Skip to content

Commit

Permalink
Merge pull request #27175 from nancymarinelli/deadChannelRecoveryBDT_V1
Browse files Browse the repository at this point in the history
Ecal DPG - Deploy new software (based on BDT) to recover isolated dead  channels…
  • Loading branch information
cmsbuild committed Aug 2, 2019
2 parents 0c4a74d + edab6bc commit ffabcf0
Show file tree
Hide file tree
Showing 10 changed files with 268 additions and 26 deletions.
1 change: 1 addition & 0 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml
Expand Up @@ -3,6 +3,7 @@
<use name="Geometry/CaloGeometry"/>
<use name="Geometry/CaloTopology"/>
<use name="Geometry/Records"/>
<use name="CommonTools/MVAUtils"/>
<use name="root"/>
<use name="rootmath"/>
<use name="rootcore"/>
Expand Down
Expand Up @@ -7,16 +7,25 @@
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include <string>

#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryNN.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"

template <typename DetIdT>
class EcalDeadChannelRecoveryAlgos {
public:
void setParameters(const edm::ParameterSet &ps);
void setCaloTopology(const CaloTopology *topology);
EcalRecHit correct(
const DetIdT id, const EcalRecHitCollection &hit_collection, std::string algo, double Sum8Cut, bool *AccFlag);
float correct(const DetIdT id,
const EcalRecHitCollection &hit_collection,
std::string algo,
double single8Cut,
double sum8Cut,
bool *accFlag);

private:
EcalDeadChannelRecoveryNN<DetIdT> nn;
EcalDeadChannelRecoveryBDTG<DetIdT> bdtg_;
};
#endif // RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryAlgos_HH
#endif
@@ -0,0 +1,52 @@
#ifndef RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryBDTG_H
#define RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryBDTG_H

#include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"

#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"

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

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

#endif
Expand Up @@ -8,30 +8,41 @@
// Feb 14 2013: Implementation of the criterion to select the "correct"
// max. cont. crystal.
//
//modified by S.Taroni, N. Marinelli 11 June 2019

#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

template <typename T>
void EcalDeadChannelRecoveryAlgos<T>::setParameters(const edm::ParameterSet &ps) {
bdtg_.setParameters(ps);
}

template <typename T>
void EcalDeadChannelRecoveryAlgos<T>::setCaloTopology(const CaloTopology *topo) {
nn.setCaloTopology(topo);
bdtg_.setCaloTopology(topo);
}

template <typename T>
EcalRecHit EcalDeadChannelRecoveryAlgos<T>::correct(
const T id, const EcalRecHitCollection &hit_collection, std::string algo, double Sum8Cut, bool *AcceptFlag) {
float EcalDeadChannelRecoveryAlgos<T>::correct(const T id,
const EcalRecHitCollection &hit_collection,
std::string algo,
double single8Cut,
double sum8Cut,
bool *acceptFlag) {
// recover as single dead channel
double NewEnergy = 0.0;

if (algo == "NeuralNetworks") {
NewEnergy = this->nn.recover(id, hit_collection, Sum8Cut, AcceptFlag);
double newEnergy = 0.0;
if (algo == "BDTG") {
*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
} else {
edm::LogError("EcalDeadChannelRecoveryAlgos") << "Invalid algorithm for dead channel recovery.";
*AcceptFlag = false;
*acceptFlag = false;
}

uint32_t flag = 0;
return EcalRecHit(id, NewEnergy, 0, flag);
return newEnergy;
}

template class EcalDeadChannelRecoveryAlgos<EBDetId>;
Expand Down
@@ -0,0 +1,149 @@
// Original Authors: S. Taroni, N. Marinelli
// University of Notre Dame - US
// Created:
//
//
//

#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h"
#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" // can I use a egammatools here?
#include "FWCore/ParameterSet/interface/FileInPath.h"
#include <iostream>

#include <iostream>
#include <ostream>
#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]));
}
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() {}

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

loadFile();
}

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) {
return 0;
}

template <>
double EcalDeadChannelRecoveryBDTG<EBDetId>::recover(
const EBDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) {
bool isCrack = false;
int cellIndex = 0.;
double neighTotEn = 0.;
float val = 0.;

//find the matrix around id
std::vector<DetId> m3x3aroundDC = EcalClusterTools::matrixDetId(topology_, id, 1);
if (m3x3aroundDC.size() < 9) {
*acceptFlag = false;
return 0;
}

// Loop over all cells in the vector "NxNaroundDC", and for each cell find it's energy
// (from the EcalRecHits collection).
for (auto const &theCells : m3x3aroundDC) {
EBDetId cell = EBDetId(theCells);
if (cell == id) {
int iEtaCentral = std::abs(cell.ieta());
int iPhiCentral = cell.iphi();

if (iEtaCentral < 2 || std::abs(iEtaCentral - 25) < 2 || std::abs(iEtaCentral - 45) < 2 ||
std::abs(iEtaCentral - 65) < 2 || iEtaCentral > 83 || (int(iPhiCentral + 0.5) % 20 == 0))
isCrack = true;
}
if (!cell.null()) {
EcalRecHitCollection::const_iterator goS_it = hit_collection.find(cell);
if (goS_it != hit_collection.end() && cell != id) {
if (goS_it->energy() < single8Cut) {
*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();
cellIndex++;
}
} else if (cell == id) { // the cell is the central one
mx_.rEn[cellIndex] = 0;
cellIndex++;
} else { //goS_it is not in the rechitcollection
*acceptFlag = false;
return 0.;
}
} else { //cell is null
*acceptFlag = false;
return 0.;
}
}
if (cellIndex > 0 && neighTotEn >= single8Cut * 8. && neighTotEn >= sum8Cut) {
bool allneighs = true;
mx_.sumE8 = neighTotEn;
for (unsigned int icell = 0; icell < 9; icell++) {
if (mx_.rEn[icell] < single8Cut && icell != 4) {
allneighs = false;
}
mx_.rEn[icell] = mx_.rEn[icell] / neighTotEn;
}
if (allneighs == true) {
// evaluate the regression
if (isCrack) {
val = exp((readerCrack->EvaluateRegression("BDTG"))[0]);
} else {
val = exp((readerNoCrack->EvaluateRegression("BDTG"))[0]);
}

*acceptFlag = true;
//return the estimated energy
return val;
} else {
*acceptFlag = false;
return 0;
}
} else {
*acceptFlag = false;
return 0;
}
}

