Skip to content

Commit 04a0040

Browse files
jikim1290alibuild
jikim1290
and
alibuild
authored
[PWGCF] adding shift correction (#13297)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent c4d9a50 commit 04a0040

File tree

1 file changed

+134
-43
lines changed

1 file changed

+134
-43
lines changed

‎PWGCF/JCorran/Tasks/jEPFlowAnalysis.cxx‎

Lines changed: 134 additions & 43 deletions
Original file line numberDiff line numberDiff line change
@@ -8,64 +8,84 @@
88
// In applying this license CERN does not waive the privileges and immunities
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
11+
1112
/// \author Maxim Virta (maxim.virta@cern.ch)
13+
/// \brief flow measurement with q-vectors
14+
/// \file jEPFlowAnalysis.cxx
1215
/// \since Jul 2024
1316

14-
#include "Framework/AnalysisTask.h"
15-
#include "Framework/RunningWorkflowInfo.h"
16-
#include "Framework/HistogramRegistry.h"
17+
#include "FlowJHistManager.h"
1718

18-
#include "Common/DataModel/EventSelection.h"
19+
#include "Common/Core/EventPlaneHelper.h"
1920
#include "Common/Core/TrackSelection.h"
20-
#include "Framework/runDataProcessing.h"
21+
#include "Common/DataModel/EventSelection.h"
22+
#include "Common/DataModel/Qvectors.h"
2123
#include "Common/DataModel/TrackSelectionTables.h"
2224

23-
#include "Common/DataModel/Qvectors.h"
24-
#include "Common/Core/EventPlaneHelper.h"
25+
#include "CCDB/BasicCCDBManager.h"
26+
#include "CCDB/CcdbApi.h"
27+
#include "Framework/AnalysisTask.h"
28+
#include "Framework/HistogramRegistry.h"
29+
#include "Framework/RunningWorkflowInfo.h"
30+
#include "Framework/runDataProcessing.h"
2531

26-
#include "FlowJHistManager.h"
27-
#include "JEPFlowAnalysis.h"
32+
#include <string>
33+
#include <vector>
2834

2935
using namespace o2;
3036
using namespace o2::framework;
3137
using namespace o2::framework::expressions;
3238
using namespace std;
3339

3440
using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Qvectors>;
35-
3641
using MyTracks = aod::Tracks;
3742

3843
struct jEPFlowAnalysis {
3944

40-
HistogramRegistry EPFlowHistograms{"EPFlow", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
41-
JEPFlowAnalysis epAnalysis;
45+
HistogramRegistry epFlowHistograms{"EPFlow", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
4246
EventPlaneHelper helperEP;
4347
FlowJHistManager histManager;
44-
Bool_t debug = kFALSE;
48+
bool debug = kFALSE;
49+
Service<o2::ccdb::BasicCCDBManager> ccdb;
50+
o2::ccdb::CcdbApi ccdbApi;
4551

4652
// Set Configurables here
4753
struct : ConfigurableGroup {
4854
Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "Minimum pT used for track selection."};
49-
Configurable<float> cfgPtMax{"cfgPtMax", 5.0f, "Maximum pT used for track selection."};
5055
Configurable<float> cfgEtaMax{"cfgEtaMax", 1.f, "Maximum eta used for track selection."};
5156
} cfgTrackCuts;
5257

5358
Configurable<bool> cfgAddEvtSel{"cfgAddEvtSel", true, "Use event selection"};
5459
Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
60+
Configurable<int> cfgnMode{"cfgnMode", 1, "the number of modulations"};
61+
62+
Configurable<bool> cfgShiftCorr{"cfgShiftCorr", false, "additional shift correction"};
63+
Configurable<std::string> cfgShiftPath{"cfgShiftPath", "Users/j/junlee/Qvector/QvecCalib/Shift", "Path for Shift"};
64+
Configurable<bool> cfgSPmethod{"cfgSPmethod", false, "flag for scalar product"};
5565

5666
Configurable<std::string> cfgDetName{"cfgDetName", "FT0C", "The name of detector to be analyzed"};
5767
Configurable<std::string> cfgRefAName{"cfgRefAName", "TPCPos", "The name of detector for reference A"};
5868
Configurable<std::string> cfgRefBName{"cfgRefBName", "TPCNeg", "The name of detector for reference B"};
5969

60-
Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (aod::track::pt < cfgTrackCuts.cfgPtMax) && (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax);
70+
ConfigurableAxis cfgAxisCent{"cfgAxisCent", {100, 0, 100}, ""};
71+
ConfigurableAxis cfgAxisPt{"cfgAxisPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 10.0, 12.0, 15.0, 30.0, 50.0, 70.0, 100.0}, ""};
72+
ConfigurableAxis cfgAxisCos{"cfgAxisCos", {102, -1.02, 1.02}, ""};
73+
74+
Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax);
6175

62-
int DetId;
63-
int RefAId;
64-
int RefBId;
76+
int detId;
77+
int refAId;
78+
int refBId;
6579
int harmInd;
6680

81+
int currentRunNumber = -999;
82+
int lastRunNumber = -999;
83+
84+
std::vector<TProfile3D*> shiftprofile{};
85+
std::string fullCCDBShiftCorrPath;
86+
6787
template <typename T>
68-
int GetDetId(const T& name)
88+
int getdetId(const T& name)
6989
{
7090
if (name.value == "FT0C") {
7191
return 0;
@@ -88,37 +108,108 @@ struct jEPFlowAnalysis {
88108

89109
void init(InitContext const&)
90110
{
91-
DetId = GetDetId(cfgDetName);
92-
RefAId = GetDetId(cfgRefAName);
93-
RefBId = GetDetId(cfgRefBName);
111+
detId = getdetId(cfgDetName);
112+
refAId = getdetId(cfgRefAName);
113+
refBId = getdetId(cfgRefBName);
114+
115+
AxisSpec axisMod{cfgnMode, 2., cfgnMode + 2.};
116+
AxisSpec axisEvtPl{360, -constants::math::PI * 1.1, constants::math::PI * 1.1};
117+
118+
AxisSpec axisCent{cfgAxisCent, "cent"};
119+
AxisSpec axisPt{cfgAxisPt, "pT"};
120+
AxisSpec axisCos{cfgAxisCos, "cos"};
121+
122+
epFlowHistograms.add("EpDet", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
123+
epFlowHistograms.add("EpRefA", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
124+
epFlowHistograms.add("EpRefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
94125

95-
epAnalysis.SetHistRegistry(&EPFlowHistograms);
96-
epAnalysis.CreateHistograms();
126+
epFlowHistograms.add("EpResDetRefA", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
127+
epFlowHistograms.add("EpResDetRefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
128+
epFlowHistograms.add("EpResRefARefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
129+
130+
epFlowHistograms.add("vncos", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisCos}});
131+
epFlowHistograms.add("vnsin", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisCos}});
97132
}
98133

99-
void process(MyCollisions::iterator const& coll, soa::Filtered<MyTracks> const& tracks)
134+
void process(MyCollisions::iterator const& coll, soa::Filtered<MyTracks> const& tracks, aod::BCsWithTimestamps const&)
100135
{
101136
if (cfgAddEvtSel && (!coll.sel8() || !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !coll.selection_bit(aod::evsel::kNoSameBunchPileup)))
102137
return;
103138

104-
Float_t cent = coll.cent();
105-
EPFlowHistograms.fill(HIST("FullCentrality"), cent);
106-
Float_t EPs[3] = {0.};
107-
for (int i = 2; i < 5; i++) { // loop over different harmonic orders
108-
harmInd = cfgnTotalSystem * 4 * (i - 2) + 3; // harmonic index to access corresponding Q-vector as all Q-vectors are in same vector
109-
EPs[0] = helperEP.GetEventPlane(coll.qvecRe()[DetId + harmInd], coll.qvecIm()[DetId + harmInd], i);
110-
EPs[1] = helperEP.GetEventPlane(coll.qvecRe()[RefAId + harmInd], coll.qvecIm()[RefAId + harmInd], i);
111-
EPs[2] = helperEP.GetEventPlane(coll.qvecRe()[RefBId + harmInd], coll.qvecIm()[RefBId + harmInd], i);
112-
113-
Float_t resNumA = helperEP.GetResolution(EPs[0], EPs[1], i);
114-
Float_t resNumB = helperEP.GetResolution(EPs[0], EPs[2], i);
115-
Float_t resDenom = helperEP.GetResolution(EPs[1], EPs[2], i);
116-
epAnalysis.FillResolutionHistograms(cent, static_cast<float>(i), resNumA, resNumB, resDenom);
117-
for (uint j = 0; j < 3; j++) { // loop over detectors used
118-
for (auto& track : tracks) {
119-
Float_t vn = TMath::Cos((i) * (track.phi() - EPs[j]));
120-
Float_t vn_sin = TMath::Sin((i) * (track.phi() - EPs[j]));
121-
epAnalysis.FillVnHistograms(i, cent, static_cast<float>(j + 1), track.pt(), vn, vn_sin);
139+
float cent = coll.cent();
140+
epFlowHistograms.fill(HIST("FullCentrality"), cent);
141+
float eps[3] = {0.};
142+
143+
if (cfgShiftCorr) {
144+
auto bc = coll.bc_as<aod::BCsWithTimestamps>();
145+
currentRunNumber = bc.runNumber();
146+
if (currentRunNumber != lastRunNumber) {
147+
shiftprofile.clear();
148+
for (int i = 0; i < cfgnMode; i++) {
149+
fullCCDBShiftCorrPath = cfgShiftPath;
150+
fullCCDBShiftCorrPath += "/v";
151+
fullCCDBShiftCorrPath += std::to_string(i + 2);
152+
auto objshift = ccdb->getForTimeStamp<TProfile3D>(fullCCDBShiftCorrPath, bc.timestamp());
153+
shiftprofile.push_back(objshift);
154+
}
155+
lastRunNumber = currentRunNumber;
156+
}
157+
}
158+
159+
for (int i = 0; i < cfgnMode; i++) { // loop over different harmonic orders
160+
harmInd = cfgnTotalSystem * 4 * (i) + 3; // harmonic index to access corresponding Q-vector as all Q-vectors are in same vector
161+
eps[0] = helperEP.GetEventPlane(coll.qvecRe()[detId + harmInd], coll.qvecIm()[detId + harmInd], i + 2);
162+
eps[1] = helperEP.GetEventPlane(coll.qvecRe()[refAId + harmInd], coll.qvecIm()[refAId + harmInd], i + 2);
163+
eps[2] = helperEP.GetEventPlane(coll.qvecRe()[refBId + harmInd], coll.qvecIm()[refBId + harmInd], i + 2);
164+
165+
auto deltapsiDet = 0.0;
166+
auto deltapsiRefA = 0.0;
167+
auto deltapsiRefB = 0.0;
168+
169+
float weight = 1.0;
170+
171+
if (cfgShiftCorr) {
172+
constexpr int kShiftBins = 10;
173+
for (int ishift = 1; ishift <= kShiftBins; ishift++) {
174+
auto coeffshiftxDet = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 0.5, ishift - 0.5));
175+
auto coeffshiftyDet = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 1.5, ishift - 0.5));
176+
auto coeffshiftxRefA = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 2.5, ishift - 0.5));
177+
auto coeffshiftyRefA = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 3.5, ishift - 0.5));
178+
auto coeffshiftxRefB = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 4.5, ishift - 0.5));
179+
auto coeffshiftyRefB = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 5.5, ishift - 0.5)); // currently only FT0C/TPCpos/TPCneg
180+
181+
deltapsiDet += ((1 / (1.0 * ishift)) * (-coeffshiftxDet * std::cos(ishift * static_cast<float>(i + 2) * eps[0]) + coeffshiftyDet * std::sin(ishift * static_cast<float>(i + 2) * eps[0])));
182+
deltapsiRefA += ((1 / (1.0 * ishift)) * (-coeffshiftxRefA * std::cos(ishift * static_cast<float>(i + 2) * eps[1]) + coeffshiftyRefA * std::sin(ishift * static_cast<float>(i + 2) * eps[1])));
183+
deltapsiRefB += ((1 / (1.0 * ishift)) * (-coeffshiftxRefB * std::cos(ishift * static_cast<float>(i + 2) * eps[2]) + coeffshiftyRefB * std::sin(ishift * static_cast<float>(i + 2) * eps[2])));
184+
}
185+
186+
eps[0] += deltapsiDet;
187+
eps[1] += deltapsiRefA;
188+
eps[2] += deltapsiRefB;
189+
}
190+
191+
if (cfgSPmethod)
192+
weight *= std::sqrt(std::pow(coll.qvecRe()[detId + harmInd], 2) + std::pow(coll.qvecIm()[detId + harmInd], 2));
193+
194+
float resNumA = helperEP.GetResolution(eps[0], eps[1], i + 2);
195+
float resNumB = helperEP.GetResolution(eps[0], eps[2], i + 2);
196+
float resDenom = helperEP.GetResolution(eps[1], eps[2], i + 2);
197+
198+
epFlowHistograms.fill(HIST("EpDet"), i + 2, cent, eps[0]);
199+
epFlowHistograms.fill(HIST("EpRefA"), i + 2, cent, eps[1]);
200+
epFlowHistograms.fill(HIST("EpRefB"), i + 2, cent, eps[2]);
201+
202+
epFlowHistograms.fill(HIST("EpResDetRefA"), i + 2, cent, resNumA);
203+
epFlowHistograms.fill(HIST("EpResDetRefB"), i + 2, cent, resNumB);
204+
epFlowHistograms.fill(HIST("EpResRefARefB"), i + 2, cent, resDenom);
205+
206+
for (int j = 0; j < cfgnMode; j++) { // loop over detectors used
207+
for (const auto& track : tracks) {
208+
float vn = std::cos((i + 2) * (track.phi() - eps[j]));
209+
float vnSin = std::sin((i + 2) * (track.phi() - eps[j]));
210+
211+
epFlowHistograms.fill(HIST("vncos"), i + 2, cent, track.pt(), vn * weight);
212+
epFlowHistograms.fill(HIST("vnsin"), i + 2, cent, track.pt(), vnSin * weight);
122213
}
123214
}
124215
}

0 commit comments

Comments
 (0)