Skip to content

Commit a2c384a

Browse files
njacazionjacazio
committed
Add first task of the LUT maker from fullsim
- Add single particle QA task - Add 2D efficiency to efficiency QA task - Add possibility to select charge depending on PDG code
1 parent 269060e commit a2c384a

5 files changed

Lines changed: 349 additions & 10 deletions

File tree

‎Analysis/ALICE3/CMakeLists.txt‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -25,6 +25,16 @@ o2_add_dpl_workflow(alice3-qa-multiplicity
2525
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
2626
COMPONENT_NAME Analysis)
2727

28+
o2_add_dpl_workflow(alice3-qa-singleparticle
29+
SOURCES src/alice3-qa-singleparticle.cxx
30+
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
31+
COMPONENT_NAME Analysis)
32+
33+
o2_add_dpl_workflow(alice3-lutmaker
34+
SOURCES src/alice3-lutmaker.cxx
35+
PUBLIC_LINK_LIBRARIES O2::Framework O2::AnalysisDataModel O2::AnalysisCore
36+
COMPONENT_NAME Analysis)
37+
2838
o2_add_dpl_workflow(alice3-pid-rich-qa
2939
SOURCES src/pidRICHqa.cxx
3040
PUBLIC_LINK_LIBRARIES O2::AnalysisDataModel O2::AnalysisCore O2::ALICE3Analysis

‎Analysis/ALICE3/src/alice3-lutmaker.cxx‎

Lines changed: 152 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,152 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
11+
/// \author Nicolo' Jacazio <nicolo.jacazio@cern.ch>, CERN
12+
/// \brief Task to extract LUTs for the fast simulation from full simulation
13+
/// \since 27/04/2021
14+
15+
// O2 includes
16+
#include "Framework/AnalysisTask.h"
17+
#include "AnalysisCore/MC.h"
18+
#include "ReconstructionDataFormats/Track.h"
19+
20+
using namespace o2;
21+
using namespace framework;
22+
using namespace framework::expressions;
23+
24+
void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
25+
{
26+
std::vector<ConfigParamSpec> options{
27+
{"lul-el", VariantType::Int, 1, {"LUT input for the Electron PDG code"}},
28+
{"lut-mu", VariantType::Int, 1, {"LUT input for the Muon PDG code"}},
29+
{"lut-pi", VariantType::Int, 1, {"LUT input for the Pion PDG code"}},
30+
{"lut-ka", VariantType::Int, 1, {"LUT input for the Kaon PDG code"}},
31+
{"lut-pr", VariantType::Int, 1, {"LUT input for the Proton PDG code"}}};
32+
std::swap(workflowOptions, options);
33+
}
34+
35+
#include "Framework/runDataProcessing.h"
36+
37+
template <o2::track::pid_constants::ID particle>
38+
struct Alice3LutMaker {
39+
static constexpr PDG_t PDGs[5] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton};
40+
static_assert(particle < 5 && "Maximum of particles reached");
41+
static constexpr int pdg = PDGs[particle];
42+
// Configurable<int> pdg{"pdg", 2212, "PDG code of the particle of interest"};
43+
Configurable<bool> selPrim{"sel-prim", false, "If true selects primaries, if not select all particles"};
44+
Configurable<float> etaMin{"eta-min", -3.f, "Lower limit in eta"};
45+
Configurable<float> etaMax{"eta-max", 3.f, "Upper limit in eta"};
46+
Configurable<float> ptMin{"pt-min", 0.f, "Lower limit in pT"};
47+
Configurable<float> ptMax{"pt-max", 100.f, "Upper limit in pT"};
48+
Configurable<int> ptBins{"pt-bins", 1000, "Number of pT bins"};
49+
Configurable<int> logPt{"log-pt", 0, "Flag to use a logarithmic pT axis"};
50+
Configurable<int> etaBins{"eta-bins", 500, "Number of eta bins"};
51+
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
52+
53+
void init(InitContext&)
54+
{
55+
const TString commonTitle = Form(" PDG %i;#it{p}_{T};#eta", pdg);
56+
AxisSpec axisPt{ptBins, ptMin, ptMax};
57+
AxisSpec axisEta{etaBins, etaMin, etaMax};
58+
59+
histos.add("pt", Form("pt PDG %i;#it{p}_{T}", pdg), kTH1F, {axisPt});
60+
histos.add("eta", Form("eta PDG %i;#eta", pdg), kTH1F, {axisEta});
61+
histos.add("CovMat_cYY", "cYY" + commonTitle, kTProfile2D, {axisPt, axisEta});
62+
histos.add("CovMat_cZY", "cZY" + commonTitle, kTProfile2D, {axisPt, axisEta});
63+
histos.add("CovMat_cZZ", "cZZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
64+
histos.add("CovMat_cSnpY", "cSnpY" + commonTitle, kTProfile2D, {axisPt, axisEta});
65+
histos.add("CovMat_cSnpZ", "cSnpZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
66+
histos.add("CovMat_cSnpSnp", "cSnpSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
67+
histos.add("CovMat_cTglY", "cTglY" + commonTitle, kTProfile2D, {axisPt, axisEta});
68+
histos.add("CovMat_cTglZ", "cTglZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
69+
histos.add("CovMat_cTglSnp", "cTglSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
70+
histos.add("CovMat_cTglTgl", "cTglTgl" + commonTitle, kTProfile2D, {axisPt, axisEta});
71+
histos.add("CovMat_c1PtY", "c1PtY" + commonTitle, kTProfile2D, {axisPt, axisEta});
72+
histos.add("CovMat_c1PtZ", "c1PtZ" + commonTitle, kTProfile2D, {axisPt, axisEta});
73+
histos.add("CovMat_c1PtSnp", "c1PtSnp" + commonTitle, kTProfile2D, {axisPt, axisEta});
74+
histos.add("CovMat_c1PtTgl", "c1PtTgl" + commonTitle, kTProfile2D, {axisPt, axisEta});
75+
histos.add("CovMat_c1Pt21Pt2", "c1Pt21Pt2" + commonTitle, kTProfile2D, {axisPt, axisEta});
76+
77+
histos.add("Efficiency", "Efficiency" + commonTitle, kTProfile2D, {axisPt, axisEta});
78+
}
79+
80+
void process(const soa::Join<aod::Tracks, aod::TracksCov, aod::McTrackLabels>& tracks,
81+
const aod::McParticles& mcParticles)
82+
{
83+
std::vector<int64_t> recoTracks(tracks.size());
84+
int ntrks = 0;
85+
86+
for (const auto& track : tracks) {
87+
const auto mcParticle = track.mcParticle();
88+
if (mcParticle.pdgCode() != pdg) {
89+
continue;
90+
}
91+
if (selPrim.value && !MC::isPhysicalPrimary(mcParticles, mcParticle)) { // Requiring is physical primary
92+
continue;
93+
}
94+
95+
recoTracks[ntrks++] = mcParticle.globalIndex();
96+
97+
histos.fill(HIST("pt"), mcParticle.pt());
98+
histos.fill(HIST("eta"), mcParticle.eta());
99+
histos.fill(HIST("CovMat_cYY"), mcParticle.pt(), mcParticle.eta(), track.cYY());
100+
histos.fill(HIST("CovMat_cZY"), mcParticle.pt(), mcParticle.eta(), track.cZY());
101+
histos.fill(HIST("CovMat_cZZ"), mcParticle.pt(), mcParticle.eta(), track.cZZ());
102+
histos.fill(HIST("CovMat_cSnpY"), mcParticle.pt(), mcParticle.eta(), track.cSnpY());
103+
histos.fill(HIST("CovMat_cSnpZ"), mcParticle.pt(), mcParticle.eta(), track.cSnpZ());
104+
histos.fill(HIST("CovMat_cSnpSnp"), mcParticle.pt(), mcParticle.eta(), track.cSnpSnp());
105+
histos.fill(HIST("CovMat_cTglY"), mcParticle.pt(), mcParticle.eta(), track.cTglY());
106+
histos.fill(HIST("CovMat_cTglZ"), mcParticle.pt(), mcParticle.eta(), track.cTglZ());
107+
histos.fill(HIST("CovMat_cTglSnp"), mcParticle.pt(), mcParticle.eta(), track.cTglSnp());
108+
histos.fill(HIST("CovMat_cTglTgl"), mcParticle.pt(), mcParticle.eta(), track.cTglTgl());
109+
histos.fill(HIST("CovMat_c1PtY"), mcParticle.pt(), mcParticle.eta(), track.c1PtY());
110+
histos.fill(HIST("CovMat_c1PtZ"), mcParticle.pt(), mcParticle.eta(), track.c1PtZ());
111+
histos.fill(HIST("CovMat_c1PtSnp"), mcParticle.pt(), mcParticle.eta(), track.c1PtSnp());
112+
histos.fill(HIST("CovMat_c1PtTgl"), mcParticle.pt(), mcParticle.eta(), track.c1PtTgl());
113+
histos.fill(HIST("CovMat_c1Pt21Pt2"), mcParticle.pt(), mcParticle.eta(), track.c1Pt21Pt2());
114+
}
115+
116+
for (const auto& mcParticle : mcParticles) {
117+
if (mcParticle.pdgCode() != pdg) {
118+
continue;
119+
}
120+
if (!MC::isPhysicalPrimary(mcParticles, mcParticle)) { // Requiring is physical primary
121+
continue;
122+
}
123+
124+
if (std::find(recoTracks.begin(), recoTracks.end(), mcParticle.globalIndex()) != recoTracks.end()) {
125+
histos.fill(HIST("Efficiency"), mcParticle.pt(), mcParticle.eta(), 1.);
126+
} else {
127+
histos.fill(HIST("Efficiency"), mcParticle.pt(), mcParticle.eta(), 0.);
128+
}
129+
}
130+
}
131+
};
132+
133+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
134+
{
135+
WorkflowSpec w;
136+
if (cfgc.options().get<int>("lul-el")) {
137+
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Electron>>(cfgc, TaskName{"alice3-lutmaker-electron"}));
138+
}
139+
if (cfgc.options().get<int>("lut-mu")) {
140+
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Muon>>(cfgc, TaskName{"alice3-lutmaker-muon"}));
141+
}
142+
if (cfgc.options().get<int>("lut-pi")) {
143+
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Pion>>(cfgc, TaskName{"alice3-lutmaker-pion"}));
144+
}
145+
if (cfgc.options().get<int>("lut-ka")) {
146+
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Kaon>>(cfgc, TaskName{"alice3-lutmaker-kaon"}));
147+
}
148+
if (cfgc.options().get<int>("lut-pr")) {
149+
w.push_back(adaptAnalysisTask<Alice3LutMaker<o2::track::PID::Proton>>(cfgc, TaskName{"alice3-lutmaker-proton"}));
150+
}
151+
return w;
152+
}

‎Analysis/ALICE3/src/alice3-qa-singleparticle.cxx‎

Lines changed: 130 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,130 @@
1+
// Copyright CERN and copyright holders of ALICE O2. This software is
2+
// distributed under the terms of the GNU General Public License v3 (GPL
3+
// Version 3), copied verbatim in the file "COPYING".
4+
//
5+
// See http://alice-o2.web.cern.ch/license for full licensing information.
6+
//
7+
// In applying this license CERN does not waive the privileges and immunities
8+
// granted to it by virtue of its status as an Intergovernmental Organization
9+
// or submit itself to any jurisdiction.
10+
/// \author Nicolo' Jacazio <nicolo.jacazio@cern.ch>, CERN
11+
12+
// O2 includes
13+
#include "Framework/AnalysisTask.h"
14+
#include "Framework/runDataProcessing.h"
15+
#include "AnalysisCore/MC.h"
16+
#include "Framework/HistogramRegistry.h"
17+
18+
using namespace o2;
19+
using namespace o2::framework;
20+
using namespace o2::framework::expressions;
21+
22+
struct Alice3SingleParticle {
23+
Configurable<int> PDG{"PDG", 2212, "PDG code of the particle of interest"};
24+
Configurable<int> IsStable{"IsStable", 0, "Flag to check stable particles"};
25+
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
26+
Configurable<int> ptBins{"pt-bins", 500, "Number of pT bins"};
27+
Configurable<float> ptMin{"pt-min", 0.f, "Lower limit in pT"};
28+
Configurable<float> ptMax{"pt-max", 5.f, "Upper limit in pT"};
29+
Configurable<int> etaBins{"eta-bins", 500, "Number of eta bins"};
30+
Configurable<float> etaMin{"eta-min", -3.f, "Lower limit in eta"};
31+
Configurable<float> etaMax{"eta-max", 3.f, "Upper limit in eta"};
32+
Configurable<float> yMin{"y-min", -3.f, "Lower limit in y"};
33+
Configurable<float> yMax{"y-max", 3.f, "Upper limit in y"};
34+
Configurable<int> prodBins{"prod-bins", 100, "Number of production vertex bins"};
35+
Configurable<float> prodMin{"prod-min", -1.f, "Lower limit in production vertex"};
36+
Configurable<float> prodMax{"prod-max", 1.f, "Upper limit in production vertex"};
37+
Configurable<float> charge{"charge", 1.f, "Particle charge to scale the reconstructed momentum"};
38+
Configurable<bool> doPrint{"doPrint", false, "Flag to print debug messages"};
39+
40+
void init(InitContext&)
41+
{
42+
const TString tit = Form("%i", PDG.value);
43+
AxisSpec axisPt{ptBins, ptMin, ptMax};
44+
AxisSpec axisEta{etaBins, etaMin, etaMax};
45+
AxisSpec axisProd{prodBins, prodMin, prodMax};
46+
47+
histos.add("particlePt", "Particle Pt " + tit + ";#it{p}_{T} (GeV/#it{c})", kTH1D, {axisPt});
48+
histos.add("prodVx", "Particle Prod. Vertex X " + tit + ";Prod. Vertex X (cm)", kTH1D, {axisProd});
49+
histos.add("prodVy", "Particle Prod. Vertex Y " + tit + ";Prod. Vertex Y (cm)", kTH1D, {axisProd});
50+
histos.add("prodVz", "Particle Prod. Vertex Z " + tit + ";Prod. Vertex Z (cm)", kTH1D, {axisProd});
51+
histos.add("prodRadius", "Particle Prod. Vertex Radius " + tit + ";Prod. Vertex Radius (cm)", kTH1D, {axisProd});
52+
histos.add("prodVxVsPt", "Particle Prod. Vertex X " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex X (cm)", kTH2D, {axisPt, axisProd});
53+
histos.add("prodVyVsPt", "Particle Prod. Vertex Y " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Y (cm)", kTH2D, {axisPt, axisProd});
54+
histos.add("prodVzVsPt", "Particle Prod. Vertex Z " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Z (cm)", kTH2D, {axisPt, axisProd});
55+
histos.add("prodRadiusVsPt", "Particle Prod. Vertex Radius " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Radius (cm)", kTH2D, {axisPt, axisProd});
56+
histos.add("prodRadius3DVsPt", "Particle Prod. Vertex Radius XYZ " + tit + ";#it{p}_{T} (GeV/#it{c});Prod. Vertex Radius XYZ (cm)", kTH2D, {axisPt, axisProd});
57+
histos.add("trackPt", "Track Pt " + tit + ";#it{p}_{T} (GeV/#it{c})", kTH1D, {axisPt});
58+
histos.add("particleEta", "Particle Eta " + tit + ";#it{#eta}", kTH1D, {axisEta});
59+
histos.add("trackEta", "Track Eta " + tit + ";#it{#eta}", kTH1D, {axisEta});
60+
histos.add("particleY", "Particle Y " + tit + ";#it{y}", kTH1D, {axisEta});
61+
histos.add("primaries", "Source for primaries " + tit + ";PDG Code", kTH1D, {{100, 0.f, 100.f}});
62+
histos.add("secondaries", "Source for secondaries " + tit + ";PDG Code", kTH1D, {{100, 0.f, 100.f}});
63+
}
64+
65+
void process(const soa::Join<o2::aod::Tracks, o2::aod::McTrackLabels>& tracks,
66+
const aod::McParticles& mcParticles)
67+
{
68+
69+
std::vector<int64_t> ParticlesOfInterest;
70+
for (const auto& mcParticle : mcParticles) {
71+
if (mcParticle.pdgCode() != PDG) {
72+
continue;
73+
}
74+
if(mcParticle.y() < yMin || mcParticle.y() > yMax) {
75+
continue;
76+
}
77+
histos.fill(HIST("particlePt"), mcParticle.pt());
78+
histos.fill(HIST("particleEta"), mcParticle.eta());
79+
histos.fill(HIST("particleY"), mcParticle.y());
80+
histos.fill(HIST("prodVx"), mcParticle.vx());
81+
histos.fill(HIST("prodVy"), mcParticle.vy());
82+
histos.fill(HIST("prodVz"), mcParticle.vz());
83+
histos.fill(HIST("prodRadius"), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy()));
84+
histos.fill(HIST("prodVxVsPt"), mcParticle.pt(), mcParticle.vx());
85+
histos.fill(HIST("prodVyVsPt"), mcParticle.pt(), mcParticle.vy());
86+
histos.fill(HIST("prodVzVsPt"), mcParticle.pt(), mcParticle.vz());
87+
histos.fill(HIST("prodRadiusVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy()));
88+
histos.fill(HIST("prodRadius3DVsPt"), mcParticle.pt(), std::sqrt(mcParticle.vx() * mcParticle.vx() + mcParticle.vy() * mcParticle.vy() + mcParticle.vz() * mcParticle.vz() ));
89+
ParticlesOfInterest.push_back(mcParticle.globalIndex());
90+
}
91+
92+
for (const auto& track : tracks) {
93+
const auto mcParticle = track.mcParticle();
94+
if (!IsStable) {
95+
if (mcParticle.mother0() < 0) {
96+
continue;
97+
}
98+
auto mother = mcParticles.iteratorAt(mcParticle.mother0());
99+
const auto ParticleIsInteresting = std::find(ParticlesOfInterest.begin(), ParticlesOfInterest.end(), mother.globalIndex()) != ParticlesOfInterest.end();
100+
if (ParticleIsInteresting && doPrint) {
101+
Printf("Track %li comes from a %i and is a %i", track.globalIndex(), mother.pdgCode(), mcParticle.pdgCode());
102+
}
103+
} else {
104+
if (mcParticle.pdgCode() != PDG) {
105+
continue;
106+
}
107+
histos.fill(HIST("trackPt"), track.pt() * charge);
108+
histos.fill(HIST("trackEta"), track.eta());
109+
if (mcParticle.mother0() < 0 && doPrint) {
110+
Printf("Track %li is a %i", track.globalIndex(), mcParticle.pdgCode());
111+
} else {
112+
auto mother = mcParticles.iteratorAt(mcParticle.mother0());
113+
if (MC::isPhysicalPrimary(mcParticles, mcParticle)) {
114+
histos.get<TH1>(HIST("primaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
115+
} else {
116+
histos.get<TH1>(HIST("secondaries"))->Fill(Form("%i", mother.pdgCode()), 1.f);
117+
}
118+
if (doPrint) {
119+
Printf("Track %li is a %i and comes from a %i and %s a primary", track.globalIndex(), mcParticle.pdgCode(), mother.pdgCode(), MC::isPhysicalPrimary(mcParticles, mcParticle) ? "is" : "is not");
120+
}
121+
}
122+
}
123+
}
124+
}
125+
};
126+
127+
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
128+
{
129+
return WorkflowSpec{adaptAnalysisTask<Alice3SingleParticle>(cfgc)};
130+
}

0 commit comments

Comments
 (0)