simulate_blazar.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation using CTA IRFs\n",
    "\n",
    "- Intrinsic spectrum assumed to be extrapolation of the Fermi-LAT spectrum (showing example of `gammapy.catalog`)\n",
    "- Add EBL absorpsion\n",
    "- Simulated CTA DL4 data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SkyCoord\n",
    "from gammapy.irf import load_cta_irfs\n",
    "from gammapy.maps import WcsGeom, MapAxis\n",
    "from gammapy.modeling.models import (\n",
    "    PowerLawSpectralModel,\n",
    "    PointSpatialModel,\n",
    "    SkyModel,\n",
    "    Models,\n",
    "    FoVBackgroundModel,\n",
    "    EBLAbsorptionNormSpectralModel\n",
    ")\n",
    "from gammapy.makers import MapDatasetMaker, SafeMaskMaker\n",
    "from gammapy.modeling import Fit\n",
    "from gammapy.data import Observation\n",
    "from gammapy.datasets import MapDataset\n",
    "\n",
    "from gammapy.catalog import SourceCatalog4FGL"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use gammapy.catalog to find the Fermi spectrum\n",
    "- see https://docs.gammapy.org/0.19/tutorials/api/catalog.html for details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'LAT 8-year point source catalog'"
      ]
     },
     "execution_count": 96,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "catalog = SourceCatalog4FGL()\n",
    "catalog.description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "## explore some other catalogs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [],
   "source": [
    "#catalog.table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [],
   "source": [
    "source = catalog[\"PKS 2155-304\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "*** Basic info ***\n",
      "\n",
      "Catalog row index (zero-based) : 5400\n",
      "Source name          : 4FGL J2158.8-3013\n",
      "Extended name        :                   \n",
      "Associations     : PKS 2155-304, PKS 2155-304, 3FGL J2158.8-3013, 3FHL J2158.8-3013, J2157-3015, 3EG J2158-3023, EGR J2200-3015\n",
      "ASSOC_PROB_BAY   : 1.000\n",
      "ASSOC_PROB_LR    : 1.000\n",
      "Class1           : bll  \n",
      "Class2           :      \n",
      "TeVCat flag      : P\n",
      "\n",
      "*** Other info ***\n",
      "\n",
      "Significance (100 MeV - 1 TeV)   : 224.639\n",
      "Npred                            : 15833.8\n",
      "\n",
      "Other flags          : 0\n",
      "\n",
      "*** Position info ***\n",
      "\n",
      "RA                   : 329.714 deg\n",
      "DEC                  : -30.225 deg\n",
      "GLON                 : 17.731 deg\n",
      "GLAT                 : -52.243 deg\n",
      "\n",
      "Semimajor (68%)      : 0.0049 deg\n",
      "Semiminor (68%)      : 0.0049 deg\n",
      "Position angle (68%) : -22.82 deg\n",
      "Semimajor (95%)      : 0.0079 deg\n",
      "Semiminor (95%)      : 0.0079 deg\n",
      "Position angle (95%) : -22.82 deg\n",
      "ROI number           : 1355\n",
      "\n",
      "*** Spectral info ***\n",
      "\n",
      "Spectrum type                                 : LogParabola      \n",
      "Detection significance (100 MeV - 1 TeV)      : 224.639\n",
      "beta                                          : 0.0410 +- 0.00412\n",
      "Significance curvature                        : 5.0\n",
      "Pivot energy                                  : 1119 MeV\n",
      "Spectral index                                : 1.767 +- 0.012\n",
      "Flux Density at pivot energy                  : 1.44e-11 +- 1.66e-13 cm-2 MeV-1 s-1\n",
      "Integral flux (1 - 100 GeV)                   : 2.05e-08 +- 2.39e-10 cm-2 s-1\n",
      "Energy flux (100 MeV - 100 GeV)               : 2.16e-10 +- 3.3e-12 erg cm-2 s-1\n",
      "\n",
      "*** Spectral points ***\n",
      "\n",
      "  e_min     e_max        flux     flux_errn   flux_errp      e2dnde     e2dnde_errn   e2dnde_errp  is_ul   flux_ul     e2dnde_ul   sqrt_ts\n",
      "   MeV       MeV     1 / (cm2 s) 1 / (cm2 s) 1 / (cm2 s) erg / (cm2 s) erg / (cm2 s) erg / (cm2 s)       1 / (cm2 s) erg / (cm2 s)        \n",
      "--------- ---------- ----------- ----------- ----------- ------------- ------------- ------------- ----- ----------- ------------- -------\n",
      "   50.000    100.000   8.811e-08   3.040e-08   3.097e-08     1.438e-11     4.960e-12     5.054e-12 False         nan           nan   2.534\n",
      "  100.000    300.000   6.859e-08   3.021e-09   3.021e-09     1.713e-11     7.543e-13     7.543e-13 False         nan           nan  25.058\n",
      "  300.000   1000.000   3.349e-08   6.226e-10   6.226e-10     2.380e-11     4.425e-13     4.425e-13 False         nan           nan  88.745\n",
      " 1000.000   3000.000   1.275e-08   2.198e-10   2.198e-10     3.124e-11     5.386e-13     5.386e-13 False         nan           nan 128.293\n",
      " 3000.000  10000.000   5.410e-09   1.249e-10   1.249e-10     3.760e-11     8.681e-13     8.681e-13 False         nan           nan 120.772\n",
      "10000.000  30000.000   1.769e-09   6.803e-11   6.803e-11     4.254e-11     1.636e-12     1.636e-12 False         nan           nan  82.057\n",
      "30000.000 300000.000   6.966e-10   4.141e-11   4.141e-11     3.530e-11     2.099e-12     2.099e-12 False         nan           nan  54.027\n",
      "*** Lightcurve info ***\n",
      "\n",
      "Lightcurve measured in the energy band: 100 MeV - 100 GeV\n",
      "\n",
      "Variability index : 656.054\n",
      "Significance peak (100 MeV - 100 GeV)    : 98.236\n",
      "Integral flux peak (100 MeV - 100 GeV)   : 1.71e-07 +- 4.01e-09 cm^-2 s^-1\n",
      "Time peak                                : 2.55e+08 s (Mission elapsed time)\n",
      "Peak interval                            : 3.63e+02 day\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(source)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5400"
      ]
     },
     "execution_count": 101,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source.row_index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<i>Row index=5400</i>\n",
       "<table id=\"table140652500143520\">\n",
       "<thead><tr><th>Source_Name</th><th>DataRelease</th><th>RAJ2000</th><th>DEJ2000</th><th>GLON</th><th>GLAT</th><th>Conf_68_SemiMajor</th><th>Conf_68_SemiMinor</th><th>Conf_68_PosAng</th><th>Conf_95_SemiMajor</th><th>Conf_95_SemiMinor</th><th>Conf_95_PosAng</th><th>ROI_num</th><th>Extended_Source_Name</th><th>Signif_Avg</th><th>Pivot_Energy</th><th>Flux1000</th><th>Unc_Flux1000</th><th>Energy_Flux100</th><th>Unc_Energy_Flux100</th><th>SpectrumType</th><th>PL_Flux_Density</th><th>Unc_PL_Flux_Density</th><th>PL_Index</th><th>Unc_PL_Index</th><th>LP_Flux_Density</th><th>Unc_LP_Flux_Density</th><th>LP_Index</th><th>Unc_LP_Index</th><th>LP_beta</th><th>Unc_LP_beta</th><th>LP_SigCurv</th><th>PLEC_Flux_Density</th><th>Unc_PLEC_Flux_Density</th><th>PLEC_Index</th><th>Unc_PLEC_Index</th><th>PLEC_Expfactor</th><th>Unc_PLEC_Expfactor</th><th>PLEC_Exp_Index</th><th>Unc_PLEC_Exp_Index</th><th>PLEC_SigCurv</th><th>Npred</th><th>Flux_Band [7]</th><th>Unc_Flux_Band [7,2]</th><th>nuFnu_Band [7]</th><th>Sqrt_TS_Band [7]</th><th>Variability_Index</th><th>Frac_Variability</th><th>Unc_Frac_Variability</th><th>Signif_Peak</th><th>Flux_Peak</th><th>Unc_Flux_Peak</th><th>Time_Peak</th><th>Peak_Interval</th><th>Flux_History [10]</th><th>Unc_Flux_History [10,2]</th><th>Sqrt_TS_History [10]</th><th>ASSOC_FGL</th><th>ASSOC_FHL</th><th>ASSOC_GAM1</th><th>ASSOC_GAM2</th><th>ASSOC_GAM3</th><th>TEVCAT_FLAG</th><th>ASSOC_TEV</th><th>CLASS1</th><th>CLASS2</th><th>ASSOC1</th><th>ASSOC2</th><th>ASSOC_PROB_BAY</th><th>ASSOC_PROB_LR</th><th>RA_Counterpart</th><th>DEC_Counterpart</th><th>Unc_Counterpart</th><th>Flags</th></tr></thead>\n",
       "<thead><tr><th></th><th></th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th>deg</th><th></th><th></th><th></th><th>MeV</th><th>1 / (cm2 s)</th><th>1 / (cm2 s)</th><th>erg / (cm2 s)</th><th>erg / (cm2 s)</th><th></th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th></th><th></th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th></th><th></th><th></th><th></th><th></th><th>1 / (cm2 MeV s)</th><th>1 / (cm2 MeV s)</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>1 / (cm2 s)</th><th>1 / (cm2 s)</th><th>erg / (cm2 s)</th><th></th><th></th><th></th><th></th><th></th><th>1 / (cm2 s)</th><th>1 / (cm2 s)</th><th>s</th><th>s</th><th>1 / (cm2 s)</th><th>1 / (cm2 s)</th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th></th><th>deg</th><th>deg</th><th>deg</th><th></th></tr></thead>\n",
       "<thead><tr><th>bytes18</th><th>int16</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>int16</th><th>bytes18</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>bytes17</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>float64</th><th>float32</th><th>float32</th><th>float32</th><th>float32</th><th>bytes18</th><th>bytes18</th><th>bytes18</th><th>bytes18</th><th>bytes18</th><th>bytes1</th><th>bytes26</th><th>bytes5</th><th>bytes5</th><th>bytes28</th><th>bytes26</th><th>float32</th><th>float32</th><th>float64</th><th>float64</th><th>float32</th><th>int16</th></tr></thead>\n",
       "<tr><td>4FGL J2158.8-3013</td><td>1</td><td>329.7141</td><td>-30.2251</td><td>17.7308</td><td>-52.2434</td><td>0.0049</td><td>0.0049</td><td>-22.820</td><td>0.0079</td><td>0.0079</td><td>-22.820</td><td>1355</td><td></td><td>224.639</td><td>1118.64</td><td>2.0495e-08</td><td>2.3869e-10</td><td>2.1621e-10</td><td>3.2972e-12</td><td>LogParabola</td><td>1.3735e-11</td><td>1.4466e-13</td><td>1.8507</td><td>0.0071</td><td>1.4365e-11</td><td>1.6555e-13</td><td>1.7668</td><td>0.0116</td><td>0.0410</td><td>0.0041</td><td>5.024</td><td>1.3951e-11</td><td>1.5212e-13</td><td>1.7437</td><td>0.0127</td><td>0.00047</td><td>0.00005</td><td>0.6667</td><td>--</td><td>5.391</td><td>15833.77</td><td>8.811208e-08 .. 6.9656847e-10</td><td>-3.0399406e-08 .. 4.141124e-11</td><td>1.4376523e-11 .. 3.529928e-11</td><td>2.5337915 .. 54.02705</td><td>656.0535</td><td>0.2555</td><td>0.0610</td><td>98.236</td><td>1.7111846e-07</td><td>4.005396e-09</td><td>255223712.0</td><td>31332582.0</td><td>1.7111846e-07 .. 1.1275729e-07</td><td>-4.005396e-09 .. 3.5797774e-09</td><td>98.23613 .. 69.82156</td><td>3FGL J2158.8-3013</td><td>3FHL J2158.8-3013</td><td>J2157-3015</td><td>3EG J2158-3023</td><td>EGR J2200-3015</td><td>P</td><td>PKS 2155-304</td><td>bll</td><td></td><td>PKS 2155-304</td><td></td><td>1.0000</td><td>0.9997</td><td>329.7169</td><td>-30.2256</td><td>0.00000</td><td>0</td></tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<Row index=5400>\n",
       "   Source_Name     DataRelease RAJ2000  DEJ2000    GLON     GLAT   Conf_68_SemiMajor Conf_68_SemiMinor Conf_68_PosAng Conf_95_SemiMajor Conf_95_SemiMinor Conf_95_PosAng ROI_num Extended_Source_Name Signif_Avg Pivot_Energy   Flux1000  Unc_Flux1000 Energy_Flux100 Unc_Energy_Flux100    SpectrumType   PL_Flux_Density Unc_PL_Flux_Density PL_Index Unc_PL_Index LP_Flux_Density Unc_LP_Flux_Density LP_Index Unc_LP_Index LP_beta  Unc_LP_beta LP_SigCurv PLEC_Flux_Density Unc_PLEC_Flux_Density PLEC_Index Unc_PLEC_Index PLEC_Expfactor Unc_PLEC_Expfactor PLEC_Exp_Index Unc_PLEC_Exp_Index PLEC_SigCurv   Npred            Flux_Band [7]              Unc_Flux_Band [7,2]               nuFnu_Band [7]           Sqrt_TS_Band [7]   Variability_Index Frac_Variability Unc_Frac_Variability Signif_Peak   Flux_Peak   Unc_Flux_Peak  Time_Peak  Peak_Interval       Flux_History [10]           Unc_Flux_History [10,2]     Sqrt_TS_History [10]     ASSOC_FGL          ASSOC_FHL          ASSOC_GAM1         ASSOC_GAM2         ASSOC_GAM3     TEVCAT_FLAG         ASSOC_TEV          CLASS1 CLASS2            ASSOC1                      ASSOC2           ASSOC_PROB_BAY ASSOC_PROB_LR RA_Counterpart DEC_Counterpart Unc_Counterpart Flags\n",
       "                                 deg      deg      deg      deg           deg               deg             deg              deg               deg             deg                                                   MeV      1 / (cm2 s) 1 / (cm2 s)  erg / (cm2 s)    erg / (cm2 s)                      1 / (cm2 MeV s)   1 / (cm2 MeV s)                         1 / (cm2 MeV s)   1 / (cm2 MeV s)                                                          1 / (cm2 MeV s)     1 / (cm2 MeV s)                                                                                                                                   1 / (cm2 s)                   1 / (cm2 s)                   erg / (cm2 s)                                                                                                    1 / (cm2 s)   1 / (cm2 s)       s            s                1 / (cm2 s)                    1 / (cm2 s)                                                                                                                                                                                                                                                                              deg             deg             deg            \n",
       "     bytes18          int16    float32  float32  float32  float32       float32           float32         float32          float32           float32         float32      int16        bytes18         float32     float32      float32     float32       float32          float32            bytes17          float32           float32       float32    float32        float32           float32       float32    float32    float32    float32    float32        float32             float32         float32      float32        float32          float32          float32          float32         float32     float32              float32                       float32                        float32                   float32             float32          float32            float32          float32      float32       float32      float64      float32               float32                        float32                   float32             bytes18            bytes18            bytes18            bytes18            bytes18          bytes1            bytes26           bytes5 bytes5           bytes28                     bytes26              float32        float32       float64         float64         float32     int16\n",
       "------------------ ----------- -------- -------- -------- -------- ----------------- ----------------- -------------- ----------------- ----------------- -------------- ------- -------------------- ---------- ------------ ----------- ------------ -------------- ------------------ ----------------- --------------- ------------------- -------- ------------ --------------- ------------------- -------- ------------ -------- ----------- ---------- ----------------- --------------------- ---------- -------------- -------------- ------------------ -------------- ------------------ ------------ ---------- ----------------------------- ------------------------------ ----------------------------- --------------------- ----------------- ---------------- -------------------- ----------- ------------- ------------- ----------- ------------- ------------------------------ ------------------------------ -------------------- ------------------ ------------------ ------------------ ------------------ ------------------ ----------- -------------------------- ------ ------ ---------------------------- -------------------------- -------------- ------------- -------------- --------------- --------------- -----\n",
       "4FGL J2158.8-3013            1 329.7141 -30.2251  17.7308 -52.2434            0.0049            0.0049        -22.820            0.0079            0.0079        -22.820    1355                         224.639      1118.64  2.0495e-08   2.3869e-10     2.1621e-10         3.2972e-12 LogParabola            1.3735e-11          1.4466e-13   1.8507       0.0071      1.4365e-11          1.6555e-13   1.7668       0.0116   0.0410      0.0041      5.024        1.3951e-11            1.5212e-13     1.7437         0.0127        0.00047            0.00005         0.6667                 --        5.391   15833.77 8.811208e-08 .. 6.9656847e-10 -3.0399406e-08 .. 4.141124e-11 1.4376523e-11 .. 3.529928e-11 2.5337915 .. 54.02705          656.0535           0.2555               0.0610      98.236 1.7111846e-07  4.005396e-09 255223712.0    31332582.0 1.7111846e-07 .. 1.1275729e-07 -4.005396e-09 .. 3.5797774e-09 98.23613 .. 69.82156 3FGL J2158.8-3013  3FHL J2158.8-3013   J2157-3015        3EG J2158-3023     EGR J2200-3015               P PKS 2155-304                bll          PKS 2155-304                                                    1.0000        0.9997       329.7169        -30.2256         0.00000     0"
      ]
     },
     "execution_count": 102,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "catalog.table[5400]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<SkyCoord (ICRS): (ra, dec) in deg\n",
       "    (329.71411133, -30.22509956)>"
      ]
     },
     "execution_count": 103,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "source.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "metadata": {},
   "outputs": [],
   "source": [
    "sky_model = source.sky_model()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SkyModel\n",
      "\n",
      "  Name                      : 4FGL J2158.8-3013\n",
      "  Datasets names            : None\n",
      "  Spectral model type       : LogParabolaSpectralModel\n",
      "  Spatial  model type       : PointSpatialModel\n",
      "  Temporal model type       : \n",
      "  Parameters:\n",
      "    amplitude               :   1.44e-11   +/- 1.7e-13 1 / (cm2 MeV s)\n",
      "    reference    (frozen)   :   1118.643       MeV         \n",
      "    alpha                   :      1.767   +/-    0.01             \n",
      "    beta                    :      0.041   +/-    0.00             \n",
      "    lon_0                   :    329.714   +/-    0.00 deg         \n",
      "    lat_0                   :    -30.225   +/-    0.00 deg         \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(sky_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:xlabel='Energy [MeV]', ylabel='e2dnde (erg / (cm2 s))'>"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaAAAAEKCAYAAABUsYHRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa5UlEQVR4nO3de7RdZXnv8e8vITSGYDiEzTiUkOYKgiik7OChSEFrECkBBeTajiAZ4WLBWmy5iK0XbLH2jBQZIBcJJ0etCQHlQCiCKLIRG2MSJCEQkAAh7IN2S0ICMVwS8vSPNTesLtZl7r3XXHNdfp8x9iDvO29PXjc8PnO+c76KCMzMzBptWN4BmJlZZ3ICMjOzXDgBmZlZLpyAzMwsF05AZmaWCycgMzPLxU55B9Aq9thjj5gwYULeYZiZtZQVK1a8GBFd5bY5AaU0YcIEli9fnncYZmYtRdJzlbb5FpyZmeXCCcjMzHLhBGRmZrlwAjIzs1w4AZmZWS6cgMzMLBdOQGZmlgsnIDNreqfesIRTb1iSdxhWZ05AZmaWCycgMzPLhROQmTW19Ru2srJ3E0uf3ciMuT2s37A175CsTvwtODOru3o+r1nZu4nXtu0A4Km+LRx9VQ8HjdutLue+5dzD6nIeGxxXQGbW1PqTT6W2tS5XQGZWd/WsLGbM7eGpvi0ADBNM7hrtyqVNuAIys6Y2b9Z0Ro4o/Kdqctdo5s2annNEVi+ugMysqY0fO+qtZz6ufNqLKyAzM8tFR1VAkiYBlwNjIuLk0na+0ZlZJa582lPmFZCk4ZJ+JemuIZzjZkl9klaX2XaMpCclrZV0abXzRMQzETG7UtvMzBqnEbfg/hpYU26DpD0l7VrSN6XMrvOBY8ocPxy4FvgYcABwuqQDJL1P0l0lP3sO9S9iZmb1k2kCkjQO+HPgpgq7HAncIWlksv8c4OrSnSLiQWBjmeMPBdYmlcwbwELghIh4NCKOK/npq8ffyczM6iPrCugq4GKg7JtjEXErcA+wUNKZwNnAKQM4/97A80Xt3qSvLEljJV0PTJN0WWm7wjEzJd24efPmAYRlZma1ZDYJQdJxQF9ErJB0VKX9IuLrkhYC1wGTI2LLQC5T7pRVrrUBOK+ku7RdesxiYHF3d/ecAcRlZmY1ZFkBHQ4cL2kdhVtjH5b03dKdJB0BHAjcDnxxgNfoBfYpao8DXhhUtGZm1lCZJaCIuCwixkXEBOA04P6I+IvifSRNA74FnAB8Cthd0lcHcJllwFRJEyXtnFznzrr8BczMLFN5v4g6CvhkRDwdETuAWcBzpTtJWgAsAfaT1CtpNkBEbAcuAO6lMNNuUUQ81rDozVqUVxi1ZqCIio9MrEh3d3csX7487zDMhmz9hq0cfVUPr23bwdQ9C99WGz92VN5hWZuStCIiustt66gvIZi1miyqlCzX1+nnLxdko//3oV3GN+9bcGbWYF5fx5qFKyCzJpbF/9P1+jrWLFwBmXUYr69jzcIVkFmHGT92FE9c8bG8wzBzBWRm1grWb9jKyt5NLH12IzPm9rB+w9a8QxoyV0Bm1rFa6V2oRsxebDRXQGZmLaAdZy/6RdSU/CKqmeWp3OzF+y46Mueoaqv2IqorIDOzFtCOsxf9DMjMrAWMHzvqrWc+7fLelisgMzPLhROQmZnlwgnIzMxy4QRkZma58CQEM7MW0S6TD/q5AjIzs1w4AZmZWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5cIJyMzMcuH3gMyGaP2GrRx9VQ+vbdvB1D0LXykeP3ZU3mGZNT0nIOs49V4FM8uVKtvtxUOzYjUTkKRu4AjgD4FXgdXAjyNiY8axmbWEdlyp0qwRKiYgSWcBnwGeBVYATwIjgQ8Cl0haDfx9RKxvQJxmdVPvqqLcSpWuXMxqq1YB7QIcHhGvltso6WBgKtAyCUjSJOByYExEnFzazjc6a1XzZk1n9v9dxjO/+z2TunZpi5UqzRpBEZHNiaWRwIPAH1BIdLdFxBcHea6bgeOAvog4sGTbMcA3gOHATRHxtRTnu6044ZS2y+nu7o7ly5cPJnwzs44laUVEdJfbVvUZUJJEjuOdz4D+PSIeq3Hd14EPR8QWSSOAhyT9MCJ+UXT+PYFXI+KVor4pEbG25FzzgWuAb5fENxy4FpgB9ALLJN1JIRldWXKOsyOir0bMZmbWINWeAX0JmAk8ACwF+ig8A9oX+FqSnD4XEavKHR+F0mpL0hyR/JSWW0cC50s6NiJekzQH+ARwbMm5HpQ0ocxlDgXWRsQzScwLgRMi4koKiXPIJM0EZk6ZMqUepzMzs0S1F1GXRcQhEfG5iPheRPw4Iu6KiLkRMRM4E9i52sklDZf0CIXkdV9ELC3eHhG3AvcACyWdCZwNnDKA+PcGni9q9yZ9leIZK+l6YJqky0rb5Y6JiMURcc6YMWMGEJaZmdVSsQKKiH8v7ZM0DBgdES8nt7Oq3tKKiDeBgyXtBtwu6cCIWF2yz9eTyuU6YHJEbClzqkpU7rJV4tkAnFfSXdo2M7MGqPkpHknfk/RuSbsAjwNPSvq7gVwkIjZRuJV3TJnzHwEcCNwODHSSQi+wT1F7HPDCAM9hZmY5SPMtuAMi4mXg48DdwHjgL2sdJKkrqXyQ9C7gI8ATJftMA74FnAB8Cthd0lcHEP8yYKqkiZJ2Bk4D7hzA8WZmlpM0CWhEMovt48AdEbGNKre5iuwF/FTSKgqJ4r6IuKtkn1HAJyPi6YjYAcwCnis9kaQFwBJgP0m9kmYDRMR24ALgXmANsCjF7DwzM2sCab4FdwOwDlgJPCjpj4CXax2UzI6bVmOfn5e0t1GoiEr3O73KOe6mUJmZmVkLqVkBRcTVEbF3RBybTK1eD3wo+9DMzKydDfhr2EkS2p5BLGZm1kG8IJ2ZmeXCCcjMzHJRNQEl7/9MLtP//uxCMjOzTlAxAUk6hcJ7O9+X9Jik4m/Mz886MDMza2/VKqDPA4dExMEUXhL9jqQTk23lPoFjZmaWWrVZcMMj4jcAEfFLSR8C7pI0jnQvopqZmVVUrQJ6pfj5T5KMjqLw2Zz3ZhyXmZm1uWoV0PmU3GqLiFeSFUgHsmSCmZnZO1RbjmFlcVvSu4v2/2GWQZmZWfur+SUESecCX6GwHHf/s58AJmUYl5mZtbk0n+L5W+C9EfFi1sGYmVnnSPMlhKeBrVkHYmZmnSVNBXQZ8B+SlgKv93dGxGcyi8rMzNpe2vWA7gceBXZkG46ZmXWKNAloe0RclHkkZmbWUdI8A/qppHMk7SVp9/6fzCMzM7O2lqYCOiP552VFfZ6GbWZmQ1IzAUXExEYEYmZmnaXmLThJfyVpt6L2/5D06Uyjso526g1LOPWGJXmHYWYZS/MMaE5EbOpvRMRLwJzMIrKOtn7DVlb2bmLpsxuZMbeH9Rv8CppZu0rzDGiYJEVEAEgaDuycbVjWCrKoUlb2buK1bYXZ/k/1beHoq3o4aNxudb3GLeceVtfzmdngpElA9wKLJF1PYfLBecA9mUZlHas/+VRqm1n7SJOALgHO4e3lGX4E3JRlUNYasqgkZszt4am+LQAME0zuGu2KxaxN1XwGFBE7IuL6iDg5Ik6KiBsi4s1GBGedZ96s6UzdczTDJSZ3jWberOl5h2RmGalYAUlaDNwI3BMR20q2TQLOAtZFxM2ZRmgdZfzYUdx30ZF5h2FmDVDtFtwc4CLgKkkbgd8BI4GJwFrgmoi4I/sQzcysHVVbEfW3wMXAxZImAHtRWJTu1xHhubFmZjYkaSYhEBHrgHWZRmJmZh0lzYuoZmZmdecEZGZmuaiYgCTdKOkTknZtZEBmZtYZqlVANwMHAXdL+omkSyQd1KC4zMyszVWbBfcL4BfAlySNBY4GPifpfcCvKLwftKgxYZqZWbtJOwtuA7Ag+UHSIcAxGcZlZmZtLlUCKhURK4AVdY7FzMw6iGfBmZlZLpyAzMwsFzVvwUk6sUz3ZuDRiOirf0hmZtYJ0jwDmg0cBvw0aR9FYXbcvpK+EhHfySg2MzNrY2luwe0A9k/WAjoJOAB4HfgAhcXqWoakSZLmSbqtXNvMzBonTQKaEBH/WdTuA/aNiI3AtgrHIGkfST+VtEbSY5L+erBBSrpZUp+k1WW2HSPpSUlrJV1a7TwR8UxEzK7UNjOzxklzC+5nku4Cbk3aJwEPStoF2FTluO3A5yLi4eRzPisk3RcRj/fvIGlP4NWIeKWob0pErC0513zgGuDbxZ2ShgPXAjOAXmCZpDuB4cCVJec428+szMyaR80EFBGflnQS8EFAFJLA9yMigA9VOe43wG+SP78iaQ2wN/B40W5HAudLOjYiXpM0B/gEcGzJuR5M1iQqdSiwNiKeAZC0EDghIq4Ejqv1dzMzs/xUTUCShgGrIuJA4PuDvUiSPKYBS4v7I+JWSROBhZJuBc6mUM2ktTfwfFG7l8KzqUpxjAX+EZgm6TIKS46/1U4SV+kxM4GZU6ZMGUBYZmZWS9UEFBE7JK2UND4i1g/mApJGU0hen42Il8tc4+tJ5XIdMDkitgzk9OXCrrRz8kmh80q6S9ulxywGFnd3d88ZQFxmZlZDmmdAewGPSfol8Pv+zog4vtaBkkZQSD7/FhE/qLDPEcCBwO3AF4ELUsTUrxfYp6g9DnhhAMebmVlO0iSgLw/mxJIEzAPWRMTcCvtMA74F/DnwLPBdSV+NiC+kvMwyYGpyG+//A6cBZwwmXjMza6ya07AjogdYB4xI/rwMeDjFuQ8H/hL4sKRHkp9jS/YZBXwyIp6OiB3ALOC50hNJWgAsAfaT1CtpdhLbdgoV073AGmBRRDyWIraGOPWGJZx6w5K8wzAza0ppPsUzBzgH2B2YTOHB//XAn1U7LiIeovwzmuJ9fl7S3kahIird7/Qq57gbuLvadczMrPmkeRH1ryhUMy8DRMRTwJ5ZBmVmZu0vTQJ6PSLe6G9I2okqM83MzMzSSJOAeiR9HniXpBkUvoiwONuwzMys3aVJQJcCvwMeBc6l8Lwl7Sw1MzOzstJ8imcHhYkB75gcYGZmNlgVKyBJiyXNTF4mLd02SdJXJJ2dbXhmZtauqlVAc4CLgKskbaRwG24kMBFYC1wTEXdkH6KZmbWjigkoIn4LXAxcnHxMdC/gVeDXEbG1MeGZmVm7SvMpHiJiHYWvIZiZmdVFmllw1uL8SSAza0ZOQGZmlotUCUjSuyTtl3UwZmbWOWomoGRF0EeAe5L2wZLuzDguq5P1G7aysncTS5/dyIy5Pazf4PkjZtYc0kxC+BJwKPAAQEQ8ksyKs4zU83nNyt5NvLZtBwBP9W3h6Kt6OGjcbnU59y3nHlaX85hZZ0pzC257RGzOPBLLRH/yqdQ2M8tLmgpotaQzgOGSpgKfAf4j27A6Wz0rixlze3iqbwsAwwSTu0a7cjGzppCmAroQeC/wOrCAwrpAn80wJqujebOmM3JE4X/myV2jmTdres4RmZkVpPkY6Vbg8uTHWsz4saPeeubjysfMmknFBCRpMVUWnouI4zOJyMzMOkK1Cuh/J/88EfifwHeT9un4szxmZjZE1T5G2gMg6YqI+NOiTYslPZh5ZFY3vvVmZs0ozSSELkmT+huSJgJd2YVkZmadIM007L8BHpD0TNKeQGFpbjMzs0FLMwvunuT9n/ckXU9ExOvZhmVmZu0u1XpAwCEUKp+dgIMkERHfziwqMzNrezUTkKTvAJMpfJD0zaQ7ACcgMzMbtDQVUDdwQERUfCfIzMxsoNLMgltN4T0gMzOzuklTAe0BPC7plxS+Bwf4SwhmZjY0adcDMjMzq6s007B7GhGImZl1lmofI32F6h8jfXcmEZmZWUeo9i24XQEkfQX4LfAdQMCZwK4Nic7MzNpWmllwH42Ib0bEKxHxckRcB5yUdWBmZtbe0iSgNyWdKWm4pGGSzuTtF1KtgvUbtrKydxNLn93IjLk9rN+wNe+QzMyaimq9XyppAvAN4HAKz4R+Dnw2ItZlHVwz2f2P9o8Zn7859f4rezfx2rYdb7VHjhj21sqktXj5BDNrF5JWRER3uW1pZsGtA06od1Dtrjj5lGubmXW6qglI0keBccCPI+K5ov6zIyJ9OdAGJnXtMqDKZMbcHp7q2wLAMMHkrtGubMzMilR8BiTpn4DLgfcB90u6sGjzBVkH1urmzZrOyBGF4Z3cNZp5s6bnHJGZWXOpVgHNBKZFxHZJXwK+J2lSRPwNhenYVsX4saPeeubjysfM7J2qzYLbKSK2A0TEJgoJ6d2SbgV2bkBsZmbWxqoloKclHdnfiIg3I2I28CSwf+aRmZlZW6uWgD4J/LK0MyK+AOyTWURmZtYRKiagiHg1Il5VwV9I+gcASeOBvRsWoZmZtaU0X0L4JnAYcHrSfgW4NrOIzMysI6RZD+gDEfHHkn4FEBEvSfIkBDMzG5I0FdA2ScNJlmaQ1AX4tX4zMxuSNAnoauB2YE9J/wg8BPxTplGZmVnbS/MtuH+TtAL4MwovoH48ItZkHpmZmbW1NM+AiIgngCcyjsXMzDpImltwZmZmdecEZGZmuXACMjOzXDgBmZlZLpyAzMwsF05AZmaWCycgMzPLhROQmZnlwgnIzMxy4QRkZma5cAIyM7NcOAGZmVkunIDMzCwXTkBmZpYLJyAzM8uFE5CZmeXCCcjMzHLhBGRmZrlwAjIzs1w4AZmZWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5cIJyMzMcuEEZGZmuXACMjOzXDgBmZlZLpyAzMwsFx2RgCRNkjRP0m3V+szMrHGaPgFJullSn6TVJf3HSHpS0lpJl1Y7R0Q8ExGza/WZmVnj7JR3ACnMB64Bvt3fIWk4cC0wA+gFlkm6ExgOXFly/NkR0deYUM3MLK2mT0AR8aCkCSXdhwJrI+IZAEkLgRMi4krguHpdW9I5wDkA48ePr9dpzcyMFrgFV8HewPNF7d6kryxJYyVdD0yTdFmlvlIRcWNEdEdEd1dXVx3DNzOzpq+AKlCZvqi0c0RsAM6r1WdmZo3TqhVQL7BPUXsc8EJOsZiZ2SC0agJaBkyVNFHSzsBpwJ05x2RmZgPQ9AlI0gJgCbCfpF5JsyNiO3ABcC+wBlgUEY/lGaeZmQ1M0z8DiojTK/TfDdzd4HDMzKxOmr4CMjOz9tT0FVAru+Xcw/IOwcysabkCMjOzXDgBmZlZLpyAzMwsF05AZmaWCycgMzPLhRNQDZJmSrpx8+bNeYdiZtZWnIBqiIjFEXHOmDFj8g7FzKytOAGZmVkunIDMzCwXiqi4jI4VkbQZeKrC5jFAuYdE5frT9O0BvDiIMAejUuz1PjbNvtX2GcoYt+r4DvT4WvvWa3zL9bfqGLfK73C5vkaN8VB/h6dGRPlnGBHhnxQ/wI0D3VauP00fsLwZ/l71PDbNvlmNcauOb73HuF7j205j3Cq/w3mOcZa/w74Fl97iQWwr15+2r1GGcu2BHJtm36zGuFXHd6DH19q3XuNbrr9Vx7hVfofTXj8Lmf0O+xZcE5K0PCK6846jXXl8s+cxzl47jLEroOZ0Y94BtDmPb/Y8xtlr+TF2BWRmZrlwBWRmZrlwAjIzs1w4AZmZWS6cgFqApI9L+pakOyQdnXc87UbS/pKul3SbpPPzjqddSdpF0gpJx+UdS7uRdJSknyW/x0flHU9aTkA5kXSzpD5Jq0v6j5H0pKS1ki4FiIj/FxFzgLOAU3MIt+UMcHzXRMR5wClAS09rbaSBjHHiEmBRY6NsXQMc3wC2ACOB3kbHOlhOQPmZDxxT3CFpOHAt8DHgAOB0SQcU7fKFZLvVNp8BjK+k44GHgJ80NsyWNp+UYyzpI8DjwH82OsgWNp/0v8M/i4iPUUjyX25wnIPmBJSTiHgQ2FjSfSiwNiKeiYg3gIXACSr4Z+CHEfFwo2NtRQMZ32T/OyPiT4AzGxtp6xrgGH8I+F/AGcAcSf5vTw0DGd+I2JFsfwn4gwaGOSQ75R2A/Td7A88XtXuBDwAXAh8BxkiaEhHX5xFcGyg7vsk98xMp/It7d+PDaitlxzgiLgCQdBbwYtF/MG1gKv0Onwh8FNgNuCaHuAbFCai5qExfRMTVwNWNDqYNVRrfB4AHGhtK2yo7xm/9IWJ+40JpS5V+h38A/KDRwQyVy+Dm0gvsU9QeB7yQUyztyOObPY9xttpqfJ2AmssyYKqkiZJ2Bk4D7sw5pnbi8c2exzhbbTW+TkA5kbQAWALsJ6lX0uyI2A5cANwLrAEWRcRjecbZqjy+2fMYZ6sTxtcfIzUzs1y4AjIzs1w4AZmZWS6cgMzMLBdOQGZmlgsnIDMzy4UTkJmZ5cIJyGyIJL0p6ZGin0trH5W9orj+MGmvk/Szkn0eKf3cf5nzPCtpv5K+qyRdLOkISY/XOodZOf4WnNnQvRoRB9fzhJJ2Sl46HIpyce0qaZ+IeF7S/inPs5DCG/dfTmIbBpwMHB4Rz0k6FrhriLFaB3IFZJaRpOL4sqSHJT0q6T1J/y7JYmPLJP1K0glJ/1mSbpW0GPiRpFGSFklaJekWSUsldUuaLelfi64zR9LclGEt4u1FDU8HFhSdZ7ikf0niWiXp3GTTAgoJqN+fAusi4rlBDYxZwgnIbOjeVXILrnjV2hcj4o+B64C/TfouB+6PiOkU1sn5F0m7JNsOA2ZFxIeBTwMvRcT7gSuAQ5J9FgLHSxqRtD8F/J+Usd5GYekJgJnA4qJts4HNSVzTKazbMzEiVgE7JB2U7HcaRYnLbLB8C85s6Krdguv/RP4K3v4P/9EUEkh/QhoJjE/+fF9E9C9C9kHgGwARsVrSquTPv5d0P3CcpDXAiIh4NGWsG4GXJJ1G4VtiW4u2HQ28X9LJSXsMMBV4lqQKkvQYhQXm/iHl9cwqcgIyy9bryT/f5O1/3wScFBFPFu8o6QPA74u7qpz3JuDzwBOkr3763UJhWeezSvoFXBgR95Y5ZgHwI6AHWBURfQO8ptk7+BacWePdC1woSQCSplXY7yHglGSfA4D39W+IiKUU1oU5g4HfDrsd+HoSR2lc5/ff2pO0b/+twYh4GtgAfG0Q1zMrywnIbOhKnwF9rcb+VwAjgFXJ9OUrKuz3TaArufV2CbAK2Fy0fRHw84h4aSDBRsQrEfHPEfFGyaabgMeBh5O4buC/3yVZALyHQgIzGzIvx2DWpCQNp/B85zVJk4GfAPv2Jw5JdwH/GhE/qXD8logY3YA4JwB3RcSBWV/L2osrILPmNQp4SNJKClXH+RHxhqTdJP2awuSHsskn8XLxi6hZkHQEhZl0L2Z1DWtfroDMzCwXroDMzCwXTkBmZpYLJyAzM8uFE5CZmeXCCcjMzHLhBGRmZrn4L+IYUO1dGapPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "source.flux_points.plot(sed_type=\"e2dnde\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Now, simulate for CTA using the above model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create a fake observation defining observation conditions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Invalid unit found in background table! Assuming (s-1 MeV-1 sr-1)\n"
     ]
    }
   ],
   "source": [
    "# Loading IRFs\n",
    "irfs = load_cta_irfs(\n",
    "    \"$GAMMAPY_DATA/cta-1dc/caldb/data/cta/1dc/bcf/South_z20_50h/irf_file.fits\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the observation parameters (typically the observation duration and the pointing position):\n",
    "livetime = 2.0 * u.hr\n",
    "pointing = source.position"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define map geometry for binned simulation\n",
    "energy_reco = MapAxis.from_edges(\n",
    "    np.logspace(-1.0, 1.0, 10), unit=\"TeV\", name=\"energy\", interp=\"log\"\n",
    ")\n",
    "geom = WcsGeom.create(\n",
    "    skydir=source.position,\n",
    "    binsz=0.02,\n",
    "    width=(3, 3),\n",
    "    frame=\"icrs\",\n",
    "    axes=[energy_reco],\n",
    ")\n",
    "# It is usually useful to have a separate binning for the true energy axis\n",
    "energy_true = MapAxis.from_edges(\n",
    "    np.logspace(-1.5, 1.5, 30), unit=\"TeV\", name=\"energy_true\", interp=\"log\"\n",
    ")\n",
    "\n",
    "empty = MapDataset.create(\n",
    "    geom, name=\"dataset-simu\", energy_axis_true=energy_true\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define model to simulate the data\n",
    "\n",
    "- intrinsic model from the 4FGL\n",
    "- EBL model (see: https://docs.gammapy.org/0.19/modeling/gallery/spectral/plot_absorbed.html?highlight=ebl)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "LogParabolaSpectralModel\n",
      "\n",
      "  type      name     value         unit        error   min max frozen link\n",
      "-------- --------- ---------- -------------- --------- --- --- ------ ----\n",
      "spectral amplitude 1.4365e-11 cm-2 MeV-1 s-1 1.655e-13 nan nan  False     \n",
      "spectral reference 1.1186e+03            MeV 0.000e+00 nan nan   True     \n",
      "spectral     alpha 1.7668e+00                1.163e-02 nan nan  False     \n",
      "spectral      beta 4.0969e-02                4.120e-03 nan nan  False     \n"
     ]
    }
   ],
   "source": [
    "int_model = sky_model.spectral_model\n",
    "print(int_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [],
   "source": [
    "redshift = 0.116\n",
    "ebl = EBLAbsorptionNormSpectralModel.read_builtin(\"dominguez\", redshift=redshift)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectral_model = ebl * int_model\n",
    "spatial_model = sky_model.spatial_model\n",
    "model_simu = SkyModel(spectral_model=spectral_model, \n",
    "                      spatial_model=spatial_model,\n",
    "                      name=\"model-simu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "bkg_model = FoVBackgroundModel(dataset_name=\"dataset-simu\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Models\n",
      "\n",
      "Component 0: SkyModel\n",
      "\n",
      "  Name                      : model-simu\n",
      "  Datasets names            : None\n",
      "  Spectral model type       : CompoundSpectralModel\n",
      "  Spatial  model type       : PointSpatialModel\n",
      "  Temporal model type       : \n",
      "  Parameters:\n",
      "    alpha_norm   (frozen)   :      1.000                   \n",
      "    redshift     (frozen)   :      0.116                   \n",
      "    amplitude               :   1.44e-11   +/- 1.7e-13 1 / (cm2 MeV s)\n",
      "    reference    (frozen)   :   1118.643       MeV         \n",
      "    alpha                   :      1.767   +/-    0.01             \n",
      "    beta                    :      0.041   +/-    0.00             \n",
      "    lon_0                   :    329.714   +/-    0.00 deg         \n",
      "    lat_0                   :    -30.225   +/-    0.00 deg         \n",
      "\n",
      "Component 1: FoVBackgroundModel\n",
      "\n",
      "  Name                      : dataset-simu-bkg\n",
      "  Datasets names            : ['dataset-simu']\n",
      "  Spectral model type       : PowerLawNormSpectralModel\n",
      "  Parameters:\n",
      "    norm                    :      1.000   +/-    0.00             \n",
      "    tilt         (frozen)   :      0.000                   \n",
      "    reference    (frozen)   :      1.000       TeV         \n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "models = Models([model_simu, bkg_model])\n",
    "print(models)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now, perform a simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Observation\n",
      "\n",
      "\tobs id            : 0 \n",
      " \ttstart            : 51544.00\n",
      "\ttstop             : 51544.08\n",
      "\tduration          : 7200.00 s\n",
      "\tpointing (icrs)   : 329.7 deg, -30.2 deg\n",
      "\n",
      "\tdeadtime fraction : 0.0%\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Create an in-memory observation\n",
    "obs = Observation.create(pointing=pointing, livetime=livetime, irfs=irfs)\n",
    "print(obs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MapDataset\n",
      "----------\n",
      "\n",
      "  Name                            : dataset-simu \n",
      "\n",
      "  Total counts                    : 0 \n",
      "  Total background counts         : 52917.67\n",
      "  Total excess counts             : -52917.67\n",
      "\n",
      "  Predicted counts                : 52917.67\n",
      "  Predicted background counts     : 52917.67\n",
      "  Predicted excess counts         : nan\n",
      "\n",
      "  Exposure min                    : 4.64e+07 m2 s\n",
      "  Exposure max                    : 3.58e+10 m2 s\n",
      "\n",
      "  Number of total bins            : 202500 \n",
      "  Number of fit bins              : 202500 \n",
      "\n",
      "  Fit statistic type              : cash\n",
      "  Fit statistic value (-2 log(L)) : nan\n",
      "\n",
      "  Number of models                : 0 \n",
      "  Number of parameters            : 0\n",
      "  Number of free parameters       : 0\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Make the MapDataset\n",
    "maker = MapDatasetMaker(selection=[\"exposure\", \"background\", \"psf\", \"edisp\"])\n",
    "\n",
    "maker_safe_mask = SafeMaskMaker(methods=[\"offset-max\"], offset_max=4.0 * u.deg)\n",
    "\n",
    "dataset = maker.run(empty, obs)\n",
    "dataset = maker_safe_mask.run(dataset, obs)\n",
    "print(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MapDataset\n",
      "----------\n",
      "\n",
      "  Name                            : dataset-simu \n",
      "\n",
      "  Total counts                    : 55515 \n",
      "  Total background counts         : 52917.67\n",
      "  Total excess counts             : 2597.33\n",
      "\n",
      "  Predicted counts                : 55665.85\n",
      "  Predicted background counts     : 52917.67\n",
      "  Predicted excess counts         : 2748.17\n",
      "\n",
      "  Exposure min                    : 4.64e+07 m2 s\n",
      "  Exposure max                    : 3.58e+10 m2 s\n",
      "\n",
      "  Number of total bins            : 202500 \n",
      "  Number of fit bins              : 202500 \n",
      "\n",
      "  Fit statistic type              : cash\n",
      "  Fit statistic value (-2 log(L)) : 152592.82\n",
      "\n",
      "  Number of models                : 2 \n",
      "  Number of parameters            : 11\n",
      "  Number of free parameters       : 6\n",
      "\n",
      "  Component 0: SkyModel\n",
      "  \n",
      "    Name                      : model-simu\n",
      "    Datasets names            : None\n",
      "    Spectral model type       : CompoundSpectralModel\n",
      "    Spatial  model type       : PointSpatialModel\n",
      "    Temporal model type       : \n",
      "    Parameters:\n",
      "      alpha_norm   (frozen)   :      1.000                   \n",
      "      redshift     (frozen)   :      0.116                   \n",
      "      amplitude               :   1.44e-11   +/- 1.7e-13 1 / (cm2 MeV s)\n",
      "      reference    (frozen)   :   1118.643       MeV         \n",
      "      alpha                   :      1.767   +/-    0.01             \n",
      "      beta                    :      0.041   +/-    0.00             \n",
      "      lon_0                   :    329.714   +/-    0.00 deg         \n",
      "      lat_0                   :    -30.225   +/-    0.00 deg         \n",
      "  \n",
      "  Component 1: FoVBackgroundModel\n",
      "  \n",
      "    Name                      : dataset-simu-bkg\n",
      "    Datasets names            : ['dataset-simu']\n",
      "    Spectral model type       : PowerLawNormSpectralModel\n",
      "    Parameters:\n",
      "      norm                    :      1.000   +/-    0.00             \n",
      "      tilt         (frozen)   :      0.000                   \n",
      "      reference    (frozen)   :      1.000       TeV         \n",
      "  \n",
      "  \n"
     ]
    }
   ],
   "source": [
    "#Add the model on the dataset and Poission fluctuate\n",
    "dataset.models = models\n",
    "dataset.fake()\n",
    "# Do a print on the dataset - there is now a counts maps\n",
    "print(dataset)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0a30bdd8019043bfad7073681e74c4a9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(SelectionSlider(continuous_update=False, description='Select energy:', layout=Layout(wid…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "dataset.excess.smooth(0.05 * u.deg).plot_interactive(\n",
    "    add_cbar=True, stretch=\"linear\"\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "WcsNDMap\n",
       "\n",
       "\tgeom  : WcsGeom \n",
       " \taxes  : ['lon', 'lat', 'energy']\n",
       "\tshape : (150, 150, 9)\n",
       "\tndim  : 3\n",
       "\tunit  : \n",
       "\tdtype : int64"
      ]
     },
     "execution_count": 134,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.counts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<WCSAxesSubplot:xlabel='Right Ascension', ylabel='Declination'>"
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVIAAAEJCAYAAADGnK/bAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAB15UlEQVR4nO29e5hlV1km/n7d6SKXIpCik3Q6UFZD03hpIYYMKmiGu1EYoygQEU3GyDgzKF5niDz+ZAR1gvP7qTwaZ8wDaABFIhKJDAgxGBiERJIYoBGm09hFETtN6K6QUFxS6e7v98feH+fd7/nWOaf6nK6u6qz3eeqpfVm3vfba3/ou71rH3B0VFRUVFUePDce7ARUVFRXrHVWQVlRUVIyJKkgrKioqxkQVpOsQZmbHuw0V6xN17BwbVEG6PrF0vBswDszsK8e7DeNgnbd/Pbd9zaIK0oqKhxYqTecYoArSioqKhyzM7BfN7FNmtsvM3mZmJ5vZjJndYGZ3tv/PGFpO5ZE22GDm4TxyAHG8AcCR9pidSxsBLLfHJ7Vpoicd3RkqrnP+I20ZgcMAHtYeL0tazR95XdLE/ZPa8rQel2dxKoPL4+tTAA4ldShKz8Zt3NiWFW3QWdyT64fR7Sdtf0DL4r7KytVn4XRc9hSAB+n6pvZ8Q/ssJ7f3DlE7D1HZh6nek6jcKIfLzr5Ea/NFWfwO+ZmO0LGWw+/4MJWnfcP5og4eS9b+HaY82TjV9vP9R87M+MGDB8dS4C666CI/cODA0HS33Xbb+9z9otJ9MzsXwIcBfKu7f83MrgXwHgDfCmDR3a80sysAnOHurxxU10mDbj6UsAnAbHu8gOYDAhqhdjKlm0qOF9EMqun2fAnAlvZ4P4AZuj6DHpbp+jKA06gtUZbm17xo7y8V0nE7lwA8gq4vovfM0YZInx1Pt+2J42VqJ9d7Cl2fovr303GUG3mmC2Utoffh8v1lKXuqfR6g+46WJA33Bbd/mZ7z2wAcoOtLlB9o3k/cO5Xq4vHAbYhyN1FZW9D0R5Q1jW4/LyXHS+35jvZ8N92bKRwzltGdlDZQO2fQjPs4jr48mco6G71+j+ecLuRfpjSBKQDbtm0bO9h14MAB3HrrPw5NZ7Zx8wjFnQTgFDN7EM3r3AfgVwE8vb1/DYCbAAwUpNW0r6ioWGcIvXzY35BS3P8VwP+LZg64G8B97v5+AGe7+91tmrsBnDWsrKqRtjiEnralWKb/rPXtoGPWfHhGj/uRP2Zt1q7iPMvDmgrXz9pAaFeZRsPaMdDVyFiLYw13WdLx9agz2juV5In2oW3vfjrmvuA+0D7j9mSa6ozkYY1ctSDVgqNMfrbQEIFGJeH6uS/4/WsZrN0GZqU9kWYBwHaqky0C7ht+zzEW9lG6aE9JC9R+3kLpuM4FdMF9trM93ofu83J/7kfXCovr/M51zI+H4YISwGYzu5XOr3b3q+Ok9X1eDGAbgC8B+Esze+nRtKYK0haO3BybQvcD2doefx2NaRXQjz0TBHx9Gl0TqeQOYGG1E8Cu9rjkIoDcn0IulPWj53Q7kH+snD7KYIHJ9zh9JvyjThZS3CYWmEjSxHnUuQW5aa/nLPgZPBFyPfrO1OVREswBfn9c7jQa4ZUJ7Egb11RYLyXpeCJT4clCbRm9ftJ3wZMv5z0g6XQMBUIYz6I3QSyi+/63YVIYSZAecPcLBtx/NoC97v5FADCzdwJ4KoAvmNk57n63mZ0D4J5hFVXTvqKiYp3B0diQw/6GYgHAd5nZqe1ChWcB+DSA6wFc2qa5FMC7hhVUNdIWhnIwhzU+1tRK2oRqgfuTNDxTqxnN5tAsgNPb432SNgJFqulx20paMD+XtmceXS0i00AivZq6UW6Wfwu6Ws8ium4HDerEcfYuwqyP8z3o709ua7SfNUpNF9iCPHDCzxj5+Z3tT/JA8pTcDlFv1FNyuagLI3MBcD2zVMYyXYs8rMWq9q/H8byZhs1lsYtsunA8HpS7cJSluN9iZu8AcDsayftPAK5G09RrzexyNEPhhcPKqoKUwAM0BvUgYammUQygLZSWP0r+cHRQ8QeiPsrMfAO6/lZ2DehHnPne1B0wJcdZ1D6L0rPAZ6FSiqAz+INn4VMShCyINR2brHxPXTPcFvZR8uTD7ILSErJ4F9Hni4Xj/WiEPJD7qnliYH8pTwrs+1TBywIzm4iUdbAk5zxOAyV/bfQl32OBXxLKHDuYHMYXpADg7q8G8Gq5/AAa7XRkVEFaUVGxzjAZjXSSqIKUkEWa2RwuaaQL6NcCzm+P59GbiTmyubm9F+nZnC+Zs9yG7H4WdFKNQjWITItg7ZQ1atUGZ+memplZ21lrDg00i/pzQIW1MGU1aPAmCypl5nSk5+gyg4NoHOXWCLpG4zMtnuvcjK6bIZ4jys4i3az1Rh9lmqdq0Vp+HOt4yFgUbL5He4CuZRV5AmqdZWNhP7o85vFQBemah9JcAur71A+UTcjb6TqbNjH4vo6ywFPTKnxa+wtpou7s4+ePQl0DQLcMNsczQZCZjFNyruBJaUqO+dm4z/k51c1RcrUwuF38Xrh+oDtJsPBm4cVCic1f9WEyIZ2vT6H3XPvo/hKAOSmH+yIwjZ5rICPaZ5NkyZ0Rz59NHnPoTexcrvrYIedZWUCvL9j9UEp7dKiCtKKiomIMRNR+7aAKUkKJL5lpN5w+ZvnMNOKlj3uQa12hTXJ0mzUnNrOz4NC0/C8FSDRQpFoYm2OhfS2gGwRjDUyDUlw2B1uyIFiYyWwuxjPvQI+7qOYnt0WfTTX0qD+LHDObQu8pUZ8j+My13U116vJXNodLEfRd6AY4OfAFOo40c2gsHe4DdkHx++M+W5Y06gpC+yxshbGbi9uiPFJ+ZnaBqDsi0k8G1Ue6psECi03QLDLNgzii1GzqB/gD00gmp59B40MDGhOLTWYVQJGfy2Vy93apv+RvVddEyUcc4H5RU17NThZqWaQ96sgEIxPAdVJgwc91KuuB+4nfK3/sy+hOGGzCZ35hbv8u9I8HLpt9oVwHCx6diErmeJS1S+5z+/m6ujx08i4pCTyRszmu/s+sbUyFK01+JVfM0aEK0oqKiooxUQXpmgRvAabk8phdNTrPGhBrOzyLs8nFpmTGQw1NjAMXupRUCeFxHXLOWlgWRY9ysqAaaxFcp7oblC+aRd1Vu1eTnN0T/GxZ4EZdE6r5BLgtmoY1fV4HP4OuaQw65rGggZZMC1e+pO6Yxe3KlnKWgjIxrrjPOE+JY5tZSoGM0K+WivKDZ5N0y+i6abKA1CKAMwe0ZXSsPdN+1ZaImtmPm9kn2r+PmNmT6N7PtxurfsrMfoGubzWzD5jZu8xsur32MDN7u5ntMbNbzGyO0l/absZ6p5ldStdv4nQZeHux/ejRa7INJ4DewI2oZpiKmXmXmf36sYTw4xU/QeyPcjWaG4gPMvLEefgg1TfGbZiiv0jHbQmmgfrZQkBHXn5+XQHFgi3KDdOc7+1v/05GvzCM9k1R3ml0Tcpl+uPn2tH+Rd/G88UWcfqn1DB+RvVFRp/PSbpFdIVTtI+PQyhGfywnf1w/pDzGkqQLcJ9lFLGsz6I8bhsL1AX0XEnRHzx+ttAfv//JmfYTWyI6MaymRroXwL9193vN7PvRLMX6TjPbCeBlAJ6C5t38rZn9b3e/E8ArAPwcgMcCeCmA/wXgcgD3uvt2M7sEwOsAvNjMZtCsULgATU/fZmbXu/u9ozTuMPJZVH1XmX9JlyvqUjvWggI8qJfkP2sdXE/cg9QXebP8Iayi/h3t8W50gz3cPv34uS/iOIQ158mWBrIWyJzQ7MOP433o9mU8y370tCH2UQP9gbhsQuM6ptH4ojV4FmmysvgZ49nVygD6d3wKZBppQMcVWwrcrign7rHlk72nzEedaaGz6C1F3oXu2OYVe0pT4/hBPPMCer5b9qOeyD7SVdNI3f0jJNRuBvDo9vhbANzs7l9190MAPgjgh9t7G9HbXDA2hL0YzWarAPAOAM9qNxz4PgA3uPtiW88NAC5q0y2itz9wRUXFuseREf5WD8fLR3o5gPe2x7sA/JaZPQrA1wD8AIBb23t/COAtAO4D8JL22rkAPg8A7n7IzO4D8Ci+3uKu9hrc/QWjNIpnTDVlgW4EXP2barZl5Hb2KS6iG83kc1C6QZH+TIOO66F5zqPLLthN6VRbzrRIBmuUocHx83A7M+1Y/cJArmGrtp71n7oQ1P+Z+WtZ6wxtjDXcgGpw7IcMRsQe9Gt7THliH3M2rsKvnr3PKEPr1/dfen6g2xfZs3CaaGemObLWyua+Yg/y8avvZTJYez7SVRekZvYMNIL0ewDA3T9tZq9Do0EuAfg4WgeHu38OwIVaRFJs6eeERv5BKkP3RfMGFDyoM4pKmDIc1IgPjk0/zsMmawiOzCfGHw+vjGLqyxS6NKtF9ATmDvRW1LDgLLUp8pcoMpqPhSQLWXZXlExmXTXD1Co2H1X4gc4zgclgUzTyAD0hlglmfk7tM50gM47qkhyX+ngQ3zcTOjH+uJ+4/qzPSoGfOM/cLhpEzNwMgWGuisn6RhlrS5AeU9PezF5uZne0f1vN7IkA3gDgYnc/GOnc/Y3ufr67X4jmHd45oNi7ADymLf8kNMt3F/l6i0ejuypP22Zm9pX4W1uvpaLixMTevXvB313rllshJvNTI5PEMdVI3f0qAFcBgJnNAngngJ9wd7YwYWZnufs9bZoXAPjuAcVej2az1Y8C+FEAH3B3N7P3Afht+unU56L5EatS2xy935vDRjNnE4Q1lQDPrmpK6fpwDuKArrNGEmblEvrXfbPJm7VLt5RTrSfOOXCgGiUHErJ13JGHtWDWtICu5p3Rd7RcDZyxCc9lZVFq7ZcS/UrdKaA8AaZkBUoBpozKFddZ81pI0sU50K+dBkMC6NdAMzM9e4/stuD+z7R4payxhsr0L16ckGnAmTuI26LjjNv15G3bcPDgwdMwNlY3Kj8Mq2na/zoaX+YftZPQIfoZgL9qfaQPAnj5kEj7GwG8xcz2oBkTlwCAuy+a2WsBfKxN9xp3H0Sh68CQm9M6UOKjmEVv4CzSfaB/gGUfOH94nCbKycxhjQZrBHyuPY+IfKTJfJ/hDvguypO5NtT81g88c11wPexTY2EXQoTbyfXvSfJEeZyehUSWTts8Xbiuvr+SCycmvxD22STDk59GrTWyzm3gPs+EelZPtqxZ3Uzq19fJkJ8P6P5USCbg+Zm57Gz5M+h69ZFOAO7+0wB+unDve1dQztdR2LHa3d8E4E1H1cCKiop1hIeoIF3rcPSbXRmY06gR9cwcZPOJNVrWOoJUngVoOKjEroFldLUWjcjzc5TaBfS2+2OzGei2mbUr1YIzbVs1X9Axa3r6zAFeJaYBEdXAsgCZuj0CQRYH+vcS4ABRlK1Qs5i1Ta5zEF+Uoe+Z+1y10Piv2jXXX0rHwSJ2D2VmOdBYA1lfxDOXtGUOVvE7Y+14MngIa6RrHY7uR761PWaTlz/iWTpWoaA+uhhUC+h+wDN0zGb+DvSWiy6j68+LepVio2ZWZqaq8Gf6DftsM/dBHEdbZtEV5vw8ujig5NvTPsg+NK6D/YAhhMIFsJ3yqFDksrifSwKXBZweB8K/XHItZAJGJ1KlhmVCWd0X/J51sUXJtcHt4L5ZRv/EBnRdEzx5cd1Rf7b8VZkS2fczPsYXpGb2BABvp0uPReOCfHN7fQ4Ng/BFwxb2rBohv6KiomJyGD9q7+7/193Pc/fzADwZwFcBXAfgCgA3uvvjAdzYng9E1UhbbERXO5hvj9U0Do2GZ/AwsThdKfCA5HpoE1HePLqIsnRvTN0MJFsuyBrEdnrG0G64Dcrx1PrZNI/nzYJKQP5LlYvoaY6RP9No9Fki/2b0tpKL9FFPaKbRzky7Uo2R31lmYUSeTDvV/t2Jrmsle//KMlBEWWwO6yKEuBbtybR4bvMgru00epbXAeTgsRDjPMrjvWaVwaB83UgzGRyTjZ2fBeCz7v45M7sYwNPb69cAuAnAKwdlroK0xRH0fyRA/644ccwDHGgGCa+SYSEX5fKgYrMs8vDgy+g5XD8LsRAqJb8a+1jVdzcj50Bv0+K4ziabmtncN6B72SIGdg1EfSFY1YRmgZ/tRBVmcrZvKn/I/PzZuvgACy+m/KgfmP2L7GpgF5Dm4ToYujdBibIU5W5Hv8DKfLzcTvXjqy97F92boTzRtp3okrFZSeAxq+M/wJNVNnkcHY6Jj/QSAG9rj89297sBwN3vNrOzhmWupn1FRcU6xEim/WYzu5X+/kNWkplNAfhBAH95tK2pGmmLTSgv8cu0QCaNh2a2mKRjrWUB3VlZzSxmBGTplOitnMSA8jOZdK9pWNsMsAai5rua+RkhXs3JLJocZS0k6XShQKb1htaVcVT1nZUCL9yXGlQpaaH8ztVtwGMmy6P9t4wej/d2KTuLuu9B1xLSYFlmOuu7YI2yZJGwFh8aa4DfOY/t0vhXK2hyGEkjPUBc9UH4fgC3u/sX2vMvmNk5rTZ6DoB7hhVQBWkLNVN4wDOyqHusBMmEikZZ2Q/GH6HSjLhdoDwlygpHhLls9V0q5Srb+o7bwEK6FFnX/CqkOH/45PZhtOj4EpWl0Wx1G2RE9QV0P+TMfObyuL1xHCarTrBqJmcTLpfH5nOwLv6uPWcBx8+ZUeaiPSzk+dm4TnYnxRjlcZ69Zy5HXQ5T6I6zjCnCVCwuT10bR4+Jm/Y/hp5ZD/RWT17Z/n/XsAKqIG1h6HeqA+WPha8H9SkG2GZ0fYwBHlyz6BeqJSHLyAZ7pM8E5hQaPxfQaBf6EWjwKJ6H68s+Fg02xbWok4NF/OFqv2RaCgsFbotOFiW+LpAv0eU+0no5kKfvltuv/Z9prixI1K++IGmyQJCOMy6Xn4uDYqUxx5NK0LV4nKhVpde5L0LAs5DPJh/uJ/U3TwaTCzaZ2akAngPgZ+jylQCuNbPL0XR5ugCIUQVpRUXFOsRkNFJ3/yqapet87SCaKP7IqIKUwFoUz9qZH0gJ+Jwnm/HjeuZjjbyshbB2ou2ApA+tic12rifznW1BvxYb5bEvV0nfGjHO6DyldsY5p+f+zLRjNblZu1FqVQbtP2UwDPJfxnV+Ns4PdBduMLKoOWvU4UpiVwEoj1Lboky9F22bR1nj0z5Q2l7Uz6Y5uwO4L3RznVJ9nIddHpPDxKP2Y6EKUkIpWMICJtuYQVcMQdJtoWvskFehpNSkAA/YbJVJHA8LsMyg95PPwRvMXBjZCqJIqwGrLNiiO2FxWaV+Yh/1TnSDHNlEFu3NXAAa1Mt4l9HGzJfH/Xc+estoge5ktYQe51f94pmw4jp4fEWeAI8npX9pUIvzZdfZl6+TmvqJ2QXE7dW2lfz/jMzfXRK8K0ddIlpRUVExAVRBuibBhNpZ9MynQTPqslwvUWb4OHPoB12GtTM1ASN9RldiLQnoahqsdW1FzwTVqHdp1YlGgwOL6N+oY0nuR/4sGh3PktGXbkZXU8sWRESdrPnH8SyVq26WRUqzH70I+Mno0b44Uq0mu7oSSuYqp8uCYNvRfbelcjKNuUSHC7BFMYXeqq/QotUlFWXuonTRrhKrQ9u0H/mvQpTG7HioGumaxRF0fVElv2hGccnoLue3x2wWQtIolacUwY2PfT+6AoqpNHE/2hnpFik/C4UZycN1TiF3YegEoSZoNmGwsAB65uMB9G9GHGAfHZvJUU+0Rf3Y/Cwl0165ulmknsueQXeC08kyzufQM/NL9KeMSpT5wDnSDfQL2ZX4HHmyiHfJz8B9xn3OE2SJOcJuB40ZZHGCyUXtgYfyxs4VFRUVE0DVSNcsNqBrGjPYtGTzimdanYFvbo8Hmb+qNTFHszTTlyLYpUi/aresnanZnu2VqpqqBsiyAFc8Q6Tj1TPzcp1dKCVOq2ru0S4OdjCYtL6I3DzN+KQzSbqoN/6r5hzl7KL8GjjM1rCrmatmNq94Uh6vRvijbaxFK3Mi0ii7JPopFpVouera0qBmye0CSnNsUAXpmsRhdM08RkYFYSEUgo8HdRa11kHFm2mo8GAhyST0EgFdSfSlCP5ycqztjHrj+ha6vij/M/oRt5mFrzIVdMLKSOAamYZc5z7jj3qBrnP5uvSUqTylDVCyNobLImtztE/zaDRdy+f+CJfQLLo7W2nkPcB1b0b+nsOnOizqrkIxY7NEfl5gkC1Fzqhv46NqpBUVFRUTQBWkaxZMcOfAE8/CbCbOyP/MhNxK1zlwE2VzWZkJryY7m3/KwVQTMO5NJXlCs+Fny3iUqqkpjzXTjtTEz3iMmXbOJie7GQLqGgGlm0ZPcyv1MderPFTVPIcFS9RMLml03H/KzGBzfg+67oho54Lk2Y5uFD4LkM1L29kiGMSuCGhAirm67Foo8Y0zJoHWMR6qRrpmwb8iqquEso0hsqgsD7gY7BmhWvPr6pmS+cW0rIxWwh8iC9zMjMwisSx8I7rOkX71BXPdJZqY+lU5DZu6avZmPkzupxAKcb4HXSE9neRR/6D6ZTNGhvqlNQJfcqFkFCN1nWjbSu4c7gM28xUscDN/f1yPtm1H11XDfcZUJnVHMDL6HAvYbPKbDGrUvqKiomIMVI10zcLR1Sh0yRvQr1GoRsizc2amAV0zV4M9TLTOTKZlNLvXA/2/Q8+BE6Cr6ag5m7WZtU1NN8jMLi1/1EAc0NVuQuvL+kOfP6AaufJas6Ccmp8alMu0TX4W5ldqu7hvoq6oP2MgDLICoq4om/uVecysUaoWG+m2o8xDPZ3yMw+4NJaywF/WZ/zO+ZlKlsr4qIJ0TWIDcmqTUnx0wMV/zlPaQIO3UeMBGh9rZiZx1HwZ3a3SNJoa56WNJbhdau7zBwIAX6c2s2tDhVVGk+EPlCPlXEcIuMxsVF9yFjVWdwi7U3QS5AlKXSN8j9OzUORn5jYxU4BdC5xWGQ9aBgvF0oIGfeZMYC8jn7zZtTCFZvVWxi5QX3jmGol7OglrW/geT5CLAM7FJFA10jUN1ogyTmZpz8qMIhMf0AK6NKdsgGs5rF2wRqWrjEBpWPNjP5p+IPpcmebItBZ+FtaU+MOPdBqU4HtcJx+ztpOt4FKNWj/QKJu1ca5TNSVtC//PhIpOJDrJMX2I69F9U7M03Ofc5ggeavvjednPmwW1WHhzG1UoqkDNJr+Mx5tNvqVgXablTwZVkFZUVFSMgWPyK6JjoQrSFocBnNIeM8WDoaaU0pLYx8gmZOavU22GfVq6DR0oHWtqGtmeo/ZzWZlPK8rK/FYa3c/cBBEBLmkcrN2X8quGwlola96snakrYGtSHrsAeBvDOI/8rNGpBs3maOTfiv5NTErR6YXkOj/LLLqap2qr6suNsvidcTv5PQ8i1CvtLHt/JVqbap9H45qYHKpGuiZhdMx+HQ5wKFcQdJ0pJyxk1RyMATYH4H5Ko74npf1EuVwOm3Jb0Nu9h01LDchwmSyYmZ+obgf2vamZxz5GricTZBqc0GBJFhSbQo+Ko+b7NLo7NmU+YF7ltJPSKBVoWeqJZ+ZJ9Xb0C7slOY+ykBzzewnBl9XDz68TTmlzEA6EAd0xV6Iv6bjM3ClRDx9n1LSSj5aF8uRM+7XnI91wvBuwUpjZyWb2j2b2cTP7lJn9Rnt9xsxuMLM72/9ntNefbmZ/elwbXVFRMWEcGeFvOMzskWb2DjP7jJl92sy+uyRLBmE9aqQPAHimuy+Z2SYAHzaz9wJ4AYAb3f1KM7sCwBUAXjlqoUzI51l0UDSetTY2YYF89mVNYze62ghrJ/vR/QmLjDWgphzXqVpwaTMWNi0X0dUoNCgS5WpbSs8ZeXagZw6ruckbZWggjFkTrLVq+7P3xP3MVsSupD62PDITXDUwPlaNLPJoEIhdO5vb43DfsOadLdzYgd5KpTDZM+J91BvXOUComir3bWAKPWrdPLrfAkOZK2xtMDujVP/kMDGN9PUA/tbdf7T9fftTAbwKK5Ql606QujtTPje1fw7gYgBPb69fA+AmNA+/DOC+YeXya2GaDA8KXSXCx7p7TtyfQ3fHI77PdexEb6OKafTMdP0Q2WTK2ABRDwu8TBCFEMp4pKXBz1Fm9ZepYI+23S9tmZH0WX5uv1KhGBppz/zHLHhYQEbekjle8vFl/sms/oyHPIV++hq/T10+CnSXe25HswF1jA0ep0o9KgkvHQMlJkU2QYbLhifpjF2grqlsGep4mIxpb2anA7gQwGUA4O7LAJbNrCRLilh3ghQAzGwjgNvQjK2r3P0WMzvb3e8GAHe/28zOao8/AuAjx6+1FRUVE4U7cGgiUfvHAvgigD8xsyehkSk/DyCVJYOwLgWpux8GcJ6ZPRLAdWa2c6VlmBn/lD1AJ6XovAY0QNcjrWIXyvxQDojMSzlsGmWzO5P7I9CQcSI1UKABhYwTqHuoZgyELIAR+bltpQhuaC0c3We+Iz9nYD/6NT3WFjUow/m4Tm4LLyjIVuSwdqqEfn6H3BagO04yrVufv6QFapt0bGrwKtLpAoOoUzXnTFvkd5axPuJj2y1lg9JlLpAlAHv37oWZfYWra63MleHISBrpZjO7lc6vdver6fwkND9m8XOtMvZ6NGb8irEuBWnA3b9kZjcBuAjAF8zsnHYGOQfAPUPyOoDT4nyD2TdeJgsS/th05yeNjGbCSwVmtuKDy4j8mYDT/DzYZ+ScBVT2sWkZTIjnHYd0lZK6GTLiPbsAlBa0R65nghTIJwUlmmt+NttZ0Gf0m0i/h84zfyu3SVc5bUWPNcA+Zt2xKfPXZn0X44vfS7ZcOGubLh3NoubLci/qjf/8bku+4CV0f9uJ3UEBZZRwfdu2bcPBgwdPwzhwAIdHEqQH3P2CAffvAnCXu9/Snr8DjSBdkSwB1mfU/sxWE4WZnQLg2QA+A+B6AJe2yS4F8K7j0sCKiopjDG800mF/w0px3w/g82b2hPbSswD8M45ClqxHjfQcANe0ftINAK5193eb2UcBXGtml6NRCl64kkI3Id83cwldMzXTblhTAbqzO6+vZpOJTcnMXC1xBVkr5oipkp2ze2rKqqkdz1niKjKi/Rm3VjWiLCABShvlBXSjkSzYF3Xw87BrITPNWTsMrZm1+AAHu1jrZdfILLrRbdYomXsMdMeMauqs0WWBPx5Xp6O7WY0GjqLszcgDnNFGHhvcNtaoA6rBZu4RbbNaUVl7x4ZPLGr/cwD+rI3Y/wuAf49WrqxElqw7QerunwDwHcn1g2hmlKPCg+iag1k0VIWSbkjBH3VGZdHBNi3p1YSKPGxa8yBmwcfl8UesvjcWamqacXtKwospPixEmdzPdZb2B4jJh9vD+xOwgM9YB9GvavZqOqAr1EtR6y3oUs74OvcRP1dJkKi/ksdQPG+YvJlQYh8zC/JghnA6ftdRtjIgdPVUiXgfx7PI3SQxxjIXQtQV6XjCYJfFROA+qo90hKL8DgCZ+b8iWbLuBOmxAq9s0oDSjvb/AeQfDtDVvHiAZUv94joLERa4SoUZJSDBQkWDGIFsZUsmfJgyxc/APto4jj7Yg+6zZh9b1mf8zDx5ZJqqBvvUL8ht4+eKY95AJp6dfcH76F62SkeFCE8mS+gFYXahfwJGcm1Z6sl4qJnwYcHK1zLhpe+MtX1+HuZBqxbLKPmidflyQANsE4FjUlH7iaEK0oqKinWGyWmkk0IVpC0c3XXnPNuHpqLmuGoamS8u00wiv2pNrB3wrJ9F/RWsObA5CuQbK0fbsqi7mtwZ5Ss0mC10Pkh7imONppd8aSV/W8l8V5M1YzqwKRymta7fj3Rcf+Y+iHTsAthN1zOKUUYFK1kOw9gIAdb25trjeSmL07K2y32+Q+rJXEuLaPqLtVDu54COqcw6GwujR+1XDVWQtjAAsaCWfXQlwZM5+pUaFMiEGl/Xj4MFhtaZBau0fvXdnV4ol/266kJgyhQHJ6L+oPhk/cGCQCccFmoKzsOTU0Zl4vSaRwN5g+hnmatGfeSZa0UDKNx/ynctPasKxoyPvEjt2oFmUud3wH0TK+O2o/teSu4RDortljSZsF1GPhnwMwWyVVqLAM7EJFA10oqKiorxUQXp2kUWwVTzlzUR1U6yAIWm09k90qjWyhoRaxfZrK9R7dBcAl+XdvLx5vY/p49niDK5fNYUNSLNeWcpT7ZHQWiErKFmgSPVFFV7Z80nYw0AeeDsdDQmcEZC53Ql90No7SXNmzGMYsXX4n+kZY16H/qDmlnUXzVtLZ/7rLQIIXPTMBMh0mWBSG1jieZ31HAHDtdg05qEo/shxDGbv3EP6PdxslDRQc2Djakoaopnpr76qDJOplKZOGrM/FBuf5j583RPzdXsOfUDydpcolUpM2HQB5v53vhcP+rMfwl0XROcZh5dQcbvTF0T2fV4xrn2fDe6u3dxuZmAjuvqw4w6M5fLNPrbnE0y/C500xcW+Oej9/7Zx63jTyc4nlhZiWCBmVGeJmfao2qkFRUVFWNhgjzSSaEK0ha8VpZNHqA7I5e0FjZnGGzaqCkViFk/22hCNS1uo66kyVwIrIGA0kQ0OGMaMLajq52r1pRF/VkL58CLal8ljqkS/7MIvN5jjZI5keqa0XeUvU91DQSUAL+E7rZ4qolyPs7D17k/WfPcn6RRN4laDtz+bK+D0CCjD3cj13x1BZ0yGLjN2XhmK077ZWKognRt4giAjXTOJjjvG5oJ1Sm5p77MjBwOSaPuBL5eOo76lByvAo6j7krFyT4KjbTyc6lQy/LzM/DHr0s/B0XgA7oENRNwUQ9/8Lwyi58r83dn91RIRRrtIxbEXBY/f2ZyR9tmkns8qeoqpex/huw5ta6SIOTFFdyWpSRtaZyzX3ha7o8NRxWkFRUVFePBK490PWAJvR9CU9O0ZHKrRpJxAkF5MvMxCxzwumsNCCjfcDY5VhI+t1HNrpJ2pAGSaMsi+rW67FmyZ86WWA4zrbmNmoe1QybED9pej+vn9rN2y8E+5Yeej+4PDoLylCLYJe2QtWgeP5kro7TIA8l1DXyqBqnBuzjO+gLoavglXvVWdIOYbMWcW2jvilCj9msXG9EdPBzpjoGgm2nwAFFzPKN8LEm6El2GB3+JEM9mYXzgbLIxlKbFdbDw2doe75J2c/tVAGSTifqYOa22sbSmv7QxcSmar+wI7gN+F9wHs3SszIbMncHHy2gI8FGGTgQZuyNz32QThvpB1QXEq7MyepO6WbjP+Z6O02z3qjn0qHExqZcYGdH/4TfW9nN/j41q2ldUVFSMgRq1X7vYgK7mF1AthLUONXMZWeBFieYlLVY1yixYwW0JDYI1tyxSrdoVl8kugJLJp+B6ONKu5Ue5O9Az+SIglmnlXL+a47qmuxTg4mBTppGHWyO0KNbIWCMuLdedQ7MbWLaUVE34bCcmfUfZbmBxzNF0/s/PxGOrtKCDNXqgP8CYLRw4gK4WzVrwMrpuJ25zphGX3BBHhSpI1yYOYXg0Un1qbGby4NFIceZj5Drio+YPBEk6zquCV1f9ZCYXH2fmvxKvuW4gX/3DfZOZmfxBzSfPkK2gUhOchVqJVqSR+sx3p3UvoytwdBNloJ+6E88zj7JgYAHFUfdoZ/znia3kc9d3ukD3lNqWCVU15Tk/5+GxoWNr0KY33H4uM5tggErIf0ggG3wlWhB/IMG7y7QQpYIEZtG/MbO2A+hqROy70vS6IQd/PJlfVzVNFvj8UZeCTapRaYAl+0Bn0N28Gcg1PxUwGa0sBH+UrYG4Ul9MS/5MGGqd/D5YQPDkw5qyWij8zBo8ygJ8zONkQaRCWd9ntpqINV39z5MH+z5Vc1W/PmukaklArnH+iWmkXqP2FRUVFeNhcj/HPDFUQdriJOTUIvVXMa2GzR/VVjMSckm7yDQFhvqook7QddZquZ2sxbJGpyYXl6npuB7WRrgNs+j2U7ZiRjVY7gNup+5vkDEEgijOmpum1zyqUasrIitLtUs1uVXDjfZnftkZua7Uoki7gPLPcSudCnIP6O6nq21XV09GzeL3n7l2WAtWF1GUm2n02t6xMKHfbDKzeQBfBnAYwCF3v8DMZgC8HY0rfB7Ai9z93kHlVEHa4hC6VCD1pwVKqzvUR8qDPwb1dnRNcx6QJSGlJvgCXefggppz3N7Mp6a+0pLwVE5lyazTMqOdvJmFClIWLCxUl5K0QL5KKQuWKH2pRL9SHmvm2uDn4uCKumO4fu4zXq6qFC01rXU8RB6+ruY8T7IZ5YonAlYUgO4z6sScBdhYiQC6CkPJL5z5hyeCyfpIn+HuB+j8CgA3uvuVZnZFe/7KQQVsGHSzoqKiYs3BJ/NzzANwMYBr2uNrAPzQsAxVI23hcs6aG2sxWWRa80TaQGghGmzgtLOFepbRW2XFK5aUfsXngyhCupIo08JZI9RgFWt6qpWz5laiDymth+vhtrGmUwr8KNNgOUmnLgjVoLIgn2pnrGmxts9jo2Rycx/re1DLJWNXqGtEy462crum0XWTLEqaLNKu779EWSptaKKUqyxYqRTBsTA5jdQBvN/MHMAfu/vVAM5297sBwN3vNrOzhhUykiA1s6cB+G8AvqnNY00d/tijbPyag6FrQqoJBvR/eDoQme/Hg18jzVoWm3JxzgORTUO+rmwApu+oGRbX2LxTRgF/1BlXUfe21A8mO+ZnUdMc6PYtuy1YKJUiyJyfy9T8pXem5WT+YxaQ+iw6mWQbjZTGzyL668zcERk7gwViyZ0Uk6/uCgbkEwu3U11D6nLKfNDZvqdR1kQFKLCSYNNmM7uVzq9uBSXjae6+rxWWN5jZZ46mSaNqpG8E8IsAbkPjlK2oqKg4Phj9x+8OuPsFA4ty39f+v8fMrgPwFABfMLNzWm30HAD3DKtoVEF6n7u/d8S06xJHkGs+rJHMINdU1aQuRZ1VO8q0sihvSa4BXU2LZ/oI6JRMq6xODoIA3Wcv8Ss5OKQamZKwS+acakccCCmR7UFpNSBU0iIz6wLINT2gn0GQBcgUXH+0Ia6zpsdlqXmduZC4Xn0uZkewCa/acsY0iPYyIyDTiHegu8+qtjXcS4vo7c8wL23eQddL/NKjx2SWiJrZaQA2uPuX2+PnAngNgOsBXArgyvb/u4aVNaog/Xsz+x8A3gnggbjo7revsO1rGqVILw/QLIIdH07J/5ldV1NeI92Z8CsRzePDyRYOlHx0kT7zC7JrYjnJo+2OsjNBwKahUrT43g70lo+qUOU6Fin9PnTfR7bpiZYzS8eD3o26QKL9myUdKF3mAtBjhk7cmfDRd76Arv90Z3t8P5VVYmpEWRnNi4X8AeTjQsueatNGWfxu5+mYx8LEVjZNhv50NoDrzAxoZOGfu/vfmtnHAFxrZpej6fIXDitoVEH6ne1/VpMdwDNHbnJFRUXFJDChTUvc/V8APCm5fhDAs1ZS1kiC1N2fsZJCM5jZxQBei8aKPgTgF9z9w+29iwC8Hs1udm9w9yvb61sBvBUNYfbH3X3JzB4G4M0AngzgIIAXu/t8m/5SAL/WVvmb7n5Ne/0mAJdFugwb0J1Rs8ABa3C6Blv5fqxRlPiZrFlw2Rrs4vxZECXMNdbcOKLPmhKvZy+Z06XgAGtQ0d7MnNZ0HKxi7Y6fZxe6m3KwyZ9pavPoamuL6AVYltB95gC/izjP3CEMNfNZU+Vnm0XXVZIFZLJrrPmxdpctEQ1tkPt5nsrOuKNcb7iGOEAUYFeRjjN1+2TuBE7H7IYssDURrMe19mb2CACvBnBhe+mDAF7j7vetoK4bAVzv7m5mTwRwLYBvNrONAK4C8BwAdwH4mJld7+7/DOAVAH4OwGMBvBTA/wJwOYB73X27mV0C4HUAXtyuRng1Gq3ZAdzWljNwRQIjPtjd6H6IbD7y9ZJQ1cHDH3sgi9Rnkc6SUFL/HAsJJX4jOeaPAOjuO6nmG9fBx2wCa/Sc3Qzsl402sECPdJk7AOj66xjcTzNofiJD2zKFnvm7C/3vj902XK76cKP9zCzg91yitgH95nFc247ub0vxuw1kwjcTzFEG0P/LqTpBZ+4QzsPHOvnqhMlt4MUK2TMvAtiWXF8x1uAS0VEJ+W9CoxW+qP27H8CfrKQid19y96BrnoYedfMpAPa4+7+4+zKAv0BDiAUaDfVI+2ftNSbLvgPAs6xxcnwfgBvcfbEVnjcAuKhNt4jKNqioOHFw+Mjwv1XEqD7Sx7n7j9D5b5jZHSutzMx+GMB/B3AWgOe1l88F8HlKdhd6Ptk/BPAWAPcBeImmd/dDZnYfgEcVyjm3TfeCoW1Dz2HPpGOddUHXA6FNZEEdDWhMJcc607Npx+mUK8hQEyzTljmIpBrUPN0rEdW5L0Ib5KAO18/9lxHtI30pqBTHFwL4OyorSwN0A2FqimYR0Qg8ZWYqa8c70f3ZDFCabJlolM3PzNohm8zsCmAtfBrdBQ2q7WdmO9B1BywlaZblOgc1+XgQA4LHmb5zbjMj2+FqLKzjjZ2/ZmbfQz7NpwH42korc/fr0ETJLkTjL302eppmJ2mb/nPouRMCpfTFckZqG8ofFdNnMlNQTS2ODrO/iU1eFgp8P/6Hm+F+dE2rzLWg5fEHG/miHjbZVIjxh5wJOKbOxLOwYOMPOTPtWahonUoNivZ/CPnkwX0J9E9epUmN07Mgi3q1nfvofuYHZndMVhY/C7d5BxqXRclMjzy6kojL1j7IBCz3f4xrpullAlMFMfvbVXiyX7T0npYk7USwTgXpfwJwTesrNTTv4LJhmczs5QBe1p7+AJFfP2RmjzOzzWg0x8dQtkejO34Vkf4uMzsJwCPa9twF4OlSzk0D2mag8cIrm3SwqLM+oH4gFh7qV4r8wbvbhf4ByR8f+yu5Pl3iCfQCLZlGqlxBDTBkGtWUHCvNKtJmgRAuR8vS9HFf7+l1FtBbJA0L+RlKlwkILZf7nIUvf/yq8WWc1Ci3FNTLfJ+D3r9ylHlccD/vR1eQs4WTBXliLJTeU1Y/91/2zjk/v9Os/QCwd+9emNlXuHnk8hsdE9r9aVIYNWp/B4Anmdnp7fn9g3N8I99VaAJJMLPtZmZtsOl8NP17EMCXADzezLYB+FcAl6BnxmcIsuxHAfwogA+0Zb4PwG+b2RltuucC+NUBbXM0vloAwEnNWtuKiopjiG3btuHgwYOnDU85AGsw2DRQkJrZS939rWb2S3IdAODuv7uCun4EwE+a2YNo3AIvboXZITP7WQDvQxNcepO7f2pAOW8E8BYz24Nm8rykbcuimb0WwMfadK9x9xKTpw9H0J352QTPVgBlFBGOTvN2eXydKS6gY87PZtIWOl5AV+thbUR9nqydlNwBJZ9Vib6iWoealuw/ZS22FOXlKL66QCLPVnQ1Ok7P2jb7O3egp9GzdqlULG1L5s7YjtwKCQ2u5P/NwJpijAv2KzIDIbMo4v1Fuu3obk5S0ij1+bPoPufRfV45PW9OogwAfk62fDKNfCysQx9pzBwPT+6tSINz99ehoSpl994D4D0jlvN1FFYauPub0DAMVgxHdyAOMw0ZSpVRKtJUcj3bbSk+yvCvRp4siKA+URZS6pfM/JVxnnE3+WNlqK9YTdgA8y0jbVxnM5/T6B6W0eZ9VPbJ6PafmtqR/37kfQF0KT4ssPag20/ZRBL3+D9TqwLT0q6M4jQt95aprHl0faTqj2XBlo0tHaO6yisLyun4iTw8QWvgTYN9mWuMhSq7ZcbGehKk7v7H7eHfufs/8L024FRRUVGxunCsL0FK+AMA549wbd2CVzZppFe1rfivEedsD06GEt151t5D9eyRe2zmZRpYmGJxzkRvbrMSxtWdUNqQJaNSqUa5FT2a0Sy1hQMiGlBS1wC3i9Oc3B7fj16wDm17t1K6UgQ6tLs59Ez+WAmWrWwaRB/ScqM81tZK7pSsTO7z3ZSHx5IGhbL3qf2YRd2j/5mdklkeWaQf6G1MkrlqlKnB73+Yy2PlWGc/fmdm3w3gqQDOFD/p6Wj8mScM+LVo1JoHHvuR2PziATor+dlkjuss6Daj3/c0Q+miTv7AS/5JoCsI1UepUX/+EDQ6DEobaQLhl408u6VsbidPONnHFuVl3Nlp9PzKW6mOh6PZHPdD7fkB9ITq7ejRxw6g5yq5Hf19xu4cjehHm7lftM28Uiszs0u8yxhjLHA5Ha8SysxkLa8UjVcfqfJteZxwupnkurZnEHc3myAnJkjXoUYafX0Sun7S+9FEzCsqKipWGQ4cXkdRe3f/IIAPmtmftuT4ExZs2vPMyTMwaydscql5dDq6UePMNGYN5H70r2DhqL9qvnHMxGo1IVmLBV0H5WFe5h50nzsLsKnWq8iCOhqsYe2cy56S48AcupoSB3cedRLwTYd6eSLdDnTf2RbKP9ce70a/aQrKE+caadfNZFhb5AAfj5NZus7BQl2cUKqT3S/cVnUVsRXEJPr9lJ616Fl03+N2SqfsCG4/j4fMQuJxmmmmY2MdaqSBr7b7kX4bei4ruPszj0mrjgN40wEVGGEmMmWGB0sIsfjIdiOPlLOA4AG2gP6oObchPkQ2H9msznyfmWuAfYL8fNq2QW4DFhwsGHjCYX+rmnPaZjZh2cyMPczuRG8P0Oef2ivniZuB7duBh93cnH/0q8AT2nvvRnd/TmZTzLfHO9o0LCQYy8n1ZTn/LgA3F/IHtqP3wajveQG9ftqBno9U2Q1cvwrJ0vvjzUg06q8CPOopTY46WWRt4/p5zKmbaDJYe/SnUTct+TMAn0GzectvoBmPHxuUoaKiouKY4dj+iuiKMapG+ih3f6OZ/TyZ+x88lg1bbRxGd3bmGTY4ghpNLQVnlLuohGagP4ihPDye9bOAhGogqi3tQT9Y043ABWsbusBAoeWzacj3ZtG4N4D+5+K2nEznrK3tRG/3mcej3XkGwNwc8NzntvmngY1TGxBhwl0fAG5t0/F7Ohn96+WBXnAsI5RzOr3Oz8Lc0ZKmxkR7DegFcwBoNBO2YkoWAZCPoUV0x6au/Y86dZyw5aBjYxRw/2VaKGu/kws2rbOoPeHB9v/dZvY8NGPz0cemSccHGwB8S3u8C91BwYRmHmzqK8r8UtPoEpIzKtAW9JvgJdOu5DJYQu4C4A9KTcsoJ845D7MKmLS/QMfqL+P645nPR89kZVM0WA5M4Qrcj54gvhe9Lcx37gQeMfuI5uSCJwO778QFFzQi96lPBd785ubW3ywA4dD/OsqTD/czg/2CKmzVlM5cMLoSLvNTR1/OULrMncDjL9rLPnMul03zTHiHjzSjL+lqqEwoRv4SIyQTqkyL4n4ZG2vMtB9VkP5mu2HJL6Phj56O5ldFKyoqKlYXPtmofbu5/K0A/tXdn99uEv92NLHJeQAvGrZB/EiC1N3f3R7eB2Dsnx1ZizB0f7CLZ91MU2SzNiNDc57MvGGNNEjLWZ3K1WNNibU7bg+7AFhr0OCQ7hLE6UK70gCEasSsmWXk9t0o73jF+ReozpPRC/DtALCpHaVTUwB2PL452XIO8L0vwCkP3g0A+OxP/Ram28ZtQa/P/7e0XzmdmWkLdJ8/e5fhpmEtMuMRax/xO5uRdGyhLCXHoUWzCc/PxW3LeMyz6C720OAn86WzOsLNNCXncRxg90GJUzo2JquR/jyAT6NnCF0B4EZ3v9LMrmjPXzmogJEEqZmdiWY7vDnO4+4/tfI2r0048hUoKlRBx6WINNA12fgDU0EC9IRo5ldSZFFazgv002IyAbsk5xlxW9upppy6LQL70RNE7EvUVVr8gW2XZwjcDzS/8AXgcZdcADzvZ5qTT70PwMuBd58HoPGfvvmtza0z0ZsUt6LrWlAzfYbuZVFrNVk1TfQHv9tsww6gy+yISTlbTTYl6VgoRb9GedwepswxdFLIhLxS87hcfn4dQzwRZWa+rvKaCBzwI5PZrM3MHo1mk/nfAhCLji5Gb0vOa9Bsxzm+IEXzu87/B81m5SfkT3bEJquBONYVRxk/UQew8gAZSkWJNPvR/9FFOtWi4jprQyxMWbtlKF1FNWe+p6tjIr8+VyawdRlr9NNm9ASctu9kdH8PKX7zez+AS85rjj/95lu/4cfGU58G4DX42s0fBwAcONDdw/UAtZM1Yt3AI+5FuwKZMMz8nZm/kMeGvid+l9PIfwOK0zEyihJvLpIFmNh3uoz+Xa60PXHMaVRrzyw0DZ5yWZk/eVyMqJBuNrNb6fxqd79a0vw+gP+K7oKjs939bgBw97vN7KxhFY0qSE9194ESuaKiomI14A4cGU2dO+DuF5RumtnzAdzj7reZ2dPHadOogvTdZvYD7XZ3JyQOo0eqVdNK6SxAf9Sd/V+qebD5VjL/eXbn6HjUFemYSsVllXxnTItRig6b81oWX+dovGpamY9PNRp2DXAa7kP1y36S0p1xR3O8bx/wLTtvBwB8/p23Yn4eWGwf9qY7em2+D8Cm9vgAPct2dLXmJfQ0UX3P2bFaDXEt0rD/POtXtSjYNObnVwaHanQ8njYnbeN6OH3kybRlfq45dPuFNXK2PJQmxxYa93OJSna0cACHJmMXPw3AD5rZD6Axik43s7cC+IKZndNqo+cAuGdYQaMK0p8H8CozewANFcrQbDJ/+uBs6x/KiYsBolw9NZfYIc9pmLoSZuXt6A4ypR5x2cEPVVqOls1mPtfPPjEOUKnJxssquR42eTXYAkqnfjVts37sS+inVgHAd6P3UT9hGfjFX+/ZdFtPbQNQbfow7b8NjR8K6PpI96P7LnUDDu6zOOY+yvo843RG8DCuZ/SnEEIlv2LmzokgFguy29vj7ejubVqaCFSYcZ3xLg/QfXXzaN4tlI6DVeremShG10gHF+P+q2h/RaPVSH/F3V/aruK8FMCV7f93DStr1Kh9trFzRUVFxXHBMaaRXgngWjO7HM38lG4kzxgoSM3sm939M+1vLPXB3W/Prq9H8KYlQJlKlDn0gf5Zm7VDNp8i/zR6GlSY7KEd7EA3cJI5/tnkDpfBcpJHKSpKV8m0IL7OGtAU+gMNrFGxdstlZcGGabmm7IjQdO6kPDd/qReEegKAvV8FHvxqc76L8v8fqnMeXS1S3R+zdC8L0LGbIjOzuW8yE5bNfEia0qqvWXT3IAiE+c59zlowa65s1fD447Gl74mfM3NHxTFbGGyh8NhQdoI+yzhwByYUtKcy/Sa0P5bp7gfR2+5hJAzTSH8ZDe3p/8vqBnDCbFqyEc3PkQJdX56a1qA0GgEeFrWMfEB3sIUQZdOKBRwoD/vYlLeZRVP5A2eTL6gwTFPKOKWgPNm1LA/XqR9uICal7CNjM3cPev18K3pEv4NtmrML7S/5a9XlEOe6sicTUBkbIzPh+TgYGZFfmRPcZ4FBPOJSOq6T26wT4Tzd4xV4pTZzH6k7osSx5n7i40X0lvuOizW2sGnoNnova/+fkCT8ioqK9Yc1uIveUNP+BYPuu/s7J9uc44dDyINKrOXw7FoyAwOsObJGw2kzDp/m5zaUOKUaeeV0GsFnLXgOvWBFKbqsGhG7LBbRW4E0j66GHuki0BP5VYNiE1k1dLTXIjL9dfRrt/OUh5+TNb0sIBR5OMDFmhu/Z+XrxrGyICL/9sJzcbuXCtfjmMcMc121nVw/a8ulvmReKdfLGj0H+1TLVxdOtlcAjxlmekyKR7oGf415qGn/7wbccwAnjCAFuvSTzJ/Hg4I35gjfW2bm6kYdmfm5JOmUrsKEaI3Cx381k7ltg8y8rDwWHkC3XSrII8LLwodXNnH7mdYTq6IyygzQFVhIjqO+kgkOus6RZRUc2cSkq3+iLdzHUV4msDVqzwwIdTlkLhimv82gS0Vi5gEvvGB31Cy67oS4voAeWyHK44mltFxW3Q8q9CMvKx+chtsyEay97UiHmvb/frUaUlFRUTEqJkF/miRGoj+Z2W8D+B13/1J7fgaAX3b3XzuGbVtVcBBQzcEAaxQZpzDjy/FmHMvoamCsAXHwSpcBbqXrrNFwftUCsmCFbpKhwYJMo9JNM5YljUahox4OhGkfAv2mo66DBx3P0zlrzUqCzzR/1pSAnitiF7rvOVserO3XACNr4eq2yOpXVwA/D2MZ3d+418Anu1pKTAi2dvRZ2FrK3i2PJbWo4pk0D/cZ98Vm9N7fpIJNvt40UsL3u/ur4sTd721XA5wwgtTQI7urQIyBczq6Pj9dZ80mUCmCjEIaoDuouazdyXU2yTej+yuepXRsWsdEUDLhM99ZtAd0j8vLXBP8XPxRhhDMPl41HVlA6ATHbhduY4lytEvSZJOfThAZtH+5To6Gl9LEs2TPtoyeL5h9r+yfBvqfmYUam+HK/NBxx/eirFJknlczlahN7IflSbnUl0eDSdOfxsWognSjmT3M3R8AADM7BcDDjl2zKioqKnKsu6g94a0AbjSzP0HzHD+FZnupEwZH0DXzeNbeT/+zwE0EYZaTdGp+ZZqOXlPtJBABmkgbM/180p4s0j+N/v00A1NSdkD5ldmzRHmsqXAaNdXjP+fX/QUyQrcGizRYw/VnQTQk17JnY41OTWE1YbU9QFcLU5eNaqcZD1OvqyanWn3UXxpTGuVnF04cK4MgYwNo/aUAm1oqrDlPBA4cenB4stXESILU3X/HzD4B4NlorODXuvv7jmnLCjCzkwF8CI1GfBKAd7j7q0u7WrdraC9z98sGlouumcIDMfPx8QAF+v1iMcBORpeik314UQf7pWLwbUf/T1fEf/3AMnNsC13nyLjSd7ag+/EGppF/LOFTVT+tgj8kFmhBaeL2ZEJF6UfcrowiFsfsy8z8kRHZLpUd4Pz6/IpM4LNprml5YtQxx+3K2ARALjgjDwsyfv884eo7y6L53JfslgC6iy14NRW3k/t4Uqa9Y/2a9kCzg/Qhd/87MzvVzB7u7l8+Vg0bgAcAPNPdl8xsE4APm9l7AbwAK9zVmrEJow3kbDZeRv/gi3v7kA8gDVRo8IopQ5lGx1pz+NoyfuEgTUWvc55Mo+CPValAqpHyx5thn7SNNTfVVvnDV02Nl3ry8/NzBVj4Rnuz98nauU5+jJIvupSG26aaJgs/fpdKhQO6m2hnQpEnP31G9XEy5Y/HXwTl5qUtC1T2bmozvz/9Fkr9d9RYg8GmDaMkMrOXAXgHgD9uL50L4K+PUZsGwhvEO9/U/jmaXa3D3XANgB9qj5fR7KxWUVFxAsDR0J+G/a0mRtVIXw7gKQBuAQB3v3OUXaOPFdofq7oNjeV7lbvfYmbprtbu/hEAHxlW5mH01torxSQz4/SampB76B5fz0zzZTkv1anugBlJzxrOXPt/N7qahi4uUHNN6+cydW02a95qDsd1pr9wG3XzlJ3oukC4Xdx/SqtivzYfZ9ZF5i8MjUwXWEDS6fU45+fhsgL8LEq0V/9pRiUCXc/cCmwV8bFG9COt+tKZ2sSLKHajH0vouppUu+V3lrlWJha1X4Ma6aiC9AF3XzYzAICZnYQu9XJV4e6HAZxnZo8EcJ2Z7RySZSjYR6r+JvYdZRw+pZioCZl94LqMkQecBjTUhRB1cFmcfxr5TyAzLWVazjXYwIIxPgp2M4TwZ9OaP14O0GWBn2hXPM8B9Pd7IPsQ9V2VgmV8rubnMrpCQSk82n6dyBgaCMxcQ+qmgeTJhLIG7tRVoDQwyPk0uvuUal7u82y5rEKFZLasVtvF38KkeKQT2th5YhhVkH7QzF4F4BQzew6A/wzgb45ds0aDu3/JzG4CcBFWuKu1NbPCN969HdOWVlRUAMDevXthZl+hS9PuvmKlbL1qpFcAuBzNL0D8DID3AHjDsWrUILS/aPpgK0RPQcMkeB2A67GCXa3bl3danG8wc47aZ0GdOAcajSFm/XDUD6MsaRBHzc9M29TAR6YRhwbBWokuFkByn7UlbnMpQLWMfo0009Z0cxcuK3NzAP39FxrZ6eju2K5BqMzM1TYHdDMOfh9zaH+xFP19pjQxFNKxmcvpp5PjCHZpf8ZzsGuFn4U1f426Z0E51m41KFp6z+wmyRZjsAmfBc/UOojrWwCcuW0bDh48+I3v7miwgt9sWjWMSn86YmZ/DeCv3f2Lx7ZJQ3EOgGtaP+kGANe6+7vN7KNY4a7WDI66TaHry8qioeqT1Khp5uPU9GqKZr5MhvoEA+pjVR5kFsGNNiwOSacmPgsO5dVmHx9PENyXIcQ4f/RHLIUEmug+90sp6q0TDpfLwlvfUbR5XuoJRKQ62qVCJBPYKiAzYRfpuM6MdaGugFEnOb6ubpKMb6wuoIxmpSv41Eee5Yl7kf5MTAbriv7Umr+vBvCzaKxfM7PDAP7A3V+zCu3rg7t/AsB3JNdXvKt1RUXF+sQkTPuVctIHlTVMI/0FNL+092/cfW9b+WMB/E8z+0V3/70xnmNNwdENCmTRUNW0+FjN3EBpDboGofj+VnTNWa6/FJxgDZe1SDUTlY3APMIAB9WU+5dxSuOcTbsAazrZqiB2B7BpmQWopiQ9k8V1S76ZJI9qbfxugNwE5z7K1p1H/h3okdPZNOYAHD8jR7mBfvdBpgUrU0LH6SAtNq6ztlkazzy2INf5vwZYeRyf3h7vRj7+x8EETfuJcdKHCdKfBPAcd//Gd+3u/2JmLwXwfgAnjCDdhHJ0lk2ubIBpVLkUgddIfelD2i35mYAdUPNOzdzMnGbf2R70uxdKwjPAEexBK3z4ugoMvs9CouRy4HbpKh01gReTY/U3az+pzzLLv0DXGfw+b0dP4EZ5Wof6zuOZOD3QXWWWsTvUNx7IyPvcR+HyUPeMtk3Hr461zO0DdN9NlMcLJUrjaqWY1M8xt3GSEif96e31a9D8ltNAQTqMkL+JhSg14Ivo/XR4RUVFxeqh5ZEO+xsFZrbRzO5Aw/K5wd1vAdDhpAM4a1g5wzTSQZPIpCaYNYFMW9B7bJpxQCm0xhIPMgs8DIqgc52sHalppGZ6FojSoJQuo8zMQQ3QZKT3MCs1MAP0B1cyk1/rZXeC9l3GvY0yMr4m0NMO9wA4vz3eJc/CWihrrtwW1VRZ089cFVoWX+fnmkX/LzHEsQbx1LWQLV9V7ZzHbEkjZC245GbRca3aambFaRBTn28SGFFQbjazW+n8ane/mhNMipM+TJA+yczuT64bmv04Thg4ugOOTVgeyAH2nYUpk9FfogygG1lWlMxxvs5mvppYatpm5aofU4VP1jY1s+M4ViLxR6UreqJOJMdxzh8ftzmLzutHzGajlqtCIspFchzgVUb8LjIBoxH+jBGg9TN0oYIes2uA/ddKk8r2ZFAqmfrs2ZzP1sdzmzN2hLoasjYP8rGOixX4SA+4+wWjlXn0nHRgiGnv7hvd/fTk7+HuXk37ioqK44IjPvxvGMzszFYTBXHSP4MeJx0YgZMOrGz3pxMaG5HzKDnKyWYSa2ox62u0GmiilxHN1dkddH0GeYCKzeetdLwP/RpExiOcRlfT4nJB6dSEV+5g5Amt4wC6GjYHPrT+kvnOGg23SbVQjdSD6t7Rnu+me1sk/+72mMnsGrhhLZL7ohQ41Dyshak7gNNzH6k7JjOZOY0uIsiCf1G2MjqiLRwI4vGcBeLiv2rn2Tp+bSdfn7QPcIIbO0+Mk14FaYuN6PcnKtgnxR9bDBwe4IF9yKPWXBYLcD6PeqLcXeiazCqUMhYAR81ZcC2gESz8PFl0e0rK0jrZ7cF5MhK8tpExNUI6XqUTwiomKRbyGgFn36NSwzL6k5rpg1wmSi2KPJm/XIWV+jhLPmIttzROSu0M6PPre87ar+MckkdpU5Euc4dMTKBOaGPnSXLSqyBt8SCAM9pj5SdmmgZ/YLoBB/u1lB84iJOaDTTWwrSsbCVOtG2Y8M7awFpc5hfluqIs1uKQ5NEPN2t/tGWWjrmerNw4zwJE6hdmMC1rO3JByG1V7TCgfuks+AS6F21R6yAbTyp8Mksh2jOKv1I1ykzbLvmytV26fJTpW9nKsqxPx4WPaLqvJqograioWHdYr5uWnPDYgJwKw6RxoGv+KJk6iy7rLJxpLUpvUhcDa0Ss0bA2klGDFKp1qJlYirpmkfUSKTvSZ6tYBmnOvNepanOZbzcsgJJpyn3Lmi7Xzf5CfeaSps1QX3JAfZTRZl63vxO9n4SOPKwpZhZLmP8lF9SW5FpG5Ru0oikwiMrGWj33f/YskQ/Ix8TRYD3/+N0JD6Y/6QeVBQSA/gFeCkplPkYekOz3i3txjZcBMpTHyquWNEDEgTP1EWoZgVKASJHxDTUgUVqlw740FQ6Zj1ehy29LrpEsQKhgQTAKp3WZ/rQeduHwZs6g6yFEhwlvbnMI/2w8LKI7HjVAFXWXqFnqu+X7GviL8+3oClX2ly8m+RfR2zx9LKzX3Z8qKioq1grW88bOJzx40xKesdm012CNaiesxWxvj5korY57dfTzeaTTtcq8mYVG/bkc1vQyTSuuZ+a6mrCloBZr1YMI5iUGBNepGl0WRFK6EZvNQGMuA82698zNohQ1rTNz4Siti9+l0tmG9aVq/ar5Zs/Jz7UP3XeWuTwCmcsio+mB0nFQag+1RfNHXXuoLax5q3XGWuukUDXSdYCSj5BNpuwnRDKzUYVCQGk3SqUp0ZJ4AxM23/fQdRX4c+3/feiaj9xm/sD5o+A6tS3q71UWQ7SNV9+AjregKwD4g82oUCoU1P83T/VnEeTIF2UNWiIZ99hkhrSl5FdmX/Yy+jdWjvQ6+XId/P73ybNk5rQKyAXJw8/JikFAx3bmgtLloqDzBXQnrJKbaRLwdfybTRUVFRVrBpX+tEahv9mUmVystbDWE9oA32NkpHWgP0pfinqyppOZR6FBlMzM3VSO8iYzDquuUmEtVtuYBVjU5GXtSLXLEomczVxda851TMl1bhvQ5fdGeUBuCnMgjDWyzExW14JyN3nFF2uQrKlxnTxmSlF53QtB2QmcJ+Pk8nuNdFk/KQ9ULQV+NxnTJPIB3fe6iAn9+B2qRrpmwb8iCpQJ3fzhs7BVX1xmDimZmc3qbKUTUF5uqKZoSSjpM+jKJE6T+d60XO0P9s+W6s/oP+EvLgmsOM6ER6QpEdJ5WS1PXjy5LKCfchTQiSSjsoVZze0sUZl08gJ6EzELL/XNa1uCwcHPlm00swMNKyDKKbljGCx81Z2hfcnvWZfcxnFGi5sYqmlfUVFRMR58QktEJ4kqSAnZLKpaWGZKTtEf0B8pZy0kq2tK8rDmouuxs3ax1qNYRo9BkBHQWVvkANFSkl/NZNVOppNjdR/ocscs8KGBN1Ca0hJRdYdkATp1rWg7s3rY/NZnV64sv2fuv6ifo9zB5uC2ZRZN3OPnyrTwaAPQBKcyHikHpiI/P89Ccp3TRxCR94ZgHnLGVCkF1MZF9ZGuUWxEL7p9O3LfDw9c/nBC8HB0kknwOmDjvw5c/cgCmR8zM8tLNC2up7QySAnd/PFlPj41RYF+IQlKy+3g9mcTQOk5mSmhVB5+ZvXrga5vofTR1ix/oESoX0RvdVK0s+QaYX+pgieizDWiz8h5ZgBsbo/n0R2bGWJc8UKQ7Jm5DBbEkT/b8WsaPeZIaSJVAX20cFT6U0VFRcV4qD7StYtD6O5bybN+ph1upfSqMZaWy2k0mmdtNqd1Rs80WtWAS+k00q7aBN8LjYk1Yg6iaF9sRe/XTjnwwYwG1gL1GUvRcWUgMAFcf2YjM205L2tdrPVF/mzhgDIAMlN6Gs2zl8xx1h6zZ1SmB5Cb+UDXIlD3xT4qm9NxmzWyr/xnbm/k5zG7Ra5njBZtf2lxxaRQBekaxQY6Vh9nZj7uRpm+xGav0qeywRSCo0Sf4o8tiwaHsOMBW6LsMIOg5KoA8g9RV2mV6Fi64iZjI+gzaNQ3iy6zy4S3OYw8GXGfzXI2RSN95qNl4cntznzC2YKGKAPo+Q25rVEmuypm5H7mY4xzNs15YimNM33/LOSzTb+1n5YlfeYOyRgJWv+kIvl1iegaxhHkwR6g+4GCjnV1UokKUwp8qODIhKTyE7lO1nrVL8htzQSMBlA4j67Yyfyi0b4d7f959As3BfeLajMs/Fmj34Nc2Ef/ZVokCnmA/pVppUASa+SskbGwKgWb+J3xpjPcr6FdDqLTlZ6l5F/OhK8Gy9gq4OdXjZbr5glevwseGyxwOT9bOpNA9ZFWVFRUjIvqI127MPRmV15rzD4+jliqNqG0kIzcrNopaxBcHpOj1URiLYS1xszsjfrDvGT6TWgnTFnhctkc43JVo5ynNmcrctQHqAT0zBesJjRroYEg92caEWt07FqIe1Eum8DqC+V28nth7RTo9kf2bkvujwCPk4C6aXRlXaZtqr+4xPRgd4JaXplFxP0XY5mpVdk+CkD3vZZoVeNgEvQnM3sMgDejea1H0Pxc8+vNbAbA29EQeeYBvMjd7x1U1qoJUjO7GMBr0TT4EIBfcPcPt/feBOD5AO5x952UZyuAtwL4MoAfd/clM3sYmod/MoCDAF7s7vNt+ksB/Fqb/Tfd/Zr2+k0ALot0JWQfD5tpunSOB6v6/zIfEV/nDz+CWzz4WHjxQM6CIzGgM/qRtp8FQGn5pAqVbPmiCkIWpCxIuJ0s+ID+yYjBfcO0Gu5n3VyF25aZ7OqKUb9qiVOaBes02KU+6yxYxH7MqDvbt7VE/9IySpMP0HVhsGuATXvdACerlyebeOcsZKP/M7Nfy1VXxNFiBT/HPAyHAPyyu99uZg8HcJuZ3QDgMgA3uvuVZnYFgCsAvHJQQRsG3ZwwbgTwJHc/D8BPAXgD3ftTNL8nrXgFgJ9r0760vXY5gHvdfTuA3wPwOgBoZ5FXA/hOAE8B8GozO6OvxIqKinWPI0eG/w2Du9/t7re3x18G8Gk02wFcDOCaNtk1AH5oWFmrppG6O09Yp6HxGce9D5nZXJJtIxoN9gh6+4pcDOC/tcfvAPCHZmYAvg/ADe6+CADtzHIRgLehmawHzmGGnPLBGoFqGRmtBeiaQzyjs9bD6Xn9eeRnc5ZJ5JmZpBF01mJLVJ7QVLNAgq7MyfLHceYCiGeIdvLihEyjjnSZxqLBHjazVSPj/KzRZ5ZGlMfvJtOiVLtX8BhgmpFaC4p4r1lgji2cjNakFCo9ZmuHNdhwP/F5poVqUDGzVCJ/gF0G+yXPpKL1AQdwaMI+0lb+fAeAWwCc7e53A42wNbOzhuVfVR+pmf0wgP8O4CwAzxshyx8CeAuA+wC8pL12LoDPA4C7HzKz+wA8iq+3uKu9Bnd/wbCKjiD/nRnINR6g/LEwLYTTcWSUB9VKPmo1SSO9+rEyHxf79fhDUJNan1cnBqAbgeZnj3JKZj5/lNkHqvWoaVwyddnU5w8ZyCcsFmrR31yP+qzjmIVq+JTV98nvQ81v9p3yO1tE1xfO0fSSjzN7Jr7H1zLwZBR1xX8W2BnrITP/Mx8pfwva3klhRDm62cxupfOr3f1qTWRm0wD+Co278f5GL1sZVlWQuvt1AK4zswvR+EufPST95wBcKJezp/QB1ysqKk4gOEYWpAfc/YJBCcxsExoh+mfu/s728hfM7JxWGz0HwD3DKjqmgtTMXg7gZe3pD7j7PuAbpvzjzGyzux9YYbF3AXgMgLvM7CQ0v6e12F5/OqV7NICbBrTNIJM3aw4B5dqxmaomnwaP4lqmQbB2lpGiueyt7fE8XWcy9hK6675Zo2Luom5Gkm3yEe3IzGnV3kqaol7PgiOhkQ8K5GRlK2tAOY8KDtCppqz1lNwhGVNCoRYGPwtjrv2/D10WSLQpq1+Rac4lcj1bBKEFs7bNTBV2gXD9g6yAjJzP71nfz969e2FmX6HL0+6+YoVnEpZ9KwPeCODT7v67dOt6AJcCuLL9/65hZR1TQeruVwG4CgDMbLuZmbu7mZ2Ppq8PHkWx8ZAfBfCjAD7Qlvk+AL9NAabnAvjVAW1zNL5aAMBGM88+3iX0SOe70PWJqc+UzeRMeGgElAnsU+ju0sQk5vmkTi0rhChQ9rFmUWYugz+YiJSzmamkc10sgCQdP3+2AUhmmvI5P/OCtF9NWBYYWZ0sBKbQ/XkWrp+j2ewOUHOVBZEyGLI8U+gtqdU8XI+6NpT+FOD+A/LnZES7eDKJSZrrVzcB16nvnMdp5hdl3/EMgHO3bcPBgwe/8d0dDVagkQ7D0wD8BIBPmtkd7bVXoRGg15rZ5Wge7YXDClpN0/5HAPykmT0I4GtoaEsOAGb2NjTa5GYzuwvAq939jYVy3gjgLWa2B82YuAQA3H3RzF4L4GNtutdE4KmiouLEwiQEaUu/LDlEn7WSsuwotOoTEmbmSrYG+mfjmMF3o593yJob5y9pFGqmciCkFFjISM16XQMEfKyaV0mjjOvsQuD2xHNlwTctj6+x+0Q191L7+bqS4zMTOoJ/cZwF+7K+LL2beC/MmlhGt2/UHZStVef2z0obWAtlk5+j7BFQY23//PZYraUssh/XM9dByTTnfg3rJmNBcN+qpsx9+YQnPxm33nrryqM5hLPN/CXDk+H3gduG+UgnhVUNNq1lsMNUzUH++PbRdf3wsg8W6A7qAA/c+PDYNZCt7wa6HwVDfVlZ1JqvL6K70mlQukzAsFkHdIUtC7gSRUfby1DfKV/ndmg/c9tYwC4WjtW0Bt3LJkWl9XCdLHBZWPGkwP7JPchXx2mdalrrxDyf1AN03UbaRzxh83sq+cx1Us4E8Xb0vpkF5FSq0vs+GkyY/TQ2qiBtYei9aB4IrNEB3UHEzv0skAJ0Z2fl56lPLNMO1NHP11kzQSFdlAdJG4I/o+noLkUZD3AR/RuvME2I82dt1jYq9Sz7qFmDjbIyvmop8MWBtxDW2UfOglzTaF9re4BG22RNUyljQBMcnC+Upe1hilFJ4IeGHu3PAndB3yr5Uk+nejIfbUzwrJXHPZ6UVbtlwT0JTNBHOjFUQVpRUbHuUAXpGsUmdKO+rEVmfjCmBal2pcclWssgv2BWnroPWINVjTDTggZRllRzYL9gPD8TrcMsLUWts2dmhDYUz6w+wsx3qNqtsggymhRD+zvz2SpKLg8tT9ucmbaZOyNzB7B2y/m3o/tu1cfK695LFk3JiphGlx0S19lNE+8rozaxZstp1F87KVRBukZxGM0OBkD/AM1MxSyIk/lIIyihZbEgjV2VlOMI9PtO+cPTj4QnggzDgiuZX4zbyR9h9qFmpjWQC1g1rXeiu9t+1mb1l6r/Mps8NMCk7ea0TD8rTYz6LBnNS4NdWVmDAoTqX+cJchldLnNgKskHNMHR+ynNjOSJsvYj33GM6WrhjmLhm/nY1eUwSQEKtEtEJ1zmuKiCtKKiYl2h+kjXCVgL41mXtUMg12oiz7AZmTXA0IAybbNk8rMpPSV51EzMTDHVotWcDCyj+3PMGnXnQMiwlTWRj58r6tqNbuBoQdJl7WKE2RvIzF9mGoSmmFGOSvmB3q92HpD8Wd2QNFoWa4T8/jXYpaZ4xijQ9kf+feiOmUzDjGO2KLh+ZTpwOzOmwCjjb1xUQbpGcRi5aQ50fUQBpSiVTFAVBNn+k5GnZFpngpipM2FuZfQf/hDYxxnCnu9lH+8Suqb9jKThCHBp8uGPSvuMP7LFJB3QFTYqyAPcr6UPltkUMQkw/Yejzhm7YhE9+pu+l1LUH5KG+2UQO4DfC7sPeAzpyrDS5BGI95KxEDgar3uLZuMiysvcKdr/7G+dFKograioqBgD1bRfw9ClFmxaZ1oD8yvDFZC5AJQTylqPmoMciCkRzTkIFsGBMFNLa7VZ6yu5IzKyu7Z/NklXcidkEXBdlcPQvTlLZjq3S62CrM2DOI0aSOQVRJlGrBrcfjRBMqC70k0tGn5/m+k6a/ucbtA2euqOQCEdu1ki/QL6NcrlJB33vQbOuN6SO2AHepo7a/1LAM7EZFAF6RqFozuQ+eWzIGJhtUzp1a8U0I+doYM9W4rJZZUI6MEMYBJ2RjRX00zdAZk7gT9cXaigKH38nF+FDX+8WXR/B4Db6bm0rqzPSr5UXUkE5JQ3XeKYLUiI+7uTunSXLX4u7iOugwn6PJYyulPJT840NWU4RLu4vEGuhUBpGXGUx23J3EEl9sM4qFH7ioqKigmgaqRrFLzWfgldk60U9Z5rj/ehP5ARUEL/oGAJa6hZIEMjqKyBqbbHKLkMdE086Dgzp2cpTWg23CdxX7fKyyK4ocVn97gvOZqfmbmqPUXbSsta+V3GxiFcZtSjHN+4zsc70Nu+kN8Ha77c7kzTjHTzlG4B3b7k919ihHDblFPLGigHTDn4qJov6FgXIWSLBTQdvwt+r5NA9ZGuYRi69I0gMbOPlAfhDHoE8hk0a42Z0MxCLdubUj9Kvq8mE+g42pj5WNk0zUjtKrh4fTSbdjvQM1lLwj4+4pJgz0x/fR7OA3SFArc9Mwnjw8024IhnABpBl/lLdZUOu2CA7rvJ+nIazRhh4cBun4wKxsJmti0788VukXwcpd+PLh0tcwHpMbuAlI7F6dQNEnm4XZxGV3BxX2bvstKfKioqKtYIqka6hnEEuebCZqJqUBxE0aV3rClkhHpdzw+6p8TwUhBEtaYsiMBQojhrrmzCzqOrabDLQvcXyCLHGiDLAhqhdZW05QycP4swcz/upuOMKK6BKebE6nvmNioZnfNzW7J9GLi/QmsrabRct2qYGa9XGRkZvzTun97+ny+Uy+NCNU0139lVkuWpGulDDBvQNaczQjR/LBzZD1pQFqmN88jDPiU1WTMhXaL1qJkNdAcyf0iZ71UFIbsguCwVFhmVKs7Dr7yPrqtPjwURf3D8nBzp5zq5rDAzs0gx0O0z7gs1azN/H9B1DWiZccwCi90upb0CgH5/eUbCZzaDska4PKbMscDmyS/z3Zf2IC3575VpERSueUqjYzlz+ZT2gFgpatR+jSMbSLqah8FaEQscXXGS+aRYEMdgzz5+xjJ6Aa4DdF3bpTStLKAUebIAka5yKmk6mn+e7rGAyiaCKC/zJbNQ5RU3rMEqVxboF9Lczig3BOQBdH/wjdMuI6ciDQrC6Koxfn6+zhsuZxNL1Ml+UV06zJNHJmR1ghzUT3Pt8Xwhj1o3Sm0qcVozn+1a1EjN7E0Ang/gHnff2V6bAfB2NN0zD+BF7n7voHI2TKg9FRUVFauC8JEO+xsRfwrgIrl2BYAb3f3xAG5szweiaqQtDtOxapGsBWX+PqAbaWdtaRAVJ7ST0LjY9Jlr/zP9in1/6lPUOrId0jMzOzMn9Tkz32Xk35+kA7puh8xMDK2TNS8ljgNdTTnzUaq2pGCtaRo9t0PkZY0qW43F+fkdRxtLvkDuS0ieaIu6bViL5Wi+WhSZtQA5zpgiYSnxmnr2JfOzZe6oSKd0KK1bmQGTpj8Bk9NI25+Gn5PLF6P30+7XoPlZ91cOKqcK0hYb0TXnebCWNi1h85lNHh3Q/PFkdB3ddGQRvQHOwo8/UBUaeo9NSK5TBVGUM4uu0OLnDKhAVL4rH5eEEh9zIEhN0+xD1IAIt585sfzBZn5hfp7Ml8wujFl0Jx9G+MbjXsnHmE1KMdmyL56Fb+ZPVJeP+h/VB8p1R7t0BVlGedIJiid8bvN29FOzos5snGr/HS1WIWp/trvfDQDufreZnTUsQxWkFRUV6w4jBps2m9mtdH61u199LNpTBWkLQx6pB3L6B2sX56NZD56ZLqwFRhmRXwMNWZ1q/sX17WgI9Zw2m/F3ohc4YQ1Eza2SOc9plD40hW6wiYNypSAOl8ksCO4P1q5LVC69XtpGj8tlTKO7GozbzPl1ZY5q3dx/2ftXrZODdLzxDOfnelhT1k1WdKWSUsuA/sARk+qVpsR1ZgHGsK74mSPPArqsBV0EEPkngRVopAeO8ueYv2Bm57Ta6DkA7hmWoQrSFocBnNwe86DmD0zNx0izG12/0hx6tBnNUzL/V2qy6TJA/mBB5ambouRj5AHP9ZdM/sij1Joom027zMzTlUmRT8tSniVPHtxmNrPZtcEulFn093NGf1pGT8CVKGpKH9N7/J4HmdzMXc3GhpYL5OyAOAf6TXPOz66Okv9S2Rkld5KyLkqUq2wiHBfH2LS/HsClAK5s/79rWIYata+oqFhXmGTU3szeBuCjAJ5gZneZ2eVoBOhzzOxOAM9pzweiaqQENoFAx5nJyTNwzOyRbje6mh8T1TNtK7TGYVoca7dK7ldzthQgUU1Dzb44Bh2z1pIRrbN7rN0xWINkt4Vqq8qIiOsZcwDo/7XRLHDFWnxoZpnbgSP4nIbfUYwL1miztfYadOJ+yYJfQP/7Y9cAa4slvmqcBzINVNNo/2dE/Rgj3DeZ5s3gYOskNdMJRu1/rHDrWSsppwpSAr9wHvA8qAcJuMz3xJubqJkNOmZ2QNSlxzxYuY4ldCOoXH6JnB1+LP6YStHeAAu+EDbZ8/CHmPlVgWZ1TLgXgO5G2VpPRgtikn7kV19eQKlJ3JZpOY/jrC3T6P9p6mGuERVQgRn07+yU5QGlibp3tP/nJS0/S5ZfN2rhPFyvugw4FrCbznXC4neWLbTYD2AbJoNjbNqvGOvOtDezk83sH83s42b2KTP7jfb6jJndYGZ3tv/PaK8/3cz+9Lg2uqKiYmKIJaLD/lYT61EjfQDAM919ycw2Afiwmb0XwAvQrEa40syuQLMaYSCJVpGZiWpK8qydaaqcVu8pOV7N9FLgZZiZFtoNKE9piaJyUjPuH9ejz8Rm9nZKp7xRrn8z5Y/ru1DuM86vLg8NnLFWni23ZE1TI+i8RFQDJ1kEn+uI/mMtmLXVUv8HQpvNIvKsnauVsIwex3gLuks82dpgrXGRrpf2eygxHbj/ot6sb0sWAJBbN+Og7v40Abi7o/feN7V/jvJqhGUA9w0rl3+zaTvyX1TkYx6EQWUJzKC7w05Galf6Dpt2KjyyXYF0xUpJsJfM7/jwo2w18VhgaB4gd2dkpiVTYZjZoJFovZa5UObQW5kU/bUo51FPgIVdaSWS5lcBnzE4OF+Ux8InBOQeSV/yK6sg30Np+P2r24Gj/rohSqRRYc++ZG4P09e4L7gOfoYd6FHrdPxklLmMHna0qIJ0AjCzjQBuQyPzrnL3W8wsXY3g7h8B8JHj19qKiopJowrSCcDdDwM4z8weCeA6M9s5JMtI0Cg80DWZWYNgEyk0nZh596Or0WYRdNXENFjDWkMWGVbWAJfJUedMo4nrnJfza7CH28UaNWs0ym7ggAjnV9M4C6qpCyHayGwITs/lR5u5z7Lgii7DVC0+oNFrNoW3I9d2gd77Z+1Q3zNre/w8JaZDaNcZiyBzL8X14MRG/2fuGO0ntrw0wBn55qWurJ95LO/HZFBN+wnD3b9kZjeh2b1lRasRzIx/pglA1zQsmUmZyRt5s6itCpgMISzDHCxFTafQNbMGRcpjZplH90MIhB8xE+5smioViz9iNhN51U7UG+kyv3AIO+5bPo5yt1OZC5KmJFTiGbL2s7AtCSwW5AylgvHEBOQCU4WL0pq4n0v+9xJlSp+f/b8M9feqSwnopyllfs3N6P5qKufXWIC+58DevXthZl+hS9Otu25FqIJ0TJjZmQAebIXoKQCeDeB1WOFqhPblnRbnG8x8E91n4ZfRdxgZlSZbLskcPKXosHaivkz1q8X/TJuM43lqS2ll1RZ0fXiRbie6H4wGeCItf9Q70FvNxZPCZnQ3R2auKn/8QHciCOxB/rHH82Q0Je5/3TRGqURxrgEirqP0zjQoVAowsh9Uubbqz+V7cayaauYLVs1RA0xA/9jlfuZJUalgkW4e/b7gjL6lE3S0ZRbAmdu24eDBg6dhDNSNnSeDcwBc0/pJNwC41t3fbWYfBXBtuzJhAcALj2cjKyoqjh2qRjom3P0TAL4juX4QK1yNwNhIx7qBBGuU2bZhkY61lWw/Tl7xomR0pQ8x8TyjsrAGHCYiR9RLPjp2E7B2wZSh3VInm5JqykY7b0ZX41HCfIDbrKYpt63ElGCNWtNxpDugtBx1r2Sr1jgPb8aRpckWMZT8ypw/tHt9H1o/m/+RLqNscZ1M62KXS4zrTPOeQ88K0b0fdGVZya8dvlhuc2mRyTioPtI1jEPId9JhQcgfhe6WxOasbgzBg71Ed2GwOa6DNttMIz5IFooMNVUjzaCVRfyxZYGzMBkXKD9PMpnJzoI480GyuV3yK4fPNAT//iSd9lmJdsPtZCG7A12uZsl8ZbNZuZsBHjNz6NG3eFLj9gD9gTB+ttJqtJILiPPHOM3oXLupneq7Dcy16XgMscLB44THBX8Lk0IVpBUVFRVjoGqkaxiGPFihJG7W4AKhdbAWm5k/fKxsgFLUX4NYrOlp4IjN/iwgkLUrI+5zmiXkz59Rl1ijLUX9OT9rKFvQ0zZLNBkNyGngg+vM+pk12Gk0gbAs6n4A3T7PAoxLcr1kuXCfH0D/87N2qM+m7VLTusRU4HejWiDnLz0bkI+zcA1Efna78Phha42121JdR4MabFqjMPQ+ygXkkVqNEquwKn18mcmlQrFkzuqHwvWHmb+AwRuOZIIpniszR9WvuoPuHWj/q4mpHz5POCwgos3sh45zjqDvoTz8USr9KfNll4QF+6iX0GUmKC+4RJni9KB77Mdk4cnl6jOrzz2jmSnYvaM+dx6P6tfkZ8l8uUB5wuM2qwLAKE2SrGCcmTzTSlE10oqKiooJoArSNYxhJHxGRtrnGTkzs+M86lKTJwtEqSmXmXyhNXDgJQt8KdQEZG2NuZsH6DrXXVoZMyg6qyvFNrfHB+i6BlsykxmSpuRm0ICS8kezqLeWx1CmRLYaTt0ZrAHydS6PsR15EC3eM4+t7Fi5q9OF/KrRZ3niHJSH33PGTmGodTMJVI10DYNNe/UlatQzjtWUz+4NGtRqPikLIOrmwVoSquom0E1UIn9AzUygJ+Rm0TVHGWrWZcwDFnALyCeFEOL7KE8WtWdmAdcTApLrZ8qYshO4rdzG7N2Uos7qjmF3BJfPQoXbqJOPvsMAr+BSQr5OcgFdkFAys3mcqbAsCTp9Fh6D2ThTKiC/80mhCtKKioqKMVEF6RpGpoXyDKxaY3Y9ytGgAtC/pR074XmJIWu0rFGpya/aIpuzrNEGWEuYRTdAoWZ+abkgP5Oak6B0qjkGShqtajCsXTE4PT8nt2VG2szHWl/G1GBtvRTQimPllQJdk5/zcEBJub78LMrgCIQ2mrk6lAGgASIug69nfcN9yc+yE93FGpyf85RYC5M07ScVtTeziwC8Hs2anDe4+5VHU04VpC14hmMzUVeZsLDRj3Va0gLdwRZpNb+mU5pVNnDVlOZIq97XjyfaV/KXZdQeLVPNXBVSbGaXTD5um05Mmcm4hG6/MlhgqiDI+i/ayP3Mgixro0bDdZKNNu9Bd/xkBHaNnOtihew5o4+y59mOfCcpdUHxuy0xHXgs8GQ7L3kGTVK6cQq3aVxMykfaLjO/Cs0P3N0F4GNmdr27//NKy9owgfZUVFRUrCqOjPA3Ap4CYI+7/4u7LwP4CzQbxK8YVSNtsQHdWZw1mmwZHmsgs+hqUKrRMdeQ8yjpOQvcKFcvc9iHKZoFuFhrmUGXK7ss9/jZWFPjNfxqEiM5Z1eHrtteStIEgq+6G3mALPJldZYCN5yHtehoV8aJVGYAa+SlwM8oz88a7Cy6QbHSO2dNL8ZIpkVmS5a5DKDnWmLNkcc8kmNlrHAfzBbazW6nkothXEzIR3ougM/T+V0AvvNoCqqCtMUjH/UozM3NpfdKJOIzR0gzKM8jRsxTwt69e7FtW/67jOcOyfvtK6jnaJ4zoM+4DYPbDaysbSUMe/6jxUr7fFB/rfQ5V9L32sK9e/fiyQP6fLWwd+/eFe89muB96DHnBuFkM7uVzq9296vp3DQDGs/BilEFaYsDBw5knbqmYWZfGXdvx+OB9dpuYP22fb22O4O7XzShou4C8Bg6fzR6jLwVofpIKyoqHqr4GIDHm9k2M5sCcAmaDeJXjKqRrm+sOy26omKtwN0PmdnPonEVbATwJnf/1NGUZUfxcykVawRmZkfzezcVDz3UsXJsUU37dQAzu8zMnq/X64dx7FDq8/WKOlaOLappfxxgZo8B8GY0jKIjaKKJry9db7O9qF2F8QV3f21bTt+qDDM7GcCHADwMzft9h7u/uk1/GYBnAPgagLsBbEKzYOVFLY9uUs+XrhYptHe12lTsWzN7E4DnA7jH3XdStqzP5wF8GcBhAIfc/QIzewKAt1O+xwL4dQBfWqVn62sT3dsI4FYA/+ruz1+t/n7Iwd3r3yr/ofkBv/Pb44ejoU5+64DrlwH48fb629v/GwF8Fs1HOwXg421aAzDdptkE4BYA39WeXwbgJe3xje3/VwH4jgk+W6ldpevHvE2D+rw9vxDA+QB2Ufq+Pm+P5wFsHvL8+wF80yo+W7FNAH4JwJ8DePdqjYGH4l817Y8D3P1ud7+9Pf4ygE8DOLd0vc12X2Rv/6erMrxBcKM3tX9s1t3f/v9i+38ZjfY6KZRWiwxaRXKs21Ts8/b8Q+hfvQr09/koeBaAz7r759rzY/5sJZjZowE8D8Ab5NZxa9OJiipIjzPMbA7Nr6LeMsp1QrYq49w270YzuwPAPQBucPdSGccCpXYV27vaGKFvB8EBvN/MbjOz/5DcvwTA28Zo3tGg1KbfB/BfsfY2SzrxcLxV4ofyH5oVdLcBeMEo1yXNC9H4GeP8JwD8gaR5JIC/B7BzFZ8pbdco7T3OfT4HMu0H5N/a/j8LjXviQro3hWaP6rNX+Zn62oTG5/tH7fWnozXt69+x+asa6XGCmW0C8FcA/szd3znseoKhqzLc/UsAbgIwqZUgo6DUromtIjlarKBvi3D3fe3/ewBch8ZlEfh+ALe7+xfGbesE2vQ0AD/YBqL+AsAzzeytq9muhxKqID0OMDMD8EYAn3b33x12vYB0VYaZnWlmj2zLOwXAswF85hg8xoraNeD6qmCFfVsq4zQze3gcA3gugF2U5MewymZ9qU3u/qvu/mh3n0PT1x9w95euZtseSqj0p+ODp6ExbT/Z+jKBJnJ6f3bd3d+jBXhhVYaZPRHANS3tZQOAa9393cf0aUZoFwBMahXJUSLtc3d/j5m9DY35u9nM7gLwand/Y1LG2QCua2QyTgLw5+7+twBgZqei2dfyZ47pU6ygTRWrh7qyqaKiomJMVNO+oqKiYkxUQVpRUVExJqograioqBgTVZBWVFRUjIkqSCsqKirGRBWkFRUVFWOiCtKKioqKMVEFacVAmNlhM7vDzHaZ2d/QqqmtZvaOEfLrL0rH9R8ys28dkvfjLVn+uGHU56x4aKMK0oph+Jq7n+fNhseLAF4ONOu73f1Hxyj3h9DsR5rCzL4Fzfi8sF36eFwwgeeseAigCtKKleCj6G3VN2dmu9rjU83sWjP7hJm93cxuMTPepf23Wu3yZjM728yeCuAHAfyPVtt9XFLXSwC8BcD727RR1ivM7J/buv6ivTZtZn9iZp9sr/9Ie/25ZvZRM7vdzP7SzKbb6/Nm9hvt9U+a2Te31/9t2547zOyfzOzh8pwnUz3/ZGbPaK9fZmbvNLO/NbM7zex3JtzvFWscVZBWjIR27f6zkG808p8B3OvuTwTwWgBPpnunAbjZ3Z+E5idQXubuH2nL+S+ttvvZpMwXo/n5jreh2QwkcAWa3dyfCOA/ttf+HwD3ufu3t9c/YGabAfwagGe7+/lofm7jl6icA+31/wngV9prvwLg5e5+HoDvRfNzHIzQxr+9bdM11vy0CwCc17b52wG82JqfNql4iKAK0ophOKXd5OMggBkANyRpvgfNVm1w910APkH3lgHEpim3odn3cyDM7N8A+KI3u8zfCOB8Mzujvf0JAH9mZi8FcKi99mwAV0V+d78XwHehcR38Q9v+S9H8/EcgttHjNv0DgN81s1cAeKS7H0IX34NGS4a7fwbA5wDsaO/d6O73ufvXAfyz1FVxgqMK0oph+FqroX0Tmo2LX56ksQH5H/TezjiHMdqOYz8G4JvbvTQ/C+B0AD/S3nseGqH5ZAC3mdlJbf26+46h+XWA89q/b3X3y+n+A9omb36k76cBnALg5jD5R3zOB+h41OesOEFQBWnFSHD3+wC8AsCvtBskMz4M4EUA0Ebiv32EIr+M5kfoOjCzDWh203+iu8+1+2leDODH2nuPcfe/R/MTGo9Es+P9+wH8LJVxBoCbATzNzLa31041sx0YADN7nLt/0t1fh8YVoIL0QwB+vE27A8AsgP87wrNWnOCogrRiZLj7P6H5KYtL5NYfATjTzD4B4JVozO/7MBh/AeC/tEEbDjZdiOang/+Vrn0IjZl+LoC3mtknAfwTgN9rfwXgNwGc0VK0Pg7gGe7+RTS/mPm2tl03o18wKn6ByvgagPcmz7mxrf/tAC5z9we0kIqHHup+pBVjow1EbXL3r7dC8UYAO7z+TnrFQwTVj1MxCZwK4O9bk98A/KcqRCseSqgaaUVFRcWYqD7SioqKijFRBWlFRUXFmKiCtKKiomJMVEFaUVFRMSaqIK2oqKgYE1WQVlRUVIyJ/x/T1HvopwusEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "dataset.counts.sum_over_axes().plot(add_cbar=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "WcsNDMap\n",
       "\n",
       "\tgeom  : WcsGeom \n",
       " \taxes  : ['lon', 'lat']\n",
       "\tshape : (150, 150)\n",
       "\tndim  : 2\n",
       "\tunit  : \n",
       "\tdtype : int64"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dataset.counts.get_image_by_coord({\"energy\": \"5 TeV\"})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [],
   "source": [
    "dataset.write(\"pks2155-simulated.fits.gz\", overwrite=True)"
   ]
  },
  {
   "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": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "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.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}