From ae2d8ad28207a09cb8dbf4e86a28a3d895b57a98 Mon Sep 17 00:00:00 2001 From: Subhadeep Roy Date: Fri, 2 Oct 2026 22:33:57 +0530 Subject: [PATCH] Updated the lambdaSpinPolarization.cxx with the pair constructions for spin correlation Updated the same-event and mixed-event pair constructions --- .../Tasks/lambdaSpinPolarization.cxx | 2305 ++++++++--------- 1 file changed, 1082 insertions(+), 1223 deletions(-) diff --git a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx index f5d26c6fbaa..67a58ab98ec 100644 --- a/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx +++ b/PWGCF/TwoParticleCorrelations/Tasks/lambdaSpinPolarization.cxx @@ -43,18 +43,27 @@ #include #include +#include #include +#include +#include #include +#include #include #include #include +#include #include +#include #include #include #include +#include +#include #include #include +#include #include using namespace o2; @@ -70,7 +79,6 @@ namespace lambdacollision DECLARE_SOA_COLUMN(Cent, cent, float); DECLARE_SOA_COLUMN(Mult, mult, float); DECLARE_SOA_COLUMN(TimeStamp, timeStamp, uint64_t); -// DECALRE_SOA_ } // namespace lambdacollision DECLARE_SOA_TABLE(LambdaCollisions, "AOD", "LAMBDACOLS", o2::soa::Index<>, @@ -172,7 +180,6 @@ using LambdaMixEventMcGenCollision = LambdaMixEventMcGenCollisions::iterator; namespace lambdamixeventtracks { -// DECLARE_SOA_INDEX_COLUMN(LambdaMixEventCollision, lambdaMixEventCollision); DECLARE_SOA_COLUMN(LambdaMixEventCollisionIdx, lambdaMixEventCollisionIdx, int); DECLARE_SOA_COLUMN(LambdaMixEventTrackIdx, lambdaMixEventTrackIdx, int); DECLARE_SOA_COLUMN(LambdaMixEventTimeStamp, lambdaMixEventTimeStamp, uint64_t); @@ -180,7 +187,6 @@ DECLARE_SOA_COLUMN(LambdaMixEventTimeStamp, lambdaMixEventTimeStamp, uint64_t); DECLARE_SOA_TABLE(LambdaMixEventTracks, "AOD", "LAMBDAMIXTRKS", o2::soa::Index<>, - // lambdamixeventtracks::LambdaMixEventCollisionId, lambdamixeventtracks::LambdaMixEventCollisionIdx, lambdamixeventtracks::LambdaMixEventTrackIdx, lambdatrack::Px, lambdatrack::Py, lambdatrack::Pz, lambdatrack::Mass, @@ -207,6 +213,21 @@ enum CollisionLabels { kTotColBeforeHasMcCollision = 1, kTotCol, kPassSelCol }; +enum EventCutFlow { kEvAll = 1, + kEvTrigger, + kEvTvx, + kEvTFBorder, + kEvItsRofBorder, + kEvItsTpcVtx, + kEvNoSameBunchPileup, + kEvGoodZvtxFT0vsPV, + kEvGoodItsLayers, + kEvCentrality, + kEvVz, + kEvOneV0, + kEvTwoV0, + kNEventCutFlow }; + enum TrackLabels { kTracksBeforeHasMcParticle = 1, kAllV0Tracks, @@ -245,23 +266,32 @@ enum RunType { kRun3 = 0, enum ParticleType { kLambda = 0, kAntiLambda }; -enum ParticlePairType { - - kLambdaAntiLambda = 0, - kAntiLambdaLambda = 1, - kLambdaLambda = 2, - kAntiLambdaAntiLambda = 3, - - kLambdaSBAntiLambda = 4, - kAntiLambdaSBLambda = 5, - kLambdaSBLambda = 6, - kAntiLambdaSBAntiLambda = 7, - - kLambdaSBSBAntiLambda = 8, - kAntiLambdaSBSBLambda = 9, - kLambdaSBSBLambda = 10, - kAntiLambdaSBSBAntiLambda = 11 -}; +enum ParticlePairType { kLambdaAntiLambda = 0, + kAntiLambdaLambda, + kLambdaLambda, + kAntiLambdaAntiLambda }; + +enum PairTier { kSameEvent = 0, + kSameEventInMixing, + kMixedEvent, + kMcGenSameEvent, + kMcGenSameEventInMixing, + kMcGenMixedEvent }; + +enum PairEventCount { kPeAll = 1, + kPeOneCand, + kPeTwoCand, + kPeUnlikeSign, + kPeLambdaLambda, + kPeAntiLambdaAntiLambda, + kNPairEventCount }; + +enum MEEventStatus { kMEEventSeen = 1, + kMEEventOutsidePools, + kMEEventWarmUp, + kMEEventMixed, + kMEEventDonated, + kNMEEventStatus }; enum ShareDauLambda { kUniqueLambda = 0, kLambdaShareDau }; @@ -284,9 +314,10 @@ enum PrmScdPairType { kPP = 0, kSP, kSS }; -// Number of daughters of the two-body decay Lambda -> p + pi static constexpr std::size_t NDaughtersTwoBody = 2; +static constexpr std::size_t NCandidatesForPair = 2; + struct LambdaTableProducer { Produces lambdaCollisionTable; @@ -294,8 +325,7 @@ struct LambdaTableProducer { Produces lambdaMCGenCollisionTable; Produces lambdaMCGenTrackTable; - Configurable cCentEstimator{"cCentEstimator", 0, - "Centrality Estimator: 0=FT0M, 1=FT0C"}; + Configurable cCentEstimator{"cCentEstimator", 0, "Centrality Estimator: 0=FT0M, 1=FT0C"}; Configurable cMinZVtx{"cMinZVtx", -10.0, "Min VtxZ (cm)"}; Configurable cMaxZVtx{"cMaxZVtx", 10.0, "Max VtxZ (cm)"}; Configurable cMinMult{"cMinMult", 0.0, "Min centrality percentile"}; @@ -303,130 +333,83 @@ struct LambdaTableProducer { Configurable cSel8Trig{"cSel8Trig", true, "Sel8 (T0A+T0C) Run3"}; Configurable cInt7Trig{"cInt7Trig", false, "kINT7 MB Trigger"}; Configurable cSel7Trig{"cSel7Trig", false, "Sel7 (V0A+V0C) Run2"}; - Configurable cTriggerTvxSel{"cTriggerTvxSel", false, - "TVX Trigger Selection"}; - Configurable cTFBorder{"cTFBorder", false, - "Timeframe Border Selection"}; - Configurable cNoItsROBorder{"cNoItsROBorder", false, - "No ITSRO Border Cut"}; - Configurable cItsTpcVtx{"cItsTpcVtx", false, - "ITS+TPC Vertex Selection"}; + Configurable cTriggerTvxSel{"cTriggerTvxSel", false, "TVX Trigger Selection"}; + Configurable cTFBorder{"cTFBorder", false, "Timeframe Border Selection"}; + Configurable cNoItsROBorder{"cNoItsROBorder", false, "No ITSRO Border Cut"}; + Configurable cItsTpcVtx{"cItsTpcVtx", false, "ITS+TPC Vertex Selection"}; Configurable cPileupReject{"cPileupReject", false, "Pileup rejection"}; - Configurable cZVtxTimeDiff{"cZVtxTimeDiff", false, - "z-vtx time diff selection"}; - Configurable cIsGoodITSLayers{"cIsGoodITSLayers", false, - "Good ITS Layers All"}; + Configurable cZVtxTimeDiff{"cZVtxTimeDiff", false, "z-vtx time diff selection"}; + Configurable cIsGoodITSLayers{"cIsGoodITSLayers", false, "Good ITS Layers All"}; - // Tracks Configurable cTrackMinPt{"cTrackMinPt", 0.15, "p_{T} minimum"}; Configurable cTrackMaxPt{"cTrackMaxPt", 999.0, "p_{T} maximum"}; Configurable cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"}; - Configurable cMinTpcCrossedRows{"cMinTpcCrossedRows", 70, - "TPC Min Crossed Rows"}; - Configurable cMinTpcCROverCls{"cMinTpcCROverCls", 0.8, - "TPC Min CR/Findable Cls"}; - Configurable cMaxTpcSharedClusters{"cMaxTpcSharedClusters", 0.4, - "TPC Max Shared Clusters"}; - Configurable cMaxChi2Tpc{"cMaxChi2Tpc", 4, "Max TPC Chi2/ndf"}; - Configurable cTpcNsigmaCut{"cTpcNsigmaCut", 3.0, - "TPC nSigma PID cut"}; - Configurable cRemoveAmbiguousTracks{"cRemoveAmbiguousTracks", false, - "Remove Ambiguous Tracks"}; - - Configurable cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.02, - "Min proton DCA to PV (cm)"}; - Configurable cMinDcaPionToPV{"cMinDcaPionToPV", 0.06, - "Min pion DCA to PV (cm)"}; - Configurable cMinV0DcaDaughters{"cMinV0DcaDaughters", 0., - "Min DCA between V0 daughters"}; - Configurable cMaxV0DcaDaughters{"cMaxV0DcaDaughters", 1., - "Max DCA between V0 daughters"}; + Configurable cMinTpcCrossedRows{"cMinTpcCrossedRows", 70, "TPC Min Crossed Rows"}; + Configurable cMinTpcCROverCls{"cMinTpcCROverCls", -999, "TPC Min CR/Findable Cls"}; + Configurable cMaxTpcSharedClusters{"cMaxTpcSharedClusters", 999, "TPC Max Shared Clusters"}; + Configurable cMaxChi2Tpc{"cMaxChi2Tpc", 999, "Max TPC Chi2/ndf"}; + Configurable cTpcNsigmaCut{"cTpcNsigmaCut", 5.0, "TPC nSigma PID cut"}; + Configurable cRemoveAmbiguousTracks{"cRemoveAmbiguousTracks", false, "Remove Ambiguous Tracks"}; + + Configurable cMinDcaProtonToPV{"cMinDcaProtonToPV", 0.05, "Min proton DCA to PV (cm)"}; + Configurable cMinDcaPionToPV{"cMinDcaPionToPV", 0.1, "Min pion DCA to PV (cm)"}; + Configurable cMinV0DcaDaughters{"cMinV0DcaDaughters", 0., "Min DCA between V0 daughters"}; + Configurable cMaxV0DcaDaughters{"cMaxV0DcaDaughters", 1.4, "Max DCA between V0 daughters"}; Configurable cMinDcaV0ToPV{"cMinDcaV0ToPV", 0.0, "Min DCA V0 to PV"}; - Configurable cMaxDcaV0ToPV{"cMaxDcaV0ToPV", 999.0, - "Max DCA V0 to PV"}; - Configurable cMinV0TransRadius{"cMinV0TransRadius", 0.5, - "Min V0 decay radius (cm)"}; - Configurable cMaxV0TransRadius{"cMaxV0TransRadius", 999.0, - "Max V0 decay radius (cm)"}; + Configurable cMaxDcaV0ToPV{"cMaxDcaV0ToPV", 0.1, "Max DCA V0 to PV"}; + Configurable cMinV0TransRadius{"cMinV0TransRadius", 0.5, "Min V0 decay radius (cm)"}; + Configurable cMaxV0TransRadius{"cMaxV0TransRadius", 999.0, "Max V0 decay radius (cm)"}; Configurable cMinV0CTau{"cMinV0CTau", 0.0, "Min cTau (cm)"}; Configurable cMaxV0CTau{"cMaxV0CTau", 30.0, "Max cTau (cm)"}; Configurable cMinV0CosPA{"cMinV0CosPA", 0.995, "Min V0 cos(PA)"}; - Configurable cKshortRejMassWindow{"cKshortRejMassWindow", 0.01, - "K0s mass rejection window"}; - Configurable cKshortRejFlag{"cKshortRejFlag", true, - "K0s mass rejection flag"}; - - // V0s kinmatic acceptance - Configurable cMinV0Mass{"cMinV0Mass", 1.10, "V0 Mass Min"}; - Configurable cMaxV0Mass{"cMaxV0Mass", 1.12, "V0 Mass Max"}; - Configurable cMinV0Pt{"cMinV0Pt", 0.8, "Minimum V0 pT"}; - Configurable cMaxV0Pt{"cMaxV0Pt", 4.2, "Minimum V0 pT"}; + Configurable cKshortRejMassWindow{"cKshortRejMassWindow", 0.001, "K0s mass rejection window"}; + Configurable cKshortRejFlag{"cKshortRejFlag", true, "K0s mass rejection flag"}; + + Configurable cMinV0Mass{"cMinV0Mass", 1.09, "V0 Mass Min"}; + Configurable cMaxV0Mass{"cMaxV0Mass", 1.14, "V0 Mass Max"}; + Configurable cMinV0Pt{"cMinV0Pt", 0.6, "Minimum V0 pT"}; + Configurable cMaxV0Pt{"cMaxV0Pt", 3.0, "Maximum V0 pT"}; Configurable cMaxV0Rap{"cMaxV0Rap", 0.5, "|rap| cut"}; Configurable cDoEtaAnalysis{"cDoEtaAnalysis", false, "Do Eta Analysis"}; - Configurable cV0TypeSelFlag{"cV0TypeSelFlag", false, - "V0 Type Selection Flag"}; - Configurable cV0TypeSelection{"cV0TypeSelection", 1, - "V0 Type Selection"}; + Configurable cV0TypeSelFlag{"cV0TypeSelFlag", true, "V0 Type Selection Flag"}; + Configurable cV0TypeSelection{"cV0TypeSelection", 1, "V0 Type Selection"}; - // V0s MC Configurable cHasMcFlag{"cHasMcFlag", true, "Has Mc Tag"}; - Configurable cSelectTrueLambda{"cSelectTrueLambda", true, - "Select True Lambda"}; - Configurable cSelMCPSV0{"cSelMCPSV0", true, - "Select Primary/Secondary V0"}; - Configurable cCheckRecoDauFlag{"cCheckRecoDauFlag", true, - "Check for reco daughter PID"}; - Configurable cGenPrimaryLambda{"cGenPrimaryLambda", true, - "Primary Generated Lambda"}; - Configurable cGenSecondaryLambda{"cGenSecondaryLambda", false, - "Secondary Generated Lambda"}; - Configurable cGenDecayChannel{"cGenDecayChannel", true, - "Gen Level Decay Channel Flag"}; + Configurable cSelectTrueLambda{"cSelectTrueLambda", false, "Select True Lambda"}; + Configurable cSelMCPSV0{"cSelMCPSV0", false, "Select Primary/Secondary V0"}; + Configurable cCheckRecoDauFlag{"cCheckRecoDauFlag", true, "Check for reco daughter PID"}; + Configurable cGenPrimaryLambda{"cGenPrimaryLambda", true, "Primary Generated Lambda"}; + Configurable cGenSecondaryLambda{"cGenSecondaryLambda", false, "Secondary Generated Lambda"}; + Configurable cGenDecayChannel{"cGenDecayChannel", false, "Gen Level Decay Channel Flag"}; Configurable cRecoMomResoFlag{"cRecoMomResoFlag", false, "Check effect of momentum space smearing on balance function"}; - // Efficiency Correction - Configurable cCorrectionFlag{"cCorrectionFlag", false, - "Correction Flag"}; - Configurable cGetEffFact{"cGetEffFact", false, - "Get Efficiency Factor Flag"}; - Configurable cGetPrimFrac{"cGetPrimFrac", false, - "Get Primary Fraction Flag"}; - Configurable cCorrFactHist{"cCorrFactHist", 0, - "Efficiency Factor Histogram"}; - Configurable cPrimFracHist{"cPrimFracHist", 0, - "Primary Fraction Histogram"}; - - // CCDB - Configurable cUrlCCDB{"cUrlCCDB", "http://ccdb-test.cern.ch:8080", "url of ccdb"}; - Configurable cPathCCDB{"cPathCCDB", "Users/y/ypatley/lambda_corr_fact", "Path for ccdb-object"}; - - // Initialize CCDB Service + Configurable cCorrectionFlag{"cCorrectionFlag", false, "Correction Flag"}; + Configurable cGetEffFact{"cGetEffFact", true, "Get Efficiency Factor Flag"}; + Configurable cGetPrimFrac{"cGetPrimFrac", false, "Get Primary Fraction Flag"}; + Configurable cCorrFactHist{"cCorrFactHist", 0, "Efficiency Factor Histogram"}; + Configurable cPrimFracHist{"cPrimFracHist", 0, "Primary Fraction Histogram"}; + + Configurable cUrlCCDB{"cUrlCCDB", "http://alice-ccdb.cern.ch", "url of ccdb"}; + Configurable cPathCCDB{"cPathCCDB", "Users/y/ypatley/SpinCorr/RecoEfficiency", "Path for ccdb-object"}; + Service ccdb{}; - // Histogram Registry. - HistogramRegistry histos{ - "histos", - {}, - OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - // initialize corr_factor objects std::vector> vCorrFactStrings = { {"hEffVsPtCentLambda", "hEffVsPtCentAntiLambda"}, {"hEffVsPtYCentLambda", "hEffVsPtYCentAntiLambda"}, {"hEffVsPtEtaCentLambda", "hEffVsPtEtaCentAntiLambda"}}; - // initialize corr_factor objects std::vector> vPrimFracStrings = { {"hPrimFracVsPtCentLambda", "hPrimFracVsPtCentAntiLambda"}, {"hPrimFracVsPtYCentLambda", "hPrimFracVsPtYCentAntiLambda"}, {"hPrimFracVsPtEtaCentLambda", "hPrimFracVsPtEtaCentAntiLambda"}}; - // Initialize Global Variables float cent = 0., mult = 0.; void init(InitContext const&) { - // Set CCDB url ccdb->setURL(cUrlCCDB.value); ccdb->setCaching(true); @@ -437,7 +420,7 @@ struct LambdaTableProducer { const AxisSpec axisVz(220, -11, 11, "V_{z} (cm)"); const AxisSpec axisPID(8000, -4000, 4000, "PdgCode"); - const AxisSpec axisV0Mass(120, 1.08, 1.20, "M_{p#pi} (GeV/#it{c}^{2})"); + const AxisSpec axisV0Mass(200, 1.08, 1.18, "M_{p#pi} (GeV/#it{c}^{2})"); const AxisSpec axisV0Pt(100., 0., 10., "p_{T} (GeV/#it{c})"); const AxisSpec axisV0Rap(48, -1.2, 1.2, "y"); const AxisSpec axisV0Eta(48, -1.2, 1.2, "#eta"); @@ -458,102 +441,77 @@ struct LambdaTableProducer { const AxisSpec axisNsigma(401, -10.025, 10.025, "n#sigma"); const AxisSpec axisdEdx(360, 20, 200, "#frac{dE}{dx}"); - // Create Histograms. - // Event histograms - histos.add("Events/h1f_collisions_info", "# of Collisions", kTH1F, - {axisCols}); - histos.add("Events/h1f_collision_posZ", "V_{z}-distribution", kTH1F, - {axisVz}); + histos.add("Events/h1f_collisions_info", "# of Collisions", kTH1F, {axisCols}); + histos.add("Events/h1f_collision_posZ", "V_{z}-distribution", kTH1F, {axisVz}); + histos.add("Events/h1f_collision_cent", "Centrality of the selected collisions", kTH1F, {axisCent}); + auto hCutFlow = histos.add("Events/hEventCutFlow", "event cut flow;;collisions", kTH1D, {{kNEventCutFlow - 1, 0.5, static_cast(kNEventCutFlow) - 0.5}}); + hCutFlow->GetXaxis()->SetBinLabel(kEvAll, "all"); + hCutFlow->GetXaxis()->SetBinLabel(kEvTrigger, "trigger (sel8)"); + hCutFlow->GetXaxis()->SetBinLabel(kEvTvx, "TVX"); + hCutFlow->GetXaxis()->SetBinLabel(kEvTFBorder, "no TF border"); + hCutFlow->GetXaxis()->SetBinLabel(kEvItsRofBorder, "no ITS ROF border"); + hCutFlow->GetXaxis()->SetBinLabel(kEvItsTpcVtx, "ITS-TPC vertex"); + hCutFlow->GetXaxis()->SetBinLabel(kEvNoSameBunchPileup, "no same-bunch pileup"); + hCutFlow->GetXaxis()->SetBinLabel(kEvGoodZvtxFT0vsPV, "good z_{vtx} FT0 vs PV"); + hCutFlow->GetXaxis()->SetBinLabel(kEvGoodItsLayers, "good ITS layers"); + hCutFlow->GetXaxis()->SetBinLabel(kEvCentrality, "centrality range"); + hCutFlow->GetXaxis()->SetBinLabel(kEvVz, "V_{z} range"); + hCutFlow->GetXaxis()->SetBinLabel(kEvOneV0, "#geq 1 selected #Lambda/#bar{#Lambda}"); + hCutFlow->GetXaxis()->SetBinLabel(kEvTwoV0, "#geq 2 selected #Lambda/#bar{#Lambda}"); histos.add("Tracks/h1f_tracks_info", "# of tracks", kTH1F, {axisTrks}); - histos.add("Tracks/h2f_armpod_before_sel", "Armenteros-Podolanski (before)", - kTH2F, {axisAlpha, axisQtarm}); - histos.add("Tracks/h2f_armpod_after_sel", "Armenteros-Podolanski (after)", - kTH2F, {axisAlpha, axisQtarm}); - 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}); - - histos.add("QA/Lambda/h2f_qt_vs_alpha", "Armenteros-Podolanski", kTH2F, - {axisAlpha, axisQtarm}); - histos.add("QA/Lambda/h1f_dca_V0_daughters", "DCA V0 daughters", kTH1F, - {axisDcaDau}); - histos.add("QA/Lambda/h1f_dca_pos_to_PV", "DCA pos-prong to PV", kTH1F, - {axisDcaProngPV}); - histos.add("QA/Lambda/h1f_dca_neg_to_PV", "DCA neg-prong to PV", kTH1F, - {axisDcaProngPV}); - histos.add("QA/Lambda/h1f_dca_V0_to_PV", "DCA V0 to PV", kTH1F, - {axisDcaV0PV}); - histos.add("QA/Lambda/h1f_V0_cospa", "cos(#theta_{PA})", kTH1F, - {axisCosPA}); - histos.add("QA/Lambda/h1f_V0_radius", "V0 decay radius", kTH1F, - {axisRadius}); + histos.add("Tracks/h2f_armpod_before_sel", "Armenteros-Podolanski (before)", kTH2F, {axisAlpha, axisQtarm}); + histos.add("Tracks/h2f_armpod_after_sel", "Armenteros-Podolanski (after)", kTH2F, {axisAlpha, axisQtarm}); + 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}); + + histos.add("QA/Lambda/h2f_qt_vs_alpha", "Armenteros-Podolanski", kTH2F, {axisAlpha, axisQtarm}); + histos.add("QA/Lambda/h1f_dca_V0_daughters", "DCA V0 daughters", kTH1F, {axisDcaDau}); + histos.add("QA/Lambda/h1f_dca_pos_to_PV", "DCA pos-prong to PV", kTH1F, {axisDcaProngPV}); + histos.add("QA/Lambda/h1f_dca_neg_to_PV", "DCA neg-prong to PV", kTH1F, {axisDcaProngPV}); + histos.add("QA/Lambda/h1f_dca_V0_to_PV", "DCA V0 to PV", kTH1F, {axisDcaV0PV}); + histos.add("QA/Lambda/h1f_V0_cospa", "cos(#theta_{PA})", kTH1F, {axisCosPA}); + histos.add("QA/Lambda/h1f_V0_radius", "V0 decay radius", kTH1F, {axisRadius}); histos.add("QA/Lambda/h1f_V0_ctau", "c#tau", kTH1F, {axisCTau}); histos.add("QA/Lambda/h1f_V0_gctau", "#gammac#tau", kTH1F, {axisGCTau}); - histos.add("QA/Lambda/h1f_pos_prong_pt", "Pos-prong p_{T}", kTH1F, - {axisTrackPt}); - histos.add("QA/Lambda/h1f_neg_prong_pt", "Neg-prong p_{T}", kTH1F, - {axisTrackPt}); - histos.add("QA/Lambda/h1f_pos_prong_eta", "Pos-prong #eta", kTH1F, - {axisV0Eta}); - histos.add("QA/Lambda/h1f_neg_prong_eta", "Neg-prong #eta", kTH1F, - {axisV0Eta}); - histos.add("QA/Lambda/h1f_pos_prong_phi", "Pos-prong #phi", kTH1F, - {axisV0Phi}); - histos.add("QA/Lambda/h1f_neg_prong_phi", "Neg-prong #phi", kTH1F, - {axisV0Phi}); - histos.add("QA/Lambda/h2f_pos_prong_dcaXY_vs_pt", "DCA vs p_{T}", kTH2F, - {axisTrackPt, axisTrackDCA}); - histos.add("QA/Lambda/h2f_neg_prong_dcaXY_vs_pt", "DCA vs p_{T}", kTH2F, - {axisTrackPt, axisTrackDCA}); - histos.add("QA/Lambda/h2f_pos_prong_dEdx_vs_p", "TPC dE/dx pos", kTH2F, - {axisMomPID, axisdEdx}); - histos.add("QA/Lambda/h2f_neg_prong_dEdx_vs_p", "TPC dE/dx neg", kTH2F, - {axisMomPID, axisdEdx}); - histos.add("QA/Lambda/h2f_pos_prong_tpc_nsigma_pr_vs_p", - "TPC n#sigma_{p} pos", kTH2F, {axisMomPID, axisNsigma}); - histos.add("QA/Lambda/h2f_neg_prong_tpc_nsigma_pr_vs_p", - "TPC n#sigma_{p} neg", kTH2F, {axisMomPID, axisNsigma}); - histos.add("QA/Lambda/h2f_pos_prong_tpc_nsigma_pi_vs_p", - "TPC n#sigma_{#pi} pos", kTH2F, {axisMomPID, axisNsigma}); - histos.add("QA/Lambda/h2f_neg_prong_tpc_nsigma_pi_vs_p", - "TPC n#sigma_{#pi} neg", kTH2F, {axisMomPID, axisNsigma}); + histos.add("QA/Lambda/h1f_pos_prong_pt", "Pos-prong p_{T}", kTH1F, {axisTrackPt}); + histos.add("QA/Lambda/h1f_neg_prong_pt", "Neg-prong p_{T}", kTH1F, {axisTrackPt}); + histos.add("QA/Lambda/h1f_pos_prong_eta", "Pos-prong #eta", kTH1F, {axisV0Eta}); + histos.add("QA/Lambda/h1f_neg_prong_eta", "Neg-prong #eta", kTH1F, {axisV0Eta}); + histos.add("QA/Lambda/h1f_pos_prong_phi", "Pos-prong #phi", kTH1F, {axisV0Phi}); + histos.add("QA/Lambda/h1f_neg_prong_phi", "Neg-prong #phi", kTH1F, {axisV0Phi}); + histos.add("QA/Lambda/h2f_pos_prong_dcaXY_vs_pt", "DCA vs p_{T}", kTH2F, {axisTrackPt, axisTrackDCA}); + histos.add("QA/Lambda/h2f_neg_prong_dcaXY_vs_pt", "DCA vs p_{T}", kTH2F, {axisTrackPt, axisTrackDCA}); + histos.add("QA/Lambda/h2f_pos_prong_dEdx_vs_p", "TPC dE/dx pos", kTH2F, {axisMomPID, axisdEdx}); + histos.add("QA/Lambda/h2f_neg_prong_dEdx_vs_p", "TPC dE/dx neg", kTH2F, {axisMomPID, axisdEdx}); + histos.add("QA/Lambda/h2f_pos_prong_tpc_nsigma_pr_vs_p", "TPC n#sigma_{p} pos", kTH2F, {axisMomPID, axisNsigma}); + histos.add("QA/Lambda/h2f_neg_prong_tpc_nsigma_pr_vs_p", "TPC n#sigma_{p} neg", kTH2F, {axisMomPID, axisNsigma}); + histos.add("QA/Lambda/h2f_pos_prong_tpc_nsigma_pi_vs_p", "TPC n#sigma_{#pi} pos", kTH2F, {axisMomPID, axisNsigma}); + histos.add("QA/Lambda/h2f_neg_prong_tpc_nsigma_pi_vs_p", "TPC n#sigma_{#pi} neg", kTH2F, {axisMomPID, axisNsigma}); histos.add("McRec/Lambda/hPt", "p_{T}", kTH1F, {axisV0Pt}); histos.add("McRec/Lambda/hEta", "#eta", kTH1F, {axisV0Eta}); histos.add("McRec/Lambda/hRap", "y", kTH1F, {axisV0Rap}); histos.add("McRec/Lambda/hPhi", "#phi", kTH1F, {axisV0Phi}); - // QA Anti-Lambda histos.addClone("QA/Lambda/", "QA/AntiLambda/"); histos.addClone("McRec/Lambda/", "McRec/AntiLambda/"); if (doprocessMCRun3) { - histos.add("Tracks/h2f_tracks_pid_before_sel", "PIDs before sel", kTH2F, - {axisPID, axisV0Pt}); - histos.add("Tracks/h2f_tracks_pid_after_sel", "PIDs after sel", kTH2F, - {axisPID, axisV0Pt}); - histos.add("Tracks/h2f_lambda_mothers_pdg", "Lambda mothers", kTH2F, - {axisPID, axisV0Pt}); - - histos.add("McGen/h1f_collision_recgen", "RecGen collisions", kTH1F, - {axisMult}); - histos.add("McGen/h1f_collisions_info", "Collisions info", kTH1F, - {axisCols}); - histos.add("McGen/h2f_collision_posZ", "V_{z} rec vs gen", kTH2F, - {axisVz, axisVz}); - histos.add("McGen/h2f_collision_cent", "Centrality rec vs gen", kTH2F, - {axisCent, axisCent}); - histos.add("McGen/h1f_lambda_daughter_PDG", "Lambda dau PDG", kTH1F, - {axisPID}); - histos.add("McGen/h1f_antilambda_daughter_PDG", "AntiLambda dau PDG", - kTH1F, {axisPID}); + histos.add("Tracks/h2f_tracks_pid_before_sel", "PIDs before sel", kTH2F, {axisPID, axisV0Pt}); + histos.add("Tracks/h2f_tracks_pid_after_sel", "PIDs after sel", kTH2F, {axisPID, axisV0Pt}); + histos.add("Tracks/h2f_lambda_mothers_pdg", "Lambda mothers", kTH2F, {axisPID, axisV0Pt}); + + histos.add("McGen/h1f_collision_recgen", "RecGen collisions", kTH1F, {axisMult}); + histos.add("McGen/h1f_collisions_info", "Collisions info", kTH1F, {axisCols}); + histos.add("McGen/h2f_collision_posZ", "V_{z} rec vs gen", kTH2F, {axisVz, axisVz}); + histos.add("McGen/h2f_collision_cent", "Centrality rec vs gen", kTH2F, {axisCent, axisCent}); + histos.add("McGen/h1f_lambda_daughter_PDG", "Lambda dau PDG", kTH1F, {axisPID}); + histos.add("McGen/h1f_antilambda_daughter_PDG", "AntiLambda dau PDG", kTH1F, {axisPID}); histos.addClone("McRec/", "McGen/"); - histos.add("McGen/Lambda/Proton/hPt", "Proton p_{T}", kTH1F, - {axisTrackPt}); + histos.add("McGen/Lambda/Proton/hPt", "Proton p_{T}", kTH1F, {axisTrackPt}); histos.add("McGen/Lambda/Proton/hEta", "Proton #eta", kTH1F, {axisV0Eta}); histos.add("McGen/Lambda/Proton/hRap", "Proton y", kTH1F, {axisV0Rap}); histos.add("McGen/Lambda/Proton/hPhi", "Proton #phi", kTH1F, {axisV0Phi}); @@ -561,120 +519,45 @@ struct LambdaTableProducer { histos.addClone("McGen/Lambda/Proton/", "McGen/AntiLambda/Proton/"); histos.addClone("McGen/Lambda/Pion/", "McGen/AntiLambda/Pion/"); - // set bin lables specific to MC - histos.get(HIST("Events/h1f_collisions_info")) - ->GetXaxis() - ->SetBinLabel(CollisionLabels::kTotColBeforeHasMcCollision, - "kTotColBeforeHasMcCollision"); - histos.get(HIST("McGen/h1f_collisions_info")) - ->GetXaxis() - ->SetBinLabel(CollisionLabels::kTotCol, "kTotCol"); - histos.get(HIST("McGen/h1f_collisions_info")) - ->GetXaxis() - ->SetBinLabel(CollisionLabels::kPassSelCol, "kPassSelCol"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kTracksBeforeHasMcParticle, - "kTracksBeforeHasMcParticle"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPrimaryLambda, "kPrimaryLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kSecondaryLambda, "kSecondaryLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kLambdaDauNotMcParticle, - "kLambdaDauNotMcParticle"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kLambdaNotPrPiMinus, - "kLambdaNotPrPiMinus"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kAntiLambdaNotAntiPrPiPlus, - "kAntiLambdaNotAntiPrPiPlus"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPassTrueLambdaSel, "kPassTrueLambdaSel"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kGenTotAccLambda, "kGenTotAccLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kGenLambdaNoDau, "kGenLambdaNoDau"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kGenLambdaToPrPi, "kGenLambdaToPrPi"); - } - - // set bin labels - histos.get(HIST("Events/h1f_collisions_info")) - ->GetXaxis() - ->SetBinLabel(CollisionLabels::kTotCol, "kTotCol"); - histos.get(HIST("Events/h1f_collisions_info")) - ->GetXaxis() - ->SetBinLabel(CollisionLabels::kPassSelCol, "kPassSelCol"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kAllV0Tracks, "kAllV0Tracks"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kV0KShortMassRej, "kV0KShortMassRej"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kNotLambdaNotAntiLambda, - "kNotLambdaNotAntiLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kV0IsBothLambdaAntiLambda, - "kV0IsBothLambdaAntiLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kNotLambdaAfterSel, "kNotLambdaAfterSel"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kV0IsLambdaOrAntiLambda, - "kV0IsLambdaOrAntiLambda"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPassV0DauTrackSel, "kPassV0DauTrackSel"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPassV0KinCuts, "kPassV0KinCuts"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPassV0TopoSel, "kPassV0TopoSel"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kAllSelPassed, "kAllSelPassed"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kEffCorrPtCent, "kEffCorrPtCent"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kEffCorrPtRapCent, "kEffCorrPtRapCent"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kNoEffCorr, "kNoEffCorr"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPFCorrPtCent, "kPFCorrPtCent"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kPFCorrPtRapCent, "kPFCorrPtRapCent"); - histos.get(HIST("Tracks/h1f_tracks_info")) - ->GetXaxis() - ->SetBinLabel(TrackLabels::kNoPFCorr, "kNoPFCorr"); + histos.get(HIST("Events/h1f_collisions_info"))->GetXaxis()->SetBinLabel(CollisionLabels::kTotColBeforeHasMcCollision, "kTotColBeforeHasMcCollision"); + histos.get(HIST("McGen/h1f_collisions_info"))->GetXaxis()->SetBinLabel(CollisionLabels::kTotCol, "kTotCol"); + histos.get(HIST("McGen/h1f_collisions_info"))->GetXaxis()->SetBinLabel(CollisionLabels::kPassSelCol, "kPassSelCol"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kTracksBeforeHasMcParticle, "kTracksBeforeHasMcParticle"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPrimaryLambda, "kPrimaryLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kSecondaryLambda, "kSecondaryLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kLambdaDauNotMcParticle, "kLambdaDauNotMcParticle"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kLambdaNotPrPiMinus, "kLambdaNotPrPiMinus"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kAntiLambdaNotAntiPrPiPlus, "kAntiLambdaNotAntiPrPiPlus"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPassTrueLambdaSel, "kPassTrueLambdaSel"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kGenTotAccLambda, "kGenTotAccLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kGenLambdaNoDau, "kGenLambdaNoDau"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kGenLambdaToPrPi, "kGenLambdaToPrPi"); + } + + histos.get(HIST("Events/h1f_collisions_info"))->GetXaxis()->SetBinLabel(CollisionLabels::kTotCol, "kTotCol"); + histos.get(HIST("Events/h1f_collisions_info"))->GetXaxis()->SetBinLabel(CollisionLabels::kPassSelCol, "kPassSelCol"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kAllV0Tracks, "kAllV0Tracks"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kV0KShortMassRej, "kV0KShortMassRej"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kNotLambdaNotAntiLambda, "kNotLambdaNotAntiLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kV0IsBothLambdaAntiLambda, "kV0IsBothLambdaAntiLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kNotLambdaAfterSel, "kNotLambdaAfterSel"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kV0IsLambdaOrAntiLambda, "kV0IsLambdaOrAntiLambda"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPassV0DauTrackSel, "kPassV0DauTrackSel"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPassV0KinCuts, "kPassV0KinCuts"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPassV0TopoSel, "kPassV0TopoSel"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kAllSelPassed, "kAllSelPassed"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kEffCorrPtCent, "kEffCorrPtCent"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kEffCorrPtRapCent, "kEffCorrPtRapCent"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kNoEffCorr, "kNoEffCorr"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPFCorrPtCent, "kPFCorrPtCent"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kPFCorrPtRapCent, "kPFCorrPtRapCent"); + histos.get(HIST("Tracks/h1f_tracks_info"))->GetXaxis()->SetBinLabel(TrackLabels::kNoPFCorr, "kNoPFCorr"); } template bool selCollision(C const& col) { - if (col.posZ() <= cMinZVtx || col.posZ() >= cMaxZVtx) { - return false; - } - + histos.fill(HIST("Events/hEventCutFlow"), kEvAll); if constexpr (run == kRun3) { if (cCentEstimator == kCentFT0M) { cent = col.centFT0M(); @@ -693,50 +576,59 @@ struct LambdaTableProducer { return false; } } + histos.fill(HIST("Events/hEventCutFlow"), kEvTrigger); - if (cent <= cMinMult || cent >= cMaxMult) { - return false; - } if (cTriggerTvxSel && !col.selection_bit(aod::evsel::kIsTriggerTVX)) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvTvx); if (cTFBorder && !col.selection_bit(aod::evsel::kNoTimeFrameBorder)) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvTFBorder); if (cNoItsROBorder && !col.selection_bit(aod::evsel::kNoITSROFrameBorder)) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvItsRofBorder); if (cItsTpcVtx && !col.selection_bit(aod::evsel::kIsVertexITSTPC)) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvItsTpcVtx); if (cPileupReject && !col.selection_bit(aod::evsel::kNoSameBunchPileup)) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvNoSameBunchPileup); if (cZVtxTimeDiff && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; } - if (cIsGoodITSLayers && - !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { + histos.fill(HIST("Events/hEventCutFlow"), kEvGoodZvtxFT0vsPV); + if (cIsGoodITSLayers && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll)) { + return false; + } + histos.fill(HIST("Events/hEventCutFlow"), kEvGoodItsLayers); + + if (cent <= cMinMult || cent >= cMaxMult) { + return false; + } + histos.fill(HIST("Events/hEventCutFlow"), kEvCentrality); + if (col.posZ() <= cMinZVtx || col.posZ() >= cMaxZVtx) { return false; } + histos.fill(HIST("Events/hEventCutFlow"), kEvVz); mult = col.multNTracksPV(); return true; } - // Kinematic Selection - bool kinCutSelection(float const& pt, float const& rap, float const& ptMin, - float const& ptMax, float const& rapMax) + bool kinCutSelection(float const& pt, float const& rap, float const& ptMin, float const& ptMax, float const& rapMax) { return pt > ptMin && pt < ptMax && rap < rapMax; } - // Track Selection template bool selTrack(T const& track) { - if (!kinCutSelection(track.pt(), std::abs(track.eta()), cTrackMinPt, - cTrackMaxPt, cTrackEtaCut)) { + if (!kinCutSelection(track.pt(), std::abs(track.eta()), cTrackMinPt, cTrackMaxPt, cTrackEtaCut)) { return false; } if (track.tpcNClsCrossedRows() <= cMinTpcCrossedRows) { @@ -754,7 +646,6 @@ struct LambdaTableProducer { return true; } - // Daughter Track Selection template bool selDaughterTracks(V const& v0, T const&, ParticleType const& v0Type) { @@ -764,8 +655,6 @@ struct LambdaTableProducer { return false; } - // Apply DCA Selection on Daughter Tracks Based on Lambda/AntiLambda - // daughters float dcaProton = 0., dcaPion = 0.; if (v0Type == kLambda) { dcaProton = std::abs(v0.dcapostopv()); @@ -783,20 +672,17 @@ struct LambdaTableProducer { template bool topoCutSelection(C const& col, V const& v0, T const&) { - if (v0.dcaV0daughters() <= cMinV0DcaDaughters || - v0.dcaV0daughters() >= cMaxV0DcaDaughters) { + if (v0.dcaV0daughters() <= cMinV0DcaDaughters || v0.dcaV0daughters() >= cMaxV0DcaDaughters) { return false; } if (v0.dcav0topv() <= cMinDcaV0ToPV || v0.dcav0topv() >= cMaxDcaV0ToPV) { return false; } - if (v0.v0radius() <= cMinV0TransRadius || - v0.v0radius() >= cMaxV0TransRadius) { + if (v0.v0radius() <= cMinV0TransRadius || v0.v0radius() >= cMaxV0TransRadius) { return false; } - float ctau = - v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * MassLambda0; + float ctau = v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * MassLambda0; if (ctau <= cMinV0CTau || ctau >= cMaxV0CTau) { return false; } @@ -817,34 +703,27 @@ struct LambdaTableProducer { tpcNSigmaPr = negtrack.tpcNSigmaPr(); tpcNSigmaPi = postrack.tpcNSigmaPi(); } - return (std::abs(tpcNSigmaPr) < cTpcNsigmaCut && - std::abs(tpcNSigmaPi) < cTpcNsigmaCut); + return (std::abs(tpcNSigmaPr) < cTpcNsigmaCut && std::abs(tpcNSigmaPi) < cTpcNsigmaCut); } template bool selLambdaMassWindow(V const& v0, T const&, ParticleType& v0type) { - // Kshort mass rejection hypothesis - if (cKshortRejFlag && - (std::abs(v0.mK0Short() - MassK0Short) <= cKshortRejMassWindow)) { + if (cKshortRejFlag && (std::abs(v0.mK0Short() - MassK0Short) <= cKshortRejMassWindow)) { histos.fill(HIST("Tracks/h1f_tracks_info"), kV0KShortMassRej); return false; } - // initialize daughter tracks auto postrack = v0.template posTrack_as(); auto negtrack = v0.template negTrack_as(); - // initialize selection flags bool lambdaFlag = false, antiLambdaFlag = false; - if ((v0.mLambda() > cMinV0Mass && v0.mLambda() < cMaxV0Mass) && - selLambdaDauWithTpcPid(postrack, negtrack)) { + if ((v0.mLambda() > cMinV0Mass && v0.mLambda() < cMaxV0Mass) && selLambdaDauWithTpcPid(postrack, negtrack)) { lambdaFlag = true; v0type = kLambda; } - if ((v0.mAntiLambda() > cMinV0Mass && v0.mAntiLambda() < cMaxV0Mass) && - selLambdaDauWithTpcPid(postrack, negtrack)) { + if ((v0.mAntiLambda() > cMinV0Mass && v0.mAntiLambda() < cMaxV0Mass) && selLambdaDauWithTpcPid(postrack, negtrack)) { antiLambdaFlag = true; v0type = kAntiLambda; } @@ -861,8 +740,7 @@ struct LambdaTableProducer { } template - bool selV0Particle(C const& col, V const& v0, T const& tracks, - ParticleType& v0Type) + bool selV0Particle(C const& col, V const& v0, T const& tracks, ParticleType& v0Type) { if (!selLambdaMassWindow(v0, tracks, v0Type)) { return false; @@ -885,7 +763,6 @@ struct LambdaTableProducer { } histos.fill(HIST("Tracks/h1f_tracks_info"), kPassV0TopoSel); - // All Selection Criterion Passed return true; } @@ -960,7 +837,6 @@ struct LambdaTableProducer { return 1.; } - // Get from CCDB auto ccdbObj = ccdb->getForTimeStamp(cPathCCDB.value, 1); if (!ccdbObj) { LOGF(warning, "CCDB OBJECT NOT FOUND"); @@ -970,7 +846,6 @@ struct LambdaTableProducer { float effCorrFact = 1., primFrac = 1.; float rap = (cDoEtaAnalysis) ? v0.eta() : v0.yLambda(); - // Get Efficiency Factor if (cGetEffFact) { auto* objEff = ccdbObj->FindObject( Form("%s", vCorrFactStrings[cCorrFactHist][part].c_str())); @@ -1021,26 +896,20 @@ struct LambdaTableProducer { template void fillLambdaQAHistos(C const& col, V const& v0, T const&) { - static constexpr std::array SubDir = { - "QA/Lambda/", "QA/AntiLambda/"}; + static constexpr std::array SubDir = {"QA/Lambda/", "QA/AntiLambda/"}; auto postrack = v0.template posTrack_as(); auto negtrack = v0.template negTrack_as(); float mass = (part == kLambda) ? v0.mLambda() : v0.mAntiLambda(); float e = RecoDecay::e(v0.px(), v0.py(), v0.pz(), mass); float gamma = e / mass; - float ctau = - v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * MassLambda0; + float ctau = v0.distovertotmom(col.posX(), col.posY(), col.posZ()) * MassLambda0; float gctau = ctau * gamma; - histos.fill(HIST(SubDir[part]) + HIST("h2f_qt_vs_alpha"), v0.alpha(), - v0.qtarm()); - histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_V0_daughters"), - v0.dcaV0daughters()); - histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_pos_to_PV"), - v0.dcapostopv()); - histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_neg_to_PV"), - v0.dcanegtopv()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_qt_vs_alpha"), v0.alpha(), v0.qtarm()); + histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_V0_daughters"), v0.dcaV0daughters()); + histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_pos_to_PV"), v0.dcapostopv()); + histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_neg_to_PV"), v0.dcanegtopv()); histos.fill(HIST(SubDir[part]) + HIST("h1f_dca_V0_to_PV"), v0.dcav0topv()); histos.fill(HIST(SubDir[part]) + HIST("h1f_V0_cospa"), v0.v0cosPA()); histos.fill(HIST(SubDir[part]) + HIST("h1f_V0_radius"), v0.v0radius()); @@ -1052,45 +921,31 @@ struct LambdaTableProducer { histos.fill(HIST(SubDir[part]) + HIST("h1f_neg_prong_eta"), negtrack.eta()); histos.fill(HIST(SubDir[part]) + HIST("h1f_pos_prong_phi"), postrack.phi()); histos.fill(HIST(SubDir[part]) + HIST("h1f_neg_prong_phi"), negtrack.phi()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_dcaXY_vs_pt"), - postrack.pt(), postrack.dcaXY()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_dcaXY_vs_pt"), - negtrack.pt(), negtrack.dcaXY()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_dEdx_vs_p"), - postrack.tpcInnerParam(), postrack.tpcSignal()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_dEdx_vs_p"), - negtrack.tpcInnerParam(), negtrack.tpcSignal()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_tpc_nsigma_pr_vs_p"), - postrack.tpcInnerParam(), postrack.tpcNSigmaPr()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_tpc_nsigma_pr_vs_p"), - negtrack.tpcInnerParam(), negtrack.tpcNSigmaPr()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_tpc_nsigma_pi_vs_p"), - postrack.tpcInnerParam(), postrack.tpcNSigmaPi()); - histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_tpc_nsigma_pi_vs_p"), - negtrack.tpcInnerParam(), negtrack.tpcNSigmaPi()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_dcaXY_vs_pt"), postrack.pt(), postrack.dcaXY()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_dcaXY_vs_pt"), negtrack.pt(), negtrack.dcaXY()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_dEdx_vs_p"), postrack.tpcInnerParam(), postrack.tpcSignal()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_dEdx_vs_p"), negtrack.tpcInnerParam(), negtrack.tpcSignal()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_tpc_nsigma_pr_vs_p"), postrack.tpcInnerParam(), postrack.tpcNSigmaPr()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_tpc_nsigma_pr_vs_p"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPr()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_pos_prong_tpc_nsigma_pi_vs_p"), postrack.tpcInnerParam(), postrack.tpcNSigmaPi()); + histos.fill(HIST(SubDir[part]) + HIST("h2f_neg_prong_tpc_nsigma_pi_vs_p"), negtrack.tpcInnerParam(), negtrack.tpcNSigmaPi()); } template - void fillKinematicHists(float const& pt, float const& eta, float const& y, - float const& phi) + void fillKinematicHists(float const& pt, float const& eta, float const& y, float const& phi) { - static constexpr std::array SubDirRG = {"McRec/", - "McGen/"}; - static constexpr std::array SubDirPart = { - "Lambda/", "AntiLambda/"}; + static constexpr std::array SubDirRG = {"McRec/", "McGen/"}; + static constexpr std::array SubDirPart = {"Lambda/", "AntiLambda/"}; histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hPt"), pt); - histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hEta"), - eta); + histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hEta"), eta); histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hRap"), y); - histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hPhi"), - phi); + histos.fill(HIST(SubDirRG[rg]) + HIST(SubDirPart[part]) + HIST("hPhi"), phi); } template - void fillLambdaRecoTables(C const& collision, B const& bc, V const& v0tracks, - T const& tracks) + void fillLambdaRecoTables(C const& collision, B const& bc, V const& v0tracks, T const& tracks) { histos.fill(HIST("Events/h1f_collisions_info"), kTotCol); @@ -1102,17 +957,16 @@ struct LambdaTableProducer { histos.fill(HIST("Events/h1f_collisions_info"), kPassSelCol); histos.fill(HIST("Events/h1f_collision_posZ"), collision.posZ()); + histos.fill(HIST("Events/h1f_collision_cent"), cent); - // Fill Collision Table - lambdaCollisionTable(cent, mult, collision.posX(), collision.posY(), - collision.posZ(), bc.timestamp()); + lambdaCollisionTable(cent, mult, collision.posX(), collision.posY(), collision.posZ(), bc.timestamp()); - // initialize v0track objects ParticleType v0Type = kLambda; PrmScdType v0PrmScdType = kPrimary; float mass = 0., corr_fact = 1.; float prPx = 0., prPy = 0., prPz = 0.; float pt = 0., eta = 0., rap = 0., phi = 0.; + std::size_t nWritten = 0; for (auto const& v0 : v0tracks) { if constexpr (dmc == kMC) { @@ -1147,8 +1001,7 @@ struct LambdaTableProducer { phi = v0.phi(); if constexpr (dmc == kMC) { - histos.fill(HIST("Tracks/h2f_tracks_pid_before_sel"), - v0.mcParticle().pdgCode(), v0.pt()); + histos.fill(HIST("Tracks/h2f_tracks_pid_before_sel"), v0.mcParticle().pdgCode(), v0.pt()); if (cSelMCPSV0) { v0PrmScdType = isPrimaryV0(v0); } @@ -1159,8 +1012,7 @@ struct LambdaTableProducer { fillLambdaMothers(v0, tracks); } histos.fill(HIST("Tracks/h1f_tracks_info"), kPassTrueLambdaSel); - histos.fill(HIST("Tracks/h2f_tracks_pid_after_sel"), - v0.mcParticle().pdgCode(), v0.pt()); + histos.fill(HIST("Tracks/h2f_tracks_pid_after_sel"), v0.mcParticle().pdgCode(), v0.pt()); if (cRecoMomResoFlag) { auto mc = v0.template mcParticle_as(); pt = mc.pt(); @@ -1168,16 +1020,14 @@ struct LambdaTableProducer { rap = mc.y(); phi = mc.phi(); float y = cDoEtaAnalysis ? eta : rap; - if (!kinCutSelection(pt, std::abs(y), cMinV0Pt, cMaxV0Pt, - cMaxV0Rap)) { + if (!kinCutSelection(pt, std::abs(y), cMinV0Pt, cMaxV0Pt, cMaxV0Rap)) { continue; } } } histos.fill(HIST("Tracks/h2f_armpod_after_sel"), v0.alpha(), v0.qtarm()); - corr_fact = (v0Type == kLambda) ? getCorrectionFactors(v0) - : getCorrectionFactors(v0); + corr_fact = (v0Type == kLambda) ? getCorrectionFactors(v0) : getCorrectionFactors(v0); if (v0Type == kLambda) { prPx = v0.pxpos(); @@ -1185,33 +1035,39 @@ struct LambdaTableProducer { prPz = v0.pzpos(); histos.fill(HIST("Tracks/h1f_lambda_pt_vs_invm"), mass, v0.pt()); fillLambdaQAHistos(collision, v0, tracks); - fillKinematicHists(v0.pt(), v0.eta(), v0.yLambda(), - v0.phi()); + fillKinematicHists(v0.pt(), v0.eta(), v0.yLambda(), v0.phi()); } else { prPx = v0.pxneg(); prPy = v0.pyneg(); prPz = v0.pzneg(); histos.fill(HIST("Tracks/h1f_antilambda_pt_vs_invm"), mass, v0.pt()); fillLambdaQAHistos(collision, v0, tracks); - fillKinematicHists(v0.pt(), v0.eta(), v0.yLambda(), - v0.phi()); + fillKinematicHists(v0.pt(), v0.eta(), v0.yLambda(), v0.phi()); } lambdaTrackTable(lambdaCollisionTable.lastIndex(), v0.px(), v0.py(), v0.pz(), pt, eta, phi, rap, mass, prPx, prPy, prPz, v0.template posTrack_as().index(), - v0.template negTrack_as().index(), v0.v0cosPA(), - v0.dcaV0daughters(), (int8_t)v0Type, v0PrmScdType, + v0.template negTrack_as().index(), + v0.v0cosPA(), + v0.dcaV0daughters(), (int8_t)v0Type, + v0PrmScdType, corr_fact); + ++nWritten; + } + + if (nWritten >= 1) { + histos.fill(HIST("Events/hEventCutFlow"), kEvOneV0); + } + if (nWritten >= NCandidatesForPair) { + histos.fill(HIST("Events/hEventCutFlow"), kEvTwoV0); } } template - void fillLambdaMcGenTables(B const& bc, C const& mcCollision, - M const& mcParticles) + void fillLambdaMcGenTables(B const& bc, C const& mcCollision, M const& mcParticles) { - lambdaMCGenCollisionTable(cent, mult, mcCollision.posX(), - mcCollision.posY(), mcCollision.posZ(), + lambdaMCGenCollisionTable(cent, mult, mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), bc.timestamp()); ParticleType v0Type = kLambda; @@ -1293,15 +1149,11 @@ struct LambdaTableProducer { histos.fill(HIST("McGen/Lambda/Pion/hEta"), vDauEta[piIdx]); histos.fill(HIST("McGen/Lambda/Pion/hRap"), vDauRap[piIdx]); histos.fill(HIST("McGen/Lambda/Pion/hPhi"), vDauPhi[piIdx]); - fillKinematicHists(mcpart.pt(), mcpart.eta(), mcpart.y(), - mcpart.phi()); + fillKinematicHists(mcpart.pt(), mcpart.eta(), mcpart.y(), mcpart.phi()); } else { - histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), - daughterPDGs[prIdx]); - histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), - daughterPDGs[piIdx]); - histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), - mcpart.pdgCode()); + histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), daughterPDGs[prIdx]); + histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), daughterPDGs[piIdx]); + histos.fill(HIST("McGen/h1f_antilambda_daughter_PDG"), mcpart.pdgCode()); histos.fill(HIST("McGen/AntiLambda/Proton/hPt"), vDauPt[prIdx]); histos.fill(HIST("McGen/AntiLambda/Proton/hEta"), vDauEta[prIdx]); histos.fill(HIST("McGen/AntiLambda/Proton/hRap"), vDauRap[prIdx]); @@ -1310,8 +1162,7 @@ struct LambdaTableProducer { histos.fill(HIST("McGen/AntiLambda/Pion/hEta"), vDauEta[piIdx]); histos.fill(HIST("McGen/AntiLambda/Pion/hRap"), vDauRap[piIdx]); histos.fill(HIST("McGen/AntiLambda/Pion/hPhi"), vDauPhi[piIdx]); - fillKinematicHists(mcpart.pt(), mcpart.eta(), - mcpart.y(), mcpart.phi()); + fillKinematicHists(mcpart.pt(), mcpart.eta(), mcpart.y(), mcpart.phi()); } lambdaMCGenTrackTable( @@ -1323,10 +1174,8 @@ struct LambdaTableProducer { } } - template - void analyzeMcRecoGen(M const& mcCollision, C const& collisions, B const&, - V const& V0s, T const& tracks, P const& mcParticles) + template + void analyzeMcRecoGen(M const& mcCollision, C const& collisions, B const&, V const& V0s, T const& tracks, P const& mcParticles) { int nRecCols = collisions.size(); if (nRecCols != 0) { @@ -1356,32 +1205,25 @@ struct LambdaTableProducer { } SliceCache cache; - Preslice> perCollision = - aod::v0data::collisionId; + Preslice> perCollision = aod::v0data::collisionId; - using CollisionsRun3 = soa::Join; - using CollisionsRun2 = - soa::Join; + using CollisionsRun3 = soa::Join; + using CollisionsRun2 = soa::Join; using Tracks = soa::Join; - using TracksRun2 = - soa::Join; + using TracksRun2 = soa::Join; using TracksMC = soa::Join; using TracksMCRun2 = soa::Join; using McV0Tracks = soa::Join; - void processDataRun3(CollisionsRun3::iterator const& collision, - aod::BCsWithTimestamps const&, aod::V0Datas const& V0s, - Tracks const& tracks) + void processDataRun3(CollisionsRun3::iterator const& collision, aod::BCsWithTimestamps const&, aod::V0Datas const& V0s, Tracks const& tracks) { auto bc = collision.bc_as(); fillLambdaRecoTables(collision, bc, V0s, tracks); } - PROCESS_SWITCH(LambdaTableProducer, processDataRun3, "Process for Run3 DATA", - true); + PROCESS_SWITCH(LambdaTableProducer, processDataRun3, "Process for Run3 DATA", true); void processMCRun3( aod::McCollisions::iterator const& mcCollision, @@ -1389,38 +1231,29 @@ struct LambdaTableProducer { aod::BCsWithTimestamps const& bcts, McV0Tracks const& V0s, TracksMC const& tracks, aod::McParticles const& mcParticles) { - analyzeMcRecoGen(mcCollision, collisions, bcts, V0s, tracks, - mcParticles); + analyzeMcRecoGen(mcCollision, collisions, bcts, V0s, tracks, mcParticles); } - PROCESS_SWITCH(LambdaTableProducer, processMCRun3, - "Process for Run3 MC RecoGen", false); + PROCESS_SWITCH(LambdaTableProducer, processMCRun3, "Process for Run3 MC RecoGen", false); }; struct LambdaTracksExtProducer { Produces lambdaTrackExtTable; - Configurable cAcceptAllLambda{"cAcceptAllLambda", false, - "Accept all lambda (ignore sharing)"}; - Configurable cRejAllLambdaShaDau{"cRejAllLambdaShaDau", true, - "Reject lambda sharing daughters"}; - Configurable cSelLambdaMassPdg{"cSelLambdaMassPdg", false, - "Select lambda closest to PDG mass"}; - Configurable cSelLambdaTScore{"cSelLambdaTScore", false, - "Select lambda by t-score"}; + Configurable cAcceptAllLambda{"cAcceptAllLambda", false, "Accept all lambda (ignore sharing)"}; + Configurable cRejAllLambdaShaDau{"cRejAllLambdaShaDau", true, "Reject lambda sharing daughters"}; + Configurable cSelLambdaMassPdg{"cSelLambdaMassPdg", false, "Select lambda closest to PDG mass"}; + Configurable cSelLambdaTScore{"cSelLambdaTScore", false, "Select lambda by t-score"}; Configurable cA{"cA", 0.6, "t-score weight: |mass - PDGmass|"}; Configurable cB{"cB", 0.6, "t-score weight: DCA daughters"}; Configurable cC{"cC", 0.6, "t-score weight: |cosPA - 1|"}; - HistogramRegistry histos{ - "histos", - {}, - OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; void init(InitContext const&) { const AxisSpec axisMult(10, 0, 10); - const AxisSpec axisMass(120, 1.08, 1.20, "M_{p#pi} (GeV/#it{c}^{2})"); + const AxisSpec axisMass(200, 1.08, 1.18, "M_{p#pi} (GeV/#it{c}^{2})"); const AxisSpec axisCPA(100, 0.995, 1.0, "cos(#theta_{PA})"); const AxisSpec axisDcaDau(75, 0., 1.5, "Daughter DCA (#sigma)"); const AxisSpec axisDEta(320, -1.6, 1.6, "#Delta#eta"); @@ -1430,48 +1263,36 @@ 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}^{Share} #Lambda#bar{#Lambda}", - kTH2D, {axisDEta, axisDPhi}); - histos.add("h2d_n2_etaphi_LaM_LaP", "#rho_{2}^{Share} #bar{#Lambda}#Lambda", - kTH2D, {axisDEta, axisDPhi}); - histos.add("h2d_n2_etaphi_LaP_LaP", "#rho_{2}^{Share} #Lambda#Lambda", - kTH2D, {axisDEta, axisDPhi}); - histos.add("h2d_n2_etaphi_LaM_LaM", - "#rho_{2}^{Share} #bar{#Lambda}#bar{#Lambda}", kTH2D, - {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaP_LaM", "#rho_{2}^{Share} #Lambda#bar{#Lambda}", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaM_LaP", "#rho_{2}^{Share} #bar{#Lambda}#Lambda", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaP_LaP", "#rho_{2}^{Share} #Lambda#Lambda", kTH2D, {axisDEta, axisDPhi}); + histos.add("h2d_n2_etaphi_LaM_LaM", "#rho_{2}^{Share} #bar{#Lambda}#bar{#Lambda}", kTH2D, {axisDEta, axisDPhi}); histos.add("Reco/h1f_lambda_invmass", "M_{#Lambda}", kTH1F, {axisMass}); histos.add("Reco/h1f_lambda_cospa", "cos(PA)", kTH1F, {axisCPA}); histos.add("Reco/h1f_lambda_dcadau", "DCA daughters", kTH1F, {axisDcaDau}); - histos.add("Reco/h1f_antilambda_invmass", "M_{#bar{#Lambda}}", kTH1F, - {axisMass}); + histos.add("Reco/h1f_antilambda_invmass", "M_{#bar{#Lambda}}", kTH1F, {axisMass}); histos.add("Reco/h1f_antilambda_cospa", "cos(PA)", kTH1F, {axisCPA}); - histos.add("Reco/h1f_antilambda_dcadau", "DCA daughters", kTH1F, - {axisDcaDau}); + histos.add("Reco/h1f_antilambda_dcadau", "DCA daughters", kTH1F, {axisDcaDau}); histos.addClone("Reco/", "SharingDau/"); } template void fillHistos(T const& track) { - static constexpr std::array SubDir = {"Reco/", - "SharingDau/"}; + static constexpr std::array SubDir = {"Reco/", "SharingDau/"}; if (track.v0Type() == kLambda) { histos.fill(HIST(SubDir[sd]) + HIST("h1f_lambda_invmass"), track.mass()); histos.fill(HIST(SubDir[sd]) + HIST("h1f_lambda_dcadau"), track.dcaDau()); histos.fill(HIST(SubDir[sd]) + HIST("h1f_lambda_cospa"), track.cosPA()); } else { - histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_invmass"), - track.mass()); - histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_dcadau"), - track.dcaDau()); - histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_cospa"), - track.cosPA()); + histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_invmass"), track.mass()); + histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_dcadau"), track.dcaDau()); + histos.fill(HIST(SubDir[sd]) + HIST("h1f_antilambda_cospa"), track.cosPA()); } } - void process(aod::LambdaCollisions::iterator const&, - aod::LambdaTracks const& tracks) + void process(aod::LambdaCollisions::iterator const&, aod::LambdaTracks const& tracks) { int nTotLambda = 0, nTotAntiLambda = 0, nSelLambda = 0, nSelAntiLambda = 0; @@ -1487,50 +1308,32 @@ struct LambdaTracksExtProducer { ++nTotAntiLambda; } - tLambda = (cA * std::abs(lambda.mass() - MassLambda0)) + - (cB * lambda.dcaDau()) + (cC * std::abs(lambda.cosPA() - 1.)); + tLambda = (cA * std::abs(lambda.mass() - MassLambda0)) + (cB * lambda.dcaDau()) + (cC * std::abs(lambda.cosPA() - 1.)); for (auto const& track : tracks) { if (lambda.index() == track.index()) { continue; } - if (lambda.posTrackId() == track.posTrackId() || - lambda.negTrackId() == track.negTrackId()) { + if (lambda.posTrackId() == track.posTrackId() || lambda.negTrackId() == track.negTrackId()) { vSharedDauLambdaIndex.push_back(track.index()); lambdaSharingDauFlag = true; 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)); + 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)); + 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 (std::abs(lambda.mass() - MassLambda0) > - std::abs(track.mass() - MassLambda0)) { + if (std::abs(lambda.mass() - MassLambda0) > std::abs(track.mass() - MassLambda0)) { lambdaMinDeltaMassFlag = false; } - tTrack = (cA * std::abs(track.mass() - MassLambda0)) + - (cB * track.dcaDau()) + (cC * std::abs(track.cosPA() - 1.)); + tTrack = (cA * std::abs(track.mass() - MassLambda0)) + (cB * track.dcaDau()) + (cC * std::abs(track.cosPA() - 1.)); if (tLambda > tTrack) { lambdaMinTScoreFlag = false; } @@ -1543,9 +1346,7 @@ struct LambdaTracksExtProducer { fillHistos(lambda); } - if (cAcceptAllLambda || (cRejAllLambdaShaDau && !lambdaSharingDauFlag) || - (cSelLambdaMassPdg && lambdaMinDeltaMassFlag) || - (cSelLambdaTScore && lambdaMinTScoreFlag)) { + if (cAcceptAllLambda || (cRejAllLambdaShaDau && !lambdaSharingDauFlag) || (cSelLambdaMassPdg && lambdaMinDeltaMassFlag) || (cSelLambdaTScore && lambdaMinTScoreFlag)) { trueLambdaFlag = true; } @@ -1557,8 +1358,7 @@ struct LambdaTracksExtProducer { } } - lambdaTrackExtTable(lambdaSharingDauFlag, vSharedDauLambdaIndex, - trueLambdaFlag); + lambdaTrackExtTable(lambdaSharingDauFlag, vSharedDauLambdaIndex, trueLambdaFlag); } if (nTotLambda != 0) { @@ -1577,299 +1377,438 @@ struct LambdaTracksExtProducer { }; struct LambdaSpinPolarization { - // Table producer Produces lambdaMixEvtCol; Produces lambdaMixEvtTrk; Produces lambdaMixEvtMGCol; Produces lambdaMixEvtMGTrk; - Configurable cNPtBins{"cNPtBins", 30, "N pT bins"}; - Configurable cMinPt{"cMinPt", 0.5f, "pT min (GeV/c)"}; - Configurable cMaxPt{"cMaxPt", 4.5f, "pT max (GeV/c)"}; - Configurable cNRapBins{"cNRapBins", 10, "N rapidity bins"}; - Configurable cMinRap{"cMinRap", -0.5f, "Rapidity min"}; - Configurable cMaxRap{"cMaxRap", 0.5f, "Rapidity max"}; - Configurable cNPhiBins{"cNPhiBins", 36, "N phi bins"}; + Configurable cMassAccMin{"cMassAccMin", 1.0956f, "Accepted mass min (GeV/c2); candidates outside the window are not paired"}; + Configurable cMassAccMax{"cMassAccMax", 1.1356f, "Accepted mass max (GeV/c2); candidates outside the window are not paired"}; + Configurable cMassHistMin{"cMassHistMin", 1.0806f, "Mass axes min (GeV/c2), below cMassAccMin"}; + Configurable cMassHistMax{"cMassHistMax", 1.1806f, "Mass axes max (GeV/c2), above cMassAccMax"}; + Configurable cNMassBins{"cNMassBins", 200, "Mass axes N bins"}; + ConfigurableAxis axisMEReplacedMass{"axisMEReplacedMass", {VARIABLE_WIDTH, 1.0956, 1.1031, 1.1081, 1.1231, 1.1281, 1.1356}, "ME: mass of the replaced leg, binned by the analysis windows LSB | gap | Sig | gap | RSB (add edges, e.g. of the 3-sigma windows, to keep more window sets)"}; + ConfigurableAxis axisDeltaR{"axisDeltaR", {VARIABLE_WIDTH, 0.0, 0.4, 0.8, 1.2, 1.8, 2.4, 3.1, 3.5}, "DeltaR bins (last bin: outside the analysis windows)"}; Configurable cNBinsCosTS{"cNBinsCosTS", 10, "N costheta* bins"}; - Configurable cNBinsDeltaR{"cNBinsDeltaR", 20, "N DeltaR bins"}; - - Configurable cMassHistMin{"cMassHistMin", 1.08f, - "Mass histogram min (GeV/c2)"}; - Configurable cMassHistMax{"cMassHistMax", 1.20f, - "Mass histogram max (GeV/c2)"}; - Configurable cNMassBins{"cNMassBins", 120, "Mass histogram N bins"}; - - Configurable cSigMinLambda{"cSigMinLambda", 1.108f, - "Signal region min (GeV/c2)"}; - Configurable cSigMaxLambda{"cSigMaxLambda", 1.123f, - "Signal region max (GeV/c2)"}; - Configurable cSbLeftMin{"cSbLeftMin", 1.080f, - "Left sideband min (GeV/c2)"}; - Configurable cSbLeftMax{"cSbLeftMax", 1.100f, - "Left sideband max (GeV/c2)"}; - Configurable cSbRightMin{"cSbRightMin", 1.135f, - "Right sideband min (GeV/c2)"}; - Configurable cSbRightMax{"cSbRightMax", 1.155f, - "Right sideband max (GeV/c2)"}; - - Configurable cInvBoostFlag{"cInvBoostFlag", true, "Inverse boost flag"}; - Configurable cDoAtlasMethod{"cDoAtlasMethod", false, - "Fill pair-boost (ATLAS) histograms"}; - Configurable cDoStarMethod{"cDoStarMethod", true, - "Fill lab-boost (STAR) histograms"}; - Configurable mixingParameter{"mixingParameter", 5, "ME pool depth"}; - Configurable cMEMode{"cMEMode", 1, - "ME mode: 0=standard, 1=kinematicConstrained"}; - - ConfigurableAxis cMultBins{"cMultBins", {VARIABLE_WIDTH, 0.f, 10.f, 30.f, 50.f, 80.f, 100.f}, "Multiplicity bins"}; - ConfigurableAxis axisCentME{"axisCentME", {VARIABLE_WIDTH, 0, 10, 30, 50, 100}, "ME centrality bins"}; - ConfigurableAxis axisVtxZME{"axisVtxZME", {VARIABLE_WIDTH, -7, -3, 0, 3, 7}, "ME vtxZ bins"}; + + Configurable cMEPoolDepth{"cMEPoolDepth", 100000, "Rolling ME pool: depth in donor events per (centrality, Vz) bin"}; + Configurable cMEPoolMinEvents{"cMEPoolMinEvents", 20, "Rolling ME pool: events a pool must hold before its bin is mixed (warm-up)"}; + Configurable cMEPoolMinCand{"cMEPoolMinCand", 2, "Rolling ME pool: only events with at least this many accepted candidates donate"}; + ConfigurableAxis axisCentME{"axisCentME", {VARIABLE_WIDTH, 0, 10, 20, 40, 100}, "ME pool centrality bins"}; + ConfigurableAxis axisVtxZME{"axisVtxZME", {VARIABLE_WIDTH, -10., -7., -4., -2., 0, 2., 4., 7., 10.}, "ME pool vtxZ bins"}; Configurable cMaxDeltaPt{"cMaxDeltaPt", 0.1f, "Kinematic ME: max |deltaPt(SE)-deltaPt(ME)| (GeV/c)"}; - Configurable cMaxDeltaPhi{"cMaxDeltaPhi", 0.1f, "Kinematic ME: max |deltaPhi(SE)-deltaPhi(ME)| (rad)"}; - Configurable cMaxDeltaRap{"cMaxDeltaRap", 0.1f, "Kinematic ME: max |deltaRap(SE)-deltaRap(ME)|"}; + Configurable cMaxDeltaPhi{"cMaxDeltaPhi", 0.02f, "Kinematic ME: max |deltaPhi(SE)-deltaPhi(ME)| (rad)"}; + Configurable cMaxDeltaRap{"cMaxDeltaRap", 0.03f, "Kinematic ME: max |deltaRap(SE)-deltaRap(ME)|"}; + Configurable cFillSEInMixing{"cFillSEInMixing", true, "Kinematic ME: fill every same-event pair the mixing processes, matched or not (SEproc)"}; + + Configurable cFillKinWSpectra{"cFillKinWSpectra", true, "Kinematic ME: fill the leg-2 spectra of SEproc and ME (weight-map input, closure)"}; + Configurable cMinPt{"cMinPt", 0.6f, "Leg spectra and ME pool grid: pT min (GeV/c)"}; + Configurable cMaxPt{"cMaxPt", 3.0f, "Leg spectra and ME pool grid: pT max (GeV/c)"}; + Configurable cMinRap{"cMinRap", -0.5f, "Leg spectra and ME pool grid: rapidity min"}; + Configurable cMaxRap{"cMaxRap", 0.5f, "Leg spectra and ME pool grid: rapidity max"}; + Configurable cNKinWPtBins{"cNKinWPtBins", 0, "Leg spectra: pT bins (0: no wider than cMaxDeltaPt)"}; + Configurable cNKinWRapBins{"cNKinWRapBins", 0, "Leg spectra: rapidity bins (0: no wider than cMaxDeltaRap)"}; + Configurable cNKinWPhiBins{"cNKinWPhiBins", 0, "Leg spectra: phi bins (0: no wider than cMaxDeltaPhi)"}; + Configurable cApplyKinWeights{"cApplyKinWeights", false, "Kinematic ME: weight the leg-2 replacements with the maps from CCDB"}; + Configurable cKinWeightCcdbUrl{"cKinWeightCcdbUrl", "http://alice-ccdb.cern.ch", "Kinematic weights: CCDB url"}; + Configurable cKinWeightCcdbPath{"cKinWeightCcdbPath", "", "Kinematic weights: CCDB path of a TList with hKinW_ (axes of KinW/*/hLeg2_)"}; + Configurable cKinWeightCcdbTimestamp{"cKinWeightCcdbTimestamp", -1, "Kinematic weights: CCDB timestamp (-1: timestamp of the first mixed collision)"}; + + Configurable cFillClosePairQA{"cFillClosePairQA", true, "Close-pair QA of the same-event pairs: SE (processDataReco) and SEproc (processDataRecoMixed)"}; + Configurable cClosePairMassWindow{"cClosePairMassWindow", 0.004f, "Close-pair QA: both legs within |m - m_Lambda| < this (GeV/c2); <= 0: all pairs"}; + ConfigurableAxis axisClosePair{"axisClosePair", {100, -0.15, 0.15}, "Close-pair QA: bins of the daughter Delta eta, Delta y and Delta phi"}; + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + Service ccdb{}; - HistogramRegistry histos{ - "histos", - {}, - OutputObjHandlingPolicy::AnalysisObject}; + uint64_t pairOrderSeed = 0; + std::array, 4> kinWeightMaps{}; + std::vector legAxesRef{}; + TAxis replacedMassAxis{}; + bool kinWeightsLoaded = false; - float cent = 0.; struct PoolTrack { - float _px, _py, _pz, _pt, _rap, _phi, _mass; - float _prPx, _prPy, _prPz; - [[nodiscard]] float px() const { return _px; } - [[nodiscard]] float py() const { return _py; } - [[nodiscard]] float pz() const { return _pz; } - [[nodiscard]] float pt() const { return _pt; } - [[nodiscard]] float rap() const { return _rap; } - [[nodiscard]] float phi() const { return _phi; } - [[nodiscard]] float mass() const { return _mass; } - [[nodiscard]] float prPx() const { return _prPx; } - [[nodiscard]] float prPy() const { return _prPy; } - [[nodiscard]] float prPz() const { return _prPz; } + float mPx, mPy, mPz, mPt, mRap, mPhi, mMass; + float mPrPx, mPrPy, mPrPz; + int8_t mSpecies; + [[nodiscard]] float px() const { return mPx; } + [[nodiscard]] float py() const { return mPy; } + [[nodiscard]] float pz() const { return mPz; } + [[nodiscard]] float pt() const { return mPt; } + [[nodiscard]] float rap() const { return mRap; } + [[nodiscard]] float phi() const { return mPhi; } + [[nodiscard]] float mass() const { return mMass; } + [[nodiscard]] float prPx() const { return mPrPx; } + [[nodiscard]] float prPy() const { return mPrPy; } + [[nodiscard]] float prPz() const { return mPrPz; } + [[nodiscard]] int8_t species() const { return mSpecies; } }; template PoolTrack toPoolTrack(T const& trk) { - return PoolTrack{trk.px(), trk.py(), trk.pz(), trk.pt(), trk.rap(), - trk.phi(), trk.mass(), trk.prPx(), trk.prPy(), trk.prPz()}; + return PoolTrack{trk.px(), trk.py(), trk.pz(), trk.pt(), trk.rap(), trk.phi(), trk.mass(), trk.prPx(), trk.prPy(), trk.prPz(), static_cast(trk.v0Type())}; } + class RollingPool + { + public: + struct Grid { + float ptMin = 0.f; + float rapMin = 0.f; + float cellPt = 1.f; + float cellRap = 1.f; + float cellPhi = TwoPI; + int nPt = 1; + int nRap = 1; + int nPhi = 1; + }; + + RollingPool() = default; + RollingPool(Grid const& g, int depth) : grid(g), maxEvents(depth) {} + + [[nodiscard]] std::size_t nEvents() const { return nPerEvent.size(); } + + void push(std::vector const& event) + { + for (auto const& cand : event) { + const int64_t seq = seqBase + static_cast(slots.size()); + slots.push_back(Slot{cand, NoSeq}); + auto [it, isNew] = cells[cand.species()].try_emplace(cellKey(cand), Cell{seq, seq}); + if (!isNew) { + slot(it->second.tail).next = seq; + it->second.tail = seq; + } + } + nPerEvent.push_back(static_cast(event.size())); + while (static_cast(nPerEvent.size()) > maxEvents) { + const int n = nPerEvent.front(); + for (int i = 0; i < n; ++i) { + retireOldest(); + } + nPerEvent.pop_front(); + } + } + + template + void findMatches(PoolTrack const& leg, F const& isMatch, std::vector& out) const + { + out.clear(); + int iPt = 0, iRap = 0, iPhi = 0; + cellCoords(leg, iPt, iRap, iPhi); + std::array phiCells{}; + int nPhiCells = 0; + for (int d = -1; d <= 1; ++d) { + const int k = ((iPhi + d) % grid.nPhi + grid.nPhi) % grid.nPhi; + if (std::find(phiCells.begin(), phiCells.begin() + nPhiCells, k) == phiCells.begin() + nPhiCells) { + phiCells[nPhiCells++] = k; + } + } + auto const& speciesCells = cells[leg.species()]; + for (int jPt = std::max(iPt - 1, 0); jPt <= std::min(iPt + 1, grid.nPt - 1); ++jPt) { + for (int jRap = std::max(iRap - 1, 0); jRap <= std::min(iRap + 1, grid.nRap - 1); ++jRap) { + for (int ip = 0; ip < nPhiCells; ++ip) { + const auto it = speciesCells.find(key(jPt, jRap, phiCells[ip])); + if (it == speciesCells.end()) { + continue; + } + for (int64_t seq = it->second.head; seq != NoSeq;) { + Slot const& sl = slot(seq); + if (isMatch(sl.cand)) { + out.push_back(&sl.cand); + } + seq = sl.next; + } + } + } + } + } + + private: + static constexpr int64_t NoSeq = -1; + struct Slot { + PoolTrack cand; + int64_t next; + }; + struct Cell { + int64_t head; + int64_t tail; + }; + + Slot& slot(int64_t seq) { return slots[static_cast(seq - seqBase)]; } + [[nodiscard]] Slot const& slot(int64_t seq) const { return slots[static_cast(seq - seqBase)]; } + + void retireOldest() + { + Slot const& old = slots.front(); + auto& speciesCells = cells[old.cand.species()]; + auto it = speciesCells.find(cellKey(old.cand)); + if (it == speciesCells.end() || it->second.head != seqBase) { + LOGF(fatal, "RollingPool: grid index out of step with the candidate queue"); + return; + } + if (old.next == NoSeq) { + speciesCells.erase(it); + } else { + it->second.head = old.next; + } + slots.pop_front(); + ++seqBase; + } + + [[nodiscard]] int64_t key(int iPt, int iRap, int iPhi) const + { + return (static_cast(iPt) * grid.nRap + iRap) * grid.nPhi + iPhi; + } + + void cellCoords(PoolTrack const& c, int& iPt, int& iRap, int& iPhi) const + { + iPt = std::clamp(static_cast(std::floor((c.pt() - grid.ptMin) / grid.cellPt)), 0, grid.nPt - 1); + iRap = std::clamp(static_cast(std::floor((c.rap() - grid.rapMin) / grid.cellRap)), 0, grid.nRap - 1); + iPhi = std::clamp(static_cast(std::floor(RecoDecay::constrainAngle(c.phi(), 0.f) / grid.cellPhi)), 0, grid.nPhi - 1); + } + + [[nodiscard]] int64_t cellKey(PoolTrack const& c) const + { + int iPt = 0, iRap = 0, iPhi = 0; + cellCoords(c, iPt, iRap, iPhi); + return key(iPt, iRap, iPhi); + } + + Grid grid{}; + int maxEvents = 0; + std::deque slots{}; + std::deque nPerEvent{}; + int64_t seqBase = 0; + std::array, 2> cells{}; + }; + + std::vector mePools{}; + std::vector mePoolsGen{}; + TAxis mePoolAxisCent{}; + TAxis mePoolAxisVz{}; + void init(InitContext const&) { - const AxisSpec axisCheck(1, 0, 1, ""); - const AxisSpec axisCent(cMultBins, "FT0M (%)"); - - const AxisSpec axisMass(cNMassBins, cMassHistMin, cMassHistMax, - "M_{#Lambda} (GeV/#it{c}^{2})"); - const AxisSpec axisPt(cNPtBins, cMinPt, cMaxPt, "p_{T} (GeV/#it{c})"); - const AxisSpec axisDRap(2 * cNRapBins, cMinRap - cMaxRap, cMaxRap - cMinRap, - "#Deltay"); - const AxisSpec axisDPhi(cNPhiBins, -PI, PI, "#Delta#varphi"); + const AxisSpec axisM1(cNMassBins, cMassHistMin, cMassHistMax, "M_{1} (GeV/#it{c}^{2})"); + const AxisSpec axisM2(cNMassBins, cMassHistMin, cMassHistMax, "M_{2} (GeV/#it{c}^{2})"); + const AxisSpec axisDR(axisDeltaR, "#DeltaR"); const AxisSpec axisCosTS(cNBinsCosTS, -1, 1, "cos(#theta*)"); - const AxisSpec axisDR(cNBinsDeltaR, 0, 3.5, "#DeltaR"); - const AxisSpec axisPosZ(200, -10, 10, "V_{z} (cm)"); - const AxisSpec axisMult(10, 0, 10, "N_{#Lambda}"); + const AxisSpec axisReplacedMass(axisMEReplacedMass, "M_{replaced leg} (GeV/#it{c}^{2})"); + const std::vector pairAxes = {axisM1, axisM2, axisDR, axisCosTS}; + const std::vector meAxes = {axisM1, axisM2, axisDR, axisCosTS, axisReplacedMass}; + + auto nBins = [](int configured, float range, float window) { + if (configured > 0) { + return configured; + } + return window > 0.f ? static_cast(std::ceil(range / window)) : 1; + }; + const AxisSpec axisLegPt(nBins(cNKinWPtBins, cMaxPt - cMinPt, cMaxDeltaPt), cMinPt, cMaxPt, "p_{T} (GeV/#it{c})"); + const AxisSpec axisLegRap(nBins(cNKinWRapBins, cMaxRap - cMinRap, cMaxDeltaRap), cMinRap, cMaxRap, "y"); + const AxisSpec axisLegPhi(nBins(cNKinWPhiBins, TwoPI, cMaxDeltaPhi), 0., TwoPI, "#varphi (rad)"); + const std::vector legAxes = {axisLegPt, axisLegRap, axisLegPhi, axisReplacedMass}; + legAxesRef.clear(); + for (auto const& spec : legAxes) { + const std::vector edges = axisEdges(spec); + legAxesRef.emplace_back(static_cast(edges.size()) - 1, edges.data()); + } + replacedMassAxis = legAxesRef.back(); + + if (!(cMassHistMin.value < cMassAccMin.value && cMassAccMax.value < cMassHistMax.value)) { + LOGF(fatal, "The mass axes [%.4f, %.4f] must be wider than the accepted mass window [%.4f, %.4f]", cMassHistMin.value, cMassHistMax.value, cMassAccMin.value, cMassAccMax.value); + } + std::vector massEdges = axisEdges(axisReplacedMass); + massEdges.push_back(cMassAccMin.value); + massEdges.push_back(cMassAccMax.value); + for (auto const& edge : massEdges) { + if (!isOnMassBinEdge(edge)) { + LOGF(warning, "Mass %.5f (accepted window or axisMEReplacedMass) is not a bin edge of the mass axes [%.4f, %.4f] / %d", edge, cMassHistMin.value, cMassHistMax.value, cNMassBins.value); + } + } + + if (doprocessDataRecoMixed || doprocessMcGenMixed) { + const std::vector centEdges = axisEdges(AxisSpec(axisCentME)); + const std::vector vzEdges = axisEdges(AxisSpec(axisVtxZME)); + mePoolAxisCent.Set(static_cast(centEdges.size()) - 1, centEdges.data()); + mePoolAxisVz.Set(static_cast(vzEdges.size()) - 1, vzEdges.data()); + + static constexpr std::array PoolSizeSteps = {1., 2., 5.}; + static constexpr double PoolSizeDecade = 10.; + const double depth = std::max(1, cMEPoolDepth.value); + std::vector poolSizeEdges = {0.}; + for (double decade = 1.; decade < depth; decade *= PoolSizeDecade) { + for (auto const& step : PoolSizeSteps) { + if (step * decade < depth) { + poolSizeEdges.push_back(step * decade); + } + } + } + poolSizeEdges.push_back(depth); + const AxisSpec axisPoolSize(poolSizeEdges, "events in the pool"); + + std::vector mixQaDirs; + if (doprocessDataRecoMixed) { + mixQaDirs.emplace_back("QA/ME/"); + } + if (doprocessMcGenMixed) { + mixQaDirs.emplace_back("McGen/QA/ME/"); + } + for (auto const& dir : mixQaDirs) { + histos.add((dir + "hMixedCentVz").c_str(), "events mixed;cent (%);V_{z} (cm)", kTH2F, {axisCentME, axisVtxZME}); + auto hStatus = histos.add((dir + "hEventStatus").c_str(), "rolling-pool mixing;;collisions", kTH1D, {{kNMEEventStatus - 1, 0.5, static_cast(kNMEEventStatus) - 0.5}}); + hStatus->GetXaxis()->SetBinLabel(kMEEventSeen, "seen"); + hStatus->GetXaxis()->SetBinLabel(kMEEventOutsidePools, "outside the pool bins"); + hStatus->GetXaxis()->SetBinLabel(kMEEventWarmUp, "pairs, pool in warm-up"); + hStatus->GetXaxis()->SetBinLabel(kMEEventMixed, "mixed"); + hStatus->GetXaxis()->SetBinLabel(kMEEventDonated, "donated to the pool"); + histos.add((dir + "hPoolEventsAtMixing").c_str(), "pool size when a collision is mixed;events in the pool;collisions", kTH1F, {axisPoolSize}); + histos.add((dir + "hMatchesLeg2").c_str(), "matches per replaced leg 2;N matches;pair type", kTH2F, {{100, -0.5, 99.5}, {4, -0.5, 3.5}}); + histos.add((dir + "hDeltaRShift").c_str(), "ME entries are filed at the #DeltaR of their seed;#DeltaR_{seed};#DeltaR_{mixed pair} - #DeltaR_{seed}", kTH2F, {{35, 0., 3.5}, {100, -0.2, 0.2}}); + histos.add((dir + "hLambdaMultVsCent").c_str(), "ME #Lambda mult;cent;N", kTH2D, {axisCentME, {50, 0, 50}}); + histos.add((dir + "hAntiLambdaMultVsCent").c_str(), "ME #bar{#Lambda} mult;cent;N", kTH2D, {axisCentME, {50, 0, 50}}); + } + + if (cMaxDeltaPt.value <= 0.f || cMaxDeltaRap.value <= 0.f || cMaxDeltaPhi.value <= 0.f) { + LOGF(fatal, "Rolling ME pool: the matching windows cMaxDeltaPt, cMaxDeltaRap and cMaxDeltaPhi must be positive"); + } + RollingPool::Grid grid; + grid.ptMin = cMinPt.value; + grid.rapMin = cMinRap.value; + grid.cellPt = cMaxDeltaPt.value; + grid.cellRap = cMaxDeltaRap.value; + grid.nPt = std::max(1, static_cast(std::ceil((cMaxPt.value - cMinPt.value) / grid.cellPt))); + grid.nRap = std::max(1, static_cast(std::ceil((cMaxRap.value - cMinRap.value) / grid.cellRap))); + grid.nPhi = std::max(1, static_cast(std::floor(TwoPI / cMaxDeltaPhi.value))); + grid.cellPhi = TwoPI / static_cast(grid.nPhi); + const auto nPools = static_cast(mePoolAxisCent.GetNbins() * mePoolAxisVz.GetNbins()); + if (doprocessDataRecoMixed) { + mePools.assign(nPools, RollingPool(grid, cMEPoolDepth.value)); + } + if (doprocessMcGenMixed) { + mePoolsGen.assign(nPools, RollingPool(grid, cMEPoolDepth.value)); + } + LOGF(info, "Rolling ME pools: %d (centrality) x %d (Vz) bins, depth %d donor events, grid %d x %d x %d (pT, y, phi)", + mePoolAxisCent.GetNbins(), mePoolAxisVz.GetNbins(), cMEPoolDepth.value, grid.nPt, grid.nRap, grid.nPhi); + } + + std::vector evDirs; + if (doprocessDataReco || doprocessDataRecoMixed) { + evDirs.emplace_back(""); + } + if (doprocessMcGen || doprocessMcGenMixed) { + evDirs.emplace_back("McGen/"); + } + const AxisSpec axisCandMass(cNMassBins, cMassHistMin, cMassHistMax, "M_{p#pi} (GeV/#it{c}^{2})"); + for (auto const& dir : evDirs) { + auto hCount = histos.add((dir + "Events/hEventCount").c_str(), "collisions of the pair analysis;;collisions", kTH1D, {{kNPairEventCount - 1, 0.5, static_cast(kNPairEventCount) - 0.5}}); + hCount->GetXaxis()->SetBinLabel(kPeAll, "all"); + hCount->GetXaxis()->SetBinLabel(kPeOneCand, "#geq 1 accepted #Lambda/#bar{#Lambda}"); + hCount->GetXaxis()->SetBinLabel(kPeTwoCand, "#geq 2 accepted (enter pairs)"); + hCount->GetXaxis()->SetBinLabel(kPeUnlikeSign, "#geq 1 #Lambda#bar{#Lambda} pair"); + hCount->GetXaxis()->SetBinLabel(kPeLambdaLambda, "#geq 1 #Lambda#Lambda pair"); + hCount->GetXaxis()->SetBinLabel(kPeAntiLambdaAntiLambda, "#geq 1 #bar{#Lambda}#bar{#Lambda} pair"); + histos.add((dir + "Events/hNCandAccepted").c_str(), "accepted candidates per collision;N_{#Lambda};N_{#bar{#Lambda}}", kTH2D, {{21, -0.5, 20.5}, {21, -0.5, 20.5}}); + histos.add((dir + "Events/hCentPairEvents").c_str(), "collisions with a pair;cent (%);collisions", kTH1D, {{100, 0., 100.}}); + for (auto const& species : {"Lambda/", "AntiLambda/"}) { + const std::string cand = dir + "QA/PairCand/" + species; + histos.add((cand + "hMassVsPt").c_str(), "accepted candidates of collisions with a pair", kTH2F, {axisCandMass, axisLegPt}); + histos.add((cand + "hRapVsPhi").c_str(), "accepted candidates of collisions with a pair", kTH2F, {axisLegRap, axisLegPhi}); + } + } + + static constexpr std::array PairTags = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; + for (auto const& tag : PairTags) { + const std::string pair{tag}; + if (doprocessDataReco) { + histos.add(("SE/hPair_" + pair).c_str(), ("SE " + pair).c_str(), kTHnSparseF, pairAxes); + } + if (doprocessDataRecoMixed) { + histos.add(("ME/hPair_" + pair).c_str(), ("ME " + pair).c_str(), kTHnSparseD, meAxes, true); + } + if (doprocessDataRecoMixed && cFillSEInMixing) { + histos.add(("SEproc/hPair_" + pair).c_str(), ("SE processed by the mixing " + pair).c_str(), kTHnSparseF, pairAxes); + } + if (doprocessDataRecoMixed && cFillKinWSpectra) { + histos.add(("KinW/SEproc/hLeg2_" + pair).c_str(), ("leg 2 of SEproc " + pair).c_str(), kTHnSparseD, legAxes); + histos.add(("KinW/ME/hLeg2_" + pair).c_str(), ("candidates replacing leg 2 " + pair).c_str(), kTHnSparseD, legAxes, true); + } + if (doprocessMcGen) { + histos.add(("McGen/SE/hPair_" + pair).c_str(), ("MC generated SE " + pair).c_str(), kTHnSparseF, pairAxes); + } + if (doprocessMcGenMixed) { + histos.add(("McGen/ME/hPair_" + pair).c_str(), ("MC generated ME " + pair).c_str(), kTHnSparseD, meAxes, true); + } + if (doprocessMcGenMixed && cFillSEInMixing) { + histos.add(("McGen/SEproc/hPair_" + pair).c_str(), ("MC generated SE processed by the mixing " + pair).c_str(), kTHnSparseF, pairAxes); + } + } + + if (cFillClosePairQA && (doprocessDataReco || doprocessDataRecoMixed)) { + const AxisSpec axisCPdEta(axisClosePair, "#Delta#eta"); + const AxisSpec axisCPdY(axisClosePair, "#Deltay"); + const AxisSpec axisCPdPhi(axisClosePair, "#Delta#varphi (rad)"); + std::vector> cpDirs; + if (doprocessDataReco) { + cpDirs.emplace_back("QA/ClosePair/", ""); + } + if (doprocessDataRecoMixed) { + cpDirs.emplace_back("QA/ClosePairSEproc/", "SEproc, "); + } + for (auto const& [dir, tag] : cpDirs) { + for (auto const& cls : {"LamLam", "ALamALam", "UnlikeSign"}) { + for (auto const& combo : {"pp", "pipi", "p1pi2", "pi1p2"}) { + const std::string name = dir + cls + "/h_" + combo; + const std::string title = tag + cls + ", daughters " + combo; + histos.add((name + "_dEta").c_str(), title.c_str(), kTH3F, {axisDR, axisCPdEta, axisCPdPhi}); + histos.add((name + "_dY").c_str(), title.c_str(), kTH3F, {axisDR, axisCPdY, axisCPdPhi}); + } + } + } + } - histos.add("QA/ME/hPoolCentVz", "ME pool;cent (%);V_{z}", kTH2F, - {axisCentME, axisVtxZME}); - histos.add("QA/ME/hLambdaMultVsCent", "ME #Lambda mult;cent;N", kTH2F, - {axisCentME, {50, 0, 50}}); - histos.add("QA/ME/hAntiLambdaMultVsCent", "ME #bar{#Lambda} mult;cent;N", - kTH2F, {axisCentME, {50, 0, 50}}); - - const std::vector massSparseAxes = { - axisCent, axisMass, axisMass, axisPt, axisPt, - axisDRap, axisDPhi, axisDR, axisCosTS}; - - histos.add("SE/Reco/Star/h2f_n2_mass_LaPLaM", - "M_{inv}: #Lambda#bar{#Lambda} inclusive", kTHnSparseF, - massSparseAxes); - histos.add("SE/Reco/Star/h2f_n2_mass_LaMLaP", - "M_{inv}: #bar{#Lambda}#Lambda inclusive", kTHnSparseF, - massSparseAxes); - histos.add("SE/Reco/Star/h2f_n2_mass_LaPLaP", - "M_{inv}: #Lambda#Lambda inclusive", kTHnSparseF, - massSparseAxes); - histos.add("SE/Reco/Star/h2f_n2_mass_LaMLaM", - "M_{inv}: #bar{#Lambda}#bar{#Lambda} inclusive", kTHnSparseF, - massSparseAxes); - - histos.add("SE/RecoBkgSigSB/Star/h2f_n2_mass_LaPLaM", - "M_{inv}: #Lambda#bar{#Lambda} Sig#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSigSB/Star/h2f_n2_mass_LaMLaP", - "M_{inv}: #bar{#Lambda}#Lambda Sig#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSigSB/Star/h2f_n2_mass_LaPLaP", - "M_{inv}: #Lambda#Lambda Sig#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSigSB/Star/h2f_n2_mass_LaMLaM", - "M_{inv}: #bar{#Lambda}#bar{#Lambda} Sig#timesSB", kTHnSparseF, - massSparseAxes); - - histos.add("SE/RecoBkgSBSB/Star/h2f_n2_mass_LaPLaM", - "M_{inv}: #Lambda#bar{#Lambda} SB#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSBSB/Star/h2f_n2_mass_LaMLaP", - "M_{inv}: #bar{#Lambda}#Lambda SB#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSBSB/Star/h2f_n2_mass_LaPLaP", - "M_{inv}: #Lambda#Lambda SB#timesSB", kTHnSparseF, - massSparseAxes); - histos.add("SE/RecoBkgSBSB/Star/h2f_n2_mass_LaMLaM", - "M_{inv}: #bar{#Lambda}#bar{#Lambda} SB#timesSB", kTHnSparseF, - massSparseAxes); - - histos.addClone("SE/Reco/Star/", "SE/Reco/Atlas/"); - histos.addClone("SE/RecoBkgSigSB/Star/", "SE/RecoBkgSigSB/Atlas/"); - histos.addClone("SE/RecoBkgSBSB/Star/", "SE/RecoBkgSBSB/Atlas/"); - - histos.addClone("SE/Reco/", "ME/Reco/"); - histos.addClone("SE/RecoBkgSigSB/", "ME/RecoBkgSigSB/"); - histos.addClone("SE/RecoBkgSBSB/", "ME/RecoBkgSBSB/"); - - histos.add("SE/RecoCorr/Star/h2f_n2_dltaR_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_dltaR_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_dltaR_LaPLaP", - "#rho_{2} #Lambda#Lambda [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_dltaR_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - - histos.add("SE/RecoCorr/Star/h2f_n2_ctheta_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_ctheta_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_ctheta_LaPLaP", - "#rho_{2} #Lambda#Lambda [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Star/h2f_n2_ctheta_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - - histos.add("SE/RecoCorr/Atlas/h2f_n2_dltaR_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_dltaR_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_dltaR_LaPLaP", - "#rho_{2} #Lambda#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_dltaR_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - - histos.add("SE/RecoCorr/Atlas/h2f_n2_ctheta_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_ctheta_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_ctheta_LaPLaP", - "#rho_{2} #Lambda#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("SE/RecoCorr/Atlas/h2f_n2_ctheta_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - - histos.add("ME/RecoCorr/Star/h2f_n2_dltaR_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_dltaR_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_dltaR_LaPLaP", - "#rho_{2} #Lambda#Lambda [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_dltaR_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - - histos.add("ME/RecoCorr/Star/h2f_n2_ctheta_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_ctheta_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_ctheta_LaPLaP", - "#rho_{2} #Lambda#Lambda [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Star/h2f_n2_ctheta_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Star]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - - histos.add("ME/RecoCorr/Atlas/h2f_n2_dltaR_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_dltaR_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_dltaR_LaPLaP", - "#rho_{2} #Lambda#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_dltaR_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDR, axisCosTS}); - - histos.add("ME/RecoCorr/Atlas/h2f_n2_ctheta_LaPLaM", - "#rho_{2} #Lambda#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_ctheta_LaMLaP", - "#rho_{2} #bar{#Lambda}#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_ctheta_LaPLaP", - "#rho_{2} #Lambda#Lambda [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - histos.add("ME/RecoCorr/Atlas/h2f_n2_ctheta_LaMLaM", - "#rho_{2} #bar{#Lambda}#bar{#Lambda} [Atlas]", kTHnSparseF, - {axisCent, axisDRap, axisDPhi, axisCosTS}); - - histos.addClone("SE/RecoCorr/Star/", "SE/RecoCorrBkgSigSB/Star/"); - histos.addClone("SE/RecoCorr/Atlas/", "SE/RecoCorrBkgSigSB/Atlas/"); - histos.addClone("SE/RecoCorr/Star/", "SE/RecoCorrBkgSBSB/Star/"); - histos.addClone("SE/RecoCorr/Atlas/", "SE/RecoCorrBkgSBSB/Atlas/"); - - histos.addClone("ME/RecoCorr/Star/", "ME/RecoCorrBkgSigSB/Star/"); - histos.addClone("ME/RecoCorr/Atlas/", "ME/RecoCorrBkgSigSB/Atlas/"); - histos.addClone("ME/RecoCorr/Star/", "ME/RecoCorrBkgSBSB/Star/"); - histos.addClone("ME/RecoCorr/Atlas/", "ME/RecoCorrBkgSBSB/Atlas/"); + if (cApplyKinWeights) { + if (cKinWeightCcdbPath.value.empty()) { + LOGF(fatal, "cApplyKinWeights is set but cKinWeightCcdbPath is empty"); + } + ccdb->setURL(cKinWeightCcdbUrl.value); + ccdb->setCaching(true); + ccdb->setFatalWhenNull(false); + } } - bool isSignal(float m) const + bool isAcceptedMass(float m) const { - return (m >= cSigMinLambda.value && m <= cSigMaxLambda.value); + return m >= cMassAccMin.value && m < cMassAccMax.value; } - bool isSideband(float m) const + + static std::vector axisEdges(AxisSpec const& axis) { - return ((m >= cSbLeftMin.value && m <= cSbLeftMax.value) || - (m >= cSbRightMin.value && m <= cSbRightMax.value)); + if (!axis.nBins.has_value()) { + return axis.binEdges; + } + const int n = axis.nBins.value(); + std::vector edges(n + 1); + for (int i = 0; i <= n; ++i) { + edges[i] = axis.binEdges[0] + (axis.binEdges[1] - axis.binEdges[0]) * i / n; + } + return edges; } - bool isInclusive(float m) const + + static constexpr double MassEdgeTolerance = 1e-3; + [[nodiscard]] bool isOnMassBinEdge(double m) const { - return (m >= cMassHistMin.value && m <= cMassHistMax.value); + const double width = (static_cast(cMassHistMax.value) - cMassHistMin.value) / cNMassBins.value; + const double k = (m - cMassHistMin.value) / width; + return std::abs(k - std::round(k)) < MassEdgeTolerance; } - void getBoostVector(std::array const& p, std::array& v, - bool inverseBoostFlag = true) + void getBoostVector(std::array const& p, std::array& v) { int n = p.size(); for (int i = 0; i < n - 1; ++i) { - if (inverseBoostFlag) { - v[i] = -p[i] / RecoDecay::e(p[0], p[1], p[2], p[3]); - } else { - v[i] = p[i] / RecoDecay::e(p[0], p[1], p[2], p[3]); - } + v[i] = -p[i] / RecoDecay::e(p[0], p[1], p[2], p[3]); } } @@ -1885,9 +1824,7 @@ struct LambdaSpinPolarization { p[2] = p[2] + gam2 * bp * b[2] + gam * b[2] * e; } - // cos of the opening angle between the two proton three-momenta. - static float cosOpeningAngle(std::array const& a, - std::array const& b) + static float cosOpeningAngle(std::array const& a, std::array const& b) { std::array n1 = {a[0], a[1], a[2]}, n2 = {b[0], b[1], b[2]}; return RecoDecay::dotProd(n1, n2) / @@ -1895,563 +1832,492 @@ struct LambdaSpinPolarization { RecoDecay::sqrtSumOfSquares(n2[0], n2[1], n2[2])); } - float cosThetaStarAtlas(std::array const& l1, - std::array const& l2, - std::array const& pr1, - std::array const& pr2) + float cosThetaStarAtlas(std::array const& l1, std::array const& l2, std::array const& pr1, std::array const& pr2) { auto l1a = l1; auto l2a = l2; auto pr1a = pr1; auto pr2a = pr2; - const float mPair = - RecoDecay::m(std::array{std::array{l1a[0], l1a[1], l1a[2]}, - std::array{l2a[0], l2a[1], l2a[2]}}, - std::array{l1a[3], l2a[3]}); - std::array llpair = {l1a[0] + l2a[0], l1a[1] + l2a[1], - l1a[2] + l2a[2], mPair}; + const float mPair = RecoDecay::m(std::array{std::array{l1a[0], l1a[1], l1a[2]}, std::array{l2a[0], l2a[1], l2a[2]}}, std::array{l1a[3], l2a[3]}); + std::array llpair = {l1a[0] + l2a[0], l1a[1] + l2a[1], l1a[2] + l2a[2], mPair}; std::array vPair{}; - getBoostVector(llpair, vPair, cInvBoostFlag); + getBoostVector(llpair, vPair); boost(l1a, vPair); boost(l2a, vPair); boost(pr1a, vPair); boost(pr2a, vPair); std::array v1p{}, v2p{}; - getBoostVector(l1a, v1p, cInvBoostFlag); - getBoostVector(l2a, v2p, cInvBoostFlag); + getBoostVector(l1a, v1p); + getBoostVector(l2a, v2p); boost(pr1a, v1p); boost(pr2a, v2p); return cosOpeningAngle(pr1a, pr2a); } - float cosThetaStarStar(std::array const& l1, - std::array const& l2, - std::array const& pr1, - std::array const& pr2) + static float deltaR(PoolTrack const& p1, PoolTrack const& p2) { - auto pr1s = pr1; - auto pr2s = pr2; - std::array v1lab{}, v2lab{}; - getBoostVector(l1, v1lab, cInvBoostFlag); - getBoostVector(l2, v2lab, cInvBoostFlag); - boost(pr1s, v1lab); - boost(pr2s, v2lab); - return cosOpeningAngle(pr1s, pr2s); + const float drap = p1.rap() - p2.rap(); + const float dphi = RecoDecay::constrainAngle(p1.phi() - p2.phi(), -PI); + return std::sqrt(drap * drap + dphi * dphi); } - template - void fillPairHistos(U const& p1, U const& p2) + template + void fillPair(PoolTrack const& p1, PoolTrack const& p2, float dR, float w, float mReplaced = 0.f) { - static constexpr std::array SubDir = { - "LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; - - constexpr bool IsSigSB = - (part_pair == kLambdaSBAntiLambda || part_pair == kAntiLambdaSBLambda || - part_pair == kLambdaSBLambda || part_pair == kAntiLambdaSBAntiLambda); - constexpr bool IsSBSB = - (part_pair == kLambdaSBSBAntiLambda || - part_pair == kAntiLambdaSBSBLambda || part_pair == kLambdaSBSBLambda || - part_pair == kAntiLambdaSBSBAntiLambda); - - constexpr int Idx = IsSigSB ? static_cast(part_pair) - 4 - : IsSBSB ? static_cast(part_pair) - 8 - : static_cast(part_pair); - - 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}; + static constexpr std::array TierDir = {"SE/hPair_", "SEproc/hPair_", "ME/hPair_", + "McGen/SE/hPair_", "McGen/SEproc/hPair_", "McGen/ME/hPair_"}; + static constexpr std::array PairTag = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; + + std::array l1 = {p1.px(), p1.py(), p1.pz(), p1.mass()}; + std::array l2 = {p2.px(), p2.py(), p2.pz(), p2.mass()}; std::array pr1 = {p1.prPx(), p1.prPy(), p1.prPz(), MassProton}; std::array pr2 = {p2.prPx(), p2.prPy(), p2.prPz(), MassProton}; - if (cDoAtlasMethod) { - float ctheta = cosThetaStarAtlas(l1, l2, pr1, pr2); - if constexpr (!IsSigSB && !IsSBSB) { - histos.fill(HIST("SE/Reco/Atlas/h2f_n2_mass_") + HIST(SubDir[Idx]), - cent, p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, - dR, ctheta); - histos.fill(HIST("SE/RecoCorr/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorr/Atlas/h2f_n2_dltaR_") + HIST(SubDir[Idx]), - cent, dR, ctheta); - } else if constexpr (IsSigSB) { - histos.fill(HIST("SE/RecoBkgSigSB/Atlas/h2f_n2_mass_") + - HIST(SubDir[Idx]), - cent, p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, - dR, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSigSB/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSigSB/Atlas/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta); - } else { - histos.fill( - HIST("SE/RecoBkgSBSB/Atlas/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSBSB/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSBSB/Atlas/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta); - } - } - - if (cDoStarMethod) { - float ctheta = cosThetaStarStar(l1, l2, pr1, pr2); - if constexpr (!IsSigSB && !IsSBSB) { - histos.fill(HIST("SE/Reco/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, - ctheta); - histos.fill(HIST("SE/RecoCorr/Star/h2f_n2_ctheta_") + HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorr/Star/h2f_n2_dltaR_") + HIST(SubDir[Idx]), - cent, dR, ctheta); - } else if constexpr (IsSigSB) { - histos.fill( - HIST("SE/RecoBkgSigSB/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSigSB/Star/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSigSB/Star/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta); - } else { - histos.fill( - HIST("SE/RecoBkgSBSB/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSBSB/Star/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta); - histos.fill(HIST("SE/RecoCorrBkgSBSB/Star/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta); - } + const float ctheta = cosThetaStarAtlas(l1, l2, pr1, pr2); + if constexpr (tier == kMixedEvent || tier == kMcGenMixedEvent) { + histos.fill(HIST(TierDir[tier]) + HIST(PairTag[part_pair]), p1.mass(), p2.mass(), dR, ctheta, mReplaced, w); + } else { + histos.fill(HIST(TierDir[tier]) + HIST(PairTag[part_pair]), p1.mass(), p2.mass(), dR, ctheta, w); } } - template - void fillPairHistosWeighted(PoolTrack const& p1, PoolTrack const& p2, - float w) + template + void fillLeg2Spectrum(PoolTrack const& leg, float mReplaced, float w) { - static constexpr std::array SubDir = { - "LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; - - constexpr bool IsSigSB = - (part_pair == kLambdaSBAntiLambda || part_pair == kAntiLambdaSBLambda || - part_pair == kLambdaSBLambda || part_pair == kAntiLambdaSBAntiLambda); - constexpr bool IsSBSB = - (part_pair == kLambdaSBSBAntiLambda || - part_pair == kAntiLambdaSBSBLambda || part_pair == kLambdaSBSBLambda || - part_pair == kAntiLambdaSBSBAntiLambda); - - constexpr int Idx = IsSigSB ? static_cast(part_pair) - 4 - : IsSBSB ? static_cast(part_pair) - 8 - : static_cast(part_pair); - - 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 pr1 = {p1.prPx(), p1.prPy(), p1.prPz(), MassProton}; - std::array pr2 = {p2.prPx(), p2.prPy(), p2.prPz(), MassProton}; + static constexpr std::array Dir = {"KinW/SEproc/hLeg2_", "KinW/ME/hLeg2_"}; + static constexpr std::array PairTag = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; + histos.fill(HIST(Dir[IsME]) + HIST(PairTag[part_pair]), leg.pt(), leg.rap(), leg.phi(), mReplaced, w); + } - if (cDoAtlasMethod) { - float ctheta = cosThetaStarAtlas(l1, l2, pr1, pr2); - - if constexpr (!IsSigSB && !IsSBSB) { - histos.fill(HIST("ME/Reco/Atlas/h2f_n2_mass_") + HIST(SubDir[Idx]), - cent, p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, - dR, ctheta, w); - histos.fill(HIST("ME/RecoCorr/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorr/Atlas/h2f_n2_dltaR_") + HIST(SubDir[Idx]), - cent, dR, ctheta, w); - } else if constexpr (IsSigSB) { - histos.fill(HIST("ME/RecoBkgSigSB/Atlas/h2f_n2_mass_") + - HIST(SubDir[Idx]), - cent, p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, - dR, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSigSB/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSigSB/Atlas/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta, w); - } else { - histos.fill( - HIST("ME/RecoBkgSBSB/Atlas/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSBSB/Atlas/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSBSB/Atlas/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta, w); - } - } - - if (cDoStarMethod) { - float ctheta = cosThetaStarStar(l1, l2, pr1, pr2); - - if constexpr (!IsSigSB && !IsSBSB) { - histos.fill(HIST("ME/Reco/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, - ctheta, w); - histos.fill(HIST("ME/RecoCorr/Star/h2f_n2_ctheta_") + HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorr/Star/h2f_n2_dltaR_") + HIST(SubDir[Idx]), - cent, dR, ctheta, w); - } else if constexpr (IsSigSB) { - histos.fill( - HIST("ME/RecoBkgSigSB/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSigSB/Star/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSigSB/Star/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta, w); - } else { - histos.fill( - HIST("ME/RecoBkgSBSB/Star/h2f_n2_mass_") + HIST(SubDir[Idx]), cent, - p1.mass(), p2.mass(), p1.pt(), p2.pt(), drap, dphi, dR, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSBSB/Star/h2f_n2_ctheta_") + - HIST(SubDir[Idx]), - cent, drap, dphi, ctheta, w); - histos.fill(HIST("ME/RecoCorrBkgSBSB/Star/h2f_n2_dltaR_") + - HIST(SubDir[Idx]), - cent, dR, ctheta, w); + [[nodiscard]] bool isSameWindow(float m1, float m2) const + { + const int b = replacedMassAxis.FindFixBin(m1); + return b >= 1 && b <= replacedMassAxis.GetNbins() && b == replacedMassAxis.FindFixBin(m2); + } + + template + void fillDaughterSeparation(std::array const& a, double ma, std::array const& b, double mb, float dR) + { + static_assert(tier == kSameEvent || tier == kSameEventInMixing, "close-pair QA is filled for SE and SEproc pairs only"); + static constexpr std::array TierDir = {"QA/ClosePair/", "QA/ClosePairSEproc/"}; + static constexpr std::array ClsDir = {"LamLam/", "ALamALam/", "UnlikeSign/"}; + static constexpr std::array ComboName = {"h_pp_", "h_pipi_", "h_p1pi2_", "h_pi1p2_"}; + const float dEta = std::asinh(a[2] / std::hypot(a[0], a[1])) - std::asinh(b[2] / std::hypot(b[0], b[1])); + const float dY = RecoDecay::y(a, ma) - RecoDecay::y(b, mb); + const float dPhi = RecoDecay::constrainAngle(std::atan2(a[1], a[0]) - std::atan2(b[1], b[0]), -PI); + histos.fill(HIST(TierDir[tier]) + HIST(ClsDir[Cls]) + HIST(ComboName[Combo]) + HIST("dEta"), dR, dEta, dPhi); + histos.fill(HIST(TierDir[tier]) + HIST(ClsDir[Cls]) + HIST(ComboName[Combo]) + HIST("dY"), dR, dY, dPhi); + } + + template + void fillClosePairQA(PoolTrack const& l1, PoolTrack const& l2, float dR) + { + constexpr int Cls = part_pair == kLambdaLambda ? 0 : (part_pair == kAntiLambdaAntiLambda ? 1 : 2); + const std::array pr1 = {l1.prPx(), l1.prPy(), l1.prPz()}; + const std::array pr2 = {l2.prPx(), l2.prPy(), l2.prPz()}; + const std::array pi1 = {l1.px() - l1.prPx(), l1.py() - l1.prPy(), l1.pz() - l1.prPz()}; + const std::array pi2 = {l2.px() - l2.prPx(), l2.py() - l2.prPy(), l2.pz() - l2.prPz()}; + fillDaughterSeparation(pr1, MassProton, pr2, MassProton, dR); + fillDaughterSeparation(pi1, MassPionCharged, pi2, MassPionCharged, dR); + fillDaughterSeparation(pr1, MassProton, pi2, MassPionCharged, dR); + fillDaughterSeparation(pi1, MassPionCharged, pr2, MassProton, dR); + } + + bool isClosePairQAMass(float m) const + { + return cClosePairMassWindow.value <= 0.f || std::abs(m - MassLambda0) < cClosePairMassWindow.value; + } + + bool isKinematicMatch(PoolTrack const& cand, PoolTrack const& leg) const + { + return std::abs(cand.pt() - leg.pt()) < cMaxDeltaPt.value && + std::abs(cand.rap() - leg.rap()) < cMaxDeltaRap.value && + std::abs(RecoDecay::constrainAngle(cand.phi() - leg.phi(), -PI)) < cMaxDeltaPhi.value; + } + + static constexpr std::size_t NLegAxes = 4; + template + float kinWeight(PoolTrack const& cand, float mReplaced) const + { + THnBase const* map = kinWeightMaps[part_pair].get(); + if (!map) { + return 1.f; + } + const std::array x = {cand.pt(), cand.rap(), cand.phi(), mReplaced}; + std::array idx{}; + for (std::size_t i = 0; i < NLegAxes; ++i) { + TAxis const* axis = map->GetAxis(i); + idx[i] = axis->FindFixBin(x[i]); + if (idx[i] < 1 || idx[i] > axis->GetNbins()) { + return 1.f; } } + const int64_t bin = map->GetBin(idx.data()); + if (bin < 0) { + return 1.f; + } + const double w = map->GetBinContent(bin); + return w > 0. ? static_cast(w) : 1.f; } - template - void analyzePairsWithMassWindow(T const& trks_1, T const& trks_2) + static constexpr double AxisEdgeTolerance = 1e-6; + [[nodiscard]] bool hasLegAxes(THnBase const& map) const { - for (auto const& trk_1 : trks_1) { - if (!isInclusive(trk_1.mass())) { - continue; + if (map.GetNdimensions() != static_cast(legAxesRef.size())) { + return false; + } + for (std::size_t i = 0; i < legAxesRef.size(); ++i) { + TAxis const* axis = map.GetAxis(i); + TAxis const& ref = legAxesRef[i]; + if (axis->GetNbins() != ref.GetNbins()) { + return false; } - bool t1sig = isSignal(trk_1.mass()); - bool t1sb = isSideband(trk_1.mass()); - - for (auto const& trk_2 : trks_2) { - if constexpr (samelambda) { - if (trk_1.index() == trk_2.index()) { - continue; - } - } - if (!isInclusive(trk_2.mass())) { - continue; + for (int b = 1; b <= ref.GetNbins() + 1; ++b) { + if (std::abs(axis->GetBinLowEdge(b) - ref.GetBinLowEdge(b)) > AxisEdgeTolerance) { + return false; } - bool t2sig = isSignal(trk_2.mass()); - bool t2sb = isSideband(trk_2.mass()); + } + } + return true; + } - fillPairHistos(trk_1, trk_2); - if ((t1sig && t2sb) || (t1sb && t2sig)) { - fillPairHistos(trk_1, trk_2); - } - if (t1sb && t2sb) { - fillPairHistos(trk_1, trk_2); - } + void loadKinWeights(uint64_t collisionTimestamp) + { + static constexpr std::array PairTag = {"LaPLaM", "LaMLaP", "LaPLaP", "LaMLaM"}; + const int64_t ts = cKinWeightCcdbTimestamp.value >= 0 ? cKinWeightCcdbTimestamp.value : static_cast(collisionTimestamp); + auto* list = ccdb->getForTimeStamp(cKinWeightCcdbPath.value, ts); + if (!list) { + LOGP(fatal, "Kinematic weights: no object at {} for timestamp {}", cKinWeightCcdbPath.value, ts); + return; + } + auto const* guard = dynamic_cast(list->FindObject("kinWeightGuardStatus")); + if (guard && !TString(guard->GetTitle()).BeginsWith("OK")) { + LOGF(fatal, "Kinematic weights at %s are flagged by their builder: %s", cKinWeightCcdbPath.value.c_str(), guard->GetTitle()); + return; + } + for (std::size_t i = 0; i < PairTag.size(); ++i) { + const std::string name = "hKinW_" + std::string(PairTag[i]); + auto* map = dynamic_cast(list->FindObject(name.c_str())); + if (!map || !hasLegAxes(*map)) { + LOGF(fatal, "Kinematic weights: %s is missing, or its axes differ from KinW/*/hLeg2_%s (pT, y, phi, window of the replaced leg) of this configuration", name.c_str(), std::string(PairTag[i]).c_str()); + return; } + kinWeightMaps[i].reset(static_cast(map->Clone())); } + auto const* form = dynamic_cast(list->FindObject("kinWeightForm")); + kinWeightsLoaded = true; + LOGF(info, "Kinematic weights loaded from %s (form: %s)", cKinWeightCcdbPath.value.c_str(), form ? form->GetTitle() : "not stated"); } - static constexpr int MEModeStandard = 0; - static constexpr int MEModeKinematic = 1; + static uint64_t splitMix(uint64_t x) + { + x += 0x9E3779B97F4A7C15ULL; + x = (x ^ (x >> 30)) * 0xBF58476D1CE4E5B9ULL; + x = (x ^ (x >> 27)) * 0x94D049BB133111EBULL; + return x ^ (x >> 31); + } - template - void analyzePairsME(T const& trks_1, T const& trks_2) + template + void setPairOrderSeed(C const& col) { + const uint64_t vz = std::bit_cast(col.posZ()); + pairOrderSeed = splitMix(((vz << 32) | static_cast(col.globalIndex())) ^ (col.timeStamp() * 0xD6E8FEB86659FD93ULL)); + } + + template + bool isCountedOrder(T const& trk_1, T const& trk_2) const + { + const int64_t g1 = trk_1.globalIndex(), g2 = trk_2.globalIndex(); + if (g1 == g2) { + return false; + } + const bool firstIsLo = g1 < g2; + const uint64_t lo = static_cast(firstIsLo ? g1 : g2); + const uint64_t hi = static_cast(firstIsLo ? g2 : g1); + const bool swapped = splitMix(pairOrderSeed ^ ((lo << 32) | hi)) & 1ULL; + return firstIsLo != swapped; + } + template + void analyzePairsSE(T const& trks_1, T const& trks_2) + { for (auto const& trk_1 : trks_1) { - if (!isInclusive(trk_1.mass())) { + if (!isAcceptedMass(trk_1.mass())) { continue; } - const bool t1sig = isSignal(trk_1.mass()); - const bool t1sb = isSideband(trk_1.mass()); - PoolTrack p1 = toPoolTrack(trk_1); + const PoolTrack p1 = toPoolTrack(trk_1); for (auto const& trk_2 : trks_2) { - if constexpr (samelambda) { - if (trk_1.index() == trk_2.index()) { - continue; - } - } - if (!isInclusive(trk_2.mass())) { + if (!isCountedOrder(trk_1, trk_2) || !isAcceptedMass(trk_2.mass())) { continue; } - const bool t2sig = isSignal(trk_2.mass()); - const bool t2sb = isSideband(trk_2.mass()); - PoolTrack p2 = toPoolTrack(trk_2); - - fillPairHistosWeighted(p1, p2, 1.0f); - - if ((t1sig && t2sb) || (t1sb && t2sig)) { - fillPairHistosWeighted(p1, p2, 1.0f); - } - - if (t1sb && t2sb) { - fillPairHistosWeighted(p1, p2, 1.0f); + const PoolTrack p2 = toPoolTrack(trk_2); + const float dR = deltaR(p1, p2); + fillPair(p1, p2, dR, 1.f); + if constexpr (tier == kSameEvent) { + if (cFillClosePairQA && isClosePairQAMass(p1.mass()) && isClosePairQAMass(p2.mass())) { + fillClosePairQA(p1, p2, dR); + } } } } } - template - void analyzePairsMEKinematic(T const& se_trks_1, T const& se_trks_2, - T const& me_pool_1, T const& me_pool_2) + template + void fillEventQA(float centrality, T const& lTrks, T const& alTrks) { - std::vector meVec1, meVec2; - meVec1.reserve(me_pool_1.size()); - meVec2.reserve(me_pool_2.size()); - for (auto const& me : me_pool_1) { - if (isInclusive(me.mass())) { - meVec1.push_back(toPoolTrack(me)); + static constexpr std::array EvDir = {"Events/", "McGen/Events/"}; + std::size_t nLambda = 0, nAntiLambda = 0; + for (auto const& trk : lTrks) { + if (isAcceptedMass(trk.mass())) { + ++nLambda; } } - for (auto const& me : me_pool_2) { - if (isInclusive(me.mass())) { - meVec2.push_back(toPoolTrack(me)); + for (auto const& trk : alTrks) { + if (isAcceptedMass(trk.mass())) { + ++nAntiLambda; } } - - if (meVec1.empty() && meVec2.empty()) { + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeAll); + histos.fill(HIST(EvDir[Gen]) + HIST("hNCandAccepted"), nLambda, nAntiLambda); + if (nLambda + nAntiLambda >= 1) { + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeOneCand); + } + if (nLambda + nAntiLambda < NCandidatesForPair) { return; } + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeTwoCand); + histos.fill(HIST(EvDir[Gen]) + HIST("hCentPairEvents"), centrality); + if (nLambda >= 1 && nAntiLambda >= 1) { + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeUnlikeSign); + } + if (nLambda >= NCandidatesForPair) { + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeLambdaLambda); + } + if (nAntiLambda >= NCandidatesForPair) { + histos.fill(HIST(EvDir[Gen]) + HIST("hEventCount"), kPeAntiLambdaAntiLambda); + } + for (auto const& trk : lTrks) { + if (isAcceptedMass(trk.mass())) { + fillPairCandQA(trk); + } + } + for (auto const& trk : alTrks) { + if (isAcceptedMass(trk.mass())) { + fillPairCandQA(trk); + } + } + } + template + void fillPairCandQA(T const& trk) + { + static constexpr std::array CandDir = {"QA/PairCand/", "McGen/QA/PairCand/"}; + static constexpr std::array Species = {"Lambda/", "AntiLambda/"}; + histos.fill(HIST(CandDir[Gen]) + HIST(Species[part]) + HIST("hMassVsPt"), trk.mass(), trk.pt()); + histos.fill(HIST(CandDir[Gen]) + HIST(Species[part]) + HIST("hRapVsPhi"), trk.rap(), trk.phi()); + } + + [[nodiscard]] int poolBin(float centrality, float vz) const + { + const int iCent = mePoolAxisCent.FindFixBin(centrality); + const int iVz = mePoolAxisVz.FindFixBin(vz); + if (iCent < 1 || iCent > mePoolAxisCent.GetNbins() || iVz < 1 || iVz > mePoolAxisVz.GetNbins()) { + return -1; + } + return (iCent - 1) * mePoolAxisVz.GetNbins() + (iVz - 1); + } + + template + void analyzePairsMEKinematic(T const& se_trks_1, T const& se_trks_2, RollingPool const& pool) + { + static constexpr std::array MixQaDir = {"QA/ME/", "McGen/QA/ME/"}; + constexpr PairTier SeTier = Gen ? kMcGenSameEventInMixing : kSameEventInMixing; + constexpr PairTier MeTier = Gen ? kMcGenMixedEvent : kMixedEvent; + std::vector matches; for (auto const& trk1 : se_trks_1) { - if (!isInclusive(trk1.mass())) { + if (!isAcceptedMass(trk1.mass())) { continue; } - const bool se1sig = isSignal(trk1.mass()); - const bool se1sb = isSideband(trk1.mass()); - PoolTrack p1 = toPoolTrack(trk1); + const PoolTrack p1 = toPoolTrack(trk1); for (auto const& trk2 : se_trks_2) { - if constexpr (samelambda) { - if (trk1.index() == trk2.index()) { - continue; - } - } - if (!isInclusive(trk2.mass())) { + if (!isCountedOrder(trk1, trk2) || !isAcceptedMass(trk2.mass())) { continue; } - const bool se2sig = isSignal(trk2.mass()); - const bool se2sb = isSideband(trk2.mass()); - PoolTrack p2 = toPoolTrack(trk2); - - { - std::vector matchA; - for (auto const& meP : meVec2) { - if (std::abs(meP.pt() - p2.pt()) < cMaxDeltaPt && - std::abs(meP.rap() - p2.rap()) < cMaxDeltaRap && - std::abs(RecoDecay::constrainAngle(meP.phi() - p2.phi(), -PI)) < - cMaxDeltaPhi) { - matchA.push_back(meP); - } - } - if (!matchA.empty()) { - const float wA = 1.0f / static_cast(matchA.size()); - for (auto const& meP2 : matchA) { - const bool me2sig = isSignal(meP2.mass()); - const bool me2sb = isSideband(meP2.mass()); - - fillPairHistosWeighted(p1, meP2, wA); - - if ((se1sig && me2sb) || (se1sb && me2sig)) { - fillPairHistosWeighted(p1, meP2, wA); - } + const PoolTrack p2 = toPoolTrack(trk2); + const float dRSeed = deltaR(p1, p2); - if (se1sb && me2sb) { - fillPairHistosWeighted(p1, meP2, wA); - } - } + if (cFillSEInMixing) { + fillPair(p1, p2, dRSeed, 1.f); + } + if constexpr (!Gen) { + if (cFillClosePairQA && isClosePairQAMass(p1.mass()) && isClosePairQAMass(p2.mass())) { + fillClosePairQA(p1, p2, dRSeed); + } + if (cFillKinWSpectra) { + fillLeg2Spectrum(p2, p2.mass(), 1.f); } } - { - std::vector matchB; - for (auto const& meP : meVec1) { - if (std::abs(meP.pt() - p1.pt()) < cMaxDeltaPt && - std::abs(meP.rap() - p1.rap()) < cMaxDeltaRap && - std::abs(RecoDecay::constrainAngle(meP.phi() - p1.phi(), -PI)) < - cMaxDeltaPhi) { - matchB.push_back(meP); + pool.findMatches(p2, [this, &p2](PoolTrack const& cand) { return isKinematicMatch(cand, p2); }, matches); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hMatchesLeg2"), matches.size(), part_pair); + for (auto const& meP2 : matches) { + float kinW = 1.f; + if constexpr (!Gen) { + kinW = kinWeight(*meP2, p2.mass()); + } + const float w = kinW / static_cast(matches.size()); + fillPair(p1, *meP2, dRSeed, w, p2.mass()); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hDeltaRShift"), dRSeed, deltaR(p1, *meP2) - dRSeed, w); + if constexpr (!Gen) { + if (cFillKinWSpectra && isSameWindow(meP2->mass(), p2.mass())) { + fillLeg2Spectrum(*meP2, p2.mass(), w); } } - if (!matchB.empty()) { - const float wB = 1.0f / static_cast(matchB.size()); - for (auto const& meP1 : matchB) { - const bool me1sig = isSignal(meP1.mass()); - const bool me1sb = isSideband(meP1.mass()); + } + } + } + } - fillPairHistosWeighted(meP1, p2, wB); + template + void mixCollision(C const& col, T const& lTrks, T const& alTrks, std::vector& pools, std::vector& eventCands) + { + static constexpr std::array MixQaDir = {"QA/ME/", "McGen/QA/ME/"}; + histos.fill(HIST(MixQaDir[Gen]) + HIST("hEventStatus"), kMEEventSeen); + const int bin = poolBin(col.cent(), col.posZ()); + if (bin < 0) { + histos.fill(HIST(MixQaDir[Gen]) + HIST("hEventStatus"), kMEEventOutsidePools); + return; + } + histos.fill(HIST(MixQaDir[Gen]) + HIST("hLambdaMultVsCent"), col.cent(), lTrks.size()); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hAntiLambdaMultVsCent"), col.cent(), alTrks.size()); - if ((me1sig && se2sb) || (me1sb && se2sig)) { - fillPairHistosWeighted(meP1, p2, wB); - } + eventCands.clear(); + for (auto const& trk : lTrks) { + if (isAcceptedMass(trk.mass())) { + eventCands.push_back(toPoolTrack(trk)); + } + } + for (auto const& trk : alTrks) { + if (isAcceptedMass(trk.mass())) { + eventCands.push_back(toPoolTrack(trk)); + } + } - if (me1sb && se2sb) { - fillPairHistosWeighted(meP1, p2, wB); - } - } + RollingPool& pool = pools[bin]; + if (eventCands.size() >= NCandidatesForPair) { + if (static_cast(pool.nEvents()) >= cMEPoolMinEvents.value) { + if constexpr (!Gen) { + if (cApplyKinWeights && !kinWeightsLoaded) { + loadKinWeights(col.timeStamp()); } } + setPairOrderSeed(col); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hEventStatus"), kMEEventMixed); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hMixedCentVz"), col.cent(), col.posZ()); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hPoolEventsAtMixing"), pool.nEvents()); + analyzePairsMEKinematic(lTrks, alTrks, pool); + analyzePairsMEKinematic(alTrks, lTrks, pool); + analyzePairsMEKinematic(lTrks, lTrks, pool); + analyzePairsMEKinematic(alTrks, alTrks, pool); + } else { + histos.fill(HIST(MixQaDir[Gen]) + HIST("hEventStatus"), kMEEventWarmUp); + } + } + if (static_cast(eventCands.size()) >= cMEPoolMinCand.value) { + pool.push(eventCands); + histos.fill(HIST(MixQaDir[Gen]) + HIST("hEventStatus"), kMEEventDonated); + } + } - } // trk2 - } // trk1 - } // analyzePairsMEKinematic - - // ========================================================================= - // O2 framework declarations - // ========================================================================= using LambdaCollisions = aod::LambdaCollisions; using LambdaMcGenCollisions = aod::LambdaMcGenCollisions; using LambdaTracks = soa::Join; - Preslice perCollisionLambda = - aod::lambdatrack::lambdaCollisionId; + 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 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); + Partition partAntiLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kAntiLambda) && (aod::lambdatrackext::trueLambdaFlag == true) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); SliceCache cachemc; - Partition partMcLambdaTracks = - (aod::lambdatrack::v0Type == (int8_t)kLambda) && - (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); + Partition partMcLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kLambda) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); - Partition partMcAntiLambdaTracks = - (aod::lambdatrack::v0Type == (int8_t)kAntiLambda) && - (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); + Partition partMcAntiLambdaTracks = (aod::lambdatrack::v0Type == (int8_t)kAntiLambda) && (aod::lambdatrack::v0PrmScd == (int8_t)kPrimary); void processDummy(LambdaCollisions::iterator const&) {} PROCESS_SWITCH(LambdaSpinPolarization, processDummy, "Dummy", false); - void processDataReco(LambdaCollisions::iterator const& collision, - LambdaTracks const&) + void processDataReco(LambdaCollisions::iterator const& collision, LambdaTracks const&) { - cent = collision.cent(); - auto lTrks = partLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); - auto alTrks = partAntiLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); - - analyzePairsWithMassWindow(lTrks, alTrks); - analyzePairsWithMassWindow(alTrks, lTrks); - analyzePairsWithMassWindow(lTrks, lTrks); - analyzePairsWithMassWindow(alTrks, alTrks); + setPairOrderSeed(collision); + auto lTrks = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto alTrks = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + if (!doprocessDataRecoMixed) { + fillEventQA(collision.cent(), lTrks, alTrks); + } + + analyzePairsSE(lTrks, alTrks); + analyzePairsSE(alTrks, lTrks); + analyzePairsSE(lTrks, lTrks); + analyzePairsSE(alTrks, alTrks); } - PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, - "SE only (data/MCReco)", false); + PROCESS_SWITCH(LambdaSpinPolarization, processDataReco, "SE only (data/MCReco)", false); - struct GetMultiplicity { - float operator()(auto const& col) const { return col.cent(); } - }; - using MixedBinning = - FlexibleBinningPolicy, - o2::aod::collision::PosZ, GetMultiplicity>; - MixedBinning binningOnVtxAndMult{ - {GetMultiplicity{}}, - {axisVtxZME, axisCentME}, - true}; - - void processDataRecoMixed(LambdaCollisions const& col, LambdaTracks const&) + void processDataRecoMixed(LambdaCollisions const& cols, LambdaTracks const&) { - for (auto const& [col1, col2] : soa::selfCombinations( - binningOnVtxAndMult, mixingParameter, -1, col, col)) { - if (col1.globalIndex() == col2.globalIndex()) { - continue; - } - cent = col1.cent(); - histos.fill(HIST("QA/ME/hPoolCentVz"), col1.cent(), col1.posZ()); - - auto lTrks1 = partLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); - auto lTrks2 = partLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); - auto alTrks1 = partAntiLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, col1.globalIndex(), cache); - auto alTrks2 = partAntiLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, col2.globalIndex(), cache); - - histos.fill(HIST("QA/ME/hLambdaMultVsCent"), col1.cent(), lTrks1.size()); - histos.fill(HIST("QA/ME/hAntiLambdaMultVsCent"), col1.cent(), - alTrks1.size()); - - if (cMEMode == MEModeStandard) { - analyzePairsME(lTrks1, alTrks2); - analyzePairsME(alTrks1, lTrks2); - analyzePairsME(lTrks1, lTrks2); - analyzePairsME(alTrks1, alTrks2); - - } else if (cMEMode == MEModeKinematic) { - analyzePairsMEKinematic(lTrks1, alTrks1, - lTrks2, alTrks2); - analyzePairsMEKinematic(alTrks1, lTrks1, - alTrks2, lTrks2); - analyzePairsMEKinematic(lTrks1, lTrks1, lTrks2, - lTrks2); - analyzePairsMEKinematic( - alTrks1, alTrks1, alTrks2, alTrks2); - } + std::vector eventCands; + for (auto const& col : cols) { + auto lTrks = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col.globalIndex(), cache); + auto alTrks = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, col.globalIndex(), cache); + fillEventQA(col.cent(), lTrks, alTrks); + mixCollision(col, lTrks, alTrks, mePools, eventCands); + } + } + PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixed, "ME: kinematic mixing against rolling pools that cross DataFrames", false); + + void processMcGen(LambdaMcGenCollisions::iterator const& collision, aod::LambdaMcGenTracks const&) + { + setPairOrderSeed(collision); + auto lTrks = partMcLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), cachemc); + auto alTrks = partMcAntiLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), cachemc); + if (!doprocessMcGenMixed) { + fillEventQA(collision.cent(), lTrks, alTrks); + } + + analyzePairsSE(lTrks, alTrks); + analyzePairsSE(alTrks, lTrks); + analyzePairsSE(lTrks, lTrks); + analyzePairsSE(alTrks, alTrks); + } + PROCESS_SWITCH(LambdaSpinPolarization, processMcGen, "MC generated: SE pairs (McGen/SE), as processDataReco", false); + + void processMcGenMixed(LambdaMcGenCollisions const& cols, aod::LambdaMcGenTracks const&) + { + std::vector eventCands; + for (auto const& col : cols) { + auto lTrks = partMcLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, col.globalIndex(), cachemc); + auto alTrks = partMcAntiLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, col.globalIndex(), cachemc); + fillEventQA(col.cent(), lTrks, alTrks); + mixCollision(col, lTrks, alTrks, mePoolsGen, eventCands); } } - PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixed, - "ME (modes 0=standard, 1=kinematicConstrained)", false); + PROCESS_SWITCH(LambdaSpinPolarization, processMcGenMixed, "MC generated: kinematic mixing against rolling pools, as processDataRecoMixed", false); - void processDataRecoMixEvent(LambdaCollisions::iterator const& collision, - LambdaTracks const&) + void processDataRecoMixEvent(LambdaCollisions::iterator const& collision, LambdaTracks const&) { - auto lTrks = partLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); - auto alTrks = partAntiLambdaTracks->sliceByCached( - aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto lTrks = partLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); + auto alTrks = partAntiLambdaTracks->sliceByCached(aod::lambdatrack::lambdaCollisionId, collision.globalIndex(), cache); if (lTrks.size() == 0 && alTrks.size() == 0) { return; } - lambdaMixEvtCol(collision.index(), collision.cent(), collision.posZ(), - collision.timeStamp()); + lambdaMixEvtCol(collision.index(), collision.cent(), collision.posZ(), collision.timeStamp()); for (auto const& track : lTrks) { lambdaMixEvtTrk(collision.index(), track.globalIndex(), track.px(), @@ -2466,18 +2332,12 @@ struct LambdaSpinPolarization { collision.timeStamp()); } } - PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixEvent, - "Mix-event table filling", false); + PROCESS_SWITCH(LambdaSpinPolarization, processDataRecoMixEvent, "Mix-event table filling", false); - void processMcGenMixEvent(LambdaMcGenCollisions::iterator const& collision, - aod::LambdaMcGenTracks const&) + void processMcGenMixEvent(LambdaMcGenCollisions::iterator const& collision, aod::LambdaMcGenTracks const&) { - auto lTrks = partMcLambdaTracks->sliceByCached( - aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), - cachemc); - auto alTrks = partMcAntiLambdaTracks->sliceByCached( - aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), - cachemc); + auto lTrks = partMcLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), cachemc); + auto alTrks = partMcAntiLambdaTracks->sliceByCached(aod::lambdamcgentrack::lambdaMcGenCollisionId, collision.globalIndex(), cachemc); if (lTrks.size() == 0 && alTrks.size() == 0) { return; @@ -2499,8 +2359,7 @@ struct LambdaSpinPolarization { track.v0Type(), collision.timeStamp()); } } - PROCESS_SWITCH(LambdaSpinPolarization, processMcGenMixEvent, - "Mix-event McGen table filling", false); + PROCESS_SWITCH(LambdaSpinPolarization, processMcGenMixEvent, "Mix-event McGen table filling", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) {