Skip to content

Commit 8115124

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

1 file changed

Lines changed: 102 additions & 69 deletions

File tree

‎PWGLF/Tasks/GlobalEventProperties/ptmultCorr.cxx‎

Lines changed: 102 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,10 @@ struct PtmultCorr {
179183
// RCT checker instance
180184
RCTFlagsChecker rctChecker;
181185

186+
static constexpr int kZeroInt{0};
187+
static constexpr float kZero{0.0f};
188+
static constexpr float kOne{1.0f};
189+
182190
void init(InitContext const&)
183191
{
184192
if (requireGoodRct) {
@@ -305,15 +313,13 @@ struct PtmultCorr {
305313
{
306314
histos.fill(HIST("EventHist"), 1);
307315

308-
if (selHasBC && !col.has_foundBC())
309-
return false;
310-
if (selHasFT0 && !col.has_foundFT0())
316+
if (!col.sel8()) {
311317
return false;
312-
318+
}
313319
// 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);
320+
// if (!col.selection_bit(o2::aod::evsel::kIsTriggerTVX))
321+
// return false;
322+
// histos.fill(HIST("EventHist"), 4);
317323

318324
// ITS ROF border
319325
if (!col.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
@@ -665,25 +671,30 @@ struct PtmultCorr {
665671
break;
666672
}
667673
} 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) {
674+
// non-primary → split weak-decay vs material by production mechanism
675+
if (mcpart.getProcess() == TMCProcess::kPDecay) {
673676
pid = kRecoWeakDecay;
677+
} else {
678+
pid = kRecoMaterial;
674679
}
680+
// also fill the union bin as a cross-check (weak + material)
681+
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), kGlobalplusITS, Rectrack.dcaXY());
682+
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), trkType, Rectrack.dcaXY());
675683
}
684+
685+
// duplicate-track (split) relabel — overrides origin, as before
676686
if (find(mclabels.begin(), mclabels.end(), Rectrack.mcParticleId()) != mclabels.end()) {
677687
pid = kRecoFake;
678688
}
679689
mclabels.push_back(Rectrack.mcParticleId());
690+
680691
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(pid), kGlobalplusITS, Rectrack.dcaXY());
681692
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), mcpart.pt(), mcpart.phi(), static_cast<double>(pid), trkType, Rectrack.dcaXY());
682693
} else {
683694
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), Rectrack.pt(), Rectrack.phi(), static_cast<double>(kRecoBkg), kGlobalplusITS, Rectrack.dcaXY());
684695
histos.fill(HIST("hPbPbRecMCdndpt"), RecCol.posZ(), RecCol.centFT0C(), Rectrack.pt(), Rectrack.phi(), static_cast<double>(kRecoBkg), trkType, Rectrack.dcaXY());
685696
}
686-
} // track (mcrec) loop
697+
}
687698
} // collision loop
688699
}
689700

@@ -716,30 +727,20 @@ struct PtmultCorr {
716727

717728
// ── CHANGE 4 ────────────────────────────────────────────────────────────
718729
// 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
730+
// Pass 1 (loose): isGenChargedTrackSelected → fills kGenAllCharged (denom).
731+
// Pass 2 (tight): + isPhysicalPrimary → fills kGenAll + species (num).
732+
// Primary fraction (pp) = kGenAll / kGenAllCharged per pT bin.
728733
// ────────────────────────────────────────────────────────────────────────
729734
for (const auto& particle : GenParticles) {
730-
// Loose check: charged + |eta| + phi (no isPhysicalPrimary requirement)
731735
if (!isGenChargedTrackSelected(particle)) {
732736
continue;
733737
}
734-
// Fill denominator: all charged particles (primaries + secondaries)
735738
if (atLeastOne) {
736739
histos.fill(HIST("hppGenMCAssoRecdndpt"), mcCollision.posZ(), particle.pt(), particle.phi(), static_cast<double>(kGenAllCharged), kNoGenpTVar);
737740
}
738-
// Tight check: physical primaries only — skip secondaries from this point on
739741
if (!particle.isPhysicalPrimary()) {
740742
continue;
741743
}
742-
// Fill primary-only histograms (numerator + species breakdown)
743744
histos.fill(HIST("hppGenMCdndpt"), mcCollision.posZ(), particle.pt(), particle.phi());
744745
if (atLeastOne) {
745746
histos.fill(HIST("hppGenMCAssoRecdndpt"), mcCollision.posZ(), particle.pt(), particle.phi(), static_cast<double>(kGenAll), kNoGenpTVar);
@@ -812,14 +813,17 @@ struct PtmultCorr {
812813
break;
813814
}
814815
} 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) {
816+
// non-primary → split weak-decay vs material by production mechanism
817+
if (mcpart.getProcess() == TMCProcess::kPDecay) {
820818
pid = kRecoWeakDecay;
819+
} else {
820+
pid = kRecoMaterial;
821821
}
822+
// union cross-check bin (weak + material)
823+
histos.fill(HIST("hppRecMCdndpt"), RecCol.posZ(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), kGlobalplusITS, Rectrack.dcaXY());
824+
histos.fill(HIST("hppRecMCdndpt"), RecCol.posZ(), mcpart.pt(), mcpart.phi(), static_cast<double>(kRecoSecondary), trkType, Rectrack.dcaXY());
822825
}
826+
// duplicate-track (split) relabel — overrides origin, as before
823827
if (find(mclabels.begin(), mclabels.end(), Rectrack.mcParticleId()) != mclabels.end()) {
824828
pid = kRecoFake;
825829
}
@@ -833,12 +837,12 @@ struct PtmultCorr {
833837
} // track (mcrec) loop
834838
} // collision loop
835839
}
836-
837-
void processEvtLossSigLossMCpp(ColMCTrueTable::iterator const& /*mcCollision*/, ColMCRecTablepp const& RecCols, TrackMCTrueTable const& GenParticles)
840+
void processEvtLossSigLossMCpp(ColMCTrueTable::iterator const& mcCollision, ColMCRecTablepp const& RecCols, TrackMCTrueTable const& GenParticles)
838841
{
839842

840843
// ── Count generated Nch in |eta| < etaRange for event loss ──────────────
841844
int nChMC = 0;
845+
int nChMCEta1 = 0;
842846
for (const auto& particle : GenParticles) {
843847
if (!particle.isPhysicalPrimary())
844848
continue;
@@ -849,15 +853,38 @@ struct PtmultCorr {
849853
continue;
850854
if (std::abs(particle.eta()) < etaRange)
851855
nChMC++;
856+
if (std::abs(particle.eta()) > kOne)
857+
continue;
858+
nChMCEta1++;
852859
}
853860

861+
if (isZvtxPosSelMC && (std::fabs(mcCollision.posZ()) > vtxRange)) {
862+
return;
863+
}
864+
865+
//---------------------------
866+
// Select only INEL > 0 generated events?
867+
//---------------------------
868+
if (applyInelgt0) {
869+
if (!(nChMCEta1 > kZeroInt)) {
870+
return;
871+
}
872+
}
873+
874+
bool atLeastOne = false;
875+
int bestCollisionIndex = -1;
876+
int biggestNContribs = -1;
877+
854878
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))
879+
880+
if (biggestNContribs < RecCol.numContrib()) {
881+
biggestNContribs = RecCol.numContrib();
882+
bestCollisionIndex = RecCol.globalIndex();
883+
}
884+
885+
if (!RecCol.sel8())
860886
continue;
887+
861888
if (!RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
862889
continue;
863890
if (!RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))
@@ -875,17 +902,7 @@ struct PtmultCorr {
875902
}
876903

