Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Ecal DPG - Deploy new software (based on BDT) to recover isolated dead channels… #27175

Merged
merged 6 commits into from Aug 2, 2019
Merged

Ecal DPG - Deploy new software (based on BDT) to recover isolated dead channels… #27175

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
b7d7959
Final (?) fix to the chi2 problem. In the end the 'source' of the dif…
nancymarinelli Jul 29, 2019
58ca071
Code formatting
nancymarinelli Jul 29, 2019
da9d251
And now we switch OFF the channel recovery until DPG request ot switc…
nancymarinelli Jul 30, 2019
ec9e102
From Silvia: Unique pointers to the tmva readers, indentation, semico…
nancymarinelli Jul 31, 2019
155b8aa
code formatting
nancymarinelli Jul 31, 2019
edab6bc
Remove algo as argument from setTopology
nancymarinelli Jul 31, 2019
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
1 change: 1 addition & 0 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml
View file
Open in desktop
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
21 changes: 15 additions & 6 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h
View file
Open in desktop
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 setCaloTopology(const CaloTopology *topology);
EcalRecHit correct(
const DetIdT id, const EcalRecHitCollection &hit_collection, std::string algo, double Sum8Cut, bool *AccFlag);
void setParameters(const edm::ParameterSet &ps);
void setCaloTopology(std::string algo, const CaloTopology *topology);
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_;
nancymarinelli marked this conversation as resolved.
Show resolved Hide resolved
};
#endif // RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryAlgos_HH
#endif
52 changes: 52 additions & 0 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h
View file
Open in desktop
@@ -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
33 changes: 22 additions & 11 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc
View file
Open in desktop
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>::setCaloTopology(const CaloTopology *topo) {
nn.setCaloTopology(topo);
void EcalDeadChannelRecoveryAlgos<T>::setParameters(const edm::ParameterSet &ps) {
bdtg_.setParameters(ps);
}

template <typename T>
EcalRecHit EcalDeadChannelRecoveryAlgos<T>::correct(
const T id, const EcalRecHitCollection &hit_collection, std::string algo, double Sum8Cut, bool *AcceptFlag) {
// recover as single dead channel
double NewEnergy = 0.0;
void EcalDeadChannelRecoveryAlgos<T>::setCaloTopology(std::string algo, const CaloTopology *topo) {
bdtg_.setCaloTopology(topo);
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm missing the reason why the first argument is needed.
It looks vestigial.
Please remove or clarify why needed.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In a way it's vestigial because it was used when the two recovery algorythms were coexisting. Although we have decided to throw away the NN, we do not want to remove the possibility for a future to include an other algo

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

but what would be the reason to keep it here, after the choice of the algo is made?
Even in the older variant of the code where both NN and BDTG were (unnecessarily) instantiated together there was still an assumption that only one of them can run.
If you really need it, the choice of the algo should be done when the EcalDeadChannelRecoveryAlgos is constructed, not in this method.
So, I do not see a reason for the algo parameter here.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The choice is made in https://github.com/nancymarinelli/cmssw/blob/deadChannelRecoveryBDT_V1/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc#L127 and the setCaloToplogy is necessary to pass the topology to the chosen algo. At this point in time, and with all the ifs related to the main validation, we do not certainly feel to make any change to how the framework was setup

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

right, I know this point of entry for the topology.
My comment still stands.
the choice of the algorithm is already made earlier (in ebDeadChannelCorrector.setParameters(ps)) there is no reason or a reasonable way to modify it here event by event.

}

if (algo == "NeuralNetworks") {
NewEnergy = this->nn.recover(id, hit_collection, Sum8Cut, AcceptFlag);
template <typename T>
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 == "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
149 changes: 149 additions & 0 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc
View file
Open in desktop
@@ -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() {}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

it looks like the two methods above are not needed (defaults will be created by the compiler anyways)

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Do they make any harm ? I personally like to have them explicit there

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The compiler is actually a good buddy of us! For example, when the destructor is compiler-generated, it is implicitly declared as noexcept which is one of the recommendations in the C++ Core Guidelines.

Edit: should have read the guidelines more carefully before posting... Just realized the same implicit noexcept goes for user-defined destructors too. Okay then this was not a good argument :)


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.
10 changes: 6 additions & 4 deletions RecoLocalCalo/EcalDeadChannelRecoveryAlgos/test/plugins/EcalDeadChannelRecoveryProducers.cc
View file
Open in desktop
Expand Up @@ -67,7 +67,8 @@ void EcalDeadChannelRecoveryProducers<DetIdT>::produce(edm::Event& evt, const ed
// create a unique_ptr to a EcalRecHitCollection, copy the RecHits into it and
// put it in the Event:
auto redCollection = std::make_unique<EcalRecHitCollection>();
deadChannelCorrector.setCaloTopology(theCaloTopology.product());
std::string dummy;
deadChannelCorrector.setCaloTopology(dummy, theCaloTopology.product());

//
// Double loop over EcalRecHit collection and "dead" cell RecHits.
Expand All @@ -80,9 +81,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
View file
Open in desktop
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
18 changes: 12 additions & 6 deletions RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc
View file
Open in desktop
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 @@ -120,12 +124,13 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,

if (flags == EcalRecHitWorkerRecover::EB_single) {
// recover as single dead channel
ebDeadChannelCorrector.setCaloTopology(caloTopology_.product());
ebDeadChannelCorrector.setCaloTopology(singleRecoveryMethod_, caloTopology_.product());

// 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 @@ -137,12 +142,13 @@ bool EcalRecHitWorkerRecover::run(const edm::Event& evt,

} else if (flags == EcalRecHitWorkerRecover::EE_single) {
// recover as single dead channel
eeDeadChannelCorrector.setCaloTopology(caloTopology_.product());
eeDeadChannelCorrector.setCaloTopology(singleRecoveryMethod_, caloTopology_.product());

// 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
1 change: 1 addition & 0 deletions RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h
View file
Open in desktop
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