Skip to content

Commit 3dc60f2

Browse files
LuzhiyonggLuzhiyongg
authored
[PWGCF] add track density correction in MC (AliceO2Group#10578)
1 parent 00f3810 commit 3dc60f2

File tree

1 file changed

+69
-0
lines changed

1 file changed

+69
-0
lines changed

‎PWGCF/Flow/Tasks/flowMc.cxx‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -41,6 +41,7 @@
4141
#include <TProfile.h>
4242
#include <TRandom3.h>
4343
#include <TPDGCode.h>
44+
#include <TF1.h>
4445

4546
using namespace o2;
4647
using namespace o2::framework;
@@ -66,6 +67,10 @@ struct FlowMc {
6667
O2_DEFINE_CONFIGURABLE(cfgIsGlobalTrack, bool, false, "Use global tracks instead of hasTPC&&hasITS")
6768
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantEnabled, bool, false, "switch of calculating flow")
6869
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
70+
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
71+
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrSlopeFactor, float, 1.0f, "A factor to scale the track density efficiency slope")
72+
Configurable<std::vector<double>> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector<double>{0.6003720411, 0.6152630970, 0.6288860646, 0.6360694031, 0.6409494798, 0.6450540203, 0.6482117301, 0.6512592056, 0.6640008690, 0.6862631416, 0.7005738691, 0.7106567432, 0.7170728333}, "parameter 0 for track density efficiency correction"};
73+
Configurable<std::vector<double>> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector<double>{-1.007592e-05, -8.932635e-06, -9.114538e-06, -1.054818e-05, -1.220212e-05, -1.312304e-05, -1.376433e-05, -1.412813e-05, -1.289562e-05, -1.050065e-05, -8.635725e-06, -7.380821e-06, -6.201250e-06}, "parameter 1 for track density efficiency correction"};
6974

