diff --git a/gammapy/data/__init__.py b/gammapy/data/__init__.py index ce34119f9b..ce3664fbc5 100644 --- a/gammapy/data/__init__.py +++ b/gammapy/data/__init__.py @@ -18,3 +18,4 @@ class InvalidDataError(Exception): from .utils import * from .observation_summary import * from .pointing import * +from .observation_stats import * diff --git a/gammapy/data/observation_stats.py b/gammapy/data/observation_stats.py new file mode 100644 index 0000000000..cecddac0fc --- /dev/null +++ b/gammapy/data/observation_stats.py @@ -0,0 +1,217 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import absolute_import, division, print_function, unicode_literals +import astropy.units as u +from ..stats import Stats +from gammapy.stats import significance_on_off + +__all__ = [ + 'ObservationStats', + 'ObservationStatsList' +] + + +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 obervation + 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 + self.gamma_rate = gamma_rate + self.bg_rate = bg_rate + + @classmethod + def from_target(cls, obs, target, bg_estimate): + """ + Auxiliary constructor from an osbervation, a target + a background estimate + + Parameters + ---------- + obs_table : `~gammapy.data.ObservationTable` + Observation index table + target_pos : `~astropy.coordinates.SkyCoord` + Target position + bg_method : str + Background estimation method + """ + n_on = cls._get_on_events(obs, target) + n_off = len(bg_estimate.off_events) + 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): + """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 + + @staticmethod + def _get_on_events(obs, target): + """Number of ON events in the region of interest (`int`) + """ + idx = target.on_region.contains(obs.events.radec) + on_events = obs.events[idx] + return len(on_events) + + @classmethod + def stack(cls, stats_list): + """Stack statistics from a list of + `~gammapy.data.ObservationStats` and returns a new instance + of `~gammapy.data.ObservationStats` + + Parameters + ---------- + stats_list : list + List of observation statistics `~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 + + a_on /= n_off + a_off /= n_off + alpha /= n_off + + # 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 + + gamma_rate /= livetime.to(u.min) + bg_rate /= livetime.to(u.min) + + total_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 total_stats + + def __add__(self, other): + """Add statistics from two observations + and returns new instance of `~gammapy.data.ObservationStats` + """ + return ObservationStats.stack(self, other) + + 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: {}\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: {}\n'.format(self.alpha) + ss += 'Bkg events in On region: {}\n'.format(self.background) + ss += 'Excess: {}\n'.format(self.excess) + ss += 'Gamma rate: {}\n'.format(self.gamma_rate) + ss += 'Bkg rate: {}\n'.format(self.bg_rate) + ss += 'Sigma: {}\n'.format(self.sigma) + + return ss + + +class ObservationStatsList(list): + """List of `~gammapy.data.ObservationList` + """ diff --git a/gammapy/data/observation_summary.py b/gammapy/data/observation_summary.py index 1cc4afac91..9b08e004d0 100644 --- a/gammapy/data/observation_summary.py +++ b/gammapy/data/observation_summary.py @@ -1,11 +1,14 @@ # 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 .observation_stats import ObservationStats, ObservationStatsList __all__ = [ 'ObservationTableSummary', + 'ObservationSummary', ] @@ -41,7 +44,7 @@ def offset(self): def plot_zenith_distribution(self, ax=None, bins=None): """Construct the zenith distribution of the observations. - + Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. @@ -73,7 +76,7 @@ def plot_zenith_distribution(self, ax=None, bins=None): def plot_offset_distribution(self, ax=None, bins=None): """Construct the offset distribution of the observations. - + Parameters ---------- ax : `~matplotlib.axes.Axes` or None, optional. @@ -125,3 +128,207 @@ def show_in_browser(self): """Make HTML file and images in tmp dir, open in browser. """ raise NotImplementedError + + +class ObservationSummary(object): + """Summary of observations. + + Class allowing to summarise informations contained in + a list of observations (`~gammapy.data.ObservationStatsList`) + + Parameters + ---------- + obs_stats : `~gammapy.data.ObservationStatsList` + List of observation statistics + """ + + 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 attributs for plotting methods + """ + cumul_obs = ObservationStatsList() + 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 + cumul_obs.append(obs) + stack = ObservationStats.stack(cumul_obs) + 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 + + 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 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 : `~matplolib.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 + + 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 : `~matplolib.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 + + 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 : `~matplolib.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 + + 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 : `~matplolib.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 + + 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 : `~matplolib.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 + + def __str__(self): + """Observation summary report (`str`) + """ + stack = ObservationStats.stack(self.obs_stats) + ss = '*** Observation summary ***\n' + ss += '{}\n'.format(stack) + return ss diff --git a/gammapy/data/tests/test_observation_stats.py b/gammapy/data/tests/test_observation_stats.py new file mode 100644 index 0000000000..9b37b3d2c5 --- /dev/null +++ b/gammapy/data/tests/test_observation_stats.py @@ -0,0 +1,71 @@ +# Licensed under a 3-clause BSD style license - see LICENSE.rst +from __future__ import absolute_import, division, print_function, unicode_literals +from numpy.testing import assert_allclose +from astropy.tests.helper import pytest +from astropy.coordinates import SkyCoord +import astropy.units as u +from ...data import DataStore, ObservationList, ObservationStats, Target +from ...utils.testing import requires_data +from ...extern.regions.shapes import CircleSkyRegion +from ...background import reflected_regions_background_estimate as refl +from ...image import ExclusionMask + + +@requires_data('gammapy-extra') +def get_obs_list(): + data_store = DataStore.from_dir( + '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/') + run_list = [23523, 23526] + obs_list = ObservationList([data_store.obs(_) for _ in run_list]) + return obs_list + + +@requires_data('gammapy-extra') +def get_obs(id): + obs_list = get_obs_list() + for run in obs_list: + if run.obs_id == id: + return run + +@pytest.fixture +def target(): + pos = SkyCoord(83.63 * u.deg, 22.01 * u.deg, frame='icrs') + on_size = 0.3 * u.deg + on_region = CircleSkyRegion(pos, on_size) + + target = Target(position=pos, + on_region=on_region, + name='Crab Nebula', + tag='crab') + return target + + +@requires_data('gammapy-extra') +@pytest.fixture +def get_mask(): + mask = ExclusionMask.read( + '$GAMMAPY_EXTRA/datasets/exclusion_masks/tevcat_exclusion.fits') + return mask + + +@requires_data('gammapy-extra') +def test_str(target): + run = get_obs(23523) + bg = refl(target.on_region, + run.pointing_radec, get_mask(), run.events) + stats = ObservationStats.from_target(run, target, bg) + text = str(stats) + assert 'Observation summary report' in text + + +@requires_data('gammapy-extra') +def test_stack(target): + obs_list = get_obs_list() + obs_stats = list() + for run in obs_list: + bg = refl(target.on_region, + run.pointing_radec, get_mask(), run.events) + obs_stats.append(ObservationStats.from_target(run, target, bg)) + sum_obs_stats = ObservationStats.stack(obs_stats) + assert_allclose(sum_obs_stats.alpha, 0.284, rtol=1.e-2) + assert_allclose(sum_obs_stats.sigma, 23.48, rtol=1.e-3) diff --git a/gammapy/data/tests/test_observation_summary.py b/gammapy/data/tests/test_observation_summary.py index e4da3d40de..aca606c600 100644 --- a/gammapy/data/tests/test_observation_summary.py +++ b/gammapy/data/tests/test_observation_summary.py @@ -3,38 +3,112 @@ from numpy.testing import assert_allclose from astropy.tests.helper import pytest from astropy.coordinates import SkyCoord -from ...data import DataStore, ObservationTableSummary +import astropy.units as u +from ...data import DataStore, ObservationTableSummary, ObservationSummary +from ...data import ObservationStats, ObservationStatsList, ObservationList +from ...data import Target from ...utils.testing import requires_data, requires_dependency +from ...extern.regions.shapes import CircleSkyRegion +from ...background import reflected_regions_background_estimate as refl +from ...image import ExclusionMask @requires_data('gammapy-extra') @pytest.fixture -def summary(): - data_store = DataStore.from_dir('$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/') +def table_summary(): + data_store = DataStore.from_dir( + '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/') target_pos = SkyCoord(83.633083, 22.0145, unit='deg') return ObservationTableSummary(data_store.obs_table, target_pos) @requires_data('gammapy-extra') -def test_str(summary): - text = str(summary) - assert ('Observation summary' in text) +def test_str(table_summary): + text = str(table_summary) + assert 'Observation summary' in text @requires_data('gammapy-extra') -def test_offset(summary): - offset = summary.offset +def test_offset(table_summary): + offset = table_summary.offset assert_allclose(offset.degree.mean(), 1., rtol=1.e-2) assert_allclose(offset.degree.std(), 0.5, rtol=1.e-2) @requires_data('gammapy-extra') @requires_dependency('matplotlib') -def test_plot_zenith(summary): - summary.plot_zenith_distribution() +def test_plot_zenith(table_summary): + table_summary.plot_zenith_distribution() @requires_data('gammapy-extra') @requires_dependency('matplotlib') -def test_plot_offset(summary): - summary.plot_offset_distribution() +def test_plot_offset(table_summary): + table_summary.plot_offset_distribution() + + +@requires_data('gammapy-extra') +@pytest.fixture +def obs_summary(): + datastore = DataStore.from_dir( + '$GAMMAPY_EXTRA/datasets/hess-crab4-hd-hap-prod2/') + run_list = [23523, 23526, 23559, 23592] + + pos = SkyCoord(83.63 * u.deg, 22.01 * u.deg, frame='icrs') + on_size = 0.3 * u.deg + on_region = CircleSkyRegion(pos, on_size) + + target = Target(position=pos, on_region=on_region, + name='Crab Nebula', tag='crab') + + mask = ExclusionMask.read( + '$GAMMAPY_EXTRA/datasets/exclusion_masks/tevcat_exclusion.fits') + + obs_list = ObservationList([datastore.obs(_) for _ in run_list]) + obs_stats = ObservationStatsList() + + for index, run in enumerate(obs_list): + + bkg = refl(on_region, run.pointing_radec, mask, run.events) + + obs_stats.append(ObservationStats.from_target(run, target, bkg)) + + summary = ObservationSummary(obs_stats) + + return summary + + +@requires_data('gammapy-extra') +@requires_dependency('matplotlib') +def test_plot_significance(obs_summary): + obs_summary.plot_significance_vs_livetime() + + +@requires_data('gammapy-extra') +@requires_dependency('matplotlib') +def test_plot_excess(obs_summary): + obs_summary.plot_excess_vs_livetime() + + +@requires_data('gammapy-extra') +@requires_dependency('matplotlib') +def test_plot_background(obs_summary): + obs_summary.plot_background_vs_livetime() + + +@requires_data('gammapy-extra') +@requires_dependency('matplotlib') +def test_plot_gamma_rate(obs_summary): + obs_summary.plot_gamma_rate() + + +@requires_data('gammapy-extra') +@requires_dependency('matplotlib') +def test_plot_background_rate(obs_summary): + obs_summary.plot_background_rate() + + +@requires_data('gammapy-extra') +def test_obs_str(obs_summary): + text = str(obs_summary) + assert 'Observation summary' in text diff --git a/gammapy/stats/data.py b/gammapy/stats/data.py index d0c173db4a..da75776de2 100644 --- a/gammapy/stats/data.py +++ b/gammapy/stats/data.py @@ -30,10 +30,10 @@ class Stats(object): # TODO: add gamma exposure def __init__(self, n_on, n_off, a_on, a_off): - self.n_on = float(n_on) - self.n_off = float(n_off) - self.a_on = float(a_on) - self.a_off = float(a_off) + self.n_on = n_on + self.n_off = n_off + self.a_on = a_on + self.a_off = a_off @property def alpha(self): @@ -63,7 +63,7 @@ def __str__(self): keys = ['n_on', 'n_off', 'a_on', 'a_off', 'alpha', 'background', 'excess'] values = [self.n_on, self.n_off, self.a_on, self.a_off, - self.alpha, self.background(), self.excess()] + self.alpha, self.background, self.excess] return '\n'.join(['%s = %s' % (k, v) for (k, v) in zip(keys, values)]) diff --git a/gammapy/stats/tests/test_data.py b/gammapy/stats/tests/test_data.py index 1978911b52..5fa8944e6c 100644 --- a/gammapy/stats/tests/test_data.py +++ b/gammapy/stats/tests/test_data.py @@ -12,6 +12,12 @@ def test_Stats(): assert_allclose(stats.background, a_on / a_off * n_off) assert_allclose(stats.excess, n_on - a_on / a_off * n_off) +def test_str(): + data = Stats(10, 10, 1, 10) + text = str(data) + assert ('alpha = 0.1' in text) + assert ('background = 1.0' in text) + assert ('excess = 9.0' in text) def test_make_stats(): pass