Skip to content

Commit 05d68bc

Browse files
Preet-BhanjanPreet Patialibuild
authored
[PWGCF] Addition New PID function (#11582)
Co-authored-by: Preet Pati <preet@preet-6.local> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 46cc605 commit 05d68bc

File tree

1 file changed

+113
-14
lines changed

1 file changed

+113
-14
lines changed

‎PWGCF/Flow/Tasks/flowPbpbPikp.cxx‎

Lines changed: 113 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -36,6 +36,7 @@
3636
#include "Common/DataModel/TrackSelectionTables.h"
3737
#include "Common/DataModel/Centrality.h"
3838
#include "Common/DataModel/PIDResponse.h"
39+
#include "Common/DataModel/PIDResponseITS.h"
3940
#include "Common/Core/trackUtilities.h"
4041
#include "Common/DataModel/Multiplicity.h"
4142
#include "CommonConstants/PhysicsConstants.h"
@@ -82,13 +83,14 @@ struct FlowPbpbPikp {
8283
O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 3.0f, "Maximal pT for ref tracks")
8384
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
8485
O2_DEFINE_CONFIGURABLE(cfgTpcCluster, int, 70, "Number of TPC clusters")
86+
O2_DEFINE_CONFIGURABLE(cfgTpcCrossRows, int, 70, "Number of TPC clusters")
8587
O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5, "Chi2 per TPC clusters")
8688
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
8789
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, true, "Fill and output NUA weights")
8890
O2_DEFINE_CONFIGURABLE(cfgOutputRunByRun, bool, true, "Fill and output NUA weights run by run")
8991
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
9092
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
91-
O2_DEFINE_CONFIGURABLE(cfgTpcNsigmaCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
93+
O2_DEFINE_CONFIGURABLE(cfgTpcCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
9294
O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.5f, "Minimum pt to use TOF N-sigma")
9395
O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 2.0f, "DCAxy range for tracks")
9496
O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "DCAz range for tracks")
@@ -97,18 +99,20 @@ struct FlowPbpbPikp {
9799
O2_DEFINE_CONFIGURABLE(cfgCutOccupancyMax, int, 2000, "Maximum occupancy cut")
98100
O2_DEFINE_CONFIGURABLE(cfgUseGlobalTrack, bool, true, "use Global track")
99101
O2_DEFINE_CONFIGURABLE(cfgITScluster, int, 5, "Number of ITS cluster")
100-
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, true, "Use track density efficiency correction")
102+
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
101103

102104
O2_DEFINE_CONFIGURABLE(cfgUseWeightPhiEtaVtxz, bool, false, "Use Phi, Eta, VertexZ dependent NUA weights")
103105
O2_DEFINE_CONFIGURABLE(cfgUseWeightPhiPtCent, bool, false, "Use Phi, Pt, Centrality dependent NUA weights")
104106
O2_DEFINE_CONFIGURABLE(cfgUseWeightPhiEtaPt, bool, true, "Use Phi, Eta, Pt dependent NUA weights")
105107
O2_DEFINE_CONFIGURABLE(cfgUseStrictPID, bool, true, "Use strict PID cuts for TPC")
106108
O2_DEFINE_CONFIGURABLE(cfgV0AT0Acut, int, 5, "V0AT0A cut")
109+
O2_DEFINE_CONFIGURABLE(cfgUseAsymmetricPID, bool, false, "Use asymmetric PID cuts");
107110

