Skip to content

Commit

Permalink
migrate EcalDeadChannelRecoveryBDTG to GBRForest
Browse files Browse the repository at this point in the history
  • Loading branch information
@guitargeek
guitargeek committed Sep 20, 2019
1 parent d920ba0 commit 49be72e
Show file tree
Hide file tree
Showing 5 changed files with 50 additions and 78 deletions.
2 changes: 1 addition & 1 deletion RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ class EcalDeadChannelRecoveryAlgos {
std::string algo,
double single8Cut,
double sum8Cut,
bool *accFlag);
bool &accFlag);

private:
EcalDeadChannelRecoveryBDTG<DetIdT> bdtg_;
Expand Down
24 changes: 4 additions & 20 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,43 +9,27 @@
#include "Geometry/CaloTopology/interface/CaloTopology.h"
#include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"

#include "CommonTools/MVAUtils/interface/TMVAZipReader.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "TMVA/Reader.h"
#include "CommonTools/MVAUtils/interface/GBRForestTools.h"

#include <string>
#include <memory>

template <typename DetIdT>
class EcalDeadChannelRecoveryBDTG {
public:
EcalDeadChannelRecoveryBDTG();
~EcalDeadChannelRecoveryBDTG();

void setParameters(const edm::ParameterSet &ps);
void setCaloTopology(const CaloTopology *topo) { topology_ = topo; }

double recover(
const DetIdT id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag);

void loadFile();
void addVariables(TMVA::Reader *reader);
const DetIdT id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool &acceptFlag);

private:
const CaloTopology *topology_;
struct XtalMatrix {
std::array<float, 9> rEn, ieta, iphi;
float sumE8;
};

XtalMatrix mx_;

edm::FileInPath bdtWeightFileNoCracks_;
edm::FileInPath bdtWeightFileCracks_;

std::unique_ptr<TMVA::Reader> readerNoCrack;
std::unique_ptr<TMVA::Reader> readerCrack;
std::unique_ptr<const GBRForest> gbrForestNoCrack_;
std::unique_ptr<const GBRForest> gbrForestCrack_;
};

