Propagator/detail/LoopProtection.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/Algebra.hpp"
0012 #include "Acts/Utilities/Logger.hpp"
0013
0014 #include <numbers>
0015
0016 namespace Acts::detail {
0017
0018 /// Estimate the loop protection limit
0019 template <typename path_aborter_t, typename propagator_state_t,
0020           typename stepper_t>
0021 void setupLoopProtection(propagator_state_t& state, const stepper_t& stepper,
0022                          path_aborter_t& pathAborter, bool releaseLimit,
0023                          const Logger& logger) {
0024   if (!state.options.loopProtection) {
0025     return;
0026   }
0027
0028   if (releaseLimit) {
0029     pathAborter.internalLimit =
0030         state.options.direction * std::numeric_limits<double>::max();
0031   }
0032
0033   // Get the field at the start position
0034   auto fieldRes =
0035       stepper.getField(state.stepping, stepper.position(state.stepping));
0036   if (!fieldRes.ok()) {
0037     // there's no great way to return the error here, so resort to warning
0038     // and not applying the loop protection in this case
0039     ACTS_WARNING("Field lookup was unsuccessful, this is very likely an error");
0040     return;
0041   }
0042   const Vector3 field = *fieldRes;
0043   const double B = field.norm();
0044   if (B == 0) {
0045     return;
0046   }
0047
0048   // Transverse component at start is taken for the loop protection
0049   const double p = stepper.absoluteMomentum(state.stepping);
0050   // Calculate the full helix path
0051   const double helixPath =
0052       state.options.direction * 2 * std::numbers::pi * p / B;
0053   // And set it as the loop limit if it overwrites the internal limit
0054   const double loopLimit = state.options.loopFraction * helixPath;
0055   const double previousLimit = pathAborter.internalLimit;
0056   if (std::abs(loopLimit) < std::abs(previousLimit)) {
0057     pathAborter.internalLimit = loopLimit;
0058
0059     ACTS_VERBOSE("Path aborter limit set to "
0060                  << loopLimit << " (full helix = " << helixPath
0061                  << ", previous limit = " << previousLimit << ")");
0062   } else {
0063     ACTS_VERBOSE("Path aborter limit not updated to "
0064                  << loopLimit << " (full helix = " << helixPath
0065                  << ", previous limit = " << previousLimit << ")");
0066   }
0067 }
0068
0069 }  // namespace Acts::detail