Add flux points computation methods

examples/example_calculate_fluxpoints.py

@@ -0,0 +1,58 @@
+"""Compute flux points
+
+This is an example script that show how to compute flux points in Gammapy. 
+TODO: Refactor and add to FluxPointsComputation class or so
+"""
+
+from gammapy.spectrum import (
+    SpectrumObservation,
+    SpectrumFit,
+    DifferentialFluxPoints,
+    SpectrumFitResult,
+    SpectrumResult
+)
+import astropy.units as u
+import numpy as np
+import copy
+import matplotlib.pyplot as plt
+
+obs = SpectrumObservation.read('$GAMMAPY_EXTRA/datasets/hess-crab4_pha/pha_obs23523.fits')
+
+fit = SpectrumFit(obs)
+fit.run()
+best_fit = copy.deepcopy(fit.result[0].fit)
+
+# Define Flux points binning
+emin = np.log10(obs.lo_threshold.to('TeV').value)
+emax = np.log10(40)
+binning = np.logspace(emin, emax, 8) * u.TeV
+
+# Fix index
+fit.model.gamma.freeze()
+
+# Fit norm in bands
+diff_flux = list()
+diff_flux_err = list()
+e_err_hi = list()
+e_err_lo = list()
+energy = list()
+for ii in range(len(binning)-1):
+    energ= np.sqrt(binning[ii] * binning[ii+1])
+    energy.append(energ)
+    e_err_hi.append(binning[ii+1] - energ)
+    e_err_lo.append(energ - binning[ii])
+    fit.fit_range = binning[[ii,ii+1]]
+    fit.run()
+    res = fit.result[0].fit
+    diff_flux.append(res.model(energ).to('cm-2 s-1 TeV-1'))
+    err = res.model_with_uncertainties(energ.to('keV').value)
+    diff_flux_err.append(err.s * u.Unit('cm-2 s-1 keV-1'))
+
+points = DifferentialFluxPoints.from_arrays(energy=energy, diff_flux=diff_flux,
+                                            diff_flux_err_hi=diff_flux_err,
+                                            diff_flux_err_lo=diff_flux_err,
+                                            energy_err_hi = e_err_hi,
+                                            energy_err_lo = e_err_lo)
+result = SpectrumResult(fit=best_fit, points=points)
+result.plot_spectrum()
+plt.show()

gammapy/spectrum/flux_point.py

