From ee9527881c92dab823bcbd0d152c86932472b7a4 Mon Sep 17 00:00:00 2001 From: Durgesh Bhatt Date: Tue, 4 Aug 2026 17:55:36 +0530 Subject: [PATCH] Add event loss, event splitting, and cut-flow QA --- PWGLF/Tasks/Resonances/deltaAnalysis.cxx | 806 ++++++++++++++++------- 1 file changed, 551 insertions(+), 255 deletions(-) diff --git a/PWGLF/Tasks/Resonances/deltaAnalysis.cxx b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx index a1697764db9..39c631fd471 100644 --- a/PWGLF/Tasks/Resonances/deltaAnalysis.cxx +++ b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx @@ -13,10 +13,13 @@ /// \brief Delta(1232) resonance analysis via proton-pion invariant mass reconstruction with advance PID and background rejection cuts /// \author Durgesh Bhatt +#include "PWGLF/DataModel/mcCentrality.h" + #include "Common/CCDB/EventSelectionParams.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -33,6 +36,7 @@ #include #include #include +#include #include #include @@ -41,6 +45,9 @@ #include #include #include +#include +#include +#include #include using namespace o2; @@ -52,6 +59,7 @@ namespace { constexpr float massProton = o2::constants::physics::MassProton; constexpr float massPion = o2::constants::physics::MassPionCharged; +constexpr float minAbsCharge = 1e-3f; } // namespace namespace delta_analysis { @@ -73,89 +81,106 @@ struct DeltaAnalysis { SliceCache cache; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - // FIX name/configurable: single space between type and name; member name == JSON key; lowerCamelCase - Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted |z-vertex| range [cm]"}; - - Configurable cfgUseNoSameBunchPileupCut{"cfgUseNoSameBunchPileupCut", false, "Apply kNoSameBunchPileup event selection"}; - - Configurable cfgUseGoodZvtxFT0vsPVCut{"cfgUseGoodZvtxFT0vsPVCut", false, "Apply kIsGoodZvtxFT0vsPV event selection"}; - Configurable applyOccupancyInTimeRangeCut{"applyOccupancyInTimeRangeCut", false, "Apply occupancy-in-time-range cut"}; - Configurable cfgOccupancyMin{"cfgOccupancyMin", 0, "Minimum track occupancy in time range"}; - Configurable cfgOccupancyMax{"cfgOccupancyMax", 9999, "Maximum track occupancy in time range"}; - - Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator: 0=FT0M 1=FT0A 2=FT0C 3=FV0A 4=NTPV"}; - Configurable cfgCentMin{"cfgCentMin", 0.f, "Minimum centrality percentile"}; - Configurable cfgCentMax{"cfgCentMax", 100.f, "Maximum centrality percentile"}; - - Configurable cfgCutPt{"cfgCutPt", 0.2f, "Minimum pT of daughter track [GeV/c]"}; - Configurable cfgCutEta{"cfgCutEta", 0.8f, "Maximum |eta| of daughter track"}; - // CHANGED: replaced single symmetric cfgCutY with independent cfgMinY / cfgMaxY - Configurable cfgMinY{"cfgMinY", -0.5f, "Minimum rapidity of reconstructed Delta"}; - Configurable cfgMaxY{"cfgMaxY", 0.5f, "Maximum rapidity of reconstructed Delta"}; - - Configurable cfgMinITSClusters{"cfgMinITSClusters", 5, "Minimum ITS clusters"}; - Configurable cfgMinTPCClusters{"cfgMinTPCClusters", 70, "Minimum TPC clusters found"}; - Configurable numberOfInvMassBins{"numberOfInvMassBins", 120, "Number of bins along the invariant mass axis"}; - Configurable cfgMinTPCCrossedRows{"cfgMinTPCCrossedRows", 70, "Minimum TPC crossed pad rows"}; - Configurable cfgMinCrossedRowsOverFindable{"cfgMinCrossedRowsOverFindable", 0.8f, "Minimum ratio crossedRows/findableClusters"}; - Configurable cfgMaxTPCSharedClusters{"cfgMaxTPCSharedClusters", 0, "Maximum TPC shared clusters"}; - Configurable cfgMaxTPCChi2NCl{"cfgMaxTPCChi2NCl", 4.0f, "Maximum TPC chi2/NCl"}; - Configurable cfgMaxITSChi2NCl{"cfgMaxITSChi2NCl", 36.0f, "Maximum ITS chi2/NCl"}; - Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require isPrimaryTrack flag"}; - Configurable requireGlobalTrackNoDCA{"requireGlobalTrackNoDCA", true, "Require isGlobalTrackWoDCA flag"}; - Configurable requirePVContributor{"requirePVContributor", true, "Require PV-contributor flag"}; - Configurable cfgLowPtTofNsigmaCut{"cfgLowPtTofNsigmaCut", 3.0f, "Low pt nSigma TOF cut for selecting pions and protons"}; - - Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "Maximum |DCAz| for all tracks [cm]"}; - Configurable> protonDCAPtBinEdges{"protonDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Proton pT bin edges for DCAxy cut [GeV/c]"}; - Configurable> protonMaxDCAxyPerPtBin{"protonMaxDCAxyPerPtBin", {0.10f, 0.08f, 0.05f, 0.05f}, "Max |DCAxy| for proton per pT bin [cm]"}; - Configurable> pionDCAPtBinEdges{"pionDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Pion pT bin edges for DCAxy cut [GeV/c]"}; - Configurable> pionMaxDCAxyPerPtBin{"pionMaxDCAxyPerPtBin", {0.20f, 0.15f, 0.10f, 0.08f}, "Max |DCAxy| for pion per pT bin [cm]"}; - - Configurable useTPCOnlyPID{"useTPCOnlyPID", false, "Use TPC-only PID (ignore TOF even if present)"}; - Configurable requireTOFForProton{"requireTOFForProton", false, "Require TOF signal for proton candidates"}; - Configurable requireTOFForPion{"requireTOFForPion", false, "Require TOF signal for pion candidates"}; - Configurable applyTOFCutWhenAvailableBelowThreshold{"applyTOFCutWhenAvailableBelowThreshold", false, "Apply TOF PID only when TOF information exists below momentum threshold"}; - - Configurable minProtonMomentum{"minProtonMomentum", 0.f, "Minimum proton momentum for TOF PID [GeV/c]"}; - Configurable minTPCNSigmaProton{"minTPCNSigmaProton", -6.0f, "Minimum (lower bound) TPC nSigma for proton"}; - Configurable minTOFNSigmaProton{"minTOFNSigmaProton", -6.0f, "Minimum (lower bound) TOF nSigma for proton"}; - Configurable minCombinedNSigmaProton{"minCombinedNSigmaProton", -6.0f, "Minimum combined nSigma for proton (asymmetric mode)"}; - Configurable maxTPCNSigmaProton{"maxTPCNSigmaProton", 3.0, "Maximum |TPC nSigma| for proton"}; - Configurable combinedNSigmaCutProton{"combinedNSigmaCutProton", 3.0, "Circular TPC+TOF combined nSigma cut for proton. Negative = asymmetric mode."}; - Configurable> protonTPCPIDMomentumBins{"protonTPCPIDMomentumBins", {0.f, 0.5f, 0.7f, 0.8f}, "Proton TPC PID momentum bin edges [GeV/c]"}; - Configurable> protonTPCNSigmaCutPerBin{"protonTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for proton per momentum bin"}; - Configurable> protonTOFPIDMomentumBins{"protonTOFPIDMomentumBins", {0.f, 999.f}, "Proton TOF PID momentum bin edges [GeV/c]"}; - Configurable> protonTOFNSigmaCutPerBin{"protonTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for proton per momentum bin"}; - - Configurable minPionMomentum{"minPionMomentum", 0.f, "Minimum pion momentum for TOF PID [GeV/c]"}; - Configurable minTPCNSigmaPion{"minTPCNSigmaPion", -6.0f, "Minimum (lower bound) TPC nSigma for pion"}; - Configurable minTOFNSigmaPion{"minTOFNSigmaPion", -6.0f, "Minimum (lower bound) TOF nSigma for pion"}; - Configurable minCombinedNSigmaPion{"minCombinedNSigmaPion", -6.0f, "Minimum combined nSigma for pion (asymmetric mode)"}; - Configurable maxTPCNSigmaPion{"maxTPCNSigmaPion", 3.0, "Maximum |TPC nSigma| for pion"}; - Configurable combinedNSigmaCutPion{"combinedNSigmaCutPion", 3.0, "Circular TPC+TOF combined nSigma cut for pion. Negative = asymmetric mode."}; - Configurable> pionTPCPIDMomentumBins{"pionTPCPIDMomentumBins", {0.f, 0.25f, 0.4f, 0.5f}, "Pion TPC PID momentum bin edges [GeV/c]"}; - Configurable> pionTPCNSigmaCutPerBin{"pionTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for pion per momentum bin"}; - Configurable> pionTOFPIDMomentumBins{"pionTOFPIDMomentumBins", {0.f, 999.f}, "Pion TOF PID momentum bin edges [GeV/c]"}; - Configurable> pionTOFNSigmaCutPerBin{"pionTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for pion per momentum bin"}; - - Configurable tpcNSigmaVetoThreshold{"tpcNSigmaVetoThreshold", 3.0f, "Reject track if TPC nSigma of a competing species is below this value"}; - Configurable tofNSigmaVetoThreshold{"tofNSigmaVetoThreshold", 3.0f, "Reject track if TOF nSigma of a competing species is below this value"}; - - Configurable applyDeepAngleCut{"applyDeepAngleCut", false, "Apply minimum opening-angle cut (removes split-track background)"}; - Configurable deepAngleCutValue{"deepAngleCutValue", 0.04, "Minimum opening angle between proton and pion [rad]"}; - - Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per signal event"}; - - Configurable enableRotationalBackground{"enableRotationalBackground", false, "Compute rotational background by rotating pion phi near PI"}; - Configurable numberOfRotations{"numberOfRotations", 10, "Number of pion-phi rotations for background"}; - Configurable rotationAngleWindow{"rotationAngleWindow", 6.f, "Pion rotated by angles within PI +/- PI/rotationAngleWindow"}; - - ConfigurableAxis cfgPtAxis{"cfgPtAxis", {VARIABLE_WIDTH, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4.0, 5.0, 7.0, 10.0}, "#it{p}_{T} (GeV/#it{c})"}; - ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f}, "Centrality (%)"}; - ConfigurableAxis cfgVtxAxis{"cfgVtxAxis", {VARIABLE_WIDTH, -12.f, -10.f, -9.f, -8.f, -7.f, -6.f, -5.f, -4.f, -3.f, -2.f, -1.f, 0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 12.f}, "Vertex z [cm]"}; - ConfigurableAxis cfgRapAxis{"cfgRapAxis", {20, -1.0, 1.0}, "Rapidity y"}; + // PDG service used only for truth-level charged-particle identification + // (see isChargedPrimaryMC() below). + Service pdgDB; + + struct : ConfigurableGroup { + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted |z-vertex| range [cm]"}; + Configurable cfgRequireRecoINELgt0{"cfgRequireRecoINELgt0", true, "Require reconstructed INEL>0 event class (at least one PV track in |eta|<1)"}; + Configurable cfgUseNoSameBunchPileupCut{"cfgUseNoSameBunchPileupCut", false, "Apply kNoSameBunchPileup event selection"}; + Configurable cfgUseGoodZvtxFT0vsPVCut{"cfgUseGoodZvtxFT0vsPVCut", false, "Apply kIsGoodZvtxFT0vsPV event selection"}; + Configurable applyOccupancyInTimeRangeCut{"applyOccupancyInTimeRangeCut", false, "Apply occupancy-in-time-range cut"}; + Configurable cfgOccupancyMin{"cfgOccupancyMin", 0, "Minimum track occupancy in time range"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 9999, "Maximum track occupancy in time range"}; + Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator: 0=FT0M 1=FT0A 2=FT0C 3=FV0A 4=NTPV"}; + Configurable cfgCentMin{"cfgCentMin", 0.f, "Minimum centrality percentile"}; + Configurable cfgCentMax{"cfgCentMax", 100.f, "Maximum centrality percentile"}; + } evSel; + + struct : ConfigurableGroup { + Configurable cfgCutPt{"cfgCutPt", 0.2f, "Minimum pT of daughter track [GeV/c]"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Maximum |eta| of daughter track"}; + Configurable cfgMinY{"cfgMinY", -0.5f, "Minimum rapidity of reconstructed Delta"}; + Configurable cfgMaxY{"cfgMaxY", 0.5f, "Maximum rapidity of reconstructed Delta"}; + Configurable cfgMinITSClusters{"cfgMinITSClusters", 5, "Minimum ITS clusters"}; + Configurable cfgMinTPCClusters{"cfgMinTPCClusters", 70, "Minimum TPC clusters found"}; + Configurable numberOfInvMassBins{"numberOfInvMassBins", 120, "Number of bins along the invariant mass axis"}; + Configurable cfgMinTPCCrossedRows{"cfgMinTPCCrossedRows", 70, "Minimum TPC crossed pad rows"}; + Configurable cfgMinCrossedRowsOverFindable{"cfgMinCrossedRowsOverFindable", 0.8f, "Minimum ratio crossedRows/findableClusters"}; + Configurable cfgMaxTPCSharedClusters{"cfgMaxTPCSharedClusters", 0, "Maximum TPC shared clusters"}; + Configurable cfgMaxTPCChi2NCl{"cfgMaxTPCChi2NCl", 4.0f, "Maximum TPC chi2/NCl"}; + Configurable cfgMaxITSChi2NCl{"cfgMaxITSChi2NCl", 36.0f, "Maximum ITS chi2/NCl"}; + Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require isPrimaryTrack flag"}; + Configurable requireGlobalTrackNoDCA{"requireGlobalTrackNoDCA", true, "Require isGlobalTrackWoDCA flag"}; + Configurable requirePVContributor{"requirePVContributor", true, "Require PV-contributor flag"}; + } trackCuts; + + struct : ConfigurableGroup { + Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "Maximum |DCAz| for all tracks [cm]"}; + Configurable> protonDCAPtBinEdges{"protonDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Proton pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> protonMaxDCAxyPerPtBin{"protonMaxDCAxyPerPtBin", {0.10f, 0.08f, 0.05f, 0.05f}, "Max |DCAxy| for proton per pT bin [cm]"}; + Configurable> pionDCAPtBinEdges{"pionDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 1000.f}, "Pion pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> pionMaxDCAxyPerPtBin{"pionMaxDCAxyPerPtBin", {0.20f, 0.15f, 0.10f, 0.08f}, "Max |DCAxy| for pion per pT bin [cm]"}; + } dcaCuts; + + struct : ConfigurableGroup { + Configurable useTPCOnlyPID{"useTPCOnlyPID", false, "Use TPC-only PID (ignore TOF even if present)"}; + Configurable applyTOFCutWhenAvailableBelowThreshold{"applyTOFCutWhenAvailableBelowThreshold", false, "Apply TOF PID only when TOF information exists below momentum threshold"}; + Configurable cfgLowPtTofNsigmaCut{"cfgLowPtTofNsigmaCut", 3.0f, "Low pt nSigma TOF cut for selecting pions and protons"}; + Configurable tpcNSigmaVetoThreshold{"tpcNSigmaVetoThreshold", 3.0f, "Reject track if TPC nSigma of a competing species is below this value"}; + Configurable tofNSigmaVetoThreshold{"tofNSigmaVetoThreshold", 3.0f, "Reject track if TOF nSigma of a competing species is below this value"}; + } pidShared; + + struct : ConfigurableGroup { + Configurable requireTOFForProton{"requireTOFForProton", false, "Require TOF signal for proton candidates"}; + Configurable minProtonMomentum{"minProtonMomentum", 0.f, "Minimum proton momentum for TOF PID [GeV/c]"}; + Configurable minTPCNSigmaProton{"minTPCNSigmaProton", -6.0f, "Minimum (lower bound) TPC nSigma for proton"}; + Configurable minTOFNSigmaProton{"minTOFNSigmaProton", -6.0f, "Minimum (lower bound) TOF nSigma for proton"}; + Configurable minCombinedNSigmaProton{"minCombinedNSigmaProton", -6.0f, "Minimum combined nSigma for proton (asymmetric mode)"}; + Configurable maxTPCNSigmaProton{"maxTPCNSigmaProton", 3.0, "Maximum |TPC nSigma| for proton"}; + Configurable combinedNSigmaCutProton{"combinedNSigmaCutProton", 3.0, "Circular TPC+TOF combined nSigma cut for proton. Negative = asymmetric mode."}; + Configurable> protonTPCPIDMomentumBins{"protonTPCPIDMomentumBins", {0.f, 0.5f, 0.7f, 0.8f}, "Proton TPC PID momentum bin edges [GeV/c]"}; + Configurable> protonTPCNSigmaCutPerBin{"protonTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for proton per momentum bin"}; + Configurable> protonTOFPIDMomentumBins{"protonTOFPIDMomentumBins", {0.f, 999.f}, "Proton TOF PID momentum bin edges [GeV/c]"}; + Configurable> protonTOFNSigmaCutPerBin{"protonTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for proton per momentum bin"}; + } protonPID; + + struct : ConfigurableGroup { + Configurable requireTOFForPion{"requireTOFForPion", false, "Require TOF signal for pion candidates"}; + Configurable minPionMomentum{"minPionMomentum", 0.f, "Minimum pion momentum for TOF PID [GeV/c]"}; + Configurable minTPCNSigmaPion{"minTPCNSigmaPion", -6.0f, "Minimum (lower bound) TPC nSigma for pion"}; + Configurable minTOFNSigmaPion{"minTOFNSigmaPion", -6.0f, "Minimum (lower bound) TOF nSigma for pion"}; + Configurable minCombinedNSigmaPion{"minCombinedNSigmaPion", -6.0f, "Minimum combined nSigma for pion (asymmetric mode)"}; + Configurable maxTPCNSigmaPion{"maxTPCNSigmaPion", 3.0, "Maximum |TPC nSigma| for pion"}; + Configurable combinedNSigmaCutPion{"combinedNSigmaCutPion", 3.0, "Circular TPC+TOF combined nSigma cut for pion. Negative = asymmetric mode."}; + Configurable> pionTPCPIDMomentumBins{"pionTPCPIDMomentumBins", {0.f, 0.25f, 0.4f, 0.5f}, "Pion TPC PID momentum bin edges [GeV/c]"}; + Configurable> pionTPCNSigmaCutPerBin{"pionTPCNSigmaCutPerBin", {5.f, 3.5f, 2.5f}, "Maximum TPC nSigma for pion per momentum bin"}; + Configurable> pionTOFPIDMomentumBins{"pionTOFPIDMomentumBins", {0.f, 999.f}, "Pion TOF PID momentum bin edges [GeV/c]"}; + Configurable> pionTOFNSigmaCutPerBin{"pionTOFNSigmaCutPerBin", {3.0f}, "Maximum TOF nSigma for pion per momentum bin"}; + } pionPID; + + struct : ConfigurableGroup { + Configurable applyDeepAngleCut{"applyDeepAngleCut", false, "Apply minimum opening-angle cut (removes split-track background)"}; + Configurable deepAngleCutValue{"deepAngleCutValue", 0.04, "Minimum opening angle between proton and pion [rad]"}; + } pairCuts; + + struct : ConfigurableGroup { + Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 5, "Number of mixed events per signal event"}; + } mixingCfg; + + struct : ConfigurableGroup { + Configurable enableRotationalBackground{"enableRotationalBackground", false, "Compute rotational background by rotating pion phi near PI"}; + Configurable numberOfRotations{"numberOfRotations", 10, "Number of pion-phi rotations for background"}; + Configurable rotationAngleWindow{"rotationAngleWindow", 6.f, "Pion rotated by angles within PI +/- PI/rotationAngleWindow"}; + } rotBkg; + + struct : ConfigurableGroup { + ConfigurableAxis cfgPtAxis{"cfgPtAxis", {VARIABLE_WIDTH, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.8, 2.0, 2.2, 2.4, 2.8, 3.2, 3.6, 4.0, 5.0, 7.0, 10.0}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f}, "Centrality (%)"}; + ConfigurableAxis cfgVtxAxis{"cfgVtxAxis", {VARIABLE_WIDTH, -12.f, -10.f, -9.f, -8.f, -7.f, -6.f, -5.f, -4.f, -3.f, -2.f, -1.f, 0.f, 1.f, 2.f, 3.f, 4.f, 5.f, 6.f, 7.f, 8.f, 9.f, 10.f, 12.f}, "Vertex z [cm]"}; + ConfigurableAxis cfgRapAxis{"cfgRapAxis", {20, -1.0, 1.0}, "Rapidity y"}; + } axes; std::vector mProtonTPCMomBins; std::vector mProtonTPCNSigCuts; @@ -181,31 +206,31 @@ struct DeltaAnalysis { << ", bins=" << binsVec.size() << ")"; } }; - chkSz(protonTPCNSigmaCutPerBin, protonTPCPIDMomentumBins, "proton TPC PID"); - chkSz(protonTOFNSigmaCutPerBin, protonTOFPIDMomentumBins, "proton TOF PID"); - chkSz(pionTPCNSigmaCutPerBin, pionTPCPIDMomentumBins, "pion TPC PID"); - chkSz(pionTOFNSigmaCutPerBin, pionTOFPIDMomentumBins, "pion TOF PID"); - chkSz(protonMaxDCAxyPerPtBin, protonDCAPtBinEdges, "proton DCA"); - chkSz(pionMaxDCAxyPerPtBin, pionDCAPtBinEdges, "pion DCA"); - - mProtonTPCMomBins = protonTPCPIDMomentumBins.value; - mProtonTPCNSigCuts = protonTPCNSigmaCutPerBin.value; - mProtonTOFMomBins = protonTOFPIDMomentumBins.value; - mProtonTOFNSigCuts = protonTOFNSigmaCutPerBin.value; - mPionTPCMomBins = pionTPCPIDMomentumBins.value; - mPionTPCNSigCuts = pionTPCNSigmaCutPerBin.value; - mPionTOFMomBins = pionTOFPIDMomentumBins.value; - mPionTOFNSigCuts = pionTOFNSigmaCutPerBin.value; - mProtonDCAPtEdges = protonDCAPtBinEdges.value; - mProtonMaxDCAxy = protonMaxDCAxyPerPtBin.value; - mPionDCAPtEdges = pionDCAPtBinEdges.value; - mPionMaxDCAxy = pionMaxDCAxyPerPtBin.value; + chkSz(protonPID.protonTPCNSigmaCutPerBin, protonPID.protonTPCPIDMomentumBins, "proton TPC PID"); + chkSz(protonPID.protonTOFNSigmaCutPerBin, protonPID.protonTOFPIDMomentumBins, "proton TOF PID"); + chkSz(pionPID.pionTPCNSigmaCutPerBin, pionPID.pionTPCPIDMomentumBins, "pion TPC PID"); + chkSz(pionPID.pionTOFNSigmaCutPerBin, pionPID.pionTOFPIDMomentumBins, "pion TOF PID"); + chkSz(dcaCuts.protonMaxDCAxyPerPtBin, dcaCuts.protonDCAPtBinEdges, "proton DCA"); + chkSz(dcaCuts.pionMaxDCAxyPerPtBin, dcaCuts.pionDCAPtBinEdges, "pion DCA"); + + mProtonTPCMomBins = protonPID.protonTPCPIDMomentumBins.value; + mProtonTPCNSigCuts = protonPID.protonTPCNSigmaCutPerBin.value; + mProtonTOFMomBins = protonPID.protonTOFPIDMomentumBins.value; + mProtonTOFNSigCuts = protonPID.protonTOFNSigmaCutPerBin.value; + mPionTPCMomBins = pionPID.pionTPCPIDMomentumBins.value; + mPionTPCNSigCuts = pionPID.pionTPCNSigmaCutPerBin.value; + mPionTOFMomBins = pionPID.pionTOFPIDMomentumBins.value; + mPionTOFNSigCuts = pionPID.pionTOFNSigmaCutPerBin.value; + mProtonDCAPtEdges = dcaCuts.protonDCAPtBinEdges.value; + mProtonMaxDCAxy = dcaCuts.protonMaxDCAxyPerPtBin.value; + mPionDCAPtEdges = dcaCuts.pionDCAPtBinEdges.value; + mPionMaxDCAxy = dcaCuts.pionMaxDCAxyPerPtBin.value; const AxisSpec ptAxis{200, 0., 10., "p_{T} (GeV/c)"}; - const AxisSpec massAxis{numberOfInvMassBins, 1.0, 8.0, "M_{inv} (GeV/#it{c}^{2})"}; - const AxisSpec centAxis{cfgCentAxis, "Centrality (%)"}; - const AxisSpec vtxAxis{cfgVtxAxis, "Vertex z [cm]"}; - const AxisSpec rapAxis{cfgRapAxis, "Rapidity y"}; + const AxisSpec massAxis{trackCuts.numberOfInvMassBins, 1.0, 8.0, "M_{inv} (GeV/#it{c}^{2})"}; + const AxisSpec centAxis{axes.cfgCentAxis, "Centrality (%)"}; + const AxisSpec vtxAxis{axes.cfgVtxAxis, "Vertex z [cm]"}; + const AxisSpec rapAxis{axes.cfgRapAxis, "Rapidity y"}; const AxisSpec nSigmaTPCaxis{100, -10., 10., "n#sigma^{TPC}"}; const AxisSpec nSigmaTOFaxis{100, -10., 10., "n#sigma^{TOF}"}; const AxisSpec momentumAxis{200, 0., 10., "p (GeV/#it{c})"}; @@ -281,7 +306,7 @@ struct DeltaAnalysis { histos.add("Analysis/hDeltaZeroInvMass", "#Delta^{0} (#rightarrow p + #pi^{-}) invariant mass", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hAntiDeltaZeroInvMass", "#bar{#Delta}^{0} (#rightarrow #bar{p} + #pi^{+}) invariant mass", kTH2F, {ptAxis, massAxis}); - if (enableRotationalBackground) { + if (rotBkg.enableRotationalBackground) { histos.add("Analysis/hDeltaPlusPlusInvMassRot", "#Delta^{++} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hAntiDeltaPlusPlusInvMassRot", "#bar{#Delta}^{++} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hDeltaZeroInvMassRot", "#Delta^{0} invariant mass - rotational background", kTH2F, {ptAxis, massAxis}); @@ -300,7 +325,7 @@ struct DeltaAnalysis { histos.add("THnSparse/hDeltaZero", "THnSparse #Delta^{0} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparse/hAntiDeltaZero", "THnSparse #bar{#Delta}^{0} same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - if (enableRotationalBackground) { + if (rotBkg.enableRotationalBackground) { histos.add("THnSparse/hDeltaPlusPlusRot", "THnSparse #Delta^{++} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparse/hAntiDeltaPlusPlusRot", "THnSparse #bar{#Delta}^{++} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparse/hDeltaZeroRot", "THnSparse #Delta^{0} rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); @@ -319,20 +344,11 @@ struct DeltaAnalysis { histos.add("THnSparse/hAntiDeltaPlusPlusMC", "THnSparse #bar{#Delta}^{++} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparse/hDeltaZeroMC", "THnSparse #Delta^{0} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparse/hAntiDeltaZeroMC", "THnSparse #bar{#Delta}^{0} MC reconstructed", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - // CHANGED: generated-particle THnSparse now carries mass at axis 0, followed by pt, cent, rap - histos.add("THnSparse/hDeltaPlusPlusGen", "THnSparse #Delta^{++} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - histos.add("THnSparse/hAntiDeltaPlusPlusGen", "THnSparse #bar{#Delta}^{++} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - histos.add("THnSparse/hDeltaZeroGen", "THnSparse #Delta^{0} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - histos.add("THnSparse/hAntiDeltaZeroGen", "THnSparse #bar{#Delta}^{0} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("Analysis/hDeltaPlusPlusInvMassMC", "#Delta^{++} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hAntiDeltaPlusPlusInvMassMC", "#bar{#Delta}^{++} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hDeltaZeroInvMassMC", "#Delta^{0} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); histos.add("Analysis/hAntiDeltaZeroInvMassMC", "#bar{#Delta}^{0} invariant mass (MC truth-matched)", kTH2F, {ptAxis, massAxis}); - histos.add("Analysis/hDeltaPlusPlusInvMassGen", "#Delta^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); - histos.add("Analysis/hAntiDeltaPlusPlusInvMassGen", "#bar{#Delta}^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); - histos.add("Analysis/hDeltaZeroInvMassGen", "#Delta^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); - histos.add("Analysis/hAntiDeltaZeroInvMassGen", "#bar{#Delta}^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); histos.add("QAChecks/hRecProtonDeltaPlusPlus", "Rec proton from #Delta^{++}", kTH1F, {ptAxis}); histos.add("QAChecks/hRecProtonAntiDeltaPlusPlus", "Rec proton from #bar{#Delta}^{++}", kTH1F, {ptAxis}); @@ -342,16 +358,7 @@ struct DeltaAnalysis { histos.add("QAChecks/hRecPionAntiDeltaPlusPlus", "Rec pion from #bar{#Delta}^{++}", kTH1F, {ptAxis}); histos.add("QAChecks/hRecPionDeltaZero", "Rec pion from #Delta^{0}", kTH1F, {ptAxis}); histos.add("QAChecks/hRecPionAntiDeltaZero", "Rec pion from #bar{#Delta}^{0}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenProtonDeltaPlusPlus", "Gen proton from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenProtonAntiDeltaPlusPlus", "Gen proton from #bar{#Delta}^{++}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenProtonDeltaZero", "Gen proton from #Delta^{0}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenProtonAntiDeltaZero", "Gen proton from #bar{#Delta}^{0}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenPionDeltaPlusPlus", "Gen pion from #Delta^{++}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenPionAntiDeltaPlusPlus", "Gen pion from #bar{#Delta}^{++}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenPionDeltaZero", "Gen pion from #Delta^{0}", kTH1F, {ptAxis}); - histos.add("QAChecks/hGenPionAntiDeltaZero", "Gen pion from #bar{#Delta}^{0}", kTH1F, {ptAxis}); - - // reconstructed-level (pre-truth-matching) invariant mass histograms for MC-reconstructed candidates + histos.add("AnalysisMCReco/hDeltaPlusPlusInvMassReco", "#Delta^{++} invariant mass - MC reconstructed (pre-truth-matching)", kTH2F, {ptAxis, massAxis}); histos.add("AnalysisMCReco/hAntiDeltaPlusPlusInvMassReco", "#bar{#Delta}^{++} invariant mass - MC reconstructed (pre-truth-matching)", kTH2F, {ptAxis, massAxis}); histos.add("AnalysisMCReco/hDeltaZeroInvMassReco", "#Delta^{0} invariant mass - MC reconstructed (pre-truth-matching)", kTH2F, {ptAxis, massAxis}); @@ -362,7 +369,6 @@ struct DeltaAnalysis { histos.add("THnSparseMCReco/hDeltaZeroReco", "THnSparse #Delta^{0} MC reconstructed (pre-truth-matching)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); histos.add("THnSparseMCReco/hAntiDeltaZeroReco", "THnSparse #bar{#Delta}^{0} MC reconstructed (pre-truth-matching)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - // detector QA for reconstructed MC tracks,no truth matching required ── histos.add("QAMC/Proton/dcaXYvsPt", "Proton DCA_{xy} vs p_{T} (MC reco, after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); histos.add("QAMC/Proton/dcaZvsPt", "Proton DCA_{z} vs p_{T} (MC reco, after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); histos.add("QAMC/Proton/tpcNSigmaVsMomentum", "Proton TPC n#sigma vs p (MC reco, after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); @@ -396,13 +402,111 @@ struct DeltaAnalysis { histos.add("QAMC/Pion/tpcClustersFoundVsPt", "Pion TPC clusters found vs p_{T} (MC reco)", kTH2F, {ptForPIDAxis, tpcClusAxis}); histos.add("QAMC/Pion/dcaXYdist", "Pion DCA_{xy} distribution (MC reco, fine bins)", kTH1F, {dcaXYaxis}); histos.add("QAMC/Pion/dcaZdist", "Pion DCA_{z} distribution (MC reco, fine bins)", kTH1F, {dcaZaxis}); + + histos.add("MCRecoEvent/hRecoEvents", "Reconstructed INEL>0 events (Nrec, MC reco)", kTH1F, {centAxis}); + } + + // ── processMCGen(): generated-level Delta spectra + QA ───────────────────────────────── + if (doprocessMCGen) { + histos.add("MCGen/hDeltaPlusPlusGenSparse", "THnSparse #Delta^{++} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/hAntiDeltaPlusPlusGenSparse", "THnSparse #bar{#Delta}^{++} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/hDeltaZeroGenSparse", "THnSparse #Delta^{0} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/hAntiDeltaZeroGenSparse", "THnSparse #bar{#Delta}^{0} generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + histos.add("MCGen/hDeltaPlusPlusInvMassGen", "#Delta^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/hAntiDeltaPlusPlusInvMassGen", "#bar{#Delta}^{++} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/hDeltaZeroInvMassGen", "#Delta^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/hAntiDeltaZeroInvMassGen", "#bar{#Delta}^{0} invariant mass - generated", kTH2F, {ptAxis, massAxis}); + + histos.add("MCGenQA/hGenProtonDeltaPlusPlus", "Gen proton from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenProtonAntiDeltaPlusPlus", "Gen proton from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenProtonDeltaZero", "Gen proton from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenProtonAntiDeltaZero", "Gen proton from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenPionDeltaPlusPlus", "Gen pion from #Delta^{++}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenPionAntiDeltaPlusPlus", "Gen pion from #bar{#Delta}^{++}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenPionDeltaZero", "Gen pion from #Delta^{0}", kTH1F, {ptAxis}); + histos.add("MCGenQA/hGenPionAntiDeltaZero", "Gen pion from #bar{#Delta}^{0}", kTH1F, {ptAxis}); + + // Generated-Delta cut flow. 0=All generated, 1=producedByGenerator, 2=truth vertex-z, + // 3=truth INEL>0, 4=PDG, 5=rapidity, 6=decay channel, 7=final generated histogram. + histos.add("MCGenQA/hGenDeltaCutFlow", "Generated #Delta cut flow", kTH1F, {{8, -0.5f, 7.5f}}); + + histos.add("MCGen/GeneratedDelta_EventAccepted/hDeltaPlusPlusInvMass", "#Delta^{++} generated, event-accepted", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaPlusPlusInvMass", "#bar{#Delta}^{++} generated, event-accepted", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hDeltaZeroInvMass", "#Delta^{0} generated, event-accepted", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaZeroInvMass", "#bar{#Delta}^{0} generated, event-accepted", kTH2F, {ptAxis, massAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hDeltaPlusPlusSparse", "THnSparse #Delta^{++} generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaPlusPlusSparse", "THnSparse #bar{#Delta}^{++} generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hDeltaZeroSparse", "THnSparse #Delta^{0} generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaZeroSparse", "THnSparse #bar{#Delta}^{0} generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + } + + // ── processEventFactor(): event-level corrections only (Event Loss / Event Splitting) ── + if (doprocessEventFactor) { + histos.add("EventFactor/hRecoEvents", "Reconstructed INEL>0 events (Nrec)", kTH1F, {centAxis}); + histos.add("EventFactor/hEventsGenAll", "Generated events passing truth |Zvtx| + truth INEL>0 cuts", kTH1F, {centAxis}); + histos.add("EventFactor/hEventsGenAccepted", "Generated events with >=1 accepted reconstructed collision", kTH1F, {centAxis}); + histos.add("EventFactor/hNRecoCollisionsPerMcCollision", "Number of reconstructed collisions per generated collision", kTH1F, {{21, -0.5f, 20.5f}}); + } + + // ── Event-level cut-flow (QA-only) histograms ─────────────────────────────────────────── + // These are pure bookkeeping: they record how many events survive each existing selection + // step (already implemented in passesEventSelectionImpl(), isTruthInelGt0(), and + // processEventFactor()) and do not themselves apply, add, or alter any selection. + { + histos.add("CutFlow/Data/hEventCutFlow", "Data event cut flow", kTH1F, {{9, -0.5f, 8.5f}}); + auto hDataCutFlow = histos.get(HIST("CutFlow/Data/hEventCutFlow")); + hDataCutFlow->GetXaxis()->SetBinLabel(1, "All"); + hDataCutFlow->GetXaxis()->SetBinLabel(2, "sel8"); + hDataCutFlow->GetXaxis()->SetBinLabel(3, "INEL>0"); + hDataCutFlow->GetXaxis()->SetBinLabel(4, "Vz"); + hDataCutFlow->GetXaxis()->SetBinLabel(5, "Occupancy"); + hDataCutFlow->GetXaxis()->SetBinLabel(6, "Centrality"); + hDataCutFlow->GetXaxis()->SetBinLabel(7, "Pileup"); + hDataCutFlow->GetXaxis()->SetBinLabel(8, "GoodZvtx"); + hDataCutFlow->GetXaxis()->SetBinLabel(9, "Accepted"); + } + + if (doprocessMC) { + histos.add("CutFlow/MCReco/hEventCutFlow", "MC reconstructed event cut flow", kTH1F, {{9, -0.5f, 8.5f}}); + auto hMcRecoCutFlow = histos.get(HIST("CutFlow/MCReco/hEventCutFlow")); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(1, "All"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(2, "sel8"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(3, "INEL>0"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(4, "Vz"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(5, "Occupancy"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(6, "Centrality"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(7, "Pileup"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(8, "GoodZvtx"); + hMcRecoCutFlow->GetXaxis()->SetBinLabel(9, "Accepted"); + } + + if (doprocessMCGen) { + histos.add("CutFlow/MCGen/hEventCutFlow", "MC generated event cut flow", kTH1F, {{5, -0.5f, 4.5f}}); + auto hMcGenCutFlow = histos.get(HIST("CutFlow/MCGen/hEventCutFlow")); + hMcGenCutFlow->GetXaxis()->SetBinLabel(1, "All"); + hMcGenCutFlow->GetXaxis()->SetBinLabel(2, "Vz"); + hMcGenCutFlow->GetXaxis()->SetBinLabel(3, "INEL>0"); + hMcGenCutFlow->GetXaxis()->SetBinLabel(4, "Centrality"); + hMcGenCutFlow->GetXaxis()->SetBinLabel(5, "Accepted"); + } + + if (doprocessEventFactor) { + histos.add("CutFlow/EventFactor/hEventAcceptedCutFlow", "Event Loss / Event Splitting: generated-event cut flow", kTH1F, {{6, -0.5f, 5.5f}}); + auto hEventFactorCutFlow = histos.get(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow")); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(1, "All"); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(2, "Vz"); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(3, "INEL>0"); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(4, "HasRecoColl"); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(5, "RecoCollAccepted"); + hEventFactorCutFlow->GetXaxis()->SetBinLabel(6, "EventAccepted"); } } // end init() template float getCentrality(CollisionType const& collision) { - switch (cfgCentralityEstimator) { + switch (evSel.cfgCentralityEstimator) { case delta_analysis::kFT0M: return collision.centFT0M(); case delta_analysis::kFT0A: @@ -418,70 +522,147 @@ struct DeltaAnalysis { } } - // passesEventSelection with optional EventSelection cuts + // INEL>0 (reconstructed): at least one PV-associated track within |eta| < 1. template - bool passesEventSelection(CollisionType const& collision) + bool isRecoInelGt0(CollisionType const& collision) + { + return collision.multNTracksPVeta1() > 0; + } + + // Truth-level charged-primary identification, used only to build the truth-level INEL>0 event + // class (isTruthInelGt0 below). + template + bool isChargedPrimaryMC(McParticleType const& mcPart) + { + if (!mcPart.isPhysicalPrimary()) { + return false; + } + auto* pdgParticle = pdgDB->GetParticle(mcPart.pdgCode()); + if (!pdgParticle) { + return false; + } + return std::abs(pdgParticle->Charge()) > minAbsCharge; + } + + // Truth-level analogue of isRecoInelGt0() above - at least one charged primary within + // |eta| < 1 among the generated particles of a given MC collision. + template + bool isTruthInelGt0(McParticleTableType const& mcParticlesInCollision) { + for (auto const& mcPart : mcParticlesInCollision) { + if (std::abs(mcPart.eta()) < 1.f && isChargedPrimaryMC(mcPart)) { + return true; + } + } + return false; + } + + struct NoCutFlowTag { + }; + struct DataCutFlowTag { + }; + struct McRecoCutFlowTag { + }; + + template + void fillEventCutFlowBin(float bin) + { + if constexpr (std::is_same_v) { + histos.fill(HIST("CutFlow/Data/hEventCutFlow"), bin); + } else if constexpr (std::is_same_v) { + histos.fill(HIST("CutFlow/MCReco/hEventCutFlow"), bin); + } + // NoCutFlowTag: intentionally does nothing. + } + + template + bool passesEventSelectionImpl(CollisionType const& collision) + { + fillEventCutFlowBin(0.f); // All collisions + if (!collision.sel8()) { return false; } - if (std::abs(collision.posZ()) > cfgCutVertex) { + fillEventCutFlowBin(1.f); // sel8 + + if (evSel.cfgRequireRecoINELgt0 && !isRecoInelGt0(collision)) { return false; } - if (applyOccupancyInTimeRangeCut) { + fillEventCutFlowBin(2.f); // Reconstructed INEL>0 + + if (std::abs(collision.posZ()) > evSel.cfgCutVertex) { + return false; + } + fillEventCutFlowBin(3.f); // |Vz| < cfgCutVertex + + if (evSel.applyOccupancyInTimeRangeCut) { const int occ = collision.trackOccupancyInTimeRange(); - if (occ < cfgOccupancyMin || occ > cfgOccupancyMax) { + if (occ < evSel.cfgOccupancyMin || occ > evSel.cfgOccupancyMax) { return false; } } + fillEventCutFlowBin(4.f); // Occupancy cut + const float cent = getCentrality(collision); - if (cent < cfgCentMin || cent > cfgCentMax) { + if (cent < evSel.cfgCentMin || cent > evSel.cfgCentMax) { return false; } + fillEventCutFlowBin(5.f); // Centrality range - if (cfgUseNoSameBunchPileupCut && + if (evSel.cfgUseNoSameBunchPileupCut && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { return false; } + fillEventCutFlowBin(6.f); // NoSameBunchPileup - if (cfgUseGoodZvtxFT0vsPVCut && + if (evSel.cfgUseGoodZvtxFT0vsPVCut && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; } + fillEventCutFlowBin(7.f); // GoodZvtxFT0vsPV + + fillEventCutFlowBin(8.f); // Final accepted event return true; } + // passesEventSelection with optional EventSelection cuts + template + bool passesEventSelection(CollisionType const& collision) + { + return passesEventSelectionImpl(collision); + } + template bool passesBasicTrackSelection(TrackType const& track) { - if (track.itsNCls() < cfgMinITSClusters) { + if (track.itsNCls() < trackCuts.cfgMinITSClusters) { return false; } - if (track.tpcNClsShared() > cfgMaxTPCSharedClusters) { + if (track.tpcNClsShared() > trackCuts.cfgMaxTPCSharedClusters) { return false; } - if (track.tpcNClsFound() < cfgMinTPCClusters) { + if (track.tpcNClsFound() < trackCuts.cfgMinTPCClusters) { return false; } - if (track.tpcNClsCrossedRows() < cfgMinTPCCrossedRows) { + if (track.tpcNClsCrossedRows() < trackCuts.cfgMinTPCCrossedRows) { return false; } - if (track.tpcNClsCrossedRows() < cfgMinCrossedRowsOverFindable * track.tpcNClsFindable()) { + if (track.tpcNClsCrossedRows() < trackCuts.cfgMinCrossedRowsOverFindable * track.tpcNClsFindable()) { return false; } - if (track.tpcChi2NCl() > cfgMaxTPCChi2NCl) { + if (track.tpcChi2NCl() > trackCuts.cfgMaxTPCChi2NCl) { return false; } - if (track.itsChi2NCl() > cfgMaxITSChi2NCl) { + if (track.itsChi2NCl() > trackCuts.cfgMaxITSChi2NCl) { return false; } - if (requirePrimaryTrack && !track.isPrimaryTrack()) { + if (trackCuts.requirePrimaryTrack && !track.isPrimaryTrack()) { return false; } - if (requireGlobalTrackNoDCA && !track.isGlobalTrackWoDCA()) { + if (trackCuts.requireGlobalTrackNoDCA && !track.isGlobalTrackWoDCA()) { return false; } - if (requirePVContributor && !track.isPVContributor()) { + if (trackCuts.requirePVContributor && !track.isPVContributor()) { return false; } return true; @@ -500,7 +681,7 @@ struct DeltaAnalysis { break; } } - return passed && (std::abs(track.dcaZ()) < cfgCutDCAz); + return passed && (std::abs(track.dcaZ()) < dcaCuts.cfgCutDCAz); } template @@ -516,7 +697,7 @@ struct DeltaAnalysis { break; } } - return passed && (std::abs(track.dcaZ()) < cfgCutDCAz); + return passed && (std::abs(track.dcaZ()) < dcaCuts.cfgCutDCAz); } template @@ -531,58 +712,54 @@ struct DeltaAnalysis { const float tofNSigPr = std::abs(track.tofNSigmaPr()); const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr; const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi; - const auto circularCutSq = static_cast(combinedNSigmaCutProton * combinedNSigmaCutProton); + const auto circularCutSq = static_cast(protonPID.combinedNSigmaCutProton * protonPID.combinedNSigmaCutProton); - const float circularVetoCutSq = tpcNSigmaVetoThreshold * tpcNSigmaVetoThreshold + tofNSigmaVetoThreshold * tofNSigmaVetoThreshold; + const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold + pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold; - if (!useTPCOnlyPID && track.hasTOF()) { - if (combinedNSigmaCutProton < 0 && totalMomentum >= minProtonMomentum) { - if (track.tofNSigmaPr() < minTOFNSigmaProton) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { + if (protonPID.combinedNSigmaCutProton < 0 && totalMomentum >= protonPID.minProtonMomentum) { + if (track.tofNSigmaPr() < protonPID.minTOFNSigmaProton) { return false; } for (int i = 0; i < nTOFBins - 1; ++i) { if (totalMomentum >= mProtonTOFMomBins[i] && totalMomentum < mProtonTOFMomBins[i + 1] && - tofNSigPr < mProtonTOFNSigCuts[i] && tofNSigPi > tofNSigmaVetoThreshold) { + tofNSigPr < mProtonTOFNSigCuts[i] && tofNSigPi > pidShared.tofNSigmaVetoThreshold) { tofPassed = true; break; } } - if (track.tpcNSigmaPr() < minCombinedNSigmaProton) { + if (track.tpcNSigmaPr() < protonPID.minCombinedNSigmaProton) { return false; } - if (tpcNSigPr < static_cast(maxTPCNSigmaProton) && tpcNSigPi > tpcNSigmaVetoThreshold) { + if (tpcNSigPr < static_cast(protonPID.maxTPCNSigmaProton) && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold) { tpcPassed = true; } - } else if (combinedNSigmaCutProton > 0 && totalMomentum >= minProtonMomentum) { + } else if (protonPID.combinedNSigmaCutProton > 0 && totalMomentum >= protonPID.minProtonMomentum) { if (combinedNSigPr < circularCutSq && combinedNSigPi > circularVetoCutSq) { tpcPassed = true; tofPassed = true; } } - // ── Part 2: low-momentum region with TOF available ────────────────── - if (totalMomentum < minProtonMomentum && tpcNSigPr < static_cast(maxTPCNSigmaProton)) { + if (totalMomentum < protonPID.minProtonMomentum && tpcNSigPr < static_cast(protonPID.maxTPCNSigmaProton)) { tpcPassed = true; - if (applyTOFCutWhenAvailableBelowThreshold) { - // TOF is available (track.hasTOF() is true in this branch): require TOF nSigma - tofPassed = (tofNSigPr < cfgLowPtTofNsigmaCut); + if (pidShared.applyTOFCutWhenAvailableBelowThreshold) { + tofPassed = (tofNSigPr < pidShared.cfgLowPtTofNsigmaCut); } else { tofPassed = true; } } } else { - // No TOF available (or useTPCOnlyPID == true) - if (totalMomentum < minProtonMomentum && tpcNSigPr < static_cast(maxTPCNSigmaProton)) { - // ── Part 2: low-momentum, no TOF → TPC-only regardless of configurable ── + if (totalMomentum < protonPID.minProtonMomentum && tpcNSigPr < static_cast(protonPID.maxTPCNSigmaProton)) { tpcPassed = true; tofPassed = true; } else { tofPassed = true; - if (track.tpcNSigmaPr() < minTPCNSigmaProton) { + if (track.tpcNSigmaPr() < protonPID.minTPCNSigmaProton) { return false; } for (int i = 0; i < nTPCBins - 1; ++i) { if (totalMomentum >= mProtonTPCMomBins[i] && totalMomentum < mProtonTPCMomBins[i + 1] && - tpcNSigPr < mProtonTPCNSigCuts[i] && tpcNSigPi > tpcNSigmaVetoThreshold) { + tpcNSigPr < mProtonTPCNSigCuts[i] && tpcNSigPi > pidShared.tpcNSigmaVetoThreshold) { tpcPassed = true; break; } @@ -604,57 +781,53 @@ struct DeltaAnalysis { const float tofNSigPr = std::abs(track.tofNSigmaPr()); const float combinedNSigPi = tpcNSigPi * tpcNSigPi + tofNSigPi * tofNSigPi; const float combinedNSigPr = tpcNSigPr * tpcNSigPr + tofNSigPr * tofNSigPr; - const auto circularCutSq = static_cast(combinedNSigmaCutPion * combinedNSigmaCutPion); - const float circularVetoCutSq = tpcNSigmaVetoThreshold * tpcNSigmaVetoThreshold + tofNSigmaVetoThreshold * tofNSigmaVetoThreshold; + const auto circularCutSq = static_cast(pionPID.combinedNSigmaCutPion * pionPID.combinedNSigmaCutPion); + const float circularVetoCutSq = pidShared.tpcNSigmaVetoThreshold * pidShared.tpcNSigmaVetoThreshold + pidShared.tofNSigmaVetoThreshold * pidShared.tofNSigmaVetoThreshold; - if (!useTPCOnlyPID && track.hasTOF()) { - if (combinedNSigmaCutPion < 0 && totalMomentum >= minPionMomentum) { - if (track.tofNSigmaPi() < minTOFNSigmaPion) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { + if (pionPID.combinedNSigmaCutPion < 0 && totalMomentum >= pionPID.minPionMomentum) { + if (track.tofNSigmaPi() < pionPID.minTOFNSigmaPion) { return false; } for (int i = 0; i < nTOFBins - 1; ++i) { if (totalMomentum >= mPionTOFMomBins[i] && totalMomentum < mPionTOFMomBins[i + 1] && - tofNSigPi < mPionTOFNSigCuts[i] && tofNSigPr > tofNSigmaVetoThreshold) { + tofNSigPi < mPionTOFNSigCuts[i] && tofNSigPr > pidShared.tofNSigmaVetoThreshold) { tofPassed = true; break; } } - if (track.tpcNSigmaPi() < minCombinedNSigmaPion) { + if (track.tpcNSigmaPi() < pionPID.minCombinedNSigmaPion) { return false; } - if (tpcNSigPi < static_cast(maxTPCNSigmaPion) && tpcNSigPr > tpcNSigmaVetoThreshold) { + if (tpcNSigPi < static_cast(pionPID.maxTPCNSigmaPion) && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold) { tpcPassed = true; } - } else if (combinedNSigmaCutPion > 0 && totalMomentum >= minPionMomentum) { + } else if (pionPID.combinedNSigmaCutPion > 0 && totalMomentum >= pionPID.minPionMomentum) { if (combinedNSigPi < circularCutSq && combinedNSigPr > circularVetoCutSq) { tpcPassed = true; tofPassed = true; } } - // ── Part 2: low-momentum region with TOF available ────────────────── - if (totalMomentum < minPionMomentum && tpcNSigPi < static_cast(maxTPCNSigmaPion)) { + if (totalMomentum < pionPID.minPionMomentum && tpcNSigPi < static_cast(pionPID.maxTPCNSigmaPion)) { tpcPassed = true; - if (applyTOFCutWhenAvailableBelowThreshold) { - // TOF is available (track.hasTOF() is true in this branch): require TOF nSigma - tofPassed = (tofNSigPi < cfgLowPtTofNsigmaCut); + if (pidShared.applyTOFCutWhenAvailableBelowThreshold) { + tofPassed = (tofNSigPi < pidShared.cfgLowPtTofNsigmaCut); } else { tofPassed = true; } } } else { - // No TOF available (or useTPCOnlyPID == true) - if (totalMomentum < minPionMomentum && tpcNSigPi < static_cast(maxTPCNSigmaPion)) { - // ── Part 2: low-momentum, no TOF → TPC-only regardless of configurable ── + if (totalMomentum < pionPID.minPionMomentum && tpcNSigPi < static_cast(pionPID.maxTPCNSigmaPion)) { tpcPassed = true; tofPassed = true; } else { tofPassed = true; - if (track.tpcNSigmaPi() < minTPCNSigmaPion) { + if (track.tpcNSigmaPi() < pionPID.minTPCNSigmaPion) { return false; } for (int i = 0; i < nTPCBins - 1; ++i) { if (totalMomentum >= mPionTPCMomBins[i] && totalMomentum < mPionTPCMomBins[i + 1] && - tpcNSigPi < mPionTPCNSigCuts[i] && tpcNSigPr > tpcNSigmaVetoThreshold) { + tpcNSigPi < mPionTPCNSigCuts[i] && tpcNSigPr > pidShared.tpcNSigmaVetoThreshold) { tpcPassed = true; break; } @@ -680,7 +853,7 @@ struct DeltaAnalysis { histos.fill(HIST("QAafter/Proton/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); histos.fill(HIST("QAafter/Proton/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); histos.fill(HIST("QAafter/Proton/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); - if (!useTPCOnlyPID && track.hasTOF()) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { const float tofNSigPr = track.tofNSigmaPr(); histos.fill(HIST("QAafter/Proton/tofNSigmaVsMomentum"), totalMomentum, tofNSigPr); histos.fill(HIST("QAafter/Proton/tofNSigmaVsPt"), pt, tofNSigPr); @@ -706,7 +879,7 @@ struct DeltaAnalysis { histos.fill(HIST("QAafter/Pion/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); histos.fill(HIST("QAafter/Pion/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); histos.fill(HIST("QAafter/Pion/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); - if (!useTPCOnlyPID && track.hasTOF()) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { const float tofNSigPi = track.tofNSigmaPi(); histos.fill(HIST("QAafter/Pion/tofNSigmaVsMomentum"), totalMomentum, tofNSigPi); histos.fill(HIST("QAafter/Pion/tofNSigmaVsPt"), pt, tofNSigPi); @@ -716,7 +889,6 @@ struct DeltaAnalysis { } } - // detector QA for reconstructed MC tracks, no truth matching required ── template void fillQAMCProton(TrackType const& track, float totalMomentum, float centralityPercent) { @@ -733,7 +905,7 @@ struct DeltaAnalysis { histos.fill(HIST("QAMC/Proton/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); histos.fill(HIST("QAMC/Proton/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); histos.fill(HIST("QAMC/Proton/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); - if (!useTPCOnlyPID && track.hasTOF()) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { const float tofNSigPr = track.tofNSigmaPr(); histos.fill(HIST("QAMC/Proton/tofNSigmaVsMomentum"), totalMomentum, tofNSigPr); histos.fill(HIST("QAMC/Proton/tofNSigmaVsPt"), pt, tofNSigPr); @@ -759,7 +931,7 @@ struct DeltaAnalysis { histos.fill(HIST("QAMC/Pion/tpcNSigmaKaonContamVsPt"), pt, track.tpcNSigmaKa()); histos.fill(HIST("QAMC/Pion/tpcCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); histos.fill(HIST("QAMC/Pion/tpcClustersFoundVsPt"), pt, track.tpcNClsFound()); - if (!useTPCOnlyPID && track.hasTOF()) { + if (!pidShared.useTPCOnlyPID && track.hasTOF()) { const float tofNSigPi = track.tofNSigmaPi(); histos.fill(HIST("QAMC/Pion/tofNSigmaVsMomentum"), totalMomentum, tofNSigPi); histos.fill(HIST("QAMC/Pion/tofNSigmaVsPt"), pt, tofNSigPi); @@ -832,7 +1004,6 @@ struct DeltaAnalysis { } } - // reconstructed-level (pre-truth-matching) invariant mass ── void fillDeltaHistogramMCReco(int protonSign, int pionSign, float pairPt, float pairMass, float centrality, float rapidity) { if (protonSign > 0) { @@ -856,19 +1027,19 @@ struct DeltaAnalysis { void fillRotationalBackground(int protonSign, int pionSign, float pxProton, float pyProton, float pzProton, float pionPhi, float pionPt, float pzPion, float centrality) { - if (numberOfRotations <= 0) { + if (rotBkg.numberOfRotations <= 0) { return; } - const float weight = 1.f / static_cast(numberOfRotations); - const float rotWindowHalf = o2::constants::math::PI / static_cast(rotationAngleWindow); + const float weight = 1.f / static_cast(rotBkg.numberOfRotations); + const float rotWindowHalf = o2::constants::math::PI / static_cast(rotBkg.rotationAngleWindow); - for (int iRot = 0; iRot < numberOfRotations; ++iRot) { + for (int iRot = 0; iRot < rotBkg.numberOfRotations; ++iRot) { float rotAngle = 0.f; - if (numberOfRotations == 1) { + if (rotBkg.numberOfRotations == 1) { rotAngle = o2::constants::math::PI; } else { rotAngle = (o2::constants::math::PI - rotWindowHalf) + - static_cast(iRot) * (2.f * rotWindowHalf / static_cast(numberOfRotations - 1)); + static_cast(iRot) * (2.f * rotWindowHalf / static_cast(rotBkg.numberOfRotations - 1)); } const float newPhi = RecoDecay::constrainAngle(pionPhi + rotAngle, 0.f); const float pxPionRot = pionPt * std::cos(newPhi); @@ -880,7 +1051,7 @@ struct DeltaAnalysis { const float rotPt = RecoDecay::pt(std::array{pxProton + pxPionRot, pyProton + pyPionRot}); const float rotY = RecoDecay::y( std::array{pxProton + pxPionRot, pyProton + pyPionRot, pzProton + pzPion}, rotMass); - if (rotY < cfgMinY || rotY > cfgMaxY) { + if (rotY < trackCuts.cfgMinY || rotY > trackCuts.cfgMaxY) { continue; } @@ -936,9 +1107,9 @@ struct DeltaAnalysis { } const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); const bool isProton = passesProtonPID(track, mom) && passesProtonDCASelection(track) && - !(requireTOFForProton && !track.hasTOF()); + !(protonPID.requireTOFForProton && !track.hasTOF()); const bool isPion = passesPionPID(track, mom) && passesPionDCASelection(track) && - !(requireTOFForPion && !track.hasTOF()); + !(pionPID.requireTOFForPion && !track.hasTOF()); if (!isProton && !isPion) { continue; } @@ -991,10 +1162,10 @@ struct DeltaAnalysis { } } - if (applyDeepAngleCut) { + if (pairCuts.applyDeepAngleCut) { const float cosAngle = std::clamp( (pxPr * pxPi + pyPr * pyPi + pzPr * pzPi) / (protonCand.mom * pionCand.mom), -1.f, 1.f); - if (std::acos(cosAngle) < static_cast(deepAngleCutValue)) { + if (std::acos(cosAngle) < static_cast(pairCuts.deepAngleCutValue)) { continue; } } @@ -1006,7 +1177,7 @@ struct DeltaAnalysis { const float pairPt = RecoDecay::pt(std::array{pxPr + pxPi, pyPr + pyPi}); const float pairY = RecoDecay::y(std::array{pxPr + pxPi, pyPr + pyPi, pzPr + pzPi}, pairMass); - if (pairY < cfgMinY || pairY > cfgMaxY) { + if (pairY < trackCuts.cfgMinY || pairY > trackCuts.cfgMaxY) { continue; } @@ -1020,7 +1191,7 @@ struct DeltaAnalysis { fillDeltaHistogramMixedEvent(protonCand.sign, pionCand.sign, pairPt, pairMass, centrality, pairY); } else { fillDeltaHistogramSameEvent(protonCand.sign, pionCand.sign, pairPt, pairMass, centrality, pairY); - if (enableRotationalBackground) { + if (rotBkg.enableRotationalBackground) { fillRotationalBackground(protonCand.sign, pionCand.sign, pxPr, pyPr, pzPr, pionCand.phi, pionCand.pt, pzPi, centrality); } @@ -1029,15 +1200,18 @@ struct DeltaAnalysis { } } - Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex; - Filter acceptanceFilter = (nabs(aod::track::eta) < cfgCutEta && nabs(aod::track::pt) > cfgCutPt); + Filter collisionFilter = nabs(aod::collision::posZ) < evSel.cfgCutVertex; + Filter acceptanceFilter = (nabs(aod::track::eta) < trackCuts.cfgCutEta && nabs(aod::track::pt) > trackCuts.cfgCutPt); - using EventCandidates = soa::Filtered>; + using EventCandidates = soa::Filtered>; using TrackCandidates = soa::Filtered>; using TrackCandidatesMC = soa::Filtered>; + using EventCandidatesMC = soa::Join; + using McCollisionsCent = soa::Join; Preslice perCol = aod::track::collisionId; Preslice perColMC = aod::track::collisionId; + Preslice perMcCollisionDelta = aod::mcparticle::mcCollisionId; // ── Event-mixing binning policies (one per centrality estimator) ──────────── using BinningTypeFT0M = ColumnBinningPolicy; @@ -1046,17 +1220,17 @@ struct DeltaAnalysis { using BinningTypeFV0A = ColumnBinningPolicy; using BinningTypeNTPV = ColumnBinningPolicy; - BinningTypeFT0M binningFT0M{{cfgVtxAxis, cfgCentAxis}, true}; - BinningTypeFT0A binningFT0A{{cfgVtxAxis, cfgCentAxis}, true}; - BinningTypeFT0C binningFT0C{{cfgVtxAxis, cfgCentAxis}, true}; - BinningTypeFV0A binningFV0A{{cfgVtxAxis, cfgCentAxis}, true}; - BinningTypeNTPV binningNTPV{{cfgVtxAxis, cfgCentAxis}, true}; + BinningTypeFT0M binningFT0M{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + BinningTypeFT0A binningFT0A{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + BinningTypeFT0C binningFT0C{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + BinningTypeFV0A binningFV0A{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + BinningTypeNTPV binningNTPV{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; - SameKindPair pairFT0M{binningFT0M, cfgNoMixedEvents, -1, &cache}; - SameKindPair pairFT0A{binningFT0A, cfgNoMixedEvents, -1, &cache}; - SameKindPair pairFT0C{binningFT0C, cfgNoMixedEvents, -1, &cache}; - SameKindPair pairFV0A{binningFV0A, cfgNoMixedEvents, -1, &cache}; - SameKindPair pairNTPV{binningNTPV, cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFT0M{binningFT0M, mixingCfg.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFT0A{binningFT0A, mixingCfg.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFT0C{binningFT0C, mixingCfg.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairFV0A{binningFV0A, mixingCfg.cfgNoMixedEvents, -1, &cache}; + SameKindPair pairNTPV{binningNTPV, mixingCfg.cfgNoMixedEvents, -1, &cache}; std::vector mProtonPool; std::vector mPionPool; @@ -1066,7 +1240,7 @@ struct DeltaAnalysis { aod::BCs const&) { for (auto const& collision : collisions) { - if (!passesEventSelection(collision)) { + if (!passesEventSelectionImpl(collision)) { continue; } const float centrality = getCentrality(collision); @@ -1116,7 +1290,6 @@ struct DeltaAnalysis { } PROCESS_SWITCH(DeltaAnalysis, processSameEvent, "Process same event", true); - // ── Shared mixed-event logic ───────────────────────────────────────────── template void runMixedEvent(PairType& mixingPair) { @@ -1135,7 +1308,7 @@ struct DeltaAnalysis { void processMixedEvent(EventCandidates const&, TrackCandidates const&) { - switch (cfgCentralityEstimator) { + switch (evSel.cfgCentralityEstimator) { case delta_analysis::kFT0M: runMixedEvent(pairFT0M); break; @@ -1158,22 +1331,30 @@ struct DeltaAnalysis { } PROCESS_SWITCH(DeltaAnalysis, processMixedEvent, "Process mixed event", true); - void processMC(soa::Join const& collisions, aod::BCs const&, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles) + // ===================================================================================== + // processMC(): RECONSTRUCTED-ONLY. Reconstructs Delta candidates, performs truth matching, + // fills reconstructed QA, fills reconstructed Delta histograms (RecoDelta - the A x eps + // numerator). Nothing generated-only is touched here. UNCHANGED from before the refactor. + // ===================================================================================== + void processMC(EventCandidatesMC const& collisions, aod::BCs const&, TrackCandidatesMC const& tracks, aod::McParticles const& mcParticles) { - constexpr float kGenCentrality = 1.f; - for (auto const& collision : collisions) { - if (!passesEventSelection(collision)) { + if (!passesEventSelectionImpl(collision)) { continue; } + const float centrality = getCentrality(collision); + + // Nrec - "the total number of selected reconstructed INEL>0 events". + histos.fill(HIST("MCRecoEvent/hRecoEvents"), centrality); + histos.fill(HIST("Event/hNcontributor"), collision.numContrib()); histos.fill(HIST("Event/hVtxZ"), collision.posZ()); const uint64_t collIdx = collision.globalIndex(); auto perColTracks = tracks.sliceBy(perColMC, collIdx); perColTracks.bindExternalIndices(&tracks); - + perColTracks.bindExternalIndices(&mcParticles); for (auto const& t0 : perColTracks) { if (!passesBasicTrackSelection(t0) || !t0.has_mcParticle()) { continue; @@ -1216,7 +1397,6 @@ struct DeltaAnalysis { } } - // fill QAMC for every reconstructed MC track passing basic+PID+DCA cuts, irrespective of truth matching. for (auto const& trackForQAMC : perColTracks) { const float momQAMC = RecoDecay::p(trackForQAMC.px(), trackForQAMC.py(), trackForQAMC.pz()); if (!passesBasicTrackSelection(trackForQAMC)) { @@ -1246,7 +1426,6 @@ struct DeltaAnalysis { continue; } - // reconstructed-level (pre-truth-matching) invariant mass. { const std::array, 2> momentaReco = { std::array{t0.px(), t0.py(), t0.pz()}, @@ -1254,7 +1433,7 @@ struct DeltaAnalysis { const float pairMassReco = RecoDecay::m(momentaReco, std::array{massProton, massPion}); const float pairPtReco = RecoDecay::pt(std::array{t0.px() + t1.px(), t0.py() + t1.py()}); const float pairYReco = RecoDecay::y(std::array{t0.px() + t1.px(), t0.py() + t1.py(), t0.pz() + t1.pz()}, pairMassReco); - if (pairYReco >= cfgMinY && pairYReco <= cfgMaxY) { + if (pairYReco >= trackCuts.cfgMinY && pairYReco <= trackCuts.cfgMaxY) { fillDeltaHistogramMCReco(t0.sign(), t1.sign(), pairPtReco, pairMassReco, centrality, pairYReco); } } @@ -1296,26 +1475,67 @@ struct DeltaAnalysis { const float pairMass = RecoDecay::m(momenta, std::array{massProton, massPion}); const float pairPt = RecoDecay::pt(std::array{t0.px() + t1.px(), t0.py() + t1.py()}); const float pairY = RecoDecay::y(std::array{t0.px() + t1.px(), t0.py() + t1.py(), t0.pz() + t1.pz()}, pairMass); - if (pairY < cfgMinY || pairY > cfgMaxY) { + if (pairY < trackCuts.cfgMinY || pairY > trackCuts.cfgMaxY) { continue; } fillDeltaHistogramMC(t0.sign(), t1.sign(), pairPt, pairMass, centrality, pairY); } } + } + PROCESS_SWITCH(DeltaAnalysis, processMC, "Process MC", false); + + void processMCGen(McCollisionsCent::iterator const& mcCollision, + aod::McParticles const& mcParticles, + soa::SmallGroups const& collisions) + { + histos.fill(HIST("CutFlow/MCGen/hEventCutFlow"), 0.f); // All MC collisions + + if (std::abs(mcCollision.posZ()) > evSel.cfgCutVertex) { + return; + } + histos.fill(HIST("CutFlow/MCGen/hEventCutFlow"), 1.f); // |Vz| < cfgCutVertex + + const bool truthInelGt0 = isTruthInelGt0(mcParticles); + if (evSel.cfgRequireRecoINELgt0 && !truthInelGt0) { + return; + } + histos.fill(HIST("CutFlow/MCGen/hEventCutFlow"), 2.f); // Truth INEL>0 + histos.fill(HIST("CutFlow/MCGen/hEventCutFlow"), 3.f); // Centrality (no cut currently applied) + histos.fill(HIST("CutFlow/MCGen/hEventCutFlow"), 4.f); // Final generated event + + bool hasAcceptedReco = false; + float genCentrality = mcCollision.centFT0M(); // fallback: MC-truth centrality proxy (see note above) + for (auto const& collision : collisions) { + if (passesEventSelection(collision)) { + hasAcceptedReco = true; + genCentrality = getCentrality(collision); // real reconstructed centrality of an accepted associated collision + break; + } + } for (auto const& mcParticle : mcParticles) { + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 0.f); // bin0: all generated particles + + if (!mcParticle.producedByGenerator()) { + continue; + } + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 1.f); // bin1: producedByGenerator + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 2.f); // bin2: truth vertex-z (event-level, already enforced above) + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 3.f); // bin3: truth INEL>0 (event-level, already enforced above) + const int pdg = mcParticle.pdgCode(); if (std::abs(pdg) != delta_analysis::PdgDeltaPlusPlus && std::abs(pdg) != delta_analysis::PdgDeltaZero) { continue; } - if (mcParticle.y() < cfgMinY || mcParticle.y() > cfgMaxY) { + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 4.f); // bin4: PDG + + if (mcParticle.y() < trackCuts.cfgMinY || mcParticle.y() > trackCuts.cfgMaxY) { continue; } + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 5.f); // bin5: rapidity + const auto daughters = mcParticle.daughters_as(); - // CHANGED: daughter loop is now used only for decay-channel validation / QA - // (hasPr, hasPi, ptPr, ptPi). It no longer accumulates daughter energy for - // the mass calculation. bool hasPr = false, hasPi = false; float ptPr = -999.f, ptPi = -999.f; for (const auto& d : daughters) { @@ -1330,7 +1550,9 @@ struct DeltaAnalysis { if (!hasPr || !hasPi) { continue; } - // CHANGED: genMass is now computed directly from the mother particle's own four-momentum (E, px, py, pz) + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 6.f); // bin6: decay channel + histos.fill(HIST("MCGenQA/hGenDeltaCutFlow"), 7.f); // bin7: final generated histogram + const float eMother = mcParticle.e(); const float pxMother = mcParticle.px(); const float pyMother = mcParticle.py(); @@ -1339,32 +1561,106 @@ struct DeltaAnalysis { const float genMass = std::sqrt(std::max(0.f, mass2)); const float genPt = mcParticle.pt(); const float genY = mcParticle.y(); + + // GeneratedDelta_All if (pdg == delta_analysis::PdgDeltaPlusPlus) { - histos.fill(HIST("Analysis/hDeltaPlusPlusInvMassGen"), genPt, genMass); - // CHANGED: THnSparse Gen fill now carries genMass as the first argument, - // matching the new {massAxis, ptAxis, centAxis, rapAxis} axis order. - histos.fill(HIST("THnSparse/hDeltaPlusPlusGen"), genMass, genPt, kGenCentrality, genY); - histos.fill(HIST("QAChecks/hGenProtonDeltaPlusPlus"), ptPr); - histos.fill(HIST("QAChecks/hGenPionDeltaPlusPlus"), ptPi); + histos.fill(HIST("MCGen/hDeltaPlusPlusInvMassGen"), genPt, genMass); + histos.fill(HIST("MCGen/hDeltaPlusPlusGenSparse"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/hGenProtonDeltaPlusPlus"), ptPr); + histos.fill(HIST("MCGenQA/hGenPionDeltaPlusPlus"), ptPi); } else if (pdg == -delta_analysis::PdgDeltaPlusPlus) { - histos.fill(HIST("Analysis/hAntiDeltaPlusPlusInvMassGen"), genPt, genMass); - histos.fill(HIST("THnSparse/hAntiDeltaPlusPlusGen"), genMass, genPt, kGenCentrality, genY); - histos.fill(HIST("QAChecks/hGenProtonAntiDeltaPlusPlus"), ptPr); - histos.fill(HIST("QAChecks/hGenPionAntiDeltaPlusPlus"), ptPi); + histos.fill(HIST("MCGen/hAntiDeltaPlusPlusInvMassGen"), genPt, genMass); + histos.fill(HIST("MCGen/hAntiDeltaPlusPlusGenSparse"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/hGenProtonAntiDeltaPlusPlus"), ptPr); + histos.fill(HIST("MCGenQA/hGenPionAntiDeltaPlusPlus"), ptPi); } else if (pdg == delta_analysis::PdgDeltaZero) { - histos.fill(HIST("Analysis/hDeltaZeroInvMassGen"), genPt, genMass); - histos.fill(HIST("THnSparse/hDeltaZeroGen"), genMass, genPt, kGenCentrality, genY); - histos.fill(HIST("QAChecks/hGenProtonDeltaZero"), ptPr); - histos.fill(HIST("QAChecks/hGenPionDeltaZero"), ptPi); + histos.fill(HIST("MCGen/hDeltaZeroInvMassGen"), genPt, genMass); + histos.fill(HIST("MCGen/hDeltaZeroGenSparse"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/hGenProtonDeltaZero"), ptPr); + histos.fill(HIST("MCGenQA/hGenPionDeltaZero"), ptPi); } else if (pdg == -delta_analysis::PdgDeltaZero) { - histos.fill(HIST("Analysis/hAntiDeltaZeroInvMassGen"), genPt, genMass); - histos.fill(HIST("THnSparse/hAntiDeltaZeroGen"), genMass, genPt, kGenCentrality, genY); - histos.fill(HIST("QAChecks/hGenProtonAntiDeltaZero"), ptPr); - histos.fill(HIST("QAChecks/hGenPionAntiDeltaZero"), ptPi); + histos.fill(HIST("MCGen/hAntiDeltaZeroInvMassGen"), genPt, genMass); + histos.fill(HIST("MCGen/hAntiDeltaZeroGenSparse"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/hGenProtonAntiDeltaZero"), ptPr); + histos.fill(HIST("MCGenQA/hGenPionAntiDeltaZero"), ptPi); + } + + // GeneratedDelta_EventAccepted (A x eps denominator) + if (hasAcceptedReco) { + if (pdg == delta_analysis::PdgDeltaPlusPlus) { + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hDeltaPlusPlusInvMass"), genPt, genMass); + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hDeltaPlusPlusSparse"), genMass, genPt, genCentrality, genY); + } else if (pdg == -delta_analysis::PdgDeltaPlusPlus) { + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaPlusPlusInvMass"), genPt, genMass); + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaPlusPlusSparse"), genMass, genPt, genCentrality, genY); + } else if (pdg == delta_analysis::PdgDeltaZero) { + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hDeltaZeroInvMass"), genPt, genMass); + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hDeltaZeroSparse"), genMass, genPt, genCentrality, genY); + } else if (pdg == -delta_analysis::PdgDeltaZero) { + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaZeroInvMass"), genPt, genMass); + histos.fill(HIST("MCGen/GeneratedDelta_EventAccepted/hAntiDeltaZeroSparse"), genMass, genPt, genCentrality, genY); + } } } } - PROCESS_SWITCH(DeltaAnalysis, processMC, "Process MC", false); + PROCESS_SWITCH(DeltaAnalysis, processMCGen, "Process MC truth: generated-level Delta spectra and QA", false); + + void processEventFactor(EventCandidatesMC const& collisions, + McCollisionsCent const& mcCollisions, + aod::McParticles const& mcParticles) + { + // ── Loop A: reconstructed collisions ────────────────────────────────────────────────── + std::unordered_map nRecoCollisionsPerMc; + std::unordered_set acceptedMcCollisionIds; + for (auto const& collision : collisions) { + if (collision.has_mcCollision()) { + ++nRecoCollisionsPerMc[collision.mcCollisionId()]; + } + if (!passesEventSelection(collision)) { + continue; + } + const float centrality = getCentrality(collision); + histos.fill(HIST("EventFactor/hRecoEvents"), centrality); + if (collision.has_mcCollision()) { + acceptedMcCollisionIds.insert(collision.mcCollisionId()); + } + } + + // ── Loop B: generated MC collisions ─────────────────────────────────────────────────── + for (auto const& mcCollision : mcCollisions) { + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 0.f); // All generated MC collisions + + if (std::abs(mcCollision.posZ()) > evSel.cfgCutVertex) { + continue; + } + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 1.f); // |Vz| < cfgCutVertex + + auto mcPartsThisColl = mcParticles.sliceBy(perMcCollisionDelta, mcCollision.globalIndex()); + if (evSel.cfgRequireRecoINELgt0 && !isTruthInelGt0(mcPartsThisColl)) { + continue; + } + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 2.f); // Truth INEL>0 + + const float truthCentrality = mcCollision.centFT0M(); // MC-truth centrality proxy (see note above) + histos.fill(HIST("EventFactor/hEventsGenAll"), truthCentrality); + + const auto itCount = nRecoCollisionsPerMc.find(mcCollision.globalIndex()); + const int nRecoTotal = (itCount != nRecoCollisionsPerMc.end()) ? itCount->second : 0; + histos.fill(HIST("EventFactor/hNRecoCollisionsPerMcCollision"), nRecoTotal); + + if (nRecoTotal > 0) { + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 3.f); // Has associated reconstructed collision + } + + const bool hasAcceptedReco = acceptedMcCollisionIds.find(mcCollision.globalIndex()) != acceptedMcCollisionIds.end(); + if (hasAcceptedReco) { + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 4.f); // Associated reco collision passes event selection + histos.fill(HIST("CutFlow/EventFactor/hEventAcceptedCutFlow"), 5.f); // Final EventAccepted + histos.fill(HIST("EventFactor/hEventsGenAccepted"), truthCentrality); + } + } + } + PROCESS_SWITCH(DeltaAnalysis, processEventFactor, "Process Event Loss / Event Splitting factors", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& context)