UnitVectors.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file UnitVectors.hpp

// This file is part of the Acts project.
//
// Copyright (C) 2020 CERN for the benefit of the Acts project
//
// This Source Code Form is subject to the terms of the Mozilla Public
// License, v. 2.0. If a copy of the MPL was not distributed with this
// file, You can obtain one at http://mozilla.org/MPL/2.0/.

#pragma once

#include <cmath>
#include <limits>
#include <utility>

#include "Acts/Utilities/Definitions.hpp"

namespace Acts {

template <typename T>
inline ActsVector<T, 3> makeDirectionUnitFromPhiEta(T phi, T eta) {
  const auto coshEta = std::cosh(eta);
  const auto tanhEta = std::tanh(eta);
  return {
      std::cos(phi) / coshEta,
      std::sin(phi) / coshEta,
      std::tanh(eta),
  };
}

template <typename T>
inline ActsVector<T, 3> makeDirectionUnitFromPhiTheta(T phi, T theta) {
  const auto cosTheta = std::cos(theta);
  const auto sinTheta = std::sin(theta);
  return {
      std::cos(phi) * sinTheta,
      std::sin(phi) * sinTheta,
      cosTheta,
  };
}

template <typename InputVector>
inline auto makeCurvilinearUnitU(
    const Eigen::MatrixBase<InputVector>& direction) {
  EIGEN_STATIC_ASSERT_FIXED_SIZE(InputVector);
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputVector);
  static_assert(3 <= InputVector::RowsAtCompileTime,
                "Direction vector must be at least three-dimensional.");

  using OutputVector = typename InputVector::PlainObject;
  using OutputScalar = typename InputVector::Scalar;

  OutputVector unitU = OutputVector::Zero();
  // explicit version of U = Z x T
  unitU[0] = -direction[1];
  unitU[1] = direction[0];
  const auto scale = unitU.template head<2>().norm();
  // if the absolute scale is tiny, the initial direction vector is aligned with
  // the z-axis. the ZxT product is ill-defined since any vector in the x-y
  // plane would be orthogonal to the direction. fix the U unit vector along the
  // x-axis to avoid this numerical instability.
  if (scale < (16 * std::numeric_limits<OutputScalar>::epsilon())) {
    unitU[0] = 1;
    unitU[1] = 0;
  } else {
    unitU.template head<2>() /= scale;
  }
  return unitU;
}

template <typename InputVector>
inline auto makeCurvilinearUnitVectors(
    const Eigen::MatrixBase<InputVector>& direction) {
  EIGEN_STATIC_ASSERT_FIXED_SIZE(InputVector);
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(InputVector);
  static_assert(3 <= InputVector::RowsAtCompileTime,
                "Direction vector must be at least three-dimensional.");

  using OutputVector = typename InputVector::PlainObject;

  std::pair<OutputVector, OutputVector> unitVectors;
  unitVectors.first = makeCurvilinearUnitU(direction);
  unitVectors.second = direction.cross(unitVectors.first);
  unitVectors.second.normalize();
  return unitVectors;
}

}  // namespace Acts