Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:28:40

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