Source code for gammapy.detect.lima

# Licensed under a 3-clause BSD style license - see LICENSE.rst
import copy
import logging
import numpy as np
from ..stats import significance, significance_on_off

__all__ = ["compute_lima_image", "compute_lima_on_off_image"]

log = logging.getLogger(__name__)


[docs]def compute_lima_image(counts, background, kernel): """Compute Li & Ma significance and flux images for known background. Parameters ---------- counts : `~gammapy.maps.WcsNDMap` Counts image background : `~gammapy.maps.WcsNDMap` Background image kernel : `astropy.convolution.Kernel2D` Convolution kernel Returns ------- images : dict Dictionary containing result maps Keys are: significance, counts, background and excess See Also -------- gammapy.stats.significance """ # Kernel is modified later make a copy here kernel = copy.deepcopy(kernel) kernel.normalize("peak") # fft convolution adds numerical noise, to ensure integer results we call # np.rint counts_conv = np.rint(counts.convolve(kernel.array).data) background_conv = background.convolve(kernel.array).data excess_conv = counts_conv - background_conv significance_conv = significance(counts_conv, background_conv, method="lima") return { "significance": counts.copy(data=significance_conv), "counts": counts.copy(data=counts_conv), "background": counts.copy(data=background_conv), "excess": counts.copy(data=excess_conv), }
[docs]def compute_lima_on_off_image(n_on, n_off, a_on, a_off, kernel): """Compute Li & Ma significance and flux images for on-off observations. Parameters ---------- n_on : `~gammapy.maps.WcsNDMap` Counts image n_off : `~gammapy.maps.WcsNDMap` Off counts image a_on : `~gammapy.maps.WcsNDMap` Relative background efficiency in the on region a_off : `~gammapy.maps.WcsNDMap` Relative background efficiency in the off region kernel : `astropy.convolution.Kernel2D` Convolution kernel Returns ------- images : dict Dictionary containing result maps Keys are: significance, n_on, background, excess, alpha See Also -------- gammapy.stats.significance_on_off """ # Kernel is modified later make a copy here kernel = copy.deepcopy(kernel) kernel.normalize("peak") # fft convolution adds numerical noise, to ensure integer results we call # np.rint n_on_conv = np.rint(n_on.convolve(kernel.array).data) a_on_conv = a_on.convolve(kernel.array).data with np.errstate(invalid="ignore", divide="ignore"): alpha_conv = a_on_conv / a_off.data significance_conv = significance_on_off( n_on_conv, n_off.data, alpha_conv, method="lima" ) with np.errstate(invalid="ignore"): background_conv = alpha_conv * n_off.data excess_conv = n_on_conv - background_conv return { "significance": n_on.copy(data=significance_conv), "n_on": n_on.copy(data=n_on_conv), "background": n_on.copy(data=background_conv), "excess": n_on.copy(data=excess_conv), "alpha": n_on.copy(data=alpha_conv), }