Skip to content

Commit 39e6b27

Browse files
ootaku303claude
ootaku303
and
claude
authored
[PWGCF] FemtoDream: add EP-rotated (B-frame) 3D q histogram (#18136)
Co-authored-by: Claude Sonnet 5 <noreply@anthropic.com>
1 parent b286a24 commit 39e6b27

3 files changed

Lines changed: 81 additions & 7 deletions

File tree

‎PWGCF/FemtoDream/Core/femtoDreamContainer.h‎

Lines changed: 55 additions & 7 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>
@@ -223,7 +224,12 @@ class FemtoDreamContainer
223224
void init_base_3Dqn(const std::string& folderName, const std::string& femtoDKout, const std::string& femtoDKside, const std::string& femtoDKlong,
224225
T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis)
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 (mStoreEProt) {
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);
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);
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
@@ -508,6 +525,16 @@ class FemtoDreamContainer
508525
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP);
509526
}
510527

528+
/// EP-rotated (B-frame) version: rotate (out, side) by the signed φ_pair − Ψ_EP so DK_y is out-of-plane (||B); DK_z unchanged.
529+
template <o2::aod::femtodreamMCparticle::MCType mc>
530+
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, const float pairPhiEPforRot)
531+
{
532+
const float DKx = femtoDKout * std::cos(pairPhiEPforRot) - femtoDKside * std::sin(pairPhiEPforRot);
533+
const float DKy = femtoDKout * std::sin(pairPhiEPforRot) + femtoDKside * std::cos(pairPhiEPforRot);
534+
const float DKz = femtoDKlong;
535+
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dEProtRmTMultPercentileQnPairphi"), DKx, DKy, DKz, mT, multPercentile, myQnBin, pairPhiEP);
536+
}
537+
511538
/// Called by setPair_3Dqn only in case of Monte Carlo truth
512539
/// Fills resolution of DKout, DKside, DKlong
513540
void setPair_3Dqn_MC(std::vector<double> k3dMC, std::vector<double> k3d, const int originPart1, const int originPart2, bool smearingByOrigin)
@@ -535,9 +562,14 @@ class FemtoDreamContainer
535562
const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo);
536563

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

539567
if (mHistogramRegistry) {
540-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
568+
if (mStoreEProt) {
569+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, pairPhiEPforRot);
570+
} else {
571+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
572+
}
541573

542574
if constexpr (isMC) {
543575
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
@@ -549,9 +581,14 @@ class FemtoDreamContainer
549581
}
550582
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
551583
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane);
584+
const float pairPhiEPMCforRot = mStoreEProt ? FemtoDreamMath::getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane) : 0.f;
552585

553586
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);
587+
if (mStoreEProt) {
588+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, pairPhiEPMCforRot);
589+
} else {
590+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
591+
}
555592
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
556593
} else {
557594
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
@@ -580,9 +617,14 @@ class FemtoDreamContainer
580617
const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo);
581618

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

584622
if (mHistogramRegistry) {
585-
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
623+
if (mStoreEProt) {
624+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, pairPhiEPforRot);
625+
} else {
626+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
627+
}
586628

587629
if constexpr (isMC) {
588630
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
@@ -594,9 +636,14 @@ class FemtoDreamContainer
594636
}
595637
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
596638
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2);
639+
const float pairPhiEPMCforRot = mStoreEProt ? FemtoDreamMath::getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2) : 0.f;
597640

598641
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);
642+
if (mStoreEProt) {
643+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, pairPhiEPMCforRot);
644+
} else {
645+
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
646+
}
600647
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
601648
} else {
602649
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
@@ -619,6 +666,7 @@ class FemtoDreamContainer
619666
int mPDGOne = 0; ///< PDG code of particle 1
620667
int mPDGTwo = 0; ///< PDG code of particle 2
621668
float mHighkstarCut = 6.;
669+
bool mStoreEProt = false; ///< Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long
622670
};
623671

624672
} // namespace o2::analysis::femtoDream

‎PWGCF/FemtoDream/Core/femtoDreamMath.h‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -373,6 +373,29 @@ class FemtoDreamMath
373373
phi_pair_onPsi = TMath::Abs(phi_pair_onPsi);
374374
return phi_pair_onPsi;
375375
}
376+
377+
/// Signed φ_pair − Ψ_EP (same as getPairPhiEP without the final |·|); used as the B-frame rotation angle.
378+
template <typename T1, typename T2>
379+
static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep)
380+
{
381+
const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
382+
const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
383+
const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2;
384+
return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep);
385+
}
386+
387+
/// Signed counterpart of the plane-calibrated (two-event-plane) getPairPhiEP.
388+
template <typename T1, typename T2>
389+
static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep1, const float Psi_ep2)
390+
{
391+
const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1);
392+
const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2);
393+
const float psidiff = Psi_ep2 - Psi_ep1;
394+
const float newPhi2 = TVector2::Phi_mpi_pi(vecpart2.Phi() - psidiff);
395+
const ROOT::Math::PtEtaPhiMVector vecpart2_calibd(vecpart2.Pt(), vecpart2.Eta(), newPhi2, vecpart2.M());
396+
const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2_calibd;
397+
return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep1);
398+
}
376399
};
377400

378401
} // 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)