# Licensed under a 3-clause BSD style license - see LICENSE.rstimportnumpyasnpimportastropy.unitsasufromastropy.tableimportTablefromgammapy.mapsimportMapAxis,RegionGeomfromgammapy.modeling.modelsimportModel,Models,SkyModelfromgammapy.utils.gaussimportGauss2DPDFfromgammapy.utils.scriptsimportmake_pathfrom.coreimportSourceCatalog,SourceCatalogObject__all__=["SourceCatalogObject1LHAASO","SourceCatalog1LHAASO",]
[docs]classSourceCatalogObject1LHAASO(SourceCatalogObject):"""One source from the 1LHAASO catalog. Catalog is represented by `~gammapy.catalog.SourceCatalog1LHAASO`. """_source_name_key="Source_Name"def_parse(self,name,which):ifwhichinself.data["Model_a"]:tag=""elifwhichinself.data["Model_b"]:tag="_b"else:raiseValueError("Invalid model component name")is_ul=Falsevalue=u.Quantity(self.data[f"{name}{tag}"])if(np.isnan(value)orvalue==0*value.unit)andf"{name}_ul{tag}"inself.data:value=self.data[f"{name}_ul{tag}"]is_ul=Truereturnvalue,is_uldef_get(self,name,which):value,_=self._parse(name,which)returnvalue
[docs]defspectral_model(self,which):"""Spectral model as a `~gammapy.modeling.models.PowerLawSpectralModel` object. * ``which="KM2A"`` - Sky model for KM2A analysis only. * ``which="WCDA"`` - Sky model for WCDA analysis only. """pars={"reference":self._get("E0",which),"amplitude":self._get("N0",which),"index":self._get("gamma",which),}errs={"amplitude":self._get("N0_err",which),"index":self._get("gamma_err",which),}model=Model.create("PowerLawSpectralModel","spectral",**pars)forname,valueinerrs.items():model.parameters[name].error=valuereturnmodel
[docs]defspatial_model(self,which):"""Spatial model as a `~gammapy.modeling.models.SpatialModel` object. * ``which="KM2A"`` - Sky model for KM2A analysis only. * ``which="WCDA"`` - Sky model for WCDA analysis only. """lat_0=self._get("DECJ2000",which)pars={"lon_0":self._get("RAJ2000",which),"lat_0":lat_0,"frame":"fk5"}pos_err=self._get("pos_err",which)scale_r95=Gauss2DPDF().containment_radius(0.95)errs={"lat_0":pos_err/scale_r95,"lon_0":pos_err/scale_r95/np.cos(lat_0),}r39,is_ul=self._parse("r39",which)ifnotis_ul:pars["sigma"]=r39model=Model.create("GaussianSpatialModel","spatial",**pars)else:model=Model.create("PointSpatialModel","spatial",**pars)forname,valueinerrs.items():model.parameters[name].error=valuereturnmodel
[docs]defsky_model(self,which="both"):"""Sky model as a `~gammapy.modeling.models.SkyModel` object. * ``which="both"`` - Create a composite template if both models are available, or, use the available one if only one is present. * ``which="KM2A"`` - Sky model for KM2A analysis if available. * ``which="WCDA"`` - Sky model for WCDA analysis if available. """ifwhich=="both":wcda=self.sky_model(which="WCDA")km2a=self.sky_model(which="KM2A")models=[mformin[wcda,km2a]ifmisnotNone]iflen(models)==2:geom=self._get_components_geom(models)mask=geom.energy_mask(energy_max=25*u.TeV)geom=geom.as_energy_truewcda_map=Models(wcda).to_template_sky_model(geom).spatial_model.mapmodel=Models(km2a).to_template_sky_model(geom,name=km2a.name)model.spatial_model.map.data[mask]=wcda_map.data[mask]model.spatial_model.filename=f"{model.name}.fits"returnmodelelse:returnmodels[0]else:_,is_ul=self._parse("N0",which)ifis_ul:returnNoneelse:returnSkyModel(spatial_model=self.spatial_model(which),spectral_model=self.spectral_model(which),name=self.name,)
[docs]classSourceCatalog1LHAASO(SourceCatalog):"""First LHAASO catalog. One source is represented by `~gammapy.catalog.SourceCatalogObject1LHAASO`. The data is from table 1 in the paper [1]_. The catalog table contains 90 rows / sources. References ---------- .. [1] 1LHAASO: The First LHAASO Catalog of Gamma-Ray Sources, https://ui.adsabs.harvard.edu/abs/2023arXiv230517030C/abstract """tag="1LHAASO""""Catalog name"""description="1LHAASO catalog from the LHAASO observatory""""Catalog description"""source_object_class=SourceCatalogObject1LHAASOdef__init__(self,filename="$GAMMAPY_DATA/catalogs/1LHAASO_catalog.fits"):table=Table.read(make_path(filename))source_name_key="Source_Name"super().__init__(table=table,source_name_key=source_name_key)
[docs]defto_models(self,which="both"):"""Create Models object from catalog. * ``which="both"`` - Use first model or create a composite template if both models are available. * ``which="KM2A"`` - Sky model for KM2A analysis if available. * ``which="WCDA"`` - Sky model for WCDA analysis if available. """models=Models()for_inself:model=_.sky_model(which)ifmodel:models.append(model)returnmodels