Remove jac_to_global from the base stepper

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

@@ -39,20 +39,13 @@ struct parameter_resetter : actor {
             // Note: How is it possible with "range"???
             const auto& mask = mask_group[index];
 
-            using frame_t = decltype(mask.local_frame());
-            using jacobian_engine = detail::jacobian_engine<frame_t>;
-
             // Reset the free vector
             stepping() = detail::bound_to_free_vector(trf3, mask,
                                                       stepping._bound_params);
 
             // Reset the path length
             stepping._s = 0;
 
-            // Reset jacobian coordinate transformation at the current surface
-            stepping._jac_to_global = jacobian_engine::bound_to_free_jacobian(
-                trf3, mask, stepping._bound_params);
-
             // Reset jacobian transport to identity matrix
             matrix_operator().set_identity(stepping._jac_transport);
         }
@@ -79,6 +72,8 @@ struct parameter_resetter : actor {
         const auto sf = navigation.get_surface();
 
         sf.template visit_mask<kernel>(sf.transform(ctx), stepping);
+
+        stepping._prev_geo_id = sf.barcode();
     }
 };
 

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

@@ -21,7 +21,28 @@ struct parameter_transporter : actor {
 
     struct state {};
 
-    /// Mask store visitor
+    /// bound to free jacobian from previous surface
+    struct bound_to_free_jacobian_kernel {
+
+        // Transformation matching this struct
+        using transform3_type = dtransform3D<algebra_t>;
+
+        template <typename mask_group_t, typename index_t,
+                  typename propagator_state_t>
+        DETRAY_HOST_DEVICE inline bound_to_free_matrix<algebra_t> operator()(
+            const mask_group_t& mask_group, const index_t& index,
+            const transform3_type& trf3, propagator_state_t& propagation) {
+
+            using frame_t = typename mask_group_t::value_type::shape::
+                template local_frame_type<algebra_t>;
+
+            using jacobian_engine_t = detail::jacobian_engine<frame_t>;
+
+            return jacobian_engine_t::bound_to_free_jacobian(
+                trf3, mask_group[index], propagation._stepping._bound_params);
+        }
+    };
+
     struct kernel {
 
         /// @name Type definitions for the struct
@@ -43,16 +64,16 @@ 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,
+            const bound_to_free_matrix<algebra_t>& bound_to_free_jacobian,
+            propagator_state_t& propagation) {
 
             using frame_t = typename mask_group_t::value_type::shape::
                 template local_frame_type<algebra_t>;
 
             using jacobian_engine_t = detail::jacobian_engine<frame_t>;
 
             using bound_matrix_t = bound_matrix<algebra_t>;
-            using bound_to_free_matrix_t =
-                typename jacobian_engine_t::bound_to_free_matrix_type;
 
             using free_matrix_t = free_matrix<algebra_t>;
             using free_to_bound_matrix_t =
@@ -85,10 +106,6 @@ struct parameter_transporter : actor {
                     .template identity<e_free_size, e_free_size>() +
                 path_correction;
 
-            // Bound to free jacobian at the departure surface
-            const bound_to_free_matrix_t& bound_to_free_jacobian =
-                stepping._jac_to_global;
-
             stepping._full_jacobian = free_to_bound_jacobian * correction_term *
                                       free_transport_jacobian *
                                       bound_to_free_jacobian;
@@ -105,6 +122,7 @@ struct parameter_transporter : actor {
     template <typename propagator_state_t>
     DETRAY_HOST_DEVICE void operator()(state& /*actor_state*/,
                                        propagator_state_t& propagation) const {
+        auto& stepping = propagation._stepping;
         const auto& navigation = propagation._navigation;
 
         // Do covariance transport when the track is on surface
@@ -113,17 +131,33 @@ struct parameter_transporter : actor {
             return;
         }
 
-        using geo_cxt_t =
-            typename propagator_state_t::detector_type::geometry_context;
+        using matrix_operator = dmatrix_operator<algebra_t>;
+        using detector_type = typename propagator_state_t::detector_type;
+        using geo_cxt_t = typename detector_type::geometry_context;
         const geo_cxt_t ctx{};
 
-        // Surface
+        // Current Surface
         const auto sf = navigation.get_surface();
 
-        sf.template visit_mask<kernel>(sf.transform(ctx), propagation);
+        bound_to_free_matrix<algebra_t> bound_to_free_jacobian =
+            matrix_operator().template zero<e_free_size, e_bound_size>();
+
+        if (!stepping._prev_geo_id.is_invalid()) {
+
+            // Previous surface
+            tracking_surface<detector_type> prev_sf{navigation.detector(),
+                                                    stepping._prev_geo_id};
+
+            bound_to_free_jacobian =
+                prev_sf.template visit_mask<bound_to_free_jacobian_kernel>(
+                    prev_sf.transform(ctx), propagation);
+        }
+
+        sf.template visit_mask<kernel>(sf.transform(ctx),
+                                       bound_to_free_jacobian, propagation);
 
         // Set surface link
-        propagation._stepping._bound_params.set_surface_link(sf.barcode());
+        stepping._bound_params.set_surface_link(sf.barcode());
     }
 };  // namespace detray
 

core/include/detray/propagator/base_stepper.hpp

@@ -11,6 +11,7 @@
 #include "detray/definitions/detail/qualifiers.hpp"
 #include "detray/definitions/pdg_particle.hpp"
 #include "detray/definitions/units.hpp"
+#include "detray/geometry/barcode.hpp"
 #include "detray/geometry/tracking_surface.hpp"
 #include "detray/propagator/actors/parameter_resetter.hpp"
 #include "detray/propagator/constrained_step.hpp"
@@ -79,9 +80,6 @@ class base_stepper {
             // A dummy barcode - should not be used
             _bound_params.set_surface_link(geometry::barcode{});
 
-            // Set the bound to free jacobian
-            _jac_to_global = cf.bound_to_free_jacobian();
-
             // Reset the path length
             _s = 0.f;
 
@@ -116,10 +114,6 @@ class base_stepper {
         free_matrix_type _jac_transport =
             matrix_operator().template identity<e_free_size, e_free_size>();
 
-        /// bound-to-free jacobian from departure surface
-        bound_to_free_matrix_type _jac_to_global =
-            matrix_operator().template zero<e_free_size, e_bound_size>();
-
         /// bound covariance
         bound_track_parameters_type _bound_params;
 
@@ -164,6 +158,9 @@ class base_stepper {
         /// The default particle hypothesis is muon
         pdg_particle<scalar_type> _ptc = muon<scalar_type>();
 
+        /// Previous barcode to calculate the bound_to_free_jacobian
+        geometry::barcode _prev_geo_id;
+
         /// Set new step constraint
         template <step::constraint type = step::constraint::e_actor>
         DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) {