108111
Configurable<std::vector<double>> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector<double>{0.7217476707, 0.7384792571, 0.7542625668, 0.7640680200, 0.7701951667, 0.7755299053, 0.7805901710, 0.7849446786, 0.7957356586, 0.8113039262, 0.8211968966, 0.8280558878, 0.8329342135}, "parameter 0 for track density efficiency correction"};
109112
Configurable<std::vector<double>> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector<double>{-2.169488e-05, -2.191913e-05, -2.295484e-05, -2.556538e-05, -2.754463e-05, -2.816832e-05, -2.846502e-05, -2.843857e-05, -2.705974e-05, -2.477018e-05, -2.321730e-05, -2.203315e-05, -2.109474e-05}, "parameter 1 for track density efficiency correction"};
110-
Configurable<std::vector<double>> cfgTofNsigmaCut{"cfgTofNsigmaCut", std::vector<double>{1.5, 1.5, 1.5}, "TOF n-sigma cut for pions, kaons, protons"};
111-
Configurable<std::vector<double>> cfgItsNsigmaCut{"cfgItsNsigmaCut", std::vector<double>{3, 2.5, 2}, "ITS n-sigma cut for pions, kaons, protons"};
113+
Configurable<std::vector<double>> cfgTofNsigmaCut{"cfgTofNsigmaCut", std::vector<double>{1.5, 1.5, 1.5, -1.5, -1.5, -1.5}, "TOF n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"};
114+
Configurable<std::vector<double>> cfgItsNsigmaCut{"cfgItsNsigmaCut", std::vector<double>{3, 2.5, 2, 3, 2.5, 2}, "ITS n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"};
115+
Configurable<std::vector<double>> cfgTpcNsigmaCut{"cfgTpcNsigmaCut", std::vector<double>{1.5, 1.5, 1.5, -1.5, -1.5, -1.5}, "TOF n-sigma cut for pions_posNsigma, kaons_posNsigma, protons_posNsigma, pions_negNsigma, kaons_negNsigma, protons_negNsigma"};
112116
Configurable<std::vector<int>> cfgUseEventCuts{"cfgUseEventCuts", std::vector<int>{1, 1, 0, 1, 0, 1, 1, 1, 1, 1, 1, 1, 0}, "Switch for various event cuts [kNoTimeFrameBorder, kNoITSROFrameBorder, kNoSameBunchPileup, kIsGoodZvtxFT0vsPV, kNoCollInTimeRangeStandard, kIsGoodITSLayersAll, kNoCollInRofStandard, kNoHighMultCollInPrevRof, Occupancy, Multiplicity correlation, T0AV0A 3 sigma cut, kIsVertexITSTPC, kTVXinTRD]"};
113117

