diff --git a/PWGLF/Tasks/Resonances/kstar0analysis.cxx b/PWGLF/Tasks/Resonances/kstar0analysis.cxx index 4c3a32afcae..ba0e75d4e1a 100644 --- a/PWGLF/Tasks/Resonances/kstar0analysis.cxx +++ b/PWGLF/Tasks/Resonances/kstar0analysis.cxx @@ -8,11 +8,16 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. + /// \file kstar0analysis.cxx /// \brief the reconstruction of k*0(892) resonance analysis using K\pi decay channel /// \author Hirak Kumar Koley +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + #include "Common/CCDB/EventSelectionParams.h" +#include "Common/CCDB/RCTSelectionFlags.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" @@ -20,6 +25,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include #include #include #include @@ -28,918 +34,1312 @@ #include #include #include +#include #include #include #include #include +#include #include +#include #include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include -#include +#include +#include #include #include -#include -#include -#include #include -#include #include -#include - using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::aod::rctsel; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::constants::physics; -struct kstar0analysis { +using LorentzVectorPtEtaPhiMass = ROOT::Math::PtEtaPhiMVector; + +enum TrackSelectionType { + AllTracks = 0, + GlobalTracks, + GlobalTracksWoPtEta, + GlobalTracksWoDCA, + QualityTracks, + InAcceptanceTracks, +}; + +enum PIDCutType { + SquareType = 1, + CircularType, +}; + +struct Kstar0analysis { + // Define slice per collision + Preslice perCollision = o2::aod::track::collisionId; SliceCache cache; - Preslice perCollision = aod::track::collisionId; + Preslice perMcCollision = o2::aod::mcparticle::mcCollisionId; + SliceCache cacheMC; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - //================================== - //|| - //|| Selection - //|| - //================================== - - // Event Selection - Configurable cfgEventSelections{"cfgEventSelections", "sel8", "Set event selection"}; - Configurable cfgEventVtxCut{"cfgEventVtxCut", 10.0, "V_z cut selection"}; - ConfigurableAxis cfgCentAxis{"cfgCentAxis", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0}, "Binning of the centrality axis"}; - Configurable cfgOccupancySel{"cfgOccupancySel", false, "Occupancy selection"}; - Configurable cfgOccupancyMax{"cfgOccupancyMax", 999999., "maximum occupancy of tracks in neighbouring collisions in a given time range"}; - Configurable cfgOccupancyMin{"cfgOccupancyMin", -100., "minimum occupancy of tracks in neighbouring collisions in a given time range"}; - - // Track Selection - // General - Configurable cfgTrackMinPt{"cfgTrackMinPt", 0.15, "set track min pT"}; - Configurable cfgTrackMaxEta{"cfgTrackMaxEta", 0.8, "set track max Eta"}; - Configurable cfgTrackMaxDCArToPVcut{"cfgTrackMaxDCArToPVcut", 0.5, "Track DCAr cut to PV Maximum"}; - Configurable cfgTrackMaxDCAzToPVcut{"cfgTrackMaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; - Configurable cfgTrackGlobalSel{"cfgTrackGlobalSel", true, "Global track selection"}; - Configurable cfgTrackPrimaryTrack{"cfgTrackPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz - Configurable cfgTrackConnectedToPV{"cfgTrackConnectedToPV", true, "PV contributor track selection"}; // PV Contriuibutor - Configurable cfgTrackGlobalWoDCATrack{"cfgTrackGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) - - // TPC - Configurable cfgTrackFindableTPCClusters{"cfgTrackFindableTPCClusters", 50, "nFindable TPC Clusters"}; - Configurable cfgTrackTPCCrossedRows{"cfgTrackTPCCrossedRows", 70, "nCrossed TPC Rows"}; - Configurable cfgTrackRowsOverFindable{"cfgTrackRowsOverFindable", 1.2, "nRowsOverFindable TPC CLusters"}; - Configurable cfgTrackTPCChi2{"cfgTrackTPCChi2", 4.0, "nTPC Chi2 per Cluster"}; - - // ITS - Configurable cfgTrackITSChi2{"cfgTrackITSChi2", 36.0, "nITS Chi2 per Cluster"}; - - // PID - Configurable cfgTrackTPCPID{"cfgTrackTPCPID", true, "Enables TPC PID"}; - Configurable cfgTrackTOFPID{"cfgTrackTOFPID", true, "Enables TOF PID"}; - Configurable cfgTrackSquarePIDCut{"cfgTrackSqurePIDCut", true, "Enables PID cut shape square switch"}; - Configurable cfgTrackCirclePIDCut{"cfgTrackCirclePIDCut", true, "Enables PID cut shape circle switch"}; - Configurable cfgTrackCircleValue{"cfgTrackCircleValue", 2, "Enables TOF TPC PID circle cut value"}; - Configurable cfgTrackTOFHard{"cfgTrackTOFHard", false, "Enables TOF Hard"}; - - Configurable cfgTrackTPCPIDnSig{"cfgTrackTPCPIDnSig", 4.0, "nTPC PID sigma"}; - Configurable cfgTrackTOFPIDnSig{"cfgTrackTOFPIDnSig", 4.0, "nTOF PID sigma"}; - Configurable cDebugLevel{"cDebugLevel", 0, "Resolution of Debug"}; - - // Mixing - ConfigurableAxis cfgBinsMixMult{"cfgBinsCent", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0}, "Binning of the centrality axis"}; - ConfigurableAxis cfgBinsMixVtx{"cfgBinsMixVtx", {VARIABLE_WIDTH, -10.0f, -5.f, 0.f, 5.f, 10.f}, "Mixing bins - z-vertex"}; - Configurable cfgMixNMixedEvents{"cfgMixNMixedEvents", 10, "Number of mixed events per event"}; - Configurable cfgVtxMixCut{"cfgVtxMixCut", 10, "Vertex Mix Cut"}; - - // MCGen - Configurable cfgForceGenReco{"cfgForceGenReco", false, "Only consider events which are reconstructed (neglect event-loss)"}; - Configurable cfgReqMcEffPID{"cfgReqMcEffPID", false, "Request McEfficiency PID"}; - - // Pair - Configurable cfgMinvNBins{"cfgMinvNBins", 300, "Number of bins for Minv axis"}; - Configurable cfgMinvMin{"cfgMinvMin", 0.60, "Minimum Minv value"}; - Configurable cfgMinvMax{"cfgMinvMax", 1.20, "Maximum Minv value"}; - - // Histogram - ConfigurableAxis binsDCAz{"binsDCAz", {40, -0.2, 0.2}, ""}; - ConfigurableAxis binsDCAxy{"binsDCAxy", {40, -0.2, 0.2}, ""}; - Configurable cfgEventCutQA{"cfgEventCutsQA", false, "Enable Event QA Hists"}; - Configurable cfgTrackCutQA{"cfgTrackCutQA", false, "Enable Track QA Hists"}; - - Configurable cfgMCHistos{"cfgMCHistos", false, "Enable MC Hists"}; - Configurable cfgMixedHistos{"cfgMixedHistos", false, "Enable Mixed Histos"}; - Configurable cfgManualEvSel{"cfgManualEvSel", false, "Enable Manual EvSel"}; - - // Main + + Service pdg; + RCTFlagsChecker rctChecker; + + struct : ConfigurableGroup { + Configurable cfgEvtZvtx{"cfgEvtZvtx", 10.0f, "Evt sel: Max. z-Vertex (cm)"}; + Configurable cfgEvtTriggerTVXSel{"cfgEvtTriggerTVXSel", true, "Evt sel: triggerTVX selection (MB)"}; + Configurable cfgEvtNoTFBorderCut{"cfgEvtNoTFBorderCut", true, "Evt sel: apply TF border cut"}; + Configurable cfgEvtNoITSROFrameBorderCut{"cfgEvtNoITSROFrameBorderCut", false, "Evt sel: apply NoITSRO border cut"}; + Configurable cfgEvtIsRCTFlagpassed{"cfgEvtIsRCTFlagpassed", false, "Evt sel: apply RCT flag selection"}; + Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; + Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; + Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", true, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; + Configurable cfgEvtSel8{"cfgEvtSel8", false, "Evt Sel 8 check for offline selection"}; + Configurable cfgEvtIsINELgt0{"cfgEvtIsINELgt0", false, "Evt sel: apply INEL>0 selection"}; + } configEvents; + + struct : ConfigurableGroup { + // Pre-selection Track cuts + Configurable trackSelection{"trackSelection", 0, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; + Configurable cMinPtcut{"cMinPtcut", 0.15f, "Minimal pT for tracks"}; + Configurable cMinTPCNClsFound{"cMinTPCNClsFound", 120, "minimum TPCNClsFound value for good track"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable cfgCutRapidity{"cfgCutRapidity", 0.5f, "rapidity range for particles"}; + Configurable cfgMinCrossedRows{"cfgMinCrossedRows", 70, "min crossed rows for good track"}; + + // DCA Selections + // DCAr to PV + Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.1f, "Track DCAr cut to PV Maximum"}; + // DCAz to PV + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 0.1f, "Track DCAz cut to PV Maximum"}; + + // Track selections + Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor + Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; + Configurable cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; + Configurable cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"}; + Configurable cTPCNClsFound{"cTPCNClsFound", false, "Switch to turn on/off TPCNClsFound cut"}; + Configurable cDCAr7SigCut{"cDCAr7SigCut", false, "Track DCAr 7 Sigma cut to PV Maximum"}; + } configTracks; + + struct : ConfigurableGroup { + /// PID Selections + Configurable pidnSigmaPreSelectionCut{"pidnSigmaPreSelectionCut", 4.0f, "pidnSigma Cut for pre-selection of tracks"}; + Configurable cByPassTOF{"cByPassTOF", false, "By pass TOF PID selection"}; // By pass TOF PID selection + Configurable cPIDcutType{"cPIDcutType", 2, "cPIDcutType = 1 for square cut, 2 for circular cut"}; // By pass TOF PID selection + + // Kaon + Configurable> kaonTPCPIDpTintv{"kaonTPCPIDpTintv", {0.5f}, "pT intervals for Kaon TPC PID cuts"}; + Configurable> kaonTPCPIDcuts{"kaonTPCPIDcuts", {2}, "nSigma list for Kaon TPC PID cuts"}; + Configurable> kaonTOFPIDpTintv{"kaonTOFPIDpTintv", {999.0f}, "pT intervals for Kaon TOF PID cuts"}; + Configurable> kaonTOFPIDcuts{"kaonTOFPIDcuts", {2}, "nSigma list for Kaon TOF PID cuts"}; + Configurable> kaonTPCTOFCombinedpTintv{"kaonTPCTOFCombinedpTintv", {999.0f}, "pT intervals for Kaon TPC-TOF PID cuts"}; + Configurable> kaonTPCTOFCombinedPIDcuts{"kaonTPCTOFCombinedPIDcuts", {2}, "nSigma list for Kaon TPC-TOF PID cuts"}; + + // Pion + Configurable> pionTPCPIDpTintv{"pionTPCPIDpTintv", {0.9f}, "pT intervals for Pion TPC PID cuts"}; + Configurable> pionTPCPIDcuts{"pionTPCPIDcuts", {2}, "nSigma list for Pion TPC PID cuts"}; + Configurable> pionTOFPIDpTintv{"pionTOFPIDpTintv", {999.0f}, "pT intervals for Pion TOF PID cuts"}; + Configurable> pionTOFPIDcuts{"pionTOFPIDcuts", {2}, "nSigma list for Pion TOF PID cuts"}; + Configurable> pionTPCTOFCombinedpTintv{"pionTPCTOFCombinedpTintv", {999.0f}, "pT intervals for Pion TPC-TOF PID cuts"}; + Configurable> pionTPCTOFCombinedPIDcuts{"pionTPCTOFCombinedPIDcuts", {2}, "nSigma list for Pion TPC-TOF PID cuts"}; + } configPID; + + struct : ConfigurableGroup { + /// Event Mixing + Configurable nEvtMixing{"nEvtMixing", 10, "Number of events to mix"}; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -10.0f, -8.0f, -6.0f, -4.0f, -2.0f, 0.0f, 2.0f, 4.0f, 6.0f, 8.0f, 10.0f}, "Mixing bins - z-vertex"}; + ConfigurableAxis cfgMultPercentileBins{"cfgMultPercentileBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "Mixing bins - multiplicity"}; + + // Rotational background + Configurable rotationalcut{"rotationalcut", 10, "Cut value (Rotation angle pi - pi/cut and pi + pi/cut)"}; + Configurable cNofRotations{"cNofRotations", 3, "Number of random rotations in the rotational background"}; + Configurable cfgRotPi{"cfgRotPi", true, "rotate Pion"}; + } configBkg; + + // Additional purity check + Configurable crejectProton{"crejectProton", false, "Switch to turn on/off proton contamination"}; + Configurable cUseOpeningAngleCut{"cUseOpeningAngleCut", false, "Kinematic Cuts for p-K pair opening angle"}; + Configurable cMinOpeningAngle{"cMinOpeningAngle", 1.4f, "Minimum opening angle between daughters"}; + Configurable cMaxOpeningAngle{"cMaxOpeningAngle", 2.4f, "Maximum opening angle between daughters"}; + Configurable cfgUseDeltaEtaPhiCuts{"cfgUseDeltaEtaPhiCuts", false, "Switch to turn on/off delta eta and delta phi cuts"}; + Configurable cfgUseDaughterEtaCutMC{"cfgUseDaughterEtaCutMC", false, "Switch to turn on/off eta cuts for daughters in MC"}; + + // MC selection cut + Configurable cUseRapcutMC{"cUseRapcutMC", true, "Use rapidity cut for MC"}; + + // cuts on mother + Configurable cfgUseCutsOnMother{"cfgUseCutsOnMother", false, "Enable additional cuts on mother"}; + Configurable cMaxPtMotherCut{"cMaxPtMotherCut", 10.0f, "Maximum pt of mother cut"}; + Configurable cMaxMinvMotherCut{"cMaxMinvMotherCut", 3.0f, "Maximum Minv of mother cut"}; + Configurable cMaxDeltaEtaCut{"cMaxDeltaEtaCut", 0.7f, "Maximum deltaEta between daughters"}; + Configurable cMaxDeltaPhiCut{"cMaxDeltaPhiCut", 1.5f, "Maximum deltaPhi between daughters"}; + + // switches + Configurable cFillMultQA{"cFillMultQA", false, "Turn on/off additional QA plots"}; + Configurable cFillTrackQA{"cFillTrackQA", false, "Turn on/off additional QA plots"}; + Configurable cFilladditionalQAeventPlots{"cFilladditionalQAeventPlots", false, "Additional QA event plots"}; + Configurable cFilladditionalMEPlots{"cFilladditionalMEPlots", false, "Additional Mixed event plots"}; + Configurable cFilldeltaEtaPhiPlots{"cFilldeltaEtaPhiPlots", false, "Enable additional cuts on daughters"}; + Configurable cFill1DQAs{"cFill1DQAs", false, "Invariant mass 1D"}; + Configurable centEstimator{"centEstimator", 0, "Select centrality estimator: 0 - FT0M, 1 - FT0A, 2 - FT0C"}; + + TRandom* rn = new TRandom(); + + // Pre-filters for efficient process + Filter zVtxFilter = (nabs(o2::aod::collision::posZ) <= configEvents.cfgEvtZvtx); + Filter collisionFilterMC = nabs(aod::mccollision::posZ) <= configEvents.cfgEvtZvtx; + + Filter acceptanceFilter = (nabs(aod::track::eta) < configTracks.cfgCutEta && nabs(aod::track::pt) > configTracks.cMinPtcut) && + (nabs(aod::track::dcaXY) < configTracks.cMaxDCArToPVcut) && (nabs(aod::track::dcaZ) < configTracks.cMaxDCAzToPVcut); + + using EventCandidates = soa::Join; + using TrackCandidates = soa::Filtered>; + + // for MC reco + using MCEventCandidates = soa::Join; + using MCTrackCandidates = soa::Filtered>; + + /// Figures + ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.12f, 0.14f, 0.16f, 0.18f, 0.2f, 0.25f, 0.3f, 0.35f, 0.4f, 0.45f, 0.5f, 0.55f, 0.6f, 0.65f, 0.7f, 0.75f, 0.8f, 0.85f, 0.9f, 0.95f, 1.0f, 1.1f, 1.2f, 1.25f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.75f, 1.8f, 1.9f, 2.0f, 2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f, 3.0f, 3.1f, 3.2f, 3.3f, 3.4f, 3.6f, 3.7f, 3.8f, 3.9f, 4.0f, 4.1f, 4.2f, 4.5f, 4.6f, 4.8f, 4.9f, 5.0f, 5.5f, 5.6f, 6.0f, 6.4f, 6.5f, 7.0f, 7.2f, 8.0f, 9.0f, 9.5f, 9.6f, 10.0f, 11.0f, 11.5f, 12.0f, 13.0f, 14.0f, 14.4f, 15.0f, 16.0f, 18.0f, 19.2f, 20.0f}, "Binning of the pT axis"}; + ConfigurableAxis binsPtQA{"binsPtQA", {VARIABLE_WIDTH, 0.0f, 0.2f, 0.4f, 0.6f, 0.8f, 1.0f, 1.2f, 1.4f, 1.6f, 1.8f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.2f, 4.4f, 4.6f, 4.8f, 5.0f, 5.2f, 5.4f, 5.6f, 5.8f, 6.0f, 6.2f, 6.4f, 6.6f, 6.8f, 7.0f, 7.2f, 7.4f, 7.6f, 7.8f, 8.0f, 8.2f, 8.4f, 8.6f, 8.8f, 9.0f, 9.2f, 9.4f, 9.6f, 9.8f, 10.0f, 10.2f, 10.4f, 10.6f, 10.8f, 11.0f, 11.2f, 11.4f, 11.6f, 11.8f, 12.0f, 12.2f, 12.4f, 12.6f, 12.8f, 13.0f, 13.2f, 13.4f, 13.6f, 13.8f, 14.0f, 14.2f, 14.4f, 14.6f, 14.8f, 15.0f, 15.2f, 15.4f, 15.6f, 15.8f, 16.0f, 16.2f, 16.4f, 16.6f, 16.8f, 17.0f, 17.2f, 17.4f, 17.6f, 17.8f, 18.0f, 18.2f, 18.4f, 18.6f, 18.8f, 19.0f, 19.2f, 19.4f, 19.6f, 19.8f, 20.0f}, "Binning of the pT axis"}; + ConfigurableAxis binsEta{"binsEta", {150, -1.5f, 1.5f}, ""}; + ConfigurableAxis binsMass{"binsMass", {300, 0.6f, 1.2f}, "Invariant Mass (GeV/#it{c}^2)"}; + ConfigurableAxis binsMult{"binsMult", {105, 0.0f, 105.0f}, "mult_{FT0M}"}; + ConfigurableAxis binsDCAz{"binsDCAz", {40, -0.2f, 0.2f}, ""}; + ConfigurableAxis binsDCAxy{"binsDCAxy", {40, -0.2f, 0.2f}, ""}; + ConfigurableAxis binsTPCXrows{"binsTPCXrows", {100, 60, 160}, ""}; + ConfigurableAxis binsnSigma{"binsnSigma", {130, -6.5f, 6.5f}, ""}; + ConfigurableAxis binsnTPCSignal{"binsnTPCSignal", {1000, 0, 1000}, ""}; + ConfigurableAxis binsEtaPhi{"binsEtaPhi", {350, -3.5f, 3.5f}, ""}; + void init(o2::framework::InitContext&) { - // HISTOGRAMS - const AxisSpec axisEta{30, -1.5, +1.5, "#eta"}; - const AxisSpec axisPhi{200, -1, +7, "#phi"}; - const AxisSpec ptAxis = {200, 0, 20.0}; - const AxisSpec pidAxis = {120, -6, 6}; - const AxisSpec minvAxis = {cfgMinvNBins, cfgMinvMin, cfgMinvMax}; - const AxisSpec axisDCAz{binsDCAz, "DCA_{z}"}; - const AxisSpec axisDCAxy{binsDCAxy, "DCA_{XY}"}; - - if (cfgEventCutQA) { - histos.add("hEvent_Cut", "Number of event after cuts", kTH1D, {{13, -0.5, 12.5}}); - histos.add("hPosZ_BC", "hPosZ_BC", kTH1F, {{300, -15.0, 15.0}}); - histos.add("hPosZ_AC", "hPosZ_AC", kTH1F, {{300, -15.0, 15.0}}); - histos.add("hcentFT0C_BC", "centFT0C_BC", kTH1F, {{110, 0.0, 110.0}}); - histos.add("hcentFT0C_AC", "centFT0C_AC", kTH1F, {{110, 0.0, 110.0}}); - } - if (cfgTrackCutQA) { - histos.add("hDCArToPv_BC", "DCArToPv_BC", kTH1F, {axisDCAxy}); - histos.add("hDCAzToPv_BC", "DCAzToPv_BC", kTH1F, {axisDCAz}); - histos.add("hIsPrim_BC", "hIsPrim_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsGood_BC", "hIsGood_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsPrimCont_BC", "hIsPrimCont_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hFindableTPCClusters_BC", "hFindableTPCClusters_BC", kTH1F, {{200, 0, 200}}); - histos.add("hFindableTPCRows_BC", "hFindableTPCRows_BC", kTH1F, {{200, 0, 200}}); - histos.add("hClustersVsRows_BC", "hClustersVsRows_BC", kTH1F, {{200, 0, 2}}); - histos.add("hTPCChi2_BC", "hTPCChi2_BC", kTH1F, {{200, 0, 100}}); - histos.add("hITSChi2_BC", "hITSChi2_BC", kTH1F, {{200, 0, 100}}); - histos.add("QA_nSigma_pion_TPC_BC", "QA_nSigma_pion_TPC_BC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_nSigma_pion_TOF_BC", "QA_nSigma_pion_TOF_BC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_pion_TPC_TOF_BC", "QA_pion_TPC_TOF_BC", {HistType::kTH2F, {pidAxis, pidAxis}}); - histos.add("QA_nSigma_kaon_TPC_BC", "QA_nSigma_kaon_TPC_BC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_nSigma_kaon_TOF_BC", "QA_nSigma_kaon_TOF_BC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_kaon_TPC_TOF_BC", "QA_kaon_TPC_TOF_BC", {HistType::kTH2F, {pidAxis, pidAxis}}); - histos.add("QA_track_pT_BC", "QA_track_pT_BC", kTH1F, {{13, 0.0, 13.0}}); - - histos.add("hDCArToPv_AC", "DCArToPv_AC", kTH1F, {axisDCAxy}); - histos.add("hDCAzToPv_AC", "DCAzToPv_AC", kTH1F, {axisDCAz}); - histos.add("hIsPrim_AC", "hIsPrim_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsGood_AC", "hIsGood_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsPrimCont_AC", "hIsPrimCont_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hFindableTPCClusters_AC", "hFindableTPCClusters_AC", kTH1F, {{200, 0, 200}}); - histos.add("hFindableTPCRows_AC", "hFindableTPCRows_AC", kTH1F, {{200, 0, 200}}); - histos.add("hClustersVsRows_AC", "hClustersVsRows_AC", kTH1F, {{200, 0, 2}}); - histos.add("hTPCChi2_AC", "hTPCChi2_AC", kTH1F, {{200, 0, 100}}); - histos.add("hITSChi2_AC", "hITSChi2_AC", kTH1F, {{200, 0, 100}}); - histos.add("QA_nSigma_pion_TPC_AC", "QA_nSigma_pion_TPC_AC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_nSigma_pion_TOF_AC", "QA_nSigma_pion_TOF_AC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_pion_TPC_TOF_AC", "QA_pion_TPC_TOF_AC", {HistType::kTH2F, {pidAxis, pidAxis}}); - histos.add("QA_nSigma_kaon_TPC_AC", "QA_nSigma_kaon_TPC_AC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_nSigma_kaon_TOF_AC", "QA_nSigma_kaon_TOF_AC", {HistType::kTH2F, {ptAxis, pidAxis}}); - histos.add("QA_kaon_TPC_TOF_AC", "QA_kaon_TPC_TOF_AC", {HistType::kTH2F, {pidAxis, pidAxis}}); - histos.add("QA_track_pT_AC", "QA_track_pT_AC", kTH1F, {{13, 0.0, 13.0}}); + rctChecker.init(configEvents.cfgEvtRCTFlagCheckerLabel, configEvents.cfgEvtRCTFlagCheckerZDCCheck, configEvents.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + + // axes + AxisSpec axisPt{binsPt, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisPtQA{binsPtQA, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisEta{binsEta, "#eta"}; + AxisSpec axisRap{binsEta, "#it{y}"}; + AxisSpec axisMassk892{binsMass, "Invariant Mass (GeV/#it{c}^2)"}; + AxisSpec axisMult{binsMult, "mult_{FT0M}"}; + AxisSpec axisDCAz{binsDCAz, "DCA_{z}"}; + AxisSpec axisDCAxy{binsDCAxy, "DCA_{XY}"}; + AxisSpec axisTPCXrow{binsTPCXrows, "#Xrows_{TPC}"}; + AxisSpec axisPIDQA{binsnSigma, "#sigma"}; + AxisSpec axisTPCSignal{binsnTPCSignal, ""}; + AxisSpec axisEtaPhi{binsEtaPhi, ""}; + AxisSpec axisPhi{350, 0, 7, "#Phi"}; + AxisSpec axisMultMix{configBkg.cfgMultPercentileBins, "Multiplicity Percentile"}; + AxisSpec axisVtxMix{configBkg.cfgVtxBins, "Vertex Z (cm)"}; + + histos.add("CollCutCounts", "No. of event after cuts", kTH1I, {{10, 0, 10}}); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(1, "All Events"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(2, "|Vz| < cut"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(3, "kIsTriggerTVX"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(5, "kNoITSROFrameBorder"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(6, "rctChecker"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(7, "sel8"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(8, "IsINELgt0"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(9, "All Passed Events"); + + histos.add("Event/posZ", "; vtx_{z} (cm); Entries", HistType::kTH1F, {{250, -12.5, 12.5}}); + histos.add("Event/centFT0M", "; FT0M Percentile; Entries", HistType::kTH1F, {{110, 0, 110}}); + + if (cFilladditionalQAeventPlots) { + // event histograms + if (doprocessData) { + histos.add("QAevent/hPairsCounterSameE", "total valid no. of pairs sameE", HistType::kTH1F, {{1, 0.5f, 1.5f}}); + histos.add("QAevent/hnTrksSameE", "n tracks per event SameE", HistType::kTH1F, {{1000, 0.0, 1000.0}}); + } + // Gen on Mixed event + if (doprocessME) { + + // Histograms for Mixed Event Pool characteristics + histos.add("QAevent/hMixPool_VtxZ", "Mixed Event Pool: Vertex Z;Vtx Z (cm);Counts", HistType::kTH1F, {axisVtxMix}); + histos.add("QAevent/hMixPool_Multiplicity", "Mixed Event Pool: Multiplicity;Multiplicity;Counts", HistType::kTH1F, {axisMultMix}); + histos.add("QAevent/hMixPool_VtxZ_vs_Multiplicity", "Mixed Event Pool: Vertex Z vs Multiplicity;Counts", HistType::kTH2F, {axisVtxMix, axisMultMix}); + + histos.add("QAevent/hPairsCounterMixedE", "total valid no. of pairs mixedE", HistType::kTH1F, {{1, 0.5f, 1.5f}}); + histos.add("QAevent/hVertexZMixedE", "Collision Vertex Z position", HistType::kTH1F, {{100, -15.0f, 15.0f}}); + histos.add("QAevent/hMultiplicityPercentMixedE", "Multiplicity percentile of collision", HistType::kTH1F, {{120, 0.0f, 120.0f}}); + histos.add("QAevent/hnTrksMixedE", "n tracks per event MixedE", HistType::kTH1F, {{1000, 0.0f, 1000.0f}}); + } } - //////////////////////////////////// - histos.add("nEvents", "nEvents", kTH1F, {{7, 0.0, 7.0}}); - histos.add("hUSS_KPi", "hUSS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hUSS_PiK", "hUSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hLSS_KPi", "hLSS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hLSS_PiK", "hLSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + if (doprocessData) { + // Track QA before cuts + // --- Track + if (cFillTrackQA) { + histos.add("QA/QAbefore/Track/TOF_TPC_Map_ka_all", "TOF + TPC Combined PID for all Kaon;{#sigma_{TOF}^{Kaon}};{#sigma_{TPC}^{Kaon}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TOF_Nsigma_ka_all", "TOF NSigma for all Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPC_Nsigma_ka_all", "TPC NSigma for all Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPConly_Nsigma_ka_all", "TPC NSigma for all Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + + histos.add("QA/QAbefore/Track/TOF_TPC_Map_pi_all", "TOF + TPC Combined PID for all Pion;{#sigma_{TOF}^{Pion}};{#sigma_{TPC}^{Pion}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TOF_Nsigma_pi_all", "TOF NSigma for all Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPC_Nsigma_pi_all", "TPC NSigma for all Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPConly_Nsigma_pi_all", "TPC NSigma for all Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Pion}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + + histos.add("QA/QAbefore/Track/dcaZ_all", "DCA_{Z} distribution of all Tracks; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm); ", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAbefore/Track/dcaXY_all", "DCA_{XY} momentum distribution of all Tracks; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAbefore/Track/TPC_CR_all", "# TPC Xrows distribution of all Tracks; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAbefore/Track/pT_all", "pT distribution of all Tracks; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAbefore/Track/eta_all", "#eta distribution of all Tracks; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + } + if (cFillMultQA) { + // Multiplicity correlation calibrations + histos.add("MultCalib/centGloPVpi", "Centrality vs Global-Tracks", kTHnSparseF, {{110, 0, 110, "Centrality"}, {500, 0, 5000, "Global Tracks"}, {500, 0, 5000, "PV tracks"}}); + histos.add("MultCalib/centGloPVka", "Centrality vs Global-Tracks", kTHnSparseF, {{110, 0, 110, "Centrality"}, {500, 0, 5000, "Global Tracks"}, {500, 0, 5000, "PV tracks"}}); + } - if (cfgMixedHistos) { - histos.add("hUSS_KPi_Mix", "hUSS_KPi_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hUSS_PiK_Mix", "hUSS_PiK_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - } - if (cfgMCHistos) { - histos.add("nEvents_Gen", "nEvents_Gen", kTH1F, {{4, 0.0, 4.0}}); - histos.add("hUSS_TrueRec", "hUSS_TrueRec", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hGen_pT_Raw", "Gen_pT_Raw (GeV/c)", kTH1F, {{800, 0., 40.}}); - histos.add("hGen_pT_GoodEv", "hGen_pT_GoodEv", kTHnSparseF, {cfgCentAxis, ptAxis}); - } + // PID QA after cuts + // --- Kaon + histos.add("QA/QAafter/Kaon/TOF_TPC_Map_ka_selected", "TOF + TPC Combined PID for selected Kaon;{#sigma_{TOF}^{Kaon}};{#sigma_{TPC}^{Kaon}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAafter/Kaon/TOF_Nsigma_ka_selected", "TOF NSigma for selected Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Kaon/TPC_Nsigma_ka_selected", "TPC NSigma for selected Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Kaon/TPC_Nsigma_ka_TPConly_selected", "TPC NSigma for selected Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Kaon/dcaZ_selected", "DCA_{Z} distribution of selected Kaons; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm); ", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAafter/Kaon/dcaXY_selected", "DCA_{XY} momentum distribution of selected Kaons; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAafter/Kaon/TPC_CR_selected", "# TPC Xrows distribution of selected Kaons; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAafter/Kaon/pT_selected", "pT distribution of selected Kaons; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAafter/Kaon/eta_selected", "#eta distribution of selected Kaons; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + histos.add("QA/QAafter/Kaon/TPC_Signal_ka_selected", "TPC Signal for selected Kaon;#it{p} (GeV/#it{c});TPC Signal (A.U.)", {HistType::kTH2F, {axisPt, axisTPCSignal}}); + histos.add("QA/QAafter/Kaon/TPCnclusterPhika_selected", "TPC ncluster vs phi for selected Kaon", kTHnSparseF, {{160, 0, 160, "TPC nCluster"}, {63, 0.0f, 6.28f, "#phi"}}); + + // --- Pion + histos.add("QA/QAafter/Pion/TOF_TPC_Map_pi_selected", "TOF + TPC Combined PID for selected Pion;{#sigma_{TOF}^{Pion}};{#sigma_{TPC}^{Pion}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAafter/Pion/TOF_Nsigma_pi_selected", "TOF NSigma for selected Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Pion/TPC_Nsigma_pi_selected", "TPC NSigma for selected Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Pion/TPC_Nsigma_pi_TPConly_selected", "TPC NSigma for selected Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Pion}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Pion/dcaZ_selected", "DCA_{Z} distribution of selected Pions; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm);", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAafter/Pion/dcaXY_selected", "DCA_{XY} momentum distribution of selected Pions; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAafter/Pion/TPC_CR_selected", "# TPC Xrows distribution of selected Pions; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAafter/Pion/pT_selected", "pT distribution of selected Pions; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAafter/Pion/eta_selected", "#eta distribution of selected Pions; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + histos.add("QA/QAafter/Pion/TPC_Signal_pi_selected", "TPC Signal for selected Pion;#it{p} (GeV/#it{c});TPC Signal (A.U.)", {HistType::kTH2F, {axisPt, axisTPCSignal}}); + histos.add("QA/QAafter/Pion/TPCnclusterPhipi_selected", "TPC ncluster vs phi for selected Pion", kTHnSparseF, {{160, 0, 160, "TPC nCluster"}, {63, 0.0f, 6.28f, "#phi"}}); + + // Mass QA 1D for quick check + if (cFill1DQAs) { + histos.add("Result/Data/k892invmass", "Invariant mass of K(892) K^{#pm}#pi^{#mp}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); + histos.add("Result/Data/antik892invmass", "Invariant mass of K(892) K^{#mp}#pi^{#pm}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); + histos.add("Result/Data/k892invmassLSPP", "Invariant mass of K(892) Like Sign Method K^{#plus}#pi^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K+ + Pi + histos.add("Result/Data/k892invmassLSMM", "Invariant mass of K(892) Like Sign Method K^{#minus}#pi^{#minus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K- + anti-Pi + } + // eta phi QA + if (cFilldeltaEtaPhiPlots) { + histos.add("QAbefore/deltaEta", "deltaEta of kaon and Pion candidates", HistType::kTH1F, {axisEtaPhi}); + histos.add("QAbefore/deltaPhi", "deltaPhi of kaon and Pion candidates", HistType::kTH1F, {axisEtaPhi}); + + histos.add("QAafter/deltaEta", "deltaEta of kaon and Pion candidates", HistType::kTH1F, {axisEtaPhi}); + histos.add("QAafter/deltaPhi", "deltaPhi of kaon and Pion candidates", HistType::kTH1F, {axisEtaPhi}); - if (cfgEventCutQA) { - std::shared_ptr hCutFlow = histos.get(HIST("hEvent_Cut")); - std::vector eventCutLabels = { - "All Events", - "sel8", - Form("|Vz| < %.1f", cfgEventVtxCut.value), - "kIsGoodZvtxFT0vsPV", - "kNoSameBunchPileup", - "kNoTimeFrameBorder", - "kNoITSROFrameBorder", - "kNoCollInTimeRangeStandard", - "kIsGoodITSLayersAll", - Form("Occupancy < %.0f", cfgOccupancyMax.value), - "All passed events"}; - - for (size_t i = 0; i < eventCutLabels.size(); ++i) { - hCutFlow->GetXaxis()->SetBinLabel(i + 1, eventCutLabels[i].c_str()); + histos.add("QAafter/PhiPiafter", "Phi of Pion candidates", HistType::kTH1F, {axisPhi}); + histos.add("QAafter/PhiKaafter", "Phi of kaon candidates", HistType::kTH1F, {axisPhi}); } + + // 3d histogram + histos.add("Result/Data/h3k892invmass", "Invariant mass of K(892) K^{#pm}#pi^{#mp}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); + histos.add("Result/Data/h3antik892invmass", "Invariant mass of K(892) K^{#mp}#pi^{#pm}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); + histos.add("Result/Data/h3k892invmassLSPP", "Invariant mass of K(892) Like Sign Method K^{#plus}#pi^{#plus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); // K+ + Pi + histos.add("Result/Data/h3k892invmassLSMM", "Invariant mass of K(892) Like Sign Method K^{#minus}#pi^{#minus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); // K- + anti-Pi } - } // end of init - //====================== - //|| For LF Analysis - //====================== - using EventCandidates = soa::Join; //, aod::CentFT0Ms, aod::CentFT0As - using EventCandidatesTrue = aod::McCollisions; - using TrackCandidates = soa::Join; - using TrackCandidatesMC = soa::Join; - - // For Mixed Event - using BinningType = ColumnBinningPolicy; - - Partition kaon = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfgTrackTPCPIDnSig); - Partition pion = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaPi) <= cfgTrackTPCPIDnSig); - Partition kaonMC = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfgTrackTPCPIDnSig); - Partition pionMC = !cfgTrackTPCPID || (nabs(aod::pidtpc::tpcNSigmaPi) <= cfgTrackTPCPIDnSig); - - double massKa = o2::constants::physics::MassKPlus; - double massPi = o2::constants::physics::MassPiMinus; - - //================================== - //|| - //|| Helper Templates - //|| - //================================== - template - void fillQA(const bool pass, const objType& obj, const int objecttype = 0) - { - if (objecttype == 1) { - if constexpr (requires { obj.posZ(); }) { - if (!pass) { - histos.fill(HIST("hPosZ_BC"), obj.posZ()); - histos.fill(HIST("hcentFT0C_BC"), obj.centFT0C()); - } else { - histos.fill(HIST("hPosZ_AC"), obj.posZ()); - histos.fill(HIST("hcentFT0C_AC"), obj.centFT0C()); - } - } - } // eventSelection histogram - - if constexpr (requires { obj.tpcCrossedRowsOverFindableCls(); }) { - if (objecttype == 3) { - if (!pass) { - histos.fill(HIST("hDCArToPv_BC"), obj.dcaXY()); - histos.fill(HIST("hDCAzToPv_BC"), obj.dcaZ()); - histos.fill(HIST("hIsPrim_BC"), obj.isPrimaryTrack()); - histos.fill(HIST("hIsGood_BC"), obj.isGlobalTrackWoDCA()); - histos.fill(HIST("hIsPrimCont_BC"), obj.isPVContributor()); - histos.fill(HIST("hFindableTPCClusters_BC"), obj.tpcNClsFindable()); - histos.fill(HIST("hFindableTPCRows_BC"), obj.tpcNClsCrossedRows()); - histos.fill(HIST("hClustersVsRows_BC"), obj.tpcCrossedRowsOverFindableCls()); - histos.fill(HIST("hTPCChi2_BC"), obj.tpcChi2NCl()); - histos.fill(HIST("hITSChi2_BC"), obj.itsChi2NCl()); - histos.fill(HIST("QA_track_pT_BC"), obj.pt()); - } else { - histos.fill(HIST("hDCArToPv_AC"), obj.dcaXY()); - histos.fill(HIST("hDCAzToPv_AC"), obj.dcaZ()); - histos.fill(HIST("hIsPrim_AC"), obj.isPrimaryTrack()); - histos.fill(HIST("hIsGood_AC"), obj.isGlobalTrackWoDCA()); - histos.fill(HIST("hIsPrimCont_AC"), obj.isPVContributor()); - histos.fill(HIST("hFindableTPCClusters_AC"), obj.tpcNClsFindable()); - histos.fill(HIST("hFindableTPCRows_AC"), obj.tpcNClsCrossedRows()); - histos.fill(HIST("hClustersVsRows_AC"), obj.tpcCrossedRowsOverFindableCls()); - histos.fill(HIST("hTPCChi2_AC"), obj.tpcChi2NCl()); - histos.fill(HIST("hITSChi2_AC"), obj.itsChi2NCl()); - histos.fill(HIST("QA_track_pT_AC"), obj.pt()); - } + if (doprocessRotational) { + if (cFill1DQAs) { + histos.add("Result/Data/k892InvMassRotation", "Invariant mass of K(892) Like Sign Method K^{#plus}#pi^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K+ + Pi + histos.add("Result/Data/antik892InvMassRotation", "Invariant mass of K(892) Like Sign Method K^{#minus}#pi^{#minus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K- + anti-Pi } - } // trackSelection - if (objecttype == 4) { - if constexpr (requires { obj.pt(); }) { - if (!pass) { - histos.fill(HIST("QA_nSigma_kaon_TPC_BC"), obj.pt(), obj.tpcNSigmaKa()); - histos.fill(HIST("QA_nSigma_kaon_TOF_BC"), obj.pt(), obj.tofNSigmaKa()); - histos.fill(HIST("QA_kaon_TPC_TOF_BC"), obj.tpcNSigmaKa(), obj.tofNSigmaKa()); - } else { - histos.fill(HIST("QA_nSigma_kaon_TPC_AC"), obj.pt(), obj.tpcNSigmaKa()); - histos.fill(HIST("QA_nSigma_kaon_TOF_AC"), obj.pt(), obj.tofNSigmaKa()); - histos.fill(HIST("QA_kaon_TPC_TOF_AC"), obj.tpcNSigmaKa(), obj.tofNSigmaKa()); - } + histos.add("Result/Data/h3k892InvMassRotation", "Invariant mass of K(892) rotation", kTHnSparseF, {axisMult, axisPt, axisMassk892}); + histos.add("Result/Data/h3antik892InvMassRotation", "Invariant mass of K(892) rotation", kTHnSparseF, {axisMult, axisPt, axisMassk892}); + } + // Mixed event histograms + if (doprocessME) { + if (cFill1DQAs) { + histos.add("Result/Data/k892invmassME_UnlikeSign", "Invariant mass of K(892) mixed event K^{#pm}#pi^{#mp}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); + histos.add("Result/Data/antik892invmassME_UnlikeSign", "Invariant mass of K(892) mixed event K^{#pm}#pi^{#mp}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); } - } // kaon pid Selection - if (objecttype == 5) { - if constexpr (requires { obj.pt(); }) { - if (!pass) { - histos.fill(HIST("QA_nSigma_pion_TPC_BC"), obj.pt(), obj.tpcNSigmaPi()); - histos.fill(HIST("QA_nSigma_pion_TOF_BC"), obj.pt(), obj.tofNSigmaPi()); - histos.fill(HIST("QA_pion_TPC_TOF_BC"), obj.tpcNSigmaPi(), obj.tofNSigmaPi()); - } else { - histos.fill(HIST("QA_nSigma_pion_TPC_AC"), obj.pt(), obj.tpcNSigmaPi()); - histos.fill(HIST("QA_nSigma_pion_TOF_AC"), obj.pt(), obj.tofNSigmaPi()); - histos.fill(HIST("QA_pion_TPC_TOF_AC"), obj.tpcNSigmaPi(), obj.tofNSigmaPi()); + histos.add("Result/Data/h3k892invmassME_UnlikeSign", "Invariant mass of K(892) mixed event K^{#pm}#pi^{#mp}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); + histos.add("Result/Data/h3antik892invmassME_UnlikeSign", "Invariant mass of K(892) mixed event K^{#pm}#pi^{#mp}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); + + if (cFilladditionalMEPlots) { + if (cFill1DQAs) { + histos.add("Result/Data/k892invmassME_LSPP", "Invariant mass of K(892) Like Sign Method K^{#plus}#pi^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K+ + Pi + histos.add("Result/Data/k892invmassME_LSMM", "Invariant mass of K(892) Like Sign Method K^{#minus}#pi^{#minus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassk892}}); // K- + anti-Pi } + histos.add("Result/Data/h3k892invmassME_LSPP", "Invariant mass of K(892) mixed event Like Sign Method K^{#plus}#pi^{#plus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); // K+ + Pi + histos.add("Result/Data/h3k892invmassME_LSMM", "Invariant mass of K(892) mixed event Like Sign Method K^{#minus}#pi^{#minus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassk892}); // K- + anti-Pi } - } // pion pid Selection - } // fill QA - - enum class objectType { MB, - MBRecParticle }; + } - template - void fillMinv(objectType type, const TrackType& trk1, const TrackType& trk2, const ROOT::Math::PxPyPzMVector& lReso, double centrality, bool IsMix, bool flip) - { - double conjugate = trk1.sign() * trk2.sign(); - switch (type) { - case objectType::MB: - if (IsMix && cfgMixedHistos) { - if (conjugate < 0) { - if (!flip) - histos.fill(HIST("hUSS_KPi_Mix"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hUSS_PiK_Mix"), centrality, lReso.Pt(), lReso.M()); - } - } else { - if (conjugate < 0) { - if (!flip) - histos.fill(HIST("hUSS_KPi"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hUSS_PiK"), centrality, lReso.Pt(), lReso.M()); - } else if (conjugate > 0) { - if (!flip) - histos.fill(HIST("hLSS_KPi"), centrality, lReso.Pt(), lReso.M()); - else - histos.fill(HIST("hLSS_PiK"), centrality, lReso.Pt(), lReso.M()); - } - } - break; + // MC QA + if (doprocessEventFactor) { + histos.add("Event/MultiplicityGenEv", "hMCEventIndices", kTH1D, {axisMult}); + histos.add("Event/MultiplicityRecoEv", "Multiplicity of Reconstructed Events", kTH1D, {axisMult}); + } - case objectType::MBRecParticle: - if (cfgMCHistos) { - if (conjugate < 0) { - histos.fill(HIST("hUSS_TrueRec"), centrality, lReso.Pt(), lReso.M()); - } - } - break; - } // switch - } // fillMinv - //====================================================================== + if (doprocessMCGen) { + histos.add("QA/MC/h2GenEtaPt_beforeanycut", " #eta-#it{p}_{T} distribution of Generated K(892); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPtQA}); + histos.add("QA/MC/h2GenPhiRapidity_beforeanycut", " #phi-y distribution of Generated K(892); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); + histos.add("QA/MC/h2GenEtaPt_afterRapcut", " #phi-#it{p}_{T} distribution of Generated K(892); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPtQA}); + histos.add("QA/MC/h2GenPhiRapidity_afterRapcut", " #phi-y distribution of Generated K(892); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); - template - std::pair eventSelection(const EventType event, const bool QA) - { - if (cfgEventCutQA && QA) { - fillQA(false, event, 1); - histos.fill(HIST("hEvent_Cut"), 0); + histos.add("Result/MC/Genk892pt", "pT distribution of True MC K(892)0", kTH2F, {axisPt, axisMult}); + histos.add("Result/MC/Genantik892pt", "pT distribution of True MC Anti-K(892)0", kTH2F, {axisPt, axisMult}); } - if (!event.sel8()) - return {false, 1}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 1); - } - if (std::abs(event.posZ()) > cfgEventVtxCut) - return {false, 2}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 2); - } - if (!event.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) - return {false, 3}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 3); - } - if (!event.selection_bit(aod::evsel::kNoSameBunchPileup)) - return {false, 4}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 4); - } - if (!event.selection_bit(aod::evsel::kNoTimeFrameBorder)) - return {false, 5}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 5); - } - if (!event.selection_bit(aod::evsel::kNoITSROFrameBorder)) - return {false, 6}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 6); - } - if (!event.selection_bit(aod::evsel::kNoCollInTimeRangeStandard)) - return {false, 7}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 7); - } - if (!event.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) - return {false, 8}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 8); - } - if (cfgOccupancySel && (event.trackOccupancyInTimeRange() > cfgOccupancyMax || event.trackOccupancyInTimeRange() < cfgOccupancyMin)) - return {false, 9}; - if (cfgEventCutQA && QA) { - histos.fill(HIST("hEvent_Cut"), 9); + if (doprocessMCRec) { + histos.add("QA/MC/h2RecoEtaPt_after", " #eta-#it{p}_{T} distribution of Reconstructed K(892); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPt}); + histos.add("QA/MC/h2RecoPhiRapidity_after", " #phi-y distribution of Reconstructed K(892); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); + + histos.add("QA/MC/trkDCAxy_pi", "DCAxy distribution of Pion track candidates", HistType::kTHnSparseF, {axisPt, axisDCAxy}); + histos.add("QA/MC/trkDCAxy_ka", "DCAxy distribution of kaon track candidates", HistType::kTHnSparseF, {axisPt, axisDCAxy}); + histos.add("QA/MC/trkDCAz_pi", "DCAz distribution of Pion track candidates", HistType::kTHnSparseF, {axisPt, axisDCAz}); + histos.add("QA/MC/trkDCAz_ka", "DCAz distribution of kaon track candidates", HistType::kTHnSparseF, {axisPt, axisDCAz}); + histos.add("QA/MC/TOF_Nsigma_pi_all", "TOF NSigma for Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TPC_Nsigma_pi_all", "TPC NSigma for Pion;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Pion}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TOF_Nsigma_ka_all", "TOF NSigma for Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TPC_Nsigma_ka_all", "TPC NSigma for Kaon;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + + histos.add("Result/MC/h3k892Recoinvmass", "Invariant mass of Reconstructed MC K(892)0", kTHnSparseF, {axisMult, axisPt, axisMassk892}); + histos.add("Result/MC/h3antik892Recoinvmass", "Invariant mass of Reconstructed MC Anti-K(892)0", kTHnSparseF, {axisMult, axisPt, axisMassk892}); } - if (cfgEventCutQA && QA) { - fillQA(true, event, 1); - histos.fill(HIST("hEvent_Cut"), 10); + if (doprocessSignalLoss) { + histos.add("Result/SignalLoss/GenTruek892pt_den", "True k892 (den)", kTH2F, {axisPt, axisMult}); + histos.add("Result/SignalLoss/GenTrueantik892pt_den", "True anti-k892 (den)", kTH2F, {axisPt, axisMult}); + + histos.add("Result/SignalLoss/GenTruek892pt_num", "True k892 (num)", kTH2F, {axisPt, axisMult}); + histos.add("Result/SignalLoss/GenTrueantik892pt_num", "True anti-k892 (num)", kTH2F, {axisPt, axisMult}); } - return {true, 11}; - }; - template - bool trackSelection(const TracksType track, const bool QA) + // Print output histograms statistics + LOG(info) << "Size of the histograms in Kstar0analysis:"; + histos.print(); + } // end of init + + float massKa = MassKaonCharged; + float massPi = MassPionCharged; + + // Centralicity estimator selection + template + float centEst(Coll collisions) { - if (cfgTrackCutQA && QA) { - fillQA(false, track, 3); + float returnValue = -999.0f; + switch (centEstimator) { + case 0: + returnValue = collisions.centFT0M(); + break; + case 1: + returnValue = collisions.centFT0A(); + break; + case 2: + returnValue = collisions.centFT0C(); + break; + default: + returnValue = collisions.centFT0M(); + break; } + return returnValue; + } + + template + bool isSelected(const Coll& collision, bool fillHist = true) + { + auto applyCut = [&](bool enabled, bool condition, int bin) { + if (!enabled) + return true; + if (!condition) + return false; + if (fillHist) + histos.fill(HIST("CollCutCounts"), bin); + return true; + }; + + if (fillHist) + histos.fill(HIST("CollCutCounts"), 0); - if (cfgTrackGlobalSel && !track.isGlobalTrack()) + if (!applyCut(true, std::abs(collision.posZ()) <= configEvents.cfgEvtZvtx, 1)) return false; - if (track.pt() < cfgTrackMinPt) + + if (!applyCut(configEvents.cfgEvtTriggerTVXSel, + collision.selection_bit(aod::evsel::kIsTriggerTVX), 2)) return false; - if (std::abs(track.eta()) > cfgTrackMaxEta) + + if (!applyCut(configEvents.cfgEvtNoTFBorderCut, + collision.selection_bit(aod::evsel::kNoTimeFrameBorder), 3)) return false; - if (std::abs(track.dcaXY()) > cfgTrackMaxDCArToPVcut) + + if (!applyCut(configEvents.cfgEvtNoITSROFrameBorderCut, + collision.selection_bit(aod::evsel::kNoITSROFrameBorder), 4)) return false; - if (std::abs(track.dcaZ()) > cfgTrackMaxDCAzToPVcut) + + if (!applyCut(configEvents.cfgEvtIsRCTFlagpassed, rctChecker(collision), 5)) return false; - if (cfgTrackPrimaryTrack && !track.isPrimaryTrack()) + + if (!applyCut(configEvents.cfgEvtSel8, collision.sel8(), 6)) return false; - if (cfgTrackGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) + + if (!applyCut(configEvents.cfgEvtIsINELgt0, collision.isInelGt0(), 7)) return false; - if (cfgTrackFindableTPCClusters > 0 && track.tpcNClsFindable() < cfgTrackFindableTPCClusters) + + if (fillHist) + histos.fill(HIST("CollCutCounts"), 8); + + return true; + } + + template + bool trackCut(const TrackType track) + { + // basic track cuts + if (configTracks.cDCAr7SigCut && std::abs(track.dcaXY()) > (0.004f + 0.013f / (track.pt()))) // 7 - Sigma cut + return false; + if (configTracks.cTPCNClsFound && (track.tpcNClsFound() < configTracks.cMinTPCNClsFound)) return false; - if (track.tpcNClsCrossedRows() < cfgTrackTPCCrossedRows) + if (track.tpcNClsCrossedRows() < configTracks.cfgMinCrossedRows) return false; - if (cfgTrackRowsOverFindable > 0 && track.tpcCrossedRowsOverFindableCls() > cfgTrackRowsOverFindable) + if (configTracks.cfgHasTOF && !track.hasTOF()) return false; - if (track.tpcChi2NCl() > cfgTrackTPCChi2) + if (configTracks.cfgPrimaryTrack && !track.isPrimaryTrack()) return false; - if (track.itsChi2NCl() > cfgTrackITSChi2) + if (configTracks.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) return false; - if (cfgTrackConnectedToPV && !track.isPVContributor()) + if (configTracks.cfgPVContributor && !track.isPVContributor()) + return false; + if (configTracks.cfgGlobalTrack && !track.isGlobalTrack()) return false; - if (cfgTrackCutQA && QA) { - fillQA(true, track, 3); - } return true; - }; + } + + // LOGF(info, "AFTER: pt: %f, hasTOF: %d, TPCSigma: %f, TOFSigma: %f, boolTPC: %d, boolTOF: %d, bool: %d", pt, candidate.hasTOF(), + // candidate.tpcNSigmaPr(), candidate.tofNSigmaPr(), tpcPIDPassed, tofPIDPassed, tpcPIDPassed || tofPIDPassed); - template - bool trackPIDKaon(const TrackPID& candidate, const bool QA) + template + bool ptDependentPidPion(const T& candidate) { - double tpcpid = 0; - double tofpid = 0; - bool tpcPIDPassed{false}, tofPIDPassed{false}; + auto vPionTPCPIDpTintv = configPID.pionTPCPIDpTintv.value; + vPionTPCPIDpTintv.insert(vPionTPCPIDpTintv.begin(), configTracks.cMinPtcut); + auto vPionTPCPIDcuts = configPID.pionTPCPIDcuts.value; + auto vPionTOFPIDpTintv = configPID.pionTOFPIDpTintv.value; + auto vPionTPCTOFCombinedpTintv = configPID.pionTPCTOFCombinedpTintv.value; + auto vPionTPCTOFCombinedPIDcuts = configPID.pionTPCTOFCombinedPIDcuts.value; + auto vPionTOFPIDcuts = configPID.pionTOFPIDcuts.value; + + float pt = candidate.pt(); + float ptSwitchToTOF = vPionTPCPIDpTintv.back(); + float tpcNsigmaPi = candidate.tpcNSigmaPi(); + float tofNsigmaPi = candidate.tofNSigmaPi(); + + bool tpcPIDPassed = false; + + // TPC PID (interval check) + for (size_t i = 0; i < vPionTPCPIDpTintv.size() - 1; ++i) { + if (pt > vPionTPCPIDpTintv[i] && pt < vPionTPCPIDpTintv[i + 1]) { + if (std::abs(tpcNsigmaPi) < vPionTPCPIDcuts[i]) + tpcPIDPassed = true; + } + } - if (cfgTrackCutQA && QA) { - fillQA(false, candidate, 4); + // TOF bypass option (for QA or MC) + if (configPID.cByPassTOF) { + return std::abs(tpcNsigmaPi) < vPionTPCPIDcuts.back(); } - // TPC - if (cfgTrackSquarePIDCut) { - if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) - tpcPIDPassed = true; - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaKa()) < cfgTrackTOFPIDnSig) { - tofPIDPassed = true; - } - } else { - if (!cfgTrackTOFHard) { - tofPIDPassed = true; - } else { - tofPIDPassed = false; + // Case 1: No TOF and pt ≤ threshold → accept only via TPC PID + if (!candidate.hasTOF() && pt <= ptSwitchToTOF) { + return tpcPIDPassed; + } + + // Case 2: No TOF but pt > threshold → reject + if (!candidate.hasTOF() && pt > ptSwitchToTOF) { + return false; + } + + // Case 3: Has TOF → use TPC + TOF (square or circular) + if (candidate.hasTOF()) { + if (configPID.cPIDcutType == SquareType) { + // Rectangular cut + for (size_t i = 0; i < vPionTOFPIDpTintv.size(); ++i) { + if (pt < vPionTOFPIDpTintv[i]) { + if (std::abs(tofNsigmaPi) < vPionTOFPIDcuts[i] && + std::abs(tpcNsigmaPi) < vPionTPCPIDcuts.back()) + return true; + } } - } - } // end of square cut - if (cfgTrackCirclePIDCut) { - if (std::abs(candidate.tpcNSigmaKa()) < cfgTrackTPCPIDnSig) - tpcpid = std::abs(candidate.tpcNSigmaKa()); - tofpid = 0; - - if (candidate.hasTOF()) { - tofpid = std::abs(candidate.tofNSigmaKa()); - } else { - if (cfgTrackTOFHard) { - tofpid = 999; + } else if (configPID.cPIDcutType == CircularType) { + // Circular cut + for (size_t i = 0; i < vPionTPCTOFCombinedpTintv.size(); ++i) { + if (pt < vPionTPCTOFCombinedpTintv[i]) { + float combinedSigma2 = + tpcNsigmaPi * tpcNsigmaPi + + tofNsigmaPi * tofNsigmaPi; + if (combinedSigma2 < vPionTPCTOFCombinedPIDcuts[i] * vPionTPCTOFCombinedPIDcuts[i]) + return true; + } } } - if (tpcpid * tpcpid + tofpid * tofpid < cfgTrackCircleValue) { - tpcPIDPassed = true; - tofPIDPassed = true; - } - } // circular cut - - // TPC & TOF - if (tpcPIDPassed && tofPIDPassed) { - if (cfgTrackCutQA && QA) { - fillQA(true, candidate, 4); - } - return true; } + return false; } - template - bool trackPIDPion(const TrackPID& candidate, const bool QA) + template + bool ptDependentPidKaon(const T& candidate) { - double tpcpid = 0; - double tofpid = 0; - bool tpcPIDPassed{false}, tofPIDPassed{false}; + auto vKaonTPCPIDpTintv = configPID.kaonTPCPIDpTintv.value; + vKaonTPCPIDpTintv.insert(vKaonTPCPIDpTintv.begin(), configTracks.cMinPtcut); + auto vKaonTPCPIDcuts = configPID.kaonTPCPIDcuts.value; + auto vKaonTOFPIDpTintv = configPID.kaonTOFPIDpTintv.value; + auto vKaonTPCTOFCombinedpTintv = configPID.kaonTPCTOFCombinedpTintv.value; + auto vKaonTPCTOFCombinedPIDcuts = configPID.kaonTPCTOFCombinedPIDcuts.value; + auto vKaonTOFPIDcuts = configPID.kaonTOFPIDcuts.value; + + float pt = candidate.pt(); + float ptSwitchToTOF = vKaonTPCPIDpTintv.back(); + float tpcNsigmaKa = candidate.tpcNSigmaKa(); + float tofNsigmaKa = candidate.tofNSigmaKa(); + + bool tpcPIDPassed = false; + + // TPC PID interval-based check + for (size_t i = 0; i < vKaonTPCPIDpTintv.size() - 1; ++i) { + if (pt > vKaonTPCPIDpTintv[i] && pt < vKaonTPCPIDpTintv[i + 1]) { + if (std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts[i]) { + tpcPIDPassed = true; + break; + } + } + } - if (cfgTrackCutQA && QA) { - fillQA(false, candidate, 5); + // TOF bypass option + if (configPID.cByPassTOF) { + return std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts.back(); } - // TPC - if (cfgTrackSquarePIDCut) { - if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) - tpcPIDPassed = true; - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaPi()) < cfgTrackTOFPIDnSig) { - tofPIDPassed = true; - } - } else { - if (!cfgTrackTOFHard) { - tofPIDPassed = true; - } else { - tofPIDPassed = false; + // Case 1: No TOF and pt ≤ ptSwitch → use TPC-only + if (!candidate.hasTOF() && pt <= ptSwitchToTOF) { + return tpcPIDPassed; + } + + // Case 2: No TOF but pt > ptSwitch → reject + if (!candidate.hasTOF() && pt > ptSwitchToTOF) { + return false; + } + + // Case 3: TOF is available → apply TPC+TOF PID logic + if (candidate.hasTOF()) { + if (configPID.cPIDcutType == SquareType) { + // Rectangular cut + for (size_t i = 0; i < vKaonTOFPIDpTintv.size(); ++i) { + if (pt < vKaonTOFPIDpTintv[i]) { + if (std::abs(tofNsigmaKa) < vKaonTOFPIDcuts[i] && + std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts.back()) { + return true; + } + } } - } - } // end of square cut - if (cfgTrackCirclePIDCut) { - if (std::abs(candidate.tpcNSigmaPi()) < cfgTrackTPCPIDnSig) - tpcpid = std::abs(candidate.tpcNSigmaPi()); - tofpid = 0; - - if (candidate.hasTOF()) { - tofpid = std::abs(candidate.tofNSigmaPi()); - } else { - if (cfgTrackTOFHard) { - tofpid = 999; + } else if (configPID.cPIDcutType == CircularType) { + // Circular cut + for (size_t i = 0; i < vKaonTPCTOFCombinedpTintv.size(); ++i) { + if (pt < vKaonTPCTOFCombinedpTintv[i]) { + float combinedSigma2 = tpcNsigmaKa * tpcNsigmaKa + + tofNsigmaKa * tofNsigmaKa; + if (combinedSigma2 < vKaonTPCTOFCombinedPIDcuts[i] * vKaonTPCTOFCombinedPIDcuts[i]) { + return true; + } + } } } - if (tpcpid * tpcpid + tofpid * tofpid < cfgTrackCircleValue) { - tpcPIDPassed = true; - tofPIDPassed = true; - } - } // circular cut - - // TPC & TOF - if (tpcPIDPassed && tofPIDPassed) { - if (cfgTrackCutQA && QA) { - fillQA(true, candidate, 5); - } - return true; } + return false; } - template - ROOT::Math::PxPyPzMVector minvReconstruction(const TracksType& trk1, const TracksType& trk2, const bool QA, const bool flip) + auto static constexpr MinPtforProtonRejection = 1.0f; + auto static constexpr MaxPtforProtonRejection = 2.0f; + auto static constexpr MaxnSigmaforProtonRejection = 2.0f; + + template + bool rejectProton(const T& candidate) { - if (!trackSelection(trk1, false) || !trackSelection(trk2, false)) - return {}; - - if (!trackPIDKaon(trk1, QA) || !trackPIDPion(trk2, QA)) - return {}; - - if (trk1.globalIndex() >= trk2.globalIndex()) - return {}; - - ROOT::Math::PxPyPzMVector lDecayDaughter1, lDecayDaughter2, lResonance; - if (!flip) { - lDecayDaughter1 = ROOT::Math::PxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2 = ROOT::Math::PxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massPi); - } else { - lDecayDaughter1 = ROOT::Math::PxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), massPi); - lDecayDaughter2 = ROOT::Math::PxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), massKa); + if (candidate.pt() > MinPtforProtonRejection && candidate.pt() < MaxPtforProtonRejection && !candidate.hasTOF() && candidate.tpcNSigmaPr() < MaxnSigmaforProtonRejection) { + return false; } - lResonance = lDecayDaughter1 + lDecayDaughter2; - - if (std::abs(lResonance.Eta()) > cfgTrackMaxEta) - return {}; - return {lResonance}; + return true; } - template - ROOT::Math::PxPyPzMVector TrueReconstruction(const TracksType& trk1, const TracksType& trk2) - { - double conjugate = trk1.sign() * trk2.sign(); - if (conjugate > 0) - return {}; + auto static constexpr MaxNok892Daughters = 2; - auto particle1 = trk1.mcParticle(); - auto particle2 = trk2.mcParticle(); + template + void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2) + { + auto centrality = centEst(collision); - if (!particle1.has_mothers() || !particle2.has_mothers()) { - return {}; + // Multiplicity correlation calibration plots + if (cFillMultQA) { + if constexpr (IsData) { + histos.fill(HIST("MultCalib/centGloPVpr"), centrality, dTracks1.size(), collision.multNTracksPV()); + histos.fill(HIST("MultCalib/centGloPVka"), centrality, dTracks2.size(), collision.multNTracksPV()); + } } - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& particle1_mom : particle1.template mothers_as()) { - mothers1.push_back(particle1_mom.globalIndex()); - mothers1PDG.push_back(particle1_mom.pdgCode()); + if (cFilladditionalQAeventPlots) { + if constexpr (IsData) { + histos.fill(HIST("QAevent/hnTrksSameE"), dTracks1.size()); + } else if constexpr (IsMix) { + histos.fill(HIST("QAevent/hVertexZMixedE"), collision.posZ()); + histos.fill(HIST("QAevent/hMultiplicityPercentMixedE"), centrality); + histos.fill(HIST("QAevent/hnTrksMixedE"), dTracks1.size()); + } } + // LOG(info) << "After pass, Collision index:" << collision.index() << "multiplicity: " << collision.centFT0M() << endl; - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& particle2_mom : particle2.template mothers_as()) { - mothers2.push_back(particle2_mom.globalIndex()); - mothers2PDG.push_back(particle2_mom.pdgCode()); - } + LorentzVectorPtEtaPhiMass lDecayDaughter1, lDecayDaughter2, lResonance, ldaughterRot, lResonanceRot; - if (mothers1PDG[0] != 313) - return {}; // mother not K*0 - if (mothers2PDG[0] != 313) - return {}; // mothers not K*0 + for (const auto& [trk1, trk2] : combinations(CombinationsFullIndexPolicy(dTracks1, dTracks2))) { + // Full index policy is needed to consider all possible combinations + if (trk1.index() == trk2.index()) + continue; // We need to run (0,1), (1,0) pairs as well. but same id pairs are not needed. - if (mothers1[0] != mothers2[0]) - return {}; // Kaon and pion not from the same K*0 + // apply the track cut + if (!trackCut(trk1) || !trackCut(trk2)) + continue; - if (std::fabs(particle1.pdgCode()) != 211 && std::fabs(particle1.pdgCode()) != 321) - return {}; - if (std::fabs(particle2.pdgCode()) != 211 && std::fabs(particle2.pdgCode()) != 321) - return {}; + //// Initialize variables + // Trk1: Pion + auto isTrk1hasTOF = trk1.hasTOF(); + auto trk1ptPi = trk1.pt(); + auto trk1etaPi = trk1.eta(); + auto trk1phiPi = trk1.phi(); + auto trk1NSigmaPiTPC = trk1.tpcNSigmaPi(); + auto trk1NSigmaPiTOF = (isTrk1hasTOF) ? trk1.tofNSigmaPi() : -999.0f; + + // Trk2: Kaon + auto isTrk2hasTOF = trk2.hasTOF(); + auto trk2ptKa = trk2.pt(); + auto trk2etaKa = trk2.eta(); + auto trk2phiKa = trk2.phi(); + auto trk2NSigmaKaTPC = trk2.tpcNSigmaKa(); + auto trk2NSigmaKaTOF = (isTrk2hasTOF) ? trk2.tofNSigmaKa() : -999.0f; + + auto deltaEta = 0.0f; + auto deltaPhi = 0.0f; + + if (cfgUseDeltaEtaPhiCuts) { + deltaEta = std::abs(trk1etaPi - trk2etaKa); + deltaPhi = std::abs(trk1phiPi - trk2phiKa); + deltaPhi = (deltaPhi > o2::constants::math::PI) ? (o2::constants::math::TwoPI - deltaPhi) : deltaPhi; + if (deltaEta >= cMaxDeltaEtaCut) + continue; + if (deltaPhi >= cMaxDeltaPhiCut) + continue; + } - double track1_mass, track2_mass; - if (std::fabs(particle1.pdgCode()) == 211) { - track1_mass = massPi; - } else { - track1_mass = massKa; - } + //// QA plots before the selection + // --- Track QA all + if constexpr (IsData) { + if (cFillTrackQA) { + histos.fill(HIST("QA/QAbefore/Track/TPC_Nsigma_pi_all"), centrality, trk1ptPi, trk1NSigmaPiTPC); + if (isTrk1hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TOF_Nsigma_pi_all"), centrality, trk1ptPi, trk1NSigmaPiTOF); + histos.fill(HIST("QA/QAbefore/Track/TOF_TPC_Map_pi_all"), trk1NSigmaPiTOF, trk1NSigmaPiTPC); + } + if (!isTrk1hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TPConly_Nsigma_pi_all"), trk1ptPi, trk1NSigmaPiTPC); + } + histos.fill(HIST("QA/QAbefore/Track/TPC_Nsigma_ka_all"), centrality, trk2ptKa, trk2NSigmaKaTPC); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TOF_Nsigma_ka_all"), centrality, trk2ptKa, trk2NSigmaKaTOF); + histos.fill(HIST("QA/QAbefore/Track/TOF_TPC_Map_ka_all"), trk2NSigmaKaTOF, trk2NSigmaKaTPC); + } + if (!isTrk2hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TPConly_Nsigma_ka_all"), trk2ptKa, trk2NSigmaKaTPC); + } - if (std::fabs(particle2.pdgCode()) == 211) { - track2_mass = massPi; - } else { - track2_mass = massKa; - } + histos.fill(HIST("QA/QAbefore/Track/dcaZ_all"), trk1ptPi, trk1.dcaZ()); + histos.fill(HIST("QA/QAbefore/Track/dcaXY_all"), trk1ptPi, trk1.dcaXY()); + histos.fill(HIST("QA/QAbefore/Track/TPC_CR_all"), trk1ptPi, trk1.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAbefore/Track/pT_all"), trk1ptPi); + histos.fill(HIST("QA/QAbefore/Track/eta_all"), trk1etaPi); + } + if (cFilldeltaEtaPhiPlots) { + histos.fill(HIST("QAbefore/deltaEta"), deltaEta); + histos.fill(HIST("QAbefore/deltaPhi"), deltaPhi); + } + } - if (track1_mass == track2_mass) { - return {}; - } + //// Apply the pid selection + if (crejectProton && rejectProton(trk2)) // to remove proton contamination from the kaon track + continue; - ROOT::Math::PxPyPzMVector lTrueDaughter1, lTrueDaughter2, lTrueReso; - lTrueDaughter1 = ROOT::Math::PxPyPzMVector(trk1.px(), trk1.py(), trk1.pz(), track1_mass); - lTrueDaughter2 = ROOT::Math::PxPyPzMVector(trk2.px(), trk2.py(), trk2.pz(), track2_mass); - lTrueReso = lTrueDaughter1 + lTrueDaughter2; + if (!ptDependentPidPion(trk1) || !ptDependentPidKaon(trk2)) + continue; - if (lTrueReso.M() < 0) - return {}; + //// QA plots after the selection + if constexpr (IsData) { // --- PID QA Pion + histos.fill(HIST("QA/QAafter/Pion/TPC_Nsigma_pi_selected"), centrality, trk1ptPi, trk1NSigmaPiTPC); + histos.fill(HIST("QA/QAafter/Pion/TPC_Signal_pi_selected"), trk1.tpcInnerParam(), trk1.tpcSignal()); + if (isTrk1hasTOF) { + histos.fill(HIST("QA/QAafter/Pion/TOF_Nsigma_pi_selected"), centrality, trk1ptPi, trk1NSigmaPiTOF); + histos.fill(HIST("QA/QAafter/Pion/TOF_TPC_Map_pi_selected"), trk1NSigmaPiTOF, trk1NSigmaPiTPC); + } + if (!isTrk1hasTOF) { + histos.fill(HIST("QA/QAafter/Pion/TPC_Nsigma_pi_TPConly_selected"), trk1ptPi, trk1NSigmaPiTPC); + } + histos.fill(HIST("QA/QAafter/Pion/dcaZ_selected"), trk1ptPi, trk1.dcaZ()); + histos.fill(HIST("QA/QAafter/Pion/dcaXY_selected"), trk1ptPi, trk1.dcaXY()); + histos.fill(HIST("QA/QAafter/Pion/TPC_CR_selected"), trk1ptPi, trk1.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAafter/Pion/pT_selected"), trk1ptPi); + histos.fill(HIST("QA/QAafter/Pion/eta_selected"), trk1etaPi); + histos.fill(HIST("QA/QAafter/Pion/TPCnclusterPhipi_selected"), trk1.tpcNClsFound(), trk1phiPi); + + // --- PID QA Kaon + histos.fill(HIST("QA/QAafter/Kaon/TPC_Nsigma_ka_selected"), centrality, trk2ptKa, trk2NSigmaKaTPC); + histos.fill(HIST("QA/QAafter/Kaon/TPC_Signal_ka_selected"), trk2.tpcInnerParam(), trk2.tpcSignal()); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/QAafter/Kaon/TOF_Nsigma_ka_selected"), centrality, trk2ptKa, trk2NSigmaKaTOF); + histos.fill(HIST("QA/QAafter/Kaon/TOF_TPC_Map_ka_selected"), trk2NSigmaKaTOF, trk2NSigmaKaTPC); + } + if (!isTrk2hasTOF) { + histos.fill(HIST("QA/QAafter/Kaon/TPC_Nsigma_ka_TPConly_selected"), trk2ptKa, trk2NSigmaKaTPC); + } + histos.fill(HIST("QA/QAafter/Kaon/dcaZ_selected"), trk2ptKa, trk2.dcaZ()); + histos.fill(HIST("QA/QAafter/Kaon/dcaXY_selected"), trk2ptKa, trk2.dcaXY()); + histos.fill(HIST("QA/QAafter/Kaon/TPC_CR_selected"), trk2ptKa, trk2.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAafter/Kaon/pT_selected"), trk2ptKa); + histos.fill(HIST("QA/QAafter/Kaon/eta_selected"), trk2etaKa); + histos.fill(HIST("QA/QAafter/Kaon/TPCnclusterPhika_selected"), trk2.tpcNClsFound(), trk2phiKa); + + if (cFilldeltaEtaPhiPlots) { + histos.fill(HIST("QAafter/PhiPrafter"), trk1phiPi); + histos.fill(HIST("QAafter/PhiKaafter"), trk2phiKa); + histos.fill(HIST("QAafter/deltaEta"), deltaEta); + histos.fill(HIST("QAafter/deltaPhi"), deltaPhi); + } + } - return {lTrueReso}; - } + //// Resonance reconstruction + lDecayDaughter1 = LorentzVectorPtEtaPhiMass(trk1ptPi, trk1etaPi, trk1phiPi, massPi); + lDecayDaughter2 = LorentzVectorPtEtaPhiMass(trk2ptKa, trk2etaKa, trk2phiKa, massKa); - template - void TrackSlicing(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix, const bool QA) - { - auto tracks1 = kaon->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); - auto tracks2 = pion->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); - auto centrality = collision1.centFT0C(); - - for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - for (bool flip : {false, true}) { - auto lReso = minvReconstruction(trk1, trk2, QA, flip); - if (lReso.M() < 0) + // Apply kinematic opening angle cut + if (cUseOpeningAngleCut) { + auto v1 = lDecayDaughter1.Vect(); + auto v2 = lDecayDaughter2.Vect(); + float alpha = std::acos(v1.Dot(v2) / (v1.R() * v2.R())); + if (alpha > cMinOpeningAngle && alpha < cMaxOpeningAngle) continue; + } - fillMinv(objectType::MB, trk1, trk2, lReso, centrality, IsMix, flip); - } // flip - } // for - } // TrackSlicing + lResonance = lDecayDaughter1 + lDecayDaughter2; - template - void TrackSlicingMC(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool IsMix, const bool QA) - { - auto tracks1 = kaonMC->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); - auto tracks2 = pionMC->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); - auto centrality = collision1.centFT0C(); + auto resonanceMass = lResonance.M(); + auto resonancePt = lResonance.Pt(); + auto resonanceRapidity = lResonance.Rapidity(); - for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) - continue; + if constexpr (IsData || IsMix) { + // Rapidity cut + if (std::abs(resonanceRapidity) > configTracks.cfgCutRapidity) + continue; + } - for (bool flip : {false, true}) { - auto lReso = minvReconstruction(trk1, trk2, QA, flip); - if (lReso.M() < 0) + if (cfgUseCutsOnMother) { + if (resonancePt >= cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (resonanceMass >= cMaxMinvMotherCut) // excluding candidates in overflow continue; + } - fillMinv(objectType::MB, trk1, trk2, lReso, centrality, IsMix, flip); - } // flip - - //============================ - //| True Reconstruction - //============================ - auto lTrueReso = TrueReconstruction(trk1, trk2); - fillMinv(objectType::MBRecParticle, trk1, trk2, lTrueReso, centrality, IsMix, false); - } // tracks - } // TrackSlicingMC - - //======================================================= - //| - //| DATA STUFF (SE) - //| - //======================================================= - int nEvents = 0; - void processDataSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& tracks) - { - if (cDebugLevel > 0) { - nEvents++; - if ((nEvents + 1) % 10000 == 0) { - std::cout << "Processed Data Events: " << nEvents << std::endl; + if (cFilladditionalQAeventPlots) { + if constexpr (IsData) { + histos.fill(HIST("QAevent/hPairsCounterSameE"), 1.0f); + } else if (IsMix) { + histos.fill(HIST("QAevent/hPairsCounterMixedE"), 1.0f); + } } - } - auto [goodEv, code] = eventSelection(collision, true); - histos.fill(HIST("nEvents"), 0.5); + //// Un-like sign pair only + if (trk1.sign() * trk2.sign() < 0) { + if constexpr (IsRot) { + for (int i = 0; i < configBkg.cNofRotations; i++) { + float theta = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / configBkg.rotationalcut, o2::constants::math::PI + o2::constants::math::PI / configBkg.rotationalcut); + if (configBkg.cfgRotPi) { + ldaughterRot = LorentzVectorPtEtaPhiMass(trk1ptPi, trk1etaPi, trk1phiPi + theta, massPi); + lResonanceRot = ldaughterRot + lDecayDaughter2; + } else { + ldaughterRot = LorentzVectorPtEtaPhiMass(trk2ptKa, trk2etaKa, trk2phiKa + theta, massKa); + lResonanceRot = lDecayDaughter1 + ldaughterRot; + } + auto resonanceRotMass = lResonanceRot.M(); + auto resonanceRotPt = lResonanceRot.Pt(); + + // Rapidity cut + if (std::abs(lResonanceRot.Rapidity()) >= configTracks.cfgCutRapidity) + continue; + + if (cfgUseCutsOnMother) { + if (resonanceRotPt >= cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (resonanceRotMass >= cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } + if (trk1.sign() < 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892InvMassRotation"), resonanceRotMass); + } + histos.fill(HIST("Result/Data/h3k892InvMassRotation"), centrality, resonanceRotPt, resonanceRotMass); + } else if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/antik892InvMassRotation"), resonanceRotMass); + } + histos.fill(HIST("Result/Data/h3antik892InvMassRotation"), centrality, resonanceRotPt, resonanceRotMass); + } + } + } + if constexpr (IsData) { + if (trk1.sign() < 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmass"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmass"), centrality, resonancePt, resonanceMass); + } else if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/antik892invmass"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3antik892invmass"), centrality, resonancePt, resonanceMass); + } + } else if (IsMix) { + if (trk1.sign() < 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmassME_UnlikeSign"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmassME_UnlikeSign"), centrality, resonancePt, resonanceMass); + } else if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/antik892invmassME_UnlikeSign"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3antik892invmassME_UnlikeSign"), centrality, resonancePt, resonanceMass); + } + } - if (!goodEv) - return; + // MC + if constexpr (IsMC) { + // now we do mc true + // ------ Temporal lambda function to prevent error in build + auto getMothersIndeces = [&](auto const& theMcParticle) { + std::vector lMothersIndeces{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother index lMother: %d", lMother.globalIndex()); + lMothersIndeces.push_back(lMother.globalIndex()); + } + return lMothersIndeces; + }; + auto getMothersPDGCodes = [&](auto const& theMcParticle) { + std::vector lMothersPDGs{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %d", lMother.pdgCode()); + lMothersPDGs.push_back(lMother.pdgCode()); + } + return lMothersPDGs; + }; + // ------ + std::vector motherstrk1 = {-1, -1}; + std::vector mothersPDGtrk1 = {-1, -1}; + + std::vector motherstrk2 = {-1, -1}; + std::vector mothersPDGtrk2 = {-1, -1}; + + // + // Get the MC particle + const auto& mctrk1 = trk1.mcParticle(); + if (mctrk1.has_mothers()) { + motherstrk1 = getMothersIndeces(mctrk1); + mothersPDGtrk1 = getMothersPDGCodes(mctrk1); + } + while (motherstrk1.size() > MaxNok892Daughters) { + motherstrk1.pop_back(); + mothersPDGtrk1.pop_back(); + } - bool INELgt0 = false; - for (const auto& track : tracks) { - if (!trackSelection(track, true)) - continue; - if (std::fabs(track.eta()) < cfgTrackMaxEta) { - INELgt0 = true; + const auto& mctrk2 = trk2.mcParticle(); + if (mctrk2.has_mothers()) { + motherstrk2 = getMothersIndeces(mctrk2); + mothersPDGtrk2 = getMothersPDGCodes(mctrk2); + } + while (motherstrk2.size() > MaxNok892Daughters) { + motherstrk2.pop_back(); + mothersPDGtrk2.pop_back(); + } + + if (std::abs(mctrk1.pdgCode()) != kPiPlus || std::abs(mctrk2.pdgCode()) != kKPlus) + continue; + + if (motherstrk1[0] != motherstrk2[0]) // Same mother + continue; + + if (std::abs(mothersPDGtrk1[0]) != Pdg::kK0Star892) + continue; + + // LOGF(info, "mother trk1 id: %d, mother trk1: %d, trk1 id: %d, trk1 pdgcode: %d, mother trk2 id: %d, mother trk2: %d, trk2 id: %d, trk2 pdgcode: %d", motherstrk1[0], mothersPDGtrk1[0], trk1.globalIndex(), mctrk1.pdgCode(), motherstrk2[0], mothersPDGtrk2[0], trk2.globalIndex(), mctrk2.pdgCode()); + + if (cUseRapcutMC && std::abs(resonanceRapidity) > configTracks.cfgCutRapidity) // rapidity cut + continue; + + histos.fill(HIST("QA/MC/h2RecoEtaPt_after"), lResonance.Eta(), resonancePt); + histos.fill(HIST("QA/MC/h2RecoPhiRapidity_after"), lResonance.Phi(), resonanceRapidity); + + // Track selection check. + histos.fill(HIST("QA/MC/trkDCAxy_pi"), trk1ptPi, trk1.dcaXY()); + histos.fill(HIST("QA/MC/trkDCAxy_ka"), trk2ptKa, trk2.dcaXY()); + histos.fill(HIST("QA/MC/trkDCAz_pi"), trk1ptPi, trk1.dcaZ()); + histos.fill(HIST("QA/MC/trkDCAz_ka"), trk2ptKa, trk2.dcaZ()); + + histos.fill(HIST("QA/MC/TPC_Nsigma_pi_all"), centrality, trk1ptPi, trk1NSigmaPiTPC); + if (isTrk1hasTOF) { + histos.fill(HIST("QA/MC/TOF_Nsigma_pi_all"), centrality, trk1ptPi, trk1NSigmaPiTOF); + } + histos.fill(HIST("QA/MC/TPC_Nsigma_ka_all"), centrality, trk2ptKa, trk2NSigmaKaTPC); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/MC/TOF_Nsigma_ka_all"), centrality, trk2ptKa, trk2NSigmaKaTOF); + } + + // MC histograms + if (mothersPDGtrk1[0] > 0) { + histos.fill(HIST("Result/MC/h3k892Recoinvmass"), centrality, resonancePt, resonanceMass); + } else { + histos.fill(HIST("Result/MC/h3antik892Recoinvmass"), centrality, resonancePt, resonanceMass); + } + } + } else { + if constexpr (IsData) { + // Like sign pair ++ + if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmassLSPP"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmassLSPP"), centrality, resonancePt, resonanceMass); + } else { // Like sign pair -- + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmassLSMM"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmassLSMM"), centrality, resonancePt, resonanceMass); + } + } else if (IsMix) { + if (cFilladditionalMEPlots) { + // Like sign pair ++ + if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmassME_LSPP"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmassME_LSPP"), centrality, resonancePt, resonanceMass); + } else { // Like sign pair -- + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/k892invmassME_LSMM"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3k892invmassME_LSMM"), centrality, resonancePt, resonanceMass); + } + } + } } } - if (!INELgt0) + } + + void processData(EventCandidates::iterator const& collision, + TrackCandidates const& tracks) + { + if (!isSelected(collision)) // Default event selection return; - histos.fill(HIST("nEvents"), 1.5); - TrackSlicing(collision, tracks, collision, tracks, false, true); + auto centrality = centEst(collision); - } // processSameEvents - PROCESS_SWITCH(kstar0analysis, processDataSameEvent, "process Data Same Event", false); + histos.fill(HIST("Event/posZ"), collision.posZ()); + histos.fill(HIST("Event/centFT0M"), centrality); - //======================================================= - //| - //| DATA STUFF (ME) - //| - //======================================================= - int nEventsMix = 0; - void processDataMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Kstar0analysis, processData, "Process Event for data without partition", false); + + void processRotational(EventCandidates::iterator const& collision, TrackCandidates const& tracks) + { + if (!isSelected(collision, false)) // Default event selection + return; + + if (!collision.isInelGt0()) // <-- + return; + + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Kstar0analysis, processRotational, "Process Rotational Background", false); + + // Processing Event Mixing + using BinningTypeVtxZT0M = ColumnBinningPolicy; + + void processME(EventCandidates const& collision, + TrackCandidates const& tracks) { auto tracksTuple = std::make_tuple(tracks); - BinningType colBinning{{cfgBinsMixVtx, cfgBinsMixMult}, true}; // true is for 'ignore overflows' (true by default) - SameKindPair pairs{colBinning, cfgMixNMixedEvents, -1, collisions, tracksTuple, &cache}; + + BinningTypeVtxZT0M colBinning{{configBkg.cfgVtxBins, configBkg.cfgMultPercentileBins}, true}; + SameKindPair pairs{colBinning, configBkg.nEvtMixing, -1, collision, tracksTuple, &cache}; // -1 is the number of the bin to skip + for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - if (cDebugLevel > 0) { - nEventsMix++; - if ((nEventsMix + 1) % 10000 == 0) { - std::cout << "Processed DATA Mixed Events : " << nEventsMix << std::endl; - } - } - auto [goodEv1, code1] = eventSelection(collision1, false); - auto [goodEv2, code2] = eventSelection(collision2, false); - bool VtxMixFlag = false; - bool CentMixFlag = false; - // bool OccupanacyMixFlag = false; - if (std::fabs(collision1.posZ() - collision2.posZ()) <= cfgVtxMixCut) // set default to maybe 10 - VtxMixFlag = true; - if (std::fabs(collision1.centFT0C() - collision2.centFT0C()) <= cfgVtxMixCut) // set default to maybe 10 - CentMixFlag = true; - - if (!goodEv1 || !goodEv2) + // LOGF(info, "Mixed event collisions: (%d, %d)", collision1.globalIndex(), collision2.globalIndex()); + + // for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) { + // LOGF(info, "Mixed event tracks pair: (%d, %d) from events (%d, %d)", t1.index(), t2.index(), collision1.index(), collision2.index()); + // } + + if (!isSelected(collision1, false)) // Default event selection continue; - if (!CentMixFlag) + + if (!isSelected(collision2, false)) // Default event selection continue; - if (!VtxMixFlag) + + if (!collision1.isInelGt0()) // <-- continue; - TrackSlicing(collision1, tracks1, collision2, tracks2, true, false); - } - } - PROCESS_SWITCH(kstar0analysis, processDataMixedEvent, "process DATA Mixed Event", false); - - //======================================================= - //| - //| MC STUFF (SE) - //| - //========================================================= - int nEventsMC = 0; - void processSameEventMC(EventCandidates::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&) - { - if (cDebugLevel > 0) { - nEventsMC++; - if ((nEventsMC + 1) % 10000 == 0) { - double histmem = histos.getSize(); - std::cout << histmem << std::endl; - std::cout << "process_SameEvent_MC: " << nEventsMC << std::endl; + if (!collision2.isInelGt0()) // <-- + continue; + + if (cFilladditionalQAeventPlots) { + // Fill histograms for the characteristics of the *mixed* events (collision1 and collision2) + // This will show the distribution of events that are actually being mixed. + if (cFill1DQAs) { + histos.fill(HIST("QAevent/hMixPool_VtxZ"), collision1.posZ()); + histos.fill(HIST("QAevent/hMixPool_Multiplicity"), collision1.centFT0M()); + } + histos.fill(HIST("QAevent/hMixPool_VtxZ_vs_Multiplicity"), collision1.posZ(), collision1.centFT0M()); + + // You might also want to fill for collision2 if you want to see both partners' distributions + // histos.fill(HIST("QAevent/hMixPool_VtxZ"), collision2.posZ()); + // histos.fill(HIST("QAevent/hMixPool_Multiplicity"), collision2.centFT0M()); + // histos.fill(HIST("QAevent/hMixPool_VtxZ_vs_Multiplicity"), collision2.posZ(), collision2.centFT0M()); } + fillHistograms(collision1, tracks1, tracks2); } - auto [goodEv, code] = eventSelection(collision, true); - - histos.fill(HIST("nEvents"), 0.5); + } + PROCESS_SWITCH(Kstar0analysis, processME, "Process EventMixing light without partition", false); - if (!goodEv) + void processMCRec(MCEventCandidates::iterator const& collision, + aod::McCollisions const&, + MCTrackCandidates const& tracks, aod::McParticles const&) + { + if (!collision.has_mcCollision()) return; - bool INELgt0 = false; - for (const auto& track : tracks) { - if (!trackSelection(track, true)) - continue; - if (std::fabs(track.eta()) < cfgTrackMaxEta) { - INELgt0 = true; - } - } - if (!INELgt0) + if (!isSelected(collision)) return; - histos.fill(HIST("nEvents"), 1.5); + auto centrality = centEst(collision); - TrackSlicingMC(collision, tracks, collision, tracks, false, true); - } // processSameEvents_MC - PROCESS_SWITCH(kstar0analysis, processSameEventMC, "process Same Event MC", false); + histos.fill(HIST("Event/posZ"), collision.posZ()); + histos.fill(HIST("Event/centFT0M"), centrality); - //======================================================= - //| - //| MC STUFF (ME) - //| - //======================================================= - int nEventsMCMix = 0; - void processMixedEventMC(EventCandidates const& collisions, TrackCandidatesMC const& tracks, aod::McParticles const&) + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Kstar0analysis, processMCRec, "Process Event for MC Rec without partition", false); + + Partition selectedMCParticles = (nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kK0Star892)); // K(892) + + void processMCGen(MCEventCandidates::iterator const& collision, aod::McCollisions const&, aod::McParticles const& mcParticles) { - auto tracksTuple = std::make_tuple(tracks); - BinningType colBinning{{cfgBinsMixVtx, cfgBinsMixMult}, true}; // true is for 'ignore overflows' (true by default) - SameKindPair pairs{colBinning, cfgMixNMixedEvents, -1, collisions, tracksTuple, &cache}; - for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - if (cDebugLevel > 0) { - nEventsMCMix++; - if ((nEventsMCMix + 1) % 10000 == 0) { - std::cout << "Processed Mixed Events: " << nEventsMCMix << std::endl; - } + if (!collision.has_mcCollision()) + return; + + bool isInAfterAllCuts = isSelected(collision, false); + + auto mcPartsAll = mcParticles.sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + // bool isTrueINELgt0 = pwglf::isINELgt0mc(mcPartsAll, pdg); + bool isTrueINELgt0 = collision.isInelGt0(); // <-- + + auto centrality = centEst(collision); + + auto mcParts = selectedMCParticles->sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + + // Not related to the real collisions + for (const auto& part : mcParts) { // loop over all K(892) particles + + std::vector daughterPDGs; + if (part.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(part.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(part.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; } - auto [goodEv1, code1] = eventSelection(collision1, false); - auto [goodEv2, code2] = eventSelection(collision2, false); - if (!goodEv1 || !goodEv2) { + + bool pass1 = std::abs(daughterPDGs[0]) == kKPlus || std::abs(daughterPDGs[1]) == kKPlus; // At least one decay to Kaon + bool pass2 = std::abs(daughterPDGs[0]) == kPiPlus || std::abs(daughterPDGs[1]) == kPiPlus; // At least one decay to Pion + + // Checking if we have both decay products + if (!pass1 || !pass2) + continue; + + // LOGF(info, "Part PDG: %d", part.pdgCode(), "DAU_ID1: %d", pass1, "DAU_ID2: %d", pass2); + + histos.fill(HIST("QA/MC/h2GenEtaPt_beforeanycut"), part.eta(), part.pt()); + histos.fill(HIST("QA/MC/h2GenPhiRapidity_beforeanycut"), part.phi(), part.y()); + + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) // rapidity cut continue; + + if (cfgUseDaughterEtaCutMC) { + for (auto const& daughters : part.daughters_as()) { + if (std::fabs(daughters.eta()) > configTracks.cfgCutEta) + continue; // eta cut for daughters + } // loop over daughters } - TrackSlicingMC(collision1, tracks1, collision2, tracks2, true, false); - } // mixing - } // processMixedEvent_MC - PROCESS_SWITCH(kstar0analysis, processMixedEventMC, "process Mixed Event MC", false); - - //====================================== - //| - //| GENERATED STUFF - //| - //====================================== - int nEventsGen = 0; - void processGen(EventCandidatesTrue::iterator const& collision, soa::SmallGroups> const& recocolls, aod::McParticles const& mcParticles) + histos.fill(HIST("QA/MC/h2GenEtaPt_afterRapcut"), part.eta(), part.pt()); + histos.fill(HIST("QA/MC/h2GenPhiRapidity_afterRapcut"), part.phi(), part.y()); + + if (isInAfterAllCuts && isTrueINELgt0) // after all event selection && INEL>0 + { + if (part.pdgCode() > 0) + histos.fill(HIST("Result/MC/Genk892pt"), part.pt(), centrality); + else + histos.fill(HIST("Result/MC/Genantik892pt"), part.pt(), centrality); + } + } + } + PROCESS_SWITCH(Kstar0analysis, processMCGen, "Process Event for MC only", false); + + void processEventFactor(MCEventCandidates const& collisions, soa::Join const& mcCollisions, aod::McParticles const& mcParticles) { - if (cDebugLevel > 0) { - ++nEventsGen; - if (nEventsGen % 10000 == 0) { - std::cout << "Processed MC (GEN) Events: " << nEventsGen << std::endl; + // Loop on reconstructed collisions + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; } + const auto& mcCollision = collision.mcCollision_as>(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheMC); + + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + bool isInAfterAllCuts = isSelected(collision, false); + + float centrality = mcCollision.centFT0M(); + + if (isTrueINELgt0 && isInAfterAllCuts) + histos.fill(HIST("Event/MultiplicityRecoEv"), centrality); } - //======================= - //|| Event & Signal loss - //======================= - if (cfgMCHistos) { - histos.fill(HIST("nEvents_Gen"), 0.5); + + // Loop on generated collisions to fill the event factor for the INEL>0 correction + for (const auto& mccolls : mcCollisions) { + float centrality = mccolls.centFT0M(); + bool inVtx10 = std::abs(mccolls.posZ()) <= configEvents.cfgEvtZvtx; + + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mccolls.globalIndex(), cacheMC); + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); // QA for Trigger efficiency + + if (inVtx10 && isTrueINELgt0) + histos.fill(HIST("Event/MultiplicityGenEv"), centrality); } + } + PROCESS_SWITCH(Kstar0analysis, processEventFactor, "Process Event factor", false); + + void processSignalLoss(MCEventCandidates const& collisions, soa::Join const& mcCollisions, aod::McParticles const& mcParticles) + { + // Loop on reconstructed collisions + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; + } + const auto& mcCollision = collision.mcCollision_as>(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheMC); + + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + bool isInAfterAllCuts = isSelected(collision, false); + bool inVtx10 = std::abs(mcCollision.posZ()) <= configEvents.cfgEvtZvtx; - for (auto& particle : mcParticles) { - if (particle.pdgCode() != 313) + float centrality = mcCollision.centFT0M(); + + // ===== NUM ===== + if (!(inVtx10 && isTrueINELgt0)) continue; - if (std::fabs(particle.eta()) > cfgTrackMaxEta) + + if (!isInAfterAllCuts) continue; - if (fabs(collision.posZ()) > cfgEventVtxCut) - break; - if (cfgMCHistos) { - histos.fill(HIST("hGen_pT_Raw"), particle.pt()); - } - } // Unreconstructed collisions(=Raw coll) for correction + for (const auto& part : particlesInCollision) { + + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) + continue; - if (recocolls.size() <= 0) { // not reconstructed - if (cfgForceGenReco) { - return; + float pt = part.pt(); + + // kaon + if (std::abs(part.pdgCode()) == Pdg::kK0Star892) { + if (part.pdgCode() > 0) + histos.fill(HIST("Result/SignalLoss/GenTruek892pt_num"), pt, centrality); + else + histos.fill(HIST("Result/SignalLoss/GenTrueantik892pt_num"), pt, centrality); + } } } - double centrality = -1; - for (auto& recocoll : recocolls) { - centrality = recocoll.centFT0C(); - auto [goodEv, code] = eventSelection(recocoll, true); - if (cfgMCHistos) { - histos.fill(HIST("nEvents_Gen"), 1.5); - } - if (!goodEv) - continue; - } // recocolls (=reconstructed collisions) - - //================= - //|| Efficiency - //================= - for (auto& particle : mcParticles) { - if (particle.pdgCode() != 313) - continue; // Not K*0 - if (std::fabs(particle.eta()) > cfgTrackMaxEta) + // Loop on generated collisions to fill the event factor for the INEL>0 correction + for (const auto& mccolls : mcCollisions) { + float centrality = mccolls.centFT0M(); + + bool inVtx10 = std::abs(mccolls.posZ()) <= configEvents.cfgEvtZvtx; + + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mccolls.globalIndex(), cacheMC); + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + + if (!(inVtx10 && isTrueINELgt0)) continue; - if (cfgMCHistos) { - histos.fill(HIST("nEvents_Gen"), 2.5); - histos.fill(HIST("hGen_pT_GoodEv"), centrality, particle.pt()); - } // cfgMCHistos - } // loop over particles - } // processMCTrue - PROCESS_SWITCH(kstar0analysis, processGen, "process Generated Particles", false); + for (const auto& part : particlesInCollision) { - void processEventsDummy(EventCandidates::iterator const&, TrackCandidates const&) - { - return; + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) + continue; + + // ========================= + // ===== K(892) ====== + // ========================= + if (std::abs(part.pdgCode()) == Pdg::kK0Star892) { + + std::vector daughterPDGs; + if (part.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(part.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(part.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } + + bool pass1 = std::abs(daughterPDGs[0]) == kKPlus || std::abs(daughterPDGs[1]) == kKPlus; // At least one decay to Kaon + bool pass2 = std::abs(daughterPDGs[0]) == kPiPlus || std::abs(daughterPDGs[1]) == kPiPlus; // At least one decay to Pion + + // Checking if we have both decay products + if (!pass1 || !pass2) + continue; + + if (part.pdgCode() > 0) + histos.fill(HIST("Result/SignalLoss/GenTruek892pt_den"), part.pt(), centrality); + else + histos.fill(HIST("Result/SignalLoss/GenTrueantik892pt_den"), part.pt(), centrality); + } + } + } } - PROCESS_SWITCH(kstar0analysis, processEventsDummy, "dummy", false); -}; // kstar0analysis + PROCESS_SWITCH(Kstar0analysis, processSignalLoss, "Process SignalLoss", false); +}; // Kstar0analysis WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; }; diff --git a/PWGLF/Tasks/Resonances/phi1020analysis.cxx b/PWGLF/Tasks/Resonances/phi1020analysis.cxx index 609a69af24d..5aa6fcd35ce 100644 --- a/PWGLF/Tasks/Resonances/phi1020analysis.cxx +++ b/PWGLF/Tasks/Resonances/phi1020analysis.cxx @@ -10,13 +10,14 @@ // or submit itself to any jurisdiction. /// \file phi1020analysis.cxx -/// \brief Reconstruction of Phi yield through track-track Minv correlations for resonance OO analysis -/// -/// -/// \author Adrian Fereydon Nassirpour #include #include #include @@ -32,192 +34,383 @@ #include #include #include +#include #include #include #include #include +#include #include +#include #include -#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include +#include +#include +#include #include -#include -#include +#include #include -#include #include -#include - using namespace o2; +using namespace o2::soa; +using namespace o2::aod; +using namespace o2::aod::rctsel; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::constants::physics; -struct phi1020analysis { +using LorentzVectorPtEtaPhiMass = ROOT::Math::PtEtaPhiMVector; +enum TrackSelectionType { + AllTracks = 0, + GlobalTracks, + GlobalTracksWoPtEta, + GlobalTracksWoDCA, + QualityTracks, + InAcceptanceTracks, +}; + +enum PIDCutType { + SquareType = 1, + CircularType, +}; + +struct Phi1020analysis { + // Define slice per collision + Preslice perCollision = o2::aod::track::collisionId; SliceCache cache; - Preslice perCollision = aod::track::collisionId; + Preslice perMcCollision = o2::aod::mcparticle::mcCollisionId; + SliceCache cacheMC; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - // Event configurables - Configurable cfg_Event_Sel{"cfg_Event_Sel", "sel8", "choose event selection"}; - Configurable cfg_Event_VtxCut{"cfg_Event_VtxCut", 10.0, "V_z cut selection"}; - Configurable cfg_Event_Timeframe{"cfg_Event_Timeframe", true, "Timeframe border cut"}; - Configurable cfg_Event_Timerange{"cfg_Event_Timerange", true, "Timerange border cut"}; - Configurable cfg_Event_Centrality{"cfg_Event_Centrality", true, "Centrality cut"}; - Configurable cfg_Event_CentralityMax{"cfg_Event_CentralityMax", 100, "CentralityMax cut"}; - Configurable cfg_Event_Pileup{"cfg_Event_Pileup", true, "Pileup border cut"}; - Configurable cfg_Event_OccupancyCut{"cfg_Event_OccupancyCut", true, "Occupancy border cut"}; - Configurable cfg_Event_MaxOccupancy{"cfg_Event_MaxOccupancy", 1, "Max TPC Occupancy"}; + Service pdg; + RCTFlagsChecker rctChecker; + + struct : ConfigurableGroup { + Configurable cfgEvtZvtx{"cfgEvtZvtx", 10.0f, "Evt sel: Max. z-Vertex (cm)"}; + Configurable cfgEvtTriggerTVXSel{"cfgEvtTriggerTVXSel", true, "Evt sel: triggerTVX selection (MB)"}; + Configurable cfgEvtNoTFBorderCut{"cfgEvtNoTFBorderCut", true, "Evt sel: apply TF border cut"}; + Configurable cfgEvtNoITSROFrameBorderCut{"cfgEvtNoITSROFrameBorderCut", false, "Evt sel: apply NoITSRO border cut"}; + Configurable cfgEvtIsRCTFlagpassed{"cfgEvtIsRCTFlagpassed", false, "Evt sel: apply RCT flag selection"}; + Configurable cfgEvtRCTFlagCheckerLabel{"cfgEvtRCTFlagCheckerLabel", "CBT_hadronPID", "Evt sel: RCT flag checker label"}; + Configurable cfgEvtRCTFlagCheckerZDCCheck{"cfgEvtRCTFlagCheckerZDCCheck", false, "Evt sel: RCT flag checker ZDC check"}; + Configurable cfgEvtRCTFlagCheckerLimitAcceptAsBad{"cfgEvtRCTFlagCheckerLimitAcceptAsBad", true, "Evt sel: RCT flag checker treat Limited Acceptance As Bad"}; + Configurable cfgEvtSel8{"cfgEvtSel8", false, "Evt Sel 8 check for offline selection"}; + Configurable cfgEvtIsINELgt0{"cfgEvtIsINELgt0", false, "Evt sel: apply INEL>0 selection"}; + } configEvents; + + struct : ConfigurableGroup { + // Pre-selection Track cuts + Configurable trackSelection{"trackSelection", 0, "Track selection: 0 -> No Cut, 1 -> kGlobalTrack, 2 -> kGlobalTrackWoPtEta, 3 -> kGlobalTrackWoDCA, 4 -> kQualityTracks, 5 -> kInAcceptanceTracks"}; + Configurable cMinPtcut{"cMinPtcut", 0.15f, "Minimal pT for tracks"}; + Configurable cMinTPCNClsFound{"cMinTPCNClsFound", 120, "minimum TPCNClsFound value for good track"}; + Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; + Configurable cfgCutRapidity{"cfgCutRapidity", 0.5f, "rapidity range for particles"}; + Configurable cfgMinCrossedRows{"cfgMinCrossedRows", 70, "min crossed rows for good track"}; + + // DCA Selections + // DCAr to PV + Configurable cMaxDCArToPVcut{"cMaxDCArToPVcut", 0.1f, "Track DCAr cut to PV Maximum"}; + // DCAz to PV + Configurable cMaxDCAzToPVcut{"cMaxDCAzToPVcut", 0.1f, "Track DCAz cut to PV Maximum"}; + + // Track selections + Configurable cfgPrimaryTrack{"cfgPrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgGlobalWoDCATrack{"cfgGlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) + Configurable cfgGlobalTrack{"cfgGlobalTrack", false, "Global track selection"}; // kGoldenChi2 | kDCAxy | kDCAz + Configurable cfgPVContributor{"cfgPVContributor", false, "PV contributor track selection"}; // PV Contriuibutor + Configurable cfgHasTOF{"cfgHasTOF", false, "Require TOF"}; + Configurable cfgUseTPCRefit{"cfgUseTPCRefit", false, "Require TPC Refit"}; + Configurable cfgUseITSRefit{"cfgUseITSRefit", false, "Require ITS Refit"}; + Configurable cTPCNClsFound{"cTPCNClsFound", false, "Switch to turn on/off TPCNClsFound cut"}; + Configurable cDCAr7SigCut{"cDCAr7SigCut", false, "Track DCAr 7 Sigma cut to PV Maximum"}; + } configTracks; + + struct : ConfigurableGroup { + /// PID Selections + Configurable pidnSigmaPreSelectionCut{"pidnSigmaPreSelectionCut", 4.0f, "pidnSigma Cut for pre-selection of tracks"}; + Configurable cByPassTOF{"cByPassTOF", false, "By pass TOF PID selection"}; // By pass TOF PID selection + Configurable cPIDcutType{"cPIDcutType", 2, "cPIDcutType = 1 for square cut, 2 for circular cut"}; // By pass TOF PID selection + + // Kaon + Configurable> kaonTPCPIDpTintv{"kaonTPCPIDpTintv", {0.5f}, "pT intervals for Kaon TPC PID cuts"}; + Configurable> kaonTPCPIDcuts{"kaonTPCPIDcuts", {2}, "nSigma list for Kaon TPC PID cuts"}; + Configurable> kaonTOFPIDpTintv{"kaonTOFPIDpTintv", {999.0f}, "pT intervals for Kaon TOF PID cuts"}; + Configurable> kaonTOFPIDcuts{"kaonTOFPIDcuts", {2}, "nSigma list for Kaon TOF PID cuts"}; + Configurable> kaonTPCTOFCombinedpTintv{"kaonTPCTOFCombinedpTintv", {999.0f}, "pT intervals for Kaon TPC-TOF PID cuts"}; + Configurable> kaonTPCTOFCombinedPIDcuts{"kaonTPCTOFCombinedPIDcuts", {2}, "nSigma list for Kaon TPC-TOF PID cuts"}; + } configPID; + + struct : ConfigurableGroup { + /// Event Mixing + Configurable nEvtMixing{"nEvtMixing", 10, "Number of events to mix"}; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -10.0f, -8.0f, -6.0f, -4.0f, -2.0f, 0.0f, 2.0f, 4.0f, 6.0f, 8.0f, 10.0f}, "Mixing bins - z-vertex"}; + ConfigurableAxis cfgMultPercentileBins{"cfgMultPercentileBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "Mixing bins - multiplicity"}; + + // Rotational background + Configurable rotationalcut{"rotationalcut", 10, "Cut value (Rotation angle pi - pi/cut and pi + pi/cut)"}; + Configurable cNofRotations{"cNofRotations", 3, "Number of random rotations in the rotational background"}; + Configurable cfgRotTrk1{"cfgRotTrk1", true, "rotate Kaon 1"}; + } configBkg; + + // Additional purity check + Configurable crejectPion{"crejectPion", false, "Switch to turn on/off pion contamination"}; + Configurable cUseOpeningAngleCut{"cUseOpeningAngleCut", false, "Kinematic Cuts for p-K pair opening angle"}; + Configurable cMinOpeningAngle{"cMinOpeningAngle", 1.4f, "Minimum opening angle between daughters"}; + Configurable cMaxOpeningAngle{"cMaxOpeningAngle", 2.4f, "Maximum opening angle between daughters"}; + Configurable cfgUseDeltaEtaPhiCuts{"cfgUseDeltaEtaPhiCuts", false, "Switch to turn on/off delta eta and delta phi cuts"}; + Configurable cfgUseDaughterEtaCutMC{"cfgUseDaughterEtaCutMC", false, "Switch to turn on/off eta cuts for daughters in MC"}; + + // MC selection cut + Configurable cUseRapcutMC{"cUseRapcutMC", true, "Use rapidity cut for MC"}; + + // cuts on mother + Configurable cfgUseCutsOnMother{"cfgUseCutsOnMother", false, "Enable additional cuts on mother"}; + Configurable cMaxPtMotherCut{"cMaxPtMotherCut", 15.0f, "Maximum pt of mother cut"}; + Configurable cMaxMinvMotherCut{"cMaxMinvMotherCut", 1.5f, "Maximum Minv of mother cut"}; + Configurable cMaxDeltaEtaCut{"cMaxDeltaEtaCut", 0.7f, "Maximum deltaEta between daughters"}; + Configurable cMaxDeltaPhiCut{"cMaxDeltaPhiCut", 1.5f, "Maximum deltaPhi between daughters"}; + + // switches + Configurable cFillMultQA{"cFillMultQA", false, "Turn on/off additional QA plots"}; + Configurable cFillTrackQA{"cFillTrackQA", false, "Turn on/off additional QA plots"}; + Configurable cFilladditionalQAeventPlots{"cFilladditionalQAeventPlots", false, "Additional QA event plots"}; + Configurable cFilladditionalMEPlots{"cFilladditionalMEPlots", false, "Additional Mixed event plots"}; + Configurable cFilldeltaEtaPhiPlots{"cFilldeltaEtaPhiPlots", false, "Enable additional cuts on daughters"}; + Configurable cFill1DQAs{"cFill1DQAs", false, "Invariant mass 1D"}; Configurable centEstimator{"centEstimator", 0, "Select centrality estimator: 0 - FT0M, 1 - FT0A, 2 - FT0C"}; - // Track configurables - Configurable cfg_Track_Sel{"cfg_Track_Sel", "globalTracks", "set track selections"}; - Configurable cfg_Track_MinPt{"cfg_Track_MinPt", 0.15, "set track min pT"}; - Configurable cfg_Track_MaxEta{"cfg_Track_MaxEta", 0.9, "set track max Eta"}; - Configurable cfg_Track_MaxDCArToPVcut{"cfg_Track_MaxDCArToPVcut", 0.5, "Track DCAr cut to PV Maximum"}; - Configurable cfg_Track_MaxDCAzToPVcut{"cfg_Track_MaxDCAzToPVcut", 2.0, "Track DCAz cut to PV Maximum"}; - Configurable cfg_Track_PrimaryTrack{"cfg_Track_PrimaryTrack", true, "Primary track selection"}; // kGoldenChi2 | kDCAxy | kDCAz - Configurable cfg_Track_ConnectedToPV{"cfg_Track_ConnectedToPV", true, "PV contributor track selection"}; // PV Contributor - Configurable cfg_Track_GlobalWoDCATrack{"cfg_Track_GlobalWoDCATrack", true, "Global track selection without DCA"}; // kQualityTracks (kTrackType | kTPCNCls | kTPCCrossedRows | kTPCCrossedRowsOverNCls | kTPCChi2NDF | kTPCRefit | kITSNCls | kITSChi2NDF | kITSRefit | kITSHits) | kInAcceptanceTracks (kPtRange | kEtaRange) - Configurable cfg_Track_nFindableTPCClusters{"cfg_Track_FindableTPCClusters", 50, "nFindable TPC Clusters"}; - Configurable cfg_Track_nTPCCrossedRows{"cfg_Track_TPCCrossedRows", 70, "nCrossed TPC Rows"}; - Configurable cfg_Track_nRowsOverFindable{"cfg_Track_RowsOverFindable", 1.2, "nRowsOverFindable TPC CLusters"}; - Configurable cfg_Track_nTPCChi2{"cfg_Track_TPCChi2", 4.0, "nTPC Chi2 per Cluster"}; - Configurable cfg_Track_nITSChi2{"cfg_Track_ITSChi2", 36.0, "nITS Chi2 per Cluster"}; - Configurable cfg_Track_TPCPID{"cfg_Track_TPCPID", true, "Enables TPC PID"}; - Configurable cfg_Track_TOFPID{"cfg_Track_TOFPID", true, "Enables TOF PID"}; - Configurable cfg_Track_Hard_TOFPID{"cfg_Track_Hard_TOFPID", true, "Enables STRICT TOF Reqruirement"}; - Configurable cfg_Track_TPCPID_nSig{"cfg_Track_TPCPID_nSig", 4, "nTPC PID sigma"}; - Configurable cfg_Track_TOFPID_nSig{"cfg_Track_TOFPID_nSig", 4, "nTOF PID sigma"}; - Configurable cfg_Track_Explicit_PID{"cfg_Track_Explicit_PID", true, "Enables explicit pid cehck"}; - Configurable cDebugLevel{"cDebugLevel", 0, "Resolution of Debug"}; - - // Pair configurables - Configurable cfg_Pair_MinvBins{"cfg_Pair_MinvBins", 300, "Number of bins for Minv axis"}; - Configurable cfg_Pair_MinvMin{"cfg_Pair_MinvMin", 0.90, "Minimum Minv value"}; - Configurable cfg_Pair_MinvMax{"cfg_Pair_MinvMax", 1.50, "Maximum Minv value"}; - - Configurable cfg_Mix_NMixedEvents{"cfg_Mix_NMixedEvents", 5, "Number of mixed events per event"}; - - // MCGen configurables - Configurable cfg_Force_GenReco{"cfg_Force_GenReco", false, "Only consider events which are reconstructed (neglect event-loss)"}; - Configurable cfg_Force_BR{"cfg_Force_BR", false, "Only consider phi->K+K-"}; - Configurable cfg_Force_Kaon_Acceptence{"cfg_Force_Kaon_Acceptence", false, "Only consider phi's whose daughters decay inside acceptence (no signal loss)"}; - - // Histogram Configurables - Configurable cfg_Event_CutQA{"cfg_Event_CutsQA", true, "Enables Track QA plots"}; - Configurable cfg_Track_CutQA{"cfg_Track_CutsQA", true, "Enables Track QA plots"}; - - // Configurables for axis - ConfigurableAxis binsDCAz{"binsDCAz", {40, -0.2, 0.2}, ""}; - ConfigurableAxis binsDCAxy{"binsDCAxy", {40, -0.2, 0.2}, ""}; - ConfigurableAxis cfg_bins_Cent{"cfg_bins_Cent", {VARIABLE_WIDTH, 0.0, 1.0, 5.0, 10.0, 20.0, 30.0, 40.0, 50.0, 60.0, 70.0, 80.0, 90.0, 100.0, 110.0}, "Binning of the centrality axis"}; - ConfigurableAxis cfg_bins_MixVtx{"cfg_bins_MixVtx", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; - ConfigurableAxis cfg_bins_MixMult{"cfg_bins_MixMult", {VARIABLE_WIDTH, 0.0f, 1.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f}, "Mixing bins - z-vertex"}; + TRandom* rn = new TRandom(); + + // Pre-filters for efficient process + Filter zVtxFilter = (nabs(o2::aod::collision::posZ) <= configEvents.cfgEvtZvtx); + Filter collisionFilterMC = nabs(aod::mccollision::posZ) <= configEvents.cfgEvtZvtx; + + Filter acceptanceFilter = (nabs(aod::track::eta) < configTracks.cfgCutEta && nabs(aod::track::pt) > configTracks.cMinPtcut) && + (nabs(aod::track::dcaXY) < configTracks.cMaxDCArToPVcut) && (nabs(aod::track::dcaZ) < configTracks.cMaxDCAzToPVcut); + + using EventCandidates = soa::Join; + using TrackCandidates = soa::Filtered>; + + // for MC reco + using MCEventCandidates = soa::Join; + using MCTrackCandidates = soa::Filtered>; + + Partition posKaon = aod::track::signed1Pt > static_cast(0); + Partition negKaon = aod::track::signed1Pt < static_cast(0); + + Partition PosKaon_MC = aod::track::signed1Pt > static_cast(0); + Partition NegKaon_MC = aod::track::signed1Pt < static_cast(0); + + /// Figures + ConfigurableAxis binsPt{"binsPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.12f, 0.14f, 0.16f, 0.18f, 0.2f, 0.25f, 0.3f, 0.35f, 0.4f, 0.45f, 0.5f, 0.55f, 0.6f, 0.65f, 0.7f, 0.75f, 0.8f, 0.85f, 0.9f, 0.95f, 1.0f, 1.1f, 1.2f, 1.25f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.75f, 1.8f, 1.9f, 2.0f, 2.1f, 2.2f, 2.3f, 2.4f, 2.5f, 2.6f, 2.7f, 2.8f, 2.9f, 3.0f, 3.1f, 3.2f, 3.3f, 3.4f, 3.6f, 3.7f, 3.8f, 3.9f, 4.0f, 4.1f, 4.2f, 4.5f, 4.6f, 4.8f, 4.9f, 5.0f, 5.5f, 5.6f, 6.0f, 6.4f, 6.5f, 7.0f, 7.2f, 8.0f, 9.0f, 9.5f, 9.6f, 10.0f, 11.0f, 11.5f, 12.0f, 13.0f, 14.0f, 14.4f, 15.0f, 16.0f, 18.0f, 19.2f, 20.0f}, "Binning of the pT axis"}; + ConfigurableAxis binsEta{"binsEta", {150, -1.5f, 1.5f}, ""}; + ConfigurableAxis binsMass{"binsMass", {300, 0.9f, 1.6f}, "Invariant Mass (GeV/#it{c}^2)"}; + ConfigurableAxis binsMult{"binsMult", {105, 0.0f, 105.0f}, "mult_{FT0M}"}; + ConfigurableAxis binsDCAz{"binsDCAz", {40, -0.2f, 0.2f}, ""}; + ConfigurableAxis binsDCAxy{"binsDCAxy", {40, -0.2f, 0.2f}, ""}; + ConfigurableAxis binsTPCXrows{"binsTPCXrows", {100, 60, 160}, ""}; + ConfigurableAxis binsnSigma{"binsnSigma", {130, -6.5f, 6.5f}, ""}; + ConfigurableAxis binsnTPCSignal{"binsnTPCSignal", {1000, 0, 1000}, ""}; + ConfigurableAxis binsEtaPhi{"binsEtaPhi", {350, -3.5f, 3.5f}, ""}; void init(o2::framework::InitContext&) { - const AxisSpec MinvAxis = {cfg_Pair_MinvBins, cfg_Pair_MinvMin, cfg_Pair_MinvMax}; - const AxisSpec PtAxis = {200, 0, 20.0}; - const AxisSpec MultAxis = {100, 0, 100}; - const AxisSpec dRAxis = {100, 0, 100}; - const AxisSpec pidAxis = {100, -5, 5}; - const AxisSpec axisDCAz{binsDCAz, "DCA_{z}"}; - const AxisSpec axisDCAxy{binsDCAxy, "DCA_{XY}"}; - - // Event QA - if (cfg_Event_CutQA) { - histos.add("hPosZ_BC", "PosZ_BC", kTH1F, {{240, -12.0, 12.0}}); - histos.add("hcentFT0_BC", "centFT0_BC", kTH1F, {{110, 0.0, 110.0}}); - histos.add("hOccupancy_BC", "Occupancy_BC", kTH1F, {{100, 0.0, 20000}}); - // - histos.add("hcentFT0_AC", "centFT0_AC", kTH1F, {{110, 0.0, 110.0}}); - histos.add("hPosZ_AC", "PosZ_AC", kTH1F, {{240, -12.0, 12.0}}); - histos.add("hOccupancy_AC", "Occupancy_AC", kTH1F, {{100, 0.0, 20000}}); + rctChecker.init(configEvents.cfgEvtRCTFlagCheckerLabel, configEvents.cfgEvtRCTFlagCheckerZDCCheck, configEvents.cfgEvtRCTFlagCheckerLimitAcceptAsBad); + + // axes + AxisSpec axisPt{binsPt, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec axisEta{binsEta, "#eta"}; + AxisSpec axisRap{binsEta, "#it{y}"}; + AxisSpec axisMassphi1020{binsMass, "Invariant Mass (GeV/#it{c}^2)"}; + AxisSpec axisMult{binsMult, "mult_{FT0M}"}; + AxisSpec axisDCAz{binsDCAz, "DCA_{z}"}; + AxisSpec axisDCAxy{binsDCAxy, "DCA_{XY}"}; + AxisSpec axisTPCXrow{binsTPCXrows, "#Xrows_{TPC}"}; + AxisSpec axisPIDQA{binsnSigma, "#sigma"}; + AxisSpec axisTPCSignal{binsnTPCSignal, ""}; + AxisSpec axisEtaPhi{binsEtaPhi, ""}; + AxisSpec axisPhi{350, 0, 7, "#Phi"}; + AxisSpec axisMultMix{configBkg.cfgMultPercentileBins, "Multiplicity Percentile"}; + AxisSpec axisVtxMix{configBkg.cfgVtxBins, "Vertex Z (cm)"}; + + histos.add("CollCutCounts", "No. of event after cuts", kTH1I, {{10, 0, 10}}); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(1, "All Events"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(2, "|Vz| < cut"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(3, "kIsTriggerTVX"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(4, "kNoTimeFrameBorder"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(5, "kNoITSROFrameBorder"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(6, "rctChecker"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(7, "sel8"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(8, "IsINELgt0"); + histos.get(HIST("CollCutCounts"))->GetXaxis()->SetBinLabel(9, "All Passed Events"); + + histos.add("Event/posZ", "; vtx_{z} (cm); Entries", HistType::kTH1F, {{250, -12.5, 12.5}}); + histos.add("Event/centFT0M", "; FT0M Percentile; Entries", HistType::kTH1F, {{110, 0, 110}}); + + if (cFilladditionalQAeventPlots) { + // event histograms + if (doprocessData) { + histos.add("QAevent/hPairsCounterSameE", "total valid no. of pairs sameE", HistType::kTH1F, {{1, 0.5f, 1.5f}}); + histos.add("QAevent/hnTrksSameE", "n tracks per event SameE", HistType::kTH1F, {{1000, 0.0, 1000.0}}); + } + // Gen on Mixed event + if (doprocessME) { + + // Histograms for Mixed Event Pool characteristics + histos.add("QAevent/hMixPool_VtxZ", "Mixed Event Pool: Vertex Z;Vtx Z (cm);Counts", HistType::kTH1F, {axisVtxMix}); + histos.add("QAevent/hMixPool_Multiplicity", "Mixed Event Pool: Multiplicity;Multiplicity;Counts", HistType::kTH1F, {axisMultMix}); + histos.add("QAevent/hMixPool_VtxZ_vs_Multiplicity", "Mixed Event Pool: Vertex Z vs Multiplicity;Counts", HistType::kTH2F, {axisVtxMix, axisMultMix}); + + histos.add("QAevent/hPairsCounterMixedE", "total valid no. of pairs mixedE", HistType::kTH1F, {{1, 0.5f, 1.5f}}); + histos.add("QAevent/hVertexZMixedE", "Collision Vertex Z position", HistType::kTH1F, {{100, -15.0f, 15.0f}}); + histos.add("QAevent/hMultiplicityPercentMixedE", "Multiplicity percentile of collision", HistType::kTH1F, {{120, 0.0f, 120.0f}}); + histos.add("QAevent/hnTrksMixedE", "n tracks per event MixedE", HistType::kTH1F, {{1000, 0.0f, 1000.0f}}); + } } - // Track QA - if (cfg_Track_CutQA) { - histos.add("hDCArToPv_BC", "DCArToPv_BC", kTH1F, {axisDCAxy}); - histos.add("hDCAzToPv_BC", "DCAzToPv_BC", kTH1F, {axisDCAz}); - histos.add("hIsPrim_BC", "hIsPrim_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsGood_BC", "hIsGood_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsPrimCont_BC", "hIsPrimCont_BC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hFindableTPCClusters_BC", "hFindableTPCClusters_BC", kTH1F, {{200, 0, 200}}); - histos.add("hFindableTPCRows_BC", "hFindableTPCRows_BC", kTH1F, {{200, 0, 200}}); - histos.add("hClustersVsRows_BC", "hClustersVsRows_BC", kTH1F, {{200, 0, 2}}); - histos.add("hTPCChi2_BC", "hTPCChi2_BC", kTH1F, {{200, 0, 100}}); - histos.add("hITSChi2_BC", "hITSChi2_BC", kTH1F, {{200, 0, 100}}); - histos.add("hTPC_nSigma_BC", "hTPC_nSigma_BC", kTH1F, {pidAxis}); - histos.add("hTOF_nSigma_BC", "hTOF_nSigma_BC", kTH1F, {pidAxis}); - histos.add("hTPC_nSigma_v_pt_BC", "hTPC_nSigma_v_pt_BC", HistType::kTHnSparseD, {pidAxis, PtAxis}); - histos.add("hTOF_nSigma_v_pt_BC", "hTOF_nSigma_v_pt_BC", HistType::kTHnSparseD, {pidAxis, PtAxis}); - // - histos.add("hDCArToPv_AC", "DCArToPv_AC", kTH1F, {axisDCAxy}); - histos.add("hDCAzToPv_AC", "DCAzToPv_AC", kTH1F, {axisDCAz}); - histos.add("hIsPrim_AC", "hIsPrim_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsGood_AC", "hIsGood_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hIsPrimCont_AC", "hIsPrimCont_AC", kTH1F, {{2, -0.5, 1.5}}); - histos.add("hFindableTPCClusters_AC", "hFindableTPCClusters_AC", kTH1F, {{200, 0, 200}}); - histos.add("hFindableTPCRows_AC", "hFindableTPCRows_AC", kTH1F, {{200, 0, 200}}); - histos.add("hClustersVsRows_AC", "hClustersVsRows_AC", kTH1F, {{200, 0, 2}}); - histos.add("hTPCChi2_AC", "hTPCChi2_AC", kTH1F, {{200, 0, 100}}); - histos.add("hITSChi2_AC", "hITSChi2_AC", kTH1F, {{200, 0, 100}}); - histos.add("hTPC_nSigma_AC", "hTPC_nSigma_AC", kTH1F, {pidAxis}); - histos.add("hTOF_nSigma_AC", "hTOF_nSigma_AC", kTH1F, {pidAxis}); - histos.add("hTPC_nSigma_v_pt_AC", "hTPC_nSigma_v_pt_AC", HistType::kTHnSparseD, {pidAxis, PtAxis}); - histos.add("hTOF_nSigma_v_pt_AC", "hTOF_nSigma_v_pt_AC", HistType::kTHnSparseD, {pidAxis, PtAxis}); + + if (doprocessData) { + // Track QA before cuts + // --- Track + if (cFillTrackQA) { + histos.add("QA/QAbefore/Track/TOF_TPC_Map_all", "TOF + TPC Combined PID for All Kaons;{#sigma_{TOF}^{Kaon}};{#sigma_{TPC}^{Kaon}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TOF_Nsigma_all", "TOF NSigma for All Kaons;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPC_Nsigma_all", "TPC NSigma for All Kaons;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/TPConly_Nsigma_all", "TPC NSigma for All Kaons;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + histos.add("QA/QAbefore/Track/dcaZ_all", "DCA_{Z} distribution of selected All Kaons; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm); ", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAbefore/Track/dcaXY_all", "DCA_{XY} momentum distribution of All Kaons; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAbefore/Track/TPC_CR_all", "# TPC Xrows distribution of All Kaons; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAbefore/Track/pT_all", "pT distribution of All Kaons; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAbefore/Track/eta_all", "#eta distribution of All Kaons; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + } + if (cFillMultQA) { + // Multiplicity correlation calibrations + histos.add("MultCalib/centGloPVTrk1", "Centrality vs Global-Tracks", kTHnSparseF, {{110, 0, 110, "Centrality"}, {500, 0, 5000, "Global Tracks"}, {500, 0, 5000, "PV tracks"}}); + histos.add("MultCalib/centGloPVTrk2", "Centrality vs Global-Tracks", kTHnSparseF, {{110, 0, 110, "Centrality"}, {500, 0, 5000, "Global Tracks"}, {500, 0, 5000, "PV tracks"}}); + } + + // PID QA after cuts + // --- Track 1 + histos.add("QA/QAafter/Trk1/TOF_TPC_Map_selected", "TOF + TPC Combined PID for selected Kaon 1;{#sigma_{TOF}^{Kaon}};{#sigma_{TPC}^{Kaon}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAafter/Trk1/TOF_Nsigma_selected", "TOF NSigma for selected Kaon 1;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk1/TPC_Nsigma_selected", "TPC NSigma for selected Kaon 1;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk1/TPC_Nsigma_TPConly_selected", "TPC NSigma for selected Kaon 1;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk1/dcaZ_selected", "DCA_{Z} distribution of selected Kaon 1; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm);", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAafter/Trk1/dcaXY_selected", "DCA_{XY} momentum distribution of selected Kaon 1; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAafter/Trk1/TPC_CR_selected", "# TPC Xrows distribution of selected Kaon 1; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAafter/Trk1/pT_selected", "pT distribution of selected Kaon 1; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAafter/Trk1/eta_selected", "#eta distribution of selected Kaon 1; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + histos.add("QA/QAafter/Trk1/TPC_Signal_selected", "TPC Signal for selected Kaon 1;#it{p} (GeV/#it{c});TPC Signal (A.U.)", {HistType::kTH2F, {axisPt, axisTPCSignal}}); + histos.add("QA/QAafter/Trk1/TPCnclusterPhipr_selected", "TPC ncluster vs phi for selected Kaon 1", kTHnSparseF, {{160, 0, 160, "TPC nCluster"}, {63, 0.0f, 6.28f, "#phi"}}); + + // --- Track 2 + histos.add("QA/QAafter/Trk2/TOF_TPC_Map_selected", "TOF + TPC Combined PID for selected Kaon 2;{#sigma_{TOF}^{Kaon}};{#sigma_{TPC}^{Kaon}}", {HistType::kTH2F, {axisPIDQA, axisPIDQA}}); + histos.add("QA/QAafter/Trk2/TOF_Nsigma_selected", "TOF NSigma for selected Kaon 2;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk2/TPC_Nsigma_selected", "TPC NSigma for selected Kaon 2;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk2/TPC_Nsigma_TPConly_selected", "TPC NSigma for selected Kaon 2;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTH2F, {axisPt, axisPIDQA}}); + histos.add("QA/QAafter/Trk2/dcaZ_selected", "DCA_{Z} distribution of selected Kaon 2; #it{p}_{T} (GeV/#it{c}); DCA_{Z} (cm); ", HistType::kTH2F, {axisPt, axisDCAz}); + histos.add("QA/QAafter/Trk2/dcaXY_selected", "DCA_{XY} momentum distribution of selected Kaon 2; #it{p}_{T} (GeV/#it{c}); DCA_{XY} (cm);", HistType::kTH2F, {axisPt, axisDCAxy}); + histos.add("QA/QAafter/Trk2/TPC_CR_selected", "# TPC Xrows distribution of selected Kaon 2; #it{p}_{T} (GeV/#it{c}); TPC X rows", HistType::kTH2F, {axisPt, axisTPCXrow}); + histos.add("QA/QAafter/Trk2/pT_selected", "pT distribution of selected Kaon 2; #it{p}_{T} (GeV/#it{c}); Counts;", {HistType::kTH1F, {axisPt}}); + histos.add("QA/QAafter/Trk2/eta_selected", "#eta distribution of selected Kaon 2; #eta; Counts;", {HistType::kTH1F, {axisEta}}); + histos.add("QA/QAafter/Trk2/TPC_Signal_selected", "TPC Signal for selected Kaon 2;#it{p} (GeV/#it{c});TPC Signal (A.U.)", {HistType::kTH2F, {axisPt, axisTPCSignal}}); + histos.add("QA/QAafter/Trk2/TPCnclusterPhika_selected", "TPC ncluster vs phi for selected Kaon 2", kTHnSparseF, {{160, 0, 160, "TPC nCluster"}, {63, 0.0f, 6.28f, "#phi"}}); + + // Mass QA 1D for quick check + if (cFill1DQAs) { + histos.add("Result/Data/phi1020invmass", "Invariant mass of #phi(1020) K^{#pm}K^{#mp}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); + histos.add("Result/Data/phi1020invmassLSPP", "Invariant mass of #phi(1020) Like Sign Method K^{#plus}K^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); // K+ + K+ + histos.add("Result/Data/phi1020invmassLSMM", "Invariant mass of #phi(1020) Like Sign Method K^{#minus}K^{#minus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); // K- + K- + } + // eta phi QA + if (cFilldeltaEtaPhiPlots) { + histos.add("QAbefore/deltaEta", "deltaEta of kaon 1 and 2 candidates", HistType::kTH1F, {axisEtaPhi}); + histos.add("QAbefore/deltaPhi", "deltaPhi of kaon 1 and 2 candidates", HistType::kTH1F, {axisEtaPhi}); + + histos.add("QAafter/deltaEta", "deltaEta of kaon 1 and 2 candidates", HistType::kTH1F, {axisEtaPhi}); + histos.add("QAafter/deltaPhi", "deltaPhi of kaon 1 and 2 candidates", HistType::kTH1F, {axisEtaPhi}); + + histos.add("QAafter/PhiafterTrk1", "Phi of kaon 1 candidates", HistType::kTH1F, {axisPhi}); + histos.add("QAafter/PhiafterTrk2", "Phi of kaon 2 candidates", HistType::kTH1F, {axisPhi}); + } + + // 3d histogram + histos.add("Result/Data/h3phi1020invmass", "Invariant mass of #phi(1020) K^{#pm}K^{#mp}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); + histos.add("Result/Data/h3phi1020invmassLSPP", "Invariant mass of #phi(1020) Like Sign Method K^{#plus}K^{#plus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); // K+ + K+ + histos.add("Result/Data/h3phi1020invmassLSMM", "Invariant mass of #phi(1020) Like Sign Method K^{#minus}K^{#minus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); // K- + K- + } + + if (doprocessRotational) { + if (cFill1DQAs) { + histos.add("Result/Data/phi1020InvMassRotation", "Invariant mass of #phi(1020) rotational Method K^{#plus}K^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); // K+ + K+ + } + histos.add("Result/Data/h3phi1020InvMassRotation", "Invariant mass of #phi(1020) rotational Method K^{#plus}K^{#plus}", kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); + } + // Mixed event histograms + if (doprocessME) { + if (cFill1DQAs) { + histos.add("Result/Data/phi1020invmassME", "Invariant mass of #phi(1020) mixed event K^{#pm}K^{#mp}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); + } + histos.add("Result/Data/h3phi1020invmassME", "Invariant mass of #phi(1020) mixed event K^{#pm}K^{#mp}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); + + if (cFilladditionalMEPlots) { + if (cFill1DQAs) { + histos.add("Result/Data/phi1020invmassME_LSPP", "Invariant mass of #phi(1020) Like Sign Method K^{#plus}K^{#plus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); // K+ + K+ + histos.add("Result/Data/phi1020invmassME_LSMM", "Invariant mass of #phi(1020) Like Sign Method K^{#minus}K^{#minus}; Invariant Mass (GeV/#it{c}^2); Counts;", {HistType::kTH1F, {axisMassphi1020}}); // K- + K- + } + histos.add("Result/Data/h3phi1020invmassME_LSPP", "Invariant mass of #phi(1020) mixed event Like Sign Method K^{#plus}K^{#plus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); // K+ + K+ + histos.add("Result/Data/h3phi1020invmassME_LSMM", "Invariant mass of #phi(1020) mixed event Like Sign Method K^{#minus}K^{#minus}", HistType::kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); // K- + K- + } } - // Data histos - histos.add("hUSS", "hUSS", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hLSS", "hLSS", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hUSS_Mix", "hUSS_Mix", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hLSS_Mix", "hLSS_Mix", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - - // MC histos - histos.add("hMC_USS", "hMC_USS", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hMC_LSS", "hMC_LSS", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hMC_USS_Mix", "hMC_USS_Mix", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hMC_LSS_Mix", "hMC_LSS_Mix", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hMC_USS_True", "hMC_USS_True", HistType::kTHnSparseD, {cfg_bins_Cent, MinvAxis, PtAxis}); - histos.add("hMC_Phi_True", "hMC_Phi_True", HistType::kTHnSparseD, {cfg_bins_Cent, PtAxis}); - - // Event Histograms - histos.add("hnEvents", "Event selection decision", kTH1I, {{10, -0.5, 9.5}}); - histos.add("hnEvents_MC", "Event selection decision", kTH1I, {{10, -0.5, 9.5}}); - histos.add("hnEvents_MC_True", "Event selection decision", kTH1I, {{10, -0.5, 9.5}}); + // MC QA + if (doprocessEventFactor) { + histos.add("Event/MultiplicityGenEv", "hMCEventIndices", kTH1D, {axisMult}); + histos.add("Event/MultiplicityRecoEv", "Multiplicity of Reconstructed Events", kTH1D, {axisMult}); + } + + if (doprocessMCGen) { + histos.add("QA/MC/h2GenEtaPt_beforeanycut", " #eta-#it{p}_{T} distribution of Generated #phi(1020); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPt}); + histos.add("QA/MC/h2GenPhiRapidity_beforeanycut", " #phi-y distribution of Generated #phi(1020); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); + histos.add("QA/MC/h2GenEtaPt_afterRapcut", " #phi-#it{p}_{T} distribution of Generated #phi(1020); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPt}); + histos.add("QA/MC/h2GenPhiRapidity_afterRapcut", " #phi-y distribution of Generated #phi(1020); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); + + histos.add("Result/MC/Genphi1020pt", "pT distribution of True MC #phi(1020)0", kTH2F, {axisPt, axisMult}); + } + + if (doprocessMCRec) { + histos.add("QA/MC/h2RecoEtaPt_after", " #eta-#it{p}_{T} distribution of Reconstructed #phi(1020); #eta; #it{p}_{T}; Counts;", HistType::kTHnSparseF, {axisEta, axisPt}); + histos.add("QA/MC/h2RecoPhiRapidity_after", " #phi-y distribution of Reconstructed #phi(1020); #phi; y; Counts;", HistType::kTHnSparseF, {axisPhi, axisRap}); + + histos.add("QA/MC/trkDCAxy_Trk1", "DCAxy distribution of kaon 1 track candidates", HistType::kTHnSparseF, {axisPt, axisDCAxy}); + histos.add("QA/MC/trkDCAxy_Trk2", "DCAxy distribution of kaon 2 track candidates", HistType::kTHnSparseF, {axisPt, axisDCAxy}); + histos.add("QA/MC/trkDCAz_Trk1", "DCAz distribution of kaon 1 track candidates", HistType::kTHnSparseF, {axisPt, axisDCAz}); + histos.add("QA/MC/trkDCAz_Trk2", "DCAz distribution of kaon 2 track candidates", HistType::kTHnSparseF, {axisPt, axisDCAz}); + histos.add("QA/MC/TOF_Nsigma_selected_Trk1", "TOF NSigma for Kaon 1;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TPC_Nsigma_selected_Trk1", "TPC NSigma for Kaon 1;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TOF_Nsigma_selected_Trk2", "TOF NSigma for Kaon 2;#it{p}_{T} (GeV/#it{c});{#sigma_{TOF}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + histos.add("QA/MC/TPC_Nsigma_selected_Trk2", "TPC NSigma for Kaon 2;#it{p}_{T} (GeV/#it{c});{#sigma_{TPC}^{Kaon}};", {HistType::kTHnSparseF, {axisMult, axisPt, axisPIDQA}}); + + histos.add("Result/MC/h3phi1020Recoinvmass", "Invariant mass of Reconstructed MC #phi(1020)0", kTHnSparseF, {axisMult, axisPt, axisMassphi1020}); + } + + if (doprocessSignalLoss) { + histos.add("Result/SignalLoss/GenTruephi1020pt_den", "True phi1020 (den)", kTH2F, {axisPt, axisMult}); + + histos.add("Result/SignalLoss/GenTruephi1020pt_num", "True phi1020 (num)", kTH2F, {axisPt, axisMult}); + } + + // Print output histograms statistics + LOG(info) << "Size of the histograms in Phi1020analysis:"; + histos.print(); } // end of init - Filter collisionFilter = nabs(aod::collision::posZ) <= cfg_Event_VtxCut; - Filter collisionFilter_MC = nabs(aod::mccollision::posZ) <= cfg_Event_VtxCut; - // Filter centralityFilter = nabs(aod::cent::centFT0C) <= cfg_Event_CentralityMax; - Filter acceptanceFilter = (nabs(aod::track::eta) < cfg_Track_MaxEta && nabs(aod::track::pt) >= cfg_Track_MinPt); - using EventCandidates = soa::Filtered>; - using EventCandidates_True = soa::Filtered; - - using TrackCandidates = soa::Filtered>; - using TrackCandidates_MC = soa::Filtered>; - - using BinningTypeVtxCent = ColumnBinningPolicy; - - Partition PosKaon_MC = - (aod::track::signed1Pt > static_cast(0)) && - (!cfg_Track_TPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfg_Track_TPCPID_nSig)); - Partition NegKaon_MC = - (aod::track::signed1Pt < static_cast(0)) && - (!cfg_Track_TPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfg_Track_TPCPID_nSig)); - Partition PosKaon = - (aod::track::signed1Pt > static_cast(0)) && - (!cfg_Track_TPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfg_Track_TPCPID_nSig)); - Partition NegKaon = - (aod::track::signed1Pt < static_cast(0)) && - (!cfg_Track_TPCPID || (nabs(aod::pidtpc::tpcNSigmaKa) <= cfg_Track_TPCPID_nSig)); - - double massKa = o2::constants::physics::MassKPlus; + float massKa = MassKaonCharged; // Centralicity estimator selection template @@ -241,452 +434,764 @@ struct phi1020analysis { return returnValue; } - //***********************************// - // First, we declare some helper functions - template - void fillQA(const bool pass, const objType& obj, const int objecttype = 0) + template + bool isSelected(const Coll& collision, bool fillHist = true) { + auto applyCut = [&](bool enabled, bool condition, int bin) { + if (!enabled) + return true; + if (!condition) + return false; + if (fillHist) + histos.fill(HIST("CollCutCounts"), bin); + return true; + }; - if (objecttype == 1) { - if constexpr (requires { obj.posZ(); }) { - if (!pass) { - histos.fill(HIST("hPosZ_BC"), obj.posZ()); - histos.fill(HIST("hcentFT0_BC"), centEst(obj)); - } else { - histos.fill(HIST("hPosZ_AC"), obj.posZ()); - histos.fill(HIST("hcentFT0_AC"), centEst(obj)); - } - } - } - if constexpr (requires { obj.tpcCrossedRowsOverFindableCls(); }) { - if (objecttype == 2) { - if (!pass) { - histos.fill(HIST("hDCArToPv_BC"), obj.dcaXY()); - histos.fill(HIST("hDCAzToPv_BC"), obj.dcaZ()); - histos.fill(HIST("hIsPrim_BC"), obj.isPrimaryTrack()); - histos.fill(HIST("hIsGood_BC"), obj.isGlobalTrackWoDCA()); - histos.fill(HIST("hIsPrimCont_BC"), obj.isPVContributor()); - histos.fill(HIST("hFindableTPCClusters_BC"), obj.tpcNClsFindable()); - histos.fill(HIST("hFindableTPCRows_BC"), obj.tpcNClsCrossedRows()); - histos.fill(HIST("hClustersVsRows_BC"), obj.tpcCrossedRowsOverFindableCls()); - histos.fill(HIST("hTPCChi2_BC"), obj.tpcChi2NCl()); - } else { - histos.fill(HIST("hDCArToPv_AC"), obj.dcaXY()); - histos.fill(HIST("hDCAzToPv_AC"), obj.dcaZ()); - histos.fill(HIST("hIsPrim_AC"), obj.isPrimaryTrack()); - histos.fill(HIST("hIsGood_AC"), obj.isGlobalTrackWoDCA()); - histos.fill(HIST("hIsPrimCont_AC"), obj.isPVContributor()); - histos.fill(HIST("hFindableTPCClusters_AC"), obj.tpcNClsFindable()); - histos.fill(HIST("hFindableTPCRows_AC"), obj.tpcNClsCrossedRows()); - histos.fill(HIST("hClustersVsRows_AC"), obj.tpcCrossedRowsOverFindableCls()); - histos.fill(HIST("hTPCChi2_AC"), obj.tpcChi2NCl()); - } - } - if (objecttype == 3) { - if (!pass) { - histos.fill(HIST("hTPC_nSigma_BC"), obj.tpcNSigmaKa()); - histos.fill(HIST("hTOF_nSigma_BC"), obj.tofNSigmaKa()); - histos.fill(HIST("hTPC_nSigma_v_pt_BC"), obj.tpcNSigmaKa(), obj.pt()); - histos.fill(HIST("hTOF_nSigma_v_pt_BC"), obj.tofNSigmaKa(), obj.pt()); - } else { - histos.fill(HIST("hTPC_nSigma_AC"), obj.tpcNSigmaKa()); - histos.fill(HIST("hTOF_nSigma_AC"), obj.tofNSigmaKa()); - histos.fill(HIST("hTPC_nSigma_v_pt_AC"), obj.tpcNSigmaKa(), obj.pt()); - histos.fill(HIST("hTOF_nSigma_v_pt_AC"), obj.tofNSigmaKa(), obj.pt()); - } - } - } - }; - //***********************************// + if (fillHist) + histos.fill(HIST("CollCutCounts"), 0); - // evsel - template - std::pair eventSelection(const EventType event, const bool QA) - { + if (!applyCut(true, std::abs(collision.posZ()) <= configEvents.cfgEvtZvtx, 1)) + return false; - if (cfg_Track_CutQA && QA) - fillQA(false, event, 1); - - if (!event.sel8()) - return {false, 1}; - if (std::abs(event.posZ()) > cfg_Event_VtxCut) - return {false, 2}; - if (cfg_Event_Timeframe && (!event.selection_bit(aod::evsel::kNoTimeFrameBorder) || !event.selection_bit(aod::evsel::kNoITSROFrameBorder))) - return {false, 3}; - if (cfg_Event_Timerange && (!event.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard))) - return {false, 4}; - if (cfg_Event_Pileup && (!event.selection_bit(aod::evsel::kNoSameBunchPileup) || !event.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) - return {false, 5}; - if (cfg_Event_Centrality && (centEst(event) > cfg_Event_CentralityMax)) - return {false, 6}; - if (cfg_Event_OccupancyCut && (event.trackOccupancyInTimeRange() > cfg_Event_MaxOccupancy)) - return {false, 7}; - - if (cfg_Track_CutQA && QA) - fillQA(true, event, 1); - - return {true, 8}; - }; - - // tracksel - template - bool trackSelection(const TrackType track, const bool QA) - { - if (cfg_Track_CutQA && QA) - fillQA(false, track, 2); + if (!applyCut(configEvents.cfgEvtTriggerTVXSel, + collision.selection_bit(aod::evsel::kIsTriggerTVX), 2)) + return false; - // basic track cuts - if (track.pt() < cfg_Track_MinPt) + if (!applyCut(configEvents.cfgEvtNoTFBorderCut, + collision.selection_bit(aod::evsel::kNoTimeFrameBorder), 3)) + return false; + + if (!applyCut(configEvents.cfgEvtNoITSROFrameBorderCut, + collision.selection_bit(aod::evsel::kNoITSROFrameBorder), 4)) + return false; + + if (!applyCut(configEvents.cfgEvtIsRCTFlagpassed, rctChecker(collision), 5)) return false; - if (std::abs(track.eta()) > cfg_Track_MaxEta) + + if (!applyCut(configEvents.cfgEvtSel8, collision.sel8(), 6)) return false; - if (std::abs(track.dcaXY()) > cfg_Track_MaxDCArToPVcut) + + if (!applyCut(configEvents.cfgEvtIsINELgt0, collision.isInelGt0(), 7)) return false; - if (std::abs(track.dcaZ()) > cfg_Track_MaxDCAzToPVcut) + + if (fillHist) + histos.fill(HIST("CollCutCounts"), 8); + + return true; + } + + template + bool trackCut(const TrackType track) + { + // basic track cuts + if (configTracks.cDCAr7SigCut && std::abs(track.dcaXY()) > (0.004f + 0.013f / (track.pt()))) // 7 - Sigma cut return false; - if (cfg_Track_PrimaryTrack && !track.isPrimaryTrack()) + if (configTracks.cTPCNClsFound && (track.tpcNClsFound() < configTracks.cMinTPCNClsFound)) return false; - if (track.tpcNClsFindable() < cfg_Track_nFindableTPCClusters) + if (track.tpcNClsCrossedRows() < configTracks.cfgMinCrossedRows) return false; - if (track.tpcNClsCrossedRows() < cfg_Track_nTPCCrossedRows) + if (configTracks.cfgHasTOF && !track.hasTOF()) return false; - if (track.tpcCrossedRowsOverFindableCls() > cfg_Track_nRowsOverFindable) + if (configTracks.cfgPrimaryTrack && !track.isPrimaryTrack()) return false; - if (track.tpcChi2NCl() > cfg_Track_nTPCChi2) + if (configTracks.cfgGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) return false; - if (track.itsChi2NCl() > cfg_Track_nITSChi2) + if (configTracks.cfgPVContributor && !track.isPVContributor()) return false; - if (cfg_Track_ConnectedToPV && !track.isPVContributor()) + if (configTracks.cfgGlobalTrack && !track.isGlobalTrack()) return false; - if (cfg_Track_CutQA && QA) - fillQA(true, track, 2); return true; - }; - - // trackpid + } - template - bool trackPIDKaon(const TrackPID& candidate, const bool QA) + template + bool ptDependentPidKaon(const T& candidate) { - bool tpcPIDPassed{false}, tofPIDPassed{false}; + auto vKaonTPCPIDpTintv = configPID.kaonTPCPIDpTintv.value; + vKaonTPCPIDpTintv.insert(vKaonTPCPIDpTintv.begin(), configTracks.cMinPtcut); + auto vKaonTPCPIDcuts = configPID.kaonTPCPIDcuts.value; + auto vKaonTOFPIDpTintv = configPID.kaonTOFPIDpTintv.value; + auto vKaonTPCTOFCombinedpTintv = configPID.kaonTPCTOFCombinedpTintv.value; + auto vKaonTPCTOFCombinedPIDcuts = configPID.kaonTPCTOFCombinedPIDcuts.value; + auto vKaonTOFPIDcuts = configPID.kaonTOFPIDcuts.value; + + float pt = candidate.pt(); + float ptSwitchToTOF = vKaonTPCPIDpTintv.back(); + float tpcNsigmaKa = candidate.tpcNSigmaKa(); + float tofNsigmaKa = candidate.tofNSigmaKa(); + + bool tpcPIDPassed = false; + + // TPC PID interval-based check + for (size_t i = 0; i < vKaonTPCPIDpTintv.size() - 1; ++i) { + if (pt > vKaonTPCPIDpTintv[i] && pt < vKaonTPCPIDpTintv[i + 1]) { + if (std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts[i]) { + tpcPIDPassed = true; + break; + } + } + } - if (cfg_Track_CutQA && QA) - fillQA(false, candidate, 3); + // TOF bypass option + if (configPID.cByPassTOF) { + return std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts.back(); + } - if (!cfg_Track_TPCPID) { - tpcPIDPassed = true; - } else { - if (std::abs(candidate.tpcNSigmaKa()) < cfg_Track_TPCPID_nSig) - tpcPIDPassed = true; + // Case 1: No TOF and pt ≤ ptSwitch → use TPC-only + if (!candidate.hasTOF() && pt <= ptSwitchToTOF) { + return tpcPIDPassed; } - if (!cfg_Track_TOFPID) { - tofPIDPassed = true; - } else { - if (candidate.hasTOF()) { - if (std::abs(candidate.tofNSigmaKa()) < cfg_Track_TOFPID_nSig) { - tofPIDPassed = true; - } - } else if (!cfg_Track_Hard_TOFPID) { - tofPIDPassed = true; - } - if (!candidate.hasTOF()) { - // std::cout << candidate.tofNSigmaKa() << std::endl; - } + // Case 2: No TOF but pt > ptSwitch → reject + if (!candidate.hasTOF() && pt > ptSwitchToTOF) { + return false; } - if (tpcPIDPassed && tofPIDPassed) { - if (cfg_Track_CutQA && QA) { - fillQA(true, candidate, 3); + + // Case 3: TOF is available → apply TPC+TOF PID logic + if (candidate.hasTOF()) { + if (configPID.cPIDcutType == SquareType) { + // Rectangular cut + for (size_t i = 0; i < vKaonTOFPIDpTintv.size(); ++i) { + if (pt < vKaonTOFPIDpTintv[i]) { + if (std::abs(tofNsigmaKa) < vKaonTOFPIDcuts[i] && + std::abs(tpcNsigmaKa) < vKaonTPCPIDcuts.back()) { + return true; + } + } + } + } else if (configPID.cPIDcutType == CircularType) { + // Circular cut + for (size_t i = 0; i < vKaonTPCTOFCombinedpTintv.size(); ++i) { + if (pt < vKaonTPCTOFCombinedpTintv[i]) { + float combinedSigma2 = tpcNsigmaKa * tpcNsigmaKa + + tofNsigmaKa * tofNsigmaKa; + if (combinedSigma2 < vKaonTPCTOFCombinedPIDcuts[i] * vKaonTPCTOFCombinedPIDcuts[i]) { + return true; + } + } + } } - return true; } + return false; } - template - void TrackSlicing(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool QA, const bool IsMix) + auto static constexpr MinPtforPionRejection = 1.0f; + auto static constexpr MaxPtforPionRejection = 2.0f; + auto static constexpr MaxnSigmaforPionRejection = 2.0f; + + template + bool rejectPion(const T& candidate) { - auto slicedtracks1 = PosKaon->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); - auto slicedtracks2 = NegKaon->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); - auto centrality = centEst(collision1); - for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(slicedtracks1, slicedtracks2))) { - auto [Minv, PhiPt] = minvReconstruction(track1, track2, QA); - if (Minv < 0) + if (candidate.pt() > MinPtforPionRejection && candidate.pt() < MaxPtforPionRejection && !candidate.hasTOF() && candidate.tpcNSigmaPi() < MaxnSigmaforPionRejection) { + return false; + } + return true; + } + + auto static constexpr MaxNophi1020Daughters = 2; + + template + void fillHistograms(const CollisionType& collision, const TracksType& dTracks1, const TracksType& dTracks2) + { + auto centrality = centEst(collision); + + // Multiplicity correlation calibration plots + if (cFillMultQA) { + if constexpr (IsData) { + histos.fill(HIST("MultCalib/centGloPVTrk1"), centrality, dTracks1.size(), collision.multNTracksPV()); + histos.fill(HIST("MultCalib/centGloPVTrk2"), centrality, dTracks2.size(), collision.multNTracksPV()); + } + } + + if (cFilladditionalQAeventPlots) { + if constexpr (IsData) { + histos.fill(HIST("QAevent/hnTrksSameE"), dTracks1.size()); + } else if constexpr (IsMix) { + histos.fill(HIST("QAevent/hVertexZMixedE"), collision.posZ()); + histos.fill(HIST("QAevent/hMultiplicityPercentMixedE"), centrality); + histos.fill(HIST("QAevent/hnTrksMixedE"), dTracks1.size()); + } + } + // LOG(info) << "After pass, Collision index:" << collision.index() << "multiplicity: " << collision.centFT0M() << endl; + + LorentzVectorPtEtaPhiMass lDecayDaughter1, lDecayDaughter2, lResonance, ldaughterRot, lResonanceRot; + + for (const auto& [trk1, trk2] : combinations(CombinationsFullIndexPolicy(dTracks1, dTracks2))) { + // Full index policy is needed to consider all possible combinations + if (trk1.index() == trk2.index()) + continue; // We need to run (0,1), (1,0) pairs as well. but same id pairs are not needed. + + // apply the track cut + if (!trackCut(trk1) || !trackCut(trk2)) continue; - double conjugate = track1.sign() * track2.sign(); - if (!IsMix) { - if (conjugate < 0) { - histos.fill(HIST("hUSS"), centrality, Minv, PhiPt); - } else if (conjugate > 0) { - histos.fill(HIST("hLSS"), centrality, Minv, PhiPt); + + //// Initialize variables + // Trk1: Proton + auto isTrk1hasTOF = trk1.hasTOF(); + auto trk1pt = trk1.pt(); + auto trk1eta = trk1.eta(); + auto trk1phi = trk1.phi(); + auto trk1NSigmaTPC = trk1.tpcNSigmaKa(); + auto trk1NSigmaTOF = (isTrk1hasTOF) ? trk1.tofNSigmaKa() : -999.0f; + + // Trk2: Kaon + auto isTrk2hasTOF = trk2.hasTOF(); + auto trk2pt = trk2.pt(); + auto trk2eta = trk2.eta(); + auto trk2phi = trk2.phi(); + auto trk2NSigmaTPC = trk2.tpcNSigmaKa(); + auto trk2NSigmaTOF = (isTrk2hasTOF) ? trk2.tofNSigmaKa() : -999.0f; + + auto deltaEta = 0.0f; + auto deltaPhi = 0.0f; + + if (cfgUseDeltaEtaPhiCuts) { + deltaEta = std::abs(trk1eta - trk2eta); + deltaPhi = std::abs(trk1phi - trk2phi); + deltaPhi = (deltaPhi > o2::constants::math::PI) ? (o2::constants::math::TwoPI - deltaPhi) : deltaPhi; + if (deltaEta >= cMaxDeltaEtaCut) + continue; + if (deltaPhi >= cMaxDeltaPhiCut) + continue; + } + + //// QA plots before the selection + // --- Track QA all + if constexpr (IsData) { + if (cFillTrackQA) { + histos.fill(HIST("QA/QAbefore/Track/TPC_Nsigma_all"), centrality, trk2pt, trk2NSigmaTPC); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TOF_Nsigma_all"), centrality, trk2pt, trk2NSigmaTOF); + histos.fill(HIST("QA/QAbefore/Track/TOF_TPC_Map_all"), trk2NSigmaTOF, trk2NSigmaTPC); + } + if (!isTrk2hasTOF) { + histos.fill(HIST("QA/QAbefore/Track/TPConly_Nsigma_all"), trk2pt, trk2NSigmaTPC); + } + + histos.fill(HIST("QA/QAbefore/Track/dcaZ_all"), trk1pt, trk1.dcaZ()); + histos.fill(HIST("QA/QAbefore/Track/dcaXY_all"), trk1pt, trk1.dcaXY()); + histos.fill(HIST("QA/QAbefore/Track/TPC_CR_all"), trk1pt, trk1.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAbefore/Track/pT_all"), trk1pt); + histos.fill(HIST("QA/QAbefore/Track/eta_all"), trk1eta); } - } else { - if (conjugate < 0) { - histos.fill(HIST("hUSS_Mix"), centrality, Minv, PhiPt); - } else if (conjugate > 0) { - histos.fill(HIST("hLSS_Mix"), centrality, Minv, PhiPt); + if (cFilldeltaEtaPhiPlots) { + histos.fill(HIST("QAbefore/deltaEta"), deltaEta); + histos.fill(HIST("QAbefore/deltaPhi"), deltaPhi); } } - } - } // TrackSlicing - template - void TrackSlicing_MC(const CollisionType& collision1, const TracksType&, const CollisionType& collision2, const TracksType&, const bool QA, const bool IsMix) - { - auto slicedtracks1 = PosKaon_MC->sliceByCached(aod::track::collisionId, collision1.globalIndex(), cache); - auto slicedtracks2 = NegKaon_MC->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); - auto centrality = centEst(collision1); - for (auto& [track1, track2] : combinations(o2::soa::CombinationsFullIndexPolicy(slicedtracks1, slicedtracks2))) { - auto [Minv, PhiPt] = minvReconstruction(track1, track2, QA); - if (Minv < 0) + //// Apply the pid selection + if (crejectPion && rejectPion(trk2)) // to remove pion contamination from the kaon track continue; - double conjugate = track1.sign() * track2.sign(); - if (!IsMix) { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS"), centrality, Minv, PhiPt); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS"), centrality, Minv, PhiPt); + + if (!ptDependentPidKaon(trk1) || !ptDependentPidKaon(trk2)) + continue; + + //// QA plots after the selection + if constexpr (IsData) { // --- PID QA Proton + histos.fill(HIST("QA/QAafter/Trk1/TPC_Nsigma_selected"), centrality, trk1pt, trk1NSigmaTPC); + histos.fill(HIST("QA/QAafter/Trk1/TPC_Signal_selected"), trk1.tpcInnerParam(), trk1.tpcSignal()); + if (isTrk1hasTOF) { + histos.fill(HIST("QA/QAafter/Trk1/TOF_Nsigma_selected"), centrality, trk1pt, trk1NSigmaTOF); + histos.fill(HIST("QA/QAafter/Trk1/TOF_TPC_Map_selected"), trk1NSigmaTOF, trk1NSigmaTPC); } - } else { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS_Mix"), centrality, Minv, PhiPt); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS_Mix"), centrality, Minv, PhiPt); + if (!isTrk1hasTOF) { + histos.fill(HIST("QA/QAafter/Trk1/TPC_Nsigma_TPConly_selected"), trk1pt, trk1NSigmaTPC); + } + histos.fill(HIST("QA/QAafter/Trk1/dcaZ_selected"), trk1pt, trk1.dcaZ()); + histos.fill(HIST("QA/QAafter/Trk1/dcaXY_selected"), trk1pt, trk1.dcaXY()); + histos.fill(HIST("QA/QAafter/Trk1/TPC_CR_selected"), trk1pt, trk1.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAafter/Trk1/pT_selected"), trk1pt); + histos.fill(HIST("QA/QAafter/Trk1/eta_selected"), trk1eta); + histos.fill(HIST("QA/QAafter/Trk1/TPCnclusterPhipr_selected"), trk1.tpcNClsFound(), trk1phi); + + // --- PID QA Kaon + histos.fill(HIST("QA/QAafter/Trk2/TPC_Nsigma_selected"), centrality, trk2pt, trk2NSigmaTPC); + histos.fill(HIST("QA/QAafter/Trk2/TPC_Signal_selected"), trk2.tpcInnerParam(), trk2.tpcSignal()); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/QAafter/Trk2/TOF_Nsigma_selected"), centrality, trk2pt, trk2NSigmaTOF); + histos.fill(HIST("QA/QAafter/Trk2/TOF_TPC_Map_selected"), trk2NSigmaTOF, trk2NSigmaTPC); + } + if (!isTrk2hasTOF) { + histos.fill(HIST("QA/QAafter/Trk2/TPC_Nsigma_TPConly_selected"), trk2pt, trk2NSigmaTPC); + } + histos.fill(HIST("QA/QAafter/Trk2/dcaZ_selected"), trk2pt, trk2.dcaZ()); + histos.fill(HIST("QA/QAafter/Trk2/dcaXY_selected"), trk2pt, trk2.dcaXY()); + histos.fill(HIST("QA/QAafter/Trk2/TPC_CR_selected"), trk2pt, trk2.tpcNClsCrossedRows()); + histos.fill(HIST("QA/QAafter/Trk2/pT_selected"), trk2pt); + histos.fill(HIST("QA/QAafter/Trk2/eta_selected"), trk2eta); + histos.fill(HIST("QA/QAafter/Trk2/TPCnclusterPhika_selected"), trk2.tpcNClsFound(), trk2phi); + + if (cFilldeltaEtaPhiPlots) { + histos.fill(HIST("QAafter/PhiafterTrk1"), trk1phi); + histos.fill(HIST("QAafter/PhiafterTrk2"), trk2phi); + histos.fill(HIST("QAafter/deltaEta"), deltaEta); + histos.fill(HIST("QAafter/deltaPhi"), deltaPhi); } } - // now we do mc true - if (!track1.has_mcParticle() || !track2.has_mcParticle()) + + //// Resonance reconstruction + lDecayDaughter1 = LorentzVectorPtEtaPhiMass(trk1pt, trk1eta, trk1phi, massKa); + lDecayDaughter2 = LorentzVectorPtEtaPhiMass(trk2pt, trk2eta, trk2phi, massKa); + + // Apply kinematic opening angle cut + if (cUseOpeningAngleCut) { + auto v1 = lDecayDaughter1.Vect(); + auto v2 = lDecayDaughter2.Vect(); + float alpha = std::acos(v1.Dot(v2) / (v1.R() * v2.R())); + if (alpha > cMinOpeningAngle && alpha < cMaxOpeningAngle) + continue; + } + + lResonance = lDecayDaughter1 + lDecayDaughter2; + + auto resonanceMass = lResonance.M(); + auto resonancePt = lResonance.Pt(); + auto resonanceRapidity = lResonance.Rapidity(); + + // Rapidity cut + if (std::abs(resonanceRapidity) > configTracks.cfgCutRapidity) continue; - auto part1 = track1.mcParticle(); - auto part2 = track2.mcParticle(); - if (std::fabs(part1.pdgCode()) != 321) - continue; // Not Kaon - if (std::fabs(part2.pdgCode()) != 321) - continue; // Not Kaon - - if (!part1.has_mothers()) - continue; // Not decaying Kaon - if (!part2.has_mothers()) - continue; // Not decaying Kaon - - std::vector mothers1{}; - std::vector mothers1PDG{}; - for (auto& part1_mom : part1.template mothers_as()) { - mothers1.push_back(part1_mom.globalIndex()); - mothers1PDG.push_back(part1_mom.pdgCode()); + + if (cfgUseCutsOnMother) { + if (resonancePt >= cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (resonanceMass >= cMaxMinvMotherCut) // excluding candidates in overflow + continue; } - std::vector mothers2{}; - std::vector mothers2PDG{}; - for (auto& part2_mom : part2.template mothers_as()) { - mothers2.push_back(part2_mom.globalIndex()); - mothers2PDG.push_back(part2_mom.pdgCode()); + if (cFilladditionalQAeventPlots) { + if constexpr (IsData) { + histos.fill(HIST("QAevent/hPairsCounterSameE"), 1.0f); + } else if (IsMix) { + histos.fill(HIST("QAevent/hPairsCounterMixedE"), 1.0f); + } } - if (mothers1PDG[0] != 333) - continue; // mother not phi - if (mothers2PDG[0] != 333) - continue; // mother not phi + //// Un-like sign pair only + if (trk1.sign() * trk2.sign() < 0) { + if constexpr (IsRot) { + for (int i = 0; i < configBkg.cNofRotations; i++) { + float theta = rn->Uniform(o2::constants::math::PI - o2::constants::math::PI / configBkg.rotationalcut, o2::constants::math::PI + o2::constants::math::PI / configBkg.rotationalcut); + if (configBkg.cfgRotTrk1) { + ldaughterRot = LorentzVectorPtEtaPhiMass(trk1pt, trk1eta, trk1phi + theta, massKa); + lResonanceRot = ldaughterRot + lDecayDaughter2; + } else { + ldaughterRot = LorentzVectorPtEtaPhiMass(trk2pt, trk2eta, trk2phi + theta, massKa); + lResonanceRot = lDecayDaughter1 + ldaughterRot; + } + auto resonanceRotMass = lResonanceRot.M(); + auto resonanceRotPt = lResonanceRot.Pt(); + + // Rapidity cut + if (std::abs(lResonanceRot.Rapidity()) >= configTracks.cfgCutRapidity) + continue; + + if (cfgUseCutsOnMother) { + if (resonanceRotPt >= cMaxPtMotherCut) // excluding candidates in overflow + continue; + if (resonanceRotMass >= cMaxMinvMotherCut) // excluding candidates in overflow + continue; + } + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020InvMassRotation"), resonanceRotMass); + } + histos.fill(HIST("Result/Data/h3phi1020InvMassRotation"), centrality, resonanceRotPt, resonanceRotMass); + } + } + if constexpr (IsData) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmass"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmass"), centrality, resonancePt, resonanceMass); + } else if (IsMix) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmassME"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmassME"), centrality, resonancePt, resonanceMass); + } - if (mothers1[0] != mothers2[0]) - continue; // Kaons not from the same phi + // MC + if constexpr (IsMC) { + // now we do mc true + // ------ Temporal lambda function to prevent error in build + auto getMothersIndeces = [&](auto const& theMcParticle) { + std::vector lMothersIndeces{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother index lMother: %d", lMother.globalIndex()); + lMothersIndeces.push_back(lMother.globalIndex()); + } + return lMothersIndeces; + }; + auto getMothersPDGCodes = [&](auto const& theMcParticle) { + std::vector lMothersPDGs{}; + for (auto const& lMother : theMcParticle.template mothers_as()) { + LOGF(debug, " mother pdgcode lMother: %d", lMother.pdgCode()); + lMothersPDGs.push_back(lMother.pdgCode()); + } + return lMothersPDGs; + }; + // ------ + std::vector motherstrk1 = {-1, -1}; + std::vector mothersPDGtrk1 = {-1, -1}; + + std::vector motherstrk2 = {-1, -1}; + std::vector mothersPDGtrk2 = {-1, -1}; + + // + // Get the MC particle + const auto& mctrk1 = trk1.mcParticle(); + if (mctrk1.has_mothers()) { + motherstrk1 = getMothersIndeces(mctrk1); + mothersPDGtrk1 = getMothersPDGCodes(mctrk1); + } + while (motherstrk1.size() > MaxNophi1020Daughters) { + motherstrk1.pop_back(); + mothersPDGtrk1.pop_back(); + } - histos.fill(HIST("hMC_USS_True"), centrality, Minv, PhiPt); - } - } // TrackSlicing + const auto& mctrk2 = trk2.mcParticle(); + if (mctrk2.has_mothers()) { + motherstrk2 = getMothersIndeces(mctrk2); + mothersPDGtrk2 = getMothersPDGCodes(mctrk2); + } + while (motherstrk2.size() > MaxNophi1020Daughters) { + motherstrk2.pop_back(); + mothersPDGtrk2.pop_back(); + } - // Invariant mass - template - std::pair minvReconstruction(const TracksType& trk1, const TracksType& trk2, const bool QA) - { - TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - //==================================================== + if (std::abs(mctrk1.pdgCode()) != kKPlus || std::abs(mctrk2.pdgCode()) != kKPlus) + continue; - if (!trackSelection(trk1, QA) || !trackSelection(trk2, false)) - return {-1.0, -1.0}; + if (motherstrk1[0] != motherstrk2[0]) // Same mother + continue; - if (cfg_Track_Explicit_PID) { - if (!trackPIDKaon(trk1, QA) || !trackPIDKaon(trk2, false)) - return {-1.0, -1.0}; - } + if (std::abs(mothersPDGtrk1[0]) != Pdg::kPhi) + continue; - if (trk1.globalIndex() >= trk2.globalIndex()) - return {-1.0, -1.0}; + // LOGF(info, "mother trk1 id: %d, mother trk1: %d, trk1 id: %d, trk1 pdgcode: %d, mother trk2 id: %d, mother trk2: %d, trk2 id: %d, trk2 pdgcode: %d", motherstrk1[0], mothersPDGtrk1[0], trk1.globalIndex(), mctrk1.pdgCode(), motherstrk2[0], mothersPDGtrk2[0], trk2.globalIndex(), mctrk2.pdgCode()); - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); - lResonance = lDecayDaughter1 + lDecayDaughter2; + histos.fill(HIST("QA/MC/h2RecoEtaPt_after"), lResonance.Eta(), resonancePt); + histos.fill(HIST("QA/MC/h2RecoPhiRapidity_after"), lResonance.Phi(), resonanceRapidity); - return {lResonance.M(), lResonance.Pt()}; - } // MinvReconstruction + // Track selection check. + histos.fill(HIST("QA/MC/trkDCAxy_Trk1"), trk1pt, trk1.dcaXY()); + histos.fill(HIST("QA/MC/trkDCAxy_Trk2"), trk2pt, trk2.dcaXY()); + histos.fill(HIST("QA/MC/trkDCAz_Trk1"), trk1pt, trk1.dcaZ()); + histos.fill(HIST("QA/MC/trkDCAz_Trk2"), trk2pt, trk2.dcaZ()); - //***************// - // DATA - //***************// + histos.fill(HIST("QA/MC/TPC_Nsigma_selected_Trk1"), centrality, trk1pt, trk1NSigmaTPC); + if (isTrk1hasTOF) { + histos.fill(HIST("QA/MC/TOF_Nsigma_selected_Trk1"), centrality, trk1pt, trk1NSigmaTOF); + } + histos.fill(HIST("QA/MC/TPC_Nsigma_selected_Trk2"), centrality, trk2pt, trk2NSigmaTPC); + if (isTrk2hasTOF) { + histos.fill(HIST("QA/MC/TOF_Nsigma_selected_Trk2"), centrality, trk2pt, trk2NSigmaTOF); + } - int nEvents = 0; - void processSameEvent(EventCandidates::iterator const& collision, TrackCandidates const& tracks) - { - if (cDebugLevel > 0) { - ++nEvents; - if (nEvents % 10000 == 0) { - std::cout << "Processed Data Events: " << nEvents << std::endl; + // MC histograms + histos.fill(HIST("Result/MC/h3phi1020Recoinvmass"), centrality, resonancePt, resonanceMass); + } + } else { + if constexpr (IsData) { + // Like sign pair ++ + if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmassLSPP"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmassLSPP"), centrality, resonancePt, resonanceMass); + } else { // Like sign pair -- + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmassLSMM"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmassLSMM"), centrality, resonancePt, resonanceMass); + } + } else if (IsMix) { + if (cFilladditionalMEPlots) { + // Like sign pair ++ + if (trk1.sign() > 0) { + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmassME_LSPP"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmassME_LSPP"), centrality, resonancePt, resonanceMass); + } else { // Like sign pair -- + if (cFill1DQAs) { + histos.fill(HIST("Result/Data/phi1020invmassME_LSMM"), resonanceMass); + } + histos.fill(HIST("Result/Data/h3phi1020invmassME_LSMM"), centrality, resonancePt, resonanceMass); + } + } + } } } + } + + void processData(EventCandidates::iterator const& collision, + TrackCandidates const& tracks) + { + if (!isSelected(collision)) // Default event selection + return; + + auto centrality = centEst(collision); + + histos.fill(HIST("Event/posZ"), collision.posZ()); + histos.fill(HIST("Event/centFT0M"), centrality); + + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Phi1020analysis, processData, "Process Event for data without partition", false); - auto [goodEv, code] = eventSelection(collision, true); - histos.fill(HIST("hnEvents"), code); - if (!goodEv) + void processRotational(EventCandidates::iterator const& collision, TrackCandidates const& tracks) + { + if (!isSelected(collision, false)) // Default event selection return; - TrackSlicing(collision, tracks, collision, tracks, true, false); - } // end of process + if (!collision.isInelGt0()) // <-- + return; - PROCESS_SWITCH(phi1020analysis, processSameEvent, "Process Same events", true); + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Phi1020analysis, processRotational, "Process Rotational Background", false); - //***************// - // DATA (MIX) - //***************// + // Processing Event Mixing + using BinningTypeVtxZT0M = ColumnBinningPolicy; - int nEvents_Mix = 0; - void processMixedEvent(EventCandidates const& collisions, TrackCandidates const& tracks) + void processME(EventCandidates const& collision, + TrackCandidates const& tracks) { auto tracksTuple = std::make_tuple(tracks); - BinningTypeVtxCent colBinning{{cfg_bins_MixVtx, cfg_bins_MixMult}, true}; - SameKindPair pairs{colBinning, cfg_Mix_NMixedEvents, -1, collisions, tracksTuple, &cache}; + + BinningTypeVtxZT0M colBinning{{configBkg.cfgVtxBins, configBkg.cfgMultPercentileBins}, true}; + SameKindPair pairs{colBinning, configBkg.nEvtMixing, -1, collision, tracksTuple, &cache}; // -1 is the number of the bin to skip for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - if (cDebugLevel > 0) { - ++nEvents_Mix; - if (nEvents_Mix % 10000 == 0) { - std::cout << "Processed Mixed Events: " << nEvents_Mix << std::endl; - } - } - auto [goodEv1, code1] = eventSelection(collision1, false); - auto [goodEv2, code2] = eventSelection(collision2, false); - if (!goodEv1 || !goodEv2) + // LOGF(info, "Mixed event collisions: (%d, %d)", collision1.globalIndex(), collision2.globalIndex()); + + // for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) { + // LOGF(info, "Mixed event tracks pair: (%d, %d) from events (%d, %d)", t1.index(), t2.index(), collision1.index(), collision2.index()); + // } + + if (!isSelected(collision1, false)) // Default event selection continue; - TrackSlicing(collision1, tracks1, collision2, tracks2, false, true); - } // mixing - } // end of process - PROCESS_SWITCH(phi1020analysis, processMixedEvent, "Process Mixed events", false); - //***************// - // RECONSTRUCTED MC - //***************// + if (!isSelected(collision2, false)) // Default event selection + continue; - int nEvents_MC = 0; - void processSameEvent_MC(EventCandidates::iterator const& collision, TrackCandidates_MC const& tracks, aod::McParticles const&) - { - if (cDebugLevel > 0) { - ++nEvents_MC; - if (nEvents_MC % 10000 == 0) { - std::cout << "Processed MC (REC) Events: " << nEvents_MC << std::endl; + if (!collision1.isInelGt0()) // <-- + continue; + + if (!collision2.isInelGt0()) // <-- + continue; + + if (cFilladditionalQAeventPlots) { + // Fill histograms for the characteristics of the *mixed* events (collision1 and collision2) + // This will show the distribution of events that are actually being mixed. + if (cFill1DQAs) { + histos.fill(HIST("QAevent/hMixPool_VtxZ"), collision1.posZ()); + histos.fill(HIST("QAevent/hMixPool_Multiplicity"), collision1.centFT0M()); + } + histos.fill(HIST("QAevent/hMixPool_VtxZ_vs_Multiplicity"), collision1.posZ(), collision1.centFT0M()); + + // You might also want to fill for collision2 if you want to see both partners' distributions + // histos.fill(HIST("QAevent/hMixPool_VtxZ"), collision2.posZ()); + // histos.fill(HIST("QAevent/hMixPool_Multiplicity"), collision2.centFT0M()); + // histos.fill(HIST("QAevent/hMixPool_VtxZ_vs_Multiplicity"), collision2.posZ(), collision2.centFT0M()); } + fillHistograms(collision1, tracks1, tracks2); } + } + PROCESS_SWITCH(Phi1020analysis, processME, "Process EventMixing light without partition", false); - auto [goodEv, code] = eventSelection(collision, true); - histos.fill(HIST("hnEvents_MC"), code); - if (!goodEv) + void processMCRec(MCEventCandidates::iterator const& collision, + aod::McCollisions const&, + MCTrackCandidates const& tracks, aod::McParticles const&) + { + if (!collision.has_mcCollision()) return; - TrackSlicing_MC(collision, tracks, collision, tracks, true, false); - } // end of process - PROCESS_SWITCH(phi1020analysis, processSameEvent_MC, "Process Same events (MC)", true); + if (!isSelected(collision)) + return; - //***************// - // RECONSTRUCTED MC (MIX) - //***************// + auto centrality = centEst(collision); - int nEvents_MC_Mix = 0; - void processMixedEvent_MC(EventCandidates const& collisions, TrackCandidates_MC const& tracks, aod::McParticles const&) + histos.fill(HIST("Event/posZ"), collision.posZ()); + histos.fill(HIST("Event/centFT0M"), centrality); + + fillHistograms(collision, tracks, tracks); + } + PROCESS_SWITCH(Phi1020analysis, processMCRec, "Process Event for MC Rec without partition", false); + + Partition selectedMCParticles = (nabs(aod::mcparticle::pdgCode) == static_cast(Pdg::kPhi)); // Lambda(1520) + + void processMCGen(MCEventCandidates::iterator const& collision, aod::McCollisions const&, aod::McParticles const& mcParticles) { - auto tracksTuple = std::make_tuple(tracks); - BinningTypeVtxCent colBinning{{cfg_bins_MixVtx, cfg_bins_MixMult}, true}; - SameKindPair pairs{colBinning, cfg_Mix_NMixedEvents, -1, collisions, tracksTuple, &cache}; - for (const auto& [collision1, tracks1, collision2, tracks2] : pairs) { - if (cDebugLevel > 0) { - ++nEvents_MC_Mix; - if (nEvents_MC_Mix % 10000 == 0) { - std::cout << "Processed Mixed Events: " << nEvents_MC_Mix << std::endl; - } + if (!collision.has_mcCollision()) + return; + + bool isInAfterAllCuts = isSelected(collision, false); + + auto mcPartsAll = mcParticles.sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + // bool isTrueINELgt0 = pwglf::isINELgt0mc(mcPartsAll, pdg); + bool isTrueINELgt0 = collision.isInelGt0(); // <-- + + auto centrality = centEst(collision); + + auto mcParts = selectedMCParticles->sliceBy(perMcCollision, collision.mcCollision().globalIndex()); + + // Not related to the real collisions + for (const auto& part : mcParts) { // loop over all Lambda(1520) particles + + std::vector daughterPDGs; + if (part.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(part.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(part.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } + + if (std::abs(daughterPDGs[0]) != PDG_t::kKPlus || std::abs(daughterPDGs[1]) != PDG_t::kKPlus) { // At least one decay to Kaon + continue; } - auto [goodEv1, code1] = eventSelection(collision1, false); - auto [goodEv2, code2] = eventSelection(collision2, false); - if (!goodEv1 || !goodEv2) + + // LOGF(info, "Part PDG: %d", part.pdgCode(), "DAU_ID1: %d", pass1, "DAU_ID2: %d", pass2); + + histos.fill(HIST("QA/MC/h2GenEtaPt_beforeanycut"), part.eta(), part.pt()); + histos.fill(HIST("QA/MC/h2GenPhiRapidity_beforeanycut"), part.phi(), part.y()); + + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) // rapidity cut continue; - TrackSlicing_MC(collision1, tracks1, collision2, tracks2, false, true); - } // mixing - } // end of process - PROCESS_SWITCH(phi1020analysis, processMixedEvent_MC, "Process Mixed events (MC)", false); - //***************// - // GENERATED MC - //***************// + if (cfgUseDaughterEtaCutMC) { + for (auto const& daughters : part.daughters_as()) { + if (std::fabs(daughters.eta()) > configTracks.cfgCutEta) + continue; // eta cut for daughters + } // loop over daughters + } + + histos.fill(HIST("QA/MC/h2GenEtaPt_afterRapcut"), part.eta(), part.pt()); + histos.fill(HIST("QA/MC/h2GenPhiRapidity_afterRapcut"), part.phi(), part.y()); - int nEvents_True = 0; - void processParticles(EventCandidates_True::iterator const& collision, soa::SmallGroups> const& recocolls, aod::McParticles const& particles) + if (isInAfterAllCuts && isTrueINELgt0) // after all event selection && INEL>0 + histos.fill(HIST("Result/MC/Genphi1020pt"), part.pt(), centrality); + } + } + PROCESS_SWITCH(Phi1020analysis, processMCGen, "Process Event for MC only", false); + + void processEventFactor(MCEventCandidates const& collisions, soa::Join const& mcCollisions, aod::McParticles const& mcParticles) { - if (cDebugLevel > 0) { - ++nEvents_True; - if (nEvents_True % 10000 == 0) { - std::cout << "Processed MC (GEN) Events: " << nEvents_True << std::endl; + // Loop on reconstructed collisions + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; } - } + const auto& mcCollision = collision.mcCollision_as>(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheMC); - if (fabs(collision.posZ()) > cfg_Event_VtxCut) - return; + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + bool isInAfterAllCuts = isSelected(collision, false); - if (recocolls.size() <= 0) { // not reconstructed - if (cfg_Force_GenReco) { - return; - } + float centrality = mcCollision.centFT0M(); + + if (isTrueINELgt0 && isInAfterAllCuts) + histos.fill(HIST("Event/MultiplicityRecoEv"), centrality); } - double centrality = -1; - for (auto& recocoll : recocolls) { // poorly reconstructed - centrality = centEst(recocoll); - auto [goodEv, code] = eventSelection(recocoll, false); - histos.fill(HIST("hnEvents_MC_True"), code); - if (!goodEv) - return; + // Loop on generated collisions to fill the event factor for the INEL>0 correction + for (const auto& mccolls : mcCollisions) { + float centrality = mccolls.centFT0M(); + bool inVtx10 = std::abs(mccolls.posZ()) <= configEvents.cfgEvtZvtx; + + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mccolls.globalIndex(), cacheMC); + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); // QA for Trigger efficiency + + if (inVtx10 && isTrueINELgt0) + histos.fill(HIST("Event/MultiplicityGenEv"), centrality); } + } + PROCESS_SWITCH(Phi1020analysis, processEventFactor, "Process Event factor", false); + + void processSignalLoss(MCEventCandidates const& collisions, soa::Join const& mcCollisions, aod::McParticles const& mcParticles) + { + // Loop on reconstructed collisions + for (const auto& collision : collisions) { + if (!collision.has_mcCollision()) { + continue; + } + const auto& mcCollision = collision.mcCollision_as>(); + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheMC); + + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + bool isInAfterAllCuts = isSelected(collision, false); + bool inVtx10 = std::abs(mcCollision.posZ()) <= configEvents.cfgEvtZvtx; - for (auto& particle : particles) { - if (particle.pdgCode() != 333) + float centrality = mcCollision.centFT0M(); + + // ===== NUM ===== + if (!(inVtx10 && isTrueINELgt0)) continue; - if (std::fabs(particle.eta()) > cfg_Track_MaxEta) + + if (!isInAfterAllCuts) continue; - if (cfg_Force_BR) { - bool baddecay = false; - for (auto& phidaughter : particle.daughters_as()) { - if (std::fabs(phidaughter.pdgCode()) != 321) { - baddecay = true; - break; - } - if (cfg_Force_Kaon_Acceptence) { - if (std::fabs(phidaughter.eta()) > cfg_Track_MaxEta) { - baddecay = true; - break; - } - } - } // loop over daughters + for (const auto& part : particlesInCollision) { + + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) + continue; - if (baddecay) + float pt = part.pt(); + + // phi + if (std::abs(part.pdgCode()) == Pdg::kPhi) + histos.fill(HIST("Result/SignalLoss/GenTruephi1020pt_num"), pt, centrality); + } + } + + // Loop on generated collisions to fill the event factor for the INEL>0 correction + for (const auto& mccolls : mcCollisions) { + float centrality = mccolls.centFT0M(); + + bool inVtx10 = std::abs(mccolls.posZ()) <= configEvents.cfgEvtZvtx; + + const auto& particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mccolls.globalIndex(), cacheMC); + bool isTrueINELgt0 = pwglf::isINELgt0mc(particlesInCollision, pdg); + + if (!(inVtx10 && isTrueINELgt0)) + continue; + + for (const auto& part : particlesInCollision) { + + if (cUseRapcutMC && std::abs(part.y()) > configTracks.cfgCutRapidity) continue; - } // enforce BR restriction - histos.fill(HIST("hMC_Phi_True"), centrality, particle.pt()); - } // loop over particles + // ========================= + // ===== PHI(1020) ====== + // ========================= + if (std::abs(part.pdgCode()) == Pdg::kPhi) { + + std::vector daughterPDGs; + if (part.has_daughters()) { + auto daughter01 = mcParticles.rawIteratorAt(part.daughtersIds()[0] - mcParticles.offset()); + auto daughter02 = mcParticles.rawIteratorAt(part.daughtersIds()[1] - mcParticles.offset()); + daughterPDGs = {daughter01.pdgCode(), daughter02.pdgCode()}; + } else { + daughterPDGs = {-1, -1}; + } - } // end of process - PROCESS_SWITCH(phi1020analysis, processParticles, "Process Particles", false); + if (std::abs(daughterPDGs[0]) != PDG_t::kKPlus || std::abs(daughterPDGs[1]) != PDG_t::kKPlus) { // At least one decay to Kaon + continue; + } + histos.fill(HIST("Result/SignalLoss/GenTruephi1020pt_den"), part.pt(), centrality); + } + } + } + } + PROCESS_SWITCH(Phi1020analysis, processSignalLoss, "Process SignalLoss", false); }; // end of main struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; };