Skip to content

[PWGCF] Addition New PID function #11582

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
victor-gonzalez merged 5 commits into from
Jun 12, 2025
Merged

[PWGCF] Addition New PID function #11582

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
64efc3b
Addition New PID function
Jun 12, 2025
4c43add
Recommiting after precommit fixes
Jun 12, 2025
3ba06f0
Recommiting after O2linter fixes
Jun 12, 2025
af3eaed
Please consider the following formatting changes
alibuild Jun 12, 2025
989463a
Merge pull request #17 from alibuild/alibot-cleanup-11582
Preet-Bhanjan Jun 12, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
127 changes: 113 additions & 14 deletions PWGCF/Flow/Tasks/flowPbpbPikp.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@
#include "Common/DataModel/TrackSelectionTables.h"
#include "Common/DataModel/Centrality.h"
#include "Common/DataModel/PIDResponse.h"
#include "Common/DataModel/PIDResponseITS.h"
#include "Common/Core/trackUtilities.h"
#include "Common/DataModel/Multiplicity.h"
#include "CommonConstants/PhysicsConstants.h"
Expand Down Expand Up @@ -82,13 +83,14 @@ struct FlowPbpbPikp {
O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 3.0f, "Maximal pT for ref tracks")
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
O2_DEFINE_CONFIGURABLE(cfgTpcCluster, int, 70, "Number of TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgTpcCrossRows, int, 70, "Number of TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5, "Chi2 per TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, true, "Fill and output NUA weights")
O2_DEFINE_CONFIGURABLE(cfgOutputRunByRun, bool, true, "Fill and output NUA weights run by run")
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
O2_DEFINE_CONFIGURABLE(cfgTpcNsigmaCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
O2_DEFINE_CONFIGURABLE(cfgTpcCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.5f, "Minimum pt to use TOF N-sigma")
O2_DEFINE_CONFIGURABLE(cfgCutDCAxy, float, 2.0f, "DCAxy range for tracks")
O2_DEFINE_CONFIGURABLE(cfgCutDCAz, float, 2.0f, "DCAz range for tracks")
Expand All @@ -97,18 +99,20 @@ struct FlowPbpbPikp {
O2_DEFINE_CONFIGURABLE(cfgCutOccupancyMax, int, 2000, "Maximum occupancy cut")
O2_DEFINE_CONFIGURABLE(cfgUseGlobalTrack, bool, true, "use Global track")
O2_DEFINE_CONFIGURABLE(cfgITScluster, int, 5, "Number of ITS cluster")
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, true, "Use track density efficiency correction")
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")

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

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"};
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"};
Configurable<std::vector<double>> cfgTofNsigmaCut{"cfgTofNsigmaCut", std::vector<double>{1.5, 1.5, 1.5}, "TOF n-sigma cut for pions, kaons, protons"};
Configurable<std::vector<double>> cfgItsNsigmaCut{"cfgItsNsigmaCut", std::vector<double>{3, 2.5, 2}, "ITS n-sigma cut for pions, kaons, protons"};
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"};
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"};
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"};
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]"};

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"};
Expand All @@ -123,9 +127,11 @@ struct FlowPbpbPikp {
ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"};
ConfigurableAxis axisParticles{"axisParticles", {3, 0, 3}, "axis for different hadrons"};
ConfigurableAxis axisTPCsignal{"axisTPCsignal", {10000, 0, 1000}, "axis for TPC signal"};
ConfigurableAxis axisTOFbeta{"axisTOFbeta", {200, 0, 2}, "axis for TOF beta"};

std::vector<double> tofNsigmaCut = cfgTofNsigmaCut;
std::vector<double> itsNsigmaCut = cfgItsNsigmaCut;
std::vector<double> tpcNsigmaCut = cfgTofNsigmaCut;

Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
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);
Expand Down Expand Up @@ -238,11 +244,15 @@ struct FlowPbpbPikp {
histos.add("c24_full_pi", "", {HistType::kTProfile, {axisMultiplicity}});
histos.add("c24_full_ka", "", {HistType::kTProfile, {axisMultiplicity}});
histos.add("c24_full_pr", "", {HistType::kTProfile, {axisMultiplicity}});

histos.add("TpcdEdx", "", {HistType::kTH2D, {axisPt, axisTPCsignal}});
histos.add("TofBeta", "", {HistType::kTH2D, {axisPt, axisTOFbeta}});

histos.add("TofTpcNsigma_before", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
histos.add("TofTpcNsigma_after", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
histos.add("partCount", "", {HistType::kTHnSparseD, {{axisParticles, axisMultiplicity, axisPt}}});

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

histos.add("hTrackCount", "Number of Tracks;; Count", {HistType::kTH1D, {{10, 0, 10}}});
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(1, "Filtered track");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(2, "Global tracks");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(3, "PV contributor");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(4, "ITS clusters");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(5, "TPC signal");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(6, "TPC clusters");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(7, "TPC crossed rows");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(8, "Pions");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(9, "Kaons");
histos.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(10, "Protons");
Comment on lines +273 to +282
Copy link
Collaborator

@victor-gonzalez victor-gonzalez Jun 12, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I would suggest using named constants and starting the axis range in -0.5
It will be much more readable and easier to evolve here and when the histogram is filled

Of course, the same for the event count histogram


if (cfgOutputNUAWeights && !cfgOutputRunByRun) {
histos.add<TH3>("NUA/hPhiEtaVtxz_ref", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, axisEta, axisVertex}});
histos.add<TH3>("NUA/hPhiEtaVtxz_ch", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, axisEta, axisVertex}});
Expand Down Expand Up @@ -386,13 +408,38 @@ struct FlowPbpbPikp {
template <typename TTrack>
bool selectionTrack(const TTrack& track)
{
if (cfgUseGlobalTrack && !(track.isGlobalTrack() && track.isPVContributor() && track.itsNCls() > cfgITScluster && track.tpcNClsFound() > cfgTpcCluster && track.hasTPC())) {
return false;
histos.fill(HIST("hTrackCount"), 0.5); // Filtered tracks
if (cfgUseGlobalTrack && !(track.isGlobalTrack())) {
return 0;
}
if (!cfgUseGlobalTrack && !(track.isPVContributor() && track.itsNCls() > cfgITScluster && track.hasTPC())) {
return false;
if (cfgUseGlobalTrack)
histos.fill(HIST("hTrackCount"), 1.5); // After global track selection

if (!(track.isPVContributor())) {
return 0;
}
return true;
histos.fill(HIST("hTrackCount"), 2.5); // After PV contributor selection

if (!(track.itsNCls() > cfgITScluster)) {
return 0;
}
histos.fill(HIST("hTrackCount"), 3.5); // After ITS cluster selection

if (!(track.hasTPC())) {
return 0;
}
histos.fill(HIST("hTrackCount"), 4.5); // If track has TPC signal

if (!(track.tpcNClsFound() > cfgTpcCluster)) {
return 0;
}
histos.fill(HIST("hTrackCount"), 5.5); // After TPC cluster selection

if (!(track.tpcNClsCrossedRows() > cfgTpcCrossRows)) {
return 0;
}
histos.fill(HIST("hTrackCount"), 6.5); // After TPC crossed rows selection
return 1;
}

template <typename TCollision, typename TTrack>
Expand All @@ -402,6 +449,7 @@ struct FlowPbpbPikp {
switch (pidIndex) {
case 1:
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt());
histos.fill(HIST("hTrackCount"), 7.5); // Pion count
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_pi_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // pion weights
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
Expand All @@ -411,6 +459,7 @@ struct FlowPbpbPikp {
break;
case 2:
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt());
histos.fill(HIST("hTrackCount"), 8.5); // Kaon count
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_ka_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // kaon weights
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
Expand All @@ -420,6 +469,7 @@ struct FlowPbpbPikp {
break;
case 3:
histos.fill(HIST("TofTpcNsigma_after"), pidIndex - 1, track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt());
histos.fill(HIST("hTrackCount"), 9.5); // Proton count
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiEtaVtxz)
histos.fill(HIST("PhiCorrected/hPhiEtaVtxz_pr_corrd"), track.phi(), track.eta(), collision.posZ(), wacc); // proton weights
if (!cfgAcceptance.value.empty() && cfgUseWeightPhiPtCent)
Expand All @@ -437,20 +487,20 @@ struct FlowPbpbPikp {
std::array<float, 3> nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
std::array<float, 3> nSigmaCombined = {std::hypot(track.tpcNSigmaPi(), track.tofNSigmaPi()), std::hypot(track.tpcNSigmaKa(), track.tofNSigmaKa()), std::hypot(track.tpcNSigmaPr(), track.tofNSigmaPr())};
int pid = -1;
float nsigma = cfgTpcNsigmaCut;
float nsigma = cfgTpcCut;

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

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

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

template <typename TTrack>
int getNsigmaPIDTpcTofAssymmetric(TTrack track)
{
// Computing Nsigma arrays for pion, kaon, and protons
std::array<float, 3> nSigmaTPC = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
std::array<float, 3> nSigmaTOF = {track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr()};
int pid = -1;

bool isPion, isKaon, isProton;
bool isTpcPion = nSigmaTPC[0] < tpcNsigmaCut[0] && nSigmaTPC[0] > tpcNsigmaCut[0 + 3];
bool isTpcKaon = nSigmaTPC[1] < tpcNsigmaCut[1] && nSigmaTPC[1] > tpcNsigmaCut[1 + 3];
bool isTpcProton = nSigmaTPC[2] < tpcNsigmaCut[2] && nSigmaTPC[2] > tpcNsigmaCut[2 + 3];

bool isTofPion = nSigmaTOF[0] < tofNsigmaCut[0] && nSigmaTOF[0] > tofNsigmaCut[0 + 3];
bool isTofKaon = nSigmaTOF[1] < tofNsigmaCut[1] && nSigmaTOF[1] > tofNsigmaCut[1 + 3];
bool isTofProton = nSigmaTOF[2] < tofNsigmaCut[2] && nSigmaTOF[2] > tofNsigmaCut[2 + 3];

if (track.pt() > cfgTofPtCut && !track.hasTOF()) {
return 0;
} else if (track.pt() > cfgTofPtCut && track.hasTOF()) {
isPion = isTofPion && isTpcPion;
isKaon = isTofKaon && isTpcKaon;
isProton = isTofProton && isTpcProton;
} else {
isPion = isTpcPion;
isKaon = isTpcKaon;
isProton = isTpcProton;
}

if ((isPion && isKaon) || (isPion && isProton) || (isKaon && isProton)) {
return 0; // more than one particle satisfy the criteria
}

if (isPion) {
pid = PIONS;
} else if (isKaon) {
pid = KAONS;
} else if (isProton) {
pid = PROTONS;
} else {
return 0; // no particle satisfies the criteria
}

return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
}

template <char... chars>
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
{
Expand Down Expand Up @@ -814,14 +910,17 @@ struct FlowPbpbPikp {
histos.fill(HIST("hEta"), track.eta());
histos.fill(HIST("hPt"), pt);

histos.fill(HIST("TpcdEdx"), pt, track.tpcSignal());
histos.fill(HIST("TofBeta"), pt, track.beta());

histos.fill(HIST("TofTpcNsigma_before"), PIONS, track.tpcNSigmaPi(), track.tofNSigmaPi(), pt);
histos.fill(HIST("TofTpcNsigma_before"), KAONS, track.tpcNSigmaKa(), track.tofNSigmaKa(), pt);
histos.fill(HIST("TofTpcNsigma_before"), PROTONS, track.tpcNSigmaPr(), track.tofNSigmaPr(), pt);

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

pidIndex = getNsigmaPIDTpcTof(track);
pidIndex = cfgUseAsymmetricPID ? getNsigmaPIDTpcTofAssymmetric(track) : getNsigmaPIDTpcTof(track);

weff = 1; // Initializing weff for each track
// NUA weights
Expand Down
Loading