Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
JacobianEngineTests.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file JacobianEngineTests.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/TrackParametrization.hpp"
13 #include "Acts/Geometry/GeometryContext.hpp"
14 #include "Acts/Propagator/detail/JacobianEngine.hpp"
15 #include "Acts/Surfaces/PlaneSurface.hpp"
16 #include "Acts/Surfaces/Surface.hpp"
17 #include "Acts/Tests/CommonHelpers/FloatComparisons.hpp"
18 
19 #include <algorithm>
20 #include <cmath>
21 #include <memory>
22 
23 namespace tt = boost::test_tools;
24 
25 namespace Acts {
26 namespace Test {
27 
29 BOOST_AUTO_TEST_CASE(jacobian_engine_helper) {
30  // (1a) Standard test with curvilinear not glazingly close to z axis
31  Vector3 direction = Vector3(7., 8., 9.).normalized();
32  FreeToBoundMatrix f2cJacobian = detail::freeToCurvilinearJacobian(direction);
33 
34  ActsScalar phi = VectorHelpers::phi(direction);
35  ActsScalar theta = VectorHelpers::theta(direction);
36  ActsScalar sinPhi = std::sin(phi);
37  ActsScalar cosPhi = std::cos(phi);
38  ActsScalar sinTheta = std::sin(theta);
39  ActsScalar cosTheta = std::cos(theta);
40 
41  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc0, eFreePos0), -sinPhi, 1e-5);
42  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc0, eFreePos1), cosPhi, 1e-5);
43  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos0), -cosPhi * cosTheta, 1e-5);
44  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos1), -sinPhi * cosTheta, 1e-5);
45  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos2), sinTheta, 1e-5);
46  CHECK_CLOSE_REL(f2cJacobian(eBoundTime, eFreeTime), 1., 1e-5);
47  CHECK_CLOSE_REL(f2cJacobian(eBoundPhi, eFreeDir0), -sinPhi / sinTheta, 1e-5);
48  CHECK_CLOSE_REL(f2cJacobian(eBoundPhi, eFreeDir1), cosPhi / sinTheta, 1e-5);
49  CHECK_CLOSE_REL(f2cJacobian(eBoundTheta, eFreeDir0), cosPhi * cosTheta, 1e-5);
50  CHECK_CLOSE_REL(f2cJacobian(eBoundTheta, eFreeDir1), sinPhi * cosTheta, 1e-5);
51  CHECK_CLOSE_REL(f2cJacobian(eBoundTheta, eFreeDir2), -sinTheta, 1e-5);
52  CHECK_CLOSE_REL(f2cJacobian(eBoundQOverP, eFreeQOverP), 1., 1e-5);
53 
54  // (1b) Test with curvilinear glazingly close to z axis
55  direction = Vector3(1., 1., 999.).normalized();
56  f2cJacobian = detail::freeToCurvilinearJacobian(direction);
57 
58  phi = VectorHelpers::phi(direction);
59  theta = VectorHelpers::theta(direction);
60  sinPhi = std::sin(phi);
61  cosPhi = std::cos(phi);
62  sinTheta = std::sin(theta);
63  cosTheta = std::cos(theta);
64 
65  const ActsScalar c = std::hypot(direction.y(), direction.z());
66  const ActsScalar invC = 1. / c;
67  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc0, eFreePos1), -direction.z() * invC,
68  1e-5);
69  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc0, eFreePos2), direction.y() * invC,
70  1e-5);
71  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos0), c, 1e-5);
72  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos1),
73  -direction.x() * direction.y() * invC, 1e-5);
74  CHECK_CLOSE_REL(f2cJacobian(eBoundLoc1, eFreePos2),
75  -direction.x() * direction.z() * invC, 1e-5);
76 
77  // (2a) Standard test with curvilinear not glazingly close to z axis
78  direction = Vector3(7., 8., 9.).normalized();
79  BoundToFreeMatrix c2fJacobian = detail::curvilinearToFreeJacobian(direction);
80 
81  phi = VectorHelpers::phi(direction);
82  theta = VectorHelpers::theta(direction);
83  sinPhi = std::sin(phi);
84  cosPhi = std::cos(phi);
85  sinTheta = std::sin(theta);
86  cosTheta = std::cos(theta);
87 
88  CHECK_CLOSE_REL(c2fJacobian(eFreePos0, eBoundLoc0), -sinPhi, 1e-5);
89  CHECK_CLOSE_REL(c2fJacobian(eFreePos0, eBoundLoc1), -cosPhi * cosTheta, 1e-5);
90  CHECK_CLOSE_REL(c2fJacobian(eFreePos1, eBoundLoc0), cosPhi, 1e-5);
91  CHECK_CLOSE_REL(c2fJacobian(eFreePos1, eBoundLoc1), -sinPhi * cosTheta, 1e-5);
92  CHECK_CLOSE_REL(c2fJacobian(eFreePos2, eBoundLoc1), sinTheta, 1e-5);
93  // Time parameter: stays as is
94  CHECK_CLOSE_REL(c2fJacobian(eFreeTime, eBoundTime), 1, 1e-5);
95  CHECK_CLOSE_REL(c2fJacobian(eFreeDir0, eBoundPhi), -sinTheta * sinPhi, 1e-5);
96  CHECK_CLOSE_REL(c2fJacobian(eFreeDir0, eBoundTheta), cosTheta * cosPhi, 1e-5);
97  CHECK_CLOSE_REL(c2fJacobian(eFreeDir1, eBoundPhi), sinTheta * cosPhi, 1e-5);
98  CHECK_CLOSE_REL(c2fJacobian(eFreeDir1, eBoundTheta), cosTheta * sinPhi, 1e-5);
99  CHECK_CLOSE_REL(c2fJacobian(eFreeDir2, eBoundTheta), -sinTheta, 1e-5);
100  // Q/P parameter: stays as is
101  CHECK_CLOSE_REL(c2fJacobian(eFreeQOverP, eBoundQOverP), 1, 1e-5);
102 
103  // (3) Test angle - directio relational jacobians
104  direction = Vector3(2., 3., 4.).normalized();
105 
106  ActsMatrix<7, 8> d2aJacobian = detail::directionToAnglesJacobian(direction);
107  ActsMatrix<8, 7> a2dJacobian = detail::anglesToDirectionJacobian(direction);
108 
109  auto roundTrip = a2dJacobian * d2aJacobian;
110  BOOST_CHECK_EQUAL(roundTrip.cols(), 8);
111  BOOST_CHECK_EQUAL(roundTrip.rows(), 8);
112 }
113 
118 BOOST_AUTO_TEST_CASE(jacobian_engine_to_bound) {
119  // Create a test context
120  GeometryContext tgContext = GeometryContext();
121 
122  // Build a start vector
123  Vector3 position{1., 2., 3.};
124  double time = 4.;
125  Vector3 direction = Vector3(5., 2., 7.).normalized();
126  double qop = 0.125;
127 
128  // Build a surface
129  auto pSurface = Surface::makeShared<PlaneSurface>(position, direction);
130 
131  // Other rotated surface
132  Vector3 odirection = Vector3(6., 2., 8.).normalized();
133  auto oSurface = Surface::makeShared<PlaneSurface>(position, odirection);
134 
135  // The free parameter vector
136  FreeVector freeParameters;
137  freeParameters << position[0], position[1], position[2], time, direction[0],
138  direction[1], direction[2], qop;
139  // And its associated bound vector
140  BoundVector boundParameters;
141  boundParameters << 0., 0., VectorHelpers::phi(direction),
142  VectorHelpers::theta(direction), qop, time;
143 
144  // Build covariance matrices for bound and free case
145  BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
146  FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
147 
148  FreeMatrix noTransportJacobian = FreeMatrix::Identity();
149  FreeMatrix realTransportJacobian = 2 * FreeMatrix::Identity();
150 
151  FreeToPathMatrix freeToPathDerivatives =
152  pSurface->freeToPathDerivative(tgContext, freeParameters);
153  BoundToFreeMatrix boundToFreeJacobian =
154  pSurface->boundToFreeJacobian(tgContext, boundParameters);
155 
156  // (1) curvilinear/bound to bound transport jacobian
157  // a) test without actual transport into the same surface
158  BoundMatrix b2bTransportJacobian = detail::boundToBoundTransportJacobian(
159  tgContext, freeParameters, boundToFreeJacobian, noTransportJacobian,
160  freeToPathDerivatives, *pSurface);
161  BoundSquareMatrix newBoundCovariance =
162  b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
163  BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance));
164  // b) test without actual transport but to a new surface
165  b2bTransportJacobian = detail::boundToBoundTransportJacobian(
166  tgContext, freeParameters, boundToFreeJacobian, noTransportJacobian,
167  freeToPathDerivatives, *oSurface);
168  newBoundCovariance =
169  b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
170  BOOST_CHECK(not boundCovariance.isApprox(newBoundCovariance));
171  // c) test to "the same" surface with transport
172  // (not really senseful, but should give a different result)
173  b2bTransportJacobian = detail::boundToBoundTransportJacobian(
174  tgContext, freeParameters, boundToFreeJacobian, realTransportJacobian,
175  freeToPathDerivatives, *pSurface);
176  newBoundCovariance =
177  b2bTransportJacobian * boundCovariance * b2bTransportJacobian.transpose();
178  BOOST_CHECK(not boundCovariance.isApprox(newBoundCovariance));
179 
180  // (2) free to bound transport jacobian
181  const ActsMatrix<7, 8>& directionToAnglesJacobian =
182  detail::directionToAnglesJacobian(direction);
183  const ActsMatrix<8, 7>& anglesToDirectionJacobian =
184  detail::anglesToDirectionJacobian(direction);
185 
186  FreeToBoundMatrix f2bTransportJacobian = detail::freeToBoundTransportJacobian(
187  tgContext, freeParameters, directionToAnglesJacobian,
188  anglesToDirectionJacobian, noTransportJacobian, freeToPathDerivatives,
189  *pSurface);
190 
191  newBoundCovariance =
192  f2bTransportJacobian * freeCovariance * f2bTransportJacobian.transpose();
193  BOOST_CHECK(not boundCovariance.isApprox(newBoundCovariance));
194 }
195 
200 BOOST_AUTO_TEST_CASE(jacobian_engine_to_curvilinear) {
201  // Create a test context
202  GeometryContext tgContext = GeometryContext();
203 
204  // Build a start vector
205  Vector3 position{1., 2., 3.};
206  double time = 4.;
207  Vector3 direction = Vector3(5., 2., 7.).normalized();
208  double qop = 0.125;
209 
210  // Build a surface, starting surface for curvilinear
211  auto pSurface = Surface::makeShared<PlaneSurface>(position, direction);
212 
213  // The free parameter vector
214  FreeVector freeParameters;
215  freeParameters << position[0], position[1], position[2], time, direction[0],
216  direction[1], direction[2], qop;
217  // And its associated bound vector
218  BoundVector boundParameters;
219  boundParameters << 0., 0., VectorHelpers::phi(direction),
220  VectorHelpers::theta(direction), qop, time;
221 
222  // Build covariance matrices for bound and free case
223  BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
224  FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
225 
226  FreeMatrix noTransportJacobian = FreeMatrix::Identity();
227 
228  FreeToPathMatrix freeToPathDerivatives =
229  pSurface->freeToPathDerivative(tgContext, freeParameters);
230  BoundToFreeMatrix boundToFreeJacobian =
231  detail::curvilinearToFreeJacobian(direction);
232 
233  // (1) curvilinear/bound to curvilinear transport jacobian
234  // a) test without actual transport into the same surface
235  BoundMatrix b2cTransportJacobian =
236  detail::boundToCurvilinearTransportJacobian(
237  direction, boundToFreeJacobian, noTransportJacobian,
238  freeToPathDerivatives);
239  BoundSquareMatrix newBoundCovariance =
240  b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose();
241  BOOST_CHECK(boundCovariance.isApprox(newBoundCovariance));
242  // b) test to another curvilinear frame at the same point (no transport)
243  b2cTransportJacobian = detail::boundToCurvilinearTransportJacobian(
244  Vector3(4., 5., 6.).normalized(), boundToFreeJacobian,
245  noTransportJacobian, freeToPathDerivatives);
246  newBoundCovariance =
247  b2cTransportJacobian * boundCovariance * b2cTransportJacobian.transpose();
248  BOOST_CHECK(not boundCovariance.isApprox(newBoundCovariance));
249 
250  // (2) free to bound transport jacobian
251  const ActsMatrix<7, 8>& directionToAnglesJacobian =
252  detail::directionToAnglesJacobian(direction);
253  const ActsMatrix<8, 7>& anglesToDirectionJacobian =
254  detail::anglesToDirectionJacobian(direction);
255 
256  FreeToBoundMatrix f2cTransportJacobian =
257  detail::freeToCurvilinearTransportJacobian(
258  direction, directionToAnglesJacobian, anglesToDirectionJacobian,
259  noTransportJacobian, freeToPathDerivatives);
260 
261  newBoundCovariance =
262  f2cTransportJacobian * freeCovariance * f2cTransportJacobian.transpose();
263  BOOST_CHECK(not boundCovariance.isApprox(newBoundCovariance));
264 }
265 
270 BOOST_AUTO_TEST_CASE(jacobian_engine_to_free) {
271  // Create a test context
272  GeometryContext tgContext = GeometryContext();
273 
274  // Build a start vector
275  Vector3 position{1., 2., 3.};
276  double time = 4.;
277  Vector3 direction = Vector3(5., 2., 7.).normalized();
278  double qop = 0.125;
279 
280  // Build a surface, starting surface for curvilinear
281  auto pSurface = Surface::makeShared<PlaneSurface>(position, direction);
282 
283  // The free parameter vector
284  FreeVector freeParameters;
285  freeParameters << position[0], position[1], position[2], time, direction[0],
286  direction[1], direction[2], qop;
287  // And its associated bound vector
288  BoundVector boundParameters;
289  boundParameters << 0., 0., VectorHelpers::phi(direction),
290  VectorHelpers::theta(direction), qop, time;
291 
292  // Build covariance matrices for bound and free case
293  BoundSquareMatrix boundCovariance = 0.025 * BoundSquareMatrix::Identity();
294  FreeSquareMatrix freeCovariance = 0.025 * FreeSquareMatrix::Identity();
295 
296  FreeMatrix noTransportJacobian = FreeMatrix::Identity();
297 
298  BoundToFreeMatrix boundToFreeJacobian =
299  pSurface->boundToFreeJacobian(tgContext, boundParameters);
300 
301  // (1) bound to free
302  BoundToFreeMatrix b2fTransportJacobian = detail::boundToFreeTransportJacobian(
303  boundToFreeJacobian, noTransportJacobian);
304 
305  FreeMatrix newFreeCovariance1 =
306  b2fTransportJacobian * boundCovariance * b2fTransportJacobian.transpose();
307  BOOST_CHECK(not newFreeCovariance1.isApprox(freeCovariance));
308 
309  // (2) curvilinear to free
310  boundToFreeJacobian = detail::curvilinearToFreeJacobian(direction);
311  BoundToFreeMatrix c2fTransportJacobian = detail::boundToFreeTransportJacobian(
312  boundToFreeJacobian, noTransportJacobian);
313 
314  FreeMatrix newFreeCovariance2 =
315  c2fTransportJacobian * boundCovariance * c2fTransportJacobian.transpose();
316  BOOST_CHECK(not newFreeCovariance2.isApprox(freeCovariance));
317  // But those should be similar/equal
318  BOOST_CHECK(newFreeCovariance1.isApprox(newFreeCovariance2));
319 }
320 
321 } // namespace Test
322 } // namespace Acts