Skip to content

Commit

Permalink
[PWGCF] Addition of nSigma PID function (#8838)
Browse files Browse the repository at this point in the history 
Co-authored-by: Preet Pati <preet@preet-2.local>
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
  • Loading branch information
@Preet-Bhanjan @alibuild
3 people authored Dec 9, 2024
1 parent 1a597d4 commit c9e57b7
Show file tree
Hide file tree
Showing 2 changed files with 94 additions and 59 deletions.
4 changes: 2 additions & 2 deletions PWGCF/Flow/Tasks/CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -39,8 +39,8 @@ o2physics_add_dpl_workflow(flow-gf
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::GFWCore
COMPONENT_NAME Analysis)

o2physics_add_dpl_workflow(flow-pbpb-pikp-task
SOURCES FlowPbPbpikp.cxx
o2physics_add_dpl_workflow(flow-pbpb-pikp
SOURCES flowPbpbPikp.cxx
PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::GFWCore
COMPONENT_NAME Analysis)

Expand Down
149 changes: 92 additions & 57 deletions PWGCF/Flow/Tasks/FlowPbPbpikp.cxx → PWGCF/Flow/Tasks/flowPbpbPikp.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,10 +9,13 @@
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.

/// \file flowPbpbPikp.cxx
/// \brief PID flow using the generic framework
/// \author Preet Bhanjan Pati <bhanjanpreet@gmail.com>

#include <CCDB/BasicCCDBManager.h>
#include <cmath>
#include <vector>
#include <iostream>
#include <utility>
#include <array>
#include <string>
Expand Down Expand Up @@ -48,13 +51,14 @@
using namespace o2;
using namespace o2::framework;
using namespace o2::framework::expressions;
using namespace std;

#define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable<TYPE> NAME{#NAME, DEFAULT, HELP};

struct GfwPidflow {
struct FlowPbpbPikp {
Service<ccdb::BasicCCDBManager> ccdb;
Configurable<int64_t> nolaterthan{"ccdb-no-later-than", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
Configurable<std::string> url{"ccdb-url", "http://ccdb-test.cern.ch:8080", "url of the ccdb repository"};
Configurable<int64_t> noLaterThan{"noLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://ccdb-test.cern.ch:8080", "url of the ccdb repository"};

O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
Expand All @@ -69,11 +73,11 @@ struct GfwPidflow {
ConfigurableAxis axisVertex{"axisVertex", {20, -10, 10}, "vertex axis for histograms"};
ConfigurableAxis axisPhi{"axisPhi", {60, 0.0, constants::math::TwoPI}, "phi axis for histograms"};
ConfigurableAxis axisEta{"axisEta", {40, -1., 1.}, "eta axis for histograms"};
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.2, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.20, 1.40, 1.60, 1.80, 2.00, 2.20, 2.40, 2.60, 2.80, 3.00, 3.50, 4.00, 5.00, 6.00, 8.00, 10.00}, "pt axis for histograms"};
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.20, 1.40, 1.60, 1.80, 2.00, 2.20, 2.40, 2.60, 2.80, 3.00, 3.50, 4.00, 5.00, 6.00, 8.00, 10.00}, "pt axis for histograms"};
ConfigurableAxis axisMultiplicity{"axisMultiplicity", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "centrality axis for histograms"};
ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {80, -5, 5}, "nsigmaTPC axis"};
ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {80, -5, 5}, "nsigmaTOF axis"};
ConfigurableAxis axisparticles{"axisparticles", {3, 0, 3}, "axis for different hadrons"};
ConfigurableAxis axisParticles{"axisParticles", {3, 0, 3}, "axis for different hadrons"};

