Acts/Propagator/RiddersPropagator.ipp

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/Propagator/RiddersPropagator.hpp"
0012
0013 #include "Acts/Definitions/TrackParametrization.hpp"
0014
0015 #include <numbers>
0016
0017 namespace Acts {
0018
0019 namespace {
0020
0021 /// @brief This function tests whether the variations on a disc as target
0022 /// surface lead to results on different sides wrt the center of the disc.
0023 /// This would lead to a flip of the phi value on the surface and therewith
0024 /// to a huge variance in that parameter. It can only occur in this
0025 /// algorithm since the ridders algorithm is unaware of the target surface.
0026 ///
0027 /// @param [in] derivatives Derivatives of a single parameter
0028 ///
0029 /// @return Boolean result whether a phi jump occurred
0030 inline bool inconsistentDerivativesOnDisc(
0031     const std::vector<BoundVector>& derivatives) {
0032   // Test each component with each other
0033   for (unsigned int i = 0; i < derivatives.size(); i++) {
0034     bool jumpedAngle = true;
0035     for (unsigned int j = 0; j < derivatives.size(); j++) {
0036       // If there is at least one with a similar angle then it seems to work
0037       // properly
0038       if (i != j && std::abs(derivatives[i](1) - derivatives[j](1)) <
0039                         std::numbers::pi / 2.) {
0040         jumpedAngle = false;
0041         break;
0042       }
0043     }
0044     // Break if a jump was detected
0045     if (jumpedAngle) {
0046       return true;
0047     }
0048   }
0049   return false;
0050 }
0051
0052 /// @brief This function fits a linear function through the final state
0053 /// parametrisations
0054 ///
0055 /// @param [in] deviations Vector of deviations
0056 /// @param [in] derivatives Vector of resulting derivatives
0057 ///
0058 /// @return Vector containing the linear fit
0059 inline BoundVector fitLinear(const std::vector<double>& deviations,
0060                              const std::vector<BoundVector>& derivatives) {
0061   assert(!derivatives.empty() && "no fit without data");
0062   if (derivatives.size() == 1) {
0063     return derivatives[0];
0064   }
0065
0066   BoundVector A = BoundVector::Zero();
0067   BoundVector C = BoundVector::Zero();
0068   double B = 0;
0069   double D = 0;
0070   const unsigned int N = deviations.size();
0071
0072   for (unsigned int i = 0; i < N; ++i) {
0073     A += deviations.at(i) * derivatives.at(i);
0074     B += deviations.at(i);
0075     C += derivatives.at(i);
0076     D += deviations.at(i) * deviations.at(i);
0077   }
0078
0079   BoundVector b = (N * A - B * C) / (N * D - B * B);
0080   BoundVector a = (C - B * b) / N;
0081
0082   return a;
0083 }
0084
0085 }  // namespace
0086
0087 template <typename propagator_t>
0088 template <typename parameters_t, typename propagator_options_t>
0089 auto RiddersPropagator<propagator_t>::propagate(
0090     const parameters_t& start, const propagator_options_t& options) const
0091     -> Result<actor_list_t_result_t<
0092         BoundTrackParameters, typename propagator_options_t::actor_list_type>> {
0093   using ThisResult = Result<actor_list_t_result_t<
0094       BoundTrackParameters, typename propagator_options_t::actor_list_type>>;
0095
0096   // Remove the covariance from our start parameters in order to skip jacobian
0097   // transport for the nominal propagation
0098   BoundTrackParameters startWithoutCov = start;
0099   startWithoutCov.covariance() = std::nullopt;
0100
0101   // Propagate the nominal parameters
0102   auto result = m_propagator.propagate(startWithoutCov, options);
0103   if (!result.ok()) {
0104     return ThisResult::failure(result.error());
0105   }
0106   // Extract results from the nominal propagation
0107   auto nominalResult = result.value();
0108   assert(nominalResult.endParameters);
0109   const auto& nominalFinalParameters = *nominalResult.endParameters;
0110
0111   BoundMatrix jacobian =
0112       wiggleAndCalculateJacobian(startWithoutCov, options, nominalResult);
0113   nominalResult.transportJacobian = jacobian;
0114
0115   if (start.covariance()) {
0116     // use nominal parameters and Ridders covariance
0117     BoundMatrix cov = jacobian * (*start.covariance()) * jacobian.transpose();
0118     // replace the covariance of the nominal result w/ the ridders covariance
0119     nominalResult.endParameters = BoundTrackParameters::createCurvilinear(
0120         nominalFinalParameters.fourPosition(options.geoContext),
0121         nominalFinalParameters.direction(), nominalFinalParameters.qOverP(),
0122         std::move(cov), nominalFinalParameters.particleHypothesis());
0123   }
0124
0125   return ThisResult::success(std::move(nominalResult));
0126 }
0127
0128 template <typename propagator_t>
0129 template <typename parameters_t, typename propagator_options_t>
0130 auto RiddersPropagator<propagator_t>::propagate(
0131     const parameters_t& start, const Surface& target,
0132     const propagator_options_t& options) const
0133     -> Result<actor_list_t_result_t<
0134         BoundTrackParameters, typename propagator_options_t::actor_list_type>> {
0135   using ThisResult = Result<actor_list_t_result_t<
0136       BoundTrackParameters, typename propagator_options_t::actor_list_type>>;
0137
0138   // Remove the covariance from our start parameters in order to skip jacobian
0139   // transport for the nominal propagation
0140   BoundTrackParameters startWithoutCov = start;
0141   startWithoutCov.covariance() = std::nullopt;
0142
0143   // Propagate the nominal parameters
0144   auto result = m_propagator.propagate(startWithoutCov, target, options);
0145   if (!result.ok()) {
0146     return ThisResult::failure(result.error());
0147   }
0148   // Extract results from the nominal propagation
0149   auto nominalResult = result.value();
0150   assert(nominalResult.endParameters);
0151   const auto& nominalFinalParameters = *nominalResult.endParameters;
0152
0153   BoundMatrix jacobian =
0154       wiggleAndCalculateJacobian(startWithoutCov, options, nominalResult);
0155   nominalResult.transportJacobian = jacobian;
0156
0157   if (start.covariance()) {
0158     // use nominal parameters and Ridders covariance
0159     BoundMatrix cov = jacobian * (*start.covariance()) * jacobian.transpose();
0160     // replace the covariance of the nominal result w/ the ridders covariance
0161     nominalResult.endParameters = BoundTrackParameters(
0162         nominalFinalParameters.referenceSurface().getSharedPtr(),
0163         nominalFinalParameters.parameters(), std::move(cov),
0164         nominalFinalParameters.particleHypothesis());
0165   }
0166
0167   return ThisResult::success(std::move(nominalResult));
0168 }
0169
0170 template <typename propagator_t>
0171 template <typename parameters_t, typename propagator_options_t>
0172 BoundMatrix RiddersPropagator<propagator_t>::wiggleAndCalculateJacobian(
0173     const parameters_t& start, const propagator_options_t& options,
0174     const actor_list_t_result_t<BoundTrackParameters,
0175                                 typename propagator_options_t::actor_list_type>&
0176         nominalResult) const {
0177   const auto& nominalFinalParameters = nominalResult.endParameters.value();
0178   // Use the curvilinear surface of the propagated parameters as target
0179   const Surface& target = nominalFinalParameters.referenceSurface();
0180
0181   // Allow larger distances to reach the target surface
0182   propagator_options_t opts = options;
0183   opts.pathLimit *= 2;
0184
0185   // Derivations of each parameter around the nominal parameters
0186   BoundMatrix jacobian = BoundMatrix::Zero();
0187
0188   // Wiggle each dimension individually
0189   for (unsigned int i = 0; i < eBoundSize; i++) {
0190     // Steps for estimating derivatives
0191     std::vector<double> deviations = options.deviationFactors;
0192     for (double& d : deviations) {
0193       d *= options.deviationScale[i];
0194     }
0195
0196     std::vector<BoundVector> derivatives =
0197         wiggleParameter(opts, start, i, target,
0198                         nominalFinalParameters.parameters(), deviations);
0199
0200     // Test if target is disc - this may lead to inconsistent results
0201     if (target.type() == Surface::Disc) {
0202       assert(!inconsistentDerivativesOnDisc(derivatives));
0203     }
0204
0205     jacobian.col(i) = fitLinear(deviations, derivatives);
0206   }
0207
0208   return jacobian;
0209 }
0210
0211 template <typename propagator_t>
0212 template <typename propagator_options_t, typename parameters_t>
0213 std::vector<BoundVector> RiddersPropagator<propagator_t>::wiggleParameter(
0214     const propagator_options_t& options, const parameters_t& start,
0215     const unsigned int param, const Surface& target, const BoundVector& nominal,
0216     const std::vector<double>& deviations) const {
0217   // Storage of the results
0218   std::vector<BoundVector> derivatives;
0219   derivatives.reserve(deviations.size());
0220   for (double h : deviations) {
0221     // Treatment for theta
0222     if (param == eBoundTheta) {
0223       const double current_theta = start.template get<eBoundTheta>();
0224       if (current_theta + h > std::numbers::pi) {
0225         h = std::numbers::pi - current_theta;
0226       }
0227       if (current_theta + h < 0) {
0228         h = -current_theta;
0229       }
0230     }
0231
0232     // Modify start parameter
0233     BoundVector values = start.parameters();
0234     values[param] += h;
0235
0236     // Propagate with updated start parameters
0237     BoundTrackParameters tp(start.referenceSurface().getSharedPtr(), values,
0238                             start.covariance(), start.particleHypothesis());
0239     const auto& r = m_propagator.propagate(tp, target, options).value();
0240     // Collect the slope
0241     derivatives.push_back((r.endParameters->parameters() - nominal) / h);
0242
0243     // Correct for a possible variation of phi around
0244     if (param == eBoundPhi) {
0245       double phi0 = nominal(eBoundPhi);
0246       double phi1 = r.endParameters->parameters()(eBoundPhi);
0247       if (std::abs(phi1 + 2. * std::numbers::pi - phi0) <
0248           std::abs(phi1 - phi0)) {
0249         derivatives.back()[eBoundPhi] =
0250             (phi1 + 2. * std::numbers::pi - phi0) / h;
0251       } else if (std::abs(phi1 - 2. * std::numbers::pi - phi0) <
0252                  std::abs(phi1 - phi0)) {
0253         derivatives.back()[eBoundPhi] =
0254             (phi1 - 2. * std::numbers::pi - phi0) / h;
0255       }
0256     }
0257   }
0258
0259   return derivatives;
0260 }
0261
0262 }  // namespace Acts