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