1 #ifndef ECALBOUNDARYINFOCALCULATOR_H_ 2 #define ECALBOUNDARYINFOCALCULATOR_H_ 47 EcalDetId hitdetid = EcalDetId(hit.
id());
53 int hitIeta = ebhitdetid.
ieta();
54 int hitIphi = ebhitdetid.
iphi();
56 for (
int ieta = -1; ieta <= 1; ieta++) {
57 for (
int iphi = -1; iphi <= 1; iphi++) {
58 if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
61 int neighbourIeta = hitIeta + ieta;
62 int neighbourIphi = hitIphi + iphi;
64 if (neighbourIphi < 1)
66 if (neighbourIphi > 360)
68 if (neighbourIeta == 0) {
69 neighbourIeta += ieta;
77 int status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
81 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
88 stati.push_back(status);
97 int hitIx = eehitdetid.
ix();
98 int hitIy = eehitdetid.
iy();
99 int hitIz = eehitdetid.
zside();
101 for (
int ix = -1; ix <= 1; ix++) {
102 for (
int iy = -1; iy <= 1; iy++) {
103 if ((ix == 0 && iy == 0) || ix * iy != 0)
106 int neighbourIx = hitIx + ix;
107 int neighbourIy = hitIy + iy;
113 int status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
116 bool present =
false;
117 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
124 stati.push_back(status);
131 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - RecHit belongs to wrong sub detector";
143 EcalDetId hitdetid = EcalDetId(hit.
id());
149 int hitIeta = ebhitdetid.
ieta();
150 int hitIphi = ebhitdetid.
iphi();
152 for (
int ieta = -1; ieta <= 1; ieta++) {
153 for (
int iphi = -1; iphi <= 1; iphi++) {
154 if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
157 int neighbourIeta = hitIeta + ieta;
158 int neighbourIphi = hitIphi + iphi;
160 if (neighbourIphi < 1)
161 neighbourIphi += 360;
162 if (neighbourIphi > 360)
163 neighbourIphi -= 360;
164 if (neighbourIeta == 0) {
165 neighbourIeta += ieta;
178 int hitIx = eehitdetid.
ix();
179 int hitIy = eehitdetid.
iy();
180 int hitIz = eehitdetid.
zside();
182 for (
int ix = -1; ix <= 1; ix++) {
183 for (
int iy = -1; iy <= 1; iy++) {
184 if ((ix == 0 && iy == 0) || ix * iy != 0)
187 int neighbourIx = hitIx + ix;
188 int neighbourIy = hitIy + iy;
197 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - RecHit belongs to wrong sub detector";
204 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"set Debug Mode!";
214 next = theNavi->
north();
218 next = theNavi->
east();
222 next = theNavi->
south();
226 next = theNavi->
west();
236 std::map<CdOrientation, CdOrientation>::iterator oIt =
oppositeDirs.find(currDirection);
239 oppDirection = oIt->second;
242 EcalDetId currDetId = theEcalNav->
pos();
249 std::map<CdOrientation, CdOrientation> turnMap =
nextDirs;
250 if (reverseOrientation)
252 std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
254 if (nIt != turnMap.end())
255 nextDirection = (*nIt).second;
257 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - no Next Direction found!?!?";
258 return nextDirection;
263 std::map<CdOrientation, CdOrientation> turnMap =
prevDirs;
264 if (reverseOrientation)
266 std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
268 if (nIt != turnMap.end())
269 nextDirection = (*nIt).second;
271 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"ERROR - no Next Direction found!?!?";
272 return nextDirection;
277 std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav(
nullptr);
285 edm::LogWarning(
"EcalBoundaryInfoCalculator") <<
"initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
330 std::vector<DetId> boundaryDetIds;
331 std::vector<int> stati;
333 double boundaryEnergy = 0;
334 double boundaryET = 0;
335 int beCellCounter = 0;
336 bool nextToBorder =
false;
338 boundaryRecHits.push_back(*hit);
340 boundaryEnergy += hit->
energy();
341 EcalDetId hitdetid = (EcalDetId) hit->
id();
342 boundaryDetIds.push_back(hitdetid);
345 double eta = cellGeom->getPosition().eta();
346 boundaryET += hit->
energy() / cosh(eta);
349 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Find Boundary RecHits...";
352 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Starting at : (" << ((
EBDetId) hitdetid).ieta() <<
"," << ((
EBDetId) hitdetid).iphi() <<
")";
354 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Starting at : (" << ((
EEDetId) hitdetid).ix() <<
"," << ((
EEDetId) hitdetid).iy() <<
")";
361 bool reverseOrientation =
false;
364 EcalDetId
start = hitdetid;
365 EcalDetId current =
start;
366 int current_status = 0;
369 bool startAlgo =
false;
373 theEcalNav->setHome(current);
376 int status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
378 stati.push_back(status);
382 currDirection =
turnLeft(currDirection, reverseOrientation);
386 <<
"EcalBoundaryInfoCalculator: No starting direction can be found: This should never happen if RecHit has a dead neighbour!";
391 currDirection =
turnRight(currDirection, reverseOrientation);
394 bool nextIsStart =
false;
395 bool atBorder =
false;
397 while (!nextIsStart) {
399 bool nextStepFound =
false;
402 while (!nextStepFound) {
404 theEcalNav->setHome(current);
407 status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
410 bool present =
false;
411 for (std::vector<int>::const_iterator
s = stati.begin();
s != stati.end(); ++
s) {
418 stati.push_back(status);
421 nextStepFound =
true;
423 currDirection =
turnRight(currDirection, reverseOrientation);
425 }
else if (next == EcalDetId(0)) {
427 currDirection =
turnLeft(currDirection, reverseOrientation);
429 }
else if (status == 0) {
430 nextStepFound =
true;
435 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
445 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Boundary path reached starting position!";
449 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
452 if ((!atBorder || status == 0) && !nextIsStart) {
453 boundaryDetIds.push_back(next);
454 if (RecHits->find(next) != RecHits->end() && status == 0) {
457 boundaryRecHits.push_back(nexthit);
458 boundaryEnergy += nexthit.
energy();
460 eta = cellGeom->getPosition().eta();
461 boundaryET += nexthit.
energy() / cosh(eta);
465 if (current_status == 0 && status == 0 && atBorder) {
467 currDirection =
turnRight(currDirection, reverseOrientation);
470 if (status == 0 && atBorder) {
472 currDirection =
turnLeft(currDirection, reverseOrientation);
476 currDirection =
turnLeft(currDirection, reverseOrientation);
479 currDirection =
turnRight(currDirection, reverseOrientation);
490 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"<<<<<<<<<<<<<<< Final Boundary object <<<<<<<<<<<<<<<";
491 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of neighbouring RecHits: " << boundaryRecHits.size();
492 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of neighbouring DetIds: " << boundaryDetIds.size();
493 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"boundary energy: " << boundaryEnergy;
494 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"boundary ET: " << boundaryET;
495 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"no of cells contributing to boundary energy: " << beCellCounter;
496 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Channel stati: ";
497 for (std::vector<int>::iterator it = stati.begin(); it != stati.end(); ++it) {
498 edm::LogInfo(
"EcalBoundaryInfoCalculator") << *it <<
" ";
504 boundInfo.
subdet = hitdetid.subdet();
505 boundInfo.
detIds = boundaryDetIds;
521 std::vector<DetId> gapDetIds;
523 double gapEnergy = 0;
525 int gapCellCounter = 0;
526 bool nextToBorder =
false;
528 gapRecHits.push_back(*hit);
530 gapEnergy += hit->
energy();
531 EcalDetId hitdetid = (EcalDetId) hit->
id();
532 gapDetIds.push_back(hitdetid);
535 double eta = cellGeom->getPosition().eta();
536 gapET += hit->
energy() / cosh(eta);
539 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Find Border RecHits...";
542 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Starting at : (" << ((
EBDetId) hitdetid).ieta() <<
"," << ((
EBDetId) hitdetid).iphi() <<
")";
544 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Starting at : (" << ((
EEDetId) hitdetid).ix() <<
"," << ((
EEDetId) hitdetid).iy() <<
")";
551 bool reverseOrientation =
false;
554 EcalDetId
start = hitdetid;
555 EcalDetId current =
start;
558 bool startAlgo =
false;
562 theEcalNav->setHome(start);
564 if (next == EcalDetId(0)) {
569 currDirection =
turnLeft(currDirection, reverseOrientation);
573 <<
"EcalBoundaryInfoCalculator: No starting direction can be found: This should never happen if RecHit is at border!";
579 currDirection =
turnLeft(currDirection, reverseOrientation);
582 bool endIsFound =
false;
583 bool startIsEnd =
false;
585 while (!endIsFound) {
587 bool nextStepFound =
false;
590 while (!nextStepFound) {
592 theEcalNav->setHome(current);
595 status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
600 }
else if (next == EcalDetId(0)) {
602 currDirection =
turnLeft(currDirection, reverseOrientation);
603 }
else if (status == 0) {
604 if (RecHits->find(next) != RecHits->end()) {
605 nextStepFound =
true;
614 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
626 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Path along gap reached starting position!";
630 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
632 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"End of gap cluster is found going left";
637 gapDetIds.push_back(next);
638 if (RecHits->find(next) != RecHits->end()) {
641 gapRecHits.push_back(nexthit);
642 gapEnergy += nexthit.
energy();
644 eta = cellGeom->getPosition().eta();
645 gapET += nexthit.
energy() / cosh(eta);
650 currDirection =
turnRight(currDirection, reverseOrientation);
655 theEcalNav->setHome(start);
658 currDirection = startDirection;
659 currDirection =
turnRight(currDirection, reverseOrientation);
666 while (!endIsFound) {
668 bool nextStepFound =
false;
671 while (!nextStepFound) {
673 theEcalNav->setHome(current);
676 status = (chit != ecalStatus->
end()) ? chit->getStatusCode() & 0x1F : -1;
681 }
else if (next == EcalDetId(0)) {
683 currDirection =
turnRight(currDirection, reverseOrientation);
684 }
else if (status == 0) {
685 if (RecHits->find(next) != RecHits->end()) {
686 nextStepFound =
true;
695 <<
"EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
704 edm::LogInfo(
"EcalBoundaryInfoCalculator")<<
"Next step: " << (EcalDetId) next <<
" Status: " << status <<
" Start: " << (EcalDetId)
start;
706 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"End of gap cluster is found going right";
711 gapDetIds.push_back(next);
712 if (RecHits->find(next) != RecHits->end()) {
715 gapRecHits.push_back(nexthit);
716 gapEnergy += nexthit.
energy();
718 eta = cellGeom->getPosition().eta();
719 gapET += nexthit.
energy() / cosh(eta);
724 currDirection =
turnLeft(currDirection, reverseOrientation);
730 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
731 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"No of RecHits along gap: " << gapRecHits.size();
732 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"No of DetIds along gap: " << gapDetIds.size();
733 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Gap energy: " << gapEnergy;
734 edm::LogInfo(
"EcalBoundaryInfoCalculator") <<
"Gap ET: " << gapET;
738 gapInfo.
subdet = hitdetid.subdet();
739 gapInfo.
detIds = gapDetIds;
744 std::vector<int> stati;
EcalBoundaryInfoCalculator()
const CaloSubdetectorGeometry * getSubdetectorGeometry(const DetId &id) const
access the subdetector geometry for the given subdetector directly
bool checkRecHitHasInvalidNeighbour(const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus) const
BoundaryInformation gapRecHits(const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry) const
CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const
std::map< CdOrientation, CdOrientation > nextDirs
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
static bool validDetId(int i, int j)
check if a valid index combination
int iphi() const
get the crystal iphi
void setHome(const T &startingPoint)
set the starting position
CdOrientation goBackOneCell(CdOrientation currDirection, EcalDetId prev, CaloNavigator< EcalDetId > *theEcalNav) const
BoundaryInformation boundaryRecHits(const edm::Handle< EcalRecHitCollection > &, const EcalRecHit *, const edm::ESHandle< CaloTopology > theCaloTopology, const edm::ESHandle< EcalChannelStatus > ecalStatus, const edm::ESHandle< CaloGeometry > geometry) const
T west() const
move the navigator west
T south() const
move the navigator south
int ieta() const
get the crystal ieta
T pos() const
get the current position
static const int ETAPHIMODE
T east() const
move the navigator east
CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const
DetId id() const
get the id
virtual std::shared_ptr< const CaloCellGeometry > getGeometry(const DetId &id) const
Get the cell geometry of a given detector id. Should return false if not found.
static bool validDetId(int crystal_ix, int crystal_iy, int iz)
std::vector< Item >::const_iterator const_iterator
std::map< CdOrientation, CdOrientation > prevDirs
const CaloSubdetectorTopology * getSubdetectorTopology(const DetId &id) const
access the subdetector Topology for the given subdetector directly
std::map< CdOrientation, CdOrientation > oppositeDirs
const_iterator find(uint32_t rawId) const
T north() const
move the navigator north
const_iterator end() const
~EcalBoundaryInfoCalculator()
bool checkRecHitHasDeadNeighbour(const EcalRecHit &hit, const edm::ESHandle< EcalChannelStatus > ecalStatus, std::vector< int > &stati) const