There was a report of possible bugs with CTA DC-1 AEFF interpolation via email.
In this notebook we have a quick look.
As far as I can see, everything is as expected and as good as it can be (given the somewhat noisy IRFs from CTA, and assuming we're not introducing smoothing in Gammapy for now).
In [1]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import astropy.units as u
from astropy.table import Table
from gammapy.irf import EffectiveAreaTable2D
Everything looks as expected to me.
Note that CTA IRFs currently are produced from diffuse photons and IRFs computed in offset bins of 0 to 1 deg, 1 to 2 deg and so on. In Gammapy, we choose to put the node for the interpolation at the bin center in offset, i.e. the first node is at 0.5 deg, the second at 1.5 deg and so on.
In [2]:
filename = '1dc/1dc/caldb/data/cta/1dc/bcf/North_z20_50h/irf_file.fits'
aeff = EffectiveAreaTable2D.read(filename)
In [3]:
# Just to compare, what the raw IRF BINTABLE HDU contains
table = Table.read(filename, hdu='EFFECTIVE AREA')
table
Out[3]:
In [4]:
print(aeff)
In [5]:
energy_axis = aeff.data.axes[0]
print(energy_axis)
energy_axis.nodes.value
Out[5]:
In [6]:
offset_axis = aeff.data.axes[1]
print(offset_axis)
print(offset_axis.nodes.value)
In [7]:
aeff.peek()
In [8]:
for offset in np.arange(0, 2, 0.5):
val = aeff.data.evaluate(offset=offset*u.deg)
plt.plot(energy_axis.nodes.value, val, label=offset)
plt.legend()
Out[8]:
In [9]:
# If no energy and offset is passed, the interpolator is called at the nodes by default
val = aeff.data.evaluate()
val2 = aeff.data.data
diff = val - val2
print(diff.max())
When extrapolating to energies below the lowest node at 0.014 TeV, negative effective area values do occur.
Running analyses that include energies so low will probably cause issues. So for now to get correct results: just put a minumum energy > 0.015 TeV.
We can discuss if we want to clip AEFF to values >= 0, or pad the data array with rows of zeros so that the interpolator does the right thing directly.
In [10]:
energy = np.logspace(-2, 2, 300) * u.TeV
val = aeff.data.evaluate(energy=energy, offset=0.5*u.deg)
In [11]:
plt.plot(energy, val, 'o')
plt.plot(energy_axis.nodes.value, aeff.data.data[:, 0], 'o')
plt.xlim(0.005, 0.03)
plt.ylim(-10000, 10000)
Out[11]:
In [12]:
plt.