From cc7c4eeebfe69fc132f0085c35b83918fbfa8f92 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Sun, 9 Nov 2025 02:03:23 +0100 Subject: [PATCH 1/3] fix bug of histogram fill for pid --- PWGLF/Tasks/Resonances/doublephimeson.cxx | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index a95af71b955..eae7465d59c 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -752,16 +752,25 @@ struct doublephimeson { // --- φ multiplicity with PID --- int phimult = 0; for (auto const& t : phitracks) { - if (t.phiMass() < minPhiMass1 || t.phiMass() > maxPhiMass1) - continue; const double kpluspt = std::hypot(t.phid1Px(), t.phid1Py()); const double kminuspt = std::hypot(t.phid2Px(), t.phid2Py()); + // PID QA before + histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), t.phid1TPC(), t.phid1TOF(), kplus1pt); + histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), t.phid1TPC(), kplus1pt); + histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), t.phid2TPC(), kminus1pt); + if (t.phiMass() < minPhiMass1 || t.phiMass() > maxPhiMass1) + continue; if (kpluspt > maxKaonPt || kminuspt > maxKaonPt) continue; if (!selectionPID(t.phid1TPC(), t.phid1TOF(), t.phid1TOFHit(), strategyPID1, kpluspt)) continue; if (!selectionPID(t.phid2TPC(), t.phid2TOF(), t.phid2TOFHit(), strategyPID2, kminuspt)) continue; + // PID QA after + histos.fill(HIST("hnsigmaTPCTOFKaon"), t1.phid1TPC(), t1.phid1TOF(), kplus1pt); + histos.fill(HIST("hnsigmaTPCKaonPlus"), t1.phid1TPC(), kplus1pt); + histos.fill(HIST("hnsigmaTPCKaonMinus"), t1.phid2TPC(), kminus1pt); + histos.fill(HIST("hPhiMass"), phi1.M(), phi1.Pt()); ++phimult; } if (phimult < 2) @@ -815,10 +824,7 @@ struct doublephimeson { const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py()); const double kminus1pt = std::hypot(t1.phid2Px(), t1.phid2Py()); - // PID QA before - histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), t1.phid1TPC(), t1.phid1TOF(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), t1.phid1TPC(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), t1.phid2TPC(), kminus1pt); + if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt) continue; @@ -837,11 +843,6 @@ struct doublephimeson { continue; if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt) continue; - // PID QA after - histos.fill(HIST("hnsigmaTPCTOFKaon"), t1.phid1TPC(), t1.phid1TOF(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonPlus"), t1.phid1TPC(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonMinus"), t1.phid2TPC(), kminus1pt); - histos.fill(HIST("hPhiMass"), phi1.M(), phi1.Pt()); const auto id1 = t1.index(); From d80793174970a767bb3035e1253783cd2993b221 Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 10 Nov 2025 09:13:33 +0100 Subject: [PATCH 2/3] Add PID histogram --- PWGLF/Tasks/Resonances/doublephimeson.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index eae7465d59c..869c92133f8 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -755,9 +755,9 @@ struct doublephimeson { const double kpluspt = std::hypot(t.phid1Px(), t.phid1Py()); const double kminuspt = std::hypot(t.phid2Px(), t.phid2Py()); // PID QA before - histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), t.phid1TPC(), t.phid1TOF(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), t.phid1TPC(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), t.phid2TPC(), kminus1pt); + histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), t.phid1TPC(), t.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), t.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), t.phid2TPC(), kminuspt); if (t.phiMass() < minPhiMass1 || t.phiMass() > maxPhiMass1) continue; if (kpluspt > maxKaonPt || kminuspt > maxKaonPt) @@ -767,10 +767,10 @@ struct doublephimeson { if (!selectionPID(t.phid2TPC(), t.phid2TOF(), t.phid2TOFHit(), strategyPID2, kminuspt)) continue; // PID QA after - histos.fill(HIST("hnsigmaTPCTOFKaon"), t1.phid1TPC(), t1.phid1TOF(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonPlus"), t1.phid1TPC(), kplus1pt); - histos.fill(HIST("hnsigmaTPCKaonMinus"), t1.phid2TPC(), kminus1pt); - histos.fill(HIST("hPhiMass"), phi1.M(), phi1.Pt()); + histos.fill(HIST("hnsigmaTPCTOFKaon"), t.phid1TPC(), t.phid1TOF(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonPlus"), t.phid1TPC(), kpluspt); + histos.fill(HIST("hnsigmaTPCKaonMinus"), t.phid2TPC(), kminuspt); + ++phimult; } if (phimult < 2) @@ -845,7 +845,7 @@ struct doublephimeson { continue; const auto id1 = t1.index(); - + histos.fill(HIST("hPhiMass"), phi1.M(), phi1.Pt()); for (auto const& t2 : phitracks) { const auto id2 = t2.index(); if (id2 <= id1) From ed1671f6647e352d6ba2f67d1cb04af3ac5b9e2a Mon Sep 17 00:00:00 2001 From: skundu692 Date: Mon, 10 Nov 2025 19:18:55 +0100 Subject: [PATCH 3/3] Improve both leg bixing --- PWGLF/Tasks/Resonances/doublephimeson.cxx | 2 - .../Strangeness/lambdaspincorrderived.cxx | 458 ++++++++++++++---- 2 files changed, 366 insertions(+), 94 deletions(-) diff --git a/PWGLF/Tasks/Resonances/doublephimeson.cxx b/PWGLF/Tasks/Resonances/doublephimeson.cxx index 869c92133f8..27409836425 100644 --- a/PWGLF/Tasks/Resonances/doublephimeson.cxx +++ b/PWGLF/Tasks/Resonances/doublephimeson.cxx @@ -824,8 +824,6 @@ struct doublephimeson { const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py()); const double kminus1pt = std::hypot(t1.phid2Px(), t1.phid2Py()); - - if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt) continue; if (!selectionPID(t1.phid1TPC(), t1.phid1TOF(), t1.phid1TOFHit(), strategyPID1, kplus1pt)) diff --git a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx index caacb8c0ed6..8accd5f7bd0 100644 --- a/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx +++ b/PWGLF/Tasks/Strangeness/lambdaspincorrderived.cxx @@ -298,6 +298,239 @@ struct lambdaspincorrderived { const ROOT::Math::PtEtaPhiMVector& daughpart1, const ROOT::Math::PtEtaPhiMVector& daughpart2, int datatype, float mixpairweight) { + // --- Mass handling --- + auto lambda1Mass = 0.0; + auto lambda2Mass = 0.0; + if (!usePDGM) { + lambda1Mass = particle1.M(); + lambda2Mass = particle2.M(); + } else { + lambda1Mass = o2::constants::physics::MassLambda; + lambda2Mass = o2::constants::physics::MassLambda; + } + 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()}; // boosting vector for pair CM + + // Step1: Boosting both Lambdas to Lambda-Lambda pair rest frame + auto lambda1CM = boostPairToCM(particle1Dummy); + auto lambda2CM = boostPairToCM(particle2Dummy); + + // Step 2: Boost Each Lambda to its Own Rest Frame + ROOT::Math::Boost boostLambda1ToCM{lambda1CM.BoostToCM()}; + ROOT::Math::Boost boostLambda2ToCM{lambda2CM.BoostToCM()}; + + // Also boost the daughter protons to the same frame + auto proton1pairCM = boostPairToCM(daughpart1); // proton1 to pair CM + auto proton2pairCM = boostPairToCM(daughpart2); // proton2 to pair CM + + // Boost protons into their respective Lambda rest frames + auto proton1LambdaRF = boostLambda1ToCM(proton1pairCM); + auto proton2LambdaRF = boostLambda2ToCM(proton2pairCM); + + // --- STAR-style Δθ (as written: dot product of proton directions in their own Λ RFs) --- + ROOT::Math::Boost boostL1_LabToRF{particle1Dummy.BoostToCM()}; // Λ1 velocity in lab + ROOT::Math::Boost boostL2_LabToRF{particle2Dummy.BoostToCM()}; // Λ2 velocity in lab + + auto p1_LRF = boostL1_LabToRF(daughpart1); + auto p2_LRF = boostL2_LabToRF(daughpart2); + + 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(); + + TVector3 k1(proton1LambdaRF.Px(), proton1LambdaRF.Py(), proton1LambdaRF.Pz()); + k1 = k1.Unit(); + TVector3 k2(proton2LambdaRF.Px(), proton2LambdaRF.Py(), proton2LambdaRF.Pz()); + k2 = k2.Unit(); + + 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; + + double cosDeltaTheta_hel = k1.Dot(k2); + if (cosDeltaTheta_hel > 1.0) + cosDeltaTheta_hel = 111.0; + if (cosDeltaTheta_hel < -1.0) + cosDeltaTheta_hel = -111.0; + + auto cosThetaDiff = -999.0; + if (cosDef == 0) { + cosThetaDiff = cosDeltaTheta_STAR_naive; + } else { + cosThetaDiff = cosDeltaTheta_hel; + } + + // --- kinematics for weights and ΔR etc. --- + double pt1 = particle1.Pt(); + double dphi1 = RecoDecay::constrainAngle(particle1.Phi(), 0.0F, harmonic); + double deta1 = particle1.Eta(); + + double pt2 = particle2.Pt(); + double dphi2 = RecoDecay::constrainAngle(particle2.Phi(), 0.0F, harmonic); + double deta2 = particle2.Eta(); + + double deta_pair = std::abs(deta1 - deta2); + double dphi_pair = RecoDecay::constrainAngle(particle1.Phi() - particle2.Phi(), 0.0F, harmonicDphi); + + double deltaR = TMath::Sqrt(deta_pair * deta_pair + dphi_pair * dphi_pair); + double deltaRap = std::abs(particle1.Rapidity() - particle2.Rapidity()); + + // ------------------------------------------------------------------ + // Per-leg, per-species efficiency weights (only for mixed events) + // ------------------------------------------------------------------ + double epsWeight1 = 1.0; + double epsWeight2 = 1.0; + + if (useweight && datatype == 1) { + // Leg 1: tag1 == 0 -> Lambda, tag1 == 1 -> AntiLambda + if (tag1 == 0) { + if (hweight1) { + epsWeight1 = hweight1->GetBinContent(hweight1->FindBin(dphi1, deta1, pt1)); + } + } else { + if (hweight4) { + epsWeight1 = hweight4->GetBinContent(hweight4->FindBin(dphi1, deta1, pt1)); + } + } + + // Leg 2: tag2 == 0 -> Lambda, tag2 == 1 -> AntiLambda + if (tag2 == 0) { + if (hweight12) { + epsWeight2 = hweight12->GetBinContent(hweight12->FindBin(dphi2, deta2, pt2)); + } + } else { + if (hweight42) { + epsWeight2 = hweight42->GetBinContent(hweight42->FindBin(dphi2, deta2, pt2)); + } + } + + // safety against zero or NaN + if (epsWeight1 <= 0.0 || !std::isfinite(epsWeight1)) { + epsWeight1 = 1.0; + } + if (epsWeight2 <= 0.0 || !std::isfinite(epsWeight2)) { + epsWeight2 = 1.0; + } + } + + // ====================================== + // SAME-EVENT (datatype == 0) + // ====================================== + if (datatype == 0) { + mixpairweight = 1.0; + histos.fill(HIST("hPtYSame"), particle1.Pt(), particle1.Rapidity(), mixpairweight); + + if (tag1 == 0 && tag2 == 0) { + histos.fill(HIST("SE_LL"), dphi1, deta1, pt1, mixpairweight); + histos.fill(HIST("SE_LL2"), dphi2, deta2, pt2, mixpairweight); + histos.fill(HIST("hSparseLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, mixpairweight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, mixpairweight); + histos.fill(HIST("hSparsePhiLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, mixpairweight); + histos.fill(HIST("hSparsePairMassLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), mixpairweight); + } + } else if (tag1 == 0 && tag2 == 1) { + histos.fill(HIST("SE_LAL"), dphi1, deta1, pt1, mixpairweight); + histos.fill(HIST("SE_LAL2"), dphi2, deta2, pt2, mixpairweight); + histos.fill(HIST("hSparseLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, mixpairweight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, mixpairweight); + histos.fill(HIST("hSparsePhiLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, mixpairweight); + histos.fill(HIST("hSparsePairMassLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), mixpairweight); + } + } else if (tag1 == 1 && tag2 == 0) { + histos.fill(HIST("hSparseAntiLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, mixpairweight); + histos.fill(HIST("SE_ALL"), dphi1, deta1, pt1, mixpairweight); + histos.fill(HIST("SE_ALL2"), dphi2, deta2, pt2, mixpairweight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapAntiLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, mixpairweight); + histos.fill(HIST("hSparsePhiAntiLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, mixpairweight); + histos.fill(HIST("hSparsePairMassAntiLambdaLambda"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), mixpairweight); + } + } else if (tag1 == 1 && tag2 == 1) { + histos.fill(HIST("hSparseAntiLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaR, mixpairweight); + histos.fill(HIST("SE_ALAL"), dphi1, deta1, pt1, mixpairweight); + histos.fill(HIST("SE_ALAL2"), dphi2, deta2, pt2, mixpairweight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapAntiLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, mixpairweight); + histos.fill(HIST("hSparsePhiAntiLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, mixpairweight); + histos.fill(HIST("hSparsePairMassAntiLambdaAntiLambda"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), mixpairweight); + } + } + } + // ====================================== + // MIXED-EVENT (datatype == 1) + // ====================================== + else if (datatype == 1) { + double weight = mixpairweight; + if (useweight) { + if (usebothweight) { + // symmetric pair weight: 1 / (ε1 * ε2) + weight = mixpairweight / (epsWeight1 * epsWeight2); + } else { + // optional single-leg mode + weight = mixpairweight / epsWeight1; + } + } + if (weight <= 0.0 || !std::isfinite(weight)) { + weight = 1.0; + } + + histos.fill(HIST("hPtYMix"), particle1.Pt(), particle1.Rapidity(), weight); + + if (tag1 == 0 && tag2 == 0) { + histos.fill(HIST("ME_LL"), dphi1, deta1, pt1, weight); + histos.fill(HIST("ME_LL2"), dphi2, deta2, pt2, weight); + histos.fill(HIST("hSparseLambdaLambdaMixed"), + particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight); + histos.fill(HIST("hSparsePhiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, weight); + histos.fill(HIST("hSparsePairMassLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), weight); + } + } else if (tag1 == 0 && tag2 == 1) { + histos.fill(HIST("ME_LAL"), dphi1, deta1, pt1, weight); + histos.fill(HIST("ME_LAL2"), dphi2, deta2, pt2, weight); + histos.fill(HIST("hSparseLambdaAntiLambdaMixed"), + particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight); + histos.fill(HIST("hSparsePhiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, weight); + histos.fill(HIST("hSparsePairMassLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), weight); + } + } else if (tag1 == 1 && tag2 == 0) { + histos.fill(HIST("ME_ALL"), dphi1, deta1, pt1, weight); + histos.fill(HIST("ME_ALL2"), dphi2, deta2, pt2, weight); + histos.fill(HIST("hSparseAntiLambdaLambdaMixed"), + particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight); + histos.fill(HIST("hSparsePhiAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, weight); + histos.fill(HIST("hSparsePairMassAntiLambdaLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), weight); + } + } else if (tag1 == 1 && tag2 == 1) { + histos.fill(HIST("ME_ALAL"), dphi1, deta1, pt1, weight); + histos.fill(HIST("ME_ALAL2"), dphi2, deta2, pt2, weight); + histos.fill(HIST("hSparseAntiLambdaAntiLambdaMixed"), + particle1.M(), particle2.M(), cosThetaDiff, deltaR, weight); + if (useAdditionalHisto) { + histos.fill(HIST("hSparseRapAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, deltaRap, weight); + histos.fill(HIST("hSparsePhiAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, dphi_pair, weight); + histos.fill(HIST("hSparsePairMassAntiLambdaAntiLambdaMixed"), particle1.M(), particle2.M(), cosThetaDiff, pairDummy.M(), weight); + } + } + } + } + + void fillHistograms2(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) + { auto lambda1Mass = 0.0; auto lambda2Mass = 0.0; @@ -386,18 +619,34 @@ struct lambdaspincorrderived { double epsWeight2 = 1.0; if (useweight && datatype == 1) { - if (tag1 == 0 && tag2 == 0) { - epsWeight1 = hweight1->GetBinContent(hweight1->FindBin(dphi1, deta1, pt1)); - epsWeight2 = hweight12->GetBinContent(hweight12->FindBin(dphi2, deta2, pt2)); - } else if (tag1 == 0 && tag2 == 1) { - epsWeight1 = hweight2->GetBinContent(hweight2->FindBin(dphi1, deta1, pt1)); - epsWeight2 = hweight22->GetBinContent(hweight22->FindBin(dphi2, deta2, pt2)); - } else if (tag1 == 1 && tag2 == 0) { - epsWeight1 = hweight3->GetBinContent(hweight3->FindBin(dphi1, deta1, pt1)); - epsWeight2 = hweight32->GetBinContent(hweight32->FindBin(dphi2, deta2, pt2)); - } else if (tag1 == 1 && tag2 == 1) { - epsWeight1 = hweight4->GetBinContent(hweight4->FindBin(dphi1, deta1, pt1)); - epsWeight2 = hweight42->GetBinContent(hweight42->FindBin(dphi2, deta2, pt2)); + // --- Leg 1: choose map only by species (tag1) --- + if (tag1 == 0) { // Λ + if (hweight1) { + epsWeight1 = hweight1->GetBinContent(hweight1->FindBin(dphi1, deta1, pt1)); + } + } else { // Λ̄ + if (hweight4) { + epsWeight1 = hweight4->GetBinContent(hweight4->FindBin(dphi1, deta1, pt1)); + } + } + + // --- Leg 2: choose map only by species (tag2) --- + if (tag2 == 0) { // Λ + if (hweight12) { + epsWeight2 = hweight12->GetBinContent(hweight12->FindBin(dphi2, deta2, pt2)); + } + } else { // Λ̄ + if (hweight42) { + epsWeight2 = hweight42->GetBinContent(hweight42->FindBin(dphi2, deta2, pt2)); + } + } + + // Safety: avoid zero/NaN + if (epsWeight1 <= 0.0 || !std::isfinite(epsWeight1)) { + epsWeight1 = 1.0; + } + if (epsWeight2 <= 0.0 || !std::isfinite(epsWeight2)) { + epsWeight2 = 1.0; } } @@ -538,16 +787,16 @@ struct lambdaspincorrderived { lambda2 = ROOT::Math::PtEtaPhiMVector(v02.lambdaPt(), v02.lambdaEta(), v02.lambdaPhi(), v02.lambdaMass()); histos.fill(HIST("deltaPhiSame"), RecoDecay::constrainAngle(v0.lambdaPhi() - v02.lambdaPhi(), 0.0F, harmonicDphi)); if (v0.v0Status() == 0 && v02.v0Status() == 0) { - fillHistograms(0, 0, lambda, lambda2, proton, proton2, 0, 1.0); + fillHistograms2(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); + fillHistograms2(0, 1, lambda, lambda2, proton, proton2, 0, 1.0); } if (v0.v0Status() == 1 && v02.v0Status() == 0) { - fillHistograms(1, 0, lambda, lambda2, proton, proton2, 0, 1.0); + fillHistograms2(1, 0, lambda, lambda2, proton, proton2, 0, 1.0); } if (v0.v0Status() == 1 && v02.v0Status() == 1) { - fillHistograms(1, 1, lambda, lambda2, proton, proton2, 0, 1.0); + fillHistograms2(1, 1, lambda, lambda2, proton, proton2, 0, 1.0); } } } @@ -610,16 +859,16 @@ struct lambdaspincorrderived { lambda2 = ROOT::Math::PtEtaPhiMVector(t2.lambdaPt(), t2.lambdaEta(), t2.lambdaPhi(), t2.lambdaMass()); histos.fill(HIST("deltaPhiMix"), RecoDecay::constrainAngle(t3.lambdaPhi() - t2.lambdaPhi(), 0.0F, harmonicDphi)); if (t3.v0Status() == 0 && t2.v0Status() == 0) { - fillHistograms(0, 0, lambda, lambda2, proton, proton2, 1, 1.0); + fillHistograms2(0, 0, lambda, lambda2, proton, proton2, 1, 1.0); } if (t3.v0Status() == 0 && t2.v0Status() == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, 1.0); + fillHistograms2(0, 1, lambda, lambda2, proton, proton2, 1, 1.0); } if (t3.v0Status() == 1 && t2.v0Status() == 0) { - fillHistograms(1, 0, lambda, lambda2, proton, proton2, 1, 1.0); + fillHistograms2(1, 0, lambda, lambda2, proton, proton2, 1, 1.0); } if (t3.v0Status() == 1 && t2.v0Status() == 1) { - fillHistograms(1, 1, lambda, lambda2, proton, proton2, 1, 1.0); + fillHistograms2(1, 1, lambda, lambda2, proton, proton2, 1, 1.0); } } } // replacement track pair @@ -690,16 +939,16 @@ struct lambdaspincorrderived { histos.fill(HIST("deltaPhiMix"), dPhi, invN); if (t3.v0Status() == 0 && t2.v0Status() == 0) { - fillHistograms(0, 0, lambda, lambda2, proton, proton2, 1, invN); + fillHistograms2(0, 0, lambda, lambda2, proton, proton2, 1, invN); } if (t3.v0Status() == 0 && t2.v0Status() == 1) { - fillHistograms(0, 1, lambda, lambda2, proton, proton2, 1, invN); + fillHistograms2(0, 1, lambda, lambda2, proton, proton2, 1, invN); } if (t3.v0Status() == 1 && t2.v0Status() == 0) { - fillHistograms(1, 0, lambda, lambda2, proton, proton2, 1, invN); + fillHistograms2(1, 0, lambda, lambda2, proton, proton2, 1, invN); } if (t3.v0Status() == 1 && t2.v0Status() == 1) { - fillHistograms(1, 1, lambda, lambda2, proton, proton2, 1, invN); + fillHistograms2(1, 1, lambda, lambda2, proton, proton2, 1, invN); } } } // end mixing-event loop @@ -798,7 +1047,7 @@ struct lambdaspincorrderived { const float dPhi = std::fabs(RecoDecay::constrainAngle(lambda.Phi(), 0.0F, harmonic) - RecoDecay::constrainAngle(lambda2.Phi(), 0.0F, harmonic)); histos.fill(HIST("deltaPhiMix"), dPhi, wBase); - fillHistograms(tX.v0Status(), t2.v0Status(), lambda, lambda2, proton, proton2, 1, wBase); + fillHistograms2(tX.v0Status(), t2.v0Status(), lambda, lambda2, proton, proton2, 1, wBase); } } } @@ -884,8 +1133,6 @@ struct lambdaspincorrderived { { return ((((((static_cast(colBin) * nStatus + statBin) * nPt + ptBin) * nEta + etaBin) * nPhi + phiBin) * nM + mBin)); } - - // ===================== Main mixing (with mass-bin + random unique sampling) ===================== void processMEV4(EventCandidates const& collisions, AllTrackCandidates const& V0s) { // Build binner from your existing configurables @@ -895,8 +1142,8 @@ struct lambdaspincorrderived { phiMix.value // φ step; φ range fixed to [0, 2π) }; - const int nCol = colBinning.getAllBinsCount(); // event-class bins (vz, centrality) - const int nStat = N_STATUS; // 2 + 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(); @@ -905,7 +1152,9 @@ struct lambdaspincorrderived { const size_t nKeys = static_cast(nCol) * nStat * nPt * nEta * nPhi * nM; std::vector> buffer(nKeys); - // ---- PASS 1: fill 6D buffer ---- + // ========================= + // PASS 1: fill 6D buffer + // ========================= for (auto const& col : collisions) { const int colBin = colBinning.getBin(std::make_tuple(col.posz(), col.cent())); auto slice = V0s.sliceBy(tracksPerCollisionV0, col.index()); @@ -918,7 +1167,6 @@ struct lambdaspincorrderived { if (status < 0 || status >= nStat) continue; - // Bin kinematics (φ already constrained via your call-site) const int ptB = mb.ptBin(t.lambdaPt()); const int etaB = mb.etaBin(t.lambdaEta()); const int phiB = mb.phiBin(RecoDecay::constrainAngle(t.lambdaPhi(), 0.0F, harmonic)); @@ -931,7 +1179,7 @@ struct lambdaspincorrderived { buffer[key].push_back(BufferCand{ .collisionIdx = static_cast(col.index()), - .rowIndex = static_cast(t.globalIndex()), // adapt accessor if needed + .rowIndex = static_cast(t.globalIndex()), .v0Status = static_cast(status), .ptBin = static_cast(ptB), .etaBin = static_cast(etaB), @@ -940,10 +1188,54 @@ struct lambdaspincorrderived { } } - // ---- PASS 2: mixing over same-event pairs ---- + // Helper: get matches for a given candidate (for mixing one leg) + auto makeMatchesFor = [&](auto const& cand, + int colBinLocal, + int64_t curColIdx) -> std::vector { + std::vector matches; + + const int status = static_cast(cand.v0Status()); + if (status < 0 || status >= nStat) + return matches; + + const int ptB = mb.ptBin(cand.lambdaPt()); + const int etaB = mb.etaBin(cand.lambdaEta()); + const int phiB = mb.phiBin(RecoDecay::constrainAngle(cand.lambdaPhi(), 0.0F, harmonic)); + const int mB = mb.massBin(cand.lambdaMass()); + if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) + return matches; + + const size_t key = linearKey(colBinLocal, status, ptB, etaB, phiB, mB, + nStat, nPt, nEta, nPhi, nM); + auto const& binVec = buffer[key]; + if (binVec.empty()) + return matches; + + matches.reserve(binVec.size()); + for (const auto& bc : binVec) { + if (bc.collisionIdx == curColIdx) + continue; // ensure different event + matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); + } + + // Optionally cap number of matches by nEvtMixing + const int cap = nEvtMixing.value; + const int n = static_cast(matches.size()); + if (cap > 0 && n > cap) { + std::shuffle(matches.begin(), matches.end(), rng); + matches.resize(cap); + } + + return matches; + }; + + // ========================= + // PASS 2: mixing over SE pairs (mix both legs) + // ========================= for (auto const& collision1 : collisions) { const int colBin = colBinning.getBin(std::make_tuple(collision1.posz(), collision1.cent())); auto poolA = V0s.sliceBy(tracksPerCollisionV0, collision1.index()); + const int64_t curColId = static_cast(collision1.index()); for (auto const& [t1, t2] : soa::combinations(o2::soa::CombinationsFullIndexPolicy(poolA, poolA))) { @@ -963,80 +1255,62 @@ struct lambdaspincorrderived { if (t1.pionIndex() == t2.protonIndex()) continue; - const int status = static_cast(t1.v0Status()); - if (status < 0 || status >= nStat) - continue; - - // Bin of t1 defines where to search (exact 6D bin) - const int ptB = mb.ptBin(t1.lambdaPt()); - const int etaB = mb.etaBin(t1.lambdaEta()); - const int phiB = mb.phiBin(RecoDecay::constrainAngle(t1.lambdaPhi(), 0.0F, harmonic)); // φ already constrained upstream - const int mB = mb.massBin(t1.lambdaMass()); - if (ptB < 0 || etaB < 0 || phiB < 0 || mB < 0) - continue; - - const size_t key = linearKey(colBin, status, ptB, etaB, phiB, mB, - nStat, nPt, nEta, nPhi, nM); - auto const& binVec = buffer[key]; - if (binVec.empty()) - continue; - - // Collect all partners from this 6D bin but different collision - std::vector matches; - matches.reserve(binVec.size()); - const int64_t curColIdx = static_cast(collision1.index()); + // matches for replacing leg 1 and leg 2 + auto matches1 = makeMatchesFor(t1, colBin, curColId); // t1 -> tX + auto matches2 = makeMatchesFor(t2, colBin, curColId); // t2 -> tY - for (const auto& bc : binVec) { - if (bc.collisionIdx == curColIdx) - continue; // ensure different event - matches.push_back(MatchRef{bc.collisionIdx, bc.rowIndex}); - } - if (matches.empty()) + const int n1 = static_cast(matches1.size()); + const int n2 = static_cast(matches2.size()); + if (n1 == 0 && n2 == 0) continue; - // ---------- YOUR PREFERRED RANDOM UNIQUE SAMPLING BLOCK ---------- - const int cap = nEvtMixing.value; - const int n = static_cast(matches.size()); - if (cap > 0 && cap < n) { - std::uniform_int_distribution dist(0, n - 1); - // pick cap unique indices - std::unordered_set chosen; - chosen.reserve(cap * 2); - while ((int)chosen.size() < cap) { - chosen.insert(dist(rng)); - } - std::vector subset; - subset.reserve(cap); - for (int idx : chosen) - subset.push_back(matches[idx]); - matches.swap(subset); - } else { - std::shuffle(matches.begin(), matches.end(), rng); - } - // ---------------------------------------------------------------- - - const float wBase = 1.0f / static_cast(matches.size()); + // Split "unit" weight between the two modes + const float w1 = (n1 > 0 ? 0.5f / static_cast(n1) : 0.f); + const float w2 = (n2 > 0 ? 0.5f / static_cast(n2) : 0.f); - // Emit mixed pairs: tX replaces t1; keep t2 - for (const auto& m : matches) { - auto tX = V0s.iteratorAt(m.rowIndex); // replace accessor if different + // --- Type A: replace leg 1 → (tX, t2) --- + for (const auto& m : matches1) { + auto tX = V0s.iteratorAt(m.rowIndex); if (!selectionV0(tX)) - continue; // optional extra guard + continue; if (!checkKinematics(t1, tX)) continue; - 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 proton1 = ROOT::Math::PtEtaPhiMVector(tX.protonPt(), tX.protonEta(), tX.protonPhi(), o2::constants::physics::MassProton); + auto lambda1 = 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 = std::fabs(RecoDecay::constrainAngle(lambda.Phi() - lambda2.Phi(), 0.0F, harmonicDphi)); - histos.fill(HIST("deltaPhiMix"), dPhi, wBase); - fillHistograms(tX.v0Status(), t2.v0Status(), lambda, lambda2, proton, proton2, 1, wBase); + const float dPhi = std::fabs(RecoDecay::constrainAngle(lambda1.Phi() - lambda2.Phi(), 0.0F, harmonicDphi)); + histos.fill(HIST("deltaPhiMix"), dPhi, w1); + + fillHistograms2(tX.v0Status(), t2.v0Status(), lambda1, lambda2, proton1, proton2, 1, w1); + } + + // --- Type B: replace leg 2 → (t1, tY) --- + for (const auto& m : matches2) { + auto tY = V0s.iteratorAt(m.rowIndex); + if (!selectionV0(tY)) + continue; + if (!checkKinematics(t2, tY)) + continue; + + auto proton1 = ROOT::Math::PtEtaPhiMVector(t1.protonPt(), t1.protonEta(), t1.protonPhi(), o2::constants::physics::MassProton); + auto lambda1 = 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 float dPhi = std::fabs(RecoDecay::constrainAngle(lambda1.Phi() - lambda2.Phi(), 0.0F, harmonicDphi)); + histos.fill(HIST("deltaPhiMix"), dPhi, w2); + + fillHistograms2(t1.v0Status(), tY.v0Status(), lambda1, lambda2, proton1, proton2, 1, w2); } } } } + PROCESS_SWITCH(lambdaspincorrderived, processMEV4, "Process data ME (5d buffer)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)