fix: Fix Fatras end of world (#2580)

Sometimes the Fatras propagation will not terminate and step with `double max` distance which results in FPEs

After these changes the FPEs during propagation disappeared

Core/include/Acts/Propagator/Navigator.hpp

@@ -274,7 +274,7 @@ class Navigator {
 
   bool targetReached(const State& state) const { return state.targetReached; }
 
-  bool endOfWorldReached(State& state) const {
+  bool endOfWorldReached(const State& state) const {
     return state.currentVolume == nullptr;
   }
 

Examples/Algorithms/Fatras/include/ActsExamples/Fatras/FatrasSimulation.hpp

@@ -78,6 +78,11 @@ class FatrasSimulation final : public IAlgorithm {
     /// have associated material.
     bool generateHitsOnPassive = false;
 
+    /// Absolute maximum step size
+    double maxStepSize = 1 * Acts::UnitConstants::m;
+    /// Absolute maximum path length
+    double pathLimit = 30 * Acts::UnitConstants::m;
+
     /// Expected average number of hits generated per particle.
     ///
     /// This is just a performance optimization hint and has no impact on the

Examples/Algorithms/Fatras/src/FatrasSimulation.cpp

@@ -157,6 +157,11 @@ struct FatrasSimulationT final : ActsExamples::detail::FatrasSimulation {
     simulation.charged.selectHitSurface.sensitive = cfg.generateHitsOnSensitive;
     simulation.charged.selectHitSurface.material = cfg.generateHitsOnMaterial;
     simulation.charged.selectHitSurface.passive = cfg.generateHitsOnPassive;
+
+    simulation.charged.maxStepSize = cfg.maxStepSize;
+    simulation.charged.pathLimit = cfg.pathLimit;
+    simulation.neutral.maxStepSize = cfg.maxStepSize;
+    simulation.neutral.pathLimit = cfg.pathLimit;
   }
   ~FatrasSimulationT() final = default;
 

Examples/Python/tests/root_file_hashes.txt

@@ -1,5 +1,5 @@
 test_pythia8__pythia8_particles.root: d9a3d4f16c01aab9984e0e848292760d44015f210878866cc94f3bb9d030034b
-test_fatras__particles_final.root: d3a900a7151623109f2e22e11beb477635eb549d5e611c5fb330a59b41054f59
+test_fatras__particles_final.root: a58450a6145f682ddaf39b6be53088074b35847e482300908c305e693db49145
 test_fatras__particles_initial.root: 0e2d9974fdd4aa5549c5c05817ec54fe6c7d9cbe7969647f227956167e31110d
 test_fatras__hits.root: d3ec06cd12e3737b5c783302ca29a532c885dc35930901d86b2f6d17041037c6
 test_geant4__particles_final.root: abf91622e96c0723651e864e050c7fdeee84b230aacb1c261ad81ff4e5df5f2c
@@ -8,17 +8,17 @@ test_geant4__hits.root: 70b592a546fd362c9341d87c9068b1f8fed657e0036e62e7262c01fa
 test_seeding__estimatedparams.root: f11b6bdc2b9b13e84f974f2dfbd9f409255f4f0eab212fe8581c352b3f2cfc27
 test_seeding__performance_seeding.root: 992f9c611d30dde0d3f3ab676bab19ada61ab6a4442828e27b65ec5e5b7a2880
 test_seeding__particles.root: 74e08ee12bdaf9f7d273369c71c742df345475c8b73bb7ce223619d0dd1f87ee
-test_seeding__particles_final.root: 3f3eeaedeac050e3dfc5fef37734b0876f4470810e93b371e5da2cb0b987bebd
+test_seeding__particles_final.root: 8737a6c9c16340839bab45ddc8c92cace98d2c30f923a71e2bd992d777f25e7c
 test_seeding__particles_initial.root: 74e08ee12bdaf9f7d273369c71c742df345475c8b73bb7ce223619d0dd1f87ee
 test_seeding_orthogonal__estimatedparams.root: 68cece5ac084cfe12a7fd7ad602ce736edb74ec8f45e0d0468923cd90ec4f95b
 test_seeding_orthogonal__performance_seeding.root: 60fbedcf5cb2b37cd8e526251940564432890d3a159d231ed819e915a904682c
 test_seeding_orthogonal__particles.root: 74e08ee12bdaf9f7d273369c71c742df345475c8b73bb7ce223619d0dd1f87ee
-test_seeding_orthogonal__particles_final.root: 3f3eeaedeac050e3dfc5fef37734b0876f4470810e93b371e5da2cb0b987bebd
+test_seeding_orthogonal__particles_final.root: 8737a6c9c16340839bab45ddc8c92cace98d2c30f923a71e2bd992d777f25e7c
 test_seeding_orthogonal__particles_initial.root: 74e08ee12bdaf9f7d273369c71c742df345475c8b73bb7ce223619d0dd1f87ee
 test_itk_seeding__estimatedparams.root: c3ac3c9894b487def066c262ce5aeae3508a0b5a68533f23b5996f2f86f3919b
 test_itk_seeding__performance_seeding.root: 78ebda54cd0f026ba4b7f316724ffd946de56a932735914baf1b7bba9505c29d
 test_itk_seeding__particles.root: 88315e93ed4cb5d40a8721502048a9d1fc100e0a7d504e25fd4502c8302f1578
-test_itk_seeding__particles_final.root: 0d38f1a4ede5fb86844c156c002ada340011e140a13dacffc1b3848f8bf34dfa
+test_itk_seeding__particles_final.root: 3ee17fd8461d2360ef3f813bc697aef3bdb8079c069b3eddd8dfb78d5015d0e6
 test_itk_seeding__particles_initial.root: 88315e93ed4cb5d40a8721502048a9d1fc100e0a7d504e25fd4502c8302f1578
 test_propagation__propagation_steps.root: 174301b25784dbb881196b658f2d7f99c0a2ea688a0129e6110fc19aa5cf8e54
 test_material_recording__geant4_material_tracks.root: e411152d370775463c22b19a351dfc7bfe40b51985e10a7c1a010aebde80715d

Fatras/include/ActsFatras/Kernel/Simulation.hpp

@@ -44,6 +44,10 @@ template <typename propagator_t, typename interactions_t,
 struct SingleParticleSimulation {
   /// How and within which geometry to propagate the particle.
   propagator_t propagator;
+  /// Absolute maximum step size
+  double maxStepSize = std::numeric_limits<double>::max();
+  /// Absolute maximum path length
+  double pathLimit = std::numeric_limits<double>::max();
   /// Decay module.
   decay_t decay;
   /// Interaction list containing the simulated interactions.
@@ -83,6 +87,8 @@ struct SingleParticleSimulation {
 
     // Construct per-call options.
     PropagatorOptions options(geoCtx, magCtx);
+    options.maxStepSize = maxStepSize;
+    options.pathLimit = pathLimit;
     // setup the interactor as part of the propagator options
     auto &actor = options.actionList.template get<Actor>();
     actor.generator = &generator;

Fatras/include/ActsFatras/Kernel/detail/SimulationActor.hpp

@@ -10,13 +10,15 @@
 
 #include "Acts/Material/ISurfaceMaterial.hpp"
 #include "Acts/Propagator/ConstrainedStep.hpp"
+#include "Acts/Propagator/StandardAborters.hpp"
 #include "Acts/Surfaces/Surface.hpp"
 #include "ActsFatras/EventData/Hit.hpp"
 #include "ActsFatras/EventData/Particle.hpp"
 #include "ActsFatras/Kernel/SimulationResult.hpp"
 
 #include <algorithm>
 #include <cassert>
+#include <cmath>
 #include <limits>
 
 namespace ActsFatras {
@@ -84,14 +86,19 @@ struct SimulationActor {
             typename navigator_t>
   void operator()(propagator_state_t &state, stepper_t &stepper,
                   navigator_t &navigator, result_type &result,
-                  const Acts::Logger & /*logger*/) const {
+                  const Acts::Logger &logger) const {
     assert(generator and "The generator pointer must be valid");
 
     // actors are called once more after the propagation terminated
     if (not result.isAlive) {
       return;
     }
 
+    if (Acts::EndOfWorldReached{}(state, stepper, navigator, logger)) {
+      result.isAlive = false;
+      return;
+    }
+
     // check if we are still on the start surface and skip if so
     if ((navigator.startSurface(state.navigation) != nullptr) &&
         (navigator.startSurface(state.navigation) ==
@@ -116,8 +123,9 @@ struct SimulationActor {
     }
 
     // decay check. needs to happen at every step, not just on surfaces.
-    if (result.properTimeLimit - result.particle.properTime() <
-        result.properTimeLimit * properTimeRelativeTolerance) {
+    if (std::isfinite(result.properTimeLimit) &&
+        (result.properTimeLimit - result.particle.properTime() <
+         result.properTimeLimit * properTimeRelativeTolerance)) {
       auto descendants = decay.run(generator, result.particle);
       for (auto &&descendant : descendants) {
         result.generatedParticles.emplace_back(std::move(descendant));
@@ -128,6 +136,7 @@ struct SimulationActor {
 
     // Regulate the step size
     if (std::isfinite(result.properTimeLimit)) {
+      assert(result.particle.mass() > 0.0 && "Particle must have mass");
       //    beta² = p²/E²
       //    gamma = 1 / sqrt(1 - beta²) = sqrt(m² + p²) / m = E / m
       //     time = proper-time * gamma

Tests/UnitTests/Fatras/Kernel/SimulationActorTests.cpp

@@ -134,13 +134,21 @@ struct MockNavigator {
     return state.targetReached;
   }
 
+  void targetReached(MockNavigatorState &state, bool reached) const {
+    state.targetReached = reached;
+  }
+
   const Acts::Surface *startSurface(const MockNavigatorState &state) const {
     return state.startSurface;
   }
 
   const Acts::Surface *currentSurface(const MockNavigatorState &state) const {
     return state.currentSurface;
   }
+
+  bool endOfWorldReached(const MockNavigatorState & /*state*/) const {
+    return false;
+  }
 };
 
 struct MockPropagatorState {