From 79d97ef9f4de7fee7ef3fd72bf427f43daaed67d Mon Sep 17 00:00:00 2001 From: Marvin Hemmer Date: Mon, 10 Aug 2026 10:18:45 +0200 Subject: [PATCH] [PWGEM] Fix potential missmatches between data and MC In skimmerGammaCalo merged the filling of the reconstruction data table and the MC info table into one function to ensure both are filled at the same time using the exact same cuts. In associateMCinfoPhoton added additionl collision check for filling the mclabels tables. Before it was only checked that the reco collision has a matching mc_collision nothing else. However, when creating the `fEventLabels` look up map for the association of the collisions isSelected() function is called. This could result in EmMcParticles poiting to garbage EMMCEventIds. --- .../TableProducer/associateMCinfoPhoton.cxx | 221 ++++++++++-------- .../TableProducer/skimmerGammaCalo.cxx | 73 +++--- 2 files changed, 162 insertions(+), 132 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx index 7eb9d556d85..94d7e81b609 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx @@ -73,6 +73,8 @@ struct AssociateMCInfoPhoton { kElectron = 0x8, }; + static constexpr uint16_t kNoMcParticleMask = (1 << 14); // isNoise bit: no genuine MC particle for this leg/track + Produces mcevents; Produces mceventlabels; Produces emmcparticles; @@ -90,6 +92,7 @@ struct AssociateMCInfoPhoton { Configurable maxPt{"maxPt", 20.f, "max pT for BinnedGenPts table"}; Configurable maxY{"maxY", 0.9f, "max |rapidity| for BinnedGenPts table"}; Configurable doStoreAllDaughters{"doStoreAllDaughters", false, "flag to enable storing of all photon, pi0, eta, eta' and omega daughters. This will increase the dervied data size!"}; + Configurable useIsSelected{"useIsSelected", true, "flag to enable the usage of isSelected from our EM event selection. Disabling this will increase the derived data size!"}; HistogramRegistry registry{"EMMCEvent"}; @@ -164,15 +167,34 @@ struct AssociateMCInfoPhoton { } } + template + bool passesCollisionGate(TCollision const& collisionIter) + { + if (!collisionIter.has_mcCollision()) { + return false; + } + if (useIsSelected.value && !collisionIter.isSelected()) { + return false; + } + return true; + } + template void selectMothersToStore(int motherId, int64_t mcParticleSize, TMCParticle& motherParticle, TMCParticle& daughterIter, TMCCollision const& mcCollisionIter, std::unordered_map& fNewLabels, std::map& fNewLabelsReversed, std::unordered_map& fEventIdx, std::unordered_map& fEventLabels, Counter& fCounter) { while (motherId > -1) { - if (motherId < mcParticleSize) { // protect against bad mother indices. why is this needed? + if (motherId < mcParticleSize) { motherParticle.setCursor(motherId); - int eventIdx = fEventLabels.find(mcCollisionIter.globalIndex())->second; - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + auto evIt = fEventLabels.find(mcCollisionIter.globalIndex()); + if (evIt == fEventLabels.end()) { + // Can't attribute this ancestor to a known event; stop walking rather than risk UB. + LOG(warning) << "selectMothersToStore: mcCollision " << mcCollisionIter.globalIndex() + << " not found in fEventLabels, aborting mother-chain walk."; + break; + } + int eventIdx = evIt->second; + if (!fNewLabels.contains(motherParticle.globalIndex())) { fNewLabels[motherParticle.globalIndex()] = fCounter.particles; fNewLabelsReversed[fCounter.particles] = motherParticle.globalIndex(); @@ -232,12 +254,7 @@ struct AssociateMCInfoPhoton { for (const auto& collision : collisions) { registry.fill(HIST("hEventCounter"), 1); - // TODO: investigate the collisions without corresponding mcCollision - if (!collision.has_mcCollision()) { - continue; - } - - if (!collision.isSelected()) { + if (!passesCollisionGate(collision)) { continue; } @@ -359,122 +376,135 @@ struct AssociateMCInfoPhoton { auto o2TrackIter = o2tracks.begin(); for (const auto& v0 : v0photons) { collisionIter.setCursor(v0.collisionId()); - if (!collisionIter.has_mcCollision()) { - continue; - } - mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - ele.setCursor(v0.negTrackId()); - pos.setCursor(v0.posTrackId()); + // check if the collision is valid AND the legs have valid MC particles + bool validMc = false; + if (passesCollisionGate(collisionIter)) { + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - o2TrackEle.setCursor(ele.trackId()); - o2TrackPos.setCursor(pos.trackId()); + ele.setCursor(v0.negTrackId()); + pos.setCursor(v0.posTrackId()); + o2TrackEle.setCursor(ele.trackId()); + o2TrackPos.setCursor(pos.trackId()); - if (!o2TrackEle.has_mcParticle() || !o2TrackPos.has_mcParticle()) { - continue; // If no MC particle is found, skip the v0 + validMc = o2TrackEle.has_mcParticle() && o2TrackPos.has_mcParticle(); } - for (const auto& leg : {pos, ele}) { // be carefull of order {pos, ele}! - o2TrackIter.setCursor(leg.trackId()); - mcParticleIter.setCursor(o2TrackIter.mcParticleId()); - // LOGF(info, "mcParticleIter.globalIndex() = %d, mcParticleIter.index() = %d", mcParticleIter.globalIndex(), mcParticleIter.index()); // these are exactly the same. - - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles); - if (isNew) { - fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex(); - fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; - fCounter.particles++; - } - v0legmclabels(iter->second, o2TrackIter.mcMask()); + for (const auto& leg : {pos, ele}) { // careful: order {pos, ele}, must stay fixed either way + if (validMc) { + o2TrackIter.setCursor(leg.trackId()); + mcParticleIter.setCursor(o2TrackIter.mcParticleId()); + + auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles); + if (isNew) { + fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex(); + auto evIt = fEventLabels.find(mcCollisionIter.globalIndex()); + if (evIt != fEventLabels.end()) { + fEventIdx[mcParticleIter.globalIndex()] = evIt->second; + } + fCounter.particles++; + } + v0legmclabels(iter->second, o2TrackIter.mcMask()); - // Next, store mother-chain of this reconstructed track. - int motherid = -999; // first mother index - if (mcParticleIter.has_mothers()) { - motherid = mcParticleIter.mothersIds()[0]; // first mother index + int motherid = -999; + if (mcParticleIter.has_mothers()) { + motherid = mcParticleIter.mothersIds()[0]; + } + selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, + fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); + } else { + v0legmclabels(-1, kNoMcParticleMask); // no MC match for this leg, but the row must still exist } - selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); - } // end of leg loop + } // end of loop over legs {pos, ele} } // end of v0 loop } if constexpr (static_cast(system & kElectron)) { auto o2TrackIter = o2tracks.begin(); - // auto emprimaryelectrons_coll = emprimaryelectrons.sliceBy(perCollisionEl, collision.globalIndex()); for (const auto& emprimaryelectron : emprimaryelectrons) { collisionIter.setCursor(emprimaryelectron.collisionId()); - if (!collisionIter.has_mcCollision()) { - continue; - } - mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - o2TrackIter.setCursor(emprimaryelectron.trackId()); - if (!o2TrackIter.has_mcParticle()) { - continue; // If no MC particle is found, skip the dilepton - } - mcParticleIter.setCursor(o2TrackIter.mcParticleId()); + int32_t mcLabel = -1; + uint16_t mcMask = kNoMcParticleMask; - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles); - if (isNew) { - fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex(); - fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; - fCounter.particles++; - } - emprimaryelectronmclabels(iter->second, o2TrackIter.mcMask()); + if (passesCollisionGate(collisionIter)) { + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); + o2TrackIter.setCursor(emprimaryelectron.trackId()); + + if (o2TrackIter.has_mcParticle()) { + mcParticleIter.setCursor(o2TrackIter.mcParticleId()); + + // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + auto [iter, isNew] = fNewLabels.try_emplace(mcParticleIter.globalIndex(), fCounter.particles); + if (isNew) { + fNewLabelsReversed[fCounter.particles] = mcParticleIter.globalIndex(); + auto evIt = fEventLabels.find(mcCollisionIter.globalIndex()); + if (evIt != fEventLabels.end()) { + fEventIdx[mcParticleIter.globalIndex()] = evIt->second; + } + fCounter.particles++; + } + mcLabel = iter->second; + mcMask = o2TrackIter.mcMask(); - // Next, store mother-chain of this reconstructed track. - int motherid = -999; // first mother index - if (mcParticleIter.has_mothers()) { - motherid = mcParticleIter.mothersIds()[0]; // first mother index + int motherid = -999; + if (mcParticleIter.has_mothers()) { + motherid = mcParticleIter.mothersIds()[0]; + } + selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, + fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); + } // if (o2TrackIter.has_mcParticle()) } - selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); + emprimaryelectronmclabels(mcLabel, mcMask); // always exactly one row per emprimaryelectron } // end of em primary electron loop } if constexpr (static_cast(system & kEMC)) { // for emc photons - // auto ememcclusters_coll = emcclusters.sliceBy(perCollisionEMC, collision.globalIndex()); for (const auto& emccluster : emcclusters) { collisionIter.setCursor(emccluster.collisionId()); - if (!collisionIter.has_mcCollision()) { - continue; - } - mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - // TODO: test - if (emccluster.mcParticleIds().size() <= 0) { - continue; - } std::vector vEmcMcParticleIds; std::vector vAmplitudes; - vEmcMcParticleIds.reserve(emccluster.mcParticleIds().size()); - vAmplitudes.reserve(emccluster.mcParticleIds().size()); - - for (size_t iCont = 0; iCont < emccluster.mcParticleIds().size(); iCont++) { - mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]); + if (passesCollisionGate(collisionIter)) { + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles); - if (isNew) { - fNewLabelsReversed[fCounter.particles] = mcPhoton.globalIndex(); - fEventIdx[mcPhoton.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; - fCounter.particles++; + if (emccluster.mcParticleIds().size() != emccluster.amplitude().size()) { + LOG(warning) << "Cluster " << emccluster.globalIndex() + << ": mcParticleIds/amplitude size mismatch (" + << emccluster.mcParticleIds().size() << " vs " << emccluster.amplitude().size() << ")"; } - vEmcMcParticleIds.emplace_back(iter->second); - vAmplitudes.emplace_back(emccluster.amplitude()[iCont]); - // ememcclustermclabels(fNewLabels.find(mcPhoton.index())->second); - - // Next, store mother-chain of this reconstructed track. - int motherid = -999; // first mother index - if (mcPhoton.has_mothers()) { - motherid = mcPhoton.mothersIds()[0]; // first mother index + size_t nCont = std::min(emccluster.mcParticleIds().size(), emccluster.amplitude().size()); + + vEmcMcParticleIds.reserve(nCont); + vAmplitudes.reserve(nCont); + + for (size_t iCont = 0; iCont < nCont; iCont++) { + mcPhoton.setCursor(emccluster.mcParticleIds()[iCont]); + // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + auto [iter, isNew] = fNewLabels.try_emplace(mcPhoton.globalIndex(), fCounter.particles); + if (isNew) { + fNewLabelsReversed[fCounter.particles] = mcPhoton.globalIndex(); + auto evIt = fEventLabels.find(mcCollisionIter.globalIndex()); + if (evIt != fEventLabels.end()) { // keep your earlier guard here too + fEventIdx[mcPhoton.globalIndex()] = evIt->second; + } + fCounter.particles++; + } + vEmcMcParticleIds.emplace_back(iter->second); + vAmplitudes.emplace_back(emccluster.amplitude()[iCont]); + + int motherid = -999; + if (mcPhoton.has_mothers()) { + motherid = mcPhoton.mothersIds()[0]; + } + selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, + fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); } - selectMothersToStore(motherid, mcParticles.size(), motherParticle, daughterIter, mcCollisionIter, fNewLabels, fNewLabelsReversed, fEventIdx, fEventLabels, fCounter); } // end of loop over mc particles of the current emc cluster + // whether or not the collision was okay, exactly one row is always written here so later the tables are ensured to have same size! ememcclustermclabels(vEmcMcParticleIds, vAmplitudes); - - } // end of em emc cluster loop + } } // Loop over the label map, create the mother/daughter relationships if these exist and write the skimmed MC stack @@ -517,7 +547,12 @@ struct AssociateMCInfoPhoton { } } - emmcparticles(fEventIdx.find(oldLabel)->second, mcParticleIter.pdgCode(), mcParticleIter.flags(), mcParticleIter.statusCode(), + auto evIt = fEventIdx.find(oldLabel); + if (evIt == fEventIdx.end()) { + LOG(warning) << "emmcparticles: no event index for particle label " << oldLabel << ", skipping."; + continue; + } + emmcparticles(evIt->second, mcParticleIter.pdgCode(), mcParticleIter.flags(), mcParticleIter.statusCode(), mothers, daughters, mcParticleIter.px(), mcParticleIter.py(), mcParticleIter.pz(), mcParticleIter.e(), mcParticleIter.vx(), mcParticleIter.vy(), mcParticleIter.vz()); diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx index 872df52551c..d749a06f9c3 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx @@ -39,6 +39,7 @@ #include #include #include +#include #include using namespace o2; @@ -84,6 +85,12 @@ struct SkimmerGammaCalo { void init(o2::framework::InitContext&) { + if ((doprocessRec || doprocessRecWithSecondaries) && + (doprocessRecMC || doprocessRecMCWithSecondaries)) { + LOG(fatal) << "processRec(WithSecondaries) and processRecMC(WithSecondaries) fill the same " + << "cluster table and must not be enabled together — this doubles MinClusters rows " + << "relative to EMCClusterMCLabels_001."; + } historeg.add("DefinitionIn", "Cluster definitions before cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); historeg.add("DefinitionOut", "Cluster definitions after cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); historeg.add("EIn", "Energy of clusters before cuts", gHistoSpecClusterE); @@ -131,6 +138,9 @@ struct SkimmerGammaCalo { template static constexpr bool HasSecondaries = !std::is_same_v; + template + static constexpr bool HasMcLabels = requires(TCluster c) { c.mcParticleIds(); }; + template void runAnalysis(TCollision const& collision, TClusters const& emcclusters, TClusterCells const& emcclustercells, TMatchedTracks const& emcmatchedtracks, TTracks const& /*tracks*/, TMatchedSecondaries const& secondaries = nullptr) { @@ -258,6 +268,16 @@ struct SkimmerGammaCalo { convertForStorage(emccluster.m02(), Observable::kM02), convertForStorage(emccluster.time(), Observable::kTime)); + if constexpr (HasMcLabels>) { + if (emccluster.mcParticleIds().size() != emccluster.amplitudeA().size()) { + LOG(warning) << "Mismatched MC label/amplitude array sizes for cluster " << emccluster.globalIndex() + << ": " << emccluster.mcParticleIds().size() << " vs " << emccluster.amplitudeA().size(); + } + std::vector mcLabels(emccluster.mcParticleIds().begin(), emccluster.mcParticleIds().end()); + std::vector amplitudes(emccluster.amplitudeA().begin(), emccluster.amplitudeA().end()); + tableEMCClusterMCLabels(mcLabels, amplitudes); + } + if (!vEta.empty()) { for (size_t iPart = 0; iPart < vEta.size(); ++iPart) { tableEmEmcMTracks(tableEmEmcClusters.lastIndex(), vEta[iPart], vPhi[iPart], vP[iPart], vPt[iPart]); @@ -302,48 +322,23 @@ struct SkimmerGammaCalo { } PROCESS_SWITCH(SkimmerGammaCalo, processRecWithSecondaries, "process reconstructed info with secondary track matching.", false); - void processMC(soa::Join::iterator const& collision, soa::Join const& emcclusters, aod::McParticles const&) + void processRecMC(soa::Join::iterator const& collision, + soa::Join const& emcclusters, + aod::EMCALClusterCells const& emcclustercells, aod::EMCALMatchedTracks const& emcmatchedtracks, + aod::FullTracks const& tracks) { - if (!collision.isSelected()) { - return; - } - - if (needEMCTrigger.value && !collision.alias_bit(kTVXinEMC)) { - return; - } - - for (const auto& emccluster : emcclusters) { + runAnalysis(collision, emcclusters, emcclustercells, emcmatchedtracks, tracks); + } + PROCESS_SWITCH(SkimmerGammaCalo, processRecMC, "process reconstructed info + MC labels at once", false); - // Definition cut - if (!(std::find(clusterDefinitions.value.begin(), clusterDefinitions.value.end(), emccluster.definition()) != clusterDefinitions.value.end())) { - continue; - } - // Energy cut - if (emccluster.energy() < minE) { - continue; - } - // timing cut - if (emccluster.time() > maxTime || emccluster.time() < minTime) { - continue; - } - // M02 cut - if (emccluster.nCells() > 1 && (emccluster.m02() > maxM02 || emccluster.m02() < minM02)) { - continue; - } - std::vector mcLabels; - std::vector amplitudes; - mcLabels.reserve(emccluster.amplitudeA().size()); - amplitudes.reserve(emccluster.amplitudeA().size()); - for (size_t iCont = 0; iCont < emccluster.amplitudeA().size(); iCont++) { - mcLabels.push_back(emccluster.mcParticleIds()[iCont]); - amplitudes.push_back(emccluster.amplitudeA()[iCont]); - } - tableEMCClusterMCLabels(mcLabels, amplitudes); - mcLabels.clear(); - amplitudes.clear(); - } + void processRecMCWithSecondaries(soa::Join::iterator const& collision, + soa::Join const& emcclusters, + aod::EMCALClusterCells const& emcclustercells, aod::EMCALMatchedTracks const& emcmatchedtracks, + aod::FullTracks const& tracks, aod::EMCMatchSecs const& emcmatchedsecondaries) + { + runAnalysis(collision, emcclusters, emcclustercells, emcmatchedtracks, tracks, emcmatchedsecondaries); } - PROCESS_SWITCH(SkimmerGammaCalo, processMC, "process MC info", false); // Run this in addition to processRec for MCs to copy the cluster mc labels from the EMCALMCClusters to the skimmed EMCClusterMCLabels table + PROCESS_SWITCH(SkimmerGammaCalo, processRecMCWithSecondaries, "process reconstructed info + MC labels at once with secondary track matching", false); void processDummy(aod::Collision const&) {