Skip to content

Commit 9fedd46

Browse files
njacazionjacazio
authored
Extend QA efficiency task (#318)
- Add track cuts to efficiency - Improve CPU efficiency
1 parent 1aec497 commit 9fedd46

1 file changed

Lines changed: 76 additions & 47 deletions

File tree

‎DPG/Tasks/qaEfficiency.cxx‎

Lines changed: 76 additions & 47 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,15 @@
1919
#include "ReconstructionDataFormats/DCA.h"
2020
#include "ReconstructionDataFormats/Track.h"
2121
#include "Common/Core/MC.h"
22+
#include "Common/Core/TrackSelection.h"
23+
#include "Common/Core/TrackSelectionDefaults.h"
2224
#include "Common/DataModel/TrackSelectionTables.h"
2325

26+
// ROOT includes
27+
#include "TPDGCode.h"
28+
#include "TEfficiency.h"
29+
#include "TList.h"
30+
2431
using namespace o2::framework;
2532

2633
void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
@@ -40,11 +47,6 @@ void customize(std::vector<o2::framework::ConfigParamSpec>& workflowOptions)
4047

4148
#include "Framework/runDataProcessing.h"
4249

43-
// ROOT includes
44-
#include "TPDGCode.h"
45-
#include "TEfficiency.h"
46-
#include "TList.h"
47-
4850
/// Task to QA the efficiency of a particular particle defined by its pdg code
4951
template <o2::track::pid_constants::ID particle>
5052
struct QaTrackingEfficiency {
@@ -82,28 +84,15 @@ struct QaTrackingEfficiency {
8284
if (pdgSign != 0 && pdgSign != 1 && pdgSign != -1) {
8385
LOG(fatal) << "Provide pdgSign as 0, 1, -1. Provided: " << pdgSign.value;
8486
}
85-
const TString tagPt = Form("%s #it{#eta} [%.2f,%.2f] #it{#varphi} [%.2f,%.2f] Prim %i",
86-
o2::track::pid_constants::sNames[particle],
87-
etaMin.value, etaMax.value,
88-
phiMin.value, phiMax.value,
89-
selPrim.value);
9087
AxisSpec axisPt{ptBins, ptMin, ptMax, "#it{p}_{T} (GeV/#it{c})"};
9188
if (logPt) {
9289
axisPt.makeLogaritmic();
9390
}
94-
95-
const TString tagEta = Form("%s #it{p}_{T} [%.2f,%.2f] #it{#varphi} [%.2f,%.2f] Prim %i",
96-
o2::track::pid_constants::sNames[particle],
97-
ptMin.value, ptMax.value,
98-
phiMin.value, phiMax.value,
99-
selPrim.value);
91+
AxisSpec axisP{ptBins, ptMin, ptMax, "#it{p} (GeV/#it{c})"};
92+
if (logPt) {
93+
axisP.makeLogaritmic();
94+
}
10095
const AxisSpec axisEta{etaBins, etaMin, etaMax, "#it{#eta}"};
101-
102-
const TString tagPhi = Form("%s #it{#eta} [%.2f,%.2f] #it{p}_{T} [%.2f,%.2f] Prim %i",
103-
o2::track::pid_constants::sNames[particle],
104-
etaMin.value, etaMax.value,
105-
ptMin.value, ptMax.value,
106-
selPrim.value);
10796
const AxisSpec axisPhi{phiBins, phiMin, phiMax, "#it{#varphi} (rad)"};
10897

10998
const AxisSpec axisSel{9, 0.5, 9.5, "Selection"};
@@ -134,40 +123,65 @@ struct QaTrackingEfficiency {
134123
histos.add("eventMultiplicity", "Event Selection", kTH1D, {{1000, 0, 5000}});
135124
histos.add("trackLength", "Track length;Track length (cm)", kTH1D, {{2000, -1000, 1000}});
136125

126+
const TString tagPt = Form("%s #it{#eta} [%.2f,%.2f] #it{#varphi} [%.2f,%.2f] Prim %i",
127+
o2::track::pid_constants::sNames[particle],
128+
etaMin.value, etaMax.value,
129+
phiMin.value, phiMax.value,
130+
selPrim.value);
137131
histos.add("pt/num", "Numerator " + tagPt, kTH1D, {axisPt});
138132
histos.add("pt/den", "Denominator " + tagPt, kTH1D, {axisPt});
139133

134+
histos.add("p/num", "Numerator " + tagPt, kTH1D, {axisP});
135+
histos.add("p/den", "Denominator " + tagPt, kTH1D, {axisP});
136+
137+
const TString tagEta = Form("%s #it{p}_{T} [%.2f,%.2f] #it{#varphi} [%.2f,%.2f] Prim %i",
138+
o2::track::pid_constants::sNames[particle],
139+
ptMin.value, ptMax.value,
140+
phiMin.value, phiMax.value,
141+
selPrim.value);
140142
histos.add("eta/num", "Numerator " + tagEta, kTH1D, {axisEta});
141143
histos.add("eta/den", "Denominator " + tagEta, kTH1D, {axisEta});
142144

145+
const TString tagPhi = Form("%s #it{#eta} [%.2f,%.2f] #it{p}_{T} [%.2f,%.2f] Prim %i",
146+
o2::track::pid_constants::sNames[particle],
147+
etaMin.value, etaMax.value,
148+
ptMin.value, ptMax.value,
149+
selPrim.value);
143150
histos.add("phi/num", "Numerator " + tagPhi, kTH1D, {axisPhi});
144151
histos.add("phi/den", "Denominator " + tagPhi, kTH1D, {axisPhi});
145152

153+
const TString tagPtEta = Form("%s #it{#varphi} [%.2f,%.2f] Prim %i",
154+
o2::track::pid_constants::sNames[particle],
155+
phiMin.value, phiMax.value,
156+
selPrim.value);
157+
histos.add("pteta/num", "Numerator " + tagPtEta, kTH2D, {axisPt, axisEta});
158+
histos.add("pteta/den", "Denominator " + tagPtEta, kTH2D, {axisPt, axisEta});
159+
146160
list.setObject(new TList);
147161
if (makeEff) {
148162
auto makeEfficiency = [&](TString effname, TString efftitle, auto templateHisto) {
149-
TAxis* axis = histos.get<TH1>(templateHisto)->GetXaxis();
163+
const TAxis* axis = histos.get<TH1>(templateHisto)->GetXaxis();
150164
if (axis->IsVariableBinSize()) {
151165
list->Add(new TEfficiency(effname, efftitle, axis->GetNbins(), axis->GetXbins()->GetArray()));
152166
} else {
153167
list->Add(new TEfficiency(effname, efftitle, axis->GetNbins(), axis->GetXmin(), axis->GetXmax()));
154168
}
155169
};
156-
auto makeEfficiency2D = [&](TString effname, TString efftitle, auto templateHistoX, auto templateHistoY) {
157-
TAxis* axisX = histos.get<TH1>(templateHistoX)->GetXaxis();
158-
TAxis* axisY = histos.get<TH1>(templateHistoY)->GetXaxis();
170+
makeEfficiency("efficiencyVsPt", "Efficiency " + tagPt + ";#it{p}_{T} (GeV/#it{c});Efficiency", HIST("pt/num"));
171+
makeEfficiency("efficiencyVsP", "Efficiency " + tagPt + ";#it{p} (GeV/#it{c});Efficiency", HIST("pt/num"));
172+
makeEfficiency("efficiencyVsEta", "Efficiency " + tagEta + ";#it{#eta};Efficiency", HIST("eta/num"));
173+
makeEfficiency("efficiencyVsPhi", "Efficiency " + tagPhi + ";#it{#varphi} (rad);Efficiency", HIST("phi/num"));
174+
175+
auto makeEfficiency2D = [&](TString effname, TString efftitle, auto templateHisto) {
176+
const TAxis* axisX = histos.get<TH2>(templateHisto)->GetXaxis();
177+
const TAxis* axisY = histos.get<TH2>(templateHisto)->GetXaxis();
159178
if (axisX->IsVariableBinSize() || axisY->IsVariableBinSize()) {
160179
list->Add(new TEfficiency(effname, efftitle, axisX->GetNbins(), axisX->GetXbins()->GetArray(), axisY->GetNbins(), axisY->GetXbins()->GetArray()));
161180
} else {
162181
list->Add(new TEfficiency(effname, efftitle, axisX->GetNbins(), axisX->GetXmin(), axisX->GetXmax(), axisY->GetNbins(), axisY->GetXmin(), axisY->GetXmax()));
163182
}
164183
};
165-
makeEfficiency("efficiencyVsPt", "Efficiency " + tagPt + ";#it{p}_{T} (GeV/#it{c});Efficiency", HIST("pt/num"));
166-
makeEfficiency("efficiencyVsP", "Efficiency " + tagPt + ";#it{p} (GeV/#it{c});Efficiency", HIST("pt/num"));
167-
makeEfficiency("efficiencyVsEta", "Efficiency " + tagEta + ";#it{#eta};Efficiency", HIST("eta/num"));
168-
makeEfficiency("efficiencyVsPhi", "Efficiency " + tagPhi + ";#it{#varphi} (rad);Efficiency", HIST("phi/num"));
169-
170-
makeEfficiency2D("efficiencyVsPtVsEta", Form("Efficiency %s #it{#varphi} [%.2f,%.2f] Prim %i;%s;%s;Efficiency", o2::track::pid_constants::sNames[particle], phiMin.value, phiMax.value, selPrim.value, "#it{p}_{T} (GeV/#it{c})", "#it{#eta}"), HIST("pt/num"), HIST("eta/num"));
184+
makeEfficiency2D("efficiencyVsPtVsEta", "Efficiency " + tagPtEta + ";#it{#varphi} (rad);Efficiency", HIST("pteta/num"));
171185
}
172186
}
173187

@@ -245,8 +259,6 @@ struct QaTrackingEfficiency {
245259
return false;
246260
};
247261

248-
std::vector<int64_t> recoTracks(tracks.size());
249-
int ntrks = 0;
250262
for (const auto& track : tracks) {
251263
const auto mcParticle = track.mcParticle();
252264
if (rejectParticle(mcParticle, HIST("trackSelection"))) {
@@ -268,10 +280,11 @@ struct QaTrackingEfficiency {
268280

269281
histos.fill(HIST("trackSelection"), 8);
270282
histos.fill(HIST("trackLength"), track.length());
283+
histos.fill(HIST("p/num"), mcParticle.p());
271284
histos.fill(HIST("pt/num"), mcParticle.pt());
272285
histos.fill(HIST("eta/num"), mcParticle.eta());
273286
histos.fill(HIST("phi/num"), mcParticle.phi());
274-
recoTracks[ntrks++] = mcParticle.globalIndex();
287+
histos.fill(HIST("pteta/num"), mcParticle.pt(), mcParticle.eta());
275288
}
276289

277290
float dNdEta = 0;
@@ -283,19 +296,30 @@ struct QaTrackingEfficiency {
283296
continue;
284297
}
285298

286-
if (makeEff) {
287-
const auto particleReconstructed = std::find(recoTracks.begin(), recoTracks.end(), mcParticle.globalIndex()) != recoTracks.end();
288-
static_cast<TEfficiency*>(list->At(0))->Fill(particleReconstructed, mcParticle.pt());
289-
static_cast<TEfficiency*>(list->At(1))->Fill(particleReconstructed, mcParticle.p());
290-
static_cast<TEfficiency*>(list->At(2))->Fill(particleReconstructed, mcParticle.eta());
291-
static_cast<TEfficiency*>(list->At(3))->Fill(particleReconstructed, mcParticle.phi());
292-
static_cast<TEfficiency*>(list->At(4))->Fill(particleReconstructed, mcParticle.pt(), mcParticle.eta());
293-
}
299+
histos.fill(HIST("p/den"), mcParticle.p());
294300
histos.fill(HIST("pt/den"), mcParticle.pt());
295301
histos.fill(HIST("eta/den"), mcParticle.eta());
296302
histos.fill(HIST("phi/den"), mcParticle.phi());
303+
histos.fill(HIST("pteta/den"), mcParticle.pt(), mcParticle.eta());
297304
}
298305
histos.fill(HIST("eventMultiplicity"), dNdEta * 0.5f / 2.f);
306+
307+
if (makeEff) {
308+
auto fillEfficiency = [&](int index, auto num, auto den) {
309+
static_cast<TEfficiency*>(list->At(index))->SetTotalHistogram(*histos.get<TH1>(den).get(), "f");
310+
static_cast<TEfficiency*>(list->At(index))->SetPassedHistogram(*histos.get<TH1>(num).get(), "f");
311+
};
312+
fillEfficiency(0, HIST("pt/num"), HIST("pt/den"));
313+
fillEfficiency(1, HIST("p/num"), HIST("p/den"));
314+
fillEfficiency(2, HIST("eta/num"), HIST("eta/den"));
315+
fillEfficiency(3, HIST("phi/num"), HIST("phi/den"));
316+
317+
auto fillEfficiency2D = [&](int index, auto num, auto den) {
318+
static_cast<TEfficiency*>(list->At(index))->SetTotalHistogram(*histos.get<TH2>(den).get(), "f");
319+
static_cast<TEfficiency*>(list->At(index))->SetPassedHistogram(*histos.get<TH2>(num).get(), "f");
320+
};
321+
fillEfficiency2D(4, HIST("pteta/num"), HIST("pteta/den"));
322+
}
299323
}
300324
};
301325

@@ -355,6 +379,7 @@ struct QaTrackingEfficiencyData {
355379
histos.get<TH1>(HIST("trackSelection"))->GetXaxis()->SetBinLabel(2, "Passed #it{p}_{T}");
356380
histos.get<TH1>(HIST("trackSelection"))->GetXaxis()->SetBinLabel(3, "Passed #it{#eta}");
357381
histos.get<TH1>(HIST("trackSelection"))->GetXaxis()->SetBinLabel(4, "Passed #it{#varphi}");
382+
histos.get<TH1>(HIST("trackSelection"))->GetXaxis()->SetBinLabel(5, "Passed quality cuts");
358383

359384
histos.add("trackLength", "Track length;Track length (cm)", kTH1D, {{2000, -1000, 1000}});
360385

@@ -397,7 +422,7 @@ struct QaTrackingEfficiencyData {
397422
}
398423

399424
void process(const o2::aod::Collision& collision,
400-
const o2::soa::Join<o2::aod::Tracks, o2::aod::TracksExtra>& tracks)
425+
const o2::soa::Join<o2::aod::Tracks, o2::aod::TracksExtra, o2::aod::TrackSelection>& tracks)
401426
{
402427

403428
histos.fill(HIST("eventSelection"), 1);
@@ -413,17 +438,21 @@ struct QaTrackingEfficiencyData {
413438
for (const auto& track : tracks) {
414439
histos.fill(HIST("trackSelection"), 1);
415440
if ((track.pt() < ptMin || track.pt() > ptMax)) { // Check pt
416-
return;
441+
continue;
417442
}
418443
histos.fill(HIST("trackSelection"), 2);
419444
if ((track.eta() < etaMin || track.eta() > etaMax)) { // Check eta
420-
return;
445+
continue;
421446
}
422447
histos.fill(HIST("trackSelection"), 3);
423448
if ((track.phi() < phiMin || track.phi() > phiMax)) { // Check phi
424-
return;
449+
continue;
425450
}
426451
histos.fill(HIST("trackSelection"), 4);
452+
if (!track.isGlobalTrack()) { // Check general cuts
453+
continue;
454+
}
455+
histos.fill(HIST("trackSelection"), 5);
427456

428457
histos.fill(HIST("trackLength"), track.length());
429458
histos.fill(HIST("pt/den"), track.pt());

0 commit comments

Comments
 (0)