Acts/Utilities/UnitVectors.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
0013 #include <cmath>
0014 #include <limits>
0015 #include <utility>
0016
0017 namespace Acts {
0018
0019 /// Construct a normalized direction vector from phi angle and pseudorapidity.
0020 ///
0021 /// @param phi is the direction angle in the x-y plane.
0022 /// @param eta is the pseudorapidity towards the z-axis.
0023 ///
0024 /// @note The input arguments intentionally use the same template type so that
0025 ///       a compile error occurs if inconsistent input types are used. Avoids
0026 ///       unexpected implicit type conversions and forces the user to
0027 ///       explicitly cast mismatched input types.
0028 template <typename T>
0029 inline Eigen::Matrix<T, 3, 1> makeDirectionFromPhiEta(T phi, T eta) {
0030   const auto coshEtaInv = 1 / std::cosh(eta);
0031   return {
0032       std::cos(phi) * coshEtaInv,
0033       std::sin(phi) * coshEtaInv,
0034       std::tanh(eta),
0035   };
0036 }
0037
0038 /// Construct a normalized direction vector from phi and theta angle.
0039 ///
0040 /// @param phi is the direction angle in radian in the x-y plane.
0041 /// @param theta is the polar angle in radian towards the z-axis.
0042 ///
0043 /// @note The input arguments intentionally use the same template type so that
0044 ///       a compile error occurs if inconsistent input types are used. Avoids
0045 ///       unexpected implicit type conversions and forces the user to
0046 ///       explicitly cast mismatched input types.
0047 template <typename T>
0048 inline Eigen::Matrix<T, 3, 1> makeDirectionFromPhiTheta(T phi, T theta) {
0049   const auto cosTheta = std::cos(theta);
0050   const auto sinTheta = std::sin(theta);
0051   return {
0052       std::cos(phi) * sinTheta,
0053       std::sin(phi) * sinTheta,
0054       cosTheta,
0055   };
0056 }
0057
0058 /// Construct a phi and theta angle from a direction vector.
0059 ///
0060 /// @param unitDir 3D vector indicating a direction
0061 ///
0062 template <typename T>
0063 inline Eigen::Matrix<T, 2, 1> makePhiThetaFromDirection(
0064     Eigen::Matrix<T, 3, 1> unitDir) {
0065   unitDir.normalize();
0066   T phi = std::atan2(unitDir[1], unitDir[0]);
0067   T theta = std::acos(unitDir[2]);
0068   return {
0069       phi,
0070       theta,
0071   };
0072 }
0073
0074 /// Construct the first curvilinear unit vector `U` for the given direction.
0075 ///
0076 /// @param direction is the input direction vector
0077 /// @returns a normalized vector in the x-y plane orthogonal to the direction.
0078 ///
0079 /// The special case of the direction vector pointing along the z-axis is
0080 /// handled by forcing the unit vector to along the x-axis.
0081 template <typename InputVector>
0082 inline auto makeCurvilinearUnitU(
0083     const Eigen::MatrixBase<InputVector>& direction) {
0084   EIGEN_STATIC_ASSERT_FIXED_SIZE(InputVector);
0085   EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputVector);
0086   static_assert(3 <= InputVector::RowsAtCompileTime,
0087                 "Direction vector must be at least three-dimensional.");
0088
0089   using OutputVector = typename InputVector::PlainObject;
0090   using OutputScalar = typename InputVector::Scalar;
0091
0092   OutputVector unitU = OutputVector::Zero();
0093   // explicit version of U = Z x T
0094   unitU[0] = -direction[1];
0095   unitU[1] = direction[0];
0096   const auto scale = unitU.template head<2>().norm();
0097   // if the absolute scale is tiny, the initial direction vector is aligned with
0098   // the z-axis. the ZxT product is ill-defined since any vector in the x-y
0099   // plane would be orthogonal to the direction. fix the U unit vector along the
0100   // x-axis to avoid this numerical instability.
0101   if (scale < (16 * std::numeric_limits<OutputScalar>::epsilon())) {
0102     unitU[0] = 1;
0103     unitU[1] = 0;
0104   } else {
0105     unitU.template head<2>() /= scale;
0106   }
0107   return unitU;
0108 }
0109
0110 /// Construct the curvilinear unit vectors `U` and `V` for the given direction.
0111 ///
0112 /// @param direction is the input direction vector
0113 /// @returns normalized unit vectors `U` and `V` orthogonal to the direction.
0114 ///
0115 /// With `T` the normalized input direction, the three vectors `U`, `V`, and
0116 /// `T` form an orthonormal basis set, i.e. they satisfy
0117 ///
0118 ///     U x V = T
0119 ///     V x T = U
0120 ///     T x U = V
0121 ///
0122 /// with the additional condition that `U` is located in the global x-y plane.
0123 template <typename InputVector>
0124 inline auto makeCurvilinearUnitVectors(
0125     const Eigen::MatrixBase<InputVector>& direction) {
0126   EIGEN_STATIC_ASSERT_FIXED_SIZE(InputVector);
0127   EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputVector);
0128   static_assert(3 <= InputVector::RowsAtCompileTime,
0129                 "Direction vector must be at least three-dimensional.");
0130
0131   using OutputVector = typename InputVector::PlainObject;
0132
0133   std::pair<OutputVector, OutputVector> unitVectors;
0134   unitVectors.first = makeCurvilinearUnitU(direction);
0135   unitVectors.second = direction.cross(unitVectors.first);
0136   unitVectors.second.normalize();
0137   return unitVectors;
0138 }
0139
0140 }  // namespace Acts