Skip to content

Commit e0a2e80

Browse files
author
Maria Calmon Behling
committed
updated CCDB files, ratio to pions, more QA plots
1 parent 83c6e4d commit e0a2e80

2 files changed

Lines changed: 152 additions & 11 deletions

File tree

‎PWGLF/TableProducer/Nuspex/particleCompositionCorrection.cxx‎

Lines changed: 76 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -70,7 +70,7 @@ struct ParticleCompositionCorrection {
7070
Configurable<float> ptMaxCut{"ptMaxCut", 10.f, "pt max cut"};
7171
Configurable<bool> enableQAHistos{"enableQAHistos", true, "enable qa histograms showing the effect of the PCC"};
7272

73-
Configurable<std::string> ccdbBasePath{"ccdbBasePath", "/Users/m/makruger/", "ccdb directory contianing the particle fraction networks"};
73+
Configurable<std::string> ccdbBasePath{"ccdbBasePath", "/Users/m/mcalmonb/", "ccdb directory contianing the particle fraction networks"};
7474
Configurable<std::string> modelPathData{"modelPathData", "PCC/data/pp", "Path to the .onnx file containing the particle fractions in data"};
7575
Configurable<std::string> modelPathMC{"modelPathMC", "PCC/pythia/pp", "Path to the .onnx file containing the particle fractions in MC"};
7676

@@ -97,14 +97,16 @@ void ParticleCompositionCorrection::init(InitContext const&)
9797
return;
9898
}
9999
if (!ccdbBasePath.value.empty()) {
100-
ccdbApi.init("http://ccdb-test.cern.ch:8080");
100+
// ccdbApi.init("http://ccdb-test.cern.ch:8080");
101+
ccdbApi.init("http://alice-ccdb.cern.ch");
101102
static const int64_t dummyTimeStamp = 2;
102-
if (!ccdbApi.retrieveBlob(ccdbBasePath.value + modelPathData.value, modelPathData.value, {}, dummyTimeStamp, false, "ParticleFractions_Data.onnx") || !ccdbApi.retrieveBlob(ccdbBasePath.value + modelPathMC.value, modelPathMC.value, {}, dummyTimeStamp, false, "ParticleFractions_MC.onnx")) {
103+
if (!ccdbApi.retrieveBlob(ccdbBasePath.value + modelPathData.value, modelPathData.value, {}, dummyTimeStamp, false, "ParticleFractions_data.onnx") || !ccdbApi.retrieveBlob(ccdbBasePath.value + modelPathMC.value, modelPathMC.value, {}, dummyTimeStamp, false, "ParticleFractions_pythia.onnx")) {
103104
LOGP(fatal, "Could not download particle fraction networks!");
104105
}
105106
}
106-
particleFractionsData.initModel(modelPathData.value + "/ParticleFractions_Data.onnx", true);
107-
particleFractionsMC.initModel(modelPathMC.value + "/ParticleFractions_MC.onnx", true);
107+
108+
particleFractionsData.initModel(modelPathData.value + "/ParticleFractions_data.onnx", true);
109+
particleFractionsMC.initModel(modelPathMC.value + "/ParticleFractions_pythia.onnx", true);
108110

109111
if (enableQAHistos) {
110112
std::vector<double> ptBinEdges = {0.15, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75,
@@ -113,19 +115,49 @@ void ParticleCompositionCorrection::init(InitContext const&)
113115
6.0, 6.5, 7.0, 8.0, 9.0, 10.0};
114116
const AxisSpec ptAxis{ptBinEdges, "#it{p}_{T} (GeV/#it{c})", "pt"};
115117

116-
histos.add("frac/data/pion", "", kTProfile, {ptAxis});
118+
const int maxMult = 100;
119+
const int nBinsMult = maxMult + 1;
120+
const AxisSpec multAxis = {nBinsMult, -0.5, nBinsMult - 0.5, "#it{N}_{ch}", "mult"};
121+
122+
// Multiplicity distributions: fiducial vs INELg0
123+
histos.add("multDist_INELg0", "", kTH1D, {multAxis});
124+
histos.add("multDist_fid", "", kTH1D, {multAxis});
125+
histos.add("multDist_fidVsINELg0", "", kTH2D, {multAxis, multAxis});
126+
127+
// Fractions data 1D
128+
// histos.add("frac/data/pion", "", kTProfile, {ptAxis});
117129
histos.add("frac/data/kaon", "", kTProfile, {ptAxis});
118130
histos.add("frac/data/proton", "", kTProfile, {ptAxis});
119131
histos.add("frac/data/sigma", "", kTProfile, {ptAxis});
132+
133+
// Fractions data 2D
134+
// histos.add("frac/data/pion_mult", "", kTProfile2D, {multAxis, ptAxis});
135+
histos.add("frac/data/kaon_mult", "", kTProfile2D, {multAxis, ptAxis});
136+
histos.add("frac/data/proton_mult", "", kTProfile2D, {multAxis, ptAxis});
137+
histos.add("frac/data/sigma_mult", "", kTProfile2D, {multAxis, ptAxis});
138+
139+
// Fractions MC
120140
histos.addClone("frac/data/", "frac/mc/");
121141

122-
histos.add("weight/pion", "", kTProfile, {ptAxis});
142+
// Weights Data / MC 1D
143+
// histos.add("weight/pion", "", kTProfile, {ptAxis});
123144
histos.add("weight/kaon", "", kTProfile, {ptAxis});
124145
histos.add("weight/proton", "", kTProfile, {ptAxis});
125146
histos.add("weight/sigma", "", kTProfile, {ptAxis});
126147

148+
// Weights Data / MC 2D
149+
// histos.add("weight/pion_mult", "", kTProfile2D, {multAxis,ptAxis});
150+
histos.add("weight/kaon_mult", "", kTProfile2D, {multAxis, ptAxis});
151+
histos.add("weight/proton_mult", "", kTProfile2D, {multAxis, ptAxis});
152+
histos.add("weight/sigma_mult", "", kTProfile2D, {multAxis, ptAxis});
153+
154+
// Weights Secondaries 1D
127155
histos.add("weight/secDec", "", kTProfile, {ptAxis});
128156
histos.add("weight/secMat", "", kTProfile, {ptAxis});
157+
158+
// Weights Secondaries 2D
159+
histos.add("weight/secDec_mult", "", kTProfile2D, {multAxis, ptAxis});
160+
histos.add("weight/secMat_mult", "", kTProfile2D, {multAxis, ptAxis});
129161
}
130162
}
131163

