cpu/propagator/jacobian_validation.cpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 // Project include(s)
0010 #include "detray/definitions/algebra.hpp"
0011 #include "detray/definitions/units.hpp"
0012 #include "detray/geometry/tracking_surface.hpp"
0013 #include "detray/navigation/caching_navigator.hpp"
0014 #include "detray/navigation/intersection/helix_intersector.hpp"
0015 #include "detray/propagator/actors/parameter_updater.hpp"
0016 #include "detray/propagator/propagator.hpp"
0017 #include "detray/propagator/rk_stepper.hpp"
0018
0019 // Detray test include(s)
0020 #include "detray/test/common/bfield.hpp"
0021 #include "detray/test/common/build_telescope_detector.hpp"
0022 #include "detray/test/common/track_generators.hpp"
0023 #include "detray/test/framework/types.hpp"
0024 #include "detray/test/utils/inspectors.hpp"
0025 #include "detray/test/utils/statistics.hpp"
0026
0027 // Vecmem include(s)
0028 #include <vecmem/memory/host_memory_resource.hpp>
0029
0030 // GTest include(s)
0031 #include <gtest/gtest.h>
0032
0033 // Boost
0034 #include "detray/options/boost_program_options.hpp"
0035
0036 // System include(s).
0037 #include <algorithm>
0038 #include <filesystem>
0039 #include <fstream>
0040 #include <iomanip>
0041 #include <iostream>
0042
0043 namespace po = boost::program_options;
0044 using namespace detray;
0045
0046 // Type declarations
0047 using test_algebra = test::algebra;
0048 using scalar = dscalar<test_algebra>;
0049 using vector3 = dvector3D<test_algebra>;
0050 using transform3_type = dtransform3D<test_algebra>;
0051 using bound_param_vector_type =
0052     bound_track_parameters<test_algebra>::parameter_vector_type;
0053 using bound_covariance_type =
0054     bound_track_parameters<test_algebra>::covariance_type;
0055 template <std::size_t ROWS, std::size_t COLS>
0056 using matrix_type = dmatrix<test_algebra, ROWS, COLS>;
0057
0058 namespace {
0059
0060 // B-field Setup
0061 // 1.996 T from averaged Bz strength of ODD field in the region of sqrt(x^2 +
0062 // y^2) < 500 mm and abs(z) < 500 mm
0063 const vector3 B_z{0.f, 0.f, 1.996f * unit<scalar>::T};
0064
0065 // Initial delta for numerical differentiaion
0066 const std::array<scalar, 5u> h_sizes_rect{1e0f, 1e0f, 2e-2f, 1e-3f, 1e-3f};
0067 const std::array<scalar, 5u> h_sizes_wire{1e0f, 1e0f, 2e-2f, 1e-3f, 1e-3f};
0068
0069 // Ridders' algorithm setup
0070 constexpr const unsigned int Nt = 100u;
0071 const std::array<scalar, 5u> safe{5.0f, 5.0f, 5.0f, 5.0f, 5.0f};
0072 const std::array<scalar, 5u> con{1.1f, 1.1f, 1.1f, 1.1f, 1.1f};
0073 constexpr const scalar big = std::numeric_limits<scalar>::max();
0074
0075 std::random_device rd;
0076 // For detector generation
0077 std::mt19937_64 mt1(rd());
0078 // For smearing initial parameter
0079 std::mt19937_64 mt2(rd());
0080
0081 // Momentum range
0082 constexpr const scalar min_mom = 0.5f * unit<scalar>::GeV;
0083 constexpr const scalar max_mom = 100.f * unit<scalar>::GeV;
0084
0085 // Detector length generator
0086 constexpr const scalar min_detector_length = 50.f * unit<scalar>::mm;
0087 constexpr const scalar max_detector_length = 500.f * unit<scalar>::mm;
0088 std::uniform_real_distribution rand_length(min_detector_length,
0089                                            max_detector_length);
0090 constexpr const scalar envelope_size = 2000.f * unit<scalar>::mm;
0091
0092 // Mask size scaler
0093 constexpr const scalar mask_scaler = 1.5f;
0094
0095 // Euler angles for the surface rotation
0096 std::uniform_real_distribution<scalar> rand_alpha(0.f,
0097                                                   2.f * constant<scalar>::pi);
0098 std::uniform_real_distribution<scalar> rand_cosbeta(constant<scalar>::inv_sqrt2,
0099                                                     1.f);
0100 std::uniform_int_distribution rand_bool(0, 1);
0101 std::uniform_real_distribution<scalar> rand_gamma(0.f,
0102                                                   2.f * constant<scalar>::pi);
0103
0104 // Particle setup = muon
0105 pdg_particle ptc = muon<scalar>();
0106
0107 // Correlation factor in the range of [-10%, 10%]
0108 constexpr scalar min_corr = -0.1f;
0109 constexpr scalar max_corr = 0.1f;
0110
0111 // Values for sampling standard deviation
0112 const std::array<scalar, 6u> stddevs_sampling{50.f * unit<scalar>::um,
0113                                               50.f * unit<scalar>::um,
0114                                               1.f * unit<scalar>::mrad,
0115                                               1.f * unit<scalar>::mrad,
0116                                               0.01f,
0117                                               1.f * unit<scalar>::ns};
0118
0119 // surface types
0120 using rect_type = rectangle2D;
0121 using wire_type = line_square;
0122
0123 }  // namespace
0124
0125 // Preprocess delta
0126 void preprocess_delta(const unsigned int i, scalar& delta,
0127                       const bound_track_parameters<test_algebra> ref_param) {
0128   if (i == e_bound_theta) {
0129     const scalar rtheta = ref_param.theta();
0130     if (rtheta < constant<scalar>::pi_2) {
0131       delta = math::min(delta, rtheta);
0132     } else if (rtheta >= constant<scalar>::pi_2) {
0133       delta = math::min(delta, constant<scalar>::pi - rtheta);
0134     }
0135   }
0136 }
0137
0138 // Containers for Ridders algorithm
0139 // Page 231 of [Numerical Recipes] 3rd edition
0140 struct ridders_derivative {
0141   std::array<std::array<std::array<scalar, Nt>, Nt>, 5u> Arr{};
0142   std::array<scalar, 5u> fac{};
0143   std::array<scalar, 5u> errt{};
0144   std::array<scalar, 5u> err{big, big, big, big, big};
0145   std::array<bool, 5u> complete{false, false, false, false, false};
0146
0147   void initialize(const bound_param_vector_type& nvec1,
0148                   const bound_param_vector_type& nvec2, const scalar delta) {
0149     for (unsigned int j = 0; j < 5u; j++) {
0150       Arr[j][0][0] = (nvec1[j] - nvec2[j]) / (2.f * delta);
0151     }
0152   }
0153
0154   void run(const bound_param_vector_type& nvec1,
0155            const bound_param_vector_type& nvec2, const scalar delta,
0156            const unsigned int p, const unsigned int i,
0157            bound_covariance_type& differentiated_jacobian) {
0158     for (unsigned int j = 0; j < 5u; j++) {
0159       Arr[j][0][p] = (nvec1[j] - nvec2[j]) / (2.f * delta);
0160     }
0161
0162     const scalar con2 = con[i] * con[i];
0163     fac[i] = con2;
0164
0165     for (unsigned int q = 1; q <= p; q++) {
0166       for (unsigned int j = 0; j < 5u; j++) {
0167         Arr[j][q][p] = (Arr[j][q - 1][p] * fac[i] - Arr[j][q - 1][p - 1]) /
0168                        (fac[i] - 1.0f);
0169         fac[i] = con2 * fac[i];
0170
0171         errt[j] = math::max(math::fabs(Arr[j][q][p] - Arr[j][q - 1][p]),
0172                             math::fabs(Arr[j][q][p] - Arr[j][q - 1][p - 1]));
0173
0174         if (errt[j] <= err[j]) {
0175           if (complete[j] == false) {
0176             err[j] = errt[j];
0177             getter::element(differentiated_jacobian, j, i) = Arr[j][q][p];
0178             /*
0179             // Please leave this for debug
0180             if (j == e_bound_theta && i == e_bound_loc0) {
0181                 std::clog << q << " " << p << " "
0182                           << getter::element(
0183                                  differentiated_jacobian, j, i)
0184                           << "  " << math::fabs(Arr[j][q][p])
0185                           << std::endl;
0186             }
0187             */
0188           }
0189         }
0190       }
0191     }
0192
0193     for (unsigned int j = 0; j < 5u; j++) {
0194       /*
0195       // Please leave this for debug
0196       if (j == e_bound_loc0 && i == e_bound_theta) {
0197           std::clog << getter::element(differentiated_jacobian, j, i)
0198                     << "  " << Arr[j][p][p] << "  "
0199                     << Arr[j][p - 1][p - 1] << "  "
0200                     << math::fabs(Arr[j][p][p] - Arr[j][p - 1][p - 1])
0201                     << "  " << safe[i] * err[j] << std::endl;
0202       }
0203       */
0204       if (math::fabs(Arr[j][p][p] - Arr[j][p - 1][p - 1]) >= safe[i] * err[j]) {
0205         complete[j] = true;
0206       }
0207     }
0208   }
0209
0210   bool finished() { return (std::ranges::count(complete, false) == 0u); }
0211 };
0212
0213 void wrap_angles(const bound_param_vector_type& ref_vector,
0214                  bound_param_vector_type& target_vector) {
0215   const scalar rphi = ref_vector.phi();
0216   const scalar tphi = target_vector.phi();
0217   scalar new_tphi = tphi;
0218
0219   if (rphi >= constant<scalar>::pi_2) {
0220     if (tphi < 0.f) {
0221       new_tphi += 2.f * constant<scalar>::pi;
0222     }
0223   } else if (rphi < -constant<scalar>::pi_2) {
0224     if (tphi >= 0.f) {
0225       new_tphi -= 2.f * constant<scalar>::pi;
0226     }
0227   }
0228
0229   target_vector.set_phi(new_tphi);
0230 }
0231
0232 scalar get_relative_difference(scalar ref_val, scalar num_val) {
0233   scalar rel_diff{0.f};
0234
0235   // If the evaluated jacovian or numerical differentiation is too small set the
0236   // relative difference to zero
0237   if (ref_val == 0.f || num_val == 0.f) {
0238     rel_diff = 0.f;
0239   } else {
0240     rel_diff = std::abs(ref_val - num_val) / std::abs(num_val);
0241   }
0242
0243   return rel_diff;
0244 }
0245
0246 // Get random initial covariance
0247 bound_covariance_type get_random_initial_covariance(const scalar ini_qop) {
0248   // Initial covariance matrix for smearing
0249   auto ini_cov = matrix::zero<bound_covariance_type>();
0250
0251   // Random correction factor
0252   std::uniform_real_distribution<scalar> rand_corr(min_corr, max_corr);
0253
0254   // Distribution for sampling standard deviations
0255   std::normal_distribution<scalar> rand_l0(0.f * unit<scalar>::um,
0256                                            stddevs_sampling[0u]);
0257   std::normal_distribution<scalar> rand_l1(0.f * unit<scalar>::um,
0258                                            stddevs_sampling[1u]);
0259   std::normal_distribution<scalar> rand_phi(0.f * unit<scalar>::mrad,
0260                                             stddevs_sampling[2u]);
0261   std::normal_distribution<scalar> rand_theta(0.f * unit<scalar>::mrad,
0262                                               stddevs_sampling[3u]);
0263   std::normal_distribution<scalar> rand_qop(0.f,
0264                                             stddevs_sampling[4u] * ini_qop);
0265   std::normal_distribution<scalar> rand_time(0.f * unit<scalar>::ns,
0266                                              stddevs_sampling[5u]);
0267
0268   std::array<scalar, 6u> stddevs{};
0269   stddevs[0] = rand_l0(mt2);
0270   stddevs[1] = rand_l1(mt2);
0271   stddevs[2] = rand_phi(mt2);
0272   stddevs[3] = rand_theta(mt2);
0273   stddevs[4] = rand_qop(mt2);
0274   stddevs[5] = rand_time(mt2);
0275
0276   for (unsigned int i = 0u; i < 6u; i++) {
0277     for (unsigned int j = 0u; j < 6u; j++) {
0278       if (i == j) {
0279         getter::element(ini_cov, i, i) = stddevs[i] * stddevs[i];
0280       } else if (i > j) {
0281         getter::element(ini_cov, i, j) =
0282             std::abs(stddevs[i]) * std::abs(stddevs[j]) * rand_corr(mt2);
0283         getter::element(ini_cov, j, i) = getter::element(ini_cov, i, j);
0284       }
0285     }
0286   }
0287
0288   return ini_cov;
0289 }
0290
0291 // Input covariance should be the diagonal matrix
0292 auto get_smeared_bound_vector(const bound_covariance_type& ini_cov,
0293                               const bound_param_vector_type& ini_vec) {
0294   using bound_vector_type = bound_vector<test_algebra>;
0295
0296   // Do the Cholesky Decomposition
0297   const bound_covariance_type L = matrix::cholesky_decomposition(ini_cov);
0298
0299   // Vecor with random elements from a normal distribution
0300   auto k = matrix::zero<bound_vector_type>();
0301   std::normal_distribution<scalar> normal_dist(0.f, 1.f);
0302   for (unsigned int i = 0u; i < 5u; i++) {
0303     // Smear the value
0304     getter::element(k, i, 0) = normal_dist(mt2);
0305   }
0306
0307   const bound_vector_type new_vec = ini_vec.vector() + L * k;
0308
0309   return bound_param_vector_type{new_vec};
0310 }
0311
0312 template <typename detector_t, typename detector_t::metadata::mask_id mask_id>
0313 std::pair<euler_rotation<test_algebra>, std::array<scalar, 3u>> tilt_surface(
0314     detector_t& det, const unsigned int sf_id, const vector3& helix_dir,
0315     const scalar alpha, const scalar beta, const scalar gamma) {
0316   const auto& sf = det.surface(sf_id);
0317   const auto& trf_link = sf.transform();
0318   auto& trf = det.transform_store().at(trf_link);
0319
0320   euler_rotation<test_algebra> euler;
0321   euler.alpha = alpha;
0322
0323   if (sf_id == 1u) {
0324     euler.beta = beta;
0325     euler.gamma = gamma;
0326   }
0327
0328   // Helix direction
0329   euler.z = helix_dir;
0330
0331   if constexpr (mask_id == detector_t::masks::id::e_rectangle2D) {
0332     // ubasis = trf.x() for bound frame
0333     euler.x = trf.x();
0334   } else if (mask_id == detector_t::masks::id::e_drift_cell) {
0335     // ubasis = vxt/|vxt| where v is trf.z() and t is helix_dir
0336     EXPECT_NEAR(vector::dot(trf.z(), helix_dir), 0.f, 1e-6f);
0337     euler.x = vector::normalize(vector::cross(trf.z(), helix_dir));
0338   }
0339
0340   EXPECT_NEAR(vector::dot(euler.x, euler.z), 0.f, 1e-6f);
0341
0342   auto [local_x, local_z] = euler();
0343
0344   if constexpr (mask_id == detector_t::masks::id::e_drift_cell) {
0345     // local_z should be the line direction
0346     local_z = vector::cross(local_z, local_x);
0347   }
0348
0349   // At the moment, we do not shift the surfaces
0350   scalar x_shift{0.f};
0351   scalar y_shift{0.f};
0352   scalar z_shift{0.f};
0353
0354   // Translation vector
0355   vector3 translation = trf.translation();
0356
0357   vector3 displacement({x_shift, y_shift, z_shift});
0358   translation = trf.translation() + displacement;
0359   transform3_type new_trf(translation, local_z, local_x, true);
0360   // @todo : Remove and use different context
0361   detray::detail::set_transform(det, new_trf, trf_link);
0362
0363   return {euler, {x_shift, y_shift, z_shift}};
0364 }
0365
0366 template <concepts::algebra algebra_t>
0367 struct bound_getter : public base_actor {
0368   // Track types
0369   using bound_track_parameters_type = bound_track_parameters<algebra_t>;
0370   using free_track_parameters_type = free_track_parameters<algebra_t>;
0371
0372   struct state {
0373     scalar m_min_path_length{detray::detail::invalid_value<scalar>()};
0374     scalar m_path_length{0.f};
0375     scalar m_abs_path_length{0.f};
0376     bound_track_parameters_type m_param_departure{};
0377     bound_track_parameters_type m_param_destination{};
0378     typename bound_track_parameters_type::covariance_type m_jacobi{};
0379     scalar m_avg_step_size{0.f};
0380     std::size_t step_count{0u};
0381     std::size_t track_ID{0u};
0382   };
0383
0384   template <typename propagator_state_t>
0385   DETRAY_HOST_DEVICE void operator()(
0386       state& actor_state, propagator_state_t& propagation,
0387       const detray::actor::parameter_transporter_result<algebra_t>& res) const {
0388     auto& navigation = propagation.navigation();
0389     auto& stepping = propagation.stepping();
0390
0391     actor_state.step_count++;
0392
0393     const scalar N = static_cast<scalar>(actor_state.step_count);
0394
0395     actor_state.m_avg_step_size =
0396         ((N - 1.f) * actor_state.m_avg_step_size + stepping.step_size()) / N;
0397
0398     // Warning for too many step counts
0399     if (actor_state.step_count > 1000000) {
0400       std::clog << "Too many step counts!" << std::endl;
0401       std::clog << "Track ID: " << actor_state.track_ID << std::endl;
0402       std::clog << "Path length: " << actor_state.m_path_length << std::endl;
0403       std::clog << "PhiI: " << actor_state.m_param_departure.phi() << std::endl;
0404       std::clog << "ThetaI: " << actor_state.m_param_departure.theta()
0405                 << std::endl;
0406       std::clog << "QopI: " << actor_state.m_param_departure.qop() << std::endl;
0407       navigation.exit();
0408       propagation.heartbeat(false);
0409     }
0410
0411     if ((navigation.is_on_sensitive() || navigation.is_on_passive()) &&
0412         navigation.geometry_identifier().index() == 0u) {
0413       actor_state.m_param_departure = res.destination_params();
0414     }
0415     // Get the bound track parameters and jacobian at the destination
0416     // surface
0417     else if ((navigation.is_on_sensitive() || navigation.is_on_passive()) &&
0418              navigation.geometry_identifier().index() == 1u) {
0419       actor_state.m_path_length = stepping.path_length();
0420       actor_state.m_abs_path_length = stepping.abs_path_length();
0421       actor_state.m_param_destination = res.destination_params();
0422       actor_state.m_jacobi =
0423           actor::parameter_transporter<algebra_t>().get_full_jacobian(
0424               propagation, actor_state.m_param_departure);
0425
0426       // Stop navigation if the destination surface is found
0427       navigation.exit();
0428       propagation.heartbeat(false);
0429     }
0430
0431     if (stepping.path_length() > actor_state.m_min_path_length) {
0432       propagation.navigation().set_no_trust();
0433     }
0434
0435     return;
0436   }
0437 };
0438
0439 /// Numerically integrate the jacobian
0440 template <typename propagator_t, typename field_t>
0441 bound_getter<test_algebra>::state evaluate_bound_param(
0442     const std::size_t trk_count, const scalar detector_length,
0443     const bound_track_parameters<test_algebra>& initial_param,
0444     const typename propagator_t::detector_type& det, const field_t& field,
0445     const scalar overstep_tolerance, const scalar path_tolerance,
0446     const scalar rk_tolerance, const scalar constraint_step,
0447     bool use_field_gradient, bool do_covariance_transport, bool do_inspect) {
0448   // Propagator is built from the stepper and navigator
0449   propagation::config cfg{};
0450   cfg.navigation.intersection.overstep_tolerance =
0451       static_cast<float>(overstep_tolerance);
0452   cfg.navigation.intersection.path_tolerance =
0453       static_cast<float>(path_tolerance);
0454   cfg.navigation.estimate_scattering_noise = false;
0455   cfg.stepping.rk_error_tol = static_cast<float>(rk_tolerance);
0456   cfg.stepping.use_eloss_gradient = true;
0457   cfg.stepping.use_field_gradient = use_field_gradient;
0458   cfg.stepping.do_covariance_transport = do_covariance_transport;
0459   propagator_t p(cfg);
0460
0461   // Actor states
0462   actor::parameter_updater_state<test_algebra> updater_state{cfg,
0463                                                              initial_param};
0464   bound_getter<test_algebra>::state bound_getter_state{};
0465   bound_getter_state.track_ID = trk_count;
0466   bound_getter_state.m_min_path_length = detector_length * 0.75f;
0467   auto actor_states = detray::tie(updater_state, bound_getter_state);
0468
0469   // Init propagator states for the reference track
0470   typename propagator_t::stepper_type::magnetic_field_type bfield_view(field);
0471   typename propagator_t::state state(initial_param, bfield_view, det);
0472
0473   // Run the propagation for the reference track
0474   state.set_particle(ptc);
0475   state.debug(do_inspect);
0476   state.stepping()
0477       .template set_constraint<detray::step::constraint::e_accuracy>(
0478           static_cast<float>(constraint_step));
0479
0480   p.propagate(state, actor_states);
0481
0482   return bound_getter_state;
0483 }
0484
0485 template <typename propagator_t, typename field_t>
0486 bound_param_vector_type get_displaced_bound_vector(
0487     const std::size_t trk_count,
0488     const bound_track_parameters<test_algebra>& ref_param,
0489     const typename propagator_t::detector_type& det,
0490     const scalar detector_length, const field_t& field,
0491     const scalar overstep_tolerance, const scalar path_tolerance,
0492     const scalar rk_tolerance, const scalar constraint_step,
0493     const unsigned int target_index, const scalar displacement) {
0494   propagation::config cfg{};
0495   cfg.navigation.intersection.overstep_tolerance =
0496       static_cast<float>(overstep_tolerance);
0497   cfg.navigation.intersection.path_tolerance =
0498       static_cast<float>(path_tolerance);
0499   cfg.navigation.estimate_scattering_noise = false;
0500   cfg.stepping.rk_error_tol = static_cast<float>(rk_tolerance);
0501   cfg.stepping.do_covariance_transport = false;
0502
0503   // Propagator is built from the stepper and navigator
0504   propagator_t p(cfg);
0505
0506   bound_track_parameters<test_algebra> dparam = ref_param;
0507   dparam[target_index] += displacement;
0508
0509   typename propagator_t::stepper_type::magnetic_field_type bfield_view(field);
0510   typename propagator_t::state dstate(dparam, bfield_view, det);
0511
0512   // Actor states
0513   actor::parameter_updater_state<test_algebra> updater_state{cfg, dparam};
0514   bound_getter<test_algebra>::state bound_getter_state{};
0515   bound_getter_state.track_ID = trk_count;
0516   bound_getter_state.m_min_path_length = detector_length * 0.75f;
0517
0518   auto actor_states = detray::tie(updater_state, bound_getter_state);
0519   dstate.set_particle(ptc);
0520   dstate.stepping()
0521       .template set_constraint<detray::step::constraint::e_accuracy>(
0522           constraint_step);
0523
0524   p.propagate(dstate, actor_states);
0525
0526   bound_param_vector_type new_vec = bound_getter_state.m_param_destination;
0527
0528   // phi needs to be wrapped w.r.t. phi of the reference vector
0529   wrap_angles(ref_param, new_vec);
0530
0531   return new_vec;
0532 }
0533
0534 /// Numerically evaluate the jacobian
0535 template <typename propagator_t, typename field_t>
0536 bound_track_parameters<test_algebra>::covariance_type directly_differentiate(
0537     const std::size_t trk_count,
0538     const bound_track_parameters<test_algebra>& ref_param,
0539     const typename propagator_t::detector_type& det,
0540     const scalar detector_length, const field_t& field,
0541     const scalar overstep_tolerance, const scalar path_tolerance,
0542     const scalar rk_tolerance, const scalar constraint_step,
0543     const std::array<scalar, 5u> hs,
0544     std::array<unsigned int, 5u>& num_iterations,
0545     std::array<bool, 25>& convergence) {
0546   // Return Jacobian
0547   bound_covariance_type differentiated_jacobian;
0548
0549   for (unsigned int i = 0u; i < 5u; i++) {
0550     scalar delta = hs[i];
0551     preprocess_delta(i, delta, ref_param);
0552
0553     const auto vec1 = get_displaced_bound_vector<propagator_t, field_t>(
0554         trk_count, ref_param, det, detector_length, field, overstep_tolerance,
0555         path_tolerance, rk_tolerance, constraint_step, i, 1.f * delta);
0556     const auto vec2 = get_displaced_bound_vector<propagator_t, field_t>(
0557         trk_count, ref_param, det, detector_length, field, overstep_tolerance,
0558         path_tolerance, rk_tolerance, constraint_step, i, -1.f * delta);
0559
0560     ridders_derivative ridder;
0561     ridder.initialize(vec1, vec2, delta);
0562
0563     for (unsigned int p = 1u; p < Nt; p++) {
0564       delta /= con[i];
0565
0566       const auto nvec1 = get_displaced_bound_vector<propagator_t, field_t>(
0567           trk_count, ref_param, det, detector_length, field, overstep_tolerance,
0568           path_tolerance, rk_tolerance, constraint_step, i, 1.f * delta);
0569       const auto nvec2 = get_displaced_bound_vector<propagator_t, field_t>(
0570           trk_count, ref_param, det, detector_length, field, overstep_tolerance,
0571           path_tolerance, rk_tolerance, constraint_step, i, -1.f * delta);
0572
0573       ridder.run(nvec1, nvec2, delta, p, i, differentiated_jacobian);
0574
0575       if (ridder.finished()) {
0576         num_iterations[i] = p;
0577         break;
0578       }
0579     }
0580
0581     // Row-major
0582     for (std::size_t j = 0u; j < 5u; j++) {
0583       convergence[i + j * 5] = ridder.complete[j];
0584     }
0585   }
0586
0587   return differentiated_jacobian;
0588 }
0589
0590 template <typename detector_t, typename detector_t::metadata::mask_id mid>
0591 bound_track_parameters<test_algebra> get_initial_parameter(
0592     const detector_t& det, const free_track_parameters<test_algebra>& vertex,
0593     const vector3& field, const scalar helix_tolerance) {
0594   // Helix from the vertex
0595   detail::helix<test_algebra> hlx(vertex, field);
0596
0597   const auto& departure_sf = det.surface(0u);
0598   const auto& trf_link = departure_sf.transform();
0599   const auto& departure_trf = det.transform_store().at(trf_link);
0600   const auto& mask_link = departure_sf.mask();
0601   const auto& departure_mask =
0602       det.mask_store().template get<mid>().at(mask_link.index());
0603
0604   using mask_t = types::get<typename detector_t::masks, mid>;
0605   helix_intersector<typename mask_t::shape, test_algebra> hlx_is{};
0606   hlx_is.run_rtsafe = false;
0607   hlx_is.convergence_tolerance = helix_tolerance;
0608   auto sfi = hlx_is(hlx, departure_sf, departure_mask, departure_trf, 0.f);
0609   EXPECT_TRUE(sfi.is_inside())
0610       << " Initial surface not found" << std::endl
0611       << " log10(Helix tolerance): " << math::log10(helix_tolerance)
0612       << " Phi: " << vector::phi(vertex.dir())
0613       << " Theta: " << vector::theta(vertex.dir())
0614       << " Mom [GeV/c]: " << vertex.p(ptc.charge()) << std::endl
0615       << sfi;
0616
0617   const auto path_length = sfi.path();
0618   // As we don't rotate or shift the initial surface anymore, the path_length
0619   // should be 0
0620   EXPECT_FLOAT_EQ(static_cast<float>(path_length), 0.f);
0621
0622   const auto pos = hlx(path_length);
0623   const auto dir = hlx.dir(path_length);
0624
0625   const free_track_parameters<test_algebra> free_par(pos, 0, dir, hlx.qop());
0626
0627   const auto bound_vec =
0628       tracking_surface{det, departure_sf}.free_to_bound_vector({}, free_par);
0629
0630   bound_track_parameters<test_algebra> ret;
0631   ret.set_surface_link(geometry::identifier{0u});
0632   ret.set_parameter_vector(bound_vec);
0633
0634   return ret;
0635 }
0636
0637 template <typename propagator_t, typename field_t>
0638 void evaluate_jacobian_difference(
0639     const std::size_t trk_count, const std::array<scalar, 3u>& euler_angles_I,
0640     const std::array<scalar, 3u>& euler_angles_F,
0641     const typename propagator_t::detector_type& det,
0642     const scalar detector_length,
0643     const bound_track_parameters<test_algebra>& track, const field_t& field,
0644     const scalar overstep_tolerance, const scalar path_tolerance,
0645     const scalar rk_tolerance, const scalar rk_tolerance_dis,
0646     const scalar constraint_step, const std::array<scalar, 5u>& hs,
0647     std::ofstream& file, scalar& ref_rel_diff, bool use_field_gradient,
0648     bool do_inspect, const bool use_precal_values = false,
0649     [[maybe_unused]] bound_covariance_type precal_diff_jacobi = {},
0650     [[maybe_unused]] std::array<unsigned int, 5u> precal_num_iterations = {},
0651     [[maybe_unused]] std::array<bool, 25u> precal_convergence = {}) {
0652   const auto phi0 = track.phi();
0653   const auto theta0 = track.theta();
0654   (void)phi0;
0655   (void)theta0;
0656
0657   auto bound_getter = evaluate_bound_param<propagator_t, field_t>(
0658       trk_count, detector_length, track, det, field, overstep_tolerance,
0659       path_tolerance, rk_tolerance, constraint_step, use_field_gradient, true,
0660       do_inspect);
0661
0662   const auto reference_param = bound_getter.m_param_departure;
0663   const auto final_param = bound_getter.m_param_destination;
0664
0665   // Sanity check
0666   ASSERT_EQ(reference_param.surface_link().index(), 0u)
0667       << " Initial surface not found " << std::endl
0668       << " log10(RK tolerance): " << math::log10(rk_tolerance)
0669       << " Path length [mm]: " << bound_getter.m_path_length
0670       << " Average step size [mm]: " << bound_getter.m_avg_step_size
0671       << " Phi: " << reference_param.phi()
0672       << " Theta: " << reference_param.theta()
0673       << " Mom [GeV/c]: " << reference_param.p(ptc.charge());
0674   ASSERT_EQ(final_param.surface_link().index(), 1u)
0675       << " Final surface not found " << std::endl
0676       << " log10(RK tolerance): " << math::log10(rk_tolerance)
0677       << " Path length [mm]: " << bound_getter.m_path_length
0678       << " Average step size [mm]: " << bound_getter.m_avg_step_size
0679       << " Phi: " << reference_param.phi()
0680       << " Theta: " << reference_param.theta()
0681       << " Mom [GeV/c]: " << reference_param.p(ptc.charge());
0682   ASSERT_TRUE(detector_length > 0.f);
0683   ASSERT_GE(bound_getter.m_path_length, 0.5f * detector_length);
0684   ASSERT_LE(bound_getter.m_path_length, 1.5f * detector_length);
0685   ASSERT_LE(bound_getter.m_path_length, max_detector_length + 200.f);
0686   ASSERT_GE(bound_getter.m_abs_path_length, bound_getter.m_path_length);
0687
0688   const auto reference_jacobian = bound_getter.m_jacobi;
0689
0690   file << trk_count << ",";
0691
0692   file << euler_angles_I[0u] << "," << euler_angles_I[1u] << ","
0693        << euler_angles_I[2u] << ",";
0694   file << euler_angles_F[0u] << "," << euler_angles_F[1u] << ","
0695        << euler_angles_F[2u] << ",";
0696
0697   file << reference_param.bound_local()[0] << ","
0698        << reference_param.bound_local()[1] << "," << reference_param.phi()
0699        << "," << reference_param.theta() << "," << reference_param.qop() << ",";
0700
0701   file << final_param.bound_local()[0] << "," << final_param.bound_local()[1]
0702        << "," << final_param.phi() << "," << final_param.theta() << ","
0703        << final_param.qop() << ",";
0704
0705   bound_covariance_type differentiated_jacobian{};
0706   std::array<unsigned int, 5u> num_iterations{};
0707   std::array<bool, 25u> convergence{};
0708
0709   if (use_precal_values) {
0710     differentiated_jacobian = precal_diff_jacobi;
0711     num_iterations = precal_num_iterations;
0712     convergence = precal_convergence;
0713   } else {
0714     differentiated_jacobian = directly_differentiate<propagator_t, field_t>(
0715         trk_count, reference_param, det, detector_length, field,
0716         overstep_tolerance, path_tolerance, rk_tolerance_dis, constraint_step,
0717         hs, num_iterations, convergence);
0718   }
0719
0720   bool total_convergence = (std::ranges::count(convergence, false) == 0);
0721
0722   // Ridders number of iterations
0723   for (unsigned int i = 0; i < 5u; i++) {
0724     file << num_iterations[i] << ",";
0725   }
0726
0727   // Convergence
0728   file << total_convergence << ",";
0729   for (unsigned int i = 0; i < 25u; i++) {
0730     file << convergence[i] << ",";
0731   }
0732
0733   // Reference track
0734   for (unsigned int i = 0; i < 5u; i++) {
0735     for (unsigned int j = 0; j < 5u; j++) {
0736       file << getter::element(reference_jacobian, i, j) << ",";
0737     }
0738   }
0739
0740   // Numerical evaluation
0741   for (unsigned int i = 0; i < 5u; i++) {
0742     for (unsigned int j = 0; j < 5u; j++) {
0743       file << getter::element(differentiated_jacobian, i, j) << ",";
0744     }
0745   }
0746
0747   // Difference between evaluation and direct jacobian
0748   for (unsigned int i = 0; i < 5u; i++) {
0749     for (unsigned int j = 0; j < 5u; j++) {
0750       const scalar ref_val = getter::element(reference_jacobian, i, j);
0751       const scalar num_val = getter::element(differentiated_jacobian, i, j);
0752       const scalar rel_diff = get_relative_difference(ref_val, num_val);
0753
0754       file << rel_diff << ",";
0755
0756       // We return dqopdqop for test
0757       if (i == 4 && j == 4) {
0758         ref_rel_diff = rel_diff;
0759       }
0760     }
0761   }
0762
0763   // Path length
0764   file << bound_getter.m_path_length << ",";
0765
0766   // Absolute path length
0767   file << bound_getter.m_abs_path_length << ",";
0768
0769   // The number of steps made
0770   file << bound_getter.step_count << ",";
0771
0772   // Average step size
0773   file << bound_getter.m_avg_step_size << ",";
0774
0775   // Log10(RK tolerance)
0776   file << math::log10(rk_tolerance) << ",";
0777
0778   // Log10(on surface tolerance)
0779   file << math::log10(path_tolerance) << ",";
0780
0781   // Overstep tolerance
0782   file << overstep_tolerance << ",";
0783
0784   // Detector length
0785   file << detector_length;
0786
0787   file << std::endl;
0788 }
0789
0790 template <typename propagator_t, typename field_t>
0791 void evaluate_covariance_transport(
0792     const std::size_t trk_count, const std::array<scalar, 3u>& euler_angles_I,
0793     const std::array<scalar, 3u>& euler_angles_F,
0794     const typename propagator_t::detector_type& det,
0795     const scalar detector_length,
0796     const bound_track_parameters<test_algebra>& track, const field_t& field,
0797     const scalar overstep_tolerance, const scalar path_tolerance,
0798     const scalar rk_tolerance, const scalar rk_tolerance_dis,
0799     const scalar constraint_step, std::ofstream& file,
0800     bool use_field_gradient) {
0801   // Copy track
0802   auto track_copy = track;
0803
0804   // Make initial covariance
0805   const bound_covariance_type ini_cov =
0806       get_random_initial_covariance(track_copy.qop());
0807
0808   track_copy.set_covariance(ini_cov);
0809
0810   auto bound_getter = evaluate_bound_param<propagator_t, field_t>(
0811       trk_count, detector_length, track_copy, det, field, overstep_tolerance,
0812       path_tolerance, rk_tolerance, constraint_step, use_field_gradient, true,
0813       false);
0814
0815   const auto reference_param = bound_getter.m_param_departure;
0816   const auto final_param = bound_getter.m_param_destination;
0817   const auto fin_cov = final_param.covariance();
0818
0819   // Sanity check
0820   ASSERT_EQ(reference_param.surface_link().index(), 0u)
0821       << " Initial surface not found " << std::endl
0822       << " log10(RK tolerance): " << math::log10(rk_tolerance)
0823       << " Path length [mm]: " << bound_getter.m_path_length
0824       << " Average step size [mm]: " << bound_getter.m_avg_step_size
0825       << " Phi: " << reference_param.phi()
0826       << " Theta: " << reference_param.theta()
0827       << " Mom [GeV/c]: " << reference_param.p(ptc.charge());
0828   ASSERT_EQ(final_param.surface_link().index(), 1u)
0829       << " Final surface not found " << std::endl
0830       << " log10(RK tolerance): " << math::log10(rk_tolerance)
0831       << " Path length [mm]: " << bound_getter.m_path_length
0832       << " Average step size [mm]: " << bound_getter.m_avg_step_size
0833       << " Phi: " << reference_param.phi()
0834       << " Theta: " << reference_param.theta()
0835       << " Mom [GeV/c]: " << reference_param.p(ptc.charge());
0836   ASSERT_TRUE(detector_length > 0.f);
0837   ASSERT_GE(bound_getter.m_path_length, 0.5f * detector_length);
0838   ASSERT_LE(bound_getter.m_path_length, 1.5f * detector_length);
0839   ASSERT_LE(bound_getter.m_path_length, max_detector_length + 200.f);
0840   ASSERT_GE(bound_getter.m_abs_path_length, bound_getter.m_path_length);
0841
0842   // Get smeared initial bound vector
0843   const bound_param_vector_type smeared_ini_vec =
0844       get_smeared_bound_vector(ini_cov, reference_param);
0845
0846   // Make smeared bound track parameter
0847   auto smeared_track = track_copy;
0848   smeared_track.set_parameter_vector(smeared_ini_vec);
0849
0850   auto smeared_bound_getter = evaluate_bound_param<propagator_t, field_t>(
0851       trk_count, detector_length, smeared_track, det, field, overstep_tolerance,
0852       path_tolerance, rk_tolerance_dis, constraint_step, use_field_gradient,
0853       false, false);
0854
0855   // Get smeared final bound vector
0856   bound_param_vector_type smeared_fin_vec =
0857       smeared_bound_getter.m_param_destination;
0858
0859   // phi needs to be wrapped w.r.t. phi of the reference vector
0860   wrap_angles(final_param, smeared_fin_vec);
0861
0862   // Get pull values
0863   std::array<scalar, 5u> pulls{};
0864
0865   bound_vector<test_algebra> diff{};
0866   for (unsigned int i = 0u; i < 5u; i++) {
0867     getter::element(diff, i, 0u) = smeared_fin_vec[i] - final_param[i];
0868     pulls[i] = getter::element(diff, i, 0u) /
0869                math::sqrt(getter::element(fin_cov, i, i));
0870   }
0871
0872   // Get Chi2
0873   const matrix_type<1u, 1u> chi2 =
0874       matrix::transpose(diff) * matrix::inverse(fin_cov) * diff;
0875   const scalar chi2_val = getter::element(chi2, 0u, 0u);
0876
0877   file << trk_count << ",";
0878
0879   file << euler_angles_I[0u] << "," << euler_angles_I[1u] << ","
0880        << euler_angles_I[2u] << ",";
0881   file << euler_angles_F[0u] << "," << euler_angles_F[1u] << ","
0882        << euler_angles_F[2u] << ",";
0883
0884   // File writing
0885   file << reference_param[e_bound_loc0] << "," << reference_param[e_bound_loc1]
0886        << "," << reference_param[e_bound_phi] << ","
0887        << reference_param[e_bound_theta] << ","
0888        << reference_param[e_bound_qoverp] << ",";
0889
0890   for (unsigned int i = 0; i < 5u; i++) {
0891     for (unsigned int j = 0; j < 5u; j++) {
0892       file << getter::element(ini_cov, i, j) << ",";
0893     }
0894   }
0895
0896   file << final_param[e_bound_loc0] << "," << final_param[e_bound_loc1] << ","
0897        << final_param[e_bound_phi] << "," << final_param[e_bound_theta] << ","
0898        << final_param[e_bound_qoverp] << ",";
0899
0900   for (unsigned int i = 0; i < 5u; i++) {
0901     for (unsigned int j = 0; j < 5u; j++) {
0902       file << getter::element(fin_cov, i, j) << ",";
0903     }
0904   }
0905
0906   file << smeared_ini_vec[e_bound_loc0] << "," << smeared_ini_vec[e_bound_loc1]
0907        << "," << smeared_ini_vec[e_bound_phi] << ","
0908        << smeared_ini_vec[e_bound_theta] << ","
0909        << smeared_ini_vec[e_bound_qoverp] << ",";
0910
0911   file << smeared_fin_vec[e_bound_loc0] << "," << smeared_fin_vec[e_bound_loc1]
0912        << "," << smeared_fin_vec[e_bound_phi] << ","
0913        << smeared_fin_vec[e_bound_theta] << ","
0914        << smeared_fin_vec[e_bound_qoverp] << ",";
0915
0916   file << pulls[0] << "," << pulls[1] << "," << pulls[2] << "," << pulls[3]
0917        << "," << pulls[4] << ",";
0918
0919   file << chi2_val << ",";
0920
0921   // Path length
0922   file << bound_getter.m_path_length << ",";
0923
0924   // Absolute path length
0925   file << bound_getter.m_abs_path_length << ",";
0926
0927   // The number of steps made
0928   file << bound_getter.step_count << ",";
0929
0930   // Average step size
0931   file << bound_getter.m_avg_step_size << ",";
0932
0933   // Log10(RK tolerance)
0934   file << math::log10(rk_tolerance) << ",";
0935
0936   // Log10(on surface tolerance)
0937   file << math::log10(path_tolerance) << ",";
0938
0939   // Overstep tolerance
0940   file << overstep_tolerance << ",";
0941
0942   // Detector length
0943   file << detector_length;
0944
0945   file << std::endl;
0946 }
0947
0948 template <typename detector_t, typename detector_t::metadata::mask_id mask_id>
0949 bound_param_vector_type get_displaced_bound_vector_helix(
0950     const bound_track_parameters<test_algebra>& track, const vector3& field,
0951     unsigned int target_index, scalar displacement, const detector_t& det,
0952     const scalar helix_tolerance) {
0953   const auto& departure_sf = det.surface(0u);
0954
0955   const auto& destination_sf = det.surface(1u);
0956   const auto& trf_link = destination_sf.transform();
0957   const auto& destination_trf = det.transform_store().at(trf_link);
0958   const auto& mask_link = destination_sf.mask();
0959   const auto& destination_mask =
0960       det.mask_store().template get<mask_id>().at(mask_link.index());
0961
0962   bound_param_vector_type dvec = track;
0963   dvec[target_index] += displacement;
0964   const auto free_vec =
0965       tracking_surface{det, departure_sf}.bound_to_free_vector({}, dvec);
0966   detail::helix<test_algebra> hlx(free_vec, field);
0967
0968   using mask_t = types::get<typename detector_t::masks, mask_id>;
0969   helix_intersector<typename mask_t::shape, test_algebra> hlx_is{};
0970   hlx_is.run_rtsafe = false;
0971   hlx_is.convergence_tolerance = helix_tolerance;
0972   auto sfi =
0973       hlx_is(hlx, destination_sf, destination_mask, destination_trf, 0.f);
0974   const auto path_length = sfi.path();
0975   const auto pos = hlx(path_length);
0976   const auto dir = hlx.dir(path_length);
0977
0978   const free_track_parameters<test_algebra> new_free_par(pos, 0, dir,
0979                                                          hlx.qop());
0980   auto new_bound_vec =
0981       tracking_surface{det, destination_sf}.free_to_bound_vector({},
0982                                                                  new_free_par);
0983
0984   // phi needs to be wrapped w.r.t. phi of the reference vector
0985   wrap_angles(dvec, new_bound_vec);
0986
0987   return new_bound_vec;
0988 }
0989
0990 template <typename detector_t, typename detector_t::metadata::mask_id mask_id>
0991 void evaluate_jacobian_difference_helix(
0992     const std::size_t trk_count, const std::array<scalar, 3u> euler_angles_I,
0993     const std::array<scalar, 3u> euler_angles_F, const detector_t& det,
0994     const scalar detector_length,
0995     const bound_track_parameters<test_algebra>& track, const vector3& field,
0996     const std::array<scalar, 5u> hs, std::ofstream& file,
0997     const scalar helix_tolerance) {
0998   const auto phi0 = track.phi();
0999   const auto theta0 = track.theta();
1000   (void)phi0;
1001   (void)theta0;
1002
1003   // Get bound to free Jacobi
1004   const auto& departure_sf = det.surface(0u);
1005   const auto bound_to_free_jacobi =
1006       tracking_surface{det, departure_sf}.bound_to_free_jacobian({}, track);
1007
1008   // Get fre vector
1009   const auto free_vec =
1010       tracking_surface{det, departure_sf}.bound_to_free_vector({}, track);
1011
1012   // Helix from the departure surface
1013   detail::helix<test_algebra> hlx(free_vec, field);
1014
1015   const auto& destination_sf = det.surface(1u);
1016   const auto& trf_link = destination_sf.transform();
1017   const auto& destination_trf = det.transform_store().at(trf_link);
1018   const auto& mask_link = destination_sf.mask();
1019   const auto& destination_mask =
1020       det.mask_store().template get<mask_id>().at(mask_link.index());
1021
1022   using mask_t = types::get<typename detector_t::masks, mask_id>;
1023   helix_intersector<typename mask_t::shape, test_algebra> hlx_is{};
1024   hlx_is.run_rtsafe = false;
1025   hlx_is.convergence_tolerance = helix_tolerance;
1026
1027   auto sfi =
1028       hlx_is(hlx, destination_sf, destination_mask, destination_trf, 0.f);
1029
1030   EXPECT_TRUE(sfi.is_inside())
1031       << " Final surface not found" << std::endl
1032       << " log10(Helix tolerance): " << math::log10(helix_tolerance)
1033       << " Phi: " << track.phi() << " Theta: " << track.theta()
1034       << " Mom [GeV/c]: " << track.p(ptc.charge());
1035
1036   const auto path_length = sfi.path();
1037
1038   // Get transport Jacobi
1039   const auto transport_jacobi = hlx.jacobian(path_length);
1040
1041   const auto pos = hlx(path_length);
1042   const auto dir = hlx.dir(path_length);
1043   const auto qop = hlx.qop();
1044
1045   // Get correction term
1046   const auto correction_term =
1047       matrix::identity<free_matrix<test_algebra>>() +
1048       tracking_surface{det, destination_sf}.path_correction(
1049           {}, pos, dir, qop * vector::cross(dir, field), 0.f);
1050
1051   const free_track_parameters<test_algebra> free_par(pos, 0.f, dir, qop);
1052
1053   // Get free to bound Jacobi
1054   const auto free_to_bound_jacobi =
1055       tracking_surface{det, destination_sf}.free_to_bound_jacobian({},
1056                                                                    free_par);
1057
1058   // Get full Jacobi
1059   const auto reference_jacobian = free_to_bound_jacobi * correction_term *
1060                                   transport_jacobi * bound_to_free_jacobi;
1061
1062   // Get bound vector
1063   const auto bound_vec =
1064       tracking_surface{det, destination_sf}.free_to_bound_vector({}, free_par);
1065
1066   /******************************
1067    *  Numerical differentiation
1068    * ****************************/
1069
1070   bound_covariance_type differentiated_jacobian{};
1071   std::array<unsigned int, 5u> num_iterations{};
1072   std::array<bool, 25u> convergence{};
1073
1074   for (unsigned int i = 0u; i < 5u; i++) {
1075     scalar delta = hs[i];
1076
1077     preprocess_delta(i, delta, track);
1078
1079     const auto vec1 = get_displaced_bound_vector_helix<detector_t, mask_id>(
1080         track, field, i, 1.f * delta, det, helix_tolerance);
1081     const auto vec2 = get_displaced_bound_vector_helix<detector_t, mask_id>(
1082         track, field, i, -1.f * delta, det, helix_tolerance);
1083
1084     ridders_derivative ridder;
1085     ridder.initialize(vec1, vec2, delta);
1086
1087     for (unsigned int p = 1u; p < Nt; p++) {
1088       delta /= con[i];
1089
1090       const auto nvec1 = get_displaced_bound_vector_helix<detector_t, mask_id>(
1091           track, field, i, 1.f * delta, det, helix_tolerance);
1092       const auto nvec2 = get_displaced_bound_vector_helix<detector_t, mask_id>(
1093           track, field, i, -1.f * delta, det, helix_tolerance);
1094
1095       ridder.run(nvec1, nvec2, delta, p, i, differentiated_jacobian);
1096
1097       if (ridder.finished()) {
1098         num_iterations[i] = p;
1099         break;
1100       }
1101     }
1102
1103     for (std::size_t j = 0u; j < 5u; j++) {
1104       convergence[i + j * 5] = ridder.complete[j];
1105     }
1106   }
1107
1108   bool total_convergence = (std::ranges::count(convergence, false) == 0);
1109
1110   file << trk_count << ",";
1111
1112   file << euler_angles_I[0u] << "," << euler_angles_I[1u] << ","
1113        << euler_angles_I[2u] << ",";
1114   file << euler_angles_F[0u] << "," << euler_angles_F[1u] << ","
1115        << euler_angles_F[2u] << ",";
1116
1117   file << track.bound_local()[0] << "," << track.bound_local()[1] << ","
1118        << track.phi() << "," << track.theta() << "," << track.qop() << ",";
1119
1120   file << bound_vec[e_bound_loc0] << "," << bound_vec[e_bound_loc1] << ","
1121        << bound_vec.phi() << "," << bound_vec.theta() << "," << bound_vec.qop()
1122        << ",";
1123
1124   // Ridders number of iterations
1125   for (unsigned int i = 0; i < 5u; i++) {
1126     file << num_iterations[i] << ",";
1127   }
1128
1129   // Convergence
1130   file << total_convergence << ",";
1131   for (unsigned int i = 0; i < 25u; i++) {
1132     file << convergence[i] << ",";
1133   }
1134
1135   // Reference track
1136   for (unsigned int i = 0; i < 5u; i++) {
1137     for (unsigned int j = 0; j < 5u; j++) {
1138       file << getter::element(reference_jacobian, i, j) << ",";
1139     }
1140   }
1141
1142   // Numerical evaluation
1143   for (unsigned int i = 0; i < 5u; i++) {
1144     for (unsigned int j = 0; j < 5u; j++) {
1145       file << getter::element(differentiated_jacobian, i, j) << ",";
1146     }
1147   }
1148
1149   // Difference between evaluation and direct jacobian
1150   for (unsigned int i = 0; i < 5u; i++) {
1151     for (unsigned int j = 0; j < 5u; j++) {
1152       const scalar ref_val = getter::element(reference_jacobian, i, j);
1153       const scalar num_val = getter::element(differentiated_jacobian, i, j);
1154       const scalar rel_diff = get_relative_difference(ref_val, num_val);
1155       file << rel_diff << ",";
1156     }
1157   }
1158
1159   // Path length
1160   file << path_length << ",";
1161
1162   // Absolute Path length (Not supported for helix)
1163   file << 0 << ",";
1164
1165   // Average step size (Doesn't exist for helix intersection)
1166   file << 0 << ",";
1167
1168   // The number of steps made (Doesn't exist for helix intersection)
1169   file << 0 << ",";
1170
1171   // Log10(RK tolerance) (Doesn't exist for helix intersection)
1172   file << 0 << ",";
1173
1174   // Log10(helix intersection tolerance)
1175   file << math::log10(helix_tolerance) << ",";
1176
1177   // Overstep tolerance (Doesn't exist for helix intersection)
1178   file << 0 << ",";
1179
1180   // Detector length
1181   file << detector_length;
1182
1183   file << std::endl;
1184 }
1185
1186 void setup_csv_header_jacobian(std::ofstream& file) {
1187   file << std::fixed << std::showpoint;
1188   file << std::setprecision(32);
1189
1190   // Track ID
1191   file << "track_ID,";
1192
1193   // Euler angles
1194   file << "alpha_I,beta_I,gamma_I,";
1195   file << "alpha_F,beta_F,gamma_F,";
1196
1197   // Initial Parameter at the departure surface
1198   file << "l0_I,l1_I,phi_I,theta_I,qop_I,";
1199
1200   // Final Parameter at the destination surface
1201   file << "l0_F,l1_F,phi_F,theta_F,qop_F,";
1202
1203   // Number of iterations to complete the numerical differentiation
1204   file << "num_iterations_l0,"
1205        << "num_iterations_l1,"
1206        << "num_iterations_phi,"
1207        << "num_iterations_theta,"
1208        << "num_iterations_qop,";
1209
1210   // Convergence
1211   file << "total_convergence,";
1212   file << "dl0dl0_C,dl0dl1_C,dl0dphi_C,dl0dtheta_C,dl0dqop_C,";
1213   file << "dl1dl0_C,dl1dl1_C,dl1dphi_C,dl1dtheta_C,dl1dqop_C,";
1214   file << "dphidl0_C,dphidl1_C,dphidphi_C,dphidtheta_C,dphidqop_C,";
1215   file << "dthetadl0_C,dthetadl1_C,dthetadphi_C,dthetadtheta_C,"
1216           "dthetadqop_C,";
1217   file << "dqopdl0_C,dqopdl1_C,dqopdphi_C,dqopdtheta_C,dqopdqop_C,";
1218
1219   // Evaluation
1220   file << "dl0dl0_E,dl0dl1_E,dl0dphi_E,dl0dtheta_E,dl0dqop_E,";
1221   file << "dl1dl0_E,dl1dl1_E,dl1dphi_E,dl1dtheta_E,dl1dqop_E,";
1222   file << "dphidl0_E,dphidl1_E,dphidphi_E,dphidtheta_E,dphidqop_E,";
1223   file << "dthetadl0_E,dthetadl1_E,dthetadphi_E,dthetadtheta_E,"
1224           "dthetadqop_E,";
1225   file << "dqopdl0_E,dqopdl1_E,dqopdphi_E,dqopdtheta_E,dqopdqop_E,";
1226
1227   // Numerical Differentiation
1228   file << "dl0dl0_D,dl0dl1_D,dl0dphi_D,dl0dtheta_D,dl0dqop_D,";
1229   file << "dl1dl0_D,dl1dl1_D,dl1dphi_D,dl1dtheta_D,dl1dqop_D,";
1230   file << "dphidl0_D,dphidl1_D,dphidphi_D,dphidtheta_D,dphidqop_D,";
1231   file << "dthetadl0_D,dthetadl1_D,dthetadphi_D,dthetadtheta_D,"
1232           "dthetadqop_D,";
1233   file << "dqopdl0_D,dqopdl1_D,dqopdphi_D,dqopdtheta_D,dqopdqop_D,";
1234
1235   // Relative Difference between Evaluation and Numerical Differentiation
1236   file << "dl0dl0_R,dl0dl1_R,dl0dphi_R,dl0dtheta_R,dl0dqop_R,";
1237   file << "dl1dl0_R,dl1dl1_R,dl1dphi_R,dl1dtheta_R,dl1dqop_R,";
1238   file << "dphidl0_R,dphidl1_R,dphidphi_R,dphidtheta_R,dphidqop_R,";
1239   file << "dthetadl0_R,dthetadl1_R,dthetadphi_R,dthetadtheta_R,"
1240           "dthetadqop_R,";
1241   file << "dqopdl0_R,dqopdl1_R,dqopdphi_R,dqopdtheta_R,dqopdqop_R,";
1242
1243   // Path length [mm]
1244   file << "path_length,";
1245
1246   // Absolute path length [mm]
1247   file << "abs_path_length,";
1248
1249   // The number of steps made
1250   file << "n_steps,";
1251
1252   // Average step size [mm]
1253   file << "average_step_size,";
1254
1255   // RK Tolerances [mm]
1256   file << "log10_rk_tolerance,";
1257
1258   // Intersection tolerance [mm]
1259   file << "log10_intersection_tolerance,";
1260
1261   // Overstep tolerance [mm]
1262   file << "overstep_tolerance,";
1263
1264   // Detector length [mm]
1265   file << "detector_length";
1266
1267   file << std::endl;
1268 }
1269
1270 void setup_csv_header_covariance(std::ofstream& file) {
1271   file << std::fixed << std::showpoint;
1272   file << std::setprecision(32);
1273
1274   // Track ID
1275   file << "track_ID,";
1276
1277   // Euler angles
1278   file << "alpha_I,beta_I,gamma_I,";
1279   file << "alpha_F,beta_F,gamma_F,";
1280
1281   // Initial parameters (vector + covariance) of reference track
1282   file << "l0_I,l1_I,phi_I,theta_I,qop_I,";
1283   file << "l0l0_I,l0l1_I,l0phi_I,l0theta_I,l0qop_I,";
1284   file << "l1l0_I,l1l1_I,l1phi_I,l1theta_I,l1qop_I,";
1285   file << "phil0_I,phil1_I,phiphi_I,phitheta_I,phiqop_I,";
1286   file << "thetal0_I,thetal1_I,thetaphi_I,thetatheta_I,thetaqop_I,";
1287   file << "qopl0_I,qopl1_I,qopphi_I,qoptheta_I,qopqop_I,";
1288
1289   // Final parameters (vector + covariance) of reference track
1290   file << "l0_F,l1_F,phi_F,theta_F,qop_F,";
1291   file << "l0l0_F,l0l1_F,l0phi_F,l0theta_F,l0qop_F,";
1292   file << "l1l0_F,l1l1_F,l1phi_F,l1theta_F,l1qop_F,";
1293   file << "phil0_F,phil1_F,phiphi_F,phitheta_F,phiqop_F,";
1294   file << "thetal0_F,thetal1_F,thetaphi_F,thetatheta_F,thetaqop_F,";
1295   file << "qopl0_F,qopl1_F,qopphi_F,qoptheta_F,qopqop_F,";
1296
1297   // Initial parameter (vector only) of smeared track
1298   file << "l0_IS,l1_IS,phi_IS,theta_IS,qop_IS,";
1299
1300   // Final parameter (vector only) of smeared track
1301   file << "l0_FS,l1_FS,phi_FS,theta_FS,qop_FS,";
1302
1303   // Pull values
1304   file << "pull_l0,pull_l1,pull_phi,pull_theta,pull_qop,";
1305
1306   // Chi2
1307   file << "chi2,";
1308
1309   // Path length [mm]
1310   file << "path_length,";
1311
1312   // Absolute path length [mm]
1313   file << "abs_path_length,";
1314
1315   // The number of steps made
1316   file << "n_steps,";
1317
1318   // Average step size [mm]
1319   file << "average_step_size,";
1320
1321   // Tolerances [mm]
1322   file << "log10_rk_tolerance,";
1323
1324   // Intersection tolerance [mm]
1325   file << "log10_intersection_tolerance,";
1326
1327   // Overstep tolerance [mm]
1328   file << "overstep_tolerance,";
1329
1330   // Detector length [mm]
1331   file << "detector_length";
1332
1333   file << std::endl;
1334 }
1335
1336 int main(int argc, char** argv) {
1337   // Options parsing
1338   po::options_description desc("\ndetray jacobian validation options");
1339   desc.add_options()("help", "produce help message");
1340   desc.add_options()("output-directory",
1341                      po::value<std::string>()->default_value(""),
1342                      "Output directory");
1343   desc.add_options()("n-tracks", po::value<std::size_t>()->default_value(100u),
1344                      "Number of tracks for generator");
1345   desc.add_options()("n-skips", po::value<std::size_t>()->default_value(0u),
1346                      "Number of skipped indices");
1347   desc.add_options()("skip-rect", po::value<bool>()->default_value(false),
1348                      "Skip rectangular telescope");
1349   desc.add_options()("skip-wire", po::value<bool>()->default_value(false),
1350                      "Skip wire telescope");
1351   desc.add_options()("log10-rk-tolerance-jac-mm",
1352                      po::value<scalar>()->default_value(-4.f),
1353                      "Set log10(rk_tolerance_jac_in_mm)");
1354   desc.add_options()("log10-rk-tolerance-dis-mm",
1355                      po::value<scalar>()->default_value(-6.f),
1356                      "Set log10(rk_tolerance_dis_in_mm)");
1357   desc.add_options()("log10-rk-tolerance-cov-mm",
1358                      po::value<scalar>()->default_value(-4.f),
1359                      "Set log10(rk_tolerance_cov_in_mm)");
1360   desc.add_options()("log10-helix-tolerance-mm",
1361                      po::value<scalar>()->default_value(-3.f),
1362                      "Set log10(helix_tolerance_in_mm)");
1363   desc.add_options()("overstep-tolerance-mm",
1364                      po::value<scalar>()->default_value(-1000.f),
1365                      "Set the overstep tolerance in mm unit");
1366   desc.add_options()("log10-on-surface-tolerance-mm",
1367                      po::value<scalar>()->default_value(-3.f),
1368                      "Set log10(path_tolerance_in_mm)");
1369   desc.add_options()("rk-tolerance-iterate-mode",
1370                      po::value<bool>()->default_value(true),
1371                      "Iterate over the rk tolerances");
1372   desc.add_options()("log10-min-rk-tolerance-mm",
1373                      po::value<scalar>()->default_value(-6.f),
1374                      "Set log10(min_rk_tolerance_in_mm)");
1375   desc.add_options()("log10-max-rk-tolerance-mm",
1376                      po::value<scalar>()->default_value(2.f),
1377                      "Set log10(max_rk_tolerance_in_mm)");
1378   desc.add_options()("mc-seed", po::value<std::size_t>()->default_value(0u),
1379                      "Monte-Carlo seed");
1380   desc.add_options()("verbose-level", po::value<int>()->default_value(1),
1381                      "Verbose level");
1382
1383   po::variables_map vm;
1384   po::store(parse_command_line(argc, argv, desc,
1385                                po::command_line_style::unix_style ^
1386                                    po::command_line_style::allow_short),
1387             vm);
1388   po::notify(vm);
1389
1390   // Help message
1391   if (vm.count("help") != 0u) {
1392     std::clog << desc << std::endl;
1393     return EXIT_FAILURE;
1394   }
1395
1396   const std::string output_directory = vm["output-directory"].as<std::string>();
1397   std::size_t n_tracks = vm["n-tracks"].as<std::size_t>();
1398   std::size_t n_skips = vm["n-skips"].as<std::size_t>();
1399   const bool skip_rect = vm["skip-rect"].as<bool>();
1400   const bool skip_wire = vm["skip-wire"].as<bool>();
1401   const scalar rk_power_jac = vm["log10-rk-tolerance-jac-mm"].as<scalar>();
1402   const scalar rk_tol_jac = std::pow(10.f, rk_power_jac) * unit<scalar>::mm;
1403   const scalar rk_power_dis = vm["log10-rk-tolerance-dis-mm"].as<scalar>();
1404   const scalar rk_tol_dis = std::pow(10.f, rk_power_dis) * unit<scalar>::mm;
1405   const scalar rk_power_cov = vm["log10-rk-tolerance-cov-mm"].as<scalar>();
1406   const scalar rk_tol_cov = std::pow(10.f, rk_power_cov) * unit<scalar>::mm;
1407   const scalar helix_power = vm["log10-helix-tolerance-mm"].as<scalar>();
1408   const scalar helix_tol = std::pow(10.f, helix_power) * unit<scalar>::mm;
1409   const scalar on_surface_power =
1410       vm["log10-on-surface-tolerance-mm"].as<scalar>() * unit<scalar>::mm;
1411   const scalar on_surface_tol = math::pow(10.f, on_surface_power);
1412   const bool rk_tolerance_iterate_mode =
1413       vm["rk-tolerance-iterate-mode"].as<bool>();
1414   const scalar log10_min_rk_tolerance =
1415       vm["log10-min-rk-tolerance-mm"].as<scalar>() * unit<scalar>::mm;
1416   const scalar log10_max_rk_tolerance =
1417       vm["log10-max-rk-tolerance-mm"].as<scalar>() * unit<scalar>::mm;
1418   const std::size_t mc_seed = vm["mc-seed"].as<std::size_t>();
1419   const int verbose_lvl = vm["verbose-level"].as<int>();
1420
1421   std::vector<scalar> log10_tols;
1422   scalar r = log10_min_rk_tolerance;
1423   while (r <= log10_max_rk_tolerance + 1e-3f) {
1424     log10_tols.push_back(r);
1425     r = r + 2.f;
1426   }
1427
1428   // Set seed for random generator
1429   mt1.seed(mc_seed);
1430
1431   // Volume material
1432   const material<scalar> volume_mat = detray::cesium_iodide_with_ded<scalar>();
1433
1434   std::string path;
1435   // Create output directory
1436   if (output_directory.empty()) {
1437     path = "";
1438   } else {
1439     std::filesystem::create_directories(output_directory);
1440     path = output_directory + "/";
1441   }
1442
1443   // Output Csv file
1444   std::ofstream helix_rect_file;
1445   std::ofstream const_rect_file;
1446   std::ofstream inhom_rect_file;
1447   std::ofstream inhom_rect_material_file;
1448   std::ofstream helix_wire_file;
1449   std::ofstream const_wire_file;
1450   std::ofstream inhom_wire_file;
1451   std::ofstream inhom_wire_material_file;
1452   helix_rect_file.open(path + "helix_rect.csv");
1453   const_rect_file.open(path + "const_rect.csv");
1454   inhom_rect_file.open(path + "inhom_rect.csv");
1455   inhom_rect_material_file.open(path + "inhom_rect_material.csv");
1456   helix_wire_file.open(path + "helix_wire.csv");
1457   const_wire_file.open(path + "const_wire.csv");
1458   inhom_wire_file.open(path + "inhom_wire.csv");
1459   inhom_wire_material_file.open(path + "inhom_wire_material.csv");
1460
1461   setup_csv_header_jacobian(helix_rect_file);
1462   setup_csv_header_jacobian(const_rect_file);
1463   setup_csv_header_jacobian(inhom_rect_file);
1464   setup_csv_header_jacobian(inhom_rect_material_file);
1465   setup_csv_header_jacobian(helix_wire_file);
1466   setup_csv_header_jacobian(const_wire_file);
1467   setup_csv_header_jacobian(inhom_wire_file);
1468   setup_csv_header_jacobian(inhom_wire_material_file);
1469
1470   std::ofstream rect_cov_transport_file;
1471   rect_cov_transport_file.open(path + "rect_cov_transport.csv");
1472   setup_csv_header_covariance(rect_cov_transport_file);
1473
1474   std::ofstream wire_cov_transport_file;
1475   wire_cov_transport_file.open(path + "wire_cov_transport.csv");
1476   setup_csv_header_covariance(wire_cov_transport_file);
1477
1478   // Output Csv file (RK tolerance iteration mode)
1479   std::vector<std::ofstream> rect_files(log10_tols.size());
1480   std::vector<std::ofstream> wire_files(log10_tols.size());
1481
1482   if (rk_tolerance_iterate_mode) {
1483     for (std::size_t i = 0u; i < log10_tols.size(); i++) {
1484       const std::string rect_name =
1485           "inhom_rect_material_" +
1486           std::to_string(static_cast<int>(log10_tols[i])) + ".csv";
1487       const std::string wire_name =
1488           "inhom_wire_material_" +
1489           std::to_string(static_cast<int>(log10_tols[i])) + ".csv";
1490       rect_files[i].open(path + rect_name);
1491       wire_files[i].open(path + wire_name);
1492
1493       setup_csv_header_jacobian(rect_files[i]);
1494       setup_csv_header_jacobian(wire_files[i]);
1495     }
1496   }
1497
1498   // Memory resource
1499   vecmem::host_memory_resource host_mr;
1500
1501   // Filter out the google test flags
1502   ::testing::InitGoogleTest(&argc, argv);
1503
1504   // Detector types
1505   using rectangle_telescope =
1506       detector<telescope_metadata<test_algebra, rect_type>>;
1507   using wire_telescope = detector<telescope_metadata<test_algebra, wire_type>>;
1508   using track_type = free_track_parameters<test_algebra>;
1509
1510   // Constant magnetic field type
1511   using const_bfield_t = bfield::const_field_t<scalar>;
1512
1513   // Magnetic field map using nearest neighbor interpolation
1514   using inhom_bfield_t = bfield::inhom_field_t<scalar>;
1515
1516   const const_bfield_t const_bfield = create_const_field<scalar>(B_z);
1517   const inhom_bfield_t inhom_bfield = create_inhom_field<scalar>();
1518
1519   // Actor chain type
1520   using actor_chain_t = actor_chain<
1521       actor::parameter_updater<test_algebra, bound_getter<test_algebra>>>;
1522
1523   // Iterate over reference (pilot) tracks for a rectangular telescope
1524   // geometry and Jacobian calculation
1525   using uniform_gen_t =
1526       detail::random_numbers<scalar, std::uniform_real_distribution<scalar>>;
1527   using trk_generator_t = random_track_generator<track_type, uniform_gen_t>;
1528   trk_generator_t::configuration trk_gen_cfg{};
1529   trk_gen_cfg.seed(mc_seed);
1530   trk_gen_cfg.n_tracks(n_tracks + n_skips);
1531   trk_gen_cfg.phi_range(-constant<scalar>::pi, constant<scalar>::pi);
1532   trk_gen_cfg.theta_range(0.f, constant<scalar>::pi);
1533   trk_gen_cfg.mom_range(min_mom, max_mom);
1534   trk_gen_cfg.origin(0.f, 0.f, 0.f);
1535   trk_gen_cfg.origin_stddev(0.f, 0.f, 0.f);
1536
1537   // Vectors for dqopdqop relative difference
1538   std::vector<std::vector<scalar>> dqopdqop_rel_diffs_rect(log10_tols.size());
1539   std::vector<std::vector<scalar>> dqopdqop_rel_diffs_wire(log10_tols.size());
1540
1541   bool do_inspect = false;
1542   if (verbose_lvl >= 4) {
1543     do_inspect = true;
1544   }
1545
1546   // Navigator types
1547   using rect_navigator_t = caching_navigator<rectangle_telescope>;
1548   using wire_navigator_t = caching_navigator<wire_telescope>;
1549
1550   // Stepper types
1551   using const_field_stepper_t =
1552       rk_stepper<const_bfield_t::view_t, test_algebra, constrained_step<scalar>,
1553                  stepper_default_policy<scalar>>;
1554   using inhom_field_stepper_t =
1555       rk_stepper<inhom_bfield_t::view_t, test_algebra, constrained_step<scalar>,
1556                  stepper_default_policy<scalar>>;
1557
1558   // Make four propagators for each case
1559   using const_field_rect_propagator_t =
1560       propagator<const_field_stepper_t, rect_navigator_t, actor_chain_t>;
1561   using inhom_field_rect_propagator_t =
1562       propagator<inhom_field_stepper_t, rect_navigator_t, actor_chain_t>;
1563
1564   using const_field_wire_propagator_t =
1565       propagator<const_field_stepper_t, wire_navigator_t, actor_chain_t>;
1566   using inhom_field_wire_propagator_t =
1567       propagator<inhom_field_stepper_t, wire_navigator_t, actor_chain_t>;
1568
1569   std::size_t track_count = 0u;
1570
1571   for (const auto track : trk_generator_t{trk_gen_cfg}) {
1572     mt2.seed(track_count);
1573
1574     // Pilot track
1575     detail::helix<test_algebra> helix_bz(track, B_z);
1576
1577     // Make a telescope geometry with rectagular surface
1578     const scalar detector_length = rand_length(mt1);
1579     const scalar constraint_step_size = detector_length * 1.25f;
1580
1581     mask<rect_type, test_algebra> rect{0u, detector_length * mask_scaler,
1582                                        detector_length * mask_scaler};
1583     mask<wire_type, test_algebra> wire{0u, detector_length * mask_scaler,
1584                                        detector_length * mask_scaler};
1585
1586     // Adjust overstep tolerance
1587     scalar overstep_tol =
1588         vm["overstep-tolerance-mm"].as<scalar>() * unit<scalar>::mm;
1589     overstep_tol = -math::min(math::fabs(overstep_tol), detector_length * 0.5f);
1590
1591     tel_det_config<test_algebra, rect_type, detail::helix> rectangle_cfg{
1592         rect, helix_bz};
1593     rectangle_cfg.envelope(envelope_size);
1594     rectangle_cfg.module_material(vacuum<scalar>{});
1595     rectangle_cfg.mat_thickness(0.f);
1596     rectangle_cfg.n_surfaces(2u);
1597     rectangle_cfg.length(detector_length);
1598     rectangle_cfg.volume_material(vacuum<scalar>{});
1599     rectangle_cfg.do_check(false);
1600
1601     auto alphaI = rand_alpha(mt1);
1602     auto alphaF = rand_alpha(mt1);
1603     auto betaF = math::acos(rand_cosbeta(mt1));
1604     if (rand_bool(mt1) == 0) {
1605       betaF = -betaF;
1606     }
1607     auto gammaF = rand_gamma(mt1);
1608
1609     std::array<scalar, 3u> euler_angles_I{alphaI, 0.f, 0.f};
1610     std::array<scalar, 3u> euler_angles_F{alphaF, betaF, gammaF};
1611
1612     // Without volume material
1613     auto [rect_det, rect_names] =
1614         build_telescope_detector<test_algebra>(host_mr, rectangle_cfg);
1615     const auto [euler_rect_initial, shift_rect_initial] =
1616         tilt_surface<decltype(rect_det),
1617                      decltype(rect_det)::masks::id::e_rectangle2D>(
1618             rect_det, 0u, helix_bz.dir(0.f), alphaI, 0.f, 0.f);
1619     const auto [euler_rect_final, shift_rect_final] =
1620         tilt_surface<decltype(rect_det),
1621                      decltype(rect_det)::masks::id::e_rectangle2D>(
1622             rect_det, 1u, helix_bz.dir(detector_length), alphaF, betaF, gammaF);
1623
1624     // With volume material
1625     rectangle_cfg.volume_material(volume_mat);
1626     auto [rect_det_w_mat, rect_names2] =
1627         build_telescope_detector<test_algebra>(host_mr, rectangle_cfg);
1628     [[maybe_unused]] const auto [euler_rect_initial2, shift_rect_initial2] =
1629         tilt_surface<decltype(rect_det_w_mat),
1630                      decltype(rect_det_w_mat)::masks::id::e_rectangle2D>(
1631             rect_det_w_mat, 0u, helix_bz.dir(0.f), alphaI, 0.f, 0.f);
1632     [[maybe_unused]] const auto [euler_rect_final2, shift_rect_final2] =
1633         tilt_surface<decltype(rect_det_w_mat),
1634                      decltype(rect_det_w_mat)::masks::id::e_rectangle2D>(
1635             rect_det_w_mat, 1u, helix_bz.dir(detector_length), alphaF, betaF,
1636             gammaF);
1637
1638     // Make a telescope geometry with wire surface
1639     tel_det_config<test_algebra, wire_type, detail::helix> wire_cfg{wire,
1640                                                                     helix_bz};
1641     wire_cfg.envelope(envelope_size);
1642     wire_cfg.module_material(vacuum<scalar>{});
1643     wire_cfg.mat_thickness(0.f);
1644     wire_cfg.n_surfaces(2u);
1645     wire_cfg.length(detector_length);
1646     wire_cfg.volume_material(vacuum<scalar>{});
1647     wire_cfg.do_check(false);
1648
1649     // Without volume material
1650     auto [wire_det, wire_names] =
1651         build_telescope_detector<test_algebra>(host_mr, wire_cfg);
1652     const auto [euler_wire_initial, shift_wire_initial] =
1653         tilt_surface<decltype(wire_det),
1654                      decltype(wire_det)::masks::id::e_drift_cell>(
1655             wire_det, 0u, helix_bz.dir(0.f), alphaI, 0.f, 0.f);
1656     const auto [euler_wire_final, shift_wire_final] =
1657         tilt_surface<decltype(wire_det),
1658                      decltype(wire_det)::masks::id::e_drift_cell>(
1659             wire_det, 1u, helix_bz.dir(detector_length), alphaF, betaF, gammaF);
1660
1661     // With volume material
1662     wire_cfg.volume_material(volume_mat);
1663     auto [wire_det_w_mat, wire_names2] =
1664         build_telescope_detector<test_algebra>(host_mr, wire_cfg);
1665     [[maybe_unused]] const auto [euler_wire_initial2, shift_wire_initial2] =
1666         tilt_surface<decltype(wire_det_w_mat),
1667                      decltype(wire_det_w_mat)::masks::id::e_drift_cell>(
1668             wire_det_w_mat, 0u, helix_bz.dir(0.f), alphaI, 0.f, 0.f);
1669     [[maybe_unused]] const auto [euler_wire_final2, shift_wire_final2] =
1670         tilt_surface<decltype(wire_det_w_mat),
1671                      decltype(wire_det_w_mat)::masks::id::e_drift_cell>(
1672             wire_det_w_mat, 1u, helix_bz.dir(detector_length), alphaF, betaF,
1673             gammaF);
1674
1675     // This IF block should locate after `tilt_surface()` calls for
1676     // debugging purpose
1677     if (track_count + 1 <= n_skips) {
1678       track_count++;
1679       continue;
1680     }
1681
1682     track_count++;
1683
1684     if (verbose_lvl >= 1) {
1685       std::clog << "[Event Property]" << std::endl;
1686       std::clog << "Track ID: " << track_count
1687                 << "  Number of processed tracks per thread: "
1688                 << track_count - n_skips << std::endl;
1689     }
1690     if (verbose_lvl >= 2) {
1691       std::clog << "[Detector Property]" << std::endl;
1692       std::clog << "Path length for the final surface: " << detector_length
1693                 << std::endl;
1694       std::clog << "Rect initial surface rotation: ("
1695                 << euler_rect_initial.alpha << " " << euler_rect_initial.beta
1696                 << " " << euler_rect_initial.gamma << ")" << std::endl;
1697       std::clog << "Rect initial surface shift: (" << shift_rect_initial[0u]
1698                 << " " << shift_rect_initial[1u] << " "
1699                 << shift_rect_initial[2u] << ")" << std::endl;
1700       std::clog << "Rect final surface rotation: (" << euler_rect_final.alpha
1701                 << " " << euler_rect_final.beta << " " << euler_rect_final.gamma
1702                 << ")" << std::endl;
1703       std::clog << "Rect final surface shift: (" << shift_rect_final[0u] << " "
1704                 << shift_rect_final[1u] << " " << shift_rect_final[2u] << ")"
1705                 << std::endl;
1706       std::clog << "Wire initial surface rotation: ("
1707                 << euler_wire_initial.alpha << " " << euler_wire_initial.beta
1708                 << " " << euler_wire_initial.gamma << ")" << std::endl;
1709       std::clog << "Wire initial surface shift: (" << shift_wire_initial[0u]
1710                 << " " << shift_wire_initial[1u] << " "
1711                 << shift_wire_initial[2u] << ")" << std::endl;
1712       std::clog << "Wire final surface rotation: (" << euler_wire_final.alpha
1713                 << " " << euler_wire_final.beta << " " << euler_wire_final.gamma
1714                 << ")" << std::endl;
1715       std::clog << "Wire final surface shift: (" << shift_wire_final[0u] << " "
1716                 << shift_wire_final[1u] << " " << shift_wire_final[2u] << ")"
1717                 << std::endl;
1718     }
1719     if (verbose_lvl >= 3) {
1720       std::clog << "[Track Property]" << std::endl;
1721       std::clog << "Phi: " << vector::phi(track.dir()) << std::endl;
1722       std::clog << "Theta: " << vector::theta(track.dir()) << std::endl;
1723       std::clog << "Mom: " << track.p(ptc.charge()) << std::endl;
1724     }
1725
1726     /**********************************
1727      * Rectangular telescope geometry
1728      **********************************/
1729
1730     if (verbose_lvl >= 3 && !skip_rect) {
1731       std::clog << "Simulating rectangular telescope..." << std::endl;
1732     }
1733
1734     // Get initial parameter
1735     const auto rect_bparam =
1736         get_initial_parameter<decltype(rect_det),
1737                               decltype(rect_det)::masks::id::e_rectangle2D>(
1738             rect_det, track, B_z, helix_tol);
1739
1740     if (!skip_rect) {
1741       auto ref_rel_diff{std::numeric_limits<scalar>::max()};
1742
1743       if (rk_tolerance_iterate_mode) {
1744         std::array<unsigned int, 5u> num_iterations{};
1745         std::array<bool, 25u> convergence{};
1746         auto differentiated_jacobian =
1747             directly_differentiate<inhom_field_rect_propagator_t>(
1748                 track_count, rect_bparam, rect_det_w_mat, detector_length,
1749                 inhom_bfield, overstep_tol, on_surface_tol, rk_tol_dis,
1750                 constraint_step_size, h_sizes_rect, num_iterations,
1751                 convergence);
1752
1753         for (std::size_t i = 0u; i < log10_tols.size(); i++) {
1754           // Rectangle Inhomogeneous field with Material
1755           evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
1756               track_count, euler_angles_I, euler_angles_F, rect_det_w_mat,
1757               detector_length, rect_bparam, inhom_bfield, overstep_tol,
1758               on_surface_tol, std::pow(10.f, log10_tols[i]), rk_tol_dis,
1759               constraint_step_size, h_sizes_rect, rect_files[i], ref_rel_diff,
1760               true, do_inspect, true, differentiated_jacobian, num_iterations,
1761               convergence);
1762
1763           dqopdqop_rel_diffs_rect[i].push_back(ref_rel_diff);
1764         }
1765       } else if (!rk_tolerance_iterate_mode) {
1766         // For helix
1767         evaluate_jacobian_difference_helix<
1768             decltype(rect_det), decltype(rect_det)::masks::id::e_rectangle2D>(
1769             track_count, euler_angles_I, euler_angles_F, rect_det,
1770             detector_length, rect_bparam, B_z, h_sizes_rect, helix_rect_file,
1771             helix_tol);
1772
1773         // Rect Const field
1774         evaluate_jacobian_difference<const_field_rect_propagator_t>(
1775             track_count, euler_angles_I, euler_angles_F, rect_det,
1776             detector_length, rect_bparam, const_bfield, overstep_tol,
1777             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1778             h_sizes_rect, const_rect_file, ref_rel_diff, true, false);
1779
1780         // Rect Inhomogeneous field
1781         evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
1782             track_count, euler_angles_I, euler_angles_F, rect_det,
1783             detector_length, rect_bparam, inhom_bfield, overstep_tol,
1784             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1785             h_sizes_rect, inhom_rect_file, ref_rel_diff, true, false);
1786
1787         // Rectangle Inhomogeneous field with Material
1788         evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
1789             track_count, euler_angles_I, euler_angles_F, rect_det_w_mat,
1790             detector_length, rect_bparam, inhom_bfield, overstep_tol,
1791             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1792             h_sizes_rect, inhom_rect_material_file, ref_rel_diff, true, false);
1793
1794         // Rectangle Inhomogeneous field with Material (Covariance
1795         // transport)
1796         evaluate_covariance_transport<inhom_field_rect_propagator_t>(
1797             track_count, euler_angles_I, euler_angles_F, rect_det_w_mat,
1798             detector_length, rect_bparam, inhom_bfield, overstep_tol,
1799             on_surface_tol, rk_tol_cov, rk_tol_dis, constraint_step_size,
1800             rect_cov_transport_file, true);
1801       }
1802     }
1803
1804     /**********************************
1805      * Wire telescope geometry
1806      **********************************/
1807
1808     if (verbose_lvl >= 3 && !skip_wire) {
1809       std::clog << "Simulating wire telescope..." << std::endl;
1810     }
1811
1812     // Get initial parameter
1813     const auto wire_bparam =
1814         get_initial_parameter<decltype(wire_det),
1815                               decltype(wire_det)::masks::id::e_drift_cell>(
1816             wire_det, track, B_z, helix_tol);
1817
1818     if (!skip_wire) {
1819       auto ref_rel_diff{std::numeric_limits<scalar>::max()};
1820
1821       if (rk_tolerance_iterate_mode) {
1822         std::array<unsigned int, 5u> num_iterations{};
1823         std::array<bool, 25u> convergence{};
1824         auto differentiated_jacobian =
1825             directly_differentiate<inhom_field_wire_propagator_t>(
1826                 track_count, wire_bparam, wire_det_w_mat, detector_length,
1827                 inhom_bfield, overstep_tol, on_surface_tol, rk_tol_dis,
1828                 constraint_step_size, h_sizes_wire, num_iterations,
1829                 convergence);
1830
1831         for (std::size_t i = 0u; i < log10_tols.size(); i++) {
1832           // Wire Inhomogeneous field with Material
1833           evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
1834               track_count, euler_angles_I, euler_angles_F, wire_det_w_mat,
1835               detector_length, wire_bparam, inhom_bfield, overstep_tol,
1836               on_surface_tol, std::pow(10.f, log10_tols[i]), rk_tol_dis,
1837               constraint_step_size, h_sizes_wire, wire_files[i], ref_rel_diff,
1838               true, do_inspect, true, differentiated_jacobian, num_iterations,
1839               convergence);
1840
1841           dqopdqop_rel_diffs_wire[i].push_back(ref_rel_diff);
1842         }
1843       } else if (!rk_tolerance_iterate_mode) {
1844         // For helix
1845         evaluate_jacobian_difference_helix<
1846             decltype(wire_det), decltype(wire_det)::masks::id::e_drift_cell>(
1847             track_count, euler_angles_I, euler_angles_F, wire_det,
1848             detector_length, wire_bparam, B_z, h_sizes_wire, helix_wire_file,
1849             helix_tol);
1850
1851         // Wire Const field
1852         evaluate_jacobian_difference<const_field_wire_propagator_t>(
1853             track_count, euler_angles_I, euler_angles_F, wire_det,
1854             detector_length, wire_bparam, const_bfield, overstep_tol,
1855             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1856             h_sizes_wire, const_wire_file, ref_rel_diff, true, false);
1857
1858         // Wire Inhomogeneous field
1859         evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
1860             track_count, euler_angles_I, euler_angles_F, wire_det,
1861             detector_length, wire_bparam, inhom_bfield, overstep_tol,
1862             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1863             h_sizes_wire, inhom_wire_file, ref_rel_diff, true, false);
1864
1865         // Wire Inhomogeneous field with Material
1866         evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
1867             track_count, euler_angles_I, euler_angles_F, wire_det_w_mat,
1868             detector_length, wire_bparam, inhom_bfield, overstep_tol,
1869             on_surface_tol, rk_tol_jac, rk_tol_dis, constraint_step_size,
1870             h_sizes_wire, inhom_wire_material_file, ref_rel_diff, true, false);
1871
1872         // Wire Inhomogeneous field with Material (Covariance transport)
1873         evaluate_covariance_transport<inhom_field_wire_propagator_t>(
1874             track_count, euler_angles_I, euler_angles_F, wire_det_w_mat,
1875             detector_length, wire_bparam, inhom_bfield, overstep_tol,
1876             on_surface_tol, rk_tol_cov, rk_tol_dis, constraint_step_size,
1877             wire_cov_transport_file, true);
1878       }
1879     }
1880   }
1881
1882   if (rk_tolerance_iterate_mode) {
1883     for (std::size_t i = 0u; i < log10_tols.size(); i++) {
1884       EXPECT_EQ(dqopdqop_rel_diffs_rect[i].size(), n_tracks);
1885       EXPECT_EQ(dqopdqop_rel_diffs_wire[i].size(), n_tracks);
1886
1887       EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_rect[i]), 1e-12f);
1888       EXPECT_LE(statistics::mean(dqopdqop_rel_diffs_rect[i]), 1e-2f);
1889       EXPECT_GE(statistics::mean(dqopdqop_rel_diffs_wire[i]), 1e-12f);
1890       EXPECT_LE(statistics::mean(dqopdqop_rel_diffs_wire[i]), 1e-2f);
1891     }
1892   }
1893
1894   // Close files
1895   helix_rect_file.close();
1896   const_rect_file.close();
1897   inhom_rect_file.close();
1898   inhom_rect_material_file.close();
1899
1900   helix_wire_file.close();
1901   const_wire_file.close();
1902   inhom_wire_file.close();
1903   inhom_wire_material_file.close();
1904
1905   rect_cov_transport_file.close();
1906   wire_cov_transport_file.close();
1907
1908   if (rk_tolerance_iterate_mode) {
1909     for (std::size_t i = 0u; i < log10_tols.size(); i++) {
1910       rect_files[i].close();
1911       wire_files[i].close();
1912     }
1913   }
1914 }