Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-08-23 03:25:25

0001 #ifndef ECALBOUNDARYINFOCALCULATOR_H_
0002 #define ECALBOUNDARYINFOCALCULATOR_H_
0003 #include <memory>
0004 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0005 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0007 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0008 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0009 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0010 #include "RecoCaloTools/Navigation/interface/EcalBarrelNavigator.h"
0011 #include "RecoCaloTools/Navigation/interface/EcalEndcapNavigator.h"
0012 #include "RecoCaloTools/Navigation/interface/EcalPreshowerNavigator.h"
0013 #include "RecoCaloTools/Navigation/interface/CaloTowerNavigator.h"
0014 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0015 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0017 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0018 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0019 #include "DataFormats/METReco/interface/BoundaryInformation.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021 #include "FWCore/Utilities/interface/Exception.h"
0022 
0023 enum CdOrientation { north, east, south, west, none };
0024 
0025 template <class EcalDetId>
0026 class EcalBoundaryInfoCalculator {
0027 public:
0028   EcalBoundaryInfoCalculator();
0029   ~EcalBoundaryInfoCalculator();
0030 
0031   BoundaryInformation boundaryRecHits(const EcalRecHitCollection&,
0032                                       const EcalRecHit*,
0033                                       const CaloTopology& theCaloTopology,
0034                                       const EcalChannelStatus& ecalStatus,
0035                                       const CaloGeometry& geometry) const;
0036 
0037   BoundaryInformation gapRecHits(const EcalRecHitCollection&,
0038                                  const EcalRecHit*,
0039                                  const CaloTopology& theCaloTopology,
0040                                  const EcalChannelStatus& ecalStatus,
0041                                  const CaloGeometry& geometry) const;
0042 
0043   bool checkRecHitHasDeadNeighbour(const EcalRecHit& hit,
0044                                    const EcalChannelStatus& ecalStatus,
0045                                    std::vector<int>& stati) const {
0046     stati.clear();
0047     EcalDetId hitdetid = EcalDetId(hit.id());
0048 
0049     if (hitdetid.subdet() == EcalBarrel) {
0050       EBDetId ebhitdetid = (EBDetId)hitdetid;
0051 
0052       int hitIeta = ebhitdetid.ieta();
0053       int hitIphi = ebhitdetid.iphi();
0054 
0055       for (int ieta = -1; ieta <= 1; ieta++) {
0056         for (int iphi = -1; iphi <= 1; iphi++) {
0057           if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
0058             //if (iphi == 0 && ieta == 0)

0059             continue;
0060           int neighbourIeta = hitIeta + ieta;
0061           int neighbourIphi = hitIphi + iphi;
0062           if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
0063             if (neighbourIphi < 1)
0064               neighbourIphi += 360;
0065             if (neighbourIphi > 360)
0066               neighbourIphi -= 360;
0067             if (neighbourIeta == 0) {
0068               neighbourIeta += ieta;
0069             }
0070           }
0071 
0072           if (EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
0073             const EBDetId detid = EBDetId(neighbourIeta, neighbourIphi, EBDetId::ETAPHIMODE);
0074             EcalChannelStatus::const_iterator chit = ecalStatus.find(detid);
0075             int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0076 
0077             if (status > 0) {
0078               bool present = false;
0079               for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
0080                 if (*s == status) {
0081                   present = true;
0082                   break;
0083                 }
0084               }
0085               if (!present)
0086                 stati.push_back(status);
0087             }
0088           }
0089         }
0090       }
0091 
0092     } else if (hitdetid.subdet() == EcalEndcap) {
0093       EEDetId eehitdetid = (EEDetId)hitdetid;
0094       int hitIx = eehitdetid.ix();
0095       int hitIy = eehitdetid.iy();
0096       int hitIz = eehitdetid.zside();
0097 
0098       for (int ix = -1; ix <= 1; ix++) {
0099         for (int iy = -1; iy <= 1; iy++) {
0100           if ((ix == 0 && iy == 0) || ix * iy != 0)
0101             //if (ix == 0 && iy == 0)

0102             continue;
0103           int neighbourIx = hitIx + ix;
0104           int neighbourIy = hitIy + iy;
0105 
0106           if (EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
0107             const EEDetId detid = EEDetId(neighbourIx, neighbourIy, hitIz, EEDetId::XYMODE);
0108             EcalChannelStatus::const_iterator chit = ecalStatus.find(detid);
0109             int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0110 
0111             if (status > 0) {
0112               bool present = false;
0113               for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
0114                 if (*s == status) {
0115                   present = true;
0116                   break;
0117                 }
0118               }
0119               if (!present)
0120                 stati.push_back(status);
0121             }
0122           }
0123         }
0124       }
0125 
0126     } else {
0127       edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
0128     }
0129 
0130     if (!stati.empty())
0131       return true;
0132     return false;
0133   }
0134 
0135   bool checkRecHitHasInvalidNeighbour(const EcalRecHit& hit, const EcalChannelStatus& ecalStatus) const {
0136     //// return true, if *direct* neighbour is invalid

0137 
0138     EcalDetId hitdetid = EcalDetId(hit.id());
0139 
0140     if (hitdetid.subdet() == EcalBarrel) {
0141       EBDetId ebhitdetid = (EBDetId)hitdetid;
0142 
0143       int hitIeta = ebhitdetid.ieta();
0144       int hitIphi = ebhitdetid.iphi();
0145 
0146       for (int ieta = -1; ieta <= 1; ieta++) {
0147         for (int iphi = -1; iphi <= 1; iphi++) {
0148           if ((iphi == 0 && ieta == 0) || iphi * ieta != 0)
0149             //if (iphi == 0 && ieta == 0)

0150             continue;
0151           int neighbourIeta = hitIeta + ieta;
0152           int neighbourIphi = hitIphi + iphi;
0153           if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
0154             if (neighbourIphi < 1)
0155               neighbourIphi += 360;
0156             if (neighbourIphi > 360)
0157               neighbourIphi -= 360;
0158             if (neighbourIeta == 0) {
0159               neighbourIeta += ieta;
0160             }
0161           }
0162 
0163           if (!EBDetId::validDetId(neighbourIeta, neighbourIphi)) {
0164             return true;
0165           }
0166         }
0167       }
0168 
0169     } else if (hitdetid.subdet() == EcalEndcap) {
0170       EEDetId eehitdetid = (EEDetId)hitdetid;
0171       int hitIx = eehitdetid.ix();
0172       int hitIy = eehitdetid.iy();
0173       int hitIz = eehitdetid.zside();
0174 
0175       for (int ix = -1; ix <= 1; ix++) {
0176         for (int iy = -1; iy <= 1; iy++) {
0177           if ((ix == 0 && iy == 0) || ix * iy != 0)
0178             //if (ix == 0 && iy == 0)

0179             continue;
0180           int neighbourIx = hitIx + ix;
0181           int neighbourIy = hitIy + iy;
0182 
0183           if (!EEDetId::validDetId(neighbourIx, neighbourIy, hitIz)) {
0184             return true;
0185           }
0186         }
0187       }
0188 
0189     } else {
0190       edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - RecHit belongs to wrong sub detector";
0191     }
0192 
0193     return false;
0194   }
0195 
0196   void setDebugMode() {
0197     edm::LogInfo("EcalBoundaryInfoCalculator") << "set Debug Mode!";
0198     debug = true;
0199   }
0200 
0201 private:
0202   EcalDetId makeStepInDirection(CdOrientation direction, const CaloNavigator<EcalDetId>* theNavi) const {
0203     EcalDetId next;
0204     switch (direction) {
0205       case north: {
0206         next = theNavi->north();
0207         break;
0208       }
0209       case east: {
0210         next = theNavi->east();
0211         break;
0212       }
0213       case south: {
0214         next = theNavi->south();
0215         break;
0216       }
0217       case west: {
0218         next = theNavi->west();
0219         break;
0220       }
0221       default:
0222         break;
0223     }
0224     return next;
0225   }
0226 
0227   CdOrientation goBackOneCell(CdOrientation currDirection, EcalDetId prev, CaloNavigator<EcalDetId>* theEcalNav) const {
0228     auto oIt = oppositeDirs.find(currDirection);
0229     CdOrientation oppDirection = none;
0230     if (oIt != oppositeDirs.end()) {
0231       oppDirection = oIt->second;
0232       theEcalNav->setHome(prev);
0233     }
0234     EcalDetId currDetId = theEcalNav->pos();
0235 
0236     return oppDirection;
0237   }
0238 
0239   CdOrientation turnRight(CdOrientation currDirection, bool reverseOrientation) const {
0240     //read nextDirection

0241     std::map<CdOrientation, CdOrientation> turnMap = nextDirs;
0242     if (reverseOrientation)
0243       turnMap = prevDirs;
0244     std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
0245     CdOrientation nextDirection = none;
0246     if (nIt != turnMap.end())
0247       nextDirection = (*nIt).second;
0248     else
0249       edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
0250     return nextDirection;
0251   }
0252 
0253   CdOrientation turnLeft(CdOrientation currDirection, bool reverseOrientation) const {
0254     //read nextDirection

0255     std::map<CdOrientation, CdOrientation> turnMap = prevDirs;
0256     if (reverseOrientation)
0257       turnMap = nextDirs;
0258     std::map<CdOrientation, CdOrientation>::iterator nIt = turnMap.find(currDirection);
0259     CdOrientation nextDirection = none;
0260     if (nIt != turnMap.end())
0261       nextDirection = (*nIt).second;
0262     else
0263       edm::LogWarning("EcalBoundaryInfoCalculator") << "ERROR - no Next Direction found!?!?";
0264     return nextDirection;
0265   }
0266 
0267   std::unique_ptr<CaloNavigator<EcalDetId>> initializeEcalNavigator(DetId startE,
0268                                                                     const CaloTopology& theCaloTopology,
0269                                                                     EcalSubdetector ecalSubDet) const {
0270     std::unique_ptr<CaloNavigator<EcalDetId>> theEcalNav;
0271     if (ecalSubDet == EcalBarrel) {
0272       theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
0273           (EBDetId)startE, (theCaloTopology.getSubdetectorTopology(DetId::Ecal, ecalSubDet)));
0274     } else if (ecalSubDet == EcalEndcap) {
0275       theEcalNav = std::make_unique<CaloNavigator<EcalDetId>>(
0276           (EEDetId)startE, (theCaloTopology.getSubdetectorTopology(DetId::Ecal, ecalSubDet)));
0277     } else {
0278       edm::LogWarning("EcalBoundaryInfoCalculator")
0279           << "initializeEcalNavigator not implemented for subDet: " << ecalSubDet;
0280     }
0281     return theEcalNav;
0282   }
0283 
0284   std::map<CdOrientation, CdOrientation> nextDirs;
0285   std::map<CdOrientation, CdOrientation> prevDirs;
0286   std::map<CdOrientation, CdOrientation> oppositeDirs;
0287   bool debug;
0288 };
0289 
0290 template <class EcalDetId>
0291 EcalBoundaryInfoCalculator<EcalDetId>::EcalBoundaryInfoCalculator() {
0292   nextDirs.clear();
0293   nextDirs[north] = east;
0294   nextDirs[east] = south;
0295   nextDirs[south] = west;
0296   nextDirs[west] = north;
0297 
0298   prevDirs.clear();
0299   prevDirs[north] = west;
0300   prevDirs[east] = north;
0301   prevDirs[south] = east;
0302   prevDirs[west] = south;
0303 
0304   oppositeDirs.clear();
0305   oppositeDirs[north] = south;
0306   oppositeDirs[south] = north;
0307   oppositeDirs[east] = west;
0308   oppositeDirs[west] = east;
0309 
0310   debug = false;
0311 }
0312 
0313 template <class EcalDetId>
0314 EcalBoundaryInfoCalculator<EcalDetId>::~EcalBoundaryInfoCalculator() {}
0315 
0316 template <class EcalDetId>
0317 BoundaryInformation EcalBoundaryInfoCalculator<EcalDetId>::boundaryRecHits(const EcalRecHitCollection& RecHits,
0318                                                                            const EcalRecHit* hit,
0319                                                                            const CaloTopology& theCaloTopology,
0320                                                                            const EcalChannelStatus& ecalStatus,
0321                                                                            const CaloGeometry& geometry) const {
0322   //initialize boundary information

0323   std::vector<EcalRecHit> boundaryRecHits;
0324   std::vector<DetId> boundaryDetIds;
0325   std::vector<int> stati;
0326 
0327   double boundaryEnergy = 0;
0328   double boundaryET = 0;
0329   int beCellCounter = 0;
0330   bool nextToBorder = false;
0331 
0332   boundaryRecHits.push_back(*hit);
0333   ++beCellCounter;
0334   boundaryEnergy += hit->energy();
0335   EcalDetId hitdetid = (EcalDetId)hit->id();
0336   boundaryDetIds.push_back(hitdetid);
0337   const CaloSubdetectorGeometry* subGeom = geometry.getSubdetectorGeometry(hitdetid);
0338   auto cellGeom = subGeom->getGeometry(hitdetid);
0339   double eta = cellGeom->getPosition().eta();
0340   boundaryET += hit->energy() / cosh(eta);
0341 
0342   if (debug) {
0343     edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Boundary RecHits...";
0344 
0345     if (hitdetid.subdet() == EcalBarrel) {
0346       edm::LogInfo("EcalBoundaryInfoCalculator")
0347           << "Starting at : (" << ((EBDetId)hitdetid).ieta() << "," << ((EBDetId)hitdetid).iphi() << ")";
0348     } else if (hitdetid.subdet() == EcalEndcap) {
0349       edm::LogInfo("EcalBoundaryInfoCalculator")
0350           << "Starting at : (" << ((EEDetId)hitdetid).ix() << "," << ((EEDetId)hitdetid).iy() << ")";
0351     }
0352   }
0353 
0354   //initialize navigator

0355   auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
0356   CdOrientation currDirection = north;
0357   bool reverseOrientation = false;
0358 
0359   EcalDetId next(0);
0360   EcalDetId start = hitdetid;
0361   EcalDetId current = start;
0362   int current_status = 0;
0363 
0364   // Search until a dead cell is ahead

0365   bool startAlgo = false;
0366   int noDirs = 0;
0367   while (!startAlgo) {
0368     next = makeStepInDirection(currDirection, theEcalNav.get());
0369     theEcalNav->setHome(current);
0370     theEcalNav->home();
0371     EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
0372     int status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0373     if (status > 0) {
0374       stati.push_back(status);
0375       startAlgo = true;
0376       break;
0377     }
0378     currDirection = turnLeft(currDirection, reverseOrientation);
0379     ++noDirs;
0380     if (noDirs > 4) {
0381       throw cms::Exception("NoStartingDirection") << "EcalBoundaryInfoCalculator: No starting direction can be found: "
0382                                                      "This should never happen if RecHit has a dead neighbour!";
0383     }
0384   }
0385 
0386   // go around dead clusters counter clock wise

0387   currDirection = turnRight(currDirection, reverseOrientation);
0388 
0389   // Search for next boundary element

0390   bool nextIsStart = false;
0391   bool atBorder = false;
0392 
0393   while (!nextIsStart) {
0394     bool nextStepFound = false;
0395     int status = 0;
0396     noDirs = 0;
0397     while (!nextStepFound) {
0398       next = makeStepInDirection(currDirection, theEcalNav.get());
0399       theEcalNav->setHome(current);
0400       theEcalNav->home();
0401       EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
0402       status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0403       if (status > 0) {
0404         // New dead cell found: update status std::vector of dead channels

0405         bool present = false;
0406         for (std::vector<int>::const_iterator s = stati.begin(); s != stati.end(); ++s) {
0407           if (*s == status) {
0408             present = true;
0409             break;
0410           }
0411         }
0412         if (!present)
0413           stati.push_back(status);
0414 
0415         if (atBorder) {
0416           nextStepFound = true;
0417         } else {
0418           currDirection = turnRight(currDirection, reverseOrientation);
0419         }
0420       } else if (next == EcalDetId(0)) {
0421         // In case the Ecal border is reached -> go along dead cells

0422         currDirection = turnLeft(currDirection, reverseOrientation);
0423         atBorder = true;
0424       } else if (status == 0) {
0425         nextStepFound = true;
0426       }
0427       ++noDirs;
0428       if (noDirs > 4) {
0429         throw cms::Exception("NoNextDirection")
0430             << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
0431       }
0432     }
0433 
0434     // make next step

0435     next = makeStepInDirection(currDirection, theEcalNav.get());
0436 
0437     if (next == start) {
0438       nextIsStart = true;
0439       if (debug)
0440         edm::LogInfo("EcalBoundaryInfoCalculator") << "Boundary path reached starting position!";
0441     }
0442 
0443     if (debug)
0444       edm::LogInfo("EcalBoundaryInfoCalculator")
0445           << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
0446 
0447     // save recHits and add energy if on the boundary (and not inside at border)

0448     if ((!atBorder || status == 0) && !nextIsStart) {
0449       boundaryDetIds.push_back(next);
0450       if (RecHits.find(next) != RecHits.end() && status == 0) {
0451         EcalRecHit nexthit = *RecHits.find(next);
0452         ++beCellCounter;
0453         boundaryRecHits.push_back(nexthit);
0454         boundaryEnergy += nexthit.energy();
0455         cellGeom = subGeom->getGeometry(hitdetid);
0456         eta = cellGeom->getPosition().eta();
0457         boundaryET += nexthit.energy() / cosh(eta);
0458       }
0459     }
0460 
0461     if (current_status == 0 && status == 0 && atBorder) {
0462       // this is for a special case, where dead cells are at border corner

0463       currDirection = turnRight(currDirection, reverseOrientation);
0464     } else {
0465       // if dead region along a border is left, turn left

0466       if (status == 0 && atBorder) {
0467         atBorder = false;
0468         currDirection = turnLeft(currDirection, reverseOrientation);
0469       }
0470       if (status == 0) {
0471         // if outside the cluster turn left to follow boundary

0472         currDirection = turnLeft(currDirection, reverseOrientation);
0473       } else {
0474         // else turn right to check if dead region can be left

0475         currDirection = turnRight(currDirection, reverseOrientation);
0476       }
0477     }
0478 
0479     // save currect position

0480     current = next;
0481     current_status = status;
0482   }
0483 
0484   if (debug) {
0485     edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Boundary object <<<<<<<<<<<<<<<";
0486     edm::LogInfo("EcalBoundaryInfoCalculator") << "no of neighbouring RecHits: " << boundaryRecHits.size();
0487     edm::LogInfo("EcalBoundaryInfoCalculator") << "no of neighbouring DetIds: " << boundaryDetIds.size();
0488     edm::LogInfo("EcalBoundaryInfoCalculator") << "boundary energy: " << boundaryEnergy;
0489     edm::LogInfo("EcalBoundaryInfoCalculator") << "boundary ET: " << boundaryET;
0490     edm::LogInfo("EcalBoundaryInfoCalculator") << "no of cells contributing to boundary energy: " << beCellCounter;
0491     edm::LogInfo("EcalBoundaryInfoCalculator") << "Channel stati: ";
0492     for (std::vector<int>::iterator it = stati.begin(); it != stati.end(); ++it) {
0493       edm::LogInfo("EcalBoundaryInfoCalculator") << *it << " ";
0494     }
0495     edm::LogInfo("EcalBoundaryInfoCalculator");
0496   }
0497 
0498   BoundaryInformation boundInfo;
0499   boundInfo.subdet = hitdetid.subdet();
0500   boundInfo.detIds = boundaryDetIds;
0501   boundInfo.recHits = boundaryRecHits;
0502   boundInfo.boundaryEnergy = boundaryEnergy;
0503   boundInfo.boundaryET = boundaryET;
0504   boundInfo.nextToBorder = nextToBorder;
0505   boundInfo.channelStatus = stati;
0506 
0507   return boundInfo;
0508 }
0509 
0510 template <class EcalDetId>
0511 BoundaryInformation EcalBoundaryInfoCalculator<EcalDetId>::gapRecHits(const EcalRecHitCollection& RecHits,
0512                                                                       const EcalRecHit* hit,
0513                                                                       const CaloTopology& theCaloTopology,
0514                                                                       const EcalChannelStatus& ecalStatus,
0515                                                                       const CaloGeometry& geometry) const {
0516   //initialize boundary information

0517   std::vector<EcalRecHit> gapRecHits;
0518   std::vector<DetId> gapDetIds;
0519 
0520   double gapEnergy = 0;
0521   double gapET = 0;
0522   int gapCellCounter = 0;
0523   bool nextToBorder = false;
0524 
0525   gapRecHits.push_back(*hit);
0526   ++gapCellCounter;
0527   gapEnergy += hit->energy();
0528   EcalDetId hitdetid = (EcalDetId)hit->id();
0529   gapDetIds.push_back(hitdetid);
0530   const CaloSubdetectorGeometry* subGeom = geometry.getSubdetectorGeometry(hitdetid);
0531   auto cellGeom = subGeom->getGeometry(hitdetid);
0532   double eta = cellGeom->getPosition().eta();
0533   gapET += hit->energy() / cosh(eta);
0534 
0535   if (debug) {
0536     edm::LogInfo("EcalBoundaryInfoCalculator") << "Find Border RecHits...";
0537 
0538     if (hitdetid.subdet() == EcalBarrel) {
0539       edm::LogInfo("EcalBoundaryInfoCalculator")
0540           << "Starting at : (" << ((EBDetId)hitdetid).ieta() << "," << ((EBDetId)hitdetid).iphi() << ")";
0541     } else if (hitdetid.subdet() == EcalEndcap) {
0542       edm::LogInfo("EcalBoundaryInfoCalculator")
0543           << "Starting at : (" << ((EEDetId)hitdetid).ix() << "," << ((EEDetId)hitdetid).iy() << ")";
0544     }
0545   }
0546 
0547   //initialize navigator

0548   auto theEcalNav = initializeEcalNavigator(hitdetid, theCaloTopology, EcalDetId::subdet());
0549   CdOrientation currDirection = north;
0550   bool reverseOrientation = false;
0551 
0552   EcalDetId next(0);
0553   EcalDetId start = hitdetid;
0554   EcalDetId current = start;
0555 
0556   // Search until a invalid cell is ahead

0557   bool startAlgo = false;
0558   int noDirs = 0;
0559   while (!startAlgo) {
0560     next = makeStepInDirection(currDirection, theEcalNav.get());
0561     theEcalNav->setHome(start);
0562     theEcalNav->home();
0563     if (next == EcalDetId(0)) {
0564       startAlgo = true;
0565       nextToBorder = true;
0566       break;
0567     }
0568     currDirection = turnLeft(currDirection, reverseOrientation);
0569     ++noDirs;
0570     if (noDirs > 4) {
0571       throw cms::Exception("NoStartingDirection") << "EcalBoundaryInfoCalculator: No starting direction can be found: "
0572                                                      "This should never happen if RecHit is at border!";
0573     }
0574   }
0575 
0576   ////////// First: go along gap left

0577   CdOrientation startDirection = currDirection;
0578   currDirection = turnLeft(currDirection, reverseOrientation);
0579 
0580   // Search for next border element

0581   bool endIsFound = false;
0582   bool startIsEnd = false;
0583 
0584   while (!endIsFound) {
0585     bool nextStepFound = false;
0586     int status = 0;
0587     noDirs = 0;
0588     while (!nextStepFound) {
0589       next = makeStepInDirection(currDirection, theEcalNav.get());
0590       theEcalNav->setHome(current);
0591       theEcalNav->home();
0592       EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
0593       status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0594       if (status > 0) {
0595         // Find dead cell along border -> end of cluster

0596         endIsFound = true;
0597         break;
0598       } else if (next == EcalDetId(0)) {
0599         // In case the Ecal border -> go along gap

0600         currDirection = turnLeft(currDirection, reverseOrientation);
0601       } else if (status == 0) {
0602         if (RecHits.find(next) != RecHits.end()) {
0603           nextStepFound = true;
0604         } else {
0605           endIsFound = true;
0606           break;
0607         }
0608       }
0609       ++noDirs;
0610       if (noDirs > 4) {
0611         throw cms::Exception("NoNextDirection")
0612             << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
0613       }
0614     }
0615 
0616     // make next step

0617     next = makeStepInDirection(currDirection, theEcalNav.get());
0618     current = next;
0619 
0620     if (next == start) {
0621       startIsEnd = true;
0622       endIsFound = true;
0623       if (debug)
0624         edm::LogInfo("EcalBoundaryInfoCalculator") << "Path along gap reached starting position!";
0625     }
0626 
0627     if (debug) {
0628       edm::LogInfo("EcalBoundaryInfoCalculator")
0629           << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
0630       if (endIsFound)
0631         edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going left";
0632     }
0633 
0634     // save recHits and add energy

0635     if (!endIsFound) {
0636       gapDetIds.push_back(next);
0637       if (RecHits.find(next) != RecHits.end()) {
0638         EcalRecHit nexthit = *RecHits.find(next);
0639         ++gapCellCounter;
0640         gapRecHits.push_back(nexthit);
0641         gapEnergy += nexthit.energy();
0642         cellGeom = subGeom->getGeometry(next);
0643         eta = cellGeom->getPosition().eta();
0644         gapET += nexthit.energy() / cosh(eta);
0645       }
0646     }
0647 
0648     // turn right to follow gap

0649     currDirection = turnRight(currDirection, reverseOrientation);
0650   }
0651 
0652   ////////// Second: go along gap right

0653   theEcalNav->setHome(start);
0654   theEcalNav->home();
0655   current = start;
0656   currDirection = startDirection;
0657   currDirection = turnRight(currDirection, reverseOrientation);
0658 
0659   // Search for next border element

0660   endIsFound = false;
0661 
0662   if (!startIsEnd) {
0663     while (!endIsFound) {
0664       bool nextStepFound = false;
0665       int status = 0;
0666       noDirs = 0;
0667       while (!nextStepFound) {
0668         next = makeStepInDirection(currDirection, theEcalNav.get());
0669         theEcalNav->setHome(current);
0670         theEcalNav->home();
0671         EcalChannelStatus::const_iterator chit = ecalStatus.find(next);
0672         status = (chit != ecalStatus.end()) ? chit->getStatusCode() & 0x1F : -1;
0673         if (status > 0) {
0674           // Find dead cell along border -> end of cluster

0675           endIsFound = true;
0676           break;
0677         } else if (next == EcalDetId(0)) {
0678           // In case the Ecal border -> go along gap

0679           currDirection = turnRight(currDirection, reverseOrientation);
0680         } else if (status == 0) {
0681           if (RecHits.find(next) != RecHits.end()) {
0682             nextStepFound = true;
0683           } else {
0684             endIsFound = true;
0685             break;
0686           }
0687         }
0688         ++noDirs;
0689         if (noDirs > 4) {
0690           throw cms::Exception("NoStartingDirection")
0691               << "EcalBoundaryInfoCalculator: No valid next direction can be found: This should never happen!";
0692         }
0693       }
0694 
0695       // make next step

0696       next = makeStepInDirection(currDirection, theEcalNav.get());
0697       current = next;
0698 
0699       if (debug) {
0700         edm::LogInfo("EcalBoundaryInfoCalculator")
0701             << "Next step: " << (EcalDetId)next << " Status: " << status << " Start: " << (EcalDetId)start;
0702         if (endIsFound)
0703           edm::LogInfo("EcalBoundaryInfoCalculator") << "End of gap cluster is found going right";
0704       }
0705 
0706       // save recHits and add energy

0707       if (!endIsFound) {
0708         gapDetIds.push_back(next);
0709         if (RecHits.find(next) != RecHits.end()) {
0710           EcalRecHit nexthit = *RecHits.find(next);
0711           ++gapCellCounter;
0712           gapRecHits.push_back(nexthit);
0713           gapEnergy += nexthit.energy();
0714           cellGeom = subGeom->getGeometry(next);
0715           eta = cellGeom->getPosition().eta();
0716           gapET += nexthit.energy() / cosh(eta);
0717         }
0718       }
0719 
0720       // turn left to follow gap

0721       currDirection = turnLeft(currDirection, reverseOrientation);
0722     }
0723   }
0724 
0725   if (debug) {
0726     edm::LogInfo("EcalBoundaryInfoCalculator") << "<<<<<<<<<<<<<<< Final Gap object <<<<<<<<<<<<<<<";
0727     edm::LogInfo("EcalBoundaryInfoCalculator") << "No of RecHits along gap: " << gapRecHits.size();
0728     edm::LogInfo("EcalBoundaryInfoCalculator") << "No of DetIds along gap: " << gapDetIds.size();
0729     edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap energy: " << gapEnergy;
0730     edm::LogInfo("EcalBoundaryInfoCalculator") << "Gap ET: " << gapET;
0731   }
0732 
0733   BoundaryInformation gapInfo;
0734   gapInfo.subdet = hitdetid.subdet();
0735   gapInfo.detIds = gapDetIds;
0736   gapInfo.recHits = gapRecHits;
0737   gapInfo.boundaryEnergy = gapEnergy;
0738   gapInfo.boundaryET = gapET;
0739   gapInfo.nextToBorder = nextToBorder;
0740   std::vector<int> stati;
0741   gapInfo.channelStatus = stati;
0742 
0743   return gapInfo;
0744 }
0745 
0746 #endif /*ECALBOUNDARYINFOCALCULATOR_H_*/