Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-07-31 22:58:41 

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 
0063     // skip bad channels
0064     if (rh.chi2() < 0)
0065       return false;
0066 
0067     if (topology.cutsFromDB()) {
0068       const auto& denseId = HCAL::detId2denseId(detId);
0069       if (denseId != HCAL::kInvalidDenseId) {
0070         threshold = topology.noiseThreshold()[denseId];
0071       } else {
0072         printf("Encountered invalid denseId for detId %u (subdetector %u)!", detId, subdet);
0073         return false;
0074       }
0075     } else {
0076       if (subdet == HcalBarrel) {
0077         threshold = params.energyThresholds()[depth - 1];
0078       } else if (subdet == HcalEndcap) {
0079         threshold = params.energyThresholds()[depth - 1 + HCAL::kMaxDepthHB];
0080       } else {
0081         printf("Rechit with detId %u has invalid subdetector %u!\n", detId, subdet);
0082         return false;
0083       }
0084     }
0085     return rh.energy() >= threshold;
0086   }
0087 
0088   template <>
0089   ALPAKA_FN_ACC bool PFRecHitProducerKernelConstruct<ECAL>::applyCuts(
0090       const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh,
0091       const ECAL::ParameterType::ConstView params,
0092       const ECAL::TopologyTypeDevice::ConstView topology) {
0093     // Reject ECAL recHits below energy threshold
0094     if (rh.energy() < params.energyThresholds()[ECAL::detId2denseId(rh.detId())])
0095       return false;
0096 
0097     // Reject ECAL recHits of bad quality
0098     if ((ECAL::checkFlag(rh.flags(), ECAL::Flags::kOutOfTime) && rh.energy() > params.cleaningThreshold()) ||
0099         ECAL::checkFlag(rh.flags(), ECAL::Flags::kTowerRecovered) || ECAL::checkFlag(rh.flags(), ECAL::Flags::kWeird) ||
0100         ECAL::checkFlag(rh.flags(), ECAL::Flags::kDiWeird))
0101       return false;
0102 
0103     return true;
0104   }
0105 
0106   template <>
0107   ALPAKA_FN_ACC void PFRecHitProducerKernelConstruct<HCAL>::constructPFRecHit(
0108       reco::PFRecHitDeviceCollection::View::element pfrh,
0109       const HCAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) {
0110     pfrh.detId() = rh.detId();
0111     pfrh.denseId() = HCAL::detId2denseId(rh.detId());
0112     pfrh.energy() = rh.energy();
0113     pfrh.time() = rh.timeM0();
0114     pfrh.depth() = HCAL::getDepth(pfrh.detId());
0115     const uint32_t subdet = getSubdet(pfrh.detId());
0116     if (subdet == HcalBarrel)
0117       pfrh.layer() = PFLayer::HCAL_BARREL1;
0118     else if (subdet == HcalEndcap)
0119       pfrh.layer() = PFLayer::HCAL_ENDCAP;
0120     else
0121       pfrh.layer() = PFLayer::NONE;
0122   }
0123 
0124   template <>
0125   ALPAKA_FN_ACC void PFRecHitProducerKernelConstruct<ECAL>::constructPFRecHit(
0126       reco::PFRecHitDeviceCollection::View::element pfrh,
0127       const ECAL::CaloRecHitSoATypeDevice::ConstView::const_element rh) {
0128     pfrh.detId() = rh.detId();
0129     pfrh.denseId() = ECAL::detId2denseId(rh.detId());
0130     pfrh.energy() = rh.energy();
0131     pfrh.time() = rh.time();
0132     pfrh.depth() = 1;
0133     const uint32_t subdet = getSubdet(pfrh.detId());
0134     if (subdet == EcalBarrel)
0135       pfrh.layer() = PFLayer::ECAL_BARREL;
0136     else if (subdet == EcalEndcap)
0137       pfrh.layer() = PFLayer::ECAL_ENDCAP;
0138     else
0139       pfrh.layer() = PFLayer::NONE;
0140   }
0141 
0142   // Kernel to associate topology information of PFRecHits
0143   template <typename CAL>
0144   struct PFRecHitProducerKernelTopology {
0145     template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0146     ALPAKA_FN_ACC void operator()(const TAcc& acc,
0147                                   const typename CAL::TopologyTypeDevice::ConstView topology,
0148                                   reco::PFRecHitDeviceCollection::View pfRecHits,
0149                                   const uint32_t* __restrict__ denseId2pfRecHit,
0150                                   uint32_t* __restrict__ num_pfRecHits) const {
0151       // First thread updates size field pfRecHits SoA
0152       if (const int32_t thread = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc)[0u]; thread == 0)
0153         pfRecHits.size() = *num_pfRecHits;
0154 
0155       // Assign position information and associate neighbours
0156       for (int32_t i : cms::alpakatools::uniform_elements(acc, *num_pfRecHits)) {
0157         const uint32_t denseId = CAL::detId2denseId(pfRecHits.detId(i));
0158 
0159         pfRecHits.x(i) = topology.positionX(denseId);
0160         pfRecHits.y(i) = topology.positionY(denseId);
0161         pfRecHits.z(i) = topology.positionZ(denseId);
0162 
0163         for (uint32_t n = 0; n < 8; n++) {
0164           pfRecHits.neighbours(i)(n) = -1;
0165           const uint32_t denseId_neighbour = topology.neighbours(denseId)(n);
0166           if (denseId_neighbour != CAL::kInvalidDenseId) {
0167             const uint32_t pfRecHit_neighbour = denseId2pfRecHit[denseId_neighbour];
0168             if (pfRecHit_neighbour != 0xffffffff)
0169               pfRecHits.neighbours(i)(n) = (int32_t)pfRecHit_neighbour;
0170           }
0171         }
0172       }
0173     }
0174   };
0175 
0176   template <typename CAL>
0177   PFRecHitProducerKernel<CAL>::PFRecHitProducerKernel(Queue& queue, const uint32_t num_recHits)
0178       : denseId2pfRecHit_(cms::alpakatools::make_device_buffer<uint32_t[]>(queue, CAL::kSize)),
0179         num_pfRecHits_(cms::alpakatools::make_device_buffer<uint32_t>(queue)),
0180         work_div_(cms::alpakatools::make_workdiv<Acc1D>(1, 1)) {
0181     alpaka::memset(queue, denseId2pfRecHit_, 0xff);  // Reset denseId -> pfRecHit index map
0182     alpaka::memset(queue, num_pfRecHits_, 0x00);     // Reset global pfRecHit counter
0183 
0184     const uint32_t items = 64;
0185     const uint32_t groups = cms::alpakatools::divide_up_by(num_recHits, items);
0186     work_div_ = cms::alpakatools::make_workdiv<Acc1D>(groups, items);
0187   }
0188 
0189   template <typename CAL>
0190   void PFRecHitProducerKernel<CAL>::processRecHits(Queue& queue,
0191                                                    const typename CAL::CaloRecHitSoATypeDevice& recHits,
0192                                                    const typename CAL::ParameterType& params,
0193                                                    const typename CAL::TopologyTypeDevice& topology,
0194                                                    reco::PFRecHitDeviceCollection& pfRecHits) {
0195     alpaka::exec<Acc1D>(queue,
0196                         work_div_,
0197                         PFRecHitProducerKernelConstruct<CAL>{},
0198                         params.view(),
0199                         topology.view(),
0200                         recHits.view(),
0201                         pfRecHits.view(),
0202                         denseId2pfRecHit_.data(),
0203                         num_pfRecHits_.data());
0204   }
0205 
0206   template <typename CAL>
0207   void PFRecHitProducerKernel<CAL>::associateTopologyInfo(Queue& queue,
0208                                                           const typename CAL::TopologyTypeDevice& topology,
0209                                                           reco::PFRecHitDeviceCollection& pfRecHits) {
0210     alpaka::exec<Acc1D>(queue,
0211                         work_div_,
0212                         PFRecHitProducerKernelTopology<CAL>{},
0213                         topology.view(),
0214                         pfRecHits.view(),
0215                         denseId2pfRecHit_.data(),
0216                         num_pfRecHits_.data());
0217   }
0218 
0219   // Instantiate templates
0220   template class PFRecHitProducerKernel<HCAL>;
0221   template class PFRecHitProducerKernel<ECAL>;
0222 }  // 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!