Skip to content

Commit 4a82301

Browse files
MaximVirtaMaxim Virtaalibuild
authored
[PWGCF] Added multiplicity study & pt correlations (#17261)
Co-authored-by: Maxim Virta <maximus@eduroam-174093473-3.dyndns.cern.ch> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent aa64e13 commit 4a82301

1 file changed

Lines changed: 103 additions & 9 deletions

File tree

‎PWGCF/GenericFramework/Tasks/flowGfwV02.cxx‎

Lines changed: 103 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -92,7 +92,7 @@ float philow = 0.0;
9292
float phiup = o2::constants::math::TwoPI;
9393
int nchbins = 300;
9494
float nchlow = 0;
95-
float nchup = 3000;
95+
float nchup = 300;
9696
std::vector<double> centbinning(90);
9797
int nBootstrap = 10;
9898
std::vector<std::pair<double, double>> etagapsPtPt;
@@ -104,6 +104,7 @@ std::vector<double> multGlobalPVCorrCutPars;
104104
} // namespace o2::analysis::gfw
105105

106106
struct FlowGfwV02 {
107+
107108
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
108109
O2_DEFINE_CONFIGURABLE(cfgMpar, int, 4, "Highest order of pt-pt correlations")
109110
O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FT0A")
@@ -124,6 +125,7 @@ struct FlowGfwV02 {
124125
O2_DEFINE_CONFIGURABLE(cfgNormalizeByCharged, bool, true, "Enable or disable the normalization by charged particles");
125126
O2_DEFINE_CONFIGURABLE(cfgConsistentEventFlag, int, 15, "Flag for consistent event selection");
126127
O2_DEFINE_CONFIGURABLE(cfgMultCut, bool, true, "Use additional event cut on mult correlations");
128+
O2_DEFINE_CONFIGURABLE(cfgUseV0, bool, false, "Use V0 analysis");
127129

128130
// Event selection cuts
129131
struct : ConfigurableGroup {
@@ -139,6 +141,16 @@ struct FlowGfwV02 {
139141
O2_DEFINE_CONFIGURABLE(cfgIsVertexITSTPC, bool, true, "kIsVertexITSTPC - Selects collisions with at least one ITS-TPC track");
140142
} cfgEventCutFlags;
141143

144+
// Event selection cuts
145+
struct : ConfigurableGroup {
146+
O2_DEFINE_CONFIGURABLE(cfgEtaSubAMin, float, -0.8, "Minimum eta for subevent A");
147+
O2_DEFINE_CONFIGURABLE(cfgEtaSubAMax, float, -0.5, "Maximum eta for subevent A");
148+
O2_DEFINE_CONFIGURABLE(cfgEtaSubBMin, float, 0.5, "Minimum eta for subevent B");
149+
O2_DEFINE_CONFIGURABLE(cfgEtaSubBMax, float, 0.8, "Maximum eta for subevent B");
150+
O2_DEFINE_CONFIGURABLE(cfgEtaSubCMin, float, -0.4, "Minimum eta for subevent C");
151+
O2_DEFINE_CONFIGURABLE(cfgEtaSubCMax, float, 0.4, "Maximum eta for subevent C");
152+
} cfgSubeventCuts;
153+
142154
struct : ConfigurableGroup {
143155
Configurable<std::vector<double>> cfgMultGlobalCutPars{"cfgMultGlobalCutPars", std::vector<double>{2272.16, -76.6932, 1.01204, -0.00631545, 1.59868e-05, 136.336, -4.97006, 0.121199, -0.0015921, 7.66197e-06}, "Global vs FT0C multiplicity cut parameter values"};
144156
Configurable<std::vector<double>> cfgMultPVCutPars{"cfgMultPVCutPars", std::vector<double>{3074.43, -106.192, 1.46176, -0.00968364, 2.61923e-05, 182.128, -7.43492, 0.193901, -0.00256715, 1.22594e-05}, "PV vs FT0C multiplicity cut parameter values"};
@@ -253,6 +265,8 @@ struct FlowGfwV02 {
253265
std::array<float, 6> itsNsigmaCut;
254266
std::array<float, 6> tpcNsigmaCut;
255267
std::array<std::unique_ptr<TH1D>, 4> hPtMid{};
268+
std::array<std::unique_ptr<TH1D>, 4> hPtForward{};
269+
std::array<std::unique_ptr<TH1D>, 4> hPtBackward{};
256270
};
257271
PIDState pidStates;
258272

@@ -275,7 +289,8 @@ struct FlowGfwV02 {
275289
PidCharged = 0,
276290
PidPions,
277291
PidKaons,
278-
PidProtons
292+
PidProtons,
293+
PidTotal
279294
};
280295
enum PiKpArrayIndex {
281296
IndPionUp = 0,
@@ -378,6 +393,24 @@ struct FlowGfwV02 {
378393
pidStates.hPtMid[PidKaons]->SetDirectory(nullptr);
379394
pidStates.hPtMid[PidProtons]->SetDirectory(nullptr);
380395

396+
pidStates.hPtForward[PidCharged] = std::make_unique<TH1D>("hPtForward_charged", "hPtForward_charged", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
397+
pidStates.hPtForward[PidPions] = std::make_unique<TH1D>("hPtForward_pions", "hPtForward_pions", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
398+
pidStates.hPtForward[PidKaons] = std::make_unique<TH1D>("hPtForward_kaons", "hPtForward_kaons", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
399+
pidStates.hPtForward[PidProtons] = std::make_unique<TH1D>("hPtForward_protons", "hPtForward_protons", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
400+
pidStates.hPtForward[PidCharged]->SetDirectory(nullptr);
401+
pidStates.hPtForward[PidPions]->SetDirectory(nullptr);
402+
pidStates.hPtForward[PidKaons]->SetDirectory(nullptr);
403+
pidStates.hPtForward[PidProtons]->SetDirectory(nullptr);
404+
405+
pidStates.hPtBackward[PidCharged] = std::make_unique<TH1D>("hPtBackward_charged", "hPtBackward_charged", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
406+
pidStates.hPtBackward[PidPions] = std::make_unique<TH1D>("hPtBackward_pions", "hPtBackward_pions", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
407+
pidStates.hPtBackward[PidKaons] = std::make_unique<TH1D>("hPtBackward_kaons", "hPtBackward_kaons", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
408+
pidStates.hPtBackward[PidProtons] = std::make_unique<TH1D>("hPtBackward_protons", "hPtBackward_protons", o2::analysis::gfw::ptbinning.size() - 1, &o2::analysis::gfw::ptbinning[0]);
409+
pidStates.hPtBackward[PidCharged]->SetDirectory(nullptr);
410+
pidStates.hPtBackward[PidPions]->SetDirectory(nullptr);
411+
pidStates.hPtBackward[PidKaons]->SetDirectory(nullptr);
412+
pidStates.hPtBackward[PidProtons]->SetDirectory(nullptr);
413+
381414
AxisSpec phiAxis = {o2::analysis::gfw::phibins, o2::analysis::gfw::philow, o2::analysis::gfw::phiup, "#phi"};
382415
AxisSpec etaAxis = {o2::analysis::gfw::etabins, -cfgTrackCuts.cfgEtaMax, cfgTrackCuts.cfgEtaMax, "#eta"};
383416
AxisSpec vtxAxis = {o2::analysis::gfw::vtxZbins, -cfgEventCuts.cfgZvtxMax, cfgEventCuts.cfgZvtxMax, "Vtx_{z} (cm)"};
@@ -398,9 +431,19 @@ struct FlowGfwV02 {
398431
AxisSpec multpvAxis = {600, 0, 600, "N_{ch} (PV)"};
399432
AxisSpec dcaZAxis = {200, -2, 2, "DCA_{z} (cm)"};
400433
AxisSpec dcaXYAxis = {200, -0.5, 0.5, "DCA_{xy} (cm)"};
434+
AxisSpec pidAxis = {4, -0.5, 3.5, "PID"}; // 0 = not identified, 1 = pion, 2 = kaon, 3 = proton
401435

402436
registry.add("v02pt", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
403-
registry.add("nchMid", "", {HistType::kTProfile2D, {ptAxis, centAxis}});
437+
registry.add("nchMid", "", {HistType::kTProfile3D, {ptAxis, centAxis, nchAxis}});
438+
registry.add("v02centmult", "", {HistType::kTProfile2D, {centAxis, nchAxis}});
439+
440+
registry.add("analysis/v0AB", "", {HistType::kTProfile3D, {pidAxis, ptAxis, centAxis}});
441+
registry.add("analysis/v0BA", "", {HistType::kTProfile3D, {pidAxis, ptAxis, centAxis}});
442+
registry.add("analysis/nchA", "", {HistType::kTProfile3D, {pidAxis, ptAxis, centAxis}});
443+
registry.add("analysis/nchB", "", {HistType::kTProfile3D, {pidAxis, ptAxis, centAxis}});
444+
registry.add("analysis/ptA", "", {HistType::kTProfile3D, {pidAxis, centAxis, nchAxis}});
445+
registry.add("analysis/ptB", "", {HistType::kTProfile3D, {pidAxis, centAxis, nchAxis}});
446+
registry.add("analysis/ptAB", "", {HistType::kTProfile3D, {pidAxis, centAxis, nchAxis}});
404447

405448
ccdb->setURL("http://alice-ccdb.cern.ch");
406449
ccdb->setCaching(true);
@@ -614,7 +657,7 @@ struct FlowGfwV02 {
614657
}
615658
if (cfgPIDEfficiency) {
616659
const std::array<std::string, 4> pidStrings = {"ch", "pi", "ka", "pr"};
617-
for (int i = 1; i < 4; i++) {
660+
for (int i = 1; i < PidTotal; i++) {
618661

619662
cfg.mEfficiency[i] = ccdb->getForTimeStamp<TH1D>(cfgEfficiency.value + pidStrings[i], timestamp);
620663
if (cfg.mEfficiency[i] == nullptr) {
@@ -806,7 +849,7 @@ struct FlowGfwV02 {
806849
}
807850

808851
template <DataType dt>
809-
void fillOutputContainers(const float& centmult, const double& rndm, const int& /*run*/ = 0)
852+
void fillOutputContainers(const float& centmult, const int& multiplicity, const double& rndm, const int& /*run*/ = 0)
810853
{
811854
for (uint l_ind = 0; l_ind < corrconfigs.size(); ++l_ind) {
812855
if (!corrconfigs.at(l_ind).pTDif) {
@@ -846,6 +889,32 @@ struct FlowGfwV02 {
846889
}
847890
}
848891
}
892+
893+
if (cfgUseV0) {
894+
double v0corrAB = 0;
895+
double v0corrBA = 0;
896+
double ptMeanForward = pidStates.hPtForward[PidCharged]->GetMean();
897+
double ptMeanBackward = pidStates.hPtBackward[PidCharged]->GetMean();
898+
double ptFractionForward = 0.;
899+
double ptFractionBackward = 0.;
900+
for (int pid = 0; pid < PidTotal; pid++) {
901+
int normIndex = (cfgNormalizeByCharged) ? PidCharged : pid;
902+
for (int i = 1; i <= fSecondAxis->GetNbins(); i++) {
903+
ptFractionForward = pidStates.hPtForward[pid]->GetBinContent(i) / pidStates.hPtForward[normIndex]->Integral();
904+
ptFractionBackward = pidStates.hPtBackward[pid]->GetBinContent(i) / pidStates.hPtBackward[normIndex]->Integral();
905+
v0corrAB = ptFractionForward * ptMeanBackward;
906+
v0corrBA = ptFractionBackward * ptMeanForward;
907+
registry.fill(HIST("analysis/v0AB"), pid, fSecondAxis->GetBinCenter(i), centmult, v0corrAB);
908+
registry.fill(HIST("analysis/v0BA"), pid, fSecondAxis->GetBinCenter(i), centmult, v0corrBA);
909+
registry.fill(HIST("analysis/nchA"), pid, fSecondAxis->GetBinCenter(i), centmult, ptFractionForward);
910+
registry.fill(HIST("analysis/nchB"), pid, fSecondAxis->GetBinCenter(i), centmult, ptFractionBackward);
911+
}
912+
registry.fill(HIST("analysis/ptA"), pid, centmult, multiplicity, ptMeanForward);
913+
registry.fill(HIST("analysis/ptB"), pid, centmult, multiplicity, ptMeanBackward);
914+
registry.fill(HIST("analysis/ptAB"), pid, centmult, multiplicity, ptMeanForward * ptMeanBackward);
915+
}
916+
}
917+
849918
// Fill the profiles for each pT bin
850919
auto dnx = fGFW->Calculate(corrconfigs.at(0), 0, kTRUE).real();
851920
if (dnx == 0)
@@ -857,9 +926,10 @@ struct FlowGfwV02 {
857926
ptFraction = pidStates.hPtMid[PidCharged]->GetBinContent(i) / pidStates.hPtMid[PidCharged]->Integral();
858927
if (std::abs(val) < 1)
859928
registry.fill(HIST("v02pt"), fSecondAxis->GetBinCenter(i), centmult, val * ptFraction, (cfgUseMultiplicityFlowWeights) ? dnx : 1.0);
860-
registry.fill(HIST("nchMid"), fSecondAxis->GetBinCenter(i), centmult, ptFraction);
929+
registry.fill(HIST("nchMid"), fSecondAxis->GetBinCenter(i), centmult, multiplicity, ptFraction);
861930
}
862931
}
932+
registry.fill(HIST("v02centmult"), centmult, multiplicity, val);
863933
return;
864934
}
865935

@@ -891,23 +961,47 @@ struct FlowGfwV02 {
891961
pidStates.hPtMid[PidPions]->Reset();
892962
pidStates.hPtMid[PidKaons]->Reset();
893963
pidStates.hPtMid[PidProtons]->Reset();
964+
pidStates.hPtBackward[PidCharged]->Reset();
965+
pidStates.hPtBackward[PidPions]->Reset();
966+
pidStates.hPtBackward[PidKaons]->Reset();
967+
pidStates.hPtBackward[PidProtons]->Reset();
968+
pidStates.hPtForward[PidCharged]->Reset();
969+
pidStates.hPtForward[PidPions]->Reset();
970+
pidStates.hPtForward[PidKaons]->Reset();
971+
pidStates.hPtForward[PidProtons]->Reset();
894972

895973
float lRandom = fRndm->Rndm();
896974

897975
// Loop over tracks and check if they are accepted
898976
AcceptedTracks acceptedTracks{0, 0, 0, 0};
899977
for (const auto& track : tracks) {
900978
processTrack(track, vtxz, xaxis.multiplicity, run, acceptedTracks);
901-
if (track.eta() > -0.4 && track.eta() < 0.4)
979+
if (track.eta() > cfgSubeventCuts.cfgEtaSubCMin && track.eta() < cfgSubeventCuts.cfgEtaSubCMax)
902980
pidStates.hPtMid[PidCharged]->Fill(track.pt(), getEfficiency(track, PidCharged));
981+
if (track.eta() > cfgSubeventCuts.cfgEtaSubAMin && track.eta() < cfgSubeventCuts.cfgEtaSubAMax) // add mean pT
982+
pidStates.hPtBackward[PidCharged]->Fill(track.pt(), getEfficiency(track, PidCharged));
983+
if (track.eta() > cfgSubeventCuts.cfgEtaSubBMin && track.eta() < cfgSubeventCuts.cfgEtaSubBMax) // add mean pT
984+
pidStates.hPtForward[PidCharged]->Fill(track.pt(), getEfficiency(track, PidCharged));
903985
// If PID is identified, fill pt spectrum for the corresponding particle
904986
int pidInd = getNsigmaPID(track);
905-
if (pidInd != -1 && track.eta() > -0.4 && track.eta() < 0.4) {
987+
if (pidInd != -1 && track.eta() > cfgSubeventCuts.cfgEtaSubCMin && track.eta() < cfgSubeventCuts.cfgEtaSubCMax) {
906988
if (cfgPIDEfficiency)
907989
pidStates.hPtMid[pidInd]->Fill(track.pt(), getEfficiency(track, pidInd));
908990
else
909991
pidStates.hPtMid[pidInd]->Fill(track.pt(), getEfficiency(track, PidCharged)); // Default to charged particles if PID efficiency is not used
910992
}
993+
if (pidInd != -1 && track.eta() > cfgSubeventCuts.cfgEtaSubAMin && track.eta() < cfgSubeventCuts.cfgEtaSubAMax) {
994+
if (cfgPIDEfficiency)
995+
pidStates.hPtBackward[pidInd]->Fill(track.pt(), getEfficiency(track, pidInd));
996+
else
997+
pidStates.hPtBackward[pidInd]->Fill(track.pt(), getEfficiency(track, PidCharged)); // Default to charged particles if PID efficiency is not used
998+
}
999+
if (pidInd != -1 && track.eta() > cfgSubeventCuts.cfgEtaSubBMin && track.eta() < cfgSubeventCuts.cfgEtaSubBMax) {
1000+
if (cfgPIDEfficiency)
1001+
pidStates.hPtForward[pidInd]->Fill(track.pt(), getEfficiency(track, pidInd));
1002+
else
1003+
pidStates.hPtForward[pidInd]->Fill(track.pt(), getEfficiency(track, PidCharged)); // Default to charged particles if PID efficiency is not used
1004+
}
9111005
}
9121006
if (cfgConsistentEventFlag & 1)
9131007
if (!acceptedTracks.nPos || !acceptedTracks.nNeg)
@@ -922,7 +1016,7 @@ struct FlowGfwV02 {
9221016
if (acceptedTracks.nPos < 2 || acceptedTracks.nMid < 2 || acceptedTracks.nNeg < 2) // o2-linter: disable=magic-number (at least two tracks in all three subevents)
9231017
return;
9241018
// Fill output containers
925-
fillOutputContainers<dt>(xaxis.centrality, lRandom, run);
1019+
fillOutputContainers<dt>(xaxis.centrality, xaxis.multiplicity, lRandom, run);
9261020
}
9271021

9281022
template <typename TTrack>

0 commit comments

Comments
 (0)