Add TablePSF and Fermi PSF

docs/irf/effective_area.rst

@@ -0,0 +1,8 @@
+Effective area
+==============
+
+`~gammapy.irf.abramowski_effective_area`
+
+`~gammapy.irf.arf_to_np` and `~gammapy.irf.np_to_arf`
+
+.. plot:: irf/plot_effective_area.py

docs/irf/energy_dispersion.rst

@@ -0,0 +1,10 @@
+Energy dispersion
+=================
+
+`~gammapy.irf.EnergyDispersion`
+
+`~gammapy.irf.gauss_energy_dispersion_matrix`
+
+`~gammapy.irf.np_to_rmf`
+
+.. plot:: irf/plot_energy_dispersion.py

docs/irf/index.rst

@@ -10,13 +10,66 @@ Introduction
 `gammapy.irf` contains functions and classes to access and model
 instrument response functions (IRFs).
 
-TODO: document
+Theory
+------
+
+For high-level gamma-ray data analysis (measuring morphology and spectra of sources)
+a canonical detector model is used, where the gamma-ray detection process is simplified
+as being fully characterized by the following three "instrument response functions":
+
+* Effective area :math:`A(p, E)` (unit: :math:`m^2`)
+* Point spread function :math:`PSF(p'|p, E)` (unit: :math:`sr^{-1}`)
+* Energy dispersion :math:`D(E'|p, E)` (unit: :math:`TeV^{-1}`)
+
+The effective area represents the gamma-ray detection efficiency,
+the PSF the angular resolution and the energy dispersion the energy resolution
+of the instrument.
+
+The full instrument response is given by
+
+.. math::
+
+   R(p', E'|p, E) = A(p, E) \times PSF(p'|p, E) \times D(E'|p, E),
+
+where :math:`p` and :math:`E` are the true gamma-ray position and energy
+and :math:`p'` and :math:`E'` are the reconstructed gamma-ray position and energy.
+
+The instrument function relates sky flux models to expected observed counts distributions via
+
+.. math::
+
+   N(p', E') = t_{obs} \int_E \int_\Omega R(p', E'|p, E) \times F(p, E) dp dE,
+
+where :math:`F`, :math:`R`, :math:`t_{obs}` and :math:`N` are the following quantities:
+
+* Sky flux model :math:`F(p, E)` (unit: :math:`m^{-2} s^{-1} TeV^{-1} sr^{-1}`)
+* Instrument response :math:`R(p', E'|p, E)` (unit: :math:`m^2 TeV^{-1} sr^{-1}`)
+* Observation time: :math:`t_{obs}` (unit: :math:`s`)
+* Expected observed counts model :math:`N(p', E')` (unit: :math:`sr^{-1} TeV^{-1}`)
+
+If you'd like to learn more about instrument response functions, have a look at the descriptions for
+`Fermi <http://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_LAT_IRFs/index.html>`__,
+for `TeV data analysis <http://inspirehep.net/record/1122589>`__ 
+and for `GammaLib <http://gammalib.sourceforge.net/user_manual/modules/obs.html#handling-the-instrument-response>`__.
+
+TODO: add an overview of what is / isn't available in Gammapy.
 
 Getting Started
 ===============
 
 TODO: document
 
+
+Using `gammapy.image`
+=====================
+
+.. toctree::
+   :maxdepth: 1
+
+   effective_area
+   psf
+   energy_dispersion
+
 Reference/API
 =============
 

docs/irf/plot_effective_area.py

@@ -0,0 +1,21 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Plot approximate effective area for HESS, HESS2 and CTA.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy.units import Quantity
+from gammapy.irf import abramowski_effective_area
+
+energy = Quantity(np.logspace(-3, 3, 100), 'TeV')
+
+for instrument in ['HESS', 'HESS2', 'CTA']:
+    a_eff = abramowski_effective_area(energy, instrument)
+    plt.plot(energy.value, a_eff.value, label=instrument)
+
+plt.loglog()
+plt.xlabel('Energy (TeV)')
+plt.ylabel('Effective Area (cm^2)')
+plt.xlim([1e-3, 1e3])
+plt.ylim([1e3, 1e12])
+plt.legend(loc='best')
+plt.show()

docs/irf/plot_energy_dispersion.py

@@ -0,0 +1,14 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Plot energy dispersion example.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+#from astropy.units import Quantity
+from gammapy import irf
+
+ebounds = np.logspace(-1, 2, 100)
+pdf_matrix = irf.gauss_energy_dispersion_matrix(ebounds, sigma=0.2)
+energy_dispersion = irf.EnergyDispersion(pdf_matrix, ebounds)
+
+energy_dispersion.plot(type='matrix')
+plt.show()

docs/irf/plot_fermi_psf.py

