{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Error Propagation using Random Gaussian Numbers\n",
    "\n",
    "### Authors: \n",
    "- Christian Michelsen (Niels Bohr Institute)\n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 06-11-2018 (latest update)\n",
    "\n",
    "***\n",
    "\n",
    "The example is based on FIRST doing the error propagation analytically, and then verifying it by running a so-called Monte-Carlo (MC) program.\n",
    "\n",
    "For more information on error propagation, see:\n",
    "-     R. J. Barlow: page 48-61 \n",
    "-     P. R. Bevington: page 36-48\n",
    "\n",
    "***\n",
    "\n",
    "DO THE FOLLOWING ANALYTICAL EXERCISE FIRST!!!\n",
    "\n",
    "1. A class of students estimate by eye, that the length of the table in Auditorium A is $L = (3.5\\pm 0.4)$m, and that the width is $W = (0.8\\pm 0.2)$m.\n",
    "\n",
    "   Assuming that there is no correlation between these two measurements, calculate ANALYTICALLY what the circumference (C), area (A), and diagonal (D) length is including (propagated) uncertainties. Repeat the calculation, given that the correlation between length and width is $\\rho(L,W) = 0.5$ - not an unreasonable number, given that they are estimated by the same (uncertain) scale."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# NOTE: Do the above analytical calculation before you continue below!\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "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",
    "from probfit import BinnedLH, Chi2Regression, Extended # Helper tool for fitting\n",
    "import sys                                             # Modules to see files and folders in directories"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../External_Functions')\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Error propagation - Simulation\n",
    "\n",
    "First we set the parameters of the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_exp = 10000           # Number of experiments\n",
    "pi = np.pi\n",
    "save_plots = False\n",
    "r = np.random\n",
    "r.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define parameters for two random numbers (Gaussianly distributed):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "mu1   =  3.5\n",
    "sig1  =  0.4\n",
    "mu2   =  0.8\n",
    "sig2  =  0.2\n",
    "rho12 =  0.5           # Correlation parameter!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "if not (-1.0 <= rho12 <= 1.0): \n",
    "    raise ValueError(f\"Correlation factor not in interval [-1,1], as it is {rho12:6.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we calculate numbers that allows the transform from uncorrelated variables `u` and `v` to correlated random numbers `x1` and `x2` below (see Barlow page 42-44 for method):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "theta = 0.5 * np.arctan( 2.0 * rho12 * sig1 * sig2 / ( np.square(sig1) - np.square(sig2) ) )\n",
    "sigu = np.sqrt( np.abs( ((sig1*np.cos(theta)**2) - (sig2*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n",
    "sigv = np.sqrt( np.abs( ((sig2*np.cos(theta)**2) - (sig1*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n",
    "\n",
    "sigu = np.sqrt( np.abs( (((sig1*np.cos(theta))**2) - (sig2*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n",
    "sigv = np.sqrt( np.abs( (((sig2*np.cos(theta))**2) - (sig1*np.sin(theta))**2 ) / ( (np.cos(theta))**2 - np.sin(theta)**2) ) )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the above is nothing more than a matrix multiplication written out! Also note that the absolute value is taken before the square root to avoid `np.sqrt(x)` with `x<0`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop over process `Nexp` times: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Gaussian:  x1 = 3.78    x2 = 0.83   ->   y = 3.15\n",
      "Gaussian:  x1 = 3.48    x2 = 0.76   ->   y = 2.66\n",
      "Gaussian:  x1 = 3.69    x2 = 0.91   ->   y = 3.35\n",
      "Gaussian:  x1 = 4.28    x2 = 0.97   ->   y = 4.14\n",
      "Gaussian:  x1 = 4.30    x2 = 0.73   ->   y = 3.12\n"
     ]
    }
   ],
   "source": [
    "x1_all = np.zeros(N_exp)\n",
    "x2_all = np.zeros_like(x1_all)\n",
    "x12_all = np.zeros((N_exp, 2))\n",
    "y_all = np.zeros_like(x1_all)\n",
    "\n",
    "for i in range(N_exp): \n",
    "\n",
    "    # Get (uncorrelated) random Gaussian numbers u and v:\n",
    "    u = r.normal(0.0, sigu)\n",
    "    v = r.normal(0.0, sigv)\n",
    "\n",
    "    # Transform into (possibly) correlated random Gaussian numbers x1 and x2:\n",
    "    x1 = mu1 + np.cos(theta)*u - np.sin(theta)*v\n",
    "    x2 = mu2 + np.sin(theta)*u + np.cos(theta)*v\n",
    "    \n",
    "    x1_all[i] = x1\n",
    "    x2_all[i] = x2\n",
    "    x12_all[i, :] = x1, x2\n",
    "    \n",
    "    # Calculate y (circumference, area or diagonal):\n",
    "    y = x1*x2                  # Just nonsense formula! Write the appropriate formula yourself!\n",
    "    y_all[i] = y\n",
    "\n",
    "    if (i < 5) :\n",
    "        print(f\"Gaussian:  x1 = {x1:4.2f}    x2 = {x2:4.2f}   ->   y ={y:5.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is written the alternative Numpy version of the above cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = r.normal(0.0, sigu, N_exp)\n",
    "v = r.normal(0.0, sigv, N_exp)\n",
    "x1_all = mu1 + np.cos(theta)*u - np.sin(theta)*v\n",
    "x2_all = mu2 + np.sin(theta)*u + np.cos(theta)*v\n",
    "x12_all = np.array([x1_all, x2_all])\n",
    "\n",
    "#y_all = x1_all - x2_all"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "Plot both input distribution and resulting 2D-histogram on screen:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAApQAAAGiCAYAAABZBFYRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl8VNX9//HXJyEJSxYQKCCLoKCALF8wUVoXEExJAEXBBbAVqIKyqBW1ol8riIhotUJFacH6CyogEYXSYlEWA1ZRFsu3CrgAggiIC0jYyXJ+f8wwnZBlJskkE4b3s495yL33zDnn3mDz8Syfa845RERERETKKircHRARERGR05sCShEREREpFwWUIiIiIlIuCihFREREpFwUUIqIiIhIuSigFBEREZFyUUApIiIiIuWigFJEREREykUBpYiIiIiUiwJKERERESmXauHugIiIiMiZpueVtdyP+/JCUtf6/xx/2zmXFpLKykgBpYiIiEgl+3FfHmvebhaSuqIbfVkvJBWVgwJKERERkUrmgHzyw92NkFFAKSIiIlLpHHkucgJKbcoRERERkXJRQClSSma20cy6hbsf4WRm15nZTjM7ZGadirjuzKxlGPrVzcy+CVFd9c3sMzOrEUTZZ8xsRCjaDYaZZZnZbeWs489m9vsSrpf4MzSz7WZ2VXn64FfXXDO7NohyHczsg1C0KRJunilvF5JPVaCAUsRPUb8kzWyImf3r5LFz7kLnXFaAepp7fyFH6rKSp4HRzrl459y/w9WJCg5cxwIZzrmjQZR9GnjIzGIrqC8h55y7wzn3WDBlzSzDzCZWRD/MrAPQEfhboLLOuf8AP5nZ1RXRF5HKlh+i/1UFCihFTkNVIFA9B9gY5j5UGDOLAwYDrwZT3jm3B/gMuKYi+xWhbgdmO+eCHWaZ7f2OiFQhCihFSsl/FNPMLjazdWaWbWZ7zeyP3mKrvP/8yTst/HMzizKzh81sh5l9Z2Yvm1mSX723eK/9aGa/P6Wd8WY238xeNbNsYIi37dVm9pOZ7TGzaf4jZN7Ru5Fm9qWZHTSzx8zsPDP7wNvfzOJG1Irrq5nFmdkhIBr4PzPbGsTzijOzp83sa+8z+vPJaeSTU9Rmdq+3nT1mNtTvu3XN7O/e/q41s4knR4vN7OQz/j/vM77J73vF1dfLzDZ5n8cuM7uvmG5fAvzknPvG+72zvP282nscb2ZbzOwWv+9kAb1LeA6vm9m3ZnbAzFaZ2YV+1zLM7HkzW+zt20dmdp7f9VTzTL8fMLNpgBXTRnUzO2pm9bzH/2tmuWaW6D1+zMym+LU50e+793uf124z+43f+eHAzcDvvM/5735N/o+Z/cfbr3lmVt37nXpm9g/v3819ZvaemRX3+yYdWOn9Xqy3fHu/9n9mZkfMrL7fc+5hnqBf5LTlcOS50HyqAgWUIuUzFZjqnEsEzgMyveev8P6ztndaeDUwxPu5EjgXiAemAZhZW+AFPL+4GwFJQONT2uoLzAdq4xmlyQPuAeoBPwd6ACNP+U5P4CKgC/A7YAbwK6Ap0A4YWMx9FdlX59xx51y8t0xH59x5RX+9gMnA+cD/AC299/WI3/WGfvd7K/C8mdXxXnseOOwtM9j7AcA5d/IZd/Q+43lB1PdX4HbnXIL3/lcU0+f2wOd+be0DfgPMNLOfAc8CG5xzL/t9ZzOeqdvi/BNoBfwM+BjPz9DfAOBRoA6wBXgcPMEZ8CbwMJ6f9Vbg0qIacM4dA9YCXb2nugI7/Mp3xRu8+TOzNOA+INXbR9+yD+fcDG9fn/I+Z//p5huBNKAF0AHP3xmAe4FvgPpAA+AhKLzQy8xqeb/7ubetE8BreP6OnjQQWO6c+95bZheQA1xQ1DMQOZ1oDaVIZFvoHVn5ycx+whPoFScHaGlm9Zxzh5xzH5ZQ9mbgj865bc65Q8CDwADzTF9fD/zdOfcv7y/VRyj8C3i1c26hcy7fOXfUObfeOfehcy7XObcd+Av/DSROeso5l+2c2wh8Crzjbf8AngCn0IaaIPoaNDMzYDhwj3Nun3PuIDAJT/B0Ug4wwTmX45x7CzgEXGBm0UB/YJxz7ohzbhMwK4hmi6zP71pbM0t0zu13zn1cTB21gYP+J5xz7wCvA8uBXhSedj3o/V6RnHMvOecOOueOA+OBjuY3Qg0scM6tcc7l4gng/sd7vhew0Tk33zmXA0wBvi3h/lcCXb0/qw7An7zH1YEU/jt67u9G4P855z51zh329i8Yf3LO7fYG3H/363MOnv8wOsf7c3ivmCntk8/L/1nPAgZ6/+4A/Bp45ZTvlfisRaTyKaAUKexa51ztkx8Kj/r5uxXP6Ntn3inZPiWUPRvPaNFJO/Dkgm3gvbbz5AXn3BHgx1O+v9P/wMzO904rfmueafBJeEaw/O31+/PRIo7jKVpJfS2N+kBNYL1fgL7Ee/6kH71B1ElHvP2q723T/74LPINiFFcfeALUXsAOM1tpZj8vpo79QEIR52fgGdnMcM6d+vNJAH4qqjIzizazyWa21fuz2u695P/z8g8S/ft86t8NR8nPYSXQDegMfAIsxfMfGl2ALUX0u1AbFPzZl6S4Pv8BzyjrO2a2zczGFvP9k8/L96ydcx956+pmZq3xjGovOuV7xT5rkdOFA/JwIflUBQooRcrBOfelc24gnmnMJ4H53mm8ov4N341nM8tJzYBcPEHeHqDJyQveNYZ1T23ulOPpeDaCtPJOuT9EMWvryqCkvpbGD3gC1wv9gvQkv2nzknzvbbOJ37mmpWy/AOfcWudcXzw/r4X8d4nCqf6D5z8UfLwjpjOAl4GRVnh3eRvg/4qpbxCeJQtX4ZmOb36y2iC6vQe/+/aO3JX0HD7AMyJ7HbDSO7LbDE8gXWi6u6g2vOX9leo3lnck9l7n3Ll4NiqNMbMeRZQ7jGcK//xTLs3CM+39a2C+dyofADNrDMTityRB5HSlKW8RAcDMfmVm9Z1z+fx3xCQfTzCUj2f94UlzgXvMrIWZxeMZUZznHU2bD1xtZr8wz0aZ8QQONhKAbOCQdyQnlHkQS+pr0LzPZSbwrHftIWbW2Mx6BvHdPDxrB8ebWU3vPd5ySrG9FHzGxfJu+LjZzJK8U8fZUGy+jTVAbW/wctLJdYC/wTMC97I3yDypK55lBEVJAI7jGXWuied5BmsxcKGZ9fNOY9+FZ51okbyj2+uBUfw3gPwAuIPiA8pMPBu92ppZTWDcKdeDfs4AZtbHzFp6g98DeNb7Fves36LwUo1X8QTEv8ITwPvrCqzwLh0QkSpCAaVI+aQBG82z83kqMMC7vvEInk0V73unersAL+FZC7YK+Ao4BtwJ4F3jeCeeDQl78Kz7+w5PEFKc+/CMfB3EE7TNK6FsaRXb1zJ4AM/054fe6d5lBL+hYjSeEb1vvf2ZS8FnMh6Y5X3GNwZR36+B7d5+3IFnrWgh3nWsGXg3h5jZRcAY4BZvoPsknuByrPd6I6AtnlHPoryMZxp5F7AJKGmt7al9+QG4Ac/mph/xbJp5P8DXVgIxeALjk8cJFL1+EufcP/GszVyB52d16malv+JZe/qTmRV3j/5a4fk5HwJWAy84594tpuwM4Ga/NZM453bi2bjkgPdOKX8z8Ocg+iBSpTmIqF3eFnzqLxGpLN5RwZ/wTGd/Fe7+VBVm9iTQ0Dk3OGDh8rdVH08w0ylQcnMzewbY6pwraQOXFMPM5gCZzrmFfudeAnY75x72O9cB+Itzrri1ryKnjY4dY90/3zp12XvZNG6yZ71zLjkklZVRuJMji4iXeXIcLscz1f00ng0V28PZp3DzTnPH4nkWKXg2QZXrlYPB8qapaR1k2XsruDsRzTk3yP/YzJoD/TglC4H3TTkKJkWqoLBNeZsnAe8aM/s/87wb+dEiysR5k+VuMU+i3+aV31ORStMXz2aY3XimDAeU4u0hkSoBzzrKw3im9J8hiFf0yenLzB7Dk+LqDxqdl0jmQrTDu6rs8g7blLd3vUwt59whM4sB/gXc7Z/Hz8xGAh2cc3eY2QDgOufcTcVUKSIiInJa6NAhxi0K0ZR3i6bfhn3KO2wjlM7jkPcwxvs5Nbrty38TGc/H87qtUKVFEREREZEQCOsaSm/KjfV4Etc+701o668x3mS7zrlcMzuAJzffD6fUMxzP2zioVavWRa1bB7XsSURERM5A69ev/8E5Vz9wyYrjKD6X1ukorAGlN/3G/5hZbWCBmbVzzn1ahnpm4Ek9QXJyslu3bl2IeyoiIiKRwsyCfRtURfaCvJC9iyL8qkQeSufcT8C7eHL6+duF9+0N3oS+SRR+HZ2IiIiIhFE4d3nX945MnnzNXCqe18j5WwSczDd3PZ63I1SN7UwiIiIiZeSAfBeaT1UQzinvRnjecBGNJ7DNdM79w8wmAOucc4vwvJ3hFTPbAuwDBoSvuyIiIiKhE0lT3mELKL0JajsVcf4Rvz8fw/PKMRERERGpovSmHBEREZFK5tAIpYiIiIiUU76LnICySuzyFhEREZHTl0YoRURERCpZpE15a4RSRERERMrljA8ox48fj5kV+tSuXbtM9U2ZMoWsrKxy92v79u2YWUjqOtNkZmaSmppKgwYNqF27Nt26deODDz4Id7dERER8HEYeUSH5VAWa8gaSkpJYsmRJgXPVqpXt0UyZMoUhQ4bQrVu3cvWpUaNGrF69mrZt25arnjPRtGnTOPfccxk1ahQ1a9bkz3/+M1deeSXr1q2jffv24e6eiIgIEFmbchRQ4gkeu3TpEu5uFBAXF1fl+nS6WLBgAXXr1vUdd+vWjaZNm/LCCy8wffr0MPZMREQkMlWNcdIq7uSI44wZM2jatCl169ZlzJgx5Ofn+8qcnCrfsWMHjz76qO94yJAhBerKyMjAzPj8889JS0ujZs2aNGrUiL/97W8AbNiwocDUe3FT3suWLaNLly7UqFGDs88+m0mTJhUqM3PmTFq3bk2NGjVo3LgxAwYM4Pjx4yF7LlWVfzAJEBsbS4sWLdi+fXt4OiQiInKKk5tyQvGpCjRC6ZWbm1vgOCoqiqio/8bbn376KQsXLuTFF18kKyuLyZMnc+WVV3L11VcDsHr1agCuu+46evfuzW233QZA/fr1i2zvV7/6FT179uS+++7jiy++4OjRowCcf/75rF69mj179tCvX78iv5uVlUV6ejrXX389jz76KJ999hkPPvgg9evXZ9iwYQCsWrWK4cOH88ADD5CWlsbevXt5/fXXOX78OHFxceV4UqefI0eOsHnzZt+zERERCT8jz0XOuJ4CSuDHH38kJiamwLmePXsWWFd59OhR5s6dS1JSEj179mTOnDksX77cF1CenJ6Oi4ujSZMmAaere/XqxaOPPgrAVVdd5Ttfs2ZNunTpUuJo2oMPPkjnzp2ZM2cOZkbPnj3Zu3cvkydP9gVNa9asoXr16kyePNn3vZtuuimIpxF5Jk+ezIkTJxg9enS4uyIiIhKRFFDi2ZSzbNmyQuf8tW7dusC5Zs2asXfv3jK3eeONN5bpe4cPH2bNmjVMmjSJvLw83/mLL76YJ554goMHD5KQkEC7du04duwYw4cP5+abb+biiy+mRo0aZe7v6SorK4snnniCqVOn0rx583B3R0REBPBMeedH0MrDyLmTcqhWrRrJyckFPq1atSpQJiEhocBxdHQ0OTk5ZW6zcePGZfre/v37yc/PZ+zYscTExPg+1113HQC7du0CIC0tjT//+c989NFHdO/enbPOOouRI0cWCEIj3datW7n++usZOnQoI0eODHd3RERECtAaSim3sqYlqlOnDmbGhAkTSEtLK3TdfxTu9ttv5/bbb+fHH39k6tSpPPbYY3Tr1q3Mo6Onk59++ok+ffrQqVMnnn/++XB3R0REJKIpoAyxhIQEDh8+XGH116pVi5SUFL766iuSk5OD+k7dunV59NFHefLJJ30jmJEsJyeH/v37U61aNebPn19ofayIiEi4OadNOREnNzeXDz/8sND55OTkUo8ktm/fnvnz59OrVy8aNWpEUlISjRo1Cvr7mzZtIjs7mz179viOq1evDvx348+kSZNIS0sjLi7Otyno3//+N2vXrmXBggUAPP3002zbto2rrrqKOnXqMGfOHPLy8sqdcP10MGLECN5//31eeeUVNm/e7DufmJioRPEiIlJl5FeR6epQUEAJHDhwgJ///OeFzn///ffUq1evVHVNnDiRYcOG0bdvXw4ePMjgwYPJyMgI+vsjR45k5cqVvuNRo0b5/uycA6BHjx68/fbbjB8/noyMDGJjY2nXrh2DBw/2le3UqRNLliwhMzOTo0ePcsEFF/D666/TqVOnUt3P6WjZsmUcP3680NR+165d9SpLERGRCmAng5RIkZyc7NatWxfuboiIiEgVZWbrnXPBrRurIK3a13BTF50Xkrp6n7sx7PejEUoRERGRShdZaygj505EREREJCw0QikiIiJSySItsbkCShEREZEwyHORs8s7ckLjMsjIyMDM6Ny5c4Hzl19+OWbGtGnTwtSz08usWbO45JJLqFOnDomJiVx22WWsWLGi1PWMHz8eMyv0+emnnwqUW7VqFZ07d6ZGjRp06NCBpUuXFlvn5s2biYmJCTpnp4iIiJTeGR1QnrRt2za++eYbwJMqaOPGjWHu0ell7969pKWlkZGRwZtvvkmrVq1IS0tj/fr1pa6rYcOGrF69usDH/7WXO3fupE+fPrRu3Zq33nqLlJQU+vbty5dffllkfffee2+h12aKiIiEm8PIIyokn6pAU95Az549WbRoESNHjuQf//gHqampZGZmhrtbp43f/e53BY579OhBVlYWs2fP5qKLLipVXXFxcb4E7kWZPn068fHxvvybV1xxBVlZWUybNo2pU6cWKPv222/z9ddf06dPHzZt2lSqfoiIiFS0fO3yjizXXHMNixYtAmDRokVcc801hcosW7aMLl26UKNGDc4++2wmTZpU4Ponn3zCDTfcQOPGjYmLi+OCCy4o9A7pIUOG0K1bN2bMmEHTpk2pW7cuY8aMIT8/P+i+5ubmkpycTP/+/Qucv+mmm+jQoQMnTpwIuq6KYmYkJCSQnZ0d8rpXrFhBamoqsbGxAERHR5OWllZoij0vL497772Xxx9/nKgo/TUXERGpSPpNC6SmpvLRRx/x/fffs2LFCtLT0wtcz8rKIj09nRYtWrBw4UIeeOABJk6cyMyZM31ltm7dSrNmzfjTn/7EO++8w8iRIxkzZgxz584tUNenn37KwoULefHFFxk+fDjPPvssixcvDrqv1apV4+WXX2bx4sW89tprALz55pssWLCAV155xRdohUNubi779u3jmWeeYdOmTQwaNKjUdezevZs6deoQGxvLpZdeyurVqwtc37JlCy1btgTwBaznnXceW7ZsKVDuL3/5C4mJifTt27eMdyMiIlJxHGjKO9LUrFmTyy+/nPvvv59OnTpx1llnFbj+4IMP0rlzZ+bMmYOZ0bNnT/bu3cvkyZMZNmwYANdeey3XXnst4HlF4qWXXsq7777La6+9xsCBA311HT16lLlz55KUlETPnj2ZM2cOy5cv972TOxht27Zl4sSJjB49mg4dOjBixAjGjRtHx44dQ/A0yubQoUO+tYrVq1cnMzOT7t27l6qO888/n2effZZ27drx3Xff8cQTT5Camsrnn39O48aNAc9rMhMTE5k3bx4DBgwgMzOTxMREjh07xokTJ4iNjeXAgQOMGzeON954I+T3KSIiEgoOq7Rd3mbWFHgZaIAnlp3hnJtqZmcB84DmwHbgRufcfjMzYCrQCzgCDHHOfVxSG1UjrK0CrrnmGmbNmlVouvvw4cOsWbOGfv36kZeXR25uLrm5uVx88cVs27aNgwcPAnDkyBHGjh1LixYtiImJISYmhr/97W/s3bu3QH2tW7cmKSnJd9ysWbNCZYIxZswY2rRpQ0pKCs2bN2fs2LFluOvQqVmzJmvXrmXp0qUMHDiQIUOGsGbNmlLVMWjQIEaNGkXXrl254YYbWLZsGXFxcfzxj38sVDY+Pp6EhATi4+MLXZswYQIXX3wxV1xxRZnvR0REJILkAvc659oCXYBRZtYWGAssd861ApZ7jwHSgVbez3BgeqAGFFB6XX311fTs2dM3ynjS/v37yc/PZ+zYsb5AMSYmhuuuuw6AXbt2AZ6NKdOnT2fMmDGsXLmStWvXkp6eTm5uboH6Tt1xHB0dTU5OTqn7GxUVxcCBAzly5AgDBgwgOjq61HWEUlRUFMnJyVx11VW89NJLpKSk8PDDD5erzrPOOotf/OIXfPrpp75zSUlJZGdn07t3b7Kzs0lPTyc7O5vq1asTGxvLjh07eP7553nooYc4dOgQhw4dIjc3l/z8fA4dOkReXl55b1VERCQk8okKyScQ59yekyOMzrmDwGagMdAXmOUtNgs4GQT1BV52Hh8Ctc2sUUltaMrbq0GDBixZsqTQ+Tp16mBmTJgwgbS0tELXmzdvDnjWMd53333ceeedvmsVuUFm9+7dPPzww3Tt2pXx48fTv39/mjVrVmHtlVbnzp0LrR8tCzPj2LFjvuOWLVuydevWAmW2bdvmW1f51Vdfcfz4cS677LJCdSUkJLB06VKuuuqqcvdLRESkPJwjlO/yrmdm6/yOZzjnZhRV0MyaA52Aj4AGzrk93kvf4pkSB0+wudPva994z+2hGAooA6hVqxYpKSl89dVXJSbHPnr0KDVq1PAd7927l/fff58LL7ywQvp166230qZNG5YvX05qaipDhgxh+fLleJY9VC7nXKF2169fT8OGDQuVPXDgAHv27KFevXrUq1evxHoPHDjABx98wK9//Wvfue7duzNr1ixycnKIiYkhPz+fJUuW0Lt3bwA6derEe++9V6CeSZMmsW3bNl588UXat29f1tsUERGpqn5wzgV8g4eZxQNvAL91zmX7/+52zjkzc2XtgALKIEyaNIm0tDTi4uJ8m2f+/e9/s3btWhYsWAB4ci9OmTKF5s2bExUVxWOPPUbdunUrpD/Tp08nKyuLDRs2EB0dzUsvvUT79u2ZOnUqv/3tbyukzZKkpKRwww030KFDBwDmzZvHu+++S0ZGRqGyCxYsYOjQoYwbN47x48cXuHb55Zdz00030bp1a/bv38/TTz9NTk4Od999t6/MiBEjeO655xgyZAi33XYbs2fPZvfu3YwePRrwTImfOjr5s5/9jO+++67IUUsREZHwMPKpvEEgM4vBE0zOds696T2918waOef2eKe0v/Oe3wU09ft6E++5YimgDEKPHj14++23GT9+vC+hdrt27Rg8eLCvzLRp0xg+fDhDhw4lKSmJe++9l40bN7Jhw4aQ9mXr1q3cf//9PPbYY1xwwQWAZ9r9qaeeYsyYMaSlpdG6deuQthlI165dmT17No8//jhmRuvWrXnjjTfo169fqeo599xzefLJJ/n++++JiYnhsssu47333vMtKwBo2rQpixcv5u677yY9PZ1WrVqxcOFCWrVqFeK7EhERqTiOkE55l8i7a/uvwGbnnP9O10XAYGCy959/8zs/2sxeAy4BDvhNjRfdhnNlHt0sl+K2sJ9Sphuem/vKe+pN59yEkupNTk5269atK6mIiIiInMHMbH0wU8QV6Zx2Ce6hNzqHpK47Wq8q8X7M7DLgPeAT4OTbVB7Cs44yE2gG7MCTNmifNwCdBqThSRs01DlXYnAVzhHKk1vYPzazBGC9mS11zp36jrz3nHN9wtA/ERERkQpTWUnJnXP/gmLn13sUUd4Bo0rTRtjSBpWwhV1ERERETiNVYg3lKVvYT/VzM/s/YDdwn3NuYyV2TURERCTkHEZ+Jb0ppzKEPaA8dQv7KZc/Bs5xzh0ys17AQjxZ20+tYzieTO5VKhejiIiISHGqynu4QyGsd1LMFnYf51y2c+6Q989vATFmVih5oXNuhnMu2TmXXL9+/Qrvt4iIiIj8V9gCyhK2sPuXaegth5ldjKe/P4aqDxkZGZgZnTsX3GV1+eWXY2ZMmzYtVE1FvFWrVtG5c2dq1KhBhw4dWLp0aZnqmTFjBm3btqVmzZqcf/75PPvss6VuKzMzk9TUVBo0aEDt2rXp1q0bH3zwQZn6IyIiUhEckO+iQvKpCsLZi0uBXwPdzWyD99PLzO4wszu8Za4HPvWuofwTMMBVQJ6jbdu28c033wDw/fffs3GjlmmWxs6dO+nTpw+tW7fmrbfeIiUlhb59+/Lll1+Wqp7XX3+d22+/nX79+vGPf/yDW265hXvvvZeZM2eWqq1p06bRuHFj/vKXv5CZmclZZ53FlVdeySeffBKyexYRESkfIy9En6ogbGsoA2xhP1lmGp48SBWqZ8+eLFq0iJEjR/KPf/yD1NRUMjMzK7rZiDF9+nTi4+N9Sd+vuOIKsrKymDZtGlOnTg1cgdf8+fO55JJLmDhxIuB5zWJWVhYLFixg2LBhQbe1YMGCAm8p6tatG02bNuWFF15g+vTpIbxzERERgTCvoawqrrnmGhYtWgTAokWLuOaaawqVWbZsGV26dKFGjRqcffbZTJo0qcD1Tz75hBtuuIHGjRsTFxfHBRdcwPPPP1+gzJAhQ+jWrRszZsygadOm1K1blzFjxpCfn0+wsrOzqVmzZoFRO/CMspoZ//znP4OuK1RWrFhBamoqsbGxAERHR5OWlsaKFStKVU9ubi6JiYkFziUkJOA/KB1MW6e+8jI2NpYWLVqwffv2UvVHRESkomjKOwKlpqby0Ucf8f3337NixQrS09MLXM/KyiI9PZ0WLVqwcOFCHnjgASZOnFggqNu6dSvNmjXjT3/6E++88w4jR45kzJgxzJ07t0Bdn376KQsXLuTFF19k+PDhPPvssyxevDjoviYmJnLdddfx6quvFjg/Z84cGjZsyC9/+csyPIHy2bJlCy1btgQ8AS/Aeeedx5YtW0pVz29+8xtWrVrF4sWLOXjwIO+88w7vvPMOd9xxh69MWdo6cuQImzdv5sILLyxVf0RERCqSprwjTM2aNbn88su5//776dSpE2eddVaB6w8++CCdO3dmzpw5mBk9e/Zk7969TJ482TcVe+2113LttdcC4Jzj0ksv5d133+WgpytpAAAgAElEQVS1115j4MCBvrqOHj3K3LlzSUpKomfPnsyZM4fly5dz9dVXB93fwYMHk5aWxtdff+1LkzR37lwGDRpEdHR0eR9HqR04cIDExETmzZvHgAEDyMzMJDExkWPHjnHixAnfaGIgvXv35q9//Sv9+vXjxIkTREdH8/zzz9O3b99ytTV58mROnDjB6NGjQ3bPIiIi8l8aofS65pprmDVrVqHp7sOHD7NmzRr69etHXl4eubm55ObmcvHFF7Nt2zYOHjwIeEbBxo4dS4sWLYiJiSEmJoa//e1v7N27t0B9rVu3JikpyXfcrFmzQmUCueqqqzj77LOZPXs2ABs2bGDTpk3ccsstZbn1kImPjychIYH4+PgyfX/lypWMGjWKRx55hKysLB5//HHuuece3nyzUEapoNvKysriiSee4JlnnqF58+Zl6peIiEioOWea8o5EV199NT179vSNMp60f/9+8vPzGTt2rC9QjImJ4brrrgNg165dAPzud79j+vTpjBkzhpUrV7J27VrS09PJzc0tUF9CQkKB4+joaHJyckrV16ioKH71q1/5Aso5c+bQvn17OnbsWKp6QiUpKYns7Gx69+5NdnY26enpZGdnU7169aBHJwHuu+8+rrvuOv73f/+Xrl278sADDzB06FDGjh1bpra2bt3K9ddfz9ChQxk5cmTI7ldERCQU8lxUSD5Vgaa8vRo0aMCSJUsKna9Tpw5mxoQJE0hLSyt0/eSo15tvvsl9993HnXfe6bt24sSJCuvv4MGDefLJJ/n444957bXXuOuuuyqsrUBatmzJ1q1bC5zbtm2bb61jsDZu3MhNN91U4Fy7du2YPn06eXl5REdHB93WTz/9RJ8+fejUqVOhzVEiIiISWgooA6hVqxYpKSl89dVXJCcnF1vu6NGj1KhRw3e8d+9e3n///QrbCNKmTRtSUlIYPXo0u3fv5uabb66QdoLRvXt3Zs2aRU5ODjExMeTn57NkyRJ69+5dqOyBAwfYs2cP9erVo169gi89atSoUaFckRs3bqRhw4a+taHBtJWTk0P//v2pVq0a8+fPJyYmpgLuWkREpOwckF9FNtSEggLKIEyaNIm0tDTi4uJ8m2f+/e9/s3btWhYsWABAjx49mDJlCs2bNycqKorHHnusUPqaUBs8eDCjR48mNTWVRo0aVWhbJRkxYgTPPfccQ4YM4bbbbmP27Nns3r27yE0wCxYsYOjQoYwbN47x48cXuHbrrbfyyCOPcO6553L55Zezdu1aZs6cyQMPPFCqtkaMGMH777/PK6+8wubNm33nExMTadu2begfgIiISKlZlZmuDgUFlEHo0aMHb7/9NuPHj/cl1G7Xrh2DBw/2lZk2bRrDhw9n6NChJCUlce+997Jx40Y2bNhQYf3q06cPo0ePDuvoJEDTpk1ZvHgxd999N+np6bRq1YqFCxfSqlWrUtXzwAMPEBMTw1//+leefPJJmjRpwrhx4/jd735XqraWLVvG8ePHufHGGwvU37VrV7Kyssp1ryIiIlKYVcCbDMMqOTnZrVu3LtzdqBQzZszgnnvu4dtvvy202UdERESKZmbrnXPFr2OrBI0urON+M/fKkNQ1qeOCsN+PRihPQ9u3b+ezzz5jwoQJDBo0SMGkiIjIaSgvgpLtRM6dnEHGjx/P1VdfTZs2bXjyySfD3R0RERE5w2mE8jSUkZFBRkZGuLshIiIiZeQw8p12eYuIiIhIOeRH0ERx5NxJOWRkZNCxY0dq1qxJkyZN6NevH+vXry9UbsqUKaXaJdytWzeGDBkSuo5WYatWraJz587UqFGDDh06sHTp0lLXMWvWLC655BLq1KlDYmIil112GStWrChQJj8/n9///vecffbZVK9enZSUFFauXFmgjJkV+bngggvKdY8iIiJStDM+oJwxYwbDhg2jX79+LF68mKlTp5KXl8eHH35YqGxpA8ozxc6dO+nTpw+tW7fmrbfeIiUlhb59+/Lll1+Wqp69e/eSlpZGRkYGb775Jq1atSItLa1AcP/HP/6RP/zhD9x///0sWLCApKQkevXqxbZt23xlVq9eXehTt27dIt90JCIiEg7OQZ6zkHyqgjN+ynvKlCnccccdjBs3zneuf//+FfraxEgzffp04uPjfTk6r7jiCrKyspg2bRpTp04Nuh7/fJPgyf+ZlZXF7NmzueiiiwB49tlnGT58OPfccw8Av/jFL2jSpAnTp0/nD3/4AwBdunQpUM+HH37Ijz/+yIABA8pzmyIiIiEVSWsoz/gRyh07dtCgQYNC52NjY31/PjllumPHDh599FHfsf90tnOOhx9+mHr16lGnTh0eeeSRUvclNzeX5ORk+vfvX+D8TTfdRIcOHapskLtixQpSU1N9zyw6Opq0tLRC09WlZWYkJCSQnZ0NwMGDB9m9ezeXXHKJr0xSUhKdOnUqceR49uzZnHPOOYUCTREREQmNMz6gbNu2LS+88AKLFi0iJyenyDInp00bNmzIrbfe6jv+/e9/7yszc+ZMJk+ezH333cerr77K8uXLWbNmTan6Uq1aNV5++WUWL17Ma6+9BsCbb77JggULeOWVVwoEuVXJli1baNmyJYAv+DvvvPPYsmVLmerLzc1l3759PPPMM2zatIlBgwYB+ALqU59DbGwsX331VZF15eXlkZmZyY033ohZ5PyXoIiInN48u7yjQvKpCs74Ke/nnnuOvn370rdvX2rVqkVaWhp33nknXbt29ZU5ObIVFxdHkyZNihzpmjJlCjfffDNjx44FoGPHjpxzzjml7k/btm2ZOHEio0ePpkOHDowYMYJx48bRsWPHMt5hxTtw4ACJiYnMmzePAQMGkJmZSWJiIseOHePEiROlCoQPHTrkS9RevXp1MjMz6d69OwB169YlISGBTZs2+crn5OTw6aef+gLZUy1dupTvvvtO090iIlLl5BE5Ax1VI6wNoy5duvDFF1/w4osv0rNnT95++226d+/OK6+8EnQdJ06c4PPPPy8QhDZp0oTzzz+/TH0aM2YMbdq0ISUlhebNm/uC1KouPj6ehIQE4uPjy1xHzZo1Wbt2LUuXLmXgwIEMGTKkwEjvr3/9a1544QXef/99fvjhB+6//34OHz5c7OjjnDlzOP/88+ncuXOZ+yQiIiIlO+MDSvCsw7v11lt544032LlzJxdddFGBTTqB7Nu3j/z8fOrUqVPg/KnHwYqKimLgwIEcOXKEAQMGEB0dXaZ6KktSUhLZ2dn07t2b7Oxs0tPTyc7Opnr16qWepo+KiiI5OZmrrrqKl156iZSUFB5++GHf9QkTJtC6dWsuu+wy6tevz/vvv8+wYcOoW7duobqOHj3KggULuOmmm8p9jyIiIqHk8GzKCcWnKlBAeYratWtzyy23sH37dvLy8oL6Tt26dYmKimL//v0Fzp96HKzdu3fz8MMP07VrV8aPH8/XX39dpnoqS8uWLdm6dWuBc9u2bfOtqyyPzp07F5jirlu3LitXrvS9z3zNmjV8++23tGvXrtB3Fy1axKFDhxRQiohIFRRZayirRi/C6Lvvvit0btu2bTRq1KjQyGBCQgKHDx8uVD4mJoY2bdoUSLC9a9cuvvjiizL16dZbb6VNmzYsX76ciy66iCFDhuCcK1NdlaF79+4sXbrUt6kpPz+fJUuW+NY++jtw4ACfffYZP/zwQ6FrRd3j+vXradiwYaHz55xzDhdccAHffvstf//73+nXr1+hMrNnz6Zdu3ZceOGFZbktERERCdIZvymnR48eXHrppfTq1Yv4+Hjee+89pk2bxkMPPVSobPv27Zk/fz69evWiUaNGJCUl0ahRIwDuvPNORo0aRZs2bejQoQOTJk0iLi6u1P2ZPn06WVlZbNiwgejoaF566SXat2/P1KlT+e1vf1vu+60II0aM4LnnnmPIkCHcdtttzJ49m927dzN69OhCZRcsWMDQoUMZN24c48ePL3AtJSWFG264gQ4dOgAwb9483n333QLvLf/Xv/7Fhg0baNeuHd9++y0TJkygefPmDB06tEBd+/btY8mSJaVauiAiIlKZ8iNoU84ZH1DeeeedzJ49mzfeeIMjR47QokULJk+ezF133VWo7MSJExk2bBh9+/bl4MGDDB482BfsDB8+nK+//pqnn36a3Nxc7rrrLqpVK93j3bp1K/fffz+PPfaY7zWBzZs356mnnmLMmDGkpaXRunXrct9zqDVt2pTFixdz9913k56eTqtWrVi4cCGtWrUqVT1du3Zl9uzZPP7445gZrVu35o033igw+hgbG8vMmTP54osvqFWrFn369OGpp54qFLy//vrr5OTkaHe3iIhIJbCqPJVaFsnJyW7dunXh7oaIiIhUUWa23jmXHM4+1GtTz1398tUhqSvj4oyw388ZP0IpIiIiEg5VZUNNKETOnYiIiIhIWGiEUkRERKSSeV69qE05IiIiIlIOkbTLW1PeIiIiIlIuCiiBjIwMOnbsSM2aNWnSpAn9+vVj/fr1hcpNmTKFrKysoOvt1q0bQ4YMCV1Hq7BVq1bRuXNnatSoQYcOHVi6dGmp68jMzCQ1NZUGDRpQu3ZtunXrxgcffFDqtoKtR0REJFz06sUIM2PGDIYNG0a/fv1YvHgxU6dOJS8vjw8//LBQ2dIGlGeKnTt30qdPH1q3bs1bb71FSkoKffv25csvvyxVPdOmTaNx48b85S9/ITMzk7POOosrr7ySTz75pFRtBVOPiIhIuEXSqxfP+DWUU6ZM4Y477ijwRpX+/ftz4sSJMPbq9DJ9+nTi4+PJyMggNjaWK664gqysLKZNm8bUqVODrmfBggXUrVvXd9ytWzeaNm3KCy+8wPTp04NuK5h6REREJHTCFtaaWVMze9fMNpnZRjO7u4gyZmZ/MrMtZvYfM+sc6n7s2LGDBg0aFDofGxvr3w/MjB07dvDoo4/6jv2ns51zPPzww9SrV486derwyCOPlLov2dnZ1KxZk5kzZxY4v23bNsyMf/7zn6WuszKsWLGC1NRU3zOLjo4mLS2NFStWlKoe/yAQPD+DFi1asH379lK1FUw9IiIiYRWi6W5NeUMucK9zri3QBRhlZm1PKZMOtPJ+hgMhH15q27YtL7zwAosWLSInJ6fIMqtXr2b16tU0bNiQW2+91Xf8+9//3ldm5syZTJ48mfvuu49XX32V5cuXs2bNmlL1JTExkeuuu45XX321wPk5c+bQsGFDfvnLX5b+BivBli1baNmyJeAJigHOO+88tmzZUq56jxw5wubNm7nwwgvL1VZR9YiIiISTw7PLOxSfqiBsAaVzbo9z7mPvnw8Cm4HGpxTrC7zsPD4EaptZo1D247nnniMvL4++fftSp04drr/+elauXFmgTJcuXejSpQtxcXE0adLEd3zeeef5ykyZMoWbb76ZsWPH0rt3b+bNm8fx48dL3Z/Bgwfz3nvv8fXXX/vOzZ07l0GDBhEdHV32G61ABw4cIDExkXnz5pGUlMTrr79OYmIix44dK9fSgcmTJ3PixAlGjx5drraKqkdERERCp0qs5DSz5kAn4KNTLjUGdvodf0PhoLNcunTpwhdffMGLL75Iz549efvtt+nevTuvvPJK0HWcOHGCzz//nK5du/rONWnShPPPP7/U/bnqqqs4++yzmT17NgAbNmxg06ZN3HLLLaWuq7LFx8eTkJBAfHx8uevKysriiSee4JlnnqF58+ZlbitQPSIiIuGiKe8QMrN44A3gt8657DLWMdzM1pnZuu+//77U309KSuLWW2/ljTfeYOfOnVx00UUFNukEsm/fPvLz86lTp06B86ceByMqKopf/epXvoByzpw5tG/fno4dO5a6rsqSlJREdnY2vXv3Jjs7m/T0dLKzs6levXqBtajB2rp1K9dffz1Dhw5l5MiRZW6rpHpERETCSWmDQsjMYvAEk7Odc28WUWQX0NTvuIn3XAHOuRnOuWTnXHL9+vXL1afatWtzyy23sH37dvLy8oL6Tt26dYmKimL//v0Fzp96HKzBgwezceNGPv74Y1577bUqPzrZsmVLtm7dWuDctm3bfGsdS+Onn36iT58+dOrUieeff77MbQWqR0REREInnLu8DfgrsNk598diii0CbvHu9u4CHHDO7QllP7777rtC57Zt20ajRo0KrVlMSEjg8OHDhcrHxMTQpk2bAmsvd+3axRdffFGmPrVp04aUlBRGjx7N7t27ufnmm8tUT2Xp3r07S5cu9W1qys/PZ8mSJXTv3r1Q2QMHDvDZZ5/xww8/FLqWk5ND//79qVatGvPnzycmJqZMbQVTj4iISLhF0ghlOPNQXgr8GvjEzDZ4zz0ENANwzv0ZeAvoBWwBjgBDQ92JHj16cOmll9KrVy/i4+N57733mDZtGg899FChsu3bt2f+/Pn06tWLRo0akZSURKNGnj1Cd955J6NGjaJNmzZ06NCBSZMmERcXV+Z+DR48mNGjR5Oamupro6oaMWIEzz33HEOGDOG2225j9uzZ7N69u8hNMAsWLGDo0KGMGzeO8ePHF6rn/fff55VXXmHz5s2+84mJibRt2zbotoKpR0REJJwcVScYDIWwBZTOuX9ByXvdnXMOGFWR/bjzzjuZPXs2b7zxBkeOHKFFixZMnjyZu+66q1DZiRMnMmzYMPr27cvBgwcZPHgwGRkZAAwfPpyvv/6ap59+mtzcXO666y6qVSv74+3Tpw+jR4+u8qOTAE2bNmXx4sXcfffdpKen06pVKxYuXEirVq1KVc+yZcs4fvw4N954Y4HzXbt29b2hKJi2gqlHREREQsc8MVvkSE5OduvWrQt3N8ptxowZ3HPPPXz77bckJCSEuzsiIiIRw8zWO+eSw9mHpNYN3C9mDAhJXUu6/ins93PGv3qxqtm+fTufffYZEyZMYNCgQQomRUREIpEjoqa8w542SAoaP348V199NW3atOHJJ58Md3dEREREAtIIZRWTkZHhW5cpIiIikelkHspIoYBSREREJAwiKaDUlLeIiIiIlItGKEVEREQqmfJQioiIiEi5uQgKKDXlLSIiIiLlohFKERERkTDIL/mFgacVBZQiIiIilcwpsbmIiIiIyH9phFJEREQkDLQpR0RERETESyOUIiIiIpVOeShFREREpJw05S0iIiIi4qURShEREZFK5oistEEKKEVEREQqm/PkoowUmvIWERERkXLRCKWIiFRpqVE3BFVuaf7rFdwTkdDSqxdFREREpMwc2uUtIiIiIuKjgFJERESk0nkSm4fiE7Als5fM7Dsz+9Tv3Hgz22VmG7yfXn7XHjSzLWb2uZn1DOZuNOUtIiIiEgaVuMs7A5gGvHzK+Wedc0/7nzCztsAA4ELgbGCZmZ3vnMsrqQGNUIqIiIhEMOfcKmBfkMX7Aq855447574CtgAXB/qSAkoRERGRMHDOQvIph9Fm9h/vlHgd77nGwE6/Mt94z5VIAaWIiIhIJXMupAFlPTNb5/cZHkQXpgPnAf8D7AGeKc/9aA2liIiIyOntB+dccmm+4Jzbe/LPZjYT+If3cBfQ1K9oE++5EmmEUkRERCQMKmuXd1HMrJHf4XXAyR3gi4ABZhZnZi2AVsCaQPVphFJERKo0vQFHIlVl7fI2s7lANzxT498A44BuZvY/eHKsbwdu9/TJbTSzTGATkAuMCrTDGxRQioiIiEQ059zAIk7/tYTyjwOPl6YNBZQiIiIiYRBJr15UQCkiIiJSyRzlTvlTpYR1U05RrwI65Xo3Mzvg91qgRyq7jyIiIiJSsnCPUGZQ9KuA/L3nnOtTOd0RERERqRyV9+bFihfWgNI5t8rMmoezDyIiIiKVzkXWGsrTIQ/lz83s/8zsn2Z2Ybg7IyIiIiIFhXvKO5CPgXOcc4fMrBewEE+CzQK8rxgaDtCsWbPK7aGISBWTGn1TwDIWFdzIyDs5r5WrnWAE05eS+gGQGnVDUG0FymkZzD0tzZsXVFsiAUXQnHeVHqF0zmU75w55//wWEGNm9YooN8M5l+ycS65fv36l91NERESktEL4Lu+wq9IBpZk1NDPz/vliPP39Mby9EhERERF/YZ3yLuZVQDEAzrk/A9cDI8wsFzgKDHCusl5UJCIiIlJxIimiCfcu76JeBeR/fRqetEIiIiIiEcOhXd4iIiIiIj5VfZe3iIiISORxgEYoRUREREQ8NEIpIme8UOUwrAy/jB1Uqe31rH5zpbZXnFDlvAxFPcH+fQmkvH+fTqe/t1I0bcoREZEqKyo2NnAhl1857UQHngizqJLL5B0+EriO6OiAZVx+EL+9Q/BcRIIWQQGlprxFREREpFw0QikiIiJS6arOW25CQQGliIiISDhoyltERERExEMjlCIiIiKVzUXWm3IUUIqIiIiEg6a8RUREREQ8NEIpIqe1YJI7hyqxc0lJsaNiAv/fqdWoEbCMO3a8VH2S01egv7tKSH4m0JS3iEjksNNnsib6Z/UCF6oWOMk3h48GLhMVgl92MTEBi7gAicuja9Usfz8Al5cXsEx+oIA+yMTngYLBULxtRwFnBNCUt4iIiIiIh0YoRURERMIhgkYoFVCKiIiIVDYHRFDaIE15i4iIiEi5aIRSREREJAycprxFREREpFwiKKDUlLeIiIiIlItGKEUk5EKRY++kUOTa+2XsoBD0REQkxCJoU44CShE5vQVISm7RgZN8R9dJCtxMCJKfu4Z1AxeqHlvi5eO14wLXEVzu7YBya5X87GKycwLWEfP94YBlQvIr1QLX4vbtD1gmqnocbx9+udjrJb0tSaS0TFPeIiIiIiIeGqEUERERqWyOiNqUo4BSREREpNJZRK2h1JS3iIiIiJSLRihFREREwkFT3iIiIiJSLgooReRMFso8kyJVTYl/v4NMH6V/R+RMo4BSRCpEoITkv4wZEJJ2omvVLPG6JSUGriQmJmCRnCZ1Srx+onbJ+SMBjtUJnBPzUJOSF+nnlny7HkGs8692JIh6AoyeVP8x8K+QesdyAzdTu0aJ16vtPRCwDo6fCFwmLy9gkfwTAepxIUrySWiS9stpLoJGKLUpR0RERETKRSOUIiIiIpXNEVFpgxRQioiIiISBXr0oIiIiIuIV1oDSzF4ys+/M7NNirpuZ/cnMtpjZf8ysc2X3UURERKRCuBB9qoBwj1BmAGklXE8HWnk/w4HpldAnERERESmFgAGlmd1pZiXnyygj59wqYF8JRfoCLzuPD4HaZtaoIvoiIiIiImUTzKacBsBaM/sYeAl42zlXWQOsjYGdfsffeM/tqaT2RaSMQpVnMr3BiJDUIyJS1UTSppyAAaVz7mEz+z3wS2AoMM3MMoG/Oue2VnQHg2Fmw/FMidOsWbMw90bkDBDgbSEWFTgVhsUGTgQeFV8rcD1NSp60yIuvHrCO4/UCl/m+c8nJz3NLzs3t6UvLwNnErzzvyxKvN4oLnOS7TY3dAcvM3XNxwDKffNW4xOuxa+IC1pGTFLhM9a9KmqgCooJYnRXMOEcQb7mx6JKTz7sgkqOLBC2C0gYFtYbSOyL5rfeTC9QB5pvZUxXYN4BdQFO/4ybec6f2b4ZzLtk5l1y/fv0K7pKIiIiI+AtmDeXdZrYeeAp4H2jvnBsBXAT0r+D+LQJu8e727gIccM5pultEREROb6Ha4V1Fps2DWUN5FtDPObfD/6RzLt/M+pSncTObC3QD6pnZN8A4IMZb/5+Bt4BewBbgCJ4pdxEREZHTXxUJBkMhmDWU40q4trk8jTvnBga47oBR5WlDRERERCqWXr0oIiIiEgZn1C5vEREREakAERRQhvtNOSIiIiJymtMIpYgUkBp1Q+BCQeTzExGRACJohFIBpYiUWlRMyf/XYTUCZ/m26oETXhNEYvMTP4sv8frBcwInLf+xY+Cu5J11osTrt1+8MmAdD5y1JWCZmQdKTtS+NycpYB1vfHdRwDJfH6gdsEzUjyUnn48+HrAKLDeI35iBkpIfPBSwirx9PwXuS4C/t+BJuP/24ZcDlhMpL3ORtYZSwwwiIiIiUi4aoRQREREJhwh69aICShEREZFw0JS3iIiIiIiHRihFREREwiCSNuUooBQREREJhwgKKDXlLSIiIiLlohFKkSogqGTiwNL810uuJ/qmUHRH5IwW7L+PJQn076oIEZaHUgGliBQUxFtwAiUud8cCZ7y2hvUDljlybp2AZbKbx5R4PadmwCqof+F3Acvc1Gx9idcvqL47YB0/5h8O3E7CthKvD97WN2Adn+w8O2CZuE2Bk88nZJd8vcaPeQHriP3xSMAyHC85aXx+EInNcfmBi5w4wTs5r5VYJhTBpEjQIiig1JS3iIiIiJSLRihFREREwiGCRigVUIqIiIiEQSStodSUt4iIiIiUiwJKERERESkXBZQiIiIiUi5aQylyGgmY0iSIlD8iZzKlBZIqpZLWUJrZS0Af4DvnXDvvubOAeUBzYDtwo3Nuv5kZMBXoBRwBhjjnPg7UhgJKkSogmCTIlfaLMJh8fjk5JV6Pqhs4f2R+fPWAZY7Xjg5cT8lpKDncLPD9HNpTO2CZV/NTSrzeqf6ugHU8dyQpYJmfjpWcH/K7T38WsI5qRyxgmaSvAv8mS9p8oMTrUQeCyDF54GDAInkHSk546fIC57sM5u9tsJSUXCpF5SY2zwCmAS/7nRsLLHfOTTazsd7jB4B0oJX3cwkw3fvPEmk4Q0RERCSCOedWAftOOd0XmOX98yzgWr/zLzuPD4HaZtYoUBsKKEVERETCwYXoUzYNnHN7vH/+Fmjg/XNjYKdfuW+850qkKW8RERGRcAjdlHc9M1vndzzDOTcj6G4458zKNwGvgFJERETk9PaDcy65lN/Za2aNnHN7vFPa33nP7wKa+pVr4j1XIk15i4iIiFQyw7MpJxSfMloEDPb+eTDwN7/zt5hHF+CA39R4sTRCKSIiIhIOlZc2aC7QDc/U+DfAOGAykNPgsVUAABGWSURBVGlmtwI7gBu9xd/CkzJoC560QUODaUMBpYiIiEgEc84NLOZSjyLKOmBUadtQQCkiIiJS2So3D2WFU0ApEklcfolJmX8ZOyhgFRYbOOF4VEJ8gAKBl2cHkxS75rdxAcu4qNgSr+fUCtyX/H0l1wFw4Pu6JV7/15F6AeuI2x+wCDGHSv4N02hf4ATetb4KnEw8an/JycQB3KHDJV7PP3Y8cB3HgygTROLyQMnGg038r6TlUqVEUECpTTkiIiIiUi4aoRQREREJhwgaoVRAKSIiIhIGkbSGMqxT3maWZmafm9kW74vJT70+xMy+N7MN3s9t4einiIiIiBQvbCOUZhYNPA+k4nlP5FozW+Sc23RK0XnOudGV3kERERGRiqQRypC4GNjinNvmnDsBvAb0DWN/RERERCqHC+GnCghnQNkY2Ol3/I333Kn6m9l/zGy+mTUt4rqIiIiIhFFV35Tzd2Cuc+64md0OzAK6n1rIzIYDwwGaNWtWuT0UqWJKysdn1WIqsSciIlKSSNqUE86AchfgP+LYxHvOxzn3o9/hi8BTRVXknJsBzABITk6OoB+PyH8Fk5A5YHJnFzgpdlBcgH/NoixwFdUDB7dxuw8ELBN9tFaJ1xO+DjwRY7mBn0te9ZL/7zI/JvA9x+4LnOQ7ev+hkgscPRawDnJyApexIH5GAdoKJrF5ZVHCcjktRVDEEs4p77VAKzNrYWaxwABgkX8BM2vkd3gNsLkS+yciIiIiQQjbCKVzLtfMRgNvA9HAS865jWY2AVjnnFsE3GVm1wC5wD5gSLj6KyIiIhJKmvL+/+3dfZBW5XnH8d9PMIgvzQtYQgNGO1FaphnXzJaKZjJGiIrVkFaDZkYLnXTW0dgJMzodMZq+6Eg6oy39IzruqIXRKFAMI9pVBNTJH5jUtSX1FUMZO0AxFExaG1942at/PIdkl93lWfec59zP3vl+Zs7sOc99POfi+Axce93nvu+KRESPpJ4jPvt2v/0lkpbUHRcAAEDLZZRQspY3AAAASmn3Ud4AAAD5aaM5JKtAQgkAAFAzF1suSCiBNtB0uh8AANoYCSXwayT6mvevuNkck5IO7Xv7qO3HTJzY/D4HDjY9ZyTGv/OLJieMa36REcQy7oPycy42m9dRkvoOHSp9H7n56/FxcARzVUp6+sDKYdv4RQgoiS5vAAAAlJHTtEGM8gYAAEApVCgBAABSoEIJAAAANFChBAAASCGjCiUJJQAAQN2CQTkAAADAL1GhBAAASCGjCiUJJdAGNvT9UyXXaTrRdPQ1vUbf+80n8Pa4o08WHvv3N71GNJkcXZI8YULz6xwsP0F6jGBi82g24fgInm2z5zYSTeMoVPWdSn0PIGd0eQMAAAAFKpQAAAApZFShJKEEAABIgC5vAAAAoECFEgAAoG4hurwBAABQUkYJJV3eAAAAKIUKJdAGms4fWWDePwDIg5XXoBwSSuDXTLOkdCTJbRzsO+p1Ljj2yubXGMkE3SOYIF2ur6Nlw6FVw7aN6LkdOsQvBQB+JaOEki5vAAAAlEKFEgAAIAFHPiVKEkoAAIC6ZTZtEF3eAAAAKIUKJQAAQAKM8gYAAEA5GSWUdHkDAACgFCqUAAAACdDlDaBSVU12XcV1qrjG0wdWlr7GWMOE5QA+tIwSSrq8AQAAUAoVSgAAgLoFXd4AAAAoK6OEMmmXt+2LbG+1vc32TUO0T7C9qmj/ke1T648SAAAAR5MsobQ9TtJ3Jc2TNFPS12zPPOK0r0v6WUR8RtLfS/rbeqMEAAContXo8q5iawcpK5SzJG2LiO0RsV/SSknzjzhnvqQVxf4aSXNsu8YYAQAAWiOimq0NpHyH8lOSdvQ73inpD4Y7JyIO2v4fSZMk7R3uolu3btV555034LMFCxbouuuu07vvvquLL7540H+zaNEiLVq0SHv37tXll18+qP3aa6/VFVdcoR07dujqq68e1H7DDTfo0ksv1datW3XNNdcMar/llls0d+5cbdmyRYsXLx7Ufscdd+icc87R5s2bdfPNNw9qX7ZsmTo6OrRx40bdfvvtg9rvvfdezZgxQ48//rjuuuuuQe0PPvigpk+frlWrVumee+4Z1L5mzRpNnjxZy5cv1/Llywe19/T06Pjjj9fdd9+t1atXD2p/7rnnJEl33nmnnnjiiQFtEydO1JNPPilJuu2227Rp06YB7ZMmTdKjjz4qSVqyZImef/75Ae3Tpk3TQw89JElavHixtmzZMqD9jDPOUHd3tySpq6tLb7zxxoD2jo4OLVu2TJJ01VVXaefOnQPaZ8+eraVLl0qSLrvsMu3bt29A+5w5c3TrrbdKkubNm6f33ntvQPsll1yiG2+8UZIGfe8kvnt89/ju8d3ju3ekdvjuoXpZDMqx3SWpS5ImTJiQOBoAAIDm2qW7ugqORKVS27Ml/VVEXFgcL5GkiFja75z1xTnP2x4v6S1JJ8dRgu7s7Ize3t7WBg8AAMYs2y9GRGfKGE6cND0+e+HgCupo/PCRG5P/eVK+Q/mCpNNtn2b7I5KulLTuiHPWSVpY7F8u6ZmjJZMAAABjQlS4tYFkXd7FO5HXS1ovaZykByLiFdt/I6k3ItZJul/Sg7a3SXpbjaQTAABgzHNf6giqk/QdyojokdRzxGff7rf/vqSv1h0XAAAARi6LQTkAAABjTpt0V1eBhBIAACCBnEZ5J116EQAAAGMfFUoAAIC6hdpmlZsqkFACAAAkQJc3AAAAUKBCCQAAkEJGFUoSSgAAgJpZdHkDAAAAv0SFEgAAoG4RjPIGAABAOXR5AwAAAAUqlAAAAClkVKEkoQQAAEiALm8AAACgQIUSAACgbiGpL58SJQklAABACvnkkySUAAAAubP9pqR3JB2SdDAiOm1/QtIqSadKelPSgoj42WiuzzuUAAAACTiq2T6EL0ZER0R0Fsc3SdoUEadL2lQcjwoJJQAAQAqHV8spu43efEkriv0Vkr4y2guRUAIAAOQvJD1t+0XbXcVnUyJid7H/lqQpo70471ACAAAkUOE8lJNt9/Y77o6I7iPO+XxE7LL9m5I22H69f2NEhD36iEgoAQAA6haqcpT33n7vRQ59u4hdxc89ttdKmiXpp7anRsRu21Ml7RltAHR5AwAAZMz2CbZPOrwv6QJJL0taJ2lhcdpCSY+N9h5UKAEAAGpmSS43oObDmCJprW2pkfs9HBFP2X5B0mrbX5f0n5IWjPYGJJQAAAAp9NVzm4jYLunMIT7fJ2lOFfegyxsAAAClUKEEAABIoMYu75ajQgkAAIBSqFACAADUrdppg5IjoQQAAKhd6WUT2wpd3gAAACiFCiUAAEACFS69mBwJJQAAQAp0eZdj+xO2N9j+SfHz48Ocd8j2lmJbV3ecAAAAaC7VO5Q3SdoUEadL2lQcD+W9iOgoti/XFx4AAEALheS+arZ2kCqhnC9pRbG/QtJXEsUBAACQRkQ1WxtIlVBOiYjdxf5baixaPpTjbPfa/qFtkk4AAIA21LJBObY3SvrkEE3f6n8QEWEPO87p0xGxy/ZvS3rG9ksR8R9D3KtLUpcknXLKKSUjBwAAqEF7FBcr0bKEMiLmDtdm+6e2p0bEbttTJe0Z5hq7ip/bbT8n6SxJgxLKiOiW1C1JnZ2dGf3vAQAAuWIt7/LWSVpY7C+U9NiRJ9j+uO0Jxf5kSedKerW2CAEAADAiqRLK70j6ku2fSJpbHMt2p+37inN+V1Kv7R9LelbSdyKChBIAAOQho0E5SSY2j4h9kuYM8XmvpD8r9jdL+mzNoQEAALReSGqTKX+qwFreAAAAKIWlFwEAAGpmRVaDckgoAQAAUsgooaTLGwAAAKVQoQQAAEghowolCSUAAEDdGOUNAAAA/AoVSgAAgAQY5Q0AAIByMkoo6fIGAABAKVQoAQAAatc+63BXgYQSAACgbqGsEkq6vAEAAFAKFUoAAIAUMpqHkoQSAAAggZymDaLLGwAAAKVQoQQAAEiBCiUAAADQQIUSAACgbiGpL58KJQklAABA7fKa2JwubwAAAJRChRIAACCFjCqUJJQAAAApZJRQ0uUNAACAUqhQAgAA1I1R3gAAACgnpMhnMW+6vAEAAFAKFUoAAIAUMhqUQ0IJAABQt8zeoaTLGwAAAKVQoQQAAEiBLm8AAACUklFCSZc3AAAASqFCCQAAULvIqkJJQgkAAFC3kNTHxOal2P6q7Vds99nuPMp5F9neanub7ZvqjBEAAAAjk+odypcl/bGkHwx3gu1xkr4raZ6kmZK+ZntmPeEBAAC0WEQ1WxtI0uUdEa9Jku2jnTZL0raI2F6cu1LSfEmvtjxAAACAVmuTZLAK7TzK+1OSdvQ73ll8BgAAgDbSsgql7Y2SPjlE07ci4rGK79Ulqas4/MD2y1VeHwNMlrQ3dRCZ4tm2Ds+2dXi2rcOzbZ0ZqQOQIqulF1uWUEbE3JKX2CVper/jacVnQ92rW1K3JNnujYhhB/qgHJ5v6/BsW4dn2zo829bh2baO7d7UMTRmDWKUdx1ekHS67dNsf0TSlZLWJY4JAAAAR0g1bdAf2d4pabakf7a9vvj8t2z3SFJEHJR0vaT1kl6TtDoiXkkRLwAAQOX6opqtDaQa5b1W0tohPv8vSRf3O+6R1PMhL99dLjo0wfNtHZ5t6/BsW4dn2zo829Zpj2eb0ShvR0Z/GAAAgLHgo+NPjtknza/kWut/fv+Lqd+3ZelFAACAukWw9GK7YqnG1rD9gO09TMdUPdvTbT9r+9ViOdJvpo4pF7aPs/0vtn9cPNu/Th1TbmyPs/1vtp9IHUtubL9p+yXbW9piRHJGbH/M9hrbr9t+zfbsZMFktFJONgklSzW21HJJF6UOIlMHJd0QETMlnS3pG3xvK/OBpPMj4kxJHZIusn124phy8001Bk2iNb4YER2puzIz9A+SnoqI35F0pvgOVyKbhFL9lmqMiP2SDi/ViJIi4geS3k4dR44iYndE/Gux/44af7GxIlQFouH/isNji609fpXPgO1pkv5Q0n2pYwFGyvZHJX1B0v2SFBH7I+LnqeKJvr5KtnaQU0LJUo0Y02yfKuksST9KG0k+ii7ZLZL2SNoQETzb6iyT9BeS2uNfs/yEpKdtv1isBodqnCbpvyX9Y/G6xn22T0gdVA5ySiiBMcv2iZIelbQ4Iv43dTy5iIhDEdGhxkpbs2z/XuqYcmD7Ekl7IuLF1LFk7PMR8Tk1XuP6hu0vpA4oE+MlfU7SPRFxlqRfSEo05qKi9yd5h7JyI16qEWgnto9VI5n8XkR8P3U8OSq6tJ4V7wJX5VxJX7b9phqvF51v+6G0IeUlInYVP/eoMW/zrLQRZWOnpJ39eivWqJFg1i+U1cTmOSWULNWIMce21XiX57WI+LvU8eTE9sm2P1bsT5T0JUmvp40qDxGxJCKmRcSpavxd+0xEXJU4rGzYPsH2SYf3JV0giVk2KhARb0naYXtG8dEcSa8mDCkb2cxDGREHbR9eqnGcpAdYqrEath+RdJ6kycWSmX8ZEfenjSob50q6WtJLxbt+knRzsUoUypkqaUUxA8QxaizfyvQ2GAumSFrb+H1T4yU9HBFPpQ0pK38u6XtF8Wm7pD9NFknk8woyK+UAAADU7DeOmRRnj7+wkmttOPBI8pVycuryBgAAQALZdHkDAACMGRFZdXmTUAIAACQQbTJCuwp0eQMAAKAUKpQAAAApZNTlzShvAACAmtl+StLkii63NyKSLtxAQgkAAIBSeIcSQJZs/77tf7d9XLHyyCus5Q0ArUGFEkC2bN8u6ThJE9VYv3dp4pAAIEsklACyVSyt9oKk9yWdExGHEocEAFmiyxtAziZJOlHSSWpUKgEALUCFEkC2bK+TtFLSaZKmRsT1iUMCgCwxDyWALNn+E0kHIuJh2+MkbbZ9fkQ8kzo2AMgNFUoAAACUwjuUAAAAKIWEEgAAAKWQUAIAAKAUEkoAAACUQkIJAACAUkgoAQAAUAoJJQAAAEohoQQAAEAp/w8uOBID3+3KugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10, 6))\n",
    "counts, xedges, yedges, im = ax.hist2d(x1_all, x2_all, bins=[60, 40], range= [[0.0, 6.0], [-1.0, 3.0]], cmin=1)\n",
    "ax.plot([0.0, 6.0], [0.0, 0.0], \"--k\")   # NOTE: This draws a line from [x1, x2], [y1, y2] with dashed line (\"--\") and in black (\"k\")\n",
    "fig.colorbar(im) # ticks=[-1, 0, 1]\n",
    "\n",
    "ax.set(title='Histogram of lengths (x) and widths (y)',\n",
    "       xlabel='x', \n",
    "       ylabel='y',\n",
    "       aspect='equal', # NOTE: This forces the x- and y-axis to have the SAME scale!!!\n",
    "      )\n",
    "\n",
    "d = {'Entries': len(x12_all),\n",
    "     'Mean x': x1_all.mean(),\n",
    "     'Mean y': x2_all.mean(),\n",
    "     'Std x': x1_all.std(ddof=1),\n",
    "     'Std y': x2_all.std(ddof=1),\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.97, text, ax, fontsize=15);\n",
    "\n",
    "fig.tight_layout()\n",
    "fig\n",
    "\n",
    "if save_plots :\n",
    "    fig.savefig(\"Dist_2Dgauss.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we bin `y_all` and fit it with a Gaussian distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian(x, N, mu, sigma):\n",
    "    return N * 1. / (sigma*np.sqrt(2*np.pi)) * np.exp(-0.5* (x-mu)**2/sigma**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlkAAAFpCAYAAACvaj13AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFT9JREFUeJzt3W+o5Nd5H/DvE9uli+2yNlKFspK6paiBpFA5LHJhQ1FlmsjaUCUQhA11nNRl80ICmwaadd4kbTHcF41SF1rTTSxiU9uKwDYWkUijah2MRfxHUhXbkuJWJGuky1qr1NlYJkuKlKcv7qx9d3fu3r9nZ+bO5wOXmTnzm5mz/Hb3fuec53dOdXcAANhbPzTrDgAA7EdCFgDAAEIWAMAAQhYAwABCFgDAAEIWAMAAQhYAwABCFgDAAEIWAMAAQhYAwACvn3UHkuSaa67pw4cPz7obAACbevLJJ/+8u6/d7Li5CFmHDx/OE088MetuAABsqqq+tZXjTBcCAAwgZAEADCBkAQAMIGQBAAwgZAEADCBkAQAMIGQBAAywaciqqhur6vNV9WxVPVNV75+0/3pVrVbV05OfO9e95oNV9XxVfbOqfmrkHwAAYB5tZTHSV5P8cnc/VVVvTvJkVT06ee43u/s/rj+4qn40ybuS/FiSH07yP6vqH3b3a3vZcQCAebbpSFZ3n+nupyb3X0nyXJJDV3jJXUke6O6/7u4/S/J8klv3orMAAItiWzVZVXU4yduSfHnSdG9Vfa2q7q+qt0zaDiV5Yd3LXsyVQxkAwL6z5ZBVVW9K8ukkH+ju7yb5SJJ/kOSWJGeS/MZ2PriqjlfVE1X1xMsvv7ydlwIAzL0thayqekPWAtYnuvszSdLdL3X3a939N0l+Kz+YElxNcuO6l98wabtId5/s7iPdfeTaazfdyBoAYKFsWvheVZXko0me6+771rVf391nJg9/Nsk3JvcfSvLJqrova4XvNyf5yp72GmALjq6cyuq58xe1HTp4II+fuH1GPQKWyVauLjya5D1Jvl5VT0/afjXJu6vqliSd5HSSX0qS7n6mqh5M8mzWrky8x5WFwCysnjuf0yvHLmo7fOLhGfUGWDabhqzu/mKSmvLUI1d4zYeSfGgX/QIAWGhWfAcAGEDIAgAYQMgCABhAyAIAGEDIAgAYQMgCABhAyAIAGEDIAgAYQMgCABhAyAIAGEDIAgAYQMgCABhAyAIAGOD1s+4AwKI4unIqq+fOX9Z+6OCBPH7i9hn0CJhnQhbAFq2eO5/TK8cuaz984uEZ9AaYd6YLAQAGELIAAAYwXQgsnGm1UeqigHkjZAELZ1pt1G7qooQ2YAQhC1h6ex3aABI1WQAAQwhZAAADmC4E9oVDBw9cNsV36OCBGfVmOrVfsFyELGBfWISgovYLlovpQgCAAYQsAIABhCwAgAHUZAHs0kZF94tQJwaMI2QB7NK0MKWgHTBdCAAwgJAFADCAkAUAMICQBQAwgJAFADCAkAUAMIAlHICZsFkysN8JWcBM7PfNkjdaoBRYHkIWwABG5AA1WQAAAwhZAAADmC4EmEJNFbBbQhbAFGqqgN0SsgBmaKMRMyEPFp+QBTBD08LUflrKApaZwncAgAGELACAAYQsAIAB1GQBS2URlmZQDA/7g5AFLJVFCCqK4WF/MF0IADCAkAUAMICQBQAwgJosgAV1dOVUVs+dv6hNgTzMj01DVlXdmOTjSa5L0klOdveHq+qtSX43yeEkp5Pc3d1/UVWV5MNJ7kzyV0l+obufGtN9gOW1eu58Tq8cu6hNgTzMj62MZL2a5Je7+6mqenOSJ6vq0SS/kOSx7l6pqhNJTiT5lSTvTHLz5OftST4yuQXYto1GawDm3aYhq7vPJDkzuf9KVT2X5FCSu5LcNjnsY0n+MGsh664kH+/uTvKlqjpYVddP3gdgW6aN1gAsgm0VvlfV4SRvS/LlJNetC07fztp0YrIWwF5Y97IXJ20AAEtjy4XvVfWmJJ9O8oHu/u5a6dWa7u6q6u18cFUdT3I8SW666abtvBRg6SzCSvXAxbYUsqrqDVkLWJ/o7s9Mml+6MA1YVdcnOTtpX01y47qX3zBpu0h3n0xyMkmOHDmyrYAGsGxcMQiLZ9PpwsnVgh9N8lx337fuqYeSvHdy/71JPreu/edrzT9J8pfqsQCAZbOVkayjSd6T5OtV9fSk7VeTrCR5sKrel+RbSe6ePPdI1pZveD5rSzj84p72GABgAWzl6sIvJqkNnn7HlOM7yT277BcAO7BR7ZbpRrj6rPgOsI9MC1MWKIXZsHchAMAAQhYAwABCFgDAAEIWAMAAQhYAwABCFgDAAEIWAMAAQhYAwABCFgDAAEIWAMAAQhYAwAD2LgTY52waDbMhZAHsczaNhtkQsoAdO7pyKqvnzl/UZoQEYI2QBezY6rnzOb1y7KK2oyunTE0BRMgC9pipKYA1ri4EABhAyAIAGMB0ITDcRksIAOxnQhYw3FaL3oWx+eMKUtg5IQuYG35xz59pV5C6kAG2Rk0WAMAARrIA2BOmFuFiQhYAe8LUIlzMdCEAwABCFgDAAKYLAUiycU0VsDNCFgBJptdUATtnuhAAYAAhCwBgACELAGAAIQsAYAAhCwBgAFcXAlvi8n6A7RGygC1xef/+cujggcu2vBGaYW8JWQBLyKbNMJ6aLACAAYQsAIABhCwAgAGELACAAYQsAIABhCwAgAGELACAAayTBUtso1XcraHElUxbyPRCO/ADQhYsiY0C1aWruE/75QnrCeGwNUIWLAnb4gBcXWqyAAAGELIAAAYQsgAABhCyAAAGELIAAAYQsgAABhCyAAAG2DRkVdX9VXW2qr6xru3Xq2q1qp6e/Ny57rkPVtXzVfXNqvqpUR0HAJhnWxnJ+p0kd0xp/83uvmXy80iSVNWPJnlXkh+bvOa/VtXr9qqzAACLYtOQ1d1fSPKdLb7fXUke6O6/7u4/S/J8klt30T8AgIW0m5qse6vqa5PpxLdM2g4leWHdMS9O2gAAlspO9y78SJL/kKQnt7+R5F9t5w2q6niS40ly00037bAbwF47dPDA1E2iDx08MIPeACyuHYWs7n7pwv2q+q0kvzd5uJrkxnWH3jBpm/YeJ5OcTJIjR470TvoB7L3HT9w+6y6whI6unMrqufMXtR06eMDfRxbajkJWVV3f3WcmD382yYUrDx9K8smqui/JDye5OclXdt1LAPa11XPnc3rl2EVt00ZUYZFsGrKq6lNJbktyTVW9mOTXktxWVbdkbbrwdJJfSpLufqaqHkzybJJXk9zT3a+N6ToAwPzaNGR197unNH/0Csd/KMmHdtMpAIBFt9PCd2CObVTfAsDVI2TBPjStvgXmhS8BLAshC4CraqtfAqYtJ+KKQxaJkAXAXJoWplxxyCLZzYrvAABsQMgCABhAyAIAGEDIAgAYQMgCABhAyAIAGMASDrBANlrE0bpBAPNHyIIFMm0RR+sGAcwn04UAAAMIWQAAAwhZAAADCFkAAAMofAdgmEMHD1x2ccahgwdm1Bu4uoQsAIaxvAjLzHQhAMAAQhYAwABCFgDAAEIWAMAAQhYAwABCFgDAAEIWAMAAQhYAwABCFgDAAFZ8B2BhbLRNj5XlmUdCFgALY1qYujR0wbwQsmDB2YAXYD4JWbDgTJMAzCeF7wAAAxjJAmApHF05ldVz5y9qUzTPSEIWAEth9dz5nF45dlGbonlGMl0IADCAkAUAMIDpQgAWmgVKmVdCFgALzQKlzCshC+bAtKueprHIKMDiELJgDky76gmAxabwHQBgACELAGAAIQsAYAAhCwBgACELAGAAVxfCVbbRJrUA7C9CFlxllmsAWA6mCwEABhCyAAAGELIAAAYQsgAABhCyAAAGELIAAAYQsgAABtg0ZFXV/VV1tqq+sa7trVX1aFX9n8ntWybtVVX/uaqer6qvVdWPj+w8AMC82spI1u8kueOSthNJHuvum5M8NnmcJO9McvPk53iSj+xNNwEAFsumIau7v5DkO5c035XkY5P7H0vyM+vaP95rvpTkYFVdv1edBQBYFDutybquu89M7n87yXWT+4eSvLDuuBcnbQAAS2XXhe/d3Ul6u6+rquNV9URVPfHyyy/vthsAAHNlpyHrpQvTgJPbs5P21SQ3rjvuhknbZbr7ZHcf6e4j11577Q67AQAwn16/w9c9lOS9SVYmt59b135vVT2Q5O1J/nLdtCIAXBWHDh7I4RMPX9a2U0dXTmX13PnL3u/xE7fv+D3Z/zYNWVX1qSS3Jbmmql5M8mtZC1cPVtX7knwryd2Twx9JcmeS55P8VZJfHNBnALiivQ4/q+fO5/TKsYvaLg1xcKlNQ1Z3v3uDp94x5dhOcs9uOwUAsOis+A4AMICQBQAwwE4L34FLKIwFYD0hC/aIwljYHzb6wgTbJWQBwDrTvjDBTqjJAgAYQMgCABhAyAIAGEBNFgy011t7ALA4hCwYyPINAMvLdCEAwABCFgDAAEIWAMAAarIAWFouTmEkIQuApeXiFEYyXQgAMICQBQAwgOlC2IGjK6eyeu78RW3qOABYT8iCHVg9dz6nV47NuhsAzDHThQAAAwhZAAADCFkAAAOoyQKAObHRRTXW81pMQhYAzIlpF9VcuiI9i8N0IQDAAEIWAMAApgsBYLDd1FpttIm1Oq35J2QBwGC7qbWaFqbUaS0G04UAAAMIWQAAA5guBIA9Mq32Kpm+gfxGtVbsH0IWAOyR7Wwer3B9/zNdCAAwgJAFADCAkAUAMICQBQAwgMJ32MRGKzUDwJUIWbCJ7VwtBAAXmC4EABhAyAIAGEDIAgAYQMgCABhAyAIAGEDIAgAYQMgCABhAyAIAGMBipACwA4cOHsjhEw9f1gYXCFkAsAOPn7h91l1gzpkuBAAYQMgCABhAyAIAGEDIAgAYQMgCABhgV1cXVtXpJK8keS3Jq919pKremuR3kxxOcjrJ3d39F7vrJgDAYtmLkax/1t23dPeRyeMTSR7r7puTPDZ5DACwVEask3VXktsm9z+W5A+T/MqAz4EdO7pyKqvnzl/WfujgAWvfALAndhuyOskfVFUn+W/dfTLJdd19ZvL8t5Nct8vPgD23eu58Tq8cu6z90tWbAWCndhuyfqK7V6vq7yZ5tKr+ZP2T3d2TAHaZqjqe5HiS3HTTTbvsBgAst2kj9EbnZ2tXIau7Vye3Z6vqs0luTfJSVV3f3Weq6vokZzd47ckkJ5PkyJEjU4MY7IWN/uMBWFQb7Zt46Qi90fnZ2nHIqqo3Jvmh7n5lcv8nk/z7JA8leW+Slcnt5/aio7BTG00NAiwqo1OLYTcjWdcl+WxVXXifT3b371fVV5M8WFXvS/KtJHfvvpsAwHZtNOIlpF0dOw5Z3f2nSf7xlPb/m+Qdu+kUzMpG/yEBLKJpYcoU4tUzYgkHWFi+3QGwV2yrAwAwgJAFADCAkAUAMICQBQAwgMJ39hULjwJcmWUdrh4hi33FwqMAV2ZZh6vHdCEAwABCFgDAAKYLWVjqrwCYZ0IWC0v9FQDzTMgCgCXnisMxhCwAWHKuOBxDyAIALrPV0a2N6mONgglZAMAU00LS0ZVTU4PXpfWxRsHWCFkAwJYYndoe62QBAAwgZAEADGC6EAAYblqBfLK/i+SFLABguI0WkN7PRfKmCwEABjCSBQDsqY3W2Fo2QhYAsKf2a43VdglZzB2rBwOwHwhZzJ1pxZH7uTASgP1J4TsAwABGstg103sAcDkhi10zvQcAlzNdCAAwgJEsFoI1VwBYNEIWC0F9FwCXmveaYCGLITYaebr0L/5G/0AAYDPzXhMsZDHEtG8R0/7ib7RhKAAsOiELAJiZrc58LCIhCwCYma3OfCwiSzgAAAwgZAEADCBkAQAMoCaLq8aCogBsxbTfFxsdN8+ELK6a/XClCADj7ZffF0IW22LxUADYGiGLbbF4KABsjcJ3AIABhCwAgAGELACAAYQsAIABhCwAgAGELACAAYQsAIABrJO1hDZaUPTSFXYtPAoAOydkLaFpC4oeXTk1dV9BC48CwM4IWSTZP/tEAcC8UJMFADCAkAUAMMCw6cKquiPJh5O8Lslvd/fKqM/aL3ZbkL6V4y4cCwCMNSRkVdXrkvyXJP88yYtJvlpVD3X3syM+b7+YVpB+aTH6RscpXAeA+TJqJOvWJM93958mSVU9kOSuJELWIArXAWC+jApZh5K8sO7xi0nePuizdmWvp+h2+7mXOnTwwNQRKgBgvlV37/2bVv1ckju6+19PHr8nydu7+951xxxPcnzy8EeSfHPPO3K5a5L8+VX4HLbOOZk/zsl8cl7mj3Myn67Gefl73X3tZgeNGslaTXLjusc3TNq+r7tPJjk56POnqqonuvvI1fxMrsw5mT/OyXxyXuaPczKf5um8jFrC4atJbq6qv19VfyvJu5I8NOizAADmzpCRrO5+taruTfI/sraEw/3d/cyIzwIAmEfD1snq7keSPDLq/Xfoqk5PsiXOyfxxTuaT8zJ/nJP5NDfnZUjhOwDAsrOtDgDAAEsRsqrqjqr6ZlU9X1UnZt0fkqq6v6rOVtU3Zt0X1lTVjVX1+ap6tqqeqar3z7pPJFX1t6vqK1X1x5Pz8u9m3SfWVNXrqup/VdXvzbovJFV1uqq+XlVPV9UTs+5PsgTThZMtfv531m3xk+TdtviZrar6p0m+l+Tj3f2PZt0fkqq6Psn13f1UVb05yZNJfsa/ldmqqkryxu7+XlW9IckXk7y/u780464tvar6N0mOJPk73f3Ts+7Psquq00mOdPfcrF22DCNZ39/ip7v/X5ILW/wwQ939hSTfmXU/+IHuPtPdT03uv5Lkuazt3sAM9ZrvTR6+YfKzv78dL4CquiHJsSS/Peu+ML+WIWRN2+LHLw64gqo6nORtSb48256QfH9a6ukkZ5M82t3Oy+z9pyT/NsnfzLojfF8n+YOqenKyq8zMLUPIArahqt6U5NNJPtDd3511f0i6+7XuviVru2fcWlWm2Geoqn46ydnufnLWfeEiP9HdP57knUnumZSlzNQyhKxNt/gB1kxqfj6d5BPd/ZlZ94eLdfe5JJ9Pcses+7Lkjib5F5MaoAeS3F5V/322XaK7Vye3Z5N8NmvlQjO1DCHLFj+wBZMC648mea6775t1f1hTVddW1cHJ/QNZu4jnT2bbq+XW3R/s7hu6+3DWfqec6u5/OeNuLbWqeuPkgp1U1RuT/GSSmV+9vu9DVne/muTCFj/PJXnQFj+zV1WfSvJHSX6kql6sqvfNuk/kaJL3ZO1b+dOTnztn3SlyfZLPV9XXsval8dHutmQAXOy6JF+sqj9O8pUkD3f378+4T/t/CQcAgFnY9yNZAACzIGQBAAwgZAEADCBkAQAMIGQBAAwgZAEADCBkAQAMIGQBAAzw/wExNQRMRrn5iwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin, xmax = 0.0, 5.0\n",
    "\n",
    "fig2, ax2 = plt.subplots(figsize=(10, 6));\n",
    "counts, bin_edges, _ = ax2.hist(y_all, 100, range=(xmin, xmax), histtype='step');\n",
    "bin_centers = (bin_edges[1:] + bin_edges[:-1])/2\n",
    "s_counts = np.sqrt(counts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the distribution of \"whatever you put into it\" (initially x1-x2), which shows what output you get and what uncertainty to expect (given by the width - think about this!). We can thus get the result by simply recording the mean and width (RMS):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Mean = 2.700,    RMS = 0.206\n"
     ]
    }
   ],
   "source": [
    "print(f\"  Mean = {y_all.mean():5.3f},    RMS = {y_all.std(ddof=1):5.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, we are in principle not even sure, if this distribution is Gaussian, so in order to check this, we fit it. Note that in the second line (defining the Minuit object to minimised) we give \"reasonable\" starting values, which is key to making fitting work (more on that later)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fit value: N = 498.48699 +/- 4.99337\n",
      "Fit value: mu = 2.70075 +/- 0.00207\n",
      "Fit value: sigma = 0.20656 +/- 0.00154\n"
     ]
    }
   ],
   "source": [
    "Chi2_object = Chi2Regression(gaussian, bin_centers[counts>0], counts[counts>0], s_counts[counts>0])\n",
    "minuit = Minuit(Chi2_object, pedantic=False, N=5000, mu=2.5, sigma=1, limit_sigma=(0, 10000), print_level=0) #   \n",
    "minuit.migrad()                          # Perform the actual fit\n",
    "\n",
    "for name in minuit.parameters:\n",
    "    print(\"Fit value: {0} = {1:.5f} +/- {2:.5f}\".format(name, minuit.values[name], minuit.errors[name]))\n",
    "\n",
    "chi2 = minuit.fval\n",
    "N_var = 3                                # Number of variables (N, mu, sigma)\n",
    "N_dof = len(counts[counts>0]) - N_var    # Number of degrees of freedom\n",
    "\n",
    "from scipy import stats\n",
    "chi2_prob =  stats.chi2.sf(chi2, N_dof) # The chi2 probability given N_DOF degrees of freedom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And plot the fit on top of the histogram:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "xaxis = np.linspace(xmin, xmax, 1000)\n",
    "yaxis = gaussian(xaxis, *minuit.args)\n",
    "ax2.plot(xaxis, yaxis);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XlcVdX6+PHPZjpMiiAgCCIShinOpDmXessik0zrlhmSZZpzWqiVWuZ47V7zpjmGppYmPxyulpYpicNXUcNZc8IBMUEBI2Revz+AnScwkekAPu/Xa79gr7X2Ws8+oj7ss85amlIKIYQQQgghRB4zUwcghBBCCCFEZSIJshBCCCGEEHeQBFkIIYQQQog7SIIshBBCCCHEHSRBFkIIIYQQ4g6SIAshhBBCCHEHSZCFEEIIIYS4gyTIQgghhBBC3EESZCGEEEIIIe5gYeoA/o6zs7Py9vY2dRhCCCGEEKIaOHjwYKJSyuVe7Sp1guzt7c2BAwdMHYYQQgghhKgGNE27WJx2MsVCCCGEEEKIO0iCLIQQQgghxB0kQRZCCCGEEOIOkiALIYQQQghxB0mQhRBCCCGEuMM9E2RN077UNO26pmnH7ihz0jTtR03TzuR/dcwv1zRNm6tp2llN045omtbqjmuC89uf0TQtuHxuRwghhBBCiNIpzhPkZUCPv5SNA35SSjUEfso/B3gaaJh/DAK+gLyEGpgEtAXaAJMKkmohhBBCCCEqk3smyEqpncDNvxT3Apbnf78cCLqj/CuV5/+AWpqmuQNPAT8qpW4qpZKAHymcdAshhBBCCGFyJZ2DXEcpFZ///TWgTv73HsDlO9pdyS+7W7kQQgghhBCVSqk/pKeUUoAqg1gA0DRtkKZpBzRNO5CQkFBW3QohhBBCCFEsJU2Qf8ufOkH+1+v55XFAvTvaeeaX3a28EKXUIqVUgFIqwMXlnltlVxmPP/44w4YNM3UYQgghhBDiHkqaIG8EClaiCAY23FH+Wv5qFo8BKflTMbYCT2qa5pj/4bwn88vKzYABA9A0rdCxZMmS++pn2bJl2NvblzqeiIgIpk+fXup+xIMhKiqK5557jjp16qBpGpGRkYXarFu3Dj8/P6ytrQkICCA6OtqoPi4ujsDAQGxtbalbty4zZ84s1MdHH32Em5sbdnZ2BAUFcf36daP6Xbt20bJlS6ytrWncuDHff/99md6nEEIIURkVZ5m3b4C9gJ+maVc0TRsIzAD+oWnaGaB7/jnAd8B54CywGHgbQCl1E5gCROcfH+eXlavu3bsTHx9vdPTr16+8hy2Sk5MTNWrUMMnYouq5desW/v7+zJo1q8j6Y8eO8eKLL/Lqq69y6NAhmjVrxjPPPENKSorepk+fPty6dYs9e/bw6aefMnHiRFatWqXXL1y4kFmzZjF//nyioqKIi4sz+vuRkJBAYGAg7du359ChQwQFBfH8889z/vz58rtxIYQQojJQSlXao3Xr1qqkgoODVWBg4F3rd+zYoQC1Zs0a5efnp2rUqKGCg4NVVlaWUkqpsLCwgrnVRkeXLl2M+gkLC1N2dnZq9+7dqkWLFspgMChfX1+VmJiolFLqqaee0q8dOnRooTiys7PV5MmTlaenp7Kzs1OdOnVSMTExRm3Wrl2rmjRpoqytrZWbm5saMGBAiV8XUbUkJCQoQO3YscOofOTIkaply5b6eXp6urK3t1eLFy9WSin1yy+/KEAdPnxYb/PGG2+oDh066OfNmzdXo0eP1s937dqlAHXmzBmllFL/+c9/lLOzs8rOztbbNGzYUL3//vtleo9CCCFERQEOqGLkoA/8TnpLliwhPDyclStXsnLlSsLDwwF46aWXiI+PZ86cOdja2upPoCMiIgr1kZ2dzZgxY5g2bRrHjh1j/Pjxet0333xDfHw87dq1K3L8KVOmsHLlSsLCwoiJiaFjx4489dRT/P777wD89ttvvPLKK7zyyiucPHmS//3vf/j5+ZXDKyGqkujoaDp16qSfGwwG2rRpw4EDB/R6BwcHmjVrprfp3Lkzhw4dQilFRkYGR48eNeqjXbt2WFpaGvXRvn17zM3N9TadOnXS64UQQojqqlonyFu2bMHe3t7oOHr0qFGbiRMn4u/vz3PPPUebNm30eZw2Nja4ubnh4OCApmm4ubnh5uaGk5NToXEyMjKYOnUqTz/9NL6+vrz++uvUrl0bAEdHR9zc3LCysip0XXp6OrNmzWLu3Ll0794dX19fpk6dilKKzZs3A3DlyhWysrJ47rnn8Pb2JiAggHHjxhXqSzxYEhIScHZ2Zv369Tg6OnLgwAGcnZ0pWPklISGB2rVrk5aWRoMGDXjvvfdwdnbm9u3bpKamcuPGDXJzc3F2dmbs2LH4+PiQnp6Oo6OjUR/Ozs7s37+fWrVqsWHDBqMxhBBCiOrKwtQBlKfOnTuzaNEiozIvLy+jc19fX/17Jycnbt68/6nRmqbRsWPH+77uzJkz3L59mz59+qBpml5++/ZtfZ5ns2bNaN++PZ06deLJJ5+kffv2vPzyy7i6ut73eKL6qVGjBl5eXtja2hZZb25uTr169f7258XFxQUvLy8sLIr+58DOzo769evj4OBQJjGLyq3DjO3EJd8uVO5Ry4bd47qaICIhhKh41TpBtrW1NUqAi/LXpCBvesr9j1PUE+Li2rRpE/Xq1TMqK3hSbWlpSVRUFHv27CEyMpL58+czc+ZMjh8/jqOj7Nb9oHJxcSExMZFu3bpx+PBhABITE2nYsKFef+PGDQwGAzt37gRgxYoV2NjYYG9vj5WVFWZmZiQmJhIaGkpoaCi5ubkkJSVRsLxiwRhNmjTRx1ixYgXVaflFUVhc8m1iZwQWKvcet9kE0QghhGlU6ykWZcHKyors7Oxy6bthw4ZYW1sTHx+Pr6+v0XHnVA4zMzM6duzIBx98QFRUFPHx8ezfv79cYhJVQ0BAAFFRUfp5RkYG+/fvJyAgQK9PSUnhyJEjepudO3fSqlUrNE3DYDDg7+9v1MfevXvJysoy6mPPnj3k5OTobaKiovR6IYQQorqq1glyRkYG165dMzpSU1Pvqw9fX18yMjLYsGEDaWlpZGZmFvvazMxMfdzMzEzS0tL0cwBra2veffddxowZw9q1azl//jyRkZEMHjyYU6dOAXkflPrkk084ePAgsbGxLFq0CEtLSxo1anRf9yGqltTUVGJiYjh27BgAZ8+eJSYmRp8C9Oabb3L06FE++eQTTpw4wZAhQzAYDPTt2xeAli1b8thjjzF06FAOHz7M6tWr+eqrrxgyZIg+xpAhQ1i4cCEREREcOnSIUaNG6XPhAfr160dmZiYjRozg5MmTTJgwgUuXLhESElLBr4YQQghRwYqz1IWpjtIu80YRy7QVLFFVsMxbQkKCfk1gYKAKDg4u1Nfo0aOVi4vL3y7zVpSCMYo6CmRnZ6uPP/5YeXt7K0tLS+Xl5aVCQkL0uE6fPq169OihnJ2dlY2NjWrevLlat25diV8XUTXc7WcnLCxMb/P//t//Uw0bNlRWVlaqdevWat++fUZ9XL58WT399NP68oDTp08vNM6kSZOUq6ursrGxUb169VLXrl0zqt+5c6dq3ry5srKyUo0aNVKbN28ul/sVlUf90E33VS6EEFUJxVzmTVMlmHNbUQICApQsKSWEEBXHe9zmwnOQ027S5OPtHJ/RxzRBCSFEGdE07aBS6p5zBav1FAshhBClkHYTVveDWQ2IMQyC78dBTvl8JkMIISoTSZCFEEIUlp0JX78IZ36ATmMJz+kM+76AzaNNHZkQQpQ7SZCFEEIUFvUpXImG3oug24eMz34TOo2BQ1/B6S2mjk4IIcqVJMhCCCGM/ZEIe/4LjYOgyfN/lj8+AZwegp8+gtycu18vhBBVnCTIQgghjO2ZC1lp8MQE43JzC+j6Plw/AcciTBObEEJUAEmQhRBC6AxkwoFleU+OXfwKN2j8PNRuCNFLKjw2IYSoKNUyQV62bBmapvH883++NRgTE4OmaXh7e5suMCGKITs7mzFjxuDn54etrS1eXl6MHTuWtLS0YvdR8Hfgr4ednZ1Ru48++gg3Nzfs7OwICgri+vXrRvW7du2iZcuWWFtb07hxY77//vsyuUdReT1ldgAyUiDgLhvCmJlBq/5w+f8g4deKDU4IISpItUyQAWxsbDhw4AC3bt0CYM2aNdSrV8/EUQlxbxkZGRw9epRp06Zx+PBhwsLC+Pbbbxk+fHix+3jppZeIj483Orp168YLL7ygt1m4cCGzZs1i/vz5REVFERcXR79+/fT6hIQEAgMDad++PYcOHSIoKIjnn3+e8+fPl+n9isqlj/nP4OAF9TvevVHzl8HMAn5ZUXGBCSFEBaq2CbKZmRmBgYFs2LABgPDwcPr0MV7k/vDhw3Tr1g1bW1vq16/PxIkTyc7+c43Pf//73zRr1gw7Oztq167Na6+9pm/1CzB58mT8/f357LPPqFOnDq6urnz22WcVc4Oi2rKzs+OHH37ghRdeoGHDhnTr1o3Ro0cTEVH8OZ82Nja4ubnpR25uLpGRkbz++ut6my+++IK33nqL3r1706pVK+bMmcO2bds4e/YsAKtWrcLKyoq5c+fSuHFjpk2bhpeXF19++WWZ37OoJFKv09HsGDT/Z96T4ruxd4WHusHx9VCJN5sSQoiSqrYJMuQ9RVuzZg3R0dG4u7tTt25dve7GjRt07dqVgIAADh8+zIoVK1i1ahX//ve/9TZJSUlMmTKFw4cP8/3333Ps2DGGDBliNMb58+c5c+YMUVFRDBkyhDFjxnDlypUKu0fxYEhOTsbR0bHE14eFheHt7U2XLl2AP59Sd+rUSW/Trl07LC0tKdi9Mjo6mvbt22Nubq636dSpE7K7ZTX261bMNAWP9CxU5VHLBu9xm/Uj9IQXpFzitelhJghUCCHKV7VOkLt06cLRo0dZsGABL730klHd559/jq+vLzNnzqRhw4Z07tyZMWPGGD0dmzJlCr169cLX15c2bdowfPjwQnMwzczM+PTTT3n44YeZMGECubm5/PLLLxVyf+LBcPXqVT7//HPGjBlTouuVUixdupSQkBA0TQPyfkHMzc3F2dmZsWPH4uPjQ3p6Oo6OjiQkJAB5UyycnZ3Zv38/tWrVYsOGDTg7O+v1ohr6dQtxqja4NS1UtXtcV2JnBOrHzPHvARrN/9hd8XEKIUQ5szB1AOXJzMyMXr16MX/+fKZNm8aKFX/Olzty5AiHDh3C3t5eL8vJySEn58+1Pbdv3860adM4efIkKSkp5OTkkJGRYTSGp6cnBoMBAIPBgK2trdE0DCFKIzU1laCgIHr06MHQoUNL1Me2bdu4ePEiwcHBRda7uLjg5eWFhUXR/xzY2dlRv359HBwcSjS+qCKy0uHcdrbntKd//i9Sf8veFTwf5R+XDpZ/bEIIUcGq9RNkgBEjRhAWFkadOnUK1T377LPExMTox9GjRzlx4gQAFy9eJDAwED8/PzZs2EBMTAxTp05F/WW+XVFJxV/bCFESt2/fpmfPnri6urJs2bIS97N48WKefPJJPD099bLatWtjZmZGYmIioaGhREZGYmFhQVJSEi4uLkBe4pyYmEiTJk04fPgwjz/+OImJiXq9qGYu7oasNH7KbVn8ax5+kmZmFyBNHgoIIaqXap8g+/r60r9//0LlTZs25dSpU/j4+ODr62t0ABw4cIDMzEzmzp1LQEAAvr6+XL16taLDFw+ojIwMnn/+eSwsLAgPD8fS0rLIdqmpqcTGxpKamlpkfUJCAhs2bDD6cB7kvdvh7+9PVFSUXrZ3716ysrIICAgAICAggD179hi9qxIVFaXXi2omNgrMLNif+0jxr2nQ5c9rhRCiGqn2CfLdDBs2jOvXr/P6669z+PBhTpw4wZIlS5g4cSIADRs2JDc3l/nz53PhwgWWL19uNEVDiPKSlZVFnz59+O2331i0aBHJyclcu3aNa9euGSWrkLc6S4MGDQgPDy+yr+XLl1OjRg169epVqG7IkCEsXLiQiIgIDh06xKhRo+jevbv+S2K/fv3IzMxkxIgRnDx5kgkTJnDp0iVCQu6yPq6o2i5EgUdr0rAu/jV1W5KqrOHCzvKLSwghTOCBTZCdnZ356aefiI+Pp0OHDrRr146wsDAaN24MQLNmzfjss8+YOXMmTZo04dtvv+XDDz80cdTiQRAXF8emTZuIiYnBx8cHd3d3/bh8+fJ99bVkyRL69euHlZVVobrBgwfz7rvvMmTIEDp27IiHhwcrV67U611dXdm0aRO7d++mRYsWrFu3joiICB566KFS36OoZDJ+h6u/gHene7e9k7kl0bl+kiALIaodrTLPlw0ICFCypJQQQpSzMz/Cqj7Qfz3ei9OInRFY7Eunvj+E9y2/hjGnoYZbOQYphBClp2naQaXUPecKVutVLIQQQhTDxT15O+PVawvsuL9L7ZpB5te8Nf0Ltua20cs9atmwe1zXMg5UCCEqhiTIQgjxoLsSDXX8wcr2vi9d9N7rMH0iCx/PhSf/fPLsPW5zWUYohBAV6oGdgyyEEALIzcmbf+z5aMmutzCAe/O8JFsIIaoJSZCFEOJBlnAKMlNLniADeLbJS7KzM8suLiGEMCFJkO/DsmXLjHbeE0KIKq/gya9nKda39gyA7HT47VjZxCSEECZWbRPk27dvM27cOLy9vbG2tsbHx4eBAwcWavf4448zbNgwE0QoxN1lZWUxbNgwnJycqFmzJiEhIfzxxx/Fvv7ixYv0798fT09PbGxsaNy4MYsWLSrU7qOPPsLNzQ07OzuCgoK4fv26Uf2VK1d46aWXcHR0xN7enscee4ykpKRS35+oRK5Eg40TOPmUvI+C5PrqobKJSQghTKzaJsgjR45k/fr1LF26lJMnT7JgwQKys7NNHZYQxTJx4kTWrl1LeHg4W7ZsITIyklGjRhX7+jNnzmBlZcXy5cs5fvw4oaGhDB061Gizm4ULFzJr1izmz59PVFQUcXFx9OvXT69PT0+nW7duJCQksHnzZo4cOcKECROK3F5dVGHxh6FuC9C0kvfhUA+sa0H8kbKLSwghTEkpVWmP1q1bq5JycnJSS5cuvWt9ly5dFFDoCAsL09ucOXNGderUSRkMBhUQEKDGjx+v7OzsShyTEMWRnZ2tnJyc1GeffaaXrVy5UhkMBpWamlrifp977jkVFBSknzdv3lyNHj1aP9+1a5cC1JkzZ5RSSoWFhalatWqpW7dulXhMUcllZSj1UW2lfvhQL6ofuqlkfS17VqmFXUrfjxBClCPggCpGDlptnyDXrFmTH374gfT09CLrIyIiiI+Pp127doSEhBAfH098fDwvvfSS3ubVV1/FxsaGgwcPEhoayueff15R4YsH2Pnz57l58yadOv25q1nnzp3JyMjg6NGjJe43OTkZR0dHAL2vO8do164dlpaWFGzOExkZSceOHZk0aRJubm40bdpU/g5UN4mnITcL3JqVvi+3ZvDbCcjJKn1fQghhYtU2Qf7888/54YcfcHFxITAwkHnz5pGSkqLXOzk54ebmhpWVFba2tri5ueHm5oaNjQ0AR48eZd++fcyZM4cmTZrQp08fXn75ZVPdjniAJCQkAHnbofft25e2bdvi7OxsVHe/tm7dyt69exkxYgQAN27cIDc3F2dnZ8aOHYuPjw/p6ek4OjrqY8THx7Nz506uXbvG5s2bGTVqFO+88w7h4eFlcJeiUriW/wtXWSTI7s0hJwMSz5S+LyGEMLFqmyAHBgZy6dIlli9fzkMPPcT06dNp1qxZsROMs2fPYm5ujp+fn17WtGnT8gpXiCK5u7vj5eVVqj5OnjxJv379+Pzzz2nRokWhehcXF7y8vArNLc7NzSU3N5clS5bQunVrBg4cyPPPP8/KlStLFY+oROKPgIUN1H6o9H255f/7eE3mIQshqr5qmyAD2Nvb07t3b+bOncuJEydIT08nLCzM1GEJ8bdcXFwASExMZO7cuaxdu5bExESjuuI6d+4c3bt3Z8yYMQwaNEgvr127NmZmZiQmJhIaGkpkZCQWFhYkJSXpYzg7O+Pp6Ymt7Z+7q/n4+BAXF1faWxSVxbWjUKcJmJmXvq/aDcHCWj6oJ4SoFqp1gnynmjVr4u7uzu+//25UbmVlVeTqFr6+vuTk5HD69Gm97NgxWeNTlD8fHx8cHR2JiorSy3bu3InBYCj0LkZqaiqxsbGkpqYW6ufSpUt07dqVgQMHMn78eKM6g8GAv7+/0Rh79+4lKyuLgIC8JbtatmxJXFwcGRkZRn16eHiUyX0KE1MqL0F2K6N3xswt8pJteYIshKgGqm2C/Nxzz7FkyRKOHj3Kr7/+ypQpUzh69ChPPfWUUTtfX1927tzJpUuXSE9PJycnB8ibTtGmTRtGjRrF8ePHiYiI4JtvvjHFrYgHjLm5OYMGDWLq1Kls376dPXv28MEHH9C/f3/s7OyM2oaHh9OgQYNC84Lj4uJ44okn6Nq1K2+//TbXrl3j2rVrRlOMhgwZwsKFC4mIiODQoUOMGjWK7t274+vrC8Arr7xCTk4Ow4cP59dff2Xt2rVEREQYLQUnqrDkS5CRAu5lMP+4gFuzvARZqbLrUwghTKDaLmjarl07Pv/8c0aPHo2mafj5+bFmzRo6duxo1O7dd9/l2LFjPPLII6SlpREWFsaAAQMAWLVqFSEhIbRq1YqmTZsybNgwPvvsMxPcjXjQfPzxx6SmpvLCCy+Qk5ND7969mTNnTrGv//HHHzl//jznz59n2bJlenn9+vWJjY0FYPDgwVy7do0hQ4bw+++/8+STT7Jw4UK9raenJ//73/8YM2YMzZo1w9PTkxkzZtC3b9+yuk1hSgVPesviA3oF3JvBwbC85FsIIaowTVXi3/QDAgJUwZJTQgghylDkjLxjQhxY/fnOhPe4zcTOCCxZn5f+D758Cl75Fu8vs0vejxBClBNN0w4qpQLu1a7aTrEQQgjxNxJOgWN9o+S41Fz8/uxbCCGqMEmQhRDiQZRwGlwalW2fNo5g75bXtxBCVGHVdg6yEEKIu8jJJuv6ryy9+hAzxm02qvKoZVO6vl384PpJQKZXCCGqLkmQhRCiGuswYztxybeNyny0q2w3ZDO4zzMMblHGiazrI3BoBVB5P98ihBD38sAmyJGRkTzxxBMkJCTo2/jezeTJkwkPD5d1kIUQVU5c8u3CH5Y7+T9Yw59zhsuSix9k/YEHiWXftxBCVJBqPQf57Nmz9O7dm1q1amFvb0+HDh04dOjQffczduxYfv7550Llixcv5tFHH6VmzZq4uLjQt29ffQktIe7m9OnTPPHEE7i6umJra0uzZs0Kbd88ePBgHnroIaytrfHw8ODtt9/m1q1b9zXO7t276dChAzVq1MDNzY2QkBCSk5P1+uTkZAYOHIibmxv29va0b9+eXbt26fU///wzTz31FC4uLtSoUYOOHTsW+fdAVEEFH6Jzfrjs+3Z5BICGZrLjohCi6qq2CXJ8fDzt27cnKSmJTZs2ER0dzT//+U/i4+Pvuy97e3tq165dqHzXrl288cYb7Nmzh23btnHz5k169OhR5M58QhSwsLDg1VdfZdu2bZw4cYKRI0fy+uuvs23bNr1N8+bNWb58OadOnWLt2rVERUXx9ttvF3uMW7duERgYSNOmTYmJiWHdunXs3r2bUaNG6W3eeecddu/ezbp164iJicHf35/AwEBSUlIA2LdvH+3atWPTpk388ssvtGnThqeffppff/217F4MYRoJp8GhHhhqlH3f+U+lfTVJkIUQVZhSqtIerVu3ViU1fPhwVbNmTZWamlpk/Y4dOxSg1qxZo/z8/FSNGjVUcHCwysrK0ttMnz5dkTeRTjVp0uSeYx46dEgBKiYmpsRxiwdTq1at1Pjx4+9aP3fuXOXu7l7s/qKjoxWgYmNj9bJJkyYpf39//bxJkyZq8uTJ+vmFCxcUoA4cOFBkn7m5ucrBwUHNmTOn2HEI06sfuqlw4RcdlFrRu/wGneWr1rzfs/z6F0KIEgIOqGLkoNX2CfKWLVt4+umnC23N+1dLliwhPDyclStXsnLlSqMte4cPH058fDxjxowp1pgFb1/XqlWr5IGLB4pSih9//JETJ07QunXrIttcvXqVb7/99q71RWnUqBGurq6sXr2anJwcEhMT2bJlC4GBf85F7dy5M1u3biUxMZGcnBy++eYbvL29ady4cZF93r59m8zMTBwdHe/vJkXlkpsDiWfKfom3O7k2kikWQogqrdomyJcuXcLT0/Oe7SZOnIi/vz/PPfccbdq0ITo6Wq+zs7PT52feS05ODhMnTqRPnz7Ur1+/VLGLB0P79u0xGAw8++yzzJkzhxdeeMGofv78+dja2uLh4YGDgwNff/11sfu2t7cnMjKSpUuXYjAYcHFxwc/Pj+nTp+tt5syZwyOPPIKLiwsGg4ElS5awdetWbGyKXuZrypQpuLq60qdPn5LdsKgcki9Cdnr5JsgujfKmWFTinVqFEOLvVNsEubh8fX31752cnLh582aJ+hk+fDg3btxg8eLFZRWaqObWrFnDwYMHmT59OhMmTGDPnj1G9f369SMmJob//e9/nD9/no8++qjYfaenp/P666/TrVs3oqOj2bp1K9HR0YwbN05vM2/ePPbv388PP/xAdHQ03bp1IzAwsMgPA65evZp58+YRERGBra1tyW9amF5C/hzy8viAXoHaDamh3YbU6+U3hhBClKNqu8ybl5cXcXH3fovPwsL4JVAleOIRGhrKli1biIqKkukVotjq1atHvXr1aNq0KcePH+eTTz7hu+++0+sdHBxwcHDg4YcfxsHBgS5dujB27Fjc3Nzu2ffq1as5d+4cu3fvxsws7/fgmTNn8vzzzzN58mTMzMwYP3484eHh/OMf/wBgwYIF1KlThzVr1vDmm2/qfW3cuJFBgwaxfv16WrVqVcavgqhwN87mfa3t+/ftSqP2Q3+OVaNO+Y0jhBDlpNo+QX7qqaf4/vvvSUtLK9dxJk+ezDfffMP27dvx8PAo17FE9WVubs4ff/zxt/VKKW7fNt7wITU1ldjYWFJTU43Kk5OT0TQNTdP0MgsLC3JycsjMzCQtLY04s3ifAAAgAElEQVSMjAw9eQbQNA1zc3OjMbZu3cqrr77KmjVr6Nq1a2lvU1QGN8+BtQPYOpXfGAXJd0EyLoQQVUy1TZDHjRuHpaUlPXv2ZM+ePZw6dYp58+axefPme1+c79q1a1y7do3U1FSys7P188zMTABmzJjBZ599xtdff421tbVe/9ckRog7hYWFsXTpUo4ePcq5c+dYvHgxX331Fb179wbg3Llz+pSLixcvEhkZyfDhw2nRogUNGjQw6is8PJwGDRoYfbgUoFu3biQlJTF27FjOnDnDvn37+OCDD+jcuTMODg44OjrStm1bPvjgA/bu3cvZs2cJDQ0lOTmZJ598EoAdO3bQu3dvZs+eTcuWLfWf74Jl4EQVdeMcOD0Ed/zyVOYcPMlQlpIgCyGqrGqbIHt4eLBr1y5q1KjB008/TevWrVmxYgXu7u7F7sPd3R13d3c+/fRTTp8+rZ8XzBVdsGABycnJdOrUSa9zd3dnzZo15XVbohqws7Pjiy++oGPHjvj7+/Of//yH6dOnM2LECABsbGw4fPgwL7zwAg8//DD9+vXD39+fTZs2FXuMpk2bsm7dOqKiomjZsiW9evWicePGRj+b4eHh+Pn5ERQURIsWLdi5cycbN26kUaO8D28tX76ctLQ03nrrLaOf75EjR5btCyIq1s3zf06BKC9m5sSqOnnJuBBCVEFaSebcVpSAgAB14MABU4chhBBVlve4zX9uNZ2VDlPd4PFxeUc52vJhd3rUuQXD9pfrOEIIcT80TTuolAq4V7tq+wRZCCHEXyRdAFTeFItydkG55Y2Xm1PuYwkhRFmTBFkIIR4UBVMeavuU+1AXlBvkZELK5XIfSwghypokyEII8aC4mZ8gV8QT5Nz8z3vIB/WEEFWQJMhCCPGguHEObGuDTfmv135BFSTI8kE9IUTVIwmyEEI8KG6er5CnxwCJ1ARDTXmCLISokqptgqxpGs7OzqSnp+tl/v7+TJ48GYDHH38cTdOwtLTEy8uLN998k6tXrxr1MWDAAH2zhTuPv645K0RVNHXqVJo2bYqdnR1169bljTfeIDEx0ajN0qVLeeSRR7CxsaFBgwZMmzbNRNGKMnHjXPkv8abT8sZKPFNB4wkhRNmptgkyQFJS0t+uSRwSEsKFCxdYtmwZp0+fpl27diQlJRm16d69O/Hx8UZHz549yzt0Icrd3r17CQ0N5eDBg0RERHDgwAFeeuklvf7nn39m0KBBvPPOO5w8eZLZs2czZcoUli1bZrqgRcllpsHvVyvsCTIATj6QFFtx4wkhRBmp1glyYGAg8+fPv2u9ra0tnp6edO3alY0bN5KYmMgXX3xh1MZgMODm5mZ0GAyG8g5diHK3adMmXn31VRo1asRjjz3GpEmT2L59u75TXnR0tP7uire3Ny+88AJt27bl4MGDJo5clMjN83lfK2AFC52jd94qFjnZFTemEEKUgWqdIL/44oucO3euWP+h16pVi7Zt2/LTTz9VQGRCVD7JycnY2NjovwB27NiRq1ev8vPPPwNw9OhRjh49ytNPP23KMEVJVeAKFjpHb8jNhltXKm5MIYQoA9U6Qba2tmbgwIHMmzevWO3d3NwKzUPesmUL9vb2+vHQQxX4n4sQFSQ1NZXp06czdOhQrK2tAXjsscdYvXo1PXv2xNLSklatWjF16lSeeeYZE0crSkRfA7kiE+QGeV9lmoUQooqxKM3FmqaNBt4AFHAUCAHcgdVAbeAg0F8plalpmgH4CmgN3ABeUkrFlmb84hg8eDBNmzbl008/vWdbTdP469bbnTt3ZtGiRfq5hUWpXjIhKp3s7GxeeeUV6taty9SpU/XyU6dOMWLECKZPn07nzp2JiYlh5MiR1K1bl+eee86EEYsSuXkO7FzBUKPixnT0zvsqCbIQooopcbanaZoHMAJorJS6rWnat8A/gWeA/yilVmuatgAYCHyR/zVJKeWrado/gZnAS3fpvsw0aNCALl268OWXX96zbXx8PB4eHkZltra2+Pr6lld4QphUbm4ur732GvHx8fz0009YWVnpdTNmzKBt27YMHToUgKZNmxITE8PMmTMlQa6KbsbmfWiuItWsC2aWcPNCxY4rhBClVNopFhaAjaZpFoAtEA90BQrWQVsOBOV/3yv/nPz6bpqmaaUcv1iGDh3KggULCj0dvtPNmzf5v//7P7p27VoRIQlhckop3njjDU6cOMHWrVupWbOmUX1ycjJmZsb/RFhYWHD79u2KDFOUlaRYcKxfsWOamUMtL3mCLISockqcICul4oDZwCXyEuMU8qZUJCulCj6yfAUoeCTrAVzOvzY7v33tv/aradogTdMOaJp2ICEhoaThGenRowe5ubmcOnXKqDwtLY0rV66wfft2nnvuOVxcXHj77bfLZEwhKrvBgwfz008/sXLlSjIzM7l27RrXrl0jMzMTyPt7s379epYvX86FCxfYuHEjS5YskWUOq6LsTLgV9+eUh4rk1EASZCFElVPiBFnTNEfyngo3AOoCdkCP0gaklFqklApQSgW4uLiUtjsAzMzMGDx4MLm5uUblYWFheHt789prr9GoUSP27t2Lo6NjmYwpRGW3aNEiLl26RNOmTXF3d9ePPXv2AHkJ9NSpU5k6dSqNGzdm5MiRDB06lA8//NDEkYv7lnIZUFCrgp8gQ15SniRTLIQQVUtpPnHWHbiglEoA0DQtAugA1NI0zSL/KbEnEJffPg6oB1zJn5LhQN6H9crFX6dTvPvuu7z77rv6eWRk5D37kA0RRHX2d1OOCvz1742oogqe4Fb0FAvIS5DTU+B2EtjIAwghRNVQmjnIl4DHNE2zzZ9L3A04AewA+uS3CQY25H+/Mf+c/Prtqjj/QwshhCid5It5X00xxUKWehNCVEGlmYO8j7wP2x0ib4k3M2AREAq8o2naWfLmGC/Nv2QpUDu//B1gXCniFkIIUVxJsXmrSdRwr/ixC5JyWclCCFGFlGpRX6XUJGDSX4rPA22KaJsO9C3NeEIIIUog6SLUqpe3qkRFK5jWIU+QhRBVSLXeSU8IIQR5UyxMMb0C8jYmsXORBFkIUaVIglwKkydPxt/f39RhCCHE30uKNc0KFgUcvSVBFkJUKdU2QR4wYACapqFpGgaDAT8/P6ZPn05OTo6pQxPinrKyshg2bBhOTk7UrFmTkJAQ/vjjj2Jff/HiRfr374+npyc2NjY0btzYaMv0Ah999BFubm7Y2dkRFBTE9evX9brIyEj971DB4ezsXCb3JyqOPWl5K0iYYgWLArLUmxCiiqm2CTJA9+7diY+P5/Tp0wwdOpT333+f2bNnmzosIe5p4sSJrF27lvDwcLZs2UJkZCSjRo0q9vVnzpzBysqK5cuXc/z4cUJDQxk6dCgrVqzQ2yxcuJBZs2Yxf/58oqKiiIuLo1+/foX6On78OPHx8cTHx3PixIkyuT9Rcepp+RsumWqKRcHYKVcgJ8t0MQghxH2o1gmywWDAzc0Nb29vRowYQbdu3Vi/fr1eHxsbi6ZpbNq0iWeeeQZbW1vc3Nz46aefALh06RLPPvssdnZ2ODo6EhISQmpqaqFx3nvvPWrWrImHhwdLliypsPsT1VNOTg6LFi3i/fffp2vXrrRv355PPvmEFStWFPspcvfu3Vm6dCndunXDx8eH4OBgnnnmGSIiIvQ2X3zxBW+99Ra9e/emVatWzJkzh23btnH27FmjvlxdXXFzc8PNzQ1XV9cyvVdR/ry0/HcFTDrFogGoXEi+ZLoYhBDiPlTrBPmvbGxsyMoq/ARj7NixBAUFceTIEZYtW4adnR0Ar7zyCsnJyezZs4eNGzeyc+fOQpsmnDp1ivj4ePbv309oaCiDBg3i2LFjFXI/ono6f/48N2/epFOnTnpZ586dycjI4OjRoyXuNzk5Wd8psqCvO8do164dlpaWHDhwwOi6Rx99lLp16xIYGCg/21WQZ0GCbMonyLW88r5KgiyEqCIeiAQ5NzeX7777ji1bttC9e/dC9S+++CKDBg3C19eXHj168Nhjj3Hs2DF2797NvHnzaN68OZ06deKjjz7iyy+/JDMzU7/WzMyM//73vzRq1IgRI0YQEBDAl19+WZG3J6qZhIS8t8SdnZ3p27cvbdu21ef+FtTdr61bt7J3715GjBgBwI0bN8jNzcXZ2ZmxY8fi4+NDeno6jo6O+hju7u4sWbKEdevW8c0335Cbm0vnzp357bffyuAuRUWppyWAoaZpd7GrVS/va8pl08UghBD3oVTrIFd2W7Zswd7enszMTDRNIzg4mMmTJxdqd+dTtAJnz57FzMzMaJWK5s2bk5mZyaVLl/D19QXAw8ODWrVq6W2aNGnCuXPnyv5mxAPJ3b30GzucPHmSfv368fnnn9OiRYtC9S4uLnh5eWFhYfzPgZ+fH35+fvp5mzZtaNCgAV999ZVsP12FeGnX86ZXaJrpgqjpAZoZJEuCLISoGqp1gty5c2cWLVqEwWCgbt26mJsXvUh+wdvOZUUz5X9EospzcXEBIDExkblz5wJw+fJlo7riOnfuHN27d2fMmDEMGjRIL69duzZmZmYkJiYSGhpKaGgoubm5JCUl3XWMgtUwYmNjS3BXwlTqaQng2NK0QZhbQo26MsVCCFFlVOspFra2tvj6+lKvXr27Jsd389BDD5Gbm2s05/Lw4cNYWlri5eWll8XFxZGcnKyfHz9+XH+6LERJ+Pj44OjoSFRUlF62c+dODAYDTZs2NWqbmppKbGxskR8evXTpEl27dmXgwIGMHz/eqM5gMODv7280xt69e8nKyiIgIKDIuLKzszl79ize3t6luDtRoZSinnbdJPOPPWrZ4D1us37sS7ZnX0wMHWZsr/BYhBDiflXrJ8il0bRpU9q1a8ewYcP4/PPPuXXrFpMmTSIkJAQrKyu9XW5uLiNHjmTChAn88MMPHDhwgGXLlpkucFHlmZubM2jQIKZOnYq/vz/W1tZ88MEH9O/fX/8AaYHw8HBCQkIICwtjwIABenlcXBxPPPEEXbt25e233+batWt63wVPiIcMGcKYMWPo2LEj3t7ejBo1iu7du+u/4H366ad4e3vTrFkz0tLSmD17Nr///juvvvpqxbwQovRSr2OtZZlkBYvd47oaF0RsgIt7ifvtdoXHIoQQ90sS5L/x9ddf8/bbb9OuXTusrKzo1asX//rXv4zaNGrUiDp16hAQEEDNmjVZunQpjRs3NlHEorr4+OOPSU1N5YUXXiAnJ4fevXszZ86cYl//448/cv78ec6fP2/0C1v9+vX1KRKDBw/m2rVrDBkyhN9//50nn3yShQsX6m2zsrJ45513+O2337C3t+fRRx8lMjKyTOZFiwqSciXvq4OnaeMAcKgHt8IxRzZrEkJUfppSytQx3FVAQID665JTQgghiun4elgbDG9FgXsz08ZycDn8bwQd0j9j94wBpo1FCPHA0jTtoFKq6LmEd6jWc5CFEOKBdisu72tleIKcv9Sbh5Zo4kCEEOLeJEEWQojqKuUKacpg2jWQCzjkfbjZUyvZWt5CCFGRJEEWQojqKuUK8crJtGsgF8h/ii1PkIUQVYEkyEIIUV2lXCFOOZs6ijyW1mBfB09JkIUQVYAkyEIIUV3diiNe1TZ1FH9yqIeHTLEQQlQBkiALIUR1lJ0Bqb9xtTIlyLW8ZIqFEKJKqLYJsqZpODs7k56erpf5+/szefJk0wVVBURHR2Npacmzzz5rVL5//366dOlCzZo1cXFxYdiwYWRkZOj1mZmZjBkzBk9PT2xtbWnRogUbNmwocRxBQUFomsZfl/mbPn06Pj4+2NjY0KhRI5YsWaLXLVu2DE3TCh1/3VyjrC1YsIAGDRqQnZ1dLv3f7bW4ly1bthAQEICdnR0eHh5MnjyZO5d1jI+P5+WXX8bV1RV7e3t69uypb2kNsGvXLjp06EDt2rWxt7enbdu2bNmypczuS5SzW1cBuEplSpDrUVe7Abm5po5ECCH+VrVNkAGSkpJYs2aNqcOoMtLS0hg4cCCtW7c2Kk9KSiIwMJCHH36YgwcPsmbNGjZu3Mi4ceP0NjNmzGDFihUsX76c48ePExQURJ8+fThz5sx9x7FkyRKSkpIKlX/11Vd8/PHHzJ49m5MnTzJ69GgGDRrEjh07AHjppZeIj483Orp168YLL7xw3zEUV05ODrNnz2bMmDFYWJT9vjt3ey3u5ezZswQFBfHMM89w5MgR5s+fz5w5c5g7d67epn///ly4cIEtW7awd+9e0tLSCAoK0pNoa2trhg0bxs6dOzl69ChBQUH06tWLEydOlNn9iXKUv0lIZXuCbNCyIfU3U0cihBB/TylVaY/WrVurkgJUz549VZs2bfSyJk2aqEmTJpW4z+pu8ODB6pNPPlHBwcEqMDBQL//uu++UpmkqNTVVL5s7d65ycHBQWVlZSimlAgMD1YABA/T67OxsZW5ursLDw+8rhrNnz6qHHnpIHTx4UAEqOjparxs6dKh6/PHHjdrXr19fzZ49u8i+4uLilLm5udqxY8d9xVAAUBcuXPjbNqtXr1YuLi4qLS2tRGP8nb97Le5l/vz5ytXV1ajsnXfeUU2aNFFKKfXHH38oTdPU5s2b9fojR47ccxwnJye1cOHC+7wTYRIx3yg1qaZ6YtwiU0fyp9NblZpUU6lL+0wdiRDiAQUcUMXIQav1E+QXX3yRc+fOcfDgQVOHUul999137Nu3j9DQ0EJ1mZmZaJqGlZWVXmZtbU1KSgoXLlwAoHPnzkRFRXHp0iWUUqxevZoaNWrQvn37YseQk5ND//79mTFjBk5OToXqO3fuzOHDhzl+/DgAO3bsIDExke7duxfZX1hYGN7e3nTp0qXYMdyvmTNnMnz4cGxsbMq033u9FveSmZlp9OcFeX9mJ0+eJD09naysLJRSGAwGo3qAX375pch4Vq1aRXJyMi1btrzveIQJpORNl6lcT5DzNgsh+ZJp4xBCiHuo1gmytbU1AwcOZN68eaYOpVJLSEjgrbfeYunSpUVOE2jTpg0Gg4F//etfZGdnc+XKFRYuXKhfC/Duu+/y8ssv4+3tjZWVFaNHj+a7777D3d292HFMmzaNunXr0qdPnyLrX3zxRaZOnUqrVq2wtLSkZ8+erFq1iubNmxdqq5Ri6dKlhISEoJXTGrBbt27lzJkzDBs2rMz7vtdrcS+dO3cmLi6OlStXkpuby8mTJ1m1ahW5ubncvHkTBwcHmjdvzty5c7l16xZ//PEHU6ZMwcLCQv8zLeDp6YnBYGDIkCFERETw6KOPlsUtivKWEgc2TqRjuHfbiuIgCbIQomqo1gkywODBg/n2229LNI/zQfHWW28RHBx81yeD7u7uLFu2jDlz5mBtbU3jxo3p1asXgJ58hoeHs2rVKtatW8eBAwcYPHgwQUFBxMbGFiuGQ4cOMW/evL/9ZSYqKoqpU6eydOlSDh06xPTp0wkODiY6OrpQ223btnHx4kWCg4OLNX4Be3t7/QBo0qSJfh4VFWXUdubMmQwaNAhHx8K7lE2bNs2or0uXip8QFOe1uJeWLVsyZ84chgwZgpWVFV26dOGVV14B/vwz++qrr7hw4QK1atXC0dERFxcXXFxcCv1CERUVRXR0NMOHD2fIkCElmlcuTCDlSuXYYvpOBntuKnv96bYQQlRaxZmHYaqjtHOQ165dq5RS6plnnlGzZ8+WOch34eDgoKysrJTBYFAGg0GZmZkpMzMzZTAYVGJiolHbq1evqtu3b6sffvhBAerixYtKqby5wP/973+N2rZq1UpNmDChWDH85z//0cc0GAzKyspKAcrKykqNGjVKKaVUly5d1JgxY4yu6927t3rllVcK9de3b1/Vo0ePYr8GBc6cOaMfgIqMjNTP75xnvG/fPmVlZaUuX75cZD83btww6qtgrnZxFOe1KK6cnBwVFxenMjMz1aJFi5SVlZXKzMwsFGtSUpJKS0tTFhYWavny5Xftr2vXrmrIkCH3FYMwkXmPKfX1P1X90E2mjsTIkQ+bKbWit6nDEEI8oCjmHOSy/9h9JTR06FBGjhxZaE6myBMdHU1OTo5+Pn78eJKTk/niiy+oVauWUduCKRPffPMNPj4+eHl5AZCcnIyZmfEbEhYWFty+fduorGCKhr29Pc7Of+7wFRwcTI8ePfTzuLg4unfvztdff03Hjh3va4yEhAQ2bNjAypUr7+t1APD19TU6r1+/Pt7e3oXazZgxg379+uHpWfQTOicnpxLNHYbivRYFUlNTSUxMxNnZWX/qfSczMzPq1q0L5P2ZdejQAUtLy0KxQt4T5ZycnL+ds21ubs4ff/xRovsSFSzlCnh3vHe7CnZFudA0WZ4gCyEqtwciQe7Rowe5ubmcOnXK1KFUSg0bNjQ6d3BwICsri0aNGull4eHheHl54erqSkREBMuXLycsLEyv79GjB7NmzcLPz48GDRqwYcMGoqOjmT59ulHfV65coUGDBgQHB7Ns2TK93NHR0WiqQsEHxurXr0+dOnX0MRYtWkT79u1p3rw5u3fvZv369SxYsMBojOXLl1OjRg19GkhZO336NBs3buTYsWPl0n9xXosC4eHhhISEEBYWxoABA4zqli9fTosWLbCzs2PJkiXs3LmTH3/8Ua/fu3cvt27dolGjRvzyyy+8++67DBgwgPr16wPw73//m7p169K8eXPMzc2JiIhg27ZtRERElMt9izKUngIZt6Cmh6kjKSROOUPKcVAKyunzAUIIUVoPRIJsZmbG4MGDee+990wdSpUVGxvL0KFDSUpKomHDhixdupTXXntNr1+wYAHjx4/ntddeIykpCV9fX1auXEnXrl3LLIaPP/4YTdMYNWoUv/32G15eXsyaNYuQkBCjdkuWLKFfv36lfsdA3bGpxp1mzZpFz549jX6BqIwOHz7M6NGjSUtLw9/fnw0bNvDEE0/o9bdv3+btt9/m8uXLuLq6EhISwscff6zXW1lZMXXqVGJjY1FK0bBhQ5YvX05QUJApbkfcj5S4vK+VbQ4yeU+QyUqDtBtg53zvC4QQwgS0uyUBlUFAQIC6393DhChPcXFx+Pj4sHPnTtq2bWvqcIQo2pkfYVUfeH0r3vNvEDsj0NQR6d6c8BGLrf4Nb+4Aj1amDkcI8YDRNO2gUirgXu2q/SoWQpSlK1eu8K9//UuSY1G5FawSUQmfIOvrMt+KM20gQgjxNx6IKRZClJW2bdtKciwqv5Q40MzA3g04bOpojOgJcv5W2EIIURlJgiyEENVNyhWoURfMK98/8bYOrqSnW/LVpiimrffSyz1q2bB7XNl9ZkEIIUqj8v3rKYQQonRuxYFD5VvBAmD3+G4w14tB7lYM6vvn3GjvcZtNGJUQQhiTOchCCFHdpFyulPOPdQ6eMsVCCFGpSYIsysTixYt59NFHqVmzJi4uLvTt27fY20xD3jJymqYVedy5lfS6devw8/PD2tqagICAQttML1myhMcffxx7e/tCWyaXlwULFtCgQQOys7PLrM+srCyGDRuGk5MTNWvWJCQk5L436IiLiyMwMBBbW1vq1q3LzJkzC7X56KOPcHNzw87OjqCgIK5fv67X/fzzzzz11FO4uLhQo0YNOnbsyM8//1zqexPlLDcXbl2tlGsg6xw8/1yKTgghKiFJkEWZ2LVrF2+88QZ79uxh27Zt3Lx5kx49ehQ7aaxXrx7x8fFGx4cffkj9+vUJCMhbjeXYsWO8+OKLvPrqqxw6dIhmzZrxzDPPkJKSovfz+++/06NHD0aMGFEu9/lXOTk5zJ49mzFjxmBhUXYzliZOnMjatWsJDw9ny5YtREZGMmrUqPvqo0+fPty6dYs9e/bw6aefMnHiRFatWqXXL1y4kFmzZjF//nyioqKIi4ujX79+ev2+ffto164dmzZt4pdffqFNmzY8/fTT/Prrr2V2n6Ic/JEAOZngUM/UkdxdTQ9IvQY5WaaORAghilac/ahNdbRu3bpE+2zv2LFDubu7q549eyoHBwf13//+VzVt2lS5uLiobdu26W0AlZCQoF8XGBiogoODSzSmMHbo0CEFqJiYmBL30bhxYzVp0iT9fOTIkaply5b6eXp6urK3t1eLFy8udO3atWtV3o93yQHqwoULf9tm9erVysXFRaWlpZVqrDtlZ2crJycn9dlnn+llK1euVAaDQaWmpharj19++UUB6vDhw3rZG2+8oTp06KCfN2/eXI0ePVo/37VrlwLUmTNniuwzNzdXOTg4qDlz5tzvLYmKdOWAUpNqKnVyk1JKqfqhm0wcUBEOLMuLMemiXlQp4xRCVDvAAVWMHLTaPkG+du0aY8eO5fXXX2fEiBHMnTuX4OBgPv30U1OH9kBITk4GoFatWiW6fvfu3Zw8edJol7zo6Gg6deqknxsMBtq0aYMpN5OZOXMmw4cPx8bGpsz6PH/+PDdv3jS6186dO5ORkcHRo0eL1Ud0dDQODg40a9bMqI9Dhw6hlNL7unOMdu3aYWlpedfX8/bt22RmZhptgy0qjw4ztuM9bjOD520EIHB5LN7jNuNRq+x+NstMwQcIZZqFEKKSqrYJsouLC507d+Yf//gHLi4uPP7443Tr1u2+5sWKksnJyWHixIn06dOH+vXrl6iPxYsX061bN6PrExIScHZ2Zv369Tg6OnLgwAGcnZ1JSEgoq9Dvy9atWzlz5gzDhg0r034L7sfZ2Zm+ffvStm1bnJ2djeqK00ft2rVJS0ujQYMGvPfeezg7O3P79m1SU1O5ceMGubm5ODs7M3bsWHx8fEhPT8fR0fGuY0yZMgVXV1f69OlTNjcqylRc8m1iZwSyoGcdADZ/8E9iZwRWzqXTauZ/gFA+qCeEqKSqbYJc8ETPxsZG/97a2prbt2+bMqwHwvDhw7lx4waLFy8u0fUpKSmsXbuW119/vcj6GjVq4OXlha2tbWnCLJK9vb1+ADRp0kQ/jwfcD1cAACAASURBVIqKMmo7c+ZMBg0aVOQT1WnTphn1denSpRLF4+7ujpeX170b3oW5uTn16tXD1dX1rm1cXFzw8vL62znUq1evZt68eURERJTL6y7KUMoVsLABWydTR3J3BU+Qb0mCLISonB64dZDzpp9Q5AoHubm5FR1OtRMaGsqWLVuIiooq8fSKVatWYTAYeP75543KXVxcSExMpFu3bhw+nLc7WGJiIg0bNix13AViYmL07xs2bMh3332Hh0fef+YFXwH279/P7t27+eqrr4rsZ/Dgwbz44ov6ed26dYsdg4uLC5B3b3PnzgXg8uXLRnXF6ePGjRsYDAZ27twJwIoVK7CxscHe3h4rKyvMzMxITEwkNDSU0NBQcnNzSUpKKjTGxo0bGTRoEOvXr6dVq1bFvg9hIreu5CWgFbSKS4kYaoC1g0yxEEJUWtX2CfK9FCRvqampellBEiJKZvLkyXzzzTds377dKJm8U3Z2NrGxsSQmJt61n8WLF/Pyyy9jbW1tVB4QEGD0FDcjI4P9+/frq1yUBV9fX/0AqF+/vn5+5zzjGTNm0K9fPzw9i15r1snJyaiv+1nhwsfHB0dHR6N73blzJwaDgaZNmxq1TU1NJTY21ujnGPJeq5SUFI4cOWLUR6tWrdA0DYPBgL+/v9EYe/fuJSsry+j13Lp1K6+++ipr1qyha9dK+Fa9KCwlrnKvgfz/27vz+Ciru///75OQhASyQQKEJJCAEW9xQ1PQUuUneitILbigvVv94opyu1erVKnW1gW3Vm2r4lKUqgVFLFSs3iooLhVBZVEEQQiQQNgTiCRAMuf3xywksmWZ5Fwz83o+HvOYmWuuXPO5MhjfOflc5wSlMRcyAO+K2YBcVFSktLQ0vfLKK5Kkf/3rX/rmm28cVxW5xo8fr8cee0wvv/yy2rdvr/LycpWXl+/T0lJaWqrCwkLdcsst+z3O/PnztWDBgv22V1x55ZVavHix7rnnHi1ZskRjxoxRUlKSRo4cGdqnvLxcCxYsCPWaL1iwQAsWLNDu3bvDdq7Lli3TjBkzdOutt4btmPXFx8dr9OjRuvfeezVr1ix98sknGjdunC6++GJ16NChwb5Tp05VYWGhpk6d2mB7v379dOKJJ+qaa67RwoULNXnyZE2aNEljxowJ7TNmzBhNmDBB06ZN0xdffKEbb7xRp59+euiXg9mzZ+vcc8/Vww8/rH79+oU+0/rT6sGDKkv39vh6WXoeLRYAPCtmA3JKSoqefPJJPfbYY8rJydHMmTN1+umnuy4rYj311FOqqKjQySefrJycnNBtypQpTTrOM888o2OOOUYnnHDCPq8dddRRmjJliiZNmqR+/fpp0aJFevPNN5Went6gjn79+unXv/61JH9Q7Nevn9atW9fkc7LWqqCgYJ/tDz74oM4++2wdccQRTT5mY/3+97/XyJEjdd5552nIkCEaNGiQHn300SYd49VXX1VqaqpOPPFE3XTTTbr77rsbzHN89dVX69e//rXGjBmjn/zkJ8rNzdWLL74Yev2FF17Qzp07ddVVVzX4TG+44YawnSfCrHa3VLUhMkaQ03NpsQDgWSbYk+tFxcXF1uUUXsAPlZWVqVevXpozZ44GDBjguhwgpGDsTJXc1ld67FjpZ3+Rjr/YdUkHN+dhadYfpNvXS4kp/vrHD3NdFYAoZ4z53Fp7yN7MmB1BBpqjtLRUDz30EOEY3hQckU338DLTQcGV/rYzigzAe2JuFgugJQYMGEA4hncFL3rz8jLTQaHFQkqlrPDNRAMA4cAIMgBEi+BFb2mRMILMYiEAvIuADADRorJUSu4kJUbAYi6p3SUZWiwAeBIBWVJBQYEefvhh12UAQMtUlkVG/7EktUuUOnaRKpl/HoD3EJAlzZs3T//7v//ruoyI99RTT6lnz55KTk7WqaeeqhUrVjT6a2tra3XzzTerT58+SklJUY8ePXTLLbdo586dDfZ7/fXX1adPH7Vv317FxcWaN29e6LXnn39exph9bj+cOzjcnnrqKRUWFqq2tjZsx9yzZ4+uvfZaderUSWlpabr00kv1/fffN+kYZWVlGjZsmFJSUtS9e3c98MAD++xz9913q1u3burQoYNGjBihjRs3hl774IMPdOaZZyo7O1upqan6yU9+og8++KDF54ZWVFkaGf3HQel5TPUGwJMIyPIvy5uSEgF/kvSwt956S9dee63uuOMOffbZZ0pJSdHZZ5+turq6Rn39rl27tHjxYt13331auHChJk6cqFdeeUXXXXddaJ+vvvpKF1xwgS666CJ98cUXOuaYY3TWWWeFFq648MILtX79+ga30047Teedd16rnLMk1dXV6eGHH9bNN9/cpNXyDuXOO+/Uq6++qqlTp+qtt97S+++/rxtvvLFJxzj//PO1fft2ffLJJ3rkkUd055136qWXXgq9PmHCBD344IN64okn9OGHH6qsrKzBPMlz587VSSedpDfeeENffvml+vfvr6FDh+rbb78N23kizLaXRkb/cVBaLi0WALzJWuvZ2wknnGBb4s4777R5eXk2KSnJ9u7d2z7++OMNXu/Tp4+VZCXZhx56aJ+vX7lypT355JNtUlKSLS4utuPGjbMdOnSw1lq7atUqK8leccUVtmPHjva+++6zp5xyis3IyLCTJk2y1lpbW1trL7/8cltYWGgTExNtfn6+veuuu2xdXV2LzsuLhg8fbs8555zQ89LSUivJvvPOO80+5h//+EebkZERen7DDTfYfv36hZ7X1NTYjh072meeeWa/X19WVmbj4+Pt7Nmzm/X+kuyqVasOus/kyZNtdna23blzZ7PeY39qa2ttp06d7GOPPRba9uKLL9qkpCRbVVXVqGN8+eWXVpJduHBhaNsVV1xhBw4cGHp+7LHH2ptuuin0/KOPPrKS7PLly/d7TJ/PZ9PT0+2jjz7a1FNCGzjytletvSvN2g//5LqUxvv3WGvv6Watz2d73vaG62oAxABJ820jMmjUjiBPmzZNDz74oJ588kktXbpUTz/9tFJTUxvs89FHH2n9+vXKy9v/qlMXX3yxjDGaN2+exo0bp7/85S/77DN48GDdfffduv322zV69Gjdfffduv/++yX5Rxfj4+M1ceJELV26VE8++aQeffRRTZgwIfwn7Ni8efN08sknh57n5uaqV69easlCLxUVFcrMzDzgeyQlJal///4HfI+JEyeqoKBAgwYNanYNh/LAAw/ouuuuU3JyctiOuXLlSm3durXBuZ5yyimhUfbGmDdvntLT03XMMcc0OMYXX3wha23oWPXf46STTlJCQsIBv5/V1dXavXt3g88E3pFjtvgfRMIqekHpedKenVL1NteVAEADURuQV61apYyMDA0ZMkQFBQUaPHiwLrnkkgb7ZGVlqVu3boqPj9/n67/++mt9/PHH+vOf/6yjjz5aw4cP189//vN99hsxYoSGDh0qSTr33HN1xhlnqKSkRJKUmJioCRMmaNCgQSosLNSwYcM0fPhw/fvf/w77+bq2adMmZWVl6fHHH1eXLl20du1aZWVladOmTc063rp16/SXv/xFN9988z7v8c9//lOZmZmaP3/+Ad/DWqvnnntOl156qYwxzT6vg3n77be1fPlyXXvttWE9bvB8srKyNHLkSA0YMEBZWVkNXmvMMTp37qydO3eqsLBQt956q7KyslRdXa2qqipt2bJFPp9PWVlZuuWWW9SrVy/V1NQoMzPzgO/xhz/8QV26dNH5558fnhNFWOVGYkAOtoPQZgHAY6I2IJ9zzjmSpKKiIl155ZX6+9//rt27dzf665cvX664uDj17ds3tK3+46Dk5OTQ6GFycrLat2+v6urq0OtPPfWUiouLlZ2drY4dO+of//iHqqqqmntantepUyf17NlTSUlJzT5GVVWVRowYoSFDhuiaa67Z5/XU1FT16NHjoH3j7777rlavXq1Ro0Y16b07duwYukn+zzz4/MMPP2yw7wMPPKDRo0fvd0T1vvvua3CsNWvWNKmOoJycHPXo0aNZXytJ8fHxys/PV5cuXQ64T3Z2tnr06HHQHurJkyfrr3/9q6ZNm0a/vkeFRpAjqQeZuZABeFTUrqTXq1cvfffdd3rvvff0/vvv6/rrr9fLL7/cpqO3U6ZM0Q033KBHHnlEgwYNUnJysu644w5t2LChzWpoK9nZ2dq8ebNuuukmXXTRRZKkzZs3Kzs7u0nHqa6u1tlnn60uXbro+eef3+97nHbaaVq4cGHoPYqK9l2F65lnntEZZ5xxwPaZA1mwYEHocVFRkd58803l5voDR/Bekj777DN9/PHHmjRp0n6Pc/XVV+uCCy4IPe/evXujawh+zzZv3qzHH39ckrR27doGrzXmGFu2bFFSUpLmzJkjSfr73/+u5ORkdezYUYmJiYqLi9PmzZt122236bbbbpPP59O2bdv2eY8ZM2Zo9OjR+uc//6njjz++0eeBttXdbJZMnJSa47qUxmsQkBv/3wgAtLaoHUGWFJpJ4ZFHHtGECRP01ltvqaamplFfW1RUJJ/Pp6+++iq0rf7jxvjoo4900kkn6dprr9XRRx+tww47TCtXrmzSMSJFcXFxgxHWsrIyrVy5UsXFxQ32q62tVUlJiTZv3rzPMXbt2qVzzjlH7dq109SpU5WQkHDQ99i1a5c+++yzfd5j06ZNmj59ui677LImn8dhhx0WuklSz549Q8/r9xmPHz9ev/zlLw8YwDt16tTgWE2Z4aJXr17KzMxscK5z5sxRUlKSjj766Ab7VlVVqaSkZJ+/ShQXF6uyslKLFi1qcIzjjz9exhglJSXpqKOOavAe//nPf7Rnz54G38+3335bF110kaZMmaLBgwc3+hzQ9rqbrf5wHB9B4x4dukhxCbRYAPCcqA3IkyZN0rPPPqslS5Zo6dKlmjJlioqKitS+fXtJ/pHK8vJylZeXq66uTjt27FB5eXmo/7Jv374aOHCgrr/+ei1evFjTp0/Xa6+91qQaDj/8cH355Zd699139e2332rs2LFaunRp2M/VC8aMGaMZM2bomWee0eLFizV69GgdccQROvXUUxvsV1paqsLCQt1yyy0Ntu/Zs0fnn3++NmzYoKeffloVFRUNPh9JuvLKK7V48WLdc889WrJkicaMGaOkpCSNHDmywbFeeOEFpaamavjw4a1yrsuWLdOMGTN06623tsrx4+PjNXr0aN17772aNWuWPvnkE40bN04XX3zxPnM6T506VYWFhZo6dWqD7f369dOJJ56oa665RgsXLtTkyZM1adIkjRkzJrTPmDFjNGHCBE2bNk1ffPGFbrzxRp1++umhXw5mz56tc889Vw8//LD69esX+jyC0+rBW7prc2T1H0tSXJyUlkOLBQDvacxUF65uLZnmbfr06XbAgAE2NTXVpqWl2dNPP90uXrw49PrEiRNDU7zVv/Xs2TO0T3Cat8TERPujH/3I/va3v7WdOnWy1u6d5u1gj3ft2mUvv/xym5GRYTMyMux1111nr7rqKjto0KBmn5eXPfHEEzY/P98mJSXZQYMG2W+//XaffYLfn1GjRu13+/5u9adae+2112xRUZFNTEy0J5xwgp07d+4+79GnTx97/fXXh/v0Qi677DI7YsSIVju+tf5/O9dcc43NyMiwqampdtSoUfud4i3473jixIn7vLZ27Vo7dOhQ2759e9utWzd7//3377PPXXfdZbt06WKTk5Pt8OHDbXl5eei1UaNG7ffz+OFnB29Y+dsia1+5xHUZTfe3odY+N4Rp3gC0CTVymjfj39ebiouLbUumCQu3u+66S2+88YY+//xz16XAkbKyMvXq1Utz5szRgAEDXJcD+FmrXb/LVtLAMdIZf3BdTdO8dqW05lMVbBivkvHDXFcDIMoZYz631hYfar8IalZre9OnT9eePXvUr18/fffdd5owYYLGjh3ruiw4VFpaqoceeohwDG/5frOSzJ7Ia7GQpPRcacc6xcnnuhIACCEgH0RNTY3uuOMOrV27Vjk5ObrqqqvCPuctIsuAAQMIx/CeSv8sJ5EZkPMkX62yRG87AO8gIB/EhRdeqAsvvNB1GQBwcMFZICJpDuSgNH+ozzX7zmwDAK5E7SwWABAzgrNApOe7raM50v2hPrTQCQB4AAEZACJdZalqbIKU0sl1JU0XaAvpTkAG4CEtCsjGmAxjzFRjzFJjzDfGmJOMMZ2MMe8YY5YH7jMD+xpjzOPGmBXGmEXGGJbkAoBwqCxVmc2SjHFdSdO1z5ASOhCQAXhKS0eQH5P0lrX2CEnHSvpG0lhJ71lriyS9F3guSUMlFQVuoyU92cL3BgBI0vYyrbcROHos+UN9ei4tFgA8pdkB2RiTLukUSc9JkrV2t7W2QtJwSS8EdntB0ojA4+GSJgXmaf5UUoYxJqfZlQMA/CpLtc5mua6i+dIIyAC8pSUjyIWSNkmaaIz50hjzrDGmg6Su1tr1gX3KJXUNPM6VtLbe15cGtjVgjBltjJlvjJkfXPYZAHAAdXukHeVar86uK2m+9Dx1N1tdVwEAIS0JyO0kHS/pSWttP0nfa287hST/esvyL0/baNbap621xdba4uzs7BaUBwAxYPs6SVZlNrIDchdTIdXucl0JAEhqWUAulVRqrZ0beD5V/sC8Idg6EbjfGHi9TFL9OYjyAtsAAM0VmAN5fSQH5OD8zdvXua0DAAKaHZCtteWS1hpj+gQ2nSZpiaQZkkYFto2SND3weIak/xeYzeJESZX1WjEAAM0RmAN5XSQH5PRgQGbMBIA3tHQlveskvWSMSZS0UtKl8ofuV4wxl0taLemCwL5vSjpL0gpJOwP7AgBaIioCcuCPi5UEZADe0KKAbK1dIKl4Py+dtp99raRrWvJ+AIAfqCyVkjNVXdPedSXNF2yxqFx78P0AoI2wkh4ARLLtZVJanusqWiYxRdtsR1osAHgGARkAIlllaWi55ki23namxQKAZxCQASCSVZbuvcgtgpXZzowgA/AMAjIARKpdVVJNRRSNIJe6LgMAJBGQASByBUdcI70HWYGAXFPhD/0A4BgBGQAiVXDWhygYQV5nO/kf0GYBwAMIyAAQqYIXtUVBD3JoJUDaLAB4AAEZACJVZalk4qTUHNeVtNg6ZfkfMIIMwAMIyAAQqbaXSR27SfEJritpsXKbKckwggzAEwjIABCpKtdGRf+xJNWqndSxK3MhA/AEAjIARKrKsqjoPw5Jz5W2M4IMwD0CMgBEImv9LRZRMoIsyX8ujCAD8AACMgBEop1bpNqaqJgDOSQtz9+DbK3rSgDEOAIyAESi4MVsUTWCnCvVVkvV21xXAiDGEZABIBKFAnIU9SCnBc6FmSwAOEZABoBIFJwvOD3fbR3hFDwX5kIG4BgBGQAiUeVaqV17KaWz60rCJ50RZADeQEAGgEhUWeZvSTDGdSXh06GLFJdAQAbgHAEZACJRZWl09R9LUlyclJZDiwUA5wjIABCJtpdFV/9xUBpzIQNwj4AMAJGmrlbasX7vrA/RJD2P1fQAOEdABoBIs2O9ZH3RNQdyUHqutH2d5KtzXQmAGEZABoBIE41zIAel5Uq+Wqlqo+tKAMQwAjIARJponAM5KDgqzoV6ABwiIANApKlc67+P1h5kianeADhFQAaASFNZJrXPkJI6uq4k/FhuGoAHEJABINJUlkbnBXqSlJwpJaTQYgHAKQIyAESa7VEckI3xjyIzggzAoXauCwAANN7A8bP0RvUqvVHaVb8dOzO0PTcj2WFVYZaexwgyAKcIyAAQQbZWbFNm+ypdfOZAXXzyMNflhE1uRrIKAoH/gXZWp8Z/p/5jZyo3I1kfjx3suDoAsYaADAARpLvZ4n8QZS0WDULw7EXSBx+o5J7/VsG4d9wVBSBm0YMMABEkx2z1P4iygNxAeq4kK+1Y57oSADGKgAwAEaS72ex/EI1zIAeF5kKmDxmAGwRkAIgg/hYLI6V1d11K60ljsRAAbhGQASCCdNcWKbWbFJ/gupTWkx4YHd9OQAbgBgEZACJId7NZSs93XUbrSuzgXymQFgsAjhCQASCCdDdbovsCvSDmQgbgEAEZACKFz6fcWArI9CADcISADACRYudmJZk90d9iIbHcNACnCMgAECkq1/rvM2IgIKfnSjUVSlaN60oAxCACMgBEiuCIaiy0WASmegutHAgAbYiADACRoiIwghwLATlwjqGVAwGgDRGQASBSVJaqyrb3T4EW7QJzIYdWDgSANkRABoBIUblW62xnyRjXlbS+1O6SDC0WAJwgIANApKgsVZnNcl1F22iXKHXsohzRYgGg7RGQASBSVK7VulgJyJKUnqccRpABOEBABoBIsHuntHOLymxn15W0nbRc5dKDDMABAjIARILAssuxN4K8VbLWdSUAYgwBGQAiQWCRkHUxNoKcYnZJ1dtcVwIgxhCQASASBOZAjpmL9KTQVG/B0XMAaCvtXBcAAGiEylLJxGmDMl1X0nbS/UtqX/7463rPt6bBS7kZyfp47GAXVQGIAQRkAIgElaVSao5qq2Pox3aafwT5uRE5Uv9hDV4qGDvTRUUAYgQtFgAQCSrXxsYS0/V17CLFtaPFAkCbIyADQCSIxYAcFy+ldfePngNAGyIgA4DX+XxSZVmoJzempOeHLlAEgLZCQAYAr/t+o+TbE3sjyJKU0SM0xR0AtBUCMgB4XbDFIFZHkHesl2p3u64EQAwhIAOA11UEpjiL1RFk6+NCPQBtioAMAF4XGkGOxYAcGDWnzQJAGyIgA4DXVZZKSWlScobrStpesK2EC/UAtCECMgB4XWVpbI4eS4HzNnvbTACgDRCQAcDrKtfEbkBulySldqPFAkCbIiADgNfF8giy5L9QjxFkAG2IgAwAXrarSqreFtsBOT2fgAygTRGQAcDLYnkO5KCMfP80b74615UAiBEEZADwsmDvbUYPt3W4lNFD8tVKO8pdVwIgRhCQAcDLtpX47zN6Oi3DqfTALwe0WQBoIy0OyMaYeGPMl8aYNwLPC40xc40xK4wxU4wxiYHtSYHnKwKvF7T0vQEg6lWskeITpY5dXVfiTnD0nJksALSRcIwg3yDpm3rPH5D0J2vtYZK2Sbo8sP1ySdsC2/8U2A8AcDAVa/z9x3Ex/Ae/4AWKjCADaCMt+olrjMmTNEzSs4HnRtJgSVMDu7wgaUTg8fDAcwVePy2wPwDgQCpWS5kx3F4hSYkpUkoWARlAm2npkMSjkm6V5As87yypwlpbG3heKik38DhX0lpJCrxeGdi/AWPMaGPMfGPM/E2bNrWwPACIcBVrYvsCvaCMHrRYAGgzzQ7IxpifStporf08jPXIWvu0tbbYWlucnZ0dzkMDQGTZVSXt3BLbF+gFZeRLFQRkAG2jJSPIAyX9zBhTImmy/K0Vj0nKMMa0C+yTJ6ks8LhMUr4kBV5Pl7SlBe8PANEt2FLACPLeEWRrXVcCIAY0OyBba39jrc2z1hZI+rmkWdbaX0qaLen8wG6jJE0PPJ4ReK7A67Os5ScdABxQKCAzgqz0HlJtjfQ9rXcAWl9rXBZ9m6RfGWNWyN9j/Fxg+3OSOge2/0rS2FZ4bwCIHsGAHOsX6Un+FguJNgsAbaLdoXc5NGvt+5LeDzxeKan/fvapkTQyHO8HADGhYrXULlnqwPUYoTaTitVS3gluawEQ9WJ4Yk0A8LiK1f5gyIyY/rmgJWayANAmCMgA4FVM8bZX+zSpfQYtFgDaBAEZALxq22oCcn0ZPfyj6gDQygjIAOBFNZVSTQUX6NWXWeD/pQEAWhkBGQC8KNhKwAjyXpk9/SPIPt+h9wWAFgjLLBYAgPC67bl/6QFJZ79YqsV2Zmh7bkayu6Jcyyzwz4VctcF1JQCiHAEZADwoZWeZlCD9a9wvpA6dXZfjDZkF/vttJS6rABADaLEAAA/KM5ulhA5SSifXpXhHZqH/noAMoJURkAHAg/LNRn/PLXMg75WeL8kQkAG0OgIyAHhQntnMBXo/1C5RSs9jqjcArY6ADAAelGc2SRlM8baPzAJGkAG0OgIyAHhN9TalmZ2MIO9PRk8CMoBWR0AGAK+pWOO/JyDvK7NA2rFeSdrtuhIAUYyADABes3WV/75Tods6vCgw1Vue2eS2DgBRjYAMAF6zLRCQg/P+Yq/A96SH2ei2DgBRjYAMAF6zdZU22zQpKdV1Jd5DQAbQBgjIAOA121Zpje3iugpv6pAlJXQgIANoVQRkAPCarSVabbu6rsKbjJEyC/wLqQBAKyEgA4CX1O6WtpdqDQH5wDJ7EpABtCoCMgB4ScUayfq02keLxQFlFvhbLKx1XQmAKEVABgAvCcxgQYvFQWQWqIPZJX2/2XUlAKIUARkAvCQwBzItFgcRnP6OFfUAtBICMgB4ybZVUkIHbVK660q8i4AMoJURkAHAS7auCgRA47oS7wouwU1ABtBKCMgA4CXbVrHE9KEkJKvcZu5dcRAAwoyADABe4fP5R0VZYvqQ1tguoX5tAAg3AjIAeEVVuVRbwwhyI5T4uklbv3NdBoAoRUAGAK8IjohmEpAPpcR2k6o2SLt2uC4FQBQiIAOAVwR7ahlBPqRVtpv/wdaVbgsBEJUIyADgFVtXSSZeSs93XYnnlQQD8hbaLACEHwEZALxi60opI1+KT3BdieeVBBdSoQ8ZQCsgIAOAV2xZIXU+zHUVEaFa7aXUHGayANAqCMgA4AXW+keQCciN16k3LRYAWgUBGQC8oGqDtLuKgNwUnXvRYgGgVRCQAcALtqzw33fu7baOSNKpt/T9Jqlmu+tKAEQZAjIAeEEoIDOC3GjBXyYYRQYQZgRkAPCCLSuk+CQpLc91JZGjUyAg04cMIMwIyADgBVu+84+IxvFjudGCC6qwWAiAMGvnugAAiGUDx89SWUW13k1coOU2V2PGzpQk5WYkO64sAiQk+0fcGUEGEGYEZABwlebOBQAAE+pJREFUqKyiWiX3nindO0qH/fhClZw+zHVJkYWZLAC0Av6WBwCuVa6RfHu4QK85mAsZQCsgIAOAa8GAR0Buus69peqtUvU215UAiCIEZABwLTjFWyfmQG6yzkX++83L3dYBIKoQkAHAtS0rpKR0qUOW60oiT/bh/vtNy9zWASCqEJABwLUtK/ytAsa4riTyZPT0zx+9+VvXlQCIIgRkAHBty3f0HzdXXLz/e0dABhBGBGQAcKi9dkmVawnILZFVREAGEFYEZABwqLdZ738Q7KVF02X3kbaVSHtqXFcCIEoQkAHAocNMqf9B9hFuC4lkWYdL1seCIQDChpX0AMChorgyycQzxVsT5WYkqyCwLPeRZoPeTJL+97HJWph2qj4eO9hxdQAiHQEZABwqMmX+GSzaJbouJaI0CMG7d0r33aEnzuiogreq3RUFIGrQYgEADh1myvw9tGi+xBQpI58L9QCEDQEZAFyp3aWeZgP9x+GQ1UfazGIhAMKDgAwArmz5Tu2Mj4AcDlmHS5tXyMjnuhIAUYCADACubFrqv89iircWyz5cqq1WrtniuhIAUYCADACubFomnzX+hS7QMoFfMg4zZY4LARANCMgA4MrmZVpju0gJya4riXxZ/gsdexOQAYQBARkAXNm0TMttrusqokOHzlJKlg4nIAMIAwIyALhQVyttXq4VBOTw6Xqk+sStcV0FgChAQAYAF7atknx7tNxHQA6bLn39I8g+ZrIA0DIEZABwITCDxXKb57iQKNL1SKWYXf5fPgCgBQjIAODCRn9A/s52d1xIFOnS13+/cYnbOgBEPAIyALiwYbGUWaidau+6kujR5Qj/tHkbCMgAWoaADAAubPha6naU6yqiS2IH/7R5G792XQmACEdABoC2tvt7act3UlcCcrgts/mMIANoMQIyALS1jd9IsgTkVrDU5ktbv5P2VLsuBUAEIyADQFvb8JX/vmtft3VEoWW+fMn6pE3LXJcCIII1OyAbY/KNMbONMUuMMV8bY24IbO9kjHnHGLM8cJ8Z2G6MMY8bY1YYYxYZY44P10kAQEQp/0pKTJUyerquJOoss/n+B8xkAaAFWjKCXCvpZmvtkZJOlHSNMeZISWMlvWetLZL0XuC5JA2VVBS4jZb0ZAveGwAi14av/aPHcfwRL9xKbDcpPsn/PQaAZmr2T2dr7Xpr7ReBxzskfSMpV9JwSS8EdntB0ojA4+GSJlm/TyVlGGNyml05AEQia/cGZIRdneKl7D4EZAAtEpbhC2NMgaR+kuZK6mqtXR94qVxS18DjXElr631ZaWDbD4812hgz3xgzf9OmTeEoDwC8o2KNtKuSKd5aU9e+BGQALdLigGyM6SjpNUk3Wmu313/NWmsl2aYcz1r7tLW22FpbnJ2d3dLyAMBbgsGNGSxaT7djpO83SjvKXVcCIEK1KCAbYxLkD8cvWWunBTZvCLZOBO43BraXScqv9+V5gW0AEDs2fCXJSF2OdF1J9Op+nP9+/UK3dQCIWO2a+4XGGCPpOUnfWGv/WO+lGZJGSRofuJ9eb/u1xpjJkgZIqqzXigEAsaF8sdSpUErq6LqSqJSbkay+T67T4iSjR194RY/X1Ya2fzx2sOPqAESKZgdkSQMlXSxpsTFmQWDb7fIH41eMMZdLWi3pgsBrb0o6S9IKSTslXdqC9waAyLThKy7Qa0WhEPzn8fpV1k796n+GSZIKxs50WBWASNPsgGyt/UiSOcDLp+1nfyvpmua+HwBEvOoKaetK6bhfuq4k+nU/Tlr9iesqAEQoJuEEgLYS7Int3s9tHbEg5zhpe5lUxWxIAJqOgAwAbWXdl/57AnLryznWf8+FegCagYAMAG1l3Zf+5aVTOrmuJPrlHOO/X/+l2zoARCQCMgC0lfULGD1uK+3TpU69pXULDr0vAPwAARkA2sLOrdK2kr1z9KL1dT9OWr/IdRUAIlBLpnkDADTSjX96Xo9K+sWbu/XJG3unHMvNSHZXVLTLOVb66jX/LycA0AQEZABoA913LpUSpJd/O1pKznRdTmzICYzWr6MPGUDT0GIBAG3gmLiVUqdehOO21P04SUYq+9x1JQAiDAEZAFqbtTo+brmU9yPXlcSW9ulS9hHS2s9cVwIgwhCQAaC1VaxRF1NBQHYh/0dS6TwZ+VxXAiCCEJABoLWVzvPf5/d3W0csyusv1VSo0JS7rgRABCEgA0BrW/uZvrdJUpe+riuJPYFfSo6PW+64EACRhIAMAK1t7Vwt8vWW4pk4qM11LpLap+t4Q0AG0HgEZABoTbt3Shu+0ue2yHUlsSkuTsotVj9GkAE0AQEZAFrTui8lX62+8BGQncnvrz6mVKrZ7roSABGCgAwArWntXEnSl77DHBcSw/L7K85YpnsD0GgEZABoTas/lrKP0Dalua4kduUP0B4b7/8sAKARCMgA0FrqaqU1n0o9B7quJLYldtBiW0hABtBoBGQAaC3rF0q7q6SCn7iuJObN9f2XVPaF/6JJADgEAjIAtJaSD/33jCA7N9f3X5Jvj1RKHzKAQyMgA0BrWf2xfx7e1K6uK4l5832HSyZOKqHNAsChEZABoDUE+49pr/CEKqVIOcfShwygUQjIANAayhdJu7YTkL2k50CpdB59yAAOiYAMAK1h5Wz/fcHJbuvAXr1Plep2S6s/cV0JAI8jIANAa1gxS+p2NP3HXtJzoBSfJH33nutKAHgcARkAwm3XDmntp1Lv01xXgvoSkqWeP5a+m+W6EgAeR0AGgHBb9aHkq5UOIyB7zmGnSZuWSpVlrisB4GEEZAAIt+/ekxI6SPknuq4EP9R7sP8+2CMOAPvRznUBABBNBo6fpZd2/ksr7OG6Ytw7oe25GckOq0JIlyOljt2kFe9J/S5yXQ0AjyIgA0AYJVSuVEHSBhUMvUUlA4a5Lgc/ZIxUdLq05F9S3R4pPsF1RQA8iBYLAAijM+Lm+x/0Geq2EBxYn7OkXZVSyUeuKwHgUYwgA0AY/Xf851K3Y6SMfNeloJ7cjGQVjJ0pSWqvXfoyKVFTJj6hZzpafTx2sOPqAHgNARkAwqVqo04wy6UjfuO6EvzAPiH4H1N1yfqF+t1GVtUDsC9aLAAgXJb9W3HGSkfQe+x5R5wlbS9VX1PiuhIAHkRABoBwWTpTa33ZUte+rivBoRw+RDJxOiN+vutKAHgQARkAwqF6m7Rytt72FftnSoC3dciSeg7UT+M+lax1XQ0AjyEgA0A4LJkh1e3W9LqBritBYx19vnrHrZfWL3RdCQCPISADQDgsflXq1FuLbaHrStBY//Uz7bbx/s8OAOohIANAS21f559T95gLJNFeETFSOukD37HSV9Mkn891NQA8hIAMAC311TRJVjp6pOtK0EQz6n4s7Vgnrf7YdSkAPISADAAtYa204GWp+/FS596uq0ETveM7QUpKkxa85LoUAB7CQiEA0BKl86SNX0tnP+a6EjRD54wMTaoaoAsWTFX/uadquzpK8q+8xwp7QOwiIANAMwwcP0tlFdV6OOEpnRmXrAGvpmrnqzOVm5HsujQ0wcdjB0vrs6QJ72jRORXSgAslKbQsNYDYREAGgGYoq6hWyV0/lh65TDruF1ry0/Ncl4TmyjnG3yLz+fNS/9HMYw2AHmQAaLaFk6XaGumES1xXgpY64RJp4xJp9SeuKwHgAQRkAGiGeNVJ/3lCyh8g5Rzruhy01NEjpZTO0iePu64EgAcQkAGgGYbFzZUq10gDb3RdCsIhMUX60ZXSt29Jm5a5rgaAYwRkAGgqa3VVu39JWX2kw4e4rgbh0v9KqV176ZM/u64EgGMEZABoqhXvqm/camng9VIcP0ajRocs6bhfSoumKFebXFcDwCF+sgNAU/h80nu/11pftnT0Ba6rQbid/CtJRje0m+a6EgAOEZABoCmWvC6VL9IjtSOldomuq0G4pedJP7pC58XPoRcZiGEEZABorNrd0qx7pK5Habrvx66rQWs5+VeqVpL03u9dVwLAEQIyADTWp3+Vtq6UTv+dLD8+o1eHLD1Z+zNp6RvSinddVwPAAX7CA0BjbFstvf+AdMRPpaL/dl0NWtkzdcOkzodJb/5a2lPjuhwAbYylpgHgIAaOn6Wyip36W8JD6h9ndfqCM1W+YKZyM5Jdl4ZWlJ2Rpl+uH6mXEu/XY3dfrT/VjpQk5WYk6+Oxgx1XB6C1EZAB4CDKKqpVct4GaeYCach4fXri/3NdEtqAPwQPll4v0Q2LpuiG0WOkHgNUMHam69IAtAFaLADgIIpMqfT27VLv06T+V7kuB21t6INSer407QqpptJ1NQDaCAEZAA6kukJPJDwmJaVKI55kUZBY1D5NOu9ZqbJMmjZacfK5rghAG+CnPQDsT90e6dVLVGDKpfMnSqldXVcEV/L7S0MfkL59S2Pb/cN1NQDaAAEZAH7IVyfNuE5aOVu3114uFZ7suiK41v9K6UdXanS7mdInf3ZdDYBWRkAGgPp8ddL0a6SF/5BOHadX6/4/1xXBK4aM1xt1J0r/N076z19dVwOgFTGLBQAE7dohTRstLXtTOvUOadCvpX8zawEC4tvpwZSbZXY+pGFv366nZn6iB2p/Lqs4pn8DogwBGQAkafNy6dVLpI3fSEMfkgaMdl0RPGjOb86Q6gZLb43V1fOe0dVH1krDn1DBHz51XRqAMCIgA4htvjpp/t9U8+btqraJumHPLZrzeq70un/kmAVBsI/4dtKwh6XsPv4pAJ/8sU6Nu1jSMNeVAQgTY611XcMBFRcX2/nz57suA0CU8K+KVx14ZnVi3Dca1+5FHRVXorlxx2nAjZOltBynNSLClC+WXrtS2vSNdPgQ6fS7pS5HuK4KwAEYYz631hYfaj8u0gMQM8oqqlVy75kqGeVTyWGPanLiPToqY4903nMa8Nv3Ccdoum5HS1fN0b17fiGVfCQ9MUD6xy+k1f+RPDwABeDgaLEAEP321Ehr5+p37Z6XHrle2rlZyugpnfWw1O8iKYE2CrRAu0S9mTpSUytO0SXt/k+jlr6tjGUzVeLrqlkJg3TZZWOknONYaAaIIG3eYmGMGSLpMUnxkp611o4/0L60WABoMl+dtK1E2vCVVP6VtHau/1Zbo102QUl9h0nHXCAVnenvJQXCbVeVtGS6tGiKfCvnKM5YKbmT1PPHUvfjpJx+UtcjpY7dCM1AG2tsi0WbBmRjTLykbyX9t6RSSfMk/Y+1dsn+9icgo1EO9W+4Uf/GW3qMRrxHS48RLefRmGNYn1S3W6qtkWp3NbzfUyPtqpSqtwVuFdLOLVJlqVSxRtpeJvlqJUl11uhbm6f/+PrqY19frUk9Xu/85qeHrg8Ik5/eP029d8zTKfGLdLxZrsK4DXtfjE+SMnv6/5rRsYuU0tl/65DlX948IcX/142EZCmhg/8+PkGKayeZeCkueKv33MRJxrg7YcDjvBqQT5L0O2vtmYHnv5Eka+39+9vfSUBe/o70yqhD7BQBISSawhRwEN/bJFWoo9bbziqzWSqzWdrePk9jLz1fyv4vKTHFdYnAXtUVUvkijXt2mu4ZlOr/a0fFaun7zf5b3a6Wv4eJ/0FI/kFgPtBr+wTrA712sOMBjVQ4SPqfl9v8bb0akM+XNMRae0Xg+cWSBlhrr623z2hJwQlI+0ha1mYF7pUlabOD90Xb4TOObny+0Y/POLrx+UY/V59xT2tt9qF28lwDnrX2aUlPu6zBGDO/Mb9dIHLxGUc3Pt/ox2cc3fh8o5/XP+O2vjqgTFJ+ved5gW0AAACAJ7R1QJ4nqcgYU2iMSZT0c0kz2rgGAAAA4IDatMXCWltrjLlW0tvyT/P2N2vt121ZQyM5bfFAm+Azjm58vtGPzzi68flGP09/xp5eahoAAABoa8xQDgAAANRDQAYAAADqISD/gDFmiDFmmTFmhTFmrOt6EF7GmL8ZYzYaY75yXQvCzxiTb4yZbYxZYoz52hhzg+uaEF7GmPbGmM+MMQsDn/HdrmtC+Blj4o0xXxpj3nBdC8LPGFNijFlsjFlgjPHkksn0INfT1KWwEXmMMadIqpI0yVp7lOt6EF7GmBxJOdbaL4wxqZI+lzSC/4ajhzHGSOpgra0yxiRI+kjSDdbaTx2XhjAyxvxKUrGkNGst68NHGWNMiaRia61nF4NhBLmh/pJWWGtXWmt3S5osabjjmhBG1to5kra6rgOtw1q73lr7ReDxDknfSMp1WxXCyfpVBZ4mBG6M9EQRY0yepGGSnnVdC2IXAbmhXElr6z0vFf9zBSKSMaZAUj9Jc91WgnAL/Pl9gaSNkt6x1vIZR5dHJd0qyee6ELQaK+n/jDGfG2NGuy5mfwjIAKKOMaajpNck3Wit3e66HoSXtbbOWnuc/Kux9jfG0C4VJYwxP5W00Vr7ueta0Kp+Yq09XtJQSdcE2h89hYDcEEthAxEu0Jf6mqSXrLXTXNeD1mOtrZA0W9IQ17UgbAZK+lmgR3WypMHGmBfdloRws9aWBe43Snpd/hZXTyEgN8RS2EAEC1zA9Zykb6y1f3RdD8LPGJNtjMkIPE6W/6LqpW6rQrhYa39jrc2z1hbI///gWdbaixyXhTAyxnQIXEQtY0wHSWdI8tzMUgTkeqy1tZKCS2F/I+kVjy6FjWYyxvxD0n8k9THGlBpjLnddE8JqoKSL5R91WhC4neW6KIRVjqTZxphF8g9qvGOtZSowIHJ0lfSRMWahpM8kzbTWvuW4pn0wzRsAAABQDyPIAAAAQD0EZAAAAKAeAjIAAABQDwEZAAAAqIeADAAAANRDQAYAAADqISADAAAA9fz/xSrbpF7iRXgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "d = {'Entries': len(y_all),\n",
    "     'Mean': y_all.mean(),\n",
    "     'Std': y_all.std(ddof=1),\n",
    "     'Chi2': chi2,\n",
    "     'NDF': N_dof,\n",
    "     'Prob': chi2_prob, \n",
    "     'N': [minuit.values['N'], minuit.errors['N']], \n",
    "     'mu': [minuit.values['mu'], minuit.errors['mu']], \n",
    "     'sigma': [minuit.values['sigma'], minuit.errors['sigma']],\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.97, text, ax2, fontsize=14);\n",
    "\n",
    "fig2.tight_layout()\n",
    "fig2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "if save_plots:\n",
    "    fig2.savefig(\"Dist_ErrorProp.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "# Questions:\n",
    "\n",
    "0. First solve the problem of obtaining the Area, Circumference & Diagonal with uncertainty ANALYTICALLY.\n",
    "\n",
    "1. Now look at the program, and assure yourself that you understand what is going on. Put in the correct expression for y in terms of x1=L and x2=W in order to calculate the circumference, area, and diagonal length, and run the program. Does the output correspond well to the results you expected from calculations to begin with?\n",
    "\n",
    "2. Imagine that you wanted to know the central value and uncertainty of y1 and y2, given the\n",
    "   same above PDFs for `x1`=$L$ and `x2`=$W$:\n",
    "   \n",
    "     `y1 = log(square(x1*tan(x2))+sqrt((x1-x2)/(cos(x2)+1.0+x1)))`\n",
    "     \n",
    "     `y2 = 1.1+sin(20*x1)`\n",
    "\n",
    "   Get the central value of y, and see if you can quickly differentiate this with\n",
    "   respect to `x1` and `x2`, and thus predict what uncertainty to expect for y using\n",
    "   the error propagation formula. It is (for once) OK to give up on the first expression :-)\n",
    "   Next, try to estimate the central value and uncertainty using random numbers\n",
    "   like above - do you trust this result more? And are the distributions Gaussian?\n",
    "\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    "3. Try to generate `x1` and `x2` with e.g. a quadratic correlation, and see that despite\n",
    "   not having any linear correlation, the result on circumference, area, and diagonal\n",
    "   length is still affected.\n",
    "\n",
    "---------------------------------------------------------------------------------- "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}