Use gammapy in core library

examples/calculate_eventdisplay_irfs.py

@@ -48,6 +48,9 @@
 )
 
 
+from gammapy.maps import MapAxis
+
+
 log = logging.getLogger("pyirf")
 
 
@@ -121,14 +124,14 @@ def main():
 
     # event display uses much finer bins for the theta cut than
     # for the sensitivity
-    theta_bins = add_overflow_bins(
-        create_bins_per_decade(10 ** (-1.9) * u.TeV, 10 ** 2.3005 * u.TeV, 50)
+    energy_axis_theta = MapAxis.from_edges(
+        add_overflow_bins(create_bins_per_decade(10 ** (-1.9) * u.TeV, 10 ** 2.3005 * u.TeV, 50)),
+        name="energy",
     )
     # same bins as event display uses
-    sensitivity_bins = add_overflow_bins(
-        create_bins_per_decade(
-            10 ** -1.9 * u.TeV, 10 ** 2.31 * u.TeV, bins_per_decade=5
-        )
+    energy_axis_sensitivity = MapAxis.from_edges(
+        add_overflow_bins(create_bins_per_decade(10 ** (-1.9) * u.TeV, 10 ** 2.3005 * u.TeV, 5)),
+        name="energy",
     )
 
     # theta cut is 68 percent containmente of the gammas
@@ -139,19 +142,19 @@ def main():
     theta_cuts_coarse = calculate_percentile_cut(
         gammas["theta"][mask_theta_cuts],
         gammas["reco_energy"][mask_theta_cuts],
-        bins=sensitivity_bins,
+        bins=energy_axis_sensitivity.edges,
         min_value=0.05 * u.deg,
         fill_value=0.32 * u.deg,
         max_value=0.32 * u.deg,
         percentile=68,
     )
 
     # interpolate to 50 bins per decade
-    theta_center = bin_center(theta_bins)
-    inter_center = bin_center(sensitivity_bins)
+    theta_center = energy_axis_theta.center
+    inter_center = energy_axis_sensitivity.center
     theta_cuts = table.QTable({
-        "low": theta_bins[:-1],
-        "high": theta_bins[1:],
+        "low": energy_axis_theta.edges_min,
+        "high": energy_axis_theta.edges_max,
         "center": theta_center,
         "cut": np.interp(np.log10(theta_center / u.TeV), np.log10(inter_center / u.TeV), theta_cuts_coarse['cut']),
     })
@@ -166,7 +169,7 @@ def main():
     sensitivity, gh_cuts = optimize_gh_cut(
         gammas,
         background,
-        reco_energy_bins=sensitivity_bins,
+        reco_energy_bins=energy_axis_sensitivity.edges,
         gh_cut_efficiencies=gh_cut_efficiencies,
         op=operator.ge,
         theta_cuts=theta_cuts,
@@ -258,8 +261,8 @@ def main():
     psf = psf_table(
         gammas[gammas["selected_gh"]],
         true_energy_bins,
-        fov_offset_bins=fov_offset_bins,
-        source_offset_bins=source_offset_bins,
+        fov_offset_axis=fov_offset_bins,
+        source_offset_axis=source_offset_bins,
     )
 
     background_rate = background_2d(
@@ -278,7 +281,7 @@ def main():
         psf, true_energy_bins, source_offset_bins, fov_offset_bins,
     ))
     hdus.append(create_rad_max_hdu(
-        theta_cuts["cut"][:, np.newaxis], theta_bins, fov_offset_bins
+        theta_cuts["cut"][:, np.newaxis], energy_axis_theta, fov_offset_bins
     ))
     hdus.append(fits.BinTableHDU(ang_res, name="ANGULAR_RESOLUTION"))
     hdus.append(fits.BinTableHDU(bias_resolution, name="ENERGY_BIAS_RESOLUTION"))

pyirf/gammapy.py

@@ -40,14 +40,14 @@ def create_effective_area_table_2d(
     Returns
     -------
     gammapy.irf.EffectiveAreaTable2D
-
     '''
     offset_axis = _create_offset_axis(fov_offset_bins)
     energy_axis_true = _create_energy_axis_true(true_energy_bins)
-
     return EffectiveAreaTable2D(
-        axes = [energy_axis_true,
-                offset_axis],
+        axes=[
+            energy_axis_true,
+            offset_axis,
+        ],
         data=effective_area,
     )
 
@@ -67,7 +67,7 @@ def create_psf_3d(
 ):
     """
     Create a ``gammapy.irf.PSF3D`` from pyirf outputs.
-    
+
     Parameters
     ----------
     psf: astropy.units.Quantity[(solid angle)^-1]
@@ -84,17 +84,19 @@ def create_psf_3d(
     Returns
     -------
     gammapy.irf.PSF3D
-    
+
     """
     offset_axis = _create_offset_axis(fov_offset_bins)
     energy_axis_true = _create_energy_axis_true(true_energy_bins)
     rad_axis = MapAxis.from_edges(source_offset_bins, name='rad')
 
     return PSF3D(
-        axes = [energy_axis_true,
-                offset_axis,
-                rad_axis],
-        data = psf
+        axes=[
+            energy_axis_true,
+            offset_axis,
+            rad_axis,
+        ],
+        data=psf
     )
 
 
@@ -109,7 +111,7 @@ def create_energy_dispersion_2d(
 ):
     """
     Create a ``gammapy.irf.EnergyDispersion2D`` from pyirf outputs.
-    
+
     Parameters
     ----------
     energy_dispersion: numpy.ndarray
@@ -126,15 +128,17 @@ def create_energy_dispersion_2d(
     Returns
     -------
     gammapy.irf.EnergyDispersion2D
-    
+
     """
     offset_axis = _create_offset_axis(fov_offset_bins)
     energy_axis_true = _create_energy_axis_true(true_energy_bins)
     migra_axis = MapAxis.from_edges(migration_bins, name="migra")
 
     return EnergyDispersion2D(
-        axes = [energy_axis_true,
-                migra_axis,
-                offset_axis],
-        data = energy_dispersion,
+        axes=[
+            energy_axis_true,
+            migra_axis,
+            offset_axis,
+        ],
+        data=energy_dispersion,
     )

pyirf/irf/effective_area.py

@@ -1,5 +1,8 @@
 import numpy as np
 import astropy.units as u
+from gammapy.maps import MapAxis
+from gammapy.irf import EffectiveAreaTable2D
+
 from ..binning import create_histogram_table
 
 
@@ -27,7 +30,12 @@ def effective_area(n_selected, n_simulated, area):
     return (n_selected / n_simulated) * area
 
 
-def effective_area_per_energy(selected_events, simulation_info, true_energy_bins):
+def effective_area_per_energy(
+    selected_events,
+    simulation_info,
+    true_energy_axis: MapAxis,
+    fov_offset_axis: MapAxis,
+):
     """
     Calculate effective area in bins of true energy.
 
@@ -37,21 +45,38 @@ def effective_area_per_energy(selected_events, simulation_info, true_energy_bins
         DL2 events table, required columns for this function: `true_energy`.
     simulation_info: pyirf.simulations.SimulatedEventsInfo
         The overall statistics of the simulated events
-    true_energy_bins: astropy.units.Quantity[energy]
+    true_energy_axis: gammapy.maps.MapAxis
         The bin edges in which to calculate effective area.
+    fov_offset_bins: astropy.units.Quantity[angle]
+        The field of view radial offset axis. This must only have a single bin.
     """
     area = np.pi * simulation_info.max_impact ** 2
 
+    if fov_offset_axis.nbin != 1:
+        raise ValueError('fov_offset_axis must only have a single bin')
+
     hist_selected = create_histogram_table(
-        selected_events, true_energy_bins, "true_energy"
+        selected_events,
+        true_energy_axis.edges,
+        "true_energy"
     )
-    hist_simulated = simulation_info.calculate_n_showers_per_energy(true_energy_bins)
+    hist_simulated = simulation_info.calculate_n_showers_per_energy(true_energy_axis.edges)
 
-    return effective_area(hist_selected["n"], hist_simulated, area)
+    aeff = effective_area(hist_selected["n"], hist_simulated, area)
+    return EffectiveAreaTable2D(
+        axes=[
+            true_energy_axis,
+            fov_offset_axis,
+        ],
+        data=aeff[:, np.newaxis],
+    )
 
 
 def effective_area_per_energy_and_fov(
-    selected_events, simulation_info, true_energy_bins, fov_offset_bins
+    selected_events,
+    simulation_info,
+    true_energy_axis: MapAxis,
+    fov_offset_axis: MapAxis,
 ):
     """
     Calculate effective area in bins of true energy and field of view offset.
@@ -64,24 +89,32 @@ def effective_area_per_energy_and_fov(
         - `true_source_fov_offset`
     simulation_info: pyirf.simulations.SimulatedEventsInfo
         The overall statistics of the simulated events
-    true_energy_bins: astropy.units.Quantity[energy]
-        The true energy bin edges in which to calculate effective area.
-    fov_offset_bins: astropy.units.Quantity[angle]
-        The field of view radial bin edges in which to calculate effective area.
+    true_energy_axis: MapAxis[energy]
+        The true energy axis in which to calculate effective area.
+    fov_offset_axis: MapAxis[angle]
+        The field of view radial offset axis in which to calculate effective area.
     """
     area = np.pi * simulation_info.max_impact ** 2
 
     hist_simulated = simulation_info.calculate_n_showers_per_energy_and_fov(
-        true_energy_bins, fov_offset_bins
+        true_energy_axis.edges,
+        fov_offset_axis.edges,
     )
 
     hist_selected, _, _ = np.histogram2d(
         selected_events["true_energy"].to_value(u.TeV),
         selected_events["true_source_fov_offset"].to_value(u.deg),
         bins=[
-            true_energy_bins.to_value(u.TeV),
-            fov_offset_bins.to_value(u.deg),
+            true_energy_axis.edges.to_value(u.TeV),
+            fov_offset_axis.edges.to_value(u.deg),
         ],
     )
 
-    return effective_area(hist_selected, hist_simulated, area)
+    aeff = effective_area(hist_selected, hist_simulated, area)
+    return EffectiveAreaTable2D(
+        axes=[
+            true_energy_axis,
+            fov_offset_axis,
+        ],
+        data=aeff,
+    )