Acts/Propagator/RiddersPropagator.ipp

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