Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:45 

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalDeadCellTriggerPrimitiveFilter
0004 // Class:      EcalDeadCellTriggerPrimitiveFilter
0005 //
0006 /**\class EcalDeadCellTriggerPrimitiveFilter EcalDeadCellTriggerPrimitiveFilter.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 dead cell regions.
0011 */
0012 //
0013 // Original Author:  Hongxuan Liu and Kenichi Hatakeyama
0014 //                   in collaboration with Konstantinos Theofilatos and Ulla Gebbert
0015 
0016 // system include files
0017 #include <memory>
0018 #include <fstream>
0019 
0020 // user include files
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/ESHandle.h"
0023 
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/stream/EDFilter.h"
0026 
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029 
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032 
0033 #include "DataFormats/Provenance/interface/ProcessHistory.h"
0034 
0035 #include "DataFormats/EcalRecHit/interface/EcalRecHit.h"
0036 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0037 #include "DataFormats/DetId/interface/DetId.h"
0038 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0039 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0040 
0041 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0042 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0043 
0044 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0045 #include "CalibCalorimetry/EcalTPGTools/interface/EcalTPGScale.h"
0046 #include "Geometry/CaloTopology/interface/EcalTrigTowerConstituentsMap.h"
0047 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0048 
0049 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0050 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0051 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0052 
0053 #include "DataFormats/HepMCCandidate/interface/GenParticle.h"
0054 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0055 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0056 
0057 #include "DataFormats/EcalDetId/interface/EcalScDetId.h"
0058 
0059 // Geometry
0060 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0061 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0062 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0063 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0064 
0065 #include "Geometry/CaloTopology/interface/CaloTowerConstituentsMap.h"
0066 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0067 
0068 #include "DataFormats/Provenance/interface/RunLumiEventNumber.h"
0069 
0070 using namespace std;
0071 
0072 class EcalDeadCellTriggerPrimitiveFilter : public edm::stream::EDFilter<> {
0073 public:
0074   explicit EcalDeadCellTriggerPrimitiveFilter(const edm::ParameterSet&);
0075   ~EcalDeadCellTriggerPrimitiveFilter() override;
0076 
0077 private:
0078   void beginStream(edm::StreamID) override;
0079   bool filter(edm::Event&, const edm::EventSetup&) override;
0080   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0081   void envSet(const edm::EventSetup&);
0082 
0083   // ----------member data ---------------------------
0084 
0085   const bool taggingMode_;
0086 
0087   const bool debug_;
0088   const int verbose_;
0089 
0090   const bool doEEfilter_;
0091 
0092   // Channel status related
0093 
0094   edm::ESHandle<EcalChannelStatus> ecalStatus;
0095   edm::ESHandle<CaloGeometry> geometry;
0096   edm::ESHandle<EcalTrigTowerConstituentsMap> ttMap_;
0097   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalStatusToken_;
0098   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0099   edm::ESGetToken<EcalTrigTowerConstituentsMap, IdealGeometryRecord> ttmapToken_;
0100 
0101   void loadEcalDigis(edm::Event& iEvent, edm::Handle<EcalTrigPrimDigiCollection>& pTPDigis);
0102   void loadEcalRecHits(edm::Event& iEvent,
0103                        edm::Handle<EcalRecHitCollection>& barrelReducedRecHitsHandle,
0104                        edm::Handle<EcalRecHitCollection>& endcapReducedRecHitsHandle);
0105 
0106   const edm::InputTag ebReducedRecHitCollection_;
0107   const edm::InputTag eeReducedRecHitCollection_;
0108   edm::EDGetTokenT<EcalRecHitCollection> ebReducedRecHitCollectionToken_;
0109   edm::EDGetTokenT<EcalRecHitCollection> eeReducedRecHitCollectionToken_;
0110 
0111   const int maskedEcalChannelStatusThreshold_;
0112 
0113   // XXX: All the following can be built at the beginning of a job
0114   // Store DetId <==> vector<double> (eta, phi, theta)
0115   std::map<DetId, std::vector<double> > EcalAllDeadChannelsValMap;
0116   // Store EB: DetId <==> vector<int> (subdet, ieta, iphi, status)
0117   // Store EE: DetId <==> vector<int> (subdet, ix, iy, iz, status)
0118   std::map<DetId, std::vector<int> > EcalAllDeadChannelsBitMap;
0119 
0120   // Store DetId <==> EcalTrigTowerDetId
0121   std::map<DetId, EcalTrigTowerDetId> EcalAllDeadChannelsTTMap;
0122 
0123   int getChannelStatusMaps();
0124 
0125   // TP filter
0126   const double etValToBeFlagged_;
0127 
0128   const edm::InputTag tpDigiCollection_;
0129   edm::EDGetTokenT<EcalTrigPrimDigiCollection> tpDigiCollectionToken_;
0130 
0131   // chnStatus > 0, then exclusive, i.e., only consider status == chnStatus
0132   // chnStatus < 0, then inclusive, i.e., consider status >= abs(chnStatus)
0133   // Return value:  + : positive zside  - : negative zside
0134   int setEvtTPstatus(const EcalTrigPrimDigiCollection&, const double& tpCntCut, const int& chnStatus, EcalTPGScale&);
0135 
0136   const bool
0137       useTTsum_;  //If set to true, the filter will compare the sum of the 5x5 tower to the provided energy threshold
0138   const bool usekTPSaturated_;  //If set to true, the filter will check the kTPSaturated flag
0139 
0140   int hasReducedRecHits_ = 0;
0141 
0142   bool useTPmethod_ = false;
0143   bool useHITmethod_ = false;
0144 
0145   edm::EDPutTokenT<bool> putToken_;
0146   EcalTPGScale::Tokens tokens_;
0147 
0148   // Only for EB since the dead front-end has one-to-one map to TT
0149   std::map<EcalTrigTowerDetId, double> accuTTetMap;
0150   std::map<EcalTrigTowerDetId, int> accuTTchnMap;
0151   std::map<EcalTrigTowerDetId, int> TTzsideMap;
0152 
0153   // For EE, the one-to-one map to dead front-end is the SuperCrystal
0154   std::map<EcalScDetId, double> accuSCetMap;
0155   std::map<EcalScDetId, int> accuSCchnMap;
0156   std::map<EcalScDetId, int> SCzsideMap;
0157 
0158   // To be used before a bug fix
0159   std::vector<DetId> avoidDuplicateVec;
0160   int setEvtRecHitstatus(const double& tpValCut,
0161                          const int& chnStatus,
0162                          const int& towerTest,
0163                          const EBRecHitCollection& HitecalEB,
0164                          const EERecHitCollection& HitecalEE);
0165 };
0166 
0167 //
0168 // constructors and destructor
0169 //
0170 EcalDeadCellTriggerPrimitiveFilter::EcalDeadCellTriggerPrimitiveFilter(const edm::ParameterSet& iConfig)
0171     : taggingMode_(iConfig.getParameter<bool>("taggingMode")),
0172       debug_(iConfig.getParameter<bool>("debug")),
0173       verbose_(iConfig.getParameter<int>("verbose")),
0174       doEEfilter_(iConfig.getUntrackedParameter<bool>("doEEfilter")),
0175       ecalStatusToken_(esConsumes<EcalChannelStatus, EcalChannelStatusRcd, edm::Transition::BeginRun>()),
0176       geomToken_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0177       ttmapToken_(esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord, edm::Transition::BeginRun>()),
0178       ebReducedRecHitCollection_(iConfig.getParameter<edm::InputTag>("ebReducedRecHitCollection")),
0179       eeReducedRecHitCollection_(iConfig.getParameter<edm::InputTag>("eeReducedRecHitCollection")),
0180       maskedEcalChannelStatusThreshold_(iConfig.getParameter<int>("maskedEcalChannelStatusThreshold")),
0181       etValToBeFlagged_(iConfig.getParameter<double>("etValToBeFlagged")),
0182       tpDigiCollection_(iConfig.getParameter<edm::InputTag>("tpDigiCollection")),
0183       tpDigiCollectionToken_(consumes<EcalTrigPrimDigiCollection>(tpDigiCollection_)),
0184       useTTsum_(iConfig.getParameter<bool>("useTTsum")),
0185       usekTPSaturated_(iConfig.getParameter<bool>("usekTPSaturated")),
0186       putToken_(produces<bool>()),
0187       tokens_(consumesCollector()) {
0188   callWhenNewProductsRegistered([this](edm::BranchDescription const& iBranch) {
0189     // If TP is available, always use TP.
0190     // In RECO file, we always have ecalTPSkim (at least from 38X for data and 39X for MC).
0191     // In AOD file, we can only have recovered rechits in the reduced rechits collection after 42X
0192     // Do NOT expect end-users provide ecalTPSkim or recovered rechits themselves!!
0193     // If they really can provide them, they must be experts to modify this code to suit their own purpose :-)
0194     if (iBranch.moduleLabel() == tpDigiCollection_.label()) {
0195       useTPmethod_ = true;
0196       //if both collections are in the job then we may already have seen the reduced collections
0197       useHITmethod_ = false;
0198     }
0199     if (iBranch.moduleLabel() == ebReducedRecHitCollection_.label() ||
0200         iBranch.moduleLabel() == eeReducedRecHitCollection_.label()) {
0201       hasReducedRecHits_++;
0202       if (not useTPmethod_ && hasReducedRecHits_ == 2) {
0203         useHITmethod_ = true;
0204         ebReducedRecHitCollectionToken_ = consumes<EcalRecHitCollection>(ebReducedRecHitCollection_);
0205         eeReducedRecHitCollectionToken_ = consumes<EcalRecHitCollection>(eeReducedRecHitCollection_);
0206       }
0207     }
0208   });
0209 }
0210 
0211 EcalDeadCellTriggerPrimitiveFilter::~EcalDeadCellTriggerPrimitiveFilter() {}
0212 
0213 void EcalDeadCellTriggerPrimitiveFilter::beginStream(edm::StreamID) {
0214   if (debug_)
0215     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0216         << "\nuseTPmethod_ : " << useTPmethod_ << "  hasReducedRecHits_ : " << hasReducedRecHits_;
0217 
0218   if (not useTPmethod_ and not useHITmethod_) {
0219     if (debug_) {
0220       edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0221           << "\nWARNING ... Cannot find either tpDigiCollection_ or reducedRecHitCollecion_ ?!";
0222       edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "  Will NOT DO ANY FILTERING !";
0223     }
0224   }
0225 
0226   if (debug_)
0227     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0228         << "useTPmethod_ : " << useTPmethod_ << "  useHITmethod_ : " << useHITmethod_;
0229 }
0230 
0231 void EcalDeadCellTriggerPrimitiveFilter::loadEcalDigis(edm::Event& iEvent,
0232                                                        edm::Handle<EcalTrigPrimDigiCollection>& pTPDigis) {
0233   iEvent.getByToken(tpDigiCollectionToken_, pTPDigis);
0234   if (!pTPDigis.isValid()) {
0235     edm::LogWarning("EcalDeadCellTriggerPrimitiveFilter")
0236         << "Can't get the product " << tpDigiCollection_.instance() << " with label " << tpDigiCollection_.label();
0237     return;
0238   }
0239 }
0240 
0241 void EcalDeadCellTriggerPrimitiveFilter::loadEcalRecHits(edm::Event& iEvent,
0242                                                          edm::Handle<EcalRecHitCollection>& barrelReducedRecHitsHandle,
0243                                                          edm::Handle<EcalRecHitCollection>& endcapReducedRecHitsHandle) {
0244   iEvent.getByToken(ebReducedRecHitCollectionToken_, barrelReducedRecHitsHandle);
0245   iEvent.getByToken(eeReducedRecHitCollectionToken_, endcapReducedRecHitsHandle);
0246 }
0247 
0248 //
0249 // static data member definitions
0250 //
0251 
0252 void EcalDeadCellTriggerPrimitiveFilter::envSet(const edm::EventSetup& iSetup) {
0253   if (debug_ && verbose_ >= 2)
0254     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "***envSet***";
0255   ecalStatus = iSetup.getHandle(ecalStatusToken_);
0256   geometry = iSetup.getHandle(geomToken_);
0257   ttMap_ = iSetup.getHandle(ttmapToken_);
0258 
0259   if (!ecalStatus.isValid())
0260     throw "Failed to get ECAL channel status!";
0261   if (!geometry.isValid())
0262     throw "Failed to get the geometry!";
0263 }
0264 
0265 // ------------ method called on each new Event  ------------
0266 bool EcalDeadCellTriggerPrimitiveFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0267   edm::RunNumber_t run = iEvent.id().run();
0268   edm::EventNumber_t event = iEvent.id().event();
0269   edm::LuminosityBlockNumber_t ls = iEvent.luminosityBlock();
0270 
0271   bool pass = true;
0272 
0273   int evtTagged = 0;
0274 
0275   if (useTPmethod_) {
0276     edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
0277     loadEcalDigis(iEvent, pTPDigis);
0278 
0279     EcalTPGScale ecalScale(tokens_, iSetup);
0280     evtTagged = setEvtTPstatus(*pTPDigis, etValToBeFlagged_, 13, ecalScale);
0281   }
0282 
0283   if (useHITmethod_) {
0284     edm::Handle<EcalRecHitCollection> barrelReducedRecHitsHandle;
0285     edm::Handle<EcalRecHitCollection> endcapReducedRecHitsHandle;
0286     loadEcalRecHits(iEvent, barrelReducedRecHitsHandle, endcapReducedRecHitsHandle);
0287     evtTagged = setEvtRecHitstatus(etValToBeFlagged_, 13, 13, *barrelReducedRecHitsHandle, *endcapReducedRecHitsHandle);
0288   }
0289 
0290   if (evtTagged) {
0291     pass = false;
0292   }
0293 
0294   if (debug_ && verbose_ >= 2) {
0295     int evtstatusABS = abs(evtTagged);
0296     printf("\nrun : %8u  event : %10llu  lumi : %4u  evtTPstatus  ABS : %d  13 : % 2d\n",
0297            run,
0298            event,
0299            ls,
0300            evtstatusABS,
0301            evtTagged);
0302   }
0303 
0304   iEvent.emplace(putToken_, pass);
0305 
0306   if (taggingMode_)
0307     return true;
0308   else
0309     return pass;
0310 }
0311 
0312 // ------------ method called once each run just before starting event loop  ------------
0313 void EcalDeadCellTriggerPrimitiveFilter::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0314   // Channel status might change for each run (data)
0315   // Event setup
0316   envSet(iSetup);
0317   getChannelStatusMaps();
0318   if (debug_ && verbose_ >= 2)
0319     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0320         << "EcalAllDeadChannelsValMap.size() : " << EcalAllDeadChannelsValMap.size()
0321         << "  EcalAllDeadChannelsBitMap.size() : " << EcalAllDeadChannelsBitMap.size();
0322   return;
0323 }
0324 
0325 int EcalDeadCellTriggerPrimitiveFilter::setEvtRecHitstatus(const double& tpValCut,
0326                                                            const int& chnStatus,
0327                                                            const int& towerTest,
0328                                                            const EBRecHitCollection& HitecalEB,
0329                                                            const EERecHitCollection& HitecalEE) {
0330   if (debug_ && verbose_ >= 2)
0331     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "***begin setEvtTPstatusRecHits***";
0332 
0333   accuTTetMap.clear();
0334   accuTTchnMap.clear();
0335   TTzsideMap.clear();
0336   accuSCetMap.clear();
0337   accuSCchnMap.clear();
0338   SCzsideMap.clear();
0339   avoidDuplicateVec.clear();
0340 
0341   /*
0342   const EBRecHitCollection HitecalEB = *(barrelRecHitsHandle.product());
0343   const EERecHitCollection HitecalEE = *(endcapRecHitsHandle.product());
0344 */
0345   int isPassCut = 0;
0346 
0347   EBRecHitCollection::const_iterator ebrechit;
0348 
0349   for (ebrechit = HitecalEB.begin(); ebrechit != HitecalEB.end(); ebrechit++) {
0350     EBDetId det = ebrechit->id();
0351 
0352     std::map<DetId, vector<double> >::iterator valItor = EcalAllDeadChannelsValMap.find(det);
0353     if (valItor == EcalAllDeadChannelsValMap.end())
0354       continue;
0355 
0356     double theta = valItor->second.back();
0357 
0358     std::map<DetId, vector<int> >::iterator bitItor = EcalAllDeadChannelsBitMap.find(det);
0359     if (bitItor == EcalAllDeadChannelsBitMap.end())
0360       continue;
0361 
0362     std::map<DetId, EcalTrigTowerDetId>::iterator ttItor = EcalAllDeadChannelsTTMap.find(det);
0363     if (ttItor == EcalAllDeadChannelsTTMap.end())
0364       continue;
0365 
0366     int status = bitItor->second.back();
0367 
0368     bool toDo = false;
0369     if (chnStatus > 0 && status == chnStatus)
0370       toDo = true;
0371     if (chnStatus < 0 && status >= abs(chnStatus))
0372       toDo = true;
0373     // This might be suitable for channels with status other than 13,
0374     // since this function is written as a general one ...
0375     if (!ebrechit->isRecovered())
0376       toDo = false;
0377     //     if( !ebrechit->checkFlag(EcalRecHit::kTowerRecovered) ) toDo = false;
0378 
0379     if (toDo) {
0380       //If we considerkTPSaturated and a recHit has a flag set, we can immediately flag the event.
0381       if (ebrechit->checkFlag(EcalRecHit::kTPSaturated) && usekTPSaturated_)
0382         return 1;
0383 
0384       EcalTrigTowerDetId ttDetId = ttItor->second;
0385       int ttzside = ttDetId.zside();
0386 
0387       std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0388       int towerTestCnt = 0;
0389       for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0390         std::map<DetId, vector<int> >::iterator bit2Itor = EcalAllDeadChannelsBitMap.find((*dit));
0391         if (bit2Itor == EcalAllDeadChannelsBitMap.end()) {
0392           towerTestCnt++;
0393           continue;
0394         }
0395         if (towerTest > 0 && bit2Itor->second.back() == towerTest)
0396           continue;
0397         if (towerTest < 0 && bit2Itor->second.back() >= abs(towerTest))
0398           continue;
0399         towerTestCnt++;
0400       }
0401       if (towerTestCnt != 0 && debug_ && verbose_ >= 2)
0402         edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0403             << "towerTestCnt : " << towerTestCnt << "  for towerTest : " << towerTest;
0404 
0405       std::vector<DetId>::iterator avoidItor;
0406       avoidItor = find(avoidDuplicateVec.begin(), avoidDuplicateVec.end(), det);
0407       if (avoidItor == avoidDuplicateVec.end()) {
0408         avoidDuplicateVec.push_back(det);
0409       } else {
0410         continue;
0411       }
0412 
0413       std::map<EcalTrigTowerDetId, double>::iterator ttetItor = accuTTetMap.find(ttDetId);
0414       if (ttetItor == accuTTetMap.end()) {
0415         accuTTetMap[ttDetId] = ebrechit->energy() * sin(theta);
0416         accuTTchnMap[ttDetId] = 1;
0417         TTzsideMap[ttDetId] = ttzside;
0418       } else {
0419         accuTTetMap[ttDetId] += ebrechit->energy() * sin(theta);
0420         accuTTchnMap[ttDetId]++;
0421       }
0422     }
0423   }  // loop over EB
0424 
0425   EERecHitCollection::const_iterator eerechit;
0426   for (eerechit = HitecalEE.begin(); eerechit != HitecalEE.end(); eerechit++) {
0427     EEDetId det = eerechit->id();
0428 
0429     std::map<DetId, vector<double> >::iterator valItor = EcalAllDeadChannelsValMap.find(det);
0430     if (valItor == EcalAllDeadChannelsValMap.end())
0431       continue;
0432 
0433     double theta = valItor->second.back();
0434 
0435     std::map<DetId, vector<int> >::iterator bitItor = EcalAllDeadChannelsBitMap.find(det);
0436     if (bitItor == EcalAllDeadChannelsBitMap.end())
0437       continue;
0438 
0439     std::map<DetId, EcalTrigTowerDetId>::iterator ttItor = EcalAllDeadChannelsTTMap.find(det);
0440     if (ttItor == EcalAllDeadChannelsTTMap.end())
0441       continue;
0442 
0443     int status = bitItor->second.back();
0444 
0445     bool toDo = false;
0446     if (chnStatus > 0 && status == chnStatus)
0447       toDo = true;
0448     if (chnStatus < 0 && status >= abs(chnStatus))
0449       toDo = true;
0450     // This might be suitable for channels with status other than 13,
0451     // since this function is written as a general one ...
0452     if (!eerechit->isRecovered())
0453       toDo = false;
0454     //     if( !eerechit->checkFlag(EcalRecHit::kTowerRecovered) ) toDo = false;
0455 
0456     if (toDo) {
0457       //If we considerkTPSaturated and a recHit has a flag set, we can immediately flag the event.
0458       if (eerechit->checkFlag(EcalRecHit::kTPSaturated) && usekTPSaturated_)
0459         return 1;
0460 
0461       // vvvv= Only for debuging or testing purpose =vvvv
0462       EcalTrigTowerDetId ttDetId = ttItor->second;
0463       //        int ttzside = ttDetId.zside();
0464 
0465       std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0466       int towerTestCnt = 0;
0467       for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0468         std::map<DetId, vector<int> >::iterator bit2Itor = EcalAllDeadChannelsBitMap.find((*dit));
0469         if (bit2Itor == EcalAllDeadChannelsBitMap.end()) {
0470           towerTestCnt++;
0471           continue;
0472         }
0473         if (towerTest > 0 && bit2Itor->second.back() == towerTest)
0474           continue;
0475         if (towerTest < 0 && bit2Itor->second.back() >= abs(towerTest))
0476           continue;
0477         towerTestCnt++;
0478       }
0479       if (towerTestCnt != 0 && debug_ && verbose_ >= 2)
0480         edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0481             << "towerTestCnt : " << towerTestCnt << "  for towerTest : " << towerTest;
0482       // ^^^^=                  END                 =^^^^
0483 
0484       EcalScDetId sc((det.ix() - 1) / 5 + 1, (det.iy() - 1) / 5 + 1, det.zside());
0485 
0486       std::vector<DetId>::iterator avoidItor;
0487       avoidItor = find(avoidDuplicateVec.begin(), avoidDuplicateVec.end(), det);
0488       if (avoidItor == avoidDuplicateVec.end()) {
0489         avoidDuplicateVec.push_back(det);
0490       } else {
0491         continue;
0492       }
0493 
0494       std::map<EcalScDetId, double>::iterator scetItor = accuSCetMap.find(sc);
0495       if (scetItor == accuSCetMap.end()) {
0496         accuSCetMap[sc] = eerechit->energy() * sin(theta);
0497         accuSCchnMap[sc] = 1;
0498         SCzsideMap[sc] = sc.zside();
0499       } else {
0500         accuSCetMap[sc] += eerechit->energy() * sin(theta);
0501         accuSCchnMap[sc]++;
0502       }
0503     }
0504   }  // loop over EE
0505 
0506   //If we are not using the TT sum, then at this point we need not do anything further, we'll pass the event
0507   if (!useTTsum_)
0508     return 0;
0509 
0510   // Checking for EB
0511   std::map<EcalTrigTowerDetId, double>::iterator ttetItor;
0512   for (ttetItor = accuTTetMap.begin(); ttetItor != accuTTetMap.end(); ttetItor++) {
0513     EcalTrigTowerDetId ttDetId = ttetItor->first;
0514 
0515     double ttetVal = ttetItor->second;
0516 
0517     std::map<EcalTrigTowerDetId, int>::iterator ttchnItor = accuTTchnMap.find(ttDetId);
0518     if (ttchnItor == accuTTchnMap.end()) {
0519       edm::LogError("EcalDeadCellTriggerPrimitiveFilter")
0520           << "\nERROR  cannot find ttDetId : " << ttDetId << " in accuTTchnMap?!";
0521     }
0522 
0523     std::map<EcalTrigTowerDetId, int>::iterator ttzsideItor = TTzsideMap.find(ttDetId);
0524     if (ttzsideItor == TTzsideMap.end()) {
0525       edm::LogError("EcalDeadCellTriggerPrimitiveFilter")
0526           << "\nERROR  cannot find ttDetId : " << ttDetId << " in TTzsideMap?!";
0527     }
0528 
0529     if (ttchnItor->second != 25 && debug_ && verbose_ >= 2)
0530       edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0531           << "WARNING ... ttchnCnt : " << ttchnItor->second << "  NOT equal  25!";
0532 
0533     if (ttetVal >= tpValCut) {
0534       isPassCut = 1;
0535       isPassCut *= ttzsideItor->second;
0536     }
0537   }
0538 
0539   // Checking for EE
0540   std::map<EcalScDetId, double>::iterator scetItor;
0541   for (scetItor = accuSCetMap.begin(); scetItor != accuSCetMap.end(); scetItor++) {
0542     EcalScDetId scDetId = scetItor->first;
0543 
0544     double scetVal = scetItor->second;
0545 
0546     std::map<EcalScDetId, int>::iterator scchnItor = accuSCchnMap.find(scDetId);
0547     if (scchnItor == accuSCchnMap.end()) {
0548       edm::LogError("EcalDeadCellTriggerPrimitiveFilter")
0549           << "\nERROR  cannot find scDetId : " << scDetId << " in accuSCchnMap?!";
0550     }
0551 
0552     std::map<EcalScDetId, int>::iterator sczsideItor = SCzsideMap.find(scDetId);
0553     if (sczsideItor == SCzsideMap.end()) {
0554       edm::LogError("EcalDeadCellTriggerPrimitiveFilter")
0555           << "\nERROR  cannot find scDetId : " << scDetId << " in SCzsideMap?!";
0556     }
0557 
0558     if (scchnItor->second != 25 && debug_ && verbose_ >= 2)
0559       edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter")
0560           << "WARNING ... scchnCnt : " << scchnItor->second << "  NOT equal  25!";
0561 
0562     if (scetVal >= tpValCut) {
0563       isPassCut = 1;
0564       isPassCut *= sczsideItor->second;
0565     }
0566   }
0567 
0568   if (debug_ && verbose_ >= 2)
0569     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "***end setEvtTPstatusRecHits***";
0570 
0571   return isPassCut;
0572 }
0573 
0574 int EcalDeadCellTriggerPrimitiveFilter::setEvtTPstatus(EcalTrigPrimDigiCollection const& tpDigis,
0575                                                        const double& tpValCut,
0576                                                        const int& chnStatus,
0577                                                        EcalTPGScale& ecalScale) {
0578   if (debug_ && verbose_ >= 2)
0579     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "***begin setEvtTPstatus***";
0580 
0581   int isPassCut = 0;
0582 
0583   std::map<DetId, std::vector<int> >::iterator bitItor;
0584   for (bitItor = EcalAllDeadChannelsBitMap.begin(); bitItor != EcalAllDeadChannelsBitMap.end(); bitItor++) {
0585     DetId maskedDetId = bitItor->first;
0586     int subdet = bitItor->second.front(), status = bitItor->second.back();
0587 
0588     // if NOT filtering on EE, skip EE subdet
0589     if (!doEEfilter_ && subdet != 1)
0590       continue;
0591 
0592     std::map<DetId, EcalTrigTowerDetId>::iterator ttItor = EcalAllDeadChannelsTTMap.find(maskedDetId);
0593     if (ttItor == EcalAllDeadChannelsTTMap.end())
0594       continue;
0595 
0596     bool toDo = false;
0597     if (chnStatus > 0 && status == chnStatus)
0598       toDo = true;
0599     if (chnStatus < 0 && status >= abs(chnStatus))
0600       toDo = true;
0601 
0602     if (toDo) {
0603       EcalTrigTowerDetId ttDetId = ttItor->second;
0604       int ttzside = ttDetId.zside();
0605 
0606       EcalTrigPrimDigiCollection::const_iterator tp = tpDigis.find(ttDetId);
0607       if (tp != tpDigis.end()) {
0608         double tpEt = ecalScale.getTPGInGeV(tp->compressedEt(), tp->id());
0609         if (tpEt >= tpValCut) {
0610           isPassCut = 1;
0611           isPassCut *= ttzside;
0612         }
0613       }
0614     }
0615   }  // loop over EB + EE
0616 
0617   if (debug_ && verbose_ >= 2)
0618     edm::LogInfo("EcalDeadCellTriggerPrimitiveFilter") << "***end setEvtTPstatus***";
0619 
0620   return isPassCut;
0621 }
0622 
0623 int EcalDeadCellTriggerPrimitiveFilter::getChannelStatusMaps() {
0624   EcalAllDeadChannelsValMap.clear();
0625   EcalAllDeadChannelsBitMap.clear();
0626 
0627   // Loop over EB ...
0628   for (int ieta = -85; ieta <= 85; ieta++) {
0629     for (int iphi = 0; iphi <= 360; iphi++) {
0630       if (!EBDetId::validDetId(ieta, iphi))
0631         continue;
0632 
0633       const EBDetId detid = EBDetId(ieta, iphi, EBDetId::ETAPHIMODE);
0634       EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
0635       // refer https://twiki.cern.ch/twiki/bin/viewauth/CMS/EcalChannelStatus
0636       int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
0637 
0638       const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
0639       auto cellGeom = subGeom->getGeometry(detid);
0640       double eta = cellGeom->getPosition().eta();
0641       double phi = cellGeom->getPosition().phi();
0642       double theta = cellGeom->getPosition().theta();
0643 
0644       if (status >= maskedEcalChannelStatusThreshold_) {
0645         std::vector<double> valVec;
0646         std::vector<int> bitVec;
0647         valVec.push_back(eta);
0648         valVec.push_back(phi);
0649         valVec.push_back(theta);
0650         bitVec.push_back(1);
0651         bitVec.push_back(ieta);
0652         bitVec.push_back(iphi);
0653         bitVec.push_back(status);
0654         EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
0655         EcalAllDeadChannelsBitMap.insert(std::make_pair(detid, bitVec));
0656       }
0657     }  // end loop iphi
0658   }  // end loop ieta
0659 
0660   // Loop over EE detid
0661   if (doEEfilter_) {
0662     for (int ix = 0; ix <= 100; ix++) {
0663       for (int iy = 0; iy <= 100; iy++) {
0664         for (int iz = -1; iz <= 1; iz++) {
0665           if (iz == 0)
0666             continue;
0667           if (!EEDetId::validDetId(ix, iy, iz))
0668             continue;
0669 
0670           const EEDetId detid = EEDetId(ix, iy, iz, EEDetId::XYMODE);
0671           EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
0672           int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
0673 
0674           const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(detid);
0675           auto cellGeom = subGeom->getGeometry(detid);
0676           double eta = cellGeom->getPosition().eta();
0677           double phi = cellGeom->getPosition().phi();
0678           double theta = cellGeom->getPosition().theta();
0679 
0680           if (status >= maskedEcalChannelStatusThreshold_) {
0681             std::vector<double> valVec;
0682             std::vector<int> bitVec;
0683             valVec.push_back(eta);
0684             valVec.push_back(phi);
0685             valVec.push_back(theta);
0686             bitVec.push_back(2);
0687             bitVec.push_back(ix);
0688             bitVec.push_back(iy);
0689             bitVec.push_back(iz);
0690             bitVec.push_back(status);
0691             EcalAllDeadChannelsValMap.insert(std::make_pair(detid, valVec));
0692             EcalAllDeadChannelsBitMap.insert(std::make_pair(detid, bitVec));
0693           }
0694         }  // end loop iz
0695       }  // end loop iy
0696     }  // end loop ix
0697   }
0698 
0699   EcalAllDeadChannelsTTMap.clear();
0700   std::map<DetId, std::vector<int> >::iterator bitItor;
0701   for (bitItor = EcalAllDeadChannelsBitMap.begin(); bitItor != EcalAllDeadChannelsBitMap.end(); bitItor++) {
0702     const DetId id = bitItor->first;
0703     EcalTrigTowerDetId ttDetId = ttMap_->towerOf(id);
0704     EcalAllDeadChannelsTTMap.insert(std::make_pair(id, ttDetId));
0705   }
0706 
0707   return 1;
0708 }
0709 
0710 //define this as a plug-in
0711 DEFINE_FWK_MODULE(EcalDeadCellTriggerPrimitiveFilter);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!