[PWGCF] Feature: Recalculation of kink/sigma pT

PWGCF/Femto/Core/femtoUtils.h

@@ -133,6 +133,96 @@ float qn(T const& col)
   return qn;
 }
 
+/// Recalculate pT for Kinks (Sigmas) using kinematic constraints
+inline float calcPtnew(float pxMother, float pyMother, float pzMother, float pxDaughter, float pyDaughter, float pzDaughter)
+{
+  // Particle masses in GeV/c^2
+  const float massPion = 0.13957f;
+  const float massNeutron = 0.93957f;
+  const float massSigmaMinus = 1.19745f;
+
+  // Calculate mother momentum and direction versor
+  float pMother = std::sqrt(pxMother * pxMother + pyMother * pyMother + pzMother * pzMother);
+  if (pMother < 1e-6f)
+    return -999.f;
+
+  float versorX = pxMother / pMother;
+  float versorY = pyMother / pMother;
+  float versorZ = pzMother / pMother;
+
+  // Calculate daughter energy
+  float ePi = std::sqrt(massPion * massPion + pxDaughter * pxDaughter + pyDaughter * pyDaughter + pzDaughter * pzDaughter);
+
+  // Scalar product of versor with daughter momentum
+  float a = versorX * pxDaughter + versorY * pyDaughter + versorZ * pzDaughter;
+
+  // Solve quadratic equation for momentum magnitude
+  float K = massSigmaMinus * massSigmaMinus + massPion * massPion - massNeutron * massNeutron;
+  float A = 4.f * (ePi * ePi - a * a);
+  float B = -4.f * a * K;
+  float C = 4.f * ePi * ePi * massSigmaMinus * massSigmaMinus - K * K;
+
+  if (std::abs(A) < 1e-6f)
+    return -999.f;
+
+  float D = B * B - 4.f * A * C;
+  if (D < 0.f)
+    return -999.f;
+
+  float sqrtD = std::sqrt(D);
+  float P1 = (-B + sqrtD) / (2.f * A);
+  float P2 = (-B - sqrtD) / (2.f * A);
+
+  // Pick physical solution: prefer P2 if positive, otherwise P1
+  if (P2 < 0.f && P1 < 0.f)
+    return -999.f;
+  if (P2 < 0.f)
+    return P1;
+
+  // Choose solution closest to original momentum
+  float p1Diff = std::abs(P1 - pMother);
+  float p2Diff = std::abs(P2 - pMother);
+  float P = (p1Diff < p2Diff) ? P1 : P2;
+
+  // Calculate pT from recalibrated momentum
+  float pxS = versorX * P;
+  float pyS = versorY * P;
+  return std::sqrt(pxS * pxS + pyS * pyS);
+}
+
+/// Helper function to calculate recalculated pT for kink particles (Sigma/SigmaPlus)
+template <typename TParticle, typename TTrackTable>
+float getRecalculatedPtForKink(const TParticle& particle, const TTrackTable& trackTable)
+{
+  // Check if particle has chaDau index
+  if constexpr (requires { particle.has_chaDau(); }) {
+    if (particle.has_chaDau()) {
+      try {
+        auto chaDaughter = trackTable.rawIteratorAt(particle.chaDauId() - trackTable.offset());
+
+        // Extract momentum components directly from dynamic columns
+        float pxDaug = chaDaughter.px();
+        float pyDaug = chaDaughter.py();
+        float pzDaug = chaDaughter.pz();
+
+        // Get momentum components from dynamic columns
+        float pxMoth = particle.px();
+        float pyMoth = particle.py();
+        float pzMoth = particle.pz();
+
+        // Recalculate pT using kinematic constraints
+        float ptRecalc = calcPtnew(pxMoth, pyMoth, pzMoth, pxDaug, pyDaug, pzDaug);
+        if (ptRecalc > 0) {
+          return ptRecalc;
+        }
+      } catch (const std::exception& e) {
+        return -1.0f;
+      }
+    }
+  }
+  return -1.0f;
+}
+
 inline bool enableTable(const char* tableName, int userSetting, o2::framework::InitContext& initContext)
 {
   if (userSetting == 1) {

PWGCF/Femto/Core/pairHistManager.h

@@ -334,12 +334,17 @@ class PairHistManager
   }
 
   template <typename T1, typename T2>
-  void setPair(const T1& particle1, const T2& particle2)
+  void setPair(const T1& particle1, const T2& particle2, float pt1Recalc = -1.0f, float pt2Recalc = -1.0f)
   {
     // pt in track table is calculated from 1/signedPt from the original track table
     // in case of He with Z=2, we have to rescale the pt with the absolute charge
-    mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * particle1.pt(), particle1.eta(), particle1.phi(), mPdgMass1);
-    mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mPdgMass2);
+
+    // Use recalculated pT if provided (for recalculated kink pT), otherwise use particle pT
+    float pt1 = (pt1Recalc > 0.0f) ? pt1Recalc : particle1.pt();
+    float pt2 = (pt2Recalc > 0.0f) ? pt2Recalc : particle2.pt();
+
+    mParticle1 = ROOT::Math::PtEtaPhiMVector(mAbsCharge1 * pt1, particle1.eta(), particle1.phi(), mPdgMass1);
+    mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * pt2, particle2.eta(), particle2.phi(), mPdgMass2);
 
     // set kT
     mKt = getKt(mParticle1, mParticle2);
@@ -363,17 +368,17 @@ class PairHistManager
   }
 
   template <typename T1, typename T2, typename T3>
-  void setPair(const T1& particle1, const T2& particle2, const T3& col)
+  void setPair(const T1& particle1, const T2& particle2, const T3& col, float pt1Recalc = -1.0f, float pt2Recalc = -1.0f)
   {
-    setPair(particle1, particle2);
+    setPair(particle1, particle2, pt1Recalc, pt2Recalc);
     mMult = col.mult();
     mCent = col.cent();
   }
 
   template <typename T1, typename T2, typename T3, typename T4>
-  void setPair(const T1& particle1, const T2& particle2, const T3& col1, const T4& col2)
+  void setPair(const T1& particle1, const T2& particle2, const T3& col1, const T4& col2, float pt1Recalc = -1.0f, float pt2Recalc = -1.0f)
   {
-    setPair(particle1, particle2);
+    setPair(particle1, particle2, pt1Recalc, pt2Recalc);
     mMult = 0.5f * (col1.mult() + col2.mult()); // if mixing with multiplicity, should be in the same mixing bin
     mCent = 0.5f * (col1.cent() + col2.cent()); // if mixing with centrality, should be in the same mixing bin
   }

PWGCF/Femto/Core/pairProcessHelpers.h

@@ -16,6 +16,7 @@
 #ifndef PWGCF_FEMTO_CORE_PAIRPROCESSHELPERS_H_
 #define PWGCF_FEMTO_CORE_PAIRPROCESSHELPERS_H_
 
+#include "PWGCF/Femto/Core/femtoUtils.h"
 #include "PWGCF/Femto/Core/modes.h"
 #include "PWGCF/Femto/DataModel/FemtoTables.h"
 
@@ -62,16 +63,19 @@ void processSameEvent(T1 const& SliceParticle,
     if (CprManager.isClosePair()) {
       continue;
     }
+    // Calculate recalculated pT for kinks (if applicable)
+    float pt1Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p1, TrackTable);
+    float pt2Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p2, TrackTable);
     // Randomize pair order if enabled
     switch (pairOrder) {
       case kOrder12:
-        PairHistManager.setPair(p1, p2, Collision);
+        PairHistManager.setPair(p1, p2, Collision, pt1Recalc, pt2Recalc);
         break;
       case kOrder21:
-        PairHistManager.setPair(p2, p1, Collision);
+        PairHistManager.setPair(p2, p1, Collision, pt2Recalc, pt1Recalc);
         break;
       default:
-        PairHistManager.setPair(p1, p2, Collision);
+        PairHistManager.setPair(p1, p2, Collision, pt1Recalc, pt2Recalc);
     }
     // fill deta-dphi histograms with kstar cutoff
     CprManager.fill(PairHistManager.getKstar());
@@ -189,7 +193,10 @@ void processSameEvent(T1 const& SliceParticle1,
     if (CprManager.isClosePair()) {
       continue;
     }
-    PairHistManager.setPair(p1, p2, Collision);
+    // Calculate recalculated pT for kinks (if applicable)
+    float pt1Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p1, TrackTable);
+    float pt2Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p2, TrackTable);
+    PairHistManager.setPair(p1, p2, Collision, pt1Recalc, pt2Recalc);
     CprManager.fill(PairHistManager.getKstar());
     if (PairHistManager.checkPairCuts()) {
       PairHistManager.template fill<mode>();
@@ -309,7 +316,10 @@ void processMixedEvent(T1 const& Collisions,
       if (CprManager.isClosePair()) {
         continue;
       }
-      PairHistManager.setPair(p1, p2, collision1, collision2);
+      // Calculate recalculated pT for kinks (if applicable)
+      float pt1Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p1, TrackTable);
+      float pt2Recalc = o2::analysis::femto::utils::getRecalculatedPtForKink(p2, TrackTable);
+      PairHistManager.setPair(p1, p2, collision1, collision2, pt1Recalc, pt2Recalc);
       CprManager.fill(PairHistManager.getKstar());
       if (PairHistManager.checkPairCuts()) {
         PairHistManager.template fill<mode>();