CMS 3D CMS Logo
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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)
 
bool checkRecHitHasDeadNeighbour (const EcalRecHit hit, const edm::ESHandle< EcalChannelStatus > ecalStatus, vector< int > &stati)
 
bool checkRecHitHasInvalidNeighbour (const EcalRecHit hit, const edm::ESHandle< EcalChannelStatus > ecalStatus)
 
 EcalBoundaryInfoCalculator ()
 
BoundaryInformation gapRecHits (const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry)
 
void setDebugMode ()
 
 ~EcalBoundaryInfoCalculator ()
 

Private Member Functions

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

Private Attributes

bool debug
 
map< CdOrientation, CdOrientation > nextDirs
 
map< CdOrientation, CdOrientation > oppositeDirs
 
map< CdOrientation, CdOrientation > prevDirs
 
CaloNavigator< EcalDetId > * theEcalNav
 

Detailed Description

template<class EcalDetId>
class EcalBoundaryInfoCalculator< EcalDetId >

Definition at line 26 of file EcalBoundaryInfoCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 307 of file EcalBoundaryInfoCalculator.h.

References debug, east, north, south, and west.

307  {
308 
309  nextDirs.clear();
310  nextDirs[north] = east;
311  nextDirs[east] = south;
312  nextDirs[south] = west;
313  nextDirs[west] = north;
314 
315  prevDirs.clear();
316  prevDirs[north] = west;
317  prevDirs[east] = north;
318  prevDirs[south] = east;
319  prevDirs[west] = south;
320 
321  oppositeDirs.clear();
322  oppositeDirs[north] = south;
323  oppositeDirs[south] = north;
324  oppositeDirs[east] = west;
325  oppositeDirs[west] = east;
326 
327  theEcalNav = 0;
328  debug = false;
329 
330 }
EcalBoundaryInfoCalculator::nextDirs
map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:299
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302
EcalBoundaryInfoCalculator::prevDirs
map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:300
EcalBoundaryInfoCalculator::oppositeDirs
map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:301
template<class EcalDetId >
EcalBoundaryInfoCalculator< EcalDetId >::~EcalBoundaryInfoCalculator ( )

Definition at line 332 of file EcalBoundaryInfoCalculator.h.

332  {
333  delete theEcalNav;
334 }
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302

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 
)

Definition at line 336 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, gather_cfg::cout, cond::rpcobimon::current, debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, CaloRecHit::energy(), PV3DBase< T, PVType, FrameType >::eta(), eta(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), EcalRecHit::id(), GetRecoTauVFromDQM_MC_cff::next, BoundaryInformation::nextToBorder, north, BoundaryInformation::recHits, alignCSCRings::s, errorMatrix2Lands_multiChannel::start, ntuplemaker::status, and BoundaryInformation::subdet.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

