Skip to content

[PWGCF] Add mc3d femtodream #15160

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
wenyaCern wants to merge 4 commits into
base: master
Choose a base branch
from
+312 −43
Open

[PWGCF] Add mc3d femtodream #15160

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
9b318ea
add MC analysis for 3D femto in femtoDream
wenyaCern Feb 25, 2026
d0bdb89
add 3Dfemto for MC
wenyaCern Feb 25, 2026
76f6d91
Please consider the following formatting changes
alibuild Feb 25, 2026
4a34d08
Merge pull request #16 from alibuild/alibot-cleanup-15160
wenyaCern Feb 25, 2026
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
56 changes: 49 additions & 7 deletions PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,6 +223,10 @@ class FemtoDreamCollisionSelection
mHistogramQn->add("Event/centVsqnVsSpher", "; cent; qn; Sphericity", kTH3F, {{10, 0, 100}, {100, 0, 1000}, {100, 0, 1}});
mHistogramQn->add("Event/qnBin", "; qnBin; entries", kTH1F, {{20, 0, 20}});
mHistogramQn->add("Event/psiEP", "; #Psi_{EP} (deg); entries", kTH1F, {{100, 0, 180}});
mHistogramQn->add("Event/epReso_FT0CTPC", "; cent; qnBin; reso_ft0c_tpc", kTH2F, {{10, 0, 100}, {10, 0, 10}});
mHistogramQn->add("Event/epReso_FT0ATPC", "; cent; qnBin; reso_ft0a_tpc", kTH2F, {{10, 0, 100}, {10, 0, 10}});
mHistogramQn->add("Event/epReso_FT0CFT0A", "; cent; qnBin; reso_ft0c_ft0a", kTH2F, {{10, 0, 100}, {10, 0, 10}});
mHistogramQn->add("Event/epReso_count", "; cent; qnBin; count", kTH2F, {{10, 0, 100}, {10, 0, 10}});

return;
}
Expand Down Expand Up @@ -330,9 +334,17 @@ class FemtoDreamCollisionSelection
/// \param col Collision
/// \return value of the qn-vector of FT0C of the event
template <typename T>
float computeqnVec(T const& col)
float computeqnVec(T const& col, int qvecMod = 0)
{
double qn = std::sqrt(col.qvecFT0CReVec()[0] * col.qvecFT0CReVec()[0] + col.qvecFT0CImVec()[0] * col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
double qn = -999.f;
if (qvecMod == 0) {
qn = std::sqrt(col.qvecFT0CReVec()[0] * col.qvecFT0CReVec()[0] + col.qvecFT0CImVec()[0] * col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
} else if (qvecMod == 1) {
qn = std::sqrt(col.qvecFT0AReVec()[0] * col.qvecFT0AReVec()[0] + col.qvecFT0AImVec()[0] * col.qvecFT0AImVec()[0]) * std::sqrt(col.sumAmplFT0A());
} else {
LOGP(error, "no selected detector of Qvec for ESE ");
return qn;
}
return qn;
}

Expand All @@ -342,15 +354,43 @@ class FemtoDreamCollisionSelection
/// \param nmode EP in which harmonic(default 2nd harmonic)
/// \return angle of the event plane (rad) of FT0C of the event
template <typename T>
float computeEP(T const& col, int nmode)
float computeEP(T const& col, int nmode, int qvecMod)
{
double EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
if (EP < 0)
double EP = -999.f;
if (qvecMod == 0) {
EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
} else if (qvecMod == 1) {
EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0AImVec()[0], col.qvecFT0AReVec()[0])));
} else if (qvecMod == 2) {
EP = ((1. / nmode) * (TMath::ATan2(col.qvecTPCallImVec()[0], col.qvecTPCallReVec()[0])));
} else {
LOGP(error, "no selected detector of Qvec for EP");
return EP;
}

if (EP < 0) {
EP += TMath::Pi();
// atan2 return in rad -pi/2-pi/2, then make it 0-pi
} // atan2 return in rad -pi/2-pi/2, then make it 0-pi
return EP;
}

