Acts/Utilities/VectorHelpers.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2016-2023 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 http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/Common.hpp"
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014 #include "Acts/Utilities/BinningType.hpp"
0015 #include "Acts/Utilities/TypeTraits.hpp"
0016
0017 #include <array>
0018 #include <limits>
0019
0020 #include "Eigen/Dense"
0021
0022 namespace Acts::VectorHelpers {
0023
0024 namespace detail {
0025 template <class T>
0026 using phi_method_t = decltype(std::declval<const T>().phi());
0027
0028 template <class T>
0029 using has_phi_method = Concepts::is_detected<phi_method_t, T>;
0030 }  // namespace detail
0031
0032 /// Calculate phi (transverse plane angle) from compatible Eigen types
0033 /// @tparam Derived Eigen derived concrete type
0034 /// @param v Any vector like Eigen type, static or dynamic
0035 /// @note Will static assert that the number of rows of @p v is at least 2, or
0036 /// in case of dynamic size, will abort execution if that is not the case.
0037 /// @return The value of the angle in the transverse plane.
0038 template <typename Derived>
0039 double phi(const Eigen::MatrixBase<Derived>& v) noexcept {
0040   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0041   if constexpr (rows != -1) {
0042     // static size, do compile time check
0043     static_assert(rows >= 2,
0044                   "Phi function not valid for vectors not at least 2D");
0045   } else {
0046     // dynamic size
0047     assert(v.rows() >= 2 &&
0048            "Phi function not valid for vectors not at least 2D");
0049   }
0050   return std::atan2(v[1], v[0]);
0051 }
0052
0053 /// Calculate phi (transverse plane angle) from anything implementing a method
0054 /// like `phi()` returning anything convertible to `double`.
0055 /// @tparam T anything that has a phi method
0056 /// @param v Any type that implements a phi method
0057 /// @return The phi value
0058 template <typename T,
0059           std::enable_if_t<detail::has_phi_method<T>::value, int> = 0>
0060 double phi(const T& v) noexcept {
0061   return v.phi();
0062 }
0063
0064 /// Calculate radius in the transverse (xy) plane of a vector
0065 /// @tparam Derived Eigen derived concrete type
0066 /// @param v Any vector like Eigen type, static or dynamic
0067 /// @note Will static assert that the number of rows of @p v is at least 2, or
0068 /// in case of dynamic size, will abort execution if that is not the case.
0069 /// @return The transverse radius value.
0070 template <typename Derived>
0071 double perp(const Eigen::MatrixBase<Derived>& v) noexcept {
0072   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0073   if constexpr (rows != -1) {
0074     // static size, do compile time check
0075     static_assert(rows >= 2,
0076                   "Perp function not valid for vectors not at least 2D");
0077   } else {
0078     // dynamic size
0079     assert(v.rows() >= 2 &&
0080            "Perp function not valid for vectors not at least 2D");
0081   }
0082   return std::hypot(v[0], v[1]);
0083 }
0084
0085 /// Calculate the theta angle (longitudinal w.r.t. z axis) of a vector
0086 /// @tparam Derived Eigen derived concrete type
0087 /// @param v Any vector like Eigen type, static or dynamic
0088 /// @note Will static assert that the number of rows of @p v is at least 3, or
0089 /// in case of dynamic size, will abort execution if that is not the case.
0090 /// @return The theta value
0091 template <typename Derived>
0092 double theta(const Eigen::MatrixBase<Derived>& v) noexcept {
0093   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0094   if constexpr (rows != -1) {
0095     // static size, do compile time check
0096     static_assert(rows == 3, "Theta function not valid for non-3D vectors.");
0097   } else {
0098     // dynamic size
0099     assert(v.rows() == 3 && "Theta function not valid for non-3D vectors.");
0100   }
0101
0102   return std::atan2(perp(v), v[2]);
0103 }
0104
0105 /// Calculate the pseudorapidity for a vector.
0106 /// @tparam Derived Eigen derived concrete type
0107 /// @param v Any vector like Eigen type, static or dynamic
0108 /// @note Will static assert that the number of rows of @p v is at least 3, or
0109 /// in case of dynamic size, will abort execution if that is not the case.
0110 /// @return The pseudorapidity value
0111 template <typename Derived>
0112 double eta(const Eigen::MatrixBase<Derived>& v) noexcept {
0113   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0114   if constexpr (rows != -1) {
0115     // static size, do compile time check
0116     static_assert(rows == 3, "Eta function not valid for non-3D vectors.");
0117   } else {
0118     // dynamic size
0119     assert(v.rows() == 3 && "Eta function not valid for non-3D vectors.");
0120   }
0121
0122   if (v[0] == 0. && v[1] == 0.) {
0123     return std::copysign(std::numeric_limits<double>::infinity(), v[2]);
0124   } else {
0125     return std::asinh(v[2] / perp(v));
0126   }
0127 }
0128
0129 /// @brief Fast evaluation of trigonomic functions.
0130 ///
0131 /// @param direction for this evaluatoin
0132 ///
0133 /// @return cos(phi), sin(phi), cos(theta), sin(theta), 1/sin(theta)
0134 inline std::array<ActsScalar, 4> evaluateTrigonomics(const Vector3& direction) {
0135   const ActsScalar x = direction(0);  // == cos(phi) * sin(theta)
0136   const ActsScalar y = direction(1);  // == sin(phi) * sin(theta)
0137   const ActsScalar z = direction(2);  // == cos(theta)
0138   // can be turned into cosine/sine
0139   const ActsScalar cosTheta = z;
0140   const ActsScalar sinTheta = std::sqrt(1 - z * z);
0141   assert(sinTheta != 0 &&
0142          "VectorHelpers: Vector is parallel to the z-axis "
0143          "which leads to division by zero");
0144   const ActsScalar invSinTheta = 1. / sinTheta;
0145   const ActsScalar cosPhi = x * invSinTheta;
0146   const ActsScalar sinPhi = y * invSinTheta;
0147
0148   return {cosPhi, sinPhi, cosTheta, sinTheta};
0149 }
0150
0151 /// Helper method to extract the binning value from a 3D vector.
0152 ///
0153 /// For this method a 3D vector is required to guarantee all potential
0154 /// binning values.
0155 inline double cast(const Vector3& position, BinningValue bval) {
0156   switch (bval) {
0157     case binX:
0158       return position[0];
0159     case binY:
0160       return position[1];
0161     case binZ:
0162       return position[2];
0163     case binR:
0164       return perp(position);
0165     case binPhi:
0166       return phi(position);
0167     case binRPhi:
0168       return perp(position) * phi(position);
0169     case binH:
0170       return theta(position);
0171     case binEta:
0172       return eta(position);
0173     case binMag:
0174       return position.norm();
0175     default:
0176       assert(false && "Invalid BinningValue enum value");
0177       return std::numeric_limits<double>::quiet_NaN();
0178   }
0179 }
0180
0181 /// @brief Calculates column-wise cross products of a matrix and a vector and
0182 /// stores the result column-wise in a matrix.
0183 ///
0184 /// @param [in] m Matrix that will be used for cross products
0185 /// @param [in] v Vector for cross products
0186 /// @return Constructed matrix
0187 inline ActsMatrix<3, 3> cross(const ActsMatrix<3, 3>& m, const Vector3& v) {
0188   ActsMatrix<3, 3> r;
0189   r.col(0) = m.col(0).cross(v);
0190   r.col(1) = m.col(1).cross(v);
0191   r.col(2) = m.col(2).cross(v);
0192
0193   return r;
0194 }
0195
0196 /// Access the three-position components in a four-position vector.
0197 inline auto position(const Vector4& pos4) {
0198   return pos4.segment<3>(ePos0);
0199 }
0200
0201 /// Access the three-position components in a free parameters vector.
0202 inline auto position(const FreeVector& params) {
0203   return params.segment<3>(eFreePos0);
0204 }
0205
0206 /// Construct a four-vector from a three-vector and scalar fourth component.
0207 template <typename vector3_t>
0208 inline auto makeVector4(const Eigen::MatrixBase<vector3_t>& vec3,
0209                         typename vector3_t::Scalar w)
0210     -> Eigen::Matrix<typename vector3_t::Scalar, 4, 1> {
0211   EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(vector3_t, 3);
0212
0213   Eigen::Matrix<typename vector3_t::Scalar, 4, 1> vec4;
0214   vec4[ePos0] = vec3[ePos0];
0215   vec4[ePos1] = vec3[ePos1];
0216   vec4[ePos2] = vec3[ePos2];
0217   vec4[eTime] = w;
0218   return vec4;
0219 }
0220
0221 /// Calculate the incident angles of a vector with in a given reference frame
0222 /// @tparam Derived Eigen derived concrete type
0223 /// @param direction The crossing direction in the global frame
0224 /// @param globalToLocal Rotation from global to local frame
0225 /// @return The angles of incidence in the two normal planes
0226 inline std::pair<double, double> incidentAngles(
0227     const Acts::Vector3& direction,
0228     const Acts::RotationMatrix3& globalToLocal) {
0229   Acts::Vector3 trfDir = globalToLocal * direction;
0230   // The angles are defined with respect to the measurement axis
0231   // i.e. "head-on" == pi/2, parallel = 0
0232   double phi = std::atan2(trfDir[2], trfDir[0]);
0233   double theta = std::atan2(trfDir[2], trfDir[1]);
0234   return {phi, theta};
0235 }
0236
0237 }  // namespace Acts::VectorHelpers