diff --git a/core/CMakeLists.txt b/core/CMakeLists.txt
index c9662bac0..fae9a8b2c 100644
--- a/core/CMakeLists.txt
+++ b/core/CMakeLists.txt
@@ -71,13 +71,13 @@ detray_add_library( detray_core core
    "include/detray/materials/material_slab.hpp"
    "include/detray/materials/predefined_materials.hpp"
    # propagator include(s)
-   "include/detray/propagator/aborters.hpp"
+   "include/detray/propagator/actors/aborters.hpp"
+   "include/detray/propagator/actors/parameter_transporter.hpp"
+   "include/detray/propagator/actors/parameter_resetter.hpp"
    "include/detray/propagator/actor_chain.hpp"
    "include/detray/propagator/base_actor.hpp"
    "include/detray/propagator/base_stepper.hpp"
    "include/detray/propagator/constrained_step.hpp"
-   "include/detray/propagator/detail/covariance_engine.hpp"
-   "include/detray/propagator/detail/jacobian_engine.hpp"
    "include/detray/propagator/line_stepper.hpp"
    "include/detray/propagator/navigation_policies.hpp"
    "include/detray/propagator/navigator.hpp"
@@ -107,7 +107,6 @@ detray_add_library( detray_core core
    "include/detray/tools/local_object_finder.hpp"
    # tracks include(s)
    "include/detray/tracks/detail/track_helper.hpp"
-   "include/detray/tracks/detail/trigonometrics.hpp"
    "include/detray/tracks/bound_track_parameters.hpp"
    "include/detray/tracks/free_track_parameters.hpp"
    "include/detray/tracks/tracks.hpp"
diff --git a/core/include/detray/coordinates/cartesian2.hpp b/core/include/detray/coordinates/cartesian2.hpp
index 849259988..4d957056c 100644
--- a/core/include/detray/coordinates/cartesian2.hpp
+++ b/core/include/detray/coordinates/cartesian2.hpp
@@ -21,29 +21,31 @@ struct cartesian2 final : public coordinate_base<cartesian2, transform3_t> {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Base type
+    // Base type
     using base_type = coordinate_base<cartesian2, transform3_t>;
-    /// Sclar type
+    // Sclar type
     using scalar_type = typename base_type::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename base_type::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename base_type::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename base_type::vector3;
-    /// Matrix actor
+    // Matrix operator
     using matrix_operator = typename base_type::matrix_operator;
-    /// Matrix size type
+    // Rotation Matrix
+    using rotation_matrix = typename base_type::rotation_matrix;
+    // Matrix size type
     using size_type = typename base_type::size_type;
-    /// 2D matrix type
+    // 2D matrix type
     template <size_type ROWS, size_type COLS>
     using matrix_type = typename base_type::template matrix_type<ROWS, COLS>;
-
-    // Trigonometrics
-    using trigonometrics = typename base_type::trigonometrics;
     // Vector types
     using bound_vector = typename base_type::bound_vector;
     using free_vector = typename base_type::free_vector;
+    // Matrix types
+    using free_to_bound_matrix = typename base_type::free_to_bound_matrix;
+    using bound_to_free_matrix = typename base_type::bound_to_free_matrix;
 
     /// @}
 
@@ -63,9 +65,9 @@ struct cartesian2 final : public coordinate_base<cartesian2, transform3_t> {
     /** This method transform from a point from global cartesian 3D frame to a
      * local 2D cartesian point */
     DETRAY_HOST_DEVICE
-    inline point2 global_to_local(const transform3_t &trf, const point3 &p,
+    inline point2 global_to_local(const transform3_t &trf3, const point3 &p,
                                   const vector3 & /*d*/) const {
-        const auto local3 = trf.point_to_local(p);
+        const auto local3 = trf3.point_to_local(p);
         return this->operator()(local3);
     }
 
@@ -73,30 +75,72 @@ struct cartesian2 final : public coordinate_base<cartesian2, transform3_t> {
      * cartesian 3D frame*/
     template <typename mask_t>
     DETRAY_HOST_DEVICE inline point3 local_to_global(
-        const transform3_t &trf, const mask_t & /*mask*/, const point2 &p,
+        const transform3_t &trf3, const mask_t & /*mask*/, const point2 &p,
         const vector3 & /*d*/) const {
-        return trf.point_to_global(point3{p[0], p[1], 0.});
+        return trf3.point_to_global(point3{p[0], p[1], 0.});
     }
 
-    DETRAY_HOST_DEVICE
-    inline matrix_type<3, 2> bound_to_free_rotation(
-        const transform3_t &trf3, const trigonometrics & /*t*/) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline vector3 normal(const transform3_t &trf3,
+                                             const mask_t & /*mask*/,
+                                             const point3 & /*pos*/,
+                                             const vector3 & /*dir*/) const {
+        vector3 ret;
+        const matrix_type<3, 1> n =
+            matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
+        ret[0] = matrix_operator().element(n, 0, 0);
+        ret[1] = matrix_operator().element(n, 1, 0);
+        ret[2] = matrix_operator().element(n, 2, 0);
+        return ret;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline rotation_matrix reference_frame(
+        const transform3_t &trf3, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        return trf3.rotation();
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_pos_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const rotation_matrix frame = reference_frame(trf3, mask, pos, dir);
 
         // Get d(x,y,z)/d(loc0, loc1)
-        return matrix_operator().template block<3, 2>(trf3.matrix(), 0, 0);
+        const matrix_type<3, 2> bound_pos_to_free_pos_derivative =
+            matrix_operator().template block<3, 2>(frame, 0, 0);
+
+        matrix_operator().template set_block(bound_to_free_jacobian,
+                                             bound_pos_to_free_pos_derivative,
+                                             e_free_pos0, e_bound_loc0);
     }
 
-    DETRAY_HOST_DEVICE
-    inline matrix_type<2, 3> free_to_bound_rotation(
-        const transform3_t &trf3, const trigonometrics & /*t*/) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_free_pos_to_bound_pos_derivative(
+        free_to_bound_matrix &free_to_bound_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
 
-        // Get transpose of transform3 matrix
-        const auto trf3T = matrix_operator().transpose(trf3);
+        const rotation_matrix frame = reference_frame(trf3, mask, pos, dir);
+        const rotation_matrix frameT = matrix_operator().transpose(frame);
 
         // Get d(loc0, loc1)/d(x,y,z)
-        return matrix_operator().template block<2, 3>(trf3T.matrix(), 0, 0);
+        const matrix_type<2, 3> free_pos_to_bound_pos_derivative =
+            matrix_operator().template block<2, 3>(frameT, 0, 0);
+
+        matrix_operator().template set_block(free_to_bound_jacobian,
+                                             free_pos_to_bound_pos_derivative,
+                                             e_bound_loc0, e_free_pos0);
     }
 
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_angle_to_free_pos_derivative(
+        bound_to_free_matrix & /*bound_to_free_jacobian*/,
+        const transform3_t & /*trf3*/, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        // Do nothing
+    }
 };  // struct cartesian2
 
 }  // namespace detray
\ No newline at end of file
diff --git a/core/include/detray/coordinates/coordinate_base.hpp b/core/include/detray/coordinates/coordinate_base.hpp
index f7d5ff981..dea9b866a 100644
--- a/core/include/detray/coordinates/coordinate_base.hpp
+++ b/core/include/detray/coordinates/coordinate_base.hpp
@@ -9,8 +9,13 @@
 
 // Project include(s).
 #include "detray/definitions/qualifiers.hpp"
+#include "detray/intersection/detail/trajectories.hpp"
+#include "detray/propagator/base_stepper.hpp"
 #include "detray/tracks/detail/track_helper.hpp"
 
+// System include(s).
+#include <climits>
+
 namespace detray {
 
 /** Coordinate base struct
@@ -21,26 +26,28 @@ struct coordinate_base {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Scalar type
+    // Scalar type
     using scalar_type = typename transform3_t::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename transform3_t::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename transform3_t::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename transform3_t::vector3;
-    /// Matrix actor
+    // Matrix operator
     using matrix_operator = typename transform3_t::matrix_actor;
-    /// Matrix size type
+    // Matrix size type
     using size_type = typename matrix_operator::size_ty;
-    /// 2D matrix type
+    // 2D matrix type
     template <size_type ROWS, size_type COLS>
     using matrix_type =
         typename matrix_operator::template matrix_type<ROWS, COLS>;
-    /// Shorthand vector/matrix types related to bound track parameters.
+    // Rotation Matrix
+    using rotation_matrix = matrix_type<3, 3>;
+    // Shorthand vector/matrix types related to bound track parameters.
     using bound_vector = matrix_type<e_bound_size, 1>;
     using bound_matrix = matrix_type<e_bound_size, e_bound_size>;
-    /// Mapping from bound track parameters.
+    // Mapping from bound track parameters.
     using bound_to_free_matrix = matrix_type<e_free_size, e_bound_size>;
     // Shorthand vector/matrix types related to free track parameters.
     using free_vector = matrix_type<e_free_size, 1>;
@@ -50,8 +57,6 @@ struct coordinate_base {
     using free_to_path_matrix = matrix_type<1, e_free_size>;
     // Track helper
     using track_helper = detail::track_helper<matrix_operator>;
-    // Trigonometrics
-    using trigonometrics = typename track_helper::trigonometrics;
 
     /// @}
 
@@ -104,75 +109,100 @@ struct coordinate_base {
         return free_vec;
     }
 
+    template <typename mask_t>
     DETRAY_HOST_DEVICE inline bound_to_free_matrix bound_to_free_jacobian(
-        const transform3_t& trf3, const bound_vector& bound_vec) {
+        const transform3_t& trf3, const mask_t& mask,
+        const bound_vector& bound_vec) const {
 
         // Declare jacobian for bound to free coordinate transform
         bound_to_free_matrix jac_to_global =
             matrix_operator().template zero<e_free_size, e_bound_size>();
 
         // Get trigonometric values
-        const auto t = track_helper().get_trigonometrics(bound_vec);
+        const scalar_type theta =
+            matrix_operator().element(bound_vec, e_bound_theta, 0);
+        const scalar_type phi =
+            matrix_operator().element(bound_vec, e_bound_phi, 0);
 
-        // Get d(x,y,z)/d(loc0, loc1)
-        const matrix_type<3, 2> bound_to_free_rotation =
-            Derived<transform3_t>().bound_to_free_rotation(trf3, t);
+        const scalar_type cos_theta = std::cos(theta);
+        const scalar_type sin_theta = std::sin(theta);
+        const scalar_type cos_phi = std::cos(phi);
+        const scalar_type sin_phi = std::sin(phi);
+
+        // Global position and direction
+        const auto free_vec = bound_to_free_vector(trf3, mask, bound_vec);
+
+        const vector3 pos = track_helper().pos(free_vec);
+        const vector3 dir = track_helper().dir(free_vec);
 
-        matrix_operator().template set_block(
-            jac_to_global, bound_to_free_rotation, e_free_pos0, e_bound_loc0);
+        // Set d(x,y,z)/d(loc0, loc1)
+        Derived<transform3_t>().set_bound_pos_to_free_pos_derivative(
+            jac_to_global, trf3, mask, pos, dir);
 
+        // Set d(bound time)/d(free time)
         matrix_operator().element(jac_to_global, e_free_time, e_bound_time) = 1;
+
+        // Set d(n_x,n_y,n_z)/d(phi, theta)
         matrix_operator().element(jac_to_global, e_free_dir0, e_bound_phi) =
-            -1 * t.sin_theta * t.sin_phi;
+            -1 * sin_theta * sin_phi;
         matrix_operator().element(jac_to_global, e_free_dir0, e_bound_theta) =
-            t.cos_theta * t.cos_phi;
+            cos_theta * cos_phi;
         matrix_operator().element(jac_to_global, e_free_dir1, e_bound_phi) =
-            t.sin_theta * t.cos_phi;
+            sin_theta * cos_phi;
         matrix_operator().element(jac_to_global, e_free_dir1, e_bound_theta) =
-            t.cos_theta * t.sin_phi;
+            cos_theta * sin_phi;
         matrix_operator().element(jac_to_global, e_free_dir2, e_bound_theta) =
-            -1 * t.sin_theta;
+            -1 * sin_theta;
         matrix_operator().element(jac_to_global, e_free_qoverp,
                                   e_bound_qoverp) = 1;
 
+        // Set d(x,y,z)/d(phi, theta)
+        Derived<transform3_t>().set_bound_angle_to_free_pos_derivative(
+            jac_to_global, trf3, mask, pos, dir);
+
         return jac_to_global;
     }
 
-    DETRAY_HOST_DEVICE
-    inline free_to_bound_matrix free_to_bound_jacobian(
-        const transform3_t& trf3, const free_vector& free_vec) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline free_to_bound_matrix free_to_bound_jacobian(
+        const transform3_t& trf3, const mask_t& mask,
+        const free_vector& free_vec) const {
 
         // Declare jacobian for bound to free coordinate transform
         free_to_bound_matrix jac_to_local =
             matrix_operator().template zero<e_bound_size, e_free_size>();
 
-        // Free direction
+        // Global position and direction
+        const vector3 pos = track_helper().pos(free_vec);
         const vector3 dir = track_helper().dir(free_vec);
 
-        // Get trigonometric values
-        const auto t = track_helper().get_trigonometrics(free_vec);
+        const scalar_type theta = getter::theta(dir);
+        const scalar_type phi = getter::phi(dir);
 
-        // Get d(loc0, loc1)/d(x,y,z)
-        const matrix_type<2, 3> free_to_bound_rotation =
-            Derived<transform3_t>().free_to_bound_rotation(trf3, t);
+        const scalar_type cos_theta = std::cos(theta);
+        const scalar_type sin_theta = std::sin(theta);
+        const scalar_type cos_phi = std::cos(phi);
+        const scalar_type sin_phi = std::sin(phi);
 
-        matrix_operator().template set_block(
-            jac_to_local, free_to_bound_rotation, e_bound_loc0, e_free_pos0);
+        // Set d(loc0, loc1)/d(x,y,z)
+        Derived<transform3_t>().set_free_pos_to_bound_pos_derivative(
+            jac_to_local, trf3, mask, pos, dir);
 
-        // Set d(Free time)/d(Bound time)
+        // Set d(free time)/d(bound time)
         matrix_operator().element(jac_to_local, e_bound_time, e_free_time) = 1;
 
-        // Set d(phi, theta)/d(free dir)
+        // Set d(phi, theta)/d(n_x, n_y, n_z)
+        // @note This codes have a serious bug when theta is equal to zero...
         matrix_operator().element(jac_to_local, e_bound_phi, e_free_dir0) =
-            -1. * t.sin_phi / t.sin_theta;
+            -1. * sin_phi / sin_theta;
         matrix_operator().element(jac_to_local, e_bound_phi, e_free_dir1) =
-            t.cos_phi / t.sin_theta;
+            cos_phi / sin_theta;
         matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir0) =
-            t.cos_phi * t.cos_theta;
+            cos_phi * cos_theta;
         matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir1) =
-            t.sin_phi * t.cos_theta;
+            sin_phi * cos_theta;
         matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir2) =
-            -1 * t.sin_theta;
+            -1 * sin_theta;
 
         // Set d(Free Qop)/d(Bound Qop)
         matrix_operator().element(jac_to_local, e_bound_qoverp, e_free_qoverp) =
@@ -180,6 +210,142 @@ struct coordinate_base {
 
         return jac_to_local;
     }
+
+    template <typename mask_t, typename stepper_state_t>
+    DETRAY_HOST_DEVICE inline free_matrix path_correction(
+        const stepper_state_t& stepping, const transform3_t& trf3,
+        const mask_t& mask) {
+
+        free_matrix path_correction =
+            matrix_operator().template zero<e_free_size, e_free_size>();
+
+        // Position and direction
+        const auto pos = stepping().pos();
+        const auto dir = stepping().dir();
+
+        // dir
+        matrix_type<1, 3> t;
+        matrix_operator().element(t, 0, 0) = dir[0];
+        matrix_operator().element(t, 0, 1) = dir[1];
+        matrix_operator().element(t, 0, 2) = dir[2];
+
+        // Surface normal vector (w)
+        matrix_type<1, 3> w;
+        const auto normal =
+            Derived<transform3_t>().normal(trf3, mask, pos, dir);
+        matrix_operator().element(w, 0, 0) = normal[0];
+        matrix_operator().element(w, 0, 1) = normal[1];
+        matrix_operator().element(w, 0, 2) = normal[2];
+
+        // w dot t
+        const scalar_type wt = vector::dot(normal, dir);
+
+        // transpose of t
+        const matrix_type<3, 1> t_T = matrix_operator().transpose(t);
+
+        // r correction term
+        const matrix_type<1, 3> r_term = -1. / wt * w;
+
+        // dr/dr0
+        const matrix_type<3, 3> drdr0 = t_T * r_term;
+
+        matrix_operator().template set_block<3, 3>(path_correction, drdr0,
+                                                   e_free_pos0, e_free_pos0);
+
+        // @note: Helical correction which doesn't exist in ACTS main
+        if constexpr (stepper_state_t::id == stepping::id::e_rk) {
+            using helix = detail::helix<transform3_t>;
+
+            // Path length
+            const auto s = stepping._s;
+
+            // helix
+            helix hlx(stepping(), &stepping._step_data.b_last);
+
+            // B field at the destination surface
+            matrix_type<1, 3> h;
+            matrix_operator().element(h, 0, 0) = hlx._h0[0];
+            matrix_operator().element(h, 0, 1) = hlx._h0[1];
+            matrix_operator().element(h, 0, 2) = hlx._h0[2];
+            // matrix_operator().set_block(h, hlx._h0, 0, 0);
+
+            // Normalized vector of h X t
+            matrix_type<1, 3> n;
+            matrix_operator().element(n, 0, 0) = hlx._n0[0];
+            matrix_operator().element(n, 0, 1) = hlx._n0[1];
+            matrix_operator().element(n, 0, 2) = hlx._n0[2];
+            // matrix_operator().set_block(n, hlx._n0, 0, 0);
+
+            // w cross h
+            matrix_type<1, 3> wh;
+            const auto _wh = vector::cross(normal, hlx._h0);
+            matrix_operator().element(wh, 0, 0) = _wh[0];
+            matrix_operator().element(wh, 0, 1) = _wh[1];
+            matrix_operator().element(wh, 0, 2) = _wh[2];
+            // matrix_operator().set_block(wh, vector::cross(normal, hlx._h0),
+            // 0, 0);
+
+            // Alpha
+            const scalar_type A = hlx._alpha;
+
+            // K
+            const scalar_type K = hlx._K;
+
+            // Ks = K*s
+            const scalar_type Ks = K * s;
+
+            // t correction term
+            matrix_type<1, 3> t_term = matrix_operator().template zero<1, 3>();
+            t_term = t_term +
+                     (Ks - std::sin(Ks)) / K * vector::dot(hlx._h0, normal) * h;
+
+            t_term = t_term + std::sin(Ks) / K * w;
+
+            t_term = t_term + (1 - std::cos(Ks)) / K * wh;
+
+            t_term = -1. / wt * t_term;
+
+            // qoverp correction term
+            const scalar_type L_term =
+                -1. / wt * A * Ks / hlx.qop() * vector::dot(normal, hlx._n0);
+
+            // dr/dt0
+            const matrix_type<3, 3> drdt0 = t_T * t_term;
+
+            // dr/dL0 (L = qoverp)
+            const matrix_type<3, 1> drdL0 = t_T * L_term;
+
+            // transpose of n
+            const matrix_type<3, 1> n_T = matrix_operator().transpose(n);
+
+            // dt/dr0
+            const scalar_type AK = A * K;
+            const matrix_type<3, 3> dtdr0 = AK * n_T * r_term;
+
+            // dt/dt0
+            const matrix_type<3, 3> dtdt0 = AK * n_T * t_term;
+
+            // dt/dL0
+            const matrix_type<3, 1> dtdL0 = AK * n_T * L_term;
+
+            matrix_operator().template set_block<3, 3>(
+                path_correction, drdt0, e_free_pos0, e_free_dir0);
+
+            matrix_operator().template set_block<3, 1>(
+                path_correction, drdL0, e_free_pos0, e_free_qoverp);
+
+            matrix_operator().template set_block<3, 3>(
+                path_correction, dtdr0, e_free_dir0, e_free_pos0);
+
+            matrix_operator().template set_block<3, 3>(
+                path_correction, dtdt0, e_free_dir0, e_free_dir0);
+
+            matrix_operator().template set_block<3, 1>(
+                path_correction, dtdL0, e_free_dir0, e_free_qoverp);
+        }
+
+        return path_correction;
+    }
 };
 
 }  // namespace detray
\ No newline at end of file
diff --git a/core/include/detray/coordinates/cylindrical2.hpp b/core/include/detray/coordinates/cylindrical2.hpp
index f08f9fde4..9d83d1429 100644
--- a/core/include/detray/coordinates/cylindrical2.hpp
+++ b/core/include/detray/coordinates/cylindrical2.hpp
@@ -24,16 +24,31 @@ struct cylindrical2 : public coordinate_base<cylindrical2, transform3_t> {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Base type
+    // Base type
     using base_type = coordinate_base<cylindrical2, transform3_t>;
-    /// Sclar type
+    // Sclar type
     using scalar_type = typename transform3_t::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename transform3_t::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename transform3_t::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename transform3_t::vector3;
+    // Matrix actor
+    using matrix_operator = typename base_type::matrix_operator;
+    // Matrix size type
+    using size_type = typename base_type::size_type;
+    // 2D matrix type
+    template <size_type ROWS, size_type COLS>
+    using matrix_type = typename base_type::template matrix_type<ROWS, COLS>;
+    // Rotation Matrix
+    using rotation_matrix = typename base_type::rotation_matrix;
+    // Vector types
+    using bound_vector = typename base_type::bound_vector;
+    using free_vector = typename base_type::free_vector;
+    // Matrix types
+    using free_to_bound_matrix = typename base_type::free_to_bound_matrix;
+    using bound_to_free_matrix = typename base_type::bound_to_free_matrix;
 
     /// @}
 
@@ -69,6 +84,90 @@ struct cylindrical2 : public coordinate_base<cylindrical2, transform3_t> {
         return trf.point_to_global(point3{x, y, z});
     }
 
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline vector3 normal(const transform3_t &trf3,
+                                             const mask_t &mask,
+                                             const point3 &pos,
+                                             const vector3 &dir) const {
+        const point2 local2 = this->global_to_local(trf3, pos, dir);
+        const scalar_type r = mask[0];
+        const scalar_type phi = local2[0] / r;
+
+        // normal vector in local coordinate
+        const vector3 local_normal{std::cos(phi), std::sin(phi), 0};
+
+        // normal vector in global coordinate
+        return trf3.rotation() * local_normal;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline rotation_matrix reference_frame(
+        const transform3_t &trf3, const mask_t &mask, const point3 &pos,
+        const vector3 &dir) const {
+
+        rotation_matrix rot = matrix_operator().template zero<3, 3>();
+
+        // y axis of the new frame is the z axis of cylindrical coordinate
+        const auto new_yaxis =
+            matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
+
+        // z axis of the new frame is the vector normal to the cylinder surface
+        const vector3 new_zaxis = normal(trf3, mask, pos, dir);
+
+        // x axis
+        const vector3 new_xaxis = vector::cross(new_yaxis, new_zaxis);
+
+        matrix_operator().element(rot, 0, 0) = new_xaxis[0];
+        matrix_operator().element(rot, 1, 0) = new_xaxis[1];
+        matrix_operator().element(rot, 2, 0) = new_xaxis[2];
+        matrix_operator().template set_block<3, 1>(rot, new_yaxis, 0, 1);
+        matrix_operator().element(rot, 0, 2) = new_zaxis[0];
+        matrix_operator().element(rot, 1, 2) = new_zaxis[1];
+        matrix_operator().element(rot, 2, 2) = new_zaxis[2];
+
+        return rot;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_pos_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+
+        // Get d(x,y,z)/d(loc0, loc1)
+        const auto bound_pos_to_free_pos_derivative =
+            matrix_operator().template block<3, 2>(frame, 0, 0);
+
+        matrix_operator().template set_block(bound_to_free_jacobian,
+                                             bound_pos_to_free_pos_derivative,
+                                             e_free_pos0, e_bound_loc0);
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_free_pos_to_bound_pos_derivative(
+        free_to_bound_matrix &free_to_bound_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+        const auto frameT = matrix_operator().transpose(frame);
+
+        // Get d(loc0, loc1)/d(x,y,z)
+        const auto free_pos_to_bound_pos_derivative =
+            matrix_operator().template block<2, 3>(frameT, 0, 0);
+
+        matrix_operator().template set_block(free_to_bound_jacobian,
+                                             free_pos_to_bound_pos_derivative,
+                                             e_bound_loc0, e_free_pos0);
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_angle_to_free_pos_derivative(
+        bound_to_free_matrix & /*bound_to_free_jacobian*/,
+        const transform3_t & /*trf3*/, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        // Do nothing
+    }
 };  // struct cylindrical2
 
 }  // namespace detray
\ No newline at end of file
diff --git a/core/include/detray/coordinates/line2.hpp b/core/include/detray/coordinates/line2.hpp
index b9a71c907..c5119ed7e 100644
--- a/core/include/detray/coordinates/line2.hpp
+++ b/core/include/detray/coordinates/line2.hpp
@@ -19,16 +19,33 @@ struct line2 : public coordinate_base<line2, transform3_t> {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Base type
+    // Base type
     using base_type = coordinate_base<line2, transform3_t>;
-    /// Sclar type
+    // Sclar type
     using scalar_type = typename base_type::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename base_type::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename base_type::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename base_type::vector3;
+    // Matrix actor
+    using matrix_operator = typename base_type::matrix_operator;
+    // Matrix size type
+    using size_type = typename base_type::size_type;
+    // 2D matrix type
+    template <size_type ROWS, size_type COLS>
+    using matrix_type = typename base_type::template matrix_type<ROWS, COLS>;
+    // Rotation Matrix
+    using rotation_matrix = typename base_type::rotation_matrix;
+    // Vector types
+    using bound_vector = typename base_type::bound_vector;
+    using free_vector = typename base_type::free_vector;
+    // Track Helper
+    using track_helper = typename base_type::track_helper;
+    // Matrix types
+    using free_to_bound_matrix = typename base_type::free_to_bound_matrix;
+    using bound_to_free_matrix = typename base_type::bound_to_free_matrix;
 
     /// @}
 
@@ -85,6 +102,136 @@ struct line2 : public coordinate_base<line2, transform3_t> {
 
         return locZ_in_global + p[0] * vector::normalize(r);
     }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline rotation_matrix reference_frame(
+        const transform3_t &trf3, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 &dir) const {
+
+        rotation_matrix rot = matrix_operator().template zero<3, 3>();
+
+        // y axis of the new frame is the z axis of line coordinate
+        const auto new_yaxis =
+            matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
+
+        // x axis of the new frame is (yaxis x track direction)
+        auto new_xaxis = vector::cross(new_yaxis, dir);
+        new_xaxis = vector::normalize(new_xaxis);
+
+        // z axis
+        const auto new_zaxis = vector::cross(new_xaxis, new_yaxis);
+
+        matrix_operator().element(rot, 0, 0) = new_xaxis[0];
+        matrix_operator().element(rot, 1, 0) = new_xaxis[1];
+        matrix_operator().element(rot, 2, 0) = new_xaxis[2];
+        matrix_operator().template set_block<3, 1>(rot, new_yaxis, 0, 1);
+        matrix_operator().element(rot, 0, 2) = new_zaxis[0];
+        matrix_operator().element(rot, 1, 2) = new_zaxis[1];
+        matrix_operator().element(rot, 2, 2) = new_zaxis[2];
+
+        return rot;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_pos_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+
+        // Get d(x,y,z)/d(loc0, loc1)
+        const auto bound_pos_to_free_pos_derivative =
+            matrix_operator().template block<3, 2>(frame, 0, 0);
+
+        matrix_operator().template set_block(bound_to_free_jacobian,
+                                             bound_pos_to_free_pos_derivative,
+                                             e_free_pos0, e_bound_loc0);
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_free_pos_to_bound_pos_derivative(
+        free_to_bound_matrix &free_to_bound_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+        const auto frameT = matrix_operator().transpose(frame);
+
+        // Get d(loc0, loc1)/d(x,y,z)
+        const auto free_pos_to_bound_pos_derivative =
+            matrix_operator().template block<2, 3>(frameT, 0, 0);
+
+        matrix_operator().template set_block(free_to_bound_jacobian,
+                                             free_pos_to_bound_pos_derivative,
+                                             e_bound_loc0, e_free_pos0);
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_angle_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        // local2
+        const auto local2 = this->global_to_local(trf3, pos, dir);
+
+        // Reference frame
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+
+        // new x_axis
+        const auto new_xaxis = getter::vector<3>(frame, 0, 0);
+
+        // new y_axis
+        const auto new_yaxis = getter::vector<3>(frame, 0, 1);
+
+        // new z_axis
+        const auto new_zaxis = getter::vector<3>(frame, 0, 2);
+
+        // the projection of direction onto ref frame normal
+        scalar_type ipdn = 1. / vector::dot(dir, new_zaxis);
+
+        // d(n_x,n_y,n_z)/dPhi
+        const auto dNdPhi = matrix_operator().template block<3, 1>(
+            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);
+
+        // d(n_x,n_y,n_z)/dTheta
+        const auto dNdTheta = matrix_operator().template block<3, 1>(
+            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);
+
+        const scalar_type C = ipdn * local2[0];
+        // and correct for the x axis components
+        vector3 phi_to_free_pos_derivative =
+            y_cross_dNdPhi - new_xaxis * vector::dot(new_xaxis, y_cross_dNdPhi);
+
+        phi_to_free_pos_derivative = C * phi_to_free_pos_derivative;
+
+        vector3 theta_to_free_pos_derivative =
+            y_cross_dNdTheta -
+            new_xaxis * vector::dot(new_xaxis, y_cross_dNdTheta);
+
+        theta_to_free_pos_derivative = C * theta_to_free_pos_derivative;
+
+        // Set the jacobian components
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos0,
+                                  e_bound_phi) = phi_to_free_pos_derivative[0];
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos1,
+                                  e_bound_phi) = phi_to_free_pos_derivative[1];
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos2,
+                                  e_bound_phi) = phi_to_free_pos_derivative[2];
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos0,
+                                  e_bound_theta) =
+            theta_to_free_pos_derivative[0];
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos1,
+                                  e_bound_theta) =
+            theta_to_free_pos_derivative[1];
+        matrix_operator().element(bound_to_free_jacobian, e_free_pos2,
+                                  e_bound_theta) =
+            theta_to_free_pos_derivative[2];
+    }
 };
 
 }  // namespace detray
