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 edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry) const
 
bool checkRecHitHasDeadNeighbour (const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus, std::vector< int > &stati) const
 
bool checkRecHitHasInvalidNeighbour (const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus) const
 
 EcalBoundaryInfoCalculator ()
 
BoundaryInformation gapRecHits (const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< 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 edm::ESHandle< 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 27 of file EcalBoundaryInfoCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 292 of file EcalBoundaryInfoCalculator.h.

References EcalBoundaryInfoCalculator< EcalDetId >::debug, east, EcalBoundaryInfoCalculator< EcalDetId >::nextDirs, north, EcalBoundaryInfoCalculator< EcalDetId >::oppositeDirs, EcalBoundaryInfoCalculator< EcalDetId >::prevDirs, south, and west.

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

Definition at line 315 of file EcalBoundaryInfoCalculator.h.

315 {}

Member Function Documentation

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

Definition at line 318 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, EcalBoundaryInfoCalculator< EcalDetId >::debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, EcalCondObjectContainer< T >::end(), edm::SortedCollection< T, SORT >::end(), EcalRecHit::energy(), PVValHelper::eta, Exception, EcalCondObjectContainer< T >::find(), edm::SortedCollection< T, SORT >::find(), CaloSubdetectorGeometry::getGeometry(), CaloGeometry::getSubdetectorGeometry(), EcalRecHit::id(), EcalBoundaryInfoCalculator< EcalDetId >::initializeEcalNavigator(), EcalBoundaryInfoCalculator< EcalDetId >::makeStepInDirection(), GetRecoTauVFromDQM_MC_cff::next, BoundaryInformation::nextToBorder, north, BoundaryInformation::recHits, alignCSCRings::s, command_line::start, mps_update::status, BoundaryInformation::subdet, EcalBoundaryInfoCalculator< EcalDetId >::turnLeft(), and EcalBoundaryInfoCalculator< EcalDetId >::turnRight().

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 44 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 136 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

136  {
138 
139  EcalDetId hitdetid = EcalDetId(hit.id());
140 
141  if (hitdetid.subdet() == EcalBarrel) {
142  EBDetId ebhitdetid = (EBDetId)hitdetid;
143 
144  int hitIeta = ebhitdetid.ieta();
145  int hitIphi = ebhitdetid.iphi();
146 
147  for (int ieta = -1; ieta <= 1; ieta++) {
148  for (int iphi = -1; iphi <= 1; iphi++) {
149  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
150  //if (iphi == 0 && ieta == 0)
151  continue;
152  int neighbourIeta = hitIeta + ieta;
153  int neighbourIphi = hitIphi + iphi;
154  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
155  if (neighbourIphi < 1)
156  neighbourIphi += 360;
157  if (neighbourIphi > 360)
158  neighbourIphi -= 360;
159  if (neighbourIeta == 0) {
160  neighbourIeta += ieta;
161  }
162  }
163 
164  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
165  return true;
166  }
167  }
168  }
169 
170  } else if (hitdetid.subdet() == EcalEndcap) {
171  EEDetId eehitdetid = (EEDetId)hitdetid;
172  int hitIx = eehitdetid.ix();
173  int hitIy = eehitdetid.iy();
174  int hitIz = eehitdetid.zside();
175 
176  for (int ix = -1; ix <= 1; ix++) {
177  for (int iy = -1; iy <= 1; iy++) {
178  if ((ix == 0 && iy == 0) || ix * iy != 0)
179  //if (ix == 0 && iy == 0)
180  continue;
181  int neighbourIx = hitIx + ix;
182  int neighbourIy = hitIy + iy;
183 
184  if (!EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
185  return true;
186  }
187  }
188  }
189 
190  } else {
191  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
192  }
193 
194  return false;
195  }
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::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:51
EEDetId::zside
int zside() const
Definition: EEDetId.h:71
EEDetId::iy
int iy() const
Definition: EEDetId.h:83
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:49
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:77
EEDetId::validDetId
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.h:248
template<class EcalDetId >
BoundaryInformation EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits ( const edm::Handle< EcalRecHitCollection > &  RecHits,
const EcalRecHit *  hit,
const edm::ESHandle< CaloTopology >  theCaloTopology,
const edm::ESHandle< EcalChannelStatus >  ecalStatus,
const edm::ESHandle< CaloGeometry >  geometry 
) const

