plugins/alpaka/CACell.h

0001 #ifndef RecoTracker_PixelSeeding_plugins_alpaka_CACell_h
0002 #define RecoTracker_PixelSeeding_plugins_alpaka_CACell_h
0003
0004 // #define ONLY_TRIPLETS_IN_HOLE
0005
0006 #include <cmath>
0007 #include <limits>
0008
0009 #include <alpaka/alpaka.hpp>
0010
0011 #include "DataFormats/TrackSoA/interface/TrackDefinitions.h"
0012 #include "DataFormats/TrackSoA/interface/TracksSoA.h"
0013 #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h"
0014 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0015 #include "HeterogeneousCore/AlpakaInterface/interface/SimpleVector.h"
0016 #include "HeterogeneousCore/AlpakaInterface/interface/VecArray.h"
0017 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0018 #include "RecoTracker/PixelSeeding/interface/CircleEq.h"
0019
0020 #include "CAStructures.h"
0021
0022 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0023   template <typename TrackerTraits>
0024   class CACellT {
0025   public:
0026     using PtrAsInt = unsigned long long;
0027
0028     static constexpr auto maxCellsPerHit = TrackerTraits::maxCellsPerHit;
0029     using OuterHitOfCellContainer = caStructures::OuterHitOfCellContainerT<TrackerTraits>;
0030     using OuterHitOfCell = caStructures::OuterHitOfCellT<TrackerTraits>;
0031     using CellNeighbors = caStructures::CellNeighborsT<TrackerTraits>;
0032     using CellTracks = caStructures::CellTracksT<TrackerTraits>;
0033     using CellNeighborsVector = caStructures::CellNeighborsVectorT<TrackerTraits>;
0034     using CellTracksVector = caStructures::CellTracksVectorT<TrackerTraits>;
0035
0036     using HitsConstView = TrackingRecHitSoAConstView<TrackerTraits>;
0037     using hindex_type = typename TrackerTraits::hindex_type;
0038     using tindex_type = typename TrackerTraits::tindex_type;
0039     static constexpr auto invalidHitId = std::numeric_limits<hindex_type>::max();
0040
0041     using TmpTuple = cms::alpakatools::VecArray<uint32_t, TrackerTraits::maxDepth>;
0042
0043     using HitContainer = typename reco::TrackSoA<TrackerTraits>::HitContainer;
0044     using Quality = ::pixelTrack::Quality;
0045     static constexpr auto bad = ::pixelTrack::Quality::bad;
0046
0047     enum class StatusBit : uint16_t { kUsed = 1, kInTrack = 2, kKilled = 1 << 15 };
0048
0049     CACellT() = default;
0050
0051     ALPAKA_FN_ACC ALPAKA_FN_INLINE void init(CellNeighborsVector& cellNeighbors,
0052                                              CellTracksVector& cellTracks,
0053                                              const HitsConstView& hh,
0054                                              int layerPairId,
0055                                              hindex_type innerHitId,
0056                                              hindex_type outerHitId) {
0057       theInnerHitId = innerHitId;
0058       theOuterHitId = outerHitId;
0059       theLayerPairId_ = layerPairId;
0060       theStatus_ = 0;
0061       theFishboneId = invalidHitId;
0062
0063       // optimization that depends on access pattern
0064       theInnerZ = hh[innerHitId].zGlobal();
0065       theInnerR = hh[innerHitId].rGlobal();
0066
0067       // link to default empty
0068       theOuterNeighbors = &cellNeighbors[0];
0069       theTracks = &cellTracks[0];
0070       ALPAKA_ASSERT_ACC(outerNeighbors().empty());
0071       ALPAKA_ASSERT_ACC(tracks().empty());
0072     }
0073
0074     template <typename TAcc>
0075     ALPAKA_FN_ACC ALPAKA_FN_INLINE __attribute__((always_inline)) int addOuterNeighbor(
0076         const TAcc& acc, typename TrackerTraits::cindex_type t, CellNeighborsVector& cellNeighbors) {
0077       // use smart cache
0078       if (outerNeighbors().empty()) {
0079         auto i = cellNeighbors.extend(acc);  // maybe wasted....
0080         if (i > 0) {
0081           cellNeighbors[i].reset();
0082           alpaka::mem_fence(acc, alpaka::memory_scope::Grid{});
0083 #ifdef ALPAKA_ACC_CPU_B_SEQ_T_SEQ_ENABLED
0084           theOuterNeighbors = &cellNeighbors[i];
0085 #else
0086           auto zero = (PtrAsInt)(&cellNeighbors[0]);
0087           alpaka::atomicCas(acc,
0088                             (PtrAsInt*)(&theOuterNeighbors),
0089                             zero,
0090                             (PtrAsInt)(&cellNeighbors[i]),
0091                             alpaka::hierarchy::Blocks{});  // if fails we cannot give "i" back...
0092 #endif
0093         } else
0094           return -1;
0095       }
0096       alpaka::mem_fence(acc, alpaka::memory_scope::Grid{});
0097       return outerNeighbors().push_back(acc, t);
0098     }
0099
0100     template <typename TAcc>
0101     ALPAKA_FN_ACC ALPAKA_FN_INLINE __attribute__((always_inline)) int addTrack(TAcc const& acc,
0102                                                                                tindex_type t,
0103                                                                                CellTracksVector& cellTracks) {
0104       if (tracks().empty()) {
0105         auto i = cellTracks.extend(acc);  // maybe wasted....
0106         if (i > 0) {
0107           cellTracks[i].reset();
0108           alpaka::mem_fence(acc, alpaka::memory_scope::Grid{});
0109 #ifdef ALPAKA_ACC_CPU_B_SEQ_T_SEQ_ENABLED
0110           theTracks = &cellTracks[i];
0111 #else
0112           auto zero = (PtrAsInt)(&cellTracks[0]);
0113           alpaka::atomicCas(acc,
0114                             (PtrAsInt*)(&theTracks),
0115                             zero,
0116                             (PtrAsInt)(&cellTracks[i]),
0117                             alpaka::hierarchy::Blocks{});  // if fails we cannot give "i" back...
0118
0119 #endif
0120         } else
0121           return -1;
0122       }
0123       alpaka::mem_fence(acc, alpaka::memory_scope::Grid{});
0124       return tracks().push_back(acc, t);
0125     }
0126
0127     ALPAKA_FN_ACC ALPAKA_FN_INLINE CellTracks& tracks() { return *theTracks; }
0128     ALPAKA_FN_ACC ALPAKA_FN_INLINE CellTracks const& tracks() const { return *theTracks; }
0129     ALPAKA_FN_ACC ALPAKA_FN_INLINE CellNeighbors& outerNeighbors() { return *theOuterNeighbors; }
0130     ALPAKA_FN_ACC ALPAKA_FN_INLINE CellNeighbors const& outerNeighbors() const { return *theOuterNeighbors; }
0131     ALPAKA_FN_ACC ALPAKA_FN_INLINE float inner_x(const HitsConstView& hh) const { return hh[theInnerHitId].xGlobal(); }
0132     ALPAKA_FN_ACC ALPAKA_FN_INLINE float outer_x(const HitsConstView& hh) const { return hh[theOuterHitId].xGlobal(); }
0133     ALPAKA_FN_ACC ALPAKA_FN_INLINE float inner_y(const HitsConstView& hh) const { return hh[theInnerHitId].yGlobal(); }
0134     ALPAKA_FN_ACC ALPAKA_FN_INLINE float outer_y(const HitsConstView& hh) const { return hh[theOuterHitId].yGlobal(); }
0135     ALPAKA_FN_ACC ALPAKA_FN_INLINE float inner_z(const HitsConstView& hh) const { return theInnerZ; }
0136     ALPAKA_FN_ACC ALPAKA_FN_INLINE float outer_z(const HitsConstView& hh) const { return hh[theOuterHitId].zGlobal(); }
0137     ALPAKA_FN_ACC ALPAKA_FN_INLINE float inner_r(const HitsConstView& hh) const { return theInnerR; }
0138     ALPAKA_FN_ACC ALPAKA_FN_INLINE float outer_r(const HitsConstView& hh) const { return hh[theOuterHitId].rGlobal(); }
0139
0140     ALPAKA_FN_ACC ALPAKA_FN_INLINE auto inner_iphi(const HitsConstView& hh) const { return hh[theInnerHitId].iphi(); }
0141     ALPAKA_FN_ACC ALPAKA_FN_INLINE auto outer_iphi(const HitsConstView& hh) const { return hh[theOuterHitId].iphi(); }
0142
0143     ALPAKA_FN_ACC ALPAKA_FN_INLINE float inner_detIndex(const HitsConstView& hh) const {
0144       return hh[theInnerHitId].detectorIndex();
0145     }
0146     ALPAKA_FN_ACC ALPAKA_FN_INLINE float outer_detIndex(const HitsConstView& hh) const {
0147       return hh[theOuterHitId].detectorIndex();
0148     }
0149
0150     constexpr unsigned int inner_hit_id() const { return theInnerHitId; }
0151     constexpr unsigned int outer_hit_id() const { return theOuterHitId; }
0152
0153     ALPAKA_FN_ACC void print_cell() const {
0154       printf("printing cell: on layerPair: %d, innerHitId: %d, outerHitId: %d \n",
0155              theLayerPairId_,
0156              theInnerHitId,
0157              theOuterHitId);
0158     }
0159
0160     ALPAKA_FN_ACC bool check_alignment(const HitsConstView& hh,
0161                                        CACellT const& otherCell,
0162                                        const float ptmin,
0163                                        const float hardCurvCut,
0164                                        const float caThetaCutBarrel,
0165                                        const float caThetaCutForward,
0166                                        const float dcaCutInnerTriplet,
0167                                        const float dcaCutOuterTriplet) const {
0168       // detIndex of the layerStart for the Phase1 Pixel Detector:
0169       // [BPX1, BPX2, BPX3, BPX4,  FP1,  FP2,  FP3,  FN1,  FN2,  FN3, LAST_VALID]
0170       // [   0,   96,  320,  672, 1184, 1296, 1408, 1520, 1632, 1744,       1856]
0171       auto ri = inner_r(hh);
0172       auto zi = inner_z(hh);
0173
0174       auto ro = outer_r(hh);
0175       auto zo = outer_z(hh);
0176
0177       auto r1 = otherCell.inner_r(hh);
0178       auto z1 = otherCell.inner_z(hh);
0179       auto isBarrel = otherCell.outer_detIndex(hh) < TrackerTraits::last_barrel_detIndex;
0180       // TODO tune CA cuts below (theta and dca)
0181       bool aligned = areAlignedRZ(r1, z1, ri, zi, ro, zo, ptmin, isBarrel ? caThetaCutBarrel : caThetaCutForward);
0182       return (aligned && dcaCut(hh,
0183                                 otherCell,
0184                                 otherCell.inner_detIndex(hh) < TrackerTraits::last_bpix1_detIndex ? dcaCutInnerTriplet
0185                                                                                                   : dcaCutOuterTriplet,
0186                                 hardCurvCut));
0187     }
0188
0189     ALPAKA_FN_ACC ALPAKA_FN_INLINE __attribute__((always_inline)) static bool areAlignedRZ(
0190         float r1, float z1, float ri, float zi, float ro, float zo, const float ptmin, const float thetaCut) {
0191       float radius_diff = std::abs(r1 - ro);
0192       float distance_13_squared = radius_diff * radius_diff + (z1 - zo) * (z1 - zo);
0193
0194       float pMin = ptmin * std::sqrt(distance_13_squared);  // this needs to be divided by
0195                                                             // radius_diff later
0196
0197       float tan_12_13_half_mul_distance_13_squared = fabs(z1 * (ri - ro) + zi * (ro - r1) + zo * (r1 - ri));
0198       return tan_12_13_half_mul_distance_13_squared * pMin <= thetaCut * distance_13_squared * radius_diff;
0199     }
0200
0201     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool dcaCut(const HitsConstView& hh,
0202                                                CACellT const& otherCell,
0203                                                const float region_origin_radius_plus_tolerance,
0204                                                const float maxCurv) const {
0205       auto x1 = otherCell.inner_x(hh);
0206       auto y1 = otherCell.inner_y(hh);
0207
0208       auto x2 = inner_x(hh);
0209       auto y2 = inner_y(hh);
0210
0211       auto x3 = outer_x(hh);
0212       auto y3 = outer_y(hh);
0213
0214       CircleEq<float> eq(x1, y1, x2, y2, x3, y3);
0215
0216       if (eq.curvature() > maxCurv)
0217         return false;
0218
0219       return std::abs(eq.dca0()) < region_origin_radius_plus_tolerance * std::abs(eq.curvature());
0220     }
0221
0222     ALPAKA_FN_ACC ALPAKA_FN_INLINE __attribute__((always_inline)) static bool dcaCutH(
0223         float x1,
0224         float y1,
0225         float x2,
0226         float y2,
0227         float x3,
0228         float y3,
0229         const float region_origin_radius_plus_tolerance,
0230         const float maxCurv) {
0231       CircleEq<float> eq(x1, y1, x2, y2, x3, y3);
0232
0233       if (eq.curvature() > maxCurv)
0234         return false;
0235
0236       return std::abs(eq.dca0()) < region_origin_radius_plus_tolerance * std::abs(eq.curvature());
0237     }
0238
0239     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool hole0(const HitsConstView& hh, CACellT const& innerCell) const {
0240       using namespace phase1PixelTopology;
0241
0242       int p = innerCell.inner_iphi(hh);
0243       if (p < 0)
0244         p += std::numeric_limits<unsigned short>::max();
0245       p = (max_ladder_bpx0 * p) / std::numeric_limits<unsigned short>::max();
0246       p %= max_ladder_bpx0;
0247       auto il = first_ladder_bpx0 + p;
0248       auto r0 = hh.averageGeometry().ladderR[il];
0249       auto ri = innerCell.inner_r(hh);
0250       auto zi = innerCell.inner_z(hh);
0251       auto ro = outer_r(hh);
0252       auto zo = outer_z(hh);
0253       auto z0 = zi + (r0 - ri) * (zo - zi) / (ro - ri);
0254       auto z_in_ladder = std::abs(z0 - hh.averageGeometry().ladderZ[il]);
0255       auto z_in_module = z_in_ladder - module_length_bpx0 * int(z_in_ladder / module_length_bpx0);
0256       auto gap = z_in_module < module_tolerance_bpx0 || z_in_module > (module_length_bpx0 - module_tolerance_bpx0);
0257       return gap;
0258     }
0259
0260     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool hole4(const HitsConstView& hh, CACellT const& innerCell) const {
0261       using namespace phase1PixelTopology;
0262
0263       int p = outer_iphi(hh);
0264       if (p < 0)
0265         p += std::numeric_limits<unsigned short>::max();
0266       p = (max_ladder_bpx4 * p) / std::numeric_limits<unsigned short>::max();
0267       p %= max_ladder_bpx4;
0268       auto il = first_ladder_bpx4 + p;
0269       auto r4 = hh.averageGeometry().ladderR[il];
0270       auto ri = innerCell.inner_r(hh);
0271       auto zi = innerCell.inner_z(hh);
0272       auto ro = outer_r(hh);
0273       auto zo = outer_z(hh);
0274       auto z4 = zo + (r4 - ro) * (zo - zi) / (ro - ri);
0275       auto z_in_ladder = std::abs(z4 - hh.averageGeometry().ladderZ[il]);
0276       auto z_in_module = z_in_ladder - module_length_bpx4 * int(z_in_ladder / module_length_bpx4);
0277       auto gap = z_in_module < module_tolerance_bpx4 || z_in_module > (module_length_bpx4 - module_tolerance_bpx4);
0278       auto holeP = z4 > hh.averageGeometry().ladderMaxZ[il] && z4 < hh.averageGeometry().endCapZ[0];
0279       auto holeN = z4 < hh.averageGeometry().ladderMinZ[il] && z4 > hh.averageGeometry().endCapZ[1];
0280       return gap || holeP || holeN;
0281     }
0282
0283     // trying to free the track building process from hardcoded layers, leaving
0284     // the visit of the graph based on the neighborhood connections between cells.
0285     template <int DEPTH, typename TAcc>
0286     ALPAKA_FN_ACC ALPAKA_FN_INLINE void find_ntuplets(TAcc const& acc,
0287                                                       const HitsConstView& hh,
0288                                                       CACellT* __restrict__ cells,
0289                                                       CellTracksVector& cellTracks,
0290                                                       HitContainer& foundNtuplets,
0291                                                       cms::alpakatools::AtomicPairCounter& apc,
0292                                                       Quality* __restrict__ quality,
0293                                                       TmpTuple& tmpNtuplet,
0294                                                       const unsigned int minHitsPerNtuplet,
0295                                                       bool startAt0) const {
0296       // the building process for a track ends if:
0297       // it has no right neighbor
0298       // it has no compatible neighbor
0299       // the ntuplets is then saved if the number of hits it contains is greater
0300       // than a threshold
0301
0302       if constexpr (DEPTH <= 0) {
0303         printf("ERROR: CACellT::find_ntuplets reached full depth!\n");
0304         ALPAKA_ASSERT_ACC(false);
0305       } else {
0306         auto doubletId = this - cells;
0307         tmpNtuplet.push_back_unsafe(doubletId);
0308         ALPAKA_ASSERT_ACC(tmpNtuplet.size() <= int(TrackerTraits::maxHitsOnTrack - 3));
0309
0310         bool last = true;
0311         for (unsigned int otherCell : outerNeighbors()) {
0312           if (cells[otherCell].isKilled())
0313             continue;  // killed by earlyFishbone
0314           last = false;
0315           cells[otherCell].template find_ntuplets<DEPTH - 1>(
0316               acc, hh, cells, cellTracks, foundNtuplets, apc, quality, tmpNtuplet, minHitsPerNtuplet, startAt0);
0317         }
0318         if (last) {  // if long enough save...
0319           if ((unsigned int)(tmpNtuplet.size()) >= minHitsPerNtuplet - 1) {
0320 #ifdef ONLY_TRIPLETS_IN_HOLE
0321             // triplets accepted only pointing to the hole
0322             if (tmpNtuplet.size() >= 3 || (startAt0 && hole4(hh, cells[tmpNtuplet[0]])) ||
0323                 ((!startAt0) && hole0(hh, cells[tmpNtuplet[0]])))
0324 #endif
0325             {
0326               hindex_type hits[TrackerTraits::maxDepth + 2];
0327               auto nh = 0U;
0328               constexpr int maxFB = 2;  // for the time being let's limit this
0329               int nfb = 0;
0330               for (auto c : tmpNtuplet) {
0331                 hits[nh++] = cells[c].theInnerHitId;
0332                 if (nfb < maxFB && cells[c].hasFishbone()) {
0333                   ++nfb;
0334                   hits[nh++] = cells[c].theFishboneId;  // Fishbone hit is always outer than inner hit
0335                 }
0336               }
0337               ALPAKA_ASSERT_ACC(nh < TrackerTraits::maxHitsOnTrack);
0338               hits[nh] = theOuterHitId;
0339               auto it = foundNtuplets.bulkFill(acc, apc, hits, nh + 1);
0340               if (it >= 0) {  // if negative is overflow....
0341                 for (auto c : tmpNtuplet)
0342                   cells[c].addTrack(acc, it, cellTracks);
0343                 quality[it] = bad;  // initialize to bad
0344               }
0345             }
0346           }
0347         }
0348         tmpNtuplet.pop_back();
0349         ALPAKA_ASSERT_ACC(tmpNtuplet.size() < int(TrackerTraits::maxHitsOnTrack - 1));
0350       }
0351     }
0352
0353     // Cell status management
0354     ALPAKA_FN_ACC ALPAKA_FN_INLINE void kill() { theStatus_ |= uint16_t(StatusBit::kKilled); }
0355     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool isKilled() const { return theStatus_ & uint16_t(StatusBit::kKilled); }
0356
0357     ALPAKA_FN_ACC ALPAKA_FN_INLINE int16_t layerPairId() const { return theLayerPairId_; }
0358
0359     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool unused() const { return 0 == (uint16_t(StatusBit::kUsed) & theStatus_); }
0360     ALPAKA_FN_ACC ALPAKA_FN_INLINE void setStatusBits(StatusBit mask) { theStatus_ |= uint16_t(mask); }
0361
0362     template <typename TAcc>
0363     ALPAKA_FN_ACC ALPAKA_FN_INLINE void setFishbone(TAcc const& acc, hindex_type id, float z, const HitsConstView& hh) {
0364       // make it deterministic: use the farther apart (in z)
0365       auto old = theFishboneId;
0366       while (old !=
0367              alpaka::atomicCas(
0368                  acc,
0369                  &theFishboneId,
0370                  old,
0371                  (invalidHitId == old || std::abs(z - theInnerZ) > std::abs(hh[old].zGlobal() - theInnerZ)) ? id : old,
0372                  alpaka::hierarchy::Blocks{}))
0373         old = theFishboneId;
0374     }
0375     ALPAKA_FN_ACC ALPAKA_FN_INLINE auto fishboneId() const { return theFishboneId; }
0376     ALPAKA_FN_ACC ALPAKA_FN_INLINE bool hasFishbone() const { return theFishboneId != invalidHitId; }
0377
0378   private:
0379     CellNeighbors* theOuterNeighbors;
0380     CellTracks* theTracks;
0381
0382     int16_t theLayerPairId_;
0383     uint16_t theStatus_;  // tbd
0384
0385     float theInnerZ;
0386     float theInnerR;
0387     hindex_type theInnerHitId;
0388     hindex_type theOuterHitId;
0389     hindex_type theFishboneId;
0390   };
0391
0392 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE
0393
0394 #endif  // RecoTracker_PixelSeeding_plugins_alpaka_CACell_h