Multiple scattering in the volume material

.gitignore

@@ -15,4 +15,4 @@
 **/*.d
 **/*.so
 **/*C_ACLiC_dict*
-**/therad_*
+**/thread_*

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;
         }

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

@@ -49,14 +49,20 @@ 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(
                 trf3, mask, stepping._bound_params.vector());
 
             // Reset jacobian transport to identity matrix
             matrix_operator().set_identity(stepping._jac_transport);
+
+            // Reset the joint covariance
+            matrix_operator().set_zero(stepping._joint_cov);
         }
     };
 

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());

utils/include/detray/simulation/random_scatterer.hpp → 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"
 

core/include/detray/propagator/base_stepper.hpp

@@ -13,13 +13,19 @@
 #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"
 
 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) {

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

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}; }
 

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