Skip to content

Commit 7b9b39d

Browse files
ynishida-styleynishida-style
committed
2 parents 32b3eff + ba02103 commit 7b9b39d

38 files changed

+2004
-1871
lines changed

‎PWGCF/Flow/Tasks/flowEventPlane.cxx‎

Lines changed: 175 additions & 40 deletions
Original file line numberDiff line numberDiff line change
@@ -80,10 +80,19 @@ struct FlowEventPlane {
8080
// Tracks
8181
Configurable<float> cTrackMinPt{"cTrackMinPt", 0.15, "p_{T} minimum"};
8282
Configurable<float> cTrackMaxPt{"cTrackMaxPt", 2.0, "p_{T} maximum"};
83+
Configurable<int> cNEtaBins{"cNEtaBins", 7, "# of eta bins"};
8384
Configurable<float> cTrackEtaCut{"cTrackEtaCut", 0.8, "Pseudorapidity cut"};
8485
Configurable<bool> cTrackGlobal{"cTrackGlobal", true, "Global Track"};
8586
Configurable<float> cTrackDcaXYCut{"cTrackDcaXYCut", 0.1, "DcaXY Cut"};
8687
Configurable<float> cTrackDcaZCut{"cTrackDcaZCut", 1., "DcaXY Cut"};
88+
Configurable<int> cNRapBins{"cNRapBins", 5, "# of y bins"};
89+
Configurable<int> cNInvMassBins{"cNInvMassBins", 500, "# of m bins"};
90+
91+
// Track PID
92+
Configurable<float> cTpcNSigmaKa{"cTpcNSigmaKa", 2., "Tpc nsigma Ka"};
93+
Configurable<float> cTpcRejCut{"cTpcRejCut", 2., "Tpc rejection cut"};
94+
Configurable<float> cTofNSigmaKa{"cTofNSigmaKa", 2., "Tof nsigma Ka"};
95+
Configurable<float> cTofRejCut{"cTofRejCut", 2., "Tof rejection cut"};
8796

8897
// Gain calibration
8998
Configurable<bool> cDoGainCalib{"cDoGainCalib", false, "Gain Calib Flag"};
@@ -131,6 +140,16 @@ struct FlowEventPlane {
131140
{"hYZNCVsCent", "hYZNCVsVx", "hYZNCVsVy", "hYZNCVsVz"}};
132141
std::map<CorrectionType, std::vector<std::vector<std::string>>> corrTypeHistNameMap = {{kFineCorr, vFineCorrHistNames}, {kCoarseCorr, vCoarseCorrHistNames}};
133142

143+
// Container for histograms
144+
struct CorrectionHistContainer {
145+
TH2F* hGainCalib;
146+
std::array<std::array<THnSparseF*, 1>, 4> vCoarseCorrHist;
147+
std::array<std::array<TProfile*, 4>, 4> vFineCorrHist;
148+
} CorrectionHistContainer;
149+
150+
// Run number
151+
int cRunNum = 0, lRunNum = 0;
152+
134153
void init(InitContext const&)
135154
{
136155
// Set CCDB url
@@ -166,10 +185,13 @@ struct FlowEventPlane {
166185
const AxisSpec axisV1{400, -4, 4, "v_{1}"};
167186

168187
const AxisSpec axisTrackPt{100, 0., 10., "p_{T} (GeV/#it{c})"};
169-
const AxisSpec axisTrackEta{16, -0.8, 0.8, "#eta"};
188+
const AxisSpec axisTrackEta{cNEtaBins, -0.8, 0.8, "#eta"};
189+
const AxisSpec axisTrackRap{cNRapBins, -0.5, 0.5, "y"};
190+
const AxisSpec axisInvMass{cNInvMassBins, 0.87, 1.12, "M_{KK} (GeV/#it{c}^{2}"};
170191

171192
const AxisSpec axisTrackDcaXY{60, -0.15, 0.15, "DCA_{XY}"};
172193
const AxisSpec axisTrackDcaZ{230, -1.15, 1.15, "DCA_{XY}"};
194+
const AxisSpec axisTrackdEdx{360, 20, 200, "#frac{dE}{dx}"};
173195

174196
// Create histograms
175197
histos.add("Event/hCent", "FT0C%", kTH1F, {axisCent});
@@ -216,13 +238,20 @@ struct FlowEventPlane {
216238
histos.add("Checks/hYaXc", "Y^{A}_{1}X^{C}_{1}", kTProfile, {axisCent});
217239
histos.add("TrackQA/hPtDcaXY", "DCA_{XY} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaXY});
218240
histos.add("TrackQA/hPtDcaZ", "DCA_{Z} vs p_{T}", kTH2F, {axisTrackPt, axisTrackDcaZ});
241+
histos.add("TrackQA/hTrackTPCdEdX", "hTrackTPCdEdX", kTH2F, {axisTrackPt, axisTrackdEdx});
242+
histos.add("TrackQA/hUSCentPtInvMass", "hUSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass});
243+
histos.add("TrackQA/hLSCentPtInvMass", "hLSCentPtInvMass", kTH3F, {axisCent, axisTrackPt, axisInvMass});
219244
histos.add("DF/hQaQc", "X^{A}_{1}X^{C}_{1} + Y^{A}_{1}Y^{C}_{1}", kTProfile, {axisCent});
220245
histos.add("DF/hAQu", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
221246
histos.add("DF/hCQu", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
222247
histos.add("DF/hAQuPos", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
223248
histos.add("DF/hCQuPos", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
224249
histos.add("DF/hAQuNeg", "u_{x}X^{A}_{1} + u_{y}Y^{A}_{1}", kTProfile2D, {axisCent, axisTrackEta});
225250
histos.add("DF/hCQuNeg", "u_{x}X^{C}_{1} + u_{y}Y^{C}_{1}", kTProfile2D, {axisCent, axisTrackEta});
251+
histos.add("DF/Reso/US/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass});
252+
histos.add("DF/Reso/US/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass});
253+
histos.add("DF/Reso/LS/hPhiQuA", "hPhiQuA", kTProfile3D, {axisCent, axisTrackRap, axisInvMass});
254+
histos.add("DF/Reso/LS/hPhiQuC", "hPhiQuC", kTProfile3D, {axisCent, axisTrackRap, axisInvMass});
226255
}
227256

228257
template <typename C>
@@ -300,45 +329,52 @@ struct FlowEventPlane {
300329

301330
void gainCalib(float const& vz, int const& runNumber, std::array<float, 4>& energy, const char* histName = "hZNAEnergy")
302331
{
303-
// Set CCDB path
304-
std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/GainCalib" + "/Run" + std::to_string(runNumber);
332+
if (cRunNum != lRunNum) {
333+
// Set CCDB path
334+
std::string ccdbPath = static_cast<std::string>(cCcdbPath) + "/GainCalib" + "/Run" + std::to_string(runNumber);
305335

306-
// Get object from CCDB
307-
auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
336+
// Get object from CCDB
337+
auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
338+
339+
// Store histogram in container
340+
CorrectionHistContainer.hGainCalib = reinterpret_cast<TH2F*>(ccdbObj->FindObject(histName));
341+
}
308342

309-
// Get gain calibration
310-
TH2F* h = reinterpret_cast<TH2F*>(ccdbObj->FindObject(histName));
311343
float v = 0.;
312344
for (int i = 0; i < static_cast<int>(energy.size()); ++i) {
313-
v = h->GetBinContent(h->FindBin(i + 0.5, vz + 0.00001));
345+
v = CorrectionHistContainer.hGainCalib->GetBinContent(CorrectionHistContainer.hGainCalib->FindBin(i + 0.5, vz + 0.00001));
314346
energy[i] *= v;
315347
}
316348
}
317349

318-
template <typename C, typename F>
319-
std::vector<float> getAvgCorrFactors(const C& ccdbObject, const F& corrType, const std::vector<float>& vCollParam)
350+
std::vector<float> getAvgCorrFactors(CorrectionType const& corrType, const std::vector<float>& vCollParam)
320351
{
321352
std::vector<float> vAvgOutput = {0., 0., 0., 0.};
322-
std::vector<std::vector<std::string>> vHistNames = corrTypeHistNameMap.at(corrType);
323353
int binarray[4];
324-
int cntrx = 0;
325-
for (auto const& x : vHistNames) {
326-
int cntry = 0;
327-
for (auto const& y : x) {
328-
if (corrType == kFineCorr) {
329-
TProfile* hp = reinterpret_cast<TProfile*>(ccdbObject->FindObject(y.c_str()));
330-
vAvgOutput[cntrx] += hp->GetBinContent(hp->GetXaxis()->FindBin(vCollParam[cntry]));
331-
} else {
332-
THnSparseF* hn = reinterpret_cast<THnSparseF*>(ccdbObject->FindObject(y.c_str()));
333-
for (int i = 0; i < static_cast<int>(vHistNames.size()); ++i) {
334-
binarray[i] = hn->GetAxis(i)->FindBin(vCollParam[i]);
335-
}
336-
vAvgOutput[cntrx] += hn->GetBinContent(hn->GetBin(binarray));
354+
if (corrType == kCoarseCorr) {
355+
int cntrx = 0;
356+
for (auto const& v : CorrectionHistContainer.vCoarseCorrHist) {
357+
for (auto const& h : v) {
358+
binarray[kCent] = h->GetAxis(kCent)->FindBin(vCollParam[kCent] + 0.0001);
359+
binarray[kVx] = h->GetAxis(kVx)->FindBin(vCollParam[kVx] + 0.0001);
360+
binarray[kVy] = h->GetAxis(kVy)->FindBin(vCollParam[kVy] + 0.0001);
361+
binarray[kVz] = h->GetAxis(kVz)->FindBin(vCollParam[kVz] + 0.0001);
362+
vAvgOutput[cntrx] += h->GetBinContent(h->GetBin(binarray));
363+
}
364+
++cntrx;
365+
}
366+
} else {
367+
int cntrx = 0;
368+
for (auto const& v : CorrectionHistContainer.vFineCorrHist) {
369+
int cntry = 0;
370+
for (auto const& h : v) {
371+
vAvgOutput[cntrx] += h->GetBinContent(h->GetXaxis()->FindBin(vCollParam[cntry] + 0.0001));
372+
++cntry;
337373
}
338-
++cntry;
374+
++cntrx;
339375
}
340-
++cntrx;
341376
}
377+
342378
return vAvgOutput;
343379
}
344380

@@ -362,20 +398,39 @@ struct FlowEventPlane {
362398
corrType = kFineCorr;
363399
}
364400

365-
// Set ccdb path
366-
ccdbPath = static_cast<std::string>(cCcdbPath) + "/CorrItr_" + std::to_string(i + 1) + "/Run" + std::to_string(runNumber);
401+
// Check current and last run number, fetch ccdb object and store corrections in container
402+
if (cRunNum != lRunNum) {
403+
// Set ccdb path
404+
ccdbPath = static_cast<std::string>(cCcdbPath) + "/CorrItr_" + std::to_string(i + 1) + "/Run" + std::to_string(runNumber);
367405

368-
// Get object from CCDB
369-
auto ccdbObj = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
406+
// Get object from CCDB
407+
auto ccdbObject = ccdbService->getForTimeStamp<TList>(ccdbPath, -1);
370408

371-
// Check CCDB Object
372-
if (!ccdbObj) {
373-
LOGF(warning, "CCDB OBJECT NOT FOUND");
374-
return;
409+
// Check CCDB Object
410+
if (!ccdbObject) {
411+
LOGF(warning, "CCDB OBJECT NOT FOUND");
412+
return;
413+
}
414+
415+
// Store histograms in Hist Container
416+
std::vector<std::vector<std::string>> vHistNames = corrTypeHistNameMap.at(corrType);
417+
int cntrx = 0;
418+
for (auto const& x : vHistNames) {
419+
int cntry = 0;
420+
for (auto const& y : x) {
421+
if (corrType == kFineCorr) {
422+
CorrectionHistContainer.vFineCorrHist[cntrx][cntry] = reinterpret_cast<TProfile*>(ccdbObject->FindObject(y.c_str()));
423+
} else {
424+
CorrectionHistContainer.vCoarseCorrHist[cntrx][cntry] = reinterpret_cast<THnSparseF*>(ccdbObject->FindObject(y.c_str()));
425+
}
426+
++cntry;
427+
}
428+
++cntrx;
429+
}
375430
}
376431

377432
// Get averages
378-
std::vector<float> vAvg = getAvgCorrFactors(ccdbObj, corrType, inputParam);
433+
std::vector<float> vAvg = getAvgCorrFactors(corrType, inputParam);
379434

380435
// Apply correction
381436
outputParam[kXa] -= vAvg[kXa];
@@ -385,6 +440,80 @@ struct FlowEventPlane {
385440
}
386441
}
387442

443+
template <typename T>
444+
bool isKaon(T const& track)
445+
{
446+
bool tofFlagKa{false}, tpcFlagKa{false};
447+
// check tof signal
448+
if (track.hasTOF()) {
449+
if (std::abs(track.tofNSigmaKa()) < cTofNSigmaKa && std::abs(track.tofNSigmaPi()) > cTofRejCut && std::abs(track.tofNSigmaPr()) > cTofRejCut) {
450+
tofFlagKa = true;
451+
}
452+
if (std::abs(track.tpcNSigmaKa()) < cTpcNSigmaKa) {
453+
tpcFlagKa = true;
454+
}
455+
} else { // select from TPC
456+
tofFlagKa = true;
457+
if (std::abs(track.tpcNSigmaKa()) < cTpcNSigmaKa && std::abs(track.tpcNSigmaPi()) > cTpcRejCut && std::abs(track.tpcNSigmaPr()) > cTpcRejCut) {
458+
tpcFlagKa = true;
459+
}
460+
}
461+
462+
if (tofFlagKa && tpcFlagKa) {
463+
return true;
464+
}
465+
466+
return false;
467+
}
468+
469+
template <typename T>
470+
void getResoFlow(T const& tracks, std::vector<float> const& vSP)
471+
{
472+
float ux = 0., uy = 0., v1a = 0., v1c = 0.;
473+
for (auto const& [track1, track2] : soa::combinations(soa::CombinationsFullIndexPolicy(tracks, tracks))) {
474+
// Discard same track
475+
if (track1.index() == track2.index()) {
476+
continue;
477+
}
478+
479+
// Select track
480+
if (!selectTrack(track1) || !selectTrack(track2)) {
481+
continue;
482+
}
483+
484+
// Identify track
485+
if (!isKaon(track1) || !isKaon(track2)) {
486+
continue;
487+
}
488+
489+
histos.fill(HIST("TrackQA/hTrackTPCdEdX"), track1.pt(), track1.tpcSignal());
490+
histos.fill(HIST("TrackQA/hTrackTPCdEdX"), track2.pt(), track2.tpcSignal());
491+
492+
// Reconstruct invariant mass
493+
float p = RecoDecay::p((track1.px() + track2.px()), (track1.py() + track2.py()), (track1.pz() + track2.pz()));
494+
float e = RecoDecay::e(track1.px(), track1.py(), track1.pz(), MassKaonCharged) + RecoDecay::e(track2.px(), track2.py(), track2.pz(), MassKaonCharged);
495+
float m = std::sqrt(RecoDecay::m2(p, e));
496+
497+
// Get directed flow
498+
std::array<float, 3> v = {track1.px() + track2.px(), track1.py() + track2.py(), track1.pz() + track2.pz()};
499+
ux = std::cos(RecoDecay::phi(v));
500+
uy = std::sin(RecoDecay::phi(v));
501+
v1a = ux * vSP[kXa] + uy * vSP[kYa];
502+
v1c = ux * vSP[kXc] + uy * vSP[kYc];
503+
504+
// Fill histograms
505+
if (track1.sign() != track2.sign()) {
506+
histos.fill(HIST("TrackQA/hUSCentPtInvMass"), cent, RecoDecay::pt(v), m);
507+
histos.fill(HIST("DF/Reso/US/hPhiQuA"), cent, RecoDecay::y(v, MassPhi), m, v1a);
508+
histos.fill(HIST("DF/Reso/US/hPhiQuC"), cent, RecoDecay::y(v, MassPhi), m, v1c);
509+
} else {
510+
histos.fill(HIST("TrackQA/hLSCentPtInvMass"), cent, RecoDecay::pt(v), m);
511+
histos.fill(HIST("DF/Reso/LS/hPhiQuA"), cent, RecoDecay::y(v, MassPhi), m, v1a);
512+
histos.fill(HIST("DF/Reso/LS/hPhiQuC"), cent, RecoDecay::y(v, MassPhi), m, v1c);
513+
}
514+
}
515+
}
516+
388517
void fillCorrHist(const std::vector<float>& vCollParam, const std::vector<float>& vSP)
389518
{
390519
histos.fill(HIST("CorrHist/hWtXZNA"), vCollParam[kCent], vCollParam[kVx], vCollParam[kVy], vCollParam[kVz], vSP[kXa]);
@@ -422,7 +551,7 @@ struct FlowEventPlane {
422551

423552
using BCsRun3 = soa::Join<aod::BCsWithTimestamps, aod::Run3MatchedToBCSparse>;
424553
using CollisionsRun3 = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0Ms, aod::CentFV0As, aod::MultsExtra>;
425-
using Tracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCPi, aod::pidTPCPr, aod::pidTOFPi, aod::pidTOFPr, aod::TrackCompColls>;
554+
using Tracks = soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA, aod::pidTPCPi, aod::pidTOFPi, aod::pidTPCKa, aod::pidTOFKa, aod::pidTPCPr, aod::pidTOFPr, aod::TrackCompColls>;
426555

427556
void process(CollisionsRun3::iterator const& collision, BCsRun3 const& /* bcs*/, aod::Zdcs const&, Tracks const& tracks)
428557
{
@@ -442,7 +571,7 @@ struct FlowEventPlane {
442571

443572
// Get bunch crossing
444573
auto bc = collision.foundBC_as<BCsRun3>();
445-
int runNumber = collision.foundBC_as<BCsRun3>().runNumber();
574+
cRunNum = collision.foundBC_as<BCsRun3>().runNumber();
446575

447576
// check zdc
448577
if (!bc.has_zdc()) {
@@ -470,8 +599,8 @@ struct FlowEventPlane {
470599

471600
// Do gain calibration
472601
if (cDoGainCalib) {
473-
gainCalib(vCollParam[kVz], runNumber, znaEnergy, "hZNASignal");
474-
gainCalib(vCollParam[kVz], runNumber, zncEnergy, "hZNCSignal");
602+
gainCalib(vCollParam[kVz], cRunNum, znaEnergy, "hZNASignal");
603+
gainCalib(vCollParam[kVz], cRunNum, zncEnergy, "hZNCSignal");
475604
}
476605

477606
// Fill zdc signal
@@ -519,7 +648,7 @@ struct FlowEventPlane {
519648

520649
// Do corrections
521650
if (cApplyRecentCorr) {
522-
applyCorrection(vCollParam, runNumber, vSP);
651+
applyCorrection(vCollParam, cRunNum, vSP);
523652
}
524653

525654
// Fill X and Y histograms for corrections after each iteration
@@ -569,6 +698,12 @@ struct FlowEventPlane {
569698
histos.fill(HIST("DF/hCQuNeg"), cent, track.eta(), v1c);
570699
}
571700
}
701+
702+
// Get resonance flow
703+
getResoFlow(tracks, vSP);
704+
705+
// Update run number
706+
lRunNum = cRunNum;
572707
}
573708
};
574709

0 commit comments

Comments
 (0)