\n",
"\n",
"**This is a fixed-text formatted version of a Jupyter notebook**\n",
"\n",
"- Try online[![Binder](https://static.mybinder.org/badge.svg)](https://mybinder.org/v2/gh/gammapy/gammapy-webpage/v0.19?urlpath=lab/tree/tutorials/analysis/2D/ring_background.ipynb)\n",
"- You may download all the notebooks in the documentation as a\n",
"[tar file](../../../_downloads/notebooks-0.19.tar).\n",
"- **Source files:**\n",
"[ring_background.ipynb](../../../_static/notebooks/ring_background.ipynb) |\n",
"[ring_background.py](../../../_static/notebooks/ring_background.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Ring background map \n",
"\n",
"## Context:\n",
"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.](https://arxiv.org/abs/astro-ph/0610959)\n",
"\n",
"## Objective:\n",
"Create an excess (gamma-ray events) map of MSH 15-52 as well as a significance map to determine how solid the signal is.\n",
"\n",
"## Proposed approach:\n",
"\n",
"The analysis workflow is roughly\n",
" - Compute the sky maps keeping each observation separately using the `Analysis` class\n",
" - Estimate the background using the `RingBackgroundMaker`\n",
" - Compute the correlated excess and significance maps using the `CorrelatedExcessMapEstimator`\n",
" \n",
"The normalised background thus obtained can be used for general modelling and fitting."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup\n",
"As usual, we'll start with some general imports..."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:56.626948Z",
"iopub.status.busy": "2021-11-22T21:06:56.626104Z",
"iopub.status.idle": "2021-11-22T21:06:57.187083Z",
"shell.execute_reply": "2021-11-22T21:06:57.187267Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import astropy.units as u\n",
"from astropy.coordinates import SkyCoord\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"from regions import CircleSkyRegion\n",
"from scipy.stats import norm\n",
"\n",
"import logging\n",
"\n",
"log = logging.getLogger(__name__)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now let's import gammapy specific classes and functions"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.189766Z",
"iopub.status.busy": "2021-11-22T21:06:57.189226Z",
"iopub.status.idle": "2021-11-22T21:06:57.381499Z",
"shell.execute_reply": "2021-11-22T21:06:57.381669Z"
}
},
"outputs": [],
"source": [
"from gammapy.analysis import Analysis, AnalysisConfig\n",
"from gammapy.makers import RingBackgroundMaker\n",
"from gammapy.estimators import ExcessMapEstimator\n",
"from gammapy.maps import Map\n",
"from gammapy.datasets import MapDatasetOnOff"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating the config file\n",
"Now, we create a config file for out analysis. You may load this from disc if you have a pre-defined config file.\n",
"\n",
"In this example, we will use a few HESS runs on the pulsar wind nebula, MSH 1552"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.384028Z",
"iopub.status.busy": "2021-11-22T21:06:57.383688Z",
"iopub.status.idle": "2021-11-22T21:06:57.385110Z",
"shell.execute_reply": "2021-11-22T21:06:57.384925Z"
}
},
"outputs": [],
"source": [
"# source_pos = SkyCoord.from_name(\"MSH 15-52\")\n",
"source_pos = SkyCoord(228.32, -59.08, unit=\"deg\")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.387876Z",
"iopub.status.busy": "2021-11-22T21:06:57.387568Z",
"iopub.status.idle": "2021-11-22T21:06:57.388800Z",
"shell.execute_reply": "2021-11-22T21:06:57.388973Z"
}
},
"outputs": [],
"source": [
"config = AnalysisConfig()\n",
"# Select observations - 2.5 degrees from the source position\n",
"config.observations.datastore = \"$GAMMAPY_DATA/hess-dl3-dr1/\"\n",
"config.observations.obs_cone = {\n",
" \"frame\": \"icrs\",\n",
" \"lon\": source_pos.ra,\n",
" \"lat\": source_pos.dec,\n",
" \"radius\": 2.5 * u.deg,\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.391807Z",
"iopub.status.busy": "2021-11-22T21:06:57.391502Z",
"iopub.status.idle": "2021-11-22T21:06:57.392629Z",
"shell.execute_reply": "2021-11-22T21:06:57.392809Z"
}
},
"outputs": [],
"source": [
"config.datasets.type = \"3d\"\n",
"config.datasets.geom.wcs.skydir = {\n",
" \"lon\": source_pos.ra,\n",
" \"lat\": source_pos.dec,\n",
" \"frame\": \"icrs\",\n",
"} # The WCS geometry - centered on MSH 15-52\n",
"config.datasets.geom.wcs.width = {\"width\": \"3 deg\", \"height\": \"3 deg\"}\n",
"config.datasets.geom.wcs.binsize = \"0.02 deg\"\n",
"\n",
"# Cutout size (for the run-wise event selection)\n",
"config.datasets.geom.selection.offset_max = 3.5 * u.deg\n",
"\n",
"# We now fix the energy axis for the counts map - (the reconstructed energy binning)\n",
"config.datasets.geom.axes.energy.min = \"0.5 TeV\"\n",
"config.datasets.geom.axes.energy.max = \"5 TeV\"\n",
"config.datasets.geom.axes.energy.nbins = 10\n",
"\n",
"# We need to extract the ring for each observation separately, hence, no stacking at this stage\n",
"config.datasets.stack = False"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.397526Z",
"iopub.status.busy": "2021-11-22T21:06:57.397186Z",
"iopub.status.idle": "2021-11-22T21:06:57.399320Z",
"shell.execute_reply": "2021-11-22T21:06:57.399505Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AnalysisConfig\n",
"\n",
" general:\n",
" log: {level: info, filename: null, filemode: null, format: null, datefmt: null}\n",
" outdir: .\n",
" observations:\n",
" datastore: $GAMMAPY_DATA/hess-dl3-dr1\n",
" obs_ids: []\n",
" obs_file: null\n",
" obs_cone: {frame: icrs, lon: 228.32 deg, lat: -59.08 deg, radius: 2.5 deg}\n",
" obs_time: {start: null, stop: null}\n",
" required_irf: [aeff, edisp, psf, bkg]\n",
" datasets:\n",
" type: 3d\n",
" stack: false\n",
" geom:\n",
" wcs:\n",
" skydir: {frame: icrs, lon: 228.32 deg, lat: -59.08 deg}\n",
" binsize: 0.02 deg\n",
" width: {width: 3.0 deg, height: 3.0 deg}\n",
" binsize_irf: 0.2 deg\n",
" selection: {offset_max: 3.5 deg}\n",
" axes:\n",
" energy: {min: 0.5 TeV, max: 5.0 TeV, nbins: 10}\n",
" energy_true: {min: 0.5 TeV, max: 20.0 TeV, nbins: 16}\n",
" map_selection: [counts, exposure, background, psf, edisp]\n",
" background:\n",
" method: null\n",
" exclusion: null\n",
" parameters: {}\n",
" safe_mask:\n",
" methods: [aeff-default]\n",
" parameters: {}\n",
" on_region: {frame: null, lon: null, lat: null, radius: null}\n",
" containment_correction: true\n",
" fit:\n",
" fit_range: {min: null, max: null}\n",
" flux_points:\n",
" energy: {min: null, max: null, nbins: null}\n",
" source: source\n",
" parameters: {selection_optional: all}\n",
" excess_map:\n",
" correlation_radius: 0.1 deg\n",
" parameters: {}\n",
" energy_edges: {min: null, max: null, nbins: null}\n",
" light_curve:\n",
" time_intervals: {start: null, stop: null}\n",
" energy_edges: {min: null, max: null, nbins: null}\n",
" source: source\n",
" parameters: {selection_optional: all}\n",
" \n"
]
}
],
"source": [
"print(config)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Getting the reduced dataset\n",
"We now use the config file to do the initial data reduction which will then be used for a ring extraction"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.401741Z",
"iopub.status.busy": "2021-11-22T21:06:57.401439Z",
"iopub.status.idle": "2021-11-22T21:06:57.439114Z",
"shell.execute_reply": "2021-11-22T21:06:57.439296Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Setting logging config: {'level': 'INFO', 'filename': None, 'filemode': None, 'format': None, 'datefmt': None}\n",
"Fetching observations.\n",
"No HDU found matching: OBS_ID = 20136, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20137, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20151, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20282, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20283, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20301, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20302, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20303, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20322, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20323, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20324, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20325, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20343, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20344, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20345, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20346, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20365, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20366, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20367, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"No HDU found matching: OBS_ID = 20368, HDU_TYPE = rad_max, HDU_CLASS = None\n",
"Number of selected observations: 20\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 33.9 ms, sys: 2.35 ms, total: 36.3 ms\n",
"Wall time: 36.2 ms\n"
]
}
],
"source": [
"%%time\n",
"# create the config\n",
"analysis = Analysis(config)\n",
"\n",
"# for this specific case,w e do not need fine bins in true energy\n",
"analysis.config.datasets.geom.axes.energy_true = (\n",
" analysis.config.datasets.geom.axes.energy\n",
")\n",
"\n",
"# `First get the required observations\n",
"analysis.get_observations()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.443745Z",
"iopub.status.busy": "2021-11-22T21:06:57.443463Z",
"iopub.status.idle": "2021-11-22T21:06:57.444856Z",
"shell.execute_reply": "2021-11-22T21:06:57.445034Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"AnalysisConfig\n",
"\n",
" general:\n",
" log: {level: INFO, filename: null, filemode: null, format: null, datefmt: null}\n",
" outdir: .\n",
" observations:\n",
" datastore: $GAMMAPY_DATA/hess-dl3-dr1\n",
" obs_ids: []\n",
" obs_file: null\n",
" obs_cone: {frame: icrs, lon: 228.32 deg, lat: -59.08 deg, radius: 2.5 deg}\n",
" obs_time: {start: null, stop: null}\n",
" required_irf: [aeff, edisp, psf, bkg]\n",
" datasets:\n",
" type: 3d\n",
" stack: false\n",
" geom:\n",
" wcs:\n",
" skydir: {frame: icrs, lon: 228.32 deg, lat: -59.08 deg}\n",
" binsize: 0.02 deg\n",
" width: {width: 3.0 deg, height: 3.0 deg}\n",
" binsize_irf: 0.2 deg\n",
" selection: {offset_max: 3.5 deg}\n",
" axes:\n",
" energy: {min: 0.5 TeV, max: 5.0 TeV, nbins: 10}\n",
" energy_true: {min: 0.5 TeV, max: 5.0 TeV, nbins: 10}\n",
" map_selection: [counts, exposure, background, psf, edisp]\n",
" background:\n",
" method: null\n",
" exclusion: null\n",
" parameters: {}\n",
" safe_mask:\n",
" methods: [aeff-default]\n",
" parameters: {}\n",
" on_region: {frame: null, lon: null, lat: null, radius: null}\n",
" containment_correction: true\n",
" fit:\n",
" fit_range: {min: null, max: null}\n",
" flux_points:\n",
" energy: {min: null, max: null, nbins: null}\n",
" source: source\n",
" parameters: {selection_optional: all}\n",
" excess_map:\n",
" correlation_radius: 0.1 deg\n",
" parameters: {}\n",
" energy_edges: {min: null, max: null, nbins: null}\n",
" light_curve:\n",
" time_intervals: {start: null, stop: null}\n",
" energy_edges: {min: null, max: null, nbins: null}\n",
" source: source\n",
" parameters: {selection_optional: all}\n",
" \n"
]
}
],
"source": [
"print(analysis.config)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:06:57.446818Z",
"iopub.status.busy": "2021-11-22T21:06:57.446528Z",
"iopub.status.idle": "2021-11-22T21:07:08.075551Z",
"shell.execute_reply": "2021-11-22T21:07:08.075786Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Creating reference dataset and makers.\n",
"Creating the background Maker.\n",
"No background maker set. Check configuration.\n",
"Start the data reduction loop.\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"CPU times: user 9.94 s, sys: 612 ms, total: 10.6 s\n",
"Wall time: 10.6 s\n"
]
}
],
"source": [
"%%time\n",
"# Data extraction\n",
"analysis.get_datasets()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Extracting the ring background\n",
"\n",
"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 `MapDatasetOnOFF` and the `RingBackgroundMaker` classes."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create exclusion mask\n",
"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."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:07:08.079301Z",
"iopub.status.busy": "2021-11-22T21:07:08.078920Z",
"iopub.status.idle": "2021-11-22T21:07:08.309633Z",
"shell.execute_reply": "2021-11-22T21:07:08.309842Z"
}
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASIAAAEHCAYAAAAK30szAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAYo0lEQVR4nO3de5RdVYHn8e/PhEceSFIEGB7Ruihq2/QQ7DSi2DQCg2m1hUF00GYmtOjMOLSgrTM6vXppI06vQbu1XS2u0QVCGhwBaZRIKxITAVsBDY/waLChLcUII5hAQRCJgd/8cXbBtVKpukXVrV236vdZ66x7zj7nnrtP3dxf9j7nnn1lm4iImp5XuwIREQmiUUhS7TpMhdlynLNJr72nCaIdkDQPeLR2PabIltoVmAqSHq9dhynUU+9pgmgHbD9B/j4RUyIftIioTrlq1liyZIn7+/t/o2xgYIBWq1WnQlMoxznzTNdjHRgY8KZNm7ZrAM2tUZnpqL+/n/Xr19euRsSMtnz58hFPoqdrFhHVJYgioroEUURUlyCKiOoSRBFRXYIoIqpLEEVEdQmiiKguQRQR1SWIIqK6BFFEVJcgiojqEkQRUV2CKCKqSxBFRHUJooioLkEUEdUliCKiugRRRFSXIIqI6hJEEVFdgigiqksQRUR1CaKIqC5BFBHVJYgioroEUURUN2VBJOlISYOSbi3Th9vWvU/SnZLukPQlSbuW8n0lrZN0haSFpWwXSZdIulfSjZL62/azUtI9ZVrZVn5N+3YRMb1MdYvoO7aXlemjAJL2A04Hlts+CJgDnFS2Px14D3AucHIpOxV42PaLgU8BZ5f99AEfAV4JHAp8RNLiqTmsiJiI6dI1mwvMkzQXmA/cX8rnAE+XSaXsOGBVmb8MOFqSgNcBa2xvtv0wsAZYUbbbDDzV9aOIiOdk7hS/3qskbaAJmg/YvtP2zyT9NXAf8ARwte2ry/afAS4EBoG3l7L9gJ8C2N4maRDYo7282FjKsH1Cdw8rIiZiKoPoZuCFtrdIej3wVeDA0n06DmgBjwBflnSy7Yts/wQ4Yth+xPY8SvmISivqUZpWl/v6+sZ5OBExXgMDA0h6vCw+Bexu213tmkk6bejkNLDQ9hYA218HdpK0BDgGGLD9kO1fA5cDrx5ltxuBpWX/c4Hdabpez5QX+/NsF287tg3sCyywvaDVaj3Ho4yITrVaLWwvABYCS8vnsLvniGyfM3RyGni6tEKQdGh57U00XbLDJM0v648G7hplt6uBoStiJwLrysF8EzhW0uLSyjq2lI1Wv8eG/hARMXXcGBxansqu2YnAuyVtozkXdFIJgRslXUbTddsG3AJ8fpT9nAdcKOlempbQSQC2N0s6C/hB2e6jtjd351AiYjIpDYLG8uXLvX79+trViJjRli9fzvr167c7nztdLt9HxCyWIIqI6hJEEVFdgigiqksQRUR1CaKIqC5BFBHVJYgioroEUURUlyCKiOoSRBFRXYIoIqpLEEVEdQmiiKguQRQR1SWIIqK6BFFEVJcgiojqEkQRUV2CKCKqSxBFRHUJooioLkEUEdUliCKiugRRRFSXIIqI6hJEEVFdzwWRpF0lfV/SBkl3SjqzlPdJWiPpnvK4uJQfKemCqpWOiFH1XBABTwJH2T4YWAaskHQY8CFgre0DgbVlOSJ6QM8FkRtbyuJOZTJwHLCqlK8Cji/zW4HBqaxjRIxPzwURgKQ5km4FHgTW2L4R2Nv2AwDlca8y/z3bZ1SrbESMaW7tCjwXtp8ClklaBHxF0kHj3Yck0bSU5gD09fVNah0jYnsDAwNIerwsPgXsbts92SIaYvsR4BpgBfBzSfsAlMcHx3iugf2AhbYXtFqt7lY2Imi1WtheACwElpbPYe91zSTtWVpCSJoHHAPcDawGVpbNVgJXjLUv248N/SEiYuqUc73PnLvtxa7ZPsAqSXNogvRS21dKuh64VNKpwH3AW2pWMiI613NBZPs24JARyjcBR099jSJionquaxYRM0+CKCKqSxBFRHUJooioLkEUEdUliCKiugRRRFSXIIqI6hJEEVFdgigiquvoFg9JhwN/CbywPEc0960d0L2qRcRs0em9ZucB7wNuohlDJCJi0nQaRIO2v9HVmkTErNVpEH1b0ieAy2kGrwfA9s1dqVVEm2Ywzc5liKne02kQvbI8Lm8rM3DU5FYnImajjoLI9mu7XZGIIeNtAY31/LSQpr+OLt9L2l3SJyWtL9PfSNq925WL2UXShENoKvcbk6fT7xF9AXgMeGuZHgXO71alImJ26fQc0Ytsv7lt+czyu2IREzKVLZX210p3bXrptEX0hKTXDC2ULzg+0Z0qxWxRs7uUrtr00mmL6N00v5yxO823qjcDp3SrUjHzTYcgkJSW0TTR6VWzW4GDJT2/LD/azUpFxOwyahBJOtn2RZL+bFg5ALY/2cW6RcQsMVaLaEF53G2EdWnTxrhNhy5Zu7b/VCvXZHYbNYhsf67Mfsv2d9vXlRPWERET1ulVs7/rsCxih6Zba6jddK7bbDDWOaJXAa8G9hx2nuj5wJxuViwiZo+xWkQ7AwtpAmu3tulR4MTxvJCkIyUNSrq1TB9uW3eGpDsk3SnpvW3l+0paJ+kKSQtL2S6SLpF0r6QbJfW3bb9S0j1lWtlWfk37dhExvYx1juha4FpJF9j+ySS83ndsv7G9QNJBwLuAQ4GtwFWS/tH2PcDpwHuAA4CTgf8DnAo8bPvFkk4Czgb+g6Q+4CM0IwQYuEnSatsPT0K9I6KLOj1H9EtJn5D09dJCWSdp3STV4beAG2z/0vY24Frg35d1c4CnyzTUiT8OWFXmLwOOVtPBfx2wxvbmEj5rgBVlu81kZMlqeuWm016p50zUaRB9EbgbaAFnAj8GfvAcXu9VkjZI+oak3y5ldwBHSNpD0nzg9cDSsu4zwOeA/wpcVMr2A34KUIJrENijvbzYWMqwfYLt9nURMY10eovHHrbPk3RGW3ft2nG+1s3AC21vkfR64KvAgbbvknQ2TQtmC7AB2AZQuoNHDNvPSP9leZTyiJjmOm0R/bo8PiDpDZIOAfYf60mSThs6OQ0stL0FwPbXgZ0kLSnL59l+he0jaLpR94yy242UFpOkucDu5TnPlBf7A/ePUjdJelTS45IeHxgYGOtwImKCBgYGKJ+5X5bPn6DzIPpYueH1/cAHgHNpftVjVLbPsb3M9jLg6aEXlXRoee1NZXmv8vgC4ATgS6PsdjUwdEXsRGCdm6/FfhM4VtJiSYuBY0vZjupmmuBaaHtBq9Ua63AiYoJarRa2F9DctbG0fA47vun1yjI7CDzXYWNPBN4taRvNECInDVUC+AdJe9C0vE4b40rXecCFku6laQmdVOq4WdJZPHvu6qO2N49WIduDz/FYImICymf/mc9fpz+wuCfNJfb+9ufYfsc4XvgzNCefR1r3++PYz6+At+xg3RdoRpOMiB7S6cnqK4DvAN8il8EjYpJ1GkTzbX+wqzWJiFmr05PVV5ZL7hERk67TIDqDJoyeKJfcHpOUURojYlJ0etVspIHRIiImxVjDgLzM9t2SXjHSets3d6daMZMMfUtjut/HlVEa6xmrRfR+msv2fzPCOgNHTXqNImLWGWsYkHeVx+f6JcaIiDGN1TU7YbT1ti+f3OpExGw0Vtfsj0ZZZyBBFB2zPW3PE+X8UF1jdc3+ZKoqEhGzV0ffI5L0V5IWtS0vlvSxrtUqZizb06r1Md3qM1t1+oXGP7T9yNBCuTs+37SOiEnRaRDNkbTL0IKkecAuo2wfEdGxTm96vQhYK+l8mpPU7+DZAewjxm06nLhOl2z66PQWj49Lug04hmZs6LNs73D0w4hO1AyjhND00mmLCOAuYJvtb0maL2k32491q2IxO7QHQrdDKeEzfXV61exdNL8h9rlStB/Nr3BERExYpyerTwMOp/mpacqvsO7VrUrF7NStS+m5RD/9ddo1e9L21qGmc/kZn7yz0RXDQ2O8XbaETu/ptEV0raQ/B+ZJ+nfAl4Gvda9aETGbdBpEHwIeAm4H/gvwdeAvulWpiHZDXatOp+g9nV6+f1rSV4Gv2n6ou1WKiNlm1BZR+Vnmv5T0C+Bu4IeSHpL04ampXkTMBmN1zd5Lc7Xs92zvYbsPeCVwuKQxf3I6IqITYwXRfwLeZntgqMD2j4CTy7qIiAkbK4h2sv2L4YXlPNFO3alSRMw2YwXR1ue4LiKiY2MF0cHlBxWHT48BvzMVFRxO0q6Svi9pg6Q7JZ1ZyvskrZF0T3lcXMqPlHRBjbpGRGdGDSLbc2w/f4RpN9u1umZPAkfZPhhYBqyQdBjNd53W2j4QWFuWI6IHdPqFxmnDjS1lcacyGTiOZ8dIWgUcX+a3AoNTWceIGJ+eCyIASXMk3Qo8CKyxfSOwt+0HAMrjXmX+e7bPqFbZiBjTeMYjmjZsPwUsKwP6f0XSQePdh5o7KQdp/gbu6+ub3EpGxHYGBgaQ9HhZ3AYssu2ebBENKQP6XwOsAH4uaR+A8vjgGM81sBRYYHtBq9XqbmUjglarhe0FwELgBeVz2HtdM0l7Dv20URnE/xia209WAyvLZiuBK8bal+3BoT9EREydcq73mXO3vdg12wdYJWkOTZBeavtKSdcDl0o6FbgPeEvNSkZE53ouiGzfBhwyQvkm4Oipr1FETFTPdc0iYuZJEEVEdQmiiKguQRQR1SWIIqK6BFFEVJcgiojqEkQRUV2CKCKqSxBFRHUJooioLkEUEdUliCKiugRRRFSXIIqI6hJEEVFdgigiqksQRUR1CaKIqC5BFBHVJYgioroEUURUlyCKiOoSRBFRXYIoIqpLEEVEdQmiiKhuSoNI0pGSbpV0p6Rr28pXSPqhpHslfaitfF9J6yRdIWlhKdtF0iVl2xsl9bdtv1LSPWVa2VZ+Tft2ETG9TFkQSVoEfBZ4k+3fBt5SyucA5wB/CLwceJukl5ennQ68BzgXOLmUnQo8bPvFwKeAs8t++oCPAK8EDgU+Imlx948sIiZqKltEbwcut30fgO0HS/mhwL22f2R7K3AxcFxZNwd4ukwqZccBq8r8ZcDRkgS8Dlhje7Pth4E1wIqy3Wbgqa4dWURMyNwpfK2XADtJugbYDfi07b8H9gN+2rbdRppWDcBngAuBQZogo31729skDQJ77GA/+5XtTujC8UTEJJnKIJoL/C5wNDAPuF7SDTzb0mlnANs/AY4Ytm5H2+9wPyMprahHSr3o6+sbvfYRMWEDAwNIerwsbgMW2XZXu2aSTisnp28F7geusv247V8A1wEH07RclrY9bf+y7Y48s72kucDuNF2vce3HtoEXAAttL2i1WuM8uogYr1arhe0FwEKgv3wOu3uOyPY5tpfZXgZ8Bfh9SXMlzafpft0F/AA4UFJL0s7AScDqUXa7Ghi6InYisK4czDeBYyUtLiepjy1lo9VvcOgPERFTx42Hh5anrGtm+y5JVwG30Zx8Ptf2HQCS/pQmNOYAX7B95yi7Og+4UNK9NC2hk8r+N0s6iybYAD5qe3N3jiYiJpPSIGgsX77c69evr12NiBlt+fLlrF+/frvzuflmdURUlyCKiOoSRBFRXYIoIqpLEEVEdQmiiKguQRQR1SWIIqK6BFFEVJcgiojqEkQRUV2CKCKqSxBFRHUJooioLkEUEdUliCKiugRRRFSXIIqI6hJEEVFdgigiqksQRUR1CaKIqC5BFBHVJYgioroEUURUlyCKiOoSRBFRXc8FkaRdJX1f0gZJd0o6s5T3SVoj6Z7yuLiUHynpgqqVjohR9VwQAU8CR9k+GFgGrJB0GPAhYK3tA4G1ZTkiekDPBZEbW8riTmUycBywqpSvAo4v81uBwamsY0SMT88FEYCkOZJuBR4E1ti+Edjb9gMA5XGvMv8922dUq2xEjGlu7Qo8F7afApZJWgR8RdJB492HJAGPUP4GfX19k1nFiBjBwMAAkh4vi9uARbbdky2iIbYfAa4BVgA/l7QPQHl8cIznGugHFtpe0Gq1ulrXiIBWq4XtBcBC4IDyOey9rpmkPUtLCEnzgGOAu4HVwMqy2UrgirH2ZfvhoT9EREydcq5309ByL3bN9gFWSZpDE6SX2r5S0vXApZJOBe4D3lKzkhHROaVB0FiyZIn7+/vH/byBgQFmWrcux9Q7eu24BgYGvGnTpu16YgmiCZL0eOnzzhg5pt4xU46r584RRcTMkyCauKfLSfOZRLUrELNLumYTJEkz7crbTDymmWqmvFcJooioLl2zDkg6RdIba9cjYqbqxe8RTQpJXwDeCDxo+6C28h8DjwFPAdtsLy+r3ippBfBz22dJWgr8PfBvgKeBz9v+dNt+5gDrgZ/ZfmMpOwV4LfAE8ADNDbsHAW+1vbWLh/uMUY77fcA7aW4gvh34E+Ck2vXt1EjHtaP3aDq8D50o/94+DcwBzrX9v0v5IuBcmjobeAfwUnrgmHbI9qycgCOAVwB3DCv/MbBkWNkpwB+X+UvK4z7AK8r8bsC/AC9ve86fAf8XuHLYft5e5teWxz8HDql53MB+wAAwryxfWupavb4TPK4R36NeOC6a8PlX4ABgZ2DD0L8vmtEl3lnmdwYW9cIxjTbN2q6Z7euAzeN4ytBQIi7Pf8D2zWX+MeAumg80kvYH3kDzv9Zwj5bHh8rjVmCXcVV+AkY57rnAPElzgfnA/aW8an07NdJxjfYeMf2P61DgXts/ctOiuRg4TtLzaUL3PADbW93ccwnT/5h2aNYG0SgMXC3pJkn/uZMnSOoHDgFuLEV/C/wPmu7AtGf7Z8Bf09wa8wAwaPvqurWaXCO8R9PdfsBP25Y3lrIDaILmfEm3SDpXUs9/obF6k6xy87ef7btm+5bHvWiaw0eMsY+FwE3ACWX5jcBny/yRtHXNpss0/LiBxcA6YE+acwtfBU6uXc/JeD9Heo96YaK5V/LctuX/CPwdsJxm+IxXlvJPA2fVru9Ep7SIhrF9f3l8EPgKTRN5RJJ2Av4B+KLty0vx4cCbyknvi4GjJF3U1UpP3DHAgO2HbP8auBx4deU6TYodvEe9YCOwtG15f5ru8kZgo5vBAAEuozk31tMSRG0kLZC029A8cCxwxw62FU0//S7bnxwqt/0/be9vu5/mqtM62yd3vfITcx9wmKT55biOpjmf0tN29B71iB8AB0pqSdqZ5t/Satv/D/ippJeW7Y4G/rlWJSfLrA0iSV8CrgdeKmljGT5kb+CfJG0Avg/8o+2rdrCLw2may0dJurVMr5+Syk/ASMdd/ne9DLiZ5tL984DPV6zmuO3g/ezJ9wjA9jbgT4Fv0vyncKntO8vq9wBflHQbzQ9I/FWVSk6ifLM6IqqbtS2iiJg+EkQRUV2CKCKqSxBFRHUJooioLkEUEdUliCKiugRRjErSU+WLgHdI+lrbj1vuK+myDp6/ZQflx0t6+RjP3VC+qFhNp8cZE5MgirE8YXuZm8HGNgOnQXNPnu0TJ7Df42nGBhqRpN+i+fd5RM27yyfhOKMDCaIYj+t5dsylfkl3lPn5ki6VdJukSyTdKGloZEsk/a/SurlB0t6SXg28CfhEaW29aITXejtwIXB12XZoX6dL+ufyWheXsoWSzpd0eyl/cyk/VtL1km6W9GVJC0v5jyWdWcpvl/SyUv4HbbeC3CJpt2HHuWvb69wi6bWl/BRJl0u6StI9kj4+yX/3ma/27f+ZpvcEbCmPc4AvAyvKcj9lyA3gA8DnyvxBNMNULC/LBv6ozH8c+IsyfwFw4iiv+y/AC2luPF7dVn4/sEuZX1Qezwb+tm2bxcAS4DpgQSn7IPDhMv9j4D1l/r9RhtsAvgYcXuYX0gwW136c7wfOL/Mvo7lZeFea0RF/BOxeln8CLK393vXSlBZRjGWepFuBTUAfsGaEbV5DM+QJtu8AbmtbtxW4sszfRPPBHpWk3wMesv0TYC3wCkmLy+rbaG74PJkm8KAZxuScoefbfhg4jKbr991S/5U0wTZkaEiQ9jp9F/ikpNNpQm4bv+k1NK00bN9NEzgvKevW2h60/Suau+FfSHQsQRRjecL2MpoP1s6Uc0TDjPaDjL92aULQ/CBBJz/Y8DbgZWVMp38Fng+8uax7A03o/C5wUxnaVpQhfIfVaY2b81vLbL/c9qlt658cXic3g9O/E5gH3DDUZevwOJ9sm+/0OKNIEEVHbA8CpwMfKIONtfsn4K0A5UrY73Swy8doBrT/DZKeRzM64b+13e9mXKfjgLeVdUttf5tmKN5FNF2oq2mGzBjax2LgBuBwSS8uZfMlvYRRSHqR7dttn03zCyzDg+g64I/Lti8BXgD8sINjjTEkiKJjtm+hGT73pGGrPgvsWcbH+SBN92mQ0V0M/Pdy0rf9ZPURND/B9LO2sutouln7ARdJuh24BfiUm4HjPwYsLl8x2AC81vZDNOduvlTqdQPbB8tw723bxxPAN0Y4zjnl9S8BTrH95PCdxPhlPKKYMDW/4baT7V+VUFkLvMS98HtaMS2kHxuTYT7w7dJlE/DuhFCMR1pEEVFdzhFFRHUJooioLkEUEdUliCKiugRRRFT3/wFrhbxIzO/7rwAAAABJRU5ErkJggg==\n",
"text/plain": [
"