File indexing completed on 2024-04-06 12:25:38
0001
0002
0003
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
0077 std::vector<DetId> m3x3aroundDC = EcalClusterTools::matrixDetId(topology_, id, 1);
0078 if (m3x3aroundDC.size() < 9) {
0079 *acceptFlag = false;
0080 return 0;
0081 }
0082
0083
0084
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) {
0109 mx_.rEn[cellIndex] = 0;
0110 cellIndex++;
0111 } else {
0112 *acceptFlag = false;
0113 return 0.;
0114 }
0115 } else {
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
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
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>;