diff --git a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h index 483f38098b7..70ed1809bb0 100644 --- a/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h +++ b/PWGEM/PhotonMeson/TableProducer/photonconversionbuilder.h @@ -26,8 +26,10 @@ #include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" #include "PWGEM/PhotonMeson/DataModel/EventTables.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/ConversionPointProfile.h" #include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" #include "PWGEM/PhotonMeson/Utils/TrackSelection.h" +#include "PWGEM/PhotonMeson/Utils/V0DedupUtilities.h" #include "Common/Core/RecoDecay.h" #include "Common/Core/TPCVDriftManager.h" @@ -39,9 +41,11 @@ #include #include +#include #include #include #include +#include #include #include #include @@ -52,6 +56,7 @@ #include #include #include +#include #include #include #include @@ -74,7 +79,9 @@ #include #include #include +#include #include +#include #include #include #include @@ -90,9 +97,20 @@ using o2::framework::HistType; using o2::framework::kTH1D; using o2::framework::kTH1F; using o2::framework::kTH2F; +using o2::framework::kTH3F; using o2::framework::kTHnSparseF; using o2::framework::Preslice; using o2::framework::Produces; +using o2::pwgem::photonmeson::cpRadiusProfile; +using o2::pwgem::photonmeson::cpScanLeg; +using o2::pwgem::photonmeson::FinderCollision; +using o2::pwgem::photonmeson::FinderCuts; +using o2::pwgem::photonmeson::FinderLeg; +using o2::pwgem::photonmeson::FinderStep; +using o2::pwgem::photonmeson::LegScan; +using o2::pwgem::photonmeson::qInv; +using o2::pwgem::photonmeson::RadiusProfile; +using o2::pwgem::photonmeson::solveConflictGroup; // NOLINTEND(google-global-names-in-headers) using MyCollisions = o2::soa::Join; @@ -135,6 +153,12 @@ struct V0CandidateHelper { float pxAtPV = 0.f; float pyAtPV = 0.f; float pzAtPV = 0.f; + float svX = 0.f; + float svY = 0.f; + float svZ = 0.f; + float purity = -999.f; + float cosPAxy = -2.f; // transverse pointing angle + float cosPA3D = -2.f; // 3D pointing angle V0CandidateHelper() = default; @@ -149,14 +173,42 @@ struct V0CandidateHelper { }; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) -using CandKey = std::tuple; +struct CandKey { + int64_t v0Id = -1; // SVertexer V0 index, -1 for finder and resolver photons + int64_t collisionId = -1; + int64_t posId = -1; + int64_t eleId = -1; + auto operator<=>(const CandKey&) const = default; +}; + +// result of a probe fit (FitOptions::probe): leg pair passed the cuts, and its dedup score +struct PairFit { + bool ok = false; + float score = 999.f; +}; + +// options for fitting a leg pair that was not built by the S-Vertexer (finder, ambiguity resolver) +struct FitOptions { + bool nonSVertexer = false; // finder PCA cut for mixed legs, side-based TPC drift correction + bool relaxAP = false; // no Armenteros-Podolanski cut + bool purityGate = false; // finder: require purity score > finderMinScore + PairFit* probe = nullptr; // only return the dedup score, do not add a candidate +}; + +// one tested pair of stored photons in the ambiguity resolver (also read by the MC diagnostics) +struct SwapRecord { + CandKey keyI{}, keyJ{}; + float q = -1.f; + bool complete = false; + bool swapped = false; +}; -// 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 + std::map> blockersByKey; // rejected candidate -> candidates that took one of its legs + std::vector snapshot; // all candidates, before deduplication + std::vector stored; // surviving candidates std::map> lossMargin; + std::vector swaps; }; template @@ -197,8 +249,36 @@ struct PhotonConversionBuilder { 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)"}; + Configurable groupMatchingScoreThreshold{"groupMatchingScoreThreshold", -1.f, "dedup mode 2, choice of the best leg-disjoint subset of a conflict group. <= 0 (master): most photons first, then lowest summed score. > 0: largest sum of max(threshold - score, 1e-3) over the photons, so a photon with a score near or above the threshold adds almost nothing and two poor candidates cannot outweigh one good one (e.g. 0.4)"}; + Configurable gapWeight{"gapWeight", 0.f, "dedup score: weight of the helix-circle gap of the legs (0 = off)"}; + Configurable gapScale{"gapScale", 1.f, "dedup score: scale of the circle gap (cm)"}; + Configurable useCosPA3D{"useCosPA3D", false, "use the 3D cosPA instead of the xy cosPA in V0PhotonCandidate"}; + Configurable resolveSplitLegs{"resolveSplitLegs", false, "after the dedup, remove photons that share a split track with a better photon - avoids duplicates for photon candidates"}; + Configurable splitLegMaxDTgl{"splitLegMaxDTgl", 0.03f, "split legs: max |tgl1 - tgl2|"}; + Configurable splitLegMaxChi2{"splitLegMaxChi2", 15.f, "split legs: max chi2 of two track pieces"}; + Configurable cfgExtendedDiag{"cfgExtendedDiag", false, "MC: truth of the stored photons, pairs and swaps"}; } dedup; + // recovery finder for not thinned data (step 2) and ambiguity resolver settings (step 5) + struct : ConfigurableGroup { + std::string prefix = "recovery_group"; + Configurable runFinder{"runFinder", false, "photon finder for free legs without ITS, candidates enter the dedup"}; + Configurable finderMinScore{"finderMinScore", 0.f, "finder: min. purity score (log-likelihood ratio true / fake)"}; + Configurable finderMaxDeltaTgl{"finderMaxDeltaTgl", 0.15f, "finder: max |tgl+ - tgl-|"}; + Configurable finderMaxCircleGap{"finderMaxCircleGap", 5.f, "finder: max |circle gap| (cm)"}; + Configurable finderMinCosXY{"finderMinCosXY", 0.9f, "finder: min. xy pointing at the touching point"}; + Configurable finderMaxDeltaZ{"finderMaxDeltaZ", 3.f, "finder: max |z+ - z-| at the touching point (cm)"}; + Configurable finderMaxDz{"finderMaxDz", 4.f, "finder: max |z of the photon line at the beam axis - z_PV| (cm)"}; + Configurable finderRMargin{"finderRMargin", 5.f, "finder: margin below the min. conversion radius (cm)"}; + Configurable finderTimeMargin{"finderTimeMargin", 100.f, "finder: margin of the leg time windows (ns)"}; + Configurable finderTimeMarginFwd{"finderTimeMarginFwd", 3000.f, "finder: extra late margin for TPC-only legs (ns)"}; + Configurable finderMaxPCAMixed{"finderMaxPCAMixed", -1.f, "finder: max PCA for an ITS and a non-ITS leg (cm), < 0 builder cuts"}; + Configurable finderFreeCollinearLegs{"finderFreeCollinearLegs", 0.f, "finder: two candidates of one collision with q_inv below this (GeV/c) count as collinear; their legs stay available to the finder (leg pairs already built are not refitted). Recovers true photons next to a close partner photon. 0 = off"}; + Configurable runSwapResolver{"runSwapResolver", false, "ambiguity resolver: for close stored photon pairs, fit the crossed leg pairing (no AP cut) and keep the pairing with the lower dedup score"}; + Configurable swapMaxQ{"swapMaxQ", 0.1f, "swap resolver: max q_inv of the two photons (GeV/c)"}; + Configurable swapMaxConvDist{"swapMaxConvDist", 20.f, "swap resolver: max distance of the conversion points (cm)"}; + } recovery; + // event-level skim for derived data struct : ConfigurableGroup { std::string prefix = "skim_group"; @@ -238,6 +318,7 @@ struct PhotonConversionBuilder { Configurable max_dcav0dau_tpc_inner_fc{"max_dcav0dau_tpc_inner_fc", 1.5, "max distance btween 2 legs to V0s with ITS hits on TPC inner FC"}; Configurable min_v0radius{"min_v0radius", 1.0, "min v0 radius"}; Configurable max_v0radius{"max_v0radius", 90.0, "max v0 radius"}; + Configurable conversionPointMethod{"conversionPointMethod", 0, "0: helix touching point of the legs (Run 2), 1: radius-profiled pair fit through the material, falls back to 0 without a valid minimum"}; Configurable margin_r_its{"margin_r_its", 3.0, "margin for r cut in cm"}; Configurable margin_r_tpc{"margin_r_tpc", 7.0, "margin for r cut in cm"}; Configurable margin_r_itstpc_tpc{"margin_r_itstpc_tpc", 7.0, "margin for r cut in cm"}; @@ -289,6 +370,8 @@ struct PhotonConversionBuilder { o2::base::MatLayerCylSet* lut = nullptr; o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; o2::aod::common::TPCVDriftManager mVDriftMgr; + uint64_t mLastCCDBTimestamp{0}; + std::set mRelaxedKeys; o2::framework::HistogramRegistry registry{ "registry", @@ -448,12 +531,32 @@ struct PhotonConversionBuilder { LOG(warning) << "deduplicationScoreWeight is smaller than zero which is not allowed"; } + if (v0cuts.conversionPointMethod == 1) { + registry.add("V0/hConvPointProfile", "conversionPointMethod 1;0 fallback to the touching point, 1 profile minimum used;candidates", kTH1F, {{2, -0.5f, 1.5f}}); + registry.add("V0/hConvPointProfileShift", "conversionPointMethod 1;R_{profile} - R_{touching point} (cm);candidates", kTH1F, {{200, -10.f, 10.f}}); + } + if (recovery.runSwapResolver && dedup.deduplicationMode == V0DeduplicationMode::KeepAll) { + LOG(fatal) << "runSwapResolver needs deduplicationMode 0, 1 or 2"; + } + if (dedup.resolveSplitLegs) { + registry.add("Dedup/hSplitLegRemoved", "split legs: photons removed per DF;photons;DFs", kTH1F, {{100, -0.5f, 99.5f}}); + } + if (recovery.runFinder) { + registry.add("Finder/hSteps", "finder;step passed;pairs", kTH1D, {{FinderStep::kNFindStep, -0.5f, static_cast(FinderStep::kNFindStep) - 0.5f}}); + registry.add("Finder/hScore", "finder;purity score;pairs", kTH1F, {{100, -20.f, 10.f}}); + registry.add("Finder/hNAdded", "finder;photons added per DF;DFs", kTH1F, {{200, 0.f, 2000.f}}); + } + if (recovery.runSwapResolver) { + registry.add("Swap/hTested", "swap resolver;0 tested, 1 swapped;pairs", kTH1D, {{2, -0.5f, 1.5f}}); + registry.add("Swap/hQ", "swap resolver;q_{inv} (GeV/c);0 tested, 1 swapped", kTH2F, {{100, 0.f, 0.1f}, {2, -0.5f, 1.5f}}); + } + 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"}; // lostByPartnerFake is a subset of lostByCrossFake, and a pair can appear in several of them. - for (const auto& nm : {"before", "bothSurvive", "bothSurviveSameColl", "oneLost", "bothLost", + for (const auto& nm : {"before", "bothSurvive", "bothSurviveSameColl", "surviveAsCrossedPairing", "halfCrossedStored", "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); @@ -478,6 +581,17 @@ struct PhotonConversionBuilder { registry.add((std::string("MCDedup/Pairs/") + nm + "/hMap").c_str(), "truth photon pairs", kTHnSparseF, {axDEtaFine, axDPhi, axQ}, true); } } + if (dedup.cfgExtendedDiag) { + const AxisSpec axPt{50, 0.f, 5.f, "p_{T,#gamma} (GeV/c)"}; + const AxisSpec axQrec{50, 0.f, 0.1f, "q_{inv} (GeV/c)"}; + registry.add("MCDedup/Stored/hClassSource", "stored photons;class;0 SVertexer, 1 finder or swap;p_{T,#gamma} (GeV/c)", kTH3F, {axClass, {2, -0.5f, 1.5f}, axPt}); + registry.add("MCDedup/Stored/hPairClass", "stored photon pairs, same collision;0 true-true, 1 crossed doublet, 2 other;q_{inv} (GeV/c)", kTH2F, {{3, -0.5f, 2.5f}, axQrec}); + registry.add("MCDedup/Copies/hChoice", "true photons with a candidate in the dedup sample;collisions with a true candidate;0 = right coll., 1 = wrong (right available), 2 = wrong (no right copy), 3 = lost", kTH2F, {{5, 0.5f, 5.5f}, {4, -0.5f, 3.5f}}); + registry.add("MCDedup/Copies/hRule", "true photons, >= 2 copies, right copy available;0 = current dedup choice, 1 = smallest |#Deltaz| of the photon line, 2 = largest xy cos(PA) (master), 3 = largest 3D cos(PA);0 = wrong, 1 = right", kTH2F, {{4, -0.5f, 3.5f}, {2, -0.5f, 1.5f}}); + registry.add("MCDedup/Copies/hPAxy", "copies of true photons with >= 2 copies;0 = wrong coll., 1 = right coll.;transverse pointing angle (rad)", kTH2F, {{2, -0.5f, 1.5f}, {200, 0.f, 0.2f}}); + registry.add("MCDedup/Copies/hPA3D", "copies of true photons with >= 2 copies;0 = wrong coll., 1 = right coll.;3D pointing angle (rad)", kTH2F, {{2, -0.5f, 1.5f}, {200, 0.f, 0.2f}}); + registry.add("MCDedup/Swap/hCategory", "swap resolver;3 #times true stored + true crossed;0 kept, 1 swapped;q_{inv} (GeV/c)", kTH3F, {{9, -0.5f, 8.5f}, {2, -0.5f, 1.5f}, axQrec}); + } } mVDriftMgr.setUseSideBasedCorrection(useSideBasedCorrection.value); @@ -523,6 +637,7 @@ struct PhotonConversionBuilder { auto timestamp = bc.timestamp(); mVDriftMgr.update(timestamp); + mLastCCDBTimestamp = timestamp; } std::pair> its_ib_Requirement = {0, {0, 1, 2}}; // no hit on 3 ITS ib layers. @@ -649,12 +764,40 @@ struct PhotonConversionBuilder { template void fillV0Table(TV0 const& v0, const bool filltable) { - using namespace o2::pwgem::photonmeson; - // Get tracks const auto& pos = v0.template posTrack_as(); const auto& ele = v0.template negTrack_as(); const auto& collision = v0.template collision_as(); // collision where this v0 belongs to. - // LOGF(info, "v0.collisionId() = %d, pos.collisionId() = %d, ele.collisionId() = %d", v0.collisionId(), pos.collisionId(), ele.collisionId()); + fitV0Candidate(collision, pos, ele, v0.globalIndex(), filltable); + } + + // the purity score (finder LLR tables) is only needed for the finder gate and the split-leg order + bool usePurity() const { return recovery.runFinder.value || dedup.resolveSplitLegs.value; } + + template + static float purityScore(TTrack const& pos, TTrack const& ele, const float cosPA3D, const float pca, const float mee, const float gap, const float psipair) + { + return o2::pwgem::photonmeson::finderPurityScore({cosPA3D, pca, mee, gap, std::fabs(pos.tgl() - ele.tgl()), psipair, + std::max(std::fabs(pos.tpcNSigmaEl()), std::fabs(ele.tpcNSigmaEl())), + std::min(std::fabs(pos.tpcNSigmaPi()), std::fabs(ele.tpcNSigmaPi()))}); + } + + // deduplication score, lower = better: master score (cosPA, PCA) + optional circle-gap term + float dedupScore(const float cosPA, const float pca, const float gap) + { + float score = getScoreV0(cosPA, pca, dedup.deduplicationScoreWeight); + if (dedup.gapWeight > 0.f) { + score += dedup.gapWeight * std::min(std::fabs(gap) / dedup.gapScale.value, 3.f); // o2-linter: disable=magic-number (cap) + } + return score; + } + + // fit of one leg pair with the V0 cuts; filltable = false adds a dedup candidate, + // true writes the tables; opt.probe only returns the score + template + void fitV0Candidate(TCollision const& collision, TTrack const& pos, TTrack const& ele, const int64_t v0Id, const bool filltable, FitOptions const& opt = {}) + { + using namespace o2::pwgem::photonmeson; + mVDriftMgr.setUseSideBasedCorrection(opt.nonSVertexer || useSideBasedCorrection.value); if (pos.sign() * ele.sign() > 0) { // reject same sign pair return; @@ -692,7 +835,9 @@ struct PhotonConversionBuilder { auto pTrack = getTrackParCov(pos); if (trackcuts.moveTPCTracks && isTPConlyTrack(pos) && !mVDriftMgr.moveTPCTrack(collision, pos, pTrack)) { - LOGP(error, "failed correction for positive tpc track"); + if (opt.probe == nullptr) { + LOGP(error, "failed correction for positive tpc track"); + } return; } auto pTrackC = pTrack; @@ -718,7 +863,9 @@ struct PhotonConversionBuilder { auto nTrack = getTrackParCov(ele); if (trackcuts.moveTPCTracks && isTPConlyTrack(ele) && !mVDriftMgr.moveTPCTrack(collision, ele, nTrack)) { - LOGP(error, "failed correction for negative tpc track"); + if (opt.probe == nullptr) { + LOGP(error, "failed correction for negative tpc track"); + } return; } auto nTrackC = nTrack; @@ -744,6 +891,19 @@ struct PhotonConversionBuilder { std::array xyz = {0.f, 0.f, 0.f}; Vtx_recalculationParCov(o2::base::Propagator::Instance(), pTrack, nTrack, xyz, matCorr); + if (v0cuts.conversionPointMethod == 1) { + LegScan scanP, scanN; + cpScanLeg(pTrack, d_bz, matCorr, scanP); + cpScanLeg(nTrack, d_bz, matCorr, scanN); + RadiusProfile prof; + cpRadiusProfile(scanP, scanN, pos.hasITS() ? 0.f : 0.25f, ele.hasITS() ? 0.f : 0.25f, prof); // o2-linter: disable=magic-number (0.5 cm extra z for legs without ITS) + const bool useProfile = prof.iMin >= 0 && prof.status != 4; // o2-linter: disable=magic-number (profile class: failed) + if (useProfile) { + registry.fill(HIST("V0/hConvPointProfileShift"), prof.rMin - std::hypot(xyz[0], xyz[1])); + xyz = prof.xyz; + } + registry.fill(HIST("V0/hConvPointProfile"), useProfile ? 1.f : 0.f); + } float rxy_tmp = RecoDecay::sqrtSumOfSquares(xyz[0], xyz[1]); if (rxy_tmp > trackcuts.maxX + v0cuts.margin_r_tpc) { return; @@ -861,7 +1021,7 @@ struct PhotonConversionBuilder { return; } - if (!filltable) { + if (!filltable && opt.probe == nullptr) { if (isITSTPCTrack(pos) && isITSTPCTrack(ele)) { registry.fill(HIST("V0/hRxy_minX_ITSTPC_ITSTPC"), std::min(pTrack.getX(), nTrack.getX()), std::min(pTrack.getX(), nTrack.getX()) - rxy); // trackiu.x() - rxy should be positive } else if (isITSonlyTrack(pos) && isITSonlyTrack(ele)) { @@ -933,7 +1093,8 @@ struct PhotonConversionBuilder { float cospaXYKF = cospaXY_KF(gammaKF_DecayVtx, KFPV); float cospaRZKF = cospaRZ_KF(gammaKF_DecayVtx, KFPV); auto centType = static_cast(mlcuts.centTypePCMMl.value); - v0photoncandidate.setPhotonCandidate(gammaKF_DecayVtx, gammaKF_PV, pos, kfp_pos_DecayVtx, ele, kfp_ele_DecayVtx, collision, cospaXYKF, cospaRZKF, cospaXYKF, psipair, phiv, centType, posdcaXY, eledcaXY, posdcaZ, eledcaZ); + // cosPA slot of the candidate (dedup score, cosPA cut): xy as in master, or 3D with dedup_group.useCosPA3D + v0photoncandidate.setPhotonCandidate(gammaKF_DecayVtx, gammaKF_PV, pos, kfp_pos_DecayVtx, ele, kfp_ele_DecayVtx, collision, dedup.useCosPA3D ? cospa_kf : cospaXYKF, cospaRZKF, cospaXYKF, psipair, phiv, centType, posdcaXY, eledcaXY, posdcaZ, eledcaZ); if (!ele.hasITS() && !pos.hasITS()) { // V0s with TPConly-TPConly if (v0cuts.max_r_itsmft_ss < rxy && rxy < trackcuts.maxX + v0cuts.margin_r_tpc) { @@ -945,6 +1106,10 @@ struct PhotonConversionBuilder { return; } } + } else if (opt.nonSVertexer && recovery.finderMaxPCAMixed > 0.f && pos.hasITS() != ele.hasITS()) { + if (v0photoncandidate.getPCA() > recovery.finderMaxPCAMixed) { + return; + } } else { // V0s with ITS hits if (rxy < v0cuts.max_r_req_its) { if (v0photoncandidate.getPCA() > v0cuts.max_dcav0dau_itsibss) { @@ -973,22 +1138,43 @@ struct PhotonConversionBuilder { return; } - if (!checkAP(v0photoncandidate.getAlpha(), v0photoncandidate.getQt(), v0cuts.max_alpha_ap, v0cuts.max_qt_ap)) { // store only photon conversions + if (!opt.relaxAP && !checkAP(v0photoncandidate.getAlpha(), v0photoncandidate.getQt(), v0cuts.max_alpha_ap, v0cuts.max_qt_ap)) { // store only photon conversions return; } - 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(), dedup.deduplicationScoreWeight); - V0CandidateHelper v0Helper(v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), + // passed all cuts: score and dedup candidate + const float gap = circleGap(o2::track::TrackAuxPar(pTrack, d_bz), o2::track::TrackAuxPar(nTrack, d_bz)); + const float score = dedupScore(v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), gap); + if (opt.probe != nullptr) { + *opt.probe = {true, score}; + return; + } + const float purity = usePurity() ? purityScore(pos, ele, cospa_kf, v0photoncandidate.getPCA(), meeSV, gap, v0photoncandidate.getPsiPair()) : -999.f; + if (opt.purityGate) { + registry.fill(HIST("Finder/hScore"), purity); + if (purity <= recovery.finderMinScore) { + return; + } + } + + pca_map[CandKey{v0Id, collision.globalIndex(), pos.globalIndex(), ele.globalIndex()}] = v0photoncandidate.getPCA(); + cospa_map[CandKey{v0Id, collision.globalIndex(), pos.globalIndex(), ele.globalIndex()}] = v0photoncandidate.getCosPA(); + + V0CandidateHelper v0Helper(v0Id, collision.globalIndex(), pos.globalIndex(), ele.globalIndex(), v0photoncandidate.getCosPA(), v0photoncandidate.getPCA(), score, meeSV); v0Helper.pxAtPV = gammaKF_PV.GetPx(); v0Helper.pyAtPV = gammaKF_PV.GetPy(); v0Helper.pzAtPV = gammaKF_PV.GetPz(); + v0Helper.svX = gammaKF_DecayVtx.GetX(); + v0Helper.svY = gammaKF_DecayVtx.GetY(); + v0Helper.svZ = gammaKF_DecayVtx.GetZ(); + v0Helper.purity = purity; + v0Helper.cosPAxy = cospaXYKF; + v0Helper.cosPA3D = cospa_kf; vecV0Dedup.emplace_back(v0Helper); if (mlcuts.applyPCMMl) { @@ -1070,7 +1256,7 @@ struct PhotonConversionBuilder { ROOT::Math::PxPyPzMVector v0_sv = vpos_sv + vele_sv; registry.fill(HIST("V0/hMeeSV_Rxy"), rxy, v0_sv.M()); - v0photonskf(collision.globalIndex(), v0.globalIndex(), v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, + v0photonskf(collision.globalIndex(), v0Id, v0legs.lastIndex() + 1, v0legs.lastIndex() + 2, gammaKF_DecayVtx.GetX(), gammaKF_DecayVtx.GetY(), gammaKF_DecayVtx.GetZ(), gammaKF_PV.GetPx(), gammaKF_PV.GetPy(), gammaKF_PV.GetPz(), v0_sv.M(), v0photoncandidate.getDcaXYToPV(), v0photoncandidate.getDcaZToPV(), @@ -1085,338 +1271,718 @@ struct PhotonConversionBuilder { } // end of fill table } - Preslice perCollision = o2::aod::v0::collisionId; - std::map, float> pca_map; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) -> pca - std::map, float> cospa_map; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) -> cospa - std::vector vecV0Dedup; // vector with all V0Candidates that is used to sort them by score (see struct V0CandidateHelper for more details of content) - std::vector> stored_v0Ids; // (pos.globalIndex(), ele.globalIndex()) - std::vector> stored_fullv0Ids; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) - std::unordered_map nv0_map; // map collisionId -> nv0 + // bookkeeping of candidates and stored photons + static CandKey keyOf(V0CandidateHelper const& c) + { + return CandKey{c.v0ID, c.colID, c.posID, c.eleID}; + } - template - void build(TCollisions const& collisions, TV0s const& v0s, TTracks const& tracks, TBCs const&, DedupDiag* diag = nullptr) + std::map candidateIndexByKey() const { - for (const auto& collision : collisions) { + std::map index; + for (int i = 0; i < static_cast(vecV0Dedup.size()); ++i) { + index[keyOf(vecV0Dedup[i])] = i; + } + return index; + } + + // replaces the stored photons after the dedup and keeps stored_v0Ids and nv0_map consistent + void setStored(std::vector keys) + { + std::stable_sort(keys.begin(), keys.end(), [](const auto& a, const auto& b) { return a.collisionId < b.collisionId; }); + stored_fullv0Ids.swap(keys); + stored_v0Ids.clear(); + for (auto& [col, n] : nv0_map) { // o2-linter: disable=const-ref-in-for-loop (counters are reset in place) + n = 0; + } + for (const auto& key : stored_fullv0Ids) { + stored_v0Ids.emplace_back(key.posId, key.eleId); + nv0_map[key.collisionId]++; + } + } + + // photon finder + // finder leg: helix circle and the time window in which the track can belong to a collision (TPC-only legs: from the drift-time asymmetry) + template + bool makeFinderLeg(TTrack const& track, TCollisions const& collisions, const int64_t bcRef, FinderLeg& leg) + { + using namespace o2::pwgem::photonmeson; + if (track.collisionId() < 0 || !track.hasTPC()) { + return false; + } + leg.id = track.globalIndex(); + leg.tgl = track.tgl(); + leg.noITS = !track.hasITS(); + leg.movable = isTPConlyTrack(track) && (track.flags() & o2::aod::track::TrackTimeAsym) != 0; + const auto bc = collisions.rawIteratorAt(track.collisionId()).template bc_as(); + leg.t = static_cast(static_cast(bc.globalBC()) - bcRef) * o2::constants::lhc::LHCBunchSpacingNS + track.trackTime(); + if (leg.movable) { + const int side = tpcDriftSide(track); + if (side == TPCSideUndetermined) { + return false; + } + leg.side = side; + leg.tMin = leg.t - track.tpcDeltaTBwd() - recovery.finderTimeMargin; + leg.tMax = leg.t + track.tpcDeltaTFwd() + recovery.finderTimeMarginFwd; + } else { + constexpr float NSigmaTime = 4.f; + const float sigma = std::max(track.trackTimeRes(), 1.f); + const float halfWidth = (track.flags() & o2::aod::track::TrackTimeResIsRange) != 0 ? sigma : NSigmaTime * sigma; + leg.side = 0; + leg.tMin = leg.t - halfWidth - recovery.finderTimeMargin; + leg.tMax = leg.t + halfWidth + recovery.finderTimeMargin; + } + leg.par = getTrackPar(track); + leg.circle = o2::track::TrackAuxPar(leg.par, d_bz); + return true; + } + + // pairs the free legs of the whole DF (not used by any candidate unless finderFreeCollinearLegs, at least one leg without ITS); a pair that passes the geometry + // (evalFinderPair) and matches a collision in z is fitted and becomes a dedup candidate with v0Id = -1 + template + void runPhotonFinder(TCollisions const& collisions, TTracks const& tracks) + { + std::set usedLegs; + std::set> builtPairs; + // legs of accepted candidates are not free; with finderFreeCollinearLegs, legs of collinear candidate pairs stay free + std::vector collinear(vecV0Dedup.size(), false); + if (recovery.finderFreeCollinearLegs > 0.f) { + std::map> candsOfColl; + for (size_t i = 0; i < vecV0Dedup.size(); ++i) { + candsOfColl[vecV0Dedup[i].colID].push_back(i); + } + for (const auto& [colId, members] : candsOfColl) { + for (size_t a = 0; a < members.size(); ++a) { + for (size_t b = a + 1; b < members.size(); ++b) { + const auto& ci = vecV0Dedup[members[a]]; + const auto& cj = vecV0Dedup[members[b]]; + if (ci.posID == cj.posID || ci.eleID == cj.eleID) { + continue; + } + if (qInv(ci.pxAtPV, ci.pyAtPV, ci.pzAtPV, cj.pxAtPV, cj.pyAtPV, cj.pzAtPV) < recovery.finderFreeCollinearLegs) { + collinear[members[a]] = true; + collinear[members[b]] = true; + } + } + } + } + } + for (size_t i = 0; i < vecV0Dedup.size(); ++i) { + const auto& c = vecV0Dedup[i]; + builtPairs.emplace(c.posID, c.eleID); + if (!collinear[i]) { + usedLegs.insert(c.posID); + usedLegs.insert(c.eleID); + } + } + std::vector colls; + int64_t bcRef = -1; + for (const auto& col : collisions) { if constexpr (isMC) { - if (!collision.has_mcCollision()) { + if (!col.has_mcCollision()) { continue; } } - - if (!collision.isSelected()) { + if (!col.isSelected()) { continue; } + const auto bc = col.template bc_as(); + if (bcRef < 0) { + bcRef = static_cast(bc.globalBC()); + } + colls.push_back({static_cast(static_cast(bc.globalBC()) - bcRef) * o2::constants::lhc::LHCBunchSpacingNS + col.collisionTime(), col.globalIndex(), col.posX(), col.posY(), col.posZ()}); + } + if (colls.empty()) { + return; + } + std::sort(colls.begin(), colls.end(), [](const FinderCollision& a, const FinderCollision& b) { return a.t < b.t; }); + const auto* vd = ccdb->getForTimeStamp("TPC/Calib/VDriftTgl", mLastCCDBTimestamp); + if (vd == nullptr) { + LOG(warning) << "photon finder: no TPC drift velocity for timestamp " << mLastCCDBTimestamp; + return; + } + const double vDrift = vd->refVDrift * vd->corrFact * 1e-3; + float xB = 0.f, yB = 0.f; + for (const auto& c : colls) { + xB += c.x / static_cast(colls.size()); + yB += c.y / static_cast(colls.size()); + } - // if constexpr (isTriggerAnalysis) { - // if (collision.triggerMask_raw() == 0) { - // continue; - // } - // } - - nv0_map[collision.globalIndex()] = 0; - - const auto& bc = collision.template foundBC_as(); - initCCDB(bc); - registry.fill(HIST("hCollisionCounter"), 1); + std::vector posLegs, eleLegs; + for (const auto& track : tracks) { + if (usedLegs.contains(track.globalIndex()) || track.pt() < trackcuts.min_pt_trackiu || !checkV0leg(track)) { + continue; + } + FinderLeg leg; + if (makeFinderLeg(track, collisions, bcRef, leg)) { + (track.sign() > 0 ? posLegs : eleLegs).push_back(leg); + } + } + std::sort(eleLegs.begin(), eleLegs.end(), [](const FinderLeg& a, const FinderLeg& b) { return a.tgl < b.tgl; }); - updateCCDB(bc); // delay update until is needed + const FinderCuts cuts{recovery.finderMaxDeltaTgl, recovery.finderMaxCircleGap, v0cuts.min_r_tpconly, v0cuts.max_r_req_its, recovery.finderRMargin, + trackcuts.maxX + v0cuts.margin_r_tpc, recovery.finderMinCosXY, recovery.finderMaxDeltaZ, recovery.finderMaxDz, d_bz}; + const FitOptions finderFit{.nonSVertexer = true, .purityGate = true}; + const size_t nBefore = vecV0Dedup.size(); + for (const auto& p : posLegs) { + for (auto itE = std::lower_bound(eleLegs.begin(), eleLegs.end(), p.tgl - recovery.finderMaxDeltaTgl, [](const FinderLeg& c, float v) { return c.tgl < v; }); + itE != eleLegs.end() && itE->tgl <= p.tgl + recovery.finderMaxDeltaTgl; ++itE) { + if ((!p.noITS && !itE->noITS) || builtPairs.contains({p.id, itE->id})) { + continue; + } + int64_t colId = -1; + const int step = evalFinderPair(p, *itE, colls, cuts, xB, yB, vDrift, colId); + if (step >= 0) { + registry.fill(HIST("Finder/hSteps"), static_cast(step)); + } + if (step == FinderStep::kFindCollision) { + fitV0Candidate(collisions.rawIteratorAt(colId), tracks.rawIteratorAt(p.id), tracks.rawIteratorAt(itE->id), -1, false, finderFit); + } + } + } + registry.fill(HIST("Finder/hNAdded"), static_cast(vecV0Dedup.size() - nBefore)); + } - vecV0Dedup.reserve(30000); // rough estimate for number of V0s per DF - const auto& v0s_per_coll = v0s.sliceBy(perCollision, collision.globalIndex()); - // LOGF(info, "n v0 = %d", v0s_per_coll.size()); - for (const auto& v0 : v0s_per_coll) { - // LOGF(info, "collision.globalIndex() = %d, v0.globalIndex() = %d, v0.posTrackId() = %d, v0.negTrackId() = %d", collision.globalIndex(), v0.globalIndex(), v0.posTrackId() , v0.negTrackId()); - fillV0Table(v0, false); - } // end of v0 loop - } // end of collision loop + // two stored photons that use two pieces of one track (same tgl, compatible parameters): keep the better one + template + void resolveSplitLegConflicts(TTracks const& tracks, DedupDiag* diag) + { + const auto indexOfKey = candidateIndexByKey(); + std::vector> order; + order.reserve(stored_fullv0Ids.size()); + for (size_t k = 0; k < stored_fullv0Ids.size(); ++k) { + const auto it = indexOfKey.find(stored_fullv0Ids[k]); + order.emplace_back(it == indexOfKey.end() ? std::numeric_limits::max() : (usePurity() ? -vecV0Dedup[it->second].purity : vecV0Dedup[it->second].score), k); + } + std::stable_sort(order.begin(), order.end(), [](const auto& a, const auto& b) { return a.first < b.first; }); + auto samePiece = [&](int64_t idA, int64_t idB) { + if (idA == idB) { + return false; + } + const auto trA = tracks.rawIteratorAt(idA); + const auto trB = tracks.rawIteratorAt(idB); + if (std::fabs(trA.tgl() - trB.tgl()) > dedup.splitLegMaxDTgl) { + return false; + } + const float chi2 = splitPieceChi2(trA, trB, d_bz); + return chi2 >= 0.f && chi2 < dedup.splitLegMaxChi2; + }; + std::vector> keptPos, keptEle; + std::vector removed(stored_fullv0Ids.size(), 0); + int nRemoved = 0; + for (const auto& [rank, k] : order) { + const auto& key = stored_fullv0Ids[k]; + int64_t owner = -1; + for (size_t m = 0; m < keptPos.size() && owner < 0; ++m) { + if (samePiece(key.posId, keptPos[m].first)) { + owner = static_cast(keptPos[m].second); + } + } + for (size_t m = 0; m < keptEle.size() && owner < 0; ++m) { + if (samePiece(key.eleId, keptEle[m].first)) { + owner = static_cast(keptEle[m].second); + } + } + if (owner >= 0) { + removed[k] = 1; + ++nRemoved; + if (diag != nullptr) { + diag->blockersByKey[key].push_back(stored_fullv0Ids[owner]); + } + continue; + } + keptPos.emplace_back(key.posId, k); + keptEle.emplace_back(key.eleId, k); + } + if (nRemoved > 0) { + std::vector kept; + kept.reserve(stored_fullv0Ids.size() - nRemoved); + for (size_t k = 0; k < stored_fullv0Ids.size(); ++k) { + if (removed[k] == 0) { + kept.push_back(stored_fullv0Ids[k]); + } + } + setStored(std::move(kept)); + } + registry.fill(HIST("Dedup/hSplitLegRemoved"), static_cast(nRemoved)); + } - stored_v0Ids.reserve(pca_map.size()); // number of photon candidates per DF - stored_fullv0Ids.reserve(pca_map.size()); // number of photon candidates per DF + // ambiguity resolver + // two stored photons of one collision, close in q and in space: fit the crossed leg pairing; if both crossed photons pass and + // their summed score is lower, swap + template + void runSwapResolver(TCollisions const& collisions, TTracks const& tracks, DedupDiag* diag) + { + const auto indexOfKey = candidateIndexByKey(); + std::map> storedOfColl; + for (const auto& key : stored_fullv0Ids) { + if (const auto it = indexOfKey.find(key); it != indexOfKey.end()) { + storedOfColl[key.collisionId].push_back(it->second); + } + } + // the AP cut is the main reason the true pairing of a crossed group was not built + auto fitPair = [&](int64_t colId, int64_t posId, int64_t eleId, PairFit* probe) { + fitV0Candidate(collisions.rawIteratorAt(colId), tracks.rawIteratorAt(posId), tracks.rawIteratorAt(eleId), -1, false, + FitOptions{.nonSVertexer = true, .relaxAP = true, .probe = probe}); + }; + auto probeFit = [&](int64_t colId, int64_t posId, int64_t eleId) { + PairFit out; + fitPair(colId, posId, eleId, &out); + return out; + }; + auto crossedPrefilter = [&](int64_t posId, int64_t eleId) { + const auto pos = tracks.rawIteratorAt(posId); + const auto ele = tracks.rawIteratorAt(eleId); + if (std::fabs(pos.tgl() - ele.tgl()) > recovery.finderMaxDeltaTgl) { + return false; + } + const o2::track::TrackAuxPar cp(getTrackPar(pos), d_bz); + const o2::track::TrackAuxPar ce(getTrackPar(ele), d_bz); + return std::fabs(circleGap(cp, ce)) < recovery.finderMaxCircleGap; + }; + std::map storedFit; + auto storedFitOf = [&](int idx) -> const PairFit& { + auto it = storedFit.find(idx); + if (it == storedFit.end()) { + const auto& c = vecV0Dedup[idx]; + it = storedFit.emplace(idx, probeFit(c.colID, c.posID, c.eleID)).first; + } + return it->second; + }; - // find minimal pca - if (dedup.deduplicationMode == V0DeduplicationMode::Pairwise) { - for (const auto& [key, value] : pca_map) { - auto v0Id = std::get<0>(key); - auto collisionId = std::get<1>(key); - auto posId = std::get<2>(key); - auto eleId = std::get<3>(key); - float v0pca = value; - float cospa = cospa_map[key]; - bool is_closest_v0 = true; - bool is_most_aligned_v0 = true; - - for (const auto& [key_tmp, value_tmp] : pca_map) { - auto v0Id_tmp = std::get<0>(key_tmp); - auto collisionId_tmp = std::get<1>(key_tmp); - auto posId_tmp = std::get<2>(key_tmp); - auto eleId_tmp = std::get<3>(key_tmp); - float v0pca_tmp = value_tmp; - float cospa_tmp = cospa_map[key_tmp]; - - if (v0Id == v0Id_tmp) { // skip exactly the same v0 + // test all close pairs of stored photons + struct Test { + SwapRecord rec; + float gain = -999.f; + }; + std::vector swaps; + for (const auto& [colId, members] : storedOfColl) { + std::vector tests; + for (size_t a = 0; a < members.size(); ++a) { + for (size_t b = a + 1; b < members.size(); ++b) { + const auto& ci = vecV0Dedup[members[a]]; + const auto& cj = vecV0Dedup[members[b]]; + if (ci.posID == cj.posID || ci.eleID == cj.eleID) { continue; } - - 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; + const float q = qInv(ci.pxAtPV, ci.pyAtPV, ci.pzAtPV, cj.pxAtPV, cj.pyAtPV, cj.pzAtPV); + if (q > recovery.swapMaxQ || std::hypot(ci.svX - cj.svX, ci.svY - cj.svY, ci.svZ - cj.svZ) > recovery.swapMaxConvDist) { + continue; } - if ((eleId == eleId_tmp || posId == posId_tmp) && v0pca > v0pca_tmp) { - if (diag != nullptr) { - diag->blockersByKey[key].push_back(key_tmp); - } - is_closest_v0 = false; - break; + const PairFit storedI = storedFitOf(members[a]); + const PairFit storedJ = storedFitOf(members[b]); + const PairFit crossedI = crossedPrefilter(ci.posID, cj.eleID) ? probeFit(colId, ci.posID, cj.eleID) : PairFit{}; + const PairFit crossedJ = crossedPrefilter(cj.posID, ci.eleID) ? probeFit(colId, cj.posID, ci.eleID) : PairFit{}; + Test t; + t.rec = {keyOf(ci), keyOf(cj), q, storedI.ok && storedJ.ok && crossedI.ok && crossedJ.ok, false}; + if (t.rec.complete) { // gain > 0: the crossed pairing has the lower (better) summed dedup score + t.gain = (storedI.score + storedJ.score) - (crossedI.score + crossedJ.score); } - } // end of pca_map tmp loop - - bool is_stored = std::find(stored_v0Ids.begin(), stored_v0Ids.end(), std::make_pair(posId, eleId)) != stored_v0Ids.end(); - if (is_closest_v0 && is_most_aligned_v0 && !is_stored) { - // auto v0 = v0s.rawIteratorAt(v0Id); - // auto collision = collisions.rawIteratorAt(collisionId); - // auto pos = tracks.rawIteratorAt(posId); - // auto ele = tracks.rawIteratorAt(eleId); - // LOGF(info, "!accept! | collision id = %d | v0id1 = %d , posid1 = %d , eleid1 = %d , pca1 = %f , cospa = %f", collisionId, v0Id, posId, eleId, v0pca, cospa); - - // fillV0Table(v0, true); - stored_v0Ids.emplace_back(std::make_pair(posId, eleId)); - stored_fullv0Ids.emplace_back(std::make_tuple(v0Id, collisionId, posId, eleId)); - nv0_map[collisionId]++; + tests.push_back(t); } - } // end of pca_map loop - // LOGF(info, "pca_map.size() = %d", pca_map.size()); - - } else if (dedup.deduplicationMode == V0DeduplicationMode::KeepAll) { - // Keep-all: Every candidate that survived the quality cuts is stored, including all 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 (dedup.deduplicationMode == V0DeduplicationMode::GroupMatching) { - const int maxEnum = std::min(dedup.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]; + } + std::vector order(tests.size()); + std::iota(order.begin(), order.end(), 0); + std::sort(order.begin(), order.end(), [&](int x, int y) { return tests[x].gain > tests[y].gain; }); + std::set used; + for (const int& k : order) { + auto& t = tests[k]; + if (t.rec.complete && t.gain > 0.f && !used.contains(t.rec.keyI) && !used.contains(t.rec.keyJ)) { + t.rec.swapped = true; + used.insert(t.rec.keyI); + used.insert(t.rec.keyJ); + swaps.push_back(t.rec); } - 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); - } + registry.fill(HIST("Swap/hTested"), 0.f); + registry.fill(HIST("Swap/hQ"), t.rec.q, t.rec.swapped ? 1.f : 0.f); + if (diag != nullptr) { + diag->swaps.push_back(t.rec); } } - std::map> groups; - for (int i = 0; i < nCand; ++i) { - groups[findRoot(i)].push_back(i); + } + if (swaps.empty()) { + return; + } + + // apply swaps: the crossed photons become candidates and replace the stored pair + const std::set storedBefore(stored_fullv0Ids.begin(), stored_fullv0Ids.end()); + std::set removedKeys; + std::vector addedKeys; + for (const auto& s : swaps) { + const int64_t colId = s.keyI.collisionId; + const int64_t posI = s.keyI.posId, eleI = s.keyI.eleId, posJ = s.keyJ.posId, eleJ = s.keyJ.eleId; + const CandKey a1 = CandKey{-1, colId, posI, eleJ}; + const CandKey a2 = CandKey{-1, colId, posJ, eleI}; + if (storedBefore.contains(a1) || storedBefore.contains(a2)) { + continue; + } + if (!pca_map.contains(a1)) { + fitPair(colId, posI, eleJ, nullptr); } + if (!pca_map.contains(a2)) { + fitPair(colId, posJ, eleI, nullptr); + } + if (!pca_map.contains(a1) || !pca_map.contains(a2)) { + continue; + } + removedKeys.insert(s.keyI); + removedKeys.insert(s.keyJ); + addedKeys.push_back(a1); + addedKeys.push_back(a2); + mRelaxedKeys.insert(a1); // refitted for the table without the AP cut, as in the resolver fit + mRelaxedKeys.insert(a2); + registry.fill(HIST("Swap/hTested"), 1.f); + if (diag != nullptr) { + diag->blockersByKey[s.keyI] = {a1, a2}; + diag->blockersByKey[s.keyJ] = {a1, a2}; + } + } + std::vector kept; + kept.reserve(stored_fullv0Ids.size()); + for (const auto& key : stored_fullv0Ids) { + if (!removedKeys.contains(key)) { + kept.push_back(key); + } + } + kept.insert(kept.end(), addedKeys.begin(), addedKeys.end()); + setStored(std::move(kept)); + } - // ---- 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); + // deduplication + Preslice perCollision = o2::aod::v0::collisionId; + std::map pca_map; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) -> pca + std::map cospa_map; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) -> cospa + std::vector vecV0Dedup; // vector with all V0Candidates that is used to sort them by score (see struct V0CandidateHelper for more details of content) + std::vector> stored_v0Ids; // (pos.globalIndex(), ele.globalIndex()) + std::vector stored_fullv0Ids; // (v0.globalIndex(), collision.globalIndex(), pos.globalIndex(), ele.globalIndex()) + std::unordered_map nv0_map; // map collisionId -> nv0 + + // mode 0: a candidate is kept if no other candidate sharing a leg has a smaller PCA and no copy in another collision has a larger cosPA + void selectPairwise(DedupDiag* diag) + { + for (const auto& [key, value] : pca_map) { + auto v0Id = key.v0Id; + auto collisionId = key.collisionId; + auto posId = key.posId; + auto eleId = key.eleId; + float v0pca = value; + float cospa = cospa_map[key]; + bool is_closest_v0 = true; + bool is_most_aligned_v0 = true; + + for (const auto& [key_tmp, value_tmp] : pca_map) { + auto collisionId_tmp = key_tmp.collisionId; + auto posId_tmp = key_tmp.posId; + auto eleId_tmp = key_tmp.eleId; + float v0pca_tmp = value_tmp; + float cospa_tmp = cospa_map[key_tmp]; + + if (key == key_tmp) { // skip exactly the same candidate + continue; + } - // Fix the order inside the group: best score first - 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; + if (collisionId != collisionId_tmp && eleId == eleId_tmp && posId == posId_tmp && cospa < cospa_tmp) { + if (diag != nullptr) { + diag->blockersByKey[key].push_back(key_tmp); } - return vecV0Dedup[a].v0ID < vecV0Dedup[b].v0ID; - }); - - if (nMembers == 1) { - selected[0] = 1; - } else if (nMembers <= maxEnum) { - - 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 - } + is_most_aligned_v0 = false; + break; + } + if ((eleId == eleId_tmp || posId == posId_tmp) && v0pca > v0pca_tmp) { + if (diag != nullptr) { + diag->blockersByKey[key].push_back(key_tmp); } - 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; - } + is_closest_v0 = false; + break; + } + } // end of pca_map tmp loop + + bool is_stored = std::find(stored_v0Ids.begin(), stored_v0Ids.end(), std::make_pair(posId, eleId)) != stored_v0Ids.end(); + if (is_closest_v0 && is_most_aligned_v0 && !is_stored) { + // auto v0 = v0s.rawIteratorAt(v0Id); + // auto collision = collisions.rawIteratorAt(collisionId); + // auto pos = tracks.rawIteratorAt(posId); + // auto ele = tracks.rawIteratorAt(eleId); + // LOGF(info, "!accept! | collision id = %d | v0id1 = %d , posid1 = %d , eleid1 = %d , pca1 = %f , cospa = %f", collisionId, v0Id, posId, eleId, v0pca, cospa); + + // fillV0Table(v0, true); + stored_v0Ids.emplace_back(std::make_pair(posId, eleId)); + stored_fullv0Ids.push_back(CandKey{v0Id, collisionId, posId, eleId}); + nv0_map[collisionId]++; + } + } // end of pca_map loop + if (recovery.runFinder || recovery.runSwapResolver) { + std::stable_sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), [](const auto& a, const auto& b) { return a.collisionId < b.collisionId; }); + } + // LOGF(info, "pca_map.size() = %d", pca_map.size()); + } - 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); - } + // mode 3: every candidate is kept, including all duplicates + void keepAll() + { + // Keep-all: Every candidate that survived the quality cuts is stored, including all 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.push_back(keyOf(v0Cand)); + nv0_map[v0Cand.colID]++; + } + std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), + [](const auto& a, const auto& b) { + return a.collisionId < b.collisionId; + }); + } + + // mode 2: candidates sharing a leg form a conflict group; per group the best leg-disjoint subset is kept + void selectByConflictGroup(DedupDiag* diag) + { + const int maxEnum = std::min(dedup.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 + 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 + 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; + } + if (vecV0Dedup[a].v0ID != vecV0Dedup[b].v0ID) { + return vecV0Dedup[a].v0ID < vecV0Dedup[b].v0ID; + } + return std::tie(vecV0Dedup[a].posID, vecV0Dedup[a].eleID, vecV0Dedup[a].colID) < std::tie(vecV0Dedup[b].posID, vecV0Dedup[b].eleID, vecV0Dedup[b].colID); + }); + + if (nMembers == 1) { + selected[0] = 1; + } else if (nMembers <= maxEnum) { + + 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 } - } else { - // Too large to enumerate: greedy by score (members is already sorted that way) - ++nGreedyFallback; - std::vector usedTracks; - usedTracks.reserve(2 * nMembers); - for (int k = 0; k < nMembers; ++k) { + } + std::vector scores(nMembers); + for (int k = 0; k < nMembers; ++k) { + scores[k] = vecV0Dedup[members[k]].score; + } + const auto sol = solveConflictGroup(legMask, scores, dedup.groupMatchingScoreThreshold.value, diag != nullptr); + for (int k = 0; k < nMembers; ++k) { + selected[k] = ((sol.bestMask & (1u << k)) != 0) ? 1 : 0; + if (diag != nullptr && selected[k] == 0 && sol.bestWithCount[k] >= 0) { 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; + diag->lossMargin[keyOf(cand)] = + std::make_pair(sol.bestCount - sol.bestWithCount[k], sol.bestWithScore[k] - sol.bestScore); } } - - // ---- store the winners, book the blockers of the losers ---------- + } else { + // Too large to enumerate: greedy by score + ++nGreedyFallback; + std::vector usedTracks; + usedTracks.reserve(2 * nMembers); 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]++; + if (std::find(usedTracks.begin(), usedTracks.end(), cand.posID) != usedTracks.end() || + std::find(usedTracks.begin(), usedTracks.end(), cand.eleID) != usedTracks.end()) { continue; } - if (diag == nullptr) { + 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.push_back(keyOf(cand)); + nv0_map[cand.colID]++; + continue; + } + if (diag == nullptr) { + continue; + } + // all winners that took one leg + const CandKey victimKey = keyOf(cand); + for (int j = 0; j < nMembers; ++j) { + if (selected[j] == 0) { continue; } - // all winners that took one leg - 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) { + const auto& winner = vecV0Dedup[members[j]]; + if (winner.posID == cand.posID || winner.eleID == cand.eleID) { + diag->blockersByKey[victimKey].push_back(keyOf(winner)); + } + } + } + } + + 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 a.collisionId < b.collisionId; + }); + } + + // mode 1: candidates sorted by score (lower = better), a candidate is kept if both legs are still free + void selectGreedy(const size_t nTracks, DedupDiag* diag) + { + // 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(nTracks, false); + std::unordered_map ownerOfLeg; // diagnostics only + stored_v0Ids.clear(); + stored_fullv0Ids.clear(); + nv0_map.clear(); + for (const auto& v0Cand : vecV0Dedup) { + const CandKey candKey = keyOf(v0Cand); + 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; } - 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 (std::find(blockers.begin(), blockers.end(), itOwner->second) == blockers.end()) { + blockers.push_back(itOwner->second); // the same winner may hold both legs } } } + continue; } - - if (nGreedyFallback > 0) { - LOGF(warning, "group matching: %d conflict group(s) larger than dedupMaxGroupSize = %d were solved greedily", nGreedyFallback, maxEnum); + 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.push_back(keyOf(v0Cand)); + nv0_map[v0Cand.colID]++; + } - // 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); - }); + // sort accepted candidates by collision ID to be compatible with table + std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), + [](const auto& a, const auto& b) { + return a.collisionId < b.collisionId; + }); + } - } 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 - stored_v0Ids.clear(); - stored_fullv0Ids.clear(); - nv0_map.clear(); - 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}) { - 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 - } - } - } + template + void build(TCollisions const& collisions, TV0s const& v0s, TTracks const& tracks, TBCs const&, DedupDiag* diag = nullptr) + { + mRelaxedKeys.clear(); + // step 1: SVertexer V0s -> dedup candidates + for (const auto& collision : collisions) { + if constexpr (isMC) { + if (!collision.has_mcCollision()) { continue; } - 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); - nv0_map[v0Cand.colID]++; } - // sort accepted candidates by collision ID to be compatible with table - std::sort(stored_fullv0Ids.begin(), stored_fullv0Ids.end(), - [](const auto& a, const auto& b) { - return std::get<1>(a) < std::get<1>(b); - }); + if (!collision.isSelected()) { + continue; + } - // End of deduplication + // if constexpr (isTriggerAnalysis) { + // if (collision.triggerMask_raw() == 0) { + // continue; + // } + // } + + nv0_map[collision.globalIndex()] = 0; + + const auto& bc = collision.template foundBC_as(); + initCCDB(bc); + registry.fill(HIST("hCollisionCounter"), 1); + + updateCCDB(bc); // delay update until is needed + + vecV0Dedup.reserve(30000); // rough estimate for number of V0s per DF + const auto& v0s_per_coll = v0s.sliceBy(perCollision, collision.globalIndex()); + // LOGF(info, "n v0 = %d", v0s_per_coll.size()); + for (const auto& v0 : v0s_per_coll) { + // LOGF(info, "collision.globalIndex() = %d, v0.globalIndex() = %d, v0.posTrackId() = %d, v0.negTrackId() = %d", collision.globalIndex(), v0.globalIndex(), v0.posTrackId() , v0.negTrackId()); + fillV0Table(v0, false); + } // end of v0 loop + } // end of collision loop + + // photon finder + if (recovery.runFinder) { + runPhotonFinder(collisions, tracks); + } + + stored_v0Ids.reserve(pca_map.size()); // number of photon candidates per DF + stored_fullv0Ids.reserve(pca_map.size()); // number of photon candidates per DF + + // step 3: deduplication + switch (dedup.deduplicationMode.value) { + case V0DeduplicationMode::Pairwise: + selectPairwise(diag); + break; + case V0DeduplicationMode::KeepAll: + keepAll(); + break; + case V0DeduplicationMode::GroupMatching: + selectByConflictGroup(diag); + break; + default: + selectGreedy(tracks.size(), diag); + break; + } + + // split legs, ambiguity resolver (on the stored photons only) + if (dedup.resolveSplitLegs) { + resolveSplitLegConflicts(tracks, diag); + } + if (recovery.runSwapResolver) { + runSwapResolver(collisions, tracks, diag); } if (diag != nullptr) { @@ -1429,8 +1995,7 @@ struct PhotonConversionBuilder { const std::set storedSetSkim(stored_fullv0Ids.begin(), stored_fullv0Ids.end()); std::unordered_map> candsPerColl; for (const auto& c : vecV0Dedup) { - if (storedSetSkim.contains(std::make_tuple(static_cast(c.v0ID), static_cast(c.colID), - static_cast(c.posID), static_cast(c.eleID)))) { + if (storedSetSkim.contains(keyOf(c))) { candsPerColl[c.colID].push_back(&c); } } @@ -1458,26 +2023,38 @@ struct PhotonConversionBuilder { } } + // refit the stored photons and write the tables (SVertexer V0s through fillV0Table, the others directly) for (const auto& fullv0Id : stored_fullv0Ids) { - auto v0Id = std::get<0>(fullv0Id); - const auto collisionId = std::get<1>(fullv0Id); + auto v0Id = fullv0Id.v0Id; + const auto collisionId = fullv0Id.collisionId; if (skimRejected.contains(collisionId)) { continue; // collision failed the skim - do not store its photons/legs } - // auto posId = std::get<2>(fullv0Id); - // auto eleId = std::get<3>(fullv0Id); + // auto posId = fullv0Id.posId; + // auto eleId = fullv0Id.eleId; // LOGF(info, "!accept! | collision id = %d | v0id = %d , posid = %d , eleid = %d", collisionId, v0Id, posId, eleId); - auto v0 = v0s.rawIteratorAt(v0Id); + if (v0Id >= 0) { + 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)) { // 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()]); + } + + fillV0Table(v0, true); + continue; + } + const auto collision = collisions.rawIteratorAt(collisionId); 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)) { // o2-linter: disable=magic-number (nonsensical default value) + if (!(collision.neeuls() >= 1 || collision.neeuls() + nv0_map[collisionId] >= 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()]); } - - fillV0Table(v0, true); + fitV0Candidate(collision, tracks.rawIteratorAt(fullv0Id.posId), tracks.rawIteratorAt(fullv0Id.eleId), v0Id, true, + FitOptions{.nonSVertexer = true, .relaxAP = mRelaxedKeys.contains(fullv0Id)}); } // end of fullv0Id loop // for (const auto& collision : collisions) { @@ -1576,8 +2153,71 @@ struct PhotonConversionBuilder { } } - template - void fillDedupTruthDiagnostics(TTracks const& tracks, TMCParticles const& mcparticles, DedupDiag const& diag) + 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 + }; + + // which collision the copies of a true photon end up in + template + void fillCopyDiagnostics(TCollisions const& collisions, TMCParticles const& mcparticles, DedupDiag const& diag, std::map const& fate) + { + auto isRightColl = [&](int64_t colId, int64_t mid) { + const auto col = collisions.rawIteratorAt(colId); + return col.has_mcCollision() && col.mcCollisionId() == mcparticles.iteratorAt(mid).mcCollisionId(); + }; + auto dzPointing = [&](V0CandidateHelper const& c) { + const auto col = collisions.rawIteratorAt(c.colID); + const float pt = std::hypot(c.pxAtPV, c.pyAtPV); + return pt > 0.f ? std::fabs(c.svZ - std::hypot(c.svX - col.posX(), c.svY - col.posY()) * c.pzAtPV / pt - col.posZ()) : 999.f; // o2-linter: disable=magic-number (no direction) + }; + for (const auto& [mid, f] : fate) { + std::set copyColls; + bool rightAvailable = false; + for (const int& ci : f.cands) { + copyColls.insert(diag.snapshot[ci].colID); + rightAvailable = rightAvailable || isRightColl(diag.snapshot[ci].colID, mid); + } + int choice = 3; // o2-linter: disable=magic-number (bin 3 = lost) + if (f.alive) { + const bool storedRight = std::any_of(f.aliveColls.begin(), f.aliveColls.end(), [&](int64_t c) { return isRightColl(c, mid); }); + choice = storedRight ? 0 : (rightAvailable ? 1 : 2); + } + const int nCopies = static_cast(copyColls.size()); + registry.fill(HIST("MCDedup/Copies/hChoice"), static_cast(std::min(nCopies, 5)), static_cast(choice)); // o2-linter: disable=magic-number (last bin) + constexpr int MinCopies = 2; + if (nCopies < MinCopies || !rightAvailable) { + continue; + } + int bestDz = -1, bestPA = -1, bestPA3D = -1; + for (const int& ci : f.cands) { + const auto& c = diag.snapshot[ci]; + const float right = isRightColl(c.colID, mid) ? 1.f : 0.f; + registry.fill(HIST("MCDedup/Copies/hPAxy"), right, std::acos(std::clamp(c.cosPAxy, -1.f, 1.f))); + registry.fill(HIST("MCDedup/Copies/hPA3D"), right, std::acos(std::clamp(c.cosPA3D, -1.f, 1.f))); + if (bestDz < 0 || dzPointing(c) < dzPointing(diag.snapshot[bestDz])) { + bestDz = ci; + } + if (bestPA < 0 || c.cosPAxy > diag.snapshot[bestPA].cosPAxy) { + bestPA = ci; + } + if (bestPA3D < 0 || c.cosPA3D > diag.snapshot[bestPA3D].cosPA3D) { + bestPA3D = ci; + } + } + if (f.alive) { // same photons for all rules + registry.fill(HIST("MCDedup/Copies/hRule"), 0.f, choice == 0 ? 1.f : 0.f); + registry.fill(HIST("MCDedup/Copies/hRule"), 1.f, isRightColl(diag.snapshot[bestDz].colID, mid) ? 1.f : 0.f); + registry.fill(HIST("MCDedup/Copies/hRule"), 2.f, isRightColl(diag.snapshot[bestPA].colID, mid) ? 1.f : 0.f); + registry.fill(HIST("MCDedup/Copies/hRule"), 3.f, isRightColl(diag.snapshot[bestPA3D].colID, mid) ? 1.f : 0.f); + } + } + } + + template + void fillDedupTruthDiagnostics(TCollisions const& collisions, TTracks const& tracks, TMCParticles const& mcparticles, DedupDiag const& diag) { const int nCand = static_cast(diag.snapshot.size()); if (nCand == 0) { @@ -1591,8 +2231,7 @@ struct PhotonConversionBuilder { 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)); + key[i] = keyOf(c); indexOfKey[key[i]] = i; cls[i] = static_cast(classifyV0Truth(tracks.rawIteratorAt(c.posID), tracks.rawIteratorAt(c.eleID), mcparticles, motherPos[i], motherEle[i])); @@ -1602,12 +2241,6 @@ struct PhotonConversionBuilder { } } - 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) { @@ -1617,7 +2250,7 @@ struct PhotonConversionBuilder { f.cands.push_back(i); if (storedSet.contains(key[i]) > 0) { f.alive = true; - f.aliveColls.insert(std::get<1>(key[i])); + f.aliveColls.insert(key[i].collisionId); } if (f.best < 0 || diag.snapshot[i].score < diag.snapshot[f.best].score) { f.best = i; @@ -1683,7 +2316,19 @@ struct PhotonConversionBuilder { } } - // ---- truth pairs, grouped by MC collision ---------------------------- + if (dedup.cfgExtendedDiag) { + fillCopyDiagnostics(collisions, mcparticles, diag, fate); + } + + // stored cross-leg fakes by collisions + std::map, std::set> storedCrossColl; + for (int i = 0; i < nCand; ++i) { + if (cls[i] == kV0CrossLegFake && storedSet.contains(key[i])) { + storedCrossColl[{motherPos[i], motherEle[i]}].insert(key[i].collisionId); + } + } + + // truth pairs, grouped by MC collision std::map> mothersByMcColl; for (const auto& [mid, f] : fate) { mothersByMcColl[mcparticles.iteratorAt(mid).mcCollisionId()].push_back(mid); @@ -1720,6 +2365,21 @@ struct PhotonConversionBuilder { } continue; } + { // the pair stored as the crossed pairing (A+B-, B+A-) in one collision or one crossed pair only + const auto itAB = storedCrossColl.find({mids[a], mids[b]}); + const auto itBA = storedCrossColl.find({mids[b], mids[a]}); + bool bothCrossed = false; + if (itAB != storedCrossColl.end() && itBA != storedCrossColl.end()) { + for (const auto& c : itAB->second) { + bothCrossed = bothCrossed || itBA->second.contains(c); + } + } + if (bothCrossed) { + registry.fill(HIST("MCDedup/Pairs/surviveAsCrossedPairing/hQ"), qTrue); + } else if (itAB != storedCrossColl.end() || itBA != storedCrossColl.end()) { + registry.fill(HIST("MCDedup/Pairs/halfCrossedStored/hQ"), qTrue); + } + } if (nLost == 1) { fillTruthPair(kOneLost, qTrue, dEta); } else { @@ -1771,6 +2431,64 @@ struct PhotonConversionBuilder { } } + // MC, cfgExtendedDiag: class of the stored photons by source, stored photon pairs at low q, outcome of the resolver tests + template + void fillExtendedDiagnostics(TTracks const& tracks, TMCParticles const& mcparticles, DedupDiag const& diag) + { + struct Stored { + CandKey key; + int cls = 0; + int64_t mPos = -1; + int64_t mEle = -1; + std::array p{}; + }; + std::map candOfKey; + for (const auto& c : diag.snapshot) { + candOfKey[keyOf(c)] = &c; + } + std::map> storedOfColl; + for (const auto& key : diag.stored) { + const auto it = candOfKey.find(key); + if (it == candOfKey.end()) { + continue; + } + Stored s{key, 0, -1, -1, {it->second->pxAtPV, it->second->pyAtPV, it->second->pzAtPV}}; + s.cls = classifyV0Truth(tracks.rawIteratorAt(key.posId), tracks.rawIteratorAt(key.eleId), mcparticles, s.mPos, s.mEle); + const float pt = std::hypot(s.p[0], s.p[1]); + registry.fill(HIST("MCDedup/Stored/hClassSource"), static_cast(s.cls), key.v0Id < 0 ? 1.f : 0.f, pt); + storedOfColl[key.collisionId].push_back(s); + } + constexpr float MaxQ = 0.1f; + for (const auto& [col, photons] : storedOfColl) { + for (size_t i = 0; i < photons.size(); ++i) { + for (size_t j = i + 1; j < photons.size(); ++j) { + const auto& a = photons[i]; + const auto& b = photons[j]; + const float q = qInv(a.p[0], a.p[1], a.p[2], b.p[0], b.p[1], b.p[2]); + if (q > MaxQ) { + continue; + } + int cat = 2; + if (a.cls == kV0True && b.cls == kV0True && a.mPos != b.mPos) { + cat = 0; + } else if (a.cls == kV0CrossLegFake && b.cls == kV0CrossLegFake && a.mPos == b.mEle && a.mEle == b.mPos) { + cat = 1; + } + registry.fill(HIST("MCDedup/Stored/hPairClass"), static_cast(cat), q); + } + } + } + for (const auto& s : diag.swaps) { + const std::array, 4> legs = {{{s.keyI.posId, s.keyI.eleId}, {s.keyJ.posId, s.keyJ.eleId}, {s.keyI.posId, s.keyJ.eleId}, {s.keyJ.posId, s.keyI.eleId}}}; + std::array isTrue{}; + for (int k = 0; k < 4; ++k) { // o2-linter: disable=magic-number (four leg pairs) + int64_t mPos = -1, mEle = -1; + isTrue[k] = classifyV0Truth(tracks.rawIteratorAt(legs[k].first), tracks.rawIteratorAt(legs[k].second), mcparticles, mPos, mEle) == kV0True ? 1 : 0; + } + registry.fill(HIST("MCDedup/Swap/hCategory"), static_cast(3 * (isTrue[0] + isTrue[1]) + isTrue[2] + isTrue[3]), s.swapped ? 1.f : 0.f, s.q); + } + } + void processRec(MyCollisions const& collisions, FilteredV0s const& v0s, MyTracksIU const& tracks, MyBCs const& bcs) { build(collisions, v0s, tracks, bcs); @@ -1788,7 +2506,10 @@ struct PhotonConversionBuilder { { DedupDiag diag; build(collisions, v0s, tracks, bcs, &diag); - fillDedupTruthDiagnostics(tracks, mcparticles, diag); + fillDedupTruthDiagnostics(collisions, tracks, mcparticles, diag); + if (dedup.cfgExtendedDiag) { + fillExtendedDiagnostics(tracks, mcparticles, diag); + } } PROCESS_SWITCH(PhotonConversionBuilder, processMC, "process reconstructed info for MC", false); diff --git a/PWGEM/PhotonMeson/Utils/ConversionPointProfile.h b/PWGEM/PhotonMeson/Utils/ConversionPointProfile.h new file mode 100644 index 00000000000..2749319def2 --- /dev/null +++ b/PWGEM/PhotonMeson/Utils/ConversionPointProfile.h @@ -0,0 +1,177 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file ConversionPointProfile.h +/// \brief radius-profiled conversion point of a photon conversion (alternative to the helix touching point) +/// \author Stefanie Mrozinski, stefanie.mrozinski@cern.ch + +#ifndef PWGEM_PHOTONMESON_UTILS_CONVERSIONPOINTPROFILE_H_ +#define PWGEM_PHOTONMESON_UTILS_CONVERSIONPOINTPROFILE_H_ + +#include "Common/Core/RecoDecay.h" + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace o2::pwgem::photonmeson +{ + +constexpr int NProfPar = 4; // (rphi, z, phi, tgl) +constexpr int NProfR = 75; +constexpr float ProfRMin = 10.f; // trial radii 10 .. 84 cm + +/// \brief chi2 of a leg pair vs the trial radius, its minimum and class +struct RadiusProfile { + std::array chi2{}; + std::array ok{}; + int status = 4; + int iMin = -1; + float rMin = -1.f; + float chi2Min = -1.f; + float halfWidth = -1.f; + std::array xyz{-999.f, -999.f, -999.f}; +}; + +/// \brief parameters and covariance of one leg at every trial radius +struct LegScan { + std::array ok{}; + std::array, NProfR> m{}; // phi of the point, z, phi of the direction, tgl + std::array, NProfR> c{}; // lower triangle of the (rphi, z, phi, tgl) covariance + std::array, NProfR> pos{}; // global point +}; + +//_______________________________________________________________________ +/// \brief propagate one leg inward from its IU to a radius and keep its parameters there +/// \param trackIU leg at its innermost update +/// \param bz magnetic field +/// \param matCorr material correction of the propagation +/// \param out output: parameters and covariance per radius, ok = false outside the reach of the leg +inline void cpScanLeg(o2::track::TrackParCov const& trackIU, const float bz, const o2::base::Propagator::MatCorrType matCorr, LegScan& out) +{ + auto trk = trackIU; + const auto iu = trk.getXYZGlo(); + const float rIU = std::hypot(iu.X(), iu.Y()); + for (int k = NProfR - 1; k >= 0; --k) { + const float r = ProfRMin + k; + out.ok[k] = false; + if (r >= rIU - 0.1f) { // o2-linter: disable=magic-number (margin) + continue; + } + float xl = 0.f; + if (!trk.getXatLabR(r, xl, bz, o2::track::DirInward)) { + break; + } + if (!o2::base::Propagator::Instance()->propagateToX(trk, xl, bz, o2::base::PropagatorImpl::MAX_SIN_PHI, o2::base::PropagatorImpl::MAX_STEP, matCorr)) { + break; + } + const auto g = trk.getXYZGlo(); + const float snp = trk.getSnp(); + const float f = 1.f / std::sqrt(std::max(1.f - snp * snp, 1e-4f)); // o2-linter: disable=magic-number (floor) + out.m[k] = {std::atan2(g.Y(), g.X()), trk.getZ(), trk.getAlpha() + std::asin(std::clamp(snp, -1.f, 1.f)), trk.getTgl()}; + out.c[k] = {trk.getSigmaY2(), trk.getSigmaZY(), trk.getSigmaZ2(), trk.getSigmaSnpY() * f, trk.getSigmaSnpZ() * f, trk.getSigmaSnp2() * f * f, trk.getSigmaTglY(), trk.getSigmaTglZ(), trk.getSigmaTglSnp() * f, trk.getSigmaTgl2()}; + out.pos[k] = {g.X(), g.Y(), g.Z()}; + out.ok[k] = true; + } +} + +inline void cpProfileFinish(LegScan const& p, LegScan const& n, RadiusProfile& out); + +//_______________________________________________________________________ +/// \brief chi2 of a common point and a common direction of two legs at every trial radius, then minimum and class +/// \param p scan of the positive leg +/// \param n scan of the negative leg +/// \param zSig2Pos extra z variance of the positive leg (cm^2, legs without ITS) +/// \param zSig2Ele extra z variance of the negative leg (cm^2, legs without ITS) +/// \param out output: profile, minimum, half width and class +inline void cpRadiusProfile(LegScan const& p, LegScan const& n, const float zSig2Pos, const float zSig2Ele, RadiusProfile& out) +{ + out = RadiusProfile{}; + for (int k = 0; k < NProfR; ++k) { + out.ok[k] = p.ok[k] && n.ok[k]; + if (!out.ok[k]) { + continue; + } + const float r = ProfRMin + k; + ROOT::Math::SVector d; + d[0] = r * RecoDecay::constrainAngle(p.m[k][0] - n.m[k][0], -o2::constants::math::PI); + d[1] = p.m[k][1] - n.m[k][1]; + d[2] = RecoDecay::constrainAngle(p.m[k][2] - n.m[k][2], -o2::constants::math::PI); + d[3] = p.m[k][3] - n.m[k][3]; + ROOT::Math::SMatrix> cov; + int idx = 0; + for (int i = 0; i < NProfPar; ++i) { + for (int j = 0; j <= i; ++j) { + cov(i, j) = p.c[k][idx] + n.c[k][idx]; + ++idx; + } + } + cov(1, 1) += zSig2Pos + zSig2Ele; + if (!cov.Invert()) { + out.ok[k] = false; + continue; + } + out.chi2[k] = static_cast(ROOT::Math::Similarity(d, cov)); + if (out.iMin < 0 || out.chi2[k] < out.chi2Min) { + out.iMin = k; + out.chi2Min = out.chi2[k]; + } + } + cpProfileFinish(p, n, out); +} +//_______________________________________________________________________ +/// \brief minimum, half width of the dchi2 < 1 region, other minima and class of a filled profile +/// \param p scan of the positive leg (point at the minimum) +/// \param n scan of the negative leg (point at the minimum) +/// \param out profile, completed in place +inline void cpProfileFinish(LegScan const& p, LegScan const& n, RadiusProfile& out) +{ + if (out.iMin < 0) { + return; + } + out.rMin = ProfRMin + out.iMin; + out.xyz = {0.5f * (p.pos[out.iMin][0] + n.pos[out.iMin][0]), 0.5f * (p.pos[out.iMin][1] + n.pos[out.iMin][1]), 0.5f * (p.pos[out.iMin][2] + n.pos[out.iMin][2])}; + // valid range, dchi2 < 1 region around the minimum, other minima + int kLo = out.iMin; + int kHi = kLo; + while (kLo - 1 >= 0 && out.ok[kLo - 1] && out.chi2[kLo - 1] - out.chi2Min < 1.f) { + --kLo; + } + while (kHi + 1 < NProfR && out.ok[kHi + 1] && out.chi2[kHi + 1] - out.chi2Min < 1.f) { + ++kHi; + } + out.halfWidth = 0.5f * (kHi - kLo) + 0.5f; + const bool edge = (out.iMin == 0 || !out.ok[out.iMin - 1]) || (out.iMin == NProfR - 1 || !out.ok[out.iMin + 1]); + bool other = false; + for (int k = 1; k < NProfR - 1 && !other; ++k) { + if (k >= kLo && k <= kHi) { + continue; + } + if (out.ok[k] && out.ok[k - 1] && out.ok[k + 1] && out.chi2[k] <= out.chi2[k - 1] && out.chi2[k] <= out.chi2[k + 1] && out.chi2[k] - out.chi2Min < 1.f) { + other = true; + } + } + out.status = edge ? 3 : (other ? 2 : (out.halfWidth > 5.f ? 1 : 0)); // o2-linter: disable=magic-number (broad: half width above 5 cm) +} + +} // namespace o2::pwgem::photonmeson + +#endif // PWGEM_PHOTONMESON_UTILS_CONVERSIONPOINTPROFILE_H_ diff --git a/PWGEM/PhotonMeson/Utils/PCMUtilities.h b/PWGEM/PhotonMeson/Utils/PCMUtilities.h index e69703a3ec5..f8f6daed325 100644 --- a/PWGEM/PhotonMeson/Utils/PCMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/PCMUtilities.h @@ -21,12 +21,18 @@ #include #include +#include #include +#include #include +#include #include #include // IWYU pragma: keep (for rotate) -#include // IWYU pragma: keep (do not replace with Math/Vector2Dfwd.h) +#include +#include +#include +#include // IWYU pragma: keep (do not replace with Math/Vector2Dfwd.h) #include #include @@ -37,12 +43,22 @@ #include //_______________________________________________________________________ +/// \brief Armenteros-Podolanski cut of photon conversions: (alpha, qT) inside the ellipse with half axes alpha_max and qt_max +/// \param alpha Armenteros alpha +/// \param qt Armenteros qT (GeV/c) +/// \param alpha_max half axis in alpha +/// \param qt_max half axis in qT (GeV/c) +/// \return true if inside the ellipse inline bool checkAP(const float alpha, const float qt, const float alpha_max = 0.95, const float qt_max = 0.05) { float ellipse = std::pow(alpha / alpha_max, 2) + std::pow(qt / qt_max, 2); return (ellipse < 1.0); } //_______________________________________________________________________ +/// \brief Armenteros alpha of a V0: (pL+ - pL-) / (pL+ + pL-), pL = momentum of a leg along the V0 momentum +/// \param pxpos px of the positive leg (GeV/c), same for pypos, pzpos +/// \param pxneg px of the negative leg (GeV/c), same for pyneg, pzneg +/// \return alpha inline float v0Alpha(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { float momTot = RecoDecay::p(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); @@ -51,6 +67,10 @@ inline float v0Alpha(float pxpos, float pypos, float pzpos, float pxneg, float p return (lQlPos - lQlNeg) / (lQlPos + lQlNeg); // longitudinal momentum asymmetry of v0 } //_______________________________________________________________________ +/// \brief Armenteros qT of a V0: momentum of the negative leg transverse to the V0 momentum +/// \param pxpos px of the positive leg (GeV/c), same for pypos, pzpos +/// \param pxneg px of the negative leg (GeV/c), same for pyneg, pzneg +/// \return qT (GeV/c) inline float v0Qt(float pxpos, float pypos, float pzpos, float pxneg, float pyneg, float pzneg) { float momTot = RecoDecay::p2(pxpos + pxneg, pypos + pyneg, pzpos + pzneg); @@ -58,6 +78,13 @@ inline float v0Qt(float pxpos, float pypos, float pzpos, float pxneg, float pyne return std::sqrt(RecoDecay::p2(pxneg, pyneg, pzneg) - dp * dp / momTot); // qt of v0 } //_______________________________________________________________________ +/// \brief conversion point of two legs (Run 2 recipe): xy = touching point of the two helix circles (centres weighted with the +/// radius of the other circle), z = radius-weighted mean of the leg z after propagating each leg to its point closest to it +/// \param prop propagator (also gives the nominal field) +/// \param trackPosInformation positive leg +/// \param trackNegInformation negative leg +/// \param xyz output: conversion point (cm) +/// \param matCorr material correction of the propagation template inline void Vtx_recalculationParCov(o2::base::Propagator* prop, const o2::track::TrackParametrizationWithError& trackPosInformation, const o2::track::TrackParametrizationWithError& trackNegInformation, std::array& xyz, o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE) { @@ -187,6 +214,12 @@ inline std::array getPropMomentumFromTrackHelix(const float s, const T } //_______________________________________________________________________ +/// \brief conversion point of two legs from their track rows (Vtx_recalculationParCov) +/// \param prop propagator (also gives the nominal field) +/// \param lTrackPos positive leg +/// \param lTrackNeg negative leg +/// \param xyz output: conversion point (cm) +/// \param matCorr material correction of the propagation template inline void Vtx_recalculation(o2::base::Propagator* prop, T1 lTrackPos, T2 lTrackNeg, std::array& xyz, o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE) { @@ -262,6 +295,11 @@ inline V0TruthClass classifyV0Truth(TTrack const& pos, TTrack const& ele, TMCPar return kV0OtherFake; } +//_______________________________________________________________________ +/// \brief electron-PID compatibility of the two legs: mean |n sigma_e| + 2 * (n sigma_e+ + n sigma_e-) in the TPC +/// \param pos positive leg +/// \param ele negative leg +/// \return kappa template float getV0Kappa(TLeg const& pos, TLeg const& ele) { @@ -269,4 +307,139 @@ float getV0Kappa(TLeg const& pos, TLeg const& ele) return kappa; } +//_______________________________________________________________________ +/// \brief gap between the helix circles of two legs in xy +/// \param cPos circle of the positive leg +/// \param cEle circle of the negative leg +/// \return distance of the centres minus the sum of the radii (cm); < 0: the circles overlap +inline float circleGap(o2::track::TrackAuxPar const& cPos, o2::track::TrackAuxPar const& cEle) +{ + return std::hypot(cEle.xC - cPos.xC, cEle.yC - cPos.yC) - (cPos.rC + cEle.rC); +} + +/// \brief return value of tpcDriftSide if the side cannot be determined +constexpr int TPCSideUndetermined = 2; + +//_______________________________________________________________________ +/// \brief TPC side in which a track drifted, from the TPC side flags, else from z, tgl and the time asymmetry +/// \param track track row (flags, z, tgl, tpcDeltaTFwd, tpcDeltaTBwd) +/// \return 1 A side, -1 C side, 0 both sides (crosses the central electrode), TPCSideUndetermined +template +inline int tpcDriftSide(TTrack const& track) +{ + constexpr float MaxCECrossingDeltaTNS = 1000.f; + const auto sides = track.flags() & (o2::aod::track::TPCSideA | o2::aod::track::TPCSideC); + if (sides == o2::aod::track::TPCSideA) { + return 1; + } + if (sides == o2::aod::track::TPCSideC) { + return -1; + } + if (sides != 0) { + return 0; + } + const float dFwd = track.tpcDeltaTFwd(); + const float dBwd = track.tpcDeltaTBwd(); + if (dFwd == dBwd && dFwd < MaxCECrossingDeltaTNS) { + return 0; + } + const bool zPositive = track.z() > 0.f; + const bool tglPositive = track.tgl() > 0.f; + if (zPositive == tglPositive) { + return tglPositive ? 1 : -1; + } + if (dBwd == 0.f && dFwd > 0.f) { + return zPositive ? 1 : -1; + } + if (dBwd > 0.f) { + return tglPositive ? 1 : -1; + } + return TPCSideUndetermined; +} + +/// \brief z and momentum of a track helix at a point of its circle (helixPointAt) +struct HelixPoint { + bool ok = false; + bool inward = false; + float z = 0.f; + std::array p{}; +}; + +//_______________________________________________________________________ +/// \brief z and momentum of a track helix at an xy point of its circle +/// \param par track parameters +/// \param xP x of the point on the circle (cm) +/// \param yP y of the point on the circle (cm) +/// \param xB x of the beam position (cm) +/// \param yB y of the beam position (cm) +/// \param bz magnetic field (kG) +/// \return ok = false for a degenerate circle; inward = the track moves towards the beam at the point (more than 10 cm from +/// the beam), then the momentum is reversed; z and momentum at the point +inline HelixPoint helixPointAt(o2::track::TrackPar const& par, const float xP, const float yP, const float xB, const float yB, const float bz) +{ + HelixPoint h; + const o2::track::TrackAuxPar c(par, bz); + std::array r{}, p{}; + par.getXYZGlo(r); + par.getPxPyPzGlo(p); + const float pt = std::hypot(p[0], p[1]); + if (!(c.rC > 0.f) || !(pt > 0.f)) { + return h; + } + const float sense = (r[0] - c.xC) * p[1] - (r[1] - c.yC) * p[0] > 0.f ? 1.f : -1.f; + float back = RecoDecay::constrainAngle(sense * (std::atan2(r[1] - c.yC, r[0] - c.xC) - std::atan2(yP - c.yC, xP - c.xC)), 0.f); + const float tx = -sense * (yP - c.yC) / c.rC; + const float ty = sense * (xP - c.xC) / c.rC; + constexpr float MinRadiusOrientation = 10.f; + h.inward = std::hypot(xP - xB, yP - yB) > MinRadiusOrientation && tx * (xP - xB) + ty * (yP - yB) < 0.f; + if (!h.inward) { + if (back > o2::constants::math::PI && (r[0] - xB) * p[0] + (r[1] - yB) * p[1] > 0.f) { + back -= o2::constants::math::TwoPI; // o2-linter: disable=two-pi-add-subtract (signed arc length to a point behind the track, not an angle range) + } + h.z = r[2] - par.getTgl() * c.rC * back; + h.p = {pt * tx, pt * ty, pt * par.getTgl()}; + } else { + h.z = r[2] + par.getTgl() * c.rC * (back > 0.f ? o2::constants::math::TwoPI - back : 0.f); + h.p = {-pt * tx, -pt * ty, -pt * par.getTgl()}; + } + h.ok = true; + return h; +} + +//_______________________________________________________________________ +/// \brief compatibility of two track pieces with one track +/// \param trA first track +/// \param trB second track +/// \param bz magnetic field (kG) +/// \return chi2 of (y, snp, tgl, q/pt) at a common x with the summed covariances (4 dof); -1 if the propagation or the inversion fails +template +inline float splitPieceChi2(TTrack const& trA, TTrack const& trB, const float bz) +{ + auto a = getTrackParCov(trA); + auto b = getTrackParCov(trB); + if (!a.rotate(b.getAlpha()) || !a.propagateTo(b.getX(), bz)) { + a = getTrackParCov(trA); + b = getTrackParCov(trB); + if (!b.rotate(a.getAlpha()) || !b.propagateTo(a.getX(), bz)) { + return -1.f; + } + } + ROOT::Math::SMatrix> c; + c(0, 0) = a.getSigmaY2() + b.getSigmaY2(); + c(1, 0) = a.getSigmaSnpY() + b.getSigmaSnpY(); + c(1, 1) = a.getSigmaSnp2() + b.getSigmaSnp2(); + c(2, 0) = a.getSigmaTglY() + b.getSigmaTglY(); + c(2, 1) = a.getSigmaTglSnp() + b.getSigmaTglSnp(); + c(2, 2) = a.getSigmaTgl2() + b.getSigmaTgl2(); + c(3, 0) = a.getSigma1PtY() + b.getSigma1PtY(); + c(3, 1) = a.getSigma1PtSnp() + b.getSigma1PtSnp(); + c(3, 2) = a.getSigma1PtTgl() + b.getSigma1PtTgl(); + c(3, 3) = a.getSigma1Pt2() + b.getSigma1Pt2(); + if (!c.Invert()) { + return -1.f; + } + const ROOT::Math::SVector d(a.getY() - b.getY(), a.getSnp() - b.getSnp(), a.getTgl() - b.getTgl(), a.getQ2Pt() - b.getQ2Pt()); + return static_cast(ROOT::Math::Similarity(c, d)); +} + #endif // PWGEM_PHOTONMESON_UTILS_PCMUTILITIES_H_ diff --git a/PWGEM/PhotonMeson/Utils/V0DedupUtilities.h b/PWGEM/PhotonMeson/Utils/V0DedupUtilities.h new file mode 100644 index 00000000000..757d86b5a86 --- /dev/null +++ b/PWGEM/PhotonMeson/Utils/V0DedupUtilities.h @@ -0,0 +1,309 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file V0DedupUtilities.h +/// \brief helpers of the photon conversion builder: conflict-group solver of the deduplication and the geometry of the photon finder +/// \author Stefanie Mrozinski, stefanie.mrozinski@cern.ch + +#ifndef PWGEM_PHOTONMESON_UTILS_V0DEDUPUTILITIES_H_ +#define PWGEM_PHOTONMESON_UTILS_V0DEDUPUTILITIES_H_ + +#include "PWGEM/PhotonMeson/Utils/PCMUtilities.h" + +#include +#include + +#include +#include +#include +#include +#include +#include + +namespace o2::pwgem::photonmeson +{ + +//_______________________________________________________________________ +/// \brief q_inv of two massless particles (photons) +/// \param px1 px of the first photon (GeV/c), same for py1, pz1 +/// \param px2 px of the second photon (GeV/c), same for py2, pz2 +/// \return sqrt(2 (|p1| |p2| - p1 . p2)) (GeV/c) +inline float qInv(const float px1, const float py1, const float pz1, const float px2, const float py2, const float pz2) +{ + return std::sqrt(std::max(0.f, 2.f * (std::hypot(px1, py1, pz1) * std::hypot(px2, py2, pz2) - px1 * px2 - py1 * py2 - pz1 * pz2))); +} + +/// \brief result of solveConflictGroup +struct ConflictGroupSolution { + uint32_t bestMask = 0; + int bestCount = 0; + float bestScore = 0.f; + std::vector bestWithCount; // per candidate: size of the best subset that contains it + std::vector bestWithScore; // per candidate: score sum of that subset +}; + +//_______________________________________________________________________ +/// \brief best leg-disjoint subset of a conflict group of the deduplication (mode 2), exact enumeration of all subsets +/// \param legMask per candidate: bit mask of its legs (group-local bits); at most 16 candidates +/// \param score per candidate: dedup score, lower = better +/// \param s0 objective: <= 0 maximum number of candidates, then minimum score sum; > 0 maximum sum of max(s0 - score, 1e-3), +/// ties by the score sum +/// \param memberBest also find, per candidate, the best subset that contains it (MC diagnostics) +/// \return mask of the best subset, its size and score sum, and the per-candidate results +inline ConflictGroupSolution solveConflictGroup(std::vector const& legMask, std::vector const& score, const float s0, const bool memberBest) +{ + const int nMembers = static_cast(legMask.size()); + const bool weighted = s0 > 0.f; + constexpr float MinWeight = 1e-3f; + constexpr float WeightTolerance = 1e-6f; + ConflictGroupSolution sol; + sol.bestWithCount.assign(nMembers, -1); + sol.bestWithScore.assign(nMembers, 0.f); + float bestWeight = -1.f; + for (uint32_t mask = 1; mask < (1u << nMembers); ++mask) { + uint64_t legs = 0; + float scoreSum = 0.f; + float weightSum = 0.f; + int count = 0; + bool disjoint = true; + for (int k = 0; k < nMembers; ++k) { + if ((mask & (1u << k)) == 0) { + continue; + } + if ((legs & legMask[k]) != 0) { + disjoint = false; + break; + } + legs |= legMask[k]; + scoreSum += score[k]; + weightSum += std::max(s0 - score[k], MinWeight); + ++count; + } + if (!disjoint) { + continue; + } + const bool better = weighted ? (weightSum > bestWeight + WeightTolerance || (std::fabs(weightSum - bestWeight) <= WeightTolerance && scoreSum < sol.bestScore)) + : (count > sol.bestCount || (count == sol.bestCount && scoreSum < sol.bestScore)); + if (better) { + sol.bestCount = count; + sol.bestScore = scoreSum; + bestWeight = weightSum; + sol.bestMask = mask; + } + if (memberBest) { + for (int k = 0; k < nMembers; ++k) { + if ((mask & (1u << k)) == 0) { + continue; + } + if (count > sol.bestWithCount[k] || (count == sol.bestWithCount[k] && scoreSum < sol.bestWithScore[k])) { + sol.bestWithCount[k] = count; + sol.bestWithScore[k] = scoreSum; + } + } + } + } + return sol; +} + +/// \brief steps of the photon finder (pairing of free legs at the touching point of their helix circles); evalFinderPair returns the last step passed +enum FinderStep { + kFindTglWindow = 0, + kFindTime, + kFindGap, + kFindRTouch, + kFindPointing, + kFindAlongMotion, + kFindZMatch, + kFindCollision, + kNFindStep +}; + +/// \brief collision for the photon finder: time (ns, relative to the first BC of the DF), index, vertex (cm) +struct FinderCollision { + double t = 0.; + int64_t id = -1; + float x = 0.f, y = 0.f, z = 0.f; +}; + +/// \brief leg for the photon finder: track parameters, helix circle, TPC drift side and the time window (ns) of the track +struct FinderLeg { + int64_t id = -1; + float tgl = 0.f; + bool noITS = false; + bool movable = false; + int side = 0; + double t = 0., tMin = 0., tMax = 0.; + o2::track::TrackPar par; + o2::track::TrackAuxPar circle; +}; + +/// \brief cuts of the photon finder (cm) and the magnetic field (kG) +struct FinderCuts { + float maxDeltaTgl = 0.f; + float maxCircleGap = 0.f; + float minRTPConly = 0.f; + float maxRReqITS = 0.f; + float rMargin = 0.f; + float maxR = 0.f; + float minCosXY = 0.f; + float maxDeltaZ = 0.f; + float maxDz = 0.f; + float bz = 0.f; +}; + +//_______________________________________________________________________ +/// \brief geometric test of one leg pair of the photon finder: tgl, time windows, circle gap, touching radius, xy pointing to the +/// beam, direction of motion, z match at the touching point (incl. the drift of TPC-only legs), collision by z +/// \param p positive leg +/// \param e negative leg +/// \param colls collisions of the DF, sorted by time +/// \param cuts cuts of the finder +/// \param xB x of the beam position (cm) +/// \param yB y of the beam position (cm) +/// \param vDrift TPC drift velocity (cm/ns) +/// \param colId output: index of the matched collision, -1 if none +/// \return last FinderStep passed, -1 if the tgl window fails; kFindCollision if a collision matches +inline int evalFinderPair(FinderLeg const& p, FinderLeg const& e, std::vector const& colls, FinderCuts const& cuts, const float xB, const float yB, const double vDrift, int64_t& colId) +{ + colId = -1; + if (std::fabs(p.tgl - e.tgl) > cuts.maxDeltaTgl) { + return -1; + } + const double tLo = std::max(p.tMin, e.tMin); + const double tHi = std::min(p.tMax, e.tMax); + if (tLo > tHi) { + return kFindTglWindow; + } + const auto& cp = p.circle; + const auto& ce = e.circle; + const float dxC = ce.xC - cp.xC; + const float dyC = ce.yC - cp.yC; + const float dC = std::hypot(dxC, dyC); + if (std::fabs(circleGap(cp, ce)) > cuts.maxCircleGap) { + return kFindTime; + } + const float xT = (cp.xC * ce.rC + ce.xC * cp.rC) / (cp.rC + ce.rC); + const float yT = (cp.yC * ce.rC + ce.yC * cp.rC) / (cp.rC + ce.rC); + const float rTouch = std::hypot(xT, yT); + const float rMin = (p.noITS && e.noITS ? cuts.minRTPConly : cuts.maxRReqITS) - cuts.rMargin; + if (rTouch > cuts.maxR || rTouch < rMin) { + return kFindGap; + } + const float dB = std::hypot(xT - xB, yT - yB); + if (!(dC > 0.f) || !(dB > 0.f) || std::fabs((xT - xB) * (-dyC / dC) + (yT - yB) * (dxC / dC)) / dB < cuts.minCosXY) { + return kFindRTouch; + } + auto nearest = [&](o2::track::TrackAuxPar const& c, float& x, float& y) { + const float d = std::max(std::hypot(xT - c.xC, yT - c.yC), 1e-6f); + x = c.xC + c.rC * (xT - c.xC) / d; + y = c.yC + c.rC * (yT - c.yC) / d; + }; + float xP = 0.f, yP = 0.f, xE = 0.f, yE = 0.f; + nearest(cp, xP, yP); + nearest(ce, xE, yE); + const auto hp = helixPointAt(p.par, xP, yP, xB, yB, cuts.bz); + const auto he = helixPointAt(e.par, xE, yE, xB, yB, cuts.bz); + if (!hp.ok || !he.ok || hp.inward || he.inward) { + return kFindPointing; + } + const double b = vDrift * (p.side - e.side); + const double a = (hp.z - p.side * vDrift * p.t) - (he.z - e.side * vDrift * e.t); + const double tBest = b == 0. ? tLo : std::clamp(-a / b, tLo, tHi); + if (std::fabs(a + b * tBest) > cuts.maxDeltaZ) { + return kFindAlongMotion; + } + const float ptG = std::hypot(hp.p[0] + he.p[0], hp.p[1] + he.p[1]); + const float tglG = ptG > 0.f ? (hp.p[2] + he.p[2]) / ptG : 0.f; + float best = std::numeric_limits::max(); + for (auto it = std::lower_bound(colls.begin(), colls.end(), tLo, [](const FinderCollision& c, double t) { return c.t < t; }); it != colls.end() && it->t <= tHi; ++it) { + const auto zPos = static_cast(hp.z + p.side * vDrift * (it->t - p.t)); + const auto zEle = static_cast(he.z + e.side * vDrift * (it->t - e.t)); + const float deltaZ = zPos - zEle; + const float dz = (zPos * ce.rC + zEle * cp.rC) / (cp.rC + ce.rC) - std::hypot(xT - it->x, yT - it->y) * tglG - it->z; + const float dist = std::hypot(dz, deltaZ); + if (dist < best && std::fabs(dz) < cuts.maxDz && std::fabs(deltaZ) < cuts.maxDeltaZ) { + best = dist; + colId = it->id; + } + } + return colId >= 0 ? kFindCollision : kFindZMatch; +} + +/// \brief input features of the finder purity score (naive-Bayes log-likelihood ratio, finderScore_v1, PbPb MC LHC26a5a) +enum FinderScoreFeature { + kScoreCosPA = 0, // 3D cosine of the pointing angle + kScorePCA, // distance of the legs at the conversion point (cm) + kScoreMee, // e+e- mass at the conversion point (GeV/c^2) + kScoreGap, // helix-circle gap of the legs (cm) + kScoreDTgl, // |tgl+ - tgl-| + kScorePsiPair, + kScoreMaxNSigmaEl, // max(|n sigma_e|) of the legs in the TPC + kScoreMinNSigmaPi, // min(|n sigma_pi|) of the legs in the TPC + kNScoreFeatures +}; + +/// \brief log-likelihood-ratio table of one feature: nBins uniform bins in [lo, hi), llr[0] underflow, llr[nBins + 1] overflow +/// (same bin lookup as TH1::FindBin) +struct FinderScoreTable { + int nBins; + double lo; + double hi; + const double* llr; +}; + +/// \brief log(p_true / p_fake) per bin of each feature (finderScore_v1) +inline constexpr std::array FinderLlrCosPA = {-0.42326231102291223, -3.0746401250262112, -2.902679468526141, -2.7737498628462491, -2.8653323764415854, -2.9977154358984874, -2.879287982624958, -2.9449559801189884, -2.8305718792936898, -2.7610493575753212, -2.6743115383729115, -2.739759837578176, -2.7445126982363175, -2.6010956990444813, -2.5377329273302194, -2.4413963521466235, -2.5007549142569676, -2.5484647707016048, -2.3339023465737707, -2.3706555897531327, -2.3350661347670365, -2.3987685637953868, -2.2562477335729185, -2.2430310349096305, -2.1429211398936303, -2.1233771903398493, -2.0022839108013391, -1.9558000169723206, -1.8729479034124974, -1.9065302845458256, -1.8542708039089619, -1.7230031706190478, -1.6490950919627059, -1.633800343346663, -1.5180730364736743, -1.389137677593709, -1.3181114339137541, -1.2615432549457208, -1.1021895055813662, -0.97913284177589832, -0.97669632321874755, -0.81070582095170196, -0.70030094035730539, -0.55130110346066707, -0.34392181466201621, -0.16444806113591265, 0.025741391801742678, 0.22392951134446418, 0.48449501950084528, 0.74029958080748082, 1.3544181524605163, 2.1416870464386246}; +inline constexpr std::array FinderLlrPca = {-0.42328563731423596, 0.74173300442294698, 0.60064005460866021, 0.48774214670698213, 0.43172212684872524, 0.34949695425445604, 0.2876013953788461, 0.13114334603545452, -0.013681195276045386, -0.10436433242742214, -0.27099684361934573, -0.38200720615175005, -0.47403001970996045, -0.57728096929145511, -0.65985610396473016, -0.75403841222733481, -0.76619973024014831, -0.76182156236311549, -0.84274340404963699, -0.93501651159467125, -0.99602554213200045, -1.029448537116822, -1.1240438482723172, -1.124554873172243, -1.2439733445571921, -1.2077440689826793, -1.3081144381420802, -1.3626265850153199, -1.3892300573742717, -1.3978783977548588, -1.4463232066948644, -0.42328563731423596}; +inline constexpr std::array FinderLlrMee = {-0.42329146929478079, 0.82940375143002754, 0.6225570701387575, 0.30541176737301567, -0.042962035242895512, -0.43745174446394713, -0.78901712806679281, -1.1834130413160895, -1.5486037315484942, -1.8446313857300292, -2.2387380284795446, -2.1811493868471543, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079, -0.42329146929478079}; +inline constexpr std::array FinderLlrGap = {-0.42329730143842487, -3.0127420263468552, -2.6803080781356927, -2.2310834132674926, -1.852433867664681, -1.4426351719772665, -1.1194240909923954, -0.72769983474283484, -0.24592785900768224, 0.40475570320106047, 1.4082550958014817, 1.1535415604659081, 0.43945228763843225, 0.14326171246634631, -0.19523294768817565, -0.39258756076294099, -0.41539680606824775, -0.42329730143842487, -0.42329730143842487, -0.42329730143842487, -0.42329730143842487, -0.42329730143842487}; +inline constexpr std::array FinderLlrDTgl = {-0.42330313374517559, 0.84513809550928687, 0.46735826574627454, 0.13082343003051436, -0.12182979016333238, -0.33846868362980387, -0.52696466906551476, -0.72450796222957359, -0.89155673129371438, -1.0310133083500654, -1.2327015991375558, -1.3726345840013925, -1.4519781151564646, -1.7004785855213946, -1.7559986132128569, -1.818203125662764, -0.42330313374517559}; +inline constexpr std::array FinderLlrPsiPair = {-0.42330196727077668, 1.5226081817845367, 0.42399589311642699, 1.6969615689293145, 1.6135799599902634, 0.44573587975283285, 0.50095693425255516, 0.43535965176674207, 0.0028547913604746104, -0.0059523726814127804, 0.61779960880440121, 0.82263503615419142, 0.81841116503800626, 0.99808371366038406, 0.57522686284035063, 0.80047346435133904, 2.521137011895664, 2.2392858597546761}; +inline constexpr std::array FinderLlrMaxNSigmaEl = {-0.42331129924865812, 0.58307082727859216, 0.52708190318359982, 0.41704446150850366, 0.16536503067868297, -0.38368187287820599, -1.4189814933894278, -3.0185512358113016, -5.190101150118414, -0.42331129924865812}; +inline constexpr std::array FinderLlrMinNSigmaPi = {-0.42330896621503916, -0.9609371208084968, -0.46431089172823747, -0.06050992822674324, 0.16077238768804142, 0.3920347340880036, 0.64255032840474002, 0.82818496003043784, 0.96447986209716363, 1.0700938177484094, 1.1082284119787249, 1.017209914244056}; + +/// \brief binning and LLR table of each feature, in the order of FinderScoreFeature +inline constexpr std::array FinderScoreTables = {{ + {50, 0.99000000953674316, 1, FinderLlrCosPA.data()}, + {30, 0, 3, FinderLlrPca.data()}, + {25, 0, 0.05000000074505806, FinderLlrMee.data()}, + {20, -5, 5, FinderLlrGap.data()}, + {15, 0, 0.15000000596046448, FinderLlrDTgl.data()}, + {16, -1.6000000238418579, 1.6000000238418579, FinderLlrPsiPair.data()}, + {8, 0, 4, FinderLlrMaxNSigmaEl.data()}, + {10, 0, 10, FinderLlrMinNSigmaPi.data()}, +}}; + +//_______________________________________________________________________ +/// \brief purity score of a finder candidate: log(p_true / p_fake) summed over the 8 features (naive Bayes) +/// \param x feature values in the order of FinderScoreFeature +/// \return score; > 0: more likely a true photon +inline float finderPurityScore(std::array const& x) +{ + float score = 0.f; + for (int i = 0; i < kNScoreFeatures; ++i) { + const auto& t = FinderScoreTables[i]; + const double v = x[i]; + int bin = 0; + if (v < t.lo) { + bin = 0; + } else if (!(v < t.hi)) { + bin = t.nBins + 1; + } else { + bin = 1 + static_cast(t.nBins * (v - t.lo) / (t.hi - t.lo)); + } + score += t.llr[bin]; + } + return score; +} + +} // namespace o2::pwgem::photonmeson + +#endif // PWGEM_PHOTONMESON_UTILS_V0DEDUPUTILITIES_H_