From ec18b06fd10d6f0559db5a0240695143dca6e144 Mon Sep 17 00:00:00 2001 From: Leonard Date: Sun, 15 Mar 2026 22:15:41 +0100 Subject: [PATCH 01/11] Add jetHadronsPID and update CMakeLists --- PWGJE/Tasks/CMakeLists.txt | 4 + PWGJE/Tasks/jetHadronsPID.cxx | 242 ++++++++++++++++++++++++++++++++++ 2 files changed, 246 insertions(+) create mode 100644 PWGJE/Tasks/jetHadronsPID.cxx diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 55ec23b2c44..5cb4350b38a 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -403,4 +403,8 @@ if(FastJet_FOUND) SOURCES bjetCentMult.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(pid-hadrons-in-jets + SOURCES jetHadronsPID.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) endif() diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx new file mode 100644 index 00000000000..48eec9da1fd --- /dev/null +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -0,0 +1,242 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. + +#include "Framework/runDataProcessing.h" +#include "PWGJE/Core/JetBkgSubUtils.h" +#include "PWGJE/Core/JetUtilities.h" +#include "PWGJE/DataModel/Jet.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/Logger.h" +#include "ReconstructionDataFormats/DCA.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" + +#include +#include +#include +#include + +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod; +using namespace o2::constants::physics; +using namespace o2::constants::math; + +// Define convenient aliases for commonly used table joins +using SelectedCollisions = soa::Join; + +// CRITICAL FIX: Replaced Proton/Deuteron/Helium PID tables with Pion PID tables +using PionTracks = soa::Join; + +struct PIDHadronsInJets { + + // Histogram registry for data + HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; + + // Parameters for ppRef analysis + Configurable isppRefAnalysis{"isppRefAnalysis", false, "Is ppRef analysis"}; + Configurable cfgAreaFrac{"cfgAreaFrac", 0.6, "fraction of jet area"}; + Configurable cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"}; + Configurable cfgMinPtTrack{"cfgMinPtTrack", 0.1, "minimum pt of tracks for jet reconstruction"}; + + // Event selection criteria + Configurable rejectITSROFBorder{"rejectITSROFBorder", true, "Reject events near the ITS ROF border"}; + Configurable rejectTFBorder{"rejectTFBorder", true, "Reject events near the TF border"}; + Configurable requireVtxITSTPC{"requireVtxITSTPC", true, "Require at least one ITS-TPC matched track"}; + Configurable rejectSameBunchPileup{"rejectSameBunchPileup", true, "Reject events with same-bunch pileup collisions"}; + Configurable requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "Require consistent FT0 vs PV z-vertex"}; + Configurable requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "Require at least one vertex track matched to TOF"}; + + // Jet selection parameters + Configurable minJetPt{"minJetPt", 10.0, "Minimum pt of the jet after bkg subtraction"}; + Configurable maxJetPt{"maxJetPt", 1e+06, "Maximum pt of the jet after bkg subtraction"}; + Configurable rJet{"rJet", 0.4, "Jet resolution parameter R"}; + Configurable zVtx{"zVtx", 10.0, "Maximum zVertex"}; + Configurable applyAreaCut{"applyAreaCut", true, "apply area cut"}; + Configurable maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"}; + Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; + + // Track quality parameters + Configurable requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"}; + Configurable minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"}; + Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 100, "minimum number of TPC crossed pad rows"}; + Configurable minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"}; + Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"}; + Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"}; + Configurable minPt{"minPt", 0.3, "minimum pt of the tracks"}; + Configurable minEta{"minEta", -0.8, "minimum eta"}; + Configurable maxEta{"maxEta", +0.8, "maximum eta"}; + Configurable maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"}; + Configurable maxDcaz{"maxDcaz", 0.05, "Maximum DCAz"}; + + Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; + + // Utility object for jet background subtraction methods + JetBkgSubUtils backgroundSub; + + // Initialize ITS PID Response object + o2::aod::ITSResponse itsResponse; + + void init(InitContext const&) + { + if (setMCDefaultItsParams) { + itsResponse.setMCDefaultParameters(); + } + + // pid of pions + registryData.add("pion_jet_tpc", "TPC Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_jet_tof", "TOF Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("pion_jet_its", "ITS Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{ITS}"}}); + } + + // ITS hit helper + template + bool hasITSHit(const TrackIts& track, int layer) + { + int ibit = layer - 1; + return (track.itsClusterMap() & (1 << ibit)); + } + + // Single-track selection for jet reconstruction + template + bool passedTrackSelectionForJetReconstruction(const JetTrack& track) + { + static constexpr int MinTpcCr = 70; + static constexpr double MaxChi2Tpc = 4.0; + static constexpr double MaxChi2Its = 36.0; + static constexpr double DcaxyMaxTrackPar0 = 0.0105; + static constexpr double DcaxyMaxTrackPar1 = 0.035; + static constexpr double DcaxyMaxTrackPar2 = 1.1; + static constexpr double DcazMaxTrack = 2.0; + + if (!track.hasITS() || !track.hasTPC()) return false; + if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false; + if (track.tpcNClsCrossedRows() < MinTpcCr) return false; + if (track.tpcChi2NCl() > MaxChi2Tpc) return false; + if (track.itsChi2NCl() > MaxChi2Its) return false; + if (std::fabs(track.eta()) > maxEta) return false; + if (track.pt() < cfgMinPtTrack) return false; + if (std::fabs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) return false; + if (std::fabs(track.dcaZ()) > DcazMaxTrack) return false; + return true; + } + + // Single-track selection for constituents + template + bool passedTrackSelection(const PionTrack& track) + { + if (requirePvContributor && !(track.isPVContributor())) return false; + if (!track.hasITS() || !track.hasTPC()) return false; + if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false; + if (track.itsNCls() < minItsNclusters) return false; + if (track.tpcNClsCrossedRows() < minTpcNcrossedRows) return false; + if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) return false; + if (track.itsChi2NCl() > maxChiSquareIts) return false; + if (track.eta() < minEta || track.eta() > maxEta) return false; + if (track.pt() < minPt) return false; + return true; + } + + // Process Data + void processData(SelectedCollisions::iterator const& collision, PionTracks const& tracks) + { + // Apply standard event selection + if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; + + // Build FastJet particles + int id(-1); + std::vector fjParticles; + for (auto const& track : tracks) { + id++; + if (!passedTrackSelectionForJetReconstruction(track)) continue; + + fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged)); + fourMomentum.set_user_index(id); + fjParticles.emplace_back(fourMomentum); + } + + if (fjParticles.empty()) return; + + // Cluster particles + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); + fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); + fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); + std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); + auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); + + // Loop over reconstructed jets + for (const auto& jet : jets) { + + if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; + if (isppRefAnalysis && std::fabs(jet.eta()) > cfgEtaJetMax) continue; + + auto jetForSub = jet; + fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); + + if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) continue; + if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) continue; + + double normalizedJetArea = jet.area() / (PI * rJet * rJet); + if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue; + if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue; + + std::vector jetConstituents = jet.constituents(); + + // loop to fill historgrams + for (const auto& particle : jetConstituents) { + + auto const& track = tracks.iteratorAt(particle.user_index()); + + // Constituent Track Selection (includes DCA checks) + if (!passedTrackSelection(track)) continue; + if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + + double pt = track.pt(); + + // Check Pion PID (+/- 3 Sigma) + double nsigmaTPCPi = track.tpcNSigmaPi(); + double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); + + // Fill TPC + if (std::abs(nsigmaTPCPi) <= 3.0) { + registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi); + } + + // Fill TOF + if (track.hasTOF()) { + double nsigmaTOFPi = track.tofNSigmaPi(); + if (std::abs(nsigmaTOFPi) <= 3.0) { + registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi); + } + } + + // Fill ITS + if (std::abs(nSigmaITSPi) <= 3.0) { + registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} \ No newline at end of file From b7d49f8291bb4f26b9bd470e48c67f6f6bced738 Mon Sep 17 00:00:00 2001 From: Leonard Date: Mon, 16 Mar 2026 12:33:09 +0100 Subject: [PATCH 02/11] small changes to jetHadronsPID --- PWGJE/Tasks/jetHadronsPID.cxx | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index 48eec9da1fd..d46c97c4690 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -36,7 +36,7 @@ using namespace o2::constants::math; // Define convenient aliases for commonly used table joins using SelectedCollisions = soa::Join; -// CRITICAL FIX: Replaced Proton/Deuteron/Helium PID tables with Pion PID tables + using PionTracks = soa::Join; struct PIDHadronsInJets { @@ -149,7 +149,7 @@ struct PIDHadronsInJets { } // Process Data - void processData(SelectedCollisions::iterator const& collision, PionTracks const& tracks) + void process(SelectedCollisions::iterator const& collision, PionTracks const& tracks) { // Apply standard event selection if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; @@ -214,12 +214,12 @@ struct PIDHadronsInJets { double nsigmaTPCPi = track.tpcNSigmaPi(); double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); - // Fill TPC + // Fill TPC histograms if (std::abs(nsigmaTPCPi) <= 3.0) { registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi); } - // Fill TOF + // Fill TOF histograms if (track.hasTOF()) { double nsigmaTOFPi = track.tofNSigmaPi(); if (std::abs(nsigmaTOFPi) <= 3.0) { @@ -227,7 +227,7 @@ struct PIDHadronsInJets { } } - // Fill ITS + // Fill ITS histograms if (std::abs(nSigmaITSPi) <= 3.0) { registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); } From db4168ca8f2dcd6930c34540c57599a079fca3e7 Mon Sep 17 00:00:00 2001 From: Leonard Date: Fri, 20 Mar 2026 21:29:40 +0100 Subject: [PATCH 03/11] fixing pid in jetHadronsPID --- PWGJE/Tasks/jetHadronsPID.cxx | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index d46c97c4690..ca5058fafa6 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -37,7 +37,7 @@ using namespace o2::constants::math; using SelectedCollisions = soa::Join; -using PionTracks = soa::Join; +using PionTracks = soa::Join; struct PIDHadronsInJets { @@ -210,24 +210,24 @@ struct PIDHadronsInJets { double pt = track.pt(); - // Check Pion PID (+/- 3 Sigma) - double nsigmaTPCPi = track.tpcNSigmaPi(); + // Check Pion PID (+/- 3 Sigma) + double nsigmaTPCPi = track.pidTPC().nSigmaPi(); double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); - // Fill TPC histograms + // Fill TPC if (std::abs(nsigmaTPCPi) <= 3.0) { registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi); } - // Fill TOF histograms + // Fill TOF if (track.hasTOF()) { - double nsigmaTOFPi = track.tofNSigmaPi(); + double nsigmaTOFPi = track.pidTOF().nSigmaPi(); if (std::abs(nsigmaTOFPi) <= 3.0) { registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi); } } - // Fill ITS histograms + // Fill ITS if (std::abs(nSigmaITSPi) <= 3.0) { registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); } From 44648230740bae20617548090675f0d98ea9e2be Mon Sep 17 00:00:00 2001 From: Leonard Date: Fri, 20 Mar 2026 21:48:31 +0100 Subject: [PATCH 04/11] more changes to pid --- PWGJE/Tasks/jetHadronsPID.cxx | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index ca5058fafa6..42106aa4f29 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -37,7 +37,7 @@ using namespace o2::constants::math; using SelectedCollisions = soa::Join; -using PionTracks = soa::Join; +using PionTracks = soa::Join; struct PIDHadronsInJets { @@ -209,11 +209,12 @@ struct PIDHadronsInJets { if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; double pt = track.pt(); + // double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); - // Check Pion PID (+/- 3 Sigma) - double nsigmaTPCPi = track.pidTPC().nSigmaPi(); - double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); + // Check Pion PID (+/- 3 Sigma) + double nsigmaTPCPi = track.tpcNSigmaPi(); + // Fill TPC if (std::abs(nsigmaTPCPi) <= 3.0) { registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi); @@ -221,16 +222,16 @@ struct PIDHadronsInJets { // Fill TOF if (track.hasTOF()) { - double nsigmaTOFPi = track.pidTOF().nSigmaPi(); + double nsigmaTOFPi = track.tofNSigmaPi(); if (std::abs(nsigmaTOFPi) <= 3.0) { registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi); } } - // Fill ITS - if (std::abs(nSigmaITSPi) <= 3.0) { - registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); - } + // // Fill ITS + // if (std::abs(nSigmaITSPi) <= 3.0) { + // registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); + // } } } } From 267a9fa2adb8e1f8641bf90cf3a6c022173a4ed9 Mon Sep 17 00:00:00 2001 From: hzalewski Date: Mon, 23 Mar 2026 10:52:47 +0100 Subject: [PATCH 05/11] add historams for hadrons in jet and ue --- PWGJE/Tasks/jetHadronsPID.cxx | 168 ++++++++++++++++++++++++++++------ 1 file changed, 139 insertions(+), 29 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index 42106aa4f29..e4aea32affa 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -24,6 +24,9 @@ #include #include +#include +#include + #include #include @@ -37,7 +40,10 @@ using namespace o2::constants::math; using SelectedCollisions = soa::Join; -using PionTracks = soa::Join; +using HadronTracks = soa::Join; struct PIDHadronsInJets { @@ -75,6 +81,7 @@ struct PIDHadronsInJets { Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"}; Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"}; Configurable minPt{"minPt", 0.3, "minimum pt of the tracks"}; + Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; // pt cut at 4.0 Configurable minEta{"minEta", -0.8, "minimum eta"}; Configurable maxEta{"maxEta", +0.8, "maximum eta"}; Configurable maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"}; @@ -94,10 +101,64 @@ struct PIDHadronsInJets { itsResponse.setMCDefaultParameters(); } - // pid of pions - registryData.add("pion_jet_tpc", "TPC Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{TPC}"}}); - registryData.add("pion_jet_tof", "TOF Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("pion_jet_its", "ITS Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 6.0, "#it{p}_{T} (GeV/#it{c})"}, {100, -3.0, 3.0, "n#sigma_{ITS}"}}); + // Pions, Kaons, Protons in Jet + registryData.add("pion_jet_tpc", "TPC Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_jet_tof", "TOF Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_jet_tpc", "TPC Kaon PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_jet_tof", "TOF Kaon PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_jet_tpc", "TPC Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_jet_tof", "TOF Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + + // Pions, Kaons, Protons Underlying Event + registryData.add("pion_ue_tpc", "TPC Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_ue_tof", "TOF Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_ue_tpc", "TPC Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_ue_tof", "TOF Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_ue_tpc", "TPC Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_ue_tof", "TOF Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + } + + // Compute two transverse directions orthogonal to vector p + void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) + { + double px = p.X(), py = p.Y(), pz = p.Z(); + double px2 = px * px, py2 = py * py, pz2 = pz * pz; + double pz4 = pz2 * pz2; + + if (px == 0 && py == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; } + if (px == 0 && py != 0) { + double ux = std::sqrt(py2 - pz4 / py2); + double uy = -pz2 / py; + u1.SetXYZ(ux, uy, pz); u2.SetXYZ(-ux, uy, pz); return; + } + if (py == 0 && px != 0) { + double ux = -pz2 / px; + double uy = std::sqrt(px2 - pz4 / px2); + u1.SetXYZ(ux, uy, pz); u2.SetXYZ(ux, -uy, pz); return; + } + + double a = px2 + py2; + double b = 2.0 * px * pz2; + double c = pz4 - py2 * py2 - px2 * py2; + double delta = b * b - 4.0 * a * c; + + if (delta < 0 || a == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; } + double u1x = (-b + std::sqrt(delta)) / (2.0 * a); + u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz); + double u2x = (-b - std::sqrt(delta)) / (2.0 * a); + u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz); + } + + double getDeltaPhi(double a1, double a2) + { + double deltaPhi(0); + double phi1 = TVector2::Phi_0_2pi(a1); + double phi2 = TVector2::Phi_0_2pi(a2); + double diff = std::fabs(phi1 - phi2); + + if (diff <= PI) deltaPhi = diff; + if (diff > PI) deltaPhi = TwoPI - diff; + return deltaPhi; } // ITS hit helper @@ -144,12 +205,12 @@ struct PIDHadronsInJets { if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) return false; if (track.itsChi2NCl() > maxChiSquareIts) return false; if (track.eta() < minEta || track.eta() > maxEta) return false; - if (track.pt() < minPt) return false; + if (track.pt() < minPt || track.pt() > maxPt) return false; return true; } // Process Data - void process(SelectedCollisions::iterator const& collision, PionTracks const& tracks) + void process(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) { // Apply standard event selection if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; @@ -174,14 +235,16 @@ struct PIDHadronsInJets { if (fjParticles.empty()) return; - // Cluster particles + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); + + if (jets.empty()) return; + auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); - // Loop over reconstructed jets for (const auto& jet : jets) { if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; @@ -197,41 +260,88 @@ struct PIDHadronsInJets { if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue; if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue; - std::vector jetConstituents = jet.constituents(); - // loop to fill historgrams + double coneRadius = std::sqrt(jet.area() / PI); + TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); + TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); + getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2); + + if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) continue; + + std::vector jetConstituents = jet.constituents(); for (const auto& particle : jetConstituents) { auto const& track = tracks.iteratorAt(particle.user_index()); - // Constituent Track Selection (includes DCA checks) if (!passedTrackSelection(track)) continue; if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; double pt = track.pt(); - // double nSigmaITSPi = static_cast(itsResponse.nSigmaITS(track)); + if (std::abs(track.tpcNSigmaPi()) <= 3.0) { + registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) { + registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); + } - // Check Pion PID (+/- 3 Sigma) - double nsigmaTPCPi = track.tpcNSigmaPi(); - - // Fill TPC - if (std::abs(nsigmaTPCPi) <= 3.0) { - registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi); + if (std::abs(track.tpcNSigmaKa()) <= 3.0) { + registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) { + registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); } + + if (std::abs(track.tpcNSigmaPr()) <= 3.0) { + registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) { + registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); + } + } + + for (auto const& track : tracks) { + + if (!passedTrackSelection(track)) continue; + if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + + double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); + double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi()); + double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1); - // Fill TOF - if (track.hasTOF()) { - double nsigmaTOFPi = track.tofNSigmaPi(); - if (std::abs(nsigmaTOFPi) <= 3.0) { - registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi); - } + double deltaEtaUe2 = track.eta() - ueAxis2.Eta(); + double deltaPhiUe2 = getDeltaPhi(track.phi(), ueAxis2.Phi()); + double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); + + double maxConeRadius = coneRadius; + if (applyAreaCut) { + maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet; } - // // Fill ITS - // if (std::abs(nSigmaITSPi) <= 3.0) { - // registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi); - // } + if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; + + double pt = track.pt(); + + if (std::abs(track.tpcNSigmaPi()) <= 3.0) { + registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) { + registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); + } + + if (std::abs(track.tpcNSigmaKa()) <= 3.0) { + registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) { + registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); + } + + if (std::abs(track.tpcNSigmaPr()) <= 3.0) { + registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr()); + } + if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) { + registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); + } } } } From afe2378ebd2f2cf58ae1b80e32f73bcee13e6124 Mon Sep 17 00:00:00 2001 From: Leonard Date: Sun, 19 Apr 2026 15:18:01 +0200 Subject: [PATCH 06/11] fixing code and adding mc process --- PWGJE/Tasks/jetHadronsPID.cxx | 388 +++++++++++++++++++++++++++++----- 1 file changed, 338 insertions(+), 50 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index e4aea32affa..fd94ed6bef3 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -29,6 +29,7 @@ #include #include +#include using namespace o2; using namespace o2::framework; @@ -36,27 +37,29 @@ using namespace o2::aod; using namespace o2::constants::physics; using namespace o2::constants::math; -// Define convenient aliases for commonly used table joins using SelectedCollisions = soa::Join; - using HadronTracks = soa::Join; + +using HadronTracksMC = soa::Join; + struct PIDHadronsInJets { - // Histogram registry for data HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - // Parameters for ppRef analysis Configurable isppRefAnalysis{"isppRefAnalysis", false, "Is ppRef analysis"}; Configurable cfgAreaFrac{"cfgAreaFrac", 0.6, "fraction of jet area"}; Configurable cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"}; Configurable cfgMinPtTrack{"cfgMinPtTrack", 0.1, "minimum pt of tracks for jet reconstruction"}; - // Event selection criteria Configurable rejectITSROFBorder{"rejectITSROFBorder", true, "Reject events near the ITS ROF border"}; Configurable rejectTFBorder{"rejectTFBorder", true, "Reject events near the TF border"}; Configurable requireVtxITSTPC{"requireVtxITSTPC", true, "Require at least one ITS-TPC matched track"}; @@ -64,7 +67,6 @@ struct PIDHadronsInJets { Configurable requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "Require consistent FT0 vs PV z-vertex"}; Configurable requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "Require at least one vertex track matched to TOF"}; - // Jet selection parameters Configurable minJetPt{"minJetPt", 10.0, "Minimum pt of the jet after bkg subtraction"}; Configurable maxJetPt{"maxJetPt", 1e+06, "Maximum pt of the jet after bkg subtraction"}; Configurable rJet{"rJet", 0.4, "Jet resolution parameter R"}; @@ -73,26 +75,22 @@ struct PIDHadronsInJets { Configurable maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"}; Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; - // Track quality parameters Configurable requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"}; Configurable minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"}; - Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 100, "minimum number of TPC crossed pad rows"}; + Configurable minTpcNcrossedRows{"minTpcNcrossedRows", 70, "minimum number of TPC crossed pad rows"}; Configurable minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"}; Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"}; Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"}; Configurable minPt{"minPt", 0.3, "minimum pt of the tracks"}; - Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; // pt cut at 4.0 + Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; // more ? Configurable minEta{"minEta", -0.8, "minimum eta"}; Configurable maxEta{"maxEta", +0.8, "maximum eta"}; - Configurable maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"}; - Configurable maxDcaz{"maxDcaz", 0.05, "Maximum DCAz"}; + Configurable maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"}; // more ? + Configurable maxDcaz{"maxDcaz", 0.05, "Maximum DCAz"};// Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; - // Utility object for jet background subtraction methods JetBkgSubUtils backgroundSub; - - // Initialize ITS PID Response object o2::aod::ITSResponse itsResponse; void init(InitContext const&) @@ -101,7 +99,32 @@ struct PIDHadronsInJets { itsResponse.setMCDefaultParameters(); } - // Pions, Kaons, Protons in Jet + registryData.add("pion_pure_tpc", "TPC Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("pion_pure_tof", "TOF Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("pion_pure_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("pion_pure_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("pion_pure_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("pion_pure_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + + registryData.add("kaon_pure_tpc", "TPC Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("kaon_pure_tof", "TOF Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_pure_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("kaon_pure_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("kaon_pure_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("kaon_pure_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + + registryData.add("proton_pure_tpc", "TPC Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); + registryData.add("proton_pure_tof", "TOF Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("proton_pure_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("proton_pure_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("proton_pure_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("proton_pure_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + + + registryData.add("z_vtx", "Z-Vertex Distribution", HistType::kTH1F, {{200, -20.0, 20.0, "Z-Vertex (cm)"}}); + registryData.add("tracks_in_jets", "Number of Tracks Inside Jets", HistType::kTH1F, {{100, 0, 100, "N_{tracks}"}}); + registryData.add("tracks_outside_jets", "Number of Tracks Outside Jets", HistType::kTH1F, {{500, 0, 500, "N_{tracks}"}}); + registryData.add("pion_jet_tpc", "TPC Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_jet_tof", "TOF Pion PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); registryData.add("kaon_jet_tpc", "TPC Kaon PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); @@ -109,16 +132,40 @@ struct PIDHadronsInJets { registryData.add("proton_jet_tpc", "TPC Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_jet_tof", "TOF Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - // Pions, Kaons, Protons Underlying Event registryData.add("pion_ue_tpc", "TPC Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_ue_tof", "TOF Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); registryData.add("kaon_ue_tpc", "TPC Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("kaon_ue_tof", "TOF Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); registryData.add("proton_ue_tpc", "TPC Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_ue_tof", "TOF Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + + registryData.add("pion_jet_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("pion_jet_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("pion_jet_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("pion_jet_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + + registryData.add("kaon_jet_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("kaon_jet_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("kaon_jet_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("kaon_jet_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + + registryData.add("proton_jet_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("proton_jet_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("proton_jet_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("proton_jet_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + + + registryData.add("mc_gen_pion_pt", "Generated Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_rec_pion_pt", "Reconstructed Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + + registryData.add("mc_gen_kaon_pt", "Generated Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_rec_kaon_pt", "Reconstructed Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + + registryData.add("mc_gen_proton_pt", "Generated Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_rec_proton_pt", "Reconstructed Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + } - // Compute two transverse directions orthogonal to vector p void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) { double px = p.X(), py = p.Y(), pz = p.Z(); @@ -161,7 +208,6 @@ struct PIDHadronsInJets { return deltaPhi; } - // ITS hit helper template bool hasITSHit(const TrackIts& track, int layer) { @@ -169,7 +215,6 @@ struct PIDHadronsInJets { return (track.itsClusterMap() & (1 << ibit)); } - // Single-track selection for jet reconstruction template bool passedTrackSelectionForJetReconstruction(const JetTrack& track) { @@ -193,7 +238,6 @@ struct PIDHadronsInJets { return true; } - // Single-track selection for constituents template bool passedTrackSelection(const PionTrack& track) { @@ -209,10 +253,9 @@ struct PIDHadronsInJets { return true; } - // Process Data - void process(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) + + void processForJets(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) { - // Apply standard event selection if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; @@ -221,7 +264,79 @@ struct PIDHadronsInJets { if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; - // Build FastJet particles + + + //pure histograms + registryData.fill(HIST("z_vtx"), collision.posZ()); + + for (auto const& track : tracks) { + if (!passedTrackSelection(track)) continue; + if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + + double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); + + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + + const double ptThreshold = 0.8; + + if (passTpcPi) { + registryData.fill(HIST("pion_pure_tpc"), pt, track.tpcNSigmaPi()); + + if (passTofPi) { + registryData.fill(HIST("pion_pure_tof"), pt, track.tofNSigmaPi()); + } + + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_pure_pt"), pt); + registryData.fill(HIST("pion_pure_eta"), eta); + registryData.fill(HIST("pion_pure_dcaxy"), dcaxy); + registryData.fill(HIST("pion_pure_dcaz"), dcaz); + } + } + + if (passTpcKa) { + registryData.fill(HIST("kaon_pure_tpc"), pt, track.tpcNSigmaKa()); + + if (passTofKa) { + registryData.fill(HIST("kaon_pure_tof"), pt, track.tofNSigmaKa()); + } + + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_pure_pt"), pt); + registryData.fill(HIST("kaon_pure_eta"), eta); + registryData.fill(HIST("kaon_pure_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_pure_dcaz"), dcaz); + } + } + + if (passTpcPr) { + registryData.fill(HIST("proton_pure_tpc"), pt, track.tpcNSigmaPr()); + + if (passTofPr) { + registryData.fill(HIST("proton_pure_tof"), pt, track.tofNSigmaPr()); + } + + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_pure_pt"), pt); + registryData.fill(HIST("proton_pure_eta"), eta); + registryData.fill(HIST("proton_pure_dcaxy"), dcaxy); + registryData.fill(HIST("proton_pure_dcaz"), dcaz); + } + } + + } + + // creating jets int id(-1); std::vector fjParticles; for (auto const& track : tracks) { @@ -235,7 +350,6 @@ struct PIDHadronsInJets { if (fjParticles.empty()) return; - fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); @@ -245,6 +359,11 @@ struct PIDHadronsInJets { auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); + + std::set tracksInJetsSet; + + // loop over jets and fill histograms + for (const auto& jet : jets) { if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; @@ -260,7 +379,6 @@ struct PIDHadronsInJets { if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue; if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue; - double coneRadius = std::sqrt(jet.area() / PI); TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); @@ -269,37 +387,83 @@ struct PIDHadronsInJets { if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) continue; std::vector jetConstituents = jet.constituents(); + + + //loop over particles in jets for (const auto& particle : jetConstituents) { - auto const& track = tracks.iteratorAt(particle.user_index()); + int trackIdx = particle.user_index(); + auto const& track = tracks.iteratorAt(trackIdx); + tracksInJetsSet.insert(trackIdx); + if (!passedTrackSelection(track)) continue; if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - if (std::abs(track.tpcNSigmaPi()) <= 3.0) { + const double ptThreshold = 0.8; + + if (passTpcPi) { registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi()); - } - if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) { - registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); + + if (passTofPi) { + registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); + } + + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_jet_pt"), pt); + registryData.fill(HIST("pion_jet_eta"), eta); + registryData.fill(HIST("pion_jet_dcaxy"), dcaxy); + registryData.fill(HIST("pion_jet_dcaz"), dcaz); + } } - if (std::abs(track.tpcNSigmaKa()) <= 3.0) { + if (passTpcKa) { registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa()); - } - if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) { - registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); + + if (passTofKa) { + registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); + } + + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_jet_pt"), pt); + registryData.fill(HIST("kaon_jet_eta"), eta); + registryData.fill(HIST("kaon_jet_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_jet_dcaz"), dcaz); + } } - if (std::abs(track.tpcNSigmaPr()) <= 3.0) { + if (passTpcPr) { registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr()); + + if (passTofPr) { + registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); + } + + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_jet_pt"), pt); + registryData.fill(HIST("proton_jet_eta"), eta); + registryData.fill(HIST("proton_jet_dcaxy"), dcaxy); + registryData.fill(HIST("proton_jet_dcaz"), dcaz); + } } - if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) { - registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); - } + } + + //loop for UE for (auto const& track : tracks) { if (!passedTrackSelection(track)) continue; @@ -318,33 +482,157 @@ struct PIDHadronsInJets { maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet; } + // Only process tracks inside the UE cone if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; double pt = track.pt(); - if (std::abs(track.tpcNSigmaPi()) <= 3.0) { + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + + const double ptThreshold = 0.8; + + if (passTpcPi) { registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi()); - } - if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) { - registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); + if (passTofPi) { + registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); + } + /* + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_ue_pt"), pt); + registryData.fill(HIST("pion_ue_eta"), eta); + registryData.fill(HIST("pion_ue_dcaxy"), dcaxy); + registryData.fill(HIST("pion_ue_dcaz"), dcaz); + } + */ } - if (std::abs(track.tpcNSigmaKa()) <= 3.0) { + if (passTpcKa) { registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa()); + + if (passTofKa) { + registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); + } + + /* + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_ue_pt"), pt); + registryData.fill(HIST("kaon_ue_eta"), eta); + registryData.fill(HIST("kaon_ue_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_ue_dcaz"), dcaz); + } + */ } - if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) { - registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); - } - if (std::abs(track.tpcNSigmaPr()) <= 3.0) { + + if (passTpcPr) { registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr()); + + if (passTofPr) { + registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); + } + + /* + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_ue_pt"), pt); + registryData.fill(HIST("proton_ue_eta"), eta); + registryData.fill(HIST("proton_ue_dcaxy"), dcaxy); + registryData.fill(HIST("proton_ue_dcaz"), dcaz); + } + */ } - if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) { - registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); + + } + } // end jet loop + + registryData.fill(HIST("tracks_in_jets"), tracksInJetsSet.size()); + + int nTracksOut = 0; + int trackIteratorIdx = 0; + for (auto const& track : tracks) { + if (passedTrackSelection(track) && tracksInJetsSet.find(trackIteratorIdx) == tracksInJetsSet.end()) { + nTracksOut++; } + trackIteratorIdx++; + } + registryData.fill(HIST("tracks_outside_jets"), nTracksOut); + } + + PROCESS_SWITCH(PIDHadronsInJets, processForJets, "Process jets", true); + + + + void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks) + { + if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; + + for (auto const& mcpart : mcParticles) { + if (!mcpart.isPhysicalPrimary()) continue; + if (std::abs(mcpart.eta()) > 0.8) continue; + + int pdg = std::abs(mcpart.pdgCode()); + double pt = mcpart.pt(); + + if (pdg == 211) { + registryData.fill(HIST("mc_gen_pion_pt"), pt); + } else if (pdg == 321) { + registryData.fill(HIST("mc_gen_kaon_pt"), pt); + } else if (pdg == 2212) { + registryData.fill(HIST("mc_gen_proton_pt"), pt); + } + } + + const double ptThreshold = 0.8; + + for (auto const& track : tracks) { + if (!passedTrackSelection(track)) continue; + if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + + double pt = track.pt(); + + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + + if (!track.has_mcParticle()) continue; + auto const& trueParticle = track.mcParticle(); + if (!trueParticle.isPhysicalPrimary()) continue; + + int pdg = std::abs(trueParticle.pdgCode()); + + if (passTpcPi && (pt < ptThreshold || passTofPi)) { + if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); + } + + if (passTpcKa && (pt < ptThreshold || passTofKa)) { + if (pdg == 321) registryData.fill(HIST("mc_rec_kaon_pt"), pt); + } + + if (passTpcPr && (pt < ptThreshold || passTofPr)) { + if (pdg == 2212) registryData.fill(HIST("mc_rec_proton_pt"), pt); } } } + PROCESS_SWITCH(PIDHadronsInJets, processMC, "Run on Monte Carlo", false); + + }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From e68116277c46c2865361949a083bde6e6bcbe973 Mon Sep 17 00:00:00 2001 From: Aleksandra Mulewicz Date: Sat, 25 Apr 2026 18:36:57 +0200 Subject: [PATCH 07/11] Update MC and jet histograms --- PWGJE/Tasks/jetHadronsPID.cxx | 78 +++++++++++++++++++++++++---------- 1 file changed, 56 insertions(+), 22 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index fd94ed6bef3..0d788485d7e 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -98,6 +98,16 @@ struct PIDHadronsInJets { if (setMCDefaultItsParams) { itsResponse.setMCDefaultParameters(); } + registryData.add("n_events", "Event counter", HistType::kTH1F, {{1, 0.5, 1.5, "N_{events}"}}); + registryData.add("n_events_raw", "All events", HistType::kTH1F, {{1, 0.5, 1.5, ""}}); + + registryData.add("jet_pt_raw", "Jet pT raw", HistType::kTH1F, {{200, 0.0, 200.0, "#it{p}_{T}^{raw} (GeV/#it{c})"}}); + registryData.add("jet_pt_subtracted", "Jet pT subtracted", HistType::kTH1F, {{200, 0.0, 200.0, "#it{p}_{T}^{sub} (GeV/#it{c})"}}); + registryData.add("jet_pt_raw_vs_sub", "Raw vs sub jet pT", HistType::kTH2F, {{200, 0, 200}, {200, 0, 200}}); + registryData.add("jet_eta", "Jet eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta_{jet}"}}); + registryData.add("jet_phi", "Jet phi", HistType::kTH1F, {{100, 0.0, TwoPI, "#phi_{jet}"}}); + registryData.add("jet_area", "Jet area", HistType::kTH1F, {{100, 0.0, 1.5, "Area"}}); + registryData.add("jet_n_constituents", "Jet multiplicity", HistType::kTH1I, {{100, 0, 100, "N_{constituents}"}}); registryData.add("pion_pure_tpc", "TPC Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_pure_tof", "TOF Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); @@ -120,7 +130,6 @@ struct PIDHadronsInJets { registryData.add("proton_pure_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); registryData.add("proton_pure_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("z_vtx", "Z-Vertex Distribution", HistType::kTH1F, {{200, -20.0, 20.0, "Z-Vertex (cm)"}}); registryData.add("tracks_in_jets", "Number of Tracks Inside Jets", HistType::kTH1F, {{100, 0, 100, "N_{tracks}"}}); registryData.add("tracks_outside_jets", "Number of Tracks Outside Jets", HistType::kTH1F, {{500, 0, 500, "N_{tracks}"}}); @@ -163,7 +172,17 @@ struct PIDHadronsInJets { registryData.add("mc_gen_proton_pt", "Generated Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); registryData.add("mc_rec_proton_pt", "Reconstructed Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - + + registryData.add("rec_pion_all", "All identified pions", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); + registryData.add("rec_kaon_all", "All identified kaons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); + registryData.add("rec_proton_all", "All identified protons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); + registryData.add("contamination_matrix_pion", "Pion Contamination Matrix;True absolute PDG;Rec #it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{3000, -0.5, 2999.5}, {120, 0.0, 4.0}}); + registryData.add("contamination_matrix_kaon", "Kaon Contamination Matrix;True absolute PDG;Rec #it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{3000, -0.5, 2999.5}, {120, 0.0, 4.0}}); + registryData.add("contamination_matrix_proton", "Proton Contamination Matrix;True absolute PDG;Rec #it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{3000, -0.5, 2999.5}, {120, 0.0, 4.0}}); + + registryData.add("mc_sec_pion_pt", "Reconstructed Pion Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_sec_kaon_pt", "Reconstructed Kaon Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_sec_proton_pt", "Reconstructed Proton Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); } void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) @@ -256,6 +275,8 @@ struct PIDHadronsInJets { void processForJets(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) { + registryData.fill(HIST("n_events_raw"), 1); + if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; @@ -264,9 +285,8 @@ struct PIDHadronsInJets { if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; - - //pure histograms + registryData.fill(HIST("n_events"), 1); registryData.fill(HIST("z_vtx"), collision.posZ()); for (auto const& track : tracks) { @@ -278,7 +298,6 @@ struct PIDHadronsInJets { double dcaxy = track.dcaXY(); double dcaz = track.dcaZ(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); @@ -336,7 +355,6 @@ struct PIDHadronsInJets { } - // creating jets int id(-1); std::vector fjParticles; for (auto const& track : tracks) { @@ -362,15 +380,20 @@ struct PIDHadronsInJets { std::set tracksInJetsSet; - // loop over jets and fill histograms - for (const auto& jet : jets) { + registryData.fill(HIST("jet_pt_raw"), jet.pt()); + registryData.fill(HIST("jet_eta"), jet.eta()); + registryData.fill(HIST("jet_phi"), jet.phi()); + registryData.fill(HIST("jet_area"), jet.area()); + registryData.fill(HIST("jet_n_constituents"), jet.constituents().size()); if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; if (isppRefAnalysis && std::fabs(jet.eta()) > cfgEtaJetMax) continue; auto jetForSub = jet; fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); + registryData.fill(HIST("jet_pt_subtracted"), jetMinusBkg.pt()); + registryData.fill(HIST("jet_pt_raw_vs_sub"), jet.pt(), jetMinusBkg.pt()); if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) continue; if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) continue; @@ -388,8 +411,6 @@ struct PIDHadronsInJets { std::vector jetConstituents = jet.constituents(); - - //loop over particles in jets for (const auto& particle : jetConstituents) { int trackIdx = particle.user_index(); @@ -462,8 +483,6 @@ struct PIDHadronsInJets { } - - //loop for UE for (auto const& track : tracks) { if (!passedTrackSelection(track)) continue; @@ -482,12 +501,10 @@ struct PIDHadronsInJets { maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet; } - // Only process tracks inside the UE cone if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; double pt = track.pt(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); @@ -549,7 +566,7 @@ struct PIDHadronsInJets { } } - } // end jet loop + } registryData.fill(HIST("tracks_in_jets"), tracksInJetsSet.size()); @@ -566,8 +583,6 @@ struct PIDHadronsInJets { PROCESS_SWITCH(PIDHadronsInJets, processForJets, "Process jets", true); - - void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks) { if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; @@ -613,20 +628,39 @@ struct PIDHadronsInJets { if (!track.has_mcParticle()) continue; auto const& trueParticle = track.mcParticle(); - if (!trueParticle.isPhysicalPrimary()) continue; - + int pdg = std::abs(trueParticle.pdgCode()); + bool isPrimary = trueParticle.isPhysicalPrimary(); if (passTpcPi && (pt < ptThreshold || passTofPi)) { - if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); + registryData.fill(HIST("rec_pion_all"), pt); //reconstructed after cuts + registryData.fill(HIST("contamination_matrix_pion"), pdg, pt); + + if (isPrimary) { + if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); // Efficiency + } else { + registryData.fill(HIST("mc_sec_pion_pt"), pt); + } } if (passTpcKa && (pt < ptThreshold || passTofKa)) { - if (pdg == 321) registryData.fill(HIST("mc_rec_kaon_pt"), pt); + registryData.fill(HIST("rec_kaon_all"), pt); + registryData.fill(HIST("contamination_matrix_kaon"), pdg, pt); + if (isPrimary) { + if (pdg == 321) registryData.fill(HIST("mc_rec_kaon_pt"), pt); + } else { + registryData.fill(HIST("mc_sec_kaon_pt"), pt); + } } if (passTpcPr && (pt < ptThreshold || passTofPr)) { - if (pdg == 2212) registryData.fill(HIST("mc_rec_proton_pt"), pt); + registryData.fill(HIST("rec_proton_all"), pt); + registryData.fill(HIST("contamination_matrix_proton"), pdg, pt); + if (isPrimary) { + if (pdg == 2212) registryData.fill(HIST("mc_rec_proton_pt"), pt); + } else { + registryData.fill(HIST("mc_sec_proton_pt"), pt); + } } } } From d989f30c7b2c67f84bcdedbb57c09b100de33e6d Mon Sep 17 00:00:00 2001 From: Leonard <107263905+MyFavoriteGitHub@users.noreply.github.com> Date: Mon, 27 Apr 2026 17:39:26 +0200 Subject: [PATCH 08/11] Rename DPL workflow from 'pid-hadrons-in-jets' to 'jet-hadrons-pid' --- PWGJE/Tasks/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index 5cb4350b38a..c822f837f8a 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -403,7 +403,7 @@ if(FastJet_FOUND) SOURCES bjetCentMult.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) - o2physics_add_dpl_workflow(pid-hadrons-in-jets + o2physics_add_dpl_workflow(jet-hadrons-pid SOURCES jetHadronsPID.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) From 25b7de036a477ee465d48f500391f2043ea3fac0 Mon Sep 17 00:00:00 2001 From: Leonard Date: Tue, 28 Apr 2026 20:20:56 +0200 Subject: [PATCH 09/11] Fixing formating and adding description --- PWGJE/Tasks/jetHadronsPID.cxx | 407 +++++++++++++++++----------------- 1 file changed, 209 insertions(+), 198 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPID.cxx index 0d788485d7e..d0e26ea03ed 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPID.cxx @@ -1,6 +1,21 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// Copyright 2019-2025 CERN and copyright holders of ALICE O2. // See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. // All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file jetHadronsPID.cxx +/// \brief Analysis of hadrons in jets +/// \author Aleksandra Mulewicz, WUT Warsaw, aleksandra.mulewicz@cern.ch +/// \author Leonard Lorenc, WUT Warsaw, leonard.lorenc@cern.ch +/// \author Hubert Zalewski, WUT Warsaw, hubert.zalewski@cern.ch +/// \author Małgorzata Janik malgorzata.janik@cern.ch +/// \author Daniela Ruggiano daniela.ruggiano@cern.ch #include "Framework/runDataProcessing.h" #include "PWGJE/Core/JetBkgSubUtils.h" @@ -29,7 +44,7 @@ #include #include -#include +#include using namespace o2; using namespace o2::framework; @@ -39,16 +54,15 @@ using namespace o2::constants::math; using SelectedCollisions = soa::Join; -using HadronTracks = soa::Join; - -using HadronTracksMC = soa::Join; struct PIDHadronsInJets { @@ -82,12 +96,12 @@ struct PIDHadronsInJets { Configurable maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"}; Configurable maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"}; Configurable minPt{"minPt", 0.3, "minimum pt of the tracks"}; - Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; // more ? + Configurable maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; Configurable minEta{"minEta", -0.8, "minimum eta"}; Configurable maxEta{"maxEta", +0.8, "maximum eta"}; - Configurable maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"}; // more ? - Configurable maxDcaz{"maxDcaz", 0.05, "Maximum DCAz"};// - + Configurable maxDcaxy{"maxDcaxy", 0.2, "Maximum DCAxy"}; + Configurable maxDcaz{"maxDcaz", 0.1, "Maximum DCAz"}; + Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; JetBkgSubUtils backgroundSub; @@ -98,6 +112,7 @@ struct PIDHadronsInJets { if (setMCDefaultItsParams) { itsResponse.setMCDefaultParameters(); } + registryData.add("n_events", "Event counter", HistType::kTH1F, {{1, 0.5, 1.5, "N_{events}"}}); registryData.add("n_events_raw", "All events", HistType::kTH1F, {{1, 0.5, 1.5, ""}}); @@ -141,38 +156,49 @@ struct PIDHadronsInJets { registryData.add("proton_jet_tpc", "TPC Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_jet_tof", "TOF Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("pion_jet_pt", "Pion pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("pion_jet_eta", "Pion Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("pion_jet_dcaxy", "Pion DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("pion_jet_dcaz", "Pion DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + + registryData.add("kaon_jet_pt", "Kaon pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("kaon_jet_eta", "Kaon Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("kaon_jet_dcaxy", "Kaon DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("kaon_jet_dcaz", "Kaon DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + + registryData.add("proton_jet_pt", "Proton pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("proton_jet_eta", "Proton Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("proton_jet_dcaxy", "Proton DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("proton_jet_dcaz", "Proton DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("pion_ue_tpc", "TPC Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_ue_tof", "TOF Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("pion_ue_pt", "Pion pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("pion_ue_eta", "Pion Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("pion_ue_dcaxy", "Pion DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("pion_ue_dcaz", "Pion DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("kaon_ue_tpc", "TPC Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("kaon_ue_tof", "TOF Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); + registryData.add("kaon_ue_pt", "Kaon pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("kaon_ue_eta", "Kaon Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("kaon_ue_dcaxy", "Kaon DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("kaon_ue_dcaz", "Kaon DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("proton_ue_tpc", "TPC Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_ue_tof", "TOF Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - - registryData.add("pion_jet_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("pion_jet_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("pion_jet_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("pion_jet_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - - registryData.add("kaon_jet_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("kaon_jet_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("kaon_jet_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("kaon_jet_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - - registryData.add("proton_jet_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("proton_jet_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("proton_jet_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("proton_jet_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - + registryData.add("proton_ue_pt", "Proton pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("proton_ue_eta", "Proton Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); + registryData.add("proton_ue_dcaxy", "Proton DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("proton_ue_dcaz", "Proton DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); registryData.add("mc_gen_pion_pt", "Generated Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); registryData.add("mc_rec_pion_pt", "Reconstructed Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_gen_kaon_pt", "Generated Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); registryData.add("mc_rec_kaon_pt", "Reconstructed Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_gen_proton_pt", "Generated Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); registryData.add("mc_rec_proton_pt", "Reconstructed Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - + registryData.add("rec_pion_all", "All identified pions", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); registryData.add("rec_kaon_all", "All identified kaons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); registryData.add("rec_proton_all", "All identified protons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); @@ -220,7 +246,7 @@ struct PIDHadronsInJets { double deltaPhi(0); double phi1 = TVector2::Phi_0_2pi(a1); double phi2 = TVector2::Phi_0_2pi(a2); - double diff = std::fabs(phi1 - phi2); + double diff = std::abs(phi1 - phi2); if (diff <= PI) deltaPhi = diff; if (diff > PI) deltaPhi = TwoPI - diff; @@ -250,10 +276,10 @@ struct PIDHadronsInJets { if (track.tpcNClsCrossedRows() < MinTpcCr) return false; if (track.tpcChi2NCl() > MaxChi2Tpc) return false; if (track.itsChi2NCl() > MaxChi2Its) return false; - if (std::fabs(track.eta()) > maxEta) return false; + if (std::abs(track.eta()) > maxEta) return false; if (track.pt() < cfgMinPtTrack) return false; - if (std::fabs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) return false; - if (std::fabs(track.dcaZ()) > DcazMaxTrack) return false; + if (std::abs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) return false; + if (std::abs(track.dcaZ()) > DcazMaxTrack) return false; return true; } @@ -276,8 +302,8 @@ struct PIDHadronsInJets { void processForJets(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) { registryData.fill(HIST("n_events_raw"), 1); - - if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; + + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) return; if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; @@ -285,74 +311,72 @@ struct PIDHadronsInJets { if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; - //pure histograms registryData.fill(HIST("n_events"), 1); registryData.fill(HIST("z_vtx"), collision.posZ()); for (auto const& track : tracks) { - if (!passedTrackSelection(track)) continue; - if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + if (!passedTrackSelection(track)) continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; - double pt = track.pt(); - double eta = track.eta(); - double dcaxy = track.dcaXY(); - double dcaz = track.dcaZ(); + double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); - bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); - bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); - bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); - bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); - bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - const double ptThreshold = 0.8; + const double ptThreshold = 0.8; - if (passTpcPi) { - registryData.fill(HIST("pion_pure_tpc"), pt, track.tpcNSigmaPi()); - - if (passTofPi) { - registryData.fill(HIST("pion_pure_tof"), pt, track.tofNSigmaPi()); - } - - if (pt < ptThreshold || passTofPi) { - registryData.fill(HIST("pion_pure_pt"), pt); - registryData.fill(HIST("pion_pure_eta"), eta); - registryData.fill(HIST("pion_pure_dcaxy"), dcaxy); - registryData.fill(HIST("pion_pure_dcaz"), dcaz); - } + if (passTpcPi) { + registryData.fill(HIST("pion_pure_tpc"), pt, track.tpcNSigmaPi()); + + if (passTofPi) { + registryData.fill(HIST("pion_pure_tof"), pt, track.tofNSigmaPi()); } - if (passTpcKa) { - registryData.fill(HIST("kaon_pure_tpc"), pt, track.tpcNSigmaKa()); - - if (passTofKa) { - registryData.fill(HIST("kaon_pure_tof"), pt, track.tofNSigmaKa()); - } - - if (pt < ptThreshold || passTofKa) { - registryData.fill(HIST("kaon_pure_pt"), pt); - registryData.fill(HIST("kaon_pure_eta"), eta); - registryData.fill(HIST("kaon_pure_dcaxy"), dcaxy); - registryData.fill(HIST("kaon_pure_dcaz"), dcaz); - } + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_pure_pt"), pt); + registryData.fill(HIST("pion_pure_eta"), eta); + registryData.fill(HIST("pion_pure_dcaxy"), dcaxy); + registryData.fill(HIST("pion_pure_dcaz"), dcaz); } + } - if (passTpcPr) { - registryData.fill(HIST("proton_pure_tpc"), pt, track.tpcNSigmaPr()); - - if (passTofPr) { - registryData.fill(HIST("proton_pure_tof"), pt, track.tofNSigmaPr()); - } - - if (pt < ptThreshold || passTofPr) { - registryData.fill(HIST("proton_pure_pt"), pt); - registryData.fill(HIST("proton_pure_eta"), eta); - registryData.fill(HIST("proton_pure_dcaxy"), dcaxy); - registryData.fill(HIST("proton_pure_dcaz"), dcaz); - } + if (passTpcKa) { + registryData.fill(HIST("kaon_pure_tpc"), pt, track.tpcNSigmaKa()); + + if (passTofKa) { + registryData.fill(HIST("kaon_pure_tof"), pt, track.tofNSigmaKa()); } + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_pure_pt"), pt); + registryData.fill(HIST("kaon_pure_eta"), eta); + registryData.fill(HIST("kaon_pure_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_pure_dcaz"), dcaz); + } + } + + if (passTpcPr) { + registryData.fill(HIST("proton_pure_tpc"), pt, track.tpcNSigmaPr()); + + if (passTofPr) { + registryData.fill(HIST("proton_pure_tof"), pt, track.tofNSigmaPr()); + } + + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_pure_pt"), pt); + registryData.fill(HIST("proton_pure_eta"), eta); + registryData.fill(HIST("proton_pure_dcaxy"), dcaxy); + registryData.fill(HIST("proton_pure_dcaz"), dcaz); + } + } } int id(-1); @@ -372,12 +396,11 @@ struct PIDHadronsInJets { fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); - + if (jets.empty()) return; auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); - std::set tracksInJetsSet; for (const auto& jet : jets) { @@ -387,14 +410,14 @@ struct PIDHadronsInJets { registryData.fill(HIST("jet_area"), jet.area()); registryData.fill(HIST("jet_n_constituents"), jet.constituents().size()); - if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; - if (isppRefAnalysis && std::fabs(jet.eta()) > cfgEtaJetMax) continue; + if (!isppRefAnalysis && ((std::abs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; + if (isppRefAnalysis && std::abs(jet.eta()) > cfgEtaJetMax) continue; auto jetForSub = jet; fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); registryData.fill(HIST("jet_pt_subtracted"), jetMinusBkg.pt()); registryData.fill(HIST("jet_pt_raw_vs_sub"), jet.pt(), jetMinusBkg.pt()); - + if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) continue; if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) continue; @@ -406,7 +429,7 @@ struct PIDHadronsInJets { TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2); - + if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) continue; std::vector jetConstituents = jet.constituents(); @@ -415,11 +438,11 @@ struct PIDHadronsInJets { int trackIdx = particle.user_index(); auto const& track = tracks.iteratorAt(trackIdx); - + tracksInJetsSet.insert(trackIdx); if (!passedTrackSelection(track)) continue; - if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; double pt = track.pt(); double eta = track.eta(); @@ -434,64 +457,63 @@ struct PIDHadronsInJets { bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - const double ptThreshold = 0.8; + const double ptThreshold = 0.8; if (passTpcPi) { - registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi()); - - if (passTofPi) { - registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); - } - - if (pt < ptThreshold || passTofPi) { - registryData.fill(HIST("pion_jet_pt"), pt); - registryData.fill(HIST("pion_jet_eta"), eta); - registryData.fill(HIST("pion_jet_dcaxy"), dcaxy); - registryData.fill(HIST("pion_jet_dcaz"), dcaz); - } + registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi()); + + if (passTofPi) { + registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi()); + } + + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_jet_pt"), pt); + registryData.fill(HIST("pion_jet_eta"), eta); + registryData.fill(HIST("pion_jet_dcaxy"), dcaxy); + registryData.fill(HIST("pion_jet_dcaz"), dcaz); + } } if (passTpcKa) { - registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa()); - - if (passTofKa) { - registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); - } - - if (pt < ptThreshold || passTofKa) { - registryData.fill(HIST("kaon_jet_pt"), pt); - registryData.fill(HIST("kaon_jet_eta"), eta); - registryData.fill(HIST("kaon_jet_dcaxy"), dcaxy); - registryData.fill(HIST("kaon_jet_dcaz"), dcaz); - } + registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa()); + + if (passTofKa) { + registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa()); + } + + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_jet_pt"), pt); + registryData.fill(HIST("kaon_jet_eta"), eta); + registryData.fill(HIST("kaon_jet_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_jet_dcaz"), dcaz); + } } if (passTpcPr) { - registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr()); - - if (passTofPr) { - registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); - } - - if (pt < ptThreshold || passTofPr) { - registryData.fill(HIST("proton_jet_pt"), pt); - registryData.fill(HIST("proton_jet_eta"), eta); - registryData.fill(HIST("proton_jet_dcaxy"), dcaxy); - registryData.fill(HIST("proton_jet_dcaz"), dcaz); - } + registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr()); + + if (passTofPr) { + registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr()); + } + + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_jet_pt"), pt); + registryData.fill(HIST("proton_jet_eta"), eta); + registryData.fill(HIST("proton_jet_dcaxy"), dcaxy); + registryData.fill(HIST("proton_jet_dcaz"), dcaz); + } } - } for (auto const& track : tracks) { if (!passedTrackSelection(track)) continue; - if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi()); double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1); - + double deltaEtaUe2 = track.eta() - ueAxis2.Eta(); double deltaPhiUe2 = getDeltaPhi(track.phi(), ueAxis2.Phi()); double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); @@ -504,6 +526,9 @@ struct PIDHadronsInJets { if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; double pt = track.pt(); + double eta = track.eta(); + double dcaxy = track.dcaXY(); + double dcaz = track.dcaZ(); bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); @@ -513,58 +538,46 @@ struct PIDHadronsInJets { bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - const double ptThreshold = 0.8; + const double ptThreshold = 0.8; if (passTpcPi) { - registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi()); - if (passTofPi) { - registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); - } - /* - if (pt < ptThreshold || passTofPi) { - registryData.fill(HIST("pion_ue_pt"), pt); - registryData.fill(HIST("pion_ue_eta"), eta); - registryData.fill(HIST("pion_ue_dcaxy"), dcaxy); - registryData.fill(HIST("pion_ue_dcaz"), dcaz); - } - */ + registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi()); + if (passTofPi) { + registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi()); + } + if (pt < ptThreshold || passTofPi) { + registryData.fill(HIST("pion_ue_pt"), pt); + registryData.fill(HIST("pion_ue_eta"), eta); + registryData.fill(HIST("pion_ue_dcaxy"), dcaxy); + registryData.fill(HIST("pion_ue_dcaz"), dcaz); + } } if (passTpcKa) { - registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa()); - - if (passTofKa) { - registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); - } - - /* - if (pt < ptThreshold || passTofKa) { - registryData.fill(HIST("kaon_ue_pt"), pt); - registryData.fill(HIST("kaon_ue_eta"), eta); - registryData.fill(HIST("kaon_ue_dcaxy"), dcaxy); - registryData.fill(HIST("kaon_ue_dcaz"), dcaz); - } - */ + registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa()); + if (passTofKa) { + registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa()); + } + if (pt < ptThreshold || passTofKa) { + registryData.fill(HIST("kaon_ue_pt"), pt); + registryData.fill(HIST("kaon_ue_eta"), eta); + registryData.fill(HIST("kaon_ue_dcaxy"), dcaxy); + registryData.fill(HIST("kaon_ue_dcaz"), dcaz); + } } - if (passTpcPr) { - registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr()); - - if (passTofPr) { - registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); - } - - /* - if (pt < ptThreshold || passTofPr) { - registryData.fill(HIST("proton_ue_pt"), pt); - registryData.fill(HIST("proton_ue_eta"), eta); - registryData.fill(HIST("proton_ue_dcaxy"), dcaxy); - registryData.fill(HIST("proton_ue_dcaz"), dcaz); - } - */ + registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr()); + if (passTofPr) { + registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr()); + } + if (pt < ptThreshold || passTofPr) { + registryData.fill(HIST("proton_ue_pt"), pt); + registryData.fill(HIST("proton_ue_eta"), eta); + registryData.fill(HIST("proton_ue_dcaxy"), dcaxy); + registryData.fill(HIST("proton_ue_dcaz"), dcaz); + } } - } } @@ -573,19 +586,19 @@ struct PIDHadronsInJets { int nTracksOut = 0; int trackIteratorIdx = 0; for (auto const& track : tracks) { - if (passedTrackSelection(track) && tracksInJetsSet.find(trackIteratorIdx) == tracksInJetsSet.end()) { - nTracksOut++; - } - trackIteratorIdx++; + if (passedTrackSelection(track) && tracksInJetsSet.find(trackIteratorIdx) == tracksInJetsSet.end()) { + nTracksOut++; + } + trackIteratorIdx++; } registryData.fill(HIST("tracks_outside_jets"), nTracksOut); } - PROCESS_SWITCH(PIDHadronsInJets, processForJets, "Process jets", true); + PROCESS_SWITCH(PIDHadronsInJets, processForJets, "Pid Analysis in and outside jets", true); void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks) { - if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return; + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) return; if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; @@ -609,11 +622,11 @@ struct PIDHadronsInJets { } } - const double ptThreshold = 0.8; + const double ptThreshold = 0.8; for (auto const& track : tracks) { if (!passedTrackSelection(track)) continue; - if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; double pt = track.pt(); @@ -626,23 +639,22 @@ struct PIDHadronsInJets { bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - if (!track.has_mcParticle()) continue; + if (!track.has_mcParticle()) continue; auto const& trueParticle = track.mcParticle(); - + int pdg = std::abs(trueParticle.pdgCode()); bool isPrimary = trueParticle.isPhysicalPrimary(); if (passTpcPi && (pt < ptThreshold || passTofPi)) { - registryData.fill(HIST("rec_pion_all"), pt); //reconstructed after cuts + registryData.fill(HIST("rec_pion_all"), pt); registryData.fill(HIST("contamination_matrix_pion"), pdg, pt); - if (isPrimary) { - if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); // Efficiency + if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); } else { registryData.fill(HIST("mc_sec_pion_pt"), pt); } } - + if (passTpcKa && (pt < ptThreshold || passTofKa)) { registryData.fill(HIST("rec_kaon_all"), pt); registryData.fill(HIST("contamination_matrix_kaon"), pdg, pt); @@ -666,7 +678,6 @@ struct PIDHadronsInJets { } PROCESS_SWITCH(PIDHadronsInJets, processMC, "Run on Monte Carlo", false); - }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From 8b4ef7ba220dee317ced639f903c8e1d4f42c028 Mon Sep 17 00:00:00 2001 From: Leonard Date: Tue, 28 Apr 2026 21:09:15 +0200 Subject: [PATCH 10/11] Changing formating in jetHadronsPid.cxx and in CMakeLists.txt --- PWGJE/Tasks/CMakeLists.txt | 2 +- .../{jetHadronsPID.cxx => jetHadronsPid.cxx} | 319 +++++++++++------- 2 files changed, 195 insertions(+), 126 deletions(-) rename PWGJE/Tasks/{jetHadronsPID.cxx => jetHadronsPid.cxx} (83%) diff --git a/PWGJE/Tasks/CMakeLists.txt b/PWGJE/Tasks/CMakeLists.txt index c822f837f8a..f20bd617e33 100644 --- a/PWGJE/Tasks/CMakeLists.txt +++ b/PWGJE/Tasks/CMakeLists.txt @@ -404,7 +404,7 @@ if(FastJet_FOUND) PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(jet-hadrons-pid - SOURCES jetHadronsPID.cxx + SOURCES jetHadronsPid.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) endif() diff --git a/PWGJE/Tasks/jetHadronsPID.cxx b/PWGJE/Tasks/jetHadronsPid.cxx similarity index 83% rename from PWGJE/Tasks/jetHadronsPID.cxx rename to PWGJE/Tasks/jetHadronsPid.cxx index d0e26ea03ed..0d5263e67ce 100644 --- a/PWGJE/Tasks/jetHadronsPID.cxx +++ b/PWGJE/Tasks/jetHadronsPid.cxx @@ -9,15 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file jetHadronsPID.cxx -/// \brief Analysis of hadrons in jets +/// \file jetHadronsPid.cxx +/// \brief Analysis of hadrons in jets /// \author Aleksandra Mulewicz, WUT Warsaw, aleksandra.mulewicz@cern.ch /// \author Leonard Lorenc, WUT Warsaw, leonard.lorenc@cern.ch /// \author Hubert Zalewski, WUT Warsaw, hubert.zalewski@cern.ch /// \author Małgorzata Janik malgorzata.janik@cern.ch /// \author Daniela Ruggiano daniela.ruggiano@cern.ch -#include "Framework/runDataProcessing.h" #include "PWGJE/Core/JetBkgSubUtils.h" #include "PWGJE/Core/JetUtilities.h" #include "PWGJE/DataModel/Jet.h" @@ -30,21 +29,22 @@ #include "Framework/AnalysisTask.h" #include "Framework/HistogramRegistry.h" #include "Framework/Logger.h" +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/DCA.h" #include "ReconstructionDataFormats/PID.h" #include "ReconstructionDataFormats/Track.h" +#include +#include + #include #include #include #include -#include -#include - #include -#include #include +#include using namespace o2; using namespace o2::framework; @@ -65,7 +65,7 @@ using HadronTracksMC = soa::Join; -struct PIDHadronsInJets { +struct jetHadronsPid { HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; @@ -112,7 +112,7 @@ struct PIDHadronsInJets { if (setMCDefaultItsParams) { itsResponse.setMCDefaultParameters(); } - + registryData.add("n_events", "Event counter", HistType::kTH1F, {{1, 0.5, 1.5, "N_{events}"}}); registryData.add("n_events_raw", "All events", HistType::kTH1F, {{1, 0.5, 1.5, ""}}); @@ -126,24 +126,24 @@ struct PIDHadronsInJets { registryData.add("pion_pure_tpc", "TPC Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_pure_tof", "TOF Pion PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("pion_pure_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("pion_pure_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("pion_pure_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - registryData.add("pion_pure_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("pion_pure_pt", "Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_pure_eta", "Pion Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_pure_dcaz", "Pion DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("pion_pure_dcaxy", "Pion DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); registryData.add("kaon_pure_tpc", "TPC Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("kaon_pure_tof", "TOF Kaon PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("kaon_pure_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("kaon_pure_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("kaon_pure_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - registryData.add("kaon_pure_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("kaon_pure_pt", "Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_pure_eta", "Kaon Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_pure_dcaz", "Kaon DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("kaon_pure_dcaxy", "Kaon DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); registryData.add("proton_pure_tpc", "TPC Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_pure_tof", "TOF Proton PID", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("proton_pure_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("proton_pure_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("proton_pure_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - registryData.add("proton_pure_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); + registryData.add("proton_pure_pt", "Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_pure_eta", "Proton Eta", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_pure_dcaz", "Proton DCAz", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + registryData.add("proton_pure_dcaxy", "Proton DCAxy", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); registryData.add("z_vtx", "Z-Vertex Distribution", HistType::kTH1F, {{200, -20.0, 20.0, "Z-Vertex (cm)"}}); registryData.add("tracks_in_jets", "Number of Tracks Inside Jets", HistType::kTH1F, {{100, 0, 100, "N_{tracks}"}}); @@ -156,48 +156,48 @@ struct PIDHadronsInJets { registryData.add("proton_jet_tpc", "TPC Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_jet_tof", "TOF Proton PID in Jets", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("pion_jet_pt", "Pion pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("pion_jet_eta", "Pion Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("pion_jet_dcaxy", "Pion DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("pion_jet_dcaz", "Pion DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("pion_jet_pt", "Pion pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_jet_eta", "Pion Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_jet_dcaxy", "Pion DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("pion_jet_dcaz", "Pion DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); - registryData.add("kaon_jet_pt", "Kaon pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("kaon_jet_eta", "Kaon Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("kaon_jet_dcaxy", "Kaon DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("kaon_jet_dcaz", "Kaon DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("kaon_jet_pt", "Kaon pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_jet_eta", "Kaon Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_jet_dcaxy", "Kaon DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("kaon_jet_dcaz", "Kaon DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); - registryData.add("proton_jet_pt", "Proton pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("proton_jet_eta", "Proton Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("proton_jet_dcaxy", "Proton DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("proton_jet_dcaz", "Proton DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("proton_jet_pt", "Proton pT in Jets", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_jet_eta", "Proton Eta in Jets", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_jet_dcaxy", "Proton DCAxy in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("proton_jet_dcaz", "Proton DCAz in Jets", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); registryData.add("pion_ue_tpc", "TPC Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("pion_ue_tof", "TOF Pion PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("pion_ue_pt", "Pion pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("pion_ue_eta", "Pion Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("pion_ue_dcaxy", "Pion DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("pion_ue_dcaz", "Pion DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("pion_ue_pt", "Pion pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("pion_ue_eta", "Pion Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("pion_ue_dcaxy", "Pion DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("pion_ue_dcaz", "Pion DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); registryData.add("kaon_ue_tpc", "TPC Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("kaon_ue_tof", "TOF Kaon PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("kaon_ue_pt", "Kaon pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("kaon_ue_eta", "Kaon Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("kaon_ue_dcaxy", "Kaon DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("kaon_ue_dcaz", "Kaon DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); + registryData.add("kaon_ue_pt", "Kaon pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("kaon_ue_eta", "Kaon Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("kaon_ue_dcaxy", "Kaon DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("kaon_ue_dcaz", "Kaon DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); registryData.add("proton_ue_tpc", "TPC Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TPC}"}}); registryData.add("proton_ue_tof", "TOF Proton PID in UE", HistType::kTH2F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}, {200, -3.0, 3.0, "n#sigma_{TOF}"}}); - registryData.add("proton_ue_pt", "Proton pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("proton_ue_eta", "Proton Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta" }}); - registryData.add("proton_ue_dcaxy", "Proton DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)" }}); - registryData.add("proton_ue_dcaz", "Proton DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)" }}); - - registryData.add("mc_gen_pion_pt", "Generated Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_rec_pion_pt", "Reconstructed Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_gen_kaon_pt", "Generated Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_rec_kaon_pt", "Reconstructed Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_gen_proton_pt", "Generated Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_rec_proton_pt", "Reconstructed Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("proton_ue_pt", "Proton pT in UE", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("proton_ue_eta", "Proton Eta in UE", HistType::kTH1F, {{100, -1.0, 1.0, "#eta"}}); + registryData.add("proton_ue_dcaxy", "Proton DCAxy in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{xy} (cm)"}}); + registryData.add("proton_ue_dcaz", "Proton DCAz in UE", HistType::kTH1F, {{100, -0.1, 0.1, "DCA_{z} (cm)"}}); + + registryData.add("mc_gen_pion_pt", "Generated Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_pion_pt", "Reconstructed Pion pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_gen_kaon_pt", "Generated Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_kaon_pt", "Reconstructed Kaon pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_gen_proton_pt", "Generated Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_rec_proton_pt", "Reconstructed Proton pT", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); registryData.add("rec_pion_all", "All identified pions", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); registryData.add("rec_kaon_all", "All identified kaons", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T}"}}); @@ -206,9 +206,9 @@ struct PIDHadronsInJets { registryData.add("contamination_matrix_kaon", "Kaon Contamination Matrix;True absolute PDG;Rec #it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{3000, -0.5, 2999.5}, {120, 0.0, 4.0}}); registryData.add("contamination_matrix_proton", "Proton Contamination Matrix;True absolute PDG;Rec #it{p}_{T} (GeV/#it{c})", HistType::kTH2F, {{3000, -0.5, 2999.5}, {120, 0.0, 4.0}}); - registryData.add("mc_sec_pion_pt", "Reconstructed Pion Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_sec_kaon_pt", "Reconstructed Kaon Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); - registryData.add("mc_sec_proton_pt", "Reconstructed Proton Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})" }}); + registryData.add("mc_sec_pion_pt", "Reconstructed Pion Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_sec_kaon_pt", "Reconstructed Kaon Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); + registryData.add("mc_sec_proton_pt", "Reconstructed Proton Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); } void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) @@ -217,16 +217,24 @@ struct PIDHadronsInJets { double px2 = px * px, py2 = py * py, pz2 = pz * pz; double pz4 = pz2 * pz2; - if (px == 0 && py == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; } + if (px == 0 && py == 0) { + u1.SetXYZ(0, 0, 0); + u2.SetXYZ(0, 0, 0); + return; + } if (px == 0 && py != 0) { double ux = std::sqrt(py2 - pz4 / py2); double uy = -pz2 / py; - u1.SetXYZ(ux, uy, pz); u2.SetXYZ(-ux, uy, pz); return; + u1.SetXYZ(ux, uy, pz); + u2.SetXYZ(-ux, uy, pz); + return; } if (py == 0 && px != 0) { double ux = -pz2 / px; double uy = std::sqrt(px2 - pz4 / px2); - u1.SetXYZ(ux, uy, pz); u2.SetXYZ(ux, -uy, pz); return; + u1.SetXYZ(ux, uy, pz); + u2.SetXYZ(ux, -uy, pz); + return; } double a = px2 + py2; @@ -234,7 +242,11 @@ struct PIDHadronsInJets { double c = pz4 - py2 * py2 - px2 * py2; double delta = b * b - 4.0 * a * c; - if (delta < 0 || a == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; } + if (delta < 0 || a == 0) { + u1.SetXYZ(0, 0, 0); + u2.SetXYZ(0, 0, 0); + return; + } double u1x = (-b + std::sqrt(delta)) / (2.0 * a); u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz); double u2x = (-b - std::sqrt(delta)) / (2.0 * a); @@ -248,8 +260,10 @@ struct PIDHadronsInJets { double phi2 = TVector2::Phi_0_2pi(a2); double diff = std::abs(phi1 - phi2); - if (diff <= PI) deltaPhi = diff; - if (diff > PI) deltaPhi = TwoPI - diff; + if (diff <= PI) + deltaPhi = diff; + if (diff > PI) + deltaPhi = TwoPI - diff; return deltaPhi; } @@ -271,52 +285,78 @@ struct PIDHadronsInJets { static constexpr double DcaxyMaxTrackPar2 = 1.1; static constexpr double DcazMaxTrack = 2.0; - if (!track.hasITS() || !track.hasTPC()) return false; - if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false; - if (track.tpcNClsCrossedRows() < MinTpcCr) return false; - if (track.tpcChi2NCl() > MaxChi2Tpc) return false; - if (track.itsChi2NCl() > MaxChi2Its) return false; - if (std::abs(track.eta()) > maxEta) return false; - if (track.pt() < cfgMinPtTrack) return false; - if (std::abs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) return false; - if (std::abs(track.dcaZ()) > DcazMaxTrack) return false; + if (!track.hasITS() || !track.hasTPC()) + return false; + if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) + return false; + if (track.tpcNClsCrossedRows() < MinTpcCr) + return false; + if (track.tpcChi2NCl() > MaxChi2Tpc) + return false; + if (track.itsChi2NCl() > MaxChi2Its) + return false; + if (std::abs(track.eta()) > maxEta) + return false; + if (track.pt() < cfgMinPtTrack) + return false; + if (std::abs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) + return false; + if (std::abs(track.dcaZ()) > DcazMaxTrack) + return false; return true; } template bool passedTrackSelection(const PionTrack& track) { - if (requirePvContributor && !(track.isPVContributor())) return false; - if (!track.hasITS() || !track.hasTPC()) return false; - if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false; - if (track.itsNCls() < minItsNclusters) return false; - if (track.tpcNClsCrossedRows() < minTpcNcrossedRows) return false; - if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) return false; - if (track.itsChi2NCl() > maxChiSquareIts) return false; - if (track.eta() < minEta || track.eta() > maxEta) return false; - if (track.pt() < minPt || track.pt() > maxPt) return false; + if (requirePvContributor && !(track.isPVContributor())) + return false; + if (!track.hasITS() || !track.hasTPC()) + return false; + if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) + return false; + if (track.itsNCls() < minItsNclusters) + return false; + if (track.tpcNClsCrossedRows() < minTpcNcrossedRows) + return false; + if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) + return false; + if (track.itsChi2NCl() > maxChiSquareIts) + return false; + if (track.eta() < minEta || track.eta() > maxEta) + return false; + if (track.pt() < minPt || track.pt() > maxPt) + return false; return true; } - void processForJets(SelectedCollisions::iterator const& collision, HadronTracks const& tracks) { registryData.fill(HIST("n_events_raw"), 1); - if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) return; - if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; - if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; - if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; - if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) return; - if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; - if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) + return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) + return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) + return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) + return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) + return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) + return; registryData.fill(HIST("n_events"), 1); registryData.fill(HIST("z_vtx"), collision.posZ()); for (auto const& track : tracks) { - if (!passedTrackSelection(track)) continue; - if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; double pt = track.pt(); double eta = track.eta(); @@ -383,21 +423,24 @@ struct PIDHadronsInJets { std::vector fjParticles; for (auto const& track : tracks) { id++; - if (!passedTrackSelectionForJetReconstruction(track)) continue; + if (!passedTrackSelectionForJetReconstruction(track)) + continue; fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged)); fourMomentum.set_user_index(id); fjParticles.emplace_back(fourMomentum); } - if (fjParticles.empty()) return; + if (fjParticles.empty()) + return; fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet); fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0)); fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef); std::vector jets = fastjet::sorted_by_pt(cs.inclusive_jets()); - if (jets.empty()) return; + if (jets.empty()) + return; auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet); @@ -410,27 +453,34 @@ struct PIDHadronsInJets { registryData.fill(HIST("jet_area"), jet.area()); registryData.fill(HIST("jet_n_constituents"), jet.constituents().size()); - if (!isppRefAnalysis && ((std::abs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue; - if (isppRefAnalysis && std::abs(jet.eta()) > cfgEtaJetMax) continue; + if (!isppRefAnalysis && ((std::abs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) + continue; + if (isppRefAnalysis && std::abs(jet.eta()) > cfgEtaJetMax) + continue; auto jetForSub = jet; fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp); registryData.fill(HIST("jet_pt_subtracted"), jetMinusBkg.pt()); registryData.fill(HIST("jet_pt_raw_vs_sub"), jet.pt(), jetMinusBkg.pt()); - if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) continue; - if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) continue; + if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) + continue; + if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) + continue; double normalizedJetArea = jet.area() / (PI * rJet * rJet); - if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue; - if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue; + if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) + continue; + if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) + continue; double coneRadius = std::sqrt(jet.area() / PI); TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0); getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2); - if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) continue; + if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) + continue; std::vector jetConstituents = jet.constituents(); @@ -441,8 +491,10 @@ struct PIDHadronsInJets { tracksInJetsSet.insert(trackIdx); - if (!passedTrackSelection(track)) continue; - if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; double pt = track.pt(); double eta = track.eta(); @@ -507,8 +559,10 @@ struct PIDHadronsInJets { for (auto const& track : tracks) { - if (!passedTrackSelection(track)) continue; - if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi()); @@ -523,7 +577,8 @@ struct PIDHadronsInJets { maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet; } - if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; + if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) + continue; double pt = track.pt(); double eta = track.eta(); @@ -594,21 +649,30 @@ struct PIDHadronsInJets { registryData.fill(HIST("tracks_outside_jets"), nTracksOut); } - PROCESS_SWITCH(PIDHadronsInJets, processForJets, "Pid Analysis in and outside jets", true); + PROCESS_SWITCH(jetHadronsPid, processForJets, "Pid Analysis in and outside jets", true); void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks) { - if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) return; - if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return; - if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return; - if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return; - if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) return; - if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return; - if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; + if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) + return; + if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) + return; + if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) + return; + if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) + return; + if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) + return; + if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) + return; + if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) + return; for (auto const& mcpart : mcParticles) { - if (!mcpart.isPhysicalPrimary()) continue; - if (std::abs(mcpart.eta()) > 0.8) continue; + if (!mcpart.isPhysicalPrimary()) + continue; + if (std::abs(mcpart.eta()) > 0.8) + continue; int pdg = std::abs(mcpart.pdgCode()); double pt = mcpart.pt(); @@ -625,8 +689,10 @@ struct PIDHadronsInJets { const double ptThreshold = 0.8; for (auto const& track : tracks) { - if (!passedTrackSelection(track)) continue; - if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) continue; + if (!passedTrackSelection(track)) + continue; + if (std::abs(track.dcaXY()) > maxDcaxy || std::abs(track.dcaZ()) > maxDcaz) + continue; double pt = track.pt(); @@ -639,17 +705,19 @@ struct PIDHadronsInJets { bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); - if (!track.has_mcParticle()) continue; + if (!track.has_mcParticle()) + continue; auto const& trueParticle = track.mcParticle(); int pdg = std::abs(trueParticle.pdgCode()); bool isPrimary = trueParticle.isPhysicalPrimary(); if (passTpcPi && (pt < ptThreshold || passTofPi)) { - registryData.fill(HIST("rec_pion_all"), pt); + registryData.fill(HIST("rec_pion_all"), pt); registryData.fill(HIST("contamination_matrix_pion"), pdg, pt); if (isPrimary) { - if (pdg == 211) registryData.fill(HIST("mc_rec_pion_pt"), pt); + if (pdg == 211) + registryData.fill(HIST("mc_rec_pion_pt"), pt); } else { registryData.fill(HIST("mc_sec_pion_pt"), pt); } @@ -659,7 +727,8 @@ struct PIDHadronsInJets { registryData.fill(HIST("rec_kaon_all"), pt); registryData.fill(HIST("contamination_matrix_kaon"), pdg, pt); if (isPrimary) { - if (pdg == 321) registryData.fill(HIST("mc_rec_kaon_pt"), pt); + if (pdg == 321) + registryData.fill(HIST("mc_rec_kaon_pt"), pt); } else { registryData.fill(HIST("mc_sec_kaon_pt"), pt); } @@ -669,18 +738,18 @@ struct PIDHadronsInJets { registryData.fill(HIST("rec_proton_all"), pt); registryData.fill(HIST("contamination_matrix_proton"), pdg, pt); if (isPrimary) { - if (pdg == 2212) registryData.fill(HIST("mc_rec_proton_pt"), pt); + if (pdg == 2212) + registryData.fill(HIST("mc_rec_proton_pt"), pt); } else { registryData.fill(HIST("mc_sec_proton_pt"), pt); } } } } - PROCESS_SWITCH(PIDHadronsInJets, processMC, "Run on Monte Carlo", false); - + PROCESS_SWITCH(jetHadronsPid, processMC, "Run on Monte Carlo", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { - return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} \ No newline at end of file + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} From 371c30d3384d013d124c1dd5f6031235d8086a8f Mon Sep 17 00:00:00 2001 From: Leonard Date: Wed, 29 Apr 2026 22:09:37 +0200 Subject: [PATCH 11/11] Changing code in jetHadronsPid.cxx according with comments on pull request --- PWGJE/Tasks/jetHadronsPid.cxx | 220 +++++++++++++++++----------------- 1 file changed, 110 insertions(+), 110 deletions(-) diff --git a/PWGJE/Tasks/jetHadronsPid.cxx b/PWGJE/Tasks/jetHadronsPid.cxx index 0d5263e67ce..80076646d9e 100644 --- a/PWGJE/Tasks/jetHadronsPid.cxx +++ b/PWGJE/Tasks/jetHadronsPid.cxx @@ -9,7 +9,7 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// \file jetHadronsPid.cxx +/// \file jetHadronsPida.cxx /// \brief Analysis of hadrons in jets /// \author Aleksandra Mulewicz, WUT Warsaw, aleksandra.mulewicz@cern.ch /// \author Leonard Lorenc, WUT Warsaw, leonard.lorenc@cern.ch @@ -21,6 +21,7 @@ #include "PWGJE/Core/JetUtilities.h" #include "PWGJE/DataModel/Jet.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseITS.h" @@ -70,7 +71,6 @@ struct jetHadronsPid { HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; Configurable isppRefAnalysis{"isppRefAnalysis", false, "Is ppRef analysis"}; - Configurable cfgAreaFrac{"cfgAreaFrac", 0.6, "fraction of jet area"}; Configurable cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"}; Configurable cfgMinPtTrack{"cfgMinPtTrack", 0.1, "minimum pt of tracks for jet reconstruction"}; @@ -86,7 +86,7 @@ struct jetHadronsPid { Configurable rJet{"rJet", 0.4, "Jet resolution parameter R"}; Configurable zVtx{"zVtx", 10.0, "Maximum zVertex"}; Configurable applyAreaCut{"applyAreaCut", true, "apply area cut"}; - Configurable maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"}; + Configurable minNormalizedJetArea{"minNormalizedJetArea", 0.6, "Minimum normalized area cut to reject fake jets"}; Configurable deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"}; Configurable requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"}; @@ -101,6 +101,7 @@ struct jetHadronsPid { Configurable maxEta{"maxEta", +0.8, "maximum eta"}; Configurable maxDcaxy{"maxDcaxy", 0.2, "Maximum DCAxy"}; Configurable maxDcaz{"maxDcaz", 0.1, "Maximum DCAz"}; + Configurable maxNSigmaPid{"maxNSigmaPid", 3.0, "Maximum nSigma for TPC and TOF PID"}; Configurable setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"}; @@ -211,67 +212,64 @@ struct jetHadronsPid { registryData.add("mc_sec_proton_pt", "Reconstructed Proton Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}}); } + // void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) + // { + // double px = p.X(), py = p.Y(), pz = p.Z(); + // double px2 = px * px, py2 = py * py, pz2 = pz * pz; + // double pz4 = pz2 * pz2; + + // if (px == 0 && py == 0) { + // u1.SetXYZ(0, 0, 0); + // u2.SetXYZ(0, 0, 0); + // return; + // } + // if (px == 0 && py != 0) { + // double ux = std::sqrt(py2 - pz4 / py2); + // double uy = -pz2 / py; + // u1.SetXYZ(ux, uy, pz); + // u2.SetXYZ(-ux, uy, pz); + // return; + // } + // if (py == 0 && px != 0) { + // double ux = -pz2 / px; + // double uy = std::sqrt(px2 - pz4 / px2); + // u1.SetXYZ(ux, uy, pz); + // u2.SetXYZ(ux, -uy, pz); + // return; + // } + + // double a = px2 + py2; + // double b = 2.0 * px * pz2; + // double c = pz4 - py2 * py2 - px2 * py2; + // double delta = b * b - 4.0 * a * c; + + // if (delta < 0 || a == 0) { + // u1.SetXYZ(0, 0, 0); + // u2.SetXYZ(0, 0, 0); + // return; + // } + // double u1x = (-b + std::sqrt(delta)) / (2.0 * a); + // u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz); + // double u2x = (-b - std::sqrt(delta)) / (2.0 * a); + // u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz); + // } + void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2) { - double px = p.X(), py = p.Y(), pz = p.Z(); - double px2 = px * px, py2 = py * py, pz2 = pz * pz; - double pz4 = pz2 * pz2; - - if (px == 0 && py == 0) { - u1.SetXYZ(0, 0, 0); - u2.SetXYZ(0, 0, 0); - return; - } - if (px == 0 && py != 0) { - double ux = std::sqrt(py2 - pz4 / py2); - double uy = -pz2 / py; - u1.SetXYZ(ux, uy, pz); - u2.SetXYZ(-ux, uy, pz); - return; - } - if (py == 0 && px != 0) { - double ux = -pz2 / px; - double uy = std::sqrt(px2 - pz4 / px2); - u1.SetXYZ(ux, uy, pz); - u2.SetXYZ(ux, -uy, pz); - return; - } - - double a = px2 + py2; - double b = 2.0 * px * pz2; - double c = pz4 - py2 * py2 - px2 * py2; - double delta = b * b - 4.0 * a * c; - - if (delta < 0 || a == 0) { + if (p.Mag2() < 1e-9) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; } - double u1x = (-b + std::sqrt(delta)) / (2.0 * a); - u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz); - double u2x = (-b - std::sqrt(delta)) / (2.0 * a); - u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz); - } - - double getDeltaPhi(double a1, double a2) - { - double deltaPhi(0); - double phi1 = TVector2::Phi_0_2pi(a1); - double phi2 = TVector2::Phi_0_2pi(a2); - double diff = std::abs(phi1 - phi2); - - if (diff <= PI) - deltaPhi = diff; - if (diff > PI) - deltaPhi = TwoPI - diff; - return deltaPhi; + u1 = p.Orthogonal(); + u2 = p.Cross(u1); } template - bool hasITSHit(const TrackIts& track, int layer) + bool hasITSLayerHit(const TrackIts& track, int layer) { int ibit = layer - 1; - return (track.itsClusterMap() & (1 << ibit)); + return (track.itsClusterMap() & (1 << ibit)) != 0; } template @@ -287,7 +285,7 @@ struct jetHadronsPid { if (!track.hasITS() || !track.hasTPC()) return false; - if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) + if ((!hasITSLayerHit(track, 1)) && (!hasITSLayerHit(track, 2)) && (!hasITSLayerHit(track, 3))) return false; if (track.tpcNClsCrossedRows() < MinTpcCr) return false; @@ -313,7 +311,7 @@ struct jetHadronsPid { return false; if (!track.hasITS() || !track.hasTPC()) return false; - if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) + if ((!hasITSLayerHit(track, 1)) && (!hasITSLayerHit(track, 2)) && (!hasITSLayerHit(track, 3))) return false; if (track.itsNCls() < minItsNclusters) return false; @@ -363,13 +361,13 @@ struct jetHadronsPid { double dcaxy = track.dcaXY(); double dcaz = track.dcaZ(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); - bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); - bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); - bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); - bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); - bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); const double ptThreshold = 0.8; @@ -419,15 +417,15 @@ struct jetHadronsPid { } } - int id(-1); + // int id(-1); std::vector fjParticles; for (auto const& track : tracks) { - id++; + // id++; if (!passedTrackSelectionForJetReconstruction(track)) continue; fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged)); - fourMomentum.set_user_index(id); + fourMomentum.set_user_index(track.index()); fjParticles.emplace_back(fourMomentum); } @@ -469,10 +467,9 @@ struct jetHadronsPid { continue; double normalizedJetArea = jet.area() / (PI * rJet * rJet); - if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) - continue; - if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) + if (applyAreaCut && normalizedJetArea < minNormalizedJetArea) { continue; + } double coneRadius = std::sqrt(jet.area() / PI); TVector3 jetAxis(jet.px(), jet.py(), jet.pz()); @@ -501,13 +498,13 @@ struct jetHadronsPid { double dcaxy = track.dcaXY(); double dcaz = track.dcaZ(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); - bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); - bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); - bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); - bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); - bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); const double ptThreshold = 0.8; @@ -565,17 +562,14 @@ struct jetHadronsPid { continue; double deltaEtaUe1 = track.eta() - ueAxis1.Eta(); - double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi()); + double deltaPhiUe1 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis1.Phi())); double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1); double deltaEtaUe2 = track.eta() - ueAxis2.Eta(); - double deltaPhiUe2 = getDeltaPhi(track.phi(), ueAxis2.Phi()); + double deltaPhiUe2 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis2.Phi())); double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2); - double maxConeRadius = coneRadius; - if (applyAreaCut) { - maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet; - } + double maxConeRadius = rJet; if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue; @@ -585,13 +579,13 @@ struct jetHadronsPid { double dcaxy = track.dcaXY(); double dcaz = track.dcaZ(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); - bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); - bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); - bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); - bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); - bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); const double ptThreshold = 0.8; @@ -651,7 +645,7 @@ struct jetHadronsPid { PROCESS_SWITCH(jetHadronsPid, processForJets, "Pid Analysis in and outside jets", true); - void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks) + void processMC(SelectedCollisions::iterator const& collision, HadronTracksMC const& tracks) { if (!collision.sel8() || std::abs(collision.posZ()) > zVtx) return; @@ -668,24 +662,6 @@ struct jetHadronsPid { if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return; - for (auto const& mcpart : mcParticles) { - if (!mcpart.isPhysicalPrimary()) - continue; - if (std::abs(mcpart.eta()) > 0.8) - continue; - - int pdg = std::abs(mcpart.pdgCode()); - double pt = mcpart.pt(); - - if (pdg == 211) { - registryData.fill(HIST("mc_gen_pion_pt"), pt); - } else if (pdg == 321) { - registryData.fill(HIST("mc_gen_kaon_pt"), pt); - } else if (pdg == 2212) { - registryData.fill(HIST("mc_gen_proton_pt"), pt); - } - } - const double ptThreshold = 0.8; for (auto const& track : tracks) { @@ -696,14 +672,14 @@ struct jetHadronsPid { double pt = track.pt(); - bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0); - bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0); + bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid); + bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid); - bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0); - bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0); + bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid); + bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid); - bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0); - bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0); + bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid); + bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid); if (!track.has_mcParticle()) continue; @@ -747,6 +723,30 @@ struct jetHadronsPid { } } PROCESS_SWITCH(jetHadronsPid, processMC, "Run on Monte Carlo", false); + void processMCTruth(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles) + { + if (std::abs(mcCollision.posZ()) > zVtx) + return; + + for (auto const& mcpart : mcParticles) { + if (!mcpart.isPhysicalPrimary()) + continue; + if (std::abs(mcpart.eta()) > maxEta) + continue; + + int pdg = std::abs(mcpart.pdgCode()); + double pt = mcpart.pt(); + + if (pdg == 211) { + registryData.fill(HIST("mc_gen_pion_pt"), pt); + } else if (pdg == 321) { + registryData.fill(HIST("mc_gen_kaon_pt"), pt); + } else if (pdg == 2212) { + registryData.fill(HIST("mc_gen_proton_pt"), pt); + } + } + } + PROCESS_SWITCH(jetHadronsPid, processMCTruth, "Run on Monte Carlo (Pure Truth)", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)