%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

"""
Explanation: CTA AEFF interpolation
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).
End of explanation
"""


filename = '1dc/1dc/caldb/data/cta/1dc/bcf/North_z20_50h/irf_file.fits'
aeff = EffectiveAreaTable2D.read(filename)

# Just to compare, what the raw IRF BINTABLE HDU contains
table = Table.read(filename, hdu='EFFECTIVE AREA')
table

print(aeff)

energy_axis = aeff.data.axes[0]
print(energy_axis)
energy_axis.nodes.value

offset_axis = aeff.data.axes[1]
print(offset_axis)
print(offset_axis.nodes.value)

"""
Explanation: Bins and nodes
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.
End of explanation
"""


aeff.peek()

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()

"""
Explanation: Peek
Let's have a quick look at AEFF, and especially at the extrapolation to offset = 0 deg.
Everything looks as expected: the AEFF values at offset = 0 deg are the linear extrapolation of the values at the two nearest nodes from offset = 0.5 deg and 1.5 deg.
End of explanation
"""


# 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())

"""
Explanation: Evaluate at nodes
We're using bilinear interpolation in energy and offset, so evaluating at nodes should give the exact same values as the data array.
End of explanation
"""


energy = np.logspace(-2, 2, 300) * u.TeV
val = aeff.data.evaluate(energy=energy, offset=0.5*u.deg)

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)

plt.

"""
Explanation: Extrapolation
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.
End of explanation
"""


