Source code for gammapy.data.obs_stats

# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import absolute_import, division, print_function, unicode_literals
from collections import OrderedDict
import numpy as np
import astropy.units as u
from ..stats import Stats, significance_on_off

__all__ = [
    'ObservationStats',
]


[docs]class ObservationStats(Stats): """Observation statistics. Class allowing to summarize observation (`~gammapy.data.DataStoreObservation`) statistics Parameters ---------- n_on : int Number of on events n_off : int Number of off events a_on : float Relative background exposure of the on region a_off : float Relative background exposure of the off region obs_id : int ID of the observation livetime : float Livetime of the observation alpha : float Normalisation between the on and the off exposure bg_rate : float Background rate (/min) gamma_rate : float Gamma rate (/min) """ def __init__(self, n_on=None, n_off=None, a_on=None, a_off=None, obs_id=None, livetime=None, alpha=None, gamma_rate=None, bg_rate=None): super(ObservationStats, self).__init__( n_on=n_on, n_off=n_off, a_on=a_on, a_off=a_off, ) self.obs_id = obs_id self.livetime = livetime self.alpha_obs = alpha or a_on / a_off self.gamma_rate = gamma_rate or n_on / livetime self.bg_rate = bg_rate or self.alpha_obs * n_off / livetime @classmethod
[docs] def from_obs(cls, obs, bg_estimate): """Create from `~gammapy.data.DataStoreObservation`. Parameters ---------- obs : `~gammapy.data.DataStoreObservation` IACT data store observation bg_estimate : `~gammapy.background.BackgroundEstimate` Background estimate """ n_on = len(bg_estimate.on_events.table) n_off = len(bg_estimate.off_events.table) a_on = bg_estimate.a_on a_off = bg_estimate.a_off obs_id = obs.obs_id livetime = obs.observation_live_time_duration alpha = a_on / a_off gamma_rate = n_on / livetime.to(u.min) bg_rate = (alpha * n_off) / livetime.to(u.min) stats = cls(n_on=n_on, n_off=n_off, a_on=a_on, a_off=a_off, obs_id=obs_id, livetime=livetime, alpha=alpha, gamma_rate=gamma_rate, bg_rate=bg_rate) return stats
@property def alpha(self): """Alpha (on / off exposure ratio) Override member function from `~gammapy.stats.Stats` to take into account weighted alpha by number of Off events """ return self.alpha_obs @property def sigma(self): """Li-Ma significance for observation statistics (`float`). """ sigma = significance_on_off( self.n_on, self.n_off, self.alpha, method='lima') return sigma @classmethod
[docs] def stack(cls, stats_list): """Stack (concatenate) list of `~gammapy.data.ObservationStats`. Parameters ---------- stats_list : list List of `~gammapy.data.ObservationStats` Returns ------- total_stats : `~gammapy.data.ObservationStats` Statistics for stacked observation """ n_on = 0 n_off = 0 a_on = 0 a_on_backup = 0 a_off = 0 a_off_backup = 0 obs_id = list() livetime = 0 alpha = 0 alpha_backup = 0 gamma_rate = 0 bg_rate = 0 for stats in stats_list: livetime += stats.livetime n_on += stats.n_on n_off += stats.n_off a_on += stats.a_on * stats.n_off a_on_backup += stats.a_on * stats.livetime.value a_off += stats.a_off * stats.n_off a_off_backup += stats.a_off * stats.livetime.value alpha += stats.alpha * stats.n_off alpha_backup += stats.alpha * stats.livetime.value obs_id.append(stats.obs_id) gamma_rate += stats.n_on - stats.alpha * stats.n_off bg_rate += stats.n_off * stats.alpha # if no off events the weighting of alpha is done # with the livetime if n_off == 0: alpha = alpha_backup / livetime.value a_on = a_on_backup / livetime.value a_off = a_off_backup / livetime.value else: a_on /= n_off a_off /= n_off alpha /= n_off gamma_rate /= livetime.to(u.min) bg_rate /= livetime.to(u.min) return cls( n_on=n_on, n_off=n_off, a_on=a_on, a_off=a_off, obs_id=obs_id, livetime=livetime, alpha=alpha, gamma_rate=gamma_rate, bg_rate=bg_rate )
[docs] def to_dict(self): """Data as an `~collections.OrderedDict`. This is useful for serialisation or putting the info in a table. """ data = OrderedDict() data['obs_id'] = self.obs_id data['livetime'] = self.livetime data['n_on'] = self.n_on data['n_off'] = self.n_off data['a_on'] = self.a_on data['a_off'] = self.a_off data['alpha'] = self.alpha data['background'] = self.background data['excess'] = self.excess data['sigma'] = self.sigma data['gamma_rate'] = self.gamma_rate data['bg_rate'] = self.bg_rate return data
def __str__(self): """Observation statistics report (`str`).""" ss = '*** Observation summary report ***\n' if type(self.obs_id) is list: obs_str = '[{0}-{1}]'.format(self.obs_id[0], self.obs_id[-1]) else: obs_str = '{}'.format(self.obs_id) ss += 'Observation Id: {}\n'.format(obs_str) ss += 'Livetime: {:.3f}\n'.format(self.livetime.to(u.h)) ss += 'On events: {}\n'.format(self.n_on) ss += 'Off events: {}\n'.format(self.n_off) ss += 'Alpha: {:.3f}\n'.format(self.alpha) ss += 'Bkg events in On region: {:.2f}\n'.format(self.background) ss += 'Excess: {:.2f}\n'.format(self.excess) ss += 'Excess / Background: {:.2f}\n'.format(np.divide(self.excess, self.background)) ss += 'Gamma rate: {:.2f}\n'.format(self.gamma_rate) ss += 'Bkg rate: {:.2f}\n'.format(self.bg_rate) ss += 'Sigma: {:.2f}\n'.format(self.sigma) return ss