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