diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.cxx index f91611d3c93..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,13 +69,13 @@ #include #include #include +#include #include #include #include #include #include #include -#include #include #include @@ -106,8 +107,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 +122,32 @@ struct V0CandidateHelper { float cosPA = -1.f; float pca = -1.f; float score = -1.f; + float mee = 0.f; // e+e- mass at the secondary vertex (GeV/c^2) 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; } +}; - bool operator>(const V0CandidateHelper& other) const - { - return score > other.score; - } +// (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) +using CandKey = std::tuple; + +// Everything the truth diagnosis need, +struct DedupDiag { + 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; }; struct PhotonConversionBuilder { @@ -164,8 +171,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, 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)"}; // single track cuts Configurable min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"}; @@ -286,6 +295,13 @@ struct PhotonConversionBuilder { maxSnp = 0.85f; // could be changed later maxStep = 2.00f; // could be changed later + 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, DedupNames[deduplicationMode.value], deduplicationScoreWeight.value); + ccdb->setURL(ccdburl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); @@ -393,6 +409,46 @@ 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"}; + // 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); + } + 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); + // 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); + // 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); + } + } + } } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -402,11 +458,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(); @@ -436,7 +492,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); } @@ -696,10 +752,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); @@ -743,7 +799,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)."; } @@ -756,7 +812,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); @@ -892,7 +948,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; } @@ -906,8 +962,13 @@ 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 +1074,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 +1137,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 +1168,235 @@ struct PhotonConversionBuilder { } } // end of pca_map loop // LOGF(info, "pca_map.size() = %d", pca_map.size()); + + } else if (deduplicationMode == V0DeduplicationMode::KeepAll) { + // Keep-all: Every candidate that survived the quality cuts is stored, INCLUDING the + // collision duplicates. + 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; + + // ---- 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]]; // path halving + x = parent[x]; + } + return x; + }; + 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); + } + } + } + std::map> groups; + for (int i = 0; i < nCand; ++i) { + groups[findRoot(i)].push_back(i); + } + + // ---- 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); + + // 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; + } + return vecV0Dedup[a].v0ID < vecV0Dedup[b].v0ID; + }); + + if (nMembers == 1) { + selected[0] = 1; // no conflict -> always kept + } else if (nMembers <= maxEnum) { + // 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}) { + 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; + // 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) { + 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; + } + // 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; + 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 { + // 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); + 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; + } + } + + // ---- 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) { + 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; + } + // 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) { + 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); + } + + // 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); + }); + } 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); - - // clear output containers + std::unordered_map ownerOfLeg; // diagnostics only 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 + 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}) { + 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); // the same winner may hold both legs + } + } + } 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 +1409,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); @@ -1156,8 +1423,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()]); @@ -1185,10 +1452,290 @@ 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; + 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, + 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); + } + } + + // ---- 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, only used for the margin plots + 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; + } + } + + // 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) { + 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); + } + } + } + // 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) { + 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); + } + + // 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) { + 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); + registry.fill(HIST("MCDedup/Photons/hVictimVsWinnerPCA"), diag.snapshot[f.best].pca, diag.snapshot[w].pca); + } + } + + // ---- 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 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) { + 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 > MaxQTruth) { + 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); + } + + // 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) { + return; + } + const auto itB = blockersOfPhoton.find(victim); + if (itB == blockersOfPhoton.end()) { + return; + } + 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; // real photon against real photon: not repairable by any algorithm + } + } + }; + 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 +1748,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..29000d78de0 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_