[PWGHF/HFC] D*-H correlation: Correcting the filling of D* and hadron table (#9145)

deependra170598 · web-flow · commit 1e27587074fe · 2025-01-07T13:39:41.000+01:00
diff --git a/PWGHF/HFC/DataModel/CorrelationTables.h b/PWGHF/HFC/DataModel/CorrelationTables.h
@@ -363,9 +363,11 @@ DECLARE_SOA_TABLE(TrkRecInfoDplus, "AOD", "TRKRECINFODPLUS", //! Tracks Reconstr
 // definition of columns and tables for Dstar-Hadron correlation pair
 namespace hf_correlation_dstar_hadron
 {
-DECLARE_SOA_INDEX_COLUMN(Collision, collision);
+DECLARE_SOA_INDEX_COLUMN(Collision, collision);      // used in pair table for indexing
+DECLARE_SOA_COLUMN(CollisionIdx, collisionIdx, int); // used in Dstar table for indexing
 // Dstar candidate properties
-DECLARE_SOA_INDEX_COLUMN(HfCandDstar, hfCandDstar);
+DECLARE_SOA_INDEX_COLUMN(HfCandDstar, hfCandDstar);      // used in pair table for indexing
+DECLARE_SOA_COLUMN(HfCandDstarIdx, hfCandDstarIdx, int); // used in Dstar table for indexing
 DECLARE_SOA_COLUMN(PhiDstar, phiDstar, float);
 DECLARE_SOA_COLUMN(EtaDstar, etaDstar, float);
 DECLARE_SOA_COLUMN(PtDstar, ptDstar, float);
@@ -410,9 +412,9 @@ DECLARE_SOA_TABLE(DstarHadronPair, "AOD", "DSTRHPAIR", // D* Hadrons pairs Infor
                   hf_correlation_dstar_hadron::DeltaM<hf_correlation_dstar_hadron::MDstar, hf_correlation_dstar_hadron::MD0>);
 
 DECLARE_SOA_TABLE(Dstar, "AOD", "DSTAR", // Only Dstar properties
-                  hf_correlation_dstar_hadron::CollisionId,
+                  hf_correlation_dstar_hadron::CollisionIdx,
                   // D* only properties
-                  hf_correlation_dstar_hadron::HfCandDstarId,
+                  hf_correlation_dstar_hadron::HfCandDstarIdx,
                   hf_correlation_dstar_hadron::PhiDstar,
                   hf_correlation_dstar_hadron::EtaDstar,
                   hf_correlation_dstar_hadron::PtDstar,
diff --git a/PWGHF/HFC/TableProducer/correlatorDstarHadrons.cxx b/PWGHF/HFC/TableProducer/correlatorDstarHadrons.cxx
@@ -13,6 +13,11 @@
 /// \author Deependra Sharma <deependra.sharma@cern.ch>, IITB
 /// \author Fabrizio Grosa <fabrizio.grosa@cern.ch>, CERN
 
+/// \brief Correlator for D* and hadrons. This task is used to produce table for D* and hadron pairs.
+
+// c++
+#include <vector>
+
 // O2
 #include "CommonConstants/PhysicsConstants.h"
 #include "Framework/ASoAHelpers.h"
@@ -26,11 +31,29 @@
 #include "PWGHF/DataModel/CandidateReconstructionTables.h"
 #include "PWGHF/DataModel/CandidateSelectionTables.h"
 #include "PWGHF/HFC/DataModel/CorrelationTables.h"
+#include "PWGHF/Utils/utilsAnalysis.h"
 
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
 
+const int nBinsPtCorrelation = 8;
+
+const double binsPtCorrelationsDefault[nBinsPtCorrelation + 1] = {0., 2., 4., 6., 8., 12., 16., 24., 100.};
+auto vecBinsPtCorrelationsDefault = std::vector<double>{binsPtCorrelationsDefault, binsPtCorrelationsDefault + nBinsPtCorrelation + 1};
+
+const double signalRegionLefBoundDefault[nBinsPtCorrelation] = {0.144, 0.144, 0.144, 0.144, 0.144, 0.144, 0.144, 0.144};
+auto vecSignalRegionLefBoundDefault = std::vector<double>{signalRegionLefBoundDefault, signalRegionLefBoundDefault + nBinsPtCorrelation};
+
+const double signalRegionRightBoundDefault[nBinsPtCorrelation] = {0.146, 0.146, 0.146, 0.146, 0.146, 0.146, 0.146, 0.146};
+auto vecSignalRegionRightBoundDefault = std::vector<double>{signalRegionRightBoundDefault, signalRegionRightBoundDefault + nBinsPtCorrelation};
+
+const double sidebandRightInnerDefault[nBinsPtCorrelation] = {0.147, 0.147, 0.147, 0.147, 0.147, 0.147, 0.147, 0.147};
+auto vecSidebandRightInnerDefault = std::vector<double>{sidebandRightInnerDefault, sidebandRightInnerDefault + nBinsPtCorrelation};
+
+const double sidebandRightOuterDefault[nBinsPtCorrelation] = {0.154, 0.154, 0.154, 0.154, 0.154, 0.154, 0.154, 0.154};
+auto vecSidebandRightOuterDefault = std::vector<double>{sidebandRightOuterDefault, sidebandRightOuterDefault + nBinsPtCorrelation};
+
 // flaging a collision if D* meson is found.
 struct HfCorrelatorDstarHadronsCollisionSelector {
   Produces<aod::DmesonSelection> collisionWDstar;
@@ -39,19 +62,19 @@ struct HfCorrelatorDstarHadronsCollisionSelector {
   Configurable<float> yCandMax{"yCandMax", 0.8, "max. cand. rapidity"};
   Configurable<float> ptCandMin{"ptCandMin", 1., "min. cand. pT"};
 
+  SliceCache cache;
+
   using DstarCandidates = soa::Join<aod::HfCandDstars, aod::HfSelDstarToD0Pi>;
   using FilteredCandidates = soa::Filtered<DstarCandidates>;
 
-  SliceCache cache;
-  Preslice<DstarCandidates> perColDstarCand = aod::hf_cand::collisionId;
-
   // candidates who passed the slection criteria defined in "CandidateSelectionTables.h"
   Filter candidateFilter = aod::hf_sel_candidate_dstar::isSelDstarToD0Pi == selectionFlagDstar;
 
+  Preslice<DstarCandidates> perColDstarCand = aod::hf_cand::collisionId;
+
   void processCollisionSelWDstar(aod::Collisions const& collisions,
                                  FilteredCandidates const& candidates)
   {
-
     for (const auto& collision : collisions) {
       bool isDstarFound = false;
       auto candidatesPerCol = candidates.sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
@@ -73,7 +96,7 @@ struct HfCorrelatorDstarHadronsCollisionSelector {
       } // candidate loop
       // LOG(info) << "processCollisionSelWDstar: isDstarFound = " << isDstarFound;
       collisionWDstar(isDstarFound); // compatible with collision table (filled collision by collision)
-    }                                // collision loop
+    } // collision loop
   }
   PROCESS_SWITCH(HfCorrelatorDstarHadronsCollisionSelector, processCollisionSelWDstar, "process only data for dstar hadron correlation", true);
 };
@@ -102,13 +125,20 @@ struct HfCorrelatorDstarHadrons {
   Configurable<float> ptAssoTrackMin{"ptAssoTrackMin", 0.5, "min Pt of Associated Track"};
   Configurable<float> ptAssoTrackMax{"ptAssoTrackMax", 50.0, "max pT of Associated Track"};
 
-  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"};
+  Configurable<std::vector<double>> binsPtCorrelations{"binsPtCorrelations", std::vector<double>{vecBinsPtCorrelationsDefault}, "pT bin limits for correlation plots"};
+  Configurable<std::vector<double>> signalRegionLefBound{"signalRegionLefBound", std::vector<double>{vecSignalRegionLefBoundDefault}, "left boundary of signal region vs pT"};
+  Configurable<std::vector<double>> signalRegionRightBound{"signalRegionRightBound", std::vector<double>{vecSignalRegionRightBoundDefault}, "right boundary of signal region vs pT"};
+  Configurable<std::vector<double>> rightSidebandOuterBoundary{"rightSidebandOuterBoundary", std::vector<double>{vecSidebandRightOuterDefault}, "right sideband outer baoundary vs pT"};
+  Configurable<std::vector<double>> rightSidebandInnerBoundary{"rightSidebandInnerBoundary", std::vector<double>{vecSidebandRightInnerDefault}, "right sideband inner boundary"};
 
-  // ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFV0M<aod::mult::MultFV0A, aod::mult::MultFV0C>> binningScheme{{binsZVtx, binsMultiplicity},true};
-  // ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFT0M<aod::mult::MultFT0A, aod::mult::MultFT0C>> binningScheme;
-  using BinningType = ColumnBinningPolicy<aod::collision::PosZ, aod::mult::MultFV0M<aod::mult::MultFV0A, aod::mult::MultFV0C>>;
-  BinningType binningScheme{{binsZVtx, binsMultiplicity}, true};
+  // Inv Mass of Dstar and D0 Candidate
+  float invMassDstarParticle;
+  float invMassD0Particle;
+  int binNumber;
+  SliceCache cache;
+
+  // 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>>;
 
   // Collision Table
   using CollisionsWDstar = soa::Join<aod::Collisions, aod::Mults, aod::DmesonSelection>;
@@ -128,12 +158,21 @@ struct HfCorrelatorDstarHadrons {
   Filter trackFilter = nabs(aod::track::eta) <= etaAbsAssoTrackMax && aod::track::pt >= ptAssoTrackMin && aod::track::pt <= ptAssoTrackMax &&
                        aod::track::dcaXY >= dcaxyAssoTrackMin && aod::track::dcaXY <= dcaxyAssoTrackMax &&
                        aod::track::dcaZ >= dcazAssoTrackMin && aod::track::dcaZ <= dcazAssoTrackMax;
-  SliceCache cache;
+
   // Preslice<DstarCandidates> perColCandidates = aod::hf_cand::collisionId;
   Preslice<FilteredCandidates> perColCandidates = aod::hf_cand::collisionId;
   // Preslice<aod::TracksWDca> perColTracks = aod::track::collisionId;
   Preslice<FilteredTracks> perColTracks = aod::track::collisionId;
 
+  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"};
+  BinningType binningScheme{{binsZVtx, binsMultiplicity}, true};
+  // Eta Phi Axes
+  ConfigurableAxis axisEta{"axisEta", {16, -1.0, 1.0}, "Eta Axis"};
+  ConfigurableAxis axisPhi{"axisPhi", {64, 0.0, 3.14}, "Phi Axis"};
+  ConfigurableAxis axisDeltaEta{"axisDeltaEta", {31, -2.0, 2.0}, "Delta Eta Axis"};
+  ConfigurableAxis axisDeltaPhi{"axisDeltaPhi", {64, -o2::constants::math::PIHalf, 3.0 * o2::constants::math::PIHalf}, "Delta Phi Axis"};
+
   HistogramRegistry registry{
     "registry",
     {{"hTriggerColCandPairCounts", "Counts of Trigger Collision, Trigger Candidates and Pair Counts", {HistType::kTH1F, {{3, 0.0, 3.0}}}}}};
@@ -145,20 +184,41 @@ struct HfCorrelatorDstarHadrons {
       LOGP(fatal, "One and only one process function must be enabled at a time.");
     }
 
+    invMassDstarParticle = -999.0;
+    invMassD0Particle = -999.0;
+    binNumber = -2;
+
     binningScheme = {{binsZVtx, binsMultiplicity}, true};
+
+    registry.add("QA/hCandsPerCol", "Candidates per Collision", {HistType::kTH1D, {{100, 0.0, 100.0}}});
+    registry.add("QA/hAssoTracksPerCol", "Tracks per Collision", {HistType::kTH1D, {{1000, 0.0, 1000.0}}});
+    registry.add("QA/hCandsVsTracksPerCol", "Candidates vs Tracks per Collision", {HistType::kTHnSparseF, {{100, 0.0, 100.0}, {1000, 0.0, 1000.0}}});
+    registry.add("QA/hCandsSignalVsTracksPerCol", "Candidates vs Tracks per Collision", {HistType::kTHnSparseF, {{100, 0.0, 100.0}, {1000, 0.0, 1000.0}}});
+    registry.add("QA/hCandsSideBandVsTracksPerCol", "Candidates vs Tracks per Collision", {HistType::kTHnSparseF, {{100, 0.0, 100.0}, {1000, 0.0, 1000.0}}});
+    // eta phi single particle distribution
+    registry.add("QA/hPhiDstarSignal", "Phi distribution of Dstar from signal region", {HistType::kTH1D, {axisPhi}});
+    registry.add("QA/hEtaDstarSignal", "Eta distribution of Dstar from signal region", {HistType::kTH1D, {axisEta}});
+    registry.add("QA/hPhiDstarSideBand", "Phi distribution of Dstar from side band region", {HistType::kTH1D, {axisPhi}});
+    registry.add("QA/hEtaDstarSideBand", "Eta distribution of Dstar from side band region", {HistType::kTH1D, {axisEta}});
+    registry.add("QA/hPhiAssoTrack", "Phi distribution of Associated Track", {HistType::kTH1D, {axisPhi}});
+    registry.add("QA/hEtaAssoTrack", "Eta distribution of Associated Track", {HistType::kTH1D, {axisEta}});
+    // delta eta phi distribution
+    registry.add("QA/hDPhiDstarAssoTrack", "Delta Phi distribution between Dstar and Associated Track", {HistType::kTH1D, {axisDeltaPhi}});
+    registry.add("QA/hDEtaDstarAssoTrack", "Delta Eta distribution between Dstar and Associated Track", {HistType::kTH1D, {axisDeltaEta}});
+    registry.add("QA/hDPhiDEtaDstarAssoTrack", "Delta Phi vs Delta Eta distribution between Dstar and Associated Track", {HistType::kTH2D, {axisDeltaPhi, axisDeltaEta}});
   }
 
   void processDataSameEvent(FilteredCollisions const& collisions, // only collisions who have altleast one D*
                             FilteredTracks const& tracks,
                             FilteredCandidates const& candidates,
                             aod::BCsWithTimestamps const&)
   {
-
     for (const auto& collision : collisions) {
       registry.fill(HIST("hTriggerColCandPairCounts"), 0); // counting trigger collision
 
       auto bc = collision.bc_as<aod::BCsWithTimestamps>();
       auto timestamp = bc.timestamp();
+      binNumber = binningScheme.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
       auto candidatesPerCol = candidates.sliceByCached(aod::hf_cand::collisionId, collision.globalIndex(), cache);
       auto tracksPerCol = tracks.sliceByCached(aod::track::collisionId, collision.globalIndex(), cache);
 
@@ -167,6 +227,64 @@ struct HfCorrelatorDstarHadrons {
       } // endif
 
       registry.fill(HIST("hTriggerColCandPairCounts"), 1); // counting number of trigger particle
+      registry.fill(HIST("QA/hCandsPerCol"), candidatesPerCol.size());
+      registry.fill(HIST("QA/hAssoTracksPerCol"), tracksPerCol.size());
+      registry.fill(HIST("QA/hCandsVsTracksPerCol"), candidatesPerCol.size(), tracksPerCol.size());
+
+      int nCandsSignal = 0;
+      int nCandsSideBand = 0;
+      // Single particle distribution for canfdidates
+      for (const auto& cand : candidatesPerCol) {
+        auto gItriggerParticle = cand.globalIndex();
+        if (cand.signSoftPi() > 0) {
+          invMassDstarParticle = cand.invMassDstar();
+          invMassD0Particle = cand.invMassD0();
+        } else {
+          invMassDstarParticle = cand.invMassAntiDstar();
+          invMassD0Particle = cand.invMassD0Bar();
+        }
+        auto ptDstar = cand.pt();
+        int corrBinPtDstar = o2::analysis::findBin(binsPtCorrelations, ptDstar);
+        auto deltaM = cand.invMassDstar() - cand.invMassD0();
+        if (deltaM > signalRegionLefBound->at(corrBinPtDstar) && deltaM < signalRegionRightBound->at(corrBinPtDstar)) {
+          // Signal Region
+          registry.fill(HIST("QA/hPhiDstarSignal"), cand.phi());
+          registry.fill(HIST("QA/hEtaDstarSignal"), cand.eta());
+          nCandsSignal++;
+        } else if (deltaM > rightSidebandInnerBoundary->at(corrBinPtDstar) && deltaM < rightSidebandOuterBoundary->at(corrBinPtDstar)) {
+          // Side Band Region
+          registry.fill(HIST("QA/hPhiDstarSideBand"), cand.phi());
+          registry.fill(HIST("QA/hEtaDstarSideBand"), cand.eta());
+          nCandsSideBand++;
+        }
+        registry.fill(HIST("QA/hCandsSignalVsTracksPerCol"), nCandsSignal, tracksPerCol.size());
+        registry.fill(HIST("QA/hCandsSideBandVsTracksPerCol"), nCandsSideBand, tracksPerCol.size());
+
+        if (enableSeparateTables) {
+          rowsDstar(collision.globalIndex(),
+                    gItriggerParticle,
+                    cand.phi(),
+                    cand.eta(),
+                    cand.pt(),
+                    invMassDstarParticle,
+                    invMassD0Particle,
+                    timestamp,
+                    binNumber);
+        }
+      } // Dstar loop
+      // Single particle distribution for tracks
+      for (const auto& track : tracksPerCol) {
+        registry.fill(HIST("QA/hPhiAssoTrack"), track.phi());
+        registry.fill(HIST("QA/hEtaAssoTrack"), track.eta());
+        if (enableSeparateTables) {
+          rowsAssoTrack(track.phi(),
+                        track.eta(),
+                        track.pt(),
+                        binNumber,
+                        collision.globalIndex(),
+                        timestamp);
+        }
+      } // Track loop
 
       // Pair creation
       for (const auto& [triggerParticle, assocParticle] : soa::combinations(soa::CombinationsFullIndexPolicy(candidatesPerCol, tracksPerCol))) {
@@ -176,7 +294,7 @@ struct HfCorrelatorDstarHadrons {
         // Track rejection based on daughter index
         if ((triggerParticle.prong0Id() == gIassocParticle) || (triggerParticle.prong1Id() == gIassocParticle) || (triggerParticle.prongPiId() == gIassocParticle)) {
           continue; // rejected pair if associated particle is same as any of daughter particle
-        }           // endif
+        } // endif
 
         // Trigger Particle Rejection
         if (triggerParticle.pt() > ptDstarMax || triggerParticle.pt() < ptDstarMin) {
@@ -188,12 +306,9 @@ struct HfCorrelatorDstarHadrons {
         } // endif
 
         registry.fill(HIST("hTriggerColCandPairCounts"), 2); // counting number of pairs
-
-        auto binNumber = binningScheme.getBin(std::make_tuple(collision.posZ(), collision.multFT0M()));
-
-        // Inv Mass of Dstar and D0 Candidate
-        float invMassDstarParticle = -999.;
-        float invMassD0Particle = -999.;
+        registry.fill(HIST("QA/hDPhiDstarAssoTrack"), triggerParticle.phi() - assocParticle.phi());
+        registry.fill(HIST("QA/hDEtaDstarAssoTrack"), triggerParticle.eta() - assocParticle.eta());
+        registry.fill(HIST("QA/hDPhiDEtaDstarAssoTrack"), triggerParticle.phi() - assocParticle.phi(), triggerParticle.eta() - assocParticle.eta());
 
         if (triggerParticle.signSoftPi() > 0) {
           invMassDstarParticle = triggerParticle.invMassDstar();
@@ -217,28 +332,7 @@ struct HfCorrelatorDstarHadrons {
                             assocParticle.pt(),
                             timestamp,
                             binNumber);
-
-        if (enableSeparateTables) {
-          rowsDstar(collision.globalIndex(),
-                    gItriggerParticle,
-                    triggerParticle.phi(),
-                    triggerParticle.eta(),
-                    triggerParticle.pt(),
-                    invMassDstarParticle,
-                    invMassD0Particle,
-                    timestamp,
-                    binNumber);
-
-          rowsAssoTrack(assocParticle.phi(),
-                        assocParticle.eta(),
-                        assocParticle.pt(),
-                        binNumber,
-                        collision.globalIndex(),
-                        timestamp);
-        }
-
       } // D-H pair loop
-
     } // collision loop
 
   } // processDataSameEvent
@@ -254,25 +348,17 @@ struct HfCorrelatorDstarHadrons {
     Pair<FilteredCollisions, FilteredCandidates, FilteredTracks, BinningType> pairData{binningScheme, 5, -1, collisions, dstarHadronTuple, &cache};
 
     for (const auto& [c1, candidatesPerCol, c2, tracksPerCol] : pairData) {
-
       auto bc = c2.bc_as<aod::BCsWithTimestamps>();
       auto timestamp = bc.timestamp();
-
       for (const auto& [triggerParticle, assocParticle] : soa::combinations(soa::CombinationsFullIndexPolicy(candidatesPerCol, tracksPerCol))) {
-
         auto gItriggerParticle = triggerParticle.globalIndex();
         auto gIassocParticle = assocParticle.globalIndex();
-
         auto yDstar = triggerParticle.y(constants::physics::MassDStar);
         if (std::abs(yDstar) > yAbsDstarMax) {
           continue;
         } // endif
 
-        int binNumber = binningScheme.getBin(std::make_tuple(c2.posZ(), c2.multFV0M()));
-
-        // Inv Mass of Dstar and D0 Candidate
-        float invMassDstarParticle = -999.;
-        float invMassD0Particle = -999.;
+        binNumber = binningScheme.getBin(std::make_tuple(c2.posZ(), c2.multFT0M()));
 
         if (triggerParticle.signSoftPi() > 0) {
           invMassDstarParticle = triggerParticle.invMassDstar();
