From 7fdeb63beffd9ac89181fbc1ee8bb56bd8188c52 Mon Sep 17 00:00:00 2001 From: Roman Lietava Date: Mon, 15 Dec 2025 15:54:29 +0100 Subject: [PATCH 1/3] fix for mccollision table --- PWGLF/Tasks/Strangeness/nonPromptCascade.cxx | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx index ec6aaffd8cf..a5ce6c86f6f 100644 --- a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx +++ b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx @@ -186,7 +186,6 @@ struct NonPromptCascadeTask { using TracksExtMC = soa::Join; using CollisionCandidatesRun3 = soa::Join; using CollisionCandidatesRun3MC = soa::Join; - using CollisionsWithLabel = soa::Join; using TracksWithLabel = soa::Join; Preslice perCollision = aod::track::collisionId; @@ -742,7 +741,7 @@ struct NonPromptCascadeTask { } PROCESS_SWITCH(NonPromptCascadeTask, processCascadesMC, "process cascades: MC analysis", false); - void processGenParticles(aod::McParticles const& mcParticles) + void processGenParticles(aod::McParticles const& mcParticles, aod::McCollisions const&) { for (const auto& p : mcParticles) { auto absCode = std::abs(p.pdgCode()); @@ -817,7 +816,7 @@ struct NonPromptCascadeTask { } return mult; } - void processdNdetaMC(CollisionsWithLabel const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks) + void processdNdetaMC(CollisionCandidatesRun3MC const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks) { // std::cout << "ProcNegMC" << std::endl; // Map: collision index -> reco multiplicity From 6ca3908370a753389ccc4e3504753c5b9ea05bd5 Mon Sep 17 00:00:00 2001 From: Roman Lietava Date: Thu, 18 Dec 2025 14:41:00 +0100 Subject: [PATCH 2/3] dev: fixing mem leaks and optimising speed --- PWGLF/Tasks/Strangeness/nonPromptCascade.cxx | 274 ++++++++++++------- 1 file changed, 170 insertions(+), 104 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx index a5ce6c86f6f..e82603390fa 100644 --- a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx +++ b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx @@ -281,7 +281,7 @@ struct NonPromptCascadeTask { int nBinsMultFV0 = cfgMaxMultFV0; AxisSpec multAxis = {nBinsMult, 0, cfgMaxMult, "Multiplicity FT0M"}; - AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FT0M"}; + AxisSpec multAxisFV0 = {nBinsMultFV0, 0, cfgMaxMultFV0, "Multiplicity FV0"}; AxisSpec nTracksAxis = {100, 0., 100., "NTracksGlobal"}; AxisSpec nTracksAxisMC = {100, 0., 100., "NTracksMC"}; std::vector centBinning; @@ -311,13 +311,13 @@ struct NonPromptCascadeTask { runsBinning.push_back(run); run++; } - AxisSpec centAxis{centBinning, "Centrality (%)"}; - // AxisSpec multAxis{multBinning, "Multiplicity"}; + AxisSpec centAxisFT0M{centBinning, "Centrality FT0M (%)"}; + AxisSpec centAxisFV0{centBinning, "Centrality FV0 (%)"}; AxisSpec trackAxisMC{trackBinning, "NTracks MC"}; AxisSpec trackAxis{trackBinning, "NTracks Global Reco"}; AxisSpec runsAxis{runsBinning, "Run Number"}; - mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxis, multAxis, centAxis, multAxisFV0, nTracksAxis, runsAxis}); + mRegistryMults.add("hCentMultsRuns", "hCentMultsRuns", HistType::kTHnSparseF, {centAxisFT0M, multAxis, centAxisFV0, multAxisFV0, nTracksAxis, runsAxis}); // // dN/deta // @@ -792,114 +792,180 @@ struct NonPromptCascadeTask { } PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false); - int getMCMult(aod::McParticles const& mcParticles, int mcCollId, std::vector& ptMCvec) - { - int mult = 0; - for (auto const& mcp : mcParticles) { - if (mcp.mcCollisionId() == mcCollId) { - // multiplicity definition: - bool accept = mcp.isPhysicalPrimary(); - accept = accept && (mcp.eta() < 0.5) && (mcp.eta() > -0.5); - int q = 0; - auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode()); - if (pdgEntry) { - q = int(std::round(pdgEntry->Charge() / 3.0)); - } else { - // LOG(warn) << "No pdg assuming neutral"; - } - accept = accept && (q != 0); - if (accept) { - ++mult; - ptMCvec.push_back(mcp.pt()); - } - } +// colls : Join +// tracks: Join +// mcCollisions: aod::McCollisions +// mcParticles : aod::McParticles + +void processdNdetaMC(CollisionCandidatesRun3MC const& colls, + aod::McCollisions const& mcCollisions, + aod::McParticles const& mcParticles, + TracksWithLabel const& tracks) +{ + //------------------------------------------------------------- + // MC mult for all MC coll + //-------------------------------------------------------------- + std::vector mcMult(mcCollisions.size(), 0); + for (auto const& mcp : mcParticles) { + int mcid = mcp.mcCollisionId(); + if (mcid < 0 || mcid >= (int)mcMult.size()) continue; + + // apply your primary/eta/charge definition here + if (!mcp.isPhysicalPrimary()) continue; + if (std::abs(mcp.eta()) > 0.5f) continue; + int q = 0; + if (auto pdg = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode())) { + q = int(std::round(pdg->Charge() / 3.0)); } - return mult; + if (q == 0) continue; + + ++mcMult[mcid]; } - void processdNdetaMC(CollisionCandidatesRun3MC const& colls, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TracksWithLabel const& tracks) - { - // std::cout << "ProcNegMC" << std::endl; - // Map: collision index -> reco multiplicity - std::vector recoMult(colls.size(), 0); - for (auto const& trk : tracks) { - int collId = trk.collisionId(); - // Here you can impose same track cuts as for MC (eta, primaries, etc.) - if (std::abs(trk.eta()) > 0.5f) { - continue; - } - ++recoMult[collId]; - } - std::vector isReco(mcParticles.size(), 0); - std::vector isRecoMult(mcParticles.size(), 0); - std::vector mcReconstructed(mcCollisions.size(), 0); - // std::cout << " reco cols with mc:" << colls.size() << " tracks:" << tracks.size() << " mccols:" << mcCollisions.size() << "mcParticles:" << mcParticles.size() << std::endl; - for (auto const& col : colls) { - int mcCollId = col.mcCollisionId(); // col.template mcCollision_as(); - int collId = col.globalIndex(); - // auto mc = col.mcCollision(); - // int mcId = mc.globalIndex(); - // std::cout << "globalIndex:" << mcId << " colID:" << mcCollId << std::endl; - std::vector mcptvec; - float mult = getMCMult(mcParticles, mcCollId, mcptvec); - mcReconstructed[mcCollId] = 1; - for (auto const& trk : tracks) { - int mcPid = trk.mcParticleId(); // depends on your label table - if (mcPid >= 0 && mcPid < (int)mcParticles.size()) { - isReco[mcPid] = 1; - isRecoMult[mcPid] = mult; - } else { - continue; - } - if (trk.collisionId() != collId) { - continue; // different event - } - auto mcPar = mcParticles.rawIteratorAt(mcPid); - // Apply same acceptance as in MC multiplicity - if (!mcPar.isPhysicalPrimary()) { - continue; - } - if (std::abs(mcPar.eta()) > 0.5f) { - continue; - } - int q = 0; - auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode()); - if (pdgEntry) { - q = int(std::round(pdgEntry->Charge() / 3.0)); - } - if (q == 0) { - continue; - } - // float multReco = recoMult[collId]; - float multReco = col.multNTracksGlobal(); - float ptReco = trk.pt(); - float ptMC = mcPar.pt(); - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco); - } - // mRegistry.fill(HIST("hNTracksMCVsTracksReco"), mult, col.multNTracksGlobal()); + // ------------------------------------------------------------ + // Build mapping: (aod::Collisions row id used by tracks.collisionId()) + // -> dense index in 'colls' (0..colls.size()-1) + // We assume col.globalIndex() refers to the original aod::Collisions row. + // ------------------------------------------------------------ + int maxCollRowId = -1; + for (auto const& trk : tracks) { + maxCollRowId = std::max(maxCollRowId, (int)trk.collisionId()); + } + std::vector collRowIdToDense(maxCollRowId + 1, -1); + + int dense = 0; + for (auto const& col : colls) { + const int collRowId = col.globalIndex(); // row id in aod::Collisions + if (collRowId >= 0 && collRowId < (int)collRowIdToDense.size()) { + collRowIdToDense[collRowId] = dense; } - // count mc particles with no reco tracks - for (auto const& mcp : mcParticles) { - int mcPidG = mcp.globalIndex(); - // std::cout << "mcPidG:" << mcPidG << std::endl; - if (!isReco[mcPidG]) { - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[mcPidG], mcp.pt()); - } + ++dense; + } + + // ------------------------------------------------------------ + // Reco multiplicity per *dense collision index in colls* + // ------------------------------------------------------------ + std::vector recoMultDense(colls.size(), 0); + for (auto const& trk : tracks) { + if (std::abs(trk.eta()) > 0.5f) { + continue; } - // count unreconstructed mc collisions - for (auto const& mc : mcCollisions) { - int gindex = mc.globalIndex(); - // std::cout << "mc globalIndex:" << gindex << std::endl; - if (!mcReconstructed[gindex]) { - std::vector mcptvec; - int mult = getMCMult(mcParticles, gindex, mcptvec); - // std::cout << "===> unreconstructed:" << mult << std::endl; - for (auto const& pt : mcptvec) { - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt); - } + const int collRowId = trk.collisionId(); + if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { + continue; + } + const int dIdx = collRowIdToDense[collRowId]; + if (dIdx >= 0) { + ++recoMultDense[dIdx]; + } + } + + // ------------------------------------------------------------ + // MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex() + // ------------------------------------------------------------ + std::vector isReco(mcParticles.size(), 0); + std::vector isRecoMult(mcParticles.size(), 0.f); + std::vector mcReconstructed(mcCollisions.size(), 0); + + // Optional cache of MC multiplicity per MC collision + std::vector mcMultCache(mcCollisions.size(), -1.f); + + // ------------------------------------------------------------ + // Single pass over tracks: fill RM for tracks whose collision is in colls + // ------------------------------------------------------------ + for (auto const& trk : tracks) { + // Accept reco track + if (std::abs(trk.eta()) > 0.5f) { + continue; + } + + // Map track's collision row id -> dense colls index + const int collRowId = trk.collisionId(); + if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { + continue; + } + const int dIdx = collRowIdToDense[collRowId]; + if (dIdx < 0) { + continue; // this track's collision is not in our 'colls' view + } + + // Get the collision row (dense index in colls view) + auto col = colls.rawIteratorAt(dIdx); + + // MC collision id (row index in aod::McCollisions) + const int mcCollId = col.mcCollisionId(); + if (mcCollId < 0 || mcCollId >= (int)mcCollisions.size()) { + continue; + } + mcReconstructed[mcCollId] = 1; + + // MC particle id (row index in aod::McParticles) + const int mcPid = trk.mcParticleId(); + if (mcPid < 0 || mcPid >= (int)mcParticles.size()) { + continue; + } + + // MC multiplicity for that MC collision (cache) + float mult = mcMultCache[mcCollId]; + if (mult < 0.f) { + std::vector tmp; + mult = mcMult[mcCollId]; + mcMultCache[mcCollId] = mult; + } + + auto mcPar = mcParticles.rawIteratorAt(mcPid); + + // Apply the same acceptance as in MC multiplicity definition + if (!mcPar.isPhysicalPrimary()) { + continue; + } + if (std::abs(mcPar.eta()) > 0.5f) { + continue; + } + + int q = 0; + if (auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode())) { + q = int(std::round(pdgEntry->Charge() / 3.0)); + } + if (q == 0) { + continue; + } + + // Mark reconstructed MC particle (now that it truly passed & matched) + isReco[mcPid] = 1; + isRecoMult[mcPid] = mult; + + const float multReco = col.multNTracksGlobal(); // or recoMultDense[dIdx] + const float ptReco = trk.pt(); + const float ptMC = mcPar.pt(); + + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco); + } + + // ------------------------------------------------------------ + // MC particles with no reco track (iterate by row index) + // ------------------------------------------------------------ + for (int pid = 0; pid < (int)mcParticles.size(); ++pid) { + if (!isReco[pid]) { + auto mcp = mcParticles.rawIteratorAt(pid); + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[pid], mcp.pt()); + } + } + + // ------------------------------------------------------------ + // Unreconstructed MC collisions (iterate by row index) + // ------------------------------------------------------------ + for (int mcid = 0; mcid < (int)mcCollisions.size(); ++mcid) { + if (!mcReconstructed[mcid]) { + std::vector mcptvec; + const int mult = mcMult[mcid]; + for (auto const& pt : mcptvec) { + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt); } } } +} + PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false); }; From ac4ad1fadc62e31d61415bfa0881495e75ba08b9 Mon Sep 17 00:00:00 2001 From: Roman Lietava Date: Thu, 18 Dec 2025 14:41:58 +0100 Subject: [PATCH 3/3] clang --- PWGLF/Tasks/Strangeness/nonPromptCascade.cxx | 300 ++++++++++--------- 1 file changed, 152 insertions(+), 148 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx index e82603390fa..4784b9e95ae 100644 --- a/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx +++ b/PWGLF/Tasks/Strangeness/nonPromptCascade.cxx @@ -792,179 +792,183 @@ struct NonPromptCascadeTask { } PROCESS_SWITCH(NonPromptCascadeTask, processCascadesData, "process cascades: Data analysis", false); -// colls : Join -// tracks: Join -// mcCollisions: aod::McCollisions -// mcParticles : aod::McParticles - -void processdNdetaMC(CollisionCandidatesRun3MC const& colls, - aod::McCollisions const& mcCollisions, - aod::McParticles const& mcParticles, - TracksWithLabel const& tracks) -{ - //------------------------------------------------------------- - // MC mult for all MC coll - //-------------------------------------------------------------- - std::vector mcMult(mcCollisions.size(), 0); - for (auto const& mcp : mcParticles) { - int mcid = mcp.mcCollisionId(); - if (mcid < 0 || mcid >= (int)mcMult.size()) continue; - - // apply your primary/eta/charge definition here - if (!mcp.isPhysicalPrimary()) continue; - if (std::abs(mcp.eta()) > 0.5f) continue; - int q = 0; - if (auto pdg = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode())) { - q = int(std::round(pdg->Charge() / 3.0)); - } - if (q == 0) continue; - - ++mcMult[mcid]; - } + // colls : Join + // tracks: Join + // mcCollisions: aod::McCollisions + // mcParticles : aod::McParticles + + void processdNdetaMC(CollisionCandidatesRun3MC const& colls, + aod::McCollisions const& mcCollisions, + aod::McParticles const& mcParticles, + TracksWithLabel const& tracks) + { + //------------------------------------------------------------- + // MC mult for all MC coll + //-------------------------------------------------------------- + std::vector mcMult(mcCollisions.size(), 0); + for (auto const& mcp : mcParticles) { + int mcid = mcp.mcCollisionId(); + if (mcid < 0 || mcid >= (int)mcMult.size()) + continue; - // ------------------------------------------------------------ - // Build mapping: (aod::Collisions row id used by tracks.collisionId()) - // -> dense index in 'colls' (0..colls.size()-1) - // We assume col.globalIndex() refers to the original aod::Collisions row. - // ------------------------------------------------------------ - int maxCollRowId = -1; - for (auto const& trk : tracks) { - maxCollRowId = std::max(maxCollRowId, (int)trk.collisionId()); - } - std::vector collRowIdToDense(maxCollRowId + 1, -1); + // apply your primary/eta/charge definition here + if (!mcp.isPhysicalPrimary()) + continue; + if (std::abs(mcp.eta()) > 0.5f) + continue; + int q = 0; + if (auto pdg = TDatabasePDG::Instance()->GetParticle(mcp.pdgCode())) { + q = int(std::round(pdg->Charge() / 3.0)); + } + if (q == 0) + continue; - int dense = 0; - for (auto const& col : colls) { - const int collRowId = col.globalIndex(); // row id in aod::Collisions - if (collRowId >= 0 && collRowId < (int)collRowIdToDense.size()) { - collRowIdToDense[collRowId] = dense; + ++mcMult[mcid]; } - ++dense; - } - // ------------------------------------------------------------ - // Reco multiplicity per *dense collision index in colls* - // ------------------------------------------------------------ - std::vector recoMultDense(colls.size(), 0); - for (auto const& trk : tracks) { - if (std::abs(trk.eta()) > 0.5f) { - continue; - } - const int collRowId = trk.collisionId(); - if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { - continue; - } - const int dIdx = collRowIdToDense[collRowId]; - if (dIdx >= 0) { - ++recoMultDense[dIdx]; + // ------------------------------------------------------------ + // Build mapping: (aod::Collisions row id used by tracks.collisionId()) + // -> dense index in 'colls' (0..colls.size()-1) + // We assume col.globalIndex() refers to the original aod::Collisions row. + // ------------------------------------------------------------ + int maxCollRowId = -1; + for (auto const& trk : tracks) { + maxCollRowId = std::max(maxCollRowId, (int)trk.collisionId()); } - } + std::vector collRowIdToDense(maxCollRowId + 1, -1); - // ------------------------------------------------------------ - // MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex() - // ------------------------------------------------------------ - std::vector isReco(mcParticles.size(), 0); - std::vector isRecoMult(mcParticles.size(), 0.f); - std::vector mcReconstructed(mcCollisions.size(), 0); - - // Optional cache of MC multiplicity per MC collision - std::vector mcMultCache(mcCollisions.size(), -1.f); - - // ------------------------------------------------------------ - // Single pass over tracks: fill RM for tracks whose collision is in colls - // ------------------------------------------------------------ - for (auto const& trk : tracks) { - // Accept reco track - if (std::abs(trk.eta()) > 0.5f) { - continue; + int dense = 0; + for (auto const& col : colls) { + const int collRowId = col.globalIndex(); // row id in aod::Collisions + if (collRowId >= 0 && collRowId < (int)collRowIdToDense.size()) { + collRowIdToDense[collRowId] = dense; + } + ++dense; } - // Map track's collision row id -> dense colls index - const int collRowId = trk.collisionId(); - if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { - continue; - } - const int dIdx = collRowIdToDense[collRowId]; - if (dIdx < 0) { - continue; // this track's collision is not in our 'colls' view + // ------------------------------------------------------------ + // Reco multiplicity per *dense collision index in colls* + // ------------------------------------------------------------ + std::vector recoMultDense(colls.size(), 0); + for (auto const& trk : tracks) { + if (std::abs(trk.eta()) > 0.5f) { + continue; + } + const int collRowId = trk.collisionId(); + if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { + continue; + } + const int dIdx = collRowIdToDense[collRowId]; + if (dIdx >= 0) { + ++recoMultDense[dIdx]; + } } - // Get the collision row (dense index in colls view) - auto col = colls.rawIteratorAt(dIdx); + // ------------------------------------------------------------ + // MC bookkeeping: index by ROW INDEX (0..size-1), not globalIndex() + // ------------------------------------------------------------ + std::vector isReco(mcParticles.size(), 0); + std::vector isRecoMult(mcParticles.size(), 0.f); + std::vector mcReconstructed(mcCollisions.size(), 0); - // MC collision id (row index in aod::McCollisions) - const int mcCollId = col.mcCollisionId(); - if (mcCollId < 0 || mcCollId >= (int)mcCollisions.size()) { - continue; - } - mcReconstructed[mcCollId] = 1; + // Optional cache of MC multiplicity per MC collision + std::vector mcMultCache(mcCollisions.size(), -1.f); - // MC particle id (row index in aod::McParticles) - const int mcPid = trk.mcParticleId(); - if (mcPid < 0 || mcPid >= (int)mcParticles.size()) { - continue; - } + // ------------------------------------------------------------ + // Single pass over tracks: fill RM for tracks whose collision is in colls + // ------------------------------------------------------------ + for (auto const& trk : tracks) { + // Accept reco track + if (std::abs(trk.eta()) > 0.5f) { + continue; + } - // MC multiplicity for that MC collision (cache) - float mult = mcMultCache[mcCollId]; - if (mult < 0.f) { - std::vector tmp; - mult = mcMult[mcCollId]; - mcMultCache[mcCollId] = mult; - } + // Map track's collision row id -> dense colls index + const int collRowId = trk.collisionId(); + if (collRowId < 0 || collRowId >= (int)collRowIdToDense.size()) { + continue; + } + const int dIdx = collRowIdToDense[collRowId]; + if (dIdx < 0) { + continue; // this track's collision is not in our 'colls' view + } - auto mcPar = mcParticles.rawIteratorAt(mcPid); + // Get the collision row (dense index in colls view) + auto col = colls.rawIteratorAt(dIdx); - // Apply the same acceptance as in MC multiplicity definition - if (!mcPar.isPhysicalPrimary()) { - continue; - } - if (std::abs(mcPar.eta()) > 0.5f) { - continue; - } + // MC collision id (row index in aod::McCollisions) + const int mcCollId = col.mcCollisionId(); + if (mcCollId < 0 || mcCollId >= (int)mcCollisions.size()) { + continue; + } + mcReconstructed[mcCollId] = 1; - int q = 0; - if (auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode())) { - q = int(std::round(pdgEntry->Charge() / 3.0)); - } - if (q == 0) { - continue; - } + // MC particle id (row index in aod::McParticles) + const int mcPid = trk.mcParticleId(); + if (mcPid < 0 || mcPid >= (int)mcParticles.size()) { + continue; + } - // Mark reconstructed MC particle (now that it truly passed & matched) - isReco[mcPid] = 1; - isRecoMult[mcPid] = mult; + // MC multiplicity for that MC collision (cache) + float mult = mcMultCache[mcCollId]; + if (mult < 0.f) { + std::vector tmp; + mult = mcMult[mcCollId]; + mcMultCache[mcCollId] = mult; + } - const float multReco = col.multNTracksGlobal(); // or recoMultDense[dIdx] - const float ptReco = trk.pt(); - const float ptMC = mcPar.pt(); + auto mcPar = mcParticles.rawIteratorAt(mcPid); - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco); - } + // Apply the same acceptance as in MC multiplicity definition + if (!mcPar.isPhysicalPrimary()) { + continue; + } + if (std::abs(mcPar.eta()) > 0.5f) { + continue; + } + + int q = 0; + if (auto pdgEntry = TDatabasePDG::Instance()->GetParticle(mcPar.pdgCode())) { + q = int(std::round(pdgEntry->Charge() / 3.0)); + } + if (q == 0) { + continue; + } + + // Mark reconstructed MC particle (now that it truly passed & matched) + isReco[mcPid] = 1; + isRecoMult[mcPid] = mult; - // ------------------------------------------------------------ - // MC particles with no reco track (iterate by row index) - // ------------------------------------------------------------ - for (int pid = 0; pid < (int)mcParticles.size(); ++pid) { - if (!isReco[pid]) { - auto mcp = mcParticles.rawIteratorAt(pid); - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[pid], mcp.pt()); + const float multReco = col.multNTracksGlobal(); // or recoMultDense[dIdx] + const float ptReco = trk.pt(); + const float ptMC = mcPar.pt(); + + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRM"), mult, multReco, ptMC, ptReco); } - } - // ------------------------------------------------------------ - // Unreconstructed MC collisions (iterate by row index) - // ------------------------------------------------------------ - for (int mcid = 0; mcid < (int)mcCollisions.size(); ++mcid) { - if (!mcReconstructed[mcid]) { - std::vector mcptvec; - const int mult = mcMult[mcid]; - for (auto const& pt : mcptvec) { - mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt); + // ------------------------------------------------------------ + // MC particles with no reco track (iterate by row index) + // ------------------------------------------------------------ + for (int pid = 0; pid < (int)mcParticles.size(); ++pid) { + if (!isReco[pid]) { + auto mcp = mcParticles.rawIteratorAt(pid); + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoTrk"), isRecoMult[pid], mcp.pt()); + } + } + + // ------------------------------------------------------------ + // Unreconstructed MC collisions (iterate by row index) + // ------------------------------------------------------------ + for (int mcid = 0; mcid < (int)mcCollisions.size(); ++mcid) { + if (!mcReconstructed[mcid]) { + std::vector mcptvec; + const int mult = mcMult[mcid]; + for (auto const& pt : mcptvec) { + mRegistrydNdeta.fill(HIST("hdNdetaRM/hdNdetaRMNotInRecoCol"), mult, pt); + } } } } -} PROCESS_SWITCH(NonPromptCascadeTask, processdNdetaMC, "process mc dN/deta", false); };