7075
ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
7176
ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
@@ -83,6 +88,12 @@ struct FlowMc {
8388
GFWWeights* mAcceptance = nullptr;
8489
bool correctionsLoaded = false;
8590

91+
std::vector<TF1*> funcEff;
92+
TH1D* hFindPtBin;
93+
TF1* funcV2;
94+
TF1* funcV3;
95+
TF1* funcV4;
96+
8697
// Connect to ccdb
8798
Service<ccdb::BasicCCDBManager> ccdb;
8899
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
@@ -134,6 +145,7 @@ struct FlowMc {
134145
histos.add<TH2>("hPtNchGlobal", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
135146
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
136147
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
148+
histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});
137149

138150
o2::framework::AxisSpec axis = axisPt;
139151
int nPtBins = axis.binEdges.size() - 1;
@@ -187,6 +199,25 @@ struct FlowMc {
187199
corrconfigsReco.push_back(fGFWReco->GetCorrelatorConfig("poiN10 refN10 | olN10 {2} refP10 {-2}", "Ch10Gap22", kTRUE));
188200
fGFWReco->CreateRegions();
189201
}
202+
203+
if (cfgTrackDensityCorrUse) {
204+
std::vector<double> pTEffBins = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.4, 1.8, 2.2, 2.6, 3.0};
205+
hFindPtBin = new TH1D("hFindPtBin", "hFindPtBin", pTEffBins.size() - 1, &pTEffBins[0]);
206+
funcEff.resize(pTEffBins.size() - 1);
207+
// LHC24g3 Eff
208+
std::vector<double> f1p0 = cfgTrackDensityP0;
209+
std::vector<double> f1p1 = cfgTrackDensityP1;
210+
for (uint ifunc = 0; ifunc < pTEffBins.size() - 1; ifunc++) {
211+
funcEff[ifunc] = new TF1(Form("funcEff%i", ifunc), "[0]+[1]*x", 0, 3000);
212+
funcEff[ifunc]->SetParameters(f1p0[ifunc], f1p1[ifunc] * cfgTrackDensityCorrSlopeFactor);
213+
}
214+
funcV2 = new TF1("funcV2", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
215+
funcV2->SetParameters(0.0186111, 0.00351907, -4.38264e-05, 1.35383e-07, -3.96266e-10);
216+
funcV3 = new TF1("funcV3", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
217+
funcV3->SetParameters(0.0174056, 0.000703329, -1.45044e-05, 1.91991e-07, -1.62137e-09);
218+
funcV4 = new TF1("funcV4", "[0]+[1]*x+[2]*x*x+[3]*x*x*x+[4]*x*x*x*x", 0, 100);
219+
funcV4->SetParameters(0.008845, 0.000259668, -3.24435e-06, 4.54837e-08, -6.01825e-10);
220+
}
190221
}
191222

192223
void loadCorrections(uint64_t timestamp)
@@ -305,6 +336,12 @@ struct FlowMc {
305336
fGFWReco->Clear();
306337
}
307338

339+
double psi2Est = 0, psi3Est = 0, psi4Est = 0;
340+
float wEPeff = 1;
341+
double v2 = 0, v3 = 0, v4 = 0;
342+
double q2x = 0, q2y = 0;
343+
double q3x = 0, q3y = 0;
344+
double q4x = 0, q4y = 0;
308345
for (auto const& mcParticle : mcParticles) {
309346
int pdgCode = std::abs(mcParticle.pdgCode());
310347
if (pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton)
@@ -318,9 +355,28 @@ struct FlowMc {
318355
for (auto const& track : tracks) {
319356
if (track.hasTPC() && track.hasITS())
320357
nChGlobal++;
358+
if (cfgTrackDensityCorrUse && cfgFlowCumulantEnabled) {
359+
bool withinPtRef = (cfgCutPtRefMin < track.pt()) && (track.pt() < cfgCutPtRefMax); // within RF pT rang
360+
if (withinPtRef) {
361+
q2x += std::cos(2 * track.phi());
362+
q2y += std::sin(2 * track.phi());
363+
q3x += std::cos(3 * track.phi());
364+
q3y += std::sin(3 * track.phi());
365+
q4x += std::cos(4 * track.phi());
366+
q4y += std::sin(4 * track.phi());
367+
}
368+
}
321369
}
322370
}
323371
}
372+
if (cfgTrackDensityCorrUse && cfgFlowCumulantEnabled) {
373+
psi2Est = std::atan2(q2y, q2x) / 2.;
374+
psi3Est = std::atan2(q3y, q3x) / 3.;
375+
psi4Est = std::atan2(q4y, q4x) / 4.;
376+
v2 = funcV2->Eval(centrality);
377+
v3 = funcV3->Eval(centrality);
378+
v4 = funcV4->Eval(centrality);
379+
}
324380

325381
for (auto const& mcParticle : mcParticles) {
326382
// focus on bulk: e, mu, pi, k, p
@@ -382,6 +438,19 @@ struct FlowMc {
382438
fWeights->fill(mcParticle.phi(), mcParticle.eta(), vtxz, mcParticle.pt(), 0, 0);
383439
if (!setCurrentParticleWeights(weff, wacc, mcParticle.phi(), mcParticle.eta(), mcParticle.pt(), vtxz))
384440
continue;
441+
if (cfgTrackDensityCorrUse && cfgFlowCumulantEnabled && withinPtRef) {
442+
double fphi = v2 * std::cos(2 * (mcParticle.phi() - psi2Est)) + v3 * std::cos(3 * (mcParticle.phi() - psi3Est)) + v4 * std::cos(4 * (mcParticle.phi() - psi4Est));
443+
fphi = (1 + 2 * fphi);
444+
int pTBinForEff = hFindPtBin->FindBin(mcParticle.pt());
445+
if (pTBinForEff >= 1 && pTBinForEff <= hFindPtBin->GetNbinsX()) {
446+
wEPeff = funcEff[pTBinForEff - 1]->Eval(fphi * nChGlobal);
447+
if (wEPeff > 0.) {
448+
wEPeff = 1. / wEPeff;
449+
weff *= wEPeff;
450+
histos.fill(HIST("hPhiWeightedTrDen"), mcParticle.phi(), wacc * wEPeff);
451+
}
452+
}
453+
}
385454

386455
if (cfgFlowCumulantEnabled) {
387456
if (withinPtRef)

0 commit comments

Comments
 (0)