{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-CSJVBCBTdi6"
   },
   "source": [
    "# Data reduction and 3D analysis\n",
    "\n",
    "\n",
    "In this tutorial, we present the basics steps for a 3D stacked analysis. The main aim is to perform a spectral and morphological analysis of a given source. We are going to analyse the data from the High Energy Stereoscopic System (H.E.S.S.) towards the SNR RX J1713.7-3946. \n",
    "\n",
    "Let's start with some basic imports:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "BoJtbHb8XNnQ"
   },
   "outputs": [],
   "source": [
    "from astropy.coordinates import SkyCoord\n",
    "import astropy.units as u\n",
    "from regions import CircleSkyRegion\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from gammapy.data import DataStore\n",
    "from gammapy.datasets import MapDataset, Datasets\n",
    "from gammapy.estimators import ExcessMapEstimator\n",
    "from gammapy.maps import MapAxis, WcsGeom, RegionGeom\n",
    "from gammapy.makers import (\n",
    "    MapDatasetMaker,\n",
    "    SafeMaskMaker,\n",
    "    FoVBackgroundMaker,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "xDA8PC4LYIjj"
   },
   "outputs": [],
   "source": [
    "data_store = DataStore.from_dir(f\"$GAMMAPY_DATA/hess-dl3-dr1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "uHQnJj-bZnEk"
   },
   "outputs": [],
   "source": [
    "target_position = SkyCoord.from_name(\"RX J1713.7-3946\").galactic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Rjm0nxG1HLuX"
   },
   "source": [
    "We can now define an observation filter to select only the relevant observations. Here we use a cone search which we define with a python dict.\n",
    "\n",
    "We then filter the ObservationTable with `~gammapy.data.ObservationTable.select_observations()`.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "wN4hVfPBG6VR"
   },
   "outputs": [],
   "source": [
    "selection = dict(\n",
    "    type=\"sky_circle\",\n",
    "    frame='galactic',\n",
    "    lon=target_position.l,\n",
    "    lat=target_position.b,\n",
    "    radius=\"5 deg\",\n",
    ")\n",
    "selected_obs_table = data_store.obs_table.select_observations(selection)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "O7Rz1406HPkt"
   },
   "outputs": [],
   "source": [
    "observations = data_store.get_observations(selected_obs_table[\"OBS_ID\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UyedN2taQOXZ"
   },
   "source": [
    "# Preparing reduced datasets geometry\n",
    "\n",
    "We define a reference geometry for our analysis.\n",
    "\n",
    "We choose a WCS based geometry. We define the center of the geometry by our target source, assign a binsize of 0.02 deg, the width of the geometry, and also define an energy axis.\n",
    "\n",
    "This geometry will be what our final dataset will look like."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "TNnZitc1QLDi"
   },
   "outputs": [],
   "source": [
    "energy_axis = MapAxis.from_energy_bounds(0.3, 10.0, 15, unit=\"TeV\")\n",
    "\n",
    "geom = WcsGeom.create(\n",
    "    skydir=target_position,\n",
    "    binsz=0.02,\n",
    "    width=(6, 6),\n",
    "    frame=\"galactic\",\n",
    "    axes=[energy_axis],\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "U-5SNAeQgX0N"
   },
   "source": [
    "Important: we have both a 'true' energy axis and 'reconstructed' energy axis. We must ensure that the true axis is both wider in energy and contains more bins than the 'reco'.\n",
    "\n",
    "![Screenshot from 2024-09-03 15-05-35.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAjgAAADdCAYAAAC2YwjuAAAABHNCSVQICAgIfAhkiAAAABl0RVh0U29mdHdhcmUAZ25vbWUtc2NyZWVuc2hvdO8Dvz4AAAApdEVYdENyZWF0aW9uIFRpbWUAbWFyLiAwMyBzZXB0LiAyMDI0IDE1OjA1OjM14kpfGAAAIABJREFUeJzt3Xl8FPX9P/DX7ibZ3BdXuMJ9XwIKcgalKoIWr59W8UCl0ipe9Sj9VkTsodVSUAsWtGJbRUVFKqhVVAg3EgKEM9xXAkkgN7l2s5/fHx92JyG7IZPs7GxmXs/HYx+Z7Fyf92c2O+98ZubzsQghBIiIiIgMxKp3AYiIiIj8jQkOERERGQ4THCIiIjKcEL0LQMaxp3t3VB45Uud9i90OW1wc7D16IHrkSLS47z5EDBigQwmJfNvZsiWqz59Hwl13oevHH+tdHCJqIrbgkOZEZSWcubm4sHEjcl5/HfsGDcKJX/0KwunUu2jNliMrC+mRkcgcN07vojRZweefY7vViuyXXtK7KERkIGzBIb+zRkejxzffeH4XlZVwnDmDknXrcP6DDyDKy3Fu0SLA6USnd9/VsaTNV/nevRDl5XoXwy/Ktm8H+DAnEfkZExzyO4vNhujRo+u8n3jvvUh69lkcSEmB8+xZnHvvPbR5+mmE9+unQymbt8qjR/Uugt8YKRYiCh68REUBZe/ZE+1mz5a/CIGiGi091HBGSgqMFAsRBQ8mOBRw0SNHeqarTp/WsSTNl5GSgqpjx/QuAhEZEBMcCjjhcHimrZGRXpcp27YNxx9+GLu7dkV6RAR2JiQgMyUF5955p9b63jiysnD6+eexd8AA7IiLQ3pkJPb07ImTjz9eb2JQtHIljtxxBzI6dEC63Y6diYk4cPXVOPvKK6guLva6zoUtW7DdYkFG+/YQDgdEZSVy58/H/uHDsTM+HjuiorBvyBDkLVgAVFd737HLhfNLliAzJQU74+ORbrcjIzkZx6ZMQenGjZ7FnDk52H/VVdjVqhUKP/8cAFCamortFovyClGuOrvLtqttW1k2hwNn//IX7OnZE+l2O3a1aQPhcMBx+rRn/bOvvuqzfrJfesmznPPcOe+hlJcjZ948HBg9GjtbtpSxdOiAo3ffjdLUVM9yxd9+K49PTIxnW2fmzKkVy8Gf/czr9nPnz5fbb9EC6eHh2N21K05Mn46KzEyfZQcA4XQib+FCHBgxotaxyX3rLb/c8O7IzkbWzJmeuHZER2PfoEHInj0b1fn5Xtc5NGECtlut8vMBoCw9HcenTkVGcrI8RklJOHbvvZeNrTH1cvbVV7HdYsGRO+4AADjz8nBi+nT5+Q8Lw+Gf/7zW8kVffYXDkydjV9u2SA8LQ0aHDjg0YQIKPvkE5Xv2eF4VmZlwlZZiR1wctlssOHTDDZetu9PPPYftFgt2tmwJUVFx2eWJGoL34FDAFX//vWc66uqra890uXDq2WeRO3++58ZTW3w8XBUVKF23DqXr1uHce++h+8qVCGnZss62C5Ytw/GHHoLrwoVa71ceOoS8Q4dwbtEidHzzTbT61a+UXZaU4Mhdd6H4kstl1VVVuLB1Ky5s3YqcefPQ7fPPET1mjNeYHNnZKP7+e5x6+mlUZmbCGh0NW3w8HFlZKN+xAydnzEDppk3o8uGHtdYTDgeO3Horir76yvOexW6H49Qp5C9divylS9Hm+efR4S9/gSU8HK7SUlijomApKoJwOGCx2xGalKSsG1L3T9p59iycOTk49dxzKKjx+LOrvByW0FCv8TRGxd69ODRpEqpOnKj1viMrCwUff4yCjz9Gy4cfRvLbb8Nit8N14QJsCQlwlZYCAGxxcbDFx3vWqxkXAJRnZODwz3/u2b4lLAzWqChUHT+Oc4sX4/ySJei0aBFaPPhg3To4fx6HJkxAWVpa7W3u2IFTO3agcPnyyybO9cn/6COcmDYNrrIyAIA1KgqwWlGekYHyjAzkLV6MHitXIvLKK+uuLARKt2yBIycHZ/70Jxl727awhIbCmZOD/A8/ROGXX6LXmjWIHDq0zupNqRcAKN+1C9X5+chMSUHF/v11F3C5cHzaNJxfskT+brEgpFUrOHNzUfzttyj+9ttai4clJ2PAiRNo+eCDyH3jDRR/9x0ubN6MqBEjvO6/urAQuRcTvJYPPwxLeLjX5YhUE0R+srtbN5EGiB1xcT6XqcjMFDsSE0UaIDK6dBEuh6PW/FPPPSfSAJEGiJNPPSWqsrLkDKdTFCxfLna2bi3SAJE5fnydbRetXi3SrFaRBohdHTuKc++/L6pOnxbVxcXiwk8/ieMPPyzSbDaR/cc/KitVV4vM8eM9+zw+bZq48NNPwnHunKg8elTkvv22Z5/bIyJE+Z49tfZZunmzZ920kBCxu3t3UbhqlRAulxBCCEdenjh0882eZQpWrKi1/tnXXvPMy375ZeE8f14IIURVVpbImTdP7O7WTeT98591Yt07aJBIA8SBlBSfdV2zbKeef17GEBYm9l5xhcgcP16c+NWv5L5OnfIsd+aVV3xuL2v2bM9yjry8WvOqzpwRO1u1EmmASI+MFNl/+IOoyMwU1cXFouLQIZH98stiR1ycOHj99Z66EUIIZ0GBZ5tZs2f73HflyZOez83u7t1F0TffeD47VdnZ4uiUKXI7NpsoWbu29srV1WL/qFGe/Rx76CFxIT1dOIuKRMWRI+Ls66+L9Ohoz/wjd93lsxzeFK5aJdIsFpEGiIM33CDKdu70zLuQni72DR0q0gCxMylJOHJza6178IYbPJ+dNItFnHziCVGVnS1nulyi4PPPPWXbe8UVfq2XM6+8IsvVsqU49tBDnr/d/SNGiMzx40Xu228LIYTIW7RIOUazZonqkhK5AadT5C9bJtJjYjx/z2defVWUbtwohBCi4vBhz9/j0Xvv9Vl/OW+8IbdvtYrKY8dU1T1RfZjgkN/4SnCqy8pEWUaGyJ4zR6THxnpOgiUbNtRariwjw3OiOPnEE173UbZ7t0iz2UQaIAq//lqZ4XSKjC5dZHLTtq2SGF3i0hPMuffeu+zJvfzAAZEeFSXSALF/1Kha82omEftHjRLVxcV11ncWFood8fEiDRCHJ0+uNc998jt8yy1e9+2L2gQnLSREHBg7VlSdOlVnOX8kOEfvu89zkir+4Qev61cXFwtXeXmt9xqa4By+/Xb5uYmJEZXHj3td5sidd4o0QOwbNqzW++f+9S8l0XvuOa/rlm7eLLaHhqpOcFyVlWJXu3Zyv0OG1EnYhRDCVV4uMjp18iSaNXkSHJtN5H/yidd9ZL/0kqf8ZTt21JrXlHpxJzhpVqtIs1rFiRkzRHVZWZ31948YIdIAcWjSJK/bz5k3z1MG4XTWmndo0iSZWNvtdf723Pb07Su3//Ofe51P1Fi8B4f8rrqoqNa9FDsiI7Fv4EBkz54NV3Exwvv1Q6+1axE9alSt9XIXLACEgDUyEu3//Gev247o3x8JF+8ZKFi2zPN+0TffeG5WbTt7NkLbtfO6fkirVrX3+fe/AwDCOndG0nPPeV0nvFcvtH78cQDAhY0bUbZjh9flWk2bBmtMTJ33bXFxnub5C9u21ZpXXVQEQF6W0lJIQgK6LV+O0A4d/L7t6vx8FHzyCQAg4c47EXPttV6Xs8bENOrygyM7G4VffAEAaPOb3yCsUyevy7WbMwcAUPbTT7Uuk+X94x8AAFtiItq//LLXdaOuvhrW2FjVZStYvhyO7GwAQIe5c71eIrSEhyNp5ky5fI3PbE1hHTog4c47vc6LvfFGz3TNz09T68XD5ULMNdcg+c03YY2IqDPbHV/UVVd53b77spurpASOvLxa81o/8QQA2RfWuX/+s866pevXo2LfPrnsjBlet0/UWExwKGDC+/ZFj//9D/1270akly/Lkov35kSnpMh7GHyIGjYMgLwh0634hx880wm3396g8lTn53uSlbibbgJsNp/LJtx2m1LOGvtqqLCLiYXz/Pla77tPDgXLluHU00/DmZuretsN0XrGDIS0aKHJtkvWr4eoqgIAT/Lp1+3/+CPgcgEA4iZO9LlceO/engTT/dlwXbiACz/9JNe98Ua/39/h/sxaY2K89v3k5v7MVh0/juqCAlX7CKuRlNb8/DSlXi7V7sUXAYvF+/6TkwHA583K7vct4eEISUysNS/2+usR3qcPACBv0SJPed3cyae9Vy/EermpnKgpeJMx+Z0tLg4DLz7+XV1UhGMPPICSH35Axb59qMrK8v5FKgQqjx8HABR/8w22+/iyrclZ40RRefAgACCkTRuvNx97U3nkiOdG5ojevetd1v0lDQAVhw83aPu1WOX/Epc+rdPxr39FyZo1cObkIHf+fOQtWICY8eORcMcdSLzzTq8tQo3hPklpwV33ABChQaeNNZ98OzB8eIPWcX82Ko8e9Ty9FtG3r2Zlc5WUIL2BN2w7CwpgS0ho+E6syv+hNT8/TamXS9X3+Wj54IMoXb8e+Z98gtgJE9Divvs8f8NlO3Yg+8UXAQAtpkyBJSyszvqtZ8zAycceQ9Xx4yj6+mv5zwRkslZw8WnA1o895jPBImosJjikCWt0tOdn148/xt4BA+A8exYnZ8xA9IgRtRIG4OKj4+4T0cCBiBw8+LL7qHkydT/GbVORENR89NtdXl+s0dHyC1gIz2Ulfwht3x79du1C1qxZyF+6FK4LF1D8v/+h+H//w6mnnkLbF16Ql86swdvYWrMe1dR/Q7lqPDbc4v77L38itNk8lz9rlS0uTrOyhbRujbgal5J8sUZHI7RtW7/uG1BfL2q0ePBBXEhPR97f/47jDzyArN/9DvZu3eDIzZXJrRCIGj4cHebO9b7+Aw8g6//+D9VFRchdsMCT4JxfsgSishLW6Gi0eOAB1eUiuhwmOKS5kJYt0fm993B44kSI8nIcvece9Nm6tdZ/e+5HW10XLiBm3Dh0fOMNVftwt3RUX3zkuCFsNe65uPSx8ku5SkuVx9b9fKIMadMGnRYvRse//Q2Fq1ah4LPPUPTll3CVliJr5kxUFxX5vCcpGNRsZaouLYX/HjyXal726Dh/vqrWD2uNS1JNeQzcF3fZQpOS0Pn99/2+/YbsG1BfL2q1efJJFK5YAVRXwxoZiQs//QRreDiiR49G4t13o+W0aT67HLBGRaHFQw8hd948FH/7LSqPHoW9SxfkLV4MQCZntkbc/0R0OcH7byEZStyNN6LVo48CAMp37kTW//1fnWXCL14mKt+zR/X27V26AJB9vvjqVO1SYV27ev7rrThwoN5la95/YO/aVXX5GsIaHY3EX/wC3T77DP0zMxHevz8AIGfu3MsmYE1S896jS+6RaAh33QPw3DDqT+E1Lh+q/WzU7Eun6swZv5XJzV22ioMHPfchBUpT6kUNV0kJDo4fD+f58+i9eTP6HzqEIRUVuKKwEL3WrUOrX//6sv0ptZ4xQ7ZCCoG8t99GyZo1qDx0SM577DHNyk7mxgSHAqbDX/8Ke69eAICcv/2tzs26sddfD0A+WVF18qSqbcekpHim3df1LyekRQtEDBwIAChcubLek3vhihVKOcePV1W2xgjr0gXtZs0CAIiqKtU3pqoRkpDgSfSqLj4xo0b0mDGeS2gFn33m17IBQPTYsZ6bg89/8IGqdUPbt0dI69YAgAtbt9a/cCNGNHd/ZkVFRYM/d/7SlHpRo2D5clSdPInw3r19Pql1OfauXRE3aRIAWdbz//oXACDm2msRrsG9UUQAExwKIGtEBLp+8IH8b08IHHvggVqtLe7/BIXDgRPTpvn8j7i6oEAmHDVOSPGTJyPk4n/r2XPmwJmT43Xdkh9/hKPG+FfuR1Orjh9Hjo97CCqPHEHOxUtmUcOHe+1NtrFERYXPON03kVqjohDSpk2tee4m/UufymoMS3i4556o4u+/95roOXNyfD7iHJqUhPjJkwEA+R9/jNJ167wuV3HgQJ0kwxoV5WlB8hWLLTbWc4/GuX/+s96n2ErXr/e0DLjF33qrZ175rl1e1zv76qsNbvmrKfa66zx1d/rZZz2PVF9KOJ0o/PJLv96/1dR6aSh378zulteKfftQXVhY72fXG/cj486zZ5H/0UfyPT4aThpigkMBFXnllWh7cTRxR1YWjj/yiGdeWMeOaH+xq/ri1atxYMwYFH75JRxnzqC6oEB+wb7wAnZ3744jt96KvIULPeta7HYkX0xCHFlZ2D9iBAo++QTOc+fgKitDWVoajk2ZgoPjx+PQxIme5Kjlgw8i+mLrz+nf/hannnhCdl1fXAxHVhbOvfceMseMgaukBBa7HZ0WLfJbXQinEwdGjcLuzp1x7p13PCcSuFwo/OILnLl4303ilCl1LgGE9+wJAKjYswe5b76Jin37cPYvf8Hp3/62UWVpebEb/8rMTJx8/HFPWYTTiYJly7BvyBBU1Hha6lIdXn9d3ovjcuHQpEnI+etfUXXyJERFBSqPHEH2rFnYN3gwMq+5plbrnCU0FPbOnQEA+R9+iKKvv0bZjh04MX16rYSq/Z/+hNCOHYHqahyaNAnZs2ejfPduVBcUwHH6NAo++wwHr7sOmWPHIvPaa2uNZ5T0zDOynyGXC4dvuaVWAlZ17BiOP/QQsl54oXFPrFks6Pzuu0BICBzZ2dg/bBjOLV6MqmPHUF1UhMrDh5G3YAH2DRiAI5Mn48S0aer3UY+m1EtDxd98sxxCQwicfeUV7O3XDzsTEpAeEYF0u12O9darlxym5MgRn9uJ/dnPPK01wuFAWHIy4i8Z74rIr/TuaZCMoyFDNQghhHA6xf6RIz29s+a9+26t2dlz5ni6ePf1OnjjjcJZWFhn0znz53t6Ovb2So+NFYUrV9YuTkGByPzZz+rd347ERK899NbsLfjckiU+Qz4+fbqnx1o3x7lzYt+wYZ71t4eFiYzkZLE9PPyyvSOX7dzpvY6sVlF59Kiqsgkhe+Q9MGaMsh2bTWQkJ4v0yEjZC+6VV9bqEfjSnoyFEKIkNdUzbIDXl83mtbfinLfe8rp8RteutYZ1qDh0SOzp3bv+4xQXJ/I//bTOPs4tWeLpJdvd6+6OuDjPdP5nn3l63VU7VIMQcrgGdy/dvl67e/SoM9SHuyfjjE6dfG676swZzzay//CHOvMbWy+enoyBeodIcFVViVPPPCPrKjpabLfb6/37qjh40Oe2ct9+W+k1+89/9rkckT+wBYcCz2ZD1w8/hL17dwDAqSefRGWNvmXavvgi+u7ciZa//CXs3bvDEhEBa1QU7L16oeXDD6PnmjXo8fXXXp9mav3kk+i7cydaPPggwpKTYQkLgy0+HpFXXYV2c+ZgwNGjnsdUPcWJj0fP1avR7fPPEX/LLQht1672ei+/jAGHD/vsobexQlq0QJ/Nm9Fl6VLETZoEW0ICHGfOyKdTxo5F8j/+gV5r13ptWYgYNAg9v/tO9sAbEwNbQgLibr4ZvdasQViNm34byhIWhh7ffot2c+YgvHdvWEJCUF1SgsihQ9Fp8WL02bIFsddfX28rR/TYseifmYmkmTMR3rev57iF9+6NVo89hn4ZGWj30kt11mv92GPo+NZbcp3wcIR26ICWv/wlen3/fa1Hn+3du6Pvrl1IXrAAMddcA1uLFrCEhiKkTRvEXHstOrz2GgYcP+61s8EWU6ei1/r1iLvpJthatIBwOBDSsiVaP/kk+u3di4Tbb0f0yJGq680tbtIkDDh8GEm//z0iBg+GNTYWFrtdtlLcfju6/Oc/6LdvH8I16CeoKfVyOa6yMmSmpCBn7ly0ffFFXHH+PIZUVGCoEBjqcmFwcTH6pKWhzW9+A1gscBUXezrv88b9N2ux29Hyl79sdMxEDWERohF31hERkeHlzJuH07/5DaJGjEDvTZvqXfbAqFG4sGkTWj78MDq9+67XZfZfeSXKtm9Hi/vvR+eLNxoTaYUtOERE5FXVxd7FPfeH+VBdWOjp0drXkBUla9eibPt2ALy5mAKDCQ4REXkVO2ECAKB81y4cmzJFjutVo9HfVVKCwi++wIFRo+A8dw6RV16JxClTvG7L/ZRi1PDhXseiI/I3XqIiIiKfzr7yCrJmzfIMpWKx22FLSICrtFT28H1R7HXXoctHH3kd1LXiwAHs7dsXEAJd/vMfJN57b8DKT+bFBIeIiOpVkZmJc+++i5K1a1F55Aiqi4thtdsR2r49ooYPR4v77vN0eujNiUcewbl33kHEoEHo89NPXgflJPI3JjhERERkOLwHh4iIiAyHCQ4REREZDhMcIiIiMhwmOERERGQ4THCIiIjIcEyV4JSXl+tdBCIiooDZuHEjvvrqKxQWFupdlIAz/GPiLpcLO3fuxPr16zFu3DgMGjRI7yIREREFzKeffoovvvgCffr0wZgxYzB8+HBEREToXSzNGTbBOXnyJFJTU7Fx40aUlJTg8ccfx9VXX613sYiIiALuv//9Lz755BMAgNVqxZVXXolRo0Zh8ODBCAkJ0bl02jBUglNaWooNGzZg/fr1OHbsmOf9xx57DKNGjdKxZERERPpatWoVli5dCgCwWCwQQsBut2PkyJEYPXo0evfuDYvFonMp/afZJzjV1dVIT0/H+vXrsX37dgghYLPZ4HQ6YbFYMH36dKSkpOhdTNKBw+EAAAghcOnH3P27t/f9sWx962q970CU0x/Lepvvazk1+21qnTW2HrXYr7+2GUxl4THW9xifP38eubm5EEJ4khx3UhMXF4cxY8Zg1KhR6NSpE5q7ZtsudezYMaxbtw7r1q1DeXk5rFar58C6k5s2bdogNTUVa9eurbWu1l/m/t6ekcrpTjqIyNxCQ0M90+4T7KWtBxaLxet7ly7b0Peaup4WZVFbvqbu153U1Pw9JCQEDocDVVVVcDgcqL44sGpz12wTnJCQEISGhsJut6O8vBwulwtA7YPnzkxrvmpq7B+Ar/f9tWxzKafaZVlO85bT/Xt9619uWzzRaVuWmgkHGdOGDRuwcOFCAIDNZoPL5YLFYsHAgQORkpKCIUOGwGaz6VxK/2n2l6iEEDhw4AA2bNiATZs2obKy0vO+xWLBPffcg5tuuknnUhIREelny5YtePPNNz2/d+rUCSkpKRg9ejSio6N1LJl2mn2CU5PT6cSuXbuwYcMGbNu2zdOqc++992LixIk6l46IiCjwfvrpJ7zxxhuIiYnBmDFjMGbMGCQnJ+tdLM0ZKsGpqby8HFu3bsX69euxf/9+3HLLLbjzzjv1LhYREVHA7N27Fz/88ANGjx6NQYMGGeoS1OUYNsGpqaCgABs3bkSLFi0wYsQIvYtDREQUEOXl5abo1M8bUyQ4REREZC6mGouKiIiIzIEJDhERERkOExwiIiIyHCY4REREZDhMcIiIiMhwmOAQERGR4TDBISIiIsNhgkPBacEC4Jpr5GvjRr1LQ02xcaNyLBcs0Ls0RGQSqkcTX71aOd/8v/8H9Ovn7yKp8OqrQEUFEB4OzJypWzF27AD++185fc01QEqKbkXR1JIlwIkTcvrpp4G4OA13dugQsHatnM7L03BHtRUWAvPny+nOnYGpUwO264BKTQXWrJHTt9wCXHGFhjvLy1OO5aBBGu6orjfeAAoK5PTs2cAlA28bwsmTwHvvyemBA4HbbtO3PFr58ksgPV1O338/0LWrvuUJBps2Ad99J6cnTgSGDfPjxnfuBFaskNPjxsmXTt5+G8jJkdO/+x1gtzdsvUYlOK+/Lqf79AmCBKeoSJ5pdU5w5syR0xaLcROc998H1q2T01Onapzg6KSwUDmWY8caO8Fxx5mcrHGCo6M33wSOHpXTs2YBRhyG5+RJ5Vjec49GCU5FBbBli5yOj9flA7NyJfDuu3J65EgmOACwebNy7BMTNUhw3BsHdE9wdu+W008/rWGCU/M/IA7yIJmlTswQpxliBBinkQQkxrNnZfM0IP+Dc7fIBZAZjqVaZqmTxsap+h4ca401jFyhapilTswQZ80YXS79yqG1ml8YjLN5M8tJzgzfP2rx2F9mPbU7MkuFqmGWOjFDnGaIEWCcRmKGGAHzxKmGWeokYC04ZqlQNcxSJ2aI0wwxAozTSMwQI2CeONUwS53wEpWOzFInZoiTl6iMxQxxmuUkZ4bvH7V47C+zntodmaVC1TBLnZghTjPECDBOIzFDjIB54lTDLHXCFhwdmaVOzBAnW3A0VlERwJ2xBUcTAT6Gbmb4/lHLLAkOW3B0ZJY6MUOcZogRCHCc4eHKtI4JjlGPZ0Bi1PEYupnhWKplljrhTcY6Mkud6BZnAL9QeSw1UPPkGGBmOJ4BT3B0YoZjqZZZ6sTcl6h0+o/CLSjrRAMBjVOn/xh5icpYzBCnWU5yZvmeVSNgx76ZnmNV92T8wgvyFRTCw+VQDZWVuhbjjjvM8Qe3cmUAd6bTf4yJieY4ls8+K19G5+7e3ciGDDHHZ3buXPkixa9/LV+aCILLkm4bNjRuveY9mngQNJsSEREZjgHOr807wSEiIiLyggkOERERGQ4THCIiIjIcJjhERERkOExwiIiIyHCY4BAREZHhqO4H5/Rp4PBhOd29O9Chg7+L1Pzk5wMZGXK6XTugZ099y6OV/fuBnBw5PWQIEBurb3m0IASQmiqnIyOBYcP0LY9WzpwBMjPldJcuQKdOAdpxgPvT2LkTKCyU0yNGAHZ7QHcfEBUVwJYtcjohARg0KAA71MHx4/IFAL17A0lJuhQjqOTmAvv2yenkZKBrV33Lo5WMDHmeBeR3cmRkw9ZT3YLz+efANdfI17Jlatc2pm3blDoxckdUf/iDEqdRO1BzOpUY77lH79Jo56uvlDj//e8A7jjAJ8fHH1fidCfnRnPmjBLjU09ptJMg6PTt/feVOP/3P12KEHR+/FGpk8WL9S6Ndn77WyXOEycavh7HovIDs9QJx6IyDo5FZRwci8q8zFInHGxTx66kg7ZO/MwMI1DzWGqACY6mzBAjYJ441TDLWFS6DLap++B1QdBsCphnELiAHnudTopmGJwRMGecRv3bNOP3j5HjVEPTOgmS8yvQ+DibdwtOEDSbAkFWJxoyQ5xmiBFgnEZihhgB88SphqZ1EiTnV4CXqHRlljoxQ5xmiBFgnEZihhgB88Sphlk4NaUuAAAZ6ElEQVTqxJyXqIKEWerELHG6/5iMHKNZjqUZ4jRDjIB54lTDLHXS2DjZguMHZqkTs8VphhgBxtncmSFGwDxxqmGWOmELjo7McvObWY69GVpweJOxcZjx+8fIcaphljox503GQcIsdWK2OM0QI8A4mzszxAiYJ041zFInvMlYR2apE7PFaYYYAcbZ3JkhRsA8caphljppbJyqx6KaOBFYs0ZOd+6sdm1j6t9fqZM2bfQti5aeeQa4+245PXCgvmXR0vffy8s2YWF6l0Q748crn9nk5ADuOMD9abz2GjBzppwO2HhbARYfrxzLgIwPp1OfKHfdBQwdKqeNOt6fWiNGKMfeyONCvvQS8OijcrpXr4avpzrBSUriIGeXio0Fxo3TuxTa69nTHF8sY8fqXQLttW4tXwEX4JOjkRNxt5CQAH3/2O1AZaVuCU6HDsY+iTdGYqI5zj19+8qXWqovUREFnM69aFITBVGHYdQEPI7UzBgnweFJ0FiCqJtwaiKeGImat2b6Hdy8ExyeBI2LJ0UiIv0Y4PxqnASHiIiI/MMA59fmneAQERERecEEh4iIiAxH9WPieXnA3r1yumNHoFs3fxep+blwAdi2TU63aAEMGKBvebRy9Chw8qSc7tcPaNVK3/JoZdMmoKoKsNmAMWP0Lo028vOBjAw53a6dcR//37cPyM2V00OGBKifmAATAkhNldORkcCwYfqWRytnzgCZmXK6a9cA998UpIqKgB075HRSEtC7t77l0Upmpjz+ADBoEJCQ0LD1VLfgpKYC11wjXwsXql3bmI4dU+rE3amYEb39thKnu3MpI7rlFhnj9dfrXRLtbN2qHMu//U3v0mjnD39Q4nT/Y2Y0DocS47336l0a7axapcT573/rXZrgkJGh1Mmf/qR3abTz+utKnGlpDV+PQzX4gVnqxGxxmiFGgHE2d2aIETBPnGqYpU44FpWOzFInZovTDDECjLO5M0OMgHniVMMsdaJLguNyqV3bmGoO5W7kOjHLsXfHaYYYgQB/MQa4P42af5tGPQEEPEad+kQxy/esGmb4fAONP/ZswfEDs9SJWU6KbMHRkE7HEjDu8QxYjO5+USorNdyJb2Y4lmqZpU54iUpHZqkTs50UeSz9RMcOw8zwtxnwBEcnZjiWapmlTgKW4ARtM6GOXUkHbZ34mdniNEOMQADi1PHEaJbPrBkuq5rlWKoRsDrReagGc16iCpKxMoKqTjRkxv/6jYqfWWNhq6M5aVonQXJ+Bcx6k3GQjJVhlv8sAnrsgyTBMerxNMvJwmw3YZohRsC4f5dqafr5DpLzK2DWFpwgYZY6YZzGYYYYAfPFaYYYAWPHqYZZ6oQ3GevILHXCOI3DDDEC5ovTDDECxo5TDbPUSWPjVD0W1fDhSjf97dqpXduYkpKUOomP17csWpo6FRg3Tk736qVnSbT18cdAebmcttn0LYtWBg1SPrNJSfqWRUvPPgvcc4+cHjRI37JoafVq2XRvt+tdEu2MH698ZjkOldSjh1InRh0bEABmzAAmT5bT/fs3fD3VCU5ionKSI8luN0eddO4sX0Y3fLjeJdBeXJw5PrM9exp3INGaxo7VuwTaa9XKHJ9ZNaKizFEn3bo1bmBv1ZeoiIiIiIIdExwiIiIyHCY4RER0eUHULwpRQzDBoeDHL1Pj4LFsvpjgUDPDBIeCH79MmzeeGIlIB8ZJcPjFSRScgqhHVCJqhGZ6fm3eCQ7/MzQunhSJiPRjgPOr6gRn82bZq6DFAtx9txZFUiFIToIFBUqdDB6sd2m0M3++Euf8+RrvTMdjO3iwEmd+vm7F0FRGhhLjzTfrXRrtPPOMEudnn+ldGu20bCljjIrSuyTa+fpr5Vg++qjepQkOR44odXLttX7eeJCcXwFg9mwlzn/9q+HrNe/BNoOEWerEjHEatftzsw1CCZgjTrP8XRo5TjXMMgBpY489x6LyA7PUCeM0DjPECJgvTjPECBg7TjXMUiccbFNHZqkTxmkcZogRMF+cZogRMHacapilTnRJcIzcJKaGWerEjHEa9UvDbJduAHPEaZa/SyPHqQYvUdWPLTh+YJY6YZzGYYYYAfPFaYYYAWPHqYZZ6oQtODoyS52YMU6jfmmYrWUDMEecZvm7NHKcarAFp35swfEDs9QJ4zQOM8QImC9OTWPUuV8UsxxLNcxSJ7zJWEdmqRPd4gzwl6kZjiePpbHUjFMzTHCCjlnqhJeodGSWOtEtTp4U/S7gl27sdvkzwMfSbJeoAON+B5nle1YNMxx3gJeodO1KOmjrxM8Yp3EEPEadekU1w7EEzBGnGWJUK2B1ovNQDY2NM0TtjpKSZLfJANCjh9q1/SxIxsoICVHqJCFBt2Jo7sorlThTUjTemY7dhD/8MHDddXI6Nla3YmgqIUE5ll266FsWLY0bp8Q5aJCuRdHUU08pw4oE5HKVDrp0UY7lwIH6liVYREUpddKunZ83HiTnVwAYPVqJc9iwhq9nEaIZ58JPPQW88Yac/uIL4JZb9C0P+c/OncrAXpMnAytW6Fseapr4eKCoCIiLAwoL9S4NNca4cUBqqpw+dgzo3FnP0lAguLPlTp2A48d1LUpjNO/RxImIiIi8YIJDREREhsMEh4iIiAyHCQ4REREZDhMcIiIiMhwmOERERGQ4THCIiOjygqhfFKKGUJ3g5ObK7hDGjQOefNL/BWquJkyQdXLrrXqXRDupqcqxf/vtAO44wF+mv/udEuexYwHddcCUlCgxPvJIAHcc4GP56adKnMuXB3TXATVtmhJnaalGO9E5wTl6VInxhRcCvvug5HQqdXL//XqXRjurVilxfvhhw9dT3ZNxZaXS15PTqXZt41q/HigrA1q21Lsk2snLU4790KEB3HGAv0z37FHiLCkJ6K4DxuFQYiwqCuCOKysDuDPg9GklTiP3A7p9u+wbE5DH1ohKS5VjadQextVyuZQ60X1kAQ2dOaPE6e5lviE4FpWfuOvFLAOeBdWx9zMzDNAY8Bh1GnqDgxEahxliVCtgdRIk51eAg23qwl0vQVUnfhbQY6/jWFRB+xn3Iw62aSxmiNMMMarFwTbr17wTHB1PgpdiguNnTHA0ZYYYAcZpJGaIUS3N68Ru12Cj6jHB0RkvURkHL1EZBy9RGYcZYlRL8zoJknMsL1HpjC04xmGGOM0QI8A4jcQMMaplljoxZwtOEGELjnGYoXXDDDECbMExEjPEqJZZ6oQtODpjC45xmCFOM8QIME4jMUOMTWHkOmELjs6Y4BiHGeI0Q4wA4zQSM8TYGDz3+MYEx094ico4zHD5xgwxArxEZSRmiLExzHbuUROn6p6M7XYgJUVO9+ypdm3jGjMGKCgAwsL0Lol2WrY0x7Hv10+JMyZG37JoxWpVYkxO1rcsWmrfXomzUyd9y6KloUOV7x7NHnzReaiGyEjlWA4aFPDdB61x44DqaiAxUe+SaKdtW+XYd+3a8PUsQjTj/99eegmYM0dOL1kCTJ2qZ2nIn44fB7p0kdPDhwNbtuhaHGqizp2BEyfkdHl50Dx+SipMnQr8619yes0aeWYlY4uPl2O5xMUBhYV6l0Y1jiZOwY8jFxsLjycRBYBxEhx+aRIFL7bYEDVfzfT8ygSHghNPiMbC40nUfFVW6l2CRmneCQ6/NI2Lx5aISF/N/HuYCQ4RERHV1czPsc07wSEiIiLyggkOERERGY7qBKe6WnZ/MG4ccM89/i9QczVlilIvTqfepdHG3r1KjH/8o96l0c7cuUqcP/2kd2m0c/31MsbbbtO7JNpZu1Y5losW6VwYDc2cqcTp7m7IaIqKlBinT9e7NMFj8mRZJzfeqHdJtLN5s3Ls589v+HqqezIWAkhNldPuftgI2LoVOHJEThu1y+zCQuXYt22r8c5CQ5WOxLp313hntWVmKnHm5wd01wGVmgpUVQGtWwdgZ1ddBSQkyOnQ0ADsUMrJUY7lsGEB223AZWQocZaW6lsWrVRVKTGWlOhblmCyYYP8noqK0rsk2jl/Xjn2ffs2fD3VCU5Qje1SszA6N5sEVb1oJKAxRkXJ3lJ1YIZjCQA2m/wZkBjfeScAO6mr5rGsrtalCAERkDjtdmVah245zPJ3qZa7XjQ57u4vCUCeY0NUpwx+0dhj37yHaiAiosCYPh1YvFhOf/MNMGGCvuUhugwmOEREdHmHD8v/4tu0US43EgUxJjhERERkOHxMnIiIiAyHCQ4REREZDhMcIiIiMhwmOERERGQ4THCIiIjIcJjgEBERUcBt2rQJX331FQoLCzXZPh8TJyIiIl0sW7YMK1asQJ8+fTB27FgMGzYMERERftl2QBMcIQQqKir8VngiIiJq3lasWIFly5YBAKxWK4YNG4YRI0Zg8ODBCGnC8BABSXByc3Oxdu1aFBUV4aGHHoKt5vgWREREZGqrVq3C0qVLa71nt9sxcuRIjB49Gr1794bFYlG1Tc0SnIqKCqxfvx6pqak4evQo+vfvj+eff75J2RgRERmXw+GAECIgLwBwXRy5Ue3PQJXR6OW+VH5+PnJyciCEgMVi8fwEgPj4eIwePRqjRo1Cp06dGvR58mu2IYTAnj17sHbtWmzdutUTRJ8+ffDss88yuSHSgcPhCMovQH/vryHxBUPcZiy3w+FQ96ElvwoNDYXFYgnIC5CXmS73c//+/XXK6f5cuVksFlitVlRXV6OiogIOhwPVKoZN90vGkZubix9//BGpqakoKiryZF4A0KNHDzz//PMICwvzx650Y5YvIpa7eZebtBOoL35vP9Xsx2q16nLCamq59TrRGq3czdWGDRuwcOHCOu9fccUVGDt2LIYOHar69pZGX6IqLy/Hpk2bsHbtWhw5cgQ2mw1Op1PZsMUCu92OTp061VuoYDth8T+N+tX8TwCo+0UAeP+S8LVsQ9f3trze6/vaptr1L/1yasr6zam+va0TDCeO0NBQEFHgbNmyBW+99Zbnn7T27dtj3LhxGDNmDGJjYxu93Sa14ISEhHhaZqqrq2Gx1L5mZrFYEB4ejsjISISEhCA0NNTzstvtnmWM8oXdHE+Q/DInIiK9pKWl4a233kJERATGjRuH0aNHo3Pnzn7Ztl9uMi4sLMSGDRuwYcMGnDx50pOFWSwWtG/fHi+++CJiYmKaXFgiIiIyhj179uD777/H6NGjMXjwYL8/Ye33p6iys7M9yc65c+cAAB07dsQLL7zAJIeIiIgAyFtdtOwXT9N+cA4ePIiNGzdi48aNiI2NxaxZs5CQkKDV7oiIiIgABKijv+rqauzevRuHDh3ChAkT2JJDREREmuJYVERERGQ4HE2ciIiIDIcJDhERERkOExwiIiIyHCY4REREZDhMcIiIiMhwmOAQ+UlhIWCxNO114IDeURARGQMTHCIiIjKcJg22SUSKyEhgyRLv87ZuBf7xDzl9993A9dd7X65tW23KRkRkNkxwiPwkLAyYOtX7vPBwJcEZNsz3ckRE5B+8REVERESGwwSHiIiIDIcJDlGQmDBBPkn11lvy9127gIkTgYQEICICeP11+f4HHyhPXW3Z4nt748bJZfr3971MdjYwcyYwYAAQEwNERwODBgGzZwP5+U2L57PPgMmT5X1FdjuQlATcdhuwerX35bdskeVt3x5wOIDKSmD+fGD4cCA+HoiKAoYMARYsAKqr69/3tm3Aww8DXbvKuktIAFJSgHfekdv2JjxcXmbcsUP+/sUXwFVXyXur4uKADRuUZXNzgd//XtZVdLSsu/79geefl/ves0d5FRQAL78sY4uMBH74of6ynzwJtGwpl3cfcyJqBEFEmvvoIyEA+Zo3z/syN9wg50+bJsS2bUJERyvrAELMnSuX+89/lPc2b/a9z5QUuUy/ft7nL10qRGSksq2oKCFiYpTfk5JkOdTKzxfi2muV7VitQiQmCmGzKe898ogQ1dW119u8WZn/9ddC9Oolp6OjhejQQQiLRZl/zz3e911dLcTTT9deNj5eiPBw5ferrxYiL6/uuna7nP/BB0IsXFi77gEhtm+Xy23bJkSLFsr7kZEyvkuXd7/ee0+InBxl+9dfX3/9zZwpl4uIEOL8efX1T0QSW3CIgkxeHnDvvUBpqWz9SEkBfvYzYMQI/+3jq6+AKVOAsjLghhuAnTvl/oqLgfR0YOhQ4OxZ4OabZXkayukEJk0CfvxRtmz84x9AURFw/rz8+corQGgosHgx8Oc/+97Oz38uW2lWrZJlOnVKtprcfLOcv3Qp8N//1l1v5kxg3jyZWjz1FJCVJVtQSkuB5cuB1q1lS9EvfuF739u3y3UtFqBHD2D8eNm61q0bUFUF3HGHjKdPH7mtCxfk76dOAbfeqmznoYeAd9+Vy7duDdx5p3x/9Wrg8GHv+3Y4lCfx7r4bSEz0XU4iugy9MywiM1DTghMSIls7lizxvlxTW3AqK4Vo107OGzJECIej7rrl5UJ06iSXef75BgR40VtvKWVbvtz7Mu7WkcjI2i0pNVtwRo0Sori47rqFhbJFBhBi8uTa8zIylJabJ57wvu/du5WWpK+/rj3P3cISEiLrZ8uWuut/+61SRm+tW+XlQrRpI+e/9lrteWlpyrq/+Y338n3yibJMerr3ZYioYdiCQxRknE55D4xWj5IvXy7vvQGAuXOBEC+dRYSHy9YQAFi2rOHb/vvf5c+UlNqtGTU98oi8z6asTLYkeTNtmryv5VJxcUpL1rZttectWCBTg8hI361D/fvLFhXAd1xOJ/DRR/Len0u5681mk/cDXSo8XLnn6ciR2vOGDlXKvmQJUF5ed/1Fi+TPkSOBwYO9l4+IGoYJDlGQiYkBnn5au+1//72yn9GjfS83bJj8efy4vMxzOadPA5mZcnriRN/L1UwO0tMvv91Ldeggf54/X/t9d1wpKfKGZF/ccfna9zXXAGPHep+XnCx/Vld7v8zkcinvt29fd/4TT8ifBQUyiarp0CFgzRo5PWOG7/ITUcMwwSEKMomJ8v4VrRw9Kn+WlMj7YXyNizV0qLJOQxIc93YB4Le/rX/MrZUrG77dS1kvfms5ncp7QshEDAC++ab+fT/zTP37dicx3qSkyCezANnKlJWlzKuoAJ59FjhxQj41ds89dde/4w4l8Vm4sPa8RYtkHElJSisTETUeezImMpmKCvmzdWvgxhsvv3x0dMOGkHBvF5BDUTRkndtvv/wyDeFwKI+ODxzYsMs7/fqp34/NJi/xTZwIrF8vbzzu3VteFtu/Xw64GhoqH0fv1q3u+iEhwK9+BcyaJW9m3rpVXgqrrATef18u88gjchtE1DRMcIhMxv1kTlKSclL153YBeRL3dQ+OFsLC5GWpCxdk/z9vvKHdvgYNAu67D/jLX2Ryk5Uln/Tq0EG2vDz1VP3J0/TpwB//KJOahQtlgvPZZ/KSW0iInE9ETcdLVETNjM2mTLtc6tfv3Vv+PHhQPvbsLz17KmXbs8d/220od1xa7/uNN2Ry8+ij8vH6vDyZrBw5IltuLtcy1KqV8pj6smUysXHfXHzbbUC7dtqWn8gsmOAQNTM1W0rcT/Wo4R7JvKIC+Pxz/5QJAGJjlRt4ly5tXPLVFO641q+XvQFrZf58+fOGGxq/DffNxhUVwGuvyTIDvLmYyJ+Y4BA1M4MHy5tlAeC777wvs2mTvMfDm+uuk53UAfKmWF9JktMJfPml7KCvodwn7gMH5PAEvhw+DKxb1/DtNsSvfy3vXXE45A3AvlqnCgqAFSvkDb2NUVYmf86aJesnK0s+8l1R0fCkbsgQYNQoOT1vnvw5cCAwZkzjykREdTHBIWpmWrdWHsN+/33gvfeUE2thIfDqq7L3XW/9rAAyOXr3XXm/R3a2bHVZvBg4dkwmM4cPyz5lBgyQY0lNm9bwst11F3DTTXJ6zhzZW/LGjcC5c/JSzpo1soffPn2Aa6+V4235S8eOwJ/+JKdXr5bJwpdfAmfOyKRm507ghReA7t3l/UGXPsXUUPfdJ39mZMj66dBB3mQcESETrFatZM/H779ffxLlTgbdY2Ox9YbIz/TuaZDIDNT0ZNyp0+W3d/SoHCvKvc3wcCGSk2UvvIAQU6cK8dxz9Y9FtWqVELGxvsdQAoTo0UOIPXvUxVpaKnsZrm+7Npscc8nlUtar2ZOxr16chRBi+nRlG97MmSPHv6pv/zfeKHtFrsndk/EDD9Qf3+7dsqdjq1XWX81xry59/frXvrfjcMgxtgAhEhKEuHCh/v0SkTpswSFqhrp0AdLSZOtK+/byEemqKtl68t13sqfcMWNq35B8qUmTZGvN738vL3vFxsr+W5KT5ePb//kPsG+f+sepo6LkJaCVK5WbZkND5YjeQ4cCzz0nL2G98opyqc2fXnxRttb88peytSYiQpapVy85wviaNcDXX8tekdVatkzGEBUl91FUJFvPhJA3Gp8+DXz6qdwXIG8e9nWJz2ZT+jt68EHZCkRE/mMRorFXoomIzMPlAlq0kJcBv/1WuanZm9WrlfmnTim9L9e0apUcPNRqlU+0ees3h4gajy04REQNUFQkkxtA9rdTn61b5c8uXbwP2QAAf/2r/HnjjUxuiLTABIeIqAESEpTH4B99FPjnP4H8fGW+yyXH4nrmGeCll+Tlt7lzvV+G274dSE2V07y5mEgbvERFRNRAx47J+5z27VPeS0yUT6SdP68MFxEbK59Eu/de79u5+27g44+BHj1kUqTFvUhEZscEh4hIBYdDJieffy4fc8/Jke/FxwN9+8oOAKdNk4/ze3PypLwkZbXKno/vvz+w5ScyCyY4REREZDi8B4eIiIgMhwkOERERGQ4THCIiIjIcJjhERERkOExwiIiIyHCY4BAREZHhMMEhIiIiw/n/7gZALhN6ezoAAAAASUVORK5CYII=)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define the true energy axis based on the above information"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "energy_axis_true = MapAxis.from_energy_bounds(\n",
    "    0.1, 20, 20, unit=\"TeV\", name=\"energy_true\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8V1sauzKQRRg"
   },
   "outputs": [],
   "source": [
    "print(geom)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E14epAwDQhz8"
   },
   "source": [
    "\n",
    "# Data reduction - FoV background method!\n",
    "\n",
    "Create the maker classes to be used\n",
    "\n",
    "The `~gammapy.datasets.MapDatasetMaker` object is initialized as well as the `~gammapy.makers.SafeMaskMaker` that carries here a maximum offset selection.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aFq9Cmv7VaGW"
   },
   "outputs": [],
   "source": [
    "stacked = MapDataset.create(\n",
    "    geom=geom, energy_axis_true=energy_axis_true, name=\"rxj-stacked\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5utUHw_GczVL"
   },
   "source": [
    "We need to mask out any sources with significant emission.\n",
    "\n",
    "You can check the http://gamma-sky.net to see if there are sources in the FoV and create an exclusion region based that."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hSSeoGIJQgI7"
   },
   "outputs": [],
   "source": [
    "circle = CircleSkyRegion(center=target_position, radius=0.7*u.deg)\n",
    "regions = [circle]\n",
    "exclusion_mask = ~geom.region_mask(regions=regions)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "nzeMvsgnNV21"
   },
   "source": [
    "## Try to plot the exclusion mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "7MZA_j9SQn2a"
   },
   "outputs": [],
   "source": [
    "exclusion_mask.plot_interactive()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "0c3hZ3wWgv7r"
   },
   "source": [
    "Create the `makers`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "Eew5DSBsQsDd"
   },
   "outputs": [],
   "source": [
    "maker = MapDatasetMaker(\n",
    "    selection=[\"counts\", \"background\", \"psf\", \"edisp\", \"exposure\"]\n",
    ")\n",
    "safe_mask_maker = SafeMaskMaker(\n",
    "    methods=[\"offset-max\", \"aeff-max\", \"bkg-peak\"], offset_max=\"2.3 deg\"\n",
    ")\n",
    "fov_bkg_maker = FoVBackgroundMaker(method=\"fit\", exclusion_mask=exclusion_mask)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Perform the data reduction loop recalling the order:\n",
    "- maker\n",
    "- safe mask\n",
    "- FoVbkg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "A_zQyQChVtUq"
   },
   "outputs": [],
   "source": [
    "for obs in observations:\n",
    "    dataset = maker.run(stacked, obs)\n",
    "    dataset = safe_mask_maker.run(dataset, obs)\n",
    "    dataset = fov_bkg_maker.run(dataset)\n",
    "    print(\n",
    "        f\"Background norm obs {obs.obs_id}: {dataset.background_model.spectral_model.norm.value:.2f}\"\n",
    "    )\n",
    "    stacked.stack(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6QyI2en2V6-6"
   },
   "outputs": [],
   "source": [
    "print(stacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2prvksOucHz9"
   },
   "source": [
    "## Compute an excess and a significance map"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "FoSTmGvu1url"
   },
   "outputs": [],
   "source": [
    "estimator = ExcessMapEstimator(\n",
    "    0.1 * u.deg, selection_optional=[], energy_edges=[0.4, 20] * u.TeV\n",
    ")\n",
    "maps = estimator.run(stacked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "PHO3nbB1cGIa"
   },
   "outputs": [],
   "source": [
    "ax1 = maps.sqrt_ts.plot(add_cbar=True)\n",
    "ax1.set_title(\"Significance map\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "pTl3Xo76hMD4"
   },
   "source": [
    "### Well done!! You have gone from DL3 data products to a DL4 dataset, which can be used for SCIENCE! 🥳\n",
    "\n",
    "\n",
    "## Lets do some science!\n",
    "\n",
    "# Exercises: 3D modeling - fit a model to our dataset\n",
    "\n",
    "\n",
    "The 3D analysis consists of a simultaneous fit of the spectral and spatial parameters of one or multiple sources. Here, we focus only on a single source, centered at the coordinates l=347.269°, b=-0.257°. To perform the fit, we firstly need to define a model: we can try with a Gaussian morphology and a powerlaw spectral shape."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import GaussianSpatialModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model = GaussianSpatialModel(\n",
    "    lon_0=target_position.l, lat_0=target_position.b, sigma=0.4*u.deg, frame='galactic',\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1 - Show on the map that this is a suitable guess for the emission extension and position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax = maps.sqrt_ts.plot(add_cbar=True)\n",
    "region = ...\n",
    "artist = ...\n",
    "ax.add_artist(artist)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2 - Create a spectral model, use a simple PowerLaw here. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import PowerLawSpectralModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = PowerLawSpectralModel(...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We limit the model position inside a 1 deg box centered on the reference position of RX J1713"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spatial_model.lon_0.min = spatial_model.lon_0.value - 0.5\n",
    "spatial_model.lon_0.max = spatial_model.lon_0.value + 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3 - Limit the `lat_0` values yourself"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import SkyModel, FoVBackgroundModel, Models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sky_model = SkyModel(\n",
    "    spectral_model=spectral_model, spatial_model=spatial_model, name=\"rxj\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we define a global ~gammapy.modeling.model.FoVBackgroundModel in order to finely adjust the level of residual CR backgroud. This should not be forgotten.\n",
    "\n",
    "Note: we assume the background is well know and can be taken from the IRF bkg model\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_model = FoVBackgroundModel(dataset_name=stacked.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we assign these models to our reduced dataset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stacked.models = Models([sky_model, bkg_model])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4 - Look at your stacked dataset. It should now contain information on the models we just applied"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(stacked)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5 - Fit the model\n",
    "\n",
    "The `~gammapy.modeling.Fit` class is orchestrating the fit, connecting the `stats` method of the dataset to the minimizer. By default, it uses `iminuit`.\n",
    "\n",
    "Its constructor takes a list of dataset as argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling import Fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit = Fit()\n",
    "result = fit.run([stacked])\n",
    "print(result)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Inspecting residuals\n",
    "\n",
    "For any fit it is useful to inspect the residual images. We have a few options on the dataset object to handle this. First we can use `.plot_residuals_spatial()` to plot a residual image, summed over all energies:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "region = CircleSkyRegion(center=target_position, radius=0.4 * u.deg)\n",
    "stacked.plot_residuals(kwargs_spectral={\"region\": region});"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6 - From here you can create the other science products\n",
    "i.e. butterfly plot, flux points, test source significance\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "colab": {
   "provenance": []
  },
  "kernelspec": {
   "display_name": "gammapy-2.1",
   "language": "python",
   "name": "gammapy-2.1"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}