diff --git a/PWGCF/Femto/Core/femtoUtils.h b/PWGCF/Femto/Core/femtoUtils.h index 5f9da9ab910..d7d0a61f9aa 100644 --- a/PWGCF/Femto/Core/femtoUtils.h +++ b/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) { diff --git a/PWGCF/Femto/DataModel/FemtoTables.h b/PWGCF/Femto/DataModel/FemtoTables.h index 66707be9e8b..4c890bdb4e6 100644 --- a/PWGCF/Femto/DataModel/FemtoTables.h +++ b/PWGCF/Femto/DataModel/FemtoTables.h @@ -529,7 +529,24 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmas_001, "FSIGMA", 1, femtobase::dynamic::Py, femtobase::dynamic::Pz, femtobase::dynamic::Theta); -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::dynamic::Pt, + femtobase::dynamic::P, + femtobase::dynamic::Px, + femtobase::dynamic::Py, + femtobase::dynamic::Pz, + femtobase::dynamic::Theta); +using FSigmas = FSigmas_002; DECLARE_SOA_TABLE_STAGED_VERSIONED(FSigmaMasks_001, "FSIGMAMASKS", 1, femtokinks::Mask); diff --git a/PWGCF/Femto/TableProducer/CMakeLists.txt b/PWGCF/Femto/TableProducer/CMakeLists.txt index fc9a5f82013..666b21a839c 100644 --- a/PWGCF/Femto/TableProducer/CMakeLists.txt +++ b/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) diff --git a/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx b/PWGCF/Femto/TableProducer/femtoProducerKinkPtConverter.cxx new file mode 100644 index 00000000000..e8642fe53a4 --- /dev/null +++ b/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 + +using namespace o2::analysis::femto; + +struct FemtoProducerKinkPtConverter { + + o2::framework::Produces producedSigmas; + + o2::framework::Configurable confUseRecalculatedPt{"confUseRecalculatedPt", true, "Use recalculated pT from kinematic constraints"}; + o2::framework::Configurable confFill1DHistos{"confFill1DHistos", true, "Fill 1D histograms"}; + o2::framework::Configurable 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(cfgc)}; + return workflow; +}