Implement optional logarithmic sampling of the photon energies

[ehlertdo]

Hey,
this is the follow-up to #474 . With this change it is possible to sample photons from tabulated photon fields where the maximum tabulated energy is much larger than the highest energies where the field density is non-zero.
I have attached a minimal CRPropa 1D simulation script that illustrates a scenario where the current version of CRPropa fails to find photons and stops with error: no photon found in sampleEps, please make sure that photon field provides photons for the interaction by adapting the energy range of the tabulated photon field.
The files for such a tabulated (CMB-like) photon field are also included and must be placed in the share folder of the local CRPropa installation.
Let me know what you think!
supplementary_files.tar.gz

[JulienDoerner]

Hey @ehlertdo,
changing the prior sampling will change the sampled photons.
This problem is small as long as the photon energy range is small, but for large primary energy the change is significant.

For example you can see in this plot the difference of both sampling techniques. Here I ran the sampleEps function for $10^4$ times on the CMB for a proton with $E = 10^{21} \ eV$.

I guess this problem can be solved by a diffrent way of evaluating the propability after the sampling. But I do not directly now how to implement this.

[JulienDoerner]

I guess an other way of solving the problem could be a change in the maximum photon energy

CRPropa3/src/PhotonBackground.cpp

Line 67 in 1ffe835

double TabularPhotonField::getMaximumPhotonEnergy(double z) const{

I would suggest to add a threshold photon density (maybe even 0) to evaluate the maximum energy. In this case the region with $n = 0$ would not be sampled.

[ehlertdo]

Good point about the difference in the sampled distribution. I had only checked that sampling from the build-in CMB and my tabulated version yields the same photon distribution after the modification. I will try to think about a solution.
In general, for the rejection sampling one can define an arbitrary proposal function which should ideally be similar to the true distribution of photons. One could consider providing a range of possible proposal distributions to choose from depending on the type of photon field, as opposed to just a flat prior. While this could speed up the sampling it would add a lot of complexity to the code.

[ehlertdo]

I implemented a photon sampling with a decreasing-power-law proposal function. Now the sampled distribution for interactions with a (I think) $10^{21}~\mathrm{eV}$ proton agrees with the default implementation in CRPropa. However, because each sample requires the evaluation of the inverse and regular CDF of the proposal function the sampling time increases by a factor of 5-10 (tested for 10k photons). Perhaps there is somem room for improvement here. The upside is that it works even in the scenario with extremely large maximum tabulated photon energies where the normal implementation fails.