From a14f25f8e95a510fb7216af3b635fc8b3acc110b Mon Sep 17 00:00:00 2001 From: Subhadeep Roy Date: Thu, 9 Apr 2026 23:53:56 +0530 Subject: [PATCH 1/4] Updated for detailed analysis --- .../Tasks/lambdaSpinPolarization.cxx | 604 ++++++++++++++---- 1 file changed, 480 insertions(+), 124 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx index e0eefb67fad..3ccfe9c3032 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx @@ -30,6 +30,8 @@ #include "Framework/runDataProcessing.h" #include +#include +#include #include #include @@ -119,11 +121,13 @@ namespace lambdatrackext DECLARE_SOA_COLUMN(LambdaSharingDaughter, lambdaSharingDaughter, bool); DECLARE_SOA_COLUMN(LambdaSharingDauIds, lambdaSharingDauIds, std::vector); DECLARE_SOA_COLUMN(TrueLambdaFlag, trueLambdaFlag, bool); +DECLARE_SOA_COLUMN(FakeFlag, fakeFlag, bool); } // namespace lambdatrackext DECLARE_SOA_TABLE(LambdaTracksExt, "AOD", "LAMBDATRACKSEXT", lambdatrackext::LambdaSharingDaughter, lambdatrackext::LambdaSharingDauIds, - lambdatrackext::TrueLambdaFlag); + lambdatrackext::TrueLambdaFlag, + lambdatrackext::FakeFlag); using LambdaTrackExt = LambdaTracksExt::iterator; @@ -239,13 +243,22 @@ enum RunType { enum ParticleType { kLambda = 0, - kAntiLambda + kAntiLambda, + kFakeLambda, + kFakeAntiLambda }; enum ParticlePairType { + // ----------------------------------------------------------------------- kLambdaAntiLambda = 0, - kLambdaLambda, - kAntiLambdaAntiLambda + kAntiLambdaLambda = 1, + kLambdaLambda = 2, + kAntiLambdaAntiLambda = 3, + // ----------------------------------------------------------------------- + kLambdaFakeAntiLambda = 4, + kAntiLambdaFakeLambda = 5, + kLambdaFakeLambda = 6, + kAntiLambdaFakeAntiLambda = 7 }; enum ShareDauLambda { @@ -315,6 +328,14 @@ struct LambdaTableProducer { Configurable cMaxChi2Tpc{"cMaxChi2Tpc", 4, "Max Chi2 Tpc"}; Configurable cTpcNsigmaCut{"cTpcNsigmaCut", 3.0, "TPC NSigma Selection Cut"}; Configurable cRemoveAmbiguousTracks{"cRemoveAmbiguousTracks", false, "Remove Ambiguous Tracks"}; + + // LS (fake) background control variables + Configurable cApplyFakeDcaCuts{"cApplyFakeDcaCuts", true, "Apply DCA-to-PV cuts on LS daughter tracks"}; + Configurable cFakeDaughterMaxDEta{"cFakeDaughterMaxDEta", -1.f, "Max delta eta between LS daughters as topology proxy (-1=off)"}; + Configurable cFakeDaughterMaxOpenAngle{"cFakeDaughterMaxOpenAngle", -1.f, "Max opening angle between LS daughters [rad] (-1=off)"}; + Configurable cAllowDualAssignment{"cAllowDualAssignment", true, "Allow both PID assignments per LS pair"}; + Configurable cMaxFakeLambdaPerEvent{"cMaxFakeLambdaPerEvent", -1, "Max fake Lambda/AntiLambda candidates per event (-1 = no limit)"}; + Configurable cRandomizeFakeSubsample{"cRandomizeFakeSubsample", false, "Randomly subsample fake candidates if over limit"}; // V0s Configurable cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.02, "Minimum Proton DCAr to PV"}; @@ -430,6 +451,9 @@ struct LambdaTableProducer { histos.add("Tracks/h1f_lambda_pt_vs_invm", "p_{T} vs M_{#Lambda}", kTH2F, {axisV0Mass, axisV0Pt}); histos.add("Tracks/h1f_antilambda_pt_vs_invm", "p_{T} vs M_{#bar{#Lambda}}", kTH2F, {axisV0Mass, axisV0Pt}); + //histogram for LS subsampling fraction (nKept/nTotal per event) + histos.add("Tracks/h1f_fake_subsample_fraction", "LS subsample fraction (nKept/nTotal)", kTH1F, {{110, 0., 1.1, "f_{subsample}"}}); + // QA Lambda histos.add("QA/Lambda/h2f_qt_vs_alpha", "Armentros-Podolanski Plot", kTH2F, {axisAlpha, axisQtarm}); histos.add("QA/Lambda/h1f_dca_V0_daughters", "DCA between V0 daughters", kTH1F, {axisDcaDau}); @@ -1016,6 +1040,45 @@ struct LambdaTableProducer { histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hPhi"), phi); } + // Basic track quality selection reused for fake-lambda (LS) daughter candidates + template + bool selTrackBasic(T const& track) + { + if (track.pt() <= cTrackMinPt || track.pt() >= cTrackMaxPt) + return false; + if (std::abs(track.eta()) >= cTrackEtaCut) + return false; + if (track.tpcNClsCrossedRows() <= cMinTpcCrossedRows) + return false; + if (track.tpcCrossedRowsOverFindableCls() < cMinTpcCROverCls) + return false; + if (track.tpcNClsShared() > cMaxTpcSharedClusters) + return false; + if (track.tpcChi2NCl() > cMaxChi2Tpc) + return false; + return true; + } + + template + bool selTrackBasicWithDCA(T const& track, bool isProton) + { + if (!selTrackBasic(track)) + return false; + if (cApplyFakeDcaCuts) { + float minDca = isProton ? (float)cMinDcaProtonToPV : (float)cMinDcaPionToPV; + if (std::abs(track.dcaXY()) < minDca) + return false; + } + return true; + } + + struct FakeCandidate { + float px, py, pz, pt, eta, phi, rap, mass; + float prPx, prPy, prPz; + int64_t id1, id2; + ParticleType type; + }; + // Reconstructed Level Tables template void fillLambdaRecoTables(C const& collision, B const& bc, V const& v0tracks, T const& tracks) @@ -1147,7 +1210,117 @@ struct LambdaTableProducer { v0.template posTrack_as().index(), v0.template negTrack_as().index(), v0.v0cosPA(), v0.dcaV0daughters(), (int8_t)v0Type, v0PrmScdType, corr_fact); } - } + + // ------------------------------------------------------------------------- + // Fake Lambda reconstruction for combinatorial background (LS pairs) + // ------------------------------------------------------------------------- + // subsample before writing to the output table. + std::vector fakeCandBuffer; + fakeCandBuffer.reserve(512); + + for (auto const& tr1 : tracks) { + if (!selTrackBasic(tr1)) { + continue; + } + + for (auto const& tr2 : tracks) { + if (tr2.index() <= tr1.index()) { + continue; + } + + if (!selTrackBasic(tr2)) { + continue; + } + + if (tr1.sign() != tr2.sign()) { + continue; + } + + if (cFakeDaughterMaxDEta > 0.f && + std::abs(tr1.eta() - tr2.eta()) > cFakeDaughterMaxDEta) { + continue; + } + if (cFakeDaughterMaxOpenAngle > 0.f) { + float cosOA = (tr1.px() * tr2.px() + tr1.py() * tr2.py() + tr1.pz() * tr2.pz()) / + (std::sqrt(tr1.px() * tr1.px() + tr1.py() * tr1.py() + tr1.pz() * tr1.pz() + 1e-10f) * + std::sqrt(tr2.px() * tr2.px() + tr2.py() * tr2.py() + tr2.pz() * tr2.pz() + 1e-10f)); + if (std::acos(std::clamp(cosOA, -1.f, 1.f)) > cFakeDaughterMaxOpenAngle) { + continue; + } + } + + float fkPx = tr1.px() + tr2.px(); + float fkPy = tr1.py() + tr2.py(); + float fkPz = tr1.pz() + tr2.pz(); + float fkPt = std::sqrt(fkPx * fkPx + fkPy * fkPy); + float fkEta = RecoDecay::eta(std::array{fkPx, fkPy, fkPz}); + float fkPhi = RecoDecay::phi(fkPx, fkPy); + + ParticleType fakeType = (tr1.sign() > 0) ? kFakeLambda : kFakeAntiLambda; + + auto bufferFakeCandidate = [&](float prPxL, float prPyL, float prPzL, float fkMass, int64_t prId, int64_t piId) { + float fkRap = RecoDecay::y(std::array{fkPx, fkPy, fkPz}, fkMass); + float fkRapEta = cDoEtaAnalysis ? std::abs(fkEta) : std::abs(fkRap); + if (!kinCutSelection(fkPt, fkRapEta, cMinV0Pt, cMaxV0Pt, cMaxV0Rap)) + return; + fakeCandBuffer.push_back({fkPx, fkPy, fkPz, fkPt, fkEta, fkPhi, fkRap, fkMass, prPxL, prPyL, prPzL, prId, piId, fakeType}); + }; + + // Assignment A: tr1 = proton candidate, tr2 = pion candidate + if (std::abs(tr1.tpcNSigmaPr()) < cTpcNsigmaCut && + std::abs(tr2.tpcNSigmaPi()) < cTpcNsigmaCut && + selTrackBasicWithDCA(tr1, true /*isProton*/) && + selTrackBasicWithDCA(tr2, false /*isPion*/)) { + std::array pPr = {tr1.px(), tr1.py(), tr1.pz()}; + std::array pPi = {tr2.px(), tr2.py(), tr2.pz()}; + float fkMassA = RecoDecay::m(std::array{pPr, pPi}, std::array{MassProton, MassPionCharged}); + if (fkMassA >= cMinV0Mass && fkMassA <= cMaxV0Mass) { + bufferFakeCandidate(tr1.px(), tr1.py(), tr1.pz(), fkMassA, tr1.index(), tr2.index()); + } + } + + // Assignment B: tr2 = proton candidate, tr1 = pion candidate + if (cAllowDualAssignment && + std::abs(tr2.tpcNSigmaPr()) < cTpcNsigmaCut && + std::abs(tr1.tpcNSigmaPi()) < cTpcNsigmaCut && + selTrackBasicWithDCA(tr2, true /*isProton*/) && + selTrackBasicWithDCA(tr1, false /*isPion*/)) { + std::array pPr = {tr2.px(), tr2.py(), tr2.pz()}; + std::array pPi = {tr1.px(), tr1.py(), tr1.pz()}; + float fkMassB = RecoDecay::m(std::array{pPr, pPi}, std::array{MassProton, MassPionCharged}); + if (fkMassB >= cMinV0Mass && fkMassB <= cMaxV0Mass) { + bufferFakeCandidate(tr2.px(), tr2.py(), tr2.pz(), fkMassB, tr2.index(), tr1.index()); + } + } + } // end tr2 loop + } // end tr1 loop + + int nFakeTotal = static_cast(fakeCandBuffer.size()); + if (cMaxFakeLambdaPerEvent > 0 && nFakeTotal > cMaxFakeLambdaPerEvent) { + if (cRandomizeFakeSubsample) { + static thread_local std::mt19937 rng{std::random_device{}()}; + std::shuffle(fakeCandBuffer.begin(), fakeCandBuffer.end(), rng); + } + fakeCandBuffer.resize(static_cast(cMaxFakeLambdaPerEvent)); + } + int nFakeKept = static_cast(fakeCandBuffer.size()); + if (nFakeTotal > 0) { + histos.fill(HIST("Tracks/h1f_fake_subsample_fraction"), static_cast(nFakeKept) / static_cast(nFakeTotal)); + } + + for (auto const& fc : fakeCandBuffer) { + if (fc.type == kFakeLambda) { + histos.fill(HIST("Tracks/h1f_lambda_pt_vs_invm"), fc.mass, fc.pt); + } else { + histos.fill(HIST("Tracks/h1f_antilambda_pt_vs_invm"), fc.mass, fc.pt); + } + lambdaTrackTable(lambdaCollisionTable.lastIndex(), + fc.px, fc.py, fc.pz, fc.pt, fc.eta, fc.phi, fc.rap, fc.mass, + fc.prPx, fc.prPy, fc.prPz, + fc.id1, fc.id2, + -999.f, -999.f, (int8_t)fc.type, kPrimary, 1.f); + } + } // end fillLambdaRecoTables // MC Generater Level Tables template @@ -1313,37 +1486,6 @@ struct LambdaTableProducer { PROCESS_SWITCH(LambdaTableProducer, processDataRun3, "Process for Run3 DATA", true); - /* void processDataRun2(CollisionsRun2::iterator const& collision, aod::V0Datas const& V0s, TracksRun2 const& tracks) - { - fillLambdaRecoTables(collision, V0s, tracks); - } - - PROCESS_SWITCH(LambdaTableProducer, processDataRun2, "Process for Run2 DATA", false); - - void processMCRecoRun3(soa::Join::iterator const& collision, aod::McCollisions const&, - McV0Tracks const& V0s, TracksMC const& tracks, aod::McParticles const&) - { - // check collision - if (!selCollision(collision)) { - return; - } - fillLambdaRecoTables(collision, V0s, tracks); - } - - PROCESS_SWITCH(LambdaTableProducer, processMCRecoRun3, "Process for Run3 McReco DATA", false); - - void processMCRecoRun2(soa::Join::iterator const& collision, aod::McCollisions const&, - McV0Tracks const& V0s, TracksMCRun2 const& tracks, aod::McParticles const&) - { - // check collision - if (!selCollision(collision)) { - return; - } - fillLambdaRecoTables(collision, V0s, tracks); - } */ - - // PROCESS_SWITCH(LambdaTableProducer, processMCRecoRun2, "Process for Run2 McReco DATA", false); - void processMCRun3(aod::McCollisions::iterator const& mcCollision, soa::SmallGroups> const& collisions, McV0Tracks const& V0s, TracksMC const& tracks, @@ -1396,9 +1538,10 @@ struct LambdaTracksExtProducer { histos.add("h1i_totantilambda_mult", "Multiplicity", kTH1I, {axisMult}); histos.add("h1i_lambda_mult", "Multiplicity", kTH1I, {axisMult}); histos.add("h1i_antilambda_mult", "Multiplicity", kTH1I, {axisMult}); - histos.add("h2d_n2_etaphi_LaP_LaM", "#rho_{2}^{SharePair}", kTH2D, {axisDEta, axisDPhi}); - histos.add("h2d_n2_etaphi_LaP_LaP", "#rho_{2}^{SharePair}", kTH2D, {axisDEta, axisDPhi}); - histos.add("h2d_n2_etaphi_LaM_LaM", "#rho_{2}^{SharePair}", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaP_LaM", "#rho_{2}^{SharePair} #Lambda#bar{#Lambda}", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaM_LaP", "#rho_{2}^{SharePair} #bar{#Lambda}#Lambda", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaP_LaP", "#rho_{2}^{SharePair} #Lambda#Lambda", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaM_LaM", "#rho_{2}^{SharePair} #bar{#Lambda}#bar{#Lambda}", kTH2D, {axisDEta, axisDPhi}); // InvMass, DcaDau and CosPA histos.add("Reco/h1f_lambda_invmass", "M_{p#pi}", kTH1F, {axisMass}); @@ -1444,7 +1587,10 @@ struct LambdaTracksExtProducer { ++nTotAntiLambda; } - tLambda = (cA * std::abs(lambda.mass() - MassLambda0)) + (cB * lambda.dcaDau()) + (cC * std::abs(lambda.cosPA() - 1.)); + bool const isFake = (lambda.v0Type() == kFakeLambda || lambda.v0Type() == kFakeAntiLambda); + if (!isFake) { + tLambda = (cA * std::abs(lambda.mass() - MassLambda0)) + (cB * lambda.dcaDau()) + (cC * std::abs(lambda.cosPA() - 1.)); + } for (auto const& track : tracks) { // check lambda index (don't analyze same lambda track !!!) @@ -1457,13 +1603,16 @@ struct LambdaTracksExtProducer { vSharedDauLambdaIndex.push_back(track.index()); lambdaSharingDauFlag = true; - // Fill DEta-DPhi Histogram - if ((lambda.v0Type() == kLambda && track.v0Type() == kAntiLambda) || (lambda.v0Type() == kAntiLambda && track.v0Type() == kLambda)) { - histos.fill(HIST("h2d_n2_etaphi_LaP_LaM"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); - } else if (lambda.v0Type() == kLambda && track.v0Type() == kLambda) { - histos.fill(HIST("h2d_n2_etaphi_LaP_LaP"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); - } else if (lambda.v0Type() == kAntiLambda && track.v0Type() == kAntiLambda) { - histos.fill(HIST("h2d_n2_etaphi_LaM_LaM"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); + if (!isFake) { + if (lambda.v0Type() == kLambda && track.v0Type() == kAntiLambda) { + histos.fill(HIST("h2d_n2_etaphi_LaP_LaM"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); + } else if (lambda.v0Type() == kAntiLambda && track.v0Type() == kLambda) { + histos.fill(HIST("h2d_n2_etaphi_LaM_LaP"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); + } else if (lambda.v0Type() == kLambda && track.v0Type() == kLambda) { + histos.fill(HIST("h2d_n2_etaphi_LaP_LaP"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); + } else if (lambda.v0Type() == kAntiLambda && track.v0Type() == kAntiLambda) { + histos.fill(HIST("h2d_n2_etaphi_LaM_LaM"), lambda.eta() - track.eta(), RecoDecay::constrainAngle((lambda.phi() - track.phi()), -PIHalf)); + } } // decision based on mass closest to PdgMass of Lambda @@ -1472,14 +1621,35 @@ struct LambdaTracksExtProducer { } // decisions based on t-score - tTrack = (cA * std::abs(track.mass() - MassLambda0)) + (cB * track.dcaDau()) + (cC * std::abs(track.cosPA() - 1.)); - if (tLambda > tTrack) { - lambdaMinTScoreFlag = false; + // Guard: only compute tTrack for real V0 candidates (fakes have sentinel cosPA/dcaDau) + bool const isFakeTrack = (track.v0Type() == kFakeLambda || track.v0Type() == kFakeAntiLambda); + if (!isFake && !isFakeTrack) { + tTrack = (cA * std::abs(track.mass() - MassLambda0)) + (cB * track.dcaDau()) + (cC * std::abs(track.cosPA() - 1.)); + if (tLambda > tTrack) { + lambdaMinTScoreFlag = false; + } } } } - // fill QA histograms + if (isFake) { + bool fakeAccepted = false; + if (cAcceptAllLambda) { + fakeAccepted = true; + } else if (cRejAllLambdaShaDau && !lambdaSharingDauFlag) { + fakeAccepted = true; + } + if (fakeAccepted) { + if (lambda.v0Type() == kFakeLambda) + ++nSelLambda; + else + ++nSelAntiLambda; + } + lambdaTrackExtTable(lambdaSharingDauFlag, vSharedDauLambdaIndex, false /*trueLambdaFlag*/, fakeAccepted /*fakeFlag*/); + continue; + } + + // fill QA histograms (real V0s only) if (lambdaSharingDauFlag) { fillHistos(lambda); } else { @@ -1505,8 +1675,8 @@ struct LambdaTracksExtProducer { } } - // fill LambdaTrackExt table - lambdaTrackExtTable(lambdaSharingDauFlag, vSharedDauLambdaIndex, trueLambdaFlag); + // fill LambdaTrackExt table (fakeFlag = false for all real V0 candidates) + lambdaTrackExtTable(lambdaSharingDauFlag, vSharedDauLambdaIndex, trueLambdaFlag, false); } // fill multiplicity histograms @@ -1536,15 +1706,21 @@ struct LambdaSpinPolarization { // Global Configurables Configurable cNPtBins{"cNPtBins", 30, "N pT Bins"}; Configurable cMinPt{"cMinPt", 0.5, "pT Min"}; - Configurable cMaxPt{"cMaxPt", 3.5, "pT Max"}; + Configurable cMaxPt{"cMaxPt", 4.5, "pT Max"}; Configurable cNRapBins{"cNRapBins", 10, "N Rapidity Bins"}; Configurable cMinRap{"cMinRap", -0.5, "Minimum Rapidity"}; Configurable cMaxRap{"cMaxRap", 0.5, "Maximum Rapidity"}; Configurable cNPhiBins{"cNPhiBins", 36, "N Phi Bins"}; - Configurable cNBinsCosTS{"cNBinsCosTS", 10, "N CosTS Bins"}; + Configurable cNBinsCosTS{"cNBinsCosTS", 20, "N CosTS Bins"}; + Configurable cNBinsDeltaR{"cNBinsDeltaR", 20, "DeltaR Bins"}; Configurable cInvBoostFlag{"cInvBoostFlag", true, "Inverse Boost Flag"}; Configurable mixingParameter{"mixingParameter", 5, "How many events are mixed"}; + Configurable cDoAtlasMethod{"cDoAtlasMethod", false, "Fill pair-boost histograms"}; + Configurable cDoStarMethod{"cDoStarMethod", true, "Fill lab-boost histograms"}; + + Configurable cMEMode{"cMEMode", 0, "ME mode: 0=standard, 1=sequential"}; + // Centrality Axis ConfigurableAxis cMultBins{"cMultBins", {VARIABLE_WIDTH, 0.0f, 10.0f, 30.0f, 50.f, 80.0f, 100.f}, "Variable Mult-Bins"}; ConfigurableAxis axisCentME{"axisCentME", {VARIABLE_WIDTH, 0, 10, 30, 50, 100}, "Mixing bins - centrality (%)"}; @@ -1563,11 +1739,11 @@ struct LambdaSpinPolarization { void init(InitContext const&) { const AxisSpec axisCheck(1, 0, 1, ""); - const AxisSpec axisPosZ(220, -11, 11, "V_{z} (cm)"); + const AxisSpec axisPosZ(220, -7, 7, "V_{z} (cm)"); const AxisSpec axisCent(cMultBins, "FT0M (%)"); const AxisSpec axisChMult(200, 0, 200, "N_{ch}"); const AxisSpec axisMult(10, 0, 10, "N_{#Lambda}"); - const AxisSpec axisMass(100, 1.06, 1.16, "M_{#Lambda} (GeV/#it{c}^{2})"); + const AxisSpec axisMass(100, 1.08, 1.2, "M_{#Lambda} (GeV/#it{c}^{2})"); const AxisSpec axisPt(cNPtBins, cMinPt, cMaxPt, "p_{T} (GeV/#it{c})"); const AxisSpec axisEta(cNRapBins, cMinRap, cMaxRap, "#eta"); const AxisSpec axisRap(cNRapBins, cMinRap, cMaxRap, "y"); @@ -1575,7 +1751,7 @@ struct LambdaSpinPolarization { const AxisSpec axisDRap(2 * cNRapBins, cMinRap - cMaxRap, cMaxRap - cMinRap, "#Deltay"); const AxisSpec axisDPhi(cNPhiBins, -PI, PI, "#Delta#varphi"); const AxisSpec axisCosTS(cNBinsCosTS, -1, 1, "cos(#theta*)"); - const AxisSpec axisDR(10, 0, 2, "#DeltaR"); + const AxisSpec axisDR(cNBinsDeltaR, 0, 4, "#DeltaR"); // Pool occupancy histos.add("QA/ME/hPoolCentVz", "ME pool occupancy;centrality (%);V_{z} (cm)", kTH2F, {axisCentME, axisVtxZME}); @@ -1584,18 +1760,36 @@ struct LambdaSpinPolarization { histos.add("QA/ME/hLambdaMultVsCent", "ME #Lambda multiplicity;centrality (%);N_{#Lambda}", kTH2F, {axisCentME, {50, 0, 50}}); histos.add("QA/ME/hAntiLambdaMultVsCent", "ME #bar{#Lambda} multiplicity;centrality (%);N_{#bar{#Lambda}}", kTH2F, {axisCentME, {50, 0, 50}}); - // inv mass vs pt for Lambda and AntiLambda + // inv mass vs pt: four separate signal pair types histos.add("Reco/h2f_n2_mass_LaPLaM", "m_{inv}^{#Lambda} vs m_{inv}^{#bar{#Lambda}}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); + histos.add("Reco/h2f_n2_mass_LaMLaP", "m_{inv}^{#bar{#Lambda}} vs m_{inv}^{#Lambda}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); histos.add("Reco/h2f_n2_mass_LaPLaP", "m_{inv}^{#Lambda} vs m_{inv}^{#Lambda}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); histos.add("Reco/h2f_n2_mass_LaMLaM", "m_{inv}^{#bar{#Lambda}} vs m_{inv}^{#bar{#Lambda}}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); - // rho2 for C2 + histos.add("RecoBkg/h2f_n2_mass_LaPFkLaM", "US-LS: m_{inv}^{#Lambda} vs m_{inv}^{fake#bar{#Lambda}}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); + histos.add("RecoBkg/h2f_n2_mass_LaMLFkLaP", "US-LS: m_{inv}^{#bar{#Lambda}} vs m_{inv}^{fake#Lambda}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); + histos.add("RecoBkg/h2f_n2_mass_LaPFkLaP", "US-LS: m_{inv}^{#Lambda} vs m_{inv}^{fake#Lambda}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); + histos.add("RecoBkg/h2f_n2_mass_LaMLFkLaM", "US-LS: m_{inv}^{#bar{#Lambda}} vs m_{inv}^{fake#bar{#Lambda}}", kTHnSparseF, {axisMass, axisMass, axisPt, axisPt}); + histos.add("RecoCorr/h2f_n2_dltaR_LaPLaM", "#rho_{2}^{#Lambda#bar{#Lambda}}", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + histos.add("RecoCorr/h2f_n2_dltaR_LaMLaP", "#rho_{2}^{#bar{#Lambda}#Lambda}", kTHnSparseF, {axisCent, axisDR, axisCosTS}); histos.add("RecoCorr/h2f_n2_dltaR_LaPLaP", "#rho_{2}^{#Lambda#Lambda}", kTHnSparseF, {axisCent, axisDR, axisCosTS}); histos.add("RecoCorr/h2f_n2_dltaR_LaMLaM", "#rho_{2}^{#bar{#Lambda}#bar{#Lambda}}", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + histos.add("RecoCorr/h2f_n2_ctheta_LaPLaM", "#rho_{2}^{#Lambda#bar{#Lambda}}", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); + histos.add("RecoCorr/h2f_n2_ctheta_LaMLaP", "#rho_{2}^{#bar{#Lambda}#Lambda}", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); histos.add("RecoCorr/h2f_n2_ctheta_LaPLaP", "#rho_{2}^{#Lambda#Lambda}", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); histos.add("RecoCorr/h2f_n2_ctheta_LaMLaM", "#rho_{2}^{#bar{#Lambda}#bar{#Lambda}}", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); + + histos.add("RecoCorrBkg/h2f_n2_dltaR_LaPFkLaM", "#rho_{2}^{#Lambda,fake#bar{#Lambda}} (US-LS)", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_dltaR_LaMLFkLaP", "#rho_{2}^{#bar{#Lambda},fake#Lambda} (US-LS)", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_dltaR_LaPFkLaP", "#rho_{2}^{#Lambda,fake#Lambda} (US-LS)", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_dltaR_LaMLFkLaM", "#rho_{2}^{#bar{#Lambda},fake#bar{#Lambda}} (US-LS)", kTHnSparseF, {axisCent, axisDR, axisCosTS}); + + histos.add("RecoCorrBkg/h2f_n2_ctheta_LaPFkLaM", "#rho_{2}^{#Lambda,fake#bar{#Lambda}} (US-LS)", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_ctheta_LaMLFkLaP", "#rho_{2}^{#bar{#Lambda},fake#Lambda} (US-LS)", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_ctheta_LaPFkLaP", "#rho_{2}^{#Lambda,fake#Lambda} (US-LS)", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); + histos.add("RecoCorrBkg/h2f_n2_ctheta_LaMLFkLaM", "#rho_{2}^{#bar{#Lambda},fake#bar{#Lambda}} (US-LS)", kTHnSparseF, {axisCent, axisDRap, axisDPhi, axisCosTS}); } void getBoostVector(std::array const& p, std::array& v, bool inverseBoostFlag = true) @@ -1626,82 +1820,179 @@ struct LambdaSpinPolarization { template void fillPairHistos(U& p1, U& p2) { - static constexpr std::string_view SubDirHist[] = {"LaPLaM", "LaPLaP", "LaMLaM"}; + static constexpr std::string_view SubDirHistUS[] = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; + static constexpr std::string_view SubDirHistBkg[] = {"LaPFkLaM", "LaMLFkLaP", "LaPFkLaP", "LaMLFkLaM"}; + + constexpr bool isBkg = (part_pair == kLambdaFakeAntiLambda || part_pair == kAntiLambdaFakeLambda || part_pair == kLambdaFakeLambda || part_pair == kAntiLambdaFakeAntiLambda); + + // Fill pair invariant mass histogram + if constexpr (!isBkg) { + histos.fill(HIST("Reco/h2f_n2_mass_") + HIST(SubDirHistUS[part_pair]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); + } else { + constexpr int bkgIdx = (int)part_pair - 4; + histos.fill(HIST("RecoBkg/h2f_n2_mass_") + HIST(SubDirHistBkg[bkgIdx]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); + } - // Fill lambda pair mass - // histos.fill(HIST("Reco/h2f_n2_mass_") + HIST(SubDirHist[part_pair]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); float drap = p1.rap() - p2.rap(); float dphi = RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PI); float dR = std::sqrt(drap * drap + dphi * dphi); std::array l1 = {p1.px(), p1.py(), p1.pz(), MassLambda0}; std::array l2 = {p2.px(), p2.py(), p2.pz(), MassLambda0}; - - std::array llpair = { - l1[0] + l2[0], - l1[1] + l2[1], - l1[2] + l2[2], - l1[3] + l2[3]}; - std::array pr1 = {p1.prPx(), p1.prPy(), p1.prPz(), MassProton}; std::array pr2 = {p2.prPx(), p2.prPy(), p2.prPz(), MassProton}; - std::array vPair; - getBoostVector(llpair, vPair, cInvBoostFlag); - - boost(l1, vPair); - boost(l2, vPair); - - boost(pr1, vPair); - boost(pr2, vPair); + if (cDoAtlasMethod) { + std::array l1_atlas = l1; + std::array l2_atlas = l2; + std::array pr1_atlas = pr1; + std::array pr2_atlas = pr2; + + std::array llpair = { + l1_atlas[0] + l2_atlas[0], + l1_atlas[1] + l2_atlas[1], + l1_atlas[2] + l2_atlas[2], + l1_atlas[3] + l2_atlas[3]}; + + std::array vPair; + getBoostVector(llpair, vPair, cInvBoostFlag); + boost(l1_atlas, vPair); + boost(l2_atlas, vPair); + boost(pr1_atlas, vPair); + boost(pr2_atlas, vPair); + + std::array v1_pair, v2_pair; + getBoostVector(l1_atlas, v1_pair, cInvBoostFlag); + getBoostVector(l2_atlas, v2_pair, cInvBoostFlag); + boost(pr1_atlas, v1_pair); + boost(pr2_atlas, v2_pair); + + std::array pr1tv_atlas = {pr1_atlas[0], pr1_atlas[1], pr1_atlas[2]}; + std::array pr2tv_atlas = {pr2_atlas[0], pr2_atlas[1], pr2_atlas[2]}; + float ctheta = + RecoDecay::dotProd(pr1tv_atlas, pr2tv_atlas) / + (RecoDecay::sqrtSumOfSquares(pr1tv_atlas[0], pr1tv_atlas[1], pr1tv_atlas[2]) * + RecoDecay::sqrtSumOfSquares(pr2tv_atlas[0], pr2tv_atlas[1], pr2tv_atlas[2])); + + if constexpr (!isBkg) { + histos.fill(HIST("RecoCorr/h2f_n2_ctheta_") + HIST(SubDirHistUS[part_pair]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorr/h2f_n2_dltaR_") + HIST(SubDirHistUS[part_pair]), cent, dR, ctheta); + } else { + constexpr int bkgIdx = (int)part_pair - 4; + histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[bkgIdx]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[bkgIdx]), cent, dR, ctheta); + } + } - std::array v1_pair, v2_pair; - getBoostVector(l1, v1_pair, cInvBoostFlag); - getBoostVector(l2, v2_pair, cInvBoostFlag); + if (cDoStarMethod) { + std::array pr1_star = pr1; + std::array pr2_star = pr2; - boost(pr1, v1_pair); - boost(pr2, v2_pair); + std::array v1_lab, v2_lab; + getBoostVector(l1, v1_lab, cInvBoostFlag); + getBoostVector(l2, v2_lab, cInvBoostFlag); + boost(pr1_star, v1_lab); + boost(pr2_star, v2_lab); - std::array pr1tv = {pr1[0], pr1[1], pr1[2]}; - std::array pr2tv = {pr2[0], pr2[1], pr2[2]}; - float ctheta = RecoDecay::dotProd(pr1tv, pr2tv) / (RecoDecay::sqrtSumOfSquares(pr1tv[0], pr1tv[1], pr1tv[2]) * RecoDecay::sqrtSumOfSquares(pr2tv[0], pr2tv[1], pr2tv[2])); + std::array pr1tv_star = {pr1_star[0], pr1_star[1], pr1_star[2]}; + std::array pr2tv_star = {pr2_star[0], pr2_star[1], pr2_star[2]}; + float ctheta = + RecoDecay::dotProd(pr1tv_star, pr2tv_star) / + (RecoDecay::sqrtSumOfSquares(pr1tv_star[0], pr1tv_star[1], pr1tv_star[2]) * + RecoDecay::sqrtSumOfSquares(pr2tv_star[0], pr2tv_star[1], pr2tv_star[2])); - histos.fill(HIST("RecoCorr/h2f_n2_ctheta_") + HIST(SubDirHist[part_pair]), cent, drap, dphi, ctheta); - histos.fill(HIST("RecoCorr/h2f_n2_dltaR_") + HIST(SubDirHist[part_pair]), cent, dR, ctheta); + if constexpr (!isBkg) { + histos.fill(HIST("RecoCorr/h2f_n2_ctheta_") + HIST(SubDirHistUS[part_pair]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorr/h2f_n2_dltaR_") + HIST(SubDirHistUS[part_pair]), cent, dR, ctheta); + } else { + constexpr int bkgIdx = (int)part_pair - 4; + histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[bkgIdx]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[bkgIdx]), cent, dR, ctheta); + } + } } - template - void analyzePairs(T const& trks_1, T const& trks_2) + template + void analyzePairsStandard(T const& trks_1, T const& trks_2) { for (auto const& trk_1 : trks_1) { for (auto const& trk_2 : trks_2) { - // check for same index for Lambda-Lambda / AntiLambda-AntiLambda - if (samelambda && ((trk_1.index() == trk_2.index()))) { - continue; + + if constexpr (samelambda) { + if (trk_1.index() == trk_2.index()) + continue; } - fillPairHistos(trk_1, trk_2); + + fillPairHistos(trk_1, trk_2); } } } template - void analyzePairsWithKinematicMatching(T const& trks_1, T const& trks_2) + void analyzePairsSequential(T const& trks_1, T const& trks_2) { + std::vector candidates; + candidates.reserve(trks_2.size()); + + for (auto const& trk : trks_2) { + candidates.push_back(&trk); + } + + if (candidates.empty()) { + return; + } + + // Randomize candidates + static thread_local std::mt19937 rng(std::random_device{}()); + std::shuffle(candidates.begin(), candidates.end(), rng); + for (auto const& trk_1 : trks_1) { - for (auto const& trk_2 : trks_2) { - if (samelambda && ((trk_1.index() == trk_2.index()))) { + + decltype(&*trks_2.begin()) currentPartner = nullptr; + + for (auto const* trk_2 : candidates) { + + if constexpr (samelambda) { + if (trk_1.index() == trk_2->index()) + continue; + } + + // Initialize + if (!currentPartner) { + currentPartner = trk_2; continue; } - // Kinematic matching - float deltaPt = std::abs(trk_1.pt() - trk_2.pt()); - float deltaPhi = std::abs(RecoDecay::constrainAngle(trk_1.phi() - trk_2.phi(), -PI)); - float deltaRap = std::abs(trk_1.rap() - trk_2.rap()); + float dPt = std::abs(trk_2->pt() - currentPartner->pt()); + float dPhi = std::abs(RecoDecay::constrainAngle(trk_2->phi() - currentPartner->phi(), -PI)); + float dRap = std::abs(trk_2->rap() - currentPartner->rap()); - if (deltaPt < cMaxDeltaPt && deltaPhi < cMaxDeltaPhi && deltaRap < cMaxDeltaRap) { - fillPairHistos(trk_1, trk_2); + if (dPt < cMaxDeltaPt && dPhi < cMaxDeltaPhi && dRap < cMaxDeltaRap) { + currentPartner = trk_2; } } + + if (currentPartner) { + fillPairHistos(trk_1, *currentPartner); + } + } + } + + template + void analyzePairsME(T const& trks_1, T const& trks_2) + { + switch (cMEMode) { + + case 0: // standard + analyzePairsStandard(trks_1, trks_2); + break; + + case 1: // sequential + analyzePairsSequential(trks_1, trks_2); + break; + + default: + LOGF(fatal, "Invalid cMEMode value!"); } } @@ -1712,27 +2003,48 @@ struct LambdaSpinPolarization { Preslice perCollisionLambda = aod::lambdatrack::lambdaCollisionId; SliceCache cache; - Partition partLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kLambda) && (aod::lambdatrackext::trueLambdaFlag == true) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); - Partition partAntiLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kAntiLambda) && (aod::lambdatrackext::trueLambdaFlag == true) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); + // Real V0 partitions: trueLambdaFlag=true (passed sharing-daughter resolution) AND fakeFlag=false + Partition partLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kLambda) && (aod::lambdatrackext::trueLambdaFlag == true) && (aod::lambdatrackext::fakeFlag == false) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); + Partition partAntiLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kAntiLambda) && (aod::lambdatrackext::trueLambdaFlag == true) && (aod::lambdatrackext::fakeFlag == false) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); + + // LS combinatorial background partitions: + // fakeFlag=true now means: passed the sharing-daughter acceptance logic (Issue 3 fix). + Partition partFakeLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kFakeLambda) && (aod::lambdatrackext::fakeFlag == true); + Partition partFakeAntiLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kFakeAntiLambda) && (aod::lambdatrackext::fakeFlag == true); void processDummy(LambdaCollisions::iterator const&) {} - PROCESS_SWITCH(LambdaSpinPolarization, processDummy, "Dummy process", true); + PROCESS_SWITCH(LambdaSpinPolarization, processDummy, "Dummy process", false); void processDataReco(LambdaCollisions::iterator const& collision, LambdaTracks const&) { - // assign centrality cent = collision.cent(); auto lambdaTracks = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); auto antiLambdaTracks = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); - // Add QA for single Lambds + auto fakeLambdaTracks = partFakeLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto fakeAntiLambdaTracks = partFakeAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); - // Analyze pairs - analyzePairs(lambdaTracks, antiLambdaTracks); - analyzePairs(lambdaTracks, lambdaTracks); - analyzePairs(antiLambdaTracks, antiLambdaTracks); + analyzePairsStandard(lambdaTracks, antiLambdaTracks); + analyzePairsStandard(antiLambdaTracks, lambdaTracks); + analyzePairsStandard(lambdaTracks, lambdaTracks); + analyzePairsStandard(antiLambdaTracks, antiLambdaTracks); + + // --------------------------------------------------------------------------- + // US-LS combinatorial background pairs + // --------------------------------------------------------------------------- + analyzePairsStandard(lambdaTracks, fakeAntiLambdaTracks); + analyzePairsStandard(fakeLambdaTracks, antiLambdaTracks); + + analyzePairsStandard(antiLambdaTracks, fakeLambdaTracks); + analyzePairsStandard(fakeAntiLambdaTracks, lambdaTracks); + + analyzePairsStandard(lambdaTracks, fakeLambdaTracks); + analyzePairsStandard(fakeLambdaTracks, lambdaTracks); + + analyzePairsStandard(antiLambdaTracks, fakeAntiLambdaTracks); + analyzePairsStandard(fakeAntiLambdaTracks, antiLambdaTracks); } PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, "Process for Data and MCReco", true); @@ -1771,33 +2083,77 @@ struct LambdaSpinPolarization { auto antiLambdaTracks_col1 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); auto antiLambdaTracks_col2 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + // Fake lambda slices (LS background) + auto fakeLambdaTracks_col2 = partFakeLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + auto fakeAntiLambdaTracks_col2 = partFakeAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + // QA: multiplicities histos.fill(HIST("QA/ME/hLambdaMultVsCent"), col1.cent(), lambdaTracks_col1.size()); histos.fill(HIST("QA/ME/hAntiLambdaMultVsCent"), col1.cent(), antiLambdaTracks_col1.size()); - // Mixed-event pairs - analyzePairsWithKinematicMatching(lambdaTracks_col1, antiLambdaTracks_col2); - analyzePairsWithKinematicMatching(antiLambdaTracks_col1, lambdaTracks_col2); - analyzePairsWithKinematicMatching(lambdaTracks_col1, lambdaTracks_col2); - analyzePairsWithKinematicMatching(antiLambdaTracks_col1, antiLambdaTracks_col2); + analyzePairsME(lambdaTracks_col1, antiLambdaTracks_col2); + analyzePairsME(antiLambdaTracks_col1, lambdaTracks_col2); + analyzePairsME(lambdaTracks_col1, lambdaTracks_col2); + analyzePairsME(antiLambdaTracks_col1, antiLambdaTracks_col2); + + analyzePairsME(lambdaTracks_col1, fakeAntiLambdaTracks_col2); + analyzePairsME(fakeLambdaTracks_col2, antiLambdaTracks_col1); + + analyzePairsME(antiLambdaTracks_col1, fakeLambdaTracks_col2); + analyzePairsME(fakeAntiLambdaTracks_col2, lambdaTracks_col1); + + analyzePairsME(lambdaTracks_col1, fakeLambdaTracks_col2); + analyzePairsME(fakeLambdaTracks_col2, lambdaTracks_col1); + + analyzePairsME(antiLambdaTracks_col1, fakeAntiLambdaTracks_col2); + analyzePairsME(fakeAntiLambdaTracks_col2, antiLambdaTracks_col1); } } PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixed, "Process for Data and MCReco for Mixed events", false); - void processDataRecoMixEvent(LambdaCollisions::iterator const& collision, LambdaTracks const& tracks) + void processDataRecoMixEvent(LambdaCollisions::iterator const& collision, LambdaTracks const&) { - // return for no lambdas in a collision - if (tracks.size() == 0) { + // Slice each particle type for this collision + auto lambdaTracks = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto antiLambdaTracks = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto fakeLambdaTracks = partFakeLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto fakeAntiLambdaTracks = partFakeAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + + // Skip collision if there are no true lambdas of any kind + if (lambdaTracks.size() == 0 && antiLambdaTracks.size() == 0 && fakeLambdaTracks.size() == 0 && fakeAntiLambdaTracks.size() == 0) { return; } - // fill collision table + // Fill collision table once per collision lambdaMixEvtCol(collision.index(), collision.cent(), collision.posZ(), collision.timeStamp()); - for (auto const& track : tracks) { - lambdaMixEvtTrk(collision.index(), track.globalIndex(), track.px(), track.py(), track.pz(), track.mass(), - track.prPx(), track.prPy(), track.prPz(), track.v0Type(), collision.timeStamp()); + for (auto const& track : lambdaTracks) { + lambdaMixEvtTrk(collision.index(), track.globalIndex(), + track.px(), track.py(), track.pz(), track.mass(), + track.prPx(), track.prPy(), track.prPz(), + track.v0Type(), collision.timeStamp()); + } + + for (auto const& track : antiLambdaTracks) { + lambdaMixEvtTrk(collision.index(), track.globalIndex(), + track.px(), track.py(), track.pz(), track.mass(), + track.prPx(), track.prPy(), track.prPz(), + track.v0Type(), collision.timeStamp()); + } + + for (auto const& track : fakeLambdaTracks) { + lambdaMixEvtTrk(collision.index(), track.globalIndex(), + track.px(), track.py(), track.pz(), track.mass(), + track.prPx(), track.prPy(), track.prPz(), + track.v0Type(), collision.timeStamp()); + } + + for (auto const& track : fakeAntiLambdaTracks) { + lambdaMixEvtTrk(collision.index(), track.globalIndex(), + track.px(), track.py(), track.pz(), track.mass(), + track.prPx(), track.prPy(), track.prPz(), + track.v0Type(), collision.timeStamp()); } } From f16e574676f7975ca659e5b9a9b1485347a505ef Mon Sep 17 00:00:00 2001 From: Subhadeep Roy Date: Fri, 10 Apr 2026 00:22:19 +0530 Subject: [PATCH 2/4] Refactor lambda spin polarization analysis code --- .../Tasks/lambdaSpinPolarization.cxx | 24 +++++++------------ 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx index 3ccfe9c3032..6141a916c2d 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx @@ -258,7 +258,7 @@ enum ParticlePairType { kLambdaFakeAntiLambda = 4, kAntiLambdaFakeLambda = 5, kLambdaFakeLambda = 6, - kAntiLambdaFakeAntiLambda = 7 + kAntiLambdaFakeAntiLambda = 7 }; enum ShareDauLambda { @@ -328,7 +328,6 @@ struct LambdaTableProducer { Configurable cMaxChi2Tpc{"cMaxChi2Tpc", 4, "Max Chi2 Tpc"}; Configurable cTpcNsigmaCut{"cTpcNsigmaCut", 3.0, "TPC NSigma Selection Cut"}; Configurable cRemoveAmbiguousTracks{"cRemoveAmbiguousTracks", false, "Remove Ambiguous Tracks"}; - // LS (fake) background control variables Configurable cApplyFakeDcaCuts{"cApplyFakeDcaCuts", true, "Apply DCA-to-PV cuts on LS daughter tracks"}; Configurable cFakeDaughterMaxDEta{"cFakeDaughterMaxDEta", -1.f, "Max delta eta between LS daughters as topology proxy (-1=off)"}; @@ -1235,9 +1234,9 @@ struct LambdaTableProducer { if (tr1.sign() != tr2.sign()) { continue; } - + if (cFakeDaughterMaxDEta > 0.f && - std::abs(tr1.eta() - tr2.eta()) > cFakeDaughterMaxDEta) { + std::abs(tr1.eta() - tr2.eta()) > cFakeDaughterMaxDEta) { continue; } if (cFakeDaughterMaxOpenAngle > 0.f) { @@ -2031,20 +2030,16 @@ struct LambdaSpinPolarization { analyzePairsStandard(lambdaTracks, lambdaTracks); analyzePairsStandard(antiLambdaTracks, antiLambdaTracks); - // --------------------------------------------------------------------------- - // US-LS combinatorial background pairs - // --------------------------------------------------------------------------- analyzePairsStandard(lambdaTracks, fakeAntiLambdaTracks); analyzePairsStandard(fakeLambdaTracks, antiLambdaTracks); analyzePairsStandard(antiLambdaTracks, fakeLambdaTracks); analyzePairsStandard(fakeAntiLambdaTracks, lambdaTracks); - analyzePairsStandard(lambdaTracks, fakeLambdaTracks); - analyzePairsStandard(fakeLambdaTracks, lambdaTracks); - - analyzePairsStandard(antiLambdaTracks, fakeAntiLambdaTracks); - analyzePairsStandard(fakeAntiLambdaTracks, antiLambdaTracks); + analyzePairsStandard(lambdaTracks, fakeLambdaTracks); + analyzePairsStandard(fakeLambdaTracks, lambdaTracks); + analyzePairsStandard(antiLambdaTracks, fakeAntiLambdaTracks); + analyzePairsStandard(fakeAntiLambdaTracks, antiLambdaTracks); } PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, "Process for Data and MCReco", true); @@ -2102,10 +2097,9 @@ struct LambdaSpinPolarization { analyzePairsME(antiLambdaTracks_col1, fakeLambdaTracks_col2); analyzePairsME(fakeAntiLambdaTracks_col2, lambdaTracks_col1); - analyzePairsME(lambdaTracks_col1, fakeLambdaTracks_col2); + analyzePairsME(lambdaTracks_col1, fakeLambdaTracks_col2); analyzePairsME(fakeLambdaTracks_col2, lambdaTracks_col1); - - analyzePairsME(antiLambdaTracks_col1, fakeAntiLambdaTracks_col2); + analyzePairsME(antiLambdaTracks_col1, fakeAntiLambdaTracks_col2); analyzePairsME(fakeAntiLambdaTracks_col2, antiLambdaTracks_col1); } } From a53d9cf3fa736832333777811df3d09e0bac5b51 Mon Sep 17 00:00:00 2001 From: Subhadeep Roy Date: Fri, 10 Apr 2026 11:17:25 +0530 Subject: [PATCH 3/4] Fix comment formatting and clean up code --- .../TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx index 6141a916c2d..5aada7cac6c 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx @@ -450,7 +450,7 @@ struct LambdaTableProducer { histos.add("Tracks/h1f_lambda_pt_vs_invm", "p_{T} vs M_{#Lambda}", kTH2F, {axisV0Mass, axisV0Pt}); histos.add("Tracks/h1f_antilambda_pt_vs_invm", "p_{T} vs M_{#bar{#Lambda}}", kTH2F, {axisV0Mass, axisV0Pt}); - //histogram for LS subsampling fraction (nKept/nTotal per event) + // histogram for LS subsampling fraction (nKept/nTotal per event) histos.add("Tracks/h1f_fake_subsample_fraction", "LS subsample fraction (nKept/nTotal)", kTH1F, {{110, 0., 1.1, "f_{subsample}"}}); // QA Lambda @@ -1236,7 +1236,7 @@ struct LambdaTableProducer { } if (cFakeDaughterMaxDEta > 0.f && - std::abs(tr1.eta() - tr2.eta()) > cFakeDaughterMaxDEta) { + std::abs(tr1.eta() - tr2.eta()) > cFakeDaughterMaxDEta) { continue; } if (cFakeDaughterMaxOpenAngle > 0.f) { From 5a70c011422d2621edf08135ddab61e967ffeb1c Mon Sep 17 00:00:00 2001 From: Subhadeep Roy Date: Fri, 10 Apr 2026 12:07:12 +0530 Subject: [PATCH 4/4] Refactor variable names and apply static_cast --- .../Tasks/lambdaSpinPolarization.cxx | 140 +++++++++--------- 1 file changed, 70 insertions(+), 70 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx index 5aada7cac6c..d079b483981 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx @@ -1064,7 +1064,7 @@ struct LambdaTableProducer { if (!selTrackBasic(track)) return false; if (cApplyFakeDcaCuts) { - float minDca = isProton ? (float)cMinDcaProtonToPV : (float)cMinDcaPionToPV; + float minDca = isProton ? static_cast(cMinDcaProtonToPV) : static_cast(cMinDcaPionToPV); if (std::abs(track.dcaXY()) < minDca) return false; } @@ -1822,14 +1822,14 @@ struct LambdaSpinPolarization { static constexpr std::string_view SubDirHistUS[] = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; static constexpr std::string_view SubDirHistBkg[] = {"LaPFkLaM", "LaMLFkLaP", "LaPFkLaP", "LaMLFkLaM"}; - constexpr bool isBkg = (part_pair == kLambdaFakeAntiLambda || part_pair == kAntiLambdaFakeLambda || part_pair == kLambdaFakeLambda || part_pair == kAntiLambdaFakeAntiLambda); + constexpr bool IsBkg = (part_pair == kLambdaFakeAntiLambda || part_pair == kAntiLambdaFakeLambda || part_pair == kLambdaFakeLambda || part_pair == kAntiLambdaFakeAntiLambda); // Fill pair invariant mass histogram - if constexpr (!isBkg) { + if constexpr (!IsBkg) { histos.fill(HIST("Reco/h2f_n2_mass_") + HIST(SubDirHistUS[part_pair]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); } else { - constexpr int bkgIdx = (int)part_pair - 4; - histos.fill(HIST("RecoBkg/h2f_n2_mass_") + HIST(SubDirHistBkg[bkgIdx]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); + constexpr int BkgIdx = static_cast(part_pair) - 4; + histos.fill(HIST("RecoBkg/h2f_n2_mass_") + HIST(SubDirHistBkg[BkgIdx]), p1.mass(), p2.mass(), p1.pt(), p2.pt()); } float drap = p1.rap() - p2.rap(); @@ -1842,71 +1842,71 @@ struct LambdaSpinPolarization { std::array pr2 = {p2.prPx(), p2.prPy(), p2.prPz(), MassProton}; if (cDoAtlasMethod) { - std::array l1_atlas = l1; - std::array l2_atlas = l2; - std::array pr1_atlas = pr1; - std::array pr2_atlas = pr2; + std::array l1Atlas = l1; + std::array l2Atlas = l2; + std::array pr1Atlas = pr1; + std::array pr2Atlas = pr2; std::array llpair = { - l1_atlas[0] + l2_atlas[0], - l1_atlas[1] + l2_atlas[1], - l1_atlas[2] + l2_atlas[2], - l1_atlas[3] + l2_atlas[3]}; + l1Atlas[0] + l2Atlas[0], + l1Atlas[1] + l2Atlas[1], + l1Atlas[2] + l2Atlas[2], + l1Atlas[3] + l2Atlas[3]}; std::array vPair; getBoostVector(llpair, vPair, cInvBoostFlag); - boost(l1_atlas, vPair); - boost(l2_atlas, vPair); - boost(pr1_atlas, vPair); - boost(pr2_atlas, vPair); - - std::array v1_pair, v2_pair; - getBoostVector(l1_atlas, v1_pair, cInvBoostFlag); - getBoostVector(l2_atlas, v2_pair, cInvBoostFlag); - boost(pr1_atlas, v1_pair); - boost(pr2_atlas, v2_pair); - - std::array pr1tv_atlas = {pr1_atlas[0], pr1_atlas[1], pr1_atlas[2]}; - std::array pr2tv_atlas = {pr2_atlas[0], pr2_atlas[1], pr2_atlas[2]}; + boost(l1Atlas, vPair); + boost(l2Atlas, vPair); + boost(pr1Atlas, vPair); + boost(pr2Atlas, vPair); + + std::array v1Pair, v2Pair; + getBoostVector(l1Atlas, v1Pair, cInvBoostFlag); + getBoostVector(l2Atlas, v2Pair, cInvBoostFlag); + boost(pr1Atlas, v1Pair); + boost(pr2Atlas, v2Pair); + + std::array pr1tvAtlas = {pr1Atlas[0], pr1Atlas[1], pr1Atlas[2]}; + std::array pr2tvAtlas = {pr2Atlas[0], pr2Atlas[1], pr2Atlas[2]}; float ctheta = - RecoDecay::dotProd(pr1tv_atlas, pr2tv_atlas) / - (RecoDecay::sqrtSumOfSquares(pr1tv_atlas[0], pr1tv_atlas[1], pr1tv_atlas[2]) * - RecoDecay::sqrtSumOfSquares(pr2tv_atlas[0], pr2tv_atlas[1], pr2tv_atlas[2])); + RecoDecay::dotProd(pr1tvAtlas, pr2tvAtlas) / + (RecoDecay::sqrtSumOfSquares(pr1tvAtlas[0], pr1tvAtlas[1], pr1tvAtlas[2]) * + RecoDecay::sqrtSumOfSquares(pr2tvAtlas[0], pr2tvAtlas[1], pr2tvAtlas[2])); - if constexpr (!isBkg) { + if constexpr (!IsBkg) { histos.fill(HIST("RecoCorr/h2f_n2_ctheta_") + HIST(SubDirHistUS[part_pair]), cent, drap, dphi, ctheta); histos.fill(HIST("RecoCorr/h2f_n2_dltaR_") + HIST(SubDirHistUS[part_pair]), cent, dR, ctheta); } else { - constexpr int bkgIdx = (int)part_pair - 4; - histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[bkgIdx]), cent, drap, dphi, ctheta); - histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[bkgIdx]), cent, dR, ctheta); + constexpr int BkgIdx = static_cast(part_pair) - 4; + histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[BkgIdx]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[BkgIdx]), cent, dR, ctheta); } } if (cDoStarMethod) { - std::array pr1_star = pr1; - std::array pr2_star = pr2; + std::array pr1Star = pr1; + std::array pr2Star = pr2; - std::array v1_lab, v2_lab; - getBoostVector(l1, v1_lab, cInvBoostFlag); - getBoostVector(l2, v2_lab, cInvBoostFlag); - boost(pr1_star, v1_lab); - boost(pr2_star, v2_lab); + std::array v1Lab, v2Lab; + getBoostVector(l1, v1Lab, cInvBoostFlag); + getBoostVector(l2, v2Lab, cInvBoostFlag); + boost(pr1Star, v1Lab); + boost(pr2Star, v2Lab); - std::array pr1tv_star = {pr1_star[0], pr1_star[1], pr1_star[2]}; - std::array pr2tv_star = {pr2_star[0], pr2_star[1], pr2_star[2]}; + std::array pr1tvStar = {pr1Star[0], pr1Star[1], pr1Star[2]}; + std::array pr2tvStar = {pr2Star[0], pr2Star[1], pr2Star[2]}; float ctheta = - RecoDecay::dotProd(pr1tv_star, pr2tv_star) / - (RecoDecay::sqrtSumOfSquares(pr1tv_star[0], pr1tv_star[1], pr1tv_star[2]) * - RecoDecay::sqrtSumOfSquares(pr2tv_star[0], pr2tv_star[1], pr2tv_star[2])); + RecoDecay::dotProd(pr1tvStar, pr2tvStar) / + (RecoDecay::sqrtSumOfSquares(pr1tvStar[0], pr1tvStar[1], pr1tvStar[2]) * + RecoDecay::sqrtSumOfSquares(pr2tvStar[0], pr2tvStar[1], pr2tvStar[2])); - if constexpr (!isBkg) { + if constexpr (!IsBkg) { histos.fill(HIST("RecoCorr/h2f_n2_ctheta_") + HIST(SubDirHistUS[part_pair]), cent, drap, dphi, ctheta); histos.fill(HIST("RecoCorr/h2f_n2_dltaR_") + HIST(SubDirHistUS[part_pair]), cent, dR, ctheta); } else { - constexpr int bkgIdx = (int)part_pair - 4; - histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[bkgIdx]), cent, drap, dphi, ctheta); - histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[bkgIdx]), cent, dR, ctheta); + constexpr int BkgIdx = static_cast(part_pair) - 4; + histos.fill(HIST("RecoCorrBkg/h2f_n2_ctheta_") + HIST(SubDirHistBkg[BkgIdx]), cent, drap, dphi, ctheta); + histos.fill(HIST("RecoCorrBkg/h2f_n2_dltaR_") + HIST(SubDirHistBkg[BkgIdx]), cent, dR, ctheta); } } } @@ -1949,7 +1949,7 @@ struct LambdaSpinPolarization { decltype(&*trks_2.begin()) currentPartner = nullptr; - for (auto const* trk_2 : candidates) { + for (auto const& trk_2 : candidates) { if constexpr (samelambda) { if (trk_1.index() == trk_2->index()) @@ -2071,36 +2071,36 @@ struct LambdaSpinPolarization { histos.fill(HIST("QA/ME/hPoolCentVz"), col1.cent(), col1.posZ()); // Lambda slices - auto lambdaTracks_col1 = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); - auto lambdaTracks_col2 = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + auto lambdaTracksCol1 = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); + auto lambdaTracksCol2 = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); // Anti-lambda slices - auto antiLambdaTracks_col1 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); - auto antiLambdaTracks_col2 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + auto antiLambdaTracksCol1 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); + auto antiLambdaTracksCol2 = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); // Fake lambda slices (LS background) - auto fakeLambdaTracks_col2 = partFakeLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); - auto fakeAntiLambdaTracks_col2 = partFakeAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + auto fakeLambdaTracksCol2 = partFakeLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); + auto fakeAntiLambdaTracksCol2 = partFakeAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); // QA: multiplicities - histos.fill(HIST("QA/ME/hLambdaMultVsCent"), col1.cent(), lambdaTracks_col1.size()); - histos.fill(HIST("QA/ME/hAntiLambdaMultVsCent"), col1.cent(), antiLambdaTracks_col1.size()); + histos.fill(HIST("QA/ME/hLambdaMultVsCent"), col1.cent(), lambdaTracksCol1.size()); + histos.fill(HIST("QA/ME/hAntiLambdaMultVsCent"), col1.cent(), antiLambdaTracksCol1.size()); - analyzePairsME(lambdaTracks_col1, antiLambdaTracks_col2); - analyzePairsME(antiLambdaTracks_col1, lambdaTracks_col2); - analyzePairsME(lambdaTracks_col1, lambdaTracks_col2); - analyzePairsME(antiLambdaTracks_col1, antiLambdaTracks_col2); + analyzePairsME(lambdaTracksCol1, antiLambdaTracksCol2); + analyzePairsME(antiLambdaTracksCol1, lambdaTracksCol2); + analyzePairsME(lambdaTracksCol1, lambdaTracksCol2); + analyzePairsME(antiLambdaTracksCol1, antiLambdaTracksCol2); - analyzePairsME(lambdaTracks_col1, fakeAntiLambdaTracks_col2); - analyzePairsME(fakeLambdaTracks_col2, antiLambdaTracks_col1); + analyzePairsME(lambdaTracksCol1, fakeAntiLambdaTracksCol2); + analyzePairsME(fakeLambdaTracksCol2, antiLambdaTracksCol1); - analyzePairsME(antiLambdaTracks_col1, fakeLambdaTracks_col2); - analyzePairsME(fakeAntiLambdaTracks_col2, lambdaTracks_col1); + analyzePairsME(antiLambdaTracksCol1, fakeLambdaTracksCol2); + analyzePairsME(fakeAntiLambdaTracksCol2, lambdaTracksCol1); - analyzePairsME(lambdaTracks_col1, fakeLambdaTracks_col2); - analyzePairsME(fakeLambdaTracks_col2, lambdaTracks_col1); - analyzePairsME(antiLambdaTracks_col1, fakeAntiLambdaTracks_col2); - analyzePairsME(fakeAntiLambdaTracks_col2, antiLambdaTracks_col1); + analyzePairsME(lambdaTracksCol1, fakeLambdaTracksCol2); + analyzePairsME(fakeLambdaTracksCol2, lambdaTracksCol1); + analyzePairsME(antiLambdaTracksCol1, fakeAntiLambdaTracksCol2); + analyzePairsME(fakeAntiLambdaTracksCol2, antiLambdaTracksCol1); } }