Skip to content

[PWGCF] Femto: add event-plane-rotated 3D q histogram #18714

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ariedel-cern merged 1 commit into from
Nov 12, 2026
+56 −2
Merged

[PWGCF] Femto: add event-plane-rotated 3D q histogram #18714

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
34 changes: 34 additions & 0 deletions PWGCF/Femto/Core/femtoMath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -255,6 +255,40 @@ inline float getPairPhiFromEventPlaneAngle(const ROOT::Math::PtEtaPhiMVector& pa
return phi;
}

/// Rotate the transverse components of the pair relative momentum from the
/// out-side frame (q_out, q_side) to the event-plane frame (q_x, q_y).
///
/// Conventions (q = p1 - p2, components in GeV/c):
/// - out-side frame (computeBertschPrattLCMS): q_out is along the pair
/// transverse momentum K_T = p1,T + p2,T (azimuth phi_pair), q_side is
/// perpendicular to K_T in the transverse plane, q_long is along the beam
/// in the LCMS.
/// - event-plane frame: x-axis along the event-plane direction Psi_EP
/// (col.eventPlaneAngle(), i.e. Psi_n of the harmonic and detector
/// configured in the collision builder), z-axis along the beam,
/// y = z cross x, i.e. at Psi_EP + pi/2.
/// q_x: in-plane component (along Psi_EP),
/// q_y: out-of-plane component (perpendicular to Psi_EP),
/// q_z: identical to q_long (not rotated, not returned here).
/// - The rotation angle is the signed pair angle
/// phi = phi_pair - Psi_EP (getPairPhiFromEventPlaneAngle), so
/// q_x = q_out cos(phi) - q_side sin(phi)
/// q_y = q_out sin(phi) + q_side cos(phi)
/// Since the longitudinal boost to the LCMS leaves the transverse plane
/// unchanged, (q_x, q_y) equal the lab transverse components of p1 - p2
/// projected onto the event-plane axes.
/// - The sign of (q_x, q_y) flips with the pair ordering and, for the
/// 2nd harmonic (Psi_2 defined modulo pi), with Psi_EP -> Psi_EP + pi;
/// only quantities invariant under (q_x, q_y) -> (-q_x, -q_y)
/// (e.g. |q_x|, |q_y|, q_x^2, q_y^2, q_x q_y) are physically meaningful.
/// \return {q_x, q_y}
inline std::tuple<float, float> rotateToEventPlane(float qOut, float qSide, float pairPhiFromEventPlaneAngle)
{
const float cosPhi = std::cos(pairPhiFromEventPlaneAngle);
const float sinPhi = std::sin(pairPhiFromEventPlaneAngle);
return {qOut * cosPhi - qSide * sinPhi, qOut * sinPhi + qSide * cosPhi};
}

// triplet-only helpers

inline ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PxPyPzEVector const& vi, ROOT::Math::PxPyPzEVector const& vj)
Expand Down
24 changes: 22 additions & 2 deletions PWGCF/Femto/Core/pairHistManager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -172,6 +172,9 @@ enum PairHist {

// event shape enginerring
kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector,
// q in the event-plane frame (see femtomath::rotateToEventPlane):
// q_x along Psi_EP (in-plane), q_y perpendicular to Psi_EP (out-of-plane), q_z = q_long (beam)
kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle,

kPairHistogramLast
};
Expand Down Expand Up @@ -256,6 +259,7 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup {
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<bool> plotEventShape{"plotEventShape", false, "(Reco/Mc) Enable 7D (q_out, q_side, q_long, mt, centrality event plane angle, qvector) histogram"};
o2::framework::Configurable<bool> plotEventShapeEpRotated{"plotEventShapeEpRotated", false, "(Reco/Mc) Enable 6D (q_x, q_y, q_z, mt, centrality, event plane angle) histogram with (q_out, q_side) rotated to the event-plane frame (q_x in-plane, q_y out-of-plane)"};
o2::framework::ConfigurableAxis eventPlaneAngle{"eventPlaneAngle", {{10, 0.f, 1.f * o2::constants::math::TwoPI}}, "event plane angle"};
o2::framework::ConfigurableAxis qvector{"qvector", {{o2::framework::VARIABLE_WIDTH, 0.50f, 68.50f, 100.50f, 126.50f, 151.50f, 176.50f, 203.50f, 232.50f, 269.50f, 322.50f, 833.50f}}, "qvector"};
o2::framework::Configurable<bool> plotSH{"plotSH", false, "(Reco) Enable spherical-harmonics decomposition of the pair momentum-difference vector"};
Expand Down Expand Up @@ -399,6 +403,7 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast>
{kTrueQsideVsQside, o2::framework::HistType::kTH2F, "hTrueQsideVsQside", "q_{side,True} vs q_{side}; q_{side,True} (GeV/#it{c}); q_{side} (GeV/#it{c})"},
{kTrueQlongVsQlong, o2::framework::HistType::kTH2F, "hTrueQlongVsQlong", "q_{long,True} vs q_{long}; q_{long,True} (GeV/#it{c}); q_{long} (GeV/#it{c})"},
{kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, o2::framework::HistType::kTHnSparseF, "hQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector", "Event shape enginering; q_{out} (GeV/#it{c}); q_{side} (GeV/#it{c}); q_{long} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); centrality (%); #varphi_{pair} - #Psi_{EP}; q-vector;"},
{kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, o2::framework::HistType::kTHnSparseF, "hQxVsQyVsQzVsMtVsCentVsEventPlaneAngle", "Event shape enginering, q in the event-plane frame; q_{x}^{EP} (in-plane) (GeV/#it{c}); q_{y}^{EP} (out-of-plane) (GeV/#it{c}); q_{z} = q_{long} (GeV/#it{c}); m_{T} (GeV/#it{c}^{2}); centrality (%); #varphi_{pair} - #Psi_{EP};"},
}};

