Skip to content

Commit 9b318ea

Browse files
wenyaCernwenyaCern
committed
add MC analysis for 3D femto in femtoDream
1 parent cedd873 commit 9b318ea

File tree

4 files changed

+226
-44
lines changed

4 files changed

+226
-44
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: 53 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,
236+
{nOriginBins, 0, nOriginBins}, {nOriginBins, 0, nOriginBins}});
237+
} else {
238+
mHistogramRegistry->add((folderName + "/relPair3dresolution").c_str(), (";" + femtoDKout + "mctruth;" + femtoDKside + "mctruth;" +femtoDKlong + "mctruth;" + femtoDKout + "_reco;" + femtoDKside + "_reco;" +femtoDKlong + "_reco;").c_str(),
239+
kTHnSparseF, {femtoDKoutAxis, femtoDKoutAxis, femtoDKsideAxis, femtoDKsideAxis, femtoDKlongAxis, femtoDKlongAxis}); }
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,27 @@ 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+
}
516+
487517
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)
518+
void setPair_3Dqn(T1 const& part1, T2 const& part2, const float multPercentile, bool IsSameSpecies, const float myQnBin, const float eventPlane, bool smearingByOrigin = false)
489519
{
490520

491521
std::vector<double> k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies);
522+
if (k3d.size()<4){
523+
LOG(error)<<"newpairfunc returned size="<< k3d.size();
524+
return;
525+
}
492526
float DKout = k3d[1];
493527
float DKside = k3d[2];
494528
float DKlong = k3d[3];
@@ -503,12 +537,17 @@ class FemtoDreamContainer
503537
if constexpr (isMC) {
504538
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
505539

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

510548
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
511549
setPair_3Dqn_base<o2::aod::femtodreamMCparticle::MCType::kTruth>(k3dMC[1], k3dMC[2], k3dMC[3], mTMC, multPercentile, myQnBin, pairPhiEPMC);
550+
setPair_3Dqn_MC(k3dMC, k3d, part1.fdMCParticle().partOriginMCTruth(), part2.fdMCParticle().partOriginMCTruth(), smearingByOrigin);
512551
} else {
513552
mHistogramRegistry->fill(HIST(mFolderSuffix[mEventType]) + HIST(o2::aod::femtodreamMCparticle::MCTypeName[o2::aod::femtodreamMCparticle::MCType::kTruth]) + HIST("/hFakePairsCounter"), 0);
514553
}
@@ -521,10 +560,14 @@ class FemtoDreamContainer
521560
}
522561

523562
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)
563+
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)
525564
{
526565

527566
std::vector<double> k3d = FemtoDreamMath::newpairfunc(part1, mMassOne, part2, mMassTwo, IsSameSpecies);
567+
if (k3d.size()<4){
568+
LOG(error)<<"newpairfunc returned size="<< k3d.size();
569+
return;
570+
}
528571
float DKout = k3d[1];
529572
float DKside = k3d[2];
530573
float DKlong = k3d[3];
@@ -539,12 +582,17 @@ class FemtoDreamContainer
539582
if constexpr (isMC) {
540583
if (part1.has_fdMCParticle() && part2.has_fdMCParticle()) {
541584

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

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

0 commit comments

Comments
 (0)