CMS 3D CMS Logo
List of all members | Public Member Functions | Private Member Functions | Private Attributes
EcalBoundaryInfoCalculator< EcalDetId > Class Template Reference

#include <EcalBoundaryInfoCalculator.h>

Public Member Functions

BoundaryInformation boundaryRecHits (const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
 
bool checkRecHitHasDeadNeighbour (const EcalRecHit &hit, const EcalChannelStatus &ecalStatus, std::vector< int > &stati) const
 
bool checkRecHitHasInvalidNeighbour (const EcalRecHit &hit, const EcalChannelStatus &ecalStatus) const
 
 EcalBoundaryInfoCalculator ()
 
BoundaryInformation gapRecHits (const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
 
void setDebugMode ()
 
 ~EcalBoundaryInfoCalculator ()
 

Private Member Functions

CdOrientation goBackOneCell (CdOrientation currDirection, EcalDetId prev, CaloNavigator< EcalDetId > *theEcalNav) const
 
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator (DetId startE, const CaloTopology &theCaloTopology, EcalSubdetector ecalSubDet) const
 
EcalDetId makeStepInDirection (CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
 
CdOrientation turnLeft (CdOrientation currDirection, bool reverseOrientation) const
 
CdOrientation turnRight (CdOrientation currDirection, bool reverseOrientation) const
 

Private Attributes

bool debug
 
std::map< CdOrientation, CdOrientation > nextDirs
 
std::map< CdOrientation, CdOrientation > oppositeDirs
 
std::map< CdOrientation, CdOrientation > prevDirs
 

Detailed Description

template<class EcalDetId>
class EcalBoundaryInfoCalculator< EcalDetId >

Definition at line 26 of file EcalBoundaryInfoCalculator.h.

Constructor & Destructor Documentation

◆ EcalBoundaryInfoCalculator()

template<class EcalDetId >
EcalBoundaryInfoCalculator< EcalDetId >::EcalBoundaryInfoCalculator ( )

Definition at line 291 of file EcalBoundaryInfoCalculator.h.

291  {
292  nextDirs.clear();
293  nextDirs[north] = east;
294  nextDirs[east] = south;
295  nextDirs[south] = west;
296  nextDirs[west] = north;
297 
298  prevDirs.clear();
299  prevDirs[north] = west;
300  prevDirs[east] = north;
301  prevDirs[south] = east;
302  prevDirs[west] = south;
303 
304  oppositeDirs.clear();
305  oppositeDirs[north] = south;
306  oppositeDirs[south] = north;
307  oppositeDirs[east] = west;
308  oppositeDirs[west] = east;
309 
310  debug = false;
311 }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:284
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:285
EcalBoundaryInfoCalculator::oppositeDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:286

◆ ~EcalBoundaryInfoCalculator()

template<class EcalDetId >
EcalBoundaryInfoCalculator< EcalDetId >::~EcalBoundaryInfoCalculator ( )

Definition at line 314 of file EcalBoundaryInfoCalculator.h.

314 {}

Member Function Documentation

◆ boundaryRecHits()

template<class EcalDetId >
BoundaryInformation EcalBoundaryInfoCalculator< EcalDetId >::boundaryRecHits ( const EcalRecHitCollection &  RecHits,
const EcalRecHit *  hit,
const CaloTopology &  theCaloTopology,
const EcalChannelStatus &  ecalStatus,
const CaloGeometry &  geometry 
) const

Definition at line 317 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

321  {
322  //initialize boundary information
323  std::vector<EcalRecHit> boundaryRecHits;
324  std::vector<DetId> boundaryDetIds;
325  std::vector<int> stati;
326 
327  double boundaryEnergy = 0;
328  double boundaryET = 0;
329  int beCellCounter = 0;
330  bool nextToBorder = false;
331 
332  boundaryRecHits.push_back(*hit);
333  ++beCellCounter;
334  boundaryEnergy += hit->energy();
335  EcalDetId hitdetid = (EcalDetId)hit->id();
336  boundaryDetIds.push_back(hitdetid);
337  const CaloSubdetectorGeometry* subGeom = geometry.getSubdetectorGeometry(hitdetid);
338  auto cellGeom = subGeom->getGeometry(hitdetid);
339  double eta = cellGeom->getPosition().eta();
340  boundaryET += hit->energy() / cosh(eta);
341 
342  if (debug) {
343  edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Boundary RecHits...";
344 
345  if (hitdetid.subdet() == EcalBarrel) {
346  edm::LogInfo("EcalBoundaryInfoCalculator")
347  << "Starting at : (" << ((EBDetId)hitdetid).ieta() << "," << ((EBDetId)hitdetid).iphi() << ")";
348  } else if (hitdetid.subdet() == EcalEndcap) {
349  edm::LogInfo("EcalBoundaryInfoCalculator")
350  << "Starting at : (" << ((EEDetId)hitdetid).ix() << "," << ((EEDetId)hitdetid).iy() << ")";
351  }
352  }
353 
354  //initialize navigator
355  auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
356  CdOrientation currDirection = north;
357  bool reverseOrientation = false;
358 
359  EcalDetId next(0);
360  EcalDetId start = hitdetid;
361  EcalDetId current = start;
362  int current_status = 0;
363 
364  // Search until a dead cell is ahead
365  bool startAlgo = false;
366  int noDirs = 0;
367  while (!startAlgo) {
368  next = makeStepInDirection(currDirection, theEcalNav.get());
369  theEcalNav->setHome(current);
370  theEcalNav->home();
371  EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
372  int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
373  if (status > 0) {
374  stati.push_back(status);
375  startAlgo = true;
376  break;
377  }
378  currDirection = turnLeft(currDirection, reverseOrientation);
379  ++noDirs;
380  if (noDirs > 4) {
381  throw cms::Exception("NoStartingDirection") << "EcalBoundaryInfoCalculator: No starting direction can be found: "
382  "This should never happen if RecHit has a dead neighbour!";
383  }
384  }
385 
386  // go around dead clusters counter clock wise
387  currDirection = turnRight(currDirection, reverseOrientation);
388 
389  // Search for next boundary element
390  bool nextIsStart = false;
391  bool atBorder = false;
392 
393  while (!nextIsStart) {
394  bool nextStepFound = false;
395  int status = 0;
396  noDirs = 0;
397  while (!nextStepFound) {
398  next = makeStepInDirection(currDirection, theEcalNav.get());
399  theEcalNav->setHome(current);
400  theEcalNav->home();
401  EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
402  status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
403  if (status > 0) {
404  // New dead cell found: update status std::vector of dead channels
405  bool present = false;
406  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
407  if (*s == status) {
408  present = true;
409  break;
410  }
411  }
412  if (!present)
413  stati.push_back(status);
414 
415  if (atBorder) {
416  nextStepFound = true;
417  } else {
418  currDirection = turnRight(currDirection, reverseOrientation);
419  }
420  } else if (next == EcalDetId(0)) {
421  // In case the Ecal border is reached -> go along dead cells
422  currDirection = turnLeft(currDirection, reverseOrientation);
423  atBorder = true;
424  } else if (status == 0) {
425  nextStepFound = true;
426  }
427  ++noDirs;
428  if (noDirs > 4) {
429  throw cms::Exception("NoNextDirection")
430  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
431  }
432  }
433 
434  // make next step
435  next = makeStepInDirection(currDirection, theEcalNav.get());
436 
437  if (next == start) {
438  nextIsStart = true;
439  if (debug)
440  edm::LogInfo("EcalBoundaryInfoCalculator") << "Boundary path reached starting position!";
441  }
442 
443  if (debug)
444  edm::LogInfo("EcalBoundaryInfoCalculator")
445  << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
446 
447  // save recHits and add energy if on the boundary (and not inside at border)
448  if ((!atBorder || status == 0) && !nextIsStart) {
449  boundaryDetIds.push_back(next);
450  if (RecHits.find(next) != RecHits.end() && status == 0) {
451  EcalRecHit nexthit = *RecHits.find(next);
452  ++beCellCounter;
453  boundaryRecHits.push_back(nexthit);
454  boundaryEnergy += nexthit.energy();
455  cellGeom = subGeom->getGeometry(hitdetid);
456  eta = cellGeom->getPosition().eta();
457  boundaryET += nexthit.energy() / cosh(eta);
458  }
459  }
460 
461  if (current_status == 0 && status == 0 && atBorder) {
462  // this is for a special case, where dead cells are at border corner
463  currDirection = turnRight(currDirection, reverseOrientation);
464  } else {
465  // if dead region along a border is left, turn left
466  if (status == 0 && atBorder) {
467  atBorder = false;
468  currDirection = turnLeft(currDirection, reverseOrientation);
469  }
470  if (status == 0) {
471  // if outside the cluster turn left to follow boundary
472  currDirection = turnLeft(currDirection, reverseOrientation);
473  } else {
474  // else turn right to check if dead region can be left
475  currDirection = turnRight(currDirection, reverseOrientation);
476  }
477  }
478 
479  // save currect position
480  current = next;
481  current_status = status;
482  }
483 
484  if (debug) {
485  edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Boundary object <<<<<<<<<<<<<<<";
486  edm::LogInfo("EcalBoundaryInfoCalculator") << "no of neighbouring RecHits: " << boundaryRecHits.size();
487  edm::LogInfo("EcalBoundaryInfoCalculator") << "no of neighbouring DetIds: " << boundaryDetIds.size();
488  edm::LogInfo("EcalBoundaryInfoCalculator") << "boundary energy: " << boundaryEnergy;
489  edm::LogInfo("EcalBoundaryInfoCalculator") << "boundary ET: " << boundaryET;
490  edm::LogInfo("EcalBoundaryInfoCalculator") << "no of cells contributing to boundary energy: " << beCellCounter;
491  edm::LogInfo("EcalBoundaryInfoCalculator") << "Channel stati: ";
492  for (std::vector<int>::iterator it = stati.begin(); it != stati.end(); ++it) {
493  edm::LogInfo("EcalBoundaryInfoCalculator") << *it << " ";
494  }
495  edm::LogInfo("EcalBoundaryInfoCalculator");
496  }
497 
498  BoundaryInformation boundInfo;
499  boundInfo.subdet = hitdetid.subdet();
500  boundInfo.detIds = boundaryDetIds;
501  boundInfo.recHits = boundaryRecHits;
502  boundInfo.boundaryEnergy = boundaryEnergy;
503  boundInfo.boundaryET = boundaryET;
504  boundInfo.nextToBorder = nextToBorder;
505  boundInfo.channelStatus = stati;
506 
507  return boundInfo;
508 }
start
Definition: start.py:1
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:239
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
command_line.start
start
Definition: command_line.py:165
BoundaryInformation::detIds
std::vector< DetId > detIds
Definition: BoundaryInformation.h:25
mps_update.status
status
Definition: mps_update.py:68
BoundaryInformation::channelStatus
std::vector< int > channelStatus
Definition: BoundaryInformation.h:26
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:53
CaloSubdetectorGeometry::getGeometry
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.
Definition: CaloSubdetectorGeometry.cc:36
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:91
EcalBoundaryInfoCalculator::boundaryRecHits
BoundaryInformation boundaryRecHits(const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
Definition: EcalBoundaryInfoCalculator.h:317
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
Definition: EcalBoundaryInfoCalculator.h:202
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:19
EcalBoundaryInfoCalculator::initializeEcalNavigator
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator(DetId startE, const CaloTopology &theCaloTopology, EcalSubdetector ecalSubDet) const
Definition: EcalBoundaryInfoCalculator.h:267
BoundaryInformation::subdet
EcalSubdetector subdet
Definition: BoundaryInformation.h:29
BoundaryInformation::recHits
std::vector< EcalRecHit > recHits
Definition: BoundaryInformation.h:23
EcalBoundaryInfoCalculator::turnLeft
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:253
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:240
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:74
EcalRecHit::energy
float energy() const
Definition: EcalRecHit.h:69

◆ checkRecHitHasDeadNeighbour()

template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasDeadNeighbour ( const EcalRecHit &  hit,
const EcalChannelStatus &  ecalStatus,
std::vector< int > &  stati 
) const
inline

Definition at line 43 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

45  {
46  stati.clear();
47  EcalDetId hitdetid = EcalDetId(hit.id());
48 
49  if (hitdetid.subdet() == EcalBarrel) {
50  EBDetId ebhitdetid = (EBDetId)hitdetid;
51 
52  int hitIeta = ebhitdetid.ieta();
53  int hitIphi = ebhitdetid.iphi();
54 
55  for (int ieta = -1; ieta <= 1; ieta++) {
56  for (int iphi = -1; iphi <= 1; iphi++) {
57  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
58  //if (iphi == 0 && ieta == 0)
59  continue;
60  int neighbourIeta = hitIeta + ieta;
61  int neighbourIphi = hitIphi + iphi;
62  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
63  if (neighbourIphi < 1)
64  neighbourIphi += 360;
65  if (neighbourIphi > 360)
66  neighbourIphi -= 360;
67  if (neighbourIeta == 0) {
68  neighbourIeta += ieta;
69  }
70  }
71 
72  if (EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
73  const EBDetId detid = EBDetId(neighbourIeta, neighbourIphi, EBDetId::ETAPHIMODE);
74  EcalChannelStatus::const_iterator chit = ecalStatus.find(detid);
75  int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
76 
77  if (status > 0) {
78  bool present = false;
79  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
80  if (*s == status) {
81  present = true;
82  break;
83  }
84  }
85  if (!present)
86  stati.push_back(status);
87  }
88  }
89  }
90  }
91 
92  } else if (hitdetid.subdet() == EcalEndcap) {
93  EEDetId eehitdetid = (EEDetId)hitdetid;
94  int hitIx = eehitdetid.ix();
95  int hitIy = eehitdetid.iy();
96  int hitIz = eehitdetid.zside();
97 
98  for (int ix = -1; ix <= 1; ix++) {
99  for (int iy = -1; iy <= 1; iy++) {
100  if ((ix == 0 && iy == 0) || ix * iy != 0)
101  //if (ix == 0 && iy == 0)
102  continue;
103  int neighbourIx = hitIx + ix;
104  int neighbourIy = hitIy + iy;
105 
106  if (EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
107  const EEDetId detid = EEDetId(neighbourIx, neighbourIy, hitIz, EEDetId::XYMODE);
108  EcalChannelStatus::const_iterator chit = ecalStatus.find(detid);
109  int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
110 
111  if (status > 0) {
112  bool present = false;
113  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
114  if (*s == status) {
115  present = true;
116  break;
117  }
118  }
119  if (!present)
120  stati.push_back(status);
121  }
122  }
123  }
124  }
125 
126  } else {
127  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
128  }
129 
130  if (!stati.empty())
131  return true;
132  return false;
133  }
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
EEDetId::XYMODE
static const int XYMODE
Definition: EEDetId.h:335
EEDetId::ix
int ix() const
Definition: EEDetId.h:77
EBDetId::validDetId
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.h:118
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
hcalRecHitTable_cff.ieta
ieta
Definition: hcalRecHitTable_cff.py:15
mps_update.status
status
Definition: mps_update.py:68
hcalRecHitTable_cff.iphi
iphi
Definition: hcalRecHitTable_cff.py:16
EBDetId::ETAPHIMODE
static const int ETAPHIMODE
Definition: EBDetId.h:158
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:53
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:91
EEDetId::zside
int zside() const
Definition: EEDetId.h:71
EEDetId::validDetId
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:19
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:74
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128
EEDetId::iy
int iy() const
Definition: EEDetId.h:83

