Skip to content

Commit 371c30d

Browse files
MyFavoriteGitHubMyFavoriteGitHub
committed
Changing code in jetHadronsPid.cxx according with comments on pull request
1 parent 8b4ef7b commit 371c30d

1 file changed

Lines changed: 110 additions & 110 deletions

File tree

‎PWGJE/Tasks/jetHadronsPid.cxx‎

Lines changed: 110 additions & 110 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,7 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12-
/// \file jetHadronsPid.cxx
12+
/// \file jetHadronsPida.cxx
1313
/// \brief Analysis of hadrons in jets
1414
/// \author Aleksandra Mulewicz, WUT Warsaw, aleksandra.mulewicz@cern.ch
1515
/// \author Leonard Lorenc, WUT Warsaw, leonard.lorenc@cern.ch
@@ -21,6 +21,7 @@
2121
#include "PWGJE/Core/JetUtilities.h"
2222
#include "PWGJE/DataModel/Jet.h"
2323

24+
#include "Common/Core/RecoDecay.h"
2425
#include "Common/Core/TrackSelection.h"
2526
#include "Common/DataModel/EventSelection.h"
2627
#include "Common/DataModel/PIDResponseITS.h"
@@ -70,7 +71,6 @@ struct jetHadronsPid {
7071
HistogramRegistry registryData{"registryData", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
7172

7273
Configurable<bool> isppRefAnalysis{"isppRefAnalysis", false, "Is ppRef analysis"};
73-
Configurable<double> cfgAreaFrac{"cfgAreaFrac", 0.6, "fraction of jet area"};
7474
Configurable<double> cfgEtaJetMax{"cfgEtaJetMax", 0.5, "max jet eta"};
7575
Configurable<double> cfgMinPtTrack{"cfgMinPtTrack", 0.1, "minimum pt of tracks for jet reconstruction"};
7676

@@ -86,7 +86,7 @@ struct jetHadronsPid {
8686
Configurable<double> rJet{"rJet", 0.4, "Jet resolution parameter R"};
8787
Configurable<double> zVtx{"zVtx", 10.0, "Maximum zVertex"};
8888
Configurable<bool> applyAreaCut{"applyAreaCut", true, "apply area cut"};
89-
Configurable<double> maxNormalizedJetArea{"maxNormalizedJetArea", 1.0, "area cut"};
89+
Configurable<double> minNormalizedJetArea{"minNormalizedJetArea", 0.6, "Minimum normalized area cut to reject fake jets"};
9090
Configurable<double> deltaEtaEdge{"deltaEtaEdge", 0.05, "eta gap from the edge"};
9191

9292
Configurable<bool> requirePvContributor{"requirePvContributor", false, "require that the track is a PV contributor"};
@@ -101,6 +101,7 @@ struct jetHadronsPid {
101101
Configurable<double> maxEta{"maxEta", +0.8, "maximum eta"};
102102
Configurable<double> maxDcaxy{"maxDcaxy", 0.2, "Maximum DCAxy"};
103103
Configurable<double> maxDcaz{"maxDcaz", 0.1, "Maximum DCAz"};
104+
Configurable<double> maxNSigmaPid{"maxNSigmaPid", 3.0, "Maximum nSigma for TPC and TOF PID"};
104105

105106
Configurable<bool> setMCDefaultItsParams{"setMCDefaultItsParams", true, "set MC default parameters"};
106107

@@ -211,67 +212,64 @@ struct jetHadronsPid {
211212
registryData.add("mc_sec_proton_pt", "Reconstructed Proton Secondaries", HistType::kTH1F, {{120, 0.0, 4.0, "#it{p}_{T} (GeV/#it{c})"}});
212213
}
213214

215+
// void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2)
216+
// {
217+
// double px = p.X(), py = p.Y(), pz = p.Z();
218+
// double px2 = px * px, py2 = py * py, pz2 = pz * pz;
219+
// double pz4 = pz2 * pz2;
220+
221+
// if (px == 0 && py == 0) {
222+
// u1.SetXYZ(0, 0, 0);
223+
// u2.SetXYZ(0, 0, 0);
224+
// return;
225+
// }
226+
// if (px == 0 && py != 0) {
227+
// double ux = std::sqrt(py2 - pz4 / py2);
228+
// double uy = -pz2 / py;
229+
// u1.SetXYZ(ux, uy, pz);
230+
// u2.SetXYZ(-ux, uy, pz);
231+
// return;
232+
// }
233+
// if (py == 0 && px != 0) {
234+
// double ux = -pz2 / px;
235+
// double uy = std::sqrt(px2 - pz4 / px2);
236+
// u1.SetXYZ(ux, uy, pz);
237+
// u2.SetXYZ(ux, -uy, pz);
238+
// return;
239+
// }
240+
241+
// double a = px2 + py2;
242+
// double b = 2.0 * px * pz2;
243+
// double c = pz4 - py2 * py2 - px2 * py2;
244+
// double delta = b * b - 4.0 * a * c;
245+
246+
// if (delta < 0 || a == 0) {
247+
// u1.SetXYZ(0, 0, 0);
248+
// u2.SetXYZ(0, 0, 0);
249+
// return;
250+
// }
251+
// double u1x = (-b + std::sqrt(delta)) / (2.0 * a);
252+
// u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz);
253+
// double u2x = (-b - std::sqrt(delta)) / (2.0 * a);
254+
// u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz);
255+
// }
256+
214257
void getPerpendicularDirections(const TVector3& p, TVector3& u1, TVector3& u2)
215258
{
216-
double px = p.X(), py = p.Y(), pz = p.Z();
217-
double px2 = px * px, py2 = py * py, pz2 = pz * pz;
218-
double pz4 = pz2 * pz2;
219-
220-
if (px == 0 && py == 0) {
221-
u1.SetXYZ(0, 0, 0);
222-
u2.SetXYZ(0, 0, 0);
223-
return;
224-
}
225-
if (px == 0 && py != 0) {
226-
double ux = std::sqrt(py2 - pz4 / py2);
227-
double uy = -pz2 / py;
228-
u1.SetXYZ(ux, uy, pz);
229-
u2.SetXYZ(-ux, uy, pz);
230-
return;
231-
}
232-
if (py == 0 && px != 0) {
233-
double ux = -pz2 / px;
234-
double uy = std::sqrt(px2 - pz4 / px2);
235-
u1.SetXYZ(ux, uy, pz);
236-
u2.SetXYZ(ux, -uy, pz);
237-
return;
238-
}
239-
240-
double a = px2 + py2;
241-
double b = 2.0 * px * pz2;
242-
double c = pz4 - py2 * py2 - px2 * py2;
243-
double delta = b * b - 4.0 * a * c;
244-
245-
if (delta < 0 || a == 0) {
259+
if (p.Mag2() < 1e-9) {
246260
u1.SetXYZ(0, 0, 0);
247261
u2.SetXYZ(0, 0, 0);
248262
return;
249263
}
250-
double u1x = (-b + std::sqrt(delta)) / (2.0 * a);
251-
u1.SetXYZ(u1x, (-pz2 - px * u1x) / py, pz);
252-
double u2x = (-b - std::sqrt(delta)) / (2.0 * a);
253-
u2.SetXYZ(u2x, (-pz2 - px * u2x) / py, pz);
254-
}
255-
256-
double getDeltaPhi(double a1, double a2)
257-
{
258-
double deltaPhi(0);
259-
double phi1 = TVector2::Phi_0_2pi(a1);
260-
double phi2 = TVector2::Phi_0_2pi(a2);
261-
double diff = std::abs(phi1 - phi2);
262-
263-
if (diff <= PI)
264-
deltaPhi = diff;
265-
if (diff > PI)
266-
deltaPhi = TwoPI - diff;
267-
return deltaPhi;
264+
u1 = p.Orthogonal();
265+
u2 = p.Cross(u1);
268266
}
269267

270268
template <typename TrackIts>
271-
bool hasITSHit(const TrackIts& track, int layer)
269+
bool hasITSLayerHit(const TrackIts& track, int layer)
272270
{
273271
int ibit = layer - 1;
274-
return (track.itsClusterMap() & (1 << ibit));
272+
return (track.itsClusterMap() & (1 << ibit)) != 0;
275273
}
276274

277275
template <typename JetTrack>
@@ -287,7 +285,7 @@ struct jetHadronsPid {
287285

288286
if (!track.hasITS() || !track.hasTPC())
289287
return false;
290-
if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3)))
288+
if ((!hasITSLayerHit(track, 1)) && (!hasITSLayerHit(track, 2)) && (!hasITSLayerHit(track, 3)))
291289
return false;
292290
if (track.tpcNClsCrossedRows() < MinTpcCr)
293291
return false;
@@ -313,7 +311,7 @@ struct jetHadronsPid {
313311
return false;
314312
if (!track.hasITS() || !track.hasTPC())
315313
return false;
316-
if ((!hasITSHit(track, 1)) && (!hasITSHit(track, 2)) && (!hasITSHit(track, 3)))
314+
if ((!hasITSLayerHit(track, 1)) && (!hasITSLayerHit(track, 2)) && (!hasITSLayerHit(track, 3)))
317315
return false;
318316
if (track.itsNCls() < minItsNclusters)
319317
return false;
@@ -363,13 +361,13 @@ struct jetHadronsPid {
363361
double dcaxy = track.dcaXY();
364362
double dcaz = track.dcaZ();
365363

366-
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0);
367-
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0);
368-
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0);
364+
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid);
365+
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid);
366+
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid);
369367

