Skip to content

Commit e846ce4

Browse files
akyadav1963akyadav1963
committed
PWGLF: fix Event_split and add secondary-origin split for primary fraction in ptmultCorr
1 parent a1f49d4 commit e846ce4

1 file changed

Lines changed: 101 additions & 69 deletions

File tree

‎PWGLF/Tasks/GlobalEventProperties/ptmultCorr.cxx‎

Lines changed: 101 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -35,6 +35,8 @@
3535
#include <Framework/runDataProcessing.h>
3636

3737
#include <TH1.h>
38+
#include <TH2.h>
39+
#include <TMCProcess.h>
3840
#include <TPDGCode.h>
3941

4042
#include <algorithm>
@@ -106,6 +108,7 @@ enum {
106108
kRecoOther,
107109
kRecoSecondary,
108110
kRecoWeakDecay,
111+
kRecoMaterial,
109112
kRecoFake,
110113
kRecoBkg,
111114
kRecTrkTypeend
@@ -140,6 +143,7 @@ struct PtmultCorr {
140143
Configurable<float> extraphicut2{"extraphicut2", 3.12661f, "Extra Phi cut 2"};
141144
Configurable<float> extraphicut3{"extraphicut3", 0.03f, "Extra Phi cut 3"};
142145
Configurable<float> extraphicut4{"extraphicut4", 6.253f, "Extra Phi cut 4"};
146+
Configurable<bool> isZvtxPosSelMC{"isZvtxPosSelMC", true, "Zvtx position selection for MC events"};
143147

144148
Configurable<int16_t> cfgMinNCrossedRows{"cfgMinNCrossedRows", 70, "Minimum TPC crossed rows"};
145149
Configurable<bool> cfgUseNclsPID{"cfgUseNclsPID", false, "Use NclsPID instead of NclsFound for the Ncls cut"};
@@ -179,6 +183,9 @@ struct PtmultCorr {
179183
// RCT checker instance
180184
RCTFlagsChecker rctChecker;
181185

186+
static constexpr int kZeroInt{0};
187+
static constexpr float kOne{1.0f};
188+
182189
void init(InitContext const&)
183190
{
184191
if (requireGoodRct) {
@@ -305,15 +312,13 @@ struct PtmultCorr {
305312
{
306313
histos.fill(HIST("EventHist"), 1);
307314

308-
if (selHasBC && !col.has_foundBC())
309-
return false;
310-
if (selHasFT0 && !col.has_foundFT0())
315+
if (!col.sel8()) {
311316
return false;
312-
317+
}
313318
// TVX trigger (replaces sel8 as primary trigger requirement)
314-
if (!col.selection_bit(o2::aod::evsel::kIsTriggerTVX))
315-
return false;
316-
histos.fill(HIST("EventHist"), 4);
319+
// if (!col.selection_bit(o2::aod::evsel::kIsTriggerTVX))
320+
// return false;
321+
// histos.fill(HIST("EventHist"), 4);
317322

318323
// ITS ROF border
319324
if (!col.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
@@ -665,25 +670,30 @@ struct PtmultCorr {
665670
break;
666671
}
667672
} else {
668-
pid = kRecoSecondary;
669-
}
670-
if (mcpart.has_mothers()) {
671-
auto mcpartMother = mcpart.template mothers_as<aod::McParticles>().front();
672-
if (mcpartMother.pdgCode() == PDG_t::kK0Short || std::abs(mcpartMother.pdgCode()) == PDG_t::kLambda0) {
673+
// non-primary → split weak-decay vs material by production mechanism
674+
if (mcpart.getProcess() == TMCProcess::kPDecay) {
673675
pid = kRecoWeakDecay;
676+
} else {
677+
pid = kRecoMaterial;
674678
}
679+
// also fill the union bin as a cross-check (weak + material)
680+
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), kGlobalplusITS, Rectrack.dcaXY());
681+
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), trkType, Rectrack.dcaXY());
675682
}
683+
684+
// duplicate-track (split) relabel — overrides origin, as before
676685
if (find(mclabels.begin(), mclabels.end(), Rectrack.mcParticleId()) != mclabels.end()) {
677686
pid = kRecoFake;
678687
}
679688
mclabels.push_back(Rectrack.mcParticleId());
689+
680690
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(pid), kGlobalplusITS, Rectrack.dcaXY());
681691
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(pid), trkType, Rectrack.dcaXY());
682692
} else {
683693
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), Rectrack.pt(), Rectrack.phi(), static_cast<double>(kRecoBkg), kGlobalplusITS, Rectrack.dcaXY());
684694
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), Rectrack.pt(), Rectrack.phi(), static_cast<double>(kRecoBkg), trkType, Rectrack.dcaXY());
685695
}
686-
} // track (mcrec) loop
696+
}
687697
} // collision loop
688698
}
689699

