detray/tracks/helix.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 // TODO: Remove this when gcc fixes their false positives.
0012 #if defined(__GNUC__) && !defined(__clang__)
0013 #pragma GCC diagnostic warning "-Wmaybe-uninitialized"
0014 #endif
0015
0016 // Project include(s).
0017 #include "detray/definitions/algebra.hpp"
0018 #include "detray/definitions/math.hpp"
0019 #include "detray/tracks/free_track_parameters.hpp"
0020 #include "detray/utils/invalid_values.hpp"
0021
0022 // System include(s).
0023 #include <ostream>
0024
0025 namespace detray::detail {
0026
0027 /// @brief describes a helical trajectory in a given B-field.
0028 ///
0029 /// Helix class for the analytical solution of track propagation in
0030 /// homogeneous B field. This Follows the notation of Eq (4.7) in
0031 /// DOI:10.1007/978-3-030-65771-0
0032 template <concepts::algebra algebra_t>
0033 class helix {
0034  public:
0035   using algebra_type = algebra_t;
0036   using scalar_type = dscalar<algebra_t>;
0037   using point3_type = dpoint3D<algebra_t>;
0038   using vector3_type = dvector3D<algebra_t>;
0039   using transform3_type = dtransform3D<algebra_t>;
0040
0041   /// Free track parameters
0042   using free_track_parameters_type = free_track_parameters<algebra_t>;
0043
0044   /// 2D Matrix type
0045   template <std::size_t ROWS, std::size_t COLS>
0046   using matrix_type = dmatrix<algebra_t, ROWS, COLS>;
0047   using free_matrix_t = free_matrix<algebra_t>;
0048
0049   DETRAY_HOST_DEVICE
0050   helix() = delete;
0051
0052   /// Parametrized constructor
0053   ///
0054   /// @param pos the the origin of the helix
0055   /// @param time the time parameter
0056   /// @param dir the initial normalized direction for the helix
0057   /// @param q the charge of the particle
0058   /// @param mag_field the magnetic field vector
0059   DETRAY_HOST_DEVICE
0060   helix(const point3_type &pos, const scalar_type time, const vector3_type &dir,
0061         const scalar_type qop, const vector3_type &mag_field)
0062       : _pos(pos),
0063         _time(time),
0064         _qop(qop),
0065         _B{vector::norm(mag_field)},
0066         _h0{vector::normalize(mag_field)},
0067         _t0{dir} {
0068     assert((math::fabs(vector::norm(_t0) - 1.f) < 1e-5f) &&
0069            "The helix direction must be normalized");
0070
0071     // Momentum
0072     const vector3_type mom =
0073         (1.f / static_cast<scalar_type>(math::fabs(qop))) * _t0;
0074
0075     // Normalized _h0 X _t0
0076     _n0 = vector::normalize(vector::cross(_h0, _t0));
0077
0078     // Magnitude of _h0 X _t0
0079     _alpha = vector::norm(vector::cross(_h0, _t0));
0080
0081     // Dot product of _h0 X _t0
0082     _delta = vector::dot(_h0, _t0);
0083
0084     // Path length scaler
0085     _K = -_qop * _B;
0086
0087     // Get longitudinal momentum parallel to B field
0088     scalar_type pz = vector::dot(mom, _h0);
0089
0090     // Get transverse momentum perpendicular to B field
0091     vector3_type pT = mom - pz * _h0;
0092
0093     // R [mm] =  pT [GeV] / B [T] in natural unit
0094     _R = vector::norm(pT) / _B;
0095
0096     // Handle the case of pT ~ 0
0097     if (vector::norm(pT) < 1e-6f) {
0098       _vz_over_vt = detail::invalid_value<scalar_type>();
0099     } else {
0100       // Get vz over vt in new coordinate
0101       _vz_over_vt = pz / vector::norm(pT);
0102     }
0103   }
0104
0105   DETRAY_HOST_DEVICE
0106   helix(const free_track_parameters_type &track, const vector3_type &mag_field)
0107       : helix(track.pos(), track.time(), track.dir(), track.qop(), mag_field) {
0108     assert(!track.is_invalid());
0109   }
0110
0111   /// @TODO Add covfie field view concept
0112   template <typename field_view_t>
0113     requires(!concepts::vector3D<field_view_t>)
0114   DETRAY_HOST_DEVICE helix(const free_track_parameters_type &track,
0115                            const field_view_t mag_field)
0116       : helix(track.pos(), track.time(), track.dir(), track.qop(),
0117               sample_field(mag_field, track.pos())) {}
0118
0119   /// @returns the radius of helix
0120   DETRAY_HOST_DEVICE
0121   scalar_type radius() const { return _R; }
0122
0123   /// @returns the position after propagating the path length of s
0124   DETRAY_HOST_DEVICE
0125   point3_type operator()(const scalar_type s) const { return this->pos(s); }
0126
0127   /// @returns the position after propagating the path length of s
0128   DETRAY_HOST_DEVICE
0129   point3_type pos(const scalar_type s) const {
0130     // Handle the case of pT ~ 0
0131     if (_vz_over_vt == detail::invalid_value<scalar_type>()) {
0132       return _pos + s * _h0;
0133     }
0134
0135     point3_type ret = _pos;
0136     ret = ret + _delta / _K * (_K * s - math::sin(_K * s)) * _h0;
0137     ret = ret + math::sin(_K * s) / _K * _t0;
0138     ret = ret + _alpha / _K * (1.f - math::cos(_K * s)) * _n0;
0139
0140     return ret;
0141   }
0142
0143   DETRAY_HOST_DEVICE
0144   point3_type pos() const { return _pos; }
0145
0146   /// @returns the tangential vector after propagating the path length of s
0147   DETRAY_HOST_DEVICE
0148   vector3_type dir(const scalar_type s) const {
0149     // Handle the case of pT ~ 0
0150     if (_vz_over_vt == detail::invalid_value<scalar_type>()) {
0151       return _t0;
0152     }
0153
0154     vector3_type ret{0.f, 0.f, 0.f};
0155
0156     ret = ret + _delta * (1 - math::cos(_K * s)) * _h0;
0157     ret = ret + math::cos(_K * s) * _t0;
0158     ret = ret + _alpha * math::sin(_K * s) * _n0;
0159
0160     return vector::normalize(ret);
0161   }
0162
0163   DETRAY_HOST_DEVICE
0164   point3_type dir() const { return _t0; }
0165
0166   DETRAY_HOST_DEVICE
0167   scalar_type time() const { return _time; }
0168
0169   DETRAY_HOST_DEVICE
0170   scalar_type qop() const { return _qop; }
0171
0172   DETRAY_HOST_DEVICE
0173   scalar_type B() const { return _B; }
0174
0175   DETRAY_HOST_DEVICE
0176   vector3_type b_field() const { return _B * _h0; }
0177
0178   /// @returns the transport jacobian after propagating the path length of s
0179   DETRAY_HOST_DEVICE
0180   free_matrix_t jacobian(const scalar_type s) const {
0181     free_matrix_t ret = matrix::zero<free_matrix_t>();
0182
0183     const matrix_type<3, 3> I33 = matrix::identity<matrix_type<3, 3>>();
0184     const matrix_type<3, 3> Z33 = matrix::zero<matrix_type<3, 3>>();
0185
0186     // Notations
0187     // r = position
0188     // t = direction
0189     // l = qoverp
0190
0191     // Get drdr
0192     auto drdr = I33;
0193     getter::set_block(ret, drdr, e_free_pos0, e_free_pos0);
0194
0195     // Get dtdr
0196     auto dtdr = Z33;
0197     getter::set_block(ret, dtdr, e_free_dir0, e_free_pos0);
0198
0199     // Get drdt
0200     auto drdt = Z33;
0201
0202     const scalar_type sin_ks = math::sin(_K * s);
0203     const scalar_type cos_ks = math::cos(_K * s);
0204     drdt = drdt + sin_ks / _K * I33;
0205
0206     matrix_type<3, 1> H0 = matrix::zero<matrix_type<3, 1>>();
0207     getter::element(H0, 0u, 0u) = _h0[0u];
0208     getter::element(H0, 1u, 0u) = _h0[1u];
0209     getter::element(H0, 2u, 0u) = _h0[2u];
0210     const matrix_type<1, 3> H0_T = matrix::transpose(H0);
0211     const matrix_type<3, 3> H0H0_T = H0 * H0_T;
0212
0213     drdt = drdt + (_K * s - sin_ks) / _K * H0H0_T;
0214
0215     drdt = drdt + (cos_ks - 1.f) / _K * matrix::column_wise_cross(I33, _h0);
0216
0217     getter::set_block(ret, drdt, e_free_pos0, e_free_dir0);
0218
0219     // Get dtdt
0220     auto dtdt = Z33;
0221     dtdt = dtdt + cos_ks * I33;
0222     dtdt = dtdt + (1 - cos_ks) * H0H0_T;
0223     dtdt = dtdt - sin_ks * matrix::column_wise_cross(I33, _h0);
0224
0225     getter::set_block(ret, dtdt, e_free_dir0, e_free_dir0);
0226
0227     // Get drdl
0228     vector3_type drdl = 1.f / _qop * (s * this->dir(s) + _pos - this->pos(s));
0229
0230     getter::set_block(ret, drdl, e_free_pos0, e_free_qoverp);
0231
0232     // Get dtdl
0233     vector3_type dtdl =
0234         -_B * s *
0235         (sin_ks * ((H0H0_T - I33) * _t0) + cos_ks * vector::cross(_h0, _t0));
0236
0237     getter::set_block(ret, dtdl, e_free_dir0, e_free_qoverp);
0238
0239     // 3x3 and 7x7 element is 1 (Maybe?)
0240     getter::element(ret, e_free_time, e_free_time) = 1.f;
0241     getter::element(ret, e_free_qoverp, e_free_qoverp) = 1.f;
0242
0243     return ret;
0244   }
0245
0246  private:
0247   /// Print
0248   DETRAY_HOST
0249   friend std::ostream &operator<<(std::ostream &os, const helix &h) {
0250     os << "helix: ";
0251     os << "ori = " << h._pos;
0252     os << ", dir = " << h._t0 << std::endl;
0253
0254     return os;
0255   }
0256
0257   /// Helper to get a field strength at a given position
0258   template <typename field_view_t>
0259   constexpr vector3_type sample_field(const field_view_t field,
0260                                       const point3_type &pos) const {
0261     const auto bvec = field.at(pos[0], pos[1], pos[2]);
0262
0263     return {bvec[0], bvec[1], bvec[2]};
0264   }
0265
0266   /// origin
0267   point3_type _pos;
0268
0269   /// time
0270   scalar_type _time;
0271
0272   /// qop
0273   scalar_type _qop;
0274
0275   /// B field strength
0276   scalar_type _B;
0277
0278   /// Normalized b field
0279   vector3_type _h0;
0280
0281   /// Normalized tangent vector
0282   vector3_type _t0;
0283
0284   /// Normalized _h0 X _t0
0285   vector3_type _n0;
0286
0287   /// Magnitude of _h0 X _t0
0288   scalar_type _alpha;
0289
0290   /// Dot product of _h0 X _t0
0291   scalar_type _delta;
0292
0293   /// Path length scaler
0294   scalar_type _K;
0295
0296   /// Radius [mm] of helix
0297   scalar_type _R;
0298
0299   /// Velocity in new z axis divided by transverse velocity
0300   scalar_type _vz_over_vt;
0301 };
0302
0303 template <concepts::algebra algebra_t, typename field_view_t>
0304   requires(!concepts::vector3D<field_view_t>)
0305 DETRAY_HOST_DEVICE helix(const free_track_parameters<algebra_t> &,
0306                          const field_view_t) -> helix<algebra_t>;
0307
0308 }  // namespace detray::detail