338  {
339 
340  //initialize boundary information
341  std::vector<EcalRecHit> boundaryRecHits;
342  std::vector<DetId> boundaryDetIds;
343  std::vector<int> stati;
344 
345  double boundaryEnergy = 0;
346  double boundaryET = 0;
347  int beCellCounter = 0;
348  bool nextToBorder = false;
349 
350  boundaryRecHits.push_back(*hit);
351  ++beCellCounter;
352  boundaryEnergy += hit->energy();
353  EcalDetId hitdetid = (EcalDetId) hit->id();
354  boundaryDetIds.push_back(hitdetid);
355  const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(hitdetid);
356  const CaloCellGeometry* cellGeom = subGeom->getGeometry(hitdetid);
357  double eta = cellGeom->getPosition().eta();
358  boundaryET += hit->energy() / cosh(eta);
359 
360  if (debug) {
361  std::cout << "Find Boundary RecHits..." << std::endl;
362 
363  if (hitdetid.subdet() == EcalBarrel) {
364  std::cout << "Starting at : (" << ((EBDetId) hitdetid).ieta() << "," << ((EBDetId) hitdetid).iphi() << ")"
365  << std::endl;
366 
367  } else if (hitdetid.subdet() == EcalEndcap) {
368  std::cout << "Starting at : (" << ((EEDetId) hitdetid).ix() << "," << ((EEDetId) hitdetid).iy() << ")"
369  << std::endl;
370  }
371  }
372 
373  //initialize navigator
374  initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
375  CdOrientation currDirection = north;
376  bool reverseOrientation = false;
377 
378  EcalDetId next(0);
379  EcalDetId start = hitdetid;
380  EcalDetId current = start;
381  int current_status = 0;
382 
383  // Search until a dead cell is ahead
384  bool startAlgo = false;
385  int noDirs = 0;
386  while (!startAlgo) {
387  next = makeStepInDirection(currDirection, theEcalNav);
388  theEcalNav->setHome(current);
389  theEcalNav->home();
390  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
391  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
392  if (status > 0) {
393  stati.push_back(status);
394  startAlgo = true;
395  break;
396  }
397  currDirection = turnLeft(currDirection, reverseOrientation);
398  ++noDirs;
399  if (noDirs > 4) {
400 
401  cout << "No starting direction can be found: This should never happen if RecHit has a dead neighbour!" << endl;
402  throw "ERROR";
403  break;
404  }
405  }
406 
407  // go around dead clusters counter clock wise
408  currDirection = turnRight(currDirection, reverseOrientation);
409 
410  // Search for next boundary element
411  bool nextIsStart = false;
412  bool atBorder = false;
413 
414  while (!nextIsStart) {
415 
416  bool nextStepFound = false;
417  int status = 0;
418  noDirs = 0;
419  while (!nextStepFound) {
420  next = makeStepInDirection(currDirection, theEcalNav);
421  theEcalNav->setHome(current);
422  theEcalNav->home();
423  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
424  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
425  if (status > 0) {
426  // New dead cell found: update status vector of dead channels
427  bool present = false;
428  for (vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
429  if (*s == status) {
430  present = true;
431  break;
432  }
433  }
434  if (!present)
435  stati.push_back(status);
436 
437  if (atBorder) {
438  nextStepFound = true;
439  } else {
440  currDirection = turnRight(currDirection, reverseOrientation);
441  }
442  } else if (next == EcalDetId(0)) {
443  // In case the Ecal border is reached -> go along dead cells
444  currDirection = turnLeft(currDirection, reverseOrientation);
445  atBorder = true;
446  } else if (status == 0) {
447  nextStepFound = true;
448  }
449  ++noDirs;
450  if (noDirs > 4) {
451  cout << "No valid next direction can be found: This should never happen!" << endl;
452  throw "ERROR";
453  break;
454  }
455  }
456 
457  // make next step
458  next = makeStepInDirection(currDirection, theEcalNav);
459 
460  if (next == start) {
461  nextIsStart = true;
462  if (debug)
463  cout << "Boundary path reached starting position!" << endl;
464  }
465 
466  if (debug)
467  cout << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start << endl;
468 
469  // save recHits and add energy if on the boundary (and not inside at border)
470  if ((!atBorder || status == 0) && !nextIsStart) {
471  boundaryDetIds.push_back(next);
472  if (RecHits->find(next) != RecHits->end() && status == 0) {
473  EcalRecHit nexthit = *RecHits->find(next);
474  ++beCellCounter;
475  boundaryRecHits.push_back(nexthit);
476  boundaryEnergy += nexthit.energy();
477  cellGeom = subGeom->getGeometry(hitdetid);
478  eta = cellGeom->getPosition().eta();
479  boundaryET += nexthit.energy() / cosh(eta);
480  }
481  }
482 
483  if (current_status == 0 && status == 0 && atBorder) {
484  // this is for a special case, where dead cells are at border corner
485  currDirection = turnRight(currDirection, reverseOrientation);
486  } else {
487  // if dead region along a border is left, turn left
488  if (status == 0 && atBorder) {
489  atBorder = false;
490  currDirection = turnLeft(currDirection, reverseOrientation);
491  }
492  if (status == 0) {
493  // if outside the cluster turn left to follow boundary
494  currDirection = turnLeft(currDirection, reverseOrientation);
495  } else {
496  // else turn right to check if dead region can be left
497  currDirection = turnRight(currDirection, reverseOrientation);
498  }
499  }
500 
501  // save currect position
502  current = next;
503  current_status = status;
504 
505  }
506 
507  if (debug) {
508  cout << "<<<<<<<<<<<<<<< Final Boundary object <<<<<<<<<<<<<<<" << endl;
509  cout << "no of neighbouring RecHits: " << boundaryRecHits.size() << endl;
510  cout << "no of neighbouring DetIds: " << boundaryDetIds.size() << endl;
511  cout << "boundary energy: " << boundaryEnergy << endl;
512  cout << "boundary ET: " << boundaryET << endl;
513  cout << "no of cells contributing to boundary energy: " << beCellCounter << endl;
514  cout << "Channel stati: ";
515  for (vector<int>::iterator it = stati.begin(); it != stati.end(); ++it) {
516  cout << *it << " ";
517  }
518  cout << endl;
519  }
520 
521  BoundaryInformation boundInfo;
522  boundInfo.subdet = hitdetid.subdet();
523  boundInfo.detIds = boundaryDetIds;
524  boundInfo.recHits = boundaryRecHits;
525  boundInfo.boundaryEnergy = boundaryEnergy;
526  boundInfo.boundaryET = boundaryET;
527  boundInfo.nextToBorder = nextToBorder;
528  boundInfo.channelStatus = stati;
529 
530  if (theEcalNav != 0) {
531  delete theEcalNav;
532  theEcalNav = 0;
533  }
534  return boundInfo;
535 }
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:136
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation)
Definition: EcalBoundaryInfoCalculator.h:249
CaloNavigator::home
void home() const
move the navigator back to the starting point
Definition: CaloNavigator.h:143
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)
Definition: EcalBoundaryInfoCalculator.h:336
EcalBoundaryInfoCalculator::initializeEcalNavigator
void initializeEcalNavigator(DetId startE, const edm::ESHandle< CaloTopology > theCaloTopology, EcalSubdetector ecalSubDet)
Definition: EcalBoundaryInfoCalculator.h:277
eta
T eta() const
Definition: Basic3DVectorLD.h:172
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.cc:46
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:43
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302
CaloRecHit::energy
float energy() const
Definition: CaloRecHit.h:19
BoundaryInformation::detIds
std::vector< DetId > detIds
Definition: BoundaryInformation.h:24
EcalBoundaryInfoCalculator::turnLeft
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation)
Definition: EcalBoundaryInfoCalculator.h:263
BoundaryInformation::channelStatus
std::vector< int > channelStatus
Definition: BoundaryInformation.h:25
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, CaloNavigator< EcalDetId > *theNavi)
Definition: EcalBoundaryInfoCalculator.h:208
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:76
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:15
BoundaryInformation::subdet
EcalSubdetector subdet
Definition: BoundaryInformation.h:28
BoundaryInformation::recHits
std::vector< EcalRecHit > recHits
Definition: BoundaryInformation.h:22
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:75
if
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GetRecoTauVFromDQM_MC_cff.next
tuple next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
errorMatrix2Lands_multiChannel.start
tuple start
Definition: errorMatrix2Lands_multiChannel.py:58
template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasDeadNeighbour ( const EcalRecHit  hit,
const edm::ESHandle< EcalChannelStatus >  ecalStatus,
vector< int > &  stati 
)
inline

