# Licensed under a 3-clause BSD style license - see LICENSE.rstimportloggingfromastropy.ioimportfitsfromgammapy.data.hdu_index_tableimportHDUIndexTablefromgammapy.utils.deprecationimportdeprecatedfromgammapy.utils.fitsimportHDULocationfromgammapy.utils.scriptsimportmake_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 infoIRF_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","FOVALIGN":"RADEC",},},"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",}defgadf_is_pointlike(header):"""Check if a GADF IRF is pointlike based on the header"""returnheader.get("HDUCLAS3")=="POINT-LIKE"
[docs]@deprecated("v1.1",alternative="load_irf_dict_from_file")defload_cta_irfs(filename):"""Load IRFs from file as written by the CTA DC1 into a dict This function has a hardcoded list of IRF types and HDU names and does not check what types of IRFs are actually present in the file. Please use `load_irf_dict_from_file` instead.. 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.backgroundimportBackground3Dfrom.edispimportEnergyDispersion2Dfrom.effective_areaimportEffectiveAreaTable2Dfrom.psfimportEnergyDependentMultiGaussPSFaeff=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")returndict(aeff=aeff,bkg=bkg,edisp=edisp,psf=psf)
classUnknownHDUClass(IOError):"""Raised when a file contains an unknown HDUCLASS"""def_get_hdu_type_and_class(header):"""Get gammapy hdu_type and class from FITS header Contains a workaround to support CTA 1DC irf file. """hdu_clas2=header.get("HDUCLAS2","")hdu_clas4=header.get("HDUCLAS4","")clas2_to_type={"rpsf":"psf","eff_area":"aeff"}hdu_type=clas2_to_type.get(hdu_clas2.lower(),hdu_clas2.lower())hdu_class=hdu_clas4.lower()ifhdu_typenotinHDUIndexTable.VALID_HDU_TYPE:raiseUnknownHDUClass(f"HDUCLAS2={hdu_clas2}, HDUCLAS4={hdu_clas4}")# workaround for CTA 1DC files with non-compliant HDUCLAS4 namesifhdu_classnotinHDUIndexTable.VALID_HDU_CLASS:hdu_class=f"{hdu_type}_{hdu_class}"ifhdu_classnotinHDUIndexTable.VALID_HDU_CLASS:raiseUnknownHDUClass(f"HDUCLAS2={hdu_clas2}, HDUCLAS4={hdu_clas4}")returnhdu_type,hdu_class
[docs]defload_irf_dict_from_file(filename):"""Load all available IRF components from given file into a dict. If multiple IRFs of the same type are present, the first encountered is returned. 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_maximportRadMax2Dfilename=make_path(filename)irf_dict={}is_pointlike=Falsewithfits.open(filename)ashdulist:forhduinhdulist:hdu_clas1=hdu.header.get("HDUCLAS1","").lower()# not an IRF componentifhdu_clas1!="response":continueis_pointlike|=hdu.header.get("HDUCLAS3")=="POINT-LIKE"try:hdu_type,hdu_class=_get_hdu_type_and_class(hdu.header)exceptUnknownHDUClassase:log.warning("File has unknown class %s",e)continueloc=HDULocation(hdu_class=hdu_class,hdu_name=hdu.name,file_dir=filename.parent,file_name=filename.name,)ifhdu_typeinirf_dict.keys():log.warning(f"more than one HDU of {hdu_type} type found")log.warning(f"loaded the {irf_dict[hdu_type].meta['EXTNAME']} HDU in the dictionary")continuedata=loc.load()irf_dict[hdu_type]=dataifis_pointlikeand"rad_max"notinirf_dict:irf_dict["rad_max"]=RadMax2D.from_irf(irf_dict["aeff"])returnirf_dict