diff --git a/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx b/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx index 61ea4cb13cf..ec76c05b8e1 100644 --- a/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx +++ b/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx @@ -21,9 +21,9 @@ #include "PWGLF/Utils/inelGt.h" #include "Common/CCDB/EventSelectionParams.h" +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -48,6 +48,7 @@ #include #include #include +#include #include #include @@ -57,7 +58,7 @@ using namespace o2::framework::expressions; struct Lambdak0sflattenicity { // Histograms are defined with HistogramRegistry - Service pdg; + Service pdg{}; HistogramRegistry rEventSelection{"eventSelection", {}, OutputObjHandlingPolicy::AnalysisObject, @@ -101,7 +102,7 @@ struct Lambdak0sflattenicity { true}; static constexpr int kNEstimators = 8; - // forward-detector segmentation, used by the channel-mapping helpers below + // forward detector segmentation static constexpr int kNChannelsPerT0Sector = 4; static constexpr int kNSectorsT0A = 24; static constexpr int kNSectorsT0C = 28; @@ -109,165 +110,184 @@ struct Lambdak0sflattenicity { static constexpr int kNFV0EtaRings = 5; static constexpr int kOuterFV0RingIndex = kNFV0EtaRings - 1; static constexpr int kNSectorsFV0OuterRing = 16; - // charge threshold used to decide whether a generated particle is charged + // TParticlePDG::Charge() is in units of e/3 static constexpr float kMinCharge = 0.01f; - // rapidity window applied to generated particles + // rapidity window for generated particles static constexpr float kMcRapidityWindow = 0.5f; - static constexpr std::string_view kHEst[kNEstimators] = { + static constexpr std::array kHEst = { "eGlobaltrack", "eFV0", "e1flatencityFV0", "eFT0", "e1flatencityFT0", "eFV0FT0C", "e1flatencityFV0FT0C", "ePtTrig"}; - static constexpr std::string_view kTEst[kNEstimators] = { + static constexpr std::array kTEst = { "GlobalTrk", "FV0", "1-flatencity_FV0", "FT0", "1-flatencityFT0", "FV0_FT0C", "1-flatencity_FV0_FT0C", "PtTrig"}; - static constexpr std::string_view kHPtEst[kNEstimators] = { + static constexpr std::array kHPtEst = { "ptVsGlobaltrack", "ptVsFV0", "ptVs1flatencityFV0", "ptVsFT0", "ptVs1flatencityFT0", "ptVsFV0FT0C", "ptVs1flatencityFV0FT0C", "pTVsPtTrig"}; - // Configurable for histograms - Configurable nBinsVz{"nBinsVz", 100, "N bins in Vz"}; - Configurable nBinsK0sMass{"nBinsK0sMass", 400, "N bins in K0sMass"}; - Configurable nBinsLambdaMass{"nBinsLambdaMass", 400, - "N bins in LambdaMass"}; - Configurable nBinsXiMass{"nBinsXiMass", 400, "N bins in XiMass"}; - - Configurable kK0sEPshiftfromMass{"kK0sEPshiftfromMass", 0.1, "distance of K0s Inv mass histogram start and end points from PDG mass"}; - Configurable kLambdaEPshiftfromMass{"kLambdaEPshiftfromMass", 0.05, "distance of Lambda Inv mass histogram start and end points from PDG mass"}; - Configurable kXiEPshiftfromMass{"kXiEPshiftfromMass", 0.05, "distance of Xi Inv mass histogram start and end points from PDG mass"}; - - Configurable nBinspT{"nBinspT", 250, "N bins in pT"}; - Configurable nBinsFlattenicity{"nBinsFlattenicity", 100, "N bins in Flattenicity"}; - - // Configurable for event selection - - Configurable applyEvSel{"applyEvSel", true, - "Apply event selection to Data and MCRec"}; - Configurable issel8{"issel8", true, - "Accept events that pass sel8 selection"}; - Configurable cutzvertex{"cutzvertex", 10.0f, - "Accepted z-vertex range (cm)"}; - Configurable isINELgt0{"isINELgt0", true, "is INEL gt 0"}; - Configurable isNoTimeFrameBorder{"isNoTimeFrameBorder", true, - "cut branch crossing at the beginning/end of TF"}; - Configurable isNoITSROFrameBorder{"isNoITSROFrameBorder", true, - "cut branch crossing at the beginning/end of ITS ROF"}; - Configurable isVertexITSTPC{"isVertexITSTPC", false, - "Is Vertex ITSTPC"}; - Configurable isNoSameBunchPileup{"isNoSameBunchPileup", false, - "Is No Same Bunch Pileup"}; - Configurable isGoodZvtxFT0vsPV{"isGoodZvtxFT0vsPV", false, - "Is Good Zvtx FT0 vs PV"}; - Configurable isTriggerTVX{"isTriggerTVX", true, - "coincidence of a signal in FT0A and FT0C"}; - - // Configurables for Flattenicity - Configurable flattenicityQA{"flattenicityQA", true, "Store Flattenicity QA plots"}; - Configurable applyCalibCh{"applyCalibCh", false, "equalize FV0"}; - Configurable applyCalibVtx{"applyCalibVtx", false, - "equalize FV0 vs vtx"}; - Configurable applyNorm{"applyNorm", false, "normalization to eta"}; - Configurable isflattenicitywithFV0{"isflattenicitywithFV0", true, - "Calculate Flattenicity with FV0"}; - Configurable isflattenicitywithFT0{"isflattenicitywithFT0", true, - "Calculate Flattenicity with FT0"}; - Configurable isflattenicitywithFV0FT0C{"isflattenicitywithFV0FT0C", true, - "Calculate Flattenicity with FV0+FT0C"}; - - // accepted values of the flattenicityforanalysis configurable + // Histogram binning + struct : ConfigurableGroup { + std::string prefix = "binning"; + Configurable nBinsVz{"nBinsVz", 100, "N bins in Vz"}; + Configurable nBinsK0sMass{"nBinsK0sMass", 400, "N bins in K0sMass"}; + Configurable nBinsLambdaMass{"nBinsLambdaMass", 400, + "N bins in LambdaMass"}; + Configurable nBinsXiMass{"nBinsXiMass", 400, "N bins in XiMass"}; + Configurable kK0sEPshiftfromMass{"kK0sEPshiftfromMass", 0.1, "distance of K0s Inv mass histogram start and end points from PDG mass"}; + Configurable kLambdaEPshiftfromMass{"kLambdaEPshiftfromMass", 0.05, "distance of Lambda Inv mass histogram start and end points from PDG mass"}; + Configurable kXiEPshiftfromMass{"kXiEPshiftfromMass", 0.05, "distance of Xi Inv mass histogram start and end points from PDG mass"}; + Configurable nBinspT{"nBinspT", 250, "N bins in pT"}; + Configurable nBinsFlattenicity{"nBinsFlattenicity", 100, "N bins in Flattenicity"}; + } binning; + + // Event selection + struct : ConfigurableGroup { + std::string prefix = "evSel"; + Configurable applyEvSel{"applyEvSel", true, + "Apply event selection to Data and MCRec"}; + Configurable issel8{"issel8", true, + "Accept events that pass sel8 selection"}; + Configurable cutzvertex{"cutzvertex", 10.0f, + "Accepted z-vertex range (cm)"}; + Configurable isINELgt0{"isINELgt0", true, "is INEL gt 0"}; + Configurable isNoTimeFrameBorder{"isNoTimeFrameBorder", true, + "cut branch crossing at the beginning/end of TF"}; + Configurable isNoITSROFrameBorder{"isNoITSROFrameBorder", true, + "cut branch crossing at the beginning/end of ITS ROF"}; + Configurable isVertexITSTPC{"isVertexITSTPC", false, + "Is Vertex ITSTPC"}; + Configurable isNoSameBunchPileup{"isNoSameBunchPileup", false, + "Is No Same Bunch Pileup"}; + Configurable isGoodZvtxFT0vsPV{"isGoodZvtxFT0vsPV", false, + "Is Good Zvtx FT0 vs PV"}; + Configurable isTriggerTVX{"isTriggerTVX", true, + "coincidence of a signal in FT0A and FT0C"}; + } evSel; + + // Flattenicity estimation + struct : ConfigurableGroup { + std::string prefix = "flatSel"; + Configurable flattenicityQA{"flattenicityQA", true, "Store Flattenicity QA plots"}; + Configurable applyCalibCh{"applyCalibCh", false, "equalize FV0"}; + Configurable applyCalibVtx{"applyCalibVtx", false, + "equalize FV0 vs vtx"}; + Configurable applyNorm{"applyNorm", false, "normalization to eta"}; + Configurable isflattenicitywithFV0{"isflattenicitywithFV0", true, + "Calculate Flattenicity with FV0"}; + Configurable isflattenicitywithFT0{"isflattenicitywithFT0", true, + "Calculate Flattenicity with FT0"}; + Configurable isflattenicitywithFV0FT0C{"isflattenicitywithFV0FT0C", true, + "Calculate Flattenicity with FV0+FT0C"}; + Configurable flattenicityforanalysis{"flattenicityforanalysis", 0, + "Which Flattenicity to be used for analysis, 0 for FV0, 1 for FT0, 2 for FV0+FT0C"}; + Configurable flattenicityforLossCorrRec{"flattenicityforLossCorrRec", true, + "Flattenicity from Rec Tracks are used for Signal and Event loss calculations"}; + } flatSel; + + // V0 selection + struct : ConfigurableGroup { + std::string prefix = "v0Sel"; + Configurable v0settingDCAv0dau{"v0settingDCAv0dau", 1, + "DCA V0 Daughters"}; + Configurable v0settingDCApostopv{"v0settingDCApostopv", 0.06, + "DCA Pos To PV"}; + Configurable v0settingDCAnegtopv{"v0settingDCAnegtopv", 0.06, + "DCA Neg To PV"}; + Configurable v0settingDCAbactopv{"v0settingDCAbactopv", 0.06, + "DCA Bchelor To PV"}; + Configurable v0settingRapidity{"v0settingRapidity", 0.5, + "V0 rapidity cut"}; + Configurable v0settingCosPAK0s{"v0settingCosPAK0s", 0.97, + "V0 CosPA for K0s"}; + Configurable v0settingRadiusK0s{"v0settingRadiusK0s", 0.5, + "v0radius for K0s"}; + Configurable v0settingcTauK0s{"v0settingcTauK0s", 20, + "v0ctau for K0s"}; + Configurable v0settingMassRejectionK0s{"v0settingMassRejectionK0s", 0.005, + "Competing Mass Rejection cut for K0s"}; + Configurable v0settingArmePodoK0s{"v0settingArmePodoK0s", 0.2, + "Armenteros-Podolanski cut for K0s"}; + Configurable v0settingCosPALambda{"v0settingCosPALambda", 0.995, + "V0 CosPA for Lambda"}; + Configurable v0settingRadiusLambda{"v0settingRadiusLambda", 0.5, + "v0radius for Lambda"}; + Configurable v0settingcTauLambda{"v0settingcTauLambda", 30, + "v0ctau for Lambda"}; + Configurable v0settingMassRejectionLambda{"v0settingMassRejectionLambda", 0.01, + "Competing Mass Rejection cut for Lambda"}; + Configurable v0settingNTPCcrossedRows{"v0settingNTPCcrossedRows", 70, + "Minimum number of TPC crossed pad rows of the V0 daughters, negative: no cut"}; + } v0Sel; + + // Daughter track acceptance and PID + struct : ConfigurableGroup { + std::string prefix = "trkPid"; + Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, + "Eta range for tracks"}; + Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.0f, "Minimum pT"}; + Configurable nSigmaTPCPion{"nSigmaTPCPion", 5, "nSigmaTPCPion"}; + Configurable nSigmaTPCProton{"nSigmaTPCProton", 5, "nSigmaTPCProton"}; + Configurable pidQAWindowK0s{"pidQAWindowK0s", 0.05, + "Half width of the K0s mass window used for the PID QA plots"}; + Configurable pidQAWindowLambda{"pidQAWindowLambda", 0.1, + "Half width of the Lambda mass window used for the PID QA plots"}; + } trkPid; + + // Cascade selection + struct : ConfigurableGroup { + std::string prefix = "cascSel"; + Configurable nTPCcrossedRows{"nTPCcrossedRows", 52, "Number of TPC crossed pad raws"}; + Configurable cascsettingDCAv0toPV{"cascsettingDCAv0toPV", 0.03, "DCA V0 To PV"}; + Configurable cascsettingDCAv0bach{"cascsettingDCAv0bach", 0.25, "DCA V0 To bachelor"}; + Configurable cascsettingDCAbaryontopv{"cascsettingDCAbaryontopv", 0.05, "DCA of the baryon daughter To PV"}; + Configurable cascsettingDCAmesontopv{"cascsettingDCAmesontopv", 0.1, "DCA of the meson daughter To PV"}; + Configurable cascsettingDCAxybaryonbach{"cascsettingDCAxybaryonbach", 0.02, "DCAxy Bachelor-Baryon To PV, below which the candidate is vetoed"}; + Configurable cascsettingCosPAbaryonbach{"cascsettingCosPAbaryonbach", 0.9999, "CosThetap Bachelor-Baryon, above which the candidate is vetoed"}; + Configurable cascsettingCosPAcascPV{"cascsettingCosPAcascPV", 0.9947, "CosThetap for Cascade to PV"}; + Configurable cascsettingCosPAv0PV{"cascsettingCosPAv0PV", 0.9876, "CosThetap for V0 to PV"}; + Configurable cascsettingv0radius{"cascsettingv0radius", 0.55, "V0 decay radius for cadcades in cm"}; + Configurable cascsettingcascradius{"cascsettingcascradius", 1.01, "Cascade decay radius for cadcades in cm"}; + Configurable cascsettingRapidity{"cascsettingRapidity", 0.5, "Cascade rapidity cut"}; + Configurable cascsettingMassRejectionLambdaXi{"cascsettingMassRejectionLambdaXi", 0.0116, "Casc Mass Rejection cut of Lambda for Xi"}; + Configurable cascsettingMassRejectioOmegaXi{"cascsettingMassRejectioOmegaXi", -1, "Casc Mass Rejection cut of Omega for Xi"}; + Configurable cascsettingproplifetime{"cascsettingproplifetime", 4.6, "Scale for lifetime cut on ctau Xi"}; + } cascSel; + + // values for flattenicityforanalysis static constexpr int kFlatFromFV0 = 0; static constexpr int kFlatFromFT0 = 1; static constexpr int kFlatFromFV0FT0C = 2; - Configurable flattenicityforanalysis{"flattenicityforanalysis", 0, - "Which Flattenicity to be used for analysis, 0 for FV0, 1 for FT0, 2 for FV0+FT0C"}; - Configurable flattenicityforLossCorrRec{"flattenicityforLossCorrRec", true, - "Flattenicity from Rec Tracks are used for Signal and Event loss calculations"}; - // Common Configurable parameters for V0 selection - Configurable v0settingDCAv0dau{"v0settingDCAv0dau", 1, - "DCA V0 Daughters"}; - Configurable v0settingDCApostopv{"v0settingDCApostopv", 0.06, - "DCA Pos To PV"}; - Configurable v0settingDCAnegtopv{"v0settingDCAnegtopv", 0.06, - "DCA Neg To PV"}; - Configurable v0settingDCAbactopv{"v0settingDCAbactopv", 0.06, - "DCA Bchelor To PV"}; - Configurable v0settingRapidity{"v0settingRapidity", 0.5, - "V0 rapidity cut"}; - - // Configurable parameters for V0 selection for KOs - Configurable v0settingCosPAK0s{"v0settingCosPAK0s", 0.97, - "V0 CosPA for K0s"}; - Configurable v0settingRadiusK0s{"v0settingRadiusK0s", 0.5, - "v0radius for K0s"}; - Configurable v0settingcTauK0s{"v0settingcTauK0s", 20, - "v0ctau for K0s"}; - Configurable v0settingMassRejectionK0s{"v0settingMassRejectionK0s", 0.005, - "Competing Mass Rejection cut for K0s"}; - Configurable v0settingArmePodoK0s{"v0settingArmePodoK0s", 0.2, - "Armenteros-Podolanski cut for K0s"}; - - // Configurable parameters for V0 selection for Lambda - Configurable v0settingCosPALambda{"v0settingCosPALambda", 0.995, - "V0 CosPA for Lambda"}; - Configurable v0settingRadiusLambda{"v0settingRadiusLambda", 0.5, - "v0radius for Lambda"}; - Configurable v0settingcTauLambda{"v0settingcTauLambda", 30, - "v0ctau for Lambda"}; - Configurable v0settingMassRejectionLambda{"v0settingMassRejectionLambda", 0.01, - "Competing Mass Rejection cut for Lambda"}; - - // Track quality of the V0 daughters - Configurable v0settingNTPCcrossedRows{"v0settingNTPCcrossedRows", 70, - "Minimum number of TPC crossed pad rows of the V0 daughters, negative: no cut"}; - - // Configurable parameters for PID selection - Configurable nSigmaTPCPion{"nSigmaTPCPion", 5, "nSigmaTPCPion"}; - Configurable nSigmaTPCProton{"nSigmaTPCProton", 5, "nSigmaTPCProton"}; - - // Mass windows around the PDG mass in which the daughter PID QA plots are filled - Configurable pidQAWindowK0s{"pidQAWindowK0s", 0.05, - "Half width of the K0s mass window used for the PID QA plots"}; - Configurable pidQAWindowLambda{"pidQAWindowLambda", 0.1, - "Half width of the Lambda mass window used for the PID QA plots"}; + // FT0M percentile class: run once over the full range for MB and once with a + // narrow window for HM. Grouped because the framework only decomposes 100 + // struct members (Framework/StructToTuple.h) + struct : ConfigurableGroup { + std::string prefix = "eventClass"; + Configurable applyCentSel{"applyCentSel", false, + "Select events in a FT0M percentile window"}; + Configurable cfgCentMin{"cfgCentMin", 0.0f, "Minimum FT0M percentile"}; + Configurable cfgCentMax{"cfgCentMax", 100.0f, "Maximum FT0M percentile"}; + Configurable nBinsCent{"nBinsCent", 100, "N bins in FT0M percentile"}; + + // keep only primaries in the MC-matched spectra, the rest go to the + // feed-down histograms + Configurable requirePrimaryMC{"requirePrimaryMC", true, + "Require isPhysicalPrimary() on the MC-matched candidate"}; + } eventClass; // Configurable v0daughter_etacut{"V0DaughterEtaCut", 0.8, // "V0DaughterEtaCut"}; // Configurable v0etacut{"v0etacut", 0.8, "v0etacut"}; - // acceptance cuts for Flattenicity correlation - Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, - "Eta range for tracks"}; - Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.0f, "Minimum pT"}; - - // Additional Cut configurables for Cascades - Configurable nTPCcrossedRows{"nTPCcrossedRows", 52, "Number of TPC crossed pad raws"}; - Configurable cascsettingDCAv0toPV{"cascsettingDCAv0toPV", 0.03, "DCA V0 To PV"}; - Configurable cascsettingDCAv0bach{"cascsettingDCAv0bach", 0.25, "DCA V0 To bachelor"}; - Configurable cascsettingDCAbaryontopv{"cascsettingDCAbaryontopv", 0.05, "DCA of the baryon daughter To PV"}; - Configurable cascsettingDCAmesontopv{"cascsettingDCAmesontopv", 0.1, "DCA of the meson daughter To PV"}; - Configurable cascsettingDCAxybaryonbach{"cascsettingDCAxybaryonbach", 0.02, "DCAxy Bachelor-Baryon To PV, below which the candidate is vetoed"}; - Configurable cascsettingCosPAbaryonbach{"cascsettingCosPAbaryonbach", 0.9999, "CosThetap Bachelor-Baryon, above which the candidate is vetoed"}; - Configurable cascsettingCosPAcascPV{"cascsettingCosPAcascPV", 0.9947, "CosThetap for Cascade to PV"}; - Configurable cascsettingCosPAv0PV{"cascsettingCosPAv0PV", 0.9876, "CosThetap for V0 to PV"}; - Configurable cascsettingv0radius{"cascsettingv0radius", 0.55, "V0 decay radius for cadcades in cm"}; - Configurable cascsettingcascradius{"cascsettingcascradius", 1.01, "Cascade decay radius for cadcades in cm"}; - Configurable cascsettingRapidity{"cascsettingRapidity", 0.5, "Cascade rapidity cut"}; - Configurable cascsettingMassRejectionLambdaXi{"cascsettingMassRejectionLambdaXi", 0.0116, "Casc Mass Rejection cut of Lambda for Xi"}; - Configurable cascsettingMassRejectioOmegaXi{"cascsettingMassRejectioOmegaXi", -1, "Casc Mass Rejection cut of Omega for Xi"}; - Configurable cascsettingproplifetime{"cascsettingproplifetime", 4.6, "Scale for lifetime cut on ctau Xi"}; - int nbin = 1; - // bin of hEventsSelected that counts the events surviving the flattenicity - // requirement; -1 when no process function applies that requirement + // hEventsSelected bin for the flattenicity requirement, -1 if not applied int nbinFlattenicity = -1; - // processGenMC also evaluates the reconstructed flattenicity; it must not fill the - // detector QA a second time when a reconstructed-level process function is running + // do not fill the detector QA twice when processGenMC runs with a rec-level process bool fillFlattenicityQAInGenMC = true; - // Vertex-dependent equalization curves for FV0, FT0A and FT0C. - // Built once in init() to avoid re-allocating a TGraph for every collision. + // vertex equalization curves, built once in init() static constexpr int kNDetVtx = 3; static constexpr int kNVtxPoints = 30; std::array gVtx; @@ -275,30 +295,29 @@ struct Lambdak0sflattenicity { void init(InitContext const&) { // Axes - // the configurables give the half width around the PDG mass - AxisSpec k0sMassAxis = {nBinsK0sMass, o2::constants::physics::MassK0Short - kK0sEPshiftfromMass, - o2::constants::physics::MassK0Short + kK0sEPshiftfromMass, + AxisSpec k0sMassAxis = {binning.nBinsK0sMass, o2::constants::physics::MassK0Short - binning.kK0sEPshiftfromMass, + o2::constants::physics::MassK0Short + binning.kK0sEPshiftfromMass, "#it{M}_{#pi^{+}#pi^{-}} [GeV/#it{c}^{2}]"}; - AxisSpec lambdaMassAxis = {nBinsLambdaMass, o2::constants::physics::MassLambda0 - kLambdaEPshiftfromMass, - o2::constants::physics::MassLambda0 + kLambdaEPshiftfromMass, + AxisSpec lambdaMassAxis = {binning.nBinsLambdaMass, o2::constants::physics::MassLambda0 - binning.kLambdaEPshiftfromMass, + o2::constants::physics::MassLambda0 + binning.kLambdaEPshiftfromMass, "#it{M}_{p#pi^{-}} [GeV/#it{c}^{2}]"}; - AxisSpec antilambdaMassAxis = {nBinsLambdaMass, o2::constants::physics::MassLambda0 - kLambdaEPshiftfromMass, - o2::constants::physics::MassLambda0 + kLambdaEPshiftfromMass, + AxisSpec antilambdaMassAxis = {binning.nBinsLambdaMass, o2::constants::physics::MassLambda0 - binning.kLambdaEPshiftfromMass, + o2::constants::physics::MassLambda0 + binning.kLambdaEPshiftfromMass, "#it{M}_{#pi^{+}#bar{p}} [GeV/#it{c}^{2}]"}; - AxisSpec xiMassAxis = {nBinsXiMass, o2::constants::physics::MassXiMinus - kXiEPshiftfromMass, - o2::constants::physics::MassXiMinus + kXiEPshiftfromMass, + AxisSpec xiMassAxis = {binning.nBinsXiMass, o2::constants::physics::MassXiMinus - binning.kXiEPshiftfromMass, + o2::constants::physics::MassXiMinus + binning.kXiEPshiftfromMass, "#it{M}_{#Lambda#pi} [GeV/#it{c}^{2}]"}; - AxisSpec vertexZAxis = {nBinsVz, -15., 15., "vrtx_{Z} [cm]"}; - AxisSpec ptAxis = {nBinspT, 0.0f, 25.0f, "#it{p}_{T} (GeV/#it{c})"}; - AxisSpec pTPCAxis = {nBinspT, 0.0f, 25.0f, "#it{p}_{TPC} (GeV/#it{c})"}; + AxisSpec vertexZAxis = {binning.nBinsVz, -15., 15., "vrtx_{Z} [cm]"}; + AxisSpec ptAxis = {binning.nBinspT, 0.0f, 25.0f, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec pTPCAxis = {binning.nBinspT, 0.0f, 25.0f, "#it{p}_{TPC} (GeV/#it{c})"}; AxisSpec decayRadiusAxis = {100, 0.0f, 100.0f, "Decay Radius (cm)"}; - AxisSpec flatAxis = {nBinsFlattenicity, 0.0f, 1.0f, "1-#rho_{ch}"}; + AxisSpec flatAxis = {binning.nBinsFlattenicity, 0.0f, 1.0f, "1-#rho_{ch}"}; + AxisSpec centAxis = {eventClass.nBinsCent, 0.0f, 100.0f, "FT0M percentile"}; - int nBinsEst[8] = {100, 500, 102, 500, 102, 500, 102, 150}; - float lowEdgeEst[8] = {-0.5, -0.5, -0.01, -0.5, -0.01, -0.5, -0.01, .0}; - float upEdgeEst[8] = {99.5, 49999.5, 1.01, 499.5, 1.01, 499.5, 1.01, 150.0}; + std::array nBinsEst = {100, 500, 102, 500, 102, 500, 102, 150}; + std::array lowEdgeEst = {-0.5, -0.5, -0.01, -0.5, -0.01, -0.5, -0.01, .0}; + std::array upEdgeEst = {99.5, 49999.5, 1.01, 499.5, 1.01, 499.5, 1.01, 150.0}; - // Vertex equalization curves, evaluated per collision in estimateFlattenicity() gVtx[0].SetName("gAmpV0"); gVtx[1].SetName("gAmpT0A"); gVtx[2].SetName("gAmpT0C"); @@ -316,35 +335,34 @@ struct Lambdak0sflattenicity { {HistType::kTH1D, {{15, 0, 15}}}); rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "all"); - if (issel8) { + if (evSel.issel8) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "sel8"); } rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "zvertex"); - if (isNoTimeFrameBorder) { + if (evSel.isNoTimeFrameBorder) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "TFBorder"); } - if (isNoITSROFrameBorder) { + if (evSel.isNoITSROFrameBorder) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "ITSROFBorder"); } - if (isVertexITSTPC) { + if (evSel.isVertexITSTPC) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "VertexITSTPC"); } - if (isNoSameBunchPileup) { + if (evSel.isNoSameBunchPileup) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "SameBunchPileup"); } - if (isGoodZvtxFT0vsPV) { + if (evSel.isGoodZvtxFT0vsPV) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "isGoodZvtxFT0vsPV"); } - if (isTriggerTVX) { + if (evSel.isTriggerTVX) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "TVX"); } - if (isINELgt0) { + if (evSel.isINELgt0) { rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "INEL>0"); } - // the process functions below reject events for which the forward detectors - // carry no usable information, so that rejection needs its own counter + // events without FV0/FT0 are rejected below, they need their own counter if (doprocessDataRun3LambdaK0s || doprocessRecMCLambdaK0s || doprocessDataRun3Cascade || doprocessRecMCRun3Cascade) { nbinFlattenicity = nbin; @@ -353,6 +371,9 @@ struct Lambdak0sflattenicity { rEventSelection.add("hFlattenicityDistribution", "hFlattenicityDistribution", {HistType::kTH1D, {flatAxis}}); + rEventSelection.add("hCentFT0M", "hCentFT0M", {HistType::kTH1D, {centAxis}}); + rEventSelection.add("hCentFT0MvsFlattenicity", "hCentFT0MvsFlattenicity", + {HistType::kTH2D, {centAxis, flatAxis}}); if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessGenMC) { rEventSelection.add("hTrueFV0amplvsFlat", "TrueFV0MvsFlat", HistType::kTH2D, {{500, -0.5, +499.5, "True Nch in FV0 region"}, flatAxis}); @@ -430,6 +451,16 @@ struct Lambdak0sflattenicity { {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); rLambda.add("Generated_MCRecoCollCheck_INELgt0_Lambda", "Generated_MCRecoCollCheck_INELgt0_Lambda", {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + rLambda.add("hMassFeedDownLambdapTFlat", "hMassFeedDownLambdapTFlat", + {HistType::kTH3D, {{lambdaMassAxis}, {ptAxis}, {flatAxis}}}); + rLambda.add("hFeedDownLambdaPtVsMotherPt", "hFeedDownLambdaPtVsMotherPt", + {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + rLambda.add("hFeedDownLambdaMotherPdg", "hFeedDownLambdaMotherPdg", + {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + rLambda.get(HIST("hFeedDownLambdaMotherPdg"))->GetXaxis()->SetBinLabel(1, "#Xi^{-}"); + rLambda.get(HIST("hFeedDownLambdaMotherPdg"))->GetXaxis()->SetBinLabel(2, "#Xi^{0}"); + rLambda.get(HIST("hFeedDownLambdaMotherPdg"))->GetXaxis()->SetBinLabel(3, "#Omega^{-}"); + rLambda.get(HIST("hFeedDownLambdaMotherPdg"))->GetXaxis()->SetBinLabel(4, "other"); } // AntiLambda reconstruction @@ -465,6 +496,16 @@ struct Lambdak0sflattenicity { {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); rAntiLambda.add("Generated_MCRecoCollCheck_INELgt0_AntiLambda", "Generated_MCRecoCollCheck_INELgt0_AntiLambda", {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + rAntiLambda.add("hMassFeedDownAntiLambdapTFlat", "hMassFeedDownAntiLambdapTFlat", + {HistType::kTH3D, {{antilambdaMassAxis}, {ptAxis}, {flatAxis}}}); + rAntiLambda.add("hFeedDownAntiLambdaPtVsMotherPt", "hFeedDownAntiLambdaPtVsMotherPt", + {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + rAntiLambda.add("hFeedDownAntiLambdaMotherPdg", "hFeedDownAntiLambdaMotherPdg", + {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + rAntiLambda.get(HIST("hFeedDownAntiLambdaMotherPdg"))->GetXaxis()->SetBinLabel(1, "#Xi^{-}"); + rAntiLambda.get(HIST("hFeedDownAntiLambdaMotherPdg"))->GetXaxis()->SetBinLabel(2, "#Xi^{0}"); + rAntiLambda.get(HIST("hFeedDownAntiLambdaMotherPdg"))->GetXaxis()->SetBinLabel(3, "#Omega^{-}"); + rAntiLambda.get(HIST("hFeedDownAntiLambdaMotherPdg"))->GetXaxis()->SetBinLabel(4, "other"); } rCommonHist.add("hArmPodoAlphavsQT", "hArmPodoAlphavsQT", @@ -502,6 +543,16 @@ struct Lambdak0sflattenicity { {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); rXi.add("Generated_MCRecoCollCheck_INELgt0_Xi", "Generated_MCRecoCollCheck_INELgt0_Xi", {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + rXi.add("hMassFeedDownXipTFlat", "hMassFeedDownXipTFlat", + {HistType::kTH3D, {{xiMassAxis}, {ptAxis}, {flatAxis}}}); + rXi.add("hFeedDownXiPtVsMotherPt", "hFeedDownXiPtVsMotherPt", + {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + rXi.add("hFeedDownXiMotherPdg", "hFeedDownXiMotherPdg", + {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + rXi.get(HIST("hFeedDownXiMotherPdg"))->GetXaxis()->SetBinLabel(1, "#Xi^{-}"); + rXi.get(HIST("hFeedDownXiMotherPdg"))->GetXaxis()->SetBinLabel(2, "#Xi^{0}"); + rXi.get(HIST("hFeedDownXiMotherPdg"))->GetXaxis()->SetBinLabel(3, "#Omega^{-}"); + rXi.get(HIST("hFeedDownXiMotherPdg"))->GetXaxis()->SetBinLabel(4, "other"); } } if (doprocessGenMC) { @@ -587,7 +638,7 @@ struct Lambdak0sflattenicity { {HistType::kTH2D, {ptAxis, flatAxis}}); } - if (flattenicityQA) { + if (flatSel.flattenicityQA) { rFlattenicity.add("hEv", "Ev", HistType::kTH1D, {{6, -0.5, 5.5, "index activated detector"}}); rFlattenicity.add("hFV0amplRing1to4", "FV01to4", HistType::kTH1D, @@ -679,15 +730,14 @@ struct Lambdak0sflattenicity { fillFlattenicityQAInGenMC = !(doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade); - // the estimator selected for the analysis has to be actually computed, - // otherwise it keeps its sentinel value and every event is rejected - if (flattenicityforanalysis == kFlatFromFV0 && !isflattenicitywithFV0 && !isflattenicitywithFV0FT0C) { + // the estimator used for the analysis has to be computed + if (flatSel.flattenicityforanalysis == kFlatFromFV0 && !flatSel.isflattenicitywithFV0 && !flatSel.isflattenicitywithFV0FT0C) { LOGF(fatal, "flattenicityforanalysis=0 (FV0) needs isflattenicitywithFV0 or isflattenicitywithFV0FT0C enabled"); } - if (flattenicityforanalysis == kFlatFromFT0 && !isflattenicitywithFT0) { + if (flatSel.flattenicityforanalysis == kFlatFromFT0 && !flatSel.isflattenicitywithFT0) { LOGF(fatal, "flattenicityforanalysis=1 (FT0) needs isflattenicitywithFT0 enabled"); } - if (flattenicityforanalysis == kFlatFromFV0FT0C && !(isflattenicitywithFV0FT0C || (isflattenicitywithFV0 && isflattenicitywithFT0))) { + if (flatSel.flattenicityforanalysis == kFlatFromFV0FT0C && !(flatSel.isflattenicitywithFV0FT0C || (flatSel.isflattenicitywithFV0 && flatSel.isflattenicitywithFT0))) { LOGF(fatal, "flattenicityforanalysis=2 (FV0+FT0C) needs isflattenicitywithFV0FT0C enabled"); } } @@ -705,8 +755,7 @@ struct Lambdak0sflattenicity { return (iCh >= 0 && iSec < kNSectorsT0C) ? iSec : -1; } - // the four inner FV0 rings hold kNChannelsPerFV0Ring channels each; everything - // beyond them belongs to the outermost ring + // four inner rings of kNChannelsPerFV0Ring channels each, the rest is the outer ring int getFV0Ring(int iCh) { const int iRing = iCh / kNChannelsPerFV0Ring; @@ -716,9 +765,7 @@ struct Lambdak0sflattenicity { return (iRing < kOuterFV0RingIndex) ? iRing : kOuterFV0RingIndex; } - // FV0 channel -> phi-ordered index; a fixed permutation of the 48 channels - // (within each inner ring the last four run backwards, and the two outer - // half-rings interleave), tabulated in kFV0PhiIndex rather than branched on. + // FV0 channel -> phi-ordered index, a fixed permutation of the 48 channels int getFV0IndexPhi(int iCh) { if (iCh < 0 || iCh >= kNCells) { @@ -750,7 +797,7 @@ struct Lambdak0sflattenicity { return flat; } // V0A signal and flatenicity calculation - static constexpr float kCalib[48] = { + static constexpr std::array kCalib = { 1.01697, 1.122, 1.03854, 1.108, 1.11634, 1.14971, 1.19321, 1.06866, 0.954675, 0.952695, 0.969853, 0.957557, 0.989784, 1.01549, 1.02182, 0.976005, 1.01865, 1.06871, 1.06264, 1.02969, 1.07378, @@ -759,39 +806,39 @@ struct Lambdak0sflattenicity { 1.00122, 1.03303, 0.887866, 0.892437, 0.906278, 0.884976, 0.864251, 0.917221, 1.10618, 1.04028, 0.893184, 0.915734, 0.892676}; // calibration T0C - static constexpr float kCalibT0C[28] = { + static constexpr std::array kCalibT0C = { 0.949829, 1.05408, 1.00681, 1.00724, 0.990663, 0.973571, 0.9855, 1.03726, 1.02526, 1.00467, 0.983008, 0.979349, 0.952352, 0.985775, 1.013, 1.01721, 0.993948, 0.996421, 0.971871, 1.02921, 0.989641, 1.01885, 1.01259, 0.929502, 1.03969, 1.02496, 1.01385, 1.01711}; // calibration T0A - static constexpr float kCalibT0A[24] = { + static constexpr std::array kCalibT0A = { 0.86041, 1.10607, 1.17724, 0.756397, 1.14954, 1.0879, 0.829438, 1.09014, 1.16515, 0.730077, 1.06722, 0.906344, 0.824167, 1.14716, 1.20692, 0.755034, 1.11734, 1.00556, 0.790522, 1.09138, 1.16225, 0.692458, 1.12428, 1.01127}; // calibration factor MFT vs vtx - static constexpr float kBiningVtxt[30] = { + static constexpr std::array kBiningVtxt = { -14.5, -13.5, -12.5, -11.5, -10.5, -9.5, -8.5, -7.5, -6.5, -5.5, -4.5, -3.5, -2.5, -1.5, -0.5, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5, 8.5, 9.5, 10.5, 11.5, 12.5, 13.5, 14.5}; // calibration factor FV0 vs vtx - static constexpr float kCalibFV0vtx[30] = { + static constexpr std::array kCalibFV0vtx = { 0.907962, 0.934607, 0.938929, 0.950987, 0.950817, 0.966362, // o2-linter: disable=pdg/explicit-mass (these are not masses but calibration values) 0.968509, 0.972741, 0.982412, 0.984872, 0.994543, 0.996003, 0.99435, 1.00266, 0.998245, 1.00584, 1.01078, 1.01003, 1.00726, 1.00872, 1.01726, 1.02015, 1.0193, 1.01106, 1.02229, 1.02104, 1.03435, 1.00822, 1.01921, 1.01736}; // calibration FT0A vs vtx - static constexpr float kCalibFT0Avtx[30] = { + static constexpr std::array kCalibFT0Avtx = { 0.924334, 0.950988, 0.959604, 0.965607, 0.970016, 0.979057, 0.978384, 0.982005, 0.992825, 0.990048, 0.998588, 0.997338, 1.00102, 1.00385, 0.99492, 1.01083, 1.00703, 1.00494, 1.00063, 1.0013, 1.00777, 1.01238, 1.01179, 1.00577, 1.01028, 1.017, 1.02975, 1.0085, 1.00856, 1.01662}; // calibration FT0C vs vtx - static constexpr float kCalibFT0Cvtx[30] = { + static constexpr std::array kCalibFT0Cvtx = { 1.02096, 1.01245, 1.02148, 1.03605, 1.03561, 1.03667, 1.04229, 1.0327, 1.03674, 1.02764, 1.01828, 1.02331, 1.01864, 1.015, 1.01197, 1.00615, 0.996845, 0.993051, @@ -799,23 +846,23 @@ struct Lambdak0sflattenicity { 0.956687, 0.932816, 0.92773, 0.914892, 0.891724, 0.872382}; static constexpr int kNeta5 = 2; // FT0C + FT0A - static constexpr float kWeigthsEta5[kNeta5] = {0.0490638, 0.010958415}; - static constexpr float kDeltaEeta5[kNeta5] = {1.1, 1.2}; + static constexpr std::array kWeigthsEta5 = {0.0490638, 0.010958415}; + static constexpr std::array kDeltaEeta5 = {1.1, 1.2}; static constexpr int kNeta6 = 2; // FT0C + FV0 - static constexpr float kWeigthsEta6[kNeta6] = {0.0490638, 0.00353962}; - static constexpr float kDeltaEeta6[kNeta6] = {1.1, 2.9}; + static constexpr std::array kWeigthsEta6 = {0.0490638, 0.00353962}; + static constexpr std::array kDeltaEeta6 = {1.1, 2.9}; static constexpr int kInnerFV0 = 32; static constexpr float kMaxEtaFV0 = 5.1; static constexpr float kMinEtaFV0 = 2.2; static constexpr float kDetaFV0 = (kMaxEtaFV0 - kMinEtaFV0) / 5.0; - // returned when the forward detectors carry no usable information + // no FV0/FT0 information static constexpr float kInvalidFlattenicity = -1.f; static constexpr int kNCells = 48; // 48 sectors in FV0 - static constexpr int kFV0PhiIndex[kNCells] = { + static constexpr std::array kFV0PhiIndex = { 0, 1, 2, 3, 7, 6, 5, 4, 8, 9, 10, 11, 15, 14, 13, 12, 16, 17, 18, 19, 23, 22, 21, 20, @@ -833,8 +880,7 @@ struct Lambdak0sflattenicity { std::array estimator{}; - // fillCounter=false evaluates the very same cuts without filling hEventsSelected, - // for the callers that only need to know whether a collision would be accepted + // fillCounter=false runs the same cuts without filling hEventsSelected template bool isEventSelected(TCollision const& collision, bool fillCounter = true) { @@ -843,17 +889,17 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } - if (issel8 && !collision.sel8()) { + if (evSel.issel8 && !collision.sel8()) { return false; } - if (issel8) { + if (evSel.issel8) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (std::abs(collision.posZ()) > cutzvertex) { + if (std::abs(collision.posZ()) > evSel.cutzvertex) { return false; } @@ -862,74 +908,78 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } - if (isNoTimeFrameBorder && + if (evSel.isNoTimeFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) { return false; } - if (isNoTimeFrameBorder) { + if (evSel.isNoTimeFrameBorder) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isNoITSROFrameBorder && + if (evSel.isNoITSROFrameBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { return false; } - if (isNoITSROFrameBorder) { + if (evSel.isNoITSROFrameBorder) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isVertexITSTPC && + if (evSel.isVertexITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { return false; } - if (isVertexITSTPC) { + if (evSel.isVertexITSTPC) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isNoSameBunchPileup && + if (evSel.isNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { return false; } - if (isNoSameBunchPileup) { + if (evSel.isNoSameBunchPileup) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isGoodZvtxFT0vsPV && + if (evSel.isGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { return false; } - if (isGoodZvtxFT0vsPV) { + if (evSel.isGoodZvtxFT0vsPV) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isTriggerTVX && + if (evSel.isTriggerTVX && !collision.selection_bit(o2::aod::evsel::kIsTriggerTVX)) { return false; } - if (isTriggerTVX) { + if (evSel.isTriggerTVX) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); } } - if (isINELgt0 && (collision.isInelGt0() == false)) { + if (eventClass.applyCentSel && (collision.centFT0M() < eventClass.cfgCentMin || collision.centFT0M() > eventClass.cfgCentMax)) { + return false; + } + + if (evSel.isINELgt0 && (collision.isInelGt0() == false)) { return false; } - if (isINELgt0) { + if (evSel.isINELgt0) { nbinev++; if (fillCounter) { rEventSelection.fill(HIST("hEventsSelected"), nbinev); @@ -940,14 +990,13 @@ struct Lambdak0sflattenicity { } // ============== Flattenicity estimation begins ===================== // - // fillQA=false evaluates the estimators without touching the QA registry, for the - // callers that would otherwise fill it twice for the same collision + // fillQA=false skips the QA registry, for callers that would fill it twice template float estimateFlattenicity(TCollision const& collision, Tracks const& tracks, bool fillQA = true) { - const bool flattenicityQAhere = flattenicityQA && fillQA; - float ampl5[kNeta5] = {0, 0}; - float ampl6[kNeta6] = {0, 0}; + const bool flattenicityQAhere = flatSel.flattenicityQA && fillQA; + std::array ampl5 = {0, 0}; + std::array ampl6 = {0, 0}; auto vtxZ = collision.posZ(); @@ -958,7 +1007,7 @@ struct Lambdak0sflattenicity { ampchannelBefore.fill(0.0); rhoLattice.fill(0); - if ((isflattenicitywithFV0 || isflattenicitywithFV0FT0C) && + if ((flatSel.isflattenicitywithFV0 || flatSel.isflattenicitywithFV0FT0C) && collision.has_foundFV0()) { auto fv0 = collision.foundFV0(); @@ -976,7 +1025,7 @@ struct Lambdak0sflattenicity { phiv0 = ((2.0 * channelv0phi) + 1 - 64.0) * constants::math::TwoPI / (32.0); } ampchannelBefore[channelv0phi] = amplCh; - if (applyCalibCh) { + if (flatSel.applyCalibCh) { amplCh *= kCalib[channelv0phi]; } sumAmpFV0 += amplCh; @@ -998,7 +1047,7 @@ struct Lambdak0sflattenicity { if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpV0vsVtxBeforeCalibration"), vtxZ, sumAmpFV0); } - if (applyCalibVtx) { + if (flatSel.applyCalibVtx) { sumAmpFV0 *= gVtx[0].Eval(vtxZ); sumAmpFV01to4Ch *= gVtx[0].Eval(vtxZ); } @@ -1008,7 +1057,7 @@ struct Lambdak0sflattenicity { } float flattenicityfv0 = 9999; - if (isflattenicitywithFV0 || isflattenicitywithFV0FT0C) { + if (flatSel.isflattenicitywithFV0 || flatSel.isflattenicitywithFV0FT0C) { flattenicityfv0 = getFlatenicity({rhoLattice.data(), rhoLattice.size()}); } @@ -1032,10 +1081,10 @@ struct Lambdak0sflattenicity { rhoLatticeT0A.fill(0); rhoLatticeT0C.fill(0); - if ((isflattenicitywithFT0 || isflattenicitywithFV0FT0C) && + if ((flatSel.isflattenicitywithFT0 || flatSel.isflattenicitywithFV0FT0C) && collision.has_foundFT0()) { auto ft0 = collision.foundFT0(); - if (isflattenicitywithFT0) { + if (flatSel.isflattenicitywithFT0) { for (std::size_t i_a = 0; i_a < ft0.amplitudeA().size(); i_a++) { float amplitude = ft0.amplitudeA()[i_a]; uint8_t channel = ft0.channelA()[i_a]; @@ -1045,7 +1094,7 @@ struct Lambdak0sflattenicity { rFlattenicity.fill(HIST("hAmpT0AVsChBeforeCalibration"), sector, amplitude); } - if (applyCalibCh) { + if (flatSel.applyCalibCh) { amplitude *= kCalibT0A[sector]; } if (flattenicityQAhere) { @@ -1069,7 +1118,7 @@ struct Lambdak0sflattenicity { rFlattenicity.fill(HIST("hAmpT0CVsChBeforeCalibration"), sector, amplitude); } - if (applyCalibCh) { + if (flatSel.applyCalibCh) { amplitude *= kCalibT0C[sector]; } if (flattenicityQAhere) { @@ -1088,7 +1137,7 @@ struct Lambdak0sflattenicity { rFlattenicity.fill(HIST("hAmpT0CvsVtxBeforeCalibration"), vtxZ, sumAmpFT0C); } - if (applyCalibVtx) { + if (flatSel.applyCalibVtx) { sumAmpFT0A *= gVtx[1].Eval(vtxZ); sumAmpFT0C *= gVtx[2].Eval(vtxZ); } @@ -1098,12 +1147,12 @@ struct Lambdak0sflattenicity { } } float flatenicityT0a = 9999; - if (isflattenicitywithFT0) { + if (flatSel.isflattenicitywithFT0) { flatenicityT0a = getFlatenicity({rhoLatticeT0A.data(), rhoLatticeT0A.size()}); } float flatenicityT0c = 9999; - if (isflattenicitywithFT0 || isflattenicitywithFV0FT0C) { + if (flatSel.isflattenicitywithFT0 || flatSel.isflattenicitywithFV0FT0C) { flatenicityT0c = getFlatenicity({rhoLatticeT0C.data(), rhoLatticeT0C.size()}); } @@ -1118,7 +1167,7 @@ struct Lambdak0sflattenicity { } if (!collision.has_foundFV0() || !collision.has_foundFT0()) { - // no forward information: flattenicity is undefined for this collision + // no FV0/FT0, flattenicity undefined return kInvalidFlattenicity; } @@ -1130,7 +1179,7 @@ struct Lambdak0sflattenicity { isOKEstimator5 = true; } if (isOKEstimator5) { - if (applyNorm) { + if (flatSel.applyNorm) { allWeights = 0; for (int i5 = 0; i5 < kNeta5; ++i5) { combinedEstimator5 += @@ -1151,7 +1200,7 @@ struct Lambdak0sflattenicity { isOKEstimator6 = true; } if (isOKEstimator6) { - if (applyNorm) { + if (flatSel.applyNorm) { allWeights = 0; for (int i6 = 0; i6 < kNeta6; ++i6) { combinedEstimator6 += @@ -1207,7 +1256,7 @@ struct Lambdak0sflattenicity { rFlattenicity.fill(HIST(kHPtEst[7]), estimator[7], pt); } - if (isflattenicitywithFV0) { + if (flatSel.isflattenicitywithFV0) { for (int iCh = 0; iCh < kNCells; ++iCh) { rFlattenicity.fill(HIST("hAmpV0VsCh"), iCh, ampchannel[iCh]); rFlattenicity.fill(HIST("hAmpV0VsChBeforeCalibration"), iCh, @@ -1224,22 +1273,50 @@ struct Lambdak0sflattenicity { rFlattenicity.fill(HIST("hFV0amplvsFlat"), sumAmpFV0, estimator[2]); } - if (flattenicityforanalysis == kFlatFromFT0) { + if (flatSel.flattenicityforanalysis == kFlatFromFT0) { finalflattenicity = estimator[4]; } - if (flattenicityforanalysis == kFlatFromFV0FT0C) { + if (flatSel.flattenicityforanalysis == kFlatFromFV0FT0C) { finalflattenicity = estimator[6]; } return finalflattenicity; } + // which weak decay produced a non-primary candidate, bins of hFeedDown*MotherPdg + static constexpr int kFdXiMinus = 0; + static constexpr int kFdXiZero = 1; + static constexpr int kFdOmegaMinus = 2; + static constexpr int kFdOther = 3; + static constexpr int kNFeedDownMothers = 4; + + template + int getFeedDownMother(TMcParticle const& mcParticle, float& motherPt) + { + motherPt = -1.f; + if (!mcParticle.has_mothers()) { + return kFdOther; + } + for (const auto& mother : mcParticle.template mothers_as()) { + motherPt = mother.pt(); + const int motherPdg = std::abs(mother.pdgCode()); + if (motherPdg == PDG_t::kXiMinus) { + return kFdXiMinus; + } + if (motherPdg == o2::constants::physics::Pdg::kXi0) { + return kFdXiZero; + } + if (motherPdg == PDG_t::kOmegaMinus) { + return kFdOmegaMinus; + } + return kFdOther; + } + return kFdOther; + } + template float estimateFlattenicityFV0MC(McParticles const& mcParticles) { rhoLatticeFV0AMC.fill(0); - float flattenicity = -1; - float etamin, etamax, minphi, maxphi, dphi; - int isegment = 0, nsectors; int multFV0 = 0; for (const auto& mcParticle : mcParticles) { @@ -1252,27 +1329,21 @@ struct Lambdak0sflattenicity { continue; } - float etap = mcParticle.eta(); - float phip = mcParticle.phi(); - isegment = 0; + const float etap = mcParticle.eta(); + const float phip = mcParticle.phi(); + int isegment = 0; for (int ieta = 0; ieta < kNFV0EtaRings; ieta++) { - etamax = kMaxEtaFV0 - ieta * kDetaFV0; - if (ieta == 0) { - etamax = kMaxEtaFV0; - } - etamin = kMaxEtaFV0 - (ieta + 1) * kDetaFV0; - if (ieta == kOuterFV0RingIndex) { - etamin = kMinEtaFV0; - } - nsectors = kNChannelsPerFV0Ring; - if (ieta == kOuterFV0RingIndex) { - nsectors = kNSectorsFV0OuterRing; - } + // the outer ring ends exactly at kMinEtaFV0, the rest follow from the ring width + const float etamax = kMaxEtaFV0 - ieta * kDetaFV0; + const float etamin = (ieta == kOuterFV0RingIndex) ? kMinEtaFV0 + : kMaxEtaFV0 - (ieta + 1) * kDetaFV0; + const int nsectors = (ieta == kOuterFV0RingIndex) ? kNSectorsFV0OuterRing + : kNChannelsPerFV0Ring; for (int iphi = 0; iphi < nsectors; iphi++) { - minphi = iphi * constants::math::TwoPI / nsectors; - maxphi = (iphi + 1) * constants::math::TwoPI / nsectors; - dphi = std::abs(maxphi - minphi); + const float minphi = iphi * constants::math::TwoPI / nsectors; + const float maxphi = (iphi + 1) * constants::math::TwoPI / nsectors; + const float dphi = std::abs(maxphi - minphi); if (etap >= etamin && etap < etamax && phip >= minphi && phip < maxphi) { rhoLatticeFV0AMC[isegment] += 1.0 / std::abs(dphi * kDetaFV0); multFV0++; @@ -1282,7 +1353,7 @@ struct Lambdak0sflattenicity { } } - flattenicity = + const float flattenicity = 1.0 - getFlatenicity({rhoLatticeFV0AMC.data(), rhoLatticeFV0AMC.size()}); rEventSelection.fill(HIST("hTrueFV0amplvsFlat"), multFV0, flattenicity); return flattenicity; @@ -1292,25 +1363,25 @@ struct Lambdak0sflattenicity { // Filters on V0s // Cannot filter on dynamic columns, so we cut on DCA to PV and DCA between // daughters only - Filter preFilterV0 = (nabs(aod::v0data::dcapostopv) > v0settingDCApostopv && - nabs(aod::v0data::dcanegtopv) > v0settingDCAnegtopv && - aod::v0data::dcaV0daughters < v0settingDCAv0dau); + Filter preFilterV0 = (nabs(aod::v0data::dcapostopv) > v0Sel.v0settingDCApostopv && + nabs(aod::v0data::dcanegtopv) > v0Sel.v0settingDCAnegtopv && + aod::v0data::dcaV0daughters < v0Sel.v0settingDCAv0dau); Filter trackFilter = - (nabs(aod::track::eta) < cfgTrkEtaCut && aod::track::pt > cfgTrkLowPtCut); + (nabs(aod::track::eta) < trkPid.cfgTrkEtaCut && aod::track::pt > trkPid.cfgTrkLowPtCut); using TrackCandidates = soa::Filtered< soa::Join>; void processDataRun3LambdaK0s( - soa::Join::iterator const& collision, soa::Filtered const& V0s, TrackCandidates const& tracks, soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/) { - if (applyEvSel && + if (evSel.applyEvSel && !(isEventSelected(collision))) { // Checking if the event passes the // selection criteria return; @@ -1328,20 +1399,22 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hVertexZ"), vtxZ); rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); for (const auto& v0 : V0s) { const auto& posDaughterTrack = v0.posTrack_as(); const auto& negDaughterTrack = v0.negTrack_as(); - if (std::abs(posDaughterTrack.eta()) > cfgTrkEtaCut || - std::abs(negDaughterTrack.eta()) > cfgTrkEtaCut || - negDaughterTrack.pt() < cfgTrkLowPtCut || - posDaughterTrack.pt() < cfgTrkLowPtCut) { + if (std::abs(posDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + std::abs(negDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + negDaughterTrack.pt() < trkPid.cfgTrkLowPtCut || + posDaughterTrack.pt() < trkPid.cfgTrkLowPtCut) { continue; } - if (posDaughterTrack.tpcNClsCrossedRows() < v0settingNTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < v0settingNTPCcrossedRows) { + if (posDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows || + negDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows) { continue; } float massK0s = v0.mK0Short(); @@ -1362,7 +1435,7 @@ struct Lambdak0sflattenicity { float v0p = std::sqrt(v0.pt() * v0.pt() + v0.pz() * v0.pz()); float ctauK0s = decaylength * o2::constants::physics::MassK0Short / v0p; - // the Lambda and the anti-Lambda share the PDG mass, so one value serves both + // same PDG mass for Lambda and AntiLambda float ctauLambda = decaylength * o2::constants::physics::MassLambda0 / v0p; float alpha = v0.alpha(); @@ -1371,16 +1444,16 @@ struct Lambdak0sflattenicity { // Cut on dynamic columns for K0s rCommonHist.fill(HIST("hArmPodoAlphavsQT"), alpha, qtarm); - if (v0.v0cosPA() >= v0settingCosPAK0s && - v0.v0radius() >= v0settingRadiusK0s && - std::abs(posDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - std::abs(negDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - ctauK0s < v0settingcTauK0s && - std::abs(v0.rapidity(0)) <= v0settingRapidity && - std::abs(massLambda - o2::constants::physics::MassLambda0) > v0settingMassRejectionK0s && + if (v0.v0cosPA() >= v0Sel.v0settingCosPAK0s && + v0.v0radius() >= v0Sel.v0settingRadiusK0s && + std::abs(posDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + std::abs(negDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + ctauK0s < v0Sel.v0settingcTauK0s && + std::abs(v0.rapidity(0)) <= v0Sel.v0settingRapidity && + std::abs(massLambda - o2::constants::physics::MassLambda0) > v0Sel.v0settingMassRejectionK0s && std::abs(massAntiLambda - o2::constants::physics::MassLambda0) > - v0settingMassRejectionK0s && - qtarm > v0settingArmePodoK0s * std::abs(alpha)) { + v0Sel.v0settingMassRejectionK0s && + qtarm > v0Sel.v0settingArmePodoK0s * std::abs(alpha)) { rKzeroShort.fill(HIST("hMassK0sSelected"), massK0s); rKzeroShort.fill(HIST("hDCAV0DaughtersK0s"), v0.dcaV0daughters()); @@ -1393,7 +1466,7 @@ struct Lambdak0sflattenicity { rKzeroShort.fill(HIST("hArmPodoAlphavsQTK0sAfterCut"), alpha, qtarm); // Filling the PID of the V0 daughters in the region of the K0s peak - if (std::abs(massK0s - o2::constants::physics::MassK0Short) < pidQAWindowK0s) { + if (std::abs(massK0s - o2::constants::physics::MassK0Short) < trkPid.pidQAWindowK0s) { rKzeroShort.fill(HIST("hNSigmaPosPionFromK0s"), posDaughterTrack.tpcNSigmaPi(), posDaughterTrack.tpcInnerParam()); @@ -1404,13 +1477,13 @@ struct Lambdak0sflattenicity { } // Cut on dynamic columns for Lambda - if (v0.v0cosPA() >= v0settingCosPALambda && - v0.v0radius() >= v0settingRadiusLambda && - std::abs(posDaughterTrack.tpcNSigmaPr()) <= nSigmaTPCProton && - std::abs(negDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - ctauLambda < v0settingcTauLambda && - std::abs(v0.rapidity(1)) <= v0settingRapidity && - std::abs(massK0s - o2::constants::physics::MassK0Short) > v0settingMassRejectionLambda) { + if (v0.v0cosPA() >= v0Sel.v0settingCosPALambda && + v0.v0radius() >= v0Sel.v0settingRadiusLambda && + std::abs(posDaughterTrack.tpcNSigmaPr()) <= trkPid.nSigmaTPCProton && + std::abs(negDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + ctauLambda < v0Sel.v0settingcTauLambda && + std::abs(v0.rapidity(1)) <= v0Sel.v0settingRapidity && + std::abs(massK0s - o2::constants::physics::MassK0Short) > v0Sel.v0settingMassRejectionLambda) { rLambda.fill(HIST("hMassLambdaSelected"), massLambda); rLambda.fill(HIST("hDCAV0DaughtersLambda"), v0.dcaV0daughters()); @@ -1422,7 +1495,7 @@ struct Lambdak0sflattenicity { rLambda.fill(HIST("hMassLambdapTFlat"), massLambda, v0.pt(), flattenicity); // Filling the PID of the V0 daughters in the region of the Lambda peak - if (std::abs(massLambda - o2::constants::physics::MassLambda0) < pidQAWindowLambda) { + if (std::abs(massLambda - o2::constants::physics::MassLambda0) < trkPid.pidQAWindowLambda) { rLambda.fill(HIST("hNSigmaPosProtonFromLambda"), posDaughterTrack.tpcNSigmaPr(), posDaughterTrack.tpcInnerParam()); @@ -1433,13 +1506,13 @@ struct Lambdak0sflattenicity { } // Cut on dynamic columns for AntiLambda - if (v0.v0cosPA() >= v0settingCosPALambda && - v0.v0radius() >= v0settingRadiusLambda && - std::abs(posDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - std::abs(negDaughterTrack.tpcNSigmaPr()) <= nSigmaTPCProton && - ctauLambda < v0settingcTauLambda && - std::abs(v0.rapidity(2)) <= v0settingRapidity && - std::abs(massK0s - o2::constants::physics::MassK0Short) > v0settingMassRejectionLambda) { + if (v0.v0cosPA() >= v0Sel.v0settingCosPALambda && + v0.v0radius() >= v0Sel.v0settingRadiusLambda && + std::abs(posDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + std::abs(negDaughterTrack.tpcNSigmaPr()) <= trkPid.nSigmaTPCProton && + ctauLambda < v0Sel.v0settingcTauLambda && + std::abs(v0.rapidity(2)) <= v0Sel.v0settingRapidity && + std::abs(massK0s - o2::constants::physics::MassK0Short) > v0Sel.v0settingMassRejectionLambda) { rAntiLambda.fill(HIST("hMassAntiLambdaSelected"), massAntiLambda); rAntiLambda.fill(HIST("hDCAV0DaughtersAntiLambda"), @@ -1453,7 +1526,7 @@ struct Lambdak0sflattenicity { rAntiLambda.fill(HIST("hMassAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); // Filling the PID of the V0 daughters in the region of the AntiLambda // peak - if (std::abs(massAntiLambda - o2::constants::physics::MassLambda0) < pidQAWindowLambda) { + if (std::abs(massAntiLambda - o2::constants::physics::MassLambda0) < trkPid.pidQAWindowLambda) { rAntiLambda.fill(HIST("hNSigmaPosPionFromAntiLambda"), posDaughterTrack.tpcNSigmaPi(), posDaughterTrack.tpcInnerParam()); @@ -1472,20 +1545,19 @@ struct Lambdak0sflattenicity { Preslice>> perCol = aod::v0data::collisionId; Preslice perColTracksMC = aod::track::collisionId; - // declaring this Preslice is what enables the mcCollisionId slice cache used by - // the sliceByCached calls below; without it they throw "Disabled cache" at run time + // needed by the sliceByCached calls below, else they throw "Disabled cache" Preslice perMCCol = aod::mcparticle::mcCollisionId; SliceCache cache1; void processRecMCLambdaK0s( - soa::Join const& collisions, soa::Filtered> const& V0s, aod::McCollisions const&, TrackCandidatesMC const& tracks, soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/, aod::McParticles const& mcParticles) { for (const auto& collision : collisions) { - if (applyEvSel && + if (evSel.applyEvSel && !(isEventSelected(collision))) { // Checking if the event passes the // selection criteria continue; @@ -1508,6 +1580,8 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hVertexZ"), vtxZ); rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); auto v0sThisCollision = V0s.sliceBy(perCol, collision.globalIndex()); const auto& mcCollision = collision.mcCollision_as(); @@ -1517,15 +1591,15 @@ struct Lambdak0sflattenicity { const auto& posDaughterTrack = v0.posTrack_as(); const auto& negDaughterTrack = v0.negTrack_as(); - if (std::abs(posDaughterTrack.eta()) > cfgTrkEtaCut || - std::abs(negDaughterTrack.eta()) > cfgTrkEtaCut || - negDaughterTrack.pt() < cfgTrkLowPtCut || - posDaughterTrack.pt() < cfgTrkLowPtCut) { + if (std::abs(posDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + std::abs(negDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + negDaughterTrack.pt() < trkPid.cfgTrkLowPtCut || + posDaughterTrack.pt() < trkPid.cfgTrkLowPtCut) { continue; } - if (posDaughterTrack.tpcNClsCrossedRows() < v0settingNTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < v0settingNTPCcrossedRows) { + if (posDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows || + negDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows) { continue; } @@ -1551,7 +1625,7 @@ struct Lambdak0sflattenicity { float v0p = std::sqrt(v0.pt() * v0.pt() + v0.pz() * v0.pz()); float ctauK0s = decaylength * o2::constants::physics::MassK0Short / v0p; - // the Lambda and the anti-Lambda share the PDG mass, so one value serves both + // same PDG mass for Lambda and AntiLambda float ctauLambda = decaylength * o2::constants::physics::MassLambda0 / v0p; float alpha = v0.alpha(); @@ -1559,18 +1633,22 @@ struct Lambdak0sflattenicity { rCommonHist.fill(HIST("hArmPodoAlphavsQT"), alpha, qtarm); auto v0mcParticle = v0.mcParticle(); + const bool isPrimaryV0 = v0mcParticle.isPhysicalPrimary(); + const bool keepForEfficiency = isPrimaryV0 || !eventClass.requirePrimaryMC; + float motherPt = -1.f; // Cut on dynamic columns for K0s - if (v0mcParticle.pdgCode() == PDG_t::kK0Short && v0.v0cosPA() >= v0settingCosPAK0s && - v0.v0radius() >= v0settingRadiusK0s && - std::abs(posDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - std::abs(negDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - ctauK0s < v0settingcTauK0s && - std::abs(v0.rapidity(0)) <= v0settingRapidity && - std::abs(massLambda - o2::constants::physics::MassLambda0) > v0settingMassRejectionK0s && + if (v0mcParticle.pdgCode() == PDG_t::kK0Short && keepForEfficiency && + v0.v0cosPA() >= v0Sel.v0settingCosPAK0s && + v0.v0radius() >= v0Sel.v0settingRadiusK0s && + std::abs(posDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + std::abs(negDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + ctauK0s < v0Sel.v0settingcTauK0s && + std::abs(v0.rapidity(0)) <= v0Sel.v0settingRapidity && + std::abs(massLambda - o2::constants::physics::MassLambda0) > v0Sel.v0settingMassRejectionK0s && std::abs(massAntiLambda - o2::constants::physics::MassLambda0) > - v0settingMassRejectionK0s && - qtarm > v0settingArmePodoK0s * std::abs(alpha)) { + v0Sel.v0settingMassRejectionK0s && + qtarm > v0Sel.v0settingArmePodoK0s * std::abs(alpha)) { rKzeroShort.fill(HIST("hMassK0sSelected"), massK0s); rKzeroShort.fill(HIST("hDCAV0DaughtersK0s"), v0.dcaV0daughters()); @@ -1583,7 +1661,7 @@ struct Lambdak0sflattenicity { rKzeroShort.fill(HIST("hArmPodoAlphavsQTK0sAfterCut"), alpha, qtarm); // Filling the PID of the V0 daughters in the region of the K0s peak - if (std::abs(massK0s - o2::constants::physics::MassK0Short) < pidQAWindowK0s) { + if (std::abs(massK0s - o2::constants::physics::MassK0Short) < trkPid.pidQAWindowK0s) { rKzeroShort.fill(HIST("hNSigmaPosPionFromK0s"), posDaughterTrack.tpcNSigmaPi(), posDaughterTrack.tpcInnerParam()); @@ -1595,63 +1673,83 @@ struct Lambdak0sflattenicity { // Cut on dynamic columns for Lambda if (v0mcParticle.pdgCode() == PDG_t::kLambda0 && - v0.v0cosPA() >= v0settingCosPALambda && - v0.v0radius() >= v0settingRadiusLambda && - std::abs(posDaughterTrack.tpcNSigmaPr()) <= nSigmaTPCProton && - std::abs(negDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - ctauLambda < v0settingcTauLambda && - std::abs(v0.rapidity(1)) <= v0settingRapidity && - std::abs(massK0s - o2::constants::physics::MassK0Short) > v0settingMassRejectionLambda) { - - rLambda.fill(HIST("hMassLambdaSelected"), massLambda); - rLambda.fill(HIST("hDCAV0DaughtersLambda"), v0.dcaV0daughters()); - rLambda.fill(HIST("hV0CosPALambda"), v0.v0cosPA()); - rLambda.fill(HIST("hrapidityLambda"), v0.rapidity(1)); - rLambda.fill(HIST("hctauLambda"), ctauLambda); - rLambda.fill(HIST("h2DdecayRadiusLambda"), v0.v0radius()); - rLambda.fill(HIST("hMassLambdapT"), massLambda, v0.pt()); - rLambda.fill(HIST("hMassLambdapTFlat"), massLambda, v0.pt(), flattenicity); - - // Filling the PID of the V0 daughters in the region of the Lambda peak - if (std::abs(massLambda - o2::constants::physics::MassLambda0) < pidQAWindowLambda) { - rLambda.fill(HIST("hNSigmaPosProtonFromLambda"), - posDaughterTrack.tpcNSigmaPr(), - posDaughterTrack.tpcInnerParam()); - rLambda.fill(HIST("hNSigmaNegPionFromLambda"), - negDaughterTrack.tpcNSigmaPi(), - negDaughterTrack.tpcInnerParam()); + v0.v0cosPA() >= v0Sel.v0settingCosPALambda && + v0.v0radius() >= v0Sel.v0settingRadiusLambda && + std::abs(posDaughterTrack.tpcNSigmaPr()) <= trkPid.nSigmaTPCProton && + std::abs(negDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + ctauLambda < v0Sel.v0settingcTauLambda && + std::abs(v0.rapidity(1)) <= v0Sel.v0settingRapidity && + std::abs(massK0s - o2::constants::physics::MassK0Short) > v0Sel.v0settingMassRejectionLambda) { + + // secondaries kept out of the efficiency numerator, booked for the feed-down + if (!isPrimaryV0) { + const int motherIndex = getFeedDownMother(v0mcParticle, motherPt); + rLambda.fill(HIST("hMassFeedDownLambdapTFlat"), massLambda, v0.pt(), flattenicity); + rLambda.fill(HIST("hFeedDownLambdaPtVsMotherPt"), v0.pt(), motherPt); + rLambda.fill(HIST("hFeedDownLambdaMotherPdg"), motherIndex); + } + if (keepForEfficiency) { + + rLambda.fill(HIST("hMassLambdaSelected"), massLambda); + rLambda.fill(HIST("hDCAV0DaughtersLambda"), v0.dcaV0daughters()); + rLambda.fill(HIST("hV0CosPALambda"), v0.v0cosPA()); + rLambda.fill(HIST("hrapidityLambda"), v0.rapidity(1)); + rLambda.fill(HIST("hctauLambda"), ctauLambda); + rLambda.fill(HIST("h2DdecayRadiusLambda"), v0.v0radius()); + rLambda.fill(HIST("hMassLambdapT"), massLambda, v0.pt()); + rLambda.fill(HIST("hMassLambdapTFlat"), massLambda, v0.pt(), flattenicity); + + // Filling the PID of the V0 daughters in the region of the Lambda peak + if (std::abs(massLambda - o2::constants::physics::MassLambda0) < trkPid.pidQAWindowLambda) { + rLambda.fill(HIST("hNSigmaPosProtonFromLambda"), + posDaughterTrack.tpcNSigmaPr(), + posDaughterTrack.tpcInnerParam()); + rLambda.fill(HIST("hNSigmaNegPionFromLambda"), + negDaughterTrack.tpcNSigmaPi(), + negDaughterTrack.tpcInnerParam()); + } } } // Cut on dynamic columns for AntiLambda if (v0mcParticle.pdgCode() == PDG_t::kLambda0Bar && - v0.v0cosPA() >= v0settingCosPALambda && - v0.v0radius() >= v0settingRadiusLambda && - std::abs(posDaughterTrack.tpcNSigmaPi()) <= nSigmaTPCPion && - std::abs(negDaughterTrack.tpcNSigmaPr()) <= nSigmaTPCProton && - ctauLambda < v0settingcTauLambda && - std::abs(v0.rapidity(2)) <= v0settingRapidity && - std::abs(massK0s - o2::constants::physics::MassK0Short) > v0settingMassRejectionLambda) { - - rAntiLambda.fill(HIST("hMassAntiLambdaSelected"), massAntiLambda); - rAntiLambda.fill(HIST("hDCAV0DaughtersAntiLambda"), - v0.dcaV0daughters()); - rAntiLambda.fill(HIST("hV0CosPAAntiLambda"), v0.v0cosPA()); - rAntiLambda.fill(HIST("hrapidityAntiLambda"), v0.rapidity(2)); - rAntiLambda.fill(HIST("hctauAntiLambda"), ctauLambda); - rAntiLambda.fill(HIST("h2DdecayRadiusAntiLambda"), v0.v0radius()); - rAntiLambda.fill(HIST("hMassAntiLambdapT"), massAntiLambda, v0.pt()); - rAntiLambda.fill(HIST("hMassAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); - - // Filling the PID of the V0 daughters in the region of the AntiLambda - // peak - if (std::abs(massAntiLambda - o2::constants::physics::MassLambda0) < pidQAWindowLambda) { - rAntiLambda.fill(HIST("hNSigmaPosPionFromAntiLambda"), - posDaughterTrack.tpcNSigmaPi(), - posDaughterTrack.tpcInnerParam()); - rAntiLambda.fill(HIST("hNSigmaNegProtonFromAntiLambda"), - negDaughterTrack.tpcNSigmaPr(), - negDaughterTrack.tpcInnerParam()); + v0.v0cosPA() >= v0Sel.v0settingCosPALambda && + v0.v0radius() >= v0Sel.v0settingRadiusLambda && + std::abs(posDaughterTrack.tpcNSigmaPi()) <= trkPid.nSigmaTPCPion && + std::abs(negDaughterTrack.tpcNSigmaPr()) <= trkPid.nSigmaTPCProton && + ctauLambda < v0Sel.v0settingcTauLambda && + std::abs(v0.rapidity(2)) <= v0Sel.v0settingRapidity && + std::abs(massK0s - o2::constants::physics::MassK0Short) > v0Sel.v0settingMassRejectionLambda) { + + // secondaries kept out of the efficiency numerator, booked for the feed-down + if (!isPrimaryV0) { + const int motherIndex = getFeedDownMother(v0mcParticle, motherPt); + rAntiLambda.fill(HIST("hMassFeedDownAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); + rAntiLambda.fill(HIST("hFeedDownAntiLambdaPtVsMotherPt"), v0.pt(), motherPt); + rAntiLambda.fill(HIST("hFeedDownAntiLambdaMotherPdg"), motherIndex); + } + if (keepForEfficiency) { + + rAntiLambda.fill(HIST("hMassAntiLambdaSelected"), massAntiLambda); + rAntiLambda.fill(HIST("hDCAV0DaughtersAntiLambda"), + v0.dcaV0daughters()); + rAntiLambda.fill(HIST("hV0CosPAAntiLambda"), v0.v0cosPA()); + rAntiLambda.fill(HIST("hrapidityAntiLambda"), v0.rapidity(2)); + rAntiLambda.fill(HIST("hctauAntiLambda"), ctauLambda); + rAntiLambda.fill(HIST("h2DdecayRadiusAntiLambda"), v0.v0radius()); + rAntiLambda.fill(HIST("hMassAntiLambdapT"), massAntiLambda, v0.pt()); + rAntiLambda.fill(HIST("hMassAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); + + // Filling the PID of the V0 daughters in the region of the AntiLambda + // peak + if (std::abs(massAntiLambda - o2::constants::physics::MassLambda0) < trkPid.pidQAWindowLambda) { + rAntiLambda.fill(HIST("hNSigmaPosPionFromAntiLambda"), + posDaughterTrack.tpcNSigmaPi(), + posDaughterTrack.tpcInnerParam()); + rAntiLambda.fill(HIST("hNSigmaNegProtonFromAntiLambda"), + negDaughterTrack.tpcNSigmaPr(), + negDaughterTrack.tpcInnerParam()); + } } } } @@ -1694,20 +1792,18 @@ struct Lambdak0sflattenicity { } } - // Filter posZFilterMC = (nabs(o2::aod::mccollision::posZ) < cutzvertex); + // Filter posZFilterMC = (nabs(o2::aod::mccollision::posZ) < evSel.cutzvertex); void processGenMC( - o2::aod::McCollision const& mcCollision, const soa::SmallGroups>& collisions, TrackCandidatesMC const& tracks, aod::FT0s const& /*ft0s*/, + o2::aod::McCollision const& mcCollision, const soa::SmallGroups>& collisions, TrackCandidatesMC const& tracks, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/, o2::aod::McParticles const& mcParticles) { - // A generated collision without a usable reconstructed counterpart has no - // reconstructed flattenicity. It keeps the sentinel value, which lands in the - // underflow of the flattenicity axis, so the event- and signal-loss counters - // stay complete while nothing contaminates the physical range. + // without a reconstructed counterpart the sentinel is kept, it falls in the + // underflow so the loss counters stay complete float flattenicity = kInvalidFlattenicity; - if (flattenicityforLossCorrRec) { + if (flatSel.flattenicityforLossCorrRec) { for (const auto& collision : collisions) { - if (applyEvSel && !isEventSelected(collision, false)) { + if (evSel.applyEvSel && !isEventSelected(collision, false)) { continue; } auto tracksThisCollision = tracks.sliceBy(perColTracksMC, collision.globalIndex()); @@ -1728,7 +1824,7 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hNEventsMCGen"), 0.5); - if (std::abs(mcCollision.posZ()) > cutzvertex) { + if (std::abs(mcCollision.posZ()) > evSel.cutzvertex) { return; } rEventSelection.fill(HIST("hNEventsMCGen"), 1.5); @@ -1794,7 +1890,7 @@ struct Lambdak0sflattenicity { // process functions; filling it here as well would double count every // collision whenever processGenMC runs alongside one of them, which is the // normal MC configuration. hNEventsMCReco already records all/passed here. - if (applyEvSel && !isEventSelected(collision, false)) { + if (evSel.applyEvSel && !isEventSelected(collision, false)) { continue; } rEventSelection.fill(HIST("hEventsSelected"), nbin - 0.5); @@ -1915,8 +2011,7 @@ struct Lambdak0sflattenicity { static constexpr float kCtauXi = 4.91; // Xi lifetime, in cm - // In #Xi^{-} #rightarrow #Lambda #pi^{-} the baryon daughter is the positive one, - // and the other way around for the #bar{#Xi}^{+}. The bachelor is always a pion. + // for Xi- the baryon daughter is the positive one, the other way around for Xi+ template bool isSelectedXi(TCollision const& collision, TCascade const& casc, TTrack const& posDaughterTrack, TTrack const& negDaughterTrack, @@ -1926,22 +2021,22 @@ struct Lambdak0sflattenicity { const auto& protonDaughter = isXiMinus ? posDaughterTrack : negDaughterTrack; const auto& pionDaughter = isXiMinus ? negDaughterTrack : posDaughterTrack; - // PID: one hypothesis per daughter, two sided - if (std::abs(protonDaughter.tpcNSigmaPr()) > nSigmaTPCProton || - std::abs(pionDaughter.tpcNSigmaPi()) > nSigmaTPCPion || - std::abs(bacDaughterTrack.tpcNSigmaPi()) > nSigmaTPCPion) { + // PID, one hypothesis per daughter + if (std::abs(protonDaughter.tpcNSigmaPr()) > trkPid.nSigmaTPCProton || + std::abs(pionDaughter.tpcNSigmaPi()) > trkPid.nSigmaTPCPion || + std::abs(bacDaughterTrack.tpcNSigmaPi()) > trkPid.nSigmaTPCPion) { return false; } // track quality - if (posDaughterTrack.tpcNClsCrossedRows() < nTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < nTPCcrossedRows || - bacDaughterTrack.tpcNClsCrossedRows() < nTPCcrossedRows) { + if (posDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows || + negDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows || + bacDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows) { return false; } - if (std::abs(posDaughterTrack.eta()) > cfgTrkEtaCut || - std::abs(negDaughterTrack.eta()) > cfgTrkEtaCut || - std::abs(bacDaughterTrack.eta()) > cfgTrkEtaCut) { + if (std::abs(posDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + std::abs(negDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || + std::abs(bacDaughterTrack.eta()) > trkPid.cfgTrkEtaCut) { return false; } @@ -1950,46 +2045,45 @@ struct Lambdak0sflattenicity { const float cosPAcasc = casc.casccosPA(collision.posX(), collision.posY(), collision.posZ()); const float cosPAv0 = casc.v0cosPA(collision.posX(), collision.posY(), collision.posZ()); - // the DCA to the PV is cut per daughter role, not per charge, so that the - // #Xi^{-} and the #bar{#Xi}^{+} are treated symmetrically + // DCA to PV per daughter role, not per charge const float dcaBaryonToPV = isXiMinus ? casc.dcapostopv() : casc.dcanegtopv(); const float dcaMesonToPV = isXiMinus ? casc.dcanegtopv() : casc.dcapostopv(); - if (std::abs(dcaBaryonToPV) < cascsettingDCAbaryontopv || std::abs(dcaMesonToPV) < cascsettingDCAmesontopv || - std::abs(casc.dcabachtopv()) < v0settingDCAbactopv || casc.dcaV0daughters() > v0settingDCAv0dau || - dcav0pv < cascsettingDCAv0toPV || casc.dcacascdaughters() > cascsettingDCAv0bach) { + if (std::abs(dcaBaryonToPV) < cascSel.cascsettingDCAbaryontopv || std::abs(dcaMesonToPV) < cascSel.cascsettingDCAmesontopv || + std::abs(casc.dcabachtopv()) < v0Sel.v0settingDCAbactopv || casc.dcaV0daughters() > v0Sel.v0settingDCAv0dau || + dcav0pv < cascSel.cascsettingDCAv0toPV || casc.dcacascdaughters() > cascSel.cascsettingDCAv0bach) { return false; } - // bachelor-baryon veto: reject the candidates where the two are compatible - if (casc.bachBaryonCosPA() > cascsettingCosPAbaryonbach || - std::abs(casc.bachBaryonDCAxyToPV()) < cascsettingDCAxybaryonbach) { + // bachelor-baryon veto + if (casc.bachBaryonCosPA() > cascSel.cascsettingCosPAbaryonbach || + std::abs(casc.bachBaryonDCAxyToPV()) < cascSel.cascsettingDCAxybaryonbach) { return false; } - if (cosPAcasc < cascsettingCosPAcascPV || cosPAv0 < cascsettingCosPAv0PV || - casc.cascradius() < cascsettingcascradius || casc.v0radius() < cascsettingv0radius) { + if (cosPAcasc < cascSel.cascsettingCosPAcascPV || cosPAv0 < cascSel.cascsettingCosPAv0PV || + casc.cascradius() < cascSel.cascsettingcascradius || casc.v0radius() < cascSel.cascsettingv0radius) { return false; } // kinematics and competing species - if (std::abs(casc.yXi()) > cascsettingRapidity || ctauXi > kCtauXi * cascsettingproplifetime) { + if (std::abs(casc.yXi()) > cascSel.cascsettingRapidity || ctauXi > kCtauXi * cascSel.cascsettingproplifetime) { return false; } - if (std::abs(casc.mLambda() - o2::constants::physics::MassLambda0) > cascsettingMassRejectionLambdaXi || - std::abs(casc.mOmega() - o2::constants::physics::MassOmegaMinus) < cascsettingMassRejectioOmegaXi) { + if (std::abs(casc.mLambda() - o2::constants::physics::MassLambda0) > cascSel.cascsettingMassRejectionLambdaXi || + std::abs(casc.mOmega() - o2::constants::physics::MassOmegaMinus) < cascSel.cascsettingMassRejectioOmegaXi) { return false; } return true; } - void processDataRun3Cascade(soa::Join::iterator const& collision, aod::CascDataExt const& Cascades, aod::V0Datas const&, DauTracks const& tracks, soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/) { - if (applyEvSel && + if (evSel.applyEvSel && !(isEventSelected(collision))) { // Checking if the event passes the // selection criteria return; @@ -2007,6 +2101,8 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hVertexZ"), vtxZ); rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); for (const auto& casc : Cascades) { @@ -2045,7 +2141,7 @@ struct Lambdak0sflattenicity { Preslice perColDauTracksMC = aod::track::collisionId; SliceCache cacheCasc; - void processRecMCRun3Cascade(soa::Join const& collisions, LabeledCascades const& Cascades, aod::V0Datas const&, LabeledDauTracks const& tracks, @@ -2053,7 +2149,7 @@ struct Lambdak0sflattenicity { aod::FV0As const& /*fv0s*/, aod::McCollisions const&, aod::McParticles const& mcParticles) { for (const auto& collision : collisions) { - if (applyEvSel && + if (evSel.applyEvSel && !(isEventSelected(collision))) { // Checking if the event passes the // selection criteria continue; @@ -2076,15 +2172,15 @@ struct Lambdak0sflattenicity { rEventSelection.fill(HIST("hVertexZ"), vtxZ); rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); auto cascsThisCollision = Cascades.sliceBy(perColCasc, collision.globalIndex()); const auto& mcCollision = collision.mcCollision_as(); for (const auto& casc : cascsThisCollision) { - // MC truth matching, mirroring what the V0 branch does: without it the - // reconstructed spectra keep the combinatorial background and any - // efficiency built from them comes out too high + // MC truth matching, else the spectra keep the combinatorial background if (!casc.has_mcParticle()) { continue; } @@ -2104,6 +2200,18 @@ struct Lambdak0sflattenicity { if (std::abs(cascMcParticle.pdgCode()) == PDG_t::kXiMinus && isSelectedXi(collision, casc, posDaughterTrack, negDaughterTrack, bacDaughterTrack, ctauXi)) { + const bool isPrimaryCasc = cascMcParticle.isPhysicalPrimary(); + if (!isPrimaryCasc) { + float motherPt = -1.f; + const int motherIndex = getFeedDownMother(cascMcParticle, motherPt); + rXi.fill(HIST("hMassFeedDownXipTFlat"), massXi, casc.pt(), flattenicity); + rXi.fill(HIST("hFeedDownXiPtVsMotherPt"), casc.pt(), motherPt); + rXi.fill(HIST("hFeedDownXiMotherPdg"), motherIndex); + } + if (!isPrimaryCasc && eventClass.requirePrimaryMC) { + continue; + } + const bool isXiMinus = casc.sign() < 0; const auto& protonDaughter = isXiMinus ? posDaughterTrack : negDaughterTrack; const auto& pionDaughter = isXiMinus ? negDaughterTrack : posDaughterTrack;