370-
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0);
371-
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0);
372-
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0);
368+
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid);
369+
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid);
370+
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid);
373371

374372
const double ptThreshold = 0.8;
375373

@@ -419,15 +417,15 @@ struct jetHadronsPid {
419417
}
420418
}
421419

422-
int id(-1);
420+
// int id(-1);
423421
std::vector<fastjet::PseudoJet> fjParticles;
424422
for (auto const& track : tracks) {
425-
id++;
423+
// id++;
426424
if (!passedTrackSelectionForJetReconstruction(track))
427425
continue;
428426

429427
fastjet::PseudoJet fourMomentum(track.px(), track.py(), track.pz(), track.energy(MassPionCharged));
430-
fourMomentum.set_user_index(id);
428+
fourMomentum.set_user_index(track.index());
431429
fjParticles.emplace_back(fourMomentum);
432430
}
433431

@@ -469,10 +467,9 @@ struct jetHadronsPid {
469467
continue;
470468

471469
double normalizedJetArea = jet.area() / (PI * rJet * rJet);
472-
if (applyAreaCut && (!isppRefAnalysis) && normalizedJetArea > maxNormalizedJetArea)
473-
continue;
474-
if (isppRefAnalysis && (jet.area() < cfgAreaFrac * PI * rJet * rJet))
470+
if (applyAreaCut && normalizedJetArea < minNormalizedJetArea) {
475471
continue;
472+
}
476473

477474
double coneRadius = std::sqrt(jet.area() / PI);
478475
TVector3 jetAxis(jet.px(), jet.py(), jet.pz());
@@ -501,13 +498,13 @@ struct jetHadronsPid {
501498
double dcaxy = track.dcaXY();
502499
double dcaz = track.dcaZ();
503500

504-
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0);
505-
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0);
506-
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0);
501+
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid);
502+
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid);
503+
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid);
507504

