diff --git a/PWGLF/Tasks/Nuspex/multiplicitypt.cxx b/PWGLF/Tasks/Nuspex/multiplicitypt.cxx index d3d64510028..ca884a420da 100644 --- a/PWGLF/Tasks/Nuspex/multiplicitypt.cxx +++ b/PWGLF/Tasks/Nuspex/multiplicitypt.cxx @@ -1,12 +1,14 @@ +#include "PWGLF/Utils/inelGt.h" + #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/Centrality.h" #include "Common/DataModel/PIDResponse.h" -#include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponseTPC.h" +#include "Common/DataModel/TrackSelectionTables.h" #include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" @@ -16,16 +18,15 @@ #include "Framework/StaticFor.h" #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/Track.h" -#include "PWGLF/Utils/inelGt.h" #include "TDatabasePDG.h" #include #include #include #include -#include #include +#include #include using namespace o2; @@ -41,7 +42,7 @@ struct multiplicitypt { // Configurables - Matching spectraTOF approach Configurable isRun3{"isRun3", true, "is Run3 dataset"}; - + // Event selection configurables (spectraTOF style) 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"}; @@ -65,7 +66,7 @@ struct multiplicitypt { // Multiplicity estimator (like spectraTOF) Configurable multiplicityEstimator{"multiplicityEstimator", 6, - "Multiplicity estimator: 0=NoMult, 1=MultFV0M, 2=MultFT0M, 3=MultFDDM, 4=MultTracklets, 5=MultTPC, 6=MultNTracksPV, 7=MultNTracksPVeta1, 8=CentFT0C, 9=CentFT0M, 10=CentFV0A"}; + "Multiplicity estimator: 0=NoMult, 1=MultFV0M, 2=MultFT0M, 3=MultFDDM, 4=MultTracklets, 5=MultTPC, 6=MultNTracksPV, 7=MultNTracksPVeta1, 8=CentFT0C, 9=CentFT0M, 10=CentFV0A"}; // Analysis switches Configurable enableDCAHistograms{"enableDCAHistograms", false, "Enable DCA histograms"}; @@ -145,7 +146,7 @@ struct multiplicitypt { // ======================================================================== // PROCESS FUNCTION DECLARATIONS - SPECTRATOF STYLE // ======================================================================== - + // Data processing void processData(CollisionTableData::iterator const& collision, TrackTableData const& tracks); @@ -166,7 +167,7 @@ struct multiplicitypt { // ======================================================================== // TRACK SELECTION FUNCTIONS - MATCHING spectraTOF // ======================================================================== - + template bool passesCutWoDCA(TrackType const& track) const { @@ -219,12 +220,12 @@ struct multiplicitypt { // ======================================================================== // PID SELECTION FUNCTIONS - TPC ONLY (OLD NON-EXCLUSIVE METHOD) // ======================================================================== - + template bool passesPIDSelection(TrackType const& track) const { float nsigmaTPC = 0.f; - + if constexpr (species == kPion) { nsigmaTPC = track.tpcNSigmaPi(); } else if constexpr (species == kKaon) { @@ -232,7 +233,7 @@ struct multiplicitypt { } else if constexpr (species == kProton) { nsigmaTPC = track.tpcNSigmaPr(); } - + // TPC-only PID (works for all pT, but better at low pT < 1 GeV/c) return (std::abs(nsigmaTPC) < cfgCutNsigma.value); } @@ -240,20 +241,20 @@ struct multiplicitypt { // ======================================================================== // EXCLUSIVE PID SELECTION - Returns best hypothesis for a track // ======================================================================== - + template int getBestPIDHypothesis(TrackType const& track) const { // Return values: -1 = no ID, 0 = pion, 1 = kaon, 2 = proton - + float nsigmaPi = std::abs(track.tpcNSigmaPi()); float nsigmaKa = std::abs(track.tpcNSigmaKa()); float nsigmaPr = std::abs(track.tpcNSigmaPr()); - + // Find the hypothesis with smallest |nσ| that passes the cut float minNSigma = 999.0f; int bestSpecies = -1; - + if (nsigmaPi < cfgCutNsigma.value && nsigmaPi < minNSigma) { minNSigma = nsigmaPi; bestSpecies = kPion; @@ -266,56 +267,76 @@ struct multiplicitypt { minNSigma = nsigmaPr; bestSpecies = kProton; } - + return bestSpecies; } // ======================================================================== // EVENT SELECTION FUNCTION - EXACT spectraTOF // ======================================================================== - + template bool isEventSelected(CollisionType const& collision) { if constexpr (fillHistograms) { ue.fill(HIST("evsel"), 1.f); - if (collision.isInelGt0()) ue.fill(HIST("evsel"), 2.f); - if (collision.isInelGt1()) ue.fill(HIST("evsel"), 3.f); + if (collision.isInelGt0()) + ue.fill(HIST("evsel"), 2.f); + if (collision.isInelGt1()) + ue.fill(HIST("evsel"), 3.f); } if (askForCustomTVX.value) { - if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) return false; + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) + return false; } else { - if (!collision.sel8()) return false; + if (!collision.sel8()) + return false; } - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 4.f); + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 4.f); - if (removeITSROFrameBorder.value && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) return false; - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 5.f); + if (removeITSROFrameBorder.value && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + return false; + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 5.f); - if (removeNoSameBunchPileup.value && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) return false; - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 6.f); + if (removeNoSameBunchPileup.value && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) + return false; + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 6.f); - if (requireIsGoodZvtxFT0vsPV.value && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) return false; - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 7.f); + if (requireIsGoodZvtxFT0vsPV.value && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) + return false; + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 7.f); - if (requireIsVertexITSTPC.value && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) return false; - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 8.f); + if (requireIsVertexITSTPC.value && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) + return false; + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 8.f); - if (removeNoTimeFrameBorder.value && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) return false; - if constexpr (fillHistograms) ue.fill(HIST("evsel"), 9.f); + if (removeNoTimeFrameBorder.value && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) + return false; + if constexpr (fillHistograms) + ue.fill(HIST("evsel"), 9.f); - if (std::abs(collision.posZ()) > cfgCutVertex.value) return false; + if (std::abs(collision.posZ()) > cfgCutVertex.value) + return false; if constexpr (fillHistograms) { ue.fill(HIST("evsel"), 13.f); - if (collision.isInelGt0()) ue.fill(HIST("evsel"), 14.f); - if (collision.isInelGt1()) ue.fill(HIST("evsel"), 15.f); + if (collision.isInelGt0()) + ue.fill(HIST("evsel"), 14.f); + if (collision.isInelGt1()) + ue.fill(HIST("evsel"), 15.f); } - if (cfgINELCut.value == 1 && !collision.isInelGt0()) return false; - if (cfgINELCut.value == 2 && !collision.isInelGt1()) return false; + if (cfgINELCut.value == 1 && !collision.isInelGt0()) + return false; + if (cfgINELCut.value == 2 && !collision.isInelGt1()) + return false; return true; } @@ -329,19 +350,19 @@ struct multiplicitypt { { auto pdgParticle = pdg->GetParticle(particle.pdgCode()); if (!pdgParticle || pdgParticle->Charge() == 0.) - return false; - + return false; + if (!particle.isPhysicalPrimary()) - return false; - + return false; + if (std::abs(particle.eta()) >= cfgCutEtaMax.value) - return false; + return false; if (particle.pt() < cfgTrkLowPtCut.value) - return false; - + return false; + if (std::abs(particle.y()) > cfgCutY.value) - return false; - + return false; + return true; } @@ -351,13 +372,17 @@ struct multiplicitypt { { int pdgCode = std::abs(particle.pdgCode()); int expectedPDG = 0; - - if constexpr (species == kPion) expectedPDG = PDGPion; - else if constexpr (species == kKaon) expectedPDG = PDGKaon; - else if constexpr (species == kProton) expectedPDG = PDGProton; - - if (pdgCode != expectedPDG) return false; - + + if constexpr (species == kPion) + expectedPDG = PDGPion; + else if constexpr (species == kKaon) + expectedPDG = PDGKaon; + else if constexpr (species == kProton) + expectedPDG = PDGProton; + + if (pdgCode != expectedPDG) + return false; + return isGoodPrimary(particle); } @@ -374,11 +399,11 @@ void multiplicitypt::init(InitContext const&) // ======================================================================== // CUSTOM TRACK CUTS INITIALIZATION - MATCHING spectraTOF // ======================================================================== - + 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.SetMinNClustersITS(min_ITS_nClusters.value); @@ -390,14 +415,14 @@ 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", {0.0, 0.1, 0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, @@ -405,14 +430,13 @@ void multiplicitypt::init(InitContext const&) 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0, 12.0, 14.0, 16.0, 18.0, 20.0, 25.0, 30.0, 40.0, 50.0}, - "pT bin limits" - }; + "pT bin limits"}; AxisSpec ptAxis = {ptBinning, "#it{p}_{T} (GeV/#it{c})"}; // ======================================================================== // HISTOGRAM REGISTRY - INCLUSIVE + PARTICLE-SPECIFIC // ======================================================================== - + // Event counting - EXACT spectraTOF approach ue.add("MC/GenRecoCollisions", "Generated and Reconstructed MC Collisions", HistType::kTH1D, {{10, 0.5, 10.5}}); auto hColl = ue.get(HIST("MC/GenRecoCollisions")); @@ -436,11 +460,11 @@ void multiplicitypt::init(InitContext const&) // ======================================================================== // INCLUSIVE CHARGED PARTICLE HISTOGRAMS // ======================================================================== - + // ALL generated primaries (before any physics selection) ue.add("Inclusive/hPtPrimGenAll", "All generated primaries (no cuts);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); - + // Generated primaries AFTER physics selection ue.add("Inclusive/hPtPrimGen", "Generated primaries (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", HistType::kTH1D, {ptAxis}); @@ -466,50 +490,50 @@ void multiplicitypt::init(InitContext const&) // ======================================================================== // PARTICLE-SPECIFIC HISTOGRAMS (Pions, Kaons, Protons) // ======================================================================== - + 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]; - + // Generated histograms 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/hPtPrimGen", name.c_str()), Form("Generated %s (after physics selection);#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), HistType::kTH1D, {ptAxis}); - + // Tracking efficiency 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}); - + // Reconstructed histograms 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}); - + // Measured spectra ue.add(Form("%s/hPtMeasured", name.c_str()), Form("Measured %s;#it{p}_{T} (GeV/#it{c});Counts", symbol.c_str()), HistType::kTH1D, {ptAxis}); - + // PID quality histograms - TPC ONLY if (enablePIDHistograms) { ue.add(Form("%s/hNsigmaTPC", name.c_str()), @@ -517,7 +541,7 @@ void multiplicitypt::init(InitContext const&) HistType::kTH2D, {ptAxis, {200, -10, 10}}); } } - + // ======================================================================== // MONITORING HISTOGRAMS // ======================================================================== @@ -535,12 +559,12 @@ void multiplicitypt::init(InitContext const&) h->GetXaxis()->SetBinLabel(13, "posZ passed"); h->GetXaxis()->SetBinLabel(14, "INEL>0 (final)"); h->GetXaxis()->SetBinLabel(15, "INEL>1 (final)"); - + 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, 2.0 * o2::constants::math::PI}}); 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 EXCLUSIVE PID for INCLUSIVE + PARTICLE-SPECIFIC (Pi, K, p) analysis"; } @@ -553,12 +577,12 @@ void multiplicitypt::processData(CollisionTableData::iterator const& collision, return; } ue.fill(HIST("hvtxZ"), collision.posZ()); - + for (const auto& track : tracks) { if (!passesTrackSelection(track)) { continue; } - + // Inclusive charged particle (always filled) ue.fill(HIST("Inclusive/hPtMeasured"), track.pt()); ue.fill(HIST("hEta"), track.eta()); @@ -566,20 +590,18 @@ void multiplicitypt::processData(CollisionTableData::iterator const& collision, // Exclusive particle identification int bestSpecies = getBestPIDHypothesis(track); - + if (bestSpecies == kPion) { ue.fill(HIST("Pion/hPtMeasured"), track.pt()); if (enablePIDHistograms) { ue.fill(HIST("Pion/hNsigmaTPC"), track.pt(), track.tpcNSigmaPi()); } - } - else if (bestSpecies == kKaon) { + } else if (bestSpecies == kKaon) { ue.fill(HIST("Kaon/hPtMeasured"), track.pt()); if (enablePIDHistograms) { ue.fill(HIST("Kaon/hNsigmaTPC"), track.pt(), track.tpcNSigmaKa()); } - } - else if (bestSpecies == kProton) { + } else if (bestSpecies == kProton) { ue.fill(HIST("Proton/hPtMeasured"), track.pt()); if (enablePIDHistograms) { ue.fill(HIST("Proton/hNsigmaTPC"), track.pt(), track.tpcNSigmaPr()); @@ -592,14 +614,14 @@ void multiplicitypt::processData(CollisionTableData::iterator const& collision, // MC PROCESSING - WITH FIXED PRIMARY FRACTION CALCULATION // ======================================================================== void multiplicitypt::processMC(TrackTableMC const& tracks, - aod::McParticles const& particles, - CollisionTableMCTrue const& mcCollisions, - CollisionTableMC const& collisions) + aod::McParticles const& particles, + CollisionTableMCTrue const& mcCollisions, + CollisionTableMC const& collisions) { LOG(info) << "=== DEBUG processMC START ==="; LOG(info) << "MC collisions: " << mcCollisions.size(); LOG(info) << "Reconstructed collisions: " << collisions.size(); - + // ======================================================================== // STEP 1: Identify which MC collisions are reconstructable // ======================================================================== @@ -608,7 +630,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, for (const auto& mcCollision : mcCollisions) { auto particlesInCollision = particles.sliceBy(perMCCol, mcCollision.globalIndex()); - + if (std::abs(mcCollision.posZ()) > cfgCutVertex.value) { continue; } @@ -618,7 +640,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, if (cfgINELCut.value == 2 && !o2::pwglf::isINELgt1mc(particlesInCollision, pdg)) { continue; } - + reconstructableMCCollisions.insert(mcCollision.globalIndex()); } @@ -636,16 +658,16 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, if (!collision.has_mcCollision()) { continue; } - + const auto& mcCollision = collision.mcCollision_as(); int64_t mcCollId = mcCollision.globalIndex(); - + if (reconstructableMCCollisions.find(mcCollId) == reconstructableMCCollisions.end()) { continue; } - + reconstructedMCCollisions.insert(mcCollId); - + if (isEventSelected(collision)) { selectedMCCollisions.insert(mcCollId); selectedCollisionIndices.insert(collision.globalIndex()); @@ -687,40 +709,42 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, for (const auto& track : tracks) { totalTracksProcessed++; - - if (!track.has_collision()) continue; - + + if (!track.has_collision()) + continue; + const auto& collision = track.collision_as(); - + if (selectedCollisionIndices.find(collision.globalIndex()) == selectedCollisionIndices.end()) { continue; } tracksFromSelectedEvents++; - - if (!passesTrackSelection(track)) continue; + + if (!passesTrackSelection(track)) + continue; tracksPassingSelection++; - + // ======================================================================== // INCLUSIVE CHARGED PARTICLE ANALYSIS // ======================================================================== - + ue.fill(HIST("Inclusive/hPtMeasured"), track.pt()); ue.fill(HIST("Inclusive/hPtAllReco"), track.pt()); ue.fill(HIST("hEta"), track.eta()); ue.fill(HIST("hPhi"), track.phi()); - + // ======================================================================== // EFFICIENCY NUMERATOR: Fill based on TRUE particle type // ======================================================================== - + if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); int pdgCode = std::abs(particle.pdgCode()); - + if (particle.isPhysicalPrimary()) { ue.fill(HIST("Inclusive/hPtNumEff"), particle.pt()); ue.fill(HIST("Inclusive/hPtPrimReco"), track.pt()); - + // Fill particle-specific efficiency numerator based on TRUE type if (pdgCode == PDGPion) { ue.fill(HIST("Pion/hPtNumEff"), particle.pt()); @@ -735,13 +759,13 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Inclusive/hPtSecReco"), track.pt()); } } - + // ======================================================================== // EXCLUSIVE PID - FIXED PRIMARY FRACTION LOGIC // ======================================================================== - + int bestSpecies = getBestPIDHypothesis(track); - + // ======================================================================== // PION CHANNEL // ======================================================================== @@ -749,14 +773,14 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Pion/hPtMeasured"), track.pt()); ue.fill(HIST("Pion/hPtAllReco"), track.pt()); particleTracksIdentified[kPion]++; - + if (enablePIDHistograms) { ue.fill(HIST("Pion/hNsigmaTPC"), track.pt(), track.tpcNSigmaPi()); } - + if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); - + // KEY FIX: Primary fraction = fraction of identified pions that are primary // This includes correctly identified pions AND misidentified kaons/protons // that happen to be primary particles @@ -769,7 +793,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, } } } - + // ======================================================================== // KAON CHANNEL // ======================================================================== @@ -777,14 +801,14 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Kaon/hPtMeasured"), track.pt()); ue.fill(HIST("Kaon/hPtAllReco"), track.pt()); particleTracksIdentified[kKaon]++; - + if (enablePIDHistograms) { ue.fill(HIST("Kaon/hNsigmaTPC"), track.pt(), track.tpcNSigmaKa()); } - + if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); - + // KEY FIX: Primary fraction of identified kaons // A misidentified pion that is primary still contributes to primary fraction if (particle.isPhysicalPrimary()) { @@ -796,7 +820,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, } } } - + // ======================================================================== // PROTON CHANNEL // ======================================================================== @@ -804,14 +828,14 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, ue.fill(HIST("Proton/hPtMeasured"), track.pt()); ue.fill(HIST("Proton/hPtAllReco"), track.pt()); particleTracksIdentified[kProton]++; - + if (enablePIDHistograms) { ue.fill(HIST("Proton/hNsigmaTPC"), track.pt(), track.tpcNSigmaPr()); } - + if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); - + // KEY FIX: Primary fraction of identified protons if (particle.isPhysicalPrimary()) { ue.fill(HIST("Proton/hPtPrimReco"), track.pt()); @@ -823,7 +847,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, } } } - + LOG(info) << "=== DEBUG TRACK COUNTING ==="; LOG(info) << "Total tracks processed: " << totalTracksProcessed; LOG(info) << "Tracks from selected events: " << tracksFromSelectedEvents; @@ -837,7 +861,7 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, LOG(info) << "Protons identified: " << particleTracksIdentified[kProton] << ", primary: " << particleTracksPrimary[kProton] << ", secondary: " << particleTracksSecondary[kProton]; - + // Calculate and log primary fractions if (particleTracksIdentified[kPion] > 0) { float pionPrimFrac = (float)particleTracksPrimary[kPion] / particleTracksIdentified[kPion]; @@ -851,20 +875,19 @@ void multiplicitypt::processMC(TrackTableMC const& tracks, float protonPrimFrac = (float)particleTracksPrimary[kProton] / particleTracksIdentified[kProton]; LOG(info) << "Proton primary fraction: " << protonPrimFrac * 100.0 << "%"; } - + LOG(info) << "=== DEBUG processMC END ==="; } - // ======================================================================== // TRUE MC PROCESSING - WITH PARTICLE-SPECIFIC SIGNAL LOSS // ======================================================================== void multiplicitypt::processTrue(CollisionTableMCTrue const& mcCollisions, - ParticleTableMC const& particles) + ParticleTableMC const& particles) { LOG(info) << "=== DEBUG processTrue START ==="; LOG(info) << "Number of MC collisions: " << mcCollisions.size(); - + int nPassPhysicsSelection = 0; int nParticlesFilledAll = 0; int nParticlesFilledAfterPS = 0; @@ -875,10 +898,10 @@ void multiplicitypt::processTrue(CollisionTableMCTrue const& mcCollisions, for (const auto& mcCollision : mcCollisions) { // Count EVERY generated event ue.fill(HIST("hEventsAllGen"), 1.0); - + ue.fill(HIST("hvtxZmc"), mcCollision.posZ()); auto particlesInCollision = particles.sliceBy(perMCCol, mcCollision.globalIndex()); - + // ======================================================================== // Fill ALL generated primaries BEFORE physics selection // ======================================================================== @@ -887,35 +910,38 @@ void multiplicitypt::processTrue(CollisionTableMCTrue const& mcCollisions, ue.fill(HIST("Inclusive/hPtPrimGenAll"), particle.pt()); nParticlesFilledAll++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Pion/hPtPrimGenAll"), particle.pt()); particleCountAll[kPion]++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Kaon/hPtPrimGenAll"), particle.pt()); particleCountAll[kKaon]++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Proton/hPtPrimGenAll"), particle.pt()); particleCountAll[kProton]++; } } - + // ======================================================================== // Apply physics selection // ======================================================================== - if (std::abs(mcCollision.posZ()) > cfgCutVertex.value) continue; - - if (cfgINELCut.value == 1 && !o2::pwglf::isINELgt0mc(particlesInCollision, pdg)) continue; - if (cfgINELCut.value == 2 && !o2::pwglf::isINELgt1mc(particlesInCollision, pdg)) continue; - + if (std::abs(mcCollision.posZ()) > cfgCutVertex.value) + continue; + + if (cfgINELCut.value == 1 && !o2::pwglf::isINELgt0mc(particlesInCollision, pdg)) + continue; + if (cfgINELCut.value == 2 && !o2::pwglf::isINELgt1mc(particlesInCollision, pdg)) + continue; + // Count physics-selected events ue.fill(HIST("hEventsPassPhysicsSelection"), 1.0); nPassPhysicsSelection++; - + // Fill primaries AFTER physics selection for (const auto& particle : particlesInCollision) { if (isGoodPrimary(particle)) { @@ -923,19 +949,19 @@ void multiplicitypt::processTrue(CollisionTableMCTrue const& mcCollisions, ue.fill(HIST("Inclusive/hPtPrimGen"), particle.pt()); nParticlesFilledAfterPS++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Pion/hPtDenEff"), particle.pt()); ue.fill(HIST("Pion/hPtPrimGen"), particle.pt()); particleCountAfterPS[kPion]++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Kaon/hPtDenEff"), particle.pt()); ue.fill(HIST("Kaon/hPtPrimGen"), particle.pt()); particleCountAfterPS[kKaon]++; } - + if (isGoodPrimarySpecies(particle)) { ue.fill(HIST("Proton/hPtDenEff"), particle.pt()); ue.fill(HIST("Proton/hPtPrimGen"), particle.pt()); @@ -949,7 +975,7 @@ void multiplicitypt::processTrue(CollisionTableMCTrue const& mcCollisions, LOG(info) << "Passing physics selection: " << nPassPhysicsSelection; LOG(info) << "Total primaries (before PS): " << nParticlesFilledAll; LOG(info) << "Total primaries (after PS): " << nParticlesFilledAfterPS; - + LOG(info) << "=== PARTICLE-SPECIFIC STATISTICS ==="; LOG(info) << "Pions - All: " << particleCountAll[kPion] << ", After PS: " << particleCountAfterPS[kPion];