diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp index 904d664eb..bef0c0b38 100644 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ b/core/include/detray/propagator/actors/parameter_resetter.hpp @@ -39,9 +39,6 @@ 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; - // Reset the free vector stepping() = detail::bound_to_free_vector(trf3, mask, stepping._bound_params); @@ -49,10 +46,6 @@ struct parameter_resetter : actor { // 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(sf.transform(ctx), stepping); + + stepping._prev_geo_id = sf.barcode(); } }; diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index 850493848..c3e6cd80c 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/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; + + template + DETRAY_HOST_DEVICE inline bound_to_free_matrix 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; + + using jacobian_engine_t = detail::jacobian_engine; + + 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,7 +64,9 @@ 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& bound_to_free_jacobian, + propagator_state_t& propagation) { using frame_t = typename mask_group_t::value_type::shape:: template local_frame_type; @@ -51,8 +74,6 @@ struct parameter_transporter : actor { using jacobian_engine_t = detail::jacobian_engine; using bound_matrix_t = bound_matrix; - using bound_to_free_matrix_t = - typename jacobian_engine_t::bound_to_free_matrix_type; using free_matrix_t = free_matrix; using free_to_bound_matrix_t = @@ -85,10 +106,6 @@ struct parameter_transporter : actor { .template identity() + 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 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; + 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(sf.transform(ctx), propagation); + bound_to_free_matrix bound_to_free_jacobian = + matrix_operator().template zero(); + + if (!stepping._prev_geo_id.is_invalid()) { + + // Previous surface + tracking_surface prev_sf{navigation.detector(), + stepping._prev_geo_id}; + + bound_to_free_jacobian = + prev_sf.template visit_mask( + prev_sf.transform(ctx), propagation); + } + + sf.template visit_mask(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 diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index 788980bf5..300b0dc91 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/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(); - /// bound-to-free jacobian from departure surface - bound_to_free_matrix_type _jac_to_global = - matrix_operator().template zero(); - /// bound covariance bound_track_parameters_type _bound_params; @@ -164,6 +158,9 @@ class base_stepper { /// The default particle hypothesis is muon pdg_particle _ptc = muon(); + /// Previous barcode to calculate the bound_to_free_jacobian + geometry::barcode _prev_geo_id; + /// Set new step constraint template DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) {