EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc

0001 //  Original Authors:  S. Taroni, N. Marinelli
0002 //  University of Notre Dame - US
0003 //  Created:
0004 //
0005 //
0006 //
0007
0008 #include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h"
0009 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h"  // can I use a egammatools here?
0010 #include "FWCore/ParameterSet/interface/FileInPath.h"
0011 #include <iostream>
0012
0013 #include <iostream>
0014 #include <memory>
0015
0016 #include <fstream>
0017 #include <ostream>
0018 #include <string>
0019
0020 template <>
0021 void EcalDeadChannelRecoveryBDTG<EBDetId>::addVariables(TMVA::Reader *reader) {
0022   for (int i = 0; i < 9; ++i) {
0023     if (i == 4)
0024       continue;
0025     reader->AddVariable("E" + std::to_string(i + 1) + "/(E1+E2+E3+E4+E6+E7+E8+E9)", &(mx_.rEn[i]));
0026   }
0027   reader->AddVariable("E1+E2+E3+E4+E6+E7+E8+E9", &(mx_.sumE8));
0028   for (int i = 0; i < 9; ++i)
0029     reader->AddVariable("iEta" + std::to_string(i + 1), &(mx_.ieta[i]));
0030   for (int i = 0; i < 9; ++i)
0031     reader->AddVariable("iPhi" + std::to_string(i + 1), &(mx_.iphi[i]));
0032 }
0033 template <>
0034 void EcalDeadChannelRecoveryBDTG<EBDetId>::loadFile() {
0035   readerNoCrack = std::make_unique<TMVA::Reader>("!Color:!Silent");
0036   readerCrack = std::make_unique<TMVA::Reader>("!Color:!Silent");
0037
0038   addVariables(readerNoCrack.get());
0039   addVariables(readerCrack.get());
0040
0041   reco::details::loadTMVAWeights(readerNoCrack.get(), "BDTG", bdtWeightFileNoCracks_.fullPath());
0042   reco::details::loadTMVAWeights(readerCrack.get(), "BDTG", bdtWeightFileCracks_.fullPath());
0043 }
0044
0045 template <typename T>
0046 EcalDeadChannelRecoveryBDTG<T>::EcalDeadChannelRecoveryBDTG() {}
0047
0048 template <typename T>
0049 EcalDeadChannelRecoveryBDTG<T>::~EcalDeadChannelRecoveryBDTG() {}
0050
0051 template <>
0052 void EcalDeadChannelRecoveryBDTG<EBDetId>::setParameters(const edm::ParameterSet &ps) {
0053   bdtWeightFileNoCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileNoCracks");
0054   bdtWeightFileCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileCracks");
0055
0056   loadFile();
0057 }
0058
0059 template <>
0060 void EcalDeadChannelRecoveryBDTG<EEDetId>::setParameters(const edm::ParameterSet &ps) {}
0061
0062 template <>
0063 double EcalDeadChannelRecoveryBDTG<EEDetId>::recover(
0064     const EEDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) {
0065   return 0;
0066 }
0067
0068 template <>
0069 double EcalDeadChannelRecoveryBDTG<EBDetId>::recover(
0070     const EBDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) {
0071   bool isCrack = false;
0072   int cellIndex = 0.;
0073   double neighTotEn = 0.;
0074   float val = 0.;
0075
0076   //find the matrix around id
0077   std::vector<DetId> m3x3aroundDC = EcalClusterTools::matrixDetId(topology_, id, 1);
0078   if (m3x3aroundDC.size() < 9) {
0079     *acceptFlag = false;
0080     return 0;
0081   }
0082
0083   //  Loop over all cells in the vector "NxNaroundDC", and for each cell find it's energy
0084   //  (from the EcalRecHits collection).
0085   for (auto const &theCells : m3x3aroundDC) {
0086     EBDetId cell = EBDetId(theCells);
0087     if (cell == id) {
0088       int iEtaCentral = std::abs(cell.ieta());
0089       int iPhiCentral = cell.iphi();
0090
0091       if (iEtaCentral < 2 || std::abs(iEtaCentral - 25) < 2 || std::abs(iEtaCentral - 45) < 2 ||
0092           std::abs(iEtaCentral - 65) < 2 || iEtaCentral > 83 || (int(iPhiCentral + 0.5) % 20 == 0))
0093         isCrack = true;
0094     }
0095     if (!cell.null()) {
0096       EcalRecHitCollection::const_iterator goS_it = hit_collection.find(cell);
0097       if (goS_it != hit_collection.end() && cell != id) {
0098         if (goS_it->energy() < single8Cut) {
0099           *acceptFlag = false;
0100           return 0.;
0101         } else {
0102           neighTotEn += goS_it->energy();
0103           mx_.rEn[cellIndex] = goS_it->energy();
0104           mx_.iphi[cellIndex] = cell.iphi();
0105           mx_.ieta[cellIndex] = cell.ieta();
0106           cellIndex++;
0107         }
0108       } else if (cell == id) {  // the cell is the central one
0109         mx_.rEn[cellIndex] = 0;
0110         cellIndex++;
0111       } else {  //goS_it is not in the rechitcollection
0112         *acceptFlag = false;
0113         return 0.;
0114       }
0115     } else {  //cell is null
0116       *acceptFlag = false;
0117       return 0.;
0118     }
0119   }
0120   if (cellIndex > 0 && neighTotEn >= single8Cut * 8. && neighTotEn >= sum8Cut) {
0121     bool allneighs = true;
0122     mx_.sumE8 = neighTotEn;
0123     for (unsigned int icell = 0; icell < 9; icell++) {
0124       if (mx_.rEn[icell] < single8Cut && icell != 4) {
0125         allneighs = false;
0126       }
0127       mx_.rEn[icell] = mx_.rEn[icell] / neighTotEn;
0128     }
0129     if (allneighs == true) {
0130       // evaluate the regression
0131       if (isCrack) {
0132         val = exp((readerCrack->EvaluateRegression("BDTG"))[0]);
0133       } else {
0134         val = exp((readerNoCrack->EvaluateRegression("BDTG"))[0]);
0135       }
0136
0137       *acceptFlag = true;
0138       //return the estimated energy
0139       return val;
0140     } else {
0141       *acceptFlag = false;
0142       return 0;
0143     }
0144   } else {
0145     *acceptFlag = false;
0146     return 0;
0147   }
0148 }
0149
0150 template class EcalDeadChannelRecoveryBDTG<EBDetId>;
0151 template class EcalDeadChannelRecoveryBDTG<EEDetId>;  //not used.