Skip to content

Commit 079e9ed

Browse files
wenyaCernalibuild
andauthored
[PWGCF] Add mc3d femtodream (#15160)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent ba6a03d commit 079e9ed

File tree

5 files changed

+312
-43
lines changed

5 files changed

+312
-43
lines changed

PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h

Lines changed: 49 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -223,6 +223,10 @@ class FemtoDreamCollisionSelection
223223
mHistogramQn->add("Event/centVsqnVsSpher", "; cent; qn; Sphericity", kTH3F, {{10, 0, 100}, {100, 0, 1000}, {100, 0, 1}});
224224
mHistogramQn->add("Event/qnBin", "; qnBin; entries", kTH1F, {{20, 0, 20}});
225225
mHistogramQn->add("Event/psiEP", "; #Psi_{EP} (deg); entries", kTH1F, {{100, 0, 180}});
226+
mHistogramQn->add("Event/epReso_FT0CTPC", "; cent; qnBin; reso_ft0c_tpc", kTH2F, {{10, 0, 100}, {10, 0, 10}});
227+
mHistogramQn->add("Event/epReso_FT0ATPC", "; cent; qnBin; reso_ft0a_tpc", kTH2F, {{10, 0, 100}, {10, 0, 10}});
228+
mHistogramQn->add("Event/epReso_FT0CFT0A", "; cent; qnBin; reso_ft0c_ft0a", kTH2F, {{10, 0, 100}, {10, 0, 10}});
229+
mHistogramQn->add("Event/epReso_count", "; cent; qnBin; count", kTH2F, {{10, 0, 100}, {10, 0, 10}});
226230

227231
return;
228232
}
@@ -330,9 +334,17 @@ class FemtoDreamCollisionSelection
330334
/// \param col Collision
331335
/// \return value of the qn-vector of FT0C of the event
332336
template <typename T>
333-
float computeqnVec(T const& col)
337+
float computeqnVec(T const& col, int qvecMod = 0)
334338
{
335-
double qn = std::sqrt(col.qvecFT0CReVec()[0] * col.qvecFT0CReVec()[0] + col.qvecFT0CImVec()[0] * col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
339+
double qn = -999.f;
340+
if (qvecMod == 0) {
341+
qn = std::sqrt(col.qvecFT0CReVec()[0] * col.qvecFT0CReVec()[0] + col.qvecFT0CImVec()[0] * col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
342+
} else if (qvecMod == 1) {
343+
qn = std::sqrt(col.qvecFT0AReVec()[0] * col.qvecFT0AReVec()[0] + col.qvecFT0AImVec()[0] * col.qvecFT0AImVec()[0]) * std::sqrt(col.sumAmplFT0A());
344+
} else {
345+
LOGP(error, "no selected detector of Qvec for ESE ");
346+
return qn;
347+
}
336348
return qn;
337349
}
338350

@@ -342,15 +354,43 @@ class FemtoDreamCollisionSelection
342354
/// \param nmode EP in which harmonic(default 2nd harmonic)
343355
/// \return angle of the event plane (rad) of FT0C of the event
344356
template <typename T>
345-
float computeEP(T const& col, int nmode)
357+
float computeEP(T const& col, int nmode, int qvecMod)
346358
{
347-
double EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
348-
if (EP < 0)
359+
double EP = -999.f;
360+
if (qvecMod == 0) {
361+
EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
362+
} else if (qvecMod == 1) {
363+
EP = ((1. / nmode) * (TMath::ATan2(col.qvecFT0AImVec()[0], col.qvecFT0AReVec()[0])));
364+
} else if (qvecMod == 2) {
365+
EP = ((1. / nmode) * (TMath::ATan2(col.qvecTPCallImVec()[0], col.qvecTPCallReVec()[0])));
366+
} else {
367+
LOGP(error, "no selected detector of Qvec for EP");
368+
return EP;
369+
}
370+
371+
if (EP < 0) {
349372
EP += TMath::Pi();
350-
// atan2 return in rad -pi/2-pi/2, then make it 0-pi
373+
} // atan2 return in rad -pi/2-pi/2, then make it 0-pi
351374
return EP;
352375
}
353376

377+
/// Compute the event plane resolution of 3 sub-events
378+
/// \tparam T type of the collision
379+
/// \param col Collision
380+
/// \param nmode EP in which harmonic(default 2nd harmonic)
381+
template <typename T>
382+
void fillEPReso(T const& col, int nmode, float centrality)
383+
{
384+
const float psi_ft0c = ((1. / nmode) * (TMath::ATan2(col.qvecFT0CImVec()[0], col.qvecFT0CReVec()[0])));
385+
const float psi_ft0a = ((1. / nmode) * (TMath::ATan2(col.qvecFT0AImVec()[0], col.qvecFT0AReVec()[0])));
386+
const float psi_tpc = ((1. / nmode) * (TMath::ATan2(col.qvecTPCallImVec()[0], col.qvecTPCallReVec()[0])));
387+
388+
mHistogramQn->fill(HIST("Event/epReso_FT0CTPC"), centrality, mQnBin + 0.f, std::cos((psi_ft0c - psi_tpc) * nmode));
389+
mHistogramQn->fill(HIST("Event/epReso_FT0ATPC"), centrality, mQnBin + 0.f, std::cos((psi_ft0a - psi_tpc) * nmode));
390+
mHistogramQn->fill(HIST("Event/epReso_FT0CFT0A"), centrality, mQnBin + 0.f, std::cos((psi_ft0c - psi_ft0a) * nmode));
391+
mHistogramQn->fill(HIST("Event/epReso_count"), centrality, mQnBin + 0.f);
392+
}
393+
354394
/// \return the 1-d qn-vector separator to 2-d
355395
std::vector<std::vector<float>> getQnBinSeparator2D(std::vector<float> flat, const int numQnBins = 10)
356396
{
@@ -413,13 +453,15 @@ class FemtoDreamCollisionSelection
413453
}
414454

415455
/// \fill event-wise informations
416-
void fillEPQA(float centrality, float fSpher, float qn, float psiEP)
456+
template <typename T>
457+
void fillEPQA(T& col, float centrality, float fSpher, float qn, float psiEP, int nmode = 2)
417458
{
418459
mHistogramQn->fill(HIST("Event/centFT0CBeforeQn"), centrality);
419460
mHistogramQn->fill(HIST("Event/centVsqn"), centrality, qn);
420461
mHistogramQn->fill(HIST("Event/centVsqnVsSpher"), centrality, qn, fSpher);
421462
mHistogramQn->fill(HIST("Event/qnBin"), mQnBin + 0.f);
422463
mHistogramQn->fill(HIST("Event/psiEP"), psiEP);
464+
fillEPReso(col, nmode, centrality);
423465
}
424466

425467
/// \todo to be implemented!

PWGCF/FemtoDream/Core/femtoDreamContainer.h

Lines changed: 52 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -223,10 +223,26 @@ class FemtoDreamContainer
223223
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});
224224
}
225225

