Skip to content

Commit 267a9fa

Browse files
hzalewskiMyFavoriteGitHub
hzalewski
authored and
MyFavoriteGitHub
committed
add historams for hadrons in jet and ue
1 parent 4464823 commit 267a9fa

1 file changed

Lines changed: 139 additions & 29 deletions

File tree

‎PWGJE/Tasks/jetHadronsPID.cxx‎

Lines changed: 139 additions & 29 deletions
Original file line numberDiff line numberDiff line change
@@ -24,6 +24,9 @@
2424
#include <fastjet/ClusterSequenceArea.hh>
2525
#include <fastjet/PseudoJet.hh>
2626

27+
#include <TVector2.h>
28+
#include <TVector3.h>
29+
2730
#include <cmath>
2831
#include <vector>
2932

@@ -37,7 +40,10 @@ using namespace o2::constants::math;
3740
using SelectedCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms>;
3841

3942

40-
using PionTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::pidTPCPi, aod::pidTOFPi>;
43+
using HadronTracks = soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA,
44+
aod::pidTPCFullPi, aod::pidTOFFullPi,
45+
aod::pidTPCFullKa, aod::pidTOFFullKa,
46+
aod::pidTPCFullPr, aod::pidTOFFullPr>;
4147

4248
struct PIDHadronsInJets {
4349

@@ -75,6 +81,7 @@ struct PIDHadronsInJets {
7581
Configurable<double> maxChiSquareTpc{"maxChiSquareTpc", 4.0, "maximum TPC chi^2/Ncls"};
7682
Configurable<double> maxChiSquareIts{"maxChiSquareIts", 36.0, "maximum ITS chi^2/Ncls"};
7783
Configurable<double> minPt{"minPt", 0.3, "minimum pt of the tracks"};
84+
Configurable<double> maxPt{"maxPt", 4.0, "maximum pt of the tracks for PID analysis"}; // pt cut at 4.0
7885
Configurable<double> minEta{"minEta", -0.8, "minimum eta"};
7986
Configurable<double> maxEta{"maxEta", +0.8, "maximum eta"};
8087
Configurable<double> maxDcaxy{"maxDcaxy", 0.05, "Maximum DCAxy"};
@@ -94,10 +101,64 @@ struct PIDHadronsInJets {
94101
itsResponse.setMCDefaultParameters();
95102
}
96103

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}"}});
104+
// Pions, Kaons, Protons in Jet
105+
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}"}});
106+
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}"}});
107+
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}"}});
108+
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}"}});
109+
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}"}});
110+
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}"}});
111+
112+
// Pions, Kaons, Protons Underlying Event
113+
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}"}});
114+
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}"}});
115+
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}"}});
116+
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}"}});
117+
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}"}});
118+
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}"}});
119+
}
120+
121+
// Compute two transverse directions orthogonal to vector p
122+
void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2)
123+
{
124+
double px = p.X(), py = p.Y(), pz = p.Z();
125+
double px2 = px * px, py2 = py * py, pz2 = pz * pz;
126+
double pz4 = pz2 * pz2;
127+
128+
if (px == 0 && py == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; }
129+
if (px == 0 && py != 0) {
130+
double ux = std::sqrt(py2 - pz4 / py2);
131+
double uy = -pz2 / py;
132+
u1.SetXYZ(ux, uy, pz); u2.SetXYZ(-ux, uy, pz); return;
133+
}
134+
if (py == 0 && px != 0) {
135+
double ux = -pz2 / px;
136+
double uy = std::sqrt(px2 - pz4 / px2);
137+
u1.SetXYZ(ux, uy, pz); u2.SetXYZ(ux, -uy, pz); return;
138+
}
139+
140+
double a = px2 + py2;
141+
double b = 2.0 * px * pz2;
142+
double c = pz4 - py2 * py2 - px2 * py2;
143+
double delta = b * b - 4.0 * a * c;
144+
145+
if (delta < 0 || a == 0) { u1.SetXYZ(0, 0, 0); u2.SetXYZ(0, 0, 0); return; }
146+
double u1x = (-b + std::sqrt(delta)) / (2.0 * a);
147+
u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz);
148+
double u2x = (-b - std::sqrt(delta)) / (2.0 * a);
149+
u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz);
150+
}
151+
152+
double getDeltaPhi(double a1, double a2)
153+
{
154+
double deltaPhi(0);
155+
double phi1 = TVector2::Phi_0_2pi(a1);
156+
double phi2 = TVector2::Phi_0_2pi(a2);
157+
double diff = std::fabs(phi1 - phi2);
158+
159+
if (diff <= PI) deltaPhi = diff;
160+
if (diff > PI) deltaPhi = TwoPI - diff;
161+
return deltaPhi;
101162
}
102163

