Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-10-31 08:16:04

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/Definitions/Algebra.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/Types.hpp"
0014 #include "Acts/Utilities/AlgebraHelpers.hpp"
0015 #include "Acts/Utilities/Enumerate.hpp"
0016 
0017 #include <cstddef>
0018 #include <ranges>
0019 
0020 #include <boost/container/static_vector.hpp>
0021 
0022 namespace Acts {
0023 
0024 /// @brief Check subspace indices for consistency
0025 ///
0026 /// Indices must be unique and within the full size of the subspace
0027 ///
0028 /// @tparam index_range_t the type of the container of indices
0029 ///
0030 /// @param indexRange the range of indices
0031 /// @param fullSize the full size of the subspace
0032 /// @param subspaceSize the size of the subspace
0033 ///
0034 /// @return true if the indices are consistent
0035 template <std::ranges::sized_range index_range_t>
0036 inline static bool checkSubspaceIndices(const index_range_t& indexRange,
0037                                         std::size_t fullSize,
0038                                         std::size_t subspaceSize) {
0039   if (subspaceSize > fullSize) {
0040     return false;
0041   }
0042   if (static_cast<std::size_t>(indexRange.size()) != subspaceSize) {
0043     return false;
0044   }
0045   for (auto it = indexRange.begin(); it != indexRange.end();) {
0046     auto index = *it;
0047     if (index >= fullSize) {
0048       return false;
0049     }
0050     ++it;
0051     if (std::find(it, indexRange.end(), index) != indexRange.end()) {
0052       return false;
0053     }
0054   }
0055   return true;
0056 }
0057 
0058 /// Serialize subspace indices to a single 64 bit integer
0059 ///
0060 /// @tparam FullSize the full size of the subspace
0061 ///
0062 /// @param indices the subspace indices
0063 ///
0064 /// @return the serialized subspace indices
0065 template <std::size_t FullSize>
0066 inline static SerializedSubspaceIndices serializeSubspaceIndices(
0067     const SubspaceIndices<FullSize>& indices)
0068   requires(FullSize <= 8)
0069 {
0070   SerializedSubspaceIndices result = 0;
0071   for (std::size_t i = 0; i < FullSize; ++i) {
0072     result |= static_cast<SerializedSubspaceIndices>(indices[i] & 0xFF)
0073               << (i * 8);
0074   }
0075   return result;
0076 }
0077 
0078 /// Deserialize subspace indices from a single 64 bit integer
0079 ///
0080 /// @tparam FullSize the full size of the subspace
0081 ///
0082 /// @param serialized the serialized subspace indices
0083 ///
0084 /// @return the subspace indices
0085 template <std::size_t FullSize>
0086 inline static SubspaceIndices<FullSize> deserializeSubspaceIndices(
0087     SerializedSubspaceIndices serialized)
0088   requires(FullSize <= 8)
0089 {
0090   SubspaceIndices<FullSize> result;
0091   for (std::size_t i = 0; i < FullSize; ++i) {
0092     result[i] = static_cast<std::uint8_t>((serialized >> (i * 8)) & 0xFF);
0093   }
0094   return result;
0095 }
0096 
0097 namespace detail {
0098 
0099 /// Helper base class for subspace operations which provides common
0100 /// functionality for fixed and variable subspace helpers
0101 ///
0102 /// @tparam Derived the derived type
0103 /// @tparam FullSize the full size of the subspace
0104 template <typename Derived, std::size_t FullSize>
0105 class SubspaceHelperBase {
0106  public:
0107   static constexpr std::size_t kFullSize = FullSize;
0108 
0109   using FullSquareMatrix = ActsSquareMatrix<kFullSize>;
0110 
0111   using size_type = std::size_t;
0112 
0113   bool empty() const { return self().empty(); }
0114   std::size_t size() const { return self().size(); }
0115 
0116   auto operator[](std::size_t i) const { return self()[i]; }
0117 
0118   auto begin() const { return self().begin(); }
0119   auto end() const { return self().end(); }
0120 
0121   bool contains(std::uint8_t index) const {
0122     return std::find(begin(), end(), index) != end();
0123   }
0124   std::size_t indexOf(std::uint8_t index) const {
0125     auto it = std::find(begin(), end(), index);
0126     return it != end() ? std::distance(begin(), it) : kFullSize;
0127   }
0128 
0129   template <typename EigenDerived>
0130   ActsVector<kFullSize> expandVector(
0131       const Eigen::DenseBase<EigenDerived>& vector) const {
0132     ActsVector<kFullSize> result = ActsVector<kFullSize>::Zero();
0133     for (auto [i, index] : enumerate(*this)) {
0134       result(index) = vector(i);
0135     }
0136     return result;
0137   }
0138 
0139   template <typename EigenDerived>
0140   FullSquareMatrix expandMatrix(
0141       const Eigen::DenseBase<EigenDerived>& matrix) const {
0142     FullSquareMatrix result = FullSquareMatrix::Zero();
0143     for (auto [i, indexI] : enumerate(*this)) {
0144       for (auto [j, indexJ] : enumerate(*this)) {
0145         result(indexI, indexJ) = matrix(i, j);
0146       }
0147     }
0148     return result;
0149   }
0150 
0151   FullSquareMatrix fullProjector() const {
0152     FullSquareMatrix result = FullSquareMatrix::Zero();
0153     for (auto [i, index] : enumerate(*this)) {
0154       result(i, index) = 1;
0155     }
0156     return result;
0157   }
0158 
0159   FullSquareMatrix fullExpander() const {
0160     FullSquareMatrix result = FullSquareMatrix::Zero();
0161     for (auto [i, index] : enumerate(*this)) {
0162       result(index, i) = 1;
0163     }
0164     return result;
0165   }
0166 
0167   ProjectorBitset projectorBitset() const
0168     requires(kFullSize <= 8)
0169   {
0170     return matrixToBitset(fullProjector()).to_ullong();
0171   }
0172 
0173  private:
0174   const Derived& self() const { return static_cast<const Derived&>(*this); }
0175 };
0176 
0177 }  // namespace detail
0178 
0179 /// Helper class for variable subspace operations
0180 ///
0181 /// @tparam FullSize the full size of the subspace
0182 /// @tparam index_t the index type
0183 template <std::size_t FullSize, typename index_t = std::uint8_t>
0184 class VariableSubspaceHelper
0185     : public detail::SubspaceHelperBase<
0186           VariableSubspaceHelper<FullSize, index_t>, FullSize> {
0187  public:
0188   /// Size of the full parameter space
0189   static constexpr std::size_t kFullSize = FullSize;
0190 
0191   /// Type alias for index type
0192   using IndexType = index_t;
0193   /// Type alias for container holding subspace indices
0194   using Container = boost::container::static_vector<IndexType, FullSize>;
0195 
0196   /// Construct a variable subspace helper with specified indices
0197   ///
0198   /// @tparam other_index_range_t Type of the index range
0199   /// @param indices Range of indices defining the subspace
0200   template <std::ranges::sized_range other_index_range_t>
0201   explicit VariableSubspaceHelper(const other_index_range_t& indices) {
0202     assert(checkSubspaceIndices(indices, kFullSize, indices.size()) &&
0203            "Invalid indices");
0204     m_indices.resize(indices.size());
0205     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0206                    [](auto index) { return static_cast<IndexType>(index); });
0207   }
0208 
0209   /// Check if the subspace is empty
0210   /// @return True if the subspace contains no indices
0211   bool empty() const { return m_indices.empty(); }
0212   /// Get the size of the subspace
0213   /// @return Number of indices in the subspace
0214   std::size_t size() const { return m_indices.size(); }
0215   /// Get the container of subspace indices
0216   /// @return Reference to the container holding the indices
0217   const Container& indices() const { return m_indices; }
0218 
0219   /// Access subspace index at position i
0220   /// @param i Position index
0221   /// @return The subspace index at position i
0222   IndexType operator[](std::size_t i) const { return m_indices[i]; }
0223 
0224   /// Get iterator to beginning of indices
0225   /// @return Iterator to the first index
0226   auto begin() const { return m_indices.begin(); }
0227   /// Get iterator to end of indices
0228   /// @return Iterator past the last index
0229   auto end() const { return m_indices.end(); }
0230 
0231  private:
0232   Container m_indices;
0233 };
0234 
0235 /// Helper class for fixed subspace operations
0236 ///
0237 /// @tparam FullSize the full size of the subspace
0238 /// @tparam SubspaceSize the size of the subspace
0239 /// @tparam index_t the index type
0240 template <std::size_t FullSize, std::size_t SubspaceSize,
0241           typename index_t = std::uint8_t>
0242 class FixedSubspaceHelper
0243     : public detail::SubspaceHelperBase<
0244           FixedSubspaceHelper<FullSize, SubspaceSize, index_t>, FullSize> {
0245  public:
0246   /// Size of the full parameter space
0247   static constexpr std::size_t kFullSize = FullSize;
0248   /// Size of the subspace parameter space
0249   static constexpr std::size_t kSubspaceSize = SubspaceSize;
0250   static_assert(kSubspaceSize <= kFullSize, "Invalid subspace size");
0251 
0252   /// Type alias for projection matrix (subspace x full)
0253   using Projector = ActsMatrix<kSubspaceSize, kFullSize>;
0254   /// Type alias for expansion matrix (full x subspace)
0255   using Expander = ActsMatrix<kFullSize, kSubspaceSize>;
0256   /// Type alias for subspace vector
0257   using Vector = ActsVector<kSubspaceSize>;
0258   /// Type alias for subspace square matrix
0259   using SquareMatrix = ActsSquareMatrix<kSubspaceSize>;
0260   /// Type alias for left application result matrix
0261   template <std::size_t K>
0262   using ApplyLeftResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0263   /// Type alias for right application result matrix
0264   template <std::size_t N>
0265   using ApplyRightResult = ActsMatrix<kSubspaceSize, kSubspaceSize>;
0266 
0267   /// Type alias for index type used to specify subspace indices
0268   using IndexType = index_t;
0269   /// Type alias for container storing subspace indices
0270   using Container = std::array<IndexType, kSubspaceSize>;
0271 
0272   /// Type alias for size type
0273   using size_type = std::size_t;
0274 
0275   /// Construct a fixed subspace helper with specified indices
0276   ///
0277   /// @tparam other_index_range_t Type of the index range
0278   /// @param indices Range of indices defining the subspace, must match SubspaceSize
0279   template <std::ranges::sized_range other_index_range_t>
0280   explicit FixedSubspaceHelper(const other_index_range_t& indices) {
0281     assert(checkSubspaceIndices(indices, kFullSize, kSubspaceSize) &&
0282            "Invalid indices");
0283     std::transform(indices.begin(), indices.end(), m_indices.begin(),
0284                    [](auto index) { return static_cast<IndexType>(index); });
0285   }
0286 
0287   /// Check if the subspace is empty
0288   /// @return True if the subspace contains no indices (always false for fixed subspaces)
0289   bool empty() const { return m_indices.empty(); }
0290   /// Get the size of the subspace
0291   /// @return Number of indices in the subspace (always SubspaceSize)
0292   std::size_t size() const { return m_indices.size(); }
0293   /// Get the container of subspace indices
0294   /// @return Reference to the array holding the indices
0295   const Container& indices() const { return m_indices; }
0296 
0297   /// Access subspace index at position i
0298   /// @param i Position index
0299   /// @return The subspace index at position i
0300   IndexType operator[](std::uint32_t i) const { return m_indices[i]; }
0301 
0302   /// Get iterator to beginning of indices
0303   /// @return Iterator to the first index
0304   auto begin() const { return m_indices.begin(); }
0305   /// Get iterator to end of indices
0306   /// @return Iterator past the last index
0307   auto end() const { return m_indices.end(); }
0308 
0309   /// Create a projection matrix from full space to subspace
0310   /// @return Matrix that projects from full space to subspace dimensions
0311   Projector projector() const {
0312     Projector result = Projector::Zero();
0313     for (auto [i, index] : enumerate(*this)) {
0314       result(i, index) = 1;
0315     }
0316     return result;
0317   }
0318 
0319   /// Create an expansion matrix from subspace to full space
0320   /// @return Matrix that expands from subspace to full space dimensions
0321   Expander expander() const {
0322     Expander result = Expander::Zero();
0323     for (auto [i, index] : enumerate(*this)) {
0324       result(index, i) = 1;
0325     }
0326     return result;
0327   }
0328 
0329   /// Apply subspace projection on the left side of a matrix
0330   /// @tparam K Number of columns in the input matrix
0331   /// @tparam Derived Eigen matrix type
0332   /// @param matrix Input matrix with rows matching full space size
0333   /// @return Matrix with subspace rows and K columns
0334   template <std::size_t K, typename Derived>
0335   ApplyLeftResult<K> applyLeftOf(
0336       const Eigen::DenseBase<Derived>& matrix) const {
0337     assert(matrix.rows() == kFullSize && "Invalid matrix size");
0338     ApplyLeftResult<K> result = ApplyLeftResult<K>::Zero();
0339     for (auto [i, indexI] : enumerate(*this)) {
0340       for (std::size_t j = 0; j < K; ++j) {
0341         result(i, j) = matrix(indexI, j);
0342       }
0343     }
0344     return result;
0345   }
0346 
0347   /// Apply subspace projection on the right side of a matrix
0348   /// @tparam N Number of rows in the input matrix
0349   /// @tparam Derived Eigen matrix type
0350   /// @param matrix Input matrix with columns matching subspace size
0351   /// @return Matrix with N rows and subspace columns
0352   template <std::size_t N, typename Derived>
0353   ApplyRightResult<N> applyRightOf(
0354       const Eigen::DenseBase<Derived>& matrix) const {
0355     assert(matrix.cols() == kSubspaceSize && "Invalid matrix size");
0356     ApplyRightResult<N> result = ApplyRightResult<N>::Zero();
0357     for (std::size_t i = 0; i < N; ++i) {
0358       for (auto [j, indexJ] : enumerate(*this)) {
0359         result(i, j) = matrix(i, indexJ);
0360       }
0361     }
0362     return result;
0363   }
0364 
0365   /// Project a full-space vector to subspace
0366   /// @tparam Derived Eigen vector type
0367   /// @param fullVector Input vector with full space dimensions
0368   /// @return Vector projected to subspace dimensions
0369   template <typename Derived>
0370   Vector projectVector(const Eigen::DenseBase<Derived>& fullVector) const {
0371     assert(fullVector.size() == kFullSize && "Invalid full vector size");
0372     Vector result = Vector::Zero();
0373     for (auto [i, index] : enumerate(*this)) {
0374       result(i) = fullVector(index);
0375     }
0376     return result;
0377   }
0378 
0379   /// Project a full-space square matrix to subspace
0380   /// @tparam Derived Eigen matrix type
0381   /// @param fullMatrix Input square matrix with full space dimensions
0382   /// @return Square matrix projected to subspace dimensions
0383   template <typename Derived>
0384   SquareMatrix projectMatrix(
0385       const Eigen::DenseBase<Derived>& fullMatrix) const {
0386     assert(fullMatrix.rows() == kFullSize && fullMatrix.cols() == kFullSize &&
0387            "Invalid full matrix size");
0388     SquareMatrix result = SquareMatrix::Zero();
0389     for (auto [i, indexI] : enumerate(*this)) {
0390       for (auto [j, indexJ] : enumerate(*this)) {
0391         result(i, j) = fullMatrix(indexI, indexJ);
0392       }
0393     }
0394     return result;
0395   }
0396 
0397  private:
0398   Container m_indices{};
0399 };
0400 
0401 /// @brief Helper type for fixed-size bound parameter subspaces
0402 /// @details Provides utilities for working with fixed-size subspaces of bound track parameters
0403 template <std::size_t SubspaceSize>
0404 using FixedBoundSubspaceHelper =
0405     FixedSubspaceHelper<Acts::eBoundSize, SubspaceSize, std::uint8_t>;
0406 
0407 /// @typedef VariableBoundSubspaceHelper
0408 /// Helper type for variable-size bound parameter subspaces.
0409 /// Provides utilities for working with variable-size subspaces of bound track
0410 /// parameters.
0411 using VariableBoundSubspaceHelper =
0412     VariableSubspaceHelper<Acts::eBoundSize, std::uint8_t>;
0413 
0414 /// Convert a projector to subspace indices
0415 ///
0416 /// @tparam kFullSize the full size of the subspace
0417 /// @tparam Derived the derived type
0418 ///
0419 /// @param projector the projector
0420 ///
0421 /// @return the subspace indices
0422 template <std::size_t kFullSize, typename Derived>
0423 SubspaceIndices<kFullSize> projectorToSubspaceIndices(
0424     const Eigen::DenseBase<Derived>& projector) {
0425   auto rows = static_cast<std::size_t>(projector.rows());
0426   auto cols = static_cast<std::size_t>(projector.cols());
0427   assert(cols == kFullSize && rows <= kFullSize && "Invalid projector size");
0428   SubspaceIndices<kFullSize> result;
0429   result.fill(kFullSize);
0430   for (std::size_t i = 0; i < rows; ++i) {
0431     for (std::size_t j = 0; j < cols; ++j) {
0432       assert((projector(i, j) == 0 || projector(i, j) == 1) &&
0433              "Invalid projector value");
0434       if (projector(i, j) == 1) {
0435         result[i] = j;
0436       }
0437     }
0438   }
0439   return result;
0440 }
0441 
0442 }  // namespace Acts