Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SeedFinder.ipp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SeedFinder.ipp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2023 CERN for the benefit of the Acts project
4 //
5 // This Source Code Form is subject to the terms of the Mozilla Public
6 // License, v. 2.0. If a copy of the MPL was not distributed with this
7 // file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 
9 #include "Acts/Seeding/CandidatesForMiddleSp.hpp"
10 
11 #include <algorithm>
12 #include <cmath>
13 #include <numeric>
14 #include <type_traits>
15 
16 namespace Acts {
17 
18 template <typename external_spacepoint_t>
19 SeedFinder<external_spacepoint_t, Acts::Cuda>::SeedFinder(
20  const Acts::SeedFinderConfig<external_spacepoint_t>& config,
21  const Acts::SeedFinderOptions& options)
22  : m_config(config), m_options(options) {
23  if (not m_config.isInInternalUnits)
24  throw std::runtime_error(
25  "SeedFinderConfig not in ACTS internal units in "
26  "Cuda/Seeding/SeedFinder");
27  if (not m_options.isInInternalUnits)
28  throw std::runtime_error(
29  "SeedFinderOptions not in ACTS internal units in "
30  "Cuda/Seeding/SeedFinder");
31  if (std::isnan(m_config.deltaRMaxTopSP)) {
32  throw std::runtime_error("Value of deltaRMaxTopSP was not initialised");
33  }
34  if (std::isnan(m_config.deltaRMinTopSP)) {
35  throw std::runtime_error("Value of deltaRMinTopSP was not initialised");
36  }
37  if (std::isnan(m_config.deltaRMaxBottomSP)) {
38  throw std::runtime_error("Value of deltaRMaxBottomSP was not initialised");
39  }
40  if (std::isnan(m_config.deltaRMinBottomSP)) {
41  throw std::runtime_error("Value of deltaRMinBottomSP was not initialised");
42  }
43 }
44 
45 // CUDA seed finding
46 template <typename external_spacepoint_t>
47 template <typename sp_range_t>
48 std::vector<Seed<external_spacepoint_t>>
49 SeedFinder<external_spacepoint_t, Acts::Cuda>::createSeedsForGroup(
50  Acts::SpacePointData& spacePointData,
51  Acts::SpacePointGrid<external_spacepoint_t>& grid,
52  const sp_range_t& bottomSPs, const std::size_t middleSPs,
53  const sp_range_t& topSPs) const {
54  std::vector<Seed<external_spacepoint_t>> outputVec;
55 
56  // Get SeedFinderConfig values
57  CudaScalar<float> deltaRMin_cuda(&m_config.deltaRMin);
58  CudaScalar<float> deltaRMax_cuda(&m_config.deltaRMax);
59  CudaScalar<float> cotThetaMax_cuda(&m_config.cotThetaMax);
60  CudaScalar<float> collisionRegionMin_cuda(&m_config.collisionRegionMin);
61  CudaScalar<float> collisionRegionMax_cuda(&m_config.collisionRegionMax);
62  CudaScalar<float> maxScatteringAngle2_cuda(&m_config.maxScatteringAngle2);
63  CudaScalar<float> sigmaScattering_cuda(&m_config.sigmaScattering);
64  CudaScalar<float> minHelixDiameter2_cuda(&m_options.minHelixDiameter2);
65  CudaScalar<float> pT2perRadius_cuda(&m_options.pT2perRadius);
66  CudaScalar<float> impactMax_cuda(&m_config.impactMax);
67  const auto seedFilterConfig = m_config.seedFilter->getSeedFilterConfig();
68  CudaScalar<float> deltaInvHelixDiameter_cuda(
69  &seedFilterConfig.deltaInvHelixDiameter);
70  CudaScalar<float> impactWeightFactor_cuda(
71  &seedFilterConfig.impactWeightFactor);
72  CudaScalar<float> filterDeltaRMin_cuda(&seedFilterConfig.deltaRMin);
73  CudaScalar<float> compatSeedWeight_cuda(&seedFilterConfig.compatSeedWeight);
74  CudaScalar<size_t> compatSeedLimit_cuda(&seedFilterConfig.compatSeedLimit);
75  CpuScalar<size_t> compatSeedLimit_cpu(&compatSeedLimit_cuda);
76  //---------------------------------
77  // Algorithm 0. Matrix Flattening
78  //---------------------------------
79 
80  std::vector<Acts::InternalSpacePoint<external_spacepoint_t>*> middleSPvec;
81  std::vector<Acts::InternalSpacePoint<external_spacepoint_t>*> bottomSPvec;
82  std::vector<Acts::InternalSpacePoint<external_spacepoint_t>*> topSPvec;
83 
84  // Get the size of spacepoints
85  int nSpM(0);
86  int nSpB(0);
87  int nSpT(0);
88 
89  {
90  auto& sp_collection = grid.at(middleSPs);
91  for (auto& sp : sp_collection) {
92  nSpM++;
93  middleSPvec.push_back(sp.get());
94  }
95  }
96  for (auto idx : bottomSPs) {
97  auto& sp_collection = grid.at(idx);
98  for (auto& sp : sp_collection) {
99  nSpB++;
100  bottomSPvec.push_back(sp.get());
101  }
102  }
103  for (std::size_t idx : topSPs) {
104  auto& sp_collection = grid.at(idx);
105  for (auto& sp : sp_collection) {
106  nSpT++;
107  topSPvec.push_back(sp.get());
108  }
109  }
110 
111  CudaScalar<int> nSpM_cuda(&nSpM);
112  CudaScalar<int> nSpB_cuda(&nSpB);
113  CudaScalar<int> nSpT_cuda(&nSpT);
114 
115  if (nSpM == 0 || nSpB == 0 || nSpT == 0)
116  return outputVec;
117 
118  // Matrix flattening
119  CpuMatrix<float> spMmat_cpu(nSpM, 6); // x y z r varR varZ
120  CpuMatrix<float> spBmat_cpu(nSpB, 6);
121  CpuMatrix<float> spTmat_cpu(nSpT, 6);
122 
123  auto fillMatrix = [](auto& mat, auto& id, auto sp) {
124  mat.set(id, 0, sp->x());
125  mat.set(id, 1, sp->y());
126  mat.set(id, 2, sp->z());
127  mat.set(id, 3, sp->radius());
128  mat.set(id, 4, sp->varianceR());
129  mat.set(id, 5, sp->varianceZ());
130  id++;
131  };
132 
133  int mIdx(0);
134  for (auto sp : middleSPvec) {
135  fillMatrix(spMmat_cpu, mIdx, sp);
136  }
137  int bIdx(0);
138  for (auto sp : bottomSPvec) {
139  fillMatrix(spBmat_cpu, bIdx, sp);
140  }
141  int tIdx(0);
142  for (auto sp : topSPvec) {
143  fillMatrix(spTmat_cpu, tIdx, sp);
144  }
145 
146  CudaMatrix<float> spMmat_cuda(nSpM, 6, &spMmat_cpu);
147  CudaMatrix<float> spBmat_cuda(nSpB, 6, &spBmat_cpu);
148  CudaMatrix<float> spTmat_cuda(nSpT, 6, &spTmat_cpu);
149  //------------------------------------
150  // Algorithm 1. Doublet Search (DS)
151  //------------------------------------
152 
153  CudaScalar<int> nSpMcomp_cuda(new int(0));
154  CudaScalar<int> nSpBcompPerSpMMax_cuda(new int(0));
155  CudaScalar<int> nSpTcompPerSpMMax_cuda(new int(0));
156  CudaVector<int> nSpBcompPerSpM_cuda(nSpM);
157  nSpBcompPerSpM_cuda.zeros();
158  CudaVector<int> nSpTcompPerSpM_cuda(nSpM);
159  nSpTcompPerSpM_cuda.zeros();
160  CudaVector<int> McompIndex_cuda(nSpM);
161  CudaMatrix<int> BcompIndex_cuda(nSpB, nSpM);
162  CudaMatrix<int> TcompIndex_cuda(nSpT, nSpM);
163  CudaMatrix<int> tmpBcompIndex_cuda(nSpB, nSpM);
164  CudaMatrix<int> tmpTcompIndex_cuda(nSpT, nSpM);
165 
166  dim3 DS_BlockSize = m_config.maxBlockSize;
167  dim3 DS_GridSize(nSpM, 1, 1);
168 
169  searchDoublet(DS_GridSize, DS_BlockSize, nSpM_cuda.get(), spMmat_cuda.get(),
170  nSpB_cuda.get(), spBmat_cuda.get(), nSpT_cuda.get(),
171  spTmat_cuda.get(), deltaRMin_cuda.get(), deltaRMax_cuda.get(),
172  cotThetaMax_cuda.get(), collisionRegionMin_cuda.get(),
173  collisionRegionMax_cuda.get(), nSpMcomp_cuda.get(),
174  nSpBcompPerSpMMax_cuda.get(), nSpTcompPerSpMMax_cuda.get(),
175  nSpBcompPerSpM_cuda.get(), nSpTcompPerSpM_cuda.get(),
176  McompIndex_cuda.get(), BcompIndex_cuda.get(),
177  tmpBcompIndex_cuda.get(), TcompIndex_cuda.get(),
178  tmpTcompIndex_cuda.get());
179 
180  CpuScalar<int> nSpMcomp_cpu(&nSpMcomp_cuda);
181  CpuScalar<int> nSpBcompPerSpMMax_cpu(&nSpBcompPerSpMMax_cuda);
182  CpuScalar<int> nSpTcompPerSpMMax_cpu(&nSpTcompPerSpMMax_cuda);
183  CpuVector<int> nSpBcompPerSpM_cpu(nSpM, &nSpBcompPerSpM_cuda);
184  CpuVector<int> nSpTcompPerSpM_cpu(nSpM, &nSpTcompPerSpM_cuda);
185  CpuVector<int> McompIndex_cpu(nSpM, &McompIndex_cuda);
186 
187  //--------------------------------------
188  // Algorithm 2. Transform coordinate
189  //--------------------------------------
190 
191  CudaMatrix<float> spMcompMat_cuda(*nSpMcomp_cpu.get(), 6);
192  CudaMatrix<float> spBcompMatPerSpM_cuda(*nSpBcompPerSpMMax_cpu.get(),
193  (*nSpMcomp_cpu.get()) * 6);
194  CudaMatrix<float> spTcompMatPerSpM_cuda(*nSpTcompPerSpMMax_cpu.get(),
195  (*nSpMcomp_cpu.get()) * 6);
196  CudaMatrix<float> circBcompMatPerSpM_cuda(*nSpBcompPerSpMMax_cpu.get(),
197  (*nSpMcomp_cpu.get()) * 6);
198  CudaMatrix<float> circTcompMatPerSpM_cuda(*nSpTcompPerSpMMax_cpu.get(),
199  (*nSpMcomp_cpu.get()) * 6);
200 
201  dim3 TC_GridSize(*nSpMcomp_cpu.get(), 1, 1);
202  dim3 TC_BlockSize = m_config.maxBlockSize;
203 
204  transformCoordinate(
205  TC_GridSize, TC_BlockSize, nSpM_cuda.get(), spMmat_cuda.get(),
206  McompIndex_cuda.get(), nSpB_cuda.get(), spBmat_cuda.get(),
207  nSpBcompPerSpMMax_cuda.get(), BcompIndex_cuda.get(), nSpT_cuda.get(),
208  spTmat_cuda.get(), nSpTcompPerSpMMax_cuda.get(), TcompIndex_cuda.get(),
209  spMcompMat_cuda.get(), spBcompMatPerSpM_cuda.get(),
210  circBcompMatPerSpM_cuda.get(), spTcompMatPerSpM_cuda.get(),
211  circTcompMatPerSpM_cuda.get());
212 
213  //------------------------------------------------------
214  // Algorithm 3. Triplet Search (TS) & Seed filtering
215  //------------------------------------------------------
216 
217  const int nTrplPerSpMLimit =
218  m_config.nAvgTrplPerSpBLimit * (*nSpBcompPerSpMMax_cpu.get());
219  CudaScalar<int> nTrplPerSpMLimit_cuda(&nTrplPerSpMLimit);
220 
221  CudaScalar<int> nTrplPerSpBLimit_cuda(&m_config.nTrplPerSpBLimit);
222  CpuScalar<int> nTrplPerSpBLimit_cpu(
223  &nTrplPerSpBLimit_cuda); // need to be USM
224 
225  CudaVector<int> nTrplPerSpM_cuda(*nSpMcomp_cpu.get());
226  nTrplPerSpM_cuda.zeros();
227  CudaMatrix<Triplet> TripletsPerSpM_cuda(nTrplPerSpMLimit,
228  *nSpMcomp_cpu.get());
229  CpuVector<int> nTrplPerSpM_cpu(*nSpMcomp_cpu.get(), true);
230  nTrplPerSpM_cpu.zeros();
231  CpuMatrix<Triplet> TripletsPerSpM_cpu(nTrplPerSpMLimit, *nSpMcomp_cpu.get(),
232  true);
233  cudaStream_t cuStream;
234  ACTS_CUDA_ERROR_CHECK(cudaStreamCreate(&cuStream));
235 
236  for (int i_m = 0; i_m <= *nSpMcomp_cpu.get(); i_m++) {
237  cudaStreamSynchronize(cuStream);
238 
239  // Search Triplet
240  if (i_m < *nSpMcomp_cpu.get()) {
241  int mIndex = *McompIndex_cpu.get(i_m);
242  int* nSpBcompPerSpM = nSpBcompPerSpM_cpu.get(mIndex);
243  int* nSpTcompPerSpM = nSpTcompPerSpM_cpu.get(mIndex);
244 
245  dim3 TS_GridSize(*nSpBcompPerSpM, 1, 1);
246  dim3 TS_BlockSize =
247  dim3(fmin(m_config.maxBlockSize, *nSpTcompPerSpM), 1, 1);
248 
249  searchTriplet(
250  TS_GridSize, TS_BlockSize, nSpTcompPerSpM_cpu.get(mIndex),
251  nSpTcompPerSpM_cuda.get(mIndex), nSpMcomp_cuda.get(),
252  spMcompMat_cuda.get(i_m, 0), nSpBcompPerSpMMax_cuda.get(),
253  BcompIndex_cuda.get(0, i_m), circBcompMatPerSpM_cuda.get(0, 6 * i_m),
254  nSpTcompPerSpMMax_cuda.get(), TcompIndex_cuda.get(0, i_m),
255  spTcompMatPerSpM_cuda.get(0, 6 * i_m),
256  circTcompMatPerSpM_cuda.get(0, 6 * i_m),
257  // Seed finder config
258  maxScatteringAngle2_cuda.get(), sigmaScattering_cuda.get(),
259  minHelixDiameter2_cuda.get(), pT2perRadius_cuda.get(),
260  impactMax_cuda.get(), nTrplPerSpMLimit_cuda.get(),
261  nTrplPerSpBLimit_cpu.get(), nTrplPerSpBLimit_cuda.get(),
262  deltaInvHelixDiameter_cuda.get(), impactWeightFactor_cuda.get(),
263  filterDeltaRMin_cuda.get(), compatSeedWeight_cuda.get(),
264  compatSeedLimit_cpu.get(), compatSeedLimit_cuda.get(),
265  // output
266  nTrplPerSpM_cuda.get(i_m), TripletsPerSpM_cuda.get(0, i_m),
267  &cuStream);
268  nTrplPerSpM_cpu.copyD2H(nTrplPerSpM_cuda.get(i_m), 1, i_m, &cuStream);
269 
270  TripletsPerSpM_cpu.copyD2H(TripletsPerSpM_cuda.get(0, i_m),
271  nTrplPerSpMLimit, nTrplPerSpMLimit * i_m,
272  &cuStream);
273  }
274 
275  if (i_m > 0) {
276  const auto m_experimentCuts = m_config.seedFilter->getExperimentCuts();
277  std::vector<typename CandidatesForMiddleSp<
278  const InternalSpacePoint<external_spacepoint_t>>::value_type>
279  candidates;
280 
281  for (int i = 0; i < *nTrplPerSpM_cpu.get(i_m - 1); i++) {
282  auto& triplet = *TripletsPerSpM_cpu.get(i, i_m - 1);
283  int mIndex = *McompIndex_cpu.get(i_m - 1);
284  int bIndex = triplet.bIndex;
285  int tIndex = triplet.tIndex;
286 
287  auto& bottomSP = *bottomSPvec[bIndex];
288  auto& middleSP = *middleSPvec[mIndex];
289  auto& topSP = *topSPvec[tIndex];
290  if (m_experimentCuts != nullptr) {
291  // add detector specific considerations on the seed weight
292  triplet.weight +=
293  m_experimentCuts->seedWeight(bottomSP, middleSP, topSP);
294  // discard seeds according to detector specific cuts (e.g.: weight)
295  if (!m_experimentCuts->singleSeedCut(triplet.weight, bottomSP,
296  middleSP, topSP)) {
297  continue;
298  }
299  }
300 
301  float Zob = 0; // It is not used in the seed filter but needs to be
302  // fixed anyway...
303 
304  candidates.emplace_back(bottomSP, middleSP, topSP, triplet.weight, Zob,
305  false);
306  }
307 
308  std::sort(candidates.begin(), candidates.end(),
309  CandidatesForMiddleSp<const InternalSpacePoint<
310  external_spacepoint_t>>::descendingByQuality);
311  std::size_t numQualitySeeds = 0; // not used but needs to be fixed
312  m_config.seedFilter->filterSeeds_1SpFixed(spacePointData, candidates,
313  numQualitySeeds,
314  std::back_inserter(outputVec));
315  }
316  }
317  ACTS_CUDA_ERROR_CHECK(cudaStreamDestroy(cuStream));
318  return outputVec;
319 }
320 } // namespace Acts