fix: Fix EigenStepper step upscaling (#3152)

The `EigenStepper` should use the same formula for up and downscaling the step size. Also according to https://cds.cern.ch/record/1113528/files/ATL-SOFT-PUB-2009-001.pdf we should use a factor 4 as a margin for rejecting the attempted step.

Core/include/Acts/Propagator/EigenStepper.ipp

@@ -158,8 +158,9 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
   // Runge-Kutta integrator state
   auto& sd = state.stepping.stepData;
 
-  double errorEstimate = 0.;
-  double h2 = 0, half_h = 0;
+  double errorEstimate = 0;
+  double h2 = 0;
+  double half_h = 0;
 
   auto pos = position(state.stepping);
   auto dir = direction(state.stepping);
@@ -178,7 +179,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
   }
 
   const auto calcStepSizeScaling = [&](const double errorEstimate_) -> double {
-    // For details about these values see ATL-SOFT-PUB-2009-001 for details
+    // For details about these values see ATL-SOFT-PUB-2009-001
     constexpr double lower = 0.25;
     constexpr double upper = 4.0;
     // This is given by the order of the Runge-Kutta method
@@ -202,6 +203,13 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
     return std::clamp(x, lower, upper);
   };
 
+  const auto isErrorTolerable = [&](const double errorEstimate_) {
+    // For details about these values see ATL-SOFT-PUB-2009-001
+    constexpr double marginFactor = 4.0;
+
+    return errorEstimate_ <= marginFactor * state.options.stepTolerance;
+  };
+
   // The following functor starts to perform a Runge-Kutta step of a certain
   // size, going up to the point where it can return an estimate of the local
   // integration error. The results are stated in the local variables above,
@@ -253,7 +261,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
     // Protect against division by zero
     errorEstimate = std::max(1e-20, errorEstimate);
 
-    return success(errorEstimate <= state.options.stepTolerance);
+    return success(isErrorTolerable(errorEstimate));
   };
 
   const double initialH =
@@ -272,7 +280,8 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
       break;
     }
 
-    h *= calcStepSizeScaling(2 * errorEstimate);
+    const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
+    h *= stepSizeScaling;
 
     // If step size becomes too small the particle remains at the initial
     // place
