Sampling of Photon Fields in the PhotoPionProduction

[LeanderSchlegel]

Dear all,

with this pull request, we want to contribute our modifications of the photon field sampling in the PhotoPionProduction (PPP), unifying the sample approach and making it available for custom fields directly accessible through the call of the PPP-interaction module by the user. In this context, we also remove the silent defaulting of all IRB fields to the IRB_Kneiske04 in the PPP.

Sampling

The new unified sampling is equally applicable to analytical as well as tabulated photon fields. The functionality is fundamentally based on the previously implemented version of CMB sampling. The new sampling has 2 steps: 1. find the maximum interaction probability of all photons in the field. For this purpose, in linear or logarithmic steps, discrete probabilities of various photon energies are checked and the maximum value is stored. 2. randomly draw photons and determine their probabilities and compare them according to the Monte-Carlo principle with the previously determined maximum probability. We countered the resulting slowdown in our general approach by improving the other code parts in the sampling class and removing redundant or unnecessary, computationally intensive calculations.

The following plot shows a comparison between the histograms gained by sampling 1e5 photons of the CMB, assuming a proton as a primary with an energy of 10^10 GeV and a redshift z=0. The blue curve was calculated using the master version and served as a reference, the red curve was calculated using our modification and calling the now directly available analytic function CMB(), the green curve finally using a tabulated version of the CMB also using our modifications. While the table data, in this case, was self-generated due to a bug in the given table, it should serve as a test of the tabular framework, the sampling, and effectively the customization ability of this approach. All curves are in very good agreement within the statistics, showing that the code mechanism works as expected and should also do so for all other tabulated files like the provided IRBs or custom fields.

Custom fields

Custom fields can simply be used by providing the field data in three tabulated files in the share-directory, the list of energies, the list of densities, and optional a list of redshifts.

Syntax to use custom fields in the PhotoPionProduction:

field = TabularPhotonField(str fieldName, bool isRedShiftdependent) sim.add(PhotoPionProduction(field,photons,neutrinos,electrons,antiNucleons))

The implementation is fundamentally based on the existing TabularPhotonField class.

Optional settings

In general, the PhotonFieldSampler is a part of the PhotoPionProduction module. The module contains two variables that can be set for sufficient sampling especially with custom fields.

The first option is 'sampleLog'. It is by default set to True to use logarithmic steps for the search for the maximum of the probability function (for example needed for IRB_Kneiske04). The logarithmic scaling can deal with large ranges of photon energies. However, if linear sampling is sufficient, for larger energy ranges this option could be faster.

The second option is the 'correctionFactor' which increases the maximum of the probability. This value is set to f = 1.6 by default, which corresponds to the previous implementation in SOPHIA. For sharp, user-defined custom fields, it could be necessary to increase this value. This will lead to a better realization of the maximum of the photon distribution, but it increases the time needed for sampling a photon. In the following, you can see an example for the effect of this parameter. Here, the CMB is sampled for a proton with E = 10^10 GeV.

Example to access the option via python:

`import crpropa as crp

cmb = crp.CMB() ppp = crp.PhotoPionProduction(cmb) ppp.photonFieldSampling.setCorrectionFactor(1.8) ppp.photonFieldSampling.setSampleLog(False) `

Best regards and thank you in advance! @reichherzerp,@JulienDoerner, LeanderSchlegel, @Froehliche-Kernschmelze

[rafaelab]

Some comments: we have the CRPropaconventions. While most of the PR already complies with it, there are small things that don't: "function(){" should be "function() {", "1eV" should be "1 eV". These things are irrelevant, but I believe it is a good practice to keep a consistent coding style. (Sorry for being petty ;) ) I won't have time to properly test the code until next week, so hopefully someone else can provide feedback before that.

[rafaelab]

Hi @LeanderSchlegel, @JulienDoerner

Thanks for the great contribution. I'm done testing the PR. I checked the interaction rates and there are no major changes introduced by the PR. The main differences in the physics are indeed in the production of secondaries, as described above.

