plugins/alpaka/CAPixelDoubletsAlgos.h

0001 #ifndef RecoTracker_PixelSeeding_plugins_alpaka_CAPixelDoubletsAlgos_h
0002 #define RecoTracker_PixelSeeding_plugins_alpaka_CAPixelDoubletsAlgos_h
0003
0004 #include <algorithm>
0005 #include <cmath>
0006 #include <cstdint>
0007 #include <cstdio>
0008 #include <limits>
0009
0010 #include <alpaka/alpaka.hpp>
0011
0012 #include "DataFormats/Math/interface/approx_atan2.h"
0013 #include "DataFormats/SiPixelClusterSoA/interface/ClusteringConstants.h"
0014 #include "DataFormats/TrackingRecHitSoA/interface/TrackingRecHitsSoA.h"
0015 #include "Geometry/CommonTopologies/interface/SimplePixelTopology.h"
0016 #include "HeterogeneousCore/AlpakaInterface/interface/VecArray.h"
0017 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0018 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0019 #include "RecoTracker/PixelSeeding/interface/CAGeometrySoA.h"
0020
0021 #include "CACell.h"
0022 #include "CAStructures.h"
0023 #include "CAHitNtupletGeneratorKernels.h"
0024
0025 // #define GPU_DEBUG
0026 // #define DOUBLETS_DEBUG
0027 // #define CA_WARNINGS
0028
0029 namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets {
0030   using namespace cms::alpakatools;
0031   using namespace ::caStructures;
0032   using namespace ::reco;
0033
0034   using HitToCell = GenericContainer;
0035
0036   template <typename TrackerTraits>
0037   using PhiBinner = PhiBinnerT<TrackerTraits>;
0038   //Move this ^ definition in CAStructures maybe
0039
0040   template <typename TAcc>
0041   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool zSizeCut(
0042       const TAcc& acc, HitsConstView hh, ::reco::CALayersSoAConstView ll, AlgoParams const& params, int i, int o) {
0043     const uint32_t mi = hh[i].detectorIndex();
0044     const auto first_forward = ll.layerStarts()[4];
0045     const auto first_bpix2 = ll.layerStarts()[1];
0046     bool innerB1 = mi < first_bpix2;
0047     bool isOuterLadder = 0 == (mi / 8) % 2;
0048     auto mes = (!innerB1) || isOuterLadder ? hh[i].clusterSizeY() : -1;
0049 #ifdef DOUBLETS_DEBUG
0050     printf("i = %d o = %d mi = %d innerB1 = %d isOuterLadder = %d first_forward = %d first_bpix2 = %d\n",
0051            i,
0052            o,
0053            mi,
0054            innerB1,
0055            isOuterLadder,
0056            first_forward,
0057            first_bpix2);
0058 #endif
0059     if (mes < 0)
0060       return false;
0061
0062     const uint32_t mo = hh[o].detectorIndex();
0063     auto so = hh[o].clusterSizeY();
0064
0065     auto dz = hh[i].zGlobal() - hh[o].zGlobal();
0066     auto dr = hh[i].rGlobal() - hh[o].rGlobal();
0067
0068     auto innerBarrel = mi < first_forward;
0069     auto onlyBarrel = mo < first_forward;
0070 #ifdef DOUBLETS_DEBUG
0071     printf("i = %d o = %d mo = %d innerB1 = %d isOuterLadder = %d \n", i, o, mo, innerBarrel, onlyBarrel);
0072 #endif
0073     if (not innerBarrel and not onlyBarrel)
0074       return false;
0075     auto dy = innerB1 ? params.maxDYsize12_ : params.maxDYsize_;
0076 #ifdef DOUBLETS_DEBUG
0077     printf("i = %d o = %d dy = %d maxDYsize12_ = %d maxDYsize_ = %d dzdrFact_ = %.2f maxDYPred_ = %d \n",
0078            i,
0079            o,
0080            dy,
0081            params.maxDYsize12_,
0082            params.maxDYsize_,
0083            params.dzdrFact_,
0084            params.maxDYPred_);
0085 #endif
0086     return onlyBarrel
0087                ? so > 0 && std::abs(so - mes) > dy
0088                : innerBarrel && std::abs(mes - int(std::abs(dz / dr) * params.dzdrFact_ + 0.5f)) > params.maxDYPred_;
0089   }
0090
0091   template <typename TAcc>
0092   ALPAKA_FN_ACC ALPAKA_FN_INLINE bool clusterCut(
0093       const TAcc& acc, HitsConstView hh, ::reco::CALayersSoAConstView ll, AlgoParams const& params, uint32_t i) {
0094     const uint32_t mi = hh[i].detectorIndex();
0095     const auto first_bpix2 = ll.layerStarts()[1];
0096     const auto first_bpix3 = ll.layerStarts()[2];
0097     bool innerB1orB2 = mi < ll.layerStarts()[2];
0098 #ifdef DOUBLETS_DEBUG
0099     printf(
0100         "i = %d mi = %d innerB1orB2 = %d innerB1 = %d innerB2 = %d minYsizeB1 = %d minYsizeB2 = %d isOuterLadder = %d "
0101         "mes = %d \n",
0102         i,
0103         mi,
0104         innerB1orB2,
0105         mi < first_bpix2,
0106         (mi >= first_bpix2) && (mi < first_bpix3),
0107         params.minYsizeB1_,
0108         params.minYsizeB2_,
0109         (0 == (mi / 8) % 2),
0110         (!(mi < first_bpix2)) || (0 == (mi / 8) % 2) ? hh[i].clusterSizeY() : -1);
0111 #endif
0112     if (!innerB1orB2)
0113       return false;
0114
0115     bool innerB1 = mi < first_bpix2;
0116
0117     bool isOuterLadder = 0 == (mi / 8) % 2;
0118     auto mes = (!innerB1) || isOuterLadder ? hh[i].clusterSizeY() : -1;
0119
0120     if (innerB1)  // B1
0121       if (mes > 0 && mes < params.minYsizeB1_)
0122         return true;  // only long cluster  (5*8)
0123     bool innerB2 = (mi >= first_bpix2) && (mi < first_bpix3);
0124     if (innerB2)  // B2 and F1
0125       if (mes > 0 && mes < params.minYsizeB2_)
0126         return true;
0127
0128     return false;
0129   }
0130
0131   template <typename TrackerTraits, typename TAcc>
0132   ALPAKA_FN_ACC ALPAKA_FN_INLINE void doubletsFromHisto(const TAcc& acc,
0133                                                         uint32_t maxNumOfDoublets,
0134                                                         CACell<TrackerTraits>* cells,
0135                                                         uint32_t* nCells,
0136                                                         HitsConstView hh,
0137                                                         ::reco::CAGraphSoAConstView cc,
0138                                                         ::reco::CALayersSoAConstView ll,
0139                                                         uint32_t const* __restrict__ offsets,
0140                                                         PhiBinner<TrackerTraits> const* phiBinner,
0141                                                         HitToCell* outerHitHisto,
0142                                                         AlgoParams const& params) {
0143     const bool doClusterCut = params.minYsizeB1_ > 0 or params.minYsizeB2_ > 0;
0144     const bool doZSizeCut = params.maxDYsize12_ > 0 or params.maxDYsize_ > 0 or params.maxDYPred_ > 0;
0145
0146     // cm (1 GeV track has 1 GeV/c / (e * 3.8T) ~ 87 cm radius in a 3.8T field)
0147     const float minRadius = params.cellPtCut_ * 87.78f;
0148     const float minRadius2T4 = 4.f * minRadius * minRadius;
0149
0150     const uint32_t nPairs = cc.metadata().size();
0151     using PhiHisto = PhiBinner<TrackerTraits>;
0152     ALPAKA_ASSERT_ACC(offsets);
0153
0154     auto layerSize = [=](uint8_t li) { return offsets[li + 1] - offsets[li]; };
0155
0156     // nPairsMax to be optimized later (originally was 64).
0157     // If it should much be bigger, consider using a block-wide parallel prefix scan,
0158     // e.g. see  https://nvlabs.github.io/cub/classcub_1_1_warp_scan.html
0159     auto& innerLayerCumulativeSize = alpaka::declareSharedVar<uint32_t[pixelTopology::maxPairs], __COUNTER__>(acc);
0160     auto& ntot = alpaka::declareSharedVar<uint32_t, __COUNTER__>(acc);
0161
0162 #ifdef DOUBLETS_DEBUG
0163     if (cms::alpakatools::once_per_grid(acc))
0164       printf(
0165           "maxNumDoublets = %d  cc.metadata().size() = %d ll.metadata().size() = %d cellZ0Cut_ = %.2f cellPtCut_ = "
0166           "%.2f doClusterCut = %d doZ0Cut = %d  doPtCut = %d doZSizeCut = %d\n",
0167           maxNumOfDoublets,
0168           cc.metadata().size(),
0169           ll.metadata().size(),
0170           params.cellZ0Cut_,
0171           params.cellPtCut_,
0172           doClusterCut,
0173           params.cellZ0Cut_ > 0,
0174           params.cellPtCut_ > 0,
0175           doZSizeCut);
0176 #endif
0177
0178     if (cms::alpakatools::once_per_block(acc)) {
0179       innerLayerCumulativeSize[0] = layerSize(cc.graph()[0][0]);
0180       for (uint32_t i = 1; i < nPairs; ++i) {
0181         innerLayerCumulativeSize[i] = innerLayerCumulativeSize[i - 1] + layerSize(cc.graph()[i][0]);
0182       }
0183       ntot = innerLayerCumulativeSize[nPairs - 1];
0184     }
0185     alpaka::syncBlockThreads(acc);
0186
0187     // declared outside the loop, as it cannot go backward
0188     uint32_t pairLayerId = 0;
0189
0190     // outermost parallel loop, using all grid elements along the slower dimension (Y or 0 in a 2D grid)
0191     for (uint32_t j : cms::alpakatools::uniform_elements_y(acc, ntot)) {
0192       // move to lower_bound ?
0193       while (j >= innerLayerCumulativeSize[pairLayerId++])
0194         ;
0195       --pairLayerId;
0196
0197       ALPAKA_ASSERT_ACC(pairLayerId < nPairs);
0198       ALPAKA_ASSERT_ACC(j < innerLayerCumulativeSize[pairLayerId]);
0199       ALPAKA_ASSERT_ACC(0 == pairLayerId || j >= innerLayerCumulativeSize[pairLayerId - 1]);
0200
0201       uint8_t inner = cc.graph()[pairLayerId][0];
0202       uint8_t outer = cc.graph()[pairLayerId][1];
0203       ALPAKA_ASSERT_ACC(outer > inner);
0204
0205       auto hoff = PhiHisto::histOff(outer);
0206       auto i = (0 == pairLayerId) ? j : j - innerLayerCumulativeSize[pairLayerId - 1];
0207       i += offsets[inner];
0208
0209       ALPAKA_ASSERT_ACC(i >= offsets[inner]);
0210       ALPAKA_ASSERT_ACC(i < offsets[inner + 1]);
0211 #ifdef DOUBLETS_DEBUG
0212       printf("pairLayerId = %d i = %d inner = %d outer = %d offsets[inner] = %d offsets[inner + 1] = %d\n",
0213              pairLayerId,
0214              i,
0215              inner,
0216              outer,
0217              offsets[inner],
0218              offsets[inner + 1]);
0219 #endif
0220       // found hit corresponding to our worker thread, now do the job
0221       if (hh[i].detectorIndex() > ll.layerStarts()[ll.metadata().size() - 1])  //TODO use cc
0222         continue;                                                              // invalid
0223
0224       /* maybe clever, not effective when zoCut is on
0225       auto bpos = (mi%8)/4;  // if barrel is 1 for z>0
0226       auto fpos = (outer>3) & (outer<7);
0227       if ( ((inner<3) & (outer>3)) && bpos!=fpos) continue;
0228       */
0229
0230       auto mez = hh[i].zGlobal();
0231
0232       if (mez < cc.minz()[pairLayerId] || mez > cc.maxz()[pairLayerId])
0233         continue;
0234 #ifdef DOUBLETS_DEBUG
0235       if (doClusterCut && outer > pixelTopology::last_barrel_layer)
0236         printf("clustCut: %d %d \n", i, clusterCut<TAcc>(acc, hh, ll, params, i));
0237 #endif
0238
0239       if (doClusterCut && outer > pixelTopology::last_barrel_layer && clusterCut<TAcc>(acc, hh, ll, params, i))
0240         continue;
0241
0242       auto mep = hh[i].iphi();
0243       auto mer = hh[i].rGlobal();
0244
0245       // all cuts: true if fails
0246       auto ptcut = [&](int j, int16_t idphi) {
0247         auto r2t4 = minRadius2T4;
0248         auto ri = mer;
0249         auto ro = hh[j].rGlobal();
0250         auto dphi = short2phi(idphi);
0251         return dphi * dphi * (r2t4 - ri * ro) > (ro - ri) * (ro - ri);
0252       };
0253       auto z0cutoff = [&](int j) {
0254         auto zo = hh[j].zGlobal();
0255         auto ro = hh[j].rGlobal();
0256         auto dr = ro - mer;
0257         return dr > cc.maxr()[pairLayerId] || dr < 0 || std::abs((mez * ro - mer * zo)) > params.cellZ0Cut_ * dr;
0258       };
0259
0260       auto iphicut = cc.phiCuts()[pairLayerId];
0261
0262       auto kl = PhiHisto::bin(int16_t(mep - iphicut));
0263       auto kh = PhiHisto::bin(int16_t(mep + iphicut));
0264       auto incr = [](auto& k) { return k = (k + 1) % PhiHisto::nbins(); };
0265
0266 #ifdef GPU_DEGBU
0267       printf("pairLayerId %d %d %.2f %.2f %.2f \n",
0268              pairLayerId,
0269              cc.phiCuts()[pairLayerId],
0270              cc.maxr()[pairLayerId],
0271              cc.maxz()[pairLayerId],
0272              cc.minz()[pairLayerId]);
0273 #endif
0274
0275       auto khh = kh;
0276       incr(khh);
0277       for (auto kk = kl; kk != khh; incr(kk)) {
0278         //#ifdef GPU_DEBUG
0279         //        if (kk != kl && kk != kh)
0280         //          nmin += phiBinner->size(kk + hoff);
0281         //#endif
0282
0283         auto const* __restrict__ p = phiBinner->begin(kk + hoff);
0284         auto const* __restrict__ e = phiBinner->end(kk + hoff);
0285         auto const maxpIndex = e - p;
0286
0287         // innermost parallel loop, using the block elements along the faster dimension (X or 1 in a 2D grid)
0288         for (uint32_t pIndex : cms::alpakatools::independent_group_elements_x(acc, maxpIndex)) {
0289           // FIXME implement alpaka::ldg and use it here? or is it const* __restrict__ enough?
0290           auto oi = p[pIndex];
0291           ALPAKA_ASSERT_ACC(oi >= offsets[outer]);
0292           ALPAKA_ASSERT_ACC(oi < offsets[outer + 1]);
0293           auto mo = hh[oi].detectorIndex();
0294
0295           // invalid
0296           if (mo > pixelClustering::maxNumModules)  //FIXME use cc?
0297             continue;
0298
0299           if (params.cellZ0Cut_ > 0. && z0cutoff(oi))
0300             continue;
0301
0302           auto mop = hh[oi].iphi();
0303           uint16_t idphi = std::min(std::abs(int16_t(mop - mep)), std::abs(int16_t(mep - mop)));
0304
0305           if (idphi > iphicut)
0306             continue;
0307 #ifdef DOUBLETS_DEBUG
0308           printf("zSizeCut: %d %d %d \n", i, oi, zSizeCut<TAcc>(acc, hh, ll, params, i, oi));
0309 #endif
0310           if (doZSizeCut && zSizeCut<TAcc>(acc, hh, ll, params, i, oi))
0311             continue;
0312
0313           if (params.cellPtCut_ > 0. && ptcut(oi, idphi))
0314             continue;
0315
0316           auto ind = alpaka::atomicAdd(acc, nCells, 1u, alpaka::hierarchy::Blocks{});
0317           if (ind >= maxNumOfDoublets) {
0318 #ifdef CA_WARNINGS
0319             printf("Warning!!!! Too many cells (limit = %d)!\n", maxNumOfDoublets);
0320 #endif
0321             alpaka::atomicSub(acc, nCells, 1u, alpaka::hierarchy::Blocks{});
0322             break;
0323           }
0324
0325           outerHitHisto->count(acc, oi - hh.offsetBPIX2());
0326           cells[ind].init(hh, pairLayerId, inner, outer, i, oi);
0327 #ifdef DOUBLETS_DEBUG
0328           printf("doublet: %d layerPair: %d inner: %d outer: %d i: %d oi: %d\n", ind, pairLayerId, inner, outer, i, oi);
0329 #endif
0330         }
0331       }
0332     }  // loop in block...
0333   }
0334
0335 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE::caPixelDoublets
0336
0337 #endif  // RecoTracker_PixelSeeding_plugins_alpaka_CAPixelDoubletsAlgos_h