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 28 of file EcalBoundaryInfoCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 297 of file EcalBoundaryInfoCalculator.h.

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

297  {
298 
299  nextDirs.clear();
300  nextDirs[north] = east;
301  nextDirs[east] = south;
302  nextDirs[south] = west;
303  nextDirs[west] = north;
304 
305  prevDirs.clear();
306  prevDirs[north] = west;
307  prevDirs[east] = north;
308  prevDirs[south] = east;
309  prevDirs[west] = south;
310 
311  oppositeDirs.clear();
312  oppositeDirs[north] = south;
313  oppositeDirs[south] = north;
314  oppositeDirs[east] = west;
315  oppositeDirs[west] = east;
316 
317  debug = false;
318 
319 }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:290
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:291
EcalBoundaryInfoCalculator::oppositeDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:292
template<class EcalDetId >
EcalBoundaryInfoCalculator< EcalDetId >::~EcalBoundaryInfoCalculator ( )

Definition at line 321 of file EcalBoundaryInfoCalculator.h.

321  {
322 }

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 324 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, EcalBoundaryInfoCalculator< EcalDetId >::debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, EcalCondObjectContainer< T >::end(), EcalRecHit::energy(), PVValHelper::eta, Exception, EcalCondObjectContainer< T >::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().

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

Definition at line 43 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

44  {
45 
46  stati.clear();
47  EcalDetId hitdetid = EcalDetId(hit.id());
48 
49  if (hitdetid.subdet() == EcalBarrel) {
50 
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 
75  const EBDetId detid = EBDetId(neighbourIeta, neighbourIphi, EBDetId::ETAPHIMODE);
76  EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
77  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
78 
79  if (status > 0) {
80  bool present = false;
81  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
82  if (*s == status) {
83  present = true;
84  break;
85  }
86  }
87  if (!present)
88  stati.push_back(status);
89  }
90  }
91  }
92  }
93 
94  } else if (hitdetid.subdet() == EcalEndcap) {
95 
96  EEDetId eehitdetid = (EEDetId) hitdetid;
97  int hitIx = eehitdetid.ix();
98  int hitIy = eehitdetid.iy();
99  int hitIz = eehitdetid.zside();
100 
101  for (int ix = -1; ix <= 1; ix++) {
102  for (int iy = -1; iy <= 1; iy++) {
103  if ((ix == 0 && iy == 0) || ix * iy != 0)
104  //if (ix == 0 && iy == 0)
105  continue;
106  int neighbourIx = hitIx + ix;
107  int neighbourIy = hitIy + iy;
108 
109  if (EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
110 
111  const EEDetId detid = EEDetId(neighbourIx, neighbourIy, hitIz, EEDetId::XYMODE);
112  EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
113  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
114 
115  if (status > 0) {
116  bool present = false;
117  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
118  if (*s == status) {
119  present = true;
120  break;
121  }
122  }
123  if (!present)
124  stati.push_back(status);
125  }
126  }
127  }
128  }
129 
130  } else {
131  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
132  }
133 
134  if (!stati.empty())
135  return true;
136  return false;
137 
138  }
EEDetId::ix
int ix() const
Definition: EEDetId.h:76
EEDetId::XYMODE
static const int XYMODE
Definition: EEDetId.h:339
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:124
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
EEDetId::zside
int zside() const
Definition: EEDetId.h:70
EEDetId::iy
int iy() const
Definition: EEDetId.h:82
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
EBDetId::ETAPHIMODE
static const int ETAPHIMODE
Definition: EBDetId.h:166
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:59
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:88
template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasInvalidNeighbour ( const EcalRecHit &  hit,
const edm::ESHandle< EcalChannelStatus >  ecalStatus 
) const
inline

Definition at line 140 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

