Noisy CTA north IRFs

Julien / Bruno mentioned via email that there are problems with CTA 1DC analyses with Gammapy, because the north IRFs have noisy AEFF and EDISP. This is a very quick look ...


In [3]:
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

In [4]:
# import os
# os.environ["CTADATA"] = "/Users/julien/Documents/WorkingDir/Tools/python/cta_1dc/1dc/"

In [5]:
import astropy.units as u 
from gammapy.data import DataStore

# Load data
data_store = DataStore.from_dir('$CTADATA/index/agn')
obs_list = data_store.obs_list([510004])

# edisp
ed = obs_list[0].edisp
emin = 0.03 * u.TeV
emax = 20 * u.TeV
e_reco =  np.logspace(np.log10(emin.value),np.log10(emax.value),40).tolist() * u.TeV
e_true=np.logspace(-1.5, 2, 500) * u.TeV
rmf = ed.to_energy_dispersion(0.*u.deg,
                              e_true=e_true,
                              e_reco=e_reco)

rmf.plot_bias()


/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py:634: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
Out[5]:
<matplotlib.axes._subplots.AxesSubplot at 0x1079f3828>

In [6]:
obs = obs_list[0]

In [9]:
obs.events.table.meta


Out[9]:
OrderedDict([('EXTNAME', 'EVENTS'),
             ('DSTYP1', 'TIME'),
             ('DSUNI1', 's'),
             ('DSVAL1', 'TABLE'),
             ('DSREF1', ':GTI'),
             ('DSTYP2', 'ENERGY'),
             ('DSUNI2', 'TeV'),
             ('DSVAL2', '0.03:50'),
             ('DSTYP3', 'POS(RA,DEC)'),
             ('DSUNI3', 'deg'),
             ('DSVAL3', 'CIRCLE(166.1196,38.2073,5)'),
             ('NDSKEYS', 3),
             ('NMCIDS', 3),
             ('MID00001', 1832),
             ('MMN00001', '2FHL_J1104.4+3812'),
             ('MID00002', 2006),
             ('MMN00002', '1FHL_J1100.6+4018'),
             ('MID00003', 1),
             ('MMN00003', 'Background model'),
             ('CREATOR', 'ctobssim (1.4.0)'),
             ('TELESCOP', 'CTA'),
             ('OBS_ID', 510004),
             ('DATE_OBS', '2021-01-01'),
             ('TIME_OBS', '14:06:51'),
             ('DATE_END', '2021-01-01'),
             ('TIME_END', '14:36:51'),
             ('TSTART', 662782080.0),
             ('TSTOP', 662783872.0),
             ('MJDREFI', 51544),
             ('MJDREFF', 0.5),
             ('TIMEUNIT', 's'),
             ('TIMESYS', 'TT'),
             ('TIMEREF', 'LOCAL'),
             ('TELAPSE', 1800.0),
             ('ONTIME', 1800.0),
             ('LIVETIME', 1764.0),
             ('DEADC', 0.98000001907),
             ('TIMEDEL', 1),
             ('OBJECT', 'AGN monitoring'),
             ('RA_OBJ', 0),
             ('DEC_OBJ', 0),
             ('RA_PNT', 166.11959839),
             ('DEC_PNT', 38.207298279),
             ('ALT_PNT', 90.0),
             ('AZ_PNT', 0),
             ('RADECSYS', 'FK5'),
             ('EQUINOX', 2000.0),
             ('CONV_DEP', 0),
             ('CONV_RA', 0),
             ('CONV_DEC', 0),
             ('OBSERVER', 'CTA Consortium'),
             ('N_TELS', 0),
             ('TELLIST', 'Baseline'),
             ('GEOLAT', 28.7619),
             ('GEOLON', 17.89),
             ('ALTITUDE', 2.2),
             ('EUNIT', 'TeV'),
             ('EVTVER', 'draft1'),
             ('CALDB', '1dc'),
             ('IRF', 'North_z20_50h')])

In [10]:
obs.aeff.peek()



In [12]:
obs.aeff.data.data.min()


Out[12]:
$0 \; \mathrm{m^{2}}$

In [13]:
obs.edisp.data.data.min()


Out[13]:
$0 \; \mathrm{}$

In [15]:
obs.edisp.plot_migration()


Out[15]:
<matplotlib.axes._subplots.AxesSubplot at 0x10847f940>

In [16]:
obs.edisp.peek()


/Users/deil/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py:634: RuntimeWarning: invalid value encountered in true_divide
  result = super().__array_ufunc__(function, method, *arrays, **kwargs)
