Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:26:44

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalDeadCellBoundaryEnergyFilter
0004 // Class:      EcalDeadCellBoundaryEnergyFilter
0005 //
0006 /**\class EcalDeadCellBoundaryEnergyFilter EcalDeadCellBoundaryEnergyFilter.cc PhysicsTools/EcalDeadCellBoundaryEnergyFilter/src/EcalDeadCellBoundaryEnergyFilter.cc
0007 
0008  Description: <one line class summary>
0009  Event filtering for anomalous ECAL events where the energy measured by ECAL is significantly biased due to energy depositions
0010  in passive or problematic detector regions. The filter will handle ECAL flags and will compute the boundary energy in the channels
0011  surrounding the problematic regions such as dead channels and gaps.
0012 
0013  // Filter Algos :
0014  // a)  "TuningMode" keep all events and save event info in a ROOT TTree for tuning/algo development
0015  // b)  "FilterMode" returns false for all events passing the AnomalousEcalVariables.isDeadEcalCluster() function (--->rejects events affected by energy deposits in Dead Cells)
0016 
0017  Implementation:
0018  <Notes on implementation>
0019  */
0020 //
0021 // Original Author:  Konstantinos Theofilatos, Ulla Gebbert and Christian Sander
0022 //         Created:  Sat Nov 14 18:43:21 CET 2009
0023 //
0024 //
0025 
0026 // system include files
0027 #include <memory>
0028 #include <fstream>
0029 
0030 // user include files
0031 #include "FWCore/Framework/interface/Frameworkfwd.h"
0032 #include "FWCore/Framework/interface/global/EDFilter.h"
0033 
0034 #include "FWCore/Framework/interface/Event.h"
0035 #include "FWCore/Framework/interface/MakerMacros.h"
0036 
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 
0039 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0040 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0041 #include "DataFormats/DetId/interface/DetId.h"
0042 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0043 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0044 
0045 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0046 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0047 
0048 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0049 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0050 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0051 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0052 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0053 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0054 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0055 
0056 #include "RecoMET/METFilters/interface/EcalBoundaryInfoCalculator.h"
0057 #include "DataFormats/METReco/interface/AnomalousECALVariables.h"
0058 
0059 class EcalDeadCellBoundaryEnergyFilter : public edm::global::EDFilter<> {
0060 public:
0061   explicit EcalDeadCellBoundaryEnergyFilter(const edm::ParameterSet&);
0062   ~EcalDeadCellBoundaryEnergyFilter() override;
0063 
0064 private:
0065   void beginJob() override;
0066   bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0067   void endJob() override;
0068 
0069   // ----------member data ---------------------------
0070   const int kMAX;
0071 
0072   edm::EDGetTokenT<EcalRecHitCollection> EBRecHitsToken_;
0073   edm::EDGetTokenT<EcalRecHitCollection> EERecHitsToken_;
0074   edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyToken_;
0075   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalStatusToken_;
0076   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0077 
0078   const std::string FilterAlgo_;
0079   const bool taggingMode_;
0080 
0081   const bool skimGap_;
0082   const bool skimDead_;
0083 
0084   const double cutBoundEnergyGapEE, cutBoundEnergyGapEB, cutBoundEnergyDeadCellsEB, cutBoundEnergyDeadCellsEE;
0085 
0086   EcalBoundaryInfoCalculator<EBDetId> ebBoundaryCalc;
0087   EcalBoundaryInfoCalculator<EEDetId> eeBoundaryCalc;
0088 
0089   double maxBoundaryEnergy_;
0090 
0091   const bool limitFilterToEB_, limitFilterToEE_;
0092   const std::vector<int> limitDeadCellToChannelStatusEB_;
0093   const std::vector<int> limitDeadCellToChannelStatusEE_;
0094 
0095   const bool enableGap_;
0096   const bool debug_;
0097 };
0098 
0099 //
0100 // static data member definitions
0101 //
0102 
0103 //
0104 // constructors and destructor
0105 //
0106 EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter(const edm::ParameterSet& iConfig)
0107     : kMAX(50)
0108       //now do what ever initialization is needed
0109       ,
0110       EBRecHitsToken_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsEB"))),
0111       EERecHitsToken_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitsEE"))),
0112       caloTopologyToken_(esConsumes()),
0113       ecalStatusToken_(esConsumes()),
0114       geometryToken_(esConsumes()),
0115       FilterAlgo_(iConfig.getUntrackedParameter<std::string>("FilterAlgo", "TuningMode")),
0116       taggingMode_(iConfig.getParameter<bool>("taggingMode")),
0117       skimGap_(iConfig.getUntrackedParameter<bool>("skimGap", false)),
0118       skimDead_(iConfig.getUntrackedParameter<bool>("skimDead", false)),
0119       cutBoundEnergyGapEE(iConfig.getUntrackedParameter<double>("cutBoundEnergyGapEE")),
0120       cutBoundEnergyGapEB(iConfig.getUntrackedParameter<double>("cutBoundEnergyGapEB")),
0121       cutBoundEnergyDeadCellsEB(iConfig.getUntrackedParameter<double>("cutBoundEnergyDeadCellsEB")),
0122       cutBoundEnergyDeadCellsEE(iConfig.getUntrackedParameter<double>("cutBoundEnergyDeadCellsEE"))
0123 
0124       ,
0125       limitFilterToEB_(iConfig.getUntrackedParameter<bool>("limitFilterToEB", false)),
0126       limitFilterToEE_(iConfig.getUntrackedParameter<bool>("limitFilterToEE", false)),
0127       limitDeadCellToChannelStatusEB_(iConfig.getParameter<std::vector<int> >("limitDeadCellToChannelStatusEB")),
0128       limitDeadCellToChannelStatusEE_(iConfig.getParameter<std::vector<int> >("limitDeadCellToChannelStatusEE"))
0129 
0130       ,
0131       enableGap_(iConfig.getUntrackedParameter<bool>("enableGap", false)),
0132       debug_(iConfig.getParameter<bool>("debug")) {
0133   maxBoundaryEnergy_ = cutBoundEnergyDeadCellsEB;
0134 
0135   if (skimGap_ && debug_)
0136     edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Skim Gap!";
0137   if (skimDead_ && debug_)
0138     edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Skim Dead!";
0139 
0140   if (debug_) {
0141     edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Constructor EcalAnomalousEvent";
0142     ebBoundaryCalc.setDebugMode();
0143     eeBoundaryCalc.setDebugMode();
0144   }
0145 
0146   produces<bool>();
0147 
0148   produces<AnomalousECALVariables>("anomalousECALVariables");
0149 }
0150 
0151 EcalDeadCellBoundaryEnergyFilter::~EcalDeadCellBoundaryEnergyFilter() {}
0152 
0153 // ------------ method called on each new Event  ------------
0154 bool EcalDeadCellBoundaryEnergyFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0155   using namespace edm;
0156 
0157   //int eventno = (int) iEvent.eventAuxiliary().event();
0158 
0159   std::vector<BoundaryInformation> v_enNeighboursGap_EB;
0160   std::vector<BoundaryInformation> v_enNeighboursGap_EE;
0161   v_enNeighboursGap_EB.reserve(50);
0162   v_enNeighboursGap_EE.reserve(50);
0163 
0164   std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EB;
0165   std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EE;
0166   v_boundaryInfoDeadCells_EB.reserve(50);
0167   v_boundaryInfoDeadCells_EE.reserve(50);
0168 
0169   // Get the Ecal RecHits
0170   auto const& EBRecHits = iEvent.get(EBRecHitsToken_);
0171   auto const& EERecHits = iEvent.get(EERecHitsToken_);
0172 
0173   auto const& theCaloTopology = iSetup.getData(caloTopologyToken_);
0174   auto const& ecalStatus = iSetup.getData(ecalStatusToken_);
0175   auto const& geometry = iSetup.getData(geometryToken_);
0176 
0177   //   int DeadChannelsCounterEB = 0;
0178   //   int DeadChannelsCounterEE = 0;
0179 
0180   int i_EBDead = 0;
0181   int i_EEDead = 0;
0182   int i_EBGap = 0;
0183   int i_EEGap = 0;
0184 
0185   std::vector<DetId> sameFlagDetIds;
0186 
0187   bool pass = true;
0188 
0189   if (!limitFilterToEE_) {
0190     if (debug_)
0191       edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "process EB";
0192 
0193     for (EcalRecHitCollection::const_iterator hit = EBRecHits.begin(); hit != EBRecHits.end(); ++hit) {
0194       bool detIdAlreadyChecked = false;
0195       DetId currDetId = (DetId)hit->id();
0196       //add limitation to channel stati
0197       EcalChannelStatus::const_iterator chit = ecalStatus.find(currDetId);
0198       int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0199       if (status != 0)
0200         continue;
0201       bool passChannelLimitation = false;
0202 
0203       // check if hit has a dead neighbour
0204       std::vector<int> deadNeighbourStati;
0205       ebBoundaryCalc.checkRecHitHasDeadNeighbour(*hit, ecalStatus, deadNeighbourStati);
0206 
0207       for (int cs = 0; cs < (int)limitDeadCellToChannelStatusEB_.size(); ++cs) {
0208         int channelAllowed = limitDeadCellToChannelStatusEB_[cs];
0209 
0210         for (std::vector<int>::iterator sit = deadNeighbourStati.begin(); sit != deadNeighbourStati.end(); ++sit) {
0211           if (channelAllowed == *sit || (channelAllowed < 0 && std::abs(channelAllowed) <= *sit)) {
0212             passChannelLimitation = true;
0213             break;
0214           }
0215         }
0216       }
0217 
0218       for (std::vector<DetId>::iterator it = sameFlagDetIds.begin(); it != sameFlagDetIds.end(); it++) {
0219         if (currDetId == *it)
0220           detIdAlreadyChecked = true;
0221       }
0222 
0223       // RecHit is at EB boundary and should be processed
0224       if (!detIdAlreadyChecked && deadNeighbourStati.empty() &&
0225           ebBoundaryCalc.checkRecHitHasInvalidNeighbour(*hit, ecalStatus)) {
0226         BoundaryInformation gapinfo =
0227             ebBoundaryCalc.gapRecHits(EBRecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0228 
0229         // get rechits along gap cluster
0230         for (std::vector<DetId>::iterator it = gapinfo.detIds.begin(); it != gapinfo.detIds.end(); it++) {
0231           sameFlagDetIds.push_back(*it);
0232         }
0233 
0234         if (gapinfo.boundaryEnergy > cutBoundEnergyGapEB) {
0235           i_EBGap++;
0236           v_enNeighboursGap_EB.push_back(gapinfo);
0237 
0238           if (debug_)
0239             edm::LogInfo("EcalDeadCellBoundaryEnergyFilter")
0240                 << "EB: gap cluster energy: " << gapinfo.boundaryEnergy << " deadCells: " << gapinfo.detIds.size();
0241         }
0242       }
0243 
0244       // RecHit is member of boundary and should be processed
0245       if (!detIdAlreadyChecked &&
0246           (passChannelLimitation || (limitDeadCellToChannelStatusEB_.empty() && !deadNeighbourStati.empty()))) {
0247         BoundaryInformation boundinfo =
0248             ebBoundaryCalc.boundaryRecHits(EBRecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0249 
0250         // get boundary of !kDead rechits arround the dead cluster
0251         for (std::vector<DetId>::iterator it = boundinfo.detIds.begin(); it != boundinfo.detIds.end(); it++) {
0252           sameFlagDetIds.push_back(*it);
0253         }
0254 
0255         if (boundinfo.boundaryEnergy > cutBoundEnergyDeadCellsEB) {
0256           i_EBDead++;
0257           v_boundaryInfoDeadCells_EB.push_back(boundinfo);
0258 
0259           if (debug_)
0260             edm::LogInfo("EcalDeadCellBoundaryEnergyFilter")
0261                 << "EB: boundary Energy dead RecHit: " << boundinfo.boundaryEnergy << " ET: " << boundinfo.boundaryET
0262                 << " deadCells: " << boundinfo.detIds.size();
0263         }
0264       }
0265     }
0266 
0267   }  //// End of EB
0268 
0269   sameFlagDetIds.clear();
0270 
0271   if (!limitFilterToEB_) {
0272     if (debug_)
0273       edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "process EE";
0274 
0275     for (EcalRecHitCollection::const_iterator hit = EERecHits.begin(); hit != EERecHits.end(); ++hit) {
0276       bool detIdAlreadyChecked = false;
0277       DetId currDetId = (DetId)hit->id();
0278       //add limitation to channel stati
0279       EcalChannelStatus::const_iterator chit = ecalStatus.find(currDetId);
0280       int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0281       if (status != 0)
0282         continue;
0283       bool passChannelLimitation = false;
0284 
0285       // check if hit has a dead neighbour
0286       std::vector<int> deadNeighbourStati;
0287       eeBoundaryCalc.checkRecHitHasDeadNeighbour(*hit, ecalStatus, deadNeighbourStati);
0288 
0289       for (int cs = 0; cs < (int)limitDeadCellToChannelStatusEE_.size(); ++cs) {
0290         int channelAllowed = limitDeadCellToChannelStatusEE_[cs];
0291 
0292         for (std::vector<int>::iterator sit = deadNeighbourStati.begin(); sit != deadNeighbourStati.end(); ++sit) {
0293           if (channelAllowed == *sit || (channelAllowed < 0 && std::abs(channelAllowed) <= *sit)) {
0294             passChannelLimitation = true;
0295             break;
0296           }
0297         }
0298       }
0299 
0300       for (std::vector<DetId>::iterator it = sameFlagDetIds.begin(); it != sameFlagDetIds.end(); it++) {
0301         if (currDetId == *it)
0302           detIdAlreadyChecked = true;
0303       }
0304 
0305       // RecHit is at EE boundary and should be processed
0306       const CaloSubdetectorGeometry* subGeom = geometry.getSubdetectorGeometry(currDetId);
0307       auto cellGeom = subGeom->getGeometry(currDetId);
0308       double eta = cellGeom->getPosition().eta();
0309 
0310       if (!detIdAlreadyChecked && deadNeighbourStati.empty() &&
0311           eeBoundaryCalc.checkRecHitHasInvalidNeighbour(*hit, ecalStatus) && std::abs(eta) < 1.6) {
0312         BoundaryInformation gapinfo =
0313             eeBoundaryCalc.gapRecHits(EERecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0314 
0315         // get rechits along gap cluster
0316         for (std::vector<DetId>::iterator it = gapinfo.detIds.begin(); it != gapinfo.detIds.end(); it++) {
0317           sameFlagDetIds.push_back(*it);
0318         }
0319 
0320         if (gapinfo.boundaryEnergy > cutBoundEnergyGapEE) {
0321           i_EEGap++;
0322           v_enNeighboursGap_EE.push_back(gapinfo);
0323 
0324           if (debug_)
0325             edm::LogInfo("EcalDeadCellBoundaryEnergyFilter")
0326                 << "EE: gap cluster energy: " << gapinfo.boundaryEnergy << " deadCells: " << gapinfo.detIds.size();
0327         }
0328       }
0329 
0330       // RecHit is member of boundary and should be processed
0331       if (!detIdAlreadyChecked &&
0332           (passChannelLimitation || (limitDeadCellToChannelStatusEE_.empty() && !deadNeighbourStati.empty()))) {
0333         BoundaryInformation boundinfo =
0334             eeBoundaryCalc.boundaryRecHits(EERecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0335 
0336         // get boundary of !kDead rechits arround the dead cluster
0337         for (std::vector<DetId>::iterator it = boundinfo.detIds.begin(); it != boundinfo.detIds.end(); it++) {
0338           sameFlagDetIds.push_back(*it);
0339         }
0340 
0341         if (boundinfo.boundaryEnergy > cutBoundEnergyDeadCellsEE) {
0342           i_EEDead++;
0343           v_boundaryInfoDeadCells_EE.push_back(boundinfo);
0344 
0345           if (debug_)
0346             edm::LogInfo("EcalDeadCellBoundaryEnergyFilter")
0347                 << "EE: boundary Energy dead RecHit: " << boundinfo.boundaryEnergy << " ET: " << boundinfo.boundaryET
0348                 << " deadCells: " << boundinfo.detIds.size();
0349         }
0350       }
0351     }
0352 
0353   }  //// End of EE
0354 
0355   sameFlagDetIds.clear();
0356 
0357   auto pAnomalousECALVariables = std::make_unique<AnomalousECALVariables>(
0358       v_enNeighboursGap_EB, v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE);
0359 
0360   bool isGap = pAnomalousECALVariables->isGapEcalCluster(cutBoundEnergyGapEB, cutBoundEnergyGapEE);
0361   bool isBoundary = pAnomalousECALVariables->isDeadEcalCluster(
0362       maxBoundaryEnergy_, limitDeadCellToChannelStatusEB_, limitDeadCellToChannelStatusEE_);
0363   pass = (!isBoundary && ((!isGap && enableGap_) || !enableGap_));
0364 
0365   iEvent.put(std::move(pAnomalousECALVariables), "anomalousECALVariables");
0366 
0367   iEvent.put(std::make_unique<bool>(pass));
0368 
0369   if (taggingMode_) {
0370     if (skimDead_ && (i_EBDead >= 1 || i_EEDead >= 1)) {
0371       return true;
0372     } else if (skimGap_ && (i_EBGap >= 1 || i_EEGap >= 1)) {
0373       return true;
0374     } else if (!skimDead_ && !skimGap_)
0375       return true;
0376     else {
0377       return false;
0378     }
0379   } else
0380     return pass;
0381 
0382   /*
0383    if (FilterAlgo_ == "TuningMode") {
0384       std::unique_ptr<AnomalousECALVariables> pAnomalousECALVariables(new AnomalousECALVariables(v_enNeighboursGap_EB,
0385             v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE));
0386       iEvent.put(std::move(pAnomalousECALVariables), "anomalousECALVariables");
0387 
0388       if (skimDead_ && (i_EBDead >= 1 || i_EEDead >= 1)) {
0389          return true;
0390       } else if (skimGap_ && (i_EBGap >= 1 || i_EEGap >= 1)) {
0391          return true;
0392       } else if (!skimDead_ && !skimGap_)
0393          return true;
0394       else {
0395          return false;
0396       }
0397    }
0398 
0399    if (FilterAlgo_ == "FilterMode") {
0400       std::unique_ptr<AnomalousECALVariables> pAnomalousECALVariables(new AnomalousECALVariables(v_enNeighboursGap_EB,
0401             v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE));
0402 
0403       bool isGap = pAnomalousECALVariables->isGapEcalCluster(cutBoundEnergyGapEB, cutBoundEnergyGapEE);
0404       bool isBoundary = pAnomalousECALVariables->isDeadEcalCluster(maxBoundaryEnergy_, limitDeadCellToChannelStatusEB_,
0405             limitDeadCellToChannelStatusEE_);
0406 
0407       bool result = (!isBoundary && ((!isGap && enableGap_) || !enableGap_));
0408       if (!result) {
0409       }
0410       return result;
0411    }
0412 */
0413 
0414   //   return true;
0415 }
0416 
0417 // ------------ method called once each job just before starting event loop  ------------
0418 void EcalDeadCellBoundaryEnergyFilter::beginJob() {}
0419 
0420 // ------------ method called once each job just after ending the event loop  ------------
0421 void EcalDeadCellBoundaryEnergyFilter::endJob() {}
0422 
0423 //define this as a plug-in
0424 DEFINE_FWK_MODULE(EcalDeadCellBoundaryEnergyFilter);