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
0024 /// @tparam Derived Eigen derived concrete type
0025 /// @param v Any vector like Eigen type, static or dynamic
0026 /// @note Will static assert that the number of rows of @p v is at least 2, or
0027 /// in case of dynamic size, will abort execution if that is not the case.
0028 /// @return The value of the angle in the transverse plane.
0029 template <typename Derived>
0030 double phi(const Eigen::MatrixBase<Derived>& v) noexcept {
0031   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0032   if constexpr (rows != -1) {
0033     // static size, do compile time check
0034     static_assert(rows >= 2,
0035                   "Phi function not valid for vectors not at least 2D");
0036   } else {
0037     // dynamic size
0038     assert(v.rows() >= 2 &&
0039            "Phi function not valid for vectors not at least 2D");
0040   }
0041   return std::atan2(v[1], v[0]);
0042 }
0043
0044 /// Calculate phi (transverse plane angle) from anything implementing a method
0045 /// like `phi()` returning anything convertible to `double`.
0046 /// @tparam T anything that has a phi method
0047 /// @param v Any type that implements a phi method
0048 /// @return The phi value
0049 template <typename T>
0050 double phi(const T& v) noexcept
0051   requires requires {
0052     { v.phi() } -> std::floating_point;
0053   }
0054 {
0055   return v.phi();
0056 }
0057
0058 /// Calculate radius in the transverse (xy) plane of a vector
0059 /// @tparam Derived Eigen derived concrete type
0060 /// @param v Any vector like Eigen type, static or dynamic
0061 /// @note Will static assert that the number of rows of @p v is at least 2, or
0062 /// in case of dynamic size, will abort execution if that is not the case.
0063 /// @return The transverse radius value.
0064 template <typename Derived>
0065 double perp(const Eigen::MatrixBase<Derived>& v) noexcept {
0066   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0067   if constexpr (rows != -1) {
0068     // static size, do compile time check
0069     static_assert(rows >= 2,
0070                   "Perp function not valid for vectors not at least 2D");
0071   } else {
0072     // dynamic size
0073     assert(v.rows() >= 2 &&
0074            "Perp function not valid for vectors not at least 2D");
0075   }
0076   return v.template head<2>().norm();
0077 }
0078
0079 /// Calculate the theta angle (longitudinal w.r.t. z axis) of a vector
0080 /// @tparam Derived Eigen derived concrete type
0081 /// @param v Any vector like Eigen type, static or dynamic
0082 /// @note Will static assert that the number of rows of @p v is at least 3, or
0083 /// in case of dynamic size, will abort execution if that is not the case.
0084 /// @return The theta value
0085 template <typename Derived>
0086 double theta(const Eigen::MatrixBase<Derived>& v) noexcept {
0087   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0088   if constexpr (rows != -1) {
0089     // static size, do compile time check
0090     static_assert(rows == 3, "Theta function not valid for non-3D vectors.");
0091   } else {
0092     // dynamic size
0093     assert(v.rows() == 3 && "Theta function not valid for non-3D vectors.");
0094   }
0095
0096   return std::atan2(perp(v), v[2]);
0097 }
0098
0099 /// Calculate the pseudorapidity for a vector.
0100 /// @tparam Derived Eigen derived concrete type
0101 /// @param v Any vector like Eigen type, static or dynamic
0102 /// @note Will static assert that the number of rows of @p v is at least 3, or
0103 /// in case of dynamic size, will abort execution if that is not the case.
0104 /// @return The pseudorapidity value
0105 template <typename Derived>
0106 double eta(const Eigen::MatrixBase<Derived>& v) noexcept {
0107   constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
0108   if constexpr (rows != -1) {
0109     // static size, do compile time check
0110     static_assert(rows == 3, "Eta function not valid for non-3D vectors.");
0111   } else {
0112     // dynamic size
0113     assert(v.rows() == 3 && "Eta function not valid for non-3D vectors.");
0114   }
0115
0116   if (v[0] == 0. && v[1] == 0.) {
0117     return std::copysign(std::numeric_limits<double>::infinity(), v[2]);
0118   } else {
0119     return std::asinh(v[2] / perp(v));
0120   }
0121 }
0122
0123 /// @brief Fast evaluation of trigonomic functions.
0124 ///
0125 /// @param direction for this evaluatoin
0126 ///
0127 /// @return cos(phi), sin(phi), cos(theta), sin(theta), 1/sin(theta)
0128 inline std::array<double, 4> evaluateTrigonomics(const Vector3& direction) {
0129   const double x = direction(0);  // == cos(phi) * sin(theta)
0130   const double y = direction(1);  // == sin(phi) * sin(theta)
0131   const double z = direction(2);  // == cos(theta)
0132   // can be turned into cosine/sine
0133   const double cosTheta = z;
0134   const double sinTheta = std::sqrt(1 - z * z);
0135   assert(sinTheta != 0 &&
0136          "VectorHelpers: Vector is parallel to the z-axis "
0137          "which leads to division by zero");
0138   const double invSinTheta = 1. / sinTheta;
0139   const double cosPhi = x * invSinTheta;
0140   const double sinPhi = y * invSinTheta;
0141
0142   return {cosPhi, sinPhi, cosTheta, sinTheta};
0143 }
0144
0145 /// Helper method to extract the binning value from a 3D vector.
0146 ///
0147 /// For this method a 3D vector is required to guarantee all potential
0148 /// axis directions to be casted from
0149 ///
0150 /// @param position is the position in global
0151 /// @param aDir is the axis direction to be extracted
0152 ///
0153 /// @return the value of the binning direction
0154 inline double cast(const Vector3& position, AxisDirection aDir) {
0155   using enum AxisDirection;
0156   switch (aDir) {
0157     case AxisX:
0158       return position[0];
0159     case AxisY:
0160       return position[1];
0161     case AxisZ:
0162       return position[2];
0163     case AxisR:
0164       return perp(position);
0165     case AxisPhi:
0166       return phi(position);
0167     case AxisRPhi:
0168       return perp(position) * phi(position);
0169     case AxisTheta:
0170       return theta(position);
0171     case AxisEta:
0172       return eta(position);
0173     case AxisMag:
0174       return position.norm();
0175     default:
0176       assert(false && "Invalid AxisDirection 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 SquareMatrix3 cross(const SquareMatrix3& m, const Vector3& v) {
0188   SquareMatrix3 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