Skip to content

Commit 403b401

Browse files
ootaku303claude
andcommitted
[PWGCF] Femto: add event-plane-rotated 3D q histogram
Add an optional 6D THnSparse (q_x, q_y, q_z, mT, centrality, phi_pair - Psi_EP) with the pair relative momentum q = p1 - p2 expressed in the event-plane frame: - right-handed frame: x-axis along Psi_EP (in-plane), y-axis at Psi_EP + pi/2 (out-of-plane), z-axis along the beam; - (q_x, q_y) are obtained by rotating (q_out, q_side) by the signed pair angle phi = phi_pair - Psi_EP: q_x = q_out cos(phi) - q_side sin(phi) q_y = q_out sin(phi) + q_side cos(phi) - q_z = q_long (LCMS, unchanged). Enabled via PairBinning.plotEventShapeEpRotated (default off); binning reuses the q_out/q_side/q_long axes. The conventions are documented in femtomath::rotateToEventPlane. Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
1 parent a1fd0bb commit 403b401

2 files changed

Lines changed: 56 additions & 2 deletions

File tree

‎PWGCF/Femto/Core/femtoMath.h‎

Lines changed: 34 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -255,6 +255,40 @@ inline float getPairPhiFromEventPlaneAngle(const ROOT::Math::PtEtaPhiMVector& pa
255255
return phi;
256256
}
257257

258+
/// Rotate the transverse components of the pair relative momentum from the
259+
/// out-side frame (q_out, q_side) to the event-plane frame (q_x, q_y).
260+
///
261+
/// Conventions (q = p1 - p2, components in GeV/c):
262+
/// - out-side frame (computeBertschPrattLCMS): q_out is along the pair
263+
/// transverse momentum K_T = p1,T + p2,T (azimuth phi_pair), q_side is
264+
/// perpendicular to K_T in the transverse plane, q_long is along the beam
265+
/// in the LCMS.
266+
/// - event-plane frame: x-axis along the event-plane direction Psi_EP
267+
/// (col.eventPlaneAngle(), i.e. Psi_n of the harmonic and detector
268+
/// configured in the collision builder), z-axis along the beam,
269+
/// y = z cross x, i.e. at Psi_EP + pi/2.
270+
/// q_x: in-plane component (along Psi_EP),
271+
/// q_y: out-of-plane component (perpendicular to Psi_EP),
272+
/// q_z: identical to q_long (not rotated, not returned here).
273+
/// - The rotation angle is the signed pair angle
274+
/// phi = phi_pair - Psi_EP (getPairPhiFromEventPlaneAngle), so
275+
/// q_x = q_out cos(phi) - q_side sin(phi)
276+
/// q_y = q_out sin(phi) + q_side cos(phi)
277+
/// Since the longitudinal boost to the LCMS leaves the transverse plane
278+
/// unchanged, (q_x, q_y) equal the lab transverse components of p1 - p2
279+
/// projected onto the event-plane axes.
280+
/// - The sign of (q_x, q_y) flips with the pair ordering and, for the
281+
/// 2nd harmonic (Psi_2 defined modulo pi), with Psi_EP -> Psi_EP + pi;
282+
/// only quantities invariant under (q_x, q_y) -> (-q_x, -q_y)
283+
/// (e.g. |q_x|, |q_y|, q_x^2, q_y^2, q_x q_y) are physically meaningful.
284+
/// \return {q_x, q_y}
285+
inline std::tuple<float, float> rotateToEventPlane(float qOut, float qSide, float pairPhiFromEventPlaneAngle)
286+
{
287+
const float cosPhi = std::cos(pairPhiFromEventPlaneAngle);
288+
const float sinPhi = std::sin(pairPhiFromEventPlaneAngle);
289+
return {qOut * cosPhi - qSide * sinPhi, qOut * sinPhi + qSide * cosPhi};
290+
}
291+
258292
// triplet-only helpers
259293

260294
inline ROOT::Math::PxPyPzEVector getqij(ROOT::Math::PxPyPzEVector const& vi, ROOT::Math::PxPyPzEVector const& vj)

‎PWGCF/Femto/Core/pairHistManager.h‎

