diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index 668dd8f7541..1e967865e42 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -43,6 +43,7 @@ #include #include +#include #include #include // for std::abs #include @@ -315,8 +316,127 @@ struct lambdaspincorrderived { HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + // ------------------------- + // Systematics control + // ------------------------- + + struct SysFromDeltaPtDeltaEtaDeltaPhi { + // Mixed -event kinematic cut + float minDeltaPt; + float maxDelatPt; + + float minDeltaEta; + float maxDeltaEta; + + float minDeltaPhi; + float maxDeltaPhi; + }; + + struct SysCuts { + // Topology cut + float maxDcaDaughters; + float minRadius; + float maxRadius; + float mincosPA; + float mindcaProton; + float mindcaPion; + float maxdcaV0ToPV; + }; + + std::vector sysCuts; // sysCuts[0] = default; sysCuts[1..] random unique + int nSysTotal = 1; + + void buildSystematicCuts() + { + + // 2/3 options per cut: index 0 = DEFAULT, index 1/2 = variations + // Fill these with the exact values you want (I used your def as index 0 + your old options as 1/2) + const std::array optDcaDaughter{1.00f, 0.90f}; + const std::array optminRadius{1.20f, 1.00f}; + const std::array optmaxRadius{23.0f, 35.0f}; + const std::array optcosPA{0.999f, 0.995f}; + + const std::array optdcaProton{0.07f, 0.05f, 0.08f}; + const std::array optdcaPion{0.2f, 0.1f, 0.3f}; + const std::array optdcaV0ToPV{1.2f, 1.0f, 1.4f}; + + // Helper: build SysCuts from chosen indices (0..2) + auto buildFromIdx = [&](int i0, int i1, int i2, int i3, + int i4, int i5, int i6) -> SysCuts { + SysCuts c{}; + c.maxDcaDaughters = optDcaDaughter[i0]; + c.minRadius = optminRadius[i1]; + c.maxRadius = optmaxRadius[i2]; + c.mincosPA = optcosPA[i3]; + + c.mindcaProton = optdcaProton[i4]; + c.mindcaPion = optdcaPion[i5]; + c.maxdcaV0ToPV = optdcaV0ToPV[i6]; + return c; + }; + sysCuts.clear(); + // sysId=0 must be strict default (all indices = 0) + SysCuts def = buildFromIdx(0, 0, 0, 0, 0, 0, 0); + sysCuts.push_back(def); + // sysId = 1: + sysCuts.push_back(buildFromIdx(1, 0, 0, 0, 0, 0, 0)); + // sysId = 2: + sysCuts.push_back(buildFromIdx(0, 1, 0, 0, 0, 0, 0)); + // sysId = 3 + sysCuts.push_back(buildFromIdx(0, 0, 1, 0, 0, 0, 0)); + // sysId = 4 + sysCuts.push_back(buildFromIdx(0, 0, 0, 1, 0, 0, 0)); + + // sysId = 5 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 1, 0, 0)); + // sysId = 6 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 2, 0, 0)); + // sysId = 7 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 0, 1, 0)); + // sysId = 8 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 0, 2, 0)); + // sysId = 9 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 0, 0, 1)); + // sysId = 10 + sysCuts.push_back(buildFromIdx(0, 0, 0, 0, 0, 0, 2)); + + nSysTotal = (int)sysCuts.size(); + + // print all configurations + LOGF(info, "========================================"); + LOGF(info, "Total systematic configurations: %d", nSysTotal); + LOGF(info, "sysId=0: DEFAULT (all cuts at default)"); + LOGF(info, "sysId=1: maxDcaDaughters %.2f -> %.2f", + optDcaDaughter[0], optDcaDaughter[1]); + LOGF(info, "sysId=2: minRadius %.2f -> %.2f", + optminRadius[0], optminRadius[1]); + LOGF(info, "sysId=3: maxRadius %.2f -> %.2f", + optmaxRadius[0], optmaxRadius[1]); + LOGF(info, "sysId=4: mincosPA %.3f -> %.3f", + optcosPA[0], optcosPA[1]); + LOGF(info, "sysId=5: mindcaProton %.2f -> %.2f (tighter)", + optdcaProton[0], optdcaProton[1]); + LOGF(info, "sysId=6: mindcaProton %.2f -> %.2f (looser)", + optdcaProton[0], optdcaProton[2]); + LOGF(info, "sysId=7: mindcaPion %.2f -> %.2f (tighter)", + optdcaPion[0], optdcaPion[1]); + LOGF(info, "sysId=8: mindcaPion %.2f -> %.2f (looser)", + optdcaPion[0], optdcaPion[2]); + LOGF(info, "sysId=9: maxdcaV0ToPV %.2f -> %.2f (tighter)", + optdcaV0ToPV[0], optdcaV0ToPV[1]); + LOGF(info, "sysId=10: maxdcaV0ToPV %.2f -> %.2f (looser)", + optdcaV0ToPV[0], optdcaV0ToPV[2]); + LOGF(info, "========================================"); + } + void init(o2::framework::InitContext&) { + buildSystematicCuts(); + + nSysTotal = (int)sysCuts.size(); // or whatever vector you fill + LOGF(info, "sysCuts.size()=%zu nSysTotal=%d", sysCuts.size(), nSysTotal); + const AxisSpec thnAxisSys{nSysTotal, -0.5f, float(nSysTotal) - 0.5f, "sysId"}; + if (fillBasicQAHistos) { histos.add("hPtRadiusV0", "V0 QA;#it{p}_{T}^{V0} (GeV/#it{c});V0 decay radius (cm)", kTH2F, {{100, 0.0, 10.0}, {120, 0.0, 45.0}}); histos.add("hPtYSame", "hPtYSame", kTH2F, {{100, 0.0, 10.0}, {200, -1.0, 1.0}}); @@ -452,6 +572,17 @@ struct lambdaspincorrderived { histos.add("hSparseAntiLambdaLambdaMixed", "hSparseAntiLambdaLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true); histos.add("hSparseAntiLambdaAntiLambdaMixed", "hSparseAntiLambdaAntiLambdaMixed", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR}, true); + // Systematic THnSparse for analysis + histos.add("hSparseLambdaLambdaSys", "hSparseLambdaLambdaSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseLambdaAntiLambdaSys", "hSparseLambdaAntiLambdaSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseAntiLambdaLambdaSys", "hSparseAntiLambdLambdaSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseAntiLambdaAntiLambdaSys", "hSparseAntiLambdaAntiLambdaSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + + histos.add("hSparseLambdaLambdaMixedSys", "hSparseLambdaLambdaMixedSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseLambdaAntiLambdaMixedSys", "hSparseLambdaAntiLambdaMixedSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseAntiLambdaLambdaMixedSys", "hSparseAntiLambdaLambdaMixedSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + histos.add("hSparseAntiLambdaAntiLambdaMixedSys", "hSparseAntiLambdaAntiLambdaMixedSys", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisR, thnAxisSys}, true); + if (fillAnalysisSparses) { histos.add("hSparseLambdaLambdaAnalysis", "hSparseLambdaLambdaAnalysis", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisDeltaR, configThnAxisDeltaRap, configThnAxisDeltaPhi}, true); histos.add("hSparseLambdaAntiLambdaAnalysis", "hSparseLambdaAntiLambdaAnalysis", HistType::kTHnSparseF, {configThnAxisInvMass, configThnAxisInvMass, configThnAxisPol, configThnAxisDeltaR, configThnAxisDeltaRap, configThnAxisDeltaPhi}, true); @@ -622,6 +753,51 @@ struct lambdaspincorrderived { return true; } + template + bool selectionV0Sys(T1 const& candidate, T2 const& SysCuts) + { + auto particle = ROOT::Math::PtEtaPhiMVector(candidate.lambdaPt(), candidate.lambdaEta(), candidate.lambdaPhi(), candidate.lambdaMass()); + if (std::abs(particle.Rapidity()) > rapidity || std::abs(particle.Eta()) > v0eta) { + return false; + } + if (candidate.lambdaMass() < MassMin || candidate.lambdaMass() > MassMax) { + return false; + } + if (candidate.v0Cospa() < SysCuts.mincosPA) { + return false; + } + if (checkDoubleStatus && candidate.doubleStatus()) { + return false; + } + if (candidate.v0Radius() > SysCuts.maxRadius) { + return false; + } + if (candidate.v0Radius() < SysCuts.minRadius) { + return false; + } + if (candidate.dcaBetweenDaughter() > SysCuts.maxDcaDaughters) { + return false; + } + + if (candidate.dcaV0ToPV() > SysCuts.maxdcaV0ToPV) { + return false; + } + + if (candidate.v0Status() == 0 && (std::abs(candidate.dcaPositive()) < SysCuts.mindcaProton || std::abs(candidate.dcaNegative()) < SysCuts.mindcaPion)) { + return false; + } + if (candidate.v0Status() == 1 && (std::abs(candidate.dcaPositive()) < SysCuts.mindcaPion || std::abs(candidate.dcaNegative()) < SysCuts.mindcaProton)) { + return false; + } + if (candidate.lambdaPt() < ptMin) { + return false; + } + if (candidate.lambdaPt() > ptMax) { + return false; + } + return true; + } + template bool checkKinematics(T1 const& c1, T2 const& c2) { @@ -1066,478 +1242,758 @@ struct lambdaspincorrderived { } } } - static inline int pairTypeCode(int tag1, int tag2) + + void fillHistogramsSys(int tag1, int tag2, + const ROOT::Math::PtEtaPhiMVector& particle1, const ROOT::Math::PtEtaPhiMVector& particle2, + const ROOT::Math::PtEtaPhiMVector& daughpart1, const ROOT::Math::PtEtaPhiMVector& daughpart2, + int datatype, float mixpairweight, int SysId, int replacedLeg = 1, int weightMapLeg = -1) { - if (tag1 == 0 && tag2 == 0) { - return 0; // LL - } else if (tag1 == 0 && tag2 == 1) { - return 1; // LAL - } else if (tag1 == 1 && tag2 == 0) { - return 2; // ALL + + auto lambda1Mass = 0.0; + auto lambda2Mass = 0.0; + if (!usePDGM) { + lambda1Mass = particle1.M(); + lambda2Mass = particle2.M(); } else { - return 3; // ALAL + lambda1Mass = o2::constants::physics::MassLambda; + lambda2Mass = o2::constants::physics::MassLambda; } - } - template - static inline bool hasSharedDaughters(const A& a, const B& b) - { - return (a.protonIndex() == b.protonIndex()) || - (a.pionIndex() == b.pionIndex()) || - (a.protonIndex() == b.pionIndex()) || - (a.pionIndex() == b.protonIndex()); - } + auto particle1Dummy = ROOT::Math::PtEtaPhiMVector(particle1.Pt(), particle1.Eta(), particle1.Phi(), lambda1Mass); + auto particle2Dummy = ROOT::Math::PtEtaPhiMVector(particle2.Pt(), particle2.Eta(), particle2.Phi(), lambda2Mass); + auto pairDummy = particle1Dummy + particle2Dummy; + ROOT::Math::Boost boostPairToCM{pairDummy.BoostToCM()}; - template - static inline bool hasSharedDaughtersMC(const A& a, const B& b) - { - return (a.protonIndexmc() == b.protonIndexmc()) || - (a.pionIndexmc() == b.pionIndexmc()) || - (a.protonIndexmc() == b.pionIndexmc()) || - (a.pionIndexmc() == b.protonIndexmc()); - } + // Step1: Boost both Lambdas to pair rest frame + auto lambda1CM = boostPairToCM(particle1Dummy); + auto lambda2CM = boostPairToCM(particle2Dummy); - ROOT::Math::PtEtaPhiMVector lambda0, proton0; - ROOT::Math::PtEtaPhiMVector lambda, proton; - ROOT::Math::PtEtaPhiMVector lambda2, proton2; + // Step2: Boost each Lambda to its own rest frame + ROOT::Math::Boost boostLambda1ToCM{lambda1CM.BoostToCM()}; + ROOT::Math::Boost boostLambda2ToCM{lambda2CM.BoostToCM()}; - Filter centralityFilter = (nabs(aod::lambdaevent::cent) < centMax && nabs(aod::lambdaevent::cent) > centMin); + // Also boost daughter protons to pair CM + auto proton1pairCM = boostPairToCM(daughpart1); + auto proton2pairCM = boostPairToCM(daughpart2); - using EventCandidates = soa::Filtered; - using AllTrackCandidates = aod::LambdaPairs; + // Then into each Lambda rest frame + auto proton1LambdaRF = boostLambda1ToCM(proton1pairCM); + auto proton2LambdaRF = boostLambda2ToCM(proton2pairCM); - void processData(EventCandidates::iterator const& collision, AllTrackCandidates const& V0s) - { - auto centrality = collision.cent(); - for (const auto& v0 : V0s) { - if (!selectionV0(v0)) { - continue; - } - if (fillBasicQAHistos) { - histos.fill(HIST("hPtRadiusV0"), v0.lambdaPt(), v0.v0Radius()); - } - if (fillBasicQAHistos) { - histos.fill(HIST("ptCent"), v0.lambdaPt(), centrality); - } - if (fillBasicQAHistos) { - histos.fill(HIST("etaCent"), v0.lambdaEta(), centrality); - } - proton = ROOT::Math::PtEtaPhiMVector(v0.protonPt(), v0.protonEta(), v0.protonPhi(), o2::constants::physics::MassProton); - lambda = ROOT::Math::PtEtaPhiMVector(v0.lambdaPt(), v0.lambdaEta(), v0.lambdaPhi(), v0.lambdaMass()); - const double phi = RecoDecay::constrainAngle(v0.lambdaPhi(), 0.0F, harmonic); - const double eta = v0.lambdaEta(); + // STAR-style alternative + ROOT::Math::Boost boostL1_LabToRF{particle1Dummy.BoostToCM()}; + ROOT::Math::Boost boostL2_LabToRF{particle2Dummy.BoostToCM()}; - if (v0.v0Status() == 0) { - histos.fill(HIST("hEtaPhiLambdaRaw"), phi, eta, getNUAWeight(0, v0.lambdaPhi(), v0.lambdaEta())); - } else { - histos.fill(HIST("hEtaPhiAntiLambdaRaw"), phi, eta, getNUAWeight(1, v0.lambdaPhi(), v0.lambdaEta())); - } - for (const auto& v02 : V0s) { - if (v02.index() <= v0.index()) { - continue; - } - if (!selectionV0(v02)) { - continue; - } - if (hasSharedDaughters(v0, v02)) - continue; - proton2 = ROOT::Math::PtEtaPhiMVector(v02.protonPt(), v02.protonEta(), v02.protonPhi(), o2::constants::physics::MassProton); - lambda2 = ROOT::Math::PtEtaPhiMVector(v02.lambdaPt(), v02.lambdaEta(), v02.lambdaPhi(), v02.lambdaMass()); - if ((v0.v0Status() == 0 && v02.v0Status() == 1) || (v0.v0Status() == 1 && v02.v0Status() == 0)) - if (fillBasicQAHistos) - histos.fill(HIST("deltaPhiSame"), RecoDecay::constrainAngle(v0.lambdaPhi() - v02.lambdaPhi(), -TMath::Pi(), harmonicDphi)); - // const int ptype = pairTypeCode(v0.v0Status(), v02.v0Status()); - if (v0.v0Status() == 0 && v02.v0Status() == 0) { - fillHistograms(0, 0, lambda, lambda2, proton, proton2, 0, 1.0); - } - if (v0.v0Status() == 0 && v02.v0Status() == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 0, 1.0); - } - if (v0.v0Status() == 1 && v02.v0Status() == 0) { - fillHistograms(0, 1, lambda2, lambda, proton2, proton, 0, 1.0); - } - if (v0.v0Status() == 1 && v02.v0Status() == 1) { - fillHistograms(1, 1, lambda, lambda2, proton, proton2, 0, 1.0); - } - } - } - } - PROCESS_SWITCH(lambdaspincorrderived, processData, "Process data", true); + auto p1_LRF = boostL1_LabToRF(daughpart1); + auto p2_LRF = boostL2_LabToRF(daughpart2); - template - void fillReplacementControlMap(int tag1, int tag2, int leg, bool isTarget, LV const& particle, float weight) - { - if (!fillReplacementQAHistos) { - return; - } - const double pt = particle.Pt(); - const double phi = RecoDecay::constrainAngle(particle.Phi(), 0.0F, harmonic); + TVector3 u1 = TVector3(p1_LRF.Px(), p1_LRF.Py(), p1_LRF.Pz()).Unit(); + TVector3 u2 = TVector3(p2_LRF.Px(), p2_LRF.Py(), p2_LRF.Pz()).Unit(); - double etaOrY = particle.Eta(); - if (userapidity) { - etaOrY = particle.Rapidity(); - } + TVector3 k1(proton1LambdaRF.Px(), proton1LambdaRF.Py(), proton1LambdaRF.Pz()); + k1 = k1.Unit(); + TVector3 k2(proton2LambdaRF.Px(), proton2LambdaRF.Py(), proton2LambdaRF.Pz()); + k2 = k2.Unit(); - if (leg == 1 && isTarget) { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("TGT_LL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("TGT_LAL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("TGT_ALL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("TGT_ALAL_leg1"), phi, etaOrY, pt, weight); - return; - } + double cosDeltaTheta_STAR_naive = u1.Dot(u2); + if (cosDeltaTheta_STAR_naive > 1.0) + cosDeltaTheta_STAR_naive = 111.0; + if (cosDeltaTheta_STAR_naive < -1.0) + cosDeltaTheta_STAR_naive = -111.0; - if (leg == 1 && !isTarget) { + double cosDeltaTheta_hel = k1.Dot(k2); + if (cosDeltaTheta_hel > 1.0) + cosDeltaTheta_hel = 111.0; + if (cosDeltaTheta_hel < -1.0) + cosDeltaTheta_hel = -111.0; - // Raw REP map: used to produce CCDB weight. - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("REP_LL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("REP_LAL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("REP_ALL_leg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("REP_ALAL_leg1"), phi, etaOrY, pt, weight); - // Correct weighted REP QA is filled in fillHistograms(), after the exact final weight is computed. - return; - } + double cosThetaDiff = (cosDef == 0) ? cosDeltaTheta_STAR_naive : cosDeltaTheta_hel; - if (leg == 2 && isTarget) { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("TGT_LL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("TGT_LAL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("TGT_ALL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("TGT_ALAL_leg2"), phi, etaOrY, pt, weight); - return; - } + double pt1 = particle1.Pt(); + double dphi1 = RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic); + double deta1 = particle1.Eta(); - if (leg == 2 && !isTarget) { + double pt2 = particle2.Pt(); + double dphi2 = RecoDecay::constrainAngle(particle2.Phi(), 0.0F, harmonic); + double deta2 = particle2.Eta(); - // Raw REP map: used to produce CCDB weight. - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("REP_LL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("REP_LAL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("REP_ALL_leg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("REP_ALAL_leg2"), phi, etaOrY, pt, weight); - // Correct weighted REP QA is filled in fillHistograms(), after the exact final weight is computed. - return; - } - } - template - void fillFixedLegControlMap(int tag1, int tag2, - int repLeg, - bool isTarget, - LV const& fixedParticle, - float weight) - { - if (!fillFixedLegQAHistos) { - return; - } - const double pt = fixedParticle.Pt(); - const double phi = RecoDecay::constrainAngle(fixedParticle.Phi(), 0.0F, harmonic); + double nuaWeight1 = getNUAWeight(tag1, particle1.Phi(), particle1.Eta()); + double nuaWeight2 = getNUAWeight(tag2, particle2.Phi(), particle2.Eta()); + const double pairNUAWeight = nuaWeight1 * nuaWeight2; - double etaOrY = fixedParticle.Eta(); - if (userapidity) { - etaOrY = fixedParticle.Rapidity(); - } + double dphi_pair = RecoDecay::constrainAngle(dphi1 - dphi2, -TMath::Pi(), harmonicDphi); + double deltaRap = std::abs(particle1.Rapidity() - particle2.Rapidity()); + double deltaR = TMath::Sqrt(deltaRap * deltaRap + dphi_pair * dphi_pair); - if (repLeg == 1) { - // leg1 is replaced, fixed leg is original leg2 - if (isTarget) { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("TGT_FIX_LL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("TGT_FIX_LAL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("TGT_FIX_ALL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("TGT_FIX_ALAL_forRepLeg1"), phi, etaOrY, pt, weight); - } else { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("SUC_FIX_LL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("SUC_FIX_LAL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("SUC_FIX_ALL_forRepLeg1"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("SUC_FIX_ALAL_forRepLeg1"), phi, etaOrY, pt, weight); - } - return; - } + // only for weight lookup; must match fillReplacementControlMap() + double yOrEta1_forWeight = deta1; + double yOrEta2_forWeight = deta2; - if (repLeg == 2) { - // leg2 is replaced, fixed leg is original leg1 - if (isTarget) { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("TGT_FIX_LL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("TGT_FIX_LAL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("TGT_FIX_ALL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("TGT_FIX_ALAL_forRepLeg2"), phi, etaOrY, pt, weight); - } else { - if (tag1 == 0 && tag2 == 0) - histos.fill(HIST("SUC_FIX_LL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 0 && tag2 == 1) - histos.fill(HIST("SUC_FIX_LAL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 0) - histos.fill(HIST("SUC_FIX_ALL_forRepLeg2"), phi, etaOrY, pt, weight); - else if (tag1 == 1 && tag2 == 1) - histos.fill(HIST("SUC_FIX_ALAL_forRepLeg2"), phi, etaOrY, pt, weight); - } - return; + if (userapidity) { + yOrEta1_forWeight = particle1.Rapidity(); + yOrEta2_forWeight = particle2.Rapidity(); } - } - // Processing Event Mixing - SliceCache cache; - using BinningType = ColumnBinningPolicy; - BinningType colBinning{{CfgVtxBins, CfgMultBins}, true}; - Preslice tracksPerCollisionV0 = aod::lambdapair::lambdaeventId; - void processMEV3(EventCandidates const& collisions, AllTrackCandidates const& V0s) - { - auto nBins = colBinning.getAllBinsCount(); - std::vector>> eventPools(nBins); + // `replacedLeg` is the position of the replaced candidate in the ordered pair + // passed to fillHistograms: 1 -> particle1, 2 -> particle2. + // `weightMapLeg` is the physical replacement branch used to select the CCDB map: + // 1 -> REP_*_leg1, 2 -> REP_*_leg2. This distinction is needed for unlike-sign + // pairs where the pair is reordered to Lambda-AntiLambda before filling. + const int replacedPos = replacedLeg; + const int ccdbMapLeg = (weightMapLeg > 0) ? weightMapLeg : replacedLeg; - for (auto& collision1 : collisions) { - const int bin = colBinning.getBin(std::make_tuple(collision1.posz(), collision1.cent())); - if (bin < 0) { - continue; - } + double epsWeightReplaced = 1.0; + double epsWeightFixed = 1.0; - auto poolA = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); + if (useweight && datatype == 1) { + const int wcat = getWeightCategory(tag1, tag2); - // if pool empty, push and continue - if (eventPools[bin].empty()) { - eventPools[bin].emplace_back(collision1.index(), std::move(poolA)); - if ((int)eventPools[bin].size() > nEvtMixing) { - eventPools[bin].pop_front(); + auto getRepEps = [&](int mapLeg, double phi, double yOrEta, double pt) -> double { + TH3D* h = nullptr; + if (mapLeg == 1) { + if (wcat == 0) + h = hweight1; + else if (wcat == 1) + h = hweight2; + else if (wcat == 2) + h = hweight4; + } else if (mapLeg == 2) { + if (wcat == 0) + h = hweight12; + else if (wcat == 1) + h = hweight22; + else if (wcat == 2) + h = hweight42; } - continue; - } + if (!h) { + return 1.0; + } + return h->GetBinContent(h->FindBin(phi, yOrEta, pt)); + }; - for (auto& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { - if (!selectionV0(t1) || !selectionV0(t2)) { - continue; + auto getFixedEps = [&](int mapLeg, double phi, double yOrEta, double pt) -> double { + TH3D* h = nullptr; + if (mapLeg == 1) { + if (wcat == 0) + h = gFixedLLRep1; + else if (wcat == 1) + h = gFixedULRep1; + else if (wcat == 2) + h = gFixedALALRep1; + } else if (mapLeg == 2) { + if (wcat == 0) + h = gFixedLLRep2; + else if (wcat == 1) + h = gFixedULRep2; + else if (wcat == 2) + h = gFixedALALRep2; } - if (t2.index() <= t1.index()) { - continue; + if (!h) { + return 1.0; } - if (hasSharedDaughters(t1, t2)) - continue; - const bool doMixLeg1 = (cfgMixLegMode.value == 0 || cfgMixLegMode.value == 2); - const bool doMixLeg2 = (cfgMixLegMode.value == 1 || cfgMixLegMode.value == 2); + return h->GetBinContent(h->FindBin(phi, yOrEta, pt)); + }; - struct PV { - AllTrackCandidates* pool; - int nRepl1 = 0; - int nRepl2 = 0; - }; + const double phiRep = (replacedPos == 2) ? dphi2 : dphi1; + const double yRep = (replacedPos == 2) ? yOrEta2_forWeight : yOrEta1_forWeight; + const double ptRep = (replacedPos == 2) ? pt2 : pt1; + epsWeightReplaced = getRepEps(ccdbMapLeg, phiRep, yRep, ptRep); - std::vector usable; - int totalRepl = 0; - int totalRepl1 = 0; - int totalRepl2 = 0; + if (cfgCcdbParam.useFixedWeight) { + const int fixedPos = (replacedPos == 2) ? 1 : 2; + const double phiFix = (fixedPos == 2) ? dphi2 : dphi1; + const double yFix = (fixedPos == 2) ? yOrEta2_forWeight : yOrEta1_forWeight; + const double ptFix = (fixedPos == 2) ? pt2 : pt1; + epsWeightFixed = getFixedEps(ccdbMapLeg, phiFix, yFix, ptFix); + } + } - int mixes = 0; - for (auto it = eventPools[bin].rbegin(); it != eventPools[bin].rend() && mixes < nEvtMixing; ++it, ++mixes) { - const int collision2idx = it->first; - auto& poolB = it->second; + if (datatype == 0) { + const double weight = pairNUAWeight; + if (tag1 == 0 && tag2 == 0) { + histos.fill(HIST("hSparseLambdaLambdaSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 0 && tag2 == 1) { + histos.fill(HIST("hSparseLambdaAntiLambdaSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 1 && tag2 == 0) { + histos.fill(HIST("hSparseAntiLambdaLambdaSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 1 && tag2 == 1) { + histos.fill(HIST("hSparseAntiLambdaAntiLambdaSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } - if (collision2idx == collision1.index()) { - continue; - } + } else if (datatype == 1) { + double weight = mixpairweight; - int nRepl1 = 0; - int nRepl2 = 0; + if (useweight) { + const double epsWeightTotal = epsWeightReplaced * epsWeightFixed; - for (auto& tX : poolB) { - if (!selectionV0(tX)) { - continue; - } + if (!std::isfinite(epsWeightTotal) || epsWeightTotal <= 0.0) { + return; + } + weight = mixpairweight / epsWeightTotal; + } - if (doMixLeg1) { - // Single-track replacement: replace a candidate only by the same species. - if (tX.v0Status() == t1.v0Status() && checkKinematics(t1, tX)) { - ++nRepl1; - } - } + // This is the pure mixing-correction weight. + // Do not include NUA here, because TGT/REP/FIX QA maps were filled without NUA. + const double weightMixingQA = weight; - if (doMixLeg2) { - // Single-track replacement: replace a candidate only by the same species. - if (tX.v0Status() == t2.v0Status() && checkKinematics(t2, tX)) { - ++nRepl2; - } - } - } + if (useweight) { + fillFinalWeightedMixingQA(tag1, tag2, ccdbMapLeg, replacedPos, particle1, particle2, weightMixingQA); + } - if (nRepl1 > 0 || nRepl2 > 0) { - usable.push_back(PV{&poolB, nRepl1, nRepl2}); - totalRepl += nRepl1 + nRepl2; - totalRepl1 += nRepl1; - totalRepl2 += nRepl2; - } - } + weight *= pairNUAWeight; + if (!std::isfinite(weight) || weight <= 0.0) { + return; + } - if (totalRepl <= 0) { - continue; - } + if (tag1 == 0 && tag2 == 0) { + histos.fill(HIST("hSparseLambdaLambdaMixedSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 0 && tag2 == 1) { + histos.fill(HIST("hSparseLambdaAntiLambdaMixedSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 1 && tag2 == 0) { + histos.fill(HIST("hSparseAntiLambdaLambdaMixedSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } else if (tag1 == 1 && tag2 == 1) { + histos.fill(HIST("hSparseAntiLambdaAntiLambdaMixedSys"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, SysId, weight); + } + } + } - const float wBase = 1.0f / static_cast(totalRepl); + static inline int pairTypeCode(int tag1, int tag2) + { + if (tag1 == 0 && tag2 == 0) { + return 0; // LL + } else if (tag1 == 0 && tag2 == 1) { + return 1; // LAL + } else if (tag1 == 1 && tag2 == 0) { + return 2; // ALL + } else { + return 3; // ALAL + } + } - // Single-track replacement target must use the same branch normalization - // as the actually used replacement candidates. Per SE pair: - // sum REP_leg1 weights = totalRepl1 / totalRepl - // sum REP_leg2 weights = totalRepl2 / totalRepl - // Therefore TGT leg1/leg2 are filled with the same weights. - if (doMixLeg1 && totalRepl1 > 0) { - fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 1, true, - ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), - static_cast(totalRepl1) * wBase); - fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 1, true, - ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), - static_cast(totalRepl1) * wBase); + template + static inline bool hasSharedDaughters(const A& a, const B& b) + { + return (a.protonIndex() == b.protonIndex()) || + (a.pionIndex() == b.pionIndex()) || + (a.protonIndex() == b.pionIndex()) || + (a.pionIndex() == b.protonIndex()); + } + + template + static inline bool hasSharedDaughtersMC(const A& a, const B& b) + { + return (a.protonIndexmc() == b.protonIndexmc()) || + (a.pionIndexmc() == b.pionIndexmc()) || + (a.protonIndexmc() == b.pionIndexmc()) || + (a.pionIndexmc() == b.protonIndexmc()); + } + + ROOT::Math::PtEtaPhiMVector lambda0, proton0; + ROOT::Math::PtEtaPhiMVector lambda, proton; + ROOT::Math::PtEtaPhiMVector lambda2, proton2; + + Filter centralityFilter = (nabs(aod::lambdaevent::cent) < centMax && nabs(aod::lambdaevent::cent) > centMin); + + using EventCandidates = soa::Filtered; + using AllTrackCandidates = aod::LambdaPairs; + + void processData(EventCandidates::iterator const& collision, AllTrackCandidates const& V0s) + { + auto centrality = collision.cent(); + for (const auto& v0 : V0s) { + if (!selectionV0(v0)) { + continue; + } + if (fillBasicQAHistos) { + histos.fill(HIST("hPtRadiusV0"), v0.lambdaPt(), v0.v0Radius()); + } + if (fillBasicQAHistos) { + histos.fill(HIST("ptCent"), v0.lambdaPt(), centrality); + } + if (fillBasicQAHistos) { + histos.fill(HIST("etaCent"), v0.lambdaEta(), centrality); + } + proton = ROOT::Math::PtEtaPhiMVector(v0.protonPt(), v0.protonEta(), v0.protonPhi(), o2::constants::physics::MassProton); + lambda = ROOT::Math::PtEtaPhiMVector(v0.lambdaPt(), v0.lambdaEta(), v0.lambdaPhi(), v0.lambdaMass()); + const double phi = RecoDecay::constrainAngle(v0.lambdaPhi(), 0.0F, harmonic); + const double eta = v0.lambdaEta(); + + if (v0.v0Status() == 0) { + histos.fill(HIST("hEtaPhiLambdaRaw"), phi, eta, getNUAWeight(0, v0.lambdaPhi(), v0.lambdaEta())); + } else { + histos.fill(HIST("hEtaPhiAntiLambdaRaw"), phi, eta, getNUAWeight(1, v0.lambdaPhi(), v0.lambdaEta())); + } + for (const auto& v02 : V0s) { + if (v02.index() <= v0.index()) { + continue; } - if (doMixLeg2 && totalRepl2 > 0) { - fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 2, true, - ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), - static_cast(totalRepl2) * wBase); - fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 2, true, - ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), - static_cast(totalRepl2) * wBase); + if (!selectionV0(v02)) { + continue; + } + if (hasSharedDaughters(v0, v02)) + continue; + proton2 = ROOT::Math::PtEtaPhiMVector(v02.protonPt(), v02.protonEta(), v02.protonPhi(), o2::constants::physics::MassProton); + lambda2 = ROOT::Math::PtEtaPhiMVector(v02.lambdaPt(), v02.lambdaEta(), v02.lambdaPhi(), v02.lambdaMass()); + if ((v0.v0Status() == 0 && v02.v0Status() == 1) || (v0.v0Status() == 1 && v02.v0Status() == 0)) + if (fillBasicQAHistos) + histos.fill(HIST("deltaPhiSame"), RecoDecay::constrainAngle(v0.lambdaPhi() - v02.lambdaPhi(), -TMath::Pi(), harmonicDphi)); + // const int ptype = pairTypeCode(v0.v0Status(), v02.v0Status()); + if (v0.v0Status() == 0 && v02.v0Status() == 0) { + fillHistograms(0, 0, lambda, lambda2, proton, proton2, 0, 1.0); + } + if (v0.v0Status() == 0 && v02.v0Status() == 1) { + fillHistograms(0, 1, lambda, lambda2, proton, proton2, 0, 1.0); + } + if (v0.v0Status() == 1 && v02.v0Status() == 0) { + fillHistograms(0, 1, lambda2, lambda, proton2, proton, 0, 1.0); + } + if (v0.v0Status() == 1 && v02.v0Status() == 1) { + fillHistograms(1, 1, lambda, lambda2, proton, proton2, 0, 1.0); } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processData, "Process data", true); - for (auto& pv : usable) { - auto& poolB = *pv.pool; + void processDataSys(EventCandidates::iterator const&, AllTrackCandidates const& V0s) + { + for (const auto& v0 : V0s) { - for (auto& tX : poolB) { - if (!selectionV0(tX)) { - continue; - } + std::vector activeSys; + activeSys.reserve((size_t)nSysTotal); - // -------- leg-1 replacement: (tX, t2) - if (doMixLeg1) { - if (tX.v0Status() == t1.v0Status() && checkKinematics(t1, tX)) { - fillReplacementControlMap(tX.v0Status(), t2.v0Status(), 1, false, - ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), - wBase); - 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()); + for (int sysId = 0; sysId < nSysTotal; ++sysId) { + const auto& sc = sysCuts[sysId]; + if (selectionV0Sys(v0, sc)) { + activeSys.push_back(sysId); + } + } - const float dPhi = RecoDecay::constrainAngle( - RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) - - RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic), - -TMath::Pi(), harmonicDphi); + if (activeSys.empty()) + continue; - if (fillBasicQAHistos) - histos.fill(HIST("deltaPhiMix"), dPhi, wBase); - fillHistograms(tX.v0Status(), t2.v0Status(), - lambda, lambda2, proton, proton2, - 1, wBase, 1); - } - } + proton = ROOT::Math::PtEtaPhiMVector(v0.protonPt(), v0.protonEta(), v0.protonPhi(), + o2::constants::physics::MassProton); + lambda = ROOT::Math::PtEtaPhiMVector(v0.lambdaPt(), v0.lambdaEta(), v0.lambdaPhi(), + v0.lambdaMass()); - // -------- leg-2 replacement: (t1, tX) - if (doMixLeg2) { - if (tX.v0Status() == t2.v0Status() && checkKinematics(t2, tX)) { - fillReplacementControlMap(t1.v0Status(), tX.v0Status(), 2, false, - ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), - wBase); - 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(tX.protonPt(), tX.protonEta(), tX.protonPhi(), - o2::constants::physics::MassProton); - auto lambda2 = ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), - tX.lambdaMass()); + for (const auto& v02 : V0s) { + if (v02.index() <= v0.index()) + continue; + if (hasSharedDaughters(v0, v02)) + continue; - const float dPhi = RecoDecay::constrainAngle( - RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) - - RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic), - -TMath::Pi(), harmonicDphi); + std::vector activePair; + activePair.reserve(activeSys.size()); - if (fillBasicQAHistos) - histos.fill(HIST("deltaPhiMix"), dPhi, wBase); - fillHistograms(t1.v0Status(), tX.v0Status(), - lambda, lambda2, proton, proton2, - 1, wBase, 2); - } - } + for (const auto& sysId : activeSys) { + const auto& sc = sysCuts[sysId]; + if (selectionV0Sys(v02, sc)) { + activePair.push_back(sysId); } } - } - // push current event into pool - auto sliced = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); - eventPools[bin].emplace_back(collision1.index(), std::move(sliced)); - if ((int)eventPools[bin].size() > nEvtMixing) { - eventPools[bin].pop_front(); + if (activePair.empty()) { + continue; + } + + proton2 = ROOT::Math::PtEtaPhiMVector(v02.protonPt(), v02.protonEta(), v02.protonPhi(), + o2::constants::physics::MassProton); + lambda2 = ROOT::Math::PtEtaPhiMVector(v02.lambdaPt(), v02.lambdaEta(), v02.lambdaPhi(), + v02.lambdaMass()); + + for (const auto& sysId : activePair) { + + if (v0.v0Status() == 0 && v02.v0Status() == 0) { + fillHistogramsSys(0, 0, lambda, lambda2, proton, proton2, 0, 1.0, sysId); + } else if (v0.v0Status() == 0 && v02.v0Status() == 1) { + fillHistogramsSys(0, 1, lambda, lambda2, proton, proton2, 0, 1.0, sysId); + } else if (v0.v0Status() == 1 && v02.v0Status() == 0) { + fillHistogramsSys(0, 1, lambda2, lambda, proton2, proton, 0, 1.0, sysId); + } else if (v0.v0Status() == 1 && v02.v0Status() == 1) { + fillHistogramsSys(1, 1, lambda, lambda2, proton, proton2, 0, 1.0, sysId); + } + } } } } - PROCESS_SWITCH(lambdaspincorrderived, processMEV3, "Process data ME (first-leg, pair-3D maps)", false); - - static constexpr int N_STATUS = 2; // v0Status ∈ {0,1} - struct MatchRef { - int64_t collisionIdx; - int64_t rowIndex; - }; + PROCESS_SWITCH(lambdaspincorrderived, processDataSys, "Process Sys analysis", true); - static inline void limitMatchesToNEvents(std::vector& matches, int nMixEvents) + template + void fillReplacementControlMap(int tag1, int tag2, int leg, bool isTarget, LV const& particle, float weight) { - if (nMixEvents <= 0 || matches.empty()) { + + if (!fillReplacementQAHistos) { return; } - std::vector kept; - kept.reserve(matches.size()); + const double pt = particle.Pt(); + const double phi = RecoDecay::constrainAngle(particle.Phi(), 0.0F, harmonic); - std::unordered_set usedEvents; - usedEvents.reserve(nMixEvents * 2); + double etaOrY = particle.Eta(); + if (userapidity) { + etaOrY = particle.Rapidity(); + } - for (const auto& m : matches) { - if (usedEvents.count(m.collisionIdx) || (int)usedEvents.size() < nMixEvents) { - kept.push_back(m); - usedEvents.insert(m.collisionIdx); - } + if (leg == 1 && isTarget) { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("TGT_LL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("TGT_LAL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("TGT_ALL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("TGT_ALAL_leg1"), phi, etaOrY, pt, weight); + return; } - matches.swap(kept); - } - struct MixBinnerR { - float ptMin, ptMax, ptStep; - float etaMin, etaMax, etaStep; - float phiMin, phiMax, phiStep; + if (leg == 1 && !isTarget) { - float mMin, mMax, mStep; - int nM_; + // Raw REP map: used to produce CCDB weight. + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("REP_LL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("REP_LAL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("REP_ALL_leg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("REP_ALAL_leg1"), phi, etaOrY, pt, weight); + // Correct weighted REP QA is filled in fillHistograms(), after the exact final weight is computed. + return; + } - std::vector rEdges; - int nR_; + if (leg == 2 && isTarget) { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("TGT_LL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("TGT_LAL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("TGT_ALL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("TGT_ALAL_leg2"), phi, etaOrY, pt, weight); + return; + } - int nPt_, nEta_, nPhi_; + if (leg == 2 && !isTarget) { - MixBinnerR(float ptMin_, float ptMax_, float ptStep_, - float etaAbsMax, float etaStep_, - float phiStep_, - float mMin_, float mMax_, int nMassBins_, - std::vector rEdges_) + // Raw REP map: used to produce CCDB weight. + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("REP_LL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("REP_LAL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("REP_ALL_leg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("REP_ALAL_leg2"), phi, etaOrY, pt, weight); + // Correct weighted REP QA is filled in fillHistograms(), after the exact final weight is computed. + return; + } + } + template + void fillFixedLegControlMap(int tag1, int tag2, + int repLeg, + bool isTarget, + LV const& fixedParticle, + float weight) + { + if (!fillFixedLegQAHistos) { + return; + } + const double pt = fixedParticle.Pt(); + const double phi = RecoDecay::constrainAngle(fixedParticle.Phi(), 0.0F, harmonic); + + double etaOrY = fixedParticle.Eta(); + if (userapidity) { + etaOrY = fixedParticle.Rapidity(); + } + + if (repLeg == 1) { + // leg1 is replaced, fixed leg is original leg2 + if (isTarget) { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("TGT_FIX_LL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("TGT_FIX_LAL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("TGT_FIX_ALL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("TGT_FIX_ALAL_forRepLeg1"), phi, etaOrY, pt, weight); + } else { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("SUC_FIX_LL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("SUC_FIX_LAL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("SUC_FIX_ALL_forRepLeg1"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("SUC_FIX_ALAL_forRepLeg1"), phi, etaOrY, pt, weight); + } + return; + } + + if (repLeg == 2) { + // leg2 is replaced, fixed leg is original leg1 + if (isTarget) { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("TGT_FIX_LL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("TGT_FIX_LAL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("TGT_FIX_ALL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("TGT_FIX_ALAL_forRepLeg2"), phi, etaOrY, pt, weight); + } else { + if (tag1 == 0 && tag2 == 0) + histos.fill(HIST("SUC_FIX_LL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 0 && tag2 == 1) + histos.fill(HIST("SUC_FIX_LAL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 0) + histos.fill(HIST("SUC_FIX_ALL_forRepLeg2"), phi, etaOrY, pt, weight); + else if (tag1 == 1 && tag2 == 1) + histos.fill(HIST("SUC_FIX_ALAL_forRepLeg2"), phi, etaOrY, pt, weight); + } + return; + } + } + // Processing Event Mixing + SliceCache cache; + using BinningType = ColumnBinningPolicy; + BinningType colBinning{{CfgVtxBins, CfgMultBins}, true}; + Preslice tracksPerCollisionV0 = aod::lambdapair::lambdaeventId; + + void processMEV3(EventCandidates const& collisions, AllTrackCandidates const& V0s) + { + auto nBins = colBinning.getAllBinsCount(); + std::vector>> eventPools(nBins); + + for (const auto& collision1 : collisions) { + const int bin = colBinning.getBin(std::make_tuple(collision1.posz(), collision1.cent())); + if (bin < 0) { + continue; + } + + auto poolA = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); + + // if pool empty, push and continue + if (eventPools[bin].empty()) { + eventPools[bin].emplace_back(collision1.index(), std::move(poolA)); + if ((int)eventPools[bin].size() > nEvtMixing) { + eventPools[bin].pop_front(); + } + continue; + } + + for (const auto& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { + if (!selectionV0(t1) || !selectionV0(t2)) { + continue; + } + if (t2.index() <= t1.index()) { + continue; + } + if (hasSharedDaughters(t1, t2)) + continue; + const bool doMixLeg1 = (cfgMixLegMode.value == 0 || cfgMixLegMode.value == 2); + const bool doMixLeg2 = (cfgMixLegMode.value == 1 || cfgMixLegMode.value == 2); + + struct PV { + AllTrackCandidates* pool; + int nRepl1 = 0; + int nRepl2 = 0; + }; + + std::vector usable; + int totalRepl = 0; + int totalRepl1 = 0; + int totalRepl2 = 0; + + int mixes = 0; + for (auto it = eventPools[bin].rbegin(); it != eventPools[bin].rend() && mixes < nEvtMixing; ++it, ++mixes) { + const int collision2idx = it->first; + auto& poolB = it->second; + + if (collision2idx == collision1.index()) { + continue; + } + + int nRepl1 = 0; + int nRepl2 = 0; + + for (const auto& tX : poolB) { + if (!selectionV0(tX)) { + continue; + } + + if (doMixLeg1) { + // Single-track replacement: replace a candidate only by the same species. + if (tX.v0Status() == t1.v0Status() && checkKinematics(t1, tX)) { + ++nRepl1; + } + } + + if (doMixLeg2) { + // Single-track replacement: replace a candidate only by the same species. + if (tX.v0Status() == t2.v0Status() && checkKinematics(t2, tX)) { + ++nRepl2; + } + } + } + + if (nRepl1 > 0 || nRepl2 > 0) { + usable.push_back(PV{&poolB, nRepl1, nRepl2}); + totalRepl += nRepl1 + nRepl2; + totalRepl1 += nRepl1; + totalRepl2 += nRepl2; + } + } + + if (totalRepl <= 0) { + continue; + } + + const float wBase = 1.0f / static_cast(totalRepl); + + // Single-track replacement target must use the same branch normalization + // as the actually used replacement candidates. Per SE pair: + // sum REP_leg1 weights = totalRepl1 / totalRepl + // sum REP_leg2 weights = totalRepl2 / totalRepl + // Therefore TGT leg1/leg2 are filled with the same weights. + if (doMixLeg1 && totalRepl1 > 0) { + fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 1, true, + ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), + static_cast(totalRepl1) * wBase); + fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 1, true, + ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), + static_cast(totalRepl1) * wBase); + } + if (doMixLeg2 && totalRepl2 > 0) { + fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 2, true, + ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), + static_cast(totalRepl2) * wBase); + fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 2, true, + ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), + static_cast(totalRepl2) * wBase); + } + + for (const auto& pv : usable) { + auto& poolB = *pv.pool; + + for (const auto& tX : poolB) { + if (!selectionV0(tX)) { + continue; + } + + // -------- leg-1 replacement: (tX, t2) + if (doMixLeg1) { + if (tX.v0Status() == t1.v0Status() && checkKinematics(t1, tX)) { + fillReplacementControlMap(tX.v0Status(), t2.v0Status(), 1, false, + ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), + wBase); + 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 = RecoDecay::constrainAngle( + RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) - + RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic), + -TMath::Pi(), harmonicDphi); + + if (fillBasicQAHistos) + histos.fill(HIST("deltaPhiMix"), dPhi, wBase); + fillHistograms(tX.v0Status(), t2.v0Status(), + lambda, lambda2, proton, proton2, + 1, wBase, 1); + } + } + + // -------- leg-2 replacement: (t1, tX) + if (doMixLeg2) { + if (tX.v0Status() == t2.v0Status() && checkKinematics(t2, tX)) { + fillReplacementControlMap(t1.v0Status(), tX.v0Status(), 2, false, + ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), + wBase); + 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(tX.protonPt(), tX.protonEta(), tX.protonPhi(), + o2::constants::physics::MassProton); + auto lambda2 = ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), + tX.lambdaMass()); + + const float dPhi = RecoDecay::constrainAngle( + RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) - + RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic), + -TMath::Pi(), harmonicDphi); + + if (fillBasicQAHistos) + histos.fill(HIST("deltaPhiMix"), dPhi, wBase); + fillHistograms(t1.v0Status(), tX.v0Status(), + lambda, lambda2, proton, proton2, + 1, wBase, 2); + } + } + } + } + } + + // push current event into pool + auto sliced = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); + eventPools[bin].emplace_back(collision1.index(), std::move(sliced)); + if ((int)eventPools[bin].size() > nEvtMixing) { + eventPools[bin].pop_front(); + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMEV3, "Process data ME (first-leg, pair-3D maps)", false); + + static constexpr int N_STATUS = 2; // v0Status ∈ {0,1} + struct MatchRef { + int64_t collisionIdx; + int64_t rowIndex; + }; + + static inline void limitMatchesToNEvents(std::vector& matches, int nMixEvents) + { + if (nMixEvents <= 0 || matches.empty()) { + return; + } + std::vector kept; + kept.reserve(matches.size()); + + std::unordered_set usedEvents; + usedEvents.reserve(nMixEvents * 2); + + for (const auto& m : matches) { + if (usedEvents.count(m.collisionIdx) || (int)usedEvents.size() < nMixEvents) { + kept.push_back(m); + usedEvents.insert(m.collisionIdx); + } + } + matches.swap(kept); + } + + struct MixBinnerR { + float ptMin, ptMax, ptStep; + float etaMin, etaMax, etaStep; + float phiMin, phiMax, phiStep; + + float mMin, mMax, mStep; + int nM_; + + std::vector rEdges; + int nR_; + + int nPt_, nEta_, nPhi_; + + MixBinnerR(float ptMin_, float ptMax_, float ptStep_, + float etaAbsMax, float etaStep_, + float phiStep_, + float mMin_, float mMax_, int nMassBins_, + std::vector rEdges_) : ptMin(ptMin_), ptMax(ptMax_), ptStep(ptStep_), @@ -1561,576 +2017,1116 @@ struct lambdaspincorrderived { 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_); - if (!(mStep > 0.f)) { - nM_ = 5; - mMin = 1.09f; - mMax = 1.14f; - mStep = (mMax - mMin) / static_cast(nM_); - } + if (!(mMax > mMin)) { + mMin = 1.09f; + mMax = 1.14f; + } + mStep = (mMax - mMin) / static_cast(nM_); + if (!(mStep > 0.f)) { + nM_ = 5; + mMin = 1.09f; + mMax = 1.14f; + mStep = (mMax - mMin) / static_cast(nM_); + } + + if (rEdges.size() < 2) { + rEdges = {3.0, 5.0, 7.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0}; + } + nR_ = static_cast(rEdges.size()) - 1; + + 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))); + } + + inline int nPt() const { return nPt_; } + inline int nEta() const { return nEta_; } + inline int nPhi() const { return nPhi_; } + inline int nM() const { return nM_; } + inline int nR() const { return nR_; } + + inline int binFromValue(float v, float vmin, float step, int nBins) const + { + 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) { + return -1; + } + if (b >= nBins) { + b = nBins - 1; + } + return b; + } + + inline int ptBin(float pt) const { return binFromValue(pt, ptMin, ptStep, nPt_); } + inline int etaBin(float eta) const { return binFromValue(eta, etaMin, etaStep, nEta_); } + inline int phiBin(float phi) const { return binFromValue(phi, phiMin, phiStep, nPhi_); } + inline int radiusBin(float r) const + { + if (!std::isfinite(r) || nR_ <= 0) { + return -1; + } + if (r < rEdges.front() || r >= rEdges.back()) { + return -1; + } + auto it = std::upper_bound(rEdges.begin(), rEdges.end(), static_cast(r)); + return static_cast(it - rEdges.begin()) - 1; + } + }; + + struct BufferCandR { + int64_t collisionIdx; + int64_t rowIndex; + uint8_t v0Status; + uint16_t ptBin, etaBin, phiBin, mBin, rBin; + }; + + struct StoredV6Candidate { + int64_t collisionIdx = -1; + int64_t globalIdx = -1; + int status = -1; + bool isDouble = false; + float cospa = 0.f; + float radius = 0.f; + float dcaPos = 0.f; + float dcaNeg = 0.f; + float dcaDau = 0.f; + float lPt = 0.f; + float lEta = 0.f; + float lPhi = 0.f; + float lMass = 0.f; + float pPt = 0.f; + float pEta = 0.f; + float pPhi = 0.f; + int64_t pIndex = -1; + int64_t piIndex = -1; + + int v0Status() const { return status; } + bool doubleStatus() const { return isDouble; } + float v0Cospa() const { return cospa; } + float v0Radius() const { return radius; } + float dcaPositive() const { return dcaPos; } + float dcaNegative() const { return dcaNeg; } + float dcaBetweenDaughter() const { return dcaDau; } + float lambdaPt() const { return lPt; } + float lambdaEta() const { return lEta; } + float lambdaPhi() const { return lPhi; } + float lambdaMass() const { return lMass; } + float protonPt() const { return pPt; } + float protonEta() const { return pEta; } + float protonPhi() const { return pPhi; } + int64_t protonIndex() const { return pIndex; } + int64_t pionIndex() const { return piIndex; } + int64_t globalIndex() const { return globalIdx; } + }; + + struct PendingV6Branch { + StoredV6Candidate target; + StoredV6Candidate fixed; + int colBin = -1; + int replacedLeg = 1; + int age = 0; + uint64_t seed = 0; + }; + + template + StoredV6Candidate storeV6Candidate(T const& t, int64_t collisionIdx) const + { + return {collisionIdx, static_cast(t.globalIndex()), static_cast(t.v0Status()), static_cast(t.doubleStatus()), + t.v0Cospa(), t.v0Radius(), t.dcaPositive(), t.dcaNegative(), t.dcaBetweenDaughter(), + t.lambdaPt(), t.lambdaEta(), t.lambdaPhi(), t.lambdaMass(), + t.protonPt(), t.protonEta(), t.protonPhi(), + static_cast(t.protonIndex()), static_cast(t.pionIndex())}; + } + + template + void fillV6MixedBranch(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) + { + const auto repProton = ROOT::Math::PtEtaPhiMVector(replacement.protonPt(), replacement.protonEta(), replacement.protonPhi(), o2::constants::physics::MassProton); + const auto repLambda = ROOT::Math::PtEtaPhiMVector(replacement.lambdaPt(), replacement.lambdaEta(), replacement.lambdaPhi(), replacement.lambdaMass()); + const auto fixedProton = ROOT::Math::PtEtaPhiMVector(fixed.protonPt(), fixed.protonEta(), fixed.protonPhi(), o2::constants::physics::MassProton); + const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(fixed.lambdaPt(), fixed.lambdaEta(), fixed.lambdaPhi(), fixed.lambdaMass()); + const int repStatus = replacement.v0Status(); + const int fixedStatus = fixed.v0Status(); + + if (replacedLeg == 1) { + fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); + fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); + if ((repStatus == 0 && fixedStatus == 1) || (repStatus == 1 && fixedStatus == 0)) { + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); + } + } + if (repStatus == 0 && fixedStatus == 1) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } else if (repStatus == 1 && fixedStatus == 0) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); + } else { + fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } + return; + } + + fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); + fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); + } + if (fixedStatus == 0 && repStatus == 1) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } else if (fixedStatus == 1 && repStatus == 0) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); + } else { + fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } + } + + struct StoredV6CandidateMC { + int64_t globalIdx = -1; + int status = -1; + bool isDouble = false; + float cospa = 0.f; + float radius = 0.f; + float dcaPos = 0.f; + float dcaNeg = 0.f; + float dcaDau = 0.f; + float lPt = 0.f; + float lEta = 0.f; + float lPhi = 0.f; + float lMass = 0.f; + float pPt = 0.f; + float pEta = 0.f; + float pPhi = 0.f; + int64_t pIndex = -1; + int64_t piIndex = -1; + + int v0Statusmc() const { return status; } + bool doubleStatusmc() const { return isDouble; } + float v0Cospamc() const { return cospa; } + float v0Radiusmc() const { return radius; } + float dcaPositivemc() const { return dcaPos; } + float dcaNegativemc() const { return dcaNeg; } + float dcaBetweenDaughtermc() const { return dcaDau; } + float lambdaPtmc() const { return lPt; } + float lambdaEtamc() const { return lEta; } + float lambdaPhimc() const { return lPhi; } + float lambdaMassmc() const { return lMass; } + float protonPtmc() const { return pPt; } + float protonEtamc() const { return pEta; } + float protonPhimc() const { return pPhi; } + int64_t protonIndexmc() const { return pIndex; } + int64_t pionIndexmc() const { return piIndex; } + int64_t globalIndex() const { return globalIdx; } + }; + + struct PendingV6BranchMC { + StoredV6CandidateMC target; + StoredV6CandidateMC fixed; + int colBin = -1; + int replacedLeg = 1; + int age = 0; + uint64_t seed = 0; + }; + + struct V6PendingState { + std::deque data; + std::deque mc; + }; + + V6PendingState v6Pending; + + template + StoredV6CandidateMC storeV6CandidateMC(T const& t) const + { + return {static_cast(t.globalIndex()), mcacc::v0Status(t), mcacc::doubleStatus(t), + mcacc::v0CosPA(t), mcacc::v0Radius(t), mcacc::dcaPos(t), mcacc::dcaNeg(t), mcacc::dcaDau(t), + mcacc::lamPt(t), mcacc::lamEta(t), mcacc::lamPhi(t), mcacc::lamMass(t), + mcacc::prPt(t), mcacc::prEta(t), mcacc::prPhi(t), + static_cast(mcacc::prIdx(t)), static_cast(mcacc::piIdx(t))}; + } + + template + void fillV6MixedBranchMC(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) + { + const auto repProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(replacement), mcacc::prEta(replacement), mcacc::prPhi(replacement), o2::constants::physics::MassProton); + const auto repLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(replacement), mcacc::lamEta(replacement), mcacc::lamPhi(replacement), mcacc::lamMass(replacement)); + const auto fixedProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(fixed), mcacc::prEta(fixed), mcacc::prPhi(fixed), o2::constants::physics::MassProton); + const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(fixed), mcacc::lamEta(fixed), mcacc::lamPhi(fixed), mcacc::lamMass(fixed)); + const int repStatus = mcacc::v0Status(replacement); + const int fixedStatus = mcacc::v0Status(fixed); + + if (replacedLeg == 1) { + fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); + fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); + } + if (repStatus == 0 && fixedStatus == 1) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } else if (repStatus == 1 && fixedStatus == 0) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); + } else { + fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + } + return; + } + + fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); + fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); + if (fillBasicQAHistos) { + histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); + } + if (fixedStatus == 0 && repStatus == 1) { + fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } else if (fixedStatus == 1 && repStatus == 0) { + fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); + } else { + fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); + } + } + + static inline size_t linearKeyR(int colBin, int statBin, + int ptBin, int etaBin, int phiBin, int mBin, int rBin, + int nStatus, int nPt, int nEta, int nPhi, int nM, int nR) + { + return (((((((static_cast(colBin) * nStatus + statBin) * nPt + ptBin) * nEta + etaBin) * nPhi + phiBin) * nM + mBin) * nR + rBin)); + } + + // ------------------------------------- + // 2) MC-only selection + kinematics cuts + // ------------------------------------- + template + bool selectionV0MC(T const& candidate) + { + auto particle = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate), + mcacc::lamEta(candidate), + mcacc::lamPhi(candidate), + mcacc::lamMass(candidate)); + if (std::abs(particle.Rapidity()) > rapidity || std::abs(particle.Eta()) > v0eta) { + return false; + } + if (mcacc::lamMass(candidate) < MassMin || mcacc::lamMass(candidate) > MassMax) { + return false; + } + if (mcacc::v0CosPA(candidate) < v0Configurations.cosPA) { + return false; + } + if (checkDoubleStatus && mcacc::doubleStatus(candidate)) { + return false; + } + if (mcacc::v0Radius(candidate) > v0Configurations.radiusMax) { + return false; + } + if (mcacc::v0Radius(candidate) < v0Configurations.radiusMin) { + return false; + } + if (mcacc::dcaDau(candidate) > v0Configurations.dcaDaughters) { + return false; + } + + if (mcacc::dcaV0ToPVMC(candidate) > v0Configurations.dcaV0ToPV) { + return false; + } + + if (mcacc::v0Status(candidate) == 0 && (std::abs(mcacc::dcaPos(candidate)) < v0Configurations.dcaProton || std::abs(mcacc::dcaNeg(candidate)) < v0Configurations.dcaPion)) { + return false; + } + if (mcacc::v0Status(candidate) == 1 && (std::abs(mcacc::dcaPos(candidate)) < v0Configurations.dcaPion || std::abs(mcacc::dcaNeg(candidate)) < v0Configurations.dcaProton)) { + return false; + } + if (mcacc::lamPt(candidate) < ptMin) { + return false; + } + if (mcacc::lamPt(candidate) > ptMax) { + return false; + } + return true; + } + + template + bool checkKinematicsMC(T1 const& candidate1, T2 const& candidate2) + { + // keep same species/status + if (mcacc::v0Status(candidate1) != mcacc::v0Status(candidate2)) { + return false; + } + + // pT window + if (std::abs(mcacc::lamPt(candidate1) - mcacc::lamPt(candidate2)) > ptMix) { + return false; + } + + // eta or rapidity window (etaMix used as Δη or Δy) + if (!userapidity) { + if (std::abs(mcacc::lamEta(candidate1) - mcacc::lamEta(candidate2)) > etaMix) { + return false; + } + } else { + const auto l1 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate1), mcacc::lamEta(candidate1), + mcacc::lamPhi(candidate1), mcacc::lamMass(candidate1)); + const auto l2 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate2), mcacc::lamEta(candidate2), + mcacc::lamPhi(candidate2), mcacc::lamMass(candidate2)); + if (std::abs(l1.Rapidity() - l2.Rapidity()) > etaMix) { + return false; + } + } + + // delta-phi window (wrapped) + const float dphi = deltaPhiMinusPiToPi((float)mcacc::lamPhi(candidate1), + (float)mcacc::lamPhi(candidate2)); + if (std::abs(dphi) > phiMix) { + return false; + } + + // mass window (optional but consistent with data) + if (std::abs(mcacc::lamMass(candidate1) - mcacc::lamMass(candidate2)) > massMix) { + return false; + } + + return true; + } + + // ----------------------------------------- + // 3) MC filter + aliases (distinct from data) + // ----------------------------------------- + Filter centralityFilterMC = (nabs(aod::lambdaeventmc::centmc) < centMax && nabs(aod::lambdaeventmc::centmc) > centMin); + + using EventCandidatesMC = soa::Filtered; + using AllTrackCandidatesMC = aod::LambdaPairmcs; + + // IMPORTANT: MC preslice uses the MC event index column + Preslice tracksPerCollisionV0mc = aod::lambdapairmc::lambdaeventmcId; + + // ----------------------------------------- + // 4) MC Same-event processing (like processData) + // ----------------------------------------- + void processMC(EventCandidatesMC::iterator const& collision, AllTrackCandidatesMC const& V0sMC) + { + const float centrality = mcacc::cent(collision); + + for (const auto& v0 : V0sMC) { + if (!selectionV0MC(v0)) { + continue; + } + if (fillBasicQAHistos) { + histos.fill(HIST("hPtRadiusV0"), mcacc::lamPt(v0), mcacc::v0Radius(v0)); + } + if (fillBasicQAHistos) { + histos.fill(HIST("ptCent"), mcacc::lamPt(v0), centrality); + } + if (fillBasicQAHistos) { + histos.fill(HIST("etaCent"), mcacc::lamEta(v0), centrality); + } + + proton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(v0), mcacc::prEta(v0), mcacc::prPhi(v0), + o2::constants::physics::MassProton); + lambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(v0), mcacc::lamEta(v0), mcacc::lamPhi(v0), + mcacc::lamMass(v0)); + if (mcacc::v0Status(v0) == 0) { + histos.fill(HIST("hEtaPhiLambdaRaw"), lambda.Phi(), lambda.Eta(), getNUAWeight(0, lambda.Phi(), lambda.Eta())); + } else { + histos.fill(HIST("hEtaPhiAntiLambdaRaw"), lambda.Phi(), lambda.Eta(), getNUAWeight(1, lambda.Phi(), lambda.Eta())); + } + + for (const auto& v02 : V0sMC) { + if (v02.index() <= v0.index()) { + continue; + } + if (!selectionV0MC(v02)) { + continue; + } + if (hasSharedDaughtersMC(v0, v02)) + continue; + proton2 = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(v02), mcacc::prEta(v02), mcacc::prPhi(v02), + o2::constants::physics::MassProton); + lambda2 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(v02), mcacc::lamEta(v02), mcacc::lamPhi(v02), + mcacc::lamMass(v02)); + + histos.fill(HIST("deltaPhiSame"), + RecoDecay::constrainAngle(mcacc::lamPhi(v0) - mcacc::lamPhi(v02), + -TMath::Pi(), harmonicDphi)); + + const int s1 = mcacc::v0Status(v0); + const int s2 = mcacc::v0Status(v02); + + if (s1 == 0 && s2 == 0) { + fillHistograms(0, 0, lambda, lambda2, proton, proton2, 0, 1.0f); + } else if (s1 == 0 && s2 == 1) { + fillHistograms(0, 1, lambda, lambda2, proton, proton2, 0, 1.0f); + } else if (s1 == 1 && s2 == 0) { + fillHistograms(0, 1, lambda2, lambda, proton2, proton, 0, 1.0f); + } else if (s1 == 1 && s2 == 1) { + fillHistograms(1, 1, lambda, lambda2, proton, proton2, 0, 1.0f); + } + } + } + } + PROCESS_SWITCH(lambdaspincorrderived, processMC, "Process MC (SE)", 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 phi with periodic wrap at 0/2pi + static inline void collectNeighborBinsPhi(int b, int nPhi, int nNeighbor, std::vector& out) + { + out.clear(); + if (nPhi <= 0 || b < 0 || b >= nPhi) { + return; + } + + if (2 * nNeighbor + 1 >= nPhi) { + out.reserve(nPhi); + for (int bb = 0; bb < nPhi; ++bb) { + out.push_back(bb); + } + return; + } + + out.reserve(2 * nNeighbor + 1); + for (int d = -nNeighbor; d <= nNeighbor; ++d) { + int bb = (b + d) % 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()); + } + + static inline void collectNeighborBinsClamp(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) { + out.push_back(bb); + } + } + } + static inline void collectNeighborBinsMass(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) { + out.push_back(bb); + } + } + + std::sort(out.begin(), out.end()); + out.erase(std::unique(out.begin(), out.end()), out.end()); + } + + static inline int getMassRegionFromEdges(float m, const std::vector& edges) + { + if (edges.size() != 4) { + return -1; + } + + if (m >= edges[0] && m < edges[1]) + return 0; // low sideband + if (m >= edges[1] && m < edges[2]) + return 1; // signal + if (m >= edges[2] && m < edges[3]) + return 2; // high sideband + + return -1; + } + static inline int getMassMixClassFromEdges(float m, const std::vector& edges) + { + // no mass separation + if (edges.size() == 2) { + if (m >= edges[0] && m < edges[1]) + return 0; + return -1; + } + + // 3 regions: SB low, signal, SB high + if (edges.size() == 4) { + if (m >= edges[0] && m < edges[1]) + return 0; // low sideband + if (m >= edges[1] && m < edges[2]) + return 1; // signal + if (m >= edges[2] && m < edges[3]) + return 0; // high sideband merged with low sideband + return -1; + } + + return -1; + } + static inline uint64_t splitmix64(uint64_t x) + { + // simple deterministic hash for reproducible shuffling + x += 0x9e3779b97f4a7c15ULL; + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + return x ^ (x >> 31); + } + + void processMEV6(EventCandidates const& collisions, AllTrackCandidates const& V0s) + { + MixBinnerR mb{ + ptMin.value, + ptMax.value, + ptMix.value, + v0etaMixBuffer.value, + etaMix.value, + phiMix.value, + MassMin.value, + MassMax.value, + cfgV5MassBins.value, + cfgMixRadiusParam.cfgMixRadiusBins.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 int nR = mb.nR(); + + const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM * nR; + std::vector> buffer(nKeys); - if (rEdges.size() < 2) { - rEdges = {3.0, 5.0, 7.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0}; + // -------- 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; } - nR_ = static_cast(rEdges.size()) - 1; - 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))); - } + auto slice = V0s.sliceBy(tracksPerCollisionV0, col.index()); - inline int nPt() const { return nPt_; } - inline int nEta() const { return nEta_; } - inline int nPhi() const { return nPhi_; } - inline int nM() const { return nM_; } - inline int nR() const { return nR_; } + for (auto const& t : slice) { + if (!selectionV0(t)) { + continue; + } - inline int binFromValue(float v, float vmin, float step, int nBins) const - { - 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) { - return -1; - } - if (b >= nBins) { - b = nBins - 1; - } - return b; - } + const int status = static_cast(t.v0Status()); + if (status < 0 || status >= nStat) { + continue; + } - inline int ptBin(float pt) const { return binFromValue(pt, ptMin, ptStep, nPt_); } - inline int etaBin(float eta) const { return binFromValue(eta, etaMin, etaStep, nEta_); } - inline int phiBin(float phi) const { return binFromValue(phi, phiMin, phiStep, nPhi_); } - inline int radiusBin(float r) const - { - if (!std::isfinite(r) || nR_ <= 0) { - return -1; - } - if (r < rEdges.front() || r >= rEdges.back()) { - return -1; + 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(), 0.0, harmonic)); + const int mB = getMassMixClassFromEdges(t.lambdaMass(), massMixEdges.value); + const int rB = mb.radiusBin(t.v0Radius()); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0 || rB < 0) { + continue; + } + + const size_t key = linearKeyR(colBin, status, ptB, etaB, phiB, mB, rB, + nStat, nPt, nEta, nPhi, nM, nR); + + buffer[key].push_back(BufferCandR{ + .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), + .rBin = static_cast(rB)}); } - auto it = std::upper_bound(rEdges.begin(), rEdges.end(), static_cast(r)); - return static_cast(it - rEdges.begin()) - 1; } - }; - - struct BufferCandR { - int64_t collisionIdx; - int64_t rowIndex; - uint8_t v0Status; - uint16_t ptBin, etaBin, phiBin, mBin, rBin; - }; - struct StoredV6Candidate { - int64_t collisionIdx = -1; - int64_t globalIdx = -1; - int status = -1; - bool isDouble = false; - float cospa = 0.f; - float radius = 0.f; - float dcaPos = 0.f; - float dcaNeg = 0.f; - float dcaDau = 0.f; - float lPt = 0.f; - float lEta = 0.f; - float lPhi = 0.f; - float lMass = 0.f; - float pPt = 0.f; - float pEta = 0.f; - float pPhi = 0.f; - int64_t pIndex = -1; - int64_t piIndex = -1; + 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); - int v0Status() const { return status; } - bool doubleStatus() const { return isDouble; } - float v0Cospa() const { return cospa; } - float v0Radius() const { return radius; } - float dcaPositive() const { return dcaPos; } - float dcaNegative() const { return dcaNeg; } - float dcaBetweenDaughter() const { return dcaDau; } - float lambdaPt() const { return lPt; } - float lambdaEta() const { return lEta; } - float lambdaPhi() const { return lPhi; } - float lambdaMass() const { return lMass; } - float protonPt() const { return pPt; } - float protonEta() const { return pEta; } - float protonPhi() const { return pPhi; } - int64_t protonIndex() const { return pIndex; } - int64_t pionIndex() const { return piIndex; } - int64_t globalIndex() const { return globalIdx; } - }; + std::vector ptBins, etaBins, phiBins; + std::vector matches1, matches2; + matches1.reserve(256); + matches2.reserve(256); - struct PendingV6Branch { - StoredV6Candidate target; - StoredV6Candidate fixed; - int colBin = -1; - int replacedLeg = 1; - int age = 0; - uint64_t seed = 0; - }; + auto collectMatchesForReplacedLeg = [&](auto const& tRep, auto const& tKeep, int colBin, int64_t curColIdx, std::vector& matches) { + matches.clear(); - template - StoredV6Candidate storeV6Candidate(T const& t, int64_t collisionIdx) const - { - return {collisionIdx, static_cast(t.globalIndex()), static_cast(t.v0Status()), static_cast(t.doubleStatus()), - t.v0Cospa(), t.v0Radius(), t.dcaPositive(), t.dcaNegative(), t.dcaBetweenDaughter(), - t.lambdaPt(), t.lambdaEta(), t.lambdaPhi(), t.lambdaMass(), - t.protonPt(), t.protonEta(), t.protonPhi(), - static_cast(t.protonIndex()), static_cast(t.pionIndex())}; - } + const int status = static_cast(tRep.v0Status()); + if (status < 0 || status >= nStat) { + return; + } - template - void fillV6MixedBranch(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) - { - const auto repProton = ROOT::Math::PtEtaPhiMVector(replacement.protonPt(), replacement.protonEta(), replacement.protonPhi(), o2::constants::physics::MassProton); - const auto repLambda = ROOT::Math::PtEtaPhiMVector(replacement.lambdaPt(), replacement.lambdaEta(), replacement.lambdaPhi(), replacement.lambdaMass()); - const auto fixedProton = ROOT::Math::PtEtaPhiMVector(fixed.protonPt(), fixed.protonEta(), fixed.protonPhi(), o2::constants::physics::MassProton); - const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(fixed.lambdaPt(), fixed.lambdaEta(), fixed.lambdaPhi(), fixed.lambdaMass()); - const int repStatus = replacement.v0Status(); - const int fixedStatus = fixed.v0Status(); + const int ptB = mb.ptBin(tRep.lambdaPt()); - if (replacedLeg == 1) { - fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); - fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); - if ((repStatus == 0 && fixedStatus == 1) || (repStatus == 1 && fixedStatus == 0)) { - if (fillBasicQAHistos) { - histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); - } + 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()); } - if (repStatus == 0 && fixedStatus == 1) { - fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); - } else if (repStatus == 1 && fixedStatus == 0) { - fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); - } else { - fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); + + const int phiB = mb.phiBin(RecoDecay::constrainAngle(tRep.lambdaPhi(), 0.0, harmonic)); + const int mB = getMassMixClassFromEdges(tRep.lambdaMass(), massMixEdges.value); + const int rB = mb.radiusBin(tRep.v0Radius()); + + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0 || rB < 0) { + return; } - return; - } + auto collectFromBins = [&](const std::vector& ptUseBins, + const std::vector& etaUseBins, + const std::vector& phiUseBins) { + for (const auto& ptUse : ptUseBins) { + for (const auto& etaUse : etaUseBins) { + for (const auto& phiUse : phiUseBins) { + const auto& vec = buffer[linearKeyR(colBin, status, ptUse, etaUse, phiUse, mB, rB, + nStat, nPt, nEta, nPhi, nM, nR)]; - fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); - fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); - if (fillBasicQAHistos) { - histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); - } - if (fixedStatus == 0 && repStatus == 1) { - fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); - } else if (fixedStatus == 1 && repStatus == 0) { - fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); - } else { - fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); - } - } + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) + continue; - struct StoredV6CandidateMC { - int64_t globalIdx = -1; - int status = -1; - bool isDouble = false; - float cospa = 0.f; - float radius = 0.f; - float dcaPos = 0.f; - float dcaNeg = 0.f; - float dcaDau = 0.f; - float lPt = 0.f; - float lEta = 0.f; - float lPhi = 0.f; - float lMass = 0.f; - float pPt = 0.f; - float pEta = 0.f; - float pPhi = 0.f; - int64_t pIndex = -1; - int64_t piIndex = -1; + auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); + + if (!selectionV0(tX)) + continue; + if (!checkKinematics(tRep, tX)) + continue; - int v0Statusmc() const { return status; } - bool doubleStatusmc() const { return isDouble; } - float v0Cospamc() const { return cospa; } - float v0Radiusmc() const { return radius; } - float dcaPositivemc() const { return dcaPos; } - float dcaNegativemc() const { return dcaNeg; } - float dcaBetweenDaughtermc() const { return dcaDau; } - float lambdaPtmc() const { return lPt; } - float lambdaEtamc() const { return lEta; } - float lambdaPhimc() const { return lPhi; } - float lambdaMassmc() const { return lMass; } - float protonPtmc() const { return pPt; } - float protonEtamc() const { return pEta; } - float protonPhimc() const { return pPhi; } - int64_t protonIndexmc() const { return pIndex; } - int64_t pionIndexmc() const { return piIndex; } - int64_t globalIndex() const { return globalIdx; } - }; + if (tX.globalIndex() == tRep.globalIndex()) + continue; + if (tX.globalIndex() == tKeep.globalIndex()) + continue; - struct PendingV6BranchMC { - StoredV6CandidateMC target; - StoredV6CandidateMC fixed; - int colBin = -1; - int replacedLeg = 1; - int age = 0; - uint64_t seed = 0; - }; + if (hasSharedDaughters(tX, tKeep)) + continue; + if (hasSharedDaughters(tX, tRep)) + continue; - struct V6PendingState { - std::deque data; - std::deque mc; - }; + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } + }; - V6PendingState v6Pending; + matches.clear(); - template - StoredV6CandidateMC storeV6CandidateMC(T const& t) const - { - return {static_cast(t.globalIndex()), mcacc::v0Status(t), mcacc::doubleStatus(t), - mcacc::v0CosPA(t), mcacc::v0Radius(t), mcacc::dcaPos(t), mcacc::dcaNeg(t), mcacc::dcaDau(t), - mcacc::lamPt(t), mcacc::lamEta(t), mcacc::lamPhi(t), mcacc::lamMass(t), - mcacc::prPt(t), mcacc::prEta(t), mcacc::prPhi(t), - static_cast(mcacc::prIdx(t)), static_cast(mcacc::piIdx(t))}; - } + // 1) exact bin first + ptBins.clear(); + etaBins.clear(); + phiBins.clear(); - template - void fillV6MixedBranchMC(TRep const& replacement, TFixed const& fixed, int replacedLeg, float controlWeight, float mixWeight) - { - const auto repProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(replacement), mcacc::prEta(replacement), mcacc::prPhi(replacement), o2::constants::physics::MassProton); - const auto repLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(replacement), mcacc::lamEta(replacement), mcacc::lamPhi(replacement), mcacc::lamMass(replacement)); - const auto fixedProton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(fixed), mcacc::prEta(fixed), mcacc::prPhi(fixed), o2::constants::physics::MassProton); - const auto fixedLambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(fixed), mcacc::lamEta(fixed), mcacc::lamPhi(fixed), mcacc::lamMass(fixed)); - const int repStatus = mcacc::v0Status(replacement); - const int fixedStatus = mcacc::v0Status(fixed); + ptBins.push_back(ptB); + etaBins.push_back(etaB); + phiBins.push_back(phiB); - if (replacedLeg == 1) { - fillReplacementControlMap(repStatus, fixedStatus, 1, false, repLambda, controlWeight); - fillFixedLegControlMap(repStatus, fixedStatus, 1, false, fixedLambda, controlWeight); - if (fillBasicQAHistos) { - histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)repLambda.Phi(), (float)fixedLambda.Phi()), mixWeight); - } - if (repStatus == 0 && fixedStatus == 1) { - fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); - } else if (repStatus == 1 && fixedStatus == 0) { - fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 1); - } else { - fillHistograms(repStatus, fixedStatus, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 1); - } - return; - } + collectFromBins(ptBins, etaBins, phiBins); - fillReplacementControlMap(fixedStatus, repStatus, 2, false, repLambda, controlWeight); - fillFixedLegControlMap(fixedStatus, repStatus, 2, false, fixedLambda, controlWeight); - if (fillBasicQAHistos) { - histos.fill(HIST("deltaPhiMix"), deltaPhiMinusPiToPi((float)fixedLambda.Phi(), (float)repLambda.Phi()), mixWeight); - } - if (fixedStatus == 0 && repStatus == 1) { - fillHistograms(0, 1, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); - } else if (fixedStatus == 1 && repStatus == 0) { - fillHistograms(0, 1, repLambda, fixedLambda, repProton, fixedProton, 1, mixWeight, 1, 2); - } else { - fillHistograms(fixedStatus, repStatus, fixedLambda, repLambda, fixedProton, repProton, 1, mixWeight, 2, 2); - } - } + // 2) if exact bin gives fewer than required matches, also search neighbors + const int targetMatches = (cfgV5MaxMatches.value > 0) ? cfgV5MaxMatches.value : 1; - static inline size_t linearKeyR(int colBin, int statBin, - int ptBin, int etaBin, int phiBin, int mBin, int rBin, - int nStatus, int nPt, int nEta, int nPhi, int nM, int nR) - { - return (((((((static_cast(colBin) * nStatus + statBin) * nPt + ptBin) * nEta + etaBin) * nPhi + phiBin) * nM + mBin) * nR + rBin)); - } + if ((int)matches.size() < targetMatches) { + std::vector ptBinsN, etaBinsN, phiBinsN; + collectNeighborBinsClamp(ptB, nPt, nN_pt, ptBinsN); + collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBinsN); + collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBinsN); - // ------------------------------------- - // 2) MC-only selection + kinematics cuts - // ------------------------------------- - template - bool selectionV0MC(T const& candidate) - { - auto particle = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate), - mcacc::lamEta(candidate), - mcacc::lamPhi(candidate), - mcacc::lamMass(candidate)); - if (std::abs(particle.Rapidity()) > rapidity || std::abs(particle.Eta()) > v0eta) { - return false; - } - if (mcacc::lamMass(candidate) < MassMin || mcacc::lamMass(candidate) > MassMax) { - return false; - } - if (mcacc::v0CosPA(candidate) < v0Configurations.cosPA) { - return false; - } - if (checkDoubleStatus && mcacc::doubleStatus(candidate)) { - return false; - } - if (mcacc::v0Radius(candidate) > v0Configurations.radiusMax) { - return false; - } - if (mcacc::v0Radius(candidate) < v0Configurations.radiusMin) { - return false; - } - if (mcacc::dcaDau(candidate) > v0Configurations.dcaDaughters) { - return false; - } + for (const auto& ptUse : ptBinsN) { + for (const auto& etaUse : etaBinsN) { + for (const auto& phiUse : phiBinsN) { + if (ptUse == ptB && etaUse == etaB && phiUse == phiB) + continue; - if (mcacc::dcaV0ToPVMC(candidate) > v0Configurations.dcaV0ToPV) { - return false; - } + const auto& vec = buffer[linearKeyR(colBin, status, ptUse, etaUse, phiUse, mB, rB, + nStat, nPt, nEta, nPhi, nM, nR)]; - if (mcacc::v0Status(candidate) == 0 && (std::abs(mcacc::dcaPos(candidate)) < v0Configurations.dcaProton || std::abs(mcacc::dcaNeg(candidate)) < v0Configurations.dcaPion)) { - return false; - } - if (mcacc::v0Status(candidate) == 1 && (std::abs(mcacc::dcaPos(candidate)) < v0Configurations.dcaPion || std::abs(mcacc::dcaNeg(candidate)) < v0Configurations.dcaProton)) { - return false; - } - if (mcacc::lamPt(candidate) < ptMin) { - return false; - } - if (mcacc::lamPt(candidate) > ptMax) { - return false; - } - return true; - } + for (auto const& bc : vec) { + if (bc.collisionIdx == curColIdx) + continue; - template - bool checkKinematicsMC(T1 const& candidate1, T2 const& candidate2) - { - // keep same species/status - if (mcacc::v0Status(candidate1) != mcacc::v0Status(candidate2)) { - return false; - } + auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); - // pT window - if (std::abs(mcacc::lamPt(candidate1) - mcacc::lamPt(candidate2)) > ptMix) { - return false; - } + if (!selectionV0(tX)) + continue; + if (!checkKinematics(tRep, tX)) + continue; - // eta or rapidity window (etaMix used as Δη or Δy) - if (!userapidity) { - if (std::abs(mcacc::lamEta(candidate1) - mcacc::lamEta(candidate2)) > etaMix) { - return false; - } - } else { - const auto l1 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate1), mcacc::lamEta(candidate1), - mcacc::lamPhi(candidate1), mcacc::lamMass(candidate1)); - const auto l2 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(candidate2), mcacc::lamEta(candidate2), - mcacc::lamPhi(candidate2), mcacc::lamMass(candidate2)); - if (std::abs(l1.Rapidity() - l2.Rapidity()) > etaMix) { - return false; + if (tX.globalIndex() == tRep.globalIndex()) + continue; + if (tX.globalIndex() == tKeep.globalIndex()) + continue; + + if (hasSharedDaughters(tX, tKeep)) + continue; + if (hasSharedDaughters(tX, tRep)) + continue; + + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + } + } + } } - } - // delta-phi window (wrapped) - const float dphi = deltaPhiMinusPiToPi((float)mcacc::lamPhi(candidate1), - (float)mcacc::lamPhi(candidate2)); - if (std::abs(dphi) > phiMix) { - return false; - } + 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()); + }; - // mass window (optional but consistent with data) - if (std::abs(mcacc::lamMass(candidate1) - mcacc::lamMass(candidate2)) > massMix) { - return false; - } + 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); + } + }; - return true; - } + const size_t pendingAtStart = v6Pending.data.size(); + size_t pendingMatched = 0; + size_t pendingExpired = 0; + size_t pendingAdded = 0; + if (!cfgV6CarryUnmatched) { + v6Pending.data.clear(); + } else { + for (auto it = v6Pending.data.begin(); it != v6Pending.data.end();) { + auto& pending = *it; + ++pending.age; + if (cfgV6MaxPendingAge.value > 0 && pending.age > cfgV6MaxPendingAge.value) { + ++pendingExpired; + it = v6Pending.data.erase(it); + continue; + } - // ----------------------------------------- - // 3) MC filter + aliases (distinct from data) - // ----------------------------------------- - Filter centralityFilterMC = (nabs(aod::lambdaeventmc::centmc) < centMax && nabs(aod::lambdaeventmc::centmc) > centMin); + auto& matches = pending.replacedLeg == 1 ? matches1 : matches2; + collectMatchesForReplacedLeg(pending.target, pending.fixed, pending.colBin, -1, matches); + limitMatchesToNEvents(matches, nEvtMixing.value); + downsampleMatches(matches, pending.seed ^ splitmix64(static_cast(pending.age))); - using EventCandidatesMC = soa::Filtered; - using AllTrackCandidatesMC = aod::LambdaPairmcs; + int nAccepted = 0; + for (auto const& m : matches) { + auto replacement = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(replacement) || !checkKinematics(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughters(replacement, pending.target) || hasSharedDaughters(replacement, pending.fixed)) { + continue; + } + ++nAccepted; + } - // IMPORTANT: MC preslice uses the MC event index column - Preslice tracksPerCollisionV0mc = aod::lambdapairmc::lambdaeventmcId; + if (nAccepted == 0) { + ++it; + continue; + } - // ----------------------------------------- - // 4) MC Same-event processing (like processData) - // ----------------------------------------- - void processMC(EventCandidatesMC::iterator const& collision, AllTrackCandidatesMC const& V0sMC) - { - const float centrality = mcacc::cent(collision); + const float controlWeight = 1.0f / static_cast(nAccepted); + const float branchNorm = cfgMixLegMode.value == 2 ? 0.5f : 1.0f; + const float mixWeight = branchNorm * controlWeight; + for (auto const& m : matches) { + auto replacement = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(replacement) || !checkKinematics(pending.target, replacement)) { + continue; + } + if (replacement.globalIndex() == pending.target.globalIndex() || replacement.globalIndex() == pending.fixed.globalIndex()) { + continue; + } + if (hasSharedDaughters(replacement, pending.target) || hasSharedDaughters(replacement, pending.fixed)) { + continue; + } + fillV6MixedBranch(replacement, pending.fixed, pending.replacedLeg, controlWeight, mixWeight); + } + ++pendingMatched; + it = v6Pending.data.erase(it); + } + } - for (const auto& v0 : V0sMC) { - if (!selectionV0MC(v0)) { + // -------- PASS 2: configurable one-leg / two-leg mixing -------- + for (auto const& col1 : collisions) { + const int colBin = colBinning.getBin(std::make_tuple(col1.posz(), col1.cent())); + if (colBin < 0) { continue; } - if (fillBasicQAHistos) { - histos.fill(HIST("hPtRadiusV0"), mcacc::lamPt(v0), mcacc::v0Radius(v0)); - } - if (fillBasicQAHistos) { - histos.fill(HIST("ptCent"), mcacc::lamPt(v0), centrality); - } - if (fillBasicQAHistos) { - histos.fill(HIST("etaCent"), mcacc::lamEta(v0), centrality); - } - proton = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(v0), mcacc::prEta(v0), mcacc::prPhi(v0), - o2::constants::physics::MassProton); - lambda = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(v0), mcacc::lamEta(v0), mcacc::lamPhi(v0), - mcacc::lamMass(v0)); - if (mcacc::v0Status(v0) == 0) { - histos.fill(HIST("hEtaPhiLambdaRaw"), lambda.Phi(), lambda.Eta(), getNUAWeight(0, lambda.Phi(), lambda.Eta())); - } else { - histos.fill(HIST("hEtaPhiAntiLambdaRaw"), lambda.Phi(), lambda.Eta(), getNUAWeight(1, lambda.Phi(), lambda.Eta())); - } + const int64_t curColIdx = static_cast(col1.index()); + auto poolA = V0s.sliceBy(tracksPerCollisionV0, col1.index()); - for (const auto& v02 : V0sMC) { - if (v02.index() <= v0.index()) { + for (auto const& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { + if (!selectionV0(t1) || !selectionV0(t2)) { continue; } - if (!selectionV0MC(v02)) { + if (t2.index() <= t1.index()) { continue; } - if (hasSharedDaughtersMC(v0, v02)) + if (hasSharedDaughters(t1, t2)) continue; - proton2 = ROOT::Math::PtEtaPhiMVector(mcacc::prPt(v02), mcacc::prEta(v02), mcacc::prPhi(v02), - o2::constants::physics::MassProton); - lambda2 = ROOT::Math::PtEtaPhiMVector(mcacc::lamPt(v02), mcacc::lamEta(v02), mcacc::lamPhi(v02), - mcacc::lamMass(v02)); - - histos.fill(HIST("deltaPhiSame"), - RecoDecay::constrainAngle(mcacc::lamPhi(v0) - mcacc::lamPhi(v02), - -TMath::Pi(), harmonicDphi)); - - const int s1 = mcacc::v0Status(v0); - const int s2 = mcacc::v0Status(v02); + const bool doMixLeg1 = (cfgMixLegMode.value == 0 || cfgMixLegMode.value == 2); + const bool doMixLeg2 = (cfgMixLegMode.value == 1 || cfgMixLegMode.value == 2); - if (s1 == 0 && s2 == 0) { - fillHistograms(0, 0, lambda, lambda2, proton, proton2, 0, 1.0f); - } else if (s1 == 0 && s2 == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 0, 1.0f); - } else if (s1 == 1 && s2 == 0) { - fillHistograms(0, 1, lambda2, lambda, proton2, proton, 0, 1.0f); - } else if (s1 == 1 && s2 == 1) { - fillHistograms(1, 1, lambda, lambda2, proton, proton2, 0, 1.0f); + // Fill TGT maps before searching for replacements. This makes TGT the + // true same-event target phase space for the selected pair, matching the + // hPtYSame definition. REP below is filled only for accepted replacements. + if (doMixLeg1) { + fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 1, true, + ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), + 1.0f); + fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 1, true, + ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), + 1.0f); + } + if (doMixLeg2) { + fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 2, true, + ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), + 1.0f); + fillFixedLegControlMap(t1.v0Status(), t2.v0Status(), 2, true, + ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), + 1.0f); } - } - } - } - PROCESS_SWITCH(lambdaspincorrderived, processMC, "Process MC (SE)", false); - static inline float phi0To2Pi(float phi) - { - // harmonic=1, min=0 => [0, 2pi) - return RecoDecay::constrainAngle(phi, 0.0f, 1); - } + if (doMixLeg1) { + collectMatchesForReplacedLeg(t1, t2, colBin, curColIdx, matches1); + limitMatchesToNEvents(matches1, nEvtMixing.value); + downsampleMatches(matches1, (uint64_t)t1.globalIndex() ^ (splitmix64((uint64_t)t2.globalIndex()) + 0x111ULL) ^ splitmix64((uint64_t)curColIdx)); + } else { + matches1.clear(); + } - 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); - } + if (doMixLeg2) { + collectMatchesForReplacedLeg(t2, t1, colBin, curColIdx, matches2); + limitMatchesToNEvents(matches2, nEvtMixing.value); + downsampleMatches(matches2, (uint64_t)t2.globalIndex() ^ (splitmix64((uint64_t)t1.globalIndex()) + 0x222ULL) ^ splitmix64((uint64_t)curColIdx)); + } else { + matches2.clear(); + } - static inline float absDeltaPhi(float phiA, float phiB) - { - return std::abs(deltaPhiMinusPiToPi(phiA, phiB)); - } + // Do not fill TGT here. TGT has already been filled above as the + // selected same-event target phase space, independent of replacement success. + int nFill1 = 0; + int nFill2 = 0; + // count actual accepted fills for leg-1 replacement + if (doMixLeg1) { + for (auto const& m : matches1) { + auto tX = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(tX)) + continue; + if (tX.v0Status() != t1.v0Status()) + continue; + if (!checkKinematics(t1, tX)) + continue; + if (tX.globalIndex() == t1.globalIndex()) + continue; + if (tX.globalIndex() == t2.globalIndex()) + continue; + if (hasSharedDaughters(tX, t2)) + continue; + if (hasSharedDaughters(tX, t1)) + continue; + ++nFill1; + } + } - // symmetric neighbors for phi with periodic wrap at 0/2pi - static inline void collectNeighborBinsPhi(int b, int nPhi, int nNeighbor, std::vector& out) - { - out.clear(); - if (nPhi <= 0 || b < 0 || b >= nPhi) { - return; - } + // count actual accepted fills for leg-2 replacement + if (doMixLeg2) { + for (auto const& m : matches2) { + auto tY = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(tY)) + continue; + if (tY.v0Status() != t2.v0Status()) + continue; + if (!checkKinematics(t2, tY)) + continue; + if (tY.globalIndex() == t2.globalIndex()) + continue; + if (tY.globalIndex() == t1.globalIndex()) + continue; + if (hasSharedDaughters(tY, t2)) + continue; + if (hasSharedDaughters(tY, t1)) + continue; + ++nFill2; + } + } - if (2 * nNeighbor + 1 >= nPhi) { - out.reserve(nPhi); - for (int bb = 0; bb < nPhi; ++bb) { - out.push_back(bb); - } - return; - } + if (cfgV6CarryUnmatched) { + const auto hasPendingSpace = [&]() { + return cfgV6MaxPendingBranches.value <= 0 || static_cast(v6Pending.data.size()) < cfgV6MaxPendingBranches.value; + }; + if (doMixLeg1 && nFill1 == 0 && hasPendingSpace()) { + v6Pending.data.push_back({storeV6Candidate(t1, curColIdx), storeV6Candidate(t2, curColIdx), colBin, 1, 0, + static_cast(t1.globalIndex()) ^ splitmix64(static_cast(t2.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAdded; + } + if (doMixLeg2 && nFill2 == 0 && hasPendingSpace()) { + v6Pending.data.push_back({storeV6Candidate(t2, curColIdx), storeV6Candidate(t1, curColIdx), colBin, 2, 0, + static_cast(t2.globalIndex()) ^ splitmix64(static_cast(t1.globalIndex())) ^ splitmix64(static_cast(curColIdx))}); + ++pendingAdded; + } + } - out.reserve(2 * nNeighbor + 1); - for (int d = -nNeighbor; d <= nNeighbor; ++d) { - int bb = (b + d) % nPhi; - if (bb < 0) { - bb += nPhi; - } - out.push_back(bb); - } + if (nFill1 <= 0 && nFill2 <= 0) { + continue; + } + // Residual-weight QA needs a leg-specific normalization: + // TGT_leg is filled once per same-event target candidate with weight 1. + // REP_leg is the average replacement distribution for that same target. + const float wSELeg1 = (nFill1 > 0) ? 1.0f / static_cast(nFill1) : 0.0f; + const float wSELeg2 = (nFill2 > 0) ? 1.0f / static_cast(nFill2) : 0.0f; - std::sort(out.begin(), out.end()); - out.erase(std::unique(out.begin(), out.end()), out.end()); - } + const int nActiveMixBranches = ((doMixLeg1 && nFill1 > 0) ? 1 : 0) + + ((doMixLeg2 && nFill2 > 0) ? 1 : 0); + float branchNorm = 1.0f; + if (cfgMixLegMode.value == 2) { + branchNorm = cfgV6CarryUnmatched ? 0.5f : 1.0f / static_cast(nActiveMixBranches); + } + const float finalMixWeightLeg1 = branchNorm * wSELeg1; + const float finalMixWeightLeg2 = branchNorm * wSELeg2; - static inline void collectNeighborBinsClamp(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) { - out.push_back(bb); - } - } - } - static inline void collectNeighborBinsMass(int b, int nBins, int nNeighbor, std::vector& out) - { - out.clear(); - out.reserve(2 * nNeighbor + 1); + // Do not fill target maps here: TGT has already been filled above, + // before searching for replacements, so that TGT represents all selected + // same-event targets rather than only targets with successful replacements. - for (int d = -nNeighbor; d <= nNeighbor; ++d) { - const int bb = b + d; - if (bb >= 0 && bb < nBins) { - out.push_back(bb); - } - } + if (doMixLeg1 && nFill1 > 0) { + for (auto const& m : matches1) { + auto tX = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(tX)) { + continue; + } + if (tX.v0Status() != t1.v0Status()) + continue; + if (!checkKinematics(t1, tX)) + continue; + if (tX.globalIndex() == t1.globalIndex()) + continue; + if (tX.globalIndex() == t2.globalIndex()) + continue; + if (hasSharedDaughters(tX, t1)) + continue; + if (hasSharedDaughters(tX, t2)) + continue; + fillReplacementControlMap(tX.v0Status(), t2.v0Status(), 1, false, + ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), + wSELeg1); + // Fixed leg for leg1 replacement is original t2. + // Fill only for successful replacement. + fillFixedLegControlMap(tX.v0Status(), t2.v0Status(), 1, false, + ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), + wSELeg1); + 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()); - std::sort(out.begin(), out.end()); - out.erase(std::unique(out.begin(), out.end()), out.end()); - } + // const int ptype = pairTypeCode(tX.v0Status(), t2.v0Status()); - static inline int getMassRegionFromEdges(float m, const std::vector& edges) - { - if (edges.size() != 4) { - return -1; - } + const float meWeight = finalMixWeightLeg1; + const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); + if ((tX.v0Status() == 0 && t2.v0Status() == 1) || (tX.v0Status() == 1 && t2.v0Status() == 0)) + if (fillBasicQAHistos) + histos.fill(HIST("deltaPhiMix"), dPhi, meWeight); + const int s1 = tX.v0Status(); + const int s2 = t2.v0Status(); - if (m >= edges[0] && m < edges[1]) - return 0; // low sideband - if (m >= edges[1] && m < edges[2]) - return 1; // signal - if (m >= edges[2] && m < edges[3]) - return 2; // high sideband + if (s1 == 0 && s2 == 1) { + fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, 1, 1); + } else if (s1 == 1 && s2 == 0) { + fillHistograms(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, 2, 1); + } else { + fillHistograms(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, 1, 1); + } + } + } - return -1; - } - static inline int getMassMixClassFromEdges(float m, const std::vector& edges) - { - // no mass separation - if (edges.size() == 2) { - if (m >= edges[0] && m < edges[1]) - return 0; - return -1; - } + if (doMixLeg2 && nFill2 > 0) { + for (auto const& m : matches2) { + auto tY = V0s.iteratorAt(static_cast(m.rowIndex)); + if (!selectionV0(tY)) { + continue; + } + if (tY.v0Status() != t2.v0Status()) + continue; + if (!checkKinematics(t2, tY)) + continue; + if (tY.globalIndex() == t2.globalIndex()) + continue; + if (tY.globalIndex() == t1.globalIndex()) + continue; + if (hasSharedDaughters(tY, t1)) + continue; + if (hasSharedDaughters(tY, t2)) + continue; + fillReplacementControlMap(t1.v0Status(), tY.v0Status(), 2, false, + ROOT::Math::PtEtaPhiMVector(tY.lambdaPt(), tY.lambdaEta(), tY.lambdaPhi(), tY.lambdaMass()), + wSELeg2); + // Fixed leg for leg2 replacement is original t1. + // Fill only for successful replacement. + fillFixedLegControlMap(t1.v0Status(), tY.v0Status(), 2, false, + ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), + wSELeg2); + 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()); - // 3 regions: SB low, signal, SB high - if (edges.size() == 4) { - if (m >= edges[0] && m < edges[1]) - return 0; // low sideband - if (m >= edges[1] && m < edges[2]) - return 1; // signal - if (m >= edges[2] && m < edges[3]) - return 0; // high sideband merged with low sideband - return -1; + // const int ptype = pairTypeCode(t1.v0Status(), tY.v0Status()); + const float meWeight = finalMixWeightLeg2; + const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); + if (fillBasicQAHistos) + histos.fill(HIST("deltaPhiMix"), dPhi, meWeight); + const int s1 = t1.v0Status(); + const int s2 = tY.v0Status(); + if (s1 == 0 && s2 == 1) { + fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, 2, 2); + } else if (s1 == 1 && s2 == 0) { + fillHistograms(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, 1, 2); + } else { + fillHistograms(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, 2, 2); + } + } + } + } + } + if (cfgV6LogPending) { + LOGF(info, "MEV6 data pending branches: carriedIn=%zu matched=%zu expired=%zu newlyPropagated=%zu carriedToNext=%zu", + pendingAtStart, pendingMatched, pendingExpired, pendingAdded, v6Pending.data.size()); } - - return -1; - } - static inline uint64_t splitmix64(uint64_t x) - { - // simple deterministic hash for reproducible shuffling - x += 0x9e3779b97f4a7c15ULL; - x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; - x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; - return x ^ (x >> 31); } + PROCESS_SWITCH(lambdaspincorrderived, processMEV6, "Process data ME v6 with radius buffer", false); - void processMEV6(EventCandidates const& collisions, AllTrackCandidates const& V0s) + void processMEV6Sys(EventCandidates const& collisions, AllTrackCandidates const& V0s) { MixBinnerR mb{ ptMin.value, @@ -2165,9 +3161,16 @@ struct lambdaspincorrderived { auto slice = V0s.sliceBy(tracksPerCollisionV0, col.index()); for (auto const& t : slice) { - if (!selectionV0(t)) { - continue; + + bool passedAny = false; + for (int sysId = 0; sysId < nSysTotal; ++sysId) { + if (selectionV0Sys(t, sysCuts[sysId])) { + passedAny = true; + break; + } } + if (!passedAny) + continue; const int status = static_cast(t.v0Status()); if (status < 0 || status >= nStat) { @@ -2240,9 +3243,9 @@ struct lambdaspincorrderived { auto collectFromBins = [&](const std::vector& ptUseBins, const std::vector& etaUseBins, const std::vector& phiUseBins) { - for (int ptUse : ptUseBins) { - for (int etaUse : etaUseBins) { - for (int phiUse : phiUseBins) { + for (const auto& ptUse : ptUseBins) { + for (const auto& etaUse : etaUseBins) { + for (const auto& phiUse : phiUseBins) { const auto& vec = buffer[linearKeyR(colBin, status, ptUse, etaUse, phiUse, mB, rB, nStat, nPt, nEta, nPhi, nM, nR)]; @@ -2252,8 +3255,16 @@ struct lambdaspincorrderived { auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); - if (!selectionV0(tX)) + bool tXPassed = false; + for (int sysId = 0; sysId < nSysTotal; ++sysId) { + if (selectionV0Sys(tX, sysCuts[sysId])) { + tXPassed = true; + break; + } + } + if (!tXPassed) continue; + if (!checkKinematics(tRep, tX)) continue; @@ -2296,9 +3307,9 @@ struct lambdaspincorrderived { collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBinsN); collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBinsN); - for (int ptUse : ptBinsN) { - for (int etaUse : etaBinsN) { - for (int phiUse : phiBinsN) { + for (const auto& ptUse : ptBinsN) { + for (const auto& etaUse : etaBinsN) { + for (const auto& phiUse : phiBinsN) { if (ptUse == ptB && etaUse == etaB && phiUse == phiB) continue; @@ -2311,8 +3322,16 @@ struct lambdaspincorrderived { auto tX = V0s.iteratorAt(static_cast(bc.rowIndex)); - if (!selectionV0(tX)) + bool tXPassed = false; + for (int sysId = 0; sysId < nSysTotal; ++sysId) { + if (selectionV0Sys(tX, sysCuts[sysId])) { + tXPassed = true; + break; + } + } + if (!tXPassed) continue; + if (!checkKinematics(tRep, tX)) continue; @@ -2430,9 +3449,16 @@ struct lambdaspincorrderived { 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; + + std::vector commonSys; + for (int sysId = 0; sysId < nSysTotal; ++sysId) { + if (selectionV0Sys(t1, sysCuts[sysId]) && selectionV0Sys(t2, sysCuts[sysId])) { + commonSys.push_back(sysId); + } } + if (commonSys.empty()) + continue; + if (t2.index() <= t1.index()) { continue; } @@ -2441,9 +3467,7 @@ struct lambdaspincorrderived { const bool doMixLeg1 = (cfgMixLegMode.value == 0 || cfgMixLegMode.value == 2); const bool doMixLeg2 = (cfgMixLegMode.value == 1 || cfgMixLegMode.value == 2); - // Fill TGT maps before searching for replacements. This makes TGT the - // true same-event target phase space for the selected pair, matching the - // hPtYSame definition. REP below is filled only for accepted replacements. + // Fill TGT maps before searching for replacements. if (doMixLeg1) { fillReplacementControlMap(t1.v0Status(), t2.v0Status(), 1, true, ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), @@ -2477,15 +3501,22 @@ struct lambdaspincorrderived { matches2.clear(); } - // Do not fill TGT here. TGT has already been filled above as the - // selected same-event target phase space, independent of replacement success. int nFill1 = 0; int nFill2 = 0; + // count actual accepted fills for leg-1 replacement if (doMixLeg1) { for (auto const& m : matches1) { auto tX = V0s.iteratorAt(static_cast(m.rowIndex)); - if (!selectionV0(tX)) + + bool hasCommon = false; + for (const auto& sysId : commonSys) { + if (selectionV0Sys(tX, sysCuts[sysId])) { + hasCommon = true; + break; + } + } + if (!hasCommon) continue; if (tX.v0Status() != t1.v0Status()) continue; @@ -2503,11 +3534,17 @@ struct lambdaspincorrderived { } } - // count actual accepted fills for leg-2 replacement if (doMixLeg2) { for (auto const& m : matches2) { auto tY = V0s.iteratorAt(static_cast(m.rowIndex)); - if (!selectionV0(tY)) + bool hasCommon = false; + for (const auto& sysId : commonSys) { + if (selectionV0Sys(tY, sysCuts[sysId])) { + hasCommon = true; + break; + } + } + if (!hasCommon) continue; if (tY.v0Status() != t2.v0Status()) continue; @@ -2544,9 +3581,7 @@ struct lambdaspincorrderived { if (nFill1 <= 0 && nFill2 <= 0) { continue; } - // Residual-weight QA needs a leg-specific normalization: - // TGT_leg is filled once per same-event target candidate with weight 1. - // REP_leg is the average replacement distribution for that same target. + const float wSELeg1 = (nFill1 > 0) ? 1.0f / static_cast(nFill1) : 0.0f; const float wSELeg2 = (nFill2 > 0) ? 1.0f / static_cast(nFill2) : 0.0f; @@ -2559,16 +3594,11 @@ struct lambdaspincorrderived { const float finalMixWeightLeg1 = branchNorm * wSELeg1; const float finalMixWeightLeg2 = branchNorm * wSELeg2; - // Do not fill target maps here: TGT has already been filled above, - // before searching for replacements, so that TGT represents all selected - // same-event targets rather than only targets with successful replacements. - if (doMixLeg1 && nFill1 > 0) { for (auto const& m : matches1) { auto tX = V0s.iteratorAt(static_cast(m.rowIndex)); - if (!selectionV0(tX)) { + if (!selectionV0(tX)) continue; - } if (tX.v0Status() != t1.v0Status()) continue; if (!checkKinematics(t1, tX)) @@ -2581,35 +3611,38 @@ struct lambdaspincorrderived { continue; if (hasSharedDaughters(tX, t2)) continue; + fillReplacementControlMap(tX.v0Status(), t2.v0Status(), 1, false, ROOT::Math::PtEtaPhiMVector(tX.lambdaPt(), tX.lambdaEta(), tX.lambdaPhi(), tX.lambdaMass()), wSELeg1); - // Fixed leg for leg1 replacement is original t2. - // Fill only for successful replacement. fillFixedLegControlMap(tX.v0Status(), t2.v0Status(), 1, false, ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()), wSELeg1); + 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()); - const float meWeight = finalMixWeightLeg1; const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); if ((tX.v0Status() == 0 && t2.v0Status() == 1) || (tX.v0Status() == 1 && t2.v0Status() == 0)) if (fillBasicQAHistos) histos.fill(HIST("deltaPhiMix"), dPhi, meWeight); + const int s1 = tX.v0Status(); const int s2 = t2.v0Status(); - if (s1 == 0 && s2 == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, 1, 1); - } else if (s1 == 1 && s2 == 0) { - fillHistograms(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, 2, 1); - } else { - fillHistograms(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, 1, 1); + for (const auto& sysId : commonSys) { + if (selectionV0Sys(tX, sysCuts[sysId])) { + if (s1 == 0 && s2 == 1) { + fillHistogramsSys(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, sysId, 1, 1); + } else if (s1 == 1 && s2 == 0) { + fillHistogramsSys(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, sysId, 2, 1); + } else { + fillHistogramsSys(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, sysId, 1, 1); + } + } } } } @@ -2617,9 +3650,8 @@ struct lambdaspincorrderived { if (doMixLeg2 && nFill2 > 0) { for (auto const& m : matches2) { auto tY = V0s.iteratorAt(static_cast(m.rowIndex)); - if (!selectionV0(tY)) { + if (!selectionV0(tY)) continue; - } if (tY.v0Status() != t2.v0Status()) continue; if (!checkKinematics(t2, tY)) @@ -2632,32 +3664,37 @@ struct lambdaspincorrderived { continue; if (hasSharedDaughters(tY, t2)) continue; + fillReplacementControlMap(t1.v0Status(), tY.v0Status(), 2, false, ROOT::Math::PtEtaPhiMVector(tY.lambdaPt(), tY.lambdaEta(), tY.lambdaPhi(), tY.lambdaMass()), wSELeg2); - // Fixed leg for leg2 replacement is original t1. - // Fill only for successful replacement. fillFixedLegControlMap(t1.v0Status(), tY.v0Status(), 2, false, ROOT::Math::PtEtaPhiMVector(t1.lambdaPt(), t1.lambdaEta(), t1.lambdaPhi(), t1.lambdaMass()), wSELeg2); + 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()); const float meWeight = finalMixWeightLeg2; const float dPhi = deltaPhiMinusPiToPi((float)lambda.Phi(), (float)lambda2.Phi()); if (fillBasicQAHistos) histos.fill(HIST("deltaPhiMix"), dPhi, meWeight); + const int s1 = t1.v0Status(); const int s2 = tY.v0Status(); - if (s1 == 0 && s2 == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, 2, 2); - } else if (s1 == 1 && s2 == 0) { - fillHistograms(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, 1, 2); - } else { - fillHistograms(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, 2, 2); + + for (const auto& sysId : commonSys) { + if (selectionV0Sys(tY, sysCuts[sysId])) { + if (s1 == 0 && s2 == 1) { + fillHistogramsSys(0, 1, lambda, lambda2, proton, proton2, 1, meWeight, sysId, 2, 2); + } else if (s1 == 1 && s2 == 0) { + fillHistogramsSys(0, 1, lambda2, lambda, proton2, proton, 1, meWeight, sysId, 1, 2); + } else { + fillHistogramsSys(s1, s2, lambda, lambda2, proton, proton2, 1, meWeight, sysId, 2, 2); + } + } } } } @@ -2668,7 +3705,7 @@ struct lambdaspincorrderived { pendingAtStart, pendingMatched, pendingExpired, pendingAdded, v6Pending.data.size()); } } - PROCESS_SWITCH(lambdaspincorrderived, processMEV6, "Process data ME v6 with radius buffer", false); + PROCESS_SWITCH(lambdaspincorrderived, processMEV6Sys, "Process data ME v6 Sys analysis", false); void processMCMEV6(EventCandidatesMC const& collisions, AllTrackCandidatesMC const& V0sMC) { @@ -2780,9 +3817,9 @@ struct lambdaspincorrderived { auto collectFromBins = [&](const std::vector& ptUseBins, const std::vector& etaUseBins, const std::vector& phiUseBins) { - for (int ptUse : ptUseBins) { - for (int etaUse : etaUseBins) { - for (int phiUse : phiUseBins) { + for (const auto& ptUse : ptUseBins) { + for (const auto& etaUse : etaUseBins) { + for (const auto& phiUse : phiUseBins) { const auto& vec = buffer[linearKeyR(colBin, status, ptUse, etaUse, phiUse, mB, rB, nStat, nPt, nEta, nPhi, nM, nR)]; @@ -2835,9 +3872,9 @@ struct lambdaspincorrderived { collectNeighborBinsClamp(etaB, nEta, nN_eta, etaBinsN); collectNeighborBinsPhi(phiB, nPhi, nN_phi, phiBinsN); - for (int ptUse : ptBinsN) { - for (int etaUse : etaBinsN) { - for (int phiUse : phiBinsN) { + for (const auto& ptUse : ptBinsN) { + for (const auto& etaUse : etaBinsN) { + for (const auto& phiUse : phiBinsN) { if (ptUse == ptB && etaUse == etaB && phiUse == phiB) { continue; }