Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:10:48

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/EventData/GenericBoundTrackParameters.hpp"
0012 #include "Acts/EventData/TrackParameters.hpp"
0013 #include "Acts/Surfaces/PlaneSurface.hpp"
0014 #include "Acts/Surfaces/Surface.hpp"
0015 
0016 #include <cmath>
0017 #include <memory>
0018 #include <type_traits>
0019 #include <utility>
0020 
0021 namespace Acts {
0022 
0023 /// This class is only a light wrapper around a surface and a vector of
0024 /// parameters. Its main purpose is to provide many constructors for the
0025 /// underlying vector. Most accessors are generated from the
0026 /// BoundTrackParameters equivalent and thus may be expensive
0027 /// @note This class holds shared ownership on its reference surface.
0028 /// @note The accessors for parameters, covariance, position, etc.
0029 /// are the weighted means of the components.
0030 /// @note If all covariances are zero, the accessor for the total
0031 /// covariance does return std::nullopt;
0032 /// TODO Add constructor from range and projector maybe?
0033 class MultiComponentBoundTrackParameters {
0034  public:
0035   using Parameters = BoundTrackParameters;
0036   using ParticleHypothesis = Parameters::ParticleHypothesis;
0037   using ParametersVector = typename Parameters::ParametersVector;
0038   using CovarianceMatrix = typename Parameters::CovarianceMatrix;
0039 
0040  private:
0041   std::vector<std::tuple<double, BoundVector, std::optional<BoundSquareMatrix>>>
0042       m_components;
0043   std::shared_ptr<const Surface> m_surface;
0044 
0045   // TODO use [[no_unique_address]] once we switch to C++20
0046   ParticleHypothesis m_particleHypothesis;
0047 
0048   /// Helper function to reduce the component vector to a single representation
0049   template <typename projector_t>
0050   auto reduce(projector_t&& proj) const {
0051     using Ret = std::decay_t<decltype(proj(std::declval<Parameters>()))>;
0052 
0053     Ret ret;
0054 
0055     if constexpr (std::is_floating_point_v<Ret>) {
0056       ret = 0.0;
0057     } else {
0058       ret = Ret::Zero();
0059     }
0060 
0061     for (auto i = 0ul; i < m_components.size(); ++i) {
0062       const auto [weight, single_pars] = (*this)[i];
0063       ret += weight * proj(single_pars);
0064     }
0065 
0066     return ret;
0067   }
0068 
0069  public:
0070   /// Construct from multiple components
0071   template <typename covariance_t>
0072   MultiComponentBoundTrackParameters(
0073       std::shared_ptr<const Surface> surface,
0074       const std::vector<std::tuple<double, BoundVector, covariance_t>>& cmps,
0075       ParticleHypothesis particleHypothesis)
0076       : m_surface(std::move(surface)),
0077         m_particleHypothesis(particleHypothesis) {
0078     static_assert(
0079         std::is_same_v<BoundSquareMatrix, covariance_t> ||
0080         std::is_same_v<std::optional<BoundSquareMatrix>, covariance_t>);
0081     if constexpr (std::is_same_v<BoundSquareMatrix, covariance_t>) {
0082       for (const auto& [weight, params, cov] : cmps) {
0083         m_components.push_back({weight, params, cov});
0084       }
0085     } else {
0086       m_components = cmps;
0087     }
0088   }
0089 
0090   /// Construct from a parameters vector on the surface and particle charge.
0091   ///
0092   /// @param surface Reference surface the parameters are defined on
0093   /// @param params Bound parameters vector
0094   /// @param particleHypothesis Particle hypothesis for these parameters
0095   /// @param cov Bound parameters covariance matrix
0096   ///
0097   /// In principle, only the charge magnitude is needed her to allow
0098   /// unambiguous extraction of the absolute momentum. The particle charge is
0099   /// required as an input here to be consistent with the other constructors
0100   /// below that that also take the charge as an input. The charge sign is
0101   /// only used in debug builds to check for consistency with the q/p
0102   /// parameter.
0103   MultiComponentBoundTrackParameters(std::shared_ptr<const Surface> surface,
0104                                      const BoundVector& params,
0105                                      std::optional<BoundSquareMatrix> cov,
0106                                      ParticleHypothesis particleHypothesis)
0107       : m_surface(std::move(surface)),
0108         m_particleHypothesis(particleHypothesis) {
0109     m_components.push_back({1., params, std::move(cov)});
0110   }
0111 
0112   /// Parameters are not default constructible due to the charge type.
0113   MultiComponentBoundTrackParameters() = delete;
0114   MultiComponentBoundTrackParameters(
0115       const MultiComponentBoundTrackParameters&) = default;
0116   MultiComponentBoundTrackParameters(MultiComponentBoundTrackParameters&&) =
0117       default;
0118   ~MultiComponentBoundTrackParameters() = default;
0119   MultiComponentBoundTrackParameters& operator=(
0120       const MultiComponentBoundTrackParameters&) = default;
0121   MultiComponentBoundTrackParameters& operator=(
0122       MultiComponentBoundTrackParameters&&) = default;
0123 
0124   /// Access the parameters
0125   const auto& components() const { return m_components; }
0126 
0127   /// Reference surface onto which the parameters are bound.
0128   const Surface& referenceSurface() const { return *m_surface; }
0129 
0130   /// Get the weight and a GenericBoundTrackParameters object for one component
0131   std::pair<double, Parameters> operator[](std::size_t i) const {
0132     return {
0133         std::get<double>(m_components[i]),
0134         Parameters(m_surface, std::get<BoundVector>(m_components[i]),
0135                    std::get<std::optional<BoundSquareMatrix>>(m_components[i]),
0136                    m_particleHypothesis)};
0137   }
0138 
0139   /// Parameters vector.
0140   ParametersVector parameters() const {
0141     return reduce([](const Parameters& p) { return p.parameters(); });
0142   }
0143 
0144   /// Optional covariance matrix.
0145   std::optional<CovarianceMatrix> covariance() const {
0146     const auto ret = reduce([](const Parameters& p) {
0147       return p.covariance() ? *p.covariance() : CovarianceMatrix::Zero();
0148     });
0149 
0150     if (ret == CovarianceMatrix::Zero()) {
0151       return std::nullopt;
0152     } else {
0153       return ret;
0154     }
0155   }
0156 
0157   /// Access a single parameter value identified by its index.
0158   ///
0159   /// @tparam kIndex Track parameter index
0160   template <BoundIndices kIndex>
0161   double get() const {
0162     return reduce([&](const Parameters& p) { return p.get<kIndex>(); });
0163   }
0164 
0165   /// Space-time position four-vector.
0166   ///
0167   /// @param[in] geoCtx Geometry context for the local-to-global
0168   /// transformation
0169   Vector4 fourPosition(const GeometryContext& geoCtx) const {
0170     return reduce([&](const Parameters& p) { return p.fourPosition(geoCtx); });
0171   }
0172 
0173   /// Spatial position three-vector.
0174   ///
0175   /// @param[in] geoCtx Geometry context for the local-to-global
0176   /// transformation
0177   Vector3 position(const GeometryContext& geoCtx) const {
0178     return reduce([&](const Parameters& p) { return p.position(geoCtx); });
0179   }
0180 
0181   /// Time coordinate.
0182   double time() const {
0183     return reduce([](const Parameters& p) { return p.time(); });
0184   }
0185 
0186   /// Unit direction three-vector, i.e. the normalized momentum
0187   /// three-vector.
0188   Vector3 direction() const {
0189     return reduce([](const Parameters& p) { return p.direction(); })
0190         .normalized();
0191   }
0192 
0193   /// Phi direction.
0194   double phi() const { return VectorHelpers::phi(direction()); }
0195 
0196   /// Theta direction.
0197   double theta() const { return VectorHelpers::theta(direction()); }
0198 
0199   /// Charge over momentum.
0200   double qOverP() const { return get<eBoundQOverP>(); }
0201 
0202   /// Absolute momentum.
0203   double absoluteMomentum() const {
0204     return reduce([](const Parameters& p) { return p.absoluteMomentum(); });
0205   }
0206 
0207   /// Transverse momentum.
0208   double transverseMomentum() const {
0209     return reduce([](const Parameters& p) { return p.transverseMomentum(); });
0210   }
0211 
0212   /// Momentum three-vector.
0213   Vector3 momentum() const {
0214     return reduce([](const Parameters& p) { return p.momentum(); });
0215   }
0216 
0217   /// Particle electric charge.
0218   double charge() const {
0219     return reduce([](const Parameters& p) { return p.charge(); });
0220   }
0221 
0222   /// Particle hypothesis.
0223   const ParticleHypothesis& particleHypothesis() const {
0224     return m_particleHypothesis;
0225   }
0226 };
0227 
0228 /// This class mimics the behaviour of the curvilinear parameters for ordinary
0229 /// track parameters. To adopt this concept, a "common surface" is constructed,
0230 /// and all parameters are projected onto this surface. The use of this is
0231 /// questionable, and if the result is reasonable depends largely on the initial
0232 /// multi component state. However, the propagator infrastructure forces the
0233 /// existence of this type
0234 /// @tparam charge_t Helper type to interpret the particle charge/momentum
0235 class MultiComponentCurvilinearTrackParameters
0236     : public MultiComponentBoundTrackParameters {
0237   using covariance_t = BoundSquareMatrix;
0238 
0239  public:
0240   using ConstructionTuple =
0241       std::tuple<double, Acts::Vector4, Acts::Vector3, double, covariance_t>;
0242 
0243  private:
0244   using Base = MultiComponentBoundTrackParameters;
0245 
0246   using BaseConstructionTuple =
0247       std::tuple<std::shared_ptr<Acts::Surface>,
0248                  std::vector<std::tuple<double, BoundVector, covariance_t>>>;
0249 
0250   /// We need this helper function in order to construct the base class properly
0251   static BaseConstructionTuple construct(
0252       const std::vector<ConstructionTuple>& curvi) {
0253     // TODO where to get a geometry context here
0254     Acts::GeometryContext gctx{};
0255 
0256     // Construct and average surface
0257     Acts::Vector3 avgPos = Acts::Vector3::Zero();
0258     Acts::Vector3 avgDir = Acts::Vector3::Zero();
0259     for (const auto& [w, pos4, dir, qop, cov] : curvi) {
0260       avgPos += w * pos4.template segment<3>(0);
0261       avgDir += w * dir;
0262     }
0263 
0264     auto s = CurvilinearSurface(avgPos, avgDir).planeSurface();
0265 
0266     std::vector<std::tuple<double, BoundVector, covariance_t>> bound;
0267     bound.reserve(curvi.size());
0268 
0269     // Project the position onto the surface, keep everything else as is
0270     for (const auto& [w, pos4, dir, qop, cov] : curvi) {
0271       Vector3 newPos = s->intersect(gctx, pos4.template segment<3>(eFreePos0),
0272                                     dir, BoundaryTolerance::Infinite())
0273                            .closest()
0274                            .position();
0275 
0276       BoundVector bv =
0277           transformFreeToCurvilinearParameters(pos4[eTime], dir, qop);
0278 
0279       // Because of the projection this should never fail
0280       bv.template segment<2>(eBoundLoc0) =
0281           *(s->globalToLocal(gctx, newPos, dir));
0282 
0283       bound.emplace_back(w, bv, cov);
0284     }
0285 
0286     return {s, bound};
0287   }
0288 
0289   /// Private constructor from a tuple
0290   MultiComponentCurvilinearTrackParameters(
0291       const BaseConstructionTuple& t, ParticleHypothesis particleHypothesis)
0292       : Base(std::get<0>(t), std::get<1>(t), particleHypothesis) {}
0293 
0294  public:
0295   MultiComponentCurvilinearTrackParameters(
0296       const std::vector<ConstructionTuple>& cmps,
0297       ParticleHypothesis particleHypothesis)
0298       : MultiComponentCurvilinearTrackParameters(construct(cmps),
0299                                                  particleHypothesis) {}
0300 };
0301 
0302 }  // namespace Acts