Definition at line 41 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

42  {
43 
44  stati.clear();
45  EcalDetId hitdetid = EcalDetId(hit.id());
46 
47  if (hitdetid.subdet() == EcalBarrel) {
48 
49  EBDetId ebhitdetid = (EBDetId) hitdetid;
50 
51  int hitIeta = ebhitdetid.ieta();
52  int hitIphi = ebhitdetid.iphi();
53 
54  for (int ieta = -1; ieta <= 1; ieta++) {
55  for (int iphi = -1; iphi <= 1; iphi++) {
56  if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
57  //if (iphi == 0 && ieta == 0)
58  continue;
59  int neighbourIeta = hitIeta + ieta;
60  int neighbourIphi = hitIphi + iphi;
61  if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
62  if (neighbourIphi < 1)
63  neighbourIphi += 360;
64  if (neighbourIphi > 360)
65  neighbourIphi -= 360;
66  if (neighbourIeta == 0) {
67  neighbourIeta += ieta;
68  }
69  }
70 
71  if (EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
72 
73  const EBDetId detid = EBDetId(neighbourIeta, neighbourIphi, EBDetId::ETAPHIMODE);
74  EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
75  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
76 
77  if (status > 0) {
78  bool present = false;
79  for (vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
80  if (*s == status) {
81  present = true;
82  break;
83  }
84  }
85  if (!present)
86  stati.push_back(status);
87  }
88  }
89  }
90  }
91 
92  } else if (hitdetid.subdet() == EcalEndcap) {
93 
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 
109  const EEDetId detid = EEDetId(neighbourIx, neighbourIy, hitIz, EEDetId::XYMODE);
110  EcalChannelStatus::const_iterator chit = ecalStatus->find(detid);
111  int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
112 
113  if (status > 0) {
114  bool present = false;
115  for (vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
116  if (*s == status) {
117  present = true;
118  break;
119  }
120  }
121  if (!present)
122  stati.push_back(status);
123  }
124  }
125  }
126  }
127 
128  } else {
129  cout << "ERROR - RecHit belongs to wrong sub detector" << endl;
130  }
131 
132  if (stati.size() > 0)
133  return true;
134  return false;
135 
136  }
EBDetId::validDetId
static bool validDetId(int i, int j)
check if a valid index combination
Definition: EBDetId.cc:59
EEDetId::ix
int ix() const
Definition: EEDetId.h:71
EEDetId::XYMODE
static const int XYMODE
Definition: EEDetId.h:316
EEDetId::validDetId
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
Definition: EEDetId.cc:562
EBDetId::iphi
int iphi() const
get the crystal iphi
Definition: EBDetId.h:46
EEDetId::zside
int zside() const
Definition: EEDetId.h:65
EEDetId::iy
int iy() const
Definition: EEDetId.h:77
EBDetId::ieta
int ieta() const
get the crystal ieta
Definition: EBDetId.h:44
EBDetId::ETAPHIMODE
static const int ETAPHIMODE
Definition: EBDetId.h:145
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:76
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:15
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::checkRecHitHasInvalidNeighbour ( const EcalRecHit  hit,
const edm::ESHandle< EcalChannelStatus >  ecalStatus 
)
inline

Definition at line 138 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 537 of file EcalBoundaryInfoCalculator.h.

References BoundaryInformation::boundaryEnergy, BoundaryInformation::boundaryET, BoundaryInformation::channelStatus, gather_cfg::cout, cond::rpcobimon::current, debug, BoundaryInformation::detIds, EcalBarrel, EcalEndcap, CaloRecHit::energy(), PV3DBase< T, PVType, FrameType >::eta(), eta(), CaloSubdetectorGeometry::getGeometry(), CaloCellGeometry::getPosition(), EcalRecHit::id(), GetRecoTauVFromDQM_MC_cff::next, BoundaryInformation::nextToBorder, north, BoundaryInformation::recHits, errorMatrix2Lands_multiChannel::start, ntuplemaker::status, and BoundaryInformation::subdet.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

539  {
540 
541  //initialize boundary information
542  std::vector<EcalRecHit> gapRecHits;
543  std::vector<DetId> gapDetIds;
544 
545  double gapEnergy = 0;
546  double gapET = 0;
547  int gapCellCounter = 0;
548  bool nextToBorder = false;
549 
550  gapRecHits.push_back(*hit);
551  ++gapCellCounter;
552  gapEnergy += hit->energy();
553  EcalDetId hitdetid = (EcalDetId) hit->id();
554  gapDetIds.push_back(hitdetid);
555  const CaloSubdetectorGeometry* subGeom = geometry->getSubdetectorGeometry(hitdetid);
556  const CaloCellGeometry* cellGeom = subGeom->getGeometry(hitdetid);
557  double eta = cellGeom->getPosition().eta();
558  gapET += hit->energy() / cosh(eta);
559 
560  if (debug) {
561  std::cout << "Find Border RecHits..." << std::endl;
562 
563  if (hitdetid.subdet() == EcalBarrel) {
564  std::cout << "Starting at : (" << ((EBDetId) hitdetid).ieta() << "," << ((EBDetId) hitdetid).iphi() << ")"
565  << std::endl;
566 
567  } else if (hitdetid.subdet() == EcalEndcap) {
568  std::cout << "Starting at : (" << ((EEDetId) hitdetid).ix() << "," << ((EEDetId) hitdetid).iy() << ")"
569  << std::endl;
570  }
571  }
572 
573  //initialize navigator
574  initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
575  CdOrientation currDirection = north;
576  bool reverseOrientation = false;
577 
578  EcalDetId next(0);
579  EcalDetId start = hitdetid;
580  EcalDetId current = start;
581 
582  // Search until a invalid cell is ahead
583  bool startAlgo = false;
584  int noDirs = 0;
585  while (!startAlgo) {
586  next = makeStepInDirection(currDirection, theEcalNav);
587  theEcalNav->setHome(start);
588  theEcalNav->home();
589  if (next == EcalDetId(0)) {
590  startAlgo = true;
591  nextToBorder = true;
592  break;
593  }
594  currDirection = turnLeft(currDirection, reverseOrientation);
595  ++noDirs;
596  if (noDirs > 4) {
597 
598  cout << "No starting direction can be found: This should never happen if RecHit is at border!" << endl;
599  throw "ERROR";
600  break;
601  }
602  }
603 
605  CdOrientation startDirection = currDirection;
606  currDirection = turnLeft(currDirection, reverseOrientation);
607 
608  // Search for next border element
609  bool endIsFound = false;
610  bool startIsEnd = false;
611 
612  while (!endIsFound) {
613 
614  bool nextStepFound = false;
615  int status = 0;
616  noDirs = 0;
617  while (!nextStepFound) {
618  next = makeStepInDirection(currDirection, theEcalNav);
619  theEcalNav->setHome(current);
620  theEcalNav->home();
621  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
622  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
623  if (status > 0) {
624  // Find dead cell along border -> end of cluster
625  endIsFound = true;
626  break;
627  } else if (next == EcalDetId(0)) {
628  // In case the Ecal border -> go along gap
629  currDirection = turnLeft(currDirection, reverseOrientation);
630  } else if (status == 0) {
631  if (RecHits->find(next) != RecHits->end()) {
632  nextStepFound = true;
633  } else {
634  endIsFound = true;
635  break;
636  }
637  }
638  ++noDirs;
639  if (noDirs > 4) {
640  cout << "No valid next direction can be found: This should never happen!" << endl;
641  throw "ERROR";
642  break;
643  }
644  }
645 
646  // make next step
647  next = makeStepInDirection(currDirection, theEcalNav);
648  current = next;
649 
650  if (next == start) {
651  startIsEnd = true;
652  endIsFound = true;
653  if (debug)
654  cout << "Path along gap reached starting position!" << endl;
655  }
656 
657  if (debug) {
658  cout << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start << endl;
659  if (endIsFound)
660  cout << "End of gap cluster is found going left" << endl;
661  }
662 
663  // save recHits and add energy
664  if (!endIsFound) {
665  gapDetIds.push_back(next);
666  if (RecHits->find(next) != RecHits->end()) {
667  EcalRecHit nexthit = *RecHits->find(next);
668  ++gapCellCounter;
669  gapRecHits.push_back(nexthit);
670  gapEnergy += nexthit.energy();
671  cellGeom = subGeom->getGeometry(next);
672  eta = cellGeom->getPosition().eta();
673  gapET += nexthit.energy() / cosh(eta);
674  }
675  }
676 
677  // turn right to follow gap
678  currDirection = turnRight(currDirection, reverseOrientation);
679 
680  }
681 
683  theEcalNav->setHome(start);
684  theEcalNav->home();
685  current = start;
686  currDirection = startDirection;
687  currDirection = turnRight(currDirection, reverseOrientation);
688 
689  // Search for next border element
690  endIsFound = false;
691 
692  if (!startIsEnd) {
693 
694  while (!endIsFound) {
695 
696  bool nextStepFound = false;
697  int status = 0;
698  noDirs = 0;
699  while (!nextStepFound) {
700  next = makeStepInDirection(currDirection, theEcalNav);
701  theEcalNav->setHome(current);
702  theEcalNav->home();
703  EcalChannelStatus::const_iterator chit = ecalStatus->find(next);
704  status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
705  if (status > 0) {
706  // Find dead cell along border -> end of cluster
707  endIsFound = true;
708  break;
709  } else if (next == EcalDetId(0)) {
710  // In case the Ecal border -> go along gap
711  currDirection = turnRight(currDirection, reverseOrientation);
712  } else if (status == 0) {
713  if (RecHits->find(next) != RecHits->end()) {
714  nextStepFound = true;
715  } else {
716  endIsFound = true;
717  break;
718  }
719  }
720  ++noDirs;
721  if (noDirs > 4) {
722  cout << "No valid next direction can be found: This should never happen!" << endl;
723  throw "ERROR";
724  break;
725  }
726  }
727 
728  // make next step
729  next = makeStepInDirection(currDirection, theEcalNav);
730  current = next;
731 
732  if (debug) {
733  cout << "Next step: " << (EcalDetId) next << " Status: " << status << " Start: " << (EcalDetId) start
734  << endl;
735  if (endIsFound)
736  cout << "End of gap cluster is found going right" << endl;
737  }
738 
739  // save recHits and add energy
740  if (!endIsFound) {
741  gapDetIds.push_back(next);
742  if (RecHits->find(next) != RecHits->end()) {
743  EcalRecHit nexthit = *RecHits->find(next);
744  ++gapCellCounter;
745  gapRecHits.push_back(nexthit);
746  gapEnergy += nexthit.energy();
747  cellGeom = subGeom->getGeometry(next);
748  eta = cellGeom->getPosition().eta();
749  gapET += nexthit.energy() / cosh(eta);
750  }
751  }
752 
753  // turn left to follow gap
754  currDirection = turnLeft(currDirection, reverseOrientation);
755 
756  }
757  }
758 
759  if (debug) {
760  cout << "<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<" << endl;
761  cout << "No of RecHits along gap: " << gapRecHits.size() << endl;
762  cout << "No of DetIds along gap: " << gapDetIds.size() << endl;
763  cout << "Gap energy: " << gapEnergy << endl;
764  cout << "Gap ET: " << gapET << endl;
765  }
766 
767  BoundaryInformation gapInfo;
768  gapInfo.subdet = hitdetid.subdet();
769  gapInfo.detIds = gapDetIds;
770  gapInfo.recHits = gapRecHits;
771  gapInfo.boundaryEnergy = gapEnergy;
772  gapInfo.boundaryET = gapET;
773  gapInfo.nextToBorder = nextToBorder;
774  vector<int> stati;
775  gapInfo.channelStatus = stati;
776 
777  if (theEcalNav != 0) {
778  delete theEcalNav;
779  theEcalNav = 0;
780  }
781  return gapInfo;
782 }
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:136
EcalBoundaryInfoCalculator::turnRight
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation)
Definition: EcalBoundaryInfoCalculator.h:249
CaloNavigator::home
void home() const
move the navigator back to the starting point
Definition: CaloNavigator.h:143
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)
Definition: EcalBoundaryInfoCalculator.h:537
EcalBoundaryInfoCalculator::initializeEcalNavigator
void initializeEcalNavigator(DetId startE, const edm::ESHandle< CaloTopology > theCaloTopology, EcalSubdetector ecalSubDet)
Definition: EcalBoundaryInfoCalculator.h:277
eta
T eta() const
Definition: Basic3DVectorLD.h:172
CaloCellGeometry::getPosition
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
Definition: CaloCellGeometry.cc:46
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:43
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302
CaloRecHit::energy
float energy() const
Definition: CaloRecHit.h:19
BoundaryInformation::detIds
std::vector< DetId > detIds
Definition: BoundaryInformation.h:24
EcalBoundaryInfoCalculator::turnLeft
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation)
Definition: EcalBoundaryInfoCalculator.h:263
BoundaryInformation::channelStatus
std::vector< int > channelStatus
Definition: BoundaryInformation.h:25
EcalBoundaryInfoCalculator::makeStepInDirection
EcalDetId makeStepInDirection(CdOrientation direction, CaloNavigator< EcalDetId > *theNavi)
Definition: EcalBoundaryInfoCalculator.h:208
EcalRecHit::id
DetId id() const
get the id
Definition: EcalRecHit.h:76
EcalCondObjectContainer::const_iterator
std::vector< Item >::const_iterator const_iterator
Definition: EcalCondObjectContainer.h:15
BoundaryInformation::subdet
EcalSubdetector subdet
Definition: BoundaryInformation.h:28
BoundaryInformation::recHits
std::vector< EcalRecHit > recHits
Definition: BoundaryInformation.h:22
PV3DBase::eta
T eta() const
Definition: PV3DBase.h:75
if
perl if(1 lt scalar(@::datatypes))
Definition: edlooper.cc:31
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
GetRecoTauVFromDQM_MC_cff.next
tuple next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
ntuplemaker.status
tuple status
Definition: ntuplemaker.py:245
errorMatrix2Lands_multiChannel.start
tuple start
Definition: errorMatrix2Lands_multiChannel.py:58
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::goBackOneCell ( CdOrientation  currDirection,
EcalDetId  prev 
)
inlineprivate

