Source code for gammapy.data.obs_summary

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
import numpy as np
import astropy.units as u
from astropy.units import Quantity
from astropy.coordinates import SkyCoord
from .obs_stats import ObservationStats

__all__ = [
    'ObservationTableSummary',
    'ObservationSummary',
]


[docs]class ObservationTableSummary(object): """Observation table summary. Class allowing to summarize informations contained in Observation index table (`~gammapy.data.ObservationTable`) Parameters ---------- obs_table : `~gammapy.data.ObservationTable` Observation index table target_pos : `~astropy.coordinates.SkyCoord` Target position """ def __init__(self, obs_table, target_pos=None): self.obs_table = obs_table self.target_pos = target_pos @property def offset(self): """Observation pointing ot target offset (`~astropy.coordinates.Angle`). """ pnt_pos = SkyCoord(self.obs_table['RA_PNT'], self.obs_table['DEC_PNT'], unit='deg') offset = pnt_pos.separation(self.target_pos) return offset def __str__(self): """Summary report (`str`). """ ss = '*** Observation summary ***\n' ss += 'Target position: {}\n'.format(self.target_pos) ss += 'Number of observations: {}\n'.format(len(self.obs_table)) livetime = Quantity(sum(self.obs_table['LIVETIME']), 'second') ss += 'Livetime: {:.2f}\n'.format(livetime.to('hour')) zenith = self.obs_table['ZEN_PNT'] ss += 'Zenith angle: (mean={:.2f}, std={:.2f})\n'.format( zenith.mean(), zenith.std()) offset = self.offset ss += 'Offset: (mean={:.2f}, std={:.2f})\n'.format( offset.mean(), offset.std()) return ss
[docs] def plot_zenith_distribution(self, ax=None, bins=None): """Construct the zenith distribution of the observations. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. If None, uses the current `~matplotlib.axes.Axes`. bins : integer or array_like, optional Binning specification, passed to `matplotlib.pyplot.hist`. By default, 30 bins from 0 deg to max zenith + 5 deg is used. Returns -------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax zenith = self.obs_table['ZEN_PNT'] if bins is None: bins = np.linspace(0, zenith.max() + 5, 30) ax.hist(zenith, bins=bins) ax.set_title('Zenith distribution') ax.set_xlabel('Zenith (Deg)') ax.set_ylabel('#Entries') return ax
[docs] def plot_offset_distribution(self, ax=None, bins=None): """Construct the offset distribution of the observations. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. If None, uses the current `~matplotlib.axes.Axes`. bins : integer or array_like, optional Binning specification, passed to `matplotlib.pyplot.hist`. By default, 10 bins from 0 deg to max offset + 0.5 deg is used. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax offset = self.offset if bins is None: bins = np.linspace(0, offset.degree.max() + 0.5, 10) ax.hist(offset.degree, bins=bins) ax.set_title('Offset distribution') ax.set_xlabel('Offset (Deg)') ax.set_ylabel('#Entries') return ax
[docs]class ObservationSummary(object): """Summary of observations. For a list of observation stats, this class can make a table and do summary printout and plots. * TODO: Data should be stored in `~astropy.table.Table`! * TODO: there should be a per-run version of the stats in addition to the cumulative version that's there now. Parameters ---------- obs_stats : list Python list of `~gammapy.data.ObservationStats` """ def __init__(self, obs_stats): self.obs_stats = obs_stats self.obs_id = np.zeros(len(self.obs_stats)) self.livetime = np.zeros(len(self.obs_stats)) * u.s self.n_on = np.zeros(len(self.obs_stats)) self.n_off = np.zeros(len(self.obs_stats)) self.alpha = np.zeros(len(self.obs_stats)) self.background = np.zeros(len(self.obs_stats)) self.excess = np.zeros(len(self.obs_stats)) self.sigma = np.zeros(len(self.obs_stats)) self.gamma_rate = np.zeros(len(self.obs_stats)) / u.min self.bg_rate = np.zeros(len(self.obs_stats)) / u.min self._init_values() def _init_values(self): """Initialise vector attributes for plotting methods.""" stats_list = [] for index, obs in enumerate(self.obs_stats): # per observation stat self.obs_id[index] = obs.obs_id # keep track of the observation self.gamma_rate[index] = obs.gamma_rate self.bg_rate[index] = obs.bg_rate # cumulative information stats_list.append(obs) stack = ObservationStats.stack(stats_list) self.livetime[index] = stack.livetime self.n_on[index] = stack.n_on self.n_off[index] = stack.n_off self.alpha[index] = stack.alpha self.background[index] = stack.background self.excess[index] = stack.excess self.sigma[index] = stack.sigma
[docs] def obs_wise_summary(self): """Observation wise summary report (`str`).""" ss = '*** Observation Wise summary ***\n' for obs in self.obs_stats: ss += '{}\n'.format(obs) return ss
def __str__(self): """Observation summary report (`str`).""" stack = ObservationStats.stack(self.obs_stats) ss = '*** Observation summary ***\n' ss += '{}\n'.format(stack) return ss
[docs] def plot_significance_vs_livetime(self, ax=None, **kwargs): """Plot significance as a function of livetime. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax ax.plot(self.livetime.to(u.h), self.sigma, "o", **kwargs) ax.set_xlabel('Livetime ({0})'.format(u.h)) ax.set_ylabel('Significance ($\sigma$)') ax.axis([0., np.amax(self.livetime.to(u.h).value) * 1.2, 0., np.amax(self.sigma) * 1.2]) ax.set_title('Significance evolution') return ax
[docs] def plot_excess_vs_livetime(self, ax=None, **kwargs): """Plot excess as a function of livetime. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax ax.plot(self.livetime.to(u.h), self.excess, "o", **kwargs) ax.set_xlabel('Livetime ({0})'.format(u.h)) ax.set_ylabel('Excess') ax.axis([0., np.amax(self.livetime.to(u.h).value) * 1.2, 0., np.amax(self.excess) * 1.2]) ax.set_title('Excess evolution') return ax
[docs] def plot_background_vs_livetime(self, ax=None, **kwargs): """Plot background as a function of livetime. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax ax.plot(self.livetime.to(u.h), self.background, "o", **kwargs) ax.set_xlabel('Livetime ({0})'.format(u.h)) ax.set_ylabel('Background') ax.axis([0., np.amax(self.livetime.to(u.h).value) * 1.2, 0., np.amax(self.background) * 1.2]) ax.set_title('Background evolution') return ax
[docs] def plot_gamma_rate(self, ax=None, **kwargs): """Plot gamma rate for each observation. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax labels = list() values = list() for index in range(len(self.gamma_rate)): labels.append(str(int(self.obs_id[index]))) values.append(index + 0.5) ax.plot(values, self.gamma_rate, "o", **kwargs) ax.set_xlabel('Observation Ids') ax.set_xticks(values) ax.set_xticklabels(labels, rotation=-22.5) ax.set_ylabel('$\gamma$ rate ({})'.format(self.gamma_rate.unit)) ax.axis([0, len(self.gamma_rate), 0., np.amax(self.gamma_rate.value) * 1.2]) ax.set_title('$\gamma$ rates') return ax
[docs] def plot_background_rate(self, ax=None, **kwargs): """Plot background rate for each observation. Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. The `~matplotlib.axes.Axes` object to be drawn on. Returns ------- ax : `~matplotlib.axes.Axes` Axis """ import matplotlib.pyplot as plt ax = plt.gca() if ax is None else ax labels = list() values = list() for index in range(len(self.bg_rate)): labels.append(str(int(self.obs_id[index]))) values.append(index + 0.5) ax.plot(values, self.bg_rate, "o", **kwargs) ax.set_xlabel('Observation Ids') ax.set_xticks(values) ax.set_xticklabels(labels, rotation=-22.5) ax.set_ylabel('Background rate ({})'.format(self.bg_rate.unit)) ax.axis([0, len(self.bg_rate), 0., np.amax(self.bg_rate.value) * 1.2]) ax.set_title('Background rates') return ax