Propagator/detail/LoopProtection.hpp

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