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 <array>
0026 #include <cassert>
0027 #include <cmath>
0028 #include <cstddef>
0029 #include <iomanip>
0030 #include <map>
0031 #include <ostream>
0032 #include <tuple>
0033 #include <vector>
0034
0035 namespace Acts::detail::Gsf {
0036
0037 /// The tolerated difference to 1 to accept weights as normalized
0038 constexpr static double s_normalizationTolerance = 1.e-4;
0039
0040 template <typename component_range_t, typename projector_t>
0041 bool weightsAreNormalized(const component_range_t &cmps,
0042                           const projector_t &proj,
0043                           double tol = s_normalizationTolerance) {
0044   double sumOfWeights = 0.0;
0045
0046   for (auto it = cmps.begin(); it != cmps.end(); ++it) {
0047     sumOfWeights += proj(*it);
0048   }
0049
0050   return std::abs(sumOfWeights - 1.0) < tol;
0051 }
0052
0053 template <typename component_range_t, typename projector_t>
0054 void normalizeWeights(component_range_t &cmps, const projector_t &proj) {
0055   double sumOfWeights = 0.0;
0056
0057   // we need decltype(auto) here to support proxy-types with reference
0058   // semantics, otherwise there is a `cannot bind ... to ...` error
0059   for (auto it = cmps.begin(); it != cmps.end(); ++it) {
0060     decltype(auto) cmp = *it;
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 it = cmps.begin(); it != cmps.end(); ++it) {
0068     decltype(auto) cmp = *it;
0069     proj(cmp) /= sumOfWeights;
0070   }
0071 }
0072
0073 // A class that prints information about the state on construction and
0074 // destruction, it also contains some assertions in the constructor and
0075 // destructor. It can be removed without change of behaviour, since it only
0076 // holds const references
0077 template <typename propagator_state_t, typename stepper_t, typename navigator_t>
0078 class ScopedGsfInfoPrinterAndChecker {
0079   const propagator_state_t &m_state;
0080   const stepper_t &m_stepper;
0081   const navigator_t &m_navigator;
0082   double m_p_initial;
0083   const Logger &m_logger;
0084
0085   const Logger &logger() const { return m_logger; }
0086
0087   void print_component_stats() const {
0088     std::size_t i = 0;
0089     for (auto cmp : m_stepper.constComponentIterable(m_state.stepping)) {
0090       auto getVector = [&](auto idx) {
0091         return cmp.pars().template segment<3>(idx).transpose();
0092       };
0093       ACTS_VERBOSE("  #" << i++ << " pos: " << getVector(eFreePos0) << ", dir: "
0094                          << getVector(eFreeDir0) << ", weight: " << cmp.weight()
0095                          << ", status: " << cmp.status()
0096                          << ", qop: " << cmp.pars()[eFreeQOverP]
0097                          << ", det(cov): " << cmp.cov().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.stepping.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(*std::min_element(tips.begin(), tips.end(),
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                    const Logger &logger) {
0317   // Clear components before adding new ones
0318   stepper.clearComponents(state.stepping);
0319
0320   // Finally loop over components
0321   for (const auto &[weight, pars, cov] : componentCache) {
0322     // Add the component to the stepper
0323     BoundTrackParameters bound(surface.getSharedPtr(), pars, cov,
0324                                stepper.particleHypothesis(state.stepping));
0325
0326     auto res = stepper.addComponent(state.stepping, std::move(bound), weight);
0327
0328     if (!res.ok()) {
0329       ACTS_ERROR("Error adding component to MultiStepper");
0330       continue;
0331     }
0332
0333     auto &cmp = *res;
0334     auto freeParams = cmp.pars();
0335     cmp.jacToGlobal() = surface.boundToFreeJacobian(
0336         state.geoContext, freeParams.template segment<3>(eFreePos0),
0337         freeParams.template segment<3>(eFreeDir0));
0338     cmp.pathAccumulated() = state.stepping.pathAccumulated;
0339     cmp.jacobian() = BoundMatrix::Identity();
0340     cmp.derivative() = FreeVector::Zero();
0341     cmp.jacTransport() = FreeMatrix::Identity();
0342   }
0343 }
0344
0345 double applyBetheHeitler(
0346     const GeometryContext &geoContext, const Surface &surface,
0347     Direction direction, const BoundTrackParameters &initialParameters,
0348     double initialWeight, const BetheHeitlerApprox &betheHeitlerApprox,
0349     std::vector<BetheHeitlerApprox::Component> &betheHeitlerCache,
0350     double weightCutoff, std::vector<GsfComponent> &componentCache,
0351     std::size_t &nInvalidBetheHeitler, double &maxPathXOverX0,
0352     const Logger &logger);
0353
0354 template <typename traj_t, typename propagator_state_t, typename stepper_t>
0355 void convoluteComponents(
0356     propagator_state_t &state, const stepper_t &stepper,
0357     const TemporaryStates<traj_t> &tmpStates,
0358     const BetheHeitlerApprox &betheHeitlerApprox,
0359     std::vector<BetheHeitlerApprox::Component> &betheHeitlerCache,
0360     double weightCutoff, std::vector<GsfComponent> &componentCache,
0361     std::size_t &nInvalidBetheHeitler, double &maxPathXOverX0,
0362     double &sumPathXOverX0, const Logger &logger) {
0363   const GeometryContext &geoContext = state.options.geoContext;
0364   const Direction direction = state.options.direction;
0365
0366   double pathXOverX0 = 0.0;
0367   auto cmps = stepper.componentIterable(state.stepping);
0368   for (auto [idx, cmp] : zip(tmpStates.tips, cmps)) {
0369     auto proxy = tmpStates.traj.getTrackState(idx);
0370     const Surface &surface = proxy.referenceSurface();
0371
0372     BoundTrackParameters bound(surface.getSharedPtr(), proxy.filtered(),
0373                                proxy.filteredCovariance(),
0374                                stepper.particleHypothesis(state.stepping));
0375
0376     pathXOverX0 += applyBetheHeitler(
0377         geoContext, surface, direction, bound, tmpStates.weights.at(idx),
0378         betheHeitlerApprox, betheHeitlerCache, weightCutoff, componentCache,
0379         nInvalidBetheHeitler, maxPathXOverX0, logger);
0380   }
0381
0382   // Store average material seen by the components
0383   // Should not be too broadly distributed
0384   sumPathXOverX0 += pathXOverX0 / tmpStates.tips.size();
0385 }
0386
0387 /// Apply the multiple scattering to the state
0388 template <typename propagator_state_t, typename stepper_t>
0389 void applyMultipleScattering(propagator_state_t &state,
0390                              const stepper_t &stepper, const Surface &surface,
0391                              const MaterialUpdateMode &updateMode,
0392                              const Logger &logger) {
0393   for (auto cmp : stepper.componentIterable(state.stepping)) {
0394     auto singleState = cmp.singleState(state);
0395     const auto &singleStepper = cmp.singleStepper(stepper);
0396
0397     detail::performMaterialInteraction(singleState, singleStepper, surface,
0398                                        updateMode, NoiseUpdateMode::addNoise,
0399                                        true, false, logger);
0400
0401     assert(singleState.stepping.cov.array().isFinite().all() &&
0402            "covariance not finite after multi scattering");
0403   }
0404 }
0405
0406 }  // namespace Acts::detail::Gsf