Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:26:43 

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

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

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

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

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

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

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

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

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

0646     currDirection = turnRight(currDirection, reverseOrientation);
0647   }
0648 
0649   ////////// Second: go along gap right

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

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

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

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

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

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

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!