Skip to content

Commit bc014f6

Browse files
njacazionjacazio
authored
[PWG-LF] add MC task (#6419)
1 parent 1d3dfcd commit bc014f6

3 files changed

Lines changed: 517 additions & 0 deletions

File tree

‎PWGLF/Tasks/QC/CMakeLists.txt‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -108,3 +108,8 @@ o2physics_add_dpl_workflow(tracked-cascade-properties
108108
SOURCES tracked_cascade_properties.cxx
109109
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
110110
COMPONENT_NAME Analysis)
111+
112+
o2physics_add_dpl_workflow(mc-particle-predictions
113+
SOURCES mcParticlePrediction.cxx
114+
PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
115+
COMPONENT_NAME Analysis)

‎PWGLF/Tasks/QC/mcParticlePrediction.cxx‎

Lines changed: 192 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,192 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
///
13+
/// \file mcParticlePrediction.cxx
14+
/// \author Nicolò Jacazio nicolo.jacazio@cern.ch
15+
/// \author Francesca Ercolessi francesca.ercolessi@cern.ch
16+
/// \brief Task to build the predictions from the models based on the generated particles
17+
///
18+
19+
#include "Framework/runDataProcessing.h"
20+
#include "Framework/AnalysisTask.h"
21+
#include "Framework/HistogramRegistry.h"
22+
#include "Framework/StaticFor.h"
23+
#include "Framework/O2DatabasePDGPlugin.h"
24+
#include "PWGLF/Utils/mcParticle.h"
25+
#include "PWGLF/Utils/inelGt.h"
26+
27+
#include "TPDGCode.h"
28+
29+
using namespace o2;
30+
using namespace o2::framework;
31+
32+
static const std::vector<std::string> parameterNames{"Enable"};
33+
static constexpr int nParameters = 1;
34+
static const int defaultParameters[o2::pwglf::PIDExtended::NIDsTot][nParameters]{{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}};
35+
bool enabledArray[o2::pwglf::PIDExtended::NIDsTot];
36+
37+
// Histograms
38+
const int FT0A = 0;
39+
const int FT0C = 1;
40+
const int FT0AC = 2;
41+
const int FV0A = 3;
42+
const int FDDA = 4;
43+
const int FDDC = 5;
44+
const int FDDAC = 6;
45+
const int ZNA = 7;
46+
const int ZNC = 8;
47+
// const int ZEM1 = 9;
48+
// const int ZEM2 = 10;
49+
// const int ZPA = 11;
50+
// const int ZPC = 12;
51+
const int nEstimators = 13;
52+
53+
const char* estimatorNames[nEstimators] = {"FT0A",
54+
"FT0C",
55+
"FT0AC",
56+
"FV0A",
57+
"FDDA",
58+
"FDDC",
59+
"FDDAC",
60+
"ZNA",
61+
"ZNC",
62+
"ZEM1",
63+
"ZEM2",
64+
"ZPA",
65+
"ZPC"};
66+
67+
std::array<std::array<std::shared_ptr<TH2>, o2::pwglf::PIDExtended::NIDsTot>, nEstimators> hpt;
68+
std::array<std::array<std::shared_ptr<TH1>, o2::pwglf::PIDExtended::NIDsTot>, nEstimators> hyield;
69+
70+
struct mcParticlePrediction {
71+
72+
// Histograms
73+
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
74+
ConfigurableAxis binsPt{"binsPt", {100, 0, 10}, "Binning of the Pt axis"};
75+
ConfigurableAxis binsMultiplicity{"binsMultiplicity", {100, 0, 1000}, "Binning of the Multiplicity axis"};
76+
Configurable<LabeledArray<int>> enabledSpecies{"enabledSpecies",
77+
{defaultParameters[0], o2::pwglf::PIDExtended::NIDsTot, nParameters, o2::pwglf::PIDExtended::arrayNames(), parameterNames},
78+
"Bethe Bloch parameters"};
79+
Configurable<bool> selectInelGt0{"selectInelGt0", true, "Select only inelastic events"};
80+
Service<o2::framework::O2DatabasePDG> pdgDB;
81+
82+
void init(o2::framework::InitContext&)
83+
{
84+
histos.add("collisions", "collisions", kTH1D, {{10, 0, 10}});
85+
auto h = histos.add<TH1>("particles", "particles", kTH1D, {{o2::pwglf::PIDExtended::NIDsTot, -0.5, -0.5 + o2::pwglf::PIDExtended::NIDsTot}});
86+
for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) {
87+
h->GetXaxis()->SetBinLabel(i + 1, o2::pwglf::PIDExtended::getName(i));
88+
}
89+
h = histos.add<TH1>("particlesPrim", "particlesPrim", kTH1D, {{o2::pwglf::PIDExtended::NIDsTot, -0.5, -0.5 + o2::pwglf::PIDExtended::NIDsTot}});
90+
for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) {
91+
h->GetXaxis()->SetBinLabel(i + 1, o2::pwglf::PIDExtended::getName(i));
92+
}
93+
94+
const AxisSpec axisPt{binsPt, "#it{p}_{T} (GeV/#it{c})"};
95+
const AxisSpec axisMultiplicity{binsMultiplicity, "Mult"};
96+
// FT0
97+
histos.add("multiplicity/FT0A", "FT0A", kTH1D, {axisMultiplicity});
98+
histos.add("multiplicity/FT0C", "FT0C", kTH1D, {axisMultiplicity});
99+
histos.add("multiplicity/FT0AC", "FT0AC", kTH1D, {axisMultiplicity});
100+
101+
// FV0
102+
histos.add("multiplicity/FV0A", "FV0A", kTH1D, {axisMultiplicity});
103+
104+
// FDD
105+
histos.add("multiplicity/FDDA", "FDDA", kTH1D, {axisMultiplicity});
106+
histos.add("multiplicity/FDDC", "FDDC", kTH1D, {axisMultiplicity});
107+
histos.add("multiplicity/FDDAC", "FDDAC", kTH1D, {axisMultiplicity});
108+
109+
// ZDC
110+
histos.add("multiplicity/ZNA", "ZNA", kTH1D, {axisMultiplicity});
111+
histos.add("multiplicity/ZNC", "ZNC", kTH1D, {axisMultiplicity});
112+
histos.add("multiplicity/ZEM1", "ZEM1", kTH1D, {axisMultiplicity});
113+
histos.add("multiplicity/ZEM2", "ZEM2", kTH1D, {axisMultiplicity});
114+
histos.add("multiplicity/ZPA", "ZPA", kTH1D, {axisMultiplicity});
115+
histos.add("multiplicity/ZPC", "ZPC", kTH1D, {axisMultiplicity});
116+
117+
// ITS
118+
histos.add("multiplicity/ITS", "ITS", kTH1D, {axisMultiplicity});
119+
120+
for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) {
121+
if (enabledSpecies->get(o2::pwglf::PIDExtended::getName(i), "Enable") != 1) {
122+
enabledArray[i] = false;
123+
continue;
124+
}
125+
enabledArray[i] = true;
126+
for (int j = 0; j < nEstimators; j++) {
127+
hpt[j][i] = histos.add<TH2>(Form("pt/%s/%s", estimatorNames[j], o2::pwglf::PIDExtended::getName(i)), o2::pwglf::PIDExtended::getName(i), kTH2D, {binsPt, axisMultiplicity});
128+
hyield[j][i] = histos.add<TH1>(Form("yield/%s/%s", estimatorNames[j], o2::pwglf::PIDExtended::getName(i)), o2::pwglf::PIDExtended::getName(i), kTH1D, {axisMultiplicity});
129+
}
130+
}
131+
}
132+
133+
int countInAcceptance(const aod::McParticles& mcParticles, const float etamin, const float etamax)
134+
{
135+
// static_assert(etamin < etamax, "etamin must be smaller than etamax");
136+
int counter = 0;
137+
for (const auto& particle : mcParticles) {
138+
if (particle.eta() > etamin && particle.eta() < etamax) {
139+
counter++;
140+
}
141+
}
142+
return counter;
143+
}
144+
145+
int countFT0A(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 3.5f, 4.9f); }
146+
int countFT0C(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -3.3f, -2.1f); }
147+
int countFV0A(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 2.2f, 5.1f); }
148+
int countFDDA(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 4.9f, 6.3f); }
149+
int countFDDC(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -7.f, -4.9f); }
150+
int countZNA(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 8.8f, 100.f); }
151+
int countZNC(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -100.f, -8.8f); }
152+
153+
void process(aod::McCollision const& mcCollision,
154+
aod::McParticles const& mcParticles)
155+
{
156+
if (selectInelGt0.value && !o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) {
157+
return;
158+
}
159+
160+
histos.fill(HIST("collisions"), 0.5);
161+
const int nFT0A = countFT0A(mcParticles);
162+
const int nFT0C = countFT0C(mcParticles);
163+
const int nFV0A = countFV0A(mcParticles);
164+
const int nFDDA = countFDDA(mcParticles);
165+
const int nFDDC = countFDDC(mcParticles);
166+
const int nZNA = countZNA(mcParticles);
167+
const int nZNC = countZNC(mcParticles);
168+
169+
for (const auto& particle : mcParticles) {
170+
particle.pdgCode();
171+
const auto id = o2::pwglf::PIDExtended::pdgToId(particle);
172+
if (id < 0) {
173+
continue;
174+
}
175+
if (!enabledArray[id]) {
176+
continue;
177+
}
178+
hpt[FT0A][id]->Fill(particle.pt(), nFT0A);
179+
hpt[FT0C][id]->Fill(particle.pt(), nFT0C);
180+
hpt[FT0AC][id]->Fill(particle.pt(), nFT0A + nFT0C);
181+
182+
hpt[FV0A][id]->Fill(particle.pt(), nFV0A);
183+
hpt[FDDA][id]->Fill(particle.pt(), nFDDA);
184+
hpt[FDDC][id]->Fill(particle.pt(), nFDDC);
185+
hpt[FDDAC][id]->Fill(particle.pt(), nFDDA + nFDDC);
186+
hpt[ZNA][id]->Fill(particle.pt(), nZNA);
187+
hpt[ZNC][id]->Fill(particle.pt(), nZNC);
188+
}
189+
}
190+
};
191+
192+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<mcParticlePrediction>(cfgc)}; }

0 commit comments

Comments
 (0)