Broken Fermi PSF example: examples/example_fermi_psf.py

[cdeil]

@ZiweiOu just mentioned at #1211 (comment) that examples/example_fermi_psf.py is currently broken:

$ python examples/example_fermi_psf.py 
Traceback (most recent call last):
  File "examples/example_fermi_psf.py", line 22, in <module>
    kernel = psf.kernel(pixel_size=pixel_size, offset_max=offset_max)
TypeError: kernel() got an unexpected keyword argument 'pixel_size'

The reason is that the TablePSF.kernel method was changed to work like this:
http://docs.gammapy.org/en/latest/api/gammapy.irf.TablePSF.html#gammapy.irf.TablePSF.kernel
and we forgot to update the example.

Also - as discussed via email: @gallanty has recently found a bug in IACTBasicImageEstimator.psf where it evaluates the PSF onto a sky image, which in his case wasn't centered properly and then morphology fitting results from Sherpa were incorrect. The bug is in IACTImageEstimator.psf where the PSF is just evaluated onto a sky image instead of calling TablePSF.kernel:

gammapy/gammapy/image/basic.py

Line 215 in c421ba9

psf_mean = mean_psf.table_psf_in_energy_band(erange, spectrum=self.spectral_model)

The proper fix here is to introduce a small PSFKernel2D class that always contains odd-sized kernel arrays and to change the TablePSF.kernel method to properly fill that kernel array (PSF centered on the central pixel), and to then call that from IACTImageEstimator.psf and if there are a few more cases that use SkyImage to represent a PSF to adapt those callers.

I'm mentioning this here so that we have a record, and because the failing examples/example_fermi_psf.py is directly related and a very nice example script one can run while coding on PSFKernel2D and TablePSF.kernel.

@gallanty - Do you have time this week or next latest to work on this? Or should I do it or try to find someone?

[gallanty]

@cdeil : Thanks for the detailed suggestions about how to fix this problem. I am busy with other things right now, but should have time from this Thursday on to try and improve the PSF calculation along the lines you suggest.

[cdeil]

@gallanty - Great!

Some more infos to help you get started on this coding:

• start by making a new git branch and check here for infos on git / Python -- everything in Gammapy is the same as in Astropy.
• make new files gammapy/irf/psf_kernel.py and gammapy/irf/tests/test_psf_kernel.py where you add the implementation and tests for the new PSFKernel2D class
• just add what you need to the class, i.e. data members array (2-dim numpy array, odd size) and pixel_size (either a Quantity or a float in deg, as you prefer). Take these two things as arguments in __init__ and just store them away, it's up to callers (e.g. TablePSF.kernel) to compute them.
• then add a to_image_hdu which returns an ImageHDU object and a write(filename) method which calls to_image_hdu and then writes to a FITS file -- i.e. make it convenient to store the PSF kernel in a format that in addition to Gammapy, also Sherpa will understand.
• Change TablePSF.kernel to return this object and try to get examples/example_fermi_psf.py to work again.
• In principle, you should also run python setup.py test -V to find and fix all callers for TablePSF.kernel and adapt them. But that quickly becomes a large and confusing task, not a simple exercise for a first Gammapy pull request. As soon as your branch is on Github and you have made a pull request, we can collaborate on it, i.e. I and other Gammapy developers can also push commits to it.

[cdeil]

The PSF kernel class was introduced and the Fermi PSF docs example adapted to the new code:
https://github.com/gammapy/gammapy/blob/master/docs/irf/plot_fermi_psf.py

Closing this old issue now.