[PWGCF] Feature: Recalculation of kink/sigma pT

PWGCF/Femto/Core/femtoUtils.h

@@ -133,6 +133,63 @@ 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);
+}
+
 inline bool enableTable(const char* tableName, int userSetting, o2::framework::InitContext& initContext)
 {
   if (userSetting == 1) {

PWGCF/Femto/DataModel/FemtoTables.h

@@ -529,7 +529,24 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmas_001, "FSIGMA", 1,
                                    femtobase::dynamic::Py<femtobase::stored::SignedPt, femtobase::stored::Phi>,
                                    femtobase::dynamic::Pz<femtobase::stored::SignedPt, femtobase::stored::Eta>,
                                    femtobase::dynamic::Theta<femtobase::stored::Eta>);
-using FSigmas = FSigmas_001;
+
+// table for basic sigma information
+DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmas_002, "FSIGMA", 2,
+                                   o2::soa::Index<>,
+                                   femtobase::stored::FColId, // use sign to differentiate between sigma minus (-1) and anti sigma minus (+1)
+                                   femtobase::stored::SignedPt,
+                                   femtobase::stored::Eta,
+                                   femtobase::stored::Phi,
+                                   femtobase::stored::Mass,
+                                   femtokinks::ChaDauId,
+                                   femtobase::dynamic::Sign<femtobase::stored::SignedPt>,
+                                   femtobase::dynamic::Pt<femtobase::stored::SignedPt>,
+                                   femtobase::dynamic::P<femtobase::stored::SignedPt, femtobase::stored::Eta>,
+                                   femtobase::dynamic::Px<femtobase::stored::SignedPt, femtobase::stored::Phi>,
+                                   femtobase::dynamic::Py<femtobase::stored::SignedPt, femtobase::stored::Phi>,
+                                   femtobase::dynamic::Pz<femtobase::stored::SignedPt, femtobase::stored::Eta>,
+                                   femtobase::dynamic::Theta<femtobase::stored::Eta>);
+using FSigmas = FSigmas_002;
 
 DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmaMasks_001, "FSIGMAMASKS", 1,
                                    femtokinks::Mask);

PWGCF/Femto/TableProducer/CMakeLists.txt

@@ -18,3 +18,8 @@ o2physics_add_dpl_workflow(femto-producer-derived-to-derived
           SOURCES ./femtoProducerDerivedToDerived.cxx
           PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
           COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(femto-producer-kink-pt-converter
+          SOURCES ./femtoProducerKinkPtConverter.cxx
+          PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore
+          COMPONENT_NAME Analysis)

PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx

@@ -0,0 +1,104 @@
+// Copyright 2019-2025 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file femtoProducerKinkPtConverter.cxx
+/// \brief Task that converts FSigmas_001 to FSigmas_002 with recalculated pT
+/// \author Henrik Fribert, TU München, henrik.fribert@tum.de
+
+#include "PWGCF/Femto/Core/femtoUtils.h"
+#include "PWGCF/Femto/DataModel/FemtoTables.h"
+
+#include "Framework/AnalysisDataModel.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/runDataProcessing.h"
+
+#include <Math/Vector4D.h>
+
+using namespace o2::analysis::femto;
+
+struct FemtoProducerKinkPtConverter {
+
+  o2::framework::Produces<o2::aod::FSigmas> producedSigmas;
+
+  o2::framework::Configurable<bool> confUseRecalculatedPt{"confUseRecalculatedPt", true, "Use recalculated pT from kinematic constraints"};
+  o2::framework::Configurable<bool> confFill1DHistos{"confFill1DHistos", true, "Fill 1D histograms"};
+  o2::framework::Configurable<bool> confFill2DHistos{"confFill2DHistos", true, "Fill 2D histograms"};
+
+  o2::framework::HistogramRegistry mHistogramRegistry{"FemtoSigmaConverter", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject};
+
+  void init(o2::framework::InitContext&)
+  {
+    if (confFill1DHistos) {
+      mHistogramRegistry.add("hPtOriginal", "Original pT;p_{T} (GeV/c);Counts", o2::framework::kTH1F, {{100, 0, 10}});
+      mHistogramRegistry.add("hPtRecalculated", "Recalculated pT;p_{T} (GeV/c);Counts", o2::framework::kTH1F, {{100, 0, 10}});
+      mHistogramRegistry.add("hRecalcSuccess", "Recalculation Success;Success (0=fail, 1=success);Counts", o2::framework::kTH1I, {{2, -0.5, 1.5}});
+    }
+    if (confFill2DHistos) {
+      mHistogramRegistry.add("hPtOriginalVsRecalculated", "Original vs Recalculated pT;p_{T,orig} (GeV/c);p_{T,recalc} (GeV/c)", o2::framework::kTH2F, {{100, 0, 10}, {100, 0, 10}});
+    }
+  }
+
+  void process(o2::aod::FSigmas_001 const& sigmasV1,
+               o2::aod::FTracks const& tracks)
+  {
+    for (const auto& sigma : sigmasV1) {
+
+      float signedPtToUse = sigma.signedPt();
+
+      if (confUseRecalculatedPt) {
+        auto chaDaughter = tracks.rawIteratorAt(sigma.chaDauId() - tracks.offset());
+
+        float pxDaug = chaDaughter.pt() * std::cos(chaDaughter.phi());
+        float pyDaug = chaDaughter.pt() * std::sin(chaDaughter.phi());
+        float pzDaug = chaDaughter.pt() * std::sinh(chaDaughter.eta());
+
+        float pxMoth = sigma.pt() * std::cos(sigma.phi());
+        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);
+
+        ROOT::Math::PtEtaPhiMVector recalcVec(ptRecalc, sigma.eta(), sigma.phi(), sigma.mass());
+        float ptFrom4Vec = recalcVec.Pt();
+
+        if (ptFrom4Vec > 0) {
+          signedPtToUse = ptFrom4Vec * utils::signum(sigma.signedPt());
+
+          if (confFill1DHistos) {
+            mHistogramRegistry.fill(HIST("hPtOriginal"), sigma.pt());
+            mHistogramRegistry.fill(HIST("hPtRecalculated"), ptFrom4Vec);
+            mHistogramRegistry.fill(HIST("hRecalcSuccess"), 1);
+          }
+          if (confFill2DHistos) {
+            mHistogramRegistry.fill(HIST("hPtOriginalVsRecalculated"), sigma.pt(), ptFrom4Vec);
+          }
+        } else {
+          if (confFill1DHistos) {
+            mHistogramRegistry.fill(HIST("hRecalcSuccess"), 0);
+          }
+        }
+      }
+
+      producedSigmas(sigma.fColId(),
+                     signedPtToUse,
+                     sigma.eta(),
+                     sigma.phi(),
+                     sigma.mass(),
+                     sigma.chaDauId());
+    }
+  }
+};
+
+o2::framework::WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc)
+{
+  o2::framework::WorkflowSpec workflow{adaptAnalysisTask<FemtoProducerKinkPtConverter>(cfgc)};
+  return workflow;
+}