{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ChiSquare Test:\n",
    "This Python macro illustrates the use of ChiSquare as a goodness-of-fit measure,\n",
    "how the ChiSquare distribution comes about, how to use the ChiSquare value (along with\n",
    "number of degrees of freedom, Ndof) as a test, and that it actually works!\n",
    "\n",
    "The program generates a certain number of datasets, each consisting of 9 points along\n",
    "a line. These are then fitted, and the result and the Chi2 of the fit is recorded along\n",
    "with the probability of the fit.\n",
    "\n",
    "## References:\n",
    "- Barlow: Chapter 6\n",
    "- Cowan: Chapter 2.7, Chapter 7\n",
    "- Bevington: Chapter 6\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   12th of November 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Initial Question/Task:\n",
    "Make sure you've read the relevant references and that you understand not only what the ChiSquare is, but also that it follows the ChiSquare distribution, and that the probability of obtaining such a ChiSquare or worse (assuming that the hypothesis/PDF is the correct one!) can be calculated from it.\n",
    "\n",
    "Executive Summary: The ChiSquare method is used both for getting a good and robust fit, and for testing the quality/probability of this fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Program settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "save_plots = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "r = np.random                         # Random generator\n",
    "r.seed(1)                            # Set a random seed (but a fixed one - more on that later.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nexp = 1000\n",
    "Npoints = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Parameters used in this problem (constant, linear coef., quadratic coef., and uncertainty in y) for generating a line (parabola) with N points and given uncertainty on points:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha0 = 3.6\n",
    "alpha1 = 0.3\n",
    "alpha2 = 0.0     # Initially, we don't include a quadratic term!\n",
    "sigmay = 1.0     # The uncertainty on each point is 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating and fitting data:\n",
    "In the following, we generate data points following a simple line. The numerical fit of the line (and any other function) is done iteratively by Minuit. The whole thing is contained in a loop, allowing us to repeat the generation and fitting, so that we can see the result on more than just a case-by-case. For this reason, we save the result of the fit (in pre-defined arrays)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Fit: a0= 3.716+-0.726  a1=0.261+-0.129  p=0.0500\n",
      "------------------------------------------------------------------\n",
      "| FCN = 8.985                   |      Ncalls=32 (32 total)      |\n",
      "| EDM = 1.49e-22 (Goal: 0.0002) |            up = 1.0            |\n",
      "------------------------------------------------------------------\n",
      "|  Valid Min.   | Valid Param.  | Above EDM | Reached call limit |\n",
      "------------------------------------------------------------------\n",
      "|     True      |     True      |   False   |       False        |\n",
      "------------------------------------------------------------------\n",
      "| Hesse failed  |   Has cov.    | Accurate  | Pos. def. | Forced |\n",
      "------------------------------------------------------------------\n",
      "|     False     |     True      |   True    |   True    | False  |\n",
      "------------------------------------------------------------------\n",
      "  Fit: a0= 3.651+-0.726  a1=0.233+-0.129  p=0.2537\n",
      "  Fit: a0= 3.932+-0.726  a1=0.282+-0.129  p=0.6989\n",
      "  Fit: a0= 3.180+-0.726  a1=0.296+-0.129  p=0.9727\n",
      "  Fit: a0= 3.543+-0.726  a1=0.343+-0.129  p=0.2892\n",
      "  Fit: a0= 4.325+-0.726  a1=0.174+-0.129  p=0.5569\n",
      "  Fit: a0= 3.943+-0.726  a1=0.328+-0.129  p=0.8759\n",
      "  Fit: a0= 4.463+-0.726  a1=0.163+-0.129  p=0.1152\n",
      "  Fit: a0= 3.602+-0.726  a1=0.313+-0.129  p=0.4949\n",
      "  Fit: a0= 3.436+-0.726  a1=0.370+-0.129  p=0.9796\n"
     ]
    }
   ],
   "source": [
    "# Arrays for storing fit results:\n",
    "array_alpha0 = np.zeros(Nexp)\n",
    "array_alpha1 = np.zeros(Nexp)\n",
    "array_Chi2 = np.zeros(Nexp)\n",
    "array_Prob = np.zeros(Nexp)\n",
    "\n",
    "# Loop, repeating the data generation and fit:\n",
    "for iexp in range(Nexp) : \n",
    "\n",
    "    # Generating data by filling values into (x,y) and associated uncertainties (here chosen to be 0 for x):\n",
    "    x = np.arange(Npoints)+1\n",
    "    ex = np.zeros_like(x)\n",
    "    y = alpha0 + alpha1*x + r.normal(0, sigmay, Npoints) + alpha2*x**2\n",
    "    # Note how we include uncertainty in y - by simply adding a Guassian number to the \"theoretical value\".\n",
    "    ey = sigmay*np.ones_like(x)\n",
    "\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Numerical fit:\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Define a fit function:\n",
    "    def fit_function(x, alpha0, alpha1):\n",
    "        return alpha0 + alpha1*x\n",
    "    \n",
    "    # Then define Chi2 calculation:\n",
    "    def chi2_owncalc(alpha0, alpha1) :\n",
    "        y_fit = fit_function(x, alpha0, alpha1)\n",
    "        chi2 = np.sum(((y - y_fit) / ey)**2)\n",
    "        return chi2\n",
    "\n",
    "    fitoutput = int(iexp == 1)        # We only want the fit output for the first fit.\n",
    "    \n",
    "    # Here we let Minuit know, what to minimise, how, and with what starting parameters:   \n",
    "    minuit = Minuit(chi2_owncalc, pedantic=False, alpha0=3.0, alpha1=0.0, print_level=fitoutput)     # Own alternative\n",
    "\n",
    "    # Perform the actual fit (and save the parameters):\n",
    "    minuit.migrad();                                                  \n",
    "    minuit_output = [minuit.fmin, minuit.params]    # Save the output parameters in case needed\n",
    "    \n",
    "    # Here we extract the fitting parameters and their errors\n",
    "    alpha0_fit = minuit.values['alpha0']\n",
    "    alpha1_fit = minuit.values['alpha1']\n",
    "    sigma_alpha0_fit = minuit.errors['alpha0']\n",
    "    sigma_alpha1_fit = minuit.errors['alpha1']\n",
    "    \n",
    "    Nvar = 2                     # Number of variables (alpha0 and alpha1)\n",
    "    Ndof_fit = Npoints - Nvar    # Number of degrees of freedom = Number of data points - Number of variables\n",
    "    \n",
    "    # In Minuit, you can just ask the fit function for the Chi2:\n",
    "    Chi2_fit = minuit.fval                          # The chi2 value\n",
    "    Prob_fit = stats.chi2.sf(Chi2_fit, Ndof_fit)    # The chi2 probability given N degrees of freedom\n",
    "    \n",
    "    # Fill the arrays with fit results (to produce plots of these at the end):\n",
    "    array_alpha0[iexp] = alpha0_fit\n",
    "    array_alpha1[iexp] = alpha1_fit\n",
    "    array_Chi2[iexp] = Chi2_fit\n",
    "    array_Prob[iexp] = Prob_fit\n",
    "    \n",
    "    # Let us see what the fit gives for the first couple of data sets:\n",
    "    if (iexp < 10) :\n",
    "        print(f\"  Fit: a0={alpha0_fit:6.3f}+-{sigma_alpha0_fit:5.3f}  a1={alpha1_fit:5.3f}+-{sigma_alpha1_fit:5.3f}  p={Prob_fit:6.4f}\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit the data and plot the result:\n",
    "Below we plot the latest fit (i.e. the last one from above), so that we can actually see what is going on (and potentially wrong!)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de2yd530n+O9DybJN32PLl9gR6UtiO7Edy6blpkbT3SSTaTuJkxppkEbpJZmJcmuaG5pJV1iYXsA7EyBYdIBiihWSne2gmhRN2mCKYCZogLbbHWArmfIlvseJI/luK45vsnyRxWf/eCmd94iSTdIkXx7y8wEIiudhjn95cEh99Z7f+/xKrTUAAEBjqOsCAABgKRGQAQCgRUAGAIAWARkAAFoEZAAAaFm9EE962mmn1dHR0YV4agAAmBc7duz4ea117aGPL0hAHh0dzcTExEI8NQAAzItSyq7DPa7FAgAAWgRkAABoEZABAKBFQAYAgBYBGQAAWgRkAABoEZABAKBFQAYAgBYBGQAAWgRkAABoEZABAKBFQAYAgBYBGQAAWgRkAABoEZABAKBFQAYAgBYBGQBIkoyPj6eUcvBjfHy865KgE6XWOu9POjY2VicmJub9eQGAhVdKyULkA1hqSik7aq1jhz7uCjIAALQIyAAALKql3s6jxQIA6KPFgsXS9WtNiwUAAMyAgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAywhI2Pj6eUcvBjfHy865IAlr1Sa533Jx0bG6sTExPz/rwAK1UpJQvx+xoOx+uNxdL1a62UsqPWOnbo464gAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAy2sG5FLKhaWUW1sfz5ZSvrgYxQEAwGJb/VrfUGu9N8nlSVJKWZXk4STfW+C6AACgE7NtsXh3kp/WWnctRDEAANC12QbkjyT59uEWSimbSikTpZSJ3bt3v/7KAACgAzMOyKWUNUmuTfKdw63XWrfUWsdqrWNr166dr/oAAGBRzeYK8q8nubnW+vhCFQMAAF2bTUD+7RyhvQIAAJaLGQXkUspwkn+R5G8WthwAAOjWax7zliS11r1JTl3gWgAAoHMm6QEAQIuADAAALQIyAAC0CMgAANAiIAMAQIuADAAALQIyAAC0CMgAANAiIAMAQIuADAAALQIyAAC0CMgAANAiIAMAQIuADAAALQIyANDYujUZHc3+JBkdbb6GFUhABhbF+Ph4SikHP8bHx7suCWjbujXZtCnZtasJB7t2NV8LyaxApdY67086NjZWJyYm5v15gcFXSslC/N5Z7uwbC250tAnFhxoZSXbuXOxqWO62bk02b87krl0ZGhlJbrwx2bhx0csopeyotY4d+vjqRa8EAFh6Hnhgdo/DXB14t2Lv3v53K5JOQvLhaLEAAJJ162b3OMzV5s3J3r39j+3d2zy+RAjIAEDzFvfwcP9jw8PN4zCfBuDdCgEZAGje2t6yJRkZyWTS9B5v2bJk3vJmGRmAdysEZACgsXFjsnNnViXNjXnCMQthAN6tEJABAF4Hx1jO0gC8W+GYN2BROa5sbuwbi8nrbW7s2+x1vWdHOubNFWQAAGgRkAEAoEVABgCAFgEZAABaBGQAAGgRkAEAoEVABgCAFgEZAABaBGQAAGgRkAEAoEVABgCgE8NJMjnZdRnTrO66AAAAVoBXXknuuivZvv3gxzNJcuedyaWXdl1dH1eQAZayrVuT0dHsT5LR0eZrgKWu1mTXruQ730n+6I+Sd74zOemk5O1vTz75yeS7303OPDP/e9I8vsSUWuu8P+nY2FidmJiY9+cFBl8pJQvxe2dZ2ro12bQp2bu399jwcLJlS7JxY3d1sez5OZ2bFb1vTz2V3HRTsm1b7wrxE080a0cfnaxfn2zYkFx9dfP5/POTUjrfs1LKjlrr2LTHl0tAHh8fzw033HDw6+uvvz7j4+OLWgPw2rr+ZThQRkebKzCHGhlJdu5c7GpYQfyczs2K2bcXX0xuvbWvVSL33deslZJcfHETgg98XHppsmbNYZ+q6z1b9gH5gK43Gnh1fkZnYWioeZvyUKUsyZtaWD78nM7Nsty3ycnk3nv7w/BttyX79jXrb3xj76rwhg3JlVfOqmWi6z07UkB2kx7AUrVu3eGvIK9bt/i1ACvDI4/0h+GbbkqefbZZO+GE5Kqrkq98pReIzz6723oXiIAMsFTdeOPhe5BvvLG7moDl49lnkx07emF427bk4YebtdWrmxvqNm7s9Q5feGHzztYKICADLFUHbsTbvDmTu3ZlaGSkCcdu0ANma9++5Pbbe0F4+/bk7rt7bVwXXJD86q/2wvDllyfHHNNtzR3SgwwsKj+jc2PfWExeb3OzZPat1uT++/tPlLjllubmuiRZu7b/RImxseTUUzsptes904MMALAc7d7d3ze8fXvyi180a8PDzY1zn/tcr294ZKS52ZcjEpABAAbF3r3JzTf39w0fOPZxaCi55JLkuut6Yfhtb2v6iZkVOwYAsBTt398bzXygXeKOO5rHk+ZK8IYNvavDV16ZHHdctzUvEwIyAEDXak0efLC/b3jHjuT555v1k09uQvD739+7OnzGGd3WvIwJyAAAi+3AaOZ23/Djjzdra9Y0o5k/8YnejXQXXKBveBEJyAAAC+nFF5vpc+0w/OMf99Yvuij5l/+yF4Yvu+yIo5lZHAIyAMB8mZxswm/7Jrr2aOazzmqC8O//fu+ItVmMZmZxCMgAAHP16KMHg/APk+SUU3qjmY8/vhnN/OUv949m1iqx5AnIAAAz8dxzycREf6vEQw81a6tX55Qk+ehHe2H4oouSVau6rJg5EpABAA7VHs184OOuu3qjmc8/P/mVX+mF4fXrMzY8nPpnf9Zt3cwLARkAWNkOjGZuh+Gbb+6NZj7ttCYE/9ZvNf3DV13V2WhmFoeADACsLLt3945YO3Du8IHRzMce2wzc+Oxne1eHR0f1Da8wAjIAsHzt3Zvcckt/GP7Zz5q1oaFmFPNv/mb/aOajjuq2ZjonIAMAy8P+/cndd/dPo7v99t5o5nXrmhD8mc/0RjMff3y3NbMkCcgAwOA5MJq53Tc8MdEbzXzSSU0I/trXeleHzzyz25oZGAIyALD0Pf309NHMjz3WrK1Zk1x+efLxj/ePZh4a6rZmBpaADAAsLS+91D+aedu2/tHMF16YvPe9vSvDl12WHH10d/Wy7AjIAEB3JieT++7rD8O33tobzXzmmc1V4d/7vd5o5pNP7rZmlj0BGQBYPI891n+ixE03Jc8806wdf3wTgL/0pd7V4XPOccQai05ABgAWxp4900czP/hgs7ZqVdMa8ZGP9MLwxRcbzcySICADAK/fvn3JHXdMH808Odmsn3decs01TRC++urmprrh4W5rhiMQkAGA2am1GbZx6GjmF15o1k89tQnCH/pQ8/mqq5pxzTAgZhSQSyknJ/lmkkuS1CSfqLX+fwtZGACwRPz859OPWPv5z5u1Y45Jrrgi+fSne60S556rb5iBNtMryP8hyQ9qrR8qpaxJ4j0RAFiOXngh70iSP/mT3o1099/frJXSjGK+9tpeGL7kEqOZWXZKrfXVv6GUE5PcluS8+lrfPGVsbKxOTEzMQ3mzV0rJDMsEOuBndG7sGwti//7knnv6RzP/6Ee90cxvelMvCB8YzXzCCd3WvIT5OZ29rveslLKj1jp26OMzuYJ8XpLdSf5TKeXtSXYk+UKt9flD/gObkmxKknXr1r3+igGA+VNr8tBD00cz79nTrJ90UtMr/NWv5gP/7t/lvz7ySHLWWd3WDB2ZyRXksST/nOSaWuu2Usp/SPJsrfV/PdL/xhVk4Ej8jM6NfWPWnn56+hFrjz7arB11VHOKxIErw1dfnbz5zQdHM3u9zY19m72u9+z1XEF+KMlDtdZtU19/N8nX5rM4AOB1eOmlpjWiHYbvuae3/pa3JO95Ty8Qv/3tRjPDq3jNgFxrfayU8mAp5cJa671J3p3kroUvDQCYZnIy+clP+sPwLbckL7/crJ9xRnNF+GMf641mPuWUbmuGATPTUyw+n2Tr1AkW9yf5+MKVBAAc9Pjj/WF4+/amfSJJjjuuCcBf+ELv6vCb3uSINXidZhSQa623JpnWnwEAzKM9e5qBG9u3906WeOCBZm3VquTSS5MPf7gXht/6VqOZYQGYpAcAXXjlleTOO/vD8J139kYzn3tu8o53NFeHr746Wb/eaGZYJAIyACy0WpNdu/rPG96xozea+Q1vaK4IX3ddbzTz2rXd1gwrmIAMwLIzPj6eG2644eDX119/fcbHxxevgCefnD6aeffuZu3AaOZPfarXKnHeefqGYQl5zXOQ58I5yMCR+BmdG/s2N4uyby+8kNx6a38Y/slPDhTQ9Am3p9FdeumSH83s9TY39m32ut6z13MOMgCQNP3B99zTH4Zvu63pJ06Sc85pQvC/+Te90cwnnthtzcCsCcgAcCQPP9wfhm+6KXnuuWbtxBObXuE/+qPe1eE3vrHbeoF5ISADQJI8+2z/aOZt25JHHmnWjjqqmT73O7/TnCixYUMznW5qNDOwvAjIAKw8L7+c3H57fxi+557mtIkkefObk3e9q3808zHHdFszsGgEZACWt1qTn/60/7zhW25JXnqpWT/99Oaq8Ec/2hvN/IY3dFsz0CkBGYDl5Yknml7hbdvy35Pk1FOTp55q1oaHmwD8+c/3rg6vW+eINaCPgAzA4Hr++d5o5gMfO3c2a0NDOStJPvSh/tHMq/3VB7w6vyUAGAyvvJLcdVd/GL7jjmT//mZ9dLQJwX/wBwdHM19+/PGpW7Z0WjYweARkYHFs3Zps3pz9SRNkbrwx2bix46JYsmpNHnig/ya6HTuSvXub9VNOacLwBz7QG818+und1gwsGwIysPC2bk02bUr27s1Qkuza1XydCMk0nnqqN5r5wI10TzzRrB19dDOa+ZOf7LVKnH++vmFgwRg1DbM0Pj6eG2644eDX119/fcbHx7sraBCMjjah+FAjI71+UV7Vsvrd9uKLzfS5A0F4+/bkvvuatVKSiy+ePpp5zZo5/aeW1b4tIvs2N/Zt9rresyONmhaQYY681mZhaKh3vmxbKc3oXl7TwL7eJieTe++dPpp5375m/eyz+8Pw2Ni8jmYe2H3rmH2bG/s2e13v2ZECshYLYOGtW3f4K8jr1i1+LSysRx6ZPpr52WebtRNOaHqFv/KVXiA+++xu6wU4DAEZWHg33niwB/mg4eHmcQbXs882N861A/FDDzVrq1c30+c2buyNZr7wQqOZl7h2C1kpRQsZK5YWC5gjr7VZmjrFYnLXrgyNjDjFYpY6f73t2zd9NPPdd/daZy64oBeEN2xILr98SYxm7nzfWFG83mav6z07UouFf8oDi2PjxmTnzqxKmhvzhOOl68Bo5m9/O/niF5Nf/uWmL/jKK5PPfCb5/veT885Lbrgh+cEPkiefbG6y+4u/SP7wD5Nf+qUlEY5h0WzdmoyO9o6x3Lq144KWvvHx8ZSpk2hKKUvunQpXkGGOvNbmxr7NzYLu2+7d/W0S27cnv/hFszY83ATj9o10IyMDc8Sa1xsLrnWM5UHDw8mWLS4EDACnWMA881qbG/s2N/O2b3v3Th/N/LOfNWtDQ8kll/SH4be9baBHM3u9seAcYznQnGIBsNLs3z99NPPtt/dGM4+MNCH4s59t+oevuCI57rhua4ZB88ADs3ucgSAgAywHtSYPPtgfhicmkuefb9YPjGZ+//t7o5nPOKPbmmE5cIzlsiQgAwyip55qAnB7NPPjjzdrRx+drF+f/Ot/3WuVuOCCgekbhoHiGMtlSUAGWOpeeilXJcmf/mkvDP/4x731iy9Ofu3XemH4ssvmPJoZmKUDN+I5xnJZcZMezJHX2tzYt9cwOdmE33arxK239kYzn3VW/3nDY2PJSSd1W/MS5vXGYvJ6Gzxu0gNYih59dPpo5meeadaOP77pFf7yl3Pd17+ev3noIaOZARaBgAywWJ57rn8087Zt/aOZL7ss+e3f7h/NvGpVkuR7X/+6cAywSARkgIWwb19yxx39Yfiuu3qjmc8/P3nnO/tHMx97bLc1A5BEQAZ4/Wpthm20T5S4+ebkxReb9dNOa64Kf/jDvSPWTj2125oBOCIBGWC2fv7z6aOZn3yyWTv22GY082c/27s6PDrqiDWAASIgA7yavXuTW27pD8P339+sDQ01o5g/+MEmCF999cCPZgZAQAbo2b8/ufvu/jD8ox/1RjOvW9cE4U9/ujea+fjju60ZgHm3fALy1q3J5s3ZnzRvZzqkG3g1tTYnSBw6mnnPnmb95JObXuE//uNe3/CZZ3ZbMwCLYnkE5K1bD455HEqameibNjVrQjKQJE8/PX0082OPNWtr1jSjmT/+8f7RzEND3dYMQCeWxyS90dEmFB9qZCTZuXPx6mBFMTFpbhZl3156qWmNOBCEt29P7r23t37RRb0gfGA089FHL2xNr5PX29zYNxaT19vgWd6T9B54YHaPA8vH5GRy333TRzO//HKzfuaZTb/w7/5ubzTzySd3WzMAS9ryCMjr1h3+CvK6dYtfC7CwHnts+mjmp59u1o4/vgnAX/xi71SJs88e6CPWxsfHc8MNNyRprk5df/31GR8f77YogGVuebRYtHqQDxoeTrZs0YPMgvFW2tzMat/27Okfzbx9e++doVWrmtaIA2OZN2xoWiemRjND4ueUxeX1NniWd4vFgRC8eXMmd+3K0MiIUyxg0LzySm8084He4bvualookmY08zXXJF/6UhOG1683mhmABbE8riC3+Ncbi8VrbW5KKamTk80NtIeOZn7hheabTj21/8rwVVc145phlvycspi83gbP8r6CDCxtTz7Z9Apv25bvJ8nppzfjmpPkmGOa0cyf/nQvEJ977kD3DQMw2ARkYH698ML00cw//WmzVkpGkuTaa3th+JJLkqOO6rJiAOgjIANzt39/cs8900czv/JKs/6mNzUh+FOfaj5fcUUuPfHE1G99q9u6AeBVCMjAzNSaPPxwLwhv29Y/mvmkk5pe4X/7b3t9w2ed1W3NADAHAjJweM88M30086OPNmtr1iSXX578/u/3WiXe/GajmQFYFgRkoJk6d+ho5nvu6a1feGHynvf0wvDb377kRzMDwFwJyLDS1Dp9NPMtt/RGM59xRnPE2sc+1muVMJoZgBVEQIbl7vHH+8Pw9u290czHHdeMZv7CF3rnDp9zjiPWAFjRBGRYTvbsaQZutMPwrl3N2qpVyaWXJh/+cC8MX3yx0cwAcAgBGQbVK68kd97Zf6rEnXf2RjOfe27yjnc0V4cPjGYeHu62ZgAYAAIyDIJamyvB7RMlduzoH828YUNy3XW9vuG1a7utGQAGlIAMS9GB0cztVondu5u1Y45JrriiN3xjw4bkvPP0DQPAPBGQoWsvvJDcemt/GP7JT5q1UpK3vjV53/t6YfjSS41mBoAFJCDDYpqcnD6a+bbbeqOZzzmnCcGf/GTz+corkxNO6LZmAFhhBGRYSIcbzfzcc83aiSc2IfirX+31Db/xjd3WCwAIyDBvnn22CcDtaXSPPNKsHXVUM5r5d3+31yrxlrcYzQwAS5CADHPx8su5Ikn+43/sH81ca7P+lrck73pX/2jmY47psmIAYIYEZHgttTY3zR0ymnlHknzuc8nppzeDNz760V6rxCmndF01ADBHAjIc6oknpo9mfuqpZm14uBnN/PnP57e+8Y18Z9eu5E1vcsQaACwjAjIr2/PP949m3ratfzTzJZckH/pQ/2jm1c2PzXe/8Y1k3boOiwcAFoKAzMrxyivJXXf1T6O7447eaObR0eSXfin5wz/sjWY+7rhOSwYAFp+AzPJUa/LAA9NHM+/d26y/4Q1NCP7gB3t9w6ef3m3NAMCSICCzPPziF9NHMz/xRLN29NHNaOYDwzc2bEjOP1/fMABwWDMKyKWUnUmeS7I/ySu11rGFLApe1YsvTh/NfN99zVopTZ/wb/xG/2jmNWu6rRkAGBizuYL8P9daf75glcDhTE4m9947fTTzvn3N+tlnNyH4E59oPo+NNRPqAADmSIsFS8sjj/SH4ZtuaibUJckJJzS9wl/5Su/q8Nlnd1svALDszDQg1yR/V0qpSf7PWuuWBayJleLZZ5sb59o30j38cLO2enUzfW7jxiYIX311cuGFRjMDAAtupgH5mlrrI6WU05P8sJRyT631n9rfUErZlGRTkqxzNiyH2rcv+dGP+q8O3313bzTzBRckv/qrvfOGL7/caGYAoBMzCsi11kemPj9RSvlekg1J/umQ79mSZEuSjI2N1Xmuk0FSa/LTn/aH4ZtvTl56qVlfu7YJwR/5SK9v+NRTu60ZAGDKawbkUspxSYZqrc9N/fm9Sf63Ba+MwbF79/TRzL/4RbM2PJxceWXyB3/Q6xseGXHEGgCwZM3kCvIZSb5XmkCzOsl/qbX+YEGrYunau3f6aOadO5u1oaFmNPN11/XC8NvednA0MwDAIHjN5FJrvT/J2xehFpaa/fsPP5p5//5mfWSkCcGf+1zz+YorkuOP77ZmAIDXyaU9GrUmDz7YC8IHRjM//3yzfvLJTQh+//t7V4fPOKPbmgGOZOvWZPPm7E+S0dHkxhubU3EAZkBAXqmeeiq56ab8/b//99nzD/+QDUnOPLC2Zk2yfn1v+MbVVzenTOgbBgbB1q3Jpk3J3r0ZSpJdu5qvEyEZmJFS6/wfODE2NlYnJibm/XlnopSShfj/NNBeemn6aOYf/7i3fuGF+fN7783v/emfNmH4ssuMZp4Br7W5sW8suNHRJhQfamSkd88ELAC/3wZPKWVHrXVs2uMC8jIzOdmE33YYvvXW3mjmM89sQvCB84bHxpKTTrJvc2DP5sa+seCGhnpnrLeV0vyOhAXi99vgOVJA1mIx6B59tP9EifZo5uOPb0Yzf/nL/aOZtUoAy9m6dYe/gmyIFTBDAvIgee656aOZH3qoWVu9ummN+OhHe2H4oouSVau6rRlgsd1448Ee5IOGh5vHAWZAQF6q9u1rjlRrnypx1129tw3PPz/5lV/pheH165Njj+22ZoCl4MCNeJs3Z3LXrgyNjDjFApgVPchLQa3Jz37WH4Zvvjl58cVm/bTTekH46qubtol5Hs08kPvWMXs2N/aNxeT1xmLyehs8epCXkt27m17h9o10Tz7ZrB17bDOa+bOf7YXi0VF9wwAAi0RAXmh79ya33NIfhu+/v1kbGmpGMX/wg/2jmY86qtuaAQBWMAF5Pu3fn9x9d/+pErff3hvNvG5dE4I//enm85VXGs0MALDECMhzVWtzgkT7RImJid5o5pNOakLw177W6xs+88xXf04AADonIM/U0083Abh9I91jjzVra9Ykl1/eG828YUMzmnloqNuaWRhbtyabN2d/0vSHuzseAJYVAflwXnopue22/r7he+/trV90UfLe9/bC8GWXJUcf3V29LJ6tWw+erzqUNMMINm1q1oRkAFgWHPM2OZncd9/00cwvv9ysHxjNfCAMj40lJ5+8MMV3yNE0MzQ6evgJXSMjyc6di13NQPJaYzF5vbGYvN4Gj2PeDnjssf4wfNNNTftE0twwNzaWfPGLvTOHjWam7YEHZvc4ADBwlndA3rNn+mjmBx9s1latalojPvIRo5mZuXXrDn8Fed26xa9lwIyPj+eGG25I0lxluf766zM+Pt5tUQBwGMunxWL//uT227Np/fps+cQneqOZJyeb9fPP7wVho5mn8bbQDLV6kA8aHk62bNGDDEuQ320sJq+3wbP8WyyeeSZZvz5bkuRv/7YJwR/6UPP5qquacc3weh0IwZs3Z3LXrgyNjDjFAgCWmeVzBTlJ/vZvc+4HPpCfTU7qG54l/+qdPXsGS5+fUxaT19vgOdIV5OV1UO+112ZnIhwDADBnyysgAwDA6yQgAwC8DuPj4ylT716XUpzQswwsrx7k6P+ZK/s2e/YMlj4/p8CrWRk9yAAA8DoJyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0LJsAvL4+HhKKUmSUkrGx8e7LQgAgIFUaq3z/qRjY2N1YmJi3p+XhVNKyUK8FpYzewZLn59T4NWUUnbUWscOfXzZXEEGAID5ICADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAEDLjANyKWVVKeWWUsr3F7IgAADo0myuIH8hyd0LVQgAACwFMwrIpZRzkvyrJN9c2HIAAKBbM72C/CdJvppk8kjfUErZVEqZKKVM7N69e16KAwCAxfaaAbmU8r4kT9Rad7za99Vat9Rax2qtY2vXrp23AgEAYDHN5AryNUmuLaXsTPKXSd5VSvmLBa0KAAA68poBudb6x7XWc2qto0k+kuTva60fW/DKAACgA85BXum2bk1GR7M/SUZHm68BAFaw1bP55lrrPyb5xwWphMW3dWuyaVOyd2/zL6Vdu5qvk2Tjxi4rAwDojCvIK9nmzcnevf2P7d3bPA4AsEIJyCvZAw/M7nEAgBVAQF7J1q2b3eMAACuAgLyS3XhjMjzc/9jwcPM4wAAbHx9PKSVJUkrJ+Ph4twUBA6XUWuf9ScfGxurExMS8Py8LYOvWZPPmTO7alaGRkSYcu0FvRkopWYifHwBgcZRSdtRaxw59fFanWLAMbdyYbNyYVULWqwIAAAfMSURBVKWk7tzZdTUAAJ3TYgEAAC0CMgAAtAjIAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAyzZIQtACxvRk2TxNhkAGDlOdKoaVeQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAIAWARkAAFoEZAAAaBGQAQCgRUAGAICW1wzIpZRjSinbSym3lVLuLKXcsBiFAQBAF1bP4HteSvKuWuueUspRSf5HKeW/11r/eYFrAwCARfeaAbnWWpPsmfryqKmPupBFAQBAV2bUg1xKWVVKuTXJE0l+WGvddpjv2VRKmSilTOzevXu+6wQAgEUxo4Bca91fa708yTlJNpRSLjnM92yptY7VWsfWrl0733UCAMCimNUpFrXWp5P8Y5JfW5BqAACgYzM5xWJtKeXkqT8fm+Q9Se5Z6MIAAKALMznF4qwkf15KWZUmUP9VrfX7C1sWAAB0YyanWPwoyfpFqAUAADpnkh4AALQIyCvc+Ph4SilJklJKxsfHuy0IAKBjpZkDMr/GxsbqxMTEvD8vAADMl1LKjlrr2KGPu4IMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAi4AMAAAtAjIAALQIyAAA0CIgAwBAS6m1zv+TlrI7ya55f+KZOS3Jzzv6bw8y+zZ79mxu7Nvc2Le5sW9zY9/mxr7NXtd7NlJrXXvogwsSkLtUSpmotY51XcegsW+zZ8/mxr7NjX2bG/s2N/Ztbuzb7C3VPdNiAQAALQIyAAC0LMeAvKXrAgaUfZs9ezY39m1u7Nvc2Le5sW9zY99mb0nu2bLrQQYAgNdjOV5BBgCAOROQAQCgZdkE5FLK/1VKeaKUckfXtQyKUsqbSin/UEq5u5RyZynlC13XNAhKKceUUraXUm6b2rcbuq5pkJRSVpVSbimlfL/rWgZFKWVnKeX2UsqtpZSJrusZBKWUk0sp3y2l3DP1O+4dXde01JVSLpx6jR34eLaU8sWu6xoEpZQvTf19cEcp5dullGO6rmkQlFK+MLVndy6119qy6UEupbwzyZ4k/7nWeknX9QyCUspZSc6qtd5cSjkhyY4kH6y13tVxaUtaKaUkOa7WuqeUclSS/5HkC7XWf+64tIFQSvlykrEkJ9Za39d1PYOglLIzyVit1QCCGSql/HmS/7fW+s1Sypokw7XWp7uua1CUUlYleTjJ1bXWrgZ/DYRSytlp/h54a631hVLKXyX5b7XW/7vbypa2UsolSf4yyYYkLyf5QZLP1Frv67SwKcvmCnKt9Z+S/KLrOgZJrfXRWuvNU39+LsndSc7utqqlrzb2TH151NTH8viX5gIrpZyT5F8l+WbXtbB8lVJOTPLOJN9Kklrry8LxrL07yU+F4xlbneTYUsrqJMNJHum4nkFwcZJ/rrXurbW+kuT/SfKbHdd00LIJyLw+pZTRJOuTbOu2ksEw1SZwa5Inkvyw1mrfZuZPknw1yWTXhQyYmuTvSik7Simbui5mAJyXZHeS/zTVzvPNUspxXRc1YD6S5NtdFzEIaq0PJ/lGkgeSPJrkmVrr33Vb1UC4I8k7SymnllKGk/xGkjd1XNNBAjIppRyf5K+TfLHW+mzX9QyCWuv+WuvlSc5JsmHqrSJeRSnlfUmeqLXu6LqWAXRNrfWKJL+e5HNTLWUc2eokVyT5s1rr+iTPJ/latyUNjqmWlGuTfKfrWgZBKeWUJB9Icm6SNyY5rpTysW6rWvpqrXcn+XqSH6Zpr7gtySudFtUiIK9wUz20f51ka631b7quZ9BMvW37j0l+reNSBsE1Sa6d6qf9yyTvKqX8RbclDYZa6yNTn59I8r00PXsc2UNJHmq9s/PdNIGZmfn1JDfXWh/vupAB8Z4kP6u17q617kvyN0l+ueOaBkKt9Vu11itqre9M0ya7JPqPEwF5RZu62exbSe6utf4fXdczKEopa0spJ0/9+dg0vxzv6baqpa/W+se11nNqraNp3r79+1qrqyyvoZRy3NRNtJlqE3hvmrcmOYJa62NJHiylXDj10LuTuPl45n472itm44Ekv1RKGZ76e/Xdae7p4TWUUk6f+rwuyXVZQq+71V0XMF9KKd9O8j8lOa2U8lCS62ut3+q2qiXvmiS/k+T2qX7aJPlfaq3/rcOaBsFZSf586i7voSR/VWt1ZBkL5Ywk32v+3s3qJP+l1vqDbksaCJ9PsnWqXeD+JB/vuJ6BMNUL+i+SfKrrWgZFrXVbKeW7SW5O0yJwS5bo+OQl6K9LKacm2Zfkc7XWp7ou6IBlc8wbAADMBy0WAADQIiADAECLgAwAAC0CMgAAtAjIAADQIiADAECLgAwAAC3/P74V9VNAFxwuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.errorbar(x, y, ey, fmt='ro', ecolor='k', elinewidth=1, capsize=2, capthick=1)\n",
    "ax.plot(x, fit_function(x, *minuit.args), '-r')\n",
    "plt.close()\n",
    "fig.tight_layout()\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha0 = 4.460+-0.726\n",
      "alpha1 = 0.220+-0.129\n",
      "Chi2   =   4.2\n",
      "Ndof   =     7\n",
      "Prob   = 0.757\n"
     ]
    }
   ],
   "source": [
    "# Print the results of the fit:\n",
    "print(\"alpha0 = {:5.3f}+-{:5.3f}\".format(alpha0_fit, sigma_alpha0_fit))\n",
    "print(\"alpha1 = {:5.3f}+-{:5.3f}\".format(alpha1_fit, sigma_alpha1_fit))\n",
    "print(\"Chi2   = {:5.1f}\".format(Chi2_fit))\n",
    "print(\"Ndof   = {:5d}\".format(Ndof_fit))\n",
    "print(\"Prob   = {:5.3f}\".format(Prob_fit))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (save_plots) :\n",
    "    fig.savefig(\"FitToLine.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Many experiments\n",
    "In case we have more than one \"experiment\" (i.e. fit), we would like to see how the fit results and the Chi2 values distribute themselves:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAALICAYAAADseNpmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5gU1Zn48e+LqIwMMKJIcAQHFTGItzgRDRujmBgiifd7VCC4ZDercaM+im6ymJhkZxNWo9FV+QUNrlFAglfMJgY0FyUqaJSgIl4QGIlM1OFiQBHO749pZhmYG9jT3TP9/TxPP91Vp6rPW8Xl9Ft1zqlIKSFJkiRJkvKrU74DkCRJkiRJJuiSJEmSJBUEE3RJkiRJkgqACbokSZIkSQXABF2SJEmSpAJggi5JkiRJUgEwQZfaWEQsiIhj8h1HPkXEKRGxNCLWRMRh27jvNRFxV7a3lSRpS7bZttlSvpmgSx9DRCyOiM9vsW5URPxx03JK6cCU0uMtfE9FRKSI6NxGoebbBOCilFJpSum5fAfTlIg4NyLejIj3I+L+iOiZ75gkSdlhm91qBd9mR0SfiHgwIt7K/FlU5DsmKVtM0KUiUAA/IvYGFuQ5hmZFxIHAbcD5QG/g78B/5zUoSVLRsc1ulY3A/wKn5TsQKdtM0KU2tvkV+4g4IiLmRsSqiHg7Iq7LbPb7zHttpkvZURHRKSK+nbmjuyIi7oyIHpt97wWZsnci4jtb1HNNREyPiLsiYhUwKlP3nIiojYjlEXFTROy02feliPhGRCyKiNURcW1E7JvZZ1VETNt8+y2OsdFYI2LniFgD7AA8HxGvNbH/DZnudKsiYl5EfLaJ7TbdtRibuWq+PCIu22KznTL1r850VazcbP9xEfFapuzFiDhls/2+CjyUUvp9SmkN8B3g1Ijo1lgskqSOxza7fbTZKaW3U0r/DTzT6B+k1I6ZoEu5dQNwQ0qpO7AvMC2z/ujMe1mmS9kcYFTmdSywD1AK3AQQEYOou7v7VaAP0AMo36Kuk4DpQBnwC2AD8C1gd+Ao4DjgG1vsMxw4HDgSuAKYmKmjLzAYOKeJ42o01pTSByml0sw2h6SU9m1i/2eAQ4GewN3AvRHRpYltydQzADgeGBcNuyyeCEzJHPeDZM5ZxmvAZ6k7X98F7oqIPpmyA4HnN22YUnoN+BDYv5k4JEkdl2124wqhzZY6LBN06eO7P3OFuzYiamm+W/R6YL+I2D2ltCal9Kdmtv0qcF1K6fXMHd2rgLOjruvb6dTd7f1jSulD4N+BtMX+c1JK96eUNqaU1qaU5qWU/pRS+iiltJi67tyf22Kf/0wprUopLQD+AvwmU/9K4FdAU5PFNBdri1JKd6WU3snE9l/AzsDAZnb5bkrp/ZTSfOAOGv4I+WNK6ZGU0gbgf4BDNqvn3pTSW5lzMhVYBByRKS4FVm5Rz0rAO+iS1HHYZneMNlvqsEzQpY/v5JRS2aYXW1/h3twY6u7IvhwRz0TEl5vZdk/gzc2W3wQ6Uzc+ek9g6aaClNLfgXe22H/p5gsRsX9EPBwRf810ofshdVfmN/f2Zp/XNrJcSuOai7VFEXFZRLwUESszP5h6NBLb5jY/tjcz9W/y180+/x3osulHR6aL4Z83+2E2eLN61gDdt6inO7C6NccgSWoXbLM7RpstdVgm6FIOpZQWpZTOAfYA/hOYHhFd2fpKOsBb1E3Uskk/4CPqGuDlwF6bCiKiBNhty+q2WL4FeBkYkOmudzUQ2380rY61WZmxa1cCZwK7Zn4wrWwhtr5b1PVWK+rZG/h/wEXAbpl6/rJZPQvY7Mp9ROxD3V2BV1r6bklSx2ObvbUCarOlDssEXcqhiDgvInqllDYCtZnVG4Aa6mYk3Wezze8BvhUR/SOilLqr51NTSh9RN07tKxHxmcwkMN+l5UarG7AKWBMRBwD/nLUDaz7WlnSj7odBDdA5Iv6dre9kb+k7EbFL1M28PhqY2op6Nv2oqgGIiNHUXY3f5BfUndPPZn6AfQ+YkVLyDrokFSHb7CbjKoQ2m8y4950zizu3MA5eajdM0KXcGg4syMySegNwdkppXaa72w+AJzJduY4EbqduPNbvgTeAdcDFAJnxZhdTN7HKcuq6Ya8APmim7suBczPb/j9a10C2VpOxtsKvqRsr9wp1Xd/WsUVXv0b8DngVmAVMSCn9pqVKUkovAv8FzKHuLsFBwBOblS8A/om6RH0FdT9Cmuv6KEnq2Gyzt1YQbXbGWuqGp0Fdb4O1rTsEqbBFSo310pHUnmSugNdS1xXujXzH01YiooK6HxM7tvJKvyRJBcU2W1JzvIMutVMR8ZVMl7GuwARgPrA4v1FJkqQt2WZLai0TdKn9Oom6iVbeou75omcnu8RIklSIbLMltYpd3CVJkiRJKgDeQZckSZIkqQB0zncAW9p9991TRUVFvsOQJKldmTdv3t9SSr1au31EfAu4kLpHGc2n7vFHfaibabon8Cxwfkrpw+a+x3ZbkqRt11S7XXAJekVFBXPnzs13GJIktSsR8eY2bFsOfBMYlFJaGxHTgLOBE4DrU0pTIuJWYAxwS3PfZbstSdK2a6rdtou7JEnFqTNQEhGdgV2oez7zMGB6pnwycHKeYpMkqSiZoEuSVGRSStXUPeppCXWJ+UpgHlC72fOKlwHl+YlQkqTiZIIuSVKRiYhdqXvsU39gT6Ar8KVGNm30US8RMTYi5kbE3JqamrYLVJKkImOCLklS8fk88EZKqSaltB6YAXwGKMt0eQfYi7pnNm8lpTQxpVSZUqrs1avV89JJkqQWmKBLklR8lgBHRsQuERHAccCLwGPA6ZltRgIP5Ck+SZKKkgm6JElFJqX0FHWTwT1L3SPWOgETgSuBSyPiVWA3YFLegpQkqQgV3GPWJElS20spjQfGb7H6deCIPIQjSZLwDrokSZIkSQXBBF2SJEmSpAJggi5JkiRJUgEwQZckSZIkqQCYoEuSJEmSVABM0CVJkiRJKgA+Zk1ShzO0ajbVtWsbLSsvK+GJccNyHJEkSbnXXHsItolSITJBl9ThVNeuZXHViEbLKsbNzHE0kiTlR3PtIdgmSoUoq13cI6IsIqZHxMsR8VJEHBURPSPi0YhYlHnfNZt1diT/9E//xLXXXpvvMCRJUgtssyVJbSHbY9BvAP43pXQAcAjwEjAOmJVSGgDMyiwXrFGjRrHTTjtRWlpa/xo/fnyr9z/mmGP42c9+tl1133rrrXznO9/Zrn1z6fvf/z4HHnggnTp14uc//3mDsmXLlnHMMcewyy678KlPfYq//OUv9WXr169nzJgxdOvWjX79+jFt2rQG+95444307t2bnj17ctVVV+XiUCRJ7ZhtdstssyWpfclagh4R3YGjgUkAKaUPU0q1wEnA5Mxmk4GTs1VnW7niiitYs2ZN/eu73/1uvkMqKPvssw833ngjhx9++FZlY8eOZfDgwbzzzjucddZZnHXWWfVl119/PfPnz2fp0qXceeedfO1rX2Pp0qUAPP3001xzzTXMnj2b+fPnM2XKlK1+DEiStCXb7ObZZktS+5LNO+j7ADXAHRHxXET8LCK6Ar1TSssBMu97bLljRIyNiLkRMbempiaLIWXXz3/+cz7zmc8watQounXrxj/8wz/wt7/9DYAf/vCHlJaW8oc//IGLLrqI0tLSrRrDiooKbrrpJg4//HC6du3KiSeeCMDDDz9MaWkpO+64I9/+9re3qnfGjBkMHjyYXXfdlRNOOIHly5fXl82cOZOBAwfSrVs3BgwYwK9//es2PAN1zj33XI477jh23nnnButXrVrFb37zG8aNG0dJSQnf+ta3ePPNN3nhhRcAuPfee/nmN79JWVkZxxxzDEcddRT33Xdffdmpp57KgQceSHl5ORdeeCFTpkxp82NRYRtaNZuKcTMbfQ2tmp3v8CQVMNvsOrbZktS+ZDNB7wx8CrglpXQY8D6t7M6eUpqYUqpMKVX26tUriyFl39NPP825557LihUrWL9+PbfffjsAV199NWvWrOGzn/0sN910E2vWrGHevHlb7X/bbbdx5513snLlSq6++moAvvzlL7NmzRq++tWvNlrfmDFjmDRpEjU1NRx22GGMHTu2vvzCCy/ke9/7HqtXr+Y3v/kN5eXlbXTkLXv11Vfp0qULpaWlfPazn2XZsmXsu+++LFy4EICFCxdywAEHcN555zFjxgwGDRq0VdkNN9zAVVdd1aBMxWvT5DaNvZqblVaSwDa7ObbZklSYspmgLwOWpZSeyixPpy5hfzsi+gBk3ldksc42MWHCBMrKyupfM2bMqC/bf//9Of744ykpKeG4447jlVde2abvHjt2LAceeCCdO3fmyCOPbHH7SZMmMWrUKIYMGULnzp25/PLLmTlzJh988AEAnTp14rXXXmPVqlX079+fwYMHb9vBZtH7779PaWkpq1evZtGiRbz33nt069aNNWvWNChfuHAh1dXVjZYtXbqURYsWNSiTJKkpttnbxzZbkgpT1hL0lNJfgaURMTCz6jjgReBBYGRm3UjggWzV2VYuv/xyamtr61+nnnpqfVnPnj3rP++0006sW7dum757wIAB27T90qVLue222+p/ePTv35+ddtqpvsvcvffey5w5c+jXrx9HHnlkgwlecq1r166sWbOGvn378te//pXDDz+c1atXU1pa2qD8mWee4eKLL260bMKECUyfPr1BmSRJTbHN3j622ZJUmLI9i/vFwC8i4gXgUOCHQBXwhYhYBHwhs9xhderU/Cnt3HnbHj3ft29fvvOd7zT48bFu3ToqKioA+MxnPsNDDz3E22+/zQEHHNDoeLhc2W+//Vi3bh3Lli0D4MMPP+S1115j4MC6azb7778/L7/8cv32L774YqvKJElqC7bZttmSVGiymqCnlP6cGUt+cErp5JTSeymld1JKx6WUBmTe381mnYXmE5/4RFaviI8ePZpbb72VZ599lpQSK1asYOrUqQBs3LiRu+66i9WrV9f/yOjevXvW6m7K+vXrWbduHSml+s8bN26ke/fuHH/88VRVVbFu3Tquv/56+vXrx0EHHQTAmWeeyY033khtbS2/+93vmDNnDqeccgoAZ5xxBjNmzGDBggVUV1czadKkBrPJSpKUbbbZttmSVGiyfQe9Q/jRj37U4JmqF154Yav3veyyy3j00Ufp06cPxxxzTKv2Of744yktLeUXv/hFfd2jRo0C4Mgjj+S6665j9OjR9OjRg8rKygYT2dx111307duX3Xffnbfeeosf/OAH23Ko2+Uf//EfKSkp4cknn2Ts2LGUlJTw+9//HqibUGf+/PnsuuuuTJkyhalTpxIRAHzrW99i8ODB9O3bl/POO49JkybRt29fAIYMGcL48eM59thjGTx4MGeeeaaNvSSpRbbZzbPNlqT2JVJK+Y6hgcrKyjR37tx8hyGpAFSMm8niqhE5K5Pas4iYl1KqzHW9tttS4WqpzbNNlPKnqXbbO+iSJEmSJBUAE3RJkiRJkgrAtk1PKkkd2NCq2VTXrm20rLyshCfGDctxRJIkSSomJuiSlFFdu7bZseuSJElSW7KLuyRJkiRJBcAEfTOPP/44EcGVV15Zv+6II46of+RIsXjsscc46KCD6N69O3369OHSSy9lw4YNrdp37dq1fP3rX6dnz56UlZXxjW98o77sRz/6EQMGDKBbt24MHjyYBx54oMG+N954I71796Znz55cddVVWT0mSVLHYpv9f7an/Vy/fj0XXHABffr0oUePHhx77LEsWLCgvry6upoTTjiBHj16sPfee3P33XfXl6WUuPzyy+nduze77bYb55xzDqtWrcr6cUlSMTJB30KvXr347W9/C8CSJUtYuXJlniPKvQMOOIBHHnmElStX8sorr7BgwQJuvvnmVu37rW99i1dffZUFCxbwt7/9rUGC3rlzZ375y1+ycuVKbrvtNs4//3xef/11AJ5++mmuueYaZs+ezfz585kyZQrTpk1rk+OTJHUMttnb335u2LCB/fbbj2eeeYb33nuPE088kZNPPrm+/KKLLqJv37787W9/4+c//zlf+9rXeOONNwCYOnUq06ZN49lnn2XJkiW8++67XHvttW12jJJUTEzQt7DTTjsxaNAgnnvuOaZPn86pp57aoHzGjBkMHjyYXXfdlRNOOIHly5fXl5199tn1V7CPO+44Fi1aVF9WUVHB+PHj2W+//ejZsye33nprzo5pW/Xp04e+ffsSEXz00Uds3LiRl156qcX91q5dy+TJk/npT39Knz596Ny5M4MHD64vv/TSSzn44IPp1KkTQ4cOZZ999mHevHkA3HvvvZx66qkceOCBlJeXc+GFFzJlypQ2O0a1f+VlJVSMm9noq7ysJN/hScoB2+ztbz+7dOnCv//7v7PXXnvRqVMnRo8ezauvvkpNTQ0Av/vd7/j617/OjjvuyLHHHsuhhx7KQw89BMCbb77JZz7zGcrLy+natStf/vKXefHFF9v0OCWpWJigN+L0009n+vTp3H///Zx00kn1659++mnGjBnDpEmTqKmp4bDDDmPs2LH15Ycddhjz58+npqaGyspKzjrrrAbf+9RTT/HCCy8wadIkrrjiCj766KOcHdO2WrJkCT169KBnz5489dRTjBo1qsV9XnnlFSKC+++/n969ezNo0CDuu+++Rrd97733eOWVV+oT+IULF3LAAQdwww03cNVVVzFo0CAWLlyYzUNSB/PEuGEsrhrR6MvZ1qXiUextdrbazzlz5rDnnnuy2267AXXd2DeXUuLVV18F6i5uvPHGGyxdupQ1a9bw8MMPM2JE4xNsSpK2jQl6I774xS/yy1/+kp133pndd9+9fv2kSZMYNWoUQ4YMoXPnzlx++eXMnDmTDz74AIArr7ySPfbYgx122IFRo0bx/PPPN/jer33ta+yyyy6MGDGC1atX89e//jWnx7Ut+vXrx8qVK3nttde49NJLGTBgQIv7rFq1ig8//JDXX3+dN998k5tvvpnzzz+/wR2LTb7+9a8zcuRIPvnJTwLw/vvvU1paytKlS1m0aBHdunVjzZo1WT8uSVLHUuxtdjbaz5UrV3LJJZdw3XXX0alT3U/DY445httuu40PPviAWbNm8ec//5m///3vAHziE59gyJAh7L333vTo0YMddtihwcUPSdL2M0FvRJcuXTjxxBMZPXp0g/VLly7ltttuo6ysjLKyMvr3789OO+3E8uXL2bBhA1dddRX77rsvZWVlDBkyhI0bNzaYXK1nz55AXZc8gHXr1uXuoLbTPvvswwEHHMBFF13U4ra77LILGzZs4LLLLqNLly4ce+yx7L///vzpT39qsN3VV1/NO++8ww033FC/rmvXrqxZs4YJEyYwffp0Vq9eTWlpadaPR5LUsRR7m/1x288PPviAU045hbPOOqtBL4Kf/vSnLFu2jL322ouqqipOPPFEunfvDsC1117L888/z4oVK1i5ciXdu3fnm9/8ZtaPTZKKkQl6E370ox9x7rnnNljXt29fvvOd71BbW1v/WrduHRUVFdx9993cd999PPbYY9TW1vLHP/4R2LqLWHuUUuKFF15ocbt99tmHiGgwg25KqcE5+MlPfsKvf/1r7r///vofPQD7778/L7/8cv3yiy++yMCBA7N0BJKkzUXEwIj482avVRHxrxHRMyIejYhFmfdd8x1raxRzm/1x2s8NGzZw7rnnst9++201ydtee+3FzJkzqamp4dFHH+WNN97gkEMOAeDZZ5/ltNNOY/fdd6e0tJTzzz+fWbNmZe+gJKmImaBvg9GjR3Prrbfy7LPPklJixYoVTJ06FYDVq1dTUlLCrrvuyurVq/nhD3+Y52i33z333MMLL7zAxo0bWbJkCRMmTOCoo45qsM2MGTPYb7/9qK6url+366678rnPfY7rrruO9evX84c//IFFixZx5JFHAnDnnXdyyy238Ktf/Ypu3bo1+L4zzjiDGTNmsGDBAqqrq5k0adJW4wElSdmRUlqYUjo0pXQocDjwd+A+YBwwK6U0AJiVWW6XiqXNbk372VibDXXDzQBuueWWrb536dKlvPvuu3zwwQf85Cc/YenSpZxyyilA3fj9++67j/fee49169Yxbdo0Bg0a1EZHKEnFxQR9Gxx55JFcd911jB49mh49elBZWVk/C/kFF1xAv3792HPPPTn44IPrk9L26J133uH000+ne/fufPrTn+bggw9mwoQJDbZZtWoVr732GuvXr2+w/vbbb2fhwoWUlZUxZswY/ud//oc999wTgPHjx7N48WL22WcfSktLKS0trf9RNGTIEMaPH8+xxx7L4MGDOfPMM03QJSk3jgNeSym9CZwETM6snwyc3OReBa5Y2uzWtJ+NtdlvvvkmkyZN4le/+hU9evSob5f/8Ic/AHWTz22aAX/q1Kn86le/qu/ifvXVV9O/f38GDhzInnvuybvvvstPf/rT3B20JHVgUWjduSorK9PcuXPzHYakAlAxbiaLq7I7M3Bz37m9ZVIhiIh5KaXK7djvduDZlNJNEVGbUirbrOy9lNJW3dwjYiwwFqBfv36Hv/nmmx8ndEmNGFo1m+ratU2Wl5eVtPjUkpbaLts2KX+aarc75yMYSZKUfxGxE3AicNW27JdSmghMhLoL620QmlT0qmvXtphcS+p4TNAl5VVzdwjKy0pyHI1UdL5E3d3ztzPLb0dEn5TS8ojoA6zIY2ySJBUdE3RJedXSHQJJbeoc4J7Nlh8ERgJVmfcH8hGUJEnFykniJEkqQhGxC/AFYMZmq6uAL0TEokxZVT5ikySpWHkHXZKkIpRS+juw2xbr3qFuVndJkpQH3kHfQkVFBSUlJZSWlrLvvvty00031ZdFBEOGDKlfvuKKK4gIHn/8cQAefPBBDjroILp168bee+/ND37wg0a/d9Prt7/9bc6OK5sef/xxBg4cSNeuXTn55JNZuXJlq/Y78MADGxx/586dufjii+vLb7zxRnr37k3Pnj256qqG8xVFBF27dq3fd+LEiVk9JklS+2Ob3bLtbbMXLlzI8OHDKSsro6KiokHZI488wqc//Wm6d+9O3759+f73v19ftmjRIoYPH07Pnj3ZY489uOCCC1pdpyTJBL1RDz30EGvWrOGee+7hyiuv5He/+119WW1tLUuWLAFg1qxZ9OrVC4A33niDCy64gJtvvplVq1bxhz/8gf3226/R7930+vznP5+7g8qSv//975xxxhmMHz+eFStWEBFbJdNNWbBgQf2xr1q1ij333JPTTjsNgKeffpprrrmG2bNnM3/+fKZMmcK0adMa7P/888/X7z927NisH5s+nqFVs6kYN7PR19Cq2fkOT1IHZZvdtI/TZnfu3Jmzzz6bH//4x1uVrVmzhv/4j/+gpqaGP/3pT9x1113cddddQN0z18855xzeeOMNFi9ezLp167jsssuyelyS1JHZxb0ZRxxxBAceeCDPPPMMn/vc5wA49dRT+eUvf8nnPvc5Bg0axNtv1018O2/ePPr378/RRx8NQL9+/ejXr1/eYm8rjz32GD169ODcc88F4PLLL+fEE0/kv//7v7fpe37729/SqVOn+vN67733cuqpp3LggQcCcOGFFzJlyhTOPPPM7B6A2kxzk735KBhJbc02e2sfp83ed9992XfffRvtObB521xeXs7w4cOZM2cO5513HocffjiHH354ffm5557Lv/3bv2XhaCSpOHgHvRnz5s1jwYIFfPKTn6xfd/LJJ/PAAw8wffr0+ru/AAcddBALFixg/PjxPP/886TUMR8Lu3DhQg444AD++Mc/Mnz4cAYOHMi7775LTU3NNn3PHXfcwXnnnUdENPjeG264gauuuopBgwaxcOHCBvscffTR9OnTh9GjR9tdTpLUgG321rLVZrfkT3/6EwcffHCjZXPmzGmyTJK0NRP0Rpx88snssssuHH300VRVVTFixP/dFdxtt93YeeedmT59Ol/84hfr1w8cOJBHHnmEJ598kiFDhtC/f38eeOCBrb63rKys/rViRft7vOz7779PaWkpNTU1LFy4kC5dugB13d1aq7a2lvvvv58LLrhgq+9dunQpixYtolu3bg2+88knn2TJkiU8++yzvPXWW3zzm9/M3kFJktot2+ymZaPNbsnNN9/MunXrGDVq1FZl8+bN44477uDaa6/NWn2S1NGZoDfi/vvvZ/Xq1fzrv/4rTz755Fblo0eP5sQTT6SkpKTB+s9//vM8+uijvPfee1xxxRWcc845vPPOOw2+t7a2tv61xx57tPmxZFvXrl1Zs2YNp5xyCm+88Qbr168HoLS0tNXfcc8993DIIYew//77b/W9EyZMYPr06axevbrBdx511FF07tyZPn36cO211/LQQw9l76AkSe2WbXbTstFmN+eRRx7hxz/+MQ8++CA777xzg7LFixdz2mmnceedd241vl+S1DQT9CbssMMOXHPNNcyZM6d+xtdNzj77bCZMmNDkviUlJXzjG9+gS5cuvPHGG20caW7tv//+vPzyy/XLL774Ij179qyfeKc17rjjjgZ3z5v63oEDBza6f6dO/rWVJP0f2+zGZaPNbsqTTz7J2LFjefjhh7cav79ixQq++MUv8v3vf5/hw4d/7LokqZiY6TRjxx135OKLL+Z73/tei9s+//zz3H777dTW1rJx40buuece1q9fz4ABA3IQae4MGzaMlStXcvfdd/P+++8zYcKErSZyq66uZr/99mPGjBlb7b9gwQKef/55zj777AbrzzjjDGbMmMGCBQuorq5m0qRJnHXWWQDMnz+f5557jg0bNvDOO+9wzTXX8JWvfKXtDlKS1O7YZm/t47TZKSXWrVvH+vXr6z9/+OGHALzwwgucccYZTJs2jcGDBzfYb+XKlQwfPpx//ud/5rzzzmvbA5SkDsgEvQVjxozhmWeeYc6cOc1uV1paytSpU9l3330pKyvjv/7rv5gxYwY9evSo3+YrX/lKg2eqTp48ua3Dz7pddtmFe++9l2uuuYZevXqxceNGqqqqGmyzfv16XnvtNVatWrXV/nfccQcjRoygZ8+eDdYPGTKE8ePHc+yxxzJ48GDOPIzzLvEAACAASURBVPPM+gS9pqaG0047je7duzNo0CB69+7NDTfc0HYHKUlql2yzG/o4bfabb75JSUkJJ5xwAkuWLKGkpITjjz8egOuvv563336b448/vv78fOlLXwLqhgY899xzfPvb325w/iRJreNj1rawePHiBstlZWWsXr0aoNFZXpctW1b/+de//nWrv7c9O+aYY3jllVeaLK+oqGhyRtzmuhlecsklXHLJJVutHzZsGK+//vq2BypJ6tBss1u2vW12c235HXfcwR133NFo2ciRIxk5cuT2BaucKy8rafFRqOVlJTwxbliOIpJkgi5JkiQVodYk3i0l8JKyyy7ukiRJkiQVABN0SZIkSZIKgAm6JEmSJEkFwARdkiRJkqQCYIIuSZIkSVIBMEGXJEmSJKkAZPUxaxGxGFgNbAA+SilVRkRPYCpQASwGzkwpvZfNeiVJkiRJau/a4g76sSmlQ1NKlZnlccCslNIAYFZmWZIkSZIkbSYXXdxPAiZnPk8GTs5BnZIkSZIktSvZTtAT8JuImBcRYzPreqeUlgNk3vfIcp2SJEmSJLV7WR2DDgxNKb0VEXsAj0bEy63ZKZPMjwXo169flkOSJEmSJKnwZfUOekrprcz7CuA+4Ajg7YjoA5B5X9HIfhNTSpUppcpevXplMyRJkiRJktqFrCXoEdE1Irpt+gwcD/wFeBAYmdlsJPBAtuqUJEmSJKmjyGYX997AfRGx6XvvTin9b0Q8A0yLiDHAEuCMLNYpSZIkSVKHkLUEPaX0OnBII+vfAY7LVj2SJElSsSsvK6Fi3MwWt5HUvmR7kjhJkiRJbeyJccPyHYKkNpCL56BLkiRJkqQWmKBLkiRJklQATNAlSSpCEVEWEdMj4uWIeCkijoqInhHxaEQsyrzvmu84JUkqJo5Bl9TmmpvIxglspLy5AfjflNLpEbETsAtwNTArpVQVEeOAccCV+QxSkqRiYoIuqc05kY1UWCKiO3A0MAogpfQh8GFEnAQck9lsMvA4JuiSJOWMXdwlSSo++wA1wB0R8VxE/CwiugK9U0rLATLve+QzSEmSio0JuiRJxacz8CnglpTSYcD71HVnb5WIGBsRcyNibk1NTVvFKElS0TFBlySp+CwDlqWUnsosT6cuYX87IvoAZN5XNLZzSmliSqkypVTZq1evnAQsSVIxMEGXJKnIpJT+CiyNiIGZVccBLwIPAiMz60YCD+QhPEmSipaTxEmSVJwuBn6RmcH9dWA0dRfup0XEGGAJcEYe45MkqeiYoEuSVIRSSn8GKhspOi7XsUiSpDp2cZckSZIkqQCYoEuSJEmSVADs4i5JH9PQqtlU165ttKy8rIQnxg3LcUSSpHxqrl3YxPZBUmNM0CXpY6quXcviqhGNllWMm5njaCRJ+dZcu7CJ7YOkxtjFXZIkSZKkAmCCLkmSJElSATBBlyRJkiSpADgGXVJRKS8raXLcX3lZSY6jkSRJkv6PCbqkouKMuZIkSSpUdnGXJEmSJKkAmKBLkiRJklQATNAlSZIkSSoAJuiSJEmSJBUAE3RJkiRJkgqACbokSZIkSQXABF2SJEmSpAJggi5JkiRJUgEwQZckSZIkqQB0zncAkiRJUrEpLyuhYtzMZsslFR8TdEmSJCnHnhg3LN8hSCpAdnGXJEmSJKkAeAddklqhua6IdkOUJHVUremKb28AKXtM0CWpFfzxIUkqRi21f80l75K2nV3cJUmSJEkqACbokiRJkiQVABN0SZIkSZIKgAm6JEmSJEkFwARdkiRJkqQCkNUEPSJ2iIjnIuLhzHL/iHgqIhZFxNSI2Cmb9UmSJEmS1FFk+zFrlwAvAd0zy/8JXJ9SmhIRtwJjgFuyXKckSZKkdmpo1Wyqa9c2We6z1lVMspagR8RewAjgB8ClERHAMODczCaTgWswQZckKe8iYjGwGtgAfJRSqoyInsBUoAJYDJyZUnovXzFKKg7VtWtZXDWiyXKfta5iks0u7j8BrgA2ZpZ3A2pTSh9llpcB5Y3tGBFjI2JuRMytqanJYkiSJKkZx6aUDk0pVWaWxwGzUkoDgFmZZUmSlCNZSdAj4svAipTSvM1XN7Jpamz/lNLElFJlSqmyV69e2QhJkiRtu5Oo6/FG5v3kPMYiSVLRyVYX96HAiRFxAtCFujHoPwHKIqJz5i76XsBbWapPkiR9PAn4TUQk4LaU0kSgd0ppOUBKaXlE7NHYjhExFhgL0K9fv1zFK0lSh5eVO+gppatSSnullCqAs4HZKaWvAo8Bp2c2Gwk8kI36JEnSxzY0pfQp4EvAv0TE0a3d0Z5vkiS1jbZ+DvqV1E0Y9yp1Y9IntXF9kiSpFVJKb2XeVwD3AUcAb0dEH4DM+4r8RShJUvHJeoKeUno8pfTlzOfXU0pHpJT2SymdkVL6INv1SZKkbRMRXSOi26bPwPHAX4AHqevxBvZ8kyQp57L9HHRJklT4egP31T0Rlc7A3Sml/42IZ4BpETEGWAKckccYJUkqOibokiQVmZTS68Ahjax/Bzgu9xFJkiRo+zHokiRJkiSpFUzQJUmSJEkqAHZxl9RqQ6tmU127ttGy8rKSHEcjSVL2NdfWbVJeVsIT44blKCJJxcQEXVKrVdeuZXHViHyHIUlSm2lNW1cxbmaOopFUbOziLkmSJElSATBBlyRJkiSpAJigS5IkSZJUAEzQJUmSJEkqACbokiRJkiQVAGdxlyRJkrRdystKWpzV3sfSSa1ngi5JkiRpu7Qm8faxdFLr2cVdkiRJkqQC4B10SWpDzXX9a67L39Cq2VTXrt3m/SRJktR+maBLUhtqLpFurstfde1aFleN2Ob9JEnNa+4CKNRdBJWkfDFBlyRJUtFo7gKoJOWbY9AlSZIkSSoAJuiSJEmSJBUAu7hLkiRJ26ClZ387jl3S9jJBlyRJkraBT9KQ1Fbs4i5JkiRJUgEwQZckSZIkqQDYxV2S2pnmxj6Wl5XY9VKSJKmdMkGXpHamuQS8uUmLJEmSVNjs4i5JkiRJUgEwQZckSZIkqQDYxV2SJElSm/G58VLrmaBL7djQqtlU165ttGx7Jwtr6TslSZK2hZOXSq1ngi61Y9W1a1lcNaLRsu2dLKy575QkSZLUdhyDLkmSJElSATBBlySpSEXEDhHxXEQ8nFnuHxFPRcSiiJgaETvlO0ZJkoqJCbokScXrEuClzZb/E7g+pTQAeA8Yk5eoJEkqUibokiQVoYjYCxgB/CyzHMAwYHpmk8nAyfmJTpKk4mSCLklScfoJcAWwMbO8G1CbUvoos7wMKG9sx4gYGxFzI2JuTU1N20cqSVKRMEGXJKnIRMSXgRUppXmbr25k09TY/imliSmlypRSZa9evdokRkmSipGPWZMkqfgMBU6MiBOALkB36u6ol0VE58xd9L2At/IYoyRJRcc76JIkFZmU0lUppb1SShXA2cDslNJXgceA0zObjQQeyFOIkiQVJRN0SZK0yZXApRHxKnVj0iflOR5JkopK1rq4R0QX4PfAzpnvnZ5SGh8R/YEpQE/gWeD8lNKH2apXkiRtv5TS48Djmc+vA0fkMx5JkopZNu+gfwAMSykdAhwKDI+II/GZqpIkSZIktShrCXqqsyazuGPmlfCZqpIkSZIktSirY9AjYoeI+DOwAngUeI1WPlNVkiRJkqRiltUEPaW0IaV0KHWPZjkC+GRjm225IiLGRsTciJhbU1OTzZAkSZIkSWoX2mQW95RSLXUTzhxJ5pmqmaJGn6maUpqYUqpMKVX26tWrLUKSJEmSJKmgZS1Bj4heEVGW+VwCfB54CZ+pKkmSJElSi7L2mDWgDzA5InagLvGfllJ6OCJeBKZExPeB5/CZqpIkSZIkbSVrCXpK6QXgsEbW+0xVSZIkSZJa0CZj0CVJkiRJ0rYxQZckSZIkqQCYoEuSJEmSVABM0CVJkiRJKgAm6JIkSZIkFQATdEmSJEmSCkA2n4MuSZIkSTk3tGo21bVrmywvLyvhiXHDchiRtH1M0CVJkiS1a9W1a1lcNaLJ8opxM3MYjbT97OIuSZIkSVIB8A66VISa6wZWXlaS42iUTeVlJU3eJbB7nyRJUmEzQZeKUEvdwNR+NZeA271PkiSpsNnFXZIkSZKkAmCCLkmSJElSATBBlyRJkiSpADgGXZLypKUJ3SRJklRcTNAlKU+cUV2SJEmbs4u7JEmSJEkFwARdkiRJkqQCYIIuSZIkSVIBcAy6JBWJlialc0y8JElSfpmgS1KRaC4BbypxlyRJUu7YxV2SJEmSpAJggi5JUpGJiC4R8XREPB8RCyLiu5n1/SPiqYhYFBFTI2KnfMcqSVIxMUGXJKn4fAAMSykdAhwKDI+II4H/BK5PKQ0A3gPG5DFGSZKKjmPQpQ6qpQnBJBWvlFIC1mQWd8y8EjAMODezfjJwDXBLruOTJKlYmaBLHZQzcktqTkTsAMwD9gNuBl4DalNKH2U2WQaUN7HvWGAsQL9+/do+WGkbDK2aTXXt2ibLvUgtqZCZoEuSVIRSShuAQyOiDLgP+GRjmzWx70RgIkBlZWWj20j5Ul27lsVVI/IdhiRtF8egS5JUxFJKtcDjwJFAWURsuni/F/BWvuKSJKkYmaBLklRkIqJX5s45EVECfB54CXgMOD2z2UjggfxEKElScbKLuyRJxacPMDkzDr0TMC2l9HBEvAhMiYjvA88Bk/IZpCRJxcYEXZKkIpNSegE4rJH1rwNH5D4iSZIEdnGXJEmSJKkgmKBLkiRJklQA7OIuSaK8rISKcTObLX9i3LAcRiRJklR8TNAlSS0m380l75IkScoOu7hLkiRJklQATNAlSZIkSSoAJuiSJEmSJBUAE3RJkiRJkgqACbokSZIkSQUgawl6RPSNiMci4qWIWBARl2TW94yIRyNiUeZ912zVKUmSJElSR5HNO+gfAZellD4JHAn8S0QMAsYBs1JKA4BZmWVJkiRJkrSZrCXoKaXlKaVnM59XAy8B5cBJwOTMZpOBk7NVpyRJkiRJHUWbjEGPiArgMOApoHdKaTnUJfHAHo1sPzYi5kbE3JqamrYISZIkSZKkgtY5218YEaXAL4F/TSmtiogW90kpTQQmAlRWVqZsxyRJkqT2b2jVbKpr1za7TXlZSY6iUa6Ul5VQMW5mi9tk4zueGDdsm+OTsimrCXpE7Ehdcv6LlNKMzOq3I6JPSml5RPQBVmSzTkmSJBWH6tq1LK4ake8wlGPZSJpb8x0tJfBSLmRzFvcAJgEvpZSu26zoQWBk5vNI4IFs1SlJkiRJUkeRzTvoQ4HzgfkR8efMuquBKmBaRIwBlgBnZLFOSZIkSZI6hKwl6CmlPwJNDTg/Llv1SJIkSZLUEbXJLO6SJEmSJGnbmKBLkiRJklQATNAlSZIkSSoAJuiSJEmSJBUAE3RJkiRJkgqACbokSZIkSQXABF2SJEmSpAJggi5JkiRJUgEwQZckSZIkqQB0zncAkmBo1Wyqa9c2WlZeVsIT44blOCJJkiRJuWaCLhWA6tq1LK4a0WhZxbiZOY5GkqT8aO6CNdRdtJakjswEXZKkIhMRfYE7gU8AG4GJKaUbIqInMBWoABYDZ6aU3stXnCo+zV2wlqRi4Bh0SZKKz0fAZSmlTwJHAv8SEYOAccCslNIAYFZmWZIk5YgJuiRJRSaltDyl9Gzm82rgJaAcOAmYnNlsMnByfiKUJKk4maBLklTEIqICOAx4CuidUloOdUk8sEcT+4yNiLkRMbempiZXoUqS1OGZoEuSVKQiohT4JfCvKaVVrd0vpTQxpVSZUqrs1atX2wUoSVKRMUGXJKkIRcSO1CXnv0gpzcisfjsi+mTK+wAr8hWfJEnFyARdkqQiExEBTAJeSildt1nRg8DIzOeRwAO5jk2SpGLmY9YkSSo+Q4HzgfkR8efMuquBKmBaRIwBlgBn5Ck+SWqXhlbNprp2bZPl5WUlPDFuWA4jUntjgi5JUpFJKf0RiCaKj8tlLJLUkVTXrmVx1YgmyyvGzcxhNGqP7OIuSZIkSVIBMEGXJEmSJKkA2MVdkiRJba6lsblQNz5XypfyspJmu6A7fly5YIIuSZKkNtfS2Fwp31pKvh0/rlywi7skSZIkSQXABF2SJEmSpAJggi5JkiRJUgEwQZckSZIkqQA4SZxU4JqbUdTZbiVJkqSOwwRdKnA+zkOSJEkqDnZxlyRJkiSpAJigS5IkSZJUAEzQJUmSJEkqAI5BlyR9LEOrZlNdu7bRsvKyEudRkCRJaiUTdEnSx1Jdu5bFVSMaLWvqCQSSJEnaml3cJUmSJEkqACbokiRJkiQVALu4S5IkSVILystKWhy6VV5WkqNo1FGZoEuSJElSC5z0VLmQtS7uEXF7RKyIiL9stq5nRDwaEYsy77tmqz5JkiRJkjqSbI5B/zkwfIt144BZKaUBwKzMsiRJkiRJ2kLWEvSU0u+Bd7dYfRIwOfN5MnBytuqTJEmSJKkjaetZ3HunlJYDZN73aOP6JEmSJElqlwpikriIGAuMBejXr1+eo5GaN7RqNtW1axstKy8rcQIRdUjNzVzrjLWSJEnZ0dYJ+tsR0SeltDwi+gArGtsopTQRmAhQWVmZ2jgm6WOprl3L4qoRjZa19OgNqb3ywpMkSVLba+su7g8CIzOfRwIPtHF9kiRJkiS1S9l8zNo9wBxgYEQsi4gxQBXwhYhYBHwhsyxJkiRJkraQtS7uKaVzmig6Llt1SJIkSZLUURXEJHFSMWhpcjmpI2ppcjnHtkuSJP0fE3QpR5qbXE7qqJpLwJ1UUZIkqaG2niROkiQVoIi4PSJWRMRfNlvXMyIejYhFmfdd8xmjJEnFxgRdkqTi9HNg+BbrxgGzUkoDgFmZZUmSlCMm6JIkFaGU0u+Bd7dYfRIwOfN5MnByToOSJKnImaBLkqRNeqeUlgNk3vdobKOIGBsRcyNibk1NTU4DlCSpIzNBlyRJ2ySlNDGlVJlSquzVq1e+w5EkqcMwQZckSZu8HRF9ADLvK/IcjyRJRcUEXZIkbfIgMDLzeSTwQB5jkSSp6PgcdEmSilBE3AMcA+weEcuA8UAVMC0ixgBLgDPyF6EKydCq2VTXrm12m/KyEp4YNyxHEUlSx2SCLklSEUopndNE0XE5DUTtQnXtWhZXjWh2m4pxM3MUjSR1XHZxlyRJkiSpAJigS5IkSZJUAEzQJUmSJEkqACbokiRJkiQVABN0SZIkSZIKgLO4S5IkSVIOlJeVNPvEAx9XKBN0SZIkScqBlpJvH1cou7hLkiRJklQAvIMuZVFz3ZbKy0pyHI1U2Fr692IXP0mSVGxM0KUsMqGQWq+5fy928ZOya2jVbKpr1zZZno2LYq0ZWytJap4JuiRJUgdXXbuWxVUjmizPxkUxL1JL0sfnGHRJkiRJkgqACbokSZIkSQXALu5qF5obO9cWk0m1VJ8kSYWipfHl0HLb5fhxqTC09G9x0zYOKem4TNDVLjQ3dq4tJpNqaayeJEmFIhttlj/2pcLQmn+LTqTasdnFXZIkSZKkAmCCLkmSJElSAbCLuyRJkiR1EK2dl8KhLYXJBF1Fy4ngpPYp15NGSpLUnrRmXgrHsRcuE3QVLSeCk9qnXE8aKUmSlCuOQZckSZIkqQCYoEuSJEmSVADs4q4OzXHmUvtUXlbSZHf15v7ttrSf49MlSVIhM0FXh+Y4c6l92t5Eurn9HJ8uSZIKnV3cJUmSJEkqAN5BlyRJkqQi0tyQsE3lDgvLDxN0SZIkSSoiLSXfDgvLHxN0ZV1LE7Nl+2rc9k4mJam4OIGcJEkqdCboyrrmJmZri6tx/qiW1BpOIKdsa+6C9CbZuPjTUj1ejJaKS2u6p7d1HZu2ae7/t9b83+Xv+K3lJEGPiOHADcAOwM9SSlW5qFeSJG072+3Wac2TQrJx8ccnkkjaXC6S2tbU0dL/by393+XF8ca1+SzuEbEDcDPwJWAQcE5EDGrreiVJ0raz3ZYkKX9y8Zi1I4BXU0qvp5Q+BKYAJ+WgXkmStO1styVJypNIKbVtBRGnA8NTShdmls8HhqSULtpsm7HA2MziQGBhG4SyO/C3Nvje9s7z0jjPS+M8L43zvDTNc9O4tjgve6eUen3cLymQdtu/N9nnOc0uz2d2eT6zy/OZfTlrt3MxBj0aWdfgqkBKaSIwsU2DiJibUqpsyzraI89L4zwvjfO8NM7z0jTPTeMK/Lzkvd0u8PPTLnlOs8vzmV2ez+zyfGZfLs9pLrq4LwP6bra8F/BWDuqVJEnbznZbkqQ8yUWC/gwwICL6R8ROwNnAgzmoV5IkbTvbbUmS8qTNu7inlD6KiIuAX1P3uJbbU0oL2rreRrRpF/p2zPPSOM9L4zwvjfO8NM1z07iCPS8F0m4X7Plpxzyn2eX5zC7PZ3Z5PrMvZ+e0zSeJkyRJkiRJLctFF3dJkiRJktQCE3RJkiRJkgpAh0/QI2J4RCyMiFcjYly+4ykUEXF7RKyIiL/kO5ZCERF9I+KxiHgpIhZExCX5jqlQRESXiHg6Ip7PnJvv5jumQhIRO0TEcxHxcL5jKRQRsTgi5kfEnyNibr7jKRQRURYR0yPi5cz/NUflO6ZCEBE9I+LRiFiUed+1kW0OjYg5mf+DXoiIs/IRayFr6TdPROwcEVMz5U9FREXuo2w/WnE+L42IFzN/H2dFxN75iLM9ae3v8og4PSJSRPiosGa05nxGxJmZv6f/n707j7Ox7v84/vpYZhoGgyhhGiRlXyZRt6VI3ZGQQmSJJMuv5U4ohLtbo52s3fZSlCjRnVa6qylrlEqIsoXsS8LM9/fHOeaeYXCGM3OdM/N+Ph7nMedcy/f6XNc5c67zua7vssbMXs/qGMNJAP/zsf5cYaX///7WTIkjO7dBN7PcwM/ATfiGjVkKtHPO/eBpYCHAzOoDh4DpzrnKXscTCsysBFDCObfCzAoAy4EW+ryAmRmQ3zl3yMzyAl8ADzrnvvY4tJBgZo8A8UBB51wzr+MJBWa2CYh3zv3hdSyhxMymAf91zk3095Cezzm3z+u4vGZmzwB7nHMJ/h9FhZ1z/U5Z5krAOefWmdll+L6jr9bx8wnkN4+Z9QSqOud6mFlboKVzThc60hHg8bwB+MY5d8TMHgAa6nieWaC/y/2/wRYAEUBv55wu8qYjwM9oeeBN4Ebn3F4zK+6c2+lJwCEuwOP5CrDSOTfOzCoC7zvn4oIdS3a/g14bWO+c+8U5dwyYCdzucUwhwTn3ObDH6zhCiXNuu3Nuhf/5QeBHoKS3UYUG53PI/zKv/5F9r+5lgJmVApoCE72ORUKbmRUE6gOTAJxzx5RcprgdmOZ/Pg1oceoCzrmfnXPr/M+3ATuBYlkWYegL5DdP6uM8G2jkvwArpzvn8XTOfeacO+J/+TVQKotjDDeB/i7/J/AMcDQrgwtDgRzP+4Axzrm9AErOzyqQ4+mAgv7nhYBtmRFIdk/QSwKbU73eghIuCYC/2l8N4BtvIwkd/mrc3+L7UfyRc07Hxucl4DEg2etAQowDPjSz5WbW3etgQkRZYBcwxV89bqKZ5fc6qBBxiXNuO/gulgLFz7awmdXGd3dtQxbEFi4C+c2Tsoxz7gSwHyiaJdGFn4z+huwK/CdTIwp/5zymZlYDKO2cU5OxcwvkM3olcKWZfWlmX5vZLVkWXfgJ5HgOATqY2RbgfaBPZgSS3RP09K4K666fnJWZRQNvAw855w54HU+ocM4lOeeq47tDUNvMcnzTCDNrBux0zi33OpYQdL1zribwd6CXv1lNTpcHqAmMc87VAA4DOaZvFDP72My+T+eRoZpt/uZIrwJdnHO6MPY/gfzm0e+iwAV8rMysA75mTs9makTh76zH1MxyAS8C/8iyiMJbIJ/RPEB5oCHQDphoZjGZHFe4CuR4tgOmOudKAbcCr/o/t0GVJ9gFhpgtQOlUr0uRSVURJHvwt69+G5jhnJvjdTyhyDm3z8wWAbcAOb2TweuB5v5OQi4CCprZa865Dh7H5Tl/FWScczvNbC6+qmOfexuV57YAW1LVPplNDkrQnXONzzTPzHaYWQnn3HZ/Ap5uNUx/M4EFwED1gXGaQH7znFxmi5nlwVdFU83d0hfQb0gzaww8ATRwzv2VRbGFq3Md0wJAZWCRv+XFpcA8M2uudujpCvR//mvn3HFgo5mtxZewL82aEMNKIMezK77fvzjnEs3sIuBiznDOOl/Z/Q76UqC8mZXxd8bTFpjncUwSovzt8CYBPzrnXvA6nlBiZsVOXnE1syigMfCTt1F5zzk3wDlXyt9BSFvgUyXnYGb5/Z384K/C3QRdzME59zuw2cwq+Cc1AnJ8J5R+84BO/uedgHdPXcB/Hp+Lr3PTt7IwtnARyG+e1Me5Nb7vLN1BT985j6e/OvYEoLna9gbkrMfUObffOXexcy7Of179Gt+xVXKevkD+598BbgAws4vxVXn/JUujDB+BHM/f8J27MbOr8d2c2RXsQLJ1gu5vX9UbWIivw683nXNrvI0qNJjZG0AiUMHMtphZV69jCgHXA/cAN5pvaKhvM2v4hDBUAvjMzFbj+wL7SO3D5CwuAb4ws1XAEmCBc+4Dj2MKFX2AGf7/perAcI/jCRUJwE1mtg5fD7oJAGYWb2YnO2C8C18ne51TfUdX9ybc0HOm3zxmNszMmvsXmwQUNbP1wCPkoBocGRXg8XwWiAbe8n8edRPoLAI8phKgAI/nQmC3mf0AfAb0dc7t9ibi0Bbg8fwHcJ//980bQOfMuMiZrYdZExEREREREQkX2foOuoiIiIiIiEi4UIIuIiIiIiIiEgKUoIuIiIiIiIiEACXoIiIiIiIiIiFACbqIiIiIiIhICFCCLiIiIiIiIhIClKCLiIiIiIiIhAAl6CIiIiIiIiIhQAm6iIiIiIiISAhQgi4iIiIiIiISApSgi4iIiIiIiIQAJegiIiIiIiIiIUAJuoiHpsYJ1QAAIABJREFUzGyNmTX0Og4vmVlLM9tsZofMrMZ5rD/VzJ46y/xDZlb2wqIUEZGcTufsCz9nZ3BbZz2/i2RXStBFMomZbTKzxqdM62xmX5x87Zyr5JxbdI5y4szMmVmeTArVa88BvZ1z0c65lafONJ//M7PvzeywmW0xs7fMrEoghfvL/cVfVl9/OQfNbKOZ9Q3yvoiISBjSOTtg5zpnO/+5+pCZbTWzF8wstwdxioQtJegiOVwI/Ii4HFhzlvkjgQeB/wOKAFcC7wBNz2NbBnQECgO3AL3NrO15lCMiIpLlwuCcDVDNORcNNALuBu47dYEQ2A+RkKUEXcRDqa/Ym1ltM1tmZgfMbIeZveBf7HP/333+K9J1zSyXmQ00s1/NbKeZTTezQqnK7eift9vMBp2ynSFmNtvMXjOzA0Bn/7YTzWyfmW03s9FmFpGqPGdmPc1snf/u8z/NrJx/nQNm9mbq5U/Zx3RjNbNIMzsE5AZWmdmGdNYtD/QC2jnnPnXO/eWcO+Kcm+GcS0i1aGEzW+CP7RszK3dK7FcAOOeecc6tcM6dcM6tBd4Frs/o+yYiIjmPztlnP2efyjn3E/BfoHKq49fPzFYDh80sj5ldbWaL/Puyxsyan1LMxWb2kX8/FpvZ5QG8VSJhTQm6SOgYCYx0zhUEygFv+qfX9/+N8VcpSwQ6+x83AGWBaGA0gJlVBMYC7YESQCGg5Cnbuh2YDcQAM4Ak4GHgYqAuvqvePU9Z5xagFlAHeAx4xb+N0vhOvu3OsF/pxupPtqP9y1RzzpVLZ91GwBbn3JIzlH1SO2Aovjvj64F/nWN5zMyAepz7ToCIiMipdM4+B/++1QNSV4Vvh68GXAy+Wm3vAR8CxYE+wAwzq5Bq+fbAP/37+q1//0WyNSXoIpnrHf9V4X1mtg/fSfhMjgNXmNnFzrlDzrmvz7Jse+AF59wvzrlDwACgrfmqjLUG3nPOfeGcOwYMBtwp6yc6595xziU75/50zi13zn3tv7O8CZgANDhlnRHOuQPOuTXA98CH/u3vB/4DnKmzmLPFei5Fge0BLDfHObfEOXcC38m7egDrDMH3HTglgGVFRCT70zn7ws7ZJ60ws734ku+JpD3PjnLObXbO/Ynv4kE0kOCcO+ac+xSYT9qLBwucc5875/4CngDqmlnpDMQiEnaUoItkrhbOuZiTD06/wp1aV3ztq38ys6Vm1uwsy14G/Jrq9a9AHuAS/7zNJ2c4544Au09Zf3PqF2Z2pZnNN7Pf/VXohuO7Wp3ajlTP/0zndTTpO1us57Ib3x2Fc/k91fMjZ4kFADPrja8telP/SV9ERETn7As7Z59U0zlX2DlXzjk30DmXfIZ9uQzYfMr8X0lbgyD1sTkE7PGvJ5JtKUEXCRHOuXXOuXb4qnmNAGabWX5Ov5IOsA1fRy0nxQIn8J2AtwOlTs4wsyh8d6LTbO6U1+OAn4Dy/up6j+OrehYMZ4v1XD4BSplZfJBiwczuBfoDjZxzW4JVroiI5Bw6Z5+31PuyDShtZqnzkVhga6rXKXfLzSwaX2ex24IUi0hIUoIuEiLMrIOZFfNfSd7nn5wE7AKS8bUFO+kN4GEzK+M/YQ0HZvmreM8GbjOz6/ydwAzl3CfuAsAB4JCZXQU8ELQdO3usZ+WcW4eviuEbZtbQzCLM7CIza2tm/TMaiJm192//ppNDr4mIiGSUztlB8Q1wGHjMzPKab4z524CZqZa51cz+5j82/wS+cc5tPr0okexDCbpI6LgFWOPvJXUk0NY5d9Rf3e1fwJf+dnF1gMnAq/h6i90IHMXXuQr+9mZ98J3gtgMHgZ3A2apyP4pvKJSDwL+BWUHcrzPGGqD/w9eZzhh8P4I2AC3xtW3LqKfw3ZlYar7edQ+Z2fjzKEdERHI2nbMvkL/NfXPg78Af+C7Id/T3/n7S68CT+Kq218LXRl4kWzPn0quJIyLZhf8K+D58VeE2eh2PiIiIpE/nbBHRHXSRbMjMbjOzfP72cM8B3wGbvI1KRERETqVztoikpgRdJHu6HV8nKtuA8viq3qm6jIiISOjROVtEUqiKu4iIiIiIiEgI0B10ERERERERkRCQx+sATnXxxRe7uLg4r8MQEREJK8uXL//DOVcsq7er87aIiEjGnem8HXIJelxcHMuWLfM6DBERkbBiZr96sV2dt0VERDLuTOdtVXEXERERERERCQFK0EVERERERERCgBJ0ERERERERkRCgBF1EREREREQkBChBFxEREREREQkBStBFREREREREQoASdBEREREREZEQoARdREREREREJAQoQRcREREREREJAUrQRUREREREREKAEnQRERERERGREKAEXURERERERCQE5PE6AAkDe/fCrFmwcCGsWQO7dkHevBAXB3XqQLNm0KgR5M7tdaQiIiIiIiJhSwm6nNmBA/DUUzBmDBw5AmXKQHw8XHIJHDsG69fDv/8NL78M5crBY4/BvfdCHn2sREREJPu6PuFTtu77M6hlloyJ4sv+Nwa1TBEJP8qkJH1ffAFt2sD27dChAzzyCFSrBmZplztyBObPh+efh/vvh/HjfUl7rVrexC2STfXo0YOSJUsyaNAgr0MREcnxtu77k00JTYNaZlz/BUEtT0TCU0Bt0M3sFjNba2brzax/OvPrm9kKMzthZq1TTa9uZolmtsbMVptZm2AGL5lk7Fho2BDy5YOvv4bp06F69dOTc/Atc9ddvuXefBN27IC6dWHUKHAuy0MXyUydO3cmIiKC6OjolMeTTz4Z8PoNGzZk4sSJ57Xt8ePHh0Vy/tRTT1GpUiVy5crF1KlT08zbsmULDRs2JF++fNSsWZPvv/8+Zd7x48fp2rUrBQoUIDY2ljfffDPNuqNGjeKSSy6hSJEiDBgwICt2RURERCTLnTNBN7PcwBjg70BFoJ2ZVTxlsd+AzsDrp0w/AnR0zlUCbgFeMrOYCw1aMtFzz0GvXvD3v8OyZVC7dmDrmcGdd8J338Ett8CDD8J998GJE5kbr0gWe+yxxzh06FDKY+jQoV6HFFLKli3LqFGjqJVOLZru3btTuXJldu/eTZs2bWjT5n/XbF988UW+++47Nm/ezPTp07n33nvZvHkzAEuWLGHIkCF8+umnfPfdd8ycOfO0BF5EREQkOwjkDnptYL1z7hfn3DFgJnB76gWcc5ucc6uB5FOm/+ycW+d/vg3YCRQLSuQSfGPGQN++vqrtc+ZAoUIZL6NIEXj3XRg0CCZNgjvugD+D20ZLJBRNnTqV6667js6dO1OgQAH+9re/8ccffwAwfPhwoqOj+e9//0vv3r2Jjo4+LYGNi4tj9OjR1KpVi/z589O8eXMA5s+fT3R0NHnz5mXgwIGnbXfOnDlUrlyZwoULc+utt7J9+/aUeQsWLKBChQoUKFCA8uXLs3Dhwkw8Aj533303jRo1IjIyMs30AwcO8OGHH9K/f3+ioqJ4+OGH+fXXX1m9ejUAb731Fv/3f/9HTEwMDRs2pG7dusydOzdlXqtWrahUqRIlS5akW7duzJw5M9P3RURERCSrBZKglwQ2p3q9xT8tQ8ysNhABbEhnXnczW2Zmy3bt2pXRoiUYPvgA/u//4Lbb4LXXfL20ny8zGDbMV1X+vfegVStfp3Ii2dySJUu4++672blzJ8ePH2fy5MkAPP744xw6dIh69eoxevRoDh06xPLly09bf8KECUyfPp39+/fz+OOPA9CsWTMOHTpE+/bt091e165dmTRpErt27aJGjRp07949ZX63bt0YNmwYBw8e5MMPP6RkyQx/dQfN+vXrueiii4iOjqZevXps2bKFcuXKsXbtWgDWrl3LVVddRYcOHZgzZw4VK1Y8bd7IkSMZMGBAmnkiIiIi2UkgCXo6DY/JUONiMysBvAp0cc4lnzrfOfeKcy7eORdfrJhusGe5jRt9d82rVIHXXw9eL+wPPACvvOJL/tu3V3V3yRaee+45YmJiUh5z5sxJmXfllVfSpEkToqKiaNSoET///HOGyu7evTuVKlUiT5481KlT55zLT5o0ic6dO3PttdeSJ08eHn30URYsWMBff/0FQK5cudiwYQMHDhygTJkyVK5cOWM7G0SHDx8mOjqagwcPsm7dOvbu3UuBAgU4dOhQmvlr165l69at6c7bvHkz69atSzNPREREJDsJJEHfApRO9boUsC3QDZhZQWABMNA593XGwpNMd+KEL3kGeOcdiI4Obvnduvl6eJ8929cTvEiYe/TRR9m3b1/Ko1WrVinzihQpkvI8IiKCo0ePZqjs8uXLZ2j5zZs3M2HChJSLBWXKlCEiIiKlmvtbb71FYmIisbGx1KlTJ02nbFktf/78HDp0iNKlS/P7779Tq1YtDh48SLT/O+fk/KVLl9KnT5905z333HPMnj07zTwRERGR7CSQBH0pUN7MyphZBNAWmBdI4f7l5wLTnXNvnX+YkmmGD4fERN/waHFxmbONRx6Bhx/2jZd+nj1Yi2QHuXKd/Ss3TwZrr5QuXZpBgwaluWBw9OhR4vz/y9dddx3vvfceO3bs4Kqrrkq3DXtWueKKKzh69ChbtmwB4NixY2zYsIEKFSoAvtoHP/30U8ryP/zwQ0DzRERERLKTcybozrkTQG9gIfAj8KZzbo2ZDTOz5gBmdo2ZbQHuBCaY2Rr/6ncB9YHOZvat/1E9U/ZEMm7tWnjqKWjbFtq1y9xtPfMMNGkCPXv6xlgXyYEuvfTSoN7F7tKlC+PHj2fFihU459i5cyezZs0CIDk5mddee42DBw+mXBgoWLBg0LZ9JsePH+fo0aM451KeJycnU7BgQZo0aUJCQgJHjx7lxRdfJDY2lipVqgBw1113MWrUKPbt28fixYtJTEykZcuWANx5553MmTOHNWvWsHXrViZNmpSmB3gRERGR7CKgcdCdc+875650zpVzzv3LP22wc26e//lS51wp51x+51xR/7BqOOdec87ldc5VT/X4NvN2RwLmnC9ZzpcPXnop87eXJw/MnOm7S3/nnaDOACVMPfPMM2nGQe/WrVvA6/7jH//go48+okSJEjRs2DCgdZo0aUJ0dDQzZsxI2Xbnzp0BqFOnDi+88AJdunShUKFCxMfHp+l87rXXXqN06dJcfPHFbNu2jX/9618Z2dXzct999xEVFcVXX31F9+7diYqK4vPPPwd8neB99913FC5cmJkzZzJr1izMfN2cPPzww1SuXJnSpUvToUMHJk2aROnSvtZV1157LU8++SQ33HADlStX5q677lKCLiIiItmSOZeh/t4yXXx8vFu2bJnXYWR/s2b57pyPGwc9emTddlev9o2t3qiRr4f3c1T5FRGRwJjZcudcfFZvV+dtyYni+i9gU0LTkC9TRELXmc7byo5yomPH4PHHoVo1uO++rN121aq+TuPefx9GjcrabYuIiIiIiIQwJeg50cSJ8Msv8PTTkDt31m+/Z09o3hweewy++y7rty8iIiIiIhKClKDnNIcPw7BhUL8+3HKLNzGYwaRJEBMD996r8dFFRERERERQgp7zjB0LO3b47p77O2fyxMUXw+jRsGwZvPiid3GIiIiIiIiECCXoOclff/mS4caN4brrvI7G15t7y5YweDD8/LPX0YiIiIiIiHhKCXpO8uqrsH079OvndSQ+ZjBmDFx0ka+zuhAbUUAktUWLFmFm9Ev1/1O7du2UYcJyih49eqQZZi4yMjJlLPNzeffdd6lbty6RkZEpQ8Wl9tJLLxEbG0t0dDQ1a9YkOTkZgGPHjtGrVy+KFy9O0aJFeeyxxwi1EUhEREREgkEJek6RlATPPgs1a/qGOAsVJUr44vr8c5gxw+toRM6qWLFifPzxxwD89ttv7N+/3+OIst748eM5dOhQyqN169bccccdAa1bqFAh+vbtS9euXU+bN3PmTJ5//nnmzp3LwYMHmTp1asrFj3HjxvHNN9+wdu1aVq5cyaxZs5ih7wsRERHJhvJ4HYBkkXnzfNXIZ83ytu15KtcnfMrWfX9i7hLmlriSkj36cOPyCA5G5s+yGErGRPFl/xuzbHsS3iIiIqhYsSIrV67ks88+o1WrViQkJKTMnzNnDoMHD2br1q3UrVuXSZMmUaJECQDatm3LZ599xvHjx6lRowbjx4+nfPnyAMTFxdGpUydmzJjBnj17GD58OD169PBkHzNi//79vPPOO6xevTqg5Rs2bAjAihUrOHLkSJp548ePZ8CAAdSqVQuAqlWrpsxbvHgxd999N4ULF6Zw4cK0a9eOOXPm0KFDh+DsSA5lZjHARKAy4IB7gbXALCAO2ATc5Zzb61GIIiIiOY4S9Jxi9GiIjYUA73Rlha37/mRTQlPfiztLQO3afEciJGRdp3Fx/Rdk2bYke2jdujWzZ8/mv//9L88880xKgr5kyRK6du3KBx98QK1atXjyySfp3r077733HgA1atRg1KhRFC1alMcff5w2bdqwYsWKlHK/+eYbVq9ezcKFC+nUqRPdunUjT57Q/oqeOXMmNWrUoFy5chdc1qpVq2jYsCHlypXj+PHjdOnShaFDhwKcVp3dOcf69esveJvCSOAD51xrM4sA8gGPA5845xLMrD/QHwiRdlEiIiLZn6q45wRr18Knn8L993sz7nkg4uN97dBfflljo0tIu/nmm3n77beJjIzk4osvTpk+adIkOnfuzLXXXkuePHl49NFHWbBgAX/99RcA/fr1o3jx4uTOnZvOnTuzatWqNOXee++95MuXj6ZNm3Lw4EF+//33LN2v8zFlyhQ6duwYlLIOHDjA+++/T2JiIl988QXTpk1j7ty5gO/O++uvv87u3bvZtGkTs2fPPu0OvGSMmRUE6gOTAJxzx5xz+4DbgWn+xaYBLbyJUEREJGcK7dszEhzjx0PevOBv93myavmFyJSq4cOHw+zZ8Mgj8OGHIVMVXyS1iy66iObNm1O9evU00zdv3syiRYuYMmVKyrSIiAi2b99O6dKlGThwIG+++Sa7d+8mOTmZ5ORkkpKSyO2/aFakSJGUdQCOHj2aRXt0fn788UdWrVrFXXfdFZTy8uXLR5cuXShevDgArVq1YtGiRbRs2ZIHHniAdevWUaVKFQoXLkyLFi347LPPgrLdHKwssAuYYmbVgOXAg8AlzrntAM657WZW3MMYRUREchwl6NndkSMwdSq0agWXXAKcUrX8PGVK1fCiRWHQIHj4YVi4EG65JfjbEAmCZ555BiBNNevSpUszaNAgBgwYcNryr776KnPnzuWzzz4jNjaW1atXU61atbDuiXzKlCncdtttxMTEBKW8cuXKpekR3zmXcnwiIiIYPXo0o0ePBuDRRx+lWrVqQdluDpYHqAn0cc59Y2Yj8VVnD4iZdQe6A8TGxmZOhCIiIjmQqrhnd2++Cfv2wQMPeB1JYHr2hLJloW9fX8/zImGiS5cujB8/nhUrVuCcY+fOncyaNQuAgwcPEhUVReHChTl48CDDhw/3ONoLk5SUxKuvvnrG6u1z5szhiiuuYOvWraetd/ToUZKSklKenzhxAoCWLVsyceJEdu/ezbZt23jnnXe44YYbADh8+DDr168nOTmZjz/+mIkTJ9K9e/fM3cnsbwuwxTn3jf/1bHwJ+w4zKwHg/7szvZWdc6845+Kdc/HFihXLkoBFRERyAiXo2d3UqVC+PNSv73UkgYmIgIQE+P57mDbt3MuLhIg6derwwgsv0KVLFwoVKkR8fDzLly8HoGPHjsTGxnLZZZdRtWpV6tSp43G0F+Y///kPSUlJ3HKGWi4HDhxgw4YNHD9+PM30V199laioKBISEnjttdeIioriqaeeAnxt9CtVqkSZMmWoVasWHTt2pGXLlgAcOXKEpk2bEh0dTY8ePRg7dix169bN3J3M5pxzvwObzayCf1Ij4AdgHtDJP60T8K4H4YmIiORYFmpVLOPj492yZcu8DiN72LQJypSBf/4TBg5MmRzXf0FQqrhnWhnOwXXXwW+/+YaGy595w64FYz9EREKBmS13zsVnYPnq+IZZiwB+Abrgu3D/JhAL/Abc6Zzbc7ZydN6WnCgzfj/oN4lIznKm87buoGdnM2b4/obbWMFm8NxzsG0bvJh1Q66JiOQkzrlv/dXUqzrnWjjn9jrndjvnGjnnyvv/njU5FxERkeBSgp5dOQfTp0ODBhAX53U0GXf99XD77b5Efe9er6MRERERERHJdErQs6slS3zVw++5x+tIzt+wYbB/Pzz/vNeRiIiIiIiIZDol6NnVjBkQGQmtW3sdyfmrWhXuugtGjoRdu7yORkREREREJFMpQc+OkpPh7bfh73+HQoW8jubCDBniG8vdP+60iIiIiIhIdpXH6wAkE3z9ta+DtXC+e37S1VdD+/YwZgw88giUKOF1RJKDxcXFsWPHDnLnzs0ll1zCww8/TO/evQEwM2rXrs033/iGlX7sscd49tln+eyzz2jYsCHz5s3jiSeeYNOmTRQpUoTu3bvzxBNPnFbuSe+88w6NGzfO+p28AH/88Qe33347P/74I2ZGvXr1GDduHCUC+L+tVKkSv/76a8rro0eP8sADD/Dyyy9z7NgxHn74Yd566y2SkpLo2rUrI0aMwMwAmDBhAk8//TS7d++mUaNGTJkyhcKFC2fafopIeLk+4VO27vszqGWWjIkKankiIicpQc+OZs/2jSd+221eRxIcgwfD66/D00/DqFFeRyM53HvvvUfjxo1ZsmQJN9xwA1WqVKFBgwYA7Nu3j99++43Y2Fg++eQTihUrBsDGjRvp2LEj8+bNo169emzevJnExMR0yw1n+fPn55VXXuGqq64CYMiQIfTq1Ys5c+acc901a9akPE9OTiYuLo477rgDgHHjxvHNN9+wdu1aDh48SL169ahatSodOnTg22+/5bHHHuOrr76iTJky3HnnnfTt25eJEydmzk6KSNjZuu9PDV8mImFDVdyzm+RkX4J+881QsKDX0QTHFVdAly4wYQJs3ep1NCIA1K5dm0qVKrF06dKUaa1ateLtt99mxYoVVKxYkYiICACWL19OmTJlqF+/PmZGbGwsbdq08Sr0TBMVFUWlSpXInTs3SUlJJCUl8eOPP2a4nI8//phcuXKlXPhYvHgxd999N4ULFyY2NpZ27dqlJP2ff/45DRo0oFKlSuTLl48+ffoEdEFAREREJBQpQc9uli6FzZuzR/X21AYMgKQk9eguIWP58uWsWbOGq6++OmVaixYtePfdd5k9e3bK3V+AKlWqsGbNGp588klWrVqFc86LkLNM1apViYqKIiEhgV69emV4/SlTptChQ4eUKuynHi/nHOvXrz/jvL1797Jnj4bvFhERkfCjKu7ZzezZkDdvpldvLxkTRVz/BRdcRsDKloW77/bdRX/8cbj44gvatsj5atGiBcnJyZgZCQkJNG36v2qTRYsWJTIyktmzZzNo0KCU6RUqVOD9999nxIgRjBgxgksvvZSRI0dy++23pyk3T57/fSX//PPPFC9ePGt2KshWr17N/v37mTBhAvXq1cvQuvv27eOdd95h1apVKdMaNmzIq6++SqdOnTh48CCzZ89OSd4bNGjAoEGD+O677yhTpgzjx48H4MiRIxQpUiR4OyVZKrPaDH/Z/8aglpnT6X0SEQk+JejZiXO+3tsbNYJM7iDJk5PngAHw2mvw0kvw1FNZv30RfJ233XDDDQwePJivvvqKPn36pJnfpUsXli1bRlRU2gtQjRs3pnHjxvz5559MmTKFdu3asXnzZooWLZpSbri3QU+tUKFCdOzYkWuuuYZNmzal6QDvbN544w2qVavGlVdemTLtgQceYN26dVSpUoXChQvTokULPvvsMwCqV6/OM888Q4sWLTh69Ci9evVi3rx5FChQIFP2S7JGZrQZvtCLynI6vU8iIsGnKu7ZyY8/wsaN0KKF15Fkjquvhlat4OWXYd8+r6ORHCx37twMGTKExMREFi1alGZe27Ztee655864blRUFD179uSiiy5i48aNmRypt5xzbNmyhb179wa8zpQpU+jYsWOaaREREYwePZpt27axZs0azIxq1aqlzO/RowcbNmxg69atVK1albi4OAqF+xCTIiIikiMpQc9O5s/3/W2ajXsqfeIJOHDAN+yaiIfy5s1Lnz59GDZs2DmXXbVqFZMnT2bfvn0kJyfzxhtvcPz4ccqXL58FkWadL774gg8//JC//vqLQ4cO8fjjj3PFFVdwcaomKVu3buWKK65ItyO3NWvWsGrVKtq2bZtm+uHDh1m/fj3Jycl8/PHHTJw4ke7du6fM/+6770hOTmb9+vUMGDCAHj16ZN5OioiIiGQiJejZyfz5UL06lCrldSSZp0YNuPVWePFFOHzY62gkh+vatStLly49bci0U0VHRzNr1izKlStHTEwMzz//PHPmzElzl/e2224jOjo65TFt2rTMDj/okpOTGTBgAMWLF6dUqVJs3bqVd999N80yx48fZ8OGDRw4cOC09adMmULTpk1Pazt+5MgRmjZtSnR0ND169GDs2LHUrVs3Zf4TTzxBwYIFue6662jWrBl9+/bNnB0UERERyWRqg55d7NkDX37pa6ed3Q0cCNdd5+sw7pFHvI5GcpBNmzaleR0TE8PBgweB03sTB9iyZUvK84ULFwZcbriqX78+y5cvP+sycXFxZ+zF/kxNA4oVK8batWvPWOa8efMCD1JEREQkhOkOenaxcKFvDPRmzbyOJPPVrQs33ADPPQd//eV1NCIiIiIiIkGhBD27mD8fihWDa67xOpKsMWAAbN8OM2Z4HYmIiIiIiEhQKEHPDk6cgP/8x9c2O8ChjMJe48ZQrZrvLnpystfRiIiIiIiIXDC1Qc8OvvkG9u7N3r23n8oMHn0U7rnHd3EiJ+27iIiIiGQL1yd8ytZ9fwa1zJIxUXzZ/8aglilZRwl6dvDhh5Arl++uck7Spg08/jg8+6z1wSrWAAAgAElEQVQSdBEREREJO1v3/cmmhOD+jo3rvyCo5UnWUhX37OCjjyA+HgoX9jqSrJU3Lzz0ECxeDEuXeh2NiIiIiIjIBdEd9DBzajWYAn8dZuXX3zCuzp08H+DVspIxUZkVXta77z4YNszXFn3WLK+jEREREREROW9K0MPMadVg3nkHXkqmz9MP0KdBA+8C80qBAnD//b4E/ZdfoGxZryMSERERERE5LwEl6GZ2CzASyA1MdM4lnDK/PvASUBVo65ybnWpeJ2Cg/+VTzrlpwQhc/D7+GPLlgzp1vI7EOw8+CC++6Hu8/LLX0YiIiIhkWMmYqKC3HVZnYSLh55wJupnlBsYANwFbgKVmNs8590OqxX4DOgOPnrJuEeBJIB5wwHL/unuDE77w0UfQoAFERnodiXcuuwzat4fJk2HIECha1OuIRERERDIkMxJpdRYmEn4C6SSuNrDeOfeLc+4YMBO4PfUCzrlNzrnVwKkDUt8MfOSc2+NPyj8CbglC3ALw22/w889w001eR+K9f/wDjhyBsWO9jkREREREROS8BJKglwQ2p3q9xT8tEAGta2bdzWyZmS3btWtXgEULH33k+6sEHSpXhr//HcaMgb/+8joaERERERGRDAskQbd0prkAyw9oXefcK865eOdcfLFixQIsWvjkE7j0UqhUyetIQsNDD8GOHfDmm15HIiIiIiIikmGBJOhbgNKpXpcCtgVY/oWsK2fjnG/87wYNwNK7DpID3XQTXH01vPSS7/iIiIiIiIiEkUAS9KVAeTMrY2YRQFtgXoDlLwSamFlhMysMNPFPkwu1YQNs2+ZL0MXHzNej+4oV8OWXXkcjIiIiIiKSIedM0J1zJ4De+BLrH4E3nXNrzGyYmTUHMLNrzGwLcCcwwczW+NfdA/wTX5K/FBjmnyYXavFi318l6Gndcw8ULuy7iy4iIiIiIhJGAhoH3Tn3PvD+KdMGp3q+FF/19fTWnQxMvoAYJT2LF0OxYr4q3fI/+fJB9+7w7LOwaRPExXkdkYiIiIiISEACStAlBC1eDPXrq/15enr1guee8/Xo/uyzXkcjIiIhoGRMVNDHhC4ZE5UpY1eLiEjOpQQ9HG3a5BsD/dFHvY4kNJUuDXfcAf/+Nzz5JERHex2RiIh4LDMS6WAn/CIiIoF0EiehRu3Pz+2hh2D/fpg+3etIREREREREAqIEPRwtXgxFikDlyl5HErrq1IFrroGRIyE52etoREREREREzkkJejj6/HOoVw9y6e07IzPfXfSff4YPPvA6GhERERERkXNSG/Qwc8nBP3xjoPfq5XUooa91a+jb13cX/dZbvY5GRETEM9cnfMrWfX8GtcySMVFBLU9ERJSgh53am9f4ntSv720g4SAiAnr2hIED4YcfoGJFryMSEQkZZrYJOAgkASecc/FmVgSYBcQBm4C7nHN7vYpRgmfrvj/ZlNDU6zBEROQcVEc6zNTc9pNvrO9q1bwOJTzcfz9ERsLLL3sdiYhIKLrBOVfdORfvf90f+MQ5Vx74xP9aREREsogS9DBTa+uPcO21kEeVHwJy8cXQrp2vN/d9+7yORkQk1N0OTPM/nwa08DAWERGRHEdZXjg5fJiKO36Brm28jiS89OkDU6fCtGnw4INeRyMiEioc8KGZOWCCc+4V4BLn3HYA59x2Myue3opm1h3oDhAbG5tV8YackjFRQR8LvWRMVKaM2S4iIuFBCXo4WbqUPC4Z6tb1OpLwUrOm75iNGeNL1tX7vYgIwPXOuW3+JPwjM/sp0BX9yfwrAPHx8S6zAgx1mZFIBzvhFxGR8KJMJZwkJvr+1qnjbRzhqHdvWLcOPvrI60hEREKCc26b/+9OYC5QG9hhZiUA/H93ehehiIhIzqMEPZx89RXri5SCokW9jiT8tG4Nl1wCo0d7HYmIiOfMLL+ZFTj5HGgCfA/MAzr5F+sEvOtNhCIiIjmTEvRw4Rx89RXLS17tdSThKSICuneHBQvgl1+8jkZExGuXAF+Y2SpgCbDAOfcBkADcZGbrgJv8r0VERCSLqA16uPj5Z9izhxW1r0JdxJ2n+++Hp5+GcePg2We9jkZExDPOuV+A08brdM7tBhplfUQikpNdn/ApW/f96XUY56ROHENfZn2WsvK9V4IeLr76CkB30C9EyZLQqhVMmgRDh/rGkxcRERERT23d9yebEpp6HcY5qRPH0JdZn6WsfO9VxT1cfPUVxMSwoWgpryMJb717w9698PrrXkciIiIiIiKShhL0cJGYCHXr4kxv2QX529+galVfZ3Eux44MJCIiIiIiIUhV3MPB/v3www9w111wxOtgwpyZ7y569+7w5ZdeRyMiIiKplIyJCnpV0pIxUUEtT0SCIzPai2eH/3cl6OFg+XLf3d7atWFRktfRhL+774bHHvPdRY+7x+toRERExE8dcInkHOHS90BWU33pcLB0qe/vNdd4G0d2kT8/dO0Kb79N8YO7vY5GREREREQEUIIeHpYsgbJloWhRryPJPh54AJKSuHvVB15HIiIiIiIiAihBDw9Ll/qqt0vwlCsHt97K3d9+AMeOeR2NiIiIiIiI2qCHvN9/h82bVb09M/TuTfEFC+Dtt6FdO6+jEREREREJSerQLesoQQ91an+eeZo04ZfCl1F29Ggl6CIiIiIiZ6AO3bKOqriHuqVLIVcuqFnT60iyn1y5eLVmU/jqK1ixwutoREREREQkh9Md9FC3dClUquTreVyC7qvrm3H481eZ37k//W598LzKKBkTpWFhRERERETkgilBD2XO+Xpwb9nS60iyrYVDm8POzrSZOpU2fd84r57y4/ovyITIREREREQkp1EV91C2cSPs2aP255mtd284ehQmTfI6EhERERERycGUoIcydRCXNSpVghtugLFjISnJ62hERERERCSHUoIeypYsgchIqFLF60iyv9694ddfYYGqq4uIiIiIiDeUoIeypUuhRg3Im9frSLK/5s2hdGl4+WWvIxERERERkRxKCXqoSkryDf1Vu7bXkeQMefLAAw/Axx/Djz96HY2IiIiIiORAStBD1c8/w+HDUKuW15HkHN26QUQEjBnjdSQiIiIiIpIDKUEPVStW+P7WrOltHDlJsWLQti1MmwYHDngdjYiIiIiI5DBK0EPVypVw0UVw1VVeR5Kz9OkDhw7B9OleRyIiIiIiIjlMHq8DkDNYsQKqVvW1jZasEx8P114Lo0dDz56QKzyvYU2dOpUuXbpQqFAhNm7cSOHChVPmnThxgrx58/Lkk08yZMiQDJebnJzMvffeG+SI/2fIkCEMHTr0nMtNmTKFhg0bUqZMGf7973/TrVu3sy4fFxdHw4YNmTp16gXFZ2bndewy06ZNm5g6dSodO3akbNmymbqt5ORkHnnkEWbNmsWOHTto3rw5L730EmXKlGHKlCl07twZOPtn5eR7fPz4cfLoO04kjZIxUcT1D/6IIiVjooJepoiEpsz4HtF3SNbRL6NQ5JzvDnqbNl5HkjP16QMdOsAnn8BNN3kdzQXZv38/I0aMICEhISjlTZ06lRMnTmRqgt6tWzduueWWlNcLFizgqaee4q233qJUqVIp08uVK8fhw4cDLnfu3LkULFjwguNLTExME0co2LRpE0OHDuVvf/tbpifos2fPZuTIkTz//PPUrVuXokWLUqJECRITEylXrlzKclnxWRHJjr7sf6PXIYhImNP3SHhTgh6KNm2CffvU/twrrVvDI4/4hlwL8wS9SZMmvPzyyzz00ENceumlXocTkFKlSqVJgH/66ScAqlevzhVXXJFm2Ywk6DVq1AhKfHXq1AlKOeHqR/8oBw899BC5UtUwyenHRURERCQYwrP+bnZ3soO4ICUUkkGRkdC9O8yfDxs3eh3NBRk4cCAA//rXv8657JIlS2jcuDHR0dHkz5+fRo0asWTJkpT5DRs2ZPHixXz55ZeYGWZGw4YNU+Zv3LiR9u3bU6xYMSIjI6levTpz584N+j6lJykpicGDB1OiRAliYmK47bbb2LJlS5pl4uLiUqpfA/z+++906tSJyy67jMjISEqUKEGzZs3YuXPnWbdlZmmqt//888+0bNmS4sWLc9FFFxEbG8udd97JiRMnzljGpk2bMDMmTJhwzriPHz/OwIEDiYuLIyIigri4OAYOHMjx48cBWLRoETfccAMAN910U8p7s2jRojNu/8MPP+TWW2+lRIkS5MuXj8qVK/P888+TlJR01n2Pi4tL2ffcuXNjZkydOjVlf042HzjXZ0VERERE0hdQgm5mt5jZWjNbb2b905kfaWaz/PO/MbM4//S8ZjbNzL4zsx/NbEBww8+mVq6E3LmhShWvI8m5evTwtT8fN87rSC5IiRIl6N27N6+88gq//vrrGZdbvXo1DRo0YO/evUydOpXp06dz4MABGjRowKpVqwAYO3YsNWrUoGrVqiQmJpKYmMjYsWMB2Lx5M9deey2rVq3ixRdfZN68edSsWZM77riDefPmpWznZCIX7PbbTz/9NOvXr2fy5MmMHDmSxMRE2rdvf9Z17rnnHhITE3n22Wf56KOPGDVqFKVKleLIkSMZ2nazZs3YunUr48aNY+HChSQkJBAZGUlycnJQ4u7UqRMJCQl07NiR+fPn06VLF0aMGEGnTp0AqFmzJmP8QwOOGjUq5b2peZYaOL/88guNGjVi8uTJLFiwgE6dOjFkyBCeeOKJs8Y7d+7clIscJ7fTtGnT05Y722dFRERERM7snFXczSw3MAa4CdgCLDWzec65H1It1hXY65y7wszaAiOANsCdQKRzroqZ5QN+MLM3nHObgr0j2cqKFVCxoq8Xd/FGyZLQqhVMnAhDhkC+fF5HdN769evHhAkTGDp0KJMnT053mWHDhhEZGcknn3xCTEwM4LsbGxcXx9ChQ5kzZw4VK1akYMGCnDhx4rTqzEOGDME5x+LFiylatCgAN998M5s3b2bw4ME0b94c8N19zp07d5qq0cFw+eWX8/rrr6e83rVrF3379mXbtm1cdtll6a6TmJjI8OHD0yTEd955Z4a2+8cff7Bu3TrefffdlH0EuPvuu4MS9/fff88bb7yRplO6Jk2akDt3bgYNGkT//v2pWrUqFStWBODqq68OqKp5jx49Up4756hXrx7Hjh3jueeeY/jw4Wd8f2rUqEHJkiWBtFXaT21qcLbPCkDevHlTPgsiIpJ5MqvTwWBTB2Qi/xNIG/TawHrn3C8AZjYTuB1InaDfDgzxP58NjDYzAxyQ38zyAFHAMUADTJ/LypVw881eRyF9+sBbb8Ebb0DXrl5Hc96KFCnCP/7xD4YOHUq/fv3SdOR10ueff06zZs1SknOAggUL0rx5c957771zbuODDz7g1ltvpVChQmmqdt9888307duXAwcOULBgQS6//PKzVv0+X6fexa3ir33y22+/nTFBv+aaa3j22WdxznHjjTdSuXJlfF9bgStatChly5alf//+7Nixg4YNG1K+fPmgxf35558D0KFDhzTLdejQgUGDBrF48WKqVq2aoZgBtm/fzpAhQ/jggw/Ytm1bmvdk586dmd5fQYECBYiOjs7w8RYRkYxRZ2Ei4SeQ21glgc2pXm/xT0t3GefcCWA/UBRfsn4Y2A78BjznnNtz6gbMrLuZLTOzZbt27crwTmQr27fD77+r/Xko+NvffEPdjR7t61k/jD388MMUKVKEwYMHpzt/z549lChR4rTpl156KXv37j1n+Tt37mT69OnkzZs3zaNv374A7N69+8J24ByKFCmS5nVkZCQAR48ePeM6s2bNonnz5jzzzDNUrVqVkiVLMmzYsICqpp9kZnz00UfEx8czYMAArrzySsqWLcu4AJtGnCvuPXt8X5envjcnE+iT8zMiOTmZ5s2bM3/+fAYOHMinn37K0qVLU6q3n+2YBUuhQoWIjo7O9O2IiIiIhJtAEvT0bnGcmq2caZnaQBJwGVAG+IeZnTYGkHPuFedcvHMuvlixYgGElI2tXOn7qx7cvWcGvXvDt9/Cl196Hc0FiY6OZsCAAbz11lt8++23p80vUqQIv//++2nTf//999OSyPQULVqU1q1bs3Tp0nQfZ7qL7aXixYszZswYtm7dyk8//UTnzp158sknmTBhQobKKVu2LNOnT2fXrl2sXLmSG2+8kZ49e/Kf//zngmM8eexPfW9Ovj7ZnCAjNmzYwLJlyxgxYgT33Xcf9erVIz4+Pkurm3fq1Ilt27Zl2fZEREREwkUgCfoWoHSq16WAU39ZpSzjr85eCNgD3A184Jw77pzbCXwJxF9o0NnayR7cq1f3Ng7xad8eYmJ8d9HDXM+ePSlZsmRKz+6pNWjQgAULFnDw4MGUaQcPHuS9996jQYMGKdMiIyP5888/T1v/lltuYfXq1VSqVIn4+PjTHifvDIeqChUqMHz4cAoXLsz3339/XmWYGdWrV+eFF14AOO9yUjt57GfOnJlm+owZMwCoX78+8L877+m9N6c62Qle3rx5U6YdP348pcxgOdNnRURERETOLJA26EuB8mZWBtgKtMWXeKc2D+gEJAKtgU+dc87MfgNuNLPXgHxAHeClYAWfLa1cCeXLQ4ECXkci4OscrmtXGDkStm2DELwTHKjIyEgGDx5M9+7dT5s3aNAg5s+fT6NGjejXrx9mxogRIzhy5EiaavEVK1Zk7NixzJo1i3LlylGgQAEqVKjAsGHDqF27NvXr16d3797ExcWxd+9evv/+e3755ZeUzul+/fVXypUrx+DBg89Y3T4r7N+/n8aNG9O+fXuuuuoq8ubNy7vvvsvevXtp0qRJwOWsXr2aBx98kDZt2nDFFVeQlJTE1KlTyZMnDzfeeOHt/ipVqkS7du0YMmQIJ06c4LrrriMxMZF//vOftGvXLqX9+ZVXXkmePHmYPHkyRYoUITIykgoVKlAgne+Rq6++mssvv5wnnniC3LlzkzdvXl588cULjvVUZ/qsgK9TwmHDhvHXX3+pozgRERGRVM55B93fprw3sBD4EXjTObfGzIaZ2cluiycBRc1sPfAIcHIotjFANPA9vkR/inNudZD3IXtZsULV20NNz56QlAQZrPocirp06ZJuJ2ZVq1Zl0aJFFCxYkE6dOnHPPfcQHR3N4sWLqVatWspy/fr1o1GjRnTr1o1rrrmG+++/H4DY2FiWLVtGtWrVePzxx7npppt44IEHWLx4cZpE1TlHUlJShtp5Z4aLLrqImjVr8u9//5vWrVvTsmVLEhMTmTFjBrfffnvA5Vx66aXExsbywgsv0Lx5c9q1a8e2bduYP38+tWrVCkqs06ZNo1+/fkyePJlbb72VSZMm0a9fP6ZNm5ayTNGiRRk9ejSrVq2iQYMGXHPNNSxfvjzd8iIiInjnnXe49NJL6dixI7169aJ+/fr073/aCJoX5EyfFfC1g09KSsKFed8OIiIiIsFmofYDKT4+3i1btszrMLyxdy8UKQIJCdCvX7qLxPVfwKaE08cdlkx2222wdCn89htERKSZpfdEREKBmS13zmV5M7LMOG/re1VEQlFmfDfp+y48ZMb7dKbzdnAHI5YLow7iQlfv3rBjB8ye7XUkIiIiIiKSTSlBDyUne9dWB3Gh56ab4Mors0VncSIiIiIiEpqUoIeSVaugRAnI6UPNhaJcuaBXL0hMhDO07RURCTdmltvMVprZfP/rMmb2jZmtM7NZZhZxrjJEREQkeJSgh5LVqyFVh1wSYjp1gvz5dRddRLKTB/F1AHvSCOBF51x5YC/Q1ZOoREREcigl6KHi+HH44QfwD5skIahQIV+S/sYb8McfXkcjInJBzKwU0BSY6H9twI3Ayc42pgEtvIlOREQkZ1KCHirWroVjx3QHPdT16gV//QUTJ3odiYjIhXoJeAw4Oe5hUWCff3hVgC1ASS8CExERyamUoIeK1f7h4XUHPbRVrAg33gjjxsGJE+deXkQkBJlZM2Cncy51pxqWzqLpjsVqZt3NbJmZLdu1a1emxCgiIpITKUEPFatX+8bXrlDB60jkXPr08Y2H/t57XkciInK+rgeam9kmYCa+qu0vATFmlse/TClgW3orO+decc7FO+fii6ljUxERkaBRgh4qVq3y3Z3Nm9frSORcmjWD2Fh1FiciYcs5N8A5V8o5Fwe0BT517v/bu/M4nev9/+OPF2MYBoPklOXMZM1SlB9akRYhpLK10Jdw2k6LCqcOp07lyBFOnVLRdk4hCaVzJEspHUK2Ibtsycg6hQbv3x+fy2SZyTXmuq7Pdc31vN9un9t1XZ/lmuf17tJnXvN5f95vdyswC7g5sFs3YLJPEUVEROKSCvRosXSpurfHioQEuPtumDkTli/3O42ISCg9BjxkZmvx7kkf7XMeERGRuKICPRrs3AnbtqlAjyV33QVJSTB8uN9JRETyxTk32znXJvB8vXOukXOumnPuFufcIb/ziYiIxBMV6NHg2ABxGsE9dpQt60259q9/UfbnvX6nERERERGRAkAFejTQCO6x6Y9/hEOHuPWbj/1OIiIiIiIiBYAK9GiwZAn87ndw9tl+J5G8qFULrr+eO76Z6s2NLiIiIiIikg8q0KOBBoiLXQ88QPmf9sC4cX4nERERERGRGKcC3W+HD0N6ugr0WHXNNawuVwWefx6c8zuNiIiIiIjEMBXoflu92userQHiYpMZYxq2hcWL4fPP/U4jIiIiIiIxTAW63zRAXMz7oE5zKFfOu4ouIiIiIiJyhlSg+23JEihSxBtwTGLSoSJFoU8fmDIF1q3zO46IiIiIiMSoBL8DxL2lS+H88yEx0e8kkh933w1DhsDIkTBihN9pRERERCTEKqYkkdpvasjfU+R4KtD9tmQJNG/udwrJr3PPhU6dYMwYePJJKF3a70QiIiIiEkJf9rvK7wgSB9TF3U+7dsHWrbr/vKB44AHIzITRo/1OIiIiIiIiMUgFup+WL/ce69XzN4eExsUXwxVXeN3cDx/2O42IiIiIiMQYFeh+Sk/3HuvU8TeHhM6DD8J338HkyX4nERERERGRGKMC3U/p6VCqFFSq5HcSCZW2bSEtDYYP9zuJiIiIiIjEGA0SF0GXDZ7J1j0Hsl+PnfIZRZLP5ab+Hwf9HhrpMcoVLgz33+9dSZ8/Hxo18juRiIiIiIjECBXoEbR1zwE2Dm7964rR3aF9+xPXSezr0QMGDYLnnoP33vM7jYiIiIiIxAh1cffLjh2wc6fuPy+ISpaEP/wBJk6Edev8TiMiIiIiIjFCBbpfNEBcwXb//ZCQAMOG+Z1ERERERERihAp0v6hAL9jOOQduuw3GjIGMDL/TiIiIiIhIDFCB7pfly6FMGa+Qk4Kpb184eBBefNHvJCIiIiIiEgNUoPslPd27em7mdxIJl/PPhzZt4IUX4Oef/U4jIiIiIiJRTgW6H5z7tUCXgu3RR+HHH+GNN/xOIiIiIiIiUU4Fuh+2b4fdu1Wgx4PLL4fGjeHvf4cjR/xOIyIiIiIiUUwFuh80QFz8MINHHoH1671p10RERERERHKhAt0Py5d7j3Xr+ptDIqN9e6hWDZ57zru9QUREREREJAcq0P2Qng5nnQVnn+13EomEwoXh4Yfh66/h88/9TiMiIiIiIlFKBbofNEBc/OnWDcqX966ii4iIiIiI5EAFeqRpBPf4lJQE994LU6f+eouDiIiIiIjIcYIq0M2spZmtMrO1ZtYvh+1FzWxcYPs8M0s9btsFZvaVmaWb2TIzKxa6+DFo61bYt08Fejy65x4oUQIGD/Y7iYiIiIiIRKHTFuhmVhh4EbgeqA10MbPaJ+3WA9jtnKsGPA/8LXBsAvAvoI9zrg7QDMgKWfpYpAHi4le5ctCnD7z7Lqxb53caERERERGJMsFcQW8ErHXOrXfO/QKMBdqdtE874M3A8wlACzMz4FpgqXNuCYBz7kfnXHxPBq0p1uLbww9DQgIMGeJ3EhERERERiTIJQexTEdh83OstQOPc9nHOHTazvUA5oAbgzGwaUB4Y65w7pTIxs15AL4AqVark9TPElvR0qFDBu5oqBULFlCRS+00Nev+/1m7BLaPHcGXCpfxQ8qzs9/iy31XhiigiIiIiIjEgmALdclh38mTOue2TAFwO/D/gZ2CGmS10zs04YUfnXgFeAWjYsGHBnihaA8QVOHkurHvXhurVmVd0CQweBpCnAl9ERERERAqmYLq4bwEqH/e6ErAtt30C952XBnYF1n/mnNvpnPsZ+Bi4KL+hY5W5o7BihQr0eJeWBl27wqhRsHOn32lERERERCRKBFOgfw1UN7M0M0sEOgNTTtpnCtAt8PxmYKZzzgHTgAvMrHigcG8KrAhN9NhTcV8GZGZqgDiB/v3hwAEYMcLvJCIiIiIiEiVOW6A75w4D9+IV2yuB8c65dDN70szaBnYbDZQzs7XAQ0C/wLG7gWF4Rf5iYJFzLm778lbfucl7oivocv75cOON8I9/eNPuiYiIiIhI3AvmHnSccx/jdU8/ft2fj3t+ELgll2P/hTfVWtyrsfM774kKdAEYMAAmToR//hOo53caERERERHxWTBd3CVEauzcBOeeCykpfkeRaHDxxXDddTBsGMWyDvqdRkTiiJkVM7P5ZrbEzNLN7C+B9WlmNs/M1pjZuMCtbSIiIhIhKtAjqPrOTbp6LicaMAAyMui85BO/k4hIfDkEXOWcuxCoD7Q0sybA34DnnXPVgd1ADx8zioiIxB0V6JFy9CjVd27WAHFyoiuvhCuuoM+8CXBQV9FFJDKcJzPwskhgccBVwITA+jeB9j7EExERiVsq0CNlwwaSDh/SFXQ51aBB/C5zF7z2mt9JRCSOmFlhM1sM7ACmA+uAPYHBYcGbKrViLsf2MrMFZrYgIyMjMoFFRETigAr0SElP9x5VoMvJmjdnXqU68OyzuoouIhHjnDvinKsPVAIaAefntFsux77inGvonGtYvnz5cG7Am7sAACAASURBVMYUERGJKyrQI+VYgV67tr85JPqYMfzyW2HbNnj1Vb/TiEiccc7tAWYDTYAUMzs2w0slYJtfuUREROKRCvRISU9nS6nyUKqU30kkCn1VpZ53P7quootIBJhZeTNLCTxPAq4GVgKzgJsDu3UDJvuTUEREJD6pQI+U9HTWnFXF7xQSrcxg0CD4/ntdRReRSDgHmGVmS4GvgenOuY+Ax4CHzGwtUA4Y7WNGERGRuJNw+l0k344cgZUrWX1ha5r7nUWiV/Pm0LSpdxW9Z09ISvI7kYgUUM65pUCDHNavx7sfXURERHygK+iRsG4dHDqkK+hyerqKLiIiIiISt1SgR0JggLjVKtDldJo1866iDx4MBw74nUZERERERCJIBXokBAr0NeVUoEsQjl1Ff+UVv5OIiIiIiEgEqUCPhPR0SEvjQGIxv5NILGjWzLsf/ZlnIDPT7zQiIiIiIhIhKtAjYflyqFPH7xQSS55+GnbsgBEj/E4iIiIiIiIRogI93LKyYNUqFeiSN5dcAjfcAM89B7t2+Z1GREREREQiQAV6uK1d6xXpKtAlr55+GvbtgyFD/E4iIiIiIiIRoAI93AIDxFG3rr85JPbUqwddu8LIkd6gcSIiIiIiUqCpQA+39HQoVAhq1fI7icSiv/zF64Hx1FN+JxERERERkTBTgR5uy5fDeedBUpLfSSQWVa0KPXvCq6/CunV+pxERERERkTBSgR5u6em6/1zy54knoEgRb350EREREREpsFSgh9Mvv8CaNSrQJX/OPRfuuw/+/W9YtszvNCIiIiIiEiYq0MNp9Wo4fFgDxEn+PfYYlCwJAwb4nURERERERMJEBXo4HRvBXVfQJb/KloX+/eGjj2D2bL/TiIiIiIhIGKhAD6fly6FwYahZ0+8kUhD88Y9QuTL07QtHj/qdRkREREREQkwFejilp0O1alC0qN9JpCBISoJnnoGFC+Hdd/1OIyIiIiIiIaYCPZw0gruEWteucNFF3r3oBw74nUZEREREREJIBXq4HDwIa9dqgDgJrUKF4O9/h02bYORIv9OIiIiIiEgIqUAPl1WrvPuEdQVdQq1ZM2jb1uvunpHhdxoREREREQkRFejhsny596gCXcLhb3+Dn36CJ5/0O4mIiIiIiISICvRwSU+HhASoXt3vJFIQ1aoFvXvDSy95vTVERERERCTmqUAPl/R0qFEDEhP9TiIF1cCBUKIEPPSQ30lERERERCQEVKCHS3q6BoiT8Dr7bK9I//hjmDrV7zQiIiIiIpJPKtDD4eefYf163X8u4XfvvV539wcfhEOH/E4jIiIiIiL5kOB3gAJp5UpwTgW6hF9iIgwfDi1b8mL7+3juwnb5eruKKUl82e+qEIUTEREREZG8UIEeDunp3qMKdImE666Dtm3p9t9/cc/oQXDuuWf8Vqn91FVeRERERMQv6uIeDunp3pXNatX8TiLxYtgwihzJgn79/E4iIiIiIiJnSFfQwyE93bsvOEHNK8GpmJKU76vXf7miI93efhv69IFLLw1RMhERERERiRRVkOGQng6XXOJ3CokhIbnvO7Mp1PoM7r8f5s2DwoXz/54iIiIiIhIx6uIeapmZsHGj7j+XyEtOhueeg4UL4ZVX/E4jIiIiIiJ5pAI91Fas8B41B7r4oXNnaNEC+veH7dv9TiMiIiIiInkQVIFuZi3NbJWZrTWzU0ahMrOiZjYusH2emaWetL2KmWWaWd/QxI5iy5d7jyrQxQ9m8M9/woED8NBDfqcREREREZE8OG2BbmaFgReB64HaQBczq33Sbj2A3c65asDzwN9O2v488J/8x40By5dDUhKkpfmdROJVjRreFfR334Xp0/1OIyIiIiIiQQrmCnojYK1zbr1z7hdgLNDupH3aAW8Gnk8AWpiZAZhZe2A9kB6ayFEuPR1q14ZCuntAfNSvH1SvDnffDQcP+p1GRERERESCEEwVWRHYfNzrLYF1Oe7jnDsM7AXKmVkJ4DHgL7/1A8ysl5ktMLMFGRkZwWaPTsuXq3u7+K9YMa+r+9q18MwzfqcREREREZEgBFOgWw7rXJD7/AV43jmX+Vs/wDn3inOuoXOuYfny5YOIFKV27YJt21SgS3S4+mq49VYYPBi+/dbvNCIiIiIichrBFOhbgMrHva4EbMttHzNLAEoDu4DGwBAz2wg8AAwws3vzmTl6pQd68WuKNYkWf/87FC8Of/gDuJP/riYiIiIiItEkmAL9a6C6maWZWSLQGZhy0j5TgG6B5zcDM53nCudcqnMuFRgOPOOceyFE2aPPsQJdV9AlWlSo4M2NPnu25kYXkWxmVtnMZpnZSjNLN7M/BtaXNbPpZrYm8FjG76wiIiLx5LQFeuCe8nuBacBKYLxzLt3MnjSztoHdRuPdc74WeAg4ZSq2uLB8OZQqBZUq+Z1E5Fc9e3pzoz/yCGza5HcaEYkOh4GHnXPnA02AewIztPQDZjjnqgMziNfzuYiIiE8SgtnJOfcx8PFJ6/583PODwC2neY9BZ5Avtixf7nVvt5xuyRfxiRm8+qrXs6N3b/j4Y31HReKcc+574PvA8/1mthJvwNd2QLPAbm8Cs/EGexUREZEI0FxgoeKcRnCX6JWW5g0W99//wptvnn5/EYkbZpYKNADmARUCxfuxIv7sXI4pOLOviIiIRBEV6KGyYwf8+KMKdIle99wDl18ODz7ozTYgInHPzJKB94EHnHP7gj2uwMy+IiIiEmVUoIfK8uXeo0Zwl2hVqBCMHg0HD2pUdxHBzIrgFef/ds5NDKz+wczOCWw/B9jhVz4REZF4FNQ96AKXDZ7J1j0Hct3efcEUBgENJ33PzulTc9ynYkpSeMKJBKtGDXjqKW/AuHffha5d/U4kIj4wM8Mb4HWlc27YcZuOzcoyOPA42Yd4IiIicUsFepC27jnAxsGtc9+h12RYchYLRt6qAbgkuj34ILz/vtfl/corNeuASHy6DLgdWGZmiwPrBuAV5uPNrAewidMMACsiIiKhpQI9VDSCu8SKwoXhrbegfn3o3h0++cTr/i4iccM59wWQ2wmrRSSziIiIyK/0W3koOAfp6RogTmJH9erw/PMwYwaMHOl3GhERERERQQV6aGzZAvv2qUCX2HLXXdCmDfTr9+sghyIiIiIi4hsV6KGgEdwlFpnBa69BqVJw221w6JDfiURERERE4poK9FBQgS6xqkIFr0hfsgQGDvQ7jYiIiIhIXFOBHgrp6XDuuVC2rN9JRPKubVvo2ROGDOGyjYtPv7+IiIiIiISFCvRQODaCu0isGj4catVi+EdDYft2v9OIiIiIiMQlFej5dfQorFihAeIktpUoAePHU/LQz3D77d73WkREREREIkoFen5t2AAHDqhAl9hXty4Dr+4Nn34Kgwf7nUZEREREJO6oQM+vZcu8R3VxlwJg3AXXQufO8MQTMGeO33FEREREROKKCvT8WrrUm65KV9ClIDCDUaMgLQ26dIGdO/1OJCIiIiISN1Sg59fSpVCtmncPr0hBUKoUjB8PGRnQtSscOeJ3IhERERGRuKACPb+WLoULLvA7hUhoXXQRvPgiTJ/udXcXEREREZGwU4GeHz/9BGvXqkCXgqlnT7jrLnj2WZg40e80IiIiIiIFngr0/Fi+HJxTgS4F1z/+AY0aQbdu8O23fqcRERERESnQVKDnx9Kl3qMKdCmoihaFCRMgKQluvBH27fM7kYiIiIhIgaUCPT+WLoWSJSE11e8kIuFTuTKMGwdr1sAdd8DRo34nEhEREREpkFSg58fSpVCvHhRSM0oB17w5/P3vMHky/OlPfqcRERERESmQVFmeKedgyRJ1b5f4cf/90Ls3DB4Mb7zhdxoRERERkQJHBfqZ2rwZ9u5VgS7xw8wbNK5FC+jVC+bM8TuRiIiIiEiBogL9TGmAOIlHRYrAe+9BWpo3aNy6dX4nEhEREREpMFSgn6ljBXq9ev7mEIm0MmXgo4+82zzatIHdu/1OJCIiIiJSIKhAP1NLl3pXEUuV8juJSORVrw4TJ8L69dCuHRw44HciEREREZGYpwL9TGmAOIl3TZvC22/DF1/ArbfCkSN+JxIRERERiWkq0M/EgQOwerUKdJGOHWH4cPjgA7jnHq/bu4iISIzp06cPTz31lN8xRERI8DtATFqxAo4eVYEuBU7FlCRS+03N41FVeazxzfxh1Cj+np7JPy7rku8MX/a7Kl/vISIisad79+688847JCYmZq97+OGH+ctf/hLU8c2aNeO2226jZ8+eef7ZL7/8cp6P8cNf//pX3n33XVauXMmYMWPo3r179rYtW7Zw2223MX/+fGrVqsVbb71F3bp1AcjKyqJPnz6MHz+eMmXKMHToUDp27Jh97MiRI3n66afJysqid+/ePPvss5H+aCISoAL9TCxe7D1eeKG/OURC7IwLY9cKunfn4bfe4uHbrvDmSz9Def8DgYiIFBSPPvoof/3rX/2OEbXOO+88Ro4cSb9+/U7Z1qtXL+rWrct//vMfRo4cSadOnUhPTwfg+eefZ9myZWzevJnFixfTpk0bLrnkEipXrsz8+fMZNGgQc+bMISUlhcsvv5wGDRqcUMCLSOSoi/uZ+OYbKFkSqlb1O4lIdDCD116DVq3gD3/w7k0XEREJkTfeeINLL72U7t27U7JkSS6//HJ27twJwDPPPENycjJz5szh3nvvJTk5mYsvvviE41NTU3nhhRe4+OKLKVGiBG3btgXgo48+Ijk5mSJFivD444+f8nMnTpxI3bp1KVOmDK1ateL777/P3jZ16lRq1qxJyZIlqV69OtOmTQtjC3i6du1KixYtKFq06Anr9+3bxyeffEK/fv1ISkriwQcf5LvvvmNpYNah9957j/vvv5+UlBSaNWvGJZdcwgcffJC9rUOHDtSpU4eKFSvSs2dPxo4dG/bPIiI5U4F+JhYtggYNoJCaTyRbkSIwYQI0bw7du3vzpYuIiITI/Pnz6dq1Kzt27CArK4sxY8YAMGDAADIzM7niiit44YUXyMzMZOHChaccP2rUKN566y327t3LgAEDAGjTpg2ZmZnceuutOf68Hj16MHr0aDIyMmjQoAG9evXK3t6zZ0+efPJJ9u/fzyeffELFihXD9MlPb+3atRQrVozk5GSuuOIKtmzZQtWqVVm1ahUAq1atolatWtx2221MnDiR2rVrn7JtxIgR9O/f/4RtIhJ5qjDz6sgRr4t7gwZ+JxGJPklJMHkyXHIJdO3qzZcuIiISpKFDh5KSkpK9TJw4MXtbjRo1uPbaa0lKSqJFixasXr06T+/dq1cv6tSpQ0JCAk2aNDnt/qNHj6Z79+40btyYhIQE+vbty9SpUzl06BAAhQoVYt26dezbt4+0tLTs+7398NNPP5GcnMz+/ftZs2YNu3fvpmTJkmRmZp6wfdWqVWzdujXHbZs3b2bNmjUnbBORyFOBnlerVnmjuF90kd9JRKJTcjJMnQr168NNN8H06X4nEhGRGNG3b1/27NmTvXTo0CF7W9myZbOfJyYmcvDgwTy9d/Xq1fO0/+bNmxk1alT2HwvS0tJITEzM7ub+3nvv8dVXX1GlShWaNGnC8uXL8/T+oVSiRAkyMzOpXLky27dv5+KLL2b//v0kJyefsP3rr7/mvvvuy3Hb0KFDmTBhwgnbRCTyVKDn1aJF3qMKdJHclS4N06ZBrVrQtq33XEREJIwKnebWw4SEvI2NXLlyZZ544okT/mBw8OBBUlNTAbj00kv58MMP+eGHH6hVq1aO97BHSrVq1Th48CBbtmwB4JdffmHdunXUrFkT8HoffPvtt9n7r1ixIqhtIhJ5GsU9r775BooV8woPEcld2bLw6adwzTVekf7++9CmzWkPO7Op3k59D03VJiISX373u9+F9Cr2nXfeSadOnbjuuuto0KABGRkZzJo1i06dOnH06FHeeecd2rVrR7FixQAoVapUyH52brKysjhy5AjOObKysjh48CCJiYmUKlWKa6+9lsGDBzN06FBGjBhBlSpVqFevHgAdO3Zk5MiRtGnThiVLlvDVV19l38N/yy23cP311/Pggw+SkpLC6NGjNc2aiI9UoOfVokXe9Gp5/CusSFwqXx5mzoTrroMbb4SxY71u778hFIW1pmoTEYlNQ4YMYfjw4dmvO3fuzGuvvRbUsQ8//DB33HEH55xzDjVr1mT27NmnPebaa69l7ty5HDp0CDNj+PDh3Hzzzbzxxhs0adKEYcOGceedd7JhwwZSUlLo3LkznTp1AuBf//oX9957L845GjduzOjRo8/oM+fFXXfdxZtvvgnA3Llz6dWrF7NmzaJZs2aMGjWK2267jTJlylCrVi3GjRuHmQHw4IMP8u2331K5cuXsIrxy5coANG7cmIEDB9K8efPsedCPfUYRiTxzzp1+J7OWwAigMPCac27wSduLAm8BFwM/Ap2ccxvN7BpgMJAI/AI84pyb+Vs/q2HDhm7BggVn8lnCKrXfVDY+cz2UKQO33gr//KffkURix969cP31MH++NwVbly5h/XGp/aaycXDrsP4MkWhjZgudcw0j/XPDcd7Wv2EREYkm4Tgv5XbePu096GZWGHgRuB6oDXQxs9on7dYD2O2cqwY8D/wtsH4ncINzrh7QDYjtyZE3bIB9+3T/uUheHbsn/fLLvT9wvfKK34lERERERKJOMIPENQLWOufWO+d+AcYC7U7apx3wZuD5BKCFmZlz7hvn3LbA+nSgWOBqe2w6NkCcplgTybuSJeHjj6FlS+jdG558EoLowSMiIiIiEi+CKdArApuPe70lsC7HfZxzh4G9QLmT9rkJ+MY5d+jkH2BmvcxsgZktyMjICDZ75H3zjXfvuY/zXIrEtOLFvXnSu3WDgQPhnnvgyBG/U4nEJTMbY2Y7zGz5cevKmtl0M1sTeCzjZ0YREZF4E0yBbjmsO/my12/uY2Z18Lq9987pBzjnXnHONXTONSxfvnwQkXyyaJFXnBeN3U4AIr4rUgRefx0eewxeegk6dYI8zmUrIiHxBtDypHX9gBnOuerAjMBrERERiZBgCvQtQOXjXlcCtuW2j5klAKWBXYHXlYAPgDucc+vyG9g3znkFuu4/F8k/Mxg8GJ5/3pt+7ZprYOdOv1OJxBXn3OcEztXHOf6WtTeB9hENJSIiEueCKdC/BqqbWZqZJQKdgSkn7TMFbxA4gJuBmc45Z2YpwFSgv3Puy1CF9sO5+zMgI0MFukgoPfCAN/XaggXQqBGsWOF3IpF4V8E59z1A4PHsnHaKmVvTJGbMnj0bM+Oxxx7LXteoUaPsacLiyciRI6lQoQJly5alf//+eT5+zpw5mNkJ09Pt3LmTG2+8kTJlynDuuecyaNCg7G19+vQhOTk5eylatGj2/OkiEnmnLdAD95TfC0wDVgLjnXPpZvakmbUN7DYaKGdma4GH+LVL3L1ANeAJM1scWHI82Ue7+ttWe08aN/Y3iEhB06kTzJ4NP/8Ml1wC//2v34lE5DRi5tY0iSnly5fn008/BWDTpk3s3bvX50SRN3/+fAYNGsTMmTNZtmwZY8eOZfz48UEff/jwYfr160etWrVOWP/444+TlZXFtm3b+N///scrr7zC1KlTAXj55ZfJzMzMXm6++WZuuummkH4uEQleMFfQcc597Jyr4Zyr6px7OrDuz865KYHnB51ztzjnqjnnGjnn1gfW/9U5V8I5V/+4ZUf4Pk741N+2yrv3/IIL/I4iUvA0bgxffw1padC6NYwcqRHeRfzxg5mdAxB4jMlztsSmxMREateuzTfffMOECRPo0KHDCdsnTpxI3bp1KVOmDK1ateL777/P3ta5c+fsq84tWrRgzZo12dtSU1MZOHAg1apVo2zZsrz88ssR+0x59d5779GhQwfq1KlDxYoV6dmzJ2PHjg36+H/84x+0bt2aChUqnLD+u+++o02bNiQlJVGlShUuueQSVuTQa23v3r1MmjSJ22+/Pd+fRUTOTFAFusCF36/2urcnJvodRaRgqlwZvvgC2raFP/4RevSAAwf8TiUSb46/Za0bMNnHLBKHbr75ZiZMmMCkSZNo1+7XWX3nz59Pjx49GD16NBkZGTRo0IBevXplb2/QoAHLli0jIyODhg0b0qlTpxPed968eSxdupTRo0fz6KOPcvjw4Yh9prxYtWoVtWrVYsSIEfTv35/atWuzatWqoI7dvn07r7/+Og899NAp2+6++24++eQTfvrpJzZu3MjChQu55pprTtlv7NixNGjQgKpVq+b7s4jImVGBHozDh6n3w1rvHlkRCZ/kZG/QuCee8EZ6v/RSWL/e71QiBZKZvQt8BdQ0sy1m1gMYDFxjZmuAawKvRSLmuuuu4/3336do0aKcddZZ2etHjx5N9+7dady4MQkJCfTt25epU6dy6JA3e+9jjz3G2WefTeHChenevTtLliw54X3/7//+j+LFi9O6dWv279/P9u3bI/q5gvXTTz+RnJzM5s2bWbNmDSVLliQzMzOoY/v27cuAAQMoVqzYKdsaNGjA3r17KV26NGlpafTs2ZP69eufst/rr7/OHXfcke/PISJnTgV6MJYvp3jWId1/LhIJhQrBk0/C1Knw3Xdw8cXw4Yd+pxIpcJxzXZxz5zjnijjnKjnnRjvnfnTOtXDOVQ88njzKu0hYFStWjLZt23LnnXeesH7z5s2MGjWKlJQUUlJSSEtLIzExke+//54jR47Qv39/qlatSkpKCo0bN+bo0aMcOXIk+/iyZcsCXjd6gINROr1niRIlyMzMZOjQoUyYMIH9+/eTnJx82uO+/PJL1q9fT+fOnXPc3rlzZy688EJ++uknvvvuO955551T7m1fuXIlS5YsoWPHjiH5LCJyZlSgB2P+fO9RBbpI5LRqBQsXwnnned3e//QniNIuiSIiEjpDhgyha9euJ6yrXLkyTzzxBHv27MleDh48SGpqKu+88w4ffPABs2bNYs+ePXzxxRcAuBgcy6RGjRp8++232a9XrFhBzZo1T3vc119/zVdffYWZYWZ89tln3HXXXTzwwAMALFq0iO7du1O0aFGqVKlC69atmTFjxgnv8frrr3PDDTeQkpIS2g8lInmiAj0Y8+bxY1IpbwArEYmctDT48kvo2ROeeQaaNoUNG/xOJSIiEXbnnXfy8ssvs2jRIpxz7Nixg3HjxgGwf/9+kpKSKFOmDPv37+eZZ57xOe2Zu+WWW5g4cSLp6els3bqV0aNHn3I//cSJE6lWrRpbt27NXvfAAw/gnMtemjZtyquvvsrw4cMBr4v722+/TVZWFtu3b2fatGnUrl07+/gjR47w9ttvq3u7SBRQgR6MefNYck4NiMO5OEV8V6wYvPoqvPMOLF8O9evDv//tdyoREYmgJk2aMGzYMO68805Kly5Nw4YNWbhwIQB33HEHVapU4dxzz+WCCy6gSZMmPqc9c40bN2bgwIE0b96cunXr0rFjx1MK9H379rFu3TqysrKCft833niDxYsXU758eerXr0/Tpk255557srf/5z//4ciRI7Rs2TJkn0VEzoxFW/efhg0bugULFvgd41f790Pp0gy7rCsPzfmX32lE4tvGjXD77d5o7127wj//CaVLn7Jbar+pbBzcOvL5RHxkZgudcw0j/XPDcd7Wv2EREYkm4Tgv5Xbe1hX001mwAJzzrqCLiL9SU2HWLG8QuXHj4IIL4JNP/E4lIiIiIhISKtBPZ+5cABarQBeJDgkJ3jRsX3wBxYvDddd596jv2eN3MhERERGRfFGBfjpz5kDduuxNKul3EhE5XpMm8M030K+fN2d63bre1GwiIiIiIjFKBfpvOXLEu4J+xRV+JxGRnBQrBs8+C/PmQZky0KYNdOnC2ft/9DuZiIiIiEieqUD/LUuWeIPEqUAXiW4NG3rjRQwaBB98wMzX+sCwYZCHEW5FRERERPymAv23zJnjPapAF4l+RYvCwIGQns78SnXg4Yfhoovg88/9TiYiIkFITU0lKSmJ5ORkqlatygsvvJC9zcxo3Lhx9utHH30UM2P27NkATJkyhXr16lGyZEl+//vf8/TTT+f4vseWTz/9NGKfK5Rmz55NzZo1KVGiBO3bt2fv3r1BHTdkyBCqV69OyZIlqVu3LpMnT87e5pyjb9++VKhQgXLlytGlSxf27duXvX3UqFGkpqZSsmRJ2rdvz+7du0P+uUTkVyrQf8ucOd6o0ZUq+Z1ERIJVtSr/d/NAmDTJ6wHTtCl06QIbNvidTERETuPDDz8kMzOTd999l8cee4zPPvsse9uePXvYtGkTADNmzKB8+fIAbNiwgTvuuIMXX3yRffv2MWfOHKpVq5bj+x5brr766sh9qBD5+eefueWWWxg4cCA7duzAzOjfv39QxyYkJPD++++zd+9eRo0axe2338769esBGDduHOPHj2fRokVs2rSJXbt28dRTTwGwePFiHn30UaZOncoPP/xAVlYWjzzySNg+o4ioQM+dc16BrqvnIrHHDNq1gxUrvBHfJ0+GWrWgb1/QX/5FRKJeo0aNqFOnDl9//XX2ug4dOvD++++zaNEiateuTWJiIgALFy4kLS2NK6+8EjOjSpUqdOrUya/oYTNr1ixKly5N165dKVGiBH379mXcuHFBHfvQQw9xwQUXUKhQIS677DLOO+88Fi5cCMB3333HpZdeSsWKFSlRogRt2rRhxYoVAHz++ec0bdqUOnXqULx4ce677z4mTpwYts8oIirQc7dmDezYoQJdJJYVL+7Nmb56Ndx6q3dfetWq3uOhQ36nExGRXCxcuJD09HTOP//87HXt27dn8uTJTJgwgZtuuil7fb169UhPT2fgwIEsWbIE55wfkcNu1apV1KpViy+++IKWLVtSs2ZNdu3aRUZGRp7eZ/fu3axevZq6desC0LlzZzZs2MDmzZvJzMzko48+onXr1gCntKVzjt27HJyrzwAAEVZJREFUd7Nr167QfCgROYUK9Nwcu29VBbpI7KtUCcaMgcWLoVEj7/70GjVg1Cj45Re/04mISED79u0pXrw4V155JYMHD84uFAHKlStH0aJFmTBhAtddd132+po1a/Lxxx8zd+5cGjduTFpa2gn3WB9735SUlOxlx44dEftMofLTTz+RnJxMRkYGq1atolixYgBkZmbm6X169+5Nt27dsv/48bvf/Y7GjRvz+9//ntKlS1O4cGF69eoFQNOmTZk9ezbLli0jMzOTl19+GfC624tIeKhAz82nn8K550LNmn4nEZFQueAC+O9/Yfp07993nz5QvTq88ooKdRGRKDBp0iT279/PAw88wNy5c0/Zfuedd9K2bVuSkpJOWH/11Vczffp0du/ezaOPPkqXLl348cdfp9ycNGkSe/bsyV7OPvvssH+WUCtRogSZmZnceOONbNiwgazATCXJyclBv8eAAQP48ccfGTFiRPa6p556iiVLlrBjxw727t1LqVKluP/++wGoX78+Q4YMoX379tSsWTN7oL6SJUuG8JOJyPFUoOfk6FGYMQOuvtq7l1VECparr4a5c71i/ZxzoHdv74r6yy/DgQN+pxMRiWuFCxdm0KBBfPXVV9mjtB/TuXNnhg4dmuuxSUlJ3H333RQrVowNBWxw0Bo1avDtt99mv16xYgVly5bNHizvdIYPH860adOYNGlS9v37AIsWLeKmm27irLPOIjk5mdtvv50ZM2Zkb+/Tpw/r1q1j69atXHDBBaSmplK6dOnQfTAROUGC3wGi0uLFsHMnXHON30lEJFzM4Lrr4NprvUJ90CD4wx/gz3+Ge++lVVY9VmQlnvZtfkvFlCS+7HdVaPKKiMSRIkWKcN999/Hkk0/SrFmz39x3yZIlLFy4kA4dOlCqVCnGjRtHVlYW1atXj0zYCLnqqqvYu3cv77zzDu3atWPo0KF07NjxhH22bt1K06ZNGTJkCB06dMhe/9Zbb/HSSy8xZ86cU65+N2jQgA8++IDbb7+dpKQkxo8fT+3atbO3L1u2jDp16rB+/Xr69+9Pnz59wvtBReKcCvScTJ/uPbZo4W8OEQk/M7j+emjZEj77DJ57DgYO5P2EoiT17gkPPugNLHcGUvtNDXFYEZH40aNHj+wr6b8lOTmZcePG8cgjj5CVlUWNGjWYOHHiCVd5b7jhBgoXLpz9+sUXX6Rbt25hyx4OxYsX57333qNXr1707NmTa665hsGDB5+wT1ZWFuvWrTthHnOAgQMHsm3bNs4777zsdQMGDMhe7rvvPmrWrMnhw4e57LLLeOmll7L3+9Of/sTMmTMpXrw4PXr00DRrImGmAj0n06dD3bpe11cRiQ9m0KyZt6SnM6XLg3R69VX45z+hVSu4+27vivtxv+CJiEjobNy48YTXKSkp7N+/Hzh1NHGALVu2ZD+fNm1a0O8by5o1a8bq1atz3Z6amppjW/1Wd/8SJUowZsyYXLdPmTIlbyFFJF90D/rJDhyAL75Q93aReFanDo+1+iNs3AiPPw4LF0Lr1lCtGgwe7E3BKCIiIiISYirQTzZjhjc/csuWficREb+dc443j/qmTTB+PKSmQv/+3rRtHTvC1Klw+LDfKUVERESkgFCBfrLJk6FUKa+bq4gIQJEicMstMGsWrFjhDSY3cya0aQMVK8JDD8GSJX6nFBEREZEYpwL9eEePwocfegNGJeZv9GYRKaDOPx9GjIBt22DSJLjsMnjhBahf31uefRbWrvU7pYiIiIjEIA0Sd7z58+GHH6BtW7+TiEi0S0yEdu285ccfYexYePttGDDAWy68kP4VLqZ5r61sKFvxjH9MKKZqu2zwTLbuyd/87poyTkRERCT8VKAfb/JkSEjwRmwWEQlWuXJwzz3esmkTvP8+TJhA70/G0JsxUK+e94e/1q2hUaM8jQQfiqnatu45wMbBrfP1HpoyTkRERCT81MX9GOdg4kRo2hRSUvxOIyKxqkoVb+70L7+EzZu97vApKd7o75deChUqwG23wbvvwq5dfqcVERERkSiiAv2YBQtg9Wro0sXvJCJSUFSqBPffD59/DhkZXlHeqhVMmwZdu0L58l7R/vjj3qBzB/LXDV1EREREYpsK9GP+/W/vntKbbvI7iYgURGXKQOfO8NZbsH07fPWVd6+6c97V9RYtvH2aN/emdvviC2/KRxERERGJG7oHHbx5jMeO9aZMUvd2EQm3woWhSRNveeop2LfPK8hnzvSWQYNg4EBITOT9s86DrJlwySXeUvHMB5wTERERkeimAh3go4+80dtvv93vJCISj0qV8rq+Hxugctcu+OwzmDsX9+5UePFFGDbM21alileo/7//Bw0aeEuZMv5lFxEREZGQiZsC/bemGfrX2EGklSzPlXMLceR/OY9UXDElKZzxRCSEKqYk5XvUcV//zZctCzfeCDfeyM2Fm7HxyWvgm2+8bvFffQVz58K4cb/un5r6a7HeoIE3H3vFimAWskj5bVNN0yYiIiJyenFToOc6zdDKlfC3JfDMM6zrr/nPRQqCAlcIJiZC48be8sAD3rqMDK9oP7YsWgQffPDrMaVKQe3aULs2PTcA/ykEdepA5cpnVLjnt001TZuIiIjI6cVNgZ6roUOhaFHo2dPvJCIiwStfHq691luO2b8flizxlpUrYcUKmDqVx3/4AWaN8fZJToZataBqVW+pVu3X5+ecA4U0dqiIiIiIX+K7QF+9Gt58E+67z/tlV0QkypxZ1/JUSE6FRtdDIzi/yC/857ryXsG+YoVXvC9YABMmwJEjvx5WrBicd55XtKelwe9//+ty/vlQvHgIP5mIiIiInCy+C/Q//9m7et6vn99JRERyFNLu+pdffuLrrCzYtAnWrft1WbvWe/z0U/j551/3nT4drr46dFlERERE5BTxW6BPm+YNsvTnP0OFCn6nERGJvCJFfu3efjLnvNHkv/vOWy66KPL5REREROJMfBbo27dDjx5Qsyb07+93GhGR6GMG5cp5i4pzERERkYgIqkA3s5bACKAw8JpzbvBJ24sCbwEXAz8CnZxzGwPb+gM9gCPA/c65aSFLfyZ+/BHatoXdu+HDD717LkVEJKxCMfVdqHIUuFH+w+B0530REREJj9MW6GZWGHgRuAbYAnxtZlOccyuO260HsNs5V83MOgN/AzqZWW2gM1AHOBf41MxqOOeOEGHFfzkA48d795tv2wbvvefNFywiImEXLUVxNPyRINoFed4XERGRMAhmPp1GwFrn3Hrn3C/AWKDdSfu0A94MPJ8AtDAzC6wf65w75JzbAKwNvF/kzJkDJUqw4vlboFMnb1C4mTPhhhsiGkNERCRGBHPeFxERkTAIpot7RWDzca+3AI1z28c5d9jM9gLlAuv/d9KxFU/+AWbWC+gVeJlpZquCSp83ZwE7+fZbuOyyMLx93PHaU0JF7Rlaas/QKzBtan/zOwEQnvb8fYje57Tn/Qict8+yvxWM71sUKTD/hqOE2jO01J6hpfYMvXCcl3I8bwdToFsO61yQ+wRzLM65V4BXgshyxsxsgXOuYTh/RjxRe4aW2jO01J6hpzYNrShvz9Oeu8N93o7y9olJatPQUnuGltoztNSeoRfJNg2mi/sWoPJxrysB23Lbx8wSgNLAriCPFRERkeihc7eIiIhPginQvwaqm1mamSXiDfo25aR9pgDdAs9vBmY651xgfWczK2pmaUB1YH5ooouIiEgYBHPeFxERkTA4bRf3wD3l9wLT8KZbGeOcSzezJ4EFzrkpwGjgbTNbi3flvHPg2HQzGw+sAA4D9/gxgntAWLvQxyG1Z2ipPUNL7Rl6atPQitr2zO28H+EYUds+MUxtGlpqz9BSe4aW2jP0Itam5l3oFhERERERERE/BdPFXURERERERETCTAW6iIiIiIiISBQo8AW6mbU0s1VmttbM+vmdJ9aZ2UYzW2Zmi81sgd95YpGZjTGzHWa2/Lh1Zc1supmtCTyW8TNjLMmlPQeZ2dbA93SxmbXyM2MsMbPKZjbLzFaaWbqZ/TGwXt/RM/Ab7anvKKc/RwcGmR0X2D7PzFIjnzJ2BNGeD5nZCjNbamYzzCzHOXjlV8H+HmlmN5uZMzNNbfUbgmlPM+sY+J6mm9k7kc4YS4L4N18lcA76JvDvPi7PNcHK6XfKk7abmY0MtPdSM7soLDkK8j3oZlYYWA1cgzdtzNdAF+fcCl+DxTAz2wg0dM7t9DtLrDKzK4FM4C3nXN3AuiHALufc4MD/YMs45x7zM2esyKU9BwGZzrmhfmaLRWZ2DnCOc26RmZUEFgLtge7oO5pnv9GeHYnz72gw52gzuxu4wDnXx8w6Azc65zr5EjjKBdmezYF5zrmfzewPQDO1Z+6C/T0y8G97KpAI3Ouc0wWMHAT5Ha0OjAeucs7tNrOznXM7fAkc5YJsz1eAb5xzL5lZbeBj51yqH3ljQU6/U560vRVwH9AKaAyMcM41DnWOgn4FvRGw1jm33jn3CzAWaOdzJolzzrnP8WY7OF474M3A8zfxfoGXIOTSnnKGnHPfO+cWBZ7vB1YCFdF39Iz8RntKcOfo4793E4AWZmYRzBhLTtuezrlZzrmfAy//hzfHveQu2N8jnwKGAAcjGS4GBdOedwEvOud2A6g4/03BtKcDSgWelwa2RTBfzAnid8p2eMW7c879D0gJ/CE+pAp6gV4R2Hzc6y3oF6P8csAnZrbQzHr5HaYAqeCc+x68X+iBs33OUxDcG+h+NEbdsc9MoDtxA2Ae+o7m20ntCfqOBnOOzt7HOXcY2AuUi0i62JPX33l6AP8Ja6LYd9o2NbMGQGXn3EeRDBajgvmO1gBqmNmXZvY/M2sZsXSxJ5j2HATcZmZbgI/xrv7KmYtIbVnQC/Sc/specPv0R8ZlzrmLgOuBewJdQUSizUtAVaA+8D3wd3/jxB4zSwbeBx5wzu3zO0+sy6E99R0N7hyt83jwgm4rM7sNaAg8F9ZEse8329TMCgHPAw9HLFFsC+Y7mgBUB5oBXYDXzCwlzLliVTDt2QV4wzlXCa9b9tuB762cmYickwr6f6AtQOXjXldCXTvyxTm3LfC4A/gAr3uN5N8Px7rIBB7VpSsfnHM/OOeOOOeOAq+i72memFkRvGLy3865iYHV+o6eoZzaU99RILhzdPY+ZpaA10VTt7TkLKjfeczsauBPQFvn3KEIZYtVp2vTkkBdYHZgjJ4mwBQNFJerYP/NT3bOZTnnNgCr8Ap2OVUw7dkD755+nHNfAcWAsyKSrmCKSG1Z0Av0r4HqZpZmZolAZ2CKz5lilpmVCAyEgpmVAK4FchzlUPJsCtAt8LwbMNnHLDHvpPuBbkTf06AF7u8dDax0zg07bpO+o2cgt/bUdxQI7hx9/PfuZmCmK8ij2+bPadsz0B17FF5xrj+ynd5vtqlzbq9z7iznXGpg4K3/4bWtBonLWTD/5icBzQHM7Cy8Lu/rI5oydgTTnpuAFgBmdj5egZ4R0ZQFyxTgjsBo7k2Avcdu/wulhFC/YTRxzh02s3uBaUBhYIxzLt3nWLGsAvBBYHyeBOAd59x//Y0Ue8zsXbyuW2cF7gkaCAwGxptZD7z/md7iX8LYkkt7NjOz+njdjjYCvX0LGHsuA24HlpnZ4sC6Aeg7eqZya88u8f4dze0cbWZPAgucc1Pw/rjxtpmtxbty3tm/xNEtyPZ8DkgG3gucyzc559r6FjrKBdmmEqQg23MacK2ZrQCOAI845370L3X0CrI9HwZeNbMH8c433fVHztzl8jtlEQDn3Mt49/G3AtYCPwN3hiWH/huJiIiIiIiI+K+gd3EXERERERERiQkq0EVERERERESigAp0ERERERERkSigAl1EREREREQkCqhAFxEREREREYkCKtBFREREREREooAKdBEREREREZEo8P8BTGeT1Q8Ze2gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x720 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if (Nexp > 1) :\n",
    "    \n",
    "    fig2, ax2 = plt.subplots(nrows=2, ncols=2, figsize=(14,10))\n",
    "    ax2 = ax2.flatten()\n",
    "\n",
    "    # Distribution of alpha0 values:\n",
    "    ax2[0].hist(array_alpha0, bins=60, range=(0.0, 6.0), histtype='step')\n",
    "    ax2[0].set_title('Histogram of alpha0')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMSE {:>10.3f}\".format(len(array_alpha0), array_alpha0.mean(), array_alpha0.std(ddof=1))\n",
    "    ax2[0].text(0.05, 0.95, string, family='monospace', transform=ax2[0].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution of alpha1 values:\n",
    "    ax2[1].hist(array_alpha1, bins=50, range=(-0.2, 0.80), histtype='step')\n",
    "    ax2[1].set_title('Histogram of alpha1')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMSE {:>10.3f}\".format(len(array_alpha1), array_alpha1.mean(), array_alpha1.std(ddof=1))\n",
    "    ax2[1].text(0.05, 0.95, string, family='monospace', transform=ax2[1].transAxes, fontsize=12, verticalalignment='top')\n",
    "\n",
    "    # Distribution (normed) of Chi2 values, along with Chi2-distribution (for Ndof=7) drawn on top:\n",
    "    ax2[2].hist(array_Chi2, bins=25, range=(0.0, 25.0), histtype='step', density=True)\n",
    "    ax2[2].set_title('Histogram of Chi2')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMSE {:>10.3f}\".format(len(array_Chi2), array_Chi2.mean(), array_Chi2.std(ddof=1))\n",
    "    ax2[2].text(0.62, 0.95, string, family='monospace', transform=ax2[2].transAxes, fontsize=12, verticalalignment='top')\n",
    "    x_axis = np.linspace(0, 20, 1000)\n",
    "    y_axis = stats.chi2.pdf(x_axis, df=7)\n",
    "    ax2[2].plot(x_axis, y_axis, 'r-')\n",
    "    ax2[2].text(0.45, 0.65, \"Note: This is not a fit!\", transform=ax2[2].transAxes, fontsize=16, verticalalignment='top')\n",
    "\n",
    "    # Distribution of Chi2 fit probabilities:\n",
    "    ax2[3].hist(array_Prob, bins=20, range=(0.0, 1.0), histtype='step')\n",
    "    ax2[3].set_title('Histogram of Prob')\n",
    "    string = \" Entries {:>6} \\n Mean {:>9.3f} \\n RMSE {:>10.3f}\".format(len(array_Prob), array_Prob.mean(), array_Prob.std(ddof=1))\n",
    "    ax2[3].text(0.3, 0.2, string, family='monospace', transform=ax2[3].transAxes, fontsize=12, verticalalignment='top')\n",
    "    \n",
    "    plt.tight_layout()\n",
    "    if (save_plots) :\n",
    "        fig2.savefig('FitResultDistributions.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Remark on $\\chi^2$ rule of thumb:\n",
    "A general rule is that the \"reduced chisquare\" (defined as $\\chi^2$/Ndof) should be \"around 1\". That is true, but the rule definitely has its limitations, since \"around\" depends on the size of Ndof! This is examplified below, considering two cases with the same reduced chisquare of 1.5, but with **very** different conclusions.\n",
    "\n",
    "The reason for this can be found in the mean and width of the $\\chi^2$ distribution (with Ndof = k), which are k and sqrt(2k). From the central limit theorem, the $\\chi^2$ should approach a Gaussian for \"large\" values of k (k > 50 is probably enough), but the **relative** width (i.e. width/mean) becomes smaller and smaller as k grows, and a reduced chisquare of e.g. 1.5 becomes more and more unlikely, as k increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "lines_to_next_cell": 2,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:\n",
      "  Low number of degrees of freedom:   Prob(  3.0,  2) = 0.2231302\n",
      "  High number of degrees of freedom:  Prob(300.0,200) = 0.0000059\n"
     ]
    }
   ],
   "source": [
    "# Numbers to illustrate the point, that using the ratio Chi2/Ndof (\"reduced chisquare\") is only a rule of thumb:\n",
    "print(\"\\n  Considering the relation between Chi2/Ndof ratio and Prob(Chi2,Ndof) vs. Ndof:\")\n",
    "print(\"  Low number of degrees of freedom:   Prob(  3.0,  2) = %9.7f\"%(stats.chi2.sf(  3.0,   2)))\n",
    "print(\"  High number of degrees of freedom:  Prob(300.0,200) = %9.7f\"%(stats.chi2.sf(300.0, 200)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The conclusion is, that while the reduced $\\chi^2$ is a good rule of thumb, it is NOT the correct way to consider a $\\chi^2$ value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "\n",
    "\n",
    "# Questions:\n",
    "\n",
    " 1) Run the code such that it does exactly one fit (initial setting), and take a look at the line fitted.\n",
    "    Does this look reasonable, and are you satisfied with the fit probability? Does the fit reproduce the input\n",
    "    numbers well (include the uncertainties in your answer)?\n",
    "\n",
    " 2) What is the chance that the input for the data could actually be from a flat distribution, i.e. without the slope?\n",
    "    [There is a correct (slower) way of doing this, and another quick but approximate way of testing this].\n",
    "\n",
    " 3) Now increase the number of experiments to e.g. 1000, and rerun the macro. Figure\n",
    "    out what you see in the window split in 2-by-2, and go through each of these to\n",
    "    see, if you understand every feature of the distributions shown, and if you are\n",
    "    happy with them! Is this what you expected? Are they also as wide as you would expect?\n",
    "    Possibly try to make some fits of these output distributions to check, if you're satisfied!\n",
    "    <br>\n",
    "    This somehow makes this the \"long\" question without any leading questions.\n",
    "    But you should try to be statistically minded enough and have an idea of what to\n",
    "    look for, at least to some degree. We will discuss this further in class.\n",
    "\n",
    " 4) Investigate if the distributions of probabilities is flat, or if it has some\n",
    "    slope to it. Do you understand why the distributions of probabilities is flat?\n",
    "    (This is typically conceptually hard to begin with, but don't worry!)\n",
    "\n",
    " 5) Find the line of code where the random points for the line are generated, and add a\n",
    "    quadratic term in x with the coefficient alpha2 = -0.04. Run the program again with Nexp=1,\n",
    "    and start by looking at the single fit in the graph. Can you see this change?\n",
    "    Now run 1000 experiments again, and see what has changed in the distributions\n",
    "    in the 2-by-2 window when including such a term. Are these changes as you would\n",
    "    have expected them?\n",
    "    <br>\n",
    "    NOTE: The reason why the quadratic term is more clear in the latter case is NOT\n",
    "    that it is a better method. It is simply because there is 1000 times more statistics!\n",
    "\n",
    " 6) With the quadratic coefficient from question 5), the linear fit of course does\n",
    "    not do very well. What changes are needed for the fit to be good again? Make\n",
    "    these changes, and see that the condition in question 4) is again met.\n",
    "\n",
    "### (Semi-)Advanced questions:\n",
    "\n",
    " 7) Try to generate data according to more advanced functions, and then fit the data with these.\n",
    "    Do you still obtain a flat distribution in probability? Do you also manage to fit these with\n",
    "    simpler/other functions and still maintain the flat distribution?\n"
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}