Skip to content

Commit 7698f3b

Browse files
jikim1290jikim1290
committed
adding sp method
1 parent d840355 commit 7698f3b

File tree

1 file changed

+289
-0
lines changed

1 file changed

+289
-0
lines changed

‎PWGCF/JCorran/jEPFlowAnalysis.cxx‎

Lines changed: 289 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,289 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
/// \author Maxim Virta (maxim.virta@cern.ch)
13+
/// \brief flow measurement with q-vectors
14+
/// \file jEPFlowAnalysis.cxx
15+
/// \since Jul 2024
16+
17+
#include "FlowJHistManager.h"
18+
19+
#include "Common/Core/EventPlaneHelper.h"
20+
#include "Common/Core/TrackSelection.h"
21+
#include "Common/DataModel/EventSelection.h"
22+
#include "Common/DataModel/Qvectors.h"
23+
#include "Common/DataModel/TrackSelectionTables.h"
24+
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"
31+
32+
#include <string>
33+
#include <vector>
34+
35+
using namespace o2;
36+
using namespace o2::framework;
37+
using namespace o2::framework::expressions;
38+
using namespace std;
39+
40+
using MyCollisions = soa::Join<aod::Collisions, aod::EvSels, aod::Qvectors>;
41+
using MyTracks = aod::Tracks;
42+
43+
struct jEPFlowAnalysis {
44+
45+
HistogramRegistry epFlowHistograms{"EPFlow", {}, OutputObjHandlingPolicy::AnalysisObject, true, true};
46+
EventPlaneHelper helperEP;
47+
FlowJHistManager histManager;
48+
bool debug = kFALSE;
49+
Service<o2::ccdb::BasicCCDBManager> ccdb;
50+
o2::ccdb::CcdbApi ccdbApi;
51+
52+
// Set Configurables here
53+
struct : ConfigurableGroup {
54+
Configurable<float> cfgPtMin{"cfgPtMin", 0.2f, "Minimum pT used for track selection."};
55+
Configurable<float> cfgEtaMax{"cfgEtaMax", 1.f, "Maximum eta used for track selection."};
56+
} cfgTrackCuts;
57+
58+
Configurable<bool> cfgAddEvtSel{"cfgAddEvtSel", true, "event selection"};
59+
Configurable<int> cfgEvtSel{"cfgEvtSel", 0, "Event selection flags\n0: Sel8\n1: Sel8+kIsGoodZvtxFT0vsPV+kNoSameBunchPileup\n2: Sel8+kIsGoodZvtxFT0vsPV+kNoSameBunchPileup+kNoCollInTimeRangeStandard\n3: Sel8+kNoSameBunchPileup"};
60+
Configurable<int> cfgMaxOccupancy{"cfgMaxOccupancy", 999999, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
61+
Configurable<int> cfgMinOccupancy{"cfgMinOccupancy", 0, "maximum occupancy of tracks in neighbouring collisions in a given time range"};
62+
63+
Configurable<int> cfgnTotalSystem{"cfgnTotalSystem", 7, "Total number of detectors in qVectorsTable"};
64+
Configurable<int> cfgnMode{"cfgnMode", 1, "the number of modulations"};
65+
66+
Configurable<bool> cfgShiftCorr{"cfgShiftCorr", false, "additional shift correction"};
67+
Configurable<std::string> cfgShiftPath{"cfgShiftPath", "Users/j/junlee/Qvector/QvecCalib/Shift", "Path for Shift"};
68+
Configurable<float> cfgVertexZ{"cfgVertexZ", 10.0, "Maximum vertex Z selection"};
69+
70+
Configurable<std::string> cfgDetName{"cfgDetName", "FT0C", "The name of detector to be analyzed"};
71+
Configurable<std::string> cfgRefAName{"cfgRefAName", "TPCPos", "The name of detector for reference A"};
72+
Configurable<std::string> cfgRefBName{"cfgRefBName", "TPCNeg", "The name of detector for reference B"};
73+
74+
ConfigurableAxis cfgAxisCent{"cfgAxisCent", {100, 0, 100}, ""};
75+
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}, ""};
76+
ConfigurableAxis cfgAxisCos{"cfgAxisCos", {102, -1.02, 1.02}, ""};
77+
ConfigurableAxis cfgAxisQvec{"cfgAxisQvec", {200, -5.0, 5.0}, ""};
78+
79+
Filter trackFilter = (aod::track::pt > cfgTrackCuts.cfgPtMin) && (nabs(aod::track::eta) < cfgTrackCuts.cfgEtaMax);
80+
81+
int detId;
82+
int refAId;
83+
int refBId;
84+
int harmInd;
85+
86+
int currentRunNumber = -999;
87+
int lastRunNumber = -999;
88+
89+
std::vector<TProfile3D*> shiftprofile{};
90+
std::string fullCCDBShiftCorrPath;
91+
92+
template <typename T>
93+
int getdetId(const T& name)
94+
{
95+
if (name.value == "FT0C") {
96+
return 0;
97+
} else if (name.value == "FT0A") {
98+
return 1;
99+
} else if (name.value == "FT0M") {
100+
return 2;
101+
} else if (name.value == "FV0A") {
102+
return 3;
103+
} else if (name.value == "TPCPos") {
104+
return 4;
105+
} else if (name.value == "TPCNeg") {
106+
return 5;
107+
} else if (name.value == "TPCTot") {
108+
return 6;
109+
} else {
110+
return 0;
111+
}
112+
}
113+
114+
void init(InitContext const&)
115+
{
116+
detId = getdetId(cfgDetName);
117+
refAId = getdetId(cfgRefAName);
118+
refBId = getdetId(cfgRefBName);
119+
120+
AxisSpec axisMod{cfgnMode, 2., cfgnMode + 2.};
121+
AxisSpec axisEvtPl{360, -constants::math::PI * 1.1, constants::math::PI * 1.1};
122+
AxisSpec axisVertex{150, -12.5, 12.5};
123+
124+
AxisSpec axisCent{cfgAxisCent, "cent"};
125+
AxisSpec axisPt{cfgAxisPt, "pT"};
126+
AxisSpec axisCos{cfgAxisCos, "cos"};
127+
AxisSpec axisQvec{cfgAxisQvec, "Qvec"};
128+
129+
epFlowHistograms.add("EpDet", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
130+
epFlowHistograms.add("EpRefA", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
131+
epFlowHistograms.add("EpRefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
132+
133+
epFlowHistograms.add("EpResDetRefA", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
134+
epFlowHistograms.add("EpResDetRefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
135+
epFlowHistograms.add("EpResRefARefB", "", {HistType::kTH3F, {axisMod, axisCent, axisEvtPl}});
136+
137+
epFlowHistograms.add("vncos", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisCos}});
138+
epFlowHistograms.add("vnsin", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisCos}});
139+
140+
epFlowHistograms.add("EpResQvecDetRefAxx", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
141+
epFlowHistograms.add("EpResQvecDetRefAxy", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
142+
epFlowHistograms.add("EpResQvecDetRefBxx", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
143+
epFlowHistograms.add("EpResQvecDetRefBxy", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
144+
epFlowHistograms.add("EpResQvecRefARefBxx", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
145+
epFlowHistograms.add("EpResQvecRefARefBxy", "", {HistType::kTH3F, {axisMod, axisCent, axisQvec}});
146+
147+
epFlowHistograms.add("SPvnxx", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisQvec}});
148+
epFlowHistograms.add("SPvnxy", "", {HistType::kTHnSparseF, {axisMod, axisCent, axisPt, axisQvec}});
149+
150+
epFlowHistograms.add("hCentrality", "", {HistType::kTH1F, {axisCent}});
151+
epFlowHistograms.add("hVertex", "", {HistType::kTH1F, {axisVertex}});
152+
}
153+
154+
void process(MyCollisions::iterator const& coll, soa::Filtered<MyTracks> const& tracks, aod::BCsWithTimestamps const&)
155+
{
156+
if (cfgAddEvtSel) {
157+
if (std::abs(coll.posZ()) > cfgVertexZ)
158+
return;
159+
switch (cfgEvtSel) {
160+
case 0: // Sel8
161+
if (!coll.sel8())
162+
return;
163+
break;
164+
case 1: // PbPb standard
165+
if (!coll.sel8() || !coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !coll.selection_bit(aod::evsel::kNoSameBunchPileup))
166+
return;
167+
break;
168+
case 2: // PbPb with pileup
169+
if (!coll.sel8() || !coll.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard) ||
170+
!coll.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV) || !coll.selection_bit(aod::evsel::kNoSameBunchPileup))
171+
return;
172+
break;
173+
case 3: // Small systems (OO, NeNe, pp)
174+
if (!coll.sel8() || !coll.selection_bit(aod::evsel::kNoSameBunchPileup))
175+
return;
176+
break;
177+
default:
178+
LOGF(warning, "Event selection flag was not found, continuing without basic event selections!\n");
179+
}
180+
// Check occupancy
181+
if (coll.trackOccupancyInTimeRange() > cfgMaxOccupancy || coll.trackOccupancyInTimeRange() < cfgMinOccupancy)
182+
return;
183+
}
184+
185+
float cent = coll.cent();
186+
epFlowHistograms.fill(HIST("hCentrality"), cent);
187+
epFlowHistograms.fill(HIST("hVertex"), coll.posZ());
188+
float eps[3] = {0.};
189+
float qx_shifted[3] = {0.};
190+
float qy_shifted[3] = {0.};
191+
192+
193+
if (cfgShiftCorr) {
194+
auto bc = coll.bc_as<aod::BCsWithTimestamps>();
195+
currentRunNumber = bc.runNumber();
196+
if (currentRunNumber != lastRunNumber) {
197+
shiftprofile.clear();
198+
for (int i = 0; i < cfgnMode; i++) {
199+
fullCCDBShiftCorrPath = cfgShiftPath;
200+
fullCCDBShiftCorrPath += "/v";
201+
fullCCDBShiftCorrPath += std::to_string(i + 2);
202+
auto objshift = ccdb->getForTimeStamp<TProfile3D>(fullCCDBShiftCorrPath, bc.timestamp());
203+
shiftprofile.push_back(objshift);
204+
}
205+
lastRunNumber = currentRunNumber;
206+
}
207+
}
208+
209+
if (coll.qvecAmp()[detId] < 1e-5 || coll.qvecAmp()[refAId] < 1e-5 || coll.qvecAmp()[refBId] < 1e-5)
210+
return;
211+
212+
for (int i = 0; i < cfgnMode; i++) { // loop over different harmonic orders
213+
harmInd = cfgnTotalSystem * 4 * (i) + 3; // harmonic index to access corresponding Q-vector as all Q-vectors are in same vector
214+
eps[0] = helperEP.GetEventPlane(coll.qvecRe()[4 * detId + harmInd], coll.qvecIm()[4 * detId + harmInd], i + 2);
215+
eps[1] = helperEP.GetEventPlane(coll.qvecRe()[4 * refAId + harmInd], coll.qvecIm()[4 * refAId + harmInd], i + 2);
216+
eps[2] = helperEP.GetEventPlane(coll.qvecRe()[4 * refBId + harmInd], coll.qvecIm()[4 * refBId + harmInd], i + 2);
217+
218+
auto deltapsiDet = 0.0;
219+
auto deltapsiRefA = 0.0;
220+
auto deltapsiRefB = 0.0;
221+
222+
float weight = 1.0;
223+
224+
if (cfgShiftCorr) {
225+
constexpr int kShiftBins = 10;
226+
for (int ishift = 1; ishift <= kShiftBins; ishift++) {
227+
auto coeffshiftxDet = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 0.5, ishift - 0.5));
228+
auto coeffshiftyDet = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 1.5, ishift - 0.5));
229+
auto coeffshiftxRefA = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 2.5, ishift - 0.5));
230+
auto coeffshiftyRefA = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 3.5, ishift - 0.5));
231+
auto coeffshiftxRefB = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 4.5, ishift - 0.5));
232+
auto coeffshiftyRefB = shiftprofile.at(i)->GetBinContent(shiftprofile.at(i)->FindBin(cent, 5.5, ishift - 0.5)); // currently only FT0C/TPCpos/TPCneg
233+
234+
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])));
235+
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])));
236+
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])));
237+
238+
}
239+
240+
eps[0] += deltapsiDet;
241+
eps[1] += deltapsiRefA;
242+
eps[2] += deltapsiRefB;
243+
244+
qx_shifted[0] = coll.qvecRe()[4 * detId + harmInd] * TMath::Cos(deltapsiDet) - coll.qvecIm()[4 * detId + harmInd] * TMath::Sin(deltapsiDet);
245+
qy_shifted[0] = coll.qvecRe()[4 * detId + harmInd] * TMath::Sin(deltapsiDet) + coll.qvecIm()[4 * detId + harmInd] * TMath::Cos(deltapsiDet);
246+
qx_shifted[1] = coll.qvecRe()[4 * refAId + harmInd] * TMath::Cos(deltapsiRefA) - coll.qvecIm()[4 * refAId + harmInd] * TMath::Sin(deltapsiRefA);
247+
qy_shifted[1] = coll.qvecRe()[4 * refAId + harmInd] * TMath::Sin(deltapsiRefA) + coll.qvecIm()[4 * refAId + harmInd] * TMath::Cos(deltapsiRefA);
248+
qx_shifted[2] = coll.qvecRe()[4 * refBId + harmInd] * TMath::Cos(deltapsiRefB) - coll.qvecIm()[4 * refBId + harmInd] * TMath::Sin(deltapsiRefB);
249+
qy_shifted[2] = coll.qvecRe()[4 * refBId + harmInd] * TMath::Sin(deltapsiRefB) + coll.qvecIm()[4 * refBId + harmInd] * TMath::Cos(deltapsiRefB);
250+
}
251+
252+
float resNumA = helperEP.GetResolution(eps[0], eps[1], i + 2);
253+
float resNumB = helperEP.GetResolution(eps[0], eps[2], i + 2);
254+
float resDenom = helperEP.GetResolution(eps[1], eps[2], i + 2);
255+
256+
epFlowHistograms.fill(HIST("EpDet"), i + 2, cent, eps[0]);
257+
epFlowHistograms.fill(HIST("EpRefA"), i + 2, cent, eps[1]);
258+
epFlowHistograms.fill(HIST("EpRefB"), i + 2, cent, eps[2]);
259+
260+
epFlowHistograms.fill(HIST("EpResDetRefA"), i + 2, cent, resNumA);
261+
epFlowHistograms.fill(HIST("EpResDetRefB"), i + 2, cent, resNumB);
262+
epFlowHistograms.fill(HIST("EpResRefARefB"), i + 2, cent, resDenom);
263+
264+
epFlowHistograms.fill(HIST("EpResQvecDetRefAxx"), i + 2, cent, qx_shifted[0] * qx_shifted[1] + qy_shifted[0] * qy_shifted[1]);
265+
epFlowHistograms.fill(HIST("EpResQvecDetRefAxy"), i + 2, cent, qx_shifted[1] * qy_shifted[0] - qx_shifted[0] * qy_shifted[1]);
266+
epFlowHistograms.fill(HIST("EpResQvecDetRefBxx"), i + 2, cent, qx_shifted[0] * qx_shifted[2] + qy_shifted[0] * qy_shifted[2]);
267+
epFlowHistograms.fill(HIST("EpResQvecDetRefBxy"), i + 2, cent, qx_shifted[2] * qy_shifted[0] - qx_shifted[0] * qy_shifted[2]);
268+
epFlowHistograms.fill(HIST("EpResQvecRefARefAxx"), i + 2, cent, qx_shifted[1] * qx_shifted[2] + qy_shifted[1] * qy_shifted[2]);
269+
epFlowHistograms.fill(HIST("EpResQvecRefARefAxy"), i + 2, cent, qx_shifted[2] * qy_shifted[1] - qx_shifted[1] * qy_shifted[2]);
270+
271+
for (const auto& track : tracks) {
272+
float vn = std::cos((i + 2) * (track.phi() - eps[0]));
273+
float vnSin = std::sin((i + 2) * (track.phi() - eps[0]));
274+
275+
epFlowHistograms.fill(HIST("vncos"), i + 2, cent, track.pt(), vn * weight);
276+
epFlowHistograms.fill(HIST("vnsin"), i + 2, cent, track.pt(), vnSin * weight);
277+
278+
epFlowHistograms.fill(HIST("SPvnxx"), i + 2, cent, track.pt(), (TMath::Cos(track.phi()) * qx_shifted[0] + TMath::Sin(track.phi()) * qy_shifted[0]) * weight);
279+
epFlowHistograms.fill(HIST("SPvnxy"), i + 2, cent, track.pt(), (TMath::Sin(track.phi()) * qx_shifted[0] - TMath::Cos(track.phi()) * qy_shifted[0]) * weight);
280+
}
281+
}
282+
}
283+
};
284+
285+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
286+
{
287+
return WorkflowSpec{
288+
adaptAnalysisTask<jEPFlowAnalysis>(cfgc)};
289+
}

0 commit comments

Comments
 (0)