detray/tracks/free_track_parameters.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/definitions/track_parametrization.hpp"
0016
0017 // System include(s)
0018 #include <ostream>
0019
0020 namespace detray {
0021
0022 template <concepts::algebra algebra_t>
0023 struct free_parameters_vector {
0024   /// @name Type definitions for the struct
0025   /// @{
0026   using algebra_type = algebra_t;
0027   using scalar_type = dscalar<algebra_t>;
0028   using point3_type = dpoint3D<algebra_t>;
0029   using vector3_type = dvector3D<algebra_t>;
0030
0031   // Shorthand vector type related to free track parameters.
0032   using vector_type = free_vector<algebra_t>;
0033
0034   /// @}
0035
0036   /// Default constructor
0037   free_parameters_vector() = default;
0038
0039   /// Construct from a 6-dim vector of parameters
0040   DETRAY_HOST_DEVICE
0041   explicit free_parameters_vector(const vector_type& vec) : m_vector(vec) {
0042     assert(!this->is_invalid());
0043   }
0044
0045   /// Construct from single parameters
0046   ///
0047   /// @param pos the global position
0048   /// @param time the time
0049   /// @param mom the global track momentum 3-vector
0050   /// @param q the particle charge
0051   DETRAY_HOST_DEVICE
0052   free_parameters_vector(const point3_type& pos, const scalar_type time,
0053                          const vector3_type& mom, const scalar_type q) {
0054     getter::set_block(m_vector, pos, e_free_pos0, 0u);
0055     getter::element(m_vector, e_free_time, 0u) = time;
0056
0057     scalar_type p = vector::norm(mom);
0058     vector3_type mom_norm = vector::normalize(mom);
0059     getter::set_block(m_vector, mom_norm, e_free_dir0, 0u);
0060     getter::element(m_vector, e_free_qoverp, 0u) = q / p;
0061
0062     assert(!this->is_invalid());
0063   }
0064
0065   /// @param rhs is the left hand side params for comparison
0066   DETRAY_HOST_DEVICE
0067   bool operator==(const free_parameters_vector& rhs) const {
0068     for (unsigned int i = 0u; i < e_free_size; i++) {
0069       if (math::fabs((*this)[i] - rhs[i]) >
0070           std::numeric_limits<scalar_type>::epsilon()) {
0071         return false;
0072       }
0073     }
0074
0075     return true;
0076   }
0077
0078   /// Convenience access to the track parameters - const
0079   DETRAY_HOST_DEVICE
0080   scalar_type operator[](std::size_t i) const {
0081     return getter::element(m_vector, static_cast<unsigned int>(i), 0u);
0082   }
0083
0084   /// Convenience access to the track parameters - non-const
0085   DETRAY_HOST_DEVICE
0086   decltype(auto) operator[](std::size_t i) {
0087     return getter::element(m_vector, static_cast<unsigned int>(i), 0u);
0088   }
0089
0090   /// @returns the global track position
0091   DETRAY_HOST_DEVICE
0092   point3_type pos() const {
0093     return {getter::element(m_vector, e_free_pos0, 0u),
0094             getter::element(m_vector, e_free_pos1, 0u),
0095             getter::element(m_vector, e_free_pos2, 0u)};
0096   }
0097
0098   /// Set the global track position
0099   DETRAY_HOST_DEVICE
0100   void set_pos(const vector3_type& pos) {
0101     assert(math::isfinite(pos[0]));
0102     assert(math::isfinite(pos[1]));
0103     assert(math::isfinite(pos[2]));
0104     getter::set_block(m_vector, pos, e_free_pos0, 0u);
0105   }
0106
0107   /// @returns the normalized, global track direction
0108   DETRAY_HOST_DEVICE
0109   vector3_type dir() const {
0110     return {getter::element(m_vector, e_free_dir0, 0u),
0111             getter::element(m_vector, e_free_dir1, 0u),
0112             getter::element(m_vector, e_free_dir2, 0u)};
0113   }
0114
0115   /// Set the global track direction
0116   /// @note Must be normalized!
0117   DETRAY_HOST_DEVICE
0118   void set_dir(const vector3_type& dir) {
0119     assert(math::isfinite(dir[0]));
0120     assert(math::isfinite(dir[1]));
0121     assert(math::isfinite(dir[2]));
0122     assert(algebra::approx_equal(vector::norm(dir), scalar_type{1},
0123                                  scalar_type{1e-5f}));
0124     getter::set_block(m_vector, dir, e_free_dir0, 0u);
0125   }
0126
0127   /// @returns the time
0128   DETRAY_HOST_DEVICE
0129   scalar_type time() const {
0130     return getter::element(m_vector, e_free_time, 0u);
0131   }
0132
0133   /// Set the time
0134   DETRAY_HOST_DEVICE
0135   void set_time(const scalar_type t) {
0136     assert(math::isfinite(t));
0137     getter::element(m_vector, e_free_time, 0u) = t;
0138   }
0139
0140   /// @returns the q/p value
0141   DETRAY_HOST_DEVICE
0142   scalar_type qop() const {
0143     return getter::element(m_vector, e_free_qoverp, 0u);
0144   }
0145
0146   /// Set the q/p value
0147   DETRAY_HOST_DEVICE
0148   void set_qop(const scalar_type qop) {
0149     assert(math::isfinite(qop));
0150     getter::element(m_vector, e_free_qoverp, 0u) = qop;
0151   }
0152
0153   /// @returns the q/p_T value
0154   DETRAY_HOST_DEVICE
0155   scalar_type qopT() const {
0156     const vector3_type dir = this->dir();
0157     assert(vector::perp(dir) != 0.f);
0158     return getter::element(m_vector, e_free_qoverp, 0u) / vector::perp(dir);
0159   }
0160
0161   /// @returns the q/p_z value
0162   DETRAY_HOST_DEVICE
0163   scalar_type qopz() const {
0164     const vector3_type dir = this->dir();
0165     return getter::element(m_vector, e_free_qoverp, 0u) / dir[2];
0166   }
0167
0168   /// @returns the absolute momentum
0169   DETRAY_HOST_DEVICE
0170   scalar_type p(const scalar_type q) const {
0171     assert(qop() != 0.f);
0172     assert(q * qop() > 0.f);
0173     return q / qop();
0174   }
0175
0176   /// @returns the global momentum 3-vector
0177   DETRAY_HOST_DEVICE
0178   vector3_type mom(const scalar_type q) const { return p(q) * dir(); }
0179
0180   /// @returns the transverse momentum
0181   DETRAY_HOST_DEVICE
0182   scalar_type pT(const scalar_type q) const {
0183     assert(this->qop() != 0.f);
0184     assert(q * qop() > 0.f);
0185     return math::fabs(q / this->qop() * vector::perp(this->dir()));
0186   }
0187
0188   /// @returns the absolute momentum z-component
0189   DETRAY_HOST_DEVICE
0190   scalar_type pz(const scalar_type q) const {
0191     assert(this->qop() != 0.f);
0192     assert(q * qop() > 0.f);
0193     return math::fabs(q / this->qop() * this->dir()[2]);
0194   }
0195
0196   /// @returns true if the parameter vector contains invalid elements or the
0197   /// direction is not normalized
0198   DETRAY_HOST_DEVICE
0199   constexpr bool is_invalid() const {
0200     bool inv_elem{false};
0201     for (std::size_t i = 0u; i < e_free_size; ++i) {
0202       inv_elem = inv_elem || !math::isfinite(getter::element(m_vector, i, 0u));
0203     }
0204     return inv_elem ||
0205            !algebra::approx_equal(vector::norm(dir()), scalar_type{1},
0206                                   scalar_type{1e-5f});
0207   }
0208
0209  private:
0210   /// Transform to a string for debugging output
0211   DETRAY_HOST
0212   friend std::ostream& operator<<(std::ostream& out_stream,
0213                                   const free_parameters_vector& fparam) {
0214     out_stream << fparam.m_vector;
0215     return out_stream;
0216   }
0217
0218   vector_type m_vector = matrix::zero<vector_type>();
0219 };
0220
0221 /// The free track parameters consist only of the parameter vector
0222 template <concepts::algebra algebra_t>
0223 using free_track_parameters = free_parameters_vector<algebra_t>;
0224
0225 }  // namespace detray