#endif
8 changes: 4 additions & 4 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,17 +29,17 @@ float EcalDeadChannelRecoveryAlgos<T>::correct(const T id,
std::string algo,
double single8Cut,
double sum8Cut,
bool *acceptFlag) {
bool &acceptFlag) {
// recover as single dead channel
double newEnergy = 0.0;
if (algo == "BDTG") {
*acceptFlag = false;
acceptFlag = false;
newEnergy = this->bdtg_.recover(id, hit_collection, single8Cut, sum8Cut, acceptFlag); //ADD here
if (newEnergy > 0.)
*acceptFlag = true; //bdtg set to 0 if there is more than one channel in the matrix that is not reponding
acceptFlag = true; //bdtg set to 0 if there is more than one channel in the matrix that is not reponding
} else {
edm::LogError("EcalDeadChannelRecoveryAlgos") << "Invalid algorithm for dead channel recovery.";
*acceptFlag = false;
acceptFlag = false;
}

return newEnergy;
Expand Down
90 changes: 39 additions & 51 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,57 +15,45 @@
#include <string>
#include <fstream>

template <>
void EcalDeadChannelRecoveryBDTG<EBDetId>::addVariables(TMVA::Reader *reader) {
for (int i = 0; i < 9; ++i) {
if (i == 4)
continue;
reader->AddVariable("E" + std::to_string(i + 1) + "/(E1+E2+E3+E4+E6+E7+E8+E9)", &(mx_.rEn[i]));
namespace {

struct XtalMatrix {
float rEn[9];
float sumE8;
float ieta[9];
float iphi[9];
};

std::vector<float> xtalMatrixToVector(XtalMatrix const &xtalMatrix) {
std::vector<float> v(xtalMatrix.rEn, xtalMatrix.rEn + 4);
v.insert(v.end(), xtalMatrix.rEn + 5, xtalMatrix.rEn + 9);
v.push_back(xtalMatrix.sumE8);
v.insert(v.end(), xtalMatrix.ieta, xtalMatrix.ieta + 9);
v.insert(v.end(), xtalMatrix.iphi, xtalMatrix.iphi + 9);
return v;
}
reader->AddVariable("E1+E2+E3+E4+E6+E7+E8+E9", &(mx_.sumE8));
for (int i = 0; i < 9; ++i)
reader->AddVariable("iEta" + std::to_string(i + 1), &(mx_.ieta[i]));
for (int i = 0; i < 9; ++i)
reader->AddVariable("iPhi" + std::to_string(i + 1), &(mx_.iphi[i]));
}
template <>
void EcalDeadChannelRecoveryBDTG<EBDetId>::loadFile() {
readerNoCrack = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:!Silent"));
readerCrack = std::unique_ptr<TMVA::Reader>(new TMVA::Reader("!Color:!Silent"));

addVariables(readerNoCrack.get());
addVariables(readerCrack.get());

reco::details::loadTMVAWeights(readerNoCrack.get(), "BDTG", bdtWeightFileNoCracks_.fullPath());
reco::details::loadTMVAWeights(readerCrack.get(), "BDTG", bdtWeightFileCracks_.fullPath());
}

template <typename T>
EcalDeadChannelRecoveryBDTG<T>::EcalDeadChannelRecoveryBDTG() {}

template <typename T>
EcalDeadChannelRecoveryBDTG<T>::~EcalDeadChannelRecoveryBDTG() {}
} // namespace

template <>
void EcalDeadChannelRecoveryBDTG<EBDetId>::setParameters(const edm::ParameterSet &ps) {
bdtWeightFileNoCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileNoCracks");
bdtWeightFileCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileCracks");

loadFile();
gbrForestNoCrack_ = createGBRForest(ps.getParameter<edm::FileInPath>("bdtWeightFileNoCracks"));
gbrForestCrack_ = createGBRForest(ps.getParameter<edm::FileInPath>("bdtWeightFileCracks"));
}

template <>
void EcalDeadChannelRecoveryBDTG<EEDetId>::setParameters(const edm::ParameterSet &ps) {}

template <>
double EcalDeadChannelRecoveryBDTG<EEDetId>::recover(
const EEDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) {
const EEDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool &acceptFlag) {
return 0;
}

template <>
double EcalDeadChannelRecoveryBDTG<EBDetId>::recover(
const EBDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) {
const EBDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool &acceptFlag) {
XtalMatrix mx;

bool isCrack = false;
int cellIndex = 0.;
double neighTotEn = 0.;
Expand All @@ -74,7 +62,7 @@ double EcalDeadChannelRecoveryBDTG<EBDetId>::recover(
//find the matrix around id
std::vector<DetId> m3x3aroundDC = EcalClusterTools::matrixDetId(topology_, id, 1);
if (m3x3aroundDC.size() < 9) {
*acceptFlag = false;
acceptFlag = false;
return 0;
}

Expand All @@ -94,53 +82,53 @@ double EcalDeadChannelRecoveryBDTG<EBDetId>::recover(
EcalRecHitCollection::const_iterator goS_it = hit_collection.find(cell);
if (goS_it != hit_collection.end() && cell != id) {
if (goS_it->energy() < single8Cut) {
*acceptFlag = false;
acceptFlag = false;
return 0.;
} else {
neighTotEn += goS_it->energy();
mx_.rEn[cellIndex] = goS_it->energy();
mx_.iphi[cellIndex] = cell.iphi();
mx_.ieta[cellIndex] = cell.ieta();
mx.rEn[cellIndex] = goS_it->energy();
mx.iphi[cellIndex] = cell.iphi();
mx.ieta[cellIndex] = cell.ieta();
cellIndex++;
}
} else if (cell == id) { // the cell is the central one
mx_.rEn[cellIndex] = 0;
mx.rEn[cellIndex] = 0;
cellIndex++;
} else { //goS_it is not in the rechitcollection
*acceptFlag = false;
acceptFlag = false;
return 0.;
}
} else { //cell is null
*acceptFlag = false;
acceptFlag = false;
return 0.;
}
}
if (cellIndex > 0 && neighTotEn >= single8Cut * 8. && neighTotEn >= sum8Cut) {
bool allneighs = true;
mx_.sumE8 = neighTotEn;
mx.sumE8 = neighTotEn;
for (unsigned int icell = 0; icell < 9; icell++) {
if (mx_.rEn[icell] < single8Cut && icell != 4) {
if (mx.rEn[icell] < single8Cut && icell != 4) {
allneighs = false;
}
mx_.rEn[icell] = mx_.rEn[icell] / neighTotEn;
mx.rEn[icell] = mx.rEn[icell] / neighTotEn;
}
if (allneighs == true) {
// evaluate the regression
if (isCrack) {
val = exp((readerCrack->EvaluateRegression("BDTG"))[0]);
val = exp(gbrForestCrack_->GetResponse(xtalMatrixToVector(mx).data()));
} else {
val = exp((readerNoCrack->EvaluateRegression("BDTG"))[0]);
val = exp(gbrForestNoCrack_->GetResponse(xtalMatrixToVector(mx).data()));
}

*acceptFlag = true;
acceptFlag = true;
//return the estimated energy
return val;
} else {
*acceptFlag = false;
acceptFlag = false;
return 0;
}
} else {
*acceptFlag = false;
acceptFlag = false;
return 0;
}
}
Expand Down
4 changes: 2 additions & 2 deletions RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -129,7 +129,7 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,
// channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
bool AcceptRecHit = true;
float ebEn = ebDeadChannelCorrector.correct(
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, AcceptRecHit);
EcalRecHit hit(detId, ebEn, 0., EcalRecHit::kDead);

if (hit.energy() != 0 and AcceptRecHit == true) {
Expand All @@ -147,7 +147,7 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,
// channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
bool AcceptRecHit = true;
float eeEn = eeDeadChannelCorrector.correct(
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, AcceptRecHit);
EcalRecHit hit(detId, eeEn, 0., EcalRecHit::kDead);
if (hit.energy() != 0 and AcceptRecHit == true) {
hit.setFlag(EcalRecHit::kNeighboursRecovered);
Expand Down

0 comments on commit 49be72e

Please sign in to comment.