CMS 3D CMS Logo

vertexFinder.dev.cc
Go to the documentation of this file.
1 #include <alpaka/alpaka.hpp>
2 
9 
10 #include "vertexFinder.h"
11 #include "clusterTracksDBSCAN.h"
12 #include "clusterTracksIterative.h"
13 #include "clusterTracksByDensity.h"
14 #include "fitVertices.h"
15 #include "sortByPt2.h"
16 #include "splitVertices.h"
17 
18 #undef PIXVERTEX_DEBUG_PRODUCE
20  namespace vertexFinder {
21  using namespace cms::alpakatools;
22  // reject outlier tracks that contribute more than this to the chi2 of the vertex fit
24  constexpr float maxChi2ForFinalFit = 5000.f;
25 
26  // split vertices with a chi2/NDoF greater than this
28 
29  template <typename TrackerTraits>
30  class LoadTracks {
31  public:
32  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
33  ALPAKA_FN_ACC void operator()(const TAcc& acc,
36  WsSoAView pws,
37  float ptMin,
38  float ptMax) const {
39  auto const* quality = tracks_view.quality();
41 
42  for (auto idx : cms::alpakatools::uniform_elements(acc, tracks_view.nTracks())) {
43  [[maybe_unused]] auto nHits = helper::nHits(tracks_view, idx);
45 
46  // initialize soa...
47  soa[idx].idv() = -1;
48 
50  continue; // no triplets
52  continue;
53 
54  auto pt = tracks_view[idx].pt();
55 
56  if (pt < ptMin)
57  continue;
58 
59  // clamp pt
60  pt = std::min<float>(pt, ptMax);
61 
62  auto& data = pws;
63  auto it = alpaka::atomicAdd(acc, &data.ntrks(), 1u, alpaka::hierarchy::Blocks{});
64  data[it].itrk() = idx;
65  data[it].zt() = reco::zip(tracks_view, idx);
66  data[it].ezt2() = tracks_view[idx].covariance()(14);
67  data[it].ptt2() = pt * pt;
68  }
69  }
70  };
71 // #define THREE_KERNELS
72 #ifndef THREE_KERNELS
74  public:
75  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
76  ALPAKA_FN_ACC void operator()(const TAcc& acc,
78  WsSoAView pws,
79  bool doSplit,
80  int minT, // min number of neighbours to be "seed"
81  float eps, // max absolute distance to cluster
82  float errmax, // max error to be "seed"
83  float chi2max // max normalized distance to cluster,
84  ) const {
85  clusterTracksByDensity(acc, pdata, pws, minT, eps, errmax, chi2max);
86  alpaka::syncBlockThreads(acc);
88  alpaka::syncBlockThreads(acc);
89  if (doSplit) {
91  alpaka::syncBlockThreads(acc);
93  alpaka::syncBlockThreads(acc);
94  }
95  sortByPt2(acc, pdata, pws);
96  }
97  };
98 #else
99  class VertexFinderKernel1 {
100  public:
101  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
102  ALPAKA_FN_ACC void operator()(const TAcc& acc,
104  WsSoAView pws,
105  int minT, // min number of neighbours to be "seed"
106  float eps, // max absolute distance to cluster
107  float errmax, // max error to be "seed"
108  float chi2max // max normalized distance to cluster,
109  ) const {
110  clusterTracksByDensity(pdata, pws, minT, eps, errmax, chi2max);
111  alpaka::syncBlockThreads(acc);
113  }
114  };
115  class VertexFinderKernel2 {
116  public:
117  template <typename TAcc, typename = std::enable_if_t<alpaka::isAccelerator<TAcc>>>
118  ALPAKA_FN_ACC void operator()(const TAcc& acc, VtxSoAView pdata, WsSoAView pws) const {
120  alpaka::syncBlockThreads(acc);
121  sortByPt2(pdata, pws);
122  }
123  };
124 #endif
125 
126  template <typename TrackerTraits>
129  float ptMin,
130  float ptMax) const {
131 #ifdef PIXVERTEX_DEBUG_PRODUCE
132  std::cout << "producing Vertices on GPU" << std::endl;
133 #endif // PIXVERTEX_DEBUG_PRODUCE
135 
136  auto soa = vertices.view();
137 
139 
140  // Initialize
141  const auto initWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1);
142  alpaka::exec<Acc1D>(queue, initWorkDiv, Init{}, soa, ws_d.view());
143 
144  // Load Tracks
145  const uint32_t blockSize = 128;
146  const uint32_t numberOfBlocks =
147  cms::alpakatools::divide_up_by(tracks_view.metadata().size() + blockSize - 1, blockSize);
148  const auto loadTracksWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(numberOfBlocks, blockSize);
149  alpaka::exec<Acc1D>(
150  queue, loadTracksWorkDiv, LoadTracks<TrackerTraits>{}, tracks_view, soa, ws_d.view(), ptMin, ptMax);
151 
152  // Running too many thread lead to problems when printf is enabled.
153  const auto finderSorterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1, 1024 - 128);
154  const auto splitterFitterWorkDiv = cms::alpakatools::make_workdiv<Acc1D>(1024, 128);
155 
156  if (oneKernel_) {
157  // implemented only for density clustesrs
158 #ifndef THREE_KERNELS
159  alpaka::exec<Acc1D>(queue,
160  finderSorterWorkDiv,
161  VertexFinderOneKernel{},
162  soa,
163  ws_d.view(),
164  doSplitting_,
165  minT,
166  eps,
167  errmax,
168  chi2max);
169 #else
170  alpaka::exec<Acc1D>(
171  queue, finderSorterWorkDiv, VertexFinderOneKernel{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
172 
173  // one block per vertex...
174  if (doSplitting_)
175  alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, soa, ws_d.view(), maxChi2ForSplit);
176  alpaka::exec<Acc1D>(queue, finderSorterWorkDiv{}, soa, ws_d.view());
177 #endif
178  } else { // five kernels
179  if (useDensity_) {
180  alpaka::exec<Acc1D>(
181  queue, finderSorterWorkDiv, ClusterTracksByDensityKernel{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
182 
183  } else if (useDBSCAN_) {
184  alpaka::exec<Acc1D>(
185  queue, finderSorterWorkDiv, ClusterTracksDBSCAN{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
186  } else if (useIterative_) {
187  alpaka::exec<Acc1D>(
188  queue, finderSorterWorkDiv, ClusterTracksIterative{}, soa, ws_d.view(), minT, eps, errmax, chi2max);
189  }
190  alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, soa, ws_d.view(), maxChi2ForFirstFit);
191 
192  // one block per vertex...
193  if (doSplitting_) {
194  alpaka::exec<Acc1D>(queue, splitterFitterWorkDiv, SplitVerticesKernel{}, soa, ws_d.view(), maxChi2ForSplit);
195 
196  alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, FitVerticesKernel{}, soa, ws_d.view(), maxChi2ForFinalFit);
197  }
198  alpaka::exec<Acc1D>(queue, finderSorterWorkDiv, SortByPt2Kernel{}, soa, ws_d.view());
199  }
200 
201  return vertices;
202  }
203 
204  template class Producer<pixelTopology::Phase1>;
205  template class Producer<pixelTopology::Phase2>;
206  template class Producer<pixelTopology::HIonPhase1>;
207  } // namespace vertexFinder
208 } // namespace ALPAKA_ACCELERATOR_NAMESPACE
ALPAKA_FN_ACC auto uniform_elements(TAcc const &acc, TArgs... args)
Definition: workdivision.h:311
Definition: helper.py:1
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView & pdata
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView WsSoAView int float eps
constexpr Idx divide_up_by(Idx value, Idx divisor)
Definition: workdivision.h:20
ALPAKA_FN_ACC void operator()(const TAcc &acc, reco::TrackSoAConstView< TrackerTraits > tracks_view, VtxSoAView soa, WsSoAView pws, float ptMin, float ptMax) const
std::conditional_t< std::is_same_v< Device, alpaka::DevCpu >, ZVertexHost, ZVertexDevice< Device > > ZVertexSoACollection
uint32_t const *__restrict__ TkSoAView< TrackerTraits > tracks_view
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView WsSoAView int float float errmax
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr float zip(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:90
constexpr uint32_t MAXTRACKS
fitVertices(pdata, pws, maxChi2ForFirstFit)
constexpr float ptMin
string quality
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView WsSoAView & pws
std::vector< Block > Blocks
Definition: Block.h:99
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView WsSoAView int minT
splitVertices(pdata, pws, maxChi2ForSplit)
ALPAKA_FN_HOST_ACC ALPAKA_FN_INLINE constexpr bool isTriplet(ConstView const &tracks, int32_t i)
Definition: TracksSoA.h:95
ALPAKA_FN_ACC void operator()(const TAcc &acc, VtxSoAView pdata, WsSoAView pws, bool doSplit, int minT, float eps, float errmax, float chi2max) const
ALPAKA_FN_ACC ALPAKA_FN_INLINE void sortByPt2(const TAcc &acc, VtxSoAView &data, WsSoAView &ws)
Definition: sortByPt2.h:26
::vertexFinder::PixelVertexWorkSpaceSoAView WsSoAView
typename reco::TrackSoA< TrackerTraits >::template Layout<>::ConstView TrackSoAConstView
Definition: TracksSoA.h:47
PortableCollection<::vertexFinder::PixelVertexWSSoALayout<> > PixelVertexWorkSpaceSoADevice
ALPAKA_FN_ACC ALPAKA_FN_INLINE void VtxSoAView WsSoAView int float float float chi2max
ZVertexSoACollection makeAsync(Queue &queue, const TkSoAConstView &tracks_view, float ptMin, float ptMax) const
TupleMultiplicity< TrackerTraits > const *__restrict__ uint32_t nHits
T1 atomicAdd(T1 *a, T2 b)
Definition: cudaCompat.h:61