.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "tutorials/analysis-2d/ring_background.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. or to run this example in your browser via Binder .. rst-class:: sphx-glr-example-title .. _sphx_glr_tutorials_analysis-2d_ring_background.py: Ring background map =================== Create an excess (gamma-ray events) and a significance map extracting a ring background. Context ------- One of the challenges of IACT analysis is accounting for the large residual hadronic emission. An excess map, assumed to be a map of only gamma-ray events, requires a good estimate of the background. However, in the absence of a solid template bkg model it is not possible to obtain reliable background model a priori. It was often found necessary in classical cherenkov astronomy to perform a local renormalization of the existing templates, usually with a ring kernel. This assumes that most of the events are background and requires to have an exclusion mask to remove regions with bright signal from the estimation. To read more about this method, see `here. `__ Objective --------- Create an excess (gamma-ray events) map of MSH 15-52 as well as a significance map to determine how solid the signal is. Proposed approach ----------------- The analysis workflow is roughly: - Compute the sky maps keeping each observation separately using the `~gammapy.analysis.Analysis` class - Estimate the background using the `~gammapy.makers.RingBackgroundMaker` - Compute the correlated excess and significance maps using the `~gammapy.estimators.ExcessMapEstimator` The normalised background thus obtained can be used for general modelling and fitting. .. GENERATED FROM PYTHON SOURCE LINES 43-48 Setup ----- As usual, we’ll start with some general imports… .. GENERATED FROM PYTHON SOURCE LINES 48-66 .. code-block:: Python import logging # %matplotlib inline import astropy.units as u from astropy.coordinates import SkyCoord from regions import CircleSkyRegion import matplotlib.pyplot as plt from gammapy.analysis import Analysis, AnalysisConfig from gammapy.datasets import MapDatasetOnOff from gammapy.estimators import ExcessMapEstimator from gammapy.makers import RingBackgroundMaker from gammapy.visualization import plot_distribution from gammapy.utils.check import check_tutorials_setup log = logging.getLogger(__name__) .. GENERATED FROM PYTHON SOURCE LINES 67-69 Check setup ----------- .. GENERATED FROM PYTHON SOURCE LINES 69-74 .. code-block:: Python check_tutorials_setup() .. rst-class:: sphx-glr-script-out .. code-block:: none System: python_executable : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/bin/python python_version : 3.9.20 machine : x86_64 system : Linux Gammapy package: version : 1.3 path : /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy Other packages: numpy : 1.26.4 scipy : 1.13.1 astropy : 5.2.2 regions : 0.8 click : 8.1.7 yaml : 6.0.2 IPython : 8.18.1 jupyterlab : not installed matplotlib : 3.9.2 pandas : not installed healpy : 1.17.3 iminuit : 2.30.1 sherpa : 4.16.1 naima : 0.10.0 emcee : 3.1.6 corner : 2.2.3 ray : 2.39.0 Gammapy environment variables: GAMMAPY_DATA : /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3 .. GENERATED FROM PYTHON SOURCE LINES 75-84 Creating the config file ------------------------ Now, we create a config file for out analysis. You may load this from disc if you have a pre-defined config file. In this example, we will use a few H.E.S.S. runs on the pulsar wind nebula, MSH 1552 .. GENERATED FROM PYTHON SOURCE LINES 84-121 .. code-block:: Python # source_pos = SkyCoord.from_name("MSH 15-52") source_pos = SkyCoord(228.32, -59.08, unit="deg") config = AnalysisConfig() # Select observations - 2.5 degrees from the source position config.observations.datastore = "$GAMMAPY_DATA/hess-dl3-dr1/" config.observations.obs_cone = { "frame": "icrs", "lon": source_pos.ra, "lat": source_pos.dec, "radius": 2.5 * u.deg, } config.datasets.type = "3d" config.datasets.geom.wcs.skydir = { "lon": source_pos.ra, "lat": source_pos.dec, "frame": "icrs", } # The WCS geometry - centered on MSH 15-52 config.datasets.geom.wcs.width = {"width": "3 deg", "height": "3 deg"} config.datasets.geom.wcs.binsize = "0.02 deg" # Cutout size (for the run-wise event selection) config.datasets.geom.selection.offset_max = 2.5 * u.deg # We now fix the energy axis for the counts map - (the reconstructed energy binning) config.datasets.geom.axes.energy.min = "0.5 TeV" config.datasets.geom.axes.energy.max = "10 TeV" config.datasets.geom.axes.energy.nbins = 10 # We need to extract the ring for each observation separately, hence, no stacking at this stage config.datasets.stack = False print(config) .. rst-class:: sphx-glr-script-out .. code-block:: none AnalysisConfig general: log: level: info filename: null filemode: null format: null datefmt: null outdir: . n_jobs: 1 datasets_file: null models_file: null observations: datastore: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3/hess-dl3-dr1 obs_ids: [] obs_file: null obs_cone: frame: icrs lon: 228.32 deg lat: -59.08 deg radius: 2.5 deg obs_time: start: null stop: null required_irf: - aeff - edisp - psf - bkg datasets: type: 3d stack: false geom: wcs: skydir: frame: icrs lon: 228.32 deg lat: -59.08 deg binsize: 0.02 deg width: width: 3.0 deg height: 3.0 deg binsize_irf: 0.2 deg selection: offset_max: 2.5 deg axes: energy: min: 0.5 TeV max: 10.0 TeV nbins: 10 energy_true: min: 0.5 TeV max: 20.0 TeV nbins: 16 map_selection: - counts - exposure - background - psf - edisp background: method: null exclusion: null parameters: {} safe_mask: methods: - aeff-default parameters: {} on_region: frame: null lon: null lat: null radius: null containment_correction: true fit: fit_range: min: null max: null flux_points: energy: min: null max: null nbins: null source: source parameters: selection_optional: all excess_map: correlation_radius: 0.1 deg parameters: {} energy_edges: min: null max: null nbins: null light_curve: time_intervals: start: null stop: null energy_edges: min: null max: null nbins: null source: source parameters: selection_optional: all metadata: creator: Gammapy 1.3 date: '2024-11-26T10:11:20.034778' origin: null .. GENERATED FROM PYTHON SOURCE LINES 122-129 Getting the reduced dataset --------------------------- We now use the config file to do the initial data reduction which will then be used for a ring extraction Create the config: .. GENERATED FROM PYTHON SOURCE LINES 129-142 .. code-block:: Python analysis = Analysis(config) # for this specific case,w e do not need fine bins in true energy analysis.config.datasets.geom.axes.energy_true = ( analysis.config.datasets.geom.axes.energy ) # First get the required observations analysis.get_observations() print(analysis.config) .. rst-class:: sphx-glr-script-out .. code-block:: none AnalysisConfig general: log: level: INFO filename: null filemode: null format: null datefmt: null outdir: . n_jobs: 1 datasets_file: null models_file: null observations: datastore: /home/runner/work/gammapy-docs/gammapy-docs/gammapy-datasets/1.3/hess-dl3-dr1 obs_ids: [] obs_file: null obs_cone: frame: icrs lon: 228.32 deg lat: -59.08 deg radius: 2.5 deg obs_time: start: null stop: null required_irf: - aeff - edisp - psf - bkg datasets: type: 3d stack: false geom: wcs: skydir: frame: icrs lon: 228.32 deg lat: -59.08 deg binsize: 0.02 deg width: width: 3.0 deg height: 3.0 deg binsize_irf: 0.2 deg selection: offset_max: 2.5 deg axes: energy: min: 0.5 TeV max: 10.0 TeV nbins: 10 energy_true: min: 0.5 TeV max: 10.0 TeV nbins: 10 map_selection: - counts - exposure - background - psf - edisp background: method: null exclusion: null parameters: {} safe_mask: methods: - aeff-default parameters: {} on_region: frame: null lon: null lat: null radius: null containment_correction: true fit: fit_range: min: null max: null flux_points: energy: min: null max: null nbins: null source: source parameters: selection_optional: all excess_map: correlation_radius: 0.1 deg parameters: {} energy_edges: min: null max: null nbins: null light_curve: time_intervals: start: null stop: null energy_edges: min: null max: null nbins: null source: source parameters: selection_optional: all metadata: creator: Gammapy 1.3 date: '2024-11-26T10:11:20.094200' origin: null .. GENERATED FROM PYTHON SOURCE LINES 143-144 Data extraction: .. GENERATED FROM PYTHON SOURCE LINES 144-147 .. code-block:: Python analysis.get_datasets() .. GENERATED FROM PYTHON SOURCE LINES 148-155 Extracting the ring background ------------------------------ Since the ring background is extracted from real off events, we need to use the Wstat statistics in this case. For this, we will use the `~gammapy.datasets.MapDatasetOnOff` and the `~gammapy.makers.RingBackgroundMaker` classes. .. GENERATED FROM PYTHON SOURCE LINES 158-165 Create exclusion mask ~~~~~~~~~~~~~~~~~~~~~ First, we need to create an exclusion mask on the known sources. In this case, we need to mask only MSH 15-52 but this depends on the sources present in our field of view. .. GENERATED FROM PYTHON SOURCE LINES 165-178 .. code-block:: Python # get the geom that we use geom = analysis.datasets[0].counts.geom energy_axis = analysis.datasets[0].counts.geom.axes["energy"] geom_image = geom.to_image().to_cube([energy_axis.squash()]) # Make the exclusion mask regions = CircleSkyRegion(center=source_pos, radius=0.4 * u.deg) exclusion_mask = ~geom_image.region_mask([regions]) exclusion_mask.sum_over_axes().plot() plt.show() .. image-sg:: /tutorials/analysis-2d/images/sphx_glr_ring_background_001.png :alt: ring background :srcset: /tutorials/analysis-2d/images/sphx_glr_ring_background_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 179-182 For the present analysis, we use a ring with an inner radius of 0.5 deg and width of 0.3 deg. .. GENERATED FROM PYTHON SOURCE LINES 182-188 .. code-block:: Python ring_maker = RingBackgroundMaker( r_in="0.5 deg", width="0.3 deg", exclusion_mask=exclusion_mask ) .. GENERATED FROM PYTHON SOURCE LINES 189-195 Create a stacked dataset ~~~~~~~~~~~~~~~~~~~~~~~~ Now, we extract the background for each dataset and then stack the maps together to create a single stacked map for further analysis .. GENERATED FROM PYTHON SOURCE LINES 195-201 .. code-block:: Python energy_axis_true = analysis.datasets[0].exposure.geom.axes["energy_true"] stacked_on_off = MapDatasetOnOff.create( geom=geom_image, energy_axis_true=energy_axis_true, name="stacked" ) .. GENERATED FROM PYTHON SOURCE LINES 202-208 .. code-block:: Python for dataset in analysis.datasets: # Ring extracting makes sense only for 2D analysis dataset_on_off = ring_maker.run(dataset.to_image()) stacked_on_off.stack(dataset_on_off) .. GENERATED FROM PYTHON SOURCE LINES 209-214 This `stacked_on_off` has `on` and `off` counts and acceptance maps which we will use in all further analysis. The `acceptance` and `acceptance_off` maps are the system acceptance of gamma-ray like events in the `on` and `off` regions respectively. .. GENERATED FROM PYTHON SOURCE LINES 214-218 .. code-block:: Python print(stacked_on_off) .. rst-class:: sphx-glr-script-out .. code-block:: none MapDatasetOnOff --------------- Name : stacked Total counts : 41803 Total background counts : 40671.90 Total excess counts : 1131.10 Predicted counts : 40672.37 Predicted background counts : 40672.37 Predicted excess counts : nan Exposure min : 1.15e+09 m2 s Exposure max : 1.33e+10 m2 s Number of total bins : 22500 Number of fit bins : 22500 Fit statistic type : wstat Fit statistic value (-2 log(L)) : 26475.50 Number of models : 0 Number of parameters : 0 Number of free parameters : 0 Total counts_off : 87621184 Acceptance : 44755 Acceptance off : 96374320 .. GENERATED FROM PYTHON SOURCE LINES 219-230 Compute correlated significance and correlated excess maps ---------------------------------------------------------- We need to convolve our maps with an appropriate smoothing kernel. The significance is computed according to the Li & Ma expression for ON and OFF Poisson measurements, see `here `__. Also, since the off counts are obtained with a ring background estimation, and are thus already correlated, we specify `correlate_off=False` to avoid over correlation. .. GENERATED FROM PYTHON SOURCE LINES 230-248 .. code-block:: Python # Using a convolution radius of 0.04 degrees estimator = ExcessMapEstimator(0.04 * u.deg, selection_optional=[], correlate_off=False) lima_maps = estimator.run(stacked_on_off) significance_map = lima_maps["sqrt_ts"] excess_map = lima_maps["npred_excess"] # We can plot the excess and significance maps fig, (ax1, ax2) = plt.subplots( figsize=(11, 4), subplot_kw={"projection": lima_maps.geom.wcs}, ncols=2 ) ax1.set_title("Significance map") significance_map.plot(ax=ax1, add_cbar=True) ax2.set_title("Excess map") excess_map.plot(ax=ax2, add_cbar=True) plt.show() .. image-sg:: /tutorials/analysis-2d/images/sphx_glr_ring_background_002.png :alt: Significance map, Excess map :srcset: /tutorials/analysis-2d/images/sphx_glr_ring_background_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 249-256 It is often important to look at the significance distribution outside the exclusion region to check that the background estimation is not contaminated by gamma-ray events. This can be the case when exclusion regions are not large enough. Typically, we expect the off distribution to be a standard normal distribution. To compute the significance distribution outside the exclusion region, we can recompute the maps after adding a `mask_fit` to our dataset. .. GENERATED FROM PYTHON SOURCE LINES 256-291 .. code-block:: Python # Mask the regions with gamma-ray emission stacked_on_off.mask_fit = exclusion_mask lima_maps2 = estimator.run(stacked_on_off) significance_map_off = lima_maps2["sqrt_ts"] kwargs_axes = {"xlabel": "Significance", "yscale": "log", "ylim": (1e-3, 1)} ax, _ = plot_distribution( significance_map, kwargs_hist={ "density": True, "alpha": 0.5, "color": "red", "label": "all bins", "bins": 51, }, kwargs_axes=kwargs_axes, ) ax, res = plot_distribution( significance_map_off, ax=ax, func="norm", kwargs_hist={ "density": True, "alpha": 0.5, "color": "blue", "label": "off bins", "bins": 51, }, kwargs_axes=kwargs_axes, ) plt.show() # sphinx_gallery_thumbnail_number = 2 .. image-sg:: /tutorials/analysis-2d/images/sphx_glr_ring_background_003.png :alt: ring background :srcset: /tutorials/analysis-2d/images/sphx_glr_ring_background_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/stats/counts_statistic.py:398: RuntimeWarning: invalid value encountered in multiply return self.alpha * self.n_off /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/stats/fit_statistics.py:207: RuntimeWarning: invalid value encountered in multiply C = alpha * (n_on + n_off) - (1 + alpha) * mu_sig /home/runner/work/gammapy-docs/gammapy-docs/gammapy/.tox/build_docs/lib/python3.9/site-packages/gammapy/stats/fit_statistics.py:208: RuntimeWarning: invalid value encountered in multiply D = np.sqrt(C**2 + 4 * alpha * (alpha + 1) * n_off * mu_sig) .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 12.858 seconds) .. _sphx_glr_download_tutorials_analysis-2d_ring_background.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v1.3?urlpath=lab/tree/notebooks/1.3/tutorials/analysis-2d/ring_background.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: ring_background.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: ring_background.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: ring_background.zip ` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_