detray/propagator/line_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/math.hpp"
0015 #include "detray/material/material.hpp"
0016 #include "detray/navigation/policies.hpp"
0017 #include "detray/propagator/base_stepper.hpp"
0018 #include "detray/utils/logging.hpp"
0019
0020 namespace detray {
0021
0022 /// Straight line stepper implementation
0023 template <concepts::algebra algebra_t,
0024           typename constraint_t = unconstrained_step<dscalar<algebra_t>>,
0025           typename policy_t = stepper_default_policy<dscalar<algebra_t>>,
0026           typename inspector_t = stepping::void_inspector>
0027 class line_stepper final
0028     : public base_stepper<algebra_t, constraint_t, policy_t, inspector_t> {
0029   using base_type =
0030       base_stepper<algebra_t, constraint_t, policy_t, inspector_t>;
0031
0032  public:
0033   static constexpr bool uses_gradient = true;
0034
0035   using algebra_type = algebra_t;
0036   using scalar_type = dscalar<algebra_t>;
0037   using vector3_type = dvector3D<algebra_t>;
0038   using transform3_type = dtransform3D<algebra_t>;
0039   using free_track_parameters_type =
0040       typename base_type::free_track_parameters_type;
0041   using bound_track_parameters_type =
0042       typename base_type::bound_track_parameters_type;
0043   template <std::size_t ROWS, std::size_t COLS>
0044   using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
0045
0046   struct state : public base_type::template state<> {
0047     static constexpr const stepping::id id = stepping::id::e_linear;
0048
0049     using base_state = typename base_type::template state<>;
0050
0051     /// Import base class constructors
0052     using base_state::base_state;
0053
0054     /// Update the track state in a straight line.
0055     DETRAY_HOST_DEVICE
0056     inline void advance_track() {
0057       auto& track = (*this)();
0058       track.set_pos(track.pos() + track.dir() * this->step_size());
0059
0060       this->update_path_lengths(this->step_size());
0061     }
0062
0063     DETRAY_HOST_DEVICE
0064     inline void advance_jacobian() {
0065       DETRAY_VERBOSE_HOST_DEVICE("Advance Jacobian");
0066
0067       // The step transport matrix in global coordinates
0068       free_matrix<algebra_t> D = matrix::identity<free_matrix<algebra_t>>();
0069
0070       // d(x,y,z)/d(n_x,n_y,n_z)
0071       matrix_type<3, 3> dxdn =
0072           this->step_size() * matrix::identity<matrix_type<3, 3>>();
0073       getter::set_block(D, dxdn, e_free_pos0, e_free_dir0);
0074
0075       /// NOTE: Let's skip the element for d(time)/d(qoverp) for the
0076       /// moment..
0077       this->transport_jacobian() = D * this->transport_jacobian();
0078     }
0079
0080     DETRAY_HOST_DEVICE
0081     constexpr vector3_type dtds() const { return {0.f, 0.f, 0.f}; }
0082
0083     DETRAY_HOST_DEVICE
0084     constexpr scalar_type dqopds(
0085         const material<scalar_type>* /*unused*/) const {
0086       return 0.f;
0087     }
0088   };
0089
0090   /// Take a step, regulared by a constrained step
0091   ///
0092   /// @param dist_to_next The straight line distance to the next surface
0093   /// @param stepping The state object of a stepper
0094   /// @param cfg The stepping configuration
0095   ///
0096   /// @returns returning the heartbeat, indicating if the stepping is alive
0097   DETRAY_HOST_DEVICE bool step(
0098       const scalar_type dist_to_next, state& stepping,
0099       const stepping::config& cfg, const bool /*unused*/ = true,
0100       const material<scalar_type>* /*unused*/ = nullptr) const {
0101     // In case of an overlap do nothing
0102     if (math::fabs(dist_to_next) <= 1.f * unit<scalar_type>::um) {
0103       stepping.run_inspector(cfg, "Step skipped (Overlap): ", dist_to_next);
0104       return true;
0105     }
0106
0107     // Straight line stepping: The distance given by the navigator is exact
0108     stepping.set_step_size(dist_to_next);
0109
0110     DETRAY_VERBOSE_HOST_DEVICE("Take step: %f mm", stepping.step_size());
0111
0112     // Check constraints
0113     if (const scalar_type max_step =
0114             stepping.constraints().template size<>(stepping.direction());
0115         math::fabs(stepping.step_size()) > math::fabs(max_step)) {
0116       // Run inspection before step size is cut
0117       stepping.run_inspector(cfg, "Before constraint: ", dist_to_next);
0118
0119       stepping.set_step_size(max_step);
0120     }
0121
0122     // Update track state
0123     stepping.advance_track();
0124
0125     // Advance jacobian transport
0126     if (cfg.do_covariance_transport) {
0127       stepping.advance_jacobian();
0128     }
0129
0130     // Count the number of steps
0131     stepping.count_trials();
0132
0133     // Run inspection if needed
0134     stepping.run_inspector(cfg, "Step complete: ", dist_to_next);
0135
0136     return true;
0137   }
0138 };
0139
0140 }  // namespace detray