files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import astropy.units as u
  from ....utils.testing import assert_quantity_allclose, requires_data, requires_dependency
  from ....maps import WcsGeom
  from .. import JFactory, DMFluxMapMaker, profiles, PrimaryFlux
  
  
  @pytest.fixture(scope="session")
  def geom():
      return WcsGeom.create(binsz=0.5, npix=10)
  
  
  @pytest.fixture(scope="session")
  def jfact():
      jfactory = JFactory(geom=geom(), profile=profiles.NFWProfile(), distance=8 * u.kpc)
      return jfactory.compute_jfactor()
  
  
  @pytest.fixture(scope="session")
  def prim_flux():
      return PrimaryFlux(mDM=1 * u.TeV, channel="W")
  
  
  @requires_dependency('scipy')
  @requires_data('gammapy-extra')
  def test_dmfluxmapmaker(jfact, prim_flux):
      x_section = 1e-26 * u.Unit("cm3 s-1")
      energy_range = [0.1, 1] * u.TeV
      dmfmm = DMFluxMapMaker(
          jfact=jfact, prim_flux=prim_flux, x_section=x_section, energy_range=energy_range
      )
      flux = dmfmm.run()
  
      actual = flux[5, 5]
      desired = 6.503327e-13 / u.cm ** 2 / u.s
      assert_quantity_allclose(actual, desired, rtol=1e-5)
  
  
  [... the rest is omitted ...]
files: gammapy/astro/darkmatter/tests/test_utils.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import astropy.units as u
  from ....utils.testing import assert_quantity_allclose, requires_data, requires_dependency
  from ....maps import WcsGeom
  from .. import JFactory, DMFluxMapMaker, profiles, PrimaryFlux
  
  
  @pytest.fixture(scope="session")
  def geom():
      return WcsGeom.create(binsz=0.5, npix=10)
  
  
  @pytest.fixture(scope="session")
  def jfact():
      jfactory = JFactory(ref_geom=geom(), profile=profiles.NFWProfile(), distance=8 * u.kpc)
      return jfactory.compute_jfactor()
  
  
  @pytest.fixture(scope="session")
  def prim_flux():
      return PrimaryFlux(mDM=1 * u.TeV, channel="W")
  
  
  @requires_dependency('scipy')
  @requires_data('gammapy-extra')
  def test_dmfluxmapmaker(jfact, prim_flux):
      x_section = 1e-26 * u.Unit("cm3 s-1")
      energy_range = [0.1, 1] * u.TeV
      dmfmm = DMFluxMapMaker(
          jfact=jfact, prim_flux=prim_flux, x_section=x_section, energy_range=energy_range
      )
      flux = dmfmm.run()
  
      actual = flux[5, 5]
      desired = 6.503327e-13 / u.cm ** 2 / u.s
      assert_quantity_allclose(actual, desired, rtol=1e-5)
  
  
  [... the rest is omitted ...]
