Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-06-18 02:20:40 

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