103164
// ITS hit helper
@@ -144,12 +205,12 @@ struct PIDHadronsInJets {
144205
if (track.tpcChi2NCl() < minChiSquareTpc || track.tpcChi2NCl() > maxChiSquareTpc) return false;
145206
if (track.itsChi2NCl() > maxChiSquareIts) return false;
146207
if (track.eta() < minEta || track.eta() > maxEta) return false;
147-
if (track.pt() < minPt) return false;
208+
if (track.pt() < minPt || track.pt() > maxPt) return false;
148209
return true;
149210
}
150211

151212
// Process Data
152-
void process(SelectedCollisions::iterator const& collision, PionTracks const& tracks)
213+
void process(SelectedCollisions::iterator const& collision, HadronTracks const& tracks)
153214
{
154215
// Apply standard event selection
155216
if (!collision.sel8() || std::fabs(collision.posZ()) > zVtx) return;
@@ -174,14 +235,16 @@ struct PIDHadronsInJets {
174235

175236
if (fjParticles.empty()) return;
176237

177-
// Cluster particles
238+
178239
fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, rJet);
179240
fastjet::AreaDefinition areaDef(fastjet::active_area, fastjet::GhostedAreaSpec(1.0));
180241
fastjet::ClusterSequenceArea cs(fjParticles, jetDef, areaDef);
181242
std::vector<fastjet::PseudoJet> jets = fastjet::sorted_by_pt(cs.inclusive_jets());
243+
244+
if (jets.empty()) return;
245+
182246
auto [rhoPerp, rhoMPerp] = jetutilities::estimateRhoPerpCone(fjParticles, jets[0], rJet);
183247

184-
// Loop over reconstructed jets
185248
for (const auto& jet : jets) {
186249

187250
if (!isppRefAnalysis && ((std::fabs(jet.eta()) + rJet) > (maxEta - deltaEtaEdge))) continue;
@@ -197,41 +260,88 @@ struct PIDHadronsInJets {
197260
if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea) continue;
198261
if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet)) continue;
199262

200-
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
201263

202-
// loop to fill historgrams
264+
double coneRadius = std::sqrt(jet.area() / PI);
265+
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
266+
TVector3 ueAxis1(0, 0, 0), ueAxis2(0, 0, 0);
267+
getPerpendicularDirections(jetAxis, ueAxis1, ueAxis2);
268+
269+
if (ueAxis1.Mag() == 0 || ueAxis2.Mag() == 0) continue;
270+
271+
std::vector<fastjet::PseudoJet> jetConstituents = jet.constituents();
203272
for (const auto& particle : jetConstituents) {
204273

205274
auto const& track = tracks.iteratorAt(particle.user_index());
206275

207-
// Constituent Track Selection (includes DCA checks)
208276
if (!passedTrackSelection(track)) continue;
209277
if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue;
210278

211279
double pt = track.pt();
212-
// double nSigmaITSPi = static_cast<double>(itsResponse.nSigmaITS<o2::track::PID::Pion>(track));
213280

281+
if (std::abs(track.tpcNSigmaPi()) <= 3.0) {
282+
registryData.fill(HIST("pion_jet_tpc"), pt, track.tpcNSigmaPi());
283+
}
284+
if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) {
285+
registryData.fill(HIST("pion_jet_tof"), pt, track.tofNSigmaPi());
286+
}
214287

