Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EstimateTrackParamsFromSeedTest.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EstimateTrackParamsFromSeedTest.cpp
1 // This file is part of the Acts project.
2 //
3 // Copyright (C) 2021 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/Direction.hpp"
13 #include "Acts/Definitions/TrackParametrization.hpp"
14 #include "Acts/Definitions/Units.hpp"
15 #include "Acts/EventData/GenericCurvilinearTrackParameters.hpp"
16 #include "Acts/EventData/TrackParameters.hpp"
17 #include "Acts/Geometry/GeometryContext.hpp"
18 #include "Acts/Geometry/GeometryIdentifier.hpp"
19 #include "Acts/Geometry/LayerCreator.hpp"
20 #include "Acts/Geometry/TrackingGeometry.hpp"
21 #include "Acts/MagneticField/ConstantBField.hpp"
22 #include "Acts/MagneticField/MagneticFieldContext.hpp"
23 #include "Acts/Propagator/EigenStepper.hpp"
24 #include "Acts/Propagator/Navigator.hpp"
25 #include "Acts/Propagator/Propagator.hpp"
26 #include "Acts/Seeding/EstimateTrackParamsFromSeed.hpp"
27 #include "Acts/Surfaces/Surface.hpp"
28 #include "Acts/Tests/CommonHelpers/CylindricalTrackingGeometry.hpp"
29 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
30 #include "Acts/Tests/CommonHelpers/MeasurementsCreator.hpp"
31 #include "Acts/Tests/CommonHelpers/TestSourceLink.hpp"
32 #include "Acts/Utilities/CalibrationContext.hpp"
33 #include "Acts/Utilities/Logger.hpp"
34 
35 #include <algorithm>
36 #include <array>
37 #include <cmath>
38 #include <iterator>
39 #include <map>
40 #include <memory>
41 #include <optional>
42 #include <ostream>
43 #include <random>
44 #include <utility>
45 #include <vector>
46 
47 #include "SpacePoint.hpp"
48 
49 namespace {
50 
51 using namespace Acts;
52 using namespace Acts::Test;
53 using namespace Acts::UnitLiterals;
54 
55 using ConstantFieldStepper = Acts::EigenStepper<>;
56 using ConstantFieldPropagator =
57  Acts::Propagator<ConstantFieldStepper, Acts::Navigator>;
58 
59 const GeometryContext geoCtx;
60 const MagneticFieldContext magCtx;
61 const CalibrationContext calCtx;
62 
63 // detector geometry
64 CylindricalTrackingGeometry geometryStore(geoCtx);
65 const auto geometry = geometryStore();
66 
67 // Two dimensional measurement with zero resolution
68 const MeasurementResolutionMap resolutions = {
69  {GeometryIdentifier(),
70  MeasurementResolution{MeasurementType::eLoc01, {0, 0}}}};
71 
72 // Construct initial track parameters.
73 CurvilinearTrackParameters makeParameters(double phi, double theta, double p,
74  double q) {
75  // create covariance matrix from reasonable standard deviations
76  Acts::BoundVector stddev;
77  stddev[Acts::eBoundLoc0] = 100_um;
78  stddev[Acts::eBoundLoc1] = 100_um;
79  stddev[Acts::eBoundTime] = 25_ns;
80  stddev[Acts::eBoundPhi] = 2_degree;
81  stddev[Acts::eBoundTheta] = 2_degree;
82  stddev[Acts::eBoundQOverP] = 1 / 100_GeV;
83  BoundSquareMatrix cov = stddev.cwiseProduct(stddev).asDiagonal();
84  // Let the particle starts from the origin
85  Vector4 mPos4(0., 0., 0., 0.);
86  return CurvilinearTrackParameters(mPos4, phi, theta, q / p, cov,
87  ParticleHypothesis::pionLike(std::abs(q)));
88 }
89 
90 std::default_random_engine rng(42);
91 
92 } // namespace
93 
94 BOOST_AUTO_TEST_CASE(trackparameters_estimation_test) {
95  // Construct a propagator with the cylinderal geometry and a constant magnetic
96  // field along z
97  Acts::Navigator navigator({
98  geometry,
99  true, // sensitive
100  true, // material
101  false // passive
102  });
103  auto field =
104  std::make_shared<Acts::ConstantBField>(Acts::Vector3(0.0, 0.0, 2._T));
105  ConstantFieldStepper stepper(std::move(field));
106 
107  ConstantFieldPropagator propagator(std::move(stepper), std::move(navigator));
108 
109  std::array<double, 2> pArray = {0.5_GeV, 1.0_GeV};
110  std::array<double, 3> phiArray = {20._degree, 0._degree - 20._degree};
111  std::array<double, 3> thetaArray = {80._degree, 90.0_degree, 100._degree};
112  std::array<double, 2> qArray = {1, -1};
113 
114  auto logger = Acts::getDefaultLogger("estimateTrackParamsFromSeed",
115  Acts::Logging::INFO);
116 
117  for (const auto& p : pArray) {
118  for (const auto& phi : phiArray) {
119  for (const auto& theta : thetaArray) {
120  for (const auto& q : qArray) {
121  BOOST_TEST_INFO("Test track with p = " << p << ", phi = " << phi
122  << ", theta = " << theta
123  << ", q = " << q);
124  auto start = makeParameters(phi, theta, p, q);
125  auto measurements = createMeasurements(propagator, geoCtx, magCtx,
126  start, resolutions, rng);
127 
128  // Create space points from different detector layers
129  std::map<GeometryIdentifier::Value, SpacePoint> spacePoints;
130  const Surface* bottomSurface = nullptr;
131  for (const auto& sl : measurements.sourceLinks) {
132  const auto geoId = sl.m_geometryId;
133  const auto& layer = geoId.layer();
134  auto it = spacePoints.find(layer);
135  // Avoid to use space point from the same layers
136  if (it != spacePoints.end()) {
137  continue;
138  }
139  const auto surface = geometry->findSurface(geoId);
140  const auto& localPos = sl.parameters;
141  Vector3 globalFakeMom(1, 1, 1);
142  Vector3 globalPos =
143  surface->localToGlobal(geoCtx, localPos, globalFakeMom);
144  // Create a space point (varianceR and varianceZ are lazily set to
145  // zero since they are not important for the test)
146  float r = std::hypot(globalPos.x(), globalPos.y());
147  spacePoints.emplace(
148  layer, SpacePoint{static_cast<float>(globalPos.x()),
149  static_cast<float>(globalPos.y()),
150  static_cast<float>(globalPos.z()), r,
151  static_cast<int>(geoId.layer()), 0., 0.});
152  if (spacePoints.size() == 1) {
153  bottomSurface = surface;
154  }
155  }
156 
157  // Check if there are at least 3 space points
158  if (spacePoints.size() < 3) {
159  BOOST_TEST_WARN("Number of space points less than 3.");
160  continue;
161  }
162 
163  // The truth track parameters at the bottom space point
164  const auto& expParams = measurements.truthParameters[0];
165  BOOST_TEST_INFO(
166  "The truth track parameters at the bottom space point: \n"
167  << expParams.transpose());
168  // The curvature of track projection on the transverse plane in unit
169  // of 1/mm
170  double rho = expParams[eBoundQOverP] * 0.3 * 2. / UnitConstants::m;
171 
172  // The space point pointers
173  std::array<const SpacePoint*, 3> spacePointPtrs{};
174  std::transform(spacePoints.begin(), std::next(spacePoints.begin(), 3),
175  spacePointPtrs.begin(),
176  [](const auto& sp) { return &sp.second; });
177 
178  // Test the partial track parameters estimator
179  auto partialParamsOpt = estimateTrackParamsFromSeed(
180  spacePointPtrs.begin(), spacePointPtrs.end(), *logger);
181  BOOST_REQUIRE(partialParamsOpt.has_value());
182  const auto& estPartialParams = partialParamsOpt.value();
183  BOOST_TEST_INFO(
184  "The estimated track parameters at the transverse plane: \n"
185  << estPartialParams.transpose());
186 
187  // The particle starting position is (0, 0, 0). Hence, d0 is zero; the
188  // phi at the point of cloest approach is exactly the phi of the truth
189  // particle
190  CHECK_CLOSE_ABS(estPartialParams[eBoundLoc0], 0., 1e-5);
191  CHECK_CLOSE_ABS(estPartialParams[eBoundPhi], phi, 1e-5);
192  CHECK_CLOSE_ABS(estPartialParams[eBoundQOverP], rho, 1e-4);
193  // The loc1, theta and time are set to zero in the estimator
194  CHECK_CLOSE_ABS(estPartialParams[eBoundLoc1], 0., 1e-10);
195  CHECK_CLOSE_ABS(estPartialParams[eBoundTheta], 0., 1e-10);
196  CHECK_CLOSE_ABS(estPartialParams[eBoundTime], 0., 1e-10);
197 
198  // Test the full track parameters estimator
199  auto fullParamsOpt = estimateTrackParamsFromSeed(
200  geoCtx, spacePointPtrs.begin(), spacePointPtrs.end(),
201  *bottomSurface, Vector3(0, 0, 2._T), 0.1_T, *logger);
202  BOOST_REQUIRE(fullParamsOpt.has_value());
203  const auto& estFullParams = fullParamsOpt.value();
204  BOOST_TEST_INFO(
205  "The estimated full track parameters at the bottom space point: "
206  "\n"
207  << estFullParams.transpose());
208 
209  CHECK_CLOSE_ABS(estFullParams[eBoundLoc0], expParams[eBoundLoc0],
210  1e-5);
211  CHECK_CLOSE_ABS(estFullParams[eBoundLoc1], expParams[eBoundLoc1],
212  1e-5);
213  // @todo Understand why the estimated phi has a limited precision
214  CHECK_CLOSE_ABS(estFullParams[eBoundPhi], expParams[eBoundPhi], 1e-1);
215  CHECK_CLOSE_ABS(estFullParams[eBoundTheta], expParams[eBoundTheta],
216  1e-2);
217  CHECK_CLOSE_ABS(estFullParams[eBoundQOverP], expParams[eBoundQOverP],
218  1e-2);
219  CHECK_CLOSE_ABS(estFullParams[eBoundTime], expParams[eBoundTime], 1.);
220  }
221  }
222  }
223  }
224 }