diff --git a/core/include/detray/navigation/policies.hpp b/core/include/detray/navigation/policies.hpp index 4e1abc430..6a66ae7c6 100644 --- a/core/include/detray/navigation/policies.hpp +++ b/core/include/detray/navigation/policies.hpp @@ -73,8 +73,7 @@ struct stepper_default_policy : actor { auto &navigation = propagation._navigation; // Not a severe change to track state expected - // Policy is called after stepsize update -> use prev. step size - if (math::fabs(stepping.prev_step_size()) < + if (math::fabs(stepping.step_size()) < math::fabs( stepping.constraints().template size<>(stepping.direction())) - pol_state.tol) { @@ -111,20 +110,19 @@ struct stepper_rk_policy : actor { const auto &stepping = propagation._stepping; auto &navigation = propagation._navigation; - // Policy is called after stepsize update -> use prev. step size - const scalar rel_correction{(stepping.prev_step_size() - navigation()) / + // How strongly did the RKN algorithm reduce the step size? + const scalar rel_correction{(stepping.step_size() - navigation()) / navigation()}; // Large correction to the stepsize - re-initialize the volume if (rel_correction > pol_state.m_threshold_no_trust) { - // Re-evaluate only next candidate navigation.set_no_trust(); } - // Medium correction - re-evaluate the current candidates + // Medium correction - re-evaluate all current candidates else if (rel_correction > pol_state.m_threshold_fair_trust) { - // Re-evaluate all candidates navigation.set_fair_trust(); } else { + // Small correction - re-evaluate only the next candidate navigation.set_high_trust(); } } 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 81e9283ae..9266b0e7c 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,11 +31,13 @@ struct parameter_resetter : actor { return; } - // Surface + // Update free params after bound params were changed by actors const auto sf = navigation.get_surface(); + stepping() = sf.bound_to_free_vector(propagation._context, + stepping.bound_params()); - sf.template visit_mask(sf.transform(propagation._context), - 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 2e11f14b1..78d0181a5 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,52 @@ struct parameter_transporter : actor { return; } - using detector_type = typename propagator_state_t::detector_type; + // Geometry context for this track + const auto& gctx = propagation._context; // Current Surface const auto sf = navigation.get_surface(); + // Bound track params of departure surface + auto& bound_params = stepping.bound_params(); + // Covariance is transported only when the previous surface is an // actual tracking surface. (i.e. This disables the covariance transport // from curvilinear frame) - if (!detail::is_invalid_value(stepping.prev_sf_index())) { + if (!bound_params.surface_link().is_invalid()) { // Previous surface - tracking_surface 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(propagation._context, - stepping.bound_params()); + const bound_to_free_matrix_t bound_to_free_jacobian = + prev_sf.bound_to_free_jacobian(gctx, bound_params); + auto vol = navigation.get_volume(); + const auto vol_mat_ptr = + vol.has_material() ? vol.material_parameters(stepping().pos()) + : nullptr; stepping.set_full_jacobian( sf.template visit_mask( - sf.transform(propagation._context), bound_to_free_jacobian, - propagation)); + sf.transform(gctx), bound_to_free_jacobian, vol_mat_ptr, + propagation._stepping)); // Calculate surface-to-surface covariance transport const bound_matrix_t new_cov = - stepping.full_jacobian() * - stepping.bound_params().covariance() * + stepping.full_jacobian() * bound_params.covariance() * matrix_operator().transpose(stepping.full_jacobian()); + stepping.bound_params().set_covariance(new_cov); } // Convert free to bound vector - stepping.bound_params().set_parameter_vector( - sf.free_to_bound_vector(propagation._context, stepping())); + bound_params.set_parameter_vector( + sf.free_to_bound_vector(gctx, stepping())); // Set surface link - stepping.bound_params().set_surface_link(sf.barcode()); - - return; + bound_params.set_surface_link(sf.barcode()); } + }; // namespace detray } // namespace detray diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index 0c4e75f53..2f30f1a14 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 08b56300f..79998df00 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. @@ -88,12 +84,11 @@ class base_stepper { const typename detector_t::geometry_context &ctx) : m_bound_params(bound_params) { - // Surface + // Departure surface const auto sf = tracking_surface{det, bound_params.surface_link()}; - sf.template visit_mask< - typename parameter_resetter::kernel>( - sf.transform(ctx), sf.index(), *this); + // Set free track parameters for stepping/navigation + 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 4984e0a0d..58c0ed813 100644 --- a/core/include/detray/propagator/propagator.hpp +++ b/core/include/detray/propagator/propagator.hpp @@ -189,16 +189,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); @@ -298,16 +298,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(), @@ -401,7 +402,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 e3e18e981..bf37c7161 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -87,20 +87,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 @@ -118,14 +132,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 @@ -144,6 +160,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; }; @@ -156,9 +175,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 185e83dd9..dd66859b1 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);