CMS 3D CMS Logo

EcalDeadCellBoundaryEnergyFilter.cc
Go to the documentation of this file.
1 // -*- C++ -*-
2 //
3 // Package: EcalDeadCellBoundaryEnergyFilter
4 // Class: EcalDeadCellBoundaryEnergyFilter
5 //
20 //
21 // Original Author: Konstantinos Theofilatos, Ulla Gebbert and Christian Sander
22 // Created: Sat Nov 14 18:43:21 CET 2009
23 //
24 //
25 
26 
27 // system include files
28 #include <memory>
29 #include <fstream>
30 
31 // user include files
34 
37 
39 
45 
48 
56 
59 
61  public:
64 
65  private:
66  void beginJob() override;
67  bool filter(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
68  void endJob() override;
69 
70  // ----------member data ---------------------------
71  const int kMAX;
72 
75 
77  const bool taggingMode_;
78 
79  const bool skimGap_;
80  const bool skimDead_;
81 
83 
86 
88 
90  const std::vector<int> limitDeadCellToChannelStatusEB_;
91  const std::vector<int> limitDeadCellToChannelStatusEE_;
92 
93  const bool enableGap_;
94  const bool debug_;
95 
96 };
97 
98 //
99 // static data member definitions
100 //
101 
102 //
103 // constructors and destructor
104 //
106  : kMAX (50)
107  //now do what ever initialization is needed
108  , EBRecHitsToken_ (consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag> ("recHitsEB")))
109  , EERecHitsToken_ (consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag> ("recHitsEE")))
110 
111  , FilterAlgo_ (iConfig.getUntrackedParameter<std::string> ("FilterAlgo", "TuningMode"))
112  , taggingMode_ (iConfig.getParameter<bool>("taggingMode"))
113  , skimGap_ (iConfig.getUntrackedParameter<bool> ("skimGap", false))
114  , skimDead_ (iConfig.getUntrackedParameter<bool> ("skimDead", false))
115  , cutBoundEnergyGapEE (iConfig.getUntrackedParameter<double> ("cutBoundEnergyGapEE"))
116  , cutBoundEnergyGapEB (iConfig.getUntrackedParameter<double> ("cutBoundEnergyGapEB"))
117  , cutBoundEnergyDeadCellsEB (iConfig.getUntrackedParameter<double> ("cutBoundEnergyDeadCellsEB"))
118  , cutBoundEnergyDeadCellsEE (iConfig.getUntrackedParameter<double> ("cutBoundEnergyDeadCellsEE"))
119 
120  , limitFilterToEB_ (iConfig.getUntrackedParameter<bool> ("limitFilterToEB", false))
121  , limitFilterToEE_ (iConfig.getUntrackedParameter<bool> ("limitFilterToEE", false))
122  , limitDeadCellToChannelStatusEB_ (iConfig.getParameter<std::vector<int> > ("limitDeadCellToChannelStatusEB"))
123  , limitDeadCellToChannelStatusEE_ (iConfig.getParameter<std::vector<int> > ("limitDeadCellToChannelStatusEE"))
124 
125  , enableGap_ (iConfig.getUntrackedParameter<bool> ("enableGap", false))
126  , debug_ (iConfig.getParameter<bool>("debug"))
127 {
128 
130 
131  if (skimGap_ && debug_ ) edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Skim Gap!";
132  if (skimDead_ && debug_ ) edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Skim Dead!";
133 
134  if( debug_ ) {
135  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "Constructor EcalAnomalousEvent";
138  }
139 
140  produces<bool>();
141 
142  produces<AnomalousECALVariables> ("anomalousECALVariables");
143 }
144 
146 }
147 
148 // ------------ method called on each new Event ------------
150  using namespace edm;
151 
152  //int eventno = (int) iEvent.eventAuxiliary().event();
153 
154  std::vector<BoundaryInformation> v_enNeighboursGap_EB;
155  std::vector<BoundaryInformation> v_enNeighboursGap_EE;
156  v_enNeighboursGap_EB.reserve(50);
157  v_enNeighboursGap_EE.reserve(50);
158 
159  std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EB;
160  std::vector<BoundaryInformation> v_boundaryInfoDeadCells_EE;
161  v_boundaryInfoDeadCells_EB.reserve(50);
162  v_boundaryInfoDeadCells_EE.reserve(50);
163 
164  // Get the Ecal RecHits
166  iEvent.getByToken(EBRecHitsToken_, EBRecHits);
168  iEvent.getByToken(EERecHitsToken_, EERecHits);
169 
170  edm::ESHandle<CaloTopology> theCaloTopology;
171  iSetup.get<CaloTopologyRecord> ().get(theCaloTopology);
172 
174  iSetup.get<EcalChannelStatusRcd> ().get(ecalStatus);
175 
177  iSetup.get<CaloGeometryRecord> ().get(geometry);
178 
179 // int DeadChannelsCounterEB = 0;
180 // int DeadChannelsCounterEE = 0;
181 
182  int i_EBDead = 0;
183  int i_EEDead = 0;
184  int i_EBGap = 0;
185  int i_EEGap = 0;
186 
187  std::vector<DetId> sameFlagDetIds;
188 
189  bool pass = true;
190 
191  if (!limitFilterToEE_) {
192 
193  if (debug_)
194  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "process EB";
195 
196  for (EcalRecHitCollection::const_iterator hit = EBRecHits->begin(); hit != EBRecHits->end(); ++hit) {
197 
198  bool detIdAlreadyChecked = false;
199  DetId currDetId = (DetId) hit->id();
200  //add limitation to channel stati
201  EcalChannelStatus::const_iterator chit = ecalStatus->find(currDetId);
202  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
203  if (status != 0)
204  continue;
205  bool passChannelLimitation = false;
206 
207  // check if hit has a dead neighbour
208  std::vector<int> deadNeighbourStati;
209  ebBoundaryCalc.checkRecHitHasDeadNeighbour(*hit, ecalStatus, deadNeighbourStati);
210 
211  for (int cs = 0; cs < (int) limitDeadCellToChannelStatusEB_.size(); ++cs) {
212  int channelAllowed = limitDeadCellToChannelStatusEB_[cs];
213 
214  for (std::vector<int>::iterator sit = deadNeighbourStati.begin(); sit != deadNeighbourStati.end(); ++sit) {
215  if (channelAllowed == *sit || (channelAllowed < 0 && std::abs(channelAllowed) <= *sit)) {
216  passChannelLimitation = true;
217  break;
218  }
219  }
220  }
221 
222  for (std::vector<DetId>::iterator it = sameFlagDetIds.begin(); it != sameFlagDetIds.end(); it++) {
223  if (currDetId == *it)
224  detIdAlreadyChecked = true;
225  }
226 
227  // RecHit is at EB boundary and should be processed
228  if (!detIdAlreadyChecked && deadNeighbourStati.empty() && ebBoundaryCalc.checkRecHitHasInvalidNeighbour(
229  *hit, ecalStatus)) {
230 
232  (const edm::Handle<EcalRecHitCollection>&) EBRecHits, (const EcalRecHit *) &(*hit), theCaloTopology,
233  ecalStatus, geometry);
234 
235  // get rechits along gap cluster
236  for (std::vector<DetId>::iterator it = gapinfo.detIds.begin(); it != gapinfo.detIds.end(); it++) {
237  sameFlagDetIds.push_back(*it);
238  }
239 
240  if (gapinfo.boundaryEnergy > cutBoundEnergyGapEB) {
241 
242  i_EBGap++;
243  v_enNeighboursGap_EB.push_back(gapinfo);
244 
245  if (debug_)
246  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "EB: gap cluster energy: " << gapinfo.boundaryEnergy << " deadCells: "
247  << gapinfo.detIds.size();
248  }
249  }
250 
251  // RecHit is member of boundary and should be processed
252  if (!detIdAlreadyChecked && (passChannelLimitation || (limitDeadCellToChannelStatusEB_.empty()
253  && !deadNeighbourStati.empty()))) {
254 
256  (const edm::Handle<EcalRecHitCollection>&) EBRecHits, (const EcalRecHit *) &(*hit), theCaloTopology,
257  ecalStatus, geometry);
258 
259  // get boundary of !kDead rechits arround the dead cluster
260  for (std::vector<DetId>::iterator it = boundinfo.detIds.begin(); it != boundinfo.detIds.end(); it++) {
261  sameFlagDetIds.push_back(*it);
262  }
263 
264  if (boundinfo.boundaryEnergy > cutBoundEnergyDeadCellsEB) {
265 
266  i_EBDead++;
267  v_boundaryInfoDeadCells_EB.push_back(boundinfo);
268 
269  if (debug_)
270  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "EB: boundary Energy dead RecHit: " << boundinfo.boundaryEnergy << " ET: "
271  << boundinfo.boundaryET << " deadCells: " << boundinfo.detIds.size();
272  }
273 
274  }
275  }
276 
277  }
278 
279  sameFlagDetIds.clear();
280 
281  if (!limitFilterToEB_) {
282 
283  if (debug_)
284  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "process EE";
285 
286  for (EcalRecHitCollection::const_iterator hit = EERecHits->begin(); hit != EERecHits->end(); ++hit) {
287 
288  bool detIdAlreadyChecked = false;
289  DetId currDetId = (DetId) hit->id();
290  //add limitation to channel stati
291  EcalChannelStatus::const_iterator chit = ecalStatus->find(currDetId);
292  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
293  if (status != 0)
294  continue;
295  bool passChannelLimitation = false;
296 
297  // check if hit has a dead neighbour
298  std::vector<int> deadNeighbourStati;
299  eeBoundaryCalc.checkRecHitHasDeadNeighbour(*hit, ecalStatus, deadNeighbourStati);
300 
301  for (int cs = 0; cs < (int) limitDeadCellToChannelStatusEE_.size(); ++cs) {
302  int channelAllowed = limitDeadCellToChannelStatusEE_[cs];
303 
304  for (std::vector<int>::iterator sit = deadNeighbourStati.begin(); sit != deadNeighbourStati.end(); ++sit) {
305  if (channelAllowed == *sit || (channelAllowed < 0 && std::abs(channelAllowed) <= *sit)) {
306  passChannelLimitation = true;
307  break;
308  }
309  }
310  }
311 
312  for (std::vector<DetId>::iterator it = sameFlagDetIds.begin(); it != sameFlagDetIds.end(); it++) {
313  if (currDetId == *it)
314  detIdAlreadyChecked = true;
315  }
316 
317  // RecHit is at EE boundary and should be processed
318  const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(currDetId);
319  auto cellGeom = subGeom->getGeometry(currDetId);
320  double eta = cellGeom->getPosition().eta();
321 
322  if (!detIdAlreadyChecked && deadNeighbourStati.empty() && eeBoundaryCalc.checkRecHitHasInvalidNeighbour(
323  *hit, ecalStatus) && std::abs(eta) < 1.6) {
324 
326  (const edm::Handle<EcalRecHitCollection>&) EERecHits, (const EcalRecHit *) &(*hit), theCaloTopology,
327  ecalStatus, geometry);
328 
329  // get rechits along gap cluster
330  for (std::vector<DetId>::iterator it = gapinfo.detIds.begin(); it != gapinfo.detIds.end(); it++) {
331  sameFlagDetIds.push_back(*it);
332  }
333 
334  if (gapinfo.boundaryEnergy > cutBoundEnergyGapEE) {
335 
336  i_EEGap++;
337  v_enNeighboursGap_EE.push_back(gapinfo);
338 
339  if (debug_)
340  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "EE: gap cluster energy: " << gapinfo.boundaryEnergy << " deadCells: "
341  << gapinfo.detIds.size();
342  }
343  }
344 
345  // RecHit is member of boundary and should be processed
346  if (!detIdAlreadyChecked && (passChannelLimitation || (limitDeadCellToChannelStatusEE_.empty()
347  && !deadNeighbourStati.empty()))) {
348 
350  (const edm::Handle<EcalRecHitCollection>&) EERecHits, (const EcalRecHit *) &(*hit), theCaloTopology,
351  ecalStatus, geometry);
352 
353  // get boundary of !kDead rechits arround the dead cluster
354  for (std::vector<DetId>::iterator it = boundinfo.detIds.begin(); it != boundinfo.detIds.end(); it++) {
355  sameFlagDetIds.push_back(*it);
356  }
357 
358  if (boundinfo.boundaryEnergy > cutBoundEnergyDeadCellsEE) {
359 
360  i_EEDead++;
361  v_boundaryInfoDeadCells_EE.push_back(boundinfo);
362 
363  if (debug_)
364  edm::LogInfo("EcalDeadCellBoundaryEnergyFilter") << "EE: boundary Energy dead RecHit: " << boundinfo.boundaryEnergy << " ET: "
365  << boundinfo.boundaryET << " deadCells: " << boundinfo.detIds.size();
366 
367  }
368 
369  }
370  }
371 
372  }
373 
374  sameFlagDetIds.clear();
375 
376  auto pAnomalousECALVariables = std::make_unique<AnomalousECALVariables>(v_enNeighboursGap_EB,
377  v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE);
378 
379 
380  bool isGap = pAnomalousECALVariables->isGapEcalCluster(cutBoundEnergyGapEB, cutBoundEnergyGapEE);
381  bool isBoundary = pAnomalousECALVariables->isDeadEcalCluster(maxBoundaryEnergy_, limitDeadCellToChannelStatusEB_,
383  pass = (!isBoundary && ((!isGap && enableGap_) || !enableGap_));
384 
385  iEvent.put(std::move(pAnomalousECALVariables), "anomalousECALVariables");
386 
387  iEvent.put(std::make_unique<bool>(pass));
388 
389  if( taggingMode_ ){
390  if (skimDead_ && (i_EBDead >= 1 || i_EEDead >= 1)) {
391  return true;
392  } else if (skimGap_ && (i_EBGap >= 1 || i_EEGap >= 1)) {
393  return true;
394  } else if (!skimDead_ && !skimGap_)
395  return true;
396  else {
397  return false;
398  }
399  }
400  else return pass;
401 
402 /*
403  if (FilterAlgo_ == "TuningMode") {
404  std::unique_ptr<AnomalousECALVariables> pAnomalousECALVariables(new AnomalousECALVariables(v_enNeighboursGap_EB,
405  v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE));
406  iEvent.put(std::move(pAnomalousECALVariables), "anomalousECALVariables");
407 
408  if (skimDead_ && (i_EBDead >= 1 || i_EEDead >= 1)) {
409  return true;
410  } else if (skimGap_ && (i_EBGap >= 1 || i_EEGap >= 1)) {
411  return true;
412  } else if (!skimDead_ && !skimGap_)
413  return true;
414  else {
415  return false;
416  }
417  }
418 
419  if (FilterAlgo_ == "FilterMode") {
420  std::unique_ptr<AnomalousECALVariables> pAnomalousECALVariables(new AnomalousECALVariables(v_enNeighboursGap_EB,
421  v_enNeighboursGap_EE, v_boundaryInfoDeadCells_EB, v_boundaryInfoDeadCells_EE));
422 
423  bool isGap = pAnomalousECALVariables->isGapEcalCluster(cutBoundEnergyGapEB, cutBoundEnergyGapEE);
424  bool isBoundary = pAnomalousECALVariables->isDeadEcalCluster(maxBoundaryEnergy_, limitDeadCellToChannelStatusEB_,
425  limitDeadCellToChannelStatusEE_);
426 
427  bool result = (!isBoundary && ((!isGap && enableGap_) || !enableGap_));
428  if (!result) {
429  }
430  return result;
431  }
432 */
433 
434 // return true;
435 }
436 
437 // ------------ method called once each job just before starting event loop ------------
439 }
440 
441 // ------------ method called once each job just after ending the event loop ------------
443 }
444 
445 //define this as a plug-in
EcalBoundaryInfoCalculator< EBDetId > ebBoundaryCalc
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
Definition: CaloGeometry.cc:49
bool checkRecHitHasInvalidNeighbour(const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus) const
BoundaryInformation gapRecHits(const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry) const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:125
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:517
edm::EDGetTokenT< EcalRecHitCollection > EERecHitsToken_
unique_ptr< ClusterSequence > cs
std::vector< EcalRecHit >::const_iterator const_iterator
EcalBoundaryInfoCalculator< EEDetId > eeBoundaryCalc
EDGetTokenT< ProductType > consumes(edm::InputTag const &tag)
int iEvent
Definition: GenABIO.cc:224
EcalDeadCellBoundaryEnergyFilter(const edm::ParameterSet &)
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
BoundaryInformation boundaryRecHits(const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry) const
std::vector< DetId > detIds
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
unsigned int id
const_iterator end() const
Definition: DetId.h:18
bool filter(edm::StreamID, edm::Event &, const edm::EventSetup &) const override
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
std::vector< Item >::const_iterator const_iterator
ESHandle< TrackerGeometry > geometry
HLT enums.
edm::EDGetTokenT< EcalRecHitCollection > EBRecHitsToken_
T get() const
Definition: EventSetup.h:71
const_iterator find(uint32_t rawId) const
const_iterator end() const
def move(src, dest)
Definition: eostools.py:511
const_iterator begin() const
bool checkRecHitHasDeadNeighbour(const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus, std::vector< int > &stati) const