Skip to content

Commit 90cf403

Browse files
ootaku303claude
ootaku303
and
claude
committed
[PWGCF] FemtoDream: add EP-rotated (B-frame) 3D q histogram
Add a storeEProt toggle to the 3D qn THnSparse in FemtoDreamContainer and femtoDreamPairTaskTrackTrack. When enabled, the out/side/long momentum-difference axes are replaced by DK_x, DK_y, DK_z rotated into the event-plane (magnetic-field) frame, so DK_y is out-of-plane (||B) for every pair. A single-component cut (e.g. on DK_x, DK_z) is applied offline on the output; this toggle only changes what is histogrammed. Co-Authored-By: Claude Sonnet 5 <noreply@anthropic.com>
1 parent bfddaae commit 90cf403

2 files changed

Lines changed: 66 additions & 10 deletions

File tree

‎PWGCF/FemtoDream/Core/femtoDreamContainer.h‎

Lines changed: 63 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -30,6 +30,7 @@
3030

3131
#include <TMath.h>
3232

33+
#include <cmath>
3334
#include <string>
3435
#include <string_view>
3536
#include <utility>
@@ -221,9 +222,14 @@ class FemtoDreamContainer
221222
/// Initialize the histograms for pairs with 3D component in divided qn bins
222223
template <typename T>
223224
void init_base_3Dqn(const std::string& folderName, const std::string& femtoDKout, const std::string& femtoDKside, const std::string& femtoDKlong,
224-
T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis)
225+
T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis, bool storeEProt)
225226
{
226-
mHistogramRegistry->add((folderName + "/relPair3dRmTMultPercentileQnPairphi").c_str(), ("; " + femtoDKout + femtoDKside + femtoDKlong + "; #it{m}_{T} (GeV/#it{c}); Centrality; qn; #varphi_{pair} - #Psi_{EP}").c_str(), o2::framework::HistType::kTHnSparseF, {femtoDKoutAxis, femtoDKsideAxis, femtoDKlongAxis, mTAxi4D, multPercentileAxis4D, qnAxis, pairPhiAxis});
227+
if (storeEProt) {
228+
// DK_x, DK_y, DK_z are the EP-rotated (B-frame) axes; DK_y is out-of-plane (||B).
229+
mHistogramRegistry->add((folderName + "/relPair3dEProtRmTMultPercentileQnPairphi").c_str(), "; DK_{x} (GeV/#it{c}); DK_{y} (GeV/#it{c}); DK_{z} (GeV/#it{c}); #it{m}_{T} (GeV/#it{c}); Centrality; qn; #varphi_{pair} - #Psi_{EP}", o2::framework::HistType::kTHnSparseF, {femtoDKoutAxis, femtoDKsideAxis, femtoDKlongAxis, mTAxi4D, multPercentileAxis4D, qnAxis, pairPhiAxis});
230+
} else {
231+
mHistogramRegistry->add((folderName + "/relPair3dRmTMultPercentileQnPairphi").c_str(), ("; " + femtoDKout + femtoDKside + femtoDKlong + "; #it{m}_{T} (GeV/#it{c}); Centrality; qn; #varphi_{pair} - #Psi_{EP}").c_str(), o2::framework::HistType::kTHnSparseF, {femtoDKoutAxis, femtoDKsideAxis, femtoDKlongAxis, mTAxi4D, multPercentileAxis4D, qnAxis, pairPhiAxis});
232+
}
227233
}
228234

229235
template <typename T>
@@ -261,15 +267,23 @@ class FemtoDreamContainer
261267
framework::AxisSpec qnAxis = {std::move(qnBins), "qn"};
262268
framework::AxisSpec pairPhiAxis = {std::move(pairPhiBins), "#varphi_{pair} - #Psi_{EP} (rad)"};
263269

