src/EventData/VectorMultiTrajectory.cpp

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