2PC Fermi catalog

gammapy/catalog/__init__.py

@@ -4,11 +4,13 @@
 from .core import SourceCatalog, SourceCatalogObject
 from .fermi import (
     SourceCatalog2FHL,
+    SourceCatalog2PC,
     SourceCatalog3FGL,
     SourceCatalog3FHL,
     SourceCatalog3PC,
     SourceCatalog4FGL,
     SourceCatalogObject2FHL,
+    SourceCatalogObject2PC,
     SourceCatalogObject3FGL,
     SourceCatalogObject3FHL,
     SourceCatalogObject3PC,
@@ -40,6 +42,7 @@
         SourceCatalog3FHL,
         SourceCatalog3PC,
         SourceCatalog3HWC,
+        SourceCatalog2PC,
         SourceCatalog1LHAASO,
     ]
 )
@@ -56,6 +59,7 @@
     "SourceCatalog3FHL",
     "SourceCatalog3HWC",
     "SourceCatalog4FGL",
+    "SourceCatalog2PC",
     "SourceCatalog3PC",
     "SourceCatalogGammaCat",
     "SourceCatalogHGPS",
@@ -68,6 +72,7 @@
     "SourceCatalogObject3FHL",
     "SourceCatalogObject3HWC",
     "SourceCatalogObject4FGL",
+    "SourceCatalogObject2PC",
     "SourceCatalogObject3PC",
     "SourceCatalogObjectGammaCat",
     "SourceCatalogObjectHGPS",

gammapy/catalog/fermi.py

@@ -28,11 +28,13 @@
     "SourceCatalog3FGL",
     "SourceCatalog3FHL",
     "SourceCatalog4FGL",
+    "SourceCatalog2PC",
     "SourceCatalog3PC",
     "SourceCatalogObject2FHL",
     "SourceCatalogObject3FGL",
     "SourceCatalogObject3FHL",
     "SourceCatalogObject4FGL",
+    "SourceCatalogObject2PC",
     "SourceCatalogObject3PC",
 ]
 
@@ -1326,6 +1328,202 @@ def spatial_model(self):
         return model
 
 
