Skip to content

Commit 944369f

Browse files
yashpatleyyashpatley
authored
[PWGCF] : Update flowEventPlane.cxx (#13578)
1 parent 85ac022 commit 944369f

File tree

1 file changed

+112
-84
lines changed

1 file changed

+112
-84
lines changed

‎PWGCF/Flow/Tasks/flowEventPlane.cxx‎

Lines changed: 112 additions & 84 deletions
Original file line numberDiff line numberDiff line change
@@ -84,8 +84,12 @@ struct FlowEventPlane {
8484
Configurable<float> cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"};
8585
Configurable<float> cTrackDcaZCut{"cTrackDcaZCut", 1., "DcaXY Cut"};
8686

87+
// Gain calibration
88+
Configurable<bool> cDoGainCalib{"cDoGainCalib", false, "Gain Calib Flag"};
89+
Configurable<bool> cUseAlphaZDC{"cUseAlphaZDC", true, "Use Alpha ZDC"};
90+
8791
// Coarse binning factor
88-
Configurable<int> cAxisCBF{"cAxisCBF", 1, "Coarse Bin Factor"};
92+
Configurable<int> cAxisCBF{"cAxisCBF", 5, "Coarse Bin Factor"};
8993

9094
// Cent Vx Vy Vz Bins
9195
Configurable<int> cAxisCentBins{"cAxisCentBins", 20, "NBins Centrality"};
@@ -97,6 +101,7 @@ struct FlowEventPlane {
97101
Configurable<float> cAxisVyMax{"cAxisVyMax", 0.006, "Vy Max"};
98102

99103
// Corrections
104+
Configurable<bool> cApplyRecentCorr{"cApplyRecentCorr", false, "Apply recentering"};
100105
Configurable<std::vector<int>> cCorrFlagVector{"cCorrFlagVector", {0, 0, 0, 0, 0, 0}, "Correction Flag"};
101106

102107
// CCDB
@@ -110,14 +115,9 @@ struct FlowEventPlane {
110115
HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject};
111116

112117
// Global objects
118+
const float zdcDenThrs = 1e-4;
113119
float cent = 0., mult = 0.;
114120
float posX = 0., posY = 0., posZ = 0.;
115-
116-
std::array<float, 4> znaXWeigthEnergy = {1., 1., 1., 1.};
117-
std::array<float, 4> znaYWeigthEnergy = {1., 1., 1., 1.};
118-
std::array<float, 4> zncXWeigthEnergy = {1., 1., 1., 1.};
119-
std::array<float, 4> zncYWeigthEnergy = {1., 1., 1., 1.};
120-
121121
std::vector<std::vector<std::string>> vCoarseCorrHistNames = {
122122
{"hXZNAVsCentVxVyVz"},
123123
{"hYZNAVsCentVxVyVz"},
@@ -128,18 +128,8 @@ struct FlowEventPlane {
128128
{"hYZNAVsCent", "hYZNAVsVx", "hYZNAVsVy", "hYZNAVsVz"},
129129
{"hXZNCVsCent", "hXZNCVsVx", "hXZNCVsVy", "hXZNCVsVz"},
130130
{"hYZNCVsCent", "hYZNCVsVx", "hYZNCVsVy", "hYZNCVsVz"}};
131-
132-
// Map for Correction Type and Histogram Names
133131
std::map<CorrectionType, std::vector<std::vector<std::string>>> corrTypeHistNameMap = {{kFineCorr, vFineCorrHistNames}, {kCoarseCorr, vCoarseCorrHistNames}};
134132

135-
// Structure to hold CCDB objects
136-
struct CcdbObjects {
137-
TList* ccdbList;
138-
TObject* obj;
139-
TProfile* hp;
140-
THnSparseF* hn;
141-
} ccdbObjects;
142-
143133
void init(InitContext const&)
144134
{
145135
// Set CCDB url
@@ -185,14 +175,12 @@ struct FlowEventPlane {
185175
histos.add("Event/hVx", "V_{x}", kTH1F, {axisVx});
186176
histos.add("Event/hVy", "V_{y}", kTH1F, {axisVy});
187177
histos.add("Event/hVz", "V_{z}", kTH1F, {axisVz});
188-
histos.add("QA/hZNASignalSector1", "ZNA Signal Sector 1", kTH1F, {axisZDCEnergy});
189-
histos.add("QA/hZNASignalSector2", "ZNA Signal Sector 2", kTH1F, {axisZDCEnergy});
190-
histos.add("QA/hZNASignalSector3", "ZNA Signal Sector 3", kTH1F, {axisZDCEnergy});
191-
histos.add("QA/hZNASignalSector4", "ZNA Signal Sector 4", kTH1F, {axisZDCEnergy});
192-
histos.add("QA/hZNCSignalSector1", "ZNC Signal Sector 1", kTH1F, {axisZDCEnergy});
193-
histos.add("QA/hZNCSignalSector2", "ZNC Signal Sector 2", kTH1F, {axisZDCEnergy});
194-
histos.add("QA/hZNCSignalSector3", "ZNC Signal Sector 3", kTH1F, {axisZDCEnergy});
195-
histos.add("QA/hZNCSignalSector4", "ZNC Signal Sector 4", kTH1F, {axisZDCEnergy});
178+
histos.add("QA/GainCalib/hZNASignal", "ZNA Signal", kTH2F, {{4, 0, 4}, {axisZDCEnergy}});
179+
histos.add("QA/GainCalib/hZNCSignal", "ZNC Signal", kTH2F, {{4, 0, 4}, {axisZDCEnergy}});
180+
histos.add("QA/hZNASignal", "ZNA Signal", kTProfile2D, {{4, 0, 4}, {axisVz}});
181+
histos.add("QA/hZNCSignal", "ZNC Signal", kTProfile2D, {{4, 0, 4}, {axisVz}});
182+
histos.add("QA/hZNAEnergyCommon", "ZNA Energy Common", kTProfile, {axisVz});
183+
histos.add("QA/hZNCEnergyCommon", "ZNC Energy Common", kTProfile, {axisVz});
196184
histos.add("CorrHist/hWtXZNA", "X^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
197185
histos.add("CorrHist/hWtYZNA", "Y^{ZNA}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
198186
histos.add("CorrHist/hWtXZNC", "X^{ZNC}_{1}", kTHnSparseF, {axisCoarseCent, axisCoarseVx, axisCoarseVy, axisCoarseVz});
@@ -221,17 +209,19 @@ struct FlowEventPlane {
221209
histos.add("Checks/hPsiSPC", "#Psi_{SP}^{C} distribution", kTH2F, {axisCent, axisPsi});
222210
histos.add("Checks/hCosPsiSPAC", "Cos(#Psi_{SP}^{A} #minus #Psi_{SP}^{C}) distribution", kTProfile, {axisCent});
223211
histos.add("Checks/hSinPsiSPAC", "Sin(#Psi_{SP}^{A} #minus #Psi_{SP}^{C}) distribution", kTProfile, {axisCent});
224-
histos.add("Checks/hXaXc", "X^{ZNC}_{1} Vs X^{ZNA}_{1}", kTProfile, {axisCent});
225-
histos.add("Checks/hYaYc", "Y^{ZNC}_{1} Vs Y^{ZNA}_{1}", kTProfile, {axisCent});
212+
histos.add("Checks/hXaXc", "X^{A}_{1}X^{C}_{1}", kTProfile, {axisCent});
213+
histos.add("Checks/hYaYc", "Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
214+
histos.add("Checks/hXaYc", "X^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
215+
histos.add("Checks/hYaXc", "Y^{A}_{1}X^{C}_{1}", kTProfile, {axisCent});
226216
histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY});
227217
histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ});
228218
histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
229-
histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
230-
histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
231-
histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
232-
histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
233-
histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
234-
histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTH3F, {axisCent, axisV1, axisTrackEta});
219+
histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
220+
histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
221+
histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
222+
histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
223+
histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
224+
histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
235225
}
236226

237227
template <typename C>
@@ -307,6 +297,23 @@ struct FlowEventPlane {
307297
return true;
308298
}
309299

300+
void gainCalib(float const& vz, int const& runNumber, std::array<float, 4>& energy, const char* histName = "hZNAEnergy")
301+
{
302+
// Set CCDB path
303+
std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/GainCalib" + "/Run" + std::to_string(runNumber);
304+
305+
// Get object from CCDB
306+
auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
307+
308+
// Get gain calibration
309+
TH2F* h = reinterpret_cast<TH2F*>(ccdbObj->FindObject(histName));
310+
float v = 0.;
311+
for (int i = 0; i < static_cast<int>(energy.size()); ++i) {
312+
v = h->GetBinContent(h->FindBin(i + 0.5, vz + 0.00001));
313+
energy[i] *= v;
314+
}
315+
}
316+
310317
template <typename C, typename F>
311318
std::vector<float> getAvgCorrFactors(const C& ccdbObject, const F& corrType, const std::vector<float>& vCollParam)
312319
{
@@ -317,18 +324,15 @@ struct FlowEventPlane {
317324
for (auto const& x : vHistNames) {
318325
int cntry = 0;
319326
for (auto const& y : x) {
320-
ccdbObjects.obj = reinterpret_cast<TObject*>(ccdbObject->FindObject(y.c_str()));
321327
if (corrType == kFineCorr) {
322-
ccdbObjects.hp = reinterpret_cast<TProfile*>(ccdbObjects.obj->Clone());
323-
vAvgOutput[cntrx] += ccdbObjects.hp->GetBinContent(ccdbObjects.hp->GetXaxis()->FindBin(vCollParam[cntry]));
324-
delete ccdbObjects.hp;
328+
TProfile* hp = reinterpret_cast<TProfile*>(ccdbObject->FindObject(y.c_str()));
329+
vAvgOutput[cntrx] += hp->GetBinContent(hp->GetXaxis()->FindBin(vCollParam[cntry]));
325330
} else {
326-
ccdbObjects.hn = reinterpret_cast<THnSparseF*>(ccdbObjects.obj->Clone());
331+
THnSparseF* hn = reinterpret_cast<THnSparseF*>(ccdbObject->FindObject(y.c_str()));
327332
for (int i = 0; i < static_cast<int>(vHistNames.size()); ++i) {
328-
binarray[i] = ccdbObjects.hn->GetAxis(i)->FindBin(vCollParam[i]);
333+
binarray[i] = hn->GetAxis(i)->FindBin(vCollParam[i]);
329334
}
330-
vAvgOutput[cntrx] += ccdbObjects.hn->GetBinContent(ccdbObjects.hn->GetBin(binarray));
331-
delete ccdbObjects.hn;
335+
vAvgOutput[cntrx] += hn->GetBinContent(hn->GetBin(binarray));
332336
}
333337
++cntry;
334338
}
@@ -337,7 +341,7 @@ struct FlowEventPlane {
337341
return vAvgOutput;
338342
}
339343

340-
void applyCorrection(const std::vector<float> inputParam, const int& runNumber, std::vector<float>& outputParam)
344+
void applyCorrection(const std::vector<float>& inputParam, const int& runNumber, std::vector<float>& outputParam)
341345
{
342346
std::vector<int> vCorrFlags = static_cast<std::vector<int>>(cCorrFlagVector);
343347
int nitr = vCorrFlags.size();
@@ -361,16 +365,16 @@ struct FlowEventPlane {
361365
ccdbPath = static_cast<std::string>(cCcdbPath) + "/CorrItr_" + std::to_string(i + 1) + "/Run" + std::to_string(runNumber);
362366

363367
// Get object from CCDB
364-
ccdbObjects.ccdbList = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
368+
auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
365369

366370
// Check CCDB Object
367-
if (!ccdbObjects.ccdbList) {
371+
if (!ccdbObj) {
368372
LOGF(warning, "CCDB OBJECT NOT FOUND");
369373
return;
370374
}
371375

372376
// Get averages
373-
std::vector<float> vAvg = getAvgCorrFactors(ccdbObjects.ccdbList, corrType, inputParam);
377+
std::vector<float> vAvg = getAvgCorrFactors(ccdbObj, corrType, inputParam);
374378

375379
// Apply correction
376380
outputParam[kXa] -= vAvg[kXa];
@@ -437,6 +441,7 @@ struct FlowEventPlane {
437441

438442
// Get bunch crossing
439443
auto bc = collision.foundBC_as<BCsRun3>();
444+
int runNumber = collision.foundBC_as<BCsRun3>().runNumber();
440445

441446
// check zdc
442447
if (!bc.has_zdc()) {
@@ -454,59 +459,82 @@ struct FlowEventPlane {
454459
return;
455460
}
456461

457-
// Fill QA histograms
458-
histos.fill(HIST("QA/hZNASignalSector1"), znaEnergy[0]);
459-
histos.fill(HIST("QA/hZNASignalSector2"), znaEnergy[1]);
460-
histos.fill(HIST("QA/hZNASignalSector3"), znaEnergy[2]);
461-
histos.fill(HIST("QA/hZNASignalSector4"), znaEnergy[3]);
462-
histos.fill(HIST("QA/hZNCSignalSector1"), zncEnergy[0]);
463-
histos.fill(HIST("QA/hZNCSignalSector2"), zncEnergy[1]);
464-
histos.fill(HIST("QA/hZNCSignalSector3"), zncEnergy[2]);
465-
histos.fill(HIST("QA/hZNCSignalSector4"), zncEnergy[3]);
466-
467-
/*auto alphaZDC = 0.395;*/
462+
// Fill gain calib histograms
463+
for (int iCh = 0; iCh < kXYAC; ++iCh) {
464+
histos.fill(HIST("QA/hZNASignal"), iCh + 0.5, vCollParam[kVz], znaEnergy[iCh]);
465+
histos.fill(HIST("QA/hZNCSignal"), iCh + 0.5, vCollParam[kVz], zncEnergy[iCh]);
466+
histos.fill(HIST("QA/hZNAEnergyCommon"), vCollParam[kVz], znaEnergyCommon);
467+
histos.fill(HIST("QA/hZNCEnergyCommon"), vCollParam[kVz], zncEnergyCommon);
468+
}
469+
470+
// Do gain calibration
471+
if (cDoGainCalib) {
472+
gainCalib(vCollParam[kVz], runNumber, znaEnergy, "hZNASignal");
473+
gainCalib(vCollParam[kVz], runNumber, zncEnergy, "hZNCSignal");
474+
}
475+
476+
// Fill zdc signal
477+
for (int iCh = 0; iCh < kXYAC; ++iCh) {
478+
histos.fill(HIST("QA/GainCalib/hZNASignal"), iCh + 0.5, znaEnergy[iCh]);
479+
histos.fill(HIST("QA/GainCalib/hZNCSignal"), iCh + 0.5, zncEnergy[iCh]);
480+
}
481+
482+
auto alphaZDC = 0.395;
468483
const double x[4] = {-1.75, 1.75, -1.75, 1.75};
469484
const double y[4] = {-1.75, -1.75, 1.75, 1.75};
470485

471486
// Calculate X and Y
472-
float znaXSumNum = 0., znaXSumDnm = 0.;
473-
float znaYSumNum = 0., znaYSumDnm = 0.;
474-
float zncXSumNum = 0., zncXSumDnm = 0.;
475-
float zncYSumNum = 0., zncYSumDnm = 0.;
487+
float znaXNum = 0., znaYNum = 0., zncXNum = 0., zncYNum = 0.;
488+
float znaDen = 0., zncDen = 0.;
489+
float znaWt = 0., zncWt = 0.;
476490

477491
// Loop over zdc sectors
478492
for (int i = 0; i < kXYAC; ++i) {
479-
znaXSumNum += znaXWeigthEnergy[i] * znaEnergy[i] * x[i];
480-
znaYSumNum += znaYWeigthEnergy[i] * znaEnergy[i] * y[i];
481-
znaXSumDnm += znaXWeigthEnergy[i] * znaEnergy[i];
482-
znaYSumDnm += znaYWeigthEnergy[i] * znaEnergy[i];
483-
zncXSumNum += zncXWeigthEnergy[i] * zncEnergy[i] * x[i];
484-
zncYSumNum += zncYWeigthEnergy[i] * zncEnergy[i] * y[i];
485-
zncXSumDnm += zncXWeigthEnergy[i] * zncEnergy[i];
486-
zncYSumDnm += zncYWeigthEnergy[i] * zncEnergy[i];
493+
if (cUseAlphaZDC) {
494+
znaWt = std::pow(znaEnergy[i], alphaZDC);
495+
zncWt = std::pow(zncEnergy[i], alphaZDC);
496+
} else {
497+
znaWt = znaEnergy[i];
498+
zncWt = zncEnergy[i];
499+
}
500+
znaXNum -= znaWt * x[i];
501+
znaYNum += znaWt * y[i];
502+
zncXNum += zncWt * x[i];
503+
zncYNum += zncWt * y[i];
504+
znaDen += znaWt;
505+
zncDen += zncWt;
506+
}
507+
508+
if (znaDen < zdcDenThrs || zncDen < zdcDenThrs) {
509+
return;
487510
}
488511

489512
// Get X and Y for A and C side ZNA
490513
std::vector<float> vSP = {0, 0, 0, 0};
491-
vSP[kXa] = znaXSumNum / znaXSumDnm;
492-
vSP[kYa] = znaYSumNum / znaYSumDnm;
493-
vSP[kXc] = zncXSumNum / zncXSumDnm;
494-
vSP[kYc] = zncYSumNum / zncYSumDnm;
514+
vSP[kXa] = znaXNum / znaDen;
515+
vSP[kYa] = znaYNum / znaDen;
516+
vSP[kXc] = zncXNum / zncDen;
517+
vSP[kYc] = zncYNum / zncDen;
495518

496519
// Do corrections
497-
int runNumber = collision.foundBC_as<BCsRun3>().runNumber();
498-
applyCorrection(vCollParam, runNumber, vSP);
520+
if (cApplyRecentCorr) {
521+
applyCorrection(vCollParam, runNumber, vSP);
522+
}
499523

500-
// Fill X and Y histograms
524+
// Fill X and Y histograms for corrections after each iteration
501525
fillCorrHist(vCollParam, vSP);
526+
527+
// Evaluate spectator plane angle and [X,Y] correlations
502528
float psiA = std::atan2(vSP[kYa], vSP[kXa]);
503-
float psiC = std::atan2(vSP[kYa], vSP[kXa]);
504-
histos.fill(HIST("Checks/hXaXc"), cent, (vSP[kXa] * vSP[kXc]));
505-
histos.fill(HIST("Checks/hYaYc"), cent, (vSP[kYa] * vSP[kYc]));
529+
float psiC = std::atan2(vSP[kYc], vSP[kXc]);
506530
histos.fill(HIST("Checks/hPsiSPA"), cent, psiA);
507531
histos.fill(HIST("Checks/hPsiSPC"), cent, psiC);
508532
histos.fill(HIST("Checks/hCosPsiSPAC"), cent, std::cos(psiA - psiC));
509533
histos.fill(HIST("Checks/hSinPsiSPAC"), cent, std::sin(psiA - psiC));
534+
histos.fill(HIST("Checks/hXaXc"), cent, (vSP[kXa] * vSP[kXc]));
535+
histos.fill(HIST("Checks/hYaYc"), cent, (vSP[kYa] * vSP[kYc]));
536+
histos.fill(HIST("Checks/hXaYc"), cent, (vSP[kXa] * vSP[kYc]));
537+
histos.fill(HIST("Checks/hYaXc"), cent, (vSP[kYa] * vSP[kXc]));
510538

511539
// Directed flow
512540
float qac = (vSP[kXa] * vSP[kXc]) + (vSP[kYa] * vSP[kYc]);
@@ -530,14 +558,14 @@ struct FlowEventPlane {
530558
v1c = ux * vSP[kXc] + uy * vSP[kYc];
531559

532560
// Fill histogram
533-
histos.fill(HIST("DF/hAQu"), cent, v1a, track.eta());
534-
histos.fill(HIST("DF/hCQu"), cent, v1c, track.eta());
561+
histos.fill(HIST("DF/hAQu"), cent, track.eta(), v1a);
562+
histos.fill(HIST("DF/hCQu"), cent, track.eta(), v1c);
535563
if (track.sign() > 0) {
536-
histos.fill(HIST("DF/hAQuPos"), cent, v1a, track.eta());
537-
histos.fill(HIST("DF/hCQuPos"), cent, v1c, track.eta());
564+
histos.fill(HIST("DF/hAQuPos"), cent, track.eta(), v1a);
565+
histos.fill(HIST("DF/hCQuPos"), cent, track.eta(), v1c);
538566
} else {
539-
histos.fill(HIST("DF/hAQuNeg"), cent, v1a, track.eta());
540-
histos.fill(HIST("DF/hCQuNeg"), cent, v1c, track.eta());
567+
histos.fill(HIST("DF/hAQuNeg"), cent, track.eta(), v1a);
568+
histos.fill(HIST("DF/hCQuNeg"), cent, track.eta(), v1c);
541569
}
542570
}
543571
}

0 commit comments

Comments
 (0)