Adapt actor chain for empty state and to produce actor state tuple if states are default constructible (#879)

Shift the update of the step size in the stepper, so that we go from a "previous step size" to a "next step size".
- Remove the volume material pointer from the state and instead get it from the navigator hand it through to the functions that need it
- Remove the "path from surface" and reimplement it in the actors that use it
- Remove the next_surface_aborter, which is superseded by the ckf_aborter in traccc (need to move the test to traccc)
- Remove the previous surface index and instead use the surface link on the bound track parameters before they are updated
- Remove the stepping direction and calculate it on the fly

Also streamlines the parameter resetter, as the jacobians are already default initialized correctly and a separate kernel is not needed, since the relevant functionality already exists on the tracking_surface

Author: niermann999

core/include/detray/navigation/policies.hpp

@@ -73,8 +73,7 @@ struct stepper_default_policy : actor {
         auto &navigation = propagation._navigation;
 
         // 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) {
@@ -111,20 +110,19 @@ struct stepper_rk_policy : actor {
         const auto &stepping = propagation._stepping;
         auto &navigation = propagation._navigation;
 
-        // Policy is called after stepsize update -> use prev. step size
-        const scalar rel_correction{(stepping.prev_step_size() - navigation()) /
+        // How strongly did the RKN algorithm reduce the step size?
+        const scalar rel_correction{(stepping.step_size() - navigation()) /
                                     navigation()};
 
         // Large correction to the stepsize - re-initialize the volume
         if (rel_correction > pol_state.m_threshold_no_trust) {
-            // Re-evaluate only next candidate
             navigation.set_no_trust();
         }
-        // Medium correction - re-evaluate the current candidates
+        // Medium correction - re-evaluate all current candidates
         else if (rel_correction > pol_state.m_threshold_fair_trust) {
-            // Re-evaluate all candidates
             navigation.set_fair_trust();
         } else {
+            // Small correction - re-evaluate only the next candidate
             navigation.set_high_trust();
         }
     }

core/include/detray/propagator/actor_chain.hpp

@@ -60,17 +60,19 @@ class actor_chain {
         // Only possible if each state is default initializable
         if constexpr ((std::default_initializable<typename actors_t::state> &&
                        ...)) {
-            tuple_t<typename actors_t::state...> states{};
-
-            return std::make_tuple(
-                states,
-                make_ref_tuple(
-                    states, std::make_index_sequence<sizeof...(actors_t)>{}));
+            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.
     ///

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,40 +19,6 @@ namespace detray {
 template <typename algebra_t>
 struct parameter_resetter : actor {
 
-    using scalar_type = dscalar<algebra_t>;
-
-    /// 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()(propagator_state_t& propagation) const {
 
@@ -65,11 +31,13 @@ struct parameter_resetter : actor {
             return;
         }
 
-        // Surface
+        // Update free params after bound params were changed by actors
         const auto sf = navigation.get_surface();
+        stepping() = sf.bound_to_free_vector(propagation._context,
+                                             stepping.bound_params());
 
-        sf.template visit_mask<kernel>(sf.transform(propagation._context),
-                                       sf.index(), stepping);
+        // Reset jacobian transport to identity matrix
+        stepping.reset_transport_jacobian();
     }
 };
 

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

@@ -21,29 +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 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>;
@@ -54,32 +52,23 @@ 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;
         }
     };
 
@@ -94,46 +83,52 @@ struct parameter_transporter : actor {
             return;
         }
 
-        using detector_type = typename propagator_state_t::detector_type;
+        // Geometry context for this track
+        const auto& gctx = propagation._context;
 
         // 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(propagation._context,
-                                               stepping.bound_params());
+            const bound_to_free_matrix_t bound_to_free_jacobian =
+                prev_sf.bound_to_free_jacobian(gctx, 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(propagation._context), bound_to_free_jacobian,
-                    propagation));
+                    sf.transform(gctx), 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(
-            sf.free_to_bound_vector(propagation._context, stepping()));
+        bound_params.set_parameter_vector(
+            sf.free_to_bound_vector(gctx, 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