From 4649a7a3c2afc73fc888e4935b419d831862348e Mon Sep 17 00:00:00 2001 From: Anton Riedel Date: Sun, 1 Nov 2026 22:29:43 +0100 Subject: [PATCH 1/3] Feat: clean up math functions --- PWGCF/Femto/Core/femtoMath.h | 283 ++++++++++++++++++++++++++ PWGCF/Femto/Core/pairHistManager.h | 256 +++-------------------- PWGCF/Femto/Core/tripletHistManager.h | 180 ++++++++-------- 3 files changed, 391 insertions(+), 328 deletions(-) create mode 100644 PWGCF/Femto/Core/femtoMath.h diff --git a/PWGCF/Femto/Core/femtoMath.h b/PWGCF/Femto/Core/femtoMath.h new file mode 100644 index 00000000000..ead624399bc --- /dev/null +++ b/PWGCF/Femto/Core/femtoMath.h @@ -0,0 +1,283 @@ +// 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 femtoMath.h +/// \brief shared kinematic math (kstar, Q3, mt, kt, spherical-harmonics kinematics, ...) +/// used by the pair and triplet histogram managers +/// \author Anton Riedel, TU München, anton.riedel@tum.de + +#ifndef PWGCF_FEMTO_CORE_FEMTOMATH_H_ +#define PWGCF_FEMTO_CORE_FEMTOMATH_H_ + +#include "PWGCF/Femto/Core/modes.h" + +#include "Common/Core/RecoDecay.h" + +#include +#include + +#include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) +#include +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include + +#include +#include +#include +#include + +namespace o2::analysis::femto::femtomath +{ + +// SH reference frame selectors for computeShKinematics(), mirrors ConfPairBinning.shFrame +constexpr int kShFrameLcmsNonIdentical = 0; +constexpr int kShFrameLcmsIdentical = 1; +constexpr int kShFramePrf = 2; + +inline float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + auto sum = (part1 + part2); + double kt = 0.5 * sum.Pt(); + return static_cast(kt); +} + +inline float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, modes::TransverseMassType mtType = modes::TransverseMassType::kAveragePdgMass) +{ + auto sum = part1 + part2; + double mt = 0; + double mass = 0; + switch (mtType) { + case modes::TransverseMassType::kAveragePdgMass: + mass = 0.5 * (part1.M() + part2.M()); + mt = std::hypot(0.5 * sum.Pt(), mass); + break; + case modes::TransverseMassType::kReducedPdgMass: + mass = 2. * (part1.M() * part2.M()) / (part1.M() + part2.M()); + mt = std::hypot(0.5 * sum.Pt(), mass); + break; + case modes::TransverseMassType::kMt4Vector: + mt = 0.5 * sum.Mt(); + break; + default: + LOG(fatal) << "Invalid transverse mass type, breaking..."; + } + return static_cast(mt); +} + +inline float getMinv(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + return static_cast((part1 + part2).M()); +} + +inline float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + // 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)); +} + +// Alpha: azimuthal opening angle between the two particles' transverse momenta, +// i.e. phi2 - phi1 wrapped to [-pi, pi). Purely a lab-frame angular variable (no +// boosting involved) -- small |alpha| means the pair sits close together in phi, +// which is where detector effects like track splitting/merging and two-track +// resolution losses tend to show up. +inline float getAlpha(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + return static_cast(RecoDecay::constrainAngle(part2.phi() - part1.phi(), -o2::constants::math::PI)); +} + +// BetaLab: lab-frame angle between the pair's relative 3-momentum (p1 - p2) and its +// total 3-momentum (p1 + p2, i.e. the kT direction). Note the (p1 - p2) here is the +// raw 3-vector difference, not the Lorentz-invariant k* -- this is a lab-frame proxy +// for how the pair's separation axis is oriented relative to its direction of motion, +// without needing to boost to the pair rest frame. +inline float getBetaLab(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + ROOT::Math::XYZVector p1(part1.Px(), part1.Py(), part1.Pz()); + ROOT::Math::XYZVector p2(part2.Px(), part2.Py(), part2.Pz()); + + ROOT::Math::XYZVector kstar = p1 - p2; + ROOT::Math::XYZVector kt = p1 + p2; + + return static_cast(std::acos(kstar.Dot(kt) / std::sqrt(kstar.Mag2() * kt.Mag2()))); +} + +inline std::tuple computeBertschPrattLCMS(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) +{ + const ROOT::Math::PxPyPzEVector p1(part1); + const ROOT::Math::PxPyPzEVector p2(part2); + const ROOT::Math::PxPyPzEVector pSum = p1 + p2; + const ROOT::Math::PxPyPzEVector pDiff = p1 - p2; + + const double tPx = pSum.Px(); + const double tPy = pSum.Py(); + const double tPz = pSum.Pz(); + const double tE = pSum.E(); + + const double tPt = std::sqrt(tPx * tPx + tPy * tPy); + const double tMt = std::sqrt(tE * tE - tPz * tPz); + + static constexpr double MinTransverseMomentum = 1e-9; + if (tPt < MinTransverseMomentum || tMt < MinTransverseMomentum) { + return {0.0f, 0.0f, 0.0f}; + } + + const double betaL = tPz / tE; + const double gammaL = tE / tMt; + + // out/side/long are linear projections of the momentum, so projecting the + // difference p1-p2 directly is equivalent to (and half the divisions of) + // computing each particle's out/side/long component and subtracting. + const double qOut = (pDiff.Px() * tPx + pDiff.Py() * tPy) / tPt; + const double qSide = (pDiff.Py() * tPx - pDiff.Px() * tPy) / tPt; + const double qLong = gammaL * (pDiff.Pz() - betaL * pDiff.E()); + + return {static_cast(qOut), static_cast(qSide), static_cast(qLong)}; +} + +// Return the bin index for value given ascending bin edges, or -1 if out of range. +// edges = {e0, e1, ..., eN} defines N bins [e0,e1), [e1,e2), ..., [e_{N-1},eN). +inline int findShBin(double value, std::vector const& edges) +{ + static constexpr std::size_t MinEdgesForOneBin = 2; + if (edges.size() < MinEdgesForOneBin || value < edges.front() || value >= edges.back()) { + return -1; + } + for (std::size_t i = 0; i < edges.size() - 1; ++i) { + if (value >= edges[i] && value < edges[i + 1]) { + return static_cast(i); + } + } + return -1; +} + +// Kinematics feeding the spherical-harmonics decomposition, ported 1:1 from +// FemtoUniverseMath::newpairfunc so results are numerically comparable. +// Returns {kv, out, side, long}, selected by shFrame: +// kShFrameLcmsNonIdentical: {kstar, fDKOut, fDKSide, fDKLong} (LCMS components of particle 1) +// kShFrameLcmsIdentical: {qinv, outLCMS, sideLCMS, longLCMS} (LCMS q-differences) +// kShFramePrf: {|qPRF|, outPRF, sidePRF, longPRF} (matches FemtoUniverse isIdenPRF=true) +inline std::tuple computeShKinematics(ROOT::Math::PtEtaPhiMVector const& part1, + ROOT::Math::PtEtaPhiMVector const& part2, + int shFrame) +{ + const ROOT::Math::PxPyPzEVector p1(part1); + const ROOT::Math::PxPyPzEVector p2(part2); + const ROOT::Math::PxPyPzEVector sum = p1 + p2; + + const double tPx = sum.Px(); + const double tPy = sum.Py(); + const double tPz = sum.Pz(); + const double tE = sum.E(); + const double tPtSq = tPx * tPx + tPy * tPy; + const double tMtSq = tE * tE - tPz * tPz; + + static constexpr double MinTransverseMomentum = 1e-9; + if (tPtSq < MinTransverseMomentum || tMtSq < MinTransverseMomentum) { + return {0.f, 0.f, 0.f, 0.f}; + } + const double tPt = std::sqrt(tPtSq); + const double tMt = std::sqrt(tMtSq); + const double tM = std::sqrt(std::max(0.0, tMtSq - tPtSq)); + + const double beta = tPz / tE; + const double gamma = tE / tMt; + + // LCMS components of particle 1 (fDKOut/Side/Long in newpairfunc) + const double e1 = p1.E(); + const double fDKOut = (p1.Px() * tPx + p1.Py() * tPy) / tPt; + const double fDKSide = (-p1.Px() * tPy + p1.Py() * tPx) / tPt; + const double fDKLong = gamma * (p1.Pz() - beta * e1); + const double fDE = gamma * (e1 - beta * p1.Pz()); + + // LCMS components of particle 2 + const double e2 = p2.E(); + const double px2 = (p2.Px() * tPx + p2.Py() * tPy) / tPt; + const double py2 = (p2.Py() * tPx - p2.Px() * tPy) / tPt; + const double pz2 = gamma * (p2.Pz() - beta * e2); + + const double outLCMS = fDKOut - px2; + const double sideLCMS = fDKSide - py2; + const double longLCMS = fDKLong - pz2; + + if (shFrame == kShFrameLcmsIdentical) { + // LCMS identical: Ylm <- LCMS q-differences, axis <- qinv + const double qinv = std::sqrt(outLCMS * outLCMS + sideLCMS * sideLCMS + longLCMS * longLCMS); + return {static_cast(qinv), + static_cast(outLCMS), static_cast(sideLCMS), static_cast(longLCMS)}; + } + if (shFrame == kShFramePrf) { + // PRF: Ylm <- PRF q-differences, axis <- |q_PRF| (matches FemtoUniverse isIdenPRF=true) + const double pE2LCMS = gamma * (e2 - beta * p2.Pz()); + const double betaOut = tPt / tMt; + const double gammaOut = tMt / tM; + const double outPRF = gammaOut * (outLCMS - betaOut * (fDE - pE2LCMS)); + const double sidePRF = sideLCMS; + const double longPRF = longLCMS; + const double qPRF = std::sqrt(outPRF * outPRF + sidePRF * sidePRF + longPRF * longPRF); + return {static_cast(qPRF), + static_cast(outPRF), static_cast(sidePRF), static_cast(longPRF)}; + } + // LCMS non-identical (default): Ylm <- LCMS components of particle 1, axis <- k* + const double betaOut = tPt / tMt; + const double gammaOut = tMt / tM; + const double fKOut = gammaOut * (fDKOut - betaOut * fDE); + const double kstar = std::sqrt(fKOut * fKOut + fDKSide * fDKSide + fDKLong * fDKLong); + return {static_cast(kstar), + static_cast(fDKOut), static_cast(fDKSide), static_cast(fDKLong)}; +} + +template +inline float getPairPhiFromEventPlaneAngle(const ROOT::Math::PtEtaPhiMVector& part1, const ROOT::Math::PtEtaPhiMVector& part2, const T1& col) +{ + const ROOT::Math::PtEtaPhiMVector trackSum = part1 + part2; + float phi = RecoDecay::constrainAngle(trackSum.Phi() - col.eventPlaneAngle()); + return phi; +} + +// triplet-only helpers + +inline ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PxPyPzEVector const& vi, ROOT::Math::PxPyPzEVector const& vj) +{ + auto trackSum = vi + vj; + auto trackDifference = vi - vj; + const double s = trackSum.M2(); + const double scaling = (s != 0.0) ? (vi.M2() - vj.M2()) / s : 0.0; + return trackDifference - scaling * trackSum; +} + +inline float getQ3(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3) +{ + // 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(std::max(0.0, -q))); +} + +} // namespace o2::analysis::femto::femtomath + +#endif // PWGCF_FEMTO_CORE_FEMTOMATH_H_ diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 0d28ce5b99e..43b14727668 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -16,6 +16,7 @@ #ifndef PWGCF_FEMTO_CORE_PAIRHISTMANAGER_H_ #define PWGCF_FEMTO_CORE_PAIRHISTMANAGER_H_ +#include "PWGCF/Femto/Core/femtoMath.h" #include "PWGCF/Femto/Core/femtoSpherHarMath.h" #include "PWGCF/Femto/Core/femtoUtils.h" #include "PWGCF/Femto/Core/histManager.h" @@ -762,14 +763,14 @@ class PairHistManager mTrueParticle2 = ROOT::Math::PtEtaPhiMVector(mcParticle2.pt(), mcParticle2.eta(), mcParticle2.phi(), mPdgMass2); // compute true kinematics - mTrueKt = getKt(mTrueParticle1, mTrueParticle2); - mTrueMt = getMt(mTrueParticle1, mTrueParticle2); - mTrueMinv = getMinv(mTrueParticle1, mTrueParticle2); - mTrueKstar = getKstar(mTrueParticle1, mTrueParticle2); + mTrueKt = femtomath::getKt(mTrueParticle1, mTrueParticle2); + mTrueMt = femtomath::getMt(mTrueParticle1, mTrueParticle2, mMtType); + mTrueMinv = femtomath::getMinv(mTrueParticle1, mTrueParticle2); + mTrueKstar = femtomath::getKstar(mTrueParticle1, mTrueParticle2); if (mPlotBertschPratt) { std::tie(mTrueQout, mTrueQside, mTrueQlong) = - computeBertschPrattLCMS(mTrueParticle1, mTrueParticle2); + femtomath::computeBertschPrattLCMS(mTrueParticle1, mTrueParticle2); } } @@ -815,10 +816,10 @@ class PairHistManager mTrueParticle1 = ROOT::Math::PtEtaPhiMVector(particle1.pt(), particle1.eta(), particle1.phi(), mPdgMass1); mTrueParticle2 = ROOT::Math::PtEtaPhiMVector(particle2.pt(), particle2.eta(), particle2.phi(), mPdgMass2); - mTrueKt = getKt(mTrueParticle1, mTrueParticle2); - mTrueMt = getMt(mTrueParticle1, mTrueParticle2); - mTrueMinv = getMinv(mTrueParticle1, mTrueParticle2); - mTrueKstar = getKstar(mTrueParticle1, mTrueParticle2); + mTrueKt = femtomath::getKt(mTrueParticle1, mTrueParticle2); + mTrueMt = femtomath::getMt(mTrueParticle1, mTrueParticle2, mMtType); + mTrueMinv = femtomath::getMinv(mTrueParticle1, mTrueParticle2); + mTrueKstar = femtomath::getKstar(mTrueParticle1, mTrueParticle2); // in the pure mc-truth path there is no reco counterpart, so the generated values are also // stored in the reco members; they are used by getKinematic(), i.e. by the kinematic histograms @@ -829,7 +830,7 @@ class PairHistManager mKstar = mTrueKstar; if (mPlotBertschPratt) { - std::tie(mTrueQout, mTrueQside, mTrueQlong) = computeBertschPrattLCMS(mTrueParticle1, mTrueParticle2); + std::tie(mTrueQout, mTrueQside, mTrueQlong) = femtomath::computeBertschPrattLCMS(mTrueParticle1, mTrueParticle2); } if (mPlotDeltaEtaDeltaPhi) { mTrueDeltaEta = particle1.eta() - particle2.eta(); @@ -1314,27 +1315,27 @@ class PairHistManager mParticle2 = ROOT::Math::PtEtaPhiMVector(mAbsCharge2 * particle2.pt(), particle2.eta(), particle2.phi(), mass2); // set kT - mKt = getKt(mParticle1, mParticle2); + mKt = femtomath::getKt(mParticle1, mParticle2); // set mT - mMt = getMt(mParticle1, mParticle2); + mMt = femtomath::getMt(mParticle1, mParticle2, mMtType); // set Minv - mMassInv = getMinv(mParticle1, mParticle2); + mMassInv = femtomath::getMinv(mParticle1, mParticle2); // set kstar - mKstar = getKstar(mParticle1, mParticle2); + mKstar = femtomath::getKstar(mParticle1, mParticle2); if (mPlotBertschPratt) { - std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); + std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2); } if (mPlotEventShape) { if (!mPlotBertschPratt) { - std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); + std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2); } if constexpr (utils::HasEventShape) { - mPairPhiFromEventPlaneAngle = getPairPhiFromEventPlaneAngle(mParticle1, mParticle2, col); + mPairPhiFromEventPlaneAngle = femtomath::getPairPhiFromEventPlaneAngle(mParticle1, mParticle2, col); mQvector = col.qvec(); } else { LOG(fatal) << "plotEventShape is enabled but the collision table has no event-shape columns. Breaking..."; @@ -1342,7 +1343,7 @@ class PairHistManager } if (mPlotSH) { - std::tie(mShKv, mShOut, mShSide, mShLong) = computeShKinematics(mParticle1, mParticle2); + std::tie(mShKv, mShOut, mShSide, mShLong) = femtomath::computeShKinematics(mParticle1, mParticle2, mShFrame); } if (mPlotDeltaEtaDeltaPhi) { @@ -1351,10 +1352,10 @@ class PairHistManager } if (mPlotAlpha || mPlotKstarVsAlphaVsMtVsMult || mPlotKstarVsAlphaVsMtVsPt1VsPt2VsMult) { - mAlpha = getAlpha(mParticle1, mParticle2); + mAlpha = femtomath::getAlpha(mParticle1, mParticle2); } if (mPlotBetaLab || mPlotKstarVsBetaLabVsMtVsMult || mPlotKstarVsBetaLabVsMtVsPt1VsPt2VsMult) { - mBetaLab = getBetaLab(mParticle1, mParticle2); + mBetaLab = femtomath::getBetaLab(mParticle1, mParticle2); } if (mPlotDalitz) { @@ -1506,8 +1507,8 @@ class PairHistManager } if (mPlotSH) { const float shCentValue = mShUseCent ? mCent : mMult; - const int iCent = findShBin(shCentValue, mShCentEdges); - const int iKt = findShBin(mKt, mShKtEdges); + const int iCent = femtomath::findShBin(shCentValue, mShCentEdges); + const int iKt = femtomath::findShBin(mKt, mShKtEdges); if (iCent >= 0 && iKt >= 0) { mYlm.doYlmUpToL(mShLMax, mShOut, mShSide, mShLong, mShYlmBuffer.data()); for (std::size_t i = 0; i < mShYlmBuffer.size(); ++i) { @@ -1529,7 +1530,7 @@ class PairHistManager if (mShPlot1D) { mShBinCount[iCent][iKt]->Fill(mShKv, 1.0); // FemtoUniverse h1D = f3d[0]: qinv (=2k*) for identical-LCMS, else k*. - const float sh1DValue = (mShFrame == ShFrameLcmsIdentical) ? (2.0f * mKstar) : mKstar; + const float sh1DValue = (mShFrame == femtomath::kShFrameLcmsIdentical) ? (2.0f * mKstar) : mKstar; mSh1D[iCent][iKt]->Fill(sh1DValue); } } @@ -1606,212 +1607,6 @@ class PairHistManager } } - float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - auto sum = (part1 + part2); - double kt = 0.5 * sum.Pt(); - return static_cast(kt); - } - - float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - auto sum = part1 + part2; - double mt = 0; - double averageMass = 0; - double reducedMass = 0; - switch (mMtType) { - case modes::TransverseMassType::kAveragePdgMass: - averageMass = 0.5 * (part1.M() + part2.M()); - mt = std::hypot(0.5 * sum.Pt(), averageMass); - break; - case modes::TransverseMassType::kReducedPdgMass: - reducedMass = 2. * (part1.M() * part2.M()) / (part1.M() + part2.M()); - mt = std::hypot(0.5 * sum.Pt(), reducedMass); - break; - case modes::TransverseMassType::kMt4Vector: - mt = 0.5 * sum.Mt(); - break; - default: - LOG(fatal) << "Invalid transverse mass type, breaking..."; - } - return static_cast(mt); - } - - float getMinv(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - return static_cast((part1 + part2).M()); - } - - float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - // 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)); - } - - float getAlpha(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - return static_cast(RecoDecay::constrainAngle(part2.phi() - part1.phi(), -o2::constants::math::PI)); - } - - float getBetaLab(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - ROOT::Math::XYZVector p1(part1.Px(), part1.Py(), part1.Pz()); - ROOT::Math::XYZVector p2(part2.Px(), part2.Py(), part2.Pz()); - - ROOT::Math::XYZVector kstar = p1 - p2; - ROOT::Math::XYZVector kt = p1 + p2; - - return static_cast(std::acos(kstar.Dot(kt) / std::sqrt(kstar.Mag2() * kt.Mag2()))); - } - - std::tuple computeBertschPrattLCMS(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - const ROOT::Math::PxPyPzEVector p1(part1); - const ROOT::Math::PxPyPzEVector p2(part2); - const ROOT::Math::PxPyPzEVector pSum = p1 + p2; - - const double tPx = pSum.Px(); - const double tPy = pSum.Py(); - const double tPz = pSum.Pz(); - const double tE = pSum.E(); - - const double tPt = std::sqrt(tPx * tPx + tPy * tPy); - const double tMt = std::sqrt(tE * tE - tPz * tPz); - - static constexpr double MinTransverseMomentum = 1e-9; - if (tPt < MinTransverseMomentum || tMt < MinTransverseMomentum) { - return {0.0f, 0.0f, 0.0f}; - } - - const double betaL = tPz / tE; - const double gammaL = tE / tMt; - - const double kout1 = (p1.Px() * tPx + p1.Py() * tPy) / tPt; - const double kside1 = (-p1.Px() * tPy + p1.Py() * tPx) / tPt; - const double klong1 = gammaL * (p1.Pz() - betaL * p1.E()); - - const double kout2 = (p2.Px() * tPx + p2.Py() * tPy) / tPt; - const double kside2 = (p2.Py() * tPx - p2.Px() * tPy) / tPt; - const double klong2 = gammaL * (p2.Pz() - betaL * p2.E()); - - auto qOut = static_cast(kout1 - kout2); - auto qSide = static_cast(kside1 - kside2); - auto qLong = static_cast(klong1 - klong2); - - return {qOut, qSide, qLong}; - } - - // Return the bin index for value given ascending bin edges, or -1 if out of range. - // edges = {e0, e1, ..., eN} defines N bins [e0,e1), [e1,e2), ..., [e_{N-1},eN). - static int findShBin(double value, std::vector const& edges) - { - static constexpr std::size_t MinEdgesForOneBin = 2; - if (edges.size() < MinEdgesForOneBin || value < edges.front() || value >= edges.back()) { - return -1; - } - for (std::size_t i = 0; i < edges.size() - 1; ++i) { - if (value >= edges[i] && value < edges[i + 1]) { - return static_cast(i); - } - } - return -1; - } - - // Kinematics feeding the spherical-harmonics decomposition, ported 1:1 from - // FemtoUniverseMath::newpairfunc so results are numerically comparable. - // Returns {kv, out, side, long}, selected by mShFrame: - // ShFrameLcmsNonIdentical: {kstar, fDKOut, fDKSide, fDKLong} (LCMS components of particle 1) - // ShFrameLcmsIdentical: {qinv, outLCMS, sideLCMS, longLCMS} (LCMS q-differences) - // ShFramePrf: {|qPRF|, outPRF, sidePRF, longPRF} (matches FemtoUniverse isIdenPRF=true) - std::tuple computeShKinematics(ROOT::Math::PtEtaPhiMVector const& part1, - ROOT::Math::PtEtaPhiMVector const& part2) - { - const ROOT::Math::PxPyPzEVector p1(part1); - const ROOT::Math::PxPyPzEVector p2(part2); - const ROOT::Math::PxPyPzEVector sum = p1 + p2; - - const double tPx = sum.Px(); - const double tPy = sum.Py(); - const double tPz = sum.Pz(); - const double tE = sum.E(); - const double tPtSq = tPx * tPx + tPy * tPy; - const double tMtSq = tE * tE - tPz * tPz; - - static constexpr double MinTransverseMomentum = 1e-9; - if (tPtSq < MinTransverseMomentum || tMtSq < MinTransverseMomentum) { - return {0.f, 0.f, 0.f, 0.f}; - } - const double tPt = std::sqrt(tPtSq); - const double tMt = std::sqrt(tMtSq); - const double tM = std::sqrt(std::max(0.0, tMtSq - tPtSq)); - - const double beta = tPz / tE; - const double gamma = tE / tMt; - - // LCMS components of particle 1 (fDKOut/Side/Long in newpairfunc) - const double e1 = p1.E(); - const double fDKOut = (p1.Px() * tPx + p1.Py() * tPy) / tPt; - const double fDKSide = (-p1.Px() * tPy + p1.Py() * tPx) / tPt; - const double fDKLong = gamma * (p1.Pz() - beta * e1); - const double fDE = gamma * (e1 - beta * p1.Pz()); - - // LCMS components of particle 2 - const double e2 = p2.E(); - const double px2 = (p2.Px() * tPx + p2.Py() * tPy) / tPt; - const double py2 = (p2.Py() * tPx - p2.Px() * tPy) / tPt; - const double pz2 = gamma * (p2.Pz() - beta * e2); - - const double outLCMS = fDKOut - px2; - const double sideLCMS = fDKSide - py2; - const double longLCMS = fDKLong - pz2; - - if (mShFrame == ShFrameLcmsIdentical) { - // LCMS identical: Ylm <- LCMS q-differences, axis <- qinv - const double qinv = std::sqrt(outLCMS * outLCMS + sideLCMS * sideLCMS + longLCMS * longLCMS); - return {static_cast(qinv), - static_cast(outLCMS), static_cast(sideLCMS), static_cast(longLCMS)}; - } - if (mShFrame == ShFramePrf) { - // PRF: Ylm <- PRF q-differences, axis <- |q_PRF| (matches FemtoUniverse isIdenPRF=true) - const double pE2LCMS = gamma * (e2 - beta * p2.Pz()); - const double betaOut = tPt / tMt; - const double gammaOut = tMt / tM; - const double outPRF = gammaOut * (outLCMS - betaOut * (fDE - pE2LCMS)); - const double sidePRF = sideLCMS; - const double longPRF = longLCMS; - const double qPRF = std::sqrt(outPRF * outPRF + sidePRF * sidePRF + longPRF * longPRF); - return {static_cast(qPRF), - static_cast(outPRF), static_cast(sidePRF), static_cast(longPRF)}; - } - // LCMS non-identical (default): Ylm <- LCMS components of particle 1, axis <- k* - const double betaOut = tPt / tMt; - const double gammaOut = tMt / tM; - const double fKOut = gammaOut * (fDKOut - betaOut * fDE); - const double kstar = std::sqrt(fKOut * fKOut + fDKSide * fDKSide + fDKLong * fDKLong); - return {static_cast(kstar), - static_cast(fDKOut), static_cast(fDKSide), static_cast(fDKLong)}; - } - - template - float getPairPhiFromEventPlaneAngle(const ROOT::Math::PtEtaPhiMVector& part1, const ROOT::Math::PtEtaPhiMVector& part2, const T1& col) - { - const ROOT::Math::PtEtaPhiMVector trackSum = part1 + part2; - float phi = RecoDecay::constrainAngle(trackSum.Phi() - col.eventPlaneAngle()); - return phi; - } - o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; bool mUsePdgMass = true; double mPdgMass1 = 0.; @@ -1929,9 +1724,6 @@ class PairHistManager bool mShUseCent = false; int mShLMax = 1; int mShFrame = 0; - static constexpr int ShFrameLcmsNonIdentical = 0; - static constexpr int ShFrameLcmsIdentical = 1; - static constexpr int ShFramePrf = 2; o2::framework::AxisSpec mShKstarSpec{{60, 0.0f, 0.3f}, "k* (GeV/#it{c})"}; // set in init() o2::framework::AxisSpec mShCentSpec{{1, 0.0f, 200.0f}, "mult/cent"}; diff --git a/PWGCF/Femto/Core/tripletHistManager.h b/PWGCF/Femto/Core/tripletHistManager.h index 1bc69b06eeb..abf5dcf3e42 100644 --- a/PWGCF/Femto/Core/tripletHistManager.h +++ b/PWGCF/Femto/Core/tripletHistManager.h @@ -16,6 +16,7 @@ #ifndef PWGCF_FEMTO_CORE_TRIPLETHISTMANAGER_H_ #define PWGCF_FEMTO_CORE_TRIPLETHISTMANAGER_H_ +#include "PWGCF/Femto/Core/femtoMath.h" #include "PWGCF/Femto/Core/femtoUtils.h" #include "PWGCF/Femto/Core/histManager.h" #include "PWGCF/Femto/Core/modes.h" @@ -67,6 +68,7 @@ enum TripletHist { kQ3VsMtVsMultVsCent, kQ3VsMtVsPt1VsPt2VsPt3VsMult, kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, + kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, // higher dimensions with mass // mc kTrueQ3VsQ3, @@ -131,6 +133,7 @@ struct ConfTripletBinning : o2::framework::ConfigurableGroup { o2::framework::Configurable plotQ3VsMtVsMultVsCent{"plotQ3VsMtVsMultVsCent", false, "Enable 4D histogram (Q3 vs Mt vs Mult Vs Cent)"}; o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMult{"plotQ3VsMtVsPt1VsPt2VsPt3VsMult", false, "Enable 6D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult)"}; o2::framework::Configurable plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent{"plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent", false, "Enable 7D histogram (Q3 vs Mt Vs Pt1 vs Pt2 vs Pt3 vs Mult vs Cent)"}; + o2::framework::Configurable plotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3{"plotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3", false, "Enable 7D histogram (Q3 vs kstar12 vs kstar13 vs kstar23 vs Pt1 vs Pt2 vs Pt3)"}; o2::framework::ConfigurableAxis q3{"q3", {{600, 0, 6}}, "q3"}; o2::framework::ConfigurableAxis kstar{"kstar", {{600, 0, 6}}, "kstar (between pairs of particles)"}; o2::framework::ConfigurableAxis mt{"mt", {{500, 0.8, 5.8}}, "mt"}; @@ -152,6 +155,12 @@ struct ConfTripletCuts : o2::framework::ConfigurableGroup { o2::framework::Configurable q3Min{"q3Min", -1, "Minimal Q3 (set to -1 to deactivate)"}; o2::framework::Configurable mtMax{"mtMax", -1, "Maximal mt (set to -1 to deactivate)"}; o2::framework::Configurable mtMin{"mtMin", -1, "Minimal mt (set to -1 to deactivate)"}; + o2::framework::Configurable kstar12Min{"kstar12Min", -1, "Minimal k* between particle 1 and particle 2 (set to -1 to deactivate)"}; + o2::framework::Configurable kstar12Max{"kstar12Max", -1, "Maximal k* between particle 1 and particle 2 (set to -1 to deactivate)"}; + o2::framework::Configurable kstar13Min{"kstar13Min", -1, "Minimal k* between particle 1 and particle 3 (set to -1 to deactivate)"}; + o2::framework::Configurable kstar13Max{"kstar13Max", -1, "Maximal k* between particle 1 and particle 3 (set to -1 to deactivate)"}; + o2::framework::Configurable kstar23Min{"kstar23Min", -1, "Minimal k* between particle 2 and particle 3 (set to -1 to deactivate)"}; + o2::framework::Configurable kstar23Max{"kstar23Max", -1, "Maximal k* between particle 2 and particle 3 (set to -1 to deactivate)"}; o2::framework::Configurable mixOnlyCommonAncestor{"mixOnlyCommonAncestor", false, "Require pair to have common ancestor (in the same event)"}; o2::framework::Configurable mixOnlyNonCommonAncestor{"mixOnlyNonCommonAncestor", false, "Require pair to have non-common ancestor (in the same event)"}; o2::framework::Configurable useMotherAsAncestor{"useMotherAsAncestor", false, "Use the first ancestor (i.e. the direct mother) instead of the partonic mother when requiring (non-)common ancestry"}; @@ -182,6 +191,7 @@ constexpr std::array, kTripletHistogramLast> // n-D with mass {kQ3VsMtVsPt1VsPt2VsPt3VsMult, o2::framework::HistType::kTHnSparseF, "hQ3VsMtVsPt1VsPt2VsPt3VsMult", "Q_{3} vs m_{T} vs p_{T,1} vs p_{T,2} vs p_{T,3} vs multiplicity; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); p_{T,1} (GeV/#it{c}) ; p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c}); Multiplicity;"}, {kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, o2::framework::HistType::kTHnSparseF, "hQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent", "Q_{3} vs m_{T} vs p_{T,1} vs p_{T,2} vs p_{T,3} vs multiplicity vs centrality; Q_{3} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); p_{T,1} (GeV/#it{c}) ; p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c}); Multiplicity; Centrality (%);"}, + {kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, o2::framework::HistType::kTHnSparseF, "hQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3", "Q_{3} vs k*_{12} vs k*_{13} vs k*_{23} vs p_{T,1} vs p_{T,2} vs p_{T,3}; Q_{3} (GeV/#it{c}); k*_{12} (GeV/#it{c}); k*_{13} (GeV/#it{c}); k*_{23} (GeV/#it{c}); p_{T,1} (GeV/#it{c}); p_{T,2} (GeV/#it{c}); p_{T,3} (GeV/#it{c});"}, {kTrueQ3VsQ3, o2::framework::HistType::kTH2F, "hTrueQ3VsQ3", "Q_{3,True} vs Q_{3}; Q_{3,true} (GeV/#it{c}); Q_{3} (GeV/#it{c})"}, {kTrueMtVsMt, o2::framework::HistType::kTH2F, "hTrueMtVsMt", "m_{T,True} vs m_{T}; m_{T,True} (GeV/#it{c}^{2}); m_{T} (GeV/#it{c}^{2})"}, {kTrueMultVsMult, o2::framework::HistType::kTH2F, "hTrueMultVsMult", "Multiplicity_{True} vs Multiplicity; Multiplicity_{True} ; Multiplicity"}, @@ -207,24 +217,25 @@ constexpr std::array, kTripletHistogramLast> }}; // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) -#define TRIPLET_HIST_ANALYSIS_MAP(conf, confMixing) \ - {kQ3, {(conf).q3}}, \ - {kMt, {(conf).mt}}, \ - {kKstar12, {(conf).kstar}}, \ - {kKstar13, {(conf).kstar}}, \ - {kKstar23, {(conf).kstar}}, \ - {kPt1VsQ3, {(conf).pt1, (conf).q3}}, \ - {kPt2VsQ3, {(conf).pt2, (conf).q3}}, \ - {kPt3VsQ3, {(conf).pt3, (conf).q3}}, \ - {kQ3VsMt, {(conf).q3, (conf).mt}}, \ - {kQ3VsMult, {(conf).q3, (conf).multiplicity}}, \ - {kQ3VsCent, {(conf).q3, (conf).centrality}}, \ - {kPt1VsPt2VsPt3, {(conf).pt1, (conf).pt2, (conf).pt3}}, \ - {kQ3VsPt1VsPt2VsPt3, {(conf).q3, (conf).pt1, (conf).pt2, (conf).pt3}}, \ - {kQ3VsMtVsMult, {(conf).q3, (conf).mt, (conf).multiplicity}}, \ - {kQ3VsMtVsMultVsCent, {(conf).q3, (conf).mt, (conf).multiplicity, (conf).centrality}}, \ - {kQ3VsMtVsPt1VsPt2VsPt3VsMult, {(conf).q3, (conf).mt, (conf).pt1, (conf).pt2, (conf).pt3, (conf).multiplicity}}, \ - {kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, {(conf).q3, (conf).mt, (conf).pt1, (conf).pt2, (conf).pt3, (conf).multiplicity, (conf).centrality}}, +#define TRIPLET_HIST_ANALYSIS_MAP(conf, confMixing) \ + {kQ3, {(conf).q3}}, \ + {kMt, {(conf).mt}}, \ + {kKstar12, {(conf).kstar}}, \ + {kKstar13, {(conf).kstar}}, \ + {kKstar23, {(conf).kstar}}, \ + {kPt1VsQ3, {(conf).pt1, (conf).q3}}, \ + {kPt2VsQ3, {(conf).pt2, (conf).q3}}, \ + {kPt3VsQ3, {(conf).pt3, (conf).q3}}, \ + {kQ3VsMt, {(conf).q3, (conf).mt}}, \ + {kQ3VsMult, {(conf).q3, (conf).multiplicity}}, \ + {kQ3VsCent, {(conf).q3, (conf).centrality}}, \ + {kPt1VsPt2VsPt3, {(conf).pt1, (conf).pt2, (conf).pt3}}, \ + {kQ3VsPt1VsPt2VsPt3, {(conf).q3, (conf).pt1, (conf).pt2, (conf).pt3}}, \ + {kQ3VsMtVsMult, {(conf).q3, (conf).mt, (conf).multiplicity}}, \ + {kQ3VsMtVsMultVsCent, {(conf).q3, (conf).mt, (conf).multiplicity, (conf).centrality}}, \ + {kQ3VsMtVsPt1VsPt2VsPt3VsMult, {(conf).q3, (conf).mt, (conf).pt1, (conf).pt2, (conf).pt3, (conf).multiplicity}}, \ + {kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, {(conf).q3, (conf).mt, (conf).pt1, (conf).pt2, (conf).pt3, (conf).multiplicity, (conf).centrality}}, \ + {kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, {(conf).q3, (conf).kstar, (conf).kstar, (conf).kstar, (conf).pt1, (conf).pt2, (conf).pt3}}, // NOLINTNEXTLINE(cppcoreguidelines-macro-usage) #define TRIPLET_HIST_MIXING_QA_MAP(confMixing) \ @@ -334,6 +345,12 @@ class TripletHistManager mQ3Max = ConfTripletCuts.q3Max.value; mMtMin = ConfTripletCuts.mtMin.value; mMtMax = ConfTripletCuts.mtMax.value; + mKstar12Min = ConfTripletCuts.kstar12Min.value; + mKstar12Max = ConfTripletCuts.kstar12Max.value; + mKstar13Min = ConfTripletCuts.kstar13Min.value; + mKstar13Max = ConfTripletCuts.kstar13Max.value; + mKstar23Min = ConfTripletCuts.kstar23Min.value; + mKstar23Max = ConfTripletCuts.kstar23Max.value; mPlotKstar = ConfTripletBinning.plotKstar.value; mPlotPt1VsPt2VsPt3 = ConfTripletBinning.plotPt1VsPt2VsPt3.value; @@ -342,6 +359,15 @@ class TripletHistManager mPlotQ3VsMtVsMultVsCent = ConfTripletBinning.plotQ3VsMtVsMultVsCent.value; mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult = ConfTripletBinning.plotQ3VsMtVsPt1VsPt2VsPt3VsMult.value; mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent = ConfTripletBinning.plotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent.value; + mPlotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3 = ConfTripletBinning.plotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3.value; + + // kstar12/13/23 are only computed on demand (they cost 3 extra femtomath::getKstar() calls per + // triplet); anything that needs the value -- the 1D plots, the 7D histogram, or a cut -- + // must be OR'ed in here, otherwise it silently stays at its default of 0 + mRequireKstar = mPlotKstar || mPlotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3 || + mKstar12Min > 0.f || mKstar12Max > 0.f || + mKstar13Min > 0.f || mKstar13Max > 0.f || + mKstar23Min > 0.f || mKstar23Max > 0.f; mPairCorrelationQa = ConfMixing.enablePairCorrelationQa.value; mEventMixingQa = ConfMixing.enableEventMixingQa.value; @@ -390,12 +416,12 @@ class TripletHistManager // set mT mMt = getMt(mParticle1, mParticle2, mParticle3); // set Q3 - mQ3 = getQ3(mParticle1, mParticle2, mParticle3); + mQ3 = femtomath::getQ3(mParticle1, mParticle2, mParticle3); - if (mPlotKstar) { - mKstar12 = getKstar(mParticle1, mParticle2); - mKstar13 = getKstar(mParticle1, mParticle3); - mKstar23 = getKstar(mParticle2, mParticle3); + if (mRequireKstar) { + mKstar12 = femtomath::getKstar(mParticle1, mParticle2); + mKstar13 = femtomath::getKstar(mParticle1, mParticle3); + mKstar23 = femtomath::getKstar(mParticle2, mParticle3); } // if one of the particles has a mass getter, we cache the value for the filling later @@ -437,7 +463,7 @@ class TripletHistManager mTrueParticle3 = ROOT::Math::PtEtaPhiMVector(particle3.pt(), particle3.eta(), particle3.phi(), mPdgMass3); mTrueMt = getMt(mTrueParticle1, mTrueParticle2, mTrueParticle3); - mTrueQ3 = getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3); + mTrueQ3 = femtomath::getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3); // in the pure mc-truth path there is no reco counterpart, so the generated values are also // stored in the reco members; they are used by getKinematic(), i.e. by the kinematic histograms @@ -445,10 +471,10 @@ class TripletHistManager mMt = mTrueMt; mQ3 = mTrueQ3; - if (mPlotKstar) { - mTrueKstar12 = getKstar(mTrueParticle1, mTrueParticle2); - mTrueKstar13 = getKstar(mTrueParticle1, mTrueParticle3); - mTrueKstar23 = getKstar(mTrueParticle2, mTrueParticle3); + if (mRequireKstar) { + mTrueKstar12 = femtomath::getKstar(mTrueParticle1, mTrueParticle2); + mTrueKstar13 = femtomath::getKstar(mTrueParticle1, mTrueParticle3); + mTrueKstar23 = femtomath::getKstar(mTrueParticle2, mTrueParticle3); mKstar12 = mTrueKstar12; mKstar13 = mTrueKstar13; mKstar23 = mTrueKstar23; @@ -493,7 +519,7 @@ class TripletHistManager // set true mT mTrueMt = getMt(mTrueParticle1, mTrueParticle2, mTrueParticle3); // set true Q3 - mTrueQ3 = getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3); + mTrueQ3 = femtomath::getQ3(mTrueParticle1, mTrueParticle2, mTrueParticle3); } template @@ -535,7 +561,13 @@ class TripletHistManager return (!(mQ3Min > 0.f) || mTrueQ3 > mQ3Min) && (!(mQ3Max > 0.f) || mTrueQ3 < mQ3Max) && (!(mMtMin > 0.f) || mTrueMt > mMtMin) && - (!(mMtMax > 0.f) || mTrueMt < mMtMax); + (!(mMtMax > 0.f) || mTrueMt < mMtMax) && + (!(mKstar12Min > 0.f) || mTrueKstar12 > mKstar12Min) && + (!(mKstar12Max > 0.f) || mTrueKstar12 < mKstar12Max) && + (!(mKstar13Min > 0.f) || mTrueKstar13 > mKstar13Min) && + (!(mKstar13Max > 0.f) || mTrueKstar13 < mKstar13Max) && + (!(mKstar23Min > 0.f) || mTrueKstar23 > mKstar23Min) && + (!(mKstar23Max > 0.f) || mTrueKstar23 < mKstar23Max); } bool checkTripletCuts() const @@ -543,7 +575,13 @@ class TripletHistManager return (!(mQ3Min > 0.f) || mQ3 > mQ3Min) && (!(mQ3Max > 0.f) || mQ3 < mQ3Max) && (!(mMtMin > 0.f) || mMt > mMtMin) && - (!(mMtMax > 0.f) || mMt < mMtMax); + (!(mMtMax > 0.f) || mMt < mMtMax) && + (!(mKstar12Min > 0.f) || mKstar12 > mKstar12Min) && + (!(mKstar12Max > 0.f) || mKstar12 < mKstar12Max) && + (!(mKstar13Min > 0.f) || mKstar13 > mKstar13Min) && + (!(mKstar13Max > 0.f) || mKstar13 < mKstar13Max) && + (!(mKstar23Min > 0.f) || mKstar23 > mKstar23Min) && + (!(mKstar23Max > 0.f) || mKstar23 < mKstar23Max); } template @@ -630,46 +668,6 @@ class TripletHistManager } private: - ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PxPyPzEVector const& vi, ROOT::Math::PxPyPzEVector const& vj) - { - auto trackSum = vi + vj; - auto trackDifference = vi - vj; - 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) - { - // 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(std::max(0.0, -q))); - } - - float getKstar(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - // 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)); - } - void initAnalysis(std::map> const& Specs) { std::string analysisDir = std::string(prefix) + std::string(AnalysisDir); @@ -711,6 +709,9 @@ class TripletHistManager if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent) { mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), getHistDesc(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), getHistType(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable), {Specs.at(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent)}); } + if (mPlotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3) { + mHistogramRegistry->add(analysisDir + getHistNameV2(kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, HistTable), getHistDesc(kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, HistTable), getHistType(kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, HistTable), {Specs.at(kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3)}); + } } void initMc(std::map> const& Specs) @@ -827,6 +828,9 @@ class TripletHistManager if (mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent) { mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent, HistTable)), mQ3, mMt, mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt(), mMult, mCent); } + if (mPlotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3) { + mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3, HistTable)), mQ3, mKstar12, mKstar13, mKstar23, mParticle1.Pt(), mParticle2.Pt(), mParticle3.Pt()); + } } void fillMc() @@ -841,33 +845,9 @@ class TripletHistManager } } - double getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) - { - auto sum = part1 + part2; - double mt = 0.; - double averageMass = 0.; - double reducedMass = 0.; - switch (mMtType) { - case modes::TransverseMassType::kAveragePdgMass: - averageMass = 0.5 * (part1.M() + part2.M()); - mt = std::hypot(0.5 * sum.Pt(), averageMass); - break; - case modes::TransverseMassType::kReducedPdgMass: - reducedMass = 2. * (part1.M() * part2.M() / (part1.M() + part2.M())); - mt = std::hypot(0.5 * sum.Pt(), reducedMass); - break; - case modes::TransverseMassType::kMt4Vector: - mt = 0.5 * sum.Mt(); - break; - default: - LOG(fatal) << "Invalid transverse mass type, breaking..."; - } - return mt; - } - float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3) { - return static_cast((getMt(part1, part2) + getMt(part2, part3) + getMt(part1, part3)) / 3.); + return (femtomath::getMt(part1, part2, mMtType) + femtomath::getMt(part2, part3, mMtType) + femtomath::getMt(part1, part3, mMtType)) / 3.f; } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; @@ -914,18 +894,26 @@ class TripletHistManager float mQ3Max = -1.f; float mMtMin = -1.f; float mMtMax = -1.f; + float mKstar12Min = -1.f; + float mKstar12Max = -1.f; + float mKstar13Min = -1.f; + float mKstar13Max = -1.f; + float mKstar23Min = -1.f; + float mKstar23Max = -1.f; // flags bool mPlot1d = true; bool mPlot2d = true; bool mPlotKstar = false; + bool mRequireKstar = false; bool mPlotPt1VsPt2VsPt3 = false; bool mPlotQ3VsPt1VsPt2VsPt3 = false; bool mPlotQ3VsMtVsMult = false; bool mPlotQ3VsMtVsMultVsCent = false; bool mPlotQ3VsMtVsPt1VsPt2VsPt3VsMult = false; bool mPlotQ3VsMtVsPt1VsPt2VsPt3VsMultVsCent = false; + bool mPlotQ3VsKstar12VsKstar13VsKstar23VsPt1VsPt2VsPt3 = false; // mixing qa bool mPairCorrelationQa = false; From 01e6d3f83c64c6470c59507af87b44023a02468d Mon Sep 17 00:00:00 2001 From: Anton Riedel Date: Sun, 1 Nov 2026 22:41:06 +0100 Subject: [PATCH 2/3] Fix: fix linter warnings --- PWGCF/Femto/Core/femtoMath.h | 16 ++++++++-------- PWGCF/Femto/Core/femtoSpherHarMath.h | 11 ++++++++--- PWGCF/Femto/Core/pairHistManager.h | 2 +- 3 files changed, 17 insertions(+), 12 deletions(-) diff --git a/PWGCF/Femto/Core/femtoMath.h b/PWGCF/Femto/Core/femtoMath.h index ead624399bc..84ef0278d9a 100644 --- a/PWGCF/Femto/Core/femtoMath.h +++ b/PWGCF/Femto/Core/femtoMath.h @@ -38,9 +38,9 @@ namespace o2::analysis::femto::femtomath { // SH reference frame selectors for computeShKinematics(), mirrors ConfPairBinning.shFrame -constexpr int kShFrameLcmsNonIdentical = 0; -constexpr int kShFrameLcmsIdentical = 1; -constexpr int kShFramePrf = 2; +constexpr int ShFrameLcmsNonIdentical = 0; +constexpr int ShFrameLcmsIdentical = 1; +constexpr int ShFramePrf = 2; inline float getKt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2) { @@ -173,9 +173,9 @@ inline int findShBin(double value, std::vector const& edges) // Kinematics feeding the spherical-harmonics decomposition, ported 1:1 from // FemtoUniverseMath::newpairfunc so results are numerically comparable. // Returns {kv, out, side, long}, selected by shFrame: -// kShFrameLcmsNonIdentical: {kstar, fDKOut, fDKSide, fDKLong} (LCMS components of particle 1) -// kShFrameLcmsIdentical: {qinv, outLCMS, sideLCMS, longLCMS} (LCMS q-differences) -// kShFramePrf: {|qPRF|, outPRF, sidePRF, longPRF} (matches FemtoUniverse isIdenPRF=true) +// ShFrameLcmsNonIdentical: {kstar, fDKOut, fDKSide, fDKLong} (LCMS components of particle 1) +// ShFrameLcmsIdentical: {qinv, outLCMS, sideLCMS, longLCMS} (LCMS q-differences) +// ShFramePrf: {|qPRF|, outPRF, sidePRF, longPRF} (matches FemtoUniverse isIdenPRF=true) inline std::tuple computeShKinematics(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, int shFrame) @@ -219,13 +219,13 @@ inline std::tuple computeShKinematics(ROOT::Math::Pt const double sideLCMS = fDKSide - py2; const double longLCMS = fDKLong - pz2; - if (shFrame == kShFrameLcmsIdentical) { + if (shFrame == ShFrameLcmsIdentical) { // LCMS identical: Ylm <- LCMS q-differences, axis <- qinv const double qinv = std::sqrt(outLCMS * outLCMS + sideLCMS * sideLCMS + longLCMS * longLCMS); return {static_cast(qinv), static_cast(outLCMS), static_cast(sideLCMS), static_cast(longLCMS)}; } - if (shFrame == kShFramePrf) { + if (shFrame == ShFramePrf) { // PRF: Ylm <- PRF q-differences, axis <- |q_PRF| (matches FemtoUniverse isIdenPRF=true) const double pE2LCMS = gamma * (e2 - beta * p2.Pz()); const double betaOut = tPt / tMt; diff --git a/PWGCF/Femto/Core/femtoSpherHarMath.h b/PWGCF/Femto/Core/femtoSpherHarMath.h index 717368f830d..65cd7e98dd5 100644 --- a/PWGCF/Femto/Core/femtoSpherHarMath.h +++ b/PWGCF/Femto/Core/femtoSpherHarMath.h @@ -34,6 +34,11 @@ class SpherHarMath static constexpr int TrigCacheSize = MaxSupportedL + 1; // sin/cos powers cache size (indices 0..5) static constexpr double SmallLength = 1e-10; // numerical guard for |r|, |z| + // thresholds for legendreUpToYlm(): the l=N block is only filled once lmax reaches N + static constexpr int MinLmaxForL3Block = 3; + static constexpr int MinLmaxForL4Block = 4; + static constexpr int MinLmaxForL5Block = 5; + /// Values of various coefficients void initializeYlms() { @@ -138,7 +143,7 @@ class SpherHarMath } // Legendre polynomials l=3 - if (lmax > 2) { // o2-linter: disable=magic-number (l index, mirrors lmax>0/1 above) + if (lmax >= MinLmaxForL3Block) { lbuf[6] = sins[3]; lbuf[7] = sins[2] * coss[1]; lbuf[8] = (5 * coss[2] - 1) * sins[1]; @@ -146,7 +151,7 @@ class SpherHarMath } // Legendre polynomials l=4 - if (lmax > 3) { // o2-linter: disable=magic-number (l index) + if (lmax >= MinLmaxForL4Block) { lbuf[10] = sins[4]; lbuf[11] = sins[3] * coss[1]; lbuf[12] = (7 * coss[2] - 1) * sins[2]; @@ -155,7 +160,7 @@ class SpherHarMath } // Legendre polynomials l=5 - if (lmax > 4) { // o2-linter: disable=magic-number (l index) + if (lmax >= MinLmaxForL5Block) { lbuf[15] = sins[5]; lbuf[16] = sins[4] * coss[1]; lbuf[17] = (9 * coss[2] - 1) * sins[3]; diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 43b14727668..6266d60cebd 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -1530,7 +1530,7 @@ class PairHistManager if (mShPlot1D) { mShBinCount[iCent][iKt]->Fill(mShKv, 1.0); // FemtoUniverse h1D = f3d[0]: qinv (=2k*) for identical-LCMS, else k*. - const float sh1DValue = (mShFrame == femtomath::kShFrameLcmsIdentical) ? (2.0f * mKstar) : mKstar; + const float sh1DValue = (mShFrame == femtomath::ShFrameLcmsIdentical) ? (2.0f * mKstar) : mKstar; mSh1D[iCent][iKt]->Fill(sh1DValue); } } From 6ffeccfc2491c36853f0ab81aedf5764bb6bf3a9 Mon Sep 17 00:00:00 2001 From: Anton Riedel Date: Sun, 1 Nov 2026 23:34:53 +0100 Subject: [PATCH 3/3] Fix: fix code check errors --- PWGCF/Femto/Core/femtoMath.h | 1 + PWGCF/Femto/Core/pairHistManager.h | 3 --- PWGCF/Femto/Core/tripletHistManager.h | 3 ++- 3 files changed, 3 insertions(+), 4 deletions(-) diff --git a/PWGCF/Femto/Core/femtoMath.h b/PWGCF/Femto/Core/femtoMath.h index 84ef0278d9a..f686b85d9b9 100644 --- a/PWGCF/Femto/Core/femtoMath.h +++ b/PWGCF/Femto/Core/femtoMath.h @@ -31,6 +31,7 @@ #include #include +#include #include #include diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 6266d60cebd..5d227b7500d 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -30,14 +30,11 @@ #include #include -#include // IWYU pragma: keep (do not replace with Math/Vector3Dfwd.h) -#include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include #include -#include #include #include #include diff --git a/PWGCF/Femto/Core/tripletHistManager.h b/PWGCF/Femto/Core/tripletHistManager.h index abf5dcf3e42..7ffc306308c 100644 --- a/PWGCF/Femto/Core/tripletHistManager.h +++ b/PWGCF/Femto/Core/tripletHistManager.h @@ -29,7 +29,6 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include -#include #include #include #include @@ -847,6 +846,8 @@ class TripletHistManager float getMt(ROOT::Math::PtEtaPhiMVector const& part1, ROOT::Math::PtEtaPhiMVector const& part2, ROOT::Math::PtEtaPhiMVector const& part3) { + // all three pairwise combinations are intentional, not a mixed-up argument order + // NOLINTNEXTLINE(readability-suspicious-call-argument) return (femtomath::getMt(part1, part2, mMtType) + femtomath::getMt(part2, part3, mMtType) + femtomath::getMt(part1, part3, mMtType)) / 3.f; }