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       auto getVector = [&](auto idx) {
0090         return cmp.pars().template segment<3>(idx).transpose();
0091       };
0092       ACTS_VERBOSE("  #" << i++ << " pos: " << getVector(eFreePos0) << ", dir: "
0093                          << getVector(eFreeDir0) << ", weight: " << cmp.weight()
0094                          << ", status: " << cmp.status()
0095                          << ", qop: " << cmp.pars()[eFreeQOverP]
0096                          << ", det(cov): " << cmp.cov().determinant());
0097     }
0098   }
0099
0100   void checks(bool onStart) const {
0101     const auto cmps = m_stepper.constComponentIterable(m_state.stepping);
0102     [[maybe_unused]] const bool allFinite =
0103         std::all_of(cmps.begin(), cmps.end(),
0104                     [](auto cmp) { return std::isfinite(cmp.weight()); });
0105     [[maybe_unused]] const bool allNormalized = weightsAreNormalized(
0106         cmps, [](const auto &cmp) { return cmp.weight(); });
0107     [[maybe_unused]] const bool zeroComponents =
0108         m_stepper.numberComponents(m_state.stepping) == 0;
0109
0110     if (onStart) {
0111       assert(!zeroComponents && "no cmps at the start");
0112       assert(allFinite && "weights not finite at the start");
0113       assert(allNormalized && "not normalized at the start");
0114     } else {
0115       assert(!zeroComponents && "no cmps at the end");
0116       assert(allFinite && "weights not finite at the end");
0117       assert(allNormalized && "not normalized at the end");
0118     }
0119   }
0120
0121  public:
0122   ScopedGsfInfoPrinterAndChecker(const propagator_state_t &state,
0123                                  const stepper_t &stepper,
0124                                  const navigator_t &navigator,
0125                                  const Logger &logger)
0126       : m_state(state),
0127         m_stepper(stepper),
0128         m_navigator(navigator),
0129         m_p_initial(stepper.absoluteMomentum(state.stepping)),
0130         m_logger{logger} {
0131     // Some initial printing
0132     checks(true);
0133     ACTS_VERBOSE("Gsf step "
0134                  << state.stepping.steps << " at mean position "
0135                  << stepper.position(state.stepping).transpose()
0136                  << " with direction "
0137                  << stepper.direction(state.stepping).transpose()
0138                  << " and momentum " << stepper.absoluteMomentum(state.stepping)
0139                  << " and charge " << stepper.charge(state.stepping));
0140     ACTS_VERBOSE("Propagation is in " << state.options.direction << " mode");
0141     print_component_stats();
0142   }
0143
0144   ~ScopedGsfInfoPrinterAndChecker() {
0145     if (m_navigator.currentSurface(m_state.navigation)) {
0146       const auto p_final = m_stepper.absoluteMomentum(m_state.stepping);
0147       ACTS_VERBOSE("Component status at end of step:");
0148       print_component_stats();
0149       ACTS_VERBOSE("Delta Momentum = " << std::setprecision(5)
0150                                        << p_final - m_p_initial);
0151     }
0152     checks(false);
0153   }
0154 };
0155
0156 double calculateDeterminant(
0157     const double *fullCalibratedCovariance,
0158     TrackStateTraits<kMeasurementSizeMax, true>::Covariance predictedCovariance,
0159     BoundSubspaceIndices projector, unsigned int calibratedSize);
0160
0161 /// Reweight the components according to `R. Frühwirth, "Track fitting
0162 /// with non-Gaussian noise"`. See also the implementation in Athena at
0163 /// PosteriorWeightsCalculator.cxx
0164 /// @note The weights are not renormalized!
0165 template <typename traj_t>
0166 void computePosteriorWeights(const traj_t &mt,
0167                              const std::vector<TrackIndexType> &tips,
0168                              std::map<TrackIndexType, double> &weights) {
0169   // Helper Function to compute detR
0170
0171   // Find minChi2, this can be used to factor some things later in the
0172   // exponentiation
0173   const auto minChi2 =
0174       mt.getTrackState(
0175             *std::ranges::min_element(tips,
0176                                       [&](const auto &a, const auto &b) {
0177                                         return mt.getTrackState(a).chi2() <
0178                                                mt.getTrackState(b).chi2();
0179                                       }))
0180           .chi2();
0181
0182   // Loop over the tips and compute new weights
0183   for (auto tip : tips) {
0184     const auto state = mt.getTrackState(tip);
0185     const double chi2 = state.chi2() - minChi2;
0186     const double detR = calculateDeterminant(
0187         state.effectiveCalibratedCovariance().data(),
0188         state.predictedCovariance(), state.projectorSubspaceIndices(),
0189         state.calibratedSize());
0190
0191     if (detR <= 0) {
0192       // If the determinant is not positive, just leave the weight as it is
0193       continue;
0194     }
0195
0196     const double factor = std::sqrt(1. / detR) * safeExp(-0.5 * chi2);
0197
0198     if (!std::isfinite(factor)) {
0199       // If something is not finite here, just leave the weight as it is
0200       continue;
0201     }
0202
0203     weights.at(tip) *= factor;
0204   }
0205 }
0206
0207 /// Enumeration type to allow templating on the state we want to project on with
0208 /// a MultiTrajectory
0209 enum class StatesType { ePredicted, eFiltered, eSmoothed };
0210
0211 inline std::ostream &operator<<(std::ostream &os, StatesType type) {
0212   constexpr static std::array names = {"predicted", "filtered", "smoothed"};
0213   os << names[static_cast<int>(type)];
0214   return os;
0215 }
0216
0217 /// @brief Projector type which maps a MultiTrajectory-Index to a tuple of
0218 /// [weight, parameters, covariance]. Therefore, it contains a MultiTrajectory
0219 /// and for now a std::map for the weights
0220 template <StatesType type, typename traj_t>
0221 struct MultiTrajectoryProjector {
0222   const traj_t &mt;
0223   const std::map<TrackIndexType, double> &weights;
0224
0225   auto operator()(TrackIndexType idx) const {
0226     const auto proxy = mt.getTrackState(idx);
0227     switch (type) {
0228       case StatesType::ePredicted:
0229         return std::tuple(weights.at(idx), proxy.predicted(),
0230                           proxy.predictedCovariance());
0231       case StatesType::eFiltered:
0232         return std::tuple(weights.at(idx), proxy.filtered(),
0233                           proxy.filteredCovariance());
0234       case StatesType::eSmoothed:
0235         return std::tuple(weights.at(idx), proxy.smoothed(),
0236                           proxy.smoothedCovariance());
0237       default:
0238         throw std::invalid_argument(
0239             "Incorrect StatesType, should be ePredicted"
0240             ", eFiltered, or eSmoothed.");
0241     }
0242   }
0243 };
0244
0245 /// Small Helper class that allows to carry a temporary value until we decide to
0246 /// update the actual value. The temporary value is deliberately only accessible
0247 /// with a mutable reference
0248 template <typename T>
0249 class Updatable {
0250   T m_tmp{};
0251   T m_val{};
0252
0253  public:
0254   Updatable() : m_tmp(0), m_val(0) {}
0255
0256   T &tmp() { return m_tmp; }
0257   void update() { m_val = m_tmp; }
0258
0259   const T &val() const { return m_val; }
0260 };
0261
0262 /// Remove components with low weights and renormalize from the component
0263 /// cache
0264 /// TODO This function does not expect normalized components, but this
0265 /// could be redundant work...
0266 void removeLowWeightComponents(std::vector<GsfComponent> &cmps,
0267                                double weightCutoff);
0268
0269 template <typename traj_t>
0270 struct TemporaryStates {
0271   traj_t traj;
0272   std::vector<TrackIndexType> tips;
0273   std::map<TrackIndexType, double> weights;
0274
0275   void clear() {
0276     traj.clear();
0277     tips.clear();
0278     weights.clear();
0279   }
0280 };
0281
0282 /// Function that updates the stepper from the MultiTrajectory
0283 template <typename traj_t, typename propagator_state_t, typename stepper_t>
0284 void updateStepper(propagator_state_t &state, const stepper_t &stepper,
0285                    const TemporaryStates<traj_t> &tmpStates,
0286                    double weightCutoff) {
0287   auto cmps = stepper.componentIterable(state.stepping);
0288   for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
0289     // we set ignored components to missed, so we can remove them after
0290     // the loop
0291     if (tmpStates.weights.at(idx) < weightCutoff) {
0292       cmp.status() = IntersectionStatus::unreachable;
0293       continue;
0294     }
0295
0296     auto proxy = tmpStates.traj.getTrackState(idx);
0297
0298     cmp.pars() = MultiTrajectoryHelpers::freeFiltered(state.geoContext, proxy);
0299     cmp.cov() = proxy.filteredCovariance();
0300     cmp.weight() = tmpStates.weights.at(idx);
0301   }
0302
0303   stepper.removeMissedComponents(state.stepping);
0304
0305   // TODO we have two normalization passes here now, this can probably be
0306   // optimized
0307   detail::Gsf::normalizeWeights(
0308       cmps, [&](auto cmp) -> double & { return cmp.weight(); });
0309 }
0310
0311 /// Function that updates the stepper from the ComponentCache
0312 template <typename propagator_state_t, typename stepper_t>
0313 void updateStepper(propagator_state_t &state, const stepper_t &stepper,
0314                    const Surface &surface,
0315                    const std::vector<GsfComponent> &componentCache) {
0316   // Clear components before adding new ones
0317   stepper.clearComponents(state.stepping);
0318
0319   // Finally loop over components
0320   for (const auto &[weight, pars, cov] : componentCache) {
0321     // Add the component to the stepper
0322     BoundTrackParameters bound(surface.getSharedPtr(), pars, cov,
0323                                stepper.particleHypothesis(state.stepping));
0324
0325     auto cmp = stepper.addComponent(state.stepping, std::move(bound), weight);
0326
0327     auto freeParams = cmp.pars();
0328     cmp.jacToGlobal() = surface.boundToFreeJacobian(
0329         state.geoContext, freeParams.template segment<3>(eFreePos0),
0330         freeParams.template segment<3>(eFreeDir0));
0331     cmp.pathAccumulated() = state.stepping.pathAccumulated;
0332     cmp.jacobian() = BoundMatrix::Identity();
0333     cmp.derivative() = FreeVector::Zero();
0334     cmp.jacTransport() = FreeMatrix::Identity();
0335   }
0336 }
0337
0338 double applyBetheHeitler(
0339     const GeometryContext &geoContext, const Surface &surface,
0340     Direction direction, const BoundTrackParameters &initialParameters,
0341     double initialWeight, const BetheHeitlerApprox &betheHeitlerApprox,
0342     std::vector<BetheHeitlerApprox::Component> &betheHeitlerCache,
0343     double weightCutoff, std::vector<GsfComponent> &componentCache,
0344     std::size_t &nInvalidBetheHeitler, double &maxPathXOverX0,
0345     const Logger &logger);
0346
0347 template <typename traj_t, typename propagator_state_t, typename stepper_t>
0348 void convoluteComponents(
0349     propagator_state_t &state, const stepper_t &stepper,
0350     const TemporaryStates<traj_t> &tmpStates,
0351     const BetheHeitlerApprox &betheHeitlerApprox,
0352     std::vector<BetheHeitlerApprox::Component> &betheHeitlerCache,
0353     double weightCutoff, std::vector<GsfComponent> &componentCache,
0354     std::size_t &nInvalidBetheHeitler, double &maxPathXOverX0,
0355     double &sumPathXOverX0, const Logger &logger) {
0356   const GeometryContext &geoContext = state.options.geoContext;
0357   const Direction direction = state.options.direction;
0358
0359   double pathXOverX0 = 0.0;
0360   auto cmps = stepper.componentIterable(state.stepping);
0361   for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
0362     auto proxy = tmpStates.traj.getTrackState(idx);
0363     const Surface &surface = proxy.referenceSurface();
0364
0365     BoundTrackParameters bound(surface.getSharedPtr(), proxy.filtered(),
0366                                proxy.filteredCovariance(),
0367                                stepper.particleHypothesis(state.stepping));
0368
0369     pathXOverX0 += applyBetheHeitler(
0370         geoContext, surface, direction, bound, tmpStates.weights.at(idx),
0371         betheHeitlerApprox, betheHeitlerCache, weightCutoff, componentCache,
0372         nInvalidBetheHeitler, maxPathXOverX0, logger);
0373   }
0374
0375   // Store average material seen by the components
0376   // Should not be too broadly distributed
0377   sumPathXOverX0 += pathXOverX0 / tmpStates.tips.size();
0378 }
0379
0380 /// Apply the multiple scattering to the state
0381 template <typename propagator_state_t, typename stepper_t>
0382 Result<void> applyMultipleScattering(propagator_state_t &state,
0383                                      const stepper_t &stepper,
0384                                      const Surface &surface,
0385                                      const MaterialUpdateMode &updateMode,
0386                                      const Logger &logger) {
0387   for (auto cmp : stepper.componentIterable(state.stepping)) {
0388     auto singleState = cmp.singleState(state);
0389     const auto &singleStepper = cmp.singleStepper(stepper);
0390
0391     const Result<detail::PointwiseMaterialEffects> materialInteractionRes =
0392         detail::performMaterialInteraction(
0393             singleState, singleStepper, surface, updateMode,
0394             NoiseUpdateMode::addNoise, true, false, logger);
0395     if (!materialInteractionRes.ok()) {
0396       ACTS_DEBUG("Error performing material interaction: "
0397                  << materialInteractionRes.error());
0398       return materialInteractionRes.error();
0399     }
0400
0401     assert(singleState.stepping.cov.array().isFinite().all() &&
0402            "covariance not finite after multi scattering");
0403   }
0404
0405   return Result<void>::success();
0406 }
0407
0408 }  // namespace Acts::detail::Gsf