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

Constructor & Destructor Documentation

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

Definition at line 294 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 318 of file EcalBoundaryInfoCalculator.h.

318  {
319 }

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 321 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, PV3DBase< T, PVType, FrameType >::eta(), Exception, EcalCondObjectContainer< T >::find(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), 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 
325  //initialize boundary information
326  std::vector<EcalRecHit> boundaryRecHits;
327  std::vector<DetId> boundaryDetIds;
328  std::vector<int> stati;
329 
330  double boundaryEnergy = 0;
331  double boundaryET = 0;
332  int beCellCounter = 0;
333  bool nextToBorder = false;
334 
335  boundaryRecHits.push_back(*hit);
336  ++beCellCounter;
337  boundaryEnergy += hit->energy();
338  EcalDetId hitdetid = (EcalDetId) hit->id();
339  boundaryDetIds.push_back(hitdetid);
340  const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(hitdetid);
341  const CaloCellGeometry* cellGeom = subGeom->getGeometry(hitdetid);
342  double eta = cellGeom->getPosition().eta();
343  boundaryET += hit->energy() / cosh(eta);
344 
345  if (debug) {
346  edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Boundary RecHits...";
347 
348  if (hitdetid.subdet() == EcalBarrel) {
349  edm::LogInfo("EcalBoundaryInfoCalculator") << "Starting at : (" << ((EBDetId) hitdetid).ieta() << "," << ((EBDetId) hitdetid).iphi() << ")";
350  } else if (hitdetid.subdet() == EcalEndcap) {
351  edm::LogInfo("EcalBoundaryInfoCalculator") << "Starting at : (" << ((EEDetId) hitdetid).ix() << "," << ((EEDetId) hitdetid).iy() << ")";
352  }
353  }
354 
355  //initialize navigator
356  auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
357  CdOrientation currDirection = north;
358  bool reverseOrientation = false;
359 
360  EcalDetId next(0);
361  EcalDetId start = hitdetid;
362  EcalDetId current = start;
363  int current_status = 0;
364 
365  // Search until a dead cell is ahead
366  bool startAlgo = false;
367  int noDirs = 0;
368  while (!startAlgo) {
369  next = makeStepInDirection(currDirection, theEcalNav.get());
370  theEcalNav->setHome(current);
371  theEcalNav->home();
372  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
373  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
374  if (status > 0) {
375  stati.push_back(status);
376  startAlgo = true;
377  break;
378  }
379  currDirection = turnLeft(currDirection, reverseOrientation);
380  ++noDirs;
381  if (noDirs > 4) {
382  throw cms::Exception("NoStartingDirection")
383  << "EcalBoundaryInfoCalculator: No starting direction can be found: This should never happen if RecHit has a dead neighbour!";
384  }
385  }
386 
387  // go around dead clusters counter clock wise
388  currDirection = turnRight(currDirection, reverseOrientation);
389 
390  // Search for next boundary element
391  bool nextIsStart = false;
392  bool atBorder = false;
393 
394  while (!nextIsStart) {
395 
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") << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start;
447 
448  // save recHits and add energy if on the boundary (and not inside at border)
449  if ((!atBorder || status == 0) && !nextIsStart) {
450  boundaryDetIds.push_back(next);
451  if (RecHits->find(next) != RecHits->end() && status == 0) {
452  EcalRecHit nexthit = *RecHits->find(next);
453  ++beCellCounter;
454  boundaryRecHits.push_back(nexthit);
455  boundaryEnergy += nexthit.energy();
456  cellGeom = subGeom->getGeometry(hitdetid);
457  eta = cellGeom->getPosition().eta();
458  boundaryET += nexthit.energy() / cosh(eta);
459  }
460  }
461 
462  if (current_status == 0 && status == 0 && atBorder) {
463  // this is for a special case, where dead cells are at border corner
464  currDirection = turnRight(currDirection, reverseOrientation);
465  } else {
466  // if dead region along a border is left, turn left
467  if (status == 0 && atBorder) {
468  atBorder = false;
469  currDirection = turnLeft(currDirection, reverseOrientation);
470  }
471  if (status == 0) {
472  // if outside the cluster turn left to follow boundary
473  currDirection = turnLeft(currDirection, reverseOrientation);
474  } else {
475  // else turn right to check if dead region can be left
476  currDirection = turnRight(currDirection, reverseOrientation);
477  }
478  }
479 
480  // save currect position
481  current = next;
482  current_status = status;
483 
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:45
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:244
CaloCellGeometry::getPosition
virtual const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.h:84
mps_update.status
status
Definition: mps_update.py:57
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
EcalBoundaryInfoCalculator::initializeEcalNavigator
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator(DetId startE, const edm::ESHandle< CaloTopology > theCaloTopology, EcalSubdetector ecalSubDet) const
Definition: EcalBoundaryInfoCalculator.h:272
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
Definition: EcalBoundaryInfoCalculator.h:207
CaloSubdetectorGeometry::getGeometry
virtual 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:44
command_line.start
start
Definition: command_line.py:165
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:321
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:258
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:77
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
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:271
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:60
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:90
template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasDeadNeighbour ( const EcalRecHit &  hit,
const edm::ESHandle< EcalChannelStatus >  ecalStatus,
std::vector< int > &  stati 
) const
inline

Definition at line 40 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

41  {
42 
43  stati.clear();
44  EcalDetId hitdetid = EcalDetId(hit.id());
45 
46  if (hitdetid.subdet() == EcalBarrel) {
47 
48  EBDetId ebhitdetid = (EBDetId) hitdetid;
49 
50  int hitIeta = ebhitdetid.ieta();
51  int hitIphi = ebhitdetid.iphi();
52 
53  for (int ieta = -1; ieta <= 1; ieta++) {
54  for (int iphi = -1; iphi <= 1; iphi++) {
55  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
56  //if (iphi == 0 && ieta == 0)
57  continue;
58  int neighbourIeta = hitIeta + ieta;
59  int neighbourIphi = hitIphi + iphi;
60  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
61  if (neighbourIphi < 1)
62  neighbourIphi += 360;
63  if (neighbourIphi > 360)
64  neighbourIphi -= 360;
65  if (neighbourIeta == 0) {
66  neighbourIeta += ieta;
67  }
68  }
69 
70  if (EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
71 
72  const EBDetId detid = EBDetId(neighbourIeta, neighbourIphi, EBDetId::ETAPHIMODE);
73  EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
74  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
75 
76  if (status > 0) {
77  bool present = false;
78  for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
79  if (*s == status) {
80  present = true;
81  break;
82  }
83  }
84  if (!present)
85  stati.push_back(status);
86  }
87  }
88  }
89  }
90 
91  } else if (hitdetid.subdet() == EcalEndcap) {
92 
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 
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.size() > 0)
132  return true;
133  return false;
134 
135  }
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:57
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:60
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:90
template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasInvalidNeighbour ( const EcalRecHit &  hit,
const edm::ESHandle< EcalChannelStatus >  ecalStatus 
) const
inline