@@ -143,6 +175,14 @@ std::tuple<float, float, float> ParticleCompositionCorrection::getWeights(aod::M
143175
}
144176
auto absPDGCode = std::abs(particle.pdgCode());
145177
// translate abs PDG code to PID variable of neural networks (0: pion, 1: kaon, 2: proton, 3: sigma)
178+
179+
// std::cout << "absPDGCode: " << absPDGCode << '\n';
180+
181+
if (absPDGCode == PDG_t::kPiPlus || absPDGCode == PDG_t::kPi0) {
182+
// std::cout << "same PDG code as Pi +- (" << PDG_t::kPiPlus <<") or Pi0 (" << PDG_t::kPi0 <<")--> skipping " << '\n';
183+
return noWeights;
184+
}
185+
146186
static const std::map<int, float> mapPID = {
147187
{PDG_t::kPiPlus, 0.f},
148188
{PDG_t::kPi0, 0.f},
@@ -175,25 +215,37 @@ std::tuple<float, float, float> ParticleCompositionCorrection::getWeights(aod::M
175215
storedWeights[particle.index()] = weights;
176216
}
177217
if (enableQAHistos && particle.isPhysicalPrimary() && std::abs(particle.eta()) < 0.8) { // o2-linter: disable=magic-number (usual range of charged-partilce measurements)
178-
if (iterMapPID->first == PDG_t::kPiPlus) {
218+
/* if (iterMapPID->first == PDG_t::kPiPlus) {
179219
histos.fill(HIST("frac/data/pion"), pt, fracData);
180220
histos.fill(HIST("frac/mc/pion"), pt, fracMC);
221+
histos.fill(HIST("frac/data/pion_mult"), dNdEta, pt, fracData);
222+
histos.fill(HIST("frac/mc/pion_mult"), dNdEta, pt, fracMC);
181223
histos.fill(HIST("weight/pion"), pt, weight);
182-
}
224+
histos.fill(HIST("weight/pion_mult"), dNdEta, pt, weight);
225+
} */
183226
if (iterMapPID->first == PDG_t::kKPlus) {
184227
histos.fill(HIST("frac/data/kaon"), pt, fracData);
185228
histos.fill(HIST("frac/mc/kaon"), pt, fracMC);
229+
histos.fill(HIST("frac/data/kaon_mult"), dNdEta, pt, fracData);
230+
histos.fill(HIST("frac/mc/kaon_mult"), dNdEta, pt, fracMC);
186231
histos.fill(HIST("weight/kaon"), pt, weight);
232+
histos.fill(HIST("weight/kaon_mult"), dNdEta, pt, weight);
187233
}
188234
if (iterMapPID->first == PDG_t::kProton) {
189235
histos.fill(HIST("frac/data/proton"), pt, fracData);
190236
histos.fill(HIST("frac/mc/proton"), pt, fracMC);
237+
histos.fill(HIST("frac/data/proton_mult"), dNdEta, pt, fracData);
238+
histos.fill(HIST("frac/mc/proton_mult"), dNdEta, pt, fracMC);
191239
histos.fill(HIST("weight/proton"), pt, weight);
240+
histos.fill(HIST("weight/proton_mult"), dNdEta, pt, weight);
192241
}
193242
if (iterMapPID->first == PDG_t::kSigmaPlus || iterMapPID->first == PDG_t::kSigmaMinus) {
194243
histos.fill(HIST("frac/data/sigma"), pt, fracData);
195244
histos.fill(HIST("frac/mc/sigma"), pt, fracMC);
245+
histos.fill(HIST("frac/data/sigma_mult"), dNdEta, pt, fracData);
246+
histos.fill(HIST("frac/mc/sigma_mult"), dNdEta, pt, fracMC);
196247
histos.fill(HIST("weight/sigma"), pt, weight);
248+
histos.fill(HIST("weight/sigma_mult"), dNdEta, pt, weight);
197249
}
198250
}
199251
return weights;
@@ -218,8 +270,10 @@ std::tuple<float, float, float> ParticleCompositionCorrection::getWeights(aod::M
218270
if (pdgParticle && pdgParticle->Charge() != 0.) {
219271
if (particle.getProcess() == TMCProcess::kPDecay) {
220272
histos.fill(HIST("weight/secDec"), particle.pt(), weight);
273+
histos.fill(HIST("weight/secDec_mult"), dNdEta, particle.pt(), weight);
221274
} else if (particle.getProcess() == TMCProcess::kPHInhelastic || particle.getProcess() == TMCProcess::kPHadronic || particle.getProcess() == TMCProcess::kPHElastic) {
222275
histos.fill(HIST("weight/secMat"), particle.pt(), weight);
276+
histos.fill(HIST("weight/secMat_mult"), dNdEta, particle.pt(), weight);
223277
}
224278
}
225279
}
@@ -232,6 +286,8 @@ void ParticleCompositionCorrection::process(aod::McCollisions::iterator const&,
232286
{
233287
// determine dNdEta of the collision
234288
float dNdEta = 0.f;
289+
float dNdEtafid = 0.f;
290+
235291
for (const auto& particle : particles) {
236292
if (!particle.isPhysicalPrimary()) {
237293
continue;
@@ -240,10 +296,19 @@ void ParticleCompositionCorrection::process(aod::McCollisions::iterator const&,
240296
if (!pdgParticle || pdgParticle->Charge() == 0.) {
241297
continue;
242298
}
243-
if (std::abs(particle.eta()) >= 0.5) { // o2-linter: disable=magic-number (particle density at mid-rapidity)
299+
if (std::abs(particle.eta()) >= 0.8) { // o2-linter: disable=magic-number (particle density at mid-rapidity)
244300
continue;
245301
}
246-
++dNdEta;
302+
if (std::abs(particle.eta()) < 0.5) { // o2-linter: disable=magic-number (particle density at mid-rapidity)
303+
++dNdEta;
304+
}
305+
++dNdEtafid;
306+
}
307+
308+
if (dNdEtafid > 0.f) {
309+
histos.fill(HIST("multDist_INELg0"), dNdEta);
310+
histos.fill(HIST("multDist_fid"), dNdEtafid);
311+
histos.fill(HIST("multDist_fidVsINELg0"), dNdEta, dNdEtafid);
247312
}
248313

249314
std::map<int32_t, std::tuple<float, float, float>> storedWeights;

‎PWGLF/Tasks/Nuspex/pccQa.cxx‎

Lines changed: 76 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -22,11 +22,27 @@
2222
#include <Framework/AnalysisDataModel.h>
2323
#include <Framework/AnalysisHelpers.h>
2424
#include <Framework/AnalysisTask.h>
25+
#include <Framework/ConfigParamSpec.h>
2526
#include <Framework/Configurable.h>
2627
#include <Framework/HistogramRegistry.h>
2728
#include <Framework/HistogramSpec.h>
2829
#include <Framework/InitContext.h>
2930
#include <Framework/O2DatabasePDGPlugin.h>
31+
32+
#include <iostream>
33+
34+
// Workaround local
35+
/* void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
36+
{
37+
workflowOptions.push_back(
38+
{"doMC",
39+
o2::framework::VariantType::Bool,
40+
false,
41+
{"Use MC info"}}
42+
);
43+
} */
44+
45+
#include "Framework/Logger.h"
3046
#include <Framework/runDataProcessing.h>
3147

3248
#include <TMCProcess.h>
@@ -43,6 +59,8 @@ struct PccQa {
4359

4460
static constexpr float MaxVtxZ = 10.f;
4561

62+
Configurable<uint32_t> maxMult{"maxMult", 100, "max multiplicity"};
63+
4664
void init(InitContext const&);
4765

4866
template <bool IS_MC, typename C, typename T>
@@ -67,34 +85,81 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
6785
return WorkflowSpec{adaptAnalysisTask<PccQa>(cfgc)};
6886
}
6987

88+
// Workaround local:
89+
/* WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
90+
{
91+
const bool doMC = cfgc.options().get<bool>("doMC");
92+
93+
LOGP(info, "doMC = {}", static_cast<bool>(doMC));
94+
95+
SetDefaultProcesses processes;
96+
processes.map = {
97+
{"processData", !doMC},
98+
{"processMC", doMC}
99+
};
100+
101+
return WorkflowSpec{
102+
adaptAnalysisTask<PccQa>(cfgc, processes)
103+
};
104+
} */
105+
70106
void PccQa::init(InitContext const&)
71107
{
108+
109+
LOGP(info, "processData = {}", static_cast<bool>(doprocessData));
110+
LOGP(info, "processMC = {}", static_cast<bool>(doprocessMC));
111+
72112
histos.add("eventCounter", "", kTH1D, {{1, 0.5, 1.5}});
73113
const AxisSpec dcaAxis = {1000, -1., 1., "#it{DCA}_{xy}", "dca"};
74114
std::vector<double> ptBinEdges = {0.2, 0.4, 0.6, 0.8, 1.0, 1.2, 1.4, 1.6, 1.8, 2.0, 2.5, 3.0};
75115
const AxisSpec ptAxis{ptBinEdges, "#it{p}_{T} (GeV/#it{c})", "pt"};
76116

117+
const int nBinsMult = 30;
118+
const AxisSpec nchAxis = {nBinsMult, 0., static_cast<double>(maxMult), "#it{N}_{ch}", "nch"};
119+
77120
histos.add("DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis});
121+
histos.add("DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
78122

79123
if (doprocessMC) {
80124
histos.add("DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis});
125+
histos.add("DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
126+
81127
histos.add("prim/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis});
128+
histos.add("prim/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
129+
130+
histos.add("prim/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
82131
histos.add("prim/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis});
132+
83133
histos.add("sec/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis});
134+
histos.add("sec/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
135+
84136
histos.add("sec/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis});
137+
histos.add("sec/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
138+
85139
histos.add("sec/dec/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis});
140+
histos.add("sec/dec/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
141+
86142
histos.add("sec/dec/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis});
143+
histos.add("sec/dec/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
144+
87145
histos.add("sec/mat/DCAxyVsPt", "", kTH2D, {ptAxis, dcaAxis});
146+
histos.add("sec/mat/DCAxyVsPtVsMult", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
147+
88148
histos.add("sec/mat/DCAxyVsPt_weighted", "", kTH2D, {ptAxis, dcaAxis});
149+
histos.add("sec/mat/DCAxyVsPtVsMult_weighted", "", kTH3D, {nchAxis, ptAxis, dcaAxis});
89150
}
90151
}
91152

92153
void PccQa::processData(CollisionTableData::iterator const& collision, TrackTableData const& tracks)
93154
{
155+
std::cout << "using processData function: " << '\n';
94156
processMeas<false>(collision, tracks);
95157
}
96158
void PccQa::processMC(CollisionTableMCTrue::iterator const&, TrackTableMC const& tracks, CollisionTableMC const& collisions, ParticleTableMC const&)
97159
{
160+
161+
std::cout << "using processMC function: " << '\n';
162+
98163
for (const auto& collision : collisions) {
99164
auto curTracks = tracks.sliceBy(perCollision, collision.globalIndex());
100165
processMeas<true>(collision, curTracks);
@@ -105,6 +170,7 @@ void PccQa::processMC(CollisionTableMCTrue::iterator const&, TrackTableMC const&
105170
template <bool IS_MC, typename C, typename T>
106171
void PccQa::processMeas(const C& collision, const T& tracks)
107172
{
173+
108174
if ((std::abs(collision.posZ()) > MaxVtxZ) || !collision.sel8()) {
109175
return;
110176
}
@@ -115,6 +181,7 @@ void PccQa::processMeas(const C& collision, const T& tracks)
115181
continue;
116182
}
117183
histos.fill(HIST("DCAxyVsPt"), track.pt(), track.dcaXY());
184+
histos.fill(HIST("DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY());
118185

119186
if constexpr (IS_MC) {
120187
if (!track.has_mcParticle()) {
@@ -123,19 +190,28 @@ void PccQa::processMeas(const C& collision, const T& tracks)
123190
const auto& particle = track.template mcParticle_as<ParticleTableMC>();
124191

125192
histos.fill(HIST("DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight());
193+
histos.fill(HIST("DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight());
126194

127195
if (particle.isPhysicalPrimary()) {
128196
histos.fill(HIST("prim/DCAxyVsPt"), track.pt(), track.dcaXY());
129197
histos.fill(HIST("prim/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight());
198+
histos.fill(HIST("prim/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY());
199+
histos.fill(HIST("prim/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight());
130200
} else {
131201
histos.fill(HIST("sec/DCAxyVsPt"), track.pt(), track.dcaXY());
132202
histos.fill(HIST("sec/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight());
203+
histos.fill(HIST("sec/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY());
204+
histos.fill(HIST("sec/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight());
133205
if (particle.getProcess() == TMCProcess::kPDecay) {
134206
histos.fill(HIST("sec/dec/DCAxyVsPt"), track.pt(), track.dcaXY());
135207
histos.fill(HIST("sec/dec/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight());
208+
histos.fill(HIST("sec/dec/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY());
209+
histos.fill(HIST("sec/dec/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight());
136210
} else if (particle.getProcess() == TMCProcess::kPHInhelastic || particle.getProcess() == TMCProcess::kPHadronic || particle.getProcess() == TMCProcess::kPHElastic) {
137211
histos.fill(HIST("sec/mat/DCAxyVsPt"), track.pt(), track.dcaXY());
138212
histos.fill(HIST("sec/mat/DCAxyVsPt_weighted"), track.pt(), track.dcaXY(), particle.pccWeight());
213+
histos.fill(HIST("sec/mat/DCAxyVsPtVsMult"), tracks.size(), track.pt(), track.dcaXY());
214+
histos.fill(HIST("sec/mat/DCAxyVsPtVsMult_weighted"), tracks.size(), track.pt(), track.dcaXY(), particle.pccWeight());
139215
}
140216
}
141217
}

0 commit comments

Comments
 (0)