215-
// Check Pion PID (+/- 3 Sigma)
216-
double nsigmaTPCPi = track.tpcNSigmaPi();
217-
218-
// Fill TPC
219-
if (std::abs(nsigmaTPCPi) <= 3.0) {
220-
registryData.fill(HIST("pion_jet_tpc"), pt, nsigmaTPCPi);
288+
if (std::abs(track.tpcNSigmaKa()) <= 3.0) {
289+
registryData.fill(HIST("kaon_jet_tpc"), pt, track.tpcNSigmaKa());
290+
}
291+
if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) {
292+
registryData.fill(HIST("kaon_jet_tof"), pt, track.tofNSigmaKa());
221293
}
294+
295+
if (std::abs(track.tpcNSigmaPr()) <= 3.0) {
296+
registryData.fill(HIST("proton_jet_tpc"), pt, track.tpcNSigmaPr());
297+
}
298+
if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) {
299+
registryData.fill(HIST("proton_jet_tof"), pt, track.tofNSigmaPr());
300+
}
301+
}
302+
303+
for (auto const& track : tracks) {
304+
305+
if (!passedTrackSelection(track)) continue;
306+
if (std::fabs(track.dcaXY()) > maxDcaxy || std::fabs(track.dcaZ()) > maxDcaz) continue;
307+
308+
double deltaEtaUe1 = track.eta() - ueAxis1.Eta();
309+
double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi());
310+
double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1);
222311

223-
// Fill TOF
224-
if (track.hasTOF()) {
225-
double nsigmaTOFPi = track.tofNSigmaPi();
226-
if (std::abs(nsigmaTOFPi) <= 3.0) {
227-
registryData.fill(HIST("pion_jet_tof"), pt, nsigmaTOFPi);
228-
}
312+
double deltaEtaUe2 = track.eta() - ueAxis2.Eta();
313+
double deltaPhiUe2 = getDeltaPhi(track.phi(), ueAxis2.Phi());
314+
double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
315+
316+
double maxConeRadius = coneRadius;
317+
if (applyAreaCut) {
318+
maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet;
229319
}
230320

231-
// // Fill ITS
232-
// if (std::abs(nSigmaITSPi) <= 3.0) {
233-
// registryData.fill(HIST("pion_jet_its"), pt, nSigmaITSPi);
234-
// }
321+
if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius) continue;
322+
323+
double pt = track.pt();
324+
325+
if (std::abs(track.tpcNSigmaPi()) <= 3.0) {
326+
registryData.fill(HIST("pion_ue_tpc"), pt, track.tpcNSigmaPi());
327+
}
328+
if (track.hasTOF() && std::abs(track.tofNSigmaPi()) <= 3.0) {
329+
registryData.fill(HIST("pion_ue_tof"), pt, track.tofNSigmaPi());
330+
}
331+
332+
if (std::abs(track.tpcNSigmaKa()) <= 3.0) {
333+
registryData.fill(HIST("kaon_ue_tpc"), pt, track.tpcNSigmaKa());
334+
}
335+
if (track.hasTOF() && std::abs(track.tofNSigmaKa()) <= 3.0) {
336+
registryData.fill(HIST("kaon_ue_tof"), pt, track.tofNSigmaKa());
337+
}
338+
339+
if (std::abs(track.tpcNSigmaPr()) <= 3.0) {
340+
registryData.fill(HIST("proton_ue_tpc"), pt, track.tpcNSigmaPr());
341+
}
342+
if (track.hasTOF() && std::abs(track.tofNSigmaPr()) <= 3.0) {
343+
registryData.fill(HIST("proton_ue_tof"), pt, track.tofNSigmaPr());
344+
}
235345
}
236346
}
237347
}

0 commit comments

Comments
 (0)