Improved Photon Field Sampling

[ehlertdo]

Hi CRPropa team,

Brief Description

Background

I was recently trying to use CRPropa with custom, tabulated photon fields relevant for AGN environments. I have processed the fields with CRPropa3-data to get the right format to use with CRPropa and, excluding some minor issues with converting the units correctly, this works generally fine. To check that everything works correctly I also created a tabulated CMB (simple black body) to compare with crp.CMB(). For programmatic reasons my maximum tabulated photon energy (4e10 eV) is very large compared to the typical CMB photon temperature (CRPropa current, epsMax_CMB=0.1 eV).

Problem

I noticed that in this case there is a problem with the sampling of photon energies when running CRPropa (e.g. for photopion) with fields where the maximum tabulated energy is much larger than any energies where the photon density is non-negligible.
This gives the following error 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.. By truncating the tabulated CMB files at lower energies, e.g. 0.1 eV as is the current default, the error can be avoided. However, ideally CRPropa would work correctly no matter the tabulated range of energies.

Solution

After some digging I noticed that CRPropa appears to use a linear sampling for the photon field energies when processing the interactions. At present, only the maximum interaction probability can be sampled logarithmically (see #368). The relevant line of code is in PhotopionProduction::sampleEps(...) where

double eps = epsMin + random.rand() * (epsMax - epsMin);

I have verified that by replacing this with a logarithmic sampling, using the inverse CDF of the reciprocal distribution,

double eps = epsMin * pow(epsMax /epsMin, random.rand());

the sampling works correctly even for the scenario mentioned above. To me it appears that this sampling is more robust and should be generally preferred. Perhaps a flag similar to PhotoPionProduction::setSampleLog(bool b) should be available to enable this behavior or the same flag could be used for both purposes.

Extra Material

Photon Field Density & Sampling

In the following I show results for three different implementations of the CMB in CRPropa to illustrate the problem and show that my suggested solution works correctly. The photon field densities of the three fields, accessed via .getMinimumPhotonEnergy(0) for each field are shown in the figure below. They are in good agreement. Also indicated are 100 random samples of the photon energies for the tabulated CMB field with large epsMax (green dotted field; field 3) for three different implementations. It is clear that the linear sampling currently implemented in CRPropa fails to properly sample the entire energy range and therefore fails to "find any photons" as the original error I received suggest. If logarithmic sampling is used, either via scipy.stats.reciproval.rvs() or by manually implementing the inverse CDF of the reciprocal distribution, the entire range is sampled and the interval with non-zero photon density is found (note that the latter two implementations are not exactly identical in the plot because I am using separate random numbers for each approach).

A Simple Propagation Exercise

The following plot shows the fraction of surviving cosmic-ray protons after propagating them for 10 Mpc in the default and tabulated CMB fields respectively with only photopion production enabled. The results for the tabulated CMB field are not in agreement with the other implementations if the vanilla CRPropa version (with linear sampling) is used (essentially no interaction; finishes in seconds while the other runs take ~5min). It can be brought into excellent agreement with the other approaches if one switches to the logarithmic sampling of photon energies (modifying the source code and reinstalling CRPropa).

From the above survival fraction one can estimate the average photopion interaction length for protons with the CMB. This is illustrated in the plot below. The results for all scenarios but the field with high epsMax and linear sampling are in good agreement with each other and literature sources such as the book by Dermer&Menon 2009 (Fig. 9.3 dash-dotted line). Note that points in the 1e3-1e4 Mpc range are due to some simple numerical approximations in the plotting process and should be at infinity. They occur because the number of particles in in a specific energy bin is lower than what would be expected for an ideal powerlaw distribution in ~50% of cases. This is not related to the problem presented here and should be ignored.

[JulienDoerner]

Hey @ehlertdo,
The sampling is not directly a linear sampling. It uses the monte carlo rejection.
In this we randomly (uniform) choose a point and decide with a second random variable if it is taken into account for the decision if the point belongs to the distribution.

CRPropa3/src/module/PhotoPionProduction.cpp

Line 427 in 1ffe835

if (random.rand() * pEpsMax < pEps)

If you now increase the sampling range, by increasing the tabulated range, to an region where no contribution to the photon field is expected, this leads to a lot of tries which are rejected. This leads to the error-message you have recieved. And you will spend a lot of time, trying to get a photon.

If you have a way to reduce the problem of testing a rejected area it could be realy helpfull.
Do you have an implementation of this which we could test? Maybe you can give us a link or suggest the changes in a pull request.

[ehlertdo]

Hi @JulienDoerner
You are right; I did not understand the current sampling approach correctly. However, in any case the current version of CRPropa has problems with the scenario mentioned above.
I created a pull request (#476) to compare the proposed changes and also included the necessary files to reproduce the problem.