◆ checkRecHitHasInvalidNeighbour()

template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasInvalidNeighbour ( const EcalRecHit &  hit,
const EcalChannelStatus &  ecalStatus 
) const
inline

Definition at line 135 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

135  {
137 
138  EcalDetId hitdetid = EcalDetId(hit.id());
139 
140  if (hitdetid.subdet() == EcalBarrel) {
141  EBDetId ebhitdetid = (EBDetId)hitdetid;
142 
143  int hitIeta = ebhitdetid.ieta();
144  int hitIphi = ebhitdetid.iphi();
145 
146  for (int ieta = -1; ieta <= 1; ieta++) {
147  for (int iphi = -1; iphi <= 1; iphi++) {
148  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
149  //if (iphi == 0 && ieta == 0)
150  continue;
151  int neighbourIeta = hitIeta + ieta;
152  int neighbourIphi = hitIphi + iphi;
153  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
154  if (neighbourIphi < 1)
155  neighbourIphi += 360;
156  if (neighbourIphi > 360)
157  neighbourIphi -= 360;
158  if (neighbourIeta == 0) {
159  neighbourIeta += ieta;
160  }
161  }
162 
163  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
164  return true;
165  }
166  }
167  }
168 
169  } else if (hitdetid.subdet() == EcalEndcap) {
170  EEDetId eehitdetid = (EEDetId)hitdetid;
171  int hitIx = eehitdetid.ix();
172  int hitIy = eehitdetid.iy();
173  int hitIz = eehitdetid.zside();
174 
175  for (int ix = -1; ix <= 1; ix++) {
176  for (int iy = -1; iy <= 1; iy++) {
177  if ((ix == 0 && iy == 0) || ix * iy != 0)
178  //if (ix == 0 && iy == 0)
179  continue;
180  int neighbourIx = hitIx + ix;
181  int neighbourIy = hitIy + iy;
182 
183  if (!EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
184  return true;
185  }
186  }
187  }
188 
189  } else {
190  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
191  }
192 
193  return false;
194  }
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
EEDetId::ix
int ix() const
Definition: EEDetId.h:77
EBDetId::validDetId
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.h:118
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
hcalRecHitTable_cff.ieta
ieta
Definition: hcalRecHitTable_cff.py:15
hcalRecHitTable_cff.iphi
iphi
Definition: hcalRecHitTable_cff.py:16
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:91
EEDetId::zside
int zside() const
Definition: EEDetId.h:71
EEDetId::validDetId
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128
EEDetId::iy
int iy() const
Definition: EEDetId.h:83

