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, EnergyDispersion2D

In [2]:
filename = '1dc/caldb/data/cta/1dc/bcf/North_z20_50h/irf_file.fits'
aeff =
edisp =,hdu='ENERGY DISPERSION')

In [3]:

NDDataArray summary info
e_true         : size =   500, min =  0.005 TeV, max = 495.450 TeV
migra          : size =   300, min =  0.005, max =  2.995
offset         : size =     6, min =  0.500 deg, max =  5.500 deg
Data           : size = 900000, min =  0.000, max = 16039.922


In [4]:

/Users/terrier/Code/anaconda2/lib/python2.7/site-packages/astropy/units/ RuntimeWarning: invalid value encountered in true_divide
  return super(Quantity, self).__truediv__(other)
/Users/terrier/Code/anaconda2/lib/python2.7/site-packages/matplotlib/ UserWarning: Power-law scaling on negative values is ill-defined, clamping to 0.
  warnings.warn("Power-law scaling on negative values is "
/Users/terrier/Code/anaconda2/lib/python2.7/site-packages/matplotlib/ RuntimeWarning: invalid value encountered in power
  np.power(resdat, gamma, resdat)

Look at effect of extrapolation towards 0. We look at bin number 260 in e_true, i.e. 2 TeV.

Negative values appear in the distribution of migration

In [5]:
for offset in (0.0*u.deg, 0.25*u.deg, 0.5*u.deg):
    val = = offset,

<matplotlib.legend.Legend at 0x1113b0890>

Now look at extracted 1D edisp (rmf)

In [6]:
rmfOK = edisp.to_energy_dispersion(offset=0.5*u.deg)
rmfbad = edisp.to_energy_dispersion(offset=0.0*u.deg)

At 0.5° it is OK

In [7]:

<matplotlib.axes._subplots.AxesSubplot at 0x1113cfa10>

At 0° this looks bad.

In [8]:

<matplotlib.axes._subplots.AxesSubplot at 0x116b6e150>

Let's look at 'predicted' number of counts.

They can get negative in the extrapolated case.

The easiest workaround is to force offset extraction at 0.5° in case it is lower.

In [9]:
model = rmfbad.e_true.log_center()**-1
plt.semilogx(rmfbad.e_reco.log_center(),rmfOK.apply(model),label='0.5 deg')

<matplotlib.legend.Legend at 0x10417c1d0>

Proposed patch to remove negative predicted counts

Simply suppress negative values produced by EnergyDispersion2D.get_response

In [10]:
from import EnergyBounds, Energy
def get_response(self, offset, e_true, e_reco=None, migra_step=5e-3):
        """Detector response R(Delta E_reco, E_true)

        Probability to reconstruct a given true energy in a given reconstructed
        energy band. In each reco bin, you integrate with a riemann sum over
        the default migra bin of your analysis.

        e_true : ``
            True energy
        e_reco : ``, None
            Reconstructed energy axis
        offset : `~astropy.coordinates.Angle`
        migra_step : float
            Integration step in migration

        rv : `~numpy.ndarray`
            Redistribution vector
        e_true = Energy(e_true)

        # Default: e_reco nodes = migra nodes * e_true nodes
        if e_reco is None:
            e_reco = EnergyBounds.from_lower_and_upper_bounds(
                self.migra.lo * e_true, self.migra.hi * e_true)
            migra = self.migra.nodes
        # Translate given e_reco binning to migra at bin center
            e_reco = EnergyBounds(e_reco)
            center = e_reco.log_centers
            migra = center / e_true

        # migration value of e_reco bounds
        migra_e_reco = e_reco / e_true

        # Define a vector of migration with mig_step step
        mrec_min = self.migra.lo[0]
        mrec_max = self.migra.hi[-1]
        mig_array = np.arange(mrec_min, mrec_max, migra_step)

        # Compute energy dispersion probability dP/dm for each element of migration array
        vals =, e_true=e_true, migra=mig_array)

        # Remove negative values likely caused by extrapolation if IRF is not defined properly
        vals[np.where(vals<0.)] = 0.0
        # Compute normalized cumulative sum to prepare integration
        tmp = np.nan_to_num(np.cumsum(vals) / np.sum(vals))

        # Determine positions (bin indices) of e_reco bounds in migration array
        pos_mig = np.digitize(migra_e_reco, mig_array) - 1
        # We ensure that no negative values are found
        pos_mig = np.maximum(pos_mig, 0)

        # We compute the difference between 2 successive bounds in e_reco
        # to get integral over reco energy bin
        integral = np.diff(tmp[pos_mig])

        return integral
EnergyDispersion2D.get_response = get_response

In [11]:
rmfOK = edisp.to_energy_dispersion(offset=0.5*u.deg)
rmfbad = edisp.to_energy_dispersion(offset=0.0*u.deg)

In [12]:
model = rmfbad.e_true.log_center()**-1
plt.semilogx(rmfbad.e_reco.log_center(),rmfOK.apply(model),label='0.5 deg')

<matplotlib.legend.Legend at 0x11225ce50>

In [ ]: