Skip to content

Commit 0ff30ed

Browse files
wenyaCernwenyaCern
committed
Add functions for flow and qn in femto dream
1 parent b178c96 commit 0ff30ed

File tree

5 files changed

+696
-6
lines changed

5 files changed

+696
-6
lines changed

‎PWGCF/DataModel/FemtoDerived.h‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,8 @@ DECLARE_SOA_COLUMN(BitMaskTrackTwo, bitmaskTrackTwo, BitMaskType); //! Bit f
5050
DECLARE_SOA_COLUMN(BitMaskTrackThree, bitmaskTrackThree, BitMaskType); //! Bit for track three
5151

5252
DECLARE_SOA_COLUMN(Downsample, downsample, bool); //! Flag for downsampling
53+
54+
DECLARE_SOA_COLUMN(QnBin, qnBin, int); //! qn bins for dividing events
5355
} // namespace femtodreamcollision
5456

5557
DECLARE_SOA_TABLE_STAGED(FDCollisions, "FDCOLLISION",
@@ -61,6 +63,10 @@ DECLARE_SOA_TABLE_STAGED(FDCollisions, "FDCOLLISION",
6163
femtodreamcollision::MagField);
6264
using FDCollision = FDCollisions::iterator;
6365

66+
DECLARE_SOA_TABLE(FDExtQnCollisions, "AOD", "FDEXTCOLLISION",
67+
femtodreamcollision::QnBin);
68+
using FDExtQnCollision = FDExtQnCollisions::iterator;
69+
6470
DECLARE_SOA_TABLE(FDColMasks, "AOD", "FDCOLMASK",
6571
femtodreamcollision::BitMaskTrackOne,
6672
femtodreamcollision::BitMaskTrackTwo,

‎PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h‎

Lines changed: 298 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,11 +18,15 @@
1818

1919
#include <string>
2020
#include <iostream>
21+
#include <vector>
2122
#include "Common/CCDB/TriggerAliases.h"
2223
#include "Common/DataModel/EventSelection.h"
2324
#include "Framework/HistogramRegistry.h"
2425
#include "Framework/Logger.h"
2526

27+
#include "Common/Core/EventPlaneHelper.h"
28+
#include "Common/DataModel/Qvectors.h"
29+
2630
using namespace o2::framework;
2731

2832
namespace o2::analysis::femtoDream
@@ -156,6 +160,36 @@ class FemtoDreamCollisionSelection
156160
return true;
157161
}
158162

163+
/// Pile-up selection of PbPb collisions
164+
/// \tparam T type of the collision
165+
/// \param col Collision
166+
/// \return whether or not the collisions fulfills the specified selections
167+
template <typename C>
168+
bool isPileUpCollisionPbPb(C const& col,
169+
bool noSameBunchPileup, bool isGoodZvtxFT0vsPV,
170+
bool isGoodITSLayersAll, bool noCollInRofStandard,
171+
bool noHighMultCollInPrevRof, bool noCollInTimeRangeStandard,
172+
bool /*isVertexITSTPC*/,
173+
int tpcOccupancyMin, int tpcOccupancyMax)
174+
{
175+
const auto occupancy = col.trackOccupancyInTimeRange();
176+
if ((occupancy < tpcOccupancyMin || occupancy > tpcOccupancyMax)) {
177+
return false;
178+
}
179+
if ((noSameBunchPileup && !col.selection_bit(aod::evsel::kNoSameBunchPileup))
180+
|| (isGoodZvtxFT0vsPV && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))
181+
|| (isGoodITSLayersAll && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll))
182+
|| (noCollInRofStandard && !col.selection_bit(aod::evsel::kNoCollInRofStandard))
183+
|| (noHighMultCollInPrevRof && !col.selection_bit(aod::evsel::kNoHighMultCollInPrevRof))
184+
|| (noCollInTimeRangeStandard && !col.selection_bit(aod::evsel::kNoCollInTimeRangeStandard))
185+
// || (isVertexITSTPC && !col.selection_bit(aod::evsel::kIsVertexITSTPC))
186+
) {
187+
return false;
188+
}
189+
190+
return true;
191+
}
192+
159193
/// Some basic QA of the event
160194
/// \tparam T type of the collision
161195
/// \param col Collision
@@ -175,6 +209,37 @@ class FemtoDreamCollisionSelection
175209
}
176210
}
177211