files: gammapy/astro/darkmatter/tests/test_utils.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """
  Utilities to compute J-factor maps
  """
  from __future__ import absolute_import, division, print_function, unicode_literals
  import astropy.units as u
  import numpy as np
  
  __all__ = [
      'JFactory',
      'DMFluxMapMaker',
  ]
  
  
  class JFactory(object):
      """Compute J-Factor maps
  
      All J-Factors are computed for annihilation. The assumend dark matter
      profiles will be centered on the center of the map
  
      Parameters
      ----------
      geom : `~gammapy.maps.WcsGeom`
          Reference geometry
      profile : `~gammapy.astro.darkmatter.profiles.DMProfile`
          Dark matter profile
      distance : `~astropy.units.Quantity`
          Distance to convert angular scale of the map
      """
  
      def __init__(self, geom, profile, distance):
          self.geom = geom
          self.profile = profile
          self.distance = distance
  
      def compute_differential_jfactor(self):
          r"""Compute differential J-Factor
  
          .. math ::
              \frac{\mathrm d J}{\mathrm d \Omega} = 
             \int_{\mathrm{LoS}} \mathrm d r \rho(r)
  
          TODO: Needs to be implemented more efficiently
          """
          separation = self.geom.separation(self.geom.center_skydir)
          rmin = separation.rad * self.distance
          rmax = self.distance
          val = [self.profile.integral(_, rmax) for _ in rmin.flatten()]
          jfact = u.Quantity(val).to('GeV2 cm-5').reshape(rmin.shape)
          return jfact / u.steradian
  
      def compute_jfactor(self):
          r"""Compute astrophysical J-Factor
  
          .. math ::
              J(\Delta\Omega) = 
             \int_{\Delta\Omega} \mathrm d \Omega^{\prime}
             \frac{\mathrm d J}{\mathrm d \Omega^{\prime}} 
  
          """
          diff_jfact = self.compute_differential_jfactor()
          jfact = diff_jfact * self.geom.to_image().solid_angle()
          return jfact
  
  
  class DMFluxMapMaker(object):
      r"""Create dark matter flux maps
  
      The gamma-ray flux is computed as follows
  
      .. math::
  
          \frac{\mathrm d \phi}{\mathrm d E \mathrm d\Omega} =
          \frac{\langle \sigma\nu \rangle}{8\pi m^2_{\mathrm{DM}}}
          \frac{\mathrm d N}{\mathrm dE} \times J(\Delta\Omega)
  
      Parameters
      ----------
      jfact : `~astropy.units.Quantity`
          J-Factor as computed by `~gammapy.astro.darkmatter.JFactory`
      prim_flux : `~gammapy.astro.darkmatter.PrimaryFlux`
          Primary gamma-ray flux
      x_section : `~astropy.units.Quantity`
          Velocity averaged annihilation cross section, $\langle \sigma\nu\rangle$
      energy_range : tuple of `~astropy.units.Quantity`
          Energy range for the map
  
      References
      ----------
      * `2011JCAP...03..051 <http://adsabs.harvard.edu/abs/2011JCAP...03..051>`_
      """
  
      def __init__(self, jfact, prim_flux, x_section, energy_range):
          self.jfact = jfact
          self.prim_flux = prim_flux
          self.x_section = x_section
          self.energy_range = energy_range
          self._flux_map = None
  
      def run(self):
          """Compute DM Flux
  
          Returns
          -------
          flux : `~astropy.units.Quantity`
              DM Flux
  
          """
          prefactor = (self.x_section / (8 * np.pi * self.prim_flux.mDM ** 2))
          int_flux = self.prim_flux.table_model.integral(
              emin=self
  
  [... the rest is omitted ...]
files: gammapy/astro/darkmatter/utils.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  """
  Utilities to compute J-factor maps
  """
  from __future__ import absolute_import, division, print_function, unicode_literals
  import astropy.units as u
  import numpy as np
  
  __all__ = [
      'JFactory',
      'DMFluxMapMaker',
  ]
  
  
  class JFactory(object):
      """Compute J-Factor maps
  
      All J-Factors are computed for annihilation. The assumend dark matter
      profiles will be centered on the center of the map
  
      Parameters
      ----------
      ref_geom : `~gammapy.maps.WcsGeom`
          Reference geometry
      profile : `~gammapy.astro.darkmatter.profiles.DMProfile`
          Dark matter profile
      distance : `~astropy.units.Quantity`
          Distance to convert angular scale of the map
      """
  
      def __init__(self, ref_geom, profile, distance):
          self.ref_geom = ref_geom
          self.profile = profile
          self.distance = distance
  
      def compute_differential_jfactor(self):
          r"""Compute differential J-Factor
  
          .. math ::
              \frac{\mathrm d J}{\mathrm d \Omega} = 
             \int_{\mathrm{LoS}} \mathrm d r \rho(r)
  
          TODO: Needs to be implemented more efficiently
          """
          separation = self.ref_geom.separation(self.ref_geom.center_skydir)
          rmin = separation.rad * self.distance
          rmax = self.distance
          val = [self.profile.integral(_, rmax) for _ in rmin.flatten()]
          jfact = u.Quantity(val).to('GeV2 cm-5').reshape(rmin.shape)
          return jfact / u.steradian
  
      def compute_jfactor(self):
          r"""Compute astrophysical J-Factor
  
          .. math ::
              J(\Delta\Omega) = 
             \int_{\Delta\Omega} \mathrm d \Omega^{\prime}
             \frac{\mathrm d J}{\mathrm d \Omega^{\prime}} 
  
          """
          diff_jfact = self.compute_differential_jfactor()
          jfact = diff_jfact * self.ref_geom.to_image().solid_angle()
          return jfact
  
  
  class DMFluxMapMaker(object):
      r"""Create dark matter flux maps
  
      The gamma-ray flux is computed as follows
  
      .. math::
  
          \frac{\mathrm d \phi}{\mathrm d E \mathrm d\Omega} =
          \frac{\langle \sigma\nu \rangle}{8\pi m^2_{\mathrm{DM}}}
          \frac{\mathrm d N}{\mathrm dE} \times J(\Delta\Omega)
  
      Parameters
      ----------
      jfact : `~astropy.units.Quantity`
          J-Factor as computed by `~gammapy.astro.darkmatter.JFactory`
      prim_flux : `~gammapy.astro.darkmatter.PrimaryFlux`
          Primary gamma-ray flux
      x_section : `~astropy.units.Quantity`
          Velocity averaged annihilation cross section, $\langle \sigma\nu\rangle$
      energy_range : tuple of `~astropy.units.Quantity`
          Energy range for the map
  
      References
      ----------
      * `2011JCAP...03..051 <http://adsabs.harvard.edu/abs/2011JCAP...03..051>`_
      """
  
      def __init__(self, jfact, prim_flux, x_section, energy_range):
          self.jfact = jfact
          self.prim_flux = prim_flux
          self.x_section = x_section
          self.energy_range = energy_range
          self._flux_map = None
  
      def run(self):
          """Compute DM Flux
  
          Returns
          -------
          flux : `~astropy.units.Quantity`
              DM Flux
  
          """
          prefactor = (self.x_section / (8 * np.pi * self.prim_flux.mDM ** 2))
          int_flux = self.prim_flux
  
  [... the rest is omitted ...]
files: gammapy/astro/darkmatter/utils.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  from astropy.coordinates import Angle
  from astropy.units import Quantity
  from ..maps import WcsNDMap
  
  __all__ = [
      'make_map_background_irf',
      'make_map_background_fov',
  ]
  
  
  def make_map_background_irf(pointing, livetime, bkg, geom, offset_max, n_integration_bins=1):
      """Compute background map from background IRFs.
  
      TODO: Call a method on bkg that returns integral over energy bin directly
      Related: https://github.com/gammapy/gammapy/pull/1342
  
      Parameters
      ----------
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      livetime : `~astropy.units.Quantity`
          Observation livetime
      bkg : `~gammapy.irf.Background3D`
          Background rate model
      geom : `~gammapy.maps.WcsGeom`
          Reference geometry
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset
      n_integration_bins : int
              Number of bins used to integrate on each energy range
  
      Returns
      -------
      background : `~gammapy.maps.WcsNDMap`
          Background predicted counts sky cube in reco energy
      """
      # Compute the expected background
      # TODO: properly transform FOV to sky coordinates
      # For now we assume the background is radially symmetric
  
      energy_axis = geom.axes[0]
      # Compute offsets of all pixels
      map_coord = geom.get_coord()
      # Retrieve energies from map coordinates
      energy_reco = map_coord[energy_axis.name] * energy_axis.unit
      # TODO: go from SkyCoord to FOV coordinates. Here assume symmetric geometry for fov_lon, fov_lat
      # Compute offset at all the pixels and energy of the Map
      fov_lon = map_coord.skycoord.separation(pointing)
      fov_lat = Angle(np.zeros_like(fov_lon), fov_lon.unit)
      data_int = Quantity(np.zeros_like(fov_lat.value), "s^-1 sr^-1")
      for ie, (e_lo, e_hi) in enumerate(zip(energy_axis.edges[0:-1], energy_axis.edges[1:])):
          data_int[ie, :, :] = bkg.integrate_on_energy_range(
              fov_lon=fov_lon[0, :, :],
              fov_lat=fov_lat[0, :, :],
              energy_range=[e_lo * energy_axis.unit, e_hi * energy_axis.unit],
              n_integration_bins=n_integration_bins,
          )
  
      d_omega = geom.solid_angle()
      data = (data_int * d_omega * livetime).to('').value
  
      # Put exposure outside offset max to zero
      # This might be more generaly dealt with a mask map
      offset = np.sqrt(fov_lon ** 2 + fov_lat ** 2)
      data[:, offset[0, :, :] >= offset_max] = 0
  
      return WcsNDMap(geom, data=data)
  
  
  def make_map_background_fov(acceptance_map, counts_map, exclusion_mask):
      """Build Normalized background map from a given acceptance map and counts map.
  
      This operation is normally performed on single observation maps.
      An exclusion map is used to avoid using regions with significant gamma-ray emission.
      All maps are assumed to follow the same WcsGeom.
  
      TODO: A model map could be used instead of an exclusion mask.
  
      Parameters
      ----------
      acceptance_map : `~gammapy.maps.WcsNDMap`
          Observation hadron acceptance map (i.e. predicted
  
  [... the rest is omitted ...]
files: gammapy/cube/background.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  from astropy.coordinates import Angle
  from astropy.units import Quantity
  from ..maps import WcsNDMap
  
  __all__ = [
      'make_map_background_irf',
      'make_map_background_fov',
  ]
  
  
  def make_map_background_irf(pointing, livetime, bkg, ref_geom, offset_max, n_integration_bins=1):
      """Compute background map from background IRFs.
  
      TODO: Call a method on bkg that returns integral over energy bin directly
      Related: https://github.com/gammapy/gammapy/pull/1342
  
      Parameters
      ----------
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      livetime : `~astropy.units.Quantity`
          Observation livetime
      bkg : `~gammapy.irf.Background3D`
          Background rate model
      ref_geom : `~gammapy.maps.WcsGeom`
          Reference geometry
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset
      n_integration_bins : int
              Number of bins used to integrate on each energy range
  
      Returns
      -------
      background : `~gammapy.maps.WcsNDMap`
          Background predicted counts sky cube in reco energy
      """
      # Compute the expected background
      # TODO: properly transform FOV to sky coordinates
      # For now we assume the background is radially symmetric
  
      energy_axis = ref_geom.axes[0]
      # Compute offsets of all pixels
      map_coord = ref_geom.get_coord()
      # Retrieve energies from map coordinates
      energy_reco = map_coord[energy_axis.name] * energy_axis.unit
      # TODO: go from SkyCoord to FOV coordinates. Here assume symmetric geometry for fov_lon, fov_lat
      # Compute offset at all the pixels and energy of the Map
      fov_lon = map_coord.skycoord.separation(pointing)
      fov_lat = Angle(np.zeros_like(fov_lon), fov_lon.unit)
      data_int = Quantity(np.zeros_like(fov_lat.value), "s^-1 sr^-1")
      for ie, (e_lo, e_hi) in enumerate(zip(energy_axis.edges[0:-1], energy_axis.edges[1:])):
          data_int[ie, :, :] = bkg.integrate_on_energy_range(
              fov_lon=fov_lon[0, :, :],
              fov_lat=fov_lat[0, :, :],
              energy_range=[e_lo * energy_axis.unit, e_hi * energy_axis.unit],
              n_integration_bins=n_integration_bins,
          )
  
      d_omega = ref_geom.solid_angle()
      data = (data_int * d_omega * livetime).to('').value
  
      # Put exposure outside offset max to zero
      # This might be more generaly dealt with a mask map
      offset = np.sqrt(fov_lon ** 2 + fov_lat ** 2)
      data[:, offset[0, :, :] >= offset_max] = 0
  
      return WcsNDMap(ref_geom, data=data)
  
  
  def make_map_background_fov(acceptance_map, counts_map, exclusion_mask):
      """Build Normalized background map from a given acceptance map and counts map.
  
      This operation is normally performed on single observation maps.
      An exclusion map is used to avoid using regions with significant gamma-ray emission.
      All maps are assumed to follow the same WcsGeom.
  
      TODO: A model map could be used instead of an exclusion mask.
  
      Parameters
      ----------
      acceptance_map : `~gammapy.maps.WcsNDMap`
  
  
  [... the rest is omitted ...]
files: gammapy/cube/background.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  from ..maps import WcsNDMap
  
  __all__ = [
      'fill_map_counts',
      'make_map_counts',
  ]
  
  
  def fill_map_counts(count_map, event_list):
      """Fill events into a counts map.
  
      The energy of the events is used for a non-spatial axis homogeneous to energy.
      The other non-spatial axis names should have an entry in the column names of the ``EventList``
  
      Parameters
      ----------
      count_map : `~gammapy.maps.Map`
          Map object, will be filled by this function.
      event_list : `~gammapy.data.EventList`
          Event list
      """
      geom = count_map.geom
  
      # Make a coordinate dictionary; skycoord is always added
      coord_dict = dict(skycoord=event_list.radec)
  
      # Now add one coordinate for each extra map axis
      for axis in geom.axes:
          if axis.type == 'energy':
              # This axis is the energy. We treat it differently because axis.name could be e.g. 'energy_reco'
              coord_dict[axis.name] = event_list.energy.to(axis.unit)
          # TODO: add proper extraction for time
          else:
              # We look for other axes name in the table column names (case insensitive)
              colnames = [_.upper() for _ in event_list.table.colnames]
              if axis.name.upper() in colnames:
                  column_name = event_list.table.colnames[colnames.index(axis.name.upper())]
                  coord_dict.update({axis.name: event_list.table[column_name].to(axis.unit)})
              else:
                  raise ValueError("Cannot find MapGeom axis {!r} in EventList".format(axis.name))
  
      count_map.fill_by_coord(coord_dict)
  
  
  def make_map_counts(events, geom, pointing, offset_max):
      """Build a WcsNDMap (space - energy) with events from an EventList.
  
      The energy of the events is used for the non-spatial axis.
  
      Parameters
      ----------
      events : `~gammapy.data.EventList`
          Event list
      geom : `~gammapy.maps.WcsGeom`
          Reference WcsGeom object used to define geometry (space - energy)
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset.
  
      Returns
      -------
      cntmap : `~gammapy.maps.WcsNDMap`
          Count cube (3D) in true energy bins
      """
      counts_map = WcsNDMap(geom)
      fill_map_counts(counts_map, events)
  
      # Compute and apply FOV offset mask
      offset = geom.separation(pointing)
      offset_mask = offset >= offset_max
      counts_map.data[:, offset_mask] = 0
  
      return counts_map
  
  
  [... the rest is omitted ...]
files: gammapy/cube/counts.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  from ..maps import WcsNDMap
  
  __all__ = [
      'fill_map_counts',
      'make_map_counts',
  ]
  
  
  def fill_map_counts(count_map, event_list):
      """Fill events into a counts map.
  
      The energy of the events is used for a non-spatial axis homogeneous to energy.
      The other non-spatial axis names should have an entry in the column names of the ``EventList``
  
      Parameters
      ----------
      count_map : `~gammapy.maps.Map`
          Map object, will be filled by this function.
      event_list : `~gammapy.data.EventList`
          Event list
      """
      geom = count_map.geom
  
      # Make a coordinate dictionary; skycoord is always added
      coord_dict = dict(skycoord=event_list.radec)
  
      # Now add one coordinate for each extra map axis
      for axis in geom.axes:
          if axis.type == 'energy':
              # This axis is the energy. We treat it differently because axis.name could be e.g. 'energy_reco'
              coord_dict[axis.name] = event_list.energy.to(axis.unit)
          # TODO: add proper extraction for time
          else:
              # We look for other axes name in the table column names (case insensitive)
              colnames = [_.upper() for _ in event_list.table.colnames]
              if axis.name.upper() in colnames:
                  column_name = event_list.table.colnames[colnames.index(axis.name.upper())]
                  coord_dict.update({axis.name: event_list.table[column_name].to(axis.unit)})
              else:
                  raise ValueError("Cannot find MapGeom axis {!r} in EventList".format(axis.name))
  
      count_map.fill_by_coord(coord_dict)
  
  
  def make_map_counts(events, ref_geom, pointing, offset_max):
      """Build a WcsNDMap (space - energy) with events from an EventList.
  
      The energy of the events is used for the non-spatial axis.
  
      Parameters
      ----------
      events : `~gammapy.data.EventList`
          Event list
      ref_geom : `~gammapy.maps.WcsGeom`
          Reference WcsGeom object used to define geometry (space - energy)
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset.
  
      Returns
      -------
      cntmap : `~gammapy.maps.WcsNDMap`
          Count cube (3D) in true energy bins
      """
      counts_map = WcsNDMap(ref_geom)
      fill_map_counts(counts_map, events)
  
      # Compute and apply FOV offset mask
      offset = ref_geom.separation(pointing)
      offset_mask = offset >= offset_max
      counts_map.data[:, offset_mask] = 0
  
      return counts_map
  
  
  [... the rest is omitted ...]
files: gammapy/cube/counts.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  from ..maps import WcsNDMap
  
  __all__ = [
      'make_map_exposure_true_energy',
  ]
  
  
  def make_map_exposure_true_energy(pointing, livetime, aeff, geom, offset_max):
      """Compute exposure WcsNDMap in true energy (i.e. not convolved by Edisp).
  
      Parameters
      ----------
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      livetime : `~astropy.units.Quantity`
          Livetime
      aeff : `~gammapy.irf.EffectiveAreaTable2D`
          Effective area table
      geom : `~gammapy.maps.WcsGeom`
          Reference WcsGeom object used to define geometry (space - energy)
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset.
  
      Returns
      -------
      expmap : `~gammapy.maps.WcsNDMap`
          Exposure cube (3D) in true energy bins
      """
      offset = geom.separation(pointing)
  
      # Retrieve energies from WcsNDMap
      # Note this would require a log_center from the geometry
      # Or even better edges, but WcsNDmap does not really allows it.
      energy = geom.axes[0].center * geom.axes[0].unit
  
      exposure = aeff.data.evaluate(offset=offset, energy=energy)
      exposure *= livetime
  
      # We check if exposure is a 3D array in case there is a single bin in energy
      # TODO: call np.atleast_3d ?
      if len(exposure.shape) < 3:
          exposure = np.expand_dims(exposure.value, 0) * exposure.unit
  
      # Put exposure outside offset max to zero
      # This might be more generaly dealt with a mask map
      exposure[:, offset >= offset_max] = 0
  
      data = exposure.to('m2 s')
  
      return WcsNDMap(geom, data)
  
  
  [... the rest is omitted ...]
files: gammapy/cube/exposure.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  from ..maps import WcsNDMap
  
  __all__ = [
      'make_map_exposure_true_energy',
  ]
  
  
  def make_map_exposure_true_energy(pointing, livetime, aeff, ref_geom, offset_max):
      """Compute exposure WcsNDMap in true energy (i.e. not convolved by Edisp).
  
      Parameters
      ----------
      pointing : `~astropy.coordinates.SkyCoord`
          Pointing direction
      livetime : `~astropy.units.Quantity`
          Livetime
      aeff : `~gammapy.irf.EffectiveAreaTable2D`
          Effective area table
      ref_geom : `~gammapy.maps.WcsGeom`
          Reference WcsGeom object used to define geometry (space - energy)
      offset_max : `~astropy.coordinates.Angle`
          Maximum field of view offset.
  
      Returns
      -------
      expmap : `~gammapy.maps.WcsNDMap`
          Exposure cube (3D) in true energy bins
      """
      offset = ref_geom.separation(pointing)
  
      # Retrieve energies from WcsNDMap
      # Note this would require a log_center from the geometry
      # Or even better edges, but WcsNDmap does not really allows it.
      energy = ref_geom.axes[0].center * ref_geom.axes[0].unit
  
      exposure = aeff.data.evaluate(offset=offset, energy=energy)
      exposure *= livetime
  
      # We check if exposure is a 3D array in case there is a single bin in energy
      # TODO: call np.atleast_3d ?
      if len(exposure.shape) < 3:
          exposure = np.expand_dims(exposure.value, 0) * exposure.unit
  
      # Put exposure outside offset max to zero
      # This might be more generaly dealt with a mask map
      exposure[:, offset >= offset_max] = 0
  
      data = exposure.to('m2 s')
  
      return WcsNDMap(ref_geom, data)
  
  
  [... the rest is omitted ...]
files: gammapy/cube/exposure.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import logging
  from astropy.utils.console import ProgressBar
  from astropy.nddata.utils import PartialOverlapError
  from astropy.coordinates import Angle
  from ..maps import WcsNDMap
  from .counts import make_map_counts
  from .exposure import make_map_exposure_true_energy
  from .background import make_map_background_irf, make_map_background_fov
  
  __all__ = [
      'MapMaker',
  ]
  
  log = logging.getLogger(__name__)
  
  
  class MapMaker(object):
      """Make all basic maps from observations.
  
      Parameters
      ----------
      geom : `~gammapy.maps.WcsGeom`
          Reference image geometry
      offset_max : `~astropy.coordinates.Angle`
          Maximum offset angle
      cutout_mode : {'trim', 'strict'}, optional
          Options for making cutouts, see :func: `~gammapy.maps.WcsNDMap.make_cutout`
          Should be left to the default value 'trim'
          unless you want only fully contained observations to be added to the map
      """
  
      def __init__(self, geom, offset_max, cutout_mode="trim"):
          self.offset_max = Angle(offset_max)
          self.geom = geom
  
          # We instantiate the end products of the MakeMaps class
          self.counts_map = WcsNDMap(self.geom)
  
          self.exposure_map = WcsNDMap(self.geom, unit="m2 s")
  
          self.background_map = WcsNDMap(self.geom)
  
          # We will need this general exclusion mask for the analysis
          self.exclusion_map = WcsNDMap(self.geom)
          self.exclusion_map.data += 1
  
          self.cutout_mode = cutout_mode
          self.maps = {}
  
      def process_obs(self, obs):
          """Process one observation and add it to the cutout image
  
          Parameters
          ----------
          obs : `~gammapy.data.DataStoreObservation`
              Observation
          """
          # First make cutout of the global image
          try:
              exclusion_mask_cutout, cutout_slices = self.exclusion_map.make_cutout(
                  obs.pointing_radec, 2 * self.offset_max, mode=self.cutout_mode
              )
          except PartialOverlapError:
              # TODO: can we silently do the right thing here? Discuss
              log.info("Observation {} not fully contained in target image. Skipping it.".format(obs.obs_id))
              return
  
          cutout_geom = exclusion_mask_cutout.geom
  
          counts_obs_map = make_map_counts(
              obs.events, cutout_geom, obs.pointing_radec, self.offset_max,
          )
  
          expo_obs_map = make_map_exposure_true_energy(
              obs.pointing_radec, obs.observation_live_time_duration,
              obs.aeff, cutout_geom, self.offset_max,
          )
  
          acceptance_obs_map = make_map_background_irf(
              obs.pointing_radec, obs.observation_live_time_duration,
              obs.bkg, cutout_geom, self.offset_max,
          )
  
          background_obs_map = make_map_background_fov(
              acceptance_obs_map, counts_obs_map, exclusion_mask_cutout,
          )
  
          self._add_cutouts(cutout_slices, counts_obs_map, expo_obs_map, background_obs_map)
  
      def _add_cutouts(self, cutout_slices, counts_obs_map, expo_obs_map, acceptance_obs_map):
          """Add current cutout to global maps."""
          self.counts_map.data[cutout_slices] += counts_obs_map.data
          self.exposure_map.data[cutout_slices] += expo_obs_map.quantity.to(self.exposure_map.unit).value
          self.background_map.data[cutout_slices]
  
  [... the rest is omitted ...]
files: gammapy/cube/make.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import logging
  from astropy.utils.console import ProgressBar
  from astropy.nddata.utils import PartialOverlapError
  from astropy.coordinates import Angle
  from ..maps import WcsNDMap
  from .counts import make_map_counts
  from .exposure import make_map_exposure_true_energy
  from .background import make_map_background_irf, make_map_background_fov
  
  __all__ = [
      'MapMaker',
  ]
  
  log = logging.getLogger(__name__)
  
  
  class MapMaker(object):
      """Make all basic maps from observations.
  
      Parameters
      ----------
      ref_geom : `~gammapy.maps.WcsGeom`
          Reference image geometry
      offset_max : `~astropy.coordinates.Angle`
          Maximum offset angle
      cutout_mode : {'trim', 'strict'}, optional
          Options for making cutouts, see :func: `~gammapy.maps.WcsNDMap.make_cutout`
          Should be left to the default value 'trim'
          unless you want only fully contained observations to be added to the map
      """
  
      def __init__(self, ref_geom, offset_max, cutout_mode="trim"):
          self.offset_max = Angle(offset_max)
          self.ref_geom = ref_geom
  
          # We instantiate the end products of the MakeMaps class
          self.counts_map = WcsNDMap(self.ref_geom)
  
          self.exposure_map = WcsNDMap(self.ref_geom, unit="m2 s")
  
          self.background_map = WcsNDMap(self.ref_geom)
  
          # We will need this general exclusion mask for the analysis
          self.exclusion_map = WcsNDMap(self.ref_geom)
          self.exclusion_map.data += 1
  
          self.cutout_mode = cutout_mode
          self.maps = {}
  
      def process_obs(self, obs):
          """Process one observation and add it to the cutout image
  
          Parameters
          ----------
          obs : `~gammapy.data.DataStoreObservation`
              Observation
          """
          # First make cutout of the global image
          try:
              exclusion_mask_cutout, cutout_slices = self.exclusion_map.make_cutout(
                  obs.pointing_radec, 2 * self.offset_max, mode=self.cutout_mode
              )
          except PartialOverlapError:
              # TODO: can we silently do the right thing here? Discuss
              log.info("Observation {} not fully contained in target image. Skipping it.".format(obs.obs_id))
              return
  
          cutout_geom = exclusion_mask_cutout.geom
  
          counts_obs_map = make_map_counts(
              obs.events, cutout_geom, obs.pointing_radec, self.offset_max,
          )
  
          expo_obs_map = make_map_exposure_true_energy(
              obs.pointing_radec, obs.observation_live_time_duration,
              obs.aeff, cutout_geom, self.offset_max,
          )
  
          acceptance_obs_map = make_map_background_irf(
              obs.pointing_radec, obs.observation_live_time_duration,
              obs.bkg, cutout_geom, self.offset_max,
          )
  
          background_obs_map = make_map_background_fov(
              acceptance_obs_map, counts_obs_map, exclusion_mask_cutout,
          )
  
          self._add_cutouts(cutout_slices, counts_obs_map, expo_obs_map, background_obs_map)
  
      def _add_cutouts(self, cutout_slices, counts_obs_map, expo_obs_map, acceptance_obs_map):
          """Add current cutout to global maps."""
          self.counts_map.data[cutout_slices] += counts_obs_map.data
          self.exposure_map.data[cutout_slices] += expo_obs_map.quantity.to(self.exposure_map.unit).value
         
  
  [... the rest is omitted ...]
files: gammapy/cube/make.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  import astropy.units as u
  from astropy.coordinates import Angle
  from ..irf import EnergyDependentTablePSF
  from ..maps import Map
  from ..cube import PSFKernel
  
  __all__ = [
      'make_psf_map',
      'PSFMap'
  ]
  
  
  def make_psf_map(psf, pointing, geom, max_offset):
      """Make a psf map for a single observation
  
      Expected axes : rad and true energy in this specific order
      The name of the rad MapAxis is expected to be 'rad'
  
      Parameters
      ----------
      psf : `~gammapy.irf.PSF3D`
          the PSF IRF
      pointing : `~astropy.coordinates.SkyCoord`
          the pointing direction
      geom : `~gammapy.maps.MapGeom`
          the map geom to be used. It provides the target geometry.
          rad and true energy axes should be given in this specific order.
      max_offset : `~astropy.coordinates.Angle`
          maximum offset w.r.t. fov center
  
      Returns
      -------
      psfmap : `~gammapy.cube.PSFMap`
          the resulting PSF map
      """
      energy_axis = geom.get_axis_by_name('energy_true')
      energy = energy_axis.center * energy_axis.unit
  
      rad_axis = geom.get_axis_by_name('theta')
      rad = Angle(rad_axis.center, unit=rad_axis.unit)
  
      # Compute separations with pointing position
      separations = pointing.separation(geom.to_image().get_coord().skycoord)
      valid = np.where(separations < max_offset)
  
      # Compute PSF values
      psf_values = psf.evaluate(offset=separations[valid], energy=energy, rad=rad)
  
      # Re-order axes to be consistent with expected geometry
      psf_values = np.transpose(psf_values, axes=(2, 0, 1))
  
      # Create Map and fill relevant entries
      psfmap = Map.from_geom(geom, unit='sr-1')
      psfmap.data[:, :, valid[0], valid[1]] += psf_values.to(psfmap.unit).value
  
      return PSFMap(psfmap)
  
  
  class PSFMap(object):
      """Class containing the Map of PSFs and allowing to interact with it.
  
      Parameters
      ----------
      psf_map : `~gammapy.maps.Map`
          the input PSF Map. Should be a Map with 2 non spatial axes.
          rad and true energy axes should be given in this specific order.
  
      Examples
      --------
  
      .. code:: python
  
          import numpy as np
          from gammapy.maps import Map, WcsGeom, MapAxis
          from gammapy.irf import EnergyDependentMultiGaussPSF
          from gammapy.cube import make_psf_map, PSFMap
          from astropy import units as u
          from astropy.coordinates import SkyCoord
  
          # Define energy axis. Note that the name is fixed.
          energy_axis = MapAxis.from_edges(np.logspace(-1., 1., 4), unit='TeV', name='energy_true')
          # Define rad axis. Again note the axis name
          rads = np.linspace(0., 0.5, 100) * u.deg
          rad_axis = MapAxis.from_edges(rads, unit='deg', name='theta')
  
          # Define parameters
          pointing = SkyCoord(0., 0., unit='deg')
          max_offset = 4 * u.deg
  
          # Create WcsGeom
          geom = WcsGeom.create(binsz=0.25*u.deg, width=10*u.deg, skydir=point
  
  [... the rest is omitted ...]
files: gammapy/cube/psf_map.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import numpy as np
  import astropy.units as u
  from astropy.coordinates import Angle
  from ..irf import EnergyDependentTablePSF
  from ..maps import Map
  from ..cube import PSFKernel
  
  __all__ = [
      'make_psf_map',
      'PSFMap'
  ]
  
  
  def make_psf_map(psf, pointing, ref_geom, max_offset):
      """Make a psf map for a single observation
  
      Expected axes : rad and true energy in this specific order
      The name of the rad MapAxis is expected to be 'rad'
  
      Parameters
      ----------
      psf : `~gammapy.irf.PSF3D`
          the PSF IRF
      pointing : `~astropy.coordinates.SkyCoord`
          the pointing direction
      ref_geom : `~gammapy.maps.MapGeom`
          the map geom to be used. It provides the target geometry.
          rad and true energy axes should be given in this specific order.
      max_offset : `~astropy.coordinates.Angle`
          maximum offset w.r.t. fov center
  
      Returns
      -------
      psfmap : `~gammapy.cube.PSFMap`
          the resulting PSF map
      """
      energy_axis = ref_geom.get_axis_by_name('energy_true')
      energy = energy_axis.center * energy_axis.unit
  
      rad_axis = ref_geom.get_axis_by_name('theta')
      rad = Angle(rad_axis.center, unit=rad_axis.unit)
  
      # Compute separations with pointing position
      separations = pointing.separation(ref_geom.to_image().get_coord().skycoord)
      valid = np.where(separations < max_offset)
  
      # Compute PSF values
      psf_values = psf.evaluate(offset=separations[valid], energy=energy, rad=rad)
  
      # Re-order axes to be consistent with expected geometry
      psf_values = np.transpose(psf_values, axes=(2, 0, 1))
  
      # Create Map and fill relevant entries
      psfmap = Map.from_geom(ref_geom, unit='sr-1')
      psfmap.data[:, :, valid[0], valid[1]] += psf_values.to(psfmap.unit).value
  
      return PSFMap(psfmap)
  
  
  class PSFMap(object):
      """Class containing the Map of PSFs and allowing to interact with it.
  
      Parameters
      ----------
      psf_map : `~gammapy.maps.Map`
          the input PSF Map. Should be a Map with 2 non spatial axes.
          rad and true energy axes should be given in this specific order.
  
      Examples
      --------
  
      .. code:: python
  
          import numpy as np
          from gammapy.maps import Map, WcsGeom, MapAxis
          from gammapy.irf import EnergyDependentMultiGaussPSF
          from gammapy.cube import make_psf_map, PSFMap
          from astropy import units as u
          from astropy.coordinates import SkyCoord
  
          # Define energy axis. Note that the name is fixed.
          energy_axis = MapAxis.from_edges(np.logspace(-1., 1., 4), unit='TeV', name='energy_true')
          # Define rad axis. Again note the axis name
          rads = np.linspace(0., 0.5, 100) * u.deg
          rad_axis = MapAxis.from_edges(rads, unit='deg', name='theta')
  
          # Define parameters
          pointing = SkyCoord(0., 0., unit='deg')
          max_offset = 4 * u.deg
  
          # Create WcsGeom
          geom = WcsGeom.create(binsz=0.25*u.deg, width=
  
  [... the rest is omitted ...]
files: gammapy/cube/psf_map.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import numpy as np
  from astropy.coordinates import SkyCoord, Angle
  from ...utils.testing import requires_data, assert_quantity_allclose
  from ...maps import Map
  from ...irf import EffectiveAreaTable2D
  from ..exposure import make_map_exposure_true_energy
  
  pytest.importorskip('scipy')
  
  
  @pytest.fixture(scope='session')
  def aeff():
      filename = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_aeff_2d_023523.fits.gz'
      return EffectiveAreaTable2D.read(filename, hdu='AEFF_2D')
  
  
  @pytest.fixture(scope='session')
  def counts_cube():
      filename = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/hess_events_simulated_023523_cntcube.fits'
      return Map.read(filename)
  
  
  @requires_data('gammapy-extra')
  def test_make_map_exposure_true_energy(aeff, counts_cube):
  
      m = make_map_exposure_true_energy(
          pointing=SkyCoord(83.633, 21.514, unit='deg'),
          livetime='1581.17 s',
          aeff=aeff,
          geom=counts_cube.geom,
          offset_max=Angle('2.2 deg'),
      )
  
      assert m.data.shape == (15, 120, 200)
      assert m.unit == 'm2 s'
      assert_quantity_allclose(np.nanmax(m.data), 4.7e8, rtol=100)
  
  
  [... the rest is omitted ...]
files: gammapy/cube/tests/test_exposure.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import numpy as np
  from astropy.coordinates import SkyCoord, Angle
  from ...utils.testing import requires_data, assert_quantity_allclose
  from ...maps import Map
  from ...irf import EffectiveAreaTable2D
  from ..exposure import make_map_exposure_true_energy
  
  pytest.importorskip('scipy')
  
  
  @pytest.fixture(scope='session')
  def aeff():
      filename = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/run023400-023599/run023523/hess_aeff_2d_023523.fits.gz'
      return EffectiveAreaTable2D.read(filename, hdu='AEFF_2D')
  
  
  @pytest.fixture(scope='session')
  def counts_cube():
      filename = '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/hess_events_simulated_023523_cntcube.fits'
      return Map.read(filename)
  
  
  @requires_data('gammapy-extra')
  def test_make_map_exposure_true_energy(aeff, counts_cube):
  
      m = make_map_exposure_true_energy(
          pointing=SkyCoord(83.633, 21.514, unit='deg'),
          livetime='1581.17 s',
          aeff=aeff,
          ref_geom=counts_cube.geom,
          offset_max=Angle('2.2 deg'),
      )
  
      assert m.data.shape == (15, 120, 200)
      assert m.unit == 'm2 s'
      assert_quantity_allclose(np.nanmax(m.data), 4.7e8, rtol=100)
  
  
  [... the rest is omitted ...]
files: gammapy/cube/tests/test_exposure.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.coordinates import SkyCoord
  from ...utils.testing import requires_data, requires_dependency
  from ...irf import EffectiveAreaTable2D, EnergyDependentMultiGaussPSF
  from ...irf.energy_dispersion import EnergyDispersion
  from ...maps import MapAxis, WcsGeom, WcsNDMap, Map
  from ...image.models import SkyGaussian
  from ...spectrum.models import PowerLaw
  from .. import (
      SkyModel,
      MapEvaluator,
      MapFit,
      make_map_exposure_true_energy,
      PSFKernel,
  )
  
  
  @pytest.fixture(scope='session')
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1., 1., 3), name="energy", unit=u.TeV)
      return WcsGeom.create(skydir=(0, 0), binsz=0.02, width=(2, 2),
                            coordsys='GAL', axes=[axis])
  
  
  @pytest.fixture(scope='session')
  def exposure(geom):
      filename = '$GAMMAPY_EXTRA/datasets/cta-1dc/caldb/data/cta//1dc/bcf/South_z20_50h/irf_file.fits'
      aeff = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA')
  
      offset_max = 3 * u.deg
  
      exposure_map = make_map_exposure_true_energy(
          pointing=SkyCoord(1, 0.5, unit='deg', frame='galactic'),
          livetime='1 hour',
          aeff=aeff,
          geom=geom,
          offset_max=offset_max,
      )
      return exposure_map
  
  
  @pytest.fixture(scope='session')
  def background(geom):
      m = Map.from_geom(geom)
      m.quantity = np.ones(m.data.shape) * 1e-5
      return m
  
  
  @pytest.fixture(scope='session')
  def edisp(geom):
      e_true = geom.get_axis_by_name('energy').edges
      return EnergyDispersion.from_diagonal_matrix(e_true=e_true)
  
  
  @pytest.fixture(scope='session')
  def psf(geom):
      filename = '$GAMMAPY_EXTRA/datasets/cta-1dc/caldb/data/cta//1dc/bcf/South_z20_50h/irf_file.fits'
      psf = EnergyDependentMultiGaussPSF.read(filename, hdu='POINT SPREAD FUNCTION')
  
      table_psf = psf.to_energy_dependent_table_psf(theta=0.5 * u.deg)
      psf_kernel = PSFKernel.from_table_psf(table_psf,
                                            geom,
                                            max_radius=0.5 * u.deg)
      return psf_kernel
  
  
  @pytest.fixture
  def sky_model():
      spatial_model = SkyGaussian(
          lon_0='0.2 deg',
          lat_0='0.1 deg',
          sigma='0.2 deg',
      )
      spectral_model = PowerLaw(
          index=3,
          amplitude='1e-11 cm-2 s-1 TeV-1',
          reference='1 TeV',
      )
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
      )
  
  
  @pytest.fixture
  def counts(sky_model, exposure, background, psf, edisp):
      evaluator = MapEvaluator(
          sky_model=sky_model,
          exposure=exposure,
          background=background,
          psf=psf,
          edisp=edisp,
      )
      npred = evaluator.compute_npred()
      return WcsNDMap(exposure.geom, npred)
  
  
  @requires_dependency('scipy')
  @
  
  [... the rest is omitted ...]
