[PWGCF] Fixed pT handling and minor fixes for Sigma-/+ in Femto framework (#18197)

Co-authored-by: Henrik Fribert <ge35huh@nidoqueen.ktas.ph.tum.de>
Co-authored-by: Henrik Fribert <ge35huh@nidoking.ktas.ph.tum.de>

Author: 3 people

PWGCF/Femto/Core/femtoUtils.h

@@ -167,14 +167,14 @@ concept HasEventShapeRow = requires(T row) {
 template <typename T>
 concept HasEventShape = HasEventShapeRow<std::decay_t<T>> || (requires { typename std::decay_t<T>::iterator; } && HasEventShapeRow<typename std::decay_t<T>::iterator>);
 
-/// Recalculate pT for Kinks (Sigmas) using kinematic constraints
-inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter)
+/// Recalculate the pT of a kink mother from its direction and the charged daughter momentum
+inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter,
+                       float massMother, float massChargedDaughter, float massNeutralDaughter)
 {
   float almost0 = 1e-6f;
-  // Particle masses in GeV/c^2
-  auto massPion = o2::constants::physics::MassPionCharged;
-  auto massNeutron = o2::constants::physics::MassNeutron;
-  auto massSigmaMinus = o2::constants::physics::MassSigmaMinus;
+  const float massPion = massChargedDaughter;
+  const float massNeutron = massNeutralDaughter;
+  const float massSigmaMinus = massMother;
 
   // Calculate mother momentum and direction versor
   float pMother = std::sqrt(pxMother * pxMother + pyMother * pyMother + pzMother * pzMother);

PWGCF/Femto/Core/kinkBuilder.h

@@ -27,6 +27,7 @@
 #include "Common/Core/RecoDecay.h"
 
 #include <CommonConstants/MathConstants.h>
+#include <CommonConstants/PhysicsConstants.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/Configurable.h>
 #include <Framework/HistogramRegistry.h>
@@ -60,34 +61,39 @@ struct ConfKinkFilters : o2::framework::ConfigurableGroup {
 
 // selections bits for all kinks
 // NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
-#define KINK_DEFAULT_BITS                                                                                                                            \
-  o2::framework::Configurable<bool> passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."}; \
-  o2::framework::Configurable<std::vector<float>> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"};                          \
-  o2::framework::Configurable<std::vector<float>> transRadMin{"transRadMin", {20.f}, "Minimum transverse radius (cm)"};                              \
-  o2::framework::Configurable<std::vector<float>> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"};                             \
-  o2::framework::Configurable<std::vector<float>> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"};        \
-  o2::framework::Configurable<std::vector<float>> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"};            \
-  o2::framework::Configurable<std::vector<float>> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"};            \
-  o2::framework::Configurable<std::vector<float>> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"};                        \
-  o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"};                         \
-  o2::framework::Configurable<std::vector<float>> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"};                                 \
-  o2::framework::Configurable<std::vector<float>> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"};                                  \
-  o2::framework::Configurable<std::vector<float>> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"};
+#define KINK_DEFAULT_BITS                                                                                                                                 \
+  o2::framework::Configurable<bool> passThrough{"passThrough", false, "If true, all Kinks are passed through. Bits for all selections are stored."};      \
+  o2::framework::Configurable<std::vector<float>> kinkTopoDcaMax{"kinkTopoDcaMax", {2.0f}, "Maximum kink topological DCA"};                               \
+  o2::framework::Configurable<std::vector<float>> transRadMin{"transRadMin", {19.6f}, "Minimum transverse radius (cm)"};                                  \
+  o2::framework::Configurable<std::vector<float>> transRadMax{"transRadMax", {100.f}, "Maximum transverse radius (cm)"};                                  \
+  o2::framework::Configurable<std::vector<float>> dauAbsEtaMax{"dauAbsEtaMax", {1.0f}, "Maximum absolute pseudorapidity for daughter track"};             \
+  o2::framework::Configurable<std::vector<float>> dauDcaPvMin{"dauDcaPvMin", {0.1f}, "Minimum DCA of daughter from primary vertex (cm)"};                 \
+  o2::framework::Configurable<std::vector<float>> mothDcaPvMax{"mothDcaPvMax", {1.0f}, "Maximum DCA of mother from primary vertex (cm)"};                 \
+  o2::framework::Configurable<std::vector<float>> qtAPMin{"qtAPMin", {0.15f}, "Minimum qT_AP for Sigma candidates"};                                      \
+  o2::framework::Configurable<std::vector<float>> qtAPMax{"qtAPMax", {0.2f}, "Maximum qT_AP for Sigma candidates"};                                       \
+  o2::framework::Configurable<std::vector<float>> cosPointingAngleMin{"cosPointingAngleMin", {0.0f}, "Minimum cosine of pointing angle"};                 \
+  o2::framework::Configurable<std::vector<float>> ptOriginalMin{"ptOriginalMin", {1.2f}, "Minimum original (not recalculated) pT of the mother (GeV/c)"}; \
+  o2::framework::Configurable<std::vector<float>> ptOriginalMax{"ptOriginalMax", {10.f}, "Maximum original (not recalculated) pT of the mother (GeV/c)"};
 
 // derived selection bits for sigma
 struct ConfSigmaBits : o2::framework::ConfigurableGroup {
   std::string prefix = std::string("SigmaBits");
   KINK_DEFAULT_BITS
+  o2::framework::Configurable<std::vector<float>> alphaAPMin{"alphaAPMin", {-1.0f}, "Minimum Alpha_AP for Sigma candidates"};
+  o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {0.0f}, "Maximum Alpha_AP for Sigma candidates"};
   o2::framework::Configurable<std::vector<float>> chaDauTpcPion{"chaDauTpcPion", {5.f}, "Maximum |nsigma_Pion| TPC for charged daughter tracks"};
 };
 
 // derived selection bits for sigma plus
 struct ConfSigmaPlusBits : o2::framework::ConfigurableGroup {
   std::string prefix = std::string("SigmaPlusBits");
   KINK_DEFAULT_BITS
+  o2::framework::Configurable<std::vector<float>> alphaAPMin{"alphaAPMin", {0.0f}, "Minimum Alpha_AP for SigmaPlus candidates"};
+  o2::framework::Configurable<std::vector<float>> alphaAPMax{"alphaAPMax", {1.0f}, "Maximum Alpha_AP for SigmaPlus candidates"};
   o2::framework::Configurable<std::vector<float>> chaDauTpcProton{"chaDauTpcProton", {5.f}, "Maximum |nsigma_Proton| TPC for charged daughter tracks"};
-  o2::framework::Configurable<std::vector<float>> chaDauTofProton{"chaDauTofProton", {5.f}, "Maximum combined |nsigma_Proton| (TPC+TOF) for charged daughter tracks"};
-  o2::framework::Configurable<float> pidThres{"pidThres", 0.75f, "Momentum threshold for using TOF/combined pid for daughter tracks (GeV/c)"};
+  o2::framework::Configurable<std::vector<float>> chaDauTofProton{"chaDauTofProton", {}, "Maximum |nsigma_Proton| TOF for charged daughter tracks"};
+  o2::framework::Configurable<bool> requireTof{"requireTof", false, "If true, TOF PID is a minimal selection. If false, TOF PID is optional"};
+  o2::framework::Configurable<bool> keepTracksWithoutTof{"keepTracksWithoutTof", true, "If true, the bit mask for the TOF selection will be true for all limits if the daughter track has no TOF"};
 };
 
 #undef KINK_DEFAULT_BITS
@@ -153,6 +159,9 @@ enum KinkSeles {
   kQtAPMax,
   kCosPointingAngleMin,
 
+  kPtOriginalMin,
+  kPtOriginalMax,
+
   kKinkSelsMax
 };
 
@@ -173,7 +182,9 @@ const std::unordered_map<KinkSeles, std::string> kinkSelectionNames = {
   {kAlphaAPMax, "alphaAPMax"},
   {kQtAPMin, "qtAPMin"},
   {kQtAPMax, "qtAPMax"},
-  {kCosPointingAngleMin, "cosPointingAngleMin"}};
+  {kCosPointingAngleMin, "cosPointingAngleMin"},
+  {kPtOriginalMin, "ptOriginalMin"},
+  {kPtOriginalMax, "ptOriginalMax"}};
 
 /// enum for all kink pre-filters (evaluated in checkFilters, before the selection bitmask)
 enum KinkFilters {
@@ -231,9 +242,9 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
     if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) {
       mMassSigmaPlusLowerLimit = filter.massMinSigmaPlus.value;
       mMassSigmaPlusUpperLimit = filter.massMaxSigmaPlus.value;
-      mPidThreshold = config.pidThres.value;
-      this->addSelection(kChaDaughTpcProton, kinkSelectionNames.at(kChaDaughTpcProton), config.chaDauTpcProton.value, limits::kAbsUpperLimit, false, false, true);
-      this->addSelection(kChaDaughTofProton, kinkSelectionNames.at(kChaDaughTofProton), config.chaDauTofProton.value, limits::kUpperLimit, false, false, true);
+      mKeepTracksWithoutTof = config.keepTracksWithoutTof.value;
+      this->addSelection(kChaDaughTpcProton, kinkSelectionNames.at(kChaDaughTpcProton), config.chaDauTpcProton.value, limits::kAbsUpperLimit, true, true, false);
+      this->addSelection(kChaDaughTofProton, kinkSelectionNames.at(kChaDaughTofProton), config.chaDauTofProton.value, limits::kAbsUpperLimit, true, config.requireTof.value, false);
     }
 
     this->addSelection(kKinkTopoDcaMax, kinkSelectionNames.at(kKinkTopoDcaMax), config.kinkTopoDcaMax.value, limits::kUpperLimit, true, true, false);