// NOLINTNEXTLINE(cppcoreguidelines-macro-usage)
Expand Down Expand Up @@ -461,6 +466,7 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast>
{kQside, {(confAnalysis).qside}}, \
{kQlong, {(confAnalysis).qlong}}, \
{kQoutVsQsideVsQlong, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong}}, \
{kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong, (confAnalysis).mt, (confAnalysis).centrality, (confAnalysis).eventPlaneAngle}}, \
{kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, {(confAnalysis).qout, (confAnalysis).qside, (confAnalysis).qlong, (confAnalysis).mt, (confAnalysis).centrality, (confAnalysis).eventPlaneAngle, (confAnalysis).qvector}},

// mixing-qa entries are independent of reco vs mc-truth status — both the reco
Expand Down Expand Up @@ -648,6 +654,7 @@ class PairHistManager
mPlotKstarVsBetaLabVsMtVsPt1VsPt2VsMult = ConfPairBinning.plotKstarVsBetaLabVsMtVsPt1VsPt2VsMult.value;
mPlotBertschPratt = ConfPairBinning.plotBertschPratt.value;
mPlotEventShape = ConfPairBinning.plotEventShape.value;
mPlotEventShapeEpRotated = ConfPairBinning.plotEventShapeEpRotated.value;

mPlotSH = ConfPairBinning.plotSH.value;
mShUseCent = ConfPairBinning.shUseCent.value;
Expand Down Expand Up @@ -1111,6 +1118,10 @@ class PairHistManager
mHistogramRegistry->add(analysisDir + getHistNameV2(kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, HistTable), getHistDesc(kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, HistTable), getHistType(kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, HistTable), {Specs.at(kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector)});
}

if (mPlotEventShapeEpRotated) {
mHistogramRegistry->add(analysisDir + getHistNameV2(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), getHistDesc(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), getHistType(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), {Specs.at(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle)});
}

if (mPlotSH) {
const int nJM = (mShLMax + 1) * (mShLMax + 1);
const int nCent = static_cast<int>(mShCentEdges.size()) - 1; // n edges -> n - 1 bins
Expand Down Expand Up @@ -1337,15 +1348,18 @@ class PairHistManager
std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2);
}

if (mPlotEventShape) {
if (mPlotEventShape || mPlotEventShapeEpRotated) {
if (!mPlotBertschPratt) {
std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2);
}
if constexpr (utils::HasEventShape<T4>) {
mPairPhiFromEventPlaneAngle = femtomath::getPairPhiFromEventPlaneAngle(mParticle1, mParticle2, col);
mQvector = col.qvec();
if (mPlotEventShapeEpRotated) {
std::tie(mQxEp, mQyEp) = femtomath::rotateToEventPlane(mQout, mQside, mPairPhiFromEventPlaneAngle);
}
} else {
LOG(fatal) << "plotEventShape is enabled but the collision table has no event-shape columns. Breaking...";
LOG(fatal) << "plotEventShape or plotEventShapeEpRotated is enabled but the collision table has no event-shape columns. Breaking...";
}
}

Expand Down Expand Up @@ -1512,6 +1526,9 @@ class PairHistManager
if (mPlotEventShape) {
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQoutVsQsideVsQlongVsMtVsCentVsEventPlaneAngleVsQvector, HistTable)), mQout, mQside, mQlong, mMt, mCent, mPairPhiFromEventPlaneAngle, mQvector);
}
if (mPlotEventShapeEpRotated) {
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable)), mQxEp, mQyEp, mQlong, mMt, mCent, mPairPhiFromEventPlaneAngle);
}
if (mPlotSH) {
const float shCentValue = mShUseCent ? mCent : mMult;
const int iCent = femtomath::findShBin(shCentValue, mShCentEdges);
Expand Down Expand Up @@ -1715,10 +1732,13 @@ class PairHistManager
bool mPlotBertschPratt = false;

bool mPlotEventShape = false;
bool mPlotEventShapeEpRotated = false;

float mQout = 0.f;
float mQside = 0.f;
float mQlong = 0.f;
float mQxEp = 0.f;
float mQyEp = 0.f;

float mDeltaEta = 0.f;
float mDeltaPhi = 0.f;
Expand Down
Loading