Skip to content

Commit d984e6f

Browse files
fixed the code check errors and added some more QA histograms
1 parent 401a700 commit d984e6f

1 file changed

Lines changed: 76 additions & 93 deletions

File tree

PWGLF/Tasks/Resonances/deltaAnalysis.cxx

Lines changed: 76 additions & 93 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@
4040
#include <Framework/OutputObjHeader.h>
4141
#include <Framework/runDataProcessing.h>
4242

43+
#include <TH1.h>
44+
4345
#include <algorithm>
4446
#include <array>
4547
#include <cmath>
@@ -83,7 +85,7 @@ struct DeltaAnalysis {
8385
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
8486
// PDG service used only for truth-level charged-particle identification
8587
// (see isChargedPrimaryMC() below).
86-
Service<o2::framework::O2DatabasePDG> pdgDB;
88+
Service<o2::framework::O2DatabasePDG> pdgDB{};
8789

8890
struct : ConfigurableGroup {
8991
Configurable<float> cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted |z-vertex| range [cm]"};
@@ -126,8 +128,6 @@ struct DeltaAnalysis {
126128

127129
struct : ConfigurableGroup {
128130
Configurable<bool> useTPCOnlyPID{"useTPCOnlyPID", false, "Use TPC-only PID (ignore TOF even if present)"};
129-
Configurable<bool> applyTOFCutWhenAvailableBelowThreshold{"applyTOFCutWhenAvailableBelowThreshold", false, "Apply TOF PID only when TOF information exists below momentum threshold"};
130-
Configurable<float> cfgLowPtTofNsigmaCut{"cfgLowPtTofNsigmaCut", 3.0f, "Low pt nSigma TOF cut for selecting pions and protons"};
131131
Configurable<float> tpcNSigmaVetoThreshold{"tpcNSigmaVetoThreshold", 3.0f, "Reject track if TPC nSigma of a competing species is below this value"};
132132
Configurable<float> tofNSigmaVetoThreshold{"tofNSigmaVetoThreshold", 3.0f, "Reject track if TOF nSigma of a competing species is below this value"};
133133
} pidShared;
@@ -229,6 +229,7 @@ struct DeltaAnalysis {
229229
const AxisSpec ptAxis{200, 0., 10., "p_{T} (GeV/c)"};
230230
const AxisSpec massAxis{trackCuts.numberOfInvMassBins, 1.0, 8.0, "M_{inv} (GeV/#it{c}^{2})"};
231231
const AxisSpec centAxis{axes.cfgCentAxis, "Centrality (%)"};
232+
const AxisSpec centDistAxis{150, 0., 105., "Centrality (%)"};
232233
const AxisSpec vtxAxis{axes.cfgVtxAxis, "Vertex z [cm]"};
233234
const AxisSpec rapAxis{axes.cfgRapAxis, "Rapidity y"};
234235
const AxisSpec nSigmaTPCaxis{100, -10., 10., "n#sigma^{TPC}"};
@@ -247,6 +248,7 @@ struct DeltaAnalysis {
247248
histos.add("Event/hNcontributor", "PV contributors; N", kTH1F, {{2001, -0.5f, 2000.5f}});
248249
histos.add("Event/hCentrality", "Centrality", kTH1F, {centAxis});
249250
histos.add("Event/hOccupancy", "Occupancy in time range", kTH1F, {occupancyAxis});
251+
histos.add("Event/centralitydistribution", "Centrality distribution (Data);vCentFT0M;Entries", kTH1F, {centDistAxis});
250252

251253
histos.add("CentQA/hCentralityVsVtxZ", "Centrality vs vertex z", kTH2F, {vtxAxis, centAxis});
252254
histos.add("CentQA/hCentralityVsOccupancy", "Centrality vs occupancy", kTH2F, {occupancyAxis, centAxis});
@@ -404,6 +406,7 @@ struct DeltaAnalysis {
404406
histos.add("QAMC/Pion/dcaZdist", "Pion DCA_{z} distribution (MC reco, fine bins)", kTH1F, {dcaZaxis});
405407

406408
histos.add("MCRecoEvent/hRecoEvents", "Reconstructed INEL>0 events (Nrec, MC reco)", kTH1F, {centAxis});
409+
histos.add("MCRecoEvent/centralitydistribution", "Centrality distribution (MC);vCentFT0M;Entries", kTH1F, {centDistAxis});
407410
}
408411

409412
// ── processMCGen(): generated-level Delta spectra + QA ─────────────────────────────────
@@ -603,11 +606,11 @@ struct DeltaAnalysis {
603606
}
604607
fillEventCutFlowBin<Tag>(4.f); // Occupancy cut
605608

606-
const float cent = getCentrality(collision);
607-
if (cent < evSel.cfgCentMin || cent > evSel.cfgCentMax) {
608-
return false;
609-
}
610-
fillEventCutFlowBin<Tag>(5.f); // Centrality range
609+
// Centrality cut removed: events are no longer rejected based on centrality
610+
// (cfgCentMin / cfgCentMax are intentionally no longer applied here). The
611+
// centrality calculation itself, getCentrality(), cfgCentralityEstimator, and
612+
// all centrality-dependent histograms/infrastructure are left untouched.
613+
fillEventCutFlowBin<Tag>(5.f); // Centrality range (bin retained for cut-flow numbering; no longer a cut)
611614

612615
if (evSel.cfgUseNoSameBunchPileupCut &&
613616
!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
@@ -700,27 +703,26 @@ struct DeltaAnalysis {
700703
return passed && (std::abs(track.dcaZ()) < dcaCuts.cfgCutDCAz);
701704
}
702705

706+
// ── PID decision flow ────────────────────────────────────────────────────────────────────────────
703707
template <typename TrackType>
704708
bool passesProtonPID(TrackType const& track, float totalMomentum)
705709
{
706-
bool tpcPassed{false}, tofPassed{false};
707-
const int nTPCBins = static_cast<int>(mProtonTPCMomBins.size());
708-
const int nTOFBins = static_cast<int>(mProtonTOFMomBins.size());
709710
const float tpcNSigPi = std::abs(track.tpcNSigmaPi());
710711
const float tpcNSigPr = std::abs(track.tpcNSigmaPr());
711712
const float tofNSigPi = std::abs(track.tofNSigmaPi());
712713
const float tofNSigPr = std::abs(track.tofNSigmaPr());
713-
const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr;
714-
const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi;
715-
const auto circularCutSq = static_cast<float>(protonPID.combinedNSigmaCutProton * protonPID.combinedNSigmaCutProton);
716714

717-
const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold + pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold;
715+
const bool useTOF = !pidShared.useTPCOnlyPID && track.hasTOF();
718716

719-
if (!pidShared.useTPCOnlyPID && track.hasTOF()) {
720-
if (protonPID.combinedNSigmaCutProton < 0 && totalMomentum >= protonPID.minProtonMomentum) {
717+
if (useTOF) {
718+
// Combined TPC+TOF PID - used whenever TOF is available, regardless of momentum.
719+
if (protonPID.combinedNSigmaCutProton < 0) {
720+
// Asymmetric mode: independent TPC and TOF nSigma windows.
721721
if (track.tofNSigmaPr() < protonPID.minTOFNSigmaProton) {
722722
return false;
723723
}
724+
bool tofPassed = false;
725+
const int nTOFBins = static_cast<int>(mProtonTOFMomBins.size());
724726
for (int i = 0; i < nTOFBins - 1; ++i) {
725727
if (totalMomentum >= mProtonTOFMomBins[i] && totalMomentum < mProtonTOFMomBins[i + 1] &&
726728
tofNSigPr < mProtonTOFNSigCuts[i] && tofNSigPi > pidShared.tofNSigmaVetoThreshold) {
@@ -731,64 +733,53 @@ struct DeltaAnalysis {
731733
if (track.tpcNSigmaPr() < protonPID.minCombinedNSigmaProton) {
732734
return false;
733735
}
734-
if (tpcNSigPr < static_cast<float>(protonPID.maxTPCNSigmaProton) && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold) {
735-
tpcPassed = true;
736-
}
737-
} else if (protonPID.combinedNSigmaCutProton > 0 && totalMomentum >= protonPID.minProtonMomentum) {
738-
if (combinedNSigPr < circularCutSq && combinedNSigPi > circularVetoCutSq) {
739-
tpcPassed = true;
740-
tofPassed = true;
741-
}
742-
}
743-
if (totalMomentum < protonPID.minProtonMomentum && tpcNSigPr < static_cast<float>(protonPID.maxTPCNSigmaProton)) {
744-
tpcPassed = true;
745-
if (pidShared.applyTOFCutWhenAvailableBelowThreshold) {
746-
tofPassed = (tofNSigPr < pidShared.cfgLowPtTofNsigmaCut);
747-
} else {
748-
tofPassed = true;
749-
}
736+
const bool tpcPassed = tpcNSigPr < static_cast<float>(protonPID.maxTPCNSigmaProton) && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold;
737+
return tofPassed && tpcPassed;
750738
}
751-
} else {
752-
if (totalMomentum < protonPID.minProtonMomentum && tpcNSigPr < static_cast<float>(protonPID.maxTPCNSigmaProton)) {
753-
tpcPassed = true;
754-
tofPassed = true;
755-
} else {
756-
tofPassed = true;
757-
if (track.tpcNSigmaPr() < protonPID.minTPCNSigmaProton) {
758-
return false;
759-
}
760-
for (int i = 0; i < nTPCBins - 1; ++i) {
761-
if (totalMomentum >= mProtonTPCMomBins[i] && totalMomentum < mProtonTPCMomBins[i + 1] &&
762-
tpcNSigPr < mProtonTPCNSigCuts[i] && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold) {
763-
tpcPassed = true;
764-
break;
765-
}
766-
}
739+
740+
// Circular mode: joint TPC+TOF nSigma radius.
741+
const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr;
742+
const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi;
743+
const auto circularCutSq = static_cast<float>(protonPID.combinedNSigmaCutProton * protonPID.combinedNSigmaCutProton);
744+
const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold +
745+
pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold;
746+
return combinedNSigPr < circularCutSq && combinedNSigPi > circularVetoCutSq;
747+
}
748+
749+
// No usable TOF: above the momentum threshold, TOF is mandatory -> reject.
750+
if (totalMomentum >= protonPID.minProtonMomentum) {
751+
return false;
752+
}
753+
754+
// TPC-only PID below the momentum threshold, using the momentum-binned TPC cuts.
755+
const int nTPCBins = static_cast<int>(mProtonTPCMomBins.size());
756+
for (int i = 0; i < nTPCBins - 1; ++i) {
757+
if (totalMomentum >= mProtonTPCMomBins[i] && totalMomentum < mProtonTPCMomBins[i + 1]) {
758+
return tpcNSigPr < mProtonTPCNSigCuts[i] && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold;
767759
}
768760
}
769-
return tpcPassed && tofPassed;
761+
return false;
770762
}
771763

772764
template <typename TrackType>
773765
bool passesPionPID(TrackType const& track, float totalMomentum)
774766
{
775-
bool tpcPassed{false}, tofPassed{false};
776-
const int nTPCBins = static_cast<int>(mPionTPCMomBins.size());
777-
const int nTOFBins = static_cast<int>(mPionTOFMomBins.size());
778767
const float tpcNSigPi = std::abs(track.tpcNSigmaPi());
779768
const float tpcNSigPr = std::abs(track.tpcNSigmaPr());
780769
const float tofNSigPi = std::abs(track.tofNSigmaPi());
781770
const float tofNSigPr = std::abs(track.tofNSigmaPr());
782-
const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi;
783-
const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr;
784-
const auto circularCutSq = static_cast<float>(pionPID.combinedNSigmaCutPion * pionPID.combinedNSigmaCutPion);
785-
const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold + pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold;
786771

787-
if (!pidShared.useTPCOnlyPID && track.hasTOF()) {
788-
if (pionPID.combinedNSigmaCutPion < 0 && totalMomentum >= pionPID.minPionMomentum) {
772+
const bool useTOF = !pidShared.useTPCOnlyPID && track.hasTOF();
773+
774+
if (useTOF) {
775+
// Combined TPC+TOF PID - used whenever TOF is available, regardless of momentum.
776+
if (pionPID.combinedNSigmaCutPion < 0) {
777+
// Asymmetric mode: independent TPC and TOF nSigma windows.
789778
if (track.tofNSigmaPi() < pionPID.minTOFNSigmaPion) {
790779
return false;
791780
}
781+
bool tofPassed = false;
782+
const int nTOFBins = static_cast<int>(mPionTOFMomBins.size());
792783
for (int i = 0; i < nTOFBins - 1; ++i) {
793784
if (totalMomentum >= mPionTOFMomBins[i] && totalMomentum < mPionTOFMomBins[i + 1] &&
794785
tofNSigPi < mPionTOFNSigCuts[i] && tofNSigPr > pidShared.tofNSigmaVetoThreshold) {
@@ -799,42 +790,32 @@ struct DeltaAnalysis {
799790
if (track.tpcNSigmaPi() < pionPID.minCombinedNSigmaPion) {
800791
return false;
801792
}
802-
if (tpcNSigPi < static_cast<float>(pionPID.maxTPCNSigmaPion) && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold) {
803-
tpcPassed = true;
804-
}
805-
} else if (pionPID.combinedNSigmaCutPion > 0 && totalMomentum >= pionPID.minPionMomentum) {
806-
if (combinedNSigPi < circularCutSq && combinedNSigPr > circularVetoCutSq) {
807-
tpcPassed = true;
808-
tofPassed = true;
809-
}
793+
const bool tpcPassed = tpcNSigPi < static_cast<float>(pionPID.maxTPCNSigmaPion) && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold;
794+
return tofPassed && tpcPassed;
810795
}
811-
if (totalMomentum < pionPID.minPionMomentum && tpcNSigPi < static_cast<float>(pionPID.maxTPCNSigmaPion)) {
812-
tpcPassed = true;
813-
if (pidShared.applyTOFCutWhenAvailableBelowThreshold) {
814-
tofPassed = (tofNSigPi < pidShared.cfgLowPtTofNsigmaCut);
815-
} else {
816-
tofPassed = true;
817-
}
818-
}
819-
} else {
820-
if (totalMomentum < pionPID.minPionMomentum && tpcNSigPi < static_cast<float>(pionPID.maxTPCNSigmaPion)) {
821-
tpcPassed = true;
822-
tofPassed = true;
823-
} else {
824-
tofPassed = true;
825-
if (track.tpcNSigmaPi() < pionPID.minTPCNSigmaPion) {
826-
return false;
827-
}
828-
for (int i = 0; i < nTPCBins - 1; ++i) {
829-
if (totalMomentum >= mPionTPCMomBins[i] && totalMomentum < mPionTPCMomBins[i + 1] &&
830-
tpcNSigPi < mPionTPCNSigCuts[i] && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold) {
831-
tpcPassed = true;
832-
break;
833-
}
834-
}
796+
797+
// Circular mode: joint TPC+TOF nSigma radius.
798+
const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi;
799+
const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr;
800+
const auto circularCutSq = static_cast<float>(pionPID.combinedNSigmaCutPion * pionPID.combinedNSigmaCutPion);
801+
const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold +
802+
pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold;
803+
return combinedNSigPi < circularCutSq && combinedNSigPr > circularVetoCutSq;
804+
}
805+
806+
// No usable TOF: above the momentum threshold, TOF is mandatory -> reject.
807+
if (totalMomentum >= pionPID.minPionMomentum) {
808+
return false;
809+
}
810+
811+
// TPC-only PID below the momentum threshold, using the momentum-binned TPC cuts.
812+
const int nTPCBins = static_cast<int>(mPionTPCMomBins.size());
813+
for (int i = 0; i < nTPCBins - 1; ++i) {
814+
if (totalMomentum >= mPionTPCMomBins[i] && totalMomentum < mPionTPCMomBins[i + 1]) {
815+
return tpcNSigPi < mPionTPCNSigCuts[i] && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold;
835816
}
836817
}
837-
return tpcPassed && tofPassed;
818+
return false;
838819
}
839820

840821
template <typename TrackType>
@@ -1252,6 +1233,7 @@ struct DeltaAnalysis {
12521233
histos.fill(HIST("Event/hVtxZ"), collision.posZ());
12531234
histos.fill(HIST("Event/hCentrality"), centrality);
12541235
histos.fill(HIST("Event/hOccupancy"), occupancy);
1236+
histos.fill(HIST("Event/centralitydistribution"), collision.centFT0M());
12551237
histos.fill(HIST("CentQA/hCentralityVsVtxZ"), collision.posZ(), centrality);
12561238
histos.fill(HIST("CentQA/hCentralityVsOccupancy"), occupancy, centrality);
12571239
histos.fill(HIST("CentQA/hEventCountVsCentrality"), centrality);
@@ -1347,6 +1329,7 @@ struct DeltaAnalysis {
13471329

13481330
// Nrec - "the total number of selected reconstructed INEL>0 events".
13491331
histos.fill(HIST("MCRecoEvent/hRecoEvents"), centrality);
1332+
histos.fill(HIST("MCRecoEvent/centralitydistribution"), collision.centFT0M());
13501333

13511334
histos.fill(HIST("Event/hNcontributor"), collision.numContrib());
13521335
histos.fill(HIST("Event/hVtxZ"), collision.posZ());
@@ -1652,7 +1635,7 @@ struct DeltaAnalysis {
16521635
histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 3.f); // Has associated reconstructed collision
16531636
}
16541637

1655-
const bool hasAcceptedReco = acceptedMcCollisionIds.find(mcCollision.globalIndex()) != acceptedMcCollisionIds.end();
1638+
const bool hasAcceptedReco = acceptedMcCollisionIds.contains(mcCollision.globalIndex());
16561639
if (hasAcceptedReco) {
16571640
histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 4.f); // Associated reco collision passes event selection
16581641
histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 5.f); // Final EventAccepted

0 commit comments

Comments
 (0)