Lines changed: 22 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -172,6 +172,9 @@ enum PairHist {
172172

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

176179
kPairHistogramLast
177180
};
@@ -256,6 +259,7 @@ struct ConfPairBinning : o2::framework::ConfigurableGroup {
256259
o2::framework::ConfigurableAxis qside{"qside", {{300, -1.5f, 1.5f}}, "q_{side} (GeV/c) in LCMS"};
257260
o2::framework::ConfigurableAxis qlong{"qlong", {{300, -1.5f, 1.5f}}, "q_{long} (GeV/c) in LCMS"};
258261
o2::framework::Configurable<bool> plotEventShape{"plotEventShape", false, "(Reco/Mc) Enable 7D (q_out, q_side, q_long, mt, centrality event plane angle, qvector) histogram"};
262+
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)"};
259263
o2::framework::ConfigurableAxis eventPlaneAngle{"eventPlaneAngle", {{10, 0.f, 1.f * o2::constants::math::TwoPI}}, "event plane angle"};
260264
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"};
261265
o2::framework::Configurable<bool> plotSH{"plotSH", false, "(Reco) Enable spherical-harmonics decomposition of the pair momentum-difference vector"};
@@ -399,6 +403,7 @@ constexpr std::array<histmanager::HistInfo<PairHist>, kPairHistogramLast>
399403
{kTrueQsideVsQside, o2::framework::HistType::kTH2F, "hTrueQsideVsQside", "q_{side,True} vs q_{side}; q_{side,True} (GeV/#it{c}); q_{side} (GeV/#it{c})"},
400404
{kTrueQlongVsQlong, o2::framework::HistType::kTH2F, "hTrueQlongVsQlong", "q_{long,True} vs q_{long}; q_{long,True} (GeV/#it{c}); q_{long} (GeV/#it{c})"},
401405
{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;"},
406+
{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};"},
402407
}};
403408

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

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

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

1121+
if (mPlotEventShapeEpRotated) {
1122+
mHistogramRegistry->add(analysisDir + getHistNameV2(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), getHistDesc(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), getHistType(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle, HistTable), {Specs.at(kQxVsQyVsQzVsMtVsCentVsEventPlaneAngle)});
1123+
}
1124+
11141125
if (mPlotSH) {
11151126
const int nJM = (mShLMax + 1) * (mShLMax + 1);
11161127
const int nCent = static_cast<int>(mShCentEdges.size()) - 1; // n edges -> n - 1 bins
@@ -1337,15 +1348,18 @@ class PairHistManager
13371348
std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2);
13381349
}
13391350

1340-
if (mPlotEventShape) {
1351+
if (mPlotEventShape || mPlotEventShapeEpRotated) {
13411352
if (!mPlotBertschPratt) {
13421353
std::tie(mQout, mQside, mQlong) = femtomath::computeBertschPrattLCMS(mParticle1, mParticle2);
13431354
}
13441355
if constexpr (utils::HasEventShape<T4>) {
13451356
mPairPhiFromEventPlaneAngle = femtomath::getPairPhiFromEventPlaneAngle(mParticle1, mParticle2, col);
13461357
mQvector = col.qvec();
1358+
if (mPlotEventShapeEpRotated) {
1359+
std::tie(mQxEp, mQyEp) = femtomath::rotateToEventPlane(mQout, mQside, mPairPhiFromEventPlaneAngle);
1360+
}
13471361
} else {
1348-
LOG(fatal) << "plotEventShape is enabled but the collision table has no event-shape columns. Breaking...";
1362+
LOG(fatal) << "plotEventShape or plotEventShapeEpRotated is enabled but the collision table has no event-shape columns. Breaking...";
13491363
}
13501364
}
13511365

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

17171734
bool mPlotEventShape = false;
1735+
bool mPlotEventShapeEpRotated = false;
17181736

17191737
float mQout = 0.f;
17201738
float mQside = 0.f;
17211739
float mQlong = 0.f;
1740+
float mQxEp = 0.f;
1741+
float mQyEp = 0.f;
17221742

17231743
float mDeltaEta = 0.f;
17241744
float mDeltaPhi = 0.f;

0 commit comments

Comments
 (0)