Skip to content

Commit 3018aba

Browse files
cnkostercnkoster
committed
Add extra NUA option and extra QA plots
1 parent e2e14a5 commit 3018aba

File tree

1 file changed

+93
-25
lines changed

1 file changed

+93
-25
lines changed

‎PWGCF/Flow/Tasks/flowSP.cxx‎

Lines changed: 93 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -85,6 +85,8 @@ struct FlowSP {
8585
O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, true, "Fill NUA weights");
8686
O2_DEFINE_CONFIGURABLE(cfgFillWeightsPOS, bool, true, "Fill NUA weights only for positive charges");
8787
O2_DEFINE_CONFIGURABLE(cfgFillWeightsNEG, bool, true, "Fill NUA weights only for negative charges");
88+
O2_DEFINE_CONFIGURABLE(cfguseNUA1D, bool, false, "Use 1D NUA weights (only phi)");
89+
O2_DEFINE_CONFIGURABLE(cfguseNUA2D, bool, true, "Use 2D NUA weights (phi and eta)");
8890
// Additional track Selections
8991
O2_DEFINE_CONFIGURABLE(cfgTrackSelsUseAdditionalTrackCut, bool, false, "Bool to enable Additional Track Cut");
9092
O2_DEFINE_CONFIGURABLE(cfgTrackSelsDoDCApt, bool, false, "Apply Pt dependent DCAz cut");
@@ -146,6 +148,7 @@ struct FlowSP {
146148
struct Config {
147149
std::vector<TH1D*> mEfficiency = {};
148150
std::vector<GFWWeights*> mAcceptance = {};
151+
std::vector<TH3D*> mAcceptance2D = {};
149152
bool correctionsLoaded = false;
150153
int lastRunNumber = 0;
151154

@@ -161,10 +164,10 @@ struct FlowSP {
161164

162165
} cfg;
163166

164-
// define output objects
165167
OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
166168
OutputObj<GFWWeights> fWeightsPOS{GFWWeights("weights_positive")};
167169
OutputObj<GFWWeights> fWeightsNEG{GFWWeights("weights_negative")};
170+
168171
HistogramRegistry registry{"registry"};
169172

170173
// Event selection cuts - Alex
@@ -247,7 +250,7 @@ struct FlowSP {
247250
AxisSpec axisDCAxy = {100, -.5, .5, "DCA_{xy} (cm)"};
248251
AxisSpec axisPhiMod = {100, 0, constants::math::PI / 9, "fmod(#varphi,#pi/9)"};
249252
AxisSpec axisPhi = {60, 0, constants::math::TwoPI, "#varphi"};
250-
AxisSpec axisEta = {64, -1.8, 1.8, "#eta"};
253+
AxisSpec axisEta = {64, -1.6,1.6, "#eta"};
251254
AxisSpec axisEtaVn = {8, -.8, .8, "#eta"};
252255
AxisSpec axisVx = {40, -0.01, 0.01, "v_{x}"};
253256
AxisSpec axisVy = {40, -0.01, 0.01, "v_{y}"};
@@ -295,14 +298,22 @@ struct FlowSP {
295298
registry.get<TH1>(HIST("hTrackCount"))->GetXaxis()->SetBinLabel(trackSel_ParticleWeights + 1, "Apply weights");
296299

297300
if (cfgFillWeights) {
298-
fWeights->setPtBins(ptbins, &ptbinning[0]);
299-
fWeights->init(true, false);
301+
if (cfguseNUA2D) {
302+
registry.add<TH3>("weights/hPhi_Eta_vz", "", kTH3D, {axisPhi, axisEta, axisVz});
303+
registry.add<TH3>("weights/hPhi_Eta_vz_positive", "", kTH3D, {axisPhi, axisEta, axisVz});
304+
registry.add<TH3>("weights/hPhi_Eta_vz_negative", "", kTH3D, {axisPhi, axisEta, axisVz});
305+
}
306+
else {
307+
// define output objects
308+
fWeights->setPtBins(ptbins, &ptbinning[0]);
309+
fWeights->init(true, false);
300310

301-
fWeightsPOS->setPtBins(ptbins, &ptbinning[0]);
302-
fWeightsPOS->init(true, false);
311+
fWeightsPOS->setPtBins(ptbins, &ptbinning[0]);
312+
fWeightsPOS->init(true, false);
303313

304-
fWeightsNEG->setPtBins(ptbins, &ptbinning[0]);
305-
fWeightsNEG->init(true, false);
314+
fWeightsNEG->setPtBins(ptbins, &ptbinning[0]);
315+
fWeightsNEG->init(true, false);
316+
}
306317
}
307318

308319
if (cfgFillEventQA) {
@@ -357,9 +368,14 @@ struct FlowSP {
357368
if (cfgFillTrackQA) {
358369
registry.add("QA/after/pt_phi", "", {HistType::kTH2D, {axisPt, axisPhiMod}});
359370
registry.add<TH1>("incl/QA/after/hPt", "", kTH1D, {axisPt});
371+
registry.add<TH1>("incl/QA/after/hPt_forward", "", kTH1D, {axisPt});
372+
registry.add<TH1>("incl/QA/after/hPt_forward_uncorrected", "", kTH1D, {axisPt});
373+
registry.add<TH1>("incl/QA/after/hPt_backward", "", kTH1D, {axisPt});
374+
registry.add<TH1>("incl/QA/after/hPt_backward_uncorrected", "", kTH1D, {axisPt});
360375
registry.add<TH1>("incl/QA/after/hPhi", "", kTH1D, {axisPhi});
361376
registry.add<TH1>("incl/QA/after/hPhi_uncorrected", "", kTH1D, {axisPhi});
362377
registry.add<TH1>("incl/QA/after/hEta", "", kTH1D, {axisEta});
378+
registry.add<TH1>("incl/QA/after/hEta_uncorrected", "", kTH1D, {axisEta});
363379
registry.add<TH3>("incl/QA/after/hPhi_Eta_vz", "", kTH3D, {axisPhi, axisEta, axisVz});
364380
registry.add<TH3>("incl/QA/after/hPhi_Eta_vz_corrected", "", kTH3D, {axisPhi, axisEta, axisVz});
365381
registry.add<TH2>("incl/QA/after/hDCAxy_pt", "", kTH2D, {axisPt, axisDCAxy});
@@ -568,6 +584,17 @@ struct FlowSP {
568584
}
569585
} // end of init
570586

587+
float getNUA2D(TH3D* hNUA, float eta, float phi, float vtxz)
588+
{
589+
int xind = hNUA->GetXaxis()->FindBin(phi);
590+
int etaind = hNUA->GetYaxis()->FindBin(eta);
591+
int vzind = hNUA->GetZaxis()->FindBin(vtxz);
592+
float weight = hNUA->GetBinContent(xind, etaind, vzind);
593+
if (weight != 0)
594+
return 1. / weight;
595+
return 1;
596+
}
597+
571598
template <typename TrackObject>
572599
int getTrackPID(TrackObject track)
573600
{
@@ -639,19 +666,34 @@ struct FlowSP {
639666

640667
int nWeights = 3;
641668

642-
if (cfgCCDB_NUA.value.empty() == false) {
643-
TList* listCorrections = ccdb->getForTimeStamp<TList>(cfgCCDB_NUA, timestamp);
644-
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights")));
645-
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_positive")));
646-
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_negative")));
647-
int sizeAcc = cfg.mAcceptance.size();
648-
if (sizeAcc < nWeights)
649-
LOGF(fatal, "Could not load acceptance weights from %s", cfgCCDB_NUA.value.c_str());
650-
else
651-
LOGF(info, "Loaded acceptance weights from %s", cfgCCDB_NUA.value.c_str());
652-
} else {
653-
LOGF(info, "cfgCCDB_NUA empty! No corrections loaded");
669+
if(cfguseNUA1D) {
670+
if (cfgCCDB_NUA.value.empty() == false) {
671+
TList* listCorrections = ccdb->getForTimeStamp<TList>(cfgCCDB_NUA, timestamp);
672+
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights")));
673+
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_positive")));
674+
cfg.mAcceptance.push_back(reinterpret_cast<GFWWeights*>(listCorrections->FindObject("weights_negative")));
675+
int sizeAcc = cfg.mAcceptance.size();
676+
if (sizeAcc < nWeights)
677+
LOGF(fatal, "Could not load acceptance weights from %s", cfgCCDB_NUA.value.c_str());
678+
else
679+
LOGF(info, "Loaded acceptance weights from %s", cfgCCDB_NUA.value.c_str());
680+
} else {
681+
LOGF(info, "cfgCCDB_NUA empty! No corrections loaded");
682+
}
683+
} else if(cfguseNUA2D) {
684+
if (cfgCCDB_NUA.value.empty() == false) {
685+
TH3D* NUA2D = ccdb->getForTimeStamp<TH3D>(cfgCCDB_NUA, timestamp);
686+
if (!NUA2D){
687+
LOGF(fatal, "Could not load acceptance weights from %s", cfgCCDB_NUA.value.c_str());
688+
} else {
689+
LOGF(info, "Loaded acceptance weights from %s", cfgCCDB_NUA.value.c_str());
690+
cfg.mAcceptance2D.push_back(NUA2D);
691+
}
692+
} else {
693+
LOGF(info, "cfgCCDB_NUA empty! No corrections loaded");
694+
}
654695
}
696+
// Get Efficiency correction
655697
if (cfgCCDB_NUE.value.empty() == false) {
656698
TList* listCorrections = ccdb->getForTimeStamp<TList>(cfgCCDB_NUE, timestamp);
657699
cfg.mEfficiency.push_back(reinterpret_cast<TH1D*>(listCorrections->FindObject("Efficiency")));
@@ -680,11 +722,20 @@ struct FlowSP {
680722
if (eff == 0)
681723
return false;
682724
weight_nue = 1. / eff;
683-
int sizeAcc = cfg.mAcceptance.size();
684-
if (sizeAcc > pID) {
685-
weight_nua = cfg.mAcceptance[pID]->getNUA(phi, eta, vtxz);
686-
} else {
687-
weight_nua = 1;
725+
726+
if(cfguseNUA1D) {
727+
int sizeAcc = cfg.mAcceptance.size();
728+
if (sizeAcc > pID) {
729+
weight_nua = cfg.mAcceptance[pID]->getNUA(phi, eta, vtxz);
730+
} else {
731+
weight_nua = 1;
732+
}
733+
} else if(cfguseNUA2D) {
734+
if(cfg.mAcceptance2D.size() > 0) {
735+
weight_nua = getNUA2D(cfg.mAcceptance2D[0], eta, phi, vtxz);
736+
} else {
737+
weight_nua = 1;
738+
}
688739
}
689740
return true;
690741
}
@@ -1004,9 +1055,17 @@ struct FlowSP {
10041055
static constexpr std::string_view Time[] = {"before/", "after/"};
10051056
// NOTE: species[kUnidentified] = "" (when no PID)
10061057
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPt"), track.pt(), wacc * weff);
1058+
if (track.eta() > 0) {
1059+
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPt_forward"), track.pt(), wacc * weff);
1060+
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPt_forward_uncorrected"), track.pt());
1061+
} else {
1062+
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPt_backward"), track.pt(), wacc * weff);
1063+
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPt_backward_uncorrected"), track.pt());
1064+
}
10071065
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPhi"), track.phi(), wacc);
10081066
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPhi_uncorrected"), track.phi());
10091067
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hEta"), track.eta(), wacc);
1068+
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hEta_uncorrected"), track.eta());
10101069
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPhi_Eta_vz"), track.phi(), track.eta(), vz);
10111070
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hPhi_Eta_vz_corrected"), track.phi(), track.eta(), vz, wacc);
10121071
registry.fill(HIST(Charge[ct]) + HIST(Species[pt]) + HIST("QA/") + HIST(Time[ft]) + HIST("hDCAxy_pt"), track.pt(), track.dcaXY(), wacc * weff);
@@ -1253,6 +1312,15 @@ struct FlowSP {
12531312
// constrain angle to 0 -> [0,0+2pi]
12541313
auto phi = RecoDecay::constrainAngle(track.phi(), 0);
12551314

1315+
if (cfguseNUA2D && cfgFillWeights) {
1316+
registry.fill(HIST("weights/hPhi_Eta_vz"), phi, track.eta(), vtxz, 1);
1317+
if(pos) {
1318+
registry.fill(HIST("weights/hPhi_Eta_vz_positive"), phi, track.eta(), vtxz, 1);
1319+
} else {
1320+
registry.fill(HIST("weights/hPhi_Eta_vz_negative"), phi, track.eta(), vtxz, 1);
1321+
}
1322+
}
1323+
12561324
// Fill NUA weights (last 0 is for Data see GFWWeights class (not a weight))
12571325
if (cfgFillWeights) {
12581326
fWeights->fill(phi, track.eta(), vtxz, track.pt(), centrality, 0);

0 commit comments

Comments
 (0)