/// Compute the event plane resolution of 3 sub-events
/// \tparam T type of the collision
/// \param col Collision
/// \param nmode EP in which harmonic(default 2nd harmonic)
template <typename T>
void fillEPReso(T const& col, int nmode, float centrality)
{
const float psi_ft0c = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
const float psi_ft0a = ((1. / nmode) * (TMath::ATan2(col.qvecFT0AImVec()[0], col.qvecFT0AReVec()[0])));
const float psi_tpc = ((1. / nmode) * (TMath::ATan2(col.qvecTPCallImVec()[0], col.qvecTPCallReVec()[0])));

mHistogramQn->fill(HIST("Event/epReso_FT0CTPC"), centrality, mQnBin + 0.f, std::cos((psi_ft0c - psi_tpc) * nmode));
mHistogramQn->fill(HIST("Event/epReso_FT0ATPC"), centrality, mQnBin + 0.f, std::cos((psi_ft0a - psi_tpc) * nmode));
mHistogramQn->fill(HIST("Event/epReso_FT0CFT0A"), centrality, mQnBin + 0.f, std::cos((psi_ft0c - psi_ft0a) * nmode));
mHistogramQn->fill(HIST("Event/epReso_count"), centrality, mQnBin + 0.f);
}

/// \return the 1-d qn-vector separator to 2-d
std::vector<std::vector<float>> getQnBinSeparator2D(std::vector<float> flat, const int numQnBins = 10)
{
Expand Down Expand Up @@ -413,13 +453,15 @@ class FemtoDreamCollisionSelection
}

/// \fill event-wise informations
void fillEPQA(float centrality, float fSpher, float qn, float psiEP)
template <typename T>
void fillEPQA(T& col, float centrality, float fSpher, float qn, float psiEP, int nmode = 2)
{
mHistogramQn->fill(HIST("Event/centFT0CBeforeQn"), centrality);
mHistogramQn->fill(HIST("Event/centVsqn"), centrality, qn);
mHistogramQn->fill(HIST("Event/centVsqnVsSpher"), centrality, qn, fSpher);
mHistogramQn->fill(HIST("Event/qnBin"), mQnBin + 0.f);
mHistogramQn->fill(HIST("Event/psiEP"), psiEP);
fillEPReso(col, nmode, centrality);
}

/// \todo to be implemented!
Expand Down
57 changes: 52 additions & 5 deletions PWGCF/FemtoDream/Core/femtoDreamContainer.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -223,10 +223,26 @@ class FemtoDreamContainer
mHistogramRegistry->add((folderName + "/relPair3dRmTMultPercentileQnPairphi").c_str(), ("; " + femtoDKout + femtoDKside + femtoDKlong + "; #it{m}_{T} (GeV/#it{c}); Centrality; qn; #varphi_{pair} - #Psi_{EP}").c_str(), kTHnSparseF, {femtoDKoutAxis, femtoDKsideAxis, femtoDKlongAxis, mTAxi4D, multPercentileAxis4D, qnAxis, pairPhiAxis});
}