◆ gapRecHits()

template<class EcalDetId >
BoundaryInformation EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits ( const EcalRecHitCollection &  RecHits,
const EcalRecHit *  hit,
const CaloTopology &  theCaloTopology,
const EcalChannelStatus &  ecalStatus,
const CaloGeometry &  geometry 
) const

Definition at line 511 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

515  {
516  //initialize boundary information
517  std::vector<EcalRecHit> gapRecHits;
518  std::vector<DetId> gapDetIds;
519 
520  double gapEnergy = 0;
521  double gapET = 0;
522  bool nextToBorder = false;
523 
524  gapRecHits.push_back(*hit);
525  gapEnergy += hit->energy();
526  EcalDetId hitdetid = (EcalDetId)hit->id();
527  gapDetIds.push_back(hitdetid);
528  const CaloSubdetectorGeometry* subGeom = geometry.getSubdetectorGeometry(hitdetid);
529  auto cellGeom = subGeom->getGeometry(hitdetid);
530  double eta = cellGeom->getPosition().eta();
531  gapET += hit->energy() / cosh(eta);
532 
533  if (debug) {
534  edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Border RecHits...";
535 
536  if (hitdetid.subdet() == EcalBarrel) {
537  edm::LogInfo("EcalBoundaryInfoCalculator")
538  << "Starting at : (" << ((EBDetId)hitdetid).ieta() << "," << ((EBDetId)hitdetid).iphi() << ")";
539  } else if (hitdetid.subdet() == EcalEndcap) {
540  edm::LogInfo("EcalBoundaryInfoCalculator")
541  << "Starting at : (" << ((EEDetId)hitdetid).ix() << "," << ((EEDetId)hitdetid).iy() << ")";
542  }
543  }
544 
545  //initialize navigator
546  auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
547  CdOrientation currDirection = north;
548  bool reverseOrientation = false;
549 
550  EcalDetId next(0);
551  EcalDetId start = hitdetid;
552  EcalDetId current = start;
553 
554  // Search until a invalid cell is ahead
555  bool startAlgo = false;
556  int noDirs = 0;
557  while (!startAlgo) {
558  next = makeStepInDirection(currDirection, theEcalNav.get());
559  theEcalNav->setHome(start);
560  theEcalNav->home();
561  if (next == EcalDetId(0)) {
562  startAlgo = true;
563  nextToBorder = true;
564  break;
565  }
566  currDirection = turnLeft(currDirection, reverseOrientation);
567  ++noDirs;
568  if (noDirs > 4) {
569  throw cms::Exception("NoStartingDirection") << "EcalBoundaryInfoCalculator: No starting direction can be found: "
570  "This should never happen if RecHit is at border!";
571  }
572  }
573 
575  CdOrientation startDirection = currDirection;
576  currDirection = turnLeft(currDirection, reverseOrientation);
577 
578  // Search for next border element
579  bool endIsFound = false;
580  bool startIsEnd = false;
581 
582  while (!endIsFound) {
583  bool nextStepFound = false;
584  int status = 0;
585  noDirs = 0;
586  while (!nextStepFound) {
587  next = makeStepInDirection(currDirection, theEcalNav.get());
588  theEcalNav->setHome(current);
589  theEcalNav->home();
590  EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
591  status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
592  if (status > 0) {
593  // Find dead cell along border -> end of cluster
594  endIsFound = true;
595  break;
596  } else if (next == EcalDetId(0)) {
597  // In case the Ecal border -> go along gap
598  currDirection = turnLeft(currDirection, reverseOrientation);
599  } else if (status == 0) {
600  if (RecHits.find(next) != RecHits.end()) {
601  nextStepFound = true;
602  } else {
603  endIsFound = true;
604  break;
605  }
606  }
607  ++noDirs;
608  if (noDirs > 4) {
609  throw cms::Exception("NoNextDirection")
610  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
611  }
612  }
613 
614  // make next step
615  next = makeStepInDirection(currDirection, theEcalNav.get());
616  current = next;
617 
618  if (next == start) {
619  startIsEnd = true;
620  endIsFound = true;
621  if (debug)
622  edm::LogInfo("EcalBoundaryInfoCalculator") << "Path along gap reached starting position!";
623  }
624 
625  if (debug) {
626  edm::LogInfo("EcalBoundaryInfoCalculator")
627  << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
628  if (endIsFound)
629  edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going left";
630  }
631 
632  // save recHits and add energy
633  if (!endIsFound) {
634  gapDetIds.push_back(next);
635  if (RecHits.find(next) != RecHits.end()) {
636  EcalRecHit nexthit = *RecHits.find(next);
637  gapRecHits.push_back(nexthit);
638  gapEnergy += nexthit.energy();
639  cellGeom = subGeom->getGeometry(next);
640  eta = cellGeom->getPosition().eta();
641  gapET += nexthit.energy() / cosh(eta);
642  }
643  }
644 
645  // turn right to follow gap
646  currDirection = turnRight(currDirection, reverseOrientation);
647  }
648 
650  theEcalNav->setHome(start);
651  theEcalNav->home();
652  current = start;
653  currDirection = startDirection;
654  currDirection = turnRight(currDirection, reverseOrientation);
655 
656  // Search for next border element
657  endIsFound = false;
658 
659  if (!startIsEnd) {
660  while (!endIsFound) {
661  bool nextStepFound = false;
662  int status = 0;
663  noDirs = 0;
664  while (!nextStepFound) {
665  next = makeStepInDirection(currDirection, theEcalNav.get());
666  theEcalNav->setHome(current);
667  theEcalNav->home();
668  EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
669  status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
670  if (status > 0) {
671  // Find dead cell along border -> end of cluster
672  endIsFound = true;
673  break;
674  } else if (next == EcalDetId(0)) {
675  // In case the Ecal border -> go along gap
676  currDirection = turnRight(currDirection, reverseOrientation);
677  } else if (status == 0) {
678  if (RecHits.find(next) != RecHits.end()) {
679  nextStepFound = true;
680  } else {
681  endIsFound = true;
682  break;
683  }
684  }
685  ++noDirs;
686  if (noDirs > 4) {
687  throw cms::Exception("NoStartingDirection")
688  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
689  }
690  }
691 
692  // make next step
693  next = makeStepInDirection(currDirection, theEcalNav.get());
694  current = next;
695 
696  if (debug) {
697  edm::LogInfo("EcalBoundaryInfoCalculator")
698  << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
699  if (endIsFound)
700  edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going right";
701  }
702 
703  // save recHits and add energy
704  if (!endIsFound) {
705  gapDetIds.push_back(next);
706  if (RecHits.find(next) != RecHits.end()) {
707  EcalRecHit nexthit = *RecHits.find(next);
708  gapRecHits.push_back(nexthit);
709  gapEnergy += nexthit.energy();
710  cellGeom = subGeom->getGeometry(next);
711  eta = cellGeom->getPosition().eta();
712  gapET += nexthit.energy() / cosh(eta);
713  }
714  }
715 
716  // turn left to follow gap
717  currDirection = turnLeft(currDirection, reverseOrientation);
718  }
719  }
720 
721  if (debug) {
722  edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
723  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of RecHits along gap: " << gapRecHits.size();
724  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of DetIds along gap: " << gapDetIds.size();
725  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap energy: " << gapEnergy;
726  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap ET: " << gapET;
727  }
728 
729  BoundaryInformation gapInfo;
730  gapInfo.subdet = hitdetid.subdet();
731  gapInfo.detIds = gapDetIds;
732  gapInfo.recHits = gapRecHits;
733  gapInfo.boundaryEnergy = gapEnergy;
734  gapInfo.boundaryET = gapET;
735  gapInfo.nextToBorder = nextToBorder;
736  std::vector<int> stati;
737  gapInfo.channelStatus = stati;
738 
739  return gapInfo;
740 }
start
Definition: start.py:1
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:239
EcalBoundaryInfoCalculator::gapRecHits
BoundaryInformation gapRecHits(const EcalRecHitCollection &, const EcalRecHit *, const CaloTopology &theCaloTopology, const EcalChannelStatus &ecalStatus, const CaloGeometry &geometry) const
Definition: EcalBoundaryInfoCalculator.h:511
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
command_line.start
start
Definition: command_line.py:165
BoundaryInformation::detIds
std::vector< DetId > detIds
Definition: BoundaryInformation.h:25
mps_update.status
status
Definition: mps_update.py:68
BoundaryInformation::channelStatus
std::vector< int > channelStatus
Definition: BoundaryInformation.h:26
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:53
CaloSubdetectorGeometry::getGeometry
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.
Definition: CaloSubdetectorGeometry.cc:36
hit::id
unsigned int id
Definition: SiStripHitEffFromCalibTree.cc:91
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:267
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
Definition: EcalBoundaryInfoCalculator.h:202
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:19
EcalBoundaryInfoCalculator::initializeEcalNavigator
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator(DetId startE, const CaloTopology &theCaloTopology, EcalSubdetector ecalSubDet) const
Definition: EcalBoundaryInfoCalculator.h:267
BoundaryInformation::subdet
EcalSubdetector subdet
Definition: BoundaryInformation.h:29
BoundaryInformation::recHits
std::vector< EcalRecHit > recHits
Definition: BoundaryInformation.h:23
EcalBoundaryInfoCalculator::turnLeft
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:253
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:240
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:74
EcalRecHit::energy
float energy() const
Definition: EcalRecHit.h:69

