PixelSeeding/plugins/gpuFishbone.h

0001 #ifndef RecoTracker_PixelSeeding_plugins_gpuFishbone_h
0002 #define RecoTracker_PixelSeeding_plugins_gpuFishbone_h
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <cstdint>
0007 #include <cstdio>
0008 #include <limits>
0009
0010 #include "DataFormats/Math/interface/approx_atan2.h"
0011 #include "HeterogeneousCore/CUDAUtilities/interface/VecArray.h"
0012 #include "HeterogeneousCore/CUDAUtilities/interface/cuda_assert.h"
0013
0014 #include "GPUCACell.h"
0015 #include "CAStructures.h"
0016
0017 namespace gpuPixelDoublets {
0018
0019   template <typename TrackerTraits>
0020   using CellNeighbors = caStructures::CellNeighborsT<TrackerTraits>;
0021   template <typename TrackerTraits>
0022   using CellTracks = caStructures::CellTracksT<TrackerTraits>;
0023   template <typename TrackerTraits>
0024   using CellNeighborsVector = caStructures::CellNeighborsVectorT<TrackerTraits>;
0025   template <typename TrackerTraits>
0026   using CellTracksVector = caStructures::CellTracksVectorT<TrackerTraits>;
0027   template <typename TrackerTraits>
0028   using OuterHitOfCell = caStructures::OuterHitOfCellT<TrackerTraits>;
0029   template <typename TrackerTraits>
0030   using HitsConstView = typename GPUCACellT<TrackerTraits>::HitsConstView;
0031
0032   template <typename TrackerTraits>
0033   __global__ void fishbone(HitsConstView<TrackerTraits> hh,
0034                            GPUCACellT<TrackerTraits>* cells,
0035                            uint32_t const* __restrict__ nCells,
0036                            OuterHitOfCell<TrackerTraits> const isOuterHitOfCellWrap,
0037                            int32_t nHits,
0038                            bool checkTrack) {
0039     constexpr auto maxCellsPerHit = GPUCACellT<TrackerTraits>::maxCellsPerHit;
0040
0041     auto const isOuterHitOfCell = isOuterHitOfCellWrap.container;
0042     int32_t offset = isOuterHitOfCellWrap.offset;
0043
0044     // x runs faster...
0045     auto firstY = threadIdx.y + blockIdx.y * blockDim.y;
0046     auto firstX = threadIdx.x;
0047
0048     float x[maxCellsPerHit], y[maxCellsPerHit], z[maxCellsPerHit], n[maxCellsPerHit];
0049     uint32_t cc[maxCellsPerHit];
0050     uint16_t d[maxCellsPerHit];
0051     uint8_t l[maxCellsPerHit];
0052
0053     for (int idy = firstY, nt = nHits - offset; idy < nt; idy += gridDim.y * blockDim.y) {
0054       auto const& vc = isOuterHitOfCell[idy];
0055       auto size = vc.size();
0056       if (size < 2)
0057         continue;
0058       // if alligned kill one of the two.
0059       // in principle one could try to relax the cut (only in r-z?) for jumping-doublets
0060       auto const& c0 = cells[vc[0]];
0061       auto xo = c0.outer_x(hh);
0062       auto yo = c0.outer_y(hh);
0063       auto zo = c0.outer_z(hh);
0064       auto sg = 0;
0065       for (int32_t ic = 0; ic < size; ++ic) {
0066         auto& ci = cells[vc[ic]];
0067         if (ci.unused())
0068           continue;  // for triplets equivalent to next
0069         if (checkTrack && ci.tracks().empty())
0070           continue;
0071         cc[sg] = vc[ic];
0072         l[sg] = ci.layerPairId();
0073         d[sg] = ci.inner_detIndex(hh);
0074         x[sg] = ci.inner_x(hh) - xo;
0075         y[sg] = ci.inner_y(hh) - yo;
0076         z[sg] = ci.inner_z(hh) - zo;
0077         n[sg] = x[sg] * x[sg] + y[sg] * y[sg] + z[sg] * z[sg];
0078         ++sg;
0079       }
0080       if (sg < 2)
0081         continue;
0082       // here we parallelize
0083       for (int32_t ic = firstX; ic < sg - 1; ic += blockDim.x) {
0084         auto& ci = cells[cc[ic]];
0085         for (auto jc = ic + 1; jc < sg; ++jc) {
0086           auto& cj = cells[cc[jc]];
0087           // must be different detectors
0088           //        if (d[ic]==d[jc]) continue;
0089           auto cos12 = x[ic] * x[jc] + y[ic] * y[jc] + z[ic] * z[jc];
0090           if (d[ic] != d[jc] && cos12 * cos12 >= 0.99999f * (n[ic] * n[jc])) {
0091             // alligned:  kill farthest (prefer consecutive layers)
0092             // if same layer prefer farthest (longer level arm) and make space for intermediate hit
0093             bool sameLayer = l[ic] == l[jc];
0094             if (n[ic] > n[jc]) {
0095               if (sameLayer) {
0096                 cj.kill();  // closest
0097                 ci.setFishbone(cj.inner_hit_id(), cj.inner_z(hh), hh);
0098               } else {
0099                 ci.kill();  // farthest
0100                 // break;  // removed to improve reproducibility. keep it for reference and tests
0101               }
0102             } else {
0103               if (!sameLayer) {
0104                 cj.kill();  // farthest
0105               } else {
0106                 ci.kill();  // closest
0107                 cj.setFishbone(ci.inner_hit_id(), ci.inner_z(hh), hh);
0108                 // break;  // removed to improve reproducibility. keep it for reference    and tests
0109               }
0110             }
0111           }
0112         }  // cj
0113       }    // ci
0114     }      // hits
0115   }
0116
0117 }  // namespace gpuPixelDoublets
0118
0119 #endif  // RecoTracker_PixelSeeding_plugins_gpuFishbone_h