Skip to content

Commit 58d02bd

Browse files
Luzhiyonggalibuild
Luzhiyongg
and
alibuild
authored
[PWGCF] add long range cumulant (#15516)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 8447267 commit 58d02bd

File tree

1 file changed

+233
-0
lines changed

1 file changed

+233
-0
lines changed

‎PWGCF/TwoParticleCorrelations/Tasks/longRangeDihadronCor.cxx‎

Lines changed: 233 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,7 @@
1717
#include "PWGCF/Core/CorrelationContainer.h"
1818
#include "PWGCF/Core/PairCuts.h"
1919
#include "PWGCF/DataModel/CorrelationsDerived.h"
20+
#include "PWGCF/GenericFramework/Core/FlowContainer.h"
2021
#include "PWGCF/GenericFramework/Core/GFW.h"
2122
#include "PWGCF/GenericFramework/Core/GFWCumulant.h"
2223
#include "PWGCF/GenericFramework/Core/GFWPowerArray.h"
@@ -156,6 +157,19 @@ struct LongRangeDihadronCor {
156157
O2_DEFINE_CONFIGURABLE(cfgMirrorFT0CDeadChannels, bool, false, "If true, mirror FT0C channels 177->145, 176->144, 178->146, 179->147, 139->115")
157158
O2_DEFINE_CONFIGURABLE(cfgRunbyRunAmplitudeFT0, bool, false, "Produce run-by-run FT0 amplitude distribution");
158159
} cfgFwdConfig;
160+
struct : ConfigurableGroup {
161+
O2_DEFINE_CONFIGURABLE(gfwOutput, bool, false, "produce cumulant results, off by default");
162+
O2_DEFINE_CONFIGURABLE(gfwNbootstrap, int, 30, "Number of subsamples")
163+
O2_DEFINE_CONFIGURABLE(gfwAcceptance, std::string, "", "CCDB path to acceptance object")
164+
O2_DEFINE_CONFIGURABLE(gfwUseNch, bool, false, "use multiplicity as x axis");
165+
ConfigurableAxis gfwAxisIndependent{"gfwAxisIndependent", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
166+
Configurable<std::vector<std::string>> gfwCorr{"gfwCorr", std::vector<std::string>{"TPC {2} FT0A {-2}"}, "User defined GFW CorrelatorConfig"};
167+
Configurable<std::vector<std::string>> gfwName{"gfwName", std::vector<std::string>{"TpcFt0A22"}, "User defined GFW Name"};
168+
GFW* fGFW = new GFW();
169+
std::vector<GFW::CorrConfig> corrconfigs;
170+
TAxis* fPtAxis;
171+
TRandom3* fRndm = new TRandom3(0);
172+
} cfgCumulantConfig;
159173

160174
SliceCache cache;
161175

@@ -198,6 +212,7 @@ struct LongRangeDihadronCor {
198212
// Corrections
199213
TH3D* mEfficiency = nullptr;
200214
TH1D* mCentralityWeight = nullptr;
215+
GFWWeights* mAcceptance = nullptr;
201216
bool correctionsLoaded = false;
202217

203218
// Define the outputs
@@ -208,6 +223,8 @@ struct LongRangeDihadronCor {
208223
OutputObj<CorrelationContainer> sameFt0aFt0c{"sameEvent_FT0A_FT0C"};
209224
OutputObj<CorrelationContainer> mixedFt0aFt0c{"mixedEvent_FT0A_FT0C"};
210225
HistogramRegistry registry{"registry"};
226+
OutputObj<FlowContainer> fFC{FlowContainer("FlowContainer")};
227+
OutputObj<FlowContainer> fFCgen{FlowContainer("FlowContainer_gen")};
211228

212229
// define global variables
213230
TRandom3* gRandom = new TRandom3();
@@ -266,6 +283,10 @@ struct LongRangeDihadronCor {
266283
kFT0CMirrorChannelEnd = 147,
267284
kFT0CMirrorChannelInnerRing = 115
268285
};
286+
enum DataType {
287+
kReco,
288+
kGen
289+
};
269290
std::array<float, 6> tofNsigmaCut;
270291
std::array<float, 6> itsNsigmaCut;
271292
std::array<float, 6> tpcNsigmaCut;
@@ -442,6 +463,69 @@ struct LongRangeDihadronCor {
442463
mixedFt0aFt0c.setObject(new CorrelationContainer("mixedEvent_FT0A_FT0C", "mixedEvent_FT0A_FT0C", corrAxisFt0aFt0c, effAxis, userAxis));
443464
}
444465

466+
if (cfgCumulantConfig.gfwOutput && doprocessSameTpcFt0a && doprocessSameTpcFt0c) {
467+
LOGF(fatal, "If you enable gfwOutput, you can only enable processSameTpcFt0a or processSameTpcFt0c, NOT both.");
468+
}
469+
if (cfgCumulantConfig.gfwOutput) {
470+
o2::framework::AxisSpec axis = axisPt;
471+
int nPtBins = axis.binEdges.size() - 1;
472+
double* ptBins = &(axis.binEdges)[0];
473+
cfgCumulantConfig.fPtAxis = new TAxis(nPtBins, ptBins);
474+
475+
std::vector<std::string> userDefineGFWCorr = cfgCumulantConfig.gfwCorr;
476+
std::vector<std::string> userDefineGFWName = cfgCumulantConfig.gfwName;
477+
TObjArray* oba = new TObjArray();
478+
if (userDefineGFWName.size() != userDefineGFWCorr.size()) {
479+
LOGF(fatal, "The GFWConfig names you provided are NOT matching with configurations. userDefineGFWName.size(): %d, userDefineGFWCorr.size(): %d", userDefineGFWName.size(), userDefineGFWCorr.size());
480+
}
481+
LOGF(info, "User adding FlowContainer Array:");
482+
if (!userDefineGFWCorr.empty() && !userDefineGFWName.empty()) {
483+
for (uint i = 0; i < userDefineGFWName.size(); i++) {
484+
if (userDefineGFWCorr.at(i).find("poi") != std::string::npos) {
485+
LOGF(info, "%d: pT-diff array %s", i, userDefineGFWName.at(i).c_str());
486+
for (auto iPt = 0; iPt < cfgCumulantConfig.fPtAxis->GetNbins(); iPt++)
487+
oba->Add(new TNamed(Form("%s_pt_%i", userDefineGFWName.at(i).c_str(), iPt + 1), Form("%s_pTDiff", userDefineGFWName.at(i).c_str())));
488+
} else {
489+
LOGF(info, "%d: %s", i, userDefineGFWName.at(i).c_str());
490+
oba->Add(new TNamed(userDefineGFWName.at(i).c_str(), userDefineGFWName.at(i).c_str()));
491+
}
492+
}
493+
}
494+
fFC->SetName("FlowContainer");
495+
fFC->SetXAxis(cfgCumulantConfig.fPtAxis);
496+
fFC->Initialize(oba, cfgCumulantConfig.gfwAxisIndependent, cfgCumulantConfig.gfwNbootstrap);
497+
fFCgen->SetName("FlowContainer_gen");
498+
fFCgen->SetXAxis(cfgCumulantConfig.fPtAxis);
499+
fFCgen->Initialize(oba, cfgCumulantConfig.gfwAxisIndependent, cfgCumulantConfig.gfwNbootstrap);
500+
501+
cfgCumulantConfig.fGFW->AddRegion("TPC", -0.8, 0.8, 1, 1);
502+
cfgCumulantConfig.fGFW->AddRegion("TPCpoi", -0.8, 0.8, 1 + cfgCumulantConfig.fPtAxis->GetNbins(), 2);
503+
cfgCumulantConfig.fGFW->AddRegion("TPCol", -0.8, 0.8, 1 + cfgCumulantConfig.fPtAxis->GetNbins(), 4);
504+
cfgCumulantConfig.fGFW->AddRegion("FT0A", 3.5, 4.9, 1, 1);
505+
cfgCumulantConfig.fGFW->AddRegion("FT0C", -3.3, -2.1, 1, 1);
506+
cfgCumulantConfig.fGFW->AddRegion("FV0", 2.2, 5.1, 1, 1);
507+
508+
if (!userDefineGFWCorr.empty() && !userDefineGFWName.empty()) {
509+
LOGF(info, "User adding GFW CorrelatorConfig:");
510+
// attentaion: here we follow the index of cfgUserDefineGFWCorr
511+
for (uint i = 0; i < userDefineGFWCorr.size(); i++) {
512+
if (i >= userDefineGFWName.size()) {
513+
LOGF(fatal, "The names you provided are more than configurations. userDefineGFWName.size(): %d > userDefineGFWCorr.size(): %d", userDefineGFWName.size(), userDefineGFWCorr.size());
514+
break;
515+
}
516+
if (userDefineGFWCorr.at(i).find("poi") != std::string::npos) {
517+
cfgCumulantConfig.corrconfigs.push_back(cfgCumulantConfig.fGFW->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kTRUE));
518+
LOGF(info, "corrconfigs.at(%d): enable pt-Diff for %s %s", cfgCumulantConfig.corrconfigs.size() - 1, userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str());
519+
} else {
520+
cfgCumulantConfig.corrconfigs.push_back(cfgCumulantConfig.fGFW->GetCorrelatorConfig(userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str(), kFALSE));
521+
LOGF(info, "corrconfigs.at(%d): %s %s", cfgCumulantConfig.corrconfigs.size() - 1, userDefineGFWCorr.at(i).c_str(), userDefineGFWName.at(i).c_str());
522+
}
523+
}
524+
}
525+
526+
cfgCumulantConfig.fGFW->CreateRegions();
527+
}
528+
445529
LOGF(info, "End of init");
446530
}
447531

@@ -627,6 +711,13 @@ struct LongRangeDihadronCor {
627711
}
628712
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight);
629713
}
714+
if (cfgCumulantConfig.gfwOutput && cfgCumulantConfig.gfwAcceptance.value.empty() == false) {
715+
mAcceptance = ccdb->getForTimeStamp<GFWWeights>(cfgCumulantConfig.gfwAcceptance, timestamp);
716+
if (mAcceptance)
717+
LOGF(info, "Loaded acceptance weights from %s (%p)", cfgCumulantConfig.gfwAcceptance.value.c_str(), (void*)mAcceptance);
718+
else
719+
LOGF(warning, "Could not load acceptance weights from %s (%p)", cfgCumulantConfig.gfwAcceptance.value.c_str(), (void*)mAcceptance);
720+
}
630721
correctionsLoaded = true;
631722
}
632723