◆ goBackOneCell()

template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::goBackOneCell ( CdOrientation  currDirection,
EcalDetId  prev,
CaloNavigator< EcalDetId > *  theEcalNav 
) const
inlineprivate

Definition at line 227 of file EcalBoundaryInfoCalculator.h.

227  {
228  auto oIt = oppositeDirs.find(currDirection);
229  CdOrientation oppDirection = none;
230  if (oIt != oppositeDirs.end()) {
231  oppDirection = oIt->second;
232  theEcalNav->setHome(prev);
233  }
234  EcalDetId currDetId = theEcalNav->pos();
235 
236  return oppDirection;
237  }
CaloNavigator::pos
T pos() const
get the current position
Definition: CaloNavigator.h:24
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:90
EcalBoundaryInfoCalculator::oppositeDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:286

◆ initializeEcalNavigator()

template<class EcalDetId>
std::unique_ptr<CaloNavigator<EcalDetId> > EcalBoundaryInfoCalculator< EcalDetId >::initializeEcalNavigator ( DetId  startE,
const CaloTopology &  theCaloTopology,
EcalSubdetector  ecalSubDet 
) const
inlineprivate

Definition at line 267 of file EcalBoundaryInfoCalculator.h.

269  {
270  std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav;
271  if (ecalSubDet == EcalBarrel) {
272  theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
273  (EBDetId)startE, (theCaloTopology.getSubdetectorTopology(DetId::Ecal, ecalSubDet)));
274  } else if (ecalSubDet == EcalEndcap) {
275  theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
276  (EEDetId)startE, (theCaloTopology.getSubdetectorTopology(DetId::Ecal, ecalSubDet)));
277  } else {
278  edm::LogWarning("EcalBoundaryInfoCalculator")
279  << "initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
280  }
281  return theEcalNav;
282  }
CaloTopology::getSubdetectorTopology
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
Definition: CaloTopology.cc:17
edm::LogWarning
Log< level::Warning, false > LogWarning
Definition: MessageLogger.h:128