@@ -716,30 +726,20 @@ struct PtmultCorr {
716726

717727
// ── CHANGE 4 ────────────────────────────────────────────────────────────
718728
// Restructured pp gen particle loop for primary fraction from MC truth.
719-
//
720-
// Pass 1 (loose): isGenChargedTrackSelected — all charged in acceptance.
721-
// → fills kGenAllCharged (denominator) when associated with a reco event.
722-
//
723-
// Pass 2 (tight): additionally require isPhysicalPrimary (via continue).
724-
// → fills kGenAll and species bins (numerator) — identical to before.
725-
//
726-
// Primary fraction (pp) = hppGenMCAssoRecdndpt[kGenAll]
727-
// / hppGenMCAssoRecdndpt[kGenAllCharged] per pT bin
729+
// Pass 1 (loose): isGenChargedTrackSelected → fills kGenAllCharged (denom).
730+
// Pass 2 (tight): + isPhysicalPrimary → fills kGenAll + species (num).
731+
// Primary fraction (pp) = kGenAll / kGenAllCharged per pT bin.
728732
// ────────────────────────────────────────────────────────────────────────
729733
for (const auto& particle : GenParticles) {
730-
// Loose check: charged + |eta| + phi (no isPhysicalPrimary requirement)
731734
if (!isGenChargedTrackSelected(particle)) {
732735
continue;
733736
}
734-
// Fill denominator: all charged particles (primaries + secondaries)
735737
if (atLeastOne) {
736738
histos.fill(HIST("hppGenMCAssoRecdndpt"), mcCollision.posZ(), particle.pt(), particle.phi(), static_cast<double>(kGenAllCharged), kNoGenpTVar);
737739
}
738-
// Tight check: physical primaries only — skip secondaries from this point on
739740
if (!particle.isPhysicalPrimary()) {
740741
continue;
741742
}
742-
// Fill primary-only histograms (numerator + species breakdown)
743743
histos.fill(HIST("hppGenMCdndpt"), mcCollision.posZ(), particle.pt(), particle.phi());
744744
if (atLeastOne) {
745745
histos.fill(HIST("hppGenMCAssoRecdndpt"), mcCollision.posZ(), particle.pt(), particle.phi(), static_cast<double>(kGenAll), kNoGenpTVar);
@@ -812,14 +812,17 @@ struct PtmultCorr {
812812
break;
813813
}
814814
} else {
815-
pid = kRecoSecondary;
816-
}
817-
if (mcpart.has_mothers()) {
818-
auto mcpartMother = mcpart.template mothers_as<aod::McParticles>().front();
819-
if (mcpartMother.pdgCode() == PDG_t::kK0Short || std::abs(mcpartMother.pdgCode()) == PDG_t::kLambda0) {
815+
// non-primary → split weak-decay vs material by production mechanism
816+
if (mcpart.getProcess() == TMCProcess::kPDecay) {
820817
pid = kRecoWeakDecay;
818+
} else {
819+
pid = kRecoMaterial;
821820
}
821+
// union cross-check bin (weak + material)
822+
histos.fill(HIST("hppRecMCdndpt"), RecCol.posZ(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), kGlobalplusITS, Rectrack.dcaXY());
823+
histos.fill(HIST("hppRecMCdndpt"), RecCol.posZ(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), trkType, Rectrack.dcaXY());
822824
}
825+
// duplicate-track (split) relabel — overrides origin, as before
823826
if (find(mclabels.begin(), mclabels.end(), Rectrack.mcParticleId()) != mclabels.end()) {
824827
pid = kRecoFake;
825828
}
@@ -833,12 +836,12 @@ struct PtmultCorr {
833836
} // track (mcrec) loop
834837
} // collision loop
835838
}
836-
837-
void processEvtLossSigLossMCpp(ColMCTrueTable::iterator const& /*mcCollision*/, ColMCRecTablepp const& RecCols, TrackMCTrueTable const& GenParticles)
839+
void processEvtLossSigLossMCpp(ColMCTrueTable::iterator const& mcCollision, ColMCRecTablepp const& RecCols, TrackMCTrueTable const& GenParticles)
838840
{
839841

840842
// ── Count generated Nch in |eta| < etaRange for event loss ──────────────
841843
int nChMC = 0;
844+
int nChMCEta1 = 0;
842845
for (const auto& particle : GenParticles) {
843846
if (!particle.isPhysicalPrimary())
844847
continue;
@@ -849,15 +852,38 @@ struct PtmultCorr {
849852
continue;
850853
if (std::abs(particle.eta()) < etaRange)
851854
nChMC++;
855+
if (std::abs(particle.eta()) > kOne)
856+
continue;
857+
nChMCEta1++;
852858
}
853859

860+
if (isZvtxPosSelMC && (std::fabs(mcCollision.posZ()) > vtxRange)) {
861+
return;
862+
}
863+
864+
//---------------------------
865+
// Select only INEL > 0 generated events?
866+
//---------------------------
867+
if (applyInelgt0) {
868+
if (!(nChMCEta1 > kZeroInt)) {
869+
return;
870+
}
871+
}
872+
873+
bool atLeastOne = false;
874+
int bestCollisionIndex = -1;
875+
int biggestNContribs = -1;
876+
854877
for (const auto& RecCol : RecCols) {
855-
if (!RecCol.has_foundBC())
856-
continue;
857-
if (!RecCol.has_foundFT0())
858-
continue;
859-
if (!RecCol.selection_bit(o2::aod::evsel::kIsTriggerTVX))
878+
879+
if (biggestNContribs < RecCol.numContrib()) {
880+
biggestNContribs = RecCol.numContrib();
881+
bestCollisionIndex = RecCol.globalIndex();
882+
}
883+
884+
if (!RecCol.sel8())
860885
continue;
886+
861887
if (!RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
862888
continue;
863889
if (!RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))
@@ -875,17 +901,7 @@ struct PtmultCorr {
875901
}
876902

877903
// ── Find best collision, apply full selection ────────────────────────────
878-
// isEventSelected() first, bestCollisionIndex filter second
879-
// matches processMCeffpp cut order
880-
bool atLeastOne = false;
881-
int bestCollisionIndex = -1;
882-
int biggestNContribs = -1;
883-
for (const auto& RecCol : RecCols) {
884-
if (biggestNContribs < RecCol.numContrib()) {
885-
biggestNContribs = RecCol.numContrib();
886-
bestCollisionIndex = RecCol.globalIndex();
887-
}
888-
}
904+
889905
for (const auto& RecCol : RecCols) {
890906
if (!isEventSelected(RecCol))
891907
continue;
@@ -924,6 +940,8 @@ struct PtmultCorr {
924940
{
925941
// ── Count generated Nch in |eta| < etaRange for event loss ──────────────
926942
int nChMC = 0;
943+
944+
int nChMCEta1 = 0;
927945
for (const auto& particle : GenParticles) {
928946
if (!particle.isPhysicalPrimary())
929947
continue;
@@ -934,17 +952,40 @@ struct PtmultCorr {
934952
continue;
935953
if (std::abs(particle.eta()) < etaRange)
936954
nChMC++;
955+
956+
if (std::abs(particle.eta()) > kOne)
957+
continue;
958+
nChMCEta1++;
959+
}
960+
961+
if (isZvtxPosSelMC && (std::fabs(mcCollision.posZ()) > vtxRange)) {
962+
return;
963+
}
964+
965+
//---------------------------
966+
// Select only INEL > 0 generated events?
967+
//---------------------------
968+
if (applyInelgt0) {
969+
if (!(nChMCEta1 > kZeroInt)) {
970+
return;
971+
}
937972
}
938973

939974
// ── Event splitting denominator ──────────────────────────────────────────
975+
bool atLeastOne = false;
976+
int bestCollisionIndex = -1;
977+
int biggestNContribs = -1;
940978

941979
for (const auto& RecCol : RecCols) {
942-
if (!RecCol.has_foundBC())
943-
continue;
944-
if (!RecCol.has_foundFT0())
945-
continue;
946-
if (!RecCol.selection_bit(o2::aod::evsel::kIsTriggerTVX))
980+
981+
if (biggestNContribs < RecCol.numContrib()) {
982+
biggestNContribs = RecCol.numContrib();
983+
bestCollisionIndex = RecCol.globalIndex();
984+
}
985+
986+
if (!RecCol.sel8())
947987
continue;
988+
948989
if (!RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
949990
continue;
950991
if (!RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))
@@ -962,24 +1003,15 @@ struct PtmultCorr {
9621003
}
9631004

9641005
// ── Find best collision, apply full selection ────────────────────────────
965-
// isEventSelected() first, bestCollisionIndex filter second
966-
// matches processMCeffPbPb cut order
967-
bool atLeastOne = false;
968-
int bestCollisionIndex = -1;
969-
int biggestNContribs = -1;
970-
for (const auto& RecCol : RecCols) {
971-
if (biggestNContribs < RecCol.numContrib()) {
972-
biggestNContribs = RecCol.numContrib();
973-
bestCollisionIndex = RecCol.globalIndex();
974-
}
975-
}
1006+
9761007
auto centrality = -999.f;
9771008
for (const auto& RecCol : RecCols) {
978-
if (!isEventSelected(RecCol))
1009+
if (RecCol.globalIndex() != bestCollisionIndex)
9791010
continue;
9801011

981-
if (RecCol.globalIndex() != bestCollisionIndex)
1012+
if (!isEventSelected(RecCol))
9821013
continue;
1014+
9831015
atLeastOne = true;
9841016
centrality = RecCol.centFT0C();
9851017
histos.fill(HIST("hCent_WRecoEvtWSelCri"), RecCol.centFT0C()); // numerator for event splitting
@@ -1003,11 +1035,11 @@ struct PtmultCorr {
10031035
if (!isGenTrackSelected(particle))
10041036
continue;
10051037

1006-
histos.fill(HIST("hPtVsNchMC_AllGen"), particle.pt(), nChMC); // denominator: all generated events
1038+
histos.fill(HIST("hPtVsNchMC_AllGen"), particle.pt(), nChMC); // denominator: all generated particles from all events
10071039
histos.fill(HIST("hgenptBeforeEvtSelPbPb"), particle.pt(), mcCollision.impactParameter());
10081040

10091041
if (atLeastOne) {
1010-
histos.fill(HIST("hPtVsNchMC_WithRecoEvt"), particle.pt(), nChMC); // numerator: ≥1 reco collision passing selection
1042+
histos.fill(HIST("hPtVsNchMC_WithRecoEvt"), particle.pt(), nChMC); // numerator: particles from events with ≥1 reco collision passing selection
10111043
histos.fill(HIST("hgenptAfterEvtSelPbPb"), particle.pt(), mcCollision.impactParameter());
10121044
}
10131045
}

0 commit comments

Comments
 (0)