@@ -247,6 +258,8 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
     this->addSelection(kQtAPMin, kinkSelectionNames.at(kQtAPMin), config.qtAPMin.value, limits::kLowerLimit, true, true, false);
     this->addSelection(kQtAPMax, kinkSelectionNames.at(kQtAPMax), config.qtAPMax.value, limits::kUpperLimit, true, true, false);
     this->addSelection(kCosPointingAngleMin, kinkSelectionNames.at(kCosPointingAngleMin), config.cosPointingAngleMin.value, limits::kLowerLimit, true, true, false);
+    this->addSelection(kPtOriginalMin, kinkSelectionNames.at(kPtOriginalMin), config.ptOriginalMin.value, limits::kLowerLimit, true, true, false);
+    this->addSelection(kPtOriginalMax, kinkSelectionNames.at(kPtOriginalMax), config.ptOriginalMax.value, limits::kUpperLimit, true, true, false);
 
     this->setupSelectionHistogram<SelectionHistName>(registry);
 
@@ -281,7 +294,7 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
   };
 
   template <typename T1, typename T2, typename T3>
-  void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& col)
+  void computeKinkKinematics(T1 const& kinkCand, T2 const& /*tracks*/, T3 const& /*col*/)
   {
     std::array<float, 3> momMother = {kinkCand.pxMoth(), kinkCand.pyMoth(), kinkCand.pzMoth()};
     float kinkMomP = RecoDecay::p(momMother);
@@ -300,12 +313,11 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
     float p2A = kinkDauP * kinkDauP;
     mQtAp = (p2V0 > 0.f) ? std::sqrt(std::max(0.f, p2A - dp * dp / p2V0)) : 0.f;
 
-    std::array<float, 3> vMother = {kinkCand.xDecVtx() - col.posX(), kinkCand.yDecVtx() - col.posY(), kinkCand.zDecVtx() - col.posZ()};
+    std::array<float, 3> vMother = {kinkCand.xDecVtx(), kinkCand.yDecVtx(), kinkCand.zDecVtx()};
     float vMotherNorm = std::sqrt(std::inner_product(vMother.begin(), vMother.end(), vMother.begin(), 0.f));
     mCosPointingAngle = (vMotherNorm > 0.f && kinkMomP > 0.f) ? (std::inner_product(momMother.begin(), momMother.end(), vMother.begin(), 0.f)) / (kinkMomP * vMotherNorm) : 0.f;
     mTransRadius = std::hypot(kinkCand.xDecVtx(), kinkCand.yDecVtx());
 
-    mKinkDauP = kinkDauP;
     mKinkDauEta = RecoDecay::eta(momDaughter);
 
     mKinkAngle = 0.f;
@@ -327,6 +339,9 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
     this->evaluateObservable(kQtAPMin, mQtAp);
     this->evaluateObservable(kQtAPMax, mQtAp);
     this->evaluateObservable(kCosPointingAngleMin, mCosPointingAngle);
+    // the stored pT is the recalculated one, the original pT is only available as selection bits
+    this->evaluateObservable(kPtOriginalMin, mKinkMotherPtOriginal);
+    this->evaluateObservable(kPtOriginalMax, mKinkMotherPtOriginal);
     this->evaluateObservable(kKinkTopoDcaMax, kinkCand.dcaKinkTopo());
 
     // Compute transRadius
@@ -345,14 +360,11 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
       this->evaluateObservable(kChaDaughTpcPion, chaDaughter.tpcNSigmaPi());
     }
     if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) {
-      if (mKinkDauP < mPidThreshold) {
-        this->evaluateObservable(kChaDaughTpcProton, chaDaughter.tpcNSigmaPr());
+      this->evaluateObservable(kChaDaughTpcProton, chaDaughter.tpcNSigmaPr());
+      if (chaDaughter.hasTOF()) {
+        this->evaluateObservable(kChaDaughTofProton, chaDaughter.tofNSigmaPr());
       } else {
-        if (chaDaughter.hasTOF()) {
-          this->evaluateObservable(kChaDaughTofProton, std::abs(chaDaughter.tofNSigmaPr()));
-        } else {
-          this->evaluateObservable(kChaDaughTofProton, 999.f);
-        }
+        this->evaluateObservable(kChaDaughTofProton, mKeepTracksWithoutTof ? 0.f : 999.f);
       }
     }
 
@@ -368,9 +380,20 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
     mKinkMotherEta = RecoDecay::eta(momMother);
     mKinkMotherPhi = RecoDecay::phi(momMother);
 
-    // Recalculate pT using kinematic constraints
-    float ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2]);
-    mKinkMotherPt = (ptRecalc > 0.f) ? ptRecalc : std::hypot(momMother[0], momMother[1]);
+    // Recalculate pT using kinematic constraints of the decay channel
+    float ptRecalc = -999.f;
+    if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) {
+      // Sigma- -> pi- n
+      ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2],
+                                  o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron);
+    }
+    if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) {
+      // Sigma+ -> p pi0
+      ptRecalc = utils::calcPtnew(momMother[0], momMother[1], momMother[2], momDaughter[0], momDaughter[1], momDaughter[2],
+                                  o2::constants::physics::MassSigmaPlus, o2::constants::physics::MassProton, o2::constants::physics::MassPionNeutral);
+    }
+    mKinkMotherPtOriginal = std::hypot(momMother[0], momMother[1]);
+    mKinkMotherPt = (ptRecalc > 0.f) ? ptRecalc : mKinkMotherPtOriginal;
   }
 
   template <typename T>
