Andreas Specovius / ctools

Commit dad5f0ebf1aea7262187922d415e9a93fda736c1

Authored by Luigi Tibaldo
1 parent b70128d1

Classical analysis moved to notebook

Inline Side-by-side
Showing 6 changed files with 314 additions and 132 deletions
doc/source/users/tutorials/hess_dr1/classical.ipynb 0 → 100644
  View file @dad5f0e
  1 +{
  2 + "cells": [
  3 + {
  4 + "cell_type": "markdown",
  5 + "metadata": {},
  6 + "source": [
  7 + "# Performing a classical On/Off analysis"
  8 + ]
  9 + },
  10 + {
  11 + "cell_type": "raw",
  12 + "metadata": {},
  13 + "source": [
  14 + "<div align=\"right\"><a href=\"../../../_downloads/classical.ipynb\" download=\"classical.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
  15 + ]
  16 + },
  17 + {
  18 + "cell_type": "markdown",
  19 + "metadata": {},
  20 + "source": [
  21 + "**In this tutorial you will learn to perform a classical On/Off analysis of the data.**\n",
  22 + "\n",
  23 + "We start by importing the gammalib, ctools, and cscripts Python modules."
  24 + ]
  25 + },
  26 + {
  27 + "cell_type": "code",
  28 + "execution_count": 4,
  29 + "metadata": {},
  30 + "outputs": [],
  31 + "source": [
  32 + "import gammalib\n",
  33 + "import ctools\n",
  34 + "import cscripts"
  35 + ]
  36 + },
  37 + {
  38 + "cell_type": "markdown",
  39 + "metadata": {},
  40 + "source": [
  41 + "We will also use matplotlib to display the results."
  42 + ]
  43 + },
  44 + {
  45 + "cell_type": "code",
  46 + "execution_count": 5,
  47 + "metadata": {},
  48 + "outputs": [],
  49 + "source": [
  50 + "%matplotlib inline\n",
  51 + "import matplotlib.pyplot as plt"
  52 + ]
  53 + },
  54 + {
  55 + "cell_type": "markdown",
  56 + "metadata": {},
  57 + "source": [
  58 + "And finally we add to our path the directory containing the example plotting scripts provided with the ctools installation."
  59 + ]
  60 + },
  61 + {
  62 + "cell_type": "code",
  63 + "execution_count": 6,
  64 + "metadata": {},
  65 + "outputs": [],
  66 + "source": [
  67 + "import sys\n",
  68 + "import os\n",
  69 + "sys.path.append(os.environ['CTOOLS']+'/share/examples/python/')"
  70 + ]
  71 + },
  72 + {
  73 + "cell_type": "markdown",
  74 + "metadata": {},
  75 + "source": [
  76 + "## Preparing the binned observations\n",
  77 + "\n",
  78 + "The first step for a classical analysis is to bin the events in bins of energy for an On (signal) region, and one or several Off (background regions).\n",
  79 + "\n",
  80 + "We will use the event selection performed in the previous tutorial."
  81 + ]
  82 + },
  83 + {
  84 + "cell_type": "code",
  85 + "execution_count": 21,
  86 + "metadata": {},
  87 + "outputs": [],
  88 + "source": [
  89 + "obsfile = \"obs_crab_selected.xml\"\n",
  90 + "\n",
  91 + "phagen = cscripts.csphagen()\n",
  92 + "phagen['inobs'] = obsfile\n",
  93 + "phagen['inmodel'] = 'NONE' # assume that the source is pointlike\n",
  94 + "phagen['ebinalg'] = 'LOG'\n",
  95 + "phagen['emin'] = 0.66\n",
  96 + "phagen['emax'] = 30.\n",
  97 + "phagen['enumbins'] = 20\n",
  98 + "phagen['coordsys'] = 'CEL'\n",
  99 + "phagen['ra'] = 83.63\n",
  100 + "phagen['dec'] = 22.01\n",
  101 + "phagen['rad'] = 0.2\n",
  102 + "phagen['bkgmethod'] = 'REFLECTED'\n",
  103 + "phagen['use_irf_bkg'] = False\n",
  104 + "phagen['maxoffset'] = 2.0\n",
  105 + "phagen['stack'] = True # treat observations separately in joint likelihood fit\n",
  106 + "\n",
  107 + "phagen.run()"
  108 + ]
  109 + },
  110 + {
  111 + "cell_type": "markdown",
  112 + "metadata": {},
  113 + "source": [
  114 + "Note that we have used the hidden parameter *use_irf_bkg* to specify that no background model should be used in the computation of the background normalisation factor, which is simply assumed to be the ratio of the solid angles of On/Off regions.\n",
  115 + "\n",
  116 + "## Likelihood fit\n",
  117 + "\n",
  118 + "The next step will be to fit a model for the source to the On data (along with a model for the background derived directly from the On/Off data). The model for the source is"
  119 + ]
  120 + },
  121 + {
  122 + "cell_type": "code",
  123 + "execution_count": 22,
  124 + "metadata": {},
  125 + "outputs": [
  126 + {
  127 + "name": "stdout",
  128 + "output_type": "stream",
  129 + "text": [
  130 + "=== GModels ===\n",
  131 + " Number of models ..........: 1\n",
  132 + " Number of parameters ......: 6\n",
  133 + "=== GModelSky ===\n",
  134 + " Name ......................: Crab\n",
  135 + " Instruments ...............: all\n",
  136 + " Instrument scale factors ..: unity\n",
  137 + " Observation identifiers ...: all\n",
  138 + " Model type ................: PointSource\n",
  139 + " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
  140 + " Number of parameters ......: 6\n",
  141 + " Number of spatial par's ...: 2\n",
  142 + " RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)\n",
  143 + " DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)\n",
  144 + " Number of spectral par's ..: 3\n",
  145 + " Prefactor ................: 5.7e-17 +/- 0 [1e-24,1e-14] ph/cm2/s/MeV (free,scale=1e-17,gradient)\n",
  146 + " Index ....................: -2.48 +/- 0 [-0,-5] (free,scale=-1,gradient)\n",
  147 + " PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
  148 + " Number of temporal par's ..: 1\n",
  149 + " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n"
  150 + ]
  151 + }
  152 + ],
  153 + "source": [
  154 + "srcmodel = os.environ['CTOOLS']+'/share/models/crab_nobkg.xml'\n",
  155 + "\n",
  156 + "print(gammalib.GModels(srcmodel))"
  157 + ]
  158 + },
  159 + {
  160 + "cell_type": "markdown",
  161 + "metadata": {},
  162 + "source": [
  163 + "We can now fit this model to the data. Note that we have to specify WSTAT as statistic to derive the background model from the data without prior information."
  164 + ]
  165 + },
  166 + {
  167 + "cell_type": "code",
  168 + "execution_count": 23,
  169 + "metadata": {},
  170 + "outputs": [],
  171 + "source": [
  172 + "like = ctools.ctlike(phagen.obs())\n",
  173 + "\n",
  174 + "like['statistic'] = 'WSTAT'\n",
  175 + "like['inmodel'] = srcmodel\n",
  176 + "\n",
  177 + "like.run()"
  178 + ]
  179 + },
  180 + {
  181 + "cell_type": "markdown",
  182 + "metadata": {},
  183 + "source": [
  184 + "We can now look at the results from the optimisation and the fitted model."
  185 + ]
  186 + },
  187 + {
  188 + "cell_type": "code",
  189 + "execution_count": 24,
  190 + "metadata": {},
  191 + "outputs": [
  192 + {
  193 + "name": "stdout",
  194 + "output_type": "stream",
  195 + "text": [
  196 + "=== GOptimizerLM ===\n",
  197 + " Optimized function value ..: 16.457\n",
  198 + " Absolute precision ........: 0.005\n",
  199 + " Acceptable value decrease .: 2\n",
  200 + " Optimization status .......: converged\n",
  201 + " Number of parameters ......: 6\n",
  202 + " Number of free parameters .: 2\n",
  203 + " Number of iterations ......: 4\n",
  204 + " Lambda ....................: 1e-07\n",
  205 + "=== GModels ===\n",
  206 + " Number of models ..........: 1\n",
  207 + " Number of parameters ......: 6\n",
  208 + "=== GModelSky ===\n",
  209 + " Name ......................: Crab\n",
  210 + " Instruments ...............: all\n",
  211 + " Instrument scale factors ..: unity\n",
  212 + " Observation identifiers ...: all\n",
  213 + " Model type ................: PointSource\n",
  214 + " Model components ..........: \"PointSource\" * \"PowerLaw\" * \"Constant\"\n",
  215 + " Number of parameters ......: 6\n",
  216 + " Number of spatial par's ...: 2\n",
  217 + " RA .......................: 83.6331 [-360,360] deg (fixed,scale=1)\n",
  218 + " DEC ......................: 22.0145 [-90,90] deg (fixed,scale=1)\n",
  219 + " Number of spectral par's ..: 3\n",
  220 + " Prefactor ................: 4.58230817013557e-17 +/- 3.14037370499303e-18 [1e-24,1e-14] ph/cm2/s/MeV (free,scale=1e-17,gradient)\n",
  221 + " Index ....................: -2.6252889331251 +/- 0.0832081186871608 [-0,-5] (free,scale=-1,gradient)\n",
  222 + " PivotEnergy ..............: 1000000 [10000,1000000000] MeV (fixed,scale=1000000,gradient)\n",
  223 + " Number of temporal par's ..: 1\n",
  224 + " Normalization ............: 1 (relative value) (fixed,scale=1,gradient)\n"
  225 + ]
  226 + }
  227 + ],
  228 + "source": [
  229 + "print(like.opt())\n",
  230 + "print(like.obs().models())"
  231 + ]
  232 + },
  233 + {
  234 + "cell_type": "markdown",
  235 + "metadata": {},
  236 + "source": [
  237 + "The fit has properly converged and the results are consistent with those obtained with the unbinned analysis.\n",
  238 + "\n",
  239 + "We will also check the spectral residuals."
  240 + ]
  241 + },
  242 + {
  243 + "cell_type": "code",
  244 + "execution_count": 26,
  245 + "metadata": {},
  246 + "outputs": [
  247 + {
  248 + "data": {
  249 + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtcAAAG9CAYAAADeCVr8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xt8lPWZ///XZYxgUOIBrK7IBL9ELQpyiKLW1ba4FlfQeqgHkqrYNbqurbrd9mtrW5J286v9tl23Xa0aW8XW4LFYSz0fUGu1KnQRPFV80KRQRQFrQCOKeP3+mAOTZJLM+Z575v18PO4Huc9XwsydK5/5fK6PuTsiIiIiIpK77YIOQERERESkXCi5FhERERHJEyXXIiIiIiJ5ouRaRERERCRPlFyLiIiIiOSJkmsRERERkTxRci0iIiIikidKrqXimdk5ZrbCzHrMbK2ZXWNmu8T2tZiZm9kXko7fPratLqiYRUREpDQpuZaKZmZfBX4AfA2oBQ4DIsBDZrZD7LC3ge+aWVUwUYqIiEhYKLmWimVmI4FW4Mvufr+7b3H3TuA0ogl2U+zQ+4EPk9ZFREREUlJyLZXsCGA4sDB5o7u/C9wH/FN8E/BtYJ6ZVRc1QhEREQkVJddSyUYB6939oxT73ojtB8DdfwusA/6lSLGJiIhICCm5lkq2HhhlZtun2LdXbH+ybwGXE23tFhEREelHybVUsqeBD4CTkzea2QjgOOCR5O3u/hDwGnBhsQIUERGRcFFyLRXL3buJDmj8HzObaWbVsfJ6dwBrgF+lOO1y4OtFC1JERERCRcm1VDR3/3/AN4EfARuBZ4DVwAx3/yDF8X8Ani1qkCIiIhIa5u5BxyAiIiIiUhbUci0iIiIikidKrkVERERE8kTJtYiIiIhInii5FhGRrJjZcDN71syeN7MXzaw16JhERIIW6gGNo0aN8rq6uqDDEBHJ2NKlS9e7++ig48iFmRkwwt3fNbNq4EngYnf/Y6rj9cwWkTBL97mdama60Kirq2PJkiVBhyEikjEz6wo6hlx5tHXm3dhqdWwZsMVGz2wRCbN0n9vqFiIiIlkzsyozWwa8BTzk7s/02d9sZkvMbMm6deuCCVJEpIiUXIuISNbcfau7TwbGAIea2UF99re7e4O7N4weHepeMCIiaVFyLSIiOXP3d4DHgJkBhyIiEqhQ97kWkfRt2bKFNWvWsHnz5qBDqSjDhw9nzJgxVFdXBx1K3pnZaGCLu79jZjsCxwA/CDgsEZFAKbkWqRBr1qxh5513pq6ujmiRByk0d2fDhg2sWbOGcePGBR1OIewF3GRmVUQ/Cb3d3X8XcEwiIoEqqeTazD4PHA/sAVzt7g8GHFLOFixYwMqVKzM+r76+njlz5hQgIqlUmzdvVmJdZGbG7rvvTrkO5HP35cCUoOMQESklBe9zbWY3mNlbZvZCn+0zzezPZvaamV0G4O6/cffzgHOA0wsdWzFkk1jncp7IYJRYF59+5iKVxcz0vq9wxWi5ng9cBfwyviH2EeLVwD8Ba4DnzOy37v5S7JBvxfaXjXnz5qV9bGurJjkTERERCaOCt1y7+xPA2302Hwq85u6r3P1D4FbgRIv6AXCfu/8p1fUquWZqR0dH4i/iuro6Ojo6gg5JJCNmxhe/+MXE+kcffcTo0aOZNWtWRtepq6tj/fr1OR8jIiKSb0H1ud4bWJ20vgaYDnyZ6GjzWjMb7+7X9j3R3duBdoCGhobQzN1++nVPJ77ueuYBnrkx2jpds9snmHjiBUSmfy6xf0KK8zs6Omhubt52ja6uxHpjY2NhghbJsxEjRvDCCy/w/vvvs+OOO/LQQw+x9957Bx2WiIhI3gRV5zpVZyR395+6+zR3vyBVYl0Oup55gCUdVyTWe95+kyUdV9D1zAO9jmtpaUm0UpsZTU1N9PT09Dqmp6eHpqamXseJlLrjjjuOe+65B4BbbrmFM888M7Hv7bff5vOf/zyTJk3isMMOY/ny5QBs2LCBY489lilTpnD++ecTnXU76uabb+bQQw9l8uTJnH/++WzdurW435CIiEiSoFqu1wD7JK2PAV5P92Qzmw3MHj9+fL7jKoiWlhagZcD9Wz/8gGdubE20ZosUWqH69acztuCMM87gu9/9LrNmzWL58uWce+65/P73v0+cP2XKFH7zm9/w6KOPctZZZ7Fs2TJaW1s58sgj+c53vsM999xDe3s7AC+//DK33XYbf/jDH6iurubCCy+ko6ODs846qyDfn4iIyFCCarl+Dqg3s3FmtgNwBvDbdE9290Xu3lxbW1uwAIPW0tKCuyeWSCSS8rhIJNLrOJFSN2nSJDo7O7nlllv453/+5177nnzyyUSf7M9+9rNs2LCB7u5unnjiCZqamgA4/vjj2XXXXQF45JFHWLp0KYcccgiTJ0/mkUceYdWqVcX9hkRERJIUvOXazG4BPg2MMrM1wDx3/4WZXQQ8AFQBN7j7i4WOJSgtLS2JFr26ujq6urr6HROJROjs7AS2tSomty5OmzaN119/nS1btiS2VVdXM23atF7HqT62pCOT6jWFcMIJJ/Af//EfPPbYY2zYsCGxPdUfiPHuTqm6Pbk7Z599Nt///vcLF6yIiEgGilEt5Ex338vdq919jLv/Irb9Xnffz93/j7u3FTqOUtHW1kZNTU2vbTU1NbS1bfsRbBo2qt95kyZNYvbs2cRb62tra5k9ezaTJk3qdVyq+tiqMiKl5txzz+U73/kOEydO7LX9qKOOSrw+H3vsMUaNGsXIkSN7bb/vvvv4+9//DsCMGTO48847eeutt4Bon+1Uf7yKiIgUS0nN0JiusPW5Thav7HH55Zfz17/+lbFjx9LW1tar4sfqXacCcNv5h2d07VT9aFVlRErRmDFjuPjii/ttb2lpYe7cuUyaNImamhpuuukmINrSfuaZZzJ16lSOPvpoxo4dC8CECRP4z//8T4499lg+/vhjqqurufrqqwfsRiUiIlJoFuZ+ug0NDb5kyZKgwxhUPOHN5GP40697mpfe2MiEvUZmdK8Jax+MDZ7MXJhfB5Kel19+mU9+8pNBh1GRUv3szWypuzcEFFIgwvDMFslVvAubfq+Wn3Sf26FsuS53J05W3V8RERGRMFJyXYLmTB/LnOljMz6vtfXBXoMnIb0BlCIiIiKSH6FMroPoc71gwYKUgwVLXVtbG83Nzb0moOk7gHIw+nhLREREJH2hTK7dfRGwqKGh4bxi3TOXxDpV9Y9C6juwcebMmSxcuBCIVhmZMWMGr7322pAl/JKritTV1fUbeCkiIiIivYUyuQ5SqoGJg7Xunn7d0wWPKW7TsFHs/MH6ftsnTZrUr2RfX33/eFCVERGRyqFPKUXyR8l1jkqpdTeXEn4tLS1DVhrp6emhqakpMVMe6EEsUsnMbB/gl8CewMdAu7v/JNioRESCFcrkOsg618kt0V3PPMCSjiu2rXd1cfa5X+J/HllJZPrnALIqqZeLl97YmHFr+YQCxSLS19q1a7nkkkt47rnnGDZsGHV1dfz3f/83++23X1rn77TTTrz77rsFjlIy8BHwVXf/k5ntDCw1s4fc/aWgAxMRCUook+sg+lwDsZbdlkGP2frhBzxzYyvP3LitP/M3/licGeNyKeGX7yoj+ohR+nJ3TjrpJM4++2xuvfVWAJYtW8abb76ZSK63bt1KVVVVkGFKBtz9DeCN2NebzOxlYG9AybWIVKxQJtdhk01ZvWzvk20Jv+i/2/4gmDZtGq+//jpbtmxJbKuurmbatGkaCClZWbx4MdXV1VxwwQWJbZMnT+axxx7jM5/5DHvttRfLli3jpZde4vOf/zyrV69m8+bNXHzxxb36/3/1q19l8eLF7Lrrrtx6662MHj06iG9H+jCzOmAK8Eyf7c1AM5CYWVNEpJwpuc5AudaQTjUQMj4Asm+Vkb4DIzUQMpxaF73IS69vzOs1J/zDSObNPnDA/S+88ALTpk1Lue/ZZ5/lhRdeYNy4cQDccMMN7Lbbbrz//vsccsghnHLKKey+++689957TJ06lR//+Md897vfpbW1lauuuiqv34dkzsx2An4NXOLuvV5Y7t4OtEN0hsYAwhMRKSol1znItYZ0qdBASAnaoYcemkisAX76059y1113AbB69WpWrlzJ7rvvznbbbcfpp58OQFNTEyeffHIg8co2ZlZNNLHucPeFQcdTitRNTqSyhDK5DnJAY7J4S2w8aYxEIuoCIaEwWAtzoRx44IHceeedKfeNGDEi8fVjjz3Gww8/zNNPP01NTQ2f/vSn2bx5c8rz4kmLBMOi/wG/AF529/8KOp5yooRcJLxCmVwHNaAxlcbGxrJIprOtMlKuXWUk/z772c/yzW9+k+uvv57zzou+dZ977jkef/zxXsd1d3ez6667UlNTwyuvvMIf//jHxL6PP/6YO++8kzPOOIMFCxZw5JFHFvV7kH4+BXwRWGFmy2Lbvunu9wYYU0VTUh4sjTkSCGlyLfmVS5WRvjTdugzEzLjrrru45JJLuOKKKxg+fDh1dXV8/vOf73XczJkzufbaa5k0aRL7778/hx12WGLfiBEjePHFF5k2bRq1tbXcdtttxf42JIm7Pwno4wMRNOZItlFyLTlXGUmWS1cZ/cVf/v7hH/6B22+/vd/2eEs2wLBhw7jvvvtSnh+vcf29732vMAGKiKQpnW5pGnNUmZRcS86SS/PFJQ9yfO211/od07eEn/7iFxERkXKwXdABSHhtGjYq63MbGxsxs8TS1NTUqysJbPuLP/k4ERGRUuDuvZZIJJLyuEgk0us4KX9quZas5VLCT0REwk/jZLYJqjyv/g9KTyiT61IpxSeqMiIiIgIqzyvbhDK5LqVSfJUs1yoj+ZpuXUREwqVcW1sbGxsTybUahipXKJNrKQ3ZVhn5jyv+lLfp1qF8H9IiIiISPkqupegG66sdT5TfeeedfvvUVzv8qqqqmDhxIu5OVVUVV111FUcccUTG1znnnHOYNWsWp556agGizM1OO+2UKBkoIlKq1DBVOKoWIiUl09HUfWtjJ69L6dlxxx1ZtmwZzz//PN///vf5xje+UfQYPvroo6LfU0REKoeSawmtgWpjK8EOh40bN7LrrrsC0clhZsyYwdSpU5k4cSJ333134rhf/vKXTJo0iYMPPpgvfvGL/a7z7W9/m3POOYePP/6Ye++9lwMOOIAjjzySr3zlK8yaNQuIDqBtbm7m2GOP5ayzzmLz5s3MnTuXiRMnMmXKFBYvXgzA/PnzueiiixLXnjVrFo899hgQbZG+/PLLOfjggznssMN48803AfjLX/7C4YcfziGHHMK3v/3tgvysREQkPNQtRAKRS5WR5Alq+tJsWGm67zJYuyK/19xzIhx3xaCHvP/++0yePJnNmzfzxhtv8OijjwIwfPhw7rrrLkaOHMn69es57LDDOOGEE3jppZdoa2vjD3/4A6NGjeLtt9/udb2vf/3rdHd3c+ONN/LBBx9w/vnn88QTTzBu3DjOPPPMXscuXbqUJ598kh133JEf//jHAKxYsYJXXnmFY489lldffXXQ2N977z0OO+ww2tra+PrXv87111/Pt771LS6++GL+9V//lbPOOourr74605+aiIiUGbVcS9GdOHlvJuw1MugwJADxbiGvvPIK999/P2eddVaiK9A3v/lNJk2axDHHHMPf/vY33nzzTR599FFOPfVURo2KTli02267Ja71ve99j3feeYfrrrsOM+OVV15h3333Zdy4cQD9kusTTjiBHXfcEYAnn3wy0Qp+wAEHEIlEhkyud9hhh0RL+LRp0xKVAP7whz8k7pWqZV1ERCpLKFuuVec63LKtMtLa+mCv+tiqjZ2DIVqYi+Hwww9n/fr1rFu3jnvvvZd169axdOlSqqurqaurY/Pmzbj7gDNzHnLIISxdupS3336b3XbbbchPKEaMGJH4eqBjt99+ez7++OPE+ubNmxNfV1dXJ2Kpqqrq1Xdbs4eKiEhcKFuu3X2RuzfX1tYGHYoEoLW1ldbWVqZNm0Z1dXWvfcm1sePLggULAopUBvPKK6+wdetWdt99d7q7u9ljjz2orq5m8eLFiT+aZsyYwe23386GDRsAenULmTlzJpdddhnHH388mzZt4oADDmDVqlWJP6xuu+22Ae991FFHJfrmv/rqq/z1r39l//33p66ujmXLlvHxxx+zevVqnn322SG/j0996lPceuutAOrvL6GkgeEi+RXKlmupTJuGjepVHzuX2tgSjHifa4i2Ht90001UVVXR2NjI7NmzaWhoYPLkyRxwwAEAHHjggVx++eUcffTRVFVVMWXKFObPn5+43he+8AU2bdrECSecwL333svPfvYzZs6cyahRozj00EMHjOPCCy/kggsuYOLEiWy//fbMnz+fYcOG8alPfYpx48YxceJEDjroIKZOnTrk9/STn/yEOXPm8JOf/IRTTjkltx+QSJENNDAc0MyCIlmyMA/2amho8CVLlhTlXvEay8lTdktxxQdA9q2PPVitTv2/bfPyyy/zyU9+MugwCurdd99lp512wt35t3/7N+rr67n00kuDDivlz97Mlrp7Q0AhBaKYz+xSkk094WxrEA91XrZdmIa6Xr7jDLNif2/6PyiedJ/barmWUElVZeS0a58CSFl9ZEJRopJScf3113PTTTfx4YcfMmXKFM4///ygQxIRkQqj5FpC48TJewcdgpS4Sy+9tCRaqiuFmd0AzALecveDgo5Hhta3lVIDw0XyT8m1hEY2VUZaWx8sUDQiAswHrgJ+GXAcJavvYMG2traS6svc1tZGc3MzPT09iW01NTW0tbUFGJVIuIWyWoiIiATP3Z8A3h7ywAqV7Syyxaze0djYSHt7e2I9EonQ3t5eUn8AiISNWq6lIsQHNmaivr6eOXPmFCAakcphZs1AM8DYsZnXtx/i2kDpDMhKZ7Bg31lk+8aeS/WObFvJGxsbE/GoK4hI7tRyLWVt07BRA+4baip1lfCLJguaIEVy4e7t7t7g7g2jR48OOpySE3+PxZempqZeXTRgW0KefFxf2baSi0j+qeVaytrqXaN1ivuW7wMSiXWqMn3ZtHTL0Kqqqpg4cSJbtmxh++235+yzz+aSSy5hu+0G/ju/s7OTp556Sp8iSMnLZrBgtn+8ZtNKnipGEck/JddS9lKV70uWaQm/sNQwzVUhBmLtuOOOLFu2DIC33nqLOXPm0N3dPegfM52dnSxYsEDJdY7MbASw2d23Bh1LpUhnsGC21Tv0iZJI6QpltxAzm21m7d3d3UGHIiXuxMl7M2GvkUGHETrF+Ih5jz32oL29nauuugp3p7Ozk3/8x39k6tSpTJ06laeeitYvv+yyy/j973/P5MmTufLKKwc8Tnozs+3MbI6Z3WNmbwGvAG+Y2Ytm9kMzq8/DPW4Bngb2N7M1ZvalXK9ZTrIZLNjW1kZNTU2vbamqd7h7ryUSiaS8XiQS6XWciBRB3zdomJZp06Z5sbS0tHhLS0vR7ieFdfPNNzvggEciEb/55pt77R/s/zt+XrHk634vvfRS2vcaasnWiBEj+m3bZZddfO3atf7ee+/5+++/7+7ur776qsff34sXL/bjjz8+cfxAx5WyVD97YIkX8PkIPA58G5gEbJe0fTfgFODXQFMhY+i75Pv/qtjvxWxlGudQz6eBzqmpqen1Pq2pqUnr3GxizDbOXO4XFsX83vR/UFzpPrfVLUQqTiaj8ft2V1i+fHni61122YUZM2YwadKkXsfks8pIqdfIzQePtaZt2bKFiy66iGXLllFVVcWrr76a8vh0jxOOcfctfTe6+9tEE+tfm1l18cOSoWRTvSP+XIifF4lECvq8yKWqSbmLP9MKTf8HpSuU3UJEMpHNaPxUVUSWL1/OokWLEuvd3d0sWrSoV8IN+asyEsTo//hf3UN9xJwvq1atoqqqij322IMrr7yST3ziEzz//PMsWbKEDz/8MOU56R5X6eKJtZn1H7Hb5xgpD8kJVWdnZ14TrHxVNYHi1vFOjj3s8vl/IIWllmuRAQxWpi9uy5YtLFy4kIULF2Z0XiqlNPq/GLO2rVu3jgsuuICLLroIM6O7u5sxY8aw3XbbcdNNN7F1a3Tc3c4778ymTZsS5w10nAxonpnVEO0O8ifgVnf/e8AxSYUKU2trqdVRl/BQy7WUvb59odIZ+HPatZU9SK5Qs7a9//77TJ48mQMPPJBjjjmGY489NlEK8cILL+Smm27isMMO49VXX2XEiBEATJo0ie23356DDz6YK6+8csDjZEAObAYeAPYBnjKzg4MNScImm+eou6u1NY+y/T+QAKTTMbtUFw1olGykM/DntGuf8tOufarXeZFIJOXgvkgkkjgmn6+TdO6XiXQGNPaFBrzkRRADGuML8GKf9f2AR4tx776LBjQW9pxinpfuAMpUz7B0lnwq9s+yWHIZxJrLQMhKlu5zWy3XUnHSbZWN18eOL3t+9hyqdhjW65iqHYax52fPSRyTT+mW5Cqk+INCQm29mU2Lr7j7q4CmSpScpPsc7Zt0qLU1f7L9hFGzeRaekmupSEMN/ElVHzsy/XM0NF6WWK/Z7RM0NF5GZPrn0rpnph93FqprhlScrwA3m9nNZvZ/zawD+EvQQUn4ZTOAshQaDcpJOv8H2XTNybdK6+6jAY0iKcyZPpY508f233H+4diN0fJ8721Y2293a+uDeY0jm5Jcg3H3inrAlYKgWuHMzGKfZD5vZpOBY4CDgMXALcnHBBKgVKRilwwUCYJarkUqxPDhw9mwYYM+ci0id2fDhg0MHz48iNsvNrMvm9lYd//A3e9x9x8AvwSmm9lNwNlBBCaVrZAlA6W/bLrmSG7Uci0Vq9IeIGPGjGHNmjWsW7cu6FAqyvDhwxkzZkwQt54JnAvcYmbjgHeA4UAV8CBwpbsvCyIwEQlOMUqtVjol1yIFEPTMjqlUV1czbty4gl1fSou7bwZ+BvwsNhPjKOB9d38n2MhEJEjqmlN46hYikqHBPjbbNGxUv225zuyoj+kkV+6+xd3fUGItIqCuOYWmlmuRPFq961Ruv+CIIY/L58yOIiIiUjoyTq7NbFdgH3dfPuTBAktuhBV3Zn7exFOhYW7+45GCOnHy3twedBAiIiISmLS6hZjZY2Y20sx2A54HbjSz/ypsaGVixZ2wdkVm56xdkV1CLoGbM31sViOz1WotIiJSHtJtua51941m9i/Aje4+z8zy2nJtZvsCl8fudWo+rx24PSfC3HvSP/7G4wsXixSdRmZLkMzsC8D97r7JzL4FTAX+093/FHBoUkI0rkMkf9JNrrc3s72A04gmwGkxsxuAWcBb7n5Q0vaZwE+IloT6ubtf4e6rgC+ZmZpspaxkMjK7b5WRdBS6yoiE3rfd/Q4zOxL4HPAj4Bpgej4unup5no/rDiV5qua6ujpVOxCRkpFutZBW4AHgNXd/LtbKnLq0QW/zidZaTTCzKuBq4DhgAnCmmU1IO2KREBpqZHaqKiPpGqjKSDYqbYraCrE19u/xwDXufjewQz4uHNTzvKOjg+bm5sR6V1cXzc3NvRJuEZGgpNty/Ya7J4ryuvuqdPpcu/sTZlbXZ/OhRJP0VQBmditwIvBSmrGIlJ3Vu04F4LbzD8/ovGxauqXi/M3MriM6/fkPzGwY+SvDmtHzfMOGDcyfP7/XtgMPPJBDDjmELVu2pEyOJ0+ezJQpU6ipqeG0004D4I477ujVzQqgp6eH888/n4cffjix7eijj+bwww9n//33Z/369fzud7/rd/2jjjqKfffdl7Vr13L//ff32z9jxgz22WcfVq9ezSOPPNJv/8yZM9lzzz1ZtWoVTzzxRL/9s2bNAmC//fbr970DnHTSSdTW1vLCCy+wZMmSxPZzzjkn8X3V1NSwbNkyli3rP+dPY2Mj1dXVPPfcc7z44ouJ8+L3iq8/9dRTvPrqq73Ora6uTvyh//jjj/OXv/yl1/7kn/nDDz/MmjVreu0fOXJk4uv777+ftWvX9tq/++67M3v2bAAWLVrEhg0ben1/yccvXLiQjRs39jp/zJgxHHPMMQDcfvvt/f7Px40bx9FHHw1E/+DasmVLr/377bcfRxxxROJ+2bz2li1bRk9PD7ff3n+oekNDAwcddBDd3d3cdddd/fYX67W37777ctRRR/X7/mbNmsWoUaP485//zNNPP93v/JEjR7Jx48Z+r7240047LaPXXl/Jr72+r8t8vPZOPvlkIPPXHsCee+7JzJnRtt9sXntDSfcB+z9pbkvH3sDqpPU1wN5mtruZXQtMMbNvDHSymTWb2RIzW6KZ5kREhnQa0U8eZ8bqXO8GfC1P1075PE8+IPmZ3Tf5ydZ7772X0XYRkWKywQYxmNnhwBHAJcCVSbtGAie5+8FD3iDacv27eJ/r2OCaz7n7v8TWvwgc6u5fzjT4hoYGT/XXViHEWwjnzZuX2YnxwYnZDGjM5BwpefHuFqnec6df9zQvvbGRCXuN7LdvMBPWPghk8bocwGAxSn6Z2VJ3bwg6jlxk+jzP1zO7rq6Orq6uftsjkQidnZ05Xz/fsnlfZfteLPZ7uNhxlvPPMlth+P7C8rMcSrrP7aFarncAdiLafWTnpGUjkG1FjzXAPknrY4DXs7yWSFk4cfLeGSfWIoMxs01mtjHFssnMNg59hbQE8jxva2ujpqam17ZSrsCTzSyrmplVJLwG7XPt7o8Dj5vZfHfv30yQneeAejMbB/wNOAPIqNSBmc0GZo8fPz5PIYkEa870scyZPjbj81pbH4z9m3uVEVVfKC/uvnMRbpPz8zwbmVTgEREptnQHNA4zs3agLvkcd//sYCeZ2S3Ap4FRZrYGmOfuvzCzi4j2AawCbnD3/j3hB+Hui4BFDQ0N52VyXqisXZFdvWvN7FhRNg0bxc4frM/q3OQqIwNVXwCUsJSB2My69cDw+DZ37z/6LkPu/lGuz/NsNTY2JpLrUuwKEhS1dvemRgMJQrrJ9R3AtcDP2VbWaUjufuYA2+8F7k33OhVnYpY9buIzQSq5LkmF+KWXS5WRlpaWQWeG7OnpoampKZHAwMDfQ7n0pytHscm/LibaZWMZcBjwNDBo40jbAbkwAAAgAElEQVS69DyXUpVLo4GScslFusn1R+5+TUEjyUDZdwtpmJtdgqyZHUWkv4uBQ4A/uvtnzOwAonMXiJSVdGr0p9NooE/yJFfpJteLzOxC4C7gg/hGd3+7IFENIdBuIZkmsGtXRKc/FykxLS0tiSojYau+IBnZ7O6bYxMEDXP3V8xs/6CDEikV+UrKReLSrXN9NtG6qE8BS2NLcWrghd2eE7Pv5iFSJGGrviAZWWNmuwC/AR4ys7tRhSbJk1KqahKPJb5EIpGUx0UikV7HieRbWi3X7j6u0IGERqnXntZAyIrz0hsbOf26/rNvDabv/NSqvlC+3P2k2JctZrYYqAXuCzCkwGmMQGVoa2ujubm51+x6qRoN+r4OwvZJXqm/jiux/3paybWZnZVqu7v/Mr/hpKfs+1xnSwMhK86Jk/ce+qA0qfpCeTKz76TYPBn4brFjEYkrRkKYbaNBukm5DK1S+68POkNj4iCz5KnOhwMzgD+5e6D9HUIxQ2MY3Hh89n3D1eIdSgO9nsMw01e5KNYMjWb21aTV4cAs4GV3P7fQ9+4r389svV4rT7FmaOzo6Cj7T/IK8T5Ip/96KmF5L6b73E63W0ivqWzNrBb4VZaxSalRi7fkoBI/8gsTd/9x8rqZ/Qj4bUDhiISCPsmTXKRbLaSvHqITEkg5UOm/ijXQzI4DbU81q2MlfuQXcjXAvkEHISLBKkRrcdj7r+dLWtVCzGyRmf02ttwD/Bm4u7ChiUihbBo2KqvzGhsbiZV0w8xoamrq1S8RtpWsSj5OgmNmK8xseWx5kejz+6dBxyUi5S+oSlRB/+5Jt+X6R0lffwR0ufuaAsSTFg1oFBnaYP3pBprZMT5jY6qxBQO1ZkvJm5X09UfAm+7+UVDBiEjlqNRKVGm1XLv748ArwM7ArsCHhQwqjXgWuXtzbW1tkGGIVJyWlhbVkQ0JM/t3M/t34JSk5XTgK7HtIiIFl5xId3Z2ln1iDel3CzkNeBb4AnAa8IyZaWYUkQqnyWdK2s6xpQH4V2Dv2HIB/Uudi4hInqTbLeRy4BB3fwvAzEYDDwN3FiowESmswSafSbU9VTZWqR/5hYG7twKY2YPAVHffFFtvAe4IMDQRkbKWbnK9XTyxjtlA+lOni0iJGWjymdOufSrja6lkVckbS++ufB8CdcGEIiJS/tJNru83sweAW2LrpwP3FiakoWlAY8gtuRFWZPmhhyatyYs508cyZ/rYjM5pbX0w9u/AAxvTLeEnRfUr4Fkzuwtw4CTgpmBDCidNPiOFpNdX+Ri09dnMxpvZp9z9a8B1wCTgYOBpoL0I8aWkAY0ht+LObRPQZGLtiuyTcslZtuX7AFauXJnHSCQT7t4GnAv8HXgHmOvu3w82KhGRoQVdUi9bQ7Vc/zfwTQB3XwgsBDCzhti+2QWNTsrXnhNh7j2ZnaNJawI1UPk+UAm/UufuS4Gl+bymmX0BaAE+CRzq7vmb11xEJMSG6jdd5+7L+26MPUTrChKRiIjkzMyejP27ycw2Ji2bzGxjHm7xAnAy8EQeriUypI6OjsTXdXV1vdZFSslQyfXwQfbtmM9ARCR/9EtI3P3I2L87u/vIpGVndx+Zh+u/7O5/zj1SkaF1dHTQ3NycWO/q6qK5uVnPNilJQyXXz5nZeX03mtmXyPNHjCKSH/olJMnM7AtmtnPs62+Z2UIzmxJ0XCKDife1jS9NTU309PT0Oqanp4empqbQ9suV8jVUn+tLgLvMrJFtyXQDsAPREedS6dauyLwv9NoV0T7Xkhfp/FKJ/xKKl8yD7EakD1YbG9Kvjy1F9W13v8PMjgQ+B/wIuBaYPtSJZvYwsGeKXZe7+93p3NzMmoFmgLFjM6tQMxRVVRCRUjRocu3ubwJHmNlngINim+9x90cLHtkgVIqvREzMcpLOPSdmf64EZqDa2FLytsb+PR64xt3vjk0kMyR3PybXm7t7O7HqUg0NDcqGJS19/3Cqq6ujq6ur33GRSET19aXkpFXn2t0XA4sLHEva3H0RsKihoaFflxUpooa5qjldAor1S2iw2ti3XxD9N1UlkXh9bAnM38zsOuAY4AdmNgxNAiYh09bWRnNzc6+uITU1NbS1tQUYlUhqesCKlJm2tjZqamp6bdMvoYp2GvAAMNPd3wF2A76W60XN7CQzWwMcDtwTm2hMpCAaGxtpb982vUYkEqG9vZ3GxsYAoxJJLd0ZGkUkJOK/bOL9qyORCG1tbQX9JaS+r6XL3XuIzVEQW38DeCMP170LuCvX64ikq7GxMfFcU1cQKWVKriVcshlAWYFTpuuXkMRZdMRrI7Cvu3/XzMYCe7r7swGHFoi+ZSoL/YeniFQedQuR8Jh4auZVRjRlellS6a2M/Ixo140zY+ubgKuDCyc4KlNZmdxdn65JUanlWsIjmwGUmjJdZLq7TzWz/wVw97+b2Q5BB1UMxSxTKSISp+RaRIqitbU143NqRo/haxd+qQDRVJQtZlYFOICZjQY+DjYkkdKnP7IkW6HsFmJms82svbu7O+hQRGQINaPHZH1uz7o1eYykYv2U6MDDPcysDXgS+P+CDak44t0B4kskEkl5XCQS6XWciEguQtlyrTrXIuGRbctzNi3d0p+7d5jZUmBGbNMXgIqcIlW1kkWkGEKZXItI5RgsyR5oX319PXPmzClUSKFgZiOBfwP2Bn5LdGDjRcAiYBlQcaP4gihTKVJO4uMYSvkTnlKoCKTkWspfNuX7oCJL+JWSTcNGsfMH67M6d+XKlXmOJpR+BfwdeBr4F6ITx+wAnOjuy4IMLEgqUylSfMVKxgeqCAQUNcFWci3lbeKp2Z23dkX0XyXXgVm961Qg9ZTqLS0tAMybN6/fPnUnSdjX3ScCmNnPgfXAWHffFGxYIiL5UaoVgZRcS3nLpnwfRFu61eIt4bYl/oW7bzWzvyixFhEpvFBWCxEpuGwmrAFNWiOl5GAz2xhbNgGT4l+b2caggxMRyVWpVgRSy7VIKrm0eIuUAHevCjoGEUlPKQzCKwelUhFIybWIlKyX3tjI6dc9PeD+VPsmFDIgqVhKfqRQSmUQXqnJ5j1XKhWB1C1ERErSiZP3ZsJeI/tt73rmgcTXv/vmSb3WRQphoOQn+Ze/SLrMrNfS1NTUq6UVtg3CSz6ukuTynktOpDs7OwP5A0Ut1yJSkuZMH8uc6WN7bevo6KD5th8m1nvefpMVt/2QL8+oTzxAW1sfLGqcUn5KtQKBSLkqt/ecWq5FpGRl08ITL9MnIlKKSnUQnuRPKJNrM5ttZu3d3d1BhyJSsvRAFsmOkh8ppra2NmpqanptC2IQXpDK7T0XyuTa3Re5e3NtbW3QoYhIAWXzwFXLdXGY2Q/N7BUzW25md5nZLkHHVChKfqSQGhsbaW9vT6xHIhHa29srejBj2N9zoUyuRaQyhf2BW2YeAg5y90nAq8A3Ao6nYJT8SKGVwiC8UhL295ySaxEJjbA/cMuJuz/o7h/FVv8IjAkynkJT8iNSXGF+zym5FpFQCfMDt4ydC9yXaoeZNZvZEjNbsm7duiKHJSJSfCrFJxJmS27Mfrr1iadmNwtlSLS2tmZ8Ts3oMXztwi8VIJpwMrOHgT1T7Lrc3e+OHXM58BGQsgCtu7cD7QANDQ2lOwJJQqGUB7GJxCm5FgmzFXfC2hWw58TMzlu7IvpvGSbXNaPH0LNuTVbnZnteuXL3Ywbbb2ZnA7OAGa6sR0QEUHItUhqybYGOJ9Zz78nsvBuPz/xeIZFty3M2Ld2VzMxmAv8XONrde4Y6XkSkUqjPtUgpiLdAZ2rPidHuHSLFdxWwM/CQmS0zs2uDDkhEpBSo5VqkVGTTAi0SEHcfH3QMIiKlSC3XIiIiIiJ5ouRaRERERCRPlFyLiEhFcXeVdBPJUEfHtmqbdXV1vdalNyXXIiIiIjKgjo4OmpubE+tdXV00NzcrwR6ABjSKiCTJtiTfvHnz8hyJiEgwzGzIY3p6emhqaqKpqSmxTZ8IRanlWkSE6OQzIiIiuVLLtUi+rV2R+SQt2cyyKHmlac9FRKL6tkDX1dXR1dXV77hIJEJnZ2eRogqPkmm5NrMRZnaTmV1vZo1BxyOSlYmnZpckazIYEREpUW1tbdTU1PTaVlNTQ1tbW0ARlbaCtlyb2Q3ALOAtdz8oaftM4CdAFfBzd78COBm4090XmdltgHrJS/g0zI0uIiIiZaKxMdrmGe9fHYlEaGtrS2wvNUH3/S50y/V8YGbyBjOrAq4GjgMmAGea2QRgDLA6dtjWAsclIiIiImlKTqQ7OztLNrEuBQVtuXb3J8ysrs/mQ4HX3H0VgJndCpwIrCGaYC+jhLqriJStbPqGTzxVLfMiIiKDCCKJ3ZttLdQQTar3BhYCp5jZNcCigU42s2YzW2JmS9atW1fYSEXKVTZ9w9eugBV3FiYeERGRMhFEtZBUxRPd3d8DhmwSc/d2oB2goaFBBRVFspFN3/BMW7lFREQqUBDJ9Rpgn6T1McDrAcQhIpnKpisJqDuJiMgQgh6EJ/kTRLeQ54B6MxtnZjsAZwC/zeQCZjbbzNq7u7sLEqCIpJBtmUF1J5Ey4e5KgERkSIUuxXcL8GlglJmtAea5+y/M7CLgAaKl+G5w9xczua67LwIWNTQ0nJfvmEVkANmWGVR3EhERyUJY/5gtdLWQMwfYfi9wbyHvLSIihWNm3yNa6elj4C3gHHdXFz8RqXihLHmnbiEiIoH7obtPcvfJwO+A7wQdkIhIKQhlcu3ui9y9uba2NuhQREQqkrtvTFodAYTz81sRkTwLolqIiIiUATNrA84CuoHPDHBMM9AMMHbs2OIFJyISkFC2XIuISOGZ2cNm9kKK5UQAd7/c3fcBOoCLUl3D3dvdvcHdG0aPHl3M8EVEAhHKlmszmw3MHj9+fNChiIiULXc/Js1DFwD3APMKGI6ISCiEsuVafa5FRIJlZvVJqycArwQVi4hIKQlly7WIiATuCjPbn2gpvi7ggoDjEREpCUquRUQkY+5+StAxiIiUolB2CxERERERKUWhTK41iYxIZXP30E6LKyIi5S2UybUGNIqIiIhIKQplci0iIiIiUoqUXIuIiIiI5ImqhYhIaVpyI6y4M7tzJ54KDXPzG4+IiEgaQtlyrQGNIhVgxZ2wdkXm561dkX1SLiIiA9Jg8vSEsuXa3RcBixoaGs4LOhYRKaA9J8LcezI758bjCxOLiIhIGkLZci0iIiIiUoqUXIuIiIiI5ImSaxERERGRPFFyLSIiIiKSJ0quRURERETyxMJYUsXMZgOzgdOBlQMcVgukW6sv3WOHOm6o/aOA9WnGVMoy+dmW8n3zcb1srqHXZuGE6bUZcffRhQimVJnZOqBrgN1BvC+GOkbvi9K7b67XDMsze6hj9NoM5r7pPbfjNQvLbQHa833sUMelsX9J0D+XYv9sS/m++bheNtfQazM8r5FSv285LUG8L4Y6Ru+L0rtvrtcMyzN7qGP02izt+5Zzt5BFBTh2qOMyuWeYBfV95vu++bheNtfQa7NwyuW1WYmCeF9ket+wKqf3Ra7XDMszO9P7hlU5vTYTQtktJKzMbIm7NwQdh0hfem2K9Kf3hZQqvTZLWzm3XJei9qADEBmAXpsi/el9IaVKr80SppZrEREREZE8Ucu1iIiIiEieKLkWEREREckTJdciIiIiInmi5FpEREREJE+UXAfIzEaY2U1mdr2ZNQYdj0icme1rZr8wszuDjkWkVOiZLaVKz+zSouQ6z8zsBjN7y8xe6LN9ppn92cxeM7PLYptPBu509/OAE4oerFSUTF6b7r7K3b8UTKQixaNntpQqPbPDS8l1/s0HZiZvMLMq4GrgOGACcKaZTQDGAKtjh20tYoxSmeaT/mtTpFLMR89sKU3z0TM7lJRc55m7PwG83WfzocBrsb8sPwRuBU4E1hB9WIP+L6TAMnxtilQEPbOlVOmZHV56OBTH3mxr7YDoA3pvYCFwipldQ4HnuRcZQMrXppntbmbXAlPM7BvBhCYSGD2zpVTpmR0C2wcdQIWwFNvc3d8D5hY7GJEkA702NwAXFDsYkRKhZ7aUKj2zQ0At18WxBtgnaX0M8HpAsYgk02tTpD+9L6RU6bUZAkqui+M5oN7MxpnZDsAZwG8DjkkE9NoUSUXvCylVem2GgJLrPDOzW4Cngf3NbI2ZfcndPwIuAh4AXgZud/cXg4xTKo9emyL96X0hpUqvzfAydw86BhERERGRsqCWaxERERGRPFFyLRXJzM4xsxVm1mNma83sGjPbJbavxcy2mNm7ScvXY/seM7PNffYdHux3IyIiIqVCybVUHDP7KvAD4GtALXAYEAEeig0QAbjN3XdKWv5f0iUu6rPv6eJ+ByIiIlKqVOdaKoqZjQRagXPd/f7Y5k4zOw1YBTQFFpyIiIiEnlqupdIcAQwnOtNagru/C9wH/FMQQYmIiEh5UHItlWYUsD5WzqivN2L7AU4zs3eSln9IOu6nSdv/VPCIRUREJDSUXEulWQ+MMrNUXaL2iu2HaO3QXZKW5BmwvpK0fWrBIxYRCQkz22pmy5KWy4KOKc7M7jSzfc3smVhsfzWzdUmx1g1w3n+a2ff6bGsws+Wxrx8xs9rCfwcSFkqupdI8DXwAnJy80cxGAMcBjwQRlIhImXjf3ScnLVfkesEBGkMyvcaBQJW7r3L36e4+GfgO0cHr8Vg7Bzj9FuD0PtvOiG0HWABckGuMUj6UXEtFcfduogMa/8fMZppZday14g5gDfCrAMMTESlLZtZpZq1m9qdYGdQDYttHmNkNZvacmf2vmZ0Y236Omd1hZouAB81sOzP7mZm9aGa/M7N7zexUM5thZncl3eefzGxhihAagbvTiPM4M3s6FudtZjYiNgPiZjObFjvGgC8At8ZOuxuYk8vPR8qLkmupOLGyet8EfgRsBJ4BVgMz3P2DIGMTEQm5Hft0C0lu8V0f60p3DfAfsW2XA4+6+yHAZ4Afxj5JBDgcONvdP0v008Y6YCLwL7F9AI8CnzSz0bH1ucCNKeL6FLB0sMDNbA/gMqK/C6YCy4GLY7tvIdpaHb/W6+7+FwB3Xw/sHJ8rQUSl+KQiufsvgF8MsK9lkPM+XaCQRETKwfuxLhepxFuUl7Kta96xwAlmFk+2hwNjY18/5O5vx74+ErjD3T8G1prZYgB3dzP7FdBkZjcSTbrPSnHvvYB1Q8R+BDABeCraOM0OwJOxfbcAj8cmFEvuEhK3LnaPd4a4h1QAJdciIiJSDPFPBreyLf8w4BR3/3PygWY2HXgvedMg170RWARsJpqAp6oG9T7RxH0wBtzv7l/su8PdO83sdeAfgZOAaX0OGR67h4i6hYiIiEhgHgC+HOvHjJlNGeC4J4FTYn2vPwF8Or4jVs3pdeBbwPwBzn8ZGD9ELE8BR5vZvrFYRphZfdL+W4CfAi+7+9r4RjPbjmgZ19VDXF8qhJJrERERyZe+fa6HqhbyPaAaWG5mL8TWU/k10UHnLwDXER0r0520vwNY7e4vDXD+PSQl5Km4+5vAl4DbzOx5osn2fkmH3A4cxLaBjHGHAk+6+9bBri+Vw9w96BhEREREBmVmO7n7u2a2O/As8Kl4C7KZXQX8b2w8TapzdwQWx87JaxJsZlcTnRvh8XxeV8Ir1Mn1qFGjvK6uLugwREQytnTp0vXuPnroI8uHntkiEmbpPrdDPaCxrq6OJUuWBB2GiEjGzKwr6BiKTc9sEQmzdJ/b6nMtIiIiIpInSq5FRERERPJEybWIiIiISJ4ouRYRERERyRMl1yIiIiIieaLkWkREREQkT5RcB6CjowMzw8yoq6ujo6Mj6JBEREQkQPG8QMJPyXWRdXR00NzcnFjv6uqiublZCbaIiIhIGQj1JDJhMdRfoj09PTQ1NdHU1JTYFuaZM0VEylH8Wa7ns4gMRi3XBbZgwYKiniciIiIiwVHLdYGtXLmSlpaWxPqVV15Jd3d3v+Nqa2u59NJLe50nIiIiIuGi5LpI5s2bB8D48eNpbm6mp6cnsa+mpoarr76axsZGAFpbWwOJUUSkEqh7h4gUkrqFFFljYyPt7e2J9UgkQnt7eyKxFhEREZHwUst1ABobG5VMi4iIZEifOkgYqOVaRERERCRPlFyLiIiISMFU2gQ5Sq5FRKQfM7vBzN4ysxcG2G9m9lMze83MlpvZ1GLHKCJSipRci4hIKvOBmYPsPw6ojy3NwDVFiEkkY9m2mlZaa6vkj5JrERHpx92fAN4e5JATgV961B+BXcxsr+JEJyJSupRci4hINvYGVietr4lt68XMms1siZktWbduXV4DUMuiSHEV+z0X1k8dSiq5NrPPm9n1Zna3mR0bdDwiIjKgVL+5+tVHc/d2d29w94bRo0cXISwRkWAVPLkeaFCMmc00sz/HBsNcBuDuv3H384BzgNMLHVsYBf3XmIhIzBpgn6T1McDrAcUiIlIyitFyPZ8+g2LMrAq4muiAmAnAmWY2IemQb8X2i4hIafotcFasashhQLe7vxF0UCIiQSt4cj3AoJhDgdfcfZW7fwjcCpwYe0j/ALjP3f+U6nqF7L9X6jo6OhJf19XV9VoXEcknM7sFeBrY38zWmNmXzOwCM7sgdsi9wCrgNeB64MKAQpUQ0qewwdP/QeEENf15qoEw04EvA8cAtWY23t2v7Xuiu7cD7QANDQ0lMf9poaZjbW1tTXy9fPlyFi1alFjv6upi7ty5LFy4kEmTJiW219fXM2fOnLzGISKVx93PHGK/A/9WpHCkRGk68uDp/6D0BJVcpxwI4+4/BX5a7GBy0bc1ua2tjcbGxpyuWV9fn9Y1tmzZwsKFC1m4cGFiW0tLS073FhGR/Mom+ck2YVKiJRK8oKqFlMVAmI6ODpqbmxPrXV1dNDc359xdQy3PIiJSTOoiIJI/QSXXzwH1ZjbOzHYAziA6OKakxR8+8aWpqYmenp5ex/T09NDU1JQ4JtuWZHfvtUQikZTHRSKRxDFqtRYREREJVjFK8fUbFOPuHwEXAQ8ALwO3u/uLGVxztpm1d3d3FyboFBYsWFC0e6XS1tZGTU1Nr201NTW0tbWldb5aJUREREQKr+B9rgcaFOPu9xIdbZ7NNRcBixoaGs7LJbZMrFy5sl/L8JVXXkmqBL+2tpZLL700sV5fX5/z/eN9sJuamoBoi3U++neLiIiISP4ENaAxtObNm5f4evz48TQ3N/fqGlJTU8PVV19dkKS3sbExq+sWYtCliIiIiPSn5DoHpdqanM8Sfhp5LiIiIpI+Jdc5amxsTCTXnZ2dgcaiEn4iIiIiwQqqWkhOghjQOJh4tY6g5buEn2aEFBERyYwKCEgoW66DGNAYFn2T/Lq6Orq6uvodF4lEEi3t8W4kmhFSREREJDehTK4lfW1tbSkHXSaX8FN3EhEREZH8CGW3EElfY2Mj7e3tifVIJEJ7e3uvZFotzyIiIiL5oeS6AiQn0p2dnSlbqTUjpIiIiEjulFxXiEwHXeY6I6SIiIhIJQplcl1q1ULKUTrdSURERESkt1Am1+6+yN2ba2trgw6lrKXTnUREREREtlG1EBlUKdTvFhEREQmLULZci4iIiIiUIiXXUlI0s5WIiIiEmZJrEREREZE8UXItIiIiIpInoUyuVYqvPHV0dCS+rqur67UuIiIiEgahrBbi7ouARQ0NDecFHYtkbsGCBaxcubLXtuXLl7No0aLEeldXF3PnzmXhwoVMmjQpsb2+vl7TtYuIiEjJCmVyLeG2cuXKtKZO37JlCwsXLmThwoWJbZpyXaQ4zGwm8BOgCvi5u1/RZ/85wA+Bv8U2XeXuPy9qkCIiJUjJtYiI9GJmVcDVwD8Ba4DnzOy37v5Sn0Nvc/eLih6giEgJC2Wfawm/lpYW3D2xRCKRlMdFIpHEMWq1FimaQ4HX3H2Vu38I3AqcGHBMIiKhoORaSkJbWxs1NTW9ttXU1NDW1hZQRCIVbW9gddL6mti2vk4xs+VmdqeZ7VOc0ERESpuSaymITCeDaWxspL29PbEeiURob2+nsbGxEOGJyOBSvXm9z/oioM7dJwEPAzelvJBZs5ktMbMl69aty3OYIiKlR8m1lIzkRLqzs1OJtUgazGxErI90Pq0BkluixwCvJx/g7hvc/YPY6vXAtFQXcvd2d29w94bRo0fnOUwRkdITyuRada5FpFKZ2XZmNsfM7jGzt4BXgDfM7EUz+6GZ1efhNs8B9WY2zsx2AM4Aftsnjr2SVk8AXs7DfUVEQi+UybW7L3L35tra2qBDkRRymQwmPnhRRAa0GPg/wDeAPd19H3ffA/hH4I/AFWbWlMsN3P0j4CLgAaJJ8+3u/qKZfdfMTogd9pVYQv888BXgnFzuKSJSLlSKT3LW2tqa+DrdyWBEJGvHuPuWvhvd/W3g18Cvzaw615u4+73AvX22fSfp628QTfBFRCSJkmvJWn19fVr9olNNBqOpzUWyE0+szWyeu7cOdoyIiBRfKLuFSGnIZRpyTWEukrN5ZvYDM7vezP7VzHYNOiAREVFyLTlKnggm3clg8t2nOtOyfyJlwoHNRPtF7wM8ZWYHBxuSiIgouZa80mQwIkXzirvPc/c73f2bRGdQvDLooEREKp2Sa8krTQYjUjTrzSxRW9rdXwVUSFpEJGBKriXvijkZTC5l/0RC7ivAzWZ2s5n9XzPrAP4SdFAiIpUulMm1JpERiCbWzc3NifWuri6am5uVYEtZs9gAA3d/HpgM3BLbtRg4M/kYEREpvlAm15pEpnK1tLQkBjA2NTXR09PTa39PTw9NTU2JY8yMBQsWBBStSEEsNrMvm9lYd//A3e9x9x8AvwSmm9lNwNkBxygiUrFU51oKohCzLNbXZzer88qVK/MciUigZgLnAreY2TjgHWA4UAU8CFzp7ssCjE9EpKKFsuVaKtOcOXMyLvvX0tJS3CBFCszdN7v7z51E6VoAABweSURBVNz9U0AEmAFMdfeIu5+nxFpEJFhKriW0ci37p/rYEnbuvsXd33D3d4KORUREopRcS2ip7J+IiJQDVb4qL0quJdSKWfZPREQk31T5qvxknFyb2a5mNqkQwYiEgbqTiFQetSxKviRXs8qk8pWER1rJtZk9ZmYjzWw34HngRjP7r8KGJlI4+kUpYWdmXzCznWNff8vMFprZ1KDjKkdqWRSRTKTbcl3r7huBk4Eb3X0acEzhwhIpHP2ilDLxbXffZGZHAp8DbgKuCTimsqCWRSmk5KpX6Va+KkR5WymcdJPr7c1sL+A04HcFjEckY0M9eJInntEvSikjW2P/Hg9c4+53AzsEGE8o6FMrKTW5Vr6S0pNuct0KPAC85u7Pmdm+gGbmkLKWamZH/WKWEvI3M7uOaKPHvWY2DA1SH1S6n1qpZVGKSZWvyk+6MzS+4e6JQYzuvkp9riUM6uvr+00kc+WVV9Ld3d3v2NraWi699NLEet+ZHQf6xQzoIShBOI3obI0/cvd3Yp8ufi3gmEpKOp9AxT+1ampqSmzrmyi3tbXR3Nzc6xMvtSwWX9/Gjba2trJ59jY2NiZeg52dncEGI7nr+xd6qgX4UzrbirUAs4H28ePHe7G0tLR4S0tL0e4nhXPzzTd7TU2NA4mlpqbGb7755sQxLS0tvfZnsnR0dAT43UlYAEs8oGdoUMu0adPy88OLib/nhtqf6ZLKzTffnNgfiUR6PS9yjTNf55Tzeek8t4OOsRLOC0OMuZyXxnXTem4P+hGimR1uZl8FRpvZvyctLUDVYOcWkrsvcvfm2traoEKQEEvnI7j6+vqsr9+3xVskn8xsk5ltTLFsMrONQcdXSvr+wsule4dq6hdXNoNKpTRVYnfKofrn7QDsRLT7yM5Jy0bg1MKGJlI4Q/2inDNnTla/mPt2QRHJN3ff2d1Hplh2dveRQcdXyjRwTKS4KrU616DJtbs/7u6twGHu3pq0/Je7q3lOKop+MUupiU3qdaiZHRVfgo6plGngWHhk07iRSratppXY2poPKmMZle7I8mFm1m5mD5rZo/GloJGJlBj9YpZSYmb/AjxBtJJTvKJTS5AxhYG6d4RTNo0b2baaVmprq+RPusn1HcD/At8iOho9vohUFP1ilhJyMXAI0OXunwGmAOvydXEzm2lmfzaz18zsshT7h5nZbbH9z5hZXb7unQ61LFaWdBo3sm01VWtretJ5z+VznEOYP3VIN7n+yN2vcfdn3X1pfCloZCIiMpjN7r4Zoomuu78C7J+PC5tZFXA1cBwwATjTzCb0OexLwN/dfTxwJfCDfNw7HWpZrExq3AhOtu+5bLtThv1Th3TrXC8yswuBu/j/27v/KLvr+s7jz5eQiENhsCH+aAI3ZBlohSLoiIfqQrsBNu0Gsxj5YWZa4aAje9bdLmp78NiVydqp9pgeapWqU5VYzwAN7LAwyoYapHA4sJqwhJhIY2iaNLMRDLFNaEYl4nv/uN97ufPz3rnzvfd+v3Nfj3O+Z+73e7+fz/c9c7/3M+/5zOd+PvCz0sGI+HFDojJrgun+WjbLiVFJpwD/C/iWpH8GDqRU94UUFw3bAyDpLmA18P2Kc1bzyjCUe4DPS1LM8MY6dOgQGzZsGHfsnHPO4W1vexvHjh2b8hfg+eefzwUXXEBHRwdXX301AHffffeUPYsf/OAH2bx5c/nYJZdcwkUXXcTZZ5/NCy+8wDe+UVxg+LrrrgNgw4YNXHzxxSxfvpznnnuOTZs2Tbr+ihUrOO2009i/f/+4ciUrV67kDW94A3v27OHRRx+dVH7VqlUAnHXWWZO+d4Arr7ySzs5OduzYwdatW8vHS9caGxujo6ODbdu2sW3btknle3p6WLBgAVu2bGHnzp2TYiztP/744/zgBz8YV3bBggXlBLX0c6iMsfJnvnnzZkZHR8eVP/nkVz4/u2nTJp577rlxzy9atIgrrrgCgJGREQ4dOjTu+6s8f3h4mCNHxk92s3TpUi699FIANm7cyNjY2Ljv74wzzuCSSy4BiknV7bffPq78zTffzPPPP89EixYtYv369dPeex/96EfHxTqxXMmOHTs499xzOXz4MPfee++47w1g165dk+69ShPvvYmvXeW999BDD00qX7r3li9fzsUXXzzp/lq1ahWnnnoqu3bt4oknnphU/uSTT+bIkSOT7j2ARx55pPwzP//889m9e/e0vfmf/exnOeeccwDKr0fp59DT08Pu3btZv349R48eZdGiRaxZs2bStSbGXrr2xOuV3uNHjhxheHh43M8Bam8brrnmGlauXAnUfu/NRq3J9fuSr5VDQQJYPqurmZlZKiLiyuRhv6SHgU7gf6dU/RJgf8X+KPD26c6JiJ9LOgwsAl6oPElSH9AHsGTJklSCO3r06KyOW3u68cYb+cxnPjMuMVq4cOGk5G6iNWvWsGHDBl566aVZlZvP5vLeuvzyy9m3bx8wPvmeyXTJbLU4MtM21DIZdla3tBckmIkXkbESZpicPkv3yUxxWusxx0VkgE9Mtc2lzoq6rwK+XLH/u8DnJpyzE1hasf8PwKKZ6k2rzS4UClMu/lIoFGoqX+97o5nl8hBjveWauSBPvddq9qJBzS432++v2e+5eq831zirqbXdrmnMtaTfm2qrpazZfFN685i12NGK7WWK46OXpVT3KHBaxf5SJg85KZ8j6XiKPedNGSroaTGzoZ4PjjV7TGy947TzMr67Wa9Bs99z9V4vM21DLRk48LmK7a+APcA9tZRt5Oaea8uaLN0nuOc600h5+XPg1cCDKdV1fNLOn0FxMbGngXMmnPOfgS8mj68FNlarN802Oy89i/XEmYfe1lqXI2eKXsRatrTirLdMHso18jWY7nrNfM+14n1QTa3tdr0Nbydwfz1l09ycXFvWZOk+cXKdbQ1Irl8L7E6xvt8BfkBxuMfHk2P/A3hX8vgEitO0Pgt8F1herc602+z5kvzMtcxcypXKVktG6k2SnVynVy6vP/88lauh3pra7Vo/0DjRGNBVZ1kzM5sjSd+j+AsE4DhgMfDJtOqPiAeAByYc+0TF459SHJttiVrmPi7NsNDb21tzvfWUma5cMT94xXRDBIBUhkJMvN6yZcvKH26rVCgU2Lt375yvZ5P5NWi+Wsdcj0i6P9m+CewC7mtsaGZWryxMom8Ntwq4ItkuB34lIj7X2pAs6+pdMGViz1y9i4NkZkxsDvk1yI9aF5FZD/xZsv0JcHFETFqxy8zSN9vVwLIyib41hqQPS/owsKZiuwb4r8lxa5F6k5+5lplLolWvehO0WlZatNr4NciumpLriHgE+HvgJIrj+l6auYSZNYuX7m07JyVbN/CfKM43vQS4keJqipYR9SQ/jZ4lIa2kfC4JWl5m4sg6vwbZVeuwkKspfmDlKuBq4DuS3tPIwKrEc4WkwcOHD7cqBDOzloiIdRGxDjgVeEtEfCQiPgK8leKUeZYR9SQ/9SZM9ZabyxABJ2it59cgm1TLv4kkPQ1cFhE/SvYXA5sj4s0Njm9G3d3dMXHJzkZZt24dALfccktTrmf51Ij7pNTLPNV79Y477mD37t3jjt16661M9YdnZ2cnN910U3m/q6uLtWvXphanzY6kJyOiew7l/x54c0T8LNl/NfB0RPxqWjGmLe02e6b3Rt7LNfNaQ0ND5Q89FgoFBgYGak7S/LNsr3J5iHEu5Wqot6Z2u9Yx168qJdaJQ7Moa2YNMjGxBlixYgULFiwYd2zBggWsWLGialnLla8D35XUL+kW4DvA11ock+WQez/T4w+TG1DzVHybJD0I3JnsX8OEKZrMrHUm9pRX64kq9bBbfkXEgKRNwDuTQ9dHxFOtjMmsnTV6WkPLjxmTa0lnAq+PiD+Q9G6KjbiAJwD/OWaWkqmGdwBs3769/PiUU05hxYoVnHfeeVXr6+npKSfXnrd0/oqIJ4EnWx2HWTuqd17ztIcqWPZUG9rx58CLABExHBEfjoibKPZa/3mjgzNrF9Ml1iMjI+X9w4cPMzIyMi7hhuLYaWsfkh5Lvr4o6UjF9qKkI62Oz8ys3VUbFrIsIrZPPBgRWyUta0hEZm2qv7+/6jnHjh1jeHiY4eHh8jH3grSXiHhn8vWkVsdi1s688qFNp1rP9QkzPPeaNAMxM7PaSbpK0knJ4z+SNCzpglbHZdauvPKhlVRLrrdI+sDEg5JuwOP8zGal2uIt/f39LV1xzXLnv0fEi5LeCfx7ijOFfLHFMZm1La98aCXVkuv/Blwv6e8k/VmyPQK8H/j9xodn1r7cC2JVvJx8/Q/AFyLiPmBhC+Mxa3ue1tCgypjriHge+A1JvwWcmxz+ZkR8u+GRmbW5UqNc7+IONu/9P0lfAi4F/jRZRMbrD5iZtVhN81xHxMPAww2OxcwmmMuUeh4yMu9dDawE1kfEv0h6I/AHLY7JzKzt1bqIjJmZZUhEjAHDFfs/BH7YuojMzAz8L0Qzs1xSUa+kTyT7p0u6sNVxmZm1OyfXZmb59JfARcB7k/0XgdtaF46ZmYGHhZhZhdJUgR6vnQtvj4i3SHoKICL+WZJnCzEzazH3XJs1wdDQUPnxsmXLxu2b1emYpOOAAJC0GPhFa0MyMzMn12YNNjQ0RF9fX3l/37599PX1OcG2ufoL4F7gdZIGgMeAP2ltSGZm5uTaLGX9/f3l1Rgl0dvby9jY2LhzxsbG6O3tLZ/T39/fmmBTUm31SUtfRAwBfwh8CjgAXIUXkTEzazkn12Yp6urqSr2slzm3SpJOlvQxSZ8HTqf4wcZXASMU5742M7MWcnJtlqK1a9eWk+HSVigUpjy3UCiMO2/t2rVNjtZy6uvA2cD3gPcDfwu8B1gdEavnWrmkX5b0LUm7k6+vnea8lyVtS7b753pdM7P5wsm1WYMNDAzQ0dEx7lhHRwcDAwMtishybnlEXBcRX6I4DV83sCoitqVU/83AQxHRBTyU7E/lJxFxfrK9K6Vrm5nlnpNrswbr6elhcHCwvF8oFBgcHKSnp6eFUVmOHSs9iIiXgX+MiBdTrH818LXk8deA/5hi3WZm856Ta7MmqEyk9+7dm8nE2tMF5sabJR1JtheB80qPJR1Jof7XJ0upl5ZUf900550gaauk/yPJCbiZWcKLyJi1sXXr1gGwfft2RkZGysf37dvH9ddfz/DwMOedd175eFdXl8eGt1hEHDfXOiRtBt4wxVMfn0U1p0fEAUnLgW9L+l5E/MMU1+oD+gBOP/30uuI1M8sTJ9dmbairq6tq7/mxY8cYHh5meHi4fCzvUwZaUURcOt1zkp6X9MaI+KGkNwI/mqaOA8nXPZL+DrgAmJRcR8QgMAjQ3d3taW/MbN7zsBCzNuTeZ5vB/cD7ksfvA+6beIKk10p6dfL4VOAdwPebFqGZWYZlJrmWtFzSVyTd0+pYzNrBbKcLnK7X2mO1551PA5dJ2g1cluwjqVvSl5Nzfg3YKulp4GHg0xHh5NrMjAYPC5H0VWAV8KOIOLfi+Ergs8BxwJcj4tMRsQe4wcm1WfMNDAzQ19c3biXJ6aYLLI3TBo/Vno8i4hCwYorjWynOq01EPA78epNDMzPLhUaPud4AfB7469IBSccBt1HsERkFtki6370eZq1TGn/d29sLFHusBwYGxo3LrmWcNnistpmZtbeGDguJiEeBH084fCHwbETsiYiXgLsozqtaE0l9yfRPWw8ePJhitGbtrdp0ge55NjMzq64VY66XAPsr9keBJZIWSfoicIGkj01XOCIGI6I7IroXL17c6FjNrEI9S7u719rMzNpJK5JrTXEsIuJQRNwYEf8mIj7V9KjMbNbmurS7JKSpmoTGaPb1zMys/bQiuR4FTqvYXwocaEEcZk1V6smdT7y0u5mZ2XitSK63AF2SzpC0ELiW4ryqZpZDeVja3czMrFkamlxLuhN4Ajhb0qikGyLi58CHgAeBZ4CNEbFzlvVeIWnw8OHD6QdtZmZmZlanhk7FFxHvneb4A8ADc6h3BBjp7u7+QL11mJmZmZmlLTMrNJqZmZmZ5Z2TazPLHc/6YWZmWeXk2szMzMwsJblMrv2BRjMzMzPLolwm1xExEhF9nZ2drQ7FzMzMzKysobOFmFm+zLdFbszMzJotlz3XZpZ/Q0ND5cfLli0btz8frmdmNhvzcRXfduXk2syabmhoiL6+vvL+vn376Ovra1jC2+zrmZlZ+3JybWYN19/fX54+TxK9vb2MjY2NO2dsbIze3t5x59Wrso5mXM/ag3sWzawWuUyuPVuImZmZmWVRLpNrzxZili3VevT6+/vL50QEhUJhyvMKhcK48+YaT7OuZ2ZmVpLL5NrM8m1gYICOjo5xxzo6OhgYGJgX17Ns8x9TZtZITq7NrOl6enoYHBws7xcKBQYHB+np6ZkX1zMzs/bl5NrMWqIysd27d2/NiW69U+rVez0zyw7/18HywMm1meWGp9QzM7Os8wqNZpZZtUyPV5pSr7e3t3zMPVtmZtnRbm1yLnuuPRWf2fx3xx13tKTsVNppHmxJV0naKekXkrpnOG+lpF2SnpV0czNjNDPLslz2XEfECDDS3d39gVbHYmaNsXv3bvr7+8cdu/XWW5nqj+rOzk5uuummcWWtbjuAdwNfmu4ESccBtwGXAaPAFkn3R8T3mxOiWe3q7TVtt95WS08uk2szax+33HJL+fGZZ55JX1/fuNUWOzo6uO2228ofUFy3bl3TY5xPIuIZqDok50Lg2YjYk5x7F7AacHJtZm0vl8NCzKw9eUq9zFgC7K/YH02OTSKpT9JWSVsPHjyYahCeOcKsuZr9nqv3eq1uG5xcm1mueEq9uZO0WdKOKbbVtVYxxbEpf5NFxGBEdEdE9+LFi+sP2swsJzwsxMyszUTEpXOsYhQ4rWJ/KXBgjnWamc0L7rk2M7PZ2gJ0STpD0kLgWuD+FsdkZpYJuUyuPRWfmVljSLpS0ihwEfBNSQ8mx39F0gMAEfFz4EPAg8AzwMaI2NmqmC1/Wj0m1vwaNFIuh4V4Kj4zs8aIiHuBe6c4fgD4nYr9B4AHmhiaZZCTs9bza5A9uUyuzczM8qLe2Q6adS0zS5eTazMzszbnpNwsPbkcc21mZmZmlkVOrs3MzMzMUuJhIWbWMs3+V7T/9W1mZo3mnmszsxkMDQ2VHy9btmzcvpmZ2UROrs3MpjE0NERfX195f9++ffT19TnBNjOzaXlYiJlZQlLVc8bGxujt7aW3t7d8zMNNzMysJJc9116h0czMzLLIKx9aLpPriBiJiL7Ozs5Wh2Jm80jpl2JpKxQKU55XKBTGnWdmZlaSy+TazKwZBgYG6OjoGHeso6ODgYGBFkVkZmZZ5+TazGwaPT09DA4OlvcLhQKDg4P09PS0MCozM8syJ9dmZjOoTKT37t3rxNrMzGbk2ULMLHc8ztnMzLLKPddmZmZmZilxcm1mZmZmlhLl+d+rkg4C+6Z5uhOodSLsWs+tdl61508FXqgxpiybzc82y9dNo7566vC92Th5ujcLEbG4EcFkVQbb7Grn+H2RvevOtc68tNnVzvG92Zrr1tZuT5zXdb5swGDa51Y7r4bnt7b659Lsn22Wr5tGffXU4XszP/dI1q87n7ZWvC+qneP3RfauO9c689JmVzvH92a2rzufh4WMNODcaufN5pp51qrvM+3rplFfPXX43myc+XJvtqNWvC9me928mk/vi7nWmZc2e7bXzav5dG+W5XpYSN5I2hoR3a2Ow2wi35tmk/l9YVnlezPb5nPPdRYNVj/FrCV8b5pN5veFZZXvzQxzz7WZmZmZWUrcc21mZmZmlhIn12ZmZmZmKXFybWZmZmaWEifXZmZmZmYpcXLdQpJOlPQ1SX8lqafV8ZiVSFou6SuS7ml1LGZZ4TbbssptdrY4uU6ZpK9K+pGkHROOr5S0S9Kzkm5ODr8buCciPgC8q+nBWluZzb0ZEXsi4obWRGrWPG6zLavcZueXk+v0bQBWVh6QdBxwG/DbwJuA90p6E7AU2J+c9nITY7T2tIHa702zdrEBt9mWTRtwm51LTq5TFhGPAj+ecPhC4NnkL8uXgLuA1cAoxcYa/FpYg83y3jRrC26zLavcZueXG4fmWMIrvR1QbKCXAMPAGklfoMHr3JtNY8p7U9IiSV8ELpD0sdaEZtYybrMtq9xm58DxrQ6gTWiKYxERR4Hrmx2MWYXp7s1DwI3NDsYsI9xmW1a5zc4B91w3xyhwWsX+UuBAi2Ixq+R702wyvy8sq3xv5oCT6+bYAnRJOkPSQuBa4P4Wx2QGvjfNpuL3hWWV780ccHKdMkl3Ak8AZ0salXRDRPwc+BDwIPAMsDEidrYyTms/vjfNJvP7wrLK92Z+KSJaHYOZmZmZ2bzgnmszMzMzs5Q4uTYzMzMzS4mTazMzMzOzlDi5NjMzMzNLiZNrMzMzM7OUOLk2MzMzM0uJk2vLNEkvS9pWsd3c6phKJN0jabmk7ySx/ZOkgxWxLpum3B9L+uSEY92StiePH5LU2fjvwMwsXW6zzTzPtWWcpH+NiF9Kuc7jk4n451LHOcAfR8SVFceuA7oj4kM1lL03Is6qOLYeOBQRn5J0A3BqRPzpXGI0M2s2t9lm7rm2nJK0V9I6Sf9X0vck/Wpy/ERJX5W0RdJTklYnx6+TdLekEeBvJb1K0l9K2inpG5IekPQeSSsk3VtxncskDU8RQg9wXw1x/rakJ5I4/0bSiclqWj+V9NbkHAFXAXclxe4D1s7l52NmliVus62dOLm2rHvNhH8xXlPx3AsR8RbgC8BHk2MfB74dEW8Dfgv4jKQTk+cuAt4XEf8OeDewDPh14P3JcwDfBn5N0uJk/3rg9iniegfw5EyBS3odcDOwIolzO/D7ydN3AtdW1HUgIv4RICJeAE6SdMpM9ZuZZZDbbGt7x7c6ALMqfhIR50/zXKl34kmKDS/A5cC7JJUa7hOA05PH34qIHyeP3wncHRG/AJ6T9DBARISkrwO9km6n2ID/3hTXfiNwsErsvwG8CXi82NHBQuCx5Lk7gUck/SHFBvvOCWUPJtf4lyrXMDPLErfZ1vacXFue/Sz5+jKv3MsC1kTErsoTJb0dOFp5aIZ6bwdGgJ9SbMynGuv3E4q/BGYiYFNE/O7EJyJir6QDwL8FrgTeOuGUE5JrmJnNF26zrS14WIjNNw8C/yUZE4ekC6Y57zFgTTKO7/XAb5aeiIgDwAHgj4AN05R/BjizSiyPA5dIWp7EcqKkrorn7wT+AngmIp4rHZT0KuBUYH+V+s3M8s5tts07Tq4t6yaO3/t0lfM/CSwAtkvakexP5X8Co8AO4EvAd4DDFc8PAfsj4vvTlP8mFY37VCLieeAG4G8kPU2x4T6r4pSNwLm88qGYkguBxyLi5ZnqNzPLILfZ1vY8FZ+1LUm/FBH/KmkR8F3gHaXeCEmfB56KiK9MU/Y1wMNJmVQbVEm3ARsj4pE06zUzyzO32ZYXHnNt7ewbyae7FwKfrGikn6Q41u8j0xWMiJ9IugVYAvxTynE95UbazGwSt9mWC+65NjMzMzNLicdcm5mZmZmlxMm1mZmZmVlKnFybmZmZmaXEybWZmZmZWUqcXJuZmZmZpeT/A+mrjRIN5nl6AAAAAElFTkSuQmCC\n",
  250 + "text/plain": [
  251 + "<matplotlib.figure.Figure at 0x110a33908>"
  252 + ]
  253 + },
  254 + "metadata": {},
  255 + "output_type": "display_data"
  256 + }
  257 + ],
  258 + "source": [
  259 + "residuals = 'residuals_classical.fits'\n",
  260 + "res = cscripts.csresspec(like.obs())\n",
  261 + "\n",
  262 + "res['statistic'] = 'WSTAT'\n",
  263 + "res['components'] = True\n",
  264 + "res['algorithm'] = 'SIGNIFICANCE'\n",
  265 + "res['outfile'] = residuals\n",
  266 + "\n",
  267 + "res.execute()\n",
  268 + "\n",
  269 + "from show_residuals import plot_residuals\n",
  270 + "plot_residuals(residuals,'',0)"
  271 + ]
  272 + }
  273 + ],
  274 + "metadata": {
  275 + "kernelspec": {
  276 + "display_name": "Python 3",
  277 + "language": "python",
  278 + "name": "python3"
  279 + },
  280 + "language_info": {
  281 + "codemirror_mode": {
  282 + "name": "ipython",
  283 + "version": 3
  284 + },
  285 + "file_extension": ".py",
  286 + "mimetype": "text/x-python",
  287 + "name": "python",
  288 + "nbconvert_exporter": "python",
  289 + "pygments_lexer": "ipython3",
  290 + "version": "3.5.4"
  291 + }
  292 + },
  293 + "nbformat": 4,
  294 + "nbformat_minor": 2
  295 +}
