METFilters/plugins/EcalBadCalibFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalBadCalibFilter
0004 // Class:      EcalBadCalibFilter
0005 //
0006 /**\class EcalBadCalibFilter EcalBadCalibFilter.cc
0007
0008  Description: <one line class summary>
0009  Event filtering to remove events with anomalous energy intercalibrations in specific ECAL channels
0010 */
0011 //
0012 // Original Authors:  D. Petyt
0013 //
0014
0015 // include files
0016
0017 #include <iostream>
0018
0019 #include "FWCore/Framework/interface/Frameworkfwd.h"
0020 #include "FWCore/Framework/interface/global/EDFilter.h"
0021 #include "FWCore/Framework/interface/Event.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0026 #include "DataFormats/DetId/interface/DetId.h"
0027 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0028 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0029 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0030
0031 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0032 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0033
0034 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0035
0036 using namespace std;
0037
0038 class EcalBadCalibFilter : public edm::global::EDFilter<> {
0039 public:
0040   explicit EcalBadCalibFilter(const edm::ParameterSet& iConfig);
0041   ~EcalBadCalibFilter() override {}
0042
0043 private:
0044   // main filter function
0045
0046   bool filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const override;
0047
0048   // input parameters
0049   // ecal rechit collection (from AOD)
0050   const edm::EDGetTokenT<EcalRecHitCollection> ecalRHSrcToken_;
0051   const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geoToken_;
0052
0053   //config parameters (defining the cuts on the bad SCs)
0054   const double ecalMin_;  // ecal rechit et threshold
0055
0056   const std::vector<unsigned int> baddetEcal_;  // DetIds of bad Ecal channels
0057
0058   const bool taggingMode_;
0059   const bool debug_;  // prints out debug info if set to true
0060 };
0061
0062 // read the parameters from the config file
0063 EcalBadCalibFilter::EcalBadCalibFilter(const edm::ParameterSet& iConfig)
0064     : ecalRHSrcToken_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("EcalRecHitSource"))),
0065       geoToken_(esConsumes()),
0066       ecalMin_(iConfig.getParameter<double>("ecalMinEt")),
0067       baddetEcal_(iConfig.getParameter<std::vector<unsigned int> >("baddetEcal")),
0068       taggingMode_(iConfig.getParameter<bool>("taggingMode")),
0069       debug_(iConfig.getParameter<bool>("debug")) {
0070   produces<bool>();
0071 }
0072
0073 bool EcalBadCalibFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0074   // load required collections
0075
0076   // Ecal rechit collection
0077   edm::Handle<EcalRecHitCollection> ecalRHs;
0078   iEvent.getByToken(ecalRHSrcToken_, ecalRHs);
0079
0080   // Calo Geometry - needed for computing E_t
0081   const CaloGeometry& geo = iSetup.getData(geoToken_);
0082
0083   // by default the event is OK
0084   bool pass = true;
0085
0086   for (const auto ecalit : baddetEcal_) {
0087     DetId ecaldet(ecalit);
0088
0089     if (ecaldet.rawId() == 0)
0090       continue;
0091
0092     // find rechit corresponding to this DetId
0093     EcalRecHitCollection::const_iterator ecalhit = ecalRHs->find(ecaldet);
0094
0095     if (ecalhit == ecalRHs->end())
0096       continue;
0097
0098     // if rechit not found, move to next DetId
0099     if (ecalhit->id().rawId() == 0 || ecalhit->id().rawId() != ecaldet.rawId()) {
0100       continue;
0101     }
0102
0103     // define energy variables
0104     float ene = 0;
0105     float et = 0;
0106
0107     // rechit has been found: obtain crystal energy
0108     ene = ecalhit->energy();
0109
0110     // compute transverse energy
0111     const GlobalPoint& posecal = geo.getPosition(ecaldet);
0112     float pf = posecal.perp() / posecal.mag();
0113     et = ene * pf;
0114
0115     // print some debug info
0116     if (debug_) {
0117       // ref: DataFormats/EcalDetId/interface/EcalSubdetector.h
0118       // EcalBarrel
0119       if (ecaldet.subdetId() == 1) {
0120         EBDetId ebdet(ecalit);
0121         int ix = ebdet.ieta();
0122         int iy = ebdet.iphi();
0123         int iz = ebdet.zside();
0124
0125         edm::LogInfo("EcalBadCalibFilter") << "DetId=" << ecaldet.rawId();
0126         edm::LogInfo("EcalBadCalibFilter") << "ieta=" << ix << " iphi=" << iy << " iz=" << iz;
0127         edm::LogInfo("EcalBadCalibFilter") << "Et=" << et << " thresh=" << ecalMin_;
0128       }
0129
0130       // EcalEndcap
0131       else if (ecaldet.subdetId() == 2) {
0132         EEDetId eedet(ecalit);
0133         int ix = eedet.ix();
0134         int iy = eedet.iy();
0135         int iz = eedet.zside();
0136
0137         edm::LogInfo("EcalBadCalibFilter") << "DetId=" << ecaldet.rawId();
0138         edm::LogInfo("EcalBadCalibFilter") << "ix=" << ix << " iy=" << iy << " iz=" << iz;
0139         edm::LogInfo("EcalBadCalibFilter") << "Et=" << et << " thresh=" << ecalMin_;
0140       }
0141     }
0142
0143     // if transverse energy is above threshold and channel has bad IC
0144     if (et > ecalMin_) {
0145       pass = false;
0146       if (debug_) {
0147         edm::LogInfo("EcalBadCalibFilter") << "DUMP EVENT" << std::endl;
0148       }
0149     }
0150   }
0151
0152   // print the decision if event is bad
0153   if (pass == false && debug_)
0154     edm::LogInfo("EcalBadCalibFilter") << "REJECT EVENT!!!";
0155
0156   iEvent.put(std::make_unique<bool>(pass));
0157
0158   return taggingMode_ || pass;
0159 }
0160
0161 #include "FWCore/Framework/interface/MakerMacros.h"
0162
0163 DEFINE_FWK_MODULE(EcalBadCalibFilter);