877904
// ── 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-
}
905+
889906
for (const auto& RecCol : RecCols) {
890907
if (!isEventSelected(RecCol))
891908
continue;
@@ -924,6 +941,8 @@ struct PtmultCorr {
924941
{
925942
// ── Count generated Nch in |eta| < etaRange for event loss ──────────────
926943
int nChMC = 0;
944+
945+
int nChMCEta1 = 0;
927946
for (const auto& particle : GenParticles) {
928947
if (!particle.isPhysicalPrimary())
929948
continue;
@@ -934,17 +953,40 @@ struct PtmultCorr {
934953
continue;
935954
if (std::abs(particle.eta()) < etaRange)
936955
nChMC++;
956+
957+
if (std::abs(particle.eta()) > kOne)
958+
continue;
959+
nChMCEta1++;
960+
}
961+
962+
if (isZvtxPosSelMC && (std::fabs(mcCollision.posZ()) > vtxRange)) {
963+
return;
964+
}
965+
966+
//---------------------------
967+
// Select only INEL > 0 generated events?
968+
//---------------------------
969+
if (applyInelgt0) {
970+
if (!(nChMCEta1 > kZeroInt)) {
971+
return;
972+
}
937973
}
938974

939975
// ── Event splitting denominator ──────────────────────────────────────────
976+
bool atLeastOne = false;
977+
int bestCollisionIndex = -1;
978+
int biggestNContribs = -1;
940979

941980
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))
981+
982+
if (biggestNContribs < RecCol.numContrib()) {
983+
biggestNContribs = RecCol.numContrib();
984+
bestCollisionIndex = RecCol.globalIndex();
985+
}
986+
987+
if (!RecCol.sel8())
947988
continue;
989+
948990
if (!RecCol.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
949991
continue;
950992
if (!RecCol.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))
@@ -962,24 +1004,15 @@ struct PtmultCorr {
9621004
}
9631005

9641006
// ── 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-
}
1007+
9761008
auto centrality = -999.f;
9771009
for (const auto& RecCol : RecCols) {
978-
if (!isEventSelected(RecCol))
1010+
if (RecCol.globalIndex() != bestCollisionIndex)
9791011
continue;
9801012

981-
if (RecCol.globalIndex() != bestCollisionIndex)
1013+
if (!isEventSelected(RecCol))
9821014
continue;
1015+
9831016
atLeastOne = true;
9841017
centrality = RecCol.centFT0C();
9851018
histos.fill(HIST("hCent_WRecoEvtWSelCri"), RecCol.centFT0C()); // numerator for event splitting
@@ -1003,11 +1036,11 @@ struct PtmultCorr {
10031036
if (!isGenTrackSelected(particle))
10041037
continue;
10051038

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

10091042
if (atLeastOne) {
1010-
histos.fill(HIST("hPtVsNchMC_WithRecoEvt"), particle.pt(), nChMC); // numerator: ≥1 reco collision passing selection
1043+
histos.fill(HIST("hPtVsNchMC_WithRecoEvt"), particle.pt(), nChMC); // numerator: particles from events with ≥1 reco collision passing selection
10111044
histos.fill(HIST("hgenptAfterEvtSelPbPb"), particle.pt(), mcCollision.impactParameter());
10121045
}
10131046
}

0 commit comments

Comments
 (0)