/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/colors.py:1238: UserWarning: Power-law scaling on negative values is ill-defined, clamping to 0.
  warnings.warn("Power-law scaling on negative values is "
---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/formatters.py in __call__(self, obj)
    339                 pass
    340             else:
--> 341                 return printer(obj)
    342             # Finally look for special method names
    343             method = get_real_method(obj, self.print_method)

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/pylabtools.py in <lambda>(fig)
    236 
    237     if 'png' in formats:
--> 238         png_formatter.for_type(Figure, lambda fig: print_figure(fig, 'png', **kwargs))
    239     if 'retina' in formats or 'png2x' in formats:
    240         png_formatter.for_type(Figure, lambda fig: retina_figure(fig, **kwargs))

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/IPython/core/pylabtools.py in print_figure(fig, fmt, bbox_inches, **kwargs)
    120 
    121     bytes_io = BytesIO()
--> 122     fig.canvas.print_figure(bytes_io, **kw)
    123     data = bytes_io.getvalue()
    124     if fmt == 'svg':

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backend_bases.py in print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, **kwargs)
   2206                     orientation=orientation,
   2207                     dryrun=True,
-> 2208                     **kwargs)
   2209                 renderer = self.figure._cachedRenderer
   2210                 bbox_inches = self.figure.get_tightbbox(renderer)

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py in print_png(self, filename_or_obj, *args, **kwargs)
    505 
    506     def print_png(self, filename_or_obj, *args, **kwargs):
--> 507         FigureCanvasAgg.draw(self)
    508         renderer = self.get_renderer()
    509         original_dpi = renderer.dpi

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/backends/backend_agg.py in draw(self)
    428             if toolbar:
    429                 toolbar.set_cursor(cursors.WAIT)
--> 430             self.figure.draw(self.renderer)
    431         finally:
    432             if toolbar:

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/figure.py in draw(self, renderer)
   1293 
   1294             mimage._draw_list_compositing_images(
-> 1295                 renderer, self, artists, self.suppressComposite)
   1296 
   1297             renderer.close_group('figure')

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    136     if not_composite or not has_images:
    137         for a in artists:
--> 138             a.draw(renderer)
    139     else:
    140         # Composite any adjacent images together

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/axes/_base.py in draw(self, renderer, inframe)
   2397             renderer.stop_rasterizing()
   2398 
-> 2399         mimage._draw_list_compositing_images(renderer, self, artists)
   2400 
   2401         renderer.close_group('axes')

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/image.py in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
    136     if not_composite or not has_images:
    137         for a in artists:
--> 138             a.draw(renderer)
    139     else:
    140         # Composite any adjacent images together

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/artist.py in draw_wrapper(artist, renderer, *args, **kwargs)
     53                 renderer.start_filter()
     54 
---> 55             return draw(artist, renderer, *args, **kwargs)
     56         finally:
     57             if artist.get_agg_filter() is not None:

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/collections.py in draw(self, renderer)
   1871                 offsets = np.column_stack([xs, ys])
   1872 
-> 1873         self.update_scalarmappable()
   1874 
   1875         if not transform.is_affine:

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/collections.py in update_scalarmappable(self)
    740             return
    741         if self._is_filled:
--> 742             self._facecolors = self.to_rgba(self._A, self._alpha)
    743         elif self._is_stroked:
    744             self._edgecolors = self.to_rgba(self._A, self._alpha)

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/cm.py in to_rgba(self, x, alpha, bytes, norm)
    273         x = ma.asarray(x)
    274         if norm:
--> 275             x = self.norm(x)
    276         rgba = self.cmap(x, alpha=alpha, bytes=bytes)
    277         return rgba

/opt/local/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/matplotlib/colors.py in __call__(self, value, clip)
   1198             resdat -= vmin
   1199             np.power(resdat, gamma, resdat)
-> 1200             resdat /= (vmax - vmin) ** gamma
   1201 
   1202             result = np.ma.array(resdat, mask=result.mask, copy=False)

~/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py in __array_ufunc__(self, function, method, *inputs, **kwargs)
    641             return result
    642 
--> 643         return self._result_as_quantity(result, unit, out)
    644 
    645     def _result_as_quantity(self, result, unit, out):

~/Library/Python/3.6/lib/python/site-packages/astropy/units/quantity.py in _result_as_quantity(self, result, unit, out)
    680         # the output is of the correct Quantity subclass, as it was passed
    681         # through check_output.
--> 682         out._set_unit(unit)
    683         return out
    684 

AttributeError: 'numpy.ndarray' object has no attribute '_set_unit'
<matplotlib.figure.Figure at 0x1088d3908>

In [ ]: