Skip to content

Commit f276b4b

Browse files
CyF1rmianaCyF1rmiana
committed
Switch to cumulant flow calculation
1 parent e0461a8 commit f276b4b

1 file changed

Lines changed: 174 additions & 8 deletions

File tree

‎PWGCF/Flow/Tasks/flowEseTask.cxx‎

Lines changed: 174 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -9,12 +9,14 @@
99
// granted to it by virtue of its status as an Intergovernmental Organization
1010
// or submit itself to any jurisdiction.
1111

12-
/// \author Junlee Kim (jikim1290@gmail.com)
12+
/// \author Junlee Kim (jikim1290@gmail.com) & Shiqi Wang(shiqi.wang@cern.ch)
1313
/// \file flowEseTask.cxx
1414
/// \brief Task for flow and event shape engineering correlation with other observation.
1515
/// \since 2023-05-15
1616
/// \version 1.0
1717

18+
#include "PWGCF/GenericFramework/Core/GFW.h"
19+
#include "PWGCF/GenericFramework/Core/GFWWeights.h"
1820
#include "PWGLF/DataModel/LFStrangenessTables.h"
1921
#include "PWGMM/Mult/DataModel/Index.h" // for Particles2Tracks table
2022

@@ -46,6 +48,7 @@
4648
#include <Math/Vector4D.h> // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h)
4749
#include <Math/Vector4Dfwd.h>
4850
#include <TF1.h>
51+
#include <TH1.h>
4952
#include <TH2.h>
5053
#include <THnSparse.h>
5154
#include <TObject.h>
@@ -146,6 +149,13 @@ struct FlowEseTask {
146149
Configurable<bool> cfgAccCor{"cfgAccCor", false, "flag to apply acceptance correction"};
147150
Configurable<std::string> cfgAccCorPath{"cfgAccCorPath", "", "path for pseudo acceptance correction"};
148151

152+
struct : ConfigurableGroup {
153+
Configurable<bool> cfgGfwEffCor{"cfgGfwEffCor", false, "flag to apply charged-track efficiency correction for GFW"};
154+
Configurable<std::string> cfgGfwEffCorPath{"cfgGfwEffCorPath", "", "CCDB path to charged-track efficiency correction for GFW"};
155+
Configurable<bool> cfgGfwAccCor{"cfgGfwAccCor", false, "flag to apply charged-track acceptance correction for GFW"};
156+
Configurable<std::string> cfgGfwAccCorPath{"cfgGfwAccCorPath", "", "CCDB path to charged-track acceptance correction for GFW"};
157+
} cfgGfwParam;
158+
149159
Configurable<bool> cfgCalcCum{"cfgCalcCum", false, "flag to calculate cumulants of cossin"};
150160
Configurable<bool> cfgCalcCum1{"cfgCalcCum1", false, "flag to calculate cumulants of coscos"};
151161

@@ -188,6 +198,8 @@ struct FlowEseTask {
188198
static constexpr int NEseCentBins = 8;
189199
static constexpr int NEseGroups = 10;
190200
static constexpr int NEseCutsPerCentBin = NEseGroups - 1;
201+
static constexpr int GFWRefMask = 1;
202+
static constexpr int GFWPoiMask = 2;
191203
static constexpr float EseCentMin = 0.0f;
192204
static constexpr float EseCentMax = 80.0f;
193205
static constexpr float EseCentBinWidth = 10.0f;
@@ -220,6 +232,21 @@ struct FlowEseTask {
220232
std::vector<TProfile3D*> shiftprofile;
221233
TProfile2D* effMap = nullptr;
222234
TProfile2D* accMap = nullptr;
235+
struct {
236+
int lastGfwCorrectionRunNumber = -999;
237+
TH1D* gfwEffMap = nullptr;
238+
GFWWeights* gfwAccWeights = nullptr;
239+
TProfile2D* histRefQn = nullptr;
240+
TProfile3D* histPoiNQn = nullptr;
241+
TProfile3D* histPoiPQn = nullptr;
242+
GFW* fGFW = new GFW();
243+
GFW::CorrConfig corrRef22;
244+
GFW::CorrConfig corrPoiN22;
245+
GFW::CorrConfig corrPoiP22;
246+
GFW::CorrConfig corrRef32;
247+
GFW::CorrConfig corrPoiN32;
248+
GFW::CorrConfig corrPoiP32;
249+
} gfwState;
223250

224251
std::string fullCCDBShiftCorrPath;
225252

@@ -362,6 +389,15 @@ struct FlowEseTask {
362389
}
363390
histos.add("QA/CentDist", "", {HistType::kTH1F, {centQaAxis}});
364391
histos.add("QA/PVzDist", "", {HistType::kTH1F, {pVzQaAxis}});
392+
histos.add("psi2/h_Ref03Gap22", "", {HistType::kTProfile, {centAxis}});
393+
histos.add("psi2/h_PoiN03Gap22", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
394+
histos.add("psi2/h_PoiP03Gap22", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
395+
histos.add("psi3/h_Ref03Gap32", "", {HistType::kTProfile, {centAxis}});
396+
histos.add("psi3/h_PoiN03Gap32", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
397+
histos.add("psi3/h_PoiP03Gap32", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
398+
gfwState.histRefQn = histos.add<TProfile2D>(Form("psi%d/h_Ref03Gap%d2VsQ%d", cfgEseHarmonic.value, cfgEseHarmonic.value, cfgEseHarmonic.value), "", HistType::kTProfile2D, {centAxis, eseGroupAxis}).get();
399+
gfwState.histPoiNQn = histos.add<TProfile3D>(Form("psi%d/h_PoiN03Gap%d2VsQ%d", cfgEseHarmonic.value, cfgEseHarmonic.value, cfgEseHarmonic.value), "", HistType::kTProfile3D, {ptAxis, centAxis, eseGroupAxis}).get();
400+
gfwState.histPoiPQn = histos.add<TProfile3D>(Form("psi%d/h_PoiP03Gap%d2VsQ%d", cfgEseHarmonic.value, cfgEseHarmonic.value, cfgEseHarmonic.value), "", HistType::kTProfile3D, {ptAxis, centAxis, eseGroupAxis}).get();
365401

366402
for (auto i = 2; i < cfgnMods + 2; i++) {
367403
histos.add(Form("psi%d/h_lambda_cos", i), "", {HistType::kTHnSparseF, {massAxis, ptAxis, cosAxis, centAxis, epAxis}});
@@ -533,6 +569,21 @@ struct FlowEseTask {
533569
fMultPVCutHigh = new TF1("fMultPVCutHigh", "[0]+[1]*x+[2]*x*x+[3]*x*x*x + 2.5*([4]+[5]*x+[6]*x*x+[7]*x*x*x+[8]*x*x*x*x)", 0, 100);
534570
fMultPVCutHigh->SetParameters(2834.66, -87.0127, 0.915126, -0.00330136, 332.513, -12.3476, 0.251663, -0.00272819, 1.12242e-05);
535571

572+
const int nPtBins = histos.get<TProfile2D>(HIST("psi2/h_PoiN03Gap22"))->GetXaxis()->GetNbins();
573+
gfwState.fGFW->AddRegion("refN03", -0.8, -0.15, 1, GFWRefMask);
574+
gfwState.fGFW->AddRegion("refP03", 0.15, 0.8, 1, GFWRefMask);
575+
gfwState.fGFW->AddRegion("poiN03", -0.8, -0.15, nPtBins + 1, GFWPoiMask);
576+
gfwState.fGFW->AddRegion("poiP03", 0.15, 0.8, nPtBins + 1, GFWPoiMask);
577+
578+
gfwState.corrRef22 = gfwState.fGFW->GetCorrelatorConfig("refN03 {2} refP03 {-2}", "Ref03Gap22", false);
579+
gfwState.corrPoiN22 = gfwState.fGFW->GetCorrelatorConfig("poiN03 {2} refP03 {-2}", "PoiN03Gap22", true);
580+
gfwState.corrPoiP22 = gfwState.fGFW->GetCorrelatorConfig("poiP03 {2} refN03 {-2}", "PoiP03Gap22", true);
581+
gfwState.corrRef32 = gfwState.fGFW->GetCorrelatorConfig("refN03 {3} refP03 {-3}", "Ref03Gap32", false);
582+
gfwState.corrPoiN32 = gfwState.fGFW->GetCorrelatorConfig("poiN03 {3} refP03 {-3}", "PoiN03Gap32", true);
583+
gfwState.corrPoiP32 = gfwState.fGFW->GetCorrelatorConfig("poiP03 {3} refN03 {-3}", "PoiP03Gap32", true);
584+
585+
gfwState.fGFW->CreateRegions();
586+
536587
ccdb->setURL(cfgCcdbParam.cfgURL);
537588
ccdbApi.init("http://alice-ccdb.cern.ch");
538589
ccdb->setCaching(true);
@@ -878,24 +929,122 @@ struct FlowEseTask {
878929
}
879930

880931
template <typename TCollision, typename V0, typename TrackType>
881-
void fillHistograms(TCollision const& collision, V0 const& V0s, TrackType const& track, int nmode, int eseGroupIndex = -1, bool fillRegular = true)
932+
void fillHistograms(TCollision const& collision, V0 const& V0s, TrackType const& track, int nmode, int eseGroupIndex = -1, bool fillRegular = true, bool fillGfw = false)
882933
{
883934
qvecDetInd = detId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
884935
qvecRefAInd = refAId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
885936
qvecRefBInd = refBId * 4 + 3 + (nmode - 2) * cfgNQvec * 4;
886937
const bool fillEse = eseGroupIndex >= 0 && nmode == cfgEseHarmonic.value;
887938

939+
THnSparse* histEseVn = nullptr;
940+
double esePlane = 0.0;
888941
if (fillEse) {
889942
const auto suffix = eseGroupSuffix(eseGroupIndex);
890-
auto* histEseVn = getEseHistogram<THnSparse>(Form("histV%d_%s", cfgEseHarmonic.value, suffix.c_str()));
943+
histEseVn = getEseHistogram<THnSparse>(Form("histV%d_%s", cfgEseHarmonic.value, suffix.c_str()));
891944
const int harmonicIndex = nmode - 2;
892-
const double esePlane = helperEP.GetEventPlane(collision.qvecFT0CReVec()[harmonicIndex], collision.qvecFT0CImVec()[harmonicIndex], nmode);
945+
esePlane = helperEP.GetEventPlane(collision.qvecFT0CReVec()[harmonicIndex], collision.qvecFT0CImVec()[harmonicIndex], nmode);
946+
}
947+
948+
TProfile2D* histPoiN22 = nullptr;
949+
if (fillGfw) {
950+
gfwState.fGFW->Clear();
951+
histPoiN22 = histos.get<TProfile2D>(HIST("psi2/h_PoiN03Gap22")).get();
952+
}
953+
if (fillEse || fillGfw) {
893954
for (const auto& trk : track) {
894955
if (!selectionTrack(trk)) {
895956
continue;
896957
}
897-
const std::array<double, 3> values = {centrality, trk.pt(), std::cos(static_cast<float>(nmode) * (trk.phi() - esePlane))};
898-
histEseVn->Fill(values.data());
958+
if (fillEse) {
959+
const std::array<double, 3> values = {centrality, trk.pt(), std::cos(static_cast<float>(nmode) * (trk.phi() - esePlane))};
960+
histEseVn->Fill(values.data());
961+
}
962+
if (!fillGfw) {
963+
continue;
964+
}
965+
double weff = 1.0;
966+
double wacc = 1.0;
967+
if (cfgGfwParam.cfgGfwEffCor && gfwState.gfwEffMap != nullptr) {
968+
const double efficiency = gfwState.gfwEffMap->GetBinContent(gfwState.gfwEffMap->FindBin(trk.pt()));
969+
if (efficiency <= 0.0) {
970+
continue;
971+
}
972+
weff = 1.0 / efficiency;
973+
}
974+
if (cfgGfwParam.cfgGfwAccCor && gfwState.gfwAccWeights != nullptr) {
975+
wacc = gfwState.gfwAccWeights->getNUA(trk.phi(), trk.eta(), collision.posZ());
976+
}
977+
const double trackWeight = weff * wacc;
978+
const int ptBin = histPoiN22->GetXaxis()->FindBin(trk.pt()) - 1;
979+
gfwState.fGFW->Fill(trk.eta(), ptBin, trk.phi(), trackWeight, GFWRefMask);
980+
gfwState.fGFW->Fill(trk.eta(), ptBin, trk.phi(), trackWeight, GFWPoiMask);
981+
}
982+
}
983+
984+
if (fillGfw) {
985+
double denominator = gfwState.fGFW->Calculate(gfwState.corrRef22, 0, true).real();
986+
if (denominator > 0.0) {
987+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrRef22, 0, false).real() / denominator;
988+
if (std::abs(correlation) < 1.0) {
989+
histos.fill(HIST("psi2/h_Ref03Gap22"), centrality, correlation, denominator);
990+
if (fillEse && cfgEseHarmonic.value == SecondHarmonic) {
991+
gfwState.histRefQn->Fill(centrality, eseGroupIndex + 0.5, correlation, denominator);
992+
}
993+
}
994+
}
995+
denominator = gfwState.fGFW->Calculate(gfwState.corrRef32, 0, true).real();
996+
if (denominator > 0.0) {
997+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrRef32, 0, false).real() / denominator;
998+
if (std::abs(correlation) < 1.0) {
999+
histos.fill(HIST("psi3/h_Ref03Gap32"), centrality, correlation, denominator);
1000+
if (fillEse && cfgEseHarmonic.value == ThirdHarmonic) {
1001+
gfwState.histRefQn->Fill(centrality, eseGroupIndex + 0.5, correlation, denominator);
1002+
}
1003+
}
1004+
}
1005+
1006+
for (int iPt = 1; iPt <= histPoiN22->GetXaxis()->GetNbins(); ++iPt) {
1007+
const double pt = histPoiN22->GetXaxis()->GetBinCenter(iPt);
1008+
denominator = gfwState.fGFW->Calculate(gfwState.corrPoiN22, iPt - 1, true).real();
1009+
if (denominator > 0.0) {
1010+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrPoiN22, iPt - 1, false).real() / denominator;
1011+
if (std::abs(correlation) < 1.0) {
1012+
histos.fill(HIST("psi2/h_PoiN03Gap22"), pt, centrality, correlation, denominator);
1013+
if (fillEse && cfgEseHarmonic.value == SecondHarmonic) {
1014+
gfwState.histPoiNQn->Fill(pt, centrality, eseGroupIndex + 0.5, correlation, denominator);
1015+
}
1016+
}
1017+
}
1018+
denominator = gfwState.fGFW->Calculate(gfwState.corrPoiP22, iPt - 1, true).real();
1019+
if (denominator > 0.0) {
1020+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrPoiP22, iPt - 1, false).real() / denominator;
1021+
if (std::abs(correlation) < 1.0) {
1022+
histos.fill(HIST("psi2/h_PoiP03Gap22"), pt, centrality, correlation, denominator);
1023+
if (fillEse && cfgEseHarmonic.value == SecondHarmonic) {
1024+
gfwState.histPoiPQn->Fill(pt, centrality, eseGroupIndex + 0.5, correlation, denominator);
1025+
}
1026+
}
1027+
}
1028+
denominator = gfwState.fGFW->Calculate(gfwState.corrPoiN32, iPt - 1, true).real();
1029+
if (denominator > 0.0) {
1030+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrPoiN32, iPt - 1, false).real() / denominator;
1031+
if (std::abs(correlation) < 1.0) {
1032+
histos.fill(HIST("psi3/h_PoiN03Gap32"), pt, centrality, correlation, denominator);
1033+
if (fillEse && cfgEseHarmonic.value == ThirdHarmonic) {
1034+
gfwState.histPoiNQn->Fill(pt, centrality, eseGroupIndex + 0.5, correlation, denominator);
1035+
}
1036+
}
1037+
}
1038+
denominator = gfwState.fGFW->Calculate(gfwState.corrPoiP32, iPt - 1, true).real();
1039+
if (denominator > 0.0) {
1040+
const double correlation = gfwState.fGFW->Calculate(gfwState.corrPoiP32, iPt - 1, false).real() / denominator;
1041+
if (std::abs(correlation) < 1.0) {
1042+
histos.fill(HIST("psi3/h_PoiP03Gap32"), pt, centrality, correlation, denominator);
1043+
if (fillEse && cfgEseHarmonic.value == ThirdHarmonic) {
1044+
gfwState.histPoiPQn->Fill(pt, centrality, eseGroupIndex + 0.5, correlation, denominator);
1045+
}
1046+
}
1047+
}
8991048
}
9001049
}
9011050

@@ -1261,12 +1410,28 @@ struct FlowEseTask {
12611410
if (cfgAccCor) {
12621411
accMap = ccdb->getForTimeStamp<TProfile2D>(cfgAccCorPath.value, bc.timestamp());
12631412
}
1413+
if (bc.runNumber() != gfwState.lastGfwCorrectionRunNumber) {
1414+
if (cfgGfwParam.cfgGfwEffCor && !cfgGfwParam.cfgGfwEffCorPath.value.empty()) {
1415+
gfwState.gfwEffMap = ccdb->getForTimeStamp<TH1D>(cfgGfwParam.cfgGfwEffCorPath.value, bc.timestamp());
1416+
if (gfwState.gfwEffMap == nullptr) {
1417+
LOGF(fatal, "Could not load charged-track efficiency histogram from %s", cfgGfwParam.cfgGfwEffCorPath.value.c_str());
1418+
}
1419+
}
1420+
if (cfgGfwParam.cfgGfwAccCor && !cfgGfwParam.cfgGfwAccCorPath.value.empty()) {
1421+
gfwState.gfwAccWeights = ccdb->getForTimeStamp<GFWWeights>(cfgGfwParam.cfgGfwAccCorPath.value, bc.timestamp());
1422+
if (gfwState.gfwAccWeights == nullptr) {
1423+
LOGF(fatal, "Could not load charged-track acceptance weights from %s", cfgGfwParam.cfgGfwAccCorPath.value.c_str());
1424+
}
1425+
}
1426+
gfwState.lastGfwCorrectionRunNumber = bc.runNumber();
1427+
}
1428+
12641429
fillEseEPQA(collision, eseGroupIndex);
12651430
if (cfgShiftCorrDef && cfgEseHarmonic.value >= cfgnMods.value + 2) {
12661431
fillShiftCorrection(collision, cfgEseHarmonic.value);
12671432
}
12681433
if (eseGroupIndex >= 0 && cfgEseHarmonic.value >= cfgnMods.value + 2) {
1269-
fillHistograms(collision, V0s, tracks, cfgEseHarmonic.value, eseGroupIndex, false);
1434+
fillHistograms(collision, V0s, tracks, cfgEseHarmonic.value, eseGroupIndex, false, true);
12701435
}
12711436
for (int i = 2; i < cfgnMods + 2; i++) {
12721437
if (cfgShiftCorrDef) {
@@ -1275,7 +1440,8 @@ struct FlowEseTask {
12751440
if (cfgQAv0) {
12761441
fillEPQA(collision, i);
12771442
}
1278-
fillHistograms(collision, V0s, tracks, i, i == cfgEseHarmonic.value ? eseGroupIndex : -1);
1443+
const bool fillGfw = eseGroupIndex >= 0 ? i == cfgEseHarmonic.value : i == SecondHarmonic;
1444+
fillHistograms(collision, V0s, tracks, i, i == cfgEseHarmonic.value ? eseGroupIndex : -1, true, fillGfw);
12791445
} // FIXME: need to fill different histograms for different harmonic
12801446
}
12811447
PROCESS_SWITCH(FlowEseTask, processData, "Process Event for data", true);

0 commit comments

Comments
 (0)