ActsExamples/EventData/Measurement.hpp

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/SubspaceHelpers.hpp"
0014 #include "Acts/EventData/Types.hpp"
0015 #include "Acts/Geometry/GeometryIdentifier.hpp"
0016 #include "Acts/Utilities/ContainerIterator.hpp"
0017 #include "ActsExamples/EventData/GeometryContainers.hpp"
0018 #include "ActsExamples/EventData/Index.hpp"
0019 #include "ActsExamples/EventData/IndexSourceLink.hpp"
0020 #include "ActsExamples/EventData/MeasurementConcept.hpp"
0021 #include "ActsExamples/EventData/SimParticle.hpp"
0022
0023 #include <cstddef>
0024 #include <iterator>
0025 #include <type_traits>
0026 #include <vector>
0027
0028 #include <boost/container/static_vector.hpp>
0029
0030 namespace ActsExamples {
0031
0032 template <typename Derived, std::size_t FullSize, bool ReadOnly>
0033 class MeasurementProxyBase;
0034 template <std::size_t FullSize, std::size_t Size, bool ReadOnly>
0035 class FixedMeasurementProxy;
0036 template <std::size_t FullSize, bool ReadOnly>
0037 class VariableMeasurementProxy;
0038
0039 template <std::size_t Size>
0040 using FixedBoundMeasurementProxy =
0041     FixedMeasurementProxy<Acts::eBoundSize, Size, false>;
0042 template <std::size_t Size>
0043 using ConstFixedBoundMeasurementProxy =
0044     FixedMeasurementProxy<Acts::eBoundSize, Size, true>;
0045 using VariableBoundMeasurementProxy =
0046     VariableMeasurementProxy<Acts::eBoundSize, false>;
0047 using ConstVariableBoundMeasurementProxy =
0048     VariableMeasurementProxy<Acts::eBoundSize, true>;
0049
0050 /// @brief A container to store and access measurements
0051 ///
0052 /// This container stores measurements of different sizes and provides
0053 /// access to them through fixed-size and variable-size proxies.
0054 ///
0055 /// The measurements are stored densely in a flat buffer and the proxies
0056 /// provide access to the individual measurements.
0057 class MeasurementContainer {
0058  public:
0059   using size_type = std::size_t;
0060   using Index = size_type;
0061   template <std::size_t Size>
0062   using FixedProxy = FixedMeasurementProxy<Acts::eBoundSize, Size, false>;
0063   template <std::size_t Size>
0064   using ConstFixedProxy = FixedMeasurementProxy<Acts::eBoundSize, Size, true>;
0065   using VariableProxy = VariableMeasurementProxy<Acts::eBoundSize, false>;
0066   using ConstVariableProxy = VariableMeasurementProxy<Acts::eBoundSize, true>;
0067   using OrderedIndices = GeometryIdMultiset<IndexSourceLink>;
0068
0069   MeasurementContainer();
0070
0071   /// @brief Get the number of measurements
0072   /// @return The number of measurements
0073   std::size_t size() const;
0074
0075   /// @brief Reserve space for a number of measurements
0076   /// @param size The number of measurements to reserve space for
0077   void reserve(std::size_t size);
0078
0079   /// @brief Add a measurement of a given size
0080   /// @param size The size of the measurement
0081   /// @param geometryId The geometry identifier of the measurement surface
0082   /// @return The index of the added measurement
0083   Index addMeasurement(std::uint8_t size, Acts::GeometryIdentifier geometryId);
0084
0085   /// @brief Get a variable-size measurement proxy
0086   /// @param index The index of the measurement
0087   /// @return The variable-size measurement proxy
0088   VariableProxy at(Index index);
0089   /// @brief Get a const variable-size measurement proxy
0090   /// @param index The index of the measurement
0091   /// @return The const variable-size measurement proxy
0092   ConstVariableProxy at(Index index) const;
0093
0094   /// @brief Get a variable-size measurement proxy
0095   /// @param index The index of the measurement
0096   /// @return The variable-size measurement proxy
0097   VariableProxy getMeasurement(Index index);
0098   /// @brief Get a const variable-size measurement proxy
0099   /// @param index The index of the measurement
0100   /// @return The const variable-size measurement proxy
0101   ConstVariableProxy getMeasurement(Index index) const;
0102
0103   /// @brief Get a fixed-size measurement proxy
0104   /// @tparam Size The size of the measurement
0105   /// @param index The index of the measurement
0106   /// @return The fixed-size measurement proxy
0107   template <std::size_t Size>
0108   FixedProxy<Size> getMeasurement(Index index) {
0109     return FixedProxy<Size>{*this, index};
0110   }
0111   /// @brief Get a const fixed-size measurement proxy
0112   /// @tparam Size The size of the measurement
0113   /// @param index The index of the measurement
0114   /// @return The const fixed-size measurement proxy
0115   template <std::size_t Size>
0116   ConstFixedProxy<Size> getMeasurement(Index index) const {
0117     return ConstFixedProxy<Size>{*this, index};
0118   }
0119
0120   /// @brief Make a measurement of a given size
0121   /// @param size The size of the measurement
0122   /// @param geometryId The geometry identifier of the measurement surface
0123   /// @return The variable-size measurement proxy
0124   VariableProxy makeMeasurement(std::uint8_t size,
0125                                 Acts::GeometryIdentifier geometryId);
0126   /// @brief Make a fixed-size measurement
0127   /// @tparam Size The size of the measurement
0128   /// @param geometryId The geometry identifier of the measurement surface
0129   /// @return The fixed-size measurement proxy
0130   template <std::size_t Size>
0131   FixedProxy<Size> makeMeasurement(Acts::GeometryIdentifier geometryId) {
0132     return getMeasurement<Size>(addMeasurement(Size, geometryId));
0133   }
0134
0135   template <MeasurementConcept OtherDerived>
0136   VariableProxy copyMeasurement(const OtherDerived& other);
0137   template <MeasurementConcept OtherDerived, std::size_t Size>
0138   FixedProxy<Size> copyMeasurement(const OtherDerived& other);
0139
0140   template <typename... Args>
0141   VariableProxy emplaceMeasurement(std::uint8_t size,
0142                                    Acts::GeometryIdentifier geometryId,
0143                                    Args&&... args);
0144
0145   template <std::size_t Size, typename... Args>
0146   FixedProxy<Size> emplaceMeasurement(Acts::GeometryIdentifier geometryId,
0147                                       Args&&... args);
0148
0149   const OrderedIndices& orderedIndices() const;
0150
0151   using iterator = Acts::ContainerIterator<MeasurementContainer, VariableProxy,
0152                                            Index, false>;
0153   using const_iterator =
0154       Acts::ContainerIterator<const MeasurementContainer, ConstVariableProxy,
0155                               Index, true>;
0156
0157   iterator begin();
0158   iterator end();
0159   const_iterator begin() const;
0160   const_iterator end() const;
0161   const_iterator cbegin() const;
0162   const_iterator cend() const;
0163
0164  public:
0165   struct MeasurementEntry {
0166     std::size_t subspaceIndexOffset{};
0167     std::size_t parameterOffset{};
0168     std::size_t covarianceOffset{};
0169     std::uint8_t size{};
0170   };
0171
0172   std::vector<MeasurementEntry> m_entries;
0173
0174   std::vector<Acts::GeometryIdentifier> m_geometryIds;
0175   std::vector<std::uint8_t> m_subspaceIndices;
0176   std::vector<double> m_parameters;
0177   std::vector<double> m_covariances;
0178
0179   OrderedIndices m_orderedIndices;
0180 };
0181
0182 /// @brief Base class for measurement proxies
0183 ///
0184 /// This class provides common functionality for fixed-size and variable-size
0185 /// measurement proxies.
0186 ///
0187 /// @tparam Derived The derived measurement proxy class
0188 /// @tparam FullSize The full size of the measurement
0189 /// @tparam ReadOnly Whether the proxy is read-only
0190 template <typename Derived, std::size_t FullSize, bool ReadOnly>
0191 class MeasurementProxyBase {
0192  public:
0193   using Index = MeasurementContainer::Index;
0194   using SubspaceIndex = std::uint8_t;
0195   using Scalar = double;
0196
0197   using FullVector = Acts::ActsVector<FullSize>;
0198   using FullSquareMatrix = Acts::ActsSquareMatrix<FullSize>;
0199
0200   using Container = std::conditional_t<ReadOnly, const MeasurementContainer,
0201                                        MeasurementContainer>;
0202
0203   MeasurementProxyBase(Container& container_, Index index_)
0204       : m_container(&container_), m_index(index_) {}
0205   template <typename OtherDerived, bool OtherReadOnly>
0206   explicit MeasurementProxyBase(
0207       const MeasurementProxyBase<OtherDerived, FullSize, OtherReadOnly>& other)
0208     requires(ReadOnly == OtherReadOnly || ReadOnly)
0209       : m_container(&other.container()), m_index(other.index()) {}
0210
0211   /// @brief Get the container of the measurement
0212   /// @return The container of the measurement
0213   Container& container() const { return *m_container; }
0214   /// @brief Get the index of the measurement
0215   /// @return The index of the measurement
0216   Index index() const { return m_index; }
0217
0218   /// @brief Get the size of the measurement
0219   /// @return The size of the measurement
0220   std::size_t size() const { return container().m_entries.at(m_index).size; }
0221
0222   /// @brief Check if the measurement contains a subspace index
0223   /// @param i The subspace index
0224   /// @return True if the measurement contains the subspace index
0225   template <typename indices_t>
0226   bool contains(indices_t i) const {
0227     return self().subspaceHelper().contains(i);
0228   }
0229
0230   /// @brief Get the index of a subspace index in the measurement
0231   /// @param i The subspace index
0232   /// @return The index of the subspace index in the measurement
0233   template <typename indices_t>
0234   std::size_t indexOf(indices_t i) const {
0235     return self().subspaceHelper().indexOf(i);
0236   }
0237
0238   /// @brief Get the geometry ID of the measurement
0239   /// @return The geometry ID
0240   Acts::GeometryIdentifier geometryId() const {
0241     return container().m_geometryIds.at(m_index);
0242   }
0243
0244   /// @brief Set the subspace indices of the measurement
0245   /// @param indices The subspace indices
0246   template <typename IndexContainer>
0247   void setSubspaceIndices(const IndexContainer& indices)
0248     requires(!ReadOnly)
0249   {
0250     assert(Acts::checkSubspaceIndices(indices, FullSize, size()) &&
0251            "Invalid indices");
0252     std::transform(indices.begin(), indices.end(),
0253                    self().subspaceIndexVector().begin(),
0254                    [](auto index) { return static_cast<Index>(index); });
0255   }
0256
0257   /// @brief Get the measurement as a full-size vector
0258   /// @return The full-size measurement vector
0259   FullVector fullParameters() const {
0260     return self().subspaceHelper().expandVector(self().parameters());
0261   }
0262
0263   /// @brief Get the covariance as a full-size square matrix
0264   /// @return The full-size covariance matrix
0265   FullSquareMatrix fullCovariance() const {
0266     return self().subspaceHelper().expandMatrix(self().covariance());
0267   }
0268
0269   /// @brief Construct the measurement from a subspace vector,
0270   /// parameters, and covariance.
0271   ///
0272   template <typename Subspace, typename ParameterDerived,
0273             typename CovarianceDerived>
0274   void fill(Subspace&& subspace,
0275             const Eigen::DenseBase<ParameterDerived>& parameters,
0276             const Eigen::DenseBase<CovarianceDerived>& covariance)
0277     requires(!ReadOnly)
0278   {
0279     self().setSubspaceIndices(std::forward<Subspace>(subspace));
0280     self().parameters() = parameters;
0281     self().covariance() = covariance;
0282   }
0283
0284   /// @brief Construct the measurement from a subspace vector,
0285   /// parameters, and covariance.
0286   ///
0287   template <MeasurementConcept OtherDerived>
0288   void fill(const OtherDerived& other)
0289     requires(!ReadOnly)
0290   {
0291     assert(size() == other.size() && "Size mismatch");
0292     fill(other.subspaceIndexVector(), other.parameters(), other.covariance());
0293   }
0294
0295   /// @brief Copy the data from another measurement
0296   /// @tparam OtherDerived The derived measurement proxy class of the other
0297   ///         measurement
0298   /// @param other The other measurement proxy
0299   template <typename OtherDerived>
0300   void copyFrom(const OtherDerived& other)
0301     requires(!ReadOnly) && requires {
0302       { this->fill(other) };
0303     }
0304   {
0305     fill(other);
0306   }
0307
0308  protected:
0309   Derived& self() { return static_cast<Derived&>(*this); }
0310   const Derived& self() const { return static_cast<const Derived&>(*this); }
0311
0312   Container* m_container;
0313   Index m_index;
0314 };
0315
0316 /// @brief Fixed-size measurement proxy
0317 ///
0318 /// This class provides access to a fixed-size measurement in a measurement
0319 /// container.
0320 ///
0321 /// @tparam FullSize The full size of the measurement
0322 /// @tparam Size The size of the measurement
0323 /// @tparam ReadOnly Whether the proxy is read-only
0324 template <std::size_t FullSize, std::size_t Size, bool ReadOnly>
0325 class FixedMeasurementProxy
0326     : public MeasurementProxyBase<
0327           FixedMeasurementProxy<FullSize, Size, ReadOnly>, FullSize, ReadOnly> {
0328  public:
0329   using Base =
0330       MeasurementProxyBase<FixedMeasurementProxy<FullSize, Size, ReadOnly>,
0331                            FullSize, ReadOnly>;
0332   using Index = typename Base::Index;
0333   using SubspaceIndex = typename Base::SubspaceIndex;
0334   using Scalar = typename Base::Scalar;
0335   using Container = typename Base::Container;
0336
0337   using SubspaceHelper = Acts::FixedSubspaceHelper<FullSize, Size>;
0338
0339   using SubspaceVector = Eigen::Matrix<SubspaceIndex, Size, 1>;
0340   using SubspaceVectorMap =
0341       std::conditional_t<ReadOnly, Eigen::Map<const SubspaceVector>,
0342                          Eigen::Map<SubspaceVector>>;
0343
0344   using ParametersVector = Eigen::Matrix<Scalar, Size, 1>;
0345   using ParametersVectorMap =
0346       std::conditional_t<ReadOnly, Eigen::Map<const ParametersVector>,
0347                          Eigen::Map<ParametersVector>>;
0348
0349   using CovarianceMatrix = Eigen::Matrix<Scalar, Size, Size>;
0350   using CovarianceMatrixMap =
0351       std::conditional_t<ReadOnly, Eigen::Map<const CovarianceMatrix>,
0352                          Eigen::Map<CovarianceMatrix>>;
0353
0354   FixedMeasurementProxy(Container& container_, Index index_)
0355       : Base(container_, index_) {
0356     assert(container().m_entries.at(index()).size == Size && "Size mismatch");
0357   }
0358   template <typename OtherDerived, bool OtherReadOnly>
0359   explicit FixedMeasurementProxy(
0360       const MeasurementProxyBase<OtherDerived, FullSize, OtherReadOnly>& other)
0361     requires(ReadOnly == OtherReadOnly || ReadOnly)
0362       : Base(other) {
0363     assert(container().m_entries.at(index()).size == Size && "Size mismatch");
0364   }
0365
0366   using Base::container;
0367   using Base::index;
0368
0369   /// @brief Get the size of the measurement
0370   /// @return The size of the measurement
0371   static constexpr std::size_t size() { return Size; }
0372
0373   /// @brief Get the subspace helper for the measurement
0374   /// @return The subspace helper
0375   SubspaceHelper subspaceHelper() const {
0376     return SubspaceHelper{subspaceIndexVector()};
0377   }
0378
0379   /// @brief Get the subspace indices of the measurement
0380   /// @return The subspace indices
0381   Acts::SubspaceIndices<Size> subspaceIndices() const {
0382     return subspaceHelper().indices();
0383   }
0384
0385   /// @brief Get the subspace index vector of the measurement
0386   /// @return The subspace index vector
0387   SubspaceVectorMap subspaceIndexVector() const {
0388     return SubspaceVectorMap{
0389         container().m_subspaceIndices.data() +
0390         container().m_entries.at(index()).subspaceIndexOffset};
0391   }
0392
0393   /// @brief Get the parameters of the measurement
0394   /// @return The parameters
0395   ParametersVectorMap parameters() const {
0396     return ParametersVectorMap{
0397         container().m_parameters.data() +
0398         container().m_entries.at(index()).parameterOffset};
0399   }
0400
0401   /// @brief Get the covariance of the measurement
0402   /// @return The covariance
0403   CovarianceMatrixMap covariance() const {
0404     return CovarianceMatrixMap{
0405         container().m_covariances.data() +
0406         container().m_entries.at(index()).covarianceOffset};
0407   }
0408 };
0409
0410 /// @brief Variable-size measurement proxy
0411 ///
0412 /// This class provides access to a variable-size measurement in a measurement
0413 /// container.
0414 ///
0415 /// @tparam FullSize The full size of the measurement
0416 /// @tparam ReadOnly Whether the proxy is read-only
0417 template <std::size_t FullSize, bool ReadOnly>
0418 class VariableMeasurementProxy
0419     : public MeasurementProxyBase<VariableMeasurementProxy<FullSize, ReadOnly>,
0420                                   FullSize, ReadOnly> {
0421  public:
0422   using Base =
0423       MeasurementProxyBase<VariableMeasurementProxy<FullSize, ReadOnly>,
0424                            FullSize, ReadOnly>;
0425   using Index = typename Base::Index;
0426   using SubspaceIndex = typename Base::SubspaceIndex;
0427   using Scalar = typename Base::Scalar;
0428   using Container = typename Base::Container;
0429
0430   using SubspaceHelper = Acts::VariableSubspaceHelper<FullSize>;
0431
0432   using SubspaceVector = Eigen::Matrix<SubspaceIndex, Eigen::Dynamic, 1>;
0433   using SubspaceVectorMap =
0434       std::conditional_t<ReadOnly, Eigen::Map<const SubspaceVector>,
0435                          Eigen::Map<SubspaceVector>>;
0436
0437   using ParametersVector = Eigen::Matrix<Scalar, Eigen::Dynamic, 1>;
0438   using ParametersVectorMap =
0439       std::conditional_t<ReadOnly, Eigen::Map<const ParametersVector>,
0440                          Eigen::Map<ParametersVector>>;
0441
0442   using CovarianceMatrix =
0443       Eigen::Matrix<Scalar, Eigen::Dynamic, Eigen::Dynamic>;
0444   using CovarianceMatrixMap =
0445       std::conditional_t<ReadOnly, Eigen::Map<const CovarianceMatrix>,
0446                          Eigen::Map<CovarianceMatrix>>;
0447
0448   VariableMeasurementProxy(Container& container_, Index index_)
0449       : Base(container_, index_) {}
0450   template <typename OtherDerived, bool OtherReadOnly>
0451   explicit VariableMeasurementProxy(
0452       const MeasurementProxyBase<OtherDerived, FullSize, OtherReadOnly>& other)
0453     requires(ReadOnly == OtherReadOnly || ReadOnly)
0454       : Base(other) {}
0455
0456   using Base::container;
0457   using Base::index;
0458
0459   /// @brief Get the subspace helper for the measurement
0460   /// @return The subspace helper
0461   SubspaceHelper subspaceHelper() const {
0462     return SubspaceHelper{subspaceIndexVector()};
0463   }
0464
0465   /// @brief Get the subspace indices of the measurement
0466   /// @return The subspace indices
0467   SubspaceVectorMap subspaceIndexVector() const {
0468     const auto size = static_cast<Eigen::Index>(this->size());
0469     return SubspaceVectorMap{
0470         container().m_subspaceIndices.data() +
0471             container().m_entries.at(index()).subspaceIndexOffset,
0472         size};
0473   }
0474
0475   /// @brief Get the parameters of the measurement
0476   /// @return The parameters
0477   ParametersVectorMap parameters() const {
0478     const auto size = static_cast<Eigen::Index>(this->size());
0479     return ParametersVectorMap{
0480         container().m_parameters.data() +
0481             container().m_entries.at(index()).parameterOffset,
0482         size};
0483   }
0484
0485   /// @brief Get the covariance of the measurement
0486   /// @return The covariance
0487   CovarianceMatrixMap covariance() const {
0488     const auto size = static_cast<Eigen::Index>(this->size());
0489     return CovarianceMatrixMap{
0490         container().m_covariances.data() +
0491             container().m_entries.at(index()).covarianceOffset,
0492         size, size};
0493   }
0494 };
0495
0496 template <MeasurementConcept OtherDerived>
0497 MeasurementContainer::VariableProxy MeasurementContainer::copyMeasurement(
0498     const OtherDerived& other) {
0499   VariableProxy meas = makeMeasurement(other.size(), other.geometryId());
0500   meas.copyFrom(other);
0501   return meas;
0502 }
0503
0504 template <MeasurementConcept OtherDerived, std::size_t Size>
0505 MeasurementContainer::FixedProxy<Size> MeasurementContainer::copyMeasurement(
0506     const OtherDerived& other) {
0507   FixedProxy<Size> meas = makeMeasurement<Size>(other.geometryId());
0508   meas.copyFrom(other);
0509   return meas;
0510 }
0511
0512 template <typename... Args>
0513 MeasurementContainer::VariableProxy MeasurementContainer::emplaceMeasurement(
0514     std::uint8_t size, Acts::GeometryIdentifier geometryId, Args&&... args) {
0515   VariableProxy meas = makeMeasurement(size, geometryId);
0516
0517   meas.fill(std::forward<Args>(args)...);
0518
0519   return meas;
0520 }
0521
0522 template <std::size_t Size, typename... Args>
0523 MeasurementContainer::FixedProxy<Size> MeasurementContainer::emplaceMeasurement(
0524     Acts::GeometryIdentifier geometryId, Args&&... args) {
0525   FixedProxy<Size> meas = makeMeasurement<Size>(geometryId);
0526
0527   meas.fill(std::forward<Args>(args)...);
0528
0529   return meas;
0530 }
0531
0532 static_assert(
0533     std::random_access_iterator<MeasurementContainer::iterator> &&
0534     std::random_access_iterator<MeasurementContainer::const_iterator>);
0535
0536 using MeasurementSimHitsMap = IndexMultimap<Index>;
0537 using MeasurementParticlesMap = IndexMultimap<SimBarcode>;
0538
0539 using SimHitMeasurementsMap = InverseMultimap<Index>;
0540 using ParticleMeasurementsMap = InverseMultimap<SimBarcode>;
0541
0542 }  // namespace ActsExamples