ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
IterativeVertexFinderTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file IterativeVertexFinderTests.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2019 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/data/test_case.hpp>
10 #include <boost/test/tools/output_test_stream.hpp>
11 #include <boost/test/unit_test.hpp>
12 
20 #include "Acts/Utilities/Units.hpp"
28 
29 namespace bdata = boost::unit_test::data;
30 using namespace Acts::UnitLiterals;
31 
32 namespace Acts {
33 namespace Test {
34 
36 using Propagator = Propagator<EigenStepper<ConstantBField>>;
37 using Linearizer = HelicalTrackLinearizer<Propagator>;
38 
39 // Create a test context
42 
43 // Vertex x/y position distribution
44 std::uniform_real_distribution<> vXYDist(-0.1_mm, 0.1_mm);
45 // Vertex z position distribution
46 std::uniform_real_distribution<> vZDist(-20_mm, 20_mm);
47 // Track d0 distribution
48 std::uniform_real_distribution<> d0Dist(-0.01_mm, 0.01_mm);
49 // Track z0 distribution
50 std::uniform_real_distribution<> z0Dist(-0.2_mm, 0.2_mm);
51 // Track pT distribution
52 std::uniform_real_distribution<> pTDist(0.4_GeV, 10_GeV);
53 // Track phi distribution
54 std::uniform_real_distribution<> phiDist(-M_PI, M_PI);
55 // Track theta distribution
56 std::uniform_real_distribution<> thetaDist(1.0, M_PI - 1.0);
57 // Track charge helper distribution
58 std::uniform_real_distribution<> qDist(-1, 1);
59 // Track IP resolution distribution
60 std::uniform_real_distribution<> resIPDist(0., 100_um);
61 // Track angular distribution
62 std::uniform_real_distribution<> resAngDist(0., 0.1);
63 // Track q/p resolution distribution
64 std::uniform_real_distribution<> resQoPDist(-0.01, 0.01);
65 // Number of vertices per test event distribution
66 std::uniform_int_distribution<> nVertexDist(1, 6);
67 // Number of tracks per vertex distribution
68 std::uniform_int_distribution<> nTracksDist(5, 15);
69 
70 // Dummy user-defined InputTrack type
71 struct InputTrack {
72  InputTrack(const BoundParameters& params) : m_parameters(params) {}
73 
74  const BoundParameters& parameters() const { return m_parameters; }
75 
76  // store e.g. link to original objects here
77 
78  private:
79  BoundParameters m_parameters;
80 };
81 
85 BOOST_AUTO_TEST_CASE(iterative_finder_test) {
86  bool debug = false;
87 
88  // Set up RNG
89  int mySeed = 31415;
90  std::mt19937 gen(mySeed);
91 
92  // Number of test events
93  unsigned int nEvents = 5; // = nTest
94 
95  for (unsigned int iEvent = 0; iEvent < nEvents; ++iEvent) {
96  // Set up constant B-Field
97  ConstantBField bField(0.0, 0.0, 1_T);
98 
99  // Set up Eigenstepper
101 
102  // Set up propagator with void navigator
103  auto propagator = std::make_shared<Propagator>(stepper);
104 
105  // Linearizer for BoundParameters type test
106  Linearizer::Config ltConfig(bField, propagator);
107  Linearizer linearizer(ltConfig);
108 
110 
111  // Set up Billoir Vertex Fitter
112  BilloirFitter::Config vertexFitterCfg;
113 
114  BilloirFitter bFitter(vertexFitterCfg);
115 
116  // Impact point estimator
118 
119  IPEstimator::Config ipEstimatorCfg(bField, propagator);
120  IPEstimator ipEstimator(ipEstimatorCfg);
121 
122  using ZScanSeedFinder = ZScanVertexFinder<BilloirFitter>;
123 
124  static_assert(VertexFinderConcept<ZScanSeedFinder>,
125  "Vertex finder does not fulfill vertex finder concept.");
126 
127  ZScanSeedFinder::Config seedFinderCfg(ipEstimator);
128 
129  ZScanSeedFinder sFinder(seedFinderCfg);
130 
131  // Vertex Finder
133 
134  static_assert(VertexFinderConcept<VertexFinder>,
135  "Vertex finder does not fulfill vertex finder concept.");
136 
137  VertexFinder::Config cfg(bFitter, linearizer, std::move(sFinder),
138  ipEstimator);
139 
140  cfg.reassignTracksAfterFirstFit = true;
141 
142  VertexFinder finder(cfg);
143 
144  // Vector to be filled with all tracks in current event
145  std::vector<std::unique_ptr<const BoundParameters>> tracks;
146 
147  std::vector<const BoundParameters*> tracksPtr;
148 
149  // Vector to be filled with truth vertices for later comparison
150  std::vector<Vertex<BoundParameters>> trueVertices;
151 
152  // start creating event with nVertices vertices
153  unsigned int nVertices = nVertexDist(gen);
154  for (unsigned int iVertex = 0; iVertex < nVertices; ++iVertex) {
155  // Number of tracks
156  unsigned int nTracks = nTracksDist(gen);
157 
158  if (debug) {
159  std::cout << "Event " << iEvent << ", Vertex " << iVertex << "/"
160  << nVertices << " with " << nTracks << " tracks."
161  << std::endl;
162  }
163  // Create perigee surface
164  std::shared_ptr<PerigeeSurface> perigeeSurface =
165  Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));
166 
167  // Create position of vertex and perigee surface
168  double x = vXYDist(gen);
169  double y = vXYDist(gen);
170  double z = vZDist(gen);
171 
172  // True vertex
173  Vertex<BoundParameters> trueV(Vector3D(x, y, z));
174  std::vector<TrackAtVertex<BoundParameters>> tracksAtTrueVtx;
175 
176  // Calculate d0 and z0 corresponding to vertex position
177  double d0_v = sqrt(x * x + y * y);
178  double z0_v = z;
179 
180  // Construct random track emerging from vicinity of vertex position
181  // Vector to store track objects used for vertex fit
182  for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
183  // Construct positive or negative charge randomly
184  double q = qDist(gen) < 0 ? -1. : 1.;
185 
186  // Construct random track parameters
187  BoundVector paramVec;
188  double z0track = z0_v + z0Dist(gen);
189  paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
190  q / pTDist(gen), 0.;
191 
192  // Resolutions
193  double res_d0 = resIPDist(gen);
194  double res_z0 = resIPDist(gen);
195  double res_ph = resAngDist(gen);
196  double res_th = resAngDist(gen);
197  double res_qp = resQoPDist(gen);
198 
199  // Fill vector of track objects with simple covariance matrix
200  Covariance covMat;
201  covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
202  0., 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
203  res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0.,
204  0., 0., 0., 0., 1.;
205  auto params = BoundParameters(geoContext, std::move(covMat), paramVec,
206  perigeeSurface);
207 
208  tracks.push_back(std::make_unique<BoundParameters>(params));
209 
210  TrackAtVertex<BoundParameters> trAtVt(0., params, tracks.back().get());
211  tracksAtTrueVtx.push_back(trAtVt);
212  }
213 
214  trueV.setTracksAtVertex(tracksAtTrueVtx);
215  trueVertices.push_back(trueV);
216 
217  } // end loop over vertices
218 
219  // shuffle list of tracks
220  std::shuffle(std::begin(tracks), std::end(tracks), gen);
221 
222  for (const auto& trk : tracks) {
223  tracksPtr.push_back(trk.get());
224  }
225 
228 
229  // find vertices
230  auto res = finder.find(tracksPtr, vertexingOptions);
231 
232  BOOST_CHECK(res.ok());
233 
234  if (!res.ok()) {
235  std::cout << res.error().message() << std::endl;
236  }
237 
238  // Retrieve vertices found by vertex finder
239  auto vertexCollection = *res;
240 
241  // check if same amount of vertices has been found with tolerance of 2
242  CHECK_CLOSE_ABS(vertexCollection.size(), nVertices, 2);
243 
244  if (debug) {
245  std::cout << "########## RESULT: ########## Event " << iEvent
246  << std::endl;
247  std::cout << "Number of true vertices: " << nVertices << std::endl;
248  std::cout << "Number of reco vertices: " << vertexCollection.size()
249  << std::endl;
250 
251  int count = 1;
252  std::cout << "----- True vertices -----" << std::endl;
253  for (const auto& vertex : trueVertices) {
254  Vector3D pos = vertex.position();
255  std::cout << count << ". True Vertex:\t Position:"
256  << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
257  << std::endl;
258  std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
259  << std::endl;
260  count++;
261  }
262  std::cout << "----- Reco vertices -----" << std::endl;
263  count = 1;
264  for (const auto& vertex : vertexCollection) {
265  Vector3D pos = vertex.position();
266  std::cout << count << ". Reco Vertex:\t Position:"
267  << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
268  << std::endl;
269  std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
270  << std::endl;
271  count++;
272  }
273  }
274 
275  // Check if all vertices have been found with close z-values
276  bool allVerticesFound = true;
277  for (const auto& trueVertex : trueVertices) {
278  SpacePointVector truePos = trueVertex.fullPosition();
279  bool currentVertexFound = false;
280  for (const auto& recoVertex : vertexCollection) {
281  SpacePointVector recoPos = recoVertex.fullPosition();
282  // check only for close z distance
283  double zDistance = std::abs(truePos[eZ] - recoPos[eZ]);
284  if (zDistance < 2_mm) {
285  currentVertexFound = true;
286  }
287  }
288  if (!currentVertexFound) {
289  allVerticesFound = false;
290  }
291  }
292 
293  // check if found vertices have compatible z values
294  BOOST_TEST(allVerticesFound);
295  }
296 }
297 
302 BOOST_AUTO_TEST_CASE(iterative_finder_test_user_track_type) {
303  bool debug = false;
304 
305  // Set up RNG
306  int mySeed = 31415;
307  std::mt19937 gen(mySeed);
308 
309  // Number of test events
310  unsigned int nEvents = 5; // = nTest
311 
312  for (unsigned int iEvent = 0; iEvent < nEvents; ++iEvent) {
313  // Set up constant B-Field
314  ConstantBField bField(0.0, 0.0, 1_T);
315 
316  // Set up Eigenstepper
318 
319  // Set up propagator with void navigator
320  auto propagator = std::make_shared<Propagator>(stepper);
321 
322  // Linearizer for user defined InputTrack type test
323  Linearizer::Config ltConfigUT(bField, propagator);
324  Linearizer linearizer(ltConfigUT);
325 
326  // Set up vertex fitter for user track type
328 
329  // Create a custom std::function to extract BoundParameters from
330  // user-defined InputTrack
331  std::function<BoundParameters(InputTrack)> extractParameters =
332  [](InputTrack params) { return params.parameters(); };
333 
334  // Set up Billoir Vertex Fitter
335  BilloirFitter::Config vertexFitterCfg;
336 
337  BilloirFitter bFitter(vertexFitterCfg, extractParameters);
338 
339  // IP Estimator
341 
342  IPEstimator::Config ipEstimatorCfg(bField, propagator);
343  IPEstimator ipEstimator(ipEstimatorCfg);
344 
345  using ZScanSeedFinder = ZScanVertexFinder<BilloirFitter>;
346  ZScanSeedFinder::Config seedFinderCfg(ipEstimator);
347 
348  ZScanSeedFinder sFinder(seedFinderCfg, extractParameters);
349 
350  // Vertex Finder
352  VertexFinder::Config cfg(bFitter, linearizer, std::move(sFinder),
353  ipEstimator);
354  cfg.reassignTracksAfterFirstFit = true;
355 
356  VertexFinder finder(cfg, extractParameters);
357 
358  // Same for user track type tracks
359  std::vector<std::unique_ptr<const InputTrack>> tracks;
360 
361  std::vector<const InputTrack*> tracksPtr;
362 
363  // Vector to be filled with truth vertices for later comparison
364  std::vector<Vertex<InputTrack>> trueVertices;
365 
366  // start creating event with nVertices vertices
367  unsigned int nVertices = nVertexDist(gen);
368  for (unsigned int iVertex = 0; iVertex < nVertices; ++iVertex) {
369  // Number of tracks
370  unsigned int nTracks = nTracksDist(gen);
371 
372  if (debug) {
373  std::cout << "Event " << iEvent << ", Vertex " << iVertex << "/"
374  << nVertices << " with " << nTracks << " tracks."
375  << std::endl;
376  }
377  // Create perigee surface
378  std::shared_ptr<PerigeeSurface> perigeeSurface =
379  Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));
380 
381  // Create position of vertex and perigee surface
382  double x = vXYDist(gen);
383  double y = vXYDist(gen);
384  double z = vZDist(gen);
385 
386  // True vertex
387  Vertex<InputTrack> trueV(Vector3D(x, y, z));
388  std::vector<TrackAtVertex<InputTrack>> tracksAtTrueVtx;
389 
390  // Calculate d0 and z0 corresponding to vertex position
391  double d0_v = sqrt(x * x + y * y);
392  double z0_v = z;
393 
394  // Construct random track emerging from vicinity of vertex position
395  // Vector to store track objects used for vertex fit
396  for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
397  // Construct positive or negative charge randomly
398  double q = qDist(gen) < 0 ? -1. : 1.;
399 
400  // Construct random track parameters
401  BoundVector paramVec;
402  double z0track = z0_v + z0Dist(gen);
403  paramVec << d0_v + d0Dist(gen), z0track, phiDist(gen), thetaDist(gen),
404  q / pTDist(gen), 0.;
405 
406  // Resolutions
407  double res_d0 = resIPDist(gen);
408  double res_z0 = resIPDist(gen);
409  double res_ph = resAngDist(gen);
410  double res_th = resAngDist(gen);
411  double res_qp = resQoPDist(gen);
412 
413  // Fill vector of track objects now for user track type
414  Covariance covMat;
415 
416  covMat << res_d0 * res_d0, 0., 0., 0., 0., 0., 0., res_z0 * res_z0, 0.,
417  0., 0., 0., 0., 0., res_ph * res_ph, 0., 0., 0., 0., 0., 0.,
418  res_th * res_th, 0., 0., 0., 0., 0., 0., res_qp * res_qp, 0., 0.,
419  0., 0., 0., 0., 1.;
420  auto paramsUT = InputTrack(BoundParameters(
421  geoContext, std::move(covMat), paramVec, perigeeSurface));
422 
423  tracks.push_back(std::make_unique<InputTrack>(paramsUT));
424 
425  auto params = extractParameters(paramsUT);
426 
427  TrackAtVertex<InputTrack> trAtVt(0., params, tracks.back().get());
428  tracksAtTrueVtx.push_back(trAtVt);
429  }
430 
431  trueV.setTracksAtVertex(tracksAtTrueVtx);
432  trueVertices.push_back(trueV);
433 
434  } // end loop over vertices
435 
436  // shuffle list of tracks
437  std::shuffle(std::begin(tracks), std::end(tracks), gen);
438 
439  for (const auto& trk : tracks) {
440  tracksPtr.push_back(trk.get());
441  }
442 
443  VertexingOptions<InputTrack> vertexingOptionsUT(geoContext,
445 
446  // find vertices
447  auto res = finder.find(tracksPtr, vertexingOptionsUT);
448 
449  BOOST_CHECK(res.ok());
450 
451  if (!res.ok()) {
452  std::cout << res.error().message() << std::endl;
453  }
454 
455  // Retrieve vertices found by vertex finder
456  auto vertexCollectionUT = *res;
457 
458  // check if same amount of vertices has been found with tolerance of 2
459  CHECK_CLOSE_ABS(vertexCollectionUT.size(), nVertices, 2);
460 
461  if (debug) {
462  std::cout << "########## RESULT: ########## Event " << iEvent
463  << std::endl;
464  std::cout << "Number of true vertices: " << nVertices << std::endl;
465  std::cout << "Number of reco vertices: " << vertexCollectionUT.size()
466  << std::endl;
467 
468  int count = 1;
469  std::cout << "----- True vertices -----" << std::endl;
470  for (const auto& vertex : trueVertices) {
471  Vector3D pos = vertex.position();
472  std::cout << count << ". True Vertex:\t Position:"
473  << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
474  << std::endl;
475  std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
476  << std::endl;
477  count++;
478  }
479  std::cout << "----- Reco vertices -----" << std::endl;
480  count = 1;
481  for (const auto& vertex : vertexCollectionUT) {
482  Vector3D pos = vertex.position();
483  std::cout << count << ". Reco Vertex:\t Position:"
484  << "(" << pos[eX] << "," << pos[eY] << "," << pos[eZ] << ")"
485  << std::endl;
486  std::cout << "Number of tracks: " << vertex.tracks().size() << std::endl
487  << std::endl;
488  count++;
489  }
490  }
491 
492  // Check if all vertices have been found with close z-values
493  bool allVerticesFound = true;
494  for (const auto& trueVertex : trueVertices) {
495  SpacePointVector truePos = trueVertex.fullPosition();
496  bool currentVertexFound = false;
497  for (const auto& recoVertex : vertexCollectionUT) {
498  SpacePointVector recoPos = recoVertex.fullPosition();
499  // check only for close z distance
500  double zDistance = std::abs(truePos[eZ] - recoPos[eZ]);
501  if (zDistance < 2_mm) {
502  currentVertexFound = true;
503  }
504  }
505  if (!currentVertexFound) {
506  allVerticesFound = false;
507  }
508  }
509 
510  // check if found vertices have compatible z values
511  BOOST_TEST(allVerticesFound);
512  }
513 }
514 
515 } // namespace Test
516 } // namespace Acts