CRPropa / CRPropa3

CRPropa is a public astrophysical simulation framework for propagating extraterrestrial ultra-high energy particles. https://crpropa.github.io/CRPropa3/
https://crpropa.desy.de
GNU General Public License v3.0
69 stars 68 forks source link

Failing example notebook for secondary #467

Closed JulienDoerner closed 7 months ago

JulienDoerner commented 7 months ago

Currently the example notebooks secondaries/photons.ipynb and secondaries/secondary_photons.ipynb raise an error:

Exception in crpropa::ModuleList::run: 
std::bad_alloc

I traced this error down to sim.add(EMInverseComptonScattering(URB_Protheroe96(),True)). With changing the photon field from Proteroe96 to Nitu21 the error does not ocurre anymore. This is also related to the failings of the crpropa-example-test.

JulienDoerner commented 7 months ago

The photon field URB_Fixsen11 gives the same propblem.

rafaelab commented 7 months ago

I think the problem has to do with the limits of the tabulated data for some models. I've tried adding some ifs to fix it, with no success.

JulienDoerner commented 7 months ago

I found the error to be related to the sampled center of mass energy. It is possible to sample slighly larger energies than the upper limit of the tabulated range $(s > 10^{23} \ \mathrm{eV}^2)$. Therefore the tabulated range in the ICSSecondariesEnergyDistribution is not sufficient enough. Increasing this upper bound of the table by a factor 2 should be enough.

I will prepare a PR to provide this fix.

rafaelab commented 7 months ago

Exactly, @JulienDoerner. A good reference would be the minimum among all background photon energies (~1e-10 eV for URB) times the maximum energy we wish to treat (~1e22 eV), squared. That would be ~1e24 eV^2.

I found the error to be related to the sampled center of mass energy. It is possible to sample slighly larger energies than the upper limit of the tabulated range (s>1023 eV2). Therefore the tabulated range in the ICSSecondariesEnergyDistribution is not sufficient enough. Increasing this upper bound of the table by a factor 2 should be enough.

I will prepare a PR to provide this fix.

JulienDoerner commented 7 months ago

If I understand the code right, the problem can occure independently of the photon field (it only depends on the form of the distribution of the CDF). In the line https://github.com/CRPropa/CRPropa3/blob/bc17418444ad4b8fb6ece32a5249065bd46b7859/src/module/EMInverseComptonScattering.cpp#L184

the maximal kinetic energy can be larger than the maximum tabulated case ($s\mathrm{kin}^\mathrm{max} = 10^{0.05} \cdot s\mathrm{tab}^\mathrm{max} \approx 1.122 \cdot s_\mathrm{tab}^\mathrm{max}$). And in the next line we also add a $m_e^2 \ c^4$, but this is much smaller than the end of the table. Therefore the factor $2$ is already quite optimistic.