diff --git a/PWGLF/Tasks/Resonances/deltaAnalysis.cxx b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx index 81f28e2b513..951f76b15c5 100644 --- a/PWGLF/Tasks/Resonances/deltaAnalysis.cxx +++ b/PWGLF/Tasks/Resonances/deltaAnalysis.cxx @@ -48,6 +48,7 @@ #include #include #include +#include #include #include #include @@ -127,7 +128,8 @@ struct DeltaAnalysis { Configurable applyTPCChi2NClCut{"applyTPCChi2NClCut", true, "Enable/disable maximum TPC chi2/NCl cut"}; Configurable applyITSChi2NClCut{"applyITSChi2NClCut", true, "Enable/disable maximum ITS chi2/NCl cut"}; Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require isPrimaryTrack flag"}; - Configurable requireGlobalTrack{"requireGlobalTrack", true, "Require global track"}; + Configurable requireGlobalTrack{"requireGlobalTrack", true, "Require global track (includes the DCA cuts of the standard global-track selection)"}; + Configurable requireGlobalTrackWoDCA{"requireGlobalTrackWoDCA", false, "Require global track without DCA cuts (isGlobalTrackWoDCA). Mutually exclusive with requireGlobalTrack"}; Configurable requirePVContributor{"requirePVContributor", true, "Require PV-contributor flag"}; } trackCuts; @@ -226,6 +228,11 @@ struct DeltaAnalysis { << ", bins=" << binsVec.size() << ")"; } }; + if (trackCuts.requireGlobalTrack && trackCuts.requireGlobalTrackWoDCA) { + LOG(fatal) << "requireGlobalTrack and requireGlobalTrackWoDCA are mutually exclusive: " + << "enable only one of them (or neither)."; + } + chkSz(protonPID.protonTPCNSigmaCutPerBin, protonPID.protonTPCPIDMomentumBins, "proton TPC PID"); chkSz(protonPID.protonTOFNSigmaCutPerBin, protonPID.protonTOFPIDMomentumBins, "proton TOF PID"); chkSz(pionPID.pionTPCNSigmaCutPerBin, pionPID.pionTPCPIDMomentumBins, "pion TPC PID"); @@ -246,15 +253,15 @@ struct DeltaAnalysis { mPionDCAPtEdges = dcaCuts.pionDCAPtBinEdges.value; mPionMaxDCAxy = dcaCuts.pionMaxDCAxyPerPtBin.value; - const AxisSpec ptAxis{50, 0., 20., "p_{T} (GeV/c)"}; + const AxisSpec ptAxis{50, 0., 12., "p_{T} (GeV/c)"}; const AxisSpec massAxis{trackCuts.numberOfInvMassBins, 1.0, 6.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, -5., 5., "n#sigma^{TPC}"}; const AxisSpec nSigmaTOFaxis{100, -5., 5., "n#sigma^{TOF}"}; - const AxisSpec momentumAxis{50, 0., 20., "p (GeV/#it{c})"}; - const AxisSpec ptForPIDAxis{50, 0., 20., "p_{T} (GeV/#it{c})"}; + const AxisSpec momentumAxis{50, 0., 12., "p (GeV/#it{c})"}; + const AxisSpec ptForPIDAxis{50, 0., 12., "p_{T} (GeV/#it{c})"}; const AxisSpec dcaXYaxis{200, -0.2f, 0.2f, "DCA_{xy} (cm)"}; const AxisSpec dcaZaxis{200, -1.2f, 1.2f, "DCA_{z} (cm)"}; const AxisSpec tpcClusAxis{200, 0, 200, "N_{clusters}^{TPC}"}; @@ -757,6 +764,9 @@ struct DeltaAnalysis { if (trackCuts.requireGlobalTrack && !track.isGlobalTrack()) { return false; } + if (trackCuts.requireGlobalTrackWoDCA && !track.isGlobalTrackWoDCA()) { + return false; + } if (trackCuts.requirePVContributor && !track.isPVContributor()) { return false; } @@ -1429,18 +1439,6 @@ struct DeltaAnalysis { using BinningTypeFV0A = ColumnBinningPolicy; using BinningTypeNTPV = ColumnBinningPolicy; - BinningTypeFT0M binningFT0M{{axes.cfgVtxAxis, axes.cfgCentAxis}}; - BinningTypeFT0A binningFT0A{{axes.cfgVtxAxis, axes.cfgCentAxis}}; - BinningTypeFT0C binningFT0C{{axes.cfgVtxAxis, axes.cfgCentAxis}}; - BinningTypeFV0A binningFV0A{{axes.cfgVtxAxis, axes.cfgCentAxis}}; - BinningTypeNTPV binningNTPV{{axes.cfgVtxAxis, axes.cfgCentAxis}}; - - 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; @@ -1524,31 +1522,6 @@ struct DeltaAnalysis { } } - void processMixedEvent(EventCandidates const&, TrackCandidates const&) - { - switch (evSel.cfgCentralityEstimator) { - case delta_analysis::kFT0M: - runMixedEvent(pairFT0M); - break; - case delta_analysis::kFT0A: - runMixedEvent(pairFT0A); - break; - case delta_analysis::kFT0C: - runMixedEvent(pairFT0C); - break; - case delta_analysis::kFV0A: - runMixedEvent(pairFV0A); - break; - case delta_analysis::kNTPV: - runMixedEvent(pairNTPV); - break; - default: - runMixedEvent(pairFT0M); - break; - } - } - PROCESS_SWITCH(DeltaAnalysis, processMixedEvent, "Process mixed event", true); - // ── MC reconstructed event mixing ──────────────────────────────────────────────────────── // Mirrors runMixedEvent()/processMixedEvent() above exactly (same event-selection call, same // "protons1+pions2" / "protons2+pions1" combinatorics, same centrality-estimator dispatch), @@ -1574,32 +1547,57 @@ struct DeltaAnalysis { } } - void processMCEventMixing(EventCandidates const&, TrackCandidates const&) + template + void runMixingWithBinning(EventCandidates const& collisions, TrackCandidates const& tracks, bool isMC) { - if (!mixingCfg.enableMCEventMixing) { - return; + BinningType binning{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + const int nMixedEvents = mixingCfg.cfgNoMixedEvents; + auto tracksTuple = std::make_tuple(tracks); + SameKindPair mixingPair{binning, nMixedEvents, -1, collisions, tracksTuple, &cache}; + if (isMC) { + runMixedEventMC(mixingPair); + } else { + runMixedEvent(mixingPair); } + } + + void dispatchMixing(EventCandidates const& collisions, TrackCandidates const& tracks, bool isMC) + { switch (evSel.cfgCentralityEstimator) { case delta_analysis::kFT0M: - runMixedEventMC(pairFT0M); + runMixingWithBinning(collisions, tracks, isMC); break; case delta_analysis::kFT0A: - runMixedEventMC(pairFT0A); + runMixingWithBinning(collisions, tracks, isMC); break; case delta_analysis::kFT0C: - runMixedEventMC(pairFT0C); + runMixingWithBinning(collisions, tracks, isMC); break; case delta_analysis::kFV0A: - runMixedEventMC(pairFV0A); + runMixingWithBinning(collisions, tracks, isMC); break; case delta_analysis::kNTPV: - runMixedEventMC(pairNTPV); + runMixingWithBinning(collisions, tracks, isMC); break; default: - runMixedEventMC(pairFT0M); + runMixingWithBinning(collisions, tracks, isMC); break; } } + + void processMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) + { + dispatchMixing(collisions, tracks, false); + } + PROCESS_SWITCH(DeltaAnalysis, processMixedEvent, "Process mixed event", true); + + void processMCEventMixing(EventCandidates const& collisions, TrackCandidates const& tracks) + { + if (!mixingCfg.enableMCEventMixing) { + return; + } + dispatchMixing(collisions, tracks, true); + } PROCESS_SWITCH(DeltaAnalysis, processMCEventMixing, "Process MC reconstructed event mixing (reuses DATA mixing pairs; no truth matching)", false); // ===================================================================================== diff --git a/PWGLF/Tasks/Resonances/lambda1520analysispo.cxx b/PWGLF/Tasks/Resonances/lambda1520analysispo.cxx index b2e09900706..142ddfd3a33 100644 --- a/PWGLF/Tasks/Resonances/lambda1520analysispo.cxx +++ b/PWGLF/Tasks/Resonances/lambda1520analysispo.cxx @@ -13,16 +13,22 @@ /// \brief Task for Lambda(1520) resonance reconstruction via proton-kaon invariant mass analysis for pO system /// /// \author Durgesh Bhatt +/// \author Narayan khadanga +#include "PWGLF/DataModel/mcCentrality.h" -#include "PWGLF/DataModel/LFResonanceTables.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" #include #include -#include -#include #include +#include #include #include #include @@ -31,1176 +37,2554 @@ #include #include #include +#include #include #include +#include +#include + +#include +#include #include +#include #include +#include +#include +#include +#include +#include +#include #include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::constants::physics; -namespace lambda_analysis +namespace +{ +constexpr float ProtonMass = o2::constants::physics::MassProton; +constexpr float KaonMass = o2::constants::physics::MassKaonCharged; +constexpr float LambdaMass = o2::constants::physics::MassLambda0; +constexpr float XiMass = o2::constants::physics::MassXiMinus; +constexpr float OmegaMass = o2::constants::physics::MassOmegaMinus; +constexpr float Lambda1520Mass = o2::constants::physics::MassLambda1520; +constexpr float MinAbsCharge = 1e-3f; +constexpr float InvalidCentrality = -999.0f; +const std::array kProtonKaonMasses{ProtonMass, KaonMass}; +} // namespace + +namespace lambda1520analysispo_analysis { -// Named PDG codes — fixes pdg/explicit-code and magic-number errors static constexpr int PdgProton{2212}; static constexpr int PdgKaon{321}; static constexpr int PdgLambda0{3122}; static constexpr int PdgXiMinus{3312}; -static constexpr int PdgXi0{3322}; static constexpr int PdgOmegaMinus{3334}; -} // namespace lambda_analysis - -struct lambda1520analysispo { - SliceCache sliceCache; - - // Preslice helpers - Preslice tracksPerResonanceCollision = aod::resodaughter::resoCollisionId; - Preslice tracksPerStandardCollision = aod::track::collisionId; - - // Pointer to MC parent particle table (used only in MC processing) - aod::ResoMCParents const* mcResonanceParentTable = nullptr; - - // ── Event-level configurables ──────────────────────────────────────────── - // FIX name/configurable: lowerCamelCase, member name == JSON string, on one - // line - Configurable applyOccupancyInTimeRangeCut{"applyOccupancyInTimeRangeCut", false, "If true, apply a cut on the number of tracks in a time window around the collision (occupancy cut)"}; - - // ── Histogram binning configurables ───────────────────────────────────── - Configurable numberOfPtBins{"numberOfPtBins", 100, "Number of bins along the transverse momentum (pT) axis"}; - Configurable numberOfInvMassBins{"numberOfInvMassBins", 120, "Number of bins along the invariant mass axis"}; - - // ── Physics configurables ──────────────────────────────────────────────── - Configurable pdgCodeLambda1520{"pdgCodeLambda1520", 3124, "PDG code of the Lambda(1520) resonance"}; - Configurable enableRotationalBackground{"enableRotationalBackground", true, "If true, compute rotational background (kaon phi rotated by ~pi) to estimate combinatorial background"}; - - // ── Track quality cuts ─────────────────────────────────────────────────── - Configurable minTrackPt{"minTrackPt", 0.15f, "Minimum transverse momentum pT of a track [GeV/c]"}; - Configurable minProtonMomentum{"minProtonMomentum", 0.f, "Minimum total momentum p for proton PID [GeV/c]"}; - Configurable minKaonMomentum{"minKaonMomentum", 0.f, "Minimum total momentum p for kaon PID [GeV/c]"}; - Configurable minPseudorapidity{"minPseudorapidity", -0.8f, "Minimum pseudorapidity eta (detector acceptance)"}; - Configurable maxPseudorapidity{"maxPseudorapidity", 0.8f, "Maximum pseudorapidity eta (detector acceptance)"}; - Configurable maxDCAz{"maxDCAz", 1.0f, "Maximum allowed distance of closest approach along z-axis (DCAz) [cm]"}; - Configurable minPairRapidity{"minPairRapidity", -0.5f, "Minimum rapidity y of the reconstructed proton-kaon pair"}; - Configurable maxPairRapidity{"maxPairRapidity", 0.5f, "Maximum rapidity y of the reconstructed proton-kaon pair"}; - - // ── TPC cluster quality cuts ───────────────────────────────────────────── - Configurable minTPCCrossedRows{"minTPCCrossedRows", 70, "Minimum number of TPC crossed pad rows (track quality)"}; - Configurable applyCrossedRowsCut{"applyCrossedRowsCut", false, "If true, require at least minTPCCrossedRows crossed rows in the TPC"}; - Configurable minTPCClustersFound{"minTPCClustersFound", 70, "Minimum number of TPC clusters found on the track"}; - Configurable applyTPCClustersCut{"applyTPCClustersCut", false, "If true, require at least minTPCClustersFound TPC clusters"}; - - // ── pT-dependent DCAxy cuts for Protons ───────────────────────────────── - Configurable> protonDCAPtBinEdges{"protonDCAPtBinEdges", {0.0f, 0.5f, 1.0f, 2.0f, 3.0f, 5.0f, 1000.0f}, "Proton pT bin edges [GeV/c] for the pT-dependent DCAxy selection"}; - Configurable> protonMaxDCAxyPerPtBin{"protonMaxDCAxyPerPtBin", {0.020f, 0.015f, 0.010f, 0.007f, 0.005f, 0.004f}, "Maximum |DCAxy| [cm] for protons in each pT bin defined by protonDCAPtBinEdges"}; - - // ── pT-dependent DCAxy cuts for Kaons ─────────────────────────────────── - Configurable> kaonDCAPtBinEdges{"kaonDCAPtBinEdges", {0.0f, 0.3f, 0.6f, 1.0f, 2.0f, 1000.0f}, "Kaon pT bin edges [GeV/c] for the pT-dependent DCAxy selection"}; - Configurable> kaonMaxDCAxyPerPtBin{"kaonMaxDCAxyPerPtBin", {0.025f, 0.018f, 0.012f, 0.008f, 0.004f}, "Maximum |DCAxy| [cm] for kaons in each pT bin defined by kaonDCAPtBinEdges"}; - - // ── Analysis mode switches ─────────────────────────────────────────────── - Configurable runQualityChecksOnly{"runQualityChecksOnly", false, "If true, only fill QA histograms and skip invariant mass computation"}; - Configurable applyDeepAngleCut{"applyDeepAngleCut", false, "If true, reject proton-kaon pairs with very small opening angle (removes split-track background)"}; - Configurable deepAngleCutValue{"deepAngleCutValue", 0.04, "Minimum allowed opening angle [rad] between proton and kaon (used if applyDeepAngleCut = true)"}; - Configurable applyKinematicPairCuts{"applyKinematicPairCuts", false, "If true, apply additional kinematic cuts on the p-K opening angle"}; - - // ── Global track selection flags ───────────────────────────────────────── - Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require track to pass the 'isPrimaryTrack' flag (kGoldenChi2 | kDCAxy | kDCAz)"}; - Configurable requireGlobalTrackNoDCA{"requireGlobalTrackNoDCA", true, "Require track to pass 'isGlobalTrackWoDCA' (quality cuts without DCA requirement)"}; - Configurable requirePVContributor{"requirePVContributor", true, "Require track to be a contributor to the primary vertex reconstruction"}; - - // ── PID configurables ──────────────────────────────────────────────────── - Configurable requireTOFForProton{"requireTOFForProton", false, "If true, only accept proton candidates that have a TOF measurement"}; - Configurable requireTOFForKaon{"requireTOFForKaon", false, "If true, only accept kaon candidates that have a TOF measurement"}; - Configurable useTPCOnlyPID{"useTPCOnlyPID", false, "If true, use only TPC nSigma for PID (ignore TOF even if available)"}; - - Configurable tpcNSigmaVetoOtherSpecies{"tpcNSigmaVetoOtherSpecies", 3.0f, "Reject track if its TPC nSigma for a different species is below this threshold (avoids misID)"}; - Configurable tpcNSigmaVetoPion{"tpcNSigmaVetoPion", 3.0f, "TPC nSigma threshold below which a track is vetoed as a pion"}; - Configurable tpcNSigmaVetoKaonForProton{"tpcNSigmaVetoKaonForProton", 3.0f, "TPC nSigma threshold: veto proton candidates that look like kaons"}; - Configurable tpcNSigmaVetoPionForKaon{"tpcNSigmaVetoPionForKaon", 3.0f, "TPC nSigma threshold: veto kaon candidates that look like pions"}; - Configurable tpcNSigmaVetoProtonForKaon{"tpcNSigmaVetoProtonForKaon", 3.0f, "TPC nSigma threshold: veto kaon candidates that look like protons"}; - - Configurable minTPCNSigmaKaon{"minTPCNSigmaKaon", -6.0f, "Minimum (most negative) TPC nSigma for kaon"}; - Configurable minTPCNSigmaProton{"minTPCNSigmaProton", -6.0f, "Minimum (most negative) TPC nSigma for proton"}; - Configurable minTOFNSigmaKaon{"minTOFNSigmaKaon", -6.0f, "Minimum (most negative) TOF nSigma for kaon"}; - Configurable minTOFNSigmaProton{"minTOFNSigmaProton", -6.0f, "Minimum (most negative) TOF nSigma for proton"}; - Configurable minCombinedTPCTOFNSigmaKaon{"minCombinedTPCTOFNSigmaKaon", -6.0f, "Minimum combined TPC+TOF nSigma for kaon (used in combined PID mode)"}; - Configurable minCombinedTPCTOFNSigmaProton{"minCombinedTPCTOFNSigmaProton", -6.0f, "Minimum combined TPC+TOF nSigma for proton (used in combined PID mode)"}; - - Configurable tpcNSigmaVetoThreshold{"tpcNSigmaVetoThreshold", 3.0f, "General TPC nSigma veto cut to reject misidentified particles when TPC bands overlap"}; - Configurable tofNSigmaVetoThreshold{"tofNSigmaVetoThreshold", 3.0f, "General TOF nSigma veto cut to reject misidentified particles"}; - - // ── Proton PID momentum-dependent TPC cuts ─────────────────────────────── - Configurable maxTPCNSigmaProton{"maxTPCNSigmaProton", 3.0, "Maximum |TPC nSigma| for proton identification (symmetric cut)"}; - Configurable combinedNSigmaCutProton{"combinedNSigmaCutProton", 3.0, "Cut on sqrt(nSigmaTPC^2 + nSigmaTOF^2) for proton. Negative value switches to asymmetric mode."}; - Configurable> protonTPCPIDMomentumBins{"protonTPCPIDMomentumBins", {0, 0.5, 0.7, 0.8}, "Momentum p bin edges [GeV/c] for momentum-dependent TPC PID cuts on protons"}; - Configurable> protonTPCNSigmaCutPerBin{"protonTPCNSigmaCutPerBin", {5., 3.5, 2.5}, "Maximum TPC nSigma for proton in each momentum bin (tighter at higher p)"}; - Configurable> protonTOFPIDMomentumBins{"protonTOFPIDMomentumBins", {0., 999.}, "Momentum p bin edges [GeV/c] for momentum-dependent TOF PID cuts on protons"}; - Configurable> protonTOFNSigmaCutPerBin{"protonTOFNSigmaCutPerBin", {3.0}, "Maximum TOF nSigma for proton in each momentum bin"}; - - // ── Kaon PID momentum-dependent TPC cuts ──────────────────────────────── - Configurable maxTPCNSigmaKaon{"maxTPCNSigmaKaon", 3.0, "Maximum |TPC nSigma| for kaon identification (symmetric cut)"}; - Configurable combinedNSigmaCutKaon{"combinedNSigmaCutKaon", 3.0, "Cut on sqrt(nSigmaTPC^2 + nSigmaTOF^2) for kaon. Negative value switches to asymmetric mode."}; - Configurable> kaonTPCPIDMomentumBins{"kaonTPCPIDMomentumBins", {0., 0.25, 0.3, 0.45}, "Momentum p bin edges [GeV/c] for momentum-dependent TPC PID cuts on kaons"}; - Configurable> kaonTPCNSigmaCutPerBin{"kaonTPCNSigmaCutPerBin", {6, 3.5, 2.5}, "Maximum TPC nSigma for kaon in each momentum bin"}; - Configurable> kaonTOFPIDMomentumBins{"kaonTOFPIDMomentumBins", {0., 999.}, "Momentum p bin edges [GeV/c] for momentum-dependent TOF PID cuts on kaons"}; - Configurable> kaonTOFNSigmaCutPerBin{"kaonTOFNSigmaCutPerBin", {3.0}, "Maximum TOF nSigma for kaon in each momentum bin"}; - - // ── Event mixing configurables ─────────────────────────────────────────── - Configurable numberOfEventsToMix{"numberOfEventsToMix", 20, "Number of events to mix with each signal event for background estimation"}; - ConfigurableAxis dcaZMixingBins{"dcaZMixingBins", {VARIABLE_WIDTH, -1.2f, -1.0f, -0.9f, -0.8f, -0.7f, -0.6f, -0.5f, -0.4f, -0.3f, -0.2f, -0.1f, 0.f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.2f}, "DCAz vertex bins used for event mixing pairing"}; - ConfigurableAxis vertexZMixingBins{"vertexZMixingBins", {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}, "Primary vertex z-position bins used for event mixing pairing"}; - ConfigurableAxis centralityMixingBins{"centralityMixingBins", {VARIABLE_WIDTH, 0.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f, 200.f}, "Centrality (FT0 %) bins used for event mixing pairing"}; - ConfigurableAxis eventPlaneMixingBins{"eventPlaneMixingBins", {VARIABLE_WIDTH, -1.5708f, -1.25664f, -0.942478f, -0.628319f, 0.f, 0.628319f, 0.942478f, 1.25664f, 1.5708f}, "Event plane angle bins used for event mixing (EP-dependent analysis)"}; - ConfigurableAxis occupancyBins{"occupancyBins", {VARIABLE_WIDTH, 0.0, 100, 500, 600, 1000, 1100, 1500, 1600, 2000, 2100, 2500, 2600, 3000, 3100, 3500, 3600, 4000, 4100, 4500, 4600, 5000, 5100, 9999}, "Track occupancy bins in the time range around the collision"}; - - // ── Rotational background configurables ───────────────────────────────── - Configurable numberOfRotations{"numberOfRotations", 10, "How many times to rotate the kaon phi for the rotational background estimate"}; - Configurable rotationAngleWindow{"rotationAngleWindow", 6.f, "The kaon is rotated by angles near PI, within a window of PI/rotationAngleWindow"}; - - // ── MC event selection flags ───────────────────────────────────────────── - Configurable mcRequireAfterAllCuts{"mcRequireAfterAllCuts", false, "MC event selection: require isInAfterAllCuts flag"}; - Configurable mcRequireINELgt0{"mcRequireINELgt0", false, "MC event selection: require at least 1 charged particle in |eta|<1 (INEL>0)"}; - Configurable mcRequireSel8{"mcRequireSel8", false, "MC event selection: require the standard Sel8 event selection"}; - Configurable mcRequireVtxWithin10cm{"mcRequireVtxWithin10cm", false, "MC event selection: require primary vertex |z| < 10 cm"}; - Configurable mcRequireTriggerTVX{"mcRequireTriggerTVX", false, "MC event selection: require the TVX (T0 vertex) trigger"}; - Configurable cEvtRecINELgt0{"cEvtRecINELgt0", false, "MC event selection: require reconstructed INEL>0 condition"}; - - // ── Histogram registry ─────────────────────────────────────────────────── - HistogramRegistry allHistograms{"allHistograms", {}, OutputObjHandlingPolicy::AnalysisObject}; - - // ============================================================ - // init() - // ============================================================ - void init(InitContext const&) - { - const AxisSpec axisCentralityPercent(110, 0, 110, "FT0 centrality (%)"); - const AxisSpec axisMomentumForPID(200, 0., 10., "p (GeV/c)"); - const AxisSpec axisPtForPID(200, 0., 10., "p_{T} (GeV/c)"); - const AxisSpec axisPt(numberOfPtBins, 0., 10., "p_{T} (GeV/c)"); - const AxisSpec axisPseudorapidity(40, -1, 1, "#eta"); - const AxisSpec axisDCAxy(240, -0.12, 0.12, "DCA_{xy} (cm)"); - const AxisSpec axisTPCNClusters(200, 0, 200, "TPC N_{clusters}"); - const AxisSpec axisTPCNSigma(401, -10.025, 10.025, "n#sigma^{TPC}"); - const AxisSpec axisTOFNSigma(401, -10.025, 10.025, "n#sigma^{TOF}"); - const AxisSpec axisTPCdEdx(380, 10, 200, "#frac{dE}{dx}"); - const AxisSpec axisVertexZ(120, -12, 12, "v_{z} (cm)"); - const AxisSpec axisEventPlaneAngle(120, -3.14, 3.14, "#theta (rad)"); - const AxisSpec axisInvariantMass(numberOfInvMassBins, 1.4, 3.0, "M_{inv} (GeV/c^{2})"); - - AxisSpec axisOccupancy = {occupancyBins, "Track occupancy [-40,100]"}; - AxisSpec axisDCAz = {dcaZMixingBins, "DCA_{z} (cm)"}; - - allHistograms.add("Event/centralityVsOccupancy", "Collision centrality vs track occupancy", kTH2F, {axisCentralityPercent, axisOccupancy}); - allHistograms.add("Event/primaryVertexZ", "Primary vertex z-position distribution", kTH1F, {{100, -15., 15.}}); - if (doprocessMix || doprocessMixDF || doprocessMixepDF) { - allHistograms.add("Event/mixingBins_centralityVsVtxZVsEventPlane", "Event mixing bin occupancy: centrality vs vtxZ vs event plane", kTH3F, {axisCentralityPercent, axisVertexZ, axisEventPlaneAngle}); - } - - allHistograms.add("QAbefore/trackEta", "Track pseudorapidity (before cuts)", kTH1F, {{50, -1.0, 1.0}}); - allHistograms.add("QAbefore/trackPt", "Track pT (before cuts)", kTH1F, {axisMomentumForPID}); - allHistograms.add("QAbefore/trackPhi", "Track azimuthal angle (before cuts)", kTH1F, {{72, 0, 6.2832}}); - allHistograms.add("QAbefore/trackEtaVsPhi", "Track eta vs phi (before cuts)", kTH2F, {axisPseudorapidity, {72, 0, 6.2832}}); - - allHistograms.add("QAbefore/Proton/tpcNSigmaVsMomentum", "TPC nSigma proton vs momentum (before PID cuts)", kTH2F, {axisMomentumForPID, axisTPCNSigma}); - allHistograms.add("QAbefore/Proton/tofNSigmaVsMomentum", "TOF nSigma proton vs momentum (before PID cuts)", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAbefore/Proton/tofNSigmaVsTPCNSigma", "TOF nSigma vs TPC nSigma proton (before PID cuts)", kTH2F, {axisTPCNSigma, axisTOFNSigma}); - - allHistograms.add("QAbefore/Kaon/tpcNSigmaVsMomentum", "TPC nSigma kaon vs momentum (before PID cuts)", kTH2F, {axisMomentumForPID, axisTPCNSigma}); - allHistograms.add("QAbefore/Kaon/tofNSigmaVsMomentum", "TOF nSigma kaon vs momentum (before PID cuts)", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAbefore/Kaon/tofNSigmaVsTPCNSigma", "TOF nSigma vs TPC nSigma kaon (before PID cuts)", kTH2F, {axisTPCNSigma, axisTOFNSigma}); - - allHistograms.add("QAafter/Proton/ptVsCentrality", "Proton pT vs centrality (after cuts)", kTH2F, {axisPtForPID, axisCentralityPercent}); - allHistograms.add("QAafter/Proton/dcaZVsPt", "Proton DCAz vs pT (after cuts)", kTH2F, {axisPtForPID, axisDCAz}); - allHistograms.add("QAafter/Proton/dcaXYVsPt", "Proton DCAxy vs pT (after cuts)", kTH2F, {axisPtForPID, axisDCAxy}); - allHistograms.add("QAafter/Proton/tpcDedxVsMomentum", "Proton TPC dE/dx signal vs momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTPCdEdx}); - allHistograms.add("QAafter/Proton/tpcNSigmaVsPt", "Proton TPC nSigma vs pT (after cuts)", kTH2F, {axisPtForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Proton/tpcNSigmaPionContamVsPt", "Proton track: TPC nSigma pion contamination check", kTH2F, {axisPtForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Proton/tpcNSigmaKaonContamVsPt", "Proton track: TPC nSigma kaon contamination check", kTH2F, {axisPtForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Proton/tpcNSigmaVsMomentum", "Proton TPC nSigma vs total momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Proton/tpcNSigmaVsCentrality", "Proton TPC nSigma vs centrality", kTH2F, {axisCentralityPercent, axisTPCNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaVsPt", "Proton TOF nSigma vs pT (after cuts)", kTH2F, {axisPtForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaVsMomentum", "Proton TOF nSigma vs total momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaVsCentrality", "Proton TOF nSigma vs centrality", kTH2F, {axisCentralityPercent, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaPionContamVsMomentum", "Proton track: TOF nSigma pion contamination check", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaKaonContamVsMomentum", "Proton track: TOF nSigma kaon contamination check", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tofNSigmaVsTPCNSigma", "Proton TOF nSigma vs TPC nSigma (after cuts)", kTH2F, {axisTPCNSigma, axisTOFNSigma}); - allHistograms.add("QAafter/Proton/tpcCrossedRowsVsPt", "Proton TPC crossed rows vs pT", kTH2F, {axisPtForPID, {200, 0, 200}}); - allHistograms.add("QAafter/Proton/tpcClustersFoundVsPt", "Proton TPC clusters found vs pT", kTH2F, {axisPtForPID, {200, 0, 200}}); - - allHistograms.add("QAafter/Kaon/ptVsCentrality", "Kaon pT vs centrality (after cuts)", kTH2F, {axisPtForPID, axisCentralityPercent}); - allHistograms.add("QAafter/Kaon/dcaZVsPt", "Kaon DCAz vs pT (after cuts)", kTH2F, {axisPtForPID, axisDCAz}); - allHistograms.add("QAafter/Kaon/dcaXYVsPt", "Kaon DCAxy vs pT (after cuts)", kTH2F, {axisPtForPID, axisDCAxy}); - allHistograms.add("QAafter/Kaon/tpcDedxVsMomentum", "Kaon TPC dE/dx signal vs momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTPCdEdx}); - allHistograms.add("QAafter/Kaon/tpcNSigmaPionContamVsPt", "Kaon track: TPC nSigma pion contamination check", kTH2F, {axisPtForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Kaon/tpcNSigmaProtonContamVsMomentum", "Kaon track: TPC nSigma proton contamination check", kTH2F, {axisMomentumForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Kaon/tpcNSigmaVsPt", "Kaon TPC nSigma vs pT (after cuts)", kTH2F, {axisPtForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Kaon/tpcNSigmaVsMomentum", "Kaon TPC nSigma vs total momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTPCNSigma}); - allHistograms.add("QAafter/Kaon/tpcNSigmaVsCentrality", "Kaon TPC nSigma vs centrality", kTH2F, {axisCentralityPercent, axisTPCNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaVsPt", "Kaon TOF nSigma vs pT (after cuts)", kTH2F, {axisPtForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaVsMomentum", "Kaon TOF nSigma vs total momentum (after cuts)", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaVsCentrality", "Kaon TOF nSigma vs centrality", kTH2F, {axisCentralityPercent, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaPionContamVsMomentum", "Kaon track: TOF nSigma pion contamination check", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaProtonContamVsMomentum", "Kaon track: TOF nSigma proton contamination check", kTH2F, {axisMomentumForPID, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tofNSigmaVsTPCNSigma", "Kaon TOF nSigma vs TPC nSigma (after cuts)", kTH2F, {axisTPCNSigma, axisTOFNSigma}); - allHistograms.add("QAafter/Kaon/tpcCrossedRowsVsPt", "Kaon TPC crossed rows vs pT", kTH2F, {axisPtForPID, {200, 0, 200}}); - allHistograms.add("QAafter/Kaon/tpcClustersFoundVsPt", "Kaon TPC clusters found vs pT", kTH2F, {axisPtForPID, {200, 0, 200}}); - - if (!doprocessMC) { - allHistograms.add("Analysis/invMass_UnlikeSign_ProtonPlusKaonMinus", "Invariant mass: p^{+}K^{-} (Lambda(1520) signal)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_UnlikeSign_ProtonMinusKaonPlus", "Invariant mass: p^{-}K^{+} (anti-Lambda(1520) signal)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_LikeSign_ProtonPlusKaonPlus", "Invariant mass: p^{+}K^{+} (like-sign background)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_LikeSign_ProtonMinusKaonMinus", "Invariant mass: p^{-}K^{-} (like-sign background)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Rotated_ProtonPlusKaonMinus", "Invariant mass: rotational background (Lambda(1520))", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Rotated_ProtonMinusKaonPlus", "Invariant mass: rotational background (anti-Lambda(1520))", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Mixed_ProtonPlusKaonMinus", "Invariant mass: mixed events p^{+}K^{-} (background)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Mixed_ProtonMinusKaonPlus", "Invariant mass: mixed events p^{-}K^{+} (background)", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Mixed_LikeSign_PlusPlus", "Invariant mass: mixed events like-sign ++", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - allHistograms.add("Analysis/invMass_Mixed_LikeSign_MinusMinus", "Invariant mass: mixed events like-sign --", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent, axisOccupancy}); - } - if (doprocessMC) { - allHistograms.add("Event/mcEventSelectionCutflow", "MC event selection cutflow (how many events pass each cut)", kTH1F, {{7, 0, 7}}); - allHistograms.add("QAChecks/protonRecoLevelPt", "Reconstructed proton pT (after PID cuts) - MC", kTH1F, {axisPtForPID}); - allHistograms.add("QAChecks/kaonRecoLevelPt", "Reconstructed kaon pT (after PID cuts) - MC", kTH1F, {axisPtForPID}); - allHistograms.add("QAChecks/protonGenLevelPt", "Generated proton pT (truth level) - MC", kTH1F, {axisPtForPID}); - allHistograms.add("QAChecks/kaonGenLevelPt", "Generated kaon pT (truth level) - MC", kTH1F, {axisPtForPID}); - allHistograms.add("Analysis/mcGenerated_Lambda1520", "Generated Lambda(1520) invariant mass vs pT vs centrality", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent}); - allHistograms.add("Analysis/mcGenerated_AntiLambda1520", "Generated anti-Lambda(1520) invariant mass vs pT vs centrality", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent}); - allHistograms.add("Analysis/mcReconstructed_Lambda1520", "Reconstructed Lambda(1520) invariant mass vs pT vs centrality - MC", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent}); - allHistograms.add("Analysis/mcReconstructed_AntiLambda1520", "Reconstructed anti-Lambda(1520) invariant mass vs pT vs centrality - MC", kTHnSparseF, {axisInvariantMass, axisPt, axisCentralityPercent}); - allHistograms.add("Analysis/mcMassResolution_Lambda1520", "Mass resolution (Mreco - Mgen) vs pT - Lambda(1520)", kTHnSparseF, {{200, -0.05, 0.05}, axisPt, axisCentralityPercent}); - allHistograms.add("Analysis/mcMassResolution_AntiLambda1520", "Mass resolution (Mreco - Mgen) vs pT - anti-Lambda(1520)", kTHnSparseF, {{200, -0.05, 0.05}, axisPt, axisCentralityPercent}); +enum CentEstimator : int { + kFT0M = 0, + kFT0A = 1, + kFT0C = 2, + kFV0A = 3 +}; + +} // namespace lambda1520analysispo_analysis + +struct Lambda1520Analysispo { + + SliceCache cache; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + Service pdgDB{}; + + struct : ConfigurableGroup { + Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted |z-vertex| range [cm]"}; + Configurable cfgUseNoITSROFrameBorderCut{"cfgUseNoITSROFrameBorderCut", false, "Apply kNoITSROFrameBorder selection"}; + Configurable cfgUseSel8{"cfgUseSel8", false, "Apply sel8 selection"}; + Configurable cfgRequireRecoINELgt0{"cfgRequireRecoINELgt0", true, "Require reconstructed INEL>0 event class (applied in ALL processes)"}; + Configurable cfgUseNoSameBunchPileupCut{"cfgUseNoSameBunchPileupCut", false, "Apply kNoSameBunchPileup event selection"}; + Configurable cfgUseGoodZvtxFT0vsPVCut{"cfgUseGoodZvtxFT0vsPVCut", false, "Apply kIsGoodZvtxFT0vsPV event selection"}; + Configurable cfgUseGoodITSLayersAll{"cfgUseGoodITSLayersAll", false, "Apply kIsGoodITSLayersAll 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"}; + Configurable cfgUseMCTruthCentrality{"cfgUseMCTruthCentrality", false, "Use MC truth centrality for generated Lambda(1520) histograms (otherwise the centrality of the accepted reconstructed collision is used when there is one)"}; + Configurable cfgApplyMCTruthINELgt0{"cfgApplyMCTruthINELgt0", true, "Apply truth-level INEL>0 selection to MC generated events (event factor, signal loss, generated QA)"}; + } evSel; + + struct : ConfigurableGroup { + Configurable cfgCutPt{"cfgCutPt", 0.15f, "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 the Lambda(1520) (reconstructed, generated and signal-loss particles)"}; + Configurable cfgMaxY{"cfgMaxY", 0.5f, "Maximum rapidity of the Lambda(1520) (reconstructed, generated and signal-loss particles)"}; + Configurable cfgMinTPCCrossedRows{"cfgMinTPCCrossedRows", 70, "Minimum TPC crossed pad rows"}; + Configurable cfgMinTPCClusters{"cfgMinTPCClusters", 120, "Minimum TPC clusters found"}; + Configurable cfgMinITSClusters{"cfgMinITSClusters", 5, "Minimum ITS clusters"}; + Configurable cfgMaxTPCSharedClusters{"cfgMaxTPCSharedClusters", 0, "Maximum TPC shared clusters"}; + Configurable cfgMinCrossedRowsOverFindable{"cfgMinCrossedRowsOverFindable", 0.8f, "Minimum ratio crossedRows/findableClusters"}; + Configurable cfgMaxTPCChi2NCl{"cfgMaxTPCChi2NCl", 4.0f, "Maximum TPC chi2/NCl"}; + Configurable cfgMaxITSChi2NCl{"cfgMaxITSChi2NCl", 36.0f, "Maximum ITS chi2/NCl"}; + Configurable applyTPCCrossedRowsCut{"applyTPCCrossedRowsCut", true, "Enable/disable minimum TPC crossed rows cut"}; + Configurable applyTPCFoundClusterCut{"applyTPCFoundClusterCut", false, "Enable/disable minimum TPC clusters found cut"}; + Configurable applyITSClusterCut{"applyITSClusterCut", false, "Enable/disable minimum ITS clusters cut"}; + Configurable applyTPCSharedClusterCut{"applyTPCSharedClusterCut", false, "Enable/disable maximum TPC shared clusters cut"}; + Configurable applyTPCCrossedRowsOverFindableCut{"applyTPCCrossedRowsOverFindableCut", false, "Enable/disable crossed rows / findable clusters cut"}; + Configurable applyTPCChi2NClCut{"applyTPCChi2NClCut", false, "Enable/disable maximum TPC chi2/NCl cut"}; + Configurable applyITSChi2NClCut{"applyITSChi2NClCut", false, "Enable/disable maximum ITS chi2/NCl cut"}; + Configurable requirePrimaryTrack{"requirePrimaryTrack", true, "Require isPrimaryTrack flag"}; + Configurable requireGlobalTrackWoDCA{"requireGlobalTrackWoDCA", true, "Require isGlobalTrackWoDCA flag"}; + Configurable requireGlobalTrack{"requireGlobalTrack", false, "Require isGlobalTrack flag"}; + Configurable requirePVContributor{"requirePVContributor", false, "Require PV-contributor flag"}; + } trackCuts; + + struct : ConfigurableGroup { + Configurable cfgCutDCAz{"cfgCutDCAz", 0.1f, "Maximum |DCAz| for all tracks [cm]"}; + Configurable> protonDCAPtBinEdges{"protonDCAPtBinEdges", {0.0f, 1000.f}, "Proton pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> protonMaxDCAxyPerPtBin{"protonMaxDCAxyPerPtBin", {0.1f}, "Max |DCAxy| for proton per pT bin [cm]"}; + Configurable> kaonDCAPtBinEdges{"kaonDCAPtBinEdges", {0.0f, 1000.f}, "Kaon pT bin edges for DCAxy cut [GeV/c]"}; + Configurable> kaonMaxDCAxyPerPtBin{"kaonMaxDCAxyPerPtBin", {0.1f}, "Max |DCAxy| for kaon per pT bin [cm]"}; + Configurable isApplyPtDepDCAxyCut{"isApplyPtDepDCAxyCut", false, "Additionally apply the continuous pT-dependent DCAxy cut"}; + Configurable dcaXYConst{"dcaXYConst", 0.004f, "Constant term of pT-dependent DCAxy cut [cm] (7-sigma: 0.004)"}; + Configurable dcaXYPtCoeff{"dcaXYPtCoeff", 0.013f, "pT-dependent coefficient of DCAxy cut [cm] (7-sigma: 0.013)"}; + Configurable dcaXYPtExponent{"dcaXYPtExponent", 1.0f, "pT exponent of DCAxy cut"}; + } dcaCuts; + + struct : ConfigurableGroup { + Configurable useTPCOnlyPIDAtAllMomenta{"useTPCOnlyPIDAtAllMomenta", false, "Bypass TOF: TPC-only PID (|nSigma| < max TPC nSigma) at all momenta"}; + Configurable applyKaonPionRejection{"applyKaonPionRejection", false, "Reject pion-like kaon candidates (no TOF, pT window, small pion nSigma)"}; + Configurable pionRejectionMinPt{"pionRejectionMinPt", 1.0f, "Pion rejection: minimum pT [GeV/c]"}; + Configurable pionRejectionMaxPt{"pionRejectionMaxPt", 2.0f, "Pion rejection: maximum pT [GeV/c]"}; + Configurable pionRejectionMaxTPCNSigma{"pionRejectionMaxTPCNSigma", 2.0f, "Pion rejection: |TPC nSigma pion| below this is rejected"}; + } pidShared; + + struct : ConfigurableGroup { + Configurable requireTOFForProton{"requireTOFForProton", false, "Require TOF signal for proton candidates"}; + Configurable maxTPCNSigmaProton{"maxTPCNSigmaProton", 2.0f, "Maximum |TPC nSigma| (used only with useTPCOnlyPIDAtAllMomenta)"}; + Configurable protonThresholdMomentum{"protonThresholdMomentum", 0.9f, "Proton momentum p [GeV/c] for tracks WITHOUT TOF: below it TPC-only PID, at/above it the track is rejected (tracks with TOF always use the TPC+TOF cut)"}; + Configurable> protonTPCPIDMomentumBins{"protonTPCPIDMomentumBins", {0.f, 0.9f}, "momentum edges p for TPC-only proton PID [GeV/c]"}; + Configurable> protonTPCNSigmaCutPerBin{"protonTPCNSigmaCutPerBin", {2.f}, "Max |TPC nSigma| per momentum bin, TPC-only proton PID"}; + Configurable> protonCombinedPIDMomentumBins{"protonCombinedPIDMomentumBins", {0.f, 999.f}, "momentum edges p for circular TPC+TOF proton PID [GeV/c]"}; + Configurable> protonCombinedNSigmaCutPerBin{"protonCombinedNSigmaCutPerBin", {2.f}, "Max combined nSigma radius per momentum bin, circular cut"}; + } protonPID; + + struct : ConfigurableGroup { + Configurable requireTOFForKaon{"requireTOFForKaon", false, "Require TOF signal for kaon candidates"}; + Configurable maxTPCNSigmaKaon{"maxTPCNSigmaKaon", 2.0f, "Maximum |TPC nSigma| (used only with useTPCOnlyPIDAtAllMomenta)"}; + Configurable kaonThresholdMomentum{"kaonThresholdMomentum", 0.5f, "Kaon momentum p [GeV/c] for tracks WITHOUT TOF: below it TPC-only PID, at/above it the track is rejected (tracks with TOF always use the TPC+TOF cut)"}; + Configurable> kaonTPCPIDMomentumBins{"kaonTPCPIDMomentumBins", {0.f, 0.5f}, "momentum edges p for TPC-only kaon PID [GeV/c]"}; + Configurable> kaonTPCNSigmaCutPerBin{"kaonTPCNSigmaCutPerBin", {2.f}, "Max |TPC nSigma| per momentum bin, TPC-only kaon PID"}; + Configurable> kaonCombinedPIDMomentumBins{"kaonCombinedPIDMomentumBins", {0.f, 999.f}, "momentum edges p for circular TPC+TOF kaon PID [GeV/c]"}; + Configurable> kaonCombinedNSigmaCutPerBin{"kaonCombinedNSigmaCutPerBin", {2.f}, "Max combined nSigma radius per momentum bin, circular cut"}; + } kaonPID; + + struct : ConfigurableGroup { + Configurable applyDaughterDeltaEtaPhiCuts{"applyDaughterDeltaEtaPhiCuts", false, "Apply |deltaEta| and |deltaPhi| cuts between daughters"}; + Configurable maxDaughterDeltaEta{"maxDaughterDeltaEta", 0.7f, "Maximum |deltaEta| between daughters"}; + Configurable maxDaughterDeltaPhi{"maxDaughterDeltaPhi", 1.5f, "Maximum |deltaPhi| between daughters"}; + Configurable applyDeepAngleCut{"applyDeepAngleCut", false, "Apply minimum opening-angle cut (removes split-track background)"}; + Configurable deepAngleCutValue{"deepAngleCutValue", 0.04f, "Minimum opening angle between proton and kaon [rad]"}; + Configurable applyOpeningAngleVeto{"applyOpeningAngleVeto", false, "Reject pairs whose opening angle lies INSIDE [veto min, veto max]"}; + Configurable openingAngleVetoMin{"openingAngleVetoMin", 1.4f, "Opening-angle veto window: lower edge [rad]"}; + Configurable openingAngleVetoMax{"openingAngleVetoMax", 2.4f, "Opening-angle veto window: upper edge [rad]"}; + Configurable applyPairPtMassUpperLimits{"applyPairPtMassUpperLimits", false, "Exclude pairs with pT or mass above the limits below (overflow bins)"}; + Configurable cfgMaxPairPt{"cfgMaxPairPt", 10.0f, "Maximum pair pT [GeV/c]"}; + Configurable cfgMaxPairMass{"cfgMaxPairMass", 3.0f, "Maximum pair mass [GeV/c^2]"}; + } pairCuts; + + struct : ConfigurableGroup { + Configurable cfgNoMixedEvents{"cfgNoMixedEvents", 10, "Number of mixed events per signal event"}; + Configurable cfgFillMixLikeSign{"cfgFillMixLikeSign", false, "Also fill like-sign pairs in mixed events"}; + } mixingCfg; + + struct : ConfigurableGroup { + Configurable enableRotationalBackground{"enableRotationalBackground", false, "Compute rotational background (unlike-sign pairs only)"}; + Configurable numberOfRotations{"numberOfRotations", 10, "Number of rotations per pair (evenly spaced, weight 1/N)"}; + Configurable rotationAngleWindow{"rotationAngleWindow", 10.f, "Rotation angles within PI +/- PI/rotationAngleWindow"}; + Configurable rotateProtonInsteadOfKaon{"rotateProtonInsteadOfKaon", true, "Rotate the proton phi (false = rotate the kaon phi)"}; + } rotBkg; + + struct : ConfigurableGroup { + Configurable cfgLambda1520PdgCode{"cfgLambda1520PdgCode", static_cast(o2::constants::physics::Pdg::kLambda1520_Py), "PDG code of the Lambda(1520) in the MC sample (change it for injected samples that use a different code). Sign is ignored: +code = Lambda(1520), -code = anti-Lambda(1520)"}; + Configurable applyGeneratedEtaCut{"applyGeneratedEtaCut", false, "Apply |eta| cut to generated particles and to truth-matched reconstructed pairs"}; + Configurable generatedEtaMax{"generatedEtaMax", 0.5f, "MC |eta| cut"}; + Configurable applyGeneratedDaughterEtaCut{"applyGeneratedDaughterEtaCut", true, "Require both generated daughters within |eta| < trackCuts.cfgCutEta"}; + Configurable applyMtScaleWeight{"applyMtScaleWeight", false, "Weight mT-scaled signal-loss histograms by mtScaledPt/pt"}; + Configurable useGeneratedMultiplicityCentrality{"useGeneratedMultiplicityCentrality", false, "Derive MC centrality from generated multiplicity"}; + Configurable> generatedMultiplicityLowerEdges{"generatedMultiplicityLowerEdges", std::vector{500, 300, 200, 120, 80, 50, 30, 10, 0}, "Generated multiplicity lower cuts"}; + Configurable> centralityPerGeneratedMultiplicityClass{"centralityPerGeneratedMultiplicityClass", std::vector{2.5f, 7.5f, 15.0f, 25.0f, 35.0f, 45.0f, 55.0f, 75.0f, 95.0f}, "Centrality bin centres"}; + Configurable generatedMultiplicityEtaMax{"generatedMultiplicityEtaMax", 0.5f, "Max eta for generated mid-rapidity multiplicity"}; + Configurable generatedMultiplicityEtaMin{"generatedMultiplicityEtaMin", -0.5f, "Min eta for generated mid-rapidity multiplicity"}; + } mcCfg; + + struct : ConfigurableGroup { + Configurable cfgFillTrackQA{"cfgFillTrackQA", true, "Fill track QA histograms (before/after selection)"}; + Configurable cfgFillPairQA{"cfgFillPairQA", false, "Fill pair QA histograms (opening angle, deltaEta, deltaPhi)"}; + Configurable cfgFillMultQA{"cfgFillMultQA", false, "Fill centrality vs track multiplicity correlation"}; + Configurable cfgFillMixQA{"cfgFillMixQA", false, "Fill mixed-event pool QA histograms"}; + } qa; + + struct : ConfigurableGroup { + Configurable referenceFrameDaughter{"referenceFrameDaughter", 0, "Reference daughter for cos(theta*): 0 = proton, 1 = kaon"}; + Configurable activateHelicityFrame{"activateHelicityFrame", false, "Compute cos(theta*) in the Helicity frame (quantization axis = Lambda(1520) flight direction)"}; + Configurable activateCollinsSoperFrame{"activateCollinsSoperFrame", false, "Compute cos(theta*) in the Collins-Soper frame (beam bisector)"}; + Configurable activateProductionFrame{"activateProductionFrame", false, "Compute cos(theta*) wrt the production-plane normal"}; + Configurable activateBeamAxisFrame{"activateBeamAxisFrame", false, "Compute cos(theta*) in the beam-axis (Gottfried-Jackson) frame"}; + Configurable activateRandomFrame{"activateRandomFrame", false, "Compute cos(theta*) wrt a random axis (null/control test, should be flat)"}; + Configurable fillRotatedBackgroundCosThetaStar{"fillRotatedBackgroundCosThetaStar", false, "Also fill rotated-background cos(theta*) histograms"}; + Configurable fillMixedEventCosThetaStar{"fillMixedEventCosThetaStar", false, "Also fill mixed-event-background cos(theta*) histograms"}; + Configurable beamEnergy{"beamEnergy", 13600.0f, "sqrt(s) beam energy in GeV (used for Collins-Soper/beam-axis/production frames)"}; + ConfigurableAxis cosThetaBins{"cosThetaBins", {20, -1.0f, 1.0f}, "cos(theta*) binning"}; + } spinAlignment; + + struct : ConfigurableGroup { + ConfigurableAxis cfgPtAxis{"cfgPtAxis", {VARIABLE_WIDTH, 0.0, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.25, 1.3, 1.4, 1.5, 1.6, 1.7, 1.75, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.1, 3.2, 3.3, 3.4, 3.6, 3.7, 3.8, 3.9, 4.0, 4.1, 4.2, 4.5, 4.6, 4.8, 4.9, 5.0, 5.5, 5.6, 6.0, 6.4, 6.5, 7.0, 7.2, 8.0, 9.0, 9.5, 9.6, 10.0, 11.0, 11.5, 12.0, 13.0, 14.0, 14.4, 15.0, 16.0, 18.0, 19.2, 20.0}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis cfgMassAxis{"cfgMassAxis", {70, 1.3f, 2.5f}, "Invariant mass (GeV/#it{c}^{2})"}; + ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.f, 5.f, 10.f, 20.f, 30.f, 40.f, 50.f, 60.f, 70.f, 80.f, 90.f, 100.f}, "Centrality (%)"}; + ConfigurableAxis cfgRapAxis{"cfgRapAxis", {20, -1.0, 1.0}, "Rapidity y"}; + ConfigurableAxis cfgVtxAxis{"cfgVtxAxis", {VARIABLE_WIDTH, -10.0f, -8.0f, -6.0f, -4.0f, -2.0f, 0.0f, 2.0f, 4.0f, 6.0f, 8.0f, 10.0f}, "Vertex z [cm] (histograms and event-mixing bins)"}; + } axes; + + struct DaughterTrackCandidate { + float px = 0.f, py = 0.f, pz = 0.f; + float pt = 0.f, eta = 0.f, phi = 0.f; + float mom = 0.f; + float dcaXY = 0.f, dcaZ = 0.f; + int sign = 0; + bool hasTOF = false; + float tpcNSigmaPr = 0.f, tpcNSigmaKa = 0.f; + float tofNSigmaPr = -999.f, tofNSigmaKa = -999.f; + int64_t globalIndex = 0; + bool hasMcLabel = false; + int mcPdgCode = 0; + int64_t mcLambda1520MotherIndex = -1; + int mcLambda1520MotherPdgCode = 0; + }; + + struct SpeciesPIDSettings { + std::vector tpcOnlyBins, tpcOnlyCuts; + std::vector combinedBins, combinedCuts; + float maxTPCNSigma = 2.f; + float thresholdMomentum = 0.f; + }; + + struct DecayKinematics { + std::array motherMom{}; + float motherE = 0.f; + std::array protonMom{}; + float protonE = 0.f; + std::array kaonMom{}; + float kaonE = 0.f; + }; + + struct AcceptedRecoCollision { + bool found = false; + float centrality = 0.f; + }; + + struct NoCutFlowTag { + }; + struct DataCutFlowTag { + }; + struct McRecoCutFlowTag { + }; + + enum class CosThetaStarSample { + Signal, + LikeSign, + Rotated, + Mixed, + LikeSignMixed, + TruthMatched + }; + + SpeciesPIDSettings mProtonPIDCuts; + SpeciesPIDSettings mKaonPIDCuts; + std::vector mProtonDCAPtEdges; + std::vector mProtonMaxDCAxy; + std::vector mKaonDCAPtEdges; + std::vector mKaonMaxDCAxy; + int mLambda1520AbsPdgCode = 0; + std::vector mProtonCandidates; + std::vector mKaonCandidates; + + float mBeam1Energy = 0.f; + std::array mBeam1Momentum{}; + float mBeam2Energy = 0.f; + std::array mBeam2Momentum{}; + std::mt19937 mRandomEngine{std::random_device{}()}; + std::uniform_real_distribution mUniform01{0., 1.}; + + template + static void setHistogramBinLabels(HistPtr const& hist, std::vector const& labels) + { + for (size_t i = 0; i < labels.size(); ++i) { + hist->GetXaxis()->SetBinLabel(static_cast(i) + 1, labels[i]); } + } - if (doprocessMCGen) { - allHistograms.add("SignalLoss/mcEventSelectionCutflow", "MC event selection cutflow (signal loss study)", kTH1F, {{7, 0, 7}}); - allHistograms.add("SignalLoss/mTScaled_fromProton", "mT-scaled Lambda(1520) pT from proton parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromAntiProton", "mT-scaled anti-Lambda(1520) pT from anti-proton parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromLambda0", "mT-scaled Lambda(1520) pT from Lambda0 parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromAntiLambda0", "mT-scaled anti-Lambda(1520) pT from anti-Lambda0 parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromXiMinus", "mT-scaled Lambda(1520) pT from Xi- parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromXiPlus", "mT-scaled anti-Lambda(1520) pT from Xi+ parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromXi0", "mT-scaled Lambda(1520) pT from Xi0 parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromAntiXi0", "mT-scaled anti-Lambda(1520) pT from anti-Xi0 parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromOmegaMinus", "mT-scaled Lambda(1520) pT from Omega- parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - allHistograms.add("SignalLoss/mTScaled_fromOmegaPlus", "mT-scaled anti-Lambda(1520) pT from Omega+ parent", kTHnSparseF, {axisPt, axisCentralityPercent}); - } - } - - // ============================================================ - // passesBasicTrackSelection() - // ============================================================ - template - bool passesBasicTrackSelection(TrackType const& track) + bool isSpinAlignmentActive() const { - if (track.pt() < minTrackPt) - return false; - if (track.eta() < minPseudorapidity || track.eta() > maxPseudorapidity) - return false; - if (requirePrimaryTrack && !track.isPrimaryTrack()) - return false; - if (requireGlobalTrackNoDCA && !track.isGlobalTrackWoDCA()) - return false; - if (requirePVContributor && !track.isPVContributor()) - return false; - if (applyCrossedRowsCut && track.tpcNClsCrossedRows() < minTPCCrossedRows) - return false; - if (applyTPCClustersCut && track.tpcNClsFound() < minTPCClustersFound) - return false; - return true; + return spinAlignment.activateHelicityFrame.value || + spinAlignment.activateCollinsSoperFrame.value || + spinAlignment.activateProductionFrame.value || + spinAlignment.activateBeamAxisFrame.value || + spinAlignment.activateRandomFrame.value; } - // ============================================================ - // passesProtonDCASelection() - // ============================================================ - template - bool passesProtonDCASelection(TrackType const& track, float trackPt) + void validateConfiguration() { - auto ptBinEdges = static_cast>(protonDCAPtBinEdges); - auto maxDCAxyValues = - static_cast>(protonMaxDCAxyPerPtBin); - int numberOfBins = static_cast(ptBinEdges.size()) - 1; + auto checkEdgesMatchCuts = [](auto const& cuts, auto const& bins, const char* label) { + const auto& cutsVec = cuts.value; + const auto& binsVec = bins.value; + if (cutsVec.size() != binsVec.size() - 1) { + LOG(fatal) << "PID bin-size mismatch for: " << label + << " (cuts=" << cutsVec.size() + << ", bins=" << binsVec.size() << ")"; + } + }; + checkEdgesMatchCuts(protonPID.protonTPCNSigmaCutPerBin, protonPID.protonTPCPIDMomentumBins, "proton TPC-only PID"); + checkEdgesMatchCuts(protonPID.protonCombinedNSigmaCutPerBin, protonPID.protonCombinedPIDMomentumBins, "proton combined PID"); + checkEdgesMatchCuts(kaonPID.kaonTPCNSigmaCutPerBin, kaonPID.kaonTPCPIDMomentumBins, "kaon TPC-only PID"); + checkEdgesMatchCuts(kaonPID.kaonCombinedNSigmaCutPerBin, kaonPID.kaonCombinedPIDMomentumBins, "kaon combined PID"); + checkEdgesMatchCuts(dcaCuts.protonMaxDCAxyPerPtBin, dcaCuts.protonDCAPtBinEdges, "proton DCAxy pT bins"); + checkEdgesMatchCuts(dcaCuts.kaonMaxDCAxyPerPtBin, dcaCuts.kaonDCAPtBinEdges, "kaon DCAxy pT bins"); + + if (mcCfg.useGeneratedMultiplicityCentrality && + mcCfg.generatedMultiplicityLowerEdges.value.size() != mcCfg.centralityPerGeneratedMultiplicityClass.value.size()) { + LOG(fatal) << "generatedMultiplicityLowerEdges and centralityPerGeneratedMultiplicityClass must have the same size"; + } - bool dcaXYPassed = false; - for (int iBin = 0; iBin < numberOfBins; iBin++) { - if (trackPt >= ptBinEdges[iBin] && - trackPt < ptBinEdges[iBin + 1] && - std::abs(track.dcaXY()) < maxDCAxyValues[iBin]) - dcaXYPassed = true; + const int nActiveFrames = static_cast(spinAlignment.activateHelicityFrame) + + static_cast(spinAlignment.activateCollinsSoperFrame) + + static_cast(spinAlignment.activateProductionFrame) + + static_cast(spinAlignment.activateBeamAxisFrame) + + static_cast(spinAlignment.activateRandomFrame); + if (nActiveFrames > 1) { + LOG(fatal) << "Multiple spin-alignment frames activated simultaneously. " + << "Please activate only ONE frame at a time to prevent data mixing in the cos(theta*) THnSparse."; + } + + if (doprocessMCGenSpinAlignment && !isSpinAlignmentActive()) { + LOG(fatal) << "processMCGenSpinAlignment requires at least one spin-alignment frame to be active"; } - if (!dcaXYPassed) - return false; - if (std::abs(track.dcaZ()) > maxDCAz) - return false; - return true; } - // ============================================================ - // passesKaonDCASelection() - // ============================================================ - template - bool passesKaonDCASelection(TrackType const& track, float trackPt) + void cacheCutSettings() { - auto ptBinEdges = static_cast>(kaonDCAPtBinEdges); - auto maxDCAxyValues = static_cast>(kaonMaxDCAxyPerPtBin); - int numberOfBins = static_cast(ptBinEdges.size()) - 1; + mLambda1520AbsPdgCode = std::abs(mcCfg.cfgLambda1520PdgCode.value); + if (mLambda1520AbsPdgCode == 0) { + LOG(fatal) << "cfgLambda1520PdgCode must be non-zero"; + } + if (doprocessMC || doprocessMCGen || doprocessSignalLoss) { + LOG(info) << "Lambda(1520) PDG code used for MC truth: " << mLambda1520AbsPdgCode; + } - bool dcaXYPassed = false; - for (int iBin = 0; iBin < numberOfBins; iBin++) { - if (trackPt >= ptBinEdges[iBin] && - trackPt < ptBinEdges[iBin + 1] && - std::abs(track.dcaXY()) < maxDCAxyValues[iBin]) - dcaXYPassed = true; + mProtonDCAPtEdges = dcaCuts.protonDCAPtBinEdges.value; + mProtonMaxDCAxy = dcaCuts.protonMaxDCAxyPerPtBin.value; + mKaonDCAPtEdges = dcaCuts.kaonDCAPtBinEdges.value; + mKaonMaxDCAxy = dcaCuts.kaonMaxDCAxyPerPtBin.value; + + mProtonPIDCuts.tpcOnlyBins = protonPID.protonTPCPIDMomentumBins.value; + mProtonPIDCuts.tpcOnlyCuts = protonPID.protonTPCNSigmaCutPerBin.value; + mProtonPIDCuts.combinedBins = protonPID.protonCombinedPIDMomentumBins.value; + mProtonPIDCuts.combinedCuts = protonPID.protonCombinedNSigmaCutPerBin.value; + mProtonPIDCuts.maxTPCNSigma = protonPID.maxTPCNSigmaProton; + mProtonPIDCuts.thresholdMomentum = protonPID.protonThresholdMomentum; + + mKaonPIDCuts.tpcOnlyBins = kaonPID.kaonTPCPIDMomentumBins.value; + mKaonPIDCuts.tpcOnlyCuts = kaonPID.kaonTPCNSigmaCutPerBin.value; + mKaonPIDCuts.combinedBins = kaonPID.kaonCombinedPIDMomentumBins.value; + mKaonPIDCuts.combinedCuts = kaonPID.kaonCombinedNSigmaCutPerBin.value; + mKaonPIDCuts.maxTPCNSigma = kaonPID.maxTPCNSigmaKaon; + mKaonPIDCuts.thresholdMomentum = kaonPID.kaonThresholdMomentum; + + auto checkTPCOnlyMomentumCoverage = [](std::vector const& edges, float threshold, const char* label) { + if (threshold > 0.f && (edges.empty() || edges.front() > 0.f || edges.back() < threshold)) { + LOG(warning) << label << ": TPC-only momentum bins do not cover [0, " << threshold + << ") GeV/c - tracks WITHOUT TOF in the uncovered range are rejected"; + } + }; + if (!pidShared.useTPCOnlyPIDAtAllMomenta) { + checkTPCOnlyMomentumCoverage(mProtonPIDCuts.tpcOnlyBins, mProtonPIDCuts.thresholdMomentum, "proton PID"); + checkTPCOnlyMomentumCoverage(mKaonPIDCuts.tpcOnlyBins, mKaonPIDCuts.thresholdMomentum, "kaon PID"); } - if (!dcaXYPassed) - return false; - if (std::abs(track.dcaZ()) > maxDCAz) - return false; - return true; } - // ============================================================ - // passesProtonPID() - // ============================================================ - template - bool passesProtonPID(TrackType const& track, float totalMomentum) + void computeBeamKinematics() { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - - auto tpcMomBins = static_cast>(protonTPCPIDMomentumBins); - auto tpcNSigCuts = - static_cast>(protonTPCNSigmaCutPerBin); - auto tofMomBins = static_cast>(protonTOFPIDMomentumBins); - auto tofNSigCuts = - static_cast>(protonTOFNSigmaCutPerBin); - int nTPCBins = static_cast(tpcMomBins.size()); - int nTOFBins = static_cast(tofMomBins.size()); - - float tpcNSigPion = std::abs(track.tpcNSigmaPi()); - float tpcNSigKaon = std::abs(track.tpcNSigmaKa()); - float tpcNSigProton = std::abs(track.tpcNSigmaPr()); - float tofNSigPion = std::abs(track.tofNSigmaPi()); - float tofNSigKaon = std::abs(track.tofNSigmaKa()); - float tofNSigProton = std::abs(track.tofNSigmaPr()); - - float combinedNSigPion = - tpcNSigPion * tpcNSigPion + tofNSigPion * tofNSigPion; - float combinedNSigKaon = - tpcNSigKaon * tpcNSigKaon + tofNSigKaon * tofNSigKaon; - float combinedNSigProton = - tpcNSigProton * tpcNSigProton + tofNSigProton * tofNSigProton; - - float circularCutSquared = - combinedNSigmaCutProton * combinedNSigmaCutProton; - float circularRejCutSquared = - tofNSigmaVetoThreshold * tpcNSigmaVetoThreshold; - - if (!useTPCOnlyPID && track.hasTOF()) { - if (track.tofNSigmaPr() < minTOFNSigmaProton) - return false; + const auto sqrtS = static_cast(spinAlignment.beamEnergy); + const double beamMomentum = std::sqrt(sqrtS * sqrtS / 4. - static_cast(ProtonMass) * static_cast(ProtonMass)); + mBeam1Energy = static_cast(sqrtS * 0.5); + mBeam1Momentum = {0.f, 0.f, static_cast(-beamMomentum)}; + mBeam2Energy = static_cast(sqrtS * 0.5); + mBeam2Momentum = {0.f, 0.f, static_cast(beamMomentum)}; + } - if (combinedNSigmaCutProton < 0) { - for (int i = 0; i < nTOFBins - 1; ++i) { - if (totalMomentum >= tofMomBins[i] && - totalMomentum < tofMomBins[i + 1] && - (tofNSigProton < tofNSigCuts[i] && - tofNSigPion > tofNSigmaVetoThreshold && - tofNSigKaon > tofNSigmaVetoThreshold)) - tofPIDPassed = true; - } - if (track.tpcNSigmaPr() < minCombinedTPCTOFNSigmaProton) - return false; - if (tpcNSigProton < maxTPCNSigmaProton && - tpcNSigPion > tpcNSigmaVetoThreshold && - tpcNSigKaon > tpcNSigmaVetoThreshold) - tpcPIDPassed = true; + void init(InitContext const&) + { + validateConfiguration(); + cacheCutSettings(); + computeBeamKinematics(); + + const AxisSpec ptAxis{axes.cfgPtAxis, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec massAxis{axes.cfgMassAxis, "M_{inv} (GeV/#it{c}^{2})"}; + const AxisSpec centAxis{axes.cfgCentAxis, "Centrality (%)"}; + const AxisSpec rapAxis{axes.cfgRapAxis, "Rapidity y"}; + const AxisSpec vtxAxis{axes.cfgVtxAxis, "Vertex z [cm]"}; + const AxisSpec cosThetaAxis{spinAlignment.cosThetaBins, "cos#theta*"}; + const AxisSpec nSigmaTPCaxis{100, -5., 5., "n#sigma^{TPC}"}; + const AxisSpec nSigmaTOFaxis{100, -5., 5., "n#sigma^{TOF}"}; + const AxisSpec momentumAxis{100, 0., 10., "p (GeV/#it{c})"}; + const AxisSpec ptForPIDAxis{100, 0., 10., "p_{T} (GeV/#it{c})"}; + const AxisSpec dcaXYaxis{200, -0.2f, 0.2f, "DCA_{xy} (cm)"}; + const AxisSpec dcaZaxis{200, -1.2f, 1.2f, "DCA_{z} (cm)"}; + const AxisSpec tpcClusAxis{200, 0, 200, "N_{clusters}^{TPC}"}; + const AxisSpec tpcRowsAxis{200, 0, 200, "Crossed rows^{TPC}"}; + const AxisSpec tpcSignalAxis{500, 0, 500, "TPC signal (a.u.)"}; + const AxisSpec occupancyAxis{200, 0, 10000, "Track occupancy [-40, 100] ns"}; + const AxisSpec etaAxis{150, -1.5, 1.5, "#eta"}; + const AxisSpec phiAxis{72, 0., o2::constants::math::TwoPI, "#varphi (rad)"}; + const AxisSpec openAngleAxis{100, 0., o2::constants::math::PI, "Opening angle (rad)"}; + const AxisSpec deltaEtaAxis{160, 0., 1.6, "|#Delta#eta|"}; + const AxisSpec deltaPhiAxis{100, 0., o2::constants::math::PI, "|#Delta#varphi| (rad)"}; + const AxisSpec massResolutionAxis{200, -0.05f, 0.05f, "M_{reco} - M_{gen} (GeV/c^{2})"}; + const std::vector eventCutLabels{"All", "Vz", "kNoITSROFrameBorder", "sel8", "INEL>0", "Occupancy", "NoSameBunchPileup", "GoodZvtxFT0vsPV", "GoodITSLayersAll", "Accepted"}; + + if (doprocessSameEvent) { + histos.add("CutFlow/Data/hEventSelection", "Data event cut flow", kTH1F, {{10, -0.5f, 9.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/Data/hEventSelection")), eventCutLabels); + + histos.add("EventQA/Data/hVertexZ", "Vertex z; z (cm)", kTH1F, {{400, -20., 20.}}); + histos.add("EventQA/Data/hPVContributors", "PV contributors; N", kTH1F, {{2001, -0.5f, 2000.5f}}); + histos.add("EventQA/Data/hCentrality", "Centrality", kTH1F, {centAxis}); + histos.add("EventQA/Data/hOccupancy", "Occupancy in time range", kTH1F, {occupancyAxis}); + histos.add("EventQA/Data/hCentralityVsVertexZ", "Centrality vs vertex z", kTH2F, {vtxAxis, centAxis}); + histos.add("EventQA/Data/hCentralityVsOccupancy", "Centrality vs occupancy", kTH2F, {occupancyAxis, centAxis}); + + if (qa.cfgFillMultQA) { + histos.add("EventQA/Data/hCentralityVsFilteredTracksVsPVTracks", "Centrality vs filtered tracks vs PV tracks", kTHnSparseF, {{110, 0, 110, "Centrality"}, {500, 0, 5000, "Filtered tracks"}, {500, 0, 5000, "PV tracks"}}); } - if ((combinedNSigmaCutProton > 0) && - (combinedNSigProton < circularCutSquared && - combinedNSigPion > circularRejCutSquared && - combinedNSigKaon > circularRejCutSquared)) { - tofPIDPassed = true; - tpcPIDPassed = true; + if (qa.cfgFillTrackQA) { + histos.add("TrackQA/Data/PreSelection/Proton/hNSigmaTPCvsP", "TPC n#sigma proton vs p (before PID)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/PreSelection/Proton/hNSigmaTOFvsP", "TOF n#sigma proton vs p (before PID)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/PreSelection/Proton/hNSigmaTOFvsNSigmaTPC", "TOF vs TPC n#sigma proton (before PID)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/PreSelection/Kaon/hNSigmaTPCvsP", "TPC n#sigma kaon vs p (before PID)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/PreSelection/Kaon/hNSigmaTOFvsP", "TOF n#sigma kaon vs p (before PID)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/PreSelection/Kaon/hNSigmaTOFvsNSigmaTPC", "TOF vs TPC n#sigma kaon (before PID)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + + histos.add("TrackQA/Data/AllTracks/hEta", "Reco dN/d#eta; #eta; Tracks", kTH1F, {{50, -1.0, 1.0}}); + histos.add("TrackQA/Data/AllTracks/hPt", "Reco pT; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/Data/AllTracks/hPhi", "Reco #varphi; #varphi (rad); Tracks", kTH1F, {phiAxis}); + histos.add("TrackQA/Data/AllTracks/hEtaPhi", "Reco #eta vs #varphi; #eta; #varphi", kTH2F, {etaAxis, phiAxis}); + + histos.add("TrackQA/Data/Selected/Proton/hDCAxyVsPt", "Proton DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/Data/Selected/Proton/hDCAzVsPt", "Proton DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTPCvsP", "Proton TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTPCvsPt", "Proton TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTPCvsCentrality", "Proton TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTOFvsP", "Proton TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTOFvsPt", "Proton TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTOFvsCentrality", "Proton TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTOFvsNSigmaTPC", "Proton TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTPCPionVsPt", "Proton track: TPC n#sigma pion contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Proton/hNSigmaTPCKaonVsPt", "Proton track: TPC n#sigma kaon contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Proton/hTPCSignalVsInnerP", "Proton TPC signal vs p (after cuts)", kTH2F, {momentumAxis, tpcSignalAxis}); + histos.add("TrackQA/Data/Selected/Proton/hTPCClustersVsPhi", "Proton TPC clusters vs #varphi", kTH2F, {tpcClusAxis, phiAxis}); + histos.add("TrackQA/Data/Selected/Proton/hTPCCrossedRowsVsPt", "Proton TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("TrackQA/Data/Selected/Proton/hTPCFoundClustersVsPt", "Proton TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("TrackQA/Data/Selected/Proton/hPt", "Proton p_{T}; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/Data/Selected/Proton/hEta", "Proton #eta; #eta; Tracks", kTH1F, {etaAxis}); + + histos.add("TrackQA/Data/Selected/Kaon/hDCAxyVsPt", "Kaon DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hDCAzVsPt", "Kaon DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsP", "Kaon TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsPt", "Kaon TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsCentrality", "Kaon TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsP", "Kaon TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsPt", "Kaon TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsCentrality", "Kaon TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsNSigmaTPC", "Kaon TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTPCPionVsPt", "Kaon track: TPC n#sigma pion contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hNSigmaTPCProtonVsPt", "Kaon track: TPC n#sigma proton contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/Data/Selected/Kaon/hTPCSignalVsInnerP", "Kaon TPC signal vs p (after cuts)", kTH2F, {momentumAxis, tpcSignalAxis}); + histos.add("TrackQA/Data/Selected/Kaon/hTPCClustersVsPhi", "Kaon TPC clusters vs #varphi", kTH2F, {tpcClusAxis, phiAxis}); + histos.add("TrackQA/Data/Selected/Kaon/hTPCCrossedRowsVsPt", "Kaon TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("TrackQA/Data/Selected/Kaon/hTPCFoundClustersVsPt", "Kaon TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("TrackQA/Data/Selected/Kaon/hPt", "Kaon p_{T}; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/Data/Selected/Kaon/hEta", "Kaon #eta; #eta; Tracks", kTH1F, {etaAxis}); } - } else { - if (!track.hasTOF() && totalMomentum >= minProtonMomentum) - return false; - tofPIDPassed = true; - if (track.tpcNSigmaPr() < minTPCNSigmaProton) - return false; - for (int i = 0; i < nTPCBins - 1; ++i) { - if (totalMomentum >= tpcMomBins[i] && - totalMomentum < tpcMomBins[i + 1] && - (tpcNSigProton < tpcNSigCuts[i] && - tpcNSigPion > tpcNSigmaVetoPion && - tpcNSigKaon > tpcNSigmaVetoKaonForProton)) - tpcPIDPassed = true; + if (qa.cfgFillPairQA) { + histos.add("PairQA/Data/hOpeningAngleBeforeCuts", "Opening angle (before pair cuts)", kTH1F, {openAngleAxis}); + histos.add("PairQA/Data/hOpeningAngleAfterCuts", "Opening angle (after pair cuts)", kTH1F, {openAngleAxis}); + histos.add("PairQA/Data/hDeltaEtaBeforeCuts", "|#Delta#eta| (before pair cuts)", kTH1F, {deltaEtaAxis}); + histos.add("PairQA/Data/hDeltaEtaAfterCuts", "|#Delta#eta| (after pair cuts)", kTH1F, {deltaEtaAxis}); + histos.add("PairQA/Data/hDeltaPhiBeforeCuts", "|#Delta#varphi| (before pair cuts)", kTH1F, {deltaPhiAxis}); + histos.add("PairQA/Data/hDeltaPhiAfterCuts", "|#Delta#varphi| (after pair cuts)", kTH1F, {deltaPhiAxis}); } - } - return (tpcPIDPassed && tofPIDPassed); - } - - // ============================================================ - // passesKaonPID() - // ============================================================ - template - bool passesKaonPID(TrackType const& track, float totalMomentum) - { - bool tpcPIDPassed{false}, tofPIDPassed{false}; - - auto tpcMomBins = static_cast>(kaonTPCPIDMomentumBins); - auto tpcNSigCuts = static_cast>(kaonTPCNSigmaCutPerBin); - auto tofMomBins = static_cast>(kaonTOFPIDMomentumBins); - auto tofNSigCuts = static_cast>(kaonTOFNSigmaCutPerBin); - int nTPCBins = static_cast(tpcMomBins.size()); - int nTOFBins = static_cast(tofMomBins.size()); - - float tpcNSigPion = std::abs(track.tpcNSigmaPi()); - float tpcNSigKaon = std::abs(track.tpcNSigmaKa()); - float tpcNSigProton = std::abs(track.tpcNSigmaPr()); - float tofNSigPion = std::abs(track.tofNSigmaPi()); - float tofNSigKaon = std::abs(track.tofNSigmaKa()); - float tofNSigProton = std::abs(track.tofNSigmaPr()); - - float combinedNSigPion = - tpcNSigPion * tpcNSigPion + tofNSigPion * tofNSigPion; - float combinedNSigKaon = - tpcNSigKaon * tpcNSigKaon + tofNSigKaon * tofNSigKaon; - float combinedNSigProton = - tpcNSigProton * tpcNSigProton + tofNSigProton * tofNSigProton; - float circularCutSquared = combinedNSigmaCutKaon * combinedNSigmaCutKaon; - float circularRejCutSquared = - tpcNSigmaVetoOtherSpecies * tofNSigmaVetoThreshold; - - if (!useTPCOnlyPID && track.hasTOF()) { - if (track.tofNSigmaKa() < minTOFNSigmaKaon) - return false; - - if (combinedNSigmaCutKaon < 0) { - for (int i = 0; i < nTOFBins - 1; ++i) { - if (totalMomentum >= tofMomBins[i] && - totalMomentum < tofMomBins[i + 1] && - (tofNSigKaon < tofNSigCuts[i] && - tofNSigPion > tofNSigmaVetoThreshold && - tofNSigProton > tofNSigmaVetoThreshold)) - tofPIDPassed = true; + if (!isSpinAlignmentActive()) { + histos.add("InvMass/Data/SameEvent/hLambda1520", "THnSparse #Lambda(1520) same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/SameEvent/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) same-event", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/SameEvent/hLikeSignProtonPlusKaonPlus", "THnSparse like-sign p K^{+}", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/SameEvent/hLikeSignProtonMinusKaonMinus", "THnSparse like-sign #bar{p} K^{-}", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + if (rotBkg.enableRotationalBackground) { + histos.add("InvMass/Data/RotatedBackground/hLambda1520", "THnSparse #Lambda(1520) rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/RotatedBackground/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) rotational background", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.get(HIST("InvMass/Data/RotatedBackground/hLambda1520"))->Sumw2(); + histos.get(HIST("InvMass/Data/RotatedBackground/hAntiLambda1520"))->Sumw2(); } - if (track.tpcNSigmaKa() < minCombinedTPCTOFNSigmaKaon) - return false; - if (tpcNSigKaon < maxTPCNSigmaKaon && - tpcNSigPion > tpcNSigmaVetoThreshold && - tpcNSigProton > tpcNSigmaVetoThreshold) - tpcPIDPassed = true; } - if ((combinedNSigmaCutKaon > 0) && - (combinedNSigKaon < circularCutSquared && - combinedNSigPion > circularRejCutSquared && - combinedNSigProton > circularRejCutSquared)) { - tofPIDPassed = true; - tpcPIDPassed = true; - } + if (isSpinAlignmentActive()) { + histos.add("CosThetaStar/Data/SameEvent/hLambda1520", "cos#theta* #Lambda(1520) unlike-sign", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/SameEvent/hAntiLambda1520", "cos#theta* #bar{#Lambda}(1520) unlike-sign", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/SameEvent/hLikeSignProtonPlusKaonPlus", "cos#theta* like-sign p K^{+}", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/SameEvent/hLikeSignProtonMinusKaonMinus", "cos#theta* like-sign #bar{p} K^{-}", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); - } else { - if (!track.hasTOF() && totalMomentum >= minKaonMomentum) - return false; - tofPIDPassed = true; - if (track.tpcNSigmaKa() < minTPCNSigmaKaon) - return false; - for (int i = 0; i < nTPCBins - 1; ++i) { - if (totalMomentum >= tpcMomBins[i] && - totalMomentum < tpcMomBins[i + 1] && - (tpcNSigKaon < tpcNSigCuts[i] && - tpcNSigPion > tpcNSigmaVetoPionForKaon && - tpcNSigProton > tpcNSigmaVetoProtonForKaon)) - tpcPIDPassed = true; + if (spinAlignment.fillRotatedBackgroundCosThetaStar && rotBkg.enableRotationalBackground) { + histos.add("CosThetaStar/Data/RotatedBackground/hLambda1520", "cos#theta* rotated bkg #Lambda(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/RotatedBackground/hAntiLambda1520", "cos#theta* rotated bkg #bar{#Lambda}(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + } } } - return (tpcPIDPassed && tofPIDPassed); - } + if (doprocessMixedEvent) { + if (!isSpinAlignmentActive()) { + histos.add("InvMass/Data/MixedEvent/hLambda1520", "THnSparse #Lambda(1520) event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/MixedEvent/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); - template - void fillTrackQAHistograms(TrackType const& track, float centralityPercent) - { - if (!passesBasicTrackSelection(track)) - return; + if (mixingCfg.cfgFillMixLikeSign) { + histos.add("InvMass/Data/MixedEvent/hLikeSignProtonPlusKaonPlus", "THnSparse like-sign p K^{+} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/Data/MixedEvent/hLikeSignProtonMinusKaonMinus", "THnSparse like-sign #bar{p} K^{-} event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + } - const float trackPt = track.pt(); - const float totalMomentum = RecoDecay::p(track.px(), track.py(), track.pz()); + if (qa.cfgFillMixQA) { + histos.add("MixingQA/Data/hPoolVertexZ", "Mixed event pool: vertex z; z (cm); Pairs", kTH1F, {vtxAxis}); + histos.add("MixingQA/Data/hPoolCentrality", "Mixed event pool: centrality; Centrality (%); Pairs", kTH1F, {centAxis}); + histos.add("MixingQA/Data/hPoolVertexZVsCentrality", "Mixed event pool: vertex z vs centrality", kTH2F, {vtxAxis, centAxis}); + histos.add("MixingQA/Data/hCollisionPairs", "Number of mixed-event collision pairs processed", kTH1F, {{1, 0.5f, 1.5f}}); + } + } - const auto tpcNSigProtonRaw = track.tpcNSigmaPr(); - allHistograms.fill(HIST("QAbefore/Proton/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigProtonRaw); - const auto tpcNSigKaonRaw = track.tpcNSigmaKa(); - allHistograms.fill(HIST("QAbefore/Kaon/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigKaonRaw); - if (track.hasTOF()) { - const auto tofNSigProtonRaw = track.tofNSigmaPr(); - allHistograms.fill(HIST("QAbefore/Proton/tofNSigmaVsMomentum"), totalMomentum, tofNSigProtonRaw); - allHistograms.fill(HIST("QAbefore/Proton/tofNSigmaVsTPCNSigma"), tpcNSigProtonRaw, tofNSigProtonRaw); - const auto tofNSigKaonRaw = track.tofNSigmaKa(); - allHistograms.fill(HIST("QAbefore/Kaon/tofNSigmaVsMomentum"), totalMomentum, tofNSigKaonRaw); - allHistograms.fill(HIST("QAbefore/Kaon/tofNSigmaVsTPCNSigma"), tpcNSigKaonRaw, tofNSigKaonRaw); - } - - const bool isProtonCandidate = - (!requireTOFForProton || track.hasTOF()) && - passesProtonPID(track, totalMomentum) && - passesProtonDCASelection(track, trackPt); - - if (isProtonCandidate) { - allHistograms.fill(HIST("QAafter/Proton/ptVsCentrality"), trackPt, centralityPercent); - allHistograms.fill(HIST("QAafter/Proton/dcaZVsPt"), trackPt, track.dcaZ()); - allHistograms.fill(HIST("QAafter/Proton/dcaXYVsPt"), trackPt, track.dcaXY()); - allHistograms.fill(HIST("QAafter/Proton/tpcDedxVsMomentum"), totalMomentum, track.tpcSignal()); - allHistograms.fill(HIST("QAafter/Proton/tpcNSigmaPionContamVsPt"), totalMomentum, track.tpcNSigmaPi()); - allHistograms.fill(HIST("QAafter/Proton/tpcNSigmaKaonContamVsPt"), totalMomentum, track.tpcNSigmaKa()); - allHistograms.fill(HIST("QAafter/Proton/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigProtonRaw); - allHistograms.fill(HIST("QAafter/Proton/tpcNSigmaVsCentrality"), centralityPercent, tpcNSigProtonRaw); - allHistograms.fill(HIST("QAafter/Proton/tpcNSigmaVsPt"), trackPt, tpcNSigProtonRaw); - allHistograms.fill(HIST("QAafter/Proton/tpcCrossedRowsVsPt"), trackPt, track.tpcNClsCrossedRows()); - allHistograms.fill(HIST("QAafter/Proton/tpcClustersFoundVsPt"), trackPt, track.tpcNClsFound()); - - if (!useTPCOnlyPID && track.hasTOF()) { - const auto tofNSigProton = track.tofNSigmaPr(); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaVsMomentum"), totalMomentum, tofNSigProton); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaVsCentrality"), centralityPercent, tofNSigProton); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaVsPt"), trackPt, tofNSigProton); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaPionContamVsMomentum"), totalMomentum, track.tofNSigmaPi()); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaKaonContamVsMomentum"), totalMomentum, track.tofNSigmaKa()); - allHistograms.fill(HIST("QAafter/Proton/tofNSigmaVsTPCNSigma"), tpcNSigProtonRaw, tofNSigProton); - } - } - - const bool isKaonCandidate = - (!requireTOFForKaon || track.hasTOF()) && - passesKaonPID(track, totalMomentum) && - passesKaonDCASelection(track, trackPt); - - if (isKaonCandidate) { - allHistograms.fill(HIST("QAafter/Kaon/ptVsCentrality"), trackPt, centralityPercent); - allHistograms.fill(HIST("QAafter/Kaon/dcaZVsPt"), trackPt, track.dcaZ()); - allHistograms.fill(HIST("QAafter/Kaon/dcaXYVsPt"), trackPt, track.dcaXY()); - allHistograms.fill(HIST("QAafter/Kaon/tpcDedxVsMomentum"), totalMomentum, track.tpcSignal()); - allHistograms.fill(HIST("QAafter/Kaon/tpcNSigmaPionContamVsPt"), totalMomentum, track.tpcNSigmaPi()); - allHistograms.fill(HIST("QAafter/Kaon/tpcNSigmaProtonContamVsMomentum"), totalMomentum, track.tpcNSigmaPr()); - allHistograms.fill(HIST("QAafter/Kaon/tpcNSigmaVsMomentum"), totalMomentum, tpcNSigKaonRaw); - allHistograms.fill(HIST("QAafter/Kaon/tpcNSigmaVsCentrality"), centralityPercent, tpcNSigKaonRaw); - allHistograms.fill(HIST("QAafter/Kaon/tpcNSigmaVsPt"), trackPt, tpcNSigKaonRaw); - allHistograms.fill(HIST("QAafter/Kaon/tpcCrossedRowsVsPt"), trackPt, track.tpcNClsCrossedRows()); - allHistograms.fill(HIST("QAafter/Kaon/tpcClustersFoundVsPt"), trackPt, track.tpcNClsFound()); - - if (!useTPCOnlyPID && track.hasTOF()) { - const auto tofNSigKaon = track.tofNSigmaKa(); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaVsMomentum"), totalMomentum, tofNSigKaon); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaVsPt"), trackPt, tofNSigKaon); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaVsCentrality"), centralityPercent, tofNSigKaon); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaPionContamVsMomentum"), totalMomentum, track.tofNSigmaPi()); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaProtonContamVsMomentum"), totalMomentum, track.tofNSigmaPr()); - allHistograms.fill(HIST("QAafter/Kaon/tofNSigmaVsTPCNSigma"), tpcNSigKaonRaw, tofNSigKaon); - } - } - } - - // ============================================================ - // fillInvariantMassHistograms() - // ============================================================ - template - void fillInvariantMassHistograms(TrackCollectionType const& protonCandidates, - TrackCollectionType const& kaonCandidates, - float centralityPercent, - int occupancyValue = 100) - { - float protonTotalMomentum = 0., kaonTotalMomentum = 0.; - - for (auto const& [protonTrack, kaonTrack] : - soa::combinations(soa::CombinationsFullIndexPolicy(protonCandidates, - kaonCandidates))) { - if (protonTrack.index() == kaonTrack.index()) - continue; - if (!passesBasicTrackSelection(protonTrack) || - !passesBasicTrackSelection(kaonTrack)) - continue; - auto pxProton = protonTrack.px(); - auto pyProton = protonTrack.py(); - auto pzProton = protonTrack.pz(); - auto pxKaon = kaonTrack.px(); - auto pyKaon = kaonTrack.py(); - auto pzKaon = kaonTrack.pz(); - // FIX root/entity: replaced TMath::Sqrt with RecoDecay::p - protonTotalMomentum = RecoDecay::p(pxProton, pyProton, pzProton); - kaonTotalMomentum = RecoDecay::p(pxKaon, pyKaon, pzKaon); - - if (requireTOFForProton && !protonTrack.hasTOF()) - continue; - if (requireTOFForKaon && !kaonTrack.hasTOF()) - continue; + if (isSpinAlignmentActive() && spinAlignment.fillMixedEventCosThetaStar) { + histos.add("CosThetaStar/Data/MixedEvent/hLambda1520", "cos#theta* mixed-event #Lambda(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/MixedEvent/hAntiLambda1520", "cos#theta* mixed-event #bar{#Lambda}(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/MixedEvent/hLikeSignProtonPlusKaonPlus", "cos#theta* mixed-event like-sign p K^{+}", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/Data/MixedEvent/hLikeSignProtonMinusKaonMinus", "cos#theta* mixed-event like-sign #bar{p} K^{-}", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + } + } - if (!passesProtonPID(protonTrack, protonTotalMomentum) || - !passesKaonPID(kaonTrack, kaonTotalMomentum)) - continue; + if (doprocessMC) { + histos.add("CutFlow/MCReco/hEventSelection", "MC reconstructed event cut flow", kTH1F, {{10, -0.5f, 9.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/MCReco/hEventSelection")), eventCutLabels); + + histos.add("EventQA/MCReco/hAcceptedEventsVsCentrality", "Reconstructed INEL>0 events (Nrec, MC reco)", kTH1F, {centAxis}); + histos.add("EventQA/MCReco/hVertexZ", "MC reco vertex z; z (cm)", kTH1F, {{400, -20., 20.}}); + histos.add("EventQA/MCReco/hPVContributors", "MC reco PV contributors; N", kTH1F, {{2001, -0.5f, 2000.5f}}); + + histos.add("InvMass/MCReco/AllPairs/hLambda1520", "THnSparse #Lambda(1520) MC reconstructed (pre-truth-matching)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/MCReco/AllPairs/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) MC reconstructed (pre-truth-matching)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + histos.add("InvMass/MCReco/TruthMatched/hLambda1520", "THnSparse #Lambda(1520) MC reconstructed (truth-matched)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/MCReco/TruthMatched/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) MC reconstructed (truth-matched)", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + histos.add("MCRecoQA/DaughterPt/hProtonFromLambda1520", "Rec proton from #Lambda(1520)", kTH1F, {ptAxis}); + histos.add("MCRecoQA/DaughterPt/hProtonFromAntiLambda1520", "Rec proton from #bar{#Lambda}(1520)", kTH1F, {ptAxis}); + histos.add("MCRecoQA/DaughterPt/hKaonFromLambda1520", "Rec kaon from #Lambda(1520)", kTH1F, {ptAxis}); + histos.add("MCRecoQA/DaughterPt/hKaonFromAntiLambda1520", "Rec kaon from #bar{#Lambda}(1520)", kTH1F, {ptAxis}); + + histos.add("TrackQA/MCReco/AllTracks/hEta", "MC Reco dN/d#eta; #eta; Tracks", kTH1F, {{50, -1.0, 1.0}}); + histos.add("TrackQA/MCReco/AllTracks/hPt", "MC Reco pT; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/MCReco/AllTracks/hPhi", "MC Reco #varphi; #varphi (rad); Tracks", kTH1F, {phiAxis}); + histos.add("TrackQA/MCReco/AllTracks/hEtaPhi", "MC Reco #eta vs #varphi; #eta; #varphi", kTH2F, {etaAxis, phiAxis}); + histos.add("PairQA/MCReco/hTruthMatchedRapidityVsPt", "Truth-matched reconstructed #Lambda(1520): y vs p_{T}", kTH2F, {ptAxis, rapAxis}); + + if (qa.cfgFillTrackQA) { + histos.add("TrackQA/MCReco/PreSelection/Proton/hNSigmaTPCvsP", "MC reco: TPC n#sigma proton vs p (before PID)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/PreSelection/Proton/hNSigmaTOFvsP", "MC reco: TOF n#sigma proton vs p (before PID)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/PreSelection/Proton/hNSigmaTOFvsNSigmaTPC", "MC reco: TOF vs TPC n#sigma proton (before PID)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTPCvsP", "MC reco: TPC n#sigma kaon vs p (before PID)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTOFvsP", "MC reco: TOF n#sigma kaon vs p (before PID)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTOFvsNSigmaTPC", "MC reco: TOF vs TPC n#sigma kaon (before PID)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + + histos.add("TrackQA/MCReco/Selected/Proton/hDCAxyVsPt", "MC reco: Proton DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hDCAzVsPt", "MC reco: Proton DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsP", "MC reco: Proton TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsPt", "MC reco: Proton TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsCentrality", "MC reco: Proton TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsP", "MC reco: Proton TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsPt", "MC reco: Proton TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsCentrality", "MC reco: Proton TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsNSigmaTPC", "MC reco: Proton TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTPCPionVsPt", "MC reco: Proton track: TPC n#sigma pion contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hNSigmaTPCKaonVsPt", "MC reco: Proton track: TPC n#sigma kaon contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hTPCSignalVsInnerP", "MC reco: Proton TPC signal vs p (after cuts)", kTH2F, {momentumAxis, tpcSignalAxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hTPCClustersVsPhi", "MC reco: Proton TPC clusters vs #varphi", kTH2F, {tpcClusAxis, phiAxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hTPCCrossedRowsVsPt", "MC reco: Proton TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hTPCFoundClustersVsPt", "MC reco: Proton TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hPt", "MC reco: Proton p_{T}; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/MCReco/Selected/Proton/hEta", "MC reco: Proton #eta; #eta; Tracks", kTH1F, {etaAxis}); + + histos.add("TrackQA/MCReco/Selected/Kaon/hDCAxyVsPt", "MC reco: Kaon DCA_{xy} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hDCAzVsPt", "MC reco: Kaon DCA_{z} vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsP", "MC reco: Kaon TPC n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsPt", "MC reco: Kaon TPC n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsCentrality", "MC reco: Kaon TPC n#sigma vs centrality", kTH2F, {centAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsP", "MC reco: Kaon TOF n#sigma vs p (after cuts)", kTH2F, {momentumAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsPt", "MC reco: Kaon TOF n#sigma vs p_{T} (after cuts)", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsCentrality", "MC reco: Kaon TOF n#sigma vs centrality", kTH2F, {centAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsNSigmaTPC", "MC reco: Kaon TOF vs TPC n#sigma (after cuts)", kTH2F, {nSigmaTPCaxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCPionVsPt", "MC reco: Kaon track: TPC n#sigma pion contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCProtonVsPt", "MC reco: Kaon track: TPC n#sigma proton contamination", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hTPCSignalVsInnerP", "MC reco: Kaon TPC signal vs p (after cuts)", kTH2F, {momentumAxis, tpcSignalAxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hTPCClustersVsPhi", "MC reco: Kaon TPC clusters vs #varphi", kTH2F, {tpcClusAxis, phiAxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hTPCCrossedRowsVsPt", "MC reco: Kaon TPC crossed rows vs p_{T}", kTH2F, {ptForPIDAxis, tpcRowsAxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hTPCFoundClustersVsPt", "MC reco: Kaon TPC clusters found vs p_{T}", kTH2F, {ptForPIDAxis, tpcClusAxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hPt", "MC reco: Kaon p_{T}; p_{T} (GeV/c); Tracks", kTH1F, {ptForPIDAxis}); + histos.add("TrackQA/MCReco/Selected/Kaon/hEta", "MC reco: Kaon #eta; #eta; Tracks", kTH1F, {etaAxis}); + } - if (!passesProtonDCASelection(protonTrack, protonTrack.pt()) || - !passesKaonDCASelection(kaonTrack, kaonTrack.pt())) - continue; + histos.add("TrackQA/MCReco/TruthMatched/Proton/hDCAxyVsPt", "Truth-matched proton DCA_{xy} vs p_{T}", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Proton/hDCAzVsPt", "Truth-matched proton DCA_{z} vs p_{T}", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Proton/hNSigmaTPCvsPt", "Truth-matched proton TPC n#sigma vs p_{T}", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Proton/hNSigmaTOFvsPt", "Truth-matched proton TOF n#sigma vs p_{T}", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Kaon/hDCAxyVsPt", "Truth-matched kaon DCA_{xy} vs p_{T}", kTH2F, {ptForPIDAxis, dcaXYaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Kaon/hDCAzVsPt", "Truth-matched kaon DCA_{z} vs p_{T}", kTH2F, {ptForPIDAxis, dcaZaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Kaon/hNSigmaTPCvsPt", "Truth-matched kaon TPC n#sigma vs p_{T}", kTH2F, {ptForPIDAxis, nSigmaTPCaxis}); + histos.add("TrackQA/MCReco/TruthMatched/Kaon/hNSigmaTOFvsPt", "Truth-matched kaon TOF n#sigma vs p_{T}", kTH2F, {ptForPIDAxis, nSigmaTOFaxis}); + + histos.add("CutFlow/MCReco/hGeneratedToReconstructed", "Generated-to-reconstructed #Lambda(1520) cut flow (efficiency QA)", kTH1F, {{19, -0.5f, 18.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/MCReco/hGeneratedToReconstructed")), + {"Generated Lambda(1520) (event accepted)", "Lambda(1520) -> p + K", "Rapidity accepted", + "Proton acceptance", "Kaon acceptance", "Both daughters acceptance", + "Proton reconstructed", "Kaon reconstructed", "Both daughters reconstructed", + "Proton track quality", "Kaon track quality", "Both track quality", + "Proton PID", "Kaon PID", "Both PID", + "Proton DCA", "Kaon DCA", "Both DCA", + "Pair cuts (final reconstructed)"}); + + histos.add("CutFlow/MCReco/hTruthMatching", "MC truth-matching cut flow (efficiency QA)", kTH1F, {{6, -0.5f, 5.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/MCReco/hTruthMatching")), + {"Reco pair", "Proton has MC label", "Kaon has MC label", "Proton is proton", "Kaon is kaon", "Common Lambda(1520) mother"}); + + histos.add("InvMass/MCReco/MassResolution/hLambda1520", "(M_{reco} - M_{gen}) truth-matched #Lambda(1520)", kTHnSparseF, {massResolutionAxis, ptAxis, centAxis}); + histos.add("InvMass/MCReco/MassResolution/hAntiLambda1520", "(M_{reco} - M_{gen}) truth-matched #bar{#Lambda}(1520)", kTHnSparseF, {massResolutionAxis, ptAxis, centAxis}); + + if (isSpinAlignmentActive()) { + histos.add("CosThetaStar/MCReco/TruthMatched/hLambda1520", "cos#theta* truth-matched #Lambda(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCReco/TruthMatched/hAntiLambda1520", "cos#theta* truth-matched #bar{#Lambda}(1520)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCReco/MassResolution/hLambda1520", "(M_{reco} - M_{gen}) vs cos#theta* #Lambda(1520)", kTHnSparseF, {massResolutionAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCReco/MassResolution/hAntiLambda1520", "(M_{reco} - M_{gen}) vs cos#theta* #bar{#Lambda}(1520)", kTHnSparseF, {massResolutionAxis, ptAxis, centAxis, cosThetaAxis}); + } + } - // FIX root/entity: replaced TMath::ACos with std::acos - if (applyDeepAngleCut && - std::acos((protonTrack.pt() * kaonTrack.pt() + pzProton * pzKaon) / - (protonTotalMomentum * kaonTotalMomentum)) < - deepAngleCutValue) - continue; + if (doprocessMCEventMixing) { + histos.add("InvMass/MCReco/MixedEvent/hLambda1520", "THnSparse #Lambda(1520) MC reconstructed event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/MCReco/MixedEvent/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) MC reconstructed event mixing", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("MixingQA/MCReco/hCollisionPairs", "Number of MC mixed-event collision pairs processed", kTH1F, {{1, 0.5f, 1.5f}}); + histos.add("MixingQA/MCReco/hCollisionPairsVsCentrality", "MC mixed-event collision pairs vs centrality", kTH1F, {centAxis}); + } - if (runQualityChecksOnly) - continue; + if (doprocessMCGen) { + histos.add("CutFlow/MCGen/hEventSelection", "MC generated event cut flow", kTH1F, {{6, -0.5f, 5.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/MCGen/hEventSelection")), + {"All", "Vz", "INEL>0", "Centrality", "Final generated event", "Associated reco event accepted"}); + + histos.add("InvMass/MCGen/Generated/hLambda1520", "THnSparse #Lambda(1520) generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/MCGen/Generated/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) generated", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + histos.add("InvMass/MCGen/GeneratedInAcceptedEvents/hLambda1520", "THnSparse #Lambda(1520) generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + histos.add("InvMass/MCGen/GeneratedInAcceptedEvents/hAntiLambda1520", "THnSparse #bar{#Lambda}(1520) generated, event-accepted", kTHnSparseF, {massAxis, ptAxis, centAxis, rapAxis}); + + histos.add("MCGenQA/DaughterPt/hProtonFromLambda1520", "Gen proton from #Lambda(1520)", kTH1F, {ptAxis}); + histos.add("MCGenQA/DaughterPt/hProtonFromAntiLambda1520", "Gen proton from #bar{#Lambda}(1520)", kTH1F, {ptAxis}); + histos.add("MCGenQA/DaughterPt/hKaonFromLambda1520", "Gen kaon from #Lambda(1520)", kTH1F, {ptAxis}); + histos.add("MCGenQA/DaughterPt/hKaonFromAntiLambda1520", "Gen kaon from #bar{#Lambda}(1520)", kTH1F, {ptAxis}); + + histos.add("CutFlow/MCGen/hGeneratedLambda1520Selection", "Generated #Lambda(1520) cut flow", kTH1F, {{8, -0.5f, 7.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection")), + {"All particles", "producedByGenerator", "Lambda(1520) PDG", "Decay p + K", "Rapidity", "Daughter eta", "Eta", "Final"}); + + histos.add("MCGenQA/Kinematics/hEtaVsPtBeforeCuts", "Generated #Lambda(1520) #eta-p_{T} (before kinematic cuts)", kTH2F, {etaAxis, ptAxis}); + histos.add("MCGenQA/Kinematics/hPhiVsRapidityBeforeCuts", "Generated #Lambda(1520) #varphi-y (before kinematic cuts)", kTH2F, {phiAxis, rapAxis}); + histos.add("MCGenQA/Kinematics/hEtaVsPtAfterCuts", "Generated #Lambda(1520) #eta-p_{T} (after kinematic cuts)", kTH2F, {etaAxis, ptAxis}); + histos.add("MCGenQA/Kinematics/hPhiVsRapidityAfterCuts", "Generated #Lambda(1520) #varphi-y (after kinematic cuts)", kTH2F, {phiAxis, rapAxis}); + histos.add("MCGenQA/Kinematics/hRapidityVsPtAccepted", "Generated #Lambda(1520) y vs p_{T} (final)", kTH2F, {ptAxis, rapAxis}); + } - std::array momProton = {pxProton, pyProton, pzProton}; - std::array momKaon = {pxKaon, pyKaon, pzKaon}; - std::array, 2> bothMomenta = {momProton, momKaon}; + if (doprocessMCGenSpinAlignment) { + histos.add("CosThetaStar/MCGen/Generated/hLambda1520", "cos#theta* generated #Lambda(1520) (spin-alignment process)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCGen/Generated/hAntiLambda1520", "cos#theta* generated #bar{#Lambda}(1520) (spin-alignment process)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCGen/GeneratedInAcceptedEvents/hLambda1520", "cos#theta* generated #Lambda(1520), event-accepted (spin-alignment process)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + histos.add("CosThetaStar/MCGen/GeneratedInAcceptedEvents/hAntiLambda1520", "cos#theta* generated #bar{#Lambda}(1520), event-accepted (spin-alignment process)", kTHnSparseF, {massAxis, ptAxis, centAxis, cosThetaAxis}); + } - float pairInvMass = RecoDecay::m(bothMomenta, std::array{MassProton, MassKaonCharged}); - float pairPt = RecoDecay::pt(std::array{pxProton + pxKaon, pyProton + pyKaon}); - float pairRapidity = RecoDecay::y(std::array{pxProton + pxKaon, pyProton + pyKaon, pzProton + pzKaon}, pairInvMass); + if (doprocessEventFactor) { + histos.add("EventLoss/hAcceptedRecoEventsVsTruthCentrality", "Accepted reconstructed events of generated INEL>0 events (Nrec, MC-truth centrality)", kTH1F, {centAxis}); + histos.add("EventLoss/hGeneratedEventsVsTruthCentrality", "Generated events passing truth |Zvtx| + truth INEL>0 cuts", kTH1F, {centAxis}); + histos.add("EventLoss/hGeneratedAcceptedEventsVsTruthCentrality", "Generated events with >=1 accepted reconstructed collision (truth centrality)", kTH1F, {centAxis}); + histos.add("EventLoss/hGeneratedAcceptedEventsVsRecoCentrality", "Generated events with >=1 accepted reconstructed collision (reconstructed centrality)", kTH1F, {centAxis}); + histos.add("EventLoss/hRecoCollisionMultiplicity", "Number of reconstructed collisions per generated collision", kTH1F, {{21, -0.5f, 20.5f}}); + histos.add("CutFlow/EventLoss/hGeneratedEventSelection", "Event Loss / Event Splitting: generated-event cut flow", kTH1F, {{6, -0.5f, 5.5f}}); + setHistogramBinLabels(histos.get(HIST("CutFlow/EventLoss/hGeneratedEventSelection")), + {"All", "Vz", "INEL>0", "HasRecoColl", "RecoCollAccepted", "EventAccepted"}); + } - if (pairRapidity < minPairRapidity || pairRapidity > maxPairRapidity) - continue; + if (doprocessSignalLoss) { + histos.add("SignalLoss/Lambda1520/hGeneratedPtVsCentrality", "Generated #Lambda(1520), generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/AntiLambda1520/hGeneratedPtVsCentrality", "Generated #bar{#Lambda}(1520), generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Lambda1520/hAcceptedEventPtVsCentrality", "Generated #Lambda(1520), accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/AntiLambda1520/hAcceptedEventPtVsCentrality", "Generated #bar{#Lambda}(1520), accepted events", kTH2F, {ptAxis, centAxis}); + + histos.add("SignalLoss/Kaon/hGeneratedTruePtVsCentrality", "Kaon, true pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Kaon/hGeneratedMtScaledPtVsCentrality", "Kaon, mT-scaled pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Kaon/hAcceptedEventTruePtVsCentrality", "Kaon, true pT, accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Kaon/hAcceptedEventMtScaledPtVsCentrality", "Kaon, mT-scaled pT, accepted events", kTH2F, {ptAxis, centAxis}); + + histos.add("SignalLoss/Proton/hGeneratedTruePtVsCentrality", "Proton, true pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Proton/hGeneratedMtScaledPtVsCentrality", "Proton, mT-scaled pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Proton/hAcceptedEventTruePtVsCentrality", "Proton, true pT, accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Proton/hAcceptedEventMtScaledPtVsCentrality", "Proton, mT-scaled pT, accepted events", kTH2F, {ptAxis, centAxis}); + + histos.add("SignalLoss/Lambda0/hGeneratedTruePtVsCentrality", "Lambda, true pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Lambda0/hGeneratedMtScaledPtVsCentrality", "Lambda, mT-scaled pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Lambda0/hAcceptedEventTruePtVsCentrality", "Lambda, true pT, accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Lambda0/hAcceptedEventMtScaledPtVsCentrality", "Lambda, mT-scaled pT, accepted events", kTH2F, {ptAxis, centAxis}); + + histos.add("SignalLoss/Xi/hGeneratedTruePtVsCentrality", "Xi, true pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Xi/hGeneratedMtScaledPtVsCentrality", "Xi, mT-scaled pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Xi/hAcceptedEventTruePtVsCentrality", "Xi, true pT, accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Xi/hAcceptedEventMtScaledPtVsCentrality", "Xi, mT-scaled pT, accepted events", kTH2F, {ptAxis, centAxis}); + + histos.add("SignalLoss/Omega/hGeneratedTruePtVsCentrality", "Omega, true pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Omega/hGeneratedMtScaledPtVsCentrality", "Omega, mT-scaled pT, generated", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Omega/hAcceptedEventTruePtVsCentrality", "Omega, true pT, accepted events", kTH2F, {ptAxis, centAxis}); + histos.add("SignalLoss/Omega/hAcceptedEventMtScaledPtVsCentrality", "Omega, mT-scaled pT, accepted events", kTH2F, {ptAxis, centAxis}); + } + } - if constexpr (!isMixedEvent && !isMCAnalysis) { - if (protonTrack.sign() * kaonTrack.sign() < 0) { - if (protonTrack.sign() > 0) - allHistograms.fill(HIST("Analysis/invMass_UnlikeSign_ProtonPlusKaonMinus"), pairInvMass, pairPt, centralityPercent, occupancyValue); - else - allHistograms.fill(HIST("Analysis/invMass_UnlikeSign_ProtonMinusKaonPlus"), pairInvMass, pairPt, centralityPercent, occupancyValue); - - if (enableRotationalBackground) { - for (int iRot = 0; iRot < numberOfRotations; iRot++) { - // FIX magic-number: named constant for the half-window - // calculation - float rotationWindowHalfWidth = - o2::constants::math::PI / rotationAngleWindow; - float rotatedKaonPhi; - if (numberOfRotations == 1) { - rotatedKaonPhi = o2::constants::math::PI; - } else { - rotatedKaonPhi = - (o2::constants::math::PI - rotationWindowHalfWidth) + - iRot * (2.f * rotationWindowHalfWidth / - (numberOfRotations - 1)); - } + template + float getCentrality(CollisionType const& collision) + { + switch (evSel.cfgCentralityEstimator) { + case lambda1520analysispo_analysis::kFT0M: + return collision.centFT0M(); + case lambda1520analysispo_analysis::kFT0A: + return collision.centFT0A(); + case lambda1520analysispo_analysis::kFT0C: + return collision.centFT0C(); + case lambda1520analysispo_analysis::kFV0A: + return collision.centFV0A(); + default: + return collision.centFT0M(); + } + } - // FIX two-pi-add-subtract: use RecoDecay::constrainAngle instead - // of manual +/- TwoPI - float newKaonPhi = RecoDecay::constrainAngle( - kaonTrack.phi() + rotatedKaonPhi, 0.f); - - float pxKaonRotated = kaonTrack.pt() * std::cos(newKaonPhi); - float pyKaonRotated = kaonTrack.pt() * std::sin(newKaonPhi); - - std::array momProtonRot = {pxProton, pyProton, - pzProton}; - std::array momKaonRot = {pxKaonRotated, pyKaonRotated, - pzKaon}; - std::array, 2> rotatedMomenta = { - momProtonRot, momKaonRot}; - - float rotatedPairInvMass = RecoDecay::m( - rotatedMomenta, std::array{MassProton, MassKaonCharged}); - float rotatedPairPt = RecoDecay::pt(std::array{ - pxProton + pxKaonRotated, pyProton + pyKaonRotated}); - float rotatedPairRapidity = RecoDecay::y( - std::array{pxProton + pxKaonRotated, pyProton + pyKaonRotated, - pzProton + pzKaon}, - MassLambda1520); - - if (rotatedPairRapidity < minPairRapidity || - rotatedPairRapidity > maxPairRapidity) - continue; - - if (protonTrack.sign() * kaonTrack.sign() < 0) { - if (protonTrack.sign() > 0) - allHistograms.fill( - HIST("Analysis/invMass_Rotated_ProtonPlusKaonMinus"), - rotatedPairInvMass, rotatedPairPt, centralityPercent, - occupancyValue); - else - allHistograms.fill( - HIST("Analysis/invMass_Rotated_ProtonMinusKaonPlus"), - rotatedPairInvMass, rotatedPairPt, centralityPercent, - occupancyValue); - } - } - } + template + bool isRecoInelGt0(CollisionType const& collision) + { + return collision.isInelGt0(); + } - } else { - if (protonTrack.sign() > 0) - allHistograms.fill( - HIST("Analysis/invMass_LikeSign_ProtonPlusKaonPlus"), - pairInvMass, pairPt, centralityPercent, occupancyValue); - else - allHistograms.fill( - HIST("Analysis/invMass_LikeSign_ProtonMinusKaonMinus"), - pairInvMass, pairPt, centralityPercent, occupancyValue); - } - } + template + bool isGeneratedChargedPrimary(McParticleType const& mcPart) + { + if (!mcPart.isPhysicalPrimary()) { + return false; + } + auto* pdgParticle = pdgDB->GetParticle(mcPart.pdgCode()); + if (!pdgParticle) { + return false; + } + return std::abs(pdgParticle->Charge()) > MinAbsCharge; + } - if constexpr (isMixedEvent) { - if (protonTrack.sign() * kaonTrack.sign() < 0) { - if (protonTrack.sign() > 0) - allHistograms.fill( - HIST("Analysis/invMass_Mixed_ProtonPlusKaonMinus"), pairInvMass, - pairPt, centralityPercent, occupancyValue); - else - allHistograms.fill( - HIST("Analysis/invMass_Mixed_ProtonMinusKaonPlus"), pairInvMass, - pairPt, centralityPercent, occupancyValue); - } else { - if (protonTrack.sign() > 0) - allHistograms.fill(HIST("Analysis/invMass_Mixed_LikeSign_PlusPlus"), - pairInvMass, pairPt, centralityPercent, - occupancyValue); - else - allHistograms.fill( - HIST("Analysis/invMass_Mixed_LikeSign_MinusMinus"), pairInvMass, - pairPt, centralityPercent, occupancyValue); - } + template + bool isTruthInelGt0(McParticleTableType const& mcParticlesInCollision) + { + for (auto const& mcPart : mcParticlesInCollision) { + if (std::abs(mcPart.eta()) < 1.f && isGeneratedChargedPrimary(mcPart)) { + return true; } + } + return false; + } - if constexpr (isMCAnalysis) { - if (protonTrack.sign() * kaonTrack.sign() < 0) { - // FIX pdg/explicit-code: use named constants from lambda_analysis - // namespace - if (std::abs(protonTrack.pdgCode()) != lambda_analysis::PdgProton || - std::abs(kaonTrack.pdgCode()) != lambda_analysis::PdgKaon) - continue; - if (protonTrack.motherId() != kaonTrack.motherId()) - continue; - if (protonTrack.motherPDG() != kaonTrack.motherPDG()) - continue; - if (protonTrack.pdgCode() == 0 || kaonTrack.pdgCode() == 0) - continue; - if (protonTrack.motherPDG() == -1 || kaonTrack.motherPDG() == -1) - continue; - if (std::abs(protonTrack.motherPDG()) != pdgCodeLambda1520) - continue; - - float trueMassOfParent = 0.; - for (auto const& parentParticle : *mcResonanceParentTable) { - if (parentParticle.mcParticleId() == protonTrack.motherId()) { - std::array parentMom = {parentParticle.px(), - parentParticle.py(), - parentParticle.pz()}; - trueMassOfParent = RecoDecay::m(parentMom, parentParticle.e()); - break; - } - } - - float massResidual = pairInvMass - trueMassOfParent; - - if (protonTrack.motherPDG() > 0) { - allHistograms.fill(HIST("Analysis/mcReconstructed_Lambda1520"), - pairInvMass, pairPt, centralityPercent); - allHistograms.fill(HIST("Analysis/mcMassResolution_Lambda1520"), - massResidual, pairPt, centralityPercent); - } else { - allHistograms.fill(HIST("Analysis/mcReconstructed_AntiLambda1520"), - pairInvMass, pairPt, centralityPercent); - allHistograms.fill(HIST("Analysis/mcMassResolution_AntiLambda1520"), - massResidual, pairPt, centralityPercent); - } - } + template + int countGeneratedChargedPrimaries(McPartsT const& partsThisMc) + { + int nChargedPrimaries = 0; + for (auto const& part : partsThisMc) { + if (!part.isPhysicalPrimary()) { + continue; + } + if (part.eta() > mcCfg.generatedMultiplicityEtaMax || part.eta() < mcCfg.generatedMultiplicityEtaMin) { + continue; + } + auto* pdgParticle = pdgDB->GetParticle(part.pdgCode()); + if (!pdgParticle || std::abs(pdgParticle->Charge()) < MinAbsCharge) { + continue; } + ++nChargedPrimaries; } + return nChargedPrimaries; } - // ── Type aliases ───────────────────────────────────────────────────────── - // FIX name/type: use UpperCamelCase for type aliases - using ResonanceCollisionsWithEP = - soa::Join; - using ResonanceMCCollisions = - soa::Join; - using ResonanceTrackTable = aod::ResoTracks; - - // ============================================================ - // processData() - // ============================================================ - void processData(ResonanceCollisionsWithEP::iterator const& collision, - ResonanceTrackTable const& tracks) + float getCentralityFromGeneratedMultiplicity(int nChargedPrimaries) { - if (cEvtRecINELgt0 && !collision.isRecINELgt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - return; + const auto& cuts = mcCfg.generatedMultiplicityLowerEdges.value; + const auto& centres = mcCfg.centralityPerGeneratedMultiplicityClass.value; + for (size_t i = 0; i < cuts.size(); ++i) { + if (nChargedPrimaries >= cuts[i]) { + return centres[i]; + } + } + return InvalidCentrality; + } - allHistograms.fill(HIST("Event/centralityVsOccupancy"), collision.cent(), - 100); - allHistograms.fill(HIST("Event/primaryVertexZ"), collision.posZ()); - for (const auto& track : tracks) { - fillTrackQAHistograms(track, collision.cent()); + template + float getTruthCentrality(McCollisionType const& mcCollision, McPartsT const& partsThisMc) + { + if (mcCfg.useGeneratedMultiplicityCentrality) { + return getCentralityFromGeneratedMultiplicity(countGeneratedChargedPrimaries(partsThisMc)); } - fillInvariantMassHistograms(tracks, tracks, collision.cent()); + return mcCollision.centFT0M(); } - PROCESS_SWITCH(lambda1520analysispo, processData, - "Process real collision data (same-event analysis)", true); - // ============================================================ - // processMC() - // FIX const-ref-in-process: table subscriptions use const& - // ============================================================ - void processMC(ResonanceMCCollisions::iterator const& collision, - soa::Join const& tracks, - aod::ResoMCParents const& mcParents) + template + void fillEventSelectionCutFlowBin(float bin) { - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 0); + if constexpr (std::is_same_v) { + histos.fill(HIST("CutFlow/Data/hEventSelection"), bin); + } else if constexpr (std::is_same_v) { + histos.fill(HIST("CutFlow/MCReco/hEventSelection"), bin); + } + } - if (mcRequireTriggerTVX && !collision.isTriggerTVX()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 1); + template + bool passesEventSelectionImpl(CollisionType const& collision) + { + fillEventSelectionCutFlowBin(0.f); - if (mcRequireVtxWithin10cm && !collision.isVtxIn10()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 2); + if (std::abs(collision.posZ()) > evSel.cfgCutVertex) { + return false; + } + fillEventSelectionCutFlowBin(1.f); - if (mcRequireINELgt0 && !collision.isINELgt0()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 3); + if (evSel.cfgUseNoITSROFrameBorderCut && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + return false; + } + fillEventSelectionCutFlowBin(2.f); - if (mcRequireSel8 && !collision.isInSel8()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 4); + if (evSel.cfgUseSel8 && !collision.sel8()) { + return false; + } + fillEventSelectionCutFlowBin(3.f); - if (cEvtRecINELgt0 && !collision.isRecINELgt0()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 5); + if (evSel.cfgRequireRecoINELgt0 && !isRecoInelGt0(collision)) { + return false; + } + fillEventSelectionCutFlowBin(4.f); - if (mcRequireAfterAllCuts && !collision.isInAfterAllCuts()) - return; - allHistograms.fill(HIST("Event/mcEventSelectionCutflow"), 6); - - auto centralityPercent = collision.cent(); - allHistograms.fill(HIST("Event/centralityVsOccupancy"), centralityPercent, - 100); - allHistograms.fill(HIST("Event/primaryVertexZ"), collision.posZ()); - - mcResonanceParentTable = &mcParents; - fillInvariantMassHistograms(tracks, tracks, centralityPercent); - - // FIX const-ref-in-for-loop: use const auto& in all range-based for loops - for (const auto& track : tracks) { - allHistograms.fill(HIST("QAbefore/trackEta"), track.eta()); - allHistograms.fill(HIST("QAbefore/trackPt"), track.pt()); - allHistograms.fill(HIST("QAbefore/trackPhi"), track.phi()); - allHistograms.fill(HIST("QAbefore/trackEtaVsPhi"), track.eta(), - track.phi()); - - // FIX pdg/explicit-code: use named PDG constants - if (std::abs(track.pdgCode()) == lambda_analysis::PdgKaon) - allHistograms.fill(HIST("QAChecks/kaonGenLevelPt"), track.pt()); - if (std::abs(track.pdgCode()) == lambda_analysis::PdgProton) - allHistograms.fill(HIST("QAChecks/protonGenLevelPt"), track.pt()); - - if (!passesBasicTrackSelection(track)) - continue; + if (evSel.applyOccupancyInTimeRangeCut) { + const int trackOccupancy = collision.trackOccupancyInTimeRange(); + if (trackOccupancy < evSel.cfgOccupancyMin || trackOccupancy > evSel.cfgOccupancyMax) { + return false; + } + } + fillEventSelectionCutFlowBin(5.f); - fillTrackQAHistograms(track, centralityPercent); + if (evSel.cfgUseNoSameBunchPileupCut && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return false; + } + fillEventSelectionCutFlowBin(6.f); - // FIX root/entity: replaced TMath::Sqrt with RecoDecay::p - float totalMom = RecoDecay::p(track.px(), track.py(), track.pz()); + if (evSel.cfgUseGoodZvtxFT0vsPVCut && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + fillEventSelectionCutFlowBin(7.f); - if (passesKaonPID(track, totalMom) && - std::abs(track.pdgCode()) == lambda_analysis::PdgKaon) - allHistograms.fill(HIST("QAChecks/kaonRecoLevelPt"), track.pt()); - if (passesProtonPID(track, totalMom) && - std::abs(track.pdgCode()) == lambda_analysis::PdgProton) - allHistograms.fill(HIST("QAChecks/protonRecoLevelPt"), track.pt()); + if (evSel.cfgUseGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) { + return false; } + fillEventSelectionCutFlowBin(8.f); - // FIX const-ref-in-for-loop: use const auto& in all range-based for loops - for (const auto& parentParticle : mcParents) { - if (std::abs(parentParticle.pdgCode()) != pdgCodeLambda1520) - continue; - if (parentParticle.y() < minPairRapidity || - parentParticle.y() > maxPairRapidity) - continue; + fillEventSelectionCutFlowBin(9.f); + return true; + } - // FIX pdg/explicit-code: use named PDG constants - bool hasProtonDaughter = (std::abs(parentParticle.daughterPDG1()) == - lambda_analysis::PdgProton || - std::abs(parentParticle.daughterPDG2()) == - lambda_analysis::PdgProton); - bool hasKaonDaughter = - (std::abs(parentParticle.daughterPDG1()) == - lambda_analysis::PdgKaon || - std::abs(parentParticle.daughterPDG2()) == lambda_analysis::PdgKaon); - if (!hasProtonDaughter || !hasKaonDaughter) - continue; + template + bool passesEventSelection(CollisionType const& collision) + { + return passesEventSelectionImpl(collision); + } + + template + AcceptedRecoCollision findFirstAcceptedRecoCollision(CollisionGroup const& collisions) + { + for (auto const& collision : collisions) { + if (passesEventSelection(collision)) { + return {true, getCentrality(collision)}; + } + } + return {}; + } + + template + bool passesBasicTrackSelection(TrackType const& track) + { + if (trackCuts.applyTPCCrossedRowsCut && + track.tpcNClsCrossedRows() < trackCuts.cfgMinTPCCrossedRows) { + return false; + } + if (trackCuts.applyTPCFoundClusterCut && + track.tpcNClsFound() < trackCuts.cfgMinTPCClusters) { + return false; + } + if (trackCuts.applyITSClusterCut && + track.itsNCls() < trackCuts.cfgMinITSClusters) { + return false; + } + if (trackCuts.applyTPCSharedClusterCut && + track.tpcNClsShared() > trackCuts.cfgMaxTPCSharedClusters) { + return false; + } + if (trackCuts.applyTPCCrossedRowsOverFindableCut && + track.tpcNClsCrossedRows() < trackCuts.cfgMinCrossedRowsOverFindable * track.tpcNClsFindable()) { + return false; + } + if (trackCuts.applyTPCChi2NClCut && + track.tpcChi2NCl() > trackCuts.cfgMaxTPCChi2NCl) { + return false; + } + if (trackCuts.applyITSChi2NClCut && + track.itsChi2NCl() > trackCuts.cfgMaxITSChi2NCl) { + return false; + } + if (trackCuts.requirePrimaryTrack && !track.isPrimaryTrack()) { + return false; + } + if (trackCuts.requireGlobalTrackWoDCA && !track.isGlobalTrackWoDCA()) { + return false; + } + if (trackCuts.requireGlobalTrack && !track.isGlobalTrack()) { + return false; + } + if (trackCuts.requirePVContributor && !track.isPVContributor()) { + return false; + } + return true; + } - std::array parentMom = { - parentParticle.px(), parentParticle.py(), parentParticle.pz()}; - float parentMass = RecoDecay::m(parentMom, parentParticle.e()); + static int findBinIndex(std::vector const& edges, float x) + { + for (size_t i = 0; i + 1 < edges.size(); ++i) { + if (x >= edges[i] && x < edges[i + 1]) { + return static_cast(i); + } + } + return -1; + } - if (parentParticle.pdgCode() > 0) - allHistograms.fill(HIST("Analysis/mcGenerated_Lambda1520"), parentMass, - parentParticle.pt(), centralityPercent); - else - allHistograms.fill(HIST("Analysis/mcGenerated_AntiLambda1520"), - parentMass, parentParticle.pt(), centralityPercent); + template + bool passesDCAxyDCAzSelection(TrackType const& track, std::vector const& ptEdges, std::vector const& maxDCAxyPerBin) + { + if (std::abs(track.dcaZ()) >= dcaCuts.cfgCutDCAz) { + return false; + } + const float pt = track.pt(); + const float dcaXY = std::abs(track.dcaXY()); + const int bin = findBinIndex(ptEdges, pt); + if (bin < 0 || dcaXY >= maxDCAxyPerBin[bin]) { + return false; + } + if (dcaCuts.isApplyPtDepDCAxyCut) { + const float dcaConst = dcaCuts.dcaXYConst; + const float dcaCoeff = dcaCuts.dcaXYPtCoeff; + const float dcaExp = dcaCuts.dcaXYPtExponent; + const float maxPtDepDCAxy = dcaConst + dcaCoeff / std::pow(pt, dcaExp); + if (dcaXY > maxPtDepDCAxy) { + return false; + } } + return true; } - PROCESS_SWITCH(lambda1520analysispo, processMC, - "Process Monte Carlo simulated events", false); - // ============================================================ - // processMCGen() - // ============================================================ - void processMCGen(ResonanceMCCollisions::iterator const& collision, - aod::ResoMCParents const& mcParents) + template + bool passesProtonDCASelection(TrackType const& track) { - float centralityPercent = collision.cent(); + return passesDCAxyDCAzSelection(track, mProtonDCAPtEdges, mProtonMaxDCAxy); + } - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 0); - if (mcRequireTriggerTVX && !collision.isTriggerTVX()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 1); - if (mcRequireVtxWithin10cm && !collision.isVtxIn10()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 2); - if (mcRequireINELgt0 && !collision.isINELgt0()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 3); - if (mcRequireSel8 && !collision.isInSel8()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 4); - if (cEvtRecINELgt0 && !collision.isRecINELgt0()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 5); - if (mcRequireAfterAllCuts && !collision.isInAfterAllCuts()) - return; - allHistograms.fill(HIST("SignalLoss/mcEventSelectionCutflow"), 6); + template + bool passesKaonDCASelection(TrackType const& track) + { + return passesDCAxyDCAzSelection(track, mKaonDCAPtEdges, mKaonMaxDCAxy); + } - // FIX const-ref-in-for-loop: use const auto& in all range-based for loops - for (const auto& parentParticle : mcParents) { - if (parentParticle.y() < minPairRapidity || - parentParticle.y() > maxPairRapidity) - continue; + bool passesMomentumDependentPID(float p, bool hasTOF, float tpcNSigma, float tofNSigma, SpeciesPIDSettings const& cuts) + { + const float absTPC = std::abs(tpcNSigma); + if (pidShared.useTPCOnlyPIDAtAllMomenta) { + return absTPC < cuts.maxTPCNSigma; + } + if (!hasTOF) { + if (p >= cuts.thresholdMomentum) { + return false; + } + const int bin = findBinIndex(cuts.tpcOnlyBins, p); + return bin >= 0 && absTPC < cuts.tpcOnlyCuts[bin]; + } + const float absTOF = std::abs(tofNSigma); + const int bin = findBinIndex(cuts.combinedBins, p); + return bin >= 0 && (absTPC * absTPC + absTOF * absTOF) < cuts.combinedCuts[bin] * cuts.combinedCuts[bin]; + } - int pdgCode = parentParticle.pdgCode(); - float parentPt = parentParticle.pt(); - double mTScaledPtSquared = -1.0; - - std::array parentMom = { - parentParticle.px(), parentParticle.py(), parentParticle.pz()}; - float parentMass = RecoDecay::m(parentMom, parentParticle.e()); - - // FIX pdg/explicit-mass: use o2::constants::physics::MassLambda1520 - // (already done via using) - auto computeMtScaledPtSquared = [&]() -> double { - return (parentPt * parentPt) + (parentMass * parentMass) - - (MassLambda1520 * MassLambda1520); - }; - - // FIX pdg/explicit-code: use named PDG constants from lambda_analysis - // namespace - if (pdgCode == lambda_analysis::PdgProton) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromProton"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == -lambda_analysis::PdgProton) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromAntiProton"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == lambda_analysis::PdgLambda0) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromLambda0"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == -lambda_analysis::PdgLambda0) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromAntiLambda0"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == lambda_analysis::PdgXiMinus) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromXiMinus"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == -lambda_analysis::PdgXiMinus) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromXiPlus"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == lambda_analysis::PdgXi0) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromXi0"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == -lambda_analysis::PdgXi0) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromAntiXi0"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == lambda_analysis::PdgOmegaMinus) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromOmegaMinus"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } else if (pdgCode == -lambda_analysis::PdgOmegaMinus) { - mTScaledPtSquared = computeMtScaledPtSquared(); - if (mTScaledPtSquared > 0) - allHistograms.fill(HIST("SignalLoss/mTScaled_fromOmegaPlus"), - std::sqrt(mTScaledPtSquared), centralityPercent); - } - } - } - PROCESS_SWITCH(lambda1520analysispo, processMCGen, - "Generator-level MC signal loss study (mT scaling)", false); - - // ── Event-mixing binning types ─────────────────────────────────────────── - // FIX name/type: UpperCamelCase type aliases - using MixingBinningVtxZAndCentrality = - ColumnBinningPolicy; - - // ============================================================ - // processMix() - // FIX const-ref-in-process: collisions argument is const& - // ============================================================ - void processMix(ResonanceCollisionsWithEP const& collisions, - ResonanceTrackTable const& tracks) - { - LOGF(debug, "Event mixing started"); - MixingBinningVtxZAndCentrality mixingBins{ - {vertexZMixingBins, centralityMixingBins}}; - auto trackPool = std::make_tuple(tracks); - SameKindPair - eventPairs{mixingBins, numberOfEventsToMix, -1, - collisions, trackPool, &sliceCache}; - - // FIX const-ref-in-for-loop: const auto& structured binding - for (const auto& [col1, tracks1, col2, tracks2] : eventPairs) { - if (cEvtRecINELgt0 && !col1.isRecINELgt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - continue; - allHistograms.fill(HIST("Event/mixingBins_centralityVsVtxZVsEventPlane"), - col1.cent(), col1.posZ(), col1.evtPl()); - fillInvariantMassHistograms(tracks1, tracks2, col1.cent()); + template + bool passesProtonPID(TrackType const& track, float totalMomentum) + { + if (protonPID.requireTOFForProton && !track.hasTOF()) { + return false; } + const bool hasTOF = track.hasTOF(); + const float tofNSigma = hasTOF ? track.tofNSigmaPr() : -999.f; + return passesMomentumDependentPID(totalMomentum, hasTOF, track.tpcNSigmaPr(), tofNSigma, mProtonPIDCuts); } - PROCESS_SWITCH(lambda1520analysispo, processMix, - "Event mixing for background estimation (standard format)", - true); - // ── Merged-DF type aliases ─────────────────────────────────────────────── - Preslice tracksPerMergedDFCollision = - aod::resodaughter::resoCollisionDFId; - using MergedDFCollisions = aod::ResoCollisionDFs; - using MergedDFTracks = aod::ResoTrackDFs; + template + bool isPionContaminatedKaonCandidate(TrackType const& track) + { + return track.pt() > pidShared.pionRejectionMinPt && + track.pt() < pidShared.pionRejectionMaxPt && + !track.hasTOF() && + std::abs(track.tpcNSigmaPi()) < pidShared.pionRejectionMaxTPCNSigma; + } - // ============================================================ - // processDatadf() - // ============================================================ - void processDatadf(MergedDFCollisions::iterator const& collision, - MergedDFTracks const& tracks) + template + bool passesKaonPID(TrackType const& track, float totalMomentum) { - if (doprocessData) - LOG(error) << "Disable processData() first when using processDatadf()!"; + if (kaonPID.requireTOFForKaon && !track.hasTOF()) { + return false; + } + if (pidShared.applyKaonPionRejection && isPionContaminatedKaonCandidate(track)) { + return false; + } + const bool hasTOF = track.hasTOF(); + const float tofNSigma = hasTOF ? track.tofNSigmaKa() : -999.f; + return passesMomentumDependentPID(totalMomentum, hasTOF, track.tpcNSigmaKa(), tofNSigma, mKaonPIDCuts); + } - if (cEvtRecINELgt0 && !collision.isRecINELgt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - return; - auto occupancyValue = 100; - if (applyOccupancyInTimeRangeCut) - occupancyValue = collision.trackOccupancyInTimeRange(); - - allHistograms.fill(HIST("Event/centralityVsOccupancy"), collision.cent(), - occupancyValue); - for (const auto& track : tracks) { - fillTrackQAHistograms(track, collision.cent()); - } - fillInvariantMassHistograms(tracks, tracks, collision.cent(), - occupancyValue); - } - PROCESS_SWITCH(lambda1520analysispo, processDatadf, - "Process real data in merged derived-data (DF) format", false); - - // ============================================================ - // processMixDF() - // FIX const-ref-in-process: collisions argument is const& - // ============================================================ - // FIX name/type: UpperCamelCase type alias - using MixingBinningDF = - ColumnBinningPolicy; - - void processMixDF(MergedDFCollisions const& collisions, - MergedDFTracks const& tracks) - { - if (doprocessMix) - LOG(fatal) << "Disable processMix() first when using processMixDF()!"; - LOGF(debug, "Event mixing (DF format) started"); - - MixingBinningDF mixingBins{{vertexZMixingBins, centralityMixingBins}}; - auto trackPool = std::make_tuple(tracks); - SameKindPair - eventPairs{mixingBins, numberOfEventsToMix, -1, - collisions, trackPool, &sliceCache}; - - // FIX const-ref-in-for-loop: const auto& structured binding - for (const auto& [col1, tracks1, col2, tracks2] : eventPairs) { - if (cEvtRecINELgt0 && !col1.isRecINELgt0()) // Check reco INELgt0 (at least one PV track in |eta| < 1) about the collision - continue; - auto occupancyValue = 100; - if (applyOccupancyInTimeRangeCut) - occupancyValue = col1.trackOccupancyInTimeRange(); - - allHistograms.fill(HIST("Event/mixingBins_centralityVsVtxZVsEventPlane"), - col1.cent(), col1.posZ(), col1.evtPl()); - fillInvariantMassHistograms(tracks1, tracks2, col1.cent(), - occupancyValue); - } - } - PROCESS_SWITCH(lambda1520analysispo, processMixDF, - "Event mixing for DF-format data", false); - - // ============================================================ - // processMixepDF() - // FIX const-ref-in-process: collisions argument is const& - // ============================================================ - // FIX name/type: UpperCamelCase type alias - using MixingBinningWithEventPlane = - ColumnBinningPolicy; - - void processMixepDF(MergedDFCollisions const& collisions, - MergedDFTracks const& tracks) - { - if (doprocessMix || doprocessMixDF) - LOG(fatal) << "Disable processMix() or processMixDF() first!"; - LOGF(debug, "Event-plane-dependent event mixing (DF format) started"); - - MixingBinningWithEventPlane mixingBins{ - {vertexZMixingBins, centralityMixingBins, eventPlaneMixingBins}}; - auto trackPool = std::make_tuple(tracks); - SameKindPair - eventPairs{mixingBins, numberOfEventsToMix, -1, - collisions, trackPool, &sliceCache}; - - // FIX const-ref-in-for-loop: const auto& structured binding - for (const auto& [col1, tracks1, col2, tracks2] : eventPairs) { - allHistograms.fill(HIST("Event/mixingBins_centralityVsVtxZVsEventPlane"), - col1.cent(), col1.posZ(), col1.evtPl()); - fillInvariantMassHistograms(tracks1, tracks2, col1.cent()); - } - } - PROCESS_SWITCH(lambda1520analysispo, processMixepDF, - "Event-plane dependent event mixing for DF-format data", - false); + template + DaughterTrackCandidate makeTrackCandidate(TrackType const& track) + { + DaughterTrackCandidate cand{}; + cand.px = track.px(); + cand.py = track.py(); + cand.pz = track.pz(); + cand.pt = track.pt(); + cand.eta = track.eta(); + cand.phi = track.phi(); + cand.mom = RecoDecay::p(track.px(), track.py(), track.pz()); + cand.dcaXY = track.dcaXY(); + cand.dcaZ = track.dcaZ(); + cand.sign = track.sign(); + cand.hasTOF = track.hasTOF(); + cand.tpcNSigmaPr = track.tpcNSigmaPr(); + cand.tpcNSigmaKa = track.tpcNSigmaKa(); + if (cand.hasTOF) { + cand.tofNSigmaPr = track.tofNSigmaPr(); + cand.tofNSigmaKa = track.tofNSigmaKa(); + } + cand.globalIndex = track.globalIndex(); + return cand; + } + + template + void attachMCTruthToCandidate(DaughterTrackCandidate& cand, TrackType const& track) + { + if (!track.has_mcParticle()) { + return; + } + auto const mcPart = track.mcParticle(); + cand.hasMcLabel = true; + cand.mcPdgCode = mcPart.pdgCode(); + if (!mcPart.has_mothers()) { + return; + } + for (auto const& mother : mcPart.template mothers_as()) { + if (std::abs(mother.pdgCode()) == mLambda1520AbsPdgCode) { + cand.mcLambda1520MotherIndex = mother.globalIndex(); + cand.mcLambda1520MotherPdgCode = mother.pdgCode(); + break; + } + } + } + + template + void buildProtonKaonCandidatePools(TrackCollection const& tracks, std::vector& protons, std::vector& kaons) + { + protons.clear(); + kaons.clear(); + for (auto const& track : tracks) { + if (!passesBasicTrackSelection(track)) { + continue; + } + const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); + const bool isProton = passesProtonPID(track, mom) && passesProtonDCASelection(track); + const bool isKaon = passesKaonPID(track, mom) && passesKaonDCASelection(track); + if (!isProton && !isKaon) { + continue; + } + DaughterTrackCandidate cand = makeTrackCandidate(track); + if constexpr (IsMC) { + attachMCTruthToCandidate(cand, track); + } + if (isProton) { + protons.push_back(cand); + } + if (isKaon) { + kaons.push_back(cand); + } + } + } + + void computeDaughterPairGeometry(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka, float& deltaEta, float& deltaPhi, float& openingAngle) + { + deltaEta = std::abs(pr.eta - ka.eta); + deltaPhi = std::abs(pr.phi - ka.phi); + if (deltaPhi > o2::constants::math::PI) { + deltaPhi = o2::constants::math::TwoPI - deltaPhi; + } + const float dot = pr.px * ka.px + pr.py * ka.py + pr.pz * ka.pz; + openingAngle = std::acos(std::clamp(dot / (pr.mom * ka.mom), -1.f, 1.f)); + } + + bool passesDaughterAngularCuts(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka) + { + if (!pairCuts.applyDaughterDeltaEtaPhiCuts && !pairCuts.applyDeepAngleCut && !pairCuts.applyOpeningAngleVeto) { + return true; + } + float deltaEta = 0.f, deltaPhi = 0.f, openingAngle = 0.f; + computeDaughterPairGeometry(pr, ka, deltaEta, deltaPhi, openingAngle); + if (pairCuts.applyDaughterDeltaEtaPhiCuts) { + const float maxDaughterDeltaEtaValue = pairCuts.maxDaughterDeltaEta; + const float maxDaughterDeltaPhiValue = pairCuts.maxDaughterDeltaPhi; + if (deltaEta >= maxDaughterDeltaEtaValue || deltaPhi >= maxDaughterDeltaPhiValue) { + return false; + } + } + if (pairCuts.applyDeepAngleCut) { + const float minOpeningAngle = pairCuts.deepAngleCutValue; + if (openingAngle < minOpeningAngle) { + return false; + } + } + if (pairCuts.applyOpeningAngleVeto) { + const float openingAngleVetoLow = pairCuts.openingAngleVetoMin; + const float openingAngleVetoHigh = pairCuts.openingAngleVetoMax; + if (openingAngle > openingAngleVetoLow && openingAngle < openingAngleVetoHigh) { + return false; + } + } + return true; + } + + bool buildProtonKaonPair(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka, float& pairMass, float& pairPt, float& pairY) + { + if (!passesDaughterAngularCuts(pr, ka)) { + return false; + } + const std::array, 2> momenta = { + std::array{pr.px, pr.py, pr.pz}, + std::array{ka.px, ka.py, ka.pz}}; + pairMass = RecoDecay::m(momenta, kProtonKaonMasses); + pairPt = RecoDecay::pt(std::array{pr.px + ka.px, pr.py + ka.py}); + pairY = RecoDecay::y(std::array{pr.px + ka.px, pr.py + ka.py, pr.pz + ka.pz}, pairMass); + + const float minRapidity = trackCuts.cfgMinY; + const float maxRapidity = trackCuts.cfgMaxY; + if (pairY < minRapidity || pairY > maxRapidity) { + return false; + } + if (pairCuts.applyPairPtMassUpperLimits) { + const float maxPairPt = pairCuts.cfgMaxPairPt; + const float maxPairMass = pairCuts.cfgMaxPairMass; + if (pairPt >= maxPairPt || pairMass >= maxPairMass) { + return false; + } + } + return true; + } + + void fillPairQABeforeCuts(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka) + { + float deltaEta = 0.f, deltaPhi = 0.f, openingAngle = 0.f; + computeDaughterPairGeometry(pr, ka, deltaEta, deltaPhi, openingAngle); + histos.fill(HIST("PairQA/Data/hOpeningAngleBeforeCuts"), openingAngle); + histos.fill(HIST("PairQA/Data/hDeltaEtaBeforeCuts"), deltaEta); + histos.fill(HIST("PairQA/Data/hDeltaPhiBeforeCuts"), deltaPhi); + } + + void fillPairQAAfterCuts(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka) + { + float deltaEta = 0.f, deltaPhi = 0.f, openingAngle = 0.f; + computeDaughterPairGeometry(pr, ka, deltaEta, deltaPhi, openingAngle); + histos.fill(HIST("PairQA/Data/hOpeningAngleAfterCuts"), openingAngle); + histos.fill(HIST("PairQA/Data/hDeltaEtaAfterCuts"), deltaEta); + histos.fill(HIST("PairQA/Data/hDeltaPhiAfterCuts"), deltaPhi); + } + + static float vecDot(const std::array& a, const std::array& b) + { + return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; + } + + static std::array vecSub(const std::array& a, const std::array& b) + { + return {a[0] - b[0], a[1] - b[1], a[2] - b[2]}; + } + + static std::array vecCross(const std::array& a, const std::array& b) + { + return {a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]}; + } + + static constexpr float VecUnitEpsilon = 1e-9f; + static std::array vecUnit(const std::array& a) + { + float mag = std::sqrt(vecDot(a, a)); + if (mag < VecUnitEpsilon) { + return {0.f, 0.f, 0.f}; + } + return {a[0] / mag, a[1] / mag, a[2] / mag}; + } + + static constexpr float BoostMomentumEpsilon = 1e-6f; + static std::array boostToRestFrame(float motherE, const std::array& motherP, + float dE, const std::array& dP) + { + float pMag2 = vecDot(motherP, motherP); + float pMag = std::sqrt(pMag2); + if (pMag < BoostMomentumEpsilon) { + return dP; + } + + std::array nHat = {motherP[0] / pMag, motherP[1] / pMag, motherP[2] / pMag}; + float beta = pMag / motherE; + float motherMass2 = motherE * motherE - pMag2; + float gamma = motherMass2 > 0.f ? motherE / std::sqrt(motherMass2) : 1.f; + + float pParallel = vecDot(dP, nHat); + std::array pPerp = {dP[0] - pParallel * nHat[0], dP[1] - pParallel * nHat[1], dP[2] - pParallel * nHat[2]}; + float pParallelBoosted = gamma * (pParallel - beta * dE); + + return {pPerp[0] + pParallelBoosted * nHat[0], pPerp[1] + pParallelBoosted * nHat[1], pPerp[2] + pParallelBoosted * nHat[2]}; + } + + static DecayKinematics makeDecayKinematics(std::array const& protonMom, std::array const& kaonMom) + { + DecayKinematics decay; + decay.protonMom = protonMom; + decay.kaonMom = kaonMom; + decay.motherMom = {protonMom[0] + kaonMom[0], protonMom[1] + kaonMom[1], protonMom[2] + kaonMom[2]}; + decay.protonE = std::sqrt(vecDot(protonMom, protonMom) + ProtonMass * ProtonMass); + decay.kaonE = std::sqrt(vecDot(kaonMom, kaonMom) + KaonMass * KaonMass); + decay.motherE = decay.protonE + decay.kaonE; + return decay; + } + + static DecayKinematics makeDecayKinematics(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka) + { + return makeDecayKinematics(std::array{pr.px, pr.py, pr.pz}, std::array{ka.px, ka.py, ka.pz}); + } + + template + static DecayKinematics makeGeneratedDecayKinematics(McParticleType const& mother, McParticleType const& proton, McParticleType const& kaon) + { + DecayKinematics decay; + decay.motherMom = {mother.px(), mother.py(), mother.pz()}; + decay.motherE = mother.e(); + decay.protonMom = {proton.px(), proton.py(), proton.pz()}; + decay.protonE = proton.e(); + decay.kaonMom = {kaon.px(), kaon.py(), kaon.pz()}; + decay.kaonE = kaon.e(); + return decay; + } + + template + static float computeGeneratedMass(McParticleType const& particle) + { + const float energy = particle.e(); + const float px = particle.px(); + const float py = particle.py(); + const float pz = particle.pz(); + return std::sqrt(std::max(0.f, energy * energy - px * px - py * py - pz * pz)); + } + + float computeCosThetaStar(DecayKinematics const& decay) + { + const std::array protonCM = boostToRestFrame(decay.motherE, decay.motherMom, decay.protonE, decay.protonMom); + const std::array kaonCM = boostToRestFrame(decay.motherE, decay.motherMom, decay.kaonE, decay.kaonMom); + const std::array referenceCM = vecUnit(spinAlignment.referenceFrameDaughter == 0 ? protonCM : kaonCM); + + if (spinAlignment.activateHelicityFrame) { + return vecDot(referenceCM, vecUnit(decay.motherMom)); + } + + if (spinAlignment.activateCollinsSoperFrame || spinAlignment.activateProductionFrame || spinAlignment.activateBeamAxisFrame) { + const std::array beam1CM = vecUnit(boostToRestFrame(decay.motherE, decay.motherMom, mBeam1Energy, mBeam1Momentum)); + const std::array beam2CM = vecUnit(boostToRestFrame(decay.motherE, decay.motherMom, mBeam2Energy, mBeam2Momentum)); + if (spinAlignment.activateCollinsSoperFrame) { + return vecDot(referenceCM, vecUnit(vecSub(beam1CM, beam2CM))); + } + if (spinAlignment.activateProductionFrame) { + return vecDot(referenceCM, vecUnit(vecCross(beam1CM, vecUnit(decay.motherMom)))); + } + return vecDot(referenceCM, beam1CM); + } + + if (spinAlignment.activateRandomFrame) { + const double randomCos = mUniform01(mRandomEngine) * 2. - 1.; + const double randomPhi = mUniform01(mRandomEngine) * o2::constants::math::TwoPI; + const double randomSin = std::sqrt(1. - randomCos * randomCos); + const std::array randomAxis = {static_cast(randomSin * std::cos(randomPhi)), + static_cast(randomSin * std::sin(randomPhi)), + static_cast(randomCos)}; + return vecDot(referenceCM, randomAxis); + } + + return 0.f; + } + + template + void fillCosThetaStarHistograms(float candidateMass, float candidatePt, DecayKinematics const& decay, float centrality, bool isParticle) + { + if (!isSpinAlignmentActive()) { + return; + } + + const float cosTheta = computeCosThetaStar(decay); + + if constexpr (Sample == CosThetaStarSample::Signal) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/Data/SameEvent/hLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/Data/SameEvent/hAntiLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } + } else if constexpr (Sample == CosThetaStarSample::LikeSign) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/Data/SameEvent/hLikeSignProtonPlusKaonPlus"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/Data/SameEvent/hLikeSignProtonMinusKaonMinus"), candidateMass, candidatePt, centrality, cosTheta); + } + } else if constexpr (Sample == CosThetaStarSample::Rotated) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/Data/RotatedBackground/hLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/Data/RotatedBackground/hAntiLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } + } else if constexpr (Sample == CosThetaStarSample::Mixed) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/Data/MixedEvent/hLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/Data/MixedEvent/hAntiLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } + } else if constexpr (Sample == CosThetaStarSample::LikeSignMixed) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/Data/MixedEvent/hLikeSignProtonPlusKaonPlus"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/Data/MixedEvent/hLikeSignProtonMinusKaonMinus"), candidateMass, candidatePt, centrality, cosTheta); + } + } else if constexpr (Sample == CosThetaStarSample::TruthMatched) { + if (isParticle) { + histos.fill(HIST("CosThetaStar/MCReco/TruthMatched/hLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/MCReco/TruthMatched/hAntiLambda1520"), candidateMass, candidatePt, centrality, cosTheta); + } + } + } + + void fillInvariantMassSameEvent(int protonSign, int kaonSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (isSpinAlignmentActive()) { + return; + } + if (protonSign * kaonSign < 0) { + if (protonSign > 0) { + histos.fill(HIST("InvMass/Data/SameEvent/hLambda1520"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/Data/SameEvent/hAntiLambda1520"), pairMass, pairPt, centrality, rapidity); + } + } else { + if (protonSign > 0) { + histos.fill(HIST("InvMass/Data/SameEvent/hLikeSignProtonPlusKaonPlus"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/Data/SameEvent/hLikeSignProtonMinusKaonMinus"), pairMass, pairPt, centrality, rapidity); + } + } + } + + void fillInvariantMassMixedEvent(int protonSign, int kaonSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (isSpinAlignmentActive()) { + return; + } + if (protonSign * kaonSign < 0) { + if (protonSign > 0) { + histos.fill(HIST("InvMass/Data/MixedEvent/hLambda1520"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/Data/MixedEvent/hAntiLambda1520"), pairMass, pairPt, centrality, rapidity); + } + } else if (mixingCfg.cfgFillMixLikeSign) { + if (protonSign > 0) { + histos.fill(HIST("InvMass/Data/MixedEvent/hLikeSignProtonPlusKaonPlus"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/Data/MixedEvent/hLikeSignProtonMinusKaonMinus"), pairMass, pairPt, centrality, rapidity); + } + } + } + + void fillInvariantMassMCRecoAllPairs(int protonSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (protonSign > 0) { + histos.fill(HIST("InvMass/MCReco/AllPairs/hLambda1520"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/MCReco/AllPairs/hAntiLambda1520"), pairMass, pairPt, centrality, rapidity); + } + } + + void fillInvariantMassMCRecoTruthMatched(int motherPdg, float pairPt, float pairMass, float centrality, float rapidity) + { + if (motherPdg > 0) { + histos.fill(HIST("InvMass/MCReco/TruthMatched/hLambda1520"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/MCReco/TruthMatched/hAntiLambda1520"), pairMass, pairPt, centrality, rapidity); + } + } + + void fillInvariantMassMCRecoMixedEvent(int protonSign, int kaonSign, float pairPt, float pairMass, float centrality, float rapidity) + { + if (protonSign * kaonSign >= 0) { + return; + } + if (protonSign > 0) { + histos.fill(HIST("InvMass/MCReco/MixedEvent/hLambda1520"), pairMass, pairPt, centrality, rapidity); + } else { + histos.fill(HIST("InvMass/MCReco/MixedEvent/hAntiLambda1520"), pairMass, pairPt, centrality, rapidity); + } + } + + void fillRotatedPairBackground(DaughterTrackCandidate const& pr, DaughterTrackCandidate const& ka, float centrality) + { + const int nRot = rotBkg.numberOfRotations; + if (nRot <= 0) { + return; + } + const float weight = 1.f / static_cast(nRot); + const float rotationHalfWindow = o2::constants::math::PI / static_cast(rotBkg.rotationAngleWindow); + const bool rotateProtonInsteadOfKaon = rotBkg.rotateProtonInsteadOfKaon; + const DaughterTrackCandidate& rotated = rotateProtonInsteadOfKaon ? pr : ka; + + for (int iRot = 0; iRot < nRot; ++iRot) { + float rotationAngle = o2::constants::math::PI; + if (nRot > 1) { + rotationAngle = (o2::constants::math::PI - rotationHalfWindow) + + static_cast(iRot) * (2.f * rotationHalfWindow / static_cast(nRot - 1)); + } + const float rotatedPhi = RecoDecay::constrainAngle(rotated.phi + rotationAngle, 0.f); + const float rotatedPx = rotated.pt * std::cos(rotatedPhi); + const float rotatedPy = rotated.pt * std::sin(rotatedPhi); + + const std::array momPr = rotateProtonInsteadOfKaon ? std::array{rotatedPx, rotatedPy, pr.pz} + : std::array{pr.px, pr.py, pr.pz}; + const std::array momKa = rotateProtonInsteadOfKaon ? std::array{ka.px, ka.py, ka.pz} + : std::array{rotatedPx, rotatedPy, ka.pz}; + const std::array, 2> rotMomenta = {momPr, momKa}; + + const float rotatedMass = RecoDecay::m(rotMomenta, kProtonKaonMasses); + const float pairPxSum = momPr[0] + momKa[0]; + const float pairPySum = momPr[1] + momKa[1]; + const float pairPzSum = momPr[2] + momKa[2]; + const float rotatedPt = RecoDecay::pt(std::array{pairPxSum, pairPySum}); + const float rotatedRapidity = RecoDecay::y(std::array{pairPxSum, pairPySum, pairPzSum}, rotatedMass); + + const float minRapidity = trackCuts.cfgMinY; + const float maxRapidity = trackCuts.cfgMaxY; + if (rotatedRapidity < minRapidity || rotatedRapidity > maxRapidity) { + continue; + } + if (pairCuts.applyPairPtMassUpperLimits) { + const float maxPairPt = pairCuts.cfgMaxPairPt; + const float maxPairMass = pairCuts.cfgMaxPairMass; + if (rotatedPt >= maxPairPt || rotatedMass >= maxPairMass) { + continue; + } + } + + if (!isSpinAlignmentActive()) { + if (pr.sign > 0) { + histos.fill(HIST("InvMass/Data/RotatedBackground/hLambda1520"), rotatedMass, rotatedPt, centrality, rotatedRapidity, weight); + } else { + histos.fill(HIST("InvMass/Data/RotatedBackground/hAntiLambda1520"), rotatedMass, rotatedPt, centrality, rotatedRapidity, weight); + } + } + + if (isSpinAlignmentActive() && spinAlignment.fillRotatedBackgroundCosThetaStar) { + fillCosThetaStarHistograms( + rotatedMass, rotatedPt, makeDecayKinematics(momPr, momKa), centrality, pr.sign > 0); + } + } + } + + template + void fillInvariantMassFromCandidatePools( + std::vector const& protonPool, + std::vector const& kaonPool, + float centrality, + bool fillPairQA = false) + { + for (auto const& pr : protonPool) { + for (auto const& ka : kaonPool) { + + if constexpr (!IsMixed) { + if (pr.globalIndex == ka.globalIndex) { + continue; + } + if (fillPairQA) { + fillPairQABeforeCuts(pr, ka); + } + } + + float pairMass = 0.f, pairPt = 0.f, pairY = 0.f; + if (!buildProtonKaonPair(pr, ka, pairMass, pairPt, pairY)) { + continue; + } + + if constexpr (!IsMixed) { + if (fillPairQA) { + fillPairQAAfterCuts(pr, ka); + } + } + + if constexpr (IsMixed) { + if constexpr (IsMC) { + fillInvariantMassMCRecoMixedEvent(pr.sign, ka.sign, pairPt, pairMass, centrality, pairY); + } else { + fillInvariantMassMixedEvent(pr.sign, ka.sign, pairPt, pairMass, centrality, pairY); + + if (isSpinAlignmentActive() && spinAlignment.fillMixedEventCosThetaStar) { + const DecayKinematics decay = makeDecayKinematics(pr, ka); + if (pr.sign * ka.sign < 0) { + fillCosThetaStarHistograms(pairMass, pairPt, decay, centrality, pr.sign > 0); + } else if (mixingCfg.cfgFillMixLikeSign) { + fillCosThetaStarHistograms(pairMass, pairPt, decay, centrality, pr.sign > 0); + } + } + } + } else { + fillInvariantMassSameEvent(pr.sign, ka.sign, pairPt, pairMass, centrality, pairY); + + if (isSpinAlignmentActive()) { + const DecayKinematics decay = makeDecayKinematics(pr, ka); + if (pr.sign * ka.sign < 0) { + fillCosThetaStarHistograms(pairMass, pairPt, decay, centrality, pr.sign > 0); + } else { + fillCosThetaStarHistograms(pairMass, pairPt, decay, centrality, pr.sign > 0); + } + } + + if (rotBkg.enableRotationalBackground && pr.sign * ka.sign < 0) { + fillRotatedPairBackground(pr, ka, centrality); + } + } + } + } + } + + template + void fillSelectedProtonQA(TrackType const& track, float totalMomentum, float centrality) + { + const float pt = track.pt(); + const float mom = totalMomentum; + const float tpcNSigPr = track.tpcNSigmaPr(); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hDCAxyVsPt"), pt, track.dcaXY()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hDCAzVsPt"), pt, track.dcaZ()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTPCvsP"), mom, tpcNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTPCvsPt"), pt, tpcNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTPCvsCentrality"), centrality, tpcNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTPCPionVsPt"), pt, track.tpcNSigmaPi()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTPCKaonVsPt"), pt, track.tpcNSigmaKa()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hTPCSignalVsInnerP"), track.tpcInnerParam(), track.tpcSignal()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hTPCClustersVsPhi"), track.tpcNClsFound(), track.phi()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hTPCCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hTPCFoundClustersVsPt"), pt, track.tpcNClsFound()); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hPt"), pt); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hEta"), track.eta()); + if (track.hasTOF()) { + const float tofNSigPr = track.tofNSigmaPr(); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTOFvsP"), mom, tofNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTOFvsPt"), pt, tofNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTOFvsCentrality"), centrality, tofNSigPr); + histos.fill(HIST("TrackQA/Data/Selected/Proton/hNSigmaTOFvsNSigmaTPC"), tpcNSigPr, tofNSigPr); + } + } + + template + void fillSelectedKaonQA(TrackType const& track, float totalMomentum, float centrality) + { + const float pt = track.pt(); + const float mom = totalMomentum; + const float tpcNSigKa = track.tpcNSigmaKa(); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hDCAxyVsPt"), pt, track.dcaXY()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hDCAzVsPt"), pt, track.dcaZ()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsP"), mom, tpcNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsPt"), pt, tpcNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTPCvsCentrality"), centrality, tpcNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTPCPionVsPt"), pt, track.tpcNSigmaPi()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTPCProtonVsPt"), pt, track.tpcNSigmaPr()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hTPCSignalVsInnerP"), track.tpcInnerParam(), track.tpcSignal()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hTPCClustersVsPhi"), track.tpcNClsFound(), track.phi()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hTPCCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hTPCFoundClustersVsPt"), pt, track.tpcNClsFound()); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hPt"), pt); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hEta"), track.eta()); + if (track.hasTOF()) { + const float tofNSigKa = track.tofNSigmaKa(); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsP"), mom, tofNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsPt"), pt, tofNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsCentrality"), centrality, tofNSigKa); + histos.fill(HIST("TrackQA/Data/Selected/Kaon/hNSigmaTOFvsNSigmaTPC"), tpcNSigKa, tofNSigKa); + } + } + + template + void fillSelectedProtonQAMC(TrackType const& track, float totalMomentum, float centrality) + { + const float pt = track.pt(); + const float mom = totalMomentum; + const float tpcNSigPr = track.tpcNSigmaPr(); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hDCAxyVsPt"), pt, track.dcaXY()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hDCAzVsPt"), pt, track.dcaZ()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsP"), mom, tpcNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsPt"), pt, tpcNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTPCvsCentrality"), centrality, tpcNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTPCPionVsPt"), pt, track.tpcNSigmaPi()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTPCKaonVsPt"), pt, track.tpcNSigmaKa()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hTPCSignalVsInnerP"), track.tpcInnerParam(), track.tpcSignal()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hTPCClustersVsPhi"), track.tpcNClsFound(), track.phi()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hTPCCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hTPCFoundClustersVsPt"), pt, track.tpcNClsFound()); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hPt"), pt); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hEta"), track.eta()); + if (track.hasTOF()) { + const float tofNSigPr = track.tofNSigmaPr(); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsP"), mom, tofNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsPt"), pt, tofNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsCentrality"), centrality, tofNSigPr); + histos.fill(HIST("TrackQA/MCReco/Selected/Proton/hNSigmaTOFvsNSigmaTPC"), tpcNSigPr, tofNSigPr); + } + } + + template + void fillSelectedKaonQAMC(TrackType const& track, float totalMomentum, float centrality) + { + const float pt = track.pt(); + const float mom = totalMomentum; + const float tpcNSigKa = track.tpcNSigmaKa(); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hDCAxyVsPt"), pt, track.dcaXY()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hDCAzVsPt"), pt, track.dcaZ()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsP"), mom, tpcNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsPt"), pt, tpcNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCvsCentrality"), centrality, tpcNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCPionVsPt"), pt, track.tpcNSigmaPi()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTPCProtonVsPt"), pt, track.tpcNSigmaPr()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hTPCSignalVsInnerP"), track.tpcInnerParam(), track.tpcSignal()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hTPCClustersVsPhi"), track.tpcNClsFound(), track.phi()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hTPCCrossedRowsVsPt"), pt, track.tpcNClsCrossedRows()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hTPCFoundClustersVsPt"), pt, track.tpcNClsFound()); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hPt"), pt); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hEta"), track.eta()); + if (track.hasTOF()) { + const float tofNSigKa = track.tofNSigmaKa(); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsP"), mom, tofNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsPt"), pt, tofNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsCentrality"), centrality, tofNSigKa); + histos.fill(HIST("TrackQA/MCReco/Selected/Kaon/hNSigmaTOFvsNSigmaTPC"), tpcNSigKa, tofNSigKa); + } + } + + template + bool findProtonKaonDaughters(McPart const& mother, McPart& proton, McPart& kaon) + { + if (!mother.has_daughters()) { + return false; + } + bool hasProtonDaughter = false, hasKaonDaughter = false; + for (auto const& d : mother.template daughters_as()) { + const int absPdgCode = std::abs(d.pdgCode()); + if (absPdgCode == lambda1520analysispo_analysis::PdgProton) { + hasProtonDaughter = true; + proton = d; + } else if (absPdgCode == lambda1520analysispo_analysis::PdgKaon) { + hasKaonDaughter = true; + kaon = d; + } + } + return hasProtonDaughter && hasKaonDaughter; + } + + struct GeneratedDaughterMatch { + bool found = false; + bool trackQuality = false; + bool pid = false; + bool dca = false; + DaughterTrackCandidate cand{}; + }; + + template + GeneratedDaughterMatch findRecoTrackOfGeneratedDaughter(TracksSlice const& tracksSlice, McPart const& mcDaughter, bool isProton) + { + GeneratedDaughterMatch result; + for (auto const& trk : tracksSlice) { + if (!trk.has_mcParticle() || trk.mcParticle().globalIndex() != mcDaughter.globalIndex()) { + continue; + } + result.found = true; + result.trackQuality = passesBasicTrackSelection(trk); + const float mom = RecoDecay::p(trk.px(), trk.py(), trk.pz()); + if (isProton) { + result.pid = passesProtonPID(trk, mom); + result.dca = passesProtonDCASelection(trk); + } else { + result.pid = passesKaonPID(trk, mom); + result.dca = passesKaonDCASelection(trk); + } + result.cand = makeTrackCandidate(trk); + break; + } + return result; + } + + template + void fillGeneratedToReconstructedCutFlow(TracksSlice const& tracksOfCollision, aod::McParticles const& mcParticles, int64_t mcCollisionId) + { + auto genParts = mcParticles.sliceBy(perMcCollision, mcCollisionId); + for (auto const& gen : genParts) { + if (!gen.producedByGenerator() || std::abs(gen.pdgCode()) != mLambda1520AbsPdgCode) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 0.f); + + auto genProton = gen; + auto genKaon = gen; + if (!findProtonKaonDaughters(gen, genProton, genKaon)) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 1.f); + + const float minRapidity = trackCuts.cfgMinY; + const float maxRapidity = trackCuts.cfgMaxY; + if (gen.y() < minRapidity || gen.y() > maxRapidity) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 2.f); + + const float maxTrackEta = trackCuts.cfgCutEta; + const float minTrackPt = trackCuts.cfgCutPt; + const bool protonInAcceptance = std::abs(genProton.eta()) < maxTrackEta && genProton.pt() > minTrackPt; + const bool kaonInAcceptance = std::abs(genKaon.eta()) < maxTrackEta && genKaon.pt() > minTrackPt; + if (protonInAcceptance) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 3.f); + } + if (kaonInAcceptance) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 4.f); + } + if (!protonInAcceptance || !kaonInAcceptance) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 5.f); + + const auto protonMatch = findRecoTrackOfGeneratedDaughter(tracksOfCollision, genProton, true); + const auto kaonMatch = findRecoTrackOfGeneratedDaughter(tracksOfCollision, genKaon, false); + if (protonMatch.found) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 6.f); + } + if (kaonMatch.found) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 7.f); + } + if (!protonMatch.found || !kaonMatch.found) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 8.f); + + if (protonMatch.trackQuality) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 9.f); + } + if (kaonMatch.trackQuality) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 10.f); + } + if (!protonMatch.trackQuality || !kaonMatch.trackQuality) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 11.f); + + if (protonMatch.pid) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 12.f); + } + if (kaonMatch.pid) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 13.f); + } + if (!protonMatch.pid || !kaonMatch.pid) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 14.f); + + if (protonMatch.dca) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 15.f); + } + if (kaonMatch.dca) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 16.f); + } + if (!protonMatch.dca || !kaonMatch.dca) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 17.f); + + float pairMass = 0.f, pairPt = 0.f, pairY = 0.f; + if (buildProtonKaonPair(protonMatch.cand, kaonMatch.cand, pairMass, pairPt, pairY)) { + histos.fill(HIST("CutFlow/MCReco/hGeneratedToReconstructed"), 18.f); + } + } + } + + static float computeMtScaledPt(float pt, float referenceMass) + { + const float mtScaledPtSquared = pt * pt + referenceMass * referenceMass - Lambda1520Mass * Lambda1520Mass; + return (mtScaledPtSquared > 0.f) ? std::sqrt(mtScaledPtSquared) : -1.f; + } + + template + void fillSignalLossProxySpecies(McPartType const& part, float centrality) + { + const int absPdgCode = std::abs(part.pdgCode()); + const float pt = part.pt(); + if (pt <= 0.f) { + return; + } + + if (absPdgCode == lambda1520analysispo_analysis::PdgKaon) { + const float mtScaledPt = computeMtScaledPt(pt, KaonMass); + const float weight = (mcCfg.applyMtScaleWeight && mtScaledPt > 0.f) ? mtScaledPt / pt : 1.f; + if constexpr (IsAcceptedEvent) { + histos.fill(HIST("SignalLoss/Kaon/hAcceptedEventTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Kaon/hAcceptedEventMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } else { + histos.fill(HIST("SignalLoss/Kaon/hGeneratedTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Kaon/hGeneratedMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } + } else if (absPdgCode == lambda1520analysispo_analysis::PdgProton) { + const float mtScaledPt = computeMtScaledPt(pt, ProtonMass); + const float weight = (mcCfg.applyMtScaleWeight && mtScaledPt > 0.f) ? mtScaledPt / pt : 1.f; + if constexpr (IsAcceptedEvent) { + histos.fill(HIST("SignalLoss/Proton/hAcceptedEventTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Proton/hAcceptedEventMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } else { + histos.fill(HIST("SignalLoss/Proton/hGeneratedTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Proton/hGeneratedMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } + } else if (absPdgCode == lambda1520analysispo_analysis::PdgLambda0) { + const float mtScaledPt = computeMtScaledPt(pt, LambdaMass); + const float weight = (mcCfg.applyMtScaleWeight && mtScaledPt > 0.f) ? mtScaledPt / pt : 1.f; + if constexpr (IsAcceptedEvent) { + histos.fill(HIST("SignalLoss/Lambda0/hAcceptedEventTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Lambda0/hAcceptedEventMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } else { + histos.fill(HIST("SignalLoss/Lambda0/hGeneratedTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Lambda0/hGeneratedMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } + } else if (absPdgCode == lambda1520analysispo_analysis::PdgXiMinus) { + const float mtScaledPt = computeMtScaledPt(pt, XiMass); + const float weight = (mcCfg.applyMtScaleWeight && mtScaledPt > 0.f) ? mtScaledPt / pt : 1.f; + if constexpr (IsAcceptedEvent) { + histos.fill(HIST("SignalLoss/Xi/hAcceptedEventTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Xi/hAcceptedEventMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } else { + histos.fill(HIST("SignalLoss/Xi/hGeneratedTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Xi/hGeneratedMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } + } else if (absPdgCode == lambda1520analysispo_analysis::PdgOmegaMinus) { + const float mtScaledPt = computeMtScaledPt(pt, OmegaMass); + const float weight = (mcCfg.applyMtScaleWeight && mtScaledPt > 0.f) ? mtScaledPt / pt : 1.f; + if constexpr (IsAcceptedEvent) { + histos.fill(HIST("SignalLoss/Omega/hAcceptedEventTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Omega/hAcceptedEventMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } else { + histos.fill(HIST("SignalLoss/Omega/hGeneratedTruePtVsCentrality"), pt, centrality); + if (mtScaledPt > 0.f) { + histos.fill(HIST("SignalLoss/Omega/hGeneratedMtScaledPtVsCentrality"), mtScaledPt, centrality, weight); + } + } + } + } + + Filter acceptanceFilter = (nabs(aod::track::eta) < trackCuts.cfgCutEta && nabs(aod::track::pt) > trackCuts.cfgCutPt); + + using EventCandidates = soa::Join; + using TrackBase = soa::Join; + 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 perMcCollision = aod::mcparticle::mcCollisionId; + + using BinningTypeFT0M = ColumnBinningPolicy; + using BinningTypeFT0A = ColumnBinningPolicy; + using BinningTypeFT0C = ColumnBinningPolicy; + using BinningTypeFV0A = ColumnBinningPolicy; + + void processSameEvent(EventCandidates const& collisions, + TrackCandidates const& tracks, + aod::BCs const&) + { + for (auto const& collision : collisions) { + if (!passesEventSelectionImpl(collision)) { + continue; + } + const float centrality = getCentrality(collision); + const int occupancy = collision.trackOccupancyInTimeRange(); + + histos.fill(HIST("EventQA/Data/hVertexZ"), collision.posZ()); + histos.fill(HIST("EventQA/Data/hPVContributors"), collision.numContrib()); + histos.fill(HIST("EventQA/Data/hCentrality"), centrality); + histos.fill(HIST("EventQA/Data/hOccupancy"), occupancy); + histos.fill(HIST("EventQA/Data/hCentralityVsVertexZ"), collision.posZ(), centrality); + histos.fill(HIST("EventQA/Data/hCentralityVsOccupancy"), occupancy, centrality); + + const int64_t collisionIndex = collision.globalIndex(); + auto tracksOfCollision = tracks.sliceBy(perCol, collisionIndex); + tracksOfCollision.bindExternalIndices(&tracks); + + if (qa.cfgFillMultQA) { + histos.fill(HIST("EventQA/Data/hCentralityVsFilteredTracksVsPVTracks"), centrality, tracksOfCollision.size(), collision.multNTracksPV()); + } + + if (qa.cfgFillTrackQA) { + for (auto const& track : tracksOfCollision) { + const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); + histos.fill(HIST("TrackQA/Data/PreSelection/Proton/hNSigmaTPCvsP"), mom, track.tpcNSigmaPr()); + histos.fill(HIST("TrackQA/Data/PreSelection/Kaon/hNSigmaTPCvsP"), mom, track.tpcNSigmaKa()); + if (track.hasTOF()) { + histos.fill(HIST("TrackQA/Data/PreSelection/Proton/hNSigmaTOFvsP"), mom, track.tofNSigmaPr()); + histos.fill(HIST("TrackQA/Data/PreSelection/Proton/hNSigmaTOFvsNSigmaTPC"), track.tpcNSigmaPr(), track.tofNSigmaPr()); + histos.fill(HIST("TrackQA/Data/PreSelection/Kaon/hNSigmaTOFvsP"), mom, track.tofNSigmaKa()); + histos.fill(HIST("TrackQA/Data/PreSelection/Kaon/hNSigmaTOFvsNSigmaTPC"), track.tpcNSigmaKa(), track.tofNSigmaKa()); + } + + histos.fill(HIST("TrackQA/Data/AllTracks/hEta"), track.eta()); + histos.fill(HIST("TrackQA/Data/AllTracks/hPt"), track.pt()); + histos.fill(HIST("TrackQA/Data/AllTracks/hPhi"), track.phi()); + histos.fill(HIST("TrackQA/Data/AllTracks/hEtaPhi"), track.eta(), track.phi()); + + if (!passesBasicTrackSelection(track)) { + continue; + } + if (passesProtonPID(track, mom) && passesProtonDCASelection(track)) { + fillSelectedProtonQA(track, mom, centrality); + } + if (passesKaonPID(track, mom) && passesKaonDCASelection(track)) { + fillSelectedKaonQA(track, mom, centrality); + } + } + } + + buildProtonKaonCandidatePools(tracksOfCollision, mProtonCandidates, mKaonCandidates); + fillInvariantMassFromCandidatePools(mProtonCandidates, mKaonCandidates, centrality, qa.cfgFillPairQA); + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processSameEvent, "Process same event (data)", true); + + template + void loopOverMixedEventPairs(PairType& mixingPair) + { + for (auto const& [c1, tracks1, c2, tracks2] : mixingPair) { + if (!passesEventSelection(c1) || !passesEventSelection(c2)) { + continue; + } + const float centrality = getCentrality(c1); + + if constexpr (IsMC) { + histos.fill(HIST("MixingQA/MCReco/hCollisionPairs"), 1.f); + histos.fill(HIST("MixingQA/MCReco/hCollisionPairsVsCentrality"), centrality); + } else { + if (qa.cfgFillMixQA) { + histos.fill(HIST("MixingQA/Data/hPoolVertexZ"), c1.posZ()); + histos.fill(HIST("MixingQA/Data/hPoolCentrality"), centrality); + histos.fill(HIST("MixingQA/Data/hPoolVertexZVsCentrality"), c1.posZ(), centrality); + histos.fill(HIST("MixingQA/Data/hCollisionPairs"), 1.f); + } + } + + std::vector protonPool1, kaonPool1, protonPool2, kaonPool2; + buildProtonKaonCandidatePools(tracks1, protonPool1, kaonPool1); + buildProtonKaonCandidatePools(tracks2, protonPool2, kaonPool2); + fillInvariantMassFromCandidatePools(protonPool1, kaonPool2, centrality); + fillInvariantMassFromCandidatePools(protonPool2, kaonPool1, centrality); + } + } + + template + void runMixingWithBinning(EventCandidates const& collisions, TrackCandidates const& tracks) + { + BinningType binning{{axes.cfgVtxAxis, axes.cfgCentAxis}, true}; + const int nMixedEvents = mixingCfg.cfgNoMixedEvents; + auto tracksTuple = std::make_tuple(tracks); + SameKindPair mixingPair{binning, nMixedEvents, -1, collisions, tracksTuple, &cache}; + loopOverMixedEventPairs(mixingPair); + } + + template + void runMixedEventForSelectedEstimator(EventCandidates const& collisions, TrackCandidates const& tracks) + { + switch (evSel.cfgCentralityEstimator) { + case lambda1520analysispo_analysis::kFT0M: + runMixingWithBinning(collisions, tracks); + break; + case lambda1520analysispo_analysis::kFT0A: + runMixingWithBinning(collisions, tracks); + break; + case lambda1520analysispo_analysis::kFT0C: + runMixingWithBinning(collisions, tracks); + break; + case lambda1520analysispo_analysis::kFV0A: + runMixingWithBinning(collisions, tracks); + break; + default: + runMixingWithBinning(collisions, tracks); + break; + } + } + + void processMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) + { + runMixedEventForSelectedEstimator(collisions, tracks); + } + PROCESS_SWITCH(Lambda1520Analysispo, processMixedEvent, "Process mixed event (data)", true); + + void processMCEventMixing(EventCandidates const& collisions, TrackCandidates const& tracks) + { + runMixedEventForSelectedEstimator(collisions, tracks); + } + PROCESS_SWITCH(Lambda1520Analysispo, processMCEventMixing, "Process MC reconstructed event mixing (reuses the mixing pairs; no truth matching)", false); + + void processMC(EventCandidatesMC const& collisions, + aod::BCs const&, + TrackCandidatesMC const& tracks, + aod::McParticles const& mcParticles) + { + std::unordered_set mcCollisionsWithEfficiencyQA; + + for (auto const& collision : collisions) { + if (!passesEventSelectionImpl(collision)) { + continue; + } + if (!collision.has_mcCollision()) { + continue; + } + const float centrality = getCentrality(collision); + + histos.fill(HIST("EventQA/MCReco/hAcceptedEventsVsCentrality"), centrality); + histos.fill(HIST("EventQA/MCReco/hVertexZ"), collision.posZ()); + histos.fill(HIST("EventQA/MCReco/hPVContributors"), collision.numContrib()); + + const int64_t collisionIndex = collision.globalIndex(); + auto tracksOfCollision = tracks.sliceBy(perColMC, collisionIndex); + tracksOfCollision.bindExternalIndices(&tracks); + tracksOfCollision.bindExternalIndices(&mcParticles); + + for (auto const& track : tracksOfCollision) { + histos.fill(HIST("TrackQA/MCReco/AllTracks/hEta"), track.eta()); + histos.fill(HIST("TrackQA/MCReco/AllTracks/hPt"), track.pt()); + histos.fill(HIST("TrackQA/MCReco/AllTracks/hPhi"), track.phi()); + histos.fill(HIST("TrackQA/MCReco/AllTracks/hEtaPhi"), track.eta(), track.phi()); + + if (qa.cfgFillTrackQA) { + const float mom = RecoDecay::p(track.px(), track.py(), track.pz()); + histos.fill(HIST("TrackQA/MCReco/PreSelection/Proton/hNSigmaTPCvsP"), mom, track.tpcNSigmaPr()); + histos.fill(HIST("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTPCvsP"), mom, track.tpcNSigmaKa()); + if (track.hasTOF()) { + histos.fill(HIST("TrackQA/MCReco/PreSelection/Proton/hNSigmaTOFvsP"), mom, track.tofNSigmaPr()); + histos.fill(HIST("TrackQA/MCReco/PreSelection/Proton/hNSigmaTOFvsNSigmaTPC"), track.tpcNSigmaPr(), track.tofNSigmaPr()); + histos.fill(HIST("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTOFvsP"), mom, track.tofNSigmaKa()); + histos.fill(HIST("TrackQA/MCReco/PreSelection/Kaon/hNSigmaTOFvsNSigmaTPC"), track.tpcNSigmaKa(), track.tofNSigmaKa()); + } + if (passesBasicTrackSelection(track)) { + if (passesProtonPID(track, mom) && passesProtonDCASelection(track)) { + fillSelectedProtonQAMC(track, mom, centrality); + } + if (passesKaonPID(track, mom) && passesKaonDCASelection(track)) { + fillSelectedKaonQAMC(track, mom, centrality); + } + } + } + } + + buildProtonKaonCandidatePools(tracksOfCollision, mProtonCandidates, mKaonCandidates); + + for (auto const& pr : mProtonCandidates) { + if (pr.mcLambda1520MotherIndex >= 0 && std::abs(pr.mcPdgCode) == lambda1520analysispo_analysis::PdgProton) { + if (pr.mcLambda1520MotherPdgCode > 0) { + histos.fill(HIST("MCRecoQA/DaughterPt/hProtonFromLambda1520"), pr.pt); + } else { + histos.fill(HIST("MCRecoQA/DaughterPt/hProtonFromAntiLambda1520"), pr.pt); + } + } + } + for (auto const& ka : mKaonCandidates) { + if (ka.mcLambda1520MotherIndex >= 0 && std::abs(ka.mcPdgCode) == lambda1520analysispo_analysis::PdgKaon) { + if (ka.mcLambda1520MotherPdgCode > 0) { + histos.fill(HIST("MCRecoQA/DaughterPt/hKaonFromLambda1520"), ka.pt); + } else { + histos.fill(HIST("MCRecoQA/DaughterPt/hKaonFromAntiLambda1520"), ka.pt); + } + } + } + + for (auto const& pr : mProtonCandidates) { + for (auto const& ka : mKaonCandidates) { + if (pr.globalIndex == ka.globalIndex) { + continue; + } + if (pr.sign * ka.sign >= 0) { + continue; + } + float pairMass = 0.f, pairPt = 0.f, pairY = 0.f; + if (!buildProtonKaonPair(pr, ka, pairMass, pairPt, pairY)) { + continue; + } + + fillInvariantMassMCRecoAllPairs(pr.sign, pairPt, pairMass, centrality, pairY); + + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 0.f); + if (pr.hasMcLabel) { + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 1.f); + } + if (ka.hasMcLabel) { + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 2.f); + } + if (!pr.hasMcLabel || !ka.hasMcLabel) { + continue; + } + if (std::abs(pr.mcPdgCode) != lambda1520analysispo_analysis::PdgProton) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 3.f); + if (std::abs(ka.mcPdgCode) != lambda1520analysispo_analysis::PdgKaon) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 4.f); + if (pr.mcLambda1520MotherIndex < 0 || pr.mcLambda1520MotherIndex != ka.mcLambda1520MotherIndex) { + continue; + } + histos.fill(HIST("CutFlow/MCReco/hTruthMatching"), 5.f); + + if (mcCfg.applyGeneratedEtaCut) { + const float pairEta = RecoDecay::eta(std::array{pr.px + ka.px, pr.py + ka.py, pr.pz + ka.pz}); + if (std::abs(pairEta) > mcCfg.generatedEtaMax) { + continue; + } + } + + const bool motherIsParticle = pr.mcLambda1520MotherPdgCode > 0; + fillInvariantMassMCRecoTruthMatched(pr.mcLambda1520MotherPdgCode, pairPt, pairMass, centrality, pairY); + + const float massResolution = pairMass - computeGeneratedMass(mcParticles.iteratorAt(pr.mcLambda1520MotherIndex)); + if (motherIsParticle) { + histos.fill(HIST("InvMass/MCReco/MassResolution/hLambda1520"), massResolution, pairPt, centrality); + } else { + histos.fill(HIST("InvMass/MCReco/MassResolution/hAntiLambda1520"), massResolution, pairPt, centrality); + } + + if (isSpinAlignmentActive()) { + const DecayKinematics decay = makeDecayKinematics(pr, ka); + fillCosThetaStarHistograms(pairMass, pairPt, decay, centrality, motherIsParticle); + + const float cosTheta = computeCosThetaStar(decay); + if (motherIsParticle) { + histos.fill(HIST("CosThetaStar/MCReco/MassResolution/hLambda1520"), massResolution, pairPt, centrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/MCReco/MassResolution/hAntiLambda1520"), massResolution, pairPt, centrality, cosTheta); + } + } + + histos.fill(HIST("PairQA/MCReco/hTruthMatchedRapidityVsPt"), pairPt, pairY); + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Proton/hDCAxyVsPt"), pr.pt, pr.dcaXY); + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Proton/hDCAzVsPt"), pr.pt, pr.dcaZ); + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Proton/hNSigmaTPCvsPt"), pr.pt, pr.tpcNSigmaPr); + if (pr.hasTOF) { + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Proton/hNSigmaTOFvsPt"), pr.pt, pr.tofNSigmaPr); + } + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Kaon/hDCAxyVsPt"), ka.pt, ka.dcaXY); + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Kaon/hDCAzVsPt"), ka.pt, ka.dcaZ); + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Kaon/hNSigmaTPCvsPt"), ka.pt, ka.tpcNSigmaKa); + if (ka.hasTOF) { + histos.fill(HIST("TrackQA/MCReco/TruthMatched/Kaon/hNSigmaTOFvsPt"), ka.pt, ka.tofNSigmaKa); + } + } + } + + if (mcCollisionsWithEfficiencyQA.insert(collision.mcCollisionId()).second) { + fillGeneratedToReconstructedCutFlow(tracksOfCollision, mcParticles, collision.mcCollisionId()); + } + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processMC, "Process MC reconstructed (truth matching)", false); + + void processMCGen(McCollisionsCent::iterator const& mcCollision, + aod::McParticles const& mcParticles, + soa::SmallGroups const& collisions) + { + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 0.f); + + if (std::abs(mcCollision.posZ()) > evSel.cfgCutVertex) { + return; + } + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 1.f); + + if (evSel.cfgApplyMCTruthINELgt0 && !isTruthInelGt0(mcParticles)) { + return; + } + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 2.f); + + float genCentrality = getTruthCentrality(mcCollision, mcParticles); + if (genCentrality == InvalidCentrality) { + return; + } + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 3.f); + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 4.f); + + const AcceptedRecoCollision acceptedReco = findFirstAcceptedRecoCollision(collisions); + const bool hasAcceptedReco = acceptedReco.found; + if (hasAcceptedReco) { + if (!evSel.cfgUseMCTruthCentrality) { + genCentrality = acceptedReco.centrality; + } + histos.fill(HIST("CutFlow/MCGen/hEventSelection"), 5.f); + } + + for (auto const& mcParticle : mcParticles) { + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 0.f); + + if (!mcParticle.producedByGenerator()) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 1.f); + + const int pdg = mcParticle.pdgCode(); + if (std::abs(pdg) != mLambda1520AbsPdgCode) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 2.f); + + auto genProton = mcParticle; + auto genKaon = mcParticle; + if (!findProtonKaonDaughters(mcParticle, genProton, genKaon)) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 3.f); + + const float genPt = mcParticle.pt(); + const float genY = mcParticle.y(); + histos.fill(HIST("MCGenQA/Kinematics/hEtaVsPtBeforeCuts"), mcParticle.eta(), genPt); + histos.fill(HIST("MCGenQA/Kinematics/hPhiVsRapidityBeforeCuts"), mcParticle.phi(), genY); + + if (genY < trackCuts.cfgMinY || genY > trackCuts.cfgMaxY) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 4.f); + + const float maxTrackEta = trackCuts.cfgCutEta; + if (mcCfg.applyGeneratedDaughterEtaCut && (std::abs(genProton.eta()) >= maxTrackEta || std::abs(genKaon.eta()) >= maxTrackEta)) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 5.f); + + if (mcCfg.applyGeneratedEtaCut && std::abs(mcParticle.eta()) > mcCfg.generatedEtaMax) { + continue; + } + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 6.f); + histos.fill(HIST("CutFlow/MCGen/hGeneratedLambda1520Selection"), 7.f); + + histos.fill(HIST("MCGenQA/Kinematics/hEtaVsPtAfterCuts"), mcParticle.eta(), genPt); + histos.fill(HIST("MCGenQA/Kinematics/hPhiVsRapidityAfterCuts"), mcParticle.phi(), genY); + histos.fill(HIST("MCGenQA/Kinematics/hRapidityVsPtAccepted"), genPt, genY); + + const float genMass = computeGeneratedMass(mcParticle); + const float ptPr = genProton.pt(); + const float ptKa = genKaon.pt(); + + if (pdg > 0) { + histos.fill(HIST("InvMass/MCGen/Generated/hLambda1520"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/DaughterPt/hProtonFromLambda1520"), ptPr); + histos.fill(HIST("MCGenQA/DaughterPt/hKaonFromLambda1520"), ptKa); + } else { + histos.fill(HIST("InvMass/MCGen/Generated/hAntiLambda1520"), genMass, genPt, genCentrality, genY); + histos.fill(HIST("MCGenQA/DaughterPt/hProtonFromAntiLambda1520"), ptPr); + histos.fill(HIST("MCGenQA/DaughterPt/hKaonFromAntiLambda1520"), ptKa); + } + + if (hasAcceptedReco) { + if (pdg > 0) { + histos.fill(HIST("InvMass/MCGen/GeneratedInAcceptedEvents/hLambda1520"), genMass, genPt, genCentrality, genY); + } else { + histos.fill(HIST("InvMass/MCGen/GeneratedInAcceptedEvents/hAntiLambda1520"), genMass, genPt, genCentrality, genY); + } + } + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processMCGen, "Process MC truth: generated-level Lambda(1520) spectra and QA", false); + + void processMCGenSpinAlignment(McCollisionsCent::iterator const& mcCollision, + aod::McParticles const& mcParticles, + soa::SmallGroups const& collisions) + { + if (std::abs(mcCollision.posZ()) > evSel.cfgCutVertex) { + return; + } + if (evSel.cfgApplyMCTruthINELgt0 && !isTruthInelGt0(mcParticles)) { + return; + } + float genCentrality = getTruthCentrality(mcCollision, mcParticles); + if (genCentrality == InvalidCentrality) { + return; + } + + const AcceptedRecoCollision acceptedReco = findFirstAcceptedRecoCollision(collisions); + const bool hasAcceptedReco = acceptedReco.found; + if (hasAcceptedReco && !evSel.cfgUseMCTruthCentrality) { + genCentrality = acceptedReco.centrality; + } + + for (auto const& mcParticle : mcParticles) { + if (!mcParticle.producedByGenerator()) { + continue; + } + const int pdg = mcParticle.pdgCode(); + if (std::abs(pdg) != mLambda1520AbsPdgCode) { + continue; + } + auto genProton = mcParticle; + auto genKaon = mcParticle; + if (!findProtonKaonDaughters(mcParticle, genProton, genKaon)) { + continue; + } + + const float genPt = mcParticle.pt(); + const float genY = mcParticle.y(); + if (genY < trackCuts.cfgMinY || genY > trackCuts.cfgMaxY) { + continue; + } + + const float maxTrackEta = trackCuts.cfgCutEta; + if (mcCfg.applyGeneratedDaughterEtaCut && (std::abs(genProton.eta()) >= maxTrackEta || std::abs(genKaon.eta()) >= maxTrackEta)) { + continue; + } + if (mcCfg.applyGeneratedEtaCut && std::abs(mcParticle.eta()) > mcCfg.generatedEtaMax) { + continue; + } + + const float genMass = computeGeneratedMass(mcParticle); + const float cosTheta = computeCosThetaStar(makeGeneratedDecayKinematics(mcParticle, genProton, genKaon)); + + if (pdg > 0) { + histos.fill(HIST("CosThetaStar/MCGen/Generated/hLambda1520"), genMass, genPt, genCentrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/MCGen/Generated/hAntiLambda1520"), genMass, genPt, genCentrality, cosTheta); + } + + if (hasAcceptedReco) { + if (pdg > 0) { + histos.fill(HIST("CosThetaStar/MCGen/GeneratedInAcceptedEvents/hLambda1520"), genMass, genPt, genCentrality, cosTheta); + } else { + histos.fill(HIST("CosThetaStar/MCGen/GeneratedInAcceptedEvents/hAntiLambda1520"), genMass, genPt, genCentrality, cosTheta); + } + } + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processMCGenSpinAlignment, "Process MC truth: generated-level Lambda(1520) spin-alignment cos(theta*) with same cuts as processMCGen", false); + + void processEventFactor(EventCandidatesMC const& collisions, + McCollisionsCent const& mcCollisions, + aod::McParticles const& mcParticles) + { + std::unordered_map nRecoCollisionsPerMc; + std::unordered_map acceptedMcCollisionRecoCent; + for (auto const& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; + } + ++nRecoCollisionsPerMc[collision.mcCollisionId()]; + if (!passesEventSelection(collision)) { + continue; + } + acceptedMcCollisionRecoCent.emplace(collision.mcCollisionId(), getCentrality(collision)); + + const auto& mcCollision = collision.mcCollision_as(); + auto mcParticlesOfRecoCollision = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); + if (evSel.cfgApplyMCTruthINELgt0 && !isTruthInelGt0(mcParticlesOfRecoCollision)) { + continue; + } + const float truthCentrality = getTruthCentrality(mcCollision, mcParticlesOfRecoCollision); + if (truthCentrality == InvalidCentrality) { + continue; + } + histos.fill(HIST("EventLoss/hAcceptedRecoEventsVsTruthCentrality"), truthCentrality); + } + + for (auto const& mcCollision : mcCollisions) { + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 0.f); + + auto mcParticlesOfGenCollision = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); + const bool isInsideVertexRange = std::abs(mcCollision.posZ()) <= evSel.cfgCutVertex; + const bool truthInelGt0 = isTruthInelGt0(mcParticlesOfGenCollision); + const float truthCentrality = getTruthCentrality(mcCollision, mcParticlesOfGenCollision); + + if (!isInsideVertexRange) { + continue; + } + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 1.f); + + if (evSel.cfgApplyMCTruthINELgt0 && !truthInelGt0) { + continue; + } + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 2.f); + + if (truthCentrality == InvalidCentrality) { + continue; + } + histos.fill(HIST("EventLoss/hGeneratedEventsVsTruthCentrality"), truthCentrality); + + const auto recoCollisionCountEntry = nRecoCollisionsPerMc.find(mcCollision.globalIndex()); + const int nRecoTotal = (recoCollisionCountEntry != nRecoCollisionsPerMc.end()) ? recoCollisionCountEntry->second : 0; + histos.fill(HIST("EventLoss/hRecoCollisionMultiplicity"), nRecoTotal); + + if (nRecoTotal > 0) { + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 3.f); + } + + const auto acceptedRecoCentralityEntry = acceptedMcCollisionRecoCent.find(mcCollision.globalIndex()); + if (acceptedRecoCentralityEntry != acceptedMcCollisionRecoCent.end()) { + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 4.f); + histos.fill(HIST("CutFlow/EventLoss/hGeneratedEventSelection"), 5.f); + histos.fill(HIST("EventLoss/hGeneratedAcceptedEventsVsTruthCentrality"), truthCentrality); + histos.fill(HIST("EventLoss/hGeneratedAcceptedEventsVsRecoCentrality"), acceptedRecoCentralityEntry->second); + } + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processEventFactor, "Process Event Loss / Event Splitting factors", false); + + void processSignalLoss(McCollisionsCent::iterator const& mcCollision, + aod::McParticles const& mcParticles, + soa::SmallGroups const& collisions) + { + if (std::abs(mcCollision.posZ()) > evSel.cfgCutVertex) { + return; + } + if (evSel.cfgApplyMCTruthINELgt0 && !isTruthInelGt0(mcParticles)) { + return; + } + const float centrality = getTruthCentrality(mcCollision, mcParticles); + if (centrality == InvalidCentrality) { + return; + } + + const bool hasAcceptedReco = findFirstAcceptedRecoCollision(collisions).found; + + const float minRapidity = trackCuts.cfgMinY; + const float maxRapidity = trackCuts.cfgMaxY; + const float maxAbsGeneratedEta = mcCfg.generatedEtaMax; + + for (auto const& part : mcParticles) { + const int absPdgCode = std::abs(part.pdgCode()); + const bool isLambda1520 = absPdgCode == mLambda1520AbsPdgCode; + const bool isProxySpecies = absPdgCode == lambda1520analysispo_analysis::PdgKaon || + absPdgCode == lambda1520analysispo_analysis::PdgProton || + absPdgCode == lambda1520analysispo_analysis::PdgLambda0 || + absPdgCode == lambda1520analysispo_analysis::PdgXiMinus || + absPdgCode == lambda1520analysispo_analysis::PdgOmegaMinus; + if (!isLambda1520 && !isProxySpecies) { + continue; + } + + if (part.y() < minRapidity || part.y() > maxRapidity) { + continue; + } + if (mcCfg.applyGeneratedEtaCut && std::abs(part.eta()) > maxAbsGeneratedEta) { + continue; + } + + if (isLambda1520) { + if (!part.producedByGenerator()) { + continue; + } + auto genProton = part; + auto genKaon = part; + if (!findProtonKaonDaughters(part, genProton, genKaon)) { + continue; + } + if (part.pdgCode() > 0) { + histos.fill(HIST("SignalLoss/Lambda1520/hGeneratedPtVsCentrality"), part.pt(), centrality); + if (hasAcceptedReco) { + histos.fill(HIST("SignalLoss/Lambda1520/hAcceptedEventPtVsCentrality"), part.pt(), centrality); + } + } else { + histos.fill(HIST("SignalLoss/AntiLambda1520/hGeneratedPtVsCentrality"), part.pt(), centrality); + if (hasAcceptedReco) { + histos.fill(HIST("SignalLoss/AntiLambda1520/hAcceptedEventPtVsCentrality"), part.pt(), centrality); + } + } + continue; + } + + if (!part.isPhysicalPrimary()) { + continue; + } + fillSignalLossProxySpecies(part, centrality); + if (hasAcceptedReco) { + fillSignalLossProxySpecies(part, centrality); + } + } + } + PROCESS_SWITCH(Lambda1520Analysispo, processSignalLoss, "Process signal loss (K, p, Lambda, Xi, Omega proxies and Lambda(1520))", false); }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(context)}; }