114118
Configurable<GFWRegions> cfgRegions{"cfgRegions", {{"refN08", "refP08", "full", "poiN", "olN", "poiP", "olP", "poi", "ol", "poiNpi", "olNpi", "poiPpi", "olPpi", "poifullpi", "olfullpi", "poiNka", "olNka", "poiPka", "olPka", "poifullka", "olfullka", "poiNpr", "olNpr", "poiPpr", "olPpr", "poifullpr", "olfullpr"}, {-0.8, 0.4, -0.8, -0.8, -0.8, 0.4, 0.4, -0.8, -0.8, -0.8, -0.8, 0.4, 0.4, -0.8, -0.8, -0.8, -0.8, 0.4, 0.4, -0.8, -0.8, -0.8, -0.8, 0.4, 0.4, -0.8, -0.8}, {-0.4, 0.8, 0.8, -0.4, -0.4, 0.8, 0.8, 0.8, 0.8, -0.4, -0.4, 0.8, 0.8, 0.8, 0.8, -0.4, -0.4, 0.8, 0.8, 0.8, 0.8, -0.4, -0.4, 0.8, 0.8, 0.8, 0.8}, {0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}, {1, 1, 1, 128, 256, 128, 256, 128, 256, 2, 16, 2, 16, 2, 16, 4, 32, 4, 32, 4, 32, 8, 64, 8, 64, 8, 64}}, "Configurations for GFW regions"};
@@ -123,9 +127,11 @@ struct FlowPbpbPikp {
123127
ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"};
124128
ConfigurableAxis axisParticles{"axisParticles", {3, 0, 3}, "axis for different hadrons"};
125129
ConfigurableAxis axisTPCsignal{"axisTPCsignal", {10000, 0, 1000}, "axis for TPC signal"};
130+
ConfigurableAxis axisTOFbeta{"axisTOFbeta", {200, 0, 2}, "axis for TOF beta"};
126131

127132
std::vector<double> tofNsigmaCut = cfgTofNsigmaCut;
128133
std::vector<double> itsNsigmaCut = cfgItsNsigmaCut;
134+
std::vector<double> tpcNsigmaCut = cfgTofNsigmaCut;
129135

130136
Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
131137
Filter trackFilter = (nabs(aod::track::dcaXY) < cfgCutDCAxy) && (nabs(aod::track::dcaZ) < cfgCutDCAz) && (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtPOIMin) && (aod::track::pt < cfgCutPtPOIMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
@@ -238,11 +244,15 @@ struct FlowPbpbPikp {
238244
histos.add("c24_full_pi", "", {HistType::kTProfile, {axisMultiplicity}});
239245
histos.add("c24_full_ka", "", {HistType::kTProfile, {axisMultiplicity}});
240246
histos.add("c24_full_pr", "", {HistType::kTProfile, {axisMultiplicity}});
247+
248+
histos.add("TpcdEdx", "", {HistType::kTH2D, {axisPt, axisTPCsignal}});
249+
histos.add("TofBeta", "", {HistType::kTH2D, {axisPt, axisTOFbeta}});
250+
241251
histos.add("TofTpcNsigma_before", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
242252
histos.add("TofTpcNsigma_after", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
243253
histos.add("partCount", "", {HistType::kTHnSparseD, {{axisParticles, axisMultiplicity, axisPt}}});
244254

245-
histos.add("hEventCount", "Number of Event;; Count", {HistType::kTH1D, {{15, 0, 15}}});
255+
histos.add("hEventCount", "Number of Events;; Count", {HistType::kTH1D, {{15, 0, 15}}});
246256
histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(1, "Filtered event");
247257
histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(2, "After sel8");
248258
histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(3, "kNoTimeFrameBorder");
@@ -259,6 +269,18 @@ struct FlowPbpbPikp {
259269
histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(14, "kIsVertexITSTPC");
260270
histos.get<TH1>(HIST("hEventCount"))->GetXaxis()->SetBinLabel(15, "kTVXinTRD");
261271

272+
histos.add("hTrackCount", "Number of Tracks;; Count", {HistType::kTH1D, {{10, 0, 10}}});
273+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "Filtered track");
274+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "Global tracks");
275+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(3, "PV contributor");
276+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(4, "ITS clusters");
277+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(5, "TPC signal");
278+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(6, "TPC clusters");
279+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(7, "TPC crossed rows");
280+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(8, "Pions");
281+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(9, "Kaons");
282+
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(10, "Protons");
283+
262284
if (cfgOutputNUAWeights && !cfgOutputRunByRun) {
263285
histos.add<TH3>("NUA/hPhiEtaVtxz_ref", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, axisEta, axisVertex}});
264286
histos.add<TH3>("NUA/hPhiEtaVtxz_ch", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, axisEta, axisVertex}});
@@ -386,13 +408,38 @@ struct FlowPbpbPikp {
386408
template <typename TTrack>
387409
bool selectionTrack(const TTrack& track)
388410
{
389-
if (cfgUseGlobalTrack && !(track.isGlobalTrack() && track.isPVContributor() && track.itsNCls() > cfgITScluster && track.tpcNClsFound() > cfgTpcCluster && track.hasTPC())) {
390-
return false;
411+
histos.fill(HIST("hTrackCount"), 0.5); // Filtered tracks
412+
if (cfgUseGlobalTrack && !(track.isGlobalTrack())) {
413+
return 0;
391414
}
392-
if (!cfgUseGlobalTrack && !(track.isPVContributor() && track.itsNCls() > cfgITScluster && track.hasTPC())) {
393-
return false;
415+
if (cfgUseGlobalTrack)
416+
histos.fill(HIST("hTrackCount"), 1.5); // After global track selection
417+
418+
if (!(track.isPVContributor())) {
419+
return 0;
394420
}
395-
return true;
421+
histos.fill(HIST("hTrackCount"), 2.5); // After PV contributor selection
422+
423+
if (!(track.itsNCls() > cfgITScluster)) {
424+
return 0;
425+
}
426+
histos.fill(HIST("hTrackCount"), 3.5); // After ITS cluster selection
427+
428+
if (!(track.hasTPC())) {
429+
return 0;
430+
}
431+
histos.fill(HIST("hTrackCount"), 4.5); // If track has TPC signal
432+
433+
if (!(track.tpcNClsFound() > cfgTpcCluster)) {
434+
return 0;
435+
}
436+
histos.fill(HIST("hTrackCount"), 5.5); // After TPC cluster selection
437+
438+
if (!(track.tpcNClsCrossedRows() > cfgTpcCrossRows)) {
439+
return 0;
440+
}
441+
histos.fill(HIST("hTrackCount"), 6.5); // After TPC crossed rows selection
442+
return 1;
396443
}
397444

398445
template <typename TCollision, typename TTrack>
@@ -402,6 +449,7 @@ struct FlowPbpbPikp {
402449
switch (pidIndex) {
403450
case 1:
404451
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt());
452+
histos.fill(HIST("hTrackCount"), 7.5); // Pion count
405453
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
406454
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_pi_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // pion weights
407455
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
@@ -411,6 +459,7 @@ struct FlowPbpbPikp {
411459
break;
412460
case 2:
413461
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt());
462+
histos.fill(HIST("hTrackCount"), 8.5); // Kaon count
414463
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
415464
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_ka_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // kaon weights
416465
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
@@ -420,6 +469,7 @@ struct FlowPbpbPikp {
420469
break;
421470
case 3:
422471
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt());
472+
histos.fill(HIST("hTrackCount"), 9.5); // Proton count
423473
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
424474
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_pr_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // proton weights
425475
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
@@ -437,20 +487,20 @@ struct FlowPbpbPikp {
437487
std::array<float, 3> nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
438488
std::array<float, 3> nSigmaCombined = {std::hypot(track.tpcNSigmaPi(), track.tofNSigmaPi()), std::hypot(track.tpcNSigmaKa(), track.tofNSigmaKa()), std::hypot(track.tpcNSigmaPr(), track.tofNSigmaPr())};
439489
int pid = -1;
440-
float nsigma = cfgTpcNsigmaCut;
490+
float nsigma = cfgTpcCut;
441491

442492
// Choose which nSigma to use
443493
std::array<float, 3> nSigmaToUse = (track.pt() > cfgTofPtCut && track.hasTOF()) ? nSigmaCombined : nSigmaTPC;
444494
if (track.pt() > cfgTofPtCut && !track.hasTOF())
445-
return -1;
495+
return 0;
446496

447497
const int numSpecies = 3;
448498
int pidCount = 0;
449499
// Select particle with the lowest nsigma
450500
for (int i = 0; i < numSpecies; ++i) {
451501
if (std::abs(nSigmaToUse[i]) < nsigma) {
452502
if (pidCount > 0 && cfgUseStrictPID)
453-
return -1; // more than one particle with low nsigma
503+
return 0; // more than one particle with low nsigma
454504

455505
pidCount++;
456506
pid = i;
@@ -461,6 +511,52 @@ struct FlowPbpbPikp {
461511
return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
462512
}
463513

514+
template <typename TTrack>
515+
int getNsigmaPIDTpcTofAssymmetric(TTrack track)
516+
{
517+
// Computing Nsigma arrays for pion, kaon, and protons
518+
std::array<float, 3> nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
519+
std::array<float, 3> nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()};
520+
int pid = -1;
521+
522+
bool isPion, isKaon, isProton;
523+
bool isTpcPion = nSigmaTPC[0] < tpcNsigmaCut[0] && nSigmaTPC[0] > tpcNsigmaCut[0 + 3];
524+
bool isTpcKaon = nSigmaTPC[1] < tpcNsigmaCut[1] && nSigmaTPC[1] > tpcNsigmaCut[1 + 3];
525+
bool isTpcProton = nSigmaTPC[2] < tpcNsigmaCut[2] && nSigmaTPC[2] > tpcNsigmaCut[2 + 3];
526+
527+
bool isTofPion = nSigmaTOF[0] < tofNsigmaCut[0] && nSigmaTOF[0] > tofNsigmaCut[0 + 3];
528+
bool isTofKaon = nSigmaTOF[1] < tofNsigmaCut[1] && nSigmaTOF[1] > tofNsigmaCut[1 + 3];
529+
bool isTofProton = nSigmaTOF[2] < tofNsigmaCut[2] && nSigmaTOF[2] > tofNsigmaCut[2 + 3];
530+
531+
if (track.pt() > cfgTofPtCut && !track.hasTOF()) {
532+
return 0;
533+
} else if (track.pt() > cfgTofPtCut && track.hasTOF()) {
534+
isPion = isTofPion && isTpcPion;
535+
isKaon = isTofKaon && isTpcKaon;
536+
isProton = isTofProton && isTpcProton;
537+
} else {
538+
isPion = isTpcPion;
539+
isKaon = isTpcKaon;
540+
isProton = isTpcProton;
541+
}
542+
543+
if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) {
544+
return 0; // more than one particle satisfy the criteria
545+
}
546+
547+
if (isPion) {
548+
pid = PIONS;
549+
} else if (isKaon) {
550+
pid = KAONS;
551+
} else if (isProton) {
552+
pid = PROTONS;
553+
} else {
554+
return 0; // no particle satisfies the criteria
555+
}
556+
557+
return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
558+
}
559+
464560
template <char... chars>
465561
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
466562
{
@@ -814,14 +910,17 @@ struct FlowPbpbPikp {
814910
histos.fill(HIST("hEta"), track.eta());
815911
histos.fill(HIST("hPt"), pt);
816912

913+
histos.fill(HIST("TpcdEdx"), pt, track.tpcSignal());
914+
histos.fill(HIST("TofBeta"), pt, track.beta());
915+
817916
histos.fill(HIST("TofTpcNsigma_before"), PIONS, track.tpcNSigmaPi(), track.tofNSigmaPi(), pt);
818917
histos.fill(HIST("TofTpcNsigma_before"), KAONS, track.tpcNSigmaKa(), track.tofNSigmaKa(), pt);
819918
histos.fill(HIST("TofTpcNsigma_before"), PROTONS, track.tpcNSigmaPr(), track.tofNSigmaPr(), pt);
820919

821920
bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range
822921
bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); // within RF pT range
823922

824-
pidIndex = getNsigmaPIDTpcTof(track);
923+
pidIndex = cfgUseAsymmetricPID ? getNsigmaPIDTpcTofAssymmetric(track) : getNsigmaPIDTpcTof(track);
825924

826925
weff = 1; // Initializing weff for each track
827926
// NUA weights

0 commit comments

Comments
 (0)