Adapt actor chain for empty state and to produce actor state tuple if states are default constructible

Author: niermann999

core/include/detray/navigation/policies.hpp

@@ -74,7 +74,7 @@ struct stepper_default_policy : actor {
 
         // Not a severe change to track state expected
         // Policy is called after stepsize update -> use prev. step size
-        if (math::fabs(stepping.prev_step_size()) <
+        if (math::fabs(stepping.step_size()) <
             math::fabs(
                 stepping.constraints().template size<>(stepping.direction())) -
                 pol_state.tol) {
@@ -112,7 +112,7 @@ struct stepper_rk_policy : actor {
         auto &navigation = propagation._navigation;
 
         // Policy is called after stepsize update -> use prev. step size
-        const scalar rel_correction{(stepping.prev_step_size() - navigation()) /
+        const scalar rel_correction{(stepping.step_size() - navigation()) /
                                     navigation()};
 
         // Large correction to the stepsize - re-initialize the volume

core/include/detray/propagator/actor_chain.hpp

@@ -1,19 +1,24 @@
 /** Detray library, part of the ACTS project (R&D line)
  *
- * (c) 2022 CERN for the benefit of the ACTS project
+ * (c) 2022-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
 
 #pragma once
 
-#include <type_traits>
-#include <utility>
-
+// Project include(s)
 #include "detray/definitions/detail/containers.hpp"
 #include "detray/definitions/detail/qualifiers.hpp"
+#include "detray/propagator/base_actor.hpp"
 #include "detray/utils/tuple_helpers.hpp"
 
+// System include(s)
+#include <concepts>
+#include <optional>
+#include <type_traits>
+#include <utility>
+
 namespace detray {
 
 /// The interface to the actors and aborters in the propagation.
@@ -44,7 +49,29 @@ class actor_chain {
     }
 
     /// @returns the actor list
-    DETRAY_HOST_DEVICE const actor_list_type &actors() const { return _actors; }
+    DETRAY_HOST_DEVICE const actor_list_type &actors() const {
+        return m_actors;
+    }
+
+    /// @returns a tuple of default constructible actor states and a
+    /// corresponding tuple of references
+    DETRAY_HOST_DEVICE
+    static constexpr auto make_actor_states() {
+        // Only possible if each state is default initializable
+        if constexpr ((std::default_initializable<typename actors_t::state> &&
+                       ...)) {
+            return tuple_t<typename actors_t::state...>{};
+        } else {
+            return std::nullopt;
+        }
+    }
+
+    /// @returns a tuple of reference for every state in the tuple @param t
+    DETRAY_HOST_DEVICE static constexpr state make_ref_tuple(
+        tuple_t<typename actors_t::state...> &t) {
+        return make_ref_tuple(t,
+                              std::make_index_sequence<sizeof...(actors_t)>{});
+    }
 
     private:
     /// Call the actors. Either single actor or composition.
@@ -58,7 +85,13 @@ class actor_chain {
                                        actor_states_t &states,
                                        propagator_state_t &p_state) const {
         if constexpr (!typename actor_t::is_comp_actor()) {
-            actr(detail::get<typename actor_t::state &>(states), p_state);
+            // No actor state defined (empty)
+            if constexpr (std::same_as<typename actor_t::state,
+                                       detray::actor::state>) {
+                actr(p_state);
+            } else {
+                actr(detail::get<typename actor_t::state &>(states), p_state);
+            }
         } else {
             actr(states, p_state);
         }
@@ -73,11 +106,19 @@ class actor_chain {
     DETRAY_HOST_DEVICE inline void run(
         actor_states_t &states, propagator_state_t &p_state,
         std::index_sequence<indices...> /*ids*/) const {
-        (run(detail::get<indices>(_actors), states, p_state), ...);
+        (run(detail::get<indices>(m_actors), states, p_state), ...);
+    }
+
+    /// @returns a tuple of reference for every state in the tuple @param t
+    template <std::size_t... indices>
+    DETRAY_HOST_DEVICE static constexpr state make_ref_tuple(
+        tuple_t<typename actors_t::state...> &t,
+        std::index_sequence<indices...> /*ids*/) {
+        return detray::tie(detail::get<indices>(t)...);
     }
 
     /// Tuple of actors
-    actor_list_type _actors = {};
+    actor_list_type m_actors = {};
 };
 
 /// Empty actor chain (placeholder)

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

@@ -95,28 +95,4 @@ struct target_aborter : actor {
     }
 };
 
-/// Aborter triggered when the next surface is reached
-struct next_surface_aborter : actor {
-    struct state {
-        // minimal step length to prevent from staying on the same surface
-        scalar min_step_length = 0.f;
-        bool success = false;
-    };
-
-    template <typename propagator_state_t>
-    DETRAY_HOST_DEVICE void operator()(state &abrt_state,
-                                       propagator_state_t &prop_state) const {
-
-        auto &navigation = prop_state._navigation;
-        auto &stepping = prop_state._stepping;
-
-        // Abort at the next sensitive surface
-        if (navigation.is_on_sensitive() &&
-            stepping.path_from_surface() > abrt_state.min_step_length) {
-            prop_state._heartbeat &= navigation.abort();
-            abrt_state.success = true;
-        }
-    }
-};
-
 }  // namespace detray

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

@@ -1,6 +1,6 @@
 /** Detray library, part of the ACTS project (R&D line)
  *
- * (c) 2022-2023 CERN for the benefit of the ACTS project
+ * (c) 2022-2024 CERN for the benefit of the ACTS project
  *
  * Mozilla Public License Version 2.0
  */
@@ -19,45 +19,8 @@ namespace detray {
 template <typename algebra_t>
 struct parameter_resetter : actor {
 
-    using scalar_type = dscalar<algebra_t>;
-
-    struct state {};
-
-    /// Mask store visitor
-    struct kernel {
-
-        // Matrix actor
-        using transform3_type = dtransform3D<algebra_t>;
-        using matrix_operator = dmatrix_operator<algebra_t>;
-
-        template <typename mask_group_t, typename index_t,
-                  typename stepper_state_t>
-        DETRAY_HOST_DEVICE inline void operator()(
-            const mask_group_t& mask_group, const index_t& index,
-            const transform3_type& trf3, const dindex sf_idx,
-            stepper_state_t& stepping) const {
-
-            // Note: How is it possible with "range"???
-            const auto& mask = mask_group[index];
-
-            // Reset the free vector
-            stepping() = detail::bound_to_free_vector(trf3, mask,
-                                                      stepping.bound_params());
-
-            // Reset the path length
-            stepping.reset_path_from_surface();
-
-            // Reset jacobian transport to identity matrix
-            stepping.reset_transport_jacobian();
-
-            // Reset the surface index
-            stepping.set_prev_sf_index(sf_idx);
-        }
-    };
-
     template <typename propagator_state_t>
-    DETRAY_HOST_DEVICE void operator()(state& /*resetter_state*/,
-                                       propagator_state_t& propagation) const {
+    DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const {
 
         const auto& navigation = propagation._navigation;
         auto& stepping = propagation._stepping;
@@ -68,14 +31,14 @@ struct parameter_resetter : actor {
             return;
         }
 
-        using geo_cxt_t =
-            typename propagator_state_t::detector_type::geometry_context;
-        const geo_cxt_t ctx{};
+        typename propagator_state_t::detector_type::geometry_context ctx{};
 
-        // Surface
+        // Update free params after bound params were changed by actors
         const auto sf = navigation.get_surface();
+        stepping() = sf.bound_to_free_vector(ctx, stepping.bound_params());
 
-        sf.template visit_mask<kernel>(sf.transform(ctx), sf.index(), stepping);
+        // Reset jacobian transport to identity matrix
+        stepping.reset_transport_jacobian();
     }
 };
 

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

@@ -21,31 +21,27 @@ struct parameter_transporter : actor {
 
     /// @name Type definitions for the struct
     /// @{
-
+    using scalar_type = dscalar<algebra_t>;
     // Transformation matching this struct
     using transform3_type = dtransform3D<algebra_t>;
-    // scalar_type
-    using scalar_type = dscalar<algebra_t>;
     // Matrix actor
     using matrix_operator = dmatrix_operator<algebra_t>;
-    // 2D matrix type
-    template <std::size_t ROWS, std::size_t COLS>
-    using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
     // bound matrix type
     using bound_matrix_t = bound_matrix<algebra_t>;
+    // Matrix type for bound to free jacobian
+    using bound_to_free_matrix_t = bound_to_free_matrix<algebra_t>;
     /// @}
 
-    struct state {};
-
     struct get_full_jacobian_kernel {
 
         template <typename mask_group_t, typename index_t,
-                  typename propagator_state_t>
+                  typename stepper_state_t>
         DETRAY_HOST_DEVICE inline bound_matrix_t operator()(
             const mask_group_t& /*mask_group*/, const index_t& /*index*/,
             const transform3_type& trf3,
-            const bound_to_free_matrix<algebra_t>& bound_to_free_jacobian,
-            const propagator_state_t& propagation) const {
+            const bound_to_free_matrix_t& bound_to_free_jacobian,
+            const material<scalar_type>* vol_mat_ptr,
+            const stepper_state_t& stepping) const {
 
             using frame_t = typename mask_group_t::value_type::shape::
                 template local_frame_type<algebra_t>;
@@ -56,38 +52,28 @@ struct parameter_transporter : actor {
             using free_to_bound_matrix_t =
                 typename jacobian_engine_t::free_to_bound_matrix_type;
 
-            // Stepper and Navigator states
-            auto& stepping = propagation._stepping;
-
-            // Free vector
-            const auto& free_params = stepping();
-
             // Free to bound jacobian at the destination surface
             const free_to_bound_matrix_t free_to_bound_jacobian =
-                jacobian_engine_t::free_to_bound_jacobian(trf3, free_params);
-
-            // Transport jacobian in free coordinate
-            const free_matrix_t& free_transport_jacobian =
-                stepping.transport_jacobian();
+                jacobian_engine_t::free_to_bound_jacobian(trf3, stepping());
 
             // Path correction factor
-            free_matrix_t path_correction = jacobian_engine_t::path_correction(
-                stepping().pos(), stepping().dir(), stepping.dtds(),
-                stepping.dqopds(), trf3);
+            const free_matrix_t path_correction =
+                jacobian_engine_t::path_correction(
+                    stepping().pos(), stepping().dir(), stepping.dtds(),
+                    stepping.dqopds(vol_mat_ptr), trf3);
 
             const free_matrix_t correction_term =
                 matrix_operator()
                     .template identity<e_free_size, e_free_size>() +
                 path_correction;
 
             return free_to_bound_jacobian * correction_term *
-                   free_transport_jacobian * bound_to_free_jacobian;
+                   stepping.transport_jacobian() * bound_to_free_jacobian;
         }
     };
 
     template <typename propagator_state_t>
-    DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
-                                       propagator_state_t& propagation) const {
+    DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const {
         auto& stepping = propagation._stepping;
         const auto& navigation = propagation._navigation;
 
@@ -97,46 +83,51 @@ struct parameter_transporter : actor {
             return;
         }
 
-        using detector_type = typename propagator_state_t::detector_type;
-        using geo_cxt_t = typename detector_type::geometry_context;
-        const geo_cxt_t ctx{};
+        typename propagator_state_t::detector_type::geometry_context ctx{};
 
         // Current Surface
         const auto sf = navigation.get_surface();
 
+        // Bound track params of departure surface
+        auto& bound_params = stepping.bound_params();
+
         // Covariance is transported only when the previous surface is an
         // actual tracking surface. (i.e. This disables the covariance transport
         // from curvilinear frame)
-        if (!detail::is_invalid_value(stepping.prev_sf_index())) {
+        if (!bound_params.surface_link().is_invalid()) {
 
             // Previous surface
-            tracking_surface<detector_type> prev_sf{navigation.detector(),
-                                                    stepping.prev_sf_index()};
+            tracking_surface prev_sf{navigation.detector(),
+                                     bound_params.surface_link()};
 
-            const bound_to_free_matrix<algebra_t> bound_to_free_jacobian =
-                prev_sf.bound_to_free_jacobian(ctx, stepping.bound_params());
+            const bound_to_free_matrix_t bound_to_free_jacobian =
+                prev_sf.bound_to_free_jacobian(ctx, bound_params);
 
+            auto vol = navigation.get_volume();
+            const auto vol_mat_ptr =
+                vol.has_material() ? vol.material_parameters(stepping().pos())
+                                   : nullptr;
             stepping.set_full_jacobian(
                 sf.template visit_mask<get_full_jacobian_kernel>(
-                    sf.transform(ctx), bound_to_free_jacobian, propagation));
+                    sf.transform(ctx), bound_to_free_jacobian, vol_mat_ptr,
+                    propagation._stepping));
 
             // Calculate surface-to-surface covariance transport
             const bound_matrix_t new_cov =
-                stepping.full_jacobian() *
-                stepping.bound_params().covariance() *
+                stepping.full_jacobian() * bound_params.covariance() *
                 matrix_operator().transpose(stepping.full_jacobian());
+
             stepping.bound_params().set_covariance(new_cov);
         }
 
         // Convert free to bound vector
-        stepping.bound_params().set_parameter_vector(
+        bound_params.set_parameter_vector(
             sf.free_to_bound_vector(ctx, stepping()));
 
         // Set surface link
-        stepping.bound_params().set_surface_link(sf.barcode());
-
-        return;
+        bound_params.set_surface_link(sf.barcode());
     }
+
 };  // namespace detray
 
 }  // namespace detray

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

@@ -35,7 +35,6 @@ struct pointwise_material_interactor : actor {
 
     struct state {
 
-        /// @TODO: Consider using the particle information in stepping::config
         /// Evaluated energy loss
         scalar_type e_loss{0.f};
         /// Evaluated projected scattering angle