Propagation can stall at border transitions

[Jean1995]

I propagate this muon

    auto debug_state = ParticleState();
    debug_state.energy = 1e9;
    debug_state.position = Cartesian3D(2.19878e+06, 0, 387704);
    debug_state.direction = Cartesian3D(-0.984808, -9.84808e-07, -0.173648);

with the following config file:

click to show config file

{
  "global":
  {
  	"do_interpolation" : true,
  	"exact_time" : true,
  	"scattering": 
  	{
  		"multiple_scattering" : "HighlandIntegral",
  		"stochastic_deflection" : ["DefaultIoniz", "DefaultBrems", "DefaultEpair", "DefaultPhotonuclear"]
  	},
  	"cuts":
  	{
  		"e_cut": "inf",
  		"v_cut": 0.05,
  		"cont_rand": true
  	},
  	"CrossSections" : {
  		"brems": {
  			"parametrization": "KelnerKokoulinPetrukhin",
  			"lpm": true
  		},
  		"epair": {
  			"parametrization": "KelnerKokoulinPetrukhin",
  			"lpm": true
  		},
  		"ioniz": {
  			"parametrization": "BetheBlochRossi"
  		},
  		"photo": {
  			"parametrization": "AbramowiczLevinLevyMaor97"
  		}
  	}
  },
  "sectors": [
  	{
  		"medium": "air",
  		"geometries": [
  			{
  				"hierarchy": 1,
  				"shape": "sphere",
  				"origin": [0, 0, -637218600],
  				"outer_radius": 1e20,
  				"inner_radius": 637413400

  			}
  		],
  		"density_distribution":
  		{
  			"type": "homogeneous",
  			"mass_density" : 0.810965e-3
  		}
  	},
  	{
  		"medium": "ice",
  		"geometries": [
  			{
  				"hierarchy": 1,
  				"shape": "sphere",
  				"origin": [0, 0, -637218600],
  				"outer_radius": 637413400,
  				"inner_radius": 637393400

  			}
  		],
  		"density_distribution":
  		{
  			"type": "homogeneous",
  			"mass_density" : 0.763776
  		}
  	},
  	{
  		"medium": "ice",
  		"geometries": [
  			{
  				"hierarchy": 1,
  				"shape": "sphere",
  				"origin": [0, 0, -637218600],
  				"outer_radius": 637393400,
  				"inner_radius": 637132400
  			}
  		],
  		"density_distribution":
  		{
  			"type": "homogeneous",
  			"mass_density" : 0.92259
  		}
  	},
  	{
  		"medium": "ice",
  		"geometries": [
  			{
  				"hierarchy": 100,
  				"shape": "cylinder",
  				"origin": [0, 0, 0],
  				"outer_radius": 80000,
  				"inner_radius": 0,
  				"height": 160000
  			}
  		],
  		"density_distribution":
  		{
  			"type": "homogeneous",
  			"mass_density" : 0.92259
  		},
  		"cuts":
  		{
  			"e_cut": 500,
  			"v_cut": 1,
  			"cont_rand": false
  		}
  	},
  	{
  		"medium": "standardrock",
  		"geometries": [
  			{
  				"hierarchy": 1,
  				"shape": "sphere",
  				"origin": [0, 0, -637218600],
  				"outer_radius": 637132400,
  				"inner_radius": 0
  			}
  		],
  		"density_distribution":
  		{
  			"type": "homogeneous",
  			"mass_density" : 2.65
  		}
  	}
  ]
}

This propagation immediately stalls and does not finish.

The reason for this appears to be the AdvanceParticle method, where the algorithm tries to find a combination of energy and distance in such a way that the particle hits a sector border. If the difference between the proposed step distance and the actual distance to the sector border is smaller than PARTICLE_POSITION_RESOLUTION (1e-3 cm), we accept this step.
In this case however, the algorithm gets stuck with a step where the difference is ~0.0044cm.

This is related to issue #267.

As a quick solution, one might want to limit the maximum number of iterations for the AdvanceParticle method. However, it is not entirely clear what do if we hit this iteration limit. Discard the random numbers? Raise an error? Just accept the current step (which, in the case of this example, would be the preferred solution).

[sudojan]

I do not completely get the problem. When the difference between the proposed distance and the distance to the sector border is smaller than the particle resolution, then the distance to propagate is set to the distance to the border, as implemented here:
https://github.com/tudo-astroparticlephysics/PROPOSAL/blob/master/src/PROPOSAL/detail/PROPOSAL/Propagator.cxx#L251
This looks fine to me. Why does the propagation get stuck, here?

[sudojan]

Nevertheless, I think it is a good idea to introduce a limit on the number of iterations, which the user can adapt.
However, the default should be reasonably high, maybe 10.000.
When reaching this limit, I propose to give a warning and terminate the propagation without raising an Error.

[Jean1995]

I do not completely get the problem. When the difference between the proposed distance and the distance to the sector border is smaller than the particle resolution, then the distance to propagate is set to the distance to the border, as implemented here: https://github.com/tudo-astroparticlephysics/PROPOSAL/blob/master/src/PROPOSAL/detail/PROPOSAL/Propagator.cxx#L251 This looks fine to me. Why does the propagation get stuck, here?

The values for the proposed distance and the distance to the sector border are 1116350.78826547 and 1116350.78387037, and the difference is smaller than the particle resolution. And with every iteration, these two values just swap. This means that the iteration is never finished.

[Jean1995]

I tried to debug the problem a bit more:

For the example here, the particle makes a propagation step of over 10 km. For this step, we calculate the intersection with a huge sphere. In the algorithm, this intersection is calculated for two very similar direction vectors (Cartesian3D(-0.984807997992707, -9.78397641566135e-07,-0.173648011383973) and Cartesian3D(-0.984807997870945, -9.78008788121459e-07, -0.173648012074526)). They are only different by 2e-6°.

I believe that this calculation is just not precise enough (numerical cancellation due to multiplication with large numbers maybe?). If I calculate the interaction points for these two different directions, I receive

x: 1099388.8151507 y: -1.09223497839943 z: 193851.905610771

x: 1099388.81961496 y: -1.0918008772515 z: 193851.905603071

and if I numerically calculate the distance of these two interactions points to our sphere, the results are numerically identical.

Therefore I guess that the propagation stalls because calculating the intersection points over a 10km step with two directions that vary only by ~1e-6° just numerically can't provide the accuracy that we require in the algorithm (1e-3 cm).

As soon as a user choses to propagate particle with large propagation steps, high energies, multiple scattering enabled and border transitions, this problem is very likely to come up. Therefore, something better than a quick fix might be necessary.

Maybe we should introduce something like a "relative accuracy", in addition to the current "absolute accuracy" (1e-3 cm) in the algorithm. So something like a combination using max(1e-3, relative_accuracy).

[sudojan]

a relative accuracy parameter, next to the absolute PARTICLE_POSITION_RESOLUTION is a good enhancement.
Furthermore, it would be helpful, if these parameters can be adjusted by the user.