Definition at line 237 of file EcalBoundaryInfoCalculator.h.

237  {
238  map<CdOrientation, CdOrientation>::iterator oIt = oppositeDirs.find(currDirection);
239  CdOrientation oppDirection=none;
240  if (oIt != oppositeDirs.end()) {
241  oppDirection = oIt->second;
242  theEcalNav->setHome(prev);
243  }
244  EcalDetId currDetId = theEcalNav->pos();
245 
246  return oppDirection;
247  }
CaloNavigator::setHome
void setHome(const T &startingPoint)
set the starting position
Definition: CaloNavigator.h:136
CaloNavigator::pos
T pos() const
get the current position
Definition: CaloNavigator.h:45
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302
EcalBoundaryInfoCalculator::oppositeDirs
map< CdOrientation, CdOrientation > oppositeDirs
Definition: EcalBoundaryInfoCalculator.h:301
template<class EcalDetId>
void EcalBoundaryInfoCalculator< EcalDetId >::initializeEcalNavigator ( DetId  startE,
const edm::ESHandle< CaloTopology >  theCaloTopology,
EcalSubdetector  ecalSubDet 
)
inlineprivate

Definition at line 277 of file EcalBoundaryInfoCalculator.h.

278  {
279  if (ecalSubDet == EcalBarrel) {
280  if (theEcalNav != 0) {
281  delete theEcalNav;
282  theEcalNav = 0;
283  }
284  theEcalNav = new CaloNavigator<EcalDetId> ((EBDetId) startE, (theCaloTopology->getSubdetectorTopology(
285  DetId::Ecal, ecalSubDet)));
286  } else if (ecalSubDet == EcalEndcap) {
287  if (theEcalNav != 0) {
288  delete theEcalNav;
289  theEcalNav = 0;
290  }
291  theEcalNav = new CaloNavigator<EcalDetId> ((EEDetId) startE, (theCaloTopology->getSubdetectorTopology(
292  DetId::Ecal, ecalSubDet)));
293  } else {
294  cout << "initializeEcalNavigator not implemented for subDet: " << ecalSubDet << endl;
295  }
296 
297  }
EcalBoundaryInfoCalculator::theEcalNav
CaloNavigator< EcalDetId > * theEcalNav
Definition: EcalBoundaryInfoCalculator.h:302
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
template<class EcalDetId>
EcalDetId EcalBoundaryInfoCalculator< EcalDetId >::makeStepInDirection ( CdOrientation  direction,
CaloNavigator< EcalDetId > *  theNavi 
)
inlineprivate

Definition at line 208 of file EcalBoundaryInfoCalculator.h.

208  {
209  EcalDetId next;
210  switch (direction) {
211  case north: {
212  //cout<<"go north"<<endl;
213  next = theNavi->north();
214  break;
215  }
216  case east: {
217  //cout<<"go east"<<endl;
218  next = theNavi->east();
219  break;
220  }
221  case south: {
222  //cout<<"go south"<<endl;
223  next = theNavi->south();
224  break;
225  }
226  case west: {
227  //cout<<"go west"<<endl;
228  next = theNavi->west();
229  break;
230  }
231  default:
232  break;
233  }
234  return next;
235  }
CaloNavigator::west
virtual T west() const
move the navigator west
Definition: CaloNavigator.h:81
CaloNavigator::north
virtual T north() const
move the navigator north
Definition: CaloNavigator.h:51
CaloNavigator::east
virtual T east() const
move the navigator east
Definition: CaloNavigator.h:71
GetRecoTauVFromDQM_MC_cff.next
tuple next
Definition: GetRecoTauVFromDQM_MC_cff.py:30
CaloNavigator::south
virtual T south() const
move the navigator south
Definition: CaloNavigator.h:61
template<class EcalDetId>
void EcalBoundaryInfoCalculator< EcalDetId >::setDebugMode ( )
inline

