EIC 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) 2020 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 <algorithm>
10 #include <cmath>
11 #include <numeric>
12 #include <type_traits>
13 
14 namespace Acts {
15 
16 template <typename external_spacepoint_t>
17 Seedfinder<external_spacepoint_t, Acts::Cuda>::Seedfinder(
18  Acts::SeedfinderConfig<external_spacepoint_t> config)
19  : m_config(std::move(config)) {
20  // calculation of scattering using the highland formula
21  // convert pT to p once theta angle is known
22  m_config.highland = 13.6 * std::sqrt(m_config.radLengthPerSeed) *
23  (1 + 0.038 * std::log(m_config.radLengthPerSeed));
24  float maxScatteringAngle = m_config.highland / m_config.minPt;
25  m_config.maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle;
26 
27  // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV and
28  // millimeter
29  // TODO: change using ACTS units
30  m_config.pTPerHelixRadius = 300. * m_config.bFieldInZ;
31  m_config.minHelixDiameter2 =
32  std::pow(m_config.minPt * 2 / m_config.pTPerHelixRadius, 2);
33  m_config.pT2perRadius =
34  std::pow(m_config.highland / m_config.pTPerHelixRadius, 2);
35 }
36 
37 // CUDA seed finding
38 template <typename external_spacepoint_t>
39 template <typename sp_range_t>
40 std::vector<Seed<external_spacepoint_t>>
41 Seedfinder<external_spacepoint_t, Acts::Cuda>::createSeedsForGroup(
42  sp_range_t bottomSPs, sp_range_t middleSPs, sp_range_t topSPs) const {
43  std::vector<Seed<external_spacepoint_t>> outputVec;
44 
45  // Get SeedfinderConfig values
46  CudaScalar<float> deltaRMin_cuda(&m_config.deltaRMin);
47  CudaScalar<float> deltaRMax_cuda(&m_config.deltaRMax);
48  CudaScalar<float> cotThetaMax_cuda(&m_config.cotThetaMax);
49  CudaScalar<float> collisionRegionMin_cuda(&m_config.collisionRegionMin);
50  CudaScalar<float> collisionRegionMax_cuda(&m_config.collisionRegionMax);
51  CudaScalar<float> maxScatteringAngle2_cuda(&m_config.maxScatteringAngle2);
52  CudaScalar<float> sigmaScattering_cuda(&m_config.sigmaScattering);
53  CudaScalar<float> minHelixDiameter2_cuda(&m_config.minHelixDiameter2);
54  CudaScalar<float> pT2perRadius_cuda(&m_config.pT2perRadius);
55  CudaScalar<float> impactMax_cuda(&m_config.impactMax);
56  const auto seedFilterConfig = m_config.seedFilter->getSeedFilterConfig();
57  CudaScalar<float> deltaInvHelixDiameter_cuda(
58  &seedFilterConfig.deltaInvHelixDiameter);
59  CudaScalar<float> impactWeightFactor_cuda(
60  &seedFilterConfig.impactWeightFactor);
61  CudaScalar<float> filterDeltaRMin_cuda(&seedFilterConfig.deltaRMin);
62  CudaScalar<float> compatSeedWeight_cuda(&seedFilterConfig.compatSeedWeight);
63  CudaScalar<size_t> compatSeedLimit_cuda(&seedFilterConfig.compatSeedLimit);
64  CpuScalar<size_t> compatSeedLimit_cpu(&compatSeedLimit_cuda);
65  //---------------------------------
66  // Algorithm 0. Matrix Flattening
67  //---------------------------------
68 
69  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*>
70  middleSPvec;
71  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*>
72  bottomSPvec;
73  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*> topSPvec;
74 
75  // Get the size of spacepoints
76  int nSpM(0);
77  int nSpB(0);
78  int nSpT(0);
79 
80  for (auto sp : middleSPs) {
81  nSpM++;
82  middleSPvec.push_back(sp);
83  }
84  for (auto sp : bottomSPs) {
85  nSpB++;
86  bottomSPvec.push_back(sp);
87  }
88  for (auto sp : topSPs) {
89  nSpT++;
90  topSPvec.push_back(sp);
91  }
92 
93  CudaScalar<int> nSpM_cuda(&nSpM);
94  CudaScalar<int> nSpB_cuda(&nSpB);
95  CudaScalar<int> nSpT_cuda(&nSpT);
96 
97  if (nSpM == 0 || nSpB == 0 || nSpT == 0)
98  return outputVec;
99 
100  // Matrix flattening
101  CpuMatrix<float> spMmat_cpu(nSpM, 6); // x y z r varR varZ
102  CpuMatrix<float> spBmat_cpu(nSpB, 6);
103  CpuMatrix<float> spTmat_cpu(nSpT, 6);
104 
105  auto fillMatrix = [](auto& mat, auto& id, auto sp) {
106  mat.set(id, 0, sp->x());
107  mat.set(id, 1, sp->y());
108  mat.set(id, 2, sp->z());
109  mat.set(id, 3, sp->radius());
110  mat.set(id, 4, sp->varianceR());
111  mat.set(id, 5, sp->varianceZ());
112  id++;
113  };
114 
115  int mIdx(0);
116  for (auto sp : middleSPs) {
117  fillMatrix(spMmat_cpu, mIdx, sp);
118  }
119  int bIdx(0);
120  for (auto sp : bottomSPs) {
121  fillMatrix(spBmat_cpu, bIdx, sp);
122  }
123  int tIdx(0);
124  for (auto sp : topSPs) {
125  fillMatrix(spTmat_cpu, tIdx, sp);
126  }
127 
128  CudaMatrix<float> spMmat_cuda(nSpM, 6, &spMmat_cpu);
129  CudaMatrix<float> spBmat_cuda(nSpB, 6, &spBmat_cpu);
130  CudaMatrix<float> spTmat_cuda(nSpT, 6, &spTmat_cpu);
131  //------------------------------------
132  // Algorithm 1. Doublet Search (DS)
133  //------------------------------------
134 
135  CudaScalar<int> nSpMcomp_cuda(new int(0));
136  CudaScalar<int> nSpBcompPerSpMMax_cuda(new int(0));
137  CudaScalar<int> nSpTcompPerSpMMax_cuda(new int(0));
138  CudaVector<int> nSpBcompPerSpM_cuda(nSpM);
139  nSpBcompPerSpM_cuda.zeros();
140  CudaVector<int> nSpTcompPerSpM_cuda(nSpM);
141  nSpTcompPerSpM_cuda.zeros();
142  CudaVector<int> McompIndex_cuda(nSpM);
143  CudaMatrix<int> BcompIndex_cuda(nSpB, nSpM);
144  CudaMatrix<int> TcompIndex_cuda(nSpT, nSpM);
145  CudaMatrix<int> tmpBcompIndex_cuda(nSpB, nSpM);
146  CudaMatrix<int> tmpTcompIndex_cuda(nSpT, nSpM);
147 
148  dim3 DS_BlockSize = m_config.maxBlockSize;
149  dim3 DS_GridSize(nSpM, 1, 1);
150 
151  searchDoublet(DS_GridSize, DS_BlockSize, nSpM_cuda.get(), spMmat_cuda.get(),
152  nSpB_cuda.get(), spBmat_cuda.get(), nSpT_cuda.get(),
153  spTmat_cuda.get(), deltaRMin_cuda.get(), deltaRMax_cuda.get(),
154  cotThetaMax_cuda.get(), collisionRegionMin_cuda.get(),
155  collisionRegionMax_cuda.get(), nSpMcomp_cuda.get(),
156  nSpBcompPerSpMMax_cuda.get(), nSpTcompPerSpMMax_cuda.get(),
157  nSpBcompPerSpM_cuda.get(), nSpTcompPerSpM_cuda.get(),
158  McompIndex_cuda.get(), BcompIndex_cuda.get(),
159  tmpBcompIndex_cuda.get(), TcompIndex_cuda.get(),
160  tmpTcompIndex_cuda.get());
161 
162  CpuScalar<int> nSpMcomp_cpu(&nSpMcomp_cuda);
163  CpuScalar<int> nSpBcompPerSpMMax_cpu(&nSpBcompPerSpMMax_cuda);
164  CpuScalar<int> nSpTcompPerSpMMax_cpu(&nSpTcompPerSpMMax_cuda);
165  CpuVector<int> nSpBcompPerSpM_cpu(nSpM, &nSpBcompPerSpM_cuda);
166  CpuVector<int> nSpTcompPerSpM_cpu(nSpM, &nSpTcompPerSpM_cuda);
167  CpuVector<int> McompIndex_cpu(nSpM, &McompIndex_cuda);
168 
169  //--------------------------------------
170  // Algorithm 2. Transform coordinate
171  //--------------------------------------
172 
173  CudaMatrix<float> spMcompMat_cuda(*nSpMcomp_cpu.get(), 6);
174  CudaMatrix<float> spBcompMatPerSpM_cuda(*nSpBcompPerSpMMax_cpu.get(),
175  (*nSpMcomp_cpu.get()) * 6);
176  CudaMatrix<float> spTcompMatPerSpM_cuda(*nSpTcompPerSpMMax_cpu.get(),
177  (*nSpMcomp_cpu.get()) * 6);
178  CudaMatrix<float> circBcompMatPerSpM_cuda(*nSpBcompPerSpMMax_cpu.get(),
179  (*nSpMcomp_cpu.get()) * 6);
180  CudaMatrix<float> circTcompMatPerSpM_cuda(*nSpTcompPerSpMMax_cpu.get(),
181  (*nSpMcomp_cpu.get()) * 6);
182 
183  dim3 TC_GridSize(*nSpMcomp_cpu.get(), 1, 1);
184  dim3 TC_BlockSize = m_config.maxBlockSize;
185 
186  transformCoordinate(
187  TC_GridSize, TC_BlockSize, nSpM_cuda.get(), spMmat_cuda.get(),
188  McompIndex_cuda.get(), nSpB_cuda.get(), spBmat_cuda.get(),
189  nSpBcompPerSpMMax_cuda.get(), BcompIndex_cuda.get(), nSpT_cuda.get(),
190  spTmat_cuda.get(), nSpTcompPerSpMMax_cuda.get(), TcompIndex_cuda.get(),
191  spMcompMat_cuda.get(), spBcompMatPerSpM_cuda.get(),
192  circBcompMatPerSpM_cuda.get(), spTcompMatPerSpM_cuda.get(),
193  circTcompMatPerSpM_cuda.get());
194 
195  //------------------------------------------------------
196  // Algorithm 3. Triplet Search (TS) & Seed filtering
197  //------------------------------------------------------
198 
199  const int nTrplPerSpMLimit =
200  m_config.nAvgTrplPerSpBLimit * (*nSpBcompPerSpMMax_cpu.get());
201  CudaScalar<int> nTrplPerSpMLimit_cuda(&nTrplPerSpMLimit);
202 
203  CudaScalar<int> nTrplPerSpBLimit_cuda(&m_config.nTrplPerSpBLimit);
204  CpuScalar<int> nTrplPerSpBLimit_cpu(
205  &nTrplPerSpBLimit_cuda); // need to be USM
206 
207  CudaVector<int> nTrplPerSpM_cuda(*nSpMcomp_cpu.get());
208  nTrplPerSpM_cuda.zeros();
209  CudaMatrix<Triplet> TripletsPerSpM_cuda(nTrplPerSpMLimit,
210  *nSpMcomp_cpu.get());
211  CpuVector<int> nTrplPerSpM_cpu(*nSpMcomp_cpu.get(), true);
212  nTrplPerSpM_cpu.zeros();
213  CpuMatrix<Triplet> TripletsPerSpM_cpu(nTrplPerSpMLimit, *nSpMcomp_cpu.get(),
214  true);
215  cudaStream_t cuStream;
216  ACTS_CUDA_ERROR_CHECK(cudaStreamCreate(&cuStream));
217 
218  for (int i_m = 0; i_m <= *nSpMcomp_cpu.get(); i_m++) {
219  cudaStreamSynchronize(cuStream);
220 
221  // Search Triplet
222  if (i_m < *nSpMcomp_cpu.get()) {
223  int mIndex = *McompIndex_cpu.get(i_m);
224  int* nSpBcompPerSpM = nSpBcompPerSpM_cpu.get(mIndex);
225  int* nSpTcompPerSpM = nSpTcompPerSpM_cpu.get(mIndex);
226 
227  dim3 TS_GridSize(*nSpBcompPerSpM, 1, 1);
228  dim3 TS_BlockSize =
229  dim3(fmin(m_config.maxBlockSize, *nSpTcompPerSpM), 1, 1);
230 
231  searchTriplet(
232  TS_GridSize, TS_BlockSize, nSpTcompPerSpM_cpu.get(mIndex),
233  nSpTcompPerSpM_cuda.get(mIndex), nSpMcomp_cuda.get(),
234  spMcompMat_cuda.get(i_m, 0), nSpBcompPerSpMMax_cuda.get(),
235  BcompIndex_cuda.get(0, i_m), circBcompMatPerSpM_cuda.get(0, 6 * i_m),
236  nSpTcompPerSpMMax_cuda.get(), TcompIndex_cuda.get(0, i_m),
237  spTcompMatPerSpM_cuda.get(0, 6 * i_m),
238  circTcompMatPerSpM_cuda.get(0, 6 * i_m),
239  // Seed finder config
240  maxScatteringAngle2_cuda.get(), sigmaScattering_cuda.get(),
241  minHelixDiameter2_cuda.get(), pT2perRadius_cuda.get(),
242  impactMax_cuda.get(), nTrplPerSpMLimit_cuda.get(),
243  nTrplPerSpBLimit_cpu.get(), nTrplPerSpBLimit_cuda.get(),
244  deltaInvHelixDiameter_cuda.get(), impactWeightFactor_cuda.get(),
245  filterDeltaRMin_cuda.get(), compatSeedWeight_cuda.get(),
246  compatSeedLimit_cpu.get(), compatSeedLimit_cuda.get(),
247  // output
248  nTrplPerSpM_cuda.get(i_m), TripletsPerSpM_cuda.get(0, i_m),
249  &cuStream);
250  nTrplPerSpM_cpu.copyD2H(nTrplPerSpM_cuda.get(i_m), 1, i_m, &cuStream);
251 
252  TripletsPerSpM_cpu.copyD2H(TripletsPerSpM_cuda.get(0, i_m),
253  nTrplPerSpMLimit, nTrplPerSpMLimit * i_m,
254  &cuStream);
255  }
256 
257  if (i_m > 0) {
258  const auto m_experimentCuts = m_config.seedFilter->getExperimentCuts();
259  std::vector<std::pair<
260  float, std::unique_ptr<const InternalSeed<external_spacepoint_t>>>>
261  seedsPerSpM;
262 
263  for (int i = 0; i < *nTrplPerSpM_cpu.get(i_m - 1); i++) {
264  auto& triplet = *TripletsPerSpM_cpu.get(i, i_m - 1);
265  int mIndex = *McompIndex_cpu.get(i_m - 1);
266  int bIndex = triplet.bIndex;
267  int tIndex = triplet.tIndex;
268 
269  auto& bottomSP = *bottomSPvec[bIndex];
270  auto& middleSP = *middleSPvec[mIndex];
271  auto& topSP = *topSPvec[tIndex];
272  if (m_experimentCuts != nullptr) {
273  // add detector specific considerations on the seed weight
274  triplet.weight +=
275  m_experimentCuts->seedWeight(bottomSP, middleSP, topSP);
276  // discard seeds according to detector specific cuts (e.g.: weight)
277  if (!m_experimentCuts->singleSeedCut(triplet.weight, bottomSP,
278  middleSP, topSP)) {
279  continue;
280  }
281  }
282 
283  float Zob = 0; // It is not used in the seed filter but needs to be
284  // fixed anyway...
285 
286  seedsPerSpM.push_back(std::make_pair(
287  triplet.weight,
288  std::make_unique<const InternalSeed<external_spacepoint_t>>(
289  bottomSP, middleSP, topSP, Zob)));
290  }
291 
292  m_config.seedFilter->filterSeeds_1SpFixed(seedsPerSpM, outputVec);
293  }
294  }
295  ACTS_CUDA_ERROR_CHECK(cudaStreamDestroy(cuStream));
296  return outputVec;
297 }
298 } // namespace Acts