diff --git a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx index 6d9e1ca4fdf..6edc236a27f 100644 --- a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx +++ b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx @@ -14,7 +14,7 @@ /// \brief Analysis to do PID with MC - Full correction factors for pions, kaons, protons #include "PWGLF/DataModel/LFParticleIdentification.h" -#include "PWGLF/DataModel/mcCentrality.h" // For McCentFT0Ms +#include "PWGLF/DataModel/mcCentrality.h" #include "PWGLF/DataModel/spectraTOF.h" #include "PWGLF/Utils/inelGt.h" @@ -71,8 +71,6 @@ struct MultiplicityPt { static constexpr int CentBinMax = 100; static constexpr int MultBinMax = 200; static constexpr int RecMultBinMax = 100; - static constexpr int DebugCountMax = 20; - static constexpr int CentMultClasses = 10; enum INELCutSelection : int { INEL = 0, @@ -83,21 +81,10 @@ struct MultiplicityPt { Configurable isRun3{"isRun3", true, "is Run3 dataset"}; Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; Configurable cfgINELCut{"cfgINELCut", 0, "INEL event selection: 0 no sel, 1 INEL>0, 2 INEL>1"}; - Configurable askForCustomTVX{"askForCustomTVX", false, "Ask for custom TVX rather than sel8"}; - Configurable removeITSROFrameBorder{"removeITSROFrameBorder", false, "Remove ITS Read-Out Frame border"}; - Configurable removeNoSameBunchPileup{"removeNoSameBunchPileup", false, "Remove no same bunch pileup"}; - Configurable requireIsGoodZvtxFT0vsPV{"requireIsGoodZvtxFT0vsPV", false, "Require good Z vertex FT0 vs PV"}; - Configurable requireIsVertexITSTPC{"requireIsVertexITSTPC", false, "Require vertex ITSTPC"}; - Configurable removeNoTimeFrameBorder{"removeNoTimeFrameBorder", false, "Remove no time frame border"}; Configurable cfgCutEtaMax{"cfgCutEtaMax", 0.8f, "Max eta range for tracks"}; Configurable cfgCutEtaMin{"cfgCutEtaMin", -0.8f, "Min eta range for tracks"}; - Configurable cfgCutY{"cfgCutY", 0.5f, "Y range for tracks"}; Configurable cfgCutNsigma{"cfgCutNsigma", 3.0f, "nsigma cut range for tracks"}; - Configurable lastRequiredTrdCluster{"lastRequiredTrdCluster", -1, "Last cluster to require in TRD"}; - Configurable requireTrdOnly{"requireTrdOnly", false, "Require only tracks from TRD"}; - Configurable requireNoTrd{"requireNoTrd", false, "Require tracks without TRD"}; - Configurable enableDCAHistograms{"enableDCAHistograms", false, "Enable DCA histograms"}; Configurable enablePIDHistograms{"enablePIDHistograms", true, "Enable PID histograms"}; Configurable useCustomTrackCuts{"useCustomTrackCuts", true, "Flag to use custom track cuts"}; Configurable itsPattern{"itsPattern", 0, "0 = Run3ITSibAny, 1 = Run3ITSallAny, 2 = Run3ITSall7Layers, 3 = Run3ITSibTwo"}; @@ -123,7 +110,6 @@ struct MultiplicityPt { Configurable phiCutHighParam1{"phiCutHighParam1", 0.16685, "First parameter for high phi cut"}; Configurable phiCutHighParam2{"phiCutHighParam2", 0.00981942, "Second parameter for high phi cut"}; - Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum constituent pT"}; Configurable cfgCutNsigmaPi{"cfgCutNsigmaPi", 3.0f, "nsigma cut for pions"}; @@ -131,7 +117,6 @@ struct MultiplicityPt { Configurable cfgCutNsigmaPr{"cfgCutNsigmaPr", 2.5f, "nsigma cut for protons"}; TrackSelection customTrackCuts; - TF1* fphiCutLow = nullptr; TF1* fphiCutHigh = nullptr; @@ -171,12 +156,10 @@ struct MultiplicityPt { float getTransformedPhi(const float phi, const int charge, const float magField) const { float transformedPhi = phi; - if (magField < 0) { + if (magField < 0) transformedPhi = o2::constants::math::TwoPI - transformedPhi; - } - if (charge < 0) { + if (charge < 0) transformedPhi = o2::constants::math::TwoPI - transformedPhi; - } transformedPhi += o2::constants::math::PI / 18.0f; transformedPhi = std::fmod(transformedPhi, o2::constants::math::PI / 9.0f); return transformedPhi; @@ -190,7 +173,6 @@ struct MultiplicityPt { if (track.pt() < pTthresholdPhiCut.value) return true; - float pt = track.pt(); float phi = track.phi(); int charge = track.sign(); @@ -198,11 +180,10 @@ struct MultiplicityPt { phi = o2::constants::math::TwoPI - phi; if (charge < 0) phi = o2::constants::math::TwoPI - phi; - phi += o2::constants::math::PI / 18.0f; phi = std::fmod(phi, o2::constants::math::PI / 9.0f); - if (phi < fphiCutHigh->Eval(pt) && phi > fphiCutLow->Eval(pt)) + if (phi < fphiCutHigh->Eval(track.pt()) && phi > fphiCutLow->Eval(track.pt())) return false; return true; } @@ -289,51 +270,24 @@ struct MultiplicityPt { return true; } - template - bool passesPIDSelection(TrackType const& track) const + template + bool passesPIDSelection(const TrackType& track, int species) const { float nsigmaTPC = 0.f; - if constexpr (species == PartPion) - nsigmaTPC = track.tpcNSigmaPi(); - else if constexpr (species == PartKaon) - nsigmaTPC = track.tpcNSigmaKa(); - else if constexpr (species == PartProton) - nsigmaTPC = track.tpcNSigmaPr(); + float cutValue = 0.f; - float cutValue = cfgCutNsigma.value; - if constexpr (species == PartPion) + if (species == PartPion) { + nsigmaTPC = track.tpcNSigmaPi(); cutValue = cfgCutNsigmaPi.value; - if constexpr (species == PartKaon) + } else if (species == PartKaon) { + nsigmaTPC = track.tpcNSigmaKa(); cutValue = cfgCutNsigmaKa.value; - if constexpr (species == PartProton) + } else if (species == PartProton) { + nsigmaTPC = track.tpcNSigmaPr(); cutValue = cfgCutNsigmaPr.value; - - return (std::abs(nsigmaTPC) < cutValue); - } - - template - int getBestPIDHypothesis(TrackType const& track) const - { - float nsigmaPi = std::abs(track.tpcNSigmaPi()); - float nsigmaKa = std::abs(track.tpcNSigmaKa()); - float nsigmaPr = std::abs(track.tpcNSigmaPr()); - - float minNSigma = 999.0f; - int bestSpecies = -1; - - if (nsigmaPi < cfgCutNsigmaPi.value && nsigmaPi < minNSigma) { - minNSigma = nsigmaPi; - bestSpecies = PartPion; - } - if (nsigmaKa < cfgCutNsigmaKa.value && nsigmaKa < minNSigma) { - minNSigma = nsigmaKa; - bestSpecies = PartKaon; } - if (nsigmaPr < cfgCutNsigmaPr.value && nsigmaPr < minNSigma) { - minNSigma = nsigmaPr; - bestSpecies = PartProton; - } - return bestSpecies; + + return std::abs(nsigmaTPC) < cutValue; } template @@ -393,15 +347,16 @@ void MultiplicityPt::init(InitContext const&) customTrackCuts.print(); } + // Axis definitions ConfigurableAxis ptBinning{"ptBinning", {VARIABLE_WIDTH, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.1, 2.2, 2.3, 2.4, 2.5, 2.6, 2.7, 2.8, 2.9, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0, 8.0, 9.0, 10.0, 11.0, 12.0, 13.0, 14.0, 15.0, 16.0, 18.0, 20.0}, "pT bin limits"}; AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; std::vector centBinningStd = {0., 1., 5., 10., 15., 20., 30., 40., 50., 60., 70., 80., 90., 100.}; + AxisSpec centAxis = {centBinningStd, "FT0M Centrality (%)"}; + std::vector centBinningFine; for (int i = 0; i <= CentBinMax; i++) centBinningFine.push_back(static_cast(i)); - - AxisSpec centAxis = {centBinningStd, "FT0M Centrality (%)"}; AxisSpec centFineAxis = {centBinningFine, "FT0M Centrality (%)"}; std::vector multBins; @@ -414,164 +369,98 @@ void MultiplicityPt::init(InitContext const&) recoMultBins.push_back(static_cast(i)); AxisSpec recoMultAxis = {recoMultBins, "N_{ch}^{reco}"}; - ue.add("Centrality/hCentRaw", "Raw FT0M Centrality (no cuts);Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); + ue.add("Centrality/hCentRaw", "Raw FT0M Centrality;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); ue.add("Centrality/hCentAfterVtx", "Centrality after vertex cut;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); - ue.add("Centrality/hCentAfterINEL", "Centrality after INEL cut;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); ue.add("Centrality/hCentAfterAll", "Centrality after all cuts;Centrality (%);Counts", HistType::kTH1D, {centFineAxis}); ue.add("Centrality/hCentVsMult", "Centrality vs Generated Multiplicity;Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centFineAxis, multAxis}); ue.add("Centrality/hMultVsCent", "Generated Multiplicity vs Centrality;N_{ch}^{gen};Centrality (%)", HistType::kTH2D, {multAxis, centFineAxis}); - ue.add("Centrality/hCentVsVz", "Centrality vs Vertex Z;Centrality (%);V_{z} (cm)", HistType::kTH2D, {centFineAxis, {40, -20, 20}}); ue.add("Centrality/hRecoMultVsCent", "Reconstructed Track Multiplicity vs Centrality;Centrality (%);N_{tracks}^{reco}", HistType::kTH2D, {centFineAxis, recoMultAxis}); - ue.add("Centrality/hGenMultPerCent", "Generated Multiplicity Distribution per Centrality Bin;Centrality (%);", HistType::kTH2D, {centFineAxis, multAxis}); - ue.add("Centrality/hVertexResVsCent", "Vertex Resolution vs Centrality;Centrality (%);V_{z} resolution (cm)", HistType::kTH2D, {centFineAxis, {100, -1, 1}}); - ue.add("INEL/hINELClass", "INEL Class for MC Collisions;INEL Class;Counts", HistType::kTH1D, {{3, 0.5, 3.5}}); - auto hINEL = ue.get(HIST("INEL/hINELClass")); - hINEL->GetXaxis()->SetBinLabel(1, "INEL0"); - hINEL->GetXaxis()->SetBinLabel(2, "INEL>0"); - hINEL->GetXaxis()->SetBinLabel(3, "INEL>1"); - - ue.add("INEL/hINELVsCent", "INEL Class vs Centrality;Centrality (%);INEL Class", HistType::kTH2D, {centFineAxis, {3, 0.5, 3.5}}); - - ue.add("CutFlow/hCutStats", "Cut Statistics;Cut Stage;Counts", HistType::kTH1D, {{6, 0.5, 6.5}}); + ue.add("CutFlow/hCutStats", "Cut Statistics;Cut Stage;Counts", HistType::kTH1D, {{5, 0.5, 5.5}}); auto hCut = ue.get(HIST("CutFlow/hCutStats")); - hCut->GetXaxis()->SetBinLabel(1, "All reco events"); + hCut->GetXaxis()->SetBinLabel(1, "All collisions"); hCut->GetXaxis()->SetBinLabel(2, "Has MC match"); hCut->GetXaxis()->SetBinLabel(3, "Has centrality"); hCut->GetXaxis()->SetBinLabel(4, "Pass vertex"); - hCut->GetXaxis()->SetBinLabel(5, "Pass INEL"); - hCut->GetXaxis()->SetBinLabel(6, "Selected"); - - ue.add("CutFlow/hCentPerCut", "Centrality Distribution at Each Cut;Cut Stage;Centrality (%)", HistType::kTH2D, {{6, 0.5, 6.5}, centFineAxis}); + hCut->GetXaxis()->SetBinLabel(5, "Selected"); - ue.add("MC/GenRecoCollisions", "Generated and Reconstructed MC Collisions", HistType::kTH1D, {{10, 0.5, 10.5}}); - auto hColl = ue.get(HIST("MC/GenRecoCollisions")); - hColl->GetXaxis()->SetBinLabel(1, "Collisions generated"); - hColl->GetXaxis()->SetBinLabel(2, "Collisions reconstructed"); - hColl->GetXaxis()->SetBinLabel(3, "INEL>0"); - hColl->GetXaxis()->SetBinLabel(4, "INEL>1"); - - ue.add("hEventLossBreakdown", "Event loss breakdown", HistType::kTH1D, {{4, 0.5, 4.5}}); + ue.add("hEventLossBreakdown", "Event loss breakdown", HistType::kTH1D, {{3, 0.5, 3.5}}); auto hLoss = ue.get(HIST("hEventLossBreakdown")); hLoss->GetXaxis()->SetBinLabel(1, "Physics selected"); hLoss->GetXaxis()->SetBinLabel(2, "Reconstructed"); hLoss->GetXaxis()->SetBinLabel(3, "Selected"); - hLoss->GetXaxis()->SetBinLabel(4, "Final efficiency"); - ue.add("MC/EventLoss/NchGenerated", "Generated charged multiplicity;N_{ch}^{gen};Counts", HistType::kTH1D, {{200, 0, 200}}); - ue.add("MC/EventLoss/NchGenerated_PhysicsSelected", "Generated charged multiplicity (physics selected);N_{ch}^{gen};Counts", HistType::kTH1D, {{200, 0, 200}}); - ue.add("MC/EventLoss/NchGenerated_Reconstructed", "Generated charged multiplicity (reconstructed);N_{ch}^{gen};Counts", HistType::kTH1D, {{200, 0, 200}}); + ue.add("MC/GenRecoCollisions", "Generated and Reconstructed MC Collisions", HistType::kTH1D, {{5, 0.5, 5.5}}); + auto hColl = ue.get(HIST("MC/GenRecoCollisions")); + hColl->GetXaxis()->SetBinLabel(1, "Collisions generated"); + hColl->GetXaxis()->SetBinLabel(2, "INEL>0"); + hColl->GetXaxis()->SetBinLabel(3, "INEL>1"); + hColl->GetXaxis()->SetBinLabel(4, "Reconstructed"); + hColl->GetXaxis()->SetBinLabel(5, "Selected"); + ue.add("MC/EventLoss/GenMultVsCent", "Generated charged particles vs FT0M centrality;FT0M Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centAxis, multAxis}); ue.add("MC/EventLoss/GenMultVsCent_Selected", "Generated vs FT0M centrality (selected events);FT0M Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centAxis, multAxis}); - ue.add("MC/EventLoss/GenMultVsCent_Rejected", "Generated vs FT0M centrality (rejected events);FT0M Centrality (%);N_{ch}^{gen}", HistType::kTH2D, {centAxis, multAxis}); - ue.add("MC/GenPtVsNch", "Generated pT vs Multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, {200, 0, 200}}); - ue.add("MC/GenPtVsNch_PhysicsSelected", "Generated pT vs Multiplicity (physics selected);#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, {200, 0, 200}}); + ue.add("Inclusive/hPtGen", "Generated primaries (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Inclusive/hPtReco", "All reconstructed tracks (track selection only);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Inclusive/hPtPrimReco", "Reconstructed primaries (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Inclusive/hPtSecReco", "Reconstructed secondaries (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimGenAll", "All generated primaries (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimGen", "Generated primaries (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimBadVertex", "Generated primaries (bad vertex);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimRecoEv", "Generated primaries (reco events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimGoodEv", "Generated primaries (good events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - - // Tracking efficiency - ue.add("Inclusive/hPtDenEff", "Tracking efficiency denominator (generated primaries in selected events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtNumEff", "Tracking efficiency numerator (reconstructed primaries matched to generated);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - - // Primary fraction - ue.add("Inclusive/hPtAllReco", "All reconstructed tracks;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtPrimReco", "Reconstructed primaries;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtSecReco", "Reconstructed secondaries;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtPrimGenAll", "All generated #pi^{#pm} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtGenINEL", "Generated #pi^{#pm} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtGenRecoEvent", "Generated #pi^{#pm} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtGen", "Generated #pi^{#pm} (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtReco", "Reconstructed #pi^{#pm} (MC matched, any status);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtPrimReco", "Reconstructed primary #pi^{#pm} (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Pion/hPtMeasured", "Measured #pi^{#pm} (PID selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - // Multiplicity-dependent - ue.add("Inclusive/hPtDenEffVsCent", "Generated primaries vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Inclusive/hPtNumEffVsCent", "Reconstructed primaries matched vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Inclusive/hPtPrimRecoVsCent", "Reconstructed primaries vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Inclusive/hPtAllRecoVsCent", "All reconstructed tracks vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Inclusive/hPtMeasuredVsCent", "All measured tracks (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); + ue.add("Pion/hPtGenRecoEvent_Mult", "Generated #pi^{#pm} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); + ue.add("Pion/hPtGenINEL_Mult", "Generated #pi^{#pm} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); - ue.add("Pion/hPtPrimGenAll", "All generated #pi^{#pm} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimGen", "Generated #pi^{#pm} (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimBadVertex", "Generated #pi^{#pm} (bad vertex);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimRecoEv", "Generated #pi^{#pm} (reco events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimGoodEv", "Generated #pi^{#pm} (good events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtDenEff", "#pi^{#pm} tracking efficiency denominator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtNumEff", "#pi^{#pm} tracking efficiency numerator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtAllReco", "All reconstructed #pi^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtPrimReco", "Reconstructed primary #pi^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtSecReco", "Reconstructed secondary #pi^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Pion/hPtDenEffVsCent", "Generated #pi^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Pion/hPtNumEffVsCent", "Reconstructed #pi^{#pm} matched vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Pion/hPtPrimRecoVsCent", "Reconstructed primary #pi^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Pion/hPtAllRecoVsCent", "All reconstructed #pi^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Pion/hPtMeasuredVsCent", "Measured #pi^{#pm} (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - if (enablePIDHistograms) + if (enablePIDHistograms) { ue.add("Pion/hNsigmaTPC", "TPC n#sigma #pi^{#pm};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); + } ue.add("Kaon/hPtPrimGenAll", "All generated K^{#pm} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtPrimGen", "Generated K^{#pm} (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtPrimBadVertex", "Generated K^{#pm} (bad vertex);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtPrimRecoEv", "Generated K^{#pm} (reco events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtPrimGoodEv", "Generated K^{#pm} (good events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtDenEff", "K^{#pm} tracking efficiency denominator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtNumEff", "K^{#pm} tracking efficiency numerator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtAllReco", "All reconstructed K^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtPrimReco", "Reconstructed primary K^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtSecReco", "Reconstructed secondary K^{#pm};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Kaon/hPtDenEffVsCent", "Generated K^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Kaon/hPtNumEffVsCent", "Reconstructed K^{#pm} matched vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Kaon/hPtPrimRecoVsCent", "Reconstructed primary K^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Kaon/hPtAllRecoVsCent", "All reconstructed K^{#pm} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Kaon/hPtMeasuredVsCent", "Measured K^{#pm} (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - if (enablePIDHistograms) + ue.add("Kaon/hPtGenINEL", "Generated K^{#pm} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Kaon/hPtGenRecoEvent", "Generated K^{#pm} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Kaon/hPtGen", "Generated K^{#pm} (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Kaon/hPtReco", "Reconstructed K^{#pm} (MC matched, any status);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Kaon/hPtPrimReco", "Reconstructed primary K^{#pm} (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Kaon/hPtMeasured", "Measured K^{#pm} (PID selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + + ue.add("Kaon/hPtGenRecoEvent_Mult", "Generated K^{#pm} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); + ue.add("Kaon/hPtGenINEL_Mult", "Generated K^{#pm} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); + + if (enablePIDHistograms) { ue.add("Kaon/hNsigmaTPC", "TPC n#sigma K^{#pm};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); + } ue.add("Proton/hPtPrimGenAll", "All generated p+#bar{p} (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtPrimGen", "Generated p+#bar{p} (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtPrimBadVertex", "Generated p+#bar{p} (bad vertex);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtPrimRecoEv", "Generated p+#bar{p} (reco events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtPrimGoodEv", "Generated p+#bar{p} (good events);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtDenEff", "p+#bar{p} tracking efficiency denominator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtNumEff", "p+#bar{p} tracking efficiency numerator;#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtAllReco", "All reconstructed p+#bar{p};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtPrimReco", "Reconstructed primary p+#bar{p};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtSecReco", "Reconstructed secondary p+#bar{p};#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - ue.add("Proton/hPtDenEffVsCent", "Generated p+#bar{p} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Proton/hPtNumEffVsCent", "Reconstructed p+#bar{p} matched vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Proton/hPtPrimRecoVsCent", "Reconstructed primary p+#bar{p} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Proton/hPtAllRecoVsCent", "All reconstructed p+#bar{p} vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - ue.add("Proton/hPtMeasuredVsCent", "Measured p+#bar{p} (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", HistType::kTH2D, {ptAxis, centAxis}); - if (enablePIDHistograms) + ue.add("Proton/hPtGenINEL", "Generated p+#bar{p} in INEL>0 events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtGenRecoEvent", "Generated p+#bar{p} in reconstructed events (|#eta|<0.8, p_{T}>0.15);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtGen", "Generated p+#bar{p} (physics selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtReco", "Reconstructed p+#bar{p} (MC matched, any status);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtPrimReco", "Reconstructed primary p+#bar{p} (MC matched);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + ue.add("Proton/hPtMeasured", "Measured p+#bar{p} (PID selected);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); + + ue.add("Proton/hPtGenRecoEvent_Mult", "Generated p+#bar{p} in reconstructed events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); + ue.add("Proton/hPtGenINEL_Mult", "Generated p+#bar{p} in INEL>0 events vs multiplicity;#it{p}_{T} (GeV/#it{c});N_{ch}^{gen}", HistType::kTH2D, {ptAxis, multAxis}); + + if (enablePIDHistograms) { ue.add("Proton/hNsigmaTPC", "TPC n#sigma p+#bar{p};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); + } + + ue.add("hEventsReco_Cent", "Reconstructed events vs centrality;FT0M Centrality (%);Counts", HistType::kTH1D, {centAxis}); + ue.add("hEventsINEL_Cent", "INEL>0 events vs centrality;FT0M Centrality (%);Counts", HistType::kTH1D, {centAxis}); ue.add("hNclFoundTPC", "Number of TPC found clusters", HistType::kTH1D, {{200, 0, 200}}); ue.add("hNclPIDTPC", "Number of TPC PID clusters", HistType::kTH1D, {{200, 0, 200}}); - ue.add("hNclFoundTPCvsPt", "TPC found clusters vs pT;#it{p}_{T} (GeV/#it{c});N_{cl,found}", HistType::kTH2D, {ptAxis, {200, 0., 200.}}); - ue.add("hNclPIDTPCvsPt", "TPC PID clusters vs pT;#it{p}_{T} (GeV/#it{c});N_{cl,PID}", HistType::kTH2D, {ptAxis, {200, 0., 200.}}); ue.add("hEta", "Track eta;#eta;Counts", HistType::kTH1D, {{20, -0.8, 0.8}}); ue.add("hPhi", "Track phi;#varphi (rad);Counts", HistType::kTH1D, {{64, 0, TwoPI}}); ue.add("hvtxZ", "Vertex Z (data);Vertex Z (cm);Events", HistType::kTH1F, {{40, -20.0, 20.0}}); - ue.add("hvtxZmc", "MC vertex Z;Vertex Z (cm);Events", HistType::kTH1F, {{40, -20.0, 20.0}}); - - ue.add("evsel", "Event selection", HistType::kTH1D, {{20, 0.5, 20.5}}); - auto hEvSel = ue.get(HIST("evsel")); - hEvSel->GetXaxis()->SetBinLabel(1, "Events read"); - hEvSel->GetXaxis()->SetBinLabel(4, "Trigger passed"); - hEvSel->GetXaxis()->SetBinLabel(5, "NoITSROFrameBorder"); - hEvSel->GetXaxis()->SetBinLabel(6, "NoSameBunchPileup"); - hEvSel->GetXaxis()->SetBinLabel(7, "IsGoodZvtxFT0vsPV"); - hEvSel->GetXaxis()->SetBinLabel(8, "IsVertexITSTPC"); - hEvSel->GetXaxis()->SetBinLabel(9, "NoTimeFrameBorder"); - hEvSel->GetXaxis()->SetBinLabel(13, "posZ passed"); - - // Phi cut monitoring - if (applyPhiCut.value) { - ue.add("PhiCut/hPtVsPhiPrimeBefore", "pT vs φ' before cut;p_{T} (GeV/c);φ'", HistType::kTH2F, {{100, 0, 10}, {100, 0, 0.4}}); - ue.add("PhiCut/hPtVsPhiPrimeAfter", "pT vs φ' after cut;p_{T} (GeV/c);φ'", HistType::kTH2F, {{100, 0, 10}, {100, 0, 0.4}}); - ue.add("PhiCut/hRejectionRate", "Track rejection rate by phi cut;p_{T} (GeV/c);Rejection Rate", HistType::kTProfile, {{100, 0, 10}}); - } - LOG(info) << "=== Initialization complete - All correction factor histograms defined ==="; + LOG(info) << "=== Initialization complete ==="; } void MultiplicityPt::processMC(TrackTableMC const& tracks, @@ -582,12 +471,12 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, aod::McCentFT0Ms const& centTable, BCsRun3 const& /*bcs*/) { - LOG(info) << "\n=== processMC START ==="; std::map mcCollisionToNch; - std::map mcCollisionVz; std::set physicsSelectedMCCollisions; std::map mcCollisionToINELClass; + std::set inel0MCCollisions; + std::map mcCollToCentFromReco; // MC collision ID -> centrality from reco ue.fill(HIST("MC/GenRecoCollisions"), 1.f, mcCollisions.size()); @@ -595,18 +484,34 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, int64_t mcCollId = mcCollision.globalIndex(); auto particlesInCollision = particles.sliceBy(perMCCol, mcCollId); + // Check if event has at least 1 particle in acceptance (INEL>0) + bool hasParticleInAcceptance = false; + for (const auto& particle : particlesInCollision) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + hasParticleInAcceptance = true; + break; + } + + if (hasParticleInAcceptance) { + inel0MCCollisions.insert(mcCollId); + } + int nGenCharged = countGeneratedChargedPrimaries(particlesInCollision, cfgCutEtaMax.value, cfgTrkLowPtCut.value); mcCollisionToNch[mcCollId] = nGenCharged; - mcCollisionVz[mcCollId] = mcCollision.posZ(); bool inel0 = o2::pwglf::isINELgt0mc(particlesInCollision, pdg); bool inel1 = o2::pwglf::isINELgt1mc(particlesInCollision, pdg); int inelClass = inel1 ? 2 : (inel0 ? 1 : 0); mcCollisionToINELClass[mcCollId] = inelClass; - ue.fill(HIST("INEL/hINELClass"), inelClass); - ue.fill(HIST("MC/EventLoss/NchGenerated"), nGenCharged); - bool physicsSelected = (std::abs(mcCollision.posZ()) <= cfgCutVertex.value); if (cfgINELCut.value == INELgt0 && !inel0) physicsSelected = false; @@ -615,120 +520,68 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, if (physicsSelected) { physicsSelectedMCCollisions.insert(mcCollId); - ue.fill(HIST("MC/EventLoss/NchGenerated_PhysicsSelected"), nGenCharged); if (inel0) - ue.fill(HIST("MC/GenRecoCollisions"), 3.f); + ue.fill(HIST("MC/GenRecoCollisions"), 2.f); if (inel1) - ue.fill(HIST("MC/GenRecoCollisions"), 4.f); - } - } - - // DEBUG: Count raw protons in MC (no cuts) - { - int nProtonsRaw = 0; - int nProtonsTotal = 0; - int nProtonsPassEta = 0; - int nProtonsPassPt = 0; - int nProtonsPassBoth = 0; - - for (const auto& mcCollision : mcCollisions) { - auto particlesInCollision = particles.sliceBy(perMCCol, mcCollision.globalIndex()); - for (const auto& particle : particlesInCollision) { - auto pdgParticle = pdg->GetParticle(particle.pdgCode()); - if (!pdgParticle || pdgParticle->Charge() == 0.) - continue; - - int pdgCode = std::abs(particle.pdgCode()); - if (pdgCode == PDG_t::kProton) { - nProtonsRaw++; // Count ALL protons regardless of status - - if (!particle.isPhysicalPrimary()) - continue; - nProtonsTotal++; - - bool passEta = (std::abs(particle.eta()) < cfgCutEtaMax.value); - bool passPt = (particle.pt() >= cfgTrkLowPtCut.value); - if (passEta) - nProtonsPassEta++; - if (passPt) - nProtonsPassPt++; - if (passEta && passPt) - nProtonsPassBoth++; - } - } - } - LOG(info) << "=== PROTON DEBUG ==="; - LOG(info) << "RAW protons in MC (any status): " << nProtonsRaw; - LOG(info) << "Physical primary protons: " << nProtonsTotal; - LOG(info) << "Protons passing eta cut (|eta|<" << cfgCutEtaMax.value << "): " << nProtonsPassEta; - LOG(info) << "Protons passing pT cut (pT>=" << cfgTrkLowPtCut.value << "): " << nProtonsPassPt; - LOG(info) << "Protons passing both cuts: " << nProtonsPassBoth; - if (nProtonsRaw == 0) { - LOG(warning) << "NO PROTONS FOUND IN MC! Check your MC generator settings."; + ue.fill(HIST("MC/GenRecoCollisions"), 3.f); } - } - - LOG(info) << "\n--- FILLING GENERATED PARTICLE SPECTRA ---"; - - for (const auto& mcCollision : mcCollisions) { - int64_t mcCollId = mcCollision.globalIndex(); - auto nchIt = mcCollisionToNch.find(mcCollId); - if (nchIt == mcCollisionToNch.end()) - continue; - int nGenCharged = nchIt->second; - bool isPhysicsSelected = (physicsSelectedMCCollisions.find(mcCollId) != physicsSelectedMCCollisions.end()); - auto particlesInCollision = particles.sliceBy(perMCCol, mcCollId); + // Fill generated particle spectra for (const auto& particle : particlesInCollision) { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); if (!pdgParticle || pdgParticle->Charge() == 0.) continue; if (!particle.isPhysicalPrimary()) continue; - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) - continue; - if (particle.pt() < cfgTrkLowPtCut.value) - continue; int pdgCode = std::abs(particle.pdgCode()); float pt = particle.pt(); - // All generated (no cuts) - SIGNAL LOSS DENOMINATOR - ue.fill(HIST("Inclusive/hPtPrimGenAll"), pt); - ue.fill(HIST("MC/GenPtVsNch"), pt, nGenCharged); - if (pdgCode == PDG_t::kPiPlus) + // Fill hPtPrimGenAll for ALL generated particles (NO CUTS) + if (pdgCode == PDG_t::kPiPlus) { ue.fill(HIST("Pion/hPtPrimGenAll"), pt); - else if (pdgCode == PDG_t::kKPlus) + } else if (pdgCode == PDG_t::kKPlus) { ue.fill(HIST("Kaon/hPtPrimGenAll"), pt); - else if (pdgCode == PDG_t::kProton) + } else if (pdgCode == PDG_t::kProton) { ue.fill(HIST("Proton/hPtPrimGenAll"), pt); + } + + // Fill hPtGenINEL for particles in INEL>0 events (with acceptance cuts) + if (hasParticleInAcceptance) { + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; - // Physics selected - SIGNAL LOSS NUMERATOR & EFFICIENCY DENOMINATOR - if (isPhysicsSelected) { - ue.fill(HIST("Inclusive/hPtPrimGen"), pt); - ue.fill(HIST("Inclusive/hPtDenEff"), pt); - ue.fill(HIST("MC/GenPtVsNch_PhysicsSelected"), pt, nGenCharged); if (pdgCode == PDG_t::kPiPlus) { - ue.fill(HIST("Pion/hPtPrimGen"), pt); - ue.fill(HIST("Pion/hPtDenEff"), pt); + ue.fill(HIST("Pion/hPtGenINEL"), pt); + ue.fill(HIST("Pion/hPtGenINEL_Mult"), pt, nGenCharged); } else if (pdgCode == PDG_t::kKPlus) { - ue.fill(HIST("Kaon/hPtPrimGen"), pt); - ue.fill(HIST("Kaon/hPtDenEff"), pt); + ue.fill(HIST("Kaon/hPtGenINEL"), pt); + ue.fill(HIST("Kaon/hPtGenINEL_Mult"), pt, nGenCharged); } else if (pdgCode == PDG_t::kProton) { - ue.fill(HIST("Proton/hPtPrimGen"), pt); - ue.fill(HIST("Proton/hPtDenEff"), pt); + ue.fill(HIST("Proton/hPtGenINEL"), pt); + ue.fill(HIST("Proton/hPtGenINEL_Mult"), pt, nGenCharged); } } - // Bad vertex for signal loss study - if (std::abs(mcCollision.posZ()) > cfgCutVertex.value) { - ue.fill(HIST("Inclusive/hPtPrimBadVertex"), pt); - if (pdgCode == PDG_t::kPiPlus) - ue.fill(HIST("Pion/hPtPrimBadVertex"), pt); - else if (pdgCode == PDG_t::kKPlus) - ue.fill(HIST("Kaon/hPtPrimBadVertex"), pt); - else if (pdgCode == PDG_t::kProton) - ue.fill(HIST("Proton/hPtPrimBadVertex"), pt); + // Apply acceptance cuts for physics-selected spectra + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + + // Fill generated spectra (physics selected only) + if (physicsSelected) { + ue.fill(HIST("Inclusive/hPtGen"), pt); + + if (pdgCode == PDG_t::kPiPlus) { + ue.fill(HIST("Pion/hPtGen"), pt); + } else if (pdgCode == PDG_t::kKPlus) { + ue.fill(HIST("Kaon/hPtGen"), pt); + } else if (pdgCode == PDG_t::kProton) { + ue.fill(HIST("Proton/hPtGen"), pt); + } } } } @@ -736,55 +589,34 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, std::map recoToMcMap; std::map recoToCentMap; - size_t iLabel = 0; - size_t nCollisions = collisions.size(); - for (const auto& label : labels) { - if (iLabel < nCollisions) { - const auto& collision = collisions.iteratorAt(iLabel); - recoToMcMap[collision.globalIndex()] = label.mcCollisionId(); - } - iLabel++; + size_t nPairs = std::min(labels.size(), collisions.size()); + for (size_t i = 0; i < nPairs; ++i) { + const auto& collision = collisions.iteratorAt(i); + const auto& label = labels.iteratorAt(i); + recoToMcMap[collision.globalIndex()] = label.mcCollisionId(); } - size_t iCent = 0; - for (const auto& cent : centTable) { - if (iCent < nCollisions) { - const auto& collision = collisions.iteratorAt(iCent); - float centValue = cent.centFT0M(); - recoToCentMap[collision.globalIndex()] = centValue; - ue.fill(HIST("Centrality/hCentRaw"), centValue); - } - iCent++; + size_t nCentPairs = std::min(centTable.size(), collisions.size()); + for (size_t i = 0; i < nCentPairs; ++i) { + const auto& collision = collisions.iteratorAt(i); + const auto& cent = centTable.iteratorAt(i); + float centValue = cent.centFT0M(); + recoToCentMap[collision.globalIndex()] = centValue; + ue.fill(HIST("Centrality/hCentRaw"), centValue); } - LOG(info) << "\n--- PROCESSING RECONSTRUCTED COLLISIONS ---"; - std::set reconstructedMCCollisions; std::set selectedMCCollisions; - // For cut statistics - std::vector centAll, centVertex, centINEL, centSelected; - int nRecoCollisions = 0; - int nSelectedEvents = 0; - int nRejectedEvents = 0; - int nNoMCMatch = 0; - int nNoCent = 0; - int nInvalidCent = 0; - int nPassVertex = 0; - int nPassINEL = 0; - int nPassAll = 0; - for (const auto& collision : collisions) { - nRecoCollisions++; - int64_t collId = collision.globalIndex(); - ue.fill(HIST("CutFlow/hCutStats"), 1); + int64_t collId = collision.globalIndex(); + + // MC matching auto mcIt = recoToMcMap.find(collId); - if (mcIt == recoToMcMap.end()) { - nNoMCMatch++; + if (mcIt == recoToMcMap.end()) continue; - } ue.fill(HIST("CutFlow/hCutStats"), 2); int64_t mcCollId = mcIt->second; @@ -796,69 +628,46 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, auto inelIt = mcCollisionToINELClass.find(mcCollId); int inelClass = (inelIt != mcCollisionToINELClass.end()) ? inelIt->second : 0; + // Centrality auto centIt = recoToCentMap.find(collId); - if (centIt == recoToCentMap.end()) { - nNoCent++; + if (centIt == recoToCentMap.end()) continue; - } ue.fill(HIST("CutFlow/hCutStats"), 3); + float cent = centIt->second; - if (cent < 0 || cent > CentBinMax) { - nInvalidCent++; + if (cent < 0 || cent > CentBinMax) continue; - } - // Store all events with valid info - centAll.push_back(cent); - ue.fill(HIST("Centrality/hCentVsMult"), cent, nGenCharged); - ue.fill(HIST("Centrality/hMultVsCent"), nGenCharged, cent); - ue.fill(HIST("Centrality/hCentVsVz"), cent, collision.posZ()); - ue.fill(HIST("INEL/hINELVsCent"), cent, inelClass); + // Store centrality for this MC collision (used later for MC truth plots) + mcCollToCentFromReco[mcCollId] = cent; - // Track cuts progressively + // Vertex cut bool passVertex = std::abs(collision.posZ()) <= cfgCutVertex.value; - if (passVertex) { - centVertex.push_back(cent); - ue.fill(HIST("Centrality/hCentAfterVtx"), cent); - ue.fill(HIST("CutFlow/hCutStats"), 4); - ue.fill(HIST("CutFlow/hCentPerCut"), 4, cent); - nPassVertex++; - } + if (!passVertex) + continue; + ue.fill(HIST("CutFlow/hCutStats"), 4); + ue.fill(HIST("Centrality/hCentAfterVtx"), cent); + // INEL cut bool passINEL = true; if (cfgINELCut.value == INELgt0 && inelClass < INELgt0) passINEL = false; if (cfgINELCut.value == INELgt1 && inelClass < INELgt1) passINEL = false; - if (passINEL) { - centINEL.push_back(cent); - ue.fill(HIST("Centrality/hCentAfterINEL"), cent); - ue.fill(HIST("CutFlow/hCutStats"), 5); - ue.fill(HIST("CutFlow/hCentPerCut"), 5, cent); - nPassINEL++; - } - - ue.fill(HIST("MC/EventLoss/GenMultVsCent"), cent, nGenCharged); - ue.fill(HIST("MC/EventLoss/NchGenerated_Reconstructed"), nGenCharged); - reconstructedMCCollisions.insert(mcCollId); - - bool passedAll = passVertex && passINEL; - if (!passedAll) { - ue.fill(HIST("MC/EventLoss/GenMultVsCent_Rejected"), cent, nGenCharged); - nRejectedEvents++; + if (!passINEL) continue; - } - // Event passed all selections - centSelected.push_back(cent); + // Event passed all cuts + ue.fill(HIST("CutFlow/hCutStats"), 5); ue.fill(HIST("Centrality/hCentAfterAll"), cent); - ue.fill(HIST("CutFlow/hCutStats"), 6); - ue.fill(HIST("CutFlow/hCentPerCut"), 6, cent); + ue.fill(HIST("MC/EventLoss/GenMultVsCent_Selected"), cent, nGenCharged); + ue.fill(HIST("Centrality/hCentVsMult"), cent, nGenCharged); + ue.fill(HIST("Centrality/hMultVsCent"), nGenCharged, cent); ue.fill(HIST("hvtxZ"), collision.posZ()); + selectedMCCollisions.insert(mcCollId); - nSelectedEvents++; - nPassAll++; + reconstructedMCCollisions.insert(mcCollId); // Get magnetic field for phi cut float magField = 0; @@ -867,160 +676,158 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, magField = getMagneticField(bc.timestamp()); } - // Process tracks in selected events int nTracksInEvent = 0; + for (const auto& track : tracks) { if (!track.has_collision()) continue; if (track.collisionId() != collId) continue; - // Fill phi cut monitoring before cut - if (applyPhiCut.value && track.pt() >= pTthresholdPhiCut.value) { - float phiPrime = getTransformedPhi(track.phi(), track.sign(), magField); - ue.fill(HIST("PhiCut/hPtVsPhiPrimeBefore"), track.pt(), phiPrime); - } - if (!passesTrackSelection(track, magField)) continue; - // Fill phi cut monitoring after cut - if (applyPhiCut.value && track.pt() >= pTthresholdPhiCut.value) { - float phiPrime = getTransformedPhi(track.phi(), track.sign(), magField); - ue.fill(HIST("PhiCut/hPtVsPhiPrimeAfter"), track.pt(), phiPrime); - } - nTracksInEvent++; - // Fill TPC cluster histograms ue.fill(HIST("hNclFoundTPC"), track.tpcNClsFound()); ue.fill(HIST("hNclPIDTPC"), track.tpcNClsPID()); - ue.fill(HIST("hNclFoundTPCvsPt"), track.pt(), track.tpcNClsFound()); - ue.fill(HIST("hNclPIDTPCvsPt"), track.pt(), track.tpcNClsPID()); - - // Inclusive histograms - ALWAYS fill for ALL tracks - ue.fill(HIST("Inclusive/hPtAllReco"), track.pt()); - ue.fill(HIST("Inclusive/hPtAllRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Inclusive/hPtMeasuredVsCent"), track.pt(), cent); ue.fill(HIST("hEta"), track.eta()); ue.fill(HIST("hPhi"), track.phi()); + ue.fill(HIST("Inclusive/hPtReco"), track.pt()); - // MC matching - FIXED: NO PID requirement for efficiency numerator if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); int pdgCode = std::abs(particle.pdgCode()); - if (particle.isPhysicalPrimary()) { - // Fill efficiency numerator for ALL primaries regardless of PID - ue.fill(HIST("Inclusive/hPtNumEff"), particle.pt()); - ue.fill(HIST("Inclusive/hPtNumEffVsCent"), particle.pt(), cent); - ue.fill(HIST("Inclusive/hPtPrimReco"), track.pt()); - ue.fill(HIST("Inclusive/hPtPrimRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Inclusive/hPtPrimRecoEv"), particle.pt()); - - // Per-species efficiency numerator - NO PID requirement! - if (pdgCode == PDG_t::kPiPlus) { - ue.fill(HIST("Pion/hPtNumEff"), particle.pt()); - ue.fill(HIST("Pion/hPtNumEffVsCent"), particle.pt(), cent); + if (pdgCode == PDG_t::kPiPlus) { + ue.fill(HIST("Pion/hPtReco"), track.pt()); + if (particle.isPhysicalPrimary()) { ue.fill(HIST("Pion/hPtPrimReco"), track.pt()); - ue.fill(HIST("Pion/hPtPrimRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Pion/hPtPrimRecoEv"), particle.pt()); - } else if (pdgCode == PDG_t::kKPlus) { - ue.fill(HIST("Kaon/hPtNumEff"), particle.pt()); - ue.fill(HIST("Kaon/hPtNumEffVsCent"), particle.pt(), cent); + ue.fill(HIST("Inclusive/hPtPrimReco"), track.pt()); + } else { + ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); + } + } else if (pdgCode == PDG_t::kKPlus) { + ue.fill(HIST("Kaon/hPtReco"), track.pt()); + if (particle.isPhysicalPrimary()) { ue.fill(HIST("Kaon/hPtPrimReco"), track.pt()); - ue.fill(HIST("Kaon/hPtPrimRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Kaon/hPtPrimRecoEv"), particle.pt()); - } else if (pdgCode == PDG_t::kProton) { - ue.fill(HIST("Proton/hPtNumEff"), particle.pt()); - ue.fill(HIST("Proton/hPtNumEffVsCent"), particle.pt(), cent); + ue.fill(HIST("Inclusive/hPtPrimReco"), track.pt()); + } else { + ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); + } + } else if (pdgCode == PDG_t::kProton) { + ue.fill(HIST("Proton/hPtReco"), track.pt()); + if (particle.isPhysicalPrimary()) { ue.fill(HIST("Proton/hPtPrimReco"), track.pt()); - ue.fill(HIST("Proton/hPtPrimRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Proton/hPtPrimRecoEv"), particle.pt()); + ue.fill(HIST("Inclusive/hPtPrimReco"), track.pt()); + } else { + ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); } - } else { - // Secondaries (non-primary particles) - ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); - if (pdgCode == PDG_t::kPiPlus) - ue.fill(HIST("Pion/hPtSecReco"), track.pt()); - else if (pdgCode == PDG_t::kKPlus) - ue.fill(HIST("Kaon/hPtSecReco"), track.pt()); - else if (pdgCode == PDG_t::kProton) - ue.fill(HIST("Proton/hPtSecReco"), track.pt()); } } - // PID selection - fill denominator for primary fraction and measured spectra - int bestSpecies = getBestPIDHypothesis(track); - if (bestSpecies == PartPion) { - // Denominator for pion primary fraction: all tracks identified as pions - ue.fill(HIST("Pion/hPtAllReco"), track.pt()); - ue.fill(HIST("Pion/hPtAllRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Pion/hPtMeasuredVsCent"), track.pt(), cent); - if (enablePIDHistograms) + if (passesPIDSelection(track, PartPion)) { + ue.fill(HIST("Pion/hPtMeasured"), track.pt()); + if (enablePIDHistograms) { ue.fill(HIST("Pion/hNsigmaTPC"), track.pt(), track.tpcNSigmaPi()); - } else if (bestSpecies == PartKaon) { - // Denominator for kaon primary fraction: all tracks identified as kaons - ue.fill(HIST("Kaon/hPtAllReco"), track.pt()); - ue.fill(HIST("Kaon/hPtAllRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Kaon/hPtMeasuredVsCent"), track.pt(), cent); - if (enablePIDHistograms) + } + } + + if (passesPIDSelection(track, PartKaon)) { + ue.fill(HIST("Kaon/hPtMeasured"), track.pt()); + if (enablePIDHistograms) { ue.fill(HIST("Kaon/hNsigmaTPC"), track.pt(), track.tpcNSigmaKa()); - } else if (bestSpecies == PartProton) { - // Denominator for proton primary fraction: all tracks identified as protons - ue.fill(HIST("Proton/hPtAllReco"), track.pt()); - ue.fill(HIST("Proton/hPtAllRecoVsCent"), track.pt(), cent); - ue.fill(HIST("Proton/hPtMeasuredVsCent"), track.pt(), cent); - if (enablePIDHistograms) + } + } + + if (passesPIDSelection(track, PartProton)) { + ue.fill(HIST("Proton/hPtMeasured"), track.pt()); + if (enablePIDHistograms) { ue.fill(HIST("Proton/hNsigmaTPC"), track.pt(), track.tpcNSigmaPr()); + } } } - // Fill event-level track multiplicity ue.fill(HIST("Centrality/hRecoMultVsCent"), cent, nTracksInEvent); } - LOG(info) << "\n=== CUT STATISTICS ==="; - LOG(info) << "Total collisions with valid info: " << centAll.size(); - LOG(info) << "Pass vertex cut: " << nPassVertex << " (" - << (centAll.size() > 0 ? 100.0 * nPassVertex / centAll.size() : 0.0) << "%)"; - LOG(info) << "Pass INEL cut: " << nPassINEL << " (" - << (centAll.size() > 0 ? 100.0 * nPassINEL / centAll.size() : 0.0) << "%)"; - LOG(info) << "Pass all cuts: " << nPassAll << " (" - << (centAll.size() > 0 ? 100.0 * nPassAll / centAll.size() : 0.0) << "%)"; - LOG(info) << "Reco collisions: " << nRecoCollisions; - LOG(info) << "Selected Events: " << nSelectedEvents; - LOG(info) << "Rejected Events: " << nRejectedEvents; - LOG(info) << "No Match: " << nNoMCMatch; - LOG(info) << "No Cent: " << nNoCent; - LOG(info) << "Invalid Cent: " << nInvalidCent; - - // Calculate mean centrality at each stage - if (!centAll.empty()) { - float meanAll = std::accumulate(centAll.begin(), centAll.end(), 0.0) / centAll.size(); - float meanVertex = centVertex.empty() ? 0 : std::accumulate(centVertex.begin(), centVertex.end(), 0.0) / centVertex.size(); - float meanINEL = centINEL.empty() ? 0 : std::accumulate(centINEL.begin(), centINEL.end(), 0.0) / centINEL.size(); - float meanSelected = centSelected.empty() ? 0 : std::accumulate(centSelected.begin(), centSelected.end(), 0.0) / centSelected.size(); - - LOG(info) << "\n=== CENTRALITY MEANS ==="; - LOG(info) << "Mean centrality (all): " << meanAll; - LOG(info) << "Mean centrality (after vertex): " << meanVertex; - LOG(info) << "Mean centrality (after INEL): " << meanINEL; - LOG(info) << "Mean centrality (selected): " << meanSelected; + for (const auto& mcCollision : mcCollisions) { + int64_t mcCollId = mcCollision.globalIndex(); + + if (reconstructedMCCollisions.find(mcCollId) == reconstructedMCCollisions.end()) + continue; + + auto centIt = mcCollToCentFromReco.find(mcCollId); + if (centIt == mcCollToCentFromReco.end()) + continue; + + auto nchIt = mcCollisionToNch.find(mcCollId); + int nGenCharged = (nchIt != mcCollisionToNch.end()) ? nchIt->second : 0; + + auto particlesInCollision = particles.sliceBy(perMCCol, mcCollId); + + for (const auto& particle : particlesInCollision) { + auto pdgParticle = pdg->GetParticle(particle.pdgCode()); + if (!pdgParticle || pdgParticle->Charge() == 0.) + continue; + if (!particle.isPhysicalPrimary()) + continue; + + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) + continue; + if (particle.pt() < cfgTrkLowPtCut.value) + continue; + + int pdgCode = std::abs(particle.pdgCode()); + float pt = particle.pt(); + + if (pdgCode == PDG_t::kPiPlus) { + ue.fill(HIST("Pion/hPtGenRecoEvent"), pt); + ue.fill(HIST("Pion/hPtGenRecoEvent_Mult"), pt, nGenCharged); + } else if (pdgCode == PDG_t::kKPlus) { + ue.fill(HIST("Kaon/hPtGenRecoEvent"), pt); + ue.fill(HIST("Kaon/hPtGenRecoEvent_Mult"), pt, nGenCharged); + } else if (pdgCode == PDG_t::kProton) { + ue.fill(HIST("Proton/hPtGenRecoEvent"), pt); + ue.fill(HIST("Proton/hPtGenRecoEvent_Mult"), pt, nGenCharged); + } + } + } + + for (const auto& mcCollId : inel0MCCollisions) { + auto centIt = mcCollToCentFromReco.find(mcCollId); + if (centIt != mcCollToCentFromReco.end()) { + float cent = centIt->second; + int nGenCharged = mcCollisionToNch[mcCollId]; + ue.fill(HIST("MC/EventLoss/GenMultVsCent"), cent, nGenCharged); + ue.fill(HIST("hEventsINEL_Cent"), cent); + } + } + + for (const auto& mcCollId : reconstructedMCCollisions) { + auto centIt = mcCollToCentFromReco.find(mcCollId); + if (centIt != mcCollToCentFromReco.end()) { + ue.fill(HIST("hEventsReco_Cent"), centIt->second); + } } ue.fill(HIST("hEventLossBreakdown"), 1.f, physicsSelectedMCCollisions.size()); ue.fill(HIST("hEventLossBreakdown"), 2.f, reconstructedMCCollisions.size()); ue.fill(HIST("hEventLossBreakdown"), 3.f, selectedMCCollisions.size()); - float efficiency = physicsSelectedMCCollisions.size() > 0 ? 100.f * selectedMCCollisions.size() / physicsSelectedMCCollisions.size() : 0; - ue.fill(HIST("hEventLossBreakdown"), 4.f, efficiency); - - LOG(info) << "\n=== FINAL EFFICIENCY ==="; - LOG(info) << "Physics selected: " << physicsSelectedMCCollisions.size(); - LOG(info) << "Reconstructed: " << reconstructedMCCollisions.size(); - LOG(info) << "Selected: " << selectedMCCollisions.size(); - LOG(info) << "Efficiency: " << efficiency << "%"; - LOG(info) << "=== processMC END ==="; + + ue.fill(HIST("MC/GenRecoCollisions"), 4.f, reconstructedMCCollisions.size()); + ue.fill(HIST("MC/GenRecoCollisions"), 5.f, selectedMCCollisions.size()); + + LOG(info) << "\n=== EVENT STATISTICS ==="; + LOG(info) << "INEL>0 MC collisions: " << inel0MCCollisions.size(); + LOG(info) << "Physics selected MC collisions: " << physicsSelectedMCCollisions.size(); + LOG(info) << "Reconstructed collisions: " << reconstructedMCCollisions.size(); + LOG(info) << "Selected collisions: " << selectedMCCollisions.size(); + + if (physicsSelectedMCCollisions.size() > 0) { + float efficiency = 100.f * selectedMCCollisions.size() / physicsSelectedMCCollisions.size(); + LOG(info) << "Final efficiency: " << efficiency << "%"; + } } void MultiplicityPt::processData(CollisionTableData::iterator const& /*collision*/,