... ...
doc/source/users/tutorials/hess_dr1/classical.rst deleted
View file @b70128d
1   -.. _hess_dr1_classical:
2   -
3   -Performing a classical On/Off analysis
4   ---------------------------------------
5   -
6   - .. admonition:: What you will learn
7   -
8   - You will learn **how to perform a classical On/Off analysis of the
9   - data**.
10   -
11   -
12   -.. code-block:: bash
13   -
14   - $ csphagen use_irf_bkg=no
15   - Input event list or observation definition XML file [obs.xml] obs_crab_selected.xml
16   - Input model definition XML file (if NONE, use point source) [NONE] crab_models.xml
17   - Source name [Crab]
18   - Algorithm for defining energy bins (FILE|LIN|LOG) [LOG]
19   - Start value for first energy bin in TeV [0.1] 0.66
20   - Stop value for last energy bin in TeV [100.0]
21   - Number of energy bins [120] 40
22   - Stack multiple observations into single PHA, ARF and RMF files? [no]
23   - Output observation definition XML file [onoff_obs.xml] onoff_obs_crab.xml
24   - Method for background estimation (REFLECTED|CUSTOM) [REFLECTED]
25   - Coordinate system (CEL - celestial, GAL - galactic) (CEL|GAL) [CEL]
26   - Right Ascension of source region centre (deg) (0-360) [83.63]
27   - Declination of source region centre (deg) (-90-90) [22.01]
28   - Radius of source region circle (deg) (0-180) [0.2]
29   -
30   -.. code-block:: xml
31   -
32   - <?xml version="1.0" encoding="UTF-8" standalone="no"?>
33   - <source_library title="source library">
34   - <source name="Crab" type="PointSource">
35   - <spectrum type="PowerLaw">
36   - <parameter name="Prefactor" value="1" scale="5.7e-18" min="0" free="1" />
37   - <parameter name="Index" value="1" scale="-2.48" min="-4" max="4" free="1" />
38   - <parameter name="PivotEnergy" value="1" scale="1000000" free="0" />
39   - </spectrum>
40   - <spatialModel type="PointSource">
41   - <parameter name="RA" value="83.633" scale="1" free="0" />
42   - <parameter name="DEC" value="22.015" scale="1" free="0" />
43   - </spatialModel>
44   - </source>
45   - </source_library>
46   -
47   -.. code-block:: bash
48   -
49   - $ ctlike statistic=wstat
50   - Input event list, counts cube or observation definition XML file [events.fits] onoff_obs_crab.xml
51   - Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_models_onoff.xml
52   - Output model definition XML file [crab_results.xml]
53   -
54   -.. code-block:: none
55   -
56   - 2018-09-13T11:49:06: +=========================================+
57   - 2018-09-13T11:49:06: | Maximum likelihood optimisation results |
58   - 2018-09-13T11:49:06: +=========================================+
59   - 2018-09-13T11:49:06: === GOptimizerLM ===
60   - 2018-09-13T11:49:06: Optimized function value ..: 60.048
61   - 2018-09-13T11:49:06: Absolute precision ........: 0.005
62   - 2018-09-13T11:49:06: Acceptable value decrease .: 2
63   - 2018-09-13T11:49:06: Optimization status .......: converged
64   - 2018-09-13T11:49:06: Number of parameters ......: 6
65   - 2018-09-13T11:49:06: Number of free parameters .: 2
66   - 2018-09-13T11:49:06: Number of iterations ......: 5
67   - 2018-09-13T11:49:06: Lambda ....................: 1e-08
68   - 2018-09-13T11:49:06: Maximum log likelihood ....: -60.048
69   - 2018-09-13T11:49:06: Observed events (Nobs) ...: 605.000
70   - 2018-09-13T11:49:06: Predicted events (Npred) ..: 604.191 (Nobs - Npred = 0.808986318662733)
71   - 2018-09-13T11:49:06: === GModels ===
72   - 2018-09-13T11:49:06: Number of models ..........: 1
73   - 2018-09-13T11:49:06: Number of parameters ......: 6
74   - 2018-09-13T11:49:06: === GModelSky ===
75   - 2018-09-13T11:49:06: Name ......................: Crab
76   - 2018-09-13T11:49:06: Instruments ...............: all
77   - 2018-09-13T11:49:06: Instrument scale factors ..: unity
78   - 2018-09-13T11:49:06: Observation identifiers ...: all
79   - 2018-09-13T11:49:06: Model type ................: PointSource
80   - 2018-09-13T11:49:06: Model components ..........: "PointSource" * "PowerLaw" * "Constant"
81   - 2018-09-13T11:49:06: Number of parameters ......: 6
82   - 2018-09-13T11:49:06: Number of spatial par's ...: 2
83   - 2018-09-13T11:49:06: RA .......................: 83.633 deg (fixed,scale=1)
84   - 2018-09-13T11:49:06: DEC ......................: 22.015 deg (fixed,scale=1)
85   - 2018-09-13T11:49:06: Number of spectral par's ..: 3
86   - 2018-09-13T11:49:06: Prefactor ................: 4.77687464316323e-17 +/- 3.12433509735449e-18 [0,infty[ ph/cm2/s/MeV (free,scale=5.7e-18,gradient)
87   - 2018-09-13T11:49:06: Index ....................: -2.60709757251519 +/- 0.0767624316833204 [9.92,-9.92] (free,scale=-2.48,gradient)
88   - 2018-09-13T11:49:06: PivotEnergy ..............: 1000000 MeV (fixed,scale=1000000,gradient)
89   - 2018-09-13T11:49:06: Number of temporal par's ..: 1
90   - 2018-09-13T11:49:06: Normalization ............: 1 (relative value) (fixed,scale=1,gradient)
91   -
92   -The intensity at 1 TeV of the Crab was fitted to
93   -:math:`(4.78 \pm 0.31) \times 10^{-11}\,{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm TeV}^{-1}`
94   -the spectral index of the power law is
95   -:math:`-2.61 \pm 0.08`.
96   -This can be compared to the values
97   -:math:`(4.64 \pm 0.26) \times 10^{-11}\,{\rm photons}\,{\rm cm}^{-2}\,{\rm s}^{-1}\,{\rm TeV}^{-1}`
98   -and
99   -:math:`-2.61 \pm 0.07`
100   -that were obtained before with the unbinned maximum likelihood fit using the
101   -energy dispersion.
102   -
103   -.. code-block:: bash
104   -
105   - $ csresspec components=yes
106   - Input event list, counts cube, or observation definition XML file [events.fits] onoff_obs_crab.xml
107   - Stack observations? [no] yes
108   - Input model definition XML file [$CTOOLS/share/models/crab.xml] crab_results_onoff.xml
109   - Residuals computation algorithm (SUB|SUBDIV|SUBDIVSQRT|SIGNIFICANCE) [SIGNIFICANCE]
110   - Output residual spectrum file [resspec.fits] resspec_onoff.fits
111   -
112   -The script will produce the FITS file ``resspec.fits`` that contains the
113   -spectral residuals for each of the model components.
114   -To display the residuals you can use the ``show_residuals.py`` script in the
115   -example folder as follows:
116   -
117   -.. code-block:: bash
118   -
119   - $CTOOLS/share/examples/python/show_residuals.py resspec.fits
120   -
121   -The figure below shows the resulting residual spectrum.
122   -Although the model provides a reasonable first order approximation of the
123   -spectral distribution of the data, there are still significant spectral
124   -residuals that cannot be explained by the model.
125   -
126   -.. figure:: classical_resspec.png
127   - :width: 600px
128   - :align: center
129   -
130   - *Residual counts spectrum after model fitting*
doc/source/users/tutorials/hess_dr1/selecting_events.ipynb
  View file @dad5f0e
... ... @@ -11,7 +11,7 @@
11 11 "cell_type": "raw",
12 12 "metadata": {},
13 13 "source": [
14   - "<div align=\"right\"><a href=\"../../../_downloads/howto_advanced_python_fitting.ipynb\" download=\"selecting_events_HESSDR1.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
  14 + "<div align=\"right\"><a href=\"../../../_downloads/selecting_events.ipynb\" download=\"selecting_events.ipynb\"><img src=\"../../../_static/download-notebook.jpg\" alt=\"Download Notebook\" height=\"40\"></a></div>"
15 15 ]
16 16 },
17 17 {
... ...
doc/source/users/tutorials/notebooks.rst
  View file @dad5f0e
... ... @@ -20,4 +20,6 @@ this page.
20 20 * :download:`How to perform a stacked analysis? <howto/howto_stacked_analysis.ipynb>`
21 21 * :download:`How to do advanced model manipulation and fitting in Python? <howto/howto_advanced_python_fitting.ipynb>`
22 22 * :download:`Data preparation for 1DC How-tos <howto/howto_prepare_1dc.ipynb>`
23   -* :download:`Selecting events for H.E.S.S. public data <hess_dr1/selecting_events.ipynb>`
  23 +* :download:`Selecting events for H.E.S.S. public data <hess_dr1/selecting_events.ipynb>`
  24 +* :download:`Classical analysis of H.E.S.S. public data <hess_dr1/classical.ipynb>`
  25 +
... ...
models/Makefile.am
  View file @dad5f0e
... ... @@ -31,6 +31,7 @@ models = 1dc_howto.xml \
31 31 crab.xml \
32 32 crab_beyond.xml \
33 33 crab_2bkg.xml \
  34 + crab_nobkg.xml \
34 35 crab_poly.xml \
35 36 crab_profile.xml \
36 37 disk.xml \
... ...
models/crab_nobkg.xml 0 → 100755
  View file @dad5f0e
  1 +<?xml version="1.0" standalone="no"?>
  2 +<source_library title="source library">
  3 + <source name="Crab" type="PointSource">
  4 + <spectrum type="PowerLaw">
  5 + <parameter name="Prefactor" scale="1e-17" value="5.7" min="1e-07" max="1000.0" free="1"/>
  6 + <parameter name="Index" scale="-1" value="2.48" min="0.0" max="+5.0" free="1"/>
  7 + <parameter name="PivotEnergy" scale="1e6" value="1." min="0.01" max="1000.0" free="0"/>
  8 + </spectrum>
  9 + <spatialModel type="PointSource">
  10 + <parameter name="RA" scale="1.0" value="83.6331" min="-360" max="360" free="0"/>
  11 + <parameter name="DEC" scale="1.0" value="22.0145" min="-90" max="90" free="0"/>
  12 + </spatialModel>
  13 + </source>
  14 +</source_library>
... ...