diff --git a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx index 1a71c04db8b..6d9e1ca4fdf 100644 --- a/PWGLF/Tasks/Nuspex/multiplicityPt.cxx +++ b/PWGLF/Tasks/Nuspex/multiplicityPt.cxx @@ -11,7 +11,7 @@ /// \author Dushmanta Sahu (dushmanta.sahu@cern.ch) /// \file multiplicityPt.cxx -/// \brief Analysis to do PID with MC +/// \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 @@ -50,7 +50,7 @@ #include #include #include -#include // For std::accumulate +#include #include #include #include @@ -66,7 +66,6 @@ using BCsRun3 = soa::Join pdg; Service ccdb; static constexpr int CentBinMax = 100; @@ -79,7 +78,6 @@ struct MultiplicityPt { INEL = 0, INELgt0 = 1, INELgt1 = 2 - }; Configurable isRun3{"isRun3", true, "is Run3 dataset"}; @@ -99,7 +97,6 @@ struct MultiplicityPt { Configurable requireTrdOnly{"requireTrdOnly", false, "Require only tracks from TRD"}; Configurable requireNoTrd{"requireNoTrd", false, "Require tracks without TRD"}; - // Analysis switches Configurable enableDCAHistograms{"enableDCAHistograms", false, "Enable DCA histograms"}; Configurable enablePIDHistograms{"enablePIDHistograms", true, "Enable PID histograms"}; Configurable useCustomTrackCuts{"useCustomTrackCuts", true, "Flag to use custom track cuts"}; @@ -119,7 +116,6 @@ struct MultiplicityPt { Configurable nClTPCFoundCut{"nClTPCFoundCut", false, "Apply TPC found clusters cut"}; Configurable nClTPCPIDCut{"nClTPCPIDCut", true, "Apply TPC clusters for PID cut"}; - // Phi cut parameters Configurable applyPhiCut{"applyPhiCut", false, "Apply phi sector cut to remove problematic TPC regions"}; Configurable pTthresholdPhiCut{"pTthresholdPhiCut", 2.0f, "pT threshold above which to apply phi cut"}; Configurable phiCutLowParam1{"phiCutLowParam1", 0.119297, "First parameter for low phi cut"}; @@ -127,58 +123,37 @@ struct MultiplicityPt { Configurable phiCutHighParam1{"phiCutHighParam1", 0.16685, "First parameter for high phi cut"}; Configurable phiCutHighParam2{"phiCutHighParam2", 0.00981942, "Second parameter for high phi cut"}; - // Basic track cuts Configurable cfgTrkEtaCut{"cfgTrkEtaCut", 0.8f, "Eta range for tracks"}; Configurable cfgTrkLowPtCut{"cfgTrkLowPtCut", 0.15f, "Minimum constituent pT"}; - // PID selection - make them configurable per particle Configurable cfgCutNsigmaPi{"cfgCutNsigmaPi", 3.0f, "nsigma cut for pions"}; Configurable cfgCutNsigmaKa{"cfgCutNsigmaKa", 2.5f, "nsigma cut for kaons"}; Configurable cfgCutNsigmaPr{"cfgCutNsigmaPr", 2.5f, "nsigma cut for protons"}; - // Custom track cuts matching spectraTOF TrackSelection customTrackCuts; - // TF1 pointers for phi cuts TF1* fphiCutLow = nullptr; TF1* fphiCutHigh = nullptr; - // Histogram Registry HistogramRegistry ue; - // Data collisions (not used but kept for completeness) using CollisionTableData = soa::Join; - - // Track tables using TrackTableData = soa::Join; using TrackTableMC = soa::Join; - - // MC particles table using ParticlesMC = aod::McParticles; - - // MC collisions table using McCollisions = aod::McCollisions; - - // Reconstructed collisions (without joins that cause size mismatch) using RecoCollisions = aod::Collisions; - // Preslice for MC particles Preslice perMCCol = aod::mcparticle::mcCollisionId; enum ParticleSpecies : int { - kPion = 0, - kKaon = 1, - kProton = 2, - kNSpecies = 3 + PartPion = 0, + PartKaon = 1, + PartProton = 2, }; - static constexpr int PDGPion = kPiPlus; - static constexpr int PDGKaon = kKPlus; - static constexpr int PDGProton = kProton; - - // Get magnetic field from CCDB int getMagneticField(uint64_t timestamp) { static o2::parameters::GRPMagField* grpo = nullptr; @@ -193,7 +168,6 @@ struct MultiplicityPt { return grpo->getNominalL3Field(); } - // Get transformed phi for phi cut (with magnetic field) float getTransformedPhi(const float phi, const int charge, const float magField) const { float transformedPhi = phi; @@ -208,36 +182,28 @@ struct MultiplicityPt { return transformedPhi; } - // Phi cut function (with magnetic field) template bool passedPhiCut(const TrackType& track, float magField) const { - if (!applyPhiCut.value) { + if (!applyPhiCut.value) return true; - } - - if (track.pt() < pTthresholdPhiCut.value) { + if (track.pt() < pTthresholdPhiCut.value) return true; - } float pt = track.pt(); float phi = track.phi(); int charge = track.sign(); - if (magField < 0) { + if (magField < 0) phi = o2::constants::math::TwoPI - phi; - } - if (charge < 0) { + 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(pt) && phi > fphiCutLow->Eval(pt)) return false; - } - return true; } @@ -249,16 +215,12 @@ struct MultiplicityPt { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); if (!pdgParticle || pdgParticle->Charge() == 0.) continue; - if (!particle.isPhysicalPrimary()) continue; - if (std::abs(particle.eta()) > etaMax) continue; - if (particle.pt() < ptMin) continue; - count++; } return count; @@ -267,17 +229,13 @@ struct MultiplicityPt { template bool passedNClTPCFoundCut(const T& trk) const { - if (!nClTPCFoundCut.value) - return true; - return trk.tpcNClsFound() >= minTPCNClsFound.value; + return !nClTPCFoundCut.value || trk.tpcNClsFound() >= minTPCNClsFound.value; } template bool passedNClTPCPIDCut(const T& trk) const { - if (!nClTPCPIDCut.value) - return true; - return trk.tpcNClsPID() >= minTPCNClsPID.value; + return !nClTPCPIDCut.value || trk.tpcNClsPID() >= minTPCNClsPID.value; } template @@ -286,12 +244,10 @@ struct MultiplicityPt { if (useCustomTrackCuts.value) { for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { if (i == static_cast(TrackSelection::TrackCuts::kDCAxy) || - i == static_cast(TrackSelection::TrackCuts::kDCAz)) { + i == static_cast(TrackSelection::TrackCuts::kDCAz)) continue; - } - if (!customTrackCuts.IsSelected(track, static_cast(i))) { + if (!customTrackCuts.IsSelected(track, static_cast(i))) return false; - } } return true; } @@ -302,9 +258,8 @@ struct MultiplicityPt { bool passesDCAxyCut(TrackType const& track) const { if (useCustomTrackCuts.value) { - if (!passesCutWoDCA(track)) { + if (!passesCutWoDCA(track)) return false; - } constexpr float DcaXYConst = 0.0105f; constexpr float DcaXYPtScale = 0.0350f; constexpr float DcaXYPtPower = 1.1f; @@ -319,26 +274,18 @@ struct MultiplicityPt { { if (track.eta() < cfgCutEtaMin.value || track.eta() > cfgCutEtaMax.value) return false; - if (track.tpcChi2NCl() < minChi2PerClusterTPC.value || track.tpcChi2NCl() > maxChi2PerClusterTPC.value) return false; - if (!passesCutWoDCA(track)) return false; - if (!passesDCAxyCut(track)) return false; - if (!passedNClTPCFoundCut(track)) return false; - if (!passedNClTPCPIDCut(track)) return false; - - // Add phi cut with magnetic field if (!passedPhiCut(track, magField)) return false; - return true; } @@ -346,21 +293,19 @@ struct MultiplicityPt { bool passesPIDSelection(TrackType const& track) const { float nsigmaTPC = 0.f; - - if constexpr (species == kPion) { + if constexpr (species == PartPion) nsigmaTPC = track.tpcNSigmaPi(); - } else if constexpr (species == kKaon) { + else if constexpr (species == PartKaon) nsigmaTPC = track.tpcNSigmaKa(); - } else if constexpr (species == kProton) { + else if constexpr (species == PartProton) nsigmaTPC = track.tpcNSigmaPr(); - } float cutValue = cfgCutNsigma.value; - if constexpr (species == kPion) + if constexpr (species == PartPion) cutValue = cfgCutNsigmaPi.value; - if constexpr (species == kKaon) + if constexpr (species == PartKaon) cutValue = cfgCutNsigmaKa.value; - if constexpr (species == kProton) + if constexpr (species == PartProton) cutValue = cfgCutNsigmaPr.value; return (std::abs(nsigmaTPC) < cutValue); @@ -378,17 +323,16 @@ struct MultiplicityPt { if (nsigmaPi < cfgCutNsigmaPi.value && nsigmaPi < minNSigma) { minNSigma = nsigmaPi; - bestSpecies = kPion; + bestSpecies = PartPion; } if (nsigmaKa < cfgCutNsigmaKa.value && nsigmaKa < minNSigma) { minNSigma = nsigmaKa; - bestSpecies = kKaon; + bestSpecies = PartKaon; } if (nsigmaPr < cfgCutNsigmaPr.value && nsigmaPr < minNSigma) { minNSigma = nsigmaPr; - bestSpecies = kProton; + bestSpecies = PartProton; } - return bestSpecies; } @@ -398,44 +342,24 @@ struct MultiplicityPt { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); if (!pdgParticle || pdgParticle->Charge() == 0.) return false; - if (!particle.isPhysicalPrimary()) return false; - if (std::abs(particle.eta()) >= cfgCutEtaMax.value) return false; if (particle.pt() < cfgTrkLowPtCut.value) return false; - return true; } - void processData(CollisionTableData::iterator const& collision, - TrackTableData const& tracks, - BCsRun3 const& bcs); + void processData(CollisionTableData::iterator const& collision, TrackTableData const& tracks, BCsRun3 const& bcs); PROCESS_SWITCH(MultiplicityPt, processData, "process data", false); - void processMC(TrackTableMC const& tracks, - aod::McParticles const& particles, - aod::McCollisions const& mcCollisions, - RecoCollisions const& collisions, - aod::McCollisionLabels const& labels, - aod::McCentFT0Ms const& centTable, - BCsRun3 const& bcs); + void processMC(TrackTableMC const& tracks, aod::McParticles const& particles, aod::McCollisions const& mcCollisions, + RecoCollisions const& collisions, aod::McCollisionLabels const& labels, aod::McCentFT0Ms const& centTable, BCsRun3 const& bcs); PROCESS_SWITCH(MultiplicityPt, processMC, "process MC", true); void init(InitContext const&); - - void endOfStream(EndOfStreamContext& /*eos*/) - { - LOG(info) << "\n=== END OF STREAM: Writing histograms to output ==="; - auto hGenMult = ue.get(HIST("MC/EventLoss/GenMultVsCent")); - if (hGenMult) { - LOG(info) << "GenMultVsCent: Entries=" << hGenMult->GetEntries() - << ", Integral=" << hGenMult->Integral(); - } - LOG(info) << "=== END OF STREAM COMPLETE ==="; - } + void endOfStream(EndOfStreamContext& /*eos*/) {} }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) @@ -446,32 +370,16 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) void MultiplicityPt::init(InitContext const&) { LOG(info) << "=================================================="; - LOG(info) << "Initializing MultiplicityPt task with full centrality diagnostics"; + LOG(info) << "Initializing MultiplicityPt task - FULL CORRECTION FACTORS"; LOG(info) << "=================================================="; - // Initialize phi cut functions if (applyPhiCut.value) { - fphiCutLow = new TF1("StandardPhiCutLow", - Form("%f/x/x+pi/18.0-%f", - phiCutLowParam1.value, phiCutLowParam2.value), - 0, 50); - fphiCutHigh = new TF1("StandardPhiCutHigh", - Form("%f/x+pi/18.0+%f", - phiCutHighParam1.value, phiCutHighParam2.value), - 0, 50); - - LOGF(info, "=== Phi Cut Parameters ==="); - LOGF(info, "Low cut: %.6f/x² + pi/18 - %.6f", - phiCutLowParam1.value, phiCutLowParam2.value); - LOGF(info, "High cut: %.6f/x + pi/18 + %.6f", - phiCutHighParam1.value, phiCutHighParam2.value); - LOGF(info, "Applied for pT > %.1f GeV/c", pTthresholdPhiCut.value); + fphiCutLow = new TF1("StandardPhiCutLow", Form("%f/x/x+pi/18.0-%f", phiCutLowParam1.value, phiCutLowParam2.value), 0, 50); + fphiCutHigh = new TF1("StandardPhiCutHigh", Form("%f/x+pi/18.0+%f", phiCutHighParam1.value, phiCutHighParam2.value), 0, 50); } if (useCustomTrackCuts.value) { - LOG(info) << "Using custom track cuts matching spectraTOF approach"; customTrackCuts = getGlobalTrackSelectionRun3ITSMatch(itsPattern.value); - customTrackCuts.SetRequireITSRefit(requireITS.value); customTrackCuts.SetRequireTPCRefit(requireTPC.value); customTrackCuts.SetRequireGoldenChi2(requireGoldenChi2.value); @@ -482,80 +390,50 @@ void MultiplicityPt::init(InitContext const&) customTrackCuts.SetMinNCrossedRowsOverFindableClustersTPC(minNCrossedRowsOverFindableClustersTPC.value); customTrackCuts.SetMaxDcaXYPtDep([](float /*pt*/) { return 10000.f; }); customTrackCuts.SetMaxDcaZ(maxDcaZ.value); - 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.}; - - // Fine centrality binning for diagnostics (100 bins, guaranteed increasing) std::vector centBinningFine; - for (int i = 0; i <= CentBinMax; i++) { + for (int i = 0; i <= CentBinMax; i++) centBinningFine.push_back(static_cast(i)); - } AxisSpec centAxis = {centBinningStd, "FT0M Centrality (%)"}; AxisSpec centFineAxis = {centBinningFine, "FT0M Centrality (%)"}; - // Multiplicity axes - properly defined std::vector multBins; - for (int i = 0; i <= MultBinMax; i++) { + for (int i = 0; i <= MultBinMax; i++) multBins.push_back(static_cast(i)); - } AxisSpec multAxis = {multBins, "N_{ch}^{gen} (|#eta|<0.8)"}; - // Reconstructed multiplicity axis - properly defined with explicit bin edges std::vector recoMultBins; - for (int i = 0; i <= RecMultBinMax; i++) { + for (int i = 0; i <= RecMultBinMax; i++) recoMultBins.push_back(static_cast(i)); - } AxisSpec recoMultAxis = {recoMultBins, "N_{ch}^{reco}"}; - // Centrality diagnostic histograms - USE FINE BINNING - ue.add("Centrality/hCentRaw", "Raw FT0M Centrality (no cuts);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}); - - // 2D correlations - USE FINE BINNING FOR DIAGNOSTICS - 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}); - - // Vertex resolution vs centrality - ue.add("Centrality/hVertexResVsCent", "Vertex Resolution vs Centrality;Centrality (%);V_{z} resolution (cm)", - HistType::kTH2D, {centFineAxis, {100, -1, 1}}); - - // INEL class distributions - ue.add("INEL/hINELClass", "INEL Class for MC Collisions;INEL Class;Counts", - HistType::kTH1D, {{3, 0.5, 3.5}}); + ue.add("Centrality/hCentRaw", "Raw FT0M Centrality (no cuts);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("INEL/hINELVsCent", "INEL Class vs Centrality;Centrality (%);INEL Class", HistType::kTH2D, {centFineAxis, {3, 0.5, 3.5}}); - // Cut flow - 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, {{6, 0.5, 6.5}}); auto hCut = ue.get(HIST("CutFlow/hCutStats")); hCut->GetXaxis()->SetBinLabel(1, "All reco events"); hCut->GetXaxis()->SetBinLabel(2, "Has MC match"); @@ -564,180 +442,136 @@ void MultiplicityPt::init(InitContext const&) 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}); + ue.add("CutFlow/hCentPerCut", "Centrality Distribution at Each Cut;Cut Stage;Centrality (%)", HistType::kTH2D, {{6, 0.5, 6.5}, centFineAxis}); - ue.add("MC/GenRecoCollisions", "Generated and Reconstructed MC Collisions", - HistType::kTH1D, {{10, 0.5, 10.5}}); + 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, {{4, 0.5, 4.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"); - // Multiplicity histograms - ue.add("MC/EventLoss/NchGenerated", "Generated charged multiplicity;N_{ch}^{gen} (|#eta|<0.8);Counts", - HistType::kTH1D, {{200, 0, 200}}); - ue.add("MC/EventLoss/NchGenerated_PhysicsSelected", "Generated charged multiplicity (physics selected);N_{ch}^{gen} (|#eta|<0.8);Counts", - HistType::kTH1D, {{200, 0, 200}}); - ue.add("MC/EventLoss/NchGenerated_Reconstructed", "Generated charged multiplicity (reconstructed);N_{ch}^{gen} (|#eta|<0.8);Counts", - HistType::kTH1D, {{200, 0, 200}}); - - // pT vs Multiplicity - 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}}); - - // Centrality vs Multiplicity correlations - USE STANDARD BINNING FOR THESE - ue.add("MC/EventLoss/GenMultVsCent", "Generated charged particles vs FT0M centrality;FT0M Centrality (%);N_{ch}^{gen} (|#eta|<0.8)", - 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}); - - // TPC cluster histograms - ue.add("hNclFoundTPC", "Number of TPC found clusters", - HistType::kTH1D, {{200, 0, 200, "N_{cl, found}"}}); - ue.add("hNclPIDTPC", "Number of TPC PID clusters", - HistType::kTH1D, {{200, 0, 200, "N_{cl, PID}"}}); - 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.}}); - - // Inclusive histograms - ue.add("Inclusive/hPtPrimGenAll", "All generated primaries (no cuts);#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/hPtPrimGen", "Generated primaries (after physics selection);#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}); - - ue.add("Inclusive/hPtNumEff", "Tracking efficiency numerator;#it{p}_{T} (GeV/#it{c});Counts", - HistType::kTH1D, {ptAxis}); - ue.add("Inclusive/hPtDenEff", "Tracking efficiency denominator;#it{p}_{T} (GeV/#it{c});Counts", - HistType::kTH1D, {ptAxis}); - - 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("Inclusive/hPtMeasuredVsCent", "All measured tracks (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%)", - HistType::kTH2D, {ptAxis, centAxis}); - - // Phi cut monitoring histograms - 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}}); - } - - // Particle-specific histograms - const std::array particleNames = {"Pion", "Kaon", "Proton"}; - const std::array particleSymbols = {"#pi^{#pm}", "K^{#pm}", "p+#bar{p}"}; - - for (int iSpecies = 0; iSpecies < kNSpecies; ++iSpecies) { - const auto& name = particleNames[iSpecies]; - const auto& symbol = particleSymbols[iSpecies]; - - ue.add(Form("%s/hPtPrimGenAll", name.c_str()), - Form("All generated %s (no cuts);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtPrimBadVertex", name.c_str()), - Form("Generated %s (bad vertex);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtPrimGen", name.c_str()), - Form("Generated %s (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtPrimRecoEv", name.c_str()), - Form("Generated %s (reco events);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtPrimGoodEv", name.c_str()), - Form("Generated %s (good events);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - - ue.add(Form("%s/hPtNumEff", name.c_str()), - Form("%s tracking efficiency numerator;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtDenEff", name.c_str()), - Form("%s tracking efficiency denominator;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - - ue.add(Form("%s/hPtAllReco", name.c_str()), - Form("All reconstructed %s;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtPrimReco", name.c_str()), - Form("Reconstructed primary %s;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - ue.add(Form("%s/hPtSecReco", name.c_str()), - Form("Reconstructed secondary %s;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), - HistType::kTH1D, {ptAxis}); - - ue.add(Form("%s/hPtMeasuredVsCent", name.c_str()), - Form("Measured %s (PID) vs centrality;#it{p}_{T} (GeV/#it{c});FT0M Centrality (%%)", symbol.c_str()), - HistType::kTH2D, {ptAxis, centAxis}); - - if (enablePIDHistograms) { - ue.add(Form("%s/hNsigmaTPC", name.c_str()), - Form("TPC n#sigma %s;#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", symbol.c_str()), - HistType::kTH2D, {ptAxis, {200, -10, 10}}); - } - } - - // Event selection histogram - constexpr int NEvSelBins = 20; - constexpr float EvSelMin = 0.5f; - constexpr float EvSelMax = 20.5f; - ue.add("evsel", "Event selection", HistType::kTH1D, {{NEvSelBins, EvSelMin, EvSelMax}}); - auto h = ue.get(HIST("evsel")); - h->GetXaxis()->SetBinLabel(1, "Events read"); - h->GetXaxis()->SetBinLabel(4, "Trigger passed"); - h->GetXaxis()->SetBinLabel(5, "NoITSROFrameBorder"); - h->GetXaxis()->SetBinLabel(6, "NoSameBunchPileup"); - h->GetXaxis()->SetBinLabel(7, "IsGoodZvtxFT0vsPV"); - h->GetXaxis()->SetBinLabel(8, "IsVertexITSTPC"); - h->GetXaxis()->SetBinLabel(9, "NoTimeFrameBorder"); - h->GetXaxis()->SetBinLabel(13, "posZ passed"); - - // Basic tracking histograms + 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/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/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}); + + // 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/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) + 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/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/hNsigmaTPC", "TPC n#sigma p+#bar{p};#it{p}_{T} (GeV/#it{c});n#sigma_{TPC}", HistType::kTH2D, {ptAxis, {200, -10, 10}}); + + 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}}); - LOG(info) << "=== Initialized MultiplicityPt task with full centrality diagnostics ==="; - LOG(info) << "Standard centrality binning: " << centBinningStd.size() - 1 << " bins (0-100%)"; - LOG(info) << "Fine centrality binning: " << centBinningFine.size() - 1 << " bins (0-100%)"; + 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) { - LOG(info) << "Phi cut ENABLED for pT > " << pTthresholdPhiCut.value << " GeV/c"; + 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}}); } -} -void MultiplicityPt::processData(CollisionTableData::iterator const& /*collision*/, - TrackTableData const& /*tracks*/, - BCsRun3 const& /*bcs*/) -{ - // Intentionally empty - data processing disabled + LOG(info) << "=== Initialization complete - All correction factor histograms defined ==="; } void MultiplicityPt::processMC(TrackTableMC const& tracks, @@ -749,228 +583,203 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, BCsRun3 const& /*bcs*/) { LOG(info) << "\n=== processMC START ==="; - LOG(info) << "Total MC collisions (generated): " << mcCollisions.size(); - LOG(info) << "Total reconstructed collisions: " << collisions.size(); - LOG(info) << "Total collision labels: " << labels.size(); - LOG(info) << "Total centrality entries: " << centTable.size(); - - LOG(info) << "\n=== CENTRALITY DEBUG - RAW DATA ==="; - LOG(info) << "First 20 centrality values from centTable:"; - int debugCount = 0; - float minCent = 999.0f, maxCent = -999.0f; - std::map centDistribution; - - for (const auto& cent : centTable) { - float c = cent.centFT0M(); - if (debugCount < DebugCountMax) { - LOG(info) << " Cent entry " << debugCount << ": " << c; - } - minCent = std::min(minCent, c); - maxCent = std::max(maxCent, c); - - int bin10 = static_cast(c / 10) * 10; - centDistribution[bin10]++; - debugCount++; - } - - LOG(info) << "Centrality range: [" << minCent << ", " << maxCent << "]"; - LOG(info) << "Distribution by 10% bins:"; - for (int i = 0; i < CentBinMax; i += 10) { - LOG(info) << " " << i << "-" << i + 10 << "%: " << centDistribution[i]; - } - - // Check if centrality is inverted (0 = peripheral, 100 = central) - // If minCent is near 0 and maxCent near 100, check correlation with multiplicity - LOG(info) << "Checking if centrality might be inverted..."; - LOG(info) << "Will check correlation with multiplicity in the next step."; std::map mcCollisionToNch; std::map mcCollisionVz; std::set physicsSelectedMCCollisions; - std::map mcCollisionToINELClass; // 0=INEL0, 1=INEL>0, 2=INEL>1 + std::map mcCollisionToINELClass; ue.fill(HIST("MC/GenRecoCollisions"), 1.f, mcCollisions.size()); - ue.fill(HIST("MC/GenRecoCollisions"), 2.f, collisions.size()); - - LOG(info) << "\n--- FIRST PASS: Building MC collision maps ---"; - - int mcWithParticles = 0; - int mcINELgt0 = 0, mcINELgt1 = 0; for (const auto& mcCollision : mcCollisions) { int64_t mcCollId = mcCollision.globalIndex(); auto particlesInCollision = particles.sliceBy(perMCCol, mcCollId); int nGenCharged = countGeneratedChargedPrimaries(particlesInCollision, cfgCutEtaMax.value, cfgTrkLowPtCut.value); - mcCollisionToNch[mcCollId] = nGenCharged; mcCollisionVz[mcCollId] = mcCollision.posZ(); - // Determine INEL class bool inel0 = o2::pwglf::isINELgt0mc(particlesInCollision, pdg); bool inel1 = o2::pwglf::isINELgt1mc(particlesInCollision, pdg); - - int inelClass = 0; - if (inel1) - inelClass = 2; - else if (inel0) - inelClass = 1; + int inelClass = inel1 ? 2 : (inel0 ? 1 : 0); mcCollisionToINELClass[mcCollId] = inelClass; ue.fill(HIST("INEL/hINELClass"), inelClass); - - if (inel0) - mcINELgt0++; - if (inel1) - mcINELgt1++; - if (nGenCharged > 0) - mcWithParticles++; - ue.fill(HIST("MC/EventLoss/NchGenerated"), nGenCharged); - // Physics selection based on vertex and INEL cuts - bool physicsSelected = true; - - if (std::abs(mcCollision.posZ()) > cfgCutVertex.value) { - physicsSelected = false; - } - - // Apply INEL cut based on configuration - if (cfgINELCut.value == INELgt0 && !inel0) { + bool physicsSelected = (std::abs(mcCollision.posZ()) <= cfgCutVertex.value); + if (cfgINELCut.value == INELgt0 && !inel0) physicsSelected = false; - } - if (cfgINELCut.value == INELgt1 && !inel1) { + if (cfgINELCut.value == INELgt1 && !inel1) physicsSelected = false; - } if (physicsSelected) { physicsSelectedMCCollisions.insert(mcCollId); ue.fill(HIST("MC/EventLoss/NchGenerated_PhysicsSelected"), nGenCharged); - - if (inel0) { + if (inel0) ue.fill(HIST("MC/GenRecoCollisions"), 3.f); - } - if (inel1) { + 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."; + } } - LOG(info) << "\n--- FIRST PASS SUMMARY ---"; - LOG(info) << "Total MC collisions processed: " << mcCollisions.size(); - LOG(info) << "MC collisions with particles: " << mcWithParticles; - LOG(info) << "INEL0: " << (mcCollisions.size() - mcINELgt0); - LOG(info) << "INEL>0: " << mcINELgt0; - LOG(info) << "INEL>1: " << mcINELgt1; - LOG(info) << "Physics-selected MC collisions: " << physicsSelectedMCCollisions.size(); + 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); + + 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) + ue.fill(HIST("Pion/hPtPrimGenAll"), pt); + else if (pdgCode == PDG_t::kKPlus) + ue.fill(HIST("Kaon/hPtPrimGenAll"), pt); + else if (pdgCode == PDG_t::kProton) + ue.fill(HIST("Proton/hPtPrimGenAll"), pt); + + // 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); + } else if (pdgCode == PDG_t::kKPlus) { + ue.fill(HIST("Kaon/hPtPrimGen"), pt); + ue.fill(HIST("Kaon/hPtDenEff"), pt); + } else if (pdgCode == PDG_t::kProton) { + ue.fill(HIST("Proton/hPtPrimGen"), pt); + ue.fill(HIST("Proton/hPtDenEff"), pt); + } + } + + // 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); + } + } + } std::map recoToMcMap; std::map recoToCentMap; - size_t nCollisions = collisions.size(); - - // Associate labels with collisions by index size_t iLabel = 0; + size_t nCollisions = collisions.size(); for (const auto& label : labels) { if (iLabel < nCollisions) { const auto& collision = collisions.iteratorAt(iLabel); - int64_t recoCollId = collision.globalIndex(); - int64_t mcCollId = label.mcCollisionId(); - recoToMcMap[recoCollId] = mcCollId; + recoToMcMap[collision.globalIndex()] = label.mcCollisionId(); } iLabel++; } - // Associate centrality with collisions by index size_t iCent = 0; for (const auto& cent : centTable) { if (iCent < nCollisions) { const auto& collision = collisions.iteratorAt(iCent); - int64_t recoCollId = collision.globalIndex(); float centValue = cent.centFT0M(); - - // Fill raw centrality histogram + recoToCentMap[collision.globalIndex()] = centValue; ue.fill(HIST("Centrality/hCentRaw"), centValue); - - recoToCentMap[recoCollId] = centValue; } iCent++; } - LOG(info) << "\n--- MAP SIZES ---"; - LOG(info) << "recoToMcMap size: " << recoToMcMap.size(); - LOG(info) << "recoToCentMap size: " << recoToCentMap.size(); - - LOG(info) << "\n=== CENTRALITY VS MULTIPLICITY DEBUG ==="; - - // Create temporary vectors to check correlation - std::vector> centMultPairs; - for (const auto& collision : collisions) { - int64_t collId = collision.globalIndex(); - - auto mcIt = recoToMcMap.find(collId); - if (mcIt == recoToMcMap.end()) - continue; - - auto centIt = recoToCentMap.find(collId); - if (centIt == recoToCentMap.end()) - continue; - - auto nchIt = mcCollisionToNch.find(mcIt->second); - if (nchIt == mcCollisionToNch.end()) - continue; - - centMultPairs.push_back({centIt->second, nchIt->second}); - } - - // Sort by centrality - std::sort(centMultPairs.begin(), centMultPairs.end()); - - LOG(info) << "Correlation between centrality and multiplicity:"; - LOG(info) << " If centrality is normal (0=central, 100=peripheral), multiplicity should decrease with centrality"; - LOG(info) << " If inverted (0=peripheral, 100=central), multiplicity should increase with centrality"; - - // Print a few samples across the range - if (centMultPairs.size() > CentMultClasses) { - for (size_t i = 0; i < centMultPairs.size(); i += centMultPairs.size() / 10) { - LOG(info) << " Cent: " << centMultPairs[i].first - << "%, Mult: " << centMultPairs[i].second; - } - } - - //=========================================================================== - // SECOND PASS: Process reconstructed collisions with detailed cut accounting - //=========================================================================== - - LOG(info) << "\n--- SECOND PASS: Processing reconstructed collisions ---"; + 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; - - // Cut counters int nPassVertex = 0; int nPassINEL = 0; int nPassAll = 0; - // For mean calculations - std::vector centAll, centVertex, centINEL, centSelected; - for (const auto& collision : collisions) { nRecoCollisions++; - int64_t collId = collision.globalIndex(); - // Fill cut flow ue.fill(HIST("CutFlow/hCutStats"), 1); - // Get MC collision ID from labels map auto mcIt = recoToMcMap.find(collId); if (mcIt == recoToMcMap.end()) { nNoMCMatch++; @@ -979,27 +788,20 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, ue.fill(HIST("CutFlow/hCutStats"), 2); int64_t mcCollId = mcIt->second; - - // Get generated multiplicity for this MC collision auto nchIt = mcCollisionToNch.find(mcCollId); - if (nchIt == mcCollisionToNch.end()) { + if (nchIt == mcCollisionToNch.end()) continue; - } - int nGenCharged = nchIt->second; - // Get INEL class auto inelIt = mcCollisionToINELClass.find(mcCollId); int inelClass = (inelIt != mcCollisionToINELClass.end()) ? inelIt->second : 0; - // Get centrality from cent map auto centIt = recoToCentMap.find(collId); if (centIt == recoToCentMap.end()) { nNoCent++; continue; } ue.fill(HIST("CutFlow/hCutStats"), 3); - float cent = centIt->second; if (cent < 0 || cent > CentBinMax) { nInvalidCent++; @@ -1023,13 +825,11 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, nPassVertex++; } - // Check INEL selection at generator level 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); @@ -1038,15 +838,11 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, nPassINEL++; } - // Fill GenMultVsCent for all reconstructed events ue.fill(HIST("MC/EventLoss/GenMultVsCent"), cent, nGenCharged); ue.fill(HIST("MC/EventLoss/NchGenerated_Reconstructed"), nGenCharged); - reconstructedMCCollisions.insert(mcCollId); - // Apply all cuts bool passedAll = passVertex && passINEL; - if (!passedAll) { ue.fill(HIST("MC/EventLoss/GenMultVsCent_Rejected"), cent, nGenCharged); nRejectedEvents++; @@ -1085,9 +881,8 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, ue.fill(HIST("PhiCut/hPtVsPhiPrimeBefore"), track.pt(), phiPrime); } - if (!passesTrackSelection(track, magField)) { + if (!passesTrackSelection(track, magField)) continue; - } // Fill phi cut monitoring after cut if (applyPhiCut.value && track.pt() >= pTthresholdPhiCut.value) { @@ -1103,65 +898,81 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, 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()); + // 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()); - if (pdgCode == PDGPion) { + // 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); ue.fill(HIST("Pion/hPtPrimReco"), track.pt()); - } else if (pdgCode == PDGKaon) { + 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("Kaon/hPtPrimReco"), track.pt()); - } else if (pdgCode == PDGProton) { + 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("Proton/hPtPrimReco"), track.pt()); + ue.fill(HIST("Proton/hPtPrimRecoVsCent"), track.pt(), cent); + ue.fill(HIST("Proton/hPtPrimRecoEv"), particle.pt()); } } else { + // Secondaries (non-primary particles) ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); - - if (pdgCode == PDGPion) { + if (pdgCode == PDG_t::kPiPlus) ue.fill(HIST("Pion/hPtSecReco"), track.pt()); - } else if (pdgCode == PDGKaon) { + else if (pdgCode == PDG_t::kKPlus) ue.fill(HIST("Kaon/hPtSecReco"), track.pt()); - } else if (pdgCode == PDGProton) { + 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 == kPion) { - ue.fill(HIST("Pion/hPtMeasuredVsCent"), track.pt(), cent); + if (bestSpecies == PartPion) { + // Denominator for pion primary fraction: all tracks identified as pions ue.fill(HIST("Pion/hPtAllReco"), track.pt()); - - if (enablePIDHistograms) { + ue.fill(HIST("Pion/hPtAllRecoVsCent"), track.pt(), cent); + ue.fill(HIST("Pion/hPtMeasuredVsCent"), track.pt(), cent); + if (enablePIDHistograms) ue.fill(HIST("Pion/hNsigmaTPC"), track.pt(), track.tpcNSigmaPi()); - } - } else if (bestSpecies == kKaon) { - ue.fill(HIST("Kaon/hPtMeasuredVsCent"), track.pt(), cent); + } else if (bestSpecies == PartKaon) { + // Denominator for kaon primary fraction: all tracks identified as kaons ue.fill(HIST("Kaon/hPtAllReco"), track.pt()); - - if (enablePIDHistograms) { + ue.fill(HIST("Kaon/hPtAllRecoVsCent"), track.pt(), cent); + ue.fill(HIST("Kaon/hPtMeasuredVsCent"), track.pt(), cent); + if (enablePIDHistograms) ue.fill(HIST("Kaon/hNsigmaTPC"), track.pt(), track.tpcNSigmaKa()); - } - } else if (bestSpecies == kProton) { - ue.fill(HIST("Proton/hPtMeasuredVsCent"), track.pt(), cent); + } else if (bestSpecies == PartProton) { + // Denominator for proton primary fraction: all tracks identified as protons ue.fill(HIST("Proton/hPtAllReco"), track.pt()); - - if (enablePIDHistograms) { + ue.fill(HIST("Proton/hPtAllRecoVsCent"), track.pt(), cent); + ue.fill(HIST("Proton/hPtMeasuredVsCent"), track.pt(), cent); + if (enablePIDHistograms) ue.fill(HIST("Proton/hNsigmaTPC"), track.pt(), track.tpcNSigmaPr()); - } } } @@ -1169,7 +980,6 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Centrality/hRecoMultVsCent"), cent, nTracksInEvent); } - // Calculate and display cut statistics LOG(info) << "\n=== CUT STATISTICS ==="; LOG(info) << "Total collisions with valid info: " << centAll.size(); LOG(info) << "Pass vertex cut: " << nPassVertex << " (" @@ -1202,7 +1012,6 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, 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); @@ -1213,3 +1022,9 @@ void MultiplicityPt::processMC(TrackTableMC const& tracks, LOG(info) << "Efficiency: " << efficiency << "%"; LOG(info) << "=== processMC END ==="; } + +void MultiplicityPt::processData(CollisionTableData::iterator const& /*collision*/, + TrackTableData const& /*tracks*/, + BCsRun3 const& /*bcs*/) +{ +}