[PWGHF] add centrality selection and calculate pT ratio in Z-h analysis

PWGHF/HFL/Tasks/taskElectronWeakBoson.cxx

@@ -17,10 +17,13 @@
 #define HomogeneousField // o2-linter: disable=name/macro (required by KFParticle)
 #endif
 
+#include "PWGHF/Core/CentralityEstimation.h"
 #include "PWGJE/DataModel/EMCALClusters.h"
 
 #include "Common/Core/RecoDecay.h"
+#include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/EventSelection.h"
+#include "Common/DataModel/Multiplicity.h"
 #include "Common/DataModel/PIDResponseTPC.h"
 #include "Common/DataModel/TrackSelectionTables.h"
 #include "EventFiltering/Zorro.h"
@@ -58,6 +61,7 @@
 using namespace o2;
 using namespace o2::framework;
 using namespace o2::framework::expressions;
+using namespace o2::hf_centrality;
 
 struct HfTaskElectronWeakBoson {
 
@@ -121,6 +125,12 @@ struct HfTaskElectronWeakBoson {
   Configurable<int> kfConstructMethod{"kfConstructMethod", 2, "KF Construct Method"};
   Configurable<int> chiSqNdfMax{"chiSqNdfMax", 10, "Chi2 Max for mass reco by KF particle"};
 
+  // Centrality estimator configuration
+  Configurable<int> centralityEstimator{"centralityEstimator", CentralityEstimator::FT0M, "Centrality estimator. See CentralityEstimator for valid values."};
+  Configurable<bool> enableCentralityAnalysis{"enableCentralityAnalysis", true, "Enable centrality-dependent analysis"};
+  Configurable<float> centralityMin{"centralityMin", -1, "minimum cut on centrality selection"};
+  Configurable<float> centralityMax{"centralityMax", 101, "maximum cut on centrality selection"};
+
   // CCDB service object
   Service<o2::ccdb::BasicCCDBManager> ccdb;
 
@@ -142,7 +152,7 @@ struct HfTaskElectronWeakBoson {
       : pt(ptr), eta(e), phi(ph), mass(m), ptchild0(ptzee0), ptchild1(ptzee1), charge(ch) {}
   };
   std::vector<HfZeeCandidate> reconstructedZ;
-
+  using CollisionsWithCent = soa::Join<aod::Collisions, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As>;
   using SelectedClusters = o2::aod::EMCALClusters;
   // PbPb
   // using TrackEle = o2::soa::Join<o2::aod::Tracks, o2::aod::FullTracks, o2::aod::TracksExtra, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCFullEl>;
@@ -186,6 +196,10 @@ struct HfTaskElectronWeakBoson {
     if (cfgSkimmedProcessing) {
       zorroSummary.setObject(zorro.getZorroSummary());
     }
+    // check centrality
+    if (centralityEstimator < CentralityEstimator::FT0A || centralityEstimator > CentralityEstimator::FV0A) {
+      LOGF(fatal, "Invalid centrality estimator: %d", static_cast<int>(centralityEstimator.value));
+    }
 
     // add configurable for CCDB path
     zorro.setBaseCCDBPath(cfgCCDBPath.value);
@@ -217,11 +231,14 @@ struct HfTaskElectronWeakBoson {
     const AxisSpec axisPtHadron{50, 0, 50, "p_{T,hadron} (GeV/c)"};
     const AxisSpec axisPtZ{150, 0, 150, "p_{T,Z} (GeV/c)"};
     const AxisSpec axisSign{2, -2, 2, "charge sign"};
+    const AxisSpec axisCentrality{10, 0, 100, "Centrality (%)"};
+    const AxisSpec axisPtRatio{200, 0, 2.0, "pt ratio for h and Z"};
 
     // create registrygrams
     registry.add("hZvtx", "Z vertex", kTH1D, {axisZvtx});
     registry.add("hEventCounterInit", "hEventCounterInit", kTH1D, {axisCounter});
     registry.add("hEventCounter", "hEventCounter", kTH1D, {axisCounter});
+    registry.add("hCentrality", "Centrality distribution", kTH1D, {axisCentrality});
     registry.add("hITSchi2", "ITS #chi^{2}", kTH1F, {axisChi2});
     registry.add("hTPCchi2", "TPC #chi^{2}", kTH1F, {axisChi2});
     registry.add("hTPCnCls", "TPC NCls", kTH1F, {axisCluster});
@@ -243,15 +260,13 @@ struct HfTaskElectronWeakBoson {
     registry.add("hEMCtime", "EMC timing", kTH1F, {axisEMCtime});
     registry.add("hIsolationEnergy", "Isolation Energy", kTH2F, {{axisE}, {axisIsoEnergy}});
     registry.add("hIsolationTrack", "Isolation Track", kTH2F, {{axisE}, {axisIsoTrack}});
-    registry.add("hInvMassZeeLs", "invariant mass for Z LS pair", kTH2F, {{axisPt}, {axisInvMassZ}});
-    registry.add("hInvMassZeeUls", "invariant mass for Z ULS pair", kTH2F, {{axisPt}, {axisInvMassZ}});
-    registry.add("hKfInvMassZeeLs", "invariant mass for Z LS pair KFp", kTH2F, {{axisPt}, {axisInvMassZ}});
-    registry.add("hKfInvMassZeeUls", "invariant mass for Z ULS pair KFp", kTH2F, {{axisPt}, {axisInvMassZ}});
+    registry.add("hInvMassZee", "invariant mass for Z ULS pair", HistType::kTHnSparseF, {axisSign, axisPt, axisInvMassZ});
+    registry.add("hKfInvMassZee", "invariant mass for Z ULS pair KFp", HistType::kTHnSparseF, {axisSign, axisPt, axisInvMassZ});
     registry.add("hTHnElectrons", "electron info", HistType::kTHnSparseF, {axisPt, axisNsigma, axisM02, axisEop, axisIsoEnergy, axisIsoTrack, axisEta, axisDedx});
     registry.add("hTHnTrMatch", "Track EMC Match", HistType::kTHnSparseF, {axisPt, axisdPhi, axisdEta});
 
     // Z-hadron correlation histograms
-    registry.add("hZHadronDphi", "Z-hadron #Delta#phi correlation", HistType::kTHnSparseF, {axisSign, axisPtZ, axisDPhiZh});
+    registry.add("hZHadronDphi", "Z-hadron #Delta#phi correlation", HistType::kTHnSparseF, {axisSign, axisPtZ, axisDPhiZh, axisPtRatio, axisPtHadron});
     registry.add("hZptSpectrum", "Z boson p_{T} spectrum", kTH2F, {{axisSign}, {axisPtZ}});
 
     // hisotgram for EMCal trigger
@@ -323,6 +338,9 @@ struct HfTaskElectronWeakBoson {
     // LOG(info) << "Invarimass cal by KF particle ";
     for (const auto& track : tracks) {
 
+      if (std::abs(track.pt() - kfpIsoEle.GetPt()) < ptMatch) {
+        continue;
+      }
       if (track.pt() < ptZeeMin) {
         continue;
       }
@@ -344,6 +362,8 @@ struct HfTaskElectronWeakBoson {
       auto child2 = RecoDecayPtEtaPhi::pVector(kfpAssEle.GetPt() * correctionPtElectron, kfpAssEle.GetEta(), kfpAssEle.GetPhi());
       double invMassEE = RecoDecay::m(std::array{child1, child2}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron});
 
+      registry.fill(HIST("hInvMassZee"), track.sign() * charge, kfpIsoEle.GetPt(), invMassEE);
+
       // reco by KFparticle
       const KFParticle* electronPairs[2] = {&kfpIsoEle, &kfpAssEle};
       KFParticle zeeKF;
@@ -362,17 +382,9 @@ struct HfTaskElectronWeakBoson {
       }
       float massZee, massZeeErr;
       zeeKF.GetMass(massZee, massZeeErr);
+      registry.fill(HIST("hKfInvMassZee"), track.sign() * charge, kfpIsoEle.GetPt(), massZee);
       // LOG(info) << "Invarimass cal by KF particle mass = " << massZee;
       // LOG(info) << "Invarimass cal by RecoDecay = " << invMassEE;
-
-      if (track.sign() * charge > 0) {
-        registry.fill(HIST("hKfInvMassZeeLs"), kfpIsoEle.GetPt(), massZee);
-        registry.fill(HIST("hInvMassZeeLs"), kfpIsoEle.GetPt(), invMassEE);
-      } else {
-        registry.fill(HIST("hKfInvMassZeeUls"), kfpIsoEle.GetPt(), massZee);
-        registry.fill(HIST("hInvMassZeeUls"), kfpIsoEle.GetPt(), invMassEE);
-      }
-
       reconstructedZ.emplace_back(
         zeeKF.GetPt(),
         zeeKF.GetEta(),
@@ -384,7 +396,8 @@ struct HfTaskElectronWeakBoson {
     }
   }
 
-  void process(soa::Filtered<aod::Collisions>::iterator const& collision,
+  // void process(soa::Filtered<aod::Collisions>::iterator const& collision,
+  void process(soa::Filtered<CollisionsWithCent>::iterator const& collision,
                aod::BCsWithTimestamps const&,
                SelectedClusters const& emcClusters,
                TrackEle const& tracks,
@@ -446,6 +459,16 @@ struct HfTaskElectronWeakBoson {
 
     registry.fill(HIST("hZvtx"), collision.posZ());
 
+    // Calculate centrality
+    if (enableCentralityAnalysis) {
+      float centrality = o2::hf_centrality::getCentralityColl(collision, centralityEstimator);
+      // LOG(info) << centrality;
+      if (centrality < centralityMin || centrality > centralityMax) {
+        return;
+      }
+      registry.fill(HIST("hCentrality"), centrality);
+    }
+
     for (const auto& track : tracks) {
 
       if (std::abs(track.eta()) > etaTrMax) {
@@ -562,8 +585,9 @@ struct HfTaskElectronWeakBoson {
             }
             // LOG(info) << "E/p" << eop;
             registry.fill(HIST("hEopNsigTPC"), match.track_as<TrackEle>().tpcNSigmaEl(), eop);
-            if (match.emcalcluster_as<SelectedClusters>().m02() < m02Min || match.emcalcluster_as<SelectedClusters>().m02() > m02Max)
+            if (match.emcalcluster_as<SelectedClusters>().m02() < m02Min || match.emcalcluster_as<SelectedClusters>().m02() > m02Max) {
               continue;
+            }
 
             if (match.track_as<TrackEle>().tpcNSigmaEl() > nsigTpcMin && match.track_as<TrackEle>().tpcNSigmaEl() < nsigTpcMax) {
               registry.fill(HIST("hEop"), match.track_as<TrackEle>().pt(), eop);
@@ -624,7 +648,8 @@ struct HfTaskElectronWeakBoson {
           }
           // calculate Z-h correlation
           double deltaPhi = RecoDecay::constrainAngle(trackAss.phi - zBoson.phi, -o2::constants::math::PIHalf);
-          registry.fill(HIST("hZHadronDphi"), zBoson.charge, zBoson.pt, deltaPhi);
+          double ptRatio = trackAss.pt / zBoson.pt;
+          registry.fill(HIST("hZHadronDphi"), zBoson.charge, zBoson.pt, deltaPhi, ptRatio, trackAss.pt);
         }
       }
     } // end of Z-hadron correlation