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 #pragma once
10 
11 // CUDA plugin include(s).
17 
18 // Acts include(s).
21 
22 // System include(s).
23 #include <cstring>
24 #include <vector>
25 
26 namespace Acts {
27 namespace Cuda {
28 
29 template <typename external_spacepoint_t>
32  const SeedFilterConfig& seedFilterConfig,
33  const TripletFilterConfig& tripletFilterConfig)
34  : m_commonConfig(std::move(commonConfig)),
35  m_seedFilterConfig(seedFilterConfig),
36  m_tripletFilterConfig(tripletFilterConfig) {
37  // calculation of scattering using the highland formula
38  // convert pT to p once theta angle is known
40  13.6 * std::sqrt(m_commonConfig.radLengthPerSeed) *
41  (1 + 0.038 * std::log(m_commonConfig.radLengthPerSeed));
42  float maxScatteringAngle = m_commonConfig.highland / m_commonConfig.minPt;
43  m_commonConfig.maxScatteringAngle2 = maxScatteringAngle * maxScatteringAngle;
44 
45  // helix radius in homogeneous magnetic field. Units are Kilotesla, MeV and
46  // millimeter
47  // TODO: change using ACTS units
53 }
54 
55 template <typename external_spacepoint_t>
56 template <typename sp_range_t>
57 std::vector<Seed<external_spacepoint_t>>
59  sp_range_t bottomSPs, sp_range_t middleSPs, sp_range_t topSPs) const {
60  // Create an empty vector, to be returned already early on when no seeds can
61  // be constructed.
62  std::vector<Seed<external_spacepoint_t>> outputVec;
63 
64  //---------------------------------
65  // Matrix Flattening
66  //---------------------------------
67 
68  // Create more convenient vectors out of the space point containers.
69  auto spVecMaker = [](sp_range_t spRange) {
70  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*> result;
71  for (auto* sp : spRange) {
72  result.push_back(sp);
73  }
74  return result;
75  };
76 
77  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*>
78  bottomSPVec(spVecMaker(bottomSPs));
79  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*>
80  middleSPVec(spVecMaker(middleSPs));
81  std::vector<const Acts::InternalSpacePoint<external_spacepoint_t>*> topSPVec(
82  spVecMaker(topSPs));
83 
84  // If either one of them is empty, we have nothing to find.
85  if ((middleSPVec.size() == 0) || (bottomSPVec.size() == 0) ||
86  (topSPVec.size() == 0)) {
87  return outputVec;
88  }
89 
90  // Create helper objects for storing information about the spacepoints on the
91  // host in single memory blobs.
92  auto bottomSPArray = make_host_array<Details::SpacePoint>(bottomSPVec.size());
93  auto middleSPArray = make_host_array<Details::SpacePoint>(middleSPVec.size());
94  auto topSPArray = make_host_array<Details::SpacePoint>(topSPVec.size());
95 
96  // Fill these memory blobs.
97  auto fillSPArray = [](Details::SpacePoint* array, const auto& spVec) {
98  for (std::size_t i = 0; i < spVec.size(); ++i) {
99  array[i].x = spVec[i]->x();
100  array[i].y = spVec[i]->y();
101  array[i].z = spVec[i]->z();
102  array[i].radius = spVec[i]->radius();
103  array[i].varianceR = spVec[i]->varianceR();
104  array[i].varianceZ = spVec[i]->varianceZ();
105  }
106  };
107  fillSPArray(bottomSPArray.get(), bottomSPVec);
108  fillSPArray(middleSPArray.get(), middleSPVec);
109  fillSPArray(topSPArray.get(), topSPVec);
110 
111  // Copy the memory blobs to the device.
112  auto bottomSPDeviceArray =
113  make_device_array<Details::SpacePoint>(bottomSPVec.size());
114  auto middleSPDeviceArray =
115  make_device_array<Details::SpacePoint>(middleSPVec.size());
116  auto topSPDeviceArray =
117  make_device_array<Details::SpacePoint>(topSPVec.size());
118  copyToDevice(bottomSPDeviceArray, bottomSPArray, bottomSPVec.size());
119  copyToDevice(middleSPDeviceArray, middleSPArray, middleSPVec.size());
120  copyToDevice(topSPDeviceArray, topSPArray, topSPVec.size());
121 
122  //---------------------------------
123  // GPU Execution
124  //---------------------------------
125 
126  // Matrices holding the counts of the viable bottom-middle and middle-top
127  // pairs.
128  auto dubletCountsHost = make_host_array<unsigned int>(middleSPVec.size());
129  memset(dubletCountsHost.get(), 0, middleSPVec.size() * sizeof(unsigned int));
130  auto middleBottomCounts = make_device_array<unsigned int>(middleSPVec.size());
131  copyToDevice(middleBottomCounts, dubletCountsHost, middleSPVec.size());
132  auto middleTopCounts = make_device_array<unsigned int>(middleSPVec.size());
133  copyToDevice(middleTopCounts, dubletCountsHost, middleSPVec.size());
134 
135  // Matrices holding the indices of the viable bottom-middle and middle-top
136  // pairs.
137  auto middleBottomDublets =
138  make_device_array<std::size_t>(middleSPVec.size() * bottomSPVec.size());
139  auto middleTopDublets =
140  make_device_array<std::size_t>(middleSPVec.size() * topSPVec.size());
141 
142  // Launch the dublet finding code.
144  m_commonConfig.maxBlockSize, bottomSPVec.size(), bottomSPDeviceArray,
145  middleSPVec.size(), middleSPDeviceArray, topSPVec.size(),
146  topSPDeviceArray, m_commonConfig.deltaRMin, m_commonConfig.deltaRMax,
147  m_commonConfig.cotThetaMax, m_commonConfig.collisionRegionMin,
148  m_commonConfig.collisionRegionMax, middleBottomCounts,
149  middleBottomDublets, middleTopCounts, middleTopDublets);
150 
151  // Count the number of dublets that we have to launch the subsequent steps
152  // for.
153  Details::DubletCounts dubletCounts =
154  Details::countDublets(m_commonConfig.maxBlockSize, middleSPVec.size(),
155  middleBottomCounts, middleTopCounts);
156 
157  // If no dublets/triplet candidates have been found, stop here.
158  if ((dubletCounts.nDublets == 0) || (dubletCounts.nTriplets == 0)) {
159  return outputVec;
160  }
161 
162  // Launch the triplet finding code on all of the previously found dublets.
163  auto tripletCandidates = Details::findTriplets(
164  m_commonConfig.maxBlockSize, dubletCounts, m_seedFilterConfig,
165  m_tripletFilterConfig, bottomSPVec.size(), bottomSPDeviceArray,
166  middleSPVec.size(), middleSPDeviceArray, topSPVec.size(),
167  topSPDeviceArray, middleBottomCounts, middleBottomDublets,
168  middleTopCounts, middleTopDublets, m_commonConfig.maxScatteringAngle2,
169  m_commonConfig.sigmaScattering, m_commonConfig.minHelixDiameter2,
170  m_commonConfig.pT2perRadius, m_commonConfig.impactMax);
171  assert(tripletCandidates.size() == middleSPVec.size());
172 
173  // Perform the final step of the filtering.
174  std::size_t middleIndex = 0;
175  auto triplet_itr = tripletCandidates.begin();
176  auto triplet_end = tripletCandidates.end();
177  for (; triplet_itr != triplet_end; ++triplet_itr, ++middleIndex) {
178  std::vector<std::pair<
179  float, std::unique_ptr<const InternalSeed<external_spacepoint_t>>>>
180  seedsPerSPM;
181  const auto& middleSP = *(middleSPVec[middleIndex]);
182  for (const Details::Triplet& triplet : *triplet_itr) {
183  assert(triplet.bottomIndex < bottomSPVec.size());
184  const auto& bottomSP = *(bottomSPVec[triplet.bottomIndex]);
185  assert(triplet.topIndex < topSPVec.size());
186  const auto& topSP = *(topSPVec[triplet.topIndex]);
187  seedsPerSPM.emplace_back(std::make_pair(
188  triplet.weight,
189  std::make_unique<const InternalSeed<external_spacepoint_t>>(
190  bottomSP, middleSP, topSP, 0)));
191  }
192  m_commonConfig.seedFilter->filterSeeds_1SpFixed(seedsPerSPM, outputVec);
193  }
194 
195  // Return the collected spacepoints.
196  return outputVec;
197 }
198 
199 } // namespace Cuda
200 } // namespace Acts