Add centrality dependence to D0–hadron correlator and correlation task

 - Introduced centrality-dependent analysis in correlator and correlation task
 - Relaxed ML selection: changed '||' to '&&' in D0/D0bar BDT condition to improve efficiency

Author: ayankrmanna

PWGHF/HFC/DataModel/CorrelationTables.h

@@ -72,6 +72,7 @@ DECLARE_SOA_COLUMN(TrackDcaXY, trackDcaXY, float);                         //! D
 DECLARE_SOA_COLUMN(TrackDcaZ, trackDcaZ, float);                           //! DCA z of the track
 DECLARE_SOA_COLUMN(TrackTPCNClsCrossedRows, trackTPCNClsCrossedRows, int); //! Number of crossed TPC Rows
 DECLARE_SOA_COLUMN(IsAutoCorrelated, isAutoCorrelated, bool);              //! Correlation Status
+DECLARE_SOA_COLUMN(Cent, cent, float);                                     //! Centrality of Collision
 DECLARE_SOA_COLUMN(TrackOrigin, trackOrigin, int);                         //! Check track origin
 DECLARE_SOA_COLUMN(IsPrompt, isPrompt, bool);                              //! Used in MC-Rec, D0 Prompt or Non-Prompt
 DECLARE_SOA_COLUMN(IsPhysicalPrimary, isPhysicalPrimary, bool);            //! Used in MC-Rec, primary associated particles