@@ -443,7 +466,7 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
   float mMassSigmaUpperLimit = 1.25f;
   float mMassSigmaPlusLowerLimit = 1.15f;
   float mMassSigmaPlusUpperLimit = 1.25f;
-  float mPidThreshold = 0.75f;
+  bool mKeepTracksWithoutTof = true;
 
   // kinematic filters
   float mPtMin = 0.f;
@@ -454,7 +477,8 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
   float mPhiMax = o2::constants::math::TwoPI;
 
   // mother kinematic
-  float mKinkMotherPt = 0.f;
+  float mKinkMotherPt = 0.f;         // recalculated pT (falls back to the original pT if the recalculation fails)
+  float mKinkMotherPtOriginal = 0.f; // pT of the mother track as reconstructed
   float mKinkMotherEta = 0.f;
   float mKinkMotherPhi = 0.f;
 
@@ -464,7 +488,6 @@ class KinkSelection : public baseselection::BaseSelection<float, o2::analysis::f
   float mCosPointingAngle = 0.f;
   float mTransRadius = 0.f;
   float mKinkDauEta = 0.f;
-  float mKinkDauP = 0.f;
   float mKinkAngle = 0.f;
 };
 
@@ -549,7 +572,8 @@ class KinkBuilder
       }
     }
 
