{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "499e39c3",
   "metadata": {},
   "source": [
    "# Signal identification and background estimation"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09090d8d",
   "metadata": {},
   "source": [
    "Signal identification is one of the first step of the analysis workflow. It consist on selecting the time period where the signal occurred, as well as other time periods before (and possibbly after) that will be use to compute the background.\n",
    "\n",
    "In bc-tool the signal identification and background calculation is performed itertively using [bayesian blocks](https://arxiv.org/abs/1304.2818). Each iteration fits the background with a polynomial, identifies peaks, and then removes them from the data use to fit the background in the next iteration. This continue until it converges. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9510dd26",
   "metadata": {},
   "source": [
    "### Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d51e5493",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gdt.missions.fermi.gbm.tte import GbmTte\n",
    "from gdt.core.binning.unbinned import bin_by_time\n",
    "from gdt.core.binning.binned import rebin_by_time\n",
    "\n",
    "from bctools.analysis import BayesianBlocksLightcurve\n",
    "\n",
    "import numpy as np\n",
    "import urllib"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69635f8e",
   "metadata": {},
   "source": [
    "### Get the light curve"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d8b11f5",
   "metadata": {},
   "source": [
    "For this example will use real GBM's data. We'll start with event data (TTE) and binned it into a light curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e517f2a7-fe1e-4a7b-9df0-08703ec733bb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Long GRB example\n",
    "filename,_ = urllib.request.urlretrieve(\"https://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/triggers/2021/bn210925800/current/glg_tte_n6_bn210925800_v00.fit\")\n",
    "\n",
    "# Short GRB example\n",
    "#filename,_ = urllib.request.urlretrieve(\"https://heasarc.gsfc.nasa.gov/FTP/fermi/data/gbm/triggers/2009/bn090120627/current/glg_tte_n5_bn090120627_v01.fit\")\n",
    "\n",
    "tte = GbmTte.open(filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "000ba72c-c3c5-450c-ae15-755cde28281f",
   "metadata": {},
   "source": [
    "Now let's select and energy range and bin the data in time to create a light curve. The algorithm correctly handles Poisson statistics, so you can bin the data very finely if you wish. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4895547e-4afb-4345-bf07-dfefa1f290a6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bin in time (32ms)\n",
    "phaii = tte.to_phaii(bin_by_time, 0.032, time_ref=0.0, energy_range = (50.0, 300.0))\n",
    "\n",
    "lc = phaii.to_lightcurve()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0376ead",
   "metadata": {},
   "source": [
    "The ``BayesianLightcurve`` class performs the Bayesian blocks calculations, as well as the background fitting and signal identification:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "0dc1060c",
   "metadata": {},
   "outputs": [],
   "source": [
    "bb_lc = BayesianBlocksLightcurve(lc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14b43ad2",
   "metadata": {},
   "source": [
    "### Identify signal"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82d3d2e7",
   "metadata": {},
   "source": [
    "The meat of the calculation is performed by ``compute_bayesian_blocks``. It iteratively fits the background, identifies the signal, and bins the light curve in Bayesian blocks. You can call it directly and the results will be cache for the other methods, or alternatively the other methods will call it for you with default parameters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "088a8807",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10.3 s, sys: 864 ms, total: 11.1 s\n",
      "Wall time: 8.76 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "bb_lc.compute_bayesian_blocks()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cfc0f8e7",
   "metadata": {},
   "source": [
    "Now we can get the start and stop times of the signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "35b58996-2412-49bc-8b78-4cef818189ea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<TimeRange: (-6.944, 236.352)>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bb_lc.signal_range"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "524396cf",
   "metadata": {},
   "source": [
    "### Plot"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b13627c2-d89b-4269-9092-fef5cfb9e0d7",
   "metadata": {},
   "source": [
    "You can use the `plot()` method for a quick view of the light curve. If you want to fine tune your plot, you can always obtain the underlaying data with the properties `lightcurve`, `bb_lightcurve` and `bkg_model`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "e32e4026-a8aa-4480-ad98-cbc87fee71c7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Axes: >"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABygklEQVR4nO3deVxU5f4H8M8ZZmFHRWVRRExxA1dywVJJVFQ009SyW1pm11yulraYldiGLW7ZtW7mT80lytIWTHJJNFNLURSXzBR3CDUWWWc7vz+QI8MiDMwwc4bP+/WaF8NzzpzzzOFhzneeVRBFUQQRERGRTChsnQEiIiIiczB4ISIiIllh8EJERESywuCFiIiIZIXBCxEREckKgxciIiKSFQYvREREJCsMXoiIiEhWlLbOQE0YjUZcu3YNHh4eEATB1tkhIiKiahBFEbdu3YK/vz8UiprXn8gyeLl27RoCAgJsnQ0iIiKqgcuXL6N58+Y1fr0sgxcPDw8AxW/e09PTxrkhS9Bq87BokT8AYPbsa1Cr3WycI6KaYVkmqlxOTg4CAgKk+3hNyTJ4KWkq8vT0ZPDiILRaJzg7Fz/39PTkBz7JFssyUdVq2+WDHXaJiIhIVmRZ80KOR6FQonPnCdJzIrliWSayPkEURdHWmTBXTk4OvLy8kJ2dzWYjIiIimbDU/ZtfC4iI7JQoitDr9TAYDLbOClG1qVQqODk5WfUcDF7ILoiiCJ0uHwCgUrly/h6SLUuVZa1Wi7S0NOTn51sye0RWJwgCmjdvDnd3d6udg8EL2QWdLh+xscUFfe7cXI7QINmyRFk2Go1ITU2Fk5MT/P39oVarGdCTLIiiiOvXr+PKlSto06aN1WpgGLwQEdkZrVYLo9GIgIAAuLq62jo7RGZp0qQJLly4AJ1OZ7XghUOliYjsVG2mTyeylbqoJeR/BhEREckKgxciIqp3Jk6ciJEjR1r9PIIg4Ntvv610+4ULFyAIApKTky1yPksfz14xeCEiIouZOHEiBEGQHt7e3oiKisLx48dtnTUTy5Ytw5o1a2ydDaohBi9ERGRRUVFRSEtLQ1paGnbt2gWlUono6GhbZ8uEl5cXGjRoYOtsUA2ZFbzExMSYRNSCIMDX11faXjbiFgQBvXr1MjlGUVERZsyYgcaNG8PNzQ0jRozAlStXLPNuSLYUCid06PAwOnR4GAqFdSc3IrImlmVAo9HA19cXvr6+6NKlC1566SVcvnwZ169fl/Z56aWXEBwcDFdXV7Rq1QqvvfYadDodgOKmD4VCgcOHD5scd/ny5QgMDETJxPCnTp3C0KFD4e7uDh8fHzz++OO4ceOGtP/XX3+N0NBQuLi4wNvbG5GRkcjLywNQvtkoISEB9913Hxo0aABvb29ER0fj3Llz0vaS5pjNmzcjIiICrq6u6Ny5Mw4cOFDl9UhLS8OQIUPg4uKCoKAgbNq06a7779mzBz169IBGo4Gfnx9efvll6PV6abvRaMS7776L1q1bQ6PRoEWLFnj77bcrPJbRaMTkyZMRHByMixcvAii+l7do0QIajQb+/v74z3/+U+V7sDdm17x07NhRiqjT0tKQkpJisr10xJ2WloYff/zRZPusWbOwZcsWxMXFYd++fcjNzUV0dDRnkKznlEpnjBmzCWPGbIJS6Wzr7FA9o9frER8fj/j4eJObRE1YqyyLooh8rd4mj9qsIpObm4sNGzagdevW8Pb2ltI9PDywZs0anDp1CsuWLcPKlSuxZMkSAEDLli0RGRmJ1atXmxxr9erV0pfktLQ09OvXD126dMHhw4eRkJCAv//+G2PHjgVQHDA8+uijeOqpp3D69GkkJiZi1KhRlb6XvLw8PP/88zh06BB27doFhUKBhx56CEaj0WS/efPmYc6cOUhOTkZwcDAeffTRKsvMa6+9htGjR+PYsWP417/+hUcffRSnT5+ucN+rV69i6NChuPfee3Hs2DF8/PHHWLVqFd566y1pn7lz5+Ldd9/Fa6+9hlOnTmHjxo3w8fEpdyytVouxY8fi8OHD2LdvHwIDA/H1119jyZIl+N///oezZ8/i22+/RWho6F3zb4/MnudFqVSa1LaUVRJxVyQ7OxurVq3CunXrEBkZCQBYv349AgICsHPnTgwePNjc7BAR1QsFOgM6vP6TTc596o3BcFVX/3YRHx8vza6al5cHPz8/xMfHmwz9fvXVV6XnLVu2xOzZs/Hll1/ixRdfBAA8/fTTmDJlChYvXgyNRoNjx44hOTkZmzdvBgB8/PHH6NatG9555x3pOP/3f/+HgIAA/Pnnn8jNzYVer8eoUaMQGBgIAHe9SY8ePdrk91WrVqFp06Y4deoUQkJCpPQ5c+Zg2LBhAIAFCxagY8eO+Ouvv9CuXbtKjz1mzBg8/fTTAIA333wTO3bswPLly7FixYpy+65YsQIBAQH46KOPIAgC2rVrh2vXruGll17C66+/jry8PCxbtgwfffQRJkwoXgD0nnvuwX333WdynNzcXAwbNgwFBQVITEyEl5cXAODSpUvw9fVFZGQkVCoVWrRogR49elSad3tlds3L2bNn4e/vj6CgIDzyyCM4f/68yfbExEQ0bdoUwcHBmDx5MjIyMqRtSUlJ0Ol0GDRokJTm7++PkJAQ7N+/v9JzFhUVIScnx+RBRET2KSIiAsnJyUhOTsZvv/2GQYMGYciQIVKzBVDcpHPffffB19cX7u7ueO2113Dp0iVp+8iRI6FUKrFlyxYAxYFJREQEWrZsCaD4frJ79264u7tLj5IA4ty5c+jcuTMGDBiA0NBQjBkzBitXrkRmZmaleT537hzGjx+PVq1awdPTE0FBQQBgkicA6NSpk/Tcz88PAEzucxXp3bt3ud8rq3k5ffo0evfubTJXSp8+fZCbm4srV67g9OnTKCoqwoABA+56zkcffRS5ubnYvn27FLgAxYFUQUEBWrVqhcmTJ2PLli21rm20BbNqXnr27InPP/8cwcHB+Pvvv/HWW28hPDwcJ0+ehLe3N4YMGYIxY8YgMDAQqampeO211/DAAw8gKSkJGo0G6enpUKvVaNiwoclxfXx8kJ6eXul5Y2NjsWDBgpq9Q5IFrTaPywOQQ7BWWXZROeHUG7apnXZRmdd3x83NDa1bt5Z+7969O7y8vLBy5Uq89dZbOHjwIB555BEsWLAAgwcPhpeXF+Li4rBo0SLpNWq1Go8//jhWr16NUaNGYePGjVi6dKm03Wg0Yvjw4Xj33XfLnd/Pzw9OTk7YsWMH9u/fj+3bt2P58uWYN28efvvtNykwKW348OEICAjAypUr4e/vD6PRiJCQEGi1WpP9VCqV9LwkwCjbtFQdlU3kJopiuW0lTV2CIMDFxaVaxx86dCjWr1+PgwcP4oEHHpDSAwICcObMGezYsQM7d+7E1KlT8f7772PPnj0m783emRW8DBkyRHoeGhqK3r1745577sHatWvx/PPPY9y4cdL2kJAQhIWFITAwEFu3bsWoUaMqPW5Ff6zS5s6di+eff176PScnBwEBAeZknYhI1gRBMKvpxp4IggCFQoGCggIAwK+//orAwEDMmzdP2qd0rUyJp59+GiEhIVixYgV0Op3JfaRbt2745ptv0LJlSyiVFV8XQRDQp08f9OnTB6+//joCAwOxZcsWk/sJANy8eROnT5/G//73P9x///0AgH379tX6fZc4ePAgnnjiCZPfu3btWuG+HTp0wDfffGNyX9y/fz88PDzQrFkzNGnSBC4uLti1a5fUFFWRZ599FiEhIRgxYgS2bt2Kfv36SdtcXFwwYsQIjBgxAtOmTUO7du2QkpKCbt26WegdW1+t/hPc3NwQGhqKs2fPVrjdz88PgYGB0nZfX19otVpkZmaa1L5kZGQgPDy80vNoNBpoNJraZJWIiOpIUVGRVJuemZmJjz76CLm5uRg+fDgAoHXr1rh06RLi4uJw7733YuvWrVLzUGnt27dHr1698NJLL+Gpp54yqXWYNm0aVq5ciUcffRQvvPACGjdujL/++gtxcXFYuXIlDh8+jF27dmHQoEFo2rQpfvvtN1y/fh3t27cvd56GDRvC29sbn376Kfz8/HDp0iW8/PLLFrsemzZtQlhYGO677z5s2LABv//+O1atWlXhvlOnTsXSpUsxY8YMTJ8+HWfOnMH8+fPx/PPPQ6FQwNnZGS+99BJefPFFqNVq9OnTB9evX8fJkycxadIkk2PNmDEDBoMB0dHR2LZtG+677z6sWbMGBoMBPXv2hKurK9atWwcXFxepX5Bc1Gqel6KiIpw+fVpq9yvr5s2buHz5srS9e/fuUKlU2LFjh7RPWloaTpw4cdfghYiI5CMhIQF+fn7w8/NDz549cejQIWzatAn9+/cHADz44IN47rnnMH36dHTp0gX79+/Ha6+9VuGxJk2aBK1Wi6eeesok3d/fH7/++isMBgMGDx6MkJAQzJw5E15eXlAoFPD09MTevXsxdOhQBAcH49VXX8WiRYtMWhBKKBQKxMXFISkpCSEhIXjuuefw/vvvW+x6LFiwAHFxcejUqRPWrl2LDRs2oEOHDhXu26xZM/z444/4/fff0blzZ0yZMgWTJk0y6eD82muvYfbs2Xj99dfRvn17jBs3rtJ+N7NmzcKCBQswdOhQ7N+/Hw0aNMDKlSvRp08fdOrUCbt27cIPP/xgMhJMDgTRjDFwc+bMwfDhw9GiRQtkZGTgrbfewp49e5CSkgJvb2/ExMRg9OjR8PPzw4ULF/DKK6/g0qVLOH36NDw8PAAUV2XFx8djzZo1aNSoEebMmYObN28iKSmp2qtP5uTkwMvLC9nZ2fD09KzZOye7wj4vZEt6vR4JCQkAiqd7qKwZojosUZYLCwuRmpqKoKAgODvX76kD3n77bcTFxZWbloPs193Kr6Xu32b9h165cgWPPvoobty4gSZNmqBXr144ePAgAgMDUVBQgJSUFHz++efIysqCn58fIiIi8OWXX0qBCwAsWbIESqUSY8eORUFBAQYMGIA1a9ZYbdlsIiKSn9zcXJw+fRrLly/Hm2++aevskJ0xK3iJi4urdJuLiwt++qnqOQicnZ2xfPlyLF++3JxTExFRPTJ9+nR88cUXGDlyZLkmIyJ5dl0nh6NQOKFNm6HScyK5Ylm2jDVr1nDhRKoUgxeyC0qlM8aP32rrbBDVGssykfVxVWkiIiKSFQYvREREJCtsNiK7oNXm4YMPmgIA5szJ4FBpki2WZSLrY/BCdkOny7d1FogsgmWZyLrYbERERESywuCFiIjqRP/+/TFr1qw6O9+aNWvQoEGDSrdfuHABgiAgOTnZqvkQBAHffvutVc9R12JiYtClSxebnZ/BCxERWczEiRMhCEK5x19//YXNmzebzJbbsmVLLF261OT1VQUcRAD7vBARkYVFRUVh9erVJmlNmjThMjAWpNPpoFKpbJ0Nm2HNCxERWZRGo4Gvr6/Jw8nJyaTZqH///rh48SKee+45qXYmMTERTz75JLKzs6W0mJgYAIBWq8WLL76IZs2awc3NDT179kRiYqLJedesWYMWLVrA1dUVDz30EG7evFmt/P7xxx8IDw+Hs7MzOnbsaHJcg8GASZMmISgoCC4uLmjbti2WLVtW7hj/93//h44dO0Kj0cDPzw/Tp0+v9HxvvPEGfHx8pOaqtLQ0DBs2DC4uLggKCsLGjRvL1UoJgoBPPvkEDz74INzc3PDWW28BAD7++GPcc889UKvVaNu2LdatWye9pqJmsaysLOlaA0BiYiIEQcCuXbsQFhYGV1dXhIeH48yZMyZ5XrhwIXx8fODh4YFJkyahsLCwWtfWWhi8kF0QBAUCA/shMLAfBIHFkuTL6mU5L6/4IYp30rTa4rSioor3NRrvpOl0xWllbz6V7WslmzdvRvPmzfHGG28gLS0NaWlpCA8Px9KlS+Hp6SmlzZkzBwDw5JNP4tdff0VcXByOHz+OMWPGICoqCmfPngUA/Pbbb3jqqacwdepUJCcnIyIiQrrBV+WFF17A7NmzcfToUYSHh2PEiBFS4GM0GtG8eXN89dVXOHXqFF5//XW88sor+Oqrr6TXf/zxx5g2bRqeeeYZpKSk4Pvvv0fr1q3LnUcURcycOROrVq3Cvn37pD4jTzzxBK5du4bExER88803+PTTT5GRkVHu9fPnz8eDDz6IlJQUPPXUU9iyZQtmzpyJ2bNn48SJE/j3v/+NJ598Ert37zbrbwEA8+bNw6JFi3D48GEolUqT9aS++uorzJ8/H2+//TYOHz4MPz8/rFixwuxzWJQoQ9nZ2SIAMTs729ZZISIHoNPpxB9++EH84YcfRJ1OZ+vsiAUFBeKpU6fEgoKC8huLwxZRzMi4k/bWW8VpTz9tuq+ra3F6auqdtCVLitPGjzfdt3Hj4vQTJ+6kffqp2XmfMGGC6OTkJLq5uUmPhx9+WBRFUezXr584c+ZMad/AwEBxyZIlJq9fvXq16OXlZZL2119/iYIgiFevXjVJHzBggDh37lxRFEXx0UcfFaOioky2jxs3rtyxSktNTRUBiAsXLpTSdDqd2Lx5c/Hdd9+t9HVTp04VR48eLf3u7+8vzps3r9L9AYibNm0S//Wvf4nt2rUTL1++LG07ffq0CEA8dOiQlHb27FkRgMm1ASDOmjXL5Ljh4eHi5MmTTdLGjBkjDh061OT9HT16VNqemZkpAhB3794tiqIo7t69WwQg7ty5U9pn69atIgCp/PXu3VucMmWKyXl69uwpdu7cucL3e7fya6n7N7/iEhGRRUVERCA5OVl6fPjhh7U63pEjRyCKIoKDg+Hu7i499uzZg3PnzgEATp8+jd69e5u8ruzvlSm9n1KpRFhYGE6fPi2lffLJJwgLC0OTJk3g7u6OlStX4tKlSwCAjIwMXLt2DQMGDLjrOZ577jkcOHAAv/zyC5o3by6lnzlzBkqlEt26dZPSWrdujYYNG5Y7RlhYmMnvp0+fRp8+fUzS+vTpY5L36urUqZP03M/PDwCk2p/aXFtrYYddIiI5yc0t/unqeifthReAWbMAZZmP9JKmBxeXO2nTpgGTJwNlO89euFB+34kTa5RFNze3CptNaspoNMLJyQlJSUnlOv26u7sDKG6SsSRBEAAUN5k899xzWLRoEXr37g0PDw+8//77+O233wAALqWv110MHDgQX3zxBX766Sc89thjUnpl+a4o3c2t/GzNJfks/bqSNIVCUe5YukqaAkt3/i15vbF0E6KdYc0L2QWtNg/vv98E77/fBFptXoX76PV6xMfHIz4+Hnq9vo5zSFQ91SnLteLmVvwofdNSq4vTNJqK91WU+qhXqYrTnJ2rt68VqdVqGAyGKtO6du0Kg8GAjIwMtG7d2uTh6+sLAOjQoQMOHjxo8rqyv1em9H56vR5JSUlo164dAOCXX35BeHg4pk6diq5du6J169ZSbQ8AeHh4oGXLlti1a9ddzzFixAhs3LgRTz/9NOLi4qT0du3aQa/X4+jRo1LaX3/9haysrCrz3b59e+zbt88kbf/+/Wjfvj2A4hFeQHGH4BI1mdOmffv2Nb621sKaF7Ib+fk3bJ0FIotgWa6eli1bYu/evXjkkUeg0WjQuHFjtGzZErm5udi1axc6d+4MV1dXBAcH47HHHsMTTzyBRYsWoWvXrrhx4wZ+/vlnhIaGYujQofjPf/6D8PBwvPfeexg5ciS2b9+OhISEauXjv//9L9q0aYP27dtjyZIlyMzMlDqstm7dGp9//jl++uknBAUFYd26dTh06BCCgoKk18fExGDKlClo2rQphgwZglu3buHXX3/FjBkzTM7z0EMPYd26dXj88cehVCrx8MMPo127doiMjMQzzzyDjz/+GCqVCrNnz4aLi0u5WpWyXnjhBYwdOxbdunXDgAED8MMPP2Dz5s3YuXMngOJaoV69emHhwoVo2bIlbty4gVdffdWcPxEAYObMmZgwYQLCwsJw3333YcOGDTh58iRatWpl9rEshTUvRERkE2+88QYuXLiAe+65R6olCA8Px5QpUzBu3Dg0adIE7733HgBg9erVeOKJJzB79my0bdsWI0aMwG+//YaAgAAAQK9evfDZZ59h+fLl6NKlC7Zv317tG/XChQvx7rvvonPnzvjll1/w3XffoXHjxgCAKVOmYNSoURg3bhx69uyJmzdvYurUqSavnzBhApYuXYoVK1agY8eOiI6OlkZBlfXwww9j7dq1ePzxx7F582YAwOeffw4fHx/07dsXDz30ECZPngwPDw84l60dK2PkyJFYtmwZ3n//fXTs2BH/+9//sHr1avTv31/a5//+7/+g0+kQFhaGmTNnVnsEVmnjxo3D66+/jpdeegndu3fHxYsX8eyzz5p9HEsSREs3FNaBnJwceHl5ITs7G56enrbODlmAVpuH2Njituu5c3MrXIlXr9dL36SioqKgLNu+T1RDlixb1SnLVSksLERqaiqCgoKqvIGR47ly5QoCAgKwc+fOKjsC26O7lV9L3b/56U9ERGRDP//8M3JzcxEaGoq0tDS8+OKLaNmyJfr27WvrrNktBi9EREQ2pNPp8Morr+D8+fPw8PBAeHg4NmzYUK+n/68KgxciIiIbGjx4MAYPHmzrbMgKgxeyC4KggL9/mPScSK5Ylomsj8EL2QWVygWTJx+ydTaIao1lmcj6+LWAiIiIZIXBCxEREckKm43ILuh0+fjvfzsAAKZNOwWVyrWKVxBZR23ne2FZJrI+Bi9kF0RRRHb2Rek5kVyxLBNZH5uNiIioTgmCgG+//bbOz9uyZUssXbq0zs9LlsfghYiILCYjIwP//ve/0aJFC2g0Gvj6+mLw4ME4cOCAtE9aWhqGDBliw1yaLyYmBl26dLHY8SZOnIiRI0dWuK2goACurq74448/7nqMNWvWoEGDBhbLk5yw2YiIiCxm9OjR0Ol0WLt2LVq1aoW///4bu3btwj///CPt4+vra8Mc2pbBYKhytegdO3YgICAA7dq1q6NcyQ9rXkhWRBEoMrAvAZE9ysrKwr59+/Duu+8iIiICgYGB6NGjB+bOnYthw4ZJ+5VtNtq/fz+6dOkCZ2dnhIWF4dtvv4UgCEhOTgYAJCYmQhAE7Nq1C2FhYXB1dUV4eDjOnDkjHePcuXN48MEH4ePjA3d3d9x7773YuXOnWflPTExEjx494ObmhgYNGqBPnz64ePEi1qxZgwULFuDYsWMQBAGCIGDNmjUAgMWLFyM0NBRubm4ICAjA1KlTkZubKx2zpHYkPj4eHTp0gEajwZNPPom1a9fiu+++k46XmJgovea7777DiBEjAADHjh1DREQEPDw84Onpie7du+Pw4cNITEzEk08+iezsbOkYMTExAIDMzEw88cQTaNiwIVxdXTFkyBCTVa5L8vTtt98iODgYzs7OGDhwIC5fvmzW9bIlBi8kG6IoYslxYPYBAY+s/J0BDFlE6RWl5UCrzav0odcXVntfna6gWvuaw93dHe7u7vj2229RVFRUrdfcunULw4cPR2hoKI4cOYI333wTL730UoX7zps3D4sWLcLhw4ehVCrx1FNPSdtyc3MxdOhQ7Ny5E0ePHsXgwYMxfPhwXLp0qVr50Ov1GDlyJPr164fjx4/jwIEDeOaZZyAIAsaNG4fZs2ejY8eOSEtLQ1paGsaNGwcAUCgU+PDDD3HixAmsXbsWP//8M1588UWTY+fn5yM2NhafffYZTp48iQ8//BBjx45FVFSUdLzw8HAAgNFoRHx8PB588EEAwGOPPYbmzZvj0KFDSEpKwssvvwyVSoXw8HAsXboUnp6e0jHmzJkDoLhJ6vDhw/j+++9x4MABiKKIoUOHQqfTmeTp7bffxtq1a/Hrr78iJycHjzzySLWulT1gsxHZBUEQ0KRJB+l5RQp0Bpy/Vbwt6VIWCnQGuKpZhMm+VKcs10ZsrHul29q0GYrx47dKv3/wQVPodPkV7hsY2A8TJyZKvy9b1hL5+TfK7Td/fvW/JCiVSqxZswaTJ0/GJ598gm7duqFfv3545JFH0KlTpwpfs2HDBgiCgJUrV8LZ2RkdOnTA1atXMXny5HL7vv322+jXrx8A4OWXX8awYcNQWFgIZ2dndO7cGZ07d5b2feutt7BlyxZ8//33mD59epV5z8nJQXZ2NqKjo3HPPfcAANq3by9td3d3h1KpLNfkNWvWLOl5UFAQ3nzzTTz77LNYsWKFlK7T6bBixQqT/Lm4uKCoqKjc8Q4ePAij0SgFM5cuXcILL7wgNSG1adNG2tfLywuCIJgc4+zZs/j+++/x66+/SsfYsGEDAgIC8O2332LMmDFSnj766CP07NkTALB27Vq0b98ev//+O3r06FHl9bI11ryQXVCpXDF16klMnXqS82KQrNX3sjx69Ghcu3YN33//PQYPHozExER069ZNamYp68yZM+jUqROcnZ2ltMpunqUDID8/PwDFHYQBIC8vDy+++CI6dOiABg0awN3dHX/88Ue1a14aNWqEiRMnSjU2y5YtQ1paWpWv2717NwYOHIhmzZrBw8MDTzzxBG7evIm8vDu1Vmq1utLgrazvvvsO0dHRUCiKb8/PP/88nn76aURGRmLhwoU4d+7cXV9/+vRpKJVKKSgBAG9vb7Rt2xanT5+W0pRKJcLCwqTf27VrhwYNGpjsY8/4tZWISEbmzs2tdJtC4WTy+5w5GZXuW3bRyJkzL9QqX6WV9KEYOHAgXn/9dTz99NOYP38+Jk6cWG5fURTL1VBV1iSsUqmk5yWvMRqNAIAXXngBP/30Ez744AO0bt0aLi4uePjhh6HVaqud79WrV+M///kPEhIS8OWXX+LVV1/Fjh070KtXrwr3v3jxIoYOHYopU6bgzTffRKNGjbBv3z5MmjTJpInGxcWl2rVw33//PWJjY6XfY2JiMH78eGzduhXbtm3D/PnzERcXh4ceeqjC11d27Sq6zhXlyRq1hdbAmhciIhlRq90qfSiVztXeV6Vyqda+ltChQweTmojS2rVrh+PHj5v0kTl8+LDZ5/jll18wceJEPPTQQwgNDYWvry8uXLhg9nG6du2KuXPnYv/+/QgJCcHGjRsBFNeeGAwGk30PHz4MvV6PRYsWoVevXggODsa1a9eqdZ6Kjnf27FlcuHABgwYNMkkPDg7Gc889h+3bt2PUqFFYvXp1pcfo0KED9Ho9fvvtNynt5s2b+PPPP02awfR6vcl1PnPmDLKysmQzwonBC9kFnS4fK1Z0xIoVHSttoyeSg/pclm/evIkHHngA69evx/Hjx5GamopNmzbhvffekzqgljV+/HgYjUY888wzOH36tFR7AphXC9C6dWts3rwZycnJOHbsmHTc6kpNTcXcuXNx4MABXLx4Edu3bze54bds2RKpqalITk7GjRs3UFRUhHvuuQd6vR7Lly/H+fPnsW7dOnzyySfVOl/Lli1x/PhxnDlzBjdu3IBOp8N3332HyMhIuLoWNzcWFBRg+vTpSExMxMWLF/Hrr7/i0KFDJnnKzc3Frl27cOPGDeTn56NNmzZ48MEHMXnyZOzbtw/Hjh3Dv/71LzRr1szkb6BSqTBjxgz89ttvOHLkCJ588kn06tVLFv1dAAYvZCdEUcT166dw/fopjiIiWavPZdnd3R09e/bEkiVL0LdvX4SEhOC1117D5MmT8dFHH1X4Gk9PT/zwww9ITk5Gly5dMG/ePLz++usAYNIPpipLlixBw4YNER4ejuHDh2Pw4MHo1q1btV9fMinc6NGjERwcjGeeeQbTp0/Hv//9bwDFfXmioqIQERGBJk2a4IsvvkCXLl2wePFivPvuuwgJCcGGDRtMmnzuZvLkyWjbti3CwsLQpEkT/Prrr/juu+9MAgwnJyfcvHkTTzzxBIKDgzF27FgMGTIECxYsAACEh4djypQpGDduHJo0aYL33nsPQHHzV/fu3REdHY3evXtDFEX8+OOPJs1urq6ueOmllzB+/Hj07t0bLi4uiIuLq/b1sjVBlOF/V05ODry8vJCdnQ1PT09bZ4csQKvNk0ZRzJ2bW2F1dU5+ITq9sUv6/dQbgznaiGrlbsOka7owY3XKclUKCwuRmpqKoKAgs27gjmLDhg3SHCYuLi5Vv8AB3LhxA35+frh8+bLVJ/Fbs2YNZs2ahaysLKsc/27l11L3b37yExGRTX3++edo1aoVmjVrhmPHjuGll17C2LFj603gAgD//PMPFi9eXK9nHzYHgxciIrKp9PR0vP7660hPT4efnx/GjBmDt99+29bZqlPBwcEIDg62dTZkg31eiIjIpl588UVcuHBBam5YsmSJ1GmVLG/ixIlWazKqKwxeiIiISFbYbER2QRAEeHkFSs+J5Iplmcj6GLyQXVCpXDFr1gVbZ4Oo1liWiayPzUZEREQkK2YFLzExMRAEweRReliXKIqIiYmBv78/XFxc0L9/f5w8edLkGEVFRZgxYwYaN24MNzc3jBgxAleuXLHMuyEiIiKHZ3bNS8eOHZGWliY9UlJSpG3vvfceFi9ejI8++giHDh2Cr68vBg4ciFu3bkn7zJo1C1u2bEFcXBz27duH3NxcREdHl1ufgeoXna4AK1fei5Ur74VOV2Dr7BDVmL2VZb1ej/j4eMTHx0Ov19s6O0QWYXbwolQq4evrKz2aNGkCoLjWZenSpZg3bx5GjRqFkJAQrF27Fvn5+dLCVtnZ2Vi1ahUWLVqEyMhIdO3aFevXr0dKSgp27txp2XdGsiKKRly7dhjXrh2GKFZ/PRIie8OybL9atmyJpUuX2jobZAFmBy9nz56Fv78/goKC8Mgjj+D8+fMAihe1Sk9PN1kNU6PRoF+/fti/fz8AICkpCTqdzmQff39/hISESPtUpKioCDk5OSYPIiKyPxMnTpS6FSiVSrRo0QLPPvssMjMzbZ01s61ZswYNGjSwdTaoAmYFLz179sTnn3+On376CStXrkR6ejrCw8Nx8+ZNpKenAwB8fHxMXuPj4yNtS09Ph1qtRsOGDSvdpyKxsbHw8vKSHgEBAeZkm4iI6lBUVBTS0tJw4cIFfPbZZ/jhhx8wdepUW2eLHIhZwcuQIUMwevRohIaGIjIyElu3bgUArF27Vtqn7LwGoihWOddBVfvMnTsX2dnZ0uPy5cvmZJuIyGwJCQnsI1JDGo0Gvr6+aN68OQYNGoRx48Zh+/bt0naDwYBJkyYhKCgILi4uaNu2LZYtWyZtT0lJgUKhwI0bNwAAmZmZUCgUGDNmjLRPbGwsevfuXWkeMjIyMHz4cLi4uCAoKAgbNmwot8/ixYsRGhoKNzc3BAQEYOrUqcjNzQUAJCYmSotDltQkxcTEAADWr1+PsLAweHh4wNfXF+PHj0dGRkatrhmZp1ZDpd3c3BAaGoqzZ89Ko47K1qBkZGRItTG+vr7QarXlqg9L71MRjUYDT09PkwcREZnS6/UVPqq73RrOnz+PhIQEqFQqKc1oNKJ58+b46quvcOrUKbz++ut45ZVX8NVXXwEAQkJC4O3tjT179gAA9u7dC29vb+zdu1c6RmJiIvr161fpeSdOnIgLFy7g559/xtdff40VK1aUCzAUCgU+/PBDnDhxAmvXrsXPP/+MF198EQAQHh6OpUuXwtPTUxqgMmfOHACAVqvFm2++iWPHjuHbb79FamoqJk6caJHrRdVTq0nqioqKcPr0adx///0ICgqCr68vduzYga5duwIo/gPv2bMH7777LgCge/fuUKlU2LFjB8aOHQsASEtLw4kTJ/Dee+/V8q0QEdVvCQkJd91e2cCI6Ohoi+YjPj4e7u7uMBgMKCwsBFBcy1FCpVJhwYIF0u9BQUHYv38/vvrqK4wdOxaCIKBv375ITEzE6NGjkZiYiAkTJmDt2rU4deoUgoODsX//fjz33HMVnv/PP//Etm3bcPDgQfTs2RMAsGrVKrRv395kv1mzZpnk4c0338Szzz6LFStWQK1Ww8vLq9yUIADw1FNPSc9btWqFDz/8ED169EBubi7c3d1rdtHILGYFL3PmzMHw4cPRokULZGRk4K233kJOTg4mTJgAQRAwa9YsvPPOO2jTpg3atGmDd955B66urhg/fjwAwMvLC5MmTcLs2bPh7e2NRo0aYc6cOVIzFNVvrq6NbZ0FIouo72U5IiICH3/8MfLz8/HZZ5/hzz//xIwZM0z2+eSTT/DZZ5/h4sWLKCgogFarRZcuXaTt/fv3x6effgoA2LNnD958802kpqZiz549yM7ORkFBAfr06VPh+U+fPg2lUomwsDAprV27duU63+7evRvvvPMOTp06hZycHOj1ehQWFiIvLw9ubm6Vvr+jR48iJiYGycnJ+Oeff2A0Fo8qu3TpEjp06GDOpaIaMit4uXLlCh599FHcuHEDTZo0Qa9evXDw4EEEBhav4/Hiiy+ioKAAU6dORWZmJnr27Int27fDw8NDOsaSJUugVCoxduxYFBQUYMCAAVizZg2cnJws+85IVtRqN7zwwnVbZ4Oo1mxZlqOiosql6fV6qcYlMjISSqX1V4Vxc3ND69atAQAffvghIiIisGDBArz55psAgK+++grPPfccFi1ahN69e8PDwwPvv/8+fvvtN+kY/fv3x8yZM/HXX3/hxIkTuP/++3Hu3Dns2bMHWVlZ6N69u8m9pTRRFAHcfW2pixcvYujQoZgyZQrefPNNNGrUCPv27cOkSZOg0+kqfV1eXh4GDRqEQYMGYf369WjSpAkuXbqEwYMHQ6vVmn2tqGbMKsVxcXF33V7SoamkU1NFnJ2dsXz5cixfvtycUxMRURWqCkyUSmWdBC9lzZ8/H0OGDMGzzz4Lf39//PLLLwgPDzcZgXTu3DmT15T0e3nrrbfQuXNneHp6ol+/foiNjUVmZuZd+7u0b98eer0ehw8fRo8ePQAAZ86cQVZWlrTP4cOHodfrsWjRIigUxd0/S/rclFCr1eUmUP3jjz9w48YNLFy4UBr5evjwYfMvCtUK1zYiIqoERxxZRv/+/dGxY0e88847AIDWrVvj8OHD+Omnn/Dnn3/itddew6FDh0xeU9LvZf369ejfvz8AoFOnTtBqtdi1a5eUVpG2bdsiKioKkydPxm+//YakpCQ8/fTTcHFxkfa55557oNfrsXz5cpw/fx7r1q3DJ598YnKcli1bIjc3F7t27cKNGzeQn5+PFi1aQK1WS6/7/vvvpRolqjsMXsgu6HQFWLOmP9as6W8XU6oT1RTLcsWef/55rFy5EpcvX8aUKVMwatQojBs3Dj179sTNmzcrnAcmIiICBoNBClQEQcD9998PALjvvvvuer7Vq1cjICAA/fr1w6hRo/DMM8+gadOm0vYuXbpg8eLFePfddxESEoINGzYgNjbW5Bjh4eGYMmUKxo0bhyZNmuC9995DkyZNsGbNGmzatAkdOnTAwoUL8cEHH9Ty6pC5BLGkcVBGcnJy4OXlhezsbA6bdhBabR5iY4t76c+dmwu1unxnuZz8QnR6Y5f0+6k3BsNVXfdV4OQ49Hp9lSN0oqKizGpqqU5ZrkphYSFSU1MRFBQEZ2dns19fWun3aO57IaqJu5VfS92/WfNCREREssIQnIjqLVEUUVTFgvYyrJw2oVQqLT6PC5GtMXghonpJFEWMW/k7jly6+/IlG6/+jq+fDa9ymRMiqjtsNiKieqlAZ8CRS1lV7pd0KQsFuiqqZ4ioTrHmhYjqvdgeItRl5snUGoC5v7O2hcgeMXghu6FSudo6C1RPqZ0AjQUn+bZUWZZ7fxuqn+qi3DJ4IbugVrvhlVfybJ0NolqzRFkuWYE5Pz/fZGI1IjkoWSbBmsv+MHgh2eD8FFRfODk5oUGDBsjIyAAAuLq6ssMwyYLRaMT169fh6upq1c9s3g2IiOyQr68vAEgBDJFcKBQKtGjRwqoBN4MXsgt6fSG++mo0AGDs2G+gVNZuVlEiW7FUWRYEAX5+fmjatOldVzkmsjdqtVpa7NJaGLyQXTAaDTh79kfpOZFcWbosOzk5WbXvAJEccZ4XIiIikhUGL0RERCQrDF6IiIhIVhi8EBERkawweCGiekmv19s6C0RUQwxeiIiISFY4VJrsglrthvnzuY4LyR/LMpH1seaFiIiIZIXBCxEREckKgxeyC3p9ITZtGoNNm8ZAry+0dXaIaoxlmcj6GLyQXTAaDTh16mucOvU1lwcgWWNZJrI+Bi9EREQkKwxeiIiISFYYvBAREZGsMHghIiIiWWHwQkRERLLC4IWIiIhkhcsDkF1QqVwxd26u9JxIrliWiayPwQvZBUEQoFa72TobRLXGskxkfWw2IiIiIllhzQvZBb2+CPHx/wYAREf/D0qlxsY5IqoZlmUi62PNC9kFo1GPY8fW4tixtTAa9bbODlGNsSwTWR+DFyIiIpIVBi9EREQkKwxeiIiISFYYvBAREZGsMHghIiIiWWHwQrKh13PkBhERcZ4XshMqlSvmzMmQnhPJFcsykfUxeCG7IAgC3Nya2DobRLXGskxkfWw2IiIiIllh8EJ2Qa8vwtat07B16zTo9UW2zg7VcxERETV+LcsykfUxeCG7YDTqcfjwChw+vIJTqpOssSwTWR+DFyIiIpKVWgUvsbGxEAQBs2bNktImTpwIQRBMHr169TJ5XVFREWbMmIHGjRvDzc0NI0aMwJUrV2qTFSKiWouKikJ0dDScnJxsnRUiuosaBy+HDh3Cp59+ik6dOpXbFhUVhbS0NOnx448/mmyfNWsWtmzZgri4OOzbtw+5ubmIjo6GwWCoaXaIiIionqhR8JKbm4vHHnsMK1euRMOGDctt12g08PX1lR6NGjWStmVnZ2PVqlVYtGgRIiMj0bVrV6xfvx4pKSnYuXNnzd8JERER1Qs1Cl6mTZuGYcOGITIyssLtiYmJaNq0KYKDgzF58mRkZGRI25KSkqDT6TBo0CApzd/fHyEhIdi/f3+FxysqKkJOTo7Jg4jIWpRKToFFZM/M/g+Ni4vDkSNHcOjQoQq3DxkyBGPGjEFgYCBSU1Px2muv4YEHHkBSUhI0Gg3S09OhVqvL1dj4+PggPT29wmPGxsZiwYIF5maViIiIHJBZwcvly5cxc+ZMbN++Hc7OzhXuM27cOOl5SEgIwsLCEBgYiK1bt2LUqFGVHlsURQiCUOG2uXPn4vnnn5d+z8nJQUBAgDlZJzunUrlg5sxU6TmRXLEsE1mfWcFLUlISMjIy0L17dynNYDBg7969+Oijj1BUVFSul76fnx8CAwNx9uxZAICvry+0Wi0yMzNNal8yMjIQHh5e4Xk1Gg00Go05WSWZEQQFGjRoaetsENUayzKR9ZnV52XAgAFISUlBcnKy9AgLC8Njjz2G5OTkCocX3rx5E5cvX4afnx8AoHv37lCpVNixY4e0T1paGk6cOFFp8EJEZI/0ej3i4+MRHx/PVc+J6pBZNS8eHh4ICQkxSXNzc4O3tzdCQkKQm5uLmJgYjB49Gn5+frhw4QJeeeUVNG7cGA899BAAwMvLC5MmTcLs2bPh7e2NRo0aYc6cOQgNDa20AzA5PoNBi1275gEABgx4G05OahvniKhyer0eCQkJFW5jWSayPot2qXdyckJKSgo+//xzZGVlwc/PDxEREfjyyy/h4eEh7bdkyRIolUqMHTsWBQUFGDBgANasWcOJoeoxg0GHAwc+AAD07x/DD3yqUwMHRsLTteJ+fOZiWSayvloHL4mJidJzFxcX/PTTT1W+xtnZGcuXL8fy5ctre3oiIiKqZ7i2EcnWtm0J7GdARFQPMXghIiIiWeE0kkRERHaudCfxqKioej8LNGteiIiqwOZJIvvC4IWIiIhkpX7XO5HdUKlc8OyzJ6TnRHLFskxkfQxeyC4IggJNm3a0dTaIao1lmcj62GxERFSGUqnEwIGc8ZvIXrHmheyCwaDFL7+8AwC4//5XOCsp1YotR2awLBNZH4MXsgsGgw579iwAAISHv8APfLK60gFNVcHNjh07oXGqXiDEskxkfWw2IiIiIllhzQsRyVrpJqLIyEjs3LnTxjkiqhlORFd9rHkhWUtI4PpGdAcDF6ov9Ho94uPjER8fX+vPQEseq64weCEiIiJZYfBCRLIml2+KRGQ5bFAjIodWl30I9Ho9tm3bZtVzEBGDF7ITSqUznn76d+k5kaUlJCSYdOjt+4B1JqFTKFRo124RevfuzbJMZCUMXsguKBROaNbsXltng6jWBMEJbm7B8Pe/FwqFk62zQ+SQ2OeFiIiIZIXBC9kFg0GLX399H7/++j4MBq2ts0NUY0ajDunpm3HgwCKWZSIrYbMR2QWDQYedO18EANx771ROqW4BpSe8AjjpVV0RRQOuXl2Nq1eBnj2nsyyTXSv7OSEXrHkhItnS6/WcmI6oHmLwQkRERLLCOmQiIrI5rutD5mDNCxEREckKQ1siIrIoa9SisGaGSmPNCxEREckKQ1eyC0qlMyZM2C09L1FYWCiNJunZp6/Vzl8fvtU5+vuzFwqFCsHB76BHjx5cHoBkJyEhQRafEfadO6o3FAontGzZ39bZoDrg6IGiIDjBwyMUgYH9uDxANcl1rhGyHcf61CAiu+DoAUpF6vPNt+zfu67PVx/KF5linxeyCwaDDr///l/8/vt/YTDobJ0dohoTRT0yMrbi8OEVLMtEVsLgheyCwaDFtm3TsW3b9ErXg3FyYhV8dbEa3naMRj0uX/4EP/00k2sb1QG9Xo/4+HjEx8dDr9fbOjvVUp08JyQkyOb92ALr2oioXtqxYycAwdbZsCtcD8u+1GVTnNywVBIRUb1UH/rOOOp7dIx3QUTkIBz1ZlMX7P3aWao5t+yCpBUNb3b0pmP7+ssSUa05+ocW1Zy939zlpCbXUi59WORQTthhl4iIiGTF/sIpInIocpmxszpYq0VkH+T/aUIOQanU4NFH46XnQPl2XSI5UChUaN36dXTr1l0qy1RzdRUwyqGppERVn4v1IchmsxHZBYVCieDgYQgOHgaFonofGloDUGQARFG0cu7IWuQ4R0dVBMEJXl73ok2bodUuy5VxxOvjiBISEir8G9XlfC71bV4Y+w0tiaow9/fiOTo2Xv0dXz8bDkGo33N21IdvW1RMTrUEclTS1FnyHOB1tjeseSG7YDDokJy8BsnJa8yeUj3pUhYKdAYr5YzIPKKox40bO3Hs2FouD0BkJQwjyS4YDFp8992TAIAOHcbAyUl11/1jexQ3FZXUvhDZC6NRj4sXl+HiRSA09JEqyzJgWpMSGRlptbyVblbQ6/WsSSDZYsklWVLXYpkjR6hyL3uzc+SOzY7w9yL7Z8n/ITk14cq1nww/BYiIyGqsHXjWZeBeVzVkVDX2eSHZ2L17t62zQCThzatijjrqRU61KdVRMhJKrrW2DF6IiBxAXQ+rtuUw7spuuGXX+5FTEGWJ4KiyIduOqFbBS2xsLARBwKxZs6Q0URQRExMDf39/uLi4oH///jh58qTJ64qKijBjxgw0btwYbm5uGDFiBK5cuVKbrBAREVE9UePg5dChQ/j000/RqVMnk/T33nsPixcvxkcffYRDhw7B19cXAwcOxK1bt6R9Zs2ahS1btiAuLg779u1Dbm4uoqOjYTBwuCsRkTnkVsNA5jP37yvXpiBz1Ch4yc3NxWOPPYaVK1eiYcOGUrooili6dCnmzZuHUaNGISQkBGvXrkV+fj42btwIAMjOzsaqVauwaNEiREZGomvXrli/fj1SUlLqxQWniimVGjz88Fd4+OGvqjWlekRERB3kyrI4W6rt1cVnjEKhQqtWL2HUqC/qzfIANbmu/H+g2qhR8DJt2jQMGzasXIe11NRUpKenY9CgQVKaRqNBv379sH//fgBAUlISdDqdyT7+/v4ICQmR9imrqKgIOTk5Jg9yLAqFEh07jkHHjmNqPaU6WY6tbzD23Emyshu2IDihYcP70L79w1Ypy/Z8TeyRrctwiYSEhEpbFxISEqr8myqVSmnWX6rBUOm4uDgcOXIEhw4dKrctPT0dAODj42OS7uPjg4sXL0r7qNVqkxqbkn1KXl9WbGwsFixYYG5WiYhqTM7zyZg7N87dAiJbzLPDIcnW4UjX1axSePnyZcycORPbt2+Hs7NzpfuVXWNGFMUq15252z5z587F888/L/2ek5ODgIAAM3JO9s5o1OP06S0AgPbtH4LRCH67pCqJooiiGnaV01qpi50oGpCVdQBffLEPDRr0xsCBg6VaGnuaZK90nqqrdBBjLk40SJZkVglKSkpCRkYGunfvLqUZDAbs3bsXH330Ec6cOQOguHbFz89P2icjI0OqjfH19YVWq0VmZqZJ7UtGRgbCw8MrPK9Go4FGUz/ajusrvb4IX389FgAwd24uFAr+venuRFHEuJW/48gl+1oiwmjU4fz5dwEAXbpssnFuiByTWX1eBgwYgJSUFCQnJ0uPsLAwPPbYY0hOTkarVq3g6+uLHTt2SK/RarXYs2ePFJh0794dKpXKZJ+0tDScOHGi0uCFiKisAp0BRy5l1fo43QK84KKqxXoTdSQhIQGFhYWyqJHkyFGyNrNqXjw8PBASEmKS5ubmBm9vbyl91qxZeOedd9CmTRu0adMG77zzDlxdXTF+/HgAgJeXFyZNmoTZs2fD29sbjRo1wpw5cxAaGmo3bXBcS8W2tm3bhsGDo22dDZKR2B5ijde7Gjqoe5XN2rXBUZS140jXj7OEW47F78ovvvgiCgoKMHXqVGRmZqJnz57Yvn07PDw8pH2WLFkCpVKJsWPHoqCgAAMGDMCaNWvg5GT/336ILMkWgXJCQoLDBeVqJ0BTw48PawYucqDX6y0eIFR1k7b1iClrnb/kf4usr9afXomJiSa/C4KAmJgYxMTEVPoaZ2dnLF++HMuXL6/t6YkcXkUftJaspawogDLnw70uA7DS5+r7QN3V1A4ZEgVXtbLWQ22Liopqda3spRYiISHBbmrK76Y2AUpV7/Fuw54treTvXlfX3F7K2d04zlcvsjg2n5E90+sdo1+FrWshqHL2dhNnX6I7uDAjERERyQq/SpPVVacGx8lJjejoz5CSklKtWUmVSnn3j6rOt205fiOXY56r47v4BKidipelqGpeGVFUolnATCgV4GzRRFbC/yyyC05OKnTs+BjS0oqraVk9Wn3mTNSmVgCl+6faa98WezP399sX7UAigKo6+KoADEQrTxHPeVs3X5Up/f+j1+vt6m8l1+C2bBMSRw7Zlv2UaLJr9fnGZc9EERj/f0k4erl6I2ZaeYp4LtTKmbKA6gRVxTcP640UclE5oXuLBkiq4Vwy53MEaI1ijUdBkeXJNXCi8ngHIrtgNOpx7tw2ZGcnw9Ozm62zU2N1HeRpjcDRy9nV3r/khlqZ6n5LN3eNFDmuqSIIAuIm98D3P/4kpUVERFT5jbtIb8CKpKMAAFHsCsD60cvdbsoGgwHx8fEAzF+NnTd7slcMXsguFC8PMAoAp1SvqbtN1KY1lGr6qIeioqJqNJ+JIAgmNSeuaqcqa1KMBh0GaooXkhWNm1AXwQtRfcPgheo1R+pgWpuJ2krU9XwSdcVR/sZEVIzBix2Re78SueffGuTYXFKX7G0ejbpki8BZLp1M63O5oOrhPC92wpFqAOojvV6P+Ph4u/wb9u3b19ZZcBgVjYIztx8JEdUegxdySCXBRHx8fK2ndKfySt/ESz+31Ddmc/9mpQMIe67hkkvNB5G9Y72+nWLzi3m2bUvAg8Pkfa3qej4OW9X2WaNJQKl0QnR08UrkcglWq3Md7Ln5xJJzMdnz+yT7xJoXIiIikhX5fk0lh+LkpMbAgUvwxx9nLD6lek1rsepzB2RbzHBcuhbIXvuRVKfZR1AocUA7BQDQkcsDkAVx5vE7WPNCdsHJSYVu3aagadNhEAR+4JN8CYISfxii8YchukZlWS5NKLbov2OvQa1cOFLww7tEKXXdB4DDaC3PUn9Da5cFW9+gavLeKuukS1RXWO4sx9afQbXF4MVM9bkpwZqMRgMuXdqLW7dS4O7ewSZ5KP33lEunT7I/omiAr+LU7ecdwBl2SS6qs/yFvWCzEdkFvb4QX3wxGH/++QqMRp1V/oGqmoslISGh0qDlbtsspWQNmuoM77anb6A1yYsjV/+LRh2GaF7BEM0rEI06W2dHVhy5XNgjRVERXNLTofnnH5O0jp9+inu+/tqGOasagxeyO1oDUFTBQ2s/92uqBjaFEtU9ZV4ePC5ehHNGhpSmKCpCpw8/hGbsWAilvmy03bABAyZPRqstW6Q0o1KJlvHx8Pv11zrNt7nY5kF2QRTvrHQ893cBetTfRQSrw5yaKScnyzRb1KY2TO7t63JmT7V0dEffvn2xd+/eau2rys2F8/XrEFq0kNKcCgvRaflyqLOz8XtMDMTbTd5tvvoK92zejHMjR+L0pEkAigOSgF27IBiNUD/yCIoaNgQAFDVoAINaDcFovHMyJyf8OX48Cr29LfROrYPBC1WLuTcfczsjF+iq/wHbvUUDuKjMvyFX9CFe9n3ZSwdqe1xmgIgsKDsb7hcvwuDigoKmTQEAisJCdFqxAuqsLBx6/XUpIGm9aRPu2bwZRdOmAVFRAACDSgX/ffuKA5KcHBQ1agQAKGzYEEWenoCiVMOKkxNOT5iAe7p2hUGjkZLPP/ggzj/0ECCYflk8+8gj1nznFsHgpQqWnvXU3kYY2WMH5AVhIpxVYqXbRwztAUFgzQzVnr2UeVtirZjlOBUUwPnmTRg0GhQ2aQKguMkm5H//gyYzE4fnzZMCEs2iRei/bBnOjxiBU5MnAwCMajX89+yBwmiEJjtbqv0obNQIhQ0aQFCpSp3MCSeeeQZ6V1fonZ2l5NSRI5E6cmS5vJ0fNQqBERHQl65BtVCtrC3U7/9asrmS4KlAWyilqRWA5i7/UzUNXMypPmdVO8mZXEaMyIFCp4MmMxNGpVKq3VBotej46afQZGUh6aWXIN4OKtp8+SVaf/ONaUCiUqH5rl3lAhLRxwdaDw+IpWtIFAqcmjwZehcX6F1cpOTUBx9E6oMPFi+yWqqp6eKwYWa9F0s1IdsDBi9WYukaDXusISGyFgaPllPdfhX1jaKwEAqDAXo3t+LfdToEb9wIzT//4Pj06VJAErxhA1p/8w1Shw/HyWeeAXC7D8nOnVAYDMUBSePGAICihg2hc3UtcyIFTj/1VHFAUqqGxDBjBnaGhpbL14Xba3TR3fEOSPZBcMIh3ZMAgA6C43w7oHqIZdm2DAapOUTQ69EyPh7O//yDPx5//E5AsnEjgr/4AqnR0Tj5738DAIxOTmj17bdQ6PU4869/Sc0+RY0awaBUQig9fYFCgT8mTIDe2dm0yWbECKQ++GC5LJVNi4qK4lxStcTgxc7drX+MI9XGKBQqnNCPvv288v4uZH8q+2avVCoRHR1d4WzFNW3WkEP/jNJlWS+KKGIlUrXlaw13vV75WgOKtAYE/LIXLpmZ+Ct6OIwlTTbfbkH7uC9w8YEBOPbMv2+/wgnt164tDkiGRqOwSXHH2Hw3dwCAMjun1PkUODNqNPQaDQqUamlqhjODh+KPocOLO7WWytvpEQ/d+UVKr16Tdr5WD73+7u+10vdfi/JUretbZnvpkaD2RL53uzpSMnEYUBwgyI21A5y66IAsp1kf6Q576phuK3N/Z8fy6hBEI8KunELiqr1IaBsOnVNxQDI+eRumHPwaCcHheOeBScCB/YDohDPLP4TGoMcUl/tx1as4IJlwRYUFhQX4869/MPtAyXUXkNtpMHQKJf6XrMJ19+J0V49BcJoViVtqV+BAqb9R0OPFP0+Wzp0KFndgl5Q/81633/zXmLw+8e6vP/BLue0DBxqgUlnhGtQSgxcymyUClrLfxkXRgMbC+dvPW4FTqpNcqQQDOrufw9U8ATfFeyDW17IsimhUkAPfWzdxpkkgDIri6zD4z/145NhPONgiFP/r+XDxrhCw/svXoDHocF+zVbji5QMAUBr0aJH9NwKy/75zXEHAtrZ9YBQUEEvdZ79v3xd7grrhb3fT+UnmD3y2XNby1S7l0kheGLxYWE0X9KuqOlwO1eW1IRp1GO78/O3nm2DPwYsjNdfZA3Mm67JHZUdwiKIO3Qyz0c0Z6BC6CQon50peKV+KoqLi5prbI2W8T51EswP7kR3YEhcjBxbvJIoYOXYinHQ6bPvfSuT7+AIAWmX+g67nkxDcVI3g3neaJLJ/bIcGHh54LVSH3GbF6S5teiHxgSCITX2wqLEolRVD7zkAgNnFJ7p9BM/bj9JptldR+R44sPhL344ddfu5XlUtdkV5rcmcWnWBn7pEDs6RhkfKjdpJvlNpuF2+DO8TJ1Do7Y2MHj2KE0URkRMnwvmff7Br5UoU+BYHJN6XLqDND98jvVcvpA++HbxAQGHjxlAWFMCtIA+G29chq0snHJsxA7datDCZEuH32HcQEREB3e7dKJlGzejTGLk+xSN5NABc1U53nUbBHlWUZ1f17ble6vi9VHX9Ktpur3NqcW0jC6pprYu5KlokMCEhoVoL+jkKpVIpfXspUReLJxJZk7UXJlTm58Prr7/gceGCSXqP+fMRMXkyXP6+0zzT+PhxdFqxAi127LizoyDAoFYDAJxv3pSSs4KDcW7UKFzt29fkuLs//hg71q1Dzj33SGl5AQG4PGgQstq1s+A7o/qGNS8WUFdBi7WV7sdiyWaqumzyqu65LNUJ2FH+9lRzdxtVVWeMRjhptTCUGrYbvHEjXNPS8MfEidLEaM1270boJ58gvVcvHJ43T9rXLS0NbunpcLlxAwU+xf1NbgUGIr1nT2SWCTIOvvkmdO7u0vwoAJDdpg2y27Qpny+5VjtZiJyGRMsprwCDF7KSyiYZk8sIFFuMbrLWOZVKJYYMicLsAz9Z5fj2rGSEoKMEmE0PH4Zrejqu3n8/dF5eAIBmP/+MTsuX43r37jj86qvSvs0SE+GWloZLgwdLwUtBkyYobNSo3ERqx6dNg1GpRE5QkJT2T0gI/gkJKZeHkqYiIlti8EJ2z8nJCVFRUbW+Acm5UyhZVmRkJJyd7aMjrSo7Gy7XryOvWTMpzfvYMQR/8QVyAwKQMm2alN7x00/hlpaGnJYt8c/t4EXn7g4nvd6kGQcAUocPh0KnQ8HtydYAIKNHD+ws6b9Sys3OnS3yXuTe+bquOEowbUsMXmohISGhWnO/lOxXdlSKo48guhs5vvfaNgvI8T2T5bhdvgz/ffug9fQ0WZPmvhdegFtaGva/846U5lRUBO+TJ+FUWGhyjOtduiCnZUuT5qGboaHYtWoVChs2NNn3wvDhVnon9qFs02/p5u6SJhD+zzkuBi9kHwQnHNU9CoBTqpN8OBUUAElJgCgCYWEAAEFwQrdzzeFy/Tqc1X8gK7QLgOJ+JW03bkTWPfeYBC/5TZvCqaAAyqIiKS27TRsceeEF5JdpojkxdWq5PBhcXFDgwnlLqH7haKNaklOHzYpG49hLBy2FQoVk/WNI1j8GhcL+ZnOkeqhMv63ArVvR6cMP4ZGaKqV5p6RA2asXUKppR6FQofexRohMKILbjUwpPTcgAJcGDcLV/v1NjvvbggXYuW4dMm4HP0DxAn/X+vZFVnBwldm09ggle1IfZ9qu6d+3otc5Uk0Ua15kSq/XW21ytJoWcEf6x6jvqrpJWLJfg6X7SJSM/ilRUYAu/e+kpAB//gn06AEEBBSn/fILMGYM0Lw5cPCg9Bq//fvR+Phx3AwNxa3bHVsLmjaF6OcH4fYInRKnnnwSEATk+ftLafl+fjg+Y0b5DJcakVOT1bRtfUMv3XxT0ZxCXN6DrIHBSz1UurbIXr61iaIRDYQrt583BysFqVaMxuKak9trsrheu4Z7Nm+GUaVCwJYtd4KXWbOAn38G1q0D/vWv4jQvL+Dvv8vVvFyJiMCNTp2Q07KllNZnyhQI06cX/3I7SBJFI/5uVnxeZ2dNbVaiIQdQlzNw16dAkcEL2QXRqMVDztNuP98EwD5GgtSWvTTLOSKnwkLce/kc3LQFAAbc2fD448BXXwH/93/AY48V71tUhMCffkKRp6fpQbp2BQoKgFJzlqBtW+DIEaBFC5Ndr1RzeL/RqMWpU8UBTZcum+DkgMsD2Bs5LppLtcPghayiJtXfleHaQZaltdyfplrytQYUmXnO/EKdyWta7P4ZTU6k4FLf/rh+e1iv69U0bNr4Mm66eAJ4+c7OKhWg1QKl+qbk+/riz0ceQb6PD0LEUuvefPBB+ZNrNMVBDSDVppC8cKSR4+Ndgaiemft7HTdkHEgEKmg88bl1A63+uYqrnk1xqaEfAKBFZho2xs2DYNRj9rTPpX0/2HUc957Yha90/liR3wUA4FHkgw4NfHHFqylCSwcZ8+cXP0rNm2JwccGft2thQux0rRZ7YqnmB9aIkLWwY0EVuKgdOQIXlRO6t2hgk3N7FdzC83vX4a2f/muS/sLedfgibh6i//hFSsty8UDznAz45P4DF+2dOU4SgsPx/v2P40CLTlLaLY0b+v37M7z2zNtwcdHcOXBgYPGDNXZUTXXV9y8qKooBnYXwv7ueq+23K7lWzVqyWUsOBEFA3OQe+P5HCy0RIIpQ5+TA9fp1ZLVuLSUHb/4GQdt/wvmoKJwdOQoA8ECX+9CwZfEcPu5zn4JBU9wHpO3VZsjJao7BbTS4p3dJU44bdvu/j/wmTfBWQzUg3E7v3QNAD7QC8C+UavYBoFbY78q3RLVVXzrgmovBC5lNrgGLNcjpWqhUKjw4zPxlFjz/+gvep04hJzBQmkZedSsXgycUj875cdMmGG/P+KrRFcE9PQ1eadeguV1p6eLTGOdHjEBRw4bQwAjD7fQL48bgwrgxxa8rdb68Du3KpZE82LKmuvQQeXaUd3xsNqpDlpzQbufOnTX+B+U/ds2YG6gkJCTYdXAj6PVQ5udLvyt0OnR9/330mT0bilKzvfoePIiOK1fC/5c7zTs6d3do3dxQ4O0NTXa2lH4lIgL7Y2Nx5nb/khKnJk/GuYcfhqHMgoA11bdvX4scx5KUSmWlTQKWXojU3hc2tRW9Xo/4+Pg6nTi0qr9FSVAVHR1d5eCDqKgoREdHV6tpqb6XAda8kH0QnJCiK25m4PIAlhEVFQWlk1PxRGupqcDo0dKEaK02b0b7tWtxadAgaeE/o1KJpocOQVVQANeMDOTenrStxcMPI+3yZWS1aXPn4IKAHevXQyzzYZzv54d8P7+6eYO3RUZG2tWINEFwgo/PQ9Jzc1S35qLkhmjpGb5L1mArLLOmkqNzcnKy+PWsaD27ulSTBW3t6f+oKmbVvHz88cfo1KkTPD094enpid69e2Pbtm3S9okTJ0IQBJNHr169TI5RVFSEGTNmoHHjxnBzc8OIESNw5coVy7wbKyjdN4I1FtajUKhwWP8UDuuf4vIANdTozB/o8NlnaFH2w2rQoOL5Ti5ckJK0np4QjEY437hxZz9BwMlnnsHhuXNR2KiRlGwYMgRJc+fi8uDBJoctG7hQMYVChebNn0Lz5izL9UFJIMmOuHXLrE+f5s2bY+HChWh9u4Pe2rVr8eCDD+Lo0aPo2LEjgOJoc/Xq1dJr1Gq1yTFmzZqFH374AXFxcfD29sbs2bMRHR2NpKQkjuypA3JZh4nuUObnQ1lQgEJvbykt7K234Jmail9i3gRQPCTY89IltPruO2R064ZLJR+kggDcd1/xRGxarfT69N69saNbNxQ1aGByrupOxFZfWKuGg+zL3WocSq9WTfbDrOBleJkl1t9++218/PHHOHjwoBS8aDQa+JZZCbVEdnY2Vq1ahXXr1kntdevXr0dAQAB27tyJwWW+2VH9IYpGuAs3bj9vjPrYHavBH3/A/coVZNx7L7ReXgCA5rt2ocvSpcjo1g2/L1gg7ev6999wzciAW9o1lAQv/7QJxrmRI5FdavQPAKBU7WjJpGt6NzfoS88qKzOlv+jYW1W3KBqh1V4HAKjVTSAI9a8s2yu5BCLVCZrre2Bd4/96g8GATZs2IS8vD71795bSExMT0bRpUzRo0AD9+vXD22+/jaZNmwIAkpKSoNPpMGjQIGl/f39/hISEYP/+/QxeHIi5i+2JRi3GOE+6/dxxlgeoiPvFi2ixfTt0Hh44+8gjUnqn5cvheekSDi5YgBvdugEACho3BgCo8vJMjnFy8mQYlUrcbNESSC5Oy2nZEqcnTaqLt2CW+jYs3WjU4sSJpwFweQBzOVJZKR1U2Lr/S03Ye57NzllKSgp69+6NwsJCuLu7Y8uWLejQoQMAYMiQIRgzZgwCAwORmpqK1157DQ888ACSkpKg0WiQnp4OtVqNhg0bmhzTx8cH6enplZ6zqKgIRaVGP+Tk5Jib7Xql9D9NfemRbu6U9zWZsr4qyrw8aG7dQl6pmsduH32IpseOIWnGf3C9U/EwY6+bmWj1/ffIad4cJ8bcCV5udOiIgoaNUOSkkvKW1rYDvvviS+hdXIFS+b3WsXiyNktN9W/Nb6Scp0J+7PmmVYLdDOo3s0to27ZtkZycjKysLHzzzTeYMGEC9uzZgw4dOmDcuHHSfiEhIQgLC0NgYCC2bt2KUaNGVXpMURTvOslUbGwsFpSqMicqy+wp7w/sR0VT1ldH8PULaHv9Io76t8WVBsWBSp8Lydjw5av407sFBj29Qtr3s9QsBF3PQMKvadiY1wUA0CQ3EM/cOxJ/eQfgywOl8tB5GtAZwC0AB0oS1bcfVFMlQ5gdoWrd3r8NOzq5NDvVB2Y3xqrVarRu3RphYWGIjY1F586dsWzZsgr39fPzQ2BgIM6ePQsA8PX1hVarRWZmpsl+GRkZ8PHxqfScc+fORXZ2tvS4fPmyudmuMUt9a0xISKj2aCVL1pY4UjVsWS4qJ7TyFKvesYaa5P6DGb9+gRf2rDVJn5u4Gst/eB/3XUiW0i57FZdfV10hUGrhvw/7PIpRj72P+Hb3S2nX3Rvh7QeexpedLddM2spThJpdK+qUNQIJR6ylMmc0Tl3WpjhKrXTpeWQ0Go1Jem2OZe9Bcq1zJ4qiSZNOaTdv3sTly5fhd3veh+7du0OlUmHHjh0YO3YsACAtLQ0nTpzAe++9V+k5NBqNyR/FnlgyOIiMjISzszOHZFeTIAh4LhTQGs0PYO7r3Bn7jh2Tfm/z7WY0O3AAf0UPx5X7iydAc79WgMH/3QC9Wg3flx4vHrkDwOdCMG5oCjAy1BX3lkxrb2yK7/p/Cb2rKxYBgDSFfam5UWCdQCsiIgK/7tkNzpBfdyoLXJRKJYYMGYLk5LrPE1XNnmfhtWWtmhxr9MzK7SuvvIIhQ4YgICAAt27dQlxcHBITE5GQkIDc3FzExMRg9OjR8PPzw4ULF/DKK6+gcePGeOih4gmbvLy8MGnSJMyePRve3t5o1KgR5syZg9DQUFlEwQaDQRZVhqX/KR255gUojic0lXxZU+blwe3qVeg8PZF/ux+KS0YG7p81CypRhGbjRikg8byeAe8zfyAzNATX+xcHL3q/prgUGYl8Pz+4iHoYlcVzdpz/13icx3gApaawd1IAKlfYohXeVe1kN4FLybdrR2iiodqR4w3RHPb4/koHZ4D9BWiWZNaV//vvv/H4448jLS0NXl5e6NSpExISEjBw4EAUFBQgJSUFn3/+ObKysuDn54eIiAh8+eWX8PDwkI6xZMkSKJVKjB07FgUFBRgwYADWrFnDzlcoDjTi4+MB1LxtVa/XyyLAsiSFVgu//fvh8vff+GvsWCkgabd2LVpu24azY8fizOOPAwAKGzaEKi8PgtEIdVYWtLc7j19+4AHcDAkxGWYsKpU4PnOm1fNv63Z0e5uhlkxFRETYVXOSNctrde4DJUGDI92Y7TEQsndmXa1Vq1ZVus3FxQU//VT1irXOzs5Yvnw5li9fbs6pydEJTjitHwbg7ssDNDl8GIqvvgL69weeeEJK77J4MQRRxKWoKGmOlDx/fxQ2agRRcacziKhSYc+HH6KgaVMYXFyk9OzgYGQHB1v4TclfVFSUQ90kaqrsN9q7USiUCAubClE0wmi8+824dCBQWVBQ37/Y2Tq4tzcMdIrxCpBdcBIVOJn7BG5p3DBGIUIwGNAjJgZuV68CKSnA7dllPS5dgmLduuLJ1m4HL0a1Gmnh4dC7ukIodaNNffBBpI4cWe5cuYGBdfKeqH5SKjUYNuy/9XbyMKK6wOCF6pQmMxNuV68iq00bGG93wm6xbRtCPv0UaHs/no+eDQAQnZzgfvkyXG7ehP7sWSl4udmpEwwLFsApPNzkuEdefrn8yeylIwhRBeS6Fo5SqbSL2hBzasPsFWtRao5XrRS5Tbdc0YeHrT9QSjjfvIlGKSkwaDT4u9QMzPf/5z9wzsrCL4sXI/v2KsVaT08Iej2aFl6FBtkQRQ8AAo7PmAG9qyt6hoZKr89u3RpiVBSgVEpT3cuNNTu12vuHoa36b9TljU4UReTn34Ber69yDqv6rKL/g4rm5KnLeXrkdg+oz+z3U44swtxp+muixbZtaPDXXzg3ahTymhWvs9Pw1Cl0W7QI/7RvbxK85LZoAYOzM5T5+VLa9a5dsfWTFUhKn4rxeExaHuB69+7FO8h4DR57U/LhXFhYaDeBLuBYo5R0unx88EHxkihcHkD+ajq8uvRIT71eb9dfKuSIV7OWHHoosigCRiNwu8Og119/oe369dC5u+PonDnSbs1370aj06dxvUsXKXi51aIFbnbsiKw2bUwOefDNNwGF6WxqBldXFCgVwO0VIvr264tDB3634hsr7243T3u6sdpDdf3dJCQk1Mm0B/Y0+obusEVTTn1cDoUYvBAAQaeDJjMThbcX0ASATsuWwe/AARz7z3+QXtK/RBTRNCkJhWXWproSEYEbXbogNyBASssNDMSBhQvLn0xRs2lgHaF92xJsUa1d+htjbYOn+j5yxpZKDzG2h0Dc0djjMhSO3AzG4KUO2fobsyonB+6XL6OwSRMU3A5UPP/6C/fNno2iRo2wa/VqaV+FXg/V7UneSuQGBOD4tGnIbdasuFbmdlv+pSFD6vaNWEFCQgKioqLs7sOnrtjjB6+jKt0vqb4OQ7fVXC0l/+dKpZJfiGSOwUspoiiiQFfcDPTAoOL5LXbsuBNwVLQSsSVXJ07YeacqvDbHdSoqRNNjx+H8z02kRt0JLDp/8j8E/LIXKROexJ8PFS+UmdXYBwqjEcr8fOjyi6QRQKceHovTDz2MXD8/GEvyoXLG2YG3R0gYa5a3ymgtfLyaqMsJ/ux5mvKy5JRXIqofGLyUUqAzoMPrZSfaKzVS4MAvKLcScS1WJ76rA4nVOm7PSynom3oESc3a4+fWPQAA3nlFSProTRgh4FHXAShSFQck/xEDMNazKb69pMAqaTVjT7w1dS0y3BsBR0qfr0Xxj2sWe0d3pYSAx12q3s8RlK7CletwWSIiW+I6tDKh0RXhnYSP8MUXc6Ey6KT0+y4kY9rBTRhw7k4H15uuXviteUfEt78fbrpCKf3D8Edw37P/h1U9HjI5doaHt13NieLurLabm7q9Dz0mMldJE2Ftj1GC/ZjIFvipXIqLygmn3hgMAOWajACgT58++PXXX6tMM5eg00F0cpI6szbfuwfdd2zHudZtcGLCxOKdRBVGrNgDVWEB/huYhlu3O8c2dQ3BOa8ctOrcBYt631m1+Ep4cWfZecUvLnU266xsXBt9+/ZFYuIuXLv8ABQC4OSkqtbr6ks/jYrmRnHEphxH+XsqFEp07jwBomiEcJelLqytuqPk7H0EG1FFGLyUIggCXNW3O9Ipyq9W7Kp2qlZaZRSFhVDl5aHo9myxEEXc99xz8LxwAYkrViDf3x8AoDHooEw5jkZOilLHFvDHxAnQu7pC9G4gpWd364Lsbl2KX2fe27UbrmonuKhUuKfVcwCKp1dnZzp5MPfG58ijH0oolRqMHLnGYu/R2v8L/F8jOWKzkRmqWz2qzM1Fgz/+MJmIrdnPP2PomDHoXHpBSkGAYDRCYTDA/coVKflGly7QbtqE49OmmRz34rBhuBoRAV2pVbodAefsqBt6vR7x8fGIj4+vs/mJKpt3wxJNF0T1SVRUFKKjo9mMfRuvQi045efDae9eND18GBlhYVJ6nxdfhMfly/htwQJc79YNAKShyZrMTJNjJD/3HHTu7ihs3FhKK2zcGMbISBTUo6pcURRhNBZJz+2FNWoIzK2iryporsk359LH5IehZYmiCJ2u+IvLsGHDIAiCQzTtyXWItyWbxVhLZT/4qVVNjU6cwMUtW+AeFiZNxtbg7Fm4vfoqOjRrZhK83GrRAsqCAjgVFEhpWW3bYvv69dB6eZkc91ZQULXz0Ldv3zqZ7t8WjMYiJCePAQAMHJgFlcqrilfcna3W0CHS6fIRG+sOAJg7Nxdqdd0tb2GL/ivWvKHbaj6Ymip9LQoLC6vYm2qDwUtZf/8NfPYZFJmZQP/+UvI9mzfD59AhGNRqKXjJDQhAnp8fbrVoYTJp29EXXijugFuKUaUqF7iYi736CSg/1Lo+1JyUvimzgykROf6nnrny84FXX4WgUkG4/34pCLnRpQsMarXU/AMARY0aYfenn5Y7RNnAhYiIiCyHwUtZgYHAxIkwtm59ZwgzgNQRI5A6YoSNM0fV5Qi1VHKrMifbYgfoijnKEHwyxeClLIUCWL0aol4PIws7EdVz9aVpsr6Re/Mrh0rbAbks4x4REWHrLFRLRR+2dXGN5fJ3tLaSTosc1um4+DeuOV47y2DwYmVV3fAjIyPrZQGujzf6+vieiYisof7dNWXMEfpxVEYQFGjQoA8AQKFw3PdZXayqLybHfj8KhRM6dHhYen43jrjMQ2VqM6Sa86vUTtl+P47QP4qfjjLiyB3PFAo17rnnZQCAUuls49zUTxV9wJUNoCxZBs25IcnpJq9UOmPMmE21OkZ1+yOUvYZlr4293/Qd+TPN3tl72agKm42IiIhIVljzQhZT0bdFufdop/LkVAtiL+T+Lbe+4d/L/jF4IbtgMBRKywM88EAWlMrazUZsbxyhjZmqR6vNs9nyAJZiy5s3AweqDgYvZDFVtV+zFqZq7KhLZBlyXUjSWhwtKGSfFyuwxJBYW93krTmcl/MbUFVKAmAiorvhHYTqrZKbpK1HOjjaN6K6xJoqcmT8bKgc/+utwGAwVPhcruQ414a1RUREYPfu3bbOBjko3rSI7o7NRjKjVCot1rRTm+M4+myxUVFRtWresvfrwyY8ckQs1/UHgxcisqiEhATW0hGRVTE0JbsgCAp4eoYBkMfyAI68VAPVjkLhhDZthkrPicjyGLwQANv3a1Eo1GjTZj4ALg9Q38m9v4dS6Yzx47faOhtEDo3BC1VbTWZWlfuNiIiI7A+DF7Kauw1j5YJs9Q+HNVNFHPULTumyznJvebyiZBcMhkIcP/4vAMADD/xt0eUBSm6a0dHR0Ov1DJjqkKPemO5Gq83DBx80BQDMmZNhl8sD1Me/CzkWBi/VUJtp7e2lY6c9Tc1fWQBhNBbZIDfWYcmapbI3Go7ksX86Xb7Zr2FAQVR9HCptITWd16Oq6kR7ny+EiIiorjF4ISIiIllhs1E9Zo1qalZ9ExGRtTF4sYLS/VzM7fNSk74pFS0waE99XMwl15757JtCtsYvD1RfyPMuQVVyhAUhrYEBBhGR/DF4IRO2GlIsCAJatOgLQRAgCOyKRfIlCAoEBvaTnhOR5TF4kbnIyEipmUXOTUUKhQaPP75Ltk1GjsLcZgfOn1OeSuWCiRMTbZ0NIofGO0U9UNInhqyDM8cSUVnsf2Rd/MStIVvWctjbtNO2vHnb8zf/+tS/pi6Xe+BNgUpjeaifbH/nI0Lx8gBLlvgBAGbOvGCXU6oTVYdWm4dly1oCYFkmshazepN9/PHH6NSpEzw9PeHp6YnevXtj27Zt0nZRFBETEwN/f3+4uLigf//+OHnypMkxioqKMGPGDDRu3Bhubm4YMWIErly5Ypl3Y4dKvpFS1fLzbyA//4ZFjlXybSw6OvqutULV3Y/IHJYsy1R7/D93PGb9FZs3b46FCxeidevWAIC1a9fiwQcfxNGjR9GxY0e89957WLx4MdasWYPg4GC89dZbGDhwIM6cOQMPDw8AwKxZs/DDDz8gLi4O3t7emD17NqKjo5GUlGQ36wDVtdLVno7ctGDvzKl+ZlU1EZHtmBW8DB8+3OT3t99+Gx9//DEOHjyIDh06YOnSpZg3bx5GjRoFoDi48fHxwcaNG/Hvf/8b2dnZWLVqFdatWyet2bN+/XoEBARg586dGDx4sIXeVu3dLaBgTQrVFjv5EhHVXI0nITAYDIiLi0NeXh569+6N1NRUpKenY9CgQdI+Go0G/fr1w/79+wEASUlJ0Ol0Jvv4+/sjJCRE2qciRUVFyMnJMXmQeUqCMQZeREQkd2Z/9UtJSUHv3r1RWFgId3d3bNmyBR06dJCCDx8fH5P9fXx8cPHiRQBAeno61Go1GjZsWG6f9PT0Ss8ZGxuLBQsWmJtVu1TRVP7W4OTkVGGzBps7iIhI7swOXtq2bYvk5GRkZWXhm2++wYQJE7Bnzx5puyAIJvuLolgurayq9pk7dy6ef/556fecnBwEBASYm3Wbqc2QWUcPNkqaT7TaPCQn1/w4jn6dLIHXiIgchdnBi1qtljrshoWF4dChQ1i2bBleeuklAMW1K35+ftL+GRkZUm2Mr68vtFotMjMzTWpfMjIyEB4eXuk5NRoNNBqNuVmtkdLzhbBfQt0RBAX8/cOk50RyxbJMZH21vjOLooiioiIEBQXB19cXO3bsQNeuXQEAWq0We/bswbvvvgsA6N69O1QqFXbs2IGxY8cCANLS0nDixAm89957tc0KyZhK5YLJkw/ZOhvllO24XZvmPtZ81A/2WpaJHIlZwcsrr7yCIUOGICAgALdu3UJcXBwSExORkJAAQRAwa9YsvPPOO2jTpg3atGmDd955B66urhg/fjwAwMvLC5MmTcLs2bPh7e2NRo0aYc6cOQgNDZVGH9m7ktoYDmkmIiKyDbOCl7///huPP/440tLS4OXlhU6dOiEhIQEDBw4EALz44osoKCjA1KlTkZmZiZ49e2L79u3SHC8AsGTJEiiVSowdOxYFBQUYMGAA1qxZI7s5XmozHbojfAN3hPcgV7z2RFTfmRW8rFq16q7bBUFATEwMYmJiKt3H2dkZy5cvx/Lly805td2zh4nm7G3NI3PodPn47387AACmTTsFlcrVxjkiqhmWZSLrk9cdzkbq4zfdun7PoigiO/ui9JxIrliWiayPwUst1MeghoiIyNYYvBBVE4NVIiL7wEkIiIiISFYYvBAREZGsMHghIiIiWWGfF7ILgiCgSZMO0nMiuWJZJrI+Bi9kF1QqV0ydetLW2SCqNZZlIutjsxERERHJCoMXO1Oy7AARERFVjM1GVsD5QMyn0+Vj5cp7AQCTJx/ilOokWyzLRNbH4IXsgiiKuH79lPScSK5Ylomsj8GLDLFmh4iI6jMGL3aIwQnJEcstEdUVdtglIiIiWWHwQkRERLLC4IWIiIhkhX1eyC4IggAvr0DpOZFcsSwTWR+DF6pzFXXsVKlcMWvWBdtkiMiCWJaJrI/NRkRERCQrDF6IiIhIVthsRHZBpyvAmjV9AQATJ+6FSuVi4xwR1QzLMpH1MXghuyCKRly7dlh6TiRXLMtE1sfgxYFwhlMiIqoP2OeFiIiIZIXBCxEREckKgxciIiKSFQYvREREJCvssEt2w9W1sa2zQGQRLMtE1sXgheyCWu2GF164butsENUayzKR9bHZiIiIiGSFwQsRERHJCpuNyC7odAXYsGEIAOCxx7ZxSnWSLZZlIutj8EJ2QRSNuHhxj/ScSK5Ylomsj8FLGZxin4iIyL6xzwsRERHJCoMXIiIikhUGL0RERCQrDF6IiIhIVthhl+yGSuVq6ywQWQTLMpF1CaIoirbOhLlycnLg5eWF7OxseHp62jo7REREVA2Wun+z2YiIiIhkhcELERERyQqDF7ILen0hNm4cho0bh0GvL7R1dohqjGWZyPrYYZfsgtFowNmzP0rPieSKZZnI+ljzQkRERLLC4IWIiIhkhcELERERyQqDFyIiIpIVBi9EREQkK7IcbVQyKXBOTo6Nc0KWotXmofD2qNKcnByo1RylQfLEskxUuZL7dm0n95fl8gBXrlxBQECArbNBRERENXD58mU0b968xq+XZfBiNBpx7do1eHh4QBAEW2enxnJychAQEIDLly9zjaZK8BpVD69T1XiNqofXqWq8RlWr7BqJoohbt27B398fCkXNe67IstlIoVDUKmKzN56envwHqAKvUfXwOlWN16h6eJ2qxmtUtYqukZeXV62Pyw67REREJCsMXoiIiEhWGLzYkEajwfz586HRaGydFbvFa1Q9vE5V4zWqHl6nqvEaVc3a10iWHXaJiIio/mLNCxEREckKgxciIiKSFQYvREREJCsMXoiIiEhWGLzUgbfffhvh4eFwdXVFgwYNKtzn0qVLGD58ONzc3NC4cWP85z//gVarNdknJSUF/fr1g4uLC5o1a4Y33nij1utD2LsVK1YgKCgIzs7O6N69O3755RdbZ6nO7N27F8OHD4e/vz8EQcC3335rsl0URcTExMDf3x8uLi7o378/Tp48abJPUVERZsyYgcaNG8PNzQ0jRozAlStX6vBdWFdsbCzuvfdeeHh4oGnTphg5ciTOnDljsk99v04ff/wxOnXqJE0W1rt3b2zbtk3aXt+vT0ViY2MhCAJmzZolpfE6ATExMRAEweTh6+srba/TaySS1b3++uvi4sWLxeeff1708vIqt12v14shISFiRESEeOTIEXHHjh2iv7+/OH36dGmf7Oxs0cfHR3zkkUfElJQU8ZtvvhE9PDzEDz74oA7fSd2Ki4sTVSqVuHLlSvHUqVPizJkzRTc3N/HixYu2zlqd+PHHH8V58+aJ33zzjQhA3LJli8n2hQsXih4eHuI333wjpqSkiOPGjRP9/PzEnJwcaZ8pU6aIzZo1E3fs2CEeOXJEjIiIEDt37izq9fo6fjfWMXjwYHH16tXiiRMnxOTkZHHYsGFiixYtxNzcXGmf+n6dvv/+e3Hr1q3imTNnxDNnzoivvPKKqFKpxBMnToiiyOtT1u+//y62bNlS7NSpkzhz5kwpnddJFOfPny927NhRTEtLkx4ZGRnS9rq8Rgxe6tDq1asrDF5+/PFHUaFQiFevXpXSvvjiC1Gj0YjZ2dmiKIriihUrRC8vL7GwsFDaJzY2VvT39xeNRqPV824LPXr0EKdMmWKS1q5dO/Hll1+2UY5sp2zwYjQaRV9fX3HhwoVSWmFhoejl5SV+8sknoiiKYlZWlqhSqcS4uDhpn6tXr4oKhUJMSEios7zXpYyMDBGAuGfPHlEUeZ0q07BhQ/Gzzz7j9Snj1q1bYps2bcQdO3aI/fr1k4IXXqdi8+fPFzt37lzhtrq+Rmw2sgMHDhxASEgI/P39pbTBgwejqKgISUlJ0j79+vUzmfBn8ODBuHbtGi5cuFDXWbY6rVaLpKQkDBo0yCR90KBB2L9/v41yZT9SU1ORnp5ucn00Gg369esnXZ+kpCTodDqTffz9/RESEuKw1zA7OxsA0KhRIwC8TmUZDAbExcUhLy8PvXv35vUpY9q0aRg2bBgiIyNN0nmd7jh79iz8/f0RFBSERx55BOfPnwdQ99eIwYsdSE9Ph4+Pj0law4YNoVarkZ6eXuk+Jb+X7ONIbty4AYPBUOF7dsT3a66Sa3C365Oeng61Wo2GDRtWuo8jEUURzz//PO677z6EhIQA4HUqkZKSAnd3d2g0GkyZMgVbtmxBhw4deH1KiYuLw5EjRxAbG1tuG69TsZ49e+Lzzz/HTz/9hJUrVyI9PR3h4eG4efNmnV8jBi81VFHHpbKPw4cPV/t4giCUSxNF0SS97D7i7c66Fb3WUVT0nh35/ZqrJtfHUa/h9OnTcfz4cXzxxRflttX369S2bVskJyfj4MGDePbZZzFhwgScOnVK2l7fr8/ly5cxc+ZMrF+/Hs7OzpXuV9+v05AhQzB69GiEhoYiMjISW7duBQCsXbtW2qeurhGDlxqaPn06Tp8+fddHybe/qvj6+paLOjMzM6HT6aQotqJ9MjIyAJSPdB1B48aN4eTkVOF7dsT3a66SHv53uz6+vr7QarXIzMysdB9HMWPGDHz//ffYvXs3mjdvLqXzOhVTq9Vo3bo1wsLCEBsbi86dO2PZsmW8PrclJSUhIyMD3bt3h1KphFKpxJ49e/Dhhx9CqVRWWstd365TWW5ubggNDcXZs2frvCwxeKmhxo0bo127dnd93C2CL6137944ceIE0tLSpLTt27dDo9Gge/fu0j579+41GT69fft2+Pv7o2XLlhZ9b/ZArVaje/fu2LFjh0n6jh07EB4ebqNc2Y+goCD4+vqaXB+tVos9e/ZI16d79+5QqVQm+6SlpeHEiRMOcw1FUcT06dOxefNm/PzzzwgKCjLZzutUMVEUUVRUxOtz24ABA5CSkoLk5GTpERYWhsceewzJyclo1aoVr1MFioqKcPr0afj5+dV9WTKrey/VyMWLF8WjR4+KCxYsEN3d3cWjR4+KR48eFW/duiWK4p2h0gMGDBCPHDki7ty5U2zevLnJUOmsrCzRx8dHfPTRR8WUlBRx8+bNoqenZ70YKr1q1Srx1KlT4qxZs0Q3NzfxwoULts5anbh165ZUVgCIixcvFo8ePSoNFV+4cKHo5eUlbt68WUxJSREfffTRCoclNm/eXNy5c6d45MgR8YEHHnCooZvPPvus6OXlJSYmJpoM38zPz5f2qe/Xae7cueLevXvF1NRU8fjx4+Irr7wiKhQKcfv27aIo8vpUpvRoI1HkdRJFUZw9e7aYmJgonj9/Xjx48KAYHR0tenh4SJ/JdXmNGLzUgQkTJogAyj12794t7XPx4kVx2LBhoouLi9ioUSNx+vTpJsOiRVEUjx8/Lt5///2iRqMRfX19xZiYGIcdJl3iv//9rxgYGCiq1WqxW7du0hDY+mD37t0VlpsJEyaIolg8NHH+/Pmir6+vqNFoxL59+4opKSkmxygoKBCnT58uNmrUSHRxcRGjo6PFS5cu2eDdWEdF1weAuHr1ammf+n6dnnrqKel/qEmTJuKAAQOkwEUUeX0qUzZ44XUSpXlbVCqV6O/vL44aNUo8efKktL0ur5Egig4+RSsRERE5FPZ5ISIiIllh8EJERESywuCFiIiIZIXBCxEREckKgxciIiKSFQYvREREJCsMXoiIiEhWGLwQERGRrDB4ISIiIllh8EJERESywuCFiIiIZIXBCxEREcnK/wM3kxtI/eZfdQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "bb_lc.plot(rebin_dt = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fddc69d",
   "metadata": {},
   "source": [
    "### Get burst duration"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b3af1b5",
   "metadata": {},
   "source": [
    "We can compute duration values like the T90. We defined it as the time length of the period between which 5% and 95% of the accumulated background-subtracted counts in the identified signal were detected."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "1282ddd8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "229.3461180481858"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bb_lc.duration(quantile = .9)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87923ba0",
   "metadata": {},
   "source": [
    "In order to estimate the error in the burst duration, bctools takes the Bayesian blocks representation of the light curve and the background fit, generates multiple poisson-fluctuated light curve samples, and repeats the duration calculation. For finely binned light curves over a long time span, this can be unfortunately quite slow.  While ideally we would need ~100 samples in orde compute the 68% cont. error, ~10 samples can give you an preliminary order of magnitude estimate in a resonable time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "2f07e710",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 2min, sys: 7.02 s, total: 2min 7s\n",
      "Wall time: 1min 32s\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([-3.21587195,  1.05369836])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "bb_lc.duration_error(.9, nsamples = 10)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:gdt-bc]",
   "language": "python",
   "name": "conda-env-gdt-bc-py"
  },
  "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.11.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}