@@ -102,7 +103,8 @@ DECLARE_SOA_TABLE(D0HadronPair, "AOD", "D0HPAIR", //! D0-Hadrons pairs Informati
                   aod::hf_correlation_d0_hadron::PtD,
                   aod::hf_correlation_d0_hadron::PtHadron,
                   aod::hf_correlation_d0_hadron::PoolBin,
-                  aod::hf_correlation_d0_hadron::IsAutoCorrelated);
+                  aod::hf_correlation_d0_hadron::IsAutoCorrelated,
+                  aod::hf_correlation_d0_hadron::Cent);
 
 DECLARE_SOA_TABLE(D0HadronRecoInfo, "AOD", "D0HRECOINFO", //! D0-Hadrons pairs Reconstructed Informations
                   aod::hf_correlation_d0_hadron::MD,

PWGHF/HFC/TableProducer/correlatorD0Hadrons.cxx

@@ -28,6 +28,7 @@
 
 #include "Common/CCDB/EventSelectionParams.h"
 #include "Common/Core/RecoDecay.h"
+#include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/EventSelection.h"
 #include "Common/DataModel/Multiplicity.h"
 #include "Common/DataModel/TrackSelectionTables.h"
@@ -76,7 +77,7 @@ double getDeltaPhi(double phiHadron, double phiD)
 // Types
 using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFT0M<aod::mult::MultFT0A, aod::mult::MultFT0C>>;
 using BinningTypeMcGen = ColumnBinningPolicy<aod::mccollision::PosZ, o2::aod::mult::MultMCFT0A>;
-using SelectedCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::Mults, aod::EvSels, aod::DmesonSelection>>;
+using SelectedCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::Mults, aod::EvSels, aod::DmesonSelection, aod::CentFT0Ms>>;
 using SelectedTracks = soa::Filtered<soa::Join<aod::TracksWDca, aod::TrackSelection, aod::TracksExtra>>;
 using SelectedCandidatesData = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0>>;
 using SelectedCandidatesDataMl = soa::Filtered<soa::Join<aod::HfCand2Prong, aod::HfSelD0, aod::HfMlD0>>;
@@ -96,10 +97,13 @@ struct HfCorrelatorD0HadronsSelection {
   Configurable<int> selectionFlagD0bar{"selectionFlagD0bar", 1, "Selection Flag for D0bar"};
   Configurable<float> yCandMax{"yCandMax", 4.0, "max. cand. rapidity"};
   Configurable<float> ptCandMin{"ptCandMin", -1., "min. cand. pT"};
+  Configurable<float> centMin{"centMin", 0., "Minimum Centrality"};
+  Configurable<float> centMax{"centMax", 100., "Maximum Centrality"};
+  Configurable<bool> useCentrality{"useCentality", false, "Flag for centrality dependent analyses"};
 
   SliceCache cache;
 
-  using SelCollisions = soa::Join<aod::Collisions, aod::EvSels>;
+  using SelCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Ms>;
 
   Preslice<aod::HfCand2Prong> perCol = aod::hf_cand::collisionId;
 
@@ -113,6 +117,7 @@ struct HfCorrelatorD0HadronsSelection {
     bool isD0Found = true;
     bool isSel8 = true;
     bool isNosameBunchPileUp = true;
+    bool isCentInRange = false;
     if (selectedD0Candidates.size() > 0) {
       auto selectedD0CandidatesGrouped = selectedD0Candidates->sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
 
@@ -130,13 +135,18 @@ struct HfCorrelatorD0HadronsSelection {
         break;
       }
     }
+    float cent = 0.;
+    if (useCentrality) {
+        cent = collision.centFT0M();
+    }
     if (useSel8) {
       isSel8 = collision.sel8();
     }
     if (selNoSameBunchPileUpColl) {
       isNosameBunchPileUp = collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup);
     }
-    isSelColl = isD0Found && isSel8 && isNosameBunchPileUp;
+    isCentInRange = (cent >= centMin && cent <=centMax);
+    isSelColl = isD0Found && isSel8 && isNosameBunchPileUp && isCentInRange;
     collisionsWithSelD0(isSelColl);
   }
   PROCESS_SWITCH(HfCorrelatorD0HadronsSelection, processD0SelectionData, "Process D0 Selection Data", false);
@@ -226,6 +236,7 @@ struct HfCorrelatorD0Hadrons {
   Configurable<bool> correlateD0WithLeadingParticle{"correlateD0WithLeadingParticle", false, "Switch for correlation of D0 mesons with leading particle only"};
   Configurable<bool> storeAutoCorrelationFlag{"storeAutoCorrelationFlag", false, "Store flag that indicates if the track is paired to its D-meson mother instead of skipping it"};
   Configurable<int> numberEventsMixed{"numberEventsMixed", 5, "Number of events mixed in ME process"};
+  Configurable<bool> useCentrality{"useCentrality", false, "Flag for centrality dependent analyses"};
 
   int leadingIndex = 0;
   double massD0{0.};
@@ -255,6 +266,7 @@ struct HfCorrelatorD0Hadrons {
   ConfigurableAxis binsMultFT0M{"binsMultFT0M", {10000, 0., 10000.}, "Multiplicity as FT0M signal amplitude"};
   ConfigurableAxis binsPosZ{"binsPosZ", {100, -10., 10.}, "primary vertex z coordinate"};
   ConfigurableAxis binsPoolBin{"binsPoolBin", {9, 0., 9.}, "PoolBin"};
+  ConfigurableAxis binsCentFt0m{"binsCentFT0m", {100, 0., 100.}, "centrality percentile (FT0M)"};
 
   BinningType corrBinning{{zPoolBins, multPoolBins}, true};
 
@@ -283,6 +295,7 @@ struct HfCorrelatorD0Hadrons {
     AxisSpec axisBdtScoreBkg = {100, 0., 1., "Bdt score background"};
     AxisSpec axisBdtScorePrompt = {100, 0., 1., "Bdt score prompt"};
     AxisSpec axisOrigin = {10, 0., 10., "Candidate origin"};
+    AxisSpec axisCent = {binsCentFt0m,"Centrality"};    
 
     // Histograms for Data
     registry.add("hPtCand", "D0, D0bar candidates", {HistType::kTH1F, {axisPtD}});
@@ -292,12 +305,14 @@ struct HfCorrelatorD0Hadrons {
     registry.add("hEta", "D0,D0bar candidates", {HistType::kTH1F, {axisEta}});
     registry.add("hPhi", "D0,D0bar candidates", {HistType::kTH1F, {axisPhi}});
     registry.add("hY", "D0,D0bar candidates", {HistType::kTH1F, {axisRapidity}});
+    registry.add("hCentFT0M","Centrality FT0M;centrality;entries", {HistType::kTH1F,{{100, 0., 100.}}});
     registry.add("hMultiplicityPreSelection", "multiplicity prior to selection;multiplicity;entries", {HistType::kTH1F, {axisMultiplicity}});
-    registry.add("hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {axisMultiplicity}});
+    registry.add("hMultiplicity", "multiplicity;multiplicity;entries", {HistType::kTH1F, {axisMultiplicity}});    
     registry.add("hMass", "D0, D0bar candidates massVsPt", {HistType::kTH2F, {{axisMassD}, {axisPtD}}});
     registry.add("hMass1D", "D0, D0bar candidates mass", {HistType::kTH1F, {axisMassD}});
     registry.add("hMassD01D", "D0 candidates mass", {HistType::kTH1F, {axisMassD}});
     registry.add("hMassD0bar1D", "D0bar candidates mass", {HistType::kTH1F, {axisMassD}});
+    registry.add("hMassD0VsPtVsCent", "D0 candidates;inv. mass (p K #pi) (GeV/#it{c}^{2});entries", {HistType::kTH3F, {{axisMassD}, {axisPtD}, {axisCent}}});
     registry.add("hMLScoresVsMassVsPtVsOrigin", "D0, D0bar candidates massVsPt", {HistType::kTHnSparseD, {{axisBdtScoreBkg}, {axisBdtScorePrompt}, {axisMassD}, {axisPtD}, {axisOrigin}}});
     // Histograms for MC Reco
     registry.add("hPtCandRec", "D0, D0bar candidates - MC reco", {HistType::kTH1F, {axisPtD}});
@@ -354,6 +369,10 @@ struct HfCorrelatorD0Hadrons {
     if (correlateD0WithLeadingParticle) {
       leadingIndex = findLeadingParticle(tracks, etaTrackMax.value);
     }
+    float cent = 0.;
+    if (useCentrality){
+      cent = collision.centFT0M();
+    }
 
     int poolBin = corrBinning.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
     registry.fill(HIST("hCollisionPoolBin"), poolBin);
@@ -404,6 +423,7 @@ struct HfCorrelatorD0Hadrons {
         registry.fill(HIST("hMass"), HfHelper::invMassD0ToPiK(candidate), candidate.pt(), efficiencyWeight);
         registry.fill(HIST("hMass1D"), HfHelper::invMassD0ToPiK(candidate), efficiencyWeight);
         registry.fill(HIST("hMassD01D"), HfHelper::invMassD0ToPiK(candidate), efficiencyWeight);
+        registry.fill(HIST("hMassD0VsPtVsCent"), HfHelper::invMassD0ToPiK(candidate), candidate.pt(), cent, efficiencyWeight);
         for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) {
           outputMlD0[iclass] = candidate.mlProbD0()[classMl->at(iclass)];
         }
@@ -413,6 +433,7 @@ struct HfCorrelatorD0Hadrons {
         registry.fill(HIST("hMass"), HfHelper::invMassD0barToKPi(candidate), candidate.pt(), efficiencyWeight);
         registry.fill(HIST("hMass1D"), HfHelper::invMassD0barToKPi(candidate), efficiencyWeight);
         registry.fill(HIST("hMassD0bar1D"), HfHelper::invMassD0barToKPi(candidate), efficiencyWeight);
+        registry.fill(HIST("hMassD0VsPtVsCent"), HfHelper::invMassD0barToKPi(candidate), candidate.pt(), cent, efficiencyWeight);
         for (unsigned int iclass = 0; iclass < classMl->size(); iclass++) {
           outputMlD0bar[iclass] = candidate.mlProbD0bar()[classMl->at(iclass)];
         }
@@ -484,11 +505,13 @@ struct HfCorrelatorD0Hadrons {
                           candidate.pt(),
                           track.pt(),
                           poolBin,
-                          correlationStatus);
+                          correlationStatus,
+                          cent);
         entryD0HadronRecoInfo(HfHelper::invMassD0ToPiK(candidate), HfHelper::invMassD0barToKPi(candidate), signalStatus);
         entryD0HadronGenInfo(false, false, 0);
         entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]);
         entryTrackRecoInfo(track.dcaXY(), track.dcaZ(), track.tpcNClsCrossedRows());
+        registry.fill(HIST("hCentFT0M"),cent );
 
       } // end inner loop (tracks)
 
@@ -614,7 +637,7 @@ struct HfCorrelatorD0Hadrons {
 
       flagD0 = candidate.flagMcMatchRec() == o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK;     // flagD0Signal 'true' if candidate matched to D0 (particle)
       flagD0bar = candidate.flagMcMatchRec() == -o2::hf_decay::hf_cand_2prong::DecayChannelMain::D0ToPiK; // flagD0Reflection 'true' if candidate, selected as D0 (particle), is matched to D0bar (antiparticle)
-
+      float cent = 100.0;  //Centrality Placeholder: will be updated later
       // ========== track loop starts here ========================
 
       for (const auto& track : tracks) {
@@ -688,7 +711,8 @@ struct HfCorrelatorD0Hadrons {
                           candidate.pt(),
                           track.pt(),
                           poolBin,
-                          correlationStatus);
+                          correlationStatus,
+                          cent);
         entryD0HadronRecoInfo(HfHelper::invMassD0ToPiK(candidate), HfHelper::invMassD0barToKPi(candidate), signalStatus);
         entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]);
         if (track.has_mcParticle()) {
@@ -729,6 +753,7 @@ struct HfCorrelatorD0Hadrons {
     bool isD0Prompt = false;
     bool isD0NonPrompt = false;
     int trackOrigin = -1;
+    float cent = 100.; //Centrality Placeholder: will be updated later
 
     for (const auto& particleTrigg : mcParticles) {
       if (std::abs(particleTrigg.pdgCode()) != Pdg::kD0) {
@@ -803,7 +828,8 @@ struct HfCorrelatorD0Hadrons {
                             particleTrigg.pt(),
                             particleAssoc.pt(),
                             poolBin,
-                            correlationStatus);
+                            correlationStatus,
+                            cent);
           entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info
           entryD0HadronGenInfo(isD0Prompt, particleAssoc.isPhysicalPrimary(), trackOrigin);
         } // end inner loop (Tracks)
@@ -853,6 +879,7 @@ struct HfCorrelatorD0Hadrons {
         invMassDstar2 = std::sqrt((eKPi + ePion) * (eKPi + ePion) - pSum2);
         std::vector<float> outputMlD0 = {-1., -1., -1.};
         std::vector<float> outputMlD0bar = {-1., -1., -1.};
+        float cent = 100.; //Centrality Placeholder: will be updated later
 
         if (candidate.isSelD0() >= selectionFlagD0) {
           if ((std::abs(invMassDstar1 - HfHelper::invMassD0ToPiK(candidate)) - softPiMass) < ptSoftPionMax) {
@@ -887,7 +914,7 @@ struct HfCorrelatorD0Hadrons {
           }
         }
         bool correlationStatus = false;
-        entryD0HadronPair(getDeltaPhi(candidate.phi(), particleAssoc.phi()), candidate.eta() - particleAssoc.eta(), candidate.pt(), particleAssoc.pt(), poolBin, correlationStatus);
+        entryD0HadronPair(getDeltaPhi(candidate.phi(), particleAssoc.phi()), candidate.eta() - particleAssoc.eta(), candidate.pt(), particleAssoc.pt(), poolBin, correlationStatus, cent);
         entryD0HadronRecoInfo(HfHelper::invMassD0ToPiK(candidate), HfHelper::invMassD0barToKPi(candidate), signalStatus);
         entryD0HadronGenInfo(false, false, 0);
         entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]);
@@ -918,6 +945,7 @@ struct HfCorrelatorD0Hadrons {
       registry.fill(HIST("hD0PoolBin"), poolBinD0);
       registry.fill(HIST("hMultFT0M"), c1.multFT0M());
       registry.fill(HIST("hZvtx"), c1.posZ());
+      float cent = 100.; //Centrality Placeholder: will be updated later
 
       for (const auto& [candidate, particleAssoc] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
 
@@ -1019,7 +1047,7 @@ struct HfCorrelatorD0Hadrons {
         } // background case D0bar
         registry.fill(HIST("hSignalStatusMERec"), signalStatus);
         bool correlationStatus = false;
-        entryD0HadronPair(getDeltaPhi(candidate.phi(), particleAssoc.phi()), candidate.eta() - particleAssoc.eta(), candidate.pt(), particleAssoc.pt(), poolBin, correlationStatus);
+        entryD0HadronPair(getDeltaPhi(candidate.phi(), particleAssoc.phi()), candidate.eta() - particleAssoc.eta(), candidate.pt(), particleAssoc.pt(), poolBin, correlationStatus, cent);
         entryD0HadronRecoInfo(HfHelper::invMassD0ToPiK(candidate), HfHelper::invMassD0barToKPi(candidate), signalStatus);
         entryD0HadronGenInfo(isD0Prompt, isPhysicalPrimary, trackOrigin);
         entryD0HadronMlInfo(outputMlD0[0], outputMlD0[1], outputMlD0[2], outputMlD0bar[0], outputMlD0bar[1], outputMlD0bar[2]);
@@ -1064,10 +1092,11 @@ struct HfCorrelatorD0Hadrons {
           if (std::abs(particleAssoc.pdgCode()) == kPiPlus && indexMotherPi >= 0 && indexMotherD0 >= 0 && indexMotherPi == indexMotherD0) {
             continue;
           }
+          float cent = 100.; //Centrality Placeholder: will be updated later
           bool correlationStatus = false;
           int trackOrigin = RecoDecay::getCharmHadronOrigin(mcParticles, particleAssoc, true);
           bool isD0Prompt = particleTrigg.originMcGen() == RecoDecay::OriginType::Prompt;
-          entryD0HadronPair(getDeltaPhi(particleAssoc.phi(), particleTrigg.phi()), particleAssoc.eta() - particleTrigg.eta(), particleTrigg.pt(), particleAssoc.pt(), poolBin, correlationStatus);
+          entryD0HadronPair(getDeltaPhi(particleAssoc.phi(), particleTrigg.phi()), particleAssoc.eta() - particleTrigg.eta(), particleTrigg.pt(), particleAssoc.pt(), poolBin, correlationStatus, cent);
           entryD0HadronRecoInfo(massD0, massD0, 0); // dummy info
           entryD0HadronGenInfo(isD0Prompt, particleAssoc.isPhysicalPrimary(), trackOrigin);
         }