◆ makeStepInDirection()

template<class EcalDetId>
EcalDetId EcalBoundaryInfoCalculator< EcalDetId >::makeStepInDirection ( CdOrientation  direction,
const CaloNavigator< EcalDetId > *  theNavi 
) const
inlineprivate

Definition at line 202 of file EcalBoundaryInfoCalculator.h.

202  {
203  EcalDetId next;
204  switch (direction) {
205  case north: {
206  next = theNavi->north();
207  break;
208  }
209  case east: {
210  next = theNavi->east();
211  break;
212  }
213  case south: {
214  next = theNavi->south();
215  break;
216  }
217  case west: {
218  next = theNavi->west();
219  break;
220  }
221  default:
222  break;
223  }
224  return next;
225  }
CaloNavigator::north
T north() const
move the navigator north
Definition: CaloNavigator.h:30
CaloNavigator::south
T south() const
move the navigator south
Definition: CaloNavigator.h:36
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
CaloNavigator::east
T east() const
move the navigator east
Definition: CaloNavigator.h:42
CaloNavigator::west
T west() const
move the navigator west
Definition: CaloNavigator.h:48

◆ setDebugMode()

template<class EcalDetId>
void EcalBoundaryInfoCalculator< EcalDetId >::setDebugMode ( )
inline

