diff --git a/core/include/detray/navigation/policies.hpp b/core/include/detray/navigation/policies.hpp index 52d45f526..4e1abc430 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.prev_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.prev_step_size() - navigation()) / navigation()}; // Large correction to the stepsize - re-initialize the volume diff --git a/core/include/detray/propagator/actors/aborters.hpp b/core/include/detray/propagator/actors/aborters.hpp index 49d9eca67..81fa2f18e 100644 --- a/core/include/detray/propagator/actors/aborters.hpp +++ b/core/include/detray/propagator/actors/aborters.hpp @@ -51,7 +51,7 @@ struct pathlimit_aborter : actor { const scalar step_limit = abrt_state.path_limit() - - math::fabs(prop_state._stepping._abs_path_length); + math::fabs(prop_state._stepping.abs_path_length()); // Check the path limit if (step_limit <= 0.f) { @@ -112,7 +112,7 @@ struct next_surface_aborter : actor { // Abort at the next sensitive surface if (navigation.is_on_sensitive() && - stepping._s > abrt_state.min_step_length) { + stepping.path_from_surface() > abrt_state.min_step_length) { prop_state._heartbeat &= navigation.abort(); abrt_state.success = true; } diff --git a/core/include/detray/propagator/actors/parameter_resetter.hpp b/core/include/detray/propagator/actors/parameter_resetter.hpp index 9aeeeb71e..7ab323941 100644 --- a/core/include/detray/propagator/actors/parameter_resetter.hpp +++ b/core/include/detray/propagator/actors/parameter_resetter.hpp @@ -42,16 +42,16 @@ struct parameter_resetter : actor { // Reset the free vector stepping() = detail::bound_to_free_vector(trf3, mask, - stepping._bound_params); + stepping.bound_params()); // Reset the path length - stepping._s = 0; + stepping.reset_path_from_surface(); // Reset jacobian transport to identity matrix - matrix_operator().set_identity(stepping._jac_transport); + stepping.reset_transport_jacobian(); // Reset the surface index - stepping._prev_sf_id = sf_idx; + stepping.set_prev_sf_index(sf_idx); } }; diff --git a/core/include/detray/propagator/actors/parameter_transporter.hpp b/core/include/detray/propagator/actors/parameter_transporter.hpp index 94e1c41d5..e9973c463 100644 --- a/core/include/detray/propagator/actors/parameter_transporter.hpp +++ b/core/include/detray/propagator/actors/parameter_transporter.hpp @@ -68,7 +68,7 @@ struct parameter_transporter : actor { // Transport jacobian in free coordinate const free_matrix_t& free_transport_jacobian = - stepping._jac_transport; + stepping.transport_jacobian(); // Path correction factor free_matrix_t path_correction = jacobian_engine_t::path_correction( @@ -107,32 +107,33 @@ struct parameter_transporter : actor { // Covariance is transported only when the previous surface is an // actual tracking surface. (i.e. This disables the covariance transport // from curvilinear frame) - if (stepping._prev_sf_id != detail::invalid_value()) { + if (!detail::is_invalid_value(stepping.prev_sf_index())) { // Previous surface tracking_surface prev_sf{navigation.detector(), - stepping._prev_sf_id}; + stepping.prev_sf_index()}; const bound_to_free_matrix bound_to_free_jacobian = - prev_sf.bound_to_free_jacobian(ctx, stepping._bound_params); + prev_sf.bound_to_free_jacobian(ctx, stepping.bound_params()); - stepping._full_jacobian = + stepping.set_full_jacobian( sf.template visit_mask( - sf.transform(ctx), bound_to_free_jacobian, propagation); + sf.transform(ctx), bound_to_free_jacobian, propagation)); // Calculate surface-to-surface covariance transport const bound_matrix_t new_cov = - stepping._full_jacobian * stepping._bound_params.covariance() * - matrix_operator().transpose(stepping._full_jacobian); - stepping._bound_params.set_covariance(new_cov); + stepping.full_jacobian() * + stepping.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( + stepping.bound_params().set_parameter_vector( sf.free_to_bound_vector(ctx, stepping())); // Set surface link - stepping._bound_params.set_surface_link(sf.barcode()); + stepping.bound_params().set_surface_link(sf.barcode()); return; } diff --git a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp index d1d973b08..58140c589 100644 --- a/core/include/detray/propagator/actors/pointwise_material_interactor.hpp +++ b/core/include/detray/propagator/actors/pointwise_material_interactor.hpp @@ -143,8 +143,8 @@ struct pointwise_material_interactor : actor { auto &stepping = prop_state._stepping; - this->update(geo_context_type{}, stepping._ptc, - stepping._bound_params, interactor_state, + this->update(geo_context_type{}, stepping.particle_hypothesis(), + stepping.bound_params(), interactor_state, static_cast(navigation.direction()), navigation.get_surface()); } diff --git a/core/include/detray/propagator/base_stepper.hpp b/core/include/detray/propagator/base_stepper.hpp index 3a2661fe5..15fde1349 100644 --- a/core/include/detray/propagator/base_stepper.hpp +++ b/core/include/detray/propagator/base_stepper.hpp @@ -43,20 +43,15 @@ template ; - using transform3_type = dtransform3D; - using matrix_operator = dmatrix_operator; - using free_track_parameters_type = free_track_parameters; using bound_track_parameters_type = bound_track_parameters; - using free_matrix_type = free_matrix; using bound_matrix_type = bound_matrix; - using bound_to_free_matrix_type = bound_to_free_matrix; + + using inspector_type = inspector_t; + using policy_type = policy_t; /// @brief State struct holding the track /// @@ -66,26 +61,23 @@ class base_stepper { /// Sets track parameters. DETRAY_HOST_DEVICE explicit state(const free_track_parameters_type &free_params) - : _track(free_params) { + : m_track(free_params) { curvilinear_frame cf(free_params); // Set bound track parameters - _bound_params.set_parameter_vector(cf.m_bound_vec); + m_bound_params.set_parameter_vector(cf.m_bound_vec); // A dummy covariance - should not be used - _bound_params.set_covariance( + m_bound_params.set_covariance( matrix_operator() .template identity()); // A dummy barcode - should not be used - _bound_params.set_surface_link(geometry::barcode{}); - - // Reset the path length - _s = 0.f; + m_bound_params.set_surface_link(geometry::barcode{}); // Reset jacobian transport to identity matrix - matrix_operator().set_identity(_jac_transport); + matrix_operator().set_identity(m_jac_transport); } /// Sets track parameters from bound track parameter. @@ -93,7 +85,7 @@ class base_stepper { DETRAY_HOST_DEVICE state( const bound_track_parameters_type &bound_params, const detector_t &det) - : _bound_params(bound_params) { + : m_bound_params(bound_params) { // Surface const auto sf = tracking_surface{det, bound_params.surface_link()}; @@ -104,135 +96,172 @@ class base_stepper { sf.transform(ctx), sf.index(), *this); } - /// free track parameter - free_track_parameters_type _track; - - /// Full jacobian - bound_matrix_type _full_jacobian = - matrix_operator().template identity(); + /// @returns free track parameters - non-const access + DETRAY_HOST_DEVICE + free_track_parameters_type &operator()() { return m_track; } - /// jacobian transport matrix - free_matrix_type _jac_transport = - matrix_operator().template identity(); + /// @returns free track parameters - const access + DETRAY_HOST_DEVICE + const free_track_parameters_type &operator()() const { return m_track; } - /// bound covariance - bound_track_parameters_type _bound_params; + /// @returns bound track parameters - const access + DETRAY_HOST_DEVICE + bound_track_parameters_type &bound_params() { return m_bound_params; } - /// @returns track parameters - const access + /// @returns bound track parameters. DETRAY_HOST_DEVICE - free_track_parameters_type &operator()() { return _track; } + const bound_track_parameters_type &bound_params() const { + return m_bound_params; + } - /// @returns track parameters. + /// Get stepping direction DETRAY_HOST_DEVICE - const free_track_parameters_type &operator()() const { return _track; } + inline step::direction direction() const { return m_direction; } - step::direction _direction{step::direction::e_forward}; + /// Set stepping direction + DETRAY_HOST_DEVICE inline void set_direction(step::direction dir) { + m_direction = dir; + } - // Stepping constraints - constraint_t _constraint = {}; + /// @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; + } - // Navigation policy state - typename policy_t::state _policy_state = {}; + /// Set next step size + DETRAY_HOST_DEVICE + inline void set_step_size(const scalar_type step) { + m_step_size = step; + } - /// The inspector type of the stepping - inspector_type _inspector; + /// @returns the current step size of this state. + DETRAY_HOST_DEVICE + inline scalar_type step_size() const { return m_step_size; } - /// Total number of step trials - std::size_t _n_total_trials{0u}; + /// Set previous step size + DETRAY_HOST_DEVICE + inline void set_prev_step_size(const scalar_type step) { + m_prev_step_size = step; + } - /// Track path length - scalar_type _path_length{0.f}; + /// @returns the previous step size of this state. + DETRAY_HOST_DEVICE + inline scalar_type prev_step_size() const { return m_prev_step_size; } - /// Absolute path length - scalar_type _abs_path_length{0.f}; + /// @returns this states path length. + DETRAY_HOST_DEVICE + inline scalar_type path_length() const { return m_path_length; } - /// Track path length from the last surface. It will be reset to 0 when - /// the track reaches a new surface - scalar_type _s{0.f}; + /// @returns this states total integration path length. + DETRAY_HOST_DEVICE + inline scalar_type abs_path_length() const { return m_abs_path_length; } - /// Current step size - scalar_type _step_size{0.f}; + /// @returns path length since last surface was encountered + DETRAY_HOST_DEVICE + inline scalar_type path_from_surface() const { return m_path_from_sf; } - /// Previous step size (DEBUG purpose only) - scalar_type _prev_step_size{0.f}; + /// Reset the path length since last surface + DETRAY_HOST_DEVICE inline void reset_path_from_surface() { + m_path_from_sf = 0.f; + } - /// The default particle hypothesis is muon - pdg_particle _ptc = muon(); + /// 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); + } - /// Previous barcode to calculate the bound_to_free_jacobian - dindex _prev_sf_id = detail::invalid_value(); + /// Updates the total number of step trials + DETRAY_HOST_DEVICE inline void count_trials() { ++m_n_total_trials; } - /// Volume material that track is passing through - const detray::material *_mat{nullptr}; + /// Set new step constraint + template + DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) { + m_constraint.template set(step_size); + } - /// Access the current volume material + /// @returns access to this states step constraints DETRAY_HOST_DEVICE - const auto &volume_material() const { - assert(_mat != nullptr); - return *_mat; - } + inline const constraint_t &constraints() const { return m_constraint; } - /// Set new step constraint + /// Remove [all] constraints template - DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) { - _constraint.template set(step_size); + DETRAY_HOST_DEVICE inline void release_step() { + m_constraint.template release(); } - /// Set new navigation direction - DETRAY_HOST_DEVICE inline void set_direction(step::direction dir) { - _direction = dir; + /// @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 access to this states step constraints + /// @returns true if the current volume contains volume material DETRAY_HOST_DEVICE - inline const constraint_t &constraints() const { return _constraint; } + bool volume_has_material() const { return (m_vol_mat != nullptr); } - /// @returns access to this states step constraints + /// Access the current volume material DETRAY_HOST_DEVICE - inline typename policy_t::state &policy_state() { - return _policy_state; + const auto &volume_material() const { + assert(m_vol_mat != nullptr); + return *m_vol_mat; } - /// @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 _n_total_trials; + /// Set new volume material access + DETRAY_HOST_DEVICE + inline void set_volume_material(const material *mat_ptr) { + m_vol_mat = mat_ptr; } - /// Updates the total number of step trials - DETRAY_HOST_DEVICE inline void count_trials() { ++_n_total_trials; } + /// Access the current particle hypothesis + DETRAY_HOST_DEVICE + const auto &particle_hypothesis() const { return m_ptc; } - /// @returns the navigation direction + /// Set new particle hypothesis DETRAY_HOST_DEVICE - inline step::direction direction() const { return _direction; } + inline void set_particle(const pdg_particle &ptc) { + m_ptc = ptc; + } - /// Remove [all] constraints - template - DETRAY_HOST_DEVICE inline void release_step() { - _constraint.template release(); + /// @returns the current transport Jacbian. + DETRAY_HOST_DEVICE + inline const free_matrix_type &transport_jacobian() const { + return m_jac_transport; } - /// Set next step size + /// @returns the current full Jacbian. DETRAY_HOST_DEVICE - inline void set_step_size(const scalar_type step) { _step_size = step; } + inline const bound_matrix_type &full_jacobian() const { + return m_full_jacobian; + } - /// @returns the current step size of this state. + /// Set new full Jacbian. DETRAY_HOST_DEVICE - inline scalar_type step_size() const { return _step_size; } + inline void set_full_jacobian(const bound_matrix_type &jac) { + m_full_jacobian = jac; + } - /// @returns this states remaining path length. + /// Reset transport Jacbian. DETRAY_HOST_DEVICE - inline scalar_type path_length() const { return _path_length; } + inline void reset_transport_jacobian() { + matrix_operator().set_identity(m_jac_transport); + } - /// @returns the current transport Jacbian. + /// @returns access to this states step constraints DETRAY_HOST_DEVICE - inline const free_matrix_type &transport_jacobian() const { - return _jac_transport; + inline typename policy_t::state &policy_state() { + return m_policy_state; } /// @returns the stepping inspector DETRAY_HOST - constexpr auto &inspector() { return _inspector; } + constexpr auto &inspector() { return m_inspector; } /// Call the stepping inspector DETRAY_HOST_DEVICE @@ -240,9 +269,71 @@ class base_stepper { [[maybe_unused]] const char *message) { if constexpr (!std::is_same_v) { - _inspector(*this, cfg, message); + m_inspector(*this, cfg, message); } } + + protected: + /// Set new transport Jacbian. + DETRAY_HOST_DEVICE + inline void set_transport_jacobian(const free_matrix_type &jac) { + m_jac_transport = jac; + } + + private: + /// Stepping direction + step::direction m_direction{step::direction::e_forward}; + + /// Total number of step trials + std::size_t m_n_total_trials{0u}; + + /// Current step size + scalar_type m_step_size{0.f}; + + /// Previous step size + scalar_type m_prev_step_size{0.f}; + + /// Track path length + scalar_type m_path_length{0.f}; + + /// Absolute path length + 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 = {}; + + /// Navigation policy state + typename policy_t::state m_policy_state = {}; + + /// The inspector type of the stepping + inspector_type m_inspector; }; }; diff --git a/core/include/detray/propagator/line_stepper.hpp b/core/include/detray/propagator/line_stepper.hpp index 59864657d..3d54c978c 100644 --- a/core/include/detray/propagator/line_stepper.hpp +++ b/core/include/detray/propagator/line_stepper.hpp @@ -1,6 +1,6 @@ /** Detray library, part of the ACTS project (R&D line) * - * (c) 2021-2022 CERN for the benefit of the ACTS project + * (c) 2021-2024 CERN for the benefit of the ACTS project * * Mozilla Public License Version 2.0 */ @@ -22,10 +22,10 @@ template { - public: using base_type = base_stepper; + public: using algebra_type = algebra_t; using scalar_type = dscalar; using vector3_type = dvector3D; @@ -49,11 +49,10 @@ class line_stepper final /// Update the track state in a straight line. DETRAY_HOST_DEVICE inline void advance_track() { - auto& track = this->_track; - track.set_pos(track.pos() + track.dir() * this->_step_size); + auto& track = (*this)(); + track.set_pos(track.pos() + track.dir() * this->step_size()); - this->_path_length += this->_step_size; - this->_s += this->_step_size; + this->update_path_lengths(this->step_size()); } DETRAY_HOST_DEVICE @@ -65,14 +64,14 @@ class line_stepper final // d(x,y,z)/d(n_x,n_y,n_z) matrix_type<3, 3> dxdn = - this->_step_size * matrix_operator().template identity<3, 3>(); + this->step_size() * matrix_operator().template identity<3, 3>(); matrix_operator().template set_block<3, 3>(D, dxdn, e_free_pos0, e_free_dir0); /// NOTE: Let's skip the element for d(time)/d(qoverp) for the /// moment.. - this->_jac_transport = D * this->_jac_transport; + this->set_transport_jacobian(D * this->transport_jacobian()); } DETRAY_HOST_DEVICE @@ -94,10 +93,10 @@ class line_stepper final const bool = true) const { // Straight line stepping: The distance given by the navigator is exact - stepping._step_size = dist_to_next; + stepping.set_step_size(dist_to_next); // Update navigation direction - const step::direction step_dir = stepping._step_size >= 0.f + const step::direction step_dir = stepping.step_size() >= 0.f ? step::direction::e_forward : step::direction::e_backward; stepping.set_direction(step_dir); @@ -105,7 +104,7 @@ class line_stepper final // Check constraints if (const scalar_type max_step = stepping.constraints().template size<>(stepping.direction()); - math::fabs(stepping._step_size) > math::fabs(max_step)) { + math::fabs(stepping.step_size()) > math::fabs(max_step)) { // Run inspection before step size is cut stepping.run_inspector(cfg, "Before constraint: "); diff --git a/core/include/detray/propagator/propagator.hpp b/core/include/detray/propagator/propagator.hpp index 75d6990ff..4151887f6 100644 --- a/core/include/detray/propagator/propagator.hpp +++ b/core/include/detray/propagator/propagator.hpp @@ -107,7 +107,7 @@ struct propagator { /// Set the particle hypothesis DETRAY_HOST_DEVICE void set_particle(const pdg_particle &ptc) { - _stepping._ptc = ptc; + _stepping.set_particle(ptc); } // Is the propagation still alive? @@ -170,8 +170,9 @@ struct propagator { // Set access to the volume material for the stepper auto vol = navigation.get_volume(); - stepping._mat = - vol.has_material() ? vol.material_parameters(track.pos()) : nullptr; + stepping.set_volume_material(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 @@ -274,9 +275,9 @@ struct propagator { // Set access to the volume material for the stepper auto vol = navigation.get_volume(); - stepping._mat = vol.has_material() - ? vol.material_parameters(track.pos()) - : nullptr; + stepping.set_volume_material( + vol.has_material() ? vol.material_parameters(track.pos()) + : nullptr); // Break automatic step size scaling by the stepper const bool reset_stepsize{navigation.is_on_surface() || @@ -374,7 +375,7 @@ struct propagator { } propagation.debug_stream << "step_size: " << std::setw(10) - << stepping._prev_step_size << std::endl; + << stepping.prev_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 46f2e4886..eab488a34 100644 --- a/core/include/detray/propagator/rk_stepper.hpp +++ b/core/include/detray/propagator/rk_stepper.hpp @@ -27,22 +27,20 @@ namespace detray { template class array_t = darray> + typename inspector_t = stepping::void_inspector> class rk_stepper final : public base_stepper { - public: using base_type = base_stepper; + public: using algebra_type = algebra_t; using scalar_type = dscalar; using point3_type = dpoint3D; using vector3_type = dvector3D; using transform3_type = dtransform3D; using matrix_operator = dmatrix_operator; - using mat_helper = matrix_helper; using free_track_parameters_type = typename base_type::free_track_parameters_type; using bound_track_parameters_type = @@ -55,36 +53,27 @@ class rk_stepper final struct state : public base_type::state { + friend rk_stepper; + static constexpr const stepping::id id = stepping::id::e_rk; DETRAY_HOST_DEVICE state(const free_track_parameters_type& t, const magnetic_field_t& mag_field) - : base_type::state(t), _magnetic_field(mag_field) {} + : base_type::state(t), m_magnetic_field(mag_field) {} template DETRAY_HOST_DEVICE state( const bound_track_parameters_type& bound_params, const magnetic_field_t& mag_field, const detector_t& det) - : base_type::state(bound_params, det), _magnetic_field(mag_field) {} + : base_type::state(bound_params, det), + m_magnetic_field(mag_field) {} - /// stepping data required for RKN4 - struct { - vector3_type b_first{0.f, 0.f, 0.f}; - vector3_type b_middle{0.f, 0.f, 0.f}; - vector3_type b_last{0.f, 0.f, 0.f}; - // t = tangential direction = dr/ds - std::array t; - // q/p - std::array qop; - // dt/ds = d^2r/ds^2 = q/p ( t X B ) - std::array dtds; - // d(q/p)/ds - std::array dqopds; - } _step_data; + /// @returns the B-field view + magnetic_field_type field() const { return m_magnetic_field; } - /// Magnetic field view - const magnetic_field_t _magnetic_field; + /// @returns access to the step data + const auto& step_data() const { return m_step_data; } /// Update the track state by Runge-Kutta-Nystrom integration. DETRAY_HOST_DEVICE @@ -133,10 +122,29 @@ class rk_stepper final [[maybe_unused]] Args&&... args) { if constexpr (!std::is_same_v) { - this->_inspector(*this, cfg, message, - std::forward(args)...); + this->inspector()(*this, cfg, message, + std::forward(args)...); } } + + private: + /// stepping data required for RKN4 + struct { + vector3_type b_first{0.f, 0.f, 0.f}; + vector3_type b_middle{0.f, 0.f, 0.f}; + vector3_type b_last{0.f, 0.f, 0.f}; + // t = tangential direction = dr/ds + std::array t; + // q/p + std::array qop; + // dt/ds = d^2r/ds^2 = q/p ( t X B ) + std::array dtds; + // d(q/p)/ds + std::array dqopds; + } m_step_data; + + /// Magnetic field view + const magnetic_field_t m_magnetic_field; }; /// Take a step, using an adaptive Runge-Kutta algorithm. diff --git a/core/include/detray/propagator/rk_stepper.ipp b/core/include/detray/propagator/rk_stepper.ipp index 8ff99fade..8a8e48b5a 100644 --- a/core/include/detray/propagator/rk_stepper.ipp +++ b/core/include/detray/propagator/rk_stepper.ipp @@ -9,16 +9,15 @@ #include "detray/geometry/tracking_volume.hpp" template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline void detray::rk_stepper::state::advance_track() { + inspector_t>::state::advance_track() { - const auto& sd = this->_step_data; - const scalar_type h{this->_step_size}; + const auto& sd = m_step_data; + const scalar_type h{this->step_size()}; const scalar_type h_6{h * static_cast(1. / 6.)}; - auto& track = this->_track; + auto& track = (*this)(); auto pos = track.pos(); auto dir = track.dir(); @@ -34,7 +33,7 @@ detray::rk_stepper_mat == nullptr)) { + if (this->volume_has_material()) { // 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]) + @@ -42,18 +41,15 @@ detray::rk_stepper_path_length += h; - this->_abs_path_length += math::fabs(h); - this->_s += h; + // Update path lengths + this->update_path_lengths(h); } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline void detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::advance_jacobian(const detray::stepping::config& cfg) { + magnetic_field_t, algebra_t, constraint_t, policy_t, + inspector_t>::state::advance_jacobian(const detray::stepping::config& cfg) { /// 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 @@ -77,10 +73,9 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< //( JacTransport = D * JacTransport ) auto D = matrix_operator().template identity(); - const auto& sd = this->_step_data; - const scalar_type h{this->_step_size}; - // const auto& mass = this->_mass; - auto& track = this->_track; + const auto& sd = m_step_data; + const scalar_type h{this->step_size()}; + auto& track = (*this)(); // Half step length const scalar_type h2{h * h}; @@ -217,6 +212,7 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< /*----------------------------------------------------------------- * Calculate the first terms of dk_n/dt1 -------------------------------------------------------------------*/ + using mat_helper = matrix_helper; // dk1/dt1 dkndt[0u] = sd.qop[0u] * mat_helper().column_wise_cross(dkndt[0u], sd.b_first); @@ -327,24 +323,24 @@ DETRAY_HOST_DEVICE inline void detray::rk_stepper< matrix_operator().set_block(D, dFdqop, 0u, 7u); matrix_operator().set_block(D, dGdqop, 4u, 7u); - this->_jac_transport = D * this->_jac_transport; + this->set_transport_jacobian(D * this->transport_jacobian()); } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::evaluate_dqopds(const std::size_t i, const scalar_type h, - const scalar_type dqopds_prev, - const detray::stepping::config& cfg) + 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) -> scalar_type { - const auto& track = this->_track; + const auto& track = (*this)(); const scalar_type qop = track.qop(); - auto& sd = this->_step_data; + auto& sd = m_step_data; - if (this->_mat == nullptr) { + if (!(this->volume_has_material())) { sd.qop[i] = qop; return 0.f; } else { @@ -368,17 +364,16 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::evaluate_dtds(const vector3_type& b_field, - const std::size_t i, const scalar_type h, - const vector3_type& dtds_prev, - const scalar_type qop) -> vector3_type { - auto& track = this->_track; + magnetic_field_t, algebra_t, constraint_t, policy_t, + inspector_t>::state::evaluate_dtds(const vector3_type& b_field, + const std::size_t i, const scalar_type h, + const vector3_type& dtds_prev, + const scalar_type qop) -> vector3_type { + auto& track = (*this)(); const auto dir = track.dir(); - auto& sd = this->_step_data; + auto& sd = m_step_data; if (i == 0u) { sd.t[i] = dir; @@ -392,11 +387,10 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::evaluate_field_gradient(const point3_type& pos) + magnetic_field_t, algebra_t, constraint_t, policy_t, + inspector_t>::state::evaluate_field_gradient(const point3_type& pos) -> matrix_type<3, 3> { matrix_type<3, 3> dBdr = matrix_operator().template zero<3, 3>(); @@ -408,7 +402,7 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< point3_type dpos1 = pos; dpos1[i] += delta; const auto bvec1_tmp = - this->_magnetic_field.at(dpos1[0], dpos1[1], dpos1[2]); + this->m_magnetic_field.at(dpos1[0], dpos1[1], dpos1[2]); vector3_type bvec1; bvec1[0u] = bvec1_tmp[0u]; bvec1[1u] = bvec1_tmp[1u]; @@ -417,7 +411,7 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< point3_type dpos2 = pos; dpos2[i] -= delta; const auto bvec2_tmp = - this->_magnetic_field.at(dpos2[0], dpos2[1], dpos2[2]); + this->m_magnetic_field.at(dpos2[0], dpos2[1], dpos2[2]); vector3_type bvec2; bvec2[0u] = bvec2_tmp[0u]; bvec2[1u] = bvec2_tmp[1u]; @@ -434,65 +428,63 @@ DETRAY_HOST_DEVICE inline auto detray::rk_stepper< } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline auto detray::rk_stepper::state::dtds() const -> vector3_type { + inspector_t>::state::dtds() const -> vector3_type { // In case there was no step before - if (this->_path_length == 0.f) { - const point3_type pos = this->_track.pos(); + if (this->path_length() == 0.f) { + const point3_type pos = (*this)().pos(); - const auto bvec_tmp = this->_magnetic_field.at(pos[0], pos[1], pos[2]); + const auto bvec_tmp = this->m_magnetic_field.at(pos[0], pos[1], pos[2]); vector3_type bvec; bvec[0u] = bvec_tmp[0u]; bvec[1u] = bvec_tmp[1u]; bvec[2u] = bvec_tmp[2u]; - return this->_track.qop() * vector::cross(this->_track.dir(), bvec); + return (*this)().qop() * vector::cross((*this)().dir(), bvec); } - return this->_step_data.dtds[3u]; + return m_step_data.dtds[3u]; } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline auto detray::rk_stepper::state::dqopds() const -> scalar_type { + inspector_t>::state::dqopds() const -> scalar_type { // In case there was no step before - if (this->_path_length == 0.f) { - return this->dqopds(this->_track.qop()); + if (this->path_length() == 0.f) { + return this->dqopds((*this)().qop()); } - return this->_step_data.dqopds[3u]; + return m_step_data.dqopds[3u]; } template class array_t> -DETRAY_HOST_DEVICE auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::dqopds(const scalar_type qop) const -> scalar_type { + typename policy_t, typename inspector_t> +DETRAY_HOST_DEVICE auto +detray::rk_stepper::state::dqopds(const scalar_type qop) const + -> scalar_type { // d(qop)ds is zero for empty space - if (this->_mat == nullptr) { + if (!(this->volume_has_material())) { return 0.f; } const auto& mat = this->volume_material(); - const scalar_type q = this->_ptc.charge(); + const scalar_type q = this->particle_hypothesis().charge(); const scalar_type p = q / qop; - const scalar_type mass = this->_ptc.mass(); + const scalar_type mass = this->particle_hypothesis().mass(); const scalar_type E = math::sqrt(p * p + mass * mass); // Compute stopping power const scalar_type stopping_power = - interaction().compute_stopping_power(mat, this->_ptc, - {mass, qop, q}); + interaction().compute_stopping_power( + mat, this->particle_hypothesis(), {mass, qop, q}); // Assert that a momentum is a positive value assert(p >= 0.f); @@ -503,22 +495,22 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< } template class array_t> -DETRAY_HOST_DEVICE auto detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, - array_t>::state::d2qopdsdqop(const scalar_type qop) const -> scalar_type { + typename policy_t, typename inspector_t> +DETRAY_HOST_DEVICE auto +detray::rk_stepper::state::d2qopdsdqop(const scalar_type qop) const + -> scalar_type { - if (this->_mat == nullptr) { + if (!(this->volume_has_material())) { return 0.f; } const auto& mat = this->volume_material(); - const scalar_type q = this->_ptc.charge(); + const scalar_type q = this->particle_hypothesis().charge(); const scalar_type p = q / qop; const scalar_type p2 = p * p; - const auto& mass = this->_ptc.mass(); + const auto& mass = this->particle_hypothesis().mass(); const scalar_type E2 = p2 + mass * mass; // Interaction object @@ -526,11 +518,12 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< // g = dE/ds = -1 * (-dE/ds) = -1 * stopping power const detail::relativistic_quantities rq(mass, qop, q); - const scalar_type g = -1.f * I.compute_stopping_power(mat, this->_ptc, rq); + const scalar_type g = + -1.f * I.compute_stopping_power(mat, this->particle_hypothesis(), rq); // dg/d(qop) = -1 * derivation of stopping power const scalar_type dgdqop = - -1.f * I.derive_stopping_power(mat, this->_ptc, rq); + -1.f * I.derive_stopping_power(mat, this->particle_hypothesis(), rq); // d(qop)/ds = - qop^3 * E * g / q^2 const scalar_type dqopds = this->dqopds(qop); @@ -541,29 +534,28 @@ DETRAY_HOST_DEVICE auto detray::rk_stepper< } template class array_t> + typename policy_t, typename inspector_t> DETRAY_HOST_DEVICE inline bool detray::rk_stepper< - magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t, array_t>:: + magnetic_field_t, algebra_t, constraint_t, policy_t, inspector_t>:: step(const scalar_type dist_to_next, detray::rk_stepper::state& stepping, + inspector_t>::state& stepping, const detray::stepping::config& cfg, const bool do_reset) const { // Get stepper and navigator states - auto& magnetic_field = stepping._magnetic_field; + auto& magnetic_field = stepping.m_magnetic_field; if (do_reset) { - stepping._step_size = dist_to_next; - } else if (stepping._step_size > 0) { - stepping._step_size = math::min(stepping._step_size, dist_to_next); + 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)); } else { - stepping._step_size = math::max(stepping._step_size, dist_to_next); + stepping.set_step_size(math::max(stepping.step_size(), dist_to_next)); } const point3_type pos = stepping().pos(); - auto& sd = stepping._step_data; + auto& sd = stepping.m_step_data; scalar_type error_estimate{0.f}; @@ -650,7 +642,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< for (unsigned int i = 0u; i < cfg.max_rk_updates; i++) { stepping.count_trials(); - error = math::max(estimate_error(stepping._step_size), + error = math::max(estimate_error(stepping.step_size()), static_cast(1e-20)); // Error is small enough @@ -661,7 +653,8 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< // Error estimate is too big // ---> Make step size smaller and estimate error again else { - stepping._step_size *= step_size_scaling(error); + stepping.set_step_size(stepping.step_size() * + step_size_scaling(error)); // Run inspection while the stepsize is getting adjusted stepping.run_inspector(cfg, "Adjust stepsize: ", i, @@ -670,7 +663,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< } // Update navigation direction - const step::direction step_dir = stepping._step_size >= 0.f + const step::direction step_dir = stepping.step_size() >= 0.f ? step::direction::e_forward : step::direction::e_backward; stepping.set_direction(step_dir); @@ -678,7 +671,7 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< // Check constraints if (const scalar_type max_step = stepping.constraints().template size<>(stepping.direction()); - math::fabs(stepping._step_size) > math::fabs(max_step)) { + math::fabs(stepping.step_size()) > math::fabs(max_step)) { // Run inspection before step size is cut stepping.run_inspector(cfg, "Before constraint: "); @@ -695,14 +688,14 @@ DETRAY_HOST_DEVICE inline bool detray::rk_stepper< } // Save the current step size - stepping._prev_step_size = stepping._step_size; + stepping.set_prev_step_size(stepping.step_size()); // Update the step size for the next step - stepping._step_size *= step_size_scaling(error); + stepping.set_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._step_size = cfg.min_stepsize; + if (stepping.step_size() < cfg.min_stepsize) { + stepping.set_step_size(cfg.min_stepsize); } // Run final inspection diff --git a/tests/include/detray/test/utils/inspectors.hpp b/tests/include/detray/test/utils/inspectors.hpp index 8e4890e5f..d07a0d476 100644 --- a/tests/include/detray/test/utils/inspectors.hpp +++ b/tests/include/detray/test/utils/inspectors.hpp @@ -445,13 +445,13 @@ struct print_inspector { << "\t\t" << step_scalor << std::endl; debug_stream << "Bfield points:" << std::endl; - const auto &f = state._step_data.b_first; + const auto &f = state.step_data().b_first; debug_stream << "\tfirst:" << tabs << f[0] << ", " << f[1] << ", " << f[2] << std::endl; - const auto &m = state._step_data.b_middle; + const auto &m = state.step_data().b_middle; debug_stream << "\tmiddle:" << tabs << m[0] << ", " << m[1] << ", " << m[2] << std::endl; - const auto &l = state._step_data.b_last; + const auto &l = state.step_data().b_last; debug_stream << "\tlast:" << tabs << l[0] << ", " << l[1] << ", " << l[2] << std::endl; diff --git a/tests/include/detray/test/utils/simulation/random_scatterer.hpp b/tests/include/detray/test/utils/simulation/random_scatterer.hpp index c1535c13e..bb272430a 100644 --- a/tests/include/detray/test/utils/simulation/random_scatterer.hpp +++ b/tests/include/detray/test/utils/simulation/random_scatterer.hpp @@ -144,8 +144,8 @@ struct random_scatterer : actor { } auto& stepping = prop_state._stepping; - const auto& ptc = stepping._ptc; - auto& bound_params = stepping._bound_params; + const auto& ptc = stepping.particle_hypothesis(); + auto& bound_params = stepping.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 +170,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)); + stepping.bound_params().set_phi(getter::phi(new_dir)); + stepping.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 1e38fd0c3..7934d2a64 100644 --- a/tests/include/detray/test/validation/material_validation_utils.hpp +++ b/tests/include/detray/test/validation/material_validation_utils.hpp @@ -191,7 +191,7 @@ struct material_tracer : detray::actor { typename parameter_transporter::state &, typename propagator_state_t::actor_chain_type:: state>) { - const auto &track_param = prop_state._stepping._bound_params; + const auto &track_param = prop_state._stepping.bound_params(); loc_pos = track_param.bound_local(); } else { const auto &track_param = prop_state._stepping(); diff --git a/tests/integration_tests/cpu/material/material_interaction.cpp b/tests/integration_tests/cpu/material/material_interaction.cpp index 04f89000b..b9ef2823e 100644 --- a/tests/integration_tests/cpu/material/material_interaction.cpp +++ b/tests/integration_tests/cpu/material/material_interaction.cpp @@ -113,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{state._stepping.bound_params().p(ptc.charge())}; // mass const auto mass = ptc.mass(); @@ -126,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(state._stepping.bound_params().covariance(), e_bound_qoverp, e_bound_qoverp)}; // Interaction object @@ -199,13 +199,13 @@ GTEST_TEST(detray_material, telescope_geometry_energy_loss) { // Propagate alt_p.propagate(alt_state, alt_actor_states); - alt_bound_param = alt_state._stepping._bound_params; + 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) { // if (alt_state._navigation.is_complete()){ const scalar altP = - alt_state._stepping._bound_params.p(ptc.charge()); + alt_state._stepping.bound_params().p(ptc.charge()); altE = std::hypot(altP, mass); break; } @@ -293,7 +293,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 = state._stepping.bound_params(); // Updated phi and theta variance if (i == 0u) { @@ -395,7 +395,7 @@ GTEST_TEST(detray_material, telescope_geometry_volume_material) { const auto eloss_approx = interaction().compute_energy_loss_bethe_bloch( - state._stepping._path_length, mat, ptc, + state._stepping.path_length(), mat, ptc, {ptc, bound_param.qop()}); const auto iniE = std::sqrt(iniP * iniP + mass * mass); diff --git a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp index 851eb9aae..28f6c75e0 100644 --- a/tests/integration_tests/cpu/propagator/jacobian_validation.cpp +++ b/tests/integration_tests/cpu/propagator/jacobian_validation.cpp @@ -407,7 +407,7 @@ 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) / + stepping.prev_step_size()) / N; // Warning for too many step counts @@ -428,23 +428,23 @@ 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 = stepping.bound_params(); } // Get the bound track parameters and jacobian at the destination // surface else if ((navigation.is_on_sensitive() || navigation.is_on_passive()) && navigation.barcode().index() == 1u) { - 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_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(); // Stop navigation if the destination surface found propagation._heartbeat &= navigation.exit(); } - if (stepping._path_length > actor_state.m_min_path_length) { + if (stepping.path_length() > actor_state.m_min_path_length) { propagation._navigation.set_no_trust(); } diff --git a/tests/integration_tests/cpu/propagator/propagator.cpp b/tests/integration_tests/cpu/propagator/propagator.cpp index 79635635b..316790076 100644 --- a/tests/integration_tests/cpu/propagator/propagator.cpp +++ b/tests/integration_tests/cpu/propagator/propagator.cpp @@ -74,44 +74,45 @@ struct helix_inspector : actor { inspector_state._nav_status.push_back(navigation.status()); // Nothing has happened yet (first call of actor chain) - if (stepping.path_length() < tol || stepping._s < tol) { + if (stepping.path_length() < tol || + stepping.path_from_surface() < tol) { return; } - if (stepping._bound_params.surface_link().is_invalid()) { + if (stepping.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(), stepping.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, stepping.bound_params()); const auto last_pos = free_params.pos(); const auto bvec = - stepping._magnetic_field.at(last_pos[0], last_pos[1], last_pos[2]); + stepping.field().at(last_pos[0], last_pos[1], last_pos[2]); const vector3 b{bvec[0], bvec[1], bvec[2]}; detail::helix hlx(free_params, &b); - const auto true_pos = hlx(stepping._s); + const auto true_pos = hlx(stepping.path_from_surface()); - const point3 relative_error{1.f / stepping._s * + const point3 relative_error{1.f / stepping.path_from_surface() * (stepping().pos() - true_pos)}; ASSERT_NEAR(getter::norm(relative_error), 0.f, tol); - auto true_J = hlx.jacobian(stepping._s); + auto true_J = hlx.jacobian(stepping.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._jac_transport, i, j), - matrix_operator().element(true_J, i, j), - stepping._s * tol * 10.f); + ASSERT_NEAR(matrix_operator().element( + stepping.transport_jacobian(), i, j), + matrix_operator().element(true_J, i, j), + stepping.path_from_surface() * tol * 10.f); } } } @@ -267,8 +268,8 @@ TEST_P(PropagatorWithRkStepper, rk4_propagator_const_bfield) { //<< lim_state.debug_stream.str() << std::endl; << lim_state._navigation.inspector().to_string() << std::endl; - ASSERT_GE(std::abs(path_limit), lim_state._stepping._abs_path_length) - << "Absolute path length: " << lim_state._stepping._abs_path_length + ASSERT_GE(std::abs(path_limit), lim_state._stepping.abs_path_length()) + << "Absolute path length: " << lim_state._stepping.abs_path_length() << ", path limit: " << path_limit << std::endl; //<< state._navigation.inspector().to_string() << std::endl; diff --git a/tests/unit_tests/cpu/detectors/telescope_detector.cpp b/tests/unit_tests/cpu/detectors/telescope_detector.cpp index e0822b919..819ebc8cb 100644 --- a/tests/unit_tests/cpu/detectors/telescope_detector.cpp +++ b/tests/unit_tests/cpu/detectors/telescope_detector.cpp @@ -233,12 +233,12 @@ GTEST_TEST(detray_detectors, telescope_detector) { navigator_x.update(stepping_x(), navigation_x, prop_cfg.navigation); // The track path lengths should match between all propagations EXPECT_NEAR( - std::abs(stepping_z1._path_length - stepping_z2._path_length) / - stepping_z1._path_length, + std::abs(stepping_z1.path_length() - stepping_z2.path_length()) / + stepping_z1.path_length(), 0.f, tol); EXPECT_NEAR( - std::abs(stepping_z1._path_length - stepping_x._path_length) / - stepping_x._path_length, + std::abs(stepping_z1.path_length() - stepping_x.path_length()) / + stepping_x.path_length(), 0.f, tol); // The track positions in z should match exactly EXPECT_NEAR(getter::norm(stepping_z1().pos() - stepping_z2().pos()) / diff --git a/tests/unit_tests/cpu/propagator/covariance_transport.cpp b/tests/unit_tests/cpu/propagator/covariance_transport.cpp index ff174bd2e..f1322caeb 100644 --- a/tests/unit_tests/cpu/propagator/covariance_transport.cpp +++ b/tests/unit_tests/cpu/propagator/covariance_transport.cpp @@ -102,7 +102,7 @@ GTEST_TEST(detray_propagator, covariance_transport) { p.propagate(propagation, detray::tie(bound_updater, rst)); // Bound state after one turn propagation - const auto& bound_param1 = propagation._stepping._bound_params; + const auto& bound_param1 = propagation._stepping.bound_params(); // const auto bound_vec0 = bound_param0; // const auto bound_vec1 = bound_param1; diff --git a/tests/unit_tests/cpu/propagator/rk_stepper.cpp b/tests/unit_tests/cpu/propagator/rk_stepper.cpp index d07d4ab4a..203e1ca37 100644 --- a/tests/unit_tests/cpu/propagator/rk_stepper.cpp +++ b/tests/unit_tests/cpu/propagator/rk_stepper.cpp @@ -248,7 +248,7 @@ TEST(detray_propagator, qop_derivative) { // RK Stepping into forward direction rk_stepper_t::state rk_state{track, hom_bfield}; - rk_state._mat = &vol_mat; + rk_state.set_volume_material(&vol_mat); for (unsigned int i_s = 0u; i_s < rk_steps; i_s++) {