diff --git a/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx b/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx index ec76c05b8e1..d476833fa64 100644 --- a/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx +++ b/PWGLF/Tasks/Strangeness/lambdak0sflattenicity.cxx @@ -27,6 +27,7 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include #include #include #include @@ -40,8 +41,9 @@ #include #include -#include +#include #include +#include #include #include @@ -100,6 +102,12 @@ struct Lambdak0sflattenicity { OutputObjHandlingPolicy::AnalysisObject, true, true}; + HistogramRegistry rCharged{ + "charged", + {}, + OutputObjHandlingPolicy::AnalysisObject, + true, + true}; static constexpr int kNEstimators = 8; // forward detector segmentation @@ -112,8 +120,6 @@ struct Lambdak0sflattenicity { static constexpr int kNSectorsFV0OuterRing = 16; // TParticlePDG::Charge() is in units of e/3 static constexpr float kMinCharge = 0.01f; - // rapidity window for generated particles - static constexpr float kMcRapidityWindow = 0.5f; static constexpr std::array kHEst = { "eGlobaltrack", "eFV0", "e1flatencityFV0", "eFT0", @@ -127,7 +133,11 @@ struct Lambdak0sflattenicity { "ptVs1flatencityFT0", "ptVsFV0FT0C", "ptVs1flatencityFV0FT0C", "pTVsPtTrig"}; - // Histogram binning + // Histogram binning. The pT and 1-rho axes are ConfigurableAxis: pass + // {nbins, lo, hi} for uniform bins or {VARIABLE_WIDTH, e0, e1, ...} for the + // edges of the published binning. The 1-rho edges are meant to be the + // percentile class boundaries measured by the processFlatDist* pass, so that + // one class is exactly one bin and no rebinning happens downstream. struct : ConfigurableGroup { std::string prefix = "binning"; Configurable nBinsVz{"nBinsVz", 100, "N bins in Vz"}; @@ -138,8 +148,21 @@ struct Lambdak0sflattenicity { 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"}; + ConfigurableAxis axisPtK0s{"axisPtK0s", {250, 0.0f, 25.0f}, "#it{p}_{T} of K0s"}; + ConfigurableAxis axisPtLambda{"axisPtLambda", {250, 0.0f, 25.0f}, "#it{p}_{T} of Lambda and AntiLambda"}; + ConfigurableAxis axisPtXi{"axisPtXi", {250, 0.0f, 25.0f}, "#it{p}_{T} of Xi"}; + ConfigurableAxis axisPtCharged{"axisPtCharged", {250, 0.0f, 25.0f}, "#it{p}_{T} of charged particles"}; + ConfigurableAxis axisPtPid{"axisPtPid", {250, 0.0f, 25.0f}, "#it{p}_{TPC} of the daughter tracks"}; + ConfigurableAxis axisFlat{"axisFlat", {100, 0.0f, 1.0f}, "1-#rho_{ch} of the spectra"}; + // detector effects move the percentiles, so the closure denominator needs + // its own class edges, taken from hFlatDistGen + ConfigurableAxis axisFlatTrue{"axisFlatTrue", {100, 0.0f, 1.0f}, "true 1-#rho_{ch}"}; + ConfigurableAxis axisFlatFine{"axisFlatFine", {2000, 0.0f, 1.0f}, "1-#rho_{ch} used to locate the percentiles"}; + ConfigurableAxis axisCent{"axisCent", {100, 0.0f, 100.0f}, "FT0M percentile"}; + ConfigurableAxis axisNch{"axisNch", {150, -0.5f, 149.5f}, "N_{ch} in the tracking acceptance"}; + ConfigurableAxis axisDcaXy{"axisDcaXy", {200, -1.0f, 1.0f}, "DCA_{xy} (cm)"}; + ConfigurableAxis axisDcaV0ToPv{"axisDcaV0ToPv", {200, 0.0f, 2.0f}, "DCA of the V0 to the PV (cm)"}; + ConfigurableAxis axisPtRes{"axisPtRes", {200, -0.5f, 0.5f}, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen})/#it{p}_{T}^{gen}"}; } binning; // Event selection @@ -170,9 +193,12 @@ struct Lambdak0sflattenicity { struct : ConfigurableGroup { std::string prefix = "flatSel"; Configurable flattenicityQA{"flattenicityQA", true, "Store Flattenicity QA plots"}; - Configurable applyCalibCh{"applyCalibCh", false, "equalize FV0"}; + // Both corrections reshape the lattice and therefore move the flattenicity. + // With both off nothing is read from CCDB at all. + Configurable applyCalibCh{"applyCalibCh", false, + "equalize the per-channel gain, from CCDB"}; Configurable applyCalibVtx{"applyCalibVtx", false, - "equalize FV0 vs vtx"}; + "equalize the per-cell z-vertex dependence, from CCDB"}; Configurable applyNorm{"applyNorm", false, "normalization to eta"}; Configurable isflattenicitywithFV0{"isflattenicitywithFV0", true, "Calculate Flattenicity with FV0"}; @@ -184,8 +210,27 @@ struct Lambdak0sflattenicity { "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"}; + // same cell convention as the detector lattice, so that gen vs rec is a + // detector effect and not a difference of definitions + Configurable genFlatDetectorLikeNorm{"genFlatDetectorLikeNorm", true, + "Weight the generator-level FV0 cells the way the detector lattice does"}; } flatSel; + // Calibration objects, read per run and only when the corresponding switch is + // on. Same layout as PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx, so + // the objects are interchangeable between the two tasks: the gain is a + // std::vector indexed by the raw detector channel, the z-vertex + // equalization a TProfile2D of (lattice cell, z_vtx). + struct : ConfigurableGroup { + std::string prefix = "ccdbConf"; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", + "url of the ccdb repository"}; + Configurable gainEqPath{"gainEqPath", "Users/s/sprasad/flattenicity/GainEq", + "CCDB directory holding FV0, FT0A and FT0C gain vectors"}; + Configurable vtxEqPath{"vtxEqPath", "Users/s/sprasad/flattenicity/ZvtxEq", + "CCDB directory holding the FV0, FT0A and FT0C z-vertex maps"}; + } ccdbConf; + // V0 selection struct : ConfigurableGroup { std::string prefix = "v0Sel"; @@ -233,6 +278,12 @@ struct Lambdak0sflattenicity { "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"}; + // negative values disable a cut, so the defaults reproduce the previous selection + Configurable minCrossedRowsOverFindable{"minCrossedRowsOverFindable", -1.f, + "Minimum TPC crossed rows over findable clusters"}; + Configurable maxTpcChi2NCl{"maxTpcChi2NCl", -1.f, "Maximum TPC chi2 per cluster"}; + Configurable maxItsChi2NCl{"maxItsChi2NCl", -1.f, "Maximum ITS chi2 per cluster"}; + Configurable minItsNCls{"minItsNCls", -1, "Minimum number of ITS clusters"}; } trkPid; // Cascade selection @@ -269,12 +320,25 @@ struct Lambdak0sflattenicity { "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"}; + // the processFlatDist* pass fills the MB and the HM 1-rho distributions + // together, so this window is applied to a histogram and not to the event + Configurable cfgCentMaxHM{"cfgCentMaxHM", 1.0f, + "Upper FT0M percentile of the high-multiplicity class"}; // 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"}; + Configurable genFlatPrimariesOnly{"genFlatPrimariesOnly", true, + "Use only primaries in the generator-level flattenicity"}; + Configurable fillChargedQA{"fillChargedQA", true, + "Fill the charged-particle histograms used for and Qpp"}; + Configurable cfgEtaChargedCut{"cfgEtaChargedCut", 0.8f, + "Eta window of the charged-particle measurement"}; + // mothers outside the Lambda window still feed it, so the matrix is + // normalised in a wider one + Configurable cfgFeedDownMotherRapidity{"cfgFeedDownMotherRapidity", 1.5f, + "Rapidity window of the generated feed-down mothers"}; } eventClass; // Configurable v0daughter_etacut{"V0DaughterEtaCut", 0.8, @@ -287,10 +351,26 @@ struct Lambdak0sflattenicity { // do not fill the detector QA twice when processGenMC runs with a rec-level process bool fillFlattenicityQAInGenMC = true; - // vertex equalization curves, built once in init() - static constexpr int kNDetVtx = 3; - static constexpr int kNVtxPoints = 30; - std::array gVtx; + // the per-event histograms are shared by every rec-level process function, + // so only the one designated in init() fills them + static constexpr int kOwnerNone = -1; + static constexpr int kOwnerFlatDistData = 0; + static constexpr int kOwnerFlatDistMC = 1; + static constexpr int kOwnerRecMCV0 = 2; + static constexpr int kOwnerDataV0 = 3; + static constexpr int kOwnerRecMCCasc = 4; + static constexpr int kOwnerDataCasc = 5; + int eventHistOwner = kOwnerNone; + + // calibration objects, refreshed when the run changes + Service ccdb{}; + int mRunNumber = -1; + std::vector gainFV0; + std::vector gainFT0A; + std::vector gainFT0C; + TProfile2D* vtxEqFV0 = nullptr; + TProfile2D* vtxEqFT0A = nullptr; + TProfile2D* vtxEqFT0C = nullptr; void init(InitContext const&) { @@ -308,24 +388,30 @@ struct Lambdak0sflattenicity { o2::constants::physics::MassXiMinus + binning.kXiEPshiftfromMass, "#it{M}_{#Lambda#pi} [GeV/#it{c}^{2}]"}; 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 ptK0sAxis = {binning.axisPtK0s, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec ptLambdaAxis = {binning.axisPtLambda, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec ptXiAxis = {binning.axisPtXi, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec ptChargedAxis = {binning.axisPtCharged, "#it{p}_{T} (GeV/#it{c})"}; + AxisSpec pTPCAxis = {binning.axisPtPid, "#it{p}_{TPC} (GeV/#it{c})"}; AxisSpec decayRadiusAxis = {100, 0.0f, 100.0f, "Decay Radius (cm)"}; - AxisSpec flatAxis = {binning.nBinsFlattenicity, 0.0f, 1.0f, "1-#rho_{ch}"}; - AxisSpec centAxis = {eventClass.nBinsCent, 0.0f, 100.0f, "FT0M percentile"}; + AxisSpec flatAxis = {binning.axisFlat, "1-#rho_{ch}"}; + AxisSpec flatTrueAxis = {binning.axisFlatTrue, "true 1-#rho_{ch}"}; + AxisSpec flatFineAxis = {binning.axisFlatFine, "1-#rho_{ch}"}; + AxisSpec centAxis = {binning.axisCent, "FT0M percentile"}; + AxisSpec nchAxis = {binning.axisNch, "#it{N}_{ch}"}; + AxisSpec dcaXyAxis = {binning.axisDcaXy, "DCA_{xy} (cm)"}; + AxisSpec dcaV0ToPvAxis = {binning.axisDcaV0ToPv, "DCA_{V0-PV} (cm)"}; + AxisSpec ptResAxis = {binning.axisPtRes, "(#it{p}_{T}^{rec} - #it{p}_{T}^{gen})/#it{p}_{T}^{gen}"}; + AxisSpec motherAxis = {kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}; 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}; - gVtx[0].SetName("gAmpV0"); - gVtx[1].SetName("gAmpT0A"); - gVtx[2].SetName("gAmpT0C"); - for (int iVtx = 0; iVtx < kNVtxPoints; ++iVtx) { - gVtx[0].SetPoint(iVtx, kBiningVtxt[iVtx], kCalibFV0vtx[iVtx]); - gVtx[1].SetPoint(iVtx, kBiningVtxt[iVtx], kCalibFT0Avtx[iVtx]); - gVtx[2].SetPoint(iVtx, kBiningVtxt[iVtx], kCalibFT0Cvtx[iVtx]); - } + ccdb->setURL(ccdbConf.ccdbUrl.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); // Histograms // Event selection @@ -364,7 +450,8 @@ struct Lambdak0sflattenicity { } // events without FV0/FT0 are rejected below, they need their own counter if (doprocessDataRun3LambdaK0s || doprocessRecMCLambdaK0s || - doprocessDataRun3Cascade || doprocessRecMCRun3Cascade) { + doprocessDataRun3Cascade || doprocessRecMCRun3Cascade || + doprocessFlatDistData || doprocessFlatDistMC) { nbinFlattenicity = nbin; rEventSelection.get(HIST("hEventsSelected"))->GetXaxis()->SetBinLabel(nbin++, "flattenicity"); } @@ -374,15 +461,69 @@ struct Lambdak0sflattenicity { rEventSelection.add("hCentFT0M", "hCentFT0M", {HistType::kTH1D, {centAxis}}); rEventSelection.add("hCentFT0MvsFlattenicity", "hCentFT0MvsFlattenicity", {HistType::kTH2D, {centAxis, flatAxis}}); - if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessGenMC) { + + // 1-rho on the fine axis, the input to the percentile boundaries of + // binning.axisFlat. Filled once per accepted collision, by the owner. + rEventSelection.add("hFlatDistRec", "hFlatDistRec", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistRecINELgt0", "hFlatDistRecINELgt0", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistRecHM", "hFlatDistRecHM", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hCentFT0MFine", "hCentFT0MFine", {HistType::kTH1D, {{1000, 0.0f, 100.0f, "FT0M percentile"}}}); + if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessGenMC || + doprocessFlatDistMC) { + // one entry per generated collision, the sample the closure denominator + // lives in: binning.axisFlatTrue has to be set from this one, or the + // numerator and the denominator do not share a class definition + rEventSelection.add("hFlatDistGen", "hFlatDistGen", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistGenINELgt0", "hFlatDistGenINELgt0", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistGenHM", "hFlatDistGenHM", {HistType::kTH1D, {flatFineAxis}}); + // same, restricted to accepted reconstructed collisions: the ratio to + // hFlatDistGen is the event selection bias on the class assignment + rEventSelection.add("hFlatDistGenInRec", "hFlatDistGenInRec", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistGenInRecINELgt0", "hFlatDistGenInRecINELgt0", {HistType::kTH1D, {flatFineAxis}}); + rEventSelection.add("hFlatDistGenInRecHM", "hFlatDistGenInRecHM", {HistType::kTH1D, {flatFineAxis}}); + } + if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessGenMC || + doprocessFlatDistMC) { rEventSelection.add("hTrueFV0amplvsFlat", "TrueFV0MvsFlat", HistType::kTH2D, - {{500, -0.5, +499.5, "True Nch in FV0 region"}, flatAxis}); + {{500, -0.5, +499.5, "True Nch in FV0 region"}, flatTrueAxis}); } - if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade) { + if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessFlatDistMC) { rEventSelection.add("hFlattenicityDistributionMCGen_Rec", "hFlattenicityDistributionMCGen_Rec", - {HistType::kTH1D, {flatAxis}}); + {HistType::kTH1D, {flatTrueAxis}}); rEventSelection.add("hFlattenicity_Corr_Gen_vs_Rec", "hFlattenicity_Corr_Gen_vs_Rec", - {HistType::kTH2D, {flatAxis, flatAxis}}); + {HistType::kTH2D, {flatTrueAxis, flatAxis}}); + // migration of the class assignment, needed pT-differentially to unfold + rEventSelection.add("hFlatGenVsRecFine", "hFlatGenVsRecFine", + {HistType::kTH2D, {flatFineAxis, flatFineAxis}}); + } + + // Charged particles in |eta| < cfgEtaChargedCut: per class is + // the scale factor of Qpp, and the DCAxy templates give the secondary + // contamination the same way the published analysis obtains it. + if (eventClass.fillChargedQA) { + rCharged.add("hNchVsFlat", "hNchVsFlat", {HistType::kTH2D, {nchAxis, flatAxis}}); + rCharged.add("hPtChVsFlat", "hPtChVsFlat", {HistType::kTH2D, {ptChargedAxis, flatAxis}}); + if (doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessFlatDistMC) { + rCharged.add("hPtChVsFlatRecPrim", "hPtChVsFlatRecPrim", + {HistType::kTH2D, {ptChargedAxis, flatAxis}}); + rCharged.add("hPtChVsFlatGenInRec", "hPtChVsFlatGenInRec", + {HistType::kTH2D, {ptChargedAxis, flatAxis}}); + rCharged.add("hNchVsFlatGenInRec", "hNchVsFlatGenInRec", {HistType::kTH2D, {nchAxis, flatAxis}}); + rCharged.add("hDcaXyPtChPrim", "hDcaXyPtChPrim", + {HistType::kTH3D, {dcaXyAxis, ptChargedAxis, flatAxis}}); + rCharged.add("hDcaXyPtChSec", "hDcaXyPtChSec", + {HistType::kTH3D, {dcaXyAxis, ptChargedAxis, flatAxis}}); + } + if (doprocessDataRun3LambdaK0s || doprocessDataRun3Cascade || doprocessFlatDistData) { + rCharged.add("hDcaXyPtCh", "hDcaXyPtCh", + {HistType::kTH3D, {dcaXyAxis, ptChargedAxis, flatAxis}}); + } + if (doprocessGenMC) { + rCharged.add("hPtChVsFlatGen", "hPtChVsFlatGen", {HistType::kTH2D, {ptChargedAxis, flatAxis}}); + rCharged.add("hNchVsFlatGen", "hNchVsFlatGen", {HistType::kTH2D, {nchAxis, flatAxis}}); + rCharged.add("hPtChVsTrueFlatGen", "hPtChVsTrueFlatGen", {HistType::kTH2D, {ptChargedAxis, flatTrueAxis}}); + rCharged.add("hNchVsTrueFlatGen", "hNchVsTrueFlatGen", {HistType::kTH2D, {nchAxis, flatTrueAxis}}); + } } if (doprocessDataRun3LambdaK0s || doprocessRecMCLambdaK0s) { @@ -408,17 +549,17 @@ struct Lambdak0sflattenicity { rKzeroShort.add("hNSigmaNegPionFromK0s", "hNSigmaNegPionFromK0s", {HistType::kTH2D, {{100, -5.f, 5.f, "n#sigma_{TPC}"}, {pTPCAxis}}}); rKzeroShort.add("hMassK0spT", "hMassK0spT", - {HistType::kTH2D, {{k0sMassAxis}, {ptAxis}}}); + {HistType::kTH2D, {{k0sMassAxis}, {ptK0sAxis}}}); rKzeroShort.add("hMassK0spTFlat", "hMassK0spTFlat", - {HistType::kTH3D, {{k0sMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{k0sMassAxis}, {ptK0sAxis}, {flatAxis}}}); rKzeroShort.add("hArmPodoAlphavsQTK0sAfterCut", "hArmPodoAlphavsQTK0sAfterCut", {HistType::kTH2D, {{200, -1, 1, "#alpha"}, {70, 0, 0.35, "Q_{T}"}}}); if (doprocessRecMCLambdaK0s) { rKzeroShort.add("Generated_MCRecoCollCheck_INEL_K0Short", "Generated_MCRecoCollCheck_INEL_K0Short", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptK0sAxis}, {flatAxis}}}); rKzeroShort.add("Generated_MCRecoCollCheck_INELgt0_K0Short", "Generated_MCRecoCollCheck_INELgt0_K0Short", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptK0sAxis}, {flatAxis}}}); } // Lambda reconstruction Mass @@ -443,24 +584,31 @@ struct Lambdak0sflattenicity { rLambda.add("h2DdecayRadiusLambda", "h2DdecayRadiusLambda", {HistType::kTH1D, {decayRadiusAxis}}); rLambda.add("hMassLambdapT", "hMassLambdapT", - {HistType::kTH2D, {{lambdaMassAxis}, {ptAxis}}}); + {HistType::kTH2D, {{lambdaMassAxis}, {ptLambdaAxis}}}); rLambda.add("hMassLambdapTFlat", "hMassLambdapTFlat", - {HistType::kTH3D, {{lambdaMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{lambdaMassAxis}, {ptLambdaAxis}, {flatAxis}}}); if (doprocessRecMCLambdaK0s) { rLambda.add("Generated_MCRecoCollCheck_INEL_Lambda", "Generated_MCRecoCollCheck_INEL_Lambda", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {flatAxis}}}); rLambda.add("Generated_MCRecoCollCheck_INELgt0_Lambda", "Generated_MCRecoCollCheck_INELgt0_Lambda", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {flatAxis}}}); rLambda.add("hMassFeedDownLambdapTFlat", "hMassFeedDownLambdapTFlat", - {HistType::kTH3D, {{lambdaMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{lambdaMassAxis}, {ptLambdaAxis}, {flatAxis}}}); rLambda.add("hFeedDownLambdaPtVsMotherPt", "hFeedDownLambdaPtVsMotherPt", - {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {ptLambdaAxis}}}); + rLambda.add("hFeedDownLambdaMatrix", "hFeedDownLambdaMatrix", + {HistType::kTHnSparseF, {{ptLambdaAxis}, {ptLambdaAxis}, {flatAxis}, {motherAxis}}}); rLambda.add("hFeedDownLambdaMotherPdg", "hFeedDownLambdaMotherPdg", - {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + {HistType::kTH1D, {{motherAxis}}}); 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"); + // generated mothers in the same events and classes as + // hFeedDownLambdaMatrix, so the matrix is a probability per mother that + // folds with the measured Xi. Xi0 is not measurable and only lives here. + rLambda.add("hGenFeedDownMotherPtFlat", "hGenFeedDownMotherPtFlat", + {HistType::kTHnSparseF, {{ptLambdaAxis}, {flatAxis}, {motherAxis}}}); } // AntiLambda reconstruction @@ -488,20 +636,22 @@ struct Lambdak0sflattenicity { rAntiLambda.add("h2DdecayRadiusAntiLambda", "h2DdecayRadiusAntiLambda", {HistType::kTH1D, {decayRadiusAxis}}); rAntiLambda.add("hMassAntiLambdapT", "hMassAntiLambdapT", - {HistType::kTH2D, {{antilambdaMassAxis}, {ptAxis}}}); + {HistType::kTH2D, {{antilambdaMassAxis}, {ptLambdaAxis}}}); rAntiLambda.add("hMassAntiLambdapTFlat", "hMassAntiLambdapTFlat", - {HistType::kTH3D, {{antilambdaMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{antilambdaMassAxis}, {ptLambdaAxis}, {flatAxis}}}); if (doprocessRecMCLambdaK0s) { rAntiLambda.add("Generated_MCRecoCollCheck_INEL_AntiLambda", "Generated_MCRecoCollCheck_INEL_AntiLambda", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {flatAxis}}}); rAntiLambda.add("Generated_MCRecoCollCheck_INELgt0_AntiLambda", "Generated_MCRecoCollCheck_INELgt0_AntiLambda", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {flatAxis}}}); rAntiLambda.add("hMassFeedDownAntiLambdapTFlat", "hMassFeedDownAntiLambdapTFlat", - {HistType::kTH3D, {{antilambdaMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{antilambdaMassAxis}, {ptLambdaAxis}, {flatAxis}}}); rAntiLambda.add("hFeedDownAntiLambdaPtVsMotherPt", "hFeedDownAntiLambdaPtVsMotherPt", - {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + {HistType::kTH2D, {{ptLambdaAxis}, {ptLambdaAxis}}}); + rAntiLambda.add("hFeedDownAntiLambdaMatrix", "hFeedDownAntiLambdaMatrix", + {HistType::kTHnSparseF, {{ptLambdaAxis}, {ptLambdaAxis}, {flatAxis}, {motherAxis}}}); rAntiLambda.add("hFeedDownAntiLambdaMotherPdg", "hFeedDownAntiLambdaMotherPdg", - {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + {HistType::kTH1D, {{motherAxis}}}); 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^{-}"); @@ -510,6 +660,42 @@ struct Lambdak0sflattenicity { rCommonHist.add("hArmPodoAlphavsQT", "hArmPodoAlphavsQT", {HistType::kTH2D, {{200, -1, 1, "#alpha"}, {70, 0, 0.35, "Q_{T}"}}}); + + // DCA of the V0 to the PV: primaries peak at zero, feed-down does not, so + // the data distribution can be fitted with the two MC templates + rLambda.add("hDcaV0ToPVLambda", "hDcaV0ToPVLambda", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + rAntiLambda.add("hDcaV0ToPVAntiLambda", "hDcaV0ToPVAntiLambda", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + if (doprocessRecMCLambdaK0s) { + rLambda.add("hDcaV0ToPVLambdaPrim", "hDcaV0ToPVLambdaPrim", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + rLambda.add("hDcaV0ToPVLambdaSec", "hDcaV0ToPVLambdaSec", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + rAntiLambda.add("hDcaV0ToPVAntiLambdaPrim", "hDcaV0ToPVAntiLambdaPrim", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + rAntiLambda.add("hDcaV0ToPVAntiLambdaSec", "hDcaV0ToPVAntiLambdaSec", + {HistType::kTH3D, {{dcaV0ToPvAxis}, {ptLambdaAxis}, {flatAxis}}}); + + rKzeroShort.add("hPtResK0s", "hPtResK0s", {HistType::kTH2D, {{ptK0sAxis}, {ptResAxis}}}); + rLambda.add("hPtResLambda", "hPtResLambda", {HistType::kTH2D, {{ptLambdaAxis}, {ptResAxis}}}); + rAntiLambda.add("hPtResAntiLambda", "hPtResAntiLambda", {HistType::kTH2D, {{ptLambdaAxis}, {ptResAxis}}}); + + // numerator classified by the measured 1-rho, denominator by the true + // one: the pair is the MC closure + rKzeroShort.add("Generated_MCRecoCollCheck_INELgt0_K0Short_TrueFlat", "Generated_MCRecoCollCheck_INELgt0_K0Short_TrueFlat", + {HistType::kTH2D, {{ptK0sAxis}, {flatTrueAxis}}}); + rLambda.add("Generated_MCRecoCollCheck_INELgt0_Lambda_TrueFlat", "Generated_MCRecoCollCheck_INELgt0_Lambda_TrueFlat", + {HistType::kTH2D, {{ptLambdaAxis}, {flatTrueAxis}}}); + rAntiLambda.add("Generated_MCRecoCollCheck_INELgt0_AntiLambda_TrueFlat", "Generated_MCRecoCollCheck_INELgt0_AntiLambda_TrueFlat", + {HistType::kTH2D, {{ptLambdaAxis}, {flatTrueAxis}}}); + rKzeroShort.add("hMassK0spTTrueFlat", "hMassK0spTTrueFlat", + {HistType::kTH3D, {{k0sMassAxis}, {ptK0sAxis}, {flatTrueAxis}}}); + rLambda.add("hMassLambdapTTrueFlat", "hMassLambdapTTrueFlat", + {HistType::kTH3D, {{lambdaMassAxis}, {ptLambdaAxis}, {flatTrueAxis}}}); + rAntiLambda.add("hMassAntiLambdapTTrueFlat", "hMassAntiLambdapTTrueFlat", + {HistType::kTH3D, {{antilambdaMassAxis}, {ptLambdaAxis}, {flatTrueAxis}}}); + } } if (doprocessRecMCRun3Cascade || doprocessDataRun3Cascade) { @@ -535,24 +721,32 @@ struct Lambdak0sflattenicity { rXi.add("hNSigmaBachPionFromXi", "hNSigmaBachPionFromXi", {HistType::kTH2D, {{100, -5.f, 5.f, "n#sigma_{TPC}"}, {pTPCAxis}}}); rXi.add("hMassXipT", "hMassXipT", - {HistType::kTH2D, {{xiMassAxis}, {ptAxis}}}); + {HistType::kTH2D, {{xiMassAxis}, {ptXiAxis}}}); rXi.add("hMassXipTFlat", "hMassXipTFlat", - {HistType::kTH3D, {{xiMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{xiMassAxis}, {ptXiAxis}, {flatAxis}}}); if (doprocessRecMCRun3Cascade) { rXi.add("Generated_MCRecoCollCheck_INEL_Xi", "Generated_MCRecoCollCheck_INEL_Xi", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptXiAxis}, {flatAxis}}}); rXi.add("Generated_MCRecoCollCheck_INELgt0_Xi", "Generated_MCRecoCollCheck_INELgt0_Xi", - {HistType::kTH2D, {{ptAxis}, {flatAxis}}}); + {HistType::kTH2D, {{ptXiAxis}, {flatAxis}}}); rXi.add("hMassFeedDownXipTFlat", "hMassFeedDownXipTFlat", - {HistType::kTH3D, {{xiMassAxis}, {ptAxis}, {flatAxis}}}); + {HistType::kTH3D, {{xiMassAxis}, {ptXiAxis}, {flatAxis}}}); rXi.add("hFeedDownXiPtVsMotherPt", "hFeedDownXiPtVsMotherPt", - {HistType::kTH2D, {{ptAxis}, {ptAxis}}}); + {HistType::kTH2D, {{ptXiAxis}, {ptXiAxis}}}); + rXi.add("hFeedDownXiMatrix", "hFeedDownXiMatrix", + {HistType::kTHnSparseF, {{ptXiAxis}, {ptXiAxis}, {flatAxis}, {motherAxis}}}); rXi.add("hFeedDownXiMotherPdg", "hFeedDownXiMotherPdg", - {HistType::kTH1D, {{kNFeedDownMothers, -0.5, kNFeedDownMothers - 0.5, "mother"}}}); + {HistType::kTH1D, {{motherAxis}}}); 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"); + + rXi.add("hPtResXi", "hPtResXi", {HistType::kTH2D, {{ptXiAxis}, {ptResAxis}}}); + rXi.add("Generated_MCRecoCollCheck_INELgt0_Xi_TrueFlat", "Generated_MCRecoCollCheck_INELgt0_Xi_TrueFlat", + {HistType::kTH2D, {{ptXiAxis}, {flatTrueAxis}}}); + rXi.add("hMassXipTTrueFlat", "hMassXipTTrueFlat", + {HistType::kTH3D, {{xiMassAxis}, {ptXiAxis}, {flatTrueAxis}}}); } } if (doprocessGenMC) { @@ -563,7 +757,7 @@ struct Lambdak0sflattenicity { {HistType::kTH1D, {vertexZAxis}}); rEventSelection.add("hFlattenicityDistributionMCGen", "hFlattenicityDistributionMCGen", - {HistType::kTH1D, {flatAxis}}); + {HistType::kTH1D, {flatTrueAxis}}); rEventSelection.add("hFlattenicityDistributionRecMCGen", "hFlattenicityDistributionRecMCGen", {HistType::kTH1D, {flatAxis}}); @@ -586,56 +780,89 @@ struct Lambdak0sflattenicity { rEventSelection.get(HIST("hNEventsMCReco"))->GetXaxis()->SetBinLabel(2, "pass ev sel"); rEventSelection.get(HIST("hNEventsMCReco"))->GetXaxis()->SetBinLabel(3, "INELgt0"); rKzeroShort.add("pGen_MCGenRecoColl_INEL_K0Short", "pGen_MCGenRecoColl_INEL_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rKzeroShort.add("Generated_MCRecoColl_INEL_K0Short", "Generated_MCRecoColl_INEL_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rKzeroShort.add("pGen_MCGenColl_INEL_K0Short", "pGen_MCGenColl_INEL_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rKzeroShort.add("pGen_MCGenRecoColl_INELgt0_K0Short", "pGen_MCGenRecoColl_INELgt0_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rKzeroShort.add("Generated_MCRecoColl_INELgt0_K0Short", "Generated_MCRecoColl_INELgt0_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rKzeroShort.add("pGen_MCGenColl_INELgt0_K0Short", "pGen_MCGenColl_INELgt0_K0Short", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptK0sAxis, flatAxis}}); rLambda.add("pGen_MCGenRecoColl_INEL_Lambda", "pGen_MCGenRecoColl_INEL_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rLambda.add("Generated_MCRecoColl_INEL_Lambda", "Generated_MCRecoColl_INEL_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rLambda.add("pGen_MCGenColl_INEL_Lambda", "pGen_MCGenColl_INEL_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rLambda.add("pGen_MCGenRecoColl_INELgt0_Lambda", "pGen_MCGenRecoColl_INELgt0_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rLambda.add("Generated_MCRecoColl_INELgt0_Lambda", "Generated_MCRecoColl_INELgt0_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rLambda.add("pGen_MCGenColl_INELgt0_Lambda", "pGen_MCGenColl_INELgt0_Lambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("pGen_MCGenRecoColl_INEL_AntiLambda", "pGen_MCGenRecoColl_INEL_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("Generated_MCRecoColl_INEL_AntiLambda", "Generated_MCRecoColl_INEL_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("pGen_MCGenColl_INEL_AntiLambda", "pGen_MCGenColl_INEL_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("pGen_MCGenRecoColl_INELgt0_AntiLambda", "pGen_MCGenRecoColl_INELgt0_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("Generated_MCRecoColl_INELgt0_AntiLambda", "Generated_MCRecoColl_INELgt0_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rAntiLambda.add("pGen_MCGenColl_INELgt0_AntiLambda", "pGen_MCGenColl_INELgt0_AntiLambda", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptLambdaAxis, flatAxis}}); rXi.add("pGen_MCGenRecoColl_INEL_Xi", "pGen_MCGenRecoColl_INEL_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); rXi.add("Generated_MCRecoColl_INEL_Xi", "Generated_MCRecoColl_INEL_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); rXi.add("pGen_MCGenColl_INEL_Xi", "pGen_MCGenColl_INEL_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); rXi.add("pGen_MCGenRecoColl_INELgt0_Xi", "pGen_MCGenRecoColl_INELgt0_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); rXi.add("Generated_MCRecoColl_INELgt0_Xi", "Generated_MCRecoColl_INELgt0_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); rXi.add("pGen_MCGenColl_INELgt0_Xi", "pGen_MCGenColl_INELgt0_Xi", - {HistType::kTH2D, {ptAxis, flatAxis}}); + {HistType::kTH2D, {ptXiAxis, flatAxis}}); + + // The same INEL>0 counters against the true 1-rho. Together with the + // reconstructed ones above they give the closure without a second job. + rEventSelection.add("hFlat_RecoColl_MC_INELgt0_TrueFlat", "hFlat_RecoColl_MC_INELgt0_TrueFlat", + {HistType::kTH1D, {flatTrueAxis}}); + rEventSelection.add("hFlat_GenRecoColl_MC_INELgt0_TrueFlat", "hFlat_GenRecoColl_MC_INELgt0_TrueFlat", + {HistType::kTH1D, {flatTrueAxis}}); + rEventSelection.add("hFlat_GenColl_MC_INELgt0_TrueFlat", "hFlat_GenColl_MC_INELgt0_TrueFlat", + {HistType::kTH1D, {flatTrueAxis}}); + rKzeroShort.add("pGen_MCGenRecoColl_INELgt0_K0Short_TrueFlat", "pGen_MCGenRecoColl_INELgt0_K0Short_TrueFlat", + {HistType::kTH2D, {ptK0sAxis, flatTrueAxis}}); + rKzeroShort.add("Generated_MCRecoColl_INELgt0_K0Short_TrueFlat", "Generated_MCRecoColl_INELgt0_K0Short_TrueFlat", + {HistType::kTH2D, {ptK0sAxis, flatTrueAxis}}); + rKzeroShort.add("pGen_MCGenColl_INELgt0_K0Short_TrueFlat", "pGen_MCGenColl_INELgt0_K0Short_TrueFlat", + {HistType::kTH2D, {ptK0sAxis, flatTrueAxis}}); + rLambda.add("pGen_MCGenRecoColl_INELgt0_Lambda_TrueFlat", "pGen_MCGenRecoColl_INELgt0_Lambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rLambda.add("Generated_MCRecoColl_INELgt0_Lambda_TrueFlat", "Generated_MCRecoColl_INELgt0_Lambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rLambda.add("pGen_MCGenColl_INELgt0_Lambda_TrueFlat", "pGen_MCGenColl_INELgt0_Lambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rAntiLambda.add("pGen_MCGenRecoColl_INELgt0_AntiLambda_TrueFlat", "pGen_MCGenRecoColl_INELgt0_AntiLambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rAntiLambda.add("Generated_MCRecoColl_INELgt0_AntiLambda_TrueFlat", "Generated_MCRecoColl_INELgt0_AntiLambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rAntiLambda.add("pGen_MCGenColl_INELgt0_AntiLambda_TrueFlat", "pGen_MCGenColl_INELgt0_AntiLambda_TrueFlat", + {HistType::kTH2D, {ptLambdaAxis, flatTrueAxis}}); + rXi.add("pGen_MCGenRecoColl_INELgt0_Xi_TrueFlat", "pGen_MCGenRecoColl_INELgt0_Xi_TrueFlat", + {HistType::kTH2D, {ptXiAxis, flatTrueAxis}}); + rXi.add("Generated_MCRecoColl_INELgt0_Xi_TrueFlat", "Generated_MCRecoColl_INELgt0_Xi_TrueFlat", + {HistType::kTH2D, {ptXiAxis, flatTrueAxis}}); + rXi.add("pGen_MCGenColl_INELgt0_Xi_TrueFlat", "pGen_MCGenColl_INELgt0_Xi_TrueFlat", + {HistType::kTH2D, {ptXiAxis, flatTrueAxis}}); } if (flatSel.flattenicityQA) { @@ -728,7 +955,33 @@ struct Lambdak0sflattenicity { LOGF(fatal, "Can not run MCGen and Data process functions together. Try one of these at a time"); } - fillFlattenicityQAInGenMC = !(doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade); + // the spectra passes come first, their event counts normalise the yields + if (doprocessRecMCLambdaK0s) { + eventHistOwner = kOwnerRecMCV0; + } else if (doprocessDataRun3LambdaK0s) { + eventHistOwner = kOwnerDataV0; + } else if (doprocessRecMCRun3Cascade) { + eventHistOwner = kOwnerRecMCCasc; + } else if (doprocessDataRun3Cascade) { + eventHistOwner = kOwnerDataCasc; + } else if (doprocessFlatDistMC) { + eventHistOwner = kOwnerFlatDistMC; + } else if (doprocessFlatDistData) { + eventHistOwner = kOwnerFlatDistData; + } + + // tied to the owner so a new process function cannot bring the double fill back + fillFlattenicityQAInGenMC = (eventHistOwner == kOwnerNone); + + // the generated 1-rho covers the FV0 only, the other estimators would + // classify rec and gen with two different observables + if ((doprocessRecMCLambdaK0s || doprocessRecMCRun3Cascade || doprocessGenMC || + doprocessFlatDistMC) && + flatSel.flattenicityforanalysis != kFlatFromFV0) { + LOGF(fatal, + "flattenicityforanalysis!=0 has no generator-level counterpart: " + "estimateFlattenicityFV0MC covers the FV0 acceptance only"); + } // the estimator used for the analysis has to be computed if (flatSel.flattenicityforanalysis == kFlatFromFV0 && !flatSel.isflattenicitywithFV0 && !flatSel.isflattenicitywithFV0FT0C) { @@ -797,53 +1050,6 @@ struct Lambdak0sflattenicity { return flat; } // V0A signal and flatenicity calculation - 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, - 1.06622, 1.15057, 1.0433, 0.83654, 0.847178, 0.890027, 0.920814, - 0.888271, 1.04662, 0.8869, 0.856348, 0.863181, 0.906312, 0.902166, - 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 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 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 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 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 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 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, - 0.985635, 0.982883, 0.981914, 0.964635, 0.967812, 0.95475, - 0.956687, 0.932816, 0.92773, 0.914892, 0.891724, 0.872382}; static constexpr int kNeta5 = 2; // FT0C + FT0A static constexpr std::array kWeigthsEta5 = {0.0490638, 0.010958415}; @@ -989,12 +1195,213 @@ struct Lambdak0sflattenicity { return true; } + // Daughter track acceptance and quality. The quality cuts are off by default + // and exist so a systematic variation can move them from the JSON. + template + bool isSelectedDaughterTrack(TTrack const& track, float minCrossedRows) + { + if (std::abs(track.eta()) > trkPid.cfgTrkEtaCut || track.pt() < trkPid.cfgTrkLowPtCut) { + return false; + } + if (track.tpcNClsCrossedRows() < minCrossedRows) { + return false; + } + if (trkPid.minCrossedRowsOverFindable > 0.f && + track.tpcCrossedRowsOverFindableCls() < trkPid.minCrossedRowsOverFindable) { + return false; + } + if (trkPid.maxTpcChi2NCl > 0.f && track.tpcChi2NCl() > trkPid.maxTpcChi2NCl) { + return false; + } + if (trkPid.maxItsChi2NCl > 0.f && track.itsChi2NCl() > trkPid.maxItsChi2NCl) { + return false; + } + if (trkPid.minItsNCls > 0 && static_cast(track.itsNCls()) < trkPid.minItsNCls) { + return false; + } + return true; + } + + // 1-rho on the fine axis for the percentile boundaries, inclusive, INEL>0 and + // in the top FT0M class + void fillFlatDistRec(float flattenicity, float centFT0M, bool isInelGt0) + { + rEventSelection.fill(HIST("hFlatDistRec"), flattenicity); + rEventSelection.fill(HIST("hCentFT0MFine"), centFT0M); + if (isInelGt0) { + rEventSelection.fill(HIST("hFlatDistRecINELgt0"), flattenicity); + if (centFT0M < eventClass.cfgCentMaxHM) { + rEventSelection.fill(HIST("hFlatDistRecHM"), flattenicity); + } + } + } + + // processGenMC only. A negative centFT0M means no accepted reconstructed + // counterpart, so the collision has no FT0M class and stays out of the HM one. + void fillFlatDistGen(float flattenicityGen, float centFT0M, bool isInelGt0) + { + rEventSelection.fill(HIST("hFlatDistGen"), flattenicityGen); + if (isInelGt0) { + rEventSelection.fill(HIST("hFlatDistGenINELgt0"), flattenicityGen); + if (centFT0M >= 0.f && centFT0M < eventClass.cfgCentMaxHM) { + rEventSelection.fill(HIST("hFlatDistGenHM"), flattenicityGen); + } + } + } + + // the same from the rec-level process functions, diagnostic only + void fillFlatDistGenInRec(float flattenicityGen, float centFT0M, bool isInelGt0) + { + rEventSelection.fill(HIST("hFlatDistGenInRec"), flattenicityGen); + if (isInelGt0) { + rEventSelection.fill(HIST("hFlatDistGenInRecINELgt0"), flattenicityGen); + if (centFT0M >= 0.f && centFT0M < eventClass.cfgCentMaxHM) { + rEventSelection.fill(HIST("hFlatDistGenInRecHM"), flattenicityGen); + } + } + } + + // Charged tracks in the tracking acceptance: per class is the + // scale factor of Qpp, and the DCAxy distribution carries the secondary + // contamination. + template + int fillChargedRec(TTracks const& tracks, float flattenicity) + { + if (!eventClass.fillChargedQA) { + return 0; + } + int nch = 0; + for (const auto& track : tracks) { + if (!track.isGlobalTrack() || std::abs(track.eta()) > eventClass.cfgEtaChargedCut) { + continue; + } + nch++; + rCharged.fill(HIST("hPtChVsFlat"), track.pt(), flattenicity); + if constexpr (isMC) { + if (!track.has_mcParticle()) { + continue; + } + const auto& mcParticle = track.mcParticle(); + if (mcParticle.isPhysicalPrimary()) { + rCharged.fill(HIST("hPtChVsFlatRecPrim"), track.pt(), flattenicity); + rCharged.fill(HIST("hDcaXyPtChPrim"), track.dcaXY(), track.pt(), flattenicity); + } else { + rCharged.fill(HIST("hDcaXyPtChSec"), track.dcaXY(), track.pt(), flattenicity); + } + } else { + rCharged.fill(HIST("hDcaXyPtCh"), track.dcaXY(), track.pt(), flattenicity); + } + } + rCharged.fill(HIST("hNchVsFlat"), nch, flattenicity); + return nch; + } + + template + int fillChargedGen(TMcParticles const& mcParticles, float flattenicity, bool trueFlat) + { + if (!eventClass.fillChargedQA) { + return 0; + } + int nch = 0; + for (const auto& mcParticle : mcParticles) { + if (!mcParticle.isPhysicalPrimary() || std::abs(mcParticle.eta()) > eventClass.cfgEtaChargedCut) { + continue; + } + auto pdgParticle = pdg->GetParticle(mcParticle.pdgCode()); + if (!(pdgParticle && std::abs(pdgParticle->Charge()) > kMinCharge)) { + continue; + } + nch++; + if (trueFlat) { + rCharged.fill(HIST("hPtChVsTrueFlatGen"), mcParticle.pt(), flattenicity); + } else { + rCharged.fill(HIST("hPtChVsFlatGen"), mcParticle.pt(), flattenicity); + } + } + if (trueFlat) { + rCharged.fill(HIST("hNchVsTrueFlatGen"), nch, flattenicity); + } else { + rCharged.fill(HIST("hNchVsFlatGen"), nch, flattenicity); + } + return nch; + } + + // ================= Calibration objects from CCDB ==================== // + // A missing or malformed object falls back to unity, so a run without + // calibration is left uncorrected instead of being equalized with somebody + // else's constants. + + std::vector fetchGainEq(const std::string& path, int run, std::size_t nChannels) + { + const auto* obj = ccdb->getForRun>(path, run); + if (!obj || obj->size() != nChannels) { + LOGF(warning, "No gain equalization of size %zu at %s for run %d, using unity", + nChannels, path.c_str(), run); + return std::vector(nChannels, 1.f); + } + return *obj; + } + + TProfile2D* fetchVtxEq(const std::string& path, int run) + { + auto* obj = ccdb->getForRun(path, run); + if (!obj) { + LOGF(warning, "No z-vertex equalization at %s for run %d, lattice left uncorrected", + path.c_str(), run); + } + return obj; + } + + // called per collision, does nothing unless a correction is switched on + template + void initCcdb(TBC const& bc) + { + const int run = bc.runNumber(); + if (run == mRunNumber) { + return; + } + mRunNumber = run; + if (flatSel.applyCalibCh) { + gainFV0 = fetchGainEq(ccdbConf.gainEqPath.value + "/FV0", run, kNCells); + gainFT0A = fetchGainEq(ccdbConf.gainEqPath.value + "/FT0A", run, kNCellsT0A); + gainFT0C = fetchGainEq(ccdbConf.gainEqPath.value + "/FT0C", run, kNCellsT0C); + } + if (flatSel.applyCalibVtx) { + vtxEqFV0 = fetchVtxEq(ccdbConf.vtxEqPath.value + "/FV0", run); + vtxEqFT0A = fetchVtxEq(ccdbConf.vtxEqPath.value + "/FT0A", run); + vtxEqFT0C = fetchVtxEq(ccdbConf.vtxEqPath.value + "/FT0C", run); + } + } + + // the gain is divided out, the convention of flattenictyPikp + static float gainEqFactor(std::vector const& gain, std::size_t channel) + { + if (channel >= gain.size() || !(gain[channel] > 0.f)) { + return 1.f; + } + return gain[channel]; + } + + // 1 for a missing map or an empty bin, so a bad bin never zeroes a cell + static float vtxEqFactor(TProfile2D const* map, int cell, float vtxZ) + { + if (!map) { + return 1.f; + } + const float factor = map->GetBinContent(map->GetXaxis()->FindBin(cell), + map->GetYaxis()->FindBin(vtxZ)); + return (factor > 0.f) ? factor : 1.f; + } + // ============== Flattenicity estimation begins ===================== // // 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 = flatSel.flattenicityQA && fillQA; + if (flatSel.applyCalibCh || flatSel.applyCalibVtx) { + initCcdb(collision.template bc_as>()); + } std::array ampl5 = {0, 0}; std::array ampl6 = {0, 0}; @@ -1002,6 +1409,8 @@ struct Lambdak0sflattenicity { float sumAmpFV0 = 0; float sumAmpFV01to4Ch = 0; + // gain-equalized but not yet vertex-equalized, the input to the z-vertex QA + float sumAmpFV0BeforeVtx = 0; ampchannel.fill(0.0); ampchannelBefore.fill(0.0); @@ -1026,7 +1435,11 @@ struct Lambdak0sflattenicity { } ampchannelBefore[channelv0phi] = amplCh; if (flatSel.applyCalibCh) { - amplCh *= kCalib[channelv0phi]; + amplCh /= gainEqFactor(gainFV0, channelv0); + } + sumAmpFV0BeforeVtx += amplCh; + if (flatSel.applyCalibVtx) { + amplCh *= vtxEqFactor(vtxEqFV0, channelv0phi, vtxZ); } sumAmpFV0 += amplCh; @@ -1045,13 +1458,7 @@ struct Lambdak0sflattenicity { } if (flattenicityQAhere) { - rFlattenicity.fill(HIST("hAmpV0vsVtxBeforeCalibration"), vtxZ, sumAmpFV0); - } - if (flatSel.applyCalibVtx) { - sumAmpFV0 *= gVtx[0].Eval(vtxZ); - sumAmpFV01to4Ch *= gVtx[0].Eval(vtxZ); - } - if (flattenicityQAhere) { + rFlattenicity.fill(HIST("hAmpV0vsVtxBeforeCalibration"), vtxZ, sumAmpFV0BeforeVtx); rFlattenicity.fill(HIST("hAmpV0vsVtx"), vtxZ, sumAmpFV0); } } @@ -1077,6 +1484,8 @@ struct Lambdak0sflattenicity { // FT0 float sumAmpFT0A = 0.f; float sumAmpFT0C = 0.f; + float sumAmpFT0ABeforeVtx = 0.f; + float sumAmpFT0CBeforeVtx = 0.f; rhoLatticeT0A.fill(0); rhoLatticeT0C.fill(0); @@ -1089,20 +1498,26 @@ struct Lambdak0sflattenicity { float amplitude = ft0.amplitudeA()[i_a]; uint8_t channel = ft0.channelA()[i_a]; int sector = getT0ASector(channel); + float amplitudeBeforeVtx = amplitude; if (sector >= 0 && sector < kNCellsT0A) { if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpT0AVsChBeforeCalibration"), sector, amplitude); } if (flatSel.applyCalibCh) { - amplitude *= kCalibT0A[sector]; + amplitude /= gainEqFactor(gainFT0A, sector); } if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpT0AVsCh"), sector, amplitude); } + amplitudeBeforeVtx = amplitude; + if (flatSel.applyCalibVtx) { + amplitude *= vtxEqFactor(vtxEqFT0A, sector, vtxZ); + } rhoLatticeT0A[sector] += amplitude; } sumAmpFT0A += amplitude; + sumAmpFT0ABeforeVtx += amplitudeBeforeVtx; if (flattenicityQAhere) { rFlattenicity.fill(HIST("hFT0A"), amplitude); } @@ -1113,35 +1528,35 @@ struct Lambdak0sflattenicity { float amplitude = ft0.amplitudeC()[i_c]; uint8_t channel = ft0.channelC()[i_c]; int sector = getT0CSector(channel); + float amplitudeBeforeVtx = amplitude; if (sector >= 0 && sector < kNCellsT0C) { if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpT0CVsChBeforeCalibration"), sector, amplitude); } if (flatSel.applyCalibCh) { - amplitude *= kCalibT0C[sector]; + amplitude /= gainEqFactor(gainFT0C, sector); } if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpT0CVsCh"), sector, amplitude); } + amplitudeBeforeVtx = amplitude; + if (flatSel.applyCalibVtx) { + amplitude *= vtxEqFactor(vtxEqFT0C, sector, vtxZ); + } rhoLatticeT0C[sector] += amplitude; } sumAmpFT0C += amplitude; + sumAmpFT0CBeforeVtx += amplitudeBeforeVtx; if (flattenicityQAhere) { rFlattenicity.fill(HIST("hFT0C"), amplitude); } } if (flattenicityQAhere) { rFlattenicity.fill(HIST("hAmpT0AvsVtxBeforeCalibration"), vtxZ, - sumAmpFT0A); + sumAmpFT0ABeforeVtx); rFlattenicity.fill(HIST("hAmpT0CvsVtxBeforeCalibration"), vtxZ, - sumAmpFT0C); - } - if (flatSel.applyCalibVtx) { - sumAmpFT0A *= gVtx[1].Eval(vtxZ); - sumAmpFT0C *= gVtx[2].Eval(vtxZ); - } - if (flattenicityQAhere) { + sumAmpFT0CBeforeVtx); rFlattenicity.fill(HIST("hAmpT0AvsVtx"), vtxZ, sumAmpFT0A); rFlattenicity.fill(HIST("hAmpT0CvsVtx"), vtxZ, sumAmpFT0C); } @@ -1289,6 +1704,21 @@ struct Lambdak0sflattenicity { static constexpr int kFdOther = 3; static constexpr int kNFeedDownMothers = 4; + // -1 for a species that does not feed the measured hadrons + static int feedDownSpecies(int pdgCode) + { + switch (std::abs(pdgCode)) { + case PDG_t::kXiMinus: + return kFdXiMinus; + case o2::constants::physics::Pdg::kXi0: + return kFdXiZero; + case PDG_t::kOmegaMinus: + return kFdOmegaMinus; + default: + return -1; + } + } + template int getFeedDownMother(TMcParticle const& mcParticle, float& motherPt) { @@ -1298,29 +1728,26 @@ struct Lambdak0sflattenicity { } 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; + const int species = feedDownSpecies(mother.pdgCode()); + return (species >= 0) ? species : kFdOther; } return kFdOther; } + // fillQA=false skips hTrueFV0amplvsFlat, for callers that would fill it twice template - float estimateFlattenicityFV0MC(McParticles const& mcParticles) + float estimateFlattenicityFV0MC(McParticles const& mcParticles, bool fillQA = true) { rhoLatticeFV0AMC.fill(0); int multFV0 = 0; for (const auto& mcParticle : mcParticles) { - if (!(mcParticle.isPhysicalPrimary() && mcParticle.pt() > 0)) { + // the measured 1-rho also sees secondaries; dropping the primary + // requirement is the handle on the leading MC non-closure source + if (eventClass.genFlatPrimariesOnly && !mcParticle.isPhysicalPrimary()) { + continue; + } + if (!(mcParticle.pt() > 0)) { continue; } @@ -1345,7 +1772,12 @@ struct Lambdak0sflattenicity { 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); + // yield per cell with the outer ring halved, as the amplitude is in + // estimateFlattenicity; the alternative normalises to the cell area + rhoLatticeFV0AMC[isegment] += + flatSel.genFlatDetectorLikeNorm + ? ((ieta == kOuterFV0RingIndex) ? 0.5f : 1.0f) + : 1.0 / std::abs(dphi * kDetaFV0); multFV0++; } isegment++; @@ -1355,7 +1787,9 @@ struct Lambdak0sflattenicity { const float flattenicity = 1.0 - getFlatenicity({rhoLatticeFV0AMC.data(), rhoLatticeFV0AMC.size()}); - rEventSelection.fill(HIST("hTrueFV0amplvsFlat"), multFV0, flattenicity); + if (fillQA) { + rEventSelection.fill(HIST("hTrueFV0amplvsFlat"), multFV0, flattenicity); + } return flattenicity; } // ====================== Flattenicity estimation ends ===================== @@ -1371,7 +1805,7 @@ struct Lambdak0sflattenicity { (nabs(aod::track::eta) < trkPid.cfgTrkEtaCut && aod::track::pt > trkPid.cfgTrkLowPtCut); using TrackCandidates = soa::Filtered< - soa::Join>; void processDataRun3LambdaK0s( @@ -1381,9 +1815,10 @@ struct Lambdak0sflattenicity { soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/) { + const bool own = (eventHistOwner == kOwnerDataV0); if (evSel.applyEvSel && - !(isEventSelected(collision))) { // Checking if the event passes the - // selection criteria + !(isEventSelected(collision, own))) { // Checking if the event passes the + // selection criteria return; } @@ -1395,26 +1830,23 @@ struct Lambdak0sflattenicity { if (flattenicity < 0.f) { return; } - rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + if (own) { + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); - rEventSelection.fill(HIST("hVertexZ"), vtxZ); - rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); - rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); - rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + rEventSelection.fill(HIST("hVertexZ"), vtxZ); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracks, flattenicity); + } for (const auto& v0 : V0s) { const auto& posDaughterTrack = v0.posTrack_as(); const auto& negDaughterTrack = v0.negTrack_as(); - 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() < v0Sel.v0settingNTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows) { + if (!isSelectedDaughterTrack(posDaughterTrack, v0Sel.v0settingNTPCcrossedRows) || + !isSelectedDaughterTrack(negDaughterTrack, v0Sel.v0settingNTPCcrossedRows)) { continue; } float massK0s = v0.mK0Short(); @@ -1493,6 +1925,7 @@ struct Lambdak0sflattenicity { rLambda.fill(HIST("h2DdecayRadiusLambda"), v0.v0radius()); rLambda.fill(HIST("hMassLambdapT"), massLambda, v0.pt()); rLambda.fill(HIST("hMassLambdapTFlat"), massLambda, v0.pt(), flattenicity); + rLambda.fill(HIST("hDcaV0ToPVLambda"), v0.dcav0topv(), 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) { @@ -1524,6 +1957,7 @@ struct Lambdak0sflattenicity { rAntiLambda.fill(HIST("hMassAntiLambdapT"), massAntiLambda, v0.pt()); rAntiLambda.fill(HIST("hMassAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); + rAntiLambda.fill(HIST("hDcaV0ToPVAntiLambda"), v0.dcav0topv(), 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) { @@ -1539,7 +1973,7 @@ struct Lambdak0sflattenicity { } using TrackCandidatesMC = - soa::Filtered>; @@ -1556,10 +1990,11 @@ struct Lambdak0sflattenicity { soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/, aod::McParticles const& mcParticles) { + const bool own = (eventHistOwner == kOwnerRecMCV0); for (const auto& collision : collisions) { if (evSel.applyEvSel && - !(isEventSelected(collision))) { // Checking if the event passes the - // selection criteria + !(isEventSelected(collision, own))) { // Checking if the event passes the + // selection criteria continue; } @@ -1576,30 +2011,36 @@ struct Lambdak0sflattenicity { if (flattenicity < 0.f) { continue; } - rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + if (own) { + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); - rEventSelection.fill(HIST("hVertexZ"), vtxZ); - rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); - rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); - rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + rEventSelection.fill(HIST("hVertexZ"), vtxZ); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracksThisCollision, flattenicity); + } auto v0sThisCollision = V0s.sliceBy(perCol, collision.globalIndex()); const auto& mcCollision = collision.mcCollision_as(); + const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache1); + const float flattenicityMCGen = estimateFlattenicityFV0MC(particlesInCollision, own); + if (own) { + rEventSelection.fill(HIST("hFlattenicityDistributionMCGen_Rec"), flattenicityMCGen); + rEventSelection.fill(HIST("hFlattenicity_Corr_Gen_vs_Rec"), flattenicityMCGen, flattenicity); + rEventSelection.fill(HIST("hFlatGenVsRecFine"), flattenicityMCGen, flattenicity); + fillFlatDistGenInRec(flattenicityMCGen, collision.centFT0M(), collision.isInelGt0()); + } + for (const auto& v0 : v0sThisCollision) { const auto& posDaughterTrack = v0.posTrack_as(); const auto& negDaughterTrack = v0.negTrack_as(); - 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() < v0Sel.v0settingNTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < v0Sel.v0settingNTPCcrossedRows) { + if (!isSelectedDaughterTrack(posDaughterTrack, v0Sel.v0settingNTPCcrossedRows) || + !isSelectedDaughterTrack(negDaughterTrack, v0Sel.v0settingNTPCcrossedRows)) { continue; } @@ -1658,6 +2099,9 @@ struct Lambdak0sflattenicity { rKzeroShort.fill(HIST("h2DdecayRadiusK0s"), v0.v0radius()); rKzeroShort.fill(HIST("hMassK0spT"), massK0s, v0.pt()); rKzeroShort.fill(HIST("hMassK0spTFlat"), massK0s, v0.pt(), flattenicity); + rKzeroShort.fill(HIST("hMassK0spTTrueFlat"), massK0s, v0.pt(), flattenicityMCGen); + rKzeroShort.fill(HIST("hPtResK0s"), v0mcParticle.pt(), + (v0.pt() - v0mcParticle.pt()) / v0mcParticle.pt()); rKzeroShort.fill(HIST("hArmPodoAlphavsQTK0sAfterCut"), alpha, qtarm); // Filling the PID of the V0 daughters in the region of the K0s peak @@ -1686,7 +2130,11 @@ struct Lambdak0sflattenicity { const int motherIndex = getFeedDownMother(v0mcParticle, motherPt); rLambda.fill(HIST("hMassFeedDownLambdapTFlat"), massLambda, v0.pt(), flattenicity); rLambda.fill(HIST("hFeedDownLambdaPtVsMotherPt"), v0.pt(), motherPt); + rLambda.fill(HIST("hFeedDownLambdaMatrix"), v0.pt(), motherPt, flattenicity, motherIndex); rLambda.fill(HIST("hFeedDownLambdaMotherPdg"), motherIndex); + rLambda.fill(HIST("hDcaV0ToPVLambdaSec"), v0.dcav0topv(), v0.pt(), flattenicity); + } else { + rLambda.fill(HIST("hDcaV0ToPVLambdaPrim"), v0.dcav0topv(), v0.pt(), flattenicity); } if (keepForEfficiency) { @@ -1698,6 +2146,10 @@ struct Lambdak0sflattenicity { rLambda.fill(HIST("h2DdecayRadiusLambda"), v0.v0radius()); rLambda.fill(HIST("hMassLambdapT"), massLambda, v0.pt()); rLambda.fill(HIST("hMassLambdapTFlat"), massLambda, v0.pt(), flattenicity); + rLambda.fill(HIST("hMassLambdapTTrueFlat"), massLambda, v0.pt(), flattenicityMCGen); + rLambda.fill(HIST("hDcaV0ToPVLambda"), v0.dcav0topv(), v0.pt(), flattenicity); + rLambda.fill(HIST("hPtResLambda"), v0mcParticle.pt(), + (v0.pt() - v0mcParticle.pt()) / v0mcParticle.pt()); // Filling the PID of the V0 daughters in the region of the Lambda peak if (std::abs(massLambda - o2::constants::physics::MassLambda0) < trkPid.pidQAWindowLambda) { @@ -1726,7 +2178,11 @@ struct Lambdak0sflattenicity { const int motherIndex = getFeedDownMother(v0mcParticle, motherPt); rAntiLambda.fill(HIST("hMassFeedDownAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); rAntiLambda.fill(HIST("hFeedDownAntiLambdaPtVsMotherPt"), v0.pt(), motherPt); + rAntiLambda.fill(HIST("hFeedDownAntiLambdaMatrix"), v0.pt(), motherPt, flattenicity, motherIndex); rAntiLambda.fill(HIST("hFeedDownAntiLambdaMotherPdg"), motherIndex); + rAntiLambda.fill(HIST("hDcaV0ToPVAntiLambdaSec"), v0.dcav0topv(), v0.pt(), flattenicity); + } else { + rAntiLambda.fill(HIST("hDcaV0ToPVAntiLambdaPrim"), v0.dcav0topv(), v0.pt(), flattenicity); } if (keepForEfficiency) { @@ -1739,6 +2195,10 @@ struct Lambdak0sflattenicity { rAntiLambda.fill(HIST("h2DdecayRadiusAntiLambda"), v0.v0radius()); rAntiLambda.fill(HIST("hMassAntiLambdapT"), massAntiLambda, v0.pt()); rAntiLambda.fill(HIST("hMassAntiLambdapTFlat"), massAntiLambda, v0.pt(), flattenicity); + rAntiLambda.fill(HIST("hMassAntiLambdapTTrueFlat"), massAntiLambda, v0.pt(), flattenicityMCGen); + rAntiLambda.fill(HIST("hDcaV0ToPVAntiLambda"), v0.dcav0topv(), v0.pt(), flattenicity); + rAntiLambda.fill(HIST("hPtResAntiLambda"), v0mcParticle.pt(), + (v0.pt() - v0mcParticle.pt()) / v0mcParticle.pt()); // Filling the PID of the V0 daughters in the region of the AntiLambda // peak @@ -1754,69 +2214,92 @@ struct Lambdak0sflattenicity { } } - const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache1); - float flattenicityMCGen = estimateFlattenicityFV0MC(particlesInCollision); - rEventSelection.fill(HIST("hFlattenicityDistributionMCGen_Rec"), flattenicityMCGen); - rEventSelection.fill(HIST("hFlattenicity_Corr_Gen_vs_Rec"), flattenicityMCGen, flattenicity); - const bool isInelGt0Rec = collision.isInelGt0(); + int nchGenInRec = 0; for (const auto& mcParticle : particlesInCollision) { if (!mcParticle.isPhysicalPrimary()) { continue; } - if (std::abs(mcParticle.y()) > kMcRapidityWindow) { + if (eventClass.fillChargedQA && std::abs(mcParticle.eta()) <= eventClass.cfgEtaChargedCut) { + auto pdgParticle = pdg->GetParticle(mcParticle.pdgCode()); + if (pdgParticle && std::abs(pdgParticle->Charge()) > kMinCharge) { + nchGenInRec++; + rCharged.fill(HIST("hPtChVsFlatGenInRec"), mcParticle.pt(), flattenicity); + } + } + + // normalisation of hFeedDownLambdaMatrix, before the Lambda rapidity cut + if (std::abs(mcParticle.y()) <= eventClass.cfgFeedDownMotherRapidity) { + const int motherSpecies = feedDownSpecies(mcParticle.pdgCode()); + if (motherSpecies >= 0) { + rLambda.fill(HIST("hGenFeedDownMotherPtFlat"), mcParticle.pt(), flattenicity, motherSpecies); + } + } + + if (std::abs(mcParticle.y()) > v0Sel.v0settingRapidity) { continue; } if (mcParticle.pdgCode() == PDG_t::kK0Short) { - rKzeroShort.fill(HIST("Generated_MCRecoCollCheck_INEL_K0Short"), mcParticle.pt(), flattenicity); // K0s + rKzeroShort.fill(HIST("Generated_MCRecoCollCheck_INEL_K0Short"), mcParticle.pt(), flattenicity); if (isInelGt0Rec) { - rKzeroShort.fill(HIST("Generated_MCRecoCollCheck_INELgt0_K0Short"), mcParticle.pt(), flattenicity); // K0s + rKzeroShort.fill(HIST("Generated_MCRecoCollCheck_INELgt0_K0Short"), mcParticle.pt(), flattenicity); + rKzeroShort.fill(HIST("Generated_MCRecoCollCheck_INELgt0_K0Short_TrueFlat"), mcParticle.pt(), flattenicityMCGen); } } if (mcParticle.pdgCode() == PDG_t::kLambda0) { - rLambda.fill(HIST("Generated_MCRecoCollCheck_INEL_Lambda"), mcParticle.pt(), flattenicity); // Lambda + rLambda.fill(HIST("Generated_MCRecoCollCheck_INEL_Lambda"), mcParticle.pt(), flattenicity); if (isInelGt0Rec) { - rLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Lambda"), mcParticle.pt(), flattenicity); // Lambda + rLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Lambda"), mcParticle.pt(), flattenicity); + rLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Lambda_TrueFlat"), mcParticle.pt(), flattenicityMCGen); } } if (mcParticle.pdgCode() == PDG_t::kLambda0Bar) { - rAntiLambda.fill(HIST("Generated_MCRecoCollCheck_INEL_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + rAntiLambda.fill(HIST("Generated_MCRecoCollCheck_INEL_AntiLambda"), mcParticle.pt(), flattenicity); if (isInelGt0Rec) { - rAntiLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + rAntiLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); + rAntiLambda.fill(HIST("Generated_MCRecoCollCheck_INELgt0_AntiLambda_TrueFlat"), mcParticle.pt(), flattenicityMCGen); } } } + if (eventClass.fillChargedQA) { + rCharged.fill(HIST("hNchVsFlatGenInRec"), nchGenInRec, flattenicity); + } } } // 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, + soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/, o2::aod::McParticles const& mcParticles) { - // without a reconstructed counterpart the sentinel is kept, it falls in the - // underflow so the loss counters stay complete - float flattenicity = kInvalidFlattenicity; - if (flatSel.flattenicityforLossCorrRec) { - for (const auto& collision : collisions) { - if (evSel.applyEvSel && !isEventSelected(collision, false)) { - continue; - } - auto tracksThisCollision = tracks.sliceBy(perColTracksMC, collision.globalIndex()); - flattenicity = estimateFlattenicity(collision, tracksThisCollision, fillFlattenicityQAInGenMC); - if (flattenicity >= 0.f) { - break; - } + // Both estimates are always computed: the reconstructed one classifies the + // numerator of the closure, the true one its denominator. Without a + // reconstructed counterpart the sentinel is kept, it falls in the underflow + // so the loss counters stay complete. + float flattenicityRec = kInvalidFlattenicity; + float centFT0M = -1.f; + for (const auto& collision : collisions) { + if (evSel.applyEvSel && !isEventSelected(collision, false)) { + continue; + } + auto tracksThisCollision = tracks.sliceBy(perColTracksMC, collision.globalIndex()); + flattenicityRec = estimateFlattenicity(collision, tracksThisCollision, fillFlattenicityQAInGenMC); + if (flattenicityRec >= 0.f) { + centFT0M = collision.centFT0M(); + break; } - rEventSelection.fill(HIST("hFlattenicityDistributionRecMCGen"), flattenicity); - } else { - flattenicity = estimateFlattenicityFV0MC(mcParticles); - rEventSelection.fill(HIST("hFlattenicityDistributionMCGen"), flattenicity); } + const float flattenicityTrue = estimateFlattenicityFV0MC(mcParticles, fillFlattenicityQAInGenMC); + rEventSelection.fill(HIST("hFlattenicityDistributionRecMCGen"), flattenicityRec); + rEventSelection.fill(HIST("hFlattenicityDistributionMCGen"), flattenicityTrue); + + // which of the two drives the loss corrections written without a suffix + const float flattenicity = flatSel.flattenicityforLossCorrRec ? flattenicityRec : flattenicityTrue; //==================================== //===== Event Loss Denominator ======= @@ -1828,16 +2311,29 @@ struct Lambdak0sflattenicity { return; } rEventSelection.fill(HIST("hNEventsMCGen"), 1.5); + + // The FT0M class has no generator-level counterpart, so with a window + // requested the collision only belongs to it through an accepted + // reconstructed one, which is what leaves centFT0M non-negative above. + // Otherwise the generated denominators stay MB while the numerator is HM. + if (eventClass.applyCentSel && centFT0M < 0.f) { + return; + } + rEventSelection.fill(HIST("hFlat_GenColl_MC"), flattenicity); - bool isINELgt0true = false; + const bool isINELgt0true = pwglf::isINELgtNmc(mcParticles, 0, pdg); - if (pwglf::isINELgtNmc(mcParticles, 0, pdg)) { - isINELgt0true = true; + if (isINELgt0true) { rEventSelection.fill(HIST("hNEventsMCGen"), 2.5); rEventSelection.fill(HIST("hFlat_GenColl_MC_INELgt0"), flattenicity); + rEventSelection.fill(HIST("hFlat_GenColl_MC_INELgt0_TrueFlat"), flattenicityTrue); } + fillFlatDistGen(flattenicityTrue, centFT0M, isINELgt0true); + fillChargedGen(mcParticles, flattenicity, false); + fillChargedGen(mcParticles, flattenicityTrue, true); + //===================================== //===== Signal Loss Denominator ======= //===================================== @@ -1847,32 +2343,35 @@ struct Lambdak0sflattenicity { if (!mcParticle.isPhysicalPrimary()) { continue; } - if (std::abs(mcParticle.y()) > kMcRapidityWindow) { - continue; - } + const bool inV0Rapidity = std::abs(mcParticle.y()) <= v0Sel.v0settingRapidity; + const bool inCascRapidity = std::abs(mcParticle.y()) <= cascSel.cascsettingRapidity; - if (mcParticle.pdgCode() == PDG_t::kK0Short) { - rKzeroShort.fill(HIST("pGen_MCGenColl_INEL_K0Short"), mcParticle.pt(), flattenicity); // K0s + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kK0Short) { + rKzeroShort.fill(HIST("pGen_MCGenColl_INEL_K0Short"), mcParticle.pt(), flattenicity); if (isINELgt0true) { - rKzeroShort.fill(HIST("pGen_MCGenColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); // K0s + rKzeroShort.fill(HIST("pGen_MCGenColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); + rKzeroShort.fill(HIST("pGen_MCGenColl_INELgt0_K0Short_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0) { - rLambda.fill(HIST("pGen_MCGenColl_INEL_Lambda"), mcParticle.pt(), flattenicity); // Lambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0) { + rLambda.fill(HIST("pGen_MCGenColl_INEL_Lambda"), mcParticle.pt(), flattenicity); if (isINELgt0true) { - rLambda.fill(HIST("pGen_MCGenColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); // Lambda + rLambda.fill(HIST("pGen_MCGenColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); + rLambda.fill(HIST("pGen_MCGenColl_INELgt0_Lambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0Bar) { - rAntiLambda.fill(HIST("pGen_MCGenColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0Bar) { + rAntiLambda.fill(HIST("pGen_MCGenColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); if (isINELgt0true) { - rAntiLambda.fill(HIST("pGen_MCGenColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + rAntiLambda.fill(HIST("pGen_MCGenColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); + rAntiLambda.fill(HIST("pGen_MCGenColl_INELgt0_AntiLambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { - rXi.fill(HIST("pGen_MCGenColl_INEL_Xi"), mcParticle.pt(), flattenicity); // Xi + if (inCascRapidity && std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { + rXi.fill(HIST("pGen_MCGenColl_INEL_Xi"), mcParticle.pt(), flattenicity); if (isINELgt0true) { - rXi.fill(HIST("pGen_MCGenColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); // Xi + rXi.fill(HIST("pGen_MCGenColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); + rXi.fill(HIST("pGen_MCGenColl_INELgt0_Xi_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } } @@ -1902,6 +2401,7 @@ struct Lambdak0sflattenicity { if (collision.isInelGt0() && isINELgt0true) { rEventSelection.fill(HIST("hNEventsMCReco"), 2.5); rEventSelection.fill(HIST("hFlat_RecoColl_MC_INELgt0"), flattenicity); + rEventSelection.fill(HIST("hFlat_RecoColl_MC_INELgt0_TrueFlat"), flattenicityTrue); recoCollIndexINELgt0++; } @@ -1915,33 +2415,35 @@ struct Lambdak0sflattenicity { if (!mcParticle.isPhysicalPrimary()) { continue; } + const bool inV0Rapidity = std::abs(mcParticle.y()) <= v0Sel.v0settingRapidity; + const bool inCascRapidity = std::abs(mcParticle.y()) <= cascSel.cascsettingRapidity; - if (std::abs(mcParticle.y()) > kMcRapidityWindow) { - continue; - } - - if (mcParticle.pdgCode() == PDG_t::kK0Short) { - rKzeroShort.fill(HIST("Generated_MCRecoColl_INEL_K0Short"), mcParticle.pt(), flattenicity); // K0s + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kK0Short) { + rKzeroShort.fill(HIST("Generated_MCRecoColl_INEL_K0Short"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rKzeroShort.fill(HIST("Generated_MCRecoColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); // K0s + rKzeroShort.fill(HIST("Generated_MCRecoColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); + rKzeroShort.fill(HIST("Generated_MCRecoColl_INELgt0_K0Short_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0) { - rLambda.fill(HIST("Generated_MCRecoColl_INEL_Lambda"), mcParticle.pt(), flattenicity); // Lambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0) { + rLambda.fill(HIST("Generated_MCRecoColl_INEL_Lambda"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rLambda.fill(HIST("Generated_MCRecoColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); // Lambda + rLambda.fill(HIST("Generated_MCRecoColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); + rLambda.fill(HIST("Generated_MCRecoColl_INELgt0_Lambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0Bar) { - rAntiLambda.fill(HIST("Generated_MCRecoColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0Bar) { + rAntiLambda.fill(HIST("Generated_MCRecoColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rAntiLambda.fill(HIST("Generated_MCRecoColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + rAntiLambda.fill(HIST("Generated_MCRecoColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); + rAntiLambda.fill(HIST("Generated_MCRecoColl_INELgt0_AntiLambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { - rXi.fill(HIST("Generated_MCRecoColl_INEL_Xi"), mcParticle.pt(), flattenicity); // Xi + if (inCascRapidity && std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { + rXi.fill(HIST("Generated_MCRecoColl_INEL_Xi"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rXi.fill(HIST("Generated_MCRecoColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); // Xi + rXi.fill(HIST("Generated_MCRecoColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); + rXi.fill(HIST("Generated_MCRecoColl_INELgt0_Xi_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } } @@ -1962,6 +2464,7 @@ struct Lambdak0sflattenicity { if (recoCollIndexINELgt0 > 0) { rEventSelection.fill(HIST("hNEventsMCGenReco"), 1.5); rEventSelection.fill(HIST("hFlat_GenRecoColl_MC_INELgt0"), flattenicity); + rEventSelection.fill(HIST("hFlat_GenRecoColl_MC_INELgt0_TrueFlat"), flattenicityTrue); } //===================================== @@ -1973,39 +2476,41 @@ struct Lambdak0sflattenicity { if (!mcParticle.isPhysicalPrimary()) { continue; } + const bool inV0Rapidity = std::abs(mcParticle.y()) <= v0Sel.v0settingRapidity; + const bool inCascRapidity = std::abs(mcParticle.y()) <= cascSel.cascsettingRapidity; - if (std::abs(mcParticle.y()) > kMcRapidityWindow) { - continue; - } - - if (mcParticle.pdgCode() == PDG_t::kK0Short) { - rKzeroShort.fill(HIST("pGen_MCGenRecoColl_INEL_K0Short"), mcParticle.pt(), flattenicity); // K0s + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kK0Short) { + rKzeroShort.fill(HIST("pGen_MCGenRecoColl_INEL_K0Short"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rKzeroShort.fill(HIST("pGen_MCGenRecoColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); // K0s + rKzeroShort.fill(HIST("pGen_MCGenRecoColl_INELgt0_K0Short"), mcParticle.pt(), flattenicity); + rKzeroShort.fill(HIST("pGen_MCGenRecoColl_INELgt0_K0Short_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0) { - rLambda.fill(HIST("pGen_MCGenRecoColl_INEL_Lambda"), mcParticle.pt(), flattenicity); // Lambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0) { + rLambda.fill(HIST("pGen_MCGenRecoColl_INEL_Lambda"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); // Lambda + rLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_Lambda"), mcParticle.pt(), flattenicity); + rLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_Lambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (mcParticle.pdgCode() == PDG_t::kLambda0Bar) { - rAntiLambda.fill(HIST("pGen_MCGenRecoColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + if (inV0Rapidity && mcParticle.pdgCode() == PDG_t::kLambda0Bar) { + rAntiLambda.fill(HIST("pGen_MCGenRecoColl_INEL_AntiLambda"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rAntiLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); // AntiLambda + rAntiLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_AntiLambda"), mcParticle.pt(), flattenicity); + rAntiLambda.fill(HIST("pGen_MCGenRecoColl_INELgt0_AntiLambda_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } - if (std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { - rXi.fill(HIST("pGen_MCGenRecoColl_INEL_Xi"), mcParticle.pt(), flattenicity); // Xi + if (inCascRapidity && std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { + rXi.fill(HIST("pGen_MCGenRecoColl_INEL_Xi"), mcParticle.pt(), flattenicity); if (recoCollIndexINELgt0 > 0) { - rXi.fill(HIST("pGen_MCGenRecoColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); // Xi + rXi.fill(HIST("pGen_MCGenRecoColl_INELgt0_Xi"), mcParticle.pt(), flattenicity); + rXi.fill(HIST("pGen_MCGenRecoColl_INELgt0_Xi_TrueFlat"), mcParticle.pt(), flattenicityTrue); } } } } // Cascade Analysis Starts here - using DauTracks = soa::Join; + using DauTracks = soa::Join; using LabeledDauTracks = soa::Join; using LabeledCascades = soa::Join; @@ -2029,14 +2534,9 @@ struct Lambdak0sflattenicity { } // track quality - if (posDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows || - negDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows || - bacDaughterTrack.tpcNClsCrossedRows() < cascSel.nTPCcrossedRows) { - return false; - } - if (std::abs(posDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || - std::abs(negDaughterTrack.eta()) > trkPid.cfgTrkEtaCut || - std::abs(bacDaughterTrack.eta()) > trkPid.cfgTrkEtaCut) { + if (!isSelectedDaughterTrack(posDaughterTrack, cascSel.nTPCcrossedRows) || + !isSelectedDaughterTrack(negDaughterTrack, cascSel.nTPCcrossedRows) || + !isSelectedDaughterTrack(bacDaughterTrack, cascSel.nTPCcrossedRows)) { return false; } @@ -2083,9 +2583,10 @@ struct Lambdak0sflattenicity { soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/) { + const bool own = (eventHistOwner == kOwnerDataCasc); if (evSel.applyEvSel && - !(isEventSelected(collision))) { // Checking if the event passes the - // selection criteria + !(isEventSelected(collision, own))) { // Checking if the event passes the + // selection criteria return; } @@ -2097,12 +2598,16 @@ struct Lambdak0sflattenicity { if (flattenicity < 0.f) { return; } - rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + if (own) { + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); - rEventSelection.fill(HIST("hVertexZ"), vtxZ); - rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); - rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); - rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + rEventSelection.fill(HIST("hVertexZ"), vtxZ); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracks, flattenicity); + } for (const auto& casc : Cascades) { @@ -2148,10 +2653,11 @@ struct Lambdak0sflattenicity { soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::FV0As const& /*fv0s*/, aod::McCollisions const&, aod::McParticles const& mcParticles) { + const bool own = (eventHistOwner == kOwnerRecMCCasc); for (const auto& collision : collisions) { if (evSel.applyEvSel && - !(isEventSelected(collision))) { // Checking if the event passes the - // selection criteria + !(isEventSelected(collision, own))) { // Checking if the event passes the + // selection criteria continue; } @@ -2168,16 +2674,29 @@ struct Lambdak0sflattenicity { if (flattenicity < 0.f) { continue; } - rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + if (own) { + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); - rEventSelection.fill(HIST("hVertexZ"), vtxZ); - rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); - rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); - rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + rEventSelection.fill(HIST("hVertexZ"), vtxZ); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracksThisCollision, flattenicity); + } auto cascsThisCollision = Cascades.sliceBy(perColCasc, collision.globalIndex()); const auto& mcCollision = collision.mcCollision_as(); + const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheCasc); + const float flattenicityMCGen = estimateFlattenicityFV0MC(particlesInCollision, own); + if (own) { + rEventSelection.fill(HIST("hFlattenicityDistributionMCGen_Rec"), flattenicityMCGen); + rEventSelection.fill(HIST("hFlattenicity_Corr_Gen_vs_Rec"), flattenicityMCGen, flattenicity); + rEventSelection.fill(HIST("hFlatGenVsRecFine"), flattenicityMCGen, flattenicity); + fillFlatDistGenInRec(flattenicityMCGen, collision.centFT0M(), collision.isInelGt0()); + } + for (const auto& casc : cascsThisCollision) { // MC truth matching, else the spectra keep the combinatorial background @@ -2206,6 +2725,7 @@ struct Lambdak0sflattenicity { const int motherIndex = getFeedDownMother(cascMcParticle, motherPt); rXi.fill(HIST("hMassFeedDownXipTFlat"), massXi, casc.pt(), flattenicity); rXi.fill(HIST("hFeedDownXiPtVsMotherPt"), casc.pt(), motherPt); + rXi.fill(HIST("hFeedDownXiMatrix"), casc.pt(), motherPt, flattenicity, motherIndex); rXi.fill(HIST("hFeedDownXiMotherPdg"), motherIndex); } if (!isPrimaryCasc && eventClass.requirePrimaryMC) { @@ -2224,33 +2744,108 @@ struct Lambdak0sflattenicity { rXi.fill(HIST("h2DdecayRadiusXi"), casc.cascradius()); rXi.fill(HIST("hMassXipT"), massXi, casc.pt()); rXi.fill(HIST("hMassXipTFlat"), massXi, casc.pt(), flattenicity); + rXi.fill(HIST("hMassXipTTrueFlat"), massXi, casc.pt(), flattenicityMCGen); + rXi.fill(HIST("hPtResXi"), cascMcParticle.pt(), + (casc.pt() - cascMcParticle.pt()) / cascMcParticle.pt()); rXi.fill(HIST("hNSigmaProtonFromXi"), protonDaughter.tpcNSigmaPr(), protonDaughter.tpcInnerParam()); rXi.fill(HIST("hNSigmaPionFromXi"), pionDaughter.tpcNSigmaPi(), pionDaughter.tpcInnerParam()); rXi.fill(HIST("hNSigmaBachPionFromXi"), bacDaughterTrack.tpcNSigmaPi(), bacDaughterTrack.tpcInnerParam()); } } - const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cacheCasc); - float flattenicityMCGen = estimateFlattenicityFV0MC(particlesInCollision); - rEventSelection.fill(HIST("hFlattenicityDistributionMCGen_Rec"), flattenicityMCGen); - rEventSelection.fill(HIST("hFlattenicity_Corr_Gen_vs_Rec"), flattenicityMCGen, flattenicity); - const bool isInelGt0Rec = collision.isInelGt0(); for (const auto& mcParticle : particlesInCollision) { - if (mcParticle.isPhysicalPrimary() && std::abs(mcParticle.y()) < kMcRapidityWindow && std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { - rXi.fill(HIST("Generated_MCRecoCollCheck_INEL_Xi"), mcParticle.pt(), flattenicity); // Xi + if (mcParticle.isPhysicalPrimary() && std::abs(mcParticle.y()) <= cascSel.cascsettingRapidity && + std::abs(mcParticle.pdgCode()) == PDG_t::kXiMinus) { + rXi.fill(HIST("Generated_MCRecoCollCheck_INEL_Xi"), mcParticle.pt(), flattenicity); if (isInelGt0Rec) { - rXi.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Xi"), mcParticle.pt(), flattenicity); // Xi + rXi.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Xi"), mcParticle.pt(), flattenicity); + rXi.fill(HIST("Generated_MCRecoCollCheck_INELgt0_Xi_TrueFlat"), mcParticle.pt(), flattenicityMCGen); } } } } } + // ================== Percentile determination pass ====================== // + // + // Event selection and flattenicity only, so a short run gives the 1-rho + // distribution the percentile boundaries are read off. The boundaries then go + // back into binning.axisFlat for the spectra pass, where one class is exactly + // one bin. + + void processFlatDistData( + soa::Join::iterator const& collision, + TrackCandidates const& tracks, + soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, + aod::FV0As const& /*fv0s*/) + { + const bool own = (eventHistOwner == kOwnerFlatDistData); + if (evSel.applyEvSel && !isEventSelected(collision, own)) { + return; + } + const float flattenicity = estimateFlattenicity(collision, tracks); + if (flattenicity < 0.f) { + return; + } + if (own) { + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + rEventSelection.fill(HIST("hVertexZ"), collision.posZ()); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracks, flattenicity); + } + } + + void processFlatDistMC( + soa::Join const& collisions, + aod::McCollisions const&, TrackCandidatesMC const& tracks, + soa::Join const& /*bcs*/, aod::FT0s const& /*ft0s*/, + aod::FV0As const& /*fv0s*/, aod::McParticles const& mcParticles) + { + const bool own = (eventHistOwner == kOwnerFlatDistMC); + for (const auto& collision : collisions) { + if (evSel.applyEvSel && !isEventSelected(collision, own)) { + continue; + } + if (!collision.has_mcCollision()) { + continue; + } + auto tracksThisCollision = tracks.sliceBy(perColTracksMC, collision.globalIndex()); + const float flattenicity = estimateFlattenicity(collision, tracksThisCollision); + if (flattenicity < 0.f) { + continue; + } + if (!own) { + continue; + } + rEventSelection.fill(HIST("hEventsSelected"), nbinFlattenicity - 0.5); + rEventSelection.fill(HIST("hVertexZ"), collision.posZ()); + rEventSelection.fill(HIST("hFlattenicityDistribution"), flattenicity); + rEventSelection.fill(HIST("hCentFT0M"), collision.centFT0M()); + rEventSelection.fill(HIST("hCentFT0MvsFlattenicity"), collision.centFT0M(), flattenicity); + fillFlatDistRec(flattenicity, collision.centFT0M(), collision.isInelGt0()); + fillChargedRec(tracksThisCollision, flattenicity); + + const auto& mcCollision = collision.mcCollision_as(); + const auto particlesInCollision = mcParticles.sliceByCached(aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache1); + const float flattenicityMCGen = estimateFlattenicityFV0MC(particlesInCollision); + rEventSelection.fill(HIST("hFlattenicityDistributionMCGen_Rec"), flattenicityMCGen); + rEventSelection.fill(HIST("hFlattenicity_Corr_Gen_vs_Rec"), flattenicityMCGen, flattenicity); + rEventSelection.fill(HIST("hFlatGenVsRecFine"), flattenicityMCGen, flattenicity); + fillFlatDistGenInRec(flattenicityMCGen, collision.centFT0M(), collision.isInelGt0()); + } + } + + PROCESS_SWITCH(Lambdak0sflattenicity, processFlatDistData, "Flattenicity distribution only, data", false); + PROCESS_SWITCH(Lambdak0sflattenicity, processFlatDistMC, "Flattenicity distribution only, MC", false); PROCESS_SWITCH(Lambdak0sflattenicity, processDataRun3LambdaK0s, "Process Run 3 Data LambdaK0s", false); PROCESS_SWITCH(Lambdak0sflattenicity, processRecMCLambdaK0s, "Process Run 3 MC reconstructed LambdaK0s", false); PROCESS_SWITCH(Lambdak0sflattenicity, processGenMC, "Process Run 3 MC generated", false); - PROCESS_SWITCH(Lambdak0sflattenicity, processDataRun3Cascade, "Process Run 3 Data Cascade", true); + PROCESS_SWITCH(Lambdak0sflattenicity, processDataRun3Cascade, "Process Run 3 Data Cascade", false); PROCESS_SWITCH(Lambdak0sflattenicity, processRecMCRun3Cascade, "Process Run 3 mc Rec Cascade", false); };