Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-12-20 03:22:20

0001 #include <algorithm>
0002 #include <memory>
0003 #include <type_traits>
0004 #include <utility>
0005 #include <variant>
0006 
0007 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0008 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0009 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0010 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0013 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0014 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0015 #include "Geometry/CaloTopology/interface/EcalBarrelTopology.h"
0016 #include "Geometry/CaloTopology/interface/EcalEndcapTopology.h"
0017 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0018 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0019 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0020 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ESProducer.h"
0021 #include "CalorimeterDefinitions.h"
0022 
0023 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0024   using namespace particleFlowRecHitProducer;
0025 
0026   template <typename CAL>
0027   class PFRecHitTopologyESProducer : public ESProducer {
0028   public:
0029     PFRecHitTopologyESProducer(edm::ParameterSet const& iConfig)
0030         : ESProducer(iConfig), cutsFromDB_(iConfig.getParameter<bool>("usePFThresholdsFromDB")) {
0031       auto cc = setWhatProduced(this);
0032       geomToken_ = cc.consumes();
0033       if constexpr (std::is_same_v<CAL, HCAL>) {
0034         hcalToken_ = cc.consumes();
0035         hcalCutsToken_ = cc.consumes();
0036       }
0037     }
0038 
0039     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0040       edm::ParameterSetDescription desc;
0041       if constexpr (std::is_same_v<CAL, HCAL>)
0042         desc.add<bool>("usePFThresholdsFromDB", true);
0043       else  // only needs to be true for HBHE
0044         desc.add<bool>("usePFThresholdsFromDB", false);
0045       descriptions.addWithDefaultLabel(desc);
0046     }
0047 
0048     std::unique_ptr<typename CAL::TopologyTypeHost> produce(const typename CAL::TopologyRecordType& iRecord) {
0049       const auto& geom = iRecord.get(geomToken_);
0050       auto product = std::make_unique<typename CAL::TopologyTypeHost>(CAL::kSize, cms::alpakatools::host());
0051       auto view = product->view();
0052       const int calEnums[2] = {CAL::kSubdetectorBarrelId, CAL::kSubdetectorEndcapId};
0053       for (const auto subdet : calEnums) {
0054         // Construct topology
0055         //  for HCAL: using dedicated record
0056         //  for ECAL: from CaloGeometry (separate for barrel and endcap)
0057         const CaloSubdetectorGeometry* geo = geom.getSubdetectorGeometry(CAL::kDetectorId, subdet);
0058         const CaloSubdetectorTopology* topo;
0059         std::variant<EcalBarrelTopology, EcalEndcapTopology> topoVar;  // need to store ECAL topology temporarily
0060         if constexpr (std::is_same_v<CAL, HCAL>)
0061           topo = &iRecord.get(hcalToken_);
0062         else if (subdet == EcalSubdetector::EcalBarrel)
0063           topo = &topoVar.emplace<EcalBarrelTopology>(geom);
0064         else
0065           topo = &topoVar.emplace<EcalEndcapTopology>(geom);
0066 
0067         // Fill product
0068         for (auto const detId : geom.getValidDetIds(CAL::kDetectorId, subdet)) {
0069           const uint32_t denseId = CAL::detId2denseId(detId);
0070           assert(denseId < CAL::kSize);
0071 
0072           // Fill SoA members with HCAL PF Thresholds from GT
0073           if constexpr (std::is_same_v<CAL, HCAL>) {
0074             view.cutsFromDB() = false;
0075             if (cutsFromDB_) {
0076               view.cutsFromDB() = true;
0077               const HcalPFCuts& pfCuts = iRecord.get(hcalCutsToken_);
0078               const HcalTopology& htopo = iRecord.get(hcalToken_);
0079               std::unique_ptr<HcalPFCuts> prod = std::make_unique<HcalPFCuts>(pfCuts);
0080               prod->setTopo(&htopo);
0081               view.noiseThreshold(denseId) = prod->getValues(detId.rawId())->noiseThreshold();
0082               view.seedThreshold(denseId) = prod->getValues(detId.rawId())->seedThreshold();
0083             }
0084           }
0085 
0086           const GlobalPoint pos = geo->getGeometry(detId)->getPosition();
0087           view.positionX(denseId) = pos.x();
0088           view.positionY(denseId) = pos.y();
0089           view.positionZ(denseId) = pos.z();
0090 
0091           for (uint32_t n = 0; n < 8; n++) {
0092             uint32_t neighDetId = getNeighbourDetId(detId, n, *topo);
0093             if (CAL::detIdInRange(neighDetId))
0094               view.neighbours(denseId)(n) = CAL::detId2denseId(neighDetId);
0095             else
0096               view.neighbours(denseId)(n) = 0xffffffff;
0097           }
0098         }
0099       }
0100 
0101       // Remove neighbours that are not backward compatible (only for HCAL)
0102       if constexpr (std::is_same_v<CAL, HCAL>) {
0103         for (const auto subdet : calEnums)
0104           for (auto const detId : geom.getValidDetIds(CAL::kDetectorId, subdet)) {
0105             const uint32_t denseId = CAL::detId2denseId(detId);
0106             for (uint32_t n = 0; n < 8; n++) {
0107               const ::reco::PFRecHitsTopologyNeighbours& neighboursOfNeighbour =
0108                   view.neighbours(view.neighbours(denseId)[n]);
0109               if (std::find(neighboursOfNeighbour.begin(), neighboursOfNeighbour.end(), denseId) ==
0110                   neighboursOfNeighbour.end())
0111                 view.neighbours(denseId)[n] = 0xffffffff;
0112             }
0113           }
0114       }
0115 
0116       // Print results (for debugging)
0117       LogDebug("PFRecHitTopologyESProducer").log([&](auto& log) {
0118         for (const auto subdet : calEnums)
0119           for (const auto detId : geom.getValidDetIds(CAL::kDetectorId, subdet)) {
0120             const uint32_t denseId = CAL::detId2denseId(detId);
0121             log.format("detId:{} denseId:{} pos:{},{},{} neighbours:{},{},{},{};{},{},{},{}\n",
0122                        (uint32_t)detId,
0123                        denseId,
0124                        view[denseId].positionX(),
0125                        view[denseId].positionY(),
0126                        view[denseId].positionZ(),
0127                        view[denseId].neighbours()(0),
0128                        view[denseId].neighbours()(1),
0129                        view[denseId].neighbours()(2),
0130                        view[denseId].neighbours()(3),
0131                        view[denseId].neighbours()(4),
0132                        view[denseId].neighbours()(5),
0133                        view[denseId].neighbours()(6),
0134                        view[denseId].neighbours()(7));
0135           }
0136       });
0137 
0138       return product;
0139     }
0140 
0141   private:
0142     edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geomToken_;
0143     edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> hcalToken_;
0144     edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0145     const bool cutsFromDB_;
0146 
0147     // specialised for HCAL/ECAL, because non-nearest neighbours are defined differently
0148     uint32_t getNeighbourDetId(const uint32_t detId, const uint32_t direction, const CaloSubdetectorTopology& topo);
0149   };
0150 
0151   template <>
0152   uint32_t PFRecHitTopologyESProducer<ECAL>::getNeighbourDetId(const uint32_t detId,
0153                                                                const uint32_t direction,
0154                                                                const CaloSubdetectorTopology& topo) {
0155     // desired order for PF: NORTH, SOUTH, EAST, WEST, NORTHEAST, SOUTHWEST, SOUTHEAST, NORTHWEST
0156     if (detId == 0)
0157       return 0;
0158 
0159     if (direction == 0)            // NORTH
0160       return topo.goNorth(detId);  // larger iphi values (except phi boundary)
0161     if (direction == 1)            // SOUTH
0162       return topo.goSouth(detId);  // smaller iphi values (except phi boundary)
0163     if (direction == 2)            // EAST
0164       return topo.goEast(detId);   // smaller ieta values
0165     if (direction == 3)            // WEST
0166       return topo.goWest(detId);   // larger ieta values
0167 
0168     if (direction == 4) {  // NORTHEAST
0169       const uint32_t NE = getNeighbourDetId(getNeighbourDetId(detId, 0, topo), 2, topo);
0170       if (NE)
0171         return NE;
0172       return getNeighbourDetId(getNeighbourDetId(detId, 2, topo), 0, topo);
0173     }
0174     if (direction == 5) {  // SOUTHWEST
0175       const uint32_t SW = getNeighbourDetId(getNeighbourDetId(detId, 1, topo), 3, topo);
0176       if (SW)
0177         return SW;
0178       return getNeighbourDetId(getNeighbourDetId(detId, 3, topo), 1, topo);
0179     }
0180     if (direction == 6) {  // SOUTHEAST
0181       const uint32_t ES = getNeighbourDetId(getNeighbourDetId(detId, 2, topo), 1, topo);
0182       if (ES)
0183         return ES;
0184       return getNeighbourDetId(getNeighbourDetId(detId, 1, topo), 2, topo);
0185     }
0186     if (direction == 7) {  // NORTHWEST
0187       const uint32_t WN = getNeighbourDetId(getNeighbourDetId(detId, 3, topo), 0, topo);
0188       if (WN)
0189         return WN;
0190       return getNeighbourDetId(getNeighbourDetId(detId, 0, topo), 3, topo);
0191     }
0192     return 0;
0193   }
0194 
0195   template <>
0196   uint32_t PFRecHitTopologyESProducer<HCAL>::getNeighbourDetId(const uint32_t detId,
0197                                                                const uint32_t direction,
0198                                                                const CaloSubdetectorTopology& topo) {
0199     // desired order for PF: NORTH, SOUTH, EAST, WEST, NORTHEAST, SOUTHWEST, SOUTHEAST, NORTHWEST
0200     if (detId == 0)
0201       return 0;
0202 
0203     if (direction == 0)            // NORTH
0204       return topo.goNorth(detId);  // larger iphi values (except phi boundary)
0205     if (direction == 1)            // SOUTH
0206       return topo.goSouth(detId);  // smaller iphi values (except phi boundary)
0207     if (direction == 2)            // EAST
0208       return topo.goEast(detId);   // smaller ieta values
0209     if (direction == 3)            // WEST
0210       return topo.goWest(detId);   // larger ieta values
0211 
0212     std::pair<uint32_t, uint32_t> directions;
0213     if (direction == 4) {  // NORTHEAST
0214       if (HCAL::getZside(detId) > 0)
0215         directions = {2, 0};  // positive eta: east -> move to smaller |ieta| (finner phi granularity) first
0216       else
0217         directions = {0, 2};      // negative eta: move in phi first then move to east (coarser phi granularity)
0218     } else if (direction == 5) {  // SOUTHWEST
0219       if (HCAL::getZside(detId) > 0)
0220         directions = {1, 3};  // positive eta: move in phi first then move to west (coarser phi granularity)
0221       else
0222         directions = {3, 1};      // negative eta: west -> move to smaller |ieta| (finner phi granularity) first
0223     } else if (direction == 6) {  // SOUTHEAST
0224       if (HCAL::getZside(detId) > 0)
0225         directions = {2, 1};  // positive eta: east -> move to smaller |ieta| (finner phi granularity) first
0226       else
0227         directions = {1, 2};      // negative eta: move in phi first then move to east (coarser phi granularity)
0228     } else if (direction == 7) {  // NORTHWEST
0229       if (HCAL::getZside(detId) > 0)
0230         directions = {0, 3};  // positive eta: move in phi first then move to west (coarser phi granularity)
0231       else
0232         directions = {3, 0};  // negative eta: west -> move to smaller |ieta| (finner phi granularity) first
0233     } else
0234       return 0;
0235     const uint32_t nn1 = getNeighbourDetId(detId, directions.first, topo);   // nearest neighbour in direction 1
0236     const uint32_t nn2 = getNeighbourDetId(detId, directions.second, topo);  // nearest neighbour in direction 2
0237     const uint32_t nnn = getNeighbourDetId(nn1, directions.second, topo);    // next-to-nearest neighbour
0238     if (nnn == nn1 || nnn == nn2)                                            // avoid duplicates
0239       return 0;
0240     return nnn;
0241   }
0242 
0243   using PFRecHitECALTopologyESProducer = PFRecHitTopologyESProducer<ECAL>;
0244   using PFRecHitHCALTopologyESProducer = PFRecHitTopologyESProducer<HCAL>;
0245 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
0246 
0247 #include "HeterogeneousCore/AlpakaCore/interface/alpaka/ModuleFactory.h"
0248 DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(PFRecHitECALTopologyESProducer);
0249 DEFINE_FWK_EVENTSETUP_ALPAKA_MODULE(PFRecHitHCALTopologyESProducer);