508-
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0);
509-
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0);
510-
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0);
505+
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid);
506+
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid);
507+
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid);
511508

512509
const double ptThreshold = 0.8;
513510

@@ -565,17 +562,14 @@ struct jetHadronsPid {
565562
continue;
566563

567564
double deltaEtaUe1 = track.eta() - ueAxis1.Eta();
568-
double deltaPhiUe1 = getDeltaPhi(track.phi(), ueAxis1.Phi());
565+
double deltaPhiUe1 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis1.Phi()));
569566
double deltaRUe1 = std::sqrt(deltaEtaUe1 * deltaEtaUe1 + deltaPhiUe1 * deltaPhiUe1);
570567

571568
double deltaEtaUe2 = track.eta() - ueAxis2.Eta();
572-
double deltaPhiUe2 = getDeltaPhi(track.phi(), ueAxis2.Phi());
569+
double deltaPhiUe2 = std::abs(RecoDecay::constrainAngle(track.phi() - ueAxis2.Phi()));
573570
double deltaRUe2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
574571

575-
double maxConeRadius = coneRadius;
576-
if (applyAreaCut) {
577-
maxConeRadius = std::sqrt(maxNormalizedJetArea) * rJet;
578-
}
572+
double maxConeRadius = rJet;
579573

580574
if (deltaRUe1 > maxConeRadius && deltaRUe2 > maxConeRadius)
581575
continue;
@@ -585,13 +579,13 @@ struct jetHadronsPid {
585579
double dcaxy = track.dcaXY();
586580
double dcaz = track.dcaZ();
587581

588-
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0);
589-
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0);
590-
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0);
582+
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid);
583+
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid);
584+
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid);
591585

