diff --git a/.gitignore b/.gitignore
index 4e3d751dc..0d3e2fbc3 100644
--- a/.gitignore
+++ b/.gitignore
@@ -15,4 +15,4 @@
 **/*.d
 **/*.so
 **/*C_ACLiC_dict*
-**/therad_*
\ No newline at end of file
+**/thread_*
\ No newline at end of file
diff --git a/core/include/detray/propagator/actors/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp
index d478f8b70..0b07d15b6 100644
--- a/core/include/detray/propagator/actors/aborters.hpp
+++ b/core/include/detray/propagator/actors/aborters.hpp
@@ -109,7 +109,7 @@ struct next_surface_aborter : actor {
 
         // Abort at the next sensitive surface
         if (navigation.is_on_sensitive() &&
-            stepping._s > abrt_state.min_step_length) {
+            stepping._path_length_per_surface > abrt_state.min_step_length) {
             prop_state._heartbeat &= navigation.abort();
             abrt_state.success = true;
         }
diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp
index 7decf874f..f33a3e837 100644
--- a/core/include/detray/propagator/actors/parameter_resetter.hpp
+++ b/core/include/detray/propagator/actors/parameter_resetter.hpp
@@ -49,7 +49,10 @@ struct parameter_resetter : actor {
                 trf3, mask, stepping._bound_params.vector()));
 
             // Reset the path length
-            stepping._s = 0;
+            stepping._path_length_per_surface = 0;
+
+            // Reset the qop at the surface
+            stepping._qop_i = stepping().qop();
 
             // Reset jacobian coordinate transformation at the current surface
             stepping._jac_to_global = local_coordinate.bound_to_free_jacobian(
@@ -57,6 +60,9 @@ struct parameter_resetter : actor {
 
             // Reset jacobian transport to identity matrix
             matrix_operator().set_identity(stepping._jac_transport);
+
+            // Reset the joint covariance
+            matrix_operator().set_zero(stepping._joint_cov);
         }
     };
 
diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp
index 4f6acbf51..fab159e64 100644
--- a/core/include/detray/propagator/actors/parameter_transporter.hpp
+++ b/core/include/detray/propagator/actors/parameter_transporter.hpp
@@ -12,13 +12,20 @@
 #include "detray/definitions/track_parametrization.hpp"
 #include "detray/geometry/surface.hpp"
 #include "detray/propagator/base_actor.hpp"