template class EcalDeadChannelRecoveryBDTG<EBDetId>;
template class EcalDeadChannelRecoveryBDTG<EEDetId>; //not used.
Expand Up @@ -80,9 +80,10 @@ void EcalDeadChannelRecoveryProducers<DetIdT>::produce(edm::Event& evt, const ed
if (it->detid() == *CheckDead) {
OverADeadRecHit = true;
bool AcceptRecHit = true;
EcalRecHit hit = deadChannelCorrector.correct(
it->detid(), *hit_collection, CorrectionMethod_, Sum8GeVThreshold_, &AcceptRecHit);

float dummy = 0;
float ebEn = deadChannelCorrector.correct(
it->detid(), *hit_collection, CorrectionMethod_, Sum8GeVThreshold_, dummy, &AcceptRecHit);
EcalRecHit hit(it->detid(), ebEn, 0., EcalRecHit::kDead);
if (hit.energy() != 0 and AcceptRecHit == true) {
hit.setFlag(EcalRecHit::kNeighboursRecovered);
} else {
Expand Down
7 changes: 7 additions & 0 deletions RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.cc
Expand Up @@ -308,6 +308,13 @@ void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descri
desc.add<edm::InputTag>("eeFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeFE"));
desc.add<edm::InputTag>("ebDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebDetId"));
desc.add<double>("singleChannelRecoveryThreshold", 8);
desc.add<double>("sum8ChannelRecoveryThreshold", 0.);
desc.add<edm::FileInPath>("bdtWeightFileNoCracks",
edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
"bdtgAllRH_8GT700MeV_noCracks_ZskimData2017_v1.xml"));
desc.add<edm::FileInPath>("bdtWeightFileCracks",
edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
"bdtgAllRH_8GT700MeV_onlyCracks_ZskimData2017_v1.xml"));
{
std::vector<std::string> temp1;
temp1.reserve(9);
Expand Down
14 changes: 10 additions & 4 deletions RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc
Expand Up @@ -24,6 +24,7 @@ EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet& ps, ed
// isolated channel recovery
singleRecoveryMethod_ = ps.getParameter<std::string>("singleChannelRecoveryMethod");
singleRecoveryThreshold_ = ps.getParameter<double>("singleChannelRecoveryThreshold");
sum8RecoveryThreshold_ = ps.getParameter<double>("sum8ChannelRecoveryThreshold");
killDeadChannels_ = ps.getParameter<bool>("killDeadChannels");
recoverEBIsolatedChannels_ = ps.getParameter<bool>("recoverEBIsolatedChannels");
recoverEEIsolatedChannels_ = ps.getParameter<bool>("recoverEEIsolatedChannels");
Expand All @@ -40,6 +41,9 @@ EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet& ps, ed

tpDigiToken_ =
c.consumes<EcalTrigPrimDigiCollection>(ps.getParameter<edm::InputTag>("triggerPrimitiveDigiCollection"));

if (recoverEBIsolatedChannels_ && singleRecoveryMethod_ == "BDTG")
ebDeadChannelCorrector.setParameters(ps);
}

void EcalRecHitWorkerRecover::set(const edm::EventSetup& es) {
Expand Down Expand Up @@ -124,8 +128,9 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,

// channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
bool AcceptRecHit = true;
EcalRecHit hit =
ebDeadChannelCorrector.correct(detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, &AcceptRecHit);
float ebEn = ebDeadChannelCorrector.correct(
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
EcalRecHit hit(detId, ebEn, 0., EcalRecHit::kDead);

if (hit.energy() != 0 and AcceptRecHit == true) {
hit.setFlag(EcalRecHit::kNeighboursRecovered);
Expand All @@ -141,8 +146,9 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,

// channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE
bool AcceptRecHit = true;
EcalRecHit hit =
eeDeadChannelCorrector.correct(detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, &AcceptRecHit);
float eeEn = eeDeadChannelCorrector.correct(
detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
EcalRecHit hit(detId, eeEn, 0., EcalRecHit::kDead);
if (hit.energy() != 0 and AcceptRecHit == true) {
hit.setFlag(EcalRecHit::kNeighboursRecovered);
} else {
Expand Down
Expand Up @@ -50,6 +50,7 @@ class EcalRecHitWorkerRecover : public EcalRecHitWorkerBaseClass {
edm::ESHandle<EcalChannelStatus> chStatus_;

double singleRecoveryThreshold_;
double sum8RecoveryThreshold_;
std::string singleRecoveryMethod_;
bool killDeadChannels_;

Expand Down

0 comments on commit ffabcf0

Please sign in to comment.