[PWGHF] Add pT histogram in MC gen for efficiency calculation

PWGHF/HFC/TableProducer/correlatorHfeHadrons.cxx

@@ -36,6 +36,8 @@
 #include <Framework/InitContext.h>
 #include <Framework/runDataProcessing.h>
 
+#include <TRandom3.h>
+
 #include <cstdint>
 #include <vector>
 
@@ -48,7 +50,7 @@ using namespace o2::aod::hf_sel_electron;
 std::vector<double> zBins{VARIABLE_WIDTH, -10.0, -2.5, 2.5, 10.0};
 std::vector<double> multBins{VARIABLE_WIDTH, 0., 200., 500.0, 5000.};
 std::vector<double> multBinsMcGen{VARIABLE_WIDTH, 0., 20., 50.0, 500.}; // In MCGen multiplicity is defined by counting primaries
-using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFV0M<aod::mult::MultFV0A, aod::mult::MultFV0C>>;
+using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFT0M<aod::mult::MultFT0A, aod::mult::MultFT0C>>;
 BinningType corrBinning{{zBins, multBins}, true};
 using BinningTypeMcGen = ColumnBinningPolicy<aod::mccollision::PosZ, o2::aod::mult::MultMCFT0A>;
 
@@ -69,10 +71,20 @@ struct HfCorrelatorHfeHadrons {
   Configurable<float> etaTrackMin{"etaTrackMin", -0.8f, "Eta range  for associated hadron tracks"};
   Configurable<float> dcaXYTrackMax{"dcaXYTrackMax", 0.5f, "DCA XY cut"};
   Configurable<float> dcaZTrackMax{"dcaZTrackMax", 1.0f, "DCA Z cut"};
+  Configurable<bool> requireEmcal{"requireEmcal", true, "Require electron to be in EMCal"};
+
+  // Sigma cut for non-EMCal electrons
+  Configurable<float> tofNSigmaEl{"tofNSigmaEl", 3.0, "Sigma cut for electrons not in EMCal"};
+  Configurable<float> tpcNsigmaElectronMin{"tpcNsigmaElectronMin", -0.5f, "min Electron TPCnsigma"};
+  Configurable<float> tpcNsigmaElectronMax{"tpcNsigmaElectronMax", 3.0f, "max Electron TPCnsigma"};
 
   // Electron hadron correlation condition
   Configurable<bool> ptCondition{"ptCondition", true, "Electron pT should be greater than associate particle pT"};
 
+  Configurable<float> eventFractionToAnalyze{"eventFractionToAnalyze", -1, "Fraction of events to analyze (use only for ME offline on very large samples)"};
+
+  TRandom3 rnd{0};
+
   SliceCache cache;
   using TableCollisions = o2::soa::Filtered<o2::soa::Join<aod::Collisions, aod::Mults, aod::EvSels>>;
   using TableCollision = TableCollisions::iterator;
@@ -87,28 +99,40 @@ struct HfCorrelatorHfeHadrons {
   Preslice<aod::Tracks> perCol = aod::track::collisionId;
   Preslice<aod::HfSelEl> perCollision = aod::hf_sel_electron::collisionId;
 
+  ConfigurableAxis binsPosZ{"binsPosZ", {100, -10., 10.}, "primary vertex z coordinate"};
   ConfigurableAxis binsDeltaEta{"binsDeltaEta", {30, -1.8, 1.8}, "#it{#Delta#eta}"};
   ConfigurableAxis binsDeltaPhi{"binsDeltaPhi", {32, -o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf}, "#it{#Delta#varphi}"};
   ConfigurableAxis binsPt{"binsPt", {50, 0.0, 50}, "#it{p_{T}}(GeV/#it{c})"};
   ConfigurableAxis binsPoolBin{"binsPoolBin", {9, 0., 9.}, "PoolBin"};
+  ConfigurableAxis binsNSigma{"binsNSigma", {30, -15., 15.}, "#it{#sigma_{TPC}}"};
 
   HistogramRegistry registry{
     "registry",
     {}};
 
   void init(InitContext&)
   {
+    AxisSpec axisPosZ = {binsPosZ, "Pos Z"};
     AxisSpec axisDeltaEta = {binsDeltaEta, "#Delta #eta = #eta_{Electron}- #eta_{Hadron}"};
     AxisSpec axisDeltaPhi = {binsDeltaPhi, "#Delta #varphi = #varphi_{Electron}- #varphi_{Hadron}"};
     AxisSpec axisPt = {binsPt, "#it{p_{T}}(GeV/#it{c})"};
     AxisSpec axisPoolBin = {binsPoolBin, "PoolBin"};
+    AxisSpec axisNSigma = {binsNSigma, "it{#sigma_{TPC}}"};
 
+    registry.add("hZvertex", "z vertex", {HistType::kTH1D, {axisPosZ}});
+    registry.add("hNevents", "No of events", {HistType::kTH1D, {{3, 1, 4}}});
     registry.add("hInclusiveEHCorrel", "Sparse for Delta phi and Delta eta Inclusive Electron with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
     registry.add("hLSEHCorrel", "Sparse for Delta phi and Delta eta Like sign Electron pair  with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
     registry.add("hULSEHCorrel", "Sparse for Delta phi and Delta eta  UnLike sign Electron pair with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
+    registry.add("hTofnSigmaVsP", " Tof nSigma info vs P; n#sigma;#it{p} (GeV#it{/c});passEMcal;", {HistType::kTH2F, {{axisNSigma}, {axisPt}}});
+    registry.add("hTpcnSigmaVsP", " TPC nSigma info vs P; n#sigma;#it{p} (GeV#it{/c});passEMcal;", {HistType::kTH2F, {{axisNSigma}, {axisPt}}});
+
     registry.add("hMCgenNonHfEHCorrel", "Sparse for Delta phi and Delta eta  for McGen Non Hf Electron  with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
     registry.add("hMCgenInclusiveEHCorrl", "Sparse for Delta phi and Delta eta  for McGen Electron pair  with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
     registry.add("hptElectron", "hptElectron", {HistType::kTH1D, {axisPt}});
+    registry.add("hptHadron", "hptHadron", {HistType::kTH1D, {axisPt}});
+    registry.add("hMCgenptHadron", "hMCgenptHadron", {HistType::kTH1D, {axisPt}});
+    registry.add("hMCgenptHadronprimary", "hMCgenptHadronprimary", {HistType::kTH1D, {axisPt}});
 
     registry.add("hMixEventInclusiveEHCorrl", "Sparse for mix event Delta phi and Delta eta Inclusive Electron with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
     registry.add("hMixEventLSEHCorrel", "Sparse for mix event Delta phi and Delta eta Like sign Electron pair with Hadron;p_{T}^{e} (GeV#it{/c});p_{T}^{h} (GeV#it{/c});#Delta#varphi;#Delta#eta;", {HistType::kTHnSparseF, {{axisPt}, {axisPt}, {axisDeltaPhi}, {axisDeltaEta}}});
@@ -147,18 +171,30 @@ struct HfCorrelatorHfeHadrons {
   {
     if (!(isRun3 ? collision.sel8() : (collision.sel7() && collision.alias_bit(kINT7))))
       return;
-    int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFV0M()));
+    int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
     auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
     int gCollisionId = collision.globalIndex();
     int64_t timeStamp = bc.timestamp();
 
+    bool skipEventTableFilling = false;
+    if (eventFractionToAnalyze > 0) {
+      if (rnd.Uniform(0, 1) > eventFractionToAnalyze) {
+        skipEventTableFilling = true;
+      }
+    }
+
+    registry.fill(HIST("hNevents"), 1);
     // Add hadron Table For Mix Event Electron Hadron correlation
-    for (const auto& hTrack : tracks) {
-      if (!selAssoHadron(hTrack)) {
-        continue;
+    if (!skipEventTableFilling) {
+      registry.fill(HIST("hZvertex"), collision.posZ());
+      for (const auto& hTrack : tracks) {
+        if (!selAssoHadron(hTrack)) {
+          continue;
+        }
+        registry.fill(HIST("hTracksBin"), poolBin);
+        registry.fill(HIST("hptHadron"), hTrack.pt());
+        entryHadron(hTrack.phi(), hTrack.eta(), hTrack.pt(), poolBin, gCollisionId, timeStamp);
       }
-      registry.fill(HIST("hTracksBin"), poolBin);
-      entryHadron(hTrack.phi(), hTrack.eta(), hTrack.pt(), poolBin, gCollisionId, timeStamp);
     }
 
     //  Construct Deta Phi between electrons and hadrons
@@ -171,16 +207,35 @@ struct HfCorrelatorHfeHadrons {
       ptElectron = eTrack.ptTrack();
       phiElectron = eTrack.phiTrack();
       etaElectron = eTrack.etaTrack();
+      bool acceptElectron = false;
 
       double deltaPhi = -999;
       double deltaEta = -999;
       double ptHadron = -999;
       double etaHadron = -999;
       double phiHadron = -999;
-
-      if (!eTrack.isEmcal()) {
-        continue;
+      if (eTrack.isEmcal()) {
+        // EMCal electron
+        if (requireEmcal) {
+          acceptElectron = true;
+        }
+      } else {
+        // Non-EMCal electron
+        if (!requireEmcal) {
+          // Apply sigma cut
+          if (std::abs(eTrack.tofNSigmaElTrack()) < tofNSigmaEl) {
+            if (eTrack.tpcNSigmaElTrack() > tpcNsigmaElectronMin && eTrack.tpcNSigmaElTrack() < tpcNsigmaElectronMax) {
+              registry.fill(HIST("hTofnSigmaVsP"), eTrack.tofNSigmaElTrack(), eTrack.ptTrack());
+              registry.fill(HIST("hTpcnSigmaVsP"), eTrack.tpcNSigmaElTrack(), eTrack.ptTrack());
+              acceptElectron = true;
+            }
+          }
+        }
       }
+
+      if (!acceptElectron)
+        continue; // skip electron if not passing criteria
+
       registry.fill(HIST("hptElectron"), ptElectron);
       int nElectronLS = 0;
       int nElectronUS = 0;
@@ -199,8 +254,10 @@ struct HfCorrelatorHfeHadrons {
         }
       }
 
-      registry.fill(HIST("hElectronBin"), poolBin);
-      entryElectron(phiElectron, etaElectron, ptElectron, nElectronLS, nElectronUS, poolBin, gCollisionId, timeStamp);
+      if (!skipEventTableFilling) {
+        registry.fill(HIST("hElectronBin"), poolBin);
+        entryElectron(phiElectron, etaElectron, ptElectron, nElectronLS, nElectronUS, poolBin, gCollisionId, timeStamp);
+      }
 
       for (const auto& hTrack : tracks) {
         // Apply Hadron cut
@@ -260,7 +317,7 @@ struct HfCorrelatorHfeHadrons {
     double ptHadronMix = -999;
     double etaHadronMix = -999;
     double phiHadronMix = -999;
-    int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFV0M()));
+    int poolBin = corrBinning.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()));
     for (const auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) {
       if (!t1.isEmcal()) {
         continue;
@@ -332,9 +389,25 @@ struct HfCorrelatorHfeHadrons {
     BinningTypeMcGen corrBinningMcGen{{zBins, multBinsMcGen}, true};
     int poolBin = corrBinningMcGen.getBin(std::make_tuple(mcCollision.posZ(), mcCollision.multMCFT0A()));
 
+    for (const auto& particleMc : mcParticles) {
+      if (particleMc.eta() < etaTrackMin || particleMc.eta() > etaTrackMax) {
+        continue;
+      }
+      if (particleMc.pt() < ptTrackMin) {
+        continue;
+      }
+
+      registry.fill(HIST("hMCgenptHadron"), particleMc.pt());
+      if (particleMc.isPhysicalPrimary()) {
+
+        registry.fill(HIST("hMCgenptHadronprimary"), particleMc.pt());
+      }
+    }
+
     double ptElectron = 0;
     double phiElectron = 0;
     double etaElectron = 0;
+
     for (const auto& electronMc : electron) {
       double ptHadron = 0;
       double phiHadron = 0;

PWGHF/HFL/TableProducer/electronSelectionWithTpcEmcal.cxx

@@ -118,6 +118,10 @@ struct HfElectronSelectionWithTpcEmcal {
   Configurable<float> m20EmcClusterElectronMin{"m20EmcClusterElectronMin", 0.0f, "min Electron  EMCal Cluster M20"};
   Configurable<float> tpcNsigmaElectronMin{"tpcNsigmaElectronMin", -0.5f, "min Electron TPCnsigma"};
   Configurable<float> tpcNsigmaElectronMax{"tpcNsigmaElectronMax", 3.0f, "max Electron TPCnsigma"};
+  Configurable<int> pdgCodeCharmMin{"pdgCodeCharmMin", 400, "Min Charm Hadron PdgCode"};
+  Configurable<int> pdgCodeCharmMax{"pdgCodeCharmMax", 600, "Max Charm Hadron PdgCode"};
+  Configurable<int> pdgCodeBeautyMin{"pdgCodeBeautyMin", 4000, "Min beauty Hadron PdgCode"};
+  Configurable<int> pdgCodeBeautyMax{"pdgCodeBeautyMax", 6000, "Max beauty Hadron PdgCode"};
 
   using TableCollisions = o2::soa::Filtered<o2::soa::Join<aod::Collisions, aod::Mults, aod::EvSels>>;
   using TableCollision = TableCollisions::iterator;
@@ -550,12 +554,9 @@ struct HfElectronSelectionWithTpcEmcal {
   void processMcGen(McGenTableCollision const& mcCollision, aod::McParticles const& mcParticles)
   {
 
-    ///// electron identification
     bool isNonHfe = false;
     for (const auto& particleMc : mcParticles) {
 
-      if (!particleMc.isPhysicalPrimary())
-        continue;
       if (!mcGensel(particleMc)) {
         continue;
       }
@@ -565,71 +566,101 @@ struct HfElectronSelectionWithTpcEmcal {
         bool isEmbEta = false;
         bool isEmbPi0 = false;
 
-        if (particleMc.has_mothers()) {
-          // Check first mother
-          auto const& mother = particleMc.mothers_first_as<aod::McParticles>();
+        // Check first mother
+        auto const& mother = particleMc.mothers_first_as<aod::McParticles>();
 
-          if (std::abs(mother.pdgCode()) == kEta || std::abs(mother.pdgCode()) == kPi0 || std::abs(mother.pdgCode()) == kGamma) {
-            registry.fill(HIST("hMcgenAllNonHfeElectron"), particleMc.pt());
-            if (mother.has_mothers()) {
-              auto const& gmother = mother.mothers_first_as<aod::McParticles>();
-              if (gmother.has_mothers()) {
-                auto const& ggmother = gmother.mothers_first_as<aod::McParticles>();
+        if (std::abs(mother.pdgCode()) == kEta || std::abs(mother.pdgCode()) == kPi0 || std::abs(mother.pdgCode()) == kGamma) {
+          registry.fill(HIST("hMcgenAllNonHfeElectron"), particleMc.pt());
+          auto const& gmother = mother.mothers_first_as<aod::McParticles>();
+          auto const& ggmother = gmother.mothers_first_as<aod::McParticles>();
+          auto const& gggmother = ggmother.mothers_first_as<aod::McParticles>();
 
-                // cases to consider: eta->e, eta->pi0->e, eta->gamma->e, eta->pi0->gamma->e, pi0->e, pi0->gamma->e
+          // cases to consider: eta->e, eta->pi0->e, eta->gamma->e, eta->pi0->gamma->e, pi0->e, pi0->gamma->e
 
-                //=================  eta->e ======================================
-                if (std::abs(mother.pdgCode()) == kEta) {
-                  isEmbEta = true;
-                }
-                //=================  eta->pi0->e ======================================
+          //=================  eta->e ======================================
+          if (std::abs(mother.pdgCode()) == kEta) {
+            if (mother.isPhysicalPrimary()) {
+              if ((std::abs(gmother.pdgCode()) >= pdgCodeCharmMin && std::abs(gmother.pdgCode()) < pdgCodeCharmMax) ||
+                  (std::abs(gmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(gmother.pdgCode()) < pdgCodeBeautyMax)) {
+                continue;
+              }
+              isEmbEta = true;
+            }
+          }
+          //=================  eta->pi0->e ======================================
 
-                if (std::abs(mother.pdgCode()) == kPi0) {
-                  isEmbPi0 = true; // pi0 -> e
+          if (std::abs(mother.pdgCode()) == kPi0) {
+            if (mother.isPhysicalPrimary()) {
+              if ((std::abs(gmother.pdgCode()) >= pdgCodeCharmMin && std::abs(gmother.pdgCode()) < pdgCodeCharmMax) ||
+                  (std::abs(gmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(gmother.pdgCode()) < pdgCodeBeautyMax)) {
+                continue;
+              }
+              isEmbPi0 = true; // pi0 -> e
+            }
 
-                  if (std::abs(gmother.pdgCode()) == kEta) {
-                    isEmbEta = true; // eta->pi0-> e
-                  }
+            if (std::abs(gmother.pdgCode()) == kEta) {
+              if (gmother.isPhysicalPrimary()) {
+                if ((std::abs(ggmother.pdgCode()) >= pdgCodeCharmMin && std::abs(ggmother.pdgCode()) < pdgCodeCharmMax) ||
+                    (std::abs(ggmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(ggmother.pdgCode()) < pdgCodeBeautyMax)) {
+                  continue;
                 }
+                isEmbEta = true; // eta->pi0-> e
+              }
+            }
+          }
 
-                /// ====================================  eta->gamma->e  and eta->pi0->gamma->e============
-                if (std::abs(mother.pdgCode()) == kGamma) {
-                  if (std::abs(gmother.pdgCode()) == kEta) {
-                    isEmbEta = true; // eta->gamma-> e
-                  }
-
-                  if (std::abs(gmother.pdgCode()) == kPi0) {
-                    isEmbPi0 = true; // pi0-> gamma-> e
+          /// ====================================  eta->gamma->e  and eta->pi0->gamma->e============
+          if (std::abs(mother.pdgCode()) == kGamma) {
+            if (std::abs(gmother.pdgCode()) == kEta) {
+              if (gmother.isPhysicalPrimary()) {
+                if ((std::abs(ggmother.pdgCode()) >= pdgCodeCharmMin && std::abs(ggmother.pdgCode()) < pdgCodeCharmMax) ||
+                    (std::abs(ggmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(ggmother.pdgCode()) < pdgCodeBeautyMax)) {
+                  continue;
+                }
+                isEmbEta = true; // eta->gamma-> e
+              }
+            }
 
-                    if (std::abs(ggmother.pdgCode()) == kEta) {
+            if (std::abs(gmother.pdgCode()) == kPi0) {
+              if (gmother.isPhysicalPrimary()) {
+                if ((std::abs(ggmother.pdgCode()) >= pdgCodeCharmMin && std::abs(ggmother.pdgCode()) < pdgCodeCharmMax) ||
+                    (std::abs(ggmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(ggmother.pdgCode()) < pdgCodeBeautyMax)) {
+                  continue;
+                }
+                isEmbPi0 = true; // pi0-> gamma-> e
+              }
 
-                      isEmbEta = true; // eta->pi0->gamma-> e
-                    }
+              if (std::abs(ggmother.pdgCode()) == kEta) {
+                if (ggmother.isPhysicalPrimary()) {
+                  if ((std::abs(gggmother.pdgCode()) >= pdgCodeCharmMin && std::abs(gggmother.pdgCode()) < pdgCodeCharmMax) ||
+                      (std::abs(gggmother.pdgCode()) >= pdgCodeBeautyMin && std::abs(gggmother.pdgCode()) < pdgCodeBeautyMax)) {
+                    continue;
                   }
+                  isEmbEta = true; // eta->pi0->gamma-> e
                 }
               }
             }
           }
-        }
-        if (isEmbPi0 || isEmbEta) {
-          registry.fill(HIST("hMcgenNonHfeElectron"), particleMc.pt());
-          isNonHfe = true;
-          if (isEmbPi0) {
+          if (isEmbPi0 || isEmbEta) {
+            registry.fill(HIST("hMcgenNonHfeElectron"), particleMc.pt());
+            isNonHfe = true;
+            if (isEmbPi0) {
 
-            registry.fill(HIST("hPi0eEmbTrkPt"), particleMc.pt());
-          }
-          if (isEmbEta) {
-            registry.fill(HIST("hEtaeEmbTrkPt"), particleMc.pt());
+              registry.fill(HIST("hPi0eEmbTrkPt"), particleMc.pt());
+            }
+            if (isEmbEta) {
+              registry.fill(HIST("hEtaeEmbTrkPt"), particleMc.pt());
+            }
           }
         }
+
         hfGenElectronSel(mcCollision.globalIndex(), particleMc.globalIndex(), particleMc.eta(), particleMc.phi(), particleMc.pt(), isNonHfe);
       }
     }
   }
 
   PROCESS_SWITCH(HfElectronSelectionWithTpcEmcal, processMcGen, "Process MC Gen mode", false);
 };
-
 WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
 {
   return WorkflowSpec{adaptAnalysisTask<HfElectronSelectionWithTpcEmcal>(cfgc)};