template <typename T>
void init_3Dqn_MC(std::string folderName, std::string femtoDKout, std::string femtoDKside, std::string femtoDKlong,
T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, bool smearingByOrigin = false)
{
mHistogramRegistry->add((folderName + "/hNoMCtruthPairsCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
mHistogramRegistry->add((folderName + "/hFakePairsCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
if (smearingByOrigin) {
const int nOriginBins = o2::aod::femtodreamMCparticle::ParticleOriginMCTruth::kNOriginMCTruthTypes;
mHistogramRegistry->add((folderName + "/relPair3dresolution").c_str(), (";" + femtoDKout + "mctruth;" + femtoDKout + "_reco;" + femtoDKside + "mctruth;" + femtoDKside + "_reco;" + femtoDKlong + "mctruth;" + femtoDKlong + "_reco;" + "MC origin particle 1; MC origin particle 2;").c_str(),
kTHnSparseF, {femtoDKoutAxis, femtoDKoutAxis, femtoDKsideAxis, femtoDKsideAxis, femtoDKlongAxis, femtoDKlongAxis, {nOriginBins, 0, nOriginBins}, {nOriginBins, 0, nOriginBins}});
} else {
mHistogramRegistry->add((folderName + "/relPair3dresolution").c_str(), (";" + femtoDKout + "mctruth;" + femtoDKside + "mctruth;" + femtoDKlong + "mctruth;" + femtoDKout + "_reco;" + femtoDKside + "_reco;" + femtoDKlong + "_reco;").c_str(),
kTHnSparseF, {femtoDKoutAxis, femtoDKoutAxis, femtoDKsideAxis, femtoDKsideAxis, femtoDKlongAxis, femtoDKlongAxis});
}
}

template <typename T>
void init_3Dqn(HistogramRegistry* registry,
T& DKoutBins, T& DKsideBins, T& DKlongBins, T& mTBins4D, T& multPercentileBins4D,
bool isMC, ConfigurableAxis qnBins = {"qnBins", {10, 0, 10}, "qn binning"}, ConfigurableAxis pairPhiBins = {"phiBins", {10, 0 - 0.05, TMath::Pi() + 0.05}, "pair phi binning"})
bool isMC, ConfigurableAxis qnBins = {"qnBins", {10, 0, 10}, "qn binning"}, ConfigurableAxis pairPhiBins = {"phiBins", {10, 0 - 0.05, TMath::Pi() + 0.05}, "pair phi binning"}, bool smearingByOrigin = false)
{
mHistogramRegistry = registry;

Expand All @@ -251,6 +267,8 @@ class FemtoDreamContainer
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);
init_3Dqn_MC(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
DKoutAxis, DKsideAxis, DKlongAxis, smearingByOrigin);
}
}

Expand Down Expand Up @@ -484,11 +502,26 @@ class FemtoDreamContainer
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
/// Fills resolution of DKout, DKside, DKlong
void setPair_3Dqn_MC(std::vector<double> k3dMC, std::vector<double> k3d, const int originPart1, const int originPart2, bool smearingByOrigin)
{
if (smearingByOrigin) {
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("_3Dqn") + HIST("/relPair3dresolution"), k3dMC[1], k3d[1], k3dMC[2], k3d[2], k3dMC[3], k3d[3], originPart1, originPart2);
} else {
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("_3Dqn") + HIST("/relPair3dresolution"), k3dMC[1], k3d[1], k3dMC[2], k3d[2], k3dMC[3], k3d[3]);
}
}

