diff --git a/examples/calculate_eventdisplay_irfs.py b/examples/calculate_eventdisplay_irfs.py index a4dd0ad54..a756017aa 100644 --- a/examples/calculate_eventdisplay_irfs.py +++ b/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,7 +142,7 @@ 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, @@ -147,11 +150,11 @@ def main(): ) # 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")) diff --git a/pyirf/gammapy.py b/pyirf/gammapy.py index 3eae90fe8..9c9a94979 100644 --- a/pyirf/gammapy.py +++ b/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, ) diff --git a/pyirf/irf/effective_area.py b/pyirf/irf/effective_area.py index daa13e767..64a7c7489 100644 --- a/pyirf/irf/effective_area.py +++ b/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, + ) diff --git a/pyirf/irf/energy_dispersion.py b/pyirf/irf/energy_dispersion.py index 9bb203999..df69257d1 100644 --- a/pyirf/irf/energy_dispersion.py +++ b/pyirf/irf/energy_dispersion.py @@ -1,12 +1,13 @@ import warnings import numpy as np import astropy.units as u -from ..binning import resample_histogram1d + +from gammapy.irf import EnergyDispersion2D +from gammapy.maps import MapAxis __all__ = [ "energy_dispersion", - "energy_dispersion_to_migration" ] @@ -28,7 +29,10 @@ def _normalize_hist(hist): def energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins, + selected_events, + true_energy_axis: MapAxis, + fov_offset_axis: MapAxis, + migration_axis: MapAxis, ): """ Calculate energy dispersion for the given DL2 event list. @@ -41,12 +45,12 @@ def energy_dispersion( selected_events: astropy.table.QTable Table of the DL2 events. Required columns: ``reco_energy``, ``true_energy``, ``true_source_fov_offset``. - true_energy_bins: astropy.units.Quantity[energy] + true_energy_aix: MapAxis[energy] Bin edges in true energy - migration_bins: astropy.units.Quantity[energy] + migration_axis: MapAxis[dimensionless] Bin edges in relative deviation, recommended range: [0.2, 5] - fov_offset_bins: astropy.units.Quantity[angle] - Bin edges in the field of view offset. + fov_offset_axis: MapAxis[angle] + Field of view offset axis. For Point-Like IRFs, only giving a single bin is appropriate. Returns @@ -68,91 +72,19 @@ def energy_dispersion( ] ), bins=[ - true_energy_bins.to_value(u.TeV), - migration_bins, - fov_offset_bins.to_value(u.deg), + true_energy_axis.edges.to_value(u.TeV), + migration_axis.edges.to_value(u.one), + fov_offset_axis.edges.to_value(u.deg), ], ) - n_events_per_energy = energy_dispersion.sum(axis=1) energy_dispersion = _normalize_hist(energy_dispersion) - return energy_dispersion - - -def energy_dispersion_to_migration( - dispersion_matrix, - disp_true_energy_edges, - disp_migration_edges, - new_true_energy_edges, - new_reco_energy_edges, -): - """ - Construct a energy migration matrix from a energy - dispersion matrix. - Depending on the new energy ranges, the sum over the first axis - can be smaller than 1. - The new true energy bins need to be a subset of the old range, - extrapolation is not supported. - New reconstruction bins outside of the old migration range are filled with - zeros. - - Parameters - ---------- - dispersion_matrix: numpy.ndarray - Energy dispersion_matrix - disp_true_energy_edges: astropy.units.Quantity[energy] - True energy edges matching the first dimension of the dispersion matrix - disp_migration_edges: numpy.ndarray - Migration edges matching the second dimension of the dispersion matrix - new_true_energy_edges: astropy.units.Quantity[energy] - True energy edges matching the first dimension of the output - new_reco_energy_edges: astropy.units.Quantity[energy] - Reco energy edges matching the second dimension of the output - - Returns: - -------- - migration_matrix: numpy.ndarray - Three-dimensional energy migration matrix. The third dimension - equals the fov offset dimension of the energy dispersion matrix. - """ - - migration_matrix = np.zeros(( - len(new_true_energy_edges) - 1, - len(new_reco_energy_edges) - 1, - dispersion_matrix.shape[2], - )) - - true_energy_interpolation = resample_histogram1d( - dispersion_matrix, - disp_true_energy_edges, - new_true_energy_edges, - axis=0, + return EnergyDispersion2D( + axes=[ + true_energy_axis, + migration_axis, + fov_offset_axis, + ], + data=energy_dispersion ) - - norm = np.sum(true_energy_interpolation, axis=1, keepdims=True) - norm[norm == 0] = 1 - true_energy_interpolation /= norm - - for idx, e_true in enumerate( - (new_true_energy_edges[1:] + new_true_energy_edges[:-1]) / 2 - ): - - # get migration for the new true energy bin - e_true_dispersion = true_energy_interpolation[idx] - - with warnings.catch_warnings(): - # silence inf/inf division warning - warnings.filterwarnings('ignore', 'invalid value encountered in true_divide') - interpolation_edges = new_reco_energy_edges / e_true - - y = resample_histogram1d( - e_true_dispersion, - disp_migration_edges, - interpolation_edges, - axis=0, - ) - - migration_matrix[idx, :, :] = y - - return migration_matrix diff --git a/pyirf/irf/psf.py b/pyirf/irf/psf.py index eaa3c6114..5f90f7e29 100644 --- a/pyirf/irf/psf.py +++ b/pyirf/irf/psf.py @@ -1,10 +1,17 @@ import numpy as np import astropy.units as u +from gammapy.irf import PSF3D +from gammapy.maps import MapAxis from ..utils import cone_solid_angle -def psf_table(events, true_energy_bins, source_offset_bins, fov_offset_bins): +def psf_table( + events, + true_energy_axis: MapAxis, + source_offset_axis: MapAxis, + fov_offset_axis: MapAxis +): """ Calculate the table based PSF (radially symmetrical bins around the true source) """ @@ -20,14 +27,22 @@ def psf_table(events, true_energy_bins, source_offset_bins, fov_offset_bins): hist, _ = np.histogramdd( array, [ - true_energy_bins.to_value(u.TeV), - fov_offset_bins.to_value(u.deg), - source_offset_bins.to_value(u.deg), + true_energy_axis.edges.to_value(u.TeV), + fov_offset_axis.edges.to_value(u.deg), + source_offset_axis.edges.to_value(u.deg), ], ) - psf = _normalize_psf(hist, source_offset_bins) - return psf + psf = _normalize_psf(hist, source_offset_axis.edges) + + return PSF3D( + axes=[ + true_energy_axis, + fov_offset_axis, + source_offset_axis, + ], + data=psf + ) def _normalize_psf(hist, source_offset_bins): diff --git a/pyirf/irf/tests/test_effective_area.py b/pyirf/irf/tests/test_effective_area.py index e08279eb1..d299740d0 100644 --- a/pyirf/irf/tests/test_effective_area.py +++ b/pyirf/irf/tests/test_effective_area.py @@ -1,6 +1,7 @@ import astropy.units as u import numpy as np from astropy.table import QTable +from gammapy.maps import MapAxis def test_effective_area(): @@ -20,7 +21,8 @@ def test_effective_area_per_energy(): from pyirf.irf import effective_area_per_energy from pyirf.simulations import SimulatedEventsInfo - true_energy_bins = [0.1, 1.0, 10.0] * u.TeV + true_energy_axis = MapAxis.from_edges([0.1, 1.0, 10.0], unit=u.TeV, name='energy_true') + fov_offset_axis = MapAxis.from_edges([0, 0.1], unit=u.deg, name='offset') selected_events = QTable( { "true_energy": np.append(np.full(1000, 0.5), np.full(10, 5)), @@ -30,30 +32,35 @@ def test_effective_area_per_energy(): # this should give 100000 events in the first bin and 10000 in the second simulation_info = SimulatedEventsInfo( n_showers=110000, - energy_min=true_energy_bins[0], - energy_max=true_energy_bins[-1], + energy_min=true_energy_axis.edges[0], + energy_max=true_energy_axis.edges[-1], max_impact=100 / np.sqrt(np.pi) * u.m, # this should give a nice round area spectral_index=-2, viewcone=0 * u.deg, ) - area = effective_area_per_energy(selected_events, simulation_info, true_energy_bins) + area = effective_area_per_energy( + selected_events, + simulation_info, + true_energy_axis, + fov_offset_axis, + ) - assert area.shape == (len(true_energy_bins) - 1,) - assert area.unit == u.m ** 2 - assert u.allclose(area, [100, 10] * u.m ** 2) + assert u.allclose(area.quantity, [[100], [10]] * u.m ** 2) def test_effective_area_energy_fov(): from pyirf.irf import effective_area_per_energy_and_fov from pyirf.simulations import SimulatedEventsInfo - true_energy_bins = [0.1, 1.0, 10.0] * u.TeV + true_energy_axis = MapAxis.from_edges([0.1, 1.0, 10.0], unit=u.TeV, name='energy_true') # choose edges so that half are in each bin in fov - fov_offset_bins = [0, np.arccos(0.98), np.arccos(0.96)] * u.rad - center_1, center_2 = 0.5 * (fov_offset_bins[:-1] + fov_offset_bins[1:]).to_value( - u.deg + fov_offset_axis = MapAxis.from_edges( + [0, np.arccos(0.98), np.arccos(0.96)], + unit=u.rad, + name='offset' ) + center_1, center_2 = fov_offset_axis.center.to_value(u.deg) selected_events = QTable( { @@ -76,18 +83,19 @@ def test_effective_area_energy_fov(): # this should give 100000 events in the first bin and 10000 in the second simulation_info = SimulatedEventsInfo( n_showers=110000, - energy_min=true_energy_bins[0], - energy_max=true_energy_bins[-1], + energy_min=true_energy_axis.edges[0], + energy_max=true_energy_axis.edges[-1], max_impact=100 / np.sqrt(np.pi) * u.m, # this should give a nice round area spectral_index=-2, - viewcone=fov_offset_bins[-1], + viewcone=fov_offset_axis.edges[-1], ) - area = effective_area_per_energy_and_fov( - selected_events, simulation_info, true_energy_bins, fov_offset_bins + area_table = effective_area_per_energy_and_fov( + selected_events, simulation_info, true_energy_axis, fov_offset_axis ) + area = area_table.quantity - assert area.shape == (len(true_energy_bins) - 1, len(fov_offset_bins) - 1) + assert area.shape == (true_energy_axis.nbin, fov_offset_axis.nbin) assert area.unit == u.m ** 2 assert u.allclose(area[:, 0], [200, 20] * u.m ** 2) assert u.allclose(area[:, 1], [100, 10] * u.m ** 2) diff --git a/pyirf/irf/tests/test_energy_dispersion.py b/pyirf/irf/tests/test_energy_dispersion.py index b7dba84a5..1887df373 100644 --- a/pyirf/irf/tests/test_energy_dispersion.py +++ b/pyirf/irf/tests/test_energy_dispersion.py @@ -1,5 +1,6 @@ import astropy.units as u import numpy as np +from gammapy.maps import MapAxis from astropy.table import QTable @@ -41,19 +42,20 @@ def test_energy_dispersion(): } ) - true_energy_bins = np.array([0.1, 1.0, 10.0, 100]) * u.TeV - fov_offset_bins = np.array([0, 1, 2]) * u.deg - migration_bins = np.linspace(0, 2, 1001) + true_energy_axis = MapAxis.from_edges([0.1, 1.0, 10.0, 100], unit=u.TeV, name='energy_true') + fov_offset_axis = MapAxis.from_edges([0, 1, 2], unit=u.deg, name='offset') + migration_axis = MapAxis.from_bounds(0.2, 5, 1000, interp='log', name='migra') - result = energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins + edisp2d = energy_dispersion( + selected_events, true_energy_axis, fov_offset_axis, migration_axis ) - assert result.shape == (3, 1000, 2) - assert np.isclose(result.sum(), 6.0) + edisp = edisp2d.quantity + assert edisp.shape == (3, 1000, 2) + assert np.isclose(edisp.sum(), 6.0) - cumulative_sum = np.cumsum(result, axis=1) - bin_centers = 0.5 * (migration_bins[1:] + migration_bins[:-1]) + cumulative_sum = np.cumsum(edisp, axis=1) + bin_centers = migration_axis.center assert np.isclose( TRUE_SIGMA_1, ( @@ -81,125 +83,3 @@ def test_energy_dispersion(): / 2, rtol=0.1, ) - - -def test_energy_dispersion_to_migration(): - from pyirf.irf import energy_dispersion - from pyirf.irf.energy_dispersion import energy_dispersion_to_migration - - np.random.seed(0) - N = 10000 - true_energy_bins = 10**np.arange(np.log10(0.2), np.log10(200), 1/10) * u.TeV - - fov_offset_bins = np.array([0, 1, 2]) * u.deg - migration_bins = np.linspace(0, 2, 101) - - true_energy = np.random.uniform( - true_energy_bins[0].value, - true_energy_bins[-1].value, - size=N - ) * u.TeV - reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) - - selected_events = QTable( - { - "reco_energy": reco_energy, - "true_energy": true_energy, - "true_source_fov_offset": np.concatenate( - [ - np.full(N // 2, 0.2), - np.full(N // 2, 1.5), - ] - ) - * u.deg, - } - ) - - dispersion_matrix = energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins - ) - - # migration matrix selecting a limited energy band with different binsizes - new_true_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV - new_reco_energy_bins = 10**np.arange(np.log10(2), np.log10(20), 1/5) * u.TeV - migration_matrix = energy_dispersion_to_migration( - dispersion_matrix, - true_energy_bins, - migration_bins, - new_true_energy_bins, - new_reco_energy_bins, - - ) - - # test dimension - assert migration_matrix.shape[0] == len(new_true_energy_bins) - 1 - assert migration_matrix.shape[1] == len(new_reco_energy_bins) - 1 - assert migration_matrix.shape[2] == dispersion_matrix.shape[2] - - # test that all migrations are included for central energies - assert np.isclose(migration_matrix.sum(axis=1).max(), 1, rtol=0.01) - - # test that migrations dont always sum to 1 (since some energies are - # not included in the matrix) - assert migration_matrix.sum() < (len(new_true_energy_bins) - 1) * (len(fov_offset_bins) - 1) - assert np.all(np.isfinite(migration_matrix)) - - -def test_overflow_bins_migration_matrix(): - from pyirf.irf import energy_dispersion - from pyirf.irf.energy_dispersion import energy_dispersion_to_migration - from pyirf.binning import add_overflow_bins - - np.random.seed(0) - N = 10000 - - # add under/overflow bins - bins = 10 ** np.arange( - np.log10(0.2), - np.log10(200), - 1 / 10, - ) * u.TeV - true_energy_bins = add_overflow_bins(bins, positive=False) - - fov_offset_bins = np.array([0, 1, 2]) * u.deg - migration_bins = np.linspace(0, 2, 101) - - true_energy = np.random.uniform( - true_energy_bins[1].value, - true_energy_bins[-2].value, - size=N - ) * u.TeV - reco_energy = true_energy * np.random.uniform(0.5, 1.5, size=N) - - selected_events = QTable( - { - "reco_energy": reco_energy, - "true_energy": true_energy, - "true_source_fov_offset": np.concatenate( - [ - np.full(N // 2, 0.2), - np.full(N // 2, 1.5), - ] - ) - * u.deg, - } - ) - - dispersion_matrix = energy_dispersion( - selected_events, true_energy_bins, fov_offset_bins, migration_bins - ) - - migration_matrix = energy_dispersion_to_migration( - dispersion_matrix, - true_energy_bins, - migration_bins, - true_energy_bins, - true_energy_bins, - ) - - # migration into underflow bin is not empty - assert np.sum(migration_matrix[:, 0, :]) > 0 - # migration from underflow bin is empty - assert np.sum(migration_matrix[0, :, :]) == 0 - - assert np.all(np.isfinite(migration_matrix)) diff --git a/pyirf/irf/tests/test_psf.py b/pyirf/irf/tests/test_psf.py index f0384dda6..e7ddf826e 100644 --- a/pyirf/irf/tests/test_psf.py +++ b/pyirf/irf/tests/test_psf.py @@ -2,6 +2,8 @@ from astropy.table import QTable import numpy as np +from gammapy.maps import MapAxis + def test_psf(): from pyirf.irf import psf_table @@ -25,25 +27,25 @@ def test_psf(): } ) - energy_bins = [0, 1.5, 3] * u.TeV - fov_bins = [0, 1] * u.deg - source_bins = np.linspace(0, 1, 201) * u.deg + energy_axis = MapAxis.from_edges([0, 1.5, 3], unit=u.TeV, name='energy_true') + fov_axis = MapAxis.from_edges([0, 1], unit=u.deg, name='offset') + source_axis = MapAxis.from_bounds(0, 1, 200, unit=u.deg, name='rad') # We return a table with one row as needed for gadf - psf = psf_table(events, energy_bins, source_bins, fov_bins) + psf = psf_table(events, energy_axis, source_axis, fov_axis) # 2 energy bins, 1 fov bin, 200 source distance bins - assert psf.shape == (2, 1, 200) - assert psf.unit == u.Unit("sr-1") + assert psf.quantity.shape == (2, 1, 200) + assert psf.quantity.unit == u.Unit("sr-1") # check that psf is normalized - bin_solid_angle = np.diff(cone_solid_angle(source_bins)) - assert np.allclose(np.sum(psf * bin_solid_angle, axis=2), 1.0) + bin_solid_angle = np.diff(cone_solid_angle(source_axis.edges)) + assert np.allclose(np.sum(psf.quantity * bin_solid_angle, axis=2), 1.0) - cumulated = np.cumsum(psf * bin_solid_angle, axis=2) + cumulated = np.cumsum(psf.quantity * bin_solid_angle, axis=2) # first energy and only fov bin - bin_centers = 0.5 * (source_bins[1:] + source_bins[:-1]) + bin_centers = source_axis.center assert u.isclose( bin_centers[np.where(cumulated[0, 0, :] >= 0.68)[0][0]], TRUE_SIGMA_1 * u.deg, diff --git a/setup.py b/setup.py index 3c35c74ae..67c08ead5 100644 --- a/setup.py +++ b/setup.py @@ -2,7 +2,6 @@ import os -gammapy = "gammapy~=0.19" extras_require = { "docs": [ @@ -15,19 +14,14 @@ "awkward1", "notebook", "tables", - gammapy, ], "tests": [ "pytest", "pytest-cov", - gammapy, "ogadf-schema~=0.2.3", "uproot~=4.0", "awkward~=1.0", ], - "gammapy": [ - gammapy, - ] } extras_require["all"] = list(set(extras_require["tests"] + extras_require["docs"])) @@ -42,6 +36,7 @@ "scipy", "tqdm", "setuptools_scm", + "gammapy~=0.19", ], include_package_data=True, extras_require=extras_require,