[PWGHF] Add finding leading particle in correlatorLcHadrons (AliceO2Group#8183)

Co-authored-by: Vít Kučera <[email protected]>

Author: 1481014945

PWGHF/HFC/DataModel/CorrelationTables.h

@@ -97,21 +97,23 @@ DECLARE_SOA_TABLE(DHadronRecoInfo, "AOD", "DHADRONRECOINFO", //! D0-Hadrons pair
 // Note: definition of columns and tables for Lc-Hadron correlation pairs
 namespace hf_correlation_lc_hadron
 {
-DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float);       //! DeltaPhi between Lc and Hadrons
-DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float);       //! DeltaEta between Lc and Hadrons
-DECLARE_SOA_COLUMN(PtLc, ptLc, float);               //! Transverse momentum of Lc
-DECLARE_SOA_COLUMN(PtHadron, ptHadron, float);       //! Transverse momentum of Hadron
-DECLARE_SOA_COLUMN(MLc, mLc, float);                 //! Invariant mass of Lc
-DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int); //! Tag for LcToPKPi/LcToPiKP
-DECLARE_SOA_COLUMN(PoolBin, poolBin, int);           //! Pool Bin for the MixedEvent
+DECLARE_SOA_COLUMN(DeltaPhi, deltaPhi, float);                //! DeltaPhi between Lc and Hadrons
+DECLARE_SOA_COLUMN(DeltaEta, deltaEta, float);                //! DeltaEta between Lc and Hadrons
+DECLARE_SOA_COLUMN(PtLc, ptLc, float);                        //! Transverse momentum of Lc
+DECLARE_SOA_COLUMN(PtHadron, ptHadron, float);                //! Transverse momentum of Hadron
+DECLARE_SOA_COLUMN(MLc, mLc, float);                          //! Invariant mass of Lc
+DECLARE_SOA_COLUMN(SignalStatus, signalStatus, int);          //! Tag for LcToPKPi/LcToPiKP
+DECLARE_SOA_COLUMN(PoolBin, poolBin, int);                    //! Pool Bin for the MixedEvent
+DECLARE_SOA_COLUMN(IsAutoCorrelated, isAutoCorrelated, bool); //! Correlation Status
 } // namespace hf_correlation_lc_hadron
 
 DECLARE_SOA_TABLE(LcHadronPair, "AOD", "LCHPAIR", //! Lc-Hadrons pairs Informations
                   aod::hf_correlation_lc_hadron::DeltaPhi,
                   aod::hf_correlation_lc_hadron::DeltaEta,
                   aod::hf_correlation_lc_hadron::PtLc,
                   aod::hf_correlation_lc_hadron::PtHadron,
-                  aod::hf_correlation_lc_hadron::PoolBin);
+                  aod::hf_correlation_lc_hadron::PoolBin,
+                  aod::hf_correlation_lc_hadron::IsAutoCorrelated);
 
 DECLARE_SOA_TABLE(LcHadronRecoInfo, "AOD", "LCHRECOINFO", //! Lc-Hadrons pairs Reconstructed Informations
                   aod::hf_correlation_lc_hadron::MLc,

PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx

@@ -15,6 +15,8 @@
 /// \author Samrangy Sadhu <[email protected]>, INFN Bari
 /// \author Swapnesh Santosh Khade <[email protected]>, IIT Indore
 
+#include <vector>
+
 #include "CommonConstants/PhysicsConstants.h"
 #include "Framework/AnalysisTask.h"
 #include "Framework/HistogramRegistry.h"
@@ -30,12 +32,14 @@
 #include "PWGHF/DataModel/CandidateReconstructionTables.h"
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/HFC/DataModel/CorrelationTables.h"
+#include "PWGHF/HFC/Utils/utilsCorrelations.h"
 
 using namespace o2;
 using namespace o2::analysis;
 using namespace o2::constants::physics;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
+using namespace o2::analysis::hf_correlations;
 
 ///
 /// Returns deltaPhi value in range [-pi/2., 3.*pi/2], typically used for correlation studies
@@ -264,44 +268,6 @@ struct HfCorrelatorD0Hadrons {
     registry.add("hCountD0TriggersGen", "D0 trigger particles - MC gen;;N of trigger D0", {HistType::kTH2F, {{1, -0.5, 0.5}, {vbins, "#it{p}_{T} (GeV/#it{c})"}}});
   }
 
-  // Find Leading Particle
-  template <typename TTracks>
-  int findLeadingParticle(TTracks const& tracks)
-  {
-    auto leadingParticle = tracks.begin();
-    for (auto const& track : tracks) {
-      if (std::abs(track.dcaXY()) >= 1. || std::abs(track.dcaZ()) >= 1.) {
-        continue;
-      }
-      if (track.pt() > leadingParticle.pt()) {
-        leadingParticle = track;
-      }
-    }
-    int leadingIndex = leadingParticle.globalIndex();
-    return leadingIndex;
-  }
-  // ======= Find Leading Particle for McGen ============
-  template <typename TMcParticles>
-  int findLeadingParticleMcGen(TMcParticles const& mcParticles)
-  {
-    auto leadingParticle = mcParticles.begin();
-    for (auto const& mcParticle : mcParticles) {
-      if (std::abs(mcParticle.eta()) > etaTrackMax) {
-        continue;
-      }
-      if (mcParticle.pt() < ptTrackMin) {
-        continue;
-      }
-      if ((std::abs(mcParticle.pdgCode()) != kElectron) && (std::abs(mcParticle.pdgCode()) != kMuonMinus) && (std::abs(mcParticle.pdgCode()) != kPiPlus) && (std::abs(mcParticle.pdgCode()) != kKPlus) && (std::abs(mcParticle.pdgCode()) != kProton)) {
-        continue;
-      }
-      if (mcParticle.pt() > leadingParticle.pt()) {
-        leadingParticle = mcParticle;
-      }
-    }
-    int leadingIndex = leadingParticle.globalIndex();
-    return leadingIndex;
-  }
   // =======  Process starts for Data, Same event ============
 
   /// D0-h correlation pair builder - for real data and data-like analysis (i.e. reco-level w/o matching request via MC truth)
@@ -315,7 +281,7 @@ struct HfCorrelatorD0Hadrons {
     }
     // find leading particle
     if (correlateD0WithLeadingParticle) {
-      leadingIndex = findLeadingParticle(tracks);
+      leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value);
     }
 
     int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