@@ -291,6 +300,8 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
 
   // When doing error propagation, update the associated Jacobian matrix
   if (state.stepping.covTransport) {
+    // using the direction before updated below
+
     // The step transport matrix in global coordinates
     FreeMatrix D;
     if (!state.stepping.extension.finalize(state, *this, navigator, h, D)) {
@@ -342,6 +353,7 @@ Acts::Result<double> Acts::EigenStepper<E, A>::step(
   (state.stepping.pars.template segment<3>(eFreeDir0)).normalize();
 
   if (state.stepping.covTransport) {
+    // using the updated direction
     state.stepping.derivative.template head<3>() =
         state.stepping.pars.template segment<3>(eFreeDir0);
     state.stepping.derivative.template segment<3>(4) = sd.k4;

Core/src/Propagator/SympyStepper.cpp

@@ -121,7 +121,7 @@ Result<double> SympyStepper::stepImpl(
   };
 
   const auto calcStepSizeScaling = [&](const double errorEstimate_) -> double {
-    // For details about these values see ATL-SOFT-PUB-2009-001 for details
+    // For details about these values see ATL-SOFT-PUB-2009-001
     constexpr double lower = 0.25;
     constexpr double upper = 4.0;
     // This is given by the order of the Runge-Kutta method
@@ -152,9 +152,11 @@ Result<double> SympyStepper::stepImpl(
 
   while (true) {
     nStepTrials++;
+
+    // For details about the factor 4 see ATL-SOFT-PUB-2009-001
     bool ok =
         rk4(pos.data(), dir.data(), t, h, qop, m, p_abs, getB, &errorEstimate,
-            state.pars.template segment<3>(eFreePos0).data(),
+            4 * stepTolerance, state.pars.template segment<3>(eFreePos0).data(),
             state.pars.template segment<3>(eFreeDir0).data(),
             state.pars.template segment<1>(eFreeTime).data(),
             state.derivative.data(),
@@ -164,7 +166,8 @@ Result<double> SympyStepper::stepImpl(
       break;
     }
 
-    h *= calcStepSizeScaling(2 * errorEstimate);
+    const double stepSizeScaling = calcStepSizeScaling(errorEstimate);
+    h *= stepSizeScaling;
 
     // If step size becomes too small the particle remains at the initial
     // place

Core/src/Propagator/codegen/sympy_stepper_math.hpp

@@ -15,8 +15,8 @@
 
 template <typename T, typename GetB>
 bool rk4(const T* p, const T* d, const T t, const T h, const T lambda,
-         const T m, const T p_abs, GetB getB, T* err, T* new_p, T* new_d,
-         T* new_time, T* path_derivatives, T* J) {
+         const T m, const T p_abs, GetB getB, T* err, const T errTol, T* new_p,
+         T* new_d, T* new_time, T* path_derivatives, T* J) {
   const auto B1 = getB(p);
   const auto x5 = std::pow(h, 2);
   const auto x0 = B1[1] * d[2];
@@ -79,7 +79,7 @@ bool rk4(const T* p, const T* d, const T t, const T h, const T lambda,
   *err =
       x5 * (std::fabs(-x29 + k2[0] + k3[0]) + std::fabs(-x30 + k2[1] + k3[1]) +
             std::fabs(-x31 + k2[2] + k3[2]));
-  if (*err > 1e-4) {
+  if (*err > errTol) {
     return false;
   }
   const auto x32 = (1.0 / 6.0) * x5;

Examples/Python/tests/root_file_hashes.txt

@@ -1,13 +1,13 @@
 test_pythia8__pythia8_particles.root: 49b89c3458a51aa9407f887be50e6bbcd9a2a0c897e6f2be5a5d6a29d1bf3505
 test_fatras__particles_simulation.root: bc970873fef0c2efd86ed5413623802353d2cd04abea72de14e8cdfc0e40076f
-test_fatras__hits.root: 821d9f3ce364bfea6c08df86b0204d7257cc1810f5b0bdfda45950e697ef5c52
+test_fatras__hits.root: 6e4beb045fa1712c4d14c280ba33c3fa13e4aff9de88d55c3e32f62ad226f724
 test_geant4__particles_simulation.root: a51f15555abefec7e9d9de31f9d856a3026af1b991db73fefaec5ae160a4af2f
 test_geant4__hits.root: adf5dcdf000a580412dc5089e17460897d6535c978eafa021584ba4281d0a1ac
-test_seeding__estimatedparams.root: 1fa2e879142059344e32ab62e3763dce14b5138b591c6e919073cb985782075f
+test_seeding__estimatedparams.root: d18718704a9ba5bf6627a0e4bf60581c640d80d66fc46dbd3e7c1d051ad99f71
 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880
 test_seeding__particles.root: 7855b021f39ad238bca098e4282667be0666f2d1630e5bcb9d51d3b5ee39fa14
 test_seeding__particles_simulation.root: 87d9c6c82422ca381a17735096b36c547eacb5dda2f26d7377614bd97a70ab1a
-test_seeding_orthogonal__estimatedparams.root: 77bf34e870bce0e44942525214f1101ace6e423a92be1d270f9991ffaf81450a
+test_seeding_orthogonal__estimatedparams.root: 974b9f45c43b81504c24439b47895a5d0a0501f205d2c377923d303fc10adc31
 test_seeding_orthogonal__performance_seeding.root: 60fbedcf5cb2b37cd8e526251940564432890d3a159d231ed819e915a904682c
 test_seeding_orthogonal__particles.root: 7855b021f39ad238bca098e4282667be0666f2d1630e5bcb9d51d3b5ee39fa14
 test_seeding_orthogonal__particles_simulation.root: 87d9c6c82422ca381a17735096b36c547eacb5dda2f26d7377614bd97a70ab1a
@@ -29,37 +29,37 @@ test_truth_tracking_kalman[odd-0.0]__performance_track_finder.root: 39aec6316cce
 test_truth_tracking_kalman[odd-1000.0]__trackstates_fitter.root: 72c79be1458c4f9c9a1661778c900f0875d257f2d391c4183a698825448919a1
 test_truth_tracking_kalman[odd-1000.0]__tracksummary_fitter.root: 3d424dec9b172f253c8c4ffbda470f678fd1081a3d36dcfea517ab0f94995ae4
 test_truth_tracking_kalman[odd-1000.0]__performance_track_finder.root: 39aec6316cceb90e314e16b02947faa691c18f57c3a851a25e547a8fc05a4593
-test_truth_tracking_gsf[generic]__trackstates_gsf.root: b843bcac04e7647aef716bd86f865ae72fa01448705ae933fd75ea5e340df95c
-test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 7de085a3b03136fc7ccf80410457a4ea30dcd702463e7c4e806b4a85a51c05ac
-test_truth_tracking_gsf[odd]__trackstates_gsf.root: ef12f504effc8dd00c378fde348291ce6d7e264c0f186b7ec82c8262f35c9f71
-test_truth_tracking_gsf[odd]__tracksummary_gsf.root: 07c72a4fe5383aeeada41ded42fc3337db67e1b75eccb9e85dfe783379869566
+test_truth_tracking_gsf[generic]__trackstates_gsf.root: 9456102b40bda129a05fbf8c27a585b45755c466e73e9a0a66519d28a93fb097
+test_truth_tracking_gsf[generic]__tracksummary_gsf.root: 618a47def602f4a9b387ed2c7eb0c33aa26ca91c0e177ca0fcef67c092d791f0
+test_truth_tracking_gsf[odd]__trackstates_gsf.root: 1b70e6e24af0c18f80d42c39e05dd20189bf4d1fcb48851f807cf723f461d7c0
+test_truth_tracking_gsf[odd]__tracksummary_gsf.root: ee806a7cddbc4bd41a74dceead569762954e8161af221ad60d637c6264e41752
 test_particle_gun__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
 test_material_mapping__material-map_tracks.root: b1138566d8d51579dce07f80166f05986942cf78c76b36875aea9b4b8b9bb957
 test_material_mapping__propagation-material.root: 2ec28ec8c10860ce0ca03b0b8abbfc9b52cef779d0113fddfe05949a8c362ec0
 test_volume_material_mapping__material-map-volume_tracks.root: 5092bbb52f15e4db7466e5a58c515794073e460af2411be2aa4ff001367928e4
 test_volume_material_mapping__propagation-volume-material.root: 89894f11a436038f6511d8b126f11ca2899d18b2b048a97c1b6c296ed85faddf
-test_digitization_example[smeared]__measurements.root: 67e17590e70c5dbc93f12070b35c085c831e9da1e5e6a223dc6e9d4c2dc7f9fa
-test_digitization_example[geometric]__measurements.root: 701b25440dae0c63a7d1a96afbf4dff98269e0c0cfa97a846cab3c146694760d
+test_digitization_example[smeared]__measurements.root: 03184f20b7b8e69f40424f76e2fd0620fb9f319d23ad69aa4bdd25c1b5cd482e
+test_digitization_example[geometric]__measurements.root: bc5d158478d7fdc4e09728cc9b4609afaa62be9dfcc1162685d03160f4d89495
 test_digitization_example_input[smeared]__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
-test_digitization_example_input[smeared]__measurements.root: 95ef20bcdc349da5bc025faf41cb76855ac087ae75df4f4bf953db5b23927aa5
+test_digitization_example_input[smeared]__measurements.root: 02a12e9949da9ba40b0120462d822f89e55562956fe3811981bdfd8648b7d946
 test_digitization_example_input[geometric]__particles.root: 5fe7dda2933ee6b9615b064d192322fe07831133cd998e5ed99a3b992b713a10
-test_digitization_example_input[geometric]__measurements.root: fa4729e28fdbbc459400dc3b7cc896c7e0cca047ab889c55c2f148d13b361637
-test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: e3828f4a3e07fcf4cd860310b4072e77aad79fd6366086e373413f67cc336960
-test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: 1d8c0517646660a82aad72f5ba5e16c70f089d120416342ceccce34fc2e3b425
+test_digitization_example_input[geometric]__measurements.root: 919abc40e61811bb1c591b6b5125594b922fd5a4bfd65d6b9f653e6a9292f662
+test_ckf_tracks_example[generic-full_seeding]__trackstates_ckf.root: f31fce48b9aaa951f05d35647a32a23d6f10dcf47fb7f0d639388082d6f4f2cb
+test_ckf_tracks_example[generic-full_seeding]__tracksummary_ckf.root: a5d11296d5a07ec44c8ce9696c1cfe140be05bc368bdef9e741a0c477ba0c755
 test_ckf_tracks_example[generic-full_seeding]__performance_seeding_trees.root: 0e0676ffafdb27112fbda50d1cf627859fa745760f98073261dcf6db3f2f991e
-test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: 9cf12f4cd8fa746a2a2309e5f73421c5b58825bc1e0454157e22fab7e81a0232
-test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: 74bc834cf220a16a96ce42a92eff3af1f0516362b68129ba98a85364237fd355
+test_ckf_tracks_example[generic-truth_estimated]__trackstates_ckf.root: 6ff08e4a2ec0acceda30d11f06cc3f4fe560418da866f57ad00a4e3a5ee22501
+test_ckf_tracks_example[generic-truth_estimated]__tracksummary_ckf.root: c5eb17a502614c2a2815797817d44ea420b2222267470f6a5341bb82d1bf01a2
 test_ckf_tracks_example[generic-truth_estimated]__performance_seeding.root: 1facb05c066221f6361b61f015cdf0918e94d9f3fce2269ec7b6a4dffeb2bc7e
-test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 0b6e072e23f4cd654e50218a84a0009c91295497b522f7a8906ed4eecf6e2c31
-test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 96ddae1bcab7b5bf79a4ba69e296420c2bed06aa7f0cc4f8d432c2271466c52d
-test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: cd96c9486eed8874752e90c82b48f684b9cb417640705d88626f9ae88c776c9e
-test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: a20288919a554b17bdf0ca9d44add8741ec2d3c2adaae57aba7e389f496b1104
+test_ckf_tracks_example[generic-truth_smeared]__trackstates_ckf.root: 89a284d518849eec0448e6cf41a2ca8d8deae3e9f8f913ca9ee3a17d422c8c0e
+test_ckf_tracks_example[generic-truth_smeared]__tracksummary_ckf.root: 472cd838f89f11498211f6f4d4412a21992fa1d76a511f3299e50cc5cb1d5ea3
+test_ckf_tracks_example[odd-full_seeding]__trackstates_ckf.root: 6ccb7b6a863744d477eb1ea420334c32d7867560019f1fdb3a651dc66b162596
+test_ckf_tracks_example[odd-full_seeding]__tracksummary_ckf.root: a39a1c798f16b0fc8a9cc2fc74e11b7550cc757db444d8518a41dbae2ae9a832
 test_ckf_tracks_example[odd-full_seeding]__performance_seeding_trees.root: 43c58577aafe07645e5660c4f43904efadf91d8cda45c5c04c248bbe0f59814f
-test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 3b63640eeefc84292e4e541069fca342dc420f01d9f008eafd965e5bddc1252d
-test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: 7ae4a98958c940ca62e2fa906a68486f9a274fae82bfa116ee4d4410ae160a75
+test_ckf_tracks_example[odd-truth_estimated]__trackstates_ckf.root: 04bbe3a9dc34a2e2106e8269188605c0b487702f3c3026dbaf974ec7619f4638
+test_ckf_tracks_example[odd-truth_estimated]__tracksummary_ckf.root: a6cbbcab9f9bd0e6614ade857de4a39b457d83c3297553b0a77cc9f5a7cb78b7
 test_ckf_tracks_example[odd-truth_estimated]__performance_seeding.root: 1a36b7017e59f1c08602ef3c2cb0483c51df248f112e3780c66594110719c575
-test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 3e860f9ddec63e9f1bebedb89c67f07dd50c937c2beb6d84dcc8f4d43506c78b
-test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: 4ae259c4a2bbd4b3801cd552cc3b1a114f1ce3660183fa7fc7518d7d04484a45
+test_ckf_tracks_example[odd-truth_smeared]__trackstates_ckf.root: 22e6d12dcbdca5bc7492c5e073cd7759b61653f9cd5eb43a26a8a3525424f62a
+test_ckf_tracks_example[odd-truth_smeared]__tracksummary_ckf.root: e2aaf2bdafb5757698ee7a1db542c4612ab35ea305dbfb7908dc204974a85ffe
 test_vertex_fitting_reading[Truth-False-100]__performance_vertexing.root: 76ef6084d758dfdfc0151ddec2170e12d73394424e3dac4ffe46f0f339ec8293
 test_vertex_fitting_reading[Iterative-False-100]__performance_vertexing.root: 60372210c830a04f95ceb78c6c68a9b0de217746ff59e8e73053750c837b57eb
 test_vertex_fitting_reading[Iterative-True-100]__performance_vertexing.root: e34f217d524a5051dbb04a811d3407df3ebe2cc4bb7f54f6bda0847dbd7b52c3
@@ -73,10 +73,10 @@ test_root_prop_step_writer[kwargsConstructor]__prop_steps.root: 02460dea7b0b56b5
 test_root_particle_writer[configPosConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
 test_root_particle_writer[configKwConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
 test_root_particle_writer[kwargsConstructor]__particles.root: 759f7c0225a30afc28c3dda6010295207254ebb0194b88f692d538b8c65d4326
-test_root_meas_writer__meas.root: 347e5815f00859edbdae5a4d3053e3c797b793975636e1b6db20e192600d29d5
-test_root_simhits_writer[configPosConstructor]__meas.root: 9610837e081f3970446ca78d4c90ce116a611285dd1f4c646548dbece45abfa2
-test_root_simhits_writer[configKwConstructor]__meas.root: 9610837e081f3970446ca78d4c90ce116a611285dd1f4c646548dbece45abfa2
-test_root_simhits_writer[kwargsConstructor]__meas.root: 9610837e081f3970446ca78d4c90ce116a611285dd1f4c646548dbece45abfa2
+test_root_meas_writer__meas.root: c92ef6936a7a20d3eb09722e9fab6aa91286439ff57b5c89c9f75c1fd156684e
+test_root_simhits_writer[configPosConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
+test_root_simhits_writer[configKwConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
+test_root_simhits_writer[kwargsConstructor]__meas.root: 66fee65a672bbf28098aac7af59ef181c2ac8331f87a379132873da86bdb1ba5
 test_root_material_writer__material.root: e3b0c44298fc1c149afbf4c8996fb92427ae41e4649b934ca495991b7852b855
 test_root_clusters_writer[configPosConstructor]__clusters.root: e842df4fe04eefff3df5f32cd1026e93286be62b8040dc700a2aff557c56dec8
 test_root_clusters_writer[configKwConstructor]__clusters.root: e842df4fe04eefff3df5f32cd1026e93286be62b8040dc700a2aff557c56dec8

Tests/UnitTests/Core/Propagator/DirectNavigatorTests.cpp

@@ -78,7 +78,7 @@ bool referenceTiming = false;
 bool oversteppingTest = false;
 double oversteppingMaxStepSize = 1_mm;
 
-/// The actual test nethod that runs the test
+/// The actual test method that runs the test
 /// can be used with several propagator types
 ///
 /// @tparam rpropagator_t is the reference propagator type

codegen/generate_sympy_stepper.py

@@ -195,7 +195,7 @@ def my_step_function_print(name_exprs, run_cse=True):
 
     lines = []
 
-    head = "template <typename T, typename GetB> bool rk4(const T* p, const T* d, const T t, const T h, const T lambda, const T m, const T p_abs, GetB getB, T* err, T* new_p, T* new_d, T* new_time, T* path_derivatives, T* J) {"
+    head = "template <typename T, typename GetB> bool rk4(const T* p, const T* d, const T t, const T h, const T lambda, const T m, const T p_abs, GetB getB, T* err, const T errTol, T* new_p, T* new_d, T* new_time, T* path_derivatives, T* J) {"
     lines.append(head)
 
     lines.append("  const auto B1 = getB(p);")
@@ -215,7 +215,7 @@ def post_expr_hook(var):
         if str(var) == "p3":
             return "const auto B3 = getB(p3);"
         if str(var) == "err":
-            return "if (*err > 1e-4) {\n  return false;\n}"
+            return "if (*err > errTol) {\n  return false;\n}"
         if str(var) == "new_time":
             return "if (J == nullptr) {\n  return true;\n}"
         if str(var) == "new_J":