Skip to content

Commit c3f1685

Browse files
Preet-BhanjanPreet Pati
Preet-Bhanjan
and
Preet Pati
authored
[PWGCF] Added weights for the particles and flow for the resonances (#10250)
Co-authored-by: Preet Pati <preet@preet-2.local>
1 parent ecf4f5a commit c3f1685

File tree

2 files changed

+243
-56
lines changed

2 files changed

+243
-56
lines changed

‎PWGCF/Flow/Tasks/flowPbpbPikp.cxx‎

Lines changed: 143 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -19,6 +19,7 @@
1919
#include <utility>
2020
#include <array>
2121
#include <string>
22+
#include <map>
2223

2324
#include "Math/Vector4D.h"
2425

@@ -75,7 +76,8 @@ struct FlowPbpbPikp {
7576
O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Use Nch for flow observables")
7677
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
7778
O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, true, "Fill NUA weights")
78-
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
79+
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, true, "Fill and output NUA weights")
80+
O2_DEFINE_CONFIGURABLE(cfgOutputRunByRun, bool, true, "Fill and output NUA weights run by run")
7981
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
8082
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
8183
O2_DEFINE_CONFIGURABLE(cfgTpcNsigmaCut, float, 2.0f, "TPC N-sigma cut for pions, kaons, protons")
@@ -107,6 +109,20 @@ struct FlowPbpbPikp {
107109
TAxis* fPtAxis;
108110
TRandom3* fRndm = new TRandom3(0);
109111

112+
std::map<int, std::vector<std::shared_ptr<TH3>>> th3sList;
113+
enum OutputSpecies {
114+
hRef = 0,
115+
hCharge,
116+
hPion,
117+
hKaon,
118+
hProton,
119+
kCount_OutputSpecies
120+
};
121+
int lastRunNumer = -1;
122+
std::vector<int> runNumbers;
123+
std::vector<GFWWeights*> mAcceptance;
124+
bool correctionsLoaded = false;
125+
110126
void init(InitContext const&)
111127
{
112128
ccdb->setURL(ccdbUrl.value);
@@ -130,6 +146,13 @@ struct FlowPbpbPikp {
130146
histos.add("c24_gap08_pr", "", {HistType::kTProfile, {axisMultiplicity}});
131147
histos.add("TofTpcNsigma", "", {HistType::kTHnSparseD, {{axisParticles, axisNsigmaTPC, axisNsigmaTOF, axisPt}}});
132148
histos.add("partCount", "", {HistType::kTHnSparseD, {{axisParticles, axisMultiplicity, axisPt}}});
149+
if (cfgOutputNUAWeights && !cfgOutputRunByRun) {
150+
histos.add<TH3>("NUA/hPhiEtaVtxz_ref", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
151+
histos.add<TH3>("NUA/hPhiEtaVtxz_ch", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
152+
histos.add<TH3>("NUA/hPhiEtaVtxz_pi", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
153+
histos.add<TH3>("NUA/hPhiEtaVtxz_ka", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
154+
histos.add<TH3>("NUA/hPhiEtaVtxz_pr", ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
155+
}
133156

134157
o2::framework::AxisSpec axis = axisPt;
135158
int nPtBins = axis.binEdges.size() - 1;
@@ -296,10 +319,10 @@ struct FlowPbpbPikp {
296319
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
297320
{
298321
double dnx, val;
299-
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
300-
if (dnx == 0)
301-
return;
302322
if (!corrconf.pTDif) {
323+
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
324+
if (dnx == 0)
325+
return;
303326
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
304327
if (std::fabs(val) < 1)
305328
histos.fill(tarName, cent, val, dnx);
@@ -319,11 +342,11 @@ struct FlowPbpbPikp {
319342
void fillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
320343
{
321344
double dnx, val;
322-
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
323-
if (dnx == 0) {
324-
return;
325-
}
326345
if (!corrconf.pTDif) {
346+
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
347+
if (dnx == 0) {
348+
return;
349+
}
327350
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
328351
if (std::fabs(val) < 1) {
329352
fFC->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
@@ -341,6 +364,80 @@ struct FlowPbpbPikp {
341364
return;
342365
}
343366

367+
void createRunByRunHistos(int runNumber)
368+
{
369+
if (cfgOutputNUAWeights) {
370+
std::vector<std::shared_ptr<TH3>> tH3s(kCount_OutputSpecies);
371+
tH3s[hRef] = histos.add<TH3>(Form("NUA/%d/hPhiEtaVtxz_ref", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
372+
tH3s[hCharge] = histos.add<TH3>(Form("NUA/%d/hPhiEtaVtxz_ch", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
373+
tH3s[hPion] = histos.add<TH3>(Form("NUA/%d/hPhiEtaVtxz_pi", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
374+
tH3s[hKaon] = histos.add<TH3>(Form("NUA/%d/hPhiEtaVtxz_ka", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
375+
tH3s[hProton] = histos.add<TH3>(Form("NUA/%d/hPhiEtaVtxz_pr", runNumber), ";#varphi;#eta;v_{z}", {HistType::kTH3D, {axisPhi, {64, -1.6, 1.6}, {40, -10, 10}}});
376+
th3sList.insert(std::make_pair(runNumber, tH3s));
377+
}
378+
}
379+
380+
void loadCorrections(aod::BCsWithTimestamps::iterator const& bc)
381+
{
382+
if (correctionsLoaded)
383+
return;
384+
if (!cfgAcceptance.value.empty()) {
385+
uint64_t timestamp = bc.timestamp();
386+
mAcceptance.clear();
387+
mAcceptance.resize(kCount_OutputSpecies);
388+
mAcceptance.push_back(ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value + "_ref", timestamp));
389+
mAcceptance.push_back(ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value + "_ch", timestamp));
390+
mAcceptance.push_back(ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value + "_pi", timestamp));
391+
mAcceptance.push_back(ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value + "_ka", timestamp));
392+
mAcceptance.push_back(ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value + "_pr", timestamp));
393+
}
394+
395+
correctionsLoaded = true;
396+
}
397+
398+
template <typename TTrack>
399+
double getAcceptance(TTrack track, const double& vtxz, int index)
400+
{ // 0 ref, 1 ch, 2 pi, 3 ka, 4 pr
401+
double wacc = 1;
402+
if (!mAcceptance.empty())
403+
wacc = mAcceptance[index]->getNUA(track.phi(), track.eta(), vtxz);
404+
return wacc;
405+
}
406+
407+
template <typename TTrack>
408+
void fillWeights(const TTrack track, const double vtxz, const int& pid_index, const int& run)
409+
{
410+
double pt = track.pt();
411+
bool withinPtPOI = (cfgCutPtPOIMin < pt) && (pt < cfgCutPtPOIMax); // within POI pT range
412+
bool withinPtRef = (cfgCutPtMin < pt) && (pt < cfgCutPtMax); // within RF pT range
413+
414+
if (cfgOutputRunByRun) {
415+
if (withinPtRef && !pid_index)
416+
th3sList[run][hRef]->Fill(track.phi(), track.eta(), vtxz); // pt-subset of charged particles for ref flow
417+
if (withinPtPOI)
418+
th3sList[run][hCharge + pid_index]->Fill(track.phi(), track.eta(), vtxz); // charged and id'ed particle weights
419+
} else {
420+
if (withinPtRef && !pid_index)
421+
histos.fill(HIST("NUA/hPhiEtaVtxz_ref"), track.phi(), track.eta(), vtxz); // pt-subset of charged particles for ref flow
422+
if (withinPtPOI) {
423+
switch (pid_index) {
424+
case 0:
425+
histos.fill(HIST("NUA/hPhiEtaVtxz_ch"), track.phi(), track.eta(), vtxz); // charged particle weights
426+
break;
427+
case 1:
428+
histos.fill(HIST("NUA/hPhiEtaVtxz_pi"), track.phi(), track.eta(), vtxz); // pion weights
429+
break;
430+
case 2:
431+
histos.fill(HIST("NUA/hPhiEtaVtxz_ka"), track.phi(), track.eta(), vtxz); // kaon weights
432+
break;
433+
case 3:
434+
histos.fill(HIST("NUA/hPhiEtaVtxz_pr"), track.phi(), track.eta(), vtxz); // proton weights
435+
break;
436+
}
437+
}
438+
}
439+
}
440+
344441
void process(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracksWithoutBayes const& tracks)
345442
{
346443
int nTot = tracks.size();
@@ -352,13 +449,23 @@ struct FlowPbpbPikp {
352449
float lRandom = fRndm->Rndm();
353450
float vtxz = collision.posZ();
354451
const auto cent = collision.centFT0C();
452+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
453+
int runNumber = bc.runNumber();
454+
if (cfgOutputRunByRun && runNumber != lastRunNumer) {
455+
lastRunNumer = runNumber;
456+
if (std::find(runNumbers.begin(), runNumbers.end(), runNumber) == runNumbers.end()) {
457+
// if run number is not in the preconfigured list, create new output histograms for this run
458+
createRunByRunHistos(runNumber);
459+
runNumbers.push_back(runNumber);
460+
}
461+
}
355462

356463
histos.fill(HIST("hVtxZ"), vtxz);
357464
histos.fill(HIST("hMult"), nTot);
358465
histos.fill(HIST("hCent"), collision.centFT0C());
359466
fGFW->Clear();
360467

361-
float weff = 1, wacc = 1;
468+
float weff = 1;
362469
int pidIndex;
363470

364471
for (auto const& track : tracks) {
@@ -376,38 +483,48 @@ struct FlowPbpbPikp {
376483

377484
// pidIndex = getBayesPIDIndex(track);
378485
pidIndex = getNsigmaPID(track);
486+
if (cfgOutputNUAWeights)
487+
fillWeights(track, vtxz, pidIndex, runNumber);
488+
489+
if (!withinPtPOI && !withinPtRef)
490+
return;
491+
double waccRef = getAcceptance(track, vtxz, 0);
492+
double waccPOI = withinPtPOI ? getAcceptance(track, vtxz, pidIndex + 1) : getAcceptance(track, vtxz, 0);
493+
if (withinPtRef && withinPtPOI && pidIndex)
494+
waccRef = waccPOI; // if particle is both (then it's overlap), override ref with POI
495+
379496
if (withinPtRef) {
380-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1);
381-
fGFW->Fill(track.eta(), 1, track.phi(), wacc * weff, 512);
497+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccRef * weff, 1);
498+
fGFW->Fill(track.eta(), 1, track.phi(), waccRef * weff, 512);
382499
}
383500
if (withinPtPOI) {
384-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 128);
385-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1024);
501+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 128);
502+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 1024);
386503
}
387504
if (withinPtPOI && withinPtRef) {
388-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 256);
389-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 2048);
505+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 256);
506+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 2048);
390507
}
391508

392509
if (pidIndex) {
393510
histos.fill(HIST("partCount"), pidIndex - 1, cent, pt);
394511
if (withinPtPOI)
395-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex));
512+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 1 << (pidIndex));
396513
if (withinPtPOI && withinPtRef)
397-
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), wacc * weff, 1 << (pidIndex + 3));
514+
fGFW->Fill(track.eta(), fPtAxis->FindBin(pt) - 1, track.phi(), waccPOI * weff, 1 << (pidIndex + 3));
398515
}
399516
} // track loop ends
400517

401518
// Filling cumulants with ROOT TProfile
402-
fillProfile(corrconfigs.at(0), HIST("c22_gap08"), cent);
403-
fillProfile(corrconfigs.at(1), HIST("c22_gap08_pi"), cent);
404-
fillProfile(corrconfigs.at(2), HIST("c22_gap08_ka"), cent);
405-
fillProfile(corrconfigs.at(3), HIST("c22_gap08_pr"), cent);
406-
fillProfile(corrconfigs.at(4), HIST("c24_full"), cent);
407-
fillProfile(corrconfigs.at(5), HIST("c24_gap08"), cent);
408-
fillProfile(corrconfigs.at(6), HIST("c24_gap08_pi"), cent);
409-
fillProfile(corrconfigs.at(7), HIST("c24_gap08_ka"), cent);
410-
fillProfile(corrconfigs.at(8), HIST("c24_gap08_pr"), cent);
519+
fillProfile(corrconfigs.at(1), HIST("c22_gap08"), cent);
520+
fillProfile(corrconfigs.at(2), HIST("c22_gap08_pi"), cent);
521+
fillProfile(corrconfigs.at(3), HIST("c22_gap08_ka"), cent);
522+
fillProfile(corrconfigs.at(4), HIST("c22_gap08_pr"), cent);
523+
fillProfile(corrconfigs.at(5), HIST("c24_full"), cent);
524+
fillProfile(corrconfigs.at(6), HIST("c24_gap08"), cent);
525+
fillProfile(corrconfigs.at(7), HIST("c24_gap08_pi"), cent);
526+
fillProfile(corrconfigs.at(8), HIST("c24_gap08_ka"), cent);
527+
fillProfile(corrconfigs.at(9), HIST("c24_gap08_pr"), cent);
411528

412529
for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {
413530
fillFC(corrconfigs.at(l_ind), cent, lRandom);

0 commit comments

Comments
 (0)