File indexing completed on 2024-09-07 04:37:45
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016 #include <memory>
0017
0018
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/ESHandle.h"
0021
0022
0023 #include "PhysicsTools/SelectorUtils/interface/JetIDSelectionFunctor.h"
0024 #include "PhysicsTools/SelectorUtils/interface/strbitset.h"
0025
0026 #include "FWCore/Framework/interface/Frameworkfwd.h"
0027 #include "FWCore/Framework/interface/stream/EDFilter.h"
0028
0029 #include "FWCore/Framework/interface/Event.h"
0030 #include "FWCore/Framework/interface/MakerMacros.h"
0031
0032 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0033
0034 #include "FWCore/Utilities/interface/Exception.h"
0035
0036 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0037 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0038 #include "DataFormats/DetId/interface/DetId.h"
0039 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0040 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0041
0042 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0043 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0044
0045 #include "DataFormats/METReco/interface/METCollection.h"
0046 #include "DataFormats/METReco/interface/CaloMET.h"
0047 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0048
0049 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0050 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0051 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0052
0053 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0054 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0055 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0056
0057
0058 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0059 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0060 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0061 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0062
0063 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0064 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
0065 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0066
0067 #include "DataFormats/METReco/interface/PFMETCollection.h"
0068 #include "DataFormats/METReco/interface/PFMET.h"
0069
0070 #include "RecoJets/JetProducers/interface/JetIDHelper.h"
0071
0072 #include "DataFormats/JetReco/interface/GenJet.h"
0073 #include "DataFormats/JetReco/interface/GenJetCollection.h"
0074 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0075 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0076 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0077
0078 #include "FWCore/Common/interface/TriggerNames.h"
0079 #include "FWCore/Framework/interface/TriggerNamesService.h"
0080 #include "DataFormats/Common/interface/TriggerResults.h"
0081
0082 #include "DataFormats/Math/interface/deltaR.h"
0083
0084 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0085
0086 class EcalDeadCellDeltaRFilter : public edm::stream::EDFilter<> {
0087 public:
0088 explicit EcalDeadCellDeltaRFilter(const edm::ParameterSet &);
0089 ~EcalDeadCellDeltaRFilter() override = default;
0090
0091 private:
0092 bool filter(edm::Event &, const edm::EventSetup &) override;
0093 void beginRun(const edm::Run &, const edm::EventSetup &) override;
0094 void endRun(const edm::Run &, const edm::EventSetup &) override;
0095 virtual void envSet(const edm::EventSetup &);
0096
0097
0098 edm::EDGetTokenT<edm::View<reco::Jet> > jetToken_;
0099 edm::Handle<edm::View<reco::Jet> > jets;
0100
0101
0102 const std::vector<double> jetSelCuts_;
0103
0104 edm::EDGetTokenT<edm::View<reco::MET> > metToken_;
0105 edm::Handle<edm::View<reco::MET> > met;
0106
0107 const bool debug_, printSkimInfo_;
0108
0109 bool isPrintedOnce;
0110
0111 void loadEventInfo(const edm::Event &iEvent, const edm::EventSetup &iSetup);
0112 void loadJets(const edm::Event &iEvent, const edm::EventSetup &iSetup);
0113 void loadMET(const edm::Event &iEvent, const edm::EventSetup &iSetup);
0114
0115 edm::RunNumber_t run;
0116 edm::EventNumber_t event;
0117 edm::LuminosityBlockNumber_t ls;
0118 bool isdata;
0119
0120 double calomet, calometPhi, tcmet, tcmetPhi, pfmet, pfmetPhi;
0121
0122
0123 edm::ESHandle<EcalChannelStatus> ecalStatus;
0124 edm::ESHandle<CaloGeometry> geometry;
0125 const EcalTrigTowerConstituentsMap *ttMap_;
0126 const edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalChannelStatusToken_;
0127 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeometryToken_;
0128 const edm::ESGetToken<EcalTrigTowerConstituentsMap, IdealGeometryRecord> ecalTrigTowerConstituentsMapToken_;
0129
0130 const int maskedEcalChannelStatusThreshold_;
0131 const int chnStatusToBeEvaluated_;
0132
0133
0134
0135 std::map<DetId, std::vector<double> > EcalAllDeadChannelsValMap;
0136
0137
0138 std::map<DetId, std::vector<int> > EcalAllDeadChannelsBitMap;
0139
0140
0141 std::map<DetId, EcalTrigTowerDetId> EcalAllDeadChannelsTTMap;
0142
0143 int getChannelStatusMaps();
0144
0145 int evtProcessedCnt, totTPFilteredCnt;
0146 double wtdEvtProcessed, wtdTPFiltered;
0147
0148 const bool isProd_;
0149 const int verbose_;
0150
0151 const bool doCracks_;
0152
0153 const std::vector<double> cracksHBHEdef_, cracksHEHFdef_;
0154
0155
0156 const std::vector<double> EcalDeadCellDeltaRFilterInput_;
0157
0158 int dPhiToMETfunc(const std::vector<reco::Jet> &jetTVec,
0159 const double &dPhiCutVal,
0160 std::vector<reco::Jet> &closeToMETjetsVec);
0161 int dRtoMaskedChnsEvtFilterFunc(const std::vector<reco::Jet> &jetTVec, const int &chnStatus, const double &dRCutVal);
0162
0163 int etaToBoundary(const std::vector<reco::Jet> &jetTVec);
0164
0165 int isCloseToBadEcalChannel(const reco::Jet &jet,
0166 const double &deltaRCut,
0167 const int &chnStatus,
0168 std::map<double, DetId> &deltaRdetIdMap);
0169
0170 const bool taggingMode_;
0171 };
0172
0173 void EcalDeadCellDeltaRFilter::loadMET(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0174 iEvent.getByToken(metToken_, met);
0175 }
0176
0177 void EcalDeadCellDeltaRFilter::loadEventInfo(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0178 run = iEvent.id().run();
0179 event = iEvent.id().event();
0180 ls = iEvent.luminosityBlock();
0181 isdata = iEvent.isRealData();
0182
0183 if (!isPrintedOnce) {
0184 if (debug_) {
0185 if (isdata)
0186 std::cout << "\nInput dataset is DATA" << std::endl << std::endl;
0187 else
0188 std::cout << "\nInput dataset is MC" << std::endl << std::endl;
0189 }
0190 isPrintedOnce = true;
0191 }
0192 }
0193
0194 void EcalDeadCellDeltaRFilter::loadJets(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0195 iEvent.getByToken(jetToken_, jets);
0196 }
0197
0198
0199
0200
0201
0202
0203
0204
0205 EcalDeadCellDeltaRFilter::EcalDeadCellDeltaRFilter(const edm::ParameterSet &iConfig)
0206 : jetToken_(consumes<edm::View<reco::Jet> >(iConfig.getParameter<edm::InputTag>("jetInputTag"))),
0207 jetSelCuts_(iConfig.getParameter<std::vector<double> >("jetSelCuts")),
0208 metToken_(consumes<edm::View<reco::MET> >(iConfig.getParameter<edm::InputTag>("metInputTag"))),
0209 debug_(iConfig.getUntrackedParameter<bool>("debug", false)),
0210 printSkimInfo_(iConfig.getUntrackedParameter<bool>("printSkimInfo", false)),
0211 ecalChannelStatusToken_(esConsumes<edm::Transition::BeginRun>()),
0212 caloGeometryToken_(esConsumes<edm::Transition::BeginRun>()),
0213 ecalTrigTowerConstituentsMapToken_(esConsumes<edm::Transition::BeginRun>()),
0214 maskedEcalChannelStatusThreshold_(iConfig.getParameter<int>("maskedEcalChannelStatusThreshold")),
0215 chnStatusToBeEvaluated_(iConfig.getParameter<int>("chnStatusToBeEvaluated")),
0216 isProd_(iConfig.getUntrackedParameter<bool>("isProd")),
0217 verbose_(iConfig.getParameter<int>("verbose")),
0218 doCracks_(iConfig.getUntrackedParameter<bool>("doCracks")),
0219 cracksHBHEdef_(iConfig.getParameter<std::vector<double> >("cracksHBHEdef")),
0220 cracksHEHFdef_(iConfig.getParameter<std::vector<double> >("cracksHEHFdef")),
0221 EcalDeadCellDeltaRFilterInput_(iConfig.getParameter<std::vector<double> >("EcalDeadCellDeltaRFilterInput")),
0222 taggingMode_(iConfig.getParameter<bool>("taggingMode")) {
0223 produces<int>("deadCellStatus");
0224 produces<int>("boundaryStatus");
0225 produces<bool>();
0226 }
0227
0228 void EcalDeadCellDeltaRFilter::envSet(const edm::EventSetup &iSetup) {
0229 if (debug_)
0230 std::cout << "***envSet***" << std::endl;
0231
0232 ttMap_ = &iSetup.getData(ecalTrigTowerConstituentsMapToken_);
0233
0234 ecalStatus = iSetup.getHandle(ecalChannelStatusToken_);
0235 geometry = iSetup.getHandle(caloGeometryToken_);
0236
0237 if (!ecalStatus.isValid())
0238 throw cms::Exception("ESDataError") << "Failed to get ECAL channel status!";
0239 if (!geometry.isValid())
0240 throw cms::Exception("ESDataError") << "Failed to get the geometry!";
0241 }
0242
0243
0244 bool EcalDeadCellDeltaRFilter::filter(edm::Event &iEvent, const edm::EventSetup &iSetup) {
0245 loadEventInfo(iEvent, iSetup);
0246 loadJets(iEvent, iSetup);
0247 loadMET(iEvent, iSetup);
0248
0249
0250
0251 bool pass = true;
0252
0253 using namespace edm;
0254
0255 std::vector<reco::Jet> seledJets;
0256
0257 for (edm::View<reco::Jet>::const_iterator ij = jets->begin(); ij != jets->end(); ij++) {
0258 if (ij->pt() > jetSelCuts_[0] && std::abs(ij->eta()) < jetSelCuts_[1]) {
0259 seledJets.push_back(reco::Jet(*ij));
0260 }
0261 }
0262
0263 if (seledJets.empty())
0264 return pass;
0265
0266 double dPhiToMET = EcalDeadCellDeltaRFilterInput_[0], dRtoDeadCell = EcalDeadCellDeltaRFilterInput_[1];
0267
0268 std::vector<reco::Jet> closeToMETjetsVec;
0269
0270 int dPhiToMETstatus = dPhiToMETfunc(seledJets, dPhiToMET, closeToMETjetsVec);
0271
0272
0273 int deadCellStatus = dRtoMaskedChnsEvtFilterFunc(closeToMETjetsVec, chnStatusToBeEvaluated_, dRtoDeadCell);
0274
0275 int boundaryStatus = etaToBoundary(closeToMETjetsVec);
0276
0277 if (debug_) {
0278 printf("\nrun : %8u event : %12llu ls : %8u dPhiToMETstatus : %d deadCellStatus : %d boundaryStatus : %d\n",
0279 run,
0280 event,
0281 ls,
0282 dPhiToMETstatus,
0283 deadCellStatus,
0284 boundaryStatus);
0285 printf("met : %6.2f metphi : % 6.3f dPhiToMET : %5.3f dRtoDeadCell : %5.3f\n",
0286 (*met)[0].pt(),
0287 (*met)[0].phi(),
0288 dPhiToMET,
0289 dRtoDeadCell);
0290 }
0291
0292 iEvent.put(std::make_unique<int>(deadCellStatus), "deadCellStatus");
0293 iEvent.put(std::make_unique<int>(boundaryStatus), "boundaryStatus");
0294
0295 if (deadCellStatus || (doCracks_ && boundaryStatus))
0296 pass = false;
0297
0298 iEvent.put(std::make_unique<bool>(pass));
0299
0300 return taggingMode_ || pass;
0301 }
0302
0303
0304 void EcalDeadCellDeltaRFilter::beginRun(const edm::Run &run, const edm::EventSetup &iSetup) {
0305 if (debug_)
0306 std::cout << "beginRun" << std::endl;
0307
0308
0309 envSet(iSetup);
0310 getChannelStatusMaps();
0311 if (debug_)
0312 std::cout << "EcalAllDeadChannelsValMap.size() : " << EcalAllDeadChannelsValMap.size()
0313 << " EcalAllDeadChannelsBitMap.size() : " << EcalAllDeadChannelsBitMap.size() << std::endl;
0314 return;
0315 }
0316
0317
0318 void EcalDeadCellDeltaRFilter::endRun(const edm::Run &run, const edm::EventSetup &iSetup) {
0319 if (debug_)
0320 std::cout << "endRun" << std::endl;
0321 return;
0322 }
0323
0324 int EcalDeadCellDeltaRFilter::etaToBoundary(const std::vector<reco::Jet> &jetTVec) {
0325 int isClose = 0;
0326
0327 int cntOrder10 = 0;
0328 for (unsigned int ij = 0; ij < jetTVec.size(); ij++) {
0329 double recoJetEta = jetTVec[ij].eta();
0330
0331 if (std::abs(recoJetEta) > cracksHBHEdef_[0] && std::abs(recoJetEta) < cracksHBHEdef_[1])
0332 isClose += (cntOrder10 * 10 + 1);
0333 if (std::abs(recoJetEta) > cracksHEHFdef_[0] && std::abs(recoJetEta) < cracksHEHFdef_[1])
0334 isClose += (cntOrder10 * 10 + 2);
0335
0336 if (isClose / pow(10, cntOrder10) >= 3)
0337 cntOrder10 = isClose / 10 + 1;
0338 }
0339
0340 return isClose;
0341 }
0342
0343
0344 int EcalDeadCellDeltaRFilter::dPhiToMETfunc(const std::vector<reco::Jet> &jetTVec,
0345 const double &dPhiCutVal,
0346 std::vector<reco::Jet> &closeToMETjetsVec) {
0347 closeToMETjetsVec.clear();
0348
0349 double minDphi = 999.0;
0350 int minIdx = -1;
0351 for (unsigned int ii = 0; ii < jetTVec.size(); ii++) {
0352 const reco::Jet &jet = jetTVec[ii];
0353
0354 double deltaPhi = std::abs(reco::deltaPhi(jet.phi(), (*met)[0].phi()));
0355 if (deltaPhi > dPhiCutVal)
0356 continue;
0357
0358 closeToMETjetsVec.push_back(jetTVec[ii]);
0359
0360 if (deltaPhi < minDphi) {
0361 minDphi = deltaPhi;
0362 minIdx = ii;
0363 }
0364 }
0365
0366 if (minIdx == -1) {
0367 }
0368
0369
0370
0371 return (int)closeToMETjetsVec.size();
0372 }
0373
0374 int EcalDeadCellDeltaRFilter::dRtoMaskedChnsEvtFilterFunc(const std::vector<reco::Jet> &jetTVec,
0375 const int &chnStatus,
0376 const double &dRCutVal) {
0377 int isClose = 0;
0378
0379 for (unsigned int ii = 0; ii < jetTVec.size(); ii++) {
0380 const reco::Jet &jet = jetTVec[ii];
0381
0382 std::map<double, DetId> dummy;
0383 int isPerJetClose = isCloseToBadEcalChannel(jet, dRCutVal, chnStatus, dummy);
0384
0385 if (isPerJetClose) {
0386 isClose++;
0387 }
0388 }
0389
0390 return isClose;
0391 }
0392
0393 int EcalDeadCellDeltaRFilter::isCloseToBadEcalChannel(const reco::Jet &jet,
0394 const double &deltaRCut,
0395 const int &chnStatus,
0396 std::map<double, DetId> &deltaRdetIdMap) {
0397 double jetEta = jet.eta(), jetPhi = jet.phi();
0398
0399 deltaRdetIdMap.clear();
0400
0401 double min_dist = 999;
0402 DetId min_detId;
0403
0404 std::map<DetId, std::vector<int> >::iterator bitItor;
0405 for (bitItor = EcalAllDeadChannelsBitMap.begin(); bitItor != EcalAllDeadChannelsBitMap.end(); bitItor++) {
0406 DetId maskedDetId = bitItor->first;
0407
0408 int status = bitItor->second.back();
0409
0410 if (chnStatus > 0 && status != chnStatus)
0411 continue;
0412 if (chnStatus < 0 && status < abs(chnStatus))
0413 continue;
0414
0415 std::map<DetId, std::vector<double> >::iterator valItor = EcalAllDeadChannelsValMap.find(maskedDetId);
0416 if (valItor == EcalAllDeadChannelsValMap.end()) {
0417 std::cout << "Error cannot find maskedDetId in EcalAllDeadChannelsValMap ?!" << std::endl;
0418 continue;
0419 }
0420
0421 double eta = (valItor->second)[0], phi = (valItor->second)[1];
0422
0423 double dist = reco::deltaR(eta, phi, jetEta, jetPhi);
0424
0425 if (min_dist > dist) {
0426 min_dist = dist;
0427 min_detId = maskedDetId;
0428 }
0429 }
0430
0431 if (min_dist > deltaRCut && deltaRCut > 0)
0432 return 0;
0433
0434 deltaRdetIdMap.insert(std::make_pair(min_dist, min_detId));
0435
0436 return 1;
0437 }
0438
0439 int EcalDeadCellDeltaRFilter::getChannelStatusMaps() {
0440 EcalAllDeadChannelsValMap.clear();
0441 EcalAllDeadChannelsBitMap.clear();
0442
0443
0444 for (int ieta = -85; ieta <= 85; ieta++) {
0445 for (int iphi = 0; iphi <= 360; iphi++) {
0446 if (!EBDetId::validDetId(ieta, iphi))
0447 continue;
0448
0449 const EBDetId detid = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
0450 EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
0451
0452 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
0453
0454 const CaloSubdetectorGeometry *subGeom = geometry->getSubdetectorGeometry(detid);
0455 auto cellGeom = subGeom->getGeometry(detid);
0456 double eta = cellGeom->getPosition().eta();
0457 double phi = cellGeom->getPosition().phi();
0458 double theta = cellGeom->getPosition().theta();
0459
0460 if (status >= maskedEcalChannelStatusThreshold_) {
0461 std::vector<double> valVec;
0462 std::vector<int> bitVec;
0463 valVec.push_back(eta);
0464 valVec.push_back(phi);
0465 valVec.push_back(theta);
0466 bitVec.push_back(1);
0467 bitVec.push_back(ieta);
0468 bitVec.push_back(iphi);
0469 bitVec.push_back(status);
0470 EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
0471 EcalAllDeadChannelsBitMap.insert(std::make_pair(detid, bitVec));
0472 }
0473 }
0474 }
0475
0476
0477 for (int ix = 0; ix <= 100; ix++) {
0478 for (int iy = 0; iy <= 100; iy++) {
0479 for (int iz = -1; iz <= 1; iz++) {
0480 if (iz == 0)
0481 continue;
0482 if (!EEDetId::validDetId(ix, iy, iz))
0483 continue;
0484
0485 const EEDetId detid = EEDetId(ix, iy, iz, EEDetId::XYMODE);
0486 EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
0487 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
0488
0489 const CaloSubdetectorGeometry *subGeom = geometry->getSubdetectorGeometry(detid);
0490 auto cellGeom = subGeom->getGeometry(detid);
0491 double eta = cellGeom->getPosition().eta();
0492 double phi = cellGeom->getPosition().phi();
0493 double theta = cellGeom->getPosition().theta();
0494
0495 if (status >= maskedEcalChannelStatusThreshold_) {
0496 std::vector<double> valVec;
0497 std::vector<int> bitVec;
0498 valVec.push_back(eta);
0499 valVec.push_back(phi);
0500 valVec.push_back(theta);
0501 bitVec.push_back(2);
0502 bitVec.push_back(ix);
0503 bitVec.push_back(iy);
0504 bitVec.push_back(iz);
0505 bitVec.push_back(status);
0506 EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
0507 EcalAllDeadChannelsBitMap.insert(std::make_pair(detid, bitVec));
0508 }
0509 }
0510 }
0511 }
0512
0513 EcalAllDeadChannelsTTMap.clear();
0514 std::map<DetId, std::vector<int> >::iterator bitItor;
0515 for (bitItor = EcalAllDeadChannelsBitMap.begin(); bitItor != EcalAllDeadChannelsBitMap.end(); bitItor++) {
0516 const DetId id = bitItor->first;
0517 EcalTrigTowerDetId ttDetId = ttMap_->towerOf(id);
0518 EcalAllDeadChannelsTTMap.insert(std::make_pair(id, ttDetId));
0519 }
0520
0521 return 1;
0522 }
0523
0524
0525 DEFINE_FWK_MODULE(EcalDeadCellDeltaRFilter);