Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-01-22 07:34:27 

0001 #include <type_traits>
0002 
0003 #include <alpaka/alpaka.hpp>
0004 
0005 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0006 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0007 
0008 #include "PFRecHitProducerKernel.h"
0009 
0010 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0011   using namespace particleFlowRecHitProducer;
0012 
0013   // Kernel to apply cuts to calorimeter hits and construct PFRecHits
0014   template <typename CAL>
0015   struct PFRecHitProducerKernelConstruct {
0016     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0017                                   const typename CAL::ParameterType::ConstView params,
0018                                   const typename CAL::TopologyTypeDevice::ConstView topology,
0019                                   const typename CAL::CaloRecHitSoATypeDevice::ConstView recHits,
0020                                   reco::PFRecHitDeviceCollection::View pfRecHits,
0021                                   uint32_t* __restrict__ denseId2pfRecHit,
0022                                   uint32_t* __restrict__ num_pfRecHits) const {
0023       // Strided loop over CaloRecHits
0024       for (int32_t i : cms::alpakatools::uniform_elements(acc, recHits.metadata().size())) {
0025         // Check energy thresholds/quality cuts (specialised for HCAL/ECAL)
0026         if (!applyCuts(recHits[i], params, topology))
0027           continue;
0028 
0029         // Use atomic operation to determine index of the PFRecHit to be constructed
0030         // The index needs to be unique and consequtive across all threads in all blocks.
0031         // This is achieved using the alpaka::hierarchy::Blocks argument.
0032         const uint32_t j = alpaka::atomicInc(acc, num_pfRecHits, 0xffffffff, alpaka::hierarchy::Blocks{});
0033 
0034         // Construct PFRecHit from CAL recHit (specialised for HCAL/ECAL)
0035         constructPFRecHit(pfRecHits[j], recHits[i]);
0036 
0037         // Fill denseId -> pfRecHit index map
0038         denseId2pfRecHit[CAL::detId2denseId(pfRecHits.detId(j))] = j;
0039       }
0040     }
0041 
0042     ALPAKA_FN_ACC static bool applyCuts(const typename CAL::CaloRecHitSoATypeDevice::ConstView::const_element rh,
0043                                         const typename CAL::ParameterType::ConstView params,
0044                                         const typename CAL::TopologyTypeDevice::ConstView topology);
0045 
0046     ALPAKA_FN_ACC static void constructPFRecHit(
0047         reco::PFRecHitDeviceCollection::View::element pfrh,
0048         const typename CAL::CaloRecHitSoATypeDevice::ConstView::const_element rh);
0049   };
0050 
0051   template <>
0052   ALPAKA_FN_ACC bool PFRecHitProducerKernelConstruct<HCAL>::applyCuts(
0053       const typename HCAL::CaloRecHitSoATypeDevice::ConstView::const_element rh,
0054       const HCAL::ParameterType::ConstView params,
0055       const HCAL::TopologyTypeDevice::ConstView topology) {
0056     // Reject HCAL recHits below enery threshold
0057     float threshold = 9999.;
0058     const uint32_t detId = rh.detId();
0059     const uint32_t depth = HCAL::getDepth(detId);
0060     const uint32_t subdet = getSubdet(detId);
0061 
0062     // skip bad channels
0063     if (rh.chi2() < 0)
0064       return false;
0065 
0066     if (topology.cutsFromDB()) {
0067       const auto& denseId = HCAL::detId2denseId(detId);
0068       if (denseId != HCAL::kInvalidDenseId) {
0069         threshold = topology.noiseThreshold()[denseId];
0070       } else {
0071         printf("Encountered invalid denseId for detId %u (subdetector %u)!", detId, subdet);
0072         return false;
0073       }
0074     } else {
0075       if (subdet == HcalBarrel) {
0076         threshold = params.energyThresholds()[depth - 1];
0077       } else if (subdet == HcalEndcap) {
0078         threshold = params.energyThresholds()[depth - 1 + HCAL::kMaxDepthHB];
0079       } else {
0080         printf("Rechit with detId %u has invalid subdetector %u!\n", detId, subdet);
0081         return false;
0082       }
0083     }
0084     return rh.energy() >= threshold;
0085   }
0086 
0087   template <>
0088   ALPAKA_FN_ACC bool PFRecHitProducerKernelConstruct<ECAL>::applyCuts(
0089       const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh,
0090       const ECAL::ParameterType::ConstView params,
0091       const ECAL::TopologyTypeDevice::ConstView topology) {
0092     // Reject ECAL recHits below energy threshold
0093     if (rh.energy() < params.energyThresholds()[ECAL::detId2denseId(rh.detId())])
0094       return false;
0095 
0096     // Reject ECAL recHits of bad quality
0097     if ((ECAL::checkFlag(rh.flags(), ECAL::Flags::kOutOfTime) && rh.energy() > params.cleaningThreshold()) ||
0098         ECAL::checkFlag(rh.flags(), ECAL::Flags::kTowerRecovered) || ECAL::checkFlag(rh.flags(), ECAL::Flags::kWeird) ||
0099         ECAL::checkFlag(rh.flags(), ECAL::Flags::kDiWeird))
0100       return false;
0101 
0102     return true;
0103   }
0104 
0105   template <>
0106   ALPAKA_FN_ACC void PFRecHitProducerKernelConstruct<HCAL>::constructPFRecHit(
0107       reco::PFRecHitDeviceCollection::View::element pfrh,
0108       const HCAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) {
0109     pfrh.detId() = rh.detId();
0110     pfrh.denseId() = HCAL::detId2denseId(rh.detId());
0111     pfrh.energy() = rh.energy();
0112     pfrh.time() = rh.timeM0();
0113     pfrh.depth() = HCAL::getDepth(pfrh.detId());
0114     const uint32_t subdet = getSubdet(pfrh.detId());
0115     if (subdet == HcalBarrel)
0116       pfrh.layer() = PFLayer::HCAL_BARREL1;
0117     else if (subdet == HcalEndcap)
0118       pfrh.layer() = PFLayer::HCAL_ENDCAP;
0119     else
0120       pfrh.layer() = PFLayer::NONE;
0121   }
0122 
0123   template <>
0124   ALPAKA_FN_ACC void PFRecHitProducerKernelConstruct<ECAL>::constructPFRecHit(
0125       reco::PFRecHitDeviceCollection::View::element pfrh,
0126       const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) {
0127     pfrh.detId() = rh.detId();
0128     pfrh.denseId() = ECAL::detId2denseId(rh.detId());
0129     pfrh.energy() = rh.energy();
0130     pfrh.time() = rh.time();
0131     pfrh.depth() = 1;
0132     const uint32_t subdet = getSubdet(pfrh.detId());
0133     if (subdet == EcalBarrel)
0134       pfrh.layer() = PFLayer::ECAL_BARREL;
0135     else if (subdet == EcalEndcap)
0136       pfrh.layer() = PFLayer::ECAL_ENDCAP;
0137     else
0138       pfrh.layer() = PFLayer::NONE;
0139   }
0140 
0141   // Kernel to associate topology information of PFRecHits
0142   template <typename CAL>
0143   struct PFRecHitProducerKernelTopology {
0144     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0145                                   const typename CAL::TopologyTypeDevice::ConstView topology,
0146                                   reco::PFRecHitDeviceCollection::View pfRecHits,
0147                                   const uint32_t* __restrict__ denseId2pfRecHit,
0148                                   uint32_t* __restrict__ num_pfRecHits) const {
0149       // First thread updates size field pfRecHits SoA
0150       if (const int32_t thread = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0u]; thread == 0)
0151         pfRecHits.size() = *num_pfRecHits;
0152 
0153       // Assign position information and associate neighbours
0154       for (int32_t i : cms::alpakatools::uniform_elements(acc, *num_pfRecHits)) {
0155         const uint32_t denseId = CAL::detId2denseId(pfRecHits.detId(i));
0156 
0157         pfRecHits.x(i) = topology.positionX(denseId);
0158         pfRecHits.y(i) = topology.positionY(denseId);
0159         pfRecHits.z(i) = topology.positionZ(denseId);
0160 
0161         for (uint32_t n = 0; n < 8; n++) {
0162           pfRecHits.neighbours(i)(n) = -1;
0163           const uint32_t denseId_neighbour = topology.neighbours(denseId)(n);
0164           if (denseId_neighbour != CAL::kInvalidDenseId) {
0165             const uint32_t pfRecHit_neighbour = denseId2pfRecHit[denseId_neighbour];
0166             if (pfRecHit_neighbour != 0xffffffff)
0167               pfRecHits.neighbours(i)(n) = (int32_t)pfRecHit_neighbour;
0168           }
0169         }
0170       }
0171     }
0172   };
0173 
0174   template <typename CAL>
0175   PFRecHitProducerKernel<CAL>::PFRecHitProducerKernel(Queue& queue, const uint32_t num_recHits)
0176       : denseId2pfRecHit_(cms::alpakatools::make_device_buffer<uint32_t[]>(queue, CAL::kSize)),
0177         num_pfRecHits_(cms::alpakatools::make_device_buffer<uint32_t>(queue)),
0178         work_div_(cms::alpakatools::make_workdiv<Acc1D>(1, 1)) {
0179     alpaka::memset(queue, denseId2pfRecHit_, 0xff);  // Reset denseId -> pfRecHit index map
0180     alpaka::memset(queue, num_pfRecHits_, 0x00);     // Reset global pfRecHit counter
0181 
0182     const uint32_t items = 64;
0183     const uint32_t groups = cms::alpakatools::divide_up_by(num_recHits, items);
0184     work_div_ = cms::alpakatools::make_workdiv<Acc1D>(groups, items);
0185   }
0186 
0187   template <typename CAL>
0188   void PFRecHitProducerKernel<CAL>::processRecHits(Queue& queue,
0189                                                    const typename CAL::CaloRecHitSoATypeDevice& recHits,
0190                                                    const typename CAL::ParameterType& params,
0191                                                    const typename CAL::TopologyTypeDevice& topology,
0192                                                    reco::PFRecHitDeviceCollection& pfRecHits) {
0193     alpaka::exec<Acc1D>(queue,
0194                         work_div_,
0195                         PFRecHitProducerKernelConstruct<CAL>{},
0196                         params.view(),
0197                         topology.view(),
0198                         recHits.view(),
0199                         pfRecHits.view(),
0200                         denseId2pfRecHit_.data(),
0201                         num_pfRecHits_.data());
0202   }
0203 
0204   template <typename CAL>
0205   void PFRecHitProducerKernel<CAL>::associateTopologyInfo(Queue& queue,
0206                                                           const typename CAL::TopologyTypeDevice& topology,
0207                                                           reco::PFRecHitDeviceCollection& pfRecHits) {
0208     alpaka::exec<Acc1D>(queue,
0209                         work_div_,
0210                         PFRecHitProducerKernelTopology<CAL>{},
0211                         topology.view(),
0212                         pfRecHits.view(),
0213                         denseId2pfRecHit_.data(),
0214                         num_pfRecHits_.data());
0215   }
0216 
0217   // Instantiate templates
0218   template class PFRecHitProducerKernel<HCAL>;
0219   template class PFRecHitProducerKernel<ECAL>;
0220 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
[ 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!