From 11bc750c8631daf9ca4bae3490631d9f0c18c55e Mon Sep 17 00:00:00 2001 From: Gyula Bencedi Date: Tue, 21 Apr 2026 11:23:55 +0200 Subject: [PATCH] Fix MC process --- .../GlobalEventProperties/flattenictyPikp.cxx | 513 +++++++++++------- 1 file changed, 303 insertions(+), 210 deletions(-) diff --git a/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx index 24d3bc33bc3..8820bc3f156 100644 --- a/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx +++ b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx @@ -122,6 +122,8 @@ static constexpr std::string_view Cstatus[] = {"preSel/", "postSel/"}; static constexpr std::string_view CstatCalib[] = {"preCalib/", "postCalib/"}; static constexpr std::string_view CdEdxMcRecPrim = "/hdEdxMcRecPrim"; static constexpr std::string_view CdEdxMcRecPrimF = "Tracks/{}/hdEdxMcRecPrim"; +static constexpr std::string_view CdEdxMcRecPrimSel = "/hdEdxMcRecPrimSel"; +static constexpr std::string_view CdEdxMcRecPrimSelF = "Tracks/{}/hdEdxMcRecPrimSel"; static constexpr std::string_view CpTvsDCAxy = "/hPtVsDCAxy"; static constexpr std::string_view CpTvsDCAxyF = "Tracks/{}/hPtVsDCAxy"; static constexpr std::string_view CpTvsDCAxyAll = "/hPtVsDCAxyAll"; @@ -219,6 +221,8 @@ std::array rhoLatticeFV0{0}; std::array fv0AmplitudeWoCalib{0}; std::array, NpartChrg> hPtGenRecEvt{}; std::array, NpartChrg> hPtGenPrimRecEvt{}; +std::array, NpartChrg> hPtGenRecEvtGtZero{}; +std::array, NpartChrg> hPtGenPrimRecEvtGtZero{}; std::array, NpartChrg> hPtEffRec{}; std::array, NpartChrg> hPtEffGen{}; std::array, NpartChrg> hPtEffRecGoodCollPrim{}; @@ -227,6 +231,10 @@ std::array, NpartChrg> hPtEffRecGoodCollMat{}; std::array, NpartChrg> hPtEffGenPrim{}; std::array, NpartChrg> hPtEffGenWeak{}; std::array, NpartChrg> hPtEffGenMat{}; +std::array, NpartChrg> hPtEffGenPrimEvtSelGen{}; +std::array, NpartChrg> hPtEffGenWeakEvtSelGen{}; +std::array, NpartChrg> hPtEffGenMatEvtSelGen{}; + std::array, NpartChrg> hDCAxyRecBadCollPrim{}; std::array, NpartChrg> hDCAxyRecBadCollWeak{}; std::array, NpartChrg> hDCAxyRecBadCollMat{}; @@ -286,12 +294,15 @@ struct FlattenictyPikp { Configurable cfgRemoveNoTimeFrameBorder{"cfgRemoveNoTimeFrameBorder", false, "Bunch crossing is far from Time Frame borders"}; Configurable cfgRemoveITSROFrameBorder{"cfgRemoveITSROFrameBorder", false, "Bunch crossing is far from ITS RO Frame border"}; Configurable cfgCutVtxZ{"cfgCutVtxZ", 10.0f, "Accepted z-vertex range"}; + Configurable useZVtxCutMC{"useZVtxCutMC", true, "use Zvtx cut in MC"}; + Configurable useINELCutMC{"useINELCutMC", true, "use INEL>0 cut in MC"}; Configurable cfgRemoveNoSameBunchPileup{"cfgRemoveNoSameBunchPileup", true, "Reject collisions in case of pileup with another collision in the same foundBC"}; Configurable cfgRequireIsGoodZvtxFT0vsPV{"cfgRequireIsGoodZvtxFT0vsPV", true, "Small difference between z-vertex from PV and from FT0"}; Configurable cfgRequireIsVertexITSTPC{"cfgRequireIsVertexITSTPC", false, "At least one ITS-TPC track (reject vertices built from ITS-only tracks)"}; Configurable cfgRequirekIsVertexTOFmatched{"cfgRequirekIsVertexTOFmatched", false, "Require kIsVertexTOFmatched: at least one of vertex contributors is matched to TOF"}; Configurable useMultMCmidrap{"useMultMCmidrap", true, "use generated Nch in ∣eta∣ < 0.8"}; Configurable cfgUseInelgt0wTVX{"cfgUseInelgt0wTVX", true, "Use INEL > 0 condition with TVX trigger, i.e. FT0A and FT0C acceptance"}; + Configurable cfgRemoveSplitVertex{"cfgRemoveSplitVertex", true, "Remove split vertices"}; } evtSelOpt; struct : ConfigurableGroup { @@ -452,6 +463,7 @@ struct FlattenictyPikp { using Colls = soa::Join; using CollsGen = soa::Join; using MCColls = soa::Join; + using CollsMCExtraMult = soa::Join; using CollsGenSgn = soa::SmallGroups>; using MyPIDTracks = soa::Join; using MyLabeledTracks = soa::Join; @@ -592,6 +604,8 @@ struct FlattenictyPikp { registryQC.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelVtxZ + 1, "Vtx-z pos"); registryQC.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelINELgt0 + 1, "INEL>0"); registryQC.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelRCTFlagChecker + 1, "RCT Flag Checker"); + // Number of tracks vs centrality + registryQC.add("Events/hNchVsCent", "Measured Nch vs Cent; centrality; Nch (|#eta|<0.8)", {kTH2F, {nChAxis, multAxis}}); // FV0 QA registryQC.add("FV0/hFV0AmplWCalib", "", {kTH2F, {channelFV0Axis, amplitudeFV0}}); registryQC.add("FV0/hFV0AmplvsVtxzWoCalib", "", {kTH2F, {vtxzAxis, amplitudeFV0Sum}}); @@ -762,7 +776,7 @@ struct FlattenictyPikp { registryMC.get(HIST("Events/hEvtMcGenColls"))->GetXaxis()->SetBinLabel(3, "Reco. coll."); registryMC.get(HIST("Events/hEvtMcGenColls"))->GetXaxis()->SetBinLabel(4, "Reco. good coll."); // - registryMC.add("Events/hNchVsCent", "Gen Nch vs Cent; mult; Gen Nch (|#eta|<0.8)", {kTH2F, {multAxis, nChAxis}}); + registryMC.add("Events/hNchGenVsCent", "Gen Nch vs Cent; mult; Gen Nch (|#eta|<0.8)", {kTH2F, {nChAxis, multAxis}}); registryMC.add("Events/hVtxZRec", "MC Rec vertex z position", kTH1F, {vtxzAxis}); registryMC.add("Events/hVtxZGen", "Generated vertex z position", kTH1F, {vtxzAxis}); registryMC.add("Events/hNchTVX", "Nch in FT0A+FT0C; Nch; status", {kTH2F, {nChAxis, {2, 0, 2}}}); @@ -797,9 +811,14 @@ struct FlattenictyPikp { const std::string strID = Form("/%s/%s", (i < Npart) ? "pos" : "neg", Pid[i % Npart]); hPtGenRecEvt[i] = registryMC.add("Tracks/hPtGenRecEvt" + strID, "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtGenPrimRecEvt[i] = registryMC.add("Tracks/hPtGenPrimRecEvt" + strID, "Gen evt w/ Nrec > 0 (primary); mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); + hPtGenRecEvtGtZero[i] = registryMC.add("Tracks/hPtGenRecEvtGtZero" + strID, "Gen evt w/ Nrec > 0; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); + hPtGenPrimRecEvtGtZero[i] = registryMC.add("Tracks/hPtGenPrimRecEvtGtZero" + strID, "Gen evt w/ Nrec > 0 (primary); mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffGenPrim[i] = registryMC.add("Tracks/hPtEffGenPrim" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffGenWeak[i] = registryMC.add("Tracks/hPtEffGenWeak" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffGenMat[i] = registryMC.add("Tracks/hPtEffGenMat" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); + hPtEffGenPrimEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenPrimEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); + hPtEffGenWeakEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenWeakEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); + hPtEffGenMatEvtSelGen[i] = registryMC.add("Tracks/hPtEffGenMatEvtSelGen" + strID, "Gen evt w/o rec Evt; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffRecGoodCollPrim[i] = registryMC.add("Tracks/hPtEffRecGoodCollPrim" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffRecGoodCollWeak[i] = registryMC.add("Tracks/hPtEffRecGoodCollWeak" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); hPtEffRecGoodCollMat[i] = registryMC.add("Tracks/hPtEffRecGoodCollMat" + strID, "; mult; flat; #it{p}_{T} (GeV/#it{c})", kTHnSparseF, {multAxis, flatAxis, ptAxis}); @@ -819,6 +838,7 @@ struct FlattenictyPikp { registryMC.add({fmt::format(CpTvsDCAxyWeakAllF.data(), CspeciesAll[i]).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); registryMC.add({fmt::format(CpTvsDCAxyMatAllF.data(), CspeciesAll[i]).c_str(), "; mult; flat; #it{p}_{T} (GeV/#it{c}); DCA_{xy} (cm)", {kTHnSparseF, {multAxis, flatAxis, ptAxis, dcaXYAxis}}}); registryMC.add({fmt::format(CdEdxMcRecPrimF.data(), CspeciesAll[i]).c_str(), "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {kTHnSparseF, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + registryMC.add({fmt::format(CdEdxMcRecPrimSelF.data(), CspeciesAll[i]).c_str(), "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {kTHnSparseF, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); } // Hash list for efficiency @@ -998,7 +1018,7 @@ struct FlattenictyPikp { } template - void filldEdx(T const& tracks, V const& v0s, C const& collision, aod::BCsWithTimestamps const& /*bcs*/) + void filldEdx(T const& tracks, V const& v0s, C const& collision, aod::BCsWithTimestamps const& bcs) { if (trkSelOpt.cfgRejectTrkAtTPCSector || v0SelOpt.cfgRejectV0sAtTPCSector || applyCalibGain || applyCalibVtx) { auto bc = collision.template bc_as(); @@ -1013,6 +1033,8 @@ struct FlattenictyPikp { const float flat = fillFlat(collision); registryData.fill(HIST("Events/hFlatVsMultEst"), flat, mult); + countTracks(tracks, collision, bcs, mult); + for (const auto& track : tracks) { float dEdx = track.tpcSignal(); if (cfgFillTrackQaHist) { @@ -1351,6 +1373,31 @@ struct FlattenictyPikp { return true; } + template + int countTracks(T const& tracks, C const& collision, aod::BCsWithTimestamps const& /*bcs*/, float mult) + { + if (trkSelOpt.cfgRejectTrkAtTPCSector || v0SelOpt.cfgRejectV0sAtTPCSector || applyCalibGain || applyCalibVtx) { + auto bc = collision.template bc_as(); + int currentRun = bc.runNumber(); + if (runNumber != currentRun) { + initCCDB(bc); + runNumber = currentRun; + } + } + + auto nTrk = 0; + for (auto const& track : tracks) { + if (!isGoodTrack(track, magField)) { + continue; + } + nTrk++; + } + if (fillHis) { + registryQC.fill(HIST("Events/hNchVsCent"), nTrk, mult); + } + return nTrk; + } + void phiMod(float& phimodn, const int& mag, const int& charge) { if (mag < Cnull) // for negative polarity field @@ -1888,9 +1935,9 @@ struct FlattenictyPikp { } } - float getMultMC(MCColls::iterator const& collision) + float getMultMC(CollsMCExtraMult::iterator const& collision) { - return getMult(collision); + return getMult(collision); } template @@ -2194,12 +2241,16 @@ struct FlattenictyPikp { constexpr int ChistIdx = id + pidSgn * Npart; // LOG(debug) << "fillMCRecTrack for pidSgn '" << pidSgn << "' and id '" << static_cast(id) << " with index " << ChistIdx; const aod::McParticles::iterator& mcParticle = track.mcParticle(); - const CollsGen::iterator& collision = track.collision_as(); + const CollsGen::iterator& collision = track.collision_as>(); + // const CollsGen::iterator& collision = track.collision_as(); if (!isChrgParticle(mcParticle.pdgCode())) { return; } - if (std::abs(mcParticle.y()) > trkSelOpt.cfgRapMax) { + if (std::abs(mcParticle.eta()) > trkSelOpt.cfgTrkEtaMax) { + return; + } + if (mcParticle.pt() < trkSelOpt.cfgTrkPtMin) { return; } if (!isPID(mcParticle)) { @@ -2230,15 +2281,15 @@ struct FlattenictyPikp { } else { hPtEffRecGoodCollPrim[ChistIdx]->Fill(mult, flat, track.pt()); hPtVsDCAxyRecGoodCollPrim[ChistIdx]->Fill(track.pt(), track.dcaXY()); + if (isDCAxyCut(track)) { + hPtEffRec[ChistIdx]->Fill(mcParticle.pt()); + } } } - if (isDCAxyCut(track)) { - hPtEffRec[ChistIdx]->Fill(mcParticle.pt()); - } } - template - void fillMCGen(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) + template + void fillMCGenRecEvt(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) { static_assert(pidSgn == CnullInt || pidSgn == ConeInt); constexpr int ChistIdx = id + pidSgn * Npart; @@ -2246,263 +2297,305 @@ struct FlattenictyPikp { if (!isPID(mcParticle)) { return; } - - if constexpr (recoEvt) { + if constexpr (isGtZeroColl) { + hPtGenRecEvtGtZero[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + if (mcParticle.isPhysicalPrimary()) { + hPtGenPrimRecEvtGtZero[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + hPtEffGen[ChistIdx]->Fill(mcParticle.pt()); + } + } else { hPtGenRecEvt[ChistIdx]->Fill(mult, flat, mcParticle.pt()); if (mcParticle.isPhysicalPrimary()) { hPtGenPrimRecEvt[ChistIdx]->Fill(mult, flat, mcParticle.pt()); } + } + } + + template + void fillMCGen(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) + { + static_assert(pidSgn == CnullInt || pidSgn == ConeInt); + constexpr int ChistIdx = id + pidSgn * Npart; + + if (!isPID(mcParticle)) { return; } - if (!mcParticle.isPhysicalPrimary()) { - if (mcParticle.getProcess() == CprocessIdWeak) { - hPtEffGenWeak[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + if constexpr (evtSel) { + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == CprocessIdWeak) { + hPtEffGenWeakEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + } else { + hPtEffGenMatEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + } } else { - hPtEffGenMat[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + hPtEffGenPrimEvtSelGen[ChistIdx]->Fill(mult, flat, mcParticle.pt()); } } else { - hPtEffGenPrim[ChistIdx]->Fill(mult, flat, mcParticle.pt()); - hPtEffGen[ChistIdx]->Fill(mcParticle.pt()); + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == CprocessIdWeak) { + hPtEffGenWeak[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + } else { + hPtEffGenMat[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + } + } else { + hPtEffGenPrim[ChistIdx]->Fill(mult, flat, mcParticle.pt()); + // hPtEffGen[ChistIdx]->Fill(mcParticle.pt()); + } } } Preslice perCollTrk = aod::track::collisionId; - PresliceUnsorted perCollMcLabel = aod::mccollisionlabel::mcCollisionId; - Preslice perCollMcPart = aod::mcparticle::mcCollisionId; - void processMC(MCColls const& mcCollisions, - CollsGen const& collisions, + void processMC(CollsMCExtraMult::iterator const& mcCollision, + soa::SmallGroups const& collisions, aod::BCsWithTimestamps const& /*bcs*/, - MyLabeledPIDTracks const& tracks, aod::FV0As const& /*fv0s*/, - aod::McParticles const& mcparticles) + aod::McParticles const& particles, + MyLabeledPIDTracks const& tracks) { - registryMC.fill(HIST("Events/hEvtGenRec"), 1.f, mcCollisions.size()); - registryMC.fill(HIST("Events/hEvtGenRec"), 2.f, collisions.size()); + LOGP(debug, "MC col {} has {} reco cols", mcCollision.globalIndex(), collisions.size()); + auto multMC = -1.; + if (evtSelOpt.useMultMCmidrap || multEst == 2) { // use generated Nch in ∣eta∣ < 0.8 + multMC = countPart(particles); + } else { + multMC = getMultMC(mcCollision); // using McCentFT0Ms + } + const float flatMC = fillFlatMC(particles); + registryMC.fill(HIST("Events/hFlatMCGen"), flatMC); - // Generated collisions - for (const auto& mcCollision : mcCollisions) { - if (mcCollision.isInelGt0()) { - registryMC.fill(HIST("Events/hEvtGenRec"), 3.f); + // Loop on rec collisions + // Obtain here: Denominator of tracking efficiency; Numerator event and signal loss + // + bool gtZeroColl = false; + auto multRecGt1 = -999; + auto flatRec = -999; + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; } - registryMC.fill(HIST("Events/hEvtMcGenColls"), 1); + registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0"), multRecGt1, flatRec); // Evt split den + if (evtSelOpt.cfgRemoveSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvt"), multRecGt1, flatRec); // Evt split num, w/ Nrec > 0 + gtZeroColl = true; + multRecGt1 = getMult(collision); + flatRec = fillFlat(collision); + } - const auto groupedColls = collisions.sliceBy(perCollMcLabel, mcCollision.globalIndex()); - const auto groupedParts = mcparticles.sliceBy(perCollMcPart, mcCollision.globalIndex()); - const float flatMC = fillFlatMC(groupedParts); - registryMC.fill(HIST("Events/hFlatMCGen"), flatMC); + if (gtZeroColl) { + registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvtSel"), multRecGt1, flatRec); // Evt split num, w/ Nrec > 0 + Evt. sel + registryMC.fill(HIST("Events/hNchGenVsCent"), multMC, multRecGt1); + registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0wRecEvtSel"), multMC, flatMC); // Evt loss num, w/ Nrec > 0 + Evt. sel + } - auto multMC = -1.; - if (evtSelOpt.useMultMCmidrap || multEst == 2) { // use generated Nch in ∣eta∣ < 0.8 - multMC = countPart(groupedParts); + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.eta()) > trkSelOpt.cfgTrkEtaMax) { + continue; + } + if (particle.pt() < trkSelOpt.cfgTrkPtMin) { + continue; + } + if (gtZeroColl) { + static_for<0, 1>([&](auto pidSgn) { + fillMCGenRecEvt(particle, multMC, flatMC); + fillMCGenRecEvt(particle, multMC, flatMC); + fillMCGenRecEvt(particle, multMC, flatMC); + }); + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTrecCollPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss num + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffGenPrimRecEvt), multRecGt1, flatRec, particle.pt()); // Tracking eff. den + } + }); } else { - multMC = getMultMC(mcCollision); // using McCentFT0Ms + static_for<0, 1>([&](auto pidSgn) { + fillMCGenRecEvt(particle, multMC, flatMC); + fillMCGenRecEvt(particle, multMC, flatMC); + fillMCGenRecEvt(particle, multMC, flatMC); + }); } + } - if (groupedColls.size() < ConeInt) { // if MC events have no rec collisions + // Loop on rec collisions + // Obtain here: Numerator of tracking efficiency; Secondary contamination correction + for (const auto& collision : collisions) { + if (trkSelOpt.cfgRejectTrkAtTPCSector || applyCalibGain || applyCalibVtx) { + auto bc = collision.bc_as(); + int currentRun = bc.runNumber(); + if (runNumber != currentRun) { + initCCDB(bc); + runNumber = currentRun; + } + } + if (!isGoodEvent(collision)) { continue; } - - auto maxNcontributors = -1; - auto bestCollIndex = -1; - for (const auto& collision : groupedColls) { - registryMC.fill(HIST("Events/hEvtMCRec"), 0.5); - if (!isGoodEvent(collision)) { + registryMC.fill(HIST("Events/hVtxZRec"), collision.posZ()); + if (evtSelOpt.cfgRemoveSplitVertex && collision.globalIndex() != mcCollision.bestCollisionIndex()) { + continue; + } + // Rec tracks; track selection w/ DCA open (for secondaries), w/ DCA close (for efficiency) + // Obtain here: DCAxy for sec contamination, MC closure + const auto& groupedTrks = tracks.sliceBy(perCollTrk, collision.globalIndex()); + int nTrk = 0; + for (const auto& track : groupedTrks) { + if (!track.has_collision()) { continue; } - registryMC.fill(HIST("Events/hEvtMCRec"), 1.5); - if (collision.isInelGt0()) { - registryMC.fill(HIST("Events/hEvtMCRec"), 2.5); + if (std::abs(track.eta()) > trkSelOpt.cfgTrkEtaMax) { + continue; } - const float multRecGt1 = getMult(collision); - const float flatRec = fillFlat(collision); - registryMC.fill(HIST("Events/hFlatMCRec"), flatRec); - if (maxNcontributors < collision.numContrib()) { - maxNcontributors = collision.numContrib(); - bestCollIndex = collision.globalIndex(); + if (track.pt() < trkSelOpt.cfgTrkPtMin) { + continue; } - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0"), multRecGt1, flatRec); // Evt split den - } - registryMC.fill(HIST("Events/hEvtMcGenColls"), 2); - // - // Rec collisions w/ largest number of PV contributors - for (const auto& collision : groupedColls) { - if (trkSelOpt.cfgRejectTrkAtTPCSector || applyCalibGain || applyCalibVtx) { - auto bc = collision.bc_as(); - int currentRun = bc.runNumber(); - if (runNumber != currentRun) { - initCCDB(bc); - runNumber = currentRun; - } + if (trkSelOpt.cfgApplyNcl && track.tpcNClsFound() < trkSelOpt.cfgNclTPCMin) { + continue; } - const float multRecGt1 = getMult(collision); - const float flatRec = fillFlat(collision); - registryMC.fill(HIST("Events/hEvtMcGenColls"), 3); - if (bestCollIndex != collision.globalIndex()) { + if (trkSelOpt.cfgApplyNclPID && track.tpcNClsPID() < trkSelOpt.cfgNclPidTPCMin) { continue; } - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvt"), multRecGt1, flatRec); // Evt split num, w/ Nrec > 0 - if (!isGoodEvent(collision)) { + float phiModn = track.phi(); + phiMod(phiModn, magField, track.sign()); + if (trkSelOpt.cfgRejectTrkAtTPCSector && (track.pt() >= trkSelOpt.cfgPhiCutPtMin && phiModn < fPhiCutHigh->Eval(track.pt()) && phiModn > fPhiCutLow->Eval(track.pt()))) { continue; } - registryMC.fill(HIST("Events/hEvtMcGenColls"), 4); - registryMC.fill(HIST("Events/hCentVsFlatRecINELgt0wRecEvtSel"), multRecGt1, flatRec); // Evt split num, w/ Nrec > 0 + Evt. sel - registryMC.fill(HIST("Events/hNchVsCent"), multRecGt1, multMC); - registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0wRecEvtSel"), multMC, flatMC); // Evt loss num, w/ Nrec > 0 + Evt. sel - // - // Rec tracks; track selection w/ DCA open (for secondaries), w/ DCA close (for efficiency) - // Obtain here: DCAxy for sec contamination, closure - const auto groupedTrks = tracks.sliceBy(perCollTrk, collision.globalIndex()); - for (const auto& track : groupedTrks) { - if (!track.has_collision()) { - continue; - } - if (std::abs(track.eta()) > trkSelOpt.cfgTrkEtaMax) { - continue; - } - if (track.pt() < trkSelOpt.cfgTrkPtMin) { - continue; - } - if (trkSelOpt.cfgApplyNcl && track.tpcNClsFound() < trkSelOpt.cfgNclTPCMin) { - continue; - } - if (trkSelOpt.cfgApplyNclPID && track.tpcNClsPID() < trkSelOpt.cfgNclPidTPCMin) { - continue; - } - float phiModn = track.phi(); - phiMod(phiModn, magField, track.sign()); - if (trkSelOpt.cfgRejectTrkAtTPCSector && (track.pt() >= trkSelOpt.cfgPhiCutPtMin && phiModn < fPhiCutHigh->Eval(track.pt()) && phiModn > fPhiCutLow->Eval(track.pt()))) { - continue; - } - if (!isDCAxyWoCut(track)) { - continue; - } - if (!track.has_mcParticle()) { - registryMC.fill(HIST("Tracks/hPtFakes"), multMC, flatMC, track.pt()); - continue; - } - auto particle = track.mcParticle_as(); - static_for<0, 1>([&](auto pidSgn) { - fillMCRecTrack(track, multMC, flatMC); - fillMCRecTrack(track, multMC, flatMC); - fillMCRecTrack(track, multMC, flatMC); - }); - static_for<0, 4>([&](auto i) { - constexpr int Cidx = i.value; - if (std::sqrt(std::pow(std::fabs(o2::aod::pidutils::tpcNSigma(track)), 2) + std::pow(std::fabs(o2::aod::pidutils::tofNSigma(track)), 2) < trkSelOpt.cfgDcaNsigmaCombinedMax)) { - if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { - if (!particle.isPhysicalPrimary()) { - if (particle.getProcess() == CprocessIdWeak) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyWeakAll), multMC, flatMC, track.pt(), track.dcaXY()); - } else { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyMatAll), multMC, flatMC, track.pt(), track.dcaXY()); - } - } else { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyPrimAll), multMC, flatMC, track.pt(), track.dcaXY()); - } - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyAll), multMC, flatMC, track.pt(), track.dcaXY()); - } - } - }); - if (isDCAxyCut(track)) { - static_for<0, 4>([&](auto i) { - constexpr int Cidx = i.value; - if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { - if (particle.isPhysicalPrimary()) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrim), track.eta(), multMC, flatMC, track.p(), track.tpcSignal()); - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffPrimRecEvt), multRecGt1, flatRec, track.pt()); // Tracking eff. num - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureRec), multMC, flatMC, track.pt()); - } - } - }); - } + if (!isDCAxyWoCut(track)) { + continue; } - - if (std::abs(mcCollision.posZ()) > evtSelOpt.cfgCutVtxZ) { + if (!track.has_mcParticle()) { + registryMC.fill(HIST("Tracks/hPtFakes"), multMC, flatMC, track.pt()); continue; } - registryMC.fill(HIST("Events/hVtxZRec"), collision.posZ()); - registryMC.fill(HIST("Events/hVtxZGen"), mcCollision.posZ()); - if (!pwglf::isINELgt0mc(groupedParts, pdg)) { + auto particle = track.mcParticle_as(); + if (collision.mcCollisionId() != particle.mcCollisionId()) { continue; } - - for (const auto& particle : groupedParts) { - if (!isChrgParticle(particle.pdgCode())) { - continue; - } - if (!particle.isPhysicalPrimary()) { - continue; - } - if (std::abs(particle.eta()) > trkSelOpt.cfgTrkEtaMax) { - continue; - } - if (particle.pt() < trkSelOpt.cfgTrkPtMin) { - continue; + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (std::abs(particle.eta()) > trkSelOpt.cfgTrkEtaMax) { + continue; + } + if (particle.pt() < trkSelOpt.cfgTrkPtMin) { + continue; + } + static_for<0, 1>([&](auto pidSgn) { // for checking purposes only: use gen Nch, gen Flat + fillMCRecTrack(track, multMC, flatMC); + fillMCRecTrack(track, multMC, flatMC); + fillMCRecTrack(track, multMC, flatMC); + }); + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + if (std::sqrt(std::pow(std::fabs(o2::aod::pidutils::tpcNSigma(track)), 2) + std::pow(std::fabs(o2::aod::pidutils::tofNSigma(track)), 2) < trkSelOpt.cfgDcaNsigmaCombinedMax)) { + if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { + if (!particle.isPhysicalPrimary()) { + if (particle.getProcess() == CprocessIdWeak) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyWeakAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); + } else { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyMatAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); + } + } else { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyPrimAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrim), track.eta(), multRecGt1, flatRec, track.p(), track.tpcSignal()); + } + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTvsDCAxyAll), multRecGt1, flatRec, track.pt(), track.dcaXY()); + } } - static_for<0, 1>([&](auto pidSgn) { - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - fillMCGen(particle, multMC, flatMC); - }); + }); + if (isDCAxyCut(track)) { static_for<0, 4>([&](auto i) { constexpr int Cidx = i.value; if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTrecCollPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss num - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffGenPrimRecEvt), multRecGt1, flatRec, particle.pt()); // Tracking eff. den + if (particle.isPhysicalPrimary()) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CdEdxMcRecPrimSel), track.eta(), multRecGt1, flatRec, track.p(), track.tpcSignal()); + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTeffPrimRecEvt), multRecGt1, flatRec, track.pt()); // Tracking eff. num + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureRec), multMC, flatMC, track.pt()); // closure + } } }); + nTrk++; } } - // Generated collisions without rec collision requirement - // Obtain here: signal loss den, event loss den, closure - // - registryMC.fill(HIST("Events/hEvtMcGen"), 0.5); - if (std::abs(mcCollision.posZ()) > evtSelOpt.cfgCutVtxZ) { + registryQC.fill(HIST("Events/hNchVsCent"), nTrk, multRecGt1); + } + static_for<0, 1>([&](auto pidSgn) { + fillEfficiency(); + fillEfficiency(); + fillEfficiency(); + }); + + // Loop on generated particles (no requirement on availaability of reconstructed collision; no event selection) + // + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { continue; } - registryMC.fill(HIST("Events/hEvtMcGen"), 1.5); - if (!pwglf::isINELgt0mc(groupedParts, pdg)) { + if (std::abs(particle.eta()) > trkSelOpt.cfgTrkEtaMax) { continue; } - registryMC.fill(HIST("Events/hEvtMcGen"), 2.5); - if (evtSelOpt.cfgUseInelgt0wTVX && !isInelGt0wTVX(groupedParts)) { // TVX trigger: FT0A + FT0C acceptance + if (particle.pt() < trkSelOpt.cfgTrkPtMin) { continue; } - registryMC.fill(HIST("Events/hEvtMcGen"), 3.5); - registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0"), multMC, flatMC); // Evt loss den + static_for<0, 1>([&](auto pidSgn) { + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + }); + } - for (const auto& mcParticle : groupedParts) { - if (!isChrgParticle(mcParticle.pdgCode())) { - continue; - } - if (std::abs(mcParticle.eta()) > trkSelOpt.cfgTrkEtaMax) { - continue; - } - if (mcParticle.pt() < trkSelOpt.cfgTrkPtMin) { - continue; - } - static_for<0, 1>([&](auto pidSgn) { - fillMCGen(mcParticle, multMC, flatMC); - fillMCGen(mcParticle, multMC, flatMC); - fillMCGen(mcParticle, multMC, flatMC); - }); - static_for<0, 4>([&](auto i) { - constexpr int Cidx = i.value; - // LOG(debug) << "fillMCGen for pidSgn '" << pidSgn << "' and id '" << static_cast(id) << " with index " << ChistIdx; - if (mcParticle.isPhysicalPrimary()) { - if (std::fabs(mcParticle.pdgCode()) == PDGs[Cidx]) { - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTgenPrimSgn), multMC, flatMC, mcParticle.pt()); // Sgn loss den - registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureGenPrim), multMC, flatMC, mcParticle.pt()); - } - } - }); + // Obtain here: Denominator of signal loss and event loss; MC closure + // + registryMC.fill(HIST("Events/hEvtMcGen"), 0.5); + if (evtSelOpt.useZVtxCutMC && std::abs(mcCollision.posZ()) > evtSelOpt.cfgCutVtxZ) { + return; + } + registryMC.fill(HIST("Events/hVtxZGen"), mcCollision.posZ()); + registryMC.fill(HIST("Events/hEvtMcGen"), 1.5); + if (evtSelOpt.useINELCutMC) { + if (!o2::pwglf::isINELgt0mc(particles, pdg)) { + return; + } + } + registryMC.fill(HIST("Events/hEvtMcGen"), 2.5); + if (evtSelOpt.cfgUseInelgt0wTVX && !isInelGt0wTVX(particles)) { // TVX trigger: FT0A + FT0C acceptance + return; + } + registryMC.fill(HIST("Events/hEvtMcGen"), 3.5); + registryMC.fill(HIST("Events/hNchVsFlatGenINELgt0"), multMC, flatMC); // Evt loss den + + for (const auto& particle : particles) { + if (!isChrgParticle(particle.pdgCode())) { + continue; + } + if (std::abs(particle.eta()) > trkSelOpt.cfgTrkEtaMax) { + continue; + } + if (particle.pt() < trkSelOpt.cfgTrkPtMin) { + continue; } static_for<0, 1>([&](auto pidSgn) { - fillEfficiency(); - fillEfficiency(); - fillEfficiency(); + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + }); + static_for<0, 4>([&](auto i) { + constexpr int Cidx = i.value; + // LOG(debug) << "fillMCGen for pidSgn '" << pidSgn << "' and id '" << static_cast(id) << " with index " << ChistIdx; + if (particle.isPhysicalPrimary()) { + if (std::fabs(particle.pdgCode()) == PDGs[Cidx]) { + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTgenPrimSgn), multMC, flatMC, particle.pt()); // Sgn loss den + registryMC.fill(HIST(Cprefix) + HIST(CspeciesAll[Cidx]) + HIST(CpTmcClosureGenPrim), multMC, flatMC, particle.pt()); // closure + } + } }); } }