212+
/// Initializes histograms for the flow calculation
213+
/// \param registry Histogram registry to be passed
214+
void initFlow(HistogramRegistry* registry, bool doQnSeparation, int mumQnBins = 10, int binPt = 100, int binEta = 32)
215+
{
216+
if (!mCutsSet) {
217+
LOGF(error, "Event selection not set - quitting!");
218+
}
219+
mReQthisEvt = new TH2D("ReQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
220+
mImQthisEvt = new TH2D("ImQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
221+
mReQ2thisEvt = new TH2D("ReQ2thisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
222+
mImQ2thisEvt = new TH2D("ImQ2thisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
223+
mMQthisEvt = new TH2D("MQthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
224+
mMQWeightthisEvt = new TH2D("MQWeightthisEvt", "", binPt, 0., 5., binEta, -0.8, 0.8);
225+
mHistogramQn = registry;
226+
mHistogramQn->add("Event/centFT0CBefore", "; cent", kTH1F, {{10, 0, 100}});
227+
mHistogramQn->add("Event/centFT0CAfter", "; cent", kTH1F, {{10, 0, 100}});
228+
mHistogramQn->add("Event/centVsqn", "; cent; qn", kTH2F, {{10, 0, 100}, {100, 0, 1000}});
229+
mHistogramQn->add("Event/centVsqnVsSpher", "; cent; qn; Sphericity", kTH3F, {{10, 0, 100}, {100, 0, 1000}, {100, 0, 1}});
230+
mHistogramQn->add("Event/qnBin", "; qnBin; entries", kTH1F, {{20, 0, 20}});
231+
232+
mHistogramQn->add<TProfile>("Event/profileC22", "; cent; c22", kTProfile, {{10, 0, 100}});
233+
mHistogramQn->add<TProfile>("Event/profileC24", "; cent; c24", kTProfile, {{10, 0, 100}});
234+
if (doQnSeparation){
235+
for (int iqn(0); iqn < mumQnBins; ++iqn) {
236+
mHistogramQn->add<TProfile>(("Qn/profileC22_" + std::to_string(iqn)).c_str(), "; cent; c22", kTProfile, {{10, 0, 100}});
237+
mHistogramQn->add(("Qn/mult_" + std::to_string(iqn)).c_str(), "; cent; c22", kTH1F, {{100, 0, 1000}});
238+
}
239+
}
240+
return;
241+
}
242+
178243
/// \todo to be implemented!
179244
/// Compute the sphericity of an event
180245
/// Important here is that the filter on tracks does not interfere here!
@@ -246,8 +311,233 @@ class FemtoDreamCollisionSelection
246311
return spt;
247312
}
248313

314+
/// \todo to be implemented!
315+
/// Compute the qn-vector(FT0C) of an event
316+
/// \tparam T type of the collision
317+
/// \param col Collision
318+
/// \return value of the qn-vector of FT0C of the event
319+
template <typename T>
320+
float computeqnVec(T const& col)
321+
{
322+
double qn = std::sqrt(col.qvecFT0CReVec()[0] * col.qvecFT0CReVec()[0] + col.qvecFT0CImVec()[0] * col.qvecFT0CImVec()[0]) * std::sqrt(col.sumAmplFT0C());
323+
return qn;
324+
}
325+
326+
/// \todo to be implemented!
327+
/// \return the 1-d qn-vector separator to 2-d
328+
std::vector<std::vector<float>> getQnBinSeparator2D(std::vector<float> flat)
329+
{
330+
constexpr size_t nBins = 11;
331+
332+
if (flat.empty() || flat.size() % nBins != 0) {
333+
LOGP(error, "ConfQnBinSeparator size = {} is not divisible by {}",
334+
flat.size(), nBins);
335+
return {{-999, -999}};
336+
}
337+
338+
size_t nCent = flat.size() / nBins;
339+
std::vector<std::vector<float>> res(nCent, std::vector<float>(nBins));
340+
341+
for (size_t i = 0; i < nCent; ++i) {
342+
for (size_t j = 0; j < nBins; ++j) {
343+
res[i][j] = flat[i * nBins + j];
344+
}
345+
}
346+
return res;
347+
}
348+
349+
/// \todo to be implemented!
350+
/// Get the bin number of qn-vector(FT0C) of an event
351+
/// \tparam T type of the collision
352+
/// \param col Collision
353+
/// \param centBinWidth centrality bin width, example: per 1%, per 10% ...
354+
/// \return bin number of qn-vector of the event
355+
template <typename T>
356+
int myqnBin(T const& col, float centMax, std::vector<float> qnBinSeparator, float fSpher, float fMult, float centBinWidth = 1.f)
357+
{
358+
auto twoDSeparator = getQnBinSeparator2D(qnBinSeparator);
359+
if (twoDSeparator.empty() || twoDSeparator[0][0] == -999.) {
360+
LOGP(warning, "ConfQnBinSeparator not set, using default fallback!");
361+
return -999; // safe fallback
362+
}
363+
364+
mHistogramQn->fill(HIST("Event/centFT0CBefore"), col.centFT0C());
365+
int qnBin = -999;
366+
float qn = computeqnVec(col);
367+
int mycentBin = static_cast<int>(col.centFT0C() / centBinWidth);
368+
if (mycentBin >= static_cast<int>(centMax / centBinWidth))
369+
return qnBin;
370+
371+
if (mycentBin > static_cast<int>(twoDSeparator.size()) -1)
372+
return qnBin;
373+
374+
for (int iqn(0); iqn < static_cast<int>(twoDSeparator[mycentBin].size()) - 1; ++iqn) {
375+
if (qn > twoDSeparator[mycentBin][iqn] && qn <= twoDSeparator[mycentBin][iqn + 1]) {
376+
qnBin = iqn;
377+
break;
378+
} else {
379+
continue;
380+
}
381+
}
382+
383+
mHistogramQn->fill(HIST("Event/centFT0CAfter"), col.centFT0C());
384+
mHistogramQn->fill(HIST("Event/centVsqn"), col.centFT0C(), qn);
385+
mHistogramQn->fill(HIST("Event/centVsqnVsSpher"), col.centFT0C(), qn, fSpher);
386+
mHistogramQn->fill(HIST("Event/qnBin"), qnBin);
387+
if (qnBin >= 0 && qnBin < 10){
388+
switch (qnBin) {
389+
case 0:
390+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("0"), fMult);
391+
break;
392+
case 1:
393+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("1"), fMult);
394+
break;
395+
case 2:
396+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("2"), fMult);
397+
break;
398+
case 3:
399+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("3"), fMult);
400+
break;
401+
case 4:
402+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("4"), fMult);
403+
break;
404+
case 5:
405+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("5"), fMult);
406+
break;
407+
case 6:
408+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("6"), fMult);
409+
break;
410+
case 7:
411+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("7"), fMult);
412+
break;
413+
case 8:
414+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("8"), fMult);
415+
break;
416+
case 9:
417+
mHistogramQn->fill(HIST("Qn/mult_") + HIST("9"), fMult);
418+
break;
419+
default:
420+
return qnBin; // invalid qn bin
421+
}
422+
}
423+
return qnBin;
424+
}
425+
426+
/// \todo to be implemented!
427+
/// Fill cumulants histo for flow calculation
428+
/// Reset hists event-by-event
429+
/// \tparam T1 type of the collision
430+
/// \tparam T2 type of the tracks
431+
/// \param tracks All tracks
432+
template <typename T1, typename T2>
433+
void fillCumulants(T1 const& col, T2 const& tracks, float fHarmonic=2.f)
434+
{
435+
int numOfTracks = col.numContrib();
436+
if (numOfTracks < 3)
437+
return ;
438+
439+
mReQthisEvt->Reset();
440+
mImQthisEvt->Reset();
441+
mReQ2thisEvt->Reset();
442+
mImQ2thisEvt->Reset();
443+
mMQthisEvt->Reset();
444+
mMQWeightthisEvt->Reset();
445+
446+
for (auto const& track : tracks) {
447+
double weight=1; // Will implement NUA&NUE correction
448+
double phi = track.phi();
449+
double pt = track.pt();
450+
double eta = track.eta();
451+
double cosnphi = weight * TMath::Cos(fHarmonic*phi);
452+
double sinnphi = weight * TMath::Sin(fHarmonic*phi);
453+
double cos2nphi = weight * TMath::Cos(2*fHarmonic*phi);
454+
double sin2nphi = weight * TMath::Sin(2*fHarmonic*phi);
455+
mReQthisEvt->Fill(pt, eta, cosnphi);
456+
mImQthisEvt->Fill(pt, eta, sinnphi);
457+
mReQ2thisEvt->Fill(pt,eta,cos2nphi);
458+
mImQ2thisEvt->Fill(pt, eta, sin2nphi);
459+
mMQthisEvt ->Fill(pt, eta);
460+
mMQWeightthisEvt ->Fill(pt, eta, weight);
461+
}
462+
return;
463+
}
464+
465+
/// \todo to be implemented!
466+
/// Do cumulants for flow calculation
467+
/// \tparam T type of the collision
468+
/// \param doQnSeparation to fill flow in divied qn bins
469+
/// \param qnBin should be <int> in 0-9
470+
/// \param fEtaGap eta gap for flow cumulant
471+
template <typename T>
472+
void doCumulants(T const& col, bool doQnSeparation = false, int qnBin = -999, int mumQnBinNum = 10, float fEtaGap = 0.3f, int binPt = 100, int binEta = 32)
473+
{
474+
if (mMQthisEvt->Integral(1, binPt, 1, binEta) < 2)
475+
return;
476+
477+
double allReQ = mReQthisEvt ->Integral(1, binPt, 1, binEta);
478+
double allImQ = mImQthisEvt ->Integral(1, binPt, 1, binEta);
479+
TComplex Q(allReQ, allImQ);
480+
TComplex QStar = TComplex::Conjugate(Q);
481+
482+
double posEtaRe = mReQthisEvt->Integral(1, binPt, mReQthisEvt->GetYaxis()->FindBin(fEtaGap+1e-6), binEta);
483+
double posEtaIm = mImQthisEvt->Integral(1, binPt, mImQthisEvt->GetYaxis()->FindBin(fEtaGap+1e-6), binEta);
484+
if (mMQthisEvt->Integral(1, binPt, mMQthisEvt->GetYaxis()->FindBin(fEtaGap+1e-6), binEta) < 2)
485+
return;
486+
float posEtaMQ = mMQWeightthisEvt->Integral(1, binPt, mMQthisEvt->GetYaxis()->FindBin(fEtaGap+1e-6), binEta);
487+
TComplex posEtaQ = TComplex(posEtaRe, posEtaIm);
488+
TComplex posEtaQStar = TComplex::Conjugate(posEtaQ);
489+
490+
double negEtaRe = mReQthisEvt->Integral(1, binPt, 1, mReQthisEvt->GetYaxis()->FindBin(-1*fEtaGap-1e-6));
491+
double negEtaIm = mImQthisEvt->Integral(1, binPt, 1, mImQthisEvt->GetYaxis()->FindBin(-1*fEtaGap-1e-6));
492+
if (mMQthisEvt->Integral(1, binPt, 1, mMQthisEvt->GetYaxis()->FindBin(-1*fEtaGap-1e-6)) < 2)
493+
return;
494+
float negEtaMQ = mMQWeightthisEvt->Integral(1, binPt, 1, mMQthisEvt->GetYaxis()->FindBin(-1*fEtaGap-1e-6));
495+
TComplex negEtaQ = TComplex(negEtaRe, negEtaIm);
496+
TComplex negEtaQStar = TComplex::Conjugate(negEtaQ);
497+
498+
mHistogramQn->get<TProfile>(HIST("Event/profileC22"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
499+
if (doQnSeparation && qnBin >= 0 && qnBin < mumQnBinNum){
500+
switch (qnBin) {
501+
case 0:
502+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("0"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
503+
break;
504+
case 1:
505+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("1"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
506+
break;
507+
case 2:
508+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("2"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
509+
break;
510+
case 3:
511+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("3"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
512+
break;
513+
case 4:
514+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("4"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
515+
break;
516+
case 5:
517+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("5"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
518+
break;
519+
case 6:
520+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("6"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
521+
break;
522+
case 7:
523+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("7"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
524+
break;
525+
case 8:
526+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("8"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
527+
break;
528+
case 9:
529+
mHistogramQn->get<TProfile>(HIST("Qn/profileC22_") + HIST("9"))->Fill(col.centFT0C(), (negEtaQ*posEtaQStar).Re()/(negEtaMQ*posEtaMQ), (negEtaMQ*posEtaMQ));
530+
break;
531+
default:
532+
return; // invalid qn bin
533+
}
534+
}
535+
return;
536+
}
537+
538+
249539
private:
250-
HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output
540+
HistogramRegistry* mHistogramRegistry = nullptr; ///< For QA output
251541
bool mCutsSet = false; ///< Protection against running without cuts
252542
bool mCheckTrigger = false; ///< Check for trigger
253543
bool mCheckOffline = false; ///< Check for offline criteria (might change)
@@ -257,6 +547,13 @@ class FemtoDreamCollisionSelection
257547
float mZvtxMax = 999.f; ///< Maximal deviation from nominal z-vertex (cm)
258548
float mMinSphericity = 0.f;
259549
float mSphericityPtmin = 0.f;
550+
HistogramRegistry* mHistogramQn = nullptr; ///< For flow cumulant output
551+
TH2D* mReQthisEvt = nullptr; ///< For flow cumulant in an event
552+
TH2D* mImQthisEvt = nullptr; ///< For flow cumulant in an event
553+
TH2D* mReQ2thisEvt = nullptr; ///< For flow cumulant in an event
554+
TH2D* mImQ2thisEvt = nullptr; ///< For flow cumulant in an event
555+
TH2D* mMQthisEvt = nullptr; ///< For flow cumulant in an event
556+
TH2D* mMQWeightthisEvt = nullptr; ///< For flow cumulant in an event
260557
};
261558
} // namespace o2::analysis::femtoDream
262559

0 commit comments

Comments
 (0)