From 90cf4035cd4bacd8d47326206a54e9e416ce9969 Mon Sep 17 00:00:00 2001 From: Takushi Omura Date: Thu, 1 Oct 2026 19:45:46 +0900 Subject: [PATCH] [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 --- PWGCF/FemtoDream/Core/femtoDreamContainer.h | 73 ++++++++++++++++--- .../Tasks/femtoDreamPairTaskTrackTrack.cxx | 3 + 2 files changed, 66 insertions(+), 10 deletions(-) diff --git a/PWGCF/FemtoDream/Core/femtoDreamContainer.h b/PWGCF/FemtoDream/Core/femtoDreamContainer.h index 98daa9eb00f..69a307d0fd5 100644 --- a/PWGCF/FemtoDream/Core/femtoDreamContainer.h +++ b/PWGCF/FemtoDream/Core/femtoDreamContainer.h @@ -30,6 +30,7 @@ #include +#include #include #include #include @@ -221,9 +222,14 @@ class FemtoDreamContainer /// Initialize the histograms for pairs with 3D component in divided qn bins template void init_base_3Dqn(const std::string& folderName, const std::string& femtoDKout, const std::string& femtoDKside, const std::string& femtoDKlong, - T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis) + T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, T& mTAxi4D, T& multPercentileAxis4D, T& qnAxis, T& pairPhiAxis, bool storeEProt) { - 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}); + if (storeEProt) { + // DK_x, DK_y, DK_z are the EP-rotated (B-frame) axes; DK_y is out-of-plane (||B). + 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}); + } else { + 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}); + } } template @@ -261,15 +267,23 @@ class FemtoDreamContainer framework::AxisSpec qnAxis = {std::move(qnBins), "qn"}; framework::AxisSpec pairPhiAxis = {std::move(pairPhiBins), "#varphi_{pair} - #Psi_{EP} (rad)"}; + // EP-rotated (B-frame) axis labels, same binning as out/side/long. + framework::AxisSpec DKxAxis = {DKoutBins, "DK_{x} (GeV/#it{c})"}; + framework::AxisSpec DKyAxis = {DKsideBins, "DK_{y} (GeV/#it{c})"}; + framework::AxisSpec DKzAxis = {DKlongBins, "DK_{z} (GeV/#it{c})"}; + framework::AxisSpec& mom1Axis = mStoreEProt ? DKxAxis : DKoutAxis; + framework::AxisSpec& mom2Axis = mStoreEProt ? DKyAxis : DKsideAxis; + framework::AxisSpec& mom3Axis = mStoreEProt ? DKzAxis : DKlongAxis; + std::string folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast("_3Dqn"); init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong, - DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis); + mom1Axis, mom2Axis, mom3Axis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis, mStoreEProt); if (isMC) { folderName = static_cast(mFolderSuffix[mEventType]) + static_cast(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast("_3Dqn"); init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong, - DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis); + mom1Axis, mom2Axis, mom3Axis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis, mStoreEProt); init_3Dqn_MC(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong, DKoutAxis, DKsideAxis, DKlongAxis, smearingByOrigin); } @@ -286,6 +300,9 @@ class FemtoDreamContainer mPDGTwo = pdg2; } + /// Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long. Call before init_3Dqn(). + void setStoreEProt(bool doStore) { mStoreEProt = doStore; } + /// Pass a pair to the container and compute all the relevant observables /// Called by setPair both in case of data/ and Monte Carlo reconstructed and for Monte Carlo truth /// \tparam T type of the femtodreamparticle @@ -501,11 +518,42 @@ class FemtoDreamContainer } } + /// Signed φ_pair − Ψ_EP (same as FemtoDreamMath::getPairPhiEP but without the final |·|); used as the B-frame rotation angle. + template + static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep) + { + const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1); + const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2); + const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2; + return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep); + } + + /// Signed counterpart of the plane-calibrated (two-event-plane) FemtoDreamMath::getPairPhiEP. + template + static float getPairPhiEPSigned(const T1& part1, const float mass1, const T2& part2, const float mass2, const float Psi_ep1, const float Psi_ep2) + { + const ROOT::Math::PtEtaPhiMVector vecpart1(part1.pt(), part1.eta(), part1.phi(), mass1); + const ROOT::Math::PtEtaPhiMVector vecpart2(part2.pt(), part2.eta(), part2.phi(), mass2); + const float psidiff = Psi_ep2 - Psi_ep1; + const float newPhi2 = TVector2::Phi_mpi_pi(vecpart2.Phi() - psidiff); + const ROOT::Math::PtEtaPhiMVector vecpart2_calibd(vecpart2.Pt(), vecpart2.Eta(), newPhi2, vecpart2.M()); + const ROOT::Math::PtEtaPhiMVector trackSum = vecpart1 + vecpart2_calibd; + return TVector2::Phi_mpi_pi(trackSum.Phi() - Psi_ep1); + } + /// Pass a pair to the container and compute all the relevant observables in divided qn bins template - 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) + 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) { - mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP); + if (storeEProt) { + // Rotate (out, side) by the signed φ_pair − Ψ_EP so DK_y is out-of-plane (||B); DK_z unchanged. + const float DKx = femtoDKout * std::cos(pairPhiEPforRot) - femtoDKside * std::sin(pairPhiEPforRot); + const float DKy = femtoDKout * std::sin(pairPhiEPforRot) + femtoDKside * std::cos(pairPhiEPforRot); + const float DKz = femtoDKlong; + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dEProtRmTMultPercentileQnPairphi"), DKx, DKy, DKz, mT, multPercentile, myQnBin, pairPhiEP); + } else { + mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP); + } } /// Called by setPair_3Dqn only in case of Monte Carlo truth @@ -535,9 +583,10 @@ class FemtoDreamContainer const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, eventPlane); + const float pairPhiEPforRot = mStoreEProt ? getPairPhiEPSigned(part1, mMassOne, part2, mMassTwo, eventPlane) : 0.f; if (mHistogramRegistry) { - setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP); + setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot); if constexpr (isMC) { if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { @@ -549,9 +598,10 @@ class FemtoDreamContainer } const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane); + const float pairPhiEPMCforRot = mStoreEProt ? getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane) : 0.f; 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 - setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC); + setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, mStoreEProt, pairPhiEPMCforRot); setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin); } else { mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); @@ -580,9 +630,10 @@ class FemtoDreamContainer const float mT = FemtoDreamMath::getmT(part1, mMassOne, part2, mMassTwo); const float pairPhiEP = FemtoDreamMath::getPairPhiEP(part1, mMassOne, part2, mMassTwo, EP1, EP2); + const float pairPhiEPforRot = mStoreEProt ? getPairPhiEPSigned(part1, mMassOne, part2, mMassTwo, EP1, EP2) : 0.f; if (mHistogramRegistry) { - setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP); + setPair_3Dqn_base(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot); if constexpr (isMC) { if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) { @@ -594,9 +645,10 @@ class FemtoDreamContainer } const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo); const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2); + const float pairPhiEPMCforRot = mStoreEProt ? getPairPhiEPSigned(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2) : 0.f; 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 - setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC); + setPair_3Dqn_base(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC, mStoreEProt, pairPhiEPMCforRot); setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin); } else { mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0); @@ -619,6 +671,7 @@ class FemtoDreamContainer int mPDGOne = 0; ///< PDG code of particle 1 int mPDGTwo = 0; ///< PDG code of particle 2 float mHighkstarCut = 6.; + bool mStoreEProt = false; ///< Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long }; } // namespace o2::analysis::femtoDream diff --git a/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx b/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx index 22e69c80e49..0c0799f38e8 100644 --- a/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx +++ b/PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx @@ -115,6 +115,7 @@ struct femtoDreamPairTaskTrackTrack { 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 <> to true)"}; ConfigurableAxis qnBins{"qnBins", {10, 0, 10}, "binning of qn interval"}; ConfigurableAxis pairPhiBins{"pairPhiBins", {12, 0., TMath::Pi()}, "binning of pair phi"}; + Configurable storeEProt{"storeEProt", false, "Store EP-rotated (B-frame) DK_x,DK_y,DK_z instead of out,side,long in the 3D qn histogram"}; } EPCal; using FilteredCollisions = soa::Filtered; @@ -333,6 +334,8 @@ struct femtoDreamPairTaskTrackTrack { } if (EPCal.do3DFemto) { + sameEventQnCont.setStoreEProt(EPCal.storeEProt); + mixedEventQnCont.setStoreEProt(EPCal.storeEProt); sameEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong, Binning4D.mT, Binning4D.multPercentile, Option.IsMC, EPCal.qnBins, EPCal.pairPhiBins, Option.SmearingByOrigin); mixedEventQnCont.init_3Dqn(&Registry, EPCal.DKout, EPCal.DKside, EPCal.DKlong,