src/alpaka/Hit.h

0001 #ifndef RecoTracker_LSTCore_src_alpaka_Hit_h
0002 #define RecoTracker_LSTCore_src_alpaka_Hit_h
0003
0004 #include "RecoTracker/LSTCore/interface/alpaka/Common.h"
0005 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0006 #include "RecoTracker/LSTCore/interface/alpaka/HitsDeviceCollection.h"
0007
0008 namespace ALPAKA_ACCELERATOR_NAMESPACE::lst {
0009
0010   template <typename TAcc>
0011   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float eta(TAcc const& acc, float x, float y, float z) {
0012     float r3 = alpaka::math::sqrt(acc, x * x + y * y + z * z);
0013     float rt = alpaka::math::sqrt(acc, x * x + y * y);
0014     float eta = ((z > 0) - (z < 0)) * alpaka::math::acosh(acc, r3 / rt);
0015     return eta;
0016   }
0017
0018   template <typename TAcc>
0019   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float phi_mpi_pi(TAcc const& acc, float x) {
0020     if (alpaka::math::abs(acc, x) <= kPi)
0021       return x;
0022
0023     constexpr float o2pi = 1.f / (2.f * kPi);
0024     float n = alpaka::math::round(acc, x * o2pi);
0025     return x - n * float(2.f * kPi);
0026   }
0027
0028   template <typename TAcc>
0029   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float phi(TAcc const& acc, float x, float y) {
0030     return phi_mpi_pi(acc, kPi + alpaka::math::atan2(acc, -y, -x));
0031   }
0032
0033   template <typename TAcc>
0034   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float deltaPhi(TAcc const& acc, float x1, float y1, float x2, float y2) {
0035     float phi1 = phi(acc, x1, y1);
0036     float phi2 = phi(acc, x2, y2);
0037     return phi_mpi_pi(acc, (phi2 - phi1));
0038   }
0039
0040   template <typename TAcc>
0041   ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE float deltaPhiChange(TAcc const& acc, float x1, float y1, float x2, float y2) {
0042     return deltaPhi(acc, x1, y1, x2 - x1, y2 - y1);
0043   }
0044
0045   ALPAKA_FN_ACC ALPAKA_FN_INLINE float calculate_dPhi(float phi1, float phi2) {
0046     // Calculate dPhi
0047     float dPhi = phi1 - phi2;
0048
0049     // Normalize dPhi to be between -pi and pi
0050     if (dPhi > kPi) {
0051       dPhi -= 2 * kPi;
0052     } else if (dPhi < -kPi) {
0053       dPhi += 2 * kPi;
0054     }
0055
0056     return dPhi;
0057   }
0058
0059   struct ModuleRangesKernel {
0060     template <typename TAcc>
0061     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0062                                   ModulesConst modules,
0063                                   HitsRanges hitsRanges,
0064                                   int nLowerModules) const {
0065       auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
0066       auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
0067
0068       for (int lowerIndex = globalThreadIdx[2]; lowerIndex < nLowerModules; lowerIndex += gridThreadExtent[2]) {
0069         uint16_t upperIndex = modules.partnerModuleIndices()[lowerIndex];
0070         if (hitsRanges.hitRanges()[lowerIndex][0] != -1 && hitsRanges.hitRanges()[upperIndex][0] != -1) {
0071           hitsRanges.hitRangesLower()[lowerIndex] = hitsRanges.hitRanges()[lowerIndex][0];
0072           hitsRanges.hitRangesUpper()[lowerIndex] = hitsRanges.hitRanges()[upperIndex][0];
0073           hitsRanges.hitRangesnLower()[lowerIndex] =
0074               hitsRanges.hitRanges()[lowerIndex][1] - hitsRanges.hitRanges()[lowerIndex][0] + 1;
0075           hitsRanges.hitRangesnUpper()[lowerIndex] =
0076               hitsRanges.hitRanges()[upperIndex][1] - hitsRanges.hitRanges()[upperIndex][0] + 1;
0077         }
0078       }
0079     }
0080   };
0081
0082   struct HitLoopKernel {
0083     template <typename TAcc>
0084     ALPAKA_FN_ACC void operator()(TAcc const& acc,
0085                                   uint16_t Endcap,          // Integer corresponding to endcap in module subdets
0086                                   uint16_t TwoS,            // Integer corresponding to TwoS in moduleType
0087                                   unsigned int nModules,    // Number of modules
0088                                   unsigned int nEndCapMap,  // Number of elements in endcap map
0089                                   EndcapGeometryDevConst endcapGeometry,
0090                                   ModulesConst modules,
0091                                   Hits hits,
0092                                   HitsRanges hitsRanges,
0093                                   unsigned int nHits) const  // Total number of hits in event
0094     {
0095       auto geoMapDetId = endcapGeometry.geoMapDetId();  // DetId's from endcap map
0096       auto geoMapPhi = endcapGeometry.geoMapPhi();      // Phi values from endcap map
0097       auto const globalThreadIdx = alpaka::getIdx<alpaka::Grid, alpaka::Threads>(acc);
0098       auto const gridThreadExtent = alpaka::getWorkDiv<alpaka::Grid, alpaka::Threads>(acc);
0099       for (unsigned int ihit = globalThreadIdx[2]; ihit < nHits; ihit += gridThreadExtent[2]) {
0100         float ihit_x = hits.xs()[ihit];
0101         float ihit_y = hits.ys()[ihit];
0102         float ihit_z = hits.zs()[ihit];
0103         int iDetId = hits.detid()[ihit];
0104
0105         hits.rts()[ihit] = alpaka::math::sqrt(acc, ihit_x * ihit_x + ihit_y * ihit_y);
0106         hits.phis()[ihit] = phi(acc, ihit_x, ihit_y);
0107         hits.etas()[ihit] =
0108             ((ihit_z > 0) - (ihit_z < 0)) *
0109             alpaka::math::acosh(
0110                 acc, alpaka::math::sqrt(acc, ihit_x * ihit_x + ihit_y * ihit_y + ihit_z * ihit_z) / hits.rts()[ihit]);
0111         auto found_pointer = std::lower_bound(modules.mapdetId(), modules.mapdetId() + nModules, iDetId);
0112         int found_index = std::distance(modules.mapdetId(), found_pointer);
0113         if (found_pointer == modules.mapdetId() + nModules)
0114           found_index = -1;
0115         uint16_t lastModuleIndex = modules.mapIdx()[found_index];
0116
0117         hits.moduleIndices()[ihit] = lastModuleIndex;
0118
0119         if (modules.subdets()[lastModuleIndex] == Endcap && modules.moduleType()[lastModuleIndex] == TwoS) {
0120           found_pointer = std::lower_bound(geoMapDetId, geoMapDetId + nEndCapMap, iDetId);
0121           found_index = std::distance(geoMapDetId, found_pointer);
0122           if (found_pointer == geoMapDetId + nEndCapMap)
0123             found_index = -1;
0124           float phi = geoMapPhi[found_index];
0125           float cos_phi = alpaka::math::cos(acc, phi);
0126           hits.highEdgeXs()[ihit] = ihit_x + 2.5f * cos_phi;
0127           hits.lowEdgeXs()[ihit] = ihit_x - 2.5f * cos_phi;
0128           float sin_phi = alpaka::math::sin(acc, phi);
0129           hits.highEdgeYs()[ihit] = ihit_y + 2.5f * sin_phi;
0130           hits.lowEdgeYs()[ihit] = ihit_y - 2.5f * sin_phi;
0131         }
0132         // Need to set initial value if index hasn't been seen before.
0133         int old = alpaka::atomicCas(acc,
0134                                     &(hitsRanges.hitRanges()[lastModuleIndex][0]),
0135                                     -1,
0136                                     static_cast<int>(ihit),
0137                                     alpaka::hierarchy::Threads{});
0138         // For subsequent visits, stores the min value.
0139         if (old != -1)
0140           alpaka::atomicMin(
0141               acc, &hitsRanges.hitRanges()[lastModuleIndex][0], static_cast<int>(ihit), alpaka::hierarchy::Threads{});
0142
0143         alpaka::atomicMax(
0144             acc, &hitsRanges.hitRanges()[lastModuleIndex][1], static_cast<int>(ihit), alpaka::hierarchy::Threads{});
0145       }
0146     }
0147   };
0148 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::lst
0149 #endif