Skip to content

Commit 770a87f

Browse files
njacazionjacazio
authored
Update FRICH QA code (#89)
- Fix slicing of MC particles - Add PID for cdeuteron - Update A3 PID TOF task with latest O2 Data Model
1 parent a2599a1 commit 770a87f

4 files changed

Lines changed: 63 additions & 79 deletions

File tree

‎ALICE3/Core/TOFResoALICE3.h‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -71,13 +71,13 @@ float TOFResoALICE3Param(const float& momentum, const float& momentumError, cons
7171
return sqrt(etexp * etexp + parameters[0] * parameters[0] + evtimereso * evtimereso);
7272
}
7373

74-
template <o2::track::PID::ID id, typename C, typename T>
75-
float TOFResoALICE3ParamTrack(const C& collision, const T& track, const Parameters& parameters)
74+
template <o2::track::PID::ID id, typename T>
75+
float TOFResoALICE3ParamTrack(const T& track, const Parameters& parameters)
7676
{
7777
const float BETA = tan(0.25f * static_cast<float>(M_PI) - 0.5f * atan(track.tgl()));
7878
const float sigmaP = sqrt(track.pt() * track.pt() * track.sigma1Pt() * track.sigma1Pt() + (BETA * BETA - 1.f) / (BETA * (BETA * BETA + 1.f)) * (track.tgl() / sqrt(track.tgl() * track.tgl() + 1.f) - 1.f) * track.sigmaTgl() * track.sigmaTgl());
7979
// const float sigmaP = std::sqrt( track.getSigma1Pt2() ) * track.pt();
80-
return TOFResoALICE3Param(track.p(), sigmaP, collision.collisionTimeRes() * 1000.f, track.length(), o2::track::pid_constants::sMasses2Z[id], parameters);
80+
return TOFResoALICE3Param(track.p(), sigmaP, track.collision().collisionTimeRes() * 1000.f, track.length(), o2::track::pid_constants::sMasses2Z[id], parameters);
8181
// return TOFResoALICE3Param(track.p(), track.sigma1Pt(), collision.collisionTimeRes() * 1000.f, track.length(), o2::track::pid_constants::sMasses[id], parameters);
8282
}
8383

‎ALICE3/TableProducer/alice3-pidTOF.cxx‎

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -84,12 +84,12 @@ struct ALICE3pidTOFTask {
8484
template <o2::track::PID::ID id>
8585
float sigma(Trks::iterator track)
8686
{
87-
return o2::pid::tof::TOFResoALICE3ParamTrack<id>(track.collision(), track, resoParameters);
87+
return o2::pid::tof::TOFResoALICE3ParamTrack<id>(track, resoParameters);
8888
}
8989
template <o2::track::PID::ID id>
9090
float nsigma(Trks::iterator track)
9191
{
92-
return (track.trackTime() - track.collision().collisionTime() * 1000.f - tof::ExpTimes<Trks::iterator, id>::GetExpectedSignal(track)) / sigma<id>(track);
92+
return (track.trackTime() * 1e3f - track.collision().collisionTime() * 1000.f - tof::ExpTimes<Trks::iterator, id>::GetExpectedSignal(track)) / sigma<id>(track);
9393
}
9494
void process(Coll const& collisions, Trks const& tracks)
9595
{
@@ -263,12 +263,13 @@ struct ALICE3pidTOFTaskQA {
263263
continue;
264264
}
265265

266+
const float tofSignal = t.trackTime() * 1e3f;
266267
// const float tof = t.tofSignal() - collisionTime_ps;
267-
const float tof = t.trackTime() - collisionTime_ps;
268+
const float tof = tofSignal - collisionTime_ps;
268269

269270
//
270271
// histos.fill(HIST("event/tofsignal"), t.p(), t.tofSignal());
271-
histos.fill(HIST("event/tofsignal"), t.p(), t.trackTime());
272+
histos.fill(HIST("event/tofsignal"), t.p(), tofSignal);
272273
histos.fill(HIST("event/pexp"), t.p(), t.tofExpMom());
273274
histos.fill(HIST("event/eta"), t.eta());
274275
histos.fill(HIST("event/length"), t.length());

‎ALICE3/Tasks/alice3-cdeuteron.cxx‎

Lines changed: 39 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -21,6 +21,7 @@
2121
#include "ReconstructionDataFormats/PID.h"
2222
#include "Common/Core/RecoDecay.h"
2323
#include "DetectorsVertexing/DCAFitterN.h"
24+
#include "Common/Core/PID/PIDResponse.h"
2425
#include "Common/Core/trackUtilities.h"
2526

2627
using namespace o2;
@@ -40,6 +41,10 @@ struct Alice3CDeuteron {
4041
Configurable<float> minDcaPion{"minDcaPion", -100, "Minimum DCA of the pion to the primary vertex"};
4142
Configurable<float> maxDca{"maxDca", 100, "Maximum track DCA to the primary vertex"};
4243
Configurable<float> minCpa{"minCpa", 0, "Minimum CPA"};
44+
Configurable<int> usePdg{"usePdg", 1, "Flag to use the PDG instead of the TOF/RICH PID"};
45+
Configurable<float> maxNsigmaDe{"maxNsigmaDe", 3.f, "Maximum Nsigma for deuteron"};
46+
Configurable<float> maxNsigmaKa{"maxNsigmaKa", 3.f, "Maximum Nsigma for kaon"};
47+
Configurable<float> maxNsigmaPi{"maxNsigmaPi", 3.f, "Maximum Nsigma for pion"};
4348
HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject};
4449
o2::vertexing::DCAFitterN<3> fitter;
4550

@@ -72,6 +77,9 @@ struct Alice3CDeuteron {
7277
const AxisSpec axisVtxY{100, -0.1, 0.1, "Vtx_{Y}"};
7378
const AxisSpec axisVtxZ{100, -0.1, 0.1, "Vtx_{Z}"};
7479
const AxisSpec axisCPA{4000, -1.1, 1.1, "CPA"};
80+
const AxisSpec axisNSigmaDe{1000, -10, 10, "N_{#sigma}^{TOF}(d)"};
81+
const AxisSpec axisNSigmaKa{1000, -10, 10, "N_{#sigma}^{TOF}(k)"};
82+
const AxisSpec axisNSigmaPi{1000, -10, 10, "N_{#sigma}^{TOF}(#pi)"};
7583
const TString tit = Form(" [%.6f, %.6f] R [%.6f, %.6f] DCA ",
7684
minRadius.value, maxRadius.value,
7785
minDca.value, maxDca.value);
@@ -121,6 +129,12 @@ struct Alice3CDeuteron {
121129
histos.add("event/mcvtxX", "mcvtxX", kTH1D, {axisVtxX});
122130
histos.add("event/mcvtxY", "mcvtxY", kTH1D, {axisVtxY});
123131
histos.add("event/mcvtxZ", "mcvtxZ", kTH1D, {axisVtxZ});
132+
histos.add("event/nsigmaDe", "nsigmaDe", kTH2D, {axisPt, axisNSigmaDe});
133+
histos.add("event/nsigmaKa", "nsigmaKa", kTH2D, {axisPt, axisNSigmaKa});
134+
histos.add("event/nsigmaPi", "nsigmaPi", kTH2D, {axisPt, axisNSigmaPi});
135+
histos.add("event/nsigmaDecut", "nsigmaDecut", kTH2D, {axisPt, axisNSigmaDe});
136+
histos.add("event/nsigmaKacut", "nsigmaKacut", kTH2D, {axisPt, axisNSigmaKa});
137+
histos.add("event/nsigmaPicut", "nsigmaPicut", kTH2D, {axisPt, axisNSigmaPi});
124138
histos.add("event/track1dcaxy", "track1dcaxy Deuteron", kTH1D, {axisDcaXY});
125139
histos.add("event/track1dcaz", "track1dcaz Deuteron", kTH1D, {axisDcaZ});
126140
histos.add("event/candperdeuteron", "candperdeuteron", kTH1D, {{1000, 0, 10000}});
@@ -170,11 +184,12 @@ struct Alice3CDeuteron {
170184

171185
void process(const soa::Join<o2::aod::Collisions, o2::aod::McCollisionLabels>::iterator& coll,
172186
const o2::aod::McCollisions& Mccoll,
173-
const soa::Join<o2::aod::Tracks, o2::aod::McTrackLabels, o2::aod::TracksExtra, o2::aod::TracksCov>& tracks,
187+
const soa::Join<o2::aod::Tracks, o2::aod::McTrackLabels, o2::aod::TracksExtra, o2::aod::TracksCov,
188+
aod::pidTOFFullPi, aod::pidTOFFullKa, aod::pidTOFFullDe>& tracks,
174189
const aod::McParticles& mcParticles)
175190
{
176-
177-
for (auto i : mcParticles) {
191+
const auto particlesInCollision = mcParticles.sliceBy(aod::mcparticle::mcCollisionId, coll.mcCollision().globalIndex());
192+
for (const auto& i : particlesInCollision) {
178193
if (i.pdgCode() != 12345) {
179194
continue;
180195
}
@@ -217,9 +232,15 @@ struct Alice3CDeuteron {
217232
for (const auto& track1 : tracks) {
218233
const auto index1 = track1.globalIndex();
219234
int ncand = 0;
220-
if (track1.mcParticle().pdgCode() != 1000010020) {
235+
histos.fill(HIST("event/nsigmaDe"), track1.pt(), track1.tofNSigmaDe());
236+
if (usePdg) {
237+
if (track1.mcParticle().pdgCode() != 1000010020) {
238+
continue;
239+
}
240+
} else if (abs(track1.tofNSigmaDe()) > maxNsigmaDe || track1.sign() < 0.f) {
221241
continue;
222242
}
243+
histos.fill(HIST("event/nsigmaDecut"), track1.pt(), track1.tofNSigmaDe());
223244
if (!getTrackPar(track1).propagateParamToDCA(collPos,
224245
magField * 10.f, &dca1, 100.)) {
225246
continue;
@@ -232,10 +253,15 @@ struct Alice3CDeuteron {
232253
if (index1 == index2) {
233254
continue;
234255
}
235-
if (track2.mcParticle().pdgCode() != -321) {
256+
histos.fill(HIST("event/nsigmaKa"), track2.pt(), track2.tofNSigmaKa());
257+
if (usePdg) {
258+
if (track2.mcParticle().pdgCode() != -321) {
259+
continue;
260+
}
261+
} else if (abs(track2.tofNSigmaKa()) > maxNsigmaKa || track2.sign() > 0.f) {
236262
continue;
237263
}
238-
264+
histos.fill(HIST("event/nsigmaKacut"), track2.pt(), track2.tofNSigmaKa());
239265
if (!getTrackPar(track2).propagateParamToDCA(collPos,
240266
magField * 10.f, &dca2, 100.)) {
241267
continue;
@@ -250,9 +276,15 @@ struct Alice3CDeuteron {
250276
if (index1 == index3) {
251277
continue;
252278
}
253-
if (track3.mcParticle().pdgCode() != 211) {
279+
histos.fill(HIST("event/nsigmaPi"), track3.pt(), track3.tofNSigmaPi());
280+
if (usePdg) {
281+
if (track3.mcParticle().pdgCode() != 211) {
282+
continue;
283+
}
284+
} else if (abs(track3.tofNSigmaPi()) > maxNsigmaPi || track3.sign() < 0.f) {
254285
continue;
255286
}
287+
histos.fill(HIST("event/nsigmaPicut"), track3.pt(), track3.tofNSigmaPi());
256288
bool iscut = false;
257289
if (abs(dca1[0]) < minDca || abs(dca1[1]) < minDca) {
258290
iscut = true;

‎ALICE3/Tasks/pidRICHqa.cxx‎

Lines changed: 16 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -44,25 +44,11 @@ DECLARE_SOA_INDEX_COLUMN(FRICH, frich);
4444

4545
DECLARE_SOA_INDEX_TABLE_USER(RICHTracksIndex, Tracks, "RICHTRK", indices::TrackId, indices::RICHId);
4646
DECLARE_SOA_INDEX_TABLE_USER(FRICHTracksIndex, Tracks, "FRICHTRK", indices::TrackId, indices::FRICHId);
47-
DECLARE_SOA_INDEX_TABLE_USER(BothRICHTracksIndex, Tracks, "FRICHTRK", indices::TrackId, indices::RICHId, indices::FRICHId);
4847
} // namespace o2::aod
4948

5049
struct richIndexBuilder { // Builder of the RICH-track index linkage
51-
Builds<o2::aod::RICHTracksIndex> ind;
52-
void init(o2::framework::InitContext&)
53-
{
54-
}
55-
};
56-
57-
struct frichIndexBuilder { // Builder of the FRICH-track index linkage
58-
Builds<o2::aod::FRICHTracksIndex> ind;
59-
void init(o2::framework::InitContext&)
60-
{
61-
}
62-
};
63-
64-
struct bothrichIndexBuilder { // Builder of the RICH and FRICH track index linkage
65-
Builds<o2::aod::BothRICHTracksIndex> ind;
50+
Builds<o2::aod::RICHTracksIndex> indB;
51+
Builds<o2::aod::FRICHTracksIndex> indF;
6652
void init(o2::framework::InitContext&)
6753
{
6854
}
@@ -91,27 +77,6 @@ struct richPidQaMc {
9177
Configurable<float> maxDelta{"maxDelta", 0.4f, "Maximum delta plotted (rad)"};
9278
Configurable<int> logAxis{"logAxis", 1, "Flag to use a log momentum axis"};
9379

94-
template <typename T>
95-
void makelogaxis(T h)
96-
{
97-
if (logAxis == 0) {
98-
return;
99-
}
100-
const int nbins = h->GetNbinsX();
101-
double binp[nbins + 1];
102-
double max = h->GetXaxis()->GetBinUpEdge(nbins);
103-
double min = h->GetXaxis()->GetBinLowEdge(1);
104-
if (min <= 0) {
105-
min = 0.00001;
106-
}
107-
double lmin = TMath::Log10(min);
108-
double ldelta = (TMath::Log10(max) - lmin) / ((double)nbins);
109-
for (int i = 0; i < nbins; i++) {
110-
binp[i] = TMath::Exp(TMath::Log(10) * (lmin + i * ldelta));
111-
}
112-
binp[nbins] = max + 1;
113-
h->GetXaxis()->Set(nbins, binp);
114-
}
11580
static constexpr int Np = 5;
11681
static constexpr std::string_view hdelta[Np] = {"delta/El", "delta/Mu", "delta/Pi", "delta/Ka", "delta/Pr"};
11782
static constexpr std::string_view hnsigma[Np] = {"nsigma/El", "nsigma/Mu", "nsigma/Pi", "nsigma/Ka", "nsigma/Pr"};
@@ -130,7 +95,10 @@ struct richPidQaMc {
13095
template <uint8_t i>
13196
void addParticleHistos()
13297
{
133-
const AxisSpec ptAxis{nBinsP, minP, maxP, "#it{p}_{T} (GeV/#it{c})"};
98+
AxisSpec ptAxis{nBinsP, minP, maxP, "#it{p}_{T} (GeV/#it{c})"};
99+
if (logAxis) {
100+
ptAxis.makeLogaritmic();
101+
}
134102
const AxisSpec nsigmaAxis{nBinsNsigma, minNsigma, maxNsigma, Form("N_{#sigma}^{RICH}(%s)", pT[pid_type])};
135103

136104
TString tit = Form("%s", pT[i]);
@@ -139,16 +107,18 @@ struct richPidQaMc {
139107
}
140108
// NSigma
141109
histos.add(hnsigmaMC[i].data(), "True " + tit, HistType::kTH2F, {ptAxis, nsigmaAxis});
142-
makelogaxis(histos.get<TH2>(HIST(hnsigmaMC[i])));
143110
histos.add(hnsigmaMCprm[i].data(), "True Primary " + tit, HistType::kTH2F, {ptAxis, nsigmaAxis});
144-
makelogaxis(histos.get<TH2>(HIST(hnsigmaMCprm[i])));
145111
histos.add(hnsigmaMCsec[i].data(), "True Secondary " + tit, HistType::kTH2F, {ptAxis, nsigmaAxis});
146-
makelogaxis(histos.get<TH2>(HIST(hnsigmaMCsec[i])));
147112
}
148113
void init(o2::framework::InitContext&)
149114
{
150-
const AxisSpec momAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"};
151-
const AxisSpec ptAxis{nBinsP, minP, maxP, "#it{p}_{T} (GeV/#it{c})"};
115+
AxisSpec momAxis{nBinsP, minP, maxP, "#it{p} (GeV/#it{c})"};
116+
AxisSpec ptAxis{nBinsP, minP, maxP, "#it{p}_{T} (GeV/#it{c})"};
117+
if (logAxis) {
118+
momAxis.makeLogaritmic();
119+
ptAxis.makeLogaritmic();
120+
}
121+
152122
const AxisSpec sigAxis{1000, 0, 0.3, "Cherenkov angle (rad)"};
153123
const AxisSpec nsigmaAxis{nBinsNsigma, minNsigma, maxNsigma, Form("N_{#sigma}^{RICH}(%s)", pT[pid_type])};
154124
const AxisSpec deltaAxis{nBinsDelta, minDelta, maxDelta, Form("#Delta(%s) (rad)", pT[pid_type])};
@@ -162,21 +132,15 @@ struct richPidQaMc {
162132
histos.add("qa/eta", ";#it{#eta}", kTH1F, {{100, -4, 4}});
163133
histos.add("qa/signalerror", ";Cherenkov angle error (rad)", kTH1F, {{100, 0, 1}});
164134
histos.add("qa/signalvsP", "", kTH2F, {momAxis, sigAxis});
165-
makelogaxis(histos.get<TH2>(HIST("qa/signalvsP")));
166135
histos.add("qa/signalvsPPrim", "", kTH2F, {momAxis, sigAxis});
167-
makelogaxis(histos.get<TH2>(HIST("qa/signalvsPPrim")));
168136
histos.add("qa/signalvsPSec", "", kTH2F, {momAxis, sigAxis});
169-
makelogaxis(histos.get<TH2>(HIST("qa/signalvsPSec")));
170137
histos.add(hdelta[pid_type].data(), "", kTH2F, {momAxis, deltaAxis});
171-
makelogaxis(histos.get<TH2>(HIST(hdelta[pid_type])));
172138
histos.add(hnsigma[pid_type].data(), "", HistType::kTH2F, {ptAxis, nsigmaAxis});
173-
makelogaxis(histos.get<TH2>(HIST(hnsigma[pid_type])));
174139
histos.add(hnsigmaprm[pid_type].data(), "Primary", HistType::kTH2F, {ptAxis, nsigmaAxis});
175-
makelogaxis(histos.get<TH2>(HIST(hnsigmaprm[pid_type])));
176140
histos.add(hnsigmasec[pid_type].data(), "Secondary", HistType::kTH2F, {ptAxis, nsigmaAxis});
177-
makelogaxis(histos.get<TH2>(HIST(hnsigmasec[pid_type])));
141+
178142
histos.add(hfrichnsigma[pid_type].data(), "", HistType::kTH2F, {ptAxis, nsigmaAxis});
179-
makelogaxis(histos.get<TH2>(HIST(hfrichnsigma[pid_type])));
143+
180144
addParticleHistos<0>();
181145
addParticleHistos<1>();
182146
addParticleHistos<2>();
@@ -202,11 +166,9 @@ struct richPidQaMc {
202166
aod::pidTOFFullEl, aod::pidTOFFullMu, aod::pidTOFFullPi,
203167
aod::pidTOFFullKa, aod::pidTOFFullPr>;
204168
using TrksfRICH = soa::Join<aod::Tracks, aod::FRICHTracksIndex, aod::TracksExtra>;
205-
using TrksbothRICH = soa::Join<aod::Tracks, aod::BothRICHTracksIndex, aod::TracksExtra>;
206169
void process(const aod::McParticles& mcParticles,
207170
const Trks& tracks,
208171
const TrksfRICH& tracksfrich,
209-
const TrksbothRICH& tracksbothrich,
210172
const aod::McTrackLabels& labels,
211173
const aod::RICHs&,
212174
const aod::FRICHs&,
@@ -337,23 +299,12 @@ struct richPidQaMc {
337299
histos.fill(HIST(hfrichnsigmasec[pid_type]), track.pt(), nsigma);
338300
}
339301
}
340-
341-
for (const auto& track : tracksbothrich) { // Barrel and Forward RICH
342-
if (!track.has_frich()) {
343-
continue;
344-
}
345-
if (!track.has_rich()) {
346-
continue;
347-
}
348-
}
349302
}
350303
};
351304

352305
WorkflowSpec defineDataProcessing(ConfigContext const& cfg)
353306
{
354-
auto workflow = WorkflowSpec{adaptAnalysisTask<richIndexBuilder>(cfg),
355-
adaptAnalysisTask<frichIndexBuilder>(cfg),
356-
adaptAnalysisTask<bothrichIndexBuilder>(cfg)};
307+
auto workflow = WorkflowSpec{adaptAnalysisTask<richIndexBuilder>(cfg)};
357308
if (cfg.options().get<int>("qa-el")) {
358309
workflow.push_back(adaptAnalysisTask<richPidQaMc<PID::Electron>>(cfg, TaskName{"pidRICH-qa-El"}));
359310
}

0 commit comments

Comments
 (0)