From 6aa2426038a061277e9f639f1b44cc54d50301e2 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Fri, 30 Jan 2026 18:52:54 +0100 Subject: [PATCH 01/11] 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 eb563fd9970..a829321f313 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -959,6 +959,7 @@ struct f1protoncorrelation { { const float maxMomPi = maxMomentumPion; const float maxMomK = maxMomentumKaon; + // const float pTofPiMin = momentumTOFPionMin; // const float pTofPiMax = momentumTOFPionMax; // const float pTofKMin = momentumTOFKaonMin; From c893cb06d99619d7e96a8a2cf7578c60bf2fb26e Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sun, 1 Feb 2026 11:04:28 +0100 Subject: [PATCH 02/11] 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 a829321f313..eb563fd9970 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -959,7 +959,6 @@ struct f1protoncorrelation { { const float maxMomPi = maxMomentumPion; const float maxMomK = maxMomentumKaon; - // const float pTofPiMin = momentumTOFPionMin; // const float pTofPiMax = momentumTOFPionMax; // const float pTofKMin = momentumTOFKaonMin; From 2eb1d315d6ddbad667a656b139d22e894602eae5 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Wed, 25 Feb 2026 19:58:49 +0100 Subject: [PATCH 03/11] 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 eb563fd9970..b2434b8fa69 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -1063,7 +1063,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 6a47f6543ca62a2da282ab82b466a013924f2c41 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 18:33:50 +0100 Subject: [PATCH 04/11] Add Event loss in phi meson task and add MC process function for spincorrelation task --- .../Strangeness/lambdaspincorrderived.cxx | 55 ++++++++++++++++--- 1 file changed, 46 insertions(+), 9 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 0ff87544f97..29ab35cc8b4 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1215,19 +1215,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; @@ -1235,17 +1252,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()) { From d72820332eca4934470a977a6b0ef003bd0eba49 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 20:29:09 +0100 Subject: [PATCH 05/11] 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 29ab35cc8b4..c5ca79f145b 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1215,36 +1215,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; @@ -1252,7 +1235,6 @@ struct lambdaspincorrderived { if (!checkKinematics(t1, tX)) { continue; } - matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); } }; From 398055cb9753662c1571125e96bf3659a75f39ac Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 2 Mar 2026 20:55:25 +0100 Subject: [PATCH 06/11] 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 b2434b8fa69..eb563fd9970 100644 --- a/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx +++ b/PWGLF/Tasks/Resonances/f1protoncorrelation.cxx @@ -1063,7 +1063,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 b1d564e55114a8b6fb588de6366bc537167cc79d Mon Sep 17 00:00:00 2001 From: skundu692 Date: Tue, 3 Mar 2026 21:19:20 +0100 Subject: [PATCH 07/11] Add new process function for mixing --- PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index c5ca79f145b..36e4b8a4917 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -204,6 +204,7 @@ struct lambdaspincorrderived { Configurable ConfWeightPathLAL2{"ConfWeightPathLAL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; Configurable ConfWeightPathALL2{"ConfWeightPathALL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; + // Mixing ///////// Configurable cfgV5NeighborPt{"cfgV5NeighborPt", 0, "v5: neighbor bins in pT (use symmetric ±N, edge-safe)"}; @@ -1800,6 +1801,7 @@ struct lambdaspincorrderived { out.erase(std::unique(out.begin(), out.end()), out.end()); } + static inline void collectNeighborBinsClamp(int b, int nBins, int nNeighbor, std::vector& out) { out.clear(); @@ -1940,7 +1942,6 @@ struct lambdaspincorrderived { continue; // same-event ordering } - // no shared daughters (same-event) if (t1.protonIndex() == t2.protonIndex()) continue; if (t1.pionIndex() == t2.pionIndex()) @@ -2184,6 +2185,7 @@ struct lambdaspincorrderived { continue; const int status = mcacc::v0Status(t1); + if (status < 0 || status >= nStat) { continue; } From a390d68bd5d1329c67d9c3fd492e82206d0069e1 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Thu, 5 Mar 2026 18:19:17 +0100 Subject: [PATCH 08/11] Add new process function for spin correlation --- .../Strangeness/lambdaspincorrderived.cxx | 29 ++++--------------- 1 file changed, 5 insertions(+), 24 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 36e4b8a4917..38c1b575da0 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1242,30 +1242,11 @@ struct lambdaspincorrderived { // 1) nominal φ-bin collectFrom(phiB); - // 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; - } + // 2) wrap only at boundaries: 0 <-> nPhi-1 + if (phiB == 0) { + collectFrom(nPhi - 1); + } else if (phiB == nPhi - 1) { + collectFrom(0); } if (matches.empty()) { From 70427e4299177684c0068571cccf08970b07330e Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sat, 7 Mar 2026 13:16:21 +0100 Subject: [PATCH 09/11] Fix event mixing sliding for spin correlation --- PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 38c1b575da0..432053c8cb0 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -1923,6 +1923,7 @@ struct lambdaspincorrderived { continue; // same-event ordering } + // no shared daughters (same-event) if (t1.protonIndex() == t2.protonIndex()) continue; if (t1.pionIndex() == t2.pionIndex()) @@ -2166,7 +2167,6 @@ struct lambdaspincorrderived { continue; const int status = mcacc::v0Status(t1); - if (status < 0 || status >= nStat) { continue; } From b14349e9dceb139fb6cfaa727d61239cdf3cd871 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sat, 7 Mar 2026 13:24:50 +0100 Subject: [PATCH 10/11] Apply clang-format --- PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx | 2 -- 1 file changed, 2 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 432053c8cb0..0ff87544f97 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -204,7 +204,6 @@ struct lambdaspincorrderived { Configurable ConfWeightPathLAL2{"ConfWeightPathLAL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; Configurable ConfWeightPathALL2{"ConfWeightPathALL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; - // Mixing ///////// Configurable cfgV5NeighborPt{"cfgV5NeighborPt", 0, "v5: neighbor bins in pT (use symmetric ±N, edge-safe)"}; @@ -1782,7 +1781,6 @@ struct lambdaspincorrderived { out.erase(std::unique(out.begin(), out.end()), out.end()); } - static inline void collectNeighborBinsClamp(int b, int nBins, int nNeighbor, std::vector& out) { out.clear(); From f6db8144a722e48b228757df5e19d28c8f33fec0 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Fri, 13 Mar 2026 21:35:08 +0100 Subject: [PATCH 11/11] Add new process function for mixing --- .../Strangeness/lambdaspincorrderived.cxx | 774 ++++++++++++++++-- 1 file changed, 724 insertions(+), 50 deletions(-) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 0ff87544f97..4e9bbcbfaa8 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -205,11 +205,12 @@ struct lambdaspincorrderived { Configurable ConfWeightPathALL2{"ConfWeightPathALL2", "Users/s/skundu/My/Object/spincorr/cent010LL", "Weight path 2"}; // Mixing ///////// + Configurable cfgV5MassBins{"cfgV5MassBins", 5, "Number of fixed mass bins for V5 mixing"}; 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 usePairKineMatch{"usePairKineMatch", true, "Require pair-level matching between (A,B) and (C,B)"}; Configurable cfgV5MaxMatches{"cfgV5MaxMatches", 50, "v5: max ME replacements per SE pair (after all cuts)"}; Configurable cfgMixSeed{"cfgMixSeed", 0xdecafbadULL, "RNG seed for downsampling matches (deterministic)"}; Configurable centMin{"centMin", 0, "Minimum Centrality"}; @@ -453,6 +454,88 @@ struct lambdaspincorrderived { return true; } + template + bool checkPairKinematics(TA const& A, TB const& B, TC const& C) + { + if (!usePairKineMatch) { + return true; + } + + const auto lA = ROOT::Math::PtEtaPhiMVector(A.lambdaPt(), A.lambdaEta(), A.lambdaPhi(), A.lambdaMass()); + const auto lB = ROOT::Math::PtEtaPhiMVector(B.lambdaPt(), B.lambdaEta(), B.lambdaPhi(), B.lambdaMass()); + const auto lC = ROOT::Math::PtEtaPhiMVector(C.lambdaPt(), C.lambdaEta(), C.lambdaPhi(), C.lambdaMass()); + + // relative pT inside the pair: |pT1 - pT2| + const float dPtAB = std::abs(A.lambdaPt() - B.lambdaPt()); + const float dPtCB = std::abs(C.lambdaPt() - B.lambdaPt()); + if (std::abs(dPtAB - dPtCB) > ptMix) { + return false; + } + + // relative longitudinal kinematics: |Δy| or |Δη| + if (userapidity) { + const float dYAB = std::abs(lA.Rapidity() - lB.Rapidity()); + const float dYCB = std::abs(lC.Rapidity() - lB.Rapidity()); + if (std::abs(dYAB - dYCB) > etaMix) { + return false; + } + } else { + const float dEtaAB = std::abs(A.lambdaEta() - B.lambdaEta()); + const float dEtaCB = std::abs(C.lambdaEta() - B.lambdaEta()); + if (std::abs(dEtaAB - dEtaCB) > etaMix) { + return false; + } + } + + // relative azimuth inside the pair: |Δφ| + const float dPhiAB = std::abs(deltaPhiMinusPiToPi((float)A.lambdaPhi(), (float)B.lambdaPhi())); + const float dPhiCB = std::abs(deltaPhiMinusPiToPi((float)C.lambdaPhi(), (float)B.lambdaPhi())); + if (std::abs(dPhiAB - dPhiCB) > phiMix) { + return false; + } + + return true; + } + template + bool checkPairKinematicsMC(TA const& A, TB const& B, TC const& C) + { + if (!usePairKineMatch) { + return true; + } + + const auto lA = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(A), mcacc::lamEta(A), mcacc::lamPhi(A), mcacc::lamMass(A)); + const auto lB = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(B), mcacc::lamEta(B), mcacc::lamPhi(B), mcacc::lamMass(B)); + const auto lC = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(C), mcacc::lamEta(C), mcacc::lamPhi(C), mcacc::lamMass(C)); + + const float dPtAB = std::abs(mcacc::lamPt(A) - mcacc::lamPt(B)); + const float dPtCB = std::abs(mcacc::lamPt(C) - mcacc::lamPt(B)); + if (std::abs(dPtAB - dPtCB) > ptMix) { + return false; + } + + if (userapidity) { + const float dYAB = std::abs(lA.Rapidity() - lB.Rapidity()); + const float dYCB = std::abs(lC.Rapidity() - lB.Rapidity()); + if (std::abs(dYAB - dYCB) > etaMix) { + return false; + } + } else { + const float dEtaAB = std::abs(mcacc::lamEta(A) - mcacc::lamEta(B)); + const float dEtaCB = std::abs(mcacc::lamEta(C) - mcacc::lamEta(B)); + if (std::abs(dEtaAB - dEtaCB) > etaMix) { + return false; + } + } + + const float dPhiAB = std::abs(deltaPhiMinusPiToPi((float)mcacc::lamPhi(A), (float)mcacc::lamPhi(B))); + const float dPhiCB = std::abs(deltaPhiMinusPiToPi((float)mcacc::lamPhi(C), (float)mcacc::lamPhi(B))); + if (std::abs(dPhiAB - dPhiCB) > phiMix) { + return false; + } + + return true; + } + void fillHistograms(int tag1, int tag2, const ROOT::Math::PtEtaPhiMVector& particle1, const ROOT::Math::PtEtaPhiMVector& particle2, const ROOT::Math::PtEtaPhiMVector& daughpart1, const ROOT::Math::PtEtaPhiMVector& daughpart2, @@ -1036,29 +1119,50 @@ struct lambdaspincorrderived { static constexpr int N_STATUS = 2; // v0Status ∈ {0,1} struct MixBinner { - // constructed from the task's configurables; φ is assumed already constrained into [0, 2π) + // constructed from the task's configurables; φ is assumed already constrained upstream float ptMin, ptMax, ptStep; float etaMin, etaMax, etaStep; float phiMin, phiMax, phiStep; - // Mass binning: [1.09, 1.14) with 50 bins (1e-3 GeV/c^2) - static constexpr float mMin = 1.09f; - static constexpr float mMax = 1.14f; - static constexpr int nM_ = 1; - static constexpr float mStep = (mMax - mMin) / nM_; + // configurable fixed mass-binning for mixing buffer + float mMin, mMax, mStep; + int nM_; int nPt_, nEta_, nPhi_; MixBinner(float ptMin_, float ptMax_, float ptStep_, float etaAbsMax, float etaStep_, - float phiStep_) - : ptMin(ptMin_), ptMax(ptMax_), ptStep(ptStep_), etaMin(-etaAbsMax), etaMax(+etaAbsMax), etaStep(etaStep_), phiMin(-static_cast(TMath::Pi())), phiMax(+static_cast(TMath::Pi())), phiStep(phiStep_) - // : ptMin(ptMin_), ptMax(ptMax_), ptStep(ptStep_), etaMin(-etaAbsMax), etaMax(+etaAbsMax), etaStep(etaStep_), phiMin(0.f), phiMax(static_cast(2.0 * TMath::Pi())), phiStep(phiStep_) + float phiStep_, + float mMin_, float mMax_, int nMassBins_) + : ptMin(ptMin_), + ptMax(ptMax_), + ptStep(ptStep_), + etaMin(-etaAbsMax), + etaMax(+etaAbsMax), + etaStep(etaStep_), + phiMin(-static_cast(TMath::Pi())), + phiMax(+static_cast(TMath::Pi())), + phiStep(phiStep_), + mMin(mMin_), + mMax(mMax_), + mStep(0.f), + nM_(std::max(1, nMassBins_)), + nPt_(0), + nEta_(0), + nPhi_(0) + // If you want phi in [0, 2pi), use: + // : ... phiMin(0.f), phiMax(static_cast(2.0 * TMath::Pi())), ... { ptStep = (ptStep > 0.f ? ptStep : 0.1f); etaStep = (etaStep > 0.f ? etaStep : 0.1f); phiStep = (phiStep > 0.f ? phiStep : 0.1f); + if (!(mMax > mMin)) { + mMin = 1.09f; + mMax = 1.14f; + } + mStep = (mMax - mMin) / static_cast(nM_); + nPt_ = std::max(1, static_cast(std::floor((ptMax - ptMin) / ptStep + 0.5f))); nEta_ = std::max(1, static_cast(std::floor((etaMax - etaMin) / etaStep + 0.5f))); nPhi_ = std::max(1, static_cast(std::ceil((phiMax - phiMin) / phiStep))); @@ -1071,14 +1175,17 @@ struct lambdaspincorrderived { inline int binFromValue(float v, float vmin, float step, int nBins) const { - if (!std::isfinite(v)) + if (!std::isfinite(v) || !std::isfinite(vmin) || !std::isfinite(step) || step <= 0.f || nBins <= 0) { return -1; + } const float x = (v - vmin) / step; int b = static_cast(std::floor(x + 1e-6f)); - if (b < 0) + if (b < 0) { return -1; - if (b >= nBins) + } + if (b >= nBins) { b = nBins - 1; // clamp exact-top edge + } return b; } @@ -1128,12 +1235,11 @@ struct lambdaspincorrderived { // ===================== Main mixing (with mass-bin + random unique sampling) ===================== void processMEV4(EventCandidates const& collisions, AllTrackCandidates const& V0s) { - // Build binner from your existing configurables MixBinner mb{ - ptMin.value, ptMax.value, ptMix.value, // pT range & step - v0etaMixBuffer.value, etaMix.value, // |eta| max & step - phiMix.value // φ step; φ range fixed to [0, 2π) - }; + ptMin.value, ptMax.value, ptMix.value, + v0etaMixBuffer.value, etaMix.value, + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.value}; const int nCol = colBinning.getAllBinsCount(); // event-class bins (vz, centrality) const int nStat = N_STATUS; // 2 @@ -1235,6 +1341,9 @@ struct lambdaspincorrderived { if (!checkKinematics(t1, tX)) { continue; } + if (!checkPairKinematics(t1, t2, tX)) { + continue; + } matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); } }; @@ -1565,11 +1674,11 @@ struct lambdaspincorrderived { // ----------------------------------------------------- void processMCMEV4(EventCandidatesMC const& collisions, AllTrackCandidatesMC const& V0sMC) { - // Same binner as in data MEV4 MixBinner mb{ ptMin.value, ptMax.value, ptMix.value, v0etaMixBuffer.value, etaMix.value, - phiMix.value}; + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.value}; const int nCol = colBinning.getAllBinsCount(); const int nStat = N_STATUS; @@ -1673,6 +1782,9 @@ struct lambdaspincorrderived { if (!checkKinematicsMC(t1, tX)) { continue; } + if (!checkPairKinematicsMC(t1, t2, tX)) { + continue; + } matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); } }; @@ -1827,12 +1939,11 @@ struct lambdaspincorrderived { void processMEV5(EventCandidates const& collisions, AllTrackCandidates const& V0s) { - // Buffer binning: use v0etaMixBuffer as the max range, etaMix as the bin step MixBinner mb{ - ptMin.value, ptMax.value, ptMix.value, // pT range & step - v0etaMixBuffer.value, etaMix.value, // |eta| (or |y|) max & step - phiMix.value // phi step - }; + ptMin.value, ptMax.value, ptMix.value, + v0etaMixBuffer.value, etaMix.value, + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.value}; const int nCol = colBinning.getAllBinsCount(); const int nStat = N_STATUS; @@ -2063,26 +2174,19 @@ struct lambdaspincorrderived { MixBinner mb{ ptMin.value, ptMax.value, ptMix.value, v0etaMixBuffer.value, etaMix.value, - phiMix.value}; + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.value}; const int nCol = colBinning.getAllBinsCount(); const int nStat = N_STATUS; const int nPt = mb.nPt(); - const int nEta = mb.nEta(); + const int nEta = mb.nEta(); // logical "nY" if userapidity=true const int nPhi = mb.nPhi(); + const int nM = mb.nM(); - const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi; + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; 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: fill buffer -------- for (auto const& col : collisions) { const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col), mcacc::cent(col))); @@ -2112,18 +2216,21 @@ struct lambdaspincorrderived { } const int phiB = mb.phiBin(RecoDecay::constrainAngle(mcacc::lamPhi(t), -TMath::Pi(), harmonic)); - if (ptB < 0 || etaB < 0 || phiB < 0) { + const int mB = mb.massBin(mcacc::lamMass(t)); + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 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}); + buffer[linearKey(colBin, status, ptB, etaB, phiB, mB, + nStat, nPt, nEta, nPhi, nM)] + .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 = static_cast(mB)}); } } @@ -2178,10 +2285,10 @@ struct lambdaspincorrderived { etaB = mb.etaBin(lv1.Rapidity()); } const int phiB = mb.phiBin(RecoDecay::constrainAngle(mcacc::lamPhi(t1), -TMath::Pi(), harmonic)); - if (ptB < 0 || etaB < 0 || phiB < 0) { + const int mB = mb.massBin(mcacc::lamMass(t1)); + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) { continue; } - collectNeighborBinsClamp(ptB, nPt, nN_pt, ptBins); collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBins); collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBins); @@ -2191,8 +2298,8 @@ struct lambdaspincorrderived { for (int ptUse : ptBins) { for (int etaUse : etaBins) { for (int phiUse : phiBins) { - auto const& vec = buffer[key5(colBin, status, ptUse, etaUse, phiUse)]; - + auto const& vec = buffer[linearKey(colBin, status, ptUse, etaUse, phiUse, mB, + nStat, nPt, nEta, nPhi, nM)]; for (auto const& bc : vec) { if (bc.collisionIdx == curColIdx) { continue; @@ -2285,6 +2392,573 @@ struct lambdaspincorrderived { } } PROCESS_SWITCH(lambdaspincorrderived, processMCMEV5, "Process MC ME v5 (paper-style)", false); + + void processMEV6(EventCandidates const& collisions, AllTrackCandidates const& V0s) + { + MixBinner mb{ + ptMin.value, ptMax.value, ptMix.value, + v0etaMixBuffer.value, etaMix.value, + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.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 int nM = mb.nM(); + + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; + std::vector> buffer(nKeys); + + // -------- PASS 1: fill buffer -------- + for (auto const& col : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(col.posz(), col.cent())); + if (colBin < 0) { + continue; + } + + 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()); + + int etaB = mb.etaBin(t.lambdaEta()); + if (userapidity) { + const auto lv = ROOT::Math::PtEtaPhiMVector(t.lambdaPt(), t.lambdaEta(), t.lambdaPhi(), t.lambdaMass()); + etaB = mb.etaBin(lv.Rapidity()); + } + + const int phiB = mb.phiBin(RecoDecay::constrainAngle(t.lambdaPhi(), -TMath::Pi(), harmonic)); + const int mB = mb.massBin(t.lambdaMass()); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) { + continue; + } + + const size_t key = linearKey(colBin, status, ptB, etaB, phiB, mB, + nStat, nPt, nEta, nPhi, nM); + + buffer[key].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 = static_cast(mB)}); + } + } + + const int nN_pt = std::max(0, cfgV5NeighborPt.value); + const int nN_eta = std::max(0, cfgV5NeighborEta.value); + const int nN_phi = std::max(0, cfgV5NeighborPhi.value); + + std::vector ptBins, etaBins, phiBins; + std::vector matches1, matches2; + matches1.reserve(256); + matches2.reserve(256); + + auto collectMatchesForReplacedLeg = [&](auto const& tRep, auto const& tKeep, int colBin, int64_t curColIdx, std::vector& matches) { + matches.clear(); + + const int status = static_cast(tRep.v0Status()); + if (status < 0 || status >= nStat) { + return; + } + + const int ptB = mb.ptBin(tRep.lambdaPt()); + + int etaB = mb.etaBin(tRep.lambdaEta()); + if (userapidity) { + const auto lv = ROOT::Math::PtEtaPhiMVector(tRep.lambdaPt(), tRep.lambdaEta(), tRep.lambdaPhi(), tRep.lambdaMass()); + etaB = mb.etaBin(lv.Rapidity()); + } + + const int phiB = mb.phiBin(RecoDecay::constrainAngle(tRep.lambdaPhi(), -TMath::Pi(), harmonic)); + const int mB = mb.massBin(tRep.lambdaMass()); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) { + return; + } + + collectNeighborBinsClamp(ptB, nPt, nN_pt, ptBins); + collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBins); + collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBins); + + for (int ptUse : ptBins) { + for (int etaUse : etaBins) { + for (int phiUse : phiBins) { + const auto& vec = buffer[linearKey(colBin, status, ptUse, etaUse, phiUse, mB, + nStat, nPt, nEta, nPhi, nM)]; + + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) { + continue; + } + + auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); + + if (!selectionV0(tX)) { + continue; + } + if (!checkKinematics(tRep, tX)) { + continue; + } + if (!checkPairKinematics(tRep, tKeep, tX)) { + continue; + } + + if (tX.globalIndex() == tRep.globalIndex()) { + continue; + } + if (tX.globalIndex() == tKeep.globalIndex()) { + continue; + } + + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } + + 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()); + }; + + auto downsampleMatches = [&](std::vector& matches, uint64_t seedBase) { + if (cfgV5MaxMatches.value > 0 && (int)matches.size() > cfgV5MaxMatches.value) { + uint64_t seed = cfgMixSeed.value ^ splitmix64(seedBase); + const int K = cfgV5MaxMatches.value; + for (int i = 0; i < K; ++i) { + seed = splitmix64(seed); + const int j = i + (int)(seed % (uint64_t)(matches.size() - i)); + std::swap(matches[i], matches[j]); + } + matches.resize(K); + } + }; + + // -------- PASS 2: two-leg mixing -------- + for (auto const& col1 : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); + if (colBin < 0) { + continue; + } + + const int64_t curColIdx = static_cast(col1.index()); + 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; + } + + if (t1.protonIndex() == t2.protonIndex()) { + continue; + } + if (t1.pionIndex() == t2.pionIndex()) { + continue; + } + if (t1.protonIndex() == t2.pionIndex()) { + continue; + } + if (t1.pionIndex() == t2.protonIndex()) { + continue; + } + + // leg 1 replaced: (t1,t2) -> (tX,t2) + collectMatchesForReplacedLeg(t1, t2, colBin, curColIdx, matches1); + + // leg 2 replaced: (t1,t2) -> (t1,tY) + collectMatchesForReplacedLeg(t2, t1, colBin, curColIdx, matches2); + + downsampleMatches(matches1, + (uint64_t)t1.globalIndex() ^ (splitmix64((uint64_t)t2.globalIndex()) + 0x111ULL) ^ splitmix64((uint64_t)curColIdx)); + downsampleMatches(matches2, + (uint64_t)t2.globalIndex() ^ (splitmix64((uint64_t)t1.globalIndex()) + 0x222ULL) ^ splitmix64((uint64_t)curColIdx)); + + const int nReuse = static_cast(matches1.size() + matches2.size()); + if (nReuse <= 0) { + continue; + } + + const float wSE = 1.0f / static_cast(nReuse); + + // replace t1 -> tX, keep t2 + for (auto const& m : matches1) { + 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 int ptype = pairTypeCode(tX.v0Status(), t2.v0Status()); + double centPairWeight = 1.0; + if (hweightCentPair) { + const int bin = hweightCentPair->FindBin(col1.cent(), ptype); + centPairWeight = hweightCentPair->GetBinContent(bin); + if (centPairWeight <= 0.0) { + centPairWeight = 1.0; + } + } + + const float meWeight = wSE * centPairWeight; + const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); + histos.fill(HIST("deltaPhiMix"), dPhi, wSE); + histos.fill(HIST("hCentPairTypeME"), col1.cent(), ptype, wSE); + + fillHistograms(tX.v0Status(), t2.v0Status(), + lambda, lambda2, proton, proton2, + 1, meWeight); + } + + // replace t2 -> tY, keep t1 + for (auto const& m : matches2) { + auto tY = V0s.iteratorAt(static_cast(m.rowIndex)); + + auto proton = ROOT::Math::PtEtaPhiMVector(t1.protonPt(), t1.protonEta(), t1.protonPhi(), + o2::constants::physics::MassProton); + auto lambda = ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), + t1.lambdaMass()); + + auto proton2 = ROOT::Math::PtEtaPhiMVector(tY.protonPt(), tY.protonEta(), tY.protonPhi(), + o2::constants::physics::MassProton); + auto lambda2 = ROOT::Math::PtEtaPhiMVector(tY.lambdaPt(), tY.lambdaEta(), tY.lambdaPhi(), + tY.lambdaMass()); + + const int ptype = pairTypeCode(t1.v0Status(), tY.v0Status()); + double centPairWeight = 1.0; + if (hweightCentPair) { + const int bin = hweightCentPair->FindBin(col1.cent(), ptype); + centPairWeight = hweightCentPair->GetBinContent(bin); + if (centPairWeight <= 0.0) { + centPairWeight = 1.0; + } + } + + const float meWeight = wSE * centPairWeight; + const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); + histos.fill(HIST("deltaPhiMix"), dPhi, wSE); + histos.fill(HIST("hCentPairTypeME"), col1.cent(), ptype, wSE); + + fillHistograms(t1.v0Status(), tY.v0Status(), + lambda, lambda2, proton, proton2, + 1, meWeight); + } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMEV6, "Process data ME v6 two-leg", false); + + void processMCMEV6(EventCandidatesMC const& collisions, AllTrackCandidatesMC const& V0sMC) + { + MixBinner mb{ + ptMin.value, ptMax.value, ptMix.value, + v0etaMixBuffer.value, etaMix.value, + phiMix.value, + MassMin.value, MassMax.value, cfgV5MassBins.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 int nM = mb.nM(); + + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; + std::vector> buffer(nKeys); + + // -------- PASS 1: fill buffer -------- + for (auto const& col : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col), mcacc::cent(col))); + if (colBin < 0) { + continue; + } + + 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)); + + int etaB = mb.etaBin(mcacc::lamEta(t)); + if (userapidity) { + const auto lv = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(t), mcacc::lamEta(t), mcacc::lamPhi(t), mcacc::lamMass(t)); + etaB = mb.etaBin(lv.Rapidity()); + } + + const int phiB = mb.phiBin(RecoDecay::constrainAngle(mcacc::lamPhi(t), -TMath::Pi(), harmonic)); + const int mB = mb.massBin(mcacc::lamMass(t)); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) { + continue; + } + + const size_t key = linearKey(colBin, status, ptB, etaB, phiB, mB, + nStat, nPt, nEta, nPhi, nM); + + buffer[key].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 = static_cast(mB)}); + } + } + + const int nN_pt = std::max(0, cfgV5NeighborPt.value); + const int nN_eta = std::max(0, cfgV5NeighborEta.value); + const int nN_phi = std::max(0, cfgV5NeighborPhi.value); + + std::vector ptBins, etaBins, phiBins; + std::vector matches1, matches2; + matches1.reserve(256); + matches2.reserve(256); + + auto collectMatchesForReplacedLeg = [&](auto const& tRep, auto const& tKeep, int colBin, int64_t curColIdx, std::vector& matches) { + matches.clear(); + + const int status = mcacc::v0Status(tRep); + if (status < 0 || status >= nStat) { + return; + } + + const int ptB = mb.ptBin(mcacc::lamPt(tRep)); + + int etaB = mb.etaBin(mcacc::lamEta(tRep)); + if (userapidity) { + const auto lv = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(tRep), mcacc::lamEta(tRep), mcacc::lamPhi(tRep), mcacc::lamMass(tRep)); + etaB = mb.etaBin(lv.Rapidity()); + } + + const int phiB = mb.phiBin(RecoDecay::constrainAngle(mcacc::lamPhi(tRep), -TMath::Pi(), harmonic)); + const int mB = mb.massBin(mcacc::lamMass(tRep)); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) { + return; + } + + collectNeighborBinsClamp(ptB, nPt, nN_pt, ptBins); + collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBins); + collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBins); + + for (int ptUse : ptBins) { + for (int etaUse : etaBins) { + for (int phiUse : phiBins) { + const auto& vec = buffer[linearKey(colBin, status, ptUse, etaUse, phiUse, mB, + nStat, nPt, nEta, nPhi, nM)]; + + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) { + continue; + } + + auto tX = V0sMC.iteratorAt(static_cast(bc.rowIndex)); + + if (!selectionV0MC(tX)) { + continue; + } + if (!checkKinematicsMC(tRep, tX)) { + continue; + } + if (!checkPairKinematicsMC(tRep, tKeep, tX)) { + continue; + } + + if (tX.globalIndex() == tRep.globalIndex()) { + continue; + } + if (tX.globalIndex() == tKeep.globalIndex()) { + continue; + } + + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } + + 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()); + }; + + auto downsampleMatches = [&](std::vector& matches, uint64_t seedBase) { + if (cfgV5MaxMatches.value > 0 && (int)matches.size() > cfgV5MaxMatches.value) { + uint64_t seed = cfgMixSeed.value ^ splitmix64(seedBase); + const int K = cfgV5MaxMatches.value; + for (int i = 0; i < K; ++i) { + seed = splitmix64(seed); + const int j = i + (int)(seed % (uint64_t)(matches.size() - i)); + std::swap(matches[i], matches[j]); + } + matches.resize(K); + } + }; + + // -------- PASS 2: two-leg mixing -------- + for (auto const& col1 : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(mcacc::posz(col1), mcacc::cent(col1))); + if (colBin < 0) { + continue; + } + + const int64_t curColIdx = static_cast(col1.index()); + 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; + } + + 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; + } + + collectMatchesForReplacedLeg(t1, t2, colBin, curColIdx, matches1); + collectMatchesForReplacedLeg(t2, t1, colBin, curColIdx, matches2); + + downsampleMatches(matches1, + (uint64_t)t1.globalIndex() ^ (splitmix64((uint64_t)t2.globalIndex()) + 0x111ULL) ^ splitmix64((uint64_t)curColIdx)); + downsampleMatches(matches2, + (uint64_t)t2.globalIndex() ^ (splitmix64((uint64_t)t1.globalIndex()) + 0x222ULL) ^ splitmix64((uint64_t)curColIdx)); + + const int nReuse = static_cast(matches1.size() + matches2.size()); + if (nReuse <= 0) { + continue; + } + + const float wSE = 1.0f / static_cast(nReuse); + + // replace t1 -> tX, keep t2 + for (auto const& m : matches1) { + 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 int ptype = pairTypeCode(mcacc::v0Status(tX), mcacc::v0Status(t2)); + double centPairWeight = 1.0; + if (hweightCentPair) { + const int bin = hweightCentPair->FindBin(mcacc::cent(col1), ptype); + centPairWeight = hweightCentPair->GetBinContent(bin); + if (centPairWeight <= 0.0) { + centPairWeight = 1.0; + } + } + + const float meWeight = wSE * centPairWeight; + const float dPhi = deltaPhiMinusPiToPi((float)lX.Phi(), (float)l2.Phi()); + histos.fill(HIST("deltaPhiMix"), dPhi, wSE); + histos.fill(HIST("hCentPairTypeME"), mcacc::cent(col1), ptype, wSE); + + fillHistograms(mcacc::v0Status(tX), mcacc::v0Status(t2), + lX, l2, pX, p2, + 1, meWeight); + } + + // replace t2 -> tY, keep t1 + for (auto const& m : matches2) { + auto tY = V0sMC.iteratorAt(static_cast(m.rowIndex)); + + auto p1 = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(t1), mcacc::prEta(t1), mcacc::prPhi(t1), + o2::constants::physics::MassProton); + auto l1 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(t1), mcacc::lamEta(t1), mcacc::lamPhi(t1), + mcacc::lamMass(t1)); + + auto pY = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(tY), mcacc::prEta(tY), mcacc::prPhi(tY), + o2::constants::physics::MassProton); + auto lY = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(tY), mcacc::lamEta(tY), mcacc::lamPhi(tY), + mcacc::lamMass(tY)); + + const int ptype = pairTypeCode(mcacc::v0Status(t1), mcacc::v0Status(tY)); + double centPairWeight = 1.0; + if (hweightCentPair) { + const int bin = hweightCentPair->FindBin(mcacc::cent(col1), ptype); + centPairWeight = hweightCentPair->GetBinContent(bin); + if (centPairWeight <= 0.0) { + centPairWeight = 1.0; + } + } + + const float meWeight = wSE * centPairWeight; + const float dPhi = deltaPhiMinusPiToPi((float)l1.Phi(), (float)lY.Phi()); + histos.fill(HIST("deltaPhiMix"), dPhi, wSE); + histos.fill(HIST("hCentPairTypeME"), mcacc::cent(col1), ptype, wSE); + + fillHistograms(mcacc::v0Status(t1), mcacc::v0Status(tY), + l1, lY, p1, pY, + 1, meWeight); + } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMCMEV6, "Process MC ME v6 two-leg", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {