1 #ifndef ECALBOUNDARYINFOCALCULATOR_H_
2 #define ECALBOUNDARYINFOCALCULATOR_H_
45 EcalDetId hitdetid = EcalDetId(hit.
id());
51 int hitIeta = ebhitdetid.
ieta();
52 int hitIphi = ebhitdetid.
iphi();
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)
59 int neighbourIeta = hitIeta + ieta;
60 int neighbourIphi = hitIphi + iphi;
62 if (neighbourIphi < 1)
64 if (neighbourIphi > 360)
66 if (neighbourIeta == 0) {
67 neighbourIeta += ieta;
75 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
79 for (vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
86 stati.push_back(status);
95 int hitIx = eehitdetid.
ix();
96 int hitIy = eehitdetid.
iy();
97 int hitIz = eehitdetid.
zside();
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)
104 int neighbourIx = hitIx + ix;
105 int neighbourIy = hitIy + iy;
111 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
114 bool present =
false;
115 for (vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
122 stati.push_back(status);
129 cout <<
"ERROR - RecHit belongs to wrong sub detector" << endl;
132 if (stati.size() > 0)
141 EcalDetId hitdetid = EcalDetId(hit.
id());
147 int hitIeta = ebhitdetid.
ieta();
148 int hitIphi = ebhitdetid.
iphi();
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)
155 int neighbourIeta = hitIeta + ieta;
156 int neighbourIphi = hitIphi + iphi;
158 if (neighbourIphi < 1)
159 neighbourIphi += 360;
160 if (neighbourIphi > 360)
161 neighbourIphi -= 360;
162 if (neighbourIeta == 0) {
163 neighbourIeta += ieta;
176 int hitIx = eehitdetid.
ix();
177 int hitIy = eehitdetid.
iy();
178 int hitIz = eehitdetid.
zside();
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)
185 int neighbourIx = hitIx + ix;
186 int neighbourIy = hitIy + iy;
195 cout <<
"ERROR - RecHit belongs to wrong sub detector" << endl;
202 cout <<
"set Debug Mode!" << endl;
213 next = theNavi->
north();
218 next = theNavi->
east();
223 next = theNavi->
south();
228 next = theNavi->
west();
238 map<CdOrientation, CdOrientation>::iterator oIt = oppositeDirs.find(currDirection);
240 if (oIt != oppositeDirs.end()) {
241 oppDirection = oIt->second;
242 theEcalNav->setHome(prev);
244 EcalDetId currDetId = theEcalNav->pos();
251 map<CdOrientation, CdOrientation> turnMap = nextDirs;
252 if (reverseOrientation)
254 map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
256 if (nIt != turnMap.end())
257 nextDirection = (*nIt).second;
259 cout <<
"ERROR - no Next Direction found!?!?" << endl;
260 return nextDirection;
265 map<CdOrientation, CdOrientation> turnMap = prevDirs;
266 if (reverseOrientation)
268 map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
270 if (nIt != turnMap.end())
271 nextDirection = (*nIt).second;
273 cout <<
"ERROR - no Next Direction found!?!?" << endl;
274 return nextDirection;
280 if (theEcalNav != 0) {
287 if (theEcalNav != 0) {
294 cout <<
"initializeEcalNavigator not implemented for subDet: " << ecalSubDet << endl;
321 oppositeDirs.clear();
341 std::vector<EcalRecHit> boundaryRecHits;
342 std::vector<DetId> boundaryDetIds;
343 std::vector<int> stati;
345 double boundaryEnergy = 0;
346 double boundaryET = 0;
347 int beCellCounter = 0;
348 bool nextToBorder =
false;
350 boundaryRecHits.push_back(*hit);
352 boundaryEnergy += hit->
energy();
353 EcalDetId hitdetid = (EcalDetId) hit->
id();
354 boundaryDetIds.push_back(hitdetid);
358 boundaryET += hit->
energy() / cosh(eta);
361 std::cout <<
"Find Boundary RecHits..." << std::endl;
374 initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
376 bool reverseOrientation =
false;
379 EcalDetId
start = hitdetid;
381 int current_status = 0;
384 bool startAlgo =
false;
387 next = makeStepInDirection(currDirection, theEcalNav);
388 theEcalNav->setHome(current);
391 int status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
393 stati.push_back(status);
397 currDirection = turnLeft(currDirection, reverseOrientation);
401 cout <<
"No starting direction can be found: This should never happen if RecHit has a dead neighbour!" << endl;
408 currDirection = turnRight(currDirection, reverseOrientation);
411 bool nextIsStart =
false;
412 bool atBorder =
false;
414 while (!nextIsStart) {
416 bool nextStepFound =
false;
419 while (!nextStepFound) {
420 next = makeStepInDirection(currDirection, theEcalNav);
421 theEcalNav->setHome(current);
424 status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
427 bool present =
false;
428 for (vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
435 stati.push_back(status);
438 nextStepFound =
true;
440 currDirection = turnRight(currDirection, reverseOrientation);
442 }
else if (next == EcalDetId(0)) {
444 currDirection = turnLeft(currDirection, reverseOrientation);
446 }
else if (status == 0) {
447 nextStepFound =
true;
451 cout <<
"No valid next direction can be found: This should never happen!" << endl;
458 next = makeStepInDirection(currDirection, theEcalNav);
463 cout <<
"Boundary path reached starting position!" << endl;
467 cout <<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId) start << endl;
470 if ((!atBorder || status == 0) && !nextIsStart) {
471 boundaryDetIds.push_back(next);
472 if (RecHits->find(next) != RecHits->end() && status == 0) {
475 boundaryRecHits.push_back(nexthit);
476 boundaryEnergy += nexthit.
energy();
479 boundaryET += nexthit.
energy() / cosh(eta);
483 if (current_status == 0 && status == 0 && atBorder) {
485 currDirection = turnRight(currDirection, reverseOrientation);
488 if (status == 0 && atBorder) {
490 currDirection = turnLeft(currDirection, reverseOrientation);
494 currDirection = turnLeft(currDirection, reverseOrientation);
497 currDirection = turnRight(currDirection, reverseOrientation);
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) {
522 boundInfo.
subdet = hitdetid.subdet();
523 boundInfo.
detIds = boundaryDetIds;
524 boundInfo.
recHits = boundaryRecHits;
530 if (theEcalNav != 0) {
542 std::vector<EcalRecHit> gapRecHits;
543 std::vector<DetId> gapDetIds;
545 double gapEnergy = 0;
547 int gapCellCounter = 0;
548 bool nextToBorder =
false;
550 gapRecHits.push_back(*hit);
552 gapEnergy += hit->
energy();
553 EcalDetId hitdetid = (EcalDetId) hit->
id();
554 gapDetIds.push_back(hitdetid);
558 gapET += hit->
energy() / cosh(eta);
561 std::cout <<
"Find Border RecHits..." << std::endl;
574 initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
576 bool reverseOrientation =
false;
579 EcalDetId
start = hitdetid;
583 bool startAlgo =
false;
586 next = makeStepInDirection(currDirection, theEcalNav);
587 theEcalNav->setHome(start);
589 if (next == EcalDetId(0)) {
594 currDirection = turnLeft(currDirection, reverseOrientation);
598 cout <<
"No starting direction can be found: This should never happen if RecHit is at border!" << endl;
606 currDirection = turnLeft(currDirection, reverseOrientation);
609 bool endIsFound =
false;
610 bool startIsEnd =
false;
612 while (!endIsFound) {
614 bool nextStepFound =
false;
617 while (!nextStepFound) {
618 next = makeStepInDirection(currDirection, theEcalNav);
619 theEcalNav->setHome(current);
622 status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
627 }
else if (next == EcalDetId(0)) {
629 currDirection = turnLeft(currDirection, reverseOrientation);
630 }
else if (status == 0) {
631 if (RecHits->find(next) != RecHits->end()) {
632 nextStepFound =
true;
640 cout <<
"No valid next direction can be found: This should never happen!" << endl;
647 next = makeStepInDirection(currDirection, theEcalNav);
654 cout <<
"Path along gap reached starting position!" << endl;
658 cout <<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId) start << endl;
660 cout <<
"End of gap cluster is found going left" << endl;
665 gapDetIds.push_back(next);
666 if (RecHits->find(next) != RecHits->end()) {
669 gapRecHits.push_back(nexthit);
670 gapEnergy += nexthit.
energy();
673 gapET += nexthit.
energy() / cosh(eta);
678 currDirection = turnRight(currDirection, reverseOrientation);
683 theEcalNav->setHome(start);
686 currDirection = startDirection;
687 currDirection = turnRight(currDirection, reverseOrientation);
694 while (!endIsFound) {
696 bool nextStepFound =
false;
699 while (!nextStepFound) {
700 next = makeStepInDirection(currDirection, theEcalNav);
701 theEcalNav->setHome(current);
704 status = (chit != ecalStatus->end()) ? chit->getStatusCode() & 0x1F : -1;
709 }
else if (next == EcalDetId(0)) {
711 currDirection = turnRight(currDirection, reverseOrientation);
712 }
else if (status == 0) {
713 if (RecHits->find(next) != RecHits->end()) {
714 nextStepFound =
true;
722 cout <<
"No valid next direction can be found: This should never happen!" << endl;
729 next = makeStepInDirection(currDirection, theEcalNav);
733 cout <<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId) start
736 cout <<
"End of gap cluster is found going right" << endl;
741 gapDetIds.push_back(next);
742 if (RecHits->find(next) != RecHits->end()) {
745 gapRecHits.push_back(nexthit);
746 gapEnergy += nexthit.
energy();
749 gapET += nexthit.
energy() / cosh(eta);
754 currDirection = turnLeft(currDirection, reverseOrientation);
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;
768 gapInfo.
subdet = hitdetid.subdet();
769 gapInfo.
detIds = gapDetIds;
777 if (theEcalNav != 0) {
EcalBoundaryInfoCalculator()
static bool validDetId(int i, int j)
check if a valid index combination
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation)
CdOrientation goBackOneCell(CdOrientation currDirection, EcalDetId prev)
BoundaryInformation gapRecHits(const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry)
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 T west() const
move the navigator west
const GlobalPoint & getPosition() const
Returns the position of reference for this cell.
virtual T north() const
move the navigator north
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
virtual const CaloCellGeometry * getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
int iphi() const
get the crystal iphi
bool checkRecHitHasInvalidNeighbour(const EcalRecHit hit, const edm::ESHandle< EcalChannelStatus > ecalStatus)
virtual T east() const
move the navigator east
map< CdOrientation, CdOrientation > nextDirs
CaloNavigator< EcalDetId > * theEcalNav
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation)
bool checkRecHitHasDeadNeighbour(const EcalRecHit hit, const edm::ESHandle< EcalChannelStatus > ecalStatus, vector< int > &stati)
map< CdOrientation, CdOrientation > prevDirs
int ieta() const
get the crystal ieta
static const int ETAPHIMODE
EcalDetId makeStepInDirection(CdOrientation direction, CaloNavigator< EcalDetId > *theNavi)
DetId id() const
get the id
std::vector< Item >::const_iterator const_iterator
ESHandle< TrackerGeometry > geometry
map< CdOrientation, CdOrientation > oppositeDirs
virtual T south() const
move the navigator south
~EcalBoundaryInfoCalculator()