File indexing completed on 2024-04-06 12:26:44
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027 #include <memory>
0028 #include <fstream>
0029
0030
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
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
0101
0102
0103
0104
0105
0106 EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter(const edm::ParameterSet& iConfig)
0107 : kMAX(50)
0108
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
0154 bool EcalDeadCellBoundaryEnergyFilter::filter(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0155 using namespace edm;
0156
0157
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
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
0178
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
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
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
0224 if (!detIdAlreadyChecked && deadNeighbourStati.empty() &&
0225 ebBoundaryCalc.checkRecHitHasInvalidNeighbour(*hit, ecalStatus)) {
0226 BoundaryInformation gapinfo =
0227 ebBoundaryCalc.gapRecHits(EBRecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0228
0229
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
0245 if (!detIdAlreadyChecked &&
0246 (passChannelLimitation || (limitDeadCellToChannelStatusEB_.empty() && !deadNeighbourStati.empty()))) {
0247 BoundaryInformation boundinfo =
0248 ebBoundaryCalc.boundaryRecHits(EBRecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0249
0250
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 }
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
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
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
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
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
0331 if (!detIdAlreadyChecked &&
0332 (passChannelLimitation || (limitDeadCellToChannelStatusEE_.empty() && !deadNeighbourStati.empty()))) {
0333 BoundaryInformation boundinfo =
0334 eeBoundaryCalc.boundaryRecHits(EERecHits, &(*hit), theCaloTopology, ecalStatus, geometry);
0335
0336
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 }
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
0384
0385
0386
0387
0388
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405
0406
0407
0408
0409
0410
0411
0412
0413
0414
0415 }
0416
0417
0418 void EcalDeadCellBoundaryEnergyFilter::beginJob() {}
0419
0420
0421 void EcalDeadCellBoundaryEnergyFilter::endJob() {}
0422
0423
0424 DEFINE_FWK_MODULE(EcalDeadCellBoundaryEnergyFilter);