custom photon fields / SOPHIA update

[MHoerbe]

Dear all,

with this pull request I want to provide the first of three parts that I wish to contribute to CRPropa. This approach follows the comments made in #218

Content:

This extension is meant to introduce variable photon fields to CRPropa. Previously, there has only been the CMB and several IRB models. This photon field extension basically follows three concepts:

The photon field may

1. have any shape, i.e. any dependence of photon energy and field density
2. have any spacial and temporal behaviour (which is part of a future pull request)
3. be modified even after CRPropa has been installed

Changelog (most severe changes):

~/src/PhotonBackground.cpp ~/include/PhotonBackground.h

There are now eight new slots PF1..PF8 that may be included in the simulation. A class has been added to draw a random photon from the photon field that takes the primary, the nucleon/gamma crossection and the photon field's shape into account. This method and all those that it depends on have been inspired by SOPHIA.

Photon fields are now contained in photon field files in ~/share/crpropa/Scaling/ . For detailed information on the format and use of these fields, there is a Wiki page in my forked repository.

~/libs/sophia/sophia_interface.f

The main point of these changes to SOPHIA is to run the code on other photon fields than the CMB and IRB Kneiske (as only these two are implemented so far in SOPHIA). SOPHIA samples photons from these fields with its own photon sampling routine which in this version has been removed and replaced by the sampling method mentioned in the section above. For this, the input and output parameters of SOPHIA changed from

sophiaevent(nature, Ein, OutPart, OutPartType, NbOutPart, z_particle, bgFlag, Zmax_IRB, idatamax, en_data, flux_data)

to

sophiaevent(nature, Ein, OutPart, OutPartType, NbOut, eps).

Furthermore, SOPHIA has recieved:

• a removal of all parts of the program that are never used under any circumstances
• a removal of all GOTO pointers which have been replaced by safer control structures
• a unified indentation to ehance the readability of the code
• a re-naming of several variables and methods from cryptic abbreviations to full identifiers
• docstrings in most functions that as far as tracable describe their functions

~/src/module/PhotonPionProduction.cpp ~/include/module/PhotonPionProduction.h

The flag "haveRedshiftDependence" and all its associated functions have been removed as now all photon fields are represented on grids that are redshift-dependent. The signature of SOPHIA has been updated. The treatment of secondary particles being produced in SOPHIA has been rewritten in a more intuitively readable way. Also, this particle treatment did not allow for more than one hadron to be added among all secondaries produced by SOPHIA. This is valid for the CMB and IRB models as not eanough energy in the center of mass system would be present to produce nucleon/anti-nucleon pairs. However, with more energetic photon fields, these particles may occur. I would like to stress that no kind of anti-nucleon can be produced in the current implementation, even if it were produced.

Tests

As photon sampling is removed from SOPHIA and now provided within CRPropa, SOPHIA's sampling routines and its CRPropa counterparts are expected to produce similar results. Therefore, the following tests have been run according to their hierarchy in the procedure of sampling photons:

Hierarchy of sampling functions (highest to lowest):

sample_eps -> prob_eps -> get_photonDensity, gaussInt -> functs -> crossection

Reference simulation setup:

SimplePropagation(10 kpc, 100 kpc)
PhotoPionProduction(field=CMB, neutrinos=True)
Source @ (100 Mpc, 0, 0)										
Injected: 100.000 monoenergetic protons with 1e12 GeV

1. crossection is the deepest layer to be tested. Note that the SOPHIA function that returns the nucleon-gamma xsection does not depend on the nucleon's energy itself but rather on the fact whether or not it deals with a proton or a neutron. The reference simulation called the crossection function roughly 500.000 times. The sampled photon's energies have been injected into the newly written C++ version of crossection producing another set of ca 500.000 crossection values. These two samples have been inserted into the scipy.stats.ks_2samp Kolmogorov-Smirnov test. A result of D=0.0 and p=1.0 corresponds to total identity in this implementation of the Kolmogorov-Smirnov test. The test achieved a result of D=0.0011, p=0.9183. Note that the functions Ef, Pl and breitwigner are only used by crossection.