A minor thing: the data file in this version is being extracted to build/data-2021-07-30/data/. Ideally it should be just build/data to keep things consistent with previous versions.

As a general comment, why isn't the photon sampling for photopion production done in the photopion production module? Isn't the current setup breaking one of CRPropa's central principles -- modularity? This new photon sampling code is specific to PhotoPionProduction. Sure, photons are being sampled, but they are being sampled based on the kinematics of photopion production that should be done in PhotoPionProduction.

In the long run it would be nice to have general photon field sampling routines in PhotonBackground to replace, for instance, the CDFs used in the EM* modules. The photon field sampling from PhotoPionProduction seems to be very particular to this interaction.

Therefore, my suggestion is to copy-paste it into the appropriate place, unless there is a strong reason for not doing so.

[LeanderSchlegel]

Hi @rafaelab,

thank you for your testing of the PR! We completely agree, Patrick has already transferred the code to the PPP-module. I looked at the directory of the extracted datafiles but did not know how to extract them directly to /data, so I tried a small change in the CMakeLists.txt 634ed5391bac90a523c46ea6df2f6fdb140d9bc0 that removes the files afterwards into this folder and deletes the other now empty one that has the version number. To make sure ccmake can not run into errors when running multiple times, I had it check if the /data folder already exists and if the version numbered folder still exists. Would that be a good solution or is there a more elegant way to do this?

[LeanderSchlegel]

Just added a small test to deterministically test the sampling by calculating the maximum probability for the CMB and a proton of a certain energy. The test is very specific and could probably be a good indicator if something in the sampling would have been altered at some point.

[reichherzerp]

These enhancements introduce getters and setters for the PPP class. Tests for the cross-language polymorphism of custom photon fields and their usage in PPP have been implemented. These tests alert if functionalities related to custom photon fields in the PPP would be restricted in the Phyton interface.

[rafaelab]

Hi @LeanderSchlegel, @reichherzerp I've successfully tested this new version. I've made some comments directly in the code. My only question is if the quantile of 0.01% is sufficient to capture all physics, or if we need something like 1e-6.

[JulienDoerner]

hi @rafaelab, I think the lower quantile is sufficient enough. Till now it is only used in the PhotoPionProduction and it is always compared to the threshold energy for the interaction. The quantile will only be used for extreme values of primary energy and blackbody tempreture, as you can see in this plot:

But I have seen a problem with the upper boundary for high tempreture fields. I have changed it. Now the maximal photon energy is scaled with the tempreture, but at least 0.1 eV. This scaling is not physicaly motivated but seems to be good for the sampling algorithm. I have compared the maximal sampled photon energy of 1e5 photons to the maximal allowed energy in the algorithm.

[lukasmerten]

Please add a short description (~one line) to the CHANGELOG.

[LeanderSchlegel]

Just added the short description :)

[LeanderSchlegel]

Hi @rafaelab, thank you for your review! @JulienDoerner and I went through the code and fixed the code style, added short comments for the quantiles based on the tests and changed the CHANGELOG to your suggestion.

[adundovi]

Thanks @reichherzerp and @LeanderSchlegel for all the hard work, and @lukasmerten and @rafaelab for reviewing and testing. I've compiled and tested the PR manually and briefly reviewed the code - it works and I don't see any issues. I would like that we add a Jupyter notebook example at some point later, but with another PR.

However, please, resolve first conflicts with the current main/master branch so we can merge it.

[LeanderSchlegel]

Thank you @adundovi for your review and testing of the code! If I see it correctly, Github says there are no conflicts left with the base branch, are there changes we still need to do?

[adundovi]

Thank you @adundovi for your review and testing of the code! If I see it correctly, Github says there are no conflicts left with the base branch, are there changes we still need to do?

Sorry, I've mixed up merge with rebase so I notice conflicts in the rebase case. Now it is merged.