@@ -63,17 +63,21 @@ def from_arrays(cls, energy, diff_flux, energy_err_hi=None,
         # Set errors to zero by default
         def_f = np.zeros(len(energy)) * diff_flux.unit
         def_e = np.zeros(len(energy)) * energy.unit
-        energy_err_hi = def_e if energy_err_hi is None else energy_err_hi
-        energy_err_lo = def_e if energy_err_lo is None else energy_err_lo
-        diff_flux_err_hi = def_f if diff_flux_err_hi is None else diff_flux_err_hi
-        diff_flux_err_lo = def_f if diff_flux_err_lo is None else diff_flux_err_lo
+        if energy_err_hi is None:
+            energy_err_hi = def_e
+        if energy_err_lo is None:
+            energy_err_lo = def_e
+        if diff_flux_err_hi is None:
+            diff_flux_err_hi = def_f
+        if diff_flux_err_lo is None:
+            diff_flux_err_lo = def_f
 
         t['ENERGY'] = energy
-        t['ENERGY_ERR_HI'] = energy_err_hi
-        t['ENERGY_ERR_LO'] = energy_err_lo
+        t['ENERGY_ERR_HI'] = Quantity(energy_err_hi)
+        t['ENERGY_ERR_LO'] = Quantity(energy_err_lo)
         t['DIFF_FLUX'] = diff_flux
-        t['DIFF_FLUX_ERR_HI'] = diff_flux_err_hi
-        t['DIFF_FLUX_ERR_LO'] = diff_flux_err_lo
+        t['DIFF_FLUX_ERR_HI'] = Quantity(diff_flux_err_hi)
+        t['DIFF_FLUX_ERR_LO'] = Quantity(diff_flux_err_lo)
         return cls(t)
 
 

gammapy/spectrum/models.py

@@ -34,58 +34,27 @@ def __str__(self):
             ss += '\n{parname} : {parval:.3g}'.format(**locals())
         return ss
 
-    def to_sherpa(self, name='default'):
-        """Return `~sherpa.models.ArithmeticModel`
-
-        Parameters
-        ----------
-        name : str, optional
-            Name of the sherpa model instance
-        """
-        import sherpa.models as m
-        if isinstance(self, PowerLaw):
-            model = m.PowLaw1D('powlaw1d.' + name)
-            model.gamma = self.parameters.index.value
-        else:
-            raise NotImplementedError
-
-        model.ref = self.parameters.reference.to('keV').value
-        model.ampl = self.parameters.amplitude.to('cm-2 s-1 keV-1').value
-
-        return model
-
-    @classmethod
-    def from_sherpa(cls, model):
-        """Create `~gammapy.spectrum.models.SpectrumModel` from
-        `~sherpa.models.ArithmeticModel`
-
-        Parameters
-        ----------
-        model : `~sherpa.models.ArithmeticModel`
-            Sherpa model
-        """
-        from . import SpectrumFit
-        pardict = dict(gamma = ['index', u.Unit('')],
-                       ref = ['reference', u.keV],
-                       ampl = ['amplitude', SpectrumFit.FLUX_FACTOR * u.Unit('cm-2 s-1 keV-1')])
-        kwargs = dict()
-
-        for par in model.pars:
-            name = par.name
-            kwargs[pardict[name][0]] =  par.val * pardict[name][1]
-
-        return cls(**kwargs)
-
     def to_dict(self):
         """Serialize to dict"""
         retval = dict()
 
         retval['name'] = self.__class__.__name__
+        retval['parameters'] = list()
         for parname, parval in self.parameters.items():
-            retval[parname] = str(parval)
+            retval['parameters'].append(dict(name=parname,
+                                             val=parval.value,
+                                             unit=str(parval.unit)))
         return retval
 
-    def plot(self, ax=None, energy_range=[0.1, 10] * u.TeV,
+    @classmethod
+    def from_dict(cls, val):
+        """Serialize from dict"""
+        kwargs = dict()
+        for _ in val['parameters']:
+            kwargs[_['name']] = _['val'] * u.Unit(_['unit'])
+        return cls(**kwargs)
+
+    def plot(self, energy_range, ax=None, 
              energy_unit='TeV', flux_unit='cm-2 s-1 TeV-1',
              energy_power=0, n_points=100, **kwargs):
         """Plot `~gammapy.spectrum.SpectralModel` 
@@ -102,9 +71,9 @@ def plot(self, ax=None, energy_range=[0.1, 10] * u.TeV,
             Unit of the energy axis
         flux_unit : str, `~astropy.units.Unit`, optional
             Unit of the flux axis
-        energy_power : int
+        energy_power : int, optional
             Power of energy to multiply flux axis with
-        n_points : int
+        n_points : int, optional
             Number of evaluation nodes
 
         Returns
@@ -150,8 +119,8 @@ class PowerLaw(SpectralModel):
     """
     def __init__(self, index, amplitude, reference):
         self.parameters = Bunch(index = index,
-                                amplitude = amplitude.to('cm-2 s-1 TeV-1'),
-                                reference = reference.to('TeV'))
+                                amplitude = amplitude,
+                                reference = reference)
 
     @staticmethod
     def evaluate(energy, index, amplitude, reference):
@@ -168,6 +137,22 @@ def integral(self, emin, emax):
 
         return prefactor * (upper - lower)
 
+    def to_sherpa(self, name='default'):
+        """Return `~sherpa.models.PowLaw1d`
+
+        Parameters
+        ----------
+        name : str, optional
+            Name of the sherpa model instance
+        """
+        import sherpa.models as m
+        model = m.PowLaw1D('powlaw1d.' + name)
+        model.gamma = self.parameters.index.value
+        model.ref = self.parameters.reference.to('keV').value
+        model.ampl = self.parameters.amplitude.to('cm-2 s-1 keV-1').value
+
+        return model
+
 
 class ExponentialCutoffPowerLaw(SpectralModel):
     """Spectral exponential cutoff power-law model.