files: gammapy/cube/tests/test_fit.py
files: # Licensed under a 3-clause BSD style license - see LICENSE.rst
  from __future__ import absolute_import, division, print_function, unicode_literals
  import pytest
  import numpy as np
  from numpy.testing import assert_allclose
  import astropy.units as u
  from astropy.coordinates import SkyCoord
  from ...utils.testing import requires_data, requires_dependency
  from ...irf import EffectiveAreaTable2D, EnergyDependentMultiGaussPSF
  from ...irf.energy_dispersion import EnergyDispersion
  from ...maps import MapAxis, WcsGeom, WcsNDMap, Map
  from ...image.models import SkyGaussian
  from ...spectrum.models import PowerLaw
  from .. import (
      SkyModel,
      MapEvaluator,
      MapFit,
      make_map_exposure_true_energy,
      PSFKernel,
  )
  
  
  @pytest.fixture(scope='session')
  def geom():
      axis = MapAxis.from_edges(np.logspace(-1., 1., 3), name="energy", unit=u.TeV)
      return WcsGeom.create(skydir=(0, 0), binsz=0.02, width=(2, 2),
                            coordsys='GAL', axes=[axis])
  
  
  @pytest.fixture(scope='session')
  def exposure(geom):
      filename = '$GAMMAPY_EXTRA/datasets/cta-1dc/caldb/data/cta//1dc/bcf/South_z20_50h/irf_file.fits'
      aeff = EffectiveAreaTable2D.read(filename, hdu='EFFECTIVE AREA')
  
      offset_max = 3 * u.deg
  
      exposure_map = make_map_exposure_true_energy(
          pointing=SkyCoord(1, 0.5, unit='deg', frame='galactic'),
          livetime='1 hour',
          aeff=aeff,
          ref_geom=geom,
          offset_max=offset_max,
      )
      return exposure_map
  
  
  @pytest.fixture(scope='session')
  def background(geom):
      m = Map.from_geom(geom)
      m.quantity = np.ones(m.data.shape) * 1e-5
      return m
  
  
  @pytest.fixture(scope='session')
  def edisp(geom):
      e_true = geom.get_axis_by_name('energy').edges
      return EnergyDispersion.from_diagonal_matrix(e_true=e_true)
  
  
  @pytest.fixture(scope='session')
  def psf(geom):
      filename = '$GAMMAPY_EXTRA/datasets/cta-1dc/caldb/data/cta//1dc/bcf/South_z20_50h/irf_file.fits'
      psf = EnergyDependentMultiGaussPSF.read(filename, hdu='POINT SPREAD FUNCTION')
  
      table_psf = psf.to_energy_dependent_table_psf(theta=0.5 * u.deg)
      psf_kernel = PSFKernel.from_table_psf(table_psf,
                                            geom,
                                            max_radius=0.5 * u.deg)
      return psf_kernel
  
  
  @pytest.fixture
  def sky_model():
      spatial_model = SkyGaussian(
          lon_0='0.2 deg',
          lat_0='0.1 deg',
          sigma='0.2 deg',
      )
      spectral_model = PowerLaw(
          index=3,
          amplitude='1e-11 cm-2 s-1 TeV-1',
          reference='1 TeV',
      )
      return SkyModel(
          spatial_model=spatial_model,
          spectral_model=spectral_model,
      )
  
  
  @pytest.fixture
  def counts(sky_model, exposure, background, psf, edisp):
      evaluator = MapEvaluator(
          sky_model=sky_model,
          exposure=exposure,
          background=background,
          psf=psf,
          edisp=edisp,
      )
      npred = evaluator.compute_npred()
      return WcsNDMap(exposure.geom, npred)
  
  
  @requires_dependency('scipy
  
  [... the rest is omitted ...]
files: gammapy/cube/tests/test_fit.py