ATLAS Offline Software
Loading...
Searching...
No Matches
MuonPerformanceAlg.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2024 CERN for the benefit of the ATLAS collaboration
3*/
4
5#include "MuonPerformanceAlg.h"
6
7#include <vector>
8
9#include "AthenaBaseComps/AthCheckMacros.h"
10#include "GaudiKernel/IToolSvc.h"
11#include "xAODMuon/Muon.h"
12#include "xAODMuon/MuonSegment.h"
13#include "xAODMuon/MuonSegmentContainer.h"
14#include "xAODTruth/TruthParticle.h"
15#include "xAODTruth/TruthParticleAuxContainer.h"
16#include "xAODTruth/TruthParticleContainer.h"
17#include "TruthUtils/MagicNumbers.h"
18#include "AthContainers/ConstAccessor.h"
19
20MuonPerformanceAlg::MuonPerformanceAlg(const std::string& name, ISvcLocator* pSvcLocator) :
21 AthAlgorithm(name, pSvcLocator), m_writeToFile(false), m_nevents(0), m_runNumber(0), m_eventNumber(0) {
22 declareProperty("writeToFile", m_writeToFile = false);
23 declareProperty("FileName", m_fileName = "MuonPerformanceAlg.txt");
24 declareProperty("ConsideredPDGs", m_pdgsToBeConsidered);
25}
26
27StatusCode MuonPerformanceAlg::initialize() {
28 // initialize cuts, please make sure the number of bins and the sizes of the cuts + string are always the same
29 unsigned int nbins = 12;
30 m_nevents = 0;
31 m_ntruth.resize(nbins);
32 m_ntruth5.resize(nbins);
33 m_ntruth10.resize(nbins);
34 m_nfound.resize(nbins);
35 m_nfound5.resize(nbins);
36 m_nfound10.resize(nbins);
37 m_nfoundr.resize(nbins);
38 m_nfoundr5.resize(nbins);
39 m_nfoundr10.resize(nbins);
40 m_nreco.resize(nbins);
41 m_nreco5.resize(nbins);
42 m_nreco10.resize(nbins);
43 m_hitCutString = {"SA 2.0 ", "CB all ", "MuidCB ", "MuGirl ", "Tag ", "Calo ",
44 "ID ", "SA(no ID)", "Tight ", "Medium ", "Loose ", "Combined "};
45
46 // add muons
47 if (m_pdgsToBeConsidered.value().empty()) {
48 m_selectedPdgs.insert(13);
49 m_selectedPdgs.insert(-13);
50 } else {
51 // add pdgs
52 for (auto pdg : m_pdgsToBeConsidered.value()) m_selectedPdgs.insert(pdg);
53 }
54 if (!m_muonsNameKey.key().empty()) ATH_CHECK(m_muonsNameKey.initialize());
55
56 ATH_CHECK(m_eventInfo.initialize());
57 ATH_CHECK(m_truthMuons.initialize());
58 ATH_CHECK(m_muonTruthParticleKey.initialize());
59 ATH_CHECK(m_truthMuonRecoMuonLinkKey.initialize());
60 ATH_CHECK(m_truthMuonTruthType.initialize());
61 ATH_CHECK(m_truthMuonTruthOrigin.initialize());
62 return StatusCode::SUCCESS;
63}
64
66
67bool MuonPerformanceAlg::passID(const xAOD::TrackParticle* tp, bool debug) const {
68 uint8_t nblh = 0x0;
69 uint8_t eblh = 0x0;
70 uint8_t nblo = 0x0;
71 uint8_t nphi = 0x0;
72 uint8_t npds = 0x0;
73 uint8_t nsctds = 0x0;
74 uint8_t nscthi = 0x0;
75 uint8_t nsctho = 0x0;
76 uint8_t npho = 0x0;
77 uint8_t ntrthi = 0x0;
78 uint8_t ntrtol = 0x0;
79
80 int tightness = 0;
81
82 bool passesIDcuts = false;
83
84 if (tp) {
85 // check blayer
86 tp->summaryValue(eblh, xAOD::expectInnermostPixelLayerHit);
87 tp->summaryValue(nblh, xAOD::numberOfInnermostPixelLayerHits);
88 tp->summaryValue(nblo, xAOD::numberOfInnermostPixelLayerOutliers);
89 if (static_cast<int>(eblh) == 0 || static_cast<int>(nblh) + static_cast<int>(nblo) > 0) tightness += 1;
90
91 // pixel hit counts
92 tp->summaryValue(nphi, xAOD::numberOfPixelHits);
93 tp->summaryValue(npds, xAOD::numberOfPixelDeadSensors);
94 if (static_cast<int>(nphi) + static_cast<int>(npds) >= 2) tightness += 2;
95
96 // sct hit counts
97 tp->summaryValue(nscthi, xAOD::numberOfSCTHits);
98 tp->summaryValue(nsctds, xAOD::numberOfSCTDeadSensors);
99 if (static_cast<int>(nscthi) + static_cast<int>(nsctds) > 5) tightness += 4;
100
101 // hole cuts
102 tp->summaryValue(npho, xAOD::numberOfPixelHoles);
103 tp->summaryValue(nsctho, xAOD::numberOfSCTHoles);
104 if (static_cast<int>(npho) + static_cast<int>(nsctho) < 2) tightness += 8;
105
106 // trt cuts
107 tp->summaryValue(ntrthi, xAOD::numberOfTRTHits);
108 tp->summaryValue(ntrtol, xAOD::numberOfTRTOutliers);
109 int ntrt = static_cast<int>(ntrthi) + static_cast<int>(ntrtol);
110 bool pass = true;
111 if (ntrt > 5 && static_cast<int>(ntrtol) / (double)ntrt > 0.9) pass = false;
112 if ((std::abs(tp->eta()) > 0.1 && std::abs(tp->eta()) < 1.9) && ntrt <= 5) pass = false;
113 if (pass) tightness += 16;
114 //
115 // B layer is not a requirement anymore
116 //
117 if (pass && tightness > 29) passesIDcuts = true;
118
119 if (passesIDcuts && debug) ATH_MSG_VERBOSE("ID trackparticle PASSED IDcuts with tightness " << tightness);
120 }
121 if (!passesIDcuts && debug) ATH_MSG_DEBUG("ID trackparticle FAILED IDcuts with tightness " << tightness);
122
123 return passesIDcuts;
124}
125
126StatusCode MuonPerformanceAlg::execute() {
127 SG::ReadHandle<xAOD::EventInfo> eventInfo(m_eventInfo);
128
129 m_runNumber = eventInfo->runNumber();
130 m_eventNumber = eventInfo->eventNumber();
131
132 SG::ReadHandle<xAOD::TruthParticleContainer> TruthMuons(m_truthMuons);
133 if (!TruthMuons.isPresent()) {
134 ATH_MSG_DEBUG("no truth muon collection present");
135 return StatusCode::SUCCESS;
136 }
137 if (!TruthMuons.isValid()) {
138 ATH_MSG_WARNING(m_truthMuons.key() << " not valid");
139 return StatusCode::FAILURE;
140 }
141 ATH_MSG_VERBOSE("Retrieved truth muons " << TruthMuons->size());
142 typedef ElementLink<xAOD::MuonContainer> MuonLink;
143
144 static const SG::ConstAccessor<int> truthTypeAcc("truthType");
145 static const SG::ConstAccessor<int> truthOriginAcc("truthOrigin");
146 static const SG::ConstAccessor<MuonLink> recoMuonLinkAcc("recoMuonLink");
147
148 m_nevents += 1;
149
150 static const SG::ConstAccessor<uint8_t> nprecLayersAcc("nprecLayers");
151
152 for (const auto truthMu : *TruthMuons) {
153 MuonLink link;
154 const int theType = truthTypeAcc(*truthMu);
155 const int theOrigin = truthOriginAcc(*truthMu);
156 ATH_MSG_VERBOSE("Truth muon: pt " << truthMu->pt() << " eta " << truthMu->eta());
157 ATH_MSG_VERBOSE("first loop: type " << theType << " origin " << theOrigin);
158 if (truthMu->pt() < 2000. || std::abs(truthMu->eta()) > 2.8) continue;
159 if (std::abs(truthMu->eta()) > 2.5)
160 ATH_MSG_VERBOSE(" SA |eta| > 2.5 muon with truth prec layers " << (int)nprecLayersAcc(*truthMu));
161
162 if (std::abs(truthMu->eta()) > 2.5 && (int)nprecLayersAcc(*truthMu) < 2) continue;
163 if (theType != 6 && theType != 7) continue;
164 if (theOrigin == 0 || theOrigin > 17) continue;
165 bool insideID = false;
166 if (std::abs(truthMu->eta()) < 2.0) insideID = true;
167 ATH_MSG_VERBOSE("Accepted Truth muon: pt " << truthMu->pt() << " eta " << truthMu->eta());
168
169 if (!insideID) m_ntruth[0] += 1;
170 if (!insideID) m_ntruth[9] += 1;
171 if (!insideID) m_ntruth[10] += 1;
172 if (insideID)
173 for (int n = 1; n < 6; n++) m_ntruth[n] += 1;
174 if (insideID)
175 for (int n = 7; n < 11; n++) m_ntruth[n] += 1;
176
177 if (truthMu->pt() > 5000.) {
178 if (!insideID) m_ntruth5[0] += 1;
179 if (!insideID) m_ntruth5[9] += 1;
180 if (!insideID) m_ntruth5[10] += 1;
181 if (insideID)
182 for (int n = 1; n < 6; n++) m_ntruth5[n] += 1;
183 if (insideID)
184 for (int n = 7; n < 11; n++) m_ntruth5[n] += 1;
185 }
186
187 if (truthMu->pt() > 10000.) {
188 if (!insideID) m_ntruth10[0] += 1;
189 if (!insideID) m_ntruth10[9] += 1;
190 if (!insideID) m_ntruth10[10] += 1;
191 if (insideID)
192 for (int n = 1; n < 6; n++) m_ntruth10[n] += 1;
193 if (insideID)
194 for (int n = 7; n < 11; n++) m_ntruth10[n] += 1;
195 }
196
197 if (recoMuonLinkAcc.isAvailable(*truthMu)) {
198 link = recoMuonLinkAcc(*truthMu);
199 ATH_MSG_VERBOSE(" link " << link.isValid());
200 if (link.isValid()) {
201 bool loose = false;
202 bool medium = false;
203 bool tight = false;
204 if (((*link)->quality() <= xAOD::Muon_v1::Loose)) loose = true;
205 if (((*link)->quality() <= xAOD::Muon_v1::Medium)) medium = true;
206 if (((*link)->quality() <= xAOD::Muon_v1::Tight)) tight = true;
207 if (!loose) ATH_MSG_WARNING("Muon did not pass Loose WP");
208 if (insideID) {
209 const xAOD::TrackParticle* tp = (*link)->primaryTrackParticle();
210 if (tp) {
211 } else
212 ATH_MSG_VERBOSE("No Track particle found on recoMuonLink");
213 // ID efficiency can only be measured on muons in the container
214
215 m_ntruth[6] += 1;
216 if (truthMu->pt() > 5000.) m_ntruth5[6] += 1;
217 if (truthMu->pt() > 10000.) m_ntruth10[6] += 1;
218
219 bool passesIDcuts = passID(tp, true);
220 ATH_MSG_VERBOSE(" all authors " << (*link)->allAuthors() << " passesIDcuts " << passesIDcuts);
221 if (passesIDcuts != (*link)->passesIDCuts()) {
222 ATH_MSG_DEBUG(" PROBLEM passedIDcuts from xAOD muon: " << (*link)->passesIDCuts()
223 << " BUT MuonPerformancAlg code gives " << passesIDcuts);
224 }
225 // Use xAOD Muon ID hit selection cuts
226 passesIDcuts = (*link)->passesIDCuts();
227 if (passesIDcuts) {
228 if (nprecLayersAcc(*truthMu) > 0) {
229 m_ntruth[11] += 1;
230 if (truthMu->pt() > 5000.) m_ntruth5[11] += 1;
231 if (truthMu->pt() > 10000.) m_ntruth10[11] += 1;
232 }
233 if (loose) {
234 m_nfound[10] += 1;
235 } else
236 print(" Muon not found by Loose ", truthMu);
237 if (medium) {
238 m_nfound[9] += 1;
239 } else if (loose)
240 print(" Muon not found by Medium ", truthMu);
241 if (tight) {
242 m_nfound[8] += 1;
243 } else if (medium)
244 print(" Muon not found by Tight ", truthMu);
245 m_nfound[6] += 1;
246 if (truthMu->pt() > 5000.) {
247 if (loose) m_nfound5[10] += 1;
248 if (medium) m_nfound5[9] += 1;
249 if (tight) m_nfound5[8] += 1;
250 m_nfound5[6] += 1;
251 }
252 if (truthMu->pt() > 10000.) {
253 if (loose) m_nfound10[10] += 1;
254 if (medium) m_nfound10[9] += 1;
255 if (tight) m_nfound10[8] += 1;
256 m_nfound10[6] += 1;
257 }
258 // CaloTag or Calolikelihood
259 if (((*link)->allAuthors() & 64 * 4) || ((*link)->allAuthors() & 64 * 8)) {
260 m_nfound[5] += 1;
261 if (truthMu->pt() > 5000.) m_nfound5[5] += 1;
262 if (truthMu->pt() > 10000.) m_nfound10[5] += 1;
263 } else
264 print(" Muon not found by CaloTag and Calolikelihood ", truthMu);
265 // MuidSA
266 if (((*link)->allAuthors() & 32)) {
267 m_nfound[7] += 1;
268 if (truthMu->pt() > 5000.) m_nfound5[7] += 1;
269 if (truthMu->pt() > 10000.) m_nfound10[7] += 1;
270 } else if (!((*link)->allAuthors() & 2))
271 print(" Muon not found by MuidSA ", truthMu);
272 // MuTag or MuTagIMO
273 if (((*link)->allAuthors() & 8) || ((*link)->allAuthors() & 16)) {
274 m_nfound[4] += 1;
275 if (truthMu->pt() > 5000.) m_nfound5[4] += 1;
276 if (truthMu->pt() > 10000.) m_nfound10[4] += 1;
277 } else
278 print(" Muon not found by MuTagIMO ", truthMu);
279 // MuGirl
280 if (((*link)->allAuthors() & 64)) {
281 m_nfound[3] += 1;
282 if (truthMu->pt() > 5000.) m_nfound5[3] += 1;
283 if (truthMu->pt() > 10000.) m_nfound10[3] += 1;
284 } else
285 print(" Muon not found by MuGirl ", truthMu);
286 // MuidCo
287 if (((*link)->allAuthors() & 2)) {
288 m_nfound[2] += 1;
289 if (truthMu->pt() > 5000.) m_nfound5[2] += 1;
290 if (truthMu->pt() > 10000.) m_nfound10[2] += 1;
291 } else
292 print(" Muon not found by MuidCo ", truthMu);
293 // MuidCo || STACO || Combined
294 if (((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4) ||
295 (*link)->muonType() == xAOD::Muon::MuonType::Combined) {
296 m_nfound[1] += 1;
297 if (truthMu->pt() > 5000.) m_nfound5[1] += 1;
298 if (truthMu->pt() > 10000.) m_nfound10[1] += 1;
299 if (nprecLayersAcc(*truthMu) > 0) {
300 m_nfound[11] += 1;
301 if (truthMu->pt() > 5000.) m_nfound5[11] += 1;
302 if (truthMu->pt() > 10000.) m_nfound10[11] += 1;
303 }
304 } else
305 print(" Combined Muon not found ", truthMu);
306 } else {
307 // SA or CB not passing ID cuts
308 if (((*link)->allAuthors() & 32) || (*link)->muonType() == xAOD::Muon::MuonType::Combined) {
309 m_nfound[7] += 1;
310 if (truthMu->pt() > 5000.) m_nfound5[7] += 1;
311 if (truthMu->pt() > 10000.) m_nfound10[7] += 1;
312 } else
313 print(" SA (CB) Muon not found ", truthMu);
314 }
315 }
316 // MuidSA
317 else if (!insideID) {
318 if ((((*link)->allAuthors() & 32)) || ((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4)) {
319 m_nfound[0] += 1;
320 if (truthMu->pt() > 5000.) m_nfound5[0] += 1;
321 if (truthMu->pt() > 10000.) m_nfound10[0] += 1;
322 } else
323 print(" Muon not found by MuidSA endcap ", truthMu);
324 if (loose) {
325 m_nfound[10] += 1;
326 if (truthMu->pt() > 5000.) m_nfound5[10] += 1;
327 if (truthMu->pt() > 10000.) m_nfound10[10] += 1;
328 } else
329 print(" Muon not found by Loose endcap ", truthMu);
330 if (medium) {
331 m_nfound[9] += 1;
332 if (truthMu->pt() > 5000.) m_nfound5[9] += 1;
333 if (truthMu->pt() > 10000.) m_nfound10[9] += 1;
334 } else if (loose)
335 print(" Muon not found by Medium endcap ", truthMu);
336 }
337 } else {
338 print(" No link Muon not found by CaloTag and Calolikelihood ", truthMu);
339 print(" No link Muon not found by MuidSA ", truthMu);
340 print(" No link Muon not found by MuTagIMO ", truthMu);
341 print(" No link Muon not found by MuGirl ", truthMu);
342 print(" No link Muon not found by MuidCo ", truthMu);
343 print(" No link Combined Muon not found ", truthMu);
344 }
345 } // end valid link
346 }
347
348 for (const auto truthMu : *TruthMuons) {
349 MuonLink link;
350 if (truthMu->pt() < 2000. || std::abs(truthMu->eta()) > 3.) continue;
351 const int theType = truthTypeAcc(*truthMu);
352 const int theOrigin = truthOriginAcc(*truthMu);
353 if (theType != 6 && theType != 7) continue;
354 if (theOrigin == 0 || theOrigin > 17) continue;
355 if (std::abs(truthMu->eta()) > 2.5 && (int)nprecLayersAcc(*truthMu) < 2) continue;
356 bool insideID = false;
357 if (recoMuonLinkAcc.isAvailable(*truthMu)) {
358 link = recoMuonLinkAcc(*truthMu);
359 if (link.isValid()) {
360 const xAOD::TrackParticle* tp = (*link)->primaryTrackParticle();
361 if (tp) {
362 if ((*link)->pt() < 2000. || std::abs((*link)->eta()) > 2.8) continue;
363 if (std::abs((*link)->eta()) < 2.0) insideID = true;
364 bool loose = false;
365 bool medium = false;
366 bool tight = false;
367 if (((*link)->quality() == xAOD::Muon_v1::Loose)) loose = true;
368 if (((*link)->quality() == xAOD::Muon_v1::Medium)) loose = true;
369 if (((*link)->quality() == xAOD::Muon_v1::Tight)) loose = true;
370 if (((*link)->quality() == xAOD::Muon_v1::Medium)) medium = true;
371 if (((*link)->quality() == xAOD::Muon_v1::Tight)) medium = true;
372 if (((*link)->quality() == xAOD::Muon_v1::Tight)) tight = true;
373 if (insideID) {
374 bool passesIDcuts = passID(tp, false);
375 // Use xAOD Muon ID hit selection cuts
376 passesIDcuts = (*link)->passesIDCuts();
377 if (passesIDcuts) {
378 if (loose) m_nfoundr[10] += 1;
379 if (medium) m_nfoundr[9] += 1;
380 if (tight) m_nfoundr[8] += 1;
381 m_nfoundr[6] += 1;
382 if (((*link)->allAuthors() & 64 * 4) || ((*link)->allAuthors() & 64 * 8)) m_nfoundr[5] += 1;
383 if (((*link)->allAuthors() & 32)) m_nfoundr[7] += 1;
384 if (((*link)->allAuthors() & 8) || ((*link)->allAuthors() & 16)) m_nfoundr[4] += 1;
385 if (((*link)->allAuthors() & 64)) m_nfoundr[3] += 1;
386 if (((*link)->allAuthors() & 2)) m_nfoundr[2] += 1;
387 if (((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4) ||
388 (*link)->muonType() == xAOD::Muon::MuonType::Combined)
389 m_nfoundr[1] += 1;
390 if ((*link)->pt() > 5000.) {
391 if (loose) m_nfoundr5[10] += 1;
392 if (medium) m_nfoundr5[9] += 1;
393 if (tight) m_nfoundr5[8] += 1;
394 m_nfoundr5[6] += 1;
395 if (((*link)->allAuthors() & 64 * 4) || ((*link)->allAuthors() & 64 * 8)) m_nfoundr5[5] += 1;
396 if (((*link)->allAuthors() & 32)) m_nfoundr5[7] += 1;
397 if (((*link)->allAuthors() & 8) || ((*link)->allAuthors() & 16)) m_nfoundr5[4] += 1;
398 if (((*link)->allAuthors() & 64)) m_nfoundr5[3] += 1;
399 if (((*link)->allAuthors() & 2)) m_nfoundr5[2] += 1;
400 if (((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4) ||
401 (*link)->muonType() == xAOD::Muon::MuonType::Combined)
402 m_nfoundr5[1] += 1;
403 }
404 if ((*link)->pt() > 10000.) {
405 if (loose) m_nfoundr10[10] += 1;
406 if (medium) m_nfoundr10[9] += 1;
407 if (tight) m_nfoundr10[8] += 1;
408 m_nfoundr10[6] += 1;
409 if (((*link)->allAuthors() & 64 * 4) || ((*link)->allAuthors() & 64 * 8)) m_nfoundr10[5] += 1;
410 if (((*link)->allAuthors() & 32)) m_nfoundr10[7] += 1;
411 if (((*link)->allAuthors() & 8) || ((*link)->allAuthors() & 16)) m_nfoundr10[4] += 1;
412 if (((*link)->allAuthors() & 64)) m_nfoundr10[3] += 1;
413 if (((*link)->allAuthors() & 2)) m_nfoundr10[2] += 1;
414 if (((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4) ||
415 (*link)->muonType() == xAOD::Muon::MuonType::Combined)
416 m_nfoundr10[1] += 1;
417 }
418 } else {
419 if (((*link)->allAuthors() & 32) || (*link)->muonType() == xAOD::Muon::MuonType::Combined) {
420 m_nfoundr[7] += 1;
421 if ((*link)->pt() > 5000.) m_nfoundr5[7] += 1;
422 if ((*link)->pt() > 10000.) m_nfoundr10[7] += 1;
423 }
424 }
425 } else if (!insideID) {
426 if ((((*link)->allAuthors() & 32)) || ((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4)) {
427 m_nfoundr[0] += 1;
428 if ((*link)->pt() > 5000.) m_nfoundr5[0] += 1;
429 if ((*link)->pt() > 10000.) m_nfoundr10[0] += 1;
430 }
431 if (loose) m_nfoundr[10] += 1;
432 if (medium) m_nfoundr[9] += 1;
433 if ((*link)->pt() > 5000.) {
434 if (loose) m_nfoundr5[10] += 1;
435 if (medium) m_nfoundr5[9] += 1;
436 }
437 if ((*link)->pt() > 10000.) {
438 if (loose) m_nfoundr10[10] += 1;
439 if (medium) m_nfoundr10[9] += 1;
440 }
441 }
442 } else
443 ATH_MSG_VERBOSE("No Track particle found on recoMuonLink");
444 } // not valid muonLink
445 } // no muonLink
446 }
447
448 SG::ReadHandle<xAOD::MuonContainer> Muons(m_muonsNameKey);
449
450 if (!Muons.isValid()) {
451 ATH_MSG_ERROR("Couldn't retrieve Muons container with key: " << m_muonsNameKey.key());
452 return StatusCode::FAILURE;
453 }
454 ATH_MSG_VERBOSE("Retrieved muons " << Muons->size());
455
456 for (const auto mu : *Muons) {
457 if (mu->pt() < 2000. || std::abs(mu->eta()) > 2.8) continue;
458 const xAOD::TrackParticle* tp = mu->primaryTrackParticle();
459
460 if (tp) {
461 ElementLink<xAOD::TruthParticleContainer> truthLink;
462 bool passesIDcuts = passID(tp, false);
463 // Use xAOD Muon ID hit selection cuts
464 passesIDcuts = mu->passesIDCuts();
465
466 bool insideID = false;
467 if (std::abs(mu->eta()) < 2.0) insideID = true;
468 static const SG::ConstAccessor<ElementLink<xAOD::TruthParticleContainer> >
469 truthParticleLinkAcc ("truthParticleLink");
470 if (truthParticleLinkAcc.isAvailable(*tp))
471 truthLink = truthParticleLinkAcc(*tp);
472 bool fake = true;
473 if (truthLink.isValid()) {
474 if (selectPdg((*truthLink)->pdgId())) fake = false;
475 }
476 bool loose = false;
477 bool medium = false;
478 bool tight = false;
479 if (mu->quality() == xAOD::Muon_v1::Loose) loose = true;
480 if (mu->quality() == xAOD::Muon_v1::Medium) loose = true;
481 if (mu->quality() == xAOD::Muon_v1::Tight) loose = true;
482 if (mu->quality() == xAOD::Muon_v1::Medium) medium = true;
483 if (mu->quality() == xAOD::Muon_v1::Tight) medium = true;
484 if (mu->quality() == xAOD::Muon_v1::Tight) tight = true;
485 if (insideID) {
486 if (passesIDcuts) {
487 if (loose) m_nreco[10] += 1;
488 if (medium) m_nreco[9] += 1;
489 if (tight) m_nreco[8] += 1;
490 m_nreco[6] += 1;
491 if ((mu->allAuthors() & 64 * 4) || (mu->allAuthors() & 64 * 8)) m_nreco[5] += 1;
492 if ((mu->allAuthors() & 32)) m_nreco[7] += 1;
493 if ((mu->allAuthors() & 8) || (mu->allAuthors() & 16)) m_nreco[4] += 1;
494 if ((mu->allAuthors() & 64)) m_nreco[3] += 1;
495 if ((mu->allAuthors() & 2)) m_nreco[2] += 1;
496 if ((mu->allAuthors() & 2) || (mu->allAuthors() & 4) || mu->muonType() == xAOD::Muon::MuonType::Combined)
497 m_nreco[1] += 1;
498 if (mu->pt() > 5000.) {
499 if (loose) m_nreco5[10] += 1;
500 if (medium) m_nreco5[9] += 1;
501 if (tight) m_nreco5[8] += 1;
502 m_nreco5[6] += 1;
503 if ((mu->allAuthors() & 64 * 4) || (mu->allAuthors() & 64 * 8)) m_nreco5[5] += 1;
504 if ((mu->allAuthors() & 32)) m_nreco5[7] += 1;
505 if ((mu->allAuthors() & 8) || (mu->allAuthors() & 16)) m_nreco5[4] += 1;
506 if ((mu->allAuthors() & 64)) m_nreco5[3] += 1;
507 if ((mu->allAuthors() & 2)) m_nreco5[2] += 1;
508 if ((mu->allAuthors() & 2) || (mu->allAuthors() & 4) || mu->muonType() == xAOD::Muon::MuonType::Combined)
509 m_nreco5[1] += 1;
510 }
511 if (mu->pt() > 10000.) {
512 if (loose) m_nreco10[10] += 1;
513 if (medium) m_nreco10[9] += 1;
514 if (tight) m_nreco10[8] += 1;
515 m_nreco10[6] += 1;
516 if ((mu->allAuthors() & 64 * 4) || (mu->allAuthors() & 64 * 8)) m_nreco10[5] += 1;
517 if ((mu->allAuthors() & 32)) m_nreco10[7] += 1;
518 if ((mu->allAuthors() & 8) || (mu->allAuthors() & 16)) m_nreco10[4] += 1;
519 if ((mu->allAuthors() & 64)) m_nreco10[3] += 1;
520 if ((mu->allAuthors() & 2)) m_nreco10[2] += 1;
521 if ((mu->allAuthors() & 2) || (mu->allAuthors() & 4) || mu->muonType() == xAOD::Muon::MuonType::Combined)
522 m_nreco10[1] += 1;
523 }
524 if (fake) {
525 if ((mu->quality() == xAOD::Muon_v1::Loose)) print(" Fake muon found by Loose ", mu);
526 if ((mu->quality() == xAOD::Muon_v1::Medium)) print(" Fake muon found by Medium ", mu);
527 if ((mu->quality() == xAOD::Muon_v1::Tight)) print(" Fake muon found by Tight ", mu);
528 if ((mu->allAuthors() & 64 * 4) || (mu->allAuthors() & 64 * 8))
529 print(" Fake muon found by CaloTag and Calolikelihood ", mu);
530 if ((mu->allAuthors() & 32)) print(" Fake muon found by MuidSA ", mu);
531 if ((mu->allAuthors() & 8) || (mu->allAuthors() & 16)) print(" Fake muon found by MuTagIMO ", mu);
532 if ((mu->allAuthors() & 64)) print(" Fake muon found by MuGirl ", mu);
533 if ((mu->allAuthors() & 2)) print(" Fake muon found by MuidCo ", mu);
534 if ((mu->allAuthors() & 2) || (mu->allAuthors() & 4) || mu->muonType() == xAOD::Muon::MuonType::Combined)
535 print(" Fake Combined muon ", mu);
536 }
537 } else {
538 if ((mu->allAuthors() & 32) || mu->muonType() == xAOD::Muon::MuonType::Combined) {
539 m_nreco[7] += 1;
540 if (fake) print(" Fake muon found by MuidSA with no ID", mu);
541 if (mu->pt() > 5000.) m_nreco5[7] += 1;
542 if (mu->pt() > 10000.) m_nreco10[7] += 1;
543 }
544 }
545 } else if (!insideID) {
546 if (loose) m_nreco[10] += 1;
547 if (medium) m_nreco[9] += 1;
548 if (loose && fake) print(" Fake muon found by Loose Endcap ", mu);
549 if (medium && fake && !loose) print(" Fake muon found by Medium Endcap ", mu);
550 if (tight && !medium && fake) print(" Fake muon found by Tight Endcap ", mu);
551 if (mu->pt() > 5000.) {
552 if (loose) m_nreco5[10] += 1;
553 if (medium) m_nreco5[9] += 1;
554 }
555 if (mu->pt() > 10000.) {
556 if (loose) m_nreco10[10] += 1;
557 if (medium) m_nreco10[9] += 1;
558 }
559 if (((mu->allAuthors() & 32)) || (mu->allAuthors() & 2) || (mu->allAuthors() & 4)) {
560 m_nreco[0] += 1;
561 if (fake) print(" Fake muon found by MuidSA Endcap ", mu);
562 if (mu->pt() > 5000.) m_nreco5[0] += 1;
563 if (mu->pt() > 10000.) m_nreco10[0] += 1;
564 }
565 }
566 }
567 }
568
569 return StatusCode::SUCCESS;
570}
571void MuonPerformanceAlg::print(const std::string& txt, const xAOD::TruthParticle* muon) {
572 // Truth: r 4251 z 3147 theta 0.91843 phi -0.781 q*p(GeV) 1.081e+02 pt(Gev) 8.589e+01
573 // uniqueID 1489 mother 23 production vertex: r 1.06363 z -102.414
574 int q = 1;
575 if (muon->pdgId() > 0) q = -1;
576 double p = sqrt(muon->e() * muon->e() - muon->m() * muon->m());
577 static const SG::ConstAccessor<uint8_t> nprecLayersAcc("nprecLayers");
578 ATH_MSG_DEBUG(txt << " run " << m_runNumber << " event " << m_eventNumber << std::endl
579 << " Truth: pdgId " << muon->pdgId() << " uniqueID " << HepMC::uniqueID(muon) << " eta " << muon->eta() << " phi "
580 << muon->phi() << " q*p (GeV) " << q * p / 1000. << " pt (GeV) " << muon->pt() / 1000. << " precisionLayers "
581 << static_cast<int>(nprecLayersAcc(*muon)));
582}
583void MuonPerformanceAlg::print(const std::string& txt, const xAOD::Muon* muon) {
584 int nprec = 0;
585 uint8_t nPrecision = 0;
586 if (muon->primaryTrackParticleLink().isValid()) {
587 muon->primaryTrackParticle()->summaryValue(nPrecision, xAOD::numberOfPrecisionLayers);
588 nprec = static_cast<int>(nPrecision);
589 }
590 double p = sqrt(muon->e() * muon->e() - muon->m() * muon->m());
591 ATH_MSG_DEBUG(txt << " run " << m_runNumber << " event " << m_eventNumber << std::endl
592 << " eta " << muon->eta() << " phi " << muon->phi() << " q*p (GeV) " << muon->charge() * p / 1000. << " pt (GeV) "
593 << muon->pt() / 1000. << " precisionLayers " << nprec << " nr segments " << muon->nMuonSegments());
594}
595
596StatusCode MuonPerformanceAlg::finalize() {
597 std::ofstream fileOutput;
598 std::string outfile = "muonPerformance_xAOD.txt";
599 fileOutput.open(outfile.c_str(), std::ios::trunc);
600 std::ostringstream sout;
601 sout.precision(4);
602
603 unsigned int width = 9;
604 unsigned int precision = 3;
605
606 sout << std::endl;
607 sout << " Summary of the xAOD Muon performance " << std::endl;
608 sout << " Muon type #Truth muons Efficiency #Truth muons Efficiency #Truth muons Efficiency" << std::endl;
609 sout << " (pt>2) (pt>5) (pt>10 GeV/c) " << std::endl;
610 for (unsigned int i = 0; i < 12; ++i) {
611 if (m_ntruth[i] == 0) {
612 sout << " " << std::endl;
613 } else {
614 sout << " " << m_hitCutString[i];
615 sout << std::setw(width) << std::setprecision(precision);
616 sout << static_cast<double>(m_ntruth[i]);
617 sout << " ";
618 sout << std::setw(width) << std::setprecision(precision);
619 sout << static_cast<double>(m_nfound[i]) / static_cast<double>(m_ntruth[i]);
620 sout << std::setw(width) << std::setprecision(precision);
621 sout << static_cast<double>(m_ntruth5[i]);
622 sout << " ";
623 sout << std::setw(width) << std::setprecision(precision);
624 if (m_ntruth5[i] != 0) {
625 sout << static_cast<double>(m_nfound5[i]) / static_cast<double>(m_ntruth5[i]);
626 } else
627 sout << 0.;
628
629 sout << std::setw(width) << std::setprecision(precision);
630 sout << static_cast<double>(m_ntruth10[i]);
631 sout << " ";
632 sout << std::setw(width) << std::setprecision(precision);
633 if (m_ntruth10[i] != 0) {
634 sout << static_cast<double>(m_nfound10[i]) / static_cast<double>(m_ntruth10[i]) << std::endl;
635 } else
636 sout << 0. << std::endl;
637 }
638 }
639 sout << " The efficiency of the ID is calculated inside |eta| < 2 for the MCP ID hit selection cuts (it uses identified - Combined, "
640 "Tagged or CaloTagged - muons with an ID track)"
641 << std::endl;
642 sout << " The efficiencies for CB all, MuidCB, MuGirl, Tag and Calo include the MCP ID cuts" << std::endl;
643 sout << " The Tight, Medium and Loose efficiencies include MCP ID cuts for muons |eta| < 2" << std::endl;
644 sout << " The SA 2.0 for |eta| >2 and SA (no ID) for |eta| < 2 doesnot include MCP ID cuts" << std::endl;
645 sout << " The Combined efficiency is defined ID tracks after MCP cuts and requiring at least 1 muon station at truth level"
646 << std::endl;
647 sout << " Fakes are calculated with the selections listed above" << std::endl;
648 sout << std::endl;
649
650 sout << " Muon type #Fake muons rate #Fake muons rate #Fake muons rate " << std::endl;
651 sout << " (all pt>0,2) (pt>5) (pt>10 GeV/c) " << std::endl;
652 for (unsigned int i = 0; i < 11; ++i) {
653 if (i == 6) continue;
654 if (m_nevents == 0) {
655 sout << " " << std::endl;
656 } else {
657 sout << " " << m_hitCutString[i];
658 sout << std::setw(width) << std::setprecision(precision);
659 sout << static_cast<double>(m_nreco[i] - m_nfoundr[i]);
660 sout << std::setw(width) << std::setprecision(precision);
661 sout << static_cast<double>(m_nreco[i] - m_nfoundr[i]) / static_cast<double>(m_nevents);
662 sout << std::setw(width) << std::setprecision(precision);
663 sout << static_cast<double>(m_nreco5[i] - m_nfoundr5[i]);
664 sout << std::setw(width) << std::setprecision(precision);
665 sout << static_cast<double>(m_nreco5[i] - m_nfoundr5[i]) / static_cast<double>(m_nevents);
666 sout << std::setw(width) << std::setprecision(precision);
667 sout << static_cast<double>(m_nreco10[i] - m_nfoundr10[i]);
668 sout << std::setw(width) << std::setprecision(precision);
669 sout << static_cast<double>(m_nreco10[i] - m_nfoundr10[i]) / static_cast<double>(m_nevents) << std::endl;
670 }
671 }
672 sout << std::endl;
673 fileOutput << sout.str() << std::endl;
674 fileOutput.close();
675 return StatusCode::SUCCESS;
676}
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition AthMsgStreamMacros.h:33
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition AthMsgStreamMacros.h:28
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition AthMsgStreamMacros.h:32
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
MuonLink
ElementLink< xAOD::MuonContainer > MuonLink
Definition BPhysHelper.cxx:21
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
debug
const bool debug
Definition MakeUncertaintyPlots.cxx:53
width
const double width
Definition TTileTripReader.cxx:24
AthAlgorithm::AthAlgorithm
AthAlgorithm(const std::string &name, ISvcLocator *pSvcLocator)
Constructor with parameters:
Definition AthAlgorithm.cxx:24
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T, V, H > &t)
Definition AthCommonDataStore.h:145
MuonPerformanceAlg::m_nfound5
std::vector< int > m_nfound5
Definition MuonPerformanceAlg.h:57
MuonPerformanceAlg::finalize
virtual StatusCode finalize() override
Definition MuonPerformanceAlg.cxx:596
MuonPerformanceAlg::m_nreco5
std::vector< int > m_nreco5
Definition MuonPerformanceAlg.h:63
MuonPerformanceAlg::MuonPerformanceAlg
MuonPerformanceAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition MuonPerformanceAlg.cxx:20
MuonPerformanceAlg::m_nreco10
std::vector< int > m_nreco10
Definition MuonPerformanceAlg.h:64
MuonPerformanceAlg::m_ntruth10
std::vector< int > m_ntruth10
Definition MuonPerformanceAlg.h:55
MuonPerformanceAlg::initialize
virtual StatusCode initialize() override
Definition MuonPerformanceAlg.cxx:27
MuonPerformanceAlg::m_selectedPdgs
std::set< int > m_selectedPdgs
Definition MuonPerformanceAlg.h:79
MuonPerformanceAlg::m_nfoundr5
std::vector< int > m_nfoundr5
Definition MuonPerformanceAlg.h:60
MuonPerformanceAlg::m_ntruth5
std::vector< int > m_ntruth5
Definition MuonPerformanceAlg.h:54
MuonPerformanceAlg::m_pdgsToBeConsidered
IntegerArrayProperty m_pdgsToBeConsidered
Definition MuonPerformanceAlg.h:78
MuonPerformanceAlg::m_nfoundr
std::vector< int > m_nfoundr
Definition MuonPerformanceAlg.h:59
MuonPerformanceAlg::m_ntruth
std::vector< int > m_ntruth
Definition MuonPerformanceAlg.h:53
MuonPerformanceAlg::m_eventInfo
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo
Definition MuonPerformanceAlg.h:66
MuonPerformanceAlg::selectPdg
bool selectPdg(int pdg) const
Definition MuonPerformanceAlg.h:80
MuonPerformanceAlg::m_nfoundr10
std::vector< int > m_nfoundr10
Definition MuonPerformanceAlg.h:61
MuonPerformanceAlg::m_truthMuonRecoMuonLinkKey
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_truthMuonRecoMuonLinkKey
Definition MuonPerformanceAlg.h:68
MuonPerformanceAlg::execute
virtual StatusCode execute() override
Definition MuonPerformanceAlg.cxx:126
MuonPerformanceAlg::m_truthMuonTruthType
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_truthMuonTruthType
Definition MuonPerformanceAlg.h:69
MuonPerformanceAlg::m_muonsNameKey
SG::ReadHandleKey< xAOD::MuonContainer > m_muonsNameKey
Definition MuonPerformanceAlg.h:49
MuonPerformanceAlg::m_nreco
std::vector< int > m_nreco
Definition MuonPerformanceAlg.h:62
MuonPerformanceAlg::m_hitCutString
std::vector< std::string > m_hitCutString
Definition MuonPerformanceAlg.h:52
MuonPerformanceAlg::m_nfound10
std::vector< int > m_nfound10
Definition MuonPerformanceAlg.h:58
MuonPerformanceAlg::m_truthMuons
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_truthMuons
Definition MuonPerformanceAlg.h:67
MuonPerformanceAlg::passID
bool passID(const xAOD::TrackParticle *, bool debug) const
Definition MuonPerformanceAlg.cxx:67
MuonPerformanceAlg::m_nfound
std::vector< int > m_nfound
Definition MuonPerformanceAlg.h:56
MuonPerformanceAlg::print
void print(const std::string &txt, const xAOD::TruthParticle *muon)
Definition MuonPerformanceAlg.cxx:571
MuonPerformanceAlg::m_truthMuonTruthOrigin
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_truthMuonTruthOrigin
Definition MuonPerformanceAlg.h:70
MuonPerformanceAlg::m_muonTruthParticleKey
SG::ReadDecorHandleKey< xAOD::MuonContainer > m_muonTruthParticleKey
Definition MuonPerformanceAlg.h:50
MuonPerformanceAlg::m_writeToFile
bool m_writeToFile
name of external file to write statistics
Definition MuonPerformanceAlg.h:42
SG::ConstAccessor
Helper class to provide constant type-safe access to aux data.
Definition ConstAccessor.h:55
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
SG::VarHandleBase::isPresent
bool isPresent() const
Is the referenced object present in SG?
Definition StoreGate/src/VarHandleBase.cxx:400
HepMC::uniqueID
int uniqueID(const T &p)
Definition MagicNumbers.h:117
Muons
Definition Muons.py:1
xAOD::TrackParticle
TrackParticle_v1 TrackParticle
Reference the current persistent version:
Definition Event/xAOD/xAODTracking/xAODTracking/TrackParticle.h:13
xAOD::TruthParticle
TruthParticle_v1 TruthParticle
Typedef to implementation.
Definition Event/xAOD/xAODTruth/xAODTruth/TruthParticle.h:15
xAOD::Muon
Muon_v1 Muon
Reference the current persistent version:
Definition Event/xAOD/xAODMuon/xAODMuon/Muon.h:13
xAOD::expectInnermostPixelLayerHit
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
Definition TrackingPrimitives.h:237
xAOD::numberOfPixelHoles
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
Definition TrackingPrimitives.h:262
xAOD::numberOfPrecisionLayers
@ numberOfPrecisionLayers
layers with at least 3 hits [unit8_t].
Definition TrackingPrimitives.h:289
xAOD::numberOfTRTHits
@ numberOfTRTHits
number of TRT hits [unit8_t].
Definition TrackingPrimitives.h:276
xAOD::numberOfSCTDeadSensors
@ numberOfSCTDeadSensors
number of dead SCT sensors crossed [unit8_t].
Definition TrackingPrimitives.h:274
xAOD::numberOfInnermostPixelLayerOutliers
@ numberOfInnermostPixelLayerOutliers
number of 0th layer barrel outliers
Definition TrackingPrimitives.h:239
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition TrackingPrimitives.h:269
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition TrackingPrimitives.h:238
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition TrackingPrimitives.h:260
xAOD::numberOfTRTOutliers
@ numberOfTRTOutliers
number of TRT outliers [unit8_t].
Definition TrackingPrimitives.h:277
xAOD::numberOfPixelDeadSensors
@ numberOfPixelDeadSensors
number of dead pixel sensors crossed [unit8_t].
Definition TrackingPrimitives.h:267
xAOD::numberOfSCTHoles
@ numberOfSCTHoles
number of SCT holes [unit8_t].
Definition TrackingPrimitives.h:271
Generated on for ATLAS Offline Software by doxygen 1.14.0