Create `MoliereInterpol` parametrization for multiple scattering

[Jean1995]

When profiling the performance of multiple scattering calculations, we already knew that Moliere takes much longer to calculate than Highland or HighlandIntegral. For applications such as CORSIKA 8, where an accurate description of large-scattering angles is mandatory, this starts to becomes a performance issue.

When profiling the execution of Moliere, I figured out that most of the time is spent in the functions f1M, f2M, F1M F2M, where an expression is evaluated using Horner's method. Looking at the manual of EGS, I found out that they use interpolation tables to store the values of these functions instead of evaluating them for every execution.

Therefore, I implemented a similar approach for PROPOSAL, using the cubicinterpolation classes.

First of all, these are validations that the interpolation works smoothly:

Validation of f1M

 

Validation of f2M

 

Validation of F1M

 

Validation of F2M

 

The results of the scattering angles stays true to the results without the interpolation. This plots shows the scattering angles for electrons in air, with a fixed initial energy of 1 GeV and a fixed final energy of 0.9 GeV.

The runtime improvement can be seen in the following plot, which shows the computing time per scattering angle.

While the actual runtime improvement always depends on the specific usecase, it is at least of a factor of two.

This PR also include some smaller improvements, namely:

• Add a direct print/cout method for multiple_scattering::ScatteringOffset
• Add a warning if Moliere::CalculateRandomAngle doesn't seem to converge (normally, only about 5 iterations are necessary)
• Add direct python bindings for Moliere, and add a python binding for the simple constructors

[sudojan]

Very good work!

the interpolation is only up to x=20, above, the functions are again evaluated analytically:

• can you give a rough estimate on how often x>20?

• regarding the validation, the transition at x=20 seems to be smooth. But why are there larger deviations for f1m at around 11 or for f2m at around 19?

• regarding the tests, one could compare histograms between Moliere and MoliereInterpol, similar to the cross-sections

[Jean1995]

• can you give a rough estimate on how often x>20?

I analyzed this for a few different use cases:

Propagation of a 1 PeV myon in ice, 500 MeV ecut:

• x<20: 428836 times
• x>20: 1516 times
• ratio of 0.4%

Propagation of a 1 PeV myon in ice, 0.05 vcut:

• x<20: 1448times
• x>20: 0 times

Propagation of a 1 PeV myon in standardrock, 500 MeV ecut:

• x<20: 233920 times
• x>20: 868 times
• ratio of 0.4%

Propagation of a 1 PeV myon in standardrock, 0.05 vcut:

• x<20: 716 times
• x>20: 0 times

Simulation of a 10 TeV electron-induced shower, with a cut of 0.5 MeV

• x<20: 68438752 times
• x>20: 359228 times
• ratio of 0.5%

So I believe it is universally at a ratio of around 0.5%

• regarding the validation, the transition at x=20 seems to be smooth. But why are there larger deviations for f1m at around 11 or for f2m at around 19?

For f1M, the numerical evaluation method changes at x=12:

https://github.com/tudo-astroparticlephysics/PROPOSAL/blob/4bd2369b782562bd880aedd3ebceba8e8373985e/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx#L194-L197

The same happens for f2M at 18.0625:

https://github.com/tudo-astroparticlephysics/PROPOSAL/blob/4bd2369b782562bd880aedd3ebceba8e8373985e/src/PROPOSAL/detail/PROPOSAL/scattering/multiple_scattering/Moliere.cxx#L233-L235

Here, the evaluation of f2M seems to be quite instable in general before we reach 18.0625.

• regarding the tests, one could compare histograms between Moliere and MoliereInterpol, similar to the cross-sections

Good point. I already thought about this, however, the results of individual calculations are not as similar as I would have thought. But I will have a closer look at this an add appropriate tests.