Definition at line 137 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

137  {
139 
140  EcalDetId hitdetid = EcalDetId(hit.id());
141 
142  if (hitdetid.subdet() == EcalBarrel) {
143 
144  EBDetId ebhitdetid = (EBDetId) hitdetid;
145 
146  int hitIeta = ebhitdetid.ieta();
147  int hitIphi = ebhitdetid.iphi();
148 
149  for (int ieta = -1; ieta <= 1; ieta++) {
150  for (int iphi = -1; iphi <= 1; iphi++) {
151  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
152  //if (iphi == 0 && ieta == 0)
153  continue;
154  int neighbourIeta = hitIeta + ieta;
155  int neighbourIphi = hitIphi + iphi;
156  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
157  if (neighbourIphi < 1)
158  neighbourIphi += 360;
159  if (neighbourIphi > 360)
160  neighbourIphi -= 360;
161  if (neighbourIeta == 0) {
162  neighbourIeta += ieta;
163  }
164  }
165 
166  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
167  return true;
168  }
169  }
170  }
171 
172  } else if (hitdetid.subdet() == EcalEndcap) {
173 
174  EEDetId eehitdetid = (EEDetId) hitdetid;
175  int hitIx = eehitdetid.ix();
176  int hitIy = eehitdetid.iy();
177  int hitIz = eehitdetid.zside();
178 
179  for (int ix = -1; ix <= 1; ix++) {
180  for (int iy = -1; iy <= 1; iy++) {
181  if ((ix == 0 && iy == 0) || ix * iy != 0)
182  //if (ix == 0 && iy == 0)
183  continue;
184  int neighbourIx = hitIx + ix;
185  int neighbourIy = hitIy + iy;
186 
187  if (!EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
188  return true;
189  }
190  }
191  }
192 
193  } else {
194  edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
195  }
196 
197  return false;
198  }
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 512 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, PV3DBase< T, PVType, FrameType >::eta(), Exception, EcalCondObjectContainer< T >::find(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), 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().