592-
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0);
593-
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0);
594-
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0);
586+
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid);
587+
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid);
588+
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid);
595589

596590
const double ptThreshold = 0.8;
597591

@@ -651,7 +645,7 @@ struct jetHadronsPid {
651645

652646
PROCESS_SWITCH(jetHadronsPid, processForJets, "Pid Analysis in and outside jets", true);
653647

654-
void processMC(SelectedCollisions::iterator const& collision, aod::McParticles const& mcParticles, HadronTracksMC const& tracks)
648+
void processMC(SelectedCollisions::iterator const& collision, HadronTracksMC const& tracks)
655649
{
656650
if (!collision.sel8() || std::abs(collision.posZ()) > zVtx)
657651
return;
@@ -668,24 +662,6 @@ struct jetHadronsPid {
668662
if (requireIsVertexTOFmatched && !collision.selection_bit(o2::aod::evsel::kIsVertexTOFmatched))
669663
return;
670664

671-
for (auto const& mcpart : mcParticles) {
672-
if (!mcpart.isPhysicalPrimary())
673-
continue;
674-
if (std::abs(mcpart.eta()) > 0.8)
675-
continue;
676-
677-
int pdg = std::abs(mcpart.pdgCode());
678-
double pt = mcpart.pt();
679-
680-
if (pdg == 211) {
681-
registryData.fill(HIST("mc_gen_pion_pt"), pt);
682-
} else if (pdg == 321) {
683-
registryData.fill(HIST("mc_gen_kaon_pt"), pt);
684-
} else if (pdg == 2212) {
685-
registryData.fill(HIST("mc_gen_proton_pt"), pt);
686-
}
687-
}
688-
689665
const double ptThreshold = 0.8;
690666

691667
for (auto const& track : tracks) {
@@ -696,14 +672,14 @@ struct jetHadronsPid {
696672

697673
double pt = track.pt();
698674

699-
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= 3.0);
700-
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= 3.0);
675+
bool passTpcPi = (std::abs(track.tpcNSigmaPi()) <= maxNSigmaPid);
676+
bool passTofPi = track.hasTOF() && (std::abs(track.tofNSigmaPi()) <= maxNSigmaPid);
701677

702-
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= 3.0);
703-
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= 3.0);
678+
bool passTpcKa = (std::abs(track.tpcNSigmaKa()) <= maxNSigmaPid);
679+
bool passTofKa = track.hasTOF() && (std::abs(track.tofNSigmaKa()) <= maxNSigmaPid);
704680

705-
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= 3.0);
706-
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= 3.0);
681+
bool passTpcPr = (std::abs(track.tpcNSigmaPr()) <= maxNSigmaPid);
682+
bool passTofPr = track.hasTOF() && (std::abs(track.tofNSigmaPr()) <= maxNSigmaPid);
707683

708684
if (!track.has_mcParticle())
709685
continue;
@@ -747,6 +723,30 @@ struct jetHadronsPid {
747723
}
748724
}
749725
PROCESS_SWITCH(jetHadronsPid, processMC, "Run on Monte Carlo", false);
726+
void processMCTruth(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& mcParticles)
727+
{
728+
if (std::abs(mcCollision.posZ()) > zVtx)
729+
return;
730+
731+
for (auto const& mcpart : mcParticles) {
732+
if (!mcpart.isPhysicalPrimary())
733+
continue;
734+
if (std::abs(mcpart.eta()) > maxEta)
735+
continue;
736+
737+
int pdg = std::abs(mcpart.pdgCode());
738+
double pt = mcpart.pt();
739+
740+
if (pdg == 211) {
741+
registryData.fill(HIST("mc_gen_pion_pt"), pt);
742+
} else if (pdg == 321) {
743+
registryData.fill(HIST("mc_gen_kaon_pt"), pt);
744+
} else if (pdg == 2212) {
745+
registryData.fill(HIST("mc_gen_proton_pt"), pt);
746+
}
747+
}
748+
}
749+
PROCESS_SWITCH(jetHadronsPid, processMCTruth, "Run on Monte Carlo (Pure Truth)", false);
750750
};
751751

752752
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)