ATLAS Offline Software
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 
10 #include "GaudiKernel/IToolSvc.h"
11 #include "xAODMuon/Muon.h"
12 #include "xAODMuon/MuonSegment.h"
19 
20 MuonPerformanceAlg::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 
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 
58  ATH_CHECK(m_muonTruthParticleKey.initialize());
62  return StatusCode::SUCCESS;
63 }
64 
66 
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 
128 
129  m_runNumber = eventInfo->runNumber();
130  m_eventNumber = eventInfo->eventNumber();
131 
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());
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 
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) {
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;
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 }
571 void 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 }
583 void 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 
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 }
MuonPerformanceAlg::m_nevents
unsigned int m_nevents
Definition: MuonPerformanceAlg.h:51
MuonPerformanceAlg::print
void print(const std::string &txt, const xAOD::TruthParticle *muon)
Definition: MuonPerformanceAlg.cxx:571
MuonPerformanceAlg::m_fileName
std::string m_fileName
Definition: MuonPerformanceAlg.h:43
xAOD::numberOfPixelHoles
@ numberOfPixelHoles
number of pixel layers on track with absence of hits [unit8_t].
Definition: TrackingPrimitives.h:261
LikeEnum::Loose
@ Loose
Definition: LikelihoodEnums.h:12
xAOD::muon
@ muon
Definition: TrackingPrimitives.h:195
MuonLink
ElementLink< xAOD::MuonContainer > MuonLink
Definition: BPhysHelper.cxx:21
MuonPerformanceAlg::m_muonsNameKey
SG::ReadHandleKey< xAOD::MuonContainer > m_muonsNameKey
Definition: MuonPerformanceAlg.h:49
AthCheckMacros.h
xAOD::uint8_t
uint8_t
Definition: Muon_v1.cxx:557
MuonPerformanceAlg::selectPdg
bool selectPdg(int pdg) const
Definition: MuonPerformanceAlg.h:80
xAOD::EventInfo_v1::eventNumber
uint64_t eventNumber() const
The current event's event number.
Muon.h
SG::ReadHandle
Definition: StoreGate/StoreGate/ReadHandle.h:70
AthCommonDataStore< AthCommonMsg< Algorithm > >::declareProperty
Gaudi::Details::PropertyBase & declareProperty(Gaudi::Property< T > &t)
Definition: AthCommonDataStore.h:145
MuonPerformanceAlg.h
TruthParticleContainer.h
MuonSegment.h
MuonPerformanceAlg::m_nfoundr10
std::vector< int > m_nfoundr10
Definition: MuonPerformanceAlg.h:61
MuonPerformanceAlg::m_nreco
std::vector< int > m_nreco
Definition: MuonPerformanceAlg.h:62
ParticleTest.tp
tp
Definition: ParticleTest.py:25
xAOD::numberOfPixelHits
@ numberOfPixelHits
these are the pixel hits, including the b-layer [unit8_t].
Definition: TrackingPrimitives.h:259
xAOD::expectInnermostPixelLayerHit
@ expectInnermostPixelLayerHit
Do we expect a 0th-layer barrel hit for this track?
Definition: TrackingPrimitives.h:236
xAOD::numberOfTRTHits
@ numberOfTRTHits
number of TRT hits [unit8_t].
Definition: TrackingPrimitives.h:275
SG::ConstAccessor< int >
MuonPerformanceAlg::m_writeToFile
bool m_writeToFile
name of external file to write statistics
Definition: MuonPerformanceAlg.h:42
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
SG::VarHandleKey::key
const std::string & key() const
Return the StoreGate ID for the referenced object.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:141
MuonPerformanceAlg::m_nfound
std::vector< int > m_nfound
Definition: MuonPerformanceAlg.h:56
SCT_CalibAlgs::nbins
@ nbins
Definition: SCT_CalibNumbers.h:10
xAOD::Muon_v1
Class describing a Muon.
Definition: Muon_v1.h:38
xAOD::EventInfo_v1::runNumber
uint32_t runNumber() const
The current event's run number.
MuonPerformanceAlg::m_selectedPdgs
std::set< int > m_selectedPdgs
Definition: MuonPerformanceAlg.h:79
MuonSegmentContainer.h
MuonPerformanceAlg::m_runNumber
int m_runNumber
Definition: MuonPerformanceAlg.h:72
MuonPerformanceAlg::m_muonTruthParticleKey
SG::ReadDecorHandleKey< xAOD::MuonContainer > m_muonTruthParticleKey
Definition: MuonPerformanceAlg.h:50
TruthParticleAuxContainer.h
MuonPerformanceAlg::m_nfoundr5
std::vector< int > m_nfoundr5
Definition: MuonPerformanceAlg.h:60
xAOD::numberOfInnermostPixelLayerOutliers
@ numberOfInnermostPixelLayerOutliers
number of 0th layer barrel outliers
Definition: TrackingPrimitives.h:238
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
lumiFormat.i
int i
Definition: lumiFormat.py:85
SG::ReadDecorHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
xAOD::TruthParticle_v1
Class describing a truth particle in the MC record.
Definition: TruthParticle_v1.h:37
HepMC::uniqueID
int uniqueID(const T &p)
Definition: MagicNumbers.h:109
jobOptions.fileOutput
fileOutput
Definition: jobOptions.py:39
MuonPerformanceAlg::m_eventNumber
int m_eventNumber
Definition: MuonPerformanceAlg.h:73
MuonPerformanceAlg::m_nreco10
std::vector< int > m_nreco10
Definition: MuonPerformanceAlg.h:64
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
MuonPerformanceAlg::m_eventInfo
SG::ReadHandleKey< xAOD::EventInfo > m_eventInfo
Definition: MuonPerformanceAlg.h:66
SG::VarHandleKey::initialize
StatusCode initialize(bool used=true)
If this object is used as a property, then this should be called during the initialize phase.
Definition: AthToolSupport/AsgDataHandles/Root/VarHandleKey.cxx:103
xAOD::double
double
Definition: CompositeParticle_v1.cxx:159
xAOD::numberOfSCTHoles
@ numberOfSCTHoles
number of SCT holes [unit8_t].
Definition: TrackingPrimitives.h:270
AthAlgorithm
Definition: AthAlgorithm.h:47
SG::ReadHandle::isValid
virtual bool isValid() override final
Can the handle be successfully dereferenced?
MuonPerformanceAlg::m_truthMuons
SG::ReadHandleKey< xAOD::TruthParticleContainer > m_truthMuons
Definition: MuonPerformanceAlg.h:67
Trk::Combined
@ Combined
Definition: TrackSummaryTool.h:32
MuonPerformanceAlg::m_ntruth
std::vector< int > m_ntruth
Definition: MuonPerformanceAlg.h:53
LikeEnum::Tight
@ Tight
Definition: LikelihoodEnums.h:15
MuonPerformanceAlg::m_nfound10
std::vector< int > m_nfound10
Definition: MuonPerformanceAlg.h:58
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:221
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
MagicNumbers.h
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_ntruth5
std::vector< int > m_ntruth5
Definition: MuonPerformanceAlg.h:54
MuonPerformanceAlg::m_ntruth10
std::vector< int > m_ntruth10
Definition: MuonPerformanceAlg.h:55
MuonPerformanceAlg::MuonPerformanceAlg
MuonPerformanceAlg(const std::string &name, ISvcLocator *pSvcLocator)
Definition: MuonPerformanceAlg.cxx:20
xAOD::numberOfTRTOutliers
@ numberOfTRTOutliers
number of TRT outliers [unit8_t].
Definition: TrackingPrimitives.h:276
Base_Fragment.width
width
Definition: Sherpa_i/share/common/Base_Fragment.py:59
Muons
Definition: Muons.py:1
LikeEnum::Medium
@ Medium
Definition: LikelihoodEnums.h:14
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
xAOD::numberOfSCTDeadSensors
@ numberOfSCTDeadSensors
number of dead SCT sensors crossed [unit8_t].
Definition: TrackingPrimitives.h:273
MuonPerformanceAlg::m_nfoundr
std::vector< int > m_nfoundr
Definition: MuonPerformanceAlg.h:59
MuonPerformanceAlg::m_truthMuonTruthOrigin
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_truthMuonTruthOrigin
Definition: MuonPerformanceAlg.h:70
extractSporadic.q
list q
Definition: extractSporadic.py:98
xAOD::numberOfSCTHits
@ numberOfSCTHits
number of hits in SCT [unit8_t].
Definition: TrackingPrimitives.h:268
SG::ConstAccessor::isAvailable
bool isAvailable(const ELT &e) const
Test to see if this variable exists in the store.
SG::VarHandleBase::isPresent
bool isPresent() const
Is the referenced object present in SG?
Definition: StoreGate/src/VarHandleBase.cxx:394
MuonPerformanceAlg::finalize
virtual StatusCode finalize() override
Definition: MuonPerformanceAlg.cxx:596
xAOD::numberOfPixelDeadSensors
@ numberOfPixelDeadSensors
number of dead pixel sensors crossed [unit8_t].
Definition: TrackingPrimitives.h:266
MuonPerformanceAlg::m_hitCutString
std::vector< std::string > m_hitCutString
Definition: MuonPerformanceAlg.h:52
xAOD::TrackParticle_v1
Class describing a TrackParticle.
Definition: TrackParticle_v1.h:43
ConstAccessor.h
Helper class to provide constant type-safe access to aux data.
TruthParticle.h
CaloNoise_fillDB.mu
mu
Definition: CaloNoise_fillDB.py:53
PrepareReferenceFile.outfile
outfile
Definition: PrepareReferenceFile.py:42
MuonPerformanceAlg::m_truthMuonTruthType
SG::ReadDecorHandleKey< xAOD::TruthParticleContainer > m_truthMuonTruthType
Definition: MuonPerformanceAlg.h:69
DataVector::size
size_type size() const noexcept
Returns the number of elements in the collection.
MuonPerformanceAlg::m_nreco5
std::vector< int > m_nreco5
Definition: MuonPerformanceAlg.h:63
MuonPerformanceAlg::m_pdgsToBeConsidered
IntegerArrayProperty m_pdgsToBeConsidered
Definition: MuonPerformanceAlg.h:78
MuonPerformanceAlg::m_nfound5
std::vector< int > m_nfound5
Definition: MuonPerformanceAlg.h:57
MuonPerformanceAlg::passID
bool passID(const xAOD::TrackParticle *, bool debug) const
Definition: MuonPerformanceAlg.cxx:67
xAOD::numberOfInnermostPixelLayerHits
@ numberOfInnermostPixelLayerHits
these are the hits in the 0th pixel barrel layer
Definition: TrackingPrimitives.h:237
MuonPerformanceAlg::initialize
virtual StatusCode initialize() override
Definition: MuonPerformanceAlg.cxx:27
xAOD::numberOfPrecisionLayers
@ numberOfPrecisionLayers
layers with at least 3 hits [unit8_t].
Definition: TrackingPrimitives.h:288