Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:37:49 

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::BeginRun>()),
0037         geomToken_(cc.esConsumes<edm::Transition::BeginRun>()) {}
0038 
0039   // Initialization
0040   void init(const edm::EventSetup& iSetup) override {
0041     bool check = theRecNumberWatcher_.check(iSetup);
0042     if (!check)
0043       return;
0044 
0045     edm::ESHandle<HcalTopology> hcalTopology = iSetup.getHandle(hcalToken_);
0046     topology_.release();
0047     topology_.reset(hcalTopology.product());
0048 
0049     // Fill a vector of valid denseid's
0050     edm::ESHandle<CaloGeometry> hGeom = iSetup.getHandle(geomToken_);
0051     const CaloGeometry& caloGeom = *hGeom;
0052 
0053     std::vector<DetId> vecHcal;
0054     std::vector<unsigned int> vDenseIdHcal;
0055     neighboursHcal_.clear();
0056     for (auto hcalSubdet : vhcalEnum_) {
0057       std::vector<DetId> vecDetIds(caloGeom.getValidDetIds(DetId::Hcal, hcalSubdet));
0058       vecHcal.insert(vecHcal.end(), vecDetIds.begin(), vecDetIds.end());
0059     }
0060     vDenseIdHcal.reserve(vecHcal.size());
0061     for (auto hDetId : vecHcal) {
0062       vDenseIdHcal.push_back(topology_.get()->detId2denseId(hDetId));
0063     }
0064     std::sort(vDenseIdHcal.begin(), vDenseIdHcal.end());
0065 
0066     // Fill a vector of cell neighbours
0067     denseIdHcalMax_ = *max_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
0068     denseIdHcalMin_ = *min_element(vDenseIdHcal.begin(), vDenseIdHcal.end());
0069     neighboursHcal_.resize(denseIdHcalMax_ - denseIdHcalMin_ + 1);
0070 
0071     for (auto denseid : vDenseIdHcal) {
0072       std::vector<DetId> neighbours(9, DetId(0));
0073 
0074       // the centre
0075       unsigned denseid_c = denseid;
0076       DetId detid_c = topology_.get()->denseId2detId(denseid_c);
0077       CaloNavigator<DET> navigator(detid_c, topology_.get());
0078 
0079       // Using enum in Geometry/CaloTopology/interface/CaloDirection.h
0080       // Order: CENTER(NONE),SOUTH,SOUTHEAST,SOUTHWEST,EAST,WEST,NORTHEAST,NORTHWEST,NORTH
0081       neighbours.at(NONE) = detid_c;
0082 
0083       neighbours.at(NORTH) = getNeighbourDetId(detid_c, 0);
0084       neighbours.at(SOUTH) = getNeighbourDetId(detid_c, 1);
0085       neighbours.at(EAST) = getNeighbourDetId(detid_c, 2);
0086       neighbours.at(WEST) = getNeighbourDetId(detid_c, 3);
0087 
0088       neighbours.at(NORTHEAST) = getNeighbourDetId(detid_c, 4);
0089       neighbours.at(SOUTHWEST) = getNeighbourDetId(detid_c, 5);
0090       neighbours.at(SOUTHEAST) = getNeighbourDetId(detid_c, 6);
0091       neighbours.at(NORTHWEST) = getNeighbourDetId(detid_c, 7);
0092 
0093       unsigned index = getIdx(denseid_c);
0094       neighboursHcal_[index] = neighbours;
0095 
0096     }  // loop over vDenseIdHcal
0097 
0098     if (debug) {
0099       backwardCompatibilityCheck(vDenseIdHcal);
0100       printNeighbourInfo(vDenseIdHcal);
0101     }
0102   }
0103 
0104   // Check neighbour
0105   void backwardCompatibilityCheck(const std::vector<unsigned int> vDenseIdHcal) {
0106     //
0107     // Check backward compatibility (does a neighbour of a channel have the channel as a neighbour?)
0108     //
0109     for (auto denseid : vDenseIdHcal) {
0110       DetId detid = topology_.get()->denseId2detId(denseid);
0111       HcalDetId hid = HcalDetId(detid);
0112       if (detid == DetId(0)) {
0113         edm::LogWarning("PFHCALDenseIdNavigator") << "Found an invalid DetId";
0114         continue;
0115       }
0116       if (!validNeighbours(denseid)) {
0117         edm::LogWarning("PFHCALDenseIdNavigator") << "The vector for neighbour information has an invalid length";
0118         continue;
0119       }
0120       std::vector<DetId> neighbours(9, DetId(0));
0121       unsigned index = getIdx(denseid);
0122       if (index >= neighboursHcal_.size()) {
0123         edm::LogWarning("PFHCALDenseIdNavigator") << "The vector for neighbour information is not found";
0124         continue;  // Skip if not found
0125       }
0126       neighbours = neighboursHcal_.at(index);
0127 
0128       //
0129       // Loop over neighbours
0130       int ineighbour = -1;
0131       for (auto neighbour : neighbours) {
0132         ineighbour++;
0133         if (neighbour == DetId(0))
0134           continue;
0135         HcalDetId hidn = HcalDetId(neighbour);
0136         std::vector<DetId> neighboursOfNeighbour(9, DetId(0));
0137         std::unordered_set<unsigned int> listOfNeighboursOfNeighbour;  // list of neighbours of neighbour
0138         unsigned denseidNeighbour = topology_.get()->detId2denseId(neighbour);
0139 
0140         if (!validNeighbours(denseidNeighbour)) {
0141           edm::LogWarning("PFHCALDenseIdNavigator")
0142               << "This neighbour does not have a valid neighbour information (another subdetector?). Ignore: "
0143               << detid.det() << " " << hid.ieta() << " " << hid.iphi() << " " << hid.depth() << " " << neighbour.det()
0144               << " " << hidn.ieta() << " " << hidn.iphi() << " " << hidn.depth();
0145           neighboursHcal_[index][ineighbour] = DetId(0);  // Not a valid neighbour. Set to DetId(0)
0146           continue;
0147         }
0148         if (getIdx(denseidNeighbour) >= neighboursHcal_.size())
0149           continue;
0150         neighboursOfNeighbour = neighboursHcal_.at(getIdx(denseidNeighbour));
0151 
0152         //
0153         // Loop over neighbours of neighbours
0154         for (auto neighbourOfNeighbour : neighboursOfNeighbour) {
0155           if (neighbourOfNeighbour == DetId(0))
0156             continue;
0157           unsigned denseidNeighbourOfNeighbour = topology_.get()->detId2denseId(neighbourOfNeighbour);
0158           if (!validNeighbours(denseidNeighbourOfNeighbour))
0159             continue;
0160           listOfNeighboursOfNeighbour.insert(denseidNeighbourOfNeighbour);
0161         }
0162 
0163         //
0164         if (listOfNeighboursOfNeighbour.find(denseid) == listOfNeighboursOfNeighbour.end()) {
0165           // This neighbour is not backward compatible.
0166           edm::LogWarning("PFHCALDenseIdNavigator")
0167               << "This neighbour does not have the original channel as its neighbour. Ignore: " << detid.det() << " "
0168               << hid.ieta() << " " << hid.iphi() << " " << hid.depth() << " " << neighbour.det() << " " << hidn.ieta()
0169               << " " << hidn.iphi() << " " << hidn.depth();
0170           neighboursHcal_[index][ineighbour] = DetId(0);
0171         }
0172 
0173       }  // loop over neighbours
0174     }  // loop over denseId_
0175   }
0176 
0177   // Print out neighbour DetId's
0178   void printNeighbourInfo(const std::vector<unsigned int> vDenseIdHcal) {
0179     // Print neighbour definitions
0180     for (auto denseid : vDenseIdHcal) {
0181       std::vector<DetId> neighbours(9, DetId(0));
0182 
0183       // the centre
0184       unsigned denseid_c = denseid;
0185       DetId detid_c = topology_.get()->denseId2detId(denseid_c);
0186       CaloNavigator<DET> navigator(detid_c, topology_.get());
0187 
0188       unsigned index = getIdx(denseid_c);
0189 
0190       neighbours = neighboursHcal_[index];
0191 
0192       const DetId N = neighbours.at(NORTH);
0193       const DetId S = neighbours.at(SOUTH);
0194       const DetId E = neighbours.at(EAST);
0195       const DetId W = neighbours.at(WEST);
0196 
0197       const DetId NE = neighbours.at(NORTHEAST);
0198       const DetId SW = neighbours.at(SOUTHWEST);
0199       const DetId SE = neighbours.at(SOUTHEAST);
0200       const DetId NW = neighbours.at(NORTHWEST);
0201 
0202       edm::LogPrint("PFHCALDenseIdNavigator")
0203           << "PFHCALDenseIdNavigator: " << HcalDetId(detid_c).ieta() << " " << HcalDetId(detid_c).iphi() << " "
0204           << HcalDetId(detid_c).depth() << " " << HcalDetId(detid_c).subdetId() << ": " << HcalDetId(N).ieta() << " "
0205           << HcalDetId(N).iphi() << " " << HcalDetId(N).depth() << " " << HcalDetId(N).subdetId() << ", "
0206           << HcalDetId(E).ieta() << " " << HcalDetId(E).iphi() << " " << HcalDetId(E).depth() << " "
0207           << HcalDetId(E).subdetId() << ", " << HcalDetId(S).ieta() << " " << HcalDetId(S).iphi() << " "
0208           << HcalDetId(S).depth() << " " << HcalDetId(S).subdetId() << ", " << HcalDetId(W).ieta() << " "
0209           << HcalDetId(W).iphi() << " " << HcalDetId(W).depth() << " " << HcalDetId(W).subdetId() << ", "
0210           << HcalDetId(NE).ieta() << " " << HcalDetId(NE).iphi() << " " << HcalDetId(NE).depth() << " "
0211           << HcalDetId(NE).subdetId() << ", " << HcalDetId(SW).ieta() << " " << HcalDetId(SW).iphi() << " "
0212           << HcalDetId(SW).depth() << " " << HcalDetId(SW).subdetId() << ", " << HcalDetId(SE).ieta() << " "
0213           << HcalDetId(SE).iphi() << " " << HcalDetId(SE).depth() << " " << HcalDetId(SE).subdetId() << ", "
0214           << HcalDetId(NW).ieta() << " " << HcalDetId(NW).iphi() << " " << HcalDetId(NW).depth() << " "
0215           << HcalDetId(NW).subdetId();
0216 
0217     }  // print ends
0218   }
0219 
0220   // Check if two DetId's have the same subdetId
0221   const bool checkSameSubDet(const DetId detId, const DetId detId2) {
0222     HcalDetId hid = HcalDetId(detId);
0223     HcalDetId hid2 = HcalDetId(detId2);
0224     return hid.subdetId() == hid2.subdetId();
0225   }
0226 
0227   //https://cmssdt.cern.ch/lxr/source/DataFormats/HcalDetId/interface/HcalDetId.h#0141
0228   static constexpr int getZside(const DetId detId) { return ((detId & HcalDetId::kHcalZsideMask2) ? (1) : (-1)); }
0229 
0230   // Obtain the neighbour's DetId
0231   const uint32_t getNeighbourDetId(const DetId detId, const uint32_t direction) {
0232     // desired order for PF: NORTH, SOUTH, EAST, WEST, NORTHEAST, SOUTHWEST, SOUTHEAST, NORTHWEST
0233     if (detId == DetId(0))
0234       return 0;
0235 
0236     if (direction == 0)                        // NORTH
0237       return topology_.get()->goNorth(detId);  // larger iphi values (except phi boundary)
0238 
0239     if (direction == 1)                        // SOUTH
0240       return topology_.get()->goSouth(detId);  // smaller iphi values (except phi boundary)
0241 
0242     // In the case of East/West, make sure we are not moving to another subdetector
0243     if (direction == 2 && checkSameSubDet(detId, topology_.get()->goEast(detId)))  // EAST
0244       return topology_.get()->goEast(detId);                                       // smaller ieta values
0245 
0246     if (direction == 3 && checkSameSubDet(detId, topology_.get()->goWest(detId)))  // WEST
0247       return topology_.get()->goWest(detId);                                       // larger ieta values
0248 
0249     // In the case of cornor cells, ensure backward compatibility.
0250     // Also make sure it's not already defined as E(ast) or W(est)
0251     if (direction == 4) {         // NORTHEAST
0252       if (getZside(detId) > 0) {  // positive eta: east -> move to smaller |ieta| (finner phi granularity) first
0253         const DetId E = getNeighbourDetId(detId, 2);
0254         const DetId NE = getNeighbourDetId(E, 0);
0255         if (getNeighbourDetId(NE, 1) == E && NE != E)
0256           return NE;  //
0257       } else {        // negative eta: move in phi first then move to east (coarser phi granularity)
0258         const DetId N = getNeighbourDetId(detId, 0);
0259         const DetId NE = getNeighbourDetId(N, 2);
0260         const DetId E = getNeighbourDetId(detId, 2);
0261         if (getNeighbourDetId(NE, 3) == N && NE != E)
0262           return NE;
0263       }
0264     }
0265     if (direction == 5) {         // SOUTHWEST
0266       if (getZside(detId) > 0) {  // positive eta: move in phi first then move to west (coarser phi granularity)
0267         const DetId S = getNeighbourDetId(detId, 1);
0268         const DetId SW = getNeighbourDetId(S, 3);
0269         const DetId W = getNeighbourDetId(detId, 3);
0270         if (getNeighbourDetId(SW, 2) == S && SW != W)
0271           return SW;
0272       } else {  // negative eta: west -> move to smaller |ieta| (finner phi granularity) first
0273         const DetId W = getNeighbourDetId(detId, 3);
0274         const DetId SW = getNeighbourDetId(W, 1);
0275         if (getNeighbourDetId(SW, 0) == W && SW != W)
0276           return SW;
0277       }
0278     }
0279     if (direction == 6) {         // SOUTHEAST
0280       if (getZside(detId) > 0) {  // positive eta: east -> move to smaller |ieta| (finner phi granularity) first
0281         const DetId E = getNeighbourDetId(detId, 2);
0282         const DetId SE = getNeighbourDetId(E, 1);
0283         if (getNeighbourDetId(SE, 0) == E && SE != E)
0284           return SE;
0285       } else {  // negative eta: move in phi first then move to east (coarser phi granularity)
0286         const DetId S = getNeighbourDetId(detId, 1);
0287         const DetId SE = getNeighbourDetId(S, 2);
0288         const DetId E = getNeighbourDetId(detId, 2);
0289         if (getNeighbourDetId(SE, 3) == S && SE != E)
0290           return SE;
0291       }
0292     }
0293     if (direction == 7) {         // NORTHWEST
0294       if (getZside(detId) > 0) {  // positive eta: move in phi first then move to west (coarser phi granularity)
0295         const DetId N = getNeighbourDetId(detId, 0);
0296         const DetId NW = getNeighbourDetId(N, 3);
0297         const DetId W = getNeighbourDetId(detId, 3);
0298         if (getNeighbourDetId(NW, 2) == N && NW != W)
0299           return NW;
0300       } else {  // negative eta: west -> move to smaller |ieta| (finner phi granularity) first
0301         const DetId W = getNeighbourDetId(detId, 3);
0302         const DetId NW = getNeighbourDetId(W, 0);
0303         if (getNeighbourDetId(NW, 1) == W && NW != W)
0304           return NW;
0305       }
0306     }
0307     return 0;
0308   }
0309 
0310   // Associate neighbour PFRecHits
0311   void associateNeighbours(reco::PFRecHit& hit,
0312                            std::unique_ptr<reco::PFRecHitCollection>& hits,
0313                            edm::RefProd<reco::PFRecHitCollection>& refProd) override {
0314     DetId detid(hit.detId());
0315     HcalDetId hid(detid);
0316     unsigned denseid = topology_.get()->detId2denseId(detid);
0317 
0318     std::vector<DetId> neighbours(9, DetId(0));
0319 
0320     if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_) {
0321       edm::LogWarning("PFRecHitHCALCachedNavigator")
0322           << " DenseId for this cell is out of the expected range." << std::endl;
0323     } else if (!validNeighbours(denseid)) {
0324       edm::LogWarning("PFRecHitHCALCachedNavigator")
0325           << " DenseId for this cell does not have the neighbour information. " << hid.ieta() << " " << hid.iphi()
0326           << " " << hid.depth() << " " << hid.subdetId();
0327     } else {
0328       unsigned index = getIdx(denseid);
0329       neighbours = neighboursHcal_.at(index);
0330     }
0331 
0332     associateNeighbour(neighbours.at(NORTH), hit, hits, refProd, 0, 1, 0);        // N
0333     associateNeighbour(neighbours.at(NORTHEAST), hit, hits, refProd, 1, 1, 0);    // NE
0334     associateNeighbour(neighbours.at(SOUTH), hit, hits, refProd, 0, -1, 0);       // S
0335     associateNeighbour(neighbours.at(SOUTHWEST), hit, hits, refProd, -1, -1, 0);  // SW
0336     associateNeighbour(neighbours.at(EAST), hit, hits, refProd, 1, 0, 0);         // E
0337     associateNeighbour(neighbours.at(SOUTHEAST), hit, hits, refProd, 1, -1, 0);   // SE
0338     associateNeighbour(neighbours.at(WEST), hit, hits, refProd, -1, 0, 0);        // W
0339     associateNeighbour(neighbours.at(NORTHWEST), hit, hits, refProd, -1, 1, 0);   // NW
0340   }
0341 
0342   // Check if we can get the valid neighbour vector for a given denseid
0343   const bool validNeighbours(const unsigned int denseid) const {
0344     if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_)
0345       return false;  // neighbour denseid is out of the expected range
0346     unsigned index = getIdx(denseid);
0347     if (neighboursHcal_.at(index).size() != 9)
0348       return false;  // the neighbour vector size should be 3x3
0349     return true;
0350   }
0351 
0352   // Get the index of neighboursHcal_ for a given denseid
0353   const unsigned int getIdx(const unsigned int denseid) const {
0354     if (denseid < denseIdHcalMin_ || denseid > denseIdHcalMax_)
0355       return unsigned(denseIdHcalMax_ - denseIdHcalMin_ +
0356                       1);  // out-of-bounce (different subdetector groups. give a dummy number.)
0357     unsigned index = denseid - denseIdHcalMin_;
0358     return index;
0359   }
0360 
0361 protected:
0362   edm::ESWatcher<HcalRecNumberingRecord> theRecNumberWatcher_;
0363   std::unique_ptr<const TOPO> topology_;
0364   std::vector<int> vhcalEnum_;
0365   std::vector<std::vector<DetId>> neighboursHcal_;
0366   unsigned int denseIdHcalMax_;
0367   unsigned int denseIdHcalMin_;
0368 
0369 private:
0370   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalToken_;
0371   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0372   const bool debug = false;
0373 };
0374 
0375 #endif
[ 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!