Replace SOPHIA by C++ version, closes #213, potentially closes #85

[MHoerbe]

Hey there, I have a working C++ v.11 version of SOPHIA. The .cpp file exists, yet I require help in embedding the code into CRPropa shall the decision be made to approve my proposal. This particularly includes the cmake files, etc.

The FORTRAN version of SOPHIA as-is is known to produce a computational bottleneck if parallelized. This causes the so-called "8-core saturation problem" in CRPropa as described in this paper, Fig. 1.

There are approaches to discretize SOPHIA as suggested in #213 but this would come at the expense of losing precision, particularly in less common parameter constellations which may not have been sampled and where interpolation is not meaningful. My approach is to replace the FORTRAN version of SOPHIA by the C++ version I have created.

I have created this version using the SOPHIA-internal stand-alone random number generator also in C++, thus being able to fully track the computational status of variables inside the two versions of SOPHIA. As a result I can show that in many parameter constellations of input parameters, the C++ and FORTRAN versions produce the exact same results, minus double precision deviations. Further improvements include:

• all variables now get properly declared as opposed to implicit declarations in FORTRAN
• all goto pointers have been removed and replaced by proper loops. The only exception is the method LUSTRF which could be rendered goto-less yet at the expense of rather unreadable code
• the amount of global variables has been reduced to improve readablity
• all input and output parameters of functions are now clearly able to be identified by a reader
• SOPHIA could be set to consider e.g. charged pions stable for things like pion propagation in CRPropa

SOPHIA as-is takes three input parameters: 1) bool onProton 2) double E_primary 3) double E_photon. The scattering angle is sampled internally from these input parameters. Each of the following parameter constellations have been tested and in the C++ version proven to reproduce the same results as in FORTRAN up to 100.000 event calls each:

• onProton = {true, false}
• log10 (E_primary / GeV) = {6, 7, 8, 9, 10, 11, 12}
• log10 (E_photon / eV) = {-4, -3, -2, -1, 0, 1, 2}

"same results" means here that if the threshold for photo-nuclear interactions is surpassed with a set of input parameters, that the same 1) number of particles with the 2) same types of particles and the 3) same 4-momenta are being produced. In some cases, there are deviations in the order of a few eV due to double precision issues.

The runtime performances of the FORTRAN and C++ versions of SOPHIA are summarized in the subsequent plot. The x-axis describes the product of E_primary x E_photon to better compare the data points to the input parameters above and almost is the CM energy of the interacting particles. Note that the higher CM energies have not been called in SOPHIA so far since the current implementations only allow for interactions with the CMB or IRB_Kneiske04 field whose photon energies are in the order of a few if not less meV.

What I suggest doing:

1. Add this version to the PhotonPionProduction, keep SOPHIA in the FORTRAN version and add a switch bool useSOPHIACplusplus or the like
2. (conduct more tests and) consider the C++ version "as good as FORTRAN" thus removing the old version of SOPHIA.

In either cases I would require help to embed the C++ file into CRPropa.

Please let me know your thoughts, Mario