Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TrackDensityVertexFinderTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TrackDensityVertexFinderTests.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 
13 #include "Acts/Definitions/Algebra.hpp"
14 #include "Acts/Definitions/Common.hpp"
15 #include "Acts/Definitions/TrackParametrization.hpp"
16 #include "Acts/Definitions/Units.hpp"
17 #include "Acts/EventData/Charge.hpp"
18 #include "Acts/EventData/GenericBoundTrackParameters.hpp"
19 #include "Acts/EventData/TrackParameters.hpp"
20 #include "Acts/Geometry/GeometryContext.hpp"
21 #include "Acts/MagneticField/MagneticFieldContext.hpp"
22 #include "Acts/Surfaces/PerigeeSurface.hpp"
23 #include "Acts/Surfaces/Surface.hpp"
24 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
25 #include "Acts/Utilities/Intersection.hpp"
26 #include "Acts/Utilities/Result.hpp"
27 #include "Acts/Utilities/UnitVectors.hpp"
28 #include "Acts/Vertexing/DummyVertexFitter.hpp"
29 #include "Acts/Vertexing/GaussianTrackDensity.hpp"
30 #include "Acts/Vertexing/TrackDensityVertexFinder.hpp"
31 #include "Acts/Vertexing/Vertex.hpp"
32 #include "Acts/Vertexing/VertexingOptions.hpp"
33 
34 #include <algorithm>
35 #include <cmath>
36 #include <functional>
37 #include <iostream>
38 #include <memory>
39 #include <random>
40 #include <string>
41 #include <system_error>
42 #include <vector>
43 
44 namespace bdata = boost::unit_test::data;
45 using namespace Acts::UnitLiterals;
46 using Acts::VectorHelpers::makeVector4;
47 
48 namespace Acts {
49 namespace Test {
50 
51 using Covariance = BoundSquareMatrix;
52 
53 // Create a test context
54 GeometryContext geoContext = GeometryContext();
55 MagneticFieldContext magFieldContext = MagneticFieldContext();
56 
62 BOOST_AUTO_TEST_CASE(track_density_finder_test) {
63  // Define some track parameter properties
64  Vector3 pos0{0, 0, 0};
65  Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
66  Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
67  Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
68  Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
69  Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
70  Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
71 
72  VertexingOptions<BoundTrackParameters> vertexingOptions(geoContext,
73  magFieldContext);
74  using Finder =
75  TrackDensityVertexFinder<DummyVertexFitter<>,
76  GaussianTrackDensity<BoundTrackParameters>>;
77  Finder finder;
78  Finder::State state;
79 
80  // Start creating some track parameters
81  Covariance covMat = Covariance::Identity();
82  std::shared_ptr<PerigeeSurface> perigeeSurface =
83  Surface::makeShared<PerigeeSurface>(pos0);
84 
85  // Test finder for some fixed track parameter values
86  auto params1a =
87  BoundTrackParameters::create(
88  perigeeSurface, geoContext, makeVector4(pos1a, 0), mom1a.normalized(),
89  1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
90  .value();
91  auto params1b =
92  BoundTrackParameters::create(
93  perigeeSurface, geoContext, makeVector4(pos1b, 0), mom1b.normalized(),
94  -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
95  .value();
96  auto params1c =
97  BoundTrackParameters::create(
98  perigeeSurface, geoContext, makeVector4(pos1c, 0), mom1c.normalized(),
99  1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
100  .value();
101 
102  // Vectors of track parameters in different orders
103  std::vector<const BoundTrackParameters*> vec1 = {&params1a, &params1b,
104  &params1c};
105  std::vector<const BoundTrackParameters*> vec2 = {&params1c, &params1a,
106  &params1b};
107 
108  auto res1 = finder.find(vec1, vertexingOptions, state);
109  auto res2 = finder.find(vec2, vertexingOptions, state);
110 
111  if (!res1.ok()) {
112  std::cout << res1.error().message() << std::endl;
113  }
114 
115  if (!res2.ok()) {
116  std::cout << res2.error().message() << std::endl;
117  }
118 
119  if (res1.ok() and res2.ok()) {
120  BOOST_CHECK(!(*res1).empty());
121  BOOST_CHECK(!(*res2).empty());
122  Vector3 result1 = (*res1).back().position();
123  Vector3 result2 = (*res2).back().position();
124  BOOST_CHECK(result1 == result2);
125  }
126 }
127 
132 BOOST_AUTO_TEST_CASE(track_density_finder_constr_test) {
133  // Define some track parameter properties
134  Vector3 pos0{0, 0, 0};
135  Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
136  Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
137  Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
138  Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
139  Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
140  Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
141 
142  // From Athena VertexSeedFinderTestAlg
143  double const expectedZResult = -13.013;
144 
145  // Create constraint for seed finding
146  Vector3 constraintPos{1.7_mm, 1.3_mm, -6_mm};
147  SquareMatrix3 constrCov = ActsSquareMatrix<3>::Identity();
148 
149  Vertex<BoundTrackParameters> constraint(constraintPos);
150  constraint.setCovariance(constrCov);
151 
152  // Finder options
153  VertexingOptions<BoundTrackParameters> vertexingOptions(
154  geoContext, magFieldContext, constraint);
155  using Finder =
156  TrackDensityVertexFinder<DummyVertexFitter<>,
157  GaussianTrackDensity<BoundTrackParameters>>;
158  Finder finder;
159  Finder::State state;
160 
161  // Start creating some track parameters
162  Covariance covMat = Covariance::Identity();
163  std::shared_ptr<PerigeeSurface> perigeeSurface =
164  Surface::makeShared<PerigeeSurface>(pos0);
165 
166  // Test finder for some fixed track parameter values
167  auto params1a =
168  BoundTrackParameters::create(
169  perigeeSurface, geoContext, makeVector4(pos1a, 0), mom1a.normalized(),
170  1_e / mom1a.norm(), covMat, ParticleHypothesis::pion())
171  .value();
172  auto params1b =
173  BoundTrackParameters::create(
174  perigeeSurface, geoContext, makeVector4(pos1b, 0), mom1b.normalized(),
175  -1_e / mom1b.norm(), covMat, ParticleHypothesis::pion())
176  .value();
177  auto params1c =
178  BoundTrackParameters::create(
179  perigeeSurface, geoContext, makeVector4(pos1c, 0), mom1c.normalized(),
180  -1_e / mom1c.norm(), covMat, ParticleHypothesis::pion())
181  .value();
182 
183  // Vector of track parameters
184  std::vector<const BoundTrackParameters*> vec1 = {&params1a, &params1b,
185  &params1c};
186 
187  auto res = finder.find(vec1, vertexingOptions, state);
188 
189  if (!res.ok()) {
190  std::cout << res.error().message() << std::endl;
191  }
192 
193  if (res.ok()) {
194  BOOST_CHECK(!(*res).empty());
195  Vector3 result = (*res).back().position();
196 
197  BOOST_CHECK(result[eX] == constraintPos[eX]);
198  BOOST_CHECK(result[eY] == constraintPos[eY]);
199  CHECK_CLOSE_ABS(result[eZ], expectedZResult, 0.001_mm);
200  }
201 }
202 
203 const double zVertexPos = 12.;
204 // x position
205 std::normal_distribution<double> xdist(1_mm, 0.1_mm);
206 // y position
207 std::normal_distribution<double> ydist(-0.7_mm, 0.1_mm);
208 // z1 position
209 std::normal_distribution<double> z1dist(zVertexPos * 1_mm, 1_mm);
210 // z2 position
211 std::normal_distribution<double> z2dist(-3_mm, 0.5_mm);
212 // Track pT distribution
213 std::uniform_real_distribution<double> pTDist(0.1_GeV, 100_GeV);
214 // Track phi distribution
215 std::uniform_real_distribution<double> phiDist(-M_PI, M_PI);
216 // Track eta distribution
217 std::uniform_real_distribution<double> etaDist(-4., 4.);
218 
223 BOOST_AUTO_TEST_CASE(track_density_finder_random_test) {
224  Covariance covMat = Covariance::Identity();
225 
226  // Perigee surface for track parameters
227  Vector3 pos0{0, 0, 0};
228  std::shared_ptr<PerigeeSurface> perigeeSurface =
229  Surface::makeShared<PerigeeSurface>(pos0);
230 
231  VertexingOptions<BoundTrackParameters> vertexingOptions(geoContext,
232  magFieldContext);
233  using Finder =
234  TrackDensityVertexFinder<DummyVertexFitter<>,
235  GaussianTrackDensity<BoundTrackParameters>>;
236  Finder finder;
237  Finder::State state;
238 
239  int mySeed = 31415;
240  std::mt19937 gen(mySeed);
241  unsigned int nTracks = 200;
242 
243  std::vector<BoundTrackParameters> trackVec;
244  trackVec.reserve(nTracks);
245 
246  // Create nTracks tracks for test case
247  for (unsigned int i = 0; i < nTracks; i++) {
248  // The position of the particle
249  Vector3 pos(xdist(gen), ydist(gen), 0);
250 
251  // Create momentum and charge of track
252  double pt = pTDist(gen);
253  double phi = phiDist(gen);
254  double eta = etaDist(gen);
255  double charge = etaDist(gen) > 0 ? 1 : -1;
256 
257  // project the position on the surface
258  Vector3 direction = makeDirectionFromPhiEta(phi, eta);
259  auto intersection =
260  perigeeSurface->intersect(geoContext, pos, direction).closest();
261  pos = intersection.position();
262 
263  // Produce most of the tracks at near z1 position,
264  // some near z2. Highest track density then expected at z1
265  pos[eZ] = ((i % 4) == 0) ? z2dist(gen) : z1dist(gen);
266 
267  trackVec.push_back(BoundTrackParameters::create(
268  perigeeSurface, geoContext, makeVector4(pos, 0),
269  direction, charge / pt, covMat,
270  ParticleHypothesis::pion())
271  .value());
272  }
273 
274  std::vector<const BoundTrackParameters*> trackPtrVec;
275  for (const auto& trk : trackVec) {
276  trackPtrVec.push_back(&trk);
277  }
278 
279  auto res3 = finder.find(trackPtrVec, vertexingOptions, state);
280  if (!res3.ok()) {
281  std::cout << res3.error().message() << std::endl;
282  }
283 
284  if (res3.ok()) {
285  BOOST_CHECK(!(*res3).empty());
286  Vector3 result = (*res3).back().position();
287  CHECK_CLOSE_ABS(result[eZ], zVertexPos, 1_mm);
288  }
289 }
290 
291 // Dummy user-defined InputTrack type
292 struct InputTrack {
293  InputTrack(const BoundTrackParameters& params) : m_parameters(params) {}
294 
295  const BoundTrackParameters& parameters() const { return m_parameters; }
296 
297  // store e.g. link to original objects here
298 
299  private:
300  BoundTrackParameters m_parameters;
301 };
302 
307 BOOST_AUTO_TEST_CASE(track_density_finder_usertrack_test) {
308  // Define some track parameter properties
309  Vector3 pos0{0, 0, 0};
310  Vector3 pos1a{1.86052_mm, -1.24035_mm, -10_mm};
311  Vector3 mom1a{400_MeV, 600_MeV, 200_MeV};
312  Vector3 pos1b{-1.24035_mm, 1.86052_mm, -3_mm};
313  Vector3 mom1b{600_MeV, 400_MeV, -200_MeV};
314  Vector3 pos1c{1.69457_mm, -0.50837_mm, -7_mm};
315  Vector3 mom1c{300_MeV, 1000_MeV, 100_MeV};
316 
317  // From Athena VertexSeedFinderTestAlg
318  double const expectedZResult = -13.013;
319 
320  // Create constraint for seed finding
321  Vector3 constraintPos{1.7_mm, 1.3_mm, -6_mm};
322  SquareMatrix3 constrCov = SquareMatrix3::Identity();
323 
324  Vertex<InputTrack> constraint(constraintPos);
325  constraint.setCovariance(constrCov);
326 
327  // Finder options
328  VertexingOptions<InputTrack> vertexingOptions(geoContext, magFieldContext,
329  constraint);
330 
331  std::function<BoundTrackParameters(InputTrack)> extractParameters =
332  [](const InputTrack& params) { return params.parameters(); };
333 
334  using Finder = TrackDensityVertexFinder<DummyVertexFitter<InputTrack>,
335  GaussianTrackDensity<InputTrack>>;
336 
337  Finder finder(extractParameters);
338  Finder::State state;
339 
340  // Start creating some track parameters
341  Covariance covMat = Covariance::Identity();
342  std::shared_ptr<PerigeeSurface> perigeeSurface =
343  Surface::makeShared<PerigeeSurface>(pos0);
344 
345  // Test finder for some fixed track parameter values
346  InputTrack params1a(BoundTrackParameters::create(perigeeSurface, geoContext,
347  makeVector4(pos1a, 0), mom1a,
348  1_e / mom1a.norm(), covMat,
349  ParticleHypothesis::pion())
350  .value());
351  InputTrack params1b(BoundTrackParameters::create(perigeeSurface, geoContext,
352  makeVector4(pos1b, 0), mom1b,
353  -1_e / mom1b.norm(), covMat,
354  ParticleHypothesis::pion())
355  .value());
356  InputTrack params1c(BoundTrackParameters::create(perigeeSurface, geoContext,
357  makeVector4(pos1c, 0), mom1c,
358  -1_e / mom1c.norm(), covMat,
359  ParticleHypothesis::pion())
360  .value());
361 
362  // Vector of track parameters
363  std::vector<const InputTrack*> vec1 = {&params1a, &params1b, &params1c};
364 
365  auto res = finder.find(vec1, vertexingOptions, state);
366 
367  if (!res.ok()) {
368  std::cout << res.error().message() << std::endl;
369  }
370 
371  if (res.ok()) {
372  BOOST_CHECK(!(*res).empty());
373  Vector3 result = (*res).back().position();
374 
375  BOOST_CHECK(result[eX] == constraintPos[eX]);
376  BOOST_CHECK(result[eY] == constraintPos[eY]);
377  CHECK_CLOSE_ABS(result[eZ], expectedZResult, 0.001_mm);
378  }
379 }
380 
381 } // namespace Test
382 } // namespace Acts