270+
// EP-rotated (B-frame) axis labels, same binning as out/side/long.
271+
framework::AxisSpec DKxAxis = {DKoutBins, "DK_{x} (GeV/#it{c})"};
272+
framework::AxisSpec DKyAxis = {DKsideBins, "DK_{y} (GeV/#it{c})"};
273+
framework::AxisSpec DKzAxis = {DKlongBins, "DK_{z} (GeV/#it{c})"};
274+
framework::AxisSpec& mom1Axis = mStoreEProt ? DKxAxis : DKoutAxis;
275+
framework::AxisSpec& mom2Axis = mStoreEProt ? DKyAxis : DKsideAxis;
276+
framework::AxisSpec& mom3Axis = mStoreEProt ? DKzAxis : DKlongAxis;
277+
264278
std::string folderName = static_cast<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast<std::string>("_3Dqn");
265279

266280
init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
267-
DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis);
281+
mom1Axis, mom2Axis, mom3Axis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis, mStoreEProt);
268282

269283
if (isMC) {
270284
folderName = static_cast<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast<std::string>("_3Dqn");
271285
init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
272-
DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis);
286+
mom1Axis, mom2Axis, mom3Axis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis, mStoreEProt);
273287
init_3Dqn_MC(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
274288
DKoutAxis, DKsideAxis, DKlongAxis, smearingByOrigin);
275289
}
@@ -286,6 +300,9 @@ class FemtoDreamContainer
286300
mPDGTwo = pdg2;
287301
}
288302

303+
/// Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long. Call before init_3Dqn().
304+
void setStoreEProt(bool doStore) { mStoreEProt = doStore; }
305+
289306
/// Pass a pair to the container and compute all the relevant observables
290307
/// Called by setPair both in case of data/ and Monte Carlo reconstructed and for Monte Carlo truth
291308
/// \tparam T type of the femtodreamparticle
@@ -501,11 +518,42 @@ class FemtoDreamContainer
501518
}
502519
}
503520

521+
/// Signed φ_pair − Ψ_EP (same as FemtoDreamMath::getPairPhiEP but without the final |·|); used as the B-frame rotation angle.
522+
template <typename T1, typename T2>
523+
static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep)
524+
{
525+
const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
526+
const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
527+
const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2;
528+
return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep);
529+
}
530+
531+
/// Signed counterpart of the plane-calibrated (two-event-plane) FemtoDreamMath::getPairPhiEP.
532+
template <typename T1, typename T2>
533+
static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep1, const float Psi_ep2)
534+
{
535+
const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
536+
const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
537+
const float psidiff = Psi_ep2 - Psi_ep1;
538+
const float newPhi2 = TVector2::Phi_mpi_pi(vecpart2.Phi() - psidiff);
539+
const ROOT::Math::PtEtaPhiMVector vecpart2_calibd(vecpart2.Pt(), vecpart2.Eta(), newPhi2, vecpart2.M());
540+
const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2_calibd;
541+
return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep1);
542+
}
543+
504544
/// Pass a pair to the container and compute all the relevant observables in divided qn bins
505545
template <o2::aod::femtodreamMCparticle::MCType mc>
506-
void setPair_3Dqn_base(const float femtoDKout, const float femtoDKside, const float femtoDKlong, const float mT, const float multPercentile, const float myQnBin, const float pairPhiEP)
546+
void setPair_3Dqn_base(const float femtoDKout, const float femtoDKside, const float femtoDKlong, const float mT, const float multPercentile, const float myQnBin, const float pairPhiEP, bool storeEProt, const float pairPhiEPforRot)
507547
{
508-
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP);
548+
if (storeEProt) {
549+
// Rotate (out, side) by the signed φ_pair − Ψ_EP so DK_y is out-of-plane (||B); DK_z unchanged.
550+
const float DKx = femtoDKout * std::cos(pairPhiEPforRot) - femtoDKside * std::sin(pairPhiEPforRot);
551+
const float DKy = femtoDKout * std::sin(pairPhiEPforRot) + femtoDKside * std::cos(pairPhiEPforRot);
552+
const float DKz = femtoDKlong;
553+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dEProtRmTMultPercentileQnPairphi"), DKx, DKy, DKz, mT, multPercentile, myQnBin, pairPhiEP);
554+
} else {
555+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP);
556+
}
509557
}
510558

