From a6c92913c70a3ec17099bfe5b45ad864dcc8428f Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Tue, 5 Nov 2024 14:48:12 +0100 Subject: [PATCH 1/2] Adapt actor chain for empty state and to produce actor state tuple if states are default constructible --- core/include/detray/navigation/policies.hpp | 4 +- .../include/detray/propagator/actor_chain.hpp | 14 +- .../detray/propagator/actors/aborters.hpp | 24 --- .../propagator/actors/parameter_resetter.hpp | 46 +----- .../actors/parameter_transporter.hpp | 72 ++++----- .../actors/pointwise_material_interactor.hpp | 1 - .../detray/propagator/base_stepper.hpp | 144 +++++------------- .../detray/propagator/line_stepper.hpp | 13 +- core/include/detray/propagator/propagator.hpp | 21 +-- core/include/detray/propagator/rk_stepper.hpp | 38 +++-- core/include/detray/propagator/rk_stepper.ipp | 125 ++++++++------- .../cpu/material/material_interaction.cpp | 50 ------ .../cpu/propagator/jacobian_validation.cpp | 6 +- .../cpu/propagator/propagator.cpp | 74 +++++---- .../unit_tests/cpu/propagator/rk_stepper.cpp | 10 +- 15 files changed, 243 insertions(+), 399 deletions(-) diff --git a/core/include/detray/navigation/policies.hpp b/core/include/detray/navigation/policies.hpp index 4e1abc430..dac2c5d32 100644 --- a/core/include/detray/navigation/policies.hpp +++ b/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 diff --git a/core/include/detray/propagator/actor_chain.hpp b/core/include/detray/propagator/actor_chain.hpp index b80ce5a4d..572880fc9 100644 --- a/core/include/detray/propagator/actor_chain.hpp +++ b/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 && ...)) { - tuple_t states{}; - - return std::make_tuple( - states, - make_ref_tuple( - states, std::make_index_sequence{})); + return tuple_t{}; } 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 &t) { + return make_ref_tuple(t, + std::make_index_sequence{}); + } + private: /// Call the actors. Either single actor or composition. /// diff --git a/core/include/detray/propagator/actors/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp index 81fa2f18e..a2de17dba 100644 --- a/core/include/detray/propagator/actors/aborters.hpp +++ b/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 - 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 diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp index 193071344..79bd20177 100644 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ b/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 struct parameter_resetter : actor { - using scalar_type = dscalar; - - /// Mask store visitor - struct kernel { - - // Matrix actor - using transform3_type = dtransform3D; - using matrix_operator = dmatrix_operator; - - template - 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 DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const { @@ -65,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(sf.transform(ctx), sf.index(), stepping); + // Reset jacobian transport to identity matrix + stepping.reset_transport_jacobian(); } }; diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index d9ae6175e..78f9c07a6 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/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; // Transformation matching this struct using transform3_type = dtransform3D; - // scalar_type - using scalar_type = dscalar; // Matrix actor using matrix_operator = dmatrix_operator; - // 2D matrix type - template - using matrix_type = dmatrix; // bound matrix type using bound_matrix_t = bound_matrix; + // Matrix type for bound to free jacobian + using bound_to_free_matrix_t = bound_to_free_matrix; /// @} struct get_full_jacobian_kernel { template + 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& bound_to_free_jacobian, - const propagator_state_t& propagation) const { + const bound_to_free_matrix_t& bound_to_free_jacobian, + const material* vol_mat_ptr, + const stepper_state_t& stepping) const { using frame_t = typename mask_group_t::value_type::shape:: template local_frame_type; @@ -54,24 +52,15 @@ 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() @@ -79,7 +68,7 @@ struct parameter_transporter : actor { 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,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 prev_sf{navigation.detector(), - stepping.prev_sf_index()}; + tracking_surface prev_sf{navigation.detector(), + bound_params.surface_link()}; - const bound_to_free_matrix 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( - 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 diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index 58140c589..dae616d22 100644 --- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp +++ b/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 diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index 8ea616767..bb39149b3 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -14,7 +14,6 @@ #include "detray/geometry/barcode.hpp" #include "detray/geometry/tracking_surface.hpp" #include "detray/materials/material.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" #include "detray/propagator/constrained_step.hpp" #include "detray/propagator/stepping_config.hpp" #include "detray/tracks/tracks.hpp" @@ -75,9 +74,6 @@ class base_stepper { // An invalid barcode - should not be used m_bound_params.set_surface_link(geometry::barcode{}); - - // Reset jacobian transport to identity matrix - matrix_operator().set_identity(m_jac_transport); } /// Sets track parameters from bound track parameter. @@ -87,13 +83,12 @@ class base_stepper { const detector_t &det) : m_bound_params(bound_params) { - // Surface + // Departure surface const auto sf = tracking_surface{det, bound_params.surface_link()}; + // Set free track parameters for stepping/navigation const typename detector_t::geometry_context ctx{}; - sf.template visit_mask< - typename parameter_resetter::kernel>( - sf.transform(ctx), sf.index(), *this); + m_track = sf.bound_to_free_vector(ctx, bound_params); } /// @returns free track parameters - non-const access @@ -108,7 +103,7 @@ class base_stepper { DETRAY_HOST_DEVICE bound_track_parameters_type &bound_params() { return m_bound_params; } - /// @returns bound track parameters. + /// @returns bound track parameters - non-const access DETRAY_HOST_DEVICE const bound_track_parameters_type &bound_params() const { return m_bound_params; @@ -116,39 +111,30 @@ class base_stepper { /// Get stepping direction DETRAY_HOST_DEVICE - inline step::direction direction() const { return m_direction; } - - /// Set stepping direction - DETRAY_HOST_DEVICE inline void set_direction(step::direction dir) { - m_direction = dir; + inline step::direction direction() const { + return m_step_size >= 0.f ? step::direction::e_forward + : step::direction::e_backward; } + /// Updates the total number of step trials + DETRAY_HOST_DEVICE inline void count_trials() { ++m_n_total_trials; } + /// @returns the total number of step trials. For steppers that don't /// use adaptive step size scaling, this is the number of steps DETRAY_HOST_DEVICE inline std::size_t n_total_trials() const { return m_n_total_trials; } - /// Set next step size + /// Set the current step size DETRAY_HOST_DEVICE inline void set_step_size(const scalar_type step) { m_step_size = step; } - /// @returns the current step size of this state. + /// @returns the step size of the current step. DETRAY_HOST_DEVICE inline scalar_type step_size() const { return m_step_size; } - /// Set previous step size - DETRAY_HOST_DEVICE - inline void set_prev_step_size(const scalar_type step) { - m_prev_step_size = step; - } - - /// @returns the previous step size of this state. - DETRAY_HOST_DEVICE - inline scalar_type prev_step_size() const { return m_prev_step_size; } - /// @returns this states path length. DETRAY_HOST_DEVICE inline scalar_type path_length() const { return m_path_length; } @@ -157,25 +143,12 @@ class base_stepper { DETRAY_HOST_DEVICE inline scalar_type abs_path_length() const { return m_abs_path_length; } - /// @returns path length since last surface was encountered - DETRAY_HOST_DEVICE - inline scalar_type path_from_surface() const { return m_path_from_sf; } - - /// Reset the path length since last surface - DETRAY_HOST_DEVICE inline void reset_path_from_surface() { - m_path_from_sf = 0.f; - } - /// Add a new segment to all path lengths (forward or backward) DETRAY_HOST_DEVICE inline void update_path_lengths(scalar_type seg) { m_path_length += seg; - m_path_from_sf += seg; m_abs_path_length += math::fabs(seg); } - /// Updates the total number of step trials - DETRAY_HOST_DEVICE inline void count_trials() { ++m_n_total_trials; } - /// Set new step constraint template DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) { @@ -192,34 +165,7 @@ class base_stepper { m_constraint.template release(); } - /// @returns the index of the previous surface for cov. transport - DETRAY_HOST_DEVICE dindex prev_sf_index() const { return m_prev_sf_id; } - - /// Set the index of the previous surface for cov. transport - DETRAY_HOST_DEVICE - void set_prev_sf_index(dindex prev_idx) { - assert(!detail::is_invalid_value(prev_idx)); - m_prev_sf_id = prev_idx; - } - - /// @returns true if the current volume contains volume material - DETRAY_HOST_DEVICE - bool volume_has_material() const { return (m_vol_mat != nullptr); } - - /// Access the current volume material - DETRAY_HOST_DEVICE - const auto &volume_material() const { - assert(m_vol_mat != nullptr); - return *m_vol_mat; - } - - /// Set new volume material access - DETRAY_HOST_DEVICE - inline void set_volume_material(const material *mat_ptr) { - m_vol_mat = mat_ptr; - } - - /// Access the current particle hypothesis + /// @returns the current particle hypothesis DETRAY_HOST_DEVICE const auto &particle_hypothesis() const { return m_ptc; } @@ -253,7 +199,7 @@ class base_stepper { matrix_operator().set_identity(m_jac_transport); } - /// @returns access to this states step constraints + /// @returns access to this states navigation policy state DETRAY_HOST_DEVICE inline typename policy_t::state &policy_state() { return m_policy_state; @@ -281,8 +227,22 @@ class base_stepper { } private: - /// Stepping direction - step::direction m_direction{step::direction::e_forward}; + /// Jacobian transport matrix + free_matrix_type m_jac_transport = + matrix_operator().template identity(); + + /// Full jacobian + bound_matrix_type m_full_jacobian = + matrix_operator().template identity(); + + /// Bound covariance + bound_track_parameters_type m_bound_params; + + /// Free track parameters + free_track_parameters_type m_track; + + /// The default particle hypothesis is 'muon' + pdg_particle m_ptc = muon(); /// Total number of step trials std::size_t m_n_total_trials{0u}; @@ -290,50 +250,20 @@ class base_stepper { /// Current step size scalar_type m_step_size{0.f}; - /// Previous step size - scalar_type m_prev_step_size{0.f}; - - /// Track path length + /// Track path length (current position along track) scalar_type m_path_length{0.f}; - /// Absolute path length + /// Absolute path length (total path length covered by the integration) scalar_type m_abs_path_length{0.f}; - /// Track path length from the last surface. It will be reset to 0 when - /// the track reaches a new surface - scalar_type m_path_from_sf{0.f}; - - /// The default particle hypothesis is muon - pdg_particle m_ptc = muon(); - - /// Volume material that track is passing through - const material *m_vol_mat{nullptr}; - - /// Previous surface index to calculate the bound_to_free_jacobian - dindex m_prev_sf_id = detail::invalid_value(); - - /// free track parameter - free_track_parameters_type m_track; - - /// Full jacobian - bound_matrix_type m_full_jacobian = - matrix_operator().template identity(); - - /// jacobian transport matrix - free_matrix_type m_jac_transport = - matrix_operator().template identity(); - - /// bound covariance - bound_track_parameters_type m_bound_params; - - /// Step size constraints - constraint_t m_constraint = {}; + /// Step size constraints (optional) + [[no_unique_address]] constraint_t m_constraint = {}; /// Navigation policy state - typename policy_t::state m_policy_state = {}; + [[no_unique_address]] typename policy_t::state m_policy_state = {}; - /// The inspector type of the stepping - inspector_type m_inspector; + /// The inspector type of the stepping (for debugging only) + [[no_unique_address]] inspector_type m_inspector; }; }; diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index dd417bbbe..df098cfd1 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -78,7 +78,9 @@ class line_stepper final constexpr vector3_type dtds() const { return {0.f, 0.f, 0.f}; } DETRAY_HOST_DEVICE - constexpr scalar_type dqopds() const { return 0.f; } + constexpr scalar_type dqopds(const material*) const { + return 0.f; + } }; /// Take a step, regulared by a constrained step @@ -90,17 +92,12 @@ class line_stepper final /// @returns returning the heartbeat, indicating if the stepping is alive DETRAY_HOST_DEVICE bool step(const scalar_type dist_to_next, state& stepping, const stepping::config& cfg, - const bool = true) const { + const bool = true, + const material* = nullptr) const { // Straight line stepping: The distance given by the navigator is exact stepping.set_step_size(dist_to_next); - // Update navigation direction - const step::direction step_dir = stepping.step_size() >= 0.f - ? step::direction::e_forward - : step::direction::e_backward; - stepping.set_direction(step_dir); - // Check constraints if (const scalar_type max_step = stepping.constraints().template size<>(stepping.direction()); diff --git a/core/include/detray/propagator/propagator.hpp b/core/include/detray/propagator/propagator.hpp index 35f9ca9c0..d66dd0e5f 100644 --- a/core/include/detray/propagator/propagator.hpp +++ b/core/include/detray/propagator/propagator.hpp @@ -172,16 +172,16 @@ struct propagator { // Set access to the volume material for the stepper auto vol = navigation.get_volume(); - stepping.set_volume_material(vol.has_material() - ? vol.material_parameters(track.pos()) - : nullptr); + const material *vol_mat_ptr = + vol.has_material() ? vol.material_parameters(track.pos()) : nullptr; // Break automatic step size scaling by the stepper when a surface // was reached and whenever the navigation is (re-)initialized const bool reset_stepsize{navigation.is_on_surface() || is_init}; // Take the step - propagation._heartbeat &= m_stepper.step( - navigation(), stepping, m_cfg.stepping, reset_stepsize); + propagation._heartbeat &= + m_stepper.step(navigation(), stepping, m_cfg.stepping, + reset_stepsize, vol_mat_ptr); // Reduce navigation trust level according to stepper update typename stepper_t::policy_type{}(stepping.policy_state(), propagation); @@ -278,16 +278,17 @@ struct propagator { // Set access to the volume material for the stepper auto vol = navigation.get_volume(); - stepping.set_volume_material( + const material *vol_mat_ptr = vol.has_material() ? vol.material_parameters(track.pos()) - : nullptr); + : nullptr; // Break automatic step size scaling by the stepper const bool reset_stepsize{navigation.is_on_surface() || is_init}; // Take the step - propagation._heartbeat &= m_stepper.step( - navigation(), stepping, m_cfg.stepping, reset_stepsize); + propagation._heartbeat &= + m_stepper.step(navigation(), stepping, m_cfg.stepping, + reset_stepsize, vol_mat_ptr); // Reduce navigation trust level according to stepper update typename stepper_t::policy_type{}(stepping.policy_state(), @@ -381,7 +382,7 @@ struct propagator { } propagation.debug_stream << "step_size: " << std::setw(10) - << stepping.prev_step_size() << std::endl; + << stepping.step_size() << std::endl; propagation.debug_stream << std::setw(10) << detail::ray(stepping()) diff --git a/core/include/detray/propagator/rk_stepper.hpp b/core/include/detray/propagator/rk_stepper.hpp index b17787f03..978cac834 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -86,20 +86,34 @@ class rk_stepper final /// @returns the B-field view magnetic_field_type field() const { return m_magnetic_field; } + /// Set the next step size + DETRAY_HOST_DEVICE + inline void set_next_step_size(const scalar_type step) { + m_next_step_size = step; + } + + /// @returns the next step size to be taken on the following step. + DETRAY_HOST_DEVICE + inline scalar_type next_step_size() const { return m_next_step_size; } + /// Update the track state by Runge-Kutta-Nystrom integration. DETRAY_HOST_DEVICE - void advance_track(const intermediate_state& sd); + void advance_track(const intermediate_state& sd, + const material* vol_mat_ptr); /// Update the jacobian transport from free propagation DETRAY_HOST_DEVICE void advance_jacobian(const stepping::config& cfg, - const intermediate_state&); + const intermediate_state&, + const material* vol_mat_ptr); /// evaulate dqopds for a given step size and material DETRAY_HOST_DEVICE detray::pair evaluate_dqopds( const std::size_t i, const scalar_type h, - const scalar_type dqopds_prev, const detray::stepping::config& cfg); + const scalar_type dqopds_prev, + const material* vol_mat_ptr, + const detray::stepping::config& cfg); /// evaulate dtds for runge kutta stepping DETRAY_HOST_DEVICE @@ -117,14 +131,16 @@ class rk_stepper final /// Evaulate d(qop)/ds DETRAY_HOST_DEVICE - scalar_type dqopds() const; + scalar_type dqopds(const material* vol_mat_ptr) const; DETRAY_HOST_DEVICE - scalar_type dqopds(const scalar_type qop) const; + scalar_type dqopds(const scalar_type qop, + const material* vol_mat_ptr) const; /// Evaulate d(d(qop)/ds)dqop DETRAY_HOST_DEVICE - scalar_type d2qopdsdqop(const scalar_type qop) const; + scalar_type d2qopdsdqop(const scalar_type qop, + const material* vol_mat_ptr) const; /// Call the stepping inspector template @@ -143,6 +159,9 @@ class rk_stepper final vector3_type m_dtds_3; scalar_type m_dqopds_3; + /// Next step size after adaptive step size scaling + scalar_type m_next_step_size{0.f}; + /// Magnetic field view const magnetic_field_t m_magnetic_field; }; @@ -155,9 +174,10 @@ class rk_stepper final /// @param do_reset whether to reset the RKN step size to "dist to next" /// /// @return returning the heartbeat, indicating if the stepping is alive - DETRAY_HOST_DEVICE bool step(const scalar_type dist_to_next, - state& stepping, const stepping::config& cfg, - bool do_reset) const; + DETRAY_HOST_DEVICE bool step( + const scalar_type dist_to_next, state& stepping, + const stepping::config& cfg, bool do_reset, + const material* vol_mat_ptr = nullptr) const; }; } // namespace detray diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 871605744..a23d04e88 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -14,8 +14,8 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t>::state:: advance_track( const detray::rk_stepper::intermediate_state& - sd) { + policy_t, inspector_t>::intermediate_state& sd, + const material* vol_mat_ptr) { const scalar_type h{this->step_size()}; const scalar_type h_6{h * static_cast(1. / 6.)}; @@ -35,7 +35,7 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< track.set_dir(dir); auto qop = track.qop(); - if (this->volume_has_material()) { + if (vol_mat_ptr != nullptr) { // Reference: Eq (82) of https://doi.org/10.1016/0029-554X(81)90063-X qop = qop + h_6 * (sd.dqopds[0u] + 2.f * (sd.dqopds[1u] + sd.dqopds[2u]) + @@ -52,7 +52,9 @@ template ::state::advance_jacobian(const detray::stepping::config& cfg, - const intermediate_state& sd) { + const intermediate_state& sd, + const material* + vol_mat_ptr) { /// The calculations are based on ATL-SOFT-PUB-2009-002. The update of the /// Jacobian matrix is requires only the calculation of eq. 17 and 18. /// Since the terms of eq. 18 are currently 0, this matrix is not needed @@ -186,21 +188,22 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< getter::element(D, e_free_qoverp, e_free_qoverp) = 1.f; } else { // Pre-calculate dqop_n/dqop1 - const scalar_type d2qop1dsdqop1 = this->d2qopdsdqop(sd.qop[0u]); + const scalar_type d2qop1dsdqop1 = + this->d2qopdsdqop(sd.qop[0u], vol_mat_ptr); dqopn_dqop[0u] = 1.f; dqopn_dqop[1u] = 1.f + half_h * d2qop1dsdqop1; const scalar_type d2qop2dsdqop1 = - this->d2qopdsdqop(sd.qop[1u]) * dqopn_dqop[1u]; + this->d2qopdsdqop(sd.qop[1u], vol_mat_ptr) * dqopn_dqop[1u]; dqopn_dqop[2u] = 1.f + half_h * d2qop2dsdqop1; const scalar_type d2qop3dsdqop1 = - this->d2qopdsdqop(sd.qop[2u]) * dqopn_dqop[2u]; + this->d2qopdsdqop(sd.qop[2u], vol_mat_ptr) * dqopn_dqop[2u]; dqopn_dqop[3u] = 1.f + h * d2qop3dsdqop1; const scalar_type d2qop4dsdqop1 = - this->d2qopdsdqop(sd.qop[3u]) * dqopn_dqop[3u]; + this->d2qopdsdqop(sd.qop[3u], vol_mat_ptr) * dqopn_dqop[3u]; /*----------------------------------------------------------------- * Calculate the first terms of d(dqop_n/ds)/dqop1 @@ -331,16 +334,16 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< template DETRAY_HOST_DEVICE inline auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, - inspector_t>::state::evaluate_dqopds(const std::size_t i, - const scalar_type h, - const scalar_type dqopds_prev, - const detray::stepping::config& cfg) - -> detray::pair { + magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t>::state:: + evaluate_dqopds(const std::size_t i, const scalar_type h, + const scalar_type dqopds_prev, + const material* vol_mat_ptr, + const detray::stepping::config& cfg) + -> detray::pair { const auto& track = (*this)(); - if (!(this->volume_has_material())) { + if (!vol_mat_ptr) { const scalar_type qop = track.qop(); return detray::make_pair(scalar_type(0.f), qop); } else if (cfg.use_mean_loss && i != 0u) { @@ -351,10 +354,10 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< // "For y we have similar formulae as for x, for y' and // \lambda similar formulae as for x'" const scalar_type qop = track.qop() + h * dqopds_prev; - return detray::make_pair(this->dqopds(qop), qop); + return detray::make_pair(this->dqopds(qop, vol_mat_ptr), qop); } else { const scalar_type qop = track.qop(); - return detray::make_pair(this->dqopds(qop), qop); + return detray::make_pair(this->dqopds(qop, vol_mat_ptr), qop); } } @@ -442,13 +445,14 @@ detray::rk_stepper -DETRAY_HOST_DEVICE inline auto -detray::rk_stepper::state::dqopds() const -> scalar_type { +DETRAY_HOST_DEVICE inline auto detray::rk_stepper< + magnetic_field_t, algebra_t, constraint_t, policy_t, + inspector_t>::state::dqopds(const material* vol_mat_ptr) const + -> scalar_type { // In case there was no step before if (this->path_length() == 0.f) { - return this->dqopds((*this)().qop()); + return this->dqopds((*this)().qop(), vol_mat_ptr); } return m_dqopds_3; @@ -456,18 +460,17 @@ detray::rk_stepper -DETRAY_HOST_DEVICE auto -detray::rk_stepper::state::dqopds(const scalar_type qop) const +DETRAY_HOST_DEVICE auto detray::rk_stepper< + magnetic_field_t, algebra_t, constraint_t, policy_t, + inspector_t>::state::dqopds(const scalar_type qop, + const material* vol_mat_ptr) const -> scalar_type { // d(qop)ds is zero for empty space - if (!(this->volume_has_material())) { + if (!vol_mat_ptr) { return 0.f; } - const auto& mat = this->volume_material(); - const scalar_type q = this->particle_hypothesis().charge(); const scalar_type p = q / qop; const scalar_type mass = this->particle_hypothesis().mass(); @@ -476,7 +479,7 @@ detray::rk_stepper().compute_stopping_power( - mat, this->particle_hypothesis(), {mass, qop, q}); + *vol_mat_ptr, this->particle_hypothesis(), {mass, qop, q}); // Assert that a momentum is a positive value assert(p >= 0.f); @@ -490,14 +493,15 @@ template DETRAY_HOST_DEVICE auto detray::rk_stepper::state::d2qopdsdqop(const scalar_type qop) const + inspector_t>::state::d2qopdsdqop(const scalar_type qop, + const material* + vol_mat_ptr) const -> scalar_type { - if (!(this->volume_has_material())) { + if (!vol_mat_ptr) { return 0.f; } - const auto& mat = this->volume_material(); const scalar_type q = this->particle_hypothesis().charge(); const scalar_type p = q / qop; const scalar_type p2 = p * p; @@ -511,14 +515,16 @@ detray::rk_stepper rq(mass, qop, q); const scalar_type g = - -1.f * I.compute_stopping_power(mat, this->particle_hypothesis(), rq); + -1.f * + I.compute_stopping_power(*vol_mat_ptr, this->particle_hypothesis(), rq); // dg/d(qop) = -1 * derivation of stopping power const scalar_type dgdqop = - -1.f * I.derive_stopping_power(mat, this->particle_hypothesis(), rq); + -1.f * + I.derive_stopping_power(*vol_mat_ptr, this->particle_hypothesis(), rq); // d(qop)/ds = - qop^3 * E * g / q^2 - const scalar_type dqopds = this->dqopds(qop); + const scalar_type dqopds = this->dqopds(qop, vol_mat_ptr); // Check Eq 3.12 of // (https://iopscience.iop.org/article/10.1088/1748-0221/4/04/P04016/meta) @@ -532,7 +538,8 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< step(const scalar_type dist_to_next, detray::rk_stepper::state& stepping, - const detray::stepping::config& cfg, const bool do_reset) const { + const detray::stepping::config& cfg, const bool do_reset, + const material* vol_mat_ptr) const { // Get stepper and navigator states auto& magnetic_field = stepping.m_magnetic_field; @@ -540,17 +547,17 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< if (do_reset) { stepping.set_step_size(dist_to_next); } else if (stepping.step_size() > 0) { - stepping.set_step_size(math::min(stepping.step_size(), dist_to_next)); + stepping.set_step_size( + math::min(stepping.next_step_size(), dist_to_next)); } else { - stepping.set_step_size(math::max(stepping.step_size(), dist_to_next)); + stepping.set_step_size( + math::max(stepping.next_step_size(), dist_to_next)); } const point3_type pos = stepping().pos(); intermediate_state sd{}; - scalar_type error_estimate{0.f}; - // First Runge-Kutta point const auto bvec = magnetic_field.at(pos[0], pos[1], pos[2]); sd.b_first[0] = bvec[0]; @@ -560,7 +567,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< // qop should be recalcuated at every point // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X detray::tie(sd.dqopds[0u], sd.qop[0u]) = - stepping.evaluate_dqopds(0u, 0.f, 0.f, cfg); + stepping.evaluate_dqopds(0u, 0.f, 0.f, vol_mat_ptr, cfg); detray::tie(sd.dtds[0u], sd.t[0u]) = stepping.evaluate_dtds( sd.b_first, 0u, 0.f, vector3_type{0.f, 0.f, 0.f}, sd.qop[0u]); @@ -580,16 +587,16 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< sd.b_middle[1] = bvec1[1]; sd.b_middle[2] = bvec1[2]; - detray::tie(sd.dqopds[1u], sd.qop[1u]) = - stepping.evaluate_dqopds(1u, half_h, sd.dqopds[0u], cfg); + detray::tie(sd.dqopds[1u], sd.qop[1u]) = stepping.evaluate_dqopds( + 1u, half_h, sd.dqopds[0u], vol_mat_ptr, cfg); detray::tie(sd.dtds[1u], sd.t[1u]) = stepping.evaluate_dtds( sd.b_middle, 1u, half_h, sd.dtds[0u], sd.qop[1u]); // Third Runge-Kutta point // qop should be recalcuated at every point // Reference: Eq (84) of https://doi.org/10.1016/0029-554X(81)90063-X - detray::tie(sd.dqopds[2u], sd.qop[2u]) = - stepping.evaluate_dqopds(2u, half_h, sd.dqopds[1u], cfg); + detray::tie(sd.dqopds[2u], sd.qop[2u]) = stepping.evaluate_dqopds( + 2u, half_h, sd.dqopds[1u], vol_mat_ptr, cfg); detray::tie(sd.dtds[2u], sd.t[2u]) = stepping.evaluate_dtds( sd.b_middle, 2u, half_h, sd.dtds[1u], sd.qop[2u]); @@ -603,7 +610,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< sd.b_last[2] = bvec2[2]; detray::tie(sd.dqopds[3u], sd.qop[3u]) = - stepping.evaluate_dqopds(3u, h, sd.dqopds[2u], cfg); + stepping.evaluate_dqopds(3u, h, sd.dqopds[2u], vol_mat_ptr, cfg); detray::tie(sd.dtds[3u], sd.t[3u]) = stepping.evaluate_dtds(sd.b_last, 3u, h, sd.dtds[2u], sd.qop[3u]); @@ -614,9 +621,8 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< const vector3_type err_vec = one_sixth * h2 * (sd.dtds[0u] - sd.dtds[1u] - sd.dtds[2u] + sd.dtds[3u]); - error_estimate = getter::norm(err_vec); - return error_estimate; + return getter::norm(err_vec); }; /// Calculate the scale factor for the stepsize adjustment using the error @@ -633,7 +639,8 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< // If the estimated error is larger than the tolerance with an additional // margin, reduce the step size and try again - for (unsigned int i = 0u; i < cfg.max_rk_updates; i++) { + const auto n_trials{cfg.max_rk_updates}; + for (unsigned int i = 0u; i < n_trials; i++) { stepping.count_trials(); error = math::max(estimate_error(stepping.step_size()), @@ -659,12 +666,6 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< stepping.m_dtds_3 = sd.dtds[3u]; stepping.m_dqopds_3 = sd.dqopds[3u]; - // Update navigation direction - const step::direction step_dir = stepping.step_size() >= 0.f - ? step::direction::e_forward - : step::direction::e_backward; - stepping.set_direction(step_dir); - // Check constraints if (const scalar_type max_step = stepping.constraints().template size<>(stepping.direction()); @@ -677,22 +678,20 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< } // Advance track state - stepping.advance_track(sd); + stepping.advance_track(sd, vol_mat_ptr); // Advance jacobian transport if (cfg.do_covariance_transport) { - stepping.advance_jacobian(cfg, sd); + stepping.advance_jacobian(cfg, sd, vol_mat_ptr); } - // Save the current step size - stepping.set_prev_step_size(stepping.step_size()); - - // Update the step size for the next step - stepping.set_step_size(stepping.step_size() * step_size_scaling(error)); + // The step size estimation fot the next step + stepping.set_next_step_size(stepping.step_size() * + step_size_scaling(error)); // Don't allow a too small step size (unless requested by the navigator) - if (stepping.step_size() < cfg.min_stepsize) { - stepping.set_step_size(cfg.min_stepsize); + if (stepping.next_step_size() < cfg.min_stepsize) { + stepping.set_next_step_size(cfg.min_stepsize); } // Run final inspection diff --git a/tests/integration_tests/cpu/material/material_interaction.cpp b/tests/integration_tests/cpu/material/material_interaction.cpp index 610a8aeb0..ab1d9e89e 100644 --- a/tests/integration_tests/cpu/material/material_interaction.cpp +++ b/tests/integration_tests/cpu/material/material_interaction.cpp @@ -165,56 +165,6 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { EXPECT_NEAR(new_var_qop, dvar_qop, 1e-10f); - /******************************** - * Test with next_surface_aborter - *********************************/ - - using alt_actor_chain_t = - actor_chain, - next_surface_aborter, interactor_t, - parameter_resetter>; - using alt_propagator_t = - propagator; - - bound_track_parameters alt_bound_param( - geometry::barcode{}.set_index(0u), bound_vector, bound_cov); - - scalar altE(0); - - unsigned int surface_count = 0; - while (surface_count < 1e4) { - surface_count++; - - // Create actor states tuples - pathlimit_aborter::state alt_aborter_state{}; - next_surface_aborter::state next_surface_aborter_state{ - 0.1f * unit::mm}; - - auto alt_actor_states = detray::tie( - alt_aborter_state, bound_updater, next_surface_aborter_state, - interactor_state, parameter_resetter_state); - - // Propagator and its state - alt_propagator_t alt_p{prop_cfg}; - alt_propagator_t::state alt_state(alt_bound_param, det); - - // Propagate - alt_p.propagate(alt_state, alt_actor_states); - - alt_bound_param = alt_state._stepping.bound_params(); - - // Terminate the propagation if the next sensitive surface was not found - if (!next_surface_aborter_state.success) { - const scalar altP = - alt_state._stepping.bound_params().p(ptc.charge()); - altE = std::hypot(altP, mass); - break; - } - } - - EXPECT_EQ(surface_count, positions.size()); - EXPECT_EQ(altE, newE); - // @todo: Validate the backward direction case as well? } diff --git a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp index 28f6c75e0..e63781fa6 100644 --- a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp +++ b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp @@ -406,9 +406,9 @@ struct bound_getter : actor { const scalar N = static_cast(actor_state.step_count); - actor_state.m_avg_step_size = ((N - 1.f) * actor_state.m_avg_step_size + - stepping.prev_step_size()) / - N; + actor_state.m_avg_step_size = + ((N - 1.f) * actor_state.m_avg_step_size + stepping.step_size()) / + N; // Warning for too many step counts if (actor_state.step_count > 1000000) { diff --git a/tests/integration_tests/cpu/propagator/propagator.cpp b/tests/integration_tests/cpu/propagator/propagator.cpp index 3e0296910..5f0d731ec 100644 --- a/tests/integration_tests/cpu/propagator/propagator.cpp +++ b/tests/integration_tests/cpu/propagator/propagator.cpp @@ -55,6 +55,7 @@ struct helix_inspector : actor { struct state { /// navigation status for every step std::vector _nav_status; + scalar path_from_surface{0.f}; }; using matrix_operator = test::matrix_operator; @@ -67,15 +68,18 @@ struct helix_inspector : actor { const auto& navigation = prop_state._navigation; const auto& stepping = prop_state._stepping; - 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{}; + // Update inspector state inspector_state._nav_status.push_back(navigation.status()); + // The propagation does not start on a surface, skipp the inital path + if (!stepping.bound_params().surface_link().is_invalid()) { + inspector_state.path_from_surface += stepping.step_size(); + } // Nothing has happened yet (first call of actor chain) if (stepping.path_length() < tol || - stepping.path_from_surface() < tol) { + inspector_state.path_from_surface < tol) { return; } @@ -98,23 +102,28 @@ struct helix_inspector : actor { detail::helix hlx(free_params, &b); - const auto true_pos = hlx(stepping.path_from_surface()); + const auto true_pos = hlx(inspector_state.path_from_surface); - const point3 relative_error{1.f / stepping.path_from_surface() * + const point3 relative_error{1.f / inspector_state.path_from_surface * (stepping().pos() - true_pos)}; ASSERT_NEAR(getter::norm(relative_error), 0.f, tol); - auto true_J = hlx.jacobian(stepping.path_from_surface()); + auto true_J = hlx.jacobian(inspector_state.path_from_surface); for (unsigned int i = 0u; i < e_free_size; i++) { for (unsigned int j = 0u; j < e_free_size; j++) { ASSERT_NEAR(matrix_operator().element( stepping.transport_jacobian(), i, j), matrix_operator().element(true_J, i, j), - stepping.path_from_surface() * tol * 10.f); + inspector_state.path_from_surface * tol * 10.f); } } + // Reset path from surface + if (navigation.is_on_sensitive() || + navigation.encountered_sf_material()) { + inspector_state.path_from_surface = 0.f; + } } }; @@ -230,44 +239,47 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { track_t lim_track(track); // Build actor states: the helix inspector can be shared - helix_inspector::state helix_insp_state{}; - pathlimit_aborter::state unlimted_aborter_state{}; - pathlimit_aborter::state pathlimit_aborter_state{path_limit}; - parameter_transporter::state transporter_state{}; - pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; + auto actor_states = actor_chain_t::make_actor_states(); + auto actor_states_lim = actor_chain_t::make_actor_states(); + auto actor_states_sync = actor_chain_t::make_actor_states(); - // Create actor states tuples - auto actor_states = - detray::tie(helix_insp_state, unlimted_aborter_state, - transporter_state, interactor_state, resetter_state); - auto sync_actor_states = - detray::tie(helix_insp_state, unlimted_aborter_state, - transporter_state, interactor_state, resetter_state); - auto lim_actor_states = - detray::tie(helix_insp_state, pathlimit_aborter_state, - transporter_state, interactor_state, resetter_state); + // Make sure the lim state is being terminated + auto& pathlimit_aborter_state = + detail::get(actor_states_lim); + pathlimit_aborter_state.set_path_limit(path_limit); // Init propagator states propagator_t::state state(track, bfield, det); + propagator_t::state sync_state(track, bfield, det); propagator_t::state lim_state(lim_track, bfield, det); + state.do_debug = true; + sync_state.do_debug = true; + lim_state.do_debug = true; + // Set step constraints state._stepping.template set_constraint( step_constr); + sync_state._stepping + .template set_constraint(step_constr); lim_state._stepping .template set_constraint(step_constr); // Propagate the entire detector - state.do_debug = true; - ASSERT_TRUE(p.propagate(state, actor_states)) + ASSERT_TRUE( + p.propagate(state, actor_chain_t::make_ref_tuple(actor_states))) //<< state.debug_stream.str() << std::endl; << state._navigation.inspector().to_string() << std::endl; + // Test propagate sync method + ASSERT_TRUE(p.propagate_sync( + sync_state, actor_chain_t::make_ref_tuple(actor_states_sync))) + //<< state.debug_stream.str() << std::endl; + << sync_state._navigation.inspector().to_string() << std::endl; + // Propagate with path limit - ASSERT_NEAR(pathlimit_aborter_state.path_limit(), path_limit, tol); - lim_state.do_debug = true; - ASSERT_FALSE(p.propagate(lim_state, lim_actor_states)) + ASSERT_FALSE(p.propagate( + lim_state, actor_chain_t::make_ref_tuple(actor_states_lim))) //<< lim_state.debug_stream.str() << std::endl; << lim_state._navigation.inspector().to_string() << std::endl; @@ -279,9 +291,9 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { // Compare the navigation status vector between propagate and // propagate_sync function const auto nav_status = - detray::get(actor_states)._nav_status; + detray::get(actor_states)._nav_status; const auto sync_nav_status = - detray::get(sync_actor_states)._nav_status; + detray::get(actor_states_sync)._nav_status; ASSERT_TRUE(nav_status.size() > 0); ASSERT_TRUE(nav_status == sync_nav_status); } diff --git a/tests/unit_tests/cpu/propagator/rk_stepper.cpp b/tests/unit_tests/cpu/propagator/rk_stepper.cpp index cc241ca77..755b31317 100644 --- a/tests/unit_tests/cpu/propagator/rk_stepper.cpp +++ b/tests/unit_tests/cpu/propagator/rk_stepper.cpp @@ -248,19 +248,17 @@ TEST(detray_propagator, qop_derivative) { // RK Stepping into forward direction rk_stepper_t::state rk_state{track, hom_bfield}; - rk_state.set_volume_material(&vol_mat); - for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) { const scalar_t qop1 = rk_state().qop(); - const scalar_t d2qopdsdqop = rk_state.d2qopdsdqop(qop1); + const scalar_t d2qopdsdqop = rk_state.d2qopdsdqop(qop1, &vol_mat); - const scalar_t dqopds1 = rk_state.dqopds(qop1); + const scalar_t dqopds1 = rk_state.dqopds(qop1, &vol_mat); - rk_stepper.step(ds, rk_state, step_cfg, true); + rk_stepper.step(ds, rk_state, step_cfg, true, &vol_mat); const scalar_t qop2 = rk_state().qop(); - const scalar_t dqopds2 = rk_state.dqopds(qop2); + const scalar_t dqopds2 = rk_state.dqopds(qop2, &vol_mat); ASSERT_TRUE(qop1 > qop2); ASSERT_NEAR((qop2 - qop1) / ds, dqopds1, 1e-4); From bb3139fe9d2c9e241ea5a41192bdddd0d031e12c Mon Sep 17 00:00:00 2001 From: Joana Niermann Date: Thu, 7 Nov 2024 19:17:32 +0100 Subject: [PATCH 2/2] Move the bound track parameters and full jacobian to parameter transporter and merge the actors that modify the bound track parameters --- .../propagator/actors/parameter_resetter.hpp | 45 ------- ..._transporter.hpp => parameter_updater.hpp} | 120 ++++++++++++++++-- .../actors/pointwise_material_interactor.hpp | 7 +- core/include/detray/propagator/base_actor.hpp | 1 + .../detray/propagator/base_stepper.hpp | 57 +-------- core/include/detray/propagator/propagator.hpp | 17 +-- tests/benchmarks/cpu/benchmark_propagator.cpp | 16 +-- .../cuda/benchmark_propagator_cuda_kernel.cu | 6 +- .../cuda/benchmark_propagator_cuda_kernel.hpp | 11 +- .../test/device/cuda/material_validation.cu | 21 ++- .../detray/test/device/propagator_test.hpp | 24 ++-- .../utils/simulation/random_scatterer.hpp | 13 +- .../validation/material_validation_utils.hpp | 105 +++++++++------ .../cpu/material/material_interaction.cpp | 30 ++--- .../cpu/propagator/jacobian_validation.cpp | 30 ++--- .../cpu/propagator/propagator.cpp | 99 +++++++++------ .../device/cuda/propagator_cuda_kernel.cu | 8 +- .../device/sycl/propagator_kernel.sycl | 10 +- .../cpu/propagator/covariance_transport.cpp | 21 +-- tutorials/src/device/cuda/propagation.hpp | 14 +- .../src/device/cuda/propagation_kernel.cu | 6 +- 21 files changed, 346 insertions(+), 315 deletions(-) delete mode 100644 core/include/detray/propagator/actors/parameter_resetter.hpp rename core/include/detray/propagator/actors/{parameter_transporter.hpp => parameter_updater.hpp} (51%) diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp deleted file mode 100644 index 79bd20177..000000000 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ /dev/null @@ -1,45 +0,0 @@ -/** Detray library, part of the ACTS project (R&D line) - * - * (c) 2022-2024 CERN for the benefit of the ACTS project - * - * Mozilla Public License Version 2.0 - */ - -#pragma once - -// Project include(s). -#include "detray/definitions/detail/qualifiers.hpp" -#include "detray/definitions/track_parametrization.hpp" -#include "detray/geometry/tracking_surface.hpp" -#include "detray/propagator/base_actor.hpp" -#include "detray/propagator/detail/jacobian_engine.hpp" - -namespace detray { - -template -struct parameter_resetter : actor { - - template - DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const { - - const auto& navigation = propagation._navigation; - auto& stepping = propagation._stepping; - - // Do covariance transport when the track is on surface - if (!(navigation.is_on_sensitive() || - navigation.encountered_sf_material())) { - return; - } - - typename propagator_state_t::detector_type::geometry_context ctx{}; - - // 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()); - - // Reset jacobian transport to identity matrix - stepping.reset_transport_jacobian(); - } -}; - -} // namespace detray diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_updater.hpp similarity index 51% rename from core/include/detray/propagator/actors/parameter_transporter.hpp rename to core/include/detray/propagator/actors/parameter_updater.hpp index 78f9c07a6..125f36892 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_updater.hpp @@ -13,6 +13,7 @@ #include "detray/geometry/tracking_surface.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/detail/jacobian_engine.hpp" +#include "detray/utils/curvilinear_frame.hpp" namespace detray { @@ -24,19 +25,81 @@ struct parameter_transporter : actor { using scalar_type = dscalar; // Transformation matching this struct using transform3_type = dtransform3D; + // The track parameters bound to the current sensitive/material surface + using bound_track_parameters_type = bound_track_parameters; // Matrix actor using matrix_operator = dmatrix_operator; // bound matrix type - using bound_matrix_t = bound_matrix; + using bound_matrix_type = bound_matrix; // Matrix type for bound to free jacobian using bound_to_free_matrix_t = bound_to_free_matrix; /// @} + struct state { + + friend parameter_transporter; + + state() = default; + + /// Start from free track parameters + DETRAY_HOST_DEVICE + explicit state(const free_track_parameters& free_params) { + + curvilinear_frame cf(free_params); + + // Set bound track parameters + m_bound_params.set_parameter_vector(cf.m_bound_vec); + + // A dummy covariance - should not be used + m_bound_params.set_covariance( + matrix_operator() + .template identity()); + + // An invalid barcode - should not be used + m_bound_params.set_surface_link(geometry::barcode{}); + } + + /// Start from bound track parameters + DETRAY_HOST_DEVICE + explicit state(const bound_track_parameters_type& bound_params) + : m_bound_params{bound_params} {} + + /// @returns bound track parameters - const access + DETRAY_HOST_DEVICE + bound_track_parameters_type& bound_params() { return m_bound_params; } + + /// @returns bound track parameters. + DETRAY_HOST_DEVICE + const bound_track_parameters_type& bound_params() const { + return m_bound_params; + } + + /// @returns the current full Jacbian. + DETRAY_HOST_DEVICE + inline const bound_matrix_type& full_jacobian() const { + return m_full_jacobian; + } + + private: + /// Set new full Jacbian. + DETRAY_HOST_DEVICE + inline void set_full_jacobian(const bound_matrix_type& jac) { + m_full_jacobian = jac; + } + + /// Full jacobian + bound_matrix_type m_full_jacobian = + matrix_operator().template identity(); + + /// bound covariance + bound_track_parameters_type m_bound_params; + }; + struct get_full_jacobian_kernel { template - DETRAY_HOST_DEVICE inline bound_matrix_t operator()( + DETRAY_HOST_DEVICE inline bound_matrix_type operator()( const mask_group_t& /*mask_group*/, const index_t& /*index*/, const transform3_type& trf3, const bound_to_free_matrix_t& bound_to_free_jacobian, @@ -73,7 +136,8 @@ struct parameter_transporter : actor { }; template - DETRAY_HOST_DEVICE void operator()(propagator_state_t& propagation) const { + DETRAY_HOST_DEVICE void operator()(state& transporter_state, + propagator_state_t& propagation) const { auto& stepping = propagation._stepping; const auto& navigation = propagation._navigation; @@ -89,7 +153,7 @@ struct parameter_transporter : actor { const auto sf = navigation.get_surface(); // Bound track params of departure surface - auto& bound_params = stepping.bound_params(); + auto& bound_params = transporter_state.bound_params(); // Covariance is transported only when the previous surface is an // actual tracking surface. (i.e. This disables the covariance transport @@ -107,17 +171,17 @@ struct parameter_transporter : actor { const auto vol_mat_ptr = vol.has_material() ? vol.material_parameters(stepping().pos()) : nullptr; - stepping.set_full_jacobian( + transporter_state.set_full_jacobian( sf.template visit_mask( 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() * bound_params.covariance() * - matrix_operator().transpose(stepping.full_jacobian()); + const bound_matrix_type new_cov = + transporter_state.full_jacobian() * bound_params.covariance() * + matrix_operator().transpose(transporter_state.full_jacobian()); - stepping.bound_params().set_covariance(new_cov); + bound_params.set_covariance(new_cov); } // Convert free to bound vector @@ -127,7 +191,43 @@ struct parameter_transporter : actor { // Set surface link bound_params.set_surface_link(sf.barcode()); } +}; + +/// Update the stepper state after the bound track parameters were updated +template +struct parameter_resetter : actor { + + template + DETRAY_HOST_DEVICE void operator()( + const parameter_transporter::state& transporter_state, + propagator_state_t& propagation) const { -}; // namespace detray + const auto& navigation = propagation._navigation; + auto& stepping = propagation._stepping; + + // Do covariance transport when the track is on surface + if (!(navigation.is_on_sensitive() || + navigation.encountered_sf_material())) { + return; + } + + typename propagator_state_t::detector_type::geometry_context ctx{}; + + // Update free params after bound params were changed by actors + const auto sf = navigation.get_surface(); + stepping() = + sf.bound_to_free_vector(ctx, transporter_state.bound_params()); + + // Reset jacobian transport to identity matrix + stepping.reset_transport_jacobian(); + } +}; + +/// Call actors that depend on the bound track parameters safely together with +/// the parameter transporter and parameter resetter +template +using parameter_updater = + composite_actor, + transporter_observers..., parameter_resetter>; } // namespace detray diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index dae616d22..031ce991a 100644 --- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp +++ b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp @@ -14,6 +14,7 @@ #include "detray/materials/detail/concepts.hpp" #include "detray/materials/detail/material_accessor.hpp" #include "detray/materials/interaction.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/ranges.hpp" @@ -127,7 +128,9 @@ struct pointwise_material_interactor : actor { template DETRAY_HOST_DEVICE inline void operator()( - state &interactor_state, propagator_state_t &prop_state) const { + state &interactor_state, + parameter_transporter::state &transporter_state, + propagator_state_t &prop_state) const { // @Todo: Make context part of propagation state using detector_type = typename propagator_state_t::detector_type; @@ -143,7 +146,7 @@ struct pointwise_material_interactor : actor { auto &stepping = prop_state._stepping; this->update(geo_context_type{}, stepping.particle_hypothesis(), - stepping.bound_params(), interactor_state, + transporter_state.bound_params(), interactor_state, static_cast(navigation.direction()), navigation.get_surface()); } diff --git a/core/include/detray/propagator/base_actor.hpp b/core/include/detray/propagator/base_actor.hpp index ccdf029b5..97a097589 100644 --- a/core/include/detray/propagator/base_actor.hpp +++ b/core/include/detray/propagator/base_actor.hpp @@ -49,6 +49,7 @@ class composite_actor final : public actor_impl_t { /// derived from another composition (final). using actor_type = actor_impl_t; using state = typename actor_type::state; + using observer_states = detray::tuple; /// Call to the implementation of the actor (the actor possibly being an /// observer itself) diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index bb39149b3..a6da74a44 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -17,7 +17,6 @@ #include "detray/propagator/constrained_step.hpp" #include "detray/propagator/stepping_config.hpp" #include "detray/tracks/tracks.hpp" -#include "detray/utils/curvilinear_frame.hpp" namespace detray { @@ -57,31 +56,16 @@ class base_stepper { /// @note It has to cast into a const track via the call operation. struct state { - /// Sets track parameters. - DETRAY_HOST_DEVICE - explicit state(const free_track_parameters_type &free_params) - : m_track(free_params) { - - curvilinear_frame cf(free_params); - - // Set bound track parameters - m_bound_params.set_parameter_vector(cf.m_bound_vec); - - // A dummy covariance - should not be used - m_bound_params.set_covariance( - matrix_operator() - .template identity()); - - // An invalid barcode - should not be used - m_bound_params.set_surface_link(geometry::barcode{}); - } + /// Construct state from free track parameters. + DETRAY_HOST_DEVICE explicit state( + const free_track_parameters_type &free_params) + : m_track{free_params} {} - /// Sets track parameters from bound track parameter. + /// Sets track parameters from bound track parameters. template DETRAY_HOST_DEVICE state( const bound_track_parameters_type &bound_params, - const detector_t &det) - : m_bound_params(bound_params) { + const detector_t &det) { // Departure surface const auto sf = tracking_surface{det, bound_params.surface_link()}; @@ -99,16 +83,6 @@ class base_stepper { DETRAY_HOST_DEVICE const free_track_parameters_type &operator()() const { return m_track; } - /// @returns bound track parameters - const access - DETRAY_HOST_DEVICE - bound_track_parameters_type &bound_params() { return m_bound_params; } - - /// @returns bound track parameters - non-const access - DETRAY_HOST_DEVICE - const bound_track_parameters_type &bound_params() const { - return m_bound_params; - } - /// Get stepping direction DETRAY_HOST_DEVICE inline step::direction direction() const { @@ -181,18 +155,6 @@ class base_stepper { return m_jac_transport; } - /// @returns the current full Jacbian. - DETRAY_HOST_DEVICE - inline const bound_matrix_type &full_jacobian() const { - return m_full_jacobian; - } - - /// Set new full Jacbian. - DETRAY_HOST_DEVICE - inline void set_full_jacobian(const bound_matrix_type &jac) { - m_full_jacobian = jac; - } - /// Reset transport Jacbian. DETRAY_HOST_DEVICE inline void reset_transport_jacobian() { @@ -231,13 +193,6 @@ class base_stepper { free_matrix_type m_jac_transport = matrix_operator().template identity(); - /// Full jacobian - bound_matrix_type m_full_jacobian = - matrix_operator().template identity(); - - /// Bound covariance - bound_track_parameters_type m_bound_params; - /// Free track parameters free_track_parameters_type m_track; diff --git a/core/include/detray/propagator/propagator.hpp b/core/include/detray/propagator/propagator.hpp index d66dd0e5f..8c3a50971 100644 --- a/core/include/detray/propagator/propagator.hpp +++ b/core/include/detray/propagator/propagator.hpp @@ -135,9 +135,9 @@ struct propagator { /// /// @note If the return value of this function is true, a propagation step /// can be taken afterwards. + template DETRAY_HOST_DEVICE void propagate_init( - state &propagation, - typename actor_chain_t::state actor_state_refs) const { + state &propagation, actor_states_t actor_state_refs) const { auto &navigation = propagation._navigation; auto &stepping = propagation._stepping; const auto &track = stepping(); @@ -163,9 +163,10 @@ struct propagator { /// /// @note If the return value of this function is true, another step can /// be taken afterwards. + template DETRAY_HOST_DEVICE bool propagate_step( state &propagation, bool is_init, - typename actor_chain_t::state actor_state_refs) const { + actor_states_t actor_state_refs) const { auto &navigation = propagation._navigation; auto &stepping = propagation._stepping; const auto &track = stepping(); @@ -223,9 +224,9 @@ struct propagator { /// @param actor_state_refs tuple containing refences to the actor states /// /// @return propagation success. - DETRAY_HOST_DEVICE bool propagate( - state &propagation, - typename actor_chain_t::state actor_state_refs) const { + template + DETRAY_HOST_DEVICE bool propagate(state &propagation, + actor_states_t actor_state_refs) const { propagate_init(propagation, actor_state_refs); bool is_init = true; @@ -259,9 +260,9 @@ struct propagator { /// @param actor_states the actor state /// /// @return propagation success. + template DETRAY_HOST_DEVICE bool propagate_sync( - state &propagation, - typename actor_chain_t::state actor_state_refs) const { + state &propagation, actor_states_t actor_state_refs) const { propagate_init(propagation, actor_state_refs); bool is_init = true; diff --git a/tests/benchmarks/cpu/benchmark_propagator.cpp b/tests/benchmarks/cpu/benchmark_propagator.cpp index ac3cc4cf6..735222992 100644 --- a/tests/benchmarks/cpu/benchmark_propagator.cpp +++ b/tests/benchmarks/cpu/benchmark_propagator.cpp @@ -16,8 +16,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/propagator.hpp" @@ -55,9 +54,11 @@ using navigator_host_type = navigator; using navigator_device_type = navigator; using field_type = bfield::const_field_t; using rk_stepper_type = rk_stepper; -using actor_chain_t = actor_chain, - pointwise_material_interactor, - parameter_resetter>; + +using parameter_updater_t = + parameter_updater>; + +using actor_chain_t = actor_chain; using propagator_host_type = propagator; using propagator_device_type = @@ -119,12 +120,11 @@ static void BM_PROPAGATOR_CPU(benchmark::State &state) { #pragma omp parallel for for (auto &track : tracks) { - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{track}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; auto actor_states = - tie(transporter_state, interactor_state, resetter_state); + detray::tie(transporter_state, interactor_state); // Create the propagator state propagator_host_type::state p_state(track, bfield, det); diff --git a/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.cu b/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.cu index 6bba06fd1..3eaeacbba 100644 --- a/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.cu +++ b/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.cu @@ -30,13 +30,11 @@ __global__ void __launch_bounds__(256, 4) propagator_benchmark_kernel( cfg.navigation.search_window = {3u, 3u}; propagator_device_type p{cfg}; - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{tracks.at(gid)}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; // Create the actor states - auto actor_states = - detray::tie(transporter_state, interactor_state, resetter_state); + auto actor_states = detray::tie(transporter_state, interactor_state); // Create the propagator state propagator_device_type::state p_state(tracks.at(gid), field_data, det); diff --git a/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.hpp b/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.hpp index ee2f505cb..a2886b502 100644 --- a/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.hpp +++ b/tests/benchmarks/cuda/benchmark_propagator_cuda_kernel.hpp @@ -15,8 +15,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/propagator.hpp" @@ -34,10 +33,10 @@ using navigator_host_type = detray::navigator; using navigator_device_type = detray::navigator; using field_type = detray::bfield::const_field_t; using rk_stepper_type = detray::rk_stepper; -using actor_chain_t = - detray::actor_chain, - detray::pointwise_material_interactor, - detray::parameter_resetter>; +using parameter_updater_t = + detray::parameter_updater>; +using actor_chain_t = detray::actor_chain; using propagator_host_type = detray::propagator; using propagator_device_type = diff --git a/tests/include/detray/test/device/cuda/material_validation.cu b/tests/include/detray/test/device/cuda/material_validation.cu index bad863fdc..aa5a90c7b 100644 --- a/tests/include/detray/test/device/cuda/material_validation.cu +++ b/tests/include/detray/test/device/cuda/material_validation.cu @@ -9,8 +9,7 @@ #include "detray/detectors/telescope_metadata.hpp" #include "detray/detectors/toy_metadata.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/line_stepper.hpp" #include "material_validation.hpp" @@ -41,11 +40,10 @@ __global__ void material_validation_kernel( // material tracer using material_tracer_t = material_validator::material_tracer; - using actor_chain_t = - actor_chain, - parameter_resetter, - pointwise_material_interactor, - material_tracer_t>; + using parameter_updater_t = + parameter_updater>; + using actor_chain_t = actor_chain; using propagator_t = propagator; detector_device_t det(det_data); @@ -66,14 +64,13 @@ __global__ void material_validation_kernel( // Create the actor states pathlimit_aborter::state aborter_state{cfg.stepping.path_limit}; - typename parameter_transporter::state transporter_state{}; - typename parameter_resetter::state resetter_state{}; + typename parameter_transporter::state transporter_state{ + tracks[trk_id]}; typename pointwise_material_interactor::state interactor_state{}; typename material_tracer_t::state mat_tracer_state{mat_steps.at(trk_id)}; - auto actor_states = - ::detray::tie(aborter_state, transporter_state, resetter_state, - interactor_state, mat_tracer_state); + auto actor_states = ::detray::tie(aborter_state, transporter_state, + interactor_state, mat_tracer_state); // Run propagation typename navigator_t::state::view_type nav_view{}; diff --git a/tests/include/detray/test/device/propagator_test.hpp b/tests/include/detray/test/device/propagator_test.hpp index 9a11c2691..53640d6ca 100644 --- a/tests/include/detray/test/device/propagator_test.hpp +++ b/tests/include/detray/test/device/propagator_test.hpp @@ -13,8 +13,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/propagator.hpp" @@ -80,16 +79,13 @@ struct propagator_test_config { // Assemble actor chain type using step_tracer_host_t = step_tracer; using step_tracer_device_t = step_tracer; -using actor_chain_host_t = - actor_chain, - pointwise_material_interactor, - parameter_resetter>; +using parameter_updater_t = + parameter_updater>; +using actor_chain_host_t = actor_chain; using actor_chain_device_t = actor_chain, - pointwise_material_interactor, - parameter_resetter>; + parameter_updater_t>; /// Precompute the tracks template > @@ -137,12 +133,10 @@ inline auto run_propagation_host(vecmem::memory_resource *mr, step_tracer_host_t::state tracer_state{*mr}; tracer_state.collect_only_on_surface(true); pathlimit_aborter::state pathlimit_state{cfg.stepping.path_limit}; - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{trk}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; - auto actor_states = - detray::tie(tracer_state, pathlimit_state, transporter_state, - interactor_state, resetter_state); + auto actor_states = detray::tie(tracer_state, pathlimit_state, + transporter_state, interactor_state); typename propagator_host_t::state state(trk, field, det); diff --git a/tests/include/detray/test/utils/simulation/random_scatterer.hpp b/tests/include/detray/test/utils/simulation/random_scatterer.hpp index 592b30863..6b866a04c 100644 --- a/tests/include/detray/test/utils/simulation/random_scatterer.hpp +++ b/tests/include/detray/test/utils/simulation/random_scatterer.hpp @@ -15,6 +15,7 @@ #include "detray/geometry/tracking_surface.hpp" #include "detray/materials/detail/concepts.hpp" #include "detray/materials/interaction.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/tracks/bound_track_parameters.hpp" #include "detray/utils/axis_rotation.hpp" @@ -130,8 +131,10 @@ struct random_scatterer : actor { }; template - DETRAY_HOST inline void operator()(state& simulator_state, - propagator_state_t& prop_state) const { + DETRAY_HOST inline void operator()( + state& simulator_state, + parameter_transporter::state& transporter_state, + propagator_state_t& prop_state) const { // @Todo: Make context part of propagation state using detector_type = typename propagator_state_t::detector_type; @@ -145,7 +148,7 @@ struct random_scatterer : actor { auto& stepping = prop_state._stepping; const auto& ptc = stepping.particle_hypothesis(); - auto& bound_params = stepping.bound_params(); + auto& bound_params = transporter_state.bound_params(); const auto sf = navigation.get_surface(); const scalar_type cos_inc_angle{ sf.cos_angle(geo_context_type{}, bound_params.dir(), @@ -170,8 +173,8 @@ struct random_scatterer : actor { simulator_state.generator); // Update Phi and Theta - stepping.bound_params().set_phi(getter::phi(new_dir)); - stepping.bound_params().set_theta(getter::theta(new_dir)); + bound_params.set_phi(getter::phi(new_dir)); + bound_params.set_theta(getter::theta(new_dir)); // Flag renavigation of the current candidate prop_state._navigation.set_high_trust(); diff --git a/tests/include/detray/test/validation/material_validation_utils.hpp b/tests/include/detray/test/validation/material_validation_utils.hpp index 7934d2a64..256c32695 100644 --- a/tests/include/detray/test/validation/material_validation_utils.hpp +++ b/tests/include/detray/test/validation/material_validation_utils.hpp @@ -13,8 +13,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/line_stepper.hpp" @@ -152,26 +151,43 @@ struct material_tracer : detray::actor { vector_t> m_mat_steps{}; }; + /// Run as observer to the parameter transporter (covariance transportsa) + template + DETRAY_HOST_DEVICE void operator()( + state &tracer, + parameter_transporter::state &transporter_state, + const propagator_state_t &prop_state) const { + + record_track_dir(tracer, prop_state); + + // Only count material if navigator encountered it + const auto &navigation = prop_state._navigation; + if (!navigation.encountered_sf_material()) { + return; + } + + // For now use default context + typename propagator_state_t::detector_type::geometry_context gctx{}; + + const auto &track_param = transporter_state.bound_params(); + dpoint2D loc_pos = track_param.bound_local(); + + record_mat_step(tracer, gctx, navigation.get_surface(), loc_pos, + track_param.dir()); + } + + /// Run in a propagation chain without covariance transport template DETRAY_HOST_DEVICE void operator()( state &tracer, const propagator_state_t &prop_state) const { using algebra_t = typename propagator_state_t::detector_type::algebra_type; - using vector3_t = dvector3D; - using point2_t = dpoint2D; - const auto &navigation = prop_state._navigation; - - // Record the initial track direction - vector3_t glob_dir = prop_state._stepping().dir(); - if (detray::detail::is_invalid_value(tracer.m_mat_record.eta) && - detray::detail::is_invalid_value(tracer.m_mat_record.phi)) { - tracer.m_mat_record.eta = getter::eta(glob_dir); - tracer.m_mat_record.phi = getter::phi(glob_dir); - } + record_track_dir(tracer, prop_state); // Only count material if navigator encountered it + const auto &navigation = prop_state._navigation; if (!navigation.encountered_sf_material()) { return; } @@ -183,21 +199,39 @@ struct material_tracer : detray::actor { const auto sf = navigation.get_surface(); // Track direction and bound position on current surface - point2_t loc_pos{}; - - // Get the local track position from the bound track parameters, - // if covariance transport is enabled in the propagation - if constexpr (detail::has_type_v< - typename parameter_transporter::state &, - typename propagator_state_t::actor_chain_type:: - state>) { - const auto &track_param = prop_state._stepping.bound_params(); - loc_pos = track_param.bound_local(); - } else { - const auto &track_param = prop_state._stepping(); - glob_dir = track_param.dir(); - loc_pos = sf.global_to_bound(gctx, track_param.pos(), glob_dir); + const auto &track_param = prop_state._stepping(); + dvector3D glob_dir = track_param.dir(); + dpoint2D loc_pos = + sf.global_to_bound(gctx, track_param.pos(), glob_dir); + + record_mat_step(tracer, gctx, sf, loc_pos, glob_dir); + } + + private: + /// Record the track direction + template + DETRAY_HOST_DEVICE inline auto record_track_dir( + state &tracer, const propagator_state_t &prop_state) const { + using algebra_t = + typename propagator_state_t::detector_type::algebra_type; + using vector3_t = dvector3D; + + // Record the initial track direction + vector3_t glob_dir = prop_state._stepping().dir(); + if (detray::detail::is_invalid_value(tracer.m_mat_record.eta) && + detray::detail::is_invalid_value(tracer.m_mat_record.phi)) { + tracer.m_mat_record.eta = getter::eta(glob_dir); + tracer.m_mat_record.phi = getter::phi(glob_dir); } + } + + /// Record the data for a material step + template + DETRAY_HOST_DEVICE inline auto record_mat_step( + state &tracer, const typename detector_t::geometry_context &gctx, + const tracking_surface sf, + const dpoint2D &loc_pos, + const dvector3D &glob_dir) const { // Fetch the material parameters and pathlength through the material const auto mat_params = sf.template visit_material( @@ -241,11 +275,10 @@ inline auto record_material( // Propagator with pathlimit aborter using material_tracer_t = material_validator::material_tracer; - using actor_chain_t = - actor_chain, - parameter_resetter, - pointwise_material_interactor, - material_tracer_t>; + using parameter_updater_t = + parameter_updater>; + using actor_chain_t = actor_chain; using propagator_t = propagator; // Propagator @@ -253,14 +286,12 @@ inline auto record_material( // Build actor and propagator states pathlimit_aborter::state pathlimit_aborter_state{cfg.stepping.path_limit}; - typename parameter_transporter::state transporter_state{}; - typename parameter_resetter::state resetter_state{}; + typename parameter_transporter::state transporter_state{track}; typename pointwise_material_interactor::state interactor_state{}; typename material_tracer_t::state mat_tracer_state{*host_mr}; - auto actor_states = - detray::tie(pathlimit_aborter_state, transporter_state, resetter_state, - interactor_state, mat_tracer_state); + auto actor_states = detray::tie(pathlimit_aborter_state, transporter_state, + interactor_state, mat_tracer_state); typename propagator_t::state propagation{track, det}; diff --git a/tests/integration_tests/cpu/material/material_interaction.cpp b/tests/integration_tests/cpu/material/material_interaction.cpp index ab1d9e89e..b350f9ec3 100644 --- a/tests/integration_tests/cpu/material/material_interaction.cpp +++ b/tests/integration_tests/cpu/material/material_interaction.cpp @@ -18,8 +18,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/line_stepper.hpp" #include "detray/propagator/propagator.hpp" @@ -73,9 +72,9 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { using navigator_t = navigator; using stepper_t = line_stepper; using interactor_t = pointwise_material_interactor; + using parameter_updater_t = parameter_updater; using actor_chain_t = - actor_chain, - interactor_t, parameter_resetter>; + actor_chain; using propagator_t = propagator; // Propagator is built from the stepper and navigator @@ -99,13 +98,12 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { geometry::barcode{}.set_index(0u), bound_vector, bound_cov); pathlimit_aborter::state aborter_state{}; - parameter_transporter::state bound_updater{}; + parameter_transporter::state bound_updater{bound_param}; interactor_t::state interactor_state{}; - parameter_resetter::state parameter_resetter_state{}; // Create actor states tuples - auto actor_states = detray::tie(aborter_state, bound_updater, - interactor_state, parameter_resetter_state); + auto actor_states = + detray::tie(aborter_state, bound_updater, interactor_state); propagator_t::state state(bound_param, det); state.do_debug = true; @@ -115,7 +113,7 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { << state.debug_stream.str() << std::endl; // new momentum - const scalar newP{state._stepping.bound_params().p(ptc.charge())}; + const scalar newP{bound_updater.bound_params().p(ptc.charge())}; // mass const auto mass = ptc.mass(); @@ -128,7 +126,7 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { // New qop variance const scalar new_var_qop{ - matrix_operator().element(state._stepping.bound_params().covariance(), + matrix_operator().element(bound_updater.bound_params().covariance(), e_bound_qoverp, e_bound_qoverp)}; // Interaction object @@ -196,9 +194,9 @@ GTEST_TEST(detray_material, telescope_geometry_scattering_angle) { using navigator_t = navigator; using stepper_t = line_stepper; using simulator_t = random_scatterer; + using parameter_updater_t = parameter_updater; using actor_chain_t = - actor_chain, - simulator_t, parameter_resetter>; + actor_chain; using propagator_t = propagator; // Propagator is built from the stepper and navigator @@ -228,16 +226,14 @@ GTEST_TEST(detray_material, telescope_geometry_scattering_angle) { for (std::size_t i = 0u; i < n_samples; i++) { pathlimit_aborter::state aborter_state{}; - parameter_transporter::state bound_updater{}; + parameter_transporter::state bound_updater{bound_param}; // Seed = sample id simulator_t::state simulator_state{i}; simulator_state.do_energy_loss = false; - parameter_resetter::state parameter_resetter_state{}; // Create actor states tuples auto actor_states = - detray::tie(aborter_state, bound_updater, simulator_state, - parameter_resetter_state); + detray::tie(aborter_state, bound_updater, simulator_state); propagator_t::state state(bound_param, det); state.do_debug = true; @@ -246,7 +242,7 @@ GTEST_TEST(detray_material, telescope_geometry_scattering_angle) { ASSERT_TRUE(p.propagate(state, actor_states)) << state.debug_stream.str() << std::endl; - const auto& final_param = state._stepping.bound_params(); + const auto& final_param = bound_updater.bound_params(); // Updated phi and theta variance if (i == 0u) { diff --git a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp index e63781fa6..ddbccbdaa 100644 --- a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp +++ b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp @@ -10,8 +10,7 @@ #include "detray/detectors/bfield.hpp" #include "detray/navigation/intersection/helix_intersector.hpp" #include "detray/navigation/navigator.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/propagator.hpp" #include "detray/propagator/rk_stepper.hpp" @@ -396,8 +395,10 @@ struct bound_getter : actor { }; template - DETRAY_HOST_DEVICE void operator()(state& actor_state, - propagator_state_t& propagation) const { + DETRAY_HOST_DEVICE void operator()( + state& actor_state, + const parameter_transporter::state& transporter_state, + propagator_state_t& propagation) const { auto& navigation = propagation._navigation; auto& stepping = propagation._stepping; @@ -428,7 +429,7 @@ struct bound_getter : actor { if ((navigation.is_on_sensitive() || navigation.is_on_passive()) && navigation.barcode().index() == 0u) { - actor_state.m_param_departure = stepping.bound_params(); + actor_state.m_param_departure = transporter_state.bound_params(); } // Get the bound track parameters and jacobian at the destination // surface @@ -437,8 +438,8 @@ struct bound_getter : actor { actor_state.m_path_length = stepping.path_length(); actor_state.m_abs_path_length = stepping.abs_path_length(); - actor_state.m_param_destination = stepping.bound_params(); - actor_state.m_jacobi = stepping.full_jacobian(); + actor_state.m_param_destination = transporter_state.bound_params(); + actor_state.m_jacobi = transporter_state.full_jacobian(); // Stop navigation if the destination surface found propagation._heartbeat &= navigation.exit(); @@ -477,9 +478,7 @@ bound_getter::state evaluate_bound_param( bound_getter::state bound_getter_state{}; bound_getter_state.track_ID = trk_count; bound_getter_state.m_min_path_length = detector_length * 0.75f; - parameter_resetter::state resetter_state{}; - auto actor_states = - detray::tie(transporter_state, bound_getter_state, resetter_state); + auto actor_states = detray::tie(transporter_state, bound_getter_state); // Init propagator states for the reference track typename propagator_t::state state(initial_param, field, det); @@ -525,13 +524,11 @@ bound_param_vector_type get_displaced_bound_vector( // Actor states parameter_transporter::state transporter_state{}; - parameter_resetter::state resetter_state{}; bound_getter::state bound_getter_state{}; bound_getter_state.track_ID = trk_count; bound_getter_state.m_min_path_length = detector_length * 0.75f; - auto actor_states = - detray::tie(transporter_state, bound_getter_state, resetter_state); + auto actor_states = detray::tie(transporter_state, bound_getter_state); dstate.set_particle(ptc); dstate._stepping .template set_constraint( @@ -1569,10 +1566,9 @@ int main(int argc, char** argv) { const inhom_bfield_t inhom_bfield = bfield::create_inhom_field(); // Actor chain type - using actor_chain_t = - actor_chain, - bound_getter, - parameter_resetter>; + using parameter_updater_t = + parameter_updater>; + using actor_chain_t = actor_chain; // Iterate over reference (pilot) tracks for a rectangular telescope // geometry and Jacobian calculation diff --git a/tests/integration_tests/cpu/propagator/propagator.cpp b/tests/integration_tests/cpu/propagator/propagator.cpp index 5f0d731ec..81eb88d8e 100644 --- a/tests/integration_tests/cpu/propagator/propagator.cpp +++ b/tests/integration_tests/cpu/propagator/propagator.cpp @@ -15,8 +15,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/base_actor.hpp" #include "detray/propagator/line_stepper.hpp" @@ -63,19 +62,18 @@ struct helix_inspector : actor { /// Check that the stepper remains on the right helical track for its pos. template DETRAY_HOST_DEVICE void operator()( - state& inspector_state, const propagator_state_t& prop_state) const { + state& inspector_state, + parameter_transporter::state& transporter_state, + const propagator_state_t& prop_state) const { const auto& navigation = prop_state._navigation; const auto& stepping = prop_state._stepping; + const auto& bound_params = transporter_state.bound_params(); typename propagator_state_t::detector_type::geometry_context ctx{}; // Update inspector state inspector_state._nav_status.push_back(navigation.status()); - // The propagation does not start on a surface, skipp the inital path - if (!stepping.bound_params().surface_link().is_invalid()) { - inspector_state.path_from_surface += stepping.step_size(); - } // Nothing has happened yet (first call of actor chain) if (stepping.path_length() < tol || @@ -83,16 +81,16 @@ struct helix_inspector : actor { return; } - if (stepping.bound_params().surface_link().is_invalid()) { + if (bound_params.surface_link().is_invalid()) { return; } // Surface - const auto sf = tracking_surface{ - navigation.detector(), stepping.bound_params().surface_link()}; + const auto sf = tracking_surface{navigation.detector(), + bound_params.surface_link()}; const free_track_parameters free_params = - sf.bound_to_free_vector(ctx, stepping.bound_params()); + sf.bound_to_free_vector(ctx, bound_params); const auto last_pos = free_params.pos(); @@ -119,6 +117,10 @@ struct helix_inspector : actor { inspector_state.path_from_surface * tol * 10.f); } } + // The propagation does not start on a surface, skipp the inital path + if (!bound_params.surface_link().is_invalid()) { + inspector_state.path_from_surface += stepping.step_size(); + } // Reset path from surface if (navigation.is_on_sensitive() || navigation.encountered_sf_material()) { @@ -209,12 +211,15 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { using policy_t = stepper_rk_policy; using stepper_t = rk_stepper; - // Include helix actor to check track position/covariance + + // Parameter update scheme: Include material interaction and helix + // inspector to check track position/covariance + using parameter_updater_t = + parameter_updater, + helix_inspector>; + using actor_chain_t = - actor_chain, - pointwise_material_interactor, - parameter_resetter>; + actor_chain; using propagator_t = propagator; // Build detector @@ -238,15 +243,30 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { // Generate second track state used for propagation with pathlimit track_t lim_track(track); - // Build actor states: the helix inspector can be shared - auto actor_states = actor_chain_t::make_actor_states(); - auto actor_states_lim = actor_chain_t::make_actor_states(); - auto actor_states_sync = actor_chain_t::make_actor_states(); + // Build actor states + pathlimit_aborter::state unlimted_aborter_state{}; + parameter_transporter::state transporter_state{track}; + pointwise_material_interactor::state interactor_state{}; + helix_inspector::state helix_insp_state{}; + auto actor_states = + detray::tie(unlimted_aborter_state, transporter_state, + interactor_state, helix_insp_state); + + parameter_transporter::state transporter_state_sync{track}; + helix_inspector::state helix_insp_state_sync{}; + auto actor_states_sync = + detray::tie(unlimted_aborter_state, transporter_state_sync, + interactor_state, helix_insp_state_sync); // Make sure the lim state is being terminated - auto& pathlimit_aborter_state = - detail::get(actor_states_lim); + pathlimit_aborter::state pathlimit_aborter_state{path_limit}; pathlimit_aborter_state.set_path_limit(path_limit); + parameter_transporter::state transporter_state_lim{ + lim_track}; + helix_inspector::state helix_insp_state_lim{}; + auto actor_states_lim = + detray::tie(pathlimit_aborter_state, transporter_state_lim, + interactor_state, helix_insp_state_lim); // Init propagator states propagator_t::state state(track, bfield, det); @@ -266,20 +286,17 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { .template set_constraint(step_constr); // Propagate the entire detector - ASSERT_TRUE( - p.propagate(state, actor_chain_t::make_ref_tuple(actor_states))) + ASSERT_TRUE(p.propagate(state, actor_states)) //<< state.debug_stream.str() << std::endl; << state._navigation.inspector().to_string() << std::endl; // Test propagate sync method - ASSERT_TRUE(p.propagate_sync( - sync_state, actor_chain_t::make_ref_tuple(actor_states_sync))) + ASSERT_TRUE(p.propagate_sync(sync_state, actor_states_sync)) //<< state.debug_stream.str() << std::endl; << sync_state._navigation.inspector().to_string() << std::endl; // Propagate with path limit - ASSERT_FALSE(p.propagate( - lim_state, actor_chain_t::make_ref_tuple(actor_states_lim))) + ASSERT_FALSE(p.propagate(lim_state, actor_states_lim)) //<< lim_state.debug_stream.str() << std::endl; << lim_state._navigation.inspector().to_string() << std::endl; @@ -291,9 +308,9 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { // Compare the navigation status vector between propagate and // propagate_sync function const auto nav_status = - detray::get(actor_states)._nav_status; + detray::get(actor_states)._nav_status; const auto sync_nav_status = - detray::get(actor_states_sync)._nav_status; + detray::get(actor_states_sync)._nav_status; ASSERT_TRUE(nav_status.size() > 0); ASSERT_TRUE(nav_status == sync_nav_status); } @@ -317,11 +334,10 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_inhom_bfield) { using policy_t = stepper_rk_policy; using stepper_t = rk_stepper; - // Include helix actor to check track position/covariance + using parameter_updater_t = + parameter_updater>; using actor_chain_t = - actor_chain, - pointwise_material_interactor, - parameter_resetter>; + actor_chain; using propagator_t = propagator; // Build detector and magnetic field @@ -343,17 +359,16 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_inhom_bfield) { // Build actor states: the helix inspector can be shared pathlimit_aborter::state unlimted_aborter_state{}; pathlimit_aborter::state pathlimit_aborter_state{path_limit}; - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{track}; + parameter_transporter::state transporter_state_lim{ + lim_track}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; // Create actor states tuples - auto actor_states = - detray::tie(unlimted_aborter_state, transporter_state, - interactor_state, resetter_state); - auto lim_actor_states = - detray::tie(pathlimit_aborter_state, transporter_state, - interactor_state, resetter_state); + auto actor_states = detray::tie(unlimted_aborter_state, + transporter_state, interactor_state); + auto lim_actor_states = detray::tie( + pathlimit_aborter_state, transporter_state_lim, interactor_state); // Init propagator states propagator_t::state state(track, bfield, det); diff --git a/tests/integration_tests/device/cuda/propagator_cuda_kernel.cu b/tests/integration_tests/device/cuda/propagator_cuda_kernel.cu index d8e4b09bb..a13320c53 100644 --- a/tests/integration_tests/device/cuda/propagator_cuda_kernel.cu +++ b/tests/integration_tests/device/cuda/propagator_cuda_kernel.cu @@ -50,14 +50,12 @@ __global__ void propagator_test_kernel( step_tracer_device_t::state tracer_state(steps.at(gid)); tracer_state.collect_only_on_surface(true); pathlimit_aborter::state aborter_state{cfg.stepping.path_limit}; - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{tracks[gid]}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; // Create the actor states - auto actor_states = - ::detray::tie(tracer_state, aborter_state, transporter_state, - interactor_state, resetter_state); + auto actor_states = ::detray::tie(tracer_state, aborter_state, + transporter_state, interactor_state); // Create the propagator state typename propagator_device_t::state state(tracks[gid], field_data, det); diff --git a/tests/integration_tests/device/sycl/propagator_kernel.sycl b/tests/integration_tests/device/sycl/propagator_kernel.sycl index b51171ca3..86041a85e 100644 --- a/tests/integration_tests/device/sycl/propagator_kernel.sycl +++ b/tests/integration_tests/device/sycl/propagator_kernel.sycl @@ -65,15 +65,15 @@ void propagator_test( step_tracer_device_t::state tracer_state(steps.at(gid)); tracer_state.collect_only_on_surface(true); pathlimit_aborter::state aborter_state{cfg.stepping.path_limit}; - parameter_transporter::state transporter_state{}; + parameter_transporter::state transporter_state{ + tracks[gid]}; pointwise_material_interactor::state interactor_state{}; - parameter_resetter::state resetter_state{}; // Create the actor states - auto actor_states = ::detray::tie( - tracer_state, aborter_state, transporter_state, - interactor_state, resetter_state); + auto actor_states = + ::detray::tie(tracer_state, aborter_state, + transporter_state, interactor_state); // Create the propagator state typename propagator_device_t::state state(tracks[gid], field_data, dev_det); diff --git a/tests/unit_tests/cpu/propagator/covariance_transport.cpp b/tests/unit_tests/cpu/propagator/covariance_transport.cpp index e9323e47c..dd0b31330 100644 --- a/tests/unit_tests/cpu/propagator/covariance_transport.cpp +++ b/tests/unit_tests/cpu/propagator/covariance_transport.cpp @@ -14,8 +14,7 @@ #include "detray/navigation/detail/ray.hpp" #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/line_stepper.hpp" #include "detray/propagator/propagator.hpp" #include "detray/tracks/tracks.hpp" @@ -68,8 +67,7 @@ GTEST_TEST(detray_propagator, covariance_transport) { using navigator_t = navigator; using cline_stepper_t = line_stepper; - using actor_chain_t = actor_chain, - parameter_resetter>; + using actor_chain_t = actor_chain>; using propagator_t = propagator; @@ -91,8 +89,7 @@ GTEST_TEST(detray_propagator, covariance_transport) { geometry::barcode{}.set_index(0u), bound_vector, bound_cov); // Actors - parameter_transporter::state bound_updater{}; - parameter_resetter::state rst{}; + parameter_transporter::state bound_updater{bound_param0}; propagation::config prop_cfg{}; prop_cfg.navigation.overstep_tolerance = -100.f * unit::um; @@ -100,10 +97,10 @@ GTEST_TEST(detray_propagator, covariance_transport) { propagator_t::state propagation(bound_param0, det); // Run propagator - p.propagate(propagation, detray::tie(bound_updater, rst)); + p.propagate(propagation, detray::tie(bound_updater)); // Bound state after one turn propagation - const auto& bound_param1 = propagation._stepping.bound_params(); + const auto& bound_param1 = bound_updater.bound_params(); // Check if the track reaches the final surface EXPECT_EQ(bound_param0.surface_link().volume(), 4095u); @@ -122,12 +119,4 @@ GTEST_TEST(detray_propagator, covariance_transport) { matrix_operator().element(bound_cov1, i, j), tol); } } - - // Check covaraince - for (unsigned int i = 0u; i < e_bound_size; i++) { - for (unsigned int j = 0u; j < e_bound_size; j++) { - EXPECT_NEAR(matrix_operator().element(bound_cov0, i, j), - matrix_operator().element(bound_cov1, i, j), tol); - } - } } diff --git a/tutorials/src/device/cuda/propagation.hpp b/tutorials/src/device/cuda/propagation.hpp index ba8a02658..1f128b96e 100644 --- a/tutorials/src/device/cuda/propagation.hpp +++ b/tutorials/src/device/cuda/propagation.hpp @@ -15,8 +15,7 @@ #include "detray/navigation/navigator.hpp" #include "detray/propagator/actor_chain.hpp" #include "detray/propagator/actors/aborters.hpp" -#include "detray/propagator/actors/parameter_resetter.hpp" -#include "detray/propagator/actors/parameter_transporter.hpp" +#include "detray/propagator/actors/parameter_updater.hpp" #include "detray/propagator/actors/pointwise_material_interactor.hpp" #include "detray/propagator/propagator.hpp" #include "detray/propagator/rk_stepper.hpp" @@ -56,11 +55,14 @@ using stepper_t = rk_stepper; // Actors + +// Add the material interaction to the bound track parameter update scheme +using parameter_updater_t = parameter_updater< + detray::tutorial::algebra_t, + pointwise_material_interactor>; +// Make actor call chain using actor_chain_t = - actor_chain, - pointwise_material_interactor, - parameter_resetter>; + actor_chain; // Propagator using propagator_t = propagator; diff --git a/tutorials/src/device/cuda/propagation_kernel.cu b/tutorials/src/device/cuda/propagation_kernel.cu index 5d0e28b5d..b246ba8c6 100644 --- a/tutorials/src/device/cuda/propagation_kernel.cu +++ b/tutorials/src/device/cuda/propagation_kernel.cu @@ -47,11 +47,9 @@ __global__ void propagation_kernel( transporter_state{}; detray::pointwise_material_interactor::state interactor_state{}; - detray::parameter_resetter::state - resetter_state{}; - auto actor_states = detray::tie(aborter_state, transporter_state, - interactor_state, resetter_state); + auto actor_states = + detray::tie(aborter_state, transporter_state, interactor_state); // Create the propagator state for the track detray::tutorial::propagator_t::state state(tracks[gid], field_data, det);