@@ -459,7 +425,7 @@ struct HfCorrelatorD0Hadrons {
     }
     // find leading particle
     if (correlateD0WithLeadingParticle) {
-      leadingIndex = findLeadingParticle(tracks);
+      leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value);
     }
     int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
     int nTracks = 0;
@@ -629,7 +595,7 @@ struct HfCorrelatorD0Hadrons {
     // MC gen level
     // find leading particle
     if (correlateD0WithLeadingParticle) {
-      leadingIndex = findLeadingParticleMcGen(mcParticles);
+      leadingIndex = findLeadingParticleMcGen(mcParticles, etaTrackMax.value, ptTrackMin.value);
     }
     for (const auto& particle1 : mcParticles) {
       // check if the particle is D0 or D0bar (for general plot filling and selection, so both cases are fine) - NOTE: decay channel is not probed!

PWGHF/HFC/TableProducer/correlatorLcHadrons.cxx

@@ -15,6 +15,8 @@
 /// \author Marianna Mazzilli <[email protected]>
 /// \author Zhen Zhang <[email protected]>
 
+#include <vector>
+
 #include "CommonConstants/PhysicsConstants.h"
 #include "Framework/AnalysisTask.h"
 #include "Framework/HistogramRegistry.h"
@@ -30,12 +32,14 @@
 #include "PWGHF/DataModel/CandidateReconstructionTables.h"
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/HFC/DataModel/CorrelationTables.h"
+#include "PWGHF/HFC/Utils/utilsCorrelations.h"
 
 using namespace o2;
 using namespace o2::analysis;
 using namespace o2::constants::physics;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
+using namespace o2::analysis::hf_correlations;
 
 ///
 /// Returns deltaPhi values in range [-pi/2., 3.*pi/2.], typically used for correlation studies
@@ -189,10 +193,14 @@ struct HfCorrelatorLcHadrons {
   ConfigurableAxis binsMultiplicity{"binsMultiplicity", {VARIABLE_WIDTH, 0.0f, 2000.0f, 6000.0f, 100000.0f}, "Mixing bins - multiplicity"};
   ConfigurableAxis binsZVtx{"binsZVtx", {VARIABLE_WIDTH, -10.0f, -2.5f, 2.5f, 10.0f}, "Mixing bins - z-vertex"};
   ConfigurableAxis binsMultiplicityMc{"binsMultiplicityMc", {VARIABLE_WIDTH, 0.0f, 20.0f, 50.0f, 500.0f}, "Mixing bins - MC multiplicity"}; // In MCGen multiplicity is defined by counting tracks
+  Configurable<bool> storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"};
+  Configurable<bool> correlateLcWithLeadingParticle{"correlateLcWithLeadingParticle", false, "Switch for correlation of Lc baryons with leading particle only"};
 
   HfHelper hfHelper;
   SliceCache cache;
   BinningType corrBinning{{binsZVtx, binsMultiplicity}, true};
+  int leadingIndex = 0;
+  bool correlationStatus = false;
 
   // Filters for ME
   Filter collisionFilter = aod::hf_selection_lc_collision::lcSel >= filterFlagLc;
@@ -263,6 +271,12 @@ struct HfCorrelatorLcHadrons {
     if (selectedLcCandidates.size() == 0) {
       return;
     }
+
+    // find leading particle
+    if (correlateLcWithLeadingParticle) {
+      leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value);
+    }
+
     int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
     int nTracks = 0;
     if (collision.numContrib() > 1) {
@@ -335,24 +349,37 @@ struct HfCorrelatorLcHadrons {
         if (std::abs(track.dcaXY()) >= dcaXYTrackMax || std::abs(track.dcaZ()) >= dcaZTrackMax) {
           continue; // Remove secondary tracks
         }
+
         // Remove Lc daughters by checking track indices
         if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex()) || (candidate.prong2Id() == track.globalIndex())) {
-          continue;
+          if (!storeAutoCorrelationFlag) {
+            continue;
+          }
+          correlationStatus = true;
         }
+
+        if (correlateLcWithLeadingParticle) {
+          if (track.globalIndex() != leadingIndex) {
+            continue;
+          }
+        }
+
         if (candidate.isSelLcToPKPi() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(track.phi(), candidate.phi()),
                             track.eta() - candidate.eta(),
                             candidate.pt(),
                             track.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPKPi(candidate), false);
         }
         if (candidate.isSelLcToPiKP() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(track.phi(), candidate.phi()),
                             track.eta() - candidate.eta(),
                             candidate.pt(),
                             track.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPiKP(candidate), false);
         }
       } // Hadron Tracks loop
