src/alpaka/Hit.h

0001 #ifndef RecoTracker_LSTCore_src_alpaka_Hit_h
0002 #define RecoTracker_LSTCore_src_alpaka_Hit_h
0003
0004 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0005 #include "HeterogeneousCore/AlpakaInterface/interface/alpakastdAlgorithm.h"
0006 #include "HeterogeneousCore/AlpakaMath/interface/deltaPhi.h"
0007
0008 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0009 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0010 #include "RecoTracker/LSTCore/interface/alpaka/HitsDeviceCollection.h"
0011
0012 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0013
0014   template <typename TAcc>
0015   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float deltaPhiChange(TAcc const& acc, float x1, float y1, float x2, float y2) {
0016     return cms::alpakatools::deltaPhi(acc, x1, y1, x2 - x1, y2 - y1);
0017   }
0018
0019   struct ModuleRangesKernel {
0020     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0021                                   ModulesConst modules,
0022                                   HitsRanges hitsRanges,
0023                                   int nLowerModules) const {
0024       for (int lowerIndex : cms::alpakatools::uniform_elements(acc, nLowerModules)) {
0025         uint16_t upperIndex = modules.partnerModuleIndices()[lowerIndex];
0026         if (hitsRanges.hitRanges()[lowerIndex][0] != -1 && hitsRanges.hitRanges()[upperIndex][0] != -1) {
0027           hitsRanges.hitRangesLower()[lowerIndex] = hitsRanges.hitRanges()[lowerIndex][0];
0028           hitsRanges.hitRangesUpper()[lowerIndex] = hitsRanges.hitRanges()[upperIndex][0];
0029           hitsRanges.hitRangesnLower()[lowerIndex] =
0030               hitsRanges.hitRanges()[lowerIndex][1] - hitsRanges.hitRanges()[lowerIndex][0] + 1;
0031           hitsRanges.hitRangesnUpper()[lowerIndex] =
0032               hitsRanges.hitRanges()[upperIndex][1] - hitsRanges.hitRanges()[upperIndex][0] + 1;
0033         }
0034       }
0035     }
0036   };
0037
0038   struct HitLoopKernel {
0039     ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0040                                   uint16_t Endcap,          // Integer corresponding to endcap in module subdets
0041                                   uint16_t TwoS,            // Integer corresponding to TwoS in moduleType
0042                                   unsigned int nModules,    // Number of modules
0043                                   unsigned int nEndCapMap,  // Number of elements in endcap map
0044                                   EndcapGeometryDevConst endcapGeometry,
0045                                   ModulesConst modules,
0046                                   HitsBaseConst hitsBase,
0047                                   HitsExtended hitsExtended,
0048                                   HitsRanges hitsRanges) const  // Total number of hits in event
0049     {
0050       auto geoMapDetId = endcapGeometry.geoMapDetId();  // DetId's from endcap map
0051       auto geoMapPhi = endcapGeometry.geoMapPhi();      // Phi values from endcap map
0052       int nHits = hitsExtended.metadata().size();
0053       ALPAKA_ASSERT_ACC(nHits == hitsBase.metadata().size());
0054       for (unsigned int ihit : cms::alpakatools::uniform_elements(acc, nHits)) {
0055         float ihit_x = hitsBase.xs()[ihit];
0056         float ihit_y = hitsBase.ys()[ihit];
0057         float ihit_z = hitsBase.zs()[ihit];
0058         int iDetId = hitsBase.detid()[ihit];
0059
0060         hitsExtended.rts()[ihit] = alpaka::math::sqrt(acc, ihit_x * ihit_x + ihit_y * ihit_y);
0061         hitsExtended.phis()[ihit] = cms::alpakatools::phi(acc, ihit_x, ihit_y);
0062         hitsExtended.etas()[ihit] =
0063             ((ihit_z > 0) - (ihit_z < 0)) *
0064             alpaka::math::acosh(acc,
0065                                 alpaka::math::sqrt(acc, ihit_x * ihit_x + ihit_y * ihit_y + ihit_z * ihit_z) /
0066                                     hitsExtended.rts()[ihit]);
0067         auto found_pointer = alpaka_std::lower_bound(modules.mapdetId(), modules.mapdetId() + nModules, iDetId);
0068         ALPAKA_ASSERT_ACC(found_pointer != modules.mapdetId() + nModules);
0069         int found_index = std::distance(modules.mapdetId(), found_pointer);
0070         uint16_t lastModuleIndex = modules.mapIdx()[found_index];
0071
0072         hitsExtended.moduleIndices()[ihit] = lastModuleIndex;
0073
0074         if (modules.subdets()[lastModuleIndex] == Endcap && modules.moduleType()[lastModuleIndex] == TwoS) {
0075           found_pointer = alpaka_std::lower_bound(geoMapDetId, geoMapDetId + nEndCapMap, iDetId);
0076           ALPAKA_ASSERT_ACC(found_pointer != geoMapDetId + nEndCapMap);
0077           found_index = std::distance(geoMapDetId, found_pointer);
0078           float phi = geoMapPhi[found_index];
0079           float cos_phi = alpaka::math::cos(acc, phi);
0080           hitsExtended.highEdgeXs()[ihit] = ihit_x + 2.5f * cos_phi;
0081           hitsExtended.lowEdgeXs()[ihit] = ihit_x - 2.5f * cos_phi;
0082           float sin_phi = alpaka::math::sin(acc, phi);
0083           hitsExtended.highEdgeYs()[ihit] = ihit_y + 2.5f * sin_phi;
0084           hitsExtended.lowEdgeYs()[ihit] = ihit_y - 2.5f * sin_phi;
0085         }
0086         // Need to set initial value if index hasn't been seen before.
0087         int old = alpaka::atomicCas(acc,
0088                                     &(hitsRanges.hitRanges()[lastModuleIndex][0]),
0089                                     -1,
0090                                     static_cast<int>(ihit),
0091                                     alpaka::hierarchy::Threads{});
0092         // For subsequent visits, stores the min value.
0093         if (old != -1)
0094           alpaka::atomicMin(
0095               acc, &hitsRanges.hitRanges()[lastModuleIndex][0], static_cast<int>(ihit), alpaka::hierarchy::Threads{});
0096
0097         alpaka::atomicMax(
0098             acc, &hitsRanges.hitRanges()[lastModuleIndex][1], static_cast<int>(ihit), alpaka::hierarchy::Threads{});
0099       }
0100     }
0101   };
0102 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
0103 #endif