Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add make_images method in MapMaker #1681

Merged
merged 7 commits into from Aug 12, 2018
Copy path View file
@@ -1,10 +1,13 @@
# 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
from ..irf import EnergyDispersion
from ..spectrum.models import PowerLaw
from ..maps import WcsNDMap, MapAxis, Map

__all__ = [
'make_map_exposure_true_energy',
'make_map_exposure_reco_energy',
]


@@ -42,3 +45,59 @@ def make_map_exposure_true_energy(pointing, livetime, aeff, geom):
exposure = (exposure * livetime).to('m2 s')

return WcsNDMap(geom, exposure.value, unit=exposure.unit)

def make_map_exposure_reco_energy(exposure_map, spectrum, edisp=None):
"""Create an exposure map in reco energy from an exposure map in true energy.
Exposure in true energy is weighted with an input spectrum and redistributed in
reco energy with the input energy dispersion.
Parameters
----------
exposure_map : `~gammapy.maps.Map`
Input exposure map in true energy. unit should have dimension of m2s
spectrum : `~gammapy.spectrum.models.SpectralModel`, default is None
Spectral model to use to weight exposure in true energy
If None is passed, a power law of photon index -2 is assumed.
edisp : `~gammapy.irf.EnergyDispersion`, default is None.
Energy Dispersion response to redistribute true energies to reco energies
If no energy dispersion is passed, a diagonal one is assumed with identical true and reco
energy bins. This is the default behaviour.
Returns
-------
exposure_image : `~gammapy.maps.Map`
Resulting image in reco energy. The unit is the same as the input map.
"""
expo_map = exposure_map.copy()
etrue_axis = expo_map.geom.get_axis_by_name("energy")
etrue_center = etrue_axis.center * etrue_axis.unit
etrue_edges = etrue_axis.edges * etrue_axis.unit
binsize = np.diff(etrue_edges)

if spectrum is None:
spectrum = PowerLaw(index=2.0)
weights = spectrum(etrue_center)*binsize
weights /= weights.sum()

for img, idx in expo_map.iter_by_image():
img *= weights[idx].value

if edisp is None:
edisp = EnergyDispersion.from_diagonal_matrix(etrue_edges, etrue_edges)

This comment has been minimized.

Copy link
@cdeil

cdeil Aug 8, 2018

Member

Isn't applying the diagonal matrix a no-op?
If so, you could just now make one and do EDISP stuff if no EDISP is given.
Faster and clearer.


# Here we do the convolution with the edisp matrix
a = np.rollaxis(expo_map.data, 0, 3)
b = np.dot(a, edisp.pdf_matrix)

This comment has been minimized.

Copy link
@cdeil

cdeil Aug 8, 2018

Member

We shouldn't duplicate this np.rollaxis, np.dot, np.rollaxis stuff, no?
Already figured out and tested here:
http://docs.gammapy.org/dev/_modules/gammapy/cube/fit.html#MapEvaluator.apply_edisp

Of course, we could also call from MapEvaluator into a function, but it shouldn't be duplicated IMO.

This comment has been minimized.

Copy link
@registerrier

registerrier Aug 9, 2018

Author Contributor

We will have other usages of an apply_edisp taking a Map and an EnergyDispersion as input; for instance, if we want a PSFKernel in reco_energy for morphological analysis.

Therefore, I propose to create a stand-alone function apply_edisp(map, edisp). It can be based on a private function _map_dot_product(array, array) that will just to the axis rotation and dot product on the data. The latter can be also used in MapEvaluator.apply_edisp().

Is this OK @cdeil, @adonath ?

This comment has been minimized.

Copy link
@cdeil

cdeil Aug 9, 2018

Member

I don't have the big picture yet how this function and the MapEvaluator relate.

This is on the line between MapMaker and MapEvaluator, could go either way.

This could also become a method on MapEvaluator, because it's doing Spectrum / npred related computations. Note that you can have @staticmethod just located on MapEvalutor, to avoid 10 standalone functions.

If you want to pursue this as a standalone function, you could also just make MapEvaluator.apply_edisp a `@staticmethod and call it here from now. Then the code duplication with the rollaxis is avoided. Or move it to a separate function and call it from both places. But I'm not super excited about ending up with ~ 10 functions that each just have 3 lines, like the ones we now have on the MapEvaluator.

@registerrier - ?

expo_ereco_data = np.rollaxis(b, 2, 0)

# Create output image by replacing true energy with reco energy axes in initial geom
geom_2D = expo_map.geom.to_image()
ereco_axis = MapAxis(edisp.e_reco.bins, unit=edisp.e_reco.unit, name="energy")

new_axes = [ereco_axis if axis == etrue_axis else axis for axis in expo_map.geom.axes]
geom_ereco = geom_2D.to_cube(axes=new_axes)

ereco_expo_map = Map.from_geom(geom=geom_ereco, unit=expo_map.unit)
ereco_expo_map.data = expo_ereco_data
return ereco_expo_map
Copy path View file
@@ -6,7 +6,7 @@
from astropy.coordinates import Angle
from ..maps import Map, WcsGeom
from .counts import fill_map_counts
from .exposure import make_map_exposure_true_energy
from .exposure import make_map_exposure_true_energy, make_map_exposure_reco_energy
from .background import make_map_background_irf, _fov_background_norm

__all__ = [
@@ -116,6 +116,32 @@ def _process_obs(self, obs, selection):
data = maps_obs[name].quantity.to(self.maps[name].unit).value
self.maps[name].fill_by_coord(coords, data)

def make_images(self, spectrum=None, edisp=None):
"""Create 2D maps by summing over energy axis.
Parameters
----------
spectrum : `~gammapy.spectrum.models.SpectralModel`, default is None
Spectral model to use to weight exposure in true energy
If None is passed, a power law of photon index -2 is assumed.
edisp : `~gammapy.irf.EnergyDispersion`, default is None.
Energy Dispersion response to redistribute true energies to reco energies
If no energy dispersion is passed, a diagonal one is assumed with identical true and reco
energy bins. This is the default behaviour.
Returns
-------
images : dict of `~gammapy.maps.Map`
the dictionary of images
"""
images = dict()
for name, map in self.maps.items():
if name == 'exposure':
expo_reco = make_map_exposure_reco_energy(map, spectrum, edisp)
images[name] = expo_reco.sum_over_axes()
images[name] = map.sum_over_axes()

This comment has been minimized.

Copy link
@cdeil

cdeil Aug 8, 2018

Member

Move this line into an else block?

I think otherwise, you're overwriting the exposure image with just the sum image like for counts and background?


return images

class MapMakerObs(object):
"""Make maps for a single IACT observation.
ProTip! Use n and p to navigate between commits in a pull request.
You can’t perform that action at this time.