226+
template <typename T>
227+
void init_3Dqn_MC(std::string folderName, std::string femtoDKout, std::string femtoDKside, std::string femtoDKlong,
228+
T& femtoDKoutAxis, T& femtoDKsideAxis, T& femtoDKlongAxis, bool smearingByOrigin = false)
229+
{
230+
mHistogramRegistry->add((folderName + "/hNoMCtruthPairsCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
231+
mHistogramRegistry->add((folderName + "/hFakePairsCounter").c_str(), "; Counter; Entries", kTH1I, {{1, 0, 1}});
232+
if (smearingByOrigin) {
233+
const int nOriginBins = o2::aod::femtodreamMCparticle::ParticleOriginMCTruth::kNOriginMCTruthTypes;
234+
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(),
235+
kTHnSparseF, {femtoDKoutAxis, femtoDKoutAxis, femtoDKsideAxis, femtoDKsideAxis, femtoDKlongAxis, femtoDKlongAxis, {nOriginBins, 0, nOriginBins}, {nOriginBins, 0, nOriginBins}});
236+
} else {
237+
mHistogramRegistry->add((folderName + "/relPair3dresolution").c_str(), (";" + femtoDKout + "mctruth;" + femtoDKside + "mctruth;" + femtoDKlong + "mctruth;" + femtoDKout + "_reco;" + femtoDKside + "_reco;" + femtoDKlong + "_reco;").c_str(),
238+
kTHnSparseF, {femtoDKoutAxis, femtoDKoutAxis, femtoDKsideAxis, femtoDKsideAxis, femtoDKlongAxis, femtoDKlongAxis});
239+
}
240+
}
241+
226242
template <typename T>
227243
void init_3Dqn(HistogramRegistry* registry,
228244
T& DKoutBins, T& DKsideBins, T& DKlongBins, T& mTBins4D, T& multPercentileBins4D,
229-
bool isMC, ConfigurableAxis qnBins = {"qnBins", {10, 0, 10}, "qn binning"}, ConfigurableAxis pairPhiBins = {"phiBins", {10, 0 - 0.05, TMath::Pi() + 0.05}, "pair phi binning"})
245+
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)
230246
{
231247
mHistogramRegistry = registry;
232248

@@ -251,6 +267,8 @@ class FemtoDreamContainer
251267
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");
252268
init_base_3Dqn(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
253269
DKoutAxis, DKsideAxis, DKlongAxis, mTAxis4D, multPercentileAxis4D, qnAxis, pairPhiAxis);
270+
init_3Dqn_MC(folderName, femtoObsDKout, femtoObsDKside, femtoObsDKlong,
271+
DKoutAxis, DKsideAxis, DKlongAxis, smearingByOrigin);
254272
}
255273
}
256274

