Addition New PID function

Author: Preet Pati

PWGCF/Flow/Tasks/flowPbpbPikp.cxx

@@ -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"
@@ -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")
@@ -97,18 +99,21 @@ 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"};
@@ -123,9 +128,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);
@@ -238,11 +245,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");
@@ -259,6 +270,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");
+
     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}});
@@ -386,13 +409,39 @@ 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>
@@ -402,6 +451,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)
@@ -411,6 +461,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)
@@ -420,6 +471,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)
@@ -437,20 +489,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;
@@ -461,6 +513,52 @@ struct FlowPbpbPikp {
     return pid + 1; // shift the pid by 1, 1 = pion, 2 = kaon, 3 = proton
   }
 
+  template <typename TTrack>
+  int getNsigmaPIDTpcTof_assymmetric(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)
   {
@@ -814,14 +912,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 ? getNsigmaPIDTpcTof_assymmetric(track) : getNsigmaPIDTpcTof(track);
 
       weff = 1; // Initializing weff for each track
       // NUA weights