ATLAS Offline Software
Loading...
Searching...
No Matches
CaloMuonLikelihoodTool.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2021 CERN for the benefit of the ATLAS collaboration
3*/
4
6
7#include <TString.h> // for Form
8
9#include <cmath>
10#include <iostream>
11#include <map>
12#include <string>
13
15#include "GaudiKernel/SystemOfUnits.h"
17#include "TFile.h"
18#include "TH1F.h"
23
25// CaloMuonLikelihoodTool constructor
27CaloMuonLikelihoodTool::CaloMuonLikelihoodTool(const std::string& type, const std::string& name, const IInterface* parent) :
28 AthAlgTool(type, name, parent), m_fileNames() {
29 declareInterface<ICaloMuonLikelihoodTool>(this);
30}
31
33// CaloMuonLikelihoodTool::initialize
36 ATH_MSG_INFO("Initializing " << name());
38 m_fileNames = {PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.A0.root", m_calibRelease.value().c_str())),
39 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.A1.root", m_calibRelease.value().c_str())),
40 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.A2.root", m_calibRelease.value().c_str())),
41 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.B0.root", m_calibRelease.value().c_str())),
42 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.B1.root", m_calibRelease.value().c_str())),
43 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.B2.root", m_calibRelease.value().c_str())),
44 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.C0.root", m_calibRelease.value().c_str())),
45 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.C1.root", m_calibRelease.value().c_str())),
46 PathResolverFindCalibFile(Form("%s/CaloMuonLikelihood.PDF.C2.root", m_calibRelease.value().c_str()))};
48 return StatusCode::SUCCESS;
49}
50
52// CaloMuonLikelihoodTool::retrieveHistograms
55 ATH_MSG_DEBUG("in CaloMuonLikelihoodTool::retrieveHistograms()");
56
57 for (int iFile = 0; iFile < 9; iFile++) {
58 m_numKeys[iFile] = 0;
59 const std::string& fileName = m_fileNames[iFile];
60
61 // --- Retrieving root files and list of keys ---
62
63 std::unique_ptr<TFile> rootFile{TFile::Open(fileName.c_str(), "READ")};
64 if (!rootFile) {
65 ATH_MSG_FATAL("Could not retrieve root file: " << fileName);
66 return StatusCode::FAILURE;
67 }
68 TList* listOfKeys = rootFile->GetListOfKeys();
69 if (!listOfKeys) {
70 ATH_MSG_FATAL("Could not retrieve key list: " << fileName);
71 return StatusCode::FAILURE;
72 }
73 ATH_MSG_DEBUG(listOfKeys->GetSize() << "histogram keys found in " << fileName);
74 if (listOfKeys->GetSize() > 22) {
75 ATH_MSG_FATAL("This exceeds the maximum allowed number of histograms");
76 return StatusCode::FAILURE;
77 }
78
79 int numKeysSignal{0}, numKeysBkg{0};
80 // --- Retrieving individual histograms ---
81 for (int iHist = 0; iHist < listOfKeys->GetSize(); iHist++) {
82 const std::string histName = listOfKeys->At(iHist)->GetName();
83 TH1F* hist = nullptr;
84 rootFile->GetObject(histName.c_str(), hist);
85 if (!hist) {
86 ATH_MSG_ERROR("cannot retrieve hist " << histName);
87 return StatusCode::FAILURE;
88 }
89 std::unique_ptr<TH1F> unique_hist{hist};
90 unique_hist->SetDirectory(nullptr);
91 bool isSignal = false;
92 size_t endOfKey = histName.find("_signal", 0);
93 if (endOfKey != std::string::npos) {
94 isSignal = true;
95 } else {
96 endOfKey = histName.find("_bgnd", 0);
97 if (endOfKey == std::string::npos) {
98 ATH_MSG_ERROR("Histogram with name: " << histName << " does not following the naming convention.");
99 return StatusCode::FAILURE;
100 }
101 }
102 const std::string key = histName.substr(0, endOfKey);
103 ATH_MSG_DEBUG("Found histogram for " << (isSignal ? "signal" : "background") << " with key " << key);
104 if (isSignal) {
105 m_TH1F_key[iFile][numKeysSignal] = key;
106 m_TH1F_sig[iFile][numKeysSignal] = std::move(unique_hist);
107 ATH_MSG_VERBOSE(m_TH1F_sig[iFile][numKeysSignal]->GetNbinsX());
108 ++numKeysSignal;
109 } else {
110 m_TH1F_bkg[iFile][numKeysBkg] = std::move(unique_hist);
111 ++numKeysBkg;
112 }
113 }
114
115 if (numKeysSignal != numKeysBkg) {
116 ATH_MSG_ERROR("The number of background histograms should be equal to the number of signal histograms.");
117 return StatusCode::FAILURE;
118 }
119 m_numKeys[iFile] = numKeysSignal;
120
121 ATH_MSG_DEBUG("Retrieved histograms from " << fileName);
122 }
123 return StatusCode::SUCCESS;
124}
125
127// CaloMuonLikelihoodTool::getLHR
130 const double dR_CUT) const {
131 ATH_MSG_DEBUG("in CaloMuonLikelihoodTool::getLHR()");
132 const EventContext& ctx = Gaudi::Hive::currentContext();
133 if (trk && ClusContainer) {
134 double eta_trk = trk->eta();
135 double p_trk = 0;
136 double qOverP = trk->qOverP();
137 if (qOverP != 0) p_trk = std::abs(1 / qOverP);
138
139 std::unique_ptr<Trk::CaloExtension> caloExt = m_caloExtensionTool->caloExtension(ctx, *trk);
140
141 if (!caloExt) return 0;
142 const Trk::TrackParameters* caloEntryInterSec = caloExt->caloEntryLayerIntersection();
143 if (!caloEntryInterSec) {
144 ATH_MSG_WARNING("getLHR() - caloEntryLayerIntersection is nullptr");
145 return 0;
146 }
147 double eta_trkAtCalo = caloEntryInterSec->eta();
148 double phi_trkAtCalo = caloEntryInterSec->parameters()[Trk::phi];
149
150 double LR = getLHR(ClusContainer, eta_trk, p_trk, eta_trkAtCalo, phi_trkAtCalo, dR_CUT);
151 return LR;
152 }
153
154 return 0;
155}
156
158// CaloMuonLikelihoodTool::getLHR
160double CaloMuonLikelihoodTool::getLHR(const xAOD::CaloClusterContainer* ClusCollection, const double eta_trk, const double p_trk,
161 const double eta_trkAtCalo, const double phi_trkAtCalo, const double dR_CUT) const {
162 const double dR_CUT2 = dR_CUT * dR_CUT;
163
164 double etot_em{0}, etot_hd{0}, etot{0};
165
166 double s[24] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
167
168 ATH_MSG_DEBUG("Loop over the CaloTopClusters...");
169 for (const xAOD::CaloCluster* theClus : *ClusCollection) {
170 double dphi = std::abs(theClus->phi() - phi_trkAtCalo);
171 if (dphi > M_PI) dphi = 2 * M_PI - dphi;
172 const double deta = std::abs(theClus->eta() - eta_trkAtCalo);
173 const double rij2 = deta * deta + dphi * dphi;
174 if (rij2 > dR_CUT2) continue;
175
176 double eemb0 = theClus->eSample(CaloSampling::PreSamplerB);
177 double eemb1 = theClus->eSample(CaloSampling::EMB1);
178 double eemb2 = theClus->eSample(CaloSampling::EMB2);
179 double eemb3 = theClus->eSample(CaloSampling::EMB3);
180 double eeme0 = theClus->eSample(CaloSampling::PreSamplerE);
181 double eeme1 = theClus->eSample(CaloSampling::EME1);
182 double eeme2 = theClus->eSample(CaloSampling::EME2);
183 double eeme3 = theClus->eSample(CaloSampling::EME3);
184 double etileg1 = theClus->eSample(CaloSampling::TileGap1);
185 double etileg2 = theClus->eSample(CaloSampling::TileGap2);
186 double etileg3 = theClus->eSample(CaloSampling::TileGap3);
187
188 s[0] += eemb0;
189 s[1] += eemb1;
190 s[2] += eemb2;
191 s[3] += eemb3;
192 s[4] += eeme0;
193 s[5] += eeme1;
194 s[6] += eeme2;
195 s[7] += eeme3;
196 s[8] += etileg1;
197 s[9] += etileg2;
198 s[10] += etileg3;
199
200 etot_em += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3;
201
202 double etileb0 = theClus->eSample(CaloSampling::TileBar0);
203 double etileb1 = theClus->eSample(CaloSampling::TileBar1);
204 double etileb2 = theClus->eSample(CaloSampling::TileBar2);
205 double etilee0 = theClus->eSample(CaloSampling::TileExt0);
206 double etilee1 = theClus->eSample(CaloSampling::TileExt1);
207 double etilee2 = theClus->eSample(CaloSampling::TileExt2);
208 double ehec0 = theClus->eSample(CaloSampling::HEC0);
209 double ehec1 = theClus->eSample(CaloSampling::HEC1);
210 double ehec2 = theClus->eSample(CaloSampling::HEC2);
211 double ehec3 = theClus->eSample(CaloSampling::HEC3);
212 double efcal0 = theClus->eSample(CaloSampling::FCAL0);
213 double efcal1 = theClus->eSample(CaloSampling::FCAL1);
214 double efcal2 = theClus->eSample(CaloSampling::FCAL2);
215
216 s[11] += etileb0;
217 s[12] += etileb1;
218 s[13] += etileb2;
219 s[14] += etilee0;
220 s[15] += etilee1;
221 s[16] += etilee2;
222 s[17] += ehec0;
223 s[18] += ehec1;
224 s[19] += ehec2;
225 s[20] += ehec3;
226 s[21] += efcal0;
227 s[22] += efcal1;
228 s[23] += efcal2;
229
230 etot_hd += etileb0 + etileb1 + etileb2 + etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
231
232 etot += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3 + etileb0 + etileb1 + etileb2 +
233 etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
234 }
235
236 ATH_MSG_VERBOSE("Values extracted from the CaloTopoClusters for a cone of: "
237 << dR_CUT << "\n"
238 << " - Energy in Calorimeter Total: " << etot << " EM: " << etot_em << " HAD: " << etot_hd
239 << "\n eemb0: " << s[0] << "\n eemb1: " << s[1] << "\n eemb2: " << s[2] << "\n eemb3: " << s[3]
240 << "\n eeme0: " << s[4] << "\n eeme1: " << s[5] << "\n eeme2: " << s[6] << "\n eeme3: " << s[7]
241 << "\n etileg1: " << s[8] << "\n etileg2: " << s[9] << "\n etileg3: " << s[10] << "\n etileb0: " << s[11]
242 << "\n etileb1: " << s[12] << "\n etileb2: " << s[13] << "\n etilee0: " << s[14] << "\n etilee1: " << s[15]
243 << "\n etilee2: " << s[16] << "\n ehec0: " << s[17] << "\n ehec1: " << s[18] << "\n ehec2: " << s[19]
244 << "\n ehec3: " << s[20] << "\n efcal0: " << s[21] << "\n efcal1: " << s[22]
245 << "\n efcal2: " << s[23]);
246
247 for (int i = 0; i < 24; ++i) s[i] /= Gaudi::Units::GeV;
248 etot_em /= Gaudi::Units::GeV;
249 etot_hd /= Gaudi::Units::GeV;
250 etot /= Gaudi::Units::GeV;
251
252 double tmp_mx(-999), tmp_mxH(-999), tmp_mxE(-999);
253 for (int i = 0; i < 24; ++i) {
254 if (s[i] > tmp_mx) { tmp_mx = s[i]; }
255 if (i < 11 && s[i] > tmp_mxE ) { tmp_mxE = s[i]; }
256 if (i >= 11 && s[i] > tmp_mxH ) { tmp_mxH = s[i]; }
257 }
258
259 double emFr{999}, mxFr{999}, mxEM{999}, mxHad{999}, EoverEtrk{999}, eemb1_wrtTotal{999}, eemb2_wrtTotal{999}, eemb3_wrtGroup{999},
260 etileb0_wrtGroup{999}, etileb1_wrtTotal{999}, etileb2_wrtGroup{999}, eeme1_wrtGroup{999}, eeme1_wrtTotal{999}, eeme2_wrtTotal{999},
261 eeme3_wrtGroup{999}, ehec0_wrtTotal{999};
262
263 if (etot > 0) {
264 mxFr = tmp_mx / etot;
265 mxEM = tmp_mxE / etot;
266 mxHad = tmp_mxH / etot;
267 emFr = etot_em / etot;
268 eemb1_wrtTotal = s[1] / etot;
269 eemb2_wrtTotal = s[2] / etot;
270 etileb1_wrtTotal = s[12] / etot;
271 eeme1_wrtTotal = s[5] / etot;
272 eeme2_wrtTotal = s[6] / etot;
273 ehec0_wrtTotal = s[17] / etot;
274 }
275
276 // p_trk in MeV
277 if (p_trk) EoverEtrk = etot / (p_trk / Gaudi::Units::GeV);
278
279 if (etot_em > 0) {
280 eemb3_wrtGroup = s[3] / etot_em;
281 eeme1_wrtGroup = s[5] / etot_em;
282 eeme3_wrtGroup = s[7] / etot_em;
283 }
284
285 if (etot_hd > 0) {
286 etileb0_wrtGroup = s[11] / etot_hd;
287 etileb2_wrtGroup = s[13] / etot_hd;
288 }
289
290 // likelihood discriminant variables
291 std::map<std::string, double> vars;
292 vars["emFr"] = emFr;
293 vars["EoverEtrk"] = EoverEtrk;
294 vars["mxFr"] = mxFr;
295 vars["mxEM"] = mxEM;
296 vars["mxHad"] = mxHad;
297 vars["eemb1_wrtTotal"] = eemb1_wrtTotal;
298 vars["eemb2_wrtTotal"] = eemb2_wrtTotal;
299 vars["eemb3_wrtGroup"] = eemb3_wrtGroup;
300 vars["etileb0_wrtGroup"] = etileb0_wrtGroup;
301 vars["etileb1_wrtTotal"] = etileb1_wrtTotal;
302 vars["etileb2_wrtGroup"] = etileb2_wrtGroup;
303 vars["eeme1_wrtGroup"] = eeme1_wrtGroup;
304 vars["eeme1_wrtTotal"] = eeme1_wrtTotal;
305 vars["eeme2_wrtTotal"] = eeme2_wrtTotal;
306 vars["eeme3_wrtGroup"] = eeme3_wrtGroup;
307 vars["ehec0_wrtTotal"] = ehec0_wrtTotal;
308
309 ATH_MSG_DEBUG("likelihood discriminant variables values: " << dR_CUT);
310 if (msgLevel(MSG::DEBUG)) {
311 std::map<std::string, double>::const_iterator iter;
312 for (iter = vars.begin(); iter != vars.end(); ++iter) ATH_MSG_DEBUG(" - " << iter->first << ": " << iter->second);
313 }
314 double LR{0}, ProbS{1}, ProbB{1};
315
316 int iFile = -1;
317 if (std::abs(eta_trk) < 1.4) {
318 if (p_trk / Gaudi::Units::GeV < 11.)
319 iFile = 0;
320 else if (p_trk / Gaudi::Units::GeV < 51.)
321 iFile = 1;
322 else
323 iFile = 2;
324 } else if (std::abs(eta_trk) >= 1.4 && std::abs(eta_trk) <= 1.6) {
325 if (p_trk / Gaudi::Units::GeV < 11.)
326 iFile = 3;
327 else if (p_trk / Gaudi::Units::GeV < 51.)
328 iFile = 4;
329 else
330 iFile = 5;
331 } else if (std::abs(eta_trk) > 1.6 && std::abs(eta_trk) < 2.5) {
332 if (p_trk / Gaudi::Units::GeV < 11.)
333 iFile = 6;
334 else if (p_trk / Gaudi::Units::GeV < 51.)
335 iFile = 7;
336 else
337 iFile = 8;
338 }
339
340 if (iFile < 0) return LR;
341
342 if (m_numKeys[iFile] == 0) {
343 ProbS = 0;
344 ProbB = 1;
345 } else {
346 for (int i = 0; i < m_numKeys[iFile]; i++) {
347 std::map<std::string, double>::const_iterator it = vars.find(m_TH1F_key[iFile][i]);
348
349 ATH_MSG_VERBOSE("getLHR "
350 << ": m_TH1F_key[" << iFile << "][" << i << "(/" << m_numKeys[iFile] << ")] = " << m_TH1F_key[iFile][i] << ", "
351 << ((it != vars.end()) ? "found" : "not found"));
352
353 if (it != vars.end()) {
354 int bin_sig = m_TH1F_sig[iFile][i]->GetXaxis()->FindFixBin(it->second);
355
356 ATH_MSG_VERBOSE("getLHR sig "
357 << ": it->first = " << it->first << ", it->second = " << it->second << ", bin_sig = " << bin_sig
358 << ", NbinsX = " << m_TH1F_sig[iFile][i]->GetNbinsX());
359
360 if (bin_sig >= 1 && bin_sig <= m_TH1F_sig[iFile][i]->GetNbinsX()) {
361 double SbinContent = m_TH1F_sig[iFile][i]->GetBinContent(bin_sig);
362 ATH_MSG_DEBUG("m_TH1F_sig Bin Content for " << it->first << ": " << SbinContent);
363 if (SbinContent)
364 ProbS *= SbinContent;
365 else {
366 ProbS *= 0.0000001;
367 ATH_MSG_DEBUG("BinContent in m_TH1F_sig was 0!" << ProbS);
368 }
369 } else {
370 ATH_MSG_DEBUG("Sbin:" << bin_sig << " and NBinMax: " << m_TH1F_sig[iFile][i]->GetNbinsX());
371 ProbS *= 1;
372 }
373 ATH_MSG_DEBUG("Temp ProbS : " << ProbS);
374 int bin_bkg = m_TH1F_bkg[iFile][i]->GetXaxis()->FindFixBin(it->second);
375
376 ATH_MSG_VERBOSE("getLHR bkg "
377 << ": it->first = " << it->first << ", it->second = " << it->second << ", bin_bkg = " << bin_bkg
378 << ", NbinsX = " << m_TH1F_bkg[iFile][i]->GetNbinsX());
379
380 if (bin_bkg >= 1 && bin_bkg <= m_TH1F_bkg[iFile][i]->GetNbinsX()) {
381 double BbinContent = m_TH1F_bkg[iFile][i]->GetBinContent(bin_bkg);
382 ATH_MSG_VERBOSE("m_TH1F_bkg Bin Content for " << it->first << ": " << BbinContent);
383 if (BbinContent)
384 ProbB *= BbinContent;
385 else {
386 ProbB *= 0.0000001;
387 ATH_MSG_DEBUG("BinContent in m_TH1F_bkg was 0! " << ProbB);
388 }
389 } else {
390 ATH_MSG_DEBUG("Bbin:" << bin_bkg << " and NBinMax: " << m_TH1F_bkg[iFile][i]->GetNbinsX());
391 ProbB *= 1;
392 }
393 ATH_MSG_DEBUG(" Temp ProbB: " << ProbB);
394 } else {
395 if (m_cnt_warn < 10) {
396 ATH_MSG_WARNING("Histogram variable <" << m_TH1F_key[iFile][i] << "> in file <" << m_fileNames[iFile]
397 << "> is not found in the calculated variable list");
398 m_cnt_warn++;
399 }
400 }
401 }
402 }
403
404 if (ProbS == 0)
405 LR = 0;
406 else
407 LR = ProbS / (ProbS + ProbB);
408
409 ATH_MSG_DEBUG("ProbS : " << ProbS);
410 ATH_MSG_DEBUG("ProbB : " << ProbB);
411 ATH_MSG_DEBUG("LR : " << LR);
412 return LR;
413}
#define M_PI
#define ATH_CHECK
Evaluate an expression and check for errors.
#define ATH_MSG_ERROR(x)
#define ATH_MSG_FATAL(x)
#define ATH_MSG_INFO(x)
#define ATH_MSG_VERBOSE(x)
#define ATH_MSG_WARNING(x)
#define ATH_MSG_DEBUG(x)
std::string PathResolverFindCalibFile(const std::string &logical_file_name)
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Gaudi::Property< std::string > m_calibRelease
ToolHandle< Trk::IParticleCaloExtensionTool > m_caloExtensionTool
std::vector< std::string > m_fileNames
std::unique_ptr< const TH1F > m_TH1F_sig[9][11]
std::unique_ptr< const TH1F > m_TH1F_bkg[9][11]
CaloMuonLikelihoodTool(const std::string &type, const std::string &name, const IInterface *parent)
virtual StatusCode initialize()
double getLHR(const xAOD::TrackParticle *trk, const xAOD::CaloClusterContainer *ClusCollection=nullptr, const double dR_CUT=0.3) const
double eta() const
Access method for pseudorapidity - from momentum.
float qOverP() const
Returns the parameter.
virtual double eta() const override final
The pseudorapidity ( ) of the particle.
@ phi
Definition ParamDefs.h:75
ParametersBase< TrackParametersDim, Charged > TrackParameters
CaloCluster_v1 CaloCluster
Define the latest version of the calorimeter cluster class.
TrackParticle_v1 TrackParticle
Reference the current persistent version:
CaloClusterContainer_v1 CaloClusterContainer
Define the latest version of the calorimeter cluster container.