Skip to content

Commit efba710

Browse files
LuzhiyonggLuzhiyongg
authored
[PWGCF] cross check v2 and efficiency (#10276)
1 parent 3296839 commit efba710

File tree

2 files changed

+173
-81
lines changed

2 files changed

+173
-81
lines changed

‎PWGCF/Flow/Tasks/flowMc.cxx‎

Lines changed: 162 additions & 80 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
/// \brief QC of synthetic flow exercise
1616

1717
#include <CCDB/BasicCCDBManager.h>
18+
#include <vector>
1819
#include <string>
1920
#include "Framework/AnalysisDataModel.h"
2021
#include "Framework/AnalysisTask.h"
@@ -23,6 +24,7 @@
2324
#include "Framework/ASoAHelpers.h"
2425
#include "Framework/RunningWorkflowInfo.h"
2526
#include "Common/DataModel/TrackSelectionTables.h"
27+
#include "Common/Core/RecoDecay.h"
2628
#include "Common/Core/TrackSelection.h"
2729
#include "Common/Core/TrackSelectionDefaults.h"
2830
#include "Common/Core/trackUtilities.h"
@@ -35,6 +37,9 @@
3537
#include "GFW.h"
3638
#include "GFWCumulant.h"
3739
#include "GFWWeights.h"
40+
#include "FlowContainer.h"
41+
#include <TProfile.h>
42+
#include <TRandom3.h>
3843
#include <TPDGCode.h>
3944

4045
using namespace o2;
@@ -51,15 +56,25 @@ struct FlowMc {
5156
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
5257
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks")
5358
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks")
59+
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
60+
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks")
5461
O2_DEFINE_CONFIGURABLE(cfgFlowAcceptance, std::string, "", "CCDB path to acceptance object")
5562
O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object")
63+
O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth")
64+
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantEnabled, bool, false, "switch of calculating flow")
65+
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
5666

5767
ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""};
68+
ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"};
5869
ConfigurableAxis axisPhi{"axisPhi", {100, 0.0f, constants::math::TwoPI}, ""};
5970
ConfigurableAxis axisNch{"axisNch", {300, 0.0f, 3000.0f}, "Nch in |eta|<0.8"};
6071

6172
ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f}, "pt axis"};
6273