2. In an analogous way, the function gaussInt has been tested, which integrates the function functs over all center-of-mass energies. the function functs is only called by gaussint, whereas functs is the only function calling crossection. The KS test yielded a result of D=9.61e-5, p=0.9999999999 between the SOPHIA output and the C++ version.

3. prob_eps calculates the non-normalized probability to encounter a photon of the given photon field. Normalization is achieved by dividing the integral of prob_eps over its entire photon energy range. prob_eps is the only function calling gaussInt and get_photonDensity. For better comparison, get_photonDensity has been replaced by an analytical expression returning the photon density of the CMB, i.e. a blackbody with a temprature of 2.73 K. The results of SOPHIA vs C++ yield D=9.61e-5, p=0.9999999999. Results on grid data (and runtimes) are discussed in a later section.

4. Finally, sample_eps calls prob_eps to sample a photon. This happens in two steps: a) calculate the normalization and the photon energy with which one most likely gets an interaction. b) monte-carlo sample a photon from the field, use norm factor and max probability for this. For comparison, the SOPHIA output has been tested vs the C++ analytical expression for the photon density and vs the grid readout. The KS test yields:

SOPHIA vs C++ (analytic):     D=0.0224, p=0.0114
SOPHIA vs C++ (grid):         D=0.0225, p=0.0110
C++ (analytic) vs C++ (grid): D=0.0001, p=0.9999

Insights: the grid version and the analytic version are able to produce nearly identical outputs, however less similarity is achieved compared to SOPHIA. In return, SOPHIA operates a dedicated sampling routine that only works for the CMB this making several assumptions about the shape of the field in beforehand. The C++ version do not require any information that is not encoded in the photon field file. One might be interested in the fact that SOPHIA uses its own random number generator whereas the C++ routines have been run with std::srand().

5. Totally finally, the sampling method is implemented into CRPropa and tested on observables. For this, the neutrino output (energies) of the reference simulation has been chosen. The KS test result of a neutrino the neutrino outputs acquired by SOPHIA vs C++ sampling have been tested in three ways fore comparison: a) using no grid i.e. the analytical density readout b) density readout of a CMB grid with 101 points of resolution c) density readout of a CMB grid with 1001 points of resolution

SOPHIA vs C++ (analytic):  1:57 | 2:20 @ D=0.0035, p=0.3514
SOPHIA vs C++ (Grid_101):  1:57 | 2:20 @ D=0.0037, p=0.2716
SOPHIA vs C++ (Grid_1001): 1:57 | 5:58 @ D=0.0035, p=0.3514

Insights: The C++ small grid / analytic versions perform roughly 20% slower than SOPHIA. A larger grid may reach the analytic version at the expense of runtime. All C++ versions perform similarly but do not reach identity with the SOPHIA version. The two points coming to my mind are: a) the C++ version of sampling has been done on CRPropa's random number generator, SOPHIA uses its own one b) SOPHIA operates a dedicated sampling method that only works on the CMB (yet another one with IRB_Kneiske) making assumptions about the CMB shape in beforehand. The C++ version does not rely on such information.

Important notes:

1. all the tests have been run on the CMB as this allows for tests against analytical expressions of the photon field and SOPHIA. The other only photon field implemented in SOPHIA is IRB_Kneiske. First tests did not yield a correlation between SOPHIA and C++. Before going deeper here, I would like to invoke your feedback on the current implementation.
2. The code does not tell CRPropa to build the photon field files stored in ~/share/crpropa/Scaling upon installation. Also, the tabulated files of the PhotoPionProduction have changed. Do you have any advice in how to include this?
3. The photon field URB_Protheroe96 could not be converted into a photon field file as there is no sufficient photon field data in the CRPropa-data repository for this. Hence, this field is currently unavailable.

Thanks in advance and best wishes, Mario