@@ -370,6 +397,12 @@ struct HfCorrelatorLcHadrons {
     if (selectedLcCandidatesMc.size() == 0) {
       return;
     }
+
+    // find leading particle
+    if (correlateLcWithLeadingParticle) {
+      leadingIndex = findLeadingParticle(tracks, dcaXYTrackMax.value, dcaZTrackMax.value);
+    }
+
     int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
     int nTracks = 0;
     if (collision.numContrib() > 1) {
@@ -472,24 +505,35 @@ struct HfCorrelatorLcHadrons {
         }
         // Removing Lc daughters by checking track indices
         if ((candidate.prong0Id() == track.globalIndex()) || (candidate.prong1Id() == track.globalIndex()) || (candidate.prong2Id() == track.globalIndex())) {
-          continue;
+          if (!storeAutoCorrelationFlag) {
+            continue;
+          }
+          correlationStatus = true;
         }
         registry.fill(HIST("hPtParticleAssocMcRec"), track.pt());
 
+        if (correlateLcWithLeadingParticle) {
+          if (track.globalIndex() != leadingIndex) {
+            continue;
+          }
+        }
+
         if (candidate.isSelLcToPKPi() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(track.phi(), candidate.phi()),
                             track.eta() - candidate.eta(),
                             candidate.pt(),
                             track.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPKPi(candidate), isLcSignal);
         }
         if (candidate.isSelLcToPiKP() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(track.phi(), candidate.phi()),
                             track.eta() - candidate.eta(),
                             candidate.pt(),
                             track.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPiKP(candidate), isLcSignal);
         }
       } // end inner loop (Tracks)