Definition at line 196 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter().

196  {
197  edm::LogInfo("EcalBoundaryInfoCalculator") << "set Debug Mode!";
198  debug = true;
199  }
edm::LogInfo
Log< level::Info, false > LogInfo
Definition: MessageLogger.h:131

◆ turnLeft()

template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnLeft ( CdOrientation  currDirection,
bool  reverseOrientation 
) const
inlineprivate

Definition at line 253 of file EcalBoundaryInfoCalculator.h.

253  {
254  //read nextDirection
255  std::map<CdOrientation, CdOrientation> turnMap = prevDirs;
256  if (reverseOrientation)
257  turnMap = nextDirs;
258  std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
259  CdOrientation nextDirection = none;
260  if (nIt != turnMap.end())
261  nextDirection = (*nIt).second;
262  else
263  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
264  return nextDirection;
265  }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:284
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:285

◆ turnRight()

template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnRight ( CdOrientation  currDirection,
bool  reverseOrientation 
) const
inlineprivate

Definition at line 239 of file EcalBoundaryInfoCalculator.h.

239  {
240  //read nextDirection
241  std::map<CdOrientation, CdOrientation> turnMap = nextDirs;
242  if (reverseOrientation)
243  turnMap = prevDirs;
244  std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
245  CdOrientation nextDirection = none;
246  if (nIt != turnMap.end())
247  nextDirection = (*nIt).second;
248  else
249  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
250  return nextDirection;
251  }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:284
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:285

Member Data Documentation

◆ debug

template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::debug
private

Definition at line 287 of file EcalBoundaryInfoCalculator.h.

Referenced by runTauIdMVA.TauIDEmbedder::loadMVA_WPs_run2_2017(), runTauIdMVA.TauIDEmbedder::runTauID(), and EcalBoundaryInfoCalculator< EBDetId >::setDebugMode().

◆ nextDirs

template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::nextDirs
private

Definition at line 284 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EBDetId >::turnLeft(), and EcalBoundaryInfoCalculator< EBDetId >::turnRight().

◆ oppositeDirs

template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::oppositeDirs
private

Definition at line 286 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EBDetId >::goBackOneCell().

◆ prevDirs

template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::prevDirs
private

Definition at line 285 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EBDetId >::turnLeft(), and EcalBoundaryInfoCalculator< EBDetId >::turnRight().

Generated for CMSSW Reference Manual by doxygen 1.8.14