Definition at line 201 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

201  {
202  cout << "set Debug Mode!" << endl;
203  debug = true;
204  }
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnLeft ( CdOrientation  currDirection,
bool  reverseOrientation 
)
inlineprivate

Definition at line 263 of file EcalBoundaryInfoCalculator.h.

263  {
264  //read nextDirection
265  map<CdOrientation, CdOrientation> turnMap = prevDirs;
266  if (reverseOrientation)
267  turnMap = nextDirs;
268  map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
269  CdOrientation nextDirection=none;
270  if (nIt != turnMap.end())
271  nextDirection = (*nIt).second;
272  else
273  cout << "ERROR - no Next Direction found!?!?" << endl;
274  return nextDirection;
275  }
EcalBoundaryInfoCalculator::nextDirs
map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:299
EcalBoundaryInfoCalculator::prevDirs
map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:300
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnRight ( CdOrientation  currDirection,
bool  reverseOrientation 
)
inlineprivate

Definition at line 249 of file EcalBoundaryInfoCalculator.h.

249  {
250  //read nextDirection
251  map<CdOrientation, CdOrientation> turnMap = nextDirs;
252  if (reverseOrientation)
253  turnMap = prevDirs;
254  map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
255  CdOrientation nextDirection=none;
256  if (nIt != turnMap.end())
257  nextDirection = (*nIt).second;
258  else
259  cout << "ERROR - no Next Direction found!?!?" << endl;
260  return nextDirection;
261  }
EcalBoundaryInfoCalculator::nextDirs
map< CdOrientation, CdOrientation > nextDirs
Definition: EcalBoundaryInfoCalculator.h:299
EcalBoundaryInfoCalculator::prevDirs
map< CdOrientation, CdOrientation > prevDirs
Definition: EcalBoundaryInfoCalculator.h:300
gather_cfg.cout
tuple cout
Definition: gather_cfg.py:121

Member Data Documentation

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

Definition at line 303 of file EcalBoundaryInfoCalculator.h.

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

Definition at line 299 of file EcalBoundaryInfoCalculator.h.

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

Definition at line 301 of file EcalBoundaryInfoCalculator.h.

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

Definition at line 300 of file EcalBoundaryInfoCalculator.h.

template<class EcalDetId>
CaloNavigator<EcalDetId>* EcalBoundaryInfoCalculator< EcalDetId >::theEcalNav
private

Definition at line 302 of file EcalBoundaryInfoCalculator.h.

Generated for CMSSW Reference Manual by doxygen 1.8.5