Rewrite EnergyDispersion class

docs/irf/plot_energy_dispersion.py

@@ -1,11 +1,11 @@
 """Plot energy dispersion example."""
 import numpy as np
 import matplotlib.pyplot as plt
-from gammapy import irf
+from gammapy.irf import EnergyDispersion
+from gammapy.spectrum import EnergyBounds
 
-ebounds = np.logspace(-1, 2, 100)
-pdf_matrix = irf.gauss_energy_dispersion_matrix(ebounds, sigma=0.2)
-energy_dispersion = irf.EnergyDispersion(pdf_matrix, ebounds)
+ebounds = EnergyBounds.equal_log_spacing(0.1, 100, 100, 'TeV')
+energy_dispersion = EnergyDispersion.from_gauss(ebounds, sigma=0.3)
 
 energy_dispersion.plot(type='matrix')
 plt.show()

gammapy/data/counts_spectrum.py

@@ -6,6 +6,7 @@
 from ..spectrum.energy import Energy, EnergyBounds
 import datetime
 from astropy.io import fits
+from gammapy.data import EventList, EventListDataset
 
 __all__ = ['CountsSpectrum']
 
@@ -32,7 +33,7 @@ class CountsSpectrum(object):
         ebounds = EnergyBounds.equal_log_spacing(1,10,10,'TeV') 
         counts = [6,3,8,4,9,5,9,5,5,1]
         spec = CountsSpectrum(counts, ebounds)
-        hdu = spec.to_pha() 
+        hdu = spec.to_fits() 
     """
 
     def __init__(self, counts, energy, livetime=None):
@@ -52,6 +53,7 @@ def __init__(self, counts, energy, livetime=None):
             self.energy = energy
 
         self._livetime = livetime
+        self.channels = np.arange(1,self.energy.nbins+1,1)
 
     @property
     def entries(self):
@@ -93,19 +95,27 @@ def from_eventlist(cls, event_list, bins):
 
         if isinstance (event_list,fits.BinTableHDU):
             event_list = EventList.read(event_list)
-        elif isinstance(event_list, gammapy.data.EventListDataSet):
+        elif isinstance(event_list, EventListDataset):
             event_list = event_list.event_list
         elif isinstance(event_list, str):
             event_list = EventList.read(event_list, hdu='EVENTS') 
 
-        energy = event_list.energy.to(bins.unit)
-        val = np.histogram(energy,bins)
+        energy = Energy(event_list.energy).to(bins.unit)
+        val, dummy = np.histogram(energy,bins.value)
         livetime = event_list.observation_live_time_duration
 
-        return cls(val, energy, livetime)
+        return cls(val, bins, livetime)
 
+    def write(self, filename, bkg=None, corr=None, rmf=None, arf=None,
+              *args, **kwargs):
+        """Write PHA to FITS file.
 
-    def to_pha(self):
+        Calls `~astropy.io.fits.HDUList.writeto`, forwarding all arguments.
+        """
+        self.to_fits(bkg=bkg, corr=corr, rmf=rmf, arf=arf).writeto(
+            filename, *args, **kwargs)
+
+    def to_fits(self, bkg=None, corr=None, rmf=None, arf=None ):
         """Output OGIP pha file
 
         Returns
@@ -119,42 +129,47 @@ def to_pha(self):
         http://heasarc.gsfc.nasa.gov/docs/heasarc/ofwg/docs/summary/
         ogip_92_007_summary.html
         """
-        col1 = fits.Column(name="CHANNEL", array=self.energy.value,
-                          format = 'E', unit = '{0}'.format(self.energy.unit))
-        col2 = fits.Column(name="COUNTS", array=self.counts, format='E')
+        col1 = fits.Column(name="CHANNEL", array=self.channels,
+                           format = 'I')
+        col2 = fits.Column(name="COUNTS", array=self.counts,
+                           format='J', unit='count')
 
         cols = fits.ColDefs([col1, col2])
 
         hdu = fits.BinTableHDU.from_columns(cols)
         header = hdu.header
 
-        # copied from np_to_pha
-        # I don't know if it should move here or stay a utility function
+        #TODO: option to store meta info in the class
         header['EXTNAME'] = 'SPECTRUM', 'name of this binary table extension'
-
-        #header['BACKFILE'] = backfile, 'Background FITS file'
-        #header['CORRFILE'] = corrfile, 'Correlation FITS file'
-        #header['RESPFILE'] = respfile, 'Redistribution matrix file (RMF)'
-        #header['ANCRFILE'] = ancrfile, 'Ancillary response file (ARF)'
+        header['TELESCOP'] = 'DUMMY', 'Telescope (mission) name'
+        header['INSTRUME'] = 'DUMMY', 'Instrument name'
+        header['FILTER'] = 'NONE', 'Instrument filter in use'
+        header['EXPOSURE'] = self.livetime.to('second').value, 'Exposure time'
+
+        header['BACKFILE'] = bkg, 'Background FITS file'
+        header['CORRFILE'] = corr, 'Correlation FITS file'
+        header['RESPFILE'] = rmf, 'Redistribution matrix file (RMF)'
+        header['ANCRFILE'] = arf, 'Ancillary response file (ARF)'
 
         header['HDUCLASS'] = 'OGIP', 'Format conforms to OGIP/GSFC spectral standards'
         header['HDUCLAS1'] = 'SPECTRUM', 'Extension contains a spectrum'
         header['HDUVERS '] = '1.2.1', 'Version number of the format'
 
         header['CHANTYPE'] = 'PHA', 'Channels assigned by detector electronics'
         header['DETCHANS'] = self.energy.nbins, 'Total number of detector channels available'
-        header['TLMIN1'] = self.energy[0].value, 'Lowest Legal channel number'
-        header['TLMAX1'] = self.energy[-1].value, 'Highest Legal channel number'
+        header['TLMIN1'] = 1, 'Lowest Legal channel number'
+        header['TLMAX1'] = self.energy.nbins, 'Highest Legal channel number'
 
         header['XFLT0001'] = 'none', 'XSPEC selection filter description'
 
         header['HDUCLAS2'] = 'NET', 'Extension contains a bkgr substr. spec.'
         header['HDUCLAS3'] = 'COUNT', 'Extension contains counts'
         header['HDUCLAS4'] = 'TYPE:I', 'Single PHA file contained'
         header['HDUVERS1'] = '1.2.1', 'Obsolete - included for backwards compatibility'
-
+
+        header['POISSERR'] = True, 'Are Poisson Distribution errors assumed'
+        header['STAT_ERR'] = 0, 'No statisitcal error was specified'
         header['SYS_ERR'] = 0, 'No systematic error was specified'
-
         header['GROUPING'] = 0, 'No grouping data has been specified'
 
         header['QUALITY '] = 0, 'No data quality information specified'

gammapy/data/event_list.py

@@ -249,8 +249,7 @@ def select_sky_cone(self, center, radius):
         mask = separation < radius
         return self[mask]
 
-    def select_sky_ring(self, center, inner_radius=None, outer_radius=None,
-                        thickness=None):
+    def select_sky_ring(self, center, inner_radius, outer_radius):
         """Select events in sky circle.
 
         Parameters
