Add maximum number of iteration steps in Propagator::AdvanceParticle (

#338)

* Add maximum number of steps in Propagator::AdvanceParticle. 
Introduce resampling of random numbers in "backscattering" case.

* If maximum number of iteration steps is reached, just propagate particle to the border with current scattering angles.

src/PROPOSAL/PROPOSAL/Constants.h

@@ -60,6 +60,11 @@ struct InterpolationSettings {
     static unsigned int NODES_RATE_INTERPOLANT;
 };
 
+// propagation settings
+struct PropagationSettings {
+    static unsigned int ADVANCE_PARTICLE_MAX_STEPS;
+};
+
 // precision parameters
 extern const double COMPUTER_PRECISION;
 extern const double HALF_PRECISION; // std::sqrt(computerPrecision);

src/PROPOSAL/detail/PROPOSAL/Constants.cxx

@@ -33,6 +33,10 @@ unsigned int InterpolationSettings::NODES_DNDX_V = 100;
 unsigned int InterpolationSettings::NODES_UTILITY = 500;
 unsigned int InterpolationSettings::NODES_RATE_INTERPOLANT = 10000;
 
+// propagation settings
+
+unsigned int PropagationSettings::ADVANCE_PARTICLE_MAX_STEPS = 200;
+
 // precision parameters
 const double PROPOSAL::COMPUTER_PRECISION = 1.e-10;
 const double PROPOSAL::HALF_PRECISION = 1.e-5; // std::sqrt(computerPrecision);

src/PROPOSAL/detail/PROPOSAL/Propagator.cxx

@@ -182,9 +182,13 @@ int Propagator::AdvanceParticle(ParticleState &state,
         return random_numbers[i++%4];
     };
 
+    int num_steps = 0; // count number of iteration steps
+    bool backscatter = false;
+
     // Iterate combinations of step lengths and scattering angles until we have
     // reached an interaction, a sector border or the maximal propagation distance
     do {
+        num_steps++;
         // Calculate grammage, energy and distance for step
         if (energy != -1 && distance == -1) {
             // Calculate grammage and distance from given energy
@@ -230,7 +234,22 @@ int Propagator::AdvanceParticle(ParticleState &state,
         // Check step
         double distance_to_border = CalculateDistanceToBorder(state.position, mean_direction, *geometry);
         bool is_inside = geometry->IsInside(state.position, mean_direction);
-        if (!is_inside) {
+
+        if (num_steps > PropagationSettings::ADVANCE_PARTICLE_MAX_STEPS) {
+            // too many iteration steps!
+            Logging::Get("proposal.propagator")->warn("Maximal number of iteration step exceeded ({}). "
+                                                      "Proposed propagation step is {} cm, while distance to border is "
+                                                      "{} cm (difference of {} cm). Initial energy particle {} MeV, "
+                                                      "particle energy after step {} MeV. Propagate to border and "
+                                                      "continue propagation.",
+                                                      PropagationSettings::ADVANCE_PARTICLE_MAX_STEPS, distance,
+                                                      distance_to_border, std::abs(distance - distance_to_border),
+                                                      state.energy, energy);
+            distance = distance_to_border;
+            grammage = density->Calculate(state.position, state.direction, distance);
+            energy = utility.EnergyDistance(state.energy, grammage);
+            advancement_type = ReachedBorder;
+        } else if (!is_inside) {
             // Special case: We are on the sector border, but scattering back outside the current sector!
             // Update sector and recalculate values
             advancement_type = InvalidStep;
@@ -242,6 +261,16 @@ int Propagator::AdvanceParticle(ParticleState &state,
             energy = energy_next_interaction;
             distance = -1;
             grammage = -1;
+
+            // if we get in this case, this might mean that we are stuck in a loop.
+            // this can happen if we backscatter in both the old and the new sector (if their medium is different).
+            // sample new random numbers to avoid this.
+            if (backscatter == true) {
+                for (auto& r: random_numbers) {
+                    r = rnd_generator();
+                }
+            }
+            backscatter = true;
         } else if (distance <= distance_to_border && distance <= max_distance && energy == energy_next_interaction) {
             // reached interaction
             advancement_type = ReachedInteraction;

src/pyPROPOSAL/detail/pybindings.cxx

@@ -237,6 +237,11 @@ PYBIND11_MODULE(proposal, m)
         .def_readwrite_static(
             "nodes_rate_interpolant", &InterpolationSettings::NODES_RATE_INTERPOLANT);
 
+    py::class_<PropagationSettings, std::shared_ptr<PropagationSettings>>(
+            m, "PropagationSettings")
+            .def_readwrite_static(
+                    "advance_particle_max_steps", &PropagationSettings::ADVANCE_PARTICLE_MAX_STEPS);
+
     /* py::class_<InterpolationDef, std::shared_ptr<InterpolationDef>>(m, */
     /*     "InterpolationDef", */
     /*     R"pbdoc( */