Skip to content

Commit 65261a2

Browse files
LuzhiyonggLuzhiyongg
authored
add long range cumulant
Added FlowContainer and GFW configurations for flow analysis.
1 parent 344bfe5 commit 65261a2

File tree

1 file changed

+234
-0
lines changed

1 file changed

+234
-0
lines changed

‎PWGCF/TwoParticleCorrelations/Tasks/longRangeDihadronCor.cxx‎

Lines changed: 234 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,70 @@ 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+
527+
cfgCumulantConfig.fGFW->CreateRegions();
528+
}
529+
445530
LOGF(info, "End of init");
446531
}
447532

@@ -627,6 +712,13 @@ struct LongRangeDihadronCor {
627712
}
628713
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgCentralityWeight.value.c_str(), (void*)mCentralityWeight);
629714
}
715+
if (cfgCumulantConfig.gfwOutput && cfgCumulantConfig.gfwAcceptance.value.empty() == false) {
716+
mAcceptance = ccdb->getForTimeStamp<GFWWeights>(cfgCumulantConfig.gfwAcceptance, timestamp);
717+
if (mAcceptance)
718+
LOGF(info, "Loaded acceptance weights from %s (%p)", cfgCumulantConfig.gfwAcceptance.value.c_str(), (void*)mAcceptance);
719+
else
720+
LOGF(warning, "Could not load acceptance weights from %s (%p)", cfgCumulantConfig.gfwAcceptance.value.c_str(), (void*)mAcceptance);
721+
}
630722
correctionsLoaded = true;
631723
}
632724

@@ -660,6 +752,80 @@ struct LongRangeDihadronCor {
660752
return true;
661753
}
662754

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

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

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

0 commit comments

Comments
 (0)