10 #include "GaudiKernel/IToolSvc.h"
21 AthAlgorithm(
name, pSvcLocator), m_writeToFile(false), m_nevents(0), m_runNumber(0), m_eventNumber(0) {
29 unsigned int nbins = 12;
43 m_hitCutString = {
"SA 2.0 ",
"CB all ",
"MuidCB ",
"MuGirl ",
"Tag ",
"Calo ",
44 "ID ",
"SA(no ID)",
"Tight ",
"Medium ",
"Loose ",
"Combined "};
62 return StatusCode::SUCCESS;
82 bool passesIDcuts =
false;
89 if (
static_cast<int>(eblh) == 0 ||
static_cast<int>(nblh) +
static_cast<int>(nblo) > 0) tightness += 1;
94 if (
static_cast<int>(nphi) +
static_cast<int>(npds) >= 2) tightness += 2;
99 if (
static_cast<int>(nscthi) +
static_cast<int>(nsctds) > 5) tightness += 4;
104 if (
static_cast<int>(npho) +
static_cast<int>(nsctho) < 2) tightness += 8;
109 int ntrt =
static_cast<int>(ntrthi) +
static_cast<int>(ntrtol);
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;
117 if (pass && tightness > 29) passesIDcuts =
true;
119 if (passesIDcuts &&
debug)
ATH_MSG_VERBOSE(
"ID trackparticle PASSED IDcuts with tightness " << tightness);
121 if (!passesIDcuts &&
debug)
ATH_MSG_DEBUG(
"ID trackparticle FAILED IDcuts with tightness " << tightness);
135 return StatusCode::SUCCESS;
139 return StatusCode::FAILURE;
152 for (
const auto truthMu : *TruthMuons) {
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));
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());
177 if (truthMu->pt() > 5000.) {
187 if (truthMu->pt() > 10000.) {
198 link = recoMuonLinkAcc(*truthMu);
216 if (truthMu->pt() > 5000.)
m_ntruth5[6] += 1;
217 if (truthMu->pt() > 10000.)
m_ntruth10[6] += 1;
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);
226 passesIDcuts = (*link)->passesIDCuts();
228 if (nprecLayersAcc(*truthMu) > 0) {
230 if (truthMu->pt() > 5000.)
m_ntruth5[11] += 1;
231 if (truthMu->pt() > 10000.)
m_ntruth10[11] += 1;
236 print(
" Muon not found by Loose ", truthMu);
240 print(
" Muon not found by Medium ", truthMu);
244 print(
" Muon not found by Tight ", truthMu);
246 if (truthMu->pt() > 5000.) {
252 if (truthMu->pt() > 10000.) {
259 if (((*link)->allAuthors() & 64 * 4) || ((*link)->allAuthors() & 64 * 8)) {
261 if (truthMu->pt() > 5000.)
m_nfound5[5] += 1;
262 if (truthMu->pt() > 10000.)
m_nfound10[5] += 1;
264 print(
" Muon not found by CaloTag and Calolikelihood ", truthMu);
266 if (((*link)->allAuthors() & 32)) {
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);
273 if (((*link)->allAuthors() & 8) || ((*link)->allAuthors() & 16)) {
275 if (truthMu->pt() > 5000.)
m_nfound5[4] += 1;
276 if (truthMu->pt() > 10000.)
m_nfound10[4] += 1;
278 print(
" Muon not found by MuTagIMO ", truthMu);
280 if (((*link)->allAuthors() & 64)) {
282 if (truthMu->pt() > 5000.)
m_nfound5[3] += 1;
283 if (truthMu->pt() > 10000.)
m_nfound10[3] += 1;
285 print(
" Muon not found by MuGirl ", truthMu);
287 if (((*link)->allAuthors() & 2)) {
289 if (truthMu->pt() > 5000.)
m_nfound5[2] += 1;
290 if (truthMu->pt() > 10000.)
m_nfound10[2] += 1;
292 print(
" Muon not found by MuidCo ", truthMu);
294 if (((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4) ||
297 if (truthMu->pt() > 5000.)
m_nfound5[1] += 1;
298 if (truthMu->pt() > 10000.)
m_nfound10[1] += 1;
299 if (nprecLayersAcc(*truthMu) > 0) {
301 if (truthMu->pt() > 5000.)
m_nfound5[11] += 1;
302 if (truthMu->pt() > 10000.)
m_nfound10[11] += 1;
305 print(
" Combined Muon not found ", truthMu);
310 if (truthMu->pt() > 5000.)
m_nfound5[7] += 1;
311 if (truthMu->pt() > 10000.)
m_nfound10[7] += 1;
313 print(
" SA (CB) Muon not found ", truthMu);
317 else if (!insideID) {
318 if ((((*link)->allAuthors() & 32)) || ((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4)) {
320 if (truthMu->pt() > 5000.)
m_nfound5[0] += 1;
321 if (truthMu->pt() > 10000.)
m_nfound10[0] += 1;
323 print(
" Muon not found by MuidSA endcap ", truthMu);
326 if (truthMu->pt() > 5000.)
m_nfound5[10] += 1;
327 if (truthMu->pt() > 10000.)
m_nfound10[10] += 1;
329 print(
" Muon not found by Loose endcap ", truthMu);
332 if (truthMu->pt() > 5000.)
m_nfound5[9] += 1;
333 if (truthMu->pt() > 10000.)
m_nfound10[9] += 1;
335 print(
" Muon not found by Medium endcap ", truthMu);
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);
348 for (
const auto truthMu : *TruthMuons) {
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;
358 link = recoMuonLinkAcc(*truthMu);
362 if ((*link)->pt() < 2000. || std::abs((*link)->eta()) > 2.8)
continue;
363 if (std::abs((*link)->eta()) < 2.0) insideID =
true;
374 bool passesIDcuts =
passID(
tp,
false);
376 passesIDcuts = (*link)->passesIDCuts();
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) ||
390 if ((*link)->pt() > 5000.) {
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) ||
404 if ((*link)->pt() > 10000.) {
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) ||
421 if ((*link)->pt() > 5000.)
m_nfoundr5[7] += 1;
425 }
else if (!insideID) {
426 if ((((*link)->allAuthors() & 32)) || ((*link)->allAuthors() & 2) || ((*link)->allAuthors() & 4)) {
428 if ((*link)->pt() > 5000.)
m_nfoundr5[0] += 1;
433 if ((*link)->pt() > 5000.) {
437 if ((*link)->pt() > 10000.) {
450 if (!
Muons.isValid()) {
452 return StatusCode::FAILURE;
457 if (
mu->pt() < 2000. || std::abs(
mu->eta()) > 2.8)
continue;
462 bool passesIDcuts =
passID(
tp,
false);
464 passesIDcuts =
mu->passesIDCuts();
466 bool insideID =
false;
467 if (std::abs(
mu->eta()) < 2.0) insideID =
true;
469 truthParticleLinkAcc (
"truthParticleLink");
471 truthLink = truthParticleLinkAcc(*
tp);
474 if (
selectPdg((*truthLink)->pdgId())) fake =
false;
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;
498 if (
mu->pt() > 5000.) {
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;
511 if (
mu->pt() > 10000.) {
516 if ((
mu->allAuthors() & 64 * 4) || (
mu->allAuthors() & 64 * 8))
m_nreco10[5] += 1;
518 if ((
mu->allAuthors() & 8) || (
mu->allAuthors() & 16))
m_nreco10[4] += 1;
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);
535 print(
" Fake Combined muon ",
mu);
540 if (fake)
print(
" Fake muon found by MuidSA with no ID",
mu);
545 }
else if (!insideID) {
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.) {
555 if (
mu->pt() > 10000.) {
559 if (((
mu->allAuthors() & 32)) || (
mu->allAuthors() & 2) || (
mu->allAuthors() & 4)) {
561 if (fake)
print(
" Fake muon found by MuidSA Endcap ",
mu);
569 return StatusCode::SUCCESS;
575 if (
muon->pdgId() > 0)
q = -1;
580 <<
muon->phi() <<
" q*p (GeV) " <<
q *
p / 1000. <<
" pt (GeV) " <<
muon->pt() / 1000. <<
" precisionLayers "
581 <<
static_cast<int>(nprecLayersAcc(*
muon)));
586 if (
muon->primaryTrackParticleLink().isValid()) {
588 nprec =
static_cast<int>(nPrecision);
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());
598 std::string
outfile =
"muonPerformance_xAOD.txt";
600 std::ostringstream sout;
603 unsigned int width = 9;
604 unsigned int precision = 3;
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) {
612 sout <<
" " << std::endl;
615 sout << std::setw(
width) << std::setprecision(precision);
616 sout << static_cast<double>(
m_ntruth[
i]);
618 sout << std::setw(
width) << std::setprecision(precision);
620 sout << std::setw(
width) << std::setprecision(precision);
623 sout << std::setw(
width) << std::setprecision(precision);
629 sout << std::setw(
width) << std::setprecision(precision);
632 sout << std::setw(
width) << std::setprecision(precision);
636 sout << 0. << std::endl;
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)"
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"
647 sout <<
" Fakes are calculated with the selections listed above" << std::endl;
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;
655 sout <<
" " << std::endl;
658 sout << std::setw(
width) << std::setprecision(precision);
660 sout << std::setw(
width) << std::setprecision(precision);
662 sout << std::setw(
width) << std::setprecision(precision);
664 sout << std::setw(
width) << std::setprecision(precision);
666 sout << std::setw(
width) << std::setprecision(precision);
668 sout << std::setw(
width) << std::setprecision(precision);
675 return StatusCode::SUCCESS;