\ No newline at end of file
diff --git a/core/include/detray/coordinates/polar2.hpp b/core/include/detray/coordinates/polar2.hpp
index abe1777c8..ec2cb38f6 100644
--- a/core/include/detray/coordinates/polar2.hpp
+++ b/core/include/detray/coordinates/polar2.hpp
@@ -22,28 +22,33 @@ struct polar2 : public coordinate_base<polar2, transform3_t> {
     /// @name Type definitions for the struct
     /// @{
 
-    /// Base type
+    // Base type
     using base_type = coordinate_base<polar2, transform3_t>;
-    /// Sclar type
+    // Sclar type
     using scalar_type = typename base_type::scalar_type;
-    /// Point in 2D space
+    // Point in 2D space
     using point2 = typename base_type::point2;
-    /// Point in 3D space
+    // Point in 3D space
     using point3 = typename base_type::point3;
-    /// Vector in 3D space
+    // Vector in 3D space
     using vector3 = typename base_type::vector3;
     /// Matrix actor
     using matrix_operator = typename base_type::matrix_operator;
     /// Matrix size type
     using size_type = typename base_type::size_type;
-    /// 2D matrix type
+    // 2D matrix type
     template <size_type ROWS, size_type COLS>
     using matrix_type = typename base_type::template matrix_type<ROWS, COLS>;
-    // Trigonometrics
-    using trigonometrics = typename base_type::trigonometrics;
+    // Rotation Matrix
+    using rotation_matrix = typename base_type::rotation_matrix;
     // Vector types
     using bound_vector = typename base_type::bound_vector;
     using free_vector = typename base_type::free_vector;
+    // Track Helper
+    using track_helper = typename base_type::track_helper;
+    // Matrix types
+    using free_to_bound_matrix = typename base_type::free_to_bound_matrix;
+    using bound_to_free_matrix = typename base_type::bound_to_free_matrix;
 
     /** This method transform from a point from 2D cartesian frame to a 2D
      * polar point */
@@ -80,18 +85,114 @@ struct polar2 : public coordinate_base<polar2, transform3_t> {
         return trf.point_to_global(point3{x, y, 0.});
     }
 
-    DETRAY_HOST_DEVICE
-    inline matrix_type<3, 2> bound_to_free_rotation(
-        const transform3_t & /*trf3*/, const trigonometrics & /*t*/) {
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline vector3 normal(const transform3_t &trf3,
+                                             const mask_t & /*mask*/,
+                                             const point3 & /*pos*/,
+                                             const vector3 & /*dir*/) const {
+        vector3 ret;
+        const auto n =
+            matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
+        ret[0] = matrix_operator().element(n, 0, 0);
+        ret[1] = matrix_operator().element(n, 1, 0);
+        ret[2] = matrix_operator().element(n, 2, 0);
+        return ret;
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline rotation_matrix reference_frame(
+        const transform3_t &trf3, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        return trf3.rotation();
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_pos_to_free_pos_derivative(
+        bound_to_free_matrix &bound_to_free_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
 
-        /*
-        matrix_type<3, 2> bound_to_free_rotation =
+        matrix_type<3, 2> bound_pos_to_free_pos_derivative =
             matrix_operator().template zero<3, 2>();
-        */
+
+        const point2 local2 = this->global_to_local(trf3, pos, dir);
+        const scalar_type lrad = local2[0];
+        const scalar_type lphi = local2[1];
+
+        const scalar_type lcos_phi = std::cos(lphi);
+        const scalar_type lsin_phi = std::sin(lphi);
+
+        // reference matrix
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+
+        // dxdu = d(x,y,z)/du
+        const matrix_type<3, 1> dxdL =
+            matrix_operator().template block<3, 1>(frame, 0, 0);
+        // dxdv = d(x,y,z)/dv
+        const matrix_type<3, 1> dydL =
+            matrix_operator().template block<3, 1>(frame, 0, 1);
+
+        const matrix_type<3, 1> col0 = dxdL * lcos_phi + dydL * lsin_phi;
+        const matrix_type<3, 1> col1 =
+            (dydL * lcos_phi - dxdL * lsin_phi) * lrad;
+
+        matrix_operator().template set_block<3, 1>(
+            bound_pos_to_free_pos_derivative, col0, e_free_pos0, e_bound_loc0);
+        matrix_operator().template set_block<3, 1>(
+            bound_pos_to_free_pos_derivative, col1, e_free_pos0, e_bound_loc1);
+
+        matrix_operator().template set_block(bound_to_free_jacobian,
+                                             bound_pos_to_free_pos_derivative,
+                                             e_free_pos0, e_bound_loc0);
     }
 
-    DETRAY_HOST_DEVICE inline matrix_type<2, 3> free_to_bound_rotation(
-        const transform3_t & /*trf3*/, const trigonometrics & /*t*/) {}
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_free_pos_to_bound_pos_derivative(
+        free_to_bound_matrix &free_to_bound_jacobian, const transform3_t &trf3,
+        const mask_t &mask, const point3 &pos, const vector3 &dir) const {
+
+        matrix_type<2, 3> free_pos_to_bound_pos_derivative =
+            matrix_operator().template zero<2, 3>();
+
+        const auto local = this->global_to_local(trf3, pos, dir);
+
+        const scalar_type lrad = local[0];
+        const scalar_type lphi = local[1];
+
+        const scalar_type lcos_phi = std::cos(lphi);
+        const scalar_type lsin_phi = std::sin(lphi);
+
+        // reference matrix
+        const auto frame = reference_frame(trf3, mask, pos, dir);
+        const auto frameT = matrix_operator().transpose(frame);
+
+        // dudG = du/d(x,y,z)
+        const matrix_type<1, 3> dudG =
+            matrix_operator().template block<1, 3>(frameT, 0, 0);
+        // dvdG = dv/d(x,y,z)
+        const matrix_type<1, 3> dvdG =
+            matrix_operator().template block<1, 3>(frameT, 1, 0);
+
+        const matrix_type<1, 3> row0 = dudG * lcos_phi + dvdG * lsin_phi;
+        const matrix_type<1, 3> row1 =
+            1. / lrad * (lcos_phi * dvdG - lsin_phi * dudG);
+
+        matrix_operator().template set_block<1, 3>(
+            free_pos_to_bound_pos_derivative, row0, e_bound_loc0, e_free_pos0);
+        matrix_operator().template set_block<1, 3>(
+            free_pos_to_bound_pos_derivative, row1, e_bound_loc1, e_free_pos0);
+
+        matrix_operator().template set_block(free_to_bound_jacobian,
+                                             free_pos_to_bound_pos_derivative,
+                                             e_bound_loc0, e_free_pos0);
+    }
+
+    template <typename mask_t>
+    DETRAY_HOST_DEVICE inline void set_bound_angle_to_free_pos_derivative(
+        bound_to_free_matrix & /*bound_to_free_jacobian*/,
+        const transform3_t & /*trf3*/, const mask_t & /*mask*/,
+        const point3 & /*pos*/, const vector3 & /*dir*/) const {
+        // Do nothing
+    }
 };
 
 }  // namespace detray
\ No newline at end of file
diff --git a/core/include/detray/intersection/detail/trajectories.hpp b/core/include/detray/intersection/detail/trajectories.hpp
index 4d5036384..bdb196bcf 100644
--- a/core/include/detray/intersection/detail/trajectories.hpp
+++ b/core/include/detray/intersection/detail/trajectories.hpp
@@ -106,6 +106,8 @@ class helix : public free_track_parameters<transform3_t> {
     // Column wise operator
     using column_wise_op = column_wise_operator<matrix_operator>;
 
+    using free_track_parameters_type::pos;
+
     DETRAY_HOST_DEVICE
     helix() = delete;
 
@@ -222,13 +224,18 @@ class helix : public free_track_parameters<transform3_t> {
             matrix_operator().template identity<3, 3>();
         const matrix_type<3, 3> Z33 = matrix_operator().template zero<3, 3>();
 
+        // Notations
+        // r = position
+        // t = direction
+        // l = qoverp
+
         // Get drdr
         auto drdr = I33;
-        matrix_operator().set_block(ret, drdr, 0, 0);
+        matrix_operator().set_block(ret, drdr, e_free_pos0, e_free_pos0);
 
         // Get dtdr
         auto dtdr = Z33;
-        matrix_operator().set_block(ret, dtdr, 4, 0);
+        matrix_operator().set_block(ret, dtdr, e_free_dir0, e_free_pos0);
 
         // Get drdt
         auto drdt = Z33;
@@ -243,7 +250,7 @@ class helix : public free_track_parameters<transform3_t> {
         drdt = drdt +
                (std::cos(_K * s) - 1) / _K * column_wise_op().cross(I33, _h0);
 
-        matrix_operator().set_block(ret, drdt, 0, 4);
+        matrix_operator().set_block(ret, drdt, e_free_pos0, e_free_dir0);
 
         // Get dtdt
         auto dtdt = Z33;
@@ -253,29 +260,28 @@ class helix : public free_track_parameters<transform3_t> {
                               matrix_operator().transpose(H0), _h0);
         dtdt = dtdt - std::sin(_K * s) * column_wise_op().cross(I33, _h0);
 
-        matrix_operator().set_block(ret, dtdt, 4, 4);
+        matrix_operator().set_block(ret, dtdt, e_free_dir0, e_free_dir0);
 
         // Get drdl
         vector3 drdl = 1 / free_track_parameters_type::qop() *
                        (s * this->dir(s) + free_track_parameters_type::pos() -
                         this->pos(s));
 
-        matrix_operator().set_block(ret, drdl, 0, 7);
+        matrix_operator().set_block(ret, drdl, e_free_pos0, e_free_qoverp);
 
         // Get dtdl
         vector3 dtdl =
             _alpha * _K * s / free_track_parameters_type::qop() * _n0;
 
-        matrix_operator().set_block(ret, dtdl, 4, 7);
+        matrix_operator().set_block(ret, dtdl, e_free_dir0, e_free_qoverp);
 
         // 3x3 and 7x7 element is 1 (Maybe?)
-        matrix_operator().element(ret, 3, 3) = 1;
-        matrix_operator().element(ret, 7, 7) = 1;
+        matrix_operator().element(ret, e_free_time, e_free_time) = 1;
+        matrix_operator().element(ret, e_free_qoverp, e_free_qoverp) = 1;
 
         return ret;
     }
 
-    private:
     /// B field
     vector3 const *_mag_field;
 
diff --git a/core/include/detray/intersection/intersection.hpp b/core/include/detray/intersection/intersection.hpp
index 2d36f8acc..0bd3c3bb7 100644
--- a/core/include/detray/intersection/intersection.hpp
+++ b/core/include/detray/intersection/intersection.hpp
@@ -59,6 +59,9 @@ struct line_plane_intersection {
     intersection::status status = intersection::status::e_missed;
     intersection::direction direction = intersection::direction::e_undefined;
 
+    // Mask index
+    dindex mask_index = dindex_invalid;
+
     // Primitive this intersection belongs to
     // TODO: rename the variable properly
     dindex index = dindex_invalid;
diff --git a/core/include/detray/intersection/intersection_kernel.hpp b/core/include/detray/intersection/intersection_kernel.hpp
index 6b603faa2..e69ff3029 100644
--- a/core/include/detray/intersection/intersection_kernel.hpp
+++ b/core/include/detray/intersection/intersection_kernel.hpp
@@ -62,6 +62,7 @@ struct intersection_initialize {
             for (auto &is : sfi) {
                 if (is.status == intersection::status::e_inside &&
                     is.path >= traj.overstep_tolerance()) {
+                    // is.mask_index = mask_index;
                     is.index = surface.volume();
                     is_container.push_back(is);
                     count++;
diff --git a/core/include/detray/masks/masks.hpp b/core/include/detray/masks/masks.hpp
index 9603b7af6..e4b55e32c 100644
--- a/core/include/detray/masks/masks.hpp
+++ b/core/include/detray/masks/masks.hpp
@@ -144,6 +144,11 @@ class mask {
                    : intersection::status::e_outside;
     }
 
+    /// @returns return local frame object
+    DETRAY_HOST_DEVICE inline constexpr local_frame_type local_frame() const {
+        return local_frame_type{};
+    }
+
     /// @returns the boundary values
     DETRAY_HOST_DEVICE
     auto values() const -> const mask_values& { return _values; }
diff --git a/core/include/detray/propagator/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp
similarity index 100%
rename from core/include/detray/propagator/aborters.hpp
rename to core/include/detray/propagator/actors/aborters.hpp
diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp
new file mode 100644
index 000000000..2c8a30e79
--- /dev/null
+++ b/core/include/detray/propagator/actors/parameter_resetter.hpp
@@ -0,0 +1,97 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// Project include(s).
+#include "detray/definitions/qualifiers.hpp"
+#include "detray/definitions/track_parametrization.hpp"
+#include "detray/propagator/base_actor.hpp"
+#include "detray/tracks/detail/track_helper.hpp"
+
+namespace detray {
+
+template <typename transform3_t>
+struct parameter_resetter : actor {
+
+    struct state {};
+
+    struct kernel {
+
+        /// @name Type definitions for the struct
+        /// @{
+
+        // Matrix actor
+        using matrix_operator = typename transform3_t::matrix_actor;
+
+        /// @}
+
+        using output_type = bool;
+
+        template <typename mask_group_t, typename index_t,
+                  typename transform_store_t, typename surface_t,
+                  typename stepper_state_t>
+        DETRAY_HOST_DEVICE inline output_type operator()(
+            const mask_group_t& mask_group, const index_t& /*index*/,
+            transform_store_t& trf_store, surface_t& surface,
+            stepper_state_t& stepping) const {
+
+            const auto& trf3 = trf_store[surface.transform()];
+
+            // Note: How is it possible with "range"???
+            const auto& mask = mask_group[surface.mask_range()];
+
+            auto local_coordinate = mask.local_frame();
+
+            // Reset the free vector
+            stepping().set_vector(local_coordinate.bound_to_free_vector(
+                trf3, mask, stepping._bound_params.vector()));
+
+            // Reset the path length
+            stepping._s = 0;
+
+            // Reset jacobian coordinate transformation at the current surface
+            stepping._jac_to_global = local_coordinate.bound_to_free_jacobian(
+                trf3, mask, stepping._bound_params.vector());
+
+            // Reset jacobian transport to identity matrix
+            matrix_operator().set_identity(stepping._jac_transport);
+
+            return true;
+        }
+    };
+
+    template <typename propagator_state_t>
+    DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
+                                       propagator_state_t& propagation) const {
+
+        auto& navigation = propagation._navigation;
+        auto& stepping = propagation._stepping;
+
+        // Do covariance transport when the track is on surface
+        if (navigation.is_on_module()) {
+
+            auto det = navigation.detector();
+            const auto& trf_store = det->transform_store();
+            const auto& mask_store = det->mask_store();
+
+            // Intersection
+            const auto& is = navigation.current();
+
+            // Surface
+            const auto& surface = det->surface_by_index(is->index);
+
+            // Set surface link
+            stepping._bound_params.set_surface_link(is->index);
+
+            mask_store.template call<kernel>(surface.mask(), trf_store, surface,
+                                             stepping);
+        }
+    }
+};
+
+}  // namespace detray
diff --git a/core/include/detray/propagator/actors/bound_to_bound_updater.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp
similarity index 78%
rename from core/include/detray/propagator/actors/bound_to_bound_updater.hpp
rename to core/include/detray/propagator/actors/parameter_transporter.hpp
index bc215e090..46c96298d 100644
--- a/core/include/detray/propagator/actors/bound_to_bound_updater.hpp
+++ b/core/include/detray/propagator/actors/parameter_transporter.hpp
@@ -16,7 +16,7 @@
 namespace detray {
 
 template <typename transform3_t>
-struct bound_to_bound_updater : actor {
+struct parameter_transporter : actor {
 
     struct state {};
 
@@ -55,11 +55,11 @@ struct bound_to_bound_updater : actor {
 
         using output_type = bool;
 
-        template <typename mask_group_t, typename surface_t,
+        template <typename mask_group_t, typename index_t, typename surface_t,
                   typename propagator_state_t>
         DETRAY_HOST_DEVICE inline output_type operator()(
-            const mask_group_t& mask_group, const surface_t& surface,
-            propagator_state_t& propagation) const {
+            const mask_group_t& mask_group, const index_t& /*index*/,
+            const surface_t& surface, propagator_state_t& propagation) const {
 
             // Stepper and Navigator states
             auto& navigation = propagation._navigation;
@@ -69,15 +69,12 @@ struct bound_to_bound_updater : actor {
             const auto& det = navigation.detector();
             const auto& transform_store = det->transform_store();
 
-            // Intersection
-            const auto& is = navigation.current();
-
             // Transform
             const auto& trf3 = transform_store[surface.transform()];
 
             // Mask
-            const auto& mask = mask_group[is->mask_index];
-            auto local_coordinate = mask.local();
+            const auto& mask = mask_group[surface.mask_range()];
+            auto local_coordinate = mask.local_frame();
 
             // Free vector
             const auto& free_vec = stepping().vector();
@@ -101,7 +98,8 @@ struct bound_to_bound_updater : actor {
             const bound_to_free_matrix& bound_to_free_jacobian =
                 stepping._jac_to_global;
 
-            // Acts version
+            // @note: (Beomki) I really don't understand why the identity matrix
+            // should be added here but it makes result better :/
             const free_matrix correction_term =
                 matrix_operator()
                     .template identity<e_free_size, e_free_size>() +
@@ -111,20 +109,6 @@ struct bound_to_bound_updater : actor {
                 free_to_bound_jacobian * correction_term *
                 free_transport_jacobian * bound_to_free_jacobian;
 
-            /*
-            const bound_matrix full_jacobian =
-                free_to_bound_jacobian *
-                (path_correction + free_transport_jacobian) *
-                bound_to_free_jacobian;
-            */
-
-            // No path correction
-            /*
-            const bound_matrix full_jacobian = free_to_bound_jacobian *
-                                               free_transport_jacobian *
-                                               bound_to_free_jacobian;
-            */
-
             const bound_matrix new_cov =
                 full_jacobian * stepping._bound_params.covariance() *
                 matrix_operator().transpose(full_jacobian);
@@ -139,24 +123,22 @@ struct bound_to_bound_updater : actor {
     template <typename propagator_state_t>
     DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
                                        propagator_state_t& propagation) const {
-
         auto& navigation = propagation._navigation;
 
         // Do covariance transport when the track is on surface
         if (navigation.is_on_module()) {
 
-            const auto& det = navigation.detector();
-            const auto& surface_container = det->surfaces();
+            auto det = navigation.detector();
             const auto& mask_store = det->mask_store();
 
             // Intersection
             const auto& is = navigation.current();
 
             // Surface
-            const auto& surface = surface_container[is->index];
+            const auto& surface = det->surface_by_index(is->index);
 
-            mask_store.template execute<kernel>(surface.mask_type(), surface,
-                                                propagation);
+            mask_store.template call<kernel>(surface.mask(), surface,
+                                             propagation);
         }
     }
 };
diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp
index 246f95a47..4bca31e55 100644
--- a/core/include/detray/propagator/base_stepper.hpp
+++ b/core/include/detray/propagator/base_stepper.hpp
@@ -10,12 +10,21 @@
 // Project include(s).
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/definitions/units.hpp"
+#include "detray/propagator/actors/parameter_resetter.hpp"
 #include "detray/propagator/constrained_step.hpp"
-#include "detray/propagator/detail/covariance_engine.hpp"
 #include "detray/tracks/tracks.hpp"
 
 namespace detray {
 
+namespace stepping {
+
+enum class id {
+    e_linear = 0,
+    e_rk = 1,
+};
+
+}  // namespace stepping
+
 /// Base stepper implementation
 template <typename transform3_t, typename constraint_t, typename policy_t>
 class base_stepper {
@@ -25,7 +34,6 @@ class base_stepper {
     using free_track_parameters_type = free_track_parameters<transform3_t>;
     using bound_track_parameters_type = bound_track_parameters<transform3_t>;
     using matrix_operator = typename transform3_t::matrix_actor;
-    using covariance_engine = detail::covariance_engine<transform3_t>;
     using track_helper = detail::track_helper<matrix_operator>;
 
     using size_type = typename transform3_type::size_type;
@@ -58,25 +66,20 @@ class base_stepper {
         state(const free_track_parameters_type &t) : _track(t) {}
 
         /// Sets track parameters from bound track parameter.
-        DETRAY_HOST_DEVICE
-        state(const bound_track_parameters_type & /*bound_params*/,
-              const transform3_type & /*trf3*/) {
-
-            // @TODO: Rewrite with kernel
-
-            /*
-            // Set the free vector
-            _track.set_vector(vector_engine().bound_to_free_vector(
-                trf3, bound_params.vector()));
-
-            // Set the bound covariance
-            _bound_covariance = bound_params.covariance();
-
-            // Reset the jacobians
-            covariance_engine().reinitialize_jacobians(
-                trf3, bound_params.vector(), _jac_to_global, _jac_transport,
-                _derivative);
-            */
+        template <typename detector_t>
+        DETRAY_HOST_DEVICE state(
+            const bound_track_parameters_type &bound_params,
+            const detector_t &det)
+            : _bound_params(bound_params) {
+
+            const auto &trf_store = det.transform_store();
+            const auto &mask_store = det.mask_store();
+            const auto &surface =
+                det.surface_by_index(bound_params.surface_link());
+
+            mask_store.template call<
+                typename parameter_resetter<transform3_t>::kernel>(
+                surface.mask(), trf_store, surface, *this);
         }
 
         /// free track parameter
@@ -86,17 +89,12 @@ class base_stepper {
         free_matrix _jac_transport =
             matrix_operator().template identity<e_free_size, e_free_size>();
 
-        /// The free parameter derivative defined at destination surface
-        free_vector _derivative =
-            matrix_operator().template zero<e_free_size, 1>();
-
         /// bound-to-free jacobian from departure surface
         bound_to_free_matrix _jac_to_global =
             matrix_operator().template zero<e_free_size, e_bound_size>();
 
         /// bound covariance
-        bound_matrix _bound_covariance =
-            matrix_operator().template zero<e_bound_size, e_bound_size>();
+        bound_track_parameters_type _bound_params;
 
         /// @returns track parameters - const access
         DETRAY_HOST_DEVICE
@@ -117,6 +115,10 @@ class base_stepper {
         /// Track path length
         scalar _path_length{0.};
 
+        /// Track path length from the last surface. It will be reset to 0 when
+        /// the track reaches a new surface
+        scalar _s{0.};
+
         /// Current step size
         scalar _step_size{0.};
 
diff --git a/core/include/detray/propagator/detail/covariance_engine.hpp b/core/include/detray/propagator/detail/covariance_engine.hpp
deleted file mode 100644
index 8bae31fbe..000000000
--- a/core/include/detray/propagator/detail/covariance_engine.hpp
+++ /dev/null
@@ -1,101 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-// Detray include(s)
-#include "detray/definitions/qualifiers.hpp"
-#include "detray/propagator/detail/jacobian_engine.hpp"
-
-namespace detray {
-
-namespace detail {
-
-// Covariance engine agnostic to surface type
-template <typename transform3_t>
-struct covariance_engine {
-
-    /// Transformation matching this struct
-    using transform3_type = transform3_t;
-    using jacobian_engine = detail::jacobian_engine<transform3_type>;
-    using track_helper = typename jacobian_engine::track_helper;
-    using matrix_operator = typename jacobian_engine::matrix_actor;
-    using size_type = typename transform3_type::size_type;
-    /// 2D matrix type
-    template <size_type ROWS, size_type COLS>
-    using matrix_type =
-        typename matrix_operator::template matrix_type<ROWS, COLS>;
-    /// Shorthand vector types related to track parameters.
-    using bound_vector = matrix_type<e_bound_size, 1>;
-    using free_vector = matrix_type<e_free_size, 1>;
-    /// Shorthand matrix types related to track parameters.
-    using bound_matrix = matrix_type<e_bound_size, e_bound_size>;
-    using free_matrix = matrix_type<e_free_size, e_free_size>;
-    /// Mapping matrix
-    using bound_to_free_matrix = matrix_type<e_free_size, e_bound_size>;
-    using free_to_bound_matrix = matrix_type<e_bound_size, e_free_size>;
-    using free_to_path_matrix = matrix_type<1, e_free_size>;
-
-    /** Function to get the full jacobian for surface-to-surface propagation
-     *
-     * @param trf3 is the transform matrix of the destination surface
-     * @param free_vec is the input free vector at the destination surface
-     * @param bound_to_free_jacobian is the coordinate transform jacobian
-     * at departure surface.
-     * @param free_transport_jacobian is the transport jacobian
-     * @param free_to_path_derivative is the derivative of free parameter w.r.t
-     * path
-     * @returns bound to bound jacobian
-     */
-    DETRAY_HOST_DEVICE inline bound_matrix bound_to_bound_jacobian(
-        const transform3_type& trf3, const free_vector& free_vec,
-        const bound_to_free_matrix& bound_to_free_jacobian,
-        const free_matrix& free_transport_jacobian,
-        const free_vector& free_to_path_derivative) const;
-
-    /** Function to update the covariance matrix for surface-to-surface
-     * propagation
-     *
-     * @param trf3 is the transform matrix of the destination surface
-     * @param bound_covariance is the covariance at the departure surface, which
-     * is to be updated for the destination surface
-     * @param free_vec is the input free vector at the destination surface
-     * @param bound_to_free_jacobian is the coordinate transform jacobian
-     * at the departure surface.
-     * @param free_transport_jacobian is the transport jacobian
-     * @param free_to_path_derivative is the derivative of free parameter w.r.t
-     * path at the destination surface
-     */
-    DETRAY_HOST_DEVICE inline void bound_to_bound_covariance_update(
-        const transform3_type& trf3, bound_matrix& bound_covariance,
-        const free_vector& free_vec,
-        const bound_to_free_matrix& bound_to_free_jacobian,
-        const free_matrix& free_transport_jacobian,
-        const free_vector& free_to_path_derivative) const;
-
-    /** Function to reinitialize the jacobian after updating covariance
-     *
-     * @param trf3 is the transform matrix of the destination surface
-     * @param bound_vec is the input bound vector at the destination surface
-     * @param bound_to_free_jacobian is the coordinate transform jacobian
-     * at the destination surface
-     * @param free_transport_jacobian is the transport jacobian
-     * @param free_to_path_derivative is the derivative of free parameter w.r.t
-     * path at the destination surface
-     */
-    DETRAY_HOST_DEVICE inline void reinitialize_jacobians(
-        const transform3_type& trf3, const bound_vector& bound_vec,
-        bound_to_free_matrix& bound_to_free_jacobian,
-        free_matrix& free_transport_jacobian,
-        free_vector& free_to_path_derivative) const;
-};
-
-}  // namespace detail
-
-}  // namespace detray
-
-#include "detray/propagator/detail/covariance_engine.ipp"
\ No newline at end of file
diff --git a/core/include/detray/propagator/detail/covariance_engine.ipp b/core/include/detray/propagator/detail/covariance_engine.ipp
deleted file mode 100644
index 5ae5cd6d4..000000000
--- a/core/include/detray/propagator/detail/covariance_engine.ipp
+++ /dev/null
@@ -1,71 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-template <typename transform3_t>
-typename detray::detail::template covariance_engine<transform3_t>::bound_matrix
-detray::detail::covariance_engine<transform3_t>::bound_to_bound_jacobian(
-    const transform3_type& trf3, const free_vector& free_vec,
-    const bound_to_free_matrix& bound_to_free_jacobian,
-    const free_matrix& free_transport_jacobian,
-    const free_vector& free_to_path_derivative) const {
-
-    // Calculate the path correction term
-    free_to_path_matrix path_correction_term =
-        jacobian_engine().free_to_path_correction(trf3, free_vec);
-
-    // Calculate the free-to-bound coordinate transform jacobian
-    free_to_bound_matrix free_to_bound_jacobian =
-        jacobian_engine().free_to_bound_coordinate(trf3, free_vec);
-
-    // Calculate the full jacobian from the local/bound parameters at the start
-    // surface to local/bound parameters at the final surface
-    // @note jac(locA->locB) = jac(gloB->locB)*(1+
-    // pathCorrectionFactor(gloB))*jacTransport(gloA->gloB) *jac(locA->gloA)
-    const auto I =
-        matrix_operator().template identity<e_free_size, e_free_size>();
-
-    return free_to_bound_jacobian *
-           (I + free_to_path_derivative * path_correction_term) *
-           free_transport_jacobian * bound_to_free_jacobian;
-}
-
-template <typename transform3_t>
-void detray::detail::covariance_engine<transform3_t>::
-    bound_to_bound_covariance_update(
-        const transform3_type& trf3, bound_matrix& bound_covariance,
-        const free_vector& free_vec,
-        const bound_to_free_matrix& bound_to_free_jacobian,
-        const free_matrix& free_transport_jacobian,
-        const free_vector& free_to_path_derivative) const {
-
-    // Get the full jacobian
-    bound_matrix full_jacobian = this->bound_to_bound_jacobian(
-        trf3, free_vec, bound_to_free_jacobian, free_transport_jacobian,
-        free_to_path_derivative);
-
-    // Update the bound covariance
-    bound_covariance = full_jacobian * bound_covariance *
-                       matrix_operator().transpose(full_jacobian);
-}
-
-template <typename transform3_t>
-void detray::detail::covariance_engine<transform3_t>::reinitialize_jacobians(
-    const transform3_type& trf3, const bound_vector& bound_vec,
-    bound_to_free_matrix& bound_to_free_jacobian,
-    free_matrix& free_transport_jacobian,
-    free_vector& free_to_path_derivative) const {
-
-    // Reset jacobian coordinate transformation at the current surface
-    bound_to_free_jacobian =
-        jacobian_engine().bound_to_free_coordinate(trf3, bound_vec);
-
-    // Reset jacobian transport to identity matrix
-    matrix_operator().set_identity(free_transport_jacobian);
-
-    // Reset derivative of position and direction to zero matrix
-    matrix_operator().set_zero(free_to_path_derivative);
-}
\ No newline at end of file
diff --git a/core/include/detray/propagator/detail/jacobian_engine.hpp b/core/include/detray/propagator/detail/jacobian_engine.hpp
deleted file mode 100644
index 52f5c8317..000000000
--- a/core/include/detray/propagator/detail/jacobian_engine.hpp
+++ /dev/null
@@ -1,80 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-// Detray include(s)
-#include "detray/definitions/qualifiers.hpp"
-#include "detray/tracks/detail/track_helper.hpp"
-
-namespace detray {
-
-namespace detail {
-
-// Jacobian engine assuming "Planar Surface with Cartesian coordinate"
-// NOTE: We may inherit jacobian engine to address different
-// types of surface, if required
-template <typename transform3_t>
-struct jacobian_engine {
-
-    /// Transformation matching this struct
-    using transform3_type = transform3_t;
-    /// scalar_type
-    using scalar_type = typename transform3_type::scalar_type;
-    /// Vector in 3D space
-    using vector3 = typename transform3_type::vector3;
-    /// Matrix operator
-    using matrix_operator = typename transform3_type::matrix_actor;
-    // Track helper
-    using track_helper = detail::track_helper<matrix_operator>;
-    using size_type = typename transform3_type::size_type;
-    /// 2D matrix type
-    template <size_type ROWS, size_type COLS>
-    using matrix_type =
-        typename matrix_operator::template matrix_type<ROWS, COLS>;
-    /// Shorthand vector types related to track parameters.
-    using bound_vector = matrix_type<e_bound_size, 1>;
-    using free_vector = matrix_type<e_free_size, 1>;
-    /// Mapping matrix
-    using bound_to_free_matrix = matrix_type<e_free_size, e_bound_size>;
-    using free_to_bound_matrix = matrix_type<e_bound_size, e_free_size>;
-    using free_to_path_matrix = matrix_type<1, e_free_size>;
-
-    /** Function to get the jacobian for bound to free coordinate transform
-     *
-     * @param trf3 is the transform matrix
-     * @param bound_vec is the input bound vector
-     * @returns bound to free jacobian
-     */
-    DETRAY_HOST_DEVICE inline bound_to_free_matrix bound_to_free_coordinate(
-        const transform3_type& trf3, const bound_vector& bound_vec) const;
-
-    /** Function to get the jacobian for free to bound coordinate transform
-     *
-     * @param trf3 is the transform matrix
-     * @param free_vec is the input free vector
-     * @returns free to bound jacobian
-     */
-    DETRAY_HOST_DEVICE inline free_to_bound_matrix free_to_bound_coordinate(
-        const transform3_type& trf3, const free_vector& free_vec) const;
-
-    /** Function to calculate the path correction term for transport jacobian,
-     * which is caused by the geometry constraint.
-     *
-     * @param trf3 is the transform matrix
-     * @param free_vec is the input free vector
-     * @returns free to path matrix
-     */
-    DETRAY_HOST_DEVICE inline free_to_path_matrix free_to_path_correction(
-        const transform3_type& trf3, const free_vector& free_vec) const;
-};
-
-}  // namespace detail
-
-}  // namespace detray
-
-#include "detray/propagator/detail/jacobian_engine.ipp"
\ No newline at end of file
diff --git a/core/include/detray/propagator/detail/jacobian_engine.ipp b/core/include/detray/propagator/detail/jacobian_engine.ipp
deleted file mode 100644
index a9a21176f..000000000
--- a/core/include/detray/propagator/detail/jacobian_engine.ipp
+++ /dev/null
@@ -1,129 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-template <typename transform3_t>
-typename detray::detail::template jacobian_engine<
-    transform3_t>::bound_to_free_matrix
-detray::detail::jacobian_engine<transform3_t>::bound_to_free_coordinate(
-    const transform3_type& trf3, const bound_vector& bound_vec) const {
-
-    // Declare jacobian for bound to free coordinate transform
-    bound_to_free_matrix jac_to_global =
-        matrix_operator().template zero<e_free_size, e_bound_size>();
-
-    // Get trigonometric values
-    const scalar_type theta =
-        matrix_operator().element(bound_vec, e_bound_theta, 0);
-    const scalar_type phi =
-        matrix_operator().element(bound_vec, e_bound_phi, 0);
-    const scalar_type cos_theta = std::cos(theta);
-    const scalar_type sin_theta = std::sin(theta);
-    const scalar_type cos_phi = std::cos(phi);
-    const scalar_type sin_phi = std::sin(phi);
-
-    // Set d(x,y,z)/d(loc0,loc1)
-    const matrix_type<3, 2> bound_to_free_rotation =
-        matrix_operator().template block<3, 2>(trf3.matrix(), 0, 0);
-
-    matrix_operator().template set_block(jac_to_global, bound_to_free_rotation,
-                                         e_free_pos0, e_bound_loc0);
-
-    matrix_operator().element(jac_to_global, e_free_time, e_bound_time) = 1;
-    matrix_operator().element(jac_to_global, e_free_dir0, e_bound_phi) =
-        -1 * sin_theta * sin_phi;
-    matrix_operator().element(jac_to_global, e_free_dir0, e_bound_theta) =
-        cos_theta * cos_phi;
-    matrix_operator().element(jac_to_global, e_free_dir1, e_bound_phi) =
-        sin_theta * cos_phi;
-    matrix_operator().element(jac_to_global, e_free_dir1, e_bound_theta) =
-        cos_theta * sin_phi;
-    matrix_operator().element(jac_to_global, e_free_dir2, e_bound_theta) =
-        -1 * sin_theta;
-    matrix_operator().element(jac_to_global, e_free_qoverp, e_bound_qoverp) = 1;
-
-    return jac_to_global;
-}
-
-template <typename transform3_t>
-typename detray::detail::template jacobian_engine<
-    transform3_t>::free_to_bound_matrix
-detray::detail::jacobian_engine<transform3_t>::free_to_bound_coordinate(
-    const transform3_type& trf3, const free_vector& free_vec) const {
-
-    // Declare jacobian for free to bound coordinate transform
-    free_to_bound_matrix jac_to_local =
-        matrix_operator().template zero<e_bound_size, e_free_size>();
-
-    // Free direction
-    const vector3 dir = track_helper().dir(free_vec);
-
-    // Get trigonometric values
-    const scalar_type theta = getter::theta(dir);
-    const scalar_type phi = getter::phi(dir);
-    const scalar_type cos_theta = std::cos(theta);
-    const scalar_type sin_theta = std::sin(theta);
-    const scalar_type inv_sin_theta = 1. / sin_theta;
-    const scalar_type cos_phi = std::cos(phi);
-    const scalar_type sin_phi = std::sin(phi);
-
-    // Set d(loc0,loc1)/d(x,y,z)
-    const matrix_type<2, 3> free_to_bound_rotation =
-        matrix_operator().template block<2, 3>(trf3.matrix_inverse(), 0, 0);
-
-    matrix_operator().template set_block(jac_to_local, free_to_bound_rotation,
-                                         e_bound_loc0, e_free_pos0);
-
-    // Set d(Free time)/d(Bound time)
-    matrix_operator().element(jac_to_local, e_bound_time, e_free_time) = 1;
-
-    // Set d(phi, theta)/d(free dir)
-    matrix_operator().element(jac_to_local, e_bound_phi, e_free_dir0) =
-        -1. * sin_phi * inv_sin_theta;
-    matrix_operator().element(jac_to_local, e_bound_phi, e_free_dir1) =
-        cos_phi * inv_sin_theta;
-    matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir0) =
-        cos_phi * cos_theta;
-    matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir1) =
-        sin_phi * cos_theta;
-    matrix_operator().element(jac_to_local, e_bound_theta, e_free_dir2) =
-        -1 * sin_theta;
-
-    // Set d(Free Qop)/d(Bound Qop)
-    matrix_operator().element(jac_to_local, e_bound_qoverp, e_free_qoverp) = 1;
-
-    return jac_to_local;
-}
-
-template <typename transform3_t>
-typename detray::detail::template jacobian_engine<
-    transform3_t>::free_to_path_matrix
-detray::detail::jacobian_engine<transform3_t>::free_to_path_correction(
-    const transform3_type& trf3, const free_vector& free_vec) const {
-
-    // Declare free to path correction
-    free_to_path_matrix free_to_path =
-        matrix_operator().template zero<1, e_free_size>();
-
-    // Free direction
-    const vector3 dir = track_helper().dir(free_vec);
-
-    // The measurement frame z axis
-    const matrix_type<3, 1> ref_z_axis =
-        matrix_operator().template block<3, 1>(trf3.matrix(), 0, 2);
-
-    // cosine angle between momentum direction and the measurement frame z axis
-    const scalar_type dz = vector::dot(ref_z_axis, dir);
-
-    // Correction term
-    const matrix_type<1, 3> correction_term =
-        -1. / dz * matrix_operator().transpose(ref_z_axis);
-
-    matrix_operator().template set_block<1, 3>(free_to_path, correction_term, 0,
-                                               e_free_pos0);
-
-    return free_to_path;
-}
\ No newline at end of file
diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp
index 89db5575c..82297b6d9 100644
--- a/core/include/detray/propagator/line_stepper.hpp
+++ b/core/include/detray/propagator/line_stepper.hpp
@@ -30,11 +30,27 @@ class line_stepper final
     using policy_type = policy_t;
     using free_track_parameters_type =
         typename base_type::free_track_parameters_type;
+    using bound_track_parameters_type =
+        typename base_type::bound_track_parameters_type;
+    using transform3_type = transform3_t;
+    using matrix_operator = typename base_type::matrix_operator;
+    using size_type = typename matrix_operator::size_ty;
+    template <size_type ROWS, size_type COLS>
+    using matrix_type =
+        typename matrix_operator::template matrix_type<ROWS, COLS>;
 
     struct state : public base_type::state {
+        static constexpr const stepping::id id = stepping::id::e_linear;
+
         DETRAY_HOST_DEVICE
         state(const free_track_parameters_type &t) : base_type::state(t) {}
 
+        template <typename detector_t>
+        DETRAY_HOST_DEVICE state(
+            const bound_track_parameters_type &bound_params,
+            const detector_t &det)
+            : base_type::state(bound_params, det) {}
+
         /// Update the track state in a straight line.
         DETRAY_HOST_DEVICE
         inline void advance_track() {
@@ -43,6 +59,25 @@ class line_stepper final
 
             this->_path_length += this->_step_size;
         }
+
+        DETRAY_HOST_DEVICE
+        inline void advance_jacobian() {
+
+            // The step transport matrix in global coordinates
+            matrix_type<e_free_size, e_free_size> D =
+                matrix_operator().template identity<e_free_size, e_free_size>();
+
+            // d(x,y,z)/d(n_x,n_y,n_z)
+            matrix_type<3, 3> dxdn =
+                this->_step_size * matrix_operator().template identity<3, 3>();
+            matrix_operator().template set_block<3, 3>(D, dxdn, e_free_pos0,
+                                                       e_free_dir0);
+
+            /// NOTE: Let's skip the element for d(time)/d(qoverp) for the
+            /// moment..
+
+            this->_jac_transport = D * this->_jac_transport;
+        }
     };
 
     /// Take a step, regulared by a constrained step
@@ -77,6 +112,9 @@ class line_stepper final
         // Update track state
         stepping.advance_track();
 
+        // Advance jacobian transport
+        stepping.advance_jacobian();
+
         // Call navigation update policy
         policy_t{}(stepping.policy_state(), propagation);
 
diff --git a/core/include/detray/propagator/navigator.hpp b/core/include/detray/propagator/navigator.hpp
index 902bc64bd..7206b2cd2 100644
--- a/core/include/detray/propagator/navigator.hpp
+++ b/core/include/detray/propagator/navigator.hpp
@@ -116,6 +116,7 @@ class navigator {
 
         public:
         using detector_type = navigator::detector_type;
+        using intersection_t = typename navigator::intersection_type;
 
         /// Default constructor
         state() = default;
@@ -141,6 +142,30 @@ class navigator {
         DETRAY_HOST_DEVICE
         scalar operator()() const { return _next->path; }
 
+        DETRAY_HOST_DEVICE
+        inline void set_unknown() { _status = navigation::status::e_unknown; }
+
+        DETRAY_HOST_DEVICE
+        inline void set_on_module() {
+            _status = navigation::status::e_on_module;
+        }
+
+        /// Updates the iterator position of the last valid candidate
+        DETRAY_HOST_DEVICE
+        inline void set_next(candidate_itr_t &&new_next) {
+            _next = std::move(new_next);
+        }
+
+        /// @returns next object that we want to reach (current target)
+        DETRAY_HOST_DEVICE
+        inline auto next() -> candidate_itr_t & { return _next; }
+
+        /// @returns next object that we want to reach (current target) - const
+        DETRAY_HOST_DEVICE
+        inline auto next() const -> const const_candidate_itr_t & {
+            return _next;
+        }
+
         /// @returns currently cached candidates - const
         DETRAY_HOST_DEVICE
         inline auto candidates() const
@@ -167,11 +192,13 @@ class navigator {
             return _next - 1;
         }
 
+        /*
         /// @returns next object that we want to reach (current target) - const
         DETRAY_HOST_DEVICE
         inline auto next() const -> const const_candidate_itr_t & {
             return _next;
         }
+        */
 
         /// @returns last valid candidate (by position in the cache) - const
         DETRAY_HOST_DEVICE
@@ -327,15 +354,18 @@ class navigator {
                    candidate.path >= track.overstep_tolerance();
         }
 
+        /*
         /// @returns next object that we want to reach (current target)
         DETRAY_HOST_DEVICE
         inline auto next() -> candidate_itr_t & { return _next; }
-
+        */
+        /*
         /// Updates the iterator position of the last valid candidate
         DETRAY_HOST_DEVICE
         inline void set_next(candidate_itr_t &&new_next) {
             _next = std::move(new_next);
         }
+        */
 
         /// Updates the iterator position of the last valid candidate
         DETRAY_HOST_DEVICE
diff --git a/core/include/detray/propagator/propagator.hpp b/core/include/detray/propagator/propagator.hpp
index 591857767..d0114b9fa 100644
--- a/core/include/detray/propagator/propagator.hpp
+++ b/core/include/detray/propagator/propagator.hpp
@@ -53,6 +53,7 @@ struct propagator {
     struct state {
 
         using detector_type = typename navigator_t::detector_type;
+        using navigator_state_type = typename navigator_t::state;
 
         /// Construct the propagation state.
         ///
@@ -78,13 +79,22 @@ struct propagator {
               _actor_states(actor_states) {}
 
         /// Construct the propagation state with bound parameter
+        DETRAY_HOST_DEVICE state(
+            const bound_track_parameters_type &param, const detector_type &det,
+            typename actor_chain_t::state actor_states = {},
+            vector_type<line_plane_intersection> &&candidates = {})
+            : _stepping(param, det),
+              _navigation(det, std::move(candidates)),
+              _actor_states(actor_states) {}
+
+        /// Construct the propagation state with bound parameter
+        template <typename field_t>
         DETRAY_HOST_DEVICE state(
             const bound_track_parameters_type &param,
-            const typename stepper_t::transform3_type &trf3,
-            const detector_type &det,
+            const field_t &magnetic_field, const detector_type &det,
             typename actor_chain_t::state actor_states = {},
             vector_type<line_plane_intersection> &&candidates = {})
-            : _stepping(param, trf3),
+            : _stepping(param, magnetic_field, det),
               _navigation(det, std::move(candidates)),
               _actor_states(actor_states) {}
 
diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp
index b0f571a64..843aa072d 100644
--- a/core/include/detray/propagator/rk_stepper.hpp
+++ b/core/include/detray/propagator/rk_stepper.hpp
@@ -11,7 +11,6 @@
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/definitions/units.hpp"
 #include "detray/propagator/base_stepper.hpp"
-#include "detray/propagator/detail/covariance_engine.hpp"
 #include "detray/propagator/navigation_policies.hpp"
 #include "detray/tracks/tracks.hpp"
 #include "detray/utils/column_wise_operator.hpp"
@@ -38,7 +37,6 @@ class rk_stepper final
     using vector2 = typename transform3_type::point2;
     using vector3 = typename transform3_type::vector3;
     using matrix_operator = typename base_type::matrix_operator;
-    using covariance_engine = typename base_type::covariance_engine;
     using column_wise_op = column_wise_operator<matrix_operator>;
 
     using free_track_parameters_type =
@@ -51,19 +49,18 @@ class rk_stepper final
 
     struct state : public base_type::state {
 
-        using field_type = magnetic_field_t;
+        static constexpr const stepping::id id = stepping::id::e_rk;
 
         DETRAY_HOST_DEVICE
         state(const free_track_parameters_type& t,
               const magnetic_field_t& mag_field)
             : base_type::state(t), _magnetic_field(mag_field) {}
 
-        DETRAY_HOST_DEVICE
-        state(const bound_track_parameters_type& bound_params,
-              const transform3_type& trf3, const field_type& mag_field)
-            : base_type::state(bound_params, trf3),
-              _magnetic_field(mag_field) {}
-
+        template <typename detector_t>
+        DETRAY_HOST_DEVICE state(
+            const bound_track_parameters_type& bound_params,
+            const magnetic_field_t& mag_field, const detector_t& det)
+            : base_type::state(bound_params, det), _magnetic_field(mag_field) {}
         /// error tolerance
         scalar _tolerance = 1e-4;
 
@@ -110,14 +107,6 @@ class rk_stepper final
      */
     template <typename propagation_state_t>
     DETRAY_HOST_DEVICE bool step(propagation_state_t& propagation);
-
-    /** Get the bound state at the surface
-     *
-     * @return returning the bound track parameter
-     */
-    template <typename propagation_state_t>
-    DETRAY_HOST_DEVICE bound_track_parameters_type bound_state(
-        propagation_state_t& /*propagation*/, const transform3_type& /*trf*/);
 };
 
 }  // namespace detray
diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp
index 73b83e5b5..88bfb0c56 100644
--- a/core/include/detray/propagator/rk_stepper.ipp
+++ b/core/include/detray/propagator/rk_stepper.ipp
@@ -11,21 +11,6 @@
 // System include(s)
 #include <cmath>
 
-template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t,
-          template <typename, std::size_t> class array_t>
-void detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                        array_t>::state::advance_derivative() {
-
-    // The derivative of position = direction
-    matrix_operator().set_block(this->_derivative, this->_track.dir(),
-                                e_free_pos0, 0);
-
-    // The derivative of direction = k4
-    matrix_operator().set_block(this->_derivative, this->_step_data.k4,
-                                e_free_dir0, 0);
-}
-
 template <typename magnetic_field_t, typename transform3_t,
           typename constraint_t, typename policy_t,
           template <typename, std::size_t> class array_t>
@@ -48,6 +33,7 @@ void detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
 
     // Update path length
     this->_path_length += h;
+    this->_s += h;
 }
 
 template <typename magnetic_field_t, typename transform3_t,
@@ -178,6 +164,7 @@ template <typename magnetic_field_t, typename transform3_t,
 template <typename propagation_state_t>
 bool detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
                         array_t>::step(propagation_state_t& propagation) {
+
     // Get stepper and navigator states
     state& stepping = propagation._stepping;
     auto& magnetic_field = stepping._magnetic_field;
@@ -286,9 +273,6 @@ bool detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
             stepping.constraints().template size<>(stepping.direction()));
     }
 
-    // Advance derivative of position and direction w.r.t path length
-    stepping.advance_derivative();
-
     // Advance track state
     stepping.advance_track();
 
@@ -300,43 +284,3 @@ bool detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
 
     return true;
 }
-
-template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t,
-          template <typename, std::size_t> class array_t>
-template <typename propagation_state_t>
-typename detray::template rk_stepper<magnetic_field_t, transform3_t,
-                                     constraint_t, policy_t,
-                                     array_t>::bound_track_parameters_type
-detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   array_t>::bound_state(propagation_state_t& /*propagation*/,
-                                         const transform3_t& /*trf3*/) {
-
-    return detray::bound_track_parameters<transform3_t>();
-
-    // @TODO: Use detector kernel
-    /*
-    auto& stepping = propagation._stepping;
-    auto& navigation = propagation._navigation;
-
-    // Get the free vector
-    const auto& free_vector = stepping().vector();
-
-    // Get the bound vector
-    const auto bound_vector =
-        vector_engine().free_to_bound_vector(trf3, free_vector);
-
-    // Get the bound covariance
-    covariance_engine().bound_to_bound_covariance_update(
-        trf3, stepping._bound_covariance, free_vector, stepping._jac_to_global,
-        stepping._jac_transport, stepping._derivative);
-
-    // Reset the jacobians in RK stepper state
-    covariance_engine().reinitialize_jacobians(
-        trf3, bound_vector, stepping._jac_to_global, stepping._jac_transport,
-        stepping._derivative);
-
-    return detray::bound_track_parameters(
-        navigation.current_object(), bound_vector, stepping._bound_covariance);
-    */
-}
diff --git a/core/include/detray/tracks/bound_track_parameters.hpp b/core/include/detray/tracks/bound_track_parameters.hpp
index ddc462593..f85762332 100644
--- a/core/include/detray/tracks/bound_track_parameters.hpp
+++ b/core/include/detray/tracks/bound_track_parameters.hpp
@@ -57,6 +57,9 @@ struct bound_track_parameters {
     DETRAY_HOST_DEVICE
     const std::size_t& surface_link() const { return m_surface_link; }
 
+    DETRAY_HOST_DEVICE
+    void set_surface_link(std::size_t link) { m_surface_link = link; }
+
     DETRAY_HOST_DEVICE
     const vector_type& vector() const { return m_vector; }
 
diff --git a/core/include/detray/tracks/detail/track_helper.hpp b/core/include/detray/tracks/detail/track_helper.hpp
index 0b84dcf66..4b618a9e1 100644
--- a/core/include/detray/tracks/detail/track_helper.hpp
+++ b/core/include/detray/tracks/detail/track_helper.hpp
@@ -10,7 +10,6 @@
 // Project include(s).
 #include "detray/definitions/qualifiers.hpp"
 #include "detray/definitions/track_parametrization.hpp"
-#include "detray/tracks/detail/trigonometrics.hpp"
 
 // System include(s).
 #include <cmath>
@@ -39,8 +38,6 @@ struct track_helper {
     using point3 = vector3;
     /// Point in 2D space
     using point2 = array_type<2>;
-    /// Trigonometrics
-    using trigonometrics = detail::trigonometrics<scalar_type>;
 
     /// Track vector types
     using bound_vector = matrix_type<e_bound_size, 1>;
@@ -90,29 +87,6 @@ struct track_helper {
         return {std::cos(phi) * sinTheta, std::sin(phi) * sinTheta,
                 std::cos(theta)};
     }
-
-    DETRAY_HOST_DEVICE
-    inline trigonometrics get_trigonometrics(
-        const bound_vector& bound_vec) const {
-        const scalar_type theta =
-            matrix_operator().element(bound_vec, e_bound_theta, 0);
-        const scalar_type phi =
-            matrix_operator().element(bound_vec, e_bound_phi, 0);
-
-        return {std::cos(theta), std::sin(theta), std::cos(phi), std::sin(phi)};
-    }
-
-    DETRAY_HOST_DEVICE
-    inline trigonometrics get_trigonometrics(
-        const free_vector& free_vec) const {
-
-        const vector3 dir = this->dir(free_vec);
-
-        const scalar_type theta = getter::theta(dir);
-        const scalar_type phi = getter::phi(dir);
-
-        return {std::cos(theta), std::sin(theta), std::cos(phi), std::sin(phi)};
-    }
 };
 
-}  // namespace detray::detail
+}  // namespace detray::detail
\ No newline at end of file
diff --git a/core/include/detray/tracks/detail/trigonometrics.hpp b/core/include/detray/tracks/detail/trigonometrics.hpp
deleted file mode 100644
index e5e7dcd7c..000000000
--- a/core/include/detray/tracks/detail/trigonometrics.hpp
+++ /dev/null
@@ -1,26 +0,0 @@
-/** Algebra plugins library, part of the ACTS project
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-// Project include(s).
-#include "detray/definitions/qualifiers.hpp"
-
-// System include(s).
-#include <cmath>
-
-namespace detray::detail {
-
-template <typename scalar_t>
-struct trigonometrics {
-    scalar_t cos_theta;
-    scalar_t sin_theta;
-    scalar_t cos_phi;
-    scalar_t sin_phi;
-};
-
-}  // namespace detray::detail
\ No newline at end of file
diff --git a/core/include/detray/tracks/free_track_parameters.hpp b/core/include/detray/tracks/free_track_parameters.hpp
index 55a2b817b..6dffb451b 100644
--- a/core/include/detray/tracks/free_track_parameters.hpp
+++ b/core/include/detray/tracks/free_track_parameters.hpp
@@ -110,8 +110,8 @@ struct free_track_parameters {
 
     DETRAY_HOST_DEVICE
     scalar_type charge() const {
-        return -1 ? matrix_operator().element(m_vector, e_free_qoverp, 0) < 0
-                  : 1.;
+        return matrix_operator().element(m_vector, e_free_qoverp, 0) < 0 ? -1.
+                                                                         : 1.;
     }
 
     DETRAY_HOST_DEVICE
diff --git a/core/include/detray/utils/invalid_values.hpp b/core/include/detray/utils/invalid_values.hpp
index 04fa39b18..cf9ea3c2f 100644
--- a/core/include/detray/utils/invalid_values.hpp
+++ b/core/include/detray/utils/invalid_values.hpp
@@ -7,6 +7,9 @@
 
 #pragma once
 
+// Project include(s).
+#include "detray/definitions/qualifiers.hpp"
+
 // System include(s)
 #include <limits>
 #include <type_traits>
diff --git a/tests/common/include/tests/common/coordinate_cartesian2.inl b/tests/common/include/tests/common/coordinate_cartesian2.inl
new file mode 100644
index 000000000..ed4db8212
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_cartesian2.inl
@@ -0,0 +1,104 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/cartesian2.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+using size_type = typename matrix_operator::size_ty;
+template <size_type ROWS, size_type COLS>
+using matrix_type = typename matrix_operator::template matrix_type<ROWS, COLS>;
+
+const scalar isclose = 1e-5;
+
+// This test cartesian2 coordinate
+TEST(coordinate, cartesian2) {
+
+    // Preparation work
+    const vector3 z = {0., 0., 1.};
+    const vector3 x = {1., 0., 0.};
+    const point3 t = {2., 3., 4.};
+    const transform3 trf(t, z, x);
+    const cartesian2<transform3> c2;
+    const point3 global1 = {4., 7., 4.};
+    const vector3 mom = {1., 2., 3.};
+    const vector3 d = vector::normalize(mom);
+    const scalar time = 0.1;
+    const scalar charge = -1.;
+    struct dummy_mask {
+    } mask;
+
+    // Global to local transformation
+    const point2 local = c2.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], 2., isclose);
+    ASSERT_NEAR(local[1], 4., isclose);
+
+    // Local to global transformation
+    const point3 global2 = c2.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+
+    // Free track parameter
+    const free_track_parameters<transform3> free_params(global1, time, mom,
+                                                        charge);
+    const auto free_vec1 = free_params.vector();
+
+    const auto bound_vec = c2.free_to_bound_vector(trf, free_vec1);
+    const auto free_vec2 = c2.bound_to_free_vector(trf, mask, bound_vec);
+
+    const matrix_operator m;
+
+    // Check if the bound vector is correct
+    ASSERT_NEAR(m.element(bound_vec, 0, 0), 2., isclose);
+    ASSERT_NEAR(m.element(bound_vec, 1, 0), 4., isclose);
+    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
+                isclose);  // atan(2)
+    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
+                isclose);  // atan(sqrt(5)/3)
+    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
+
+    // Check if the same free vector is obtained
+    for (int i = 0; i < 8; i++) {
+        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
+                    isclose);
+    }
+
+    // Normal vector
+    const vector3 n =
+        c2.normal(trf, mask, free_params.pos(), free_params.dir());
+    ASSERT_EQ(n, vector3({0., 0., 1.}));
+
+    // Test Jacobian transformation
+    const matrix_type<6, 6> J =
+        c2.free_to_bound_jacobian(trf, mask, free_vec1) *
+        c2.bound_to_free_jacobian(trf, mask, bound_vec);
+
+    for (std::size_t i = 0; i < 6; i++) {
+        for (std::size_t j = 0; j < 6; j++) {
+            if (i == j) {
+                EXPECT_NEAR(m.element(J, i, j), 1., isclose);
+            } else {
+                EXPECT_NEAR(m.element(J, i, j), 0., isclose);
+            }
+        }
+    }
+}
diff --git a/tests/common/include/tests/common/coordinate_cartesian3.inl b/tests/common/include/tests/common/coordinate_cartesian3.inl
new file mode 100644
index 000000000..b15fe1a52
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_cartesian3.inl
@@ -0,0 +1,54 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/cartesian3.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+
+const scalar isclose = 1e-5;
+
+// This test cartesian3 coordinate
+TEST(coordinate, cartesian3) {
+
+    // Preparation work
+    const vector3 z = {0., 0., 1.};
+    const vector3 x = {1., 0., 0.};
+    const point3 t = {2., 3., 4.};
+    const transform3 trf(t, z, x);
+    const cartesian3<transform3> c3;
+    const point3 global1 = {4., 7., 5.};
+    const vector3 mom = {1., 2., 3.};
+    const vector3 d = vector::normalize(mom);
+    struct dummy_mask {
+    } mask;
+
+    // Global to local transformation
+    const point3 local = c3.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], 2., isclose);
+    ASSERT_NEAR(local[1], 4., isclose);
+    ASSERT_NEAR(local[2], 1., isclose);
+
+    // Local to global transformation
+    const point3 global2 = c3.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+}
diff --git a/tests/common/include/tests/common/coordinate_cylindrical2.inl b/tests/common/include/tests/common/coordinate_cylindrical2.inl
new file mode 100644
index 000000000..027042cf3
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_cylindrical2.inl
@@ -0,0 +1,113 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/cylindrical2.hpp"
+#include "detray/masks/masks.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+// System include(s).
+#include <climits>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+using size_type = typename matrix_operator::size_ty;
+template <size_type ROWS, size_type COLS>
+using matrix_type = typename matrix_operator::template matrix_type<ROWS, COLS>;
+
+const scalar isclose = 1e-5;
+
+// This test cylindrical2 coordinate
+TEST(coordinate, cylindrical2) {
+
+    // Preparation work
+    const vector3 z = {0., 0., 1.};
+    const vector3 x = {1., 0., 0.};
+    const point3 t = {2., 3., 4.};
+    const transform3 trf(t, z, x);
+    const cylindrical2<transform3> c2;
+    // Global position on surface
+    const point3 global1 = {scalar{3.4142136}, scalar{4.4142136}, scalar{9.}};
+    const vector3 mom = {1., 2., 3.};
+    const vector3 d = vector::normalize(mom);
+    const scalar time = 0.1;
+    const scalar charge = -1.;
+
+    const scalar r = 2.;
+    const scalar hz = std::numeric_limits<scalar>::infinity();
+    mask<cylinder2D<>> mask{0UL, r, -hz, hz};
+
+    // Global to local transformation
+    const point2 local = c2.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], r * scalar{M_PI_4}, isclose);
+    ASSERT_NEAR(local[1], 5., isclose);
+
+    // Local to global transformation
+    const point3 global2 = c2.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+
+    // Free track parameter
+    const free_track_parameters<transform3> free_params(global1, time, mom,
+                                                        charge);
+    const auto free_vec1 = free_params.vector();
+
+    const auto bound_vec = c2.free_to_bound_vector(trf, free_vec1);
+    const auto free_vec2 = c2.bound_to_free_vector(trf, mask, bound_vec);
+
+    const matrix_operator m;
+
+    // Check if the bound vector is correct
+    ASSERT_NEAR(m.element(bound_vec, 0, 0), r * scalar{M_PI_4}, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 1, 0), 5, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
+                isclose);  // atan(2)
+    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
+                isclose);  // atan(sqrt(5)/3)
+    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
+
+    // Check if the same free vector is obtained
+    for (int i = 0; i < 8; i++) {
+        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
+                    isclose);
+    }
+
+    // Normal vector
+    const vector3 n =
+        c2.normal(trf, mask, free_params.pos(), free_params.dir());
+    ASSERT_NEAR(n[0], 1. / std::sqrt(2), isclose);
+    ASSERT_NEAR(n[1], 1. / std::sqrt(2), isclose);
+    ASSERT_NEAR(n[2], 0., isclose);
+
+    // Test Jacobian transformation
+    const matrix_type<6, 6> J =
+        c2.free_to_bound_jacobian(trf, mask, free_vec1) *
+        c2.bound_to_free_jacobian(trf, mask, bound_vec);
+
+    for (std::size_t i = 0; i < 6; i++) {
+        for (std::size_t j = 0; j < 6; j++) {
+            if (i == j) {
+                EXPECT_NEAR(m.element(J, i, j), 1., isclose);
+            } else {
+                EXPECT_NEAR(m.element(J, i, j), 0., isclose);
+            }
+        }
+    }
+}
diff --git a/tests/common/include/tests/common/coordinate_cylindrical3.inl b/tests/common/include/tests/common/coordinate_cylindrical3.inl
new file mode 100644
index 000000000..de3bf4d7d
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_cylindrical3.inl
@@ -0,0 +1,54 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/cylindrical3.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+
+const scalar isclose = 1e-5;
+
+// This test cylindrical3 coordinate
+TEST(coordinate, cylindrical3) {
+
+    // Preparation work
+    const vector3 z = {0., 0., 1.};
+    const vector3 x = {1., 0., 0.};
+    const point3 t = {2., 3., 4.};
+    const transform3 trf(t, z, x);
+    const cylindrical3<transform3> c3;
+    const point3 global1 = {scalar{3.4142136}, scalar{4.4142136}, scalar{9.}};
+    const vector3 mom = {1., 2., 3.};
+    const vector3 d = vector::normalize(mom);
+    struct dummy_mask {
+    } mask;
+
+    // Global to local transformation
+    const point3 local = c3.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], 2, isclose);
+    ASSERT_NEAR(local[1], M_PI_4, isclose);
+    ASSERT_NEAR(local[2], 5., isclose);
+
+    // Local to global transformation
+    const point3 global2 = c3.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+}
diff --git a/tests/common/include/tests/common/coordinate_line2.inl b/tests/common/include/tests/common/coordinate_line2.inl
new file mode 100644
index 000000000..bac23620f
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_line2.inl
@@ -0,0 +1,149 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/line2.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+using size_type = typename matrix_operator::size_ty;
+template <size_type ROWS, size_type COLS>
+using matrix_type = typename matrix_operator::template matrix_type<ROWS, COLS>;
+
+const scalar isclose = 1e-5;
+
+TEST(coordinate, line2_case1) {
+
+    // Preparation work
+    vector3 z = {1., 1., 1.};
+    z = vector::normalize(z);
+    vector3 x = {1., 0., -1.};
+    x = vector::normalize(x);
+    const point3 t = {0., 0., 0.};
+    const transform3 trf(t, z, x);
+    const line2<transform3> l2;
+    const point3 global1 = {1, 1.5, 0.5};
+    const vector3 mom = {0., 1., 1.};
+    const vector3 d = vector::normalize(mom);
+    const scalar time = 0.1;
+    const scalar charge = -1.;
+    struct dummy_mask {
+    } mask;
+
+    // Global to local transformation
+    const point2 local = l2.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], -1. / sqrt(2), isclose);
+    ASSERT_NEAR(local[1], sqrt(3), isclose);
+
+    // Local to global transformation
+    const point3 global2 = l2.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+
+    // Free track parameter
+    const free_track_parameters<transform3> free_params(global1, time, mom,
+                                                        charge);
+    const auto free_vec1 = free_params.vector();
+
+    const auto bound_vec = l2.free_to_bound_vector(trf, free_vec1);
+    const auto free_vec2 = l2.bound_to_free_vector(trf, mask, bound_vec);
+
+    const matrix_operator m;
+
+    // Check if the bound vector is correct
+    ASSERT_NEAR(m.element(bound_vec, 0, 0), -1. / sqrt(2), isclose);
+    ASSERT_NEAR(m.element(bound_vec, 1, 0), sqrt(3), isclose);
+    ASSERT_NEAR(m.element(bound_vec, 2, 0), M_PI_2, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 3, 0), M_PI_4, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1. / sqrt(2), isclose);
+    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
+
+    // Check if the same free vector is obtained
+    for (int i = 0; i < 8; i++) {
+        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
+                    isclose);
+    }
+
+    // Test Jacobian transformation
+    const matrix_type<6, 6> J =
+        l2.free_to_bound_jacobian(trf, mask, free_vec1) *
+        l2.bound_to_free_jacobian(trf, mask, bound_vec);
+
+    for (std::size_t i = 0; i < 6; i++) {
+        for (std::size_t j = 0; j < 6; j++) {
+            if (i == j) {
+                EXPECT_NEAR(m.element(J, i, j), 1., isclose);
+            } else {
+                EXPECT_NEAR(m.element(J, i, j), 0., isclose);
+            }
+        }
+    }
+}
+
+TEST(coordinate, line2_case2) {
+
+    // Preparation work
+    vector3 z = {1., 2., 3.};
+    z = vector::normalize(z);
+    vector3 x = {2., -4., 2.};
+    x = vector::normalize(x);
+    const point3 t = {0., 0., 0.};
+    const transform3 trf(t, z, x);
+    const line2<transform3> l2;
+    const point2 local1 = {1, 2};
+    const vector3 mom = {1., 6., -2.};
+    const vector3 d = vector::normalize(mom);
+    const scalar time = 0.1;
+    const scalar charge = -1.;
+    struct dummy_mask {
+    } mask;
+
+    // local to global transformation
+    const point3 global = l2.local_to_global(trf, mask, local1, d);
+
+    // global to local transform
+    const point2 local2 = l2.global_to_local(trf, global, d);
+
+    // Check if the same local position is obtained
+    ASSERT_NEAR(local1[0], local2[0], isclose);
+    ASSERT_NEAR(local1[1], local2[1], isclose);
+
+    // Free track parameter
+    const free_track_parameters<transform3> free_params(global, time, mom,
+                                                        charge);
+    const auto free_vec = free_params.vector();
+    const auto bound_vec = l2.free_to_bound_vector(trf, free_vec);
+
+    // Test Jacobian transformation
+    const matrix_type<6, 6> J = l2.free_to_bound_jacobian(trf, mask, free_vec) *
+                                l2.bound_to_free_jacobian(trf, mask, bound_vec);
+
+    const matrix_operator m;
+
+    for (std::size_t i = 0; i < 6; i++) {
+        for (std::size_t j = 0; j < 6; j++) {
+            if (i == j) {
+                EXPECT_NEAR(m.element(J, i, j), 1., isclose);
+            } else {
+                EXPECT_NEAR(m.element(J, i, j), 0., isclose);
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/tests/common/include/tests/common/coordinate_polar2.inl b/tests/common/include/tests/common/coordinate_polar2.inl
new file mode 100644
index 000000000..ead064cab
--- /dev/null
+++ b/tests/common/include/tests/common/coordinate_polar2.inl
@@ -0,0 +1,105 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2020-2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/coordinates/polar2.hpp"
+#include "detray/tracks/tracks.hpp"
+
+// GTest include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using point2 = __plugin::point2<scalar>;
+using point3 = __plugin::point3<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+using matrix_operator = typename transform3::matrix_actor;
+using size_type = typename matrix_operator::size_ty;
+template <size_type ROWS, size_type COLS>
+using matrix_type = typename matrix_operator::template matrix_type<ROWS, COLS>;
+
+const scalar isclose = 1e-5;
+
+// This test polar2 coordinate
+TEST(coordinate, polar2) {
+
+    // Preparation work
+    const vector3 z = {0., 0., 1.};
+    const vector3 x = {1., 0., 0.};
+    const point3 t = {2., 3., 4.};
+    const transform3 trf(t, z, x);
+    const polar2<transform3> p2;
+    const point3 global1 = {4., 7., 4.};
+    const vector3 mom = {1., 2., 3.};
+    const vector3 d = vector::normalize(mom);
+    const scalar time = 0.1;
+    const scalar charge = -1.;
+    struct dummy_mask {
+    } mask;
+
+    // Global to local transformation
+    const point2 local = p2.global_to_local(trf, global1, d);
+
+    // Check if the local position is correct
+    ASSERT_NEAR(local[0], std::sqrt(20.), isclose);
+    ASSERT_NEAR(local[1], atan2(4., 2.), isclose);
+
+    // Local to global transformation
+    const point3 global2 = p2.local_to_global(trf, mask, local, d);
+
+    // Check if the same global position is obtained
+    ASSERT_NEAR(global1[0], global2[0], isclose);
+    ASSERT_NEAR(global1[1], global2[1], isclose);
+    ASSERT_NEAR(global1[2], global2[2], isclose);
+
+    // Free track parameter
+    const free_track_parameters<transform3> free_params(global1, time, mom,
+                                                        charge);
+    const auto free_vec1 = free_params.vector();
+
+    const auto bound_vec = p2.free_to_bound_vector(trf, free_vec1);
+    const auto free_vec2 = p2.bound_to_free_vector(trf, mask, bound_vec);
+
+    const matrix_operator m;
+
+    // Check if the bound vector is correct
+    ASSERT_NEAR(m.element(bound_vec, 0, 0), std::sqrt(20.), isclose);
+    ASSERT_NEAR(m.element(bound_vec, 1, 0), atan2(4., 2.), isclose);
+    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
+                isclose);  // atan(2)
+    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
+                isclose);  // atan(sqrt(5)/3)
+    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
+    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
+
+    // Check if the same free vector is obtained
+    for (int i = 0; i < 8; i++) {
+        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
+                    isclose);
+    }
+
+    // Normal vector
+    const vector3 n =
+        p2.normal(trf, mask, free_params.pos(), free_params.dir());
+    ASSERT_EQ(n, vector3({0., 0., 1.}));
+
+    // Test Jacobian transformation
+    const matrix_type<6, 6> J =
+        p2.free_to_bound_jacobian(trf, mask, free_vec1) *
+        p2.bound_to_free_jacobian(trf, mask, bound_vec);
+
+    for (std::size_t i = 0; i < 6; i++) {
+        for (std::size_t j = 0; j < 6; j++) {
+
+            if (i == j) {
+                EXPECT_NEAR(m.element(J, i, j), 1., isclose);
+            } else {
+                EXPECT_NEAR(m.element(J, i, j), 0., isclose);
+            }
+        }
+    }
+}
\ No newline at end of file
diff --git a/tests/common/include/tests/common/coordinates.inl b/tests/common/include/tests/common/coordinates.inl
deleted file mode 100644
index ca8e91d37..000000000
--- a/tests/common/include/tests/common/coordinates.inl
+++ /dev/null
@@ -1,328 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2020-2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-// Project include(s).
-#include "detray/coordinates/coordinates.hpp"
-#include "detray/tracks/tracks.hpp"
-
-// GTest include(s).
-#include <gtest/gtest.h>
-
-using namespace detray;
-using point2 = __plugin::point2<scalar>;
-using point3 = __plugin::point3<scalar>;
-using vector3 = __plugin::vector3<scalar>;
-using transform3 = __plugin::transform3<detray::scalar>;
-using matrix_operator = typename transform3::matrix_actor;
-
-const scalar isclose = 1e-5;
-
-// This test cartesian2 coordinate
-TEST(test_host_basics, cartesian2) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const cartesian2<transform3> c2;
-    const point3 global1 = {4., 7., 4.};
-    const vector3 mom = {1., 2., 3.};
-    const vector3 d = vector::normalize(mom);
-    const scalar time = 0.1;
-    const scalar charge = -1.;
-    struct dummy_mask {
-    } mask;
-
-    // Global to local transformation
-    const point2 local = c2.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], 2., isclose);
-    ASSERT_NEAR(local[1], 4., isclose);
-
-    // Local to global transformation
-    const point3 global2 = c2.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-
-    // Free track parameter
-    const free_track_parameters<transform3> free_params(global1, time, mom,
-                                                        charge);
-    const auto free_vec1 = free_params.vector();
-
-    const auto bound_vec = c2.free_to_bound_vector(trf, free_vec1);
-    const auto free_vec2 = c2.bound_to_free_vector(trf, mask, bound_vec);
-
-    const matrix_operator m;
-
-    // Check if the bound vector is correct
-    ASSERT_NEAR(m.element(bound_vec, 0, 0), 2., isclose);
-    ASSERT_NEAR(m.element(bound_vec, 1, 0), 4., isclose);
-    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
-                isclose);  // atan(2)
-    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
-                isclose);  // atan(sqrt(5)/3)
-    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
-
-    // Check if the same free vector is obtained
-    for (int i = 0; i < 8; i++) {
-        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
-                    isclose);
-    }
-}
-
-// This test cartesian3 coordinate
-TEST(test_host_basics, cartesian3) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const cartesian3<transform3> c3;
-    const point3 global1 = {4., 7., 5.};
-    const vector3 mom = {1., 2., 3.};
-    const vector3 d = vector::normalize(mom);
-    struct dummy_mask {
-    } mask;
-
-    // Global to local transformation
-    const point3 local = c3.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], 2., isclose);
-    ASSERT_NEAR(local[1], 4., isclose);
-    ASSERT_NEAR(local[2], 1., isclose);
-
-    // Local to global transformation
-    const point3 global2 = c3.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-}
-
-// This test polar2 coordinate
-TEST(test_host_basics, polar2) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const polar2<transform3> p2;
-    const point3 global1 = {4., 7., 4.};
-    const vector3 mom = {1., 2., 3.};
-    const vector3 d = vector::normalize(mom);
-    const scalar time = 0.1;
-    const scalar charge = -1.;
-    struct dummy_mask {
-    } mask;
-
-    // Global to local transformation
-    const point2 local = p2.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], std::sqrt(20.), isclose);
-    ASSERT_NEAR(local[1], atan2(4., 2.), isclose);
-
-    // Local to global transformation
-    const point3 global2 = p2.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-
-    // Free track parameter
-    const free_track_parameters<transform3> free_params(global1, time, mom,
-                                                        charge);
-    const auto free_vec1 = free_params.vector();
-
-    const auto bound_vec = p2.free_to_bound_vector(trf, free_vec1);
-    const auto free_vec2 = p2.bound_to_free_vector(trf, mask, bound_vec);
-
-    const matrix_operator m;
-
-    // Check if the bound vector is correct
-    ASSERT_NEAR(m.element(bound_vec, 0, 0), std::sqrt(20.), isclose);
-    ASSERT_NEAR(m.element(bound_vec, 1, 0), atan2(4., 2.), isclose);
-    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
-                isclose);  // atan(2)
-    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
-                isclose);  // atan(sqrt(5)/3)
-    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
-
-    // Check if the same free vector is obtained
-    for (int i = 0; i < 8; i++) {
-        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
-                    isclose);
-    }
-}
-
-// This test cylindrical2 coordinate
-TEST(test_host_basics, cylindrical2) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const cylindrical2<transform3> c2;
-    // Global position on surface
-    const point3 global1 = {scalar{3.4142136}, scalar{4.4142136}, scalar{9.}};
-    const vector3 mom = {1., 2., 3.};
-    const vector3 d = vector::normalize(mom);
-    const scalar time = 0.1;
-    const scalar charge = -1.;
-
-    // Define cylinder mask
-    struct cylinder_mask {
-        scalar r = 0.;
-        scalar operator[](dindex) const { return r; }
-    };
-
-    const scalar r = 2.;
-    const cylinder_mask mask{r};
-
-    // Global to local transformation
-    const point2 local = c2.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], r * scalar{M_PI_4}, isclose);
-    ASSERT_NEAR(local[1], 5., isclose);
-
-    // Local to global transformation
-    const point3 global2 = c2.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-
-    // Free track parameter
-    const free_track_parameters<transform3> free_params(global1, time, mom,
-                                                        charge);
-    const auto free_vec1 = free_params.vector();
-
-    const auto bound_vec = c2.free_to_bound_vector(trf, free_vec1);
-    const auto free_vec2 = c2.bound_to_free_vector(trf, mask, bound_vec);
-
-    const matrix_operator m;
-
-    // Check if the bound vector is correct
-    ASSERT_NEAR(m.element(bound_vec, 0, 0), r * scalar{M_PI_4}, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 1, 0), 5, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 2, 0), 1.1071487,
-                isclose);  // atan(2)
-    ASSERT_NEAR(m.element(bound_vec, 3, 0), 0.64052231,
-                isclose);  // atan(sqrt(5)/3)
-    ASSERT_NEAR(m.element(bound_vec, 4, 0), -1 / 3.7416574, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
-
-    // Check if the same free vector is obtained
-    for (int i = 0; i < 8; i++) {
-        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
-                    isclose);
-    }
-}
-
-// This test cylindrical2 coordinate
-TEST(test_host_basics, cylindrical3) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const cylindrical3<transform3> c3;
-    const point3 global1 = {scalar{3.4142136}, scalar{4.4142136}, scalar{9.}};
-    const vector3 mom = {1., 2., 3.};
-    const vector3 d = vector::normalize(mom);
-    struct dummy_mask {
-    } mask;
-
-    // Global to local transformation
-    const point3 local = c3.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], 2, isclose);
-    ASSERT_NEAR(local[1], M_PI_4, isclose);
-    ASSERT_NEAR(local[2], 5., isclose);
-
-    // Local to global transformation
-    const point3 global2 = c3.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-}
-
-// This test line2 coordinate
-TEST(test_host_basics, line2) {
-
-    // Preparation work
-    const vector3 z = {0., 0., 1.};
-    const vector3 x = {1., 0., 0.};
-    const point3 t = {2., 3., 4.};
-    const transform3 trf(t, z, x);
-    const line2<transform3> l2;
-    const point3 global1 = {3., 3., 9.};
-    const vector3 mom = {0., 2., 0.};
-    const vector3 d = vector::normalize(mom);
-    const scalar time = 0.1;
-    const scalar charge = -1.;
-    struct dummy_mask {
-    } mask;
-
-    // Global to local transformation
-    const point2 local = l2.global_to_local(trf, global1, d);
-
-    // Check if the local position is correct
-    ASSERT_NEAR(local[0], -1., isclose);
-    ASSERT_NEAR(local[1], 5., isclose);
-
-    // Local to global transformation
-    const point3 global2 = l2.local_to_global(trf, mask, local, d);
-
-    // Check if the same global position is obtained
-    ASSERT_NEAR(global1[0], global2[0], isclose);
-    ASSERT_NEAR(global1[1], global2[1], isclose);
-    ASSERT_NEAR(global1[2], global2[2], isclose);
-
-    // Free track parameter
-    const free_track_parameters<transform3> free_params(global1, time, mom,
-                                                        charge);
-    const auto free_vec1 = free_params.vector();
-
-    const auto bound_vec = l2.free_to_bound_vector(trf, free_vec1);
-    const auto free_vec2 = l2.bound_to_free_vector(trf, mask, bound_vec);
-
-    const matrix_operator m;
-
-    // Check if the bound vector is correct
-    ASSERT_NEAR(m.element(bound_vec, 0, 0), -1., isclose);
-    ASSERT_NEAR(m.element(bound_vec, 1, 0), 5, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 2, 0), M_PI_2, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 3, 0), M_PI_2, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 4, 0), -0.5, isclose);
-    ASSERT_NEAR(m.element(bound_vec, 5, 0), 0.1, isclose);
-
-    // Check if the same free vector is obtained
-    for (int i = 0; i < 8; i++) {
-        ASSERT_NEAR(m.element(free_vec1, i, 0), m.element(free_vec2, i, 0),
-                    isclose);
-    }
-}
\ No newline at end of file
diff --git a/tests/common/include/tests/common/covariance_engine.inl b/tests/common/include/tests/common/covariance_engine.inl
deleted file mode 100644
index f970b20d1..000000000
--- a/tests/common/include/tests/common/covariance_engine.inl
+++ /dev/null
@@ -1,16 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-// detray include(s)
-#include "detray/definitions/units.hpp"
-#include "detray/propagator/detail/covariance_engine.hpp"
-#include "detray/tracks/tracks.hpp"
-
-// google-test include(s)
-#include <gtest/gtest.h>
-
-TEST(covariance_engine, jacobian_coordinate) {}
\ No newline at end of file
diff --git a/tests/common/include/tests/common/line_stepper.inl b/tests/common/include/tests/common/line_stepper.inl
new file mode 100644
index 000000000..60fe76f95
--- /dev/null
+++ b/tests/common/include/tests/common/line_stepper.inl
@@ -0,0 +1,116 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/definitions/units.hpp"
+#include "detray/propagator/actor_chain.hpp"
+#include "detray/propagator/actors/parameter_resetter.hpp"
+#include "detray/propagator/actors/parameter_transporter.hpp"
+#include "detray/propagator/line_stepper.hpp"
+#include "detray/propagator/navigator.hpp"
+#include "detray/propagator/propagator.hpp"
+#include "tests/common/tools/create_telescope_detector.hpp"
+
+// Covfie include(s).
+#include <covfie/core/field.hpp>
+#include <covfie/core/field_view.hpp>
+
+// Vecmem include(s).
+#include <vecmem/memory/host_memory_resource.hpp>
+
+// google-test include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+using matrix_operator = standard_matrix_operator<scalar>;
+using transform3 = __plugin::transform3<detray::scalar>;
+
+constexpr scalar epsilon = 1e-6;
+
+TEST(line_stepper, covariance_transport) {
+
+    vecmem::host_memory_resource host_mr;
+
+    // Use rectangular surfaces
+    constexpr bool unbounded = true;
+
+    // Create telescope detector with a single plane
+    detail::ray<transform3> traj{{0, 0, 0}, 0, {1, 0, 0}, -1};
+    std::vector<scalar> positions = {0., 10., 20., 30., 40., 50., 60.};
+
+    const auto det =
+        create_telescope_detector<unbounded>(host_mr, positions, traj);
+
+    using navigator_t = navigator<decltype(det)>;
+    using cline_stepper_t = line_stepper<transform3, constrained_step<>>;
+    using actor_chain_t = actor_chain<dtuple, parameter_transporter<transform3>,
+                                      parameter_resetter<transform3>>;
+    using propagator_t =
+        propagator<cline_stepper_t, navigator_t, actor_chain_t>;
+
+    // Bound vector
+    typename bound_track_parameters<transform3>::vector_type bound_vector;
+    getter::element(bound_vector, e_bound_loc0, 0) = 0.;
+    getter::element(bound_vector, e_bound_loc1, 0) = 0.;
+    getter::element(bound_vector, e_bound_phi, 0) = 0.;
+    getter::element(bound_vector, e_bound_theta, 0) = M_PI / 4.;
+    getter::element(bound_vector, e_bound_qoverp, 0) = -1. / 10.;
+    getter::element(bound_vector, e_bound_time, 0) = 0.;
+
+    // Bound covariance
+    typename bound_track_parameters<transform3>::covariance_type bound_cov =
+        matrix_operator().template identity<e_bound_size, e_bound_size>();
+
+    // Note: Set angle error as ZERO, to constrain the loc0 and loc1 divergence
+    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 0.;
+    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
+
+    // Bound track parameter
+    const bound_track_parameters<transform3> bound_param0(0, bound_vector,
+                                                          bound_cov);
+
+    // Actors
+    parameter_transporter<transform3>::state bound_updater{};
+    parameter_resetter<transform3>::state rst{};
+    actor_chain_t::state actor_states = std::tie(bound_updater, rst);
+
+    propagator_t p({}, {});
+    propagator_t::state propagation(bound_param0, det, actor_states);
+
+    // Run propagator
+    p.propagate(propagation);
+
+    // Bound state after one turn propagation
+    const auto &bound_param1 = propagation._stepping._bound_params;
+
+    // const auto bound_vec0 = bound_param0.vector();
+    // const auto bound_vec1 = bound_param1.vector();
+
+    // Check if the track reaches the final surface
+    EXPECT_EQ(bound_param0.surface_link(), 0);
+    EXPECT_EQ(bound_param1.surface_link(), 5);
+
+    const auto bound_cov0 = bound_param0.covariance();
+    const auto bound_cov1 = bound_param1.covariance();
+
+    // Check covaraince
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
+                        matrix_operator().element(bound_cov1, i, j), epsilon);
+        }
+    }
+
+    // Check covaraince
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+
+            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
+                        matrix_operator().element(bound_cov1, i, j), epsilon);
+        }
+    }
+}
diff --git a/tests/common/include/tests/common/path_correction.inl b/tests/common/include/tests/common/path_correction.inl
new file mode 100644
index 000000000..08ed4e066
--- /dev/null
+++ b/tests/common/include/tests/common/path_correction.inl
@@ -0,0 +1,524 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+// Project include(s).
+#include "detray/definitions/units.hpp"
+#include "detray/materials/predefined_materials.hpp"
+#include "detray/propagator/actor_chain.hpp"
+#include "detray/propagator/actors/parameter_resetter.hpp"
+#include "detray/propagator/actors/parameter_transporter.hpp"
+#include "detray/propagator/line_stepper.hpp"
+#include "detray/propagator/navigator.hpp"
+#include "detray/propagator/propagator.hpp"
+#include "detray/propagator/rk_stepper.hpp"
+#include "detray/tracks/tracks.hpp"
+#include "tests/common/tools/detector_metadata.hpp"
+
+// Covfie include(s).
+#include <covfie/core/field.hpp>
+#include <covfie/core/field_view.hpp>
+
+// Vecmem include(s).
+#include <vecmem/memory/host_memory_resource.hpp>
+
+// google-test include(s).
+#include <gtest/gtest.h>
+
+using namespace detray;
+
+struct surface_targeter : actor {
+
+    struct state {
+        scalar _path;
+        dindex _target_surface_index = dindex_invalid;
+    };
+
+    /// Enforces thepath limit on a stepper state
+    ///
+    /// @param abrt_state contains the path limit
+    /// @param prop_state state of the propagation
+    template <typename propagator_state_t>
+    DETRAY_HOST_DEVICE void operator()(state &actor_state,
+                                       propagator_state_t &prop_state) const {
+
+        prop_state._heartbeat = true;
+        auto &navigation = prop_state._navigation;
+        auto &stepping = prop_state._stepping;
+        navigation.set_full_trust();
+
+        scalar residual = actor_state._path - stepping.path_length();
+        stepping.set_constraint(residual);
+
+        typename propagator_state_t::navigator_state_type::intersection_t is;
+        is.index = actor_state._target_surface_index;
+        is.path = residual;
+        auto &candidates = navigation.candidates();
+        candidates.clear();
+        candidates.push_back(is);
+        navigation.set_next(candidates.begin());
+        navigation.set_unknown();
+
+        if (residual < std::abs(1e-6)) {
+            prop_state._heartbeat = false;
+            candidates.push_back({});
+            navigation.next()++;
+            navigation.set_on_module();
+        }
+    }
+};
+
+// Type Definitions
+using vector2 = __plugin::vector2<scalar>;
+using vector3 = __plugin::vector3<scalar>;
+using matrix_operator = standard_matrix_operator<scalar>;
+using registry_type = detector_registry::default_detector;
+using detector_type =
+    detector<registry_type, covfie::field, std::array, std::tuple,
+             vecmem::vector, vecmem::jagged_vector>;
+using mask_container = typename detector_type::mask_container;
+using material_container = typename detector_type::material_container;
+using transform3 = typename detector_type::transform3;
+using surface_type = typename detector_type::surface_container::value_type;
+using mask_link_type = typename surface_type::mask_link;
+using material_link_type = typename surface_type::material_link;
+using mag_field_t = detector_type::bfield_type;
+using navigator_t = navigator<detector_type>;
+using crk_stepper_t =
+    rk_stepper<mag_field_t::view_t, transform3, constrained_step<>>;
+using actor_chain_t =
+    actor_chain<dtuple, surface_targeter, parameter_transporter<transform3>,
+                parameter_resetter<transform3>>;
+using propagator_t = propagator<crk_stepper_t, navigator_t, actor_chain_t>;
+
+namespace env {
+
+constexpr scalar epsilon = 1e-3;
+constexpr scalar rk_tolerance = 1e-8;
+
+// B field
+const vector3 B{0, 0, 1. * unit_constants::T};
+mag_field_t mag_field(mag_field_t::backend_t::configuration_t{B[0], B[1],
+                                                              B[2]});
+
+// memory resource
+vecmem::host_memory_resource resource;
+
+// Context
+const typename detector_type::context ctx{};
+
+}  // namespace env
+
+// This tests the path correction of cartesian coordinate
+TEST(path_correction, cartesian) {
+
+    // Create a detector
+    detector_type det(env::resource, std::move(env::mag_field));
+
+    // Mask and material ID
+    constexpr registry_type::mask_ids mask_id =
+        registry_type::mask_ids::e_rectangle2;
+    constexpr registry_type::material_ids material_id =
+        registry_type::material_ids::e_slab;
+
+    // Add a volume
+    det.new_volume({0., 0., 0., 0., -M_PI, M_PI});
+
+    typename detector_type::surface_container surfaces(&env::resource);
+    typename detector_type::transform_container transforms(env::resource);
+    typename detector_type::mask_container masks(env::resource);
+    typename detector_type::material_container materials(env::resource);
+
+    // Add a surface
+    const auto trf_index = transforms.size(env::ctx);
+    mask_link_type mask_link{mask_id, masks.template size<mask_id>()};
+    material_link_type material_link{material_id,
+                                     materials.template size<material_id>()};
+    surfaces.emplace_back(trf_index, mask_link, material_link, 0,
+                          dindex_invalid, false);
+
+    // Add a transform
+    const vector3 t{0, 0, 0};
+    const vector3 z{1, 0, 0};
+    const vector3 x{0, 1, 0};
+    transforms.emplace_back(env::ctx, t, z, x);
+
+    // Add a mask
+    const scalar hx = 100. * unit_constants::mm;
+    const scalar hy = 100. * unit_constants::mm;
+    masks.template add_value<mask_id>(0UL, hx, hy);
+
+    // Add a material
+    const material<scalar> mat = silicon<scalar>();
+    const scalar thickness = 2 * unit_constants::mm;
+    materials.template add_value<material_id>(mat, thickness);
+
+    typename detector_type::volume_type &vol = det.volume_by_index(0);
+    det.add_objects_per_volume(env::ctx, vol, surfaces, masks, materials,
+                               transforms);
+
+    // Generate track starting point
+    vector2 local{2, 3};
+    vector3 mom{1 * unit_constants::MeV, 0., 0.};
+
+    scalar time = 0.;
+    scalar q = -1.;
+
+    // bound vector
+    typename bound_track_parameters<transform3>::vector_type bound_vector;
+    getter::element(bound_vector, e_bound_loc0, 0) = local[0];
+    getter::element(bound_vector, e_bound_loc1, 0) = local[1];
+    getter::element(bound_vector, e_bound_phi, 0) = getter::phi(mom);
+    getter::element(bound_vector, e_bound_theta, 0) = getter::theta(mom);
+    getter::element(bound_vector, e_bound_qoverp, 0) = q / getter::norm(mom);
+    getter::element(bound_vector, e_bound_time, 0) = time;
+
+    // bound covariance
+    typename bound_track_parameters<transform3>::covariance_type bound_cov =
+        matrix_operator().template zero<e_bound_size, e_bound_size>();
+    getter::element(bound_cov, e_bound_loc0, e_bound_loc0) = 1.;
+    getter::element(bound_cov, e_bound_loc1, e_bound_loc1) = 1.;
+    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 1.;
+    // Note: Set theta error as ZERO, to constrain the loc1 divergence
+    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
+    getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 1.;
+    getter::element(bound_cov, e_bound_time, e_bound_time) = 1.;
+
+    // bound track parameter
+    const bound_track_parameters<transform3> bound_param0(0, bound_vector,
+                                                          bound_cov);
+
+    // Path length per turn
+    scalar S = 2. * getter::perp(mom) / getter::norm(env::B) * M_PI;
+
+    // Actors
+    surface_targeter::state targeter{S, 0};
+    parameter_transporter<transform3>::state bound_updater{};
+    parameter_resetter<transform3>::state rst{};
+
+    actor_chain_t::state actor_states = std::tie(targeter, bound_updater, rst);
+
+    propagator_t p({}, {});
+    propagator_t::state propagation(bound_param0, env::mag_field, det,
+                                    actor_states);
+
+    crk_stepper_t::state &crk_state = propagation._stepping;
+
+    // Decrease tolerance down to 1e-8
+    crk_state._tolerance = env::rk_tolerance;
+
+    // Run propagator
+    p.propagate(propagation);
+
+    EXPECT_NEAR(crk_state.path_length(), targeter._path, env::epsilon);
+
+    // Bound state after one turn propagation
+    const auto bound_param1 = crk_state._bound_params;
+    const auto bound_vec0 = bound_param0.vector();
+    const auto bound_vec1 = bound_param1.vector();
+    const auto bound_cov0 = bound_param0.covariance();
+    const auto bound_cov1 = bound_param1.covariance();
+
+    // Check if the bound state stays the same after one turn propagation
+
+    // vector
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        EXPECT_NEAR(matrix_operator().element(bound_vec0, i, 0),
+                    matrix_operator().element(bound_vec1, i, 0), env::epsilon);
+    }
+
+    // covaraince
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
+                        matrix_operator().element(bound_cov1, i, j),
+                        env::epsilon);
+        }
+    }
+
+    /*
+    // Print the matrix elements
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            printf("%f ", matrix_operator().element(bound_cov0, i, j));
+        }
+        printf("\n");
+    }
+    printf("\n");
+
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            printf("%f ", matrix_operator().element(bound_cov1, i, j));
+        }
+        printf("\n");
+    }
+    printf("\n");
+    */
+}
+
+TEST(path_correction, polar) {
+
+    // Create a detector
+    detector_type det(env::resource);
+
+    // Mask and material ID
+    constexpr registry_type::mask_ids mask_id =
+        registry_type::mask_ids::e_ring2;
+    constexpr registry_type::material_ids material_id =
+        registry_type::material_ids::e_slab;
+
+    // Add a volume
+    det.new_volume({0., 0., 0., 0., -M_PI, M_PI});
+
+    typename detector_type::surface_container surfaces(&env::resource);
+    typename detector_type::transform_container transforms(env::resource);
+    typename detector_type::mask_container masks(env::resource);
+    typename detector_type::material_container materials(env::resource);
+
+    // Add a surface
+    const auto trf_index = transforms.size(env::ctx);
+    mask_link_type mask_link{mask_id, masks.template size<mask_id>()};
+    material_link_type material_link{material_id,
+                                     materials.template size<material_id>()};
+    surfaces.emplace_back(trf_index, mask_link, material_link, 0,
+                          dindex_invalid, false);
+
+    // Add a transform
+    const vector3 t{0, 0, 0};
+    const vector3 z{1, 0, 0};
+    const vector3 x{0, 1, 0};
+    transforms.emplace_back(env::ctx, t, z, x);
+
+    // Add a mask
+    const scalar r_low = 0. * unit_constants::mm;
+    const scalar r_high = 100. * unit_constants::mm;
+    masks.template add_value<mask_id>(0UL, r_low, r_high);
+
+    // Add a material
+    const material<scalar> mat = silicon<scalar>();
+    const scalar thickness = 2 * unit_constants::mm;
+    materials.template add_value<material_id>(mat, thickness);
+
+    typename detector_type::volume_type &vol = det.volume_by_index(0);
+    det.add_objects_per_volume(env::ctx, vol, surfaces, masks, materials,
+                               transforms);
+
+    // Generate track starting point
+    vector2 local{2, M_PI / 6.};
+    vector3 mom{1 * unit_constants::MeV, 0., 0.};
+
+    scalar time = 0.;
+    scalar q = -1.;
+
+    // bound vector
+    typename bound_track_parameters<transform3>::vector_type bound_vector;
+    getter::element(bound_vector, e_bound_loc0, 0) = local[0];
+    getter::element(bound_vector, e_bound_loc1, 0) = local[1];
+    getter::element(bound_vector, e_bound_phi, 0) = getter::phi(mom);
+    getter::element(bound_vector, e_bound_theta, 0) = getter::theta(mom);
+    getter::element(bound_vector, e_bound_qoverp, 0) = q / getter::norm(mom);
+    getter::element(bound_vector, e_bound_time, 0) = time;
+
+    // bound covariance
+    typename bound_track_parameters<transform3>::covariance_type bound_cov =
+        matrix_operator().template zero<e_bound_size, e_bound_size>();
+    getter::element(bound_cov, e_bound_loc0, e_bound_loc0) = 1.;
+    getter::element(bound_cov, e_bound_loc1, e_bound_loc1) = 1.;
+    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 1.;
+    // Note: Set theta error as ZERO, to constrain the loc1 divergence
+    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
+    getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 1.;
+    getter::element(bound_cov, e_bound_time, e_bound_time) = 1.;
+
+    // bound track parameter
+    const bound_track_parameters<transform3> bound_param0(0, bound_vector,
+                                                          bound_cov);
+
+    // Path length per turn
+    scalar S = 2. * getter::perp(mom) / getter::norm(env::B) * M_PI;
+
+    // Actors
+    surface_targeter::state targeter{S, 0};
+    parameter_transporter<transform3>::state bound_updater{};
+    parameter_resetter<transform3>::state rst{};
+
+    actor_chain_t::state actor_states = std::tie(targeter, bound_updater, rst);
+
+    propagator_t p({}, {});
+    propagator_t::state propagation(bound_param0, env::mag_field, det,
+                                    actor_states);
+
+    crk_stepper_t::state &crk_state = propagation._stepping;
+
+    // RK stepper and its state
+    EXPECT_FLOAT_EQ(crk_state().pos()[0], 0);
+    EXPECT_FLOAT_EQ(crk_state().pos()[1], 2 * std::cos(M_PI / 6.));
+    EXPECT_FLOAT_EQ(crk_state().pos()[2], 2 * std::sin(M_PI / 6.));
+
+    // Decrease tolerance down to 1e-8
+    crk_state._tolerance = env::rk_tolerance;
+
+    // Run propagator
+    p.propagate(propagation);
+
+    EXPECT_NEAR(crk_state.path_length(), targeter._path, env::epsilon);
+
+    // Bound state after one turn propagation
+    const auto bound_param1 = crk_state._bound_params;
+    const auto bound_vec0 = bound_param0.vector();
+    const auto bound_vec1 = bound_param1.vector();
+    const auto bound_cov0 = bound_param0.covariance();
+    const auto bound_cov1 = bound_param1.covariance();
+
+    // Check if the bound state stays the same after one turn propagation
+
+    // vector
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        EXPECT_NEAR(matrix_operator().element(bound_vec0, i, 0),
+                    matrix_operator().element(bound_vec1, i, 0), env::epsilon);
+    }
+
+    // covaraince
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
+                        matrix_operator().element(bound_cov1, i, j),
+                        env::epsilon);
+        }
+    }
+}
+
+TEST(path_correction, cylindrical) {
+
+    // Create a detector
+    detector_type det(env::resource);
+
+    // Mask and material ID
+    constexpr registry_type::mask_ids mask_id =
+        registry_type::mask_ids::e_cylinder2;
+    constexpr registry_type::material_ids material_id =
+        registry_type::material_ids::e_slab;
+
+    // Add a volume
+    det.new_volume({0., 0., 0., 0., -M_PI, M_PI});
+
+    typename detector_type::surface_container surfaces(&env::resource);
+    typename detector_type::transform_container transforms(env::resource);
+    typename detector_type::mask_container masks(env::resource);
+    typename detector_type::material_container materials(env::resource);
+
+    // Add a surface
+    const auto trf_index = transforms.size(env::ctx);
+    mask_link_type mask_link{mask_id, masks.template size<mask_id>()};
+    material_link_type material_link{material_id,
+                                     materials.template size<material_id>()};
+    surfaces.emplace_back(trf_index, mask_link, material_link, 0,
+                          dindex_invalid, false);
+
+    // Add a transform
+    const vector3 t{-50 * unit_constants::mm, 0, 0};
+    const vector3 z{0, 0, 1};
+    const vector3 x{1, 0, 0};
+    transforms.emplace_back(env::ctx, t, z, x);
+
+    // Add a mask
+    const scalar r = 50 * unit_constants::mm;
+    const scalar half_length_1 = 1000. * unit_constants::mm;
+    const scalar half_length_2 = 1000. * unit_constants::mm;
+    masks.template add_value<mask_id>(0UL, r, half_length_1, half_length_2);
+
+    // Add a material
+    const material<scalar> mat = silicon<scalar>();
+    const scalar thickness = 2 * unit_constants::mm;
+    materials.template add_value<material_id>(mat, thickness);
+
+    typename detector_type::volume_type &vol = det.volume_by_index(0);
+    det.add_objects_per_volume(env::ctx, vol, surfaces, masks, materials,
+                               transforms);
+
+    // Generate track starting point
+    vector2 local{0., 0.};
+    vector3 mom{1 * unit_constants::MeV, 0., 0.};
+    scalar time = 0.;
+    scalar q = -1.;
+
+    // bound vector
+    typename bound_track_parameters<transform3>::vector_type bound_vector;
+    getter::element(bound_vector, e_bound_loc0, 0) = local[0];
+    getter::element(bound_vector, e_bound_loc1, 0) = local[1];
+    getter::element(bound_vector, e_bound_phi, 0) = getter::phi(mom);
+    getter::element(bound_vector, e_bound_theta, 0) = getter::theta(mom);
+    getter::element(bound_vector, e_bound_qoverp, 0) = q / getter::norm(mom);
+    getter::element(bound_vector, e_bound_time, 0) = time;
+
+    // bound covariance
+    typename bound_track_parameters<transform3>::covariance_type bound_cov =
+        matrix_operator().template zero<e_bound_size, e_bound_size>();
+    getter::element(bound_cov, e_bound_loc0, e_bound_loc0) = 1.;
+    getter::element(bound_cov, e_bound_loc1, e_bound_loc1) = 1.;
+    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 1.;
+    // Note: Set theta error as ZERO, to constrain the loc1 divergence
+    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
+    getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 1.;
+    getter::element(bound_cov, e_bound_time, e_bound_time) = 1.;
+
+    // bound track parameter
+    const bound_track_parameters<transform3> bound_param0(0, bound_vector,
+                                                          bound_cov);
+
+    // Path length per turn
+    scalar S = 2. * getter::perp(mom) / getter::norm(env::B) * M_PI;
+
+    // Actors
+    surface_targeter::state targeter{S, 0};
+    parameter_transporter<transform3>::state bound_updater{};
+    parameter_resetter<transform3>::state rst{};
+
+    actor_chain_t::state actor_states = std::tie(targeter, bound_updater, rst);
+
+    propagator_t p({}, {});
+    propagator_t::state propagation(bound_param0, env::mag_field, det,
+                                    actor_states);
+
+    crk_stepper_t::state &crk_state = propagation._stepping;
+
+    // RK stepper and its state
+    EXPECT_FLOAT_EQ(crk_state().pos()[0], 0);
+    EXPECT_FLOAT_EQ(crk_state().pos()[1], 0);  // y
+    EXPECT_FLOAT_EQ(crk_state().pos()[2], 0);  // z
+
+    // Decrease tolerance down to 1e-8
+    crk_state._tolerance = env::rk_tolerance;
+
+    // Run propagator
+    p.propagate(propagation);
+
+    EXPECT_NEAR(crk_state.path_length(), targeter._path, env::epsilon);
+
+    // Bound state after one turn propagation
+    const auto bound_param1 = crk_state._bound_params;
+    const auto bound_vec0 = bound_param0.vector();
+    const auto bound_vec1 = bound_param1.vector();
+    const auto bound_cov0 = bound_param0.covariance();
+    const auto bound_cov1 = bound_param1.covariance();
+
+    // Check if the bound state stays the same after one turn propagation
+
+    // vector
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        EXPECT_NEAR(matrix_operator().element(bound_vec0, i, 0),
+                    matrix_operator().element(bound_vec1, i, 0), env::epsilon);
+    }
+
+    // covaraince
+    for (std::size_t i = 0; i < e_bound_size; i++) {
+        for (std::size_t j = 0; j < e_bound_size; j++) {
+            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
+                        matrix_operator().element(bound_cov1, i, j),
+                        env::epsilon);
+        }
+    }
+}
\ No newline at end of file
diff --git a/tests/common/include/tests/common/test_telescope_detector.inl b/tests/common/include/tests/common/test_telescope_detector.inl
index 8a059dc8e..d5c3b8d56 100644
--- a/tests/common/include/tests/common/test_telescope_detector.inl
+++ b/tests/common/include/tests/common/test_telescope_detector.inl
@@ -33,9 +33,8 @@ struct prop_state {
 
     stepping_t _stepping;
     navigation_t _navigation;
-    using field_type = typename stepping_t::field_type;
 
-    template <typename track_t>
+    template <typename track_t, typename field_type>
     prop_state(const track_t &t_in, const field_type &field,
                const typename navigation_t::detector_type &det)
         : _stepping(t_in, field), _navigation(det) {}
diff --git a/tests/common/include/tests/common/tools_guided_navigator.inl b/tests/common/include/tests/common/tools_guided_navigator.inl
index c0505bdea..239a2604b 100644
--- a/tests/common/include/tests/common/tools_guided_navigator.inl
+++ b/tests/common/include/tests/common/tools_guided_navigator.inl
@@ -10,8 +10,8 @@
 #include <vecmem/memory/host_memory_resource.hpp>
 
 #include "detray/definitions/units.hpp"
-#include "detray/propagator/aborters.hpp"
 #include "detray/propagator/actor_chain.hpp"
+#include "detray/propagator/actors/aborters.hpp"
 #include "detray/propagator/navigation_policies.hpp"
 #include "detray/propagator/navigator.hpp"
 #include "detray/propagator/propagator.hpp"
diff --git a/tests/common/include/tests/common/tools_propagator.inl b/tests/common/include/tests/common/tools_propagator.inl
index bd70f749a..facc51060 100644
--- a/tests/common/include/tests/common/tools_propagator.inl
+++ b/tests/common/include/tests/common/tools_propagator.inl
@@ -11,15 +11,14 @@
 
 #include "detray/definitions/units.hpp"
 #include "detray/intersection/detail/trajectories.hpp"
-#include "detray/propagator/aborters.hpp"
 #include "detray/propagator/actor_chain.hpp"
+#include "detray/propagator/actors/aborters.hpp"
 #include "detray/propagator/base_actor.hpp"
 #include "detray/propagator/line_stepper.hpp"
 #include "detray/propagator/navigator.hpp"
 #include "detray/propagator/propagator.hpp"
 #include "detray/propagator/rk_stepper.hpp"
 #include "detray/tracks/tracks.hpp"
-#include "tests/common/tools/create_telescope_detector.hpp"
 #include "tests/common/tools/create_toy_geometry.hpp"
 #include "tests/common/tools/inspectors.hpp"
 #include "tests/common/tools/track_generators.hpp"
diff --git a/tests/common/include/tests/common/tools_stepper.inl b/tests/common/include/tests/common/tools_stepper.inl
index 8b87c4912..abc1c8a4e 100644
--- a/tests/common/include/tests/common/tools_stepper.inl
+++ b/tests/common/include/tests/common/tools_stepper.inl
@@ -5,7 +5,7 @@
  * Mozilla Public License Version 2.0
  */
 
-// detray include(s)
+// Project include(s)
 #include "detray/core/transform_store.hpp"
 #include "detray/core/type_registry.hpp"
 #include "detray/definitions/units.hpp"
@@ -234,134 +234,3 @@ TEST(ALGEBRA_PLUGIN, rk_stepper) {
         EXPECT_NEAR(getter::norm(backward_relative_error), 0, epsilon);
     }
 }
-
-// This tests the covariance transport in rk stepper
-TEST(ALGEBRA_PLUGIN, covariance_transport) {
-
-    /*
-
-    // test surface
-    const vector3 u{0, 1, 0};
-    const vector3 w{1, 0, 0};
-    const vector3 t{0, 0, 0};
-    const transform3 trf(t, w, u);
-
-    // Generate track starting point
-    vector3 local{2, 3, 0};
-    vector3 mom{0.02, 0., 0.};
-    scalar time = 0.;
-    scalar q = -1.;
-
-    // bound vector
-    typename bound_track_parameters<transform3>::vector_type bound_vector;
-    getter::element(bound_vector, e_bound_loc0, 0) = local[0];
-    getter::element(bound_vector, e_bound_loc1, 0) = local[1];
-    getter::element(bound_vector, e_bound_phi, 0) = getter::phi(mom);
-    getter::element(bound_vector, e_bound_theta, 0) = getter::theta(mom);
-    getter::element(bound_vector, e_bound_qoverp, 0) = q / getter::norm(mom);
-    getter::element(bound_vector, e_bound_time, 0) = time;
-
-    // bound covariance
-    typename bound_track_parameters<transform3>::covariance_type bound_cov =
-        matrix_operator().template zero<e_bound_size, e_bound_size>();
-    getter::element(bound_cov, e_bound_loc0, e_bound_loc0) = 1.;
-    getter::element(bound_cov, e_bound_loc1, e_bound_loc1) = 1.;
-    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 1.;
-
-    // Note: Set theta error as ZERO, to constrain the loc1 divergence
-    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
-    getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 1.;
-    getter::element(bound_cov, e_bound_time, e_bound_time) = 1.;
-
-    // B field
-    vector3 B{0, 0, 1. * unit_constants::T};
-    mag_field_t mag_field(B);
-
-    // bound track parameter
-    const bound_track_parameters<transform3> bound_param0(0, bound_vector,
-                                                          bound_cov);
-
-    prop_state<crk_stepper_t::state, nav_state> propagation{
-        crk_stepper_t::state(bound_param0, trf, mag_field), nav_state{}};
-    crk_stepper_t::state &crk_state = propagation._stepping;
-    nav_state &n_state = propagation._navigation;
-
-    // Decrease tolerance down to 1e-8
-    crk_state._tolerance = 1e-8;
-
-    // RK stepper and its state
-    crk_stepper_t crk_stepper;
-
-    ASSERT_FLOAT_EQ(crk_state().pos()[0], 0);
-    ASSERT_FLOAT_EQ(crk_state().pos()[1], local[0]);
-    ASSERT_FLOAT_EQ(crk_state().pos()[2], local[1]);
-    ASSERT_NEAR(crk_state().dir()[0], 1, epsilon);
-    ASSERT_NEAR(crk_state().dir()[1], 0, epsilon);
-    ASSERT_NEAR(crk_state().dir()[2], 0, epsilon);
-
-    // helix trajectory
-    detail::helix helix(crk_state(), &B);
-
-    // Path length per turn
-    scalar S = 2. * getter::norm(mom) / getter::norm(B) * M_PI;
-
-    // Run stepper for one turn
-    unsigned int max_steps = 1e4;
-    for (unsigned int i = 0; i < max_steps; i++) {
-
-        crk_state.set_constraint(S - crk_state.path_length());
-
-        n_state._step_size = S;
-
-        crk_stepper.step(propagation);
-
-        if (std::abs(S - crk_state.path_length()) < 1e-6) {
-            break;
-        }
-
-        // Make sure that we didn't reach the end of for loop
-        ASSERT_TRUE(i < max_steps - 1);
-    }
-
-    // Transport jacobian check
-
-    auto jac_transport = crk_state._jac_transport;
-    auto true_J = helix.jacobian(crk_state.path_length());
-
-    for (std::size_t i = 0; i < e_free_size; i++) {
-        for (std::size_t j = 0; j < e_free_size; j++) {
-            EXPECT_NEAR(matrix_operator().element(jac_transport, i, j),
-                        matrix_operator().element(true_J, i, j),
-                        crk_state.path_length() * epsilon);
-        }
-    }
-
-    // Bound parameters check
-
-    // Bound state after one turn propagation
-    const auto bound_param1 = crk_stepper.bound_state(propagation, trf);
-
-    const auto bound_vec0 = bound_param0.vector();
-    const auto bound_vec1 = bound_param1.vector();
-
-    const auto bound_cov0 = bound_param0.covariance();
-    const auto bound_cov1 = bound_param1.covariance();
-
-    // Check if the bound state stays the same after one turn propagation
-
-    // vector
-    for (size_type i = 0; i < e_bound_size; i++) {
-        EXPECT_NEAR(matrix_operator().element(bound_vec0, i, 0),
-                    matrix_operator().element(bound_vec1, i, 0), epsilon);
-    }
-
-    // covaraince
-    for (size_type i = 0; i < e_bound_size; i++) {
-        for (size_type j = 0; j < e_bound_size; j++) {
-            EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j),
-                        matrix_operator().element(bound_cov1, i, j),
-                        crk_state.path_length() * epsilon);
-        }
-    }
-    */
-}
diff --git a/tests/unit_tests/array/CMakeLists.txt b/tests/unit_tests/array/CMakeLists.txt
index 840bbf20a..99ca64bfd 100644
--- a/tests/unit_tests/array/CMakeLists.txt
+++ b/tests/unit_tests/array/CMakeLists.txt
@@ -9,9 +9,13 @@ detray_add_test( array
    "array_actor_chain.cpp"
    "array_annulus2D.cpp"
    "array_container.cpp"
-   "array_coordinates.cpp"
+   "array_coordinate_cartesian2.cpp"
+   "array_coordinate_cartesian3.cpp"
+   "array_coordinate_cylindrical2.cpp"
+   "array_coordinate_cylindrical3.cpp"
+   "array_coordinate_line2.cpp"
+   "array_coordinate_polar2.cpp"
    "array_core.cpp"
-   "array_covariance_engine.cpp"
    "array_cylinder_intersection.cpp"
    "array_cylinder.cpp"
    "array_detector.cpp"
@@ -26,13 +30,15 @@ detray_add_test( array
    "array_guided_navigator.cpp"
    "array_helix_trajectory.cpp"
    "array_intersection_kernel.cpp"
-   "array_particle_gun.cpp"
+   "array_line_intersection.cpp"
+   "array_line_stepper.cpp"
+   "array_line.cpp"
    "array_mask_store.cpp"
    "array_materials.cpp"
    "array_navigator.cpp"
+   "array_particle_gun.cpp"
+   "array_path_correction.cpp"
    "array_planar_intersection.cpp" 
-   "array_line_intersection.cpp"
-   "array_line.cpp"
    "array_propagator.cpp"
    "array_rectangle2D.cpp"
    "array_ring2D.cpp"
diff --git a/tests/unit_tests/array/array_coordinate_cartesian2.cpp b/tests/unit_tests/array/array_coordinate_cartesian2.cpp
new file mode 100644
index 000000000..f6a1ca99c
--- /dev/null
+++ b/tests/unit_tests/array/array_coordinate_cartesian2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/coordinate_cartesian2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_coordinate_cartesian3.cpp b/tests/unit_tests/array/array_coordinate_cartesian3.cpp
new file mode 100644
index 000000000..07da64154
--- /dev/null
+++ b/tests/unit_tests/array/array_coordinate_cartesian3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/coordinate_cartesian3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_coordinate_cylindrical2.cpp b/tests/unit_tests/array/array_coordinate_cylindrical2.cpp
new file mode 100644
index 000000000..0b981467b
--- /dev/null
+++ b/tests/unit_tests/array/array_coordinate_cylindrical2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/coordinate_cylindrical2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_coordinate_cylindrical3.cpp b/tests/unit_tests/array/array_coordinate_cylindrical3.cpp
new file mode 100644
index 000000000..31195c45a
--- /dev/null
+++ b/tests/unit_tests/array/array_coordinate_cylindrical3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/coordinate_cylindrical3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_covariance_engine.cpp b/tests/unit_tests/array/array_coordinate_line2.cpp
similarity index 82%
rename from tests/unit_tests/array/array_covariance_engine.cpp
rename to tests/unit_tests/array/array_coordinate_line2.cpp
index 1bf23ee84..b87112ab6 100644
--- a/tests/unit_tests/array/array_covariance_engine.cpp
+++ b/tests/unit_tests/array/array_coordinate_line2.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/array_definitions.hpp"
-#include "tests/common/covariance_engine.inl"
+#include "tests/common/coordinate_line2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_coordinate_polar2.cpp b/tests/unit_tests/array/array_coordinate_polar2.cpp
new file mode 100644
index 000000000..2d06429e5
--- /dev/null
+++ b/tests/unit_tests/array/array_coordinate_polar2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/coordinate_polar2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/array/array_coordinates.cpp b/tests/unit_tests/array/array_line_stepper.cpp
similarity index 83%
rename from tests/unit_tests/array/array_coordinates.cpp
rename to tests/unit_tests/array/array_line_stepper.cpp
index 4d3e112ed..d55d144fa 100644
--- a/tests/unit_tests/array/array_coordinates.cpp
+++ b/tests/unit_tests/array/array_line_stepper.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/array_definitions.hpp"
-#include "tests/common/coordinates.inl"
+#include "tests/common/line_stepper.inl"
diff --git a/tests/unit_tests/array/array_path_correction.cpp b/tests/unit_tests/array/array_path_correction.cpp
new file mode 100644
index 000000000..9344e0349
--- /dev/null
+++ b/tests/unit_tests/array/array_path_correction.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/array_definitions.hpp"
+#include "tests/common/path_correction.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/cuda/CMakeLists.txt b/tests/unit_tests/cuda/CMakeLists.txt
index d1007e6e1..1e7b61431 100644
--- a/tests/unit_tests/cuda/CMakeLists.txt
+++ b/tests/unit_tests/cuda/CMakeLists.txt
@@ -42,8 +42,6 @@ foreach(algebra ${algebras})
       "detector_cuda_kernel.cu"
       "navigator_cuda.cpp" "navigator_cuda_kernel.hpp"
       "navigator_cuda_kernel.cu"
-      "rk_stepper_cuda.cpp" "rk_stepper_cuda_kernel.hpp"
-      "rk_stepper_cuda_kernel.cu"
       "sf_finders_grid_cuda.cpp" "sf_finders_grid_cuda_kernel.hpp"
       "sf_finders_grid_cuda_kernel.cu"
       "propagator_cuda.cpp" "propagator_cuda_kernel.hpp"
diff --git a/tests/unit_tests/cuda/propagator_cuda.cpp b/tests/unit_tests/cuda/propagator_cuda.cpp
index d39a1aaa1..46cadbbbe 100644
--- a/tests/unit_tests/cuda/propagator_cuda.cpp
+++ b/tests/unit_tests/cuda/propagator_cuda.cpp
@@ -151,7 +151,6 @@ TEST_P(CudaPropagatorWithRkStepper, propagator) {
             // ASSERT_NEAR((host_pl - device_pl) / host_pl, 0, is_close);
 
             ASSERT_EQ(host_positions[i].size(), device_positions[i].size());
-
             ASSERT_NEAR(host_pl, device_pl, host_pl * is_close);
 
             auto& host_pos = host_positions[i][j];
diff --git a/tests/unit_tests/cuda/propagator_cuda_kernel.hpp b/tests/unit_tests/cuda/propagator_cuda_kernel.hpp
index efe73b075..00b3c52cd 100644
--- a/tests/unit_tests/cuda/propagator_cuda_kernel.hpp
+++ b/tests/unit_tests/cuda/propagator_cuda_kernel.hpp
@@ -21,8 +21,8 @@
 #include <covfie/core/field_view.hpp>
 
 #include "detray/definitions/units.hpp"
-#include "detray/propagator/aborters.hpp"
 #include "detray/propagator/actor_chain.hpp"
+#include "detray/propagator/actors/aborters.hpp"
 #include "detray/propagator/base_actor.hpp"
 #include "detray/propagator/navigator.hpp"
 #include "detray/propagator/propagator.hpp"
@@ -63,7 +63,7 @@ constexpr unsigned int theta_steps{10};
 constexpr unsigned int phi_steps{10};
 
 constexpr scalar rk_tolerance{1e-4};
-constexpr scalar overstep_tolerance{-7 * unit_constants::um};
+constexpr scalar overstep_tolerance{-3 * unit_constants::um};
 constexpr scalar constrainted_step_size{2. * unit_constants::mm};
 constexpr scalar is_close{1e-4};
 constexpr scalar path_limit{2 * unit_constants::m};
diff --git a/tests/unit_tests/cuda/rk_stepper_cuda.cpp b/tests/unit_tests/cuda/rk_stepper_cuda.cpp
deleted file mode 100644
index e047f64c6..000000000
--- a/tests/unit_tests/cuda/rk_stepper_cuda.cpp
+++ /dev/null
@@ -1,127 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#include <gtest/gtest.h>
-
-#include <vecmem/memory/cuda/managed_memory_resource.hpp>
-#include <vecmem/memory/host_memory_resource.hpp>
-
-#include "rk_stepper_cuda_kernel.hpp"
-
-TEST(rk_stepper_cuda, bound_state) {
-
-    // VecMem memory resource(s)
-    vecmem::cuda::managed_memory_resource mng_mr;
-
-    // test surface
-    const vector3 u{0, 1, 0};
-    const vector3 w{1, 0, 0};
-    const vector3 t{0, 0, 0};
-    const transform3 trf(t, w, u);
-
-    // Generate track starting point
-    vector3 local{2, 3, 0};
-    vector3 mom{0.02, 0., 0.};
-    scalar time = 0.;
-    scalar q = -1.;
-
-    // bound vector
-    typename bound_track_parameters<transform3>::vector_type bound_vector;
-    getter::element(bound_vector, e_bound_loc0, 0) = local[0];
-    getter::element(bound_vector, e_bound_loc1, 0) = local[1];
-    getter::element(bound_vector, e_bound_phi, 0) = getter::phi(mom);
-    getter::element(bound_vector, e_bound_theta, 0) = getter::theta(mom);
-    getter::element(bound_vector, e_bound_qoverp, 0) = q / getter::norm(mom);
-    getter::element(bound_vector, e_bound_time, 0) = time;
-
-    // bound covariance
-    typename bound_track_parameters<transform3>::covariance_type bound_cov =
-        matrix_operator().template zero<e_bound_size, e_bound_size>();
-    getter::element(bound_cov, e_bound_loc0, e_bound_loc0) = 1.;
-    getter::element(bound_cov, e_bound_loc1, e_bound_loc1) = 1.;
-    getter::element(bound_cov, e_bound_phi, e_bound_phi) = 1.;
-
-    // Note: Set theta error as ZERO, to constrain the loc1 divergence
-    getter::element(bound_cov, e_bound_theta, e_bound_theta) = 0.;
-    getter::element(bound_cov, e_bound_qoverp, e_bound_qoverp) = 1.;
-
-    // bound track parameter
-    const bound_track_parameters<transform3> in_param(0, bound_vector,
-                                                      bound_cov);
-    const vector3 B{0, 0, 1. * unit_constants::T};
-
-    /**
-     * Get CPU bound parameter after one turn
-     */
-
-    mag_field_t mag_field(
-        mag_field_t::backend_t::configuration_t{B[0], B[1], B[2]});
-    prop_state<crk_stepper_t::state, nav_state> propagation{
-        crk_stepper_t::state(in_param, trf, mag_field), nav_state{}};
-    crk_stepper_t::state &crk_state = propagation._stepping;
-    nav_state &n_state = propagation._navigation;
-
-    // Decrease tolerance down to 1e-8
-    crk_state.set_tolerance(rk_tolerance);
-
-    // RK stepper and its state
-    crk_stepper_t crk_stepper;
-
-    // Path length per turn
-    scalar S = 2. * std::fabs(1. / in_param.qop()) / getter::norm(B) * M_PI;
-
-    // Run stepper for half turn
-    unsigned int max_steps = 1e4;
-
-    for (unsigned int i = 0; i < max_steps; i++) {
-
-        crk_state.set_constraint(S - crk_state.path_length());
-
-        n_state._step_size = S;
-
-        crk_stepper.step(propagation);
-
-        if (std::abs(S - crk_state.path_length()) < 1e-6) {
-            break;
-        }
-
-        // Make sure that we didn't reach the end of for loop
-        ASSERT_TRUE(i < max_steps - 1);
-    }
-
-    // Bound state after one turn propagation
-    const auto out_param_cpu = crk_stepper.bound_state(propagation, trf);
-
-    /**
-     * Get CUDA bound parameter after one turn
-     */
-    vecmem::vector<bound_track_parameters<transform3>> out_param_cuda(1,
-                                                                      &mng_mr);
-
-    bound_state_test(vecmem::get_data(out_param_cuda), in_param, B, trf);
-
-    /**
-     * Compare CPU and CUDA
-     */
-    const auto bvec_cpu = out_param_cpu.vector();
-    const auto bcov_cpu = out_param_cpu.covariance();
-
-    const auto bvec_cuda = out_param_cuda[0].vector();
-    const auto bcov_cuda = out_param_cuda[0].covariance();
-
-    for (std::size_t i = 0; i < e_bound_size; i++) {
-        EXPECT_NEAR(matrix_operator().element(bvec_cpu, i, 0),
-                    matrix_operator().element(bvec_cuda, i, 0), epsilon);
-    }
-
-    for (std::size_t i = 0; i < e_bound_size; i++) {
-        for (std::size_t j = 0; j < e_bound_size; j++) {
-            EXPECT_NEAR(matrix_operator().element(bcov_cpu, i, j),
-                        matrix_operator().element(bcov_cuda, i, j), epsilon);
-        }
-    }
-}
diff --git a/tests/unit_tests/cuda/rk_stepper_cuda_kernel.cu b/tests/unit_tests/cuda/rk_stepper_cuda_kernel.cu
deleted file mode 100644
index 7d8c91e60..000000000
--- a/tests/unit_tests/cuda/rk_stepper_cuda_kernel.cu
+++ /dev/null
@@ -1,88 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#include <covfie/core/backend/primitive/constant.hpp>
-#include <covfie/core/field.hpp>
-#include <covfie/core/field_view.hpp>
-#include <covfie/core/vector.hpp>
-#include <vecmem/containers/device_vector.hpp>
-
-#include "detray/definitions/cuda_definitions.hpp"
-#include "rk_stepper_cuda_kernel.hpp"
-
-namespace {
-using field_t = covfie::field<
-    covfie::backend::constant<covfie::vector::vector_d<detray::scalar, 3>,
-                              covfie::vector::vector_d<detray::scalar, 3>>>;
-}
-
-namespace detray {
-
-__global__ void bound_state_test_kernel(
-    vecmem::data::vector_view<bound_track_parameters<transform3>> out_param,
-    const bound_track_parameters<transform3> in_param, const field_t::view_t B,
-    const transform3 trf) {
-
-    vecmem::device_vector<bound_track_parameters<transform3>> out_param_cuda(
-        out_param);
-
-    prop_state<crk_stepper_t::state, nav_state> propagation{
-        crk_stepper_t::state(in_param, trf, B), nav_state{}};
-    crk_stepper_t::state &crk_state = propagation._stepping;
-    nav_state &n_state = propagation._navigation;
-
-    // Decrease tolerance down to 1e-8
-    crk_state.set_tolerance(rk_tolerance);
-
-    // RK stepper and its state
-    crk_stepper_t crk_stepper;
-
-    vector3 Bv{B.at(0.f, 0.f, 0.f)[0], B.at(0.f, 0.f, 0.f)[1],
-               B.at(0.f, 0.f, 0.f)[2]};
-
-    // Path length per turn
-    scalar S = 2. * std::fabs(1. / in_param.qop()) / getter::norm(Bv) * M_PI;
-
-    // Run stepper for one turn
-    unsigned int max_steps = 1e4;
-    for (unsigned int i = 0; i < max_steps; i++) {
-
-        crk_state.set_constraint(S - crk_state.path_length());
-
-        n_state._step_size = S;
-
-        crk_stepper.step(propagation);
-
-        if (std::abs(S - crk_state.path_length()) < 1e-6) {
-            break;
-        }
-    }
-
-    // Bound state after one turn propagation
-    out_param_cuda[0] = crk_stepper.bound_state(propagation, trf);
-}
-
-void bound_state_test(
-    vecmem::data::vector_view<bound_track_parameters<transform3>> out_param,
-    const bound_track_parameters<transform3> in_param, const vector3 B,
-    const transform3 trf) {
-
-    constexpr int thread_dim = 1;
-    constexpr int block_dim = 1;
-
-    field_t f(field_t::backend_t::configuration_t{B[0], B[1], B[2]});
-
-    // run the test kernel
-    bound_state_test_kernel<<<block_dim, thread_dim>>>(out_param, in_param, f,
-                                                       trf);
-
-    // cuda error check
-    DETRAY_CUDA_ERROR_CHECK(cudaGetLastError());
-    DETRAY_CUDA_ERROR_CHECK(cudaDeviceSynchronize());
-}
-
-}  // namespace detray
diff --git a/tests/unit_tests/cuda/rk_stepper_cuda_kernel.hpp b/tests/unit_tests/cuda/rk_stepper_cuda_kernel.hpp
deleted file mode 100644
index 2cca70132..000000000
--- a/tests/unit_tests/cuda/rk_stepper_cuda_kernel.hpp
+++ /dev/null
@@ -1,97 +0,0 @@
-/** Detray library, part of the ACTS project (R&D line)
- *
- * (c) 2022 CERN for the benefit of the ACTS project
- *
- * Mozilla Public License Version 2.0
- */
-
-#pragma once
-
-#if defined(array)
-#include "detray/plugins/algebra/array_definitions.hpp"
-#elif defined(eigen)
-#include "detray/plugins/algebra/eigen_definitions.hpp"
-#elif defined(smatrix)
-#include "detray/plugins/algebra/smatrix_definitions.hpp"
-#elif defined(vc_array)
-#include "detray/plugins/algebra/vc_array_definitions.hpp"
-#endif
-
-#include <covfie/core/backend/primitive/constant.hpp>
-#include <covfie/core/field.hpp>
-#include <covfie/core/vector.hpp>
-
-#include "detray/definitions/qualifiers.hpp"
-#include "detray/definitions/units.hpp"
-#include "detray/propagator/rk_stepper.hpp"
-#include "detray/tracks/tracks.hpp"
-
-using namespace detray;
-
-namespace {
-
-// type definitions
-using vector3 = __plugin::vector3<scalar>;
-using point3 = __plugin::point3<scalar>;
-using transform3 = __plugin::transform3<scalar>;
-using mag_field_t = covfie::field<covfie::backend::constant<
-    covfie::vector::vector_d<scalar, 3>, covfie::vector::vector_d<scalar, 3>>>;
-using rk_stepper_t = rk_stepper<mag_field_t::view_t, transform3>;
-using crk_stepper_t =
-    rk_stepper<mag_field_t::view_t, transform3, constrained_step<>>;
-using matrix_operator = standard_matrix_operator<scalar>;
-
-// geomery navigation configurations
-constexpr unsigned int theta_steps = 100;
-constexpr unsigned int phi_steps = 100;
-constexpr unsigned int rk_steps = 100;
-
-constexpr scalar epsilon = 1e-3;
-constexpr scalar rk_tolerance = 1e-4;
-
-// dummy navigation struct
-struct nav_state {
-
-    DETRAY_HOST_DEVICE
-    scalar operator()() const { return _step_size; }
-
-    DETRAY_HOST_DEVICE
-    inline auto current_object() const -> dindex { return dindex_invalid; }
-
-    DETRAY_HOST_DEVICE
-    inline void set_full_trust() {}
-
-    DETRAY_HOST_DEVICE
-    inline void set_high_trust() {}
-
-    DETRAY_HOST_DEVICE
-    inline void set_fair_trust() {}
-
-    DETRAY_HOST_DEVICE
-    inline void set_no_trust() {}
-
-    DETRAY_HOST_DEVICE
-    inline bool abort() { return false; }
-
-    scalar _step_size = 1. * unit_constants::mm;
-};
-
-// dummy propagator state
-template <typename stepping_t, typename navigation_t>
-struct prop_state {
-    stepping_t _stepping;
-    navigation_t _navigation;
-};
-
-}  // anonymous namespace
-
-namespace detray {
-
-// Test function for Runge-Kutta stepper bound state
-// This test investigates only one track
-void bound_state_test(
-    vecmem::data::vector_view<bound_track_parameters<transform3>> out_param,
-    const bound_track_parameters<transform3> in_param, const vector3 B,
-    const transform3 trf);
-
-}  // namespace detray
diff --git a/tests/unit_tests/eigen/CMakeLists.txt b/tests/unit_tests/eigen/CMakeLists.txt
index 9de27f73e..eeefeab17 100644
--- a/tests/unit_tests/eigen/CMakeLists.txt
+++ b/tests/unit_tests/eigen/CMakeLists.txt
@@ -9,9 +9,13 @@ detray_add_test( eigen
    "eigen_actor_chain.cpp"
    "eigen_annulus2D.cpp"
    "eigen_container.cpp"
-   "eigen_coordinates.cpp"
+   "eigen_coordinate_cartesian2.cpp"
+   "eigen_coordinate_cartesian3.cpp"
+   "eigen_coordinate_cylindrical2.cpp"
+   "eigen_coordinate_cylindrical3.cpp"
+   "eigen_coordinate_line2.cpp"
+   "eigen_coordinate_polar2.cpp"
    "eigen_core.cpp"
-   "eigen_covariance_engine.cpp"
    "eigen_cylinder_intersection.cpp" 
    "eigen_cylinder.cpp"
    "eigen_detector.cpp"
@@ -26,13 +30,15 @@ detray_add_test( eigen
    "eigen_guided_navigator.cpp"
    "eigen_helix_trajectory.cpp"
    "eigen_intersection_kernel.cpp"
+   "eigen_line_intersection.cpp"
+   "eigen_line_stepper.cpp"
+   "eigen_line.cpp"
    "eigen_mask_store.cpp"
    "eigen_materials.cpp"
    "eigen_navigator.cpp"
    "eigen_particle_gun.cpp"
+   "eigen_path_correction.cpp"
    "eigen_planar_intersection.cpp"
-   "eigen_line_intersection.cpp"
-   "eigen_line.cpp"
    "eigen_propagator.cpp"
    "eigen_rectangle2D.cpp"
    "eigen_ring2D.cpp"
diff --git a/tests/unit_tests/eigen/eigen_coordinate_cartesian2.cpp b/tests/unit_tests/eigen/eigen_coordinate_cartesian2.cpp
new file mode 100644
index 000000000..8ff588882
--- /dev/null
+++ b/tests/unit_tests/eigen/eigen_coordinate_cartesian2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/eigen_definitions.hpp"
+#include "tests/common/coordinate_cartesian2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_coordinate_cartesian3.cpp b/tests/unit_tests/eigen/eigen_coordinate_cartesian3.cpp
new file mode 100644
index 000000000..ff767ecfe
--- /dev/null
+++ b/tests/unit_tests/eigen/eigen_coordinate_cartesian3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/eigen_definitions.hpp"
+#include "tests/common/coordinate_cartesian3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_coordinate_cylindrical2.cpp b/tests/unit_tests/eigen/eigen_coordinate_cylindrical2.cpp
new file mode 100644
index 000000000..aa4f0ae3c
--- /dev/null
+++ b/tests/unit_tests/eigen/eigen_coordinate_cylindrical2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/eigen_definitions.hpp"
+#include "tests/common/coordinate_cylindrical2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_coordinate_cylindrical3.cpp b/tests/unit_tests/eigen/eigen_coordinate_cylindrical3.cpp
new file mode 100644
index 000000000..7a67e3857
--- /dev/null
+++ b/tests/unit_tests/eigen/eigen_coordinate_cylindrical3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/eigen_definitions.hpp"
+#include "tests/common/coordinate_cylindrical3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_covariance_engine.cpp b/tests/unit_tests/eigen/eigen_coordinate_line2.cpp
similarity index 82%
rename from tests/unit_tests/eigen/eigen_covariance_engine.cpp
rename to tests/unit_tests/eigen/eigen_coordinate_line2.cpp
index 724ae9ca6..017406573 100644
--- a/tests/unit_tests/eigen/eigen_covariance_engine.cpp
+++ b/tests/unit_tests/eigen/eigen_coordinate_line2.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/eigen_definitions.hpp"
-#include "tests/common/covariance_engine.inl"
+#include "tests/common/coordinate_line2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_coordinate_polar2.cpp b/tests/unit_tests/eigen/eigen_coordinate_polar2.cpp
new file mode 100644
index 000000000..f148e848f
--- /dev/null
+++ b/tests/unit_tests/eigen/eigen_coordinate_polar2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/eigen_definitions.hpp"
+#include "tests/common/coordinate_polar2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/eigen/eigen_line_stepper.cpp b/tests/unit_tests/eigen/eigen_line_stepper.cpp
new file mode 100644
index 000000000..e69de29bb
diff --git a/tests/unit_tests/eigen/eigen_coordinates.cpp b/tests/unit_tests/eigen/eigen_path_correction.cpp
similarity index 84%
rename from tests/unit_tests/eigen/eigen_coordinates.cpp
rename to tests/unit_tests/eigen/eigen_path_correction.cpp
index a4a961581..5592a0628 100644
--- a/tests/unit_tests/eigen/eigen_coordinates.cpp
+++ b/tests/unit_tests/eigen/eigen_path_correction.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/eigen_definitions.hpp"
-#include "tests/common/coordinates.inl"
\ No newline at end of file
+#include "tests/common/line_stepper.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/CMakeLists.txt b/tests/unit_tests/smatrix/CMakeLists.txt
index 36eb63c7e..39ea56fb9 100644
--- a/tests/unit_tests/smatrix/CMakeLists.txt
+++ b/tests/unit_tests/smatrix/CMakeLists.txt
@@ -9,9 +9,13 @@ detray_add_test( smatrix
    "smatrix_actor_chain.cpp"
    "smatrix_annulus2D.cpp"
    "smatrix_container.cpp"
-   "smatrix_coordinates.cpp"
+   "smatrix_coordinate_cartesian2.cpp"
+   "smatrix_coordinate_cartesian3.cpp"
+   "smatrix_coordinate_cylindrical2.cpp"
+   "smatrix_coordinate_cylindrical3.cpp"
+   "smatrix_coordinate_line2.cpp"
+   "smatrix_coordinate_polar2.cpp"
    "smatrix_core.cpp" 
-   "smatrix_covariance_engine.cpp"
    "smatrix_cylinder_intersection.cpp"
    "smatrix_cylinder.cpp"
    "smatrix_detector.cpp"
@@ -26,12 +30,14 @@ detray_add_test( smatrix
    "smatrix_guided_navigator.cpp"
    "smatrix_helix_trajectory.cpp"
    "smatrix_intersection_kernel.cpp"
+   "smatrix_line_intersection.cpp"
+   "smatrix_line_stepper.cpp"
+   "smatrix_line.cpp"
    "smatrix_mask_store.cpp"
    "smatrix_materials.cpp"
    "smatrix_navigator.cpp"
    "smatrix_particle_gun.cpp"
-   "smatrix_line_intersection.cpp"
-   "smatrix_line.cpp"
+   "smatrix_path_correction.cpp"
    "smatrix_planar_intersection.cpp"
    "smatrix_propagator.cpp"
    "smatrix_rectangle2D.cpp"
diff --git a/tests/unit_tests/smatrix/smatrix_coordinate_cartesian2.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_cartesian2.cpp
new file mode 100644
index 000000000..671a55218
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_cartesian2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/coordinate_cartesian2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_coordinate_cartesian3.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_cartesian3.cpp
new file mode 100644
index 000000000..d3aba71ae
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_cartesian3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/coordinate_cartesian3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical2.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical2.cpp
new file mode 100644
index 000000000..769044bd9
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/coordinate_cylindrical2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical3.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical3.cpp
new file mode 100644
index 000000000..4385c315b
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_cylindrical3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/coordinate_cylindrical3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_covariance_engine.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_line2.cpp
similarity index 82%
rename from tests/unit_tests/smatrix/smatrix_covariance_engine.cpp
rename to tests/unit_tests/smatrix/smatrix_coordinate_line2.cpp
index 86658c813..26d8a9170 100644
--- a/tests/unit_tests/smatrix/smatrix_covariance_engine.cpp
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_line2.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/smatrix_definitions.hpp"
-#include "tests/common/covariance_engine.inl"
+#include "tests/common/coordinate_line2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_coordinate_polar2.cpp b/tests/unit_tests/smatrix/smatrix_coordinate_polar2.cpp
new file mode 100644
index 000000000..cd1563a69
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_coordinate_polar2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/coordinate_polar2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/smatrix/smatrix_coordinates.cpp b/tests/unit_tests/smatrix/smatrix_line_stepper.cpp
similarity index 84%
rename from tests/unit_tests/smatrix/smatrix_coordinates.cpp
rename to tests/unit_tests/smatrix/smatrix_line_stepper.cpp
index b35385cc4..ca05895f4 100644
--- a/tests/unit_tests/smatrix/smatrix_coordinates.cpp
+++ b/tests/unit_tests/smatrix/smatrix_line_stepper.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/smatrix_definitions.hpp"
-#include "tests/common/coordinates.inl"
\ No newline at end of file
+#include "tests/common/line_stepper.inl"
diff --git a/tests/unit_tests/smatrix/smatrix_path_correction.cpp b/tests/unit_tests/smatrix/smatrix_path_correction.cpp
new file mode 100644
index 000000000..fe2905e40
--- /dev/null
+++ b/tests/unit_tests/smatrix/smatrix_path_correction.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/smatrix_definitions.hpp"
+#include "tests/common/path_correction.inl"
diff --git a/tests/unit_tests/vc/CMakeLists.txt b/tests/unit_tests/vc/CMakeLists.txt
index b4f8be4d8..7049c9d37 100644
--- a/tests/unit_tests/vc/CMakeLists.txt
+++ b/tests/unit_tests/vc/CMakeLists.txt
@@ -9,9 +9,13 @@ detray_add_test( vc_array
    "vc_array_actor_chain.cpp"
    "vc_array_annulus2D.cpp"
    "vc_array_container.cpp"
-   "vc_array_coordinates.cpp"
+   "vc_array_coordinate_cartesian2.cpp"
+   "vc_array_coordinate_cartesian3.cpp"
+   "vc_array_coordinate_cylindrical2.cpp"
+   "vc_array_coordinate_cylindrical3.cpp"
+   "vc_array_coordinate_line2.cpp"
+   "vc_array_coordinate_polar2.cpp"
    "vc_array_core.cpp" 
-   "vc_array_covariance_engine.cpp"
    "vc_array_cylinder_intersection.cpp" 
    "vc_array_cylinder.cpp"
    "vc_array_detector.cpp"
@@ -27,11 +31,13 @@ detray_add_test( vc_array
    "vc_array_helix_trajectory.cpp"
    "vc_array_intersection_kernel.cpp"
    "vc_array_line_intersection.cpp"
+   "vc_array_line_stepper.cpp"
    "vc_array_line.cpp"
    "vc_array_mask_store.cpp" 
    "vc_array_materials.cpp"
    "vc_array_navigator.cpp"
    "vc_array_particle_gun.cpp"
+   "vc_array_path_correction.cpp"
    "vc_array_planar_intersection.cpp"
    "vc_array_propagator.cpp"
    "vc_array_rectangle2D.cpp"
diff --git a/tests/unit_tests/vc/vc_array_coordinate_cartesian2.cpp b/tests/unit_tests/vc/vc_array_coordinate_cartesian2.cpp
new file mode 100644
index 000000000..b7930f78d
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_coordinate_cartesian2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/coordinate_cartesian2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_coordinate_cartesian3.cpp b/tests/unit_tests/vc/vc_array_coordinate_cartesian3.cpp
new file mode 100644
index 000000000..7c07bb0b9
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_coordinate_cartesian3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/coordinate_cartesian3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_coordinate_cylindrical2.cpp b/tests/unit_tests/vc/vc_array_coordinate_cylindrical2.cpp
new file mode 100644
index 000000000..cfb3612c1
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_coordinate_cylindrical2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/coordinate_cylindrical2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_coordinate_cylindrical3.cpp b/tests/unit_tests/vc/vc_array_coordinate_cylindrical3.cpp
new file mode 100644
index 000000000..8535f838b
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_coordinate_cylindrical3.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/coordinate_cylindrical3.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_covariance_engine.cpp b/tests/unit_tests/vc/vc_array_coordinate_line2.cpp
similarity index 82%
rename from tests/unit_tests/vc/vc_array_covariance_engine.cpp
rename to tests/unit_tests/vc/vc_array_coordinate_line2.cpp
index f28ae08e7..971485af7 100644
--- a/tests/unit_tests/vc/vc_array_covariance_engine.cpp
+++ b/tests/unit_tests/vc/vc_array_coordinate_line2.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/vc_array_definitions.hpp"
-#include "tests/common/covariance_engine.inl"
\ No newline at end of file
+#include "tests/common/coordinate_line2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_coordinate_polar2.cpp b/tests/unit_tests/vc/vc_array_coordinate_polar2.cpp
new file mode 100644
index 000000000..3cfc52775
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_coordinate_polar2.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/coordinate_polar2.inl"
\ No newline at end of file
diff --git a/tests/unit_tests/vc/vc_array_coordinates.cpp b/tests/unit_tests/vc/vc_array_line_stepper.cpp
similarity index 84%
rename from tests/unit_tests/vc/vc_array_coordinates.cpp
rename to tests/unit_tests/vc/vc_array_line_stepper.cpp
index c9471ae18..eeb43ed24 100644
--- a/tests/unit_tests/vc/vc_array_coordinates.cpp
+++ b/tests/unit_tests/vc/vc_array_line_stepper.cpp
@@ -6,4 +6,4 @@
  */
 
 #include "detray/plugins/algebra/vc_array_definitions.hpp"
-#include "tests/common/coordinates.inl"
\ No newline at end of file
+#include "tests/common/line_stepper.inl"
diff --git a/tests/unit_tests/vc/vc_array_path_correction.cpp b/tests/unit_tests/vc/vc_array_path_correction.cpp
new file mode 100644
index 000000000..fd89ef298
--- /dev/null
+++ b/tests/unit_tests/vc/vc_array_path_correction.cpp
@@ -0,0 +1,9 @@
+/** Detray library, part of the ACTS project (R&D line)
+ *
+ * (c) 2022 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#include "detray/plugins/algebra/vc_array_definitions.hpp"
+#include "tests/common/path_correction.inl"
\ No newline at end of file