PFClusterProducer/interface/PFECALHashNavigator.h

0001 #ifndef RecoParticleFlow_PFClusterProducer_PFECALHashNavigator_h
0002 #define RecoParticleFlow_PFClusterProducer_PFECALHashNavigator_h
0003
0004 #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h"
0005 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0006 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0007 #include "FWCore/Framework/interface/ConsumesCollector.h"
0008
0009 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0012
0013 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0014 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0015 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0016 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0017
0018 static const CaloDirection orderedDir[8] = {SOUTHWEST, SOUTH, SOUTHEAST, WEST, EAST, NORTHWEST, NORTH, NORTHEAST};
0019
0020 class PFECALHashNavigator : public PFRecHitNavigatorBase {
0021 public:
0022   PFECALHashNavigator() {}
0023
0024   PFECALHashNavigator(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc)
0025       : PFRecHitNavigatorBase(iConfig, cc),
0026         neighbourmapcalculated_(false),
0027         crossBarrelEndcapBorder_(iConfig.getParameter<bool>("crossBarrelEndcapBorder")),
0028         geomToken_(cc.esConsumes<edm::Transition::BeginRun>()) {}
0029
0030   void init(const edm::EventSetup& iSetup) override {
0031     edm::ESHandle<CaloGeometry> geoHandle = iSetup.getHandle(geomToken_);
0032
0033     const CaloSubdetectorGeometry* ebTmp = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalBarrel);
0034     const CaloSubdetectorGeometry* eeTmp = geoHandle->getSubdetectorGeometry(DetId::Ecal, EcalEndcap);
0035
0036     barrelGeometry_ = dynamic_cast<const EcalBarrelGeometry*>(ebTmp);
0037     endcapGeometry_ = dynamic_cast<const EcalEndcapGeometry*>(eeTmp);
0038
0039     // get the ecalBarrel topology
0040     barrelTopology_ = std::make_unique<EcalBarrelTopology>(*geoHandle);
0041     endcapTopology_ = std::make_unique<EcalEndcapTopology>(*geoHandle);
0042
0043     ecalNeighbArray(*barrelGeometry_, *barrelTopology_, *endcapGeometry_, *endcapTopology_);
0044   }
0045
0046   void associateNeighbours(reco::PFRecHit& rh,
0047                            std::unique_ptr<reco::PFRecHitCollection>& hits,
0048                            edm::RefProd<reco::PFRecHitCollection>& refprod) override {
0049     DetId center(rh.detId());
0050
0051     DetId north = move(center, NORTH);
0052     DetId northeast = move(center, NORTHEAST);
0053     DetId northwest = move(center, NORTHWEST);
0054     DetId south = move(center, SOUTH);
0055     DetId southeast = move(center, SOUTHEAST);
0056     DetId southwest = move(center, SOUTHWEST);
0057     DetId east = move(center, EAST);
0058     DetId west = move(center, WEST);
0059
0060     associateNeighbour(north, rh, hits, refprod, 0, 1, 0);
0061     associateNeighbour(northeast, rh, hits, refprod, 1, 1, 0);
0062     associateNeighbour(south, rh, hits, refprod, 0, -1, 0);
0063     associateNeighbour(southwest, rh, hits, refprod, -1, -1, 0);
0064     associateNeighbour(east, rh, hits, refprod, 1, 0, 0);
0065     associateNeighbour(southeast, rh, hits, refprod, 1, -1, 0);
0066     associateNeighbour(west, rh, hits, refprod, -1, 0, 0);
0067     associateNeighbour(northwest, rh, hits, refprod, -1, 1, 0);
0068   }
0069
0070 protected:
0071   void ecalNeighbArray(const EcalBarrelGeometry& barrelGeom,
0072                        const CaloSubdetectorTopology& barrelTopo,
0073                        const EcalEndcapGeometry& endcapGeom,
0074                        const CaloSubdetectorTopology& endcapTopo) {
0075     const unsigned nbarrel = 62000;
0076     // Barrel first. The hashed index runs from 0 to 61199
0077     neighboursEB_.resize(nbarrel);
0078
0079     //std::cout << " Building the array of neighbours (barrel) " ;
0080
0081     const std::vector<DetId>& vec(barrelGeom.getValidDetIds(DetId::Ecal, EcalBarrel));
0082     unsigned size = vec.size();
0083     for (unsigned ic = 0; ic < size; ++ic) {
0084       // We get the 9 cells in a square.
0085       std::vector<DetId> neighbours(barrelTopo.getWindow(vec[ic], 3, 3));
0086       //      std::cout << " Cell " << EBDetId(vec[ic]) << std::endl;
0087       unsigned nneighbours = neighbours.size();
0088
0089       unsigned hashedindex = EBDetId(vec[ic]).hashedIndex();
0090       if (hashedindex >= nbarrel) {
0091         LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0092       }
0093
0094       // If there are 9 cells, it is easy, and this order is know:
0095       //      6  7  8
0096       //      3  4  5
0097       //      0  1  2   (0 = SOUTHWEST)
0098
0099       if (nneighbours == 9) {
0100         neighboursEB_[hashedindex].reserve(8);
0101         for (unsigned in = 0; in < nneighbours; ++in) {
0102           // remove the centre
0103           if (neighbours[in] != vec[ic]) {
0104             neighboursEB_[hashedindex].push_back(neighbours[in]);
0105             //          std::cout << " Neighbour " << in << " " << EBDetId(neighbours[in]) << std::endl;
0106           }
0107         }
0108       } else {
0109         DetId central(vec[ic]);
0110         neighboursEB_[hashedindex].resize(8, DetId(0));
0111         for (unsigned idir = 0; idir < 8; ++idir) {
0112           DetId testid = central;
0113           bool status = stdmove(testid, orderedDir[idir], barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0114           if (status)
0115             neighboursEB_[hashedindex][idir] = testid;
0116         }
0117       }
0118     }
0119
0120     // Moved to the endcap
0121
0122     //  std::cout << " done " << size << std::endl;
0123     //  std::cout << " Building the array of neighbours (endcap) " ;
0124
0125     //  vec.clear();
0126     const std::vector<DetId>& vecee = endcapGeom.getValidDetIds(DetId::Ecal, EcalEndcap);
0127     size = vecee.size();
0128     // There are some holes in the hashedIndex for the EE. Hence the array is bigger than the number
0129     // of crystals
0130     const unsigned nendcap = 19960;
0131
0132     neighboursEE_.resize(nendcap);
0133     for (unsigned ic = 0; ic < size; ++ic) {
0134       // We get the 9 cells in a square.
0135       std::vector<DetId> neighbours(endcapTopo.getWindow(vecee[ic], 3, 3));
0136       unsigned nneighbours = neighbours.size();
0137       // remove the centre
0138       unsigned hashedindex = EEDetId(vecee[ic]).hashedIndex();
0139
0140       if (hashedindex >= nendcap) {
0141         LogDebug("CaloGeometryTools") << " Array overflow " << std::endl;
0142       }
0143
0144       if (nneighbours == 9) {
0145         neighboursEE_[hashedindex].reserve(8);
0146         for (unsigned in = 0; in < nneighbours; ++in) {
0147           // remove the centre
0148           if (neighbours[in] != vecee[ic]) {
0149             neighboursEE_[hashedindex].push_back(neighbours[in]);
0150           }
0151         }
0152       } else {
0153         DetId central(vecee[ic]);
0154         neighboursEE_[hashedindex].resize(8, DetId(0));
0155         for (unsigned idir = 0; idir < 8; ++idir) {
0156           DetId testid = central;
0157           bool status = stdmove(testid, orderedDir[idir], barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0158
0159           if (status)
0160             neighboursEE_[hashedindex][idir] = testid;
0161         }
0162       }
0163     }
0164     //  std::cout << " done " << size <<std::endl;
0165     neighbourmapcalculated_ = true;
0166   }
0167
0168   bool stdsimplemove(DetId& cell,
0169                      const CaloDirection& dir,
0170                      const CaloSubdetectorTopology& barrelTopo,
0171                      const CaloSubdetectorTopology& endcapTopo,
0172                      const EcalBarrelGeometry& barrelGeom,
0173                      const EcalEndcapGeometry& endcapGeom) const {
0174     std::vector<DetId> neighbours;
0175
0176     // BARREL CASE
0177     if (cell.subdetId() == EcalBarrel) {
0178       EBDetId ebDetId = cell;
0179
0180       neighbours = barrelTopo.getNeighbours(ebDetId, dir);
0181
0182       // first try to move according to the standard navigation
0183       if (!neighbours.empty() && !neighbours[0].null()) {
0184         cell = neighbours[0];
0185         return true;
0186       }
0187
0188       // failed.
0189
0190       if (crossBarrelEndcapBorder_) {
0191         // are we on the outer ring ?
0192         const int ietaAbs(ebDetId.ietaAbs());  // abs value of ieta
0193         if (EBDetId::MAX_IETA == ietaAbs) {
0194           // get ee nbrs for for end of barrel crystals
0195
0196           // yes we are
0197           const EcalBarrelGeometry::OrderedListOfEEDetId& ol(*barrelGeom.getClosestEndcapCells(ebDetId));
0198
0199           // take closest neighbour on the other side, that is in the barrel.
0200           cell = *(ol.begin());
0201           return true;
0202         }
0203       }
0204     }
0205
0206     // ENDCAP CASE
0207     else if (cell.subdetId() == EcalEndcap) {
0208       EEDetId eeDetId = cell;
0209
0210       neighbours = endcapTopo.getNeighbours(eeDetId, dir);
0211
0212       if (!neighbours.empty() && !neighbours[0].null()) {
0213         cell = neighbours[0];
0214         return true;
0215       }
0216
0217       // failed.
0218
0219       if (crossBarrelEndcapBorder_) {
0220         // are we on the outer ring ?
0221         const int iphi(eeDetId.iPhiOuterRing());
0222         if (iphi != 0) {
0223           // yes we are
0224           const EcalEndcapGeometry::OrderedListOfEBDetId& ol(*endcapGeom.getClosestBarrelCells(eeDetId));
0225
0226           // take closest neighbour on the other side, that is in the barrel.
0227           cell = *(ol.begin());
0228           return true;
0229         }
0230       }
0231     }
0232
0233     // everything failed
0234     cell = DetId(0);
0235     return false;
0236   }
0237
0238   bool stdmove(DetId& cell,
0239                const CaloDirection& dir,
0240                const CaloSubdetectorTopology& barrelTopo,
0241                const CaloSubdetectorTopology& endcapTopo,
0242                const EcalBarrelGeometry& barrelGeom,
0243                const EcalEndcapGeometry& endcapGeom)
0244
0245       const {
0246     bool result;
0247
0248     if (dir == NORTH) {
0249       result = stdsimplemove(cell, NORTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0250       return result;
0251     } else if (dir == SOUTH) {
0252       result = stdsimplemove(cell, SOUTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0253       return result;
0254     } else if (dir == EAST) {
0255       result = stdsimplemove(cell, EAST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0256       return result;
0257     } else if (dir == WEST) {
0258       result = stdsimplemove(cell, WEST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0259       return result;
0260     }
0261
0262     // One has to try both paths
0263     else if (dir == NORTHEAST) {
0264       result = stdsimplemove(cell, NORTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0265       if (result)
0266         return stdsimplemove(cell, EAST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0267       else {
0268         result = stdsimplemove(cell, EAST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0269         if (result)
0270           return stdsimplemove(cell, NORTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0271         else
0272           return false;
0273       }
0274     } else if (dir == NORTHWEST) {
0275       result = stdsimplemove(cell, NORTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0276       if (result)
0277         return stdsimplemove(cell, WEST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0278       else {
0279         result = stdsimplemove(cell, WEST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0280         if (result)
0281           return stdsimplemove(cell, NORTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0282         else
0283           return false;
0284       }
0285     } else if (dir == SOUTHEAST) {
0286       result = stdsimplemove(cell, SOUTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0287       if (result)
0288         return stdsimplemove(cell, EAST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0289       else {
0290         result = stdsimplemove(cell, EAST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0291         if (result)
0292           return stdsimplemove(cell, SOUTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0293         else
0294           return false;
0295       }
0296     } else if (dir == SOUTHWEST) {
0297       result = stdsimplemove(cell, SOUTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0298       if (result)
0299         return stdsimplemove(cell, WEST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0300       else {
0301         result = stdsimplemove(cell, SOUTH, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0302         if (result)
0303           return stdsimplemove(cell, WEST, barrelTopo, endcapTopo, barrelGeom, endcapGeom);
0304         else
0305           return false;
0306       }
0307     }
0308     cell = DetId(0);
0309     return false;
0310   }
0311
0312   DetId move(DetId cell, const CaloDirection& dir) const {
0313     DetId originalcell = cell;
0314     if (dir == NONE || cell == DetId(0))
0315       return false;
0316
0317     // Conversion CaloDirection and index in the table
0318     // CaloDirection :NONE,SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST, NORTHEAST,NORTHWEST,NORTH
0319     // Table : SOUTHWEST,SOUTH,SOUTHEAST,WEST,EAST,NORTHWEST,NORTH, NORTHEAST
0320     static const int calodirections[9] = {-1, 1, 2, 0, 4, 3, 7, 5, 6};
0321
0322     assert(neighbourmapcalculated_);
0323
0324     DetId result = (originalcell.subdetId() == EcalBarrel)
0325                        ? neighboursEB_[EBDetId(originalcell).hashedIndex()][calodirections[dir]]
0326                        : neighboursEE_[EEDetId(originalcell).hashedIndex()][calodirections[dir]];
0327     return result;
0328   }
0329
0330   std::unique_ptr<EcalEndcapTopology> endcapTopology_;
0331   std::unique_ptr<EcalBarrelTopology> barrelTopology_;
0332
0333   const EcalEndcapGeometry* endcapGeometry_;
0334   const EcalBarrelGeometry* barrelGeometry_;
0335
0336   /// for each ecal barrel rechit, keep track of the neighbours
0337   std::vector<std::vector<DetId> > neighboursEB_;
0338
0339   /// for each ecal endcap rechit, keep track of the neighbours
0340   std::vector<std::vector<DetId> > neighboursEE_;
0341
0342   /// set to true in ecalNeighbArray
0343   bool neighbourmapcalculated_;
0344
0345   /// if true, navigation will cross the barrel-endcap border
0346   bool crossBarrelEndcapBorder_;
0347
0348 private:
0349   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0350 };
0351
0352 #endif