511559
/// Called by setPair_3Dqn only in case of Monte Carlo truth
@@ -535,9 +583,10 @@ class FemtoDreamContainer
535583
const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo);
536584

537585
const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, eventPlane);
586+
const float pairPhiEPforRot = mStoreEProt ? getPairPhiEPSigned(part1, mMassOne, part2, mMassTwo, eventPlane) : 0.f;
538587

539588
if (mHistogramRegistry) {
540-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
589+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot);
541590

542591
if constexpr (isMC) {
543592
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
@@ -549,9 +598,10 @@ class FemtoDreamContainer
549598
}
550599
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
551600
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane);
601+
const float pairPhiEPMCforRot = mStoreEProt ? getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane) : 0.f;
552602

553603
if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates
554-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
604+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, mStoreEProt, pairPhiEPMCforRot);
555605
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
556606
} else {
557607
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
@@ -580,9 +630,10 @@ class FemtoDreamContainer
580630
const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo);
581631

582632
const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, EP1, EP2);
633+
const float pairPhiEPforRot = mStoreEProt ? getPairPhiEPSigned(part1, mMassOne, part2, mMassTwo, EP1, EP2) : 0.f;
583634

584635
if (mHistogramRegistry) {
585-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
636+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot);
586637

587638
if constexpr (isMC) {
588639
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
@@ -594,9 +645,10 @@ class FemtoDreamContainer
594645
}
595646
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
596647
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2);
648+
const float pairPhiEPMCforRot = mStoreEProt ? getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2) : 0.f;
597649

598650
if (std::abs(part1.fdMCParticle().pdgMCTruth()) == mPDGOne && std::abs(part2.fdMCParticle().pdgMCTruth()) == mPDGTwo) { // Note: all pair-histogramms are filled with MC truth information ONLY in case of non-fake candidates
599-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
651+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, mStoreEProt, pairPhiEPMCforRot);
600652
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
601653
} else {
602654
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
@@ -619,6 +671,7 @@ class FemtoDreamContainer
619671
int mPDGOne = 0; ///< PDG code of particle 1
620672
int mPDGTwo = 0; ///< PDG code of particle 2
621673
float mHighkstarCut = 6.;
674+
bool mStoreEProt = false; ///< Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long
622675
};
623676

624677
} // namespace o2::analysis::femtoDream

‎PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,7 @@ struct femtoDreamPairTaskTrackTrack {
115115
ConfigurableAxis DKlong{"DKlong", {500, -2., 2.}, "binning DKlong for the 3-D femtoscopy plot: R_long(LCMS) vs mT vs multiplicity percentile vs qnBin vs pait phi wrt EP (set <<do3DFemto>> to true)"};
116116
ConfigurableAxis qnBins{"qnBins", {10, 0, 10}, "binning of qn interval"};
117117
ConfigurableAxis pairPhiBins{"pairPhiBins", {12, 0., TMath::Pi()}, "binning of pair phi"};
118+
Configurable<bool> storeEProt{"storeEProt", false, "Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long in the 3D qn histogram"};
118119
} EPCal;
119120

120121
using FilteredCollisions = soa::Filtered<FDCollisions>;
@@ -333,6 +334,8 @@ struct femtoDreamPairTaskTrackTrack {
333334
}
334335

335336
if (EPCal.do3DFemto) {
337+
sameEventQnCont.setStoreEProt(EPCal.storeEProt);
338+
mixedEventQnCont.setStoreEProt(EPCal.storeEProt);
336339
sameEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong,
337340
Binning4D.mT, Binning4D.multPercentile, Option.IsMC, EPCal.qnBins, EPCal.pairPhiBins, Option.SmearingByOrigin);
338341
mixedEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong,

0 commit comments

Comments
 (0)