Skip to content

[PWGCF] FemtoDream: add EP-rotated (B-frame) 3D q histogram #18136

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

Open
ootaku303 wants to merge 1 commit into
base: master
Choose a base branch
from
+66 −10
Open

[PWGCF] FemtoDream: add EP-rotated (B-frame) 3D q histogram #18136

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
73 changes: 63 additions & 10 deletions PWGCF/FemtoDream/Core/femtoDreamContainer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,6 +30,7 @@

#include <TMath.h>

#include <cmath>
#include <string>
#include <string_view>
#include <utility>
Expand Down Expand Up @@ -221,9 +222,14 @@ class FemtoDreamContainer
/// Initialize the histograms for pairs with 3D component in divided qn bins
template <typename T>
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 <typename T>
Expand Down Expand Up @@ -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<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kRecon]) + static_cast<std::string>("_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<std::string>(mFolderSuffix[mEventType]) + static_cast<std::string>(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + static_cast<std::string>("_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);
}
Expand All @@ -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
Expand Down Expand Up @@ -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 <typename T1, typename T2>
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 <typename T1, typename T2>
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 <o2::aod::femtodreamMCparticle::MCType mc>
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
Expand Down Expand Up @@ -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<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot);

if constexpr (isMC) {
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
Expand All @@ -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<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(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);
Expand Down Expand Up @@ -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<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP);
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kRecon>(DKout, DKside, DKlong, mT, multPercentile, myQnBin, pairPhiEP, mStoreEProt, pairPhiEPforRot);

if constexpr (isMC) {
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
Expand All @@ -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<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(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);
Expand All @@ -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
Expand Down
3 changes: 3 additions & 0 deletions PWGCF/FemtoDream/Tasks/femtoDreamPairTaskTrackTrack.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 <<do3DFemto>> to true)"};
ConfigurableAxis qnBins{"qnBins", {10, 0, 10}, "binning of qn interval"};
ConfigurableAxis pairPhiBins{"pairPhiBins", {12, 0., TMath::Pi()}, "binning of pair phi"};
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"};
} EPCal;

using FilteredCollisions = soa::Filtered<FDCollisions>;
Expand Down Expand Up @@ -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,
Expand Down
Loading