Filter collisionFilter = nabs(aod::collision::posZ) < cfgCutVertex;
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);
Expand All @@ -86,14 +90,15 @@ struct GfwPidflow {
TAxis* fPtAxis;
TRandom3* fRndm = new TRandom3(0);

using aodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>>;
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>>;
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::FT0Mults, aod::MultZeqs, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs>>;
// 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>>;
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::pidTPCFullPi, aod::pidTPCFullKa, aod::pidTPCFullPr, aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullPr>>;

void init(InitContext const&)
{
ccdb->setURL(url.value);
ccdb->setURL(ccdbUrl.value);
ccdb->setCaching(true);
ccdb->setCreatedNotAfter(nolaterthan.value);
ccdb->setCreatedNotAfter(noLaterThan.value);

histos.add("hPhi", "", {HistType::kTH1D, {axisPhi}});
histos.add("hEta", "", {HistType::kTH1D, {axisEta}});
Expand All @@ -106,28 +111,28 @@ struct GfwPidflow {
histos.add("c22_gap08_ka", "", {HistType::kTProfile, {axisMultiplicity}});
histos.add("c22_gap08_pr", "", {HistType::kTProfile, {axisMultiplicity}});
histos.add("c24_full", "", {HistType::kTProfile, {axisMultiplicity}});
histos.add("TofTpcNsigma", "", {HistType::kTHnSparseD, {{axisparticles, axisNsigmaTPC, axisNsigmaTOF}}});

histos.add("TofTpcNsigma", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
histos.add("partCount", "", {HistType::kTHnSparseD, {{axisParticles, axisMultiplicity, axisPt}}});
o2::framework::AxisSpec axis = axisPt;
int nPtBins = axis.binEdges.size() - 1;
double* PtBins = &(axis.binEdges)[0];
fPtAxis = new TAxis(nPtBins, PtBins);
double* ptBins = &(axis.binEdges)[0];
fPtAxis = new TAxis(nPtBins, ptBins);

TObjArray* oba = new TObjArray();
oba->Add(new TNamed("Ch08Gap22", "Ch08Gap22"));
for (Int_t i = 0; i < fPtAxis->GetNbins(); i++)
for (int i = 0; i < fPtAxis->GetNbins(); i++)
oba->Add(new TNamed(Form("Ch08Gap22_pt_%i", i + 1), "Ch08Gap22_pTDiff"));
oba->Add(new TNamed("Pi08Gap22", "Pi08Gap22"));
for (Int_t i = 0; i < fPtAxis->GetNbins(); i++)
for (int i = 0; i < fPtAxis->GetNbins(); i++)
oba->Add(new TNamed(Form("Pi08Gap22_pt_%i", i + 1), "Pi08Gap22_pTDiff"));
oba->Add(new TNamed("Ka08Gap22", "Ka08Gap22"));
for (Int_t i = 0; i < fPtAxis->GetNbins(); i++)
for (int i = 0; i < fPtAxis->GetNbins(); i++)
oba->Add(new TNamed(Form("Ka08Gap22_pt_%i", i + 1), "Ka08Gap22_pTDiff"));
oba->Add(new TNamed("Pr08Gap22", "Pr08Gap22"));
for (Int_t i = 0; i < fPtAxis->GetNbins(); i++)
for (int i = 0; i < fPtAxis->GetNbins(); i++)
oba->Add(new TNamed(Form("Pr08Gap22_pt_%i", i + 1), "Pr08Gap22_pTDiff"));
oba->Add(new TNamed("ChFull24", "ChFull24"));
for (Int_t i = 0; i < fPtAxis->GetNbins(); i++)
for (int i = 0; i < fPtAxis->GetNbins(); i++)
oba->Add(new TNamed(Form("ChFull24_pt_%i", i + 1), "ChFull24_pTDiff"));

fFC->SetName("FlowContainer");
Expand Down Expand Up @@ -169,8 +174,36 @@ struct GfwPidflow {
fGFW->CreateRegions();
}

enum Particles {
PIONS,
KAONS,
PROTONS
};

template <typename TTrack>
std::pair<int, int> GetBayesID(TTrack track)
int getNsigmaPID(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> nSigmaCombined = {std::hypot(track.tpcNSigmaPi(), track.tofNSigmaPi()), std::hypot(track.tpcNSigmaKa(), track.tofNSigmaKa()), std::hypot(track.tpcNSigmaPr(), track.tofNSigmaPr())};
int pid = -1;
float nsigma = 3.0;

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

// Select particle with the lowest nsigma
for (int i = 0; i < 3; ++i) {
if (std::abs(nSigmaToUse[i]) < nsigma) {
pid = i;
nsigma = std::abs(nSigmaToUse[i]);
}
}
return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
}

/*template <typename TTrack>
std::pair<int, int> getBayesID(TTrack track)
{
std::array<int, 3> bayesprobs = {static_cast<int>(track.bayesPi()), static_cast<int>(track.bayesKa()), static_cast<int>(track.bayesPr())};
int bayesid = -1;
Expand All @@ -186,50 +219,50 @@ struct GfwPidflow {
}
template <typename TTrack>
int GetBayesPIDIndex(TTrack track)
int getBayesPIDIndex(TTrack track)
{
int maxProb[3] = {80, 80, 80};
int pidID = -1;
std::pair<int, int> idprob = GetBayesID(track);
if (idprob.first == 0 || idprob.first == 1 || idprob.first == 2) { // 0 = pion, 1 = kaon, 2 = proton
std::pair<int, int> idprob = getBayesID(track);
if (idprob.first == PIONS || idprob.first == KAONS || idprob.first == PROTONS) { // 0 = pion, 1 = kaon, 2 = proton
pidID = idprob.first;
float nsigmaTPC[3] = {track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr()};
if (idprob.second > maxProb[pidID]) {
if (abs(nsigmaTPC[pidID]) > 3)
if (std::fabs(nsigmaTPC[pidID]) > 3)
return 0;
return pidID + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
} else {
return 0;
}
}
return 0;
}
}*/

template <char... chars>
void FillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
{
double dnx, val;
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
if (dnx == 0)
return;
if (!corrconf.pTDif) {
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
if (TMath::Abs(val) < 1)
if (std::fabs(val) < 1)
histos.fill(tarName, cent, val, dnx);
return;
}
for (Int_t i = 1; i <= fPtAxis->GetNbins(); i++) {
for (int i = 1; i <= fPtAxis->GetNbins(); i++) {
dnx = fGFW->Calculate(corrconf, i - 1, kTRUE).real();
if (dnx == 0)
continue;
val = fGFW->Calculate(corrconf, i - 1, kFALSE).real() / dnx;
if (TMath::Abs(val) < 1)
if (std::fabs(val) < 1)
histos.fill(tarName, fPtAxis->GetBinCenter(i), val, dnx);
}
return;
}

void FillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
void fillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
{
double dnx, val;
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
Expand All @@ -238,84 +271,86 @@ struct GfwPidflow {
}
if (!corrconf.pTDif) {
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
if (TMath::Abs(val) < 1) {
if (std::fabs(val) < 1) {
fFC->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
}
return;
}
for (Int_t i = 1; i <= fPtAxis->GetNbins(); i++) {
for (int i = 1; i <= fPtAxis->GetNbins(); i++) {
dnx = fGFW->Calculate(corrconf, i - 1, kTRUE).real();
if (dnx == 0)
continue;
val = fGFW->Calculate(corrconf, i - 1, kFALSE).real() / dnx;
if (TMath::Abs(val) < 1)
if (std::fabs(val) < 1)
fFC->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
}
return;
}

void process(aodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, aodTracks const& tracks)
void process(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracks const& tracks)
{
int Ntot = tracks.size();
if (Ntot < 1)
int nTot = tracks.size();
if (nTot < 1)
return;
if (!collision.sel8())
return;
float l_Random = fRndm->Rndm();
float lRandom = fRndm->Rndm();

float vtxz = collision.posZ();
histos.fill(HIST("hVtxZ"), vtxz);
histos.fill(HIST("hMult"), Ntot);
histos.fill(HIST("hMult"), nTot);
histos.fill(HIST("hCent"), collision.centFT0C());
fGFW->Clear();
const auto cent = collision.centFT0C();
float weff = 1, wacc = 1;
int pidIndex;
for (auto& track : tracks) {
for (auto const& track : tracks) {
double pt = track.pt();
histos.fill(HIST("hPhi"), track.phi());
histos.fill(HIST("hEta"), track.eta());
histos.fill(HIST("hPt"), pt);

histos.fill(HIST("TofTpcNsigma"), 0, track.tpcNSigmaPi(), track.tofNSigmaPi());
histos.fill(HIST("TofTpcNsigma"), 1, track.tpcNSigmaKa(), track.tofNSigmaKa());
histos.fill(HIST("TofTpcNsigma"), 2, track.tpcNSigmaPr(), track.tofNSigmaPr());
histos.fill(HIST("TofTpcNsigma"), PIONS, track.tpcNSigmaPi(), track.tofNSigmaPi(), pt);
histos.fill(HIST("TofTpcNsigma"), KAONS, track.tpcNSigmaKa(), track.tofNSigmaKa(), pt);
histos.fill(HIST("TofTpcNsigma"), 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
bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range
bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); // within RF pT range

pidIndex = GetBayesPIDIndex(track);
if (WithinPtRef)
// pidIndex = getBayesPIDIndex(track);
pidIndex = getNsigmaPID(track);
if (withinPtRef)
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1);
if (WithinPtPOI)
if (withinPtPOI)
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 128);
if (WithinPtPOI && WithinPtRef)
if (withinPtPOI && withinPtRef)
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 256);
fGFW->Fill(track.eta(), 1, track.phi(), wacc * weff, 512);

if (pidIndex) {
if (WithinPtPOI)
histos.fill(HIST("partCount"), pidIndex - 1, cent, pt);
if (withinPtPOI)
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex));
if (WithinPtPOI && WithinPtRef)
if (withinPtPOI && withinPtRef)
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex + 3));
}
}

// Filling c22 with ROOT TProfile
FillProfile(corrconfigs.at(0), HIST("c22_gap08"), cent);
FillProfile(corrconfigs.at(1), HIST("c22_gap08_pi"), cent);
FillProfile(corrconfigs.at(2), HIST("c22_gap08_ka"), cent);
FillProfile(corrconfigs.at(3), HIST("c22_gap08_pr"), cent);
FillProfile(corrconfigs.at(4), HIST("c24_full"), cent);
fillProfile(corrconfigs.at(0), HIST("c22_gap08"), cent);
fillProfile(corrconfigs.at(1), HIST("c22_gap08_pi"), cent);
fillProfile(corrconfigs.at(2), HIST("c22_gap08_ka"), cent);
fillProfile(corrconfigs.at(3), HIST("c22_gap08_pr"), cent);
fillProfile(corrconfigs.at(4), HIST("c24_full"), cent);

for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {
FillFC(corrconfigs.at(l_ind), cent, l_Random);
fillFC(corrconfigs.at(l_ind), cent, lRandom);
}

} // end of process
};

WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
{
return WorkflowSpec{adaptAnalysisTask<GfwPidflow>(cfgc)};
return WorkflowSpec{adaptAnalysisTask<FlowPbpbPikp>(cfgc)};
}

0 comments on commit c9e57b7

Please sign in to comment.