From 97d2afed290c4b3205e2c8355281bbb4183c319c Mon Sep 17 00:00:00 2001 From: skundu692 Date: Fri, 30 Jan 2026 18:52:54 +0100 Subject: [PATCH 1/7] Add process function for systematic --- PWGLF/Tasks/Resonances/f1protoncorrelation.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx index 63568887fab..2627cb00e01 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -952,6 +952,7 @@ struct f1protoncorrelation { { const float maxMomPi = maxMomentumPion; const float maxMomK = maxMomentumKaon; + // const float pTofPiMin = momentumTOFPionMin; // const float pTofPiMax = momentumTOFPionMax; // const float pTofKMin = momentumTOFKaonMin; From b422afdf6aeaf9aea216d516fcca5ca4ea731688 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sun, 1 Feb 2026 11:04:28 +0100 Subject: [PATCH 2/7] Fix clang issue --- PWGLF/Tasks/Resonances/f1protoncorrelation.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx index 2627cb00e01..63568887fab 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -952,7 +952,6 @@ struct f1protoncorrelation { { const float maxMomPi = maxMomentumPion; const float maxMomK = maxMomentumKaon; - // const float pTofPiMin = momentumTOFPionMin; // const float pTofPiMax = momentumTOFPionMax; // const float pTofKMin = momentumTOFKaonMin; From 0b44af2ab0c0367212361720e2a6d09fe751d7db Mon Sep 17 00:00:00 2001 From: skundu692 Date: Wed, 25 Feb 2026 19:58:49 +0100 Subject: [PATCH 3/7] Old process function --- PWGLF/Tasks/Resonances/f1protoncorrelation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx index 63568887fab..e96cc342fdc 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -1056,7 +1056,7 @@ struct f1protoncorrelation { histos.fill(HIST("hPhaseSpaceProtonKaonSame"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, protontrack.protonCharge(), kaonCharge, bz, bz), relative_momentum); // Phase Space Proton kaon if (pionCharge == protontrack.protonCharge()) histos.fill(HIST("hPhaseSpaceProtonPionSame"), Proton.Eta() - Pion.Eta(), PhiAtSpecificRadiiTPC(Proton, Pion, protontrack.protonCharge(), pionCharge, bz, bz), relative_momentum); // Phase Space Proton Pion - histos.fill(HIST("h2SameEventf1pptCorrelation"), F1.M(), relative_momentum, Proton.Pt()); + histos.fill(HIST("h2SameEventf1pptCorrelation"), F1.M(), relative_momentum, Proton.Pt()); } activePair.push_back(sysId); } From c89738745961637926210dacfef5cfa95b417550 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 18:33:50 +0100 Subject: [PATCH 4/7] Add Event loss in phi meson task and add MC process function for spincorrelation task --- .../Strangeness/lambdaspincorrderived.cxx | 56 ++++++++++++++++--- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 51e13a779f1..7a6e94d7d9e 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1132,19 +1132,36 @@ struct lambdaspincorrderived { if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) continue; - // Collect partners from nominal key, plus wrapped neighbor only for φ-edge bins std::vector matches; - matches.reserve(128); // or keep binVec.size() if you prefer + const int maxKeep = maxMatchesPerPair.value; // default 25 + matches.reserve(std::max(64, maxKeep > 0 ? maxKeep : 64)); + const int64_t curColIdx = static_cast(collision1.index()); + std::unordered_set seenRow; + seenRow.reserve(static_cast(std::max(256, 4 * (maxKeep > 0 ? maxKeep : 64)))); - auto collectFrom = [&](int phiBinUse) { - const size_t keyUse = linearKey(colBin, status, ptB, etaB, phiBinUse, mB, + auto collectFrom = [&](int ptUse, int etaUse, int phiUse) { + if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { + return; // early stop + } + + const size_t keyUse = linearKey(colBin, status, ptUse, etaUse, phiUse, mB, nStat, nPt, nEta, nPhi, nM); auto const& vec = buffer[keyUse]; + for (const auto& bc : vec) { + if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { + break; + } if (bc.collisionIdx == curColIdx) { continue; // must be from different event } + + // dedupe first + if (!seenRow.insert(bc.rowIndex).second) { + continue; + } + auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); if (!selectionV0(tX)) { continue; @@ -1152,17 +1169,37 @@ struct lambdaspincorrderived { if (!checkKinematics(t1, tX)) { continue; } + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); } }; // 1) nominal φ-bin collectFrom(phiB); - // 2) wrap only at boundaries: 0 <-> nPhi-1 - if (phiB == 0) { - collectFrom(nPhi - 1); - } else if (phiB == nPhi - 1) { - collectFrom(0); + // scan pt±1, eta±1, phi±1 (wrapped) + for (int dpt = -1; dpt <= 1; ++dpt) { + const int ptUse = ptB + dpt; + if (ptUse < 0 || ptUse >= nPt) { + continue; + } + for (int deta = -1; deta <= 1; ++deta) { + const int etaUse = etaB + deta; + if (etaUse < 0 || etaUse >= nEta) { + continue; + } + for (int phiUse : phiBins) { + collectFrom(ptUse, etaUse, phiUse); + if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { + break; + } + } + if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { + break; + } + } + if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { + break; + } } if (matches.empty()) { @@ -1455,6 +1492,7 @@ struct lambdaspincorrderived { // enable it PROCESS_SWITCH(lambdaspincorrderived, processMCMEV3, "Process MC ME (MEV3)", false); + // ----------------------------------------------------- // 5) MC Event Mixing using your MEV4 6D-buffer approach // ----------------------------------------------------- From f92acf0e7948090ccb8e34a706e42fbba36224d5 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 20:29:09 +0100 Subject: [PATCH 5/7] Add new MC process function for mixing --- .../Strangeness/lambdaspincorrderived.cxx | 26 +++---------------- 1 file changed, 4 insertions(+), 22 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 7a6e94d7d9e..93191f8a934 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1132,36 +1132,19 @@ struct lambdaspincorrderived { if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) continue; + // Collect partners from nominal key, plus wrapped neighbor only for φ-edge bins std::vector matches; - const int maxKeep = maxMatchesPerPair.value; // default 25 - matches.reserve(std::max(64, maxKeep > 0 ? maxKeep : 64)); - + matches.reserve(128); // or keep binVec.size() if you prefer const int64_t curColIdx = static_cast(collision1.index()); - std::unordered_set seenRow; - seenRow.reserve(static_cast(std::max(256, 4 * (maxKeep > 0 ? maxKeep : 64)))); - - auto collectFrom = [&](int ptUse, int etaUse, int phiUse) { - if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { - return; // early stop - } - const size_t keyUse = linearKey(colBin, status, ptUse, etaUse, phiUse, mB, + auto collectFrom = [&](int phiBinUse) { + const size_t keyUse = linearKey(colBin, status, ptB, etaB, phiBinUse, mB, nStat, nPt, nEta, nPhi, nM); auto const& vec = buffer[keyUse]; - for (const auto& bc : vec) { - if (maxKeep > 0 && static_cast(matches.size()) >= maxKeep) { - break; - } if (bc.collisionIdx == curColIdx) { continue; // must be from different event } - - // dedupe first - if (!seenRow.insert(bc.rowIndex).second) { - continue; - } - auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); if (!selectionV0(tX)) { continue; @@ -1169,7 +1152,6 @@ struct lambdaspincorrderived { if (!checkKinematics(t1, tX)) { continue; } - matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); } }; From c04115bbf6737ec9fa291869470cef682ef62a99 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 20:55:25 +0100 Subject: [PATCH 6/7] Fix clang error --- PWGLF/Tasks/Resonances/f1protoncorrelation.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx index e96cc342fdc..63568887fab 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -1056,7 +1056,7 @@ struct f1protoncorrelation { histos.fill(HIST("hPhaseSpaceProtonKaonSame"), Proton.Eta() - Kaon.Eta(), PhiAtSpecificRadiiTPC(Proton, Kaon, protontrack.protonCharge(), kaonCharge, bz, bz), relative_momentum); // Phase Space Proton kaon if (pionCharge == protontrack.protonCharge()) histos.fill(HIST("hPhaseSpaceProtonPionSame"), Proton.Eta() - Pion.Eta(), PhiAtSpecificRadiiTPC(Proton, Pion, protontrack.protonCharge(), pionCharge, bz, bz), relative_momentum); // Phase Space Proton Pion - histos.fill(HIST("h2SameEventf1pptCorrelation"), F1.M(), relative_momentum, Proton.Pt()); + histos.fill(HIST("h2SameEventf1pptCorrelation"), F1.M(), relative_momentum, Proton.Pt()); } activePair.push_back(sysId); } From 4cfeb6aef858ec640fec40f4735740bc08579fbf Mon Sep 17 00:00:00 2001 From: skundu692 Date: Tue, 3 Mar 2026 21:19:20 +0100 Subject: [PATCH 7/7] Add new process function for mixing --- .../Strangeness/lambdaspincorrderived.cxx | 424 +++++++++++++++++- 1 file changed, 417 insertions(+), 7 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 93191f8a934..a4171e50263 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -203,6 +203,20 @@ struct lambdaspincorrderived { Configurable ConfWeightPathALL2{"ConfWeightPathALL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; // event sel///////// + + Configurable cfgV5NeighborPt{"cfgV5NeighborPt", 0, "v5: neighbor bins in pT (use symmetric ±N, edge-safe)"}; + Configurable cfgV5NeighborEta{"cfgV5NeighborEta", 0, "v5: neighbor bins in eta (use symmetric ±N, edge-safe)"}; + Configurable cfgV5NeighborPhi{"cfgV5NeighborPhi", 0, "v5: neighbor bins in phi (use symmetric ±N, periodic wrap)"}; + + Configurable cfgV5MaxAbsDeltaPairPt{"cfgV5MaxAbsDeltaPairPt", 0.1f, "v5: |pT(pair_ME)-pT(pair_SE)| max"}; + Configurable cfgV5MaxAbsDeltaPairEta{"cfgV5MaxAbsDeltaPairEta", 0.05f, "v5: |eta(pair_ME)-eta(pair_SE)| max"}; + Configurable cfgV5MaxAbsDeltaPairPhi{"cfgV5MaxAbsDeltaPairPhi", 0.03f, "v5: |dphi(pair_ME,pair_SE)| max (periodic)"}; + + // optional cap to keep it fast + Configurable cfgV5MaxMatches{"cfgV5MaxMatches", 50, "v5: max ME replacements per SE pair (after all cuts)"}; + + Configurable cfgMaxMatchesPerSE{"cfgMaxMatchesPerSE", 50, "Cap number of ME matches per SE pair (0=all)"}; + Configurable cfgMixSeed{"cfgMixSeed", 0xdecafbadULL, "RNG seed for downsampling matches (deterministic)"}; Configurable maxMatchesPerPair{"maxMatchesPerPair", 25, "Max mixed candidates per (t1,t2)"}; Configurable centMin{"centMin", 0, "Minimum Centrality"}; Configurable centMax{"centMax", 80, "Maximum Centrality"}; @@ -416,9 +430,6 @@ struct lambdaspincorrderived { if (std::abs(RecoDecay::constrainAngle(RecoDecay::constrainAngle(candidate1.lambdaPhi(), 0.f, harmonic) - RecoDecay::constrainAngle(candidate2.lambdaPhi(), 0.f, harmonic), -TMath::Pi(), 1)) > phiMix) { return false; } - /*if (std::abs(RecoDecay::constrainAngle(candidate1.lambdaPhi(), 0.0F, harmonic) - RecoDecay::constrainAngle(candidate2.lambdaPhi(), 0.0F, harmonic)) > phiMix) { - return false; - }*/ if (std::abs(candidate1.lambdaMass() - candidate2.lambdaMass()) > massMix) { return false; } @@ -1276,10 +1287,7 @@ struct lambdaspincorrderived { if (std::abs(mcacc::lamEta(c1) - mcacc::lamEta(c2)) > etaMix) { return false; } - if (std::abs(RecoDecay::constrainAngle( - RecoDecay::constrainAngle(mcacc::lamPhi(c1), 0.f, harmonic) - - RecoDecay::constrainAngle(mcacc::lamPhi(c2), 0.f, harmonic), - -TMath::Pi(), 1)) > phiMix) { + if (std::abs(RecoDecay::constrainAngle(RecoDecay::constrainAngle(mcacc::lamPhi(c1), 0.f, harmonic) - RecoDecay::constrainAngle(mcacc::lamPhi(c2), 0.f, harmonic), -TMath::Pi(), 1)) > phiMix) { return false; } if (std::abs(mcacc::lamMass(c1) - mcacc::lamMass(c2)) > massMix) { @@ -1646,6 +1654,408 @@ struct lambdaspincorrderived { } } PROCESS_SWITCH(lambdaspincorrderived, processMCMEV4, "Process MC ME (5d buffer)", false); + + static inline float phi0To2Pi(float phi) + { + // harmonic=1, min=0 => [0, 2pi) + return RecoDecay::constrainAngle(phi, 0.0f, 1); + } + + static inline float deltaPhiMinusPiToPi(float phiA, float phiB) + { + // returns in [-pi, pi) + const float d = phi0To2Pi(phiA) - phi0To2Pi(phiB); + return RecoDecay::constrainAngle(d, -TMath::Pi(), 1); + } + + static inline float absDeltaPhi(float phiA, float phiB) + { + return std::abs(deltaPhiMinusPiToPi(phiA, phiB)); + } + + // symmetric neighbors for continuous mixing (pt/eta): include bin, ±1, ±2..., edge-safe + static inline void collectNeighborBins1D(int b, int nBins, int nNeighbor, std::vector& out) + { + out.clear(); + out.reserve(2 * nNeighbor + 1); + for (int d = -nNeighbor; d <= nNeighbor; ++d) { + const int bb = b + d; + if (bb < 0 || bb >= nBins) { + continue; + } + out.push_back(bb); + } + std::sort(out.begin(), out.end()); + out.erase(std::unique(out.begin(), out.end()), out.end()); + } + + // symmetric neighbors for phi: periodic wrap + static inline void collectNeighborBinsPhi(int b, int nPhi, int nNeighbor, std::vector& out) + { + out.clear(); + out.reserve(2 * nNeighbor + 1); + for (int d = -nNeighbor; d <= nNeighbor; ++d) { + int bb = b + d; + bb %= nPhi; + if (bb < 0) { + bb += nPhi; + } + out.push_back(bb); + } + std::sort(out.begin(), out.end()); + out.erase(std::unique(out.begin(), out.end()), out.end()); + } + + template + static inline bool passPairClosenessV5(LVec const& l1SE, LVec const& l2SE, + LVec const& l1ME, LVec const& l2ME, + float maxAbsDeltaAbsDPt, + float maxAbsDeltaAbsDEta, + float maxAbsDeltaAbsDPhi) + { + const float dPtSE = std::abs(static_cast(l1SE.Pt() - l2SE.Pt())); + const float dEtaSE = std::abs(static_cast(l1SE.Eta() - l2SE.Eta())); + const float dPhiSE = std::abs(RecoDecay::constrainAngle( + RecoDecay::constrainAngle(static_cast(l1SE.Phi()), 0.f, 1) - + RecoDecay::constrainAngle(static_cast(l2SE.Phi()), 0.f, 1), + -TMath::Pi(), 1)); + + const float dPtME = std::abs(static_cast(l1ME.Pt() - l2ME.Pt())); + const float dEtaME = std::abs(static_cast(l1ME.Eta() - l2ME.Eta())); + const float dPhiME = std::abs(RecoDecay::constrainAngle( + RecoDecay::constrainAngle(static_cast(l1ME.Phi()), 0.f, 1) - + RecoDecay::constrainAngle(static_cast(l2ME.Phi()), 0.f, 1), + -TMath::Pi(), 1)); + + if (std::abs(dPtME - dPtSE) > maxAbsDeltaAbsDPt) + return false; + if (std::abs(dEtaME - dEtaSE) > maxAbsDeltaAbsDEta) + return false; + if (std::abs(dPhiME - dPhiSE) > maxAbsDeltaAbsDPhi) + return false; + return true; + } + + void processMEV5(EventCandidates const& collisions, AllTrackCandidates const& V0s) + { + MixBinner mb{ptMin.value, ptMax.value, ptMix.value, + v0eta.value, etaMix.value, + phiMix.value}; + + const int nCol = colBinning.getAllBinsCount(); + const int nStat = N_STATUS; + const int nPt = mb.nPt(); + const int nEta = mb.nEta(); + const int nPhi = mb.nPhi(); + + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi; + std::vector> buffer(nKeys); + + auto key5 = [&](int colBin, int status, int ptB, int etaB, int phiB) -> size_t { + // linearKey without mass + // ((((col* nStat + status)*nPt + ptB)*nEta + etaB)*nPhi + phiB) + size_t k = static_cast(colBin); + k = k * static_cast(nStat) + static_cast(status); + k = k * static_cast(nPt) + static_cast(ptB); + k = k * static_cast(nEta) + static_cast(etaB); + k = k * static_cast(nPhi) + static_cast(phiB); + return k; + }; + + // ---- PASS 1: fill buffer ---- + for (auto const& col : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(col.posz(), col.cent())); + auto slice = V0s.sliceBy(tracksPerCollisionV0, col.index()); + + for (auto const& t : slice) { + if (!selectionV0(t)) + continue; + + const int status = static_cast(t.v0Status()); + if (status < 0 || status >= nStat) + continue; + + const int ptB = mb.ptBin(t.lambdaPt()); + const int etaB = mb.etaBin(t.lambdaEta()); + const int phiB = mb.phiBin(phi0To2Pi(t.lambdaPhi())); + if (ptB < 0 || etaB < 0 || phiB < 0) + continue; + + buffer[key5(colBin, status, ptB, etaB, phiB)].push_back(BufferCand{ + .collisionIdx = static_cast(col.index()), + .rowIndex = static_cast(t.globalIndex()), // keep your choice + .v0Status = static_cast(status), + .ptBin = static_cast(ptB), + .etaBin = static_cast(etaB), + .phiBin = static_cast(phiB), + .mBin = 0}); + } + } + + std::vector ptBins, etaBins, phiBins; + std::vector matches; + matches.reserve(256); + + // ---- PASS 2: for each SE pair, find ME leg-1 replacements ---- + for (auto const& col1 : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); + auto poolA = V0s.sliceBy(tracksPerCollisionV0, col1.index()); + + for (auto const& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { + if (!selectionV0(t1) || !selectionV0(t2)) + continue; + if (t2.index() <= t1.index()) + continue; + + // no shared daughters + if (t1.protonIndex() == t2.protonIndex()) + continue; + if (t1.pionIndex() == t2.pionIndex()) + continue; + if (t1.protonIndex() == t2.pionIndex()) + continue; + if (t1.pionIndex() == t2.protonIndex()) + continue; + + const int status = static_cast(t1.v0Status()); + if (status < 0 || status >= nStat) + continue; + + const int ptB0 = mb.ptBin(t1.lambdaPt()); + const int etaB0 = mb.etaBin(t1.lambdaEta()); + const int phiB0 = mb.phiBin(phi0To2Pi(t1.lambdaPhi())); + if (ptB0 < 0 || etaB0 < 0 || phiB0 < 0) + continue; + + collectNeighborBins1D(ptB0, nPt, cfgV5NeighborPt.value, ptBins); + collectNeighborBins1D(etaB0, nEta, cfgV5NeighborEta.value, etaBins); + collectNeighborBinsPhi(phiB0, nPhi, cfgV5NeighborPhi.value, phiBins); + + matches.clear(); + const int64_t curColIdx = static_cast(col1.index()); + + // SE pair 4-vectors (Lambda only; consistent with your fillHistograms) + const auto l1SE = ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()); + const auto l2SE = ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()); + + for (int ptB : ptBins) { + for (int etaB : etaBins) { + for (int phiB : phiBins) { + auto const& vec = buffer[key5(colBin, status, ptB, etaB, phiB)]; + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) + continue; // different event + + auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); + + if (tX.globalIndex() == t2.globalIndex()) + continue; + + if (!selectionV0(tX)) + continue; + const auto l1ME = ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()); + const auto l2ME = l2SE; // fixed second leg from SE + + if (!passPairClosenessV5(l1SE, l2SE, l1ME, l2ME, + cfgV5MaxAbsDeltaPairPt.value, + cfgV5MaxAbsDeltaPairEta.value, + cfgV5MaxAbsDeltaPairPhi.value)) { + continue; + } + + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } + + if (matches.empty()) + continue; + + // dedupe + std::sort(matches.begin(), matches.end(), + [](auto const& a, auto const& b) { return std::tie(a.collisionIdx, a.rowIndex) < std::tie(b.collisionIdx, b.rowIndex); }); + matches.erase(std::unique(matches.begin(), matches.end(), + [](auto const& a, auto const& b) { return a.collisionIdx == b.collisionIdx && a.rowIndex == b.rowIndex; }), + matches.end()); + if (matches.empty()) + continue; + + // optional cap + if (cfgV5MaxMatches.value > 0 && (int)matches.size() > cfgV5MaxMatches.value) { + matches.resize(cfgV5MaxMatches.value); + } + + const float wBase = 1.0f / static_cast(matches.size()); + + for (auto const& m : matches) { + auto tX = V0s.iteratorAt(static_cast(m.rowIndex)); + + auto proton = ROOT::Math::PtEtaPhiMVector(tX.protonPt(), tX.protonEta(), tX.protonPhi(), o2::constants::physics::MassProton); + auto lambda = ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()); + auto proton2 = ROOT::Math::PtEtaPhiMVector(t2.protonPt(), t2.protonEta(), t2.protonPhi(), o2::constants::physics::MassProton); + auto lambda2 = ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()); + + const float dPhi = deltaPhiMinusPiToPi(static_cast(lambda.Phi()), static_cast(lambda2.Phi())); + histos.fill(HIST("deltaPhiMix"), dPhi, wBase); + fillHistograms(tX.v0Status(), t2.v0Status(), lambda, lambda2, proton, proton2, /*datatype=*/1, /*mixpairweight=*/wBase); + } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMEV5, "Process data ME v5 (paper-style)", false); + + void processMCMEV5(EventCandidatesMC const& collisions, AllTrackCandidatesMC const& V0sMC) + { + MixBinner mb{ptMin.value, ptMax.value, ptMix.value, + v0eta.value, etaMix.value, + phiMix.value}; + + const int nCol = colBinning.getAllBinsCount(); + const int nStat = N_STATUS; + const int nPt = mb.nPt(); + const int nEta = mb.nEta(); + const int nPhi = mb.nPhi(); + + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi; + std::vector> buffer(nKeys); + + auto key5 = [&](int colBin, int status, int ptB, int etaB, int phiB) -> size_t { + size_t k = static_cast(colBin); + k = k * static_cast(nStat) + static_cast(status); + k = k * static_cast(nPt) + static_cast(ptB); + k = k * static_cast(nEta) + static_cast(etaB); + k = k * static_cast(nPhi) + static_cast(phiB); + return k; + }; + + // PASS 1: buffer fill + for (auto const& col : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col), mcacc::cent(col))); + auto slice = V0sMC.sliceBy(tracksPerCollisionV0mc, col.index()); + + for (auto const& t : slice) { + if (!selectionV0MC(t)) + continue; + + const int status = mcacc::v0Status(t); + if (status < 0 || status >= nStat) + continue; + + const int ptB = mb.ptBin(mcacc::lamPt(t)); + const int etaB = mb.etaBin(mcacc::lamEta(t)); + const int phiB = mb.phiBin(phi0To2Pi(mcacc::lamPhi(t))); + if (ptB < 0 || etaB < 0 || phiB < 0) + continue; + + buffer[key5(colBin, status, ptB, etaB, phiB)].push_back(BufferCand{ + .collisionIdx = static_cast(col.index()), + .rowIndex = static_cast(t.globalIndex()), + .v0Status = static_cast(status), + .ptBin = static_cast(ptB), + .etaBin = static_cast(etaB), + .phiBin = static_cast(phiB), + .mBin = 0}); + } + } + + std::vector ptBins, etaBins, phiBins; + std::vector matches; + matches.reserve(256); + + // PASS 2: build ME + for (auto const& col1 : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col1), mcacc::cent(col1))); + auto poolA = V0sMC.sliceBy(tracksPerCollisionV0mc, col1.index()); + + for (auto const& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { + if (!selectionV0MC(t1) || !selectionV0MC(t2)) + continue; + if (t2.index() <= t1.index()) + continue; + + // no shared daughters + if (mcacc::prIdx(t1) == mcacc::prIdx(t2)) + continue; + if (mcacc::piIdx(t1) == mcacc::piIdx(t2)) + continue; + if (mcacc::prIdx(t1) == mcacc::piIdx(t2)) + continue; + if (mcacc::piIdx(t1) == mcacc::prIdx(t2)) + continue; + + const int status = mcacc::v0Status(t1); + if (status < 0 || status >= nStat) + continue; + + const int ptB0 = mb.ptBin(mcacc::lamPt(t1)); + const int etaB0 = mb.etaBin(mcacc::lamEta(t1)); + const int phiB0 = mb.phiBin(phi0To2Pi(mcacc::lamPhi(t1))); + if (ptB0 < 0 || etaB0 < 0 || phiB0 < 0) + continue; + + collectNeighborBins1D(ptB0, nPt, cfgV5NeighborPt.value, ptBins); + collectNeighborBins1D(etaB0, nEta, cfgV5NeighborEta.value, etaBins); + collectNeighborBinsPhi(phiB0, nPhi, cfgV5NeighborPhi.value, phiBins); + + matches.clear(); + const int64_t curColIdx = static_cast(col1.index()); + + const auto l1SE = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(t1), mcacc::lamEta(t1), mcacc::lamPhi(t1), mcacc::lamMass(t1)); + const auto l2SE = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(t2), mcacc::lamEta(t2), mcacc::lamPhi(t2), mcacc::lamMass(t2)); + + for (int ptB : ptBins) { + for (int etaB : etaBins) { + for (int phiB : phiBins) { + auto const& vec = buffer[key5(colBin, status, ptB, etaB, phiB)]; + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) + continue; + auto tX = V0sMC.iteratorAt(static_cast(bc.rowIndex)); + + if (tX.globalIndex() == t2.globalIndex()) + continue; + + if (!selectionV0MC(tX)) + continue; + const auto l1ME = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(tX), mcacc::lamEta(tX), mcacc::lamPhi(tX), mcacc::lamMass(tX)); + const auto l2ME = l2SE; + if (!passPairClosenessV5(l1SE, l2SE, l1ME, l2ME, cfgV5MaxAbsDeltaPairPt.value, cfgV5MaxAbsDeltaPairEta.value, cfgV5MaxAbsDeltaPairPhi.value)) { + continue; + } + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } + if (matches.empty()) + continue; + std::sort(matches.begin(), matches.end(), + [](auto const& a, auto const& b) { return std::tie(a.collisionIdx, a.rowIndex) < std::tie(b.collisionIdx, b.rowIndex); }); + matches.erase(std::unique(matches.begin(), matches.end(), + [](auto const& a, auto const& b) { return a.collisionIdx == b.collisionIdx && a.rowIndex == b.rowIndex; }), + matches.end()); + if (matches.empty()) + continue; + + if (cfgV5MaxMatches.value > 0 && (int)matches.size() > cfgV5MaxMatches.value) { + matches.resize(cfgV5MaxMatches.value); + } + const float wBase = 1.0f / static_cast(matches.size()); + for (auto const& m : matches) { + auto tX = V0sMC.iteratorAt(static_cast(m.rowIndex)); + auto pX = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(tX), mcacc::prEta(tX), mcacc::prPhi(tX), o2::constants::physics::MassProton); + auto lX = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(tX), mcacc::lamEta(tX), mcacc::lamPhi(tX), mcacc::lamMass(tX)); + auto p2 = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(t2), mcacc::prEta(t2), mcacc::prPhi(t2), o2::constants::physics::MassProton); + auto l2 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(t2), mcacc::lamEta(t2), mcacc::lamPhi(t2), mcacc::lamMass(t2)); + const float dPhi = deltaPhiMinusPiToPi(static_cast(lX.Phi()), static_cast(l2.Phi())); + histos.fill(HIST("deltaPhiMix"), dPhi, wBase); + fillHistograms(mcacc::v0Status(tX), mcacc::v0Status(t2), lX, l2, pX, p2, /*datatype=*/1, /*mixpairweight=*/wBase); + } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMCMEV5, "Process MC ME v5 (paper-style)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {