Acts/EventData/VectorMultiTrajectory.hpp

0001 // This file is part of the Acts project.
0002 //
0003 // Copyright (C) 2022 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 http://mozilla.org/MPL/2.0/.
0008
0009 #pragma once
0010
0011 #include "Acts/Definitions/TrackParametrization.hpp"
0012 #include "Acts/EventData/MultiTrajectory.hpp"
0013 #include "Acts/EventData/MultiTrajectoryBackendConcept.hpp"
0014 #include "Acts/EventData/SourceLink.hpp"
0015 #include "Acts/EventData/TrackStatePropMask.hpp"
0016 #include "Acts/EventData/detail/DynamicColumn.hpp"
0017 #include "Acts/EventData/detail/DynamicKeyIterator.hpp"
0018 #include "Acts/Utilities/Concepts.hpp"
0019 #include "Acts/Utilities/HashedString.hpp"
0020 #include "Acts/Utilities/Helpers.hpp"
0021 #include "Acts/Utilities/ThrowAssert.hpp"
0022
0023 #include <any>
0024 #include <cassert>
0025 #include <cstddef>
0026 #include <iosfwd>
0027 #include <memory>
0028 #include <optional>
0029 #include <stdexcept>
0030 #include <string>
0031 #include <type_traits>
0032 #include <unordered_map>
0033 #include <utility>
0034 #include <vector>
0035
0036 #include <boost/histogram.hpp>
0037
0038 namespace Acts {
0039 class Surface;
0040 template <typename T>
0041 struct IsReadOnlyMultiTrajectory;
0042
0043 namespace detail_vmt {
0044
0045 using MultiTrajectoryTraits::IndexType;
0046 constexpr auto kInvalid = MultiTrajectoryTraits::kInvalid;
0047 constexpr auto MeasurementSizeMax = MultiTrajectoryTraits::MeasurementSizeMax;
0048
0049 class VectorMultiTrajectoryBase {
0050  public:
0051   struct Statistics {
0052     using axis_t = boost::histogram::axis::variant<
0053         boost::histogram::axis::category<std::string>,
0054         boost::histogram::axis::category<>>;
0055
0056     using axes_t = std::vector<axis_t>;
0057     using hist_t = boost::histogram::histogram<axes_t>;
0058
0059     hist_t hist;
0060
0061     void toStream(std::ostream& os, std::size_t n = 1);
0062   };
0063
0064   template <typename T>
0065   Statistics statistics(T& instance) const {
0066     using namespace boost::histogram;
0067     using cat = axis::category<std::string>;
0068
0069     Statistics::axes_t axes;
0070     axes.emplace_back(cat({
0071         "count",
0072         "index",
0073         "parPred",
0074         "covPred",
0075         "parFilt",
0076         "covFilt",
0077         "parSmth",
0078         "covSmth",
0079         "meas",
0080         "measOffset",
0081         "measCov",
0082         "measCovOffset",
0083         "jac",
0084         "sourceLinks",
0085         "projectors",
0086     }));
0087
0088     axes.emplace_back(axis::category<>({0, 1}));
0089
0090     auto h = make_histogram(axes);
0091
0092     for (IndexType i = 0; i < instance.size(); i++) {
0093       auto ts = instance.getTrackState(i);
0094
0095       bool isMeas = ts.typeFlags().test(TrackStateFlag::MeasurementFlag);
0096
0097       h("count", isMeas);
0098
0099       h("index", isMeas, weight(sizeof(IndexData)));
0100
0101       using scalar = typename decltype(ts.predicted())::Scalar;
0102       std::size_t par_size = eBoundSize * sizeof(scalar);
0103       std::size_t cov_size = eBoundSize * eBoundSize * sizeof(scalar);
0104
0105       const IndexData& index = m_index[i];
0106       if (ts.hasPredicted() &&
0107           ACTS_CHECK_BIT(index.allocMask, TrackStatePropMask::Predicted)) {
0108         h("parPred", isMeas, weight(par_size));
0109         h("covPred", isMeas, weight(cov_size));
0110       }
0111       if (ts.hasFiltered() &&
0112           ACTS_CHECK_BIT(index.allocMask, TrackStatePropMask::Filtered)) {
0113         h("parFilt", isMeas, weight(par_size));
0114         h("covFilt", isMeas, weight(cov_size));
0115       }
0116       if (ts.hasSmoothed() &&
0117           ACTS_CHECK_BIT(index.allocMask, TrackStatePropMask::Smoothed)) {
0118         h("parSmth", isMeas, weight(par_size));
0119         h("covSmth", isMeas, weight(cov_size));
0120       }
0121       h("sourceLinks", isMeas, weight(sizeof(SourceLink)));
0122       h("measOffset", isMeas,
0123         weight(sizeof(decltype(m_measOffset)::value_type)));
0124       h("measCovOffset", isMeas,
0125         weight(sizeof(decltype(m_measCovOffset)::value_type)));
0126       if (ts.hasCalibrated() &&
0127           ACTS_CHECK_BIT(index.allocMask, TrackStatePropMask::Calibrated)) {
0128         std::size_t meas_size = ts.calibratedSize() * sizeof(scalar);
0129         std::size_t meas_cov_size =
0130             ts.calibratedSize() * ts.calibratedSize() * sizeof(scalar);
0131
0132         h("meas", isMeas, weight(meas_size));
0133         h("measCov", isMeas, weight(meas_cov_size));
0134         h("sourceLinks", isMeas, weight(sizeof(const SourceLink)));
0135         h("projectors", isMeas, weight(sizeof(ProjectorBitset)));
0136       }
0137
0138       if (ts.hasJacobian() &&
0139           ACTS_CHECK_BIT(index.allocMask, TrackStatePropMask::Jacobian)) {
0140         h("jac", isMeas, weight(cov_size));
0141       }
0142     }
0143
0144     return Statistics{h};
0145   }
0146
0147  protected:
0148   struct IndexData {
0149     IndexType ipredicted = kInvalid;
0150     IndexType ifiltered = kInvalid;
0151     IndexType ismoothed = kInvalid;
0152     IndexType ijacobian = kInvalid;
0153     IndexType iprojector = kInvalid;
0154
0155     float chi2 = 0;
0156     double pathLength = 0;
0157     TrackStateType::raw_type typeFlags{};
0158
0159     IndexType iuncalibrated = kInvalid;
0160     IndexType icalibratedsourcelink = kInvalid;
0161     IndexType measdim = kInvalid;
0162
0163     TrackStatePropMask allocMask = TrackStatePropMask::None;
0164   };
0165
0166   VectorMultiTrajectoryBase() = default;
0167
0168   VectorMultiTrajectoryBase(const VectorMultiTrajectoryBase& other)
0169       : m_index{other.m_index},
0170         m_previous{other.m_previous},
0171         m_next{other.m_next},
0172         m_params{other.m_params},
0173         m_cov{other.m_cov},
0174         m_meas{other.m_meas},
0175         m_measOffset{other.m_measOffset},
0176         m_measCov{other.m_measCov},
0177         m_measCovOffset{other.m_measCovOffset},
0178         m_jac{other.m_jac},
0179         m_sourceLinks{other.m_sourceLinks},
0180         m_projectors{other.m_projectors},
0181         m_referenceSurfaces{other.m_referenceSurfaces} {
0182     for (const auto& [key, value] : other.m_dynamic) {
0183       m_dynamic.insert({key, value->clone()});
0184     }
0185     m_dynamicKeys = other.m_dynamicKeys;
0186   };
0187
0188   VectorMultiTrajectoryBase(VectorMultiTrajectoryBase&& other) = default;
0189
0190   // BEGIN INTERFACE HELPER
0191   template <typename T>
0192   static constexpr bool has_impl(T& instance, HashedString key,
0193                                  IndexType istate) {
0194     using namespace Acts::HashedStringLiteral;
0195     switch (key) {
0196       case "predicted"_hash:
0197         return instance.m_index[istate].ipredicted != kInvalid;
0198       case "filtered"_hash:
0199         return instance.m_index[istate].ifiltered != kInvalid;
0200       case "smoothed"_hash:
0201         return instance.m_index[istate].ismoothed != kInvalid;
0202       case "calibrated"_hash:
0203         return instance.m_measOffset[istate] != kInvalid;
0204       case "calibratedCov"_hash:
0205         return instance.m_measCovOffset[istate] != kInvalid;
0206       case "jacobian"_hash:
0207         return instance.m_index[istate].ijacobian != kInvalid;
0208       case "projector"_hash:
0209         return instance.m_index[istate].iprojector != kInvalid;
0210       case "uncalibratedSourceLink"_hash:
0211         return instance.m_sourceLinks[instance.m_index[istate].iuncalibrated]
0212             .has_value();
0213       case "previous"_hash:
0214       case "next"_hash:
0215       case "referenceSurface"_hash:
0216       case "measdim"_hash:
0217       case "chi2"_hash:
0218       case "pathLength"_hash:
0219       case "typeFlags"_hash:
0220         return true;
0221       default:
0222         return instance.m_dynamic.find(key) != instance.m_dynamic.end();
0223     }
0224   }
0225
0226   template <bool EnsureConst, typename T>
0227   static std::any component_impl(T& instance, HashedString key,
0228                                  IndexType istate) {
0229     if constexpr (EnsureConst) {
0230       static_assert(std::is_const_v<std::remove_reference_t<T>>,
0231                     "Is not const");
0232     }
0233     using namespace Acts::HashedStringLiteral;
0234     switch (key) {
0235       case "previous"_hash:
0236         return &instance.m_previous[istate];
0237       case "next"_hash:
0238         return &instance.m_next[istate];
0239       case "predicted"_hash:
0240         return &instance.m_index[istate].ipredicted;
0241       case "filtered"_hash:
0242         return &instance.m_index[istate].ifiltered;
0243       case "smoothed"_hash:
0244         return &instance.m_index[istate].ismoothed;
0245       case "projector"_hash:
0246         return &instance.m_projectors[instance.m_index[istate].iprojector];
0247       case "measdim"_hash:
0248         return &instance.m_index[istate].measdim;
0249       case "chi2"_hash:
0250         return &instance.m_index[istate].chi2;
0251       case "pathLength"_hash:
0252         return &instance.m_index[istate].pathLength;
0253       case "typeFlags"_hash:
0254         return &instance.m_index[istate].typeFlags;
0255       default:
0256         auto it = instance.m_dynamic.find(key);
0257         if (it == instance.m_dynamic.end()) {
0258           throw std::runtime_error("Unable to handle this component");
0259         }
0260         std::conditional_t<EnsureConst, const detail::DynamicColumnBase*,
0261                            detail::DynamicColumnBase*>
0262             col = it->second.get();
0263         assert(col && "Dynamic column is null");
0264         return col->get(istate);
0265     }
0266   }
0267
0268   template <typename T>
0269   static bool hasColumn_impl(T& instance, HashedString key) {
0270     using namespace Acts::HashedStringLiteral;
0271     switch (key) {
0272       case "predicted"_hash:
0273       case "filtered"_hash:
0274       case "smoothed"_hash:
0275       case "calibrated"_hash:
0276       case "calibratedCov"_hash:
0277       case "jacobian"_hash:
0278       case "projector"_hash:
0279       case "previous"_hash:
0280       case "next"_hash:
0281       case "uncalibratedSourceLink"_hash:
0282       case "referenceSurface"_hash:
0283       case "measdim"_hash:
0284       case "chi2"_hash:
0285       case "pathLength"_hash:
0286       case "typeFlags"_hash:
0287         return true;
0288       default:
0289         return instance.m_dynamic.find(key) != instance.m_dynamic.end();
0290     }
0291   }
0292
0293  public:
0294   detail::DynamicKeyRange<detail::DynamicColumnBase> dynamicKeys_impl() const {
0295     return {m_dynamic.begin(), m_dynamic.end()};
0296   }
0297
0298   // END INTERFACE HELPER
0299
0300  public:
0301   IndexType calibratedSize_impl(IndexType istate) const {
0302     return m_index[istate].measdim;
0303   }
0304
0305   SourceLink getUncalibratedSourceLink_impl(IndexType istate) const {
0306     return m_sourceLinks[m_index[istate].iuncalibrated].value();
0307   }
0308
0309   const Surface* referenceSurface_impl(IndexType istate) const {
0310     return m_referenceSurfaces[istate].get();
0311   }
0312
0313  protected:
0314   /// index to map track states to the corresponding
0315   std::vector<IndexData> m_index;
0316   std::vector<IndexType> m_previous;
0317   std::vector<IndexType> m_next;
0318   std::vector<typename detail_lt::Types<eBoundSize>::Coefficients> m_params;
0319   std::vector<typename detail_lt::Types<eBoundSize>::Covariance> m_cov;
0320
0321   std::vector<double> m_meas;
0322   std::vector<MultiTrajectoryTraits::IndexType> m_measOffset;
0323   std::vector<double> m_measCov;
0324   std::vector<MultiTrajectoryTraits::IndexType> m_measCovOffset;
0325
0326   std::vector<typename detail_lt::Types<eBoundSize>::Covariance> m_jac;
0327   std::vector<std::optional<SourceLink>> m_sourceLinks;
0328   std::vector<ProjectorBitset> m_projectors;
0329
0330   // owning vector of shared pointers to surfaces
0331   //
0332   // This might be problematic when appending a large number of surfaces
0333   // trackstates, because vector has to reallocated and thus copy. This might
0334   // be handled in a smart way by moving but not sure.
0335   std::vector<std::shared_ptr<const Surface>> m_referenceSurfaces;
0336
0337   std::vector<HashedString> m_dynamicKeys;
0338   std::unordered_map<HashedString, std::unique_ptr<detail::DynamicColumnBase>>
0339       m_dynamic;
0340 };
0341
0342 }  // namespace detail_vmt
0343
0344 class VectorMultiTrajectory;
0345
0346 template <>
0347 struct IsReadOnlyMultiTrajectory<VectorMultiTrajectory> : std::false_type {};
0348
0349 class VectorMultiTrajectory final
0350     : public detail_vmt::VectorMultiTrajectoryBase,
0351       public MultiTrajectory<VectorMultiTrajectory> {
0352 #ifndef DOXYGEN
0353   friend MultiTrajectory<VectorMultiTrajectory>;
0354 #endif
0355
0356  public:
0357   VectorMultiTrajectory() = default;
0358   VectorMultiTrajectory(const VectorMultiTrajectory& other)
0359       : VectorMultiTrajectoryBase{other} {}
0360
0361   VectorMultiTrajectory(VectorMultiTrajectory&& other)
0362       : VectorMultiTrajectoryBase{std::move(other)} {}
0363
0364   Statistics statistics() const {
0365     return detail_vmt::VectorMultiTrajectoryBase::statistics(*this);
0366   }
0367
0368   // BEGIN INTERFACE
0369   TrackStateProxy::Parameters parameters_impl(IndexType parIdx) {
0370     return TrackStateProxy::Parameters{m_params[parIdx].data()};
0371   }
0372
0373   ConstTrackStateProxy::Parameters parameters_impl(IndexType parIdx) const {
0374     return ConstTrackStateProxy::Parameters{m_params[parIdx].data()};
0375   }
0376
0377   TrackStateProxy::Covariance covariance_impl(IndexType parIdx) {
0378     return TrackStateProxy::Covariance{m_cov[parIdx].data()};
0379   }
0380
0381   ConstTrackStateProxy::Covariance covariance_impl(IndexType parIdx) const {
0382     return ConstTrackStateProxy::Covariance{m_cov[parIdx].data()};
0383   }
0384
0385   TrackStateProxy::Covariance jacobian_impl(IndexType istate) {
0386     IndexType jacIdx = m_index[istate].ijacobian;
0387     return TrackStateProxy::Covariance{m_jac[jacIdx].data()};
0388   }
0389
0390   ConstTrackStateProxy::Covariance jacobian_impl(IndexType istate) const {
0391     IndexType jacIdx = m_index[istate].ijacobian;
0392     return ConstTrackStateProxy::Covariance{m_jac[jacIdx].data()};
0393   }
0394
0395   template <std::size_t measdim>
0396   TrackStateProxy::Measurement<measdim> measurement_impl(IndexType istate) {
0397     IndexType offset = m_measOffset[istate];
0398     return TrackStateProxy::Measurement<measdim>{&m_meas[offset]};
0399   }
0400
0401   template <std::size_t measdim>
0402   ConstTrackStateProxy::Measurement<measdim> measurement_impl(
0403       IndexType istate) const {
0404     IndexType offset = m_measOffset[istate];
0405     return ConstTrackStateProxy::Measurement<measdim>{&m_meas[offset]};
0406   }
0407
0408   template <std::size_t measdim>
0409   TrackStateProxy::MeasurementCovariance<measdim> measurementCovariance_impl(
0410       IndexType istate) {
0411     IndexType offset = m_measCovOffset[istate];
0412     return TrackStateProxy::MeasurementCovariance<measdim>{&m_measCov[offset]};
0413   }
0414
0415   template <std::size_t measdim>
0416   ConstTrackStateProxy::MeasurementCovariance<measdim>
0417   measurementCovariance_impl(IndexType istate) const {
0418     IndexType offset = m_measCovOffset[istate];
0419     return ConstTrackStateProxy::MeasurementCovariance<measdim>{
0420         &m_measCov[offset]};
0421   }
0422
0423   IndexType addTrackState_impl(
0424       TrackStatePropMask mask = TrackStatePropMask::All,
0425       IndexType iprevious = kInvalid);
0426
0427   void addTrackStateComponents_impl(IndexType istate, TrackStatePropMask mask);
0428
0429   void reserve(std::size_t n);
0430
0431   void shareFrom_impl(IndexType iself, IndexType iother,
0432                       TrackStatePropMask shareSource,
0433                       TrackStatePropMask shareTarget);
0434
0435   void unset_impl(TrackStatePropMask target, IndexType istate);
0436
0437   bool has_impl(HashedString key, IndexType istate) const {
0438     return detail_vmt::VectorMultiTrajectoryBase::has_impl(*this, key, istate);
0439   }
0440
0441   IndexType size_impl() const {
0442     return m_index.size();
0443   }
0444
0445   void clear_impl();
0446
0447   std::any component_impl(HashedString key, IndexType istate) {
0448     return detail_vmt::VectorMultiTrajectoryBase::component_impl<false>(
0449         *this, key, istate);
0450   }
0451
0452   std::any component_impl(HashedString key, IndexType istate) const {
0453     return detail_vmt::VectorMultiTrajectoryBase::component_impl<true>(
0454         *this, key, istate);
0455   }
0456
0457   template <typename T>
0458   void addColumn_impl(const std::string& key) {
0459     Acts::HashedString hashedKey = hashString(key);
0460     m_dynamic.insert({hashedKey, std::make_unique<detail::DynamicColumn<T>>()});
0461   }
0462
0463   bool hasColumn_impl(HashedString key) const {
0464     return detail_vmt::VectorMultiTrajectoryBase::hasColumn_impl(*this, key);
0465   }
0466
0467   void allocateCalibrated_impl(IndexType istate, std::size_t measdim) {
0468     throw_assert(measdim > 0 && measdim <= eBoundSize,
0469                  "Invalid measurement dimension detected");
0470
0471     if (m_measOffset[istate] != kInvalid &&
0472         m_measCovOffset[istate] != kInvalid &&
0473         m_index[istate].measdim == measdim) {
0474       return;
0475     }
0476
0477     m_index[istate].measdim = measdim;
0478
0479     m_measOffset[istate] = static_cast<IndexType>(m_meas.size());
0480     m_meas.resize(m_meas.size() + measdim);
0481
0482     m_measCovOffset[istate] = static_cast<IndexType>(m_measCov.size());
0483     m_measCov.resize(m_measCov.size() + measdim * measdim);
0484   }
0485
0486   void setUncalibratedSourceLink_impl(IndexType istate, SourceLink sourceLink) {
0487     m_sourceLinks[m_index[istate].iuncalibrated] = std::move(sourceLink);
0488   }
0489
0490   void setReferenceSurface_impl(IndexType istate,
0491                                 std::shared_ptr<const Surface> surface) {
0492     m_referenceSurfaces[istate] = std::move(surface);
0493   }
0494
0495   void copyDynamicFrom_impl(IndexType dstIdx, HashedString key,
0496                             const std::any& srcPtr);
0497
0498   // END INTERFACE
0499 };
0500
0501 ACTS_STATIC_CHECK_CONCEPT(MutableMultiTrajectoryBackend, VectorMultiTrajectory);
0502
0503 class ConstVectorMultiTrajectory;
0504
0505 template <>
0506 struct IsReadOnlyMultiTrajectory<ConstVectorMultiTrajectory> : std::true_type {
0507 };
0508
0509 class ConstVectorMultiTrajectory final
0510     : public detail_vmt::VectorMultiTrajectoryBase,
0511       public MultiTrajectory<ConstVectorMultiTrajectory> {
0512 #ifndef DOXYGEN
0513   friend MultiTrajectory<ConstVectorMultiTrajectory>;
0514 #endif
0515
0516  public:
0517   ConstVectorMultiTrajectory() = default;
0518
0519   ConstVectorMultiTrajectory(const ConstVectorMultiTrajectory& other)
0520       : VectorMultiTrajectoryBase{other} {}
0521
0522   ConstVectorMultiTrajectory(const VectorMultiTrajectory& other)
0523       : VectorMultiTrajectoryBase{other} {}
0524
0525   ConstVectorMultiTrajectory(VectorMultiTrajectory&& other)
0526       : VectorMultiTrajectoryBase{std::move(other)} {}
0527
0528   ConstVectorMultiTrajectory(ConstVectorMultiTrajectory&&) = default;
0529
0530   Statistics statistics() const {
0531     return detail_vmt::VectorMultiTrajectoryBase::statistics(*this);
0532   }
0533
0534   // BEGIN INTERFACE
0535
0536   ConstTrackStateProxy::Parameters parameters_impl(IndexType parIdx) const {
0537     return ConstTrackStateProxy::Parameters{m_params[parIdx].data()};
0538   }
0539
0540   ConstTrackStateProxy::Covariance covariance_impl(IndexType parIdx) const {
0541     return ConstTrackStateProxy::Covariance{m_cov[parIdx].data()};
0542   }
0543
0544   ConstTrackStateProxy::Covariance jacobian_impl(IndexType istate) const {
0545     IndexType jacIdx = m_index[istate].ijacobian;
0546     return ConstTrackStateProxy::Covariance{m_jac[jacIdx].data()};
0547   }
0548
0549   template <std::size_t measdim>
0550   ConstTrackStateProxy::Measurement<measdim> measurement_impl(
0551       IndexType istate) const {
0552     IndexType offset = m_measOffset[istate];
0553     return ConstTrackStateProxy::Measurement<measdim>{&m_meas[offset]};
0554   }
0555
0556   template <std::size_t measdim>
0557   ConstTrackStateProxy::MeasurementCovariance<measdim>
0558   measurementCovariance_impl(IndexType istate) const {
0559     IndexType offset = m_measCovOffset[istate];
0560     return ConstTrackStateProxy::MeasurementCovariance<measdim>{
0561         &m_measCov[offset]};
0562   }
0563
0564   bool has_impl(HashedString key, IndexType istate) const {
0565     return detail_vmt::VectorMultiTrajectoryBase::has_impl(*this, key, istate);
0566   }
0567
0568   IndexType size_impl() const {
0569     return m_index.size();
0570   }
0571
0572   std::any component_impl(HashedString key, IndexType istate) const {
0573     return detail_vmt::VectorMultiTrajectoryBase::component_impl<true>(
0574         *this, key, istate);
0575   }
0576
0577   bool hasColumn_impl(HashedString key) const {
0578     return detail_vmt::VectorMultiTrajectoryBase::hasColumn_impl(*this, key);
0579   }
0580
0581   // END INTERFACE
0582 };
0583
0584 ACTS_STATIC_CHECK_CONCEPT(ConstMultiTrajectoryBackend,
0585                           ConstVectorMultiTrajectory);
0586
0587 }  // namespace Acts