From b8076e42d14e36d5497037cbd9ed09b4b4e427b2 Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Tue, 11 Aug 2026 13:29:21 +0200 Subject: [PATCH 1/3] add new deduplication modes + MC information --- .../TableProducer/photonconversionbuilder.cxx | 584 ++++++++++++++++-- PWGEM/PhotonMeson/Utils/PCMUtilities.h | 50 ++ 2 files changed, 595 insertions(+), 39 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index f91611d3c93..eb05287e955 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -74,7 +74,6 @@ #include #include #include -#include #include #include @@ -106,8 +105,10 @@ enum TrackPropMode { }; enum V0DeduplicationMode { - Pairwise = 0, // old algorithm - GreedyMatching = 1 // new algorithm + Pairwise = 0, // old algorithm, V0 with best PCA + cosPA at the same time in a group wins + GreedyMatching = 1, // new algorithm, best score wins (score combines PCA and cosPA -- see detailed implementation in PCMUtilities.cxx) + GroupMatching = 2, // new algorithm, best combination of leg-disjoint V0s wins, not the single best V0 + KeepAll = 3 // no deduplication, keep all V0s }; // Struct needed for deduplication mode 1. Used to keep track of v0 properties and a score based on pca and cosPA @@ -119,28 +120,29 @@ struct V0CandidateHelper { float cosPA = -1.f; float pca = -1.f; float score = -1.f; + float mee = 0.f; // e+e- mass V0CandidateHelper() = default; - V0CandidateHelper(int v0, int col, int pos, int ele, float cpa, float p, float s) - : v0ID(v0), colID(col), posID(pos), eleID(ele), cosPA(cpa), pca(p), score(s) + V0CandidateHelper(int v0, int col, int pos, int ele, float cpa, float p, float s, float m = 0.f) + : v0ID(v0), colID(col), posID(pos), eleID(ele), cosPA(cpa), pca(p), score(s), mee(m) { } - bool operator==(const V0CandidateHelper& other) const - { - return score == other.score; - } + bool operator==(const V0CandidateHelper& other) const { return score == other.score; } + bool operator<(const V0CandidateHelper& other) const { return score < other.score; } + bool operator>(const V0CandidateHelper& other) const { return score > other.score; } +}; - bool operator<(const V0CandidateHelper& other) const - { - return score < other.score; - } +// (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) +using CandKey = std::tuple; - bool operator>(const V0CandidateHelper& other) const - { - return score > other.score; - } +// Everything the truth diagnosis need, +struct DedupDiag { + std::map> blockersByKey; // rejected candidate + std::vector snapshot; // all candidates, before deduplication + std::vector stored; // stored candidates after deduplication + std::map> lossMargin; }; struct PhotonConversionBuilder { @@ -164,8 +166,10 @@ struct PhotonConversionBuilder { Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; Configurable useMatCorrType{"useMatCorrType", 0, "0: none, 1: TGeo, 2: LUT"}; Configurable modeTrackPropagation{"modeTrackPropagation", 0, "0: use real track propagation, including material, 1: use fast approximation using only geometry, 2: Use real track propagation and make comparison to fast propagation (only for debugging and testing)"}; - Configurable deduplicationMode{"deduplicationMode", 0, "0: Pairwise deduplication, 1: Based on Greedy matching (best score wins)"}; + Configurable deduplicationMode{"deduplicationMode", 0, "0: Pairwise deduplication, 1: Based on Greedy matching (best score wins), 2: Based on Group matching (crossed pairs are both kept), 3: Keep all V0s"}; Configurable deduplicationScoreWeight{"deduplicationScoreWeight", 0.5f, "0.:only pca goes into the score, 1: only cosPA goes itno score, any number in between is a mix of pca and cosPA"}; + Configurable cfgDedupTruthMaps{"cfgDedupTruthMaps", true, "fill the fine (dEta, dPhi, q) truth maps in MC"}; + Configurable dedupMaxGroupSize{"dedupMaxGroupSize", 12, "group matching (mode 2): conflict groups up to this size are solved exactly, larger ones greedily (capped at 16)"}; // single track cuts Configurable min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"}; @@ -286,6 +290,13 @@ struct PhotonConversionBuilder { maxSnp = 0.85f; // could be changed later maxStep = 2.00f; // could be changed later + static constexpr std::array kDedupNames = {"pairwise", "greedy matching", "group matching", "keep all"}; + if (deduplicationMode < 0 || deduplicationMode >= static_cast(kDedupNames.size())) { + LOG(fatal) << "unknown deduplicationMode " << deduplicationMode.value; + } + LOGF(info, "photon-conversion-builder: deduplicationMode = %d (%s), score weight = %.2f", + deduplicationMode.value, kDedupNames[deduplicationMode.value], deduplicationScoreWeight.value); + ccdb->setURL(ccdburl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -393,6 +404,34 @@ struct PhotonConversionBuilder { if (deduplicationScoreWeight < 0.f) { // o2-linter: disable=magic-number (score has to be above zero) LOG(warning) << "deduplicationScoreWeight is smaller than zero which is not allowed"; } + + if (doprocessMC) { + const AxisSpec axQ{60, 0.f, 0.3f, "q_{inv}^{true} (GeV/c)"}; + const AxisSpec axDEta{80, -1.6f, 1.6f, "#Delta#eta_{#gamma#gamma}^{true}"}; + const AxisSpec axClass{3, -0.5f, 2.5f, "0 = true, 1 = cross-leg fake, 2 = other fake"}; + + for (const auto& nm : {"before", "bothSurvive", "bothSurviveSameColl", "oneLost", "bothLost", "lostByCrossFake", "lostByPartnerFake", "lostByOtherFake", "lostByTrue"}) { + const std::string dir = std::string("MCDedup/Pairs/") + nm + "/"; + registry.add((dir + "hQ").c_str(), "truth photon pairs with >=1 true V0 candidate before deduplication;q_{inv}^{true} (GeV/c);pairs", kTH1F, {axQ}, true); + registry.add((dir + "hDEta").c_str(), "truth photon pairs with >=1 true V0 candidate before deduplication;#Delta#eta^{true};pairs", kTH1F, {axDEta}, true); + } + registry.add("MCDedup/Candidates/hClassStage", "candidate composition;class;0 = before, 1 = stored", kTH2F, {axClass, {2, -0.5f, 1.5f}}, true); + + const AxisSpec axBlocker{4, -0.5f, 3.5f, "0 = true, 1 = cross-leg fake, 2 = other fake, 3 = no blocker (cut, not dedup)"}; + + registry.add("MCDedup/Photons/hFate", "fate of the true photons;0 = alive, 1 = lost (blocked), 2 = lost (no blocker);photons", kTH1F, {{3, -0.5f, 2.5f}}, true); + registry.add("MCDedup/Photons/hBlockerClass", "blockers of a killed true photon (ENTRIES = blockers, not photons);class of the blocker;blocker entries", kTH1F, {axBlocker}, true); + registry.add("MCDedup/Photons/hKillMargin", "mode 0/1 only;S_{victim} - S_{winner};kills", kTH1F, {{200, 0.f, 1.f}}, true); + registry.add("MCDedup/Photons/hVictimVsWinnerPCA", "PCA of the two;PCA_{victim} (cm);PCA_{winner} (cm)", kTH2F, {{60, 0.f, 3.f}, {60, 0.f, 3.f}}, true); + registry.add("MCDedup/Photons/hLostMargin", "why the true photon lost (mode 2);cardinality deficit;#Sigma S(best set with it) - #Sigma S(winner)", kTH2F, {{5, -0.5f, 4.5f}, {100, 0.f, 2.f}}, true); + if (cfgDedupTruthMaps) { + const AxisSpec axDEtaFine{640, -1.6f, 1.6f, "#Delta#eta^{true}"}; + const AxisSpec axDPhi{144, -o2::constants::math::PI, o2::constants::math::PI, "#Delta#varphi^{true} (rad)"}; + for (const auto& nm : {"before", "bothSurvive", "bothLost", "lostByPartnerFake"}) { + registry.add((std::string("MCDedup/Pairs/") + nm + "/hMap").c_str(), "truth photon pairs", kTHnSparseF, {axDEtaFine, axDPhi, axQ}, true); + } + } + } } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -906,8 +945,12 @@ struct PhotonConversionBuilder { pca_map[std::make_tuple(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex())] = v0photoncandidate.getPCA(); cospa_map[std::make_tuple(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex())] = v0photoncandidate.getCosPA(); + ROOT::Math::PxPyPzMVector vposDedup(kfp_pos_DecayVtx.GetPx(), kfp_pos_DecayVtx.GetPy(), kfp_pos_DecayVtx.GetPz(), o2::constants::physics::MassElectron); + ROOT::Math::PxPyPzMVector veleDedup(kfp_ele_DecayVtx.GetPx(), kfp_ele_DecayVtx.GetPy(), kfp_ele_DecayVtx.GetPz(), o2::constants::physics::MassElectron); + const auto meeSV = static_cast((vposDedup + veleDedup).M()); + const auto score = getScoreV0(v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), deduplicationScoreWeight); - V0CandidateHelper v0Helper(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), score); + V0CandidateHelper v0Helper(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), score, meeSV); vecV0Dedup.emplace_back(v0Helper); if (applyPCMMl) { @@ -1013,7 +1056,7 @@ struct PhotonConversionBuilder { std::unordered_map nv0_map; // map collisionId -> nv0 template - void build(TCollisions const& collisions, TV0s const& v0s, TTracks const& tracks, TBCs const&) + void build(TCollisions const& collisions, TV0s const& v0s, TTracks const& tracks, TBCs const&, DedupDiag* diag = nullptr) { for (const auto& collision : collisions) { if constexpr (isMC) { @@ -1076,14 +1119,17 @@ struct PhotonConversionBuilder { continue; } - if (collisionId != collisionId_tmp && eleId == eleId_tmp && posId == posId_tmp && cospa < cospa_tmp) { // same ele and pos, but attached to different collision - // LOGF(info, "!reject! | collision id = %d | posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , cospa1 = %f , cospa2 = %f", collisionId, posId, eleId, posId_tmp, eleId_tmp, cospa, cospa_tmp); + if (collisionId != collisionId_tmp && eleId == eleId_tmp && posId == posId_tmp && cospa < cospa_tmp) { + if (diag != nullptr) { + diag->blockersByKey[key].push_back(key_tmp); + } is_most_aligned_v0 = false; break; } - if ((eleId == eleId_tmp || posId == posId_tmp) && v0pca > v0pca_tmp) { - // LOGF(info, "!reject! | collision id = %d | posid1 = %d , eleid1 = %d , posid2 = %d , eleid2 = %d , pca1 = %f , pca2 = %f", collisionId, posId, eleId, posId_tmp, eleId_tmp, v0pca, v0pca_tmp); + if (diag != nullptr) { + diag->blockersByKey[key].push_back(key_tmp); + } is_closest_v0 = false; break; } @@ -1104,37 +1150,225 @@ struct PhotonConversionBuilder { } } // end of pca_map loop // LOGF(info, "pca_map.size() = %d", pca_map.size()); + + } else if (deduplicationMode == V0DeduplicationMode::KeepAll) { + // Every candidate that survived the quality cuts is stored + stored_v0Ids.clear(); + stored_fullv0Ids.clear(); + nv0_map.clear(); + for (const auto& v0Cand : vecV0Dedup) { + stored_v0Ids.emplace_back(v0Cand.posID, v0Cand.eleID); + stored_fullv0Ids.emplace_back(v0Cand.v0ID, v0Cand.colID, v0Cand.posID, v0Cand.eleID); + nv0_map[v0Cand.colID]++; + } + std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), + [](const auto& a, const auto& b) { + return std::get<1>(a) < std::get<1>(b); + }); + } else if (deduplicationMode == V0DeduplicationMode::GroupMatching) { + const int maxEnum = std::min(dedupMaxGroupSize.value, 16); // o2-linter: disable=magic-number (subset check of photons) + const int nCand = static_cast(vecV0Dedup.size()); + stored_v0Ids.clear(); + stored_fullv0Ids.clear(); + nv0_map.clear(); + int nGreedyFallback = 0; + + // Groups are build of V0s which share/are connected through a shared daughter track. + std::vector parent(nCand); + for (int i = 0; i < nCand; ++i) { + parent[i] = i; + } + auto findRoot = [&parent](int x) { + while (parent[x] != x) { + parent[x] = parent[parent[x]]; + x = parent[x]; + } + return x; + }; + // Remember the first candidate using each track. If another candidate uses the same track, merge their groups. + std::unordered_map firstCandOfTrack; + for (int i = 0; i < nCand; ++i) { + for (const int trackId : {vecV0Dedup[i].posID, vecV0Dedup[i].eleID}) { + auto [it, inserted] = firstCandOfTrack.try_emplace(trackId, i); + if (!inserted) { + parent[findRoot(i)] = findRoot(it->second); + } + } + } + // Collect all candidates belonging to the same group (i.e. sharing a daughter track) and sort them by score. Then select the best subset of candidates that do not share any daughter tracks. + std::map> groups; + for (int i = 0; i < nCand; ++i) { + groups[findRoot(i)].push_back(i); + } + + for (auto& [root, members] : groups) { + const int nMembers = static_cast(members.size()); + std::vector selected(nMembers, 0); + // Process candidates from best to worst score. + std::sort(members.begin(), members.end(), [this](int a, int b) { + if (vecV0Dedup[a].score != vecV0Dedup[b].score) { + return vecV0Dedup[a].score < vecV0Dedup[b].score; + } + return vecV0Dedup[a].v0ID < vecV0Dedup[b].v0ID; + }); + + if (nMembers == 1) { + // A single candidate has no conflict and can always be kept. + selected[0] = 1; + } else if (nMembers <= maxEnum) { + // For small groups, enumerate all subsets and find the optimal set of candidates that do not share any daughter tracks. + std::unordered_map legBit; + std::vector legMask(nMembers, 0); + for (int k = 0; k < nMembers; ++k) { + for (const int trackId : {vecV0Dedup[members[k]].posID, vecV0Dedup[members[k]].eleID}) { + const int bit = legBit.try_emplace(trackId, static_cast(legBit.size())).first->second; + legMask[k] |= (1ull << bit); // at most 2 * 16 = 32 legs, fits into 64 bit + } + } + uint32_t bestMask = 0; + int bestCount = 0; + float bestScore = 0.f; + std::vector bestWithCount(nMembers, -1); + std::vector bestWithScore(nMembers, 0.f); + for (uint32_t mask = 1; mask < (1u << nMembers); ++mask) { + uint64_t usedLegsLocal = 0; + float scoreSum = 0.f; + int count = 0; + bool disjoint = true; + for (int k = 0; k < nMembers; ++k) { + if ((mask & (1u << k)) == 0) { + continue; + } + if ((usedLegsLocal & legMask[k]) != 0) { + disjoint = false; + break; + } + usedLegsLocal |= legMask[k]; + scoreSum += vecV0Dedup[members[k]].score; + ++count; + } + if (!disjoint) { + continue; + } + if (count > bestCount || (count == bestCount && scoreSum < bestScore)) { + bestCount = count; + bestScore = scoreSum; + bestMask = mask; + } + if (diag != nullptr) { + for (int k = 0; k < nMembers; ++k) { + if ((mask & (1u << k)) == 0) { + continue; + } + if (count > bestWithCount[k] || (count == bestWithCount[k] && scoreSum < bestWithScore[k])) { + bestWithCount[k] = count; + bestWithScore[k] = scoreSum; + } + } + } + } + for (int k = 0; k < nMembers; ++k) { + selected[k] = ((bestMask & (1u << k)) != 0) ? 1 : 0; + if (diag != nullptr && selected[k] == 0 && bestWithCount[k] >= 0) { + const auto& cand = vecV0Dedup[members[k]]; + diag->lossMargin[std::make_tuple(static_cast(cand.v0ID), static_cast(cand.colID), + static_cast(cand.posID), static_cast(cand.eleID))] = + std::make_pair(bestCount - bestWithCount[k], bestWithScore[k] - bestScore); + } + } + } else { + + ++nGreedyFallback; + std::vector usedTracks; + usedTracks.reserve(2 * nMembers); + for (int k = 0; k < nMembers; ++k) { + const auto& cand = vecV0Dedup[members[k]]; + if (std::find(usedTracks.begin(), usedTracks.end(), cand.posID) != usedTracks.end() || + std::find(usedTracks.begin(), usedTracks.end(), cand.eleID) != usedTracks.end()) { + continue; + } + usedTracks.push_back(cand.posID); + usedTracks.push_back(cand.eleID); + selected[k] = 1; + } + } + + for (int k = 0; k < nMembers; ++k) { + const auto& cand = vecV0Dedup[members[k]]; + if (selected[k] != 0) { + stored_v0Ids.emplace_back(cand.posID, cand.eleID); + stored_fullv0Ids.emplace_back(cand.v0ID, cand.colID, cand.posID, cand.eleID); + nv0_map[cand.colID]++; + continue; + } + if (diag == nullptr) { + continue; + } + const CandKey victimKey = std::make_tuple(static_cast(cand.v0ID), static_cast(cand.colID), + static_cast(cand.posID), static_cast(cand.eleID)); + for (int j = 0; j < nMembers; ++j) { + if (selected[j] == 0) { + continue; + } + const auto& winner = vecV0Dedup[members[j]]; + if (winner.posID == cand.posID || winner.eleID == cand.eleID) { + diag->blockersByKey[victimKey].emplace_back(static_cast(winner.v0ID), static_cast(winner.colID), + static_cast(winner.posID), static_cast(winner.eleID)); + } + } + } + } + + if (nGreedyFallback > 0) { + LOGF(warning, "group matching: %d conflict group(s) larger than dedupMaxGroupSize = %d were solved greedily", nGreedyFallback, maxEnum); + } + + std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), + [](const auto& a, const auto& b) { + return std::get<1>(a) < std::get<1>(b); + }); + } else { // Sort best candidates first, depending on score std::sort(vecV0Dedup.begin(), vecV0Dedup.end()); // container to keep track of which electron and positron tracks have already been used std::vector usedLegs(tracks.size(), false); - + std::unordered_map ownerOfLeg; // diagnostics only // clear output containers stored_v0Ids.clear(); stored_fullv0Ids.clear(); nv0_map.clear(); - // Loop over all v0s, starting with the one with the best score // If a v0 contains a track that is already in another (better) V0 candidate, skip V0 - for (const auto& v0Cand : vecV0Dedup) { - // Skip if one of the tracks is already used + for (const auto& v0Cand : vecV0Dedup) + // Skip if one of the tracks is already used + { + const CandKey candKey = std::make_tuple(static_cast(v0Cand.v0ID), static_cast(v0Cand.colID), + static_cast(v0Cand.posID), static_cast(v0Cand.eleID)); if (usedLegs[v0Cand.posID] || usedLegs[v0Cand.eleID]) { + if (diag != nullptr) { + auto& blockers = diag->blockersByKey[candKey]; + for (const int legId : {v0Cand.posID, v0Cand.eleID}) { + const auto itOwner = ownerOfLeg.find(legId); + if (itOwner == ownerOfLeg.end()) { + continue; + } + if (std::find(blockers.begin(), blockers.end(), itOwner->second) == blockers.end()) { + blockers.push_back(itOwner->second); + } + } + } continue; } - // Accept candidate usedLegs[v0Cand.posID] = true; usedLegs[v0Cand.eleID] = true; - + if (diag != nullptr) { + ownerOfLeg[v0Cand.posID] = candKey; + ownerOfLeg[v0Cand.eleID] = candKey; + } stored_v0Ids.emplace_back(v0Cand.posID, v0Cand.eleID); - - stored_fullv0Ids.emplace_back( - v0Cand.v0ID, - v0Cand.colID, - v0Cand.posID, - v0Cand.eleID); - + stored_fullv0Ids.emplace_back(v0Cand.v0ID, v0Cand.colID, v0Cand.posID, v0Cand.eleID); nv0_map[v0Cand.colID]++; } @@ -1147,6 +1381,11 @@ struct PhotonConversionBuilder { // End of deduplication } + if (diag != nullptr) { + diag->snapshot = vecV0Dedup; + diag->stored = stored_fullv0Ids; + } + for (const auto& fullv0Id : stored_fullv0Ids) { auto v0Id = std::get<0>(fullv0Id); // auto collisionId = std::get<1>(fullv0Id); @@ -1189,6 +1428,270 @@ struct PhotonConversionBuilder { Filter v0Filter = o2::aod::v0::v0Type > (uint8_t)0; using FilteredV0s = soa::Filtered; + // Pair-level outcome categories. The outcome (survival or loss) and cause are listed. + enum PairFate { + kBefore = 0, + kBothSurvive, + kOneLost, + kBothLost, + kLostByCrossFake, + kLostByPartnerFake, + kLostByOtherFake, + kLostByTrue + }; + + void fillTruthPair(int fate, float q, float dEta) + { + switch (fate) { + case kBefore: + registry.fill(HIST("MCDedup/Pairs/before/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/before/hDEta"), dEta); + break; + case kBothSurvive: + registry.fill(HIST("MCDedup/Pairs/bothSurvive/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/bothSurvive/hDEta"), dEta); + break; + case kOneLost: + registry.fill(HIST("MCDedup/Pairs/oneLost/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/oneLost/hDEta"), dEta); + break; + case kBothLost: + registry.fill(HIST("MCDedup/Pairs/bothLost/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/bothLost/hDEta"), dEta); + break; + case kLostByCrossFake: + registry.fill(HIST("MCDedup/Pairs/lostByCrossFake/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/lostByCrossFake/hDEta"), dEta); + break; + case kLostByPartnerFake: + registry.fill(HIST("MCDedup/Pairs/lostByPartnerFake/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/lostByPartnerFake/hDEta"), dEta); + break; + case kLostByOtherFake: + registry.fill(HIST("MCDedup/Pairs/lostByOtherFake/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/lostByOtherFake/hDEta"), dEta); + break; + case kLostByTrue: + registry.fill(HIST("MCDedup/Pairs/lostByTrue/hQ"), q); + registry.fill(HIST("MCDedup/Pairs/lostByTrue/hDEta"), dEta); + break; + default: + break; + } + } + + void fillTruthPairMap(int fate, float dEta, float dPhi, float q) + { + if (!cfgDedupTruthMaps) { + return; + } + switch (fate) { + case kBefore: + registry.fill(HIST("MCDedup/Pairs/before/hMap"), dEta, dPhi, q); + break; + case kBothSurvive: + registry.fill(HIST("MCDedup/Pairs/bothSurvive/hMap"), dEta, dPhi, q); + break; + case kBothLost: + registry.fill(HIST("MCDedup/Pairs/bothLost/hMap"), dEta, dPhi, q); + break; + case kLostByPartnerFake: + registry.fill(HIST("MCDedup/Pairs/lostByPartnerFake/hMap"), dEta, dPhi, q); + break; + default: + break; + } + } + + template + void fillDedupTruthDiagnostics(TTracks const& tracks, TMCParticles const& mcparticles, DedupDiag const& diag) + { + const int nCand = static_cast(diag.snapshot.size()); + if (nCand == 0) { + return; + } + const std::set storedSet(diag.stored.begin(), diag.stored.end()); + + std::vector key(nCand); + std::vector cls(nCand, static_cast(kV0OtherFake)); + std::vector motherPos(nCand, -1), motherEle(nCand, -1); + std::map indexOfKey; + for (int i = 0; i < nCand; ++i) { + const auto& c = diag.snapshot[i]; + key[i] = std::make_tuple(static_cast(c.v0ID), static_cast(c.colID), + static_cast(c.posID), static_cast(c.eleID)); + indexOfKey[key[i]] = i; + cls[i] = static_cast(classifyV0Truth(tracks.rawIteratorAt(c.posID), tracks.rawIteratorAt(c.eleID), + mcparticles, motherPos[i], motherEle[i])); + registry.fill(HIST("MCDedup/Candidates/hClassStage"), static_cast(cls[i]), 0.f); + if (storedSet.contains(key[i]) > 0) { + registry.fill(HIST("MCDedup/Candidates/hClassStage"), static_cast(cls[i]), 1.f); + } + } + + struct PhotonFate { + bool alive = false; + int best = -1; // best-score candidate + std::vector cands; // all true candidates of this photon + std::set aliveColls; // collisions with a surviving candidate + }; + std::map fate; + for (int i = 0; i < nCand; ++i) { + if (cls[i] != kV0True) { + continue; + } + auto& f = fate[motherPos[i]]; + f.cands.push_back(i); + if (storedSet.contains(key[i]) > 0) { + f.alive = true; + f.aliveColls.insert(std::get<1>(key[i])); + } + if (f.best < 0 || diag.snapshot[i].score < diag.snapshot[f.best].score) { + f.best = i; + } + } + std::map> blockersOfPhoton; + for (const auto& [mid, f] : fate) { + if (f.alive) { + registry.fill(HIST("MCDedup/Photons/hFate"), 0.f); + continue; + } + auto& blockers = blockersOfPhoton[mid]; + for (const int ci : f.cands) { + const auto itBlocker = diag.blockersByKey.find(key[ci]); + if (itBlocker == diag.blockersByKey.end()) { + continue; + } + for (const auto& bKey : itBlocker->second) { + const auto itIndex = indexOfKey.find(bKey); + if (itIndex == indexOfKey.end()) { + continue; + } + if (std::find(blockers.begin(), blockers.end(), itIndex->second) == blockers.end()) { + blockers.push_back(itIndex->second); + } + } + } + + int bestDeficit = -1; + float bestExcess = 0.f; + for (const int ci : f.cands) { + const auto itMargin = diag.lossMargin.find(key[ci]); + if (itMargin == diag.lossMargin.end()) { + continue; + } + if (bestDeficit < 0 || itMargin->second.first < bestDeficit || + (itMargin->second.first == bestDeficit && itMargin->second.second < bestExcess)) { + bestDeficit = itMargin->second.first; + bestExcess = itMargin->second.second; + } + } + if (bestDeficit >= 0) { + registry.fill(HIST("MCDedup/Photons/hLostMargin"), static_cast(bestDeficit), bestExcess); + } + + registry.fill(HIST("MCDedup/Photons/hFate"), blockers.empty() ? 2.f : 1.f); + + if (blockers.empty()) { + registry.fill(HIST("MCDedup/Photons/hBlockerClass"), 3.f); // o2-linter: disable=magic-number (bin 3 = no blocker) + continue; + } + for (const int w : blockers) { + registry.fill(HIST("MCDedup/Photons/hBlockerClass"), static_cast(cls[w])); + } + + if (f.best >= 0 && deduplicationMode.value < static_cast(V0DeduplicationMode::GroupMatching)) { + const int w = blockers.front(); + registry.fill(HIST("MCDedup/Photons/hKillMargin"), diag.snapshot[f.best].score - diag.snapshot[w].score); + registry.fill(HIST("MCDedup/Photons/hVictimVsWinnerPCA"), diag.snapshot[f.best].pca, diag.snapshot[w].pca); + } + } + + std::map> mothersByMcColl; + for (const auto& [mid, f] : fate) { + mothersByMcColl[mcparticles.iteratorAt(mid).mcCollisionId()].push_back(mid); + } + constexpr float kMaxQTruth = 0.3f; + for (const auto& [mcCol, mids] : mothersByMcColl) { + for (size_t a = 0; a < mids.size(); ++a) { + for (size_t b = a + 1; b < mids.size(); ++b) { + const auto g1 = mcparticles.iteratorAt(mids[a]); + const auto g2 = mcparticles.iteratorAt(mids[b]); + const float e1 = std::hypot(g1.px(), g1.py(), g1.pz()); + const float e2 = std::hypot(g2.px(), g2.py(), g2.pz()); + const float qTrue = std::sqrt(std::max(0.f, 2.f * (e1 * e2 - g1.px() * g2.px() - g1.py() * g2.py() - g1.pz() * g2.pz()))); + if (qTrue > kMaxQTruth) { + continue; + } + const float dEta = g1.eta() - g2.eta(); + const float dPhi = RecoDecay::constrainAngle(static_cast(g1.phi() - g2.phi()), -o2::constants::math::PI); + fillTruthPair(kBefore, qTrue, dEta); + fillTruthPairMap(kBefore, dEta, dPhi, qTrue); + + const auto& f1 = fate.at(mids[a]); + const auto& f2 = fate.at(mids[b]); + const int nLost = (f1.alive ? 0 : 1) + (f2.alive ? 0 : 1); + if (nLost == 0) { + fillTruthPair(kBothSurvive, qTrue, dEta); + fillTruthPairMap(kBothSurvive, dEta, dPhi, qTrue); + bool sameColl = false; + for (const auto& c1 : f1.aliveColls) { + sameColl = sameColl || (f2.aliveColls.count(c1) > 0); + } + if (sameColl) { + registry.fill(HIST("MCDedup/Pairs/bothSurviveSameColl/hQ"), qTrue); + } + continue; + } + if (nLost == 1) { + fillTruthPair(kOneLost, qTrue, dEta); + } else { + fillTruthPair(kBothLost, qTrue, dEta); + fillTruthPairMap(kBothLost, dEta, dPhi, qTrue); + } + + bool byCross = false, byPartner = false, byOther = false, byTrue = false; + auto attribute = [&](int64_t victim, int64_t partner) { + if (fate.at(victim).alive) { + return; + } + const auto itB = blockersOfPhoton.find(victim); + if (itB == blockersOfPhoton.end()) { + return; + } + for (const int w : itB->second) { + if (cls[w] == kV0CrossLegFake) { + byCross = true; + if (motherPos[w] == partner || motherEle[w] == partner) { + byPartner = true; + } + } else if (cls[w] == kV0OtherFake) { + byOther = true; + } else if (cls[w] == kV0True) { + byTrue = true; + } + } + }; + attribute(mids[a], mids[b]); + attribute(mids[b], mids[a]); + if (byCross) { + fillTruthPair(kLostByCrossFake, qTrue, dEta); + } + if (byPartner) { + fillTruthPair(kLostByPartnerFake, qTrue, dEta); + fillTruthPairMap(kLostByPartnerFake, dEta, dPhi, qTrue); + } + if (byOther) { + fillTruthPair(kLostByOtherFake, qTrue, dEta); + } + if (byTrue) { + fillTruthPair(kLostByTrue, qTrue, dEta); + } + } + } + } + } + void processRec(MyCollisions const& collisions, FilteredV0s const& v0s, MyTracksIU const& tracks, aod::BCsWithTimestamps const& bcs) { build(collisions, v0s, tracks, bcs); @@ -1201,9 +1704,12 @@ struct PhotonConversionBuilder { // } // PROCESS_SWITCH(PhotonConversionBuilder, processRec_SWT, "process reconstructed info for data", false); - void processMC(MyCollisionsMC const& collisions, FilteredV0s const& v0s, MyTracksIUMC const& tracks, aod::BCsWithTimestamps const& bcs) + void processMC(MyCollisionsMC const& collisions, FilteredV0s const& v0s, MyTracksIUMC const& tracks, + aod::BCsWithTimestamps const& bcs, aod::McParticles const& mcparticles) { - build(collisions, v0s, tracks, bcs); + DedupDiag diag; + build(collisions, v0s, tracks, bcs, &diag); + fillDedupTruthDiagnostics(tracks, mcparticles, diag); } PROCESS_SWITCH(PhotonConversionBuilder, processMC, "process reconstructed info for MC", false); diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index 7005cbcb66a..f7453543576 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -29,11 +29,13 @@ #include // IWYU pragma: keep (for rotate) #include // IWYU pragma: keep (do not replace with Math/Vector2Dfwd.h) #include +#include #include #include #include +#include //_______________________________________________________________________ inline bool checkAP(const float alpha, const float qt, const float alpha_max = 0.95, const float qt_max = 0.05) @@ -212,5 +214,53 @@ inline float getScoreV0(float cosPA, float pca, float weight) float score = wCos * cosScore + wPca * pcaScore; return score; } +//_______________________________________________________________________ +/// \brief truth class of a V0 candidate, from the mothers of its two legs +enum V0TruthClass { + kV0True = 0, // both legs come from the SAME true photon + kV0CrossLegFake = 1, // legs come from TWO different true photons + kV0OtherFake = 2 // at least one leg is not a photon conversion leg +}; +//_______________________________________________________________________ +/// \brief global index of the photon a track descends from +/// \param track track with an MC label +/// \param mcparticles the full McParticles table +/// \return index of the mother if it is a photon, -1 otherwise +template +inline int64_t getPhotonMotherId(TTrack const& track, TMCParticles const& mcparticles) +{ + if (!track.has_mcParticle()) { + return -1; + } + const auto mcp = mcparticles.iteratorAt(track.mcParticleId()); + if (!mcp.has_mothers()) { + return -1; + } + const auto& mothers = mcp.mothersIds(); + if (mothers.empty() || mothers[0] < 0) { + return -1; + } + return (mcparticles.iteratorAt(mothers[0]).pdgCode() == PDG_t::kGamma) ? static_cast(mothers[0]) : -1; +} + +//_______________________________________________________________________ +/// \brief classify a V0 candidate by the photon mothers of its two legs +/// \param motherPos photon mother of the positive leg, -1 if none +/// \param motherEle photon mother of the negative leg, -1 if none +/// \return kV0True, kV0CrossLegFake or kV0OtherFake +template +inline V0TruthClass classifyV0Truth(TTrack const& pos, TTrack const& ele, TMCParticles const& mcparticles, + int64_t& motherPos, int64_t& motherEle) +{ + motherPos = getPhotonMotherId(pos, mcparticles); + motherEle = getPhotonMotherId(ele, mcparticles); + if (motherPos >= 0 && motherPos == motherEle) { + return kV0True; + } + if (motherPos >= 0 && motherEle >= 0) { + return kV0CrossLegFake; + } + return kV0OtherFake; +} #endif // PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_ From 4990761025c23b0b1aa94223cfceb104217e99fe Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Tue, 11 Aug 2026 13:58:07 +0200 Subject: [PATCH 2/3] fix o2-linter errors --- .../TableProducer/photonconversionbuilder.cxx | 154 +++++++++++------- 1 file changed, 98 insertions(+), 56 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index eb05287e955..9462c9ddf07 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -120,7 +120,7 @@ struct V0CandidateHelper { float cosPA = -1.f; float pca = -1.f; float score = -1.f; - float mee = 0.f; // e+e- mass + float mee = 0.f; // e+e- mass at the secondary vertex (GeV/c^2) V0CandidateHelper() = default; @@ -139,9 +139,12 @@ using CandKey = std::tuple; // Everything the truth diagnosis need, struct DedupDiag { - std::map> blockersByKey; // rejected candidate - std::vector snapshot; // all candidates, before deduplication - std::vector stored; // stored candidates after deduplication + std::map> blockersByKey; // rejected candidate -> ALL candidates that took one of its legs + std::vector snapshot; // ALL candidates, before deduplication + std::vector stored; // the survivors + // group matching only: rejected candidate -> (how many candidates the best solution CONTAINING + // it would have cost, by how much its best solution was worse in total score). + // deficit == 0 means an equally large alternative existed and the score alone decided. std::map> lossMargin; }; @@ -290,12 +293,12 @@ struct PhotonConversionBuilder { maxSnp = 0.85f; // could be changed later maxStep = 2.00f; // could be changed later - static constexpr std::array kDedupNames = {"pairwise", "greedy matching", "group matching", "keep all"}; - if (deduplicationMode < 0 || deduplicationMode >= static_cast(kDedupNames.size())) { + static constexpr std::array DedupNames = {"pairwise", "greedy matching", "group matching", "keep all"}; + if (deduplicationMode < 0 || deduplicationMode >= static_cast(DedupNames.size())) { LOG(fatal) << "unknown deduplicationMode " << deduplicationMode.value; } LOGF(info, "photon-conversion-builder: deduplicationMode = %d (%s), score weight = %.2f", - deduplicationMode.value, kDedupNames[deduplicationMode.value], deduplicationScoreWeight.value); + deduplicationMode.value, DedupNames[deduplicationMode.value], deduplicationScoreWeight.value); ccdb->setURL(ccdburl); ccdb->setCaching(true); @@ -409,8 +412,12 @@ struct PhotonConversionBuilder { const AxisSpec axQ{60, 0.f, 0.3f, "q_{inv}^{true} (GeV/c)"}; const AxisSpec axDEta{80, -1.6f, 1.6f, "#Delta#eta_{#gamma#gamma}^{true}"}; const AxisSpec axClass{3, -0.5f, 2.5f, "0 = true, 1 = cross-leg fake, 2 = other fake"}; - - for (const auto& nm : {"before", "bothSurvive", "bothSurviveSameColl", "oneLost", "bothLost", "lostByCrossFake", "lostByPartnerFake", "lostByOtherFake", "lostByTrue"}) { + // One folder per pair fate, with the SAME histogram names inside. The post-processing then + // only varies the folder, and everything belonging to one fate sits together. + // NB: the "lost..." folders are diagnostic selections, NOT an exclusive decomposition - + // lostByPartnerFake is a subset of lostByCrossFake, and a pair can appear in several of them. + for (const auto& nm : {"before", "bothSurvive", "bothSurviveSameColl", "oneLost", "bothLost", + "lostByCrossFake", "lostByPartnerFake", "lostByOtherFake", "lostByTrue"}) { const std::string dir = std::string("MCDedup/Pairs/") + nm + "/"; registry.add((dir + "hQ").c_str(), "truth photon pairs with >=1 true V0 candidate before deduplication;q_{inv}^{true} (GeV/c);pairs", kTH1F, {axQ}, true); registry.add((dir + "hDEta").c_str(), "truth photon pairs with >=1 true V0 candidate before deduplication;#Delta#eta^{true};pairs", kTH1F, {axDEta}, true); @@ -421,13 +428,21 @@ struct PhotonConversionBuilder { registry.add("MCDedup/Photons/hFate", "fate of the true photons;0 = alive, 1 = lost (blocked), 2 = lost (no blocker);photons", kTH1F, {{3, -0.5f, 2.5f}}, true); registry.add("MCDedup/Photons/hBlockerClass", "blockers of a killed true photon (ENTRIES = blockers, not photons);class of the blocker;blocker entries", kTH1F, {axBlocker}, true); + // only meaningful for mode 0/1, where the smaller score always wins. Mode 2 compares SETS and + // may deliberately keep the worse-scoring candidates, so it is not filled there. registry.add("MCDedup/Photons/hKillMargin", "mode 0/1 only;S_{victim} - S_{winner};kills", kTH1F, {{200, 0.f, 1.f}}, true); registry.add("MCDedup/Photons/hVictimVsWinnerPCA", "PCA of the two;PCA_{victim} (cm);PCA_{winner} (cm)", kTH2F, {{60, 0.f, 3.f}, {60, 0.f, 3.f}}, true); - registry.add("MCDedup/Photons/hLostMargin", "why the true photon lost (mode 2);cardinality deficit;#Sigma S(best set with it) - #Sigma S(winner)", kTH2F, {{5, -0.5f, 4.5f}, {100, 0.f, 2.f}}, true); + // mode 2 only: why exactly was this true photon not part of the winning set? + // x = 0 -> an equally large solution containing it existed, only the total score decided. + // That is a pure ranking problem and the part that a better discriminant can win back. + // x > 0 -> keeping it would have cost x candidates. Structural, no arbitration can repair it. + registry.add("MCDedup/Photons/hLostMargin", "why the true photon lost (mode 2);cardinality deficit;#Sigma S(best set with it) - #Sigma S(winner)", + kTH2F, {{5, -0.5f, 4.5f}, {100, 0.f, 2.f}}, true); if (cfgDedupTruthMaps) { const AxisSpec axDEtaFine{640, -1.6f, 1.6f, "#Delta#eta^{true}"}; const AxisSpec axDPhi{144, -o2::constants::math::PI, o2::constants::math::PI, "#Delta#varphi^{true} (rad)"}; for (const auto& nm : {"before", "bothSurvive", "bothLost", "lostByPartnerFake"}) { + // the fine map lands in the same folder as hQ and hDEta of that fate registry.add((std::string("MCDedup/Pairs/") + nm + "/hMap").c_str(), "truth photon pairs", kTHnSparseF, {axDEtaFine, axDPhi, axQ}, true); } } @@ -441,11 +456,11 @@ struct PhotonConversionBuilder { } // In case override, don't proceed, please - no CCDB access required - if (d_bz_input > -990) { + if (d_bz_input > -990) { // o2-linter: disable=magic-number (override value) d_bz = d_bz_input; o2::parameters::GRPMagField grpmag; - if (std::fabs(d_bz) > 1e-5) { - grpmag.setL3Current(30000.f / (d_bz / 5.0f)); + if (std::fabs(d_bz) > 1e-5) { // o2-linter: disable=magic-number (override value) + grpmag.setL3Current(30000.f / (d_bz / 5.0f)); // o2-linter: disable=magic-number (override value) } o2::base::Propagator::initFieldFromGRP(&grpmag); mRunNumber = bc.runNumber(); @@ -475,7 +490,7 @@ struct PhotonConversionBuilder { } mRunNumber = bc.runNumber(); - if (useMatCorrType == 2) { + if (useMatCorrType == 2) { // o2-linter: disable=magic-number (material budget correction) // setMatLUT only after magfield has been initalized (setMatLUT has implicit and problematic init field call if not) o2::base::Propagator::Instance()->setMatLUT(lut); } @@ -735,10 +750,10 @@ struct PhotonConversionBuilder { auto phiHelix = RecoDecay::constrainAngle(std::atan2(diffY, diffX) - o2::constants::math::PI / 2.); // Electron - float arcLenghtEle = helixPosEle.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosEle.rC) * helixPosEle.rC : o2::constants::math::PI / 2.2 * helixPosEle.rC; // This assumes that the photon momentum vector is a tangent of the circle + float arcLenghtEle = helixPosEle.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosEle.rC) * helixPosEle.rC : o2::constants::math::PI / 2.2 * helixPosEle.rC; // This assumes that the photon momentum vector is a tangent of the circle // o2-linter: disable=magic-number (geometrical assumption for propagation) auto propTrackEle = getPropMomentumFromTrackHelix(arcLenghtEle, ele, helixPosEle, d_bz / 10., phiHelix - ele.phi()); // Positron - float arcLenghtPos = helixPosPos.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosPos.rC) * helixPosPos.rC : o2::constants::math::PI / 2.2 * helixPosPos.rC; // This assumes that the photon momentum vector is a tangent of the circle + float arcLenghtPos = helixPosPos.rC * 0.9 > propV0LegsRadius ? std::asin(propV0LegsRadius / helixPosPos.rC) * helixPosPos.rC : o2::constants::math::PI / 2.2 * helixPosPos.rC; // This assumes that the photon momentum vector is a tangent of the circle // o2-linter: disable=magic-number (geometrical assumption for propagation) auto propTrackPos = getPropMomentumFromTrackHelix(arcLenghtPos, pos, helixPosPos, d_bz / 10., phiHelix - pos.phi()); phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(propTrackPos[0], propTrackPos[1], propTrackPos[2], propTrackEle[0], propTrackEle[1], propTrackEle[2], pos.sign(), ele.sign(), d_bz); @@ -782,7 +797,7 @@ struct PhotonConversionBuilder { } } - if (phiv == 999.f || psipair == 999.f) { + if (phiv == 999.f || psipair == 999.f) { // o2-linter: disable=magic-number (nonsensical default value) LOG(debug) << "Propagation failed for all radii (" << propV0LegsRadius << ", 30, 10 cm). Using default values for phiv and psipair (999.f)."; } @@ -795,7 +810,7 @@ struct PhotonConversionBuilder { KFParticle gammaKF; gammaKF.SetConstructMethod(2); gammaKF.Construct(GammaDaughters.data(), 2); - if (kfMassConstrain > -0.1) { + if (kfMassConstrain > -0.1) { // o2-linter: disable=magic-number (nonsensical default value) gammaKF.SetNonlinearMassConstraint(kfMassConstrain); } KFPVertex kfpVertex = createKFPVertexFromCollision(collision); @@ -931,7 +946,7 @@ struct PhotonConversionBuilder { return; } - if (v0photoncandidate.getChi2NDF() > 6e+3) { // protection for uint16. + if (v0photoncandidate.getChi2NDF() > 6e+3) { // protection for uint16. // o2-linter: disable=magic-number (protection for uint16) return; } @@ -950,7 +965,8 @@ struct PhotonConversionBuilder { const auto meeSV = static_cast((vposDedup + veleDedup).M()); const auto score = getScoreV0(v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), deduplicationScoreWeight); - V0CandidateHelper v0Helper(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), score, meeSV); + V0CandidateHelper v0Helper(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), + v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), score, meeSV); vecV0Dedup.emplace_back(v0Helper); if (applyPCMMl) { @@ -1152,7 +1168,8 @@ struct PhotonConversionBuilder { // LOGF(info, "pca_map.size() = %d", pca_map.size()); } else if (deduplicationMode == V0DeduplicationMode::KeepAll) { - // Every candidate that survived the quality cuts is stored + // Keep-all: Every candidate that survived the quality cuts is stored, INCLUDING the + // collision duplicates. stored_v0Ids.clear(); stored_fullv0Ids.clear(); nv0_map.clear(); @@ -1173,38 +1190,39 @@ struct PhotonConversionBuilder { nv0_map.clear(); int nGreedyFallback = 0; - // Groups are build of V0s which share/are connected through a shared daughter track. + // ---- conflict groups: candidates connected through a shared leg ----- std::vector parent(nCand); for (int i = 0; i < nCand; ++i) { parent[i] = i; } auto findRoot = [&parent](int x) { while (parent[x] != x) { - parent[x] = parent[parent[x]]; + parent[x] = parent[parent[x]]; // path halving x = parent[x]; } return x; }; - // Remember the first candidate using each track. If another candidate uses the same track, merge their groups. std::unordered_map firstCandOfTrack; for (int i = 0; i < nCand; ++i) { - for (const int trackId : {vecV0Dedup[i].posID, vecV0Dedup[i].eleID}) { + for (const int& trackId : {vecV0Dedup[i].posID, vecV0Dedup[i].eleID}) { auto [it, inserted] = firstCandOfTrack.try_emplace(trackId, i); if (!inserted) { parent[findRoot(i)] = findRoot(it->second); } } } - // Collect all candidates belonging to the same group (i.e. sharing a daughter track) and sort them by score. Then select the best subset of candidates that do not share any daughter tracks. std::map> groups; for (int i = 0; i < nCand; ++i) { groups[findRoot(i)].push_back(i); } - for (auto& [root, members] : groups) { + // ---- decide each group on its own ---------------------------------- + for (auto& [root, members] : groups) { // o2-linter: disable=const-ref-in-for-loop (members is sorted in place below) const int nMembers = static_cast(members.size()); std::vector selected(nMembers, 0); - // Process candidates from best to worst score. + + // Fix the order inside the group: best score first, ties by v0ID. This makes the + // result independent of the order in which the candidates ended up in vecV0Dedup. std::sort(members.begin(), members.end(), [this](int a, int b) { if (vecV0Dedup[a].score != vecV0Dedup[b].score) { return vecV0Dedup[a].score < vecV0Dedup[b].score; @@ -1213,14 +1231,14 @@ struct PhotonConversionBuilder { }); if (nMembers == 1) { - // A single candidate has no conflict and can always be kept. - selected[0] = 1; + selected[0] = 1; // no conflict -> always kept } else if (nMembers <= maxEnum) { - // For small groups, enumerate all subsets and find the optimal set of candidates that do not share any daughter tracks. + // Exact. Instead of comparing track IDs for every subset, each leg of the group gets + // a local bit number; "share a leg" is then a single bitwise AND. std::unordered_map legBit; std::vector legMask(nMembers, 0); for (int k = 0; k < nMembers; ++k) { - for (const int trackId : {vecV0Dedup[members[k]].posID, vecV0Dedup[members[k]].eleID}) { + for (const int& trackId : {vecV0Dedup[members[k]].posID, vecV0Dedup[members[k]].eleID}) { const int bit = legBit.try_emplace(trackId, static_cast(legBit.size())).first->second; legMask[k] |= (1ull << bit); // at most 2 * 16 = 32 legs, fits into 64 bit } @@ -1228,6 +1246,9 @@ struct PhotonConversionBuilder { uint32_t bestMask = 0; int bestCount = 0; float bestScore = 0.f; + // per candidate: the best solution that CONTAINS it. Comparing that against the overall + // best tells afterwards whether a rejected candidate lost on cardinality (structural, + // no algorithm can help) or on total score in a tie (a ranking problem, i.e. fixable). std::vector bestWithCount(nMembers, -1); std::vector bestWithScore(nMembers, 0.f); for (uint32_t mask = 1; mask < (1u << nMembers); ++mask) { @@ -1250,6 +1271,9 @@ struct PhotonConversionBuilder { if (!disjoint) { continue; } + // lexicographic: cardinality first, then total score. Strict comparisons mean the + // lowest mask wins a true tie, and since members is sorted by score those are the + // better candidates. if (count > bestCount || (count == bestCount && scoreSum < bestScore)) { bestCount = count; bestScore = scoreSum; @@ -1277,7 +1301,8 @@ struct PhotonConversionBuilder { } } } else { - + // Too large to enumerate: greedy by score (members is already sorted that way). This + // gives a maximal, not necessarily a maximum matching, so it is counted below. ++nGreedyFallback; std::vector usedTracks; usedTracks.reserve(2 * nMembers); @@ -1293,6 +1318,7 @@ struct PhotonConversionBuilder { } } + // ---- store the winners, book the blockers of the losers ---------- for (int k = 0; k < nMembers; ++k) { const auto& cand = vecV0Dedup[members[k]]; if (selected[k] != 0) { @@ -1304,6 +1330,9 @@ struct PhotonConversionBuilder { if (diag == nullptr) { continue; } + // Record ALL winners that took a leg of this candidate. A cross-leg fake is typically + // blocked by TWO different true photons (one per leg); storing only the first would + // make the later MC cause-of-loss attribution depend on the loop order. const CandKey victimKey = std::make_tuple(static_cast(cand.v0ID), static_cast(cand.colID), static_cast(cand.posID), static_cast(cand.eleID)); for (int j = 0; j < nMembers; ++j) { @@ -1323,6 +1352,7 @@ struct PhotonConversionBuilder { LOGF(warning, "group matching: %d conflict group(s) larger than dedupMaxGroupSize = %d were solved greedily", nGreedyFallback, maxEnum); } + // the table expects candidates ordered by collision std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), [](const auto& a, const auto& b) { return std::get<1>(a) < std::get<1>(b); @@ -1334,33 +1364,29 @@ struct PhotonConversionBuilder { // container to keep track of which electron and positron tracks have already been used std::vector usedLegs(tracks.size(), false); std::unordered_map ownerOfLeg; // diagnostics only - // clear output containers stored_v0Ids.clear(); stored_fullv0Ids.clear(); nv0_map.clear(); - // Loop over all v0s, starting with the one with the best score - // If a v0 contains a track that is already in another (better) V0 candidate, skip V0 - for (const auto& v0Cand : vecV0Dedup) - // Skip if one of the tracks is already used - { + for (const auto& v0Cand : vecV0Dedup) { const CandKey candKey = std::make_tuple(static_cast(v0Cand.v0ID), static_cast(v0Cand.colID), static_cast(v0Cand.posID), static_cast(v0Cand.eleID)); if (usedLegs[v0Cand.posID] || usedLegs[v0Cand.eleID]) { if (diag != nullptr) { + // both legs, not just the first one found: a cross-leg fake is usually blocked by + // two different winners, and that is exactly the case being measured here auto& blockers = diag->blockersByKey[candKey]; - for (const int legId : {v0Cand.posID, v0Cand.eleID}) { + for (const int& legId : {v0Cand.posID, v0Cand.eleID}) { const auto itOwner = ownerOfLeg.find(legId); if (itOwner == ownerOfLeg.end()) { continue; } if (std::find(blockers.begin(), blockers.end(), itOwner->second) == blockers.end()) { - blockers.push_back(itOwner->second); + blockers.push_back(itOwner->second); // the same winner may hold both legs } } } continue; } - // Accept candidate usedLegs[v0Cand.posID] = true; usedLegs[v0Cand.eleID] = true; if (diag != nullptr) { @@ -1395,8 +1421,8 @@ struct PhotonConversionBuilder { auto v0 = v0s.rawIteratorAt(v0Id); if constexpr (enableFilter) { - auto collision_tmp = v0.template collision_as(); // collision where this v0 belongs. - if (!(collision_tmp.neeuls() >= 1 || collision_tmp.neeuls() + nv0_map[collision_tmp.globalIndex()] >= 2)) { + auto collision_tmp = v0.template collision_as(); // collision where this v0 belongs. + if (!(collision_tmp.neeuls() >= 1 || collision_tmp.neeuls() + nv0_map[collision_tmp.globalIndex()] >= 2)) { // o2-linter: disable=magic-number (nonsensical default value) continue; } // LOGF(info, "collision_tmp.globalIndex() = %d, collision_tmp.neeuls() = %d, nv0_map = %d", collision_tmp.globalIndex(), collision_tmp.neeuls(), nv0_map[collision_tmp.globalIndex()]); @@ -1424,16 +1450,17 @@ struct PhotonConversionBuilder { vecV0Dedup.clear(); } // end of build - //! type of V0. 0: built solely for cascades (does not pass standard V0 cuts), 1: standard 2, 3: photon-like with TPC-only use. Regular analysis should always use type 1 or 3. Filter v0Filter = o2::aod::v0::v0Type > (uint8_t)0; using FilteredV0s = soa::Filtered; - // Pair-level outcome categories. The outcome (survival or loss) and cause are listed. enum PairFate { kBefore = 0, kBothSurvive, kOneLost, kBothLost, + // NB: the four "lostBy..." entries are diagnostic selections, NOT an exclusive decomposition. + // kLostByPartnerFake is a SUBSET of kLostByCrossFake, and cross + other + true do not add up + // to oneLost + bothLost (a pair can lose two photons to two different blockers). kLostByCrossFake, kLostByPartnerFake, kLostByOtherFake, @@ -1529,9 +1556,12 @@ struct PhotonConversionBuilder { } } + // ---- per true photon: which photon passed and which one was removed + // IMPORTANT: a photon is lost only when ALL of its candidates died, so all of them are asked + // for blockers, not just the best one - different candidates have different opponents. struct PhotonFate { bool alive = false; - int best = -1; // best-score candidate + int best = -1; // best-score candidate, only used for the margin plots std::vector cands; // all true candidates of this photon std::set aliveColls; // collisions with a surviving candidate }; @@ -1550,6 +1580,9 @@ struct PhotonConversionBuilder { f.best = i; } } + + // motherId -> indices of ALL candidates that took a leg from this photon. Empty means the + // photon disappeared without anybody claiming a leg: then the mode cut, it did not dedup. std::map> blockersOfPhoton; for (const auto& [mid, f] : fate) { if (f.alive) { @@ -1557,7 +1590,7 @@ struct PhotonConversionBuilder { continue; } auto& blockers = blockersOfPhoton[mid]; - for (const int ci : f.cands) { + for (const int& ci : f.cands) { const auto itBlocker = diag.blockersByKey.find(key[ci]); if (itBlocker == diag.blockersByKey.end()) { continue; @@ -1572,10 +1605,11 @@ struct PhotonConversionBuilder { } } } - + // Why was it not in the winning set? Take the most favourable of its candidates: the one + // that could have been kept at the smallest cost in cardinality. int bestDeficit = -1; float bestExcess = 0.f; - for (const int ci : f.cands) { + for (const int& ci : f.cands) { const auto itMargin = diag.lossMargin.find(key[ci]); if (itMargin == diag.lossMargin.end()) { continue; @@ -1590,16 +1624,19 @@ struct PhotonConversionBuilder { registry.fill(HIST("MCDedup/Photons/hLostMargin"), static_cast(bestDeficit), bestExcess); } + // ONE entry per photon, so the photon-level survival rate is comparable between modes. + // hVictimBlockerClass below counts one entry per BLOCKER and a victim can have two. registry.fill(HIST("MCDedup/Photons/hFate"), blockers.empty() ? 2.f : 1.f); if (blockers.empty()) { registry.fill(HIST("MCDedup/Photons/hBlockerClass"), 3.f); // o2-linter: disable=magic-number (bin 3 = no blocker) continue; } - for (const int w : blockers) { + for (const int& w : blockers) { registry.fill(HIST("MCDedup/Photons/hBlockerClass"), static_cast(cls[w])); } - + // margin and PCA comparison relate the best candidate to its blocker - that is only the + // quantity that decided in mode 0/1. Mode 2 compares sets, so it is skipped there. if (f.best >= 0 && deduplicationMode.value < static_cast(V0DeduplicationMode::GroupMatching)) { const int w = blockers.front(); registry.fill(HIST("MCDedup/Photons/hKillMargin"), diag.snapshot[f.best].score - diag.snapshot[w].score); @@ -1607,11 +1644,12 @@ struct PhotonConversionBuilder { } } + // ---- truth pairs, grouped by MC collision ---------------------------- std::map> mothersByMcColl; for (const auto& [mid, f] : fate) { mothersByMcColl[mcparticles.iteratorAt(mid).mcCollisionId()].push_back(mid); } - constexpr float kMaxQTruth = 0.3f; + constexpr float MaxQTruth = 0.3f; for (const auto& [mcCol, mids] : mothersByMcColl) { for (size_t a = 0; a < mids.size(); ++a) { for (size_t b = a + 1; b < mids.size(); ++b) { @@ -1620,7 +1658,7 @@ struct PhotonConversionBuilder { const float e1 = std::hypot(g1.px(), g1.py(), g1.pz()); const float e2 = std::hypot(g2.px(), g2.py(), g2.pz()); const float qTrue = std::sqrt(std::max(0.f, 2.f * (e1 * e2 - g1.px() * g2.px() - g1.py() * g2.py() - g1.pz() * g2.pz()))); - if (qTrue > kMaxQTruth) { + if (qTrue > MaxQTruth) { continue; } const float dEta = g1.eta() - g2.eta(); @@ -1650,6 +1688,8 @@ struct PhotonConversionBuilder { fillTruthPairMap(kBothLost, dEta, dPhi, qTrue); } + // Who killed them? These four flags are DIAGNOSES, not a decomposition: partnerFake is + // a subset of crossFake, and one pair can end up in several categories at once. bool byCross = false, byPartner = false, byOther = false, byTrue = false; auto attribute = [&](int64_t victim, int64_t partner) { if (fate.at(victim).alive) { @@ -1659,16 +1699,18 @@ struct PhotonConversionBuilder { if (itB == blockersOfPhoton.end()) { return; } - for (const int w : itB->second) { + for (const int& w : itB->second) { if (cls[w] == kV0CrossLegFake) { byCross = true; + // partner fake = a cross-leg fake using a leg of the OTHER photon of this very + // pair -> the double-kill mechanism if (motherPos[w] == partner || motherEle[w] == partner) { byPartner = true; } } else if (cls[w] == kV0OtherFake) { byOther = true; } else if (cls[w] == kV0True) { - byTrue = true; + byTrue = true; // real photon against real photon: not repairable by any algorithm } } }; From 82953138a4653d8d8eaf2c70bf1c13c58e337324 Mon Sep 17 00:00:00 2001 From: Stefanie Mrozinski Date: Tue, 11 Aug 2026 14:39:40 +0200 Subject: [PATCH 3/3] fix requested changes + missing header for code-checker --- PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx | 6 ++++-- PWGEM/PhotonMeson/Utils/PCMUtilities.h | 4 ++-- 2 files changed, 6 insertions(+), 4 deletions(-) diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index 9462c9ddf07..7535143e3d4 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx @@ -47,6 +47,7 @@ #include #include #include +#include #include #include #include @@ -68,6 +69,7 @@ #include #include #include +#include #include #include #include @@ -169,7 +171,7 @@ struct PhotonConversionBuilder { Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; Configurable useMatCorrType{"useMatCorrType", 0, "0: none, 1: TGeo, 2: LUT"}; Configurable modeTrackPropagation{"modeTrackPropagation", 0, "0: use real track propagation, including material, 1: use fast approximation using only geometry, 2: Use real track propagation and make comparison to fast propagation (only for debugging and testing)"}; - Configurable deduplicationMode{"deduplicationMode", 0, "0: Pairwise deduplication, 1: Based on Greedy matching (best score wins), 2: Based on Group matching (crossed pairs are both kept), 3: Keep all V0s"}; + Configurable deduplicationMode{"deduplicationMode", 0, "0: Pairwise deduplication, 1: Based on Greedy matching (best score wins), 2: Based on Group matching (crossed pairs are both kept), 3: Keep all V0s, our default in the config is mode 0, however if a wrong configuration is used, the default will be mode 1 (Greedy matching) to avoid crashes"}; Configurable deduplicationScoreWeight{"deduplicationScoreWeight", 0.5f, "0.:only pca goes into the score, 1: only cosPA goes itno score, any number in between is a mix of pca and cosPA"}; Configurable cfgDedupTruthMaps{"cfgDedupTruthMaps", true, "fill the fine (dEta, dPhi, q) truth maps in MC"}; Configurable dedupMaxGroupSize{"dedupMaxGroupSize", 12, "group matching (mode 2): conflict groups up to this size are solved exactly, larger ones greedily (capped at 16)"}; @@ -1530,7 +1532,7 @@ struct PhotonConversionBuilder { } } - template + template void fillDedupTruthDiagnostics(TTracks const& tracks, TMCParticles const& mcparticles, DedupDiag const& diag) { const int nCand = static_cast(diag.snapshot.size()); diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index f7453543576..29000d78de0 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -226,7 +226,7 @@ enum V0TruthClass { /// \param track track with an MC label /// \param mcparticles the full McParticles table /// \return index of the mother if it is a photon, -1 otherwise -template +template inline int64_t getPhotonMotherId(TTrack const& track, TMCParticles const& mcparticles) { if (!track.has_mcParticle()) { @@ -248,7 +248,7 @@ inline int64_t getPhotonMotherId(TTrack const& track, TMCParticles const& mcpart /// \param motherPos photon mother of the positive leg, -1 if none /// \param motherEle photon mother of the negative leg, -1 if none /// \return kV0True, kV0CrossLegFake or kV0OtherFake -template +template inline V0TruthClass classifyV0Truth(TTrack const& pos, TTrack const& ele, TMCParticles const& mcparticles, int64_t& motherPos, int64_t& motherEle) {