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 "RecoVertex/PixelVertexFinding/interface/PixelVertexWorkSpaceLayout.h"
0007 #include "RecoVertex/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       ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0032                                     ::reco::TrackSoAConstView tracks_view,
0033                                     VtxSoAView data,
0034                                     TrkSoAView trkdata,
0035                                     WsSoAView ws,
0036                                     float ptMin,
0037                                     float ptMax) const {
0038         auto const* quality = tracks_view.quality();
0039
0040         for (auto idx : cms::alpakatools::uniform_elements(acc, tracks_view.nTracks())) {
0041           [[maybe_unused]] auto nHits = reco::nHits(tracks_view, idx);
0042           ALPAKA_ASSERT_ACC(nHits >= 3);
0043
0044           // initialize the track data
0045           trkdata[idx].idv() = -1;
0046
0047           // do not use triplets
0048           if (reco::isTriplet(tracks_view, idx))
0049             continue;
0050
0051           // use only "high purity" track
0052           if (quality[idx] < ::pixelTrack::Quality::highPurity)
0053             continue;
0054
0055           auto pt = tracks_view[idx].pt();
0056           // pT min cut
0057           if (pt < ptMin)
0058             continue;
0059
0060           // clamp pT to the pTmax
0061           pt = std::min<float>(pt, ptMax);
0062
0063           // load the track data into the workspace
0064           auto it = alpaka::atomicAdd(acc, &ws.ntrks(), 1u, alpaka::hierarchy::Blocks{});
0065           ws[it].itrk() = idx;
0066           ws[it].zt() = reco::zip(tracks_view, idx);
0067           ws[it].ezt2() = tracks_view[idx].covariance()(14);
0068           ws[it].ptt2() = pt * pt;
0069         }
0070       }
0071     };
0072
0073 // #define THREE_KERNELS
0074 #ifndef THREE_KERNELS
0075     class VertexFinderOneKernel {
0076     public:
0077       ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0078                                     VtxSoAView data,
0079                                     TrkSoAView trkdata,
0080                                     WsSoAView ws,
0081                                     bool doSplit,
0082                                     int minT,      // min number of neighbours to be "seed"
0083                                     float eps,     // max absolute distance to cluster
0084                                     float errmax,  // max error to be "seed"
0085                                     float chi2max  // max normalized distance to cluster,
0086       ) const {
0087         clusterTracksByDensity(acc, data, trkdata, ws, minT, eps, errmax, chi2max);
0088         alpaka::syncBlockThreads(acc);
0089         fitVertices(acc, data, trkdata, ws, maxChi2ForFirstFit);
0090         alpaka::syncBlockThreads(acc);
0091         if (doSplit) {
0092           splitVertices(acc, data, trkdata, ws, maxChi2ForSplit);
0093           alpaka::syncBlockThreads(acc);
0094           fitVertices(acc, data, trkdata, ws, maxChi2ForFinalFit);
0095           alpaka::syncBlockThreads(acc);
0096         }
0097         sortByPt2(acc, data, trkdata, ws);
0098       }
0099     };
0100 #else
0101     class VertexFinderKernel1 {
0102     public:
0103       ALPAKA_FN_ACC void operator()(Acc1D const& acc,
0104                                     VtxSoAView data,
0105                                     WsSoAView ws,
0106                                     int minT,      // min number of neighbours to be "seed"
0107                                     float eps,     // max absolute distance to cluster
0108                                     float errmax,  // max error to be "seed"
0109                                     float chi2max  // max normalized distance to cluster,
0110       ) const {
0111         clusterTracksByDensity(acc, data, ws, minT, eps, errmax, chi2max);
0112         alpaka::syncBlockThreads(acc);
0113         fitVertices(acc, data, ws, maxChi2ForFirstFit);
0114       }
0115     };
0116
0117     class VertexFinderKernel2 {
0118     public:
0119       ALPAKA_FN_ACC void operator()(Acc1D const& acc, VtxSoAView data, WsSoAView ws) const {
0120         fitVertices(acc, data, ws, maxChi2ForFinalFit);
0121         alpaka::syncBlockThreads(acc);
0122         sortByPt2(data, ws);
0123       }
0124     };
0125 #endif
0126
0127     template <typename TrackerTraits>
0128     ZVertexSoACollection Producer<TrackerTraits>::makeAsync(
0129         Queue& queue, ::reco::TrackSoAConstView const& tracks_view, int maxVertices, float ptMin, float ptMax) const {
0130 #ifdef PIXVERTEX_DEBUG_PRODUCE
0131       std::cout << "producing Vertices on GPU" << std::endl;
0132 #endif  // PIXVERTEX_DEBUG_PRODUCE
0133       const auto maxTracks = tracks_view.metadata().size();
0134       ZVertexSoACollection vertices({{maxVertices, maxTracks}}, queue);
0135       auto data = vertices.view();
0136       auto trkdata = vertices.view<reco::ZVertexTracksSoA>();
0137
0138       PixelVertexWorkSpaceSoADevice workspace(maxTracks, queue);
0139       auto ws = workspace.view();
0140
0141       // Initialize
0142       const auto initWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1);
0143       alpaka::exec<Acc1D>(queue, initWorkDiv, Init{}, data, ws);
0144
0145       // Load Tracks
0146       const uint32_t blockSize = 128;
0147       const uint32_t numberOfBlocks =
0148           cms::alpakatools::divide_up_by(tracks_view.metadata().size() + blockSize - 1, blockSize);
0149       const auto loadTracksWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(numberOfBlocks, blockSize);
0150       alpaka::exec<Acc1D>(
0151           queue, loadTracksWorkDiv, LoadTracks<TrackerTraits>{}, tracks_view, data, trkdata, ws, ptMin, ptMax);
0152
0153       // Running too many thread lead to problems when printf is enabled.
0154       const auto finderSorterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1024 - 128);
0155       const auto splitterFitterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1024, 128);
0156
0157       if (oneKernel_) {
0158         // implemented only for density clustesrs
0159 #ifndef THREE_KERNELS
0160         alpaka::exec<Acc1D>(queue,
0161                             finderSorterWorkDiv,
0162                             VertexFinderOneKernel{},
0163                             data,
0164                             trkdata,
0165                             ws,
0166                             doSplitting_,
0167                             minT,
0168                             eps,
0169                             errmax,
0170                             chi2max);
0171 #else
0172         alpaka::exec<Acc1D>(
0173             queue, finderSorterWorkDiv, VertexFinderOneKernel{}, data, trkdata, ws, minT, eps, errmax, chi2max);
0174
0175         // one block per vertex...
0176         if (doSplitting_)
0177           alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, data, trkdata, ws, maxChi2ForSplit);
0178         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv{}, data, ws);
0179 #endif
0180       } else {  // five kernels
0181         if (useDensity_) {
0182           alpaka::exec<Acc1D>(
0183               queue, finderSorterWorkDiv, ClusterTracksByDensityKernel{}, data, trkdata, ws, minT, eps, errmax, chi2max);
0184
0185         } else if (useDBSCAN_) {
0186           alpaka::exec<Acc1D>(
0187               queue, finderSorterWorkDiv, ClusterTracksDBSCAN{}, data, trkdata, ws, minT, eps, errmax, chi2max);
0188         } else if (useIterative_) {
0189           alpaka::exec<Acc1D>(
0190               queue, finderSorterWorkDiv, ClusterTracksIterative{}, data, trkdata, ws, minT, eps, errmax, chi2max);
0191         }
0192         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, data, trkdata, ws, maxChi2ForFirstFit);
0193
0194         // one block per vertex...
0195         if (doSplitting_) {
0196           alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, data, trkdata, ws, maxChi2ForSplit);
0197
0198           alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, data, trkdata, ws, maxChi2ForFinalFit);
0199         }
0200         alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, SortByPt2Kernel{}, data, trkdata, ws);
0201       }
0202
0203       return vertices;
0204     }
0205
0206     template class Producer<pixelTopology::Phase1>;
0207     template class Producer<pixelTopology::Phase2>;
0208     template class Producer<pixelTopology::HIonPhase1>;
0209   }  // namespace vertexFinder
0210 }  // namespace ALPAKA_ACCELERATOR_NAMESPACE