MultiComponentTrackParameters.hpp

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

// This file is part of the Acts project.
//
// Copyright (C) 2021 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/EventData/GenericBoundTrackParameters.hpp"
#include "Acts/EventData/TrackParameters.hpp"
#include "Acts/Surfaces/PlaneSurface.hpp"
#include "Acts/Surfaces/Surface.hpp"

#include <cmath>
#include <memory>
#include <type_traits>
#include <utility>

namespace Acts {

class MultiComponentBoundTrackParameters {
 public:
  using Parameters = BoundTrackParameters;
  using ParticleHypothesis = Parameters::ParticleHypothesis;
  using Scalar = typename Parameters::Scalar;
  using ParametersVector = typename Parameters::ParametersVector;
  using CovarianceMatrix = typename Parameters::CovarianceMatrix;

 private:
  std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
      m_components;
  std::shared_ptr<const Surface> m_surface;

  // TODO use [[no_unique_address]] once we switch to C++20
  ParticleHypothesis m_particleHypothesis;

  template <typename projector_t>
  auto reduce(projector_t&& proj) const {
    using Ret = std::decay_t<decltype(proj(std::declval<Parameters>()))>;

    Ret ret;

    if constexpr (std::is_floating_point_v<Ret>) {
      ret = 0.0;
    } else {
      ret = Ret::Zero();
    }

    for (auto i = 0ul; i < m_components.size(); ++i) {
      const auto [weight, single_pars] = (*this)[i];
      ret += weight * proj(single_pars);
    }

    return ret;
  }

 public:
  template <typename covariance_t>
  MultiComponentBoundTrackParameters(
      std::shared_ptr<const Surface> surface,
      const std::vector<std::tuple<double, BoundVector, covariance_t>>& cmps,
      ParticleHypothesis particleHypothesis)
      : m_surface(std::move(surface)),
        m_particleHypothesis(particleHypothesis) {
    static_assert(
        std::is_same_v<BoundSquareMatrix, covariance_t> ||
        std::is_same_v<std::optional<BoundSquareMatrix>, covariance_t>);
    if constexpr (std::is_same_v<BoundSquareMatrix, covariance_t>) {
      for (const auto& [weight, params, cov] : cmps) {
        m_components.push_back({weight, params, cov});
      }
    } else {
      m_components = cmps;
    }
  }

  MultiComponentBoundTrackParameters(std::shared_ptr<const Surface> surface,
                                     const BoundVector& params,
                                     std::optional<BoundSquareMatrix> cov,
                                     ParticleHypothesis particleHypothesis)
      : m_surface(std::move(surface)),
        m_particleHypothesis(particleHypothesis) {
    m_components.push_back({1., params, std::move(cov)});
  }

  MultiComponentBoundTrackParameters() = delete;
  MultiComponentBoundTrackParameters(
      const MultiComponentBoundTrackParameters&) = default;
  MultiComponentBoundTrackParameters(MultiComponentBoundTrackParameters&&) =
      default;
  ~MultiComponentBoundTrackParameters() = default;
  MultiComponentBoundTrackParameters& operator=(
      const MultiComponentBoundTrackParameters&) = default;
  MultiComponentBoundTrackParameters& operator=(
      MultiComponentBoundTrackParameters&&) = default;

  const auto& components() const { return m_components; }

  const Surface& referenceSurface() const { return *m_surface; }

  std::pair<double, Parameters> operator[](std::size_t i) const {
    return std::make_pair(
        std::get<double>(m_components[i]),
        Parameters(m_surface, std::get<BoundVector>(m_components[i]),
                   std::get<std::optional<BoundSquareMatrix>>(m_components[i]),
                   m_particleHypothesis));
  }

  ParametersVector parameters() const {
    return reduce([](const Parameters& p) { return p.parameters(); });
  }

  std::optional<CovarianceMatrix> covariance() const {
    const auto ret = reduce([](const Parameters& p) {
      return p.covariance() ? *p.covariance() : CovarianceMatrix::Zero();
    });

    if (ret == CovarianceMatrix::Zero()) {
      return std::nullopt;
    } else {
      return ret;
    }
  }

  template <BoundIndices kIndex>
  Scalar get() const {
    return reduce([&](const Parameters& p) { return p.get<kIndex>(); });
  }

  Vector4 fourPosition(const GeometryContext& geoCtx) const {
    return reduce([&](const Parameters& p) { return p.fourPosition(geoCtx); });
  }

  Vector3 position(const GeometryContext& geoCtx) const {
    return reduce([&](const Parameters& p) { return p.position(geoCtx); });
  }

  Scalar time() const {
    return reduce([](const Parameters& p) { return p.time(); });
  }

  Vector3 direction() const {
    return reduce([](const Parameters& p) { return p.direction(); })
        .normalized();
  }

  Scalar phi() const { return VectorHelpers::phi(direction()); }

  Scalar theta() const { return VectorHelpers::theta(direction()); }

  Scalar qOverP() const { return get<eBoundQOverP>(); }

  Scalar absoluteMomentum() const {
    return reduce([](const Parameters& p) { return p.absoluteMomentum(); });
  }

  Scalar transverseMomentum() const {
    return reduce([](const Parameters& p) { return p.transverseMomentum(); });
  }

  Vector3 momentum() const {
    return reduce([](const Parameters& p) { return p.momentum(); });
  }

  Scalar charge() const {
    return reduce([](const Parameters& p) { return p.charge(); });
  }

  const ParticleHypothesis& particleHypothesis() const {
    return m_particleHypothesis;
  }
};

class MultiComponentCurvilinearTrackParameters
    : public MultiComponentBoundTrackParameters {
  using covariance_t = BoundSquareMatrix;

 public:
  using ConstructionTuple = std::tuple<double, Acts::Vector4, Acts::Vector3,
                                       ActsScalar, covariance_t>;

 private:
  using Base = MultiComponentBoundTrackParameters;

  using BaseConstructionTuple =
      std::tuple<std::shared_ptr<Acts::Surface>,
                 std::vector<std::tuple<double, BoundVector, covariance_t>>>;

  static BaseConstructionTuple construct(
      const std::vector<ConstructionTuple>& curvi) {
    // TODO where to get a geometry context here
    Acts::GeometryContext gctx{};

    // Construct and average surface
    Acts::Vector3 avgPos = Acts::Vector3::Zero();
    Acts::Vector3 avgDir = Acts::Vector3::Zero();
    for (const auto& [w, pos4, dir, qop, cov] : curvi) {
      avgPos += w * pos4.template segment<3>(0);
      avgDir += w * dir;
    }

    auto s = Surface::makeShared<PlaneSurface>(avgPos, avgDir);

    std::vector<std::tuple<double, BoundVector, covariance_t>> bound;
    bound.reserve(curvi.size());

    // Project the position onto the surface, keep everything else as is
    for (const auto& [w, pos4, dir, qop, cov] : curvi) {
      Vector3 newPos =
          s->intersect(gctx, pos4.template segment<3>(eFreePos0), dir, false)
              .closest()
              .position();

      BoundVector bv =
          detail::transformFreeToCurvilinearParameters(pos4[eTime], dir, qop);

      // Because of the projection this should never fail
      bv.template segment<2>(eBoundLoc0) =
          *(s->globalToLocal(gctx, newPos, dir));

      bound.emplace_back(w, bv, cov);
    }

    return {s, bound};
  }

  MultiComponentCurvilinearTrackParameters(
      const BaseConstructionTuple& t, ParticleHypothesis particleHypothesis)
      : Base(std::get<0>(t), std::get<1>(t), particleHypothesis) {}

 public:
  MultiComponentCurvilinearTrackParameters(
      const std::vector<ConstructionTuple>& cmps,
      ParticleHypothesis particleHypothesis)
      : MultiComponentCurvilinearTrackParameters(construct(cmps),
                                                 particleHypothesis) {}
};

}  // namespace Acts