detray/propagator/base_stepper.hpp

0001 // This file is part of the ACTS project.
0002 //
0003 // Copyright (C) 2016 CERN for the benefit of the ACTS project
0004 //
0005 // This Source Code Form is subject to the terms of the Mozilla Public
0006 // License, v. 2.0. If a copy of the MPL was not distributed with this
0007 // file, You can obtain one at https://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 // Project include(s).
0012 #include "detray/definitions/algebra.hpp"
0013 #include "detray/definitions/detail/qualifiers.hpp"
0014 #include "detray/definitions/pdg_particle.hpp"
0015 #include "detray/definitions/units.hpp"
0016 #include "detray/geometry/identifier.hpp"
0017 #include "detray/geometry/tracking_surface.hpp"
0018 #include "detray/propagator/constrained_step.hpp"
0019 #include "detray/propagator/detail/print_stepper_state.hpp"
0020 #include "detray/propagator/stepping_config.hpp"
0021 #include "detray/propagator/transport_jacobian.hpp"
0022 #include "detray/tracks/tracks.hpp"
0023 #include "detray/utils/curvilinear_frame.hpp"
0024 #include "detray/utils/logging.hpp"
0025
0026 namespace detray {
0027
0028 namespace stepping {
0029
0030 /// A void inpector that does nothing.
0031 ///
0032 /// Inspectors can be plugged in to understand the current stepper state.
0033 struct void_inspector {
0034   template <typename state_t>
0035   DETRAY_HOST_DEVICE constexpr void operator()(const state_t & /*ignored*/,
0036                                                const char * /*ignored*/) const {
0037     /*Do nothing*/
0038   }
0039 };
0040
0041 }  // namespace stepping
0042
0043 /// Base stepper implementation
0044 template <concepts::algebra algebra_t, typename constraint_t, typename policy_t,
0045           typename inspector_t = stepping::void_inspector>
0046 class base_stepper {
0047  public:
0048   using scalar_type = dscalar<algebra_t>;
0049
0050   using free_track_parameters_type = free_track_parameters<algebra_t>;
0051   using bound_track_parameters_type = bound_track_parameters<algebra_t>;
0052   using free_matrix_type = free_matrix<algebra_t>;
0053   using bound_matrix_type = bound_matrix<algebra_t>;
0054   using magnetic_field_type = void;
0055
0056   using inspector_type = inspector_t;
0057   using policy_type = policy_t;
0058
0059   /// @brief State struct holding the track
0060   ///
0061   /// @note It has to cast into a const track via the call operation.
0062   template <typename internal_jacobian_t = free_matrix_type>
0063   struct state {
0064     using internal_jacobian_type = internal_jacobian_t;
0065
0066     /// Sets track parameters.
0067     DETRAY_HOST_DEVICE
0068     explicit state(const free_track_parameters_type &free_params)
0069         : m_track(free_params) {
0070       assert(!m_track.is_invalid());
0071
0072       // HACK: When the overload resolution for the transport Jacobian
0073       // type is resolved, turn this into a default member
0074       // initialization.
0075       if constexpr (concepts::transport_jacobian<internal_jacobian_type>) {
0076         m_jac_transport = internal_jacobian_type::identity();
0077       } else {
0078         m_jac_transport = matrix::identity<internal_jacobian_type>();
0079       }
0080     }
0081
0082     /// Sets track parameters from bound track parameter.
0083     template <typename detector_t>
0084     DETRAY_HOST_DEVICE state(const bound_track_parameters_type &bound_params,
0085                              const detector_t &det,
0086                              const typename detector_t::geometry_context &ctx) {
0087       assert(!bound_params.is_invalid());
0088       assert(!bound_params.surface_link().is_invalid());
0089
0090       // Departure surface
0091       const auto sf = tracking_surface{det, bound_params.surface_link()};
0092
0093       // Set free track parameters for stepping/navigation
0094       m_track = sf.bound_to_free_vector(ctx, bound_params);
0095
0096       assert(!m_track.is_invalid());
0097
0098       // HACK: When the overload resolution for the transport Jacobian
0099       // type is resolved, turn this into a default member
0100       // initialization.
0101       if constexpr (concepts::transport_jacobian<internal_jacobian_type>) {
0102         m_jac_transport = internal_jacobian_type::identity();
0103       } else {
0104         m_jac_transport = matrix::identity<internal_jacobian_type>();
0105       }
0106     }
0107
0108     /// @returns free track parameters - non-const access
0109     DETRAY_HOST_DEVICE
0110     free_track_parameters_type &operator()() { return m_track; }
0111
0112     /// @returns free track parameters - const access
0113     DETRAY_HOST_DEVICE
0114     const free_track_parameters_type &operator()() const { return m_track; }
0115
0116     /// Get stepping direction
0117     DETRAY_HOST_DEVICE
0118     inline step::direction direction() const {
0119       return m_step_size >= 0.f ? step::direction::e_forward
0120                                 : step::direction::e_backward;
0121     }
0122
0123     /// Updates the total number of step trials
0124     DETRAY_HOST_DEVICE inline void count_trials() { ++m_n_total_trials; }
0125
0126     /// @returns the total number of step trials. For steppers that don't
0127     /// use adaptive step size scaling, this is the number of steps
0128     DETRAY_HOST_DEVICE inline std::size_t n_total_trials() const {
0129       return m_n_total_trials;
0130     }
0131
0132     /// Set the current step size
0133     DETRAY_HOST_DEVICE
0134     inline void set_step_size(const scalar_type step) {
0135       assert(math::fabs(step) >= 1e-4f * unit<scalar_type>::mm);
0136       m_step_size = step;
0137     }
0138
0139     /// @returns the step size of the current step.
0140     DETRAY_HOST_DEVICE
0141     inline scalar_type step_size() const { return m_step_size; }
0142
0143     /// @returns this states path length.
0144     DETRAY_HOST_DEVICE
0145     inline scalar_type path_length() const { return m_path_length; }
0146
0147     /// @returns this states total integration path length.
0148     DETRAY_HOST_DEVICE
0149     inline scalar_type abs_path_length() const { return m_abs_path_length; }
0150
0151     /// Add a new segment to all path lengths (forward or backward)
0152     DETRAY_HOST_DEVICE inline void update_path_lengths(scalar_type seg) {
0153       m_path_length += seg;
0154       m_abs_path_length += math::fabs(seg);
0155     }
0156
0157     /// Set new step constraint
0158     template <step::constraint type = step::constraint::e_actor>
0159     DETRAY_HOST_DEVICE inline void set_constraint(scalar_type step_size) {
0160       m_constraint.template set<type>(step_size);
0161     }
0162
0163     /// @returns access to this states step constraints
0164     DETRAY_HOST_DEVICE
0165     inline const constraint_t &constraints() const { return m_constraint; }
0166
0167     /// Remove [all] constraints
0168     template <step::constraint type = step::constraint::e_actor>
0169     DETRAY_HOST_DEVICE inline void release_step() {
0170       m_constraint.template release<type>();
0171     }
0172
0173     /// @returns the current particle hypothesis
0174     DETRAY_HOST_DEVICE
0175     const auto &particle_hypothesis() const { return m_ptc; }
0176
0177     /// Set new particle hypothesis
0178     DETRAY_HOST_DEVICE
0179     inline void set_particle(const pdg_particle<scalar_type> &ptc) {
0180       m_ptc = ptc;
0181     }
0182
0183     /// @returns the current transport Jacbian - const
0184     /// @note converts to the full 8x8 matrix in case the reduced storage
0185     DETRAY_HOST_DEVICE
0186     inline free_matrix_type transport_jacobian() const
0187       requires(!std::same_as<free_matrix_type, internal_jacobian_type>)
0188     {
0189       return m_jac_transport.operator dmatrix<algebra_t, 8, 8>();
0190     }
0191
0192     /// @returns the current transport Jacbian - const
0193     DETRAY_HOST_DEVICE
0194     inline const free_matrix_type &transport_jacobian() const
0195       requires std::same_as<free_matrix_type, internal_jacobian_type>
0196     {
0197       return m_jac_transport;
0198     }
0199
0200     /// @returns the current transport Jacbian.
0201     DETRAY_HOST_DEVICE
0202     inline free_matrix_type &transport_jacobian()
0203       requires std::same_as<free_matrix_type, internal_jacobian_type>
0204     {
0205       return m_jac_transport;
0206     }
0207
0208     /// @returns the current transport Jacbian.
0209     DETRAY_HOST_DEVICE
0210     inline internal_jacobian_type &internal_transport_jacobian() {
0211       return m_jac_transport;
0212     }
0213
0214     /// @returns the current transport Jacbian.
0215     DETRAY_HOST_DEVICE
0216     inline const internal_jacobian_type &internal_transport_jacobian() const {
0217       return m_jac_transport;
0218     }
0219
0220     /// Reset transport Jacbian.
0221     DETRAY_HOST_DEVICE
0222     inline void reset_transport_jacobian() {
0223       // HACK: When the overload resolution for the transport Jacobian
0224       // type is resolved, remove this conditional logic.
0225       if constexpr (concepts::transport_jacobian<internal_jacobian_type>) {
0226         m_jac_transport = internal_jacobian_type::identity();
0227       } else {
0228         m_jac_transport = matrix::identity<internal_jacobian_type>();
0229       }
0230     }
0231
0232     /// @returns access to this states navigation policy state
0233     DETRAY_HOST_DEVICE
0234     inline typename policy_t::state &policy_state() { return m_policy_state; }
0235
0236     /// @returns the stepping inspector
0237     DETRAY_HOST
0238     constexpr auto &inspector() { return m_inspector; }
0239
0240     /// Call the stepping inspector
0241     DETRAY_HOST_DEVICE
0242     inline void run_inspector([[maybe_unused]] const stepping::config &cfg,
0243                               [[maybe_unused]] const char *message,
0244                               [[maybe_unused]] const scalar_type dist) {
0245       if constexpr (!std::is_same_v<inspector_t, stepping::void_inspector>) {
0246         m_inspector(*this, cfg, message, dist);
0247       }
0248
0249       DETRAY_DEBUG_HOST("" << message << "\n"
0250                            << detray::stepping::print_state(*this, dist));
0251     }
0252
0253    private:
0254     /// Jacobian transport matrix
0255     internal_jacobian_type m_jac_transport;
0256
0257     /// Free track parameters
0258     free_track_parameters_type m_track;
0259
0260     /// The default particle hypothesis is 'muon'
0261     pdg_particle<scalar_type> m_ptc = muon<scalar_type>();
0262
0263     /// Total number of step trials
0264     std::size_t m_n_total_trials{0u};
0265
0266     /// Current step size
0267     scalar_type m_step_size{0.f};
0268
0269     /// Track path length (current position along track)
0270     scalar_type m_path_length{0.f};
0271
0272     /// Absolute path length (total path length covered by the integration)
0273     scalar_type m_abs_path_length{0.f};
0274
0275     /// Step size constraints (optional)
0276     DETRAY_NO_UNIQUE_ADDRESS constraint_t m_constraint = {};
0277
0278     /// Navigation policy state
0279     DETRAY_NO_UNIQUE_ADDRESS typename policy_t::state m_policy_state = {};
0280
0281     /// The inspector type of the stepping (for debugging only)
0282     DETRAY_NO_UNIQUE_ADDRESS inspector_type m_inspector;
0283   };
0284 };
0285
0286 }  // namespace detray