[PWGCF] Fixes in femto frame work (#18578)

Author: ariedel-cern

PWGCF/Femto/Core/closePairRejection.h

@@ -17,6 +17,7 @@
 #define PWGCF_FEMTO_CORE_CLOSEPAIRREJECTION_H_
 
 #include "PWGCF/Femto/Core/histManager.h"
+#include "PWGCF/Femto/Core/modes.h"
 
 #include "Common/Core/RecoDecay.h"
 
@@ -61,7 +62,8 @@ enum CprHist {
 
 // template configurable group for Cpr
 // DefaultCutAverage: default of cutAverage (e.g. false for pairs of generated particles, where no detector effects are present)
-template <auto& Prefix, bool DefaultCutAverage = true>
+// DefaultKinematicVariable: default of kinematicVariable (e.g. Q3 for triplets, kstar for pairs)
+template <auto& Prefix, bool DefaultCutAverage = true, modes::KinematicVariable DefaultKinematicVariable = modes::KinematicVariable::kKstar>
 struct ConfCpr : o2::framework::ConfigurableGroup {
   std::string prefix = std::string(Prefix);
   o2::framework::Configurable<bool> cutAverage{"cutAverage", DefaultCutAverage, "Apply CPR if the average deta-dphistar is below the configured values"};
@@ -71,6 +73,7 @@ struct ConfCpr : o2::framework::ConfigurableGroup {
   o2::framework::Configurable<bool> plotAverage{"plotAverage", true, "Plot average deta dphi distribution"};
   o2::framework::Configurable<bool> plotAngularCorrelation{"plotAngularCorrelation", false, "Plot angular correlation of particles (eta1 vs eta2 & phi1 vs phi2)"};
   o2::framework::Configurable<bool> plotKinematic{"plotKinematic", true, "Plot kinematic (kstar/Q3) distribution of blocked pairs/triplets"};
+  o2::framework::Configurable<int> kinematicVariable{"kinematicVariable", static_cast<int>(DefaultKinematicVariable), "Type of kinematic variable (0->kstar, 1->kT, 2->mT, 3->Q3) reported by the kinematic plot"};
   o2::framework::Configurable<float> detaMax{"detaMax", 0.01f, "Maximium deta"};
   o2::framework::Configurable<float> dphistarMax{"dphistarMax", 0.01f, "Maximum dphistar"};
   o2::framework::Configurable<float> detaCenter{"detaCenter", 0.f, "Center of deta cut"};
@@ -241,6 +244,7 @@ class CloseTrackRejection
 
     mPlotAngularCorrelation = confCpr.plotAngularCorrelation.value;
     mPlotKinematic = confCpr.plotKinematic.value;
+    mKinematicVariable = static_cast<modes::KinematicVariable>(confCpr.kinematicVariable.value);
 
     if (confCpr.seed.value >= 0) {
       uint64_t randomSeed = 0;
@@ -259,7 +263,8 @@ class CloseTrackRejection
     mHistogramRegistry = registry;
 
     if (mPlotKinematic) {
-      mHistogramRegistry->add(std::string(prefix) + getHistNameV2(kKinematic, HistTable), getHistDesc(kKinematic, HistTable), getHistType(kKinematic, HistTable), {specs.at(kKinematic)});
+      std::string kinematicDesc = histmanager::replaceFirstAxisTitle(getHistDesc(kKinematic, HistTable), modes::kinematicVariableToAxisTitle(mKinematicVariable));
+      mHistogramRegistry->add(std::string(prefix) + getHistNameV2(kKinematic, HistTable), kinematicDesc.c_str(), getHistType(kKinematic, HistTable), {specs.at(kKinematic)});
     }
 
     if (mPlotAverage) {
@@ -301,9 +306,9 @@ class CloseTrackRejection
     bool isClose = mUseEpisoidalCut ? isClosePairEpisoidalCut() : isClosePairRectengularCut();
 
     if (isClose) {
-      fillBlocked(pairHistManager.getKinematic());
+      fillBlocked(pairHistManager.getKinematic(mKinematicVariable));
     } else {
-      fillPassed(pairHistManager.getKinematic());
+      fillPassed(pairHistManager.getKinematic(mKinematicVariable));
     }
 
     return isClose;
@@ -493,6 +498,7 @@ class CloseTrackRejection
   bool mPlotAverage = false;
   bool mPlotAngularCorrelation = false;
   bool mPlotKinematic = false;
+  modes::KinematicVariable mKinematicVariable = modes::KinematicVariable::kKstar;
 
   float mKinematicMin = -1.f;
   float mKinematicMax = -1.f;
@@ -531,10 +537,10 @@ template <auto& prefix>
 class ClosePairRejectionTrackTrack
 {
  public:
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr,
+            T1 const& confCpr,
             int absChargeTrack1,
             int absChargeTrack2)
   {
@@ -597,10 +603,10 @@ template <auto& prefixTrackV0>
 class ClosePairRejectionTrackV0 // can also be used for any particle type that has pos/neg daughters, like resonances
 {
  public:
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr,
+            T1 const& confCpr,
             int absChargeTrack)
   {
     mCtr.init(registry, specs, confCpr, absChargeTrack, 1);
@@ -719,10 +725,10 @@ template <auto& prefix>
 class ClosePairRejectionTrackKink
 {
  public:
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr,
+            T1 const& confCpr,
             int absChargeTrack)
   {
     mCtr.init(registry, specs, confCpr, absChargeTrack, 1);
@@ -753,10 +759,10 @@ template <auto& prefix>
 class ClosePairRejectionTrackTwoPhotonResonance
 {
  public:
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* /*registry*/,
             std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& /*specs*/,
-            T const& /*confCpr*/,
+            T1 const& /*confCpr*/,
             int /*absChargeTrack*/)
   {
   }
@@ -774,10 +780,10 @@ template <auto& prefix>
 class ClosePairRejectionMcParticleMcParticle
 {
  public:
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr)
+            T1 const& confCpr)
   {
     mCtr.init(registry, specs, confCpr, 1, 1);
     mCtr.setMagField(confCpr.magField.value);

PWGCF/Femto/Core/closeTripletRejection.h

@@ -17,6 +17,7 @@
 #define PWGCF_FEMTO_CORE_CLOSETRIPLETREJECTION_H_
 
 #include "PWGCF/Femto/Core/closePairRejection.h"
+#include "PWGCF/Femto/Core/modes.h"
 
 #include <Framework/HistogramRegistry.h>
 #include <Framework/HistogramSpec.h>
@@ -31,11 +32,12 @@ constexpr const char PrefixCtrTrackTrackV0[] = "CtrTrackTrackV0";
 constexpr const char PrefixCtrTrackTrackCascade[] = "CtrTrackTrackCascade";
 constexpr const char PrefixCtrMcParticleMcParticleMcParticle[] = "CtrMcParticleMcParticleMcParticle";
 
-using ConfCtrTrackTrackTrack = closepairrejection::ConfCpr<PrefixCtrTrackTrackTrack>;
-using ConfCtrTrackTrackV0 = closepairrejection::ConfCpr<PrefixCtrTrackTrackV0>;
-using ConfCtrTrackTrackCascade = closepairrejection::ConfCpr<PrefixCtrTrackTrackCascade>;
+// Q3 (not kstar) is the natural kinematic variable for triplets, so it is the default here
+using ConfCtrTrackTrackTrack = closepairrejection::ConfCpr<PrefixCtrTrackTrackTrack, true, modes::KinematicVariable::kQ3>;
+using ConfCtrTrackTrackV0 = closepairrejection::ConfCpr<PrefixCtrTrackTrackV0, true, modes::KinematicVariable::kQ3>;
+using ConfCtrTrackTrackCascade = closepairrejection::ConfCpr<PrefixCtrTrackTrackCascade, true, modes::KinematicVariable::kQ3>;
 // generated particles, no cut by default (close triplet rejection corrects a detector effect)
-using ConfCtrMcParticleMcParticleMcParticle = closepairrejection::ConfCpr<PrefixCtrMcParticleMcParticleMcParticle, false>;
+using ConfCtrMcParticleMcParticleMcParticle = closepairrejection::ConfCpr<PrefixCtrMcParticleMcParticleMcParticle, false, modes::KinematicVariable::kQ3>;
 
 // directory names
 constexpr char PrefixTrack1Track2Se[] = "CPR_Track1Track2/SE/";
@@ -81,10 +83,10 @@ class CloseTripletRejectionTrackTrackTrack
   CloseTripletRejectionTrackTrackTrack() = default;
   ~CloseTripletRejectionTrackTrackTrack() = default;
 
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<closepairrejection::CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr,
+            T1 const& confCpr,
             int absChargeTrack1,
             int absChargeTrack2,
             int absChargeTrack3)
@@ -128,10 +130,10 @@ class CloseTripletRejectionTrackTrackV0
   CloseTripletRejectionTrackTrackV0() = default;
   ~CloseTripletRejectionTrackTrackV0() = default;
 
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<closepairrejection::CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr,
+            T1 const& confCpr,
             int absChargeTrack1,
             int absChargeTrack2)
   {
@@ -225,10 +227,10 @@ class CloseTripletRejectionMcParticleMcParticleMcParticle
   CloseTripletRejectionMcParticleMcParticleMcParticle() = default;
   ~CloseTripletRejectionMcParticleMcParticleMcParticle() = default;
 
-  template <typename T>
+  template <typename T1>
   void init(o2::framework::HistogramRegistry* registry,
             std::map<closepairrejection::CprHist, std::vector<o2::framework::AxisSpec>> const& specs,
-            T const& confCpr)
+            T1 const& confCpr)
   {
     mCtr12.init(registry, specs, confCpr);
     mCtr23.init(registry, specs, confCpr);

PWGCF/Femto/Core/histManager.h

@@ -85,5 +85,19 @@ constexpr const char* getHistDesc(EnumType variable, const ArrayType& array)
 
   return it != array.end() ? it->histdesc.data() : "";
 }
+
+// Replace the x-axis segment of a "title; x axis; y axis[; z axis...]" histogram description
+// with newTitle. Used to patch a static HistTable description at runtime once the actual axis
+// quantity (e.g. which momentum type, which kinematic variable) is resolved from a Configurable.
+inline std::string replaceFirstAxisTitle(std::string_view desc, std::string_view newTitle)
+{
+  std::string result(desc);
+  const auto first = result.find(';');
+  const auto second = (first == std::string::npos) ? std::string::npos : result.find(';', first + 1);
+  if (second != std::string::npos) {
+    result.replace(first + 1, second - first - 1, " " + std::string(newTitle) + " ");
+  }
+  return result;
+}
 } // namespace o2::analysis::femto::histmanager
 #endif // PWGCF_FEMTO_CORE_HISTMANAGER_H_

PWGCF/Femto/Core/mcBuilder.h

@@ -673,25 +673,11 @@ class McBuilder
     if (!mcParticle.has_mothers()) {
       return -1;
     }
-    auto motherIds = mcParticle.mothersIds();
-    // adapted these checks from MCUtils in PWGEM
-    const int defaultMotherSize = 2;
-    std::vector<int> allMotherIds;
-    if (motherIds.size() == defaultMotherSize && motherIds[1] > motherIds[0]) {
-      // keep for now, might be needed later
-      // && ((80 < std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 90) || (100 < std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 110))) {
-      // Pythia: mother range
-      for (int i = motherIds[0]; i <= motherIds[1]; i++) {
-        allMotherIds.push_back(i);
-      }
-    } else {
-      // Otherwise just use them as given
-      for (const int& id : motherIds) {
-        allMotherIds.push_back(id);
-      }
-    }
+    // RecoDecay::getOrderedMotherIds() implements the Pythia mother-id-range expansion (gated on the particle's gen status code)
+    std::vector<int64_t> allMotherIds;
+    RecoDecay::getOrderedMotherIds(mcParticle, allMotherIds, /*searchUpToQuark=*/true);
     // Loop over all mothers
-    for (const int& i : allMotherIds) {
+    for (const int64_t& i : allMotherIds) {
 
       if (i < 0 || i >= mcParticles.size()) {
         continue;
@@ -722,19 +708,12 @@ class McBuilder
       if (!current.has_mothers()) {
         break;
       }
-      auto motherIds = current.mothersIds();
-      int nextIndex = -1;
-      const int defaultMotherSize = 2;
-      if (motherIds.size() == defaultMotherSize && motherIds[1] > motherIds[0]) {
-        nextIndex = motherIds[0];
-      } else {
-        for (const int& id : motherIds) {
-          if (id >= 0 && id < mcParticles.size()) {
-            nextIndex = id;
-            break;
-          }
-        }
-      }
+      // RecoDecay::getOrderedMotherIds() implements the Pythia mother-id-range expansion (gated on
+      // the particle's gen status code). This function only ever follows a single branch upward,
+      // so just the first ordered id is used.
+      std::vector<int64_t> motherIdsOrdered;
+      RecoDecay::getOrderedMotherIds(current, motherIdsOrdered, /*searchUpToQuark=*/true);
+      int64_t nextIndex = motherIdsOrdered.empty() ? -1 : motherIdsOrdered.front();
 
       if (nextIndex < 0 || nextIndex >= mcParticles.size()) {
         break;

PWGCF/Femto/Core/modes.h

@@ -282,5 +282,22 @@ enum class KinematicVariable : o2::analysis::femto::datatypes::KinematicVariable
   kKinematicVariableLast = 4
 };
 
+/// axis title for the kinematic variable plotted by CPR / pair- and triplet-cleaner QA histograms
+constexpr const char* kinematicVariableToAxisTitle(KinematicVariable type)
+{
+  switch (type) {
+    case KinematicVariable::kKstar:
+      return "k* (GeV/#it{c})";
+    case KinematicVariable::kKt:
+      return "k_{T} (GeV/#it{c})";
+    case KinematicVariable::kMt:
+      return "m_{T} (GeV/#it{c}^{2})";
+    case KinematicVariable::kQ3:
+      return "Q_{3} (GeV/#it{c})";
+    default:
+      return "kinematic variable (GeV/#it{c})";
+  }
+}
+
 }; // namespace o2::analysis::femto::modes
 #endif // PWGCF_FEMTO_CORE_MODES_H_

PWGCF/Femto/Core/pairBuilder.h

@@ -91,9 +91,7 @@ class PairTrackTrackBuilder
             std::map<T12, std::vector<o2::framework::AxisSpec>> const& trackHistSpec2,
             std::map<T13, std::vector<o2::framework::AxisSpec>> const& pairHistSpec,
             std::map<T14, std::vector<o2::framework::AxisSpec>> const& cprHistSpec,
-            std::map<T15, std::vector<o2::framework::AxisSpec>> const& pairCleanerHistSpec,
-            int puritySlotMask1 = 1,
-            int puritySlotMask2 = 1)
+            std::map<T15, std::vector<o2::framework::AxisSpec>> const& pairCleanerHistSpec)
   {
 
     // check if correlate the same tracks or not
@@ -108,8 +106,8 @@ class PairTrackTrackBuilder
     if (mSameSpecies) {
       mTrackCleaner1.init(confTrackCleaner1);
 
-      mTrackHistManagerSe1.template init<modeSe>(registry, trackHistSpec1, confTrackSelection1, puritySlotMask1);
-      mTrackHistManagerMe1.template init<modeMe>(registry, trackHistSpec1, confTrackSelection1, puritySlotMask1);
+      mTrackHistManagerSe1.template init<modeSe>(registry, trackHistSpec1, confTrackSelection1);
+      mTrackHistManagerMe1.template init<modeMe>(registry, trackHistSpec1, confTrackSelection1);
 
       mPairHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value);
       mPairHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value);
@@ -121,10 +119,10 @@ class PairTrackTrackBuilder
     } else {
       mTrackCleaner1.init(confTrackCleaner1);
       mTrackCleaner2.init(confTrackCleaner2);
-      mTrackHistManagerSe1.template init<modeSe>(registry, trackHistSpec1, confTrackSelection1, puritySlotMask1);
-      mTrackHistManagerSe2.template init<modeSe>(registry, trackHistSpec2, confTrackSelection2, puritySlotMask2);
-      mTrackHistManagerMe1.template init<modeMe>(registry, trackHistSpec1, confTrackSelection1, puritySlotMask1);
-      mTrackHistManagerMe2.template init<modeMe>(registry, trackHistSpec2, confTrackSelection2, puritySlotMask2);
+      mTrackHistManagerSe1.template init<modeSe>(registry, trackHistSpec1, confTrackSelection1);
+      mTrackHistManagerSe2.template init<modeSe>(registry, trackHistSpec2, confTrackSelection2);
+      mTrackHistManagerMe1.template init<modeMe>(registry, trackHistSpec1, confTrackSelection1);
+      mTrackHistManagerMe2.template init<modeMe>(registry, trackHistSpec2, confTrackSelection2);
 
       mPairHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection2.pdgCodeAbs.value);
       mPairHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value);