[PWGCF] Addition of nSigma PID function

PWGCF/Flow/Tasks/CMakeLists.txt

@@ -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)
 

PWGCF/Flow/Tasks/FlowPbPbpikp.cxx → PWGCF/Flow/Tasks/flowPbpbPikp.cxx

@@ -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>
@@ -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")
@@ -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);
@@ -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}});
@@ -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");
@@ -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;
@@ -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();
@@ -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)};
 }