Skip to content

Commit e4d35b1

Browse files
Preet-BhanjanPreet Patialibuild
authored
[PWGCF] PID flow for reconstructed particles (#9233)
Co-authored-by: Preet Pati <preet@preet-2.local> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 309df8d commit e4d35b1

File tree

1 file changed

+115
-26
lines changed

1 file changed

+115
-26
lines changed

‎PWGCF/Flow/Tasks/flowPbpbPikp.cxx‎

Lines changed: 115 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -11,7 +11,7 @@
1111

1212
/// \file flowPbpbPikp.cxx
1313
/// \brief PID flow using the generic framework
14-
/// \author Preet Bhanjan Pati <bhanjanpreet@gmail.com>
14+
/// \author Preet Bhanjan Pati <preet.bhanjan.pati@cern.ch>
1515

1616
#include <CCDB/BasicCCDBManager.h>
1717
#include <cmath>
@@ -20,6 +20,8 @@
2020
#include <array>
2121
#include <string>
2222

23+
#include "Math/Vector4D.h"
24+
2325
#include "Framework/runDataProcessing.h"
2426
#include "Framework/AnalysisTask.h"
2527
#include "Framework/ASoAHelpers.h"
@@ -69,6 +71,8 @@ struct FlowPbpbPikp {
6971
O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5, "Chi2 per TPC clusters")
7072
O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables")
7173
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
74+
O2_DEFINE_CONFIGURABLE(cfgTpcNsigmaCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
75+
O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 1.8f, "Minimum pt to use TOF N-sigma")
7276

7377
ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"};
7478
ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
@@ -78,10 +82,20 @@ struct FlowPbpbPikp {
7882
ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"};
7983
ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"};
8084
ConfigurableAxis axisParticles{"axisParticles", {3, 0, 3}, "axis for different hadrons"};
85+
ConfigurableAxis axisPhiMass{"axisPhiMass", {10000, 0, 2}, "axis for invariant mass distibution for Phi"};
86+
ConfigurableAxis axisTPCsignal{"axisTPCsignal", {10000, 0, 1000}, "axis for TPC signal"};
8187

8288
Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
8389
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtPOIMin) && (aod::track::pt < cfgCutPtPOIMax) && ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
8490

91+
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>>;
92+
// using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::pidBayes, aod::pidBayesPi, aod::pidBayesKa, aod::pidBayesPr, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
93+
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
94+
95+
SliceCache cache;
96+
Partition<AodTracks> posTracks = aod::track::signed1Pt > 0.0f;
97+
Partition<AodTracks> negTracks = aod::track::signed1Pt < 0.0f;
98+
8599
OutputObj<FlowContainer> fFC{FlowContainer("FlowContainer")};
86100
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
87101

@@ -90,10 +104,6 @@ struct FlowPbpbPikp {
90104
TAxis* fPtAxis;
91105
TRandom3* fRndm = new TRandom3(0);
92106

93-
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>>;
94-
// using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::pidBayes, aod::pidBayesPi, aod::pidBayesKa, aod::pidBayesPr, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
95-
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;
96-
97107
void init(InitContext const&)
98108
{
99109
ccdb->setURL(ccdbUrl.value);
@@ -106,13 +116,17 @@ struct FlowPbpbPikp {
106116
histos.add("hMult", "", {HistType::kTH1D, {{3000, 0.5, 3000.5}}});
107117
histos.add("hCent", "", {HistType::kTH1D, {{90, 0, 90}}});
108118
histos.add("hPt", "", {HistType::kTH1D, {axisPt}});
119+
histos.add("hPhiMass", "", {HistType::kTH1D, {axisPhiMass}});
109120
histos.add("c22_gap08", "", {HistType::kTProfile, {axisMultiplicity}});
110121
histos.add("c22_gap08_pi", "", {HistType::kTProfile, {axisMultiplicity}});
111122
histos.add("c22_gap08_ka", "", {HistType::kTProfile, {axisMultiplicity}});
112123
histos.add("c22_gap08_pr", "", {HistType::kTProfile, {axisMultiplicity}});
113124
histos.add("c24_full", "", {HistType::kTProfile, {axisMultiplicity}});
125+
histos.add("KplusTPC", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}});
126+
histos.add("KminusTPC", "", {HistType::kTH2D, {{axisPt, axisTPCsignal}}});
114127
histos.add("TofTpcNsigma", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
115128
histos.add("partCount", "", {HistType::kTHnSparseD, {{axisParticles, axisMultiplicity, axisPt}}});
129+
116130
o2::framework::AxisSpec axis = axisPt;
117131
int nPtBins = axis.binEdges.size() - 1;
118132
double* ptBins = &(axis.binEdges)[0];
@@ -143,6 +157,8 @@ struct FlowPbpbPikp {
143157
fGFW->AddRegion("refN08", -0.8, -0.4, 1, 1);
144158
fGFW->AddRegion("refP08", 0.4, 0.8, 1, 1);
145159
fGFW->AddRegion("full", -0.8, 0.8, 1, 512);
160+
fGFW->AddRegion("poi", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 1024);
161+
fGFW->AddRegion("ol", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 2048);
146162

147163
// charged parts
148164
fGFW->AddRegion("poiN", -0.8, -0.4, 1 + fPtAxis->GetNbins(), 128);
@@ -171,6 +187,7 @@ struct FlowPbpbPikp {
171187
corrconfigs.push_back(fGFW->GetCorrelatorConfig("poiNpr refN08 | olNpr {2} refP08 {-2}", "Pr08Gap22", kTRUE));
172188

173189
corrconfigs.push_back(fGFW->GetCorrelatorConfig("full {2 2 -2 -2}", "ChFull24", kFALSE));
190+
corrconfigs.push_back(fGFW->GetCorrelatorConfig("poi full | ol {2 2 -2 -2}", "ChFull24", kTRUE));
174191
fGFW->CreateRegions();
175192
}
176193

@@ -180,17 +197,28 @@ struct FlowPbpbPikp {
180197
PROTONS
181198
};
182199

200+
template <typename TTrack>
201+
bool isFakeKaon(TTrack track)
202+
{
203+
const auto pglobal = track.p();
204+
const auto ptpc = track.tpcInnerParam();
205+
if (std::abs(pglobal - ptpc) > 0.1) {
206+
return true;
207+
}
208+
return false;
209+
}
210+
183211
template <typename TTrack>
184212
int getNsigmaPID(TTrack track)
185213
{
186214
// Computing Nsigma arrays for pion, kaon, and protons
187215
std::array<float, 3> nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
188216
std::array<float, 3> nSigmaCombined = {std::hypot(track.tpcNSigmaPi(), track.tofNSigmaPi()), std::hypot(track.tpcNSigmaKa(), track.tofNSigmaKa()), std::hypot(track.tpcNSigmaPr(), track.tofNSigmaPr())};
189217
int pid = -1;
190-
float nsigma = 3.0;
218+
float nsigma = cfgTpcNsigmaCut;
191219

192220
// Choose which nSigma to use
193-
std::array<float, 3> nSigmaToUse = (track.pt() > 0.4 && track.hasTOF()) ? nSigmaCombined : nSigmaTPC;
221+
std::array<float, 3> nSigmaToUse = (track.pt() > cfgTofPtCut && track.hasTOF()) ? nSigmaCombined : nSigmaTPC;
194222

195223
// Select particle with the lowest nsigma
196224
for (int i = 0; i < 3; ++i) {
@@ -238,6 +266,24 @@ struct FlowPbpbPikp {
238266
return 0;
239267
}*/
240268

269+
template <typename TTrack, typename vector, char... chars>
270+
void resurrectParticle(TTrack trackplus, TTrack trackminus, vector plusdaug, vector minusdaug, vector mom, double plusmass, double minusmass, const ConstStr<chars...>& hist)
271+
{
272+
for (auto const& [partplus, partminus] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(trackplus, trackminus))) {
273+
if (getNsigmaPID(partplus) != 2)
274+
continue;
275+
if (getNsigmaPID(partminus) != 2)
276+
continue;
277+
278+
plusdaug = ROOT::Math::PxPyPzMVector(partplus.px(), partplus.py(), partplus.pz(), plusmass);
279+
minusdaug = ROOT::Math::PxPyPzMVector(partminus.px(), partminus.py(), partminus.pz(), minusmass);
280+
mom = plusdaug + minusdaug;
281+
282+
histos.fill(hist, mom.M());
283+
}
284+
return;
285+
}
286+
241287
template <char... chars>
242288
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
243289
{
@@ -287,6 +333,10 @@ struct FlowPbpbPikp {
287333
return;
288334
}
289335

336+
ROOT::Math::PxPyPzMVector Phimom, kplusdaug, kminusdaug;
337+
double massKplus = o2::constants::physics::MassKPlus;
338+
double massKminus = o2::constants::physics::MassKMinus;
339+
290340
void process(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracks const& tracks)
291341
{
292342
int nTot = tracks.size();
@@ -302,39 +352,78 @@ struct FlowPbpbPikp {
302352
histos.fill(HIST("hCent"), collision.centFT0C());
303353
fGFW->Clear();
304354
const auto cent = collision.centFT0C();
355+
356+
auto posSlicedTracks = posTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
357+
auto negSlicedTracks = negTracks->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
358+
305359
float weff = 1, wacc = 1;
306360
int pidIndex;
307-
for (auto const& track : tracks) {
308-
double pt = track.pt();
309-
histos.fill(HIST("hPhi"), track.phi());
310-
histos.fill(HIST("hEta"), track.eta());
361+
362+
// resurrectParticle(posSlicedTracks, negSlicedTracks, kplusdaug, kminusdaug, Phimom, massKplus, massKminus, HIST("hPhiMass"));
363+
364+
for (auto const& trackA : posSlicedTracks) {
365+
if (getNsigmaPID(trackA) != 2)
366+
continue;
367+
if (isFakeKaon(trackA))
368+
continue;
369+
auto trackAID = trackA.globalIndex();
370+
371+
for (auto const& trackB : negSlicedTracks) {
372+
auto trackBID = trackB.globalIndex();
373+
if (getNsigmaPID(trackB) != 2)
374+
continue;
375+
if (isFakeKaon(trackB))
376+
continue;
377+
if (trackAID == trackBID)
378+
continue;
379+
380+
histos.fill(HIST("KplusTPC"), trackA.pt(), trackA.tpcSignal());
381+
histos.fill(HIST("KminusTPC"), trackB.pt(), trackB.tpcSignal());
382+
383+
kplusdaug = ROOT::Math::PxPyPzMVector(trackA.px(), trackA.py(), trackA.pz(), massKplus);
384+
kminusdaug = ROOT::Math::PxPyPzMVector(trackB.px(), trackB.py(), trackB.pz(), massKminus);
385+
Phimom = kplusdaug + kminusdaug;
386+
387+
histos.fill(HIST("hPhiMass"), Phimom.M());
388+
}
389+
}
390+
391+
for (auto const& track1 : tracks) {
392+
double pt = track1.pt();
393+
histos.fill(HIST("hPhi"), track1.phi());
394+
histos.fill(HIST("hEta"), track1.eta());
311395
histos.fill(HIST("hPt"), pt);
312396

313-
histos.fill(HIST("TofTpcNsigma"), PIONS, track.tpcNSigmaPi(), track.tofNSigmaPi(), pt);
314-
histos.fill(HIST("TofTpcNsigma"), KAONS, track.tpcNSigmaKa(), track.tofNSigmaKa(), pt);
315-
histos.fill(HIST("TofTpcNsigma"), PROTONS, track.tpcNSigmaPr(), track.tofNSigmaPr(), pt);
397+
histos.fill(HIST("TofTpcNsigma"), PIONS, track1.tpcNSigmaPi(), track1.tofNSigmaPi(), pt);
398+
histos.fill(HIST("TofTpcNsigma"), KAONS, track1.tpcNSigmaKa(), track1.tofNSigmaKa(), pt);
399+
histos.fill(HIST("TofTpcNsigma"), PROTONS, track1.tpcNSigmaPr(), track1.tofNSigmaPr(), pt);
316400

317401
bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range
318402
bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); // within RF pT range
319403

320-
// pidIndex = getBayesPIDIndex(track);
321-
pidIndex = getNsigmaPID(track);
322-
if (withinPtRef)
323-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1);
324-
if (withinPtPOI)
325-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 128);
326-
if (withinPtPOI && withinPtRef)
327-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 256);
328-
fGFW->Fill(track.eta(), 1, track.phi(), wacc * weff, 512);
404+
// pidIndex = getBayesPIDIndex(track1);
405+
pidIndex = getNsigmaPID(track1);
406+
if (withinPtRef) {
407+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 1);
408+
fGFW->Fill(track1.eta(), 1, track1.phi(), wacc * weff, 512);
409+
}
410+
if (withinPtPOI) {
411+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 128);
412+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 1024);
413+
}
414+
if (withinPtPOI && withinPtRef) {
415+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 256);
416+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 2048);
417+
}
329418

330419
if (pidIndex) {
331420
histos.fill(HIST("partCount"), pidIndex - 1, cent, pt);
332421
if (withinPtPOI)
333-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex));
422+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 1 << (pidIndex));
334423
if (withinPtPOI && withinPtRef)
335-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex + 3));
424+
fGFW->Fill(track1.eta(), fPtAxis->FindBin(pt) - 1, track1.phi(), wacc * weff, 1 << (pidIndex + 3));
336425
}
337-
}
426+
} // track1 loop ends
338427

339428
// Filling c22 with ROOT TProfile
340429
fillProfile(corrconfigs.at(0), HIST("c22_gap08"), cent);

0 commit comments

Comments
 (0)