@@ -484,11 +502,26 @@ class FemtoDreamContainer
484502
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[mc]) + HIST("_3Dqn") + HIST("/relPair3dRmTMultPercentileQnPairphi"), femtoDKout, femtoDKside, femtoDKlong, mT, multPercentile, myQnBin, pairPhiEP);
485503
}
486504

505+
/// Called by setPair_3Dqn only in case of Monte Carlo truth
506+
/// Fills resolution of DKout, DKside, DKlong
507+
void setPair_3Dqn_MC(std::vector<double> k3dMC, std::vector<double> k3d, const int originPart1, const int originPart2, bool smearingByOrigin)
508+
{
509+
if (smearingByOrigin) {
510+
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);
511+
} else {
512+
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]);
513+
}
514+
}
515+
487516
template <bool isMC, typename T1, typename T2>
488-
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane)
517+
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane, bool smearingByOrigin = false)
489518
{
490519

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

506-
std::vector<double> k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
539+
std::vector<double> k3dMC = FemtoDreamMath::newpairfuncMC(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
540+
if (k3dMC.size() < 4) {
541+
LOG(error) << "newpairfunc returned size=" << k3d.size();
542+
return;
543+
}
507544
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
508545
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, eventPlane);
509546

510547
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
511548
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
549+
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
512550
} else {
513551
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
514552
}
@@ -521,10 +559,14 @@ class FemtoDreamContainer
521559
}
522560

523561
template <bool isMC, typename T1, typename T2>
524-
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float EP1, const float EP2)
562+
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)
525563
{
526564

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

542-
std::vector<double> k3dMC = FemtoDreamMath::newpairfunc(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
584+
std::vector<double> k3dMC = FemtoDreamMath::newpairfuncMC(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, IsSameSpecies);
585+
if (k3dMC.size() < 4) {
586+
LOG(error) << "newpairfunc returned size=" << k3d.size();
587+
return;
588+
}
543589
const float mTMC = FemtoDreamMath::getmT(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo);
544590
const float pairPhiEPMC = FemtoDreamMath::getPairPhiEP(part1.fdMCParticle(), mMassOne, part2.fdMCParticle(), mMassTwo, EP1, EP2);
545591

546592
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
547593
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
594+
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
548595
} else {
549596
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
550597
}

PWGCF/FemtoDream/Core/femtoDreamMath.h

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -249,6 +249,98 @@ class FemtoDreamMath
249249
return vect;
250250
}
251251

