1
|
|
|
# Licensed under a 3-clause BSD style license - see LICENSE.rst |
2
|
|
|
import pytest |
3
|
|
|
import numpy as np |
4
|
|
|
from numpy.testing import assert_allclose, assert_equal |
5
|
|
|
import astropy.units as u |
6
|
|
|
from gammapy.irf import EffectiveAreaTable2D |
7
|
|
|
from gammapy.maps import MapAxis |
8
|
|
|
from gammapy.utils.testing import ( |
9
|
|
|
assert_quantity_allclose, |
10
|
|
|
mpl_plot_check, |
11
|
|
|
requires_data, |
12
|
|
|
requires_dependency, |
13
|
|
|
) |
14
|
|
|
|
15
|
|
|
|
16
|
|
|
@pytest.fixture(scope="session") |
17
|
|
|
def aeff(): |
18
|
|
|
filename = "$GAMMAPY_DATA/hess-dl3-dr1/data/hess_dl3_dr1_obs_id_023523.fits.gz" |
19
|
|
|
return EffectiveAreaTable2D.read(filename, hdu="AEFF") |
20
|
|
|
|
21
|
|
|
|
22
|
|
|
@requires_data() |
23
|
|
|
def test_basic(aeff): |
24
|
|
|
assert aeff.axes["energy_true"].nbin == 96 |
25
|
|
|
assert aeff.axes["offset"].nbin == 6 |
26
|
|
|
assert aeff.data.shape == (96, 6) |
27
|
|
|
|
28
|
|
|
assert aeff.axes["energy_true"].unit == "TeV" |
29
|
|
|
assert aeff.axes["offset"].unit == "deg" |
30
|
|
|
assert aeff.unit == "m2" |
31
|
|
|
|
32
|
|
|
assert_quantity_allclose(aeff.meta["HI_THRES"], 100, rtol=1e-3) |
33
|
|
|
assert_quantity_allclose(aeff.meta["LO_THRES"], 0.870964, rtol=1e-3) |
34
|
|
|
|
35
|
|
|
test_val = aeff.evaluate(energy_true="14 TeV", offset="0.2 deg") |
36
|
|
|
assert_allclose(test_val.value, 683177.5, rtol=1e-3) |
37
|
|
|
|
38
|
|
|
|
39
|
|
|
def test_from_parametrization(): |
40
|
|
|
# Log center of this is 100 GeV |
41
|
|
|
area_ref = 1.65469579e07 * u.cm ** 2 |
42
|
|
|
|
43
|
|
|
axis = MapAxis.from_energy_edges([80, 125] * u.GeV, name="energy_true") |
44
|
|
|
area = EffectiveAreaTable2D.from_parametrization(axis, "HESS") |
45
|
|
|
|
46
|
|
|
assert_allclose(area.quantity, area_ref) |
47
|
|
|
assert area.unit == area_ref.unit |
48
|
|
|
|
49
|
|
|
# Log center of this is 0.1, 2 TeV |
50
|
|
|
area_ref = [1.65469579e07, 1.46451957e09] * u.cm * u.cm |
51
|
|
|
|
52
|
|
|
axis = MapAxis.from_energy_edges([0.08, 0.125, 32] * u.TeV, name="energy_true") |
53
|
|
|
area = EffectiveAreaTable2D.from_parametrization(axis, "HESS") |
54
|
|
|
assert_allclose(area.quantity[:, 0], area_ref) |
55
|
|
|
assert area.unit == area_ref.unit |
56
|
|
|
assert area.meta["TELESCOP"] == "HESS" |
57
|
|
|
|
58
|
|
|
|
59
|
|
|
@requires_dependency("matplotlib") |
60
|
|
|
@requires_data() |
61
|
|
|
def test_plot(aeff): |
62
|
|
|
with mpl_plot_check(): |
63
|
|
|
aeff.plot() |
64
|
|
|
|
65
|
|
|
with mpl_plot_check(): |
66
|
|
|
aeff.plot_energy_dependence() |
67
|
|
|
|
68
|
|
|
with mpl_plot_check(): |
69
|
|
|
aeff.plot_offset_dependence() |
70
|
|
|
|
71
|
|
|
|
72
|
|
|
def test_to_table(): |
73
|
|
|
energy_axis_true = MapAxis.from_energy_bounds( |
74
|
|
|
"1 TeV", "10 TeV", nbin=10, name="energy_true" |
75
|
|
|
) |
76
|
|
|
|
77
|
|
|
offset_axis = MapAxis.from_bounds(0, 1, nbin=4, name="offset", unit="deg") |
78
|
|
|
|
79
|
|
|
aeff = EffectiveAreaTable2D( |
80
|
|
|
axes=[energy_axis_true, offset_axis], data=1, unit="cm2" |
81
|
|
|
) |
82
|
|
|
hdu = aeff.to_table_hdu() |
83
|
|
|
assert_equal( |
84
|
|
|
hdu.data["ENERG_LO"][0], aeff.axes["energy_true"].edges[:-1].value |
85
|
|
|
) |
86
|
|
|
assert hdu.header["TUNIT1"] == aeff.axes["energy_true"].unit |
87
|
|
|
|
88
|
|
|
|
89
|
|
|
def test_wrong_axis_order(): |
90
|
|
|
energy_axis_true = MapAxis.from_energy_bounds( |
91
|
|
|
"1 TeV", "10 TeV", nbin=10, name="energy_true" |
92
|
|
|
) |
93
|
|
|
|
94
|
|
|
offset = np.linspace(0, 1, 4) * u.deg |
95
|
|
|
offset_axis = MapAxis.from_nodes(offset, name="offset") |
96
|
|
|
|
97
|
|
|
data = np.ones(shape=(offset_axis.nbin, energy_axis_true.nbin)) |
98
|
|
|
|
99
|
|
|
with pytest.raises(ValueError): |
100
|
|
|
EffectiveAreaTable2D( |
101
|
|
|
axes=[energy_axis_true, offset_axis], data=data, unit="cm2" |
102
|
|
|
) |
103
|
|
|
|
104
|
|
|
|
105
|
|
|
@requires_data("gammapy-data") |
106
|
|
|
def test_aeff2d_pointlike(): |
107
|
|
|
filename = "$GAMMAPY_DATA/joint-crab/dl3/magic/run_05029748_DL3.fits" |
108
|
|
|
|
109
|
|
|
aeff = EffectiveAreaTable2D.read(filename) |
110
|
|
|
hdu = aeff.to_table_hdu() |
111
|
|
|
|
112
|
|
|
assert aeff.is_pointlike |
113
|
|
|
assert hdu.header["HDUCLAS3"] == "POINT-LIKE" |
114
|
|
|
|