Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-05-18 09:15:19

0001 #ifndef RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h
0002 #define RecoParticleFlow_PFClusterProducer_PFHCALDenseIdNavigator_h
0003 
0004 #include "FWCore/Framework/interface/ESWatcher.h"
0005 
0006 #include "RecoParticleFlow/PFClusterProducer/interface/PFRecHitNavigatorBase.h"
0007 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0008 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0009 #include "FWCore/Framework/interface/ConsumesCollector.h"
0010 
0011 #include "RecoCaloTools/Navigation/interface/CaloNavigator.h"
0012 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0013 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0014 #include "DataFormats/EcalDetId/interface/ESDetId.h"
0015 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0016 
0017 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0018 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0019 #include "Geometry/CaloTopology/interface/EcalPreshowerTopology.h"
0020 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0021 
0022 #include "Geometry/CaloTopology/interface/CaloTowerTopology.h"
0023 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0024 
0025 template <typename DET, typename TOPO, bool ownsTopo = true>
0026 class PFHCALDenseIdNavigator : public PFRecHitNavigatorBase {
0027 public:
0028   ~PFHCALDenseIdNavigator() override {
0029     if (!ownsTopo) {
0030       topology_.release();
0031     }
0032   }
0033 
0034   PFHCALDenseIdNavigator(const edm::ParameterSet& iConfig, edm::ConsumesCollector& cc)
0035       : vhcalEnum_(iConfig.getParameter<std::vector<int>>("hcalEnums")),
0036         hcalToken_(cc.esConsumes<edm::Transition::BeginLuminosityBlock>()),
0037         geomToken_(cc.esConsumes<edm::Transition::BeginLuminosityBlock>()) {}
0038 
0039   void init(const edm::EventSetup& iSetup) override {
0040     bool check = theRecNumberWatcher_.check(iSetup);
0041     if (!check)
0042       return;
0043 
0044     edm::ESHandle<HcalTopology> hcalTopology = iSetup.getHandle(hcalToken_);
0045     topology_.release();
0046     topology_.reset(hcalTopology.product());
0047 
0048     // Fill a vector of valid denseid's
0049     edm::ESHandle<CaloGeometry> hGeom = iSetup.getHandle(geomToken_);
0050     const CaloGeometry& caloGeom = *hGeom;
0051 
0052     std::vector<DetId> vecHcal;
0053     std::vector<unsigned int> vDenseIdHcal;
0054     neighboursHcal_.clear();
0055     for (auto hcalSubdet : vhcalEnum_) {
0056       std::vector<DetId> vecDetIds(caloGeom.getValidDetIds(DetId::Hcal, hcalSubdet));
0057       vecHcal.insert(vecHcal.end(), vecDetIds.begin(), vecDetIds.end());
0058     }
0059     vDenseIdHcal.reserve(vecHcal.size());
0060     for (auto hDetId : vecHcal) {
0061       vDenseIdHcal.push_back(topology_.get()->detId2denseId(hDetId));
0062     }
0063     std::sort(vDenseIdHcal.begin(), vDenseIdHcal.end());
0064 
0065     // Fill a vector of cell neighbours
0066     denseIdHcalMax_ = *max_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
0067     denseIdHcalMin_ = *min_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
0068     neighboursHcal_.resize(denseIdHcalMax_ - denseIdHcalMin_ + 1);
0069 
0070     for (auto denseid : vDenseIdHcal) {
0071       DetId N(0);
0072       DetId E(0);
0073       DetId S(0);
0074       DetId W(0);
0075       DetId NW(0);
0076       DetId NE(0);
0077       DetId SW(0);
0078       DetId SE(0);
0079       std::vector<DetId> neighbours(9, DetId(0));
0080 
0081       // the centre
0082       unsigned denseid_c = denseid;
0083       DetId detid_c = topology_.get()->denseId2detId(denseid_c);
0084       CaloNavigator<DET> navigator(detid_c, topology_.get());
0085 
0086       // Using enum in Geometry/CaloTopology/interface/CaloDirection.h
0087       // Order: CENTER(NONE),SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST,NORTHEAST,NORTHWEST,NORTH
0088       neighbours.at(NONE) = detid_c;
0089 
0090       navigator.home();
0091       N = navigator.north();
0092       neighbours.at(NORTH) = N;
0093       if (N != DetId(0)) {
0094         NE = navigator.east();
0095       } else {
0096         navigator.home();
0097         E = navigator.east();
0098         NE = navigator.north();
0099       }
0100       neighbours.at(NORTHEAST) = NE;
0101 
0102       navigator.home();
0103       S = navigator.south();
0104       neighbours.at(SOUTH) = S;
0105       if (S != DetId(0)) {
0106         SW = navigator.west();
0107       } else {
0108         navigator.home();
0109         W = navigator.west();
0110         SW = navigator.south();
0111       }
0112       neighbours.at(SOUTHWEST) = SW;
0113 
0114       navigator.home();
0115       E = navigator.east();
0116       neighbours.at(EAST) = E;
0117       if (E != DetId(0)) {
0118         SE = navigator.south();
0119       } else {
0120         navigator.home();
0121         S = navigator.south();
0122         SE = navigator.east();
0123       }
0124       neighbours.at(SOUTHEAST) = SE;
0125 
0126       navigator.home();
0127       W = navigator.west();
0128       neighbours.at(WEST) = W;
0129       if (W != DetId(0)) {
0130         NW = navigator.north();
0131       } else {
0132         navigator.home();
0133         N = navigator.north();
0134         NW = navigator.west();
0135       }
0136       neighbours.at(NORTHWEST) = NW;
0137 
0138       unsigned index = getIdx(denseid_c);
0139       neighboursHcal_[index] = neighbours;
0140     }
0141   }
0142 
0143   void associateNeighbours(reco::PFRecHit& hit,
0144                            std::unique_ptr<reco::PFRecHitCollection>& hits,
0145                            edm::RefProd<reco::PFRecHitCollection>& refProd) override {
0146     DetId detid(hit.detId());
0147     unsigned denseid = topology_.get()->detId2denseId(detid);
0148 
0149     std::vector<DetId> neighbours(9, DetId(0));
0150 
0151     if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_) {
0152       edm::LogWarning("PFRecHitHCALCachedNavigator") << " DenseId for this cell is out of the range." << std::endl;
0153     } else if (!validNeighbours(denseid)) {
0154       edm::LogWarning("PFRecHitHCALCachedNavigator")
0155           << " DenseId for this cell does not have the neighbour information." << std::endl;
0156     } else {
0157       unsigned index = getIdx(denseid);
0158       neighbours = neighboursHcal_.at(index);
0159     }
0160 
0161     associateNeighbour(neighbours.at(NORTH), hit, hits, refProd, 0, 1, 0);        // N
0162     associateNeighbour(neighbours.at(NORTHEAST), hit, hits, refProd, 1, 1, 0);    // NE
0163     associateNeighbour(neighbours.at(SOUTH), hit, hits, refProd, 0, -1, 0);       // S
0164     associateNeighbour(neighbours.at(SOUTHWEST), hit, hits, refProd, -1, -1, 0);  // SW
0165     associateNeighbour(neighbours.at(EAST), hit, hits, refProd, 1, 0, 0);         // E
0166     associateNeighbour(neighbours.at(SOUTHEAST), hit, hits, refProd, 1, -1, 0);   // SE
0167     associateNeighbour(neighbours.at(WEST), hit, hits, refProd, -1, 0, 0);        // W
0168     associateNeighbour(neighbours.at(NORTHWEST), hit, hits, refProd, -1, 1, 0);   // NW
0169   }
0170 
0171   bool validNeighbours(const unsigned int denseid) const {
0172     bool ok = true;
0173     unsigned index = getIdx(denseid);
0174     if (neighboursHcal_.at(index).size() != 9)
0175       ok = false;  // the neighbour vector size should be 3x3
0176     return ok;
0177   }
0178 
0179   unsigned int getIdx(const unsigned int denseid) const {
0180     unsigned index = denseid - denseIdHcalMin_;
0181     return index;
0182   }
0183 
0184 protected:
0185   edm::ESWatcher<HcalRecNumberingRecord> theRecNumberWatcher_;
0186   std::unique_ptr<const TOPO> topology_;
0187   std::vector<int> vhcalEnum_;
0188   std::vector<std::vector<DetId>> neighboursHcal_;
0189   unsigned int denseIdHcalMax_;
0190   unsigned int denseIdHcalMin_;
0191 
0192 private:
0193   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalToken_;
0194   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0195 };
0196 
0197 #endif