# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
import astropy.units as u
from astropy.table import Table
from gammapy.irf import EDispKernelMap, PSFMap
from gammapy.maps import Map
from .core import Maker
from .utils import (
make_edisp_kernel_map,
make_edisp_map,
make_map_background_irf,
make_map_exposure_true_energy,
make_psf_map,
)
__all__ = ["MapDatasetMaker"]
log = logging.getLogger(__name__)
[docs]class MapDatasetMaker(Maker):
"""Make maps for a single IACT observation.
Parameters
----------
selection : list
List of str, selecting which maps to make.
Available: 'counts', 'exposure', 'background', 'psf', 'edisp'
By default, all maps are made.
background_oversampling : int
Background evaluation oversampling factor in energy.
background_interp_missing_data: bool
Interpolate missing values in background 3d map.
Default is True, have to be set to True for CTA IRF.
"""
tag = "MapDatasetMaker"
available_selection = ["counts", "exposure", "background", "psf", "edisp"]
def __init__(
self,
selection=None,
background_oversampling=None,
background_interp_missing_data=True,
):
self.background_oversampling = background_oversampling
self.background_interp_missing_data = background_interp_missing_data
if selection is None:
selection = self.available_selection
selection = set(selection)
if not selection.issubset(self.available_selection):
difference = selection.difference(self.available_selection)
raise ValueError(f"{difference} is not a valid method.")
self.selection = selection
[docs] @staticmethod
def make_counts(geom, observation):
"""Make counts map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference map geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
counts : `~gammapy.maps.Map`
Counts map.
"""
counts = Map.from_geom(geom)
counts.fill_events(observation.events)
return counts
[docs] @staticmethod
def make_exposure(geom, observation, use_region_center=True):
"""Make exposure map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference map geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
exposure : `~gammapy.maps.Map`
Exposure map.
"""
if isinstance(observation.aeff, Map):
return observation.aeff.interp_to_geom(
geom=geom,
)
return make_map_exposure_true_energy(
pointing=observation.pointing_radec,
livetime=observation.observation_live_time_duration,
aeff=observation.aeff,
geom=geom,
use_region_center=use_region_center,
)
[docs] @staticmethod
def make_exposure_irf(geom, observation, use_region_center=True):
"""Make exposure map with irf geometry.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
exposure : `~gammapy.maps.Map`
Exposure map.
"""
return make_map_exposure_true_energy(
pointing=observation.pointing_radec,
livetime=observation.observation_live_time_duration,
aeff=observation.aeff,
geom=geom,
use_region_center=use_region_center,
)
[docs] def make_background(self, geom, observation):
"""Make background map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
background : `~gammapy.maps.Map`
Background map.
"""
bkg = observation.bkg
if isinstance(bkg, Map):
return bkg.interp_to_geom(geom=geom, preserve_counts=True)
bkg_coordsys = observation.bkg.meta.get("FOVALIGN", "RADEC")
if bkg_coordsys == "ALTAZ":
pointing = observation.fixed_pointing_info
elif bkg_coordsys == "RADEC":
pointing = observation.pointing_radec
else:
raise ValueError(
f"Invalid background coordinate system: {bkg_coordsys!r}\n"
"Options: ALTAZ, RADEC"
)
use_region_center = getattr(self, "use_region_center", True)
if self.background_interp_missing_data:
bkg.interp_missing_data(axis_name="energy")
return make_map_background_irf(
pointing=pointing,
ontime=observation.observation_time_duration,
bkg=bkg,
geom=geom,
oversampling=self.background_oversampling,
use_region_center=use_region_center,
)
[docs] def make_edisp(self, geom, observation):
"""Make energy dispersion map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
edisp : `~gammapy.irf.EDispMap`
Edisp map.
"""
exposure = self.make_exposure_irf(geom.squash(axis_name="migra"), observation)
use_region_center = getattr(self, "use_region_center", True)
return make_edisp_map(
edisp=observation.edisp,
pointing=observation.pointing_radec,
geom=geom,
exposure_map=exposure,
use_region_center=use_region_center,
)
[docs] def make_edisp_kernel(self, geom, observation):
"""Make energy dispersion kernel map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference geom. Must contain "energy" and "energy_true" axes in that order.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
edisp : `~gammapy.irf.EDispKernelMap`
EdispKernel map.
"""
if isinstance(observation.edisp, EDispKernelMap):
exposure = None
interp_map = observation.edisp.edisp_map.interp_to_geom(geom)
return EDispKernelMap(edisp_kernel_map=interp_map, exposure_map=exposure)
exposure = self.make_exposure_irf(geom.squash(axis_name="energy"), observation)
use_region_center = getattr(self, "use_region_center", True)
return make_edisp_kernel_map(
edisp=observation.edisp,
pointing=observation.pointing_radec,
geom=geom,
exposure_map=exposure,
use_region_center=use_region_center,
)
[docs] def make_psf(self, geom, observation):
"""Make psf map.
Parameters
----------
geom : `~gammapy.maps.Geom`
Reference geom.
observation : `~gammapy.data.Observation`
Observation container.
Returns
-------
psf : `~gammapy.irf.PSFMap`
Psf map.
"""
psf = observation.psf
if isinstance(psf, PSFMap):
return PSFMap(psf.psf_map.interp_to_geom(geom))
exposure = self.make_exposure_irf(geom.squash(axis_name="rad"), observation)
return make_psf_map(
psf=psf,
pointing=observation.pointing_radec,
geom=geom,
exposure_map=exposure,
)
[docs] def run(self, dataset, observation):
"""Make map dataset.
Parameters
----------
dataset : `~gammapy.datasets.MapDataset`
Reference dataset.
observation : `~gammapy.data.Observation`
Observation
Returns
-------
dataset : `~gammapy.datasets.MapDataset`
Map dataset.
"""
kwargs = {"gti": observation.gti}
kwargs["meta_table"] = self.make_meta_table(observation)
mask_safe = Map.from_geom(dataset.counts.geom, dtype=bool)
mask_safe.data[...] = True
kwargs["mask_safe"] = mask_safe
if "counts" in self.selection:
counts = self.make_counts(dataset.counts.geom, observation)
else:
counts = Map.from_geom(dataset.counts.geom, data=0)
kwargs["counts"] = counts
if "exposure" in self.selection:
exposure = self.make_exposure(dataset.exposure.geom, observation)
kwargs["exposure"] = exposure
if "background" in self.selection:
kwargs["background"] = self.make_background(
dataset.counts.geom, observation
)
if "psf" in self.selection:
psf = self.make_psf(dataset.psf.psf_map.geom, observation)
kwargs["psf"] = psf
if "edisp" in self.selection:
if dataset.edisp.edisp_map.geom.axes[0].name.upper() == "MIGRA":
edisp = self.make_edisp(dataset.edisp.edisp_map.geom, observation)
else:
edisp = self.make_edisp_kernel(
dataset.edisp.edisp_map.geom, observation
)
kwargs["edisp"] = edisp
return dataset.__class__(name=dataset.name, **kwargs)