/
test_group_observations.py
53 lines (40 loc) · 1.92 KB
/
test_group_observations.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
from astropy.tests.helper import assert_quantity_allclose
from gammapy.data import ObservationTable
from gammapy.datasets import gammapy_extra
from gammapy.spectrum import SpectrumFit, SpectrumObservationGrouping, group_obs_table
obs_table_file = gammapy_extra.filename(
'datasets/hess-crab4_pha/observation_table.fits')
obs_table = ObservationTable.read(obs_table_file)
fit = SpectrumFit.from_observation_table(obs_table)
fit.model = 'PL'
fit.energy_threshold_low = '1 TeV'
fit.energy_threshold_high = '10 TeV'
fit.run(method='sherpa')
#Use each obs in one group
obs_table1 = group_obs_table(obs_table, eff_range=[90, 95], n_eff_bin=5)
obs_table1.write('grouped_table1_debug.fits', overwrite=True)
grouping = SpectrumObservationGrouping(obs_table1)
grouping.run()
fit_band2 = SpectrumFit.from_observation_table(grouping.stacked_obs_table)
fit_band2.model = 'PL'
fit_band2.energy_threshold_low = '1 TeV'
fit_band2.energy_threshold_high = '10 TeV'
fit_band2.run(method='sherpa')
assert_quantity_allclose(fit.result.parameters["index"],
fit_band2.result.parameters["index"], rtol=1e-5)
assert_quantity_allclose(fit.result.parameters["norm"],
fit_band2.result.parameters["norm"], rtol=1e-5)
# Put all runs in one group
obs_table2 = group_obs_table(obs_table, n_eff_bin=1, n_off_bin=1, n_zen_bin=1)
obs_table2.write('grouped_table2_debug.fits', overwrite=True)
grouping = SpectrumObservationGrouping(obs_table2)
grouping.run()
fit_band3 = SpectrumFit.from_observation_table(grouping.stacked_obs_table)
fit_band3.model = 'PL'
fit_band3.energy_threshold_low = '100 GeV'
fit_band3.energy_threshold_high = '10 TeV'
fit_band3.run(method='sherpa')
assert_quantity_allclose(fit.result.parameters["index"],
fit_band3.result.parameters["index"], rtol=1e0)
assert_quantity_allclose(fit_band3.result.parameters["norm"],
fit.result.parameters["norm"], rtol=1e0)