74+
// Cent vs IP
75+
TH1D* mCentVsIPTruth = nullptr;
76+
bool centVsIPTruthLoaded = false;
77+
6378
// Corrections
6479
TH1D* mEfficiency = nullptr;
6580
GFWWeights* mAcceptance = nullptr;
@@ -71,6 +86,14 @@ struct FlowMc {
7186
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
7287

7388
OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
89+
OutputObj<FlowContainer> fFCTrue{FlowContainer("FlowContainerTrue")};
90+
OutputObj<FlowContainer> fFCReco{FlowContainer("FlowContainerReco")};
91+
GFW* fGFWTrue = new GFW();
92+
GFW* fGFWReco = new GFW();
93+
TAxis* fPtAxis;
94+
std::vector<GFW::CorrConfig> corrconfigsTruth;
95+
std::vector<GFW::CorrConfig> corrconfigsReco;
96+
TRandom3* fRndm = new TRandom3(0);
7497

7598
void init(InitContext&)
7699
{
@@ -103,15 +126,61 @@ struct FlowMc {
103126
histos.add<TH2>("hEPVsPhi", "hEPVsPhi;Event Plane Angle; #varphi", HistType::kTH2D, {axisPhi, axisPhi});
104127
histos.add<TH2>("hPtNchGenerated", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
105128
histos.add<TH2>("hPtNchGlobal", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
129+
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
130+
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
106131

107-
if (cfgOutputNUAWeights) {
108-
o2::framework::AxisSpec axis = axisPt;
109-
int nPtBins = axis.binEdges.size() - 1;
110-
double* ptBins = &(axis.binEdges)[0];
132+
o2::framework::AxisSpec axis = axisPt;
133+
int nPtBins = axis.binEdges.size() - 1;
134+
double* ptBins = &(axis.binEdges)[0];
135+
fPtAxis = new TAxis(nPtBins, ptBins);
111136

137+
if (cfgOutputNUAWeights) {
112138
fWeights->setPtBins(nPtBins, ptBins);
113139
fWeights->init(true, false);
114140
}
141+
142+
if (cfgFlowCumulantEnabled) {
143+
TObjArray* oba = new TObjArray();
144+
oba->Add(new TNamed("ChFull22", "ChFull22"));
145+
for (auto i = 0; i < fPtAxis->GetNbins(); i++)
146+
oba->Add(new TNamed(Form("ChFull22_pt_%i", i + 1), "ChFull22_pTDiff"));
147+
oba->Add(new TNamed("Ch10Gap22", "Ch10Gap22"));
148+
for (auto i = 0; i < fPtAxis->GetNbins(); i++)
149+
oba->Add(new TNamed(Form("Ch10Gap22_pt_%i", i + 1), "Ch10Gap22_pTDiff"));
150+
fFCTrue->SetName("FlowContainerTrue");
151+
fFCTrue->SetXAxis(fPtAxis);
152+
fFCTrue->Initialize(oba, axisCentrality, cfgFlowCumulantNbootstrap);
153+
fFCReco->SetName("FlowContainerReco");
154+
fFCReco->SetXAxis(fPtAxis);
155+
fFCReco->Initialize(oba, axisCentrality, cfgFlowCumulantNbootstrap);
156+
delete oba;
157+
158+
fGFWTrue->AddRegion("full", -0.8, 0.8, 1, 1);
159+
fGFWTrue->AddRegion("refN10", -0.8, -0.5, 1, 1);
160+
fGFWTrue->AddRegion("refP10", 0.5, 0.8, 1, 1);
161+
fGFWTrue->AddRegion("poiN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 2);
162+
fGFWTrue->AddRegion("poifull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 2);
163+
fGFWTrue->AddRegion("olN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 4);
164+
fGFWTrue->AddRegion("olfull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 4);
165+
corrconfigsTruth.push_back(fGFWTrue->GetCorrelatorConfig("full {2 -2}", "ChFull22", kFALSE));
166+
corrconfigsTruth.push_back(fGFWTrue->GetCorrelatorConfig("poifull full | olfull {2 -2}", "ChFull22", kTRUE));
167+
corrconfigsTruth.push_back(fGFWTrue->GetCorrelatorConfig("refN10 {2} refP10 {-2}", "Ch10Gap22", kFALSE));
168+
corrconfigsTruth.push_back(fGFWTrue->GetCorrelatorConfig("poiN10 refN10 | olN10 {2} refP10 {-2}", "Ch10Gap22", kTRUE));
169+
fGFWTrue->CreateRegions();
170+
171+
fGFWReco->AddRegion("full", -0.8, 0.8, 1, 1);
172+
fGFWReco->AddRegion("refN10", -0.8, -0.5, 1, 1);
173+
fGFWReco->AddRegion("refP10", 0.5, 0.8, 1, 1);
174+
fGFWReco->AddRegion("poiN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 2);
175+
fGFWReco->AddRegion("poifull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 2);
176+
fGFWReco->AddRegion("olN10", -0.8, -0.5, 1 + fPtAxis->GetNbins(), 4);
177+
fGFWReco->AddRegion("olfull", -0.8, 0.8, 1 + fPtAxis->GetNbins(), 4);
178+
corrconfigsReco.push_back(fGFWReco->GetCorrelatorConfig("full {2 -2}", "ChFull22", kFALSE));
179+
corrconfigsReco.push_back(fGFWReco->GetCorrelatorConfig("poifull full | olfull {2 -2}", "ChFull22", kTRUE));
180+
corrconfigsReco.push_back(fGFWReco->GetCorrelatorConfig("refN10 {2} refP10 {-2}", "Ch10Gap22", kFALSE));
181+
corrconfigsReco.push_back(fGFWReco->GetCorrelatorConfig("poiN10 refN10 | olN10 {2} refP10 {-2}", "Ch10Gap22", kTRUE));
182+
fGFWReco->CreateRegions();
183+
}
115184
}
116185

117186
void loadCorrections(uint64_t timestamp)
@@ -152,6 +221,54 @@ struct FlowMc {
152221
return true;
153222
}
154223

224+
void fillFC(GFW* fGFW, bool isMCTruth, const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
225+
{
226+
double dnx, val;
227+
dnx = fGFW->Calculate(corrconf, 0, kTRUE).real();
228+
if (!corrconf.pTDif) {
229+
if (dnx == 0)
230+
return;
231+
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
232+
if (std::fabs(val) < 1) {
233+
if (isMCTruth)
234+
fFCTrue->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
235+
else
236+
fFCReco->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
237+
}
238+
return;
239+
}
240+
for (auto i = 1; i <= fPtAxis->GetNbins(); i++) {
241+
dnx = fGFW->Calculate(corrconf, i - 1, kTRUE).real();
242+
if (dnx == 0)
243+
continue;
244+
val = fGFW->Calculate(corrconf, i - 1, kFALSE).real() / dnx;
245+
if (std::fabs(val) < 1) {
246+
if (isMCTruth)
247+
fFCTrue->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
248+
else
249+
fFCReco->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
250+
}
251+
}
252+
return;
253+
}
254+
255+
void loadCentVsIPTruth(uint64_t timestamp)
256+
{
257+
if (centVsIPTruthLoaded)
258+
return;
259+
if (cfgCentVsIPTruth.value.empty() == false) {
260+
mCentVsIPTruth = ccdb->getForTimeStamp<TH1D>(cfgCentVsIPTruth, timestamp);
261+
if (mCentVsIPTruth)
262+
LOGF(info, "Loaded CentVsIPTruth weights from %s (%p)", cfgCentVsIPTruth.value.c_str(), (void*)mCentVsIPTruth);
263+
else
264+
LOGF(fatal, "Failed to load CentVsIPTruth weights from %s", cfgCentVsIPTruth.value.c_str());
265+
266+
centVsIPTruthLoaded = true;
267+
} else {
268+
LOGF(fatal, "when calculate flow, Cent Vs IP distribution must be provided");
269+
}
270+
}
271+
155272
using RecoTracks = soa::Join<aod::TracksIU, aod::TracksExtra>;
156273

157274
void process(aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join<aod::McParticles, aod::ParticlesToTracks> const& mcParticles, RecoTracks const&)
@@ -160,11 +277,12 @@ struct FlowMc {
160277
float imp = mcCollision.impactParameter();
161278
float evPhi = mcCollision.eventPlaneAngle();
162279
float vtxz = mcCollision.posZ();
163-
if (evPhi < 0)
164-
evPhi += constants::math::TwoPI;
280+
evPhi = RecoDecay::constrainAngle(evPhi);
165281

166282
int64_t nCh = 0;
167283
int64_t nChGlobal = 0;
284+
float centrality = 0;
285+
float lRandom = fRndm->Rndm();
168286
float weff = 1.;
169287
float wacc = 1.;
170288
auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
@@ -174,6 +292,12 @@ struct FlowMc {
174292
// event within range
175293
histos.fill(HIST("hImpactParameter"), imp);
176294
histos.fill(HIST("hEventPlaneAngle"), evPhi);
295+
if (cfgFlowCumulantEnabled) {
296+
loadCentVsIPTruth(bc.timestamp());
297+
centrality = mCentVsIPTruth->GetBinContent(mCentVsIPTruth->GetXaxis()->FindBin(imp));
298+
fGFWTrue->Clear();
299+
fGFWReco->Clear();
300+
}
177301

178302
for (auto const& mcParticle : mcParticles) {
179303
int pdgCode = std::abs(mcParticle.pdgCode());
@@ -204,13 +328,11 @@ struct FlowMc {
204328
continue;
205329

206330
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
207-
if (deltaPhi < 0)
208-
deltaPhi += constants::math::TwoPI;
209-
if (deltaPhi > constants::math::TwoPI)
210-
deltaPhi -= constants::math::TwoPI;
331+
deltaPhi = RecoDecay::constrainAngle(deltaPhi);
211332
histos.fill(HIST("hPtVsPhiGenerated"), deltaPhi, mcParticle.pt());
212333
histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt());
213334
histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal);
335+
histos.fill(HIST("hPtMCGen"), mcParticle.pt());
214336

215337
nCh++;
216338

@@ -241,10 +363,30 @@ struct FlowMc {
241363
}
242364

243365
bool withinPtRef = (cfgCutPtRefMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtRefMax); // within RF pT range
366+
bool withinPtPOI = (cfgCutPtPOIMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtPOIMax); // within POI pT range
244367
if (cfgOutputNUAWeights && withinPtRef)
245368
fWeights->fill(mcParticle.phi(), mcParticle.eta(), vtxz, mcParticle.pt(), 0, 0);
246369
if (!setCurrentParticleWeights(weff, wacc, mcParticle.phi(), mcParticle.eta(), mcParticle.pt(), vtxz))
247370
continue;
371+
372+
if (cfgFlowCumulantEnabled) {
373+
if (withinPtRef)
374+
fGFWTrue->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 1);
375+
if (withinPtPOI)
376+
fGFWTrue->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 2);
377+
if (withinPtPOI && withinPtRef)
378+
fGFWTrue->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 4);
379+
380+
if (validGlobal) {
381+
if (withinPtRef)
382+
fGFWReco->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 1);
383+
if (withinPtPOI)
384+
fGFWReco->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 2);
385+
if (withinPtPOI && withinPtRef)
386+
fGFWReco->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), wacc * weff, 4);
387+
}
388+
}
389+
248390
if (withinPtRef) {
249391
histos.fill(HIST("hEPVsPhiMC"), evPhi, mcParticle.phi());
250392
}
@@ -260,6 +402,7 @@ struct FlowMc {
260402
histos.fill(HIST("hPtVsPhiGlobal"), deltaPhi, mcParticle.pt(), wacc * weff);
261403
histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
262404
histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal);
405+
histos.fill(HIST("hPtMCGlobal"), mcParticle.pt());
263406
}
264407
// if any track present, fill
265408
if (validTrack)
@@ -271,79 +414,18 @@ struct FlowMc {
271414
if (validITSABTrack)
272415
histos.fill(HIST("hBVsPtVsPhiITSABTrack"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
273416
}
274-
}
275-
histos.fill(HIST("hNchVsImpactParameter"), imp, nCh);
276-
}
277-
278-
using LabeledCascades = soa::Join<aod::CascDataExt, aod::McCascLabels>;
279-
280-
void processCascade(aod::McParticle const& mcParticle, soa::SmallGroups<LabeledCascades> const& cascades, RecoTracks const&, aod::McCollisions const&)
281-
{
282-
auto mcCollision = mcParticle.mcCollision();
283-
float imp = mcCollision.impactParameter();
284-
285-
int pdgCode = std::abs(mcParticle.pdgCode());
286-
if (pdgCode != PDG_t::kXiMinus && pdgCode != PDG_t::kOmegaMinus)
287-
return;
288-
289-
if (!mcParticle.isPhysicalPrimary())
290-
return;
291-
if (std::fabs(mcParticle.eta()) > 0.8)
292-
return;
293-
294-
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
295-
if (deltaPhi < 0)
296-
deltaPhi += constants::math::TwoPI;
297-
if (deltaPhi > constants::math::TwoPI)
298-
deltaPhi -= constants::math::TwoPI;
299-
if (pdgCode == PDG_t::kXiMinus)
300-
histos.fill(HIST("hBVsPtVsPhiGeneratedXi"), imp, deltaPhi, mcParticle.pt());
301-
if (pdgCode == PDG_t::kOmegaMinus)
302-
histos.fill(HIST("hBVsPtVsPhiGeneratedOmega"), imp, deltaPhi, mcParticle.pt());
303-
304-
if (cascades.size() > 0) {
305-
if (pdgCode == PDG_t::kXiMinus)
306-
histos.fill(HIST("hBVsPtVsPhiGlobalXi"), imp, deltaPhi, mcParticle.pt());
307-
if (pdgCode == PDG_t::kOmegaMinus)
308-
histos.fill(HIST("hBVsPtVsPhiGlobalOmega"), imp, deltaPhi, mcParticle.pt());
309-
}
310-
}
311-
PROCESS_SWITCH(FlowMc, processCascade, "Process cascades", true);
312417

313-
using LabeledV0s = soa::Join<aod::V0Datas, aod::McV0Labels>;
314-
315-
void processV0s(aod::McParticle const& mcParticle, soa::SmallGroups<LabeledV0s> const& v0s, RecoTracks const&, aod::McCollisions const&)
316-
{
317-
auto mcCollision = mcParticle.mcCollision();
318-
float imp = mcCollision.impactParameter();
319-
320-
int pdgCode = std::abs(mcParticle.pdgCode());
321-
if (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0)
322-
return;
323-
324-
if (!mcParticle.isPhysicalPrimary())
325-
return;
326-
if (std::fabs(mcParticle.eta()) > 0.8)
327-
return;
328-
329-
float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle();
330-
if (deltaPhi < 0)
331-
deltaPhi += constants::math::TwoPI;
332-
if (deltaPhi > constants::math::TwoPI)
333-
deltaPhi -= constants::math::TwoPI;
334-
if (pdgCode == PDG_t::kK0Short)
335-
histos.fill(HIST("hBVsPtVsPhiGeneratedK0Short"), imp, deltaPhi, mcParticle.pt());
336-
if (pdgCode == PDG_t::kLambda0)
337-
histos.fill(HIST("hBVsPtVsPhiGeneratedLambda"), imp, deltaPhi, mcParticle.pt());
338-
339-
if (v0s.size() > 0) {
340-
if (pdgCode == PDG_t::kK0Short)
341-
histos.fill(HIST("hBVsPtVsPhiGlobalK0Short"), imp, deltaPhi, mcParticle.pt());
342-
if (pdgCode == PDG_t::kLambda0)
343-
histos.fill(HIST("hBVsPtVsPhiGlobalLambda"), imp, deltaPhi, mcParticle.pt());
418+
if (cfgFlowCumulantEnabled) {
419+
for (uint l_ind = 0; l_ind < corrconfigsTruth.size(); l_ind++) {
420+
fillFC(fGFWTrue, true, corrconfigsTruth.at(l_ind), centrality, lRandom);
421+
}
422+
for (uint l_ind = 0; l_ind < corrconfigsReco.size(); l_ind++) {
423+
fillFC(fGFWReco, false, corrconfigsReco.at(l_ind), centrality, lRandom);
424+
}
425+
}
344426
}
427+
histos.fill(HIST("hNchVsImpactParameter"), imp, nCh);
345428
}
346-
PROCESS_SWITCH(FlowMc, processV0s, "Process V0s", true);
347429
};
348430

349431
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)