-    if (mProduceSigmas || mProduceSigmaMasks || mProduceSigmaExtras || mProduceSigmaPlus || mProduceSigmaPlusMasks || mProduceSigmaPlusExtras) {
+    if (mProduceSigmas || mProduceLiteSigmas || mProduceSigmaMasks || mProduceSigmaExtras ||
+        mProduceSigmaPlus || mProduceLiteSigmaPlus || mProduceSigmaPlusMasks || mProduceSigmaPlusExtras) {
       mFillAnyTable = true;
     } else {
       LOG(info) << "No tables configured, Selection object will not be configured...";

PWGCF/Femto/Core/kinkHistManager.h

@@ -284,25 +284,27 @@ class KinkHistManager
             std::map<trackhistmanager::TrackHist, std::vector<o2::framework::AxisSpec>> const& ChaDauSpecs)
   {
     mHistogramRegistry = registry;
-    mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value;
+    mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value);
 
     int chaDauPdgCodeAbs = 0;
     int chaDauCharge = 0;
     const int absCharge = 1;
 
-    if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) {
+    if (mPdgCode == PDG_t::kSigmaMinus) {
       if (ConfKinkSelection.sign.value < 0) {
         chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus);
         chaDauCharge = -1;
       } else {
+        mPdgCode = -1 * mPdgCode; // anti-Sigma- is positively charged and has negative pdg code
         chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus);
         chaDauCharge = 1;
       }
-    } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) {
+    } else if (mPdgCode == PDG_t::kSigmaPlus) {
       if (ConfKinkSelection.sign.value > 0) {
         chaDauPdgCodeAbs = std::abs(PDG_t::kProton);
         chaDauCharge = 1;
       } else {
+        mPdgCode = -1 * mPdgCode; // anti-Sigma+ is negatively charged and has negative pdg code
         chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar);
         chaDauCharge = -1;
       }
@@ -332,26 +334,28 @@ class KinkHistManager
             T3 const& ConfChaDauBinningQa)
   {
     mHistogramRegistry = registry;
-    mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value;
+    mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value);
     this->enableOptionalHistograms(ConfKinkBinningQa);
 
     int chaDauPdgCodeAbs = 0;
     int chaDauCharge = 0;
     const int absCharge = 1;
 
-    if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) {
+    if (mPdgCode == PDG_t::kSigmaMinus) {
       if (ConfKinkSelection.sign.value < 0) {
         chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus);
         chaDauCharge = -1;
       } else {
+        mPdgCode = -1 * mPdgCode; // anti-Sigma- is positively charged and has negative pdg code
         chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus);
         chaDauCharge = 1;
       }
-    } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) {
+    } else if (mPdgCode == PDG_t::kSigmaPlus) {
       if (ConfKinkSelection.sign.value > 0) {
         chaDauPdgCodeAbs = std::abs(PDG_t::kProton);
         chaDauCharge = 1;
       } else {
+        mPdgCode = -1 * mPdgCode; // anti-Sigma+ is negatively charged and has negative pdg code
         chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar);
         chaDauCharge = -1;
       }

PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx

@@ -16,6 +16,7 @@
 #include "PWGCF/Femto/Core/femtoUtils.h"
 #include "PWGCF/Femto/DataModel/FemtoTables.h"
 
+#include <CommonConstants/PhysicsConstants.h>
 #include <Framework/AnalysisDataModel.h>
 #include <Framework/AnalysisHelpers.h>
 #include <Framework/AnalysisTask.h>
@@ -73,7 +74,9 @@ struct FemtoProducerKinkPtConverter {
         float pyMoth = sigma.pt() * std::sin(sigma.phi());
         float pzMoth = sigma.pt() * std::sinh(sigma.eta());
 
-        float ptRecalc = utils::calcPtnew(pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug);
+        // Sigma- -> pi- n
+        float ptRecalc = utils::calcPtnew(pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug,
+                                          o2::constants::physics::MassSigmaMinus, o2::constants::physics::MassPionCharged, o2::constants::physics::MassNeutron);
 
         ROOT::Math::PtEtaPhiMVector recalcVec(ptRecalc, sigma.eta(), sigma.phi(), sigma.mass());
         float ptFrom4Vec = recalcVec.Pt();