Dead channels 11 0 x

RecoLocalCalo/EcalDeadChannelRecoveryAlgos/BuildFile.xml

@@ -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"/>

RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryAlgos.h

@@ -7,16 +7,23 @@
 #include "DataFormats/EcalDetId/interface/EBDetId.h"
 #include <string>
 
+#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
+
 #include "RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryNN.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);
+template <typename DetIdT> class EcalDeadChannelRecoveryAlgos {
+ public:
+  void setParameters(const edm::ParameterSet&ps);
+  void setCaloTopology(std::string algo, const CaloTopology *topology);
+  EcalRecHit correct(const DetIdT id,
+                     const EcalRecHitCollection &hit_collection,
+                     std::string algo, double single8Cut, double sum8Cut, bool *accFlag);
 
-private:
-  EcalDeadChannelRecoveryNN<DetIdT> nn;
+ private:
+  EcalDeadChannelRecoveryNN<DetIdT> nn_;
+  EcalDeadChannelRecoveryBDTG<DetIdT> bdtg_;
 };
-#endif  // RecoLocalCalo_EcalDeadChannelRecoveryAlgos_EcalDeadChannelRecoveryAlgos_HH
+#endif

RecoLocalCalo/EcalDeadChannelRecoveryAlgos/interface/EcalDeadChannelRecoveryBDTG.h

@@ -0,0 +1,54 @@
+#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_;
+
+  TMVA::Reader *readerNoCrack;
+  TMVA::Reader *readerCrack;
+
+};
+
+
+
+#endif

RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryAlgos.cc

@@ -8,30 +8,52 @@
 //      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>
+void EcalDeadChannelRecoveryAlgos<T>::setCaloTopology( std::string algo, 
+    const CaloTopology *topo) {
+  if (algo=="BDTG"){
+    bdtg_.setCaloTopology(topo);
+  }else{
+    nn_.setCaloTopology(topo);
+  }
+
+}
+
+
 template <typename T>
 EcalRecHit EcalDeadChannelRecoveryAlgos<T>::correct(
-    const T id, const EcalRecHitCollection &hit_collection, std::string algo, double Sum8Cut, bool *AcceptFlag) {
+    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;
-
+  double newEnergy = 0.0;
   if (algo == "NeuralNetworks") {
-    NewEnergy = this->nn.recover(id, hit_collection, Sum8Cut, AcceptFlag);
-  } else {
-    edm::LogError("EcalDeadChannelRecoveryAlgos") << "Invalid algorithm for dead channel recovery.";
-    *AcceptFlag = false;
+
+     newEnergy = this->nn_.recover(id, hit_collection, sum8Cut, acceptFlag);
+  }else 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;
   }
 
   uint32_t flag = 0;
-  return EcalRecHit(id, NewEnergy, 0, flag);
+  return EcalRecHit(id, newEnergy, 0, flag);
+
 }
 
 template class EcalDeadChannelRecoveryAlgos<EBDetId>;

RecoLocalCalo/EcalDeadChannelRecoveryAlgos/src/EcalDeadChannelRecoveryBDTG.cc

@@ -0,0 +1,171 @@
+//  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 = 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 <typename T> EcalDeadChannelRecoveryBDTG<T>::EcalDeadChannelRecoveryBDTG () {
+}
+
+
+template <typename T>EcalDeadChannelRecoveryBDTG<T>::~EcalDeadChannelRecoveryBDTG() {
+}
+
+template <>
+void EcalDeadChannelRecoveryBDTG<EBDetId>::setParameters(const edm::ParameterSet&ps)
+{
+
+  bdtWeightFileNoCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileNoCracks");
+  bdtWeightFileCracks_ = ps.getParameter<edm::FileInPath>("bdtWeightFileCracks");
+
+  this->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);
+  // do not apply the recovery if the matrix is not 3x3
+  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() ){
+	//keep the en, iphi, ieta of xtals of the matrix
+	if ( 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 { // 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.