Improve gammapy.spectrum.cosmic_ray_flux

[cdeil]

We have gammapy.spectrum.cosmic_ray_flux which needs to be removed or improved.

It's not well-tested, see coverage report.

This is a bit related to the case of crab spectrum and the API question there, see #1788 .

Can SpectralModel represent fluxes per solid angle, i.e. CR surface brightness in m^-2 s^-1 TeV^-1 sr^-1?

If yes, we could return SpectralModel objects here and use the existing classes. If no, then we could keep as-is and evaluate directly and just improve docs and tests.

But the first question is whether to remove completely, or whether to keep and polish.
-> I'm +1 to keep and polish and am volunteering to do this for v0.13 in this case, because I think the CR electron and proton spectrum is quite commonly used by gamma-ray people (e.g. I used the electron spectrum once to do a rough estimate if HESS can see the Fermi bubbles, that's when I added this).

@adonath @registerrier - Thoughts?

[adonath]

I agree the cosmic_ray_flux method might be useful. But it should be adapted to return an instance of SpectralModel. Technically all our spectral models can represent flux per solid angle, so I think this is not a blocker.

[cdeil]

Looks like some methods do work, and some don't for spectral models per solid angle:

>>> import astropy.units as u
>>> from gammapy.spectrum.models import PowerLaw
>>> model = PowerLaw(amplitude="1 cm-2 s-1 TeV-1 sr-1")
>>> model(1 * u.TeV)
<Quantity 1. 1 / (cm2 s sr TeV)>
>>> model.integral(1*u.TeV, 10*u.TeV)
<Quantity 0.9 1 / (cm2 s sr)>
>>> model.plot(energy_range=[1, 10]*u.TeV)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/Users/deil/work/code/gammapy/gammapy/spectrum/models.py", line 280, in plot
    flux = self(energy).to(flux_unit)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py", line 671, in to
    return self._new_view(self._to_value(unit, equivalencies), unit)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/quantity.py", line 643, in _to_value
    equivalencies=equivalencies)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 989, in to
    return self._get_converter(other, equivalencies=equivalencies)(value)
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 920, in _get_converter
    raise exc
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 906, in _get_converter
    self, other, self._normalize_equivalencies(equivalencies))
  File "/Users/deil/software/anaconda3/envs/gammapy-dev/lib/python3.7/site-packages/astropy/units/core.py", line 890, in _apply_equivalencies
    unit_str, other_str))
astropy.units.core.UnitConversionError: '1 / (cm2 s sr TeV)' and '1 / (cm2 s TeV)' are not convertible

So I'm not sure returning a SpectralModel for this is a good idea.

@adonath - Thoughts?

[adonath]

@cdeil The .plot() method uses default units for SEDs. I don't see a problem when this fails, either users have to set the flux_unit argument or we have to change the code to take the default unit from the amplitude parameter or similar. I'm still in favour of returning a SpectralModel object...

[registerrier]

I am also +1 for that solution.
I am not sure that there is a solid use case for CR spectra as SpectralModel, since we would not fold them. But having SpectralModel able to deal with intensities seems useful e.g. Diffuse Galactic and Extragalactic backgrounds.

[cdeil]

OK, I'll try to change to return a SpectralModel for this and open a PR soon.

[cdeil]

Actually it's unlikely that I'll get to this before going on vacation tomorrow. So I'm un-assigning myself here and will ask for volunteers for this in the weekly dev call tomorrow.

To summarise the outcome of the discussion above: the task here is to make a PR that updates gammapy/spectrum/cosmic_ray.py and gammapy/spectrum/tests/test_cosmic_ray.py.

The function cosmic_ray_flux should be renamed to cosmic_ray_spectrum and rewritten to return a spectral model, using the exising Gammapy spectral model classes (see here) and deleting the current _power_law and _log_normal helper functions.

The tests should execute cosmic_ray_spectrum once for each particle type, and then call evaluate or integral for the returned spectrum for one energy or energy range, and assert on the resulting value, possibly using pytest.parametrize as in gammapy/spectrum/tests/test_crab.py
To execute the tests locally and to check code coverage, the following command can be used:

pytest -v gammapy/spectrum/tests/test_cosmic_ray.py --cov=gammapy.spectrum.cosmic_ray --cov-report=html

If anyone is interested in this task, please leave a comment here any time.

[cdeil]

This was done by @JouvinLea in #2282

Following the suggestion in #1788 (comment) we now have create_crab_spectral_model

I will do a follow-up commit in master to change cosmic_ray_spectrum -> create_cosmic_ray_spectral_model so that we get a more uniform API.

[cdeil]

Done in b09efaf