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