src/EventData/VectorMultiTrajectory.cpp

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 #include "Acts/EventData/VectorMultiTrajectory.hpp"
0010
0011 #include "Acts/EventData/MultiTrajectory.hpp"
0012 #include "Acts/EventData/TrackStatePropMask.hpp"
0013 #include "Acts/Utilities/Helpers.hpp"
0014
0015 #include <iomanip>
0016 #include <ostream>
0017 #include <type_traits>
0018
0019 #include <boost/histogram.hpp>
0020 #include <boost/histogram/axis/category.hpp>
0021 #include <boost/histogram/make_histogram.hpp>
0022
0023 namespace Acts {
0024
0025 auto VectorMultiTrajectory::addTrackState_impl(TrackStatePropMask mask,
0026                                                IndexType iprevious)
0027     -> IndexType {
0028   using PropMask = TrackStatePropMask;
0029
0030   m_index.emplace_back();
0031   IndexData& p = m_index.back();
0032   IndexType index = m_index.size() - 1;
0033   m_previous.emplace_back(iprevious);
0034   m_next.emplace_back(kInvalid);
0035
0036   p.allocMask = mask;
0037
0038   // always set, but can be null
0039   m_referenceSurfaces.emplace_back(nullptr);
0040
0041   assert(m_params.size() == m_cov.size());
0042
0043   if (ACTS_CHECK_BIT(mask, PropMask::Predicted)) {
0044     m_params.emplace_back();
0045     m_cov.emplace_back();
0046     p.ipredicted = m_params.size() - 1;
0047   }
0048
0049   if (ACTS_CHECK_BIT(mask, PropMask::Filtered)) {
0050     m_params.emplace_back();
0051     m_cov.emplace_back();
0052     p.ifiltered = m_params.size() - 1;
0053   }
0054
0055   if (ACTS_CHECK_BIT(mask, PropMask::Smoothed)) {
0056     m_params.emplace_back();
0057     m_cov.emplace_back();
0058     p.ismoothed = m_params.size() - 1;
0059   }
0060
0061   assert(m_params.size() == m_cov.size());
0062
0063   if (ACTS_CHECK_BIT(mask, PropMask::Jacobian)) {
0064     m_jac.emplace_back();
0065     p.ijacobian = m_jac.size() - 1;
0066   }
0067
0068   m_sourceLinks.emplace_back(std::nullopt);
0069   p.iuncalibrated = m_sourceLinks.size() - 1;
0070
0071   m_measOffset.push_back(kInvalid);
0072   m_measCovOffset.push_back(kInvalid);
0073
0074   if (ACTS_CHECK_BIT(mask, PropMask::Calibrated)) {
0075     m_sourceLinks.emplace_back(std::nullopt);
0076     p.icalibratedsourcelink = m_sourceLinks.size() - 1;
0077
0078     m_projectors.emplace_back();
0079     p.iprojector = m_projectors.size() - 1;
0080   }
0081
0082   // dynamic columns
0083   for (auto& [key, vec] : m_dynamic) {
0084     vec->add();
0085   }
0086
0087   return index;
0088 }
0089
0090 void VectorMultiTrajectory::addTrackStateComponents_impl(
0091     IndexType istate, TrackStatePropMask mask) {
0092   using PropMask = TrackStatePropMask;
0093
0094   IndexData& p = m_index[istate];
0095   PropMask currentMask = p.allocMask;
0096
0097   assert(m_params.size() == m_cov.size());
0098
0099   if (ACTS_CHECK_BIT(mask, PropMask::Predicted) &&
0100       !ACTS_CHECK_BIT(currentMask, PropMask::Predicted)) {
0101     m_params.emplace_back();
0102     m_cov.emplace_back();
0103     p.ipredicted = m_params.size() - 1;
0104   }
0105
0106   if (ACTS_CHECK_BIT(mask, PropMask::Filtered) &&
0107       !ACTS_CHECK_BIT(currentMask, PropMask::Filtered)) {
0108     m_params.emplace_back();
0109     m_cov.emplace_back();
0110     p.ifiltered = m_params.size() - 1;
0111   }
0112
0113   if (ACTS_CHECK_BIT(mask, PropMask::Smoothed) &&
0114       !ACTS_CHECK_BIT(currentMask, PropMask::Smoothed)) {
0115     m_params.emplace_back();
0116     m_cov.emplace_back();
0117     p.ismoothed = m_params.size() - 1;
0118   }
0119
0120   assert(m_params.size() == m_cov.size());
0121
0122   if (ACTS_CHECK_BIT(mask, PropMask::Jacobian) &&
0123       !ACTS_CHECK_BIT(currentMask, PropMask::Jacobian)) {
0124     m_jac.emplace_back();
0125     p.ijacobian = m_jac.size() - 1;
0126   }
0127
0128   if (ACTS_CHECK_BIT(mask, PropMask::Calibrated) &&
0129       !ACTS_CHECK_BIT(currentMask, PropMask::Calibrated)) {
0130     m_sourceLinks.emplace_back(std::nullopt);
0131     p.icalibratedsourcelink = m_sourceLinks.size() - 1;
0132
0133     m_projectors.emplace_back();
0134     p.iprojector = m_projectors.size() - 1;
0135   }
0136
0137   p.allocMask |= mask;
0138 }
0139
0140 void VectorMultiTrajectory::shareFrom_impl(IndexType iself, IndexType iother,
0141                                            TrackStatePropMask shareSource,
0142                                            TrackStatePropMask shareTarget) {
0143   // auto other = getTrackState(iother);
0144   IndexData& self = m_index[iself];
0145   IndexData& other = m_index[iother];
0146
0147   assert(ACTS_CHECK_BIT(getTrackState(iother).getMask(), shareSource) &&
0148          "Source has incompatible allocation");
0149
0150   using PM = TrackStatePropMask;
0151
0152   IndexType sourceIndex{kInvalid};
0153   switch (shareSource) {
0154     case PM::Predicted:
0155       sourceIndex = other.ipredicted;
0156       break;
0157     case PM::Filtered:
0158       sourceIndex = other.ifiltered;
0159       break;
0160     case PM::Smoothed:
0161       sourceIndex = other.ismoothed;
0162       break;
0163     case PM::Jacobian:
0164       sourceIndex = other.ijacobian;
0165       break;
0166     default:
0167       throw std::domain_error{"Unable to share this component"};
0168   }
0169
0170   assert(sourceIndex != kInvalid);
0171
0172   switch (shareTarget) {
0173     case PM::Predicted:
0174       assert(shareSource != PM::Jacobian);
0175       self.ipredicted = sourceIndex;
0176       break;
0177     case PM::Filtered:
0178       assert(shareSource != PM::Jacobian);
0179       self.ifiltered = sourceIndex;
0180       break;
0181     case PM::Smoothed:
0182       assert(shareSource != PM::Jacobian);
0183       self.ismoothed = sourceIndex;
0184       break;
0185     case PM::Jacobian:
0186       assert(shareSource == PM::Jacobian);
0187       self.ijacobian = sourceIndex;
0188       break;
0189     default:
0190       throw std::domain_error{"Unable to share this component"};
0191   }
0192 }
0193
0194 void VectorMultiTrajectory::unset_impl(TrackStatePropMask target,
0195                                        IndexType istate) {
0196   using PM = TrackStatePropMask;
0197
0198   switch (target) {
0199     case PM::Predicted:
0200       m_index[istate].ipredicted = kInvalid;
0201       break;
0202     case PM::Filtered:
0203       m_index[istate].ifiltered = kInvalid;
0204       break;
0205     case PM::Smoothed:
0206       m_index[istate].ismoothed = kInvalid;
0207       break;
0208     case PM::Jacobian:
0209       m_index[istate].ijacobian = kInvalid;
0210       break;
0211     case PM::Calibrated:
0212       m_measOffset[istate] = kInvalid;
0213       m_measCovOffset[istate] = kInvalid;
0214       break;
0215     default:
0216       throw std::domain_error{"Unable to unset this component"};
0217   }
0218 }
0219
0220 void VectorMultiTrajectory::clear_impl() {
0221   m_index.clear();
0222   m_previous.clear();
0223   m_next.clear();
0224   m_params.clear();
0225   m_cov.clear();
0226   m_meas.clear();
0227   m_measOffset.clear();
0228   m_measCov.clear();
0229   m_measCovOffset.clear();
0230   m_jac.clear();
0231   m_sourceLinks.clear();
0232   m_projectors.clear();
0233   m_referenceSurfaces.clear();
0234   for (auto& [key, vec] : m_dynamic) {
0235     vec->clear();
0236   }
0237 }
0238
0239 void detail_vmt::VectorMultiTrajectoryBase::Statistics::toStream(
0240     std::ostream& os, std::size_t n) {
0241   using namespace boost::histogram;
0242   using cat = axis::category<std::string>;
0243
0244   auto& h = hist;
0245
0246   auto column_axis = axis::get<cat>(h.axis(0));
0247   auto type_axis = axis::get<axis::category<>>(h.axis(1));
0248
0249   auto p = [&](const auto& key, const auto v, const std::string suffix = "") {
0250     os << std::setw(20) << key << ": ";
0251     if constexpr (std::is_same_v<std::decay_t<decltype(v)>, double>) {
0252       os << std::fixed << std::setw(8) << std::setprecision(2) << v / n
0253          << suffix;
0254     } else {
0255       os << std::fixed << std::setw(8) << double(v) / n << suffix;
0256     }
0257     os << std::endl;
0258   };
0259
0260   for (int t = 0; t < type_axis.size(); t++) {
0261     os << (type_axis.bin(t) == 1 ? "meas" : "other") << ":" << std::endl;
0262     double total = 0;
0263     for (int c = 0; c < column_axis.size(); c++) {
0264       std::string key = column_axis.bin(c);
0265       auto v = h.at(c, t);
0266       if (key == "count") {
0267         p(key, static_cast<std::size_t>(v));
0268         continue;
0269       }
0270       p(key, v / 1024 / 1024, "M");
0271       total += v;
0272     }
0273     p("total", total / 1024 / 1024, "M");
0274   }
0275 }
0276
0277 void VectorMultiTrajectory::reserve(std::size_t n) {
0278   m_index.reserve(n);
0279   m_previous.reserve(n);
0280   m_next.reserve(n);
0281   m_params.reserve(n * 2);
0282   m_cov.reserve(n * 2);
0283   m_meas.reserve(n * 2);
0284   m_measOffset.reserve(n);
0285   m_measCov.reserve(n * 2 * 2);
0286   m_measCovOffset.reserve(n);
0287   m_jac.reserve(n);
0288   m_sourceLinks.reserve(n);
0289   m_projectors.reserve(n);
0290   m_referenceSurfaces.reserve(n);
0291
0292   for (auto& [key, vec] : m_dynamic) {
0293     vec->reserve(n);
0294   }
0295 }
0296
0297 void VectorMultiTrajectory::copyDynamicFrom_impl(IndexType dstIdx,
0298                                                  HashedString key,
0299                                                  const std::any& srcPtr) {
0300   auto it = m_dynamic.find(key);
0301   if (it == m_dynamic.end()) {
0302     throw std::invalid_argument{
0303         "Destination container does not have matching dynamic column"};
0304   }
0305
0306   it->second->copyFrom(dstIdx, srcPtr);
0307 }
0308
0309 }  // namespace Acts