15 #include "GaudiKernel/SystemOfUnits.h"
29 declareInterface<ICaloMuonLikelihoodTool>(
this);
48 return StatusCode::SUCCESS;
55 ATH_MSG_DEBUG(
"in CaloMuonLikelihoodTool::retrieveHistograms()");
57 for (
int iFile = 0; iFile < 9; iFile++) {
66 return StatusCode::FAILURE;
71 return StatusCode::FAILURE;
75 ATH_MSG_FATAL(
"This exceeds the maximum allowed number of histograms");
76 return StatusCode::FAILURE;
79 int numKeysSignal{0}, numKeysBkg{0};
81 for (
int iHist = 0; iHist <
listOfKeys->GetSize(); iHist++) {
87 return StatusCode::FAILURE;
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) {
96 endOfKey =
histName.find(
"_bgnd", 0);
97 if (endOfKey == std::string::npos) {
99 return StatusCode::FAILURE;
102 const std::string
key =
histName.substr(0, endOfKey);
103 ATH_MSG_DEBUG(
"Found histogram for " << (isSignal ?
"signal" :
"background") <<
" with key " <<
key);
106 m_TH1F_sig[iFile][numKeysSignal] = std::move(unique_hist);
110 m_TH1F_bkg[iFile][numKeysBkg] = std::move(unique_hist);
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;
123 return StatusCode::SUCCESS;
130 const double dR_CUT)
const {
132 const EventContext& ctx = Gaudi::Hive::currentContext();
133 if (trk && ClusContainer) {
134 double eta_trk = trk->
eta();
139 std::unique_ptr<Trk::CaloExtension> caloExt =
m_caloExtensionTool->caloExtension(ctx, *trk);
141 if (!caloExt)
return 0;
143 if (!caloEntryInterSec) {
147 double eta_trkAtCalo = caloEntryInterSec->
eta();
148 double phi_trkAtCalo = caloEntryInterSec->parameters()[
Trk::phi];
150 double LR =
getLHR(ClusContainer, eta_trk, p_trk, eta_trkAtCalo, phi_trkAtCalo, dR_CUT);
161 const double eta_trkAtCalo,
const double phi_trkAtCalo,
const double dR_CUT)
const {
162 const double dR_CUT2 = dR_CUT * dR_CUT;
164 double etot_em{0}, etot_hd{0}, etot{0};
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};
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;
200 etot_em += eemb0 + eemb1 + eemb2 + eemb3 + eeme0 + eeme1 + eeme2 + eeme3 + etileg1 + etileg2 + etileg3;
230 etot_hd += etileb0 + etileb1 + etileb2 + etilee0 + etilee1 + etilee2 + ehec0 + ehec1 + ehec2 + ehec3 + efcal0 + efcal1 + efcal2;
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;
236 ATH_MSG_VERBOSE(
"Values extracted from the CaloTopoClusters for a cone of: "
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]);
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]; }
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};
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;
280 eemb3_wrtGroup =
s[3] / etot_em;
281 eeme1_wrtGroup =
s[5] / etot_em;
282 eeme3_wrtGroup =
s[7] / etot_em;
286 etileb0_wrtGroup =
s[11] / etot_hd;
287 etileb2_wrtGroup =
s[13] / etot_hd;
291 std::map<std::string, double> vars;
293 vars[
"EoverEtrk"] = EoverEtrk;
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;
309 ATH_MSG_DEBUG(
"likelihood discriminant variables values: " << dR_CUT);
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);
314 double LR{0}, ProbS{1}, ProbB{1};
317 if (std::abs(eta_trk) < 1.4) {
324 }
else if (std::abs(eta_trk) >= 1.4 && std::abs(eta_trk) <= 1.6) {
331 }
else if (std::abs(eta_trk) > 1.6 && std::abs(eta_trk) < 2.5) {
340 if (iFile < 0)
return LR;
347 std::map<std::string, double>::const_iterator
it = vars.find(
m_TH1F_key[iFile][
i]);
350 <<
": m_TH1F_key[" << iFile <<
"][" <<
i <<
"(/" <<
m_numKeys[iFile] <<
")] = " <<
m_TH1F_key[iFile][
i] <<
", "
351 << ((
it != vars.end()) ?
"found" :
"not found"));
353 if (
it != vars.end()) {
354 int bin_sig =
m_TH1F_sig[iFile][
i]->GetXaxis()->FindFixBin(
it->second);
357 <<
": it->first = " <<
it->first <<
", it->second = " <<
it->second <<
", bin_sig = " << bin_sig
358 <<
", NbinsX = " <<
m_TH1F_sig[iFile][
i]->GetNbinsX());
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);
364 ProbS *= SbinContent;
374 int bin_bkg =
m_TH1F_bkg[iFile][
i]->GetXaxis()->FindFixBin(
it->second);
377 <<
": it->first = " <<
it->first <<
", it->second = " <<
it->second <<
", bin_bkg = " << bin_bkg
378 <<
", NbinsX = " <<
m_TH1F_bkg[iFile][
i]->GetNbinsX());
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);
384 ProbB *= BbinContent;
397 <<
"> is not found in the calculated variable list");
407 LR = ProbS / (ProbS + ProbB);