From 09ae77160f71eac54cb171be0497f32772afdf29 Mon Sep 17 00:00:00 2001 From: nancy Date: Mon, 29 Jul 2019 12:21:41 +0200 Subject: [PATCH] Cherry pick files from PR27175 for Backport --- .../BuildFile.xml | 1 + .../interface/EcalDeadChannelRecoveryAlgos.h | 27 +- .../interface/EcalDeadChannelRecoveryBDTG.h | 52 ++ .../src/EcalDeadChannelRecoveryAlgos.cc | 38 +- .../src/EcalDeadChannelRecoveryBDTG.cc | 150 ++++ .../EcalDeadChannelRecoveryProducers.cc | 34 +- .../plugins/EcalRecHitProducer.cc | 547 ++++++------ .../plugins/EcalRecHitWorkerRecover.cc | 783 +++++++++--------- .../plugins/EcalRecHitWorkerRecover.h | 127 ++- .../EcalRecProducers/python/ecalRecHit_cfi.py | 13 +- 10 files changed, 986 insertions(+), 786 deletions(-) create mode 100644 RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h create mode 100644 RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml index 72df18a54ebf1..eab37ca28a8e6 100644 --- a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml @@ -3,6 +3,7 @@ + diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h index 15986fbf8c5a8..a94f55508e3ca 100644 --- a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h @@ -7,16 +7,25 @@ #include "DataFormats/EcalDetId/interface/EBDetId.h" #include -#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryNN.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" -template class EcalDeadChannelRecoveryAlgos { - public: +#include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h" + +#include "FWCore/ParameterSet/interface/ParameterSet.h" + +template +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 nn; +private: + EcalDeadChannelRecoveryBDTG bdtg_; }; -#endif // RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryAlgos_HH +#endif diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h new file mode 100644 index 0000000000000..58faf836b6594 --- /dev/null +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h @@ -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 +#include + +template +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 rEn, ieta, iphi; + float sumE8; + }; + + XtalMatrix mx_; + + edm::FileInPath bdtWeightFileNoCracks_; + edm::FileInPath bdtWeightFileCracks_; + + TMVA::Reader *readerNoCrack; + TMVA::Reader *readerCrack; +}; + +#endif diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc index d35a2a249cb9b..6c09ffe12f528 100644 --- a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc @@ -8,33 +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 -void EcalDeadChannelRecoveryAlgos::setCaloTopology( - const CaloTopology *topo) { - nn.setCaloTopology(topo); +void EcalDeadChannelRecoveryAlgos::setParameters(const edm::ParameterSet &ps) { + bdtg_.setParameters(ps); } template -EcalRecHit EcalDeadChannelRecoveryAlgos::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::setCaloTopology(const CaloTopology *topo) { + bdtg_.setCaloTopology(topo); +} - if (algo == "NeuralNetworks") { - NewEnergy = this->nn.recover(id, hit_collection, Sum8Cut, AcceptFlag); +template +float EcalDeadChannelRecoveryAlgos::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; + edm::LogError("EcalDeadChannelRecoveryAlgos") << "Invalid algorithm for dead channel recovery."; + *acceptFlag = false; } - uint32_t flag = 0; - return EcalRecHit(id, NewEnergy, 0, flag); + return newEnergy; } template class EcalDeadChannelRecoveryAlgos; diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc new file mode 100644 index 0000000000000..4a194135fcee3 --- /dev/null +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc @@ -0,0 +1,150 @@ +// 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 + +#include +#include +#include +#include + +template <> +void EcalDeadChannelRecoveryBDTG::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::loadFile() { + readerNoCrack = new TMVA::Reader("!Color:!Silent"); + readerCrack = new TMVA::Reader("!Color:!Silent"); + + this->addVariables(readerNoCrack); + this->addVariables(readerCrack); + + reco::details::loadTMVAWeights(readerNoCrack, "BDTG", bdtWeightFileNoCracks_.fullPath()); + reco::details::loadTMVAWeights(readerCrack, "BDTG", bdtWeightFileCracks_.fullPath()); +} + +template +EcalDeadChannelRecoveryBDTG::EcalDeadChannelRecoveryBDTG() {} + +template +EcalDeadChannelRecoveryBDTG::~EcalDeadChannelRecoveryBDTG() {} + +template <> +void EcalDeadChannelRecoveryBDTG::setParameters(const edm::ParameterSet &ps) { + bdtWeightFileNoCracks_ = ps.getParameter("bdtWeightFileNoCracks"); + bdtWeightFileCracks_ = ps.getParameter("bdtWeightFileCracks"); + + this->loadFile(); +} + +template <> +void EcalDeadChannelRecoveryBDTG::setParameters(const edm::ParameterSet &ps) {} + +template <> +double EcalDeadChannelRecoveryBDTG::recover( + const EEDetId id, const EcalRecHitCollection &hit_collection, double single8Cut, double sum8Cut, bool *acceptFlag) { + return 0; +} + +template <> +double EcalDeadChannelRecoveryBDTG::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 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; +template class EcalDeadChannelRecoveryBDTG; //not used. diff --git a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/test/plugins/EcalDeadChannelRecoveryProducers.cc b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/test/plugins/EcalDeadChannelRecoveryProducers.cc index 818bc2ff80fe6..3233f3441b9f5 100644 --- a/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/test/plugins/EcalDeadChannelRecoveryProducers.cc +++ b/RecoLocalCalo/EcalDeadChannelRecoveryAlgos/test/plugins/EcalDeadChannelRecoveryProducers.cc @@ -28,8 +28,7 @@ using namespace cms; using namespace std; template -EcalDeadChannelRecoveryProducers::EcalDeadChannelRecoveryProducers( - const edm::ParameterSet& ps) { +EcalDeadChannelRecoveryProducers::EcalDeadChannelRecoveryProducers(const edm::ParameterSet& ps) { //now do what ever other initialization is needed CorrectDeadCells_ = ps.getParameter("CorrectDeadCells"); CorrectionMethod_ = ps.getParameter("CorrectionMethod"); @@ -37,8 +36,7 @@ EcalDeadChannelRecoveryProducers::EcalDeadChannelRecoveryProducers( DeadChannelFileName_ = ps.getParameter("DeadChannelsFile"); Sum8GeVThreshold_ = ps.getParameter("Sum8GeVThreshold"); - hitToken_ = - consumes(ps.getParameter("hitTag")); + hitToken_ = consumes(ps.getParameter("hitTag")); produces(reducedHitCollection_); } @@ -50,8 +48,7 @@ EcalDeadChannelRecoveryProducers::~EcalDeadChannelRecoveryProducers() { // ------------ method called to produce the data ------------ template -void EcalDeadChannelRecoveryProducers::produce( - edm::Event& evt, const edm::EventSetup& iSetup) { +void EcalDeadChannelRecoveryProducers::produce(edm::Event& evt, const edm::EventSetup& iSetup) { using namespace edm; edm::ESHandle theCaloTopology; @@ -76,18 +73,17 @@ void EcalDeadChannelRecoveryProducers::produce( // Double loop over EcalRecHit collection and "dead" cell RecHits. // If we step into a "dead" cell call "deadChannelCorrector::correct()" // - for (EcalRecHitCollection::const_iterator it = hit_collection->begin(); - it != hit_collection->end(); ++it) { + for (EcalRecHitCollection::const_iterator it = hit_collection->begin(); it != hit_collection->end(); ++it) { std::vector::const_iterator CheckDead = ChannelsDeadID.begin(); bool OverADeadRecHit = false; while (CheckDead != ChannelsDeadID.end()) { 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 { @@ -109,7 +105,8 @@ void EcalDeadChannelRecoveryProducers::produce( } // method called once each job just before starting event loop ------------ -template <> void EcalDeadChannelRecoveryProducers::beginJob() { +template <> +void EcalDeadChannelRecoveryProducers::beginJob() { //Open the DeadChannel file, read it. FILE* DeadCha; printf("Dead Channels FILE: %s\n", DeadChannelFileName_.c_str()); @@ -120,7 +117,6 @@ template <> void EcalDeadChannelRecoveryProducers::beginJob() { int iphi = -10000; while (fileStatus != EOF) { - fileStatus = fscanf(DeadCha, "%d %d\n", &ieta, &iphi); // Problem reading Dead Channels file @@ -138,7 +134,8 @@ template <> void EcalDeadChannelRecoveryProducers::beginJob() { fclose(DeadCha); } -template <> void EcalDeadChannelRecoveryProducers::beginJob() { +template <> +void EcalDeadChannelRecoveryProducers::beginJob() { //Open the DeadChannel file, read it. FILE* DeadCha; printf("Dead Channels FILE: %s\n", DeadChannelFileName_.c_str()); @@ -149,7 +146,6 @@ template <> void EcalDeadChannelRecoveryProducers::beginJob() { int iy = -10000; int iz = -10000; while (fileStatus != EOF) { - fileStatus = fscanf(DeadCha, "%d %d %d\n", &ix, &iy, &iz); // Problem reading Dead Channels file @@ -172,10 +168,8 @@ template <> void EcalDeadChannelRecoveryProducers::beginJob() { template void EcalDeadChannelRecoveryProducers::endJob() {} -typedef class EcalDeadChannelRecoveryProducers - EBDeadChannelRecoveryProducers; -typedef class EcalDeadChannelRecoveryProducers - EEDeadChannelRecoveryProducers; +typedef class EcalDeadChannelRecoveryProducers EBDeadChannelRecoveryProducers; +typedef class EcalDeadChannelRecoveryProducers EEDeadChannelRecoveryProducers; DEFINE_FWK_MODULE(EBDeadChannelRecoveryProducers); DEFINE_FWK_MODULE(EEDeadChannelRecoveryProducers); diff --git a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.cc b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.cc index 0d7646d84d1a8..0bf9119382d63 100644 --- a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.cc +++ b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.cc @@ -27,298 +27,294 @@ #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" #include "FWCore/ParameterSet/interface/ParameterSetDescription.h" -EcalRecHitProducer::EcalRecHitProducer(const edm::ParameterSet& ps) -{ - - ebRechitCollection_ = ps.getParameter("EBrechitCollection"); - eeRechitCollection_ = ps.getParameter("EErechitCollection"); - - recoverEBIsolatedChannels_ = ps.getParameter("recoverEBIsolatedChannels"); - recoverEEIsolatedChannels_ = ps.getParameter("recoverEEIsolatedChannels"); - recoverEBVFE_ = ps.getParameter("recoverEBVFE"); - recoverEEVFE_ = ps.getParameter("recoverEEVFE"); - recoverEBFE_ = ps.getParameter("recoverEBFE"); - recoverEEFE_ = ps.getParameter("recoverEEFE"); - killDeadChannels_ = ps.getParameter("killDeadChannels"); - - - produces< EBRecHitCollection >(ebRechitCollection_); - produces< EERecHitCollection >(eeRechitCollection_); - - - ebUncalibRecHitToken_ = - consumes( ps.getParameter("EBuncalibRecHitCollection")); - - eeUncalibRecHitToken_ = - consumes( ps.getParameter("EEuncalibRecHitCollection")); - - ebDetIdToBeRecoveredToken_ = - consumes>(ps.getParameter("ebDetIdToBeRecovered")); - - eeDetIdToBeRecoveredToken_= - consumes>(ps.getParameter("eeDetIdToBeRecovered")); - - ebFEToBeRecoveredToken_ = consumes>(ps.getParameter("ebFEToBeRecovered")); - - eeFEToBeRecoveredToken_= consumes>( ps.getParameter("eeFEToBeRecovered")) ; - - std::string componentType = ps.getParameter("algo"); - edm::ConsumesCollector c{consumesCollector()}; - worker_ = std::unique_ptr{EcalRecHitWorkerFactory::get()->create(componentType, ps, c)}; - - // to recover problematic channels - componentType = ps.getParameter("algoRecover"); - workerRecover_ = std::unique_ptr{EcalRecHitWorkerFactory::get()->create(componentType, ps, c)}; - - edm::ParameterSet cleaningPs = - ps.getParameter("cleaningConfig"); - cleaningAlgo_ = std::make_unique(cleaningPs); +EcalRecHitProducer::EcalRecHitProducer(const edm::ParameterSet& ps) { + ebRechitCollection_ = ps.getParameter("EBrechitCollection"); + eeRechitCollection_ = ps.getParameter("EErechitCollection"); + + recoverEBIsolatedChannels_ = ps.getParameter("recoverEBIsolatedChannels"); + recoverEEIsolatedChannels_ = ps.getParameter("recoverEEIsolatedChannels"); + recoverEBVFE_ = ps.getParameter("recoverEBVFE"); + recoverEEVFE_ = ps.getParameter("recoverEEVFE"); + recoverEBFE_ = ps.getParameter("recoverEBFE"); + recoverEEFE_ = ps.getParameter("recoverEEFE"); + killDeadChannels_ = ps.getParameter("killDeadChannels"); + + produces(ebRechitCollection_); + produces(eeRechitCollection_); + + ebUncalibRecHitToken_ = + consumes(ps.getParameter("EBuncalibRecHitCollection")); + + eeUncalibRecHitToken_ = + consumes(ps.getParameter("EEuncalibRecHitCollection")); + + ebDetIdToBeRecoveredToken_ = consumes>(ps.getParameter("ebDetIdToBeRecovered")); + + eeDetIdToBeRecoveredToken_ = consumes>(ps.getParameter("eeDetIdToBeRecovered")); + + ebFEToBeRecoveredToken_ = consumes>(ps.getParameter("ebFEToBeRecovered")); + + eeFEToBeRecoveredToken_ = consumes>(ps.getParameter("eeFEToBeRecovered")); + + std::string componentType = ps.getParameter("algo"); + edm::ConsumesCollector c{consumesCollector()}; + worker_ = EcalRecHitWorkerFactory::get()->create(componentType, ps, c); + + // to recover problematic channels + componentType = ps.getParameter("algoRecover"); + workerRecover_ = EcalRecHitWorkerFactory::get()->create(componentType, ps, c); + + edm::ParameterSet cleaningPs = ps.getParameter("cleaningConfig"); + cleaningAlgo_ = std::make_unique(cleaningPs); } EcalRecHitProducer::~EcalRecHitProducer() = default; -void -EcalRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) -{ - using namespace edm; - - Handle< EBUncalibratedRecHitCollection > pEBUncalibRecHits; - Handle< EEUncalibratedRecHitCollection > pEEUncalibRecHits; - - const EBUncalibratedRecHitCollection* ebUncalibRecHits = nullptr; - const EEUncalibratedRecHitCollection* eeUncalibRecHits = nullptr; - - // get the barrel uncalib rechit collection - - evt.getByToken( ebUncalibRecHitToken_, pEBUncalibRecHits); - ebUncalibRecHits = pEBUncalibRecHits.product(); - LogDebug("EcalRecHitDebug") << "total # EB uncalibrated rechits: " << ebUncalibRecHits->size(); - - - - evt.getByToken( eeUncalibRecHitToken_, pEEUncalibRecHits); - eeUncalibRecHits = pEEUncalibRecHits.product(); // get a ptr to the product - LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits: " << eeUncalibRecHits->size(); - - // collection of rechits to put in the event - auto ebRecHits = std::make_unique(); - auto eeRecHits = std::make_unique(); - - worker_->set(es); - - if ( recoverEBIsolatedChannels_ || recoverEEIsolatedChannels_ - || recoverEBFE_ || recoverEEFE_ - || recoverEBVFE_ || recoverEEVFE_ - || killDeadChannels_ ) { - workerRecover_->set(es); - } +void EcalRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) { + using namespace edm; - if (ebUncalibRecHits) - { - // loop over uncalibrated rechits to make calibrated ones - for(EBUncalibratedRecHitCollection::const_iterator it = ebUncalibRecHits->begin(); it != ebUncalibRecHits->end(); ++it) { - worker_->run(evt, *it, *ebRecHits); - } - } + Handle pEBUncalibRecHits; + Handle pEEUncalibRecHits; - if (eeUncalibRecHits) - { - // loop over uncalibrated rechits to make calibrated ones - for(EEUncalibratedRecHitCollection::const_iterator it = eeUncalibRecHits->begin(); it != eeUncalibRecHits->end(); ++it) { - worker_->run(evt, *it, *eeRecHits); - } - } + const EBUncalibratedRecHitCollection* ebUncalibRecHits = nullptr; + const EEUncalibratedRecHitCollection* eeUncalibRecHits = nullptr; - // sort collections before attempting recovery, to avoid insertion of double recHits - ebRecHits->sort(); - eeRecHits->sort(); - - if ( recoverEBIsolatedChannels_ || recoverEBFE_ || killDeadChannels_ ) - { - edm::Handle< std::set > pEBDetId; - const std::set * detIds = nullptr; - evt.getByToken( ebDetIdToBeRecoveredToken_, pEBDetId); - detIds = pEBDetId.product(); - - - if ( detIds ) { - edm::ESHandle chStatus; - es.get().get(chStatus); - for( std::set::const_iterator it = detIds->begin(); it != detIds->end(); ++it ) { - // get channel status map to treat dead VFE separately - EcalChannelStatusMap::const_iterator chit = chStatus->find( *it ); - EcalChannelStatusCode chStatusCode; - if ( chit != chStatus->end() ) { - chStatusCode = *chit; - } else { - edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " - << (*it).rawId() - << "! something wrong with EcalChannelStatus in your DB? "; - } - EcalUncalibratedRecHit urh; - if ( chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE ) { // dead VFE (from DB info) - // uses the EcalUncalibratedRecHit to pass the DetId info - urh = EcalUncalibratedRecHit( *it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_VFE ); - if ( recoverEBVFE_ || killDeadChannels_ ) workerRecover_->run( evt, urh, *ebRecHits ); - } else { - // uses the EcalUncalibratedRecHit to pass the DetId info - urh = EcalUncalibratedRecHit( *it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_single ); - if ( recoverEBIsolatedChannels_ || killDeadChannels_ ) workerRecover_->run( evt, urh, *ebRecHits ); - } - - } - } - } + // get the barrel uncalib rechit collection - if ( recoverEEIsolatedChannels_ || recoverEEVFE_ || killDeadChannels_ ) - { - edm::Handle< std::set > pEEDetId; - const std::set * detIds = nullptr; - - evt.getByToken( eeDetIdToBeRecoveredToken_, pEEDetId); - detIds = pEEDetId.product(); - - if ( detIds ) { - edm::ESHandle chStatus; - es.get().get(chStatus); - for( std::set::const_iterator it = detIds->begin(); it != detIds->end(); ++it ) { - // get channel status map to treat dead VFE separately - EcalChannelStatusMap::const_iterator chit = chStatus->find( *it ); - EcalChannelStatusCode chStatusCode; - if ( chit != chStatus->end() ) { - chStatusCode = *chit; - } else { - edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " - << (*it).rawId() - << "! something wrong with EcalChannelStatus in your DB? "; - } - EcalUncalibratedRecHit urh; - if ( chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) { // dead VFE (from DB info) - // uses the EcalUncalibratedRecHit to pass the DetId info - urh = EcalUncalibratedRecHit( *it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_VFE ); - if ( recoverEEVFE_ || killDeadChannels_ ) workerRecover_->run( evt, urh, *eeRecHits ); - } else { - // uses the EcalUncalibratedRecHit to pass the DetId info - urh = EcalUncalibratedRecHit( *it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_single ); - if ( recoverEEIsolatedChannels_ || killDeadChannels_ ) workerRecover_->run( evt, urh, *eeRecHits ); - } - } - } + evt.getByToken(ebUncalibRecHitToken_, pEBUncalibRecHits); + ebUncalibRecHits = pEBUncalibRecHits.product(); + LogDebug("EcalRecHitDebug") << "total # EB uncalibrated rechits: " << ebUncalibRecHits->size(); + + evt.getByToken(eeUncalibRecHitToken_, pEEUncalibRecHits); + eeUncalibRecHits = pEEUncalibRecHits.product(); // get a ptr to the product + LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits: " << eeUncalibRecHits->size(); + + // collection of rechits to put in the event + auto ebRecHits = std::make_unique(); + auto eeRecHits = std::make_unique(); + + worker_->set(es); + + if (recoverEBIsolatedChannels_ || recoverEEIsolatedChannels_ || recoverEBFE_ || recoverEEFE_ || recoverEBVFE_ || + recoverEEVFE_ || killDeadChannels_) { + workerRecover_->set(es); + } + + if (ebUncalibRecHits) { + // loop over uncalibrated rechits to make calibrated ones + for (EBUncalibratedRecHitCollection::const_iterator it = ebUncalibRecHits->begin(); it != ebUncalibRecHits->end(); + ++it) { + worker_->run(evt, *it, *ebRecHits); + } + } + + if (eeUncalibRecHits) { + // loop over uncalibrated rechits to make calibrated ones + for (EEUncalibratedRecHitCollection::const_iterator it = eeUncalibRecHits->begin(); it != eeUncalibRecHits->end(); + ++it) { + worker_->run(evt, *it, *eeRecHits); + } + } + + // sort collections before attempting recovery, to avoid insertion of double recHits + ebRecHits->sort(); + eeRecHits->sort(); + + if (recoverEBIsolatedChannels_ || recoverEBFE_ || killDeadChannels_) { + edm::Handle> pEBDetId; + const std::set* detIds = nullptr; + evt.getByToken(ebDetIdToBeRecoveredToken_, pEBDetId); + detIds = pEBDetId.product(); + + if (detIds) { + edm::ESHandle chStatus; + es.get().get(chStatus); + for (std::set::const_iterator it = detIds->begin(); it != detIds->end(); ++it) { + // get channel status map to treat dead VFE separately + EcalChannelStatusMap::const_iterator chit = chStatus->find(*it); + EcalChannelStatusCode chStatusCode; + if (chit != chStatus->end()) { + chStatusCode = *chit; + } else { + edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << (*it).rawId() + << "! something wrong with EcalChannelStatus in your DB? "; + } + EcalUncalibratedRecHit urh; + if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) { // dead VFE (from DB info) + // uses the EcalUncalibratedRecHit to pass the DetId info + urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_VFE); + if (recoverEBVFE_ || killDeadChannels_) + workerRecover_->run(evt, urh, *ebRecHits); + } else { + // uses the EcalUncalibratedRecHit to pass the DetId info + urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_single); + if (recoverEBIsolatedChannels_ || killDeadChannels_) + workerRecover_->run(evt, urh, *ebRecHits); } + } + } + } - if ( recoverEBFE_ || killDeadChannels_ ) - { - edm::Handle< std::set > pEBFEId; - const std::set * ttIds = nullptr; - - evt.getByToken( ebFEToBeRecoveredToken_, pEBFEId); - ttIds = pEBFEId.product(); - - if ( ttIds ) { - for( std::set::const_iterator it = ttIds->begin(); it != ttIds->end(); ++it ) { - // uses the EcalUncalibratedRecHit to pass the DetId info - int ieta = (((*it).ietaAbs()-1)*5+1)*(*it).zside(); // from EcalTrigTowerConstituentsMap - int iphi = (((*it).iphi()-1)*5+11)%360; // from EcalTrigTowerConstituentsMap - if( iphi <= 0 ) iphi += 360; // from EcalTrigTowerConstituentsMap - EcalUncalibratedRecHit urh( EBDetId(ieta, iphi, EBDetId::ETAPHIMODE), 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_FE ); - workerRecover_->run( evt, urh, *ebRecHits ); - } - } + if (recoverEEIsolatedChannels_ || recoverEEVFE_ || killDeadChannels_) { + edm::Handle> pEEDetId; + const std::set* detIds = nullptr; + + evt.getByToken(eeDetIdToBeRecoveredToken_, pEEDetId); + detIds = pEEDetId.product(); + + if (detIds) { + edm::ESHandle chStatus; + es.get().get(chStatus); + for (std::set::const_iterator it = detIds->begin(); it != detIds->end(); ++it) { + // get channel status map to treat dead VFE separately + EcalChannelStatusMap::const_iterator chit = chStatus->find(*it); + EcalChannelStatusCode chStatusCode; + if (chit != chStatus->end()) { + chStatusCode = *chit; + } else { + edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << (*it).rawId() + << "! something wrong with EcalChannelStatus in your DB? "; } + EcalUncalibratedRecHit urh; + if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) { // dead VFE (from DB info) + // uses the EcalUncalibratedRecHit to pass the DetId info + urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_VFE); + if (recoverEEVFE_ || killDeadChannels_) + workerRecover_->run(evt, urh, *eeRecHits); + } else { + // uses the EcalUncalibratedRecHit to pass the DetId info + urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_single); + if (recoverEEIsolatedChannels_ || killDeadChannels_) + workerRecover_->run(evt, urh, *eeRecHits); + } + } + } + } - if ( recoverEEFE_ || killDeadChannels_ ) - { - edm::Handle< std::set > pEEFEId; - const std::set * scIds = nullptr; - - evt.getByToken( eeFEToBeRecoveredToken_, pEEFEId); - scIds = pEEFEId.product(); - - - if ( scIds ) { - for( std::set::const_iterator it = scIds->begin(); it != scIds->end(); ++it ) { - // uses the EcalUncalibratedRecHit to pass the DetId info - if (EEDetId::validDetId( ((*it).ix()-1)*5+1, ((*it).iy()-1)*5+1, (*it).zside() )) { - EcalUncalibratedRecHit urh( EEDetId( ((*it).ix()-1)*5+1, ((*it).iy()-1)*5+1, (*it).zside() ), 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_FE ); - workerRecover_->run( evt, urh, *eeRecHits ); - } - } - } + if (recoverEBFE_ || killDeadChannels_) { + edm::Handle> pEBFEId; + const std::set* ttIds = nullptr; + + evt.getByToken(ebFEToBeRecoveredToken_, pEBFEId); + ttIds = pEBFEId.product(); + + if (ttIds) { + for (std::set::const_iterator it = ttIds->begin(); it != ttIds->end(); ++it) { + // uses the EcalUncalibratedRecHit to pass the DetId info + int ieta = (((*it).ietaAbs() - 1) * 5 + 1) * (*it).zside(); // from EcalTrigTowerConstituentsMap + int iphi = (((*it).iphi() - 1) * 5 + 11) % 360; // from EcalTrigTowerConstituentsMap + if (iphi <= 0) + iphi += 360; // from EcalTrigTowerConstituentsMap + EcalUncalibratedRecHit urh( + EBDetId(ieta, iphi, EBDetId::ETAPHIMODE), 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_FE); + workerRecover_->run(evt, urh, *ebRecHits); + } + } + } + + if (recoverEEFE_ || killDeadChannels_) { + edm::Handle> pEEFEId; + const std::set* scIds = nullptr; + + evt.getByToken(eeFEToBeRecoveredToken_, pEEFEId); + scIds = pEEFEId.product(); + + if (scIds) { + for (std::set::const_iterator it = scIds->begin(); it != scIds->end(); ++it) { + // uses the EcalUncalibratedRecHit to pass the DetId info + if (EEDetId::validDetId(((*it).ix() - 1) * 5 + 1, ((*it).iy() - 1) * 5 + 1, (*it).zside())) { + EcalUncalibratedRecHit urh(EEDetId(((*it).ix() - 1) * 5 + 1, ((*it).iy() - 1) * 5 + 1, (*it).zside()), + 0, + 0, + 0, + 0, + EcalRecHitWorkerBaseClass::EE_FE); + workerRecover_->run(evt, urh, *eeRecHits); } + } + } + } - // without re-sorting, find (used below in cleaning) will lead - // to undefined results - ebRecHits->sort(); - eeRecHits->sort(); - - // apply spike cleaning - if (cleaningAlgo_){ - cleaningAlgo_->setFlags(*ebRecHits); - cleaningAlgo_->setFlags(*eeRecHits); - } + // without re-sorting, find (used below in cleaning) will lead + // to undefined results + ebRecHits->sort(); + eeRecHits->sort(); + // apply spike cleaning + if (cleaningAlgo_) { + cleaningAlgo_->setFlags(*ebRecHits); + cleaningAlgo_->setFlags(*eeRecHits); + } - // put the collection of recunstructed hits in the event - LogInfo("EcalRecHitInfo") << "total # EB calibrated rechits: " << ebRecHits->size(); - LogInfo("EcalRecHitInfo") << "total # EE calibrated rechits: " << eeRecHits->size(); + // put the collection of recunstructed hits in the event + LogInfo("EcalRecHitInfo") << "total # EB calibrated rechits: " << ebRecHits->size(); + LogInfo("EcalRecHitInfo") << "total # EE calibrated rechits: " << eeRecHits->size(); - evt.put(std::move(ebRecHits), ebRechitCollection_); - evt.put(std::move(eeRecHits), eeRechitCollection_); + evt.put(std::move(ebRecHits), ebRechitCollection_); + evt.put(std::move(eeRecHits), eeRechitCollection_); } void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { edm::ParameterSetDescription desc; - desc.add("recoverEEVFE",false); - desc.add("EErechitCollection","EcalRecHitsEE"); - desc.add("recoverEBIsolatedChannels",false); - desc.add("recoverEBVFE",false); - desc.add("laserCorrection",true); - desc.add("EBLaserMIN",0.5); - desc.add("killDeadChannels",true); + desc.add("recoverEEVFE", false); + desc.add("EErechitCollection", "EcalRecHitsEE"); + desc.add("recoverEBIsolatedChannels", false); + desc.add("recoverEBVFE", false); + desc.add("laserCorrection", true); + desc.add("EBLaserMIN", 0.5); + desc.add("killDeadChannels", true); { std::vector temp1; temp1.reserve(3); temp1.push_back(14); temp1.push_back(78); temp1.push_back(142); - desc.add >("dbStatusToBeExcludedEB",temp1); + desc.add>("dbStatusToBeExcludedEB", temp1); } - desc.add("EEuncalibRecHitCollection",edm::InputTag("ecalMultiFitUncalibRecHit","EcalUncalibRecHitsEE")); + desc.add("EEuncalibRecHitCollection", + edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEE")); { std::vector temp1; temp1.reserve(3); temp1.push_back(14); temp1.push_back(78); temp1.push_back(142); - desc.add >("dbStatusToBeExcludedEE",temp1); + desc.add>("dbStatusToBeExcludedEE", temp1); } - desc.add("EELaserMIN",0.5); - desc.add("ebFEToBeRecovered",edm::InputTag("ecalDetIdToBeRecovered","ebFE")); + desc.add("EELaserMIN", 0.5); + desc.add("ebFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebFE")); { edm::ParameterSetDescription psd0; - psd0.add("e6e2thresh",0.04); - psd0.add("tightenCrack_e6e2_double",3); - psd0.add("e4e1Threshold_endcap",0.3); - psd0.add("tightenCrack_e4e1_single",3); - psd0.add("tightenCrack_e1_double",2); - psd0.add("cThreshold_barrel",4); - psd0.add("e4e1Threshold_barrel",0.08); - psd0.add("tightenCrack_e1_single",2); - psd0.add("e4e1_b_barrel",-0.024); - psd0.add("e4e1_a_barrel",0.04); - psd0.add("ignoreOutOfTimeThresh",1000000000.0); - psd0.add("cThreshold_endcap",15); - psd0.add("e4e1_b_endcap",-0.0125); - psd0.add("e4e1_a_endcap",0.02); - psd0.add("cThreshold_double",10); - desc.add("cleaningConfig",psd0); + psd0.add("e6e2thresh", 0.04); + psd0.add("tightenCrack_e6e2_double", 3); + psd0.add("e4e1Threshold_endcap", 0.3); + psd0.add("tightenCrack_e4e1_single", 3); + psd0.add("tightenCrack_e1_double", 2); + psd0.add("cThreshold_barrel", 4); + psd0.add("e4e1Threshold_barrel", 0.08); + psd0.add("tightenCrack_e1_single", 2); + psd0.add("e4e1_b_barrel", -0.024); + psd0.add("e4e1_a_barrel", 0.04); + psd0.add("ignoreOutOfTimeThresh", 1000000000.0); + psd0.add("cThreshold_endcap", 15); + psd0.add("e4e1_b_endcap", -0.0125); + psd0.add("e4e1_a_endcap", 0.02); + psd0.add("cThreshold_double", 10); + desc.add("cleaningConfig", psd0); } - desc.add("logWarningEtThreshold_EE_FE",50); - desc.add("eeDetIdToBeRecovered",edm::InputTag("ecalDetIdToBeRecovered","eeDetId")); - desc.add("recoverEBFE",true); - desc.add("eeFEToBeRecovered",edm::InputTag("ecalDetIdToBeRecovered","eeFE")); - desc.add("ebDetIdToBeRecovered",edm::InputTag("ecalDetIdToBeRecovered","ebDetId")); - desc.add("singleChannelRecoveryThreshold",8); + desc.add("logWarningEtThreshold_EE_FE", 50); + desc.add("eeDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeDetId")); + desc.add("recoverEBFE", true); + desc.add("eeFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeFE")); + desc.add("ebDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebDetId")); + desc.add("singleChannelRecoveryThreshold", 8); + desc.add("sum8ChannelRecoveryThreshold", 0.); + desc.add("bdtWeightFileNoCracks", + edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/" + "bdtgAllRH_8GT700MeV_noCracks_ZskimData2017_v1.xml")); + desc.add("bdtWeightFileCracks", + edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/" + "bdtgAllRH_8GT700MeV_onlyCracks_ZskimData2017_v1.xml")); { std::vector temp1; temp1.reserve(9); @@ -331,13 +327,13 @@ void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descri temp1.push_back("kDeadVFE"); temp1.push_back("kDeadFE"); temp1.push_back("kNoDataNoTP"); - desc.add >("ChannelStatusToBeExcluded",temp1); + desc.add>("ChannelStatusToBeExcluded", temp1); } - desc.add("EBrechitCollection","EcalRecHitsEB"); - desc.add("triggerPrimitiveDigiCollection",edm::InputTag("ecalDigis","EcalTriggerPrimitives")); - desc.add("recoverEEFE",true); - desc.add("singleChannelRecoveryMethod","NeuralNetworks"); - desc.add("EBLaserMAX",3.0); + desc.add("EBrechitCollection", "EcalRecHitsEB"); + desc.add("triggerPrimitiveDigiCollection", edm::InputTag("ecalDigis", "EcalTriggerPrimitives")); + desc.add("recoverEEFE", true); + desc.add("singleChannelRecoveryMethod", "NeuralNetworks"); + desc.add("EBLaserMAX", 3.0); { edm::ParameterSetDescription psd0; { @@ -347,7 +343,7 @@ void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descri temp2.push_back("kDAC"); temp2.push_back("kNoLaser"); temp2.push_back("kNoisy"); - psd0.add >("kGood",temp2); + psd0.add>("kGood", temp2); } { std::vector temp2; @@ -355,13 +351,13 @@ void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descri temp2.push_back("kFixedG0"); temp2.push_back("kNonRespondingIsolated"); temp2.push_back("kDeadVFE"); - psd0.add >("kNeighboursRecovered",temp2); + psd0.add>("kNeighboursRecovered", temp2); } { std::vector temp2; temp2.reserve(1); temp2.push_back("kNoDataNoTP"); - psd0.add >("kDead",temp2); + psd0.add>("kDead", temp2); } { std::vector temp2; @@ -369,25 +365,26 @@ void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descri temp2.push_back("kNNoisy"); temp2.push_back("kFixedG6"); temp2.push_back("kFixedG1"); - psd0.add >("kNoisy",temp2); + psd0.add>("kNoisy", temp2); } { std::vector temp2; temp2.reserve(1); temp2.push_back("kDeadFE"); - psd0.add >("kTowerRecovered",temp2); + psd0.add>("kTowerRecovered", temp2); } - desc.add("flagsMapDBReco",psd0); + desc.add("flagsMapDBReco", psd0); } - desc.add("EBuncalibRecHitCollection",edm::InputTag("ecalMultiFitUncalibRecHit","EcalUncalibRecHitsEB")); - desc.add("algoRecover","EcalRecHitWorkerRecover"); - desc.add("algo","EcalRecHitWorkerSimple"); - desc.add("EELaserMAX",8.0); - desc.add("logWarningEtThreshold_EB_FE",50); - desc.add("recoverEEIsolatedChannels",false); - desc.add("skipTimeCalib",false); - descriptions.add("ecalRecHit",desc); + desc.add("EBuncalibRecHitCollection", + edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEB")); + desc.add("algoRecover", "EcalRecHitWorkerRecover"); + desc.add("algo", "EcalRecHitWorkerSimple"); + desc.add("EELaserMAX", 8.0); + desc.add("logWarningEtThreshold_EB_FE", 50); + desc.add("recoverEEIsolatedChannels", false); + desc.add("skipTimeCalib", false); + descriptions.add("ecalRecHit", desc); } #include "FWCore/Framework/interface/MakerMacros.h" -DEFINE_FWK_MODULE( EcalRecHitProducer ); +DEFINE_FWK_MODULE(EcalRecHitProducer); diff --git a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc index 26c8fa1a4d104..cf3dc55d616ad 100644 --- a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc +++ b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.cc @@ -18,263 +18,268 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Framework/interface/EDProducer.h" - -EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet&ps, edm::ConsumesCollector& c) : - EcalRecHitWorkerBaseClass(ps,c) -{ - rechitMaker_ = std::make_unique(); - // isolated channel recovery - singleRecoveryMethod_ = ps.getParameter("singleChannelRecoveryMethod"); - singleRecoveryThreshold_ = ps.getParameter("singleChannelRecoveryThreshold"); - killDeadChannels_ = ps.getParameter("killDeadChannels"); - recoverEBIsolatedChannels_ = ps.getParameter("recoverEBIsolatedChannels"); - recoverEEIsolatedChannels_ = ps.getParameter("recoverEEIsolatedChannels"); - recoverEBVFE_ = ps.getParameter("recoverEBVFE"); - recoverEEVFE_ = ps.getParameter("recoverEEVFE"); - recoverEBFE_ = ps.getParameter("recoverEBFE"); - recoverEEFE_ = ps.getParameter("recoverEEFE"); - - dbStatusToBeExcludedEE_ = ps.getParameter >("dbStatusToBeExcludedEE"); - dbStatusToBeExcludedEB_ = ps.getParameter >("dbStatusToBeExcludedEB"); - - logWarningEtThreshold_EB_FE_ = ps.getParameter("logWarningEtThreshold_EB_FE"); - logWarningEtThreshold_EE_FE_ = ps.getParameter("logWarningEtThreshold_EE_FE"); - - tpDigiToken_ = - c.consumes(ps.getParameter("triggerPrimitiveDigiCollection")); +EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet& ps, edm::ConsumesCollector& c) + : EcalRecHitWorkerBaseClass(ps, c) { + rechitMaker_ = std::make_unique(); + // isolated channel recovery + singleRecoveryMethod_ = ps.getParameter("singleChannelRecoveryMethod"); + singleRecoveryThreshold_ = ps.getParameter("singleChannelRecoveryThreshold"); + sum8RecoveryThreshold_ = ps.getParameter("sum8ChannelRecoveryThreshold"); + killDeadChannels_ = ps.getParameter("killDeadChannels"); + recoverEBIsolatedChannels_ = ps.getParameter("recoverEBIsolatedChannels"); + recoverEEIsolatedChannels_ = ps.getParameter("recoverEEIsolatedChannels"); + recoverEBVFE_ = ps.getParameter("recoverEBVFE"); + recoverEEVFE_ = ps.getParameter("recoverEEVFE"); + recoverEBFE_ = ps.getParameter("recoverEBFE"); + recoverEEFE_ = ps.getParameter("recoverEEFE"); + + dbStatusToBeExcludedEE_ = ps.getParameter >("dbStatusToBeExcludedEE"); + dbStatusToBeExcludedEB_ = ps.getParameter >("dbStatusToBeExcludedEB"); + + logWarningEtThreshold_EB_FE_ = ps.getParameter("logWarningEtThreshold_EB_FE"); + logWarningEtThreshold_EE_FE_ = ps.getParameter("logWarningEtThreshold_EE_FE"); + + tpDigiToken_ = + c.consumes(ps.getParameter("triggerPrimitiveDigiCollection")); + + if (recoverEBIsolatedChannels_ && singleRecoveryMethod_ == "BDTG") + ebDeadChannelCorrector.setParameters(ps); } - -void EcalRecHitWorkerRecover::set(const edm::EventSetup& es) -{ - - es.get().get(laser); - es.get().get(caloTopology_); - ecalScale_.setEventSetup( es ); - es.get().get(pEcalMapping_); - ecalMapping_ = pEcalMapping_.product(); - // geometry... - es.get().get("EcalBarrel",pEBGeom_); - es.get().get(caloGeometry_); - es.get().get(chStatus_); - geo_ = caloGeometry_.product(); - ebGeom_ = pEBGeom_.product(); - es.get().get(ttMap_); - recoveredDetIds_EB_.clear(); - recoveredDetIds_EE_.clear(); - tpgscale_.setEventSetup(es); +void EcalRecHitWorkerRecover::set(const edm::EventSetup& es) { + es.get().get(laser); + es.get().get(caloTopology_); + ecalScale_.setEventSetup(es); + es.get().get(pEcalMapping_); + ecalMapping_ = pEcalMapping_.product(); + // geometry... + es.get().get("EcalBarrel", pEBGeom_); + es.get().get(caloGeometry_); + es.get().get(chStatus_); + geo_ = caloGeometry_.product(); + ebGeom_ = pEBGeom_.product(); + es.get().get(ttMap_); + recoveredDetIds_EB_.clear(); + recoveredDetIds_EE_.clear(); + tpgscale_.setEventSetup(es); } - -bool -EcalRecHitWorkerRecover::run( const edm::Event & evt, - const EcalUncalibratedRecHit& uncalibRH, - EcalRecHitCollection & result ) -{ - DetId detId=uncalibRH.id(); - uint32_t flags = (0xF & uncalibRH.flags()); - - // get laser coefficient - //float lasercalib = laser->getLaserCorrection( detId, evt.time()); - - // killDeadChannels_ = true, means explicitely kill dead channels even if the recovered energies are computed in the code - // if you don't want to store the recovered energies in the rechit you can produce LogWarnings if logWarningEtThreshold_EB(EE)_FE>0 - // logWarningEtThreshold_EB(EE)_FE_<0 will not compute the recovered energies at all (faster) - - if ( killDeadChannels_ ) { - if ( (flags == EcalRecHitWorkerRecover::EB_single && !recoverEBIsolatedChannels_) - || (flags == EcalRecHitWorkerRecover::EE_single && !recoverEEIsolatedChannels_) - || (flags == EcalRecHitWorkerRecover::EB_VFE && !recoverEBVFE_) - || (flags == EcalRecHitWorkerRecover::EE_VFE && !recoverEEVFE_) - ) { - EcalRecHit hit( detId, 0., 0., EcalRecHit::kDead ); - hit.setFlag( EcalRecHit::kDead) ; - insertRecHit( hit, result); // insert trivial rechit with kDead flag - return true; - } - if ( flags == EcalRecHitWorkerRecover::EB_FE && !recoverEBFE_) { - EcalTrigTowerDetId ttDetId( ((EBDetId)detId).tower() ); - std::vector vid = ttMap_->constituentsOf( ttDetId ); - for ( std::vector::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) { - EcalRecHit hit( (*dit), 0., 0., EcalRecHit::kDead ); - hit.setFlag( EcalRecHit::kDead ) ; - insertRecHit( hit, result ); // insert trivial rechit with kDead flag - } - if(logWarningEtThreshold_EB_FE_<0)return true; // if you don't want log warning just return true - } - if ( flags == EcalRecHitWorkerRecover::EE_FE && !recoverEEFE_) { - EEDetId id( detId ); - EcalScDetId sc( 1+(id.ix()-1)/5, 1+(id.iy()-1)/5, id.zside() ); - std::vector eeC; - for(int dx=1; dx<=5; ++dx){ - for(int dy=1; dy<=5; ++dy){ - int ix = (sc.ix()-1)*5 + dx; - int iy = (sc.iy()-1)*5 + dy; - int iz = sc.zside(); - if(EEDetId::validDetId(ix, iy, iz)){ - eeC.push_back(EEDetId(ix, iy, iz)); - } - } - } - for ( size_t i = 0; i < eeC.size(); ++i ) { - EcalRecHit hit( eeC[i], 0., 0., EcalRecHit::kDead ); - hit.setFlag( EcalRecHit::kDead ) ; - insertRecHit( hit, result ); // insert trivial rechit with kDead flag - } - if(logWarningEtThreshold_EE_FE_<0) return true; // if you don't want log warning just return true - } +bool EcalRecHitWorkerRecover::run(const edm::Event& evt, + const EcalUncalibratedRecHit& uncalibRH, + EcalRecHitCollection& result) { + DetId detId = uncalibRH.id(); + uint32_t flags = (0xF & uncalibRH.flags()); + + // get laser coefficient + //float lasercalib = laser->getLaserCorrection( detId, evt.time()); + + // killDeadChannels_ = true, means explicitely kill dead channels even if the recovered energies are computed in the code + // if you don't want to store the recovered energies in the rechit you can produce LogWarnings if logWarningEtThreshold_EB(EE)_FE>0 + // logWarningEtThreshold_EB(EE)_FE_<0 will not compute the recovered energies at all (faster) + + if (killDeadChannels_) { + if ((flags == EcalRecHitWorkerRecover::EB_single && !recoverEBIsolatedChannels_) || + (flags == EcalRecHitWorkerRecover::EE_single && !recoverEEIsolatedChannels_) || + (flags == EcalRecHitWorkerRecover::EB_VFE && !recoverEBVFE_) || + (flags == EcalRecHitWorkerRecover::EE_VFE && !recoverEEVFE_)) { + EcalRecHit hit(detId, 0., 0., EcalRecHit::kDead); + hit.setFlag(EcalRecHit::kDead); + insertRecHit(hit, result); // insert trivial rechit with kDead flag + return true; + } + if (flags == EcalRecHitWorkerRecover::EB_FE && !recoverEBFE_) { + EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower()); + std::vector vid = ttMap_->constituentsOf(ttDetId); + for (std::vector::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) { + EcalRecHit hit((*dit), 0., 0., EcalRecHit::kDead); + hit.setFlag(EcalRecHit::kDead); + insertRecHit(hit, result); // insert trivial rechit with kDead flag + } + if (logWarningEtThreshold_EB_FE_ < 0) + return true; // if you don't want log warning just return true + } + if (flags == EcalRecHitWorkerRecover::EE_FE && !recoverEEFE_) { + EEDetId id(detId); + EcalScDetId sc(1 + (id.ix() - 1) / 5, 1 + (id.iy() - 1) / 5, id.zside()); + std::vector eeC; + for (int dx = 1; dx <= 5; ++dx) { + for (int dy = 1; dy <= 5; ++dy) { + int ix = (sc.ix() - 1) * 5 + dx; + int iy = (sc.iy() - 1) * 5 + dy; + int iz = sc.zside(); + if (EEDetId::validDetId(ix, iy, iz)) { + eeC.push_back(EEDetId(ix, iy, iz)); + } } + } + for (size_t i = 0; i < eeC.size(); ++i) { + EcalRecHit hit(eeC[i], 0., 0., EcalRecHit::kDead); + hit.setFlag(EcalRecHit::kDead); + insertRecHit(hit, result); // insert trivial rechit with kDead flag + } + if (logWarningEtThreshold_EE_FE_ < 0) + return true; // if you don't want log warning just return true + } + } - if ( flags == EcalRecHitWorkerRecover::EB_single ) { - // recover as single dead channel - ebDeadChannelCorrector.setCaloTopology(caloTopology_.product()); - - // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE - bool AcceptRecHit=true; - EcalRecHit hit = ebDeadChannelCorrector.correct( detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, &AcceptRecHit); - - if ( hit.energy() != 0 and AcceptRecHit == true ) { - hit.setFlag( EcalRecHit::kNeighboursRecovered ) ; - } else { - // recovery failed - hit.setFlag( EcalRecHit::kDead ) ; - } - insertRecHit( hit, result ); - - } else if ( flags == EcalRecHitWorkerRecover::EE_single ) { - // recover as single dead channel - eeDeadChannelCorrector.setCaloTopology(caloTopology_.product()); - - // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE - bool AcceptRecHit=true; - EcalRecHit hit = eeDeadChannelCorrector.correct( detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, &AcceptRecHit); - if ( hit.energy() != 0 and AcceptRecHit == true ) { - hit.setFlag( EcalRecHit::kNeighboursRecovered ) ; - } else { - // recovery failed - hit.setFlag( EcalRecHit::kDead ) ; - } - insertRecHit( hit, result ); - - } else if ( flags == EcalRecHitWorkerRecover::EB_VFE ) { - // recover as dead VFE - EcalRecHit hit( detId, 0., 0.); - hit.setFlag( EcalRecHit::kDead ) ; - // recovery not implemented - insertRecHit( hit, result ); - } else if ( flags == EcalRecHitWorkerRecover::EB_FE ) { - // recover as dead TT - - EcalTrigTowerDetId ttDetId( ((EBDetId)detId).tower() ); - edm::Handle pTPDigis; - evt.getByToken(tpDigiToken_, pTPDigis); - const EcalTrigPrimDigiCollection * tpDigis = nullptr; - tpDigis = pTPDigis.product(); - - EcalTrigPrimDigiCollection::const_iterator tp = tpDigis->find( ttDetId ); - // recover the whole trigger tower - if ( tp != tpDigis->end() ) { - //std::vector vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) ); - std::vector vid = ttMap_->constituentsOf( ttDetId ); - float tpEt = ecalScale_.getTPGInGeV( tp->compressedEt(), tp->id() ); - float tpEtThreshEB = logWarningEtThreshold_EB_FE_; - if(tpEt>tpEtThreshEB){ - edm::LogWarning("EnergyInDeadEB_FE")<<"TP energy in the dead TT = "<::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) { - if (alreadyInserted(*dit)) continue; - float theta = ebGeom_->getGeometry(*dit)->getPosition().theta(); - float tpEt = ecalScale_.getTPGInGeV( tp->compressedEt(), tp->id() ); - if(checkChannelStatus(*dit, dbStatusToBeExcludedEB_)){ - EcalRecHit hit( *dit, tpEt /((float)vid.size()) / sin(theta), 0.); - hit.setFlag( EcalRecHit::kTowerRecovered ) ; - if ( tp->compressedEt() == 0xFF ) hit.setFlag( EcalRecHit::kTPSaturated ); - if ( tp->sFGVB() ) hit.setFlag( EcalRecHit::kL1SpikeFlag ); - insertRecHit( hit, result ); - } - } - } else { - // tp not found => recovery failed - std::vector vid = ttMap_->constituentsOf( ttDetId ); - for ( std::vector::const_iterator dit = vid.begin(); dit != vid.end(); ++dit ) { - if (alreadyInserted(*dit)) continue; - EcalRecHit hit( *dit,0., 0. ); - hit.setFlag( EcalRecHit::kDead ) ; - insertRecHit( hit, result ); - } - } - } - } else if ( flags == EcalRecHitWorkerRecover::EE_FE ) { - // Structure for recovery: - // ** SC --> EEDetId constituents (eeC) --> associated Trigger Towers (aTT) --> EEDetId constituents (aTTC) - // ** energy for a SC EEDetId = [ sum_aTT(energy) - sum_aTTC(energy) ] / N_eeC - // .. i.e. the total energy of the TTs covering the SC minus - // .. the energy of the recHits in the TTs but not in the SC - //std::vector vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) ); - // due to lack of implementation of the EcalTrigTowerDetId ix,iy methods in EE we compute Et recovered energies (in EB we compute E) - - EEDetId eeId( detId ); - EcalScDetId sc( (eeId.ix()-1)/5+1, (eeId.iy()-1)/5+1, eeId.zside() ); - std::set eeC; - for(int dx=1; dx<=5; ++dx){ - for(int dy=1; dy<=5; ++dy){ - int ix = (sc.ix()-1)*5 + dx; - int iy = (sc.iy()-1)*5 + dy; - int iz = sc.zside(); - if(EEDetId::validDetId(ix, iy, iz)){ - EEDetId id(ix, iy, iz); - if (checkChannelStatus(id,dbStatusToBeExcludedEE_)){ - eeC.insert(id); - } // check status - } - } - } - - edm::Handle pTPDigis; - evt.getByToken(tpDigiToken_, pTPDigis); - const EcalTrigPrimDigiCollection * tpDigis = nullptr; - tpDigis = pTPDigis.product(); - - // associated trigger towers - std::set aTT; - for ( std::set::const_iterator it = eeC.begin(); it!=eeC.end(); ++it ) { - aTT.insert( ttMap_->towerOf( *it ) ); - } - // associated trigger towers: total energy - float totE = 0; - // associated trigger towers: EEDetId constituents - std::set aTTC; - bool atLeastOneTPSaturated = false; - for ( std::set::const_iterator it = aTT.begin(); it != aTT.end(); ++it ) { - // add the energy of this trigger tower - EcalTrigPrimDigiCollection::const_iterator itTP = tpDigis->find( *it ); - if ( itTP != tpDigis->end() ) { - - std::vector v = ttMap_->constituentsOf( *it ); - - // from the constituents, remove dead channels - std::vector::iterator ttcons = v.begin(); - while (ttcons != v.end()){ - if (!checkChannelStatus(*ttcons,dbStatusToBeExcludedEE_)){ - ttcons=v.erase(ttcons); - } else { - ++ttcons; - } - }// while - - if ( itTP->compressedEt() == 0xFF ){ // In the case of a saturated trigger tower, a fraction - atLeastOneTPSaturated = true; //of the saturated energy is put in: number of xtals in dead region/total xtals in TT *63.75 - - //Alternative recovery algorithm that I will now investigate. - //Estimate energy sums the energy in the working channels, then decides how much energy - //to put here depending on that. Duncan 20101203 - - totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v); - - /* + if (flags == EcalRecHitWorkerRecover::EB_single) { + // recover as single dead channel + ebDeadChannelCorrector.setCaloTopology(caloTopology_.product()); + + // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE + bool AcceptRecHit = true; + 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); + } else { + // recovery failed + hit.setFlag(EcalRecHit::kDead); + } + insertRecHit(hit, result); + + } else if (flags == EcalRecHitWorkerRecover::EE_single) { + // recover as single dead channel + eeDeadChannelCorrector.setCaloTopology(caloTopology_.product()); + + // channel recovery. Accepted new RecHit has the flag AcceptRecHit=TRUE + bool AcceptRecHit = true; + 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 { + // recovery failed + hit.setFlag(EcalRecHit::kDead); + } + insertRecHit(hit, result); + + } else if (flags == EcalRecHitWorkerRecover::EB_VFE) { + // recover as dead VFE + EcalRecHit hit(detId, 0., 0.); + hit.setFlag(EcalRecHit::kDead); + // recovery not implemented + insertRecHit(hit, result); + } else if (flags == EcalRecHitWorkerRecover::EB_FE) { + // recover as dead TT + + EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower()); + edm::Handle pTPDigis; + evt.getByToken(tpDigiToken_, pTPDigis); + const EcalTrigPrimDigiCollection* tpDigis = nullptr; + tpDigis = pTPDigis.product(); + + EcalTrigPrimDigiCollection::const_iterator tp = tpDigis->find(ttDetId); + // recover the whole trigger tower + if (tp != tpDigis->end()) { + //std::vector vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) ); + std::vector vid = ttMap_->constituentsOf(ttDetId); + float tpEt = ecalScale_.getTPGInGeV(tp->compressedEt(), tp->id()); + float tpEtThreshEB = logWarningEtThreshold_EB_FE_; + if (tpEt > tpEtThreshEB) { + edm::LogWarning("EnergyInDeadEB_FE") << "TP energy in the dead TT = " << tpEt << " at " << ttDetId; + } + if (!killDeadChannels_ || recoverEBFE_) { + // democratic energy sharing + + for (std::vector::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) { + if (alreadyInserted(*dit)) + continue; + float theta = ebGeom_->getGeometry(*dit)->getPosition().theta(); + float tpEt = ecalScale_.getTPGInGeV(tp->compressedEt(), tp->id()); + if (checkChannelStatus(*dit, dbStatusToBeExcludedEB_)) { + EcalRecHit hit(*dit, tpEt / ((float)vid.size()) / sin(theta), 0.); + hit.setFlag(EcalRecHit::kTowerRecovered); + if (tp->compressedEt() == 0xFF) + hit.setFlag(EcalRecHit::kTPSaturated); + if (tp->sFGVB()) + hit.setFlag(EcalRecHit::kL1SpikeFlag); + insertRecHit(hit, result); + } + } + } else { + // tp not found => recovery failed + std::vector vid = ttMap_->constituentsOf(ttDetId); + for (std::vector::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) { + if (alreadyInserted(*dit)) + continue; + EcalRecHit hit(*dit, 0., 0.); + hit.setFlag(EcalRecHit::kDead); + insertRecHit(hit, result); + } + } + } + } else if (flags == EcalRecHitWorkerRecover::EE_FE) { + // Structure for recovery: + // ** SC --> EEDetId constituents (eeC) --> associated Trigger Towers (aTT) --> EEDetId constituents (aTTC) + // ** energy for a SC EEDetId = [ sum_aTT(energy) - sum_aTTC(energy) ] / N_eeC + // .. i.e. the total energy of the TTs covering the SC minus + // .. the energy of the recHits in the TTs but not in the SC + //std::vector vid = ecalMapping_->dccTowerConstituents( ecalMapping_->DCCid( ttDetId ), ecalMapping_->iTT( ttDetId ) ); + // due to lack of implementation of the EcalTrigTowerDetId ix,iy methods in EE we compute Et recovered energies (in EB we compute E) + + EEDetId eeId(detId); + EcalScDetId sc((eeId.ix() - 1) / 5 + 1, (eeId.iy() - 1) / 5 + 1, eeId.zside()); + std::set eeC; + for (int dx = 1; dx <= 5; ++dx) { + for (int dy = 1; dy <= 5; ++dy) { + int ix = (sc.ix() - 1) * 5 + dx; + int iy = (sc.iy() - 1) * 5 + dy; + int iz = sc.zside(); + if (EEDetId::validDetId(ix, iy, iz)) { + EEDetId id(ix, iy, iz); + if (checkChannelStatus(id, dbStatusToBeExcludedEE_)) { + eeC.insert(id); + } // check status + } + } + } + + edm::Handle pTPDigis; + evt.getByToken(tpDigiToken_, pTPDigis); + const EcalTrigPrimDigiCollection* tpDigis = nullptr; + tpDigis = pTPDigis.product(); + + // associated trigger towers + std::set aTT; + for (std::set::const_iterator it = eeC.begin(); it != eeC.end(); ++it) { + aTT.insert(ttMap_->towerOf(*it)); + } + // associated trigger towers: total energy + float totE = 0; + // associated trigger towers: EEDetId constituents + std::set aTTC; + bool atLeastOneTPSaturated = false; + for (std::set::const_iterator it = aTT.begin(); it != aTT.end(); ++it) { + // add the energy of this trigger tower + EcalTrigPrimDigiCollection::const_iterator itTP = tpDigis->find(*it); + if (itTP != tpDigis->end()) { + std::vector v = ttMap_->constituentsOf(*it); + + // from the constituents, remove dead channels + std::vector::iterator ttcons = v.begin(); + while (ttcons != v.end()) { + if (!checkChannelStatus(*ttcons, dbStatusToBeExcludedEE_)) { + ttcons = v.erase(ttcons); + } else { + ++ttcons; + } + } // while + + if (itTP->compressedEt() == 0xFF) { // In the case of a saturated trigger tower, a fraction + atLeastOneTPSaturated = + true; //of the saturated energy is put in: number of xtals in dead region/total xtals in TT *63.75 + + //Alternative recovery algorithm that I will now investigate. + //Estimate energy sums the energy in the working channels, then decides how much energy + //to put here depending on that. Duncan 20101203 + + totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v); + + /* These commented out lines use 64GeV*fraction of the TT overlapping the dead FE @@ -286,171 +291,155 @@ EcalRecHitWorkerRecover::run( const edm::Event & evt, } //std::cout << count << ", " << v.size() << std::endl; totE+=((float)count/(float)v.size())* ((it->ietaAbs()>26)?2*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ):ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() ));*/ - } - else {totE += ((it->ietaAbs()>26)?2:1)*ecalScale_.getTPGInGeV( itTP->compressedEt(), itTP->id() );} - - - // get the trigger tower constituents - - if (itTP->compressedEt() == 0){ // If there's no energy in TT, the constituents are removed from the recovery. - for (size_t i = 0 ; i < v.size(); ++i) - eeC.erase(v[i]); - } - else if (itTP->compressedEt()!=0xFF){ //If it's saturated the energy has already been determined, so we do not want to subtract any channels - for ( size_t j = 0; j < v.size(); ++j ) { - aTTC.insert( v[j] ); - } - } - - } - } - // remove crystals of dead SC - // (this step is not needed if sure that SC crystals are not - // in the recHit collection) - - for ( std::set::const_iterator it = eeC.begin(); it != eeC.end(); ++it ) { - aTTC.erase(*it); - } - // compute the total energy for the dead SC - const EcalRecHitCollection * hits = &result; - for ( std::set::const_iterator it = aTTC.begin(); it != aTTC.end(); ++it ) { - EcalRecHitCollection::const_iterator jt = hits->find( *it ); - if ( jt != hits->end() ) { - float energy = jt->energy(); // Correct conversion to Et - float eta = geo_->getPosition(jt->id()).eta(); - float pf = 1.0/cosh(eta); - // use Et instead of E, consistent with the Et estimation of the associated TT - totE -= energy*pf; - } - } - - - float scEt = totE; - float scEtThreshEE = logWarningEtThreshold_EE_FE_; - if(scEt>scEtThreshEE){ - edm::LogWarning("EnergyInDeadEE_FE")<<"TP energy in the dead TT = "<::const_iterator it = eeC.begin(); it != eeC.end(); ++it ) { - - float eta = geo_->getPosition(*it).eta(); //Convert back to E from Et for the recovered hits - float pf = 1.0/cosh(eta); - EcalRecHit hit( *it, totE / ((float)eeC.size()*pf), 0); - - if (atLeastOneTPSaturated) hit.setFlag(EcalRecHit::kTPSaturated ); - hit.setFlag(EcalRecHit::kTowerRecovered); - insertRecHit( hit, result ); - - }// for - }// if + } else { + totE += ((it->ietaAbs() > 26) ? 2 : 1) * ecalScale_.getTPGInGeV(itTP->compressedEt(), itTP->id()); } - return true; -} -float EcalRecHitWorkerRecover::estimateEnergy(int ieta, EcalRecHitCollection* hits, const std::set& sId, const std::vector& vId ) { - - float xtalE=0; - int count = 0; - for (std::vector::const_iterator vIdit = vId.begin(); vIdit != vId.end(); ++ vIdit){ - std::set::const_iterator sIdit = sId.find(*vIdit); - if (sIdit==sId.end()){ - float energy = hits->find(*vIdit)->energy(); - float eta = geo_->getPosition(*vIdit).eta(); - float pf = 1.0/cosh(eta); - xtalE += energy*pf; - count++; - } - } - - if (count==0) { // If there are no overlapping crystals return saturated value. - - double etsat = tpgscale_.getTPGInGeV(0xFF, - ttMap_->towerOf(*vId.begin())); // get saturation value for the first - // constituent, for the others it's the same - - return etsat/cosh(ieta)*(ieta>26?2:1); // account for duplicated TT in EE for ieta>26 - } - else return xtalE*((vId.size()/(float)count) - 1)*(ieta>26?2:1); - - + // get the trigger tower constituents + + if (itTP->compressedEt() == 0) { // If there's no energy in TT, the constituents are removed from the recovery. + for (size_t i = 0; i < v.size(); ++i) + eeC.erase(v[i]); + } else if (itTP->compressedEt() != 0xFF) { + //If it's saturated the energy has already been determined, so we do not want to subtract any channels + for (size_t j = 0; j < v.size(); ++j) { + aTTC.insert(v[j]); + } + } + } + } + // remove crystals of dead SC + // (this step is not needed if sure that SC crystals are not + // in the recHit collection) + + for (std::set::const_iterator it = eeC.begin(); it != eeC.end(); ++it) { + aTTC.erase(*it); + } + // compute the total energy for the dead SC + const EcalRecHitCollection* hits = &result; + for (std::set::const_iterator it = aTTC.begin(); it != aTTC.end(); ++it) { + EcalRecHitCollection::const_iterator jt = hits->find(*it); + if (jt != hits->end()) { + float energy = jt->energy(); // Correct conversion to Et + float eta = geo_->getPosition(jt->id()).eta(); + float pf = 1.0 / cosh(eta); + // use Et instead of E, consistent with the Et estimation of the associated TT + totE -= energy * pf; + } + } + + float scEt = totE; + float scEtThreshEE = logWarningEtThreshold_EE_FE_; + if (scEt > scEtThreshEE) { + edm::LogWarning("EnergyInDeadEE_FE") << "TP energy in the dead TT = " << scEt << " at " << sc; + } + + // assign the energy to the SC crystals + if (!killDeadChannels_ || recoverEEFE_) { // if eeC is empty, i.e. there are no hits + // in the tower, nothing is returned. No negative values from noise. + for (std::set::const_iterator it = eeC.begin(); it != eeC.end(); ++it) { + float eta = geo_->getPosition(*it).eta(); //Convert back to E from Et for the recovered hits + float pf = 1.0 / cosh(eta); + EcalRecHit hit(*it, totE / ((float)eeC.size() * pf), 0); + + if (atLeastOneTPSaturated) + hit.setFlag(EcalRecHit::kTPSaturated); + hit.setFlag(EcalRecHit::kTowerRecovered); + insertRecHit(hit, result); + + } // for + } // if + } + return true; } +float EcalRecHitWorkerRecover::estimateEnergy(int ieta, + EcalRecHitCollection* hits, + const std::set& sId, + const std::vector& vId) { + float xtalE = 0; + int count = 0; + for (std::vector::const_iterator vIdit = vId.begin(); vIdit != vId.end(); ++vIdit) { + std::set::const_iterator sIdit = sId.find(*vIdit); + if (sIdit == sId.end()) { + float energy = hits->find(*vIdit)->energy(); + float eta = geo_->getPosition(*vIdit).eta(); + float pf = 1.0 / cosh(eta); + xtalE += energy * pf; + count++; + } + } -void EcalRecHitWorkerRecover::insertRecHit( const EcalRecHit &hit, EcalRecHitCollection &collection ) -{ - // skip already inserted DetId's and raise a log warning - if ( alreadyInserted( hit.id() ) ) { - edm::LogWarning("EcalRecHitWorkerRecover") << "DetId already recovered! Skipping..."; - return; - } - EcalRecHitCollection::iterator it = collection.find( hit.id() ); - if ( it == collection.end() ) { - // insert the hit in the collection - collection.push_back( hit ); - } else { - // overwrite existing recHit - *it = hit; - } - if ( hit.id().subdetId() == EcalBarrel ) { - recoveredDetIds_EB_.insert( hit.id() ); - } else if ( hit.id().subdetId() == EcalEndcap ) { - recoveredDetIds_EE_.insert( hit.id() ); - } else { - edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << hit.id().rawId(); - } + if (count == 0) { // If there are no overlapping crystals return saturated value. + + double etsat = tpgscale_.getTPGInGeV(0xFF, + ttMap_->towerOf(*vId.begin())); // get saturation value for the first + // constituent, for the others it's the same + + return etsat / cosh(ieta) * (ieta > 26 ? 2 : 1); // account for duplicated TT in EE for ieta>26 + } else + return xtalE * ((vId.size() / (float)count) - 1) * (ieta > 26 ? 2 : 1); } +void EcalRecHitWorkerRecover::insertRecHit(const EcalRecHit& hit, EcalRecHitCollection& collection) { + // skip already inserted DetId's and raise a log warning + if (alreadyInserted(hit.id())) { + edm::LogWarning("EcalRecHitWorkerRecover") << "DetId already recovered! Skipping..."; + return; + } + EcalRecHitCollection::iterator it = collection.find(hit.id()); + if (it == collection.end()) { + // insert the hit in the collection + collection.push_back(hit); + } else { + // overwrite existing recHit + *it = hit; + } + if (hit.id().subdetId() == EcalBarrel) { + recoveredDetIds_EB_.insert(hit.id()); + } else if (hit.id().subdetId() == EcalEndcap) { + recoveredDetIds_EE_.insert(hit.id()); + } else { + edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << hit.id().rawId(); + } +} -bool EcalRecHitWorkerRecover::alreadyInserted( const DetId & id ) -{ - bool res = false; - if ( id.subdetId() == EcalBarrel ) { - res = ( recoveredDetIds_EB_.find( id ) != recoveredDetIds_EB_.end() ); - } else if ( id.subdetId() == EcalEndcap ) { - res = ( recoveredDetIds_EE_.find( id ) != recoveredDetIds_EE_.end() ); - } else { - edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << id.rawId(); - } - return res; +bool EcalRecHitWorkerRecover::alreadyInserted(const DetId& id) { + bool res = false; + if (id.subdetId() == EcalBarrel) { + res = (recoveredDetIds_EB_.find(id) != recoveredDetIds_EB_.end()); + } else if (id.subdetId() == EcalEndcap) { + res = (recoveredDetIds_EE_.find(id) != recoveredDetIds_EE_.end()); + } else { + edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << id.rawId(); + } + return res; } // In the future, this will be used to calibrate the TT energy. There is a dependance on // eta at lower energies that can be corrected for here after more validation. -float EcalRecHitWorkerRecover::recCheckCalib(float eTT, int ieta){ - - return eTT; - -} +float EcalRecHitWorkerRecover::recCheckCalib(float eTT, int ieta) { return eTT; } // return false is the channel has status in the list of statusestoexclude // true otherwise (channel ok) // Careful: this function works on raw (encoded) channel statuses -bool EcalRecHitWorkerRecover::checkChannelStatus(const DetId& id, - const std::vector& statusestoexclude){ - - - if (!chStatus_.isValid()) - edm::LogError("ObjectNotFound") << "Channel Status not set"; - - - EcalChannelStatus::const_iterator chIt = chStatus_->find( id ); +bool EcalRecHitWorkerRecover::checkChannelStatus(const DetId& id, const std::vector& statusestoexclude) { + if (!chStatus_.isValid()) + edm::LogError("ObjectNotFound") << "Channel Status not set"; + + EcalChannelStatus::const_iterator chIt = chStatus_->find(id); uint16_t dbStatus = 0; - if ( chIt != chStatus_->end() ) { + if (chIt != chStatus_->end()) { dbStatus = chIt->getEncodedStatusCode(); } else { - edm::LogError("ObjectNotFound") << "No channel status found for xtal " - << id.rawId() - << "! something wrong with EcalChannelStatus in your DB? "; + edm::LogError("ObjectNotFound") << "No channel status found for xtal " << id.rawId() + << "! something wrong with EcalChannelStatus in your DB? "; } - - for (std::vector::const_iterator status = statusestoexclude.begin(); - status!= statusestoexclude.end(); ++status){ - - if ( *status == dbStatus) return false; - + + for (std::vector::const_iterator status = statusestoexclude.begin(); status != statusestoexclude.end(); + ++status) { + if (*status == dbStatus) + return false; } return true; @@ -458,4 +447,4 @@ bool EcalRecHitWorkerRecover::checkChannelStatus(const DetId& id, #include "FWCore/Framework/interface/MakerMacros.h" #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h" -DEFINE_EDM_PLUGIN( EcalRecHitWorkerFactory, EcalRecHitWorkerRecover, "EcalRecHitWorkerRecover" ); +DEFINE_EDM_PLUGIN(EcalRecHitWorkerFactory, EcalRecHitWorkerRecover, "EcalRecHitWorkerRecover"); diff --git a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h index 474a13878e8d2..a73d209792f7e 100644 --- a/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h +++ b/RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h @@ -28,72 +28,67 @@ #include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h" class EcalRecHitWorkerRecover : public EcalRecHitWorkerBaseClass { - public: - EcalRecHitWorkerRecover(const edm::ParameterSet&, edm::ConsumesCollector& c); - ~EcalRecHitWorkerRecover() override {}; - - void set(const edm::EventSetup& es) override; - bool run(const edm::Event& evt, const EcalUncalibratedRecHit& uncalibRH, EcalRecHitCollection & result) override; - - protected: - - void insertRecHit( const EcalRecHit &hit, EcalRecHitCollection &collection ); - float recCheckCalib(float energy, int ieta); - bool alreadyInserted( const DetId & id ); - float estimateEnergy(int ieta, EcalRecHitCollection* hits, - const std::set& sId, - const std::vector& vId); - bool checkChannelStatus(const DetId& id, - const std::vector& statusestoexclude); - - edm::ESHandle laser; - - // isolated dead channels - edm::ESHandle caloTopology_; - edm::ESHandle caloGeometry_; - edm::ESHandle chStatus_; - - - double singleRecoveryThreshold_; - std::string singleRecoveryMethod_; - bool killDeadChannels_; - - bool recoverEBIsolatedChannels_; - bool recoverEEIsolatedChannels_; - bool recoverEBVFE_; - bool recoverEEVFE_; - bool recoverEBFE_; - bool recoverEEFE_; - - // list of channel statuses for which recovery in EE should - // not be attempted - std::vector dbStatusToBeExcludedEE_; - std::vector dbStatusToBeExcludedEB_; - - // dead FE - EcalTPGScale ecalScale_; - edm::EDGetTokenT tpDigiToken_; - edm::ESHandle< EcalElectronicsMapping > pEcalMapping_; - const EcalElectronicsMapping *ecalMapping_; - double logWarningEtThreshold_EB_FE_; - double logWarningEtThreshold_EE_FE_; - - edm::ESHandle ttMap_; - - edm::ESHandle pEBGeom_; - const CaloSubdetectorGeometry * ebGeom_; - const CaloGeometry* geo_; - - std::unique_ptr rechitMaker_; - - std::set recoveredDetIds_EB_; - std::set recoveredDetIds_EE_; - - EcalTPGScale tpgscale_; - - EcalDeadChannelRecoveryAlgos ebDeadChannelCorrector; - EcalDeadChannelRecoveryAlgos eeDeadChannelCorrector; - +public: + EcalRecHitWorkerRecover(const edm::ParameterSet&, edm::ConsumesCollector& c); + ~EcalRecHitWorkerRecover() override{}; + + void set(const edm::EventSetup& es) override; + bool run(const edm::Event& evt, const EcalUncalibratedRecHit& uncalibRH, EcalRecHitCollection& result) override; + +protected: + void insertRecHit(const EcalRecHit& hit, EcalRecHitCollection& collection); + float recCheckCalib(float energy, int ieta); + bool alreadyInserted(const DetId& id); + float estimateEnergy(int ieta, EcalRecHitCollection* hits, const std::set& sId, const std::vector& vId); + bool checkChannelStatus(const DetId& id, const std::vector& statusestoexclude); + + edm::ESHandle laser; + + // isolated dead channels + edm::ESHandle caloTopology_; + edm::ESHandle caloGeometry_; + edm::ESHandle chStatus_; + + double singleRecoveryThreshold_; + double sum8RecoveryThreshold_; + std::string singleRecoveryMethod_; + bool killDeadChannels_; + + bool recoverEBIsolatedChannels_; + bool recoverEEIsolatedChannels_; + bool recoverEBVFE_; + bool recoverEEVFE_; + bool recoverEBFE_; + bool recoverEEFE_; + + // list of channel statuses for which recovery in EE should + // not be attempted + std::vector dbStatusToBeExcludedEE_; + std::vector dbStatusToBeExcludedEB_; + + // dead FE + EcalTPGScale ecalScale_; + edm::EDGetTokenT tpDigiToken_; + edm::ESHandle pEcalMapping_; + const EcalElectronicsMapping* ecalMapping_; + double logWarningEtThreshold_EB_FE_; + double logWarningEtThreshold_EE_FE_; + + edm::ESHandle ttMap_; + + edm::ESHandle pEBGeom_; + const CaloSubdetectorGeometry* ebGeom_; + const CaloGeometry* geo_; + + std::unique_ptr rechitMaker_; + + std::set recoveredDetIds_EB_; + std::set recoveredDetIds_EE_; + + EcalTPGScale tpgscale_; + + EcalDeadChannelRecoveryAlgos ebDeadChannelCorrector; + EcalDeadChannelRecoveryAlgos eeDeadChannelCorrector; }; #endif diff --git a/RecoLocalCalo/EcalRecProducers/python/ecalRecHit_cfi.py b/RecoLocalCalo/EcalRecProducers/python/ecalRecHit_cfi.py index 21e91ceadfe0c..cab0a08e6ad4a 100644 --- a/RecoLocalCalo/EcalRecProducers/python/ecalRecHit_cfi.py +++ b/RecoLocalCalo/EcalRecProducers/python/ecalRecHit_cfi.py @@ -49,7 +49,7 @@ # for channel recovery algoRecover = cms.string("EcalRecHitWorkerRecover"), - recoverEBIsolatedChannels = cms.bool(False), + recoverEBIsolatedChannels = cms.bool(True),##default is false recoverEEIsolatedChannels = cms.bool(False), recoverEBVFE = cms.bool(False), recoverEEVFE = cms.bool(False), @@ -76,8 +76,11 @@ eeDetIdToBeRecovered = cms.InputTag("ecalDetIdToBeRecovered:eeDetId"), ebFEToBeRecovered = cms.InputTag("ecalDetIdToBeRecovered:ebFE"), eeFEToBeRecovered = cms.InputTag("ecalDetIdToBeRecovered:eeFE"), - singleChannelRecoveryMethod = cms.string("NeuralNetworks"), - singleChannelRecoveryThreshold = cms.double(8), + singleChannelRecoveryMethod = cms.string("BDTG"), + singleChannelRecoveryThreshold = cms.double(0.70), #Threshold in GeV + sum8ChannelRecoveryThreshold = cms.double(0.), #Threshold in GeV + bdtWeightFileNoCracks = cms.FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/bdtgAllRH_8GT700MeV_noCracks_ZskimData2017_v1.xml"), + bdtWeightFileCracks = cms.FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/bdtgAllRH_8GT700MeV_onlyCracks_ZskimData2017_v1.xml"), triggerPrimitiveDigiCollection = cms.InputTag("ecalDigis:EcalTriggerPrimitives"), cleaningConfig=cleaningAlgoConfig, @@ -88,4 +91,6 @@ fastSim.toModify(ecalRecHit, killDeadChannels = False, recoverEBFE = False, - recoverEEFE = False) + recoverEEFE = False, + recoverEBIsolatedChannels = False + )