58 ATH_MSG_DEBUG(
"in CaloMuonLikelihoodTool::retrieveHistograms()");
60 for (
int iFile = 0; iFile < 9; iFile++) {
66 std::unique_ptr<TFile> rootFile{TFile::Open(fileName.c_str(),
"READ")};
69 return StatusCode::FAILURE;
71 TList* listOfKeys = rootFile->GetListOfKeys();
74 return StatusCode::FAILURE;
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;
82 int numKeysSignal{0}, numKeysBkg{0};
84 for (
int iHist = 0; iHist < listOfKeys->GetSize(); iHist++) {
85 const std::string histName = listOfKeys->At(iHist)->GetName();
87 rootFile->GetObject(histName.c_str(), hist);
90 return StatusCode::FAILURE;
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) {
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;
105 const std::string key = histName.substr(0, endOfKey);
106 ATH_MSG_DEBUG(
"Found histogram for " << (isSignal ?
"signal" :
"background") <<
" with key " << key);
109 m_TH1F_sig[iFile][numKeysSignal] = std::move(unique_hist);
113 m_TH1F_bkg[iFile][numKeysBkg] = std::move(unique_hist);
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;
126 return StatusCode::SUCCESS;
164 const double eta_trkAtCalo,
const double phi_trkAtCalo,
const double dR_CUT)
const {
165 const double dR_CUT2 = dR_CUT * dR_CUT;
167 double etot_em{0}, etot_hd{0}, etot{0};
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};
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;
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);
203 etot_em += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3;
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);
233 etot_hd += etileb0 + etileb1 + etileb2 + etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
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;
239 ATH_MSG_VERBOSE(
"Values extracted from the CaloTopoClusters for a cone of: "
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]);
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;
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]; }
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};
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;
280 if (p_trk) EoverEtrk = etot / (p_trk / Gaudi::Units::GeV);
283 eemb3_wrtGroup = s[3] / etot_em;
284 eeme1_wrtGroup = s[5] / etot_em;
285 eeme3_wrtGroup = s[7] / etot_em;
289 etileb0_wrtGroup = s[11] / etot_hd;
290 etileb2_wrtGroup = s[13] / etot_hd;
294 std::map<std::string, double> vars;
296 vars[
"EoverEtrk"] = EoverEtrk;
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;
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);
317 double LR{0}, ProbS{1}, ProbB{1};
320 if (std::abs(eta_trk) < 1.4) {
321 if (p_trk / Gaudi::Units::GeV < 11.)
323 else if (p_trk / Gaudi::Units::GeV < 51.)
327 }
else if (std::abs(eta_trk) >= 1.4 && std::abs(eta_trk) <= 1.6) {
328 if (p_trk / Gaudi::Units::GeV < 11.)
330 else if (p_trk / Gaudi::Units::GeV < 51.)
334 }
else if (std::abs(eta_trk) > 1.6 && std::abs(eta_trk) < 2.5) {
335 if (p_trk / Gaudi::Units::GeV < 11.)
337 else if (p_trk / Gaudi::Units::GeV < 51.)
343 if (iFile < 0)
return LR;
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]);
353 <<
": m_TH1F_key[" << iFile <<
"][" << i <<
"(/" <<
m_numKeys[iFile] <<
")] = " <<
m_TH1F_key[iFile][i] <<
", "
354 << ((it != vars.end()) ?
"found" :
"not found"));
356 if (it != vars.end()) {
357 int bin_sig =
m_TH1F_sig[iFile][i]->GetXaxis()->FindFixBin(it->second);
360 <<
": it->first = " << it->first <<
", it->second = " << it->second <<
", bin_sig = " << bin_sig
361 <<
", NbinsX = " <<
m_TH1F_sig[iFile][i]->GetNbinsX());
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);
367 ProbS *= SbinContent;
377 int bin_bkg =
m_TH1F_bkg[iFile][i]->GetXaxis()->FindFixBin(it->second);
380 <<
": it->first = " << it->first <<
", it->second = " << it->second <<
", bin_bkg = " << bin_bkg
381 <<
", NbinsX = " <<
m_TH1F_bkg[iFile][i]->GetNbinsX());
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);
387 ProbB *= BbinContent;
400 <<
"> is not found in the calculated variable list");
410 LR = ProbS / (ProbS + ProbB);