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