From 54bce21652c1870deb43f67b5f5860dae767c41b Mon Sep 17 00:00:00 2001 From: arkaprabha-saha Date: Fri, 20 Feb 2026 11:20:05 +0100 Subject: [PATCH 1/2] Added the MC process for yield correction of light nuclei production and changes variable names --- PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx | 231 ++++++++++++++------------ 1 file changed, 124 insertions(+), 107 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx b/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx index a0076bff27a..4f236b7653b 100644 --- a/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx +++ b/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx @@ -14,162 +14,179 @@ /// \author Arkaprabha Saha #include "Common/Core/PID/TPCPIDResponse.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/PIDResponseTOF.h" +#include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" - -#include "MathUtils/BetheBlochAleph.h" +#include "DataFormatsTPC/BetheBlochAleph.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" - +#include #include - #include #include #include using namespace o2; using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; + using CollisionWithEvSel = soa::Join; using TotalTracks = soa::Join; +using CollisionMC = soa::Join; +using TracksMC = soa::Join; namespace { -static const std::vector particleName{"d"}; +static const int kPdgCodeHe3 = 1000020030; +static const float kMaxEta = 0.8f; +static const float kMaxGenRapidity = 0.5f; +static const float kMaxVertexZ = 10.f; +static const int kMinTpcCrossedRows = 70; static const double kBetheBlochDefault[6]{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}; -static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; -static const float maxEtaCut = 0.8f; -static const int minTpcCrossedRowsCut = 70; -static const float maxVertexZCut = 10.f; -} // namespace +static const std::vector kParamNamesBB{"p0", "p1", "p2", "p3", "p4", "resolution"}; +static const std::vector kParticleNamesBB{"He3"}; +} struct AntiNucleiTask { - // Histogram registry: for holding histograms - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - - // Configurable track cuts - Configurable trackNclusTPCcut{"trackNclusTPCcut", 70.0f, "min number of TPC clusters"}; - Configurable trackNclusITScut{"trackNclusITScut", 4.0f, "min number of ITS clusters"}; - Configurable maxChi2TPC{"maxChi2TPC", 4.0f, "max chi2 per cluster TPC"}; - Configurable minChi2TPC{"minChi2TPC", 0.0f, "min chi2 per cluster TPC"}; - Configurable chi2ITS{"chi2ITS", 36.0f, "max chi2 per cluster ITS"}; - Configurable trackDCAz{"trackDCAz", 0.1f, "maxDCAz"}; - Configurable trackDCAxy{"trackDCAxy", 0.1f, "maxDCAxy"}; - Configurable tpcNSigmaCut{"tpcNSigmaCut", 3.0f, "tpcNSigmaCut"}; - Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault, 1, 6, particleName, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for deuteron"}; + HistogramRegistry histo{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + Configurable mMinClsTPC{"minClsTPC", 70.0f, ""}; + Configurable mMinClsITS{"minClsITS", 4.0f, ""}; + Configurable mMaxChi2TPC{"maxChi2TPC", 4.0f, ""}; + Configurable mMinChi2TPC{"minChi2TPC", 0.0f, ""}; + Configurable mMaxChi2ITS{"maxChi2ITS", 36.0f, ""}; + Configurable mMaxDCAZ{"maxDCAZ", 0.02f, ""}; + Configurable mMaxDCAXY{"maxDCAXY", 0.02f, ""}; + Configurable mMaxNsigmaTPC{"maxNsigmaTPC", 3.0f, ""}; + Configurable> mBetheBlochParams{"betheBlochParams", {kBetheBlochDefault, 1, 6, kParticleNamesBB, kParamNamesBB}, ""}; void init(InitContext const&) - { // Defining the Histogram Axes - ConfigurableAxis etaAxis{"etaAxis", {16, -0.8, +0.8}, "#eta"}; - ConfigurableAxis phiAxis{"phiAxis", {70, 0.f, 7.f}, "#phi"}; - ConfigurableAxis zVtxAxis{"zVtxAxis", {100, -20.f, 20.f}, "Primary Vertex z (cm)"}; - ConfigurableAxis nSigmaAxis{"nSigmaAxis", {50, -5.f, 5.f}, "N_{#sigma}"}; - ConfigurableAxis ptAxis{"ptAxis", {200, -10.0f, 10.0f}, "p_{T} (GeV/c)"}; - ConfigurableAxis centAxis{"centAxis", {100, 0, 100.0f}, "Centrality"}; - ConfigurableAxis momAxis{"momAxis", {5.e2, 0.f, 5.f}, "momentum axis binning"}; - ConfigurableAxis tpcAxis{"tpcAxis", {4.e2, 0.f, 4.e3f}, "tpc signal axis binning"}; - - // Creating histograms - histos.add("RawzVtx", "RawzVtx", kTH1F, {{zVtxAxis, "Primary Vertex z (cm)"}}); - histos.add("zVtx", "zVtx", kTH1F, {{zVtxAxis, "Primary Vertex z (cm)"}}); - histos.add("RawEta", "RawEta", kTH1F, {{etaAxis, "#eta"}}); - histos.add("Eta", "Eta", kTH1F, {{etaAxis, "#eta"}}); - histos.add("RawPhi", "RawPhi", kTH1F, {{phiAxis, "#phi (rad)"}}); - histos.add("Phi", "Phi", kTH1F, {{phiAxis, "#phi (rad)"}}); - histos.add("RawPt", "RawPt", kTH1F, {{ptAxis, "#it{p}_{T} (GeV/#it{c})"}}); - histos.add("Pt", "Pt", kTH1F, {{ptAxis, "#it{p}_{T} (GeV/#it{c})"}}); - histos.add("TpcSignal", "TpcSignal", kTH2F, {{momAxis, "#it{p}_{TPC} (GeV/#it{c})"}, {tpcAxis, "d#it{E}/d#it{x}_{TPC}"}}); - histos.add("RawtpcNSigma", "RawtpcNSigma", kTH3F, {{centAxis, "Centrality"}, {ptAxis, "#it{p}_{T} (GeV/#it{c})"}, {nSigmaAxis, "N_{#sigma}"}}); - histos.add("tpcNSigma", "tpcNSigma", kTH3F, {{centAxis, "Centrality"}, {ptAxis, "#it{p}_{T} (GeV/#it{c})"}, {nSigmaAxis, "N_{#sigma}"}}); - histos.add("RawtofNSigma", "RawtofNSigma", kTH3F, {{centAxis, "Centrality"}, {ptAxis, "#it{p}_{T} (GeV/#it{c})"}, {nSigmaAxis, "N_{#sigma}"}}); - histos.add("tofNSigma", "tofNSigma", kTH3F, {{centAxis, "Centrality"}, {ptAxis, "#it{p}_{T} (GeV/#it{c})"}, {nSigmaAxis, "N_{#sigma}"}}); + { + ConfigurableAxis axisEta{"eta", {16, -0.8, +0.8}, "#eta"}; + ConfigurableAxis axisPhi{"phi", {70, 0.f, 7.f}, "#phi (rad)"}; + ConfigurableAxis axisZVtx{"zVtx", {100, -20.f, 20.f}, "Vertex z (cm)"}; + ConfigurableAxis axisNSigma{"nSigma", {50, -5.f, 5.f}, "n#sigma"}; + ConfigurableAxis axisPt{"pt", {200, -10.0f, 10.0f}, "#it{p}_{T} (GeV/#it{c})"}; + ConfigurableAxis axisMomTPC{"momTPC", {5.e2, 0.f, 5.f}, "#it{p}_{TPC} (GeV/#it{c})"}; + ConfigurableAxis axisSignalTPC{"signalTPC", {4.e2, 0.f, 4.e3f}, "d#it{E}/d#it{x}_{TPC} (a.u.)"}; + ConfigurableAxis axisCent{"centAxis", {100, 0, 100.0f}, "Centrality"}; + + histo.add("ZVtx_Raw", "Raw Z-vertex", kTH1F, {{axisZVtx}}); + histo.add("ZVtx", "Z-vertex", kTH1F, {{axisZVtx}}); + histo.add("Eta_Raw", "Raw Eta", kTH1F, {{axisEta}}); + histo.add("Eta", "Eta", kTH1F, {{axisEta}}); + histo.add("Phi_Raw", "Raw Phi", kTH1F, {{axisPhi}}); + histo.add("Phi", "Phi", kTH1F, {{axisPhi}}); + histo.add("Pt_Raw", "Raw Pt", kTH1F, {{axisPt}}); + histo.add("Pt", "Pt", kTH1F, {{axisPt}}); + histo.add("TPCSignal", "TPC dEdx vs Momentum", kTH2F, {{axisMomTPC}, {axisSignalTPC}}); + histo.add("NSigmaTPC_Raw", "Raw TPC nSigma", kTH3F, {{axisCent}, {axisPt}, {axisNSigma}}); + histo.add("NSigmaTOF_Raw", "Raw TOF nSigma", kTH3F, {{axisCent}, {axisPt}, {axisNSigma}}); + histo.add("NSigmaTPC", "TPC nSigma", kTH3F, {{axisCent}, {axisPt}, {axisNSigma}}); + histo.add("NSigmaTOF", "TOF nSigma", kTH3F, {{axisCent}, {axisPt}, {axisNSigma}}); + histo.add("GenPt", "Generated He3 (|y|<0.5)", kTH1F, {{axisPt}}); + histo.add("RecPt", "Reconstructed He3 (|#eta|<0.8)", kTH1F, {{axisPt}}); } - // Function to apply track cuts template - bool isGoodTrack(const T& track) + bool passTrackSelection(const T& track) { - if (track.eta() > maxEtaCut) - return false; - if (track.tpcNClsFound() < trackNclusTPCcut) - return false; - if (track.tpcNClsCrossedRows() < minTpcCrossedRowsCut) - return false; - if (track.itsNCls() < trackNclusITScut) - return false; - if (track.tpcChi2NCl() > maxChi2TPC) - return false; - if (track.tpcChi2NCl() < minChi2TPC) - return false; - if (track.itsChi2NCl() > chi2ITS) - return false; - if (std::abs(track.dcaXY()) > trackDCAxy) - return false; - if (std::abs(track.dcaZ()) > trackDCAz) - return false; - + if (std::abs(track.eta()) > kMaxEta) return false; + if (track.tpcNClsFound() < mMinClsTPC) return false; + if (track.tpcNClsCrossedRows() < kMinTpcCrossedRows) return false; + if (track.itsNCls() < mMinClsITS) return false; + if (track.tpcChi2NCl() > mMaxChi2TPC || track.tpcChi2NCl() < mMinChi2TPC) return false; + if (track.itsChi2NCl() > mMaxChi2ITS) return false; + if (std::abs(track.dcaXY()) > mMaxDCAXY || std::abs(track.dcaZ()) > mMaxDCAZ) return false; return true; } - // The process function + static float getSignedPtMC(auto const& particle) + { + return (particle.pdgCode() > 0) ? particle.pt() : -particle.pt(); + } + + static bool isPrimaryHe3(auto const& particle) + { + return std::abs(particle.pdgCode()) == kPdgCodeHe3 && particle.isPhysicalPrimary(); + } + void process(CollisionWithEvSel::iterator const& collision, TotalTracks const& tracks) { - // Event Selection - if (std::abs(collision.posZ()) > maxVertexZCut) { - return; - } + if (std::abs(collision.posZ()) > kMaxVertexZ) return; + histo.fill(HIST("ZVtx_Raw"), collision.posZ()); - // Filling the z-vertex histogram before the event selection cuts. - histos.fill(HIST("RawzVtx"), collision.posZ()); + if (!collision.sel8()) return; + histo.fill(HIST("ZVtx"), collision.posZ()); - // Applying the built-in O2 event selection (sel8). - if (!collision.sel8()) { - return; + for (const auto& track : tracks) { + double expectedSignal = tpc::BetheBlochAleph( + static_cast(track.tpcInnerParam()), + mBetheBlochParams->get("p0"), mBetheBlochParams->get("p1"), + mBetheBlochParams->get("p2"), mBetheBlochParams->get("p3"), + mBetheBlochParams->get("p4")); + + float nSigmaTPC = static_cast((track.tpcSignal() - expectedSignal) / (expectedSignal * mBetheBlochParams->get("resolution"))); + float signedPt = (track.sign() > 0) ? 2 * track.pt() : -2 * track.pt(); + + histo.fill(HIST("Eta_Raw"), track.eta()); + histo.fill(HIST("Phi_Raw"), track.phi()); + histo.fill(HIST("Pt_Raw"), signedPt); + histo.fill(HIST("NSigmaTPC_Raw"), collision.centFT0C(), signedPt, nSigmaTPC); + histo.fill(HIST("NSigmaTOF_Raw"), collision.centFT0C(), signedPt, track.tofNSigmaDe()); + + if (passTrackSelection(track)) { + histo.fill(HIST("Eta"), track.eta()); + histo.fill(HIST("Phi"), track.phi()); + histo.fill(HIST("Pt"), signedPt); + histo.fill(HIST("NSigmaTPC"), collision.centFT0C(), signedPt, nSigmaTPC); + histo.fill(HIST("TPCSignal"), track.tpcInnerParam(), track.tpcSignal()); + + if (std::abs(nSigmaTPC) < mMaxNsigmaTPC) { + histo.fill(HIST("NSigmaTOF"), collision.centFT0C(), signedPt, track.tofNSigmaDe()); + } + } } + } - // Filling the z-vertex histogram after the event selection cuts. - histos.fill(HIST("zVtx"), collision.posZ()); + void processMC(CollisionMC const& collisions, aod::McCollisions const&, TracksMC const& tracks, aod::McParticles const& particlesMC) + { + for (const auto& collision : collisions) { + if (std::abs(collision.posZ()) > kMaxVertexZ || !collision.sel8()) continue; + + for (const auto& particle : particlesMC) { + if (particle.mcCollisionId() == collision.mcCollisionId() && isPrimaryHe3(particle)) { + if (std::abs(particle.y()) < kMaxGenRapidity) { + histo.fill(HIST("GenPt"), getSignedPtMC(particle)); + } + } + } - // Track Selection - for (const auto& track : tracks) { + for (const auto& track : tracks) { + if (track.collisionId() != collision.index() || !passTrackSelection(track) || track.mcParticleId() == -1) continue; - double expBethe{common::BetheBlochAleph(static_cast(track.tpcInnerParam()), cfgBetheBlochParams->get("p0"), cfgBetheBlochParams->get("p1"), cfgBetheBlochParams->get("p2"), cfgBetheBlochParams->get("p3"), cfgBetheBlochParams->get("p4"))}; - double expSigma{expBethe * cfgBetheBlochParams->get("resolution")}; - float tpcNSigma = static_cast((track.tpcSignal() - expBethe) / expSigma); - - float pt = track.sign() > 0 ? 2 * track.pt() : -2 * track.pt(); - // Filling histograms with track data before applying any cuts. - histos.fill(HIST("RawEta"), track.eta()); - histos.fill(HIST("RawPhi"), track.phi()); - histos.fill(HIST("RawPt"), pt); - histos.fill(HIST("RawtpcNSigma"), collision.centFT0C(), pt, tpcNSigma); - histos.fill(HIST("RawtofNSigma"), collision.centFT0C(), pt, track.tofNSigmaDe()); - - // If the track is good, fill the "after cuts" histograms. - if (isGoodTrack(track)) { - histos.fill(HIST("Eta"), track.eta()); - histos.fill(HIST("Phi"), track.phi()); - histos.fill(HIST("Pt"), pt); - histos.fill(HIST("tpcNSigma"), collision.centFT0C(), pt, tpcNSigma); - histos.fill(HIST("TpcSignal"), track.tpcInnerParam(), track.tpcSignal()); - - if (std::abs(tpcNSigma) < tpcNSigmaCut) { - histos.fill(HIST("tofNSigma"), collision.centFT0C(), pt, track.tofNSigmaDe()); + auto particle = particlesMC.iteratorAt(track.mcParticleId()); + if (isPrimaryHe3(particle)) { + if (std::abs(track.rapidity(MassHelium3)) < kMaxGenRapidity) { + histo.fill(HIST("RecPt"), getSignedPtMC(particle)); + } } } } } + PROCESS_SWITCH(AntiNucleiTask, processMC, "processMC", true); PROCESS_SWITCH(AntiNucleiTask, process, "process", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + return WorkflowSpec{adaptAnalysisTask(cfgc)}; } From 3a88785cc0e56ff01cf48169b78e96a176f3e302 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Fri, 20 Feb 2026 10:24:46 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx | 38 ++++++++++++++++++--------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx b/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx index 4f236b7653b..385a2dd4f5d 100644 --- a/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx +++ b/PWGLF/Tasks/Nuspex/AntiNucleiTask.cxx @@ -21,13 +21,16 @@ #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" + #include "DataFormatsTPC/BetheBlochAleph.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" + #include #include + #include #include #include @@ -52,7 +55,7 @@ static const int kMinTpcCrossedRows = 70; static const double kBetheBlochDefault[6]{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}; static const std::vector kParamNamesBB{"p0", "p1", "p2", "p3", "p4", "resolution"}; static const std::vector kParticleNamesBB{"He3"}; -} +} // namespace struct AntiNucleiTask { HistogramRegistry histo{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -98,13 +101,20 @@ struct AntiNucleiTask { template bool passTrackSelection(const T& track) { - if (std::abs(track.eta()) > kMaxEta) return false; - if (track.tpcNClsFound() < mMinClsTPC) return false; - if (track.tpcNClsCrossedRows() < kMinTpcCrossedRows) return false; - if (track.itsNCls() < mMinClsITS) return false; - if (track.tpcChi2NCl() > mMaxChi2TPC || track.tpcChi2NCl() < mMinChi2TPC) return false; - if (track.itsChi2NCl() > mMaxChi2ITS) return false; - if (std::abs(track.dcaXY()) > mMaxDCAXY || std::abs(track.dcaZ()) > mMaxDCAZ) return false; + if (std::abs(track.eta()) > kMaxEta) + return false; + if (track.tpcNClsFound() < mMinClsTPC) + return false; + if (track.tpcNClsCrossedRows() < kMinTpcCrossedRows) + return false; + if (track.itsNCls() < mMinClsITS) + return false; + if (track.tpcChi2NCl() > mMaxChi2TPC || track.tpcChi2NCl() < mMinChi2TPC) + return false; + if (track.itsChi2NCl() > mMaxChi2ITS) + return false; + if (std::abs(track.dcaXY()) > mMaxDCAXY || std::abs(track.dcaZ()) > mMaxDCAZ) + return false; return true; } @@ -120,10 +130,12 @@ struct AntiNucleiTask { void process(CollisionWithEvSel::iterator const& collision, TotalTracks const& tracks) { - if (std::abs(collision.posZ()) > kMaxVertexZ) return; + if (std::abs(collision.posZ()) > kMaxVertexZ) + return; histo.fill(HIST("ZVtx_Raw"), collision.posZ()); - if (!collision.sel8()) return; + if (!collision.sel8()) + return; histo.fill(HIST("ZVtx"), collision.posZ()); for (const auto& track : tracks) { @@ -159,7 +171,8 @@ struct AntiNucleiTask { void processMC(CollisionMC const& collisions, aod::McCollisions const&, TracksMC const& tracks, aod::McParticles const& particlesMC) { for (const auto& collision : collisions) { - if (std::abs(collision.posZ()) > kMaxVertexZ || !collision.sel8()) continue; + if (std::abs(collision.posZ()) > kMaxVertexZ || !collision.sel8()) + continue; for (const auto& particle : particlesMC) { if (particle.mcCollisionId() == collision.mcCollisionId() && isPrimaryHe3(particle)) { @@ -170,7 +183,8 @@ struct AntiNucleiTask { } for (const auto& track : tracks) { - if (track.collisionId() != collision.index() || !passTrackSelection(track) || track.mcParticleId() == -1) continue; + if (track.collisionId() != collision.index() || !passTrackSelection(track) || track.mcParticleId() == -1) + continue; auto particle = particlesMC.iteratorAt(track.mcParticleId()); if (isPrimaryHe3(particle)) {