Helpers.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2016-2018 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/.

// AlgebraHelper.h, Acts project

#pragma once

#include "Acts/Utilities/BinningType.hpp"
#include "Acts/Utilities/Definitions.hpp"
#include "Acts/Utilities/ParameterDefinitions.hpp"
#include "Acts/Utilities/TypeTraits.hpp"

#include <bitset>
#include <cassert>
#include <cmath>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <memory>
#include <string>
#include <vector>

#define ACTS_CHECK_BIT(value, mask) ((value & mask) == mask)

namespace Acts {

namespace VectorHelpers {
namespace detail {
template <class T>
using phi_method_t = decltype(std::declval<const T>().phi());

template <class T>
using has_phi_method = Concepts ::is_detected<phi_method_t, T>;

}  // namespace detail

template <typename Derived>
double phi(const Eigen::MatrixBase<Derived>& v) noexcept {
  constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
  if constexpr (rows != -1) {
    // static size, do compile time check
    static_assert(rows >= 2,
                  "Phi function not valid for vectors not at least 2D");
  } else {
    // dynamic size
    if (v.rows() < 2) {
      std::cerr << "Phi function not valid for vectors not at least 2D"
                << std::endl;
      std::abort();
    }
  }

  return std::atan2(v[1], v[0]);
}

template <typename T,
          std::enable_if_t<detail::has_phi_method<T>::value, int> = 0>
double phi(const T& v) noexcept {
  return v.phi();
}

template <typename Derived>
double perp(const Eigen::MatrixBase<Derived>& v) noexcept {
  constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
  if constexpr (rows != -1) {
    // static size, do compile time check
    static_assert(rows >= 2,
                  "Perp function not valid for vectors not at least 2D");
  } else {
    // dynamic size
    if (v.rows() < 2) {
      std::cerr << "Perp function not valid for vectors not at least 2D"
                << std::endl;
      std::abort();
    }
  }
  return std::sqrt(v[0] * v[0] + v[1] * v[1]);
}

template <typename Derived>
double theta(const Eigen::MatrixBase<Derived>& v) noexcept {
  constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
  if constexpr (rows != -1) {
    // static size, do compile time check
    static_assert(rows >= 3, "Theta function not valid for non-3D vectors.");
  } else {
    // dynamic size
    if (v.rows() < 3) {
      std::cerr << "Theta function not valid for non-3D vectors." << std::endl;
      std::abort();
    }
  }

  return std::atan2(std::sqrt(v[0] * v[0] + v[1] * v[1]), v[2]);
}

template <typename Derived>
double eta(const Eigen::MatrixBase<Derived>& v) noexcept {
  constexpr int rows = Eigen::MatrixBase<Derived>::RowsAtCompileTime;
  if constexpr (rows != -1) {
    // static size, do compile time check
    static_assert(rows >= 3, "Eta function not valid for non-3D vectors.");
  } else {
    // dynamic size
    if (v.rows() < 3) {
      std::cerr << "Eta function not valid for non-3D vectors." << std::endl;
      std::abort();
    }
  }

  return std::atanh(v[2] / v.norm());
}

inline double cast(const Vector3D& position, BinningValue bval) {
  switch (bval) {
    case binX:
      return position[0];
    case binY:
      return position[1];
    case binZ:
      return position[2];
    case binR:
      return perp(position);
    case binPhi:
      return phi(position);
    case binRPhi:
      return perp(position) * phi(position);
    case binH:
      return theta(position);
    case binEta:
      return eta(position);
    case binMag:
      return position.norm();
    default:
      assert(false and "Invalid BinningValue enum value");
      return std::numeric_limits<double>::quiet_NaN();
  }
}

inline ActsMatrixD<3, 3> cross(const ActsMatrixD<3, 3>& m, const Vector3D& v) {
  ActsMatrixD<3, 3> r;
  r.col(0) = m.col(0).cross(v);
  r.col(1) = m.col(1).cross(v);
  r.col(2) = m.col(2).cross(v);

  return r;
}

inline auto position(const Vector4D& pos4) {
  return pos4.segment<3>(ePos0);
}

inline auto position(const FreeVector& params) {
  return params.segment<3>(eFreePos0);
}

template <typename vector3_t>
inline auto makeVector4(const Eigen::MatrixBase<vector3_t>& vec3,
                        typename vector3_t::Scalar w)
    -> Eigen::Matrix<typename vector3_t::Scalar, 4, 1> {
  EIGEN_STATIC_ASSERT_VECTOR_SPECIFIC_SIZE(vector3_t, 3);

  Eigen::Matrix<typename vector3_t::Scalar, 4, 1> vec4;
  vec4[ePos0] = vec3[ePos0];
  vec4[ePos1] = vec3[ePos1];
  vec4[ePos2] = vec3[ePos2];
  vec4[eTime] = w;
  return vec4;
}

}  // namespace VectorHelpers

namespace detail {

inline double roundWithPrecision(double val, int precision) {
  if (val < 0 && std::abs(val) * std::pow(10, precision) < 1.) {
    return -val;
  }
  return val;
}
}  // namespace detail

inline std::string toString(const ActsMatrixXd& matrix, int precision = 4,
                            const std::string& offset = "") {
  std::ostringstream sout;

  sout << std::setiosflags(std::ios::fixed) << std::setprecision(precision);
  if (matrix.cols() == 1) {
    sout << "(";
    for (int i = 0; i < matrix.rows(); ++i) {
      double val = detail::roundWithPrecision(matrix(i, 0), precision);
      sout << val;
      if (i != matrix.rows() - 1) {
        sout << ", ";
      }
    }
    sout << ")";
  } else {
    for (int i = 0; i < matrix.rows(); ++i) {
      for (int j = 0; j < matrix.cols(); ++j) {
        if (j == 0) {
          sout << "(";
        }
        double val = detail::roundWithPrecision(matrix(i, j), precision);
        sout << val;
        if (j == matrix.cols() - 1) {
          sout << ")";
        } else {
          sout << ", ";
        }
      }
      if (i != matrix.rows() -
                   1) {  // make the end line and the offset in the next line
        sout << std::endl;
        sout << offset;
      }
    }
  }
  return sout.str();
}

inline std::string toString(const Acts::Translation3D& translation,
                            int precision = 4) {
  Acts::Vector3D trans;
  trans[0] = translation.x();
  trans[1] = translation.y();
  trans[2] = translation.z();
  return toString(trans, precision);
}

inline std::string toString(const Acts::Transform3D& transform,
                            int precision = 4, const std::string& offset = "") {
  std::ostringstream sout;
  sout << "Translation : " << toString(transform.translation(), precision)
       << std::endl;
  std::string rotationOffset = offset + "              ";
  sout << offset << "Rotation    : "
       << toString(transform.rotation(), precision + 2, rotationOffset);
  return sout.str();
}

template <typename T>
std::vector<T*> unpack_shared_vector(
    const std::vector<std::shared_ptr<T>>& items) {
  std::vector<T*> rawPtrs;
  rawPtrs.reserve(items.size());
  for (const std::shared_ptr<T>& item : items) {
    rawPtrs.push_back(item.get());
  }
  return rawPtrs;
}

template <typename T>
std::vector<const T*> unpack_shared_vector(
    const std::vector<std::shared_ptr<const T>>& items) {
  std::vector<const T*> rawPtrs;
  rawPtrs.reserve(items.size());
  for (const std::shared_ptr<const T>& item : items) {
    rawPtrs.push_back(item.get());
  }
  return rawPtrs;
}

template <template <size_t> class Callable, size_t N, size_t NMAX,
          typename... Args>
decltype(Callable<N>::invoke(std::declval<Args>()...)) template_switch(
    size_t v, Args&&... args) {
  if (v == N) {
    return Callable<N>::invoke(std::forward<Args>(args)...);
  }
  if constexpr (N < NMAX) {
    return template_switch<Callable, N + 1, NMAX>(v,
                                                  std::forward<Args>(args)...);
  }
  std::cerr << "template_switch<Fn, " << N << ", " << NMAX << ">(v=" << v
            << ") is not valid (v > NMAX)" << std::endl;
  std::abort();
}

template <typename MatrixType>
MatrixType bitsetToMatrix(const std::bitset<MatrixType::RowsAtCompileTime *
                                            MatrixType::ColsAtCompileTime>
                              bs) {
  constexpr int rows = MatrixType::RowsAtCompileTime;
  constexpr int cols = MatrixType::ColsAtCompileTime;

  static_assert(rows != -1 && cols != -1,
                "bitsetToMatrix does not support dynamic matrices");

  MatrixType m;
  auto* p = m.data();
  for (size_t i = 0; i < rows * cols; i++) {
    p[i] = bs[rows * cols - 1 - i];
  }
  return m;
}

template <typename Derived>
auto matrixToBitset(const Eigen::PlainObjectBase<Derived>& m) {
  using MatrixType = Eigen::PlainObjectBase<Derived>;
  constexpr size_t rows = MatrixType::RowsAtCompileTime;
  constexpr size_t cols = MatrixType::ColsAtCompileTime;

  std::bitset<rows * cols> res;

  auto* p = m.data();
  for (size_t i = 0; i < rows * cols; i++) {
    res[rows * cols - 1 - i] = p[i];
  }

  return res;
}

}  // namespace Acts