plugins/alpaka/PFRecHitProducerKernel.dev.cc

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