140  {
142 
143  EcalDetId hitdetid = EcalDetId(hit.id());
144 
145  if (hitdetid.subdet() == EcalBarrel) {
146 
147  EBDetId ebhitdetid = (EBDetId) hitdetid;
148 
149  int hitIeta = ebhitdetid.ieta();
150  int hitIphi = ebhitdetid.iphi();
151 
152  for (int ieta = -1; ieta <= 1; ieta++) {
153  for (int iphi = -1; iphi <= 1; iphi++) {
154  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
155  //if (iphi == 0 && ieta == 0)
156  continue;
157  int neighbourIeta = hitIeta + ieta;
158  int neighbourIphi = hitIphi + iphi;
159  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
160  if (neighbourIphi < 1)
161  neighbourIphi += 360;
162  if (neighbourIphi > 360)
163  neighbourIphi -= 360;
164  if (neighbourIeta == 0) {
165  neighbourIeta += ieta;
166  }
167  }
168 
169  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
170  return true;
171  }
172  }
173  }
174 
175  } else if (hitdetid.subdet() == EcalEndcap) {
176 
177  EEDetId eehitdetid = (EEDetId) hitdetid;
178  int hitIx = eehitdetid.ix();
179  int hitIy = eehitdetid.iy();
180  int hitIz = eehitdetid.zside();
181 
182  for (int ix = -1; ix <= 1; ix++) {
183  for (int iy = -1; iy <= 1; iy++) {
184  if ((ix == 0 && iy == 0) || ix * iy != 0)
185  //if (ix == 0 && iy == 0)
186  continue;
187  int neighbourIx = hitIx + ix;
188  int neighbourIy = hitIy + iy;
189 
190  if (!EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
191  return true;
192  }
193  }
194  }
195 
196  } else {
197  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
198  }
199 
200  return false;
201  }
EEDetId::ix
int ix() const
Definition: EEDetId.h:76
EBDetId::validDetId
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.h:124
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:53
EEDetId::zside
int zside() const
Definition: EEDetId.h:70
EEDetId::iy
int iy() const
Definition: EEDetId.h:82
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:51
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 515 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, EcalBoundaryInfoCalculator< EcalDetId >::debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, EcalCondObjectContainer< T >::end(), EcalRecHit::energy(), PVValHelper::eta, Exception, EcalCondObjectContainer< T >::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 
519  //initialize boundary information
520  std::vector<EcalRecHit> gapRecHits;
521  std::vector<DetId> gapDetIds;
522 
523  double gapEnergy = 0;
524  double gapET = 0;
525  int gapCellCounter = 0;
526  bool nextToBorder = false;
527 
528  gapRecHits.push_back(*hit);
529  ++gapCellCounter;
530  gapEnergy += hit->energy();
531  EcalDetId hitdetid = (EcalDetId) hit->id();
532  gapDetIds.push_back(hitdetid);
533  const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(hitdetid);
534  auto cellGeom = subGeom->getGeometry(hitdetid);
535  double eta = cellGeom->getPosition().eta();
536  gapET += hit->energy() / cosh(eta);
537 
538  if (debug) {
539  edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Border RecHits...";
540 
541  if (hitdetid.subdet() == EcalBarrel) {
542  edm::LogInfo("EcalBoundaryInfoCalculator") << "Starting at : (" << ((EBDetId) hitdetid).ieta() << "," << ((EBDetId) hitdetid).iphi() << ")";
543  } else if (hitdetid.subdet() == EcalEndcap) {
544  edm::LogInfo("EcalBoundaryInfoCalculator") << "Starting at : (" << ((EEDetId) hitdetid).ix() << "," << ((EEDetId) hitdetid).iy() << ")";
545  }
546  }
547 
548  //initialize navigator
549  auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
550  CdOrientation currDirection = north;
551  bool reverseOrientation = false;
552 
553  EcalDetId next(0);
554  EcalDetId start = hitdetid;
555  EcalDetId current = start;
556 
557  // Search until a invalid cell is ahead
558  bool startAlgo = false;
559  int noDirs = 0;
560  while (!startAlgo) {
561  next = makeStepInDirection(currDirection, theEcalNav.get());
562  theEcalNav->setHome(start);
563  theEcalNav->home();
564  if (next == EcalDetId(0)) {
565  startAlgo = true;
566  nextToBorder = true;
567  break;
568  }
569  currDirection = turnLeft(currDirection, reverseOrientation);
570  ++noDirs;
571  if (noDirs > 4) {
572  throw cms::Exception("NoStartingDirection")
573  << "EcalBoundaryInfoCalculator: No starting direction can be found: This should never happen if RecHit is at border!";
574  }
575  }
576 
578  CdOrientation startDirection = currDirection;
579  currDirection = turnLeft(currDirection, reverseOrientation);
580 
581  // Search for next border element
582  bool endIsFound = false;
583  bool startIsEnd = false;
584 
585  while (!endIsFound) {
586 
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") << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start;
631  if (endIsFound)
632  edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going left";
633  }
634 
635  // save recHits and add energy
636  if (!endIsFound) {
637  gapDetIds.push_back(next);
638  if (RecHits->find(next) != RecHits->end()) {
639  EcalRecHit nexthit = *RecHits->find(next);
640  ++gapCellCounter;
641  gapRecHits.push_back(nexthit);
642  gapEnergy += nexthit.energy();
643  cellGeom = subGeom->getGeometry(next);
644  eta = cellGeom->getPosition().eta();
645  gapET += nexthit.energy() / cosh(eta);
646  }
647  }
648 
649  // turn right to follow gap
650  currDirection = turnRight(currDirection, reverseOrientation);
651 
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 
666  while (!endIsFound) {
667 
668  bool nextStepFound = false;
669  int status = 0;
670  noDirs = 0;
671  while (!nextStepFound) {
672  next = makeStepInDirection(currDirection, theEcalNav.get());
673  theEcalNav->setHome(current);
674  theEcalNav->home();
675  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
676  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
677  if (status > 0) {
678  // Find dead cell along border -> end of cluster
679  endIsFound = true;
680  break;
681  } else if (next == EcalDetId(0)) {
682  // In case the Ecal border -> go along gap
683  currDirection = turnRight(currDirection, reverseOrientation);
684  } else if (status == 0) {
685  if (RecHits->find(next) != RecHits->end()) {
686  nextStepFound = true;
687  } else {
688  endIsFound = true;
689  break;
690  }
691  }
692  ++noDirs;
693  if (noDirs > 4) {
694  throw cms::Exception("NoStartingDirection")
695  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
696  }
697  }
698 
699  // make next step
700  next = makeStepInDirection(currDirection, theEcalNav.get());
701  current = next;
702 
703  if (debug) {
704  edm::LogInfo("EcalBoundaryInfoCalculator")<< "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start;
705  if (endIsFound)
706  edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going right";
707  }
708 
709  // save recHits and add energy
710  if (!endIsFound) {
711  gapDetIds.push_back(next);
712  if (RecHits->find(next) != RecHits->end()) {
713  EcalRecHit nexthit = *RecHits->find(next);
714  ++gapCellCounter;
715  gapRecHits.push_back(nexthit);
716  gapEnergy += nexthit.energy();
717  cellGeom = subGeom->getGeometry(next);
718  eta = cellGeom->getPosition().eta();
719  gapET += nexthit.energy() / cosh(eta);
720  }
721  }
722 
723  // turn left to follow gap
724  currDirection = turnLeft(currDirection, reverseOrientation);
725 
726  }
727  }
728 
729  if (debug) {
730  edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
731  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of RecHits along gap: " << gapRecHits.size();
732  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of DetIds along gap: " << gapDetIds.size();
733  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap energy: " << gapEnergy;
734  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap ET: " << gapET;
735  }
736 
737  BoundaryInformation gapInfo;
738  gapInfo.subdet = hitdetid.subdet();
739  gapInfo.detIds = gapDetIds;
740  gapInfo.recHits = gapRecHits;
741  gapInfo.boundaryEnergy = gapEnergy;
742  gapInfo.boundaryET = gapET;
743  gapInfo.nextToBorder = nextToBorder;
744  std::vector<int> stati;
745  gapInfo.channelStatus = stati;
746 
747  return gapInfo;
748 }
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:49
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:515
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:247
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:275
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
Definition: EcalBoundaryInfoCalculator.h:210
command_line.start
start
Definition: command_line.py:166
BoundaryInformation::detIds
std::vector< DetId > detIds
Definition: BoundaryInformation.h:26
EcalRecHit::energy
float energy() const
Definition: EcalRecHit.h:68
BoundaryInformation::channelStatus
std::vector< int > channelStatus
Definition: BoundaryInformation.h:27
edm::SortedCollection::end
const_iterator end() const
Definition: SortedCollection.h:310
EcalBoundaryInfoCalculator::turnLeft
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:261
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:46
EcalCondObjectContainer< EcalChannelStatusCode >::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:17
BoundaryInformation::subdet
EcalSubdetector subdet
Definition: BoundaryInformation.h:30
BoundaryInformation::recHits
std::vector< EcalRecHit > recHits
Definition: BoundaryInformation.h:24
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:271
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:59
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:88
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::goBackOneCell ( CdOrientation  currDirection,
EcalDetId  prev,
CaloNavigator< EcalDetId > *  theEcalNav 
) const
inlineprivate

