plugins/alpaka/vertexFinder.dev.cc

0001 #include <alpaka/alpaka.hpp>
0002
0003 #include "DataFormats/TrackSoA/interface/alpaka/TrackUtilities.h"
0004 #include "HeterogeneousCore/AlpakaInterface/interface/config.h"
0005 #include "HeterogeneousCore/AlpakaInterface/interface/workdivision.h"
0006 #include "RecoTracker/PixelVertexFinding/interface/PixelVertexWorkSpaceLayout.h"
0007 #include "RecoTracker/PixelVertexFinding/plugins/alpaka/PixelVertexWorkSpaceSoADeviceAlpaka.h"
0008
0009 #include "vertexFinder.h"
0010 #include "clusterTracksDBSCAN.h"
0011 #include "clusterTracksIterative.h"
0012 #include "clusterTracksByDensity.h"
0013 #include "fitVertices.h"
0014 #include "sortByPt2.h"
0015 #include "splitVertices.h"
0016
0017 #undef PIXVERTEX_DEBUG_PRODUCE
0018 namespace ALPAKA_ACCELERATOR_NAMESPACE {
0019   namespace vertexFinder {
0020     using namespace cms::alpakatools;
0021     // reject outlier tracks that contribute more than this to the chi2 of the vertex fit
0022     constexpr float maxChi2ForFirstFit = 50.f;
0023     constexpr float maxChi2ForFinalFit = 5000.f;
0024
0025     // split vertices with a chi2/NDoF greater than this
0026     constexpr float maxChi2ForSplit = 9.f;
0027
0028     template <typename TrackerTraits>
0029     class LoadTracks {
0030     public:
0031       template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0032       ALPAKA_FN_ACC void operator()(const TAcc& acc,
0033                                     reco::TrackSoAConstView<TrackerTraits> tracks_view,
0034                                     VtxSoAView soa,
0035                                     WsSoAView pws,
0036                                     float ptMin,
0037                                     float ptMax) const {
0038         auto const* quality = tracks_view.quality();
0039         using helper = TracksUtilities<TrackerTraits>;
0040
0041         for (auto idx : cms::alpakatools::uniform_elements(acc, tracks_view.nTracks())) {
0042           [[maybe_unused]] auto nHits = helper::nHits(tracks_view, idx);
0043           ALPAKA_ASSERT_ACC(nHits >= 3);
0044
0045           // initialize soa...
0046           soa[idx].idv() = -1;
0047
0048           if (reco::isTriplet(tracks_view, idx))
0049             continue;  // no triplets
0050           if (quality[idx] < ::pixelTrack::Quality::highPurity)
0051             continue;
0052
0053           auto pt = tracks_view[idx].pt();
0054
0055           if (pt < ptMin)
0056             continue;
0057
0058           // clamp pt
0059           pt = std::min<float>(pt, ptMax);
0060
0061           auto& data = pws;
0062           auto it = alpaka::atomicAdd(acc, &data.ntrks(), 1u, alpaka::hierarchy::Blocks{});
0063           data[it].itrk() = idx;
0064           data[it].zt() = reco::zip(tracks_view, idx);
0065           data[it].ezt2() = tracks_view[idx].covariance()(14);
0066           data[it].ptt2() = pt * pt;
0067         }
0068       }
0069     };
0070 // #define THREE_KERNELS
0071 #ifndef THREE_KERNELS
0072     class VertexFinderOneKernel {
0073     public:
0074       template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0075       ALPAKA_FN_ACC void operator()(const TAcc& acc,
0076                                     VtxSoAView pdata,
0077                                     WsSoAView pws,
0078                                     bool doSplit,
0079                                     int minT,      // min number of neighbours to be "seed"
0080                                     float eps,     // max absolute distance to cluster
0081                                     float errmax,  // max error to be "seed"
0082                                     float chi2max  // max normalized distance to cluster,
0083       ) const {
0084         clusterTracksByDensity(acc, pdata, pws, minT, eps, errmax, chi2max);
0085         alpaka::syncBlockThreads(acc);
0086         fitVertices(acc, pdata, pws, maxChi2ForFirstFit);
0087         alpaka::syncBlockThreads(acc);
0088         if (doSplit) {
0089           splitVertices(acc, pdata, pws, maxChi2ForSplit);
0090           alpaka::syncBlockThreads(acc);
0091           fitVertices(acc, pdata, pws, maxChi2ForFinalFit);
0092           alpaka::syncBlockThreads(acc);
0093         }
0094         sortByPt2(acc, pdata, pws);
0095       }
0096     };
0097 #else
0098     class VertexFinderKernel1 {
0099     public:
0100       template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0101       ALPAKA_FN_ACC void operator()(const TAcc& acc,
0102                                     VtxSoAView pdata,
0103                                     WsSoAView pws,
0104                                     int minT,      // min number of neighbours to be "seed"
0105                                     float eps,     // max absolute distance to cluster
0106                                     float errmax,  // max error to be "seed"
0107                                     float chi2max  // max normalized distance to cluster,
0108       ) const {
0109         clusterTracksByDensity(pdata, pws, minT, eps, errmax, chi2max);
0110         alpaka::syncBlockThreads(acc);
0111         fitVertices(pdata, pws, maxChi2ForFirstFit);
0112       }
0113     };
0114     class VertexFinderKernel2 {
0115     public:
0116       template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
0117       ALPAKA_FN_ACC void operator()(const TAcc& acc, VtxSoAView pdata, WsSoAView pws) const {
0118         fitVertices(pdata, pws, maxChi2ForFinalFit);
0119         alpaka::syncBlockThreads(acc);
0120         sortByPt2(pdata, pws);
0121       }
0122     };
0123 #endif
0124
0125     template <typename TrackerTraits>
0126     ZVertexSoACollection Producer<TrackerTraits>::makeAsync(Queue& queue,
0127                                                             const reco::TrackSoAConstView<TrackerTraits>& tracks_view,
0128                                                             float ptMin,
0129                                                             float ptMax) const {
0130 #ifdef PIXVERTEX_DEBUG_PRODUCE
0131       std::cout << "producing Vertices on GPU" << std::endl;
0132 #endif  // PIXVERTEX_DEBUG_PRODUCE
0133       ZVertexSoACollection vertices(queue);
0134
0135       auto soa = vertices.view();
0136
0137       auto ws_d = PixelVertexWorkSpaceSoADevice(::zVertex::MAXTRACKS, queue);
0138
0139       // Initialize
0140       const auto initWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1);
0141       alpaka::exec<Acc1D>(queue, initWorkDiv, Init{}, soa, ws_d.view());
0142
0143       // Load Tracks
0144       const uint32_t blockSize = 128;
0145       const uint32_t numberOfBlocks =
0146           cms::alpakatools::divide_up_by(tracks_view.metadata().size() + blockSize - 1, blockSize);
0147       const auto loadTracksWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(numberOfBlocks, blockSize);
0148       alpaka::exec<Acc1D>(
0149           queue, loadTracksWorkDiv, LoadTracks<TrackerTraits>{}, tracks_view, soa, ws_d.view(), ptMin, ptMax);
0150
0151       // Running too many thread lead to problems when printf is enabled.
0152       const auto finderSorterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1024 - 128);
0153       const auto splitterFitterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1024, 128);
0154
0155       if (oneKernel_) {
0156         // implemented only for density clustesrs
0157 #ifndef THREE_KERNELS
0158         alpaka::exec<Acc1D>(queue,
0159                             finderSorterWorkDiv,
0160                             VertexFinderOneKernel{},
0161                             soa,
0162                             ws_d.view(),
0163                             doSplitting_,
0164                             minT,
0165                             eps,
0166                             errmax,
0167                             chi2max);
0168 #else
0169         alpaka::exec<Acc1D>(
0170             queue, finderSorterWorkDiv, VertexFinderOneKernel{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
0171
0172         // one block per vertex...
0173         if (doSplitting_)
0174           alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, soa, ws_d.view(), maxChi2ForSplit);
0175         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv{}, soa, ws_d.view());
0176 #endif
0177       } else {  // five kernels
0178         if (useDensity_) {
0179           alpaka::exec<Acc1D>(
0180               queue, finderSorterWorkDiv, ClusterTracksByDensityKernel{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
0181
0182         } else if (useDBSCAN_) {
0183           alpaka::exec<Acc1D>(
0184               queue, finderSorterWorkDiv, ClusterTracksDBSCAN{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
0185         } else if (useIterative_) {
0186           alpaka::exec<Acc1D>(
0187               queue, finderSorterWorkDiv, ClusterTracksIterative{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
0188         }
0189         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, soa, ws_d.view(), maxChi2ForFirstFit);
0190
0191         // one block per vertex...
0192         if (doSplitting_) {
0193           alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, soa, ws_d.view(), maxChi2ForSplit);
0194
0195           alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, soa, ws_d.view(), maxChi2ForFinalFit);
0196         }
0197         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, SortByPt2Kernel{}, soa, ws_d.view());
0198       }
0199
0200       return vertices;
0201     }
0202
0203     template class Producer<pixelTopology::Phase1>;
0204     template class Producer<pixelTopology::Phase2>;
0205     template class Producer<pixelTopology::HIonPhase1>;
0206   }  // namespace vertexFinder
0207 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE