Merge pull request #1035 from cdeil/fp

Some cleanup of FluxPoints code and tests

docs/spectrum/index.rst

@@ -78,7 +78,6 @@ Content
    :maxdepth: 1
 
    fitting
-   flux_points
    energy_group
    plotting_fermi_spectra
 

gammapy/catalog/fermi.py

@@ -17,7 +17,7 @@
 from ..utils.table import table_standardise_units_inplace
 from ..image import SkyImage
 from ..image.models import Delta2D, Template2D
-from ..spectrum import FluxPoints, compute_flux_points_dnde
+from ..spectrum import FluxPoints
 from ..spectrum.models import (
     PowerLaw,
     PowerLaw2,
@@ -453,9 +453,9 @@ def flux_points(self):
 
         flux_points = FluxPoints(table)
 
-        flux_points_dnde = compute_flux_points_dnde(
-            flux_points, model=self.spectral_model)
-        return flux_points_dnde
+        # TODO: change this and leave it up to the caller to convert to dnde
+        # See https://github.com/gammapy/gammapy/issues/1034
+        return flux_points.to_sed_type('dnde', model=self.spectral_model)
 
     @property
     def spectral_model(self):
@@ -526,9 +526,9 @@ def flux_points(self):
 
         flux_points = FluxPoints(table)
 
-        flux_points_dnde = compute_flux_points_dnde(
-            flux_points, model=self.spectral_model)
-        return flux_points_dnde
+        # TODO: change this and leave it up to the caller to convert to dnde
+        # See https://github.com/gammapy/gammapy/issues/1034
+        return flux_points.to_sed_type('dnde', model=self.spectral_model)
 
     @property
     def spectral_model(self):

gammapy/scripts/tests/test_spectrum_pipe.py

@@ -30,16 +30,16 @@ def get_config():
 @requires_dependency('scipy')
 @requires_dependency('sherpa')
 def test_spectrum_analysis_iact(tmpdir):
-    conf = get_config()
-    conf['outdir'] = tmpdir
+    config = get_config()
+    config['outdir'] = tmpdir
 
-    ana = SpectrumAnalysisIACT(
-        observations=obs_list(),
-        config=conf
-    )
+    analysis = SpectrumAnalysisIACT(observations=obs_list(), config=config)
+    analysis.run()
+    flux_points = analysis.flux_point_estimator.flux_points
+
+    print(flux_points)
+    print(flux_points.table)
 
-    assert 'IACT' in str(ana)
-    ana.run()
-    actual = ana.flux_point_estimator.flux_points.table[0]['dnde']
-    desired = 7.208219787928114e-08
+    actual = flux_points.table['dnde'].quantity[0]
+    desired = 7.208219780210792e-12 * u.Unit('cm-2 s-1 TeV-1')
     assert_quantity_allclose(actual, desired)

gammapy/spectrum/__init__.py

@@ -9,7 +9,6 @@
 from .diffuse import *
 from .energy_group import *
 from .flux_point import *
-from .sherpa_chi2asym import *
 from .utils import *
 from .extract import *
 from .simulation import *

gammapy/spectrum/crab.py

@@ -1,17 +1,4 @@
 # Licensed under a 3-clause BSD style license - see LICENSE.rst
-"""Published Crab nebula reference spectra.
-
-The Crab is often used as a standard candle in gamma-ray astronomy.
-Statements like "this source has a flux of 10 % Crab" or
-"our sensitivity is 2 % Crab" are common.
-
-Here we provide a reference of what the Crab flux actually is according
-to several publications:
-
-* 'HEGRA' : http://adsabs.harvard.edu/abs/2000ApJ...539..317A
-* 'HESS PL' and 'HESS ECPL' : http://adsabs.harvard.edu/abs/2006A%26A...457..899A
-* 'Meyer' :  http://adsabs.harvard.edu/abs/2010A%26A...523A...2M, Appendix D
-"""
 from __future__ import absolute_import, division, print_function, unicode_literals
 import numpy as np
 from astropy import units as u
@@ -23,39 +10,27 @@
 ]
 
 # HESS publication: 2006A&A...457..899A
-#'int_flux' = 2.26e-11
 hess_pl = {'amplitude': 3.45e-11 * u.Unit('1 / (cm2 s TeV)'),
            'index': 2.63,
-           'reference' : 1 * u.TeV}
+           'reference': 1 * u.TeV}
 
-# Note that for ecpl, the diff_flux is not
-# the differential flux at 1 TeV, that you
-# get by multiplying with exp(-e / cutoff)
-# int_flux: 2.27e-11
 hess_ecpl = {'amplitude': 3.76e-11 * u.Unit('1 / (cm2 s TeV)'),
              'index': 2.39,
              'lambda_': 1 / (14.3 * u.TeV),
              'reference': 1 * u.TeV}
 
 # HEGRA publication : 2004ApJ...614..897A
-# int_flux': 1.75e-11
 hegra = {'amplitude': 2.83e-11 * u.Unit('1 / (cm2 s TeV)'),
          'index': 2.62,
          'reference': 1 * u.TeV}
 
-# Meyer et al. publication: 2010arXiv1008.4524M
-# diff_flux and index were calculated numerically
-# by hand at 1 TeV as a finite differene
-meyer = {'diff_flux': 3.3457e-11 * u.Unit('1 / (cm2 s TeV)'),
-         'index': 2.5362,
-         'int_flux': 2.0744e-11}
 
-
-#TODO: make this a general LogPolynomial spectral model
 class MeyerCrabModel(SpectralModel):
+    """Meyer 2010 log polynomial Crab spectral model.
+    
+    See 2010A%26A...523A...2M, Appendix D.
     """
-    Log polynomial model as used by 2010arXiv1008.4524M.
-    """
+
     def __init__(self):
         coefficients = np.array([-0.00449161, 0, 0.0473174, -0.179475,
                                  -0.53616, -10.2708])
@@ -73,15 +48,16 @@ def evaluate(energy, coefficients):
 
 
 class CrabSpectrum(object):
-    """
-    Crab spectral model.
+    """Crab nebula spectral model.
+
+    The Crab nebula is often used as a standard candle in gamma-ray astronomy.
+    Fluxes and sensitivities are often quoted relative to the Crab spectrum.
 
     The following references are available:
 
-    * 'meyer', 2010arXiv1008.4524M
-    * 'hegra', 2004ApJ...614..897A
-    * 'hess_pl', 2006A&A...457..899A
-    * 'hess_ecpl', 2006A&A...457..899A
+    * 'meyer', http://adsabs.harvard.edu/abs/2010A%26A...523A...2M, Appendix D
+    * 'hegra', http://adsabs.harvard.edu/abs/2000ApJ...539..317A
+    * 'hess_pl' and 'hess_ecpl': http://adsabs.harvard.edu/abs/2006A%26A...457..899A
 
     Parameters
     ----------
@@ -90,44 +66,42 @@ class CrabSpectrum(object):
 
     Examples
     --------
-    Plot the 'hess_ecpl' reference crab spectrum between 1 TeV and 100 TeV:
-
-        >>> from gammapy.spectrum import CrabSpectrum
-        >>> import astropy.units as u
-        >>> crab_hess_ecpl = CrabSpectrum('hess_ecpl')
-        >>> crab_hess_ecpl.model.plot([1, 100] * u.TeV)
-
-    Use a reference crab spectrum as unit to measure flux:
-
-        >>> from astropy import units as u
-        >>> from gammapy.spectrum import CrabSpectrum
-        >>> from gammapy.spectrum.models import PowerLaw
-        >>>
-        >>> # define power law model
-        >>> amplitude = 1e-12 * u.Unit('1 / (cm2 s TeV)')
-        >>> index = 2.3
-        >>> reference = 1 * u.TeV
-        >>> pwl = PowerLaw(index, amplitude, reference)
-        >>>
-        >>> # define crab reference
+    Let's first import what we need::
+
+        import astropy.units as u
+        from gammapy.spectrum import CrabSpectrum
+        from gammapy.spectrum.models import PowerLaw
+
+    Plot the 'hess_ecpl' reference Crab spectrum between 1 TeV and 100 TeV::
+
+        crab_hess_ecpl = CrabSpectrum('hess_ecpl')
+        crab_hess_ecpl.model.plot([1, 100] * u.TeV)
+
+    Use a reference crab spectrum as unit to measure a differential flux (at 10 TeV)::
+
+        >>> pwl = PowerLaw(index=2.3, amplitude=1e-12 * u.Unit('1 / (cm2 s TeV)'), reference=1 * u.TeV)
         >>> crab = CrabSpectrum('hess_pl').model
-        >>>
-        >>> # compute flux at 10 TeV in crab units
         >>> energy = 10 * u.TeV
-        >>> flux_cu = (pwl(energy) / crab(energy)).to('%')
-        >>> print(flux_cu)
+        >>> dnde_cu = (pwl(energy) / crab(energy)).to('%')
+        >>> print(dnde_cu)
         6.19699156377 %
 
-    And the same for integral fluxes:
+    And the same for integral fluxes (between 1 and 10 TeV)::
 
         >>> # compute integral flux in crab units
         >>> emin, emax = [1, 10] * u.TeV
         >>> flux_int_cu = (pwl.integral(emin, emax) / crab.integral(emin, emax)).to('%')
         >>> print(flux_int_cu)
         3.5350582166 %
-
     """
+
+    references = [
+        'meyer', 'hegra', 'hess_pl', 'hess_ecpl',
+    ]
+    """Available references (see class docstring)."""
+
     def __init__(self, reference='meyer'):
+
         if reference == 'meyer':
             model = MeyerCrabModel()
         elif reference == 'hegra':
@@ -137,7 +111,8 @@ def __init__(self, reference='meyer'):
         elif reference == 'hess_ecpl':
             model = ExponentialCutoffPowerLaw(**hess_ecpl)
         else:
-            raise ValueError("Unknown reference, choose one of the following:"
-                             "'meyer', 'hegra', 'hess_pl' or 'hess_ecpl'")
+            fmt = 'Invalid reference: {!r}. Choices: {!r}'
+            raise ValueError(fmt.format(reference, self.references))
+
         self.model = model
         self.reference = reference

gammapy/spectrum/extract.py

@@ -18,7 +18,7 @@
 
 
 class SpectrumExtraction(object):
-    """Creating input data to 1D spectrum fitting
+    """Creating input data to 1D spectrum fitting.
 
     This class is responsible for extracting a
     `~gammapy.spectrum.SpectrumObservation` from a