diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 079f217c5f3..9c533c8bf4d 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -32,6 +32,7 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include #include #include @@ -903,15 +904,20 @@ class PairHistManager float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) { - // compute pair momentum - auto sum = part1 + part2; - // Boost particle 1 to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) - // make a copy of particle 1 - auto particle1Prf = ROOT::Math::PtEtaPhiMVector(part1); - // get lorentz boost into pair rest frame - ROOT::Math::Boost boostPrf(sum.BoostToCM()); - // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k* - return static_cast(boostPrf(particle1Prf).P()); + // Use Cartesian 4-vectors: addition/M2() become pure arithmetic + const ROOT::Math::PxPyPzEVector p1(part1); + const ROOT::Math::PxPyPzEVector p2(part2); + + // Mandelstam s = (p1 + p2)^2 + const double s = (p1 + p2).M2(); + const double m1sq = p1.M2(); + const double m2sq = p2.M2(); + + // Källen function λ(s, m1^2, m2^2) = (s - m1^2 - m2^2)² - 4*m1^2*m2^2 + const double kallen = (s - m1sq - m2sq) * (s - m1sq - m2sq) - 4.0 * m1sq * m2sq; + + // k* = 0.5 * sqrt(λ/s) + return static_cast(0.5 * std::sqrt(std::max(0.0, kallen) / s)); } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; diff --git a/PWGCF/Femto/Core/tripletHistManager.h b/PWGCF/Femto/Core/tripletHistManager.h index 6f9f2323f31..8e6674f0c10 100644 --- a/PWGCF/Femto/Core/tripletHistManager.h +++ b/PWGCF/Femto/Core/tripletHistManager.h @@ -383,23 +383,24 @@ class TripletHistManager float getQ3() const { return mQ3; } private: - ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PtEtaPhiMVector const& pi, ROOT::Math::PtEtaPhiMVector const& pj) + ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PxPyPzEVector const& vi, ROOT::Math::PxPyPzEVector const& vj) { - // Convert to PxPyPzEVector to get proper Lorentz dot product - ROOT::Math::PxPyPzEVector vi(pi); - ROOT::Math::PxPyPzEVector vj(pj); - auto trackSum = vi + vj; auto trackDifference = vi - vj; - double scaling = trackDifference.Dot(trackSum) / trackSum.Dot(trackSum); + const double s = trackSum.M2(); + const double scaling = (s != 0.0) ? (vi.M2() - vj.M2()) / s : 0.0; return trackDifference - scaling * trackSum; } float getQ3(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3) { - auto q12 = getqij(part1, part2); - auto q23 = getqij(part2, part3); - auto q31 = getqij(part3, part1); + // upfront conversion to PxPyPzEVector + const ROOT::Math::PxPyPzEVector p1(part1); + const ROOT::Math::PxPyPzEVector p2(part2); + const ROOT::Math::PxPyPzEVector p3(part3); + auto q12 = getqij(p1, p2); + auto q23 = getqij(p2, p3); + auto q31 = getqij(p3, p1); double q = q12.M2() + q23.M2() + q31.M2(); return static_cast(std::sqrt(-q)); } diff --git a/PWGCF/Femto/Tasks/femtoPairEfficiency.cxx b/PWGCF/Femto/Tasks/femtoPairEfficiency.cxx index 10179a63fba..6705f731599 100644 --- a/PWGCF/Femto/Tasks/femtoPairEfficiency.cxx +++ b/PWGCF/Femto/Tasks/femtoPairEfficiency.cxx @@ -431,7 +431,7 @@ struct FemtoPairEfficiency { ROOT::Math::PtEtaPhiMVector particle1; ROOT::Math::PtEtaPhiMVector particle2; - for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsUpperIndexPolicy(tracks, tracks))) { + for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(tracks, tracks))) { bool order1 = checkTrack(p1, TrackSel1) && checkTrackPid(p1, TrackSel1) && checkTrack(p2, TrackSel2) && checkTrackPid(p2, TrackSel2);