CalculateResiduals.hpp

Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CalculateResiduals.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 "Acts/Definitions/TrackParametrization.hpp"
#include "Acts/Utilities/detail/periodic.hpp"

#include <cassert>

#include <Eigen/Core>

namespace Acts {
namespace detail {

template <typename index_container_t, typename measured_t, typename residuals_t>
inline void calculateResiduals(BoundIndices size,
                               const index_container_t& indices,
                               const BoundVector& reference,
                               const Eigen::MatrixBase<measured_t>& measured,
                               Eigen::MatrixBase<residuals_t>& residuals) {
  using OutputScalar = typename residuals_t::Scalar;

  EIGEN_STATIC_ASSERT_VECTOR_ONLY(measured_t);
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(residuals_t);
  assert((size <= eBoundSize) and "Measured subspace is too large");
  assert((size <= measured.size()) and "Inconsistent measured size");
  assert((size <= residuals.size()) and "Inconsistent residuals size");

  for (size_t i = 0; i < size; ++i) {
    size_t fullIndex = indices[i];
    // this is neither elegant nor smart but it is the simplest solution.
    //
    // only phi must be handled specially here. the theta limits can only be
    // correctly handled if phi is updated, too. since we can not ensure that
    // both are available, it is probably less error-prone to treat theta as a
    // regular, unrestricted parameter.
    if (fullIndex == eBoundPhi) {
      residuals[i] = difference_periodic<OutputScalar>(
          measured[i], reference[fullIndex],
          static_cast<OutputScalar>(2 * M_PI));
    } else {
      residuals[i] = measured[i] - reference[fullIndex];
    }
  }
}

template <typename index_container_t, typename measured_t, typename residuals_t>
inline void calculateResiduals(FreeIndices size,
                               const index_container_t& indices,
                               const FreeVector& reference,
                               const Eigen::MatrixBase<measured_t>& measured,
                               Eigen::MatrixBase<residuals_t>& residuals) {
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(measured_t);
  EIGEN_STATIC_ASSERT_VECTOR_ONLY(residuals_t);
  assert((size <= eFreeSize) and "Measured subspace is too large");
  assert((size <= measured.size()) and "Inconsistent measured size");
  assert((size <= residuals.size()) and "Inconsistent residuals size");

  for (size_t i = 0; i < size; ++i) {
    // all free parameters are unrestricted. no need to call parameter traits
    residuals[i] = measured[i] - reference[indices[i]];
  }
}

}  // namespace detail
}  // namespace Acts