Definition at line 235 of file EcalBoundaryInfoCalculator.h.

235  {
236  std::map<CdOrientation, CdOrientation>::iterator oIt = oppositeDirs.find(currDirection);
237  CdOrientation oppDirection=none;
238  if (oIt != oppositeDirs.end()) {
239  oppDirection = oIt->second;
240  theEcalNav->setHome(prev);
241  }
242  EcalDetId currDetId = theEcalNav->pos();
243 
244  return oppDirection;
245  }
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:106
CaloNavigator::pos
T pos() const
get the current position
Definition: CaloNavigator.h:32
EcalBoundaryInfoCalculator::oppositeDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:292
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 275 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 210 of file EcalBoundaryInfoCalculator.h.

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

210  {
211  EcalDetId next;
212  switch (direction) {
213  case north: {
214  next = theNavi->north();
215  break;
216  }
217  case east: {
218  next = theNavi->east();
219  break;
220  }
221  case south: {
222  next = theNavi->south();
223  break;
224  }
225  case west: {
226  next = theNavi->west();
227  break;
228  }
229  default:
230  break;
231  }
232  return next;
233  }
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:31
CaloNavigator::west
T west() const
move the navigator west
Definition: CaloNavigator.h:59
CaloNavigator::south
T south() const
move the navigator south
Definition: CaloNavigator.h:45
CaloNavigator::east
T east() const
move the navigator east
Definition: CaloNavigator.h:52
CaloNavigator::north
T north() const
move the navigator north
Definition: CaloNavigator.h:38
template<class EcalDetId>
void EcalBoundaryInfoCalculator< EcalDetId >::setDebugMode ( )
inline

