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, std::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
 
std::map< CdOrientation,
CdOrientation
nextDirs
 
std::map< CdOrientation,
CdOrientation
oppositeDirs
 
std::map< CdOrientation,
CdOrientation
prevDirs
 
CaloNavigator< EcalDetId > * theEcalNav
 

Detailed Description

template<class EcalDetId>
class EcalBoundaryInfoCalculator< EcalDetId >

Definition at line 24 of file EcalBoundaryInfoCalculator.h.

Constructor & Destructor Documentation

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

Definition at line 305 of file EcalBoundaryInfoCalculator.h.

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

305  {
306 
307  nextDirs.clear();
308  nextDirs[north] = east;
309  nextDirs[east] = south;
310  nextDirs[south] = west;
311  nextDirs[west] = north;
312 
313  prevDirs.clear();
314  prevDirs[north] = west;
315  prevDirs[east] = north;
316  prevDirs[south] = east;
317  prevDirs[west] = south;
318 
319  oppositeDirs.clear();
324 
325  theEcalNav = 0;
326  debug = false;
327 
328 }
std::map< CdOrientation, CdOrientation > nextDirs
CaloNavigator< EcalDetId > * theEcalNav
std::map< CdOrientation, CdOrientation > prevDirs
std::map< CdOrientation, CdOrientation > oppositeDirs
template<class EcalDetId >
EcalBoundaryInfoCalculator< EcalDetId >::~EcalBoundaryInfoCalculator ( )

Definition at line 330 of file EcalBoundaryInfoCalculator.h.

330  {
331  delete theEcalNav;
332 }
CaloNavigator< EcalDetId > * theEcalNav

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

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

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 39 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

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

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 275 of file EcalBoundaryInfoCalculator.h.

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

Definition at line 206 of file EcalBoundaryInfoCalculator.h.

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

Definition at line 199 of file EcalBoundaryInfoCalculator.h.

Referenced by EcalDeadCellBoundaryEnergyFilter::filter().

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

Definition at line 261 of file EcalBoundaryInfoCalculator.h.

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  std::cout << "ERROR - no Next Direction found!?!?" << std::endl;
272  return nextDirection;
273  }
std::map< CdOrientation, CdOrientation > nextDirs
std::map< CdOrientation, CdOrientation > prevDirs
tuple cout
Definition: gather_cfg.py:121
template<class EcalDetId>
CdOrientation EcalBoundaryInfoCalculator< EcalDetId >::turnRight ( CdOrientation  currDirection,
bool  reverseOrientation 
)
inlineprivate

Definition at line 247 of file EcalBoundaryInfoCalculator.h.

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  std::cout << "ERROR - no Next Direction found!?!?" << std::endl;
258  return nextDirection;
259  }
std::map< CdOrientation, CdOrientation > nextDirs
std::map< CdOrientation, CdOrientation > prevDirs
tuple cout
Definition: gather_cfg.py:121

Member Data Documentation

template<class EcalDetId>
bool EcalBoundaryInfoCalculator< EcalDetId >::debug
private
template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::nextDirs
private
template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::oppositeDirs
private
template<class EcalDetId>
std::map<CdOrientation, CdOrientation> EcalBoundaryInfoCalculator< EcalDetId >::prevDirs
private
template<class EcalDetId>
CaloNavigator<EcalDetId>* EcalBoundaryInfoCalculator< EcalDetId >::theEcalNav
private