Skip to content

Commit 46ba359

Browse files
wenyaCernwenyaCern
committed
Add functions for femto in qn bins
1 parent 5903b01 commit 46ba359

File tree

5 files changed

+689
-6
lines changed

5 files changed

+689
-6
lines changed

‎PWGCF/DataModel/FemtoDerived.h‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -50,6 +50,9 @@ 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); //! Bin of qn-vector of the event
55+
5356
} // namespace femtodreamcollision
5457

5558
DECLARE_SOA_TABLE_STAGED(FDCollisions, "FDCOLLISION",
@@ -69,6 +72,10 @@ DECLARE_SOA_TABLE(FDColMasks, "AOD", "FDCOLMASK",
6972
DECLARE_SOA_TABLE(FDDownSample, "AOD", "FDDOWNSAMPLE",
7073
femtodreamcollision::Downsample);
7174

75+
DECLARE_SOA_TABLE(FDExtQnCollisions, "AOD", "FDEXTQNCOLLISION",
76+
femtodreamcollision::QnBin);
77+
using FDExtCollision = FDExtQnCollisions::iterator;
78+
7279
namespace femtodreamMCcollision
7380
{
7481
DECLARE_SOA_COLUMN(MultMCgenPartEta08, multMCgenPartEta08, int); //! Multiplicity of the event as given by the generator in |eta|<0.8

‎PWGCF/FemtoDream/Core/femtoDreamCollisionSelection.h‎

Lines changed: 311 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,41 @@ 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+
{
170+
int tpcOccupancyMin = 0;
171+
int tpcOccupancyMax = 1000;
172+
const auto occupancy = col.trackOccupancyInTimeRange();
173+
if ((occupancy < tpcOccupancyMin || occupancy > tpcOccupancyMax)) {
174+
return false;
175+
}
176+
177+
bool noSameBunchPileup = true;
178+
bool isGoodZvtxFT0vsPV = true;
179+
bool isGoodITSLayersAll = true;
180+
bool noCollInRofStandard = true;
181+
bool noHighMultCollInPrevRof = true;
182+
bool noCollInTimeRangeStandard = true;
183+
// bool isVertexITSTPC = true;
184+
if ((noSameBunchPileup && !col.selection_bit(aod::evsel::kNoSameBunchPileup))
185+
|| (isGoodZvtxFT0vsPV && !col.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV))
186+
|| (isGoodITSLayersAll && !col.selection_bit(aod::evsel::kIsGoodITSLayersAll))
187+
|| (noCollInRofStandard && !col.selection_bit(aod::evsel::kNoCollInRofStandard))
188+
|| (noHighMultCollInPrevRof && !col.selection_bit(aod::evsel::kNoHighMultCollInPrevRof))
189+
|| (noCollInTimeRangeStandard && !col.selection_bit(aod::evsel::kNoCollInTimeRangeStandard))
190+
// || (isVertexITSTPC && !col.selection_bit(aod::evsel::kIsVertexITSTPC))
191+
) {
192+
return false;
193+
}
194+
195+
return true;
196+
}
197+
159198
/// Some basic QA of the event
160199
/// \tparam T type of the collision
161200
/// \param col Collision
@@ -175,6 +214,37 @@ class FemtoDreamCollisionSelection
175214
}
176215
}
177216

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

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

0 commit comments

Comments
 (0)