514  {
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  const CaloCellGeometry* 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") << "Starting at : (" << ((EBDetId) hitdetid).ieta() << "," << ((EBDetId) hitdetid).iphi() << ")";
540  } else if (hitdetid.subdet() == EcalEndcap) {
541  edm::LogInfo("EcalBoundaryInfoCalculator") << "Starting at : (" << ((EEDetId) hitdetid).ix() << "," << ((EEDetId) hitdetid).iy() << ")";
542  }
543  }
544 
545  //initialize navigator
546  auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
547  CdOrientation currDirection = north;
548  bool reverseOrientation = false;
549 
550  EcalDetId next(0);
551  EcalDetId start = hitdetid;
552  EcalDetId current = start;
553 
554  // Search until a invalid cell is ahead
555  bool startAlgo = false;
556  int noDirs = 0;
557  while (!startAlgo) {
558  next = makeStepInDirection(currDirection, theEcalNav.get());
559  theEcalNav->setHome(start);
560  theEcalNav->home();
561  if (next == EcalDetId(0)) {
562  startAlgo = true;
563  nextToBorder = true;
564  break;
565  }
566  currDirection = turnLeft(currDirection, reverseOrientation);
567  ++noDirs;
568  if (noDirs > 4) {
569  throw cms::Exception("NoStartingDirection")
570  << "EcalBoundaryInfoCalculator: No starting direction can be found: This should never happen if RecHit is at border!";
571  }
572  }
573 
575  CdOrientation startDirection = currDirection;
576  currDirection = turnLeft(currDirection, reverseOrientation);
577 
578  // Search for next border element
579  bool endIsFound = false;
580  bool startIsEnd = false;
581 
582  while (!endIsFound) {
583 
584  bool nextStepFound = false;
585  int status = 0;
586  noDirs = 0;
587  while (!nextStepFound) {
588  next = makeStepInDirection(currDirection, theEcalNav.get());
589  theEcalNav->setHome(current);
590  theEcalNav->home();
591  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
592  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
593  if (status > 0) {
594  // Find dead cell along border -> end of cluster
595  endIsFound = true;
596  break;
597  } else if (next == EcalDetId(0)) {
598  // In case the Ecal border -> go along gap
599  currDirection = turnLeft(currDirection, reverseOrientation);
600  } else if (status == 0) {
601  if (RecHits->find(next) != RecHits->end()) {
602  nextStepFound = true;
603  } else {
604  endIsFound = true;
605  break;
606  }
607  }
608  ++noDirs;
609  if (noDirs > 4) {
610  throw cms::Exception("NoNextDirection")
611  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
612  }
613  }
614 
615  // make next step
616  next = makeStepInDirection(currDirection, theEcalNav.get());
617  current = next;
618 
619  if (next == start) {
620  startIsEnd = true;
621  endIsFound = true;
622  if (debug)
623  edm::LogInfo("EcalBoundaryInfoCalculator") << "Path along gap reached starting position!";
624  }
625 
626  if (debug) {
627  edm::LogInfo("EcalBoundaryInfoCalculator") << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start;
628  if (endIsFound)
629  edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going left";
630  }
631 
632  // save recHits and add energy
633  if (!endIsFound) {
634  gapDetIds.push_back(next);
635  if (RecHits->find(next) != RecHits->end()) {
636  EcalRecHit nexthit = *RecHits->find(next);
637  ++gapCellCounter;
638  gapRecHits.push_back(nexthit);
639  gapEnergy += nexthit.energy();
640  cellGeom = subGeom->getGeometry(next);
641  eta = cellGeom->getPosition().eta();
642  gapET += nexthit.energy() / cosh(eta);
643  }
644  }
645 
646  // turn right to follow gap
647  currDirection = turnRight(currDirection, reverseOrientation);
648 
649  }
650 
652  theEcalNav->setHome(start);
653  theEcalNav->home();
654  current = start;
655  currDirection = startDirection;
656  currDirection = turnRight(currDirection, reverseOrientation);
657 
658  // Search for next border element
659  endIsFound = false;
660 
661  if (!startIsEnd) {
662 
663  while (!endIsFound) {
664 
665  bool nextStepFound = false;
666  int status = 0;
667  noDirs = 0;
668  while (!nextStepFound) {
669  next = makeStepInDirection(currDirection, theEcalNav.get());
670  theEcalNav->setHome(current);
671  theEcalNav->home();
672  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
673  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
674  if (status > 0) {
675  // Find dead cell along border -> end of cluster
676  endIsFound = true;
677  break;
678  } else if (next == EcalDetId(0)) {
679  // In case the Ecal border -> go along gap
680  currDirection = turnRight(currDirection, reverseOrientation);
681  } else if (status == 0) {
682  if (RecHits->find(next) != RecHits->end()) {
683  nextStepFound = true;
684  } else {
685  endIsFound = true;
686  break;
687  }
688  }
689  ++noDirs;
690  if (noDirs > 4) {
691  throw cms::Exception("NoStartingDirection")
692  << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
693  }
694  }
695 
696  // make next step
697  next = makeStepInDirection(currDirection, theEcalNav.get());
698  current = next;
699 
700  if (debug) {
701  edm::LogInfo("EcalBoundaryInfoCalculator")<< "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 
726  if (debug) {
727  edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
728  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of RecHits along gap: " << gapRecHits.size();
729  edm::LogInfo("EcalBoundaryInfoCalculator") << "No of DetIds along gap: " << gapDetIds.size();
730  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap energy: " << gapEnergy;
731  edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap ET: " << gapET;
732  }
733 
734  BoundaryInformation gapInfo;
735  gapInfo.subdet = hitdetid.subdet();
736  gapInfo.detIds = gapDetIds;
737  gapInfo.recHits = gapRecHits;
738  gapInfo.boundaryEnergy = gapEnergy;
739  gapInfo.boundaryET = gapET;
740  gapInfo.nextToBorder = nextToBorder;
741  std::vector<int> stati;
742  gapInfo.channelStatus = stati;
743 
744  return gapInfo;
745 }
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:45
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:512
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
Definition: EcalBoundaryInfoCalculator.h:244
CaloCellGeometry::getPosition
virtual const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.h:84
mps_update.status
status
Definition: mps_update.py:57
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
EcalBoundaryInfoCalculator::initializeEcalNavigator
std::unique_ptr< CaloNavigator< EcalDetId > > initializeEcalNavigator(DetId startE, const edm::ESHandle< CaloTopology > theCaloTopology, EcalSubdetector ecalSubDet) const
Definition: EcalBoundaryInfoCalculator.h:272
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator< EcalDetId > *theNavi) const
Definition: EcalBoundaryInfoCalculator.h:207
CaloSubdetectorGeometry::getGeometry
virtual 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:44
command_line.start
start
Definition: command_line.py:165
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:258
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:77
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
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:76
edm::SortedCollection::find
iterator find(key_type k)
Definition: SortedCollection.h:271
EcalCondObjectContainer::find
const_iterator find(uint32_t rawId) const
Definition: EcalCondObjectContainer.h:60
EcalCondObjectContainer::end
const_iterator end() const
Definition: EcalCondObjectContainer.h:90
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::goBackOneCell ( CdOrientation  currDirection,
EcalDetId  prev,
CaloNavigator< EcalDetId > *  theEcalNav 
) const
inlineprivate

Definition at line 232 of file EcalBoundaryInfoCalculator.h.

232  {
233  std::map<CdOrientation, CdOrientation>::iterator oIt = oppositeDirs.find(currDirection);
234  CdOrientation oppDirection=none;
235  if (oIt != oppositeDirs.end()) {
236  oppDirection = oIt->second;
237  theEcalNav->setHome(prev);
238  }
239  EcalDetId currDetId = theEcalNav->pos();
240 
241  return oppDirection;
242  }
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:289
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 272 of file EcalBoundaryInfoCalculator.h.

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

273  {
274  std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav(nullptr);
275  if (ecalSubDet == EcalBarrel) {
276  theEcalNav.reset(new CaloNavigator<EcalDetId> ((EBDetId) startE, (theCaloTopology->getSubdetectorTopology(
277  DetId::Ecal, ecalSubDet))));
278  } else if (ecalSubDet == EcalEndcap) {
279  theEcalNav.reset(new CaloNavigator<EcalDetId> ((EEDetId) startE, (theCaloTopology->getSubdetectorTopology(
280  DetId::Ecal, ecalSubDet))));
281  } else {
282  edm::LogWarning("EcalBoundaryInfoCalculator") << "initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
283  }
284  return theEcalNav;
285  }
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 207 of file EcalBoundaryInfoCalculator.h.

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

207  {
208  EcalDetId next;
209  switch (direction) {
210  case north: {
211  next = theNavi->north();
212  break;
213  }
214  case east: {
215  next = theNavi->east();
216  break;
217  }
218  case south: {
219  next = theNavi->south();
220  break;
221  }
222  case west: {
223  next = theNavi->west();
224  break;
225  }
226  default:
227  break;
228  }
229  return next;
230  }
GetRecoTauVFromDQM_MC_cff.next
next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
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 200 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::EcalDeadCellBoundaryEnergyFilter().

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

Definition at line 258 of file EcalBoundaryInfoCalculator.h.

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

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

Definition at line 244 of file EcalBoundaryInfoCalculator.h.

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

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

Member Data Documentation

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

Definition at line 290 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 287 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 289 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 288 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