\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/api/mask_maps.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",
"[mask_maps.ipynb](../../_static/notebooks/mask_maps.ipynb) |\n",
"[mask_maps.py](../../_static/notebooks/mask_maps.py)\n",
"
\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Mask maps\n",
"\n",
"## Prerequisites\n",
"\n",
"- Understanding of basic analyses in 1D or 3D.\n",
"- Usage of `~regions` and catalogs, see the [catalog notebook](catalog.ipynb). \n",
"\n",
"## Context\n",
"\n",
"There are two main categories of masks in Gammapy for different use cases.\n",
"- Fitting often requires to ignore some parts of a reduced dataset, e.g. to restrict the fit to a specific energy range or to ignore parts of the region of interest that the user does not want to model, or both. \n",
"Gammapy's `Datasets` therefore contain a `mask_fit` sharing the same geometry as the data (i.e. `counts`).\n",
"- During data reduction, some background makers will normalize the background model template\n",
"on the data themselves. To limit contamination by real photons, one has to exclude parts of the \n",
"field-of-view where signal is expected to be large. To do so, one needs to provide an exclusion mask. \n",
"The latter can be provided in a different geometry as it will be reprojected by the `Makers`.\n",
"\n",
"We explain in more details these two types of masks below:\n",
"\n",
"### Masks for fitting: `mask_fit`\n",
"\n",
"The region of interest used for the fit can defined through the dataset `mask_fit` attribute.\n",
"The `mask_fit` is a map containing boolean values where pixels used in the fit are stored as True. \n",
"\n",
"A spectral fit (1D or 3D) can be restricted to a specific energy range where e.g. the background\n",
"is well estimated or where the number of counts is large enough.\n",
"Similarly, 2D and 3D analyses usually require to work with a wider map than the region of interest so sources laying outside but reconstructed inside because of the PSF are correctly taken into account. Then the `mask_fit` have to include a margin that take into account the PSF width. We will show an example in the boundary mask sub-section.\n",
"\n",
"The `mask_fit` also can be used to exclude sources or complex regions for which we don't have good enough models. In that case the masking is an extra security, it is prefereable to include the available models even if the sources are masked and frozen.\n",
"\n",
"Note that a dataset contains also a `mask_safe` attribute that is created and filled during \n",
"data reduction. It is not to be modified directly by users. The `mask_safe` is defined only from \n",
"the options passed to the `~gammapy.makers.SafeMaskMaker` (More details [here](https://docs.gammapy.org/dev/makers/index.html#safe-data-range-handling)).\n",
"\n",
"### Exclusion masks\n",
"\n",
"Background templates stored in the DL3 IRF are often not reliable enough to be used without some corrections. A set of common techniques to perform background or normalisation from the data is implemented in gammapy: reflected regions for 1D spectrum analysis, field-of-view (FoV) background or ring background for 2D and 3D analyses.\n",
"\n",
"To avoid contamination of the background estimate from gamma-ray bright regions these methods require to exclude those regions from the data used for the estimation. To do so, we use exclusion masks. They are maps containing boolean values where excluded pixels are stored as False. \n",
"\n",
"\n",
"## Proposed approach\n",
"\n",
"Even if the use cases for exclusion masks and fit masks are different, the way to create these masks is exactly the same, so in the following we show how to work with masks in general:\n",
"- Creating masks from scratch\n",
"- Combining multiple masks\n",
"- Extending and reducing an existing mask\n",
"- Reading and writing masks\n",
"\n",
"\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.240600Z",
"iopub.status.busy": "2021-11-22T21:09:28.235980Z",
"iopub.status.idle": "2021-11-22T21:09:28.410758Z",
"shell.execute_reply": "2021-11-22T21:09:28.410934Z"
}
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.413301Z",
"iopub.status.busy": "2021-11-22T21:09:28.412930Z",
"iopub.status.idle": "2021-11-22T21:09:28.932760Z",
"shell.execute_reply": "2021-11-22T21:09:28.932955Z"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"from astropy.coordinates import SkyCoord, Angle\n",
"import astropy.units as u\n",
"from regions import CircleSkyRegion, Regions\n",
"from gammapy.maps import Map, WcsGeom, MapAxis\n",
"from gammapy.catalog import CATALOG_REGISTRY\n",
"from gammapy.datasets import Datasets\n",
"from gammapy.estimators import ExcessMapEstimator\n",
"from gammapy.modeling.models import FoVBackgroundModel"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Creating a mask for fitting\n",
"\n",
"One can build a `mask_fit` to restrict the energy range of pixels used to\n",
"fit a `Dataset`. The mask being a `Map` it needs to use the same geometry\n",
"(i.e. a `Geom` object) as the `Dataset` it will be applied to. \n",
"\n",
"We show here how to proceed on a `MapDataset` taken from Fermi data used in the 3FHL catalog. The dataset is already in the form of a `Datasets` object. We read it from disk.\n"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.935100Z",
"iopub.status.busy": "2021-11-22T21:09:28.934813Z",
"iopub.status.idle": "2021-11-22T21:09:28.946171Z",
"shell.execute_reply": "2021-11-22T21:09:28.946358Z"
}
},
"outputs": [],
"source": [
"filename = \"$GAMMAPY_DATA/fermi-3fhl-crab/Fermi-LAT-3FHL_datasets.yaml\"\n",
"datasets = Datasets.read(filename=filename)\n",
"dataset = datasets[\"Fermi-LAT\"]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can check the default energy range of the dataset. In the absence of a `mask_fit` it is equal to the safe energy range."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.948216Z",
"iopub.status.busy": "2021-11-22T21:09:28.947909Z",
"iopub.status.idle": "2021-11-22T21:09:28.989634Z",
"shell.execute_reply": "2021-11-22T21:09:28.989871Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"HDU 'MASK_FIT' not found\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fit range : (WcsNDMap\n",
"\n",
"\tgeom : WcsGeom \n",
" \taxes : ['lon', 'lat']\n",
"\tshape : (50, 40)\n",
"\tndim : 2\n",
"\tunit : TeV\n",
"\tdtype : float64\n",
", WcsNDMap\n",
"\n",
"\tgeom : WcsGeom \n",
" \taxes : ['lon', 'lat']\n",
"\tshape : (50, 40)\n",
"\tndim : 2\n",
"\tunit : TeV\n",
"\tdtype : float64\n",
")\n"
]
}
],
"source": [
"print(f\"Fit range : {dataset.energy_range}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Create a mask in energy\n",
"\n",
"We show first how to use a simple helper function `~gammapy.maps.Geom.energy_range()`.\n",
"\n",
"We obtain the `Geom` that is stored on the `counts` map inside the `Dataset`\n",
"and we can directly create the `Map`.\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.992000Z",
"iopub.status.busy": "2021-11-22T21:09:28.991696Z",
"iopub.status.idle": "2021-11-22T21:09:28.992913Z",
"shell.execute_reply": "2021-11-22T21:09:28.993092Z"
}
},
"outputs": [],
"source": [
"mask_energy = dataset.counts.geom.energy_mask(10 * u.GeV, 700 * u.GeV)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now set the dataset `mask_fit` attribute. \n",
"\n",
"And we check that the total fit range has changed accordingly. The bin edges closest to requested range provide the actual fit range. "
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:28.995454Z",
"iopub.status.busy": "2021-11-22T21:09:28.995129Z",
"iopub.status.idle": "2021-11-22T21:09:28.996787Z",
"shell.execute_reply": "2021-11-22T21:09:28.997060Z"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Fit range : (WcsNDMap\n",
"\n",
"\tgeom : WcsGeom \n",
" \taxes : ['lon', 'lat']\n",
"\tshape : (50, 40)\n",
"\tndim : 2\n",
"\tunit : TeV\n",
"\tdtype : float64\n",
", WcsNDMap\n",
"\n",
"\tgeom : WcsGeom \n",
" \taxes : ['lon', 'lat']\n",
"\tshape : (50, 40)\n",
"\tndim : 2\n",
"\tunit : TeV\n",
"\tdtype : float64\n",
")\n"
]
}
],
"source": [
"dataset.mask_fit = mask_energy\n",
"print(f\"Fit range : {dataset.energy_range}\")"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Mask some sky regions\n",
"\n",
"One might also exclude some specific part of the sky for the fit. For instance, if one wants not to model a specific source in the region of interest, or if one want to reduce the region of interest in the dataset `Geom`.\n",
"\n",
"In the following we restrict the fit region to a square around the Crab nebula. **Note**: the dataset geometry is aligned on the galactic frame, we use the same frame to define the box to ensure a correct alignment. \n",
"We can now create the map. We use the `WcsGeom.region_mask` method putting all pixels outside the regions to False (because we only want to consider pixels inside the region. For convenience we can directly pass a ds9 region string to the method:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:29.007283Z",
"iopub.status.busy": "2021-11-22T21:09:29.006912Z",
"iopub.status.idle": "2021-11-22T21:09:29.011255Z",
"shell.execute_reply": "2021-11-22T21:09:29.011483Z"
}
},
"outputs": [],
"source": [
"regions = \"galactic;box(184.55, -5.78, 3.0, 3.0)\"\n",
"mask_map = dataset.counts.geom.region_mask(regions)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now combine this mask with the energy mask using the logical and operator"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:29.013566Z",
"iopub.status.busy": "2021-11-22T21:09:29.013279Z",
"iopub.status.idle": "2021-11-22T21:09:29.014635Z",
"shell.execute_reply": "2021-11-22T21:09:29.014807Z"
}
},
"outputs": [],
"source": [
"dataset.mask_fit &= mask_map"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's check the result and plot the full mask."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"execution": {
"iopub.execute_input": "2021-11-22T21:09:29.017125Z",
"iopub.status.busy": "2021-11-22T21:09:29.016836Z",
"iopub.status.idle": "2021-11-22T21:09:29.289022Z",
"shell.execute_reply": "2021-11-22T21:09:29.289267Z"
}
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/Users/adonath/software/mambaforge/envs/gammapy-dev/lib/python3.9/site-packages/astropy/visualization/wcsaxes/core.py:211: MatplotlibDeprecationWarning: Passing parameters norm and vmin/vmax simultaneously is deprecated since 3.3 and will become an error two minor releases later. Please pass vmin/vmax directly to the norm when creating it.\n",
" return super().imshow(X, *args, origin=origin, **kwargs)\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA0MAAACyCAYAAACA5hTOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAbnklEQVR4nO3deZRtZX3m8e/DReEyqoBxilxMjBDTgiLOKKJGMXEFGhMvjizHuDDd6nJOOssYE1uUaEewTWfQG6cYBY2aOMQBRQUVCQgkcWhxtqNBGe4FB+DXf5xdsG9Z59Q5p2qfU+fu72etvWoP797v+9Z+6t56aw8nVYUkSZIk9c1u826AJEmSJM2DgyFJkiRJveRgSJIkSVIvORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSLy3UYCjJ15Ncm2R7azp93u1aTZItST6e5Jok/57koSPKJskrk1zeTKcmyQrljm59D3YkqWXflzuusM99m7L7rrDtX5I8y/4NbbfZu6nchjk3E/bvwU3ZK5N8fUS5DdO/Zv8+ZG/hzs2E/Xt+kkuSXJ3ksiTP3+j9a/bvQ/YW7txM2L8PLGvnT5NcvJH71+zfh+wt1LlJcuskb0/y3Qz+rf50knuPKJ8syu8UVbUwE/B14KEzqGf3dT7eucCfAZuBE4ErgIOGlH0G8CXgDsDtgX8FfneV428Bapx2N8c+edm6XwN+Ahxg/8zeIp2bCft3L+AJwNOBr495fLM3m+wt3LmZsH8vAO4B7A7cBfgGsHUj969H2Vu4czNJ/1bY92zgDzdy//qSvUU7N8CdgOcCtwU2Mfj3+j+BfYaUX5jfKToN2XpPo344gJOBTwGvBn4EXAYc19q+P/DXwPeA7wAvBza19v008Brgh822A4D3AVcBn2/WfaopfwZw2rL63wc8e4V2/UpzYvZtrTtnWCCAzwBPby0/BThvkvCs0teXAB9btv+pwFlTnpNdun9mb+Oem0n71yrzUKb8hdvsdZO9RTs30/avVfbPgddt1P71MXuLcm7W0r+m3dcDh2zU/vU1e4t0bpYd7yrgyCHbFuZ3ioW6TW4M92YwUjyQwTfkr1uX5LYB1wG/DNwd+HXgqcv2/Rpwa+BPGPwA7ABuAzypmWgd66QkuwEkORB4CPD2Fdp0V+BrVXV1a91FzfqV3LXZPk7ZYUb19c3A0UuXHZs+PBb42wnraLd3V+7fuMzeTe3bqNlbD2avm+yth4XIXnNujgYunbBOs9dx9hbk3KzlZ+uJwDlVddmEdZq97v/dW7hzk+QI4ObAV4cUWZzfKdZrdDiLicFfCrYzuOy4ND2tbhrtf7VVdi8Go8vbAL/AYLS+ubX9JODjrX2/2dq2CfgZcJfWuhv/UtAs/xvwsGb+WcA/DWnzE1g2Embww/emIeWvBw5tLd+56UfGGUmv1tdm+SPAS5r5hzG4zHmzKc/JLt0/s7dxz82k/WuVmerqg9nrLnuLdm6m7V9T7o8Y/FKwx0btXx+ztyjnZo39+yrLbivaaP3rcfYW5tw0x9oPuBh48YgyC/M7xSJeGTq+qm7Rmv6yte3/Lc1U1TXN7D7AwcDNgO8luSLJFcBfMPirwJJvteYPYnAivjVkOwxGq49v5h/PYIS6ku0MQtO2H3D1CmVXKr8fsL2aMz2Gcfq6jcFfIWDww/u2qvrZSgdb9oDf48Zo71KbF6J/EzJ7o2307K2V2esue2u1ENlrHup9IvAbVfWTVfrUZvY6zt5GOTcdZu8BDAYL71qlP8uZve6ztyHOzRjZWyq3mcGtgudV1StGtG9D/U4xyu6T7rCgvsVgdHlgVV03pEz75PyAwWW5OwBfbtb94rLybwEuSXI4cBjwniHHvRS4U5J966ZLp4cDbxtR/nDgc62yk1yyH6evZwGvT/Jg4L8Cxww7WFUdt0p9C92/GTB7O5tn9tbK7HWXvbXa8NlL8mTgRcADq+rbqxx/ObPXYfY20rnp8N+9JzF4lmL7KuWWM3vd/7u3Ic7NGNkjyR4M+v8dBi9IGGVD/U4x0lovlc1yYowH6patK+CXm/l/AP4Xg5HpbsAvAQ8ase87GAR4L+BQ4JsrlPln4IvA36zS7vMYPOi3J3ACo99s87sMLsneHrgdg+BM9PaNUX1t7fPG5vt56Tqcl126f2Zv456bCfu3W1PuOAZvjNoTuPlG7l+Psrdw52bC/j2OwV+yD5vg+GZvNtlbuHMzSf+a8pubMscuQv/6kr1FOzcMrtK8j8FgaJy3vS3M7xRrDuwsp6az1zK49LY0vXtEwNs/HPsD/xv4NnAl8C80r88csu9BwD9y09tFXgl8dFmZxzd1PHiMk3t20/Yv0foBZ/Cw5vbWchg8DPjDZjqVEfdXDgnP0L629jmm2eeF63Bedun+mb2Ne24m7N9Sve3p7I3cvx5lb+HOzYT9u4zBcwntc/iGjdy/HmVv4c7NJP1r1p3E4I8MI/893yj960v2Fu3cAA9qjnPNsvNy9JCfrYX5nSLNgbSKJK8EblNVT2qteyCDy6dbquqGuTVOuzSzp3kxe5oXs6d5MXv9s4gvUJiJJIcmuVsG7sXg/ejvbm2/GfDfgb/yB0PryexpXsye5sXsaV7MnhwMDbcvgwezdgB/D5zG4N5FkhzG4B7P2wKvnU/ztAsze5oXs6d5MXuaF7PXc94mJ0mSJKmXvDIkSZIkqZccDEmSJEnqpYX70NUk3tenNbnVrW5Vl19++cR/CDB7WqtpswfmT2tj9jQvZk/zNE7+vDKk3jnkkEMy7zaon8ye5sXsaV7MnuZpnPzNfDCU5Kgk1yd5dGvd1iQXJHn2rNsjSZIkqZ9mOhhKsonBJ/t+aNmmrcBRwH2S7DPLNkmSJEnqp1lfGfo94Ezg+8vWL13Cqta8JEmSJHVmZoOhJLcHTgDesMLms4DzgfOr6upZtUmSJElSf83ybXKvBV5YVdcnO1/8qaptwLYZtkWSJElSz3U6GEpyCvC0ZnF/4O+agdCBwCOTXFdV7xnjODs6a6Q0gtnTPJk/zYvZ07yYPc1ap4OhqjoDOGP5+iRvAt4/zkCoOc7erX1957xmxuxpnsyf5sXsaV7MnmbNzxmSJEmS1EuzfGboRlV18jzqlSRJkqQlXhmSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSL83lBQqSxlPlW0W7svzDn/XzzF93zN9oZq87Zk/amVeGJEmSJPXSWIOhJAcneWgzvznJvt02S5IkSZK6tepgKMnTgHcBf9GsugPwng7bJEmSJEmdG+fK0CnA/YGrAKrqK8Ctu2yUJEmSJHVtnMHQT6rqp0sLSXYHfLJRkiRJ0kIbZzD0iSQvATYneRjwTuB93TZLkiRJkro1zmDoRcAPgIuBZwD/BPxBl42SJEmSpK6t+jlDVXUD8JfNJEmSJEm7hKGDoSQXM+LZoKq6WyctkiRJkqQZGHVl6Debr6c0X9/cfH0ccE1nLZIkSZKkGRg6GKqqbwAkuX9V3b+16UVJPg28rOvGSZIkSVJXxnmBwt5JHrC0kOR+wN7dNUmSJEmSurfqCxSApwB/k2T/ZvkK4MmdtUiSJEmSZmCct8l9ATg8yX5AqurK7pslSZIkSd1adTCU5A+XLQNQVT4zJEmSJGlhjfPM0I7WdD1wHLBl0oqSHJPkyiQXNtMftrZtTXJBkmdPelxJkiRJmsY4t8md1l5O8mrgvVPWd05V/eYK67cCRwFvTbJPVW2f8viSJEmSNJZxrgwttxdwp3VuR5qv1ZqXJEmSpM6M88zQxQwGKQCbgIOAP56yvvsmuQj4LvC8qrq0WX8WcD7wlqq6espjS5IkSdLYxnm1dvu2tuuA/6iq66ao6wLg4KranuSRwHuAOwNU1TZg27Adk+yYoj5pzcye5sn8aV7MnubF7GnWxrlN7uVV9Y1m+k5VXZfkzeMcPMkpSy9MAG58Fqiq/gm4WZIDxzlOVe29NI1TXlovZk/zZP40L2ZP82L2NGvjDIbu2l5Isjtw5DgHr6ozquqIqjoCuCHNe7mT3Kup+/LJmitJkiRJ62PobXJJXgy8BNic5Kql1cBPgf8zRV2PBp6Z5DrgWmBrVdUq+0iSJElSJ4YOhqrqFcArkryiql681oqq6nTg9LUeR5IkSZLWw6grQ4dW1b8D70xyj+Xbq+qCTlsmSZIkSR0a9Ta55wJPB05bYVsBx3bSIkmSJEmagVG3yT29mT2uqn7c3pZkz05bJUmSJEkdG+dtcp8Zc50kSZIkLYxRzwzdBrg9g7fJ3Z3Bm+QA9gP2mkHbJEmSJKkzo54ZejhwMnAH4M9a669m8MptSZIkSVpYo54Z2gZsS3JiVZ05wzZJkiRJUudGXRkCoKrOTPIbwF2BPVvrX9ZlwyRJkiSpS6u+QCHJG4DHAL/H4Lmh3wYO7rhdkiRJktSpcd4md7+qeiLwo6r6I+C+wC922yxJkiRJ6tY4g6Frm6/XJLkd8DPgkO6aJEmSJEndW/WZIeD9SW4BvAq4ACjgL7tslCRJkiR1bZwXKPxxM3tmkvczeInCoZ22SpIkSZI6Ns5tcjeqqp9U1ZXAOztqjyRJkiTNxESDoZasayskSZIkacamHQzVurZCkiRJkmZs6DNDSd7HyoOeAAd01iJJkiRJmoFRL1B49ZTbJEmSJGnDGzoYqqpPzLIhkiRJkjRL0z4zJEmSJEkLbaaDoSTHJLkwyaVJPtFavzXJBUmePcv2SJIkSeqvVT90db0kuQXweuARVfXNJLdubd4KHAW8Nck+VbV9Vu2SJEmS1E+rXhlK8s/NQGZp+ZZJPjRFXY8FzqqqbwJU1ffb1TRfCz/DSJIkSdIMjHOb3IFVdcXSQlX9CLj18OJD/QpwyyRnJ/lCkie2tp0FnA+cX1VXT3FsSZIkSZrIOLfJ3ZDkjktXdJIczHQfuro7cCTwEGAzcG6S86rqy1W1Ddg2bMckO6aoT1ozs6d5Mn+aF7OneTF7mrVxBkO/D3yq9cKDBwJPH+fgSU4BntYs/j3wwaraAexI8kngcODLqx2nqvZuHXOagZg0FbOneTJ/mhezp3kxe5q1VQdDVfXBJPcA7sPgeZ7nVNV/jnPwqjoDOAMgyWHA6Ul2B24O3Bt4zbQNlyRJkqS1GDoYSnJoVf17MxAC+G7z9Y7NbXMXTFJRVf1bkg8CXwRuAP6qqi6ZqtWSJEmStEajrgw9l8HtcKetsK2AYyetrKpeBbxq0v0kSZIkab0NHQxV1dJzQcdV1Y/b25Ls2WmrJEmSJKlj47xa+zNjrpMkSZKkhTHqmaHbALcHNie5Ozd9GOp+wF4zaJskSZIkdWbUM0MPB04G7sDguaGlwdBVwEu6bZYkSZIkdWvUM0PbgG1JTqyqM2fYJkmSJEnq3DjPDB2Z5BZLC0lumeTl3TVJkiRJkro3zmDouKq6Ymmhqn4EPLKzFkmSJEnSDIwzGNqUZI+lhSSbgT1GlJckSZKkDW/UCxSWvAX4aJI3Mviw1ScD2zptlSRJkiR1bNXBUFWdmuRi4CEM3ij3x1X1oc5bJkmSJEkdGufKEFX1AeADHbdFkiRJkmZm1WeGktwnyeeTbE/y0yTXJ7lqFo2TJEmSpK6M8wKF04GTgK8Am4GnAq/rslGSJEmS1LVxb5P7apJNVXU98MYkn+m4XZIkSZLUqXEGQ9ckuTlwYZJTge8Be3fbLEmSJEnq1ji3yT0B2AQ8C9gB/CJwYpeNkiRJkqSujfNq7W80s9cCf9RtcyRJkiRpNoYOhprPFqph26vqbp20SJIkSZJmYNSVod+cWSskSZIkacaGDoZat8etiyTPBx7Xqvcw4KCq+mGSrcALgL+tqteuZ72SJEmStJKZfehqVb2qqo6oqiOAFwOfqKofNpu3AkcB90myz6THliRJkqRJzetDV08C3t5aTvO1WvOSJEmS1JlxBkNU1VeBTVV1fVW9EXjwtBUm2Qt4BHBma/VZwPnA+VV19bTHliRJkqRxzeNDVx8FfLp1ixxVtQ3YtoZjSpIkSdJExv3Q1d2Y4kNXk5yS5MJmul2zeis73yI3znF2LE2T7CetldnTPJk/zYvZ07yYPc3aJB+6+mMm/NDVqjoDOGNpOcn+wIOAx094nBuvRCUZ+tlH0noze5on86d5MXuaF7OnWRt6ZSjJbyU5pbX82SRfa6ZHT1nfCcCHq8rRviRJkqS5GnWb3AuA97aW92Dw+utjgGdOU1lVvamqtk6zryRJkiStp1G3yd28qr7VWv5UVV0OXJ5kLS9QkCRJkqS5G3Vl6Jbthap6VmvxoG6aI0mSJEmzMWow9NkkT1u+MskzgM911yRJkiRJ6t6o2+SeA7wnyWOBC5p1RzJ4duj4jtslSZIkSZ0aOhiqqu8D90tyLHDXZvU/VtXHZtIySZIkSerQOJ8z9DHAAZA0B0nm3QT1mPnTvJg9SbMy6pkhSZIkSdplORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYEiSJElSLzkYkiRJktRLDoYkSZIk9ZKDIUmSJEm95GBIkiRJUi85GJIkSZLUSw6GJEmSJPWSgyFJkiRJveRgSJIkSVIv7T7vBkzqgAMOYMuWLVx22WUccsghM63bOneZ+mqa/Zay1xxjl/4e9aXORckezO/fvj7kYB51mr3x7OrnpQ91mr2NWV9f6hwnf6maOqNzlWRHVe1tnYtf5zz6uFZ9+B71oU6zt/Hq60udZm9j1tmHPs6rzrXow/eoD32cV52r8TY5SZIkSb3kYEiSJElSLy3sbXKSJEmStBZeGZIkSZLUSw6GJEmSJPXShhoMJfmbJN9Pcklr3RFJzktyYZLzk9yrWf+4Zt3SdEOSI5ptxzRlT13POpttd0tybpJLk1ycZM8u60yyJcm1rX6+obVPZ/1stt8xyfYkz5u0zgn7eK9W/y5KcsI0fVyrWefP7Jm9Kdts9szeupl19iats9m2pvz1IXtT9HPu+TN7Zm9e2fs5VbVhJuCBwD2AS1rrPgwc18w/Ejh7hf3+C/C11vI7gM3AacCh61Ung89l+iJweLN8ALCp4zq3tMstO04ndba2nwm8E3jepHVO2Me9gN2b+dsC328tj93HRcuf2TN7Zs/s9S1788hfH7K3iPkze2ZvXtlbPm2oK0NV9Ungh8tXA/s18/sD311h15OAt7eWd2v2uwHIOtb568AXq+qiZt/Lq+r6juscpbM6kxwPfA24dJo6J6mvqq6pquua9Xs25Saqbz3MOn9mz+zd2Dizt7zOUczeOurD/7t9yN6kdW6E/Jm9n6tvFLPXpVmNusadWDYyBg4Dvgl8C/gOcPAK+/xf4Ndayw8HLgBOW886gWcDbwY+1Bz/BTOocwuwA/gX4BPA0TOoc2/gXGAf4KXs/JeCseuc5FwC92bwg7gdOGHaPi5a/sye2TN7Zq9v2ZtH/vqQvUXMn9kze/PK3k5tn2VlU34z/xw4sZn/HeAjy8rfG7h4FnUCzwMuAw5kcJnvXOAhHde5B3BAM39kE6r9Oq7z1cDvNPM7/XB0eS6b9YcBnwP27EP+zJ7ZM3tmr2/Zm0f++pC9Rcyf2TN788reTu2YZ+VjfjOvhBs/DynAVcvKvwZ4ySzqBLYCb2qV+x/A82fRz1a5s4F7dtzPc4CvN9MVDC59PmuGffz4tH1ctPyZPbNn9sxe37I3j/z1IXuLmD+zZ/bmlb32tKGeGRriu8CDmvljga8sbUiyG/DbwN/NqM4PAXdLsleS3Zsy/9plnUkOSrKpmb8TcGcG93Z2VmdVHV1VW6pqC/Ba4E+r6vSu6ktySPP9JMnBwF0Y/GBuBLPOn9kze0vMHmZvTvrw/24fsje0zg2cP7OH2Zu5eY7EVhgdvh34HvAz4NvAU4AHAF8ALgI+CxzZKn8McN6M63w8g/scLwFO7bpO4MSmvosY3Ef5qFn0s7XfS5nisumEfXxC08cLmz4e34f8mT2zZ/bMXt+yN4/89SF7i5g/s2f25pW95dPS5StJkiRJ6pVFuE1OkiRJktadgyFJkiRJveRgSJIkSVIvORiSJEmS1EsOhiRJkiT1koMhSZIkSb3kYGiEJL+Q5G1JvpbkC0nOTXLCKvtsSXLJlPWdnOR2reW/SvKrY+57TJL3T1PvmMe/XZJ3NfNHJHnkFMd4aZLnrX/rdj1mb6fjm70ZMns7Hd/szZDZ2+n4Zm+GzN5Ox+9d9hwMDZEkwHuAT1bVnarqSGArcIcOqz0ZuPGHo6qeWlXr9Unva1JV362qRzeLRwAT/3BoPGZvZ2Zvdszezsze7Ji9nZm92TF7O+tj9hwMDXcs8NOqesPSiqr6RlW9Dm78i8A5SS5opvstP8CoMklekOTiJBcl+Z9JHg3cE3hrkguTbE5ydpJ7NuUf0RzjoiQfHbcTSU5q6rkkyStb67cn+ZPmeOcl+YVm/S81y59P8rIk21t9uSTJzYGXAY9p2vmY5X8BaMptaeZ/P8mXknwEuEurzC8l+WDzF5hzkhw6bp96wOyZvXkxe2ZvXsye2ZsXs9f37FWV0woT8N+A14zYvhewZzN/Z+D8Zn4LcMkqZY4DPgPs1Szfqvl6NnDPVh1nM/iBOQj4FnBIu/yy9hwDvH/ZutsB32z23x34GHB8s62ARzXzpwJ/0My/Hzipmf9dYPsK/ToZOL1Vz0uB57WWL2nKHwlc3Hwf9gO+ulQO+Chw52b+3sDH5n3ON8pk9sye2TN7fZvMntkze2ZvXtPuaCxJzgAewOCvB0cBNwNOT3IEcD3wKyvsNqzMQ4E3VtU1AFX1w1Wqvw+Dy7eXjVl+yVHA2VX1g6YPbwUeyOBy8E8Z/CAAfAF4WDN/X+D4Zv5twKvHrGslRwPvXupnkvc2X/cB7ge8M8lS2T3WUM8uzexNxeytA7M3FbO3DszeVMzeOjB7U1no7DkYGu5S4MSlhao6JcmBwPnNqucA/wEczuB2wx+vcIxhZcJgpD6uScu39xvmZ9UM0xn84K4lC9ex8y2Xe7bmV2r3bsAVVXXEGurclZm98Zm99WX2xmf21pfZG5/ZW19mb3y7ZPZ8Zmi4jwF7Jnlma91erfn9ge9V1Q3AE4BNKxxjWJkPA09OshdAkls1668G9l3hOOcCD0pyyLLyq/lss9+BSTYBJwGfWGWf87jpH4WtQ8osb+fXgXs0bbsHcEiz/pPACc39sPsCjwKoqquAy5L8drNPkhw+Zp/6wOyZvXkxe2ZvXsye2ZsXs9fz7DkYGqIZRR/PIFyXJfkcsA14YVPk9cCTkpzH4HLojhUOs2KZqvog8F7g/CQXAksPo70JeEOaB+pabfkB8HTgrCQXAe8Y0uyHJPn20sTgPs4XAx8HLgIuqKp/WKXrzwae2/T3tsCVK5T5OPCrTTsfA5wJ3KrpyzOBLzftvqBp64VNmXNax3gc8JSmP5cCv7VKu3rD7Jm9eTF7Zm9ezJ7ZmxezZ/Zy05UzCZq/XlxbVZVkK4OH6zZccLXrMXuaF7OneTF7mhezdxOfGdJyRzJ4CDDAFcCT59sc9YjZ07yYPc2L2dO8mL2GV4YkSZIk9ZLPDEmSJEnqJQdDkiRJknrJwZAkSZKkXnIwJEmSJKmXHAxJkiRJ6iUHQ5IkSZJ66f8DLHV0HnNRUJAAAAAASUVORK5CYII=\n",
"text/plain": [
"