Source code for gammapy.catalog.lhaaso

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import numpy as np
import astropy.units as u
from astropy.table import Table
from gammapy.maps import MapAxis, RegionGeom
from gammapy.modeling.models import Model, Models, SkyModel
from gammapy.utils.gauss import Gauss2DPDF
from gammapy.utils.scripts import make_path
from .core import SourceCatalog, SourceCatalogObject

__all__ = [
    "SourceCatalogObject1LHAASO",
    "SourceCatalog1LHAASO",
]


[docs]class SourceCatalogObject1LHAASO(SourceCatalogObject): """One source from the 1LHAASO catalog. Catalog is represented by `~gammapy.catalog.SourceCatalog1LHAASO`. """ _source_name_key = "Source_Name" def _parse(self, name, which): if which in self.data["Model_a"]: tag = "" elif which in self.data["Model_b"]: tag = "_b" else: raise ValueError("Invalid model component name") is_ul = False value = u.Quantity(self.data[f"{name}{tag}"]) if ( np.isnan(value) or value == 0 * value.unit ) and f"{name}_ul{tag}" in self.data: value = self.data[f"{name}_ul{tag}"] is_ul = True return value, is_ul def _get(self, name, which): value, _ = self._parse(name, which) return value
[docs] def spectral_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) for name, value in errs.items(): model.parameters[name].error = value return model
[docs] def spatial_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) if not is_ul: pars["sigma"] = r39 model = Model.create("GaussianSpatialModel", "spatial", **pars) else: model = Model.create("PointSpatialModel", "spatial", **pars) for name, value in errs.items(): model.parameters[name].error = value return model
@staticmethod def _get_components_geom(models): energy_axis = MapAxis.from_energy_bounds( "0.1 TeV", "2000 TeV", nbin=10, per_decade=True, name="energy" ) regions = [m.spatial_model.evaluation_region for m in models] geom = RegionGeom.from_regions( regions, binsz_wcs="0.05 deg", axes=[energy_axis] ) return geom.to_wcs_geom()
[docs] def sky_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. """ if which == "both": wcda = self.sky_model(which="WCDA") km2a = self.sky_model(which="KM2A") models = [m for m in [wcda, km2a] if m is not None] if len(models) == 2: geom = self._get_components_geom(models) mask = geom.energy_mask(energy_max=25 * u.TeV) geom = geom.as_energy_true wcda_map = Models(wcda).to_template_sky_model(geom).spatial_model.map model = 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" return model else: return models[0] else: _, is_ul = self._parse("N0", which) if is_ul: return None else: return SkyModel( spatial_model=self.spatial_model(which), spectral_model=self.spectral_model(which), name=self.name, )
[docs]class SourceCatalog1LHAASO(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 = SourceCatalogObject1LHAASO def __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] def to_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 _ in self: model = _.sky_model(which) if model: models.append(model) return models