+class SourceCatalogObject2PC(SourceCatalogObjectFermiPCBase):
+    @property
+    def _auxiliary_filename(self):
+        return make_path(
+            "$GAMMAPY_DATA/catalogs/fermi/2PC_auxiliary/PSR"
+            + self.name
+            + "_2PC_data.fits.gz"
+        )
+
+    def _info_more(self):
+        d = self.data
+        ss = "\n*** Other info ***\n\n"
+        ss += "{:<20s} : {:s}\n".format("Binary", d["Binary"])
+        return ss
+
+    def _info_pulsar(self):
+        d = self.data
+        ss = "\n*** Pulsar info ***\n\n"
+        ss += "{:<20s} : {:.3f}\n".format("Period", d["Period"])
+        ss += "{:<20s} : {:.3e}\n".format("P_Dot", d["P_Dot"])
+        ss += "{:<20s} : {:.3e}\n".format("E_Dot", d["E_Dot"])
+        ss += "{:<20s} : {}\n".format("Type", d["Type"])
+        return ss
+
+    def _info_spectral_fit(self):
+        d = self.data_spectral
+        ss = "\n*** Spectral info ***\n\n"
+        if d is None:
+            ss += "No spectral info available.\n"
+            return ss
+        ss += "{:<20s} : {}\n".format("On peak", d["On_Peak"])
+        ss += "{:<20s} : {:.0f}\n".format("TS DC", d["TS_DC"])
+        ss += "{:<20s} : {:.0f}\n".format("TS cutoff", d["TS_Cutoff"])
+        ss += "{:<20s} : {:.0f}\n".format("TS b free", d["TS_bfree"])
+
+        indentation = " " * 4
+        fmt_e = "{}{:<20s} : {:.3e} +- {:.3e}\n"
+        fmt_f = "{}{:<20s} : {:.3f} +- {:.3f}\n"
+
+        if not isinstance(d["PLEC1_Prefactor"], np.ma.core.MaskedConstant):
+
+            ss += "\n{}* PLSuperExpCutoff b = 1 *\n\n".format(indentation)
+            ss += fmt_e.format(
+                indentation, "Amplitude", d["PLEC1_Prefactor"], d["Unc_PLEC1_Prefactor"]
+            )
+            ss += fmt_f.format(
+                indentation,
+                "Index 1",
+                d["PLEC1_Photon_Index"],
+                d["Unc_PLEC1_Photon_Index"],
+            )
+            ss += "{}{:<20s} : {:.3f}\n".format(indentation, "Index 2", 1)
+            ss += "{}{:<20s} : {:.3f}\n".format(
+                indentation, "Reference", d["PLEC1_Scale"]
+            )
+            ss += fmt_f.format(
+                indentation, "Ecut", d["PLEC1_Cutoff"], d["Unc_PLEC1_Cutoff"]
+            )
+
+        if not isinstance(d["PLEC_Prefactor"], np.ma.core.MaskedConstant):
+
+            ss += "\n{}* PLSuperExpCutoff b free *\n\n".format(indentation)
+            ss += fmt_e.format(
+                indentation, "Amplitude", d["PLEC_Prefactor"], d["Unc_PLEC_Prefactor"]
+            )
+            ss += fmt_f.format(
+                indentation,
+                "Index 1",
+                d["PLEC_Photon_Index"],
+                d["Unc_PLEC_Photon_Index"],
+            )
+            ss += fmt_f.format(
+                indentation,
+                "Index 2",
+            )
+            d["PLEC_Exponential_Index"],
+            d["Unc_PLEC_Exponential_Index"],
+            ss += "{}{:<20s} : {:.3f}\n".format(
+                indentation, "Reference", d["PLEC_Scale"]
+            )
+            ss += fmt_f.format(
+                indentation, "Ecut", d["PLEC_Cutoff"], d["Unc_PLEC_Cutoff"]
+            )
+
+        if not isinstance(d["PL_Prefactor"], np.ma.core.MaskedConstant):
+
+            ss += "\n{}* PowerLaw *\n\n".format(indentation)
+            ss += fmt_e.format(
+                indentation, "Amplitude", d["PL_Prefactor"], d["Unc_PL_Prefactor"]
+            )
+            ss += fmt_f.format(
+                indentation, "Index", d["PL_Photon_Index"], d["Unc_PL_Photon_Index"]
+            )
+            ss += "{}{:<20s} : {:.3f}\n".format(indentation, "Reference", d["PL_Scale"])
+
+        return ss
+
+    def _info_phasogram(self):
+        d = self.data
+        ss = "\n*** Phasogram info ***\n\n"
+        ss += "{:<20s} : {:d}\n".format("Number of peaks", d["Num_Peaks"])
+        ss += "{:<20s} : {:.3f}\n".format("Peak separation", d["Peak_Sep"])
+        return ss
+
+    def spectral_model(self):
+        d = self.data_spectral
+        if d is None:
+            log.warning(f"No spectral model available for source {self.name}")
+            return None
+        if d["TS_Cutoff"] < 9:
+            tag = "PowerLawSpectralModel"
+            pars = {
+                "reference": d["PL_Scale"],
+                "amplitude": d["PL_Prefactor"],
+                "index": d["PL_Photon_Index"],
+            }
+            errs = {
+                "amplitude": d["Unc_PL_Prefactor"],
+                "index": d["Unc_PL_Photon_Index"],
+            }
+        elif d["TS_bfree"] >= 9:
+            tag = "SuperExpCutoffPowerLaw3FGLSpectralModel"
+            pars = {
+                "index_1": d["PLEC_Photon_Index"],
+                "index_2": d["PLEC_Exponential_Index"],
+                "amplitude": d["PLEC_Prefactor"],
+                "reference": d["PLEC_Scale"],
+                "ecut": d["PLEC_Cutoff"],
+            }
+            errs = {
+                "index_1": d["Unc_PLEC_Photon_Index"],
+                "index_2": d["Unc_PLEC_Exponential_Index"],
+                "amplitude": d["Unc_PLEC_Prefactor"],
+                "ecut": d["Unc_PLEC_Cutoff"],
+            }
+        elif d["TS_bfree"] < 9:
+            tag = "SuperExpCutoffPowerLaw3FGLSpectralModel"
+            pars = {
+                "index_1": d["PLEC1_Photon_Index"],
+                "index_2": 1,
+                "amplitude": d["PLEC1_Prefactor"],
+                "reference": d["PLEC1_Scale"],
+                "ecut": d["PLEC1_Cutoff"],
+            }
+            errs = {
+                "index_1": d["Unc_PLEC1_Photon_Index"],
+                "amplitude": d["Unc_PLEC1_Prefactor"],
+                "ecut": d["Unc_PLEC1_Cutoff"],
+            }
+        else:
+            log.warning(f"No spectral model available for source {self.name}")
+            return None
+
+        model = Model.create(tag, "spectral", **pars)
+
+        for name, value in errs.items():
+            model.parameters[name].error = value
+
+        return model
+
+    @property
+    def flux_points_table(self):
+        """Flux points (`~astropy.table.Table`)."""
+        try:
+            fp_data = Table.read(self._auxiliary_filename, hdu="PULSAR_SED")
+        except (KeyError, FileNotFoundError):
+            log.warning(f"No flux points available for source {self.name}")
+            return None
+        table = Table()
+
+        table["e_min"] = fp_data["Energy_Min"]
+        table["e_max"] = fp_data["Energy_Max"]
+        table["e_ref"] = fp_data["Center_Energy"]
+
+        table["flux"] = fp_data["PhotonFlux"]
+        table["flux_err"] = fp_data["Unc_PhotonFlux"]
+
+        table["e2dnde"] = fp_data["EnergyFlux"]
+        table["e2dnde_err"] = fp_data["Unc_EnergyFlux"]
+
+        is_ul = np.where(table["e2dnde_err"] == 0, True, False)
+        table["is_ul"] = is_ul
+
+        table["flux_ul"] = np.nan * table["flux_err"].unit
+        flux_ul = compute_flux_points_ul(table["flux"], table["flux_err"])
+        table["flux_ul"][is_ul] = flux_ul[is_ul]
+
+        table["e2dnde_ul"] = np.nan * table["e2dnde"].unit
+        e2dnde_ul = compute_flux_points_ul(table["e2dnde"], table["e2dnde_err"])
+        table["e2dnde_ul"][is_ul] = e2dnde_ul[is_ul]
+
+        table = Table()
+
+        return table
+
+
 class SourceCatalogObject3PC(SourceCatalogObjectFermiPCBase):
 
     asso = ["assoc_new"]