@@ -259,38 +258,21 @@ def select_sky_ring(self, center, inner_radius=None, outer_radius=None,
             Sky ring center
         inner_radius : `~astropy.coordinates.Angle`
             Sky ring inner radius
-        outer_radius : `~astropy.coordinates.Angle`, None
+        outer_radius : `~astropy.coordinates.Angle`
             Sky ring outer radius
-        inner_radius : `~astropy.coordinates.Angle`
-            Sky ring thinkness
 
         Returns
         -------
         event_list : `EventList`
             Copy of event list with selection applied.
         """
 
-        if outer_radius is None:
-            if thickness is None or inner_radius is None:
-                raise ValueError("Not enough information specified")
-            else:
-                outer_radius = inner_radius + thickness
-
-        else:
-            if inner_radius is None:
-                if thickness is None:
-                    raise ValueError("Not enough information specified")
-                else:
-                    inner_radius = outer_radius + thickness
-
-            else:
-                thickness = outer_radius - inner_radius
-
         position = self.radec
         separation = center.separation(position)
-        mask1 = inner_radius > separation
-        mask2 = outer_radius < separation
+        mask1 = inner_radius < separation
+        mask2 = separation < outer_radius
         mask = mask1 * mask2
+
         return self[mask]
 
     def select_sky_box(self, lon_lim, lat_lim, frame='icrs'):

gammapy/hspec/wstat.py

@@ -10,6 +10,7 @@ class spectral_data(object):
     """
 
     def __init__(self, data, bkg):
+
         self.bkg = bkg
         self.data = data
         self.ONexpo = data.exposure
@@ -61,6 +62,7 @@ def __init__(self, dataids):
         for dataid in dataids:
             spec = spectral_data(sau.get_data(dataid), sau.get_bkg(dataid))
             self.datasets.append(spec)
+
         self.bkg = np.concatenate([a.fit_bkg for a in self.datasets])
         self.alpha = np.concatenate([a.fit_alpha for a in self.datasets])
         self.Ttot = np.concatenate([a.fit_Ttot for a in self.datasets])

gammapy/irf/effective_area_table.py

@@ -198,12 +198,14 @@ def to_fits(self, header=None, energy_unit='TeV', effarea_unit='m2', **kwargs):
         hdu.add_datasum()
         return fits.HDUList([prim_hdu, hdu])
 
-    def write(self, filename, *args, **kwargs):
+    def write(self, filename, energy_unit='TeV', effarea_unit='m2',
+              *args, **kwargs):
         """Write ARF to FITS file.
 
         Calls `~astropy.io.fits.HDUList.writeto`, forwarding all arguments.
         """
-        self.to_fits().writeto(filename, *args, **kwargs)
+        self.to_fits(energy_unit=energy_unit, effarea_unit=effarea_unit).writeto(
+            filename, *args, **kwargs)
 
     @classmethod
     def read(cls, filename):
@@ -482,8 +484,8 @@ def to_effective_area_table(self, offset, energy_lo=None, energy_hi=None):
         """
 
         if energy_lo is None and energy_hi is None:
-            energy_lo = self.energy_lo
-            energy_hi = self.energy_hi
+            energy_lo = self.energ_lo
+            energy_hi = self.energ_hi
         elif energy_lo is None or energy_hi is None:
             raise ValueError("Only 1 energy vector given, need 2")
         if not isinstance(energy_lo, Quantity) or not isinstance(energy_hi, Quantity):