Source code for gammapy.irf.io

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import logging
from astropy.io import fits
from gammapy.data.hdu_index_table import HDUIndexTable
from gammapy.utils.fits import HDULocation
from gammapy.utils.scripts import make_path

__all__ = ["load_cta_irfs", "load_irf_dict_from_file"]

log = logging.getLogger(__name__)


IRF_DL3_AXES_SPECIFICATION = {
    "THETA": {"name": "offset", "interp": "lin"},
    "ENERG": {"name": "energy_true", "interp": "log"},
    "ETRUE": {"name": "energy_true", "interp": "log"},
    "RAD": {"name": "rad", "interp": "lin"},
    "DETX": {"name": "fov_lon", "interp": "lin"},
    "DETY": {"name": "fov_lat", "interp": "lin"},
    "MIGRA": {"name": "migra", "interp": "lin"},
}

COMMON_HEADERS = {
    "HDUCLASS": "GADF",
    "HDUDOC": "https://github.com/open-gamma-ray-astro/gamma-astro-data-formats",
    "HDUVERS": "0.2",
}

COMMON_IRF_HEADERS = {
    **COMMON_HEADERS,
    "HDUCLAS1": "RESPONSE",
}


# The key is the class tag.
# TODO: extend the info here with the minimal header info
IRF_DL3_HDU_SPECIFICATION = {
    "bkg_3d": {
        "extname": "BACKGROUND",
        "column_name": "BKG",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "BKG",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "BKG_3D",
        },
    },
    "bkg_2d": {
        "extname": "BACKGROUND",
        "column_name": "BKG",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "BKG",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "BKG_2D",
        },
    },
    "edisp_2d": {
        "extname": "ENERGY DISPERSION",
        "column_name": "MATRIX",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "EDISP",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "EDISP_2D",
        },
    },
    "psf_table": {
        "extname": "PSF_2D_TABLE",
        "column_name": "RPSF",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "RPSF",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "PSF_TABLE",
        },
    },
    "psf_3gauss": {
        "extname": "PSF_2D_GAUSS",
        "column_name": {
            "sigma_1": "SIGMA_1",
            "sigma_2": "SIGMA_2",
            "sigma_3": "SIGMA_3",
            "scale": "SCALE",
            "ampl_2": "AMPL_2",
            "ampl_3": "AMPL_3",
        },
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "RPSF",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "PSF_3GAUSS",
        },
    },
    "psf_king": {
        "extname": "PSF_2D_KING",
        "column_name": {
            "sigma": "SIGMA",
            "gamma": "GAMMA",
        },
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "RPSF",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "PSF_KING",
        },
    },
    "aeff_2d": {
        "extname": "EFFECTIVE AREA",
        "column_name": "EFFAREA",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "EFF_AREA",
            "HDUCLAS3": "FULL-ENCLOSURE",  # added here to have HDUCLASN in order
            "HDUCLAS4": "AEFF_2D",
        },
    },
    "rad_max_2d": {
        "extname": "RAD_MAX",
        "column_name": "RAD_MAX",
        "mandatory_keywords": {
            **COMMON_IRF_HEADERS,
            "HDUCLAS2": "RAD_MAX",
            "HDUCLAS3": "POINT-LIKE",
            "HDUCLAS4": "RAD_MAX_2D",
        },
    },
}


IRF_MAP_HDU_SPECIFICATION = {
    "edisp_kernel_map": "edisp",
    "edisp_map": "edisp",
    "psf_map": "psf",
    "psf_map_reco": "psf",
}


def gadf_is_pointlike(header):
    """Check if a GADF IRF is pointlike based on the header"""
    return header.get("HDUCLAS3") == "POINT-LIKE"


[docs]def load_cta_irfs(filename): """load CTA instrument response function and return a dictionary container. The IRF format should be compliant with the one discussed at http://gamma-astro-data-formats.readthedocs.io/en/latest/irfs/. The various IRFs are accessible with the following keys: - 'aeff' is a `~gammapy.irf.EffectiveAreaTable2D` - 'edisp' is a `~gammapy.irf.EnergyDispersion2D` - 'psf' is a `~gammapy.irf.EnergyDependentMultiGaussPSF` - 'bkg' is a `~gammapy.irf.Background3D` Parameters ---------- filename : str the input filename. Default is Returns ------- cta_irf : dict the IRF dictionary Examples -------- Access the CTA 1DC IRFs stored in the gammapy datasets >>> from gammapy.irf import load_cta_irfs >>> filename = "$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits" >>> cta_irf = load_cta_irfs(filename) >>> print(cta_irf['aeff']) EffectiveAreaTable2D -------------------- <BLANKLINE> axes : ['energy_true', 'offset'] shape : (42, 6) ndim : 2 unit : m2 dtype : >f4 <BLANKLINE> """ from .background import Background3D from .edisp import EnergyDispersion2D from .effective_area import EffectiveAreaTable2D from .psf import EnergyDependentMultiGaussPSF aeff = EffectiveAreaTable2D.read(filename, hdu="EFFECTIVE AREA") bkg = Background3D.read(filename, hdu="BACKGROUND") edisp = EnergyDispersion2D.read(filename, hdu="ENERGY DISPERSION") psf = EnergyDependentMultiGaussPSF.read(filename, hdu="POINT SPREAD FUNCTION") return dict(aeff=aeff, bkg=bkg, edisp=edisp, psf=psf)
[docs]def load_irf_dict_from_file(filename): """Open a fits file and generate a dictionary containing the Gammapy objects corresponding to the IRF components stored Parameters ---------- filename : str, Path path to the file containing the IRF components, if EVENTS and GTI HDUs are included in the file, they are ignored Returns ------- irf_dict : dict of `~gammapy.irf.IRF` dictionary with instances of the Gammapy objects corresponding to the IRF components """ from .rad_max import RadMax2D filename = make_path(filename) hdulist = fits.open(make_path(filename)) irf_dict = {} is_pointlike = False for hdu in hdulist: hdu_class = hdu.header.get("HDUCLAS1", "").lower() if hdu_class == "response": hdu_class = hdu.header.get("HDUCLAS4", "").lower() is_pointlike |= hdu.header["HDUCLAS3"] == "POINT-LIKE" loc = HDULocation( hdu_class=hdu_class, hdu_name=hdu.name, file_dir=filename.parent, file_name=filename.name, ) for name in HDUIndexTable.VALID_HDU_TYPE: if name in hdu_class: if name in irf_dict.keys(): log.warning(f"more than one HDU of {name} type found") log.warning( f"loaded the {irf_dict[name].meta['EXTNAME']} HDU in the dictionary" ) continue data = loc.load() # TODO: maybe introduce IRF.type attribute... irf_dict[name] = data else: # not an IRF component continue if is_pointlike and "rad_max" not in irf_dict: irf_dict["rad_max"] = RadMax2D.from_irf(irf_dict["aeff"]) return irf_dict