@@ -660,6 +751,80 @@ struct LongRangeDihadronCor {
660751
return true;
661752
}
662753

754+
bool getAcceptanceWeight(float& weight_nua, float phi, float eta, float vtxz)
755+
{
756+
if (mAcceptance)
757+
weight_nua = mAcceptance->getNUA(phi, eta, vtxz);
758+
else
759+
weight_nua = 1;
760+
return true;
761+
}
762+
763+
template <char... chars>
764+
void fillProfile(const GFW::CorrConfig& corrconf, const ConstStr<chars...>& tarName, const double& cent)
765+
{
766+
double dnx, val;
767+
dnx = cfgCumulantConfig.fGFW->Calculate(corrconf, 0, kTRUE).real();
768+
if (dnx == 0)
769+
return;
770+
if (!corrconf.pTDif) {
771+
val = cfgCumulantConfig.fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
772+
if (std::fabs(val) < 1)
773+
registry.fill(tarName, cent, val, dnx);
774+
return;
775+
}
776+
return;
777+
}
778+
779+
template <DataType dt>
780+
void fillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
781+
{
782+
double dnx, val;
783+
dnx = cfgCumulantConfig.fGFW->Calculate(corrconf, 0, kTRUE).real();
784+
if (!corrconf.pTDif) {
785+
if (dnx == 0)
786+
return;
787+
val = cfgCumulantConfig.fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
788+
if (std::fabs(val) < 1) {
789+
(dt == kGen) ? fFCgen->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm) : fFC->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
790+
}
791+
return;
792+
}
793+
for (auto i = 1; i <= cfgCumulantConfig.fPtAxis->GetNbins(); i++) {
794+
dnx = cfgCumulantConfig.fGFW->Calculate(corrconf, i - 1, kTRUE).real();
795+
if (dnx == 0)
796+
continue;
797+
val = cfgCumulantConfig.fGFW->Calculate(corrconf, i - 1, kFALSE).real() / dnx;
798+
if (std::fabs(val) < 1) {
799+
(dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm) : fFC->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
800+
}
801+
}
802+
return;
803+
}
804+
805+
template <typename TFT0s>
806+
void fillGFWFT0(TFT0s const& ft0, int corType)
807+
{
808+
std::size_t channelSize = 0;
809+
if (corType == kFT0C) {
810+
channelSize = ft0.channelC().size();
811+
} else if (corType == kFT0A) {
812+
channelSize = ft0.channelA().size();
813+
} else {
814+
LOGF(fatal, "Cor Index %d out of range", corType);
815+
}
816+
for (std::size_t iCh = 0; iCh < channelSize; iCh++) {
817+
int chanelid = 0;
818+
float ampl = 0.;
819+
getChannel(ft0, iCh, chanelid, ampl, corType, MixedEvent + 1);
820+
auto phi = getPhiFT0(chanelid, corType);
821+
auto eta = getEtaFT0(chanelid, corType);
822+
for (float ihit = 0; ihit < ampl; ihit++) {
823+
cfgCumulantConfig.fGFW->Fill(eta, 1, phi, 1., 1);
824+
}
825+
}
826+
}
827+
663828
// fill multiple histograms
664829
template <typename TCollision, typename TTracks>
665830
void fillYield(TCollision collision, TTracks tracks) // function to fill the yield and etaphi histograms.
@@ -1021,6 +1186,40 @@ struct LongRangeDihadronCor {
10211186
registry.fill(HIST("Nch"), tracks.size());
10221187
registry.fill(HIST("zVtx"), collision.posZ());
10231188

1189+
if (cfgCumulantConfig.gfwOutput) {
1190+
cfgCumulantConfig.fGFW->Clear();
1191+
float lRandom = cfgCumulantConfig.fRndm->Rndm();
1192+
float weff = 1, wacc = 1;
1193+
float independent = cent;
1194+
if (cfgCumulantConfig.gfwUseNch)
1195+
independent = static_cast<float>(tracks.size());
1196+
1197+
// fill TPC Q-vector
1198+
for (const auto& track : tracks) {
1199+
if (!trackSelected(track))
1200+
continue;
1201+
bool withinPtPOI = (0.2 < track.pt()) && (track.pt() < 10.0); // o2-linter: disable=magic-number (within POI pT range)
1202+
bool withinPtRef = (0.2 < track.pt()) && (track.pt() < 3.0); // o2-linter: disable=magic-number (within RF pT range)
1203+
getAcceptanceWeight(wacc, track.phi(), track.eta(), collision.posZ());
1204+
if (!getEfficiencyCorrection(weff, track.eta(), track.pt(), collision.posZ()))
1205+
continue;
1206+
if (withinPtRef)
1207+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 1);
1208+
if (withinPtPOI)
1209+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 2);
1210+
if (withinPtPOI && withinPtRef)
1211+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 4);
1212+
}
1213+
// fill FT0 Q-vector
1214+
const auto& ft0 = collision.foundFT0();
1215+
fillGFWFT0(ft0, kFT0A);
1216+
1217+
// Filling Flow Container
1218+
for (uint l_ind = 0; l_ind < cfgCumulantConfig.corrconfigs.size(); l_ind++) {
1219+
fillFC<kReco>(cfgCumulantConfig.corrconfigs.at(l_ind), independent, lRandom);
1220+
}
1221+
}
1222+
10241223
if (cfgSelCollByNch && (tracks.size() < cfgCutMultMin || tracks.size() >= cfgCutMultMax)) {
10251224
return;
10261225
}
@@ -1141,6 +1340,40 @@ struct LongRangeDihadronCor {
11411340
registry.fill(HIST("Nch"), tracks.size());
11421341
registry.fill(HIST("zVtx"), collision.posZ());
11431342

1343+
if (cfgCumulantConfig.gfwOutput) {
1344+
cfgCumulantConfig.fGFW->Clear();
1345+
float lRandom = cfgCumulantConfig.fRndm->Rndm();
1346+
float weff = 1, wacc = 1;
1347+
float independent = cent;
1348+
if (cfgCumulantConfig.gfwUseNch)
1349+
independent = static_cast<float>(tracks.size());
1350+
1351+
// fill TPC Q-vector
1352+
for (const auto& track : tracks) {
1353+
if (!trackSelected(track))
1354+
continue;
1355+
bool withinPtPOI = (0.2 < track.pt()) && (track.pt() < 10.0); // o2-linter: disable=magic-number (within POI pT range)
1356+
bool withinPtRef = (0.2 < track.pt()) && (track.pt() < 3.0); // o2-linter: disable=magic-number (within RF pT range)
1357+
getAcceptanceWeight(wacc, track.phi(), track.eta(), collision.posZ());
1358+
if (!getEfficiencyCorrection(weff, track.eta(), track.pt(), collision.posZ()))
1359+
continue;
1360+
if (withinPtRef)
1361+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 1);
1362+
if (withinPtPOI)
1363+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 2);
1364+
if (withinPtPOI && withinPtRef)
1365+
cfgCumulantConfig.fGFW->Fill(track.eta(), cfgCumulantConfig.fPtAxis->FindBin(track.pt()) - 1, track.phi(), wacc * weff, 4);
1366+
}
1367+
// fill FT0 Q-vector
1368+
const auto& ft0 = collision.foundFT0();
1369+
fillGFWFT0(ft0, kFT0C);
1370+
1371+
// Filling Flow Container
1372+
for (uint l_ind = 0; l_ind < cfgCumulantConfig.corrconfigs.size(); l_ind++) {
1373+
fillFC<kReco>(cfgCumulantConfig.corrconfigs.at(l_ind), independent, lRandom);
1374+
}
1375+
}
1376+
11441377
if (cfgSelCollByNch && (tracks.size() < cfgCutMultMin || tracks.size() >= cfgCutMultMax)) {
11451378
return;
11461379
}

0 commit comments

Comments
 (0)