@@ -506,6 +550,11 @@ struct HfCorrelatorLcHadrons {
     int counterLcHadron = 0;
     registry.fill(HIST("hMcEvtCount"), 0);
 
+    // find leading particle
+    if (correlateLcWithLeadingParticle) {
+      leadingIndex = findLeadingParticleMcGen(mcParticles, etaTrackMax.value, ptTrackMin.value);
+    }
+
     auto getTracksSize = [&mcParticles](aod::McCollision const& /*collision*/) {
       int nTracks = 0;
       for (const auto& track : mcParticles) {
@@ -556,15 +605,26 @@ struct HfCorrelatorLcHadrons {
           }
 
           if ((std::abs(particleAssoc.pdgCode()) != kElectron) && (std::abs(particleAssoc.pdgCode()) != kMuonMinus) && (std::abs(particleAssoc.pdgCode()) != kPiPlus) && (std::abs(particle.pdgCode()) != kKPlus) && (std::abs(particleAssoc.pdgCode()) != kProton)) {
-            continue;
+            if (!storeAutoCorrelationFlag) {
+              continue;
+            }
+            correlationStatus = true;
           }
+
+          if (correlateLcWithLeadingParticle) {
+            if (particleAssoc.globalIndex() != leadingIndex) {
+              continue;
+            }
+          }
+
           int poolBin = corrBinningMcGen.getBin(std::make_tuple(mcCollision.posZ(), getTracksSize(mcCollision)));
           registry.fill(HIST("hPtParticleAssocMcGen"), particleAssoc.pt());
           entryLcHadronPair(getDeltaPhi(particleAssoc.phi(), particle.phi()),
                             particleAssoc.eta() - particle.eta(),
                             particle.pt(),
                             particleAssoc.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(MassLambdaCPlus, true);
         } // end inner loop
       }
@@ -600,15 +660,17 @@ struct HfCorrelatorLcHadrons {
                             t1.eta() - t2.eta(),
                             t1.pt(),
                             t2.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPKPi(t1), false);
         }
         if (t1.isSelLcToPiKP() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(t1.phi(), t2.phi()),
                             t1.eta() - t2.eta(),
                             t1.pt(),
                             t2.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPiKP(t1), false);
         }
       }
@@ -634,15 +696,17 @@ struct HfCorrelatorLcHadrons {
                             t1.eta() - t2.eta(),
                             t1.pt(),
                             t2.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPKPi(t1), false);
         }
         if (t1.isSelLcToPiKP() >= selectionFlagLc) {
           entryLcHadronPair(getDeltaPhi(t1.phi(), t2.phi()),
                             t1.eta() - t2.eta(),
                             t1.pt(),
                             t2.pt(),
-                            poolBin);
+                            poolBin,
+                            correlationStatus);
           entryLcHadronRecoInfo(hfHelper.invMassLcToPiKP(t1), false);
         }
       }
@@ -696,7 +760,8 @@ struct HfCorrelatorLcHadrons {
                           t1.eta() - t2.eta(),
                           t1.pt(),
                           t2.pt(),
-                          poolBin);
+                          poolBin,
+                          correlationStatus);
       }
     }
   }