template <bool isMC, typename T1, typename T2>
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane)
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane, bool smearingByOrigin = false)
{

std::vector<double> k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies);
if (k3d.size() < 4) {
LOG(error) << "newpairfunc returned size=" << k3d.size();
return;
}
float DKout = k3d[1];
float DKside = k3d[2];
float DKlong = k3d[3];
Expand All @@ -503,12 +536,17 @@ class FemtoDreamContainer
if constexpr (isMC) {
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {

std::vector<double> k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
std::vector<double> k3dMC = FemtoDreamMath::newpairfuncMC(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
if (k3dMC.size() < 4) {
LOG(error) << "newpairfunc returned size=" << k3d.size();
return;
}
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane);

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_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 @@ -521,10 +559,14 @@ class FemtoDreamContainer
}

template <bool isMC, typename T1, typename T2>
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float EP1, const float EP2)
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float EP1, const float EP2, bool smearingByOrigin = false)
{

std::vector<double> k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies);
if (k3d.size() < 4) {
LOG(error) << "newpairfunc returned size=" << k3d.size();
return;
}
float DKout = k3d[1];
float DKside = k3d[2];
float DKlong = k3d[3];
Expand All @@ -539,12 +581,17 @@ class FemtoDreamContainer
if constexpr (isMC) {
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {

std::vector<double> k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
std::vector<double> k3dMC = FemtoDreamMath::newpairfuncMC(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
if (k3dMC.size() < 4) {
LOG(error) << "newpairfunc returned size=" << k3d.size();
return;
}
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2);

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_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
92 changes: 92 additions & 0 deletions PWGCF/FemtoDream/Core/femtoDreamMath.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -249,6 +249,98 @@ class FemtoDreamMath
return vect;
}

/// Compute the 3d components of the pair momentum in LCMS and PRF
/// Copy from femto universe
/// \tparam T type of tracks
/// \param part1 Particle 1
/// \param mass1 Mass of particle 1
/// \param part2 Particle 2
/// \param mass2 Mass of particle 2
/// \param isiden Identical or non-identical particle pair
template <typename T>
static std::vector<double> newpairfuncMC(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden)
{
const double trkPx1 = part1.pt() * std::cos(part1.phi());
const double trkPy1 = part1.pt() * std::sin(part1.phi());
const double trkPz1 = part1.pt() * std::sinh(part1.eta());

const double trkPx2 = part2.pt() * std::cos(part2.phi());
const double trkPy2 = part2.pt() * std::sin(part2.phi());
const double trkPz2 = part2.pt() * std::sinh(part2.eta());

const double e1 = std::sqrt(std::pow(trkPx1, 2) + std::pow(trkPy1, 2) + std::pow(trkPz1, 2) + std::pow(mass1, 2));
const double e2 = std::sqrt(std::pow(trkPx2, 2) + std::pow(trkPy2, 2) + std::pow(trkPz2, 2) + std::pow(mass2, 2));

const ROOT::Math::PxPyPzEVector vecpart1(trkPx1, trkPy1, trkPz1, e1);
const ROOT::Math::PxPyPzEVector vecpart2(trkPx2, trkPy2, trkPz2, e2);
const ROOT::Math::PxPyPzEVector trackSum = vecpart1 + vecpart2;

std::vector<double> vect;

const double tPx = trackSum.px();
const double tPy = trackSum.py();
const double tPz = trackSum.pz();
const double tE = trackSum.E();

const double tPtSq = (tPx * tPx + tPy * tPy);
const double tMtSq = (tE * tE - tPz * tPz);
const double tM = std::sqrt(tMtSq - tPtSq);
const double tMt = std::sqrt(tMtSq);
const double tPt = std::sqrt(tPtSq);

// Boost to LCMS

const double beta_LCMS = tPz / tE;
const double gamma_LCMS = tE / tMt;

const double fDKOut = (trkPx1 * tPx + trkPy1 * tPy) / tPt;
const double fDKSide = (-trkPx1 * tPy + trkPy1 * tPx) / tPt;
const double fDKLong = gamma_LCMS * (trkPz1 - beta_LCMS * e1);
const double fDE = gamma_LCMS * (e1 - beta_LCMS * trkPz1);

const double px1LCMS = fDKOut;
const double py1LCMS = fDKSide;
const double pz1LCMS = fDKLong;
const double pE1LCMS = fDE;

const double px2LCMS = (trkPx2 * tPx + trkPy2 * tPy) / tPt;
const double py2LCMS = (trkPy2 * tPx - trkPx2 * tPy) / tPt;
const double pz2LCMS = gamma_LCMS * (trkPz2 - beta_LCMS * e2);
const double pE2LCMS = gamma_LCMS * (e2 - beta_LCMS * trkPz2);

const double fDKOutLCMS = px1LCMS - px2LCMS;
const double fDKSideLCMS = py1LCMS - py2LCMS;
const double fDKLongLCMS = pz1LCMS - pz2LCMS;

// Boost to PRF

const double betaOut = tPt / tMt;
const double gammaOut = tMt / tM;

const double fDKOutPRF = gammaOut * (fDKOutLCMS - betaOut * (pE1LCMS - pE2LCMS));
const double fDKSidePRF = fDKSideLCMS;
const double fDKLongPRF = fDKLongLCMS;
const double fKOut = gammaOut * (fDKOut - betaOut * fDE);

const double qlcms = std::sqrt(fDKOutLCMS * fDKOutLCMS + fDKSideLCMS * fDKSideLCMS + fDKLongLCMS * fDKLongLCMS);
const double qinv = std::sqrt(fDKOutPRF * fDKOutPRF + fDKSidePRF * fDKSidePRF + fDKLongPRF * fDKLongPRF);
const double kstar = std::sqrt(fKOut * fKOut + fDKSide * fDKSide + fDKLong * fDKLong);

if (isiden) {
vect.push_back(qinv);
vect.push_back(fDKOutLCMS);
vect.push_back(fDKSideLCMS);
vect.push_back(fDKLongLCMS);
vect.push_back(qlcms);
} else {
vect.push_back(kstar);
vect.push_back(fDKOut);
vect.push_back(fDKSide);
vect.push_back(fDKLong);
}
return vect;
}

/// Compute the phi angular of a pair with respect to the event plane
/// \tparam T type of tracks
/// \param part1 Particle 1
Expand Down
Loading
Loading