CMS 3D CMS Logo

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