Skip to content

Commit ec18b06

Browse files
MyFavoriteGitHubMyFavoriteGitHub
committed
Add jetHadronsPID and update CMakeLists
1 parent 3ebe857 commit ec18b06

2 files changed

Lines changed: 246 additions & 0 deletions

File tree

‎PWGJE/Tasks/CMakeLists.txt‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -403,4 +403,8 @@ if(FastJet_FOUND)
403403
SOURCES bjetCentMult.cxx
404404
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::PWGJECore O2Physics::AnalysisCore
405405
COMPONENT_NAME Analysis)
406+
o2physics_add_dpl_workflow(pid-hadrons-in-jets
407+
SOURCES jetHadronsPID.cxx
408+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGJECore FastJet::FastJet FastJet::Contrib O2Physics::EventFilteringUtils
409+
COMPONENT_NAME Analysis)
406410
endif()

‎PWGJE/Tasks/jetHadronsPID.cxx‎

Lines changed: 242 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,242 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
5+
#include "Framework/runDataProcessing.h"
6+
#include "PWGJE/Core/JetBkgSubUtils.h"
7+
#include "PWGJE/Core/JetUtilities.h"
8+
#include "PWGJE/DataModel/Jet.h"
9+
10+
#include "Common/Core/TrackSelection.h"
11+
#include "Common/DataModel/EventSelection.h"
12+
#include "Common/DataModel/PIDResponseITS.h"
13+
#include "Common/DataModel/TrackSelectionTables.h"
14+
15+
#include "Framework/AnalysisTask.h"
16+
#include "Framework/HistogramRegistry.h"
17+
#include "Framework/Logger.h"
18+
#include "ReconstructionDataFormats/DCA.h"
19+
#include "ReconstructionDataFormats/PID.h"
20+
#include "ReconstructionDataFormats/Track.h"
21+
22+
#include <fastjet/AreaDefinition.hh>
23+
#include <fastjet/ClusterSequence.hh>
24+
#include <fastjet/ClusterSequenceArea.hh>
25+
#include <fastjet/PseudoJet.hh>
26+
27+
#include <cmath>
28+
#include <vector>
29+
30+
using namespace o2;
31+
using namespace o2::framework;
32+
using namespace o2::aod;
33+
using namespace o2::constants::physics;
34+
using namespace o2::constants::math;
35+
36+
// Define convenient aliases for commonly used table joins
37+
using SelectedCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms>;
38+
39+
// CRITICAL FIX: Replaced Proton/Deuteron/Helium PID tables with Pion PID tables
40+
using PionTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCPi, aod::pidTOFPi>;
41+
42+
struct PIDHadronsInJets {
43+
44+
// Histogram registry for data
45+
HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
46+
47+
// Parameters for ppRef analysis
48+
Configurable<bool> isppRefAnalysis{"isppRefAnalysis", false, "Is ppRef analysis"};
49+
Configurable<double> cfgAreaFrac{"cfgAreaFrac", 0.6, "fraction of jet area"};
50+
Configurable<double> cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"};
51+
Configurable<double> cfgMinPtTrack{"cfgMinPtTrack", 0.1, "minimum pt of tracks for jet reconstruction"};
52+
53+
// Event selection criteria
54+
Configurable<bool> rejectITSROFBorder{"rejectITSROFBorder", true, "Reject events near the ITS ROF border"};
55+
Configurable<bool> rejectTFBorder{"rejectTFBorder", true, "Reject events near the TF border"};
56+
Configurable<bool> requireVtxITSTPC{"requireVtxITSTPC", true, "Require at least one ITS-TPC matched track"};
57+
Configurable<bool> rejectSameBunchPileup{"rejectSameBunchPileup", true, "Reject events with same-bunch pileup collisions"};
58+
Configurable<bool> requireIsGoodZvtxFT0VsPV{"requireIsGoodZvtxFT0VsPV", true, "Require consistent FT0 vs PV z-vertex"};
59+
Configurable<bool> requireIsVertexTOFmatched{"requireIsVertexTOFmatched", false, "Require at least one vertex track matched to TOF"};
60+
61+
// Jet selection parameters
62+
Configurable<double> minJetPt{"minJetPt", 10.0, "Minimum pt of the jet after bkg subtraction"};
63+
Configurable<double> maxJetPt{"maxJetPt", 1e+06, "Maximum pt of the jet after bkg subtraction"};
64+
Configurable<double> rJet{"rJet", 0.4, "Jet resolution parameter R"};
65+
Configurable<double> zVtx{"zVtx", 10.0, "Maximum zVertex"};
66+
Configurable<bool> applyAreaCut{"applyAreaCut", true, "apply area cut"};
67+
Configurable<double> maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"};
68+
Configurable<double> deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"};
69+
70+
// Track quality parameters
71+
Configurable<bool> requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"};
72+
Configurable<int> minItsNclusters{"minItsNclusters", 5, "minimum number of ITS clusters"};
73+
Configurable<int> minTpcNcrossedRows{"minTpcNcrossedRows", 100, "minimum number of TPC crossed pad rows"};
74+
Configurable<double> minChiSquareTpc{"minChiSquareTpc", 0.0, "minimum TPC chi^2/Ncls"};
75+
Configurable<double> maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"};
76+
Configurable<double> maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"};
77+
Configurable<double> minPt{"minPt", 0.3, "minimum pt of the tracks"};
78+
Configurable<double> minEta{"minEta", -0.8, "minimum eta"};
79+
Configurable<double> maxEta{"maxEta", +0.8, "maximum eta"};
80+
Configurable<double> maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"};
81+
Configurable<double> maxDcaz{"maxDcaz", 0.05, "Maximum DCAz"};
82+
83+
Configurable<bool> setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"};
84+
85+
// Utility object for jet background subtraction methods
86+
JetBkgSubUtils backgroundSub;
87+
88+
// Initialize ITS PID Response object
89+
o2::aod::ITSResponse itsResponse;
90+
91+
void init(InitContext const&)
92+
{
93+
if (setMCDefaultItsParams) {
94+
itsResponse.setMCDefaultParameters();
95+
}
96+
97+
// pid of pions
98+
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}"}});
99+
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}"}});
100+
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}"}});
101+
}
102+
103+
// ITS hit helper
104+
template <typename TrackIts>
105+
bool hasITSHit(const TrackIts& track, int layer)
106+
{
107+
int ibit = layer - 1;
108+
return (track.itsClusterMap() & (1 << ibit));
109+
}
110+
111+
// Single-track selection for jet reconstruction
112+
template <typename JetTrack>
113+
bool passedTrackSelectionForJetReconstruction(const JetTrack& track)
114+
{
115+
static constexpr int MinTpcCr = 70;
116+
static constexpr double MaxChi2Tpc = 4.0;
117+
static constexpr double MaxChi2Its = 36.0;
118+
static constexpr double DcaxyMaxTrackPar0 = 0.0105;
119+
static constexpr double DcaxyMaxTrackPar1 = 0.035;
120+
static constexpr double DcaxyMaxTrackPar2 = 1.1;
121+
static constexpr double DcazMaxTrack = 2.0;
122+
123+
if (!track.hasITS() || !track.hasTPC()) return false;
124+
if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false;
125+
if (track.tpcNClsCrossedRows() < MinTpcCr) return false;
126+
if (track.tpcChi2NCl() > MaxChi2Tpc) return false;
127+
if (track.itsChi2NCl() > MaxChi2Its) return false;
128+
if (std::fabs(track.eta()) > maxEta) return false;
129+
if (track.pt() < cfgMinPtTrack) return false;
130+
if (std::fabs(track.dcaXY()) > (DcaxyMaxTrackPar0 + DcaxyMaxTrackPar1 / std::pow(track.pt(), DcaxyMaxTrackPar2))) return false;
131+
if (std::fabs(track.dcaZ()) > DcazMaxTrack) return false;
132+
return true;
133+
}
134+
135+
// Single-track selection for constituents
136+
template <typename PionTrack>
137+
bool passedTrackSelection(const PionTrack& track)
138+
{
139+
if (requirePvContributor && !(track.isPVContributor())) return false;
140+
if (!track.hasITS() || !track.hasTPC()) return false;
141+
if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3))) return false;
142+
if (track.itsNCls() < minItsNclusters) return false;
143+
if (track.tpcNClsCrossedRows() < minTpcNcrossedRows) return false;
144+
if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) return false;
145+
if (track.itsChi2NCl() > maxChiSquareIts) return false;
146+
if (track.eta() < minEta || track.eta() > maxEta) return false;
147+
if (track.pt() < minPt) return false;
148+
return true;
149+
}
150+
151+
// Process Data
152+
void processData(SelectedCollisions::iterator const& collision, PionTracks const& tracks)
153+
{
154+
// Apply standard event selection
155+
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return;
156+
if (rejectITSROFBorder && !collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) return;
157+
if (rejectTFBorder && !collision.selection_bit(o2::aod::evsel::kNoTimeFrameBorder)) return;
158+
if (requireVtxITSTPC && !collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) return;
159+
if (rejectSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) return;
160+
if (requireIsGoodZvtxFT0VsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) return;
161+
if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched)) return;
162+
163+
// Build FastJet particles
164+
int id(-1);
165+
std::vector<fastjet::PseudoJet> fjParticles;
166+
for (auto const& track : tracks) {
167+
id++;
168+
if (!passedTrackSelectionForJetReconstruction(track)) continue;
169+
170+
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
171+
fourMomentum.set_user_index(id);
172+
fjParticles.emplace_back(fourMomentum);
173+
}
174+
175+
if (fjParticles.empty()) return;
176+
177+
// Cluster particles
178+
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
179+
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
180+
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
181+
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
182+
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
183+
184+
// Loop over reconstructed jets
185+
for (const auto& jet : jets) {
186+
187+
if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue;
188+
if (isppRefAnalysis && std::fabs(jet.eta()) > cfgEtaJetMax) continue;
189+
190+
auto jetForSub = jet;
191+
fastjet::PseudoJet jetMinusBkg = backgroundSub.doRhoAreaSub(jetForSub, rhoPerp, rhoMPerp);
192+
193+
if (isppRefAnalysis && (jet.pt() < minJetPt || jet.pt() > maxJetPt)) continue;
194+
if (!isppRefAnalysis && (jetMinusBkg.pt() < minJetPt || jetMinusBkg.pt() > maxJetPt)) continue;
195+
196+
double normalizedJetArea = jet.area() / (PI * rJet * rJet);
197+
if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue;
198+
if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue;
199+
200+
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
201+
202+
// loop to fill historgrams
203+
for (const auto& particle : jetConstituents) {
204+
205+
auto const& track = tracks.iteratorAt(particle.user_index());
206+
207+
// Constituent Track Selection (includes DCA checks)
208+
if (!passedTrackSelection(track)) continue;
209+
if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue;
210+
211+
double pt = track.pt();
212+
213+
// Check Pion PID (+/- 3 Sigma)
214+
double nsigmaTPCPi = track.tpcNSigmaPi();
215+
double nSigmaITSPi = static_cast<double>(itsResponse.nSigmaITS<o2::track::PID::Pion>(track));
216+
217+
// Fill TPC
218+
if (std::abs(nsigmaTPCPi) <= 3.0) {
219+
registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi);
220+
}
221+
222+
// Fill TOF
223+
if (track.hasTOF()) {
224+
double nsigmaTOFPi = track.tofNSigmaPi();
225+
if (std::abs(nsigmaTOFPi) <= 3.0) {
226+
registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi);
227+
}
228+
}
229+
230+
// Fill ITS
231+
if (std::abs(nSigmaITSPi) <= 3.0) {
232+
registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi);
233+
}
234+
}
235+
}
236+
}
237+
};
238+
239+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
240+
{
241+
return WorkflowSpec{adaptAnalysisTask<PIDHadronsInJets>(cfgc)};
242+
}

0 commit comments

Comments
 (0)