TrackFitting/detail/GsfUtils.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/Direction.hpp"
0012 #include "Acts/Definitions/TrackParametrization.hpp"
0013 #include "Acts/EventData/BoundTrackParameters.hpp"
0014 #include "Acts/EventData/MultiTrajectoryHelpers.hpp"
0015 #include "Acts/EventData/Types.hpp"
0016 #include "Acts/Propagator/detail/PointwiseMaterialInteraction.hpp"
0017 #include "Acts/Surfaces/Surface.hpp"
0018 #include "Acts/TrackFitting/BetheHeitlerApprox.hpp"
0019 #include "Acts/TrackFitting/GsfComponent.hpp"
0020 #include "Acts/Utilities/AlgebraHelpers.hpp"
0021 #include "Acts/Utilities/Intersection.hpp"
0022 #include "Acts/Utilities/Logger.hpp"
0023 #include "Acts/Utilities/Zip.hpp"
0024
0025 #include <algorithm>
0026 #include <array>
0027 #include <cassert>
0028 #include <cmath>
0029 #include <cstddef>
0030 #include <iomanip>
0031 #include <map>
0032 #include <ostream>
0033 #include <tuple>
0034 #include <vector>
0035
0036 namespace Acts::detail::Gsf {
0037
0038 /// The tolerated difference to 1 to accept weights as normalized
0039 constexpr static double s_normalizationTolerance = 1.e-4;
0040
0041 template <typename component_range_t, typename projector_t>
0042 bool weightsAreNormalized(const component_range_t &cmps,
0043                           const projector_t &proj,
0044                           double tol = s_normalizationTolerance) {
0045   double sumOfWeights = 0.0;
0046
0047   for (auto &&cmp : cmps) {
0048     sumOfWeights += proj(cmp);
0049   }
0050
0051   return std::abs(sumOfWeights - 1.0) < tol;
0052 }
0053
0054 template <typename component_range_t, typename projector_t>
0055 void normalizeWeights(component_range_t &cmps, const projector_t &proj) {
0056   double sumOfWeights = 0.0;
0057
0058   // we need decltype(auto) here to support proxy-types with reference
0059   // semantics, otherwise there is a `cannot bind ... to ...` error
0060   for (auto &&cmp : cmps) {
0061     assert(std::isfinite(proj(cmp)) && "weight not finite in normalization");
0062     sumOfWeights += proj(cmp);
0063   }
0064
0065   assert(sumOfWeights > 0 && "sum of weights is not > 0");
0066
0067   for (auto &&cmp : cmps) {
0068     proj(cmp) /= sumOfWeights;
0069   }
0070 }
0071
0072 // A class that prints information about the state on construction and
0073 // destruction, it also contains some assertions in the constructor and
0074 // destructor. It can be removed without change of behaviour, since it only
0075 // holds const references
0076 template <typename propagator_state_t, typename stepper_t, typename navigator_t>
0077 class ScopedGsfInfoPrinterAndChecker {
0078   const propagator_state_t &m_state;
0079   const stepper_t &m_stepper;
0080   const navigator_t &m_navigator;
0081   double m_p_initial;
0082   const Logger &m_logger;
0083
0084   const Logger &logger() const { return m_logger; }
0085
0086   void print_component_stats() const {
0087     std::size_t i = 0;
0088     for (auto cmp : m_stepper.constComponentIterable(m_state.stepping)) {
0089       const auto &singleStepper = cmp.singleStepper(m_stepper);
0090       ACTS_VERBOSE(
0091           "  #" << i++
0092                 << " pos: " << singleStepper.position(cmp.state()).transpose()
0093                 << ", dir: " << singleStepper.direction(cmp.state()).transpose()
0094                 << ", weight: " << cmp.weight() << ", status: " << cmp.status()
0095                 << ", qop: " << singleStepper.qOverP(cmp.state())
0096                 << ", det(cov): "
0097                 << singleStepper.covariance(cmp.state()).value().determinant());
0098     }
0099   }
0100
0101   void checks(bool onStart) const {
0102     const auto cmps = m_stepper.constComponentIterable(m_state.stepping);
0103     [[maybe_unused]] const bool allFinite =
0104         std::all_of(cmps.begin(), cmps.end(),
0105                     [](auto cmp) { return std::isfinite(cmp.weight()); });
0106     [[maybe_unused]] const bool allNormalized = weightsAreNormalized(
0107         cmps, [](const auto &cmp) { return cmp.weight(); });
0108     [[maybe_unused]] const bool zeroComponents =
0109         m_stepper.numberComponents(m_state.stepping) == 0;
0110
0111     if (onStart) {
0112       assert(!zeroComponents && "no cmps at the start");
0113       assert(allFinite && "weights not finite at the start");
0114       assert(allNormalized && "not normalized at the start");
0115     } else {
0116       assert(!zeroComponents && "no cmps at the end");
0117       assert(allFinite && "weights not finite at the end");
0118       assert(allNormalized && "not normalized at the end");
0119     }
0120   }
0121
0122  public:
0123   ScopedGsfInfoPrinterAndChecker(const propagator_state_t &state,
0124                                  const stepper_t &stepper,
0125                                  const navigator_t &navigator,
0126                                  const Logger &logger)
0127       : m_state(state),
0128         m_stepper(stepper),
0129         m_navigator(navigator),
0130         m_p_initial(stepper.absoluteMomentum(state.stepping)),
0131         m_logger{logger} {
0132     // Some initial printing
0133     checks(true);
0134     ACTS_VERBOSE("Gsf step "
0135                  << state.steps << " at mean position "
0136                  << stepper.position(state.stepping).transpose()
0137                  << " with direction "
0138                  << stepper.direction(state.stepping).transpose()
0139                  << " and momentum " << stepper.absoluteMomentum(state.stepping)
0140                  << " and charge " << stepper.charge(state.stepping));
0141     ACTS_VERBOSE("Propagation is in " << state.options.direction << " mode");
0142     print_component_stats();
0143   }
0144
0145   ~ScopedGsfInfoPrinterAndChecker() {
0146     if (m_navigator.currentSurface(m_state.navigation)) {
0147       const auto p_final = m_stepper.absoluteMomentum(m_state.stepping);
0148       ACTS_VERBOSE("Component status at end of step:");
0149       print_component_stats();
0150       ACTS_VERBOSE("Delta Momentum = " << std::setprecision(5)
0151                                        << p_final - m_p_initial);
0152     }
0153     checks(false);
0154   }
0155 };
0156
0157 double calculateDeterminant(
0158     const double *fullCalibratedCovariance,
0159     TrackStateTraits<kMeasurementSizeMax, true>::Covariance predictedCovariance,
0160     BoundSubspaceIndices projector, unsigned int calibratedSize);
0161
0162 /// Reweight the components according to `R. Frühwirth, "Track fitting
0163 /// with non-Gaussian noise"`. See also the implementation in Athena at
0164 /// PosteriorWeightsCalculator.cxx
0165 /// @note The weights are not renormalized!
0166 template <typename traj_t>
0167 void computePosteriorWeights(const traj_t &mt,
0168                              const std::vector<TrackIndexType> &tips,
0169                              std::map<TrackIndexType, double> &weights) {
0170   // Helper Function to compute detR
0171
0172   // Find minChi2, this can be used to factor some things later in the
0173   // exponentiation
0174   const auto minChi2 =
0175       mt.getTrackState(
0176             *std::ranges::min_element(tips,
0177                                       [&](const auto &a, const auto &b) {
0178                                         return mt.getTrackState(a).chi2() <
0179                                                mt.getTrackState(b).chi2();
0180                                       }))
0181           .chi2();
0182
0183   //! [posterior weights]
0184   // Loop over the tips and compute new weights
0185   for (auto tip : tips) {
0186     const auto state = mt.getTrackState(tip);
0187     const double chi2 = state.chi2() - minChi2;
0188     const double detR = calculateDeterminant(
0189         state.effectiveCalibratedCovariance().data(),
0190         state.predictedCovariance(), state.projectorSubspaceIndices(),
0191         state.calibratedSize());
0192
0193     if (detR <= 0) {
0194       // If the determinant is not positive, just leave the weight as it is
0195       continue;
0196     }
0197
0198     const double factor = std::sqrt(1. / detR) * safeExp(-0.5 * chi2);
0199
0200     if (!std::isfinite(factor)) {
0201       // If something is not finite here, just leave the weight as it is
0202       continue;
0203     }
0204
0205     weights.at(tip) *= factor;
0206   }
0207   //! [posterior weights]
0208 }
0209
0210 /// Enumeration type to allow templating on the state we want to project on with
0211 /// a MultiTrajectory
0212 enum class StatesType { ePredicted, eFiltered, eSmoothed };
0213
0214 inline std::ostream &operator<<(std::ostream &os, StatesType type) {
0215   constexpr static std::array names = {"predicted", "filtered", "smoothed"};
0216   os << names[static_cast<int>(type)];
0217   return os;
0218 }
0219
0220 /// @brief Projector type which maps a MultiTrajectory-Index to a tuple of
0221 /// [weight, parameters, covariance]. Therefore, it contains a MultiTrajectory
0222 /// and for now a std::map for the weights
0223 template <StatesType type, typename traj_t>
0224 struct MultiTrajectoryProjector {
0225   const traj_t &mt;
0226   const std::map<TrackIndexType, double> &weights;
0227
0228   auto operator()(TrackIndexType idx) const {
0229     const auto proxy = mt.getTrackState(idx);
0230     switch (type) {
0231       case StatesType::ePredicted:
0232         return std::tuple(weights.at(idx), proxy.predicted(),
0233                           proxy.predictedCovariance());
0234       case StatesType::eFiltered:
0235         return std::tuple(weights.at(idx), proxy.filtered(),
0236                           proxy.filteredCovariance());
0237       case StatesType::eSmoothed:
0238         return std::tuple(weights.at(idx), proxy.smoothed(),
0239                           proxy.smoothedCovariance());
0240       default:
0241         throw std::invalid_argument(
0242             "Incorrect StatesType, should be ePredicted"
0243             ", eFiltered, or eSmoothed.");
0244     }
0245   }
0246 };
0247
0248 /// Small Helper class that allows to carry a temporary value until we decide to
0249 /// update the actual value. The temporary value is deliberately only accessible
0250 /// with a mutable reference
0251 template <typename T>
0252 class Updatable {
0253   T m_tmp{};
0254   T m_val{};
0255
0256  public:
0257   Updatable() : m_tmp(0), m_val(0) {}
0258
0259   T &tmp() { return m_tmp; }
0260   void update() { m_val = m_tmp; }
0261
0262   const T &val() const { return m_val; }
0263 };
0264
0265 /// Remove components with low weights and renormalize from the component
0266 /// cache
0267 /// TODO This function does not expect normalized components, but this
0268 /// could be redundant work...
0269 void removeLowWeightComponents(std::vector<GsfComponent> &cmps,
0270                                double weightCutoff);
0271
0272 template <typename traj_t>
0273 struct TemporaryStates {
0274   traj_t traj;
0275   std::vector<TrackIndexType> tips;
0276   std::map<TrackIndexType, double> weights;
0277
0278   void clear() {
0279     traj.clear();
0280     tips.clear();
0281     weights.clear();
0282   }
0283 };
0284
0285 /// Function that updates the stepper from the MultiTrajectory
0286 template <typename traj_t, typename propagator_state_t, typename stepper_t>
0287 void updateStepper(propagator_state_t &state, const stepper_t &stepper,
0288                    const TemporaryStates<traj_t> &tmpStates,
0289                    double weightCutoff) {
0290   auto cmps = stepper.componentIterable(state.stepping);
0291   for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
0292     // we set ignored components to missed, so we can remove them after
0293     // the loop
0294     if (tmpStates.weights.at(idx) < weightCutoff) {
0295       cmp.status() = IntersectionStatus::unreachable;
0296       continue;
0297     }
0298
0299     auto proxy = tmpStates.traj.getTrackState(idx);
0300
0301     // The same update the KF and the single-stepper CKF apply after a filter
0302     // step, so a component ends up in the same state a single stepper would
0303     cmp.update(
0304         stepper, MultiTrajectoryHelpers::freeFiltered(state.geoContext, proxy),
0305         proxy.filtered(), proxy.filteredCovariance(), proxy.referenceSurface());
0306     cmp.weight() = tmpStates.weights.at(idx);
0307   }
0308
0309   stepper.removeMissedComponents(state.stepping);
0310
0311   // TODO we have two normalization passes here now, this can probably be
0312   // optimized
0313   detail::Gsf::normalizeWeights(
0314       cmps, [&](auto cmp) -> double & { return cmp.weight(); });
0315 }
0316
0317 /// Function that updates the stepper from the ComponentCache
0318 template <typename propagator_state_t, typename stepper_t>
0319 void updateStepper(propagator_state_t &state, const stepper_t &stepper,
0320                    const Surface &surface,
0321                    const std::vector<GsfComponent> &componentCache) {
0322   // Clear components before adding new ones
0323   stepper.clearComponents(state.stepping);
0324
0325   // Finally loop over components
0326   for (const auto &[weight, pars, cov] : componentCache) {
0327     // Add the component to the stepper
0328     BoundTrackParameters bound(surface.getSharedPtr(), pars, cov,
0329                                stepper.particleHypothesis(state.stepping));
0330
0331     stepper.addComponent(state.stepping, std::move(bound), weight);
0332   }
0333 }
0334
0335 double applyBetheHeitler(
0336     const GeometryContext &geoContext, const Surface &surface,
0337     Direction direction, const BoundTrackParameters &initialParameters,
0338     double initialWeight, const BetheHeitlerApprox &betheHeitlerApprox,
0339     std::vector<BetheHeitlerApprox::Component> &betheHeitlerCache,
0340     double weightCutoff, std::vector<GsfComponent> &componentCache,
0341     std::size_t &nInvalidBetheHeitler, double &maxPathXOverX0,
0342     const Logger &logger);
0343
0344 /// Apply the multiple scattering to the state
0345 template <typename propagator_state_t, typename stepper_t>
0346 Result<void> applyMultipleScattering(propagator_state_t &state,
0347                                      const stepper_t &stepper,
0348                                      const Surface &surface,
0349                                      const MaterialUpdateMode &updateMode,
0350                                      const Logger &logger) {
0351   for (auto cmp : stepper.componentIterable(state.stepping)) {
0352     auto singleState = cmp.singleState(state);
0353     const auto &singleStepper = cmp.singleStepper(stepper);
0354
0355     const Result<detail::PointwiseMaterialEffects> materialInteractionRes =
0356         detail::performMaterialInteraction(
0357             singleState, singleStepper, surface, updateMode,
0358             NoiseUpdateMode::addNoise, true, false, logger);
0359     if (!materialInteractionRes.ok()) {
0360       ACTS_DEBUG("Error performing material interaction: "
0361                  << materialInteractionRes.error());
0362       return materialInteractionRes.error();
0363     }
0364
0365     assert(singleStepper.covariance(singleState.stepping)
0366                .value()
0367                .array()
0368                .isFinite()
0369                .all() &&
0370            "covariance not finite after multi scattering");
0371   }
0372
0373   return Result<void>::success();
0374 }
0375
0376 }  // namespace Acts::detail::Gsf