252+
/// Compute the 3d components of the pair momentum in LCMS and PRF
253+
/// Copy from femto universe
254+
/// \tparam T type of tracks
255+
/// \param part1 Particle 1
256+
/// \param mass1 Mass of particle 1
257+
/// \param part2 Particle 2
258+
/// \param mass2 Mass of particle 2
259+
/// \param isiden Identical or non-identical particle pair
260+
template <typename T>
261+
static std::vector<double> newpairfuncMC(const T& part1, const float mass1, const T& part2, const float mass2, bool isiden)
262+
{
263+
const double trkPx1 = part1.pt() * std::cos(part1.phi());
264+
const double trkPy1 = part1.pt() * std::sin(part1.phi());
265+
const double trkPz1 = part1.pt() * std::sinh(part1.eta());
266+
267+
const double trkPx2 = part2.pt() * std::cos(part2.phi());
268+
const double trkPy2 = part2.pt() * std::sin(part2.phi());
269+
const double trkPz2 = part2.pt() * std::sinh(part2.eta());
270+
271+
const double e1 = std::sqrt(std::pow(trkPx1, 2) + std::pow(trkPy1, 2) + std::pow(trkPz1, 2) + std::pow(mass1, 2));
272+
const double e2 = std::sqrt(std::pow(trkPx2, 2) + std::pow(trkPy2, 2) + std::pow(trkPz2, 2) + std::pow(mass2, 2));
273+
274+
const ROOT::Math::PxPyPzEVector vecpart1(trkPx1, trkPy1, trkPz1, e1);
275+
const ROOT::Math::PxPyPzEVector vecpart2(trkPx2, trkPy2, trkPz2, e2);
276+
const ROOT::Math::PxPyPzEVector trackSum = vecpart1 + vecpart2;
277+
278+
std::vector<double> vect;
279+
280+
const double tPx = trackSum.px();
281+
const double tPy = trackSum.py();
282+
const double tPz = trackSum.pz();
283+
const double tE = trackSum.E();
284+
285+
const double tPtSq = (tPx * tPx + tPy * tPy);
286+
const double tMtSq = (tE * tE - tPz * tPz);
287+
const double tM = std::sqrt(tMtSq - tPtSq);
288+
const double tMt = std::sqrt(tMtSq);
289+
const double tPt = std::sqrt(tPtSq);
290+
291+
// Boost to LCMS
292+
293+
const double beta_LCMS = tPz / tE;
294+
const double gamma_LCMS = tE / tMt;
295+
296+
const double fDKOut = (trkPx1 * tPx + trkPy1 * tPy) / tPt;
297+
const double fDKSide = (-trkPx1 * tPy + trkPy1 * tPx) / tPt;
298+
const double fDKLong = gamma_LCMS * (trkPz1 - beta_LCMS * e1);
299+
const double fDE = gamma_LCMS * (e1 - beta_LCMS * trkPz1);
300+
301+
const double px1LCMS = fDKOut;
302+
const double py1LCMS = fDKSide;
303+
const double pz1LCMS = fDKLong;
304+
const double pE1LCMS = fDE;
305+
306+
const double px2LCMS = (trkPx2 * tPx + trkPy2 * tPy) / tPt;
307+
const double py2LCMS = (trkPy2 * tPx - trkPx2 * tPy) / tPt;
308+
const double pz2LCMS = gamma_LCMS * (trkPz2 - beta_LCMS * e2);
309+
const double pE2LCMS = gamma_LCMS * (e2 - beta_LCMS * trkPz2);
310+
311+
const double fDKOutLCMS = px1LCMS - px2LCMS;
312+
const double fDKSideLCMS = py1LCMS - py2LCMS;
313+
const double fDKLongLCMS = pz1LCMS - pz2LCMS;
314+
315+
// Boost to PRF
316+
317+
const double betaOut = tPt / tMt;
318+
const double gammaOut = tMt / tM;
319+
320+
const double fDKOutPRF = gammaOut * (fDKOutLCMS - betaOut * (pE1LCMS - pE2LCMS));
321+
const double fDKSidePRF = fDKSideLCMS;
322+
const double fDKLongPRF = fDKLongLCMS;
323+
const double fKOut = gammaOut * (fDKOut - betaOut * fDE);
324+
325+
const double qlcms = std::sqrt(fDKOutLCMS * fDKOutLCMS + fDKSideLCMS * fDKSideLCMS + fDKLongLCMS * fDKLongLCMS);
326+
const double qinv = std::sqrt(fDKOutPRF * fDKOutPRF + fDKSidePRF * fDKSidePRF + fDKLongPRF * fDKLongPRF);
327+
const double kstar = std::sqrt(fKOut * fKOut + fDKSide * fDKSide + fDKLong * fDKLong);
328+
329+
if (isiden) {
330+
vect.push_back(qinv);
331+
vect.push_back(fDKOutLCMS);
332+
vect.push_back(fDKSideLCMS);
333+
vect.push_back(fDKLongLCMS);
334+
vect.push_back(qlcms);
335+
} else {
336+
vect.push_back(kstar);
337+
vect.push_back(fDKOut);
338+
vect.push_back(fDKSide);
339+
vect.push_back(fDKLong);
340+
}
341+
return vect;
342+
}
343+
252344
/// Compute the phi angular of a pair with respect to the event plane
253345
/// \tparam T type of tracks
254346
/// \param part1 Particle 1

0 commit comments

Comments
 (0)