{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "48171309-57ba-411b-bcd3-fafaf43b029b",
   "metadata": {},
   "source": [
    "# Interface with gammapy\n",
    "\n",
    "This notebook demonstrates how you can easily set up a model for `gammapy` that interfaces with the `OptDepth` class from `ebltable`.\n",
    "For this notebook to run, you need to have `gammapy` installed. I tested the notebook with `gammapy` version 1.0."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c502df4c-bb0f-426f-85bd-c1d5a2b4e232",
   "metadata": {},
   "source": [
    "We start with the necessary `gammapy` imports."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "d6502940-a447-495e-afd6-3aec7414f010",
   "metadata": {},
   "outputs": [],
   "source": [
    "from gammapy.modeling.models import TemplateSpectralModel, SpectralModel\n",
    "from gammapy.modeling import Parameter"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b9809923-4395-42d6-aa99-69fc7e8b33d0",
   "metadata": {},
   "source": [
    "And continue with addtional imports. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "907a72c1-21a2-4bc9-a9af-7737b39340e2",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from astropy import units as u\n",
    "from ebltable.tau_from_model import OptDepth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "801940be-426b-4fef-9a48-d81c7a5906b8",
   "metadata": {},
   "source": [
    "Then we pick an EBL model in initialize the `OptDepth` class. For illustration we also set up a `gammapy` template spectrum. For this, it's however not possible to change the EBL normalization parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bcc4a3e0-10b7-4e9f-acfb-6f36a82508fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "ebl_model = 'finke2022'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1ae78746-3b3f-4509-975f-fde8eddf127b",
   "metadata": {},
   "outputs": [],
   "source": [
    "tau =  OptDepth.readmodel(model=ebl_model)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3c23d882-3997-4573-8728-643c0ef35e86",
   "metadata": {},
   "outputs": [],
   "source": [
    "z = 0.2\n",
    "ETeV = np.logspace(-1,1.5,50)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d035f19f-8397-4dd9-be2d-b69282707b6a",
   "metadata": {},
   "outputs": [],
   "source": [
    "atten = np.exp(-1. * tau.opt_depth(z,ETeV))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "5e74b5a1-91f0-460c-bb24-188a64f25de1",
   "metadata": {},
   "outputs": [],
   "source": [
    "template = TemplateSpectralModel(energy=ETeV * u.TeV, values=atten)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edb43045-9d54-412a-af3e-a12b6295c25b",
   "metadata": {},
   "source": [
    "Now we build our own spectral model class. The class has the possibility to be initiated from a built-in model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "939d0c6a-55d7-4b20-9da9-22ffb9cb1ddb",
   "metadata": {},
   "outputs": [],
   "source": [
    "class OptDepthSpectralModel(SpectralModel):\n",
    "    \n",
    "    r\"\"\"Gamma-ray absorption models build from the `~ebltable` package. \n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "    optdepth : `~ebltable.tau_from_model.OptDepth`\n",
    "        optical depth class\n",
    "    redshift : float\n",
    "        Redshift of the absorption model\n",
    "    alpha_norm: float\n",
    "        Norm of the EBL model\n",
    "    \"\"\"\n",
    "\n",
    "    tag = [\"OptDepthSpectralModel\", \"tau-norm\"]\n",
    "    alpha_norm = Parameter(\"alpha_norm\", 1.0, frozen=True)\n",
    "    redshift = Parameter(\"redshift\", 0.1, frozen=True)\n",
    "\n",
    "    def __init__(self, optdepth, redshift, alpha_norm):\n",
    "        self.optdepth = optdepth\n",
    "        self._cashed_redshift = None\n",
    "        self._cashed_tau = None\n",
    "        \n",
    "        super().__init__(redshift=redshift, alpha_norm=alpha_norm)\n",
    "\n",
    "    @classmethod\n",
    "    def read_model(cls, model, redshift, alpha_norm=1.):\n",
    "        \"\"\"\n",
    "        Initiate the class from a built-in model \n",
    "        \n",
    "        Parameters\n",
    "        ----------\n",
    "        model: str\n",
    "            model name implemented in `ebltable.tau_from_model.OptDepth`\n",
    "        \n",
    "        redshift: float\n",
    "            source redshift\n",
    "        \"\"\"\n",
    "        optdepth = OptDepth.readmodel(model=model)\n",
    "        \n",
    "        return cls(\n",
    "            optdepth=optdepth,\n",
    "            redshift=redshift,\n",
    "            alpha_norm=alpha_norm,\n",
    "        )\n",
    "\n",
    "    def evaluate(self, energy, redshift, alpha_norm):\n",
    "        \"\"\"Evaluate model for energy and parameter value.\"\"\"\n",
    "        if self._cashed_redshift is None \\\n",
    "            or not self._cashed_redshift == redshift:\n",
    "            \n",
    "            try: \n",
    "                self._cashed_tau = self.optdepth.opt_depth(redshift.value, energy.to(u.TeV).value)\n",
    "            except AttributeError:\n",
    "                self._cashed_tau = self.optdepth.opt_depth(redshift, energy.to(u.TeV).value)\n",
    "        \n",
    "        absorption = np.exp(-self._cashed_tau * alpha_norm)\n",
    "        return absorption"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dec1f077-8773-48cd-a359-7920b556574c",
   "metadata": {},
   "source": [
    "Initialize the class from a built-in model is done in the following way:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a1097c21-ed55-4228-8d9b-4b107f0b531d",
   "metadata": {},
   "outputs": [],
   "source": [
    "optdepth = OptDepthSpectralModel.read_model(ebl_model, z)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2849fe99-9b79-42a8-b3fe-c86073a99f9e",
   "metadata": {},
   "source": [
    "Here are the parameters of the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "454f811a-d19f-42c6-90f6-31a05ba4b407",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><i>Table length=2</i>\n",
       "<table id=\"table5563331392\" class=\"table-striped table-bordered table-condensed\">\n",
       "<thead><tr><th>type</th><th>name</th><th>value</th><th>unit</th><th>error</th><th>min</th><th>max</th><th>frozen</th><th>is_norm</th><th>link</th></tr></thead>\n",
       "<thead><tr><th>str8</th><th>str10</th><th>float64</th><th>str1</th><th>int64</th><th>float64</th><th>float64</th><th>bool</th><th>bool</th><th>str1</th></tr></thead>\n",
       "<tr><td>spectral</td><td>alpha_norm</td><td>1.0000e+00</td><td></td><td>0.000e+00</td><td>nan</td><td>nan</td><td>True</td><td>False</td><td></td></tr>\n",
       "<tr><td>spectral</td><td>redshift</td><td>2.0000e-01</td><td></td><td>0.000e+00</td><td>nan</td><td>nan</td><td>True</td><td>False</td><td></td></tr>\n",
       "</table></div>"
      ],
      "text/plain": [
       "<Table length=2>\n",
       "  type      name      value    unit   error     min     max   frozen is_norm link\n",
       "  str8     str10     float64   str1   int64   float64 float64  bool    bool  str1\n",
       "-------- ---------- ---------- ---- --------- ------- ------- ------ ------- ----\n",
       "spectral alpha_norm 1.0000e+00      0.000e+00     nan     nan   True   False     \n",
       "spectral   redshift 2.0000e-01      0.000e+00     nan     nan   True   False     "
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "optdepth.parameters.to_table()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d770bcf-ef60-4fa4-97d9-dac646c479ec",
   "metadata": {},
   "source": [
    "Lastly, we plot everything for comparison. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "4bbd288c-b198-4637-b0e0-8d81d3814c06",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x14b9fd340>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGhCAYAAACZCkVQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB/zElEQVR4nO3dd3gU5drH8e/sJpveG6QQQg+EGpDepAmKgA2xgWJBRUAsqKhgxcN7QGyo4AHUox4UxYoKUgTpVRBCDyQkIY30nt3n/SOwGhIggZDJbu7Pde0lOzM785vdcffOM8/MoymlFEIIIYQQNs6gdwAhhBBCiJogRY0QQggh7IIUNUIIIYSwC1LUCCGEEMIuSFEjhBBCCLsgRY0QQggh7IIUNUIIIYSwCw56B6gtFouFxMREPDw80DRN7zhCCCGEqAKlFDk5OQQHB2MwXLwtpt4UNYmJiYSFhekdQwghhBCXIT4+ntDQ0IsuU2+KGg8PD6DsTfH09NQ5jRBCCCGqIjs7m7CwMOvv+MXUm6Lm3CknT09PKWqEEEIIG1OVriPSUVgIIYQQdkGKGiGEEELYBSlqhBBCCGEXbKqo+fHHH2nZsiXNmzfno48+0juOEEIIIeoQm+koXFpaytSpU1m7di2enp506tSJm266CV9fX72jCSGEEKIOsJmWmm3bttGmTRtCQkLw8PBg2LBh/Prrr3rHEkIIIUQdUWtFzfr16xk+fDjBwcFomsa3335bYZn58+cTERGBs7Mz0dHRbNiwwTovMTGRkJAQ6/PQ0FASEhJqI7oQQgghbECtnX7Ky8ujffv23Hvvvdx8880V5i9dupQpU6Ywf/58evbsyYcffsjQoUM5cOAAjRo1QilV4TUXu2a9qKiIoqIi6/Ps7Oya2ZHzFBSb+WXfcVJPfYKzgwsmRxecTG64mNxwNrni6uSBs7sfLh4N8HR2xcvJBWdH01XJIoQQQtRntVbUDB06lKFDh15w/ty5cxk/fjz3338/APPmzePXX3/l/fffZ9asWYSEhJRrmTl16hRdu3a94PpmzZrFSy+9VHM7cAFpuUW8vmw1QzK+Id8JCkyQ7wT5zlrZf/8xTRnKijAHpXBSitdPF9Cq2IFizcQGN0e+91A0MIdR6jETNycH3ExGCtKewWR0wsXBHRdHN1wcPXB18sTL1R8f30b4BEUS5O5DgKvHJcfEEEIIIexZnegoXFxczM6dO3nmmWfKTR88eDCbNm0C4JprruGvv/4iISEBT09PVqxYwYsvvnjBdT777LNMnTrV+vzcbZZrmoNR4xo/Z+5cbrnksvkmyHWBHBfIcdHwi8on2K0AFGzL9MIv2QMnp0R+KDkNgCOFOEeeLHux+eyjEMg5u8K4v9etlIanxYKrRTEhw4HIEg+KHb045uLEVpciAlya07jFE/i4mvB1c6QkezdBviGEBTXFwcGxJt8SIYQQQhd1oqhJS0vDbDYTFBRUbnpQUBCnT5f9wDs4ODBnzhz69++PxWLh6aefxs/P74LrdHJywsnJ6armBmjo5cJb991ISvphLLk5mHNyseTmYsnJwZxb9m919jSYa3HZIzALQKHun83B8BBKC/MJ/eE7Jn+/kdSOobSb1Ia8IjMFOan0flmR56qR6wbZrpDlpshwV6R6Qpy3kXgvKHBRaJoix6iRY9RoaDlNm+ITUAyxmhu/ePjROjOJ+csGWXMHtJxGoUHDqBTeZoWXxYC7xYSH5oqTqQFm/140cA8g1COQ5o6KVo3a4ufd4Kq/n0IIIcTlqhNFzTnn95FRSpWbduONN3LjjTfWdqxLcvD1JXjW6xecbykuxpKbizkrC3Nm5tlHFh7X9sfo5QVAVoZGZrKZNn360Kd7YwBKkj04mmPGPQeCKl8zAJq7O5agAEq8nMn3dMLQrzs7PVwpzTtDfu4RBpScwMkxGJeWAWTkFZOdl00BCtAwaxrpDhrpABSffWRC+kHOTiyzG2ak5BBd4EW2UyBH3bzY5WYhxKst7do+RkNvZ0J9XHByMF7ZmymEEEJcpjpR1Pj7+2M0Gq2tMuekpKRUaL2xRQaTCYOvLw4XuaeO1/Ab8Bp+Q7lpDj7eNP7qK0rTUilNTaU0La3svymplCQlUpqQiDkrC5Wbi5abiwkwAY0nPo9LVBsAmn/yKV0XHsL71s4ETLoGKCsWi09EUurvy6m045w8fZDTGbGk5SSQWZhCgjIQa/Agz5xBsSUDkyGZfIOBcHM+EZYMKDjBYQc3vvfwo/XpE8zeVrZeTYPWTabhphzwxQM/R3+C3MMJD2hNZMQ1NAtti8EoRY8QQoiro04UNSaTiejoaFatWsWoUaOs01etWsWIESN0TKYvzWTCpW3URZex5OVRkpRESWLi2UcSpkZ/9x0qPhFLaWoqqrTUOq00JYXjQ4eBgwOmxuG0aRVJdGQrnFr2xTmyFQ6VnNZLSo0jNyWRP1MTKEqPIzdrL72LjmHSgmkW6E5iZgHG0tPEmTTKOv9klj3yj8LJ1XASHC2ACsXbIZgQt3C6kEvL4A50iRyAr7ftF69CCCH0panKrpW+CnJzczl69CgAHTt2ZO7cufTv3x9fX18aNWrE0qVLufvuu/nggw/o3r07CxYsYOHChezfv5/w8PAr3n52djZeXl5kZWXh6el5xeuzFeacHIpjYzF6e2Nq1AiAgr17OXnPWFRhYaWvcQgIwKlVK5xbtcI5shVOkZGYGje+6CX0SilOZ2awY9+3xKXGcDonlvTiFM6QQ5qhmFQHDcsFXr8wKZmIQndSncLY5xXISS8PosIH0zv6ZrxcpBOzEELUZ9X5/a61ombdunX079+/wvSxY8eyZMkSoOzme7NnzyYpKYmoqCjefPNN+vTpUyPbr69FzYUoi4XS06cpOnKEwpiDFB48SNHBgxSfPAmVHBJGLy+a/PIzDj4+1tdr1biEPC8/hz8TDrEjNYGY9GOczjyAY8EOkhxK+CohiQZmMwBv+XjxkbcXHbI82ZD4HA29nGkdCD78mzDPFrQO60G3tkPwcPOukfdBCCFE3VYnixq9SVFTNZa8PAoPH6bo0KGzxU4MRTEHMXp50Wz979bWmlOPTaL45EkCpz2Ne8+eV7TNrPRkko7vI/vUAf5IW88uywlMeVGsSbkVgObOWzkdsdy6vFEpQko0QpUXYS5NaBkxhH4dhhPgLp+rEELYGylqKiFFzeVTxcWUJCVhOnsaUCnFkd59MKelEf7557h26ghA7u+/k73iZ9x6dMe1W3ccgwKvaLtZBSUcTs7hz5hV/HnqQ06rMyQ4lJBjrNhCpJSGgzmQhg7BdNdKadWwGwOjR+PrK311hBDClklRUwkpampWSUoKBX/+iXufPhjO3g8o6cUZZH75pXUZU9OmuHXvXlbkdOmC0cPjirdrMZvZf3wH2w/+ytG03ZwqOcV+o5Fih7xyyzkqxebYUyQZQ0nxaseRho0IadKVHu2GYTJd/fsXCSGEqBlS1FRCipqrL3/XbnLXriVv82YK9+8v3zfHaMQlKgr3fn1xv3YATi2aX7TjcXUdSk3kt2M7ORL7C5nZ2zGZc1mQ8vewGrcFNyDGyUTb5G5ovuPpGOZDa79MIht6ERHSqsZyCCGEqFlS1FRCipraZc7MJG/bNvI2byZ/0+ayDsj/4BgaiseAa/EaORLnyMirkiHtdDyn9m0gL3Yz76k/OGoqhmOPkFLaGIDeAe+xxz+e0BJFhMWPFt7t6NlmJNGt+sn9dIQQoo6QoqYSUtToqyQxkdw//iB3zVryNm1CFRcD0GDGi/iMGQOUdVIGMLi5XZUMxcWFHEsvZndcJrvjMjiT8hjb3DMqLOdbaqGp2R2voJH0a34tQ5p3lJHVhRBCJ1LUVEKKmrrDkp9P7saN5K5eQ8CUyTg2KBtTKuOLL0ie9QY+d91F0NNP1UqW+KQj/LbzfxxI3kysOYHjJjMl5w/XYXEiwuxJR0d/ejYdzrXX3Iaj9MsRQohaUZ3f7zpxR2FRvxhcXfEcNAjPQYPKTS/Yuw9VXIyD39/DSZhz8yjYtRO37t3RHGv+RnxhDZtz7w0vWJ9n5Z7ht21L2Za4nY2lhWRZjqAZCzlhSOUEqfjt3ELvVU8Q49KW9IadyYlozsBrxuDqfHVal4QQQlSdtNSIOkMpRdHhIzgE+FvHycr6/nsSn56G0ccHz6HX4Xn99bh07FitG/9diVKzmZVHdrN7xzscyYthbEYy/QszAVjv4syjDQIJLwJfh//QvYkfPZr60bqBO46O8veCEELUBDn9VAkpamxTxv+WkvrOO5jT/x4y3CG4IV43DMf7lputQz/UFovZzImYHaTs+41taStZ6pZGk1wvNiROB8BIMeHNnyfU7EoLl1b0aDmCPp1G4OAgwz0IIcTlkKKmElLU2C5VWkrelq1k//gjOatWWTsUA7h274b3LbfgMWgQBlPtd+YtLS1hf9xJdiUZ2XwsndSE7zga/G25ZbzMFlqUeNDSI4o+7cbQtXV/DLXU0iSEELZOippKSFFjHyyFheSuW0fmsq/J27jRei8co7c3XiNG4H3brTg1bapbvuLiItbsWMbWoz9wuOAIh01FFBrKdzzWzF4EOUbRLbA9NzduS4eWvXRKK4QQdZ8UNZWQosb+FJ9KIOubb8j85htKT5+2Tg999x08Bg7UMdnf8vJzWLX1c7af+JUjJcc56KhQBku5ZUJKLLyY1xZTk/40jr4O/+ArH5VeCCHshRQ1lZCixn4ps5ncDRvI/GoZ+du303ztGuu9bgr27sXBzw/HkBCdU5bJLMjj6/1/sPrEZtIzfyPFIZPIomI+T0q2LvNYQCjKFEiflk/Tv+O1BHjI5eNCiPpLippKSFFTP1jy8qwFjVKK2JGjKDp6lNC35tWZ1pt/Sk5PYPfOH/A+vhf/1C0EmWPp3TgEs6bhfuQBkkqb0jzQnZ7+6wnyyGdI57toEtZG79hCCFFr5D41ot76592ILdnZGH180JyccO3SxTq9JDERh4CAq3Lfm+oK8gvhusETrM9TUuKZuOUjDqfuZl9AR5KSsjmSkouv83K+ViXMX/MjjYqhGYG08r+GAZ1G0yK8g347IIQQdYi01Ai7V5qaikNAgPV57OjRlJ5OxueOO/C+7VYcfHx0THdxGXnFbI09w29bxnNQxXGykgu8gkvAwaUfnRt0YUSrXnQKaVL7QYUQ4iqR00+VkKJGQFmBc3zUTZjT0gDQnJ3xvvlm/O67t870u7mY4/H7WbXzMw6kbuU4KZx0VKjzhnXwKzXQwuzOtV596Nb+Vho1by8DdAohbJYUNZWQokacYykuJnvFCs58/AlFMTFlE41GvG64Hr/778epeXN9A1ZD/Omj/PjXatacOUls7j6KDCdBK/tf+t3TKfQtKCQbN37zaM4uX386hl/HkJ734u4kZ56FELZBippKSFEjzqeUIn/LFtIWLCB/8xbrdPdrr8Xvgftx7dhRx3SXJzUrlVVbPmZvwnpuSk6nbeFhXLRi3vbxYqG3F+2z3diU+AItG3jSKdSFIMt/6dLyOrq0HiB3PRZC1ElS1FRCihpxMQX79pG+8CNyVq2y3tDPtXNn/B56ELdevdDOO8VjK0qKizhxYBu/7f+Cnbm7seRFsSZlJADNnHeQHLEMAHezhcalJkKNwTT370SPNjfQOqKznLYSQuhOippKSFEjqqLoeCzpi/5D1nffQ0kJAP6PPELApMd0TlZzkrML2XUyg137PmNX7qecdDRTZKhYtLmUmnByaE1Tz1Z0DWzJtcFNaRHevtYGExVCCJCiplJS1IjqKElO5sziJWQuW0bjr77EKSICAHNODgY3N7v6Yc8vzOOP3d+z58RaTuQcIl47Q7yjwlxJ61RAqYU3UwPJ94nCMbQdjuEtiGzRDZNJbhAohLg6pKiphBQ14nJYCgowuLhYnyc88SRFhw/RYOZMXKOjdUx2dWVkpbIydi8bTv3FoYwDFBX9SZZDLs2LS1iW+PeQFKODgzjuaCI6bzi+IXfQMsiDRu5naOznSkRIKx33QAhhL+Tme0LUkH8WNOacHHL/+ANLVla56fbIxyuA0R0GMLrDAOu0jKxU/or5na0+J+D0PjyzYjjlkEehQWNfSgAJp+IA6B3wHnv84/EptRBiNhFk8CfUvSktgrvQtc11BPnV/UvnhRC2SVpqhKgGc3Y2ub+vx2v4DdZpaQsWYgoLxWPIELs6LVUVxcVFbNv/G6kqiiMpBRxOzkHlPMU2j/QLvsah1AMXQyOCXBrR1tmN9p4BRLccQOOQlrWYXAhhK+T0UyWkqBFXQ0liIkeHXAclJTg1b4b/Y4/hMWiQzV4tVVNSMxLZtv83DiVs5VTOUU5bUkk0FpHuUHnR1ye/gFeSizhtCifbownrfAw0Coiid/S9hPr51vv3U4j6TIqaSkhRI64Gc04OZz75hDNLPsaSkwOAc+vWBEyZjFvv3vJjfJ745JNsTD7JntOHOJJxDEP2H2RoGQzLy2VqRiYA6QYD/cJD0ZSi+NDzOJp8aBroTmvnj3FyzCbcrw1RjXvToWVPubeOEPWAFDWVkKJGXE3m7GzSFy8m4+NPsOTnA+DSsSMBkyfj1q2rzunqvrycTJKO/0XmyX0knd7F1+pPClQJO+Jew2wp+4rq1Pg5jrhYrK9xtihCSg00UF4EuzQiqMUd9GrcjsiAEAz17DSgEPZMippKSFEjakPpmTOkf/QfMj77DFVUBIBrt24ETJqEayfbu0Ox3opKzZxMz+dIci7bdz9JfNFxUrRcEh0UxZXcWwcAswuNzUYa40xn9070aDWCRpGdcXJ2rd3wQogaIUVNJaSoEbWpJCWF9A8XkPHll9ab+Ln17UPApEm4tGmjczrbV1iUz66Dv7Mv9g9OZhwg1pxNjKZRakxB0/7+SpuWnsFd2TkUKyNbnMP5n58nTT3bMrDnq7Rq6IGTg9wxWYi6ToqaSkhRI/RQkpBA2gcfkPnNcjCbAQj//DNcO3XSOZl9yirMZ+OJAxzb9xlx2Yfom5FLz9zj+JDDd+5uPB/gR4sCAztPvI6DQaNFkAeRrm8Q7t2cfh1up13zbnrvghDiPFLUVEKKGqGn4pMnSX3vPYqPHafxV19aL/225OVhcHPTOZ19UxYLp+OPsHXf92xK/p3iYl9+T72LjPwSPAxnoOVs67JBJRaaWXxo4dWOXq1vonPr/jL+lRA6k6KmElLUiLrAUlyMwWQq+3deHkevuw73Hj0Ieu45jF5eOqerP5RSJGQWsPPwXtbvm8lxSyKxJkuFoSE8Sh1wMXWmQ0An7mw7hE4hTXRKLET9JUVNJaSoEXVN9i+/kDDlcRwbNaLpTz+iOcrlyXpKTk9g1fbP2Ju4ntjSUxwzlVJyXpHTuBgijeHc0WYCbaOvw+ggN2UX4mqToqYSUtSIuqhg3z4subm4de8OgCouJvW9+fjcMQbHoCCd09VvmTlpfBezjTWndhOTuYNCw1GUBp5mM7/HJZCNJ8e8e5LYtBPdetwpwz8IcZVIUVMJKWqELcj43/84PfMlNJMJ79tuw++BB3AMCtQ7lgBOJh7m+43vk5eyh0dSDuKp5aOA60KDSXUw0qtwDL2i72FYu4a4O0kLjhA1RYqaSkhRI2xB/q5dpM59k/wdOwDQnJzwHn0bfvffj2OgFDd1RUlxEYe2/Uri3uXMMW0n3Qi5h5+nUHng7GhgVOO1tA0N4Jb+kzCZnPSOK4RNk6KmElLUCFuhlCJ/yxZS33mXgl27gLLixuf22/F74H4c/P11Tij+yWI2s+vwdrae9ufrXac4nppDq6bPkmDS8C+10JnGjOj0ML063nDplQkhKpCiphJS1Ahbo5Qib9Mm0t55l4I9ewDQnJ3xGTMGv/vH4+Dnp29AUYFSiu1HT/LpunvZaUwhx/j3cA3NijSucYvmrgEvENZArqISoqqkqKmEFDXCVimlyPtjI6nvvkPhn3sB0Fxc8Bk9Gt/77pXTUnVUTl4mn6+azcaUVew1FVgvFzdZNBq5XMfzvR8iOqSpzimFqPukqKmEFDXC1imlyNuwgdR33qVw3z4AjP7+NF+zGu3svW9E3XQ07i/+9/u/2FK4h5NnPyqlDIRq7ZnaYiiDu4/RN6AQdZgUNZWQokbYi7KWmz9Im/8+rt26Ejh5snV6aWIijiFyaXFdZTGbWbjjV5YcWEKuIcY6vVuBiYmtnqNdn1HWu00LIcpIUVMJKWqEvVFKQWmp9aZ9uRs3Ev/Ag3jffBMNX3lF53TiUn6I2c7yDVPZZcrg/sxsJmZmEWtoTGrbB2gzZCxurh56RxSiTqjO77f8SSCEjdI0rdxdiPO3bQeLBc3FRcdUoqqGR3Zh0YMbWNL9A1qZBpCvnIiwnKDw8Cvc+EU33v5qKmZLvfibU4gaIy01QtiRgv37cQgIsHYezt+xg5R58/C7917c+/eXUxt1WNaZVA78MI/P8r7mdzdHOmb4kKxe5akhLRkQGYh23pANQtQXcvqpElLUiPro1GOTyFm1CgBTeDi+947Da8QIDNKaU2dl5qQxb/lE1pwYQlx+WXE6KOwo17d2ZmT/h3ROJ0Ttk6KmElLUiPqoJDmZjP/+l4ylX2LJzgbA6O2Nzx1j8LnjDrmRXx2WlV/C+78fY/HGWNqFPMdB11I6FjlzX6fn6Nd5lN7xhKg1UtRUQooaUZ9Z8vLI/PobznzyCSWnTgGgmUx43jgc37vvwbllC50TiguJT0vnjW9uZqMpDbOmoSlFj2JvJvR9gw4te+kdT4irToqaSkhRIwSo0lJyflvNmcWLKfjzT+t0l87R+IwZg+egQXLPmzpq058/89HWmWx3ygfAZFEMpSXPjl4iV0oJuyZFTSWkqBGivPxduzmzZAk5q1eD2QyAQ2AgTX/9Rfrc1GEr/viETw68yX6nUgC8S9yZ0GUWd7bvp28wIa4SKWoqIUWNEJUrSU4mc+mXZHz1JS7t2hP23rvWeYUHDuDUqpVcNVXHWMxm5n87jc+zVpJjLPsKD3Pox4fXzyDMW/pJCfsiRU0lpKgR4uJUSQnmrCxr5+HiU6c4NmgwpsaNifh6GQZXV50TivPFZ6YxYcXLxJWsBcDTrHG750AeHfV/GIxGndMJUTPk5ntCiGrTHB3LXQ1VdOgQBldXHBs2LFfQ5O/ajSop0SOiOE+Ytz8/3fE2z3R4B58Sd7KNigV5q1jxZn8Sju/XO54QtU5aaoQQF2TOzcOcmYEpNBQoO1V1tP+1GL298Ro+HK+bRuHcsqXOKQVAXl42b3x5L845u5l+Jp1C5cjuJg/RecyLOJqc9I4nxGWT00+VkKJGiCuXt2UrCU89iTk1zTrNuXVrvEaNwvOG63Hw8dExnQCIP7qPrK8mElW0hwQHI9P8w5jcey5d2gzQO5oQl0WKmkpIUSNEzVClpeT+8QdZ3ywnZ+1aOHcqytER9x498Lx+GO7XXovR3V3foPWYsljY8f37fJT0NptcTbTOd+DWa37g5k4hMtyCsDlS1FRCihohal5pRgbZP/5E5vJvKDoQY52uOTnh3qdPWYHTt69cIq6TvYc3MXvtJI4n3EVicXOub9eQ10e2xcvV8dIvFqKOsMuiJj4+nrvvvpuUlBQcHBx44YUXuPXWW6v8eilqhLi6io4dI3vFz2SvWEFxbKx1utdNNxH8+ms6JqvfzBbFB78f481Vhym1KAY1/IDBbYZw68CJekcTokrssqhJSkoiOTmZDh06kJKSQqdOnTh06BBubm5Ver0UNULUDqUURQcPkr1iBdk/raDBjBdx79sXgMKYGM4sWYLn9dfj3qePzknrlz/jM5nz1Rvs9vsJTSmuM4cz884vcXWu2neoEHqxy0u6GzZsSIcOHQAIDAzE19eXM2fO6BtKCFGBpmk4R0YS+MQTNF39G269e1vnZf/8C1nffU/m19+Ue40lP7+2Y9Y77cO8mXPvJHoX+aA0jZ8d4rjj0+5s2bdS72hC1JgaK2rWr1/P8OHDCQ4ORtM0vv322wrLzJ8/n4iICJydnYmOjmbDhg2Xta0dO3ZgsVgICwu7wtRCiKtJ07RydyP2GDQI37H34DVqpHVa0fHjHO7WnfhHHiXz228xZ2XpkLR+CPAJZv6D65noNRwPs4VjJsVjO6by1leT9Y4mRI1wqKkV5eXl0b59e+69915uvvnmCvOXLl3KlClTmD9/Pj179uTDDz9k6NChHDhwgEaNGgEQHR1NUVFRhdeuXLmS4OBgANLT07nnnnv46KOPaiq6EKKWuLSNwqVtVLlpeX9sRBUXk7tmDblr1pDk4IBb1654DB6Mx8ABOPj56ZTWfj008nV6Hh3Bq2smsN+plI/y13D4o4G8cfcPeDhJp25hu65KnxpN01i+fDkjR460TuvatSudOnXi/ffft06LjIxk5MiRzJo1q0rrLSoqYtCgQTzwwAPcfffdl1z2nwVSdnY2YWFh0qdGiDpGKUXR4SPkrFxJzsqVFB058vdMgwHXzp3LCpxBA3EMCtIvqB0qLS3hxU9v40ftCErTcDJHsGjY27Rr0FjvaEJY1bk+NcXFxezcuZPBgweXmz548GA2bdpUpXUopRg3bhzXXnvtJQsagFmzZuHl5WV9yKkqIeomTdNwbtmCgMcm0uSH72ny8woCHn8c5zZtwGIhf9s2kl99laN9+3Hi9jGkL1pM8akEvWPbBQcHR16/dzmPhU0EswtFxljuXDGGJdt+1DuaEJelVoqatLQ0zGYzQef9lRUUFMTp06ertI6NGzeydOlSvv32Wzp06ECHDh3Yt2/fBZd/9tlnycrKsj7i4+OvaB+EELXDKSIC/4ceJOLrZTT97TcCp03DpWNHAAr27CFl9mzOfPyxzintywMDJvDRwM9wLA0FYy5vHniGmR/fjsVs1juaENVSY31qquL8O1kqpap8d8tevXphsViqvC0nJyecnGS8EyFsmSk0BL97x+F37zhKklPI+W0VOStX4Tnk71bf/O3bSXnrLbxvvgXvf3RAFtXTtVFz1tzxFZM+v5XdDqf40bKPrm8No9f4L/Dw8tU7nhBVUistNf7+/hiNxgqtMikpKRVab4QQojKOQYH43nkn4R8vwbVzZ+v0rJ9+omDHTvK3b7dOU0phzs7WI6ZN83ZxZ8m4H7nbsRvT0zIZmr2FjLd6cTJmp97RhKiSWilqTCYT0dHRrFq1qtz0VatW0aNHj9qIIISwU/4TJhA4bRo+t/19h/HCv/7iSM9enHrsMXI3bEDJaZQqMxiNPH3HQlr3/5Rk/GhkSSBl+XA+XP6s3tGEuKQaO/2Um5vL0aNHrc9jY2PZs2cPvr6+NGrUiKlTp3L33XfTuXNnunfvzoIFC4iLi2PChAk1FUEIUQ85NmiA373jyk3L27IFVVJCzqrfyFn1G47BwXjfditeo27CMShQn6A2pmXna0kPW8/mxWN4ISCN5Owfyfy8iKfHzJFBMUWdVWOXdK9bt47+/ftXmD527FiWLFkClN18b/bs2SQlJREVFcWbb75Jn1q6VboMkyBE/VJ46DCZXy8j69vvsJw7FWU04nFtf7xvuw23Hj3QjEZ9Q9qAwqJ8nvl0BIdUEqdin+b66A68MjIKR6PN3JBe2Di7HPvpSklRI0T9ZCksJOfXX8lY+iUFu3ZZpzsGB+N96y143XSztN5UwX/Wbue1lSlYFPRs6su8W5sR4B2gdyxRD0hRUwkpaoQQRUeOkPHlV2R994/WG0dHvEeOwP/hh3E8e+dyUbnVMck89sVuOnl8QIb3MV7r8z4dWvbSO5awc1LUVEKKGiHEOdbWm/8tpWD3bgAaL1uGS1QbnZPVfTuPHeepdcNJdTDgU2rhmZbTGNbrHr1jCTsmRU0lpKgRQlQmf9cu8jZtJmDio9ZpmV9/jXNUW5xbttAxWd2198gWpq97gBMmcLIoJvjdxP03vqx3LGGnpKiphBQ1QoiqKElO5tjAQaiSEiK+/w7nFlLYVCY5PYGpy25kr3MxmlLcauzA9Ds+xiCdr0UNq3NjPwkhhK1QJaW4DxyAa5cu5QoaS16ejqnqniC/EP4zbgN9i/1RmsaXlj+ZsmgwxSUlekcT9Zi01AghRCVUcTGayQRAaUYGx28YjtcNN+D/2ESM7u46p6s7LGYzr35+D8vMf6I0jSBDd76//R1cHWWYGlEzpKVGCCGu0LmCBiDnl18wp6dz5uOPOTZ0KFk//EA9+XvwkgxGIy/e/Rljg+5DKQPJls1c9/kDZBXm6x1N1ENS1AghxCX4jBlD2MIFOIY3wpyaRuJTTxN39z0UHjqsd7Q644mhU3mo5csoiwMZ7ObeTweSlpGkdyxRz0hRI4QQVeDeuzdNfviBgCmT0Zydyd+xg9ibbiJ51izMOTl6x6sTHus+gifazcbJAkdMOUz46jrOnEnWO5aoR6SoEUKIKjKYTPhPmEDTn37EY9AgMJs58/EnHBs2jKzvv5dTUsC90YOY0fRxPMwWbstJI+O960hPPqV3LFFPSEdhIYS4TLkb/iD51VcpPnkSALe+fQh+4w0cfHx0Tqa/PXt+I/Tb8fiTSZwhBKf7fiAotKnesYQNko7CQghRC9x79yLih+8JePxxNJOJvN/XEztiJPnbt+sdTXcdOgyk4M7vOY0/HiTx3A/D2b5/td6xhJ2TokYIIa6AwWTC/6EHafzVl5iaNKE0JYXiUwl6x6oTwpq3h/t+4Xn/YLa5Gnl14+McSsrSO5awY1LUCCFEDXBu2ZKIr76kwSsv4z1qpHV6PTnDf0ENGjVn0vWf0brQgayEMYxeuJWYpGy9Ywk7JUWNEELUEIObGz633mp9XpqRwYmbbyH3j406ptJfy/AOfHj3VrwCepOZX8JdH23lUGKa3rGEHZKiRgghrpL0Dz6g8MABkmfNQpWW6h1HV96uJj4Z35U2wZ6EaD8yaUU/6WMjapwUNUIIcZUEPP44PnfeScicf6M5OOgdR3deLo58fG80hsC1JDpqPLNlEnsPb9I7lrAjUtQIIcRVYnB2psELz+PcqpV1WvbPP1N48KCOqfTl7+HK7CGfEFyiSHEw8OT6B4k5vlPvWMJOSFEjhBC1JG/LVhKeepqTd91dry/7bhHegX/3+4igEgtJjhpPrB7L8fj9escSdkCKGiGEqCXOrSNx6dAeS24ucePvJ2d1/e1T0rZZN2b3mI9/qYV4k8bkX24nPumI3rGEjZOiRgghaonR05NGH32E+7XXooqLOfXYJDK//lrvWLrp1Lovr3eZg0+phRMmuP/Huzmdk6F3LGHDpKgRQohaZHB2JvTtt/C66SawWEia/jxpCxfW2/vZdG93HS+3fxUXs5FEUx43fnUf6fkyQKi4PFLUCCFELdMcHGj42qv4PXA/AKlz5pLyr9koi0XnZPro13kUL3RfAGZnCoxHuWHpfWTmnNE7lrBBUtQIIYQONE0j8IknCJw2DYAzS5aQ9OyzqJISnZPpY3jkNbzQ5U2wmMg1HOTJL67HYjbrHUvYGClqhBBCR373jiP4X2+A0UjWd98TP3EiloICvWPp4ra2vXgm4mEclGKrUy7zPrlH70jCxkhRI4QQOvMaMYLQ995Fc3Ym7/f1xI2/H0tent6xdHFnv/u507kXd2blMOnkj+xc8R+9IwkbIkWNEELUAR79+tFo0SIMnp4U7NpF1vff6x1JN0/e/gH9HK/DAWizdRqHd63TO5KwEVLUCCFEHeHaqSONFi7Af9JjeN9+u95xdNXlgXf506UrRq2ELzc8yF9Ht+odSdgATdWT6wizs7Px8vIiKysLT09PveMIIUSVqJISMBrRDPXvb9CcrDO88mlffvaApkUan961BQ9XV71jiVpWnd/v+vd/iRBC2AhLQQHxjzxK8qw36uV9bDy8fLl94Hwalii01F48+fUBLJb69z6IqpOiRggh6qj8bdvI27CBzGXLKImL0zuOLjq16s0r/X5lf/5wft2fzL9XHtI7kqjDpKgRQog6yr1vXxrMnEmj//wHU3i43nF007VZCG/c3BaAbzetYP43z+ucSNRVDnoHEEIIcWE+t48u99ySn4+hHvYrualTKEdPbOT7zP/wUbZG8NowRvZ/SO9Yoo6RlhohhLARhQcOcGzIdWSvXKl3FF1MvfEmmpe4U6JpzIl9m10HN+gdSdQxUtQIIYSNyPz2W0pTU0l84klyN27UO06tc3Bw5N+jv6NJMWQaDby6/lEKi/L1jiXqEClqhBDCRgRNm4bHkCGokhJOTXyMwpgYvSPVOj/vBrzS90PczRaOOCle+fxOvSOJOkSKGiGEsBGa0Ujw/83GrUd3VEEBpyY+RmlGht6xal27Fj24w2MwAD9pR/h546c6JxJ1hRQ1QghhQwwmEyFvvoljWBglCQkkPvEEqrRU71i17tGb/k2XQjfMmsacmLlkFtTPsbJEeVLUCCGEjTF6eRH67rtoLi7kbdpMyptv6h2p1hmMRmYM/wS3UkeSHUsZ991MvSOJOkCKGiGEsEHOLVsQ/PprAJz5zyKyV6zQOVHtCw9uwd1RLwNwrOgXFm7/VedEQm9S1AghhI3yHDoUv/vHA5A4/XkKD9W/u+0+2vUGGpsGAvDR3unEnz6ucyKhJylqhBDChgU8/jhuPXpYOw6bMzP1jlTrlox4Gf8SR/Idinj5+7v0jiN0JEWNEELYMM1oJGTuHBxDQymJjyfhyadQZrPesWqVn6sHj7d4FKNSBBYnsfWnBXpHEjqRokYIIWyc0dub0HffQXN2RnN2QhUX6x2p1t3YZzwvlvbgtbQztNo+k9TEE3pHEjrQVD0Zzz47OxsvLy+ysrLw9PTUO44QQtS4wsOHcWrWDM1QP/9eLSku4sS/etDcfJQ/naNp+9QqDEaj3rHEFarO73f9PPKFEMIOObdoYS1olFKUnjmjc6La5WhywnTrQk4ZnPjQK47Z/3tQ70iilklRI4QQdsZSUEDik09x4vYxmLOy9I5Tq8JbdeLTiKFscHXhh6ItHDh5Qu9IohZJUSOEEHbGUlhIwe7dlCQkkLdlq95xat1TYz6iV34Ajifv5MUVidSTXhYC6VMjhBB2qWD/figtxaV9e72j6CIxs4ABc36noMTMe3d04vp2DfWOJC6T9KkRQoh6zqVNm3pb0AAEe7vwUN8mACxZ+Tk5eZn6BhK1QooaIYSwc0VHjpA0c2a9u3/Ng32aMKDhOxz2X8zcrx/WO46oBVLUCCGEHbMUFHBy7Dgy/7eU9EWL9I5Tq1xNDrRp0BSLpvFL6V5iEw7qHUlcZVLUCCGEHTO4uBD4xBMApL79DoUHDuicqHZNHDWXxsWQazQwb8VEveOIq0yKGiGEsHNeN43CY9AgKCkh4amnsRQW6h2p1phMTtwZfi8A6x1Ps3nvLzonEleTFDVCCGHnNE2jwcsv4RAQQPGxY6T8e47ekWrV7YOn0qHQiVJN44PNL+gdR1xFUtQIIUQ94ODjQ8PXXwcg47//JXfDBp0T1a6Hus7EqBS7nAv56rd39Y4jrhIpaoQQop5w790Ln7vuAiDxuecozcjQOVHt6dXhBnqXBALw7okvKa1nV4LVF1LUCCFEPRL45BOYmjbFnJrG6RdfrFd325089D2wmDjjmMHLa/+rdxxxFdhcUZOfn094eDhPPvmk3lGEEMLmGJydCfm/2eDoSM6q38j6ZrnekWpNs9BIuvneDsC3Jz8iIz9X50SiptlcUfPaa6/RtWtXvWMIIYTNcm7dmoBJjwGQ/NprFMfF6Zyo9vx7yKNopT4oh0zekBvy2R2bKmqOHDnCwYMHGTZsmN5RhBDCpvnddx+unTtjyc8n6bnpKItF70i1wsvZlfv8+wOwtnQXh0/u0TeQqFE1VtSsX7+e4cOHExwcjKZpfPvttxWWmT9/PhERETg7OxMdHc2Gava+f/LJJ5k1a1YNJRZCiPpLMxpp+MYbGNzdcW7TBlVaqnekWjPphhdpXqTRoriEgz+8qnccUYMcampFeXl5tG/fnnvvvZebb765wvylS5cyZcoU5s+fT8+ePfnwww8ZOnQoBw4coFGjRgBER0dTVFRU4bUrV65k+/bttGjRghYtWrBp06ZL5ikqKiq3ruzs7CvYOyGEsD+m0BCarlqJg4+P3lFqlcFo5Jk2M+my6l4s6jdiD2wnonUXvWOJGqCpq9D1XdM0li9fzsiRI63TunbtSqdOnXj//fet0yIjIxk5cmSVWl+effZZ/vvf/2I0GsnNzaWkpIQnnniCF198sdLlZ86cyUsvvVRhelWGLhdCiPpIWSxoBpvqlXBFdv3fcDrlrWevc2faPbNa7zjiArKzs/Hy8qrS73etHL3FxcXs3LmTwYMHl5s+ePDgKrW6AMyaNYv4+HhOnDjBv//9bx544IELFjRQVgRlZWVZH/Hx8Ve0D0IIYc+Kjh3j5F13k71qld5Rak3QzW9w2uDIFqfD/LZlqd5xRA2osdNPF5OWlobZbCYoKKjc9KCgIE6fPn1Vtunk5ISTk9NVWbcQQtibrB9+oGDXLlIzMvAYMKBetNiENGnDg0HN2eycS++97zGw22i9I4krVCtFzTmappV7rpSqMK0qxo0bV0OJhBBCAPg//DDmzEz8H3ywXhQ05wxtMY6Uw++RkdmUrPwSvFwd9Y4krkCtHLn+/v4YjcYKrTIpKSkVWm+EEELUPoOTEw1nzsQxOFjvKLVqZL8HKSx4l80Zt/DVTummYOtqpagxmUxER0ez6rxztatWraJHjx61EUEIIUQ15G3bhqWSq1HtjaZp3NO9MQCfbjmJxVJ/ho2wRzVW1OTm5rJnzx727NkDQGxsLHv27CHu7J0qp06dykcffcSiRYuIiYnh8ccfJy4ujgkTJtRUBCGEEDUg5d//Ju6esaR98IHeUWrFyI7BBLtmEMabLJURvG1ajfWp2bFjB/3797c+nzp1KgBjx45lyZIljB49mvT0dF5++WWSkpKIiopixYoVhIeH11QEIYQQNcC5XTsA0hd+hOd1Q3Fu2ULnRFeXq8mBbsGfsMoxCS32Y8bwmN6RxGW6KvepqYuqc527EELUd/ETJ5L722qc27ej8eefoxmNeke6qjbv/ZUHdz+JphRLrnmPTq376h1JnFXn7lMjhBDCtjR44QUM7u4U/rmXjM8+1zvOVde93RCiihxRmsbnm9/QO464TFLUCCGEqMAxKIjAJ58EIGXePEoSEnROdPUNaDAcgC1aHFm5Z3ROIy6HFDVCCCEq5X3brbh0jkbl53P65Vf0jnPV3XPdcwSWWsgyGli8YobeccRlkKJGCCFEpTSDgYYvvwyOjuT+/ju5v/+ud6SrymRyoqdDawB+z7TvfbVXUtQIIYS4IKcmTfC9+24Akme9gSou1jnR1TXu2pk4KsVRJ8Xyvev0jiOqSYoaIYQQF+X/yMMY/fwoPnGCM//9TO84V1WTsDYEGroB8MHeL3VOI6pLihohhBAXZXR3J3Dq4wCkzZ9PaVqazomurvs6jgMgoWQzx88k6xtGVIsUNUIIIS7Ja9QonKOisOTmkvrWW3rHuapuadMDkzkMzVDK+ytf0zuOqAYpaoQQQlySZjAQNP053Pv1w2/8eL3jXFUGg4FbvaIA2J27muJi+x8Dy15IUSOEEKJKXDt2JOyD9zE1bqx3lKvuoYFP4WG2YNYUa9b8R+84ooqkqBFCCHFZzLl5eke4any8Ani8uBMr4xMI3vON3nFEFUlRI4QQolrMOTkkPv88x4cPx5Jnv4VNj8HPYFQa7Qq3c+roX3rHEVUgRY0QQohq0Uwm8rduozQpidz16/WOc9WENGnDPpfOmIFtv83WO46oAge9AwghhLAtBicnGr7yCprJhGunjnrHuapSOt3C0PgEcrWd9M9KxccrQO9I4iKkpUYIIUS1uXXravcFDUDffvdiUQYsGvyw+Su944hLkKJGCCHEFSk+lUDhgQN6x7gqTCYnBnpPJ/3Qy+zK7K53HHEJUtQIIYS4bDnr1nF82DASp01DlZbqHeeqGND5OkpwZs3BFErMFr3jiIuQokYIIcRlc+3QAYOrK0VHjpLxv6V6x7kqOjXywc/NRHZhMZsPn9Q7jrgIKWqEEEJcNqO3NwGTJwGQ+s47mLOydE5U84wGjeuCv6NFs2f5ZuODescRFyFFjRBCiCvifdttODVvhiUri/T/LNI7zlUREeBLkqPGXyRiMZv1jiMuQIoaIYQQV0QzGgmYMgWAM59+SklKir6BroJRfR7ByaJIdjSwbudyveOIC5CiRgghxBVzv/ZaXNq3RxUUkP7BB3rHqXF+3g2ILHEBYPX+z3VOIy5EihohhBBXTNM0AqZOBSDjy68ojo/XOVHN6+AZDcD+0qM6JxEXIkWNEEKIGuHW9RrcevaE0lJS335H7zg1bnj3R9CU4phJsf/YDr3jiEpIUSOEEKLGBDz+OADZP/5I4aHDOqepWS3C29G82AjAD1vt7xSbPZCiRgghRI1xiWqDx3XXgVKkzpund5wa18qlPQCrCu2vM7Q9kKJGCCFEjQqYNAmMRnLXriV/1y6949SogT2fBiDZEMfpnAyd04jzSVEjhBCiRjk1icD7plE4tWgBSukdp0b1bdwaQ2kAmmZm8a6VescR53HQO4AQQgj7E/TMM2guLmgG+/rb2WAw0NKjKzEFP7L72OfQd7TekcQ/2NfRJoQQok4wuLnZXUFzzpiAcADiOUp+YZ7OacQ/2ecRJ4QQok6wFBSQ/tFH5Pz2m95Rasz1Pe/Fy2zB22Jm8+av9I4j/kGKGiGEEFdNxmefkfLvOST/3/+hSkr0jlMjTCYnZuREsuJUEm4xG/WOI/5BihohhBBXjfftY3Bu1w7/CQ+DHZ2O8o8ciQY0Sl2Hslj0jiPOsp8jTAghRJ1jdHcj4suleI8aiWY06h2nxrTsMZxC5UiASmH/vg16xxFnSVEjhBCi1ig7ucTb1d2LuQEt6RMeyhe73tI7jjhLihohhBBXnSotJePLL4kdOQpzTo7ecWqEyb8tuQYDx4tP6B1FnCVFjRBCiKtP0zjz8ScUHTpExmef652mRtzUayqhJ29gS+wMTmcV6h1HIEWNEEKIWqAZjfhPeAiAM0uWYMnP1znRlWsS2gIX/xtQOLAqJlnvOAIpaoQQQtQSz6FDcQwLw5yZScaXX+odp0YMat0AgJX7k3ROIkCKGiGEELVEc3DA74H7ATizaDGW4mKdE125fk2d6R38GolqPKfTTukdp96TokYIIUSt8Ro5EoegIEpTUsj6Zrneca5Yy+AGpLhmk+JoYPmG9/SOU+9JUSOEEKLWGEwm/MbfB0D6woU2f5dhg9FIFKEA7Er9Xec0QooaIYQQtcr71lsx+vpSkpBA1k8/6R3nivVqOhKA/Q5ZMsClzqSoEUIIUasMLi74jhsHQPqChTY/zMC5AS5zjAa+X79A7zj1mhQ1Qgghap3PHWMweHpSfPw4OStX6R3niphMTrQp9QFgZ7z9jEZui6SoEUIIUeuM7u743nUnAGkLPrT54ROae0YBcNKcoHOS+k2KGiGEELrwuftuNFdXig7EkLd+vd5xrsg1La4H4IijIr+kSOc09ZcUNUIIIXTh4OOD3/3j8X/kEZzbtdM7zhXp0X4YmF0oNVhYdWSP3nHqLQe9AwghhKi/Ah55RO8INcLBaMTT0Ixs9rE6dhsjWnfVO1K9JC01QgghRA1o7lXWr+Zo2hadk9RfUtQIIYTQXe7GjZy48y4K9uzRO8plG+RZdgVUcckefYPUY1LUCCGE0F32Tyso2LmT9P/8R+8ol21ghxsxKEWyo4EDx3foHadekj41QgghdOf3wP0YPTzwve9evaNctiC/EJ5LNdK7KJ6kg3uhSWe9I9U70lIjhBBCd04REQQ9+wyOQUF6R7ki4S4dCS41U3xym95R6iUpaoQQQtQ5qrRU7wiXxaFR2VVPvum7dU5SP0lRI4QQos4oPHSYuAcfJPn11/WOcln8W/fiP14evOebTlbuGb3j1DtS1AghhKgzzBkZ5K3fQNa332HOzdU7TrU1imjNx55erHVzYe32r/SOU+9IUSOEEKLOcO16DaZmTbHk55P17Xd6x6k2g9FI55JI2qdEkFgQqnecesemiprY2Fj69+9P69atadu2LXl5eXpHEkIIUYM0TcNnzBgAMj7/3CYHumwZ+S/+SH+IP9Ntu9OzLbKpombcuHG8/PLLHDhwgN9//x0nJye9IwkhhKhhXiNGYnBzo/j4cfK32N7deaPDy27Ct/Nkpk0WZbbMZoqa/fv34+joSO/evQHw9fXFwUFusyOEEPbG6O6G14gRQFlrja2JCvEi3PkwzRyWsO/oHr3j1Cs1VtSsX7+e4cOHExwcjKZpfPvttxWWmT9/PhERETg7OxMdHc2GDRuqvP4jR47g7u7OjTfeSKdOnXjdRnvGCyGEuDSfO8pOQeWsXkNJYqLOaarH2dGIf/AS9jf8k3V7PtY7Tr1SY0VNXl4e7du359133610/tKlS5kyZQrTp09n9+7d9O7dm6FDhxIXF2ddJjo6mqioqAqPxMRESkpK2LBhA++99x6bN29m1apVrFq16oJ5ioqKyM7OLvcQQghhG5yaNcO1a1ewWMhY+qXecaqtkaGsP83RzL06J6lfaqyoGTp0KK+++io33XRTpfPnzp3L+PHjuf/++4mMjGTevHmEhYXx/vvvW5fZuXMnf/31V4VHcHAwoaGhdOnShbCwMJycnBg2bBh7LjLw2axZs/Dy8rI+wsLCampXhRBC1AKfO+8AIPOrr7AUF+ucpnpa+UUDcIJUnZPUL7XSp6a4uJidO3cyePDgctMHDx7Mpk2bqrSOLl26kJycTEZGBhaLhfXr1xMZGXnB5Z999lmysrKsj/j4+CvaByGEELXL49prcWjQAPOZM+T8+qvecaqlT/tbADjpqEhOT9A5Tf1RK0VNWloaZrOZoPPG9AgKCuL06dNVWoeDgwOvv/46ffr0oV27djRv3pwbbrjhgss7OTnh6elZ7iGEEMJ2aA4O+Iy+DYCMz2yrw3Bkk2iCSixYNI21O23v9JmtqtWrnzRNK/dcKVVh2sUMHTqUffv28ddffzF37tyajieEEKKO8b71VnB0pGDPHgr279c7TrU0tpT9Mf1Xwh86J6k/auWaaH9/f4xGY4VWmZSUlAqtN0IIIcQ5Dv7+eF1/PZaiQgw2dm+yCNeWbDXv5ETxCb2j1Bu10lJjMpmIjo6ucLXSqlWr6NGjR21EEEIIYaMaznqd0DffxKlZM72jVEuHVjcDsNfRQKnZrHOa+qHGiprc3Fz27NljvSIpNjaWPXv2WC/Znjp1Kh999BGLFi0iJiaGxx9/nLi4OCZMmFBTEYQQQtih6nRTqEv6tx+MsjiijIVsOHlA7zj1Qo2dftqxYwf9+/e3Pp86dSoAY8eOZcmSJYwePZr09HRefvllkpKSiIqKYsWKFYSHh9dUBCGEEHas6Phxsn/8Ef+JE9EMdf+G+K6OTrgRQT6H+fXoFvo3aat3JLunqXoyMEV2djZeXl5kZWXJlVBCCGFjVHExR3r3wZyVRdiHH+Det6/ekark3i+msKN4NdGljVgy/ie949ik6vx+1/1SVwghRL2nmUx43XIz7v37Y/T10ztOlfXxKMuarE7qnKR+kBEhhRBC2ITAJ5+0uf41/TvexrKfP6NDYRFn0hLx9Q/WO5Jdk5YaIYQQNsHWChqAxiEt+SBB47W0M8TvlfvVXG1S1AghhLApxacSSJk3D0tBgd5RquS0V3sA8o9VbVggcfmkqBFCCGEzlFLEjb+P9A8+JPsn2+h4q0KvwQwUntmudxS7J0WNEEIIm6FpGj633gpA5jfLdU5TRc1a0zM8lKcDsygsytc7jV2TokYIIYRN8Rw+HDSNgl27KElM1DvOJXWKGoCmwIBiy/71esexa1LUCCGEsCmOQUG4du4MQPbPP+uc5tIcHBzpWPokqYdeJ66old5x7JoUNUIIIWyO5/XDAMj+aYXOSaqmZZNeWHBg58kMvaPYNSlqhBBC2ByPIUPAaKTwwAGKYmP1jnNJ0eE+AOyKk6LmapKiRgghhM1x8PHBrUcPwDZaa9qFeNA75FVcfSZxKHa33nHslhQ1QgghbJL1FNSKFdT1YQw9XZ3IdM4l0aTx+95lesexW1LUCCGEsEkeAweimUwUHz9O0cGDese5pHD8ATiUukPnJPZLihohhBA2yejubh2t2xZuxNfEszUAiZZUnZPYLylqhBBC2CzP668HIMsGTkG1DOkCQJKxSOck9kuKGiGEEDbLvV9fDK6ulCYmUfjnn3rHuajOkYMASHcwEH/6uM5p7JMUNUIIIWyWwdmZoBdeIPyz/+Lcrp3ecS4qyC+EgFILALsOrtI5jX1y0DuAEEIIcSW8R43UO0KVNTC7kOpQxF/pxxmhdxg7JC01QgghRC0p8hoIwM4SnYPYKSlqhBBC2Lyi48dJeuklUv79b72jXFRzn6YAnC44qXMS+yRFjRBCCJtXmppG5hf/I+PLr7AUF+sd54I6Nigb0DJXndI5iX2SokYIIYTNc+0cjc9ddxEydy6a0ah3nAvqGdIMAGXMIiktTuc09keKGiGEEDZPMxpp8Px03Hv1rNNFTahfQ/zPXgG144BcAVXT5OonIYQQohbdlOtB4+KTuLrJnYVrmrTUCCGEsBuFhw6R/K/Z5Kxbp3eUC7rGMYrhefm4pZ7QO4rdkaJGCCGE3che8TNnFi8m86s6PBJ2QCQALllHdQ5if6SoEUIIYTc8hw0DIG/9eszZ2TqnqZxLaCu2ODux0yFe7yh2R4oaIYQQdsO5ZQucmjdDlZSQs+o3veNUyrtxcx5oGMSbASZSziTqHceuSFEjhBDCrpwbuTt7xQqdk1SuUcMWNC80EJVr4tApuQlfTZKiRgghhF3xHDoUgLwtWyhNT9c5TeWcDR+yOf5lkkvC9I5iV6SoEUIIYVdM4eE4t20LZjPZv/6qd5xKNQ90B+BISo7OSeyLFDVCCCHszrkOw3X1FFSzIA/AwsnEY3pHsStS1AghhLA7noMHAVCwazfmrCyd01TkVfwHDVs8wzH1jN5R7IoUNUIIIeyOY0gIpiZNwGIhb/MWveNU0LZJR3KNBpIdID3ztN5x7IYUNUIIIeySe+9eAORt/EPnJBVFhLbG22xBaRo7YtboHcduyNhP5zGbzZSUlOgdQwhRhzk6OmKsw4MmijJuvXpz5uNPyN3wB0opNE3TO1I5DUsdyTSaOXRqG0O4Q+84dkGKmrOUUpw+fZrMzEy9owghbIC3tzcNGjSocz+U4m+uXTqjOTlRevo0xUeP4tS8ud6RygnSfIghjfjsI3pHsRtS1Jx1rqAJDAzE1dVVvqiEEJVSSpGfn09KSgoADRs21DmRuBCDszPBs17H1KQJpmbN9I5TQbBrOJSmkVyaoncUuyFFDWWnnM4VNH5+fnrHEULUcS4uLgCkpKQQGBgop6LqsHOXdtdFEQHtIWknycYCvaPYDekoDNY+NK6urjonEULYinPfF9IHT1yuTq0GApDooJGZn6tzGvsgRc0/yCknIURVyfeF7chZs4aEp58mb+s2vaOU0yy0DZhdQYNN8Qf1jmMXpKgRQghh13LXriP7+x/IWV23Ru02GAy4EAzAjkQpamqC9KkRQghh1zyvH4bR2xuPQQP1jlJBgHMj4kqOEpuyB7hL7zg2T1pqbNy4ceMYOXJkrW5zyZIleHt71+o2hRDicrl160bgE1NxaddO7ygVdLNkA6Cy1uucxD5IUSOEEELopGlAFACnDfk6J7EPUtRcgFKK/OJSXR5KqcvK3K9fPyZNmsTTTz+Nr68vDRo0YObMmeWW0TSN999/n6FDh+Li4kJERARfffWVdf66devQNK3cTQj37NmDpmmcOHGCdevWce+995KVlYWmaWiaVmEbQghR11gKCshZt46ML7/UO0o5PSKHsigpmU8SkigqlMLmSkmfmgsoKDHT+sVfddn2gZeH4Gq6vI/m448/ZurUqWzdupXNmzczbtw4evbsyaBBg6zLvPDCC7zxxhu89dZbfPrpp4wZM4aoqCgiIyMvuf4ePXowb948XnzxRQ4dOgSAu7v7ZWUVQojaUnT0GKcmPIzBzQ3vUaPQHB31jgRAeGgrfAuNeFJE7LF9RLTpqnckmyYtNXamXbt2zJgxg+bNm3PPPffQuXNnVq9eXW6ZW2+9lfvvv58WLVrwyiuv0LlzZ955550qrd9kMuHl5YWmaTRo0IAGDRpIUSOEqPOc27TG6OODJS+Pgj179I5jpRkMJDqGA3DmxD6d09g+aam5ABdHIwdeHqLbti9Xu/M6wjVs2NB6O/dzunfvXuH5njr0P7kQQtQ0zWDArWdPsn/8kdw/NuLapYvekax2ejVghTkJQ/Iaorlf7zg2TYqaC9A07bJPAenJ8bwmVU3TsFgsl3zduRuJGQxljXf/7Ncjd0wVQtgDt15lRU3ehg3w+BS941gd93LnS4snnYuO6R3F5snpp3poy5YtFZ63atUKgICAAACSkpKs889vxTGZTJjN5qsbUgghaph7z54AFB44QGl6us5p/hYZOoAOmT44FHTQO4rNk6KmHvrqq69YtGgRhw8fZsaMGWzbto2JEycC0KxZM8LCwpg5cyaHDx/mp59+Ys6cOeVe37hxY3Jzc1m9ejVpaWnk50uPfSFE3ecQEIBT67ILIvI2btQ5zd96dLydDUnTWJs8guLSS7esiwuToqYeeumll/jf//5Hu3bt+Pjjj/nss89o3bo1UHb66osvvuDgwYO0b9+ef/3rX7z66qvlXt+jRw8mTJjA6NGjCQgIYPbs2XrshhBCVJt7z14A5P7xh85J/tbQyxl3JwdKLYqT6Xl6x7Fpmrrcm6LYmOzsbLy8vMjKysLT07PcvMLCQmJjY4mIiMDZ2VmnhLVD0zSWL19e63chFsLe1KfvDXuSt20bcfeMxejrS/M/NqAZ6sbf9re98wPZZ7Zwf/+e3NJnmN5x6pSL/X6fr258mkIIIUQtcO3QAYOrK+YzZyiMidE7jpWv02wSGn/PrmML9Y5i06SoEUIIUW9oJhOuZ29rkbehDp2Ccg4FILk46RJLiouRoqaeUUrJqSchRL3m3qvsKqi8OtSvJsKvLQDJBulTcyWkqBFCCFGvuPUq6yycv2cP5txcndOUadesDwAJjoqCIilsLpdNFTVvvvkmbdq0oXXr1kyaNOmyB34UQghRf5nCwvC6+SaCnnoS6sjvSFTTrrhYLJRqGrsO/q53HJtlM0VNamoq7777Ljt37mTfvn3s3Lmzwk3khBBCiKoIfu01fMeOxejhoXcUABwcHAkuLRsiZ/+JunMPHVtjM0UNQGlpKYWFhZSUlFBSUkJgYKDekYQQQogaEaTKLleOO3NQ5yS2q8aKmvXr1zN8+HCCg4PRNI1vv/22wjLz58+33tMhOjqaDRs2VHn9AQEBPPnkkzRq1Ijg4GAGDhxI06ZNayq+EEKIeqbk9Gkyly2jJCFB7ygAeHp0AGCPZtI3iA2rsaImLy+P9u3b8+6771Y6f+nSpUyZMoXp06eze/duevfuzdChQ4mLi7MuEx0dTVRUVIVHYmIiGRkZ/Pjjj5w4cYKEhAQ2bdrE+vXrL5inqKiI7Ozscg8hhBDinMRnnyXp+RfIWb1a7ygAhIUPAiDRkqFzEttVY0XN0KFDefXVV7npppsqnT937lzGjx/P/fffT2RkJPPmzSMsLIz333/fuszOnTv566+/KjyCg4P57bffaNasGb6+vri4uHD99ddftE/NrFmz8PLysj7CwsJqalfFJYwbN04uG68h/fr1Y8qUKXrHAGDmzJl06NChysufOHECTdMqDIgqRF3h3qcvzu3bYfTx0TsKAF1Dy4arKdaSKSot0TmNbaqVPjXFxcXs3LmTwYMHl5s+ePBgNm3aVKV1hIWFsWnTJgoLCzGbzaxbt46WLVtecPlnn32WrKws6yM+Pv6K9qEu0jTtoo9x48bpHbFK7OnHTwo6IWyH77ixRCxditfw4XpHAaBjcASaxQHNUMrOo9v1jmOTHGpjI2lpaZjNZoKCgspNDwoK4vTp01VaR7du3Rg2bBgdO3bEYDAwYMAAbrzxxgsu7+TkhJOT0xXlruuSkv6+8+TSpUt58cUXOXTokHWai4uLHrFEFZSUlODo6Kh3DCHqNU3T9I5QjsnBgYjSUo6b4MDhX+jRqofekWxOrV79dP4BpJSq1kH12muvERMTw/79+3n77bdr54Aszrvwo6SwGssWVG3ZamjQoIH14eXlhaZp5aatX7+e6OhonJ2dadKkCS+99BKlpaXW12uaxocffsgNN9yAq6srkZGRbN68maNHj9KvXz/c3Nzo3r07x44ds77m3CmIDz/8kLCwMFxdXbn11lvJzMy8YM5ffvmFXr164e3tjZ+fHzfccEO5dUZERADQsWNHNE2jX79+1nmLFy8mMjISZ2dnWrVqxfz58y/6nixbtoy2bdvi4uKCn58fAwcOJC+v7H0914ry0ksvERgYiKenJw899BDFxcXW1yulmD17Nk2aNMHFxYX27duzbNmyctvYv38/119/PZ6ennh4eNC7d2+OHTvGzJkz+fjjj/nuu++srWXr1q2ztkR9+eWX9OvXD2dnZ/773/+Snp7OmDFjCA0NxdXVlbZt2/LFF19cdP/Od+7zWLRoEY0aNcLd3Z2HH34Ys9nM7NmzadCgAYGBgbz22mvlXhcXF8eIESNwd3fH09OT2267jeTk5HLLvPHGGwQFBeHh4cH48eMpLDzveL+Mz0eIusicm0dJSoreMQAIVB44WSycyTh26YVFBbXSUuPv74/RaKzQKpOSklKh9abOeT34wvOaD4Y7v/r7+f81g5L8ypcN7wX3/vT383ltIT+94nIzsy4v53l+/fVX7rrrLt5++23rj+6DDz4IwIwZM6zLvfLKK8ydO5e5c+cybdo07rjjDpo0acKzzz5Lo0aNuO+++5g4cSI///yz9TVHjx7lyy+/5IcffiA7O5vx48fz6KOP8tlnn1WaJS8vj6lTp9K2bVvy8vJ48cUXGTVqFHv27MFgMLBt2zauueYafvvtN9q0aYPJVNbzf+HChcyYMYN3332Xjh07snv3bh544AHc3NwYO3Zshe0kJSUxZswYZs+ezahRo8jJyWHDhg3lbtK4evVqnJ2dWbt2LSdOnODee+/F39/f+qP//PPP88033/D+++/TvHlz1q9fz1133UVAQAB9+/YlISGBPn360K9fP9asWYOnpycbN26ktLSUJ598kpiYGLKzs1m8eDEAvr6+JCYmAjBt2jTmzJnD4sWLcXJyorCwkOjoaKZNm4anpyc//fQTd999N02aNKFr165V/qyPHTvGzz//zC+//MKxY8e45ZZbiI2NpUWLFvz+++9s2rSJ++67jwEDBtCtWzfrUBlubm78/vvvlJaW8sgjjzB69GjWrVsHwJdffsmMGTN477336N27N59++ilvv/02TZo0sW63up+PEHXRmU8+IXnWG3iNHEnwrNf1jsMYp558ELuA3Z6RekexTeoqANTy5cvLTbvmmmvUww8/XG5aZGSkeuaZZ65GhAqysrIUoLKysirMKygoUAcOHFAFBQUVXzjD88KP/95SftlXG1x42UXDyi/7r4jKl7tMixcvVl5eXtbnvXv3Vq+//nq5ZT799FPVsGFD63NAPf/889bnmzdvVoD6z3/+Y532xRdfKGdn57/fjhkzlNFoVPHx8dZpP//8szIYDCopKUkppdTYsWPViBEjLpg1JSVFAWrfvn1KKaViY2MVoHbv3l1uubCwMPX555+Xm/bKK6+o7t27V7renTt3KkCdOHGi0vljx45Vvr6+Ki8vzzrt/fffV+7u7spsNqvc3Fzl7OysNm3aVO5148ePV2PGjFFKKfXss8+qiIgIVVxcfMFtnL/v5/Zv3rx5lb7mn4YNG6aeeOIJ6/O+ffuqyZMnX3D5GTNmKFdXV5WdnW2dNmTIENW4cWNlNput01q2bKlmzZqllFJq5cqVymg0qri4OOv8/fv3K0Bt27ZNKaVU9+7d1YQJE8ptq2vXrqp9+/bW55f6fC70udqLi35vCJuRvWaNOtCylTp6/fV6R1FKKbXr10+VmuGpjrzcUe8odcbFfr/PV2MtNbm5uRw9etT6PDY2lj179uDr60ujRo2YOnUqd999N507d6Z79+4sWLCAuLg4JkyYUFMRro7nEi88TzOWf/7U0cqXA9DOO9M3Zd/lZ6qCnTt3sn379nKnHcxmM4WFheTn5+Pq6gpAu3btrPPPtZq1bdu23LTCwkKys7Px9Cy7MVSjRo0IDQ21LtO9e3csFguHDh2iQYMGFbIcO3aMF154gS1btpCWlobFYgHKToFERUVVmj81NZX4+HjGjx/PAw88YJ1eWlqKl5dXpa9p3749AwYMoG3btgwZMoTBgwdzyy234POPKxvat29v3fdz2XNzc4mPjyclJYXCwkIGDRpUbr3FxcV07NgRgD179tC7d+/L6g/TuXPncs/NZjNvvPEGS5cuJSEhgaKiIoqKinBzc6vWehs3bozHP+6KGhQUhNFoxGAwlJuWcrZ5PSYmhrCwsHJXBLZu3Rpvb29iYmLo0qULMTExFf7f7N69O2vXrgUu7/MRoi5yOft9V3zsOObcPIzu1fv/r6b5hLUCINBctf6morwaK2p27NhB//79rc+nTp0KwNixY1myZAmjR48mPT2dl19+maSkJKKiolixYgXh4eE1FeHqMFXjAL9ay14Gi8XCSy+9VOkl9s7OztZ///PH+VwfpcqmnStEKnNumQv1cRo+fDhhYWEsXLiQ4OBgLBYLUVFR5fqyVJYfyk5xnH8qxmg0VvYSjEYjq1atYtOmTaxcuZJ33nmH6dOns3XrVmu/nYvtw7lt/vTTT4SEhJSbf67T+ZV0vj6/WJkzZw5vvvkm8+bNo23btri5uTFlypSLvi+VOb/A0jSt0mnn9k9doC/bhaZX5nI+HyHqIgd/fxyCG1KamETh/v24db1G1zw+wY15PNCfBAcH5qWeJDigjv9G1jE1VtT069fvkgNMPvLIIzzyyCM1tUlxEZ06deLQoUM0a9asxtcdFxdHYmIiwcFl/Y02b96MwWCgRYsWFZZNT08nJiaGDz/8kN69ewPwxx9/lFvmXB8as9lsnRYUFERISAjHjx/nzjvvrHI2TdPo2bMnPXv25MUXXyQ8PJzly5dbi+w///yTgoICa3GyZcsW3N3dCQ0NxcfHBycnJ+Li4ujbt2+l62/Xrh0ff/zxBa9eMplM5fbjYjZs2MCIESO46667gLJC4ciRI0RGXt1z6a1btyYuLo74+Hhra82BAwfIysqybjsyMpItW7Zwzz33WF/3z/tCXe7nI0Rd5BLVlpzEJAr/2qd7UePl5c92Z2eyjAYOxm6XoqaaaqWjsKh9L774IjfccANhYWHceuutGAwG9u7dy759+3j11VevaN3Ozs6MHTuWf//732RnZzNp0iRuu+22Sk89+fj44Ofnx4IFC2jYsCFxcXE888wz5ZYJDAzExcWFX375hdDQUJydnfHy8mLmzJlMmjQJT09Phg4dSlFRETt27CAjI8NapPzT1q1bWb16NYMHDyYwMJCtW7eSmpparkgoLi5m/PjxPP/885w8eZIZM2YwceJEDAYDHh4ePPnkkzz++ONYLBZ69epFdnY2mzZtwt3dnbFjxzJx4kTeeecdbr/9dp599lm8vLzYsmUL11xzDS1btqRx48b8+uuvHDp0CD8/v4ueimnWrBlff/01mzZtwsfHh7lz53L69OmrXtQMHDiQdu3aceeddzJv3jxrR+G+fftaT5FNnjyZsWPH0rlzZ3r16sVnn33G/v37y3UUru7nI0Rd5dKuLTkrV1Kw9+p2C6iqtnk9Sck2UBjVXO8oNsemBrQUVTdkyBB+/PFHVq1aRZcuXejWrRtz586tkdN9zZo146abbmLYsGEMHjyYqKioC17KazAY+N///sfOnTuJiori8ccf5//+7//KLePg4MDbb7/Nhx9+SHBwMCNGjADg/vvv56OPPmLJkiW0bduWvn37smTJkgueSvL09GT9+vUMGzaMFi1a8PzzzzNnzhyGDh1qXWbAgAE0b96cPn36cNtttzF8+HBmzpxpnf/KK6/w4osvMmvWLCIjIxkyZAg//PCDdZt+fn6sWbOG3Nxc+vbtS3R0NAsXLrS22jzwwAO0bNmSzp07ExAQwMaNFx5t94UXXqBTp04MGTKEfv360aBBg1q5cd+5sdl8fHzo06cPAwcOpEmTJixdutS6zOjRo3nxxReZNm0a0dHRnDx5kocffrjceqr7+QhRVzlHlfWrKdxXN4oai/eD7MwZTmqR9E+rLk1d6pyRncjOzsbLy4usrCxrh9dzCgsLiY2NtQ62KS5s5syZfPvttzZ5999x48aRmZlZ6WCrQlSXfG/YD3NuLoe7XANK0XzjHzj4+ema5+UfDrBoYywP9mnCc8Pk0u6L/X6fT1pqhBBC1GtGd3dMZ0+tFtSB1poGTvFc4/UVWQkL9Y5ic6SoEUIIUe+du7S7sA70q9EKNhITvJMj6je9o9gcKWpEtcycOdMmTz0BLFmyRE49CSEq5dy27J5ZBX/pX9SEB5adckpzqNqVlOJvUtQIIYSo9/7ZUqN3V9PIxmVXIWYYDWTmpOmaxdZIUSOEEKLec2rVChwdUSUlmNP0LSRCApvgevYGlweOb9c1i62RokYIIUS9ZzCZaPrzClps34ZDQIC+WYxG/EvLfp5jk/Q/HWZLpKgRQgghAFNoKJqhbvws+ljKhmZJzLjImIKigrrx6QkhhBDCytdQNkhtWsFFBlUWFUhRI4QQQgDmnBwSnnqa4yNGokpLdc1i8CrrLPyX0VfXHLZGihphMxo3bsy8efP0jqGrfv36MWXKFL1jAGWX93fo0KHKy584cQJN02z2lgDC/hlcXclds4aiQ4coOnZc1ywBgR0BSDXn6JrD1khRYwfi4+MZP348wcHBmEwmwsPDmTx5Munp6dVe17lxgf5pyZIlaJqGpmkYjUZ8fHzo2rUrL7/8MllZWTW0F+W35+3tXePr1cO4ceNqZTwnIcSV04xGgp5/nrCFCzCFheqapblfIwAKlVzSXR1S1Ni448eP07lzZw4fPswXX3zB0aNH+eCDD1i9ejXdu3fnzJkzNbIdT09PkpKSOHXqFJs2beLBBx/kk08+oUOHDiQmyjnfK1VSUqJ3BCEE4D1qJO69e2NwddU1R9uAs0WVMZP8wjxds9gSKWouIb8kv9qPUsvf52JLLaXkl+RTWFpYpfVW16OPPorJZGLlypX07duXRo0aMXToUH777TcSEhKYPn26ddnGjRvzyiuvcMcdd+Du7k5wcDDvvPNOufkAo0aNQtM063Moa8Fp0KABDRs2JDIykvHjx7Np0yZyc3N5+umnrcsppZg9ezZNmjTBxcWF9u3bs2zZMuv8devWoWkaP/30E+3bt8fZ2ZmuXbuy7+x4K+vWrePee+8lKyvL2jr0z1G08/Pzue+++/Dw8KBRo0YsWLDgou/PsmXLaNu2LS4uLvj5+TFw4EDy8sq+IM61orz00ksEBgbi6enJQw89RHFxcZX3B2D//v1cf/31eHp64uHhQe/evTl27BgzZ87k448/5rvvvrPuy7p166ynYb788kv69euHs7Mz//3vf0lPT2fMmDGEhobi6upK27Zt+eKLLy5xBJR37pTQokWLaNSoEe7u7jz88MOYzWZmz55NgwYNCAwM5LXXXiv3uri4OEaMGIG7uzuenp7cdtttJCcnl1vmjTfeICgoCA8PD8aPH09hYfljGmDx4sVERkbi7OxMq1atLjh6uxDi4pr7h+CoFEqDg8d36B3Hdqh6IisrSwEqKyurwryCggJ14MABVVBQUGFe1JKoaj9+if3F+vpfYn9RUUui1Lifx5Vbb+8velf62upIT09Xmqap119/vdL5DzzwgPLx8VEWi0UppVR4eLjy8PBQs2bNUocOHVJvv/22MhqNauXKlUoppVJSUhSgFi9erJKSklRKSopSSqnFixcrLy+vSrcxefJk5eHhoUpLS5VSSj333HOqVatW6pdfflHHjh1TixcvVk5OTmrdunVKKaXWrl2rABUZGalWrlyp9u7dq2644QbVuHFjVVxcrIqKitS8efOUp6enSkpKUklJSSonJ8ea39fXV7333nvqyJEjatasWcpgMKiYmJhKsyUmJioHBwc1d+5cFRsbq/bu3avee+896/rGjh2r3N3d1ejRo9Vff/2lfvzxRxUQEKCee+456zoutT+nTp1Svr6+6qabblLbt29Xhw4dUosWLVIHDx5UOTk56rbbblPXXXeddV+KiopUbGysAlTjxo3V119/rY4fP64SEhLUqVOn1P/93/+p3bt3q2PHjlk/ny1btljz9O3bV02ePPmCx8SMGTOUu7u7uuWWW9T+/fvV999/r0wmkxoyZIh67LHH1MGDB9WiRYsUoDZv3qyUUspisaiOHTuqXr16qR07dqgtW7aoTp06qb59+1rXu3TpUmUymdTChQvVwYMH1fTp05WHh4dq3769dZkFCxaohg0bWvfp66+/Vr6+vmrJkiVKKWXd7927d18wvy252PeGsF0Ws1llrVypkv89R5l1/mxvWNBG3bCgpfppzQJdc+jtYr/f53PQsZ4SV+jIkSMopYiMrHxo+sjISDIyMkhNTSUwMBCAnj178swzzwDQokULNm7cyJtvvsmgQYMIOHvDKW9vbxo0aFClDK1atSInJ4f09HTc3NyYO3cua9asoXv37gA0adKEP/74gw8//JC+fftaXzdjxgwGDRoEwMcff0xoaCjLly/ntttuw8vLy9oydL5hw4bxyCOPADBt2jTefPNN1q1bR6tWrSosm5SURGlpKTfddBPh4eEAtD17K/RzTCYTixYtwtXVlTZt2vDyyy/z1FNP8corr1BQUHDJ/Xnvvffw8vLif//7H46Ojtb39RwXFxeKiooq3ZcpU6Zw0003lZv25JNPWv/92GOP8csvv/DVV1/RtWvXC30EFVgsFhYtWoSHhwetW7emf//+HDp0iBUrVmAwGGjZsiX/+te/WLduHd26deO3335j7969xMbGEhYWBsCnn35KmzZt2L59O126dGHevHncd9993H///QC8+uqr/Pbbb+Vaa1555RXmzJlj3aeIiAgOHDjAhx9+yNixY6ucXwhdaRqnX3oZc1oa7tf2x7VjR92ivJbmR7uiPWwPKNItg62RouYStt6xtdqvMRlN1n8PaDSArXdsxaCVP9P3y82/XHG2S1Fnxy/RNM067dyP8z+fX8kVRf/cxoEDBygsLLQWK+cUFxfT8bwvhn/m8PX1pWXLlsTExFxye+3atbP++1zhk5KSUumy7du3Z8CAAbRt25YhQ4YwePBgbrnlFnx8fMot4/qPc+fdu3cnNzeX+Ph4UlJSLrk/e/bsoXfv3taCpjo6d+5c7rnZbOaNN95g6dKlJCQkUFRURFFREW5ubtVab+PGjfHw8LA+DwoKwmg0YvjHTcWCgoKs71tMTAxhYWHWggagdevWeHt7ExMTQ5cuXYiJiWHChAnlttO9e3fWrl0LQGpqqrXD+gMPPGBdprS0FC8vr2rlF0JPmqbhEhVF7rp1FO77S9eiptA1BIr2UHomTrcMtkaKmktwdbyyzmIOBgccDBXf5itdL0CzZs2sxURlV9gcPHgQHx8f/P39L7qefxY91RUTE4Onpyd+fn4cP152CeRPP/1ESEhIueWcnJwuua6q5Di/eNA0DcvZMVLOZzQaWbVqFZs2bWLlypW88847TJ8+na1btxIREXHJLOfWe7H9cXFxuWTmCzm/WJkzZw5vvvkm8+bNo23btri5uTFlypRyfXyqorL36GLvm1Kq0vf+QtMrc25dCxcurNCqZDQaq5xdiLrAuV1bcteto2CfvkMUmD3DIAOM2fG65rAl0lHYhvn5+TFo0CDmz59PQUFBuXmnT5/ms88+Y/To0eV+mLZs2VJuuS1btpQ7dePo6IjZXLXh7lNSUvj8888ZOXIkBoOB1q1b4+TkRFxcHM2aNSv3+GcrwPk5MjIyOHz4sDWHyWSqcoZL0TSNnj178tJLL7F7925MJhPLly+3zv/zzz/LvXdbtmzB3d2d0NDQKu1Pu3bt2LBhwwWvXqrOvmzYsIERI0Zw11130b59e5o0acKRI0euYO+rpnXr1sTFxREf//cX54EDB8jKyrKe2oyMjKz02DknKCiIkJAQjh8/XuG9ulQBKURdYx2xW+ei5qiHkdHBQSxwPqRrDlsiLTU27t1336VHjx4MGTKEV199lYiICPbv389TTz1FSEhIhatcNm7cyOzZsxk5ciSrVq3iq6++4qeffrLOb9y4MatXr6Znz544OTlZT9UopTh9+jRKKTIzM9m8eTOvv/46Xl5evPHGGwB4eHjw5JNP8vjjj2OxWOjVqxfZ2dls2rQJd3f3cv0qXn75Zfz8/AgKCmL69On4+/tbW5saN25Mbm4uq1evtp4ecr2Myyu3bt3K6tWrGTx4MIGBgWzdupXU1NRyfZCKi4sZP348zz//PCdPnmTGjBlMnDgRg8FQpf2ZOHEi77zzDrfffjvPPvssXl5ebNmyhWuuuYaWLVvSuHFjfv31Vw4dOoSfn99FT8U0a9aMr7/+mk2bNuHj48PcuXM5ffr0BftM1ZSBAwfSrl077rzzTubNm0dpaSmPPPIIffv2tZ4imzx5MmPHjqVz58706tWLzz77jP3799OkSRPrembOnMmkSZPw9PRk6NChFBUVsWPHDjIyMpg6depV3QchapJzVBQAxSdOYM7OxujpqU8On2AO5DkRXFK91tr6TFpqbFzz5s3ZsWMHTZs2ZfTo0TRt2pQHH3yQ/v37s3nzZnx9y99i+4knnmDnzp107NjR2rFzyJAh1vlz5sxh1apVhIWFlesHk52dTcOGDQkJCaF79+7Wzp+7d++mYcOG1uVeeeUVXnzxRWbNmkVkZCRDhgzhhx9+qPDX+htvvMHkyZOJjo4mKSmJ77//HpOprC9Sjx49mDBhAqNHjyYgIIDZs2df1nvj6enJ+vXrGTZsGC1atOD5559nzpw5DB061LrMgAEDaN68OX369OG2225j+PDh5S4hv9T++Pn5sWbNGnJzc+nbty/R0dEsXLjQerrngQceoGXLlnTu3JmAgAA2btx4wbwvvPACnTp1YsiQIfTr148GDRrUyo37zt1w0cfHhz59+jBw4ECaNGnC0qVLrcuMHj2aF198kWnTphEdHc3Jkyd5+OGHy63n/vvv56OPPmLJkiW0bduWvn37smTJEmmpETbHwccHx7OtsYV//aVbjrbNr6PFqZ6ohDFYLEq3HLZEU+d6etq57OxsvLy8yMrKwvO8qruwsJDY2FgiIiJwdnbWKeHV17hxY6ZMmaLrbfbXrVtH//79ycjI0P2uwePGjSMzM7PCHZSFqIr68r1RXyVMfYLsFSsImDIF/wkP6ZKhxGyh5fM/Y1Gw7bkBBHrWz+PsYr/f55OWGiGEEOI8zmf71RT8pV+/GkejgYZeZRcjxGcUXGJpAVLUCCGEEBW4tDvbWXivvp2FO3isobf/++w/+I2uOWyFdBSuR06cOKF3BPr160ddOeO5ZMkSvSMIIeoo58hIMBgoTUmhJDkZx6AgXXKUOK5lj3s2TU6vBMbrksGWSEuNEEIIcR6DqytOzZoB+l7a7eNQdrHHmeLKbzIqypOiRgghhKiE89lTUAX79LsCyt8lGIAMS45uGWyJnH4SQgghKuHWvTuWrCycW7a49MJXSbBPM0jeRIZBxn+qCilqhBBCiEp4XX89Xtdfr2uGJsFtIRnSHCxYzGYMMuzIRcnpJyGEEKKOat3kGgDyDQbik4/pnKbuk6JGCCGEuAClFMWnTlGSmKjL9r3cffEpLRsw9uDJnbpksCVS1AghhBAXkDpnDscGDiJdx1tAKBoA8JcMAXVJUtQIIYQQF+DUogU4OmLJy9Mtg4NjOADHc5J1y2ArpKgRNqNfv366jltVF9S196C6eepafiEuxWPIEFru3EHwa6/plsHfuaylJik3SbcMtkKKGjsRHx/P+PHjCQ4OxmQyER4ezuTJk0lPT6/Wei70ozNu3Dg0TUPTNBwdHQkKCmLQoEEsWrQIi8VSQ3tx6Ry2yJ72xZ7NmjWLLl264OHhQWBgICNHjuTQoUN6xxI6Mzg5YTCZdM3QxHR2IMu8PbrmsAVS1NiB48eP07lzZw4fPswXX3zB0aNH+eCDD1i9ejXdu3fnzJkzNbKd6667jqSkJE6cOMHPP/9M//79mTx5MjfccAOlpaU1so36qrhYTpbr7ffff+fRRx9ly5YtrFq1itLSUgYPHkyejqcdhACIPHsVt9miT2dlWyJFzSVY8vOr/VD/+IFXpaVl0wsLq7Tey/Hoo49iMplYuXIlffv2pVGjRgwdOpTffvuNhIQEpk+fDpS1GEycOJGJEyfi7e2Nn58fzz//vHUspnHjxvH777/z1ltvWVtl/jlelJOTEw0aNCAkJIROnTrx3HPP8d133/Hzzz9bx1FSSjF79myaNGmCi4sL7du3Z9myZeXyXkkOi8XC008/ja+vLw0aNGDmzJkXfW+WLVtG27ZtcXFxwc/Pj4EDB5b7kbpUlqrsj8Vi4V//+hfNmjXDycmJRo0a8drZpuoL7cu57U6dOhV/f38GDRoEwC+//EKvXr2sWW644QaOHaveZZz9+vXjscceY8qUKfj4+BAUFMSCBQvIy8vj3nvvxcPDg6ZNm/Lzzz9bX1NUVMSkSZMIDAzE2dmZXr16sX379nLrzcvL45577sHd3Z2GDRsyZ86ccvOr8l5drm3bttGvXz9cXFxo1aoV27dvZ8GCBdx44401sn4oe+/HjRtHmzZtaN++PYsXLyYuLo6dO+WKk/ouc9kyYm+6mbQPF+iy/YiGUYSXlBAuf/xcmqonsrKyFKCysrIqzCsoKFAHDhxQBQUFFeYdaNmq2o+sn3/+e7s//6wOtGylTtx1d7n1HurWvdLXVld6errSNE29/vrrlc5/4IEHlI+Pj7JYLKpv377K3d1dTZ48WR08eFD997//Va6urmrBggVKKaUyMzNV9+7d1QMPPKCSkpJUUlKSKi0tVUopNXbsWDVixIhKt9G+fXs1dOhQpZRSzz33nGrVqpX65Zdf1LFjx9TixYuVk5OTWrdunXX5y83Rt29f5enpqWbOnKkOHz6sPv74Y6Vpmlq5cmWluRITE5WDg4OaO3euio2NVXv37lXvvfeeysnJqXKWquzP008/rXx8fNSSJUvU0aNH1YYNG9TChQsvuS/u7u7qqaeeUgcPHlQxMTFKKaWWLVumvv76a3X48GG1e/duNXz4cNW2bVtlNputeSdPnnzhA+LsMh4eHuqVV15Rhw8fVq+88ooyGAxq6NChasGCBerw4cPq4YcfVn5+fiovL08ppdSkSZNUcHCwWrFihdq/f78aO3as8vHxUenp6db1Pvzwwyo0NFStXLlS7d27V91www3W9646n/2l8p9v8+bNytnZWc2aNUsdPnxYjRo1Sg0bNkw1a9ZM7dq1q8Lyr732mnJzc7voY/369Zfc7pEjRxSg9u3bV+n8i31vCPuS9p9F6kDLVip+8hRdtp+dma7UDE+lZniq3OwMXTLo6WK/3+eTokbZdlGzZcsWBajly5dXOn/u3LkKUMnJyapv374qMjJSWSwW6/xp06apyMhI6/ML/ehcrKgZPXq0ioyMVLm5ucrZ2Vlt2rSp3Pzx48erMWPGlNvG5eTo27ev6tWrV7lpXbp0UdOmTas0186dOxWgTpw4Uen8S2Wpyv5kZ2crJycnaxFzoW1Uti8dOnS44GvOSUlJKffDWtWi5p/vU2lpqXJzc1N33/33MZiUlKQAtXnzZpWbm6scHR3VZ599Zp1fXFysgoOD1ezZs5VSSuXk5CiTyaT+97//WZdJT09XLi4uavLkydX67Ktb1HTv3l3deeed1udLly5VBoNBjRo1qtLl09PT1ZEjRy76yM/Pv+g2LRaLGj58eIXj7Z+kqKk/steuVQdatlLHbhyhW4asGQ2VmuGpYg9s1y2DXqpT1MgwCZfQclf1m561f3Qq8xg4sGwdhvJn+pqt/u2Ks1WFOnsaRdM0ALp162b9N0D37t2ZM2cOZrMZ42XeflsphaZpHDhwgMLCQuuplHOKi4vp2LFjuWmXm6Ndu3blnjds2JCUlMpHr23fvj0DBgygbdu2DBkyhMGDB3PLLbfg4+NTpSx//fXXJfcnJiaGoqIiBgwYcMHMF9K5c+cK044dO8YLL7zAli1bSEtLs3bCjouLIyoqqsrr/uf7ZDQa8fPzo23bttZpQUFBAKSkpHDs2DFKSkro2bOndb6joyPXXHMNMTEx1lzFxcV0797duoyvry8tW7YEqNZnXx2nTp1i8+bN/N///Z91mslkQinFSy+9VOlrfH198fX1vextAkycOJG9e/fyxx9/XNF6hH1wiogAoPjECZTFgmao/Z4bqcZAPM2xZCYeg8iK3x2ijBQ1l2Bwdb2i12sODmgOFd/mK13vOc2aNbMWFCNHjqww/+DBg/j4+ODv718j26tMTEwMERER1h/gn376iZCQkHLLODk51ci2HB0dyz3XNO2CV18ZjUZWrVrFpk2bWLlyJe+88w7Tp09n69atRJz9kqqKi+2Pi4tLNffgb25ubhWmDR8+nLCwMBYuXEhwcDAWi4WoqKhqdySu7H3657RzRZzFYqlQ+J5zrlg99++LuVqf/bmi6p8F4KFDh7jmmmvKFWn/9Prrr/P6669fdL0///wzvXv3rnTeY489xvfff8/69esJDQ29zOTCnjiGhKA5OqKKiihJTMIUGnLpF9WwD/1c2OocwrDkVXRgdK1v31ZIR2Eb5+fnx6BBg5g/fz4FBQXl5p0+fZrPPvuM0aNHW3+ctmzZUm6ZLVu20Lx5c2vriMlkwmw2V3n7a9asYd++fdx88820bt0aJycn4uLiaNasWblHWFhYhe3WZI4L0TSNnj178tJLL7F7925MJhPLly+vUpaq7E/z5s1xcXFh9erVF8xQ1X1JT08nJiaG559/ngEDBhAZGUlGRsZl7HX1NGvWDJPJVK5VoqSkhB07dhAZGWldxtHRsdx7lZGRweHDhwGq9dlXR1ZWVrmWuzNnzjB79uyLFkoTJkxgz549F31U1kqmlGLixIl88803rFmzplqFr7BvmoMDjuGNACiOjdUlQ5GjG2eMRtIK5Aqoi5GWGjvw7rvv0qNHD4YMGcKrr75KREQE+/fv56mnniIkJMR6JQ6U3c9m6tSpPPTQQ+zatYt33nmn3FUsjRs3ZuvWrZw4cQJ3d3d8fX0xnG1qLSoq4vTp05jNZpKTk/nll1+YNWsWN9xwA/fccw9Go5Enn3ySxx9/HIvFQq9evcjOzmbTpk24u7szduzYK8pRXVu3bmX16tUMHjyYwMBAtm7dSmpqqvWH+lJZPDw8Lrk/zs7OTJs2jaeffhqTyUTPnj1JTU1l//79jB8/vlr74uPjg5+fHwsWLKBhw4bExcXxzDPPVHu/q8vNzY2HH36Yp556Cl9fXxo1asTs2bPJz8+37oO7uzvjx4/nqaeews/Pj6CgIKZPn249NqryXl2ODh06YDabmT17NrfeeiuTJ08mPDycmJgYTp48SXh4eIXXXO7pp0cffZTPP/+c7777Dg8PD06fPg2Al5fXFbXICfvgFBFB8dFjZUVN7161vv1WwRPYv30/JU26X3rh+uxqdu6pSy63o7CtOHHihBo3bpxq0KCBcnR0VGFhYeqxxx5TaWlp1mX69u2rHnnkETVhwgTl6empfHx81DPPPFOuk+yhQ4dUt27dlIuLiwJUbGysUqqsozCgAOXg4KACAgLUwIED1aJFi6xX5ihV1sHyrbfeUi1btlSOjo4qICBADRkyRP3+++9XnKOyTqYjRoxQY8eOrfQ9OXDggBoyZIgKCAhQTk5OqkWLFuqdd94pt8ylslRlf8xms3r11VdVeHi4cnR0VI0aNSp3NVpV90UppVatWqUiIyOVk5OTateunVq3bl25juBV7Sh8/jLh4eHqzTffLDftn+stKChQjz32mPL391dOTk6qZ8+eatu2beWWz8nJUXfddZdydXVVQUFBavbs2eW2VdXP/p/ZFi9erC71NfTyyy8rPz8/5ezsrMaOHatSU1NVp06dVLNmzS76uuo6d3yf/1i8eHGly9vD94aouuQ5c9WBlq1U4syZumx/xd5EFT7tRzXyvT902b6eqtNRWFPqEifL7UR2djZeXl5kZWXh6elZbl5hYSGxsbFERETg7OysU8Krr1+/fnTo0IF58+ZJjjqYpT6aOXMm69atY926dXpHqbb68r0hymQu/5akZ5/FtVs3wpcsrvXt7z2VyY3vbiTQw4lt0wfW+vb1dLHf7/PJ6SchhG5+/fVX3nrrLb1jCHFJTk3OXgF1/Lgu2w9yg97+74NjFjl5nfFw89YlR10nRY0QQjebN2/WO4IQVWI623G8NCUFc24eRveKVy9eTQGeHsT4naDIoBFzfAfXtK1frTVVJUVNPVJXmvjrSg6oW1mEEHWX0dMTo58f5vR0ik+cwCWqTa1u32A0EmCGUwY4nrRXipoLkEu6hRBCiCqw3oQvVp9TUL7mslsZJJw5osv2bYG01AghhBBV4HP33XiNHIFLx066bN/b4AGkk5afoMv2bYEUNUIIIUQVeA4ZrOv2/UwBQDpnSs/omqMuk9NP/1BPrm4XQtQA+b4QtS3Arezu3Bnk6pyk7pKihr/HycnPz9c5iRDCVpz7vjh/nC1hv5TZTN62bWR8+SXqAmPOXU1hfmUDyJ4xlNT6tm2FnH6ibOBDb29v62jPrq6uFQb3E0IIKGuhyc/PJyUlBW9v78se3V7YIKWIH38/qqQEtx49a31gyybN+kDcuyQ7OFBcWoqpksGS6zt5R85q0KABgLWwEUKIi/H29rZ+b4j6QXNwwLVbN9BAFRfX+vZbNWyOUgY0zcKB1FN0aNi41jPUdXWyqBk1ahTr1q1jwIABLFu2rNy8H3/8kSeeeAKLxcK0adO4//77a2SbmqbRsGFDAgMDKSmRpj0hxIU5OjpKC0091WjhAt22bXJwwGj2weKQzl/JsVLUVKJOFjWTJk3ivvvu4+OPPy43vbS0lKlTp7J27Vo8PT3p1KkTN91002WNyHshRqNRvqyEEELUSUHKgSTg9Ml10KG/3nHqnDrZUbh///54eHhUmL5t2zbatGlDSEgIHh4eDBs2jF9//VWHhEIIIeozi04XljSxFAGQmfGXLtuv66pd1Kxfv57hw4cTHByMpml8++23FZaZP3++deTa6OhoNmzYUBNZSUxMJCTk745ZoaGhJCTITYiEEELUjqIjRzjcvQfHrhuqy/YbOTSka0Eh/gUFumy/rqt2UZOXl0f79u159913K52/dOlSpkyZwvTp09m9eze9e/dm6NChxMXFWZeJjo4mKiqqwiMxMfGi267svhBylZIQQoja4hAYiDkjwzqwZW0b6NOLj06ncG1WYa1v2xZUu0/N0KFDGTr0whXq3LlzGT9+vLUD77x58/j11195//33mTVrFgA7d+68rLAhISHlWmZOnTpF165dK122qKiIoqIi6/OsrCwAsrOzL2vbQgghBJpGgbc35vR00v76C5fWkbW6eYtbENlFCmNxQr35PTu3n1W64aW6AoBavny59XlRUZEyGo3qm2++KbfcpEmTVJ8+faq17rVr16qbb7653LSSkhLVrFkzderUKZWdna2aNWum0tLSKn39jBkzFCAPechDHvKQhzzs4BEfH3/J2qFGr35KS0vDbDYTFBRUbnpQUBCnT5+u8nqGDBnCrl27yMvLIzQ0lOXLl9OlSxccHByYM2cO/fv3x2Kx8PTTT+Pn51fpOp599lmmTp1qfW6xWDhz5gx+fn7lTll16dKF7du3XzLTxZbLzs4mLCyM+Ph4PD09q7yftqCq748tbr8m1n2566ju66qzfFWWra/HM+h7TNvr8Xw5r62J796qLmPPx7S9fkf/c71KKXJycggODr7k667KJd3n93NRSlWr78vFrmi68cYbufHGGy+5DicnJ5ycnMpN8/b2rrCc0Wis0kFeleU8PT3t7n+Yqr4/trj9mlj35a6juq+rzvJVWba+Hs+g7zFtr8fz5by2Jr97q7ouezym7fU7+vz1enl5Vel1NXpJt7+/P0ajsUKrTEpKSoXWm7ri0UcfrdHl7I3e+301t18T677cdVT3ddVZvirL6v256knPfbfX4/lyXluT371yPNvf9i93vdrZvjGX92JNY/ny5YwcOdI6rWvXrkRHRzN//nzrtNatWzNixAhrR2F7k52djZeXF1lZWXb3V4Cof+R4FvZGjun6o9qnn3Jzczl69Kj1eWxsLHv27MHX15dGjRoxdepU7r77bjp37kz37t1ZsGABcXFxTJgwoUaD1yVOTk7MmDGjwukuIWyRHM/C3sgxXX9Uu6Vm3bp19O9f8dbMY8eOZcmSJUDZzfdmz55NUlISUVFRvPnmm/Tp06dGAgshhBBCVOaKTj8JIYQQQtQVdXLsJyGEEEKI6pKiRgghhBB2QYoaIYQQQtgFKWqEEEIIYRekqKllo0aNwsfHh1tuuUXvKEJclh9//JGWLVvSvHlzPvroI73jCHFF5DvZvsjVT7Vs7dq15Obm8vHHH7Ns2TK94whRLaWlpbRu3Zq1a9fi6elJp06d2Lp1K76+vnpHE+KyyHeyfZGWmlrWv39/PDw89I4hxGXZtm0bbdq0ISQkBA8PD4YNG3bRsdqEqOvkO9m+SFHzD+vXr2f48OEEBwejaRrffvtthWXmz59PREQEzs7OREdHs2HDhtoPKsRlutJjPDExkZCQEOvz0NBQEhISaiO6EBXId7Y4nxQ1/5CXl0f79u159913K52/dOlSpkyZwvTp09m9eze9e/dm6NChxMXFWZeJjo4mKiqqwiMxMbG2dkOIC7rSY7yys9Wapl3VzEJcSE18Zws7o0SlALV8+fJy06655ho1YcKEctNatWqlnnnmmWqte+3atermm2++0ohCXJHLOcY3btyoRo4caZ03adIk9dlnn131rEJcypV8Z8t3sv2QlpoqKi4uZufOnQwePLjc9MGDB7Np0yadUglRc6pyjF9zzTX89ddfJCQkkJOTw4oVKxgyZIgecYW4KPnOrp+qPUp3fZWWlobZbCYoKKjc9KCgIE6fPl3l9QwZMoRdu3aRl5dHaGgoy5cvp0uXLjUdV4hqq8ox7uDgwJw5c+jfvz8Wi4Wnn34aPz8/PeIKcVFV/c6W72T7IkVNNZ3ff0ApVa0+BXKliKjrLnWM33jjjdx44421HUuIy3Kp41m+k+2LnH6qIn9/f4xGY4VWmZSUlAp/CQhhi+QYF/ZEjuf6SYqaKjKZTERHR7Nq1apy01etWkWPHj10SiVEzZFjXNgTOZ7rJzn99A+5ubkcPXrU+jw2NpY9e/bg6+tLo0aNmDp1KnfffTedO3eme/fuLFiwgLi4OCZMmKBjaiGqTo5xYU/keBYV6Hz1VZ2ydu1aBVR4jB071rrMe++9p8LDw5XJZFKdOnVSv//+u36BhagmOcaFPZHjWZxPxn4SQgghhF2QPjVCCCGEsAtS1AghhBDCLkhRI4QQQgi7IEWNEEIIIeyCFDVCCCGEsAtS1AghhBDCLkhRI4QQQgi7IEWNEEIIIeyCFDVCCCGEsAtS1AghhBDCLkhRI4QQQgi7IEWNEEIIIezC/wOZ0bxpMQOfdAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(ETeV, atten, label=\"Input\")\n",
    "plt.loglog(ETeV, template(ETeV * u.TeV), ls='--', label=\"Template spectral model\")\n",
    "plt.loglog(ETeV, optdepth(ETeV * u.TeV), ls='-.', label=\"OptDepth spectral model\")\n",
    "plt.loglog(ETeV, optdepth.evaluate(ETeV * u.TeV, z, 2.), ls='-.', label=r\"OptDepth spectral model, $\\alpha = 2$\")\n",
    "plt.ylim(1e-10, 2.)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00c540c3-4844-49ee-ab87-309f29729ef3",
   "metadata": {},
   "source": [
    "Now, you can use the `OptDepthSpectralModel` within a `gammapy` `CompoundSpectralModel` to fit the EBL normalization `alpha_norm` (or the source redshift) together with the intrinsic spectral parameters to your favorite blazar observation."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:gammapy-1.0]",
   "language": "python",
   "name": "conda-env-gammapy-1.0-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.9.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}