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;
 
  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);