From c5c763b57b3413a82792ea9de33f95fd2a0b16d4 Mon Sep 17 00:00:00 2001 From: Sushanta Tripathy Date: Thu, 5 Nov 2026 13:36:54 +0100 Subject: [PATCH 1/2] Enabling multiplicity dependent studies --- PWGLF/Tasks/Resonances/f0phiproxy.cxx | 269 ++++++++++++++++++--- PWGLF/Tasks/Resonances/lambdastarproxy.cxx | 60 ++++- 2 files changed, 292 insertions(+), 37 deletions(-) diff --git a/PWGLF/Tasks/Resonances/f0phiproxy.cxx b/PWGLF/Tasks/Resonances/f0phiproxy.cxx index b83cdee0f7c..3d1aec4712f 100644 --- a/PWGLF/Tasks/Resonances/f0phiproxy.cxx +++ b/PWGLF/Tasks/Resonances/f0phiproxy.cxx @@ -13,6 +13,7 @@ /// \brief Dipion-based half-momentum kaon proxies with kaon/phi and pion/K*0 bachelors, in data and MC. /// \author Sushanta Tripathy +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponseTOF.h" @@ -34,7 +35,9 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include +#include #include #include #include @@ -46,10 +49,12 @@ #include #include #include +#include #include #include #include #include +#include #include #include @@ -57,7 +62,7 @@ using namespace o2; using namespace o2::framework; struct F0phiproxy { - using Collisions = soa::Join; + using Collisions = soa::Join; using Tracks = soa::Join; @@ -72,6 +77,10 @@ struct F0phiproxy { Configurable vertexMax{"vertexMax", 10.f, "Maximum absolute vertex z (cm)"}; Configurable requireSel8{"requireSel8", true, "Require sel8"}; Configurable requireINELgt0{"requireINELgt0", true, "Require at least one PV track in |eta|<1; MC gen uses charged physical primaries"}; + Configurable requireNoSameBunchPileup{"requireNoSameBunchPileup", false, "Require kNoSameBunchPileup, as in Lambda-star trigger mode 4"}; + Configurable requireGoodZvtxFT0vsPV{"requireGoodZvtxFT0vsPV", false, "Require kIsGoodZvtxFT0vsPV, as in Lambda-star trigger mode 4"}; + Configurable useFT0M{"useFT0M", true, "Use calibrated FT0M percentiles for reconstructed data/MC; false retains legacy NPV axes and mixing"}; + ConfigurableAxis ft0mAxis{"ft0mAxis", {0., 0., 1., 5., 10., 15., 20., 30., 40., 50., 70., 100.}, "FT0M percentile classes; first zero selects variable-width bins; must cover [0,100]"}; Configurable trackPtMin{"trackPtMin", 0.15f, "Minimum track pT (GeV/c)"}; Configurable trackEtaMax{"trackEtaMax", 0.8f, "Maximum track |eta|"}; Configurable requireGlobalTrack{"requireGlobalTrack", true, "Require global track selection"}; @@ -88,6 +97,7 @@ struct F0phiproxy { Configurable mixDepth{"mixDepth", 5, "Previous events per mixing bin; zero disables ME booking and mixing"}; Configurable mixZWidth{"mixZWidth", 2.f, "Vertex bin width (cm)"}; Configurable mixMultWidth{"mixMultWidth", 20.f, "PV multiplicity bin width"}; + Configurable mixFT0MWidth{"mixFT0MWidth", 1.f, "FT0M mixing width in percentage points; pools also respect analysis-class boundaries"}; Configurable enablePhiProxy{"enablePhiProxy", true, "Dipion/2 as on-shell K plus bachelor K"}; Configurable enableKstarProxy{"enableKstarProxy", false, "Dipion/2 as on-shell K plus bachelor pion; use the dedicated Kstar configuration"}; @@ -119,7 +129,11 @@ struct F0phiproxy { static constexpr int IonCodeThreshold = 1000000000; static constexpr int IonZDivisor = 10000; static constexpr int IonZModulo = 1000; + static constexpr float FT0MPercentileMax = 100.f; + static constexpr size_t MinAxisSpecificationSize = 3; + static constexpr double MaxUniformFT0MBins = 10000.; int resolvedF0Code = 0; + std::vector ft0mEdges{}; struct TruthInfo { int64_t id = -1; @@ -141,12 +155,28 @@ struct F0phiproxy { std::vector pions{}, kaons{}; }; + using MixingKey = std::tuple; + enum EventStage : int { Received = 1, + VertexAccepted, + Sel8Accepted, + InelAccepted, + QualityAccepted, + ActivityAccepted, + NEventStages = ActivityAccepted }; + std::shared_ptr eventSelection{}, eventStages{}, selectedEvents{}, ft0mStatusQA{}, mixPartners{}; + std::shared_ptr eventStagesByActivity{}, activityVsNPV{}, activityVsVertex{}, mixPartnersByActivity{}; + std::array, 4> pidQA{}; + std::array, 4> pidQAByActivity{}; + std::array, NMixClasses> dipionSpectra{}; std::array, NMixClasses>, NChannels> compactProxy{}, fullProxy{}; std::array, NCharges>, NMixClasses>, NChannels> controls{}; std::array, NTruthClasses>, NChannels> recoComponents{}; std::array, NChannels> recoTrueParent{}, recoPtResponse{}, proxyPtResponse{}; std::shared_ptr recoTrueF0{}; + std::array, NTruthClasses>, NChannels> recoComponentsByActivity{}; + std::array, NChannels> recoTrueParentByActivity{}, recoPtResponseByActivity{}, proxyPtResponseByActivity{}; + std::shared_ptr recoTrueF0ByActivity{}; enum GenSpectrum : size_t { GenDipion, GenF0, GenPhi, @@ -166,18 +196,97 @@ struct F0phiproxy { static double bachelorMass(size_t channel) { return channel == Phi ? constants::physics::MassKaonCharged : constants::physics::MassPionCharged; } static int bachelorPdg(size_t channel) { return channel == Phi ? PdgKaon : PdgPion; } + std::string recoName(std::string const& name) const + { + return useFT0M.value ? std::string("FT0M/").append(name) : name; + } + + static int ft0mStatus(float percentile) + { + if (!std::isfinite(percentile)) { + return 1; + } + if (percentile < 0.f) { + return 2; + } + return percentile < FT0MPercentileMax ? 0 : 3; + } + + static std::vector resolveFT0MEdges(std::vector const& specification) + { + if (specification.size() < MinAxisSpecificationSize || + !std::all_of(specification.begin(), specification.end(), [](double x) { return std::isfinite(x); })) { + throw std::runtime_error("Invalid FT0M axis specification"); + } + std::vector edges; + if (specification.front() == 0.) { + edges.assign(specification.begin() + 1, specification.end()); + } else { + const double bins = specification.front(); + if (specification.size() != MinAxisSpecificationSize || bins < 1. || bins > MaxUniformFT0MBins || std::floor(bins) != bins) { + throw std::runtime_error("Invalid uniform FT0M axis specification"); + } + for (int i = 0; i <= static_cast(bins); ++i) { + edges.push_back(specification[1] + (specification[2] - specification[1]) * i / bins); + } + } + if (edges.front() != 0. || edges.back() != FT0MPercentileMax || + std::adjacent_find(edges.begin(), edges.end(), [](double a, double b) { return a >= b; }) != edges.end()) { + throw std::runtime_error("FT0M classes must have increasing edges covering [0,100]"); + } + return edges; + } + + AxisSpec activityAxis(bool compact = false) const + { + if (useFT0M.value) { + return AxisSpec{ft0mEdges, "FT0M percentile (MB>0), calibrated"}; + } + return compact ? AxisSpec{proxyMultiplicityAxis, "N_{PV}(|#eta|<1), grouped"} : AxisSpec{multiplicityBins.value, 0., 200., "N_{PV}(|#eta|<1)"}; + } + + int ft0mClass(float percentile) const + { + if (ft0mStatus(percentile) != 0) { + return -1; + } + return static_cast(std::upper_bound(ft0mEdges.begin(), ft0mEdges.end(), percentile) - ft0mEdges.begin()) - 1; + } + + MixingKey mixingKey(int run, float vertex, float mult) const + { + const int activityClass = useFT0M.value ? ft0mClass(mult) : 0; + if (useFT0M.value && activityClass < 0) { + throw std::runtime_error("Invalid FT0M event cannot enter a mixing pool"); + } + const float width = useFT0M.value ? mixFT0MWidth.value : mixMultWidth.value; + auto bin = [](float value, float binWidth) { + const double index = std::floor(value / binWidth); + if (!std::isfinite(index) || index < std::numeric_limits::min() || index > std::numeric_limits::max()) { + throw std::runtime_error("Mixing coordinate cannot be represented by the configured bin width"); + } + return static_cast(index); + }; + return {run, bin(vertex, mixZWidth.value), activityClass, bin(mult, width)}; + } + void init(InitContext const&) { if (doprocess && doprocessMCReco) { throw std::runtime_error("Enable either process (data) or processMCReco, never both"); } - if (mixDepth.value < 0 || mixZWidth.value <= 0.f || mixMultWidth.value <= 0.f || + if (mixDepth.value < 0 || !std::isfinite(mixZWidth.value) || mixZWidth.value <= 0.f || + !std::isfinite(mixMultWidth.value) || mixMultWidth.value <= 0.f || + !std::isfinite(mixFT0MWidth.value) || mixFT0MWidth.value <= 0.f || vertexMax.value <= 0.f || trackPtMin.value < 0.f || trackEtaMax.value <= 0.f || pionTPC.value <= 0.f || pionTOF.value <= 0.f || kaonTPC.value <= 0.f || kaonTOF.value <= 0.f || dipionMassMin.value >= dipionMassMax.value || dipionYMax.value <= 0.f || proxyYMax.value <= 0.f || multiplicityBins.value <= 0 || proxyDipionMassBins.value <= 0 || kstarMassBins.value <= 0 || mcF0PdgCode.value < 0) { throw std::runtime_error("Invalid F0PhiProxy configuration"); } + if (useFT0M.value) { + ft0mEdges = resolveFT0MEdges(ft0mAxis.value); + } if (doprocess || doprocessMCReco) { if ((enablePhiProxy.value || enableKstarProxy.value) && !enableCompactProxy.value && !enableFullProxySparse.value) { throw std::runtime_error("Enable compact or full proxy output when a proxy channel is enabled"); @@ -202,28 +311,57 @@ struct F0phiproxy { const AxisSpec massPi{300, 0.2, 1.7, "m_{#pi#pi} (GeV/c^{2})"}; const AxisSpec compactMassPi{proxyDipionMassBins.value, 0.2, 1.7, "m_{#pi#pi} (GeV/c^{2})"}; const AxisSpec pt{100, 0., 10., "p_{T} (GeV/c)"}; - const AxisSpec mult{multiplicityBins.value, 0., 200., "N_{PV}(|#eta|<1)"}; - const AxisSpec compactMult{proxyMultiplicityAxis, "N_{PV}(|#eta|<1), grouped"}; + const AxisSpec mult = activityAxis(); + const AxisSpec compactMult = activityAxis(true); const AxisSpec compactPt{proxyPtAxis, "proxy pair p_{T} (GeV/c), grouped"}; const AxisSpec charge{3, -0.5, 2.5, "dipion: 0=unlike, 1=++, 2=--"}; const AxisSpec bachelorCharge{2, -1.5, 1.5, "bachelor charge"}; const AxisSpec tof{2, -0.5, 1.5, "bachelor has TOF"}; - histos.add("hEventSelection", "Received; selected", HistType::kTH1D, {{2, 0.5, 2.5}}); + eventSelection = histos.add(recoName("hEventSelection").c_str(), "Received; selected including activity validity", HistType::kTH1D, {{2, 0.5, 2.5}}); + eventSelection->GetXaxis()->SetBinLabel(1, "Received"); + eventSelection->GetXaxis()->SetBinLabel(2, "Selected"); + const AxisSpec stages{NEventStages, 0.5, NEventStages + 0.5, "Event selection stage"}; + eventStages = histos.add(recoName("hEventSelectionSteps").c_str(), "Reconstructed event cut flow", HistType::kTH1D, {stages}); + eventStagesByActivity = histos.add(recoName("hEventSelectionVsMultiplicity").c_str(), "Valid-activity collisions;activity;selection stage", HistType::kTH2D, {compactMult, stages}); + const std::array stageNames{"Received", "Vertex", "sel8 or disabled", "INEL>0 or disabled", "Quality flags or disabled", "Valid activity selected"}; + for (int i = 0; i < NEventStages; ++i) { + eventStages->GetXaxis()->SetBinLabel(i + 1, stageNames[i]); + eventStagesByActivity->GetYaxis()->SetBinLabel(i + 1, stageNames[i]); + } + selectedEvents = histos.add(recoName("hSelectedEvents").c_str(), "Selected collisions counted once, including zero-candidate events;activity;Events", HistType::kTH1D, {compactMult}); + activityVsNPV = histos.add(recoName("hMultiplicityVsNPV").c_str(), "Selected collisions;activity;N_{PV}(|#eta|<1), overflow retained", HistType::kTH2D, {compactMult, {400, 0., 400.}}); + activityVsVertex = histos.add(recoName("hMultiplicityVsVertex").c_str(), "Selected collisions;activity;vertex z (cm)", HistType::kTH2D, {compactMult, {100, -vertexMax.value, vertexMax.value}}); + auto estimator = histos.add("hMultiplicityEstimator", "Configured reconstructed estimator;mode;Workers", HistType::kTH1D, {{2, -0.5, 1.5}}); + estimator->GetXaxis()->SetBinLabel(1, "NPV"); + estimator->GetXaxis()->SetBinLabel(2, "FT0M"); + estimator->Fill(useFT0M.value ? 1. : 0.); + if (useFT0M.value) { + ft0mStatusQA = histos.add("FT0M/hPercentileStatus", "After base event selection;FT0M status;Events", HistType::kTH1D, {{4, -0.5, 3.5}}); + const std::array statusNames{"Valid [0,100)", "Nonfinite", "Below zero", "At or above 100"}; + for (size_t i = 0; i < statusNames.size(); ++i) { + ft0mStatusQA->GetXaxis()->SetBinLabel(static_cast(i + 1), statusNames[i]); + } + } if (mixingEnabled()) { - histos.add("hMixPartners", "Mixing partners;N", HistType::kTH1D, {{101, -0.5, 100.5}}); + mixPartnersByActivity = histos.add(recoName("hMixPartnersVsMultiplicity").c_str(), "Available previous events;activity;mixing partners", HistType::kTH2D, {compactMult, {101, -0.5, 100.5}}); + mixPartners = histos.add(recoName("hMixPartners").c_str(), "Mixing partners;N", HistType::kTH1D, {{101, -0.5, 100.5}}); } if (enablePidQA.value) { - histos.add("hPionTPC", "Quality-selected tracks;pT;nSigmaTPC pion", HistType::kTH2F, {pt, {100, -10., 10.}}); - histos.add("hPionTOF", "TOF-matched quality tracks;pT;nSigmaTOF pion", HistType::kTH2F, {pt, {100, -10., 10.}}); - histos.add("hKaonTPC", "Quality-selected tracks;pT;nSigmaTPC kaon", HistType::kTH2F, {pt, {100, -10., 10.}}); - histos.add("hKaonTOF", "TOF-matched quality tracks;pT;nSigmaTOF kaon", HistType::kTH2F, {pt, {100, -10., 10.}}); + const std::array names{"hPionTPC", "hPionTOF", "hKaonTPC", "hKaonTOF"}; + for (size_t i = 0; i < names.size(); ++i) { + pidQA[i] = histos.add(recoName(names[i]).c_str(), "Quality-selected tracks, before PID;pT;nSigma", HistType::kTH2F, {pt, {100, -10., 10.}}); + if (useFT0M.value) { + const std::string name = std::string(names[i]).append("ByMultiplicity"); + pidQAByActivity[i] = histos.add(recoName(name).c_str(), "Quality-selected tracks, before PID;pT;nSigma;FT0M percentile", HistType::kTH3F, {pt, {100, -10., 10.}, compactMult}); + } + } } const std::array dipionNames{"hDipionSE", "hDipionME"}; for (size_t mix = 0; mix < dipionNames.size(); ++mix) { if (mix != 0 && !mixingEnabled()) { continue; } - dipionSpectra[mix] = histos.add(dipionNames[mix], "Dipion SE or raw ME", HistType::kTHnSparseF, {massPi, pt, mult, charge}); + dipionSpectra[mix] = histos.add(recoName(dipionNames[mix]).c_str(), "Dipion SE or raw ME", HistType::kTHnSparseF, {massPi, pt, mult, charge}); } const std::array, NChannels> compactNames{{{"hProxyCompactSE", "hProxyCompactME"}, {"hKstarProxyCompactSE", "hKstarProxyCompactME"}}}; const std::array, NChannels> fullNames{{{"hProxySE", "hProxyME"}, {"hKstarProxySE", "hKstarProxyME"}}}; @@ -238,15 +376,15 @@ struct F0phiproxy { if (channelEnabled(channel)) { // The measured dipion mass is retained; a mass window is not an f0 tag. if (enableCompactProxy.value) { - compactProxy[channel][mix] = histos.add(compactNames[channel][mix], "Grouped-pT/multiplicity proxy; no constituent charge is measured", HistType::kTHnSparseF, {mass, compactPt, compactMassPi, compactMult, charge, bachelorCharge, tof}); + compactProxy[channel][mix] = histos.add(recoName(compactNames[channel][mix]).c_str(), "Grouped-pT/multiplicity proxy; no constituent charge is measured", HistType::kTHnSparseF, {mass, compactPt, compactMassPi, compactMult, charge, bachelorCharge, tof}); } if (enableFullProxySparse.value) { - fullProxy[channel][mix] = histos.add(fullNames[channel][mix], "Diagnostic full proxy", HistType::kTHnSparseF, {mass, pt, massPi, mult, charge, bachelorCharge, tof}); + fullProxy[channel][mix] = histos.add(recoName(fullNames[channel][mix]).c_str(), "Diagnostic full proxy", HistType::kTHnSparseF, {mass, pt, massPi, mult, charge, bachelorCharge, tof}); } } if (controlEnabled(channel)) { for (size_t sign = 0; sign < NCharges; ++sign) { - controls[channel][mix][sign] = histos.add(controlNames[channel][mix][sign], "Standard control SE or raw ME", HistType::kTHnSparseF, {mass, pt, mult}); + controls[channel][mix][sign] = histos.add(recoName(controlNames[channel][mix][sign]).c_str(), "Standard control SE or raw ME", HistType::kTHnSparseF, {mass, pt, mult}); } } } @@ -263,27 +401,40 @@ struct F0phiproxy { const std::array names{"hMCGenDipion", "hMCGenTrueF0", "hMCGenPhiKK", "hMCGenTruePhi", "hMCGenKstarKPi", "hMCGenTrueKstar", "hMCGenProxyPhi", "hMCGenProxyPhiTrueF0", "hMCGenProxyKstar", "hMCGenProxyKstarTrueF0"}; for (size_t i = 0; i < names.size(); ++i) { const AxisSpec axis = i <= GenF0 ? pionMass : (i == GenKstar || i == GenKstarTrue || i >= GenProxyKstar ? mass[Kstar] : mass[Phi]); - generated[i] = histos.add(names[i], "Generated unweighted fiducial spectrum", HistType::kTH2D, {axis, pt}, true); + generated[i] = histos.add(names[i], "Generated unweighted fiducial spectrum; no detector FT0M selection", HistType::kTH2D, {axis, pt}, true); } } if (doprocessMCReco) { - recoTrueF0 = histos.add("hMCRecoTrueF0", "Selected common-parent f0 -> pi pi", HistType::kTH2D, {pionMass, pt}, true); + recoTrueF0 = histos.add(recoName("hMCRecoTrueF0").c_str(), "Selected common-parent f0 -> pi pi", HistType::kTH2D, {pionMass, pt}, true); + if (useFT0M.value) { + recoTrueF0ByActivity = histos.add("FT0M/hMCRecoTrueF0ByMultiplicity", "Selected common-parent f0;dipion mass;pair pT;FT0M percentile", HistType::kTHnSparseD, {pionMass, pt, activityAxis(true)}, true); + } const std::array channelNames{"Phi", "Kstar"}; const std::array classNames{"TrueF0", "Combinatorial", "MisID", "Unmatched"}; for (size_t channel = 0; channel < NChannels; ++channel) { if (channelEnabled(channel)) { for (size_t category = 0; category < classNames.size(); ++category) { - const std::string name = std::string("hMCRecoProxy") + channelNames[channel] + classNames[category]; - recoComponents[channel][category] = histos.add(name.c_str(), "Exclusive proxy truth component; not a coalescence ancestry tag", HistType::kTH2D, {mass[channel], pt}, true); + const std::string name = std::string("hMCRecoProxy").append(channelNames[channel]).append(classNames[category]); + recoComponents[channel][category] = histos.add(recoName(name).c_str(), "Exclusive proxy truth component; not a coalescence ancestry tag", HistType::kTH2D, {mass[channel], pt}, true); + if (useFT0M.value) { + recoComponentsByActivity[channel][category] = histos.add(recoName(std::string(name).append("ByMultiplicity")).c_str(), "Exclusive proxy truth component; not coalescence ancestry", HistType::kTHnSparseD, {mass[channel], pt, activityAxis(true)}, true); + } + } + const std::string name = std::string("hMCRecoProxy").append(channelNames[channel]).append("PtResponse"); + proxyPtResponse[channel] = histos.add(recoName(name).c_str(), "Primary true f0 daughters and true bachelor;generated proxy pair pT;reconstructed proxy pair pT", HistType::kTH2D, {pt, pt}, true); + if (useFT0M.value) { + proxyPtResponseByActivity[channel] = histos.add(recoName(std::string(name).append("ByMultiplicity")).c_str(), "Selected true-f0 proxy response;generated pair pT;reconstructed pair pT;FT0M percentile", HistType::kTHnSparseD, {pt, pt, activityAxis(true)}, true); } - const std::string name = std::string("hMCRecoProxy") + channelNames[channel] + "PtResponse"; - proxyPtResponse[channel] = histos.add(name.c_str(), "Primary true f0 daughters and true bachelor;generated proxy pair pT;reconstructed proxy pair pT", HistType::kTH2D, {pt, pt}, true); } if (controlEnabled(channel)) { - const std::string truthName = std::string("hMCRecoTrue") + channelNames[channel]; - const std::string responseName = std::string("hMCReco") + channelNames[channel] + "PtResponse"; - recoTrueParent[channel] = histos.add(truthName.c_str(), "Selected exact common-parent standard control", HistType::kTH2D, {mass[channel], pt}, true); - recoPtResponse[channel] = histos.add(responseName.c_str(), "Selected standard control;generated pair pT;reconstructed pair pT", HistType::kTH2D, {pt, pt}, true); + const std::string truthName = std::string("hMCRecoTrue").append(channelNames[channel]); + const std::string responseName = std::string("hMCReco").append(channelNames[channel]).append("PtResponse"); + recoTrueParent[channel] = histos.add(recoName(truthName).c_str(), "Selected exact common-parent standard control", HistType::kTH2D, {mass[channel], pt}, true); + recoPtResponse[channel] = histos.add(recoName(responseName).c_str(), "Selected standard control;generated pair pT;reconstructed pair pT", HistType::kTH2D, {pt, pt}, true); + if (useFT0M.value) { + recoTrueParentByActivity[channel] = histos.add(recoName(std::string(truthName).append("ByMultiplicity")).c_str(), "Selected exact common-parent standard control", HistType::kTHnSparseD, {mass[channel], pt, activityAxis(true)}, true); + recoPtResponseByActivity[channel] = histos.add(recoName(std::string(responseName).append("ByMultiplicity")).c_str(), "Selected standard response;generated pair pT;reconstructed pair pT;FT0M percentile", HistType::kTHnSparseD, {pt, pt, activityAxis(true)}, true); + } } } } @@ -398,11 +549,13 @@ struct F0phiproxy { if (mc && !mixed) { const size_t truthClass = proxyTruthClass(first, second, bachelor, channel); recoComponents[channel][truthClass]->Fill(pair.M(), pair.Pt()); + fillSparse(recoComponentsByActivity[channel][truthClass], {pair.M(), pair.Pt(), mult}); if (truthClass == 0 && first.truth.primary && second.truth.primary && bachelor.truth.primary && bachelor.truth.id != first.truth.id && bachelor.truth.id != second.truth.id) { const double px = ProxyMomentumScale * (first.truth.px + second.truth.px) + bachelor.truth.px; const double py = ProxyMomentumScale * (first.truth.py + second.truth.py) + bachelor.truth.py; proxyPtResponse[channel]->Fill(std::hypot(px, py), pair.Pt()); + fillSparse(proxyPtResponseByActivity[channel], {std::hypot(px, py), pair.Pt(), mult}); } } } @@ -422,6 +575,8 @@ struct F0phiproxy { if (mc && !mixed && first.sign != second.sign && identified(first, PdgKaon) && identified(second, bachelorPdg(channel)) && commonParent(first.truth, second.truth, parent)) { recoTrueParent[channel]->Fill(pair.M(), pair.Pt()); recoPtResponse[channel]->Fill(truthPairPt(first.truth, second.truth), pair.Pt()); + fillSparse(recoTrueParentByActivity[channel], {pair.M(), pair.Pt(), mult}); + fillSparse(recoPtResponseByActivity[channel], {truthPairPt(first.truth, second.truth), pair.Pt(), mult}); } } void fillControls(Event const& event, std::deque const& pool, float mult, bool mc) @@ -467,21 +622,56 @@ struct F0phiproxy { void processReco(TCollisions const& collisions, TTracks const& tracks, TSlice const& slice) { // Pools reset at each dataframe; the key also separates runs, z and multiplicity. - std::map>, std::deque> pools; + std::map> pools; const std::deque emptyPool; for (auto const& collision : collisions) { - histos.fill(HIST("hEventSelection"), 1.); - if (std::abs(collision.posZ()) >= vertexMax.value || (requireSel8.value && !collision.sel8()) || - (requireINELgt0.value && collision.multNTracksPVeta1() < 1)) { + const float npv = collision.multNTracksPVeta1(); + const float mult = useFT0M.value ? collision.centFT0M() : npv; + const bool validActivity = useFT0M.value ? ft0mStatus(mult) == 0 : std::isfinite(mult) && mult >= 0.f; + auto stage = [&](int step) { + eventStages->Fill(step); + if (validActivity) { + eventStagesByActivity->Fill(mult, step); + } + }; + eventSelection->Fill(1.); + stage(Received); + if (!std::isfinite(collision.posZ()) || std::abs(collision.posZ()) >= vertexMax.value) { + continue; + } + stage(VertexAccepted); + if (requireSel8.value && !collision.sel8()) { + continue; + } + stage(Sel8Accepted); + if (requireINELgt0.value && npv < 1.f) { continue; } - histos.fill(HIST("hEventSelection"), 2.); - const float mult = collision.multNTracksPVeta1(); + stage(InelAccepted); + if ((requireNoSameBunchPileup.value && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) || + (requireGoodZvtxFT0vsPV.value && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))) { + continue; + } + stage(QualityAccepted); + if (useFT0M.value) { + // Missing calibration (105) and non-INEL markers (100.5) are not percentiles. + // Never substitute NPV in an FT0M histogram or pool. + ft0mStatusQA->Fill(ft0mStatus(mult)); + } + if (!validActivity) { + continue; + } + stage(ActivityAccepted); + eventSelection->Fill(2.); + selectedEvents->Fill(mult); + activityVsNPV->Fill(mult, npv); + activityVsVertex->Fill(mult, collision.posZ()); std::deque* activePool = nullptr; if (mixingEnabled()) { const int run = collision.template bc_as().runNumber(); - activePool = &pools[{run, {static_cast(std::floor(collision.posZ() / mixZWidth.value)), static_cast(std::floor(mult / mixMultWidth.value))}}]; - histos.fill(HIST("hMixPartners"), activePool->size()); + activePool = &pools[mixingKey(run, collision.posZ(), mult)]; + mixPartners->Fill(activePool->size()); + mixPartnersByActivity->Fill(mult, activePool->size()); } auto const& pool = activePool ? *activePool : emptyPool; Event event; @@ -493,11 +683,17 @@ struct F0phiproxy { } const bool tof = track.hasTOF(); if (enablePidQA.value) { - histos.fill(HIST("hPionTPC"), track.pt(), track.tpcNSigmaPi()); - histos.fill(HIST("hKaonTPC"), track.pt(), track.tpcNSigmaKa()); + auto fillPID = [&](size_t index, float value) { + pidQA[index]->Fill(track.pt(), value); + if (pidQAByActivity[index]) { + pidQAByActivity[index]->Fill(track.pt(), value, mult); + } + }; + fillPID(0, track.tpcNSigmaPi()); + fillPID(2, track.tpcNSigmaKa()); if (tof) { - histos.fill(HIST("hPionTOF"), track.pt(), track.tofNSigmaPi()); - histos.fill(HIST("hKaonTOF"), track.pt(), track.tofNSigmaKa()); + fillPID(1, track.tofNSigmaPi()); + fillPID(3, track.tofNSigmaKa()); } } TrackCandidate candidate{.id = track.globalIndex(), .sign = track.sign(), .px = track.px(), .py = track.py(), .pz = track.pz(), .hasTOF = tof, .truth = {}}; @@ -543,6 +739,7 @@ struct F0phiproxy { fillSparse(dipionSpectra[0], {pair.M(), pair.Pt(), mult, static_cast(chargeCategory(first.sign, second.sign))}); if (mc && trueF0(first, second)) { recoTrueF0->Fill(pair.M(), pair.Pt()); + fillSparse(recoTrueF0ByActivity, {pair.M(), pair.Pt(), mult}); } if (pair.M() < dipionMassMin.value || pair.M() >= dipionMassMax.value) { continue; diff --git a/PWGLF/Tasks/Resonances/lambdastarproxy.cxx b/PWGLF/Tasks/Resonances/lambdastarproxy.cxx index 2ff5ad926a9..625403f7c0e 100644 --- a/PWGLF/Tasks/Resonances/lambdastarproxy.cxx +++ b/PWGLF/Tasks/Resonances/lambdastarproxy.cxx @@ -74,6 +74,16 @@ struct Lambdastarproxy { static constexpr float CutPtMinDefault = 0.5f; static constexpr float CutEtaMaxDefault = 0.8f; + static constexpr int FT0MPercentileBins = 100; + static constexpr float FT0MPercentileMax = 100.f; + enum class FT0MEventStatus : int { + Valid = 1, + Negative, + OutOfRange, + NonFinite, + Unavailable + }; + static constexpr float CutMCPtMinDefault = 0.1f; static constexpr float CutMCEtaMaxDefault = 0.8f; static constexpr int FillMCTruthDefault = 1; @@ -352,6 +362,37 @@ struct Lambdastarproxy { return fallback; } + template + void fillEventMultiplicityCounters(TCollision const& collision, float eventMult) + { + // Use the exact estimator of the pair maps, including the legacy track-count + // fallback. Nonfinite values have no percentile and are counted in overflow. + histos.fill(HIST("hEventMultiplicity"), std::isfinite(eventMult) ? eventMult : FT0MPercentileMax); + + // Keep the raw FT0M denominator separate: fallback track counts must never + // acquire an FT0M percentile label. Missing/nonfinite FT0M goes to overflow; + // finite negative and >=100 values retain their underflow/overflow meaning. + float ft0m = FT0MPercentileMax; + auto status = FT0MEventStatus::Unavailable; + if constexpr (requires { collision.centFT0M(); }) { + const float value = collision.centFT0M(); + if (!std::isfinite(value)) { + status = FT0MEventStatus::NonFinite; + } else { + ft0m = value; + if (value < 0.f) { + status = FT0MEventStatus::Negative; + } else if (value >= FT0MPercentileMax) { + status = FT0MEventStatus::OutOfRange; + } else { + status = FT0MEventStatus::Valid; + } + } + } + histos.fill(HIST("hEventFT0M"), ft0m); + histos.fill(HIST("hEventFT0MStatus"), static_cast(status)); + } + int chargeFromPdg(int pdg) const { const int apdg = std::abs(pdg); @@ -529,13 +570,28 @@ struct Lambdastarproxy { AxisSpec dcaXYAxis{200, -0.2, 0.2, "DCA_{xy} (cm)"}; AxisSpec dcaZAxis{200, -0.2, 0.2, "DCA_{z} (cm)"}; - // Count reconstructed collisions received by the data process, before any + // Count reconstructed collisions received by the data or reco-MC process, before any // task event cuts and after all of them, independently of track candidates. histos.add("hEventSelection", "Proxy event selection;Selection stage;Events", HistType::kTH1D, {AxisSpec{2, 0.5, 2.5, "Selection stage"}}); histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(1, "Before event selection"); histos.get(HIST("hEventSelection"))->GetXaxis()->SetBinLabel(2, "After event selection"); + // Small event-level denominators, independent of candidate and QA switches. + // Double storage and Sumw2 preserve large merged event counts and errors. + histos.add("hEventFT0M", "Selected reconstructed events;Raw FT0M percentile;Events", + HistType::kTH1D, {AxisSpec{FT0MPercentileBins, 0., FT0MPercentileMax, "Raw FT0M percentile"}}, true); + histos.add("hEventMultiplicity", "Selected reconstructed events (same estimator as pair maps);Multiplicity/centrality;Events", + HistType::kTH1D, {centAxis}, true); + histos.add("hEventFT0MStatus", "Selected reconstructed events;Raw FT0M status;Events", + HistType::kTH1D, {AxisSpec{static_cast(FT0MEventStatus::Unavailable), 0.5, static_cast(FT0MEventStatus::Unavailable) + 0.5, "Raw FT0M status"}}, true); + auto* ft0mStatusAxis = histos.get(HIST("hEventFT0MStatus"))->GetXaxis(); + ft0mStatusAxis->SetBinLabel(static_cast(FT0MEventStatus::Valid), "Finite FT0M [0,100)"); + ft0mStatusAxis->SetBinLabel(static_cast(FT0MEventStatus::Negative), "Negative FT0M"); + ft0mStatusAxis->SetBinLabel(static_cast(FT0MEventStatus::OutOfRange), "FT0M >= 100"); + ft0mStatusAxis->SetBinLabel(static_cast(FT0MEventStatus::NonFinite), "Nonfinite FT0M"); + ft0mStatusAxis->SetBinLabel(static_cast(FT0MEventStatus::Unavailable), "FT0M unavailable"); + // Invariant-mass spectra histos.add("hInvMassPKUnlike", "pK invariant mass (unlike-sign);M_{pK} (GeV/c^{2});Counts", @@ -1945,6 +2001,8 @@ struct Lambdastarproxy { // Compute event multiplicity (FT0M or fallback) const float eventMult = eventMultiplicityFT0MOrFallback(collision, eventMultFallback); + // Exactly once after event cuts, before PID probes or candidate requirements. + fillEventMultiplicityCounters(collision, eventMult); if (lstarEnableProxyPIDProbes.value != 0) { if constexpr (IsMC) { fillProxyPIDProbes(tracks, eventMult, collision.mcCollisionId()); From 0c520f69cf066b8459c5b3f12571ff4e034de7aa Mon Sep 17 00:00:00 2001 From: Sushanta Tripathy Date: Thu, 5 Nov 2026 13:50:15 +0100 Subject: [PATCH 2/2] Fix static_cast for event selection stage axis --- PWGLF/Tasks/Resonances/f0phiproxy.cxx | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/PWGLF/Tasks/Resonances/f0phiproxy.cxx b/PWGLF/Tasks/Resonances/f0phiproxy.cxx index 3d1aec4712f..1792f20da65 100644 --- a/PWGLF/Tasks/Resonances/f0phiproxy.cxx +++ b/PWGLF/Tasks/Resonances/f0phiproxy.cxx @@ -13,6 +13,7 @@ /// \brief Dipion-based half-momentum kaon proxies with kaon/phi and pion/K*0 bachelors, in data and MC. /// \author Sushanta Tripathy +#include "Common/CCDB/EventSelectionParams.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" @@ -320,7 +321,7 @@ struct F0phiproxy { eventSelection = histos.add(recoName("hEventSelection").c_str(), "Received; selected including activity validity", HistType::kTH1D, {{2, 0.5, 2.5}}); eventSelection->GetXaxis()->SetBinLabel(1, "Received"); eventSelection->GetXaxis()->SetBinLabel(2, "Selected"); - const AxisSpec stages{NEventStages, 0.5, NEventStages + 0.5, "Event selection stage"}; + const AxisSpec stages{NEventStages, 0.5, static_cast(NEventStages) + 0.5, "Event selection stage"}; eventStages = histos.add(recoName("hEventSelectionSteps").c_str(), "Reconstructed event cut flow", HistType::kTH1D, {stages}); eventStagesByActivity = histos.add(recoName("hEventSelectionVsMultiplicity").c_str(), "Valid-activity collisions;activity;selection stage", HistType::kTH2D, {compactMult, stages}); const std::array stageNames{"Received", "Vertex", "sel8 or disabled", "INEL>0 or disabled", "Quality flags or disabled", "Valid activity selected"};