@@ -1686,6 +1884,56 @@ def __init__(self, filename="$GAMMAPY_DATA/catalogs/fermi/gll_psch_v13.fit.gz"):
         self.energy_bounds_table = Table.read(filename, hdu="EnergyBounds")
 
 
+class SourceCatalog2PC(SourceCatalog):
+    """Fermi-LAT 2nd pulsar catalog.
+
+    - https://ui.adsabs.harvard.edu/abs/2013ApJS..208...17A
+    - https://fermi.gsfc.nasa.gov/ssc/data/access/lat/2nd_PSR_catalog/
+
+    One source is represented by `~gammapy.catalog.SourceCatalogObject2PC`.
+
+    TODO  : Fix the UnitsWarning here ...
+    """
+
+    tag = "2PC"
+    description = "LAT 2nd pulsar catalog"
+    source_object_class = SourceCatalogObject2PC
+
+    def __init__(self, filename="$GAMMAPY_DATA/catalogs/fermi/2PC_catalog_v04.fits.gz"):
+
+        filename = make_path(filename)
+
+        with warnings.catch_warnings():
+            warnings.simplefilter("ignore", u.UnitsWarning)
+            table_psr = Table.read(filename, hdu="PULSAR_CATALOG")
+            table_spectral = Table.read(filename, hdu="SPECTRAL")
+            table_off_peak = Table.read(filename, hdu="OFF_PEAK")
+
+        table_standardise_units_inplace(table_psr)
+        table_standardise_units_inplace(table_spectral)
+        table_standardise_units_inplace(table_off_peak)
+
+        source_name_key = "PSR_Name"
+
+        super().__init__(table=table_psr, source_name_key=source_name_key)
+
+        self.source_object_class._source_name_key = source_name_key
+
+        self.off_peak_table = table_off_peak
+        self.spectral_table = table_spectral
+
+    def to_models(self, **kwargs):
+        models = Models()
+        for m in self:
+            sky_model = m.sky_model()
+            if sky_model is not None:
+                models.append(sky_model)
+        return models
+
+    def _get_name_spectral(self, data):
+        return f"{data[self._source_name_key].strip()}"
+
+
 class SourceCatalog3PC(SourceCatalog):
     """Fermi-LAT 3rd pulsar catalog.
 

gammapy/catalog/tests/data/2pc_j0034-0534.txt

@@ -0,0 +1,75 @@
+
+*** Basic info ***
+
+Catalog row index (zero-based) : 3
+Source name          : J0034-0534
+
+*** Other info ***
+
+Binary               : Y
+
+*** Pulsar info ***
+
+Period               : 1.880 ms
+P_Dot                : 4.980e-21
+E_Dot                : 1.720e+34 erg / s
+Type                 : MSP
+
+*** Position info ***
+
+RA                   : 8.591 deg
+DEC                  : -5.577 deg
+GLON                 : 111.492 deg
+GLAT                 : -68.069 deg
+
+*** Spectral info ***
+
+On peak              : N
+TS DC                : 563
+TS cutoff            : 45
+TS b free            : 0
+
+    * PLSuperExpCutoff b = 1 *
+
+    Amplitude            : 8.691e-12 1 / (MeV s cm2) +- 1.208e-12 1 / (MeV s cm2)
+    Index 1              : 1.436 +- 0.169
+    Index 2              : 1.000
+    Reference            : 755.700 MeV
+    Ecut                 : 1831.000 MeV +- 416.700 MeV
+
+    * PLSuperExpCutoff b free *
+
+    Amplitude            : 0.000e+00 1 / (MeV s cm2) +- 0.000e+00 1 / (MeV s cm2)
+    Index 1              : -- +- --
+    Index 2              : -- +- --
+    Reference            : 0.000 MeV
+    Ecut                 : 0.000 MeV +- 0.000 MeV
+
+    * PowerLaw *
+
+    Amplitude            : 2.435e-12 1 / (MeV s cm2) +- 1.520e-13 1 / (MeV s cm2)
+    Index                : 2.219 +- 0.049
+    Reference            : 1000.000 MeV
+
+*** Spectral points ***
+
+e_min   e_max  e_ref        flux           flux_err         e2dnde      e2dnde_err  is_ul     flux_ul        e2dnde_ul  
+ GeV     GeV    GeV   ph / (GeV s cm2) ph / (GeV s cm2) erg / (s cm2) erg / (s cm2)       ph / (GeV s cm2) erg / (s cm2)
+------ ------- ------ ---------------- ---------------- ------------- ------------- ----- ---------------- -------------
+ 0.100   0.178  0.128        1.747e-07        0.000e+00     4.586e-12     0.000e+00  True        1.747e-07     4.586e-12
+ 0.178   0.316  0.228        5.343e-08        1.041e-08     4.436e-12     8.638e-13 False              nan           nan
+ 0.316   0.562  0.405        1.886e-08        2.588e-09     4.950e-12     6.794e-13 False              nan           nan
+ 0.562   1.000  0.720        5.558e-09        6.914e-10     4.614e-12     5.740e-13 False              nan           nan
+ 1.000   1.778  1.280        1.927e-09        2.320e-10     5.059e-12     6.090e-13 False              nan           nan
+ 1.778   3.162  2.276        5.261e-10        7.922e-11     4.368e-12     6.576e-13 False              nan           nan
+ 3.162   5.623  4.048        8.004e-11        2.215e-11     2.101e-12     5.816e-13 False              nan           nan
+ 5.623  10.000  7.199        1.634e-11        0.000e+00     1.356e-12     0.000e+00  True        1.634e-11     1.356e-12
+10.000  17.783 12.801        3.678e-12        0.000e+00     9.655e-13     0.000e+00  True        3.678e-12     9.655e-13
+17.783  31.623 22.764        1.158e-12        0.000e+00     9.616e-13     0.000e+00  True        1.158e-12     9.616e-13
+31.623  56.234 40.482        6.298e-13        0.000e+00     1.653e-12     0.000e+00  True        6.298e-13     1.653e-12
+56.234 100.000 71.988        9.865e-13        0.000e+00     8.190e-12     0.000e+00  True        9.865e-13     8.190e-12
+
+*** Phasogram info ***
+
+Number of peaks      : 2
+Peak separation      : 0.285