fixes after rebase

core/include/detray/geometry/detail/shape_utils.hpp

@@ -18,7 +18,7 @@ namespace detray::detail {
 /// using a small angle approximation
 template <typename scalar_t>
 constexpr scalar_t phi_tolerance(scalar_t tol, scalar_t radius) {
-    return tol / radius;
+    return radius > 0.f ? tol / radius : tol;
 }
 
 }  // namespace detray::detail

core/include/detray/navigation/detail/helix.hpp

@@ -173,10 +173,10 @@ class helix {
     DETRAY_HOST_DEVICE
     free_matrix_t jacobian(const scalar_type s) const {
 
-        auto ret = matrix::zero<free_matrix_t>();
+        free_matrix_t ret = matrix::zero<free_matrix_t>();
 
-        const auto I33 = matrix::identity<matrix_type<3, 3>>();
-        const auto Z33 = matrix::zero<matrix_type<3, 3>>();
+        const matrix_type<3, 3> I33 = matrix::identity<matrix_type<3, 3>>();
+        const matrix_type<3, 3> Z33 = matrix::zero<matrix_type<3, 3>>();
 
         // Notations
         // r = position
@@ -198,7 +198,7 @@ class helix {
         const scalar cos_ks = math::cos(_K * s);
         drdt = drdt + sin_ks / _K * I33;
 
-        auto H0 = matrix::zero<matrix_type<3, 1>>();
+        matrix_type<3, 1> H0 = matrix::zero<matrix_type<3, 1>>();
         getter::element(H0, 0u, 0u) = _h0[0u];
         getter::element(H0, 1u, 0u) = _h0[1u];
         getter::element(H0, 2u, 0u) = _h0[2u];

core/include/detray/propagator/actors/parameter_transporter.hpp

@@ -60,7 +60,7 @@ struct parameter_transporter : actor {
                     stepping().pos(), stepping().dir(), stepping.dtds(),
                     stepping.dqopds(vol_mat_ptr), trf3);
 
-            const auto correction_term =
+            const free_matrix_t correction_term =
                 matrix::identity<free_matrix_t>() + path_correction;
 
             return free_to_bound_jacobian * correction_term *

core/include/detray/propagator/detail/jacobian_cartesian.hpp

@@ -31,24 +31,27 @@ struct jacobian<cartesian2D<algebra_t>> {
     using vector3_type = dvector3D<algebra_t>;
 
     // Rotation Matrix
+    template <std::size_t ROWS, std::size_t COLS>
+    using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
     using rotation_matrix = dmatrix<algebra_t, 3, 3>;
     using bound_to_free_matrix_type = bound_to_free_matrix<algebra_t>;
     using free_to_bound_matrix_type = free_to_bound_matrix<algebra_t>;
     using free_to_path_matrix_type = free_to_path_matrix<algebra_t>;
     /// @}
 
     DETRAY_HOST_DEVICE
-    static inline auto reference_frame(const transform3_type &trf3,
-                                       const point3_type & /*pos*/,
-                                       const vector3_type & /*dir*/) {
+    static inline rotation_matrix reference_frame(
+        const transform3_type &trf3, const point3_type & /*pos*/,
+        const vector3_type & /*dir*/) {
         return trf3.rotation();
     }
 
     DETRAY_HOST_DEVICE static inline free_to_path_matrix_type path_derivative(
         const transform3_type &trf3, const point3_type & /*pos*/,
         const vector3_type &dir, const vector3_type & /*dtds*/) {
 
-        auto derivative = matrix::zero<free_to_path_matrix_type>();
+        free_to_path_matrix_type derivative =
+            matrix::zero<free_to_path_matrix_type>();
 
         const vector3_type normal = coordinate_frame::normal(trf3);
 
@@ -69,9 +72,12 @@ struct jacobian<cartesian2D<algebra_t>> {
 
         const rotation_matrix frame = reference_frame(trf3, pos, dir);
 
-        // getter::block<3, 2>(frame, 0u, 0u) is d(x,y,z)/d(loc0, loc1)
+        // Get d(x,y,z)/d(loc0, loc1)
+        const matrix_type<3, 2> bound_pos_to_free_pos_derivative =
+            getter::block<3, 2>(frame, 0u, 0u);
+
         getter::set_block(bound_to_free_jacobian,
-                          getter::block<3, 2>(frame, 0u, 0u), e_free_pos0,
+                          bound_pos_to_free_pos_derivative, e_free_pos0,
                           e_bound_loc0);
     }
 
@@ -81,11 +87,15 @@ struct jacobian<cartesian2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        const auto frameT = matrix::transpose(reference_frame(trf3, pos, dir));
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frameT = matrix::transpose(frame);
+
+        // Get d(loc0, loc1)/d(x,y,z)
+        const matrix_type<2, 3> free_pos_to_bound_pos_derivative =
+            getter::block<2, 3>(frameT, 0, 0);
 
-        // getter::block<2, 3>(frameT, 0, 0) is Get d(loc0, loc1)/d(x,y,z)
         getter::set_block(free_to_bound_jacobian,
-                          getter::block<2, 3>(frameT, 0u, 0u), e_bound_loc0,
+                          free_pos_to_bound_pos_derivative, e_bound_loc0,
                           e_free_pos0);
     }
 

core/include/detray/propagator/detail/jacobian_cylindrical.hpp

@@ -32,6 +32,8 @@ struct jacobian<cylindrical2D<algebra_t>> {
     using vector3_type = dvector3D<algebra_t>;
 
     // Rotation Matrix
+    template <std::size_t ROWS, std::size_t COLS>
+    using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
     using rotation_matrix = dmatrix<algebra_t, 3, 3>;
     using bound_to_free_matrix_type = bound_to_free_matrix<algebra_t>;
     using free_to_bound_matrix_type = free_to_bound_matrix<algebra_t>;
@@ -42,10 +44,10 @@ struct jacobian<cylindrical2D<algebra_t>> {
                                                   const point3_type &pos,
                                                   const vector3_type &dir) {
 
-        auto rot = matrix::zero<rotation_matrix>();
+        rotation_matrix rot = matrix::zero<rotation_matrix>();
 
         // y axis of the new frame is the z axis of cylindrical coordinate
-        const auto &new_yaxis = trf3.z();
+        const vector3_type new_yaxis = trf3.z();
 
         // z axis of the new frame is the vector normal to the cylinder surface
         const point3_type local =
@@ -55,9 +57,13 @@ struct jacobian<cylindrical2D<algebra_t>> {
         // x axis
         const vector3_type new_xaxis = vector::cross(new_yaxis, new_zaxis);
 
-        getter::set_block(rot, new_xaxis, 0u, 0u);
+        getter::element(rot, 0u, 0u) = new_xaxis[0];
+        getter::element(rot, 1u, 0u) = new_xaxis[1];
+        getter::element(rot, 2u, 0u) = new_xaxis[2];
         getter::set_block(rot, new_yaxis, 0u, 1u);
-        getter::set_block(rot, new_zaxis, 0u, 2u);
+        getter::element(rot, 0u, 2u) = new_zaxis[0];
+        getter::element(rot, 1u, 2u) = new_zaxis[1];
+        getter::element(rot, 2u, 2u) = new_zaxis[2];
 
         return rot;
     }
@@ -66,13 +72,15 @@ struct jacobian<cylindrical2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir, const vector3_type & /*dtds*/) {
 
-        auto derivative = matrix::zero<free_to_path_matrix_type>();
+        free_to_path_matrix_type derivative =
+            matrix::zero<free_to_path_matrix_type>();
 
         const point3_type local =
             coordinate_frame::global_to_local_3D(trf3, pos, dir);
         const vector3_type normal = coordinate_frame::normal(trf3, local);
 
-        const vector3_type pos_term = -1.f / vector::dot(normal, dir) * normal;
+        const vector3_type pos_term =
+            (-1.f / vector::dot(normal, dir)) * normal;
 
         getter::element(derivative, 0u, e_free_pos0) = pos_term[0];
         getter::element(derivative, 0u, e_free_pos1) = pos_term[1];
@@ -87,11 +95,14 @@ struct jacobian<cylindrical2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        const auto frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
+
+        // Get d(x,y,z)/d(loc0, loc1)
+        const matrix_type<3, 2> bound_pos_to_free_pos_derivative =
+            getter::block<3, 2>(frame, 0u, 0u);
 
-        // getter::block<3, 2>(frame, 0u, 0u) is d(x,y,z)/d(loc0, loc1)
         getter::set_block(bound_to_free_jacobian,
-                          getter::block<3, 2>(frame, 0u, 0u), e_free_pos0,
+                          bound_pos_to_free_pos_derivative, e_free_pos0,
                           e_bound_loc0);
     }
 
@@ -101,11 +112,15 @@ struct jacobian<cylindrical2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        const auto frameT = matrix::transpose(reference_frame(trf3, pos, dir));
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frameT = matrix::transpose(frame);
+
+        // Get d(loc0, loc1)/d(x,y,z)
+        const matrix_type<2, 3> free_pos_to_bound_pos_derivative =
+            getter::block<2, 3>(frameT, 0u, 0u);
 
-        // getter::block<2, 3>(frameT, 0, 0) is Get d(loc0, loc1)/d(x,y,z)
         getter::set_block(free_to_bound_jacobian,
-                          getter::block<2, 3>(frameT, 0u, 0u), e_bound_loc0,
+                          free_pos_to_bound_pos_derivative, e_bound_loc0,
                           e_free_pos0);
     }
 

core/include/detray/propagator/detail/jacobian_engine.hpp

@@ -49,7 +49,8 @@ requires std::is_object_v<typename frame_t::loc_point> struct jacobian_engine {
         const bound_parameters_vector<algebra_type>& bound_vec) {
 
         // Declare jacobian for bound to free coordinate transform
-        auto jac_to_global = matrix::zero<bound_to_free_matrix_type>();
+        bound_to_free_matrix_type jac_to_global =
+            matrix::zero<bound_to_free_matrix_type>();
 
         // Get trigonometric values
         const scalar_type theta{bound_vec.theta()};
@@ -99,7 +100,8 @@ requires std::is_object_v<typename frame_t::loc_point> struct jacobian_engine {
         const free_track_parameters<algebra_type>& free_params) {
 
         // Declare jacobian for bound to free coordinate transform
-        auto jac_to_local = matrix::zero<free_to_bound_matrix_type>();
+        free_to_bound_matrix_type jac_to_local =
+            matrix::zero<free_to_bound_matrix_type>();
 
         // Global position and direction
         const vector3_type pos = free_params.pos();
@@ -146,7 +148,8 @@ requires std::is_object_v<typename frame_t::loc_point> struct jacobian_engine {
         free_to_path_matrix_type path_derivative =
             jacobian_t::path_derivative(trf3, pos, dir, dtds);
 
-        auto derivative = matrix::zero<path_to_free_matrix_type>();
+        path_to_free_matrix_type derivative =
+            matrix::zero<path_to_free_matrix_type>();
         getter::element(derivative, e_free_pos0, 0u) = dir[0];
         getter::element(derivative, e_free_pos1, 0u) = dir[1];
         getter::element(derivative, e_free_pos2, 0u) = dir[2];

core/include/detray/propagator/detail/jacobian_line.hpp

@@ -43,10 +43,10 @@ struct jacobian<line2D<algebra_t>> {
                                                   const point3_type & /*pos*/,
                                                   const vector3_type &dir) {
 
-        auto rot = matrix::zero<rotation_matrix>();
+        rotation_matrix rot = matrix::zero<rotation_matrix>();
 
         // y axis of the new frame is the z axis of line coordinate
-        const vector3_type &new_yaxis = trf3.z();
+        const vector3_type new_yaxis = trf3.z();
 
         // x axis of the new frame is (yaxis x track direction)
         const vector3_type new_xaxis =
@@ -55,9 +55,13 @@ struct jacobian<line2D<algebra_t>> {
         // z axis
         const vector3_type new_zaxis = vector::cross(new_xaxis, new_yaxis);
 
-        getter::set_block(rot, new_xaxis, 0u, 0u);
+        getter::element(rot, 0u, 0u) = new_xaxis[0];
+        getter::element(rot, 1u, 0u) = new_xaxis[1];
+        getter::element(rot, 2u, 0u) = new_xaxis[2];
         getter::set_block(rot, new_yaxis, 0u, 1u);
-        getter::set_block(rot, new_zaxis, 0u, 2u);
+        getter::element(rot, 0u, 2u) = new_zaxis[0];
+        getter::element(rot, 1u, 2u) = new_zaxis[1];
+        getter::element(rot, 2u, 2u) = new_zaxis[2];
 
         return rot;
     }
@@ -66,7 +70,8 @@ struct jacobian<line2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir, const vector3_type &dtds) {
 
-        auto derivative = matrix::zero<free_to_path_matrix_type>();
+        free_to_path_matrix_type derivative =
+            matrix::zero<free_to_path_matrix_type>();
 
         // The vector between position and center
         const point3_type center = trf3.translation();
@@ -82,7 +87,7 @@ struct jacobian<line2D<algebra_t>> {
         // Cosine of angle between momentum direction and local frame z axis
         const scalar_type dz = vector::dot(local_zaxis, dir);
 
-        // Local x axis component of pc:
+        // local x axis component of pc:
         const vector3_type pc_x = pc - pz * local_zaxis;
 
         const scalar_type norm =
@@ -107,11 +112,14 @@ struct jacobian<line2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        const auto frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
+
+        // Get d(x,y,z)/d(loc0, loc1)
+        const auto bound_pos_to_free_pos_derivative =
+            getter::block<3, 2>(frame, 0u, 0u);
 
-        // getter::block<3, 2>(frame, 0u, 0u) is d(x,y,z)/d(loc0, loc1)
         getter::set_block(bound_to_free_jacobian,
-                          getter::block<3, 2>(frame, 0u, 0u), e_free_pos0,
+                          bound_pos_to_free_pos_derivative, e_free_pos0,
                           e_bound_loc0);
     }
 
@@ -121,11 +129,15 @@ struct jacobian<line2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        const auto frameT = matrix::transpose(reference_frame(trf3, pos, dir));
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frameT = matrix::transpose(frame);
+
+        // Get d(loc0, loc1)/d(x,y,z)
+        const auto free_pos_to_bound_pos_derivative =
+            getter::block<2, 3>(frameT, 0u, 0u);
 
-        // getter::block<2, 3>(frameT, 0, 0) is Get d(loc0, loc1)/d(x,y,z)
         getter::set_block(free_to_bound_jacobian,
-                          getter::block<2, 3>(frameT, 0u, 0u), e_bound_loc0,
+                          free_pos_to_bound_pos_derivative, e_bound_loc0,
                           e_free_pos0);
     }
 
@@ -135,30 +147,29 @@ struct jacobian<line2D<algebra_t>> {
         const transform3_type &trf3, const point3_type &pos,
         const vector3_type &dir) {
 
-        // Local 2D position
         const auto local2 = coordinate_frame::global_to_local(trf3, pos, dir);
 
         // Reference frame
-        const auto frame = reference_frame(trf3, pos, dir);
+        const rotation_matrix frame = reference_frame(trf3, pos, dir);
 
-        // New x, y and z axes
-        const auto new_xaxis = getter::vector<3>(frame, 0u, 0u);
-        const auto new_yaxis = getter::vector<3>(frame, 0u, 1u);
-        const auto new_zaxis = getter::vector<3>(frame, 0u, 2u);
+        // New x, y and z-axis
+        const vector3_type new_xaxis = getter::vector<3>(frame, 0u, 0u);
+        const vector3_type new_yaxis = getter::vector<3>(frame, 0u, 1u);
+        const vector3_type new_zaxis = getter::vector<3>(frame, 0u, 2u);
 
         // The projection of direction onto ref frame normal
         const scalar_type ipdn{1.f / vector::dot(dir, new_zaxis)};
 
         // d(n_x,n_y,n_z)/dPhi
-        const auto dNdPhi =
+        const vector3_type dNdPhi =
             getter::vector<3>(bound_to_free_jacobian, e_free_dir0, e_bound_phi);
 
         // Get new_yaxis X d(n_x,n_y,n_z)/dPhi
-        auto y_cross_dNdPhi = vector::cross(new_yaxis, dNdPhi);
+        const vector3_type y_cross_dNdPhi = vector::cross(new_yaxis, dNdPhi);
 
         // d(n_x,n_y,n_z)/dTheta
-        const auto dNdTheta = getter::vector<3>(bound_to_free_jacobian,
-                                                e_free_dir0, e_bound_theta);
+        const vector3_type dNdTheta = getter::vector<3>(
+            bound_to_free_jacobian, e_free_dir0, e_bound_theta);
 
         // Build the cross product of d(D)/d(eBoundPhi) components with y axis
         auto y_cross_dNdTheta = vector::cross(new_yaxis, dNdTheta);
@@ -177,10 +188,18 @@ struct jacobian<line2D<algebra_t>> {
         theta_to_free_pos_derivative = C * theta_to_free_pos_derivative;
 
         // Set the jacobian components
-        getter::set_block(bound_to_free_jacobian, phi_to_free_pos_derivative,
-                          e_free_pos0, e_bound_phi);
-        getter::set_block(bound_to_free_jacobian, theta_to_free_pos_derivative,
-                          e_free_pos0, e_bound_theta);
+        getter::element(bound_to_free_jacobian, e_free_pos0, e_bound_phi) =
+            phi_to_free_pos_derivative[0];
+        getter::element(bound_to_free_jacobian, e_free_pos1, e_bound_phi) =
+            phi_to_free_pos_derivative[1];
+        getter::element(bound_to_free_jacobian, e_free_pos2, e_bound_phi) =
+            phi_to_free_pos_derivative[2];
+        getter::element(bound_to_free_jacobian, e_free_pos0, e_bound_theta) =
+            theta_to_free_pos_derivative[0];
+        getter::element(bound_to_free_jacobian, e_free_pos1, e_bound_theta) =
+            theta_to_free_pos_derivative[1];
+        getter::element(bound_to_free_jacobian, e_free_pos2, e_bound_theta) =
+            theta_to_free_pos_derivative[2];
     }
 };