Definition at line 513 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, EcalBoundaryInfoCalculator< EcalDetId >::debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, EcalCondObjectContainer< T >::end(), edm::SortedCollection< T, SORT >::end(), EcalRecHit::energy(), PVValHelper::eta, Exception, EcalCondObjectContainer< T >::find(), edm::SortedCollection< T, SORT >::find(), CaloSubdetectorGeometry::getGeometry(), CaloGeometry::getSubdetectorGeometry(), EcalRecHit::id(), EcalBoundaryInfoCalculator< EcalDetId >::initializeEcalNavigator(), EcalBoundaryInfoCalculator< EcalDetId >::makeStepInDirection(), GetRecoTauVFromDQM_MC_cff::next, BoundaryInformation::nextToBorder, north, BoundaryInformation::recHits, command_line::start, mps_update::status, BoundaryInformation::subdet, EcalBoundaryInfoCalculator< EcalDetId >::turnLeft(), and EcalBoundaryInfoCalculator< EcalDetId >::turnRight().

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 228 of file EcalBoundaryInfoCalculator.h.

228  {
229  std::map<CdOrientation, CdOrientation>::iterator oIt = oppositeDirs.find(currDirection);
230  CdOrientation oppDirection = none;
231  if (oIt != oppositeDirs.end()) {
232  oppDirection = oIt->second;
233  theEcalNav->setHome(prev);
234  }
235  EcalDetId currDetId = theEcalNav->pos();
236 
237  return oppDirection;
238  }
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:91
CaloNavigator::pos
T pos() const
get the current position
Definition: CaloNavigator.h:25
EcalBoundaryInfoCalculator::oppositeDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:287
template<class EcalDetId>
std::unique_ptr<CaloNavigator<EcalDetId> > EcalBoundaryInfoCalculator< EcalDetId >::initializeEcalNavigator ( DetId  startE,
const edm::ESHandle< CaloTopology >  theCaloTopology,
EcalSubdetector  ecalSubDet 
) const
inlineprivate

Definition at line 268 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EcalDetId >::boundaryRecHits(), and EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits().

270  {
271  std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav(nullptr);
272  if (ecalSubDet == EcalBarrel) {
273  theEcalNav.reset(new CaloNavigator<EcalDetId>(
274  (EBDetId)startE, (theCaloTopology->getSubdetectorTopology(DetId::Ecal, ecalSubDet))));
275  } else if (ecalSubDet == EcalEndcap) {
276  theEcalNav.reset(new CaloNavigator<EcalDetId>(
277  (EEDetId)startE, (theCaloTopology->getSubdetectorTopology(DetId::Ecal, ecalSubDet))));
278  } else {
279  edm::LogWarning("EcalBoundaryInfoCalculator")
280  << "initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
281  }
282  return theEcalNav;
283  }
CaloTopology::getSubdetectorTopology
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
Definition: CaloTopology.cc:17
template<class EcalDetId>
EcalDetId EcalBoundaryInfoCalculator< EcalDetId >::makeStepInDirection ( CdOrientation  direction,
const CaloNavigator< EcalDetId > *  theNavi 
) const
inlineprivate

Definition at line 203 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EcalDetId >::boundaryRecHits(), and EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits().

203  {
204  EcalDetId next;
205  switch (direction) {
206  case north: {
207  next = theNavi->north();
208  break;
209  }
210  case east: {
211  next = theNavi->east();
212  break;
213  }
214  case south: {
215  next = theNavi->south();
216  break;
217  }
218  case west: {
219  next = theNavi->west();
220  break;
221  }
222  default:
223  break;
224  }
225  return next;
226  }
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
CaloNavigator::west
T west() const
move the navigator west
Definition: CaloNavigator.h:49
CaloNavigator::south
T south() const
move the navigator south
Definition: CaloNavigator.h:37
CaloNavigator::east
T east() const
move the navigator east
Definition: CaloNavigator.h:43
CaloNavigator::north
T north() const
move the navigator north
Definition: CaloNavigator.h:31
template<class EcalDetId>
void EcalBoundaryInfoCalculator< EcalDetId >::setDebugMode ( )
inline

Definition at line 197 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter().

197  {
198  edm::LogInfo("EcalBoundaryInfoCalculator") << "set Debug Mode!";
199  debug = true;
200  }
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnLeft ( CdOrientation  currDirection,
bool  reverseOrientation 
) const
inlineprivate

Definition at line 254 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EcalDetId >::boundaryRecHits(), and EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits().

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

Definition at line 240 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EcalDetId >::boundaryRecHits(), and EcalBoundaryInfoCalculator< EcalDetId >::gapRecHits().

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

Member Data Documentation

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

Definition at line 288 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 285 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 287 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalBoundaryInfoCalculator< EcalDetId >::EcalBoundaryInfoCalculator(), and EcalBoundaryInfoCalculator< EBDetId >::goBackOneCell().

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

Definition at line 286 of file EcalBoundaryInfoCalculator.h.

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

Generated for CMSSW Reference Manual by doxygen 1.8.11