From 40093f2a3d6771c5de05294af2c7a296d1c2d3a5 Mon Sep 17 00:00:00 2001 From: ptak Date: Thu, 9 Jul 2026 09:34:43 +0200 Subject: [PATCH] [PWGCF] Add spherical harmonics decomposition to unified Femto framework Signed-off-by: ptak --- PWGCF/Femto/Core/femtoSpherHarMath.h | 225 +++++++++++++++++++++++++++ PWGCF/Femto/Core/pairHistManager.h | 212 +++++++++++++++++++++++++ 2 files changed, 437 insertions(+) create mode 100644 PWGCF/Femto/Core/femtoSpherHarMath.h diff --git a/PWGCF/Femto/Core/femtoSpherHarMath.h b/PWGCF/Femto/Core/femtoSpherHarMath.h new file mode 100644 index 00000000000..717368f830d --- /dev/null +++ b/PWGCF/Femto/Core/femtoSpherHarMath.h @@ -0,0 +1,225 @@ +// 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 femtoSpherHarMath.h +/// \brief Container for the calculation of spherical harmonics components +/// \author Pritam Chakraborty, WUT Warsaw, pritam.chakraborty@pw.edu.pl + +#ifndef PWGCF_FEMTO_CORE_FEMTOSPHERHARMATH_H_ +#define PWGCF_FEMTO_CORE_FEMTOSPHERHARMATH_H_ + +#include + +#include +#include +#include + +namespace o2::analysis::femto +{ + +/// \class SpherHarMath +/// \brief Container for math calculations of quantities related to pairs +class SpherHarMath +{ + public: + static constexpr int MaxSupportedL = 5; // analytic Ylm implemented up to l=5 + 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| + + /// Values of various coefficients + void initializeYlms() + { + double oneoversqrtpi = 1.0 / std::sqrt(o2::constants::math::PI); + + // l=0 prefactors + fgPrefactors[0] = 0.5 * oneoversqrtpi; + + // l=1 prefactors + fgPrefactors[1] = 0.5 * std::sqrt(3.0 / 2.0) * oneoversqrtpi; + fgPrefactors[2] = 0.5 * std::sqrt(3.0) * oneoversqrtpi; + fgPrefactors[3] = -fgPrefactors[1]; + + // l=2 prefactors + fgPrefactors[4] = 0.25 * std::sqrt(15.0 / 2.0) * oneoversqrtpi; + fgPrefactors[5] = 0.5 * std::sqrt(15.0 / 2.0) * oneoversqrtpi; + fgPrefactors[6] = 0.25 * std::sqrt(5.0) * oneoversqrtpi; + fgPrefactors[7] = -fgPrefactors[5]; + fgPrefactors[8] = fgPrefactors[4]; + + // l=3 prefactors + fgPrefactors[9] = 0.125 * std::sqrt(35.0) * oneoversqrtpi; + fgPrefactors[10] = 0.25 * std::sqrt(105.0 / 2.0) * oneoversqrtpi; + fgPrefactors[11] = 0.125 * std::sqrt(21.0) * oneoversqrtpi; + fgPrefactors[12] = 0.25 * std::sqrt(7.0) * oneoversqrtpi; + fgPrefactors[13] = -fgPrefactors[11]; + fgPrefactors[14] = fgPrefactors[10]; + fgPrefactors[15] = -fgPrefactors[9]; + + // l=4 prefactors + fgPrefactors[16] = 3.0 / 16.0 * std::sqrt(35.0 / 2.0) * oneoversqrtpi; + fgPrefactors[17] = 3.0 / 8.0 * std::sqrt(35.0) * oneoversqrtpi; + fgPrefactors[18] = 3.0 / 8.0 * std::sqrt(5.0 / 2.0) * oneoversqrtpi; + fgPrefactors[19] = 3.0 / 8.0 * std::sqrt(5.0) * oneoversqrtpi; + fgPrefactors[20] = 3.0 / 16.0 * oneoversqrtpi; + fgPrefactors[21] = -fgPrefactors[19]; + fgPrefactors[22] = fgPrefactors[18]; + fgPrefactors[23] = -fgPrefactors[17]; + fgPrefactors[24] = fgPrefactors[16]; + + // l=5 prefactors + fgPrefactors[25] = 3.0 / 32.0 * std::sqrt(77.0) * oneoversqrtpi; + fgPrefactors[26] = 3.0 / 16.0 * std::sqrt(385.0 / 2.0) * oneoversqrtpi; + fgPrefactors[27] = 1.0 / 32.0 * std::sqrt(385.0) * oneoversqrtpi; + fgPrefactors[28] = 1.0 / 8.0 * std::sqrt(1155.0 / 2.0) * oneoversqrtpi; + fgPrefactors[29] = 1.0 / 16.0 * std::sqrt(165.0 / 2.0) * oneoversqrtpi; + fgPrefactors[30] = 1.0 / 16.0 * std::sqrt(11.0) * oneoversqrtpi; + fgPrefactors[31] = -fgPrefactors[29]; + fgPrefactors[32] = fgPrefactors[28]; + fgPrefactors[33] = -fgPrefactors[27]; + fgPrefactors[34] = fgPrefactors[26]; + fgPrefactors[35] = -fgPrefactors[25]; + + fgPrefshift[0] = 0; + fgPrefshift[1] = 2; + fgPrefshift[2] = 6; + fgPrefshift[3] = 12; + fgPrefshift[4] = 20; + fgPrefshift[5] = 30; + + fgPlmshift[0] = 0; + fgPlmshift[1] = 2; + fgPlmshift[2] = 5; + fgPlmshift[3] = 9; + fgPlmshift[4] = 14; + fgPlmshift[5] = 20; + } + + /// Function to calculate Legendre Polynomials + /// \param lmax Maximum value of L component + /// \param ctheta Value of theta + /// \param lbuf values of coefficients + void legendreUpToYlm(int lmax, double ctheta, std::array& lbuf) + { + // Calculate a set of legendre polynomials up to a given l + // with spherical input + std::array sins{}; + std::array coss{}; + sins[0] = 0.0; + coss[0] = 1.0; + sins[1] = std::sqrt(1 - ctheta * ctheta); + coss[1] = ctheta; + for (int iter = 2; iter < TrigCacheSize; iter++) { + sins[iter] = sins[iter - 1] * sins[1]; + coss[iter] = coss[iter - 1] * coss[1]; + } + + // Legendre polynomials l=0 + lbuf[0] = 1.0; + + // Legendre polynomials l=1 + if (lmax > 0) { + lbuf[1] = sins[1]; + lbuf[2] = coss[1]; + } + + // Legendre polynomials l=2 + if (lmax > 1) { + lbuf[3] = sins[2]; + lbuf[4] = sins[1] * coss[1]; + lbuf[5] = 3 * coss[2] - 1; + } + + // Legendre polynomials l=3 + if (lmax > 2) { // o2-linter: disable=magic-number (l index, mirrors lmax>0/1 above) + lbuf[6] = sins[3]; + lbuf[7] = sins[2] * coss[1]; + lbuf[8] = (5 * coss[2] - 1) * sins[1]; + lbuf[9] = 5 * coss[3] - 3 * coss[1]; + } + + // Legendre polynomials l=4 + if (lmax > 3) { // o2-linter: disable=magic-number (l index) + lbuf[10] = sins[4]; + lbuf[11] = sins[3] * coss[1]; + lbuf[12] = (7 * coss[2] - 1) * sins[2]; + lbuf[13] = (7 * coss[3] - 3 * coss[1]) * sins[1]; + lbuf[14] = 35 * coss[4] - 30 * coss[2] + 3; + } + + // Legendre polynomials l=5 + if (lmax > 4) { // o2-linter: disable=magic-number (l index) + lbuf[15] = sins[5]; + lbuf[16] = sins[4] * coss[1]; + lbuf[17] = (9 * coss[2] - 1) * sins[3]; + lbuf[18] = (3 * coss[3] - 1 * coss[1]) * sins[2]; + lbuf[19] = (21 * coss[4] - 14 * coss[2] + 1) * sins[1]; + lbuf[20] = 63 * coss[5] - 70 * coss[3] + 15 * coss[1]; + } + } + + /// Function to calculate a set of Ylms up to a given l with cartesian input + void doYlmUpToL(int lmax, double x, double y, double z, std::complex* ylms) + { + double ctheta = 0.0; + double phi = 0.0; + + double r = std::sqrt(x * x + y * y + z * z); + if (r < SmallLength || std::fabs(z) < SmallLength) { + ctheta = 0.0; + } else { + ctheta = z / r; + } + phi = std::atan2(y, x); + doYlmUpToL(lmax, ctheta, phi, ylms); + } + + /// Function to calculate a set of Ylms up to a given l with spherical input + void doYlmUpToL(int lmax, double ctheta, double phi, std::complex* ylms) + { + int lcur = 0; + double lpol = 0.0; + + std::array coss{}; + std::array sins{}; + + std::array lbuf{}; + legendreUpToYlm(lmax, ctheta, lbuf); + initializeYlms(); + + for (int iter = 1; iter <= lmax; iter++) { + coss[iter - 1] = std::cos(iter * phi); + sins[iter - 1] = std::sin(iter * phi); + } + + ylms[lcur++] = fgPrefactors[0] * lbuf[0] * std::complex(1, 0); + + for (int il = 1; il <= lmax; il++) { + // First im = 0 + ylms[lcur + il] = fgPrefactors[fgPrefshift[il]] * lbuf[static_cast(fgPlmshift[il])] * std::complex(1.0, 0.0); + // Im != 0 + for (int im = 1; im <= il; im++) { + lpol = lbuf[static_cast(fgPlmshift[il]) - im]; + ylms[lcur + il - im] = fgPrefactors[fgPrefshift[il] - im] * lpol * std::complex(coss[im - 1], -sins[im - 1]); + ylms[lcur + il + im] = fgPrefactors[fgPrefshift[il] + im] * lpol * std::complex(coss[im - 1], sins[im - 1]); + } + lcur += 2 * il + 1; + } + } + + private: + std::array fgPrefactors{}; + std::array fgPrefshift{}; + std::array fgPlmshift{}; +}; + +} // namespace o2::analysis::femto + +#endif // PWGCF_FEMTO_CORE_FEMTOSPHERHARMATH_H_ diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 3750f4f0e8c..a513a158f66 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/femtoSpherHarMath.h" #include "PWGCF/Femto/Core/femtoUtils.h" #include "PWGCF/Femto/Core/histManager.h" #include "PWGCF/Femto/Core/modes.h" @@ -30,12 +31,16 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include #include #include #include +#include +#include #include #include +#include #include #include #include @@ -209,6 +214,12 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup { o2::framework::ConfigurableAxis qout{"qout", {{300, -1.5f, 1.5f}}, "q_{out} (GeV/c) in LCMS"}; o2::framework::ConfigurableAxis qside{"qside", {{300, -1.5f, 1.5f}}, "q_{side} (GeV/c) in LCMS"}; o2::framework::ConfigurableAxis qlong{"qlong", {{300, -1.5f, 1.5f}}, "q_{long} (GeV/c) in LCMS"}; + o2::framework::Configurable plotSH{"plotSH", false, "(Reco) Enable spherical-harmonics decomposition of the pair momentum-difference vector"}; + o2::framework::Configurable shLMax{"shLMax", 2, "Maximum l for SH decomposition (0..5). FemtoUniverse hard-codes 1."}; + o2::framework::Configurable shFrame{"shFrame", 1, "SH reference frame/variable: 0=LCMS non-identical (k*), 1=LCMS identical (qinv, FemtoUniverse default), 2=PRF (q_PRF, matches FemtoUniverse isIdenPRF=true)"}; + o2::framework::ConfigurableAxis shKstar{"shKstar", {{60, 0.0f, 0.3f}}, "k*/qinv binning for SH histograms"}; + o2::framework::ConfigurableAxis shCentBins{"shCentBins", {o2::framework::VARIABLE_WIDTH, 0.0f, 200.0f}, "SH: multiplicity/centrality bin edges (like FemtoUniverse confMultKstarBins)"}; + o2::framework::ConfigurableAxis shKtBins{"shKtBins", {o2::framework::VARIABLE_WIDTH, 0.1f, 0.2f, 0.3f, 0.4f}, "SH: kT bin edges (like FemtoUniverse confKtKstarBins)"}; }; struct ConfPairCuts : o2::framework::ConfigurableGroup { @@ -528,6 +539,17 @@ class PairHistManager mPlotDeltaEtaDeltaPhi = ConfPairBinning.plotDeltaEtaDeltaPhi.value; mPlotBertschPratt = ConfPairBinning.plotBertschPratt.value; + mPlotSH = ConfPairBinning.plotSH.value; + mShLMax = ConfPairBinning.shLMax.value; + mShFrame = ConfPairBinning.shFrame.value; + if (mPlotSH) { + mShKstarSpec = {ConfPairBinning.shKstar, "k* (GeV/#it{c})"}; + mYlm.initializeYlms(); + // copy bin edges, stripping the leading VARIABLE_WIDTH (0) marker + mShCentEdges.assign(ConfPairBinning.shCentBins.value.begin() + 1, ConfPairBinning.shCentBins.value.end()); + mShKtEdges.assign(ConfPairBinning.shKtBins.value.begin() + 1, ConfPairBinning.shKtBins.value.end()); + } + // transverse mass type mMtType = static_cast(ConfPairBinning.transverseMassType.value); @@ -637,6 +659,10 @@ class PairHistManager std::tie(mQout, mQside, mQlong) = computeBertschPrattLCMS(mParticle1, mParticle2); } + if (mPlotSH) { + std::tie(mShKv, mShOut, mShSide, mShLong) = computeShKinematics(mParticle1, mParticle2); + } + if (mPlotDeltaEtaDeltaPhi) { mDeltaEta = particle1.eta() - particle2.eta(); mDeltaPhi = RecoDecay::constrainAngle(particle1.phi() - particle2.phi(), -o2::constants::math::PIHalf); @@ -941,6 +967,65 @@ class PairHistManager mHistogramRegistry->add(analysisDir + getHistNameV2(kQlong, HistTable), getHistDesc(kQlong, HistTable), getHistType(kQlong, HistTable), {Specs.at(kQlong)}); mHistogramRegistry->add(analysisDir + getHistNameV2(kQoutQsideQlong, HistTable), getHistDesc(kQoutQsideQlong, HistTable), getHistType(kQoutQsideQlong, HistTable), {Specs.at(kQoutQsideQlong)}); } + if (mPlotSH) { + const int nJM = (mShLMax + 1) * (mShLMax + 1); + const int nCent = static_cast(mShCentEdges.size()) - 1; // n edges -> n - 1 bins + const int nKt = static_cast(mShKtEdges.size()) - 1; + mShYlmBuffer.assign(nJM, {}); + + mShReal.resize(nCent); + mShImag.resize(nCent); + for (int iCent = 0; iCent < nCent; ++iCent) { + mShReal[iCent].resize(nKt); + mShImag[iCent].resize(nKt); + // folder name: mult_{low}_{high} + const std::string centFolder = "mult_" + std::to_string(static_cast(mShCentEdges[iCent])) + + "_" + std::to_string(static_cast(mShCentEdges[iCent + 1])); + for (int iKt = 0; iKt < nKt; ++iKt) { + mShReal[iCent][iKt].resize(nJM); + mShImag[iCent][iKt].resize(nJM); + // folder name: kT_{low*100}_{high*100} + std::string ktFolder = "kT_"; + ktFolder += std::to_string(static_cast(mShKtEdges[iKt] * 100.0)); + ktFolder += "_"; + ktFolder += std::to_string(static_cast(mShKtEdges[iKt + 1] * 100.0)); + std::string dir = std::string(prefix) + std::string(AnalysisDir) + "SH/"; + dir += centFolder; + dir += "/"; + dir += ktFolder; + dir += "/"; + + int ihist = 0; + for (int l = 0; l <= mShLMax; ++l) { + for (int m = -l; m <= l; ++m) { + std::string lm = std::to_string(l); + lm += (m < 0) ? std::to_string(l - m) : std::to_string(m); + std::string nameRe = dir; + nameRe += "ReYlm"; + nameRe += lm; + std::string nameIm = dir; + nameIm += "ImYlm"; + nameIm += lm; + // shared "Y_{l}^{m}" suffix for both titles + std::string ylmLabel = "Y_{"; + ylmLabel += std::to_string(l); + ylmLabel += "}^{"; + ylmLabel += std::to_string(m); + ylmLabel += "}"; + std::string titleRe = "Re "; + titleRe += ylmLabel; + titleRe += "; k* (GeV/#it{c}); Re[A_{l}^{m}]"; + std::string titleIm = "Im "; + titleIm += ylmLabel; + titleIm += "; k* (GeV/#it{c}); Im[A_{l}^{m}]"; + mShReal[iCent][iKt][ihist] = mHistogramRegistry->add(nameRe.c_str(), titleRe.c_str(), o2::framework::kTH1D, {mShKstarSpec}); + mShImag[iCent][iKt][ihist] = mHistogramRegistry->add(nameIm.c_str(), titleIm.c_str(), o2::framework::kTH1D, {mShKstarSpec}); + ++ihist; + } + } + } + } + } } // reco-vs-truth correlation histograms (kReco and kMc both set) @@ -1120,6 +1205,17 @@ class PairHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQlong, HistTable)), mQlong); mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQoutQsideQlong, HistTable)), mQout, mQside, mQlong); } + if (mPlotSH) { + const int iCent = findShBin(mMult, mShCentEdges); + const int iKt = 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) { + mShReal[iCent][iKt][i]->Fill(mShKv, std::real(mShYlmBuffer[i])); + mShImag[iCent][iKt][i]->Fill(mShKv, -std::imag(mShYlmBuffer[i])); + } + } + } } // reco-vs-truth correlation fill (kReco and kMc both set) @@ -1283,6 +1379,97 @@ class PairHistManager 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)}; + } + o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; bool mUsePdgMass = true; double mPdgMass1 = 0.; @@ -1374,6 +1561,31 @@ class PairHistManager float mDeltaEta = 0.f; float mDeltaPhi = 0.f; + // Spherical harmonics + bool mPlotSH = 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() + + // kinematics computed in setPair(): axis value + 3 components feeding Ylm + float mShKv = 0.f; // kstar (non-identical) or qinv (identical) + float mShOut = 0.f; + float mShSide = 0.f; + float mShLong = 0.f; + + // SH histograms binned in [iCent][iKt][ihist]; ihist = l*(l+1)+m + std::vector>>> mShReal; + std::vector>>> mShImag; + std::vector> mShYlmBuffer; // reused, allocated once + std::vector mShCentEdges; + std::vector mShKtEdges; + + o2::analysis::femto::SpherHarMath mYlm{}; + // qa bool mPairCorrelationQa = false; bool mEventMixingQa = false;