Definition at line 203 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter().

203  {
204  edm::LogInfo("EcalBoundaryInfoCalculator") << "set Debug Mode!";
205  debug = true;
206  }
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnLeft ( CdOrientation  currDirection,
bool  reverseOrientation 
) const
inlineprivate

Definition at line 261 of file EcalBoundaryInfoCalculator.h.

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

261  {
262  //read nextDirection
263  std::map<CdOrientation, CdOrientation> turnMap = prevDirs;
264  if (reverseOrientation)
265  turnMap = nextDirs;
266  std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
267  CdOrientation nextDirection=none;
268  if (nIt != turnMap.end())
269  nextDirection = (*nIt).second;
270  else
271  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
272  return nextDirection;
273  }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:290
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:291
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnRight ( CdOrientation  currDirection,
bool  reverseOrientation 
) const
inlineprivate

Definition at line 247 of file EcalBoundaryInfoCalculator.h.

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

247  {
248  //read nextDirection
249  std::map<CdOrientation, CdOrientation> turnMap = nextDirs;
250  if (reverseOrientation)
251  turnMap = prevDirs;
252  std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
253  CdOrientation nextDirection=none;
254  if (nIt != turnMap.end())
255  nextDirection = (*nIt).second;
256  else
257  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
258  return nextDirection;
259  }
EcalBoundaryInfoCalculator::nextDirs
std::map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:290
EcalBoundaryInfoCalculator::prevDirs
std::map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:291

Member Data Documentation

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

Definition at line 293 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 290 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 292 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 291 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