+#include "detray/propagator/propagation_config.hpp"
+
+#include <iostream>
 
 namespace detray {
 
 template <typename transform3_t>
 struct parameter_transporter : actor {
 
-    struct state {};
+    struct state {
+        // @TODO remove once we can call stepping config...
+        bool do_scatter = false;
+        bool do_covariance_transport = true;
+    };
 
     /// Mask store visitor
     struct kernel {
@@ -58,7 +65,8 @@ struct parameter_transporter : actor {
                   typename propagator_state_t>
         DETRAY_HOST_DEVICE inline void operator()(
             const mask_group_t& mask_group, const index_t& index,
-            const transform3_type& trf3, propagator_state_t& propagation) {
+            const transform3_type& trf3, propagator_state_t& propagation,
+            state& actor_state) {
 
             // Stepper and Navigator states
             auto& stepping = propagation._stepping;
@@ -119,13 +127,34 @@ struct parameter_transporter : actor {
                           matrix_operator().transpose(stepping._full_jacobian);
             }
 
+            //(void)actor_state;
+            // Calculate multiple scattering term
+            // @TODO: Take stepping::config instead
+            const auto joint_cov = stepping.calculate_ms_covariance(
+                actor_state.do_scatter, actor_state.do_covariance_transport);
+
+            // Add the multiple scattering term
+            // new_cov = new_cov + stepping._joint_cov;
+            new_cov = new_cov + joint_cov;
+
+            /*            
+            std::cout << "Joint Cov" << std::endl;
+            for (int i = 0; i < 6; i++) {
+                for (int j = 0; j < 6; j++) {
+                    std::cout << getter::element(new_cov, i, j) << " ";
+                }
+                std::cout << std::endl;
+            }
+            std::cout << std::endl;
+            */
+
             // Calculate surface-to-surface covariance transport
             stepping._bound_params.set_covariance(new_cov);
         }
     };
 
     template <typename propagator_state_t>
-    DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
+    DETRAY_HOST_DEVICE void operator()(state& actor_state,
                                        propagator_state_t& propagation) const {
         const auto& navigation = propagation._navigation;
 
@@ -141,7 +170,8 @@ struct parameter_transporter : actor {
         // Surface
         const auto sf = navigation.get_surface();
 
-        sf.template visit_mask<kernel>(sf.transform(ctx), propagation);
+        sf.template visit_mask<kernel>(sf.transform(ctx), propagation,
+                                       actor_state);
 
         // Set surface link
         propagation._stepping._bound_params.set_surface_link(sf.barcode());
diff --git a/utils/include/detray/simulation/random_scatterer.hpp b/core/include/detray/propagator/actors/random_scatterer.hpp
similarity index 98%
rename from utils/include/detray/simulation/random_scatterer.hpp
rename to core/include/detray/propagator/actors/random_scatterer.hpp
index 9f2ef0277..17f78f26b 100644
--- a/utils/include/detray/simulation/random_scatterer.hpp
+++ b/core/include/detray/propagator/actors/random_scatterer.hpp
@@ -14,10 +14,10 @@
 #include "detray/geometry/surface.hpp"
 #include "detray/materials/interaction.hpp"
 #include "detray/propagator/base_actor.hpp"
-#include "detray/simulation/landau_distribution.hpp"
-#include "detray/simulation/scattering_helper.hpp"
+#include "detray/propagator/detail/scattering_helper.hpp"
 #include "detray/tracks/bound_track_parameters.hpp"
 #include "detray/utils/axis_rotation.hpp"
+#include "detray/utils/landau_distribution.hpp"
 #include "detray/utils/ranges.hpp"
 #include "detray/utils/unit_vectors.hpp"
 
diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp
index 95596e1ed..3672344f3 100644
--- a/core/include/detray/propagator/base_stepper.hpp
+++ b/core/include/detray/propagator/base_stepper.hpp
@@ -13,6 +13,7 @@
 #include "detray/geometry/surface.hpp"
 #include "detray/propagator/actors/parameter_resetter.hpp"
 #include "detray/propagator/constrained_step.hpp"
+#include "detray/propagator/detail/random_device.hpp"
 #include "detray/propagator/stepping_config.hpp"
 #include "detray/tracks/tracks.hpp"
 
@@ -20,6 +21,11 @@ namespace detray {
 
 namespace stepping {
 
+/// A void random device that does nothing.
+struct void_random_device {
+    void set_seed(const uint_fast64_t /*sd*/) { return; }
+};
+
 /// A void inpector that does nothing.
 ///
 /// Inspectors can be plugged in to understand the current stepper state.
@@ -33,6 +39,7 @@ struct void_inspector {
 
 /// Base stepper implementation
 template <typename transform3_t, typename constraint_t, typename policy_t,
+          typename random_device_t = stepping::void_random_device,
           typename inspector_t = stepping::void_inspector>
 class base_stepper {
 
@@ -69,6 +76,7 @@ class base_stepper {
     struct state {
 
         /// Sets track parameters.
+        /// @TODO: Remove the free track parameter input if possible
         DETRAY_HOST_DEVICE
         state(const free_track_parameters_type &t) : _track(t) {}
 
@@ -86,9 +94,12 @@ class base_stepper {
             sf.template visit_mask<
                 typename parameter_resetter<transform3_t>::kernel>(
                 sf.transform(ctx), *this);
+
+            _qop_i = bound_params.qop();
         }
 
         /// free track parameter
+        /// @TODO: Remove the covariance if possible
         free_track_parameters_type _track;
 
         /// Full jacobian
@@ -106,6 +117,10 @@ class base_stepper {
         /// bound covariance
         bound_track_parameters_type _bound_params;
 
+        /// joint covariance (for multiple scattering)
+        bound_matrix _joint_cov =
+            matrix_operator().template zero<e_bound_size, e_bound_size>();
+
         /// @returns track parameters - const access
         DETRAY_HOST_DEVICE
         free_track_parameters_type &operator()() { return _track; }
@@ -130,16 +145,23 @@ class base_stepper {
 
         /// Track path length from the last surface. It will be reset to 0 when
         /// the track reaches a new surface
-        scalar _s{0.};
+        scalar _path_length_per_surface{0.};
 
         /// Current step size
         scalar _step_size{0.};
 
         /// The particle mass
         scalar_type _mass{105.7f * unit<scalar_type>::MeV};
+
         /// The particle pdg
         int _pdg = 13;  // default muon
 
+        /// qop at the initial surface
+        scalar_type _qop_i{0.f};
+
+        /// Random device
+        random_device_t rand_device;
+
         /// Set new step constraint
         template <step::constraint type = step::constraint::e_actor>
         DETRAY_HOST_DEVICE inline void set_constraint(scalar step_size) {
diff --git a/core/include/detray/propagator/detail/random_device.hpp b/core/include/detray/propagator/detail/random_device.hpp
new file mode 100644
index 000000000..0fbba7c63
--- /dev/null
+++ b/core/include/detray/propagator/detail/random_device.hpp
@@ -0,0 +1,26 @@
+/** Detray plugins library, part of the ACTS project
+ *
+ * (c) 2024 CERN for the benefit of the ACTS project
+ *
+ * Mozilla Public License Version 2.0
+ */
+
+#pragma once
+
+// System include(s).
+#include <random>
+
+namespace detray {
+
+namespace stepping {
+
+struct default_random_device {
+    /// Random generator
+    std::random_device rd{};
+    std::mt19937_64 generator{rd()};
+    void set_seed(const uint_fast64_t sd) { generator.seed(sd); }
+};
+
+}  // namespace stepping
+
+}  // namespace detray
\ No newline at end of file
diff --git a/utils/include/detray/simulation/scattering_helper.hpp b/core/include/detray/propagator/detail/scattering_helper.hpp
similarity index 100%
rename from utils/include/detray/simulation/scattering_helper.hpp
rename to core/include/detray/propagator/detail/scattering_helper.hpp
diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp
index c4b3840a5..497eed101 100644
--- a/core/include/detray/propagator/line_stepper.hpp
+++ b/core/include/detray/propagator/line_stepper.hpp
@@ -12,24 +12,27 @@
 #include "detray/definitions/detail/qualifiers.hpp"
 #include "detray/navigation/policies.hpp"
 #include "detray/propagator/base_stepper.hpp"
+#include "detray/propagator/detail/random_device.hpp"
 
 namespace detray {
 
 /// Straight line stepper implementation
 template <typename transform3_t, typename constraint_t = unconstrained_step,
           typename policy_t = stepper_default_policy,
+          typename random_device_t = stepping::void_random_device,
           typename inspector_t = stepping::void_inspector>
 class line_stepper final
     : public base_stepper<transform3_t, constraint_t, policy_t, inspector_t> {
 
     public:
-    using base_type =
-        base_stepper<transform3_t, constraint_t, policy_t, inspector_t>;
+    using base_type = base_stepper<transform3_t, constraint_t, policy_t,
+                                   random_device_t, inspector_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 bound_matrix = typename bound_track_parameters_type::covariance_type;
     using matrix_operator = typename base_type::matrix_operator;
     using size_type = typename matrix_operator::size_ty;
     using vector3 = typename line_stepper::transform3_type::vector3;
@@ -56,7 +59,7 @@ class line_stepper final
             track.set_pos(track.pos() + track.dir() * this->_step_size);
 
             this->_path_length += this->_step_size;
-            this->_s += this->_step_size;
+            this->_path_length_per_surface += this->_step_size;
         }
 
         DETRAY_HOST_DEVICE
@@ -78,6 +81,12 @@ class line_stepper final
             this->_jac_transport = D * this->_jac_transport;
         }
 
+        DETRAY_HOST_DEVICE
+        inline bound_matrix calculate_ms_covariance(
+            bool /*do_scatter*/, bool /*do_covariance_transport*/) const {
+            return matrix_operator().template zero<6, 6>();
+        }
+
         DETRAY_HOST_DEVICE
         inline vector3 dtds() const { return vector3{0.f, 0.f, 0.f}; }
 
diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp
index 68ca04bc6..dbaa2c8cc 100644
--- a/core/include/detray/propagator/rk_stepper.hpp
+++ b/core/include/detray/propagator/rk_stepper.hpp
@@ -28,14 +28,15 @@ namespace detray {
 template <typename magnetic_field_t, typename transform3_t,
           typename constraint_t = unconstrained_step,
           typename policy_t = stepper_rk_policy,
+          typename random_device_t = stepping::void_random_device,
           typename inspector_t = stepping::void_inspector,
           template <typename, std::size_t> class array_t = darray>
 class rk_stepper final
     : public base_stepper<transform3_t, constraint_t, policy_t, inspector_t> {
 
     public:
-    using base_type =
-        base_stepper<transform3_t, constraint_t, policy_t, inspector_t>;
+    using base_type = base_stepper<transform3_t, constraint_t, policy_t,
+                                   random_device_t, inspector_t>;
 
     using transform3_type = transform3_t;
     using scalar_type = typename transform3_type::scalar_type;
@@ -48,11 +49,13 @@ class rk_stepper final
         typename base_type::free_track_parameters_type;
     using bound_track_parameters_type =
         typename base_type::bound_track_parameters_type;
+    using bound_matrix = typename bound_track_parameters_type::covariance_type;
     using magnetic_field_type = magnetic_field_t;
     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>;
+    using interaction_type = interaction<scalar_type>;
 
     DETRAY_HOST_DEVICE
     rk_stepper() {}
@@ -95,19 +98,25 @@ class rk_stepper final
 
         /// Update the track state by Runge-Kutta-Nystrom integration.
         DETRAY_HOST_DEVICE
-        inline void advance_track();
+        inline void advance_track(
+            const stepping::config<scalar_type>& cfg = {});
 
         /// Update the jacobian transport from free propagation
         DETRAY_HOST_DEVICE
         inline void advance_jacobian(
             const stepping::config<scalar_type>& cfg = {});
 
+        /// Update the covariance from multiple scattering
+        /// @TODO Take stepping::config
+        DETRAY_HOST_DEVICE
+        inline bound_matrix calculate_ms_covariance(
+            bool do_scatter, bool do_covariance_transport) const;
+
         /// evaulate dqopds for a given step size and material
         DETRAY_HOST_DEVICE
         inline scalar_type evaluate_dqopds(
             const std::size_t i, const typename transform3_t::scalar_type h,
-            const scalar dqopds_prev,
-            const detray::stepping::config<scalar_type>& cfg);
+            const scalar dqopds_prev);
 
         /// evaulate dtds for runge kutta stepping
         DETRAY_HOST_DEVICE
diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp
index ae279a4bc..d48b80650 100644
--- a/core/include/detray/propagator/rk_stepper.ipp
+++ b/core/include/detray/propagator/rk_stepper.ipp
@@ -6,14 +6,18 @@
  */
 
 // Project include(s).
+#include "detray/definitions/detail/macros.hpp"
 #include "detray/geometry/detector_volume.hpp"
+#include "detray/materials/interaction.hpp"
+#include "detray/propagator/actors/random_scatterer.hpp"
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE void
 detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   inspector_t, array_t>::state::advance_track() {
+                   random_device_t, inspector_t, array_t>::state::
+    advance_track(const detray::stepping::config<scalar_type>& cfg) {
 
     const auto& sd = this->_step_data;
     const scalar_type h{this->_step_size};
@@ -34,26 +38,73 @@ detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
     track.set_dir(dir);
 
     auto qop = track.qop();
-    if (!(this->_mat == vacuum<scalar_type>())) {
+    const auto mat = this->_mat;
+    if (!(mat == vacuum<scalar_type>())) {
+
+        // only for simulation with valid random device
+        if constexpr (!std::is_same<random_device_t,
+                                    stepping::void_random_device>::value) {
+
+            // Scatter if it is for the simulation
+            if (cfg.do_scatter) {
+                const scalar_type projected_scattering_angle =
+                    interaction_type().compute_multiple_scattering_theta0(
+                        h / mat.X0(), this->_pdg, this->_mass, qop,
+                        track.charge());
+
+                dir = random_scatterer<transform3_t>().scatter(
+                    dir, projected_scattering_angle,
+                    this->rand_device.generator);
+            }
+        }
+
         // Reference: Eq (82) of https://doi.org/10.1016/0029-554X(81)90063-X
-        qop =
-            qop + h_6 * (sd.dqopds[0u] + 2.f * (sd.dqopds[1u] + sd.dqopds[2u]) +
+        if (!cfg.do_scatter) {
+            qop = qop +
+                  h_6 * (sd.dqopds[0u] + 2.f * (sd.dqopds[1u] + sd.dqopds[2u]) +
                          sd.dqopds[3u]);
+        } else {
+            // only for simulation with valid random device
+            if constexpr (!std::is_same<random_device_t,
+                                        stepping::void_random_device>::value) {
+                // Apply mean probable energy loss instead
+                const scalar_type e_loss_mpv =
+                    interaction_type().compute_energy_loss_landau(
+                        h, mat, this->_pdg, this->_mass, qop, track.charge());
+
+                const scalar_type e_loss_sigma =
+                    interaction_type().compute_energy_loss_landau_sigma(
+                        h, mat, this->_pdg, this->_mass, qop, track.charge());
+
+                // Get the new momentum
+                const auto new_mom = random_scatterer<transform3_t>().attenuate(
+                    e_loss_mpv, e_loss_sigma, this->_mass, track.charge() / qop,
+                    this->rand_device.generator);
+
+                qop = track.charge() / new_mom;
+            }
+        }
     }
     track.set_qop(qop);
 
     // Update path length
     this->_path_length += h;
-    this->_s += h;
+    this->_path_length_per_surface += h;
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE void
 detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   inspector_t, array_t>::state::
+                   random_device_t, inspector_t, array_t>::state::
     advance_jacobian(const detray::stepping::config<scalar_type>& cfg) {
+
+    // Immediately exit if it is for simulation
+    if (cfg.do_scatter || !cfg.do_covariance_transport) {
+        return;
+    }
+
     /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the
     /// Jacobian matrix is requires only the calculation of eq. 17 and 18.
     /// Since the terms of eq. 18 are currently 0, this matrix is not needed
@@ -389,15 +440,126 @@ detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
-DETRAY_HOST_DEVICE auto
-detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   inspector_t, array_t>::state::
-    evaluate_dqopds(const std::size_t i,
-                    const typename transform3_t::scalar_type h,
-                    const scalar dqopds_prev,
-                    const detray::stepping::config<scalar_type>& cfg) ->
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
+DETRAY_HOST_DEVICE auto detray::rk_stepper<
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t,
+    array_t>::state::calculate_ms_covariance(bool do_scatter,
+                                             bool do_covariance_transport) const
+    -> bound_matrix {
+
+    // Joint covariance
+    bound_matrix joint_cov = matrix_operator().template zero<6, 6>();
+
+    // Return if there is no material
+    if (this->_mat == vacuum<scalar_type>()) {
+        return joint_cov;
+    }
+
+    // Return if it is for simulation
+    if (do_scatter || do_covariance_transport) {
+        return joint_cov;
+    }
+
+    // Implement the thick scatterer method of Eq (4.103 - 111) of "Pattern
+    // Recognition, Tracking and Vertex Reconstruction in Particle Detectors"
+
+    // Variance per unit length of the projected scattering angle
+    const scalar s1 = interaction_type().compute_multiple_scattering_theta0(
+        1.f / this->_mat.X0(), this->_pdg, this->_mass, this->_qop_i,
+        this->_track.charge());
+    const scalar s2 = interaction_type().compute_multiple_scattering_theta0(
+        1.f / this->_mat.X0(), this->_pdg, this->_mass, this->_track.qop(),
+        this->_track.charge());
+    const scalar variance = s1 * s2;
+
+    // 4X4 Identity matrix
+    matrix_type<4, 4> E = matrix_operator().template zero<4, 4>();
+
+    const scalar L = this->_path_length_per_surface;
+    const scalar half_L2 = L * L / 2.f;
+    const scalar third_L3 = L * L * L / 3.f;
+
+    // Set E (Eq. 4.107). Note that we flip the index of (0,1) <-> (2,3)
+    getter::element(E, 0, 0) = third_L3;
+    getter::element(E, 1, 1) = third_L3;
+    getter::element(E, 2, 2) = L;
+    getter::element(E, 3, 3) = L;
+    getter::element(E, 0, 2) = half_L2;
+    getter::element(E, 1, 3) = half_L2;
+    getter::element(E, 2, 0) = half_L2;
+    getter::element(E, 3, 1) = half_L2;
+    E = variance * E;
+
+    const matrix_type<3, 3> C2G =
+        mat_helper().curvilinear_to_global(this->_track.dir());
+    // Take the bottom-right 2x2 block
+    const matrix_type<2, 2> block =
+        matrix_operator().template block<2, 2>(C2G, 0, 0);
+
+    // Jacobian for (Eq. 4.108).  Note that we flip the index of (0,1) <-> (2,3)
+    matrix_type<4, 4> T = matrix_operator().template identity<4, 4>();
+    matrix_operator().template set_block<2, 2>(T, block, 2, 2);
+
+    E = T * E * matrix_operator().transpose(T);
+
+    // Jacobian d(phi,theta)/d(tx,ty). Note that it is different from
+    // (Eq. 4.109) due to the different coordinate system
+    matrix_type<4, 4> J = matrix_operator().template identity<4, 4>();
+
+    const auto t = this->_track.dir();
+    const scalar_type phi = getter::phi(t);
+    const scalar_type theta = getter::theta(t);
+    const scalar_type sin_phi = math::sin(phi);
+    const scalar_type cos_phi = math::cos(phi);
+    const scalar_type sin_theta = math::sin(theta);
+    const scalar_type cos_theta = math::cos(theta);
+
+    // -sin(phi)/sin(theta)
+    getter::element(J, 2, 2) = -sin_phi / sin_theta;
+    getter::element(J, 2, 3) = -cos_phi / sin_theta;
+    getter::element(J, 3, 2) = -cos_phi * cos_theta;
+    getter::element(J, 3, 3) = -sin_phi * cos_theta;
+
+    E = J * E * matrix_operator().transpose(J);
+
+    // Set the joint covariance
+    matrix_operator().template set_block<4, 4>(joint_cov, E, 0, 0);
+
+    // Add qop variance
+    scalar_type sigma_qop =
+        interaction_type().compute_energy_loss_landau_sigma_QOverP(
+            L, this->_mat, this->_pdg, this->_mass, this->_track.qop(),
+            this->_track.charge());
+
+    // Flip the sign if the step length is negative
+    scalar_type var_qop = sigma_qop * sigma_qop;
+    if (L < 0) {
+        var_qop = -var_qop;
+    }
+
+    /*
+    scalar_type qop1 = this->_qop_i;
+    scalar_type qop2 = this->_track.qop();
+
+    std::cout << L << " " << sigma_qop << " " << -1/qop1 << " " << -1/qop2 << std::endl;
+    */
+
+    getter::element(joint_cov, e_bound_qoverp, e_bound_qoverp) += var_qop;
+
+    return joint_cov;
+}
+
+template <typename magnetic_field_t, typename transform3_t,
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
+DETRAY_HOST_DEVICE auto detray::rk_stepper<
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t,
+    array_t>::state::evaluate_dqopds(const std::size_t i,
+                                     const typename transform3_t::scalar_type h,
+                                     const scalar dqopds_prev) ->
     typename transform3_t::scalar_type {
 
     const auto& track = this->_track;
@@ -409,29 +571,31 @@ detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
         return 0.f;
     } else {
 
-        if (cfg.use_mean_loss) {
-            if (i == 0u) {
-                sd.qop[i] = qop;
-            } else {
-
-                // qop_n is calculated recursively like the direction of
-                // evaluate_dtds.
-                //
-                // https://doi.org/10.1016/0029-554X(81)90063-X says:
-                // "For y  we  have  similar  formulae  as  for x, for y' and
-                // \lambda similar  formulae as for  x'"
-                sd.qop[i] = qop + h * dqopds_prev;
-            }
+        if (i == 0u) {
+            sd.qop[i] = qop;
+        }
+        // For reconstruction with mean energy loss
+        else {
+
+            // qop_n is calculated recursively like the direction of
+            // evaluate_dtds.
+            //
+            // https://doi.org/10.1016/0029-554X(81)90063-X says:
+            // "For y  we  have  similar  formulae  as  for x, for y' and
+            // \lambda similar  formulae as for  x'"
+            sd.qop[i] = qop + h * dqopds_prev;
         }
+
         return this->dqopds(sd.qop[i]);
     }
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE auto detray::rk_stepper<
-    magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t,
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t,
     array_t>::state::evaluate_dtds(const vector3& b_field, const std::size_t i,
                                    const typename transform3_t::scalar_type h,
                                    const vector3& dtds_prev,
@@ -453,11 +617,11 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper<
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE auto detray::rk_stepper<
-    magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t,
-    array_t>::state::evaluate_field_gradient(const vector3& pos)
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t, array_t>::state::evaluate_field_gradient(const vector3& pos)
     -> matrix_type<3, 3> {
 
     matrix_type<3, 3> dBdr = matrix_operator().template zero<3, 3>();
@@ -495,22 +659,22 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper<
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE auto
 detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   inspector_t, array_t>::state::dqopds() const ->
-    typename transform3_t::scalar_type {
+                   random_device_t, inspector_t, array_t>::state::dqopds() const
+    -> typename transform3_t::scalar_type {
     return this->_step_data.dqopds[3u];
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
-DETRAY_HOST_DEVICE auto
-detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
-                   inspector_t, array_t>::state::dqopds(const scalar_type qop)
-    const -> typename transform3_t::scalar_type {
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
+DETRAY_HOST_DEVICE auto detray::rk_stepper<
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t, array_t>::state::dqopds(const scalar_type qop) const ->
+    typename transform3_t::scalar_type {
 
     const auto& mat = this->_mat;
 
@@ -539,11 +703,11 @@ detray::rk_stepper<magnetic_field_t, transform3_t, constraint_t, policy_t,
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 DETRAY_HOST_DEVICE auto detray::rk_stepper<
-    magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t,
-    array_t>::state::d2qopdsdqop(const scalar_type qop) const ->
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t, array_t>::state::d2qopdsdqop(const scalar_type qop) const ->
     typename transform3_t::scalar_type {
 
     using scalar_t = typename transform3_t::scalar_type;
@@ -582,13 +746,14 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper<
 }
 
 template <typename magnetic_field_t, typename transform3_t,
-          typename constraint_t, typename policy_t, typename inspector_t,
-          template <typename, std::size_t> class array_t>
+          typename constraint_t, typename policy_t, typename random_device_t,
+          typename inspector_t, template <typename, std::size_t> class array_t>
 template <typename propagation_state_t>
 DETRAY_HOST_DEVICE bool detray::rk_stepper<
-    magnetic_field_t, transform3_t, constraint_t, policy_t, inspector_t,
-    array_t>::step(propagation_state_t& propagation,
-                   const detray::stepping::config<scalar_type>& cfg) {
+    magnetic_field_t, transform3_t, constraint_t, policy_t, random_device_t,
+    inspector_t, array_t>::step(propagation_state_t& propagation,
+                                const detray::stepping::config<scalar_type>&
+                                    cfg) {
 
     // Get stepper and navigator states
     state& stepping = propagation._stepping;
@@ -611,7 +776,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper<
 
     // qop should be recalcuated at every point
     // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X
-    sd.dqopds[0u] = stepping.evaluate_dqopds(0u, 0.f, 0.f, cfg);
+    sd.dqopds[0u] = stepping.evaluate_dqopds(0u, 0.f, 0.f);
     sd.dtds[0u] = stepping.evaluate_dtds(sd.b_first, 0u, 0.f,
                                          vector3{0.f, 0.f, 0.f}, sd.qop[0u]);
 
@@ -630,16 +795,14 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper<
         sd.b_middle[1] = bvec1[1];
         sd.b_middle[2] = bvec1[2];
 
-        sd.dqopds[1u] =
-            stepping.evaluate_dqopds(1u, half_h, sd.dqopds[0u], cfg);
+        sd.dqopds[1u] = stepping.evaluate_dqopds(1u, half_h, sd.dqopds[0u]);
         sd.dtds[1u] = stepping.evaluate_dtds(sd.b_middle, 1u, half_h,
                                              sd.dtds[0u], sd.qop[1u]);
 
         // Third Runge-Kutta point
         // qop should be recalcuated at every point
         // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X
-        sd.dqopds[2u] =
-            stepping.evaluate_dqopds(2u, half_h, sd.dqopds[1u], cfg);
+        sd.dqopds[2u] = stepping.evaluate_dqopds(2u, half_h, sd.dqopds[1u]);
         sd.dtds[2u] = stepping.evaluate_dtds(sd.b_middle, 2u, half_h,
                                              sd.dtds[1u], sd.qop[2u]);
 
@@ -652,7 +815,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper<
         sd.b_last[1] = bvec2[1];
         sd.b_last[2] = bvec2[2];
 
-        sd.dqopds[3u] = stepping.evaluate_dqopds(3u, h, sd.dqopds[2u], cfg);
+        sd.dqopds[3u] = stepping.evaluate_dqopds(3u, h, sd.dqopds[2u]);
         sd.dtds[3u] =
             stepping.evaluate_dtds(sd.b_last, 3u, h, sd.dtds[2u], sd.qop[3u]);
 
@@ -728,7 +891,7 @@ DETRAY_HOST_DEVICE bool detray::rk_stepper<
     }
 
     // Advance track state
-    stepping.advance_track();
+    stepping.advance_track(cfg);
 
     // Advance jacobian transport
     stepping.advance_jacobian(cfg);
diff --git a/core/include/detray/propagator/stepping_config.hpp b/core/include/detray/propagator/stepping_config.hpp
index 579f1e291..85e28c2e5 100644
--- a/core/include/detray/propagator/stepping_config.hpp
+++ b/core/include/detray/propagator/stepping_config.hpp
@@ -12,6 +12,7 @@
 
 // System include(s).
 #include <limits>
+#include <random>
 
 namespace detray::stepping {
 
@@ -34,13 +35,22 @@ struct config {
     scalar_t path_limit{5.f * unit<scalar_t>::m};
     /// Maximum number of Runge-Kutta step trials
     std::size_t max_rk_updates{10000u};
+    /*
+    /// Scatter trajectory for simulation
+    bool do_scattering{false};
     /// Use mean energy loss (Bethe)
     /// if false, most probable energy loss (Landau) will be used
     bool use_mean_loss{true};
+    */
     /// Use eloss gradient in error propagation
     bool use_eloss_gradient{false};
     /// Use b field gradient in error propagation
     bool use_field_gradient{false};
+
+    /// Do the scattering (for simulation)
+    bool do_scatter{false};
+    /// Do covariance transport (for Reconstruction)
+    bool do_covariance_transport{true};
 };
 
 }  // namespace detray::stepping
diff --git a/core/include/detray/utils/axis_rotation.hpp b/core/include/detray/utils/axis_rotation.hpp
index 5e40817ce..c45b24373 100644
--- a/core/include/detray/utils/axis_rotation.hpp
+++ b/core/include/detray/utils/axis_rotation.hpp
@@ -41,6 +41,8 @@ struct axis_rotation {
 
         scalar_type cos_theta{math::cos(theta)};
 
+        // Rodrigues' rotation_formula
+        // https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
         matrix_type<3, 3> I = matrix_operator().template identity<3, 3>();
         matrix_type<3, 3> axis_cross = mat_helper().cross_matrix(U);
         matrix_type<3, 3> axis_outer = mat_helper().outer_product(U, U);
diff --git a/utils/include/detray/simulation/landau_distribution.hpp b/core/include/detray/utils/landau_distribution.hpp
similarity index 100%
rename from utils/include/detray/simulation/landau_distribution.hpp
rename to core/include/detray/utils/landau_distribution.hpp
diff --git a/core/include/detray/utils/matrix_helper.hpp b/core/include/detray/utils/matrix_helper.hpp
index 5038ab01e..b7c3c787c 100644
--- a/core/include/detray/utils/matrix_helper.hpp
+++ b/core/include/detray/utils/matrix_helper.hpp
@@ -67,6 +67,12 @@ struct matrix_helper {
     }
 
     /// Cross product matrix
+    ///
+    /// For a given vector (v) = (x,y,z), return the following matrix:
+    ///
+    /// [ 0  -z   y ]
+    /// [ z   0  -x ]
+    /// [-y   x   0 ]
     DETRAY_HOST_DEVICE
     inline matrix_type<3, 3> cross_matrix(const vector3& v) const {
         matrix_type<3, 3> ret;
@@ -84,6 +90,8 @@ struct matrix_helper {
     }
 
     /// Outer product operation
+    ///
+    /// @return 3x3 matrix M = v * v^T where v is a 3-elements column matrix
     DETRAY_HOST_DEVICE
     inline matrix_type<3, 3> outer_product(const vector3& v1,
                                            const vector3& v2) const {
@@ -101,6 +109,8 @@ struct matrix_helper {
     }
 
     /// Cholesky decompose
+    ///
+    /// @return the lower triangle matrix (L) that satisfies M = L * L^T
     template <size_type N>
     DETRAY_HOST_DEVICE inline matrix_type<N, N> cholesky_decompose(
         const matrix_type<N, N>& mat) const {
@@ -127,6 +137,93 @@ struct matrix_helper {
 
         return L;
     }
+
+    /// @return the rotation matrix (R) that converts the local curvilinear
+    /// cartesian to the global cartesian
+    ///
+    /// The local z-axis (w) is defined as the given vector. If we say that u
+    /// and v are the local x and y axis, R is defined as follows:
+    ///
+    ///     [ u₁  v₁  w₁ ]
+    /// R = [ u₂  v₂  w₂ ]
+    ///     [ u₃  v₃  w₃ ]
+    ///
+    /// One way to obtain R is solving Rodrigues' rotation_formula (This is the
+    /// formula implemented in axis_rotation.hpp):
+    ///
+    /// cos_α I + sin_α cross_matrix(k) + (1 - cos_α) * outer_product(k)
+    ///
+    /// k: normalization of (0,0,1) X w
+    /// α: the angle required to rotate (0,0,1) to w around the axis of k
+    /// cross_matrix: function defined above
+    /// outer_product: function defined_above
+    ///
+    ///
+    /// In a nutshell, the rest of R elements is calculated as:
+    ///
+    /// u₁ = (w₁ w₁ w₃ + w₂ w₂) / ( 1 - w₃ w₃ )
+    /// u₂ = -w₁w₂/(1 + w₃)
+    /// u₃ = -w₁
+    /// v₁ = -w₁w₂/(1 + w₃)
+    /// v₂ = (w₁ w₁ + w₂ w₂ w₃) / ( 1 - w₃ w₃ )
+    /// v₃ = -w₂
+    DETRAY_HOST_DEVICE
+    inline matrix_type<3, 3> curvilinear_to_global(const vector3& w) const {
+
+        matrix_type<3, 3> R = matrix_operator().template zero<3, 3>();
+
+        const scalar_type w1w2 = w[0] * w[1];
+
+        // Set u
+        getter::element(R, 0u, 0u) =
+            (w[0] * w[0] * w[2] + w[1] * w[1]) / (1 - w[2] * w[2]);
+        getter::element(R, 1u, 0u) = -w1w2 / (1 + w[2]);
+        getter::element(R, 2u, 0u) = -w[0];
+
+        // Set v
+        getter::element(R, 0u, 1u) = getter::element(R, 1u, 0u);
+        getter::element(R, 1u, 1u) =
+            (w[0] * w[0] + w[1] * w[1] * w[2]) / (1 - w[2] * w[2]);
+        getter::element(R, 2u, 1u) = -w[1];
+
+        // Set w
+        getter::element(R, 0u, 2u) = w[0u];
+        getter::element(R, 1u, 2u) = w[1u];
+        getter::element(R, 2u, 2u) = w[2u];
+
+        return R;
+    }
+
+    /// @return the rotation matrix (R) that converts the global cartesian to
+    /// the local curvilinear cartesian
+    ///
+    /// It is the transpose of curvlinear_to_global
+    DETRAY_HOST_DEVICE
+    inline matrix_type<3, 3> global_to_curvilinear(const vector3& w) {
+
+        matrix_type<3, 3> R = matrix_operator().template zero<3, 3>();
+
+        const scalar_type w1w2 = w[0] * w[1];
+
+        // Set u
+        getter::element(R, 0u, 0u) =
+            (w[0] * w[0] * w[2] + w[1] * w[1]) / (1 - w[2] * w[2]);
+        getter::element(R, 1u, 0u) = -w1w2 / (1 + w[2]);
+        getter::element(R, 2u, 0u) = w[0];
+
+        // Set v
+        getter::element(R, 0u, 1u) = getter::element(R, 1u, 0u);
+        getter::element(R, 1u, 1u) =
+            (w[0] * w[0] + w[1] * w[1] * w[2]) / (1 - w[2] * w[2]);
+        getter::element(R, 2u, 1u) = w[1];
+
+        // Set w
+        getter::element(R, 0u, 2u) = -w[0u];
+        getter::element(R, 1u, 2u) = -w[1u];
+        getter::element(R, 2u, 2u) = w[2u];
+
+        return R;
+    }
 };
 
 }  // namespace detray
diff --git a/tests/include/detray/test/helix_navigation.hpp b/tests/include/detray/test/helix_navigation.hpp
index bf8c1d7f8..3c785c55b 100644
--- a/tests/include/detray/test/helix_navigation.hpp
+++ b/tests/include/detray/test/helix_navigation.hpp
@@ -107,9 +107,10 @@ class helix_navigation : public test::fixture_base<> {
         using navigator_t = navigator<detector_t, inspector_t, intersection_t>;
         // Runge-Kutta stepper
         using bfield_t = bfield::const_field_t;
-        using stepper_t = rk_stepper<typename bfield_t::view_t, transform3_t,
-                                     unconstrained_step, stepper_rk_policy,
-                                     stepping::print_inspector>;
+        using stepper_t =
+            rk_stepper<typename bfield_t::view_t, transform3_t,
+                       unconstrained_step, stepper_rk_policy,
+                       stepping::void_random_device, stepping::print_inspector>;
         // Propagator with pathlimit aborter
         using actor_chain_t = actor_chain<dtuple, pathlimit_aborter>;
         using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
diff --git a/tests/include/detray/test/straight_line_navigation.hpp b/tests/include/detray/test/straight_line_navigation.hpp
index dc59d200b..d49c77fc8 100644
--- a/tests/include/detray/test/straight_line_navigation.hpp
+++ b/tests/include/detray/test/straight_line_navigation.hpp
@@ -110,7 +110,8 @@ class straight_line_navigation : public test::fixture_base<> {
         //  Line stepper
         using stepper_t =
             line_stepper<transform3_t, unconstrained_step,
-                         stepper_default_policy, stepping::print_inspector>;
+                         stepper_default_policy, stepping::void_random_device,
+                         stepping::print_inspector>;
         // Propagator with pathlimit aborter
         using actor_chain_t = actor_chain<dtuple, pathlimit_aborter>;
         using propagator_t = propagator<stepper_t, navigator_t, actor_chain_t>;
diff --git a/tests/integration_tests/cpu/material/material_interaction.cpp b/tests/integration_tests/cpu/material/material_interaction.cpp
index f5ba6f374..b2252cd36 100644
--- a/tests/integration_tests/cpu/material/material_interaction.cpp
+++ b/tests/integration_tests/cpu/material/material_interaction.cpp
@@ -22,10 +22,10 @@
 #include "detray/propagator/actors/parameter_resetter.hpp"
 #include "detray/propagator/actors/parameter_transporter.hpp"
 #include "detray/propagator/actors/pointwise_material_interactor.hpp"
+#include "detray/propagator/actors/random_scatterer.hpp"
 #include "detray/propagator/line_stepper.hpp"
 #include "detray/propagator/propagator.hpp"
 #include "detray/propagator/rk_stepper.hpp"
-#include "detray/simulation/random_scatterer.hpp"
 #include "detray/test/types.hpp"
 #include "detray/utils/inspectors.hpp"
 #include "detray/utils/statistics.hpp"
diff --git a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp
index 4faa1681c..06dd8c267 100644
--- a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp
+++ b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp
@@ -13,6 +13,7 @@
 #include "detray/navigation/navigator.hpp"
 #include "detray/propagator/actors/parameter_resetter.hpp"
 #include "detray/propagator/actors/parameter_transporter.hpp"
+#include "detray/propagator/detail/random_device.hpp"
 #include "detray/propagator/propagator.hpp"
 #include "detray/propagator/rk_stepper.hpp"
 #include "detray/simulation/event_generator/track_generators.hpp"
@@ -340,12 +341,13 @@ struct bound_getter : actor {
 /// Numerically integrate the jacobian
 template <typename propagator_t, typename field_t>
 bound_getter<transform3_type>::state evaluate_bound_param(
-    const scalar detector_length,
+    const std::size_t mc_seed, const scalar detector_length,
     const bound_track_parameters<transform3_type>& initial_param,
     const typename propagator_t::detector_type& det, const field_t& field,
     const scalar overstep_tolerance, const scalar on_surface_tolerance,
     const scalar rk_tolerance, const scalar constraint_step,
-    bool use_field_gradient, bool do_inspect = false) {
+    bool use_field_gradient, bool do_scatter, bool do_covariance_transport,
+    bool do_inspect = false) {
 
     // Propagator is built from the stepper and navigator
     propagation::config<scalar> cfg{};
@@ -354,10 +356,13 @@ bound_getter<transform3_type>::state evaluate_bound_param(
     cfg.stepping.rk_error_tol = rk_tolerance;
     cfg.stepping.use_eloss_gradient = true;
     cfg.stepping.use_field_gradient = use_field_gradient;
+    cfg.stepping.do_scatter = do_scatter;
+    cfg.stepping.do_covariance_transport = do_covariance_transport;
     propagator_t p(cfg);
 
     // Actor states
-    parameter_transporter<transform3_type>::state transporter_state{};
+    parameter_transporter<transform3_type>::state transporter_state{
+        cfg.stepping.do_scatter, cfg.stepping.do_covariance_transport};
     bound_getter<transform3_type>::state bound_getter_state{};
     bound_getter_state.m_min_path_length = detector_length * 0.75f;
     parameter_resetter<transform3_type>::state resetter_state{};
@@ -366,6 +371,7 @@ bound_getter<transform3_type>::state evaluate_bound_param(
 
     // Init propagator states for the reference track
     typename propagator_t::state state(initial_param, field, det);
+    state._stepping.rand_device.set_seed(mc_seed);
 
     // Run the propagation for the reference track
     state.do_debug = do_inspect;
@@ -395,6 +401,9 @@ get_displaced_bound_vector(
     cfg.navigation.overstep_tolerance = overstep_tolerance;
     cfg.navigation.on_surface_tolerance = on_surface_tolerance;
     cfg.stepping.rk_error_tol = rk_tolerance;
+    // we don't have to scatering or covariance transport
+    cfg.stepping.do_scatter = false;
+    cfg.stepping.do_covariance_transport = false;
 
     // Propagator is built from the stepper and navigator
     propagator_t p(cfg);
@@ -408,7 +417,8 @@ get_displaced_bound_vector(
     typename propagator_t::state dstate(dparam, field, det);
 
     // Actor states
-    parameter_transporter<transform3_type>::state transporter_state{};
+    parameter_transporter<transform3_type>::state transporter_state{
+        cfg.stepping.do_scatter, cfg.stepping.do_covariance_transport};
     parameter_resetter<transform3_type>::state resetter_state{};
     bound_getter<transform3_type>::state bound_getter_state{};
     bound_getter_state.m_min_path_length = detector_length * 0.75f;
@@ -609,8 +619,8 @@ bound_track_parameters<transform3_type> get_initial_parameter(
 
 template <typename propagator_t, typename field_t>
 void evaluate_jacobian_difference(
-    const unsigned int trk_count, typename propagator_t::detector_type& det,
-    const scalar detector_length,
+    const std::size_t mc_seed, const unsigned int trk_count,
+    typename propagator_t::detector_type& det, const scalar detector_length,
     const bound_track_parameters<transform3_type>& track, const field_t& field,
     const material<scalar> volume_mat, const scalar overstep_tolerance,
     const scalar on_surface_tolerance, const scalar rk_tolerance,
@@ -625,10 +635,15 @@ void evaluate_jacobian_difference(
 
     det.volumes()[0u].set_material(volume_mat);
 
+    const bool do_scatter = false;
+    const bool do_covariance_transport = true;
+
+    // No scattering for jacobian comparsion study
+    // But perform covariance transport
     auto bound_getter = evaluate_bound_param<propagator_t, field_t>(
-        detector_length, track, det, field, overstep_tolerance,
+        mc_seed, detector_length, track, det, field, overstep_tolerance,
         on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient,
-        do_inspect);
+        do_scatter, do_covariance_transport, do_inspect);
 
     const auto reference_param = bound_getter.m_param_departure;
     const auto final_param = bound_getter.m_param_destination;
@@ -735,6 +750,7 @@ void evaluate_jacobian_difference(
 
 template <typename propagator_t, typename field_t>
 void evaluate_covariance_transport(
+    const bool include_multiple_scattering, const std::size_t mc_seed,
     const unsigned int trk_count, typename propagator_t::detector_type& det,
     const scalar detector_length,
     const bound_track_parameters<transform3_type>& track, const field_t& field,
@@ -754,10 +770,14 @@ void evaluate_covariance_transport(
 
     track_copy.set_covariance(ini_cov);
 
+    bool do_scatter = false;
+    bool do_covariance_transport = true;
+
+    // Reference track (Without Scattering but with covariance transport)
     auto bound_getter = evaluate_bound_param<propagator_t, field_t>(
-        detector_length, track_copy, det, field, overstep_tolerance,
-        on_surface_tolerance, rk_tolerance, constraint_step,
-        use_field_gradient);
+        mc_seed, detector_length, track_copy, det, field, overstep_tolerance,
+        on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient,
+        do_scatter, do_covariance_transport);
 
     const auto reference_param = bound_getter.m_param_departure;
     const auto ini_vec = reference_param.vector();
@@ -793,10 +813,18 @@ void evaluate_covariance_transport(
     auto smeared_track = track_copy;
     smeared_track.set_vector(smeared_ini_vec);
 
+    if (include_multiple_scattering) {
+        do_scatter = true;
+    } else {
+        do_scatter = false;
+    }
+    do_covariance_transport = false;
+
+    // Smeared track (With scattering and no covariance transport)
     auto smeared_bound_getter = evaluate_bound_param<propagator_t, field_t>(
-        detector_length, smeared_track, det, field, overstep_tolerance,
-        on_surface_tolerance, rk_tolerance, constraint_step,
-        use_field_gradient);
+        mc_seed, detector_length, smeared_track, det, field, overstep_tolerance,
+        on_surface_tolerance, rk_tolerance, constraint_step, use_field_gradient,
+        do_scatter, do_covariance_transport);
 
     // Get smeared final bound vector
     bound_vector_type smeared_fin_vec =
@@ -1324,6 +1352,9 @@ int main(int argc, char** argv) {
                        "Set log10(max_rk_tolerance_in_mm)");
     desc.add_options()("mc-seed", po::value<std::size_t>()->default_value(0u),
                        "Monte-Carlo seed");
+    desc.add_options()("include-multiple-scattering",
+                       po::value<bool>()->default_value(false),
+                       "Include multiple scattering");
     desc.add_options()("verbose-level", po::value<int>()->default_value(1),
                        "Verbose level");
 
@@ -1360,6 +1391,8 @@ int main(int argc, char** argv) {
     const scalar log10_max_rk_tolerance =
         vm["log10-max-rk-tolerance-mm"].as<scalar>() * unit<scalar>::mm;
     const std::size_t mc_seed = vm["mc-seed"].as<std::size_t>();
+    const bool include_multiple_scattering =
+        vm["include-multiple-scattering"].as<bool>();
     const int verbose_lvl = vm["verbose-level"].as<int>();
 
     std::vector<scalar> log10_tols;
@@ -1506,10 +1539,10 @@ int main(int argc, char** argv) {
     // Stepper types
     using const_field_stepper_t =
         rk_stepper<const_bfield_t::view_t, transform3_type, constrained_step<>,
-                   stepper_default_policy>;
+                   stepper_default_policy, stepping::default_random_device>;
     using inhom_field_stepper_t =
         rk_stepper<inhom_bfield_t::view_t, transform3_type, constrained_step<>,
-                   stepper_default_policy>;
+                   stepper_default_policy, stepping::default_random_device>;
 
     // Make four propagators for each case
     using const_field_rect_propagator_t =
@@ -1653,11 +1686,11 @@ int main(int argc, char** argv) {
 
                     // Rectangle Inhomogeneous field with Material
                     evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
-                        track_count, rect_det, detector_length, rect_bparam,
-                        inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                        std::pow(10.f, log10_tols[i]), constraint_step_size,
-                        h_sizes_rect, rect_files[i], ref_rel_diff, true,
-                        do_inspect);
+                        mc_seed, track_count, rect_det, detector_length,
+                        rect_bparam, inhom_bfield, volume_mat, overstep_tol,
+                        on_surface_tol, std::pow(10.f, log10_tols[i]),
+                        constraint_step_size, h_sizes_rect, rect_files[i],
+                        ref_rel_diff, true, do_inspect);
 
                     dqopdqop_rel_diffs_rect[i].push_back(ref_rel_diff);
                 }
@@ -1672,8 +1705,8 @@ int main(int argc, char** argv) {
 
                 // Rect Const field
                 evaluate_jacobian_difference<const_field_rect_propagator_t>(
-                    track_count, rect_det, detector_length, rect_bparam,
-                    const_bfield, vacuum<scalar>(), overstep_tol,
+                    mc_seed, track_count, rect_det, detector_length,
+                    rect_bparam, const_bfield, vacuum<scalar>(), overstep_tol,
                     on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect,
                     const_rect_file, ref_rel_diff);
 
@@ -1688,24 +1721,25 @@ int main(int argc, char** argv) {
 
                 // Rect Inhomogeneous field
                 evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
-                    track_count, rect_det, detector_length, rect_bparam,
-                    inhom_bfield, vacuum<scalar>(), overstep_tol,
+                    mc_seed, track_count, rect_det, detector_length,
+                    rect_bparam, inhom_bfield, vacuum<scalar>(), overstep_tol,
                     on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect,
                     inhom_rect_file, ref_rel_diff);
 
                 // Rectangle Inhomogeneous field with Material
                 evaluate_jacobian_difference<inhom_field_rect_propagator_t>(
-                    track_count, rect_det, detector_length, rect_bparam,
-                    inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                    rk_tol, constraint_step_size, h_sizes_rect,
+                    mc_seed, track_count, rect_det, detector_length,
+                    rect_bparam, inhom_bfield, volume_mat, overstep_tol,
+                    on_surface_tol, rk_tol, constraint_step_size, h_sizes_rect,
                     inhom_rect_material_file, ref_rel_diff);
 
                 // Rectangle Inhomogeneous field with Material (Covariance
                 // transport)
                 evaluate_covariance_transport<inhom_field_rect_propagator_t>(
-                    track_count, rect_det, detector_length, rect_bparam,
-                    inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                    rk_tol, constraint_step_size, rect_cov_transport_file);
+                    include_multiple_scattering, mc_seed, track_count, rect_det,
+                    detector_length, rect_bparam, inhom_bfield, volume_mat,
+                    overstep_tol, on_surface_tol, rk_tol, constraint_step_size,
+                    rect_cov_transport_file);
             }
         }
 
@@ -1731,11 +1765,11 @@ int main(int argc, char** argv) {
                 for (std::size_t i = 0u; i < log10_tols.size(); i++) {
                     // Wire Inhomogeneous field with Material
                     evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
-                        track_count, wire_det, detector_length, wire_bparam,
-                        inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                        std::pow(10.f, log10_tols[i]), constraint_step_size,
-                        h_sizes_wire, wire_files[i], ref_rel_diff, true,
-                        do_inspect);
+                        mc_seed, track_count, wire_det, detector_length,
+                        wire_bparam, inhom_bfield, volume_mat, overstep_tol,
+                        on_surface_tol, std::pow(10.f, log10_tols[i]),
+                        constraint_step_size, h_sizes_wire, wire_files[i],
+                        ref_rel_diff, true, do_inspect);
 
                     dqopdqop_rel_diffs_wire[i].push_back(ref_rel_diff);
                 }
@@ -1750,8 +1784,8 @@ int main(int argc, char** argv) {
 
                 // Wire Const field
                 evaluate_jacobian_difference<const_field_wire_propagator_t>(
-                    track_count, wire_det, detector_length, wire_bparam,
-                    const_bfield, vacuum<scalar>(), overstep_tol,
+                    mc_seed, track_count, wire_det, detector_length,
+                    wire_bparam, const_bfield, vacuum<scalar>(), overstep_tol,
                     on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire,
                     const_wire_file, ref_rel_diff);
 
@@ -1766,23 +1800,24 @@ int main(int argc, char** argv) {
 
                 // Wire Inhomogeneous field
                 evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
-                    track_count, wire_det, detector_length, wire_bparam,
-                    inhom_bfield, vacuum<scalar>(), overstep_tol,
+                    mc_seed, track_count, wire_det, detector_length,
+                    wire_bparam, inhom_bfield, vacuum<scalar>(), overstep_tol,
                     on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire,
                     inhom_wire_file, ref_rel_diff);
 
                 // Wire Inhomogeneous field with Material
                 evaluate_jacobian_difference<inhom_field_wire_propagator_t>(
-                    track_count, wire_det, detector_length, wire_bparam,
-                    inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                    rk_tol, constraint_step_size, h_sizes_wire,
+                    mc_seed, track_count, wire_det, detector_length,
+                    wire_bparam, inhom_bfield, volume_mat, overstep_tol,
+                    on_surface_tol, rk_tol, constraint_step_size, h_sizes_wire,
                     inhom_wire_material_file, ref_rel_diff);
 
                 // Wire Inhomogeneous field with Material (Covariance transport)
                 evaluate_covariance_transport<inhom_field_wire_propagator_t>(
-                    track_count, wire_det, detector_length, wire_bparam,
-                    inhom_bfield, volume_mat, overstep_tol, on_surface_tol,
-                    rk_tol, constraint_step_size, wire_cov_transport_file);
+                    include_multiple_scattering, mc_seed, track_count, wire_det,
+                    detector_length, wire_bparam, inhom_bfield, volume_mat,
+                    overstep_tol, on_surface_tol, rk_tol, constraint_step_size,
+                    wire_cov_transport_file);
             }
         }
     }
diff --git a/tests/integration_tests/cpu/propagator/propagator.cpp b/tests/integration_tests/cpu/propagator/propagator.cpp
index 43bf28e15..0e0c1f871 100644
--- a/tests/integration_tests/cpu/propagator/propagator.cpp
+++ b/tests/integration_tests/cpu/propagator/propagator.cpp
@@ -71,7 +71,8 @@ struct helix_inspector : actor {
         }
 
         // Nothing has happened yet (first call of actor chain)
-        if (stepping.path_length() < tol || stepping._s < tol) {
+        if (stepping.path_length() < tol ||
+            stepping._path_length_per_surface < tol) {
             return;
         }
 
@@ -94,21 +95,21 @@ struct helix_inspector : actor {
 
         detail::helix<transform3> hlx(free_params, &b);
 
-        const auto true_pos = hlx(stepping._s);
+        const auto true_pos = hlx(stepping._path_length_per_surface);
 
-        const point3 relative_error{1.f / stepping._s *
+        const point3 relative_error{1.f / stepping._path_length_per_surface *
                                     (stepping().pos() - true_pos)};
 
         ASSERT_NEAR(getter::norm(relative_error), 0.f, tol);
 
-        auto true_J = hlx.jacobian(stepping._s);
+        auto true_J = hlx.jacobian(stepping._path_length_per_surface);
 
         for (unsigned int i = 0u; i < e_free_size; i++) {
             for (unsigned int j = 0u; j < e_free_size; j++) {
                 ASSERT_NEAR(
                     matrix_operator().element(stepping._jac_transport, i, j),
                     matrix_operator().element(true_J, i, j),
-                    stepping._s * tol * 10.f);
+                    stepping._path_length_per_surface * tol * 10.f);
             }
         }
     }
diff --git a/tests/scripts/run_jacobian_validation.sh b/tests/scripts/run_jacobian_validation.sh
index d3542471d..a1954943d 100644
--- a/tests/scripts/run_jacobian_validation.sh
+++ b/tests/scripts/run_jacobian_validation.sh
@@ -16,6 +16,9 @@ n_tracks_per_thread=1000
 log10_min_rk_tol=-6
 log10_max_rk_tol=2
 
+# log10 rk tolerance for covariance transport
+log10_rk_tol=-4
+
 # Helix intersector tolerance [log10 in mm]
 log10_helix_tol=-3
 # Surface tolerance [log10 in mm]
@@ -24,13 +27,19 @@ log10_on_surface_tol=-3
 # Monte-Carlo seed
 mc_seed=0
 
+# include multiple scattering
+include_multiple_scattering=false
+
+# Skip the first phase
+skip_first_phase=false
+
 # Skip the second phase
 skip_second_phase=false
 
 # Verbose level
 verbose_level=1
 
-while getopts "hd:n:t:p:q:i:s:r:v:" arg; do
+while getopts "hd:n:t:p:q:c:i:s:m:f:r:v:" arg; do
     case $arg in
         h)
             echo ""
@@ -42,8 +51,11 @@ while getopts "hd:n:t:p:q:i:s:r:v:" arg; do
             echo "-t <Number of tracks per thread>"
             echo "-p <log10(min_rk_error_tolerance_in_mm)>"
             echo "-q <log10(max_rk_error_tolerance_in_mm)>"
+            echo "-c <log10(rk_error_tolerance_in_mm_for_covariance_transport)>"
             echo "-i <log10(intersection_tolerance_in_mm)>"
             echo "-s <Monte-Carlo seed>"
+            echo "-m <Include myltiple scattering>"
+            echo "-f <Skip the first phase>"
             echo "-r <Skip the second phase>"
             echo "-v <Verbose level>"
             echo ""
@@ -69,6 +81,10 @@ while getopts "hd:n:t:p:q:i:s:r:v:" arg; do
             log10_max_rk_tol=$OPTARG
             echo "log10(max_rk_error_tolerance_in_mm): ${log10_max_rk_tol}"
         ;;
+        c)
+            log10_rk_tol=$OPTARG
+            echo "log10(rk_error_tolerance_in_mm_for_covariance_transport): ${log10_rk_tol}"
+        ;;
         i)
             log10_helix_tol=$OPTARG
             log10_on_surface_tol=$OPTARG
@@ -78,6 +94,14 @@ while getopts "hd:n:t:p:q:i:s:r:v:" arg; do
             mc_seed=$OPTARG
             echo "Monte-Carlo seed: ${mc_seed}"
         ;;
+        m)
+            include_multiple_scattering=$OPTARG
+            echo "Include multiple scattering in the second phase: ${include_multiple_scattering}"
+        ;;
+        f)
+            skip_first_phase=$OPTARG
+            echo "Skip the first phase: ${skip_first_phase}"
+        ;;
         r)
             skip_second_phase=$OPTARG
             echo "Skip the second phase: ${skip_second_phase}"
@@ -100,35 +124,40 @@ fi
 # RK tolerance iteration #
 ##########################
 
-echo "Starting rk toleracne iteration..."
-
-# Remove the old directories
-for (( i=0; i < ${n_threads}; ++i ))
-do
-    rm -rf ${PWD}/thread_${i}
-done
-rm -rf ${PWD}/merged
-
-for (( i=0; i < ${n_threads}; ++i ))
-do
-    n_skips=`expr ${i} \* ${n_tracks_per_thread}`
+if [ "$skip_first_phase" = false ] ; then
     
-    command_rk_tolerance="${dir}/detray_integration_test_jacobian_validation \
-    --output-directory=thread_${i} \
-    --rk-tolerance-iterate-mode=true \
-    --n-tracks=${n_tracks_per_thread} \
-    --n-skips=${n_skips} \
-    --log10-min-rk-tolerance=${log10_min_rk_tol} \
-    --log10-max-rk-tolerance=${log10_max_rk_tol} \
-    --log10-helix-tolerance=${log10_helix_tol} \
-    --log10-on-surface-tolerance=${log10_on_surface_tol} \
-    --mc-seed=${mc_seed} \
-    --verbose-level=${verbose_level}"
-    ${command_rk_tolerance} &
-done
-wait
-
-echo "Finished rk toleracne iteration"
+    echo "Starting rk toleracne iteration..."
+    
+    # Remove the old directories
+    for (( i=0; i < ${n_threads}; ++i ))
+    do
+        rm -rf ${PWD}/thread_${i}
+    done
+    rm -rf ${PWD}/merged
+    
+    for (( i=0; i < ${n_threads}; ++i ))
+    do
+        n_skips=`expr ${i} \* ${n_tracks_per_thread}`
+        
+        command_rk_tolerance="${dir}/detray_integration_test_jacobian_validation \
+        --output-directory=thread_${i} \
+        --rk-tolerance-iterate-mode=true \
+        --n-tracks=${n_tracks_per_thread} \
+        --n-skips=${n_skips} \
+        --log10-min-rk-tolerance=${log10_min_rk_tol} \
+        --log10-max-rk-tolerance=${log10_max_rk_tol} \
+        --log10-helix-tolerance=${log10_helix_tol} \
+        --log10-on-surface-tolerance=${log10_on_surface_tol} \
+        --mc-seed=${mc_seed} \
+        --include-multiple-scattering=${include_multiple_scattering} \
+        --verbose-level=${verbose_level}"
+        ${command_rk_tolerance} &
+    done
+    wait
+    
+    echo "Finished rk toleracne iteration"
+    
+fi
 
 #####################################
 # Jacobi validation & Cov transport #
@@ -147,10 +176,11 @@ if [ "$skip_second_phase" = false ] ; then
         --rk-tolerance-iterate-mode=false \
         --n-tracks=${n_tracks_per_thread} \
         --n-skips=${n_skips} \
-        --log10-rk-tolerance=${log10_min_rk_tol} \
+        --log10-rk-tolerance=${log10_rk_tol} \
         --log10-helix-tolerance=${log10_helix_tol} \
         --log10-on-surface-tolerance=${log10_on_surface_tol} \
-        --mc-seed=${mc_seed}
+        --mc-seed=${mc_seed} \
+        --include-multiple-scattering=${include_multiple_scattering} \
         --verbose-level=${verbose_level}"
         ${command_jacobi_validation} &
     done
@@ -201,7 +231,7 @@ echo "Finished merging Csv files"
 
 cd ${output_dir}
 
-if [ "$skip_second_phase" = false ] ; then
+if [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = false ]; then
     
     # Run rk_tolerance_comparision.C
     root -q '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')'
@@ -211,8 +241,14 @@ if [ "$skip_second_phase" = false ] ; then
     
     # Run covariance_validation.C
     root -q -l ../../../tests/validation/root/covariance_validation.C+O
-else 
-
+    
+    elif [ "$skip_first_phase" = true ] && [ "$skip_second_phase" = false ]; then
+    
+    # Run covariance_validation.C
+    root ../../../tests/validation/root/covariance_validation.C+O
+    
+    elif [ "$skip_first_phase" = false ] && [ "$skip_second_phase" = true ]; then
+    
     # Run rk_tolerance_comparision.C
     root '../../../tests/validation/root/rk_tolerance_comparison.C+O('${log10_min_rk_tol}','${log10_max_rk_tol}')'
 fi
diff --git a/tests/unit_tests/cpu/material/energy_loss.cpp b/tests/unit_tests/cpu/material/energy_loss.cpp
index e6c1850de..9b644298e 100644
--- a/tests/unit_tests/cpu/material/energy_loss.cpp
+++ b/tests/unit_tests/cpu/material/energy_loss.cpp
@@ -10,9 +10,9 @@
 #include "detray/materials/material.hpp"
 #include "detray/materials/material_slab.hpp"
 #include "detray/materials/predefined_materials.hpp"
-#include "detray/simulation/landau_distribution.hpp"
-#include "detray/simulation/random_scatterer.hpp"
+#include "detray/propagator/actors/random_scatterer.hpp"
 #include "detray/test/types.hpp"
+#include "detray/utils/landau_distribution.hpp"
 #include "detray/utils/statistics.hpp"
 
 // GTest include(s).
diff --git a/tests/unit_tests/cpu/simulation/landau_sampling.cpp b/tests/unit_tests/cpu/simulation/landau_sampling.cpp
index 5734c5f8a..37d1da866 100644
--- a/tests/unit_tests/cpu/simulation/landau_sampling.cpp
+++ b/tests/unit_tests/cpu/simulation/landau_sampling.cpp
@@ -6,7 +6,7 @@
  */
 
 // Project include(s).
-#include "detray/simulation/landau_distribution.hpp"
+#include "detray/utils/landau_distribution.hpp"
 
 // GTest include(s)
 #include <gtest/gtest.h>
diff --git a/tests/unit_tests/cpu/simulation/scattering.cpp b/tests/unit_tests/cpu/simulation/scattering.cpp
index 72040b1e5..912aa7928 100644
--- a/tests/unit_tests/cpu/simulation/scattering.cpp
+++ b/tests/unit_tests/cpu/simulation/scattering.cpp
@@ -7,8 +7,8 @@
 
 // Project include(s).
 #include "detray/propagator/actors/pointwise_material_interactor.hpp"
-#include "detray/simulation/random_scatterer.hpp"
-#include "detray/simulation/scattering_helper.hpp"
+#include "detray/propagator/actors/random_scatterer.hpp"
+#include "detray/propagator/detail/scattering_helper.hpp"
 #include "detray/test/types.hpp"
 #include "detray/tracks/tracks.hpp"
 #include "detray/utils/statistics.hpp"
diff --git a/tests/unit_tests/cpu/utils/matrix_helper.cpp b/tests/unit_tests/cpu/utils/matrix_helper.cpp
index c749523ec..71862987c 100644
--- a/tests/unit_tests/cpu/utils/matrix_helper.cpp
+++ b/tests/unit_tests/cpu/utils/matrix_helper.cpp
@@ -165,4 +165,47 @@ GTEST_TEST(detray_utils, cholesky_decomposition) {
                             static_cast<float>(getter::element(B, i, j)));
         }
     }
-}
\ No newline at end of file
+}
+
+GTEST_TEST(detray_utils, curvilinear) {
+
+    const vector3 t1 = vector::normalize(vector3{1.f, 2.f, 3.f});
+    const vector3 z{0.f, 0.f, 1.f};
+
+    // Convert the (0,0,1) in local coordinate to global cooridnate
+    const matrix_type<3, 3> R1 =
+        matrix_helper<matrix_operator>().curvilinear_to_global(t1);
+    const vector3 t2 = R1 * z;
+
+    // The converted vector should be eqaul to the local-z axis in global
+    // coordinate
+    EXPECT_NEAR(static_cast<float>(t1[0]), static_cast<float>(t2[0]),
+                tolerance);
+    EXPECT_NEAR(static_cast<float>(t1[1]), static_cast<float>(t2[1]),
+                tolerance);
+    EXPECT_NEAR(static_cast<float>(t1[2]), static_cast<float>(t2[2]),
+                tolerance);
+
+    // Convert the local-z axis in global coordinate to local coordinate
+    const matrix_type<3, 3> R2 =
+        matrix_helper<matrix_operator>().global_to_curvilinear(t1);
+    const vector3 t3 = R2 * t1;
+
+    // The converted vector should be equal to unit z vector
+    EXPECT_NEAR(static_cast<float>(z[0]), static_cast<float>(t3[0]), tolerance);
+    EXPECT_NEAR(static_cast<float>(z[1]), static_cast<float>(t3[1]), tolerance);
+    EXPECT_NEAR(static_cast<float>(z[2]), static_cast<float>(t3[2]), tolerance);
+
+    // Round trip test
+    const matrix_type<3, 3> I = R1 * R2;
+
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 0, 0)), 1.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 0, 1)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 0, 2)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 1, 0)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 1, 1)), 1.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 1, 2)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 2, 0)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 2, 1)), 0.f, tolerance);
+    EXPECT_NEAR(static_cast<float>(getter::element(I, 2, 2)), 1.f, tolerance);
+}
diff --git a/tests/validation/root/covariance_validation.C b/tests/validation/root/covariance_validation.C
index 4cb858bda..6481d09b8 100644
--- a/tests/validation/root/covariance_validation.C
+++ b/tests/validation/root/covariance_validation.C
@@ -30,6 +30,9 @@
 #include <vector>
 
 namespace {
+double pull_min = -5.;
+double pull_max = 5.;
+double bin_width = 0.1;
 double x_pos = 0.16f;
 double title_x = x_pos;
 double title_y = 0.8f;
@@ -124,8 +127,8 @@ std::pair<double, double> fit_pval(TH1D* h_pval) {
 }
 
 void set_yaxis_title(TH1D* h, const double text_size) {
-    double bin_width = h->GetBinWidth(0u);
-    std::string str = std::to_string(bin_width);
+    double bw = h->GetBinWidth(0u);
+    std::string str = std::to_string(bw);
     str.erase(str.find_last_not_of('0') + 1, std::string::npos);
     str.erase(str.find_last_not_of('.') + 1, std::string::npos);
     std::string y_axis_title = "Counts / (" + str + ")";
@@ -322,9 +325,7 @@ void read_tree(TTree* t, const std::string& tag, const std::string& title) {
     const std::array<float, 2> cdim1{700, 500};
     const std::array<float, 2> cdim2{700, 500};
 
-    float pull_min = -5.f;
-    float pull_max = 5.f;
-    int n_bins = 100;
+    int n_bins = (pull_max - pull_min) / bin_width;
 
     double pull_l0;
     double pull_l1;
diff --git a/tests/validation/root/rk_tolerance_comparison.C b/tests/validation/root/rk_tolerance_comparison.C
index b6709de03..d5744f4d4 100644
--- a/tests/validation/root/rk_tolerance_comparison.C
+++ b/tests/validation/root/rk_tolerance_comparison.C
@@ -217,7 +217,7 @@ void draw_mean_step_size(const std::string header_title,
     gr->GetXaxis()->SetTitle("log_{10}(#tau [mm])");
     gr->GetYaxis()->SetTitle("log_{10}(Mean of average step sizes [mm])");
     gr->GetXaxis()->SetLimits(x_vec.front() - 0.5, x_vec.back() + 0.5);
-    gr->GetYaxis()->SetRangeUser(means.front() - 0.2, means.back() + 0.95);
+    gr->GetYaxis()->SetRangeUser(means.front() - 0.2, means.back() + 1.00);
     gr->GetXaxis()->SetLabelSize(label_font_size);
     gr->GetYaxis()->SetLabelSize(label_font_size);
     gr->GetXaxis()->SetTitleSize(title_font_size);