Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SingleSeedVertexFinderTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SingleSeedVertexFinderTests.cpp
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 <boost/test/unit_test.hpp>
10 
11 #include "Acts/Definitions/Algebra.hpp"
12 #include "Acts/Definitions/Units.hpp"
13 #include "Acts/Vertexing/SingleSeedVertexFinder.hpp"
14 
15 #include <cmath>
16 #include <iostream>
17 #include <random>
18 #include <vector>
19 
21 struct SpacePoint4SSVFT {
22  SpacePoint4SSVFT(double x, double y, double z) : m_x(x), m_y(y), m_z(z) {}
23  double m_x;
24  double m_y;
25  double m_z;
26  double x() const { return m_x; }
27  double y() const { return m_y; }
28  double z() const { return m_z; }
29  double r() const { return std::sqrt(m_x * m_x + m_y * m_y); }
30 };
31 
37 double getRndDouble(std::mt19937& gen, double from, double to) {
38  return gen() / 4294967296. * (to - from) + from;
39 }
40 
46 int getRndInt(std::mt19937& gen, int from, int to) {
47  return (int)(gen() / 4294967296. * (to - from) + from);
48 }
49 
53 std::vector<double> makePlaneFromTriplet(SpacePoint4SSVFT aa,
54  SpacePoint4SSVFT bb,
55  SpacePoint4SSVFT cc) {
56  Acts::Vector3 a{aa.x(), aa.y(), aa.z()};
57  Acts::Vector3 b{bb.x(), bb.y(), bb.z()};
58  Acts::Vector3 c{cc.x(), cc.y(), cc.z()};
59 
60  Acts::Vector3 ba = b - a, ca = c - a;
61 
62  Acts::Vector3 abg = ba.cross(ca).normalized();
63  double delta = -1. * abg.dot(a);
64 
65  // plane (alpha*x + beta*y + gamma*z + delta = 0)
66  return {abg[0], abg[1], abg[2], delta};
67 }
68 
70 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_small_planes_test) {
71  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
72  singleSeedVtxCfg.maxPhideviation = 0.2;
73  singleSeedVtxCfg.maxXYdeviation = 0.1;
74  singleSeedVtxCfg.maxXYZdeviation = 0.1;
75  singleSeedVtxCfg.minTheta = 0.5;
76  singleSeedVtxCfg.rMinNear = 6.f * Acts::UnitConstants::mm;
77  singleSeedVtxCfg.rMaxNear = 14.9f * Acts::UnitConstants::mm;
78  singleSeedVtxCfg.rMinMiddle = 15.f * Acts::UnitConstants::mm;
79  singleSeedVtxCfg.rMaxMiddle = 24.9f * Acts::UnitConstants::mm;
80  singleSeedVtxCfg.rMinFar = 25.f * Acts::UnitConstants::mm;
81  singleSeedVtxCfg.rMaxFar = 44.9f * Acts::UnitConstants::mm;
82  singleSeedVtxCfg.numPhiSlices = 10;
83  singleSeedVtxCfg.useFracPhiSlices = 1.0;
84  singleSeedVtxCfg.numZSlices = 10;
85  singleSeedVtxCfg.useFracZSlices = 1.0;
86  singleSeedVtxCfg.maxAbsZ = 50. * Acts::UnitConstants::mm;
87  singleSeedVtxCfg.maxZPosition = 20.f * Acts::UnitConstants::mm;
88  singleSeedVtxCfg.maxRPosition = 5.f * Acts::UnitConstants::mm;
89  singleSeedVtxCfg.minimalizeWRT = "planes";
90  singleSeedVtxCfg.maxIterations = 3;
91  singleSeedVtxCfg.removeFraction = 0.3;
92  singleSeedVtxCfg.minVtxShift = 0.3f * Acts::UnitConstants::mm;
93  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
94  singleSeedVtxCfg);
95 
96  double vtxX = 1.;
97  double vtxY = -1.;
98  double vtxZ = 10;
99 
100  std::vector<std::vector<double>> positions = {
101  {10.8, 0., 7.}, {20.9, 0.7, 4.},
102  {30.5, 1.9, 1.}, // plane ((-0.25,-0.25,-0.9), 9.0)
103  {2.1, 10.6, 15.2}, {2.7, 19.36666, 19.5},
104  {4.5, 34., 25.4}, // plane ((0.8, -0.3, 0.5), -6.1)
105  {-6.25, -7.9, 10.5}, {-12.8, -15.08, 11.},
106  {-20., -22., 11.5}, // plane ((-0.02,-0.05,-0.98), 9.77)
107  {-6., 8., 2.4}, {-12., 15., -3.4},
108  {-17., 21., -8.4}, // plane ((0.1,0.5,0.5), -4.6)
109  {7.8, 7.8, 16.0}, {14.8, 15., 17.},
110  {22.8, 23.3, 18.}, // this will be removed after iteration
111  {-5.1, 8.1, -10.}, {-1.3, 9., -10.4},
112  {3.1, 10.1, -11.1}}; // this will not form a triplet
113 
114  std::vector<SpacePoint4SSVFT> inputSpacepoints;
115  for (auto pos : positions) {
116  inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
117  }
118 
119  auto t1 = std::chrono::high_resolution_clock::now();
120  auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
121  auto t2 = std::chrono::high_resolution_clock::now();
122 
123  bool vtxFound = false;
124  if (vtx.ok()) {
125  std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
126  << " ms at x = " << vtx.value()[0] << "mm, y = " << vtx.value()[1]
127  << "mm, z = " << vtx.value()[2] << "mm" << std::endl;
128  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
129  << "mm, z = " << vtxZ << "mm" << std::endl;
130 
131  if (std::abs(vtxX - vtx.value()[0]) < 1e-3 &&
132  std::abs(vtxY - vtx.value()[1]) < 1e-3 &&
133  std::abs(vtxZ - vtx.value()[2]) < 1e-3) {
134  vtxFound = true;
135  }
136  } else {
137  std::cout << "Not found a vertex in the event after "
138  << (t2 - t1).count() / 1e6 << " ms" << std::endl;
139  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
140  << "mm, z = " << vtxZ << "mm" << std::endl;
141  }
142 
143  // check if found vertex has compatible position
144  BOOST_CHECK(vtxFound);
145 }
146 
148 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_small_rays_test) {
149  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
150  singleSeedVtxCfg.maxPhideviation = 0.2;
151  singleSeedVtxCfg.maxXYdeviation = 0.1;
152  singleSeedVtxCfg.maxXYZdeviation = 0.1;
153  singleSeedVtxCfg.minTheta = 0.5;
154  singleSeedVtxCfg.rMinNear = 6.f * Acts::UnitConstants::mm;
155  singleSeedVtxCfg.rMaxNear = 14.9f * Acts::UnitConstants::mm;
156  singleSeedVtxCfg.rMinMiddle = 15.f * Acts::UnitConstants::mm;
157  singleSeedVtxCfg.rMaxMiddle = 24.9f * Acts::UnitConstants::mm;
158  singleSeedVtxCfg.rMinFar = 25.f * Acts::UnitConstants::mm;
159  singleSeedVtxCfg.rMaxFar = 44.9f * Acts::UnitConstants::mm;
160  singleSeedVtxCfg.numPhiSlices = 10;
161  singleSeedVtxCfg.useFracPhiSlices = 1.0;
162  singleSeedVtxCfg.numZSlices = 10;
163  singleSeedVtxCfg.useFracZSlices = 1.0;
164  singleSeedVtxCfg.maxAbsZ = 50. * Acts::UnitConstants::mm;
165  singleSeedVtxCfg.maxZPosition = 20.f * Acts::UnitConstants::mm;
166  singleSeedVtxCfg.maxRPosition = 5.f * Acts::UnitConstants::mm;
167  singleSeedVtxCfg.minimalizeWRT = "rays";
168  singleSeedVtxCfg.maxIterations = 1;
169  singleSeedVtxCfg.removeFraction = 0.3;
170  singleSeedVtxCfg.minVtxShift = 0.3f * Acts::UnitConstants::mm;
171  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
172  singleSeedVtxCfg);
173 
174  double vtxX = 1.;
175  double vtxY = -1.;
176  double vtxZ = 10;
177 
178  std::vector<std::vector<double>> positions = {
179  {11., 0., 7.}, {21., 1., 4.},
180  {31., 2., 1.}, // this ray is Ok
181  {2., 9., 15.}, {3., 19., 20.},
182  {4.5, 34., 27.5}, // this ray is Ok
183  {-6., -8., 10.5}, {-13., -15., 11.},
184  {-20., -22., 11.5}, // this ray is Ok
185  {7., 7., 16.0}, {14., 14., 17.},
186  {21., 21., 18.}, // this will be removed after iteration
187  {-5., 8., -10.}, {-1.5, 9., -10.5},
188  {3., 10., -11.}}; // this will not form a triplet
189 
190  std::vector<SpacePoint4SSVFT> inputSpacepoints;
191  for (auto pos : positions) {
192  inputSpacepoints.emplace_back(pos[0], pos[1], pos[2]);
193  }
194 
195  auto t1 = std::chrono::high_resolution_clock::now();
196  auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
197  auto t2 = std::chrono::high_resolution_clock::now();
198 
199  bool vtxFound = false;
200  if (vtx.ok()) {
201  std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
202  << " ms at x = " << vtx.value()[0] << "mm, y = " << vtx.value()[1]
203  << "mm, z = " << vtx.value()[2] << "mm" << std::endl;
204  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
205  << "mm, z = " << vtxZ << "mm" << std::endl;
206 
207  if (std::abs(vtxX - vtx.value()[0]) < 1e-3 &&
208  std::abs(vtxY - vtx.value()[1]) < 1e-3 &&
209  std::abs(vtxZ - vtx.value()[2]) < 1e-3) {
210  vtxFound = true;
211  }
212  } else {
213  std::cout << "Not found a vertex in the event after "
214  << (t2 - t1).count() / 1e6 << " ms" << std::endl;
215  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
216  << "mm, z = " << vtxZ << "mm" << std::endl;
217  }
218 
219  // check if found vertex has compatible position
220  BOOST_CHECK(vtxFound);
221 }
222 
224 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_planes_test) {
225  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
226  singleSeedVtxCfg.maxPhideviation = 0.1;
227  singleSeedVtxCfg.maxXYdeviation = 0.1;
228  singleSeedVtxCfg.maxXYZdeviation = 0.1;
229  singleSeedVtxCfg.minTheta = 0.5;
230  singleSeedVtxCfg.rMinNear = 8.f * Acts::UnitConstants::mm;
231  singleSeedVtxCfg.rMaxNear = 12.f * Acts::UnitConstants::mm;
232  singleSeedVtxCfg.rMinMiddle = 18.f * Acts::UnitConstants::mm;
233  singleSeedVtxCfg.rMaxMiddle = 22.f * Acts::UnitConstants::mm;
234  singleSeedVtxCfg.rMinFar = 28.f * Acts::UnitConstants::mm;
235  singleSeedVtxCfg.rMaxFar = 32.f * Acts::UnitConstants::mm;
236  singleSeedVtxCfg.numPhiSlices = 60;
237  singleSeedVtxCfg.useFracPhiSlices = 0.8;
238  singleSeedVtxCfg.numZSlices = 150;
239  singleSeedVtxCfg.useFracZSlices = 0.8;
240  singleSeedVtxCfg.maxAbsZ = 75. * Acts::UnitConstants::mm;
241  singleSeedVtxCfg.maxZPosition = 15.f * Acts::UnitConstants::mm;
242  singleSeedVtxCfg.maxRPosition = 2.5f * Acts::UnitConstants::mm;
243  singleSeedVtxCfg.minimalizeWRT = "planes";
244  singleSeedVtxCfg.maxIterations = 20;
245  singleSeedVtxCfg.removeFraction = 0.1;
246  singleSeedVtxCfg.minVtxShift = 0.05f * Acts::UnitConstants::mm;
247  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
248  singleSeedVtxCfg);
249 
250  std::mt19937 gen(299792458);
251 
252  int vtxFound = 0;
253  int nEvents = 5;
254  for (int event = 0; event < nEvents; event++) {
255  double vtxX = getRndDouble(gen, -1., 1.);
256  double vtxY = getRndDouble(gen, -1., 1.);
257  double vtxZ = getRndDouble(gen, -10., 10.);
258 
259  std::vector<SpacePoint4SSVFT> inputSpacepoints;
260 
261  // make straight lines originating from the given vertex
262  int nTracks = getRndInt(gen, 200, 400);
263  for (int track = 0; track < nTracks; ++track) {
264  // initial position of the track
265  double posX = vtxX;
266  double posY = vtxY;
267  double posZ = vtxZ;
268 
269  // initial direction of the track
270  double dirX = getRndDouble(gen, -1., 1.);
271  double dirY = getRndDouble(gen, -1., 1.);
272  double dirZ = getRndDouble(gen, -1., 1.);
273  // rotation of the track
274  double theta = getRndDouble(gen, 0.03, 0.09) *
275  (getRndDouble(gen, -1., 1.) > 0 ? 1 : -1);
276  double sgn = std::copysign(1., dirY);
277 
278  for (int i = 1; i <= 3; ++i) {
279  if (i != 1) {
280  // rotate direction, not for the first time
281  double dirXtmp = cos(theta) * dirX - sin(theta) * dirY;
282  double dirYtmp = sin(theta) * dirX + cos(theta) * dirY;
283  dirX = dirXtmp;
284  dirY = dirYtmp;
285  }
286 
287  // double sgn=std::copysign(1.,dirY);
288  double dirR2 = dirX * dirX + dirY * dirY;
289  double D = posX * (posY + dirY) - posY * (posX + dirX);
290  // add some smearing to the layers
291  // layers are (9-11), (19-21), and (29-31)
292  double r = i * 10. + getRndDouble(gen, -1., 1.);
293  // intersection of the layer and the straigh line
294  double x1 =
295  (D * dirY + sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) / dirR2;
296  double y1 =
297  (-D * dirX + std::fabs(dirY) * std::sqrt(r * r * dirR2 - D * D)) /
298  dirR2;
299  // how many units from the vertex to the intersection
300  double zDist = std::fabs((x1 - posX) / dirX);
301 
302  // position of the new spacepoint
303  posX = x1;
304  posY = y1;
305  // use the same amount of units for distance in Z
306  posZ += zDist * dirZ;
307 
308  if (i == 3) {
309  // move z position, so the vertex will be part of the plane
310  auto abgd = makePlaneFromTriplet({vtxX, vtxY, vtxZ},
311  inputSpacepoints.rbegin()[1],
312  inputSpacepoints.rbegin()[0]);
313  posZ = -1. * (abgd[0] * posX + abgd[1] * posY + abgd[3]) / abgd[2];
314  }
315 
316  inputSpacepoints.emplace_back(posX, posY, posZ);
317  }
318  }
319 
320  auto t1 = std::chrono::high_resolution_clock::now();
321  auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
322  auto t2 = std::chrono::high_resolution_clock::now();
323 
324  if (vtx.ok()) {
325  std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
326  << " ms at x = " << vtx.value()[0]
327  << "mm, y = " << vtx.value()[1] << "mm, z = " << vtx.value()[2]
328  << "mm" << std::endl;
329  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
330  << "mm, z = " << vtxZ << "mm" << std::endl;
331  std::cout << "Difference is in x = " << std::abs(vtxX - vtx.value()[0])
332  << "mm, y = " << std::abs(vtxY - vtx.value()[1])
333  << "mm, z = " << std::abs(vtxZ - vtx.value()[2]) << "mm"
334  << std::endl;
335  if (std::abs(vtxX - vtx.value()[0]) < 2.0 &&
336  std::abs(vtxY - vtx.value()[1]) < 2.0 &&
337  std::abs(vtxZ - vtx.value()[2]) < 0.3) {
338  ++vtxFound;
339  }
340  } else {
341  std::cout << "Not found a vertex in the event after "
342  << (t2 - t1).count() / 1e6 << " ms" << std::endl;
343  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
344  << "mm, z = " << vtxZ << "mm" << std::endl;
345  }
346  }
347 
348  std::cout << "Found " << vtxFound << " out of " << nEvents << " vertices"
349  << std::endl;
350 
351  // check if all vertices have compatible positions
352  BOOST_CHECK(vtxFound == nEvents);
353 }
354 
356 BOOST_AUTO_TEST_CASE(single_seed_vertex_finder_full_rays_test) {
357  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT>::Config singleSeedVtxCfg;
358  singleSeedVtxCfg.maxPhideviation = 0.1;
359  singleSeedVtxCfg.maxXYdeviation = 0.1;
360  singleSeedVtxCfg.maxXYZdeviation = 0.1;
361  singleSeedVtxCfg.minTheta = 0.5;
362  singleSeedVtxCfg.rMinNear = 8.f * Acts::UnitConstants::mm;
363  singleSeedVtxCfg.rMaxNear = 12.f * Acts::UnitConstants::mm;
364  singleSeedVtxCfg.rMinMiddle = 18.f * Acts::UnitConstants::mm;
365  singleSeedVtxCfg.rMaxMiddle = 22.f * Acts::UnitConstants::mm;
366  singleSeedVtxCfg.rMinFar = 28.f * Acts::UnitConstants::mm;
367  singleSeedVtxCfg.rMaxFar = 32.f * Acts::UnitConstants::mm;
368  singleSeedVtxCfg.numPhiSlices = 60;
369  singleSeedVtxCfg.useFracPhiSlices = 0.8;
370  singleSeedVtxCfg.numZSlices = 150;
371  singleSeedVtxCfg.useFracZSlices = 0.8;
372  singleSeedVtxCfg.maxAbsZ = 75. * Acts::UnitConstants::mm;
373  singleSeedVtxCfg.maxZPosition = 15.f * Acts::UnitConstants::mm;
374  singleSeedVtxCfg.maxRPosition = 2.5f * Acts::UnitConstants::mm;
375  singleSeedVtxCfg.minimalizeWRT = "rays";
376  singleSeedVtxCfg.maxIterations = 10;
377  singleSeedVtxCfg.removeFraction = 0.1;
378  singleSeedVtxCfg.minVtxShift = 0.05f * Acts::UnitConstants::mm;
379  Acts::SingleSeedVertexFinder<SpacePoint4SSVFT> SingleSeedVertexFinder(
380  singleSeedVtxCfg);
381 
382  std::mt19937 gen(299792458);
383 
384  int vtxFound = 0;
385  int nEvents = 5;
386  for (int event = 0; event < nEvents; event++) {
387  double vtxX = getRndDouble(gen, -1., 1.);
388  double vtxY = getRndDouble(gen, -1., 1.);
389  double vtxZ = getRndDouble(gen, -10., 10.);
390 
391  std::vector<SpacePoint4SSVFT> inputSpacepoints;
392 
393  // make straight lines originating from the given vertex
394  int nTracks = getRndInt(gen, 200, 400);
395  for (int track = 0; track < nTracks; ++track) {
396  // direction of the track
397  double dirX = getRndDouble(gen, -1., 1.);
398  double dirY = getRndDouble(gen, -1., 1.);
399  double dirZ = getRndDouble(gen, -1., 1.);
400  // use upper or lower intersection?
401  int part = (getRndDouble(gen, -1., 1.) > 0 ? 1 : -1);
402 
403  for (int rIndx = 1; rIndx <= 3; rIndx += 1) {
404  double sgn = std::copysign(1., dirY);
405  double dirR2 = dirX * dirX + dirY * dirY;
406  double D = vtxX * (vtxY + dirY) - vtxY * (vtxX + dirX);
407  // add some smearing to the layers
408  // layers are (9-11), (19-21), and (29-31)
409  double r = rIndx * 10 + getRndDouble(gen, -1., 1.);
410  // intersection of the layer and the straigh line
411  double x1 =
412  (D * dirY + part * sgn * dirX * std::sqrt(r * r * dirR2 - D * D)) /
413  dirR2;
414  double y1 = (-D * dirX + part * std::fabs(dirY) *
415  std::sqrt(r * r * dirR2 - D * D)) /
416  dirR2;
417  // how many units from the vertex to the intersection
418  double zDist = std::fabs((x1 - vtxX) / dirX);
419  // use the same amount of units for distance in Z
420  inputSpacepoints.emplace_back(x1, y1, zDist * dirZ + vtxZ);
421  }
422  }
423 
424  auto t1 = std::chrono::high_resolution_clock::now();
425  auto vtx = SingleSeedVertexFinder.findVertex(inputSpacepoints);
426  auto t2 = std::chrono::high_resolution_clock::now();
427 
428  if (vtx.ok()) {
429  std::cout << "Found a vertex in the event in " << (t2 - t1).count() / 1e6
430  << " ms at x = " << vtx.value()[0]
431  << "mm, y = " << vtx.value()[1] << "mm, z = " << vtx.value()[2]
432  << "mm" << std::endl;
433  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
434  << "mm, z = " << vtxZ << "mm" << std::endl;
435  std::cout << "Difference is in x = " << std::abs(vtxX - vtx.value()[0])
436  << "mm, y = " << std::abs(vtxY - vtx.value()[1])
437  << "mm, z = " << std::abs(vtxZ - vtx.value()[2]) << "mm"
438  << std::endl;
439  if (std::abs(vtxX - vtx.value()[0]) < 0.3 &&
440  std::abs(vtxY - vtx.value()[1]) < 0.3 &&
441  std::abs(vtxZ - vtx.value()[2]) < 0.3) {
442  ++vtxFound;
443  }
444  } else {
445  std::cout << "Not found a vertex in the event after "
446  << (t2 - t1).count() / 1e6 << " ms" << std::endl;
447  std::cout << "Truth vertex was at x = " << vtxX << "mm, y = " << vtxY
448  << "mm, z = " << vtxZ << "mm" << std::endl;
449  }
450  }
451 
452  std::cout << "Found " << vtxFound << " out of " << nEvents << " vertices"
453  << std::endl;
454 
455  // check if all vertices have compatible positions
456  BOOST_CHECK(vtxFound == nEvents);
457 }