ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
ImpactPointEstimatorTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file ImpactPointEstimatorTests.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 
19 #include "Acts/Utilities/Units.hpp"
22 
23 using namespace Acts::UnitLiterals;
24 
25 namespace Acts {
26 namespace Test {
27 
29 using Propagator = Propagator<EigenStepper<ConstantBField>>;
30 // Create a test context
33 
34 // Track d0 distribution
35 std::uniform_real_distribution<> d0Dist(-0.01_mm, 0.01_mm);
36 // Track z0 distribution
37 std::uniform_real_distribution<> z0Dist(-0.2_mm, 0.2_mm);
38 // Track pT distribution
39 std::uniform_real_distribution<> pTDist(0.4_GeV, 10._GeV);
40 // Track phi distribution
41 std::uniform_real_distribution<> phiDist(-M_PI, M_PI);
42 // Track theta distribution
43 std::uniform_real_distribution<> thetaDist(1.0, M_PI - 1.0);
44 // Track IP resolution distribution
45 std::uniform_real_distribution<> resIPDist(0., 100._um);
46 // Track angular distribution
47 std::uniform_real_distribution<> resAngDist(0., 0.1);
48 // Track q/p resolution distribution
49 std::uniform_real_distribution<> resQoPDist(-0.1, 0.1);
50 // Track charge helper distribution
51 std::uniform_real_distribution<> qDist(-1, 1);
52 // Vertex x/y position distribution
53 std::uniform_real_distribution<> vXYDist(-0.1_mm, 0.1_mm);
54 // Vertex z position distribution
55 std::uniform_real_distribution<> vZDist(-20_mm, 20_mm);
56 // Number of tracks distritbution
57 std::uniform_int_distribution<> nTracksDist(3, 10);
58 
61 BOOST_AUTO_TEST_CASE(impactpoint_estimator_params_distance_test) {
62  // Debug mode
63  bool debugMode = false;
64  // Number of tests
65  unsigned int nTests = 10;
66 
67  // Set up RNG
68  int mySeed = 31415;
69  std::mt19937 gen(mySeed);
70 
71  // Set up constant B-Field
72  ConstantBField bField(Vector3D(0., 0., 2._T));
73 
74  // Set up Eigenstepper
76 
77  // Set up propagator with void navigator
78  auto propagator = std::make_shared<Propagator>(stepper);
79 
80  // Set up the ImpactPointEstimator
82  bField, propagator);
83 
85 
86  // Reference position
87  Vector3D refPosition(0., 0., 0.);
88 
89  // Start running tests
90  for (unsigned int i = 0; i < nTests; i++) {
91  // Create a track
92  // Resolutions
93  double resD0 = resIPDist(gen);
94  double resZ0 = resIPDist(gen);
95  double resPh = resAngDist(gen);
96  double resTh = resAngDist(gen);
97  double resQp = resQoPDist(gen);
98 
99  // Covariance matrix
100  Covariance covMat;
101  covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
102  0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
103  0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
104 
105  // The charge
106  double q = qDist(gen) < 0 ? -1. : 1.;
107 
108  // Impact parameters (IP)
109  double d0 = d0Dist(gen);
110  double z0 = z0Dist(gen);
111 
112  if (debugMode) {
113  std::cout << "IP: (" << d0 << "," << z0 << ")" << std::endl;
114  }
115 
116  // The track parameters
118  paramVec << d0, z0, phiDist(gen), thetaDist(gen), q / pTDist(gen), 0.;
119 
120  // Corresponding surface
121  std::shared_ptr<PerigeeSurface> perigeeSurface =
122  Surface::makeShared<PerigeeSurface>(refPosition);
123 
124  // Creating the track
125  BoundParameters myTrack(geoContext, std::move(covMat), paramVec,
126  perigeeSurface);
127 
128  // Distance in transverse plane
129  double transverseDist = std::sqrt(std::pow(d0, 2) + std::pow(z0, 2));
130 
131  // Estimate 3D distance
132  auto distanceRes =
133  ipEstimator.calculate3dDistance(geoContext, myTrack, refPosition);
134  BOOST_CHECK(distanceRes.ok());
135 
136  double distance = *distanceRes;
137 
138  BOOST_CHECK(distance < transverseDist);
139 
140  if (debugMode) {
141  std::cout << std::setprecision(10)
142  << "Distance in transverse plane: " << transverseDist
143  << std::endl;
144  std::cout << std::setprecision(10) << "Distance in 3D: " << distance
145  << std::endl;
146  }
147 
148  auto res = ipEstimator.estimate3DImpactParameters(
149  geoContext, magFieldContext, myTrack, refPosition);
150 
151  BOOST_CHECK(res.ok());
152 
153  BoundParameters trackAtIP3d = std::move(**res);
154 
155  const auto& myTrackParams = myTrack.parameters();
156  const auto& trackIP3dParams = trackAtIP3d.parameters();
157 
158  // d0 and z0 should have changed
159  BOOST_CHECK_NE(myTrackParams[ParID_t::eLOC_D0],
160  trackIP3dParams[ParID_t::eLOC_D0]);
161  BOOST_CHECK_NE(myTrackParams[ParID_t::eLOC_Z0],
162  trackIP3dParams[ParID_t::eLOC_Z0]);
163  // Theta along helix and q/p shoud remain the same
164  CHECK_CLOSE_REL(myTrackParams[ParID_t::eTHETA],
165  trackIP3dParams[ParID_t::eTHETA], 1e-5);
166  CHECK_CLOSE_REL(myTrackParams[ParID_t::eQOP],
167  trackIP3dParams[ParID_t::eQOP], 1e-5);
168 
169  if (debugMode) {
170  std::cout << std::setprecision(10) << "Old track parameters: \n"
171  << myTrackParams << std::endl;
172  std::cout << std::setprecision(10) << "Parameters at IP3d: \n"
173  << trackIP3dParams << std::endl;
174  }
175  } // end for loop tests
176 }
177 
180 BOOST_AUTO_TEST_CASE(impactpoint_estimator_compatibility_test) {
181  // Debug mode
182  bool debugMode = false;
183  // Number of tests
184  unsigned int nTests = 10;
185 
186  // Set up RNG
187  int mySeed = 31415;
188  std::mt19937 gen(mySeed);
189 
190  // Set up constant B-Field
191  ConstantBField bField(Vector3D(0., 0., 2.) * units::_T);
192 
193  // Set up Eigenstepper
195 
196  // Set up propagator with void navigator
197  auto propagator = std::make_shared<Propagator>(stepper);
198 
199  // Set up the ImpactPointEstimator
201  bField, propagator);
202 
204 
205  // Reference position
206  Vector3D refPosition(0., 0., 0.);
207 
208  // Lists to store distances and comp values
209  std::vector<double> distancesList;
210  std::vector<double> compatibilityList;
211 
212  // Generate covariance matrix values once
213  // for all cov matrices below
214  // Resolutions
215  double resD0 = resIPDist(gen);
216  double resZ0 = resIPDist(gen);
217  double resPh = resAngDist(gen);
218  double resTh = resAngDist(gen);
219  double resQp = resQoPDist(gen);
220 
221  // Start running tests
222  for (unsigned int i = 0; i < nTests; i++) {
223  // Create a track
224 
225  // Covariance matrix
226  Covariance covMat;
227  covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
228  0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
229  0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
230 
231  // The charge
232  double q = qDist(gen) < 0 ? -1. : 1.;
233 
234  // Impact parameters (IP)
235  double d0 = d0Dist(gen);
236  double z0 = z0Dist(gen);
237 
238  // The track parameters
240  paramVec << d0, z0, phiDist(gen), thetaDist(gen), q / pTDist(gen), 0.;
241 
242  // Corresponding surface
243  std::shared_ptr<PerigeeSurface> perigeeSurface =
244  Surface::makeShared<PerigeeSurface>(refPosition);
245 
246  // Creating the track
247  BoundParameters myTrack(geoContext, std::move(covMat), paramVec,
248  perigeeSurface);
249 
250  // Estimate 3D distance
251  auto distanceRes =
252  ipEstimator.calculate3dDistance(geoContext, myTrack, refPosition);
253  BOOST_CHECK(distanceRes.ok());
254 
255  distancesList.push_back(*distanceRes);
256 
257  auto res = ipEstimator.estimate3DImpactParameters(
258  geoContext, magFieldContext, myTrack, refPosition);
259 
260  BOOST_CHECK(res.ok());
261 
262  BoundParameters params = std::move(**res);
263 
264  auto compRes =
265  ipEstimator.get3dVertexCompatibility(geoContext, &params, refPosition);
266 
267  BOOST_CHECK(compRes.ok());
268 
269  compatibilityList.push_back(*compRes);
270 
271  } // end create tracks loop
272 
273  // Now test for all above constructed tracks
274  // if distances and compatibility values are
275  // compatible with one another
276  for (unsigned int i = 0; i < nTests; i++) {
277  for (unsigned int j = i + 1; j < nTests; j++) {
278  double relDiffComp =
279  (compatibilityList[i] - compatibilityList[j]) / compatibilityList[i];
280 
281  double relDiffDist =
282  (distancesList[i] - distancesList[j]) / distancesList[i];
283 
284  if (debugMode) {
285  std::cout << "Comparing track " << i << " with track " << j
286  << std::endl;
287  std::cout << "\t" << i << ": Comp.: " << compatibilityList[i]
288  << ", dist.: " << distancesList[i] << std::endl;
289  std::cout << "\t" << j << ": Comp.: " << compatibilityList[j]
290  << ", dist.: " << distancesList[j] << std::endl;
291  std::cout << "\t Rel.diff.: Comp(1-2)/1: " << relDiffComp
292  << ", Dist(1-2)/1: " << relDiffDist << std::endl;
293  }
294 
295  // Relative differences of compatibility values and distances
296  // should have the same sign, i.e. the following product
297  // should always be positive
298  double res = relDiffComp * relDiffDist;
299 
300  BOOST_CHECK(res > 0.);
301  }
302  }
303 }
304 
308 BOOST_AUTO_TEST_CASE(impactpoint_estimator_athena_test) {
309  // Set up constant B-Field
310  ConstantBField bField(Vector3D(0., 0., 1.9971546939_T));
311 
312  // Set up Eigenstepper
314 
315  // Set up propagator with void navigator
316  auto propagator = std::make_shared<Propagator>(stepper);
317 
318  // Set up the ImpactPointEstimator
320  bField, propagator);
321 
323 
324  // Use same values as in Athena unit test
325  Vector3D pos1(2_mm, 1_mm, -10_mm);
326  Vector3D mom1(400_MeV, 600_MeV, 200_MeV);
327  Vector3D vtxPos(1.2_mm, 0.8_mm, -7_mm);
328 
329  // Start creating some track parameters
330  Covariance covMat = Covariance::Identity();
331  std::shared_ptr<PerigeeSurface> perigeeSurface =
332  Surface::makeShared<PerigeeSurface>(pos1);
333 
334  // Some fixed track parameter values
335  BoundParameters params1(geoContext, covMat, pos1, mom1, 1, 0, perigeeSurface);
336 
337  auto res1 = ipEstimator.calculate3dDistance(geoContext, params1, vtxPos);
338  BOOST_CHECK(res1.ok());
339  double distance = (*res1);
340 
341  // Desired result from Athena unit test
342  const double result = 3.10391_mm;
343  CHECK_CLOSE_ABS(distance, result, 0.00001_mm);
344 
345  auto res2 = ipEstimator.estimate3DImpactParameters(
346  geoContext, magFieldContext, params1, vtxPos);
347  BOOST_CHECK(res2.ok());
348  BoundParameters endParams = std::move(**res2);
349  Vector3D surfaceCenter = endParams.referenceSurface().center(geoContext);
350 
351  BOOST_CHECK_EQUAL(surfaceCenter, vtxPos);
352 }
353 
357 BOOST_AUTO_TEST_CASE(impactpoint_estimator_parameter_estimation_test) {
358  // Number of tracks to test with
359  unsigned int nTracks = 10;
360 
361  // Set up RNG
362  int mySeed = 31415;
363  std::mt19937 gen(mySeed);
364 
365  // Set up constant B-Field
366  ConstantBField bField(0.0, 0.0, 1_T);
367 
368  // Set up Eigenstepper
370 
371  // Set up propagator with void navigator
372  auto propagator = std::make_shared<Propagator>(stepper);
373 
374  // Create perigee surface
375  std::shared_ptr<PerigeeSurface> perigeeSurface =
376  Surface::makeShared<PerigeeSurface>(Vector3D(0., 0., 0.));
377 
378  // Create position of vertex and perigee surface
379  double x = vXYDist(gen);
380  double y = vXYDist(gen);
381  double z = vZDist(gen);
382 
383  SpacePointVector vertexPosition(x, y, z, 0.);
384 
385  // Constraint for vertex fit
386  Vertex<BoundParameters> myConstraint;
387  // Some abitrary values
388  SpacePointSymMatrix myCovMat = SpacePointSymMatrix::Zero();
389  myCovMat(0, 0) = 30.;
390  myCovMat(1, 1) = 30.;
391  myCovMat(2, 2) = 30.;
392  myCovMat(3, 3) = 30.;
393  myConstraint.setFullCovariance(std::move(myCovMat));
394  myConstraint.setFullPosition(vertexPosition);
395 
396  // Calculate d0 and z0 corresponding to vertex position
397  double d0_v = std::sqrt(x * x + y * y);
398  double z0_v = z;
399 
400  // Set up the ImpactPointEstimator
402  bField, propagator);
404 
405  // Construct random track emerging from vicinity of vertex position
406  // Vector to store track objects used for vertex fit
407  for (unsigned int iTrack = 0; iTrack < nTracks; iTrack++) {
408  // Construct positive or negative charge randomly
409  double q = qDist(gen) < 0 ? -1. : 1.;
410 
411  // Construct random track parameters
412  BoundVector paramVec;
413  paramVec << d0_v + d0Dist(gen), z0_v + z0Dist(gen), phiDist(gen),
414  thetaDist(gen), q / pTDist(gen), 0.;
415 
416  // Resolutions
417  double resD0 = resIPDist(gen);
418  double resZ0 = resIPDist(gen);
419  double resPh = resAngDist(gen);
420  double resTh = resAngDist(gen);
421  double resQp = resQoPDist(gen);
422 
423  // Fill vector of track objects with simple covariance matrix
424  Covariance covMat;
425  covMat << resD0 * resD0, 0., 0., 0., 0., 0., 0., resZ0 * resZ0, 0., 0., 0.,
426  0., 0., 0., resPh * resPh, 0., 0., 0., 0., 0., 0., resTh * resTh, 0.,
427  0., 0., 0., 0., 0., resQp * resQp, 0., 0., 0., 0., 0., 0., 1.;
428 
429  BoundParameters track = BoundParameters(geoContext, std::move(covMat),
430  paramVec, perigeeSurface);
431 
432  // Check if IP are retrieved
433  ImpactParametersAndSigma output =
434  ipEstimator
435  .estimateImpactParameters(track, myConstraint, geoContext,
437  .value();
438  BOOST_CHECK_NE(output.IPd0, 0.);
439  BOOST_CHECK_NE(output.IPz0, 0.);
440  }
441 }
442 
443 } // namespace Test
444 } // namespace Acts