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:92
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:125
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:68

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

◆ 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:122

◆ 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:125

◆ 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