@@ -0,0 +1,32 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Plot Fermi PSF.
+"""
+import numpy as np
+import matplotlib.pyplot as plt
+from astropy.coordinates import Angle
+from astropy.units import Quantity
+from gammapy.datasets import FermiGalacticCenter
+from gammapy.irf import EnergyDependentTablePSF
+
+filename = FermiGalacticCenter.filenames()['psf']
+fermi_psf = EnergyDependentTablePSF.read(filename)
+
+energies = Quantity([1], 'GeV')
+for energy in energies:
+    print('0')
+    psf = fermi_psf.table_psf_at_energy(energy=energy)
+    print('1')
+    psf.normalize()
+    print('2')
+    print('3', psf.eval(Angle(0.1, 'deg')))
+    kernel = psf.kernel(pixel_size=Angle(0.1, 'deg'))
+    print('4')
+
+    print(np.nansum(kernel.value))
+    plt.imshow(np.log(kernel.value))
+    #psf.plot_psf_vs_theta()
+
+#plt.xlim(1e-2, 10)
+#plt.gca().set_xscale('linear')
+#plt.gca().set_yscale('linear')
+plt.show()

docs/irf/psf.rst

@@ -0,0 +1,9 @@
+Point spread function
+=====================
+
+The `~gammapy.irf.TablePSF` and `~gammapy.irf.EnergyDependentTablePSF` classes
+represent radially-symmetric PSFs where the PSF is given at a number of offets:
+
+.. plot:: irf/plot_fermi_psf.py
+
+TODO: discuss noise when PSF is estimated from real or simulated data.

examples/fermi_psf.py

@@ -0,0 +1,39 @@
+# Licensed under a 3-clause BSD style license - see LICENSE.rst
+"""Compute Fermi PSF image for a given energy band.
+
+The input should 
+"""
+import numpy as np
+from astropy.coordinates import Angle
+from astropy.units import Quantity
+from astropy.io import fits
+from gammapy.datasets import FermiGalacticCenter
+from gammapy.irf import EnergyDependentTablePSF
+
+# Parameters
+filename = FermiGalacticCenter.filenames()['psf']
+pixel_size = Angle(0.1, 'deg')
+offset_max = Angle(2, 'deg')
+energy = Quantity(10, 'GeV')
+energy_band = Quantity([10, 500], 'GeV')
+outfile = 'fermi_psf_image.fits'
+
+# Compute PSF image
+fermi_psf = EnergyDependentTablePSF.read(filename)
+#psf = fermi_psf.table_psf_at_energy(energy=energy)
+psf = fermi_psf.table_psf_in_energy_band(energy_band=energy_band, spectral_index=2.5)
+psf.normalize()
+kernel = psf.kernel(pixel_size=pixel_size, offset_max=offset_max)
+kernel_image = kernel.value
+
+kernel_image_integral = kernel_image.sum() * pixel_size.to('radian').value ** 2
+print('Kernel image integral: {0}'.format(kernel_image_integral))
+print('shape: {0}'.format(kernel_image.shape))
+
+#import IPython; IPython.embed()
+# Print some info and save to FITS file
+#print(fermi_psf.info())
+
+print(psf.info())
+print('Writing {0}'.format(outfile))
+fits.writeto(outfile, data=kernel_image, clobber=True)

gammapy/irf/__init__.py

@@ -2,4 +2,6 @@
 """Instrument response function (IRF) functionality.
 """
 from .effective_area import *
+from .psf_core import *
+from .psf_table import *
 from .energy_dispersion import *

gammapy/irf/effective_area.py

@@ -18,14 +18,14 @@ def abramowski_effective_area(energy, instrument='HESS'):
 
     Parameters
     ----------
-    energy : array_like
-        Energy in TeV
-    instruments : {'HESS', 'HESS2', 'CTA'}
-        Name of the instrument
+    energy : `~astropy.units.Quantity`
+        Energy
+    instrument : {'HESS', 'HESS2', 'CTA'}
+        Instrument name
 
     Returns
     -------
-    effective_area : array
+    effective_area : `~astropy.units.Quantity`
         Effective area in cm^2
     """
     # Put the parameters g in a dictionary.
@@ -36,18 +36,21 @@ def abramowski_effective_area(energy, instrument='HESS'):
             'HESS2': [2.05e9, 0.0891, 1e5],
             'CTA': [1.71e11, 0.0891, 1e5]}
 
-    energy = Quantity(energy, 'TeV').to('MeV').value
+    if not isinstance(energy, Quantity):
+        raise ValueError("energy must be a Quantity object.")
+
+    energy = energy.to('MeV').value
 
     if not instrument in pars.keys():
         ss = 'Unknown instrument: {0}\n'.format(instrument)
         ss += 'Valid instruments: HESS, HESS2, CTA'
         raise ValueError(ss)
-   
+
     g1 = pars[instrument][0]
     g2 = pars[instrument][1]
     g3 = -pars[instrument][2]
     value = g1 * energy ** (-g2) * np.exp(g3 / energy)
-    return value
+    return Quantity(value, 'cm^2')
 
 
 def np_to_arf(effective_area, energy_bounds, telescope='DUMMY',

gammapy/morphology/psf.py → gammapy/irf/psf_core.py

@@ -19,9 +19,9 @@
 import numpy as np
 from numpy import log, exp
 from astropy.convolution import Gaussian2DKernel
-from .utils import read_json
-from .gauss import Gauss2DPDF, MultiGauss2D
 from ..utils.const import sigma_to_fwhm, fwhm_to_sigma
+from ..morphology import read_json
+from ..morphology import Gauss2DPDF, MultiGauss2D
 
 __all__ = ['GaussPSF',
            'HESSMultiGaussPSF',