{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Central Limit Theorem (CLT)\n",
    "\n",
    "### Authors: \n",
    "- Christian Michelsen (Niels Bohr Institute)\n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "- 15-11-2018 (latest update)\n",
    "\n",
    "***\n",
    "\n",
    "The aim of this notebook is to illustrate the Central Limit Theorem (CLT) through concrete examples.\n",
    "\n",
    "When you add random numbers from different distributions, but with the same variance (or standard deviation), together and plot the distribution of these sums, you end up with a Gaussian distribution, as dictated by the CLT. \n",
    "The example also illustrates how widths (and therefore uncertainties) are added in quadrature, as one has to divide the sum by the square root of the number of random numbers that went into the sum in order to get a Gaussian of unit width (when using random numbers of unit width, i.e. RMS $= \\sigma = 1$).\n",
    "\n",
    "For more information on the Central Limit Theorem, see:\n",
    "- **R. Barlow**: page 49 (and page 45 for Uniform distribution)\n",
    "- **G. Cowan**: page 33\n",
    "- __[Wikipedia: \"Central limit theorem\"](http://en.wikipedia.org/wiki/Central_limit_theorem)__\n",
    "***\n",
    "\n",
    "First, we import the modules we want to use:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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": "markdown",
   "metadata": {},
   "source": [
    "Here we set the random seed for the random number generator (RNG). This ensures reproducability (the same results every time the notebook is restarted). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.random                         # Random generator\n",
    "r.seed(42)                            # Set a random seed (but a fixed one - more on that later.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here we set the parameters for the experiement. We are going to play around with these more later on. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_experiments = 1000                  # Number of sums produced\n",
    "N_uniform     = 10                    # Number of uniform numbers used in sum\n",
    "N_exponential = 0                     # Number of exponential numbers used in sum\n",
    "N_cauchy      = 0                     # Number of cauchy numbers used in sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And finally we set up some global constants, like $\\pi$, and the bool flags about the program:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "pi = 3.14159265358979323846264338328  # Selfexplanatory!!!\n",
    "pi = np.pi                            # Another way of doing it - surely better!\n",
    "\n",
    "verbose = True                        # Print some numbers or not?\n",
    "N_verbose = 10                        # If so, how many?\n",
    "save_plots = False                    # Save the plots produced to file(s)?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Initial functions\n",
    "We load the external functions from the `AppStat2018/External_Functions` directory. The first line adds the relevant path to the places it will go look for functions, and the second line imports the function wanted. If you get an error here, check that you have copied \"External Functions\" to your directory, and that you refer to the correct path."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "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": [
    "## Loop over process:\n",
    "\n",
    "Having loaded everything that we need, we can start the actual program. We start out by initializing a counter to count how many of the produced sums that fall outside some range ($\\pm 3\\sigma$) and some zero-filled numpy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "N3_sigma = 0        # Counter for the number of produced sums, that fall outside +-3 sigma\n",
    "\n",
    "x_uniform = np.zeros((N_uniform, N_experiments))\n",
    "x_exponential = np.zeros((N_exponential, N_experiments))\n",
    "x_cauchy = np.zeros((N_cauchy, N_experiments))\n",
    "x_sum = np.zeros((N_experiments))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Note on code:\n",
    "\n",
    "The loops are here written out explicitly, in order to assure that you understand what is going on in them. Commented out below is the much shorter and faster Python version of the code, which you should make yourself familiar with immediately (yes, right now!). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   Uniform:     -0.4346\n",
      "   Uniform:      1.5613\n",
      "   Uniform:      0.8037\n",
      "   Uniform:      0.3418\n",
      "   Uniform:     -1.1916\n",
      "   Uniform:     -1.1917\n",
      "   Uniform:     -1.5308\n",
      "   Uniform:      1.2685\n",
      "   Uniform:      0.3503\n",
      "   Uniform:      0.7208\n",
      "At iexp :  500\n",
      "At iexp :  1000\n"
     ]
    }
   ],
   "source": [
    "for iexp in range(N_experiments): \n",
    "\n",
    "    if ((iexp+1) % 500 == 0): \n",
    "        print(\"At iexp : \", iexp+1)                       # Show progress!\n",
    "    sum_value = 0.0                                       # sum_value is the number we are going to add random numbers to!\n",
    "                                                          # According to the CLT, it should be Gaussianly distributed.\n",
    " \n",
    "    # Generating uniform numbers (with mean 0, and RMS of 1):\n",
    "    for i in range(N_uniform): \n",
    "        x = np.sqrt(12.0) * (r.uniform() - 0.5)           # Uniform between +-sqrt(3). Why? Possibly see references above!\n",
    "        sum_value += x\n",
    "        x_uniform[i, iexp] = x\n",
    "        if (verbose and iexp == 0 and i < N_verbose):\n",
    "            print(f\"   Uniform:     {x:7.4f}\")\n",
    "    \n",
    "    # Alternative Numpy way of doing the above loop (however, without the printing)\n",
    "    #x_uniform[:, iexp] = np.random.uniform(-0.5, 0.5, size=N_uniform) * np.sqrt(12)   \n",
    "    #sum_value += x_uniform[:, iexp].sum()\n",
    "    \n",
    "\n",
    "    # Generating exponential numbers (with mean 0, and RMS of 1):\n",
    "    for i in range(N_exponential): \n",
    "        x = r.exponential() - 1.0                       # Exponential starting at -1. Why?\n",
    "        # x = -np.log(r.uniform()) - 1.0                # Alternative way to produce x exponentially (we'll get to why this works!).\n",
    "        sum_value += x\n",
    "        \n",
    "        x_exponential[i, iexp] = x\n",
    "        if (verbose and iexp == 0 and i < N_verbose) : \n",
    "            print(f\"   Exponential: {x:7.4f}\")\n",
    "    \n",
    "    # Alternative Numpy way of doing the above loop\n",
    "    #x_exponential[:, iexp] = np.random.exponential(size=N_uniform) - 1\n",
    "    #sum_value += x_exponential[:, iexp].sum()\n",
    "\n",
    "    # Generating numbers according to a Cauchy distribution (1 / (1 + x^2)):\n",
    "    for i in range(N_cauchy): \n",
    "        x = np.tan(pi * (r.uniform() - 0.5))              # Cauchy with mean 0... now produced in a \"complicated\" way!\n",
    "        # x = r.standard_cauchy()                         # Alternative way to produce x according to Cauchy PDF.\n",
    "        sum_value += x\n",
    "        x_cauchy[i, iexp] = x\n",
    "        if (verbose and iexp == 0 and i < N_verbose):\n",
    "            print(f\"   Cauchy:      {x:7.4f}\")\n",
    "            \n",
    "    # Alternative Numpy way of doing the above loop\n",
    "    #x_cauchy[:, iexp] = np.random.standard_cauchy(size=N_uniform)\n",
    "    #sum_value += x_cauchy[:, iexp].sum()\n",
    "\n",
    "    N_total = N_uniform + N_exponential + N_cauchy\n",
    "    sum_value = sum_value / np.sqrt(N_total)     # Ask yourself, why I divide by sqrt(N)?\n",
    "    x_sum[iexp] = sum_value\n",
    "    if not (-3.0 < sum_value < 3.0):\n",
    "        N3_sigma += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have created three 2-dimensional arrays of shape (N_pdf, N_experiments), e.g. (10, 1000). We now flatten the arrays to get 1D arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_uniform = x_uniform.flatten()\n",
    "x_exponential = x_exponential.flatten()\n",
    "x_cauchy = x_cauchy.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw the input distributions:\n",
    "\n",
    "It is always important to visualize out data to see if our code produced the expected results, find outliers and just generally get a better understand of it.\n",
    "\n",
    "We first define the number of bins and the ranges of the different distributions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_bins = 100\n",
    "x_ranges = [(-2.5, 2.5), (-1.5, 5.5), (-6, 6)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now loop over our three distributions and plot them (if they are not empty). In each subplot we plot histograms of the distributions with the given number of bins and ranges as defined above and with the mean, standard deviation and truncated standard deviation for each distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3Xl0V9W9///nG6gVCBAQiIpgVAYvQp1yq7byNWi/dUAL4lT8amulpVqtrZfWYu3Podd20UmuturVW63aquCAw636lVst1nqtLV4HiBUMFAWkkDKFMAjC/v3BJ/kmDJKQk+GTPB9rZfH57H2G9wmuvXi5z9knUkpIkiRJkhqvQ0sXIEmSJElthQFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsBSZiLi3yPi/6v1/dKIWBYRVRGxT0vWJkmSJDWH8D1YqhYRCRiUUiqv1XY9MDCldEEDj/UxoBI4NqX0RqaFSpIkSa2UM1hqKkXA3kBZQ3eMbfxvU5IkSXnHf8Sq3iKiNCIWR8TEiFgeEUsj4ku1+u+JiBsjYjAwN9e8OiKez/V/KiL+EhFrcn9+qta+MyPiBxHxErAeODjXdmNE/HfuNsP/jIh9IuL+iKjMHaO4+X4DkiRJ0kczYKmh9gV6AP2A8cCtEdGz9gYppXnAYbmvhSmlEyOiF/AUcAuwD3AT8NR2z2ZdCEwAugHv5to+n2vvBxwCvAz8CugF/BW4LusLlCRJkvaUAUsNtRn4fkppc0rpaaAKGFKP/UYB76SUfp1S+jCl9CDwNnBGrW3uSSmV5fo359p+lVKan1JaAzwDzE8p/S6l9CHwMHBkZlcmSZIkNZIBS7VtAT62XdvH2Baqqq3IhZtq64GCehx7f/7frFS1d9k2M1Vt0U72W1br84adfK/PuSVJkqRmYcBSbe8Bxdu1HcSOwWhPvA8cuF3bAGBJre8uaSlJkqS8ZsBSbdOA70XEARHRISI+w7Zb+B7J4NhPA4Mj4vyI6BQR5wFDgd9mcGxJkiSpVejU0gWoVfl+7uePQE9gPvB/UkpzGnvglNKKiDgduBm4HSgHTk8p/aOxx5YkSZJaC180LEmSJEkZ8RZBSZIkScqIAUtSqxQRd+deaL3TW1Rjm1siojwi3oyIo5q7Rkn5z7FGUtYMWJJaq3uAUz6i/1RgUO5nAtue7ZOkhroHxxpJGTJgSWqVUkp/AFZ+xCajgfvSNn8CCiNiv+apTlJb4VgjKWutYhXB3r17p+Li4pYuQ1IjvPrqq/9IKfVpxlP2o+7LqRfn2pZuv2FETGDb/3mma9euRx966KHNUqCk7DnWSGouezretIqAVVxczKxZs1q6DEmNEBFZvJC6SaSU7gTuBCgpKUmON1L+cqyR1Fz2dLzxFkFJ+WoJ0L/W9wNybZKUJccaSQ1iwJKUr54EvpBb4etYYE1KaYdbdiSpkRxrJDVIq7hFUJK2FxEPAqVA74hYDFwHfAwgpfTvwNPAaUA5sB74UstUKimfOdZIypoBS1KrlFIat5v+BFzWTOVIaqMcayRlzVsEJUmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAagaHHXYYM2fObOkyJEmSJDWxdhOwiouL6dy5MwUFBTU/v/vd7+q9f0RQXl6+R+cuKyujtLR0j/bVzj3xxBMcd9xxfPzjH+eiiy6q07d582bGjx9Pt27dGDBgAA899FCd/ltuuYWioiJ69uzJ1VdfXadv5syZDBkyhK5duzJmzBjWrFlT07d48WJKS0vp0qULRx11FHPmzGmy65MkSVJ+ajcBC+A///M/qaqqqvn5zGc+09IlaQ/16NGDb3/724wfP36HvilTpjB79mwWLVrEfffdx8UXX8yiRYsAeOWVV7j++ut5/vnnmTNnDlOnTq0JYOvXr+ecc87huuuuY/ny5UREnQA2YcIEhg0bxooVKzjvvPM477zzmudiJUmSlDfaVcDalYsuuojx48dz/PHH061bNy6++OKavlNPPZWCggIADj/8cAoKCvjmN79Z079w4UIigocffpgDDzyQbt268aMf/QiA73znOxQUFNChQ4cdZsu2bNnCDTfcQHFxMUVFRVxxxRVs3ry5pv+HP/wh++23H927d6ekpIQlS5Y05a8g75SWljJ27Fh69eq1Q9/DDz/MFVdcQWFhIaWlpRx33HE89thjADzyyCOMHTuWww47jH79+vHlL3+ZqVOnAvD73/+eHj16cP7559O1a1e+9a1vMW3aNAAqKyuZMWMGkyZNonPnzlx55ZW8++67vPnmm8130ZIkSWr1DFg5M2bM4IEHHqCsrIyHH36YP//5zwA888wzVFVVAfDGG29QVVXFv/3bv+2w/0MPPcSrr77KsmXL+OxnPwvAj370I6qqqhgwYMAO299000088cQTvPTSS5SXl/PWW29xyy23ADB37lxuvPFG/vu//5vVq1dz22238fGPf7ypLr3NmTt3LoceeigXXHABU6dOZejQocydO7dO380338zEiRN32vfHP/6Rz372swwcOJCVK1dSUVFBeXk5e++9NwUFBYwYMYIFCxZwyCGH1OwrSZIkQT0CVkT0j4jfR8RbEVEWEd/ItfeKiP+KiHdyf/bMtUdE3BIR5RHxZkQc1dQXUV9jxoyhsLCw5mf58uU1fWeccQYDBgxgwIABDB8+nHnz5jXo2Ndeey29e/emS5cuHHnkkbvd/pe//CXXXXcd/fr1o1u3blx22WVMnz4d2Pa815YtW3j77bfZsmULn/zkJ+ndu3fDLrYdW7duHQUFBcyePZslS5bQrVu3mpBc3Td//nzeeeednfb9/e9/56233qoJtVVVVTV9a9eupaysjFWrVtXZV5IkSYL6zWB9CExMKQ0FjgUui4ihwCTguZTSIOC53HeAU4FBuZ8JwO2ZV72HHn/8cVavXl3z07dv35q+2rea7bXXXmzcuLFBxx40aFCDtl+0aBEXXnhhTdj74he/SEVFBQCDBw/mjjvu4MYbb6RXr16cc845VFZWNuj47VnXrl2pqqrijTfeYOLEiaxdu7bmNs/qvltuuYUnn3xyp31nn302ixcvJqUEQEFBQU1f//79WblyJccdd1ydfSVJkiSoR8BKKS1NKf1P7vNa4K9AP2A0cG9us3uBMbnPo4H70jZ/AgojYr/MK29mEfGR/Z06dWrQ8fr378+zzz5bE/YqKyvrzJpddNFFvPTSS8yfP5958+Zx991371Hd7dHgwYN5++23a76/9dZbDBkyZI/6evXqRZ8+fRg4cCAbNmxg8eLFAGzatIn58+fX7CtJkiRBA5/Biohi4EjgFaAopbQ01/V3oCj3uR+wqNZui3Nt2x9rQkTMiohZ1TM3rdm+++6b6bLcF198Mddeey3vv/8+KSXmzZvHjBkzAJg/fz7PP/88mzZtolOnTmzdupXu3btndu62YMuWLWzcuJEtW7bUfP7www8BOPfcc7nllltYvXo1L7zwAi+//DJnnnkmAOeccw7Tp0+nrKyMJUuWcNddd9WsBnjiiSeyZs0aHnjgAdatW8dPf/pTzj33XAC6d+/OySefzOTJk9m4cSNTpkypuZ1UkiRJqlbvgBURBcCjwDdTSnXuV0vb7qVKDTlxSunOlFJJSqmkT58+Ddl1j51xxhl13oN177337n6nnB/+8Idcfvnl9OvXb4d3J+3Mli1bas7z3nvv1Zz7rrvuAmDixImMGDGC448/nh49ejB27NiaZ8I2bdrEpEmT2GeffRg8eDDHHXccF1544Z5ddBv161//ms6dOzN58mR+85vf0LlzZ2688UYArrzySoYNG0b//v254IILuOuuu+jfvz8AxxxzDNdddx0jR45k2LBhnHvuuTUBq0uXLjz88MNcf/319OnTh61btzJ58uSac95xxx3Mnj2bnj17MnXqVKZNm7bbmU1JkiS1L1H9nMlHbhTxMeC3wLMppZtybXOB0pTS0twtgDNTSkMi4o7c5we3325Xxy8pKUmzZs3K4HIktZSIeDWlVNLSdeyO442U3xxrJDWXPR1v6rOKYAB3AX+tDlc5TwJfzH3+IvBErfYv5FYTPBZY81HhSpIkSZLaivqszPBp4EJgdkS8nmv7LjAZeCgixgPvAufm+p4GTgPKgfXAlzKtWJIkSZJaqd0GrJTSH4FdPWhy0k62T8BljaxLkiRJkvJOg1YRlCRJkiTtmgFLkiRJkjLSLgLW9ddfT0TwzDPPALB69Wo+/vGPU1pa2rKFaY8tXryY0tJSunTpwlFHHdWgd5TdcsstFBUV0bNnzzpL7r/zzjuccsop9OrVi759+/KFL3yBNWvWZHJOSZIktQ/tImABDBo0iEceeQSAJ598kgMPPLCFK1JjTJgwgWHDhrFixQrOO++8mndZ7c4rr7zC9ddfz/PPP8+cOXOYOnUqDz30EACVlZWMGzeOv/3tbyxcuJCNGzcyceLERp9TkiRJ7Ue7CVj//M//zGuvvcaHH37Io48+yllnnVXTt2XLFm644QaKi4spKiriiiuuYPPmzQD87W9/48QTT6Rnz550796dc845h9WrVwOwcOFCIoIpU6bQp08fDjroIF555ZUWub72pLKykhkzZjBp0iQ6d+7MlVdeybvvvsubb765230feeQRxo4dy2GHHUa/fv348pe/zNSpUwE4+uij+eIXv0iPHj3o0qUL559/Pi+//HKjzylJkqT2o90ErIjghBNO4Mknn2TZsmUccsghNX033XQTTzzxBC+99BLl5eW89dZb3HLLLQB88MEHfOUrX2HRokUsXryYFStWcMMNN9Q5dmVlJUuXLmX06NE79Cl75eXl7L333hQUFDBixAgWLFjAIYccwty5c3e779y5czn00EO5+eabmThxIkOHDt3lfi+//DKf+MQnGn1OSZIktR/tJmABnH322XzjG99g1KhRddp/+ctfct1119GvXz+6devGZZddxvTp0wE49NBDGTduHAUFBTUzWG+88Uad/S+99FI6derEaaedxrx585rtetqrdevWUVBQwNq1aykrK2PVqlV069aNqqqqeu87f/583nnnnV3u9+qrr/KrX/2Kf/3Xf230OSVJktR+1OdFw23Gpz71Kf7pn/6Jc845hz/+8Y817YsWLeLCCy+kQ4dteXPr1q3su+++ACxfvpwrrriCF198kXXr1rFp0yaOPvroOsft1asXAHvttRcbN25spqtp+9577z2GDh0KwIEHHkhZWRkAXbt2paqqiv79+7Ny5UoA1q5dS0FBwW6PWb1v9QzlY489tsN+Cxcu5KyzzuK+++5j4MCBjT6nJEmS2o92NYMVEcyYMYNDDz20Tnv//v159tlnWb16NatXr6aysrJmJurqq68mIvjrX//K6tWrmTx5MtvepaymNmDAAKqqqqiqqqoJVwADBw5kw4YNLF68GIBNmzYxf/58hgwZsttjDh48mLfffrvm+1tvvVVnv+XLl3PyySdz4403csopp2RyTkmSJLUf7Spg7crFF1/Mtddey/vvv09KiXnz5jFjxgxg2yxF9+7d6dq1K3/729+44447Wrhade/enZNPPpnJkyezceNGpkyZwoABAxg+fHid7UpLS/nOd75Tp+2cc85h+vTplJWVsWTJEu66666a1QDXrFnDKaecwqWXXsoFF1ywR+eUJElS+2bAAiZOnMiIESM4/vjj6dGjB2PHjmX58uUAXHfddcyaNYsePXpw7rnncsYZZ7RwtQK44447mD17Nj179mTq1KlMmzaNiKizzcKFC1m2bFmdtmOOOYbrrruOkSNHMmzYMM4999yagPX444/z2muv8b3vfY+CgoKan4acU5IkSe1btIbb3UpKStKsWbNaugxJjRARr6aUSlq6jt1xvJHym2ONpOayp+ONM1iSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYDXT99dfvsIS3JEmSJEE7ClhlZWWcdNJJFBYW0rdvX84666w6/TNnzuSAAw5ooerUUIsXL6a0tJQuXbpw1FFHMWfOnHrtt2nTJi677DKKiorYZ599uOqqq+q8OPqOO+6guLiYbt26MWbMGFatWlXTt2HDBr761a/Sq1cvCgsL+drXvpb5dUmSJCm/tZuANWbMGEaNGsWKFSuYN28eJ598ckuXpEaYMGECw4YNY8WKFZx33nk177Landtvv51XXnmFt99+m9dee41p06Zx//33A/D6669z1VVX8dRTT7Fs2TI2b97Mt7/97Zp9r7zySsrLyykrK2PFihUGLEmSJO2gXQSsf/zjH5SXlzN+/Hg6duxIYWEhEyZMqOkvKCjg1FNP5f333695ueyjjz4KwNatW/nWt75F7969GT58OO+9915LXYZyKisrmTFjBpMmTaJz585ceeWVvPvuu7z55pu73feFF17g/PPPp2fPngwYMIBx48Yxffp0AP7whz9wwgkncNhhh9GlSxe+/vWv1/Rt2LCBe++9l5///Ofst99+dOzYkWHDhjXpdUqSJCn/tIuA1atXLwYMGMCECROYOXMmH3zwQZ3+qqoqnnnmGfbff3+qqqqoqqqquYXwoYce4oknnqCsrIzp06fzxBNPtMQlqJby8nL23ntvCgoKGDFiBAsWLOCQQw5h7ty5u913Zy/WLi8v32XfqlWrWLlyJfPmzSMiePzxxykqKmLo0KE89thjjb8YSZIktSntImB16NCB5557jk6dOnHmmWdSVFTEDTfcUK99n3rqKS644AKKiooYNGgQY8aMaeJqtTvr1q2joKCAtWvXUlZWxqpVq+jWrRtVVVW73be0tJQHHniAFStWsHDhQh5++GHWr18PwAknnMDMmTOZPXs2VVVV3H777QCsX7+eyspKNm3axIIFC3j33Xe59dZbufDCC1m6dGmTXqskSZLyS7sIWAADBw7k/vvvZ8WKFTz66KPcfPPNPP3007vdr6Kigr59+9Z8LyoqasoyVQ9du3alqqqK/v37s3LlSo477jjWrl1LQUHBbve99NJLOfbYYxk+fDijRo1izJgxdO/eHYAjjjiCH//4x4wZM4YhQ4ZwzDHHANCtWze6dOnCli1bmDhxInvvvTcjR45k8ODB/OlPf2rSa5UkSVJ+aTcBq1qHDh046aSTKC0tpaysrE77zvTt25fly5fXfF+2bFmT16iPNnDgQDZs2MDixYuBbSsDzp8/nyFDhux237322otf/OIXvP/++5SVlRERHH744TX9l1xyCfPnz2fJkiV84hOfoLi4mB49enDwwQcTEUREnePt7LZCSZIktV/tImBt2bKFa6+9tuYf5G+88QYvvvgiRx99dM02++67LxUVFTsEqNNPP53f/OY3LFu2jHfeeYfHH3+8WWvXjrp3787JJ5/M5MmT2bhxI1OmTGHAgAEMHz68znalpaV85zvfqdO2bt06ysvL2bp1K7/73e/45S9/WWfBk9mzZ7N161bKy8u5+uqrueSSSwDo2bMnJ5xwAjfddBObN2/mxRdfZN68eRx77LFNf8GSJEnKG51auoDm0KFDB8rLyznmmGNYs2YNRUVFfO973+PEE0+s2Wbw4MF8+ctf5vDDD6dTp07cdtttfO5zn+Pss8/mz3/+M4cddhj7778/o0ePZtOmTS14NYJt76u64IIL6NmzJ4ceeijTpk3bYXZp4cKFFBcX12lbv349o0aNYtGiRey///7cdtttHHfccTX911xzDc8//zxdunRh/PjxdZZpv/vuu7nooosoLCykX79+/PrXv2b//fdv0uuUJElSfonWcItTSUlJmjVrVkuXIakRIuLVlFJJS9exO443Un5zrJHUXPZ0vGkXtwhKkiRJUnMwYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSA1UDXX389F1xwQUuX0apFBOXl5S1dhiRJktTs2k3AKisr46STTqKwsJC+ffty1lln1emfOXMmBxxwQObnzbdAli/haPHixZSWltKlSxeOOuoo5syZU6/9Nm3axGWXXUZRURH77LMPV111VZ2XBd9xxx0UFxfTrVs3xowZw6pVq2r6NmzYwFe/+lV69epFYWEhX/va1zK/LkmSJOW3dhOwxowZw6hRo1ixYgXz5s3j5JNPbumS1AgTJkxg2LBhrFixgvPOO4/zzjuvXvvdfvvtvPLKK7z99tu89tprTJs2jfvvvx+A119/nauuuoqnnnqKZcuWsXnz5jrvwbryyispLy+nrKyMFStWGLAkSZK0g3YRsP7xj39QXl7O+PHj6dixI4WFhUyYMKGmv6CggFNPPZX333+fgoICCgoKePTRRwHYunUr3/rWt+jduzfDhw/nvffeq9c5X3zxRQoKCvjhD3/ItGnTao67fPnymm0uuugivvnNbzJ27FgKCgo48MADqays5J577uH444+v2a72LNjChQuJCKZMmUKfPn046KCDeOWVV+qc+7bbbmPIkCF0796dESNG8Pe//x2AZ555hiOOOIJu3bpRVFTENddcUzN7c+qpp1JQUADA4YcfTkFBAd/85jdrjjl79mxKS0vp2bMnRx55JC+//HJN3+uvv86RRx5Jjx49+MEPflCv309jVFZWMmPGDCZNmkTnzp258soreffdd3nzzTd3u+8LL7zA+eefT8+ePRkwYADjxo1j+vTpAPzhD3/ghBNO4LDDDqNLly58/etfr+nbsGED9957Lz//+c/Zb7/96NixI8OGDWvS65QkSVL+aRcBq1evXgwYMIAJEyYwc+ZMPvjggzr9VVVVPPPMM+y///5UVVVRVVVVcwvhQw89xBNPPEFZWRnTp0/niSeeqNc5R4wYQVVVFd/97nc577zzao7bt2/fOtvdd999fOlLX2LNmjU8/vjjfOxjH6vX8SsrK1m6dCmjR4/mhhtuqGl/8MEH+clPfsKjjz7KqlWruOGGG9i8eTOwLSz+/Oc/Z+XKlbz88sv8x3/8B08++SSwLXxVVVUB8MYbb1BVVcW//du/AbB27Vo++9nPcv755/OPf/yDH/zgB4wdO5b169eTUuL888/nC1/4AhUVFaxYsaJe9TdGeXk5e++9NwUFBYwYMYIFCxZwyCGHMHfu3N3uu7MXa1ffErmzvlWrVrFy5UrmzZtHRPD4449TVFTE0KFDeeyxxxp/MZIkSWpTdhuwIuLuiFgeEXNqtU2LiNdzPwsj4vVce3FEbKjV9+9NWXx9dejQgeeee45OnTpx5plnUlRUVCeUfJSnnnqKCy64gKKiIgYNGsSYMWMyrW3kyJGcccYZdOzYkSOPPJLOnTvXa79LL72UTp06cdpppzFv3rya9l/96ldMnDiRYcOG0bFjR0488UT69+8PwKhRoxgxYgQf+9jHOPjggxk5ciRvvPHGbs/129/+ln333ZcJEybQsWNHTjvtNPr27ctLL73EggULmD9/Ppdffjl77bUXV1555Z79Ihpg3bp1FBQUsHbtWsrKyli1ahXdunWrCYgfpbS0lAceeIAVK1awcOFCHn74YdavXw/ACSecwMyZM5k9ezZVVVXcfvvtAKxfv57Kyko2bdrEggULePfdd7n11lu58MILWbp0aZNeqyRJkvJLfWaw7gFOqd2QUjovpXRESukI4FFgeq3u+dV9KaVLsiu1cQYOHMj999/PihUrePTRR7n55pt5+umnd7tfRUVFnVmnoqKiTOsaNGjQHu3Xq1cvAPbaay82btxY075o0SKKi4t3us+f//xnRo4cSZ8+fSgsLOSxxx6rmd36KIsWLaKsrIzCwsKan3feeYelS5dSUVFBYWFhzcxb1r+fnenatStVVVX079+flStXctxxx7F27dqaWxw/yqWXXsqxxx7L8OHDGTVqFGPGjKF79+4AHHHEEfz4xz9mzJgxDBkyhGOOOQaAbt260aVLF7Zs2cLEiRPZe++9GTlyJIMHD+ZPf/pTk16rJEmS8stuA1ZK6Q/Ayp31RUQA5wIPZlxXk+nQoQMnnXQSpaWllJWV1Wnfmb59+9Z5bmrZsmUNPt9H6dSp0w5te++9Nx9++GHN98rKynqfr3///ixcuHCnfePGjWP06NG8//77rF69mtNPP32H2+K2/ZXueMyRI0eyevXqmp/169fzhS98gb59+7J69eqaoNbQ38+eGDhwIBs2bGDx4sXAtpUB58+fz5AhQ3a771577cUvfvEL3n//fcrKyogIDj/88Jr+Sy65hPnz57NkyRI+8YlPUFxcTI8ePTj44IOJiB1+Pzu7rVCSJEntV2OfwRoBLEspvVOr7aCIeC0iXoiIEbvaMSImRMSsiJhVUVHRyDI+2pYtW7j22mtr/kH+xhtv8OKLL3L00UfXbLPvvvtSUVGxQ0A4/fTT+c1vfsOyZct45513ePzxxxt07n333Ze33367TmDanUGDBjF37lwqKytZv349zzzzTL33veiii/jZz37GnDlz2LJlCy+88AKLFi0Ctj1L1bt3bzp16sTMmTN59tlnd1rv9kuejxo1ijlz5vDwww/z4Ycfsm7dOqZPn86qVas46KCDGDhwILfeeiubN29mypQp9a51T3Xv3p2TTz6ZyZMns3HjRqZMmcKAAQMYPnx4ne1KS0v5zne+U6dt3bp1lJeXs3XrVn73u9/xy1/+ss6CJ7Nnz2br1q2Ul5dz9dVXc8kl2yZhe/bsyQknnMBNN93E5s2befHFF5k3bx7HHntsk1+vJEmS8kdjA9Y46s5eLQUGpJSOBP4FeCAiuu9sx5TSnSmlkpRSSZ8+fRpZxkfr0KED5eXlHHPMMRQUFDB27Fi+973vceKJJ9ZsM3jwYL785S9z+OGHc8ABB9Qs/nD22WczevRoDjvsMM466yxGjx7doHOfe+65dO/enX79+nHAAQdQnzB59NFHc/bZZ3P44Ydz5plnctRRR9X7fOeffz7/8i//wtixYyksLOSaa66puX3vtttu45prrqF79+784he/4LTTTtth/x/+8Idcfvnl9OvXj6uvvhrYFmieeeYZ/v3f/52+fftSXFzMb37zGzp06EBE8MADD3DPPffQu3dvmvrvstodd9zB7Nmz6dmzJ1OnTmXatGl70neJAAAeqklEQVQ7zC4tXLhwh8C8fv16Ro0aRUFBAZdccgm33XYbxx13XE1/9e/nU5/6FKeffnqdZdrvvvtu5s6dS2FhIePHj+fXv/41+++/f9NeqCRJkvJK1OcWp4goBn6bUhpWq60TsAQ4OqW0eBf7zQS+lVKa9VHHLykpSbNmfeQmklq5iHg1pVTS0nXsjuONlN8cayQ1lz0dbxozg/UZ4O3a4Soi+kREx9zng4FBwIJGnEOSJEmS8kZ9lml/EHgZGBIRiyNifK7r8+y4uMX/At7MLdv+CHBJSmmnC2Tku+oXCe/sp3fv3i1dniRJkqQWsOMSdttJKY3bRftFO2l7lG3Ltrd51S8SliRJkqRqjV3kQpIkSZKUY8CS1GpFxCkRMTciyiNi0k76B0TE73OvhngzInZcGlOSdsOxRlKWDFiSWqXcgjm3AqcCQ4FxETF0u82+BzyUezXE54HbmrdKSfnOsUZS1gxYklqrTwLlKaUFKaVNwFRg+xfRJaD6XXs9gPebsT5JbYNjjaRMGbAktVb9gEW1vi/OtdV2PXBBRCwGnga+vrMDRcSEiJgVEbPq87JvSe2KY42kTBmwJOWzccA9KaUDgNOAX0fEDuNaSunOlFJJSqmkT58+zV6kpLznWCOp3na7TLt29OnJz7Nk9Qb6FXbmpUkntnQ5eaH6dwb4e1N9LQH61/p+QK6ttvHAKQAppZcjYm+gN7C8WSqU1BY41kjKlDNYe2DJ6g0snDyqJjBo96p/Z/7e1AB/AQZFxEERsRfbHix/crtt3gNOAoiIfwL2BrwvR1JDONZIypQBS1KrlFL6ELgceBb4K9tW8CqLiO9HxOdym00EvhIRbwAPAhellFLLVCwpHznWSMqatwhKarVSSk+z7YHy2m3X1vr8FvDp5q5LUtviWCMpSwasjPhcliRJkiRvEcyIz2VJkiRJcgZLaiWcBZUkScp/zmBJrYSzoJIkSfnPgCVJkiRJGfEWQamBvJVPkiRJu+IMltRA3sonSZKkXXEGS7tVPWMD5MWsjTNMkiRJainOYGm3qmds8mXWxhkmSZIktRQDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcRl2ptR1suduxy5JEmS1Lo4g9WMsl7u3OXIJUmSpNbFgCVJkiRJGTFgSZIkSVJGfAarhX3Uc1lZP7Ol5uUzcpIkSe2PAauFVT9HBVA86al692XFENB0qv/+murvTpIkSa2Ptwi2cy6UIUmSJGXHgCVJkiRJGTFgSZIkSVJGfAZLrZILfEiSJCkfGbDUKjXHAh+SJElS1nYbsCLibuB0YHlKaViu7XrgK0BFbrPvppSezvVdDYwHtgBXpJSebYK6m01DV9nrV9iZ4klPtYlZl9rXXt/t8v2aG6L677r6c3u6dkmSJO1cfWaw7gF+Ady3XfuUlNJPazdExFDg88BhwP7A7yJicEppSwa1toiGLrVd/Y/stjDrUnsWqT7btYVrbojagaq9XbskSZJ2brcBK6X0h4gorufxRgNTU0ofAH+LiHLgk8DLe1yhGqS9zia1Vf59SpIk5ZfGrCJ4eUS8GRF3R0TPXFs/YFGtbRbn2nYQERMiYlZEzKqoqNjZJtoD+fheq09Pfp7iSU9RPOkpPj35+ZYup1XJx79PSZKk9mxPA9btwCHAEcBS4GcNPUBK6c6UUklKqaRPnz57WIbaguoQYZCQJElSvtujVQRTSsuqP0fEfwC/zX1dAvSvtekBuTZpB/VdRCPrc3mrnSRJkprKHs1gRcR+tb6eCczJfX4S+HxEfDwiDgIGAX9uXIlqq6pnrpoj8HirnSRJkppDfZZpfxAoBXpHxGLgOqA0Io4AErAQ+CpASqksIh4C3gI+BC7L5xUEJUmSJKkh6rOK4LidNN/1Edv/APhBY4qSJEmSpHy0R89gqeVUP0sENMuzS9trSy9SliRJkrJmwMozO3v5b3Xoqf68K1ks9NDQFylnsZCFoU6SJEn5woDVBtQ3dFSHs/qGoyzsLBA2VENDnSRJktRSDFiNUN+Zo4Yer6lnarKuu7lsX7ezWZIkSWptDFiNkPU/8JtrpiZfg0ntup3NkiRJUmtkwNIeyfdZMGfAJEmS1BQMWNoj+RpOfJ5LkiRJTcmApTpaehl4SZIkKZ8ZsFRHFqv+SZIkSe2VAasZZPEuKH20rGfesnhn2J7K1+fbJEmSZMBqFg2dFaq9EMOeai+LOdQOQlnOvLXEO8OqteW/L0mSpLbOgNUKZfEP7PaymIO3NEqSJKk1MWApM2191sxbPSVJkrQ7Bixlpq3PmjlbJkmSpN0xYDWR1rJQQWupo7YsnjHLR9svxNEWZ/kkSZLaOwNWE2kt/3huLXXU1hprag61Z8Da6iyfJElSe2fAUqvSXLNb28/stdfQJ0mSpGwZsNSqNFfQqX0eZ5MkSZKUlQ4tXYDU0qpnsz49+fmWLkWSJEl5zhksAdnemtcaF9b4KG199UNJkiQ1HwOWgGxvzfN5JkmSJLVX3iKovJTvS717W6IkSVLb5AyW8lK+z5J5W6IkSVLbZMBSm5Vvz4JJkiQp/xmw1Gbl+yyXJEmS8o/PYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJWXAZdclSZIELnIhZcJl1yVJkgQGrHrLesnvfH9RriRJkqQdGbDqKeslv11CXJIkSWp7fAZLakHVM5m7m81sj894RcQpETE3IsojYtIutjk3It6KiLKIeKC5a5TUNjjeSMqSM1hSC6rvTGZ7e8YrIjoCtwL/G1gM/CUinkwpvVVrm0HA1cCnU0qrIqJvy1QrKZ853kjKmgFLylDWz+q1Y58EylNKCwAiYiowGnir1jZfAW5NKa0CSCktb/YqJbUFjjeSMmXAyhMuipEffLYuM/2ARbW+LwaO2W6bwQAR8RLQEbg+pfR/d3awiJgATAAYMGBA5sVKymuZjTeONZKgHgErIu4GTgeWp5SG5dp+ApwBbALmA19KKa2OiGLgr8Dc3O5/Sild0gR1tzv+w13aQSdgEFAKHAD8ISKGp5RWb79hSulO4E6AkpKS1JxFSmoT6jXeONZIgvotcnEPcMp2bf8FDEspfQKYx7b7kqvNTykdkfsxXClv1HfBiT09rrOPDbIE6F/r+wG5ttoWA0+mlDanlP7GtrFoUDPVJ6ntcLyRlKndzmCllP6Qm5mq3Taj1tc/AWdnW5bU/JpqltDZxz3yF2BQRBzEtn/ofB44f7ttHgfGAb+KiN5su4VnQbNWKaktcLyRlKkslmm/GHim1veDIuK1iHghIkbsaqeImBARsyJiVkVFRQZlSGorUkofApcDz7LttuOHUkplEfH9iPhcbrNngRUR8Rbwe+DbKaUVLVOxpHzleCMpa41a5CIirgE+BO7PNS0FBqSUVkTE0cDjEXFYSqly+329T1nSR0kpPQ08vV3btbU+J+Bfcj+StMccbyRlaY9nsCLiIrYtfvF/cgMPKaUPqv+PTkrpVbYtgDE4gzolSZIkqdXbo4AVEacAVwGfSymtr9XeJ/fCPiLiYLY9AOo9ypIkSZLahfos0/4g25Yl7R0Ri4Hr2LZq4MeB/4oI+H/Lsf8v4PsRsRnYClySUlrZRLVLkiRJUqtSn1UEx+2k+a5dbPso8Ghji5Jai+ol1qs/S5IkSR+lUYtcSG2dS6xLkiSpIbJYpl2SJEmShAFLkiRJkjJjwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMlKvgBURd0fE8oiYU6utV0T8V0S8k/uzZ649IuKWiCiPiDcj4qimKl6SJEmSWpP6zmDdA5yyXdsk4LmU0iDgudx3gFOBQbmfCcDtjS9TkiRJklq/egWslNIfgJXbNY8G7s19vhcYU6v9vrTNn4DCiNgvi2IlSZIkqTVrzDNYRSmlpbnPfweKcp/7AYtqbbc41yZJkiRJbVomi1yklBKQGrJPREyIiFkRMauioiKLMiRJkiSpRTUmYC2rvvUv9+fyXPsSoH+t7Q7ItdWRUrozpVSSUirp06dPI8qQJEmSpNahMQHrSeCLuc9fBJ6o1f6F3GqCxwJrat1KKEmSJEltVqf6bBQRDwKlQO+IWAxcB0wGHoqI8cC7wLm5zZ8GTgPKgfXAlzKuWZIkSZJapXoFrJTSuF10nbSTbRNwWWOKkiRJkqR8lMkiF5IkSZIkA5YkSZIkZcaAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSWq1IuKUiJgbEeURMekjtjsrIlJElDRnfZLaBscaSVkyYElqlSKiI3ArcCowFBgXEUN3sl034BvAK81boaS2wLFGUtYMWJJaq08C5SmlBSmlTcBUYPROtvtX4EfAxuYsTlKb4VgjKVMGLEmtVT9gUa3vi3NtNSLiKKB/SumpjzpQREyIiFkRMauioiL7SiXlM8caSZkyYEnKSxHRAbgJmLi7bVNKd6aUSlJKJX369Gn64iS1GY41khrKgCWptVoC9K/1/YBcW7VuwDBgZkQsBI4FnvThc0kN5FgjKVOdWrqA1urTk59nyeoN9Cvs3NKlSO3VX4BBEXEQ2/6x83ng/OrOlNIaoHf194iYCXwrpTSrmeuUlN8cayRlyoC1C0tWb2Dh5FEtXYbUbqWUPoyIy4FngY7A3Smlsoj4PjArpfRky1YoqS1wrJGUNQOWpFYrpfQ08PR2bdfuYtvS5qhJUtvjWCMpSz6DJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGXKZdyiP9CjtTPOmpms8vTTqxhSuSJElSbQYsKY/UDlTVQUuSJEmtxx4HrIgYAkyr1XQwcC1QCHwFqMi1fzf3Aj9JkiRJatP2OGCllOYCRwBEREdgCfAY8CVgSkrpp5lUKEmSJEl5IqtFLk4C5qeU3s3oeJIkSZKUd7IKWJ8HHqz1/fKIeDMi7o6InjvbISImRMSsiJhVUVGxs00kSZIkKa80OmBFxF7A54CHc023A4ew7fbBpcDPdrZfSunOlFJJSqmkT58+jS1DkiRJklpcFjNYpwL/k1JaBpBSWpZS2pJS2gr8B/DJDM4hSZIkSa1eFgFrHLVuD4yI/Wr1nQnMyeAckiRJktTqNeo9WBHRFfjfwFdrNf84Io4AErBwuz5JkiRJarMaFbBSSuuAfbZru7BRFUmSJElSnspqFUFJkiRJavcMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZaRTYw8QEQuBtcAW4MOUUklE9AKmAcXAQuDclNKqxp5LkiRJklqzrGawRqaUjkgpleS+TwKeSykNAp7LfZckSZKkNq2pbhEcDdyb+3wvMKaJziNJkiRJrUYWASsBMyLi1YiYkGsrSiktzX3+O1CUwXkkSZIkqVVr9DNYwPEppSUR0Rf4r4h4u3ZnSilFRNp+p1wYmwAwYMCADMqQJEmSpJbV6BmslNKS3J/LgceATwLLImI/gNyfy3ey350ppZKUUkmfPn0aW4YkSZIktbhGBayI6BoR3ao/A58F5gBPAl/MbfZF4InGnEeSJEmS8kFjbxEsAh6LiOpjPZBS+r8R8RfgoYgYD7wLnNvI80iSJElSq9eogJVSWgAcvpP2FcBJjTm2JEmSJOWbplqmXZIkSZLaHQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDlqRWKSJOiYi5EVEeEZN20v8vEfFWRLwZEc9FxIEtUaek/Od4IylLBixJrU5EdARuBU4FhgLjImLodpu9BpSklD4BPAL8uHmrlNQWON5IypoBS1Jr9EmgPKW0IKW0CZgKjK69QUrp9yml9bmvfwIOaOYaJbUNjjeSMmXAktQa9QMW1fq+ONe2K+OBZ3bVGRETImJWRMyqqKjIqERJbURm441jjSQwYEnKcxFxAVAC/GRX26SU7kwplaSUSvr06dN8xUlqU3Y33jjWSALo1NIFSNJOLAH61/p+QK6tjoj4DHANcEJK6YNmqk1S2+J4IylTzmBJao3+AgyKiIMiYi/g88CTtTeIiCOBO4DPpZSWt0CNktoGxxtJmTJgSWp1UkofApcDzwJ/BR5KKZVFxPcj4nO5zX4CFAAPR8TrEfHkLg4nSbvkeCMpa94iKKlVSik9DTy9Xdu1tT5/ptmLktQmOd5IypIzWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZWSPA1ZE9I+I30fEWxFRFhHfyLVfHxFLIuL13M9p2ZUrSZIkSa1Xp0bs+yEwMaX0PxHRDXg1Iv4r1zclpfTTxpcnSZIkSfljjwNWSmkpsDT3eW1E/BXol1VhkiRJkpRvMnkGKyKKgSOBV3JNl0fEmxFxd0T03MU+EyJiVkTMqqioyKIMSZIkSWpRjQ5YEVEAPAp8M6VUCdwOHAIcwbYZrp/tbL+U0p0ppZKUUkmfPn0aW4YkSZIktbhGBayI+BjbwtX9KaXpACmlZSmlLSmlrcB/AJ9sfJmSJEmS1Po1ZhXBAO4C/ppSuqlW+361NjsTmLPn5UmSJElS/mjMKoKfBi4EZkfE67m27wLjIuIIIAELga82qkJJkiRJyhONWUXwj0DspOvpPS9HkiRJkvJXJqsISpIkSZIMWJIkSZKUGQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZabKAFRGnRMTciCiPiElNdR5JbdfuxpGI+HhETMv1vxIRxc1fpaR851gjKUtNErAioiNwK3AqMBQYFxFDm+Jcktqmeo4j44FVKaWBwBTgR81bpaR851gjKWtNNYP1SaA8pbQgpbQJmAqMbqJzSWqb6jOOjAbuzX1+BDgpIqIZa5SU/xxrJGWqUxMdtx+wqNb3xcAxtTeIiAnAhNzXqoiY20S1APQG/tHQnaJ1/f+pPbqGViTf64dWeA0N/G+0qes/MOPj7XYcqb1NSunDiFgD7MN217ndePNBRMzJuNbm1Or+O2ygfK8f8v8a8r3+IRkfz7Fm5/L9vxPI/2uw/pa3R+NNUwWs3Ur/f3t3FypXdYZx/P9omopotRhE0UgUIhjsRUIQBamWiIQgBvGDCKJiKCRFUXohQm5Er6TohSBoQGkrtpVqqwe0WPwIAWnSBD9iNVWSVDQ2qGgbBTFVfHuxl2FIzjmzZs7MXnudeX6wYc85c5jn3bPmndn7rNk7YjOwuY3HkrQzIla28VjjUnsNteeH+muoPf9c9Pab2reD85dXew3zIX/pDDNxr+mW2mtw/vKG7TfjmiL4EbC45/aZ6WdmZrly+sjh+0haAJwEfNZKOjObL9xrzGykxrWDtQNYKulsSQuBdcDUmB7LzOannD4yBdyU1q8BXo6IaDGjmdXPvcbMRmosUwTT/ORbgReAY4HHIuLtcTxWplamIo5Z7TXUnh/qr6Gq/DP1EUn3ADsjYgp4FHhc0h7gc5oPRv1UtR2m4fzl1V6D8/dwr5lR7fmh/hqcv7yhapAPwJiZmZmZmY3G2C40bGZmZmZmNmm8g2VmZmZmZjYiE7ODJelXkv4paZekP0s6uXSmQUm6VtLbkr6TVM1pLyWtlvSupD2S7iqdZxCSHpP0Sa3XMpG0WNIrkt5JY+f20pna0G/MSfqhpCfT77dLWtJ+ytll1PDL9LzukvSSpFFfh2xOcl/3kq6WFF3raTn5JV3X89r6XdsZ+8kYQ2el/vB6GkdrSuScSb/+q8aDqb5dkla0nTHlqLrfuNeUV3u/ca+ZRkRMxAJcDixI6/cB95XONEQN59Fc8GwLsLJ0nszMxwJ7gXOAhcCbwLLSuQbI/1NgBfCP0lmGzH86sCKtnwi8V9P2H7LmvmMO+AXwcFpfBzxZOvcQNfwMOD6tb+xSDbmv+zQmtwLbutTTMrf/UuB14Mfp9qmlcw9Rw2ZgY1pfBrxfOvcR+Wbtv8Aa4C+AgAuB7R3dzp3tN+415Zfa+417zfTLxPwHKyL+GhHfppvbaK5zUZWI2B0R75bOMaALgD0RsS8i/gf8AVhbOFO2iNhKc8aoKkXEgYh4La1/CewGziibauxyxtxa4Ddp/SlglSS1mLGfvjVExCsR8VW62bWelvu6v5fmgNfXbYbLkJP/58BDEfEfgIj4pOWM/eTUEMCP0vpJwL9bzNdXRv9dC/w2GtuAkyWd3k66w2rvN+415dXeb9xrpjExO1hHuIVmT9TG7wzgw57b+5n/H/A7KU1LWQ5sL5tk7HLG3OH7pAMvB4FTWkmXZ9DXzXq61dP65k9TLBZHxHNtBsuUs/3PBc6V9KqkbZJWt5YuT04NdwM3SNoPPA/c1k60kenC+0vt/ca9prza+417zTTGch2sUiS9CJw2za82RcSz6T6bgG+BJ9rMliunBrNBSToBeBq4IyK+KJ3HRkfSDcBK4JLSWXJJOgZ4ALi5cJS5WEAzbedSmiP6WyX9JCL+WzTVYK4Hfh0R90u6iOY6T+dHxHelg1n3uNcUVXu/mbheM692sCListl+L+lm4ApgVaRJlV3Tr4YKfQQs7rl9ZvqZtUTSD2h2rp6IiD+VztOCnDH3/X32S1pAM2Xhs3biZcl63Ui6DNgEXBIRh1rKlqNf/hOB84EtaabUacCUpCsjYmdrKWeWs/3308zD/wb4l6T3aD4A7WgnYl85NawHVgNExN8kHQcsAro0/Wg2XXh/qb3fuNeUV3u/ca+ZxsRMEUz/Tr0TuLJnLrGN3w5gqaSzJS2k+YLvVOFMEyPN838U2B0RD5TO05KcMTcF3JTWrwFe7thBl741SFoOPELT07r2JjVr/og4GBGLImJJRCyh+V5Hlz7w5IyhZ2iOJiNpEc0Unn1thuwjp4YPgFUAks4DjgM+bTXl3EwBN6YzfF0IHIyIAy1nqL3fuNeUV3u/ca+ZzlzOulHTAuyhmT/5RloeLp1piBquojmKcQj4GHihdKbM3Gtozl63l2aqY/FMA2T/PXAA+CZt+/WlMw2Y/2KaL5fu6hn7a0rnaqHuo8YccA/NGys0zf2PqS/8HTindOYhangx9YHvn9ep0pkHyX/EfbfQvTN79dv+opl69A7wFrCudOYhalgGvEpz1q83gMtLZz4i/1H9F9gAbOh5Dh5K9b1VagzV3m/ca8ovtfcb95qjF6U/NDMzMzMzszmamCmCZmZmZmZm4+YdLDMzMzMzsxHxDpaZmZmZmdmIeAfLzMzMzMxsRLyDZWZmZmZmNiLewTIzMzMzMxsR72CZmZmZmZmNyP8BocXZ/scoiG8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(12,6))\n",
    "\n",
    "x_all = [x_uniform, x_exponential, x_cauchy]\n",
    "titles = ['Uniform', 'Exponential', 'Cauchy']\n",
    "\n",
    "for ax_i, x, title, x_range in zip(ax, x_all, titles, x_ranges):\n",
    "    \n",
    "    if len(x) > 0:\n",
    "        \n",
    "        ax_i.hist(x, bins=N_bins, range=x_range, histtype='step')\n",
    "        ymax = ax_i.get_ylim()[1]*1.2\n",
    "        ax_i.set(title=title, ylim=(0, ymax), xlim=x_range)\n",
    "\n",
    "        d = {'Entries': len(x),\n",
    "             'Mean': x.mean(),\n",
    "             'Std': x.std(ddof=1),\n",
    "             'Std_truncated': x[(x_range[0]<x) & (x<x_range[1])].std(ddof=1),\n",
    "        }\n",
    "\n",
    "        text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "        add_text_to_ax(0.1, 0.97, text, ax_i, fontsize=12)\n",
    "\n",
    "fig.tight_layout()\n",
    "\n",
    "if save_plots:\n",
    "    fig.savefig('Plot_CentralLimit_Input.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Outside $\\pm 3 \\sigma$\n",
    "\n",
    "Below we write how many of the experiments had a result outside the range [-3,3], i.e. beyond $\\pm 3 \\sigma$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Number of experiments beyond 3 sigma / total:    0 / 1000 = 0.0000\n"
     ]
    }
   ],
   "source": [
    "print(f\"\\nNumber of experiments beyond 3 sigma / total: {N3_sigma:4d} / {N_experiments:d} = {N3_sigma/N_experiments:6.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- __Question__: *Do the numbers above make sense? What did you expect it to be?*\n",
    "***\n",
    "## Draw output plots with corresponding fits to the screen:\n",
    "\n",
    "Having visualized the individual distributions above, we now visualize the sum of the different distributions, `x_sum`. We first plot a histogram of it in the range from `xmin` to `xmax`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAGDCAYAAADtZ0xmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xu4ZGdZJ+zfQ4DQEKDFhIidxA7nicrJBmEyqCTioOE0jsPAIBMwYzwiOMxIgzjyjTpf6ygHBWcIgkRBgQ85aQQ5BDxECSScQgiYGDuSTiAJpk0iLUnI8/1Rq8lO9z7U6uzatXf3fV/XvnatVWvVemqtqr1/9da71lvdHQAAYHp3mHcBAACw0QjRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDRtEVV1YVd8z7zrmqar+XVV9oapuqKqHz7ue26OqDq+qz1bVfZZZZmdVfe9a1rWWxh7Pqvqeqrp8RrVsraquqjsO0++pqlNX6bEfW1WfXzC9qsd11n8bquohVfXXs3p82KiEaFgHFvunWlXPrqq/2jvd3d/a3R9e4XFuEwQOQr+e5Ke7+4ju/sS8i7mdTk/yF919ZZJU1Ruq6pfnXNNaW/Z4Dq/l+8+hrnT393f3mSstN02N3f2X3f2g1ahrsdfJNH8bbo/u/nSS3VX1pFltAzYiIRqY2joI59+S5MI517BafjzJ78+7iDk7mI7notbBe2a1vCnJj827CFhPhGjYIBa2VlfVo6rqvKq6rqq+VFUvGxb7i+H37uEr8sdU1R2q6iVVdVlVXVVVv1dV91zwuP95uO/LVfUL+2znpVX1tqp6Y1Vdl+TZw7b/pqp2V9WVVfWqqrrzgsfrqvrJqrq4qq6vql+qqvtV1V8P9b514fL7PMdFax26PtyQ5LAkn6qqv1tk3X9dVddU1bHD9EOr6tqqevAK+/WFVbVrqPXzVXXyMP82LX77diUY9tN/r6pPV9U/V9XrqurooRvA9VX1gar6hiW2eVyS+yY5d5g+Pckzk/zccNz+eMHiDxu28U9V9ZaqusuCx3liVX1yOBZ/XVUPWWJ7VVUvH/bpdVV1QVV923Dfh6vqvyxY9jbfgMzxeO59LX9q2Cf/ccF9Lxge78qqes6C+YdX1a9X1T/U5H3xf6tq0xK1HTYse01VXZrklH3u//p+qar7V9WfD8fgmqp6y1I17n2dDK+rLyb53X1fO4NH1qQ7z7VV9bt7j+u++3/BMbj/Uq+Tuu179vCqekVVXTH8vKKqDh/u21vbUvvvB4aarh/eE/9tQRkfTnLy3scChGjYqF6Z5JXdfY8k90vy1mH+dw2/Nw9fkf9NkmcPP4/LJLgdkeRVSVJVJyT57Uz+Md8nyT2TbNlnW09J8rYkmzNpjfpakp9NcmSSxyQ5OclP7rPOv03yHUkeneTnkpyR5IeTHJvk25I8Y4nntWit3f3V7j5iWOah3X2/fVfs7r9O8pokZw7B6Y1JfqG7P7fEtlJVD0ry00ke2d13H+reudTyi/j3SR6f5IFJnpTkPUlenOSoTP6+/swS6317kku7++ah9jMy2be/Nhy3hV+bPy3JE5Icn+Qhmeyf1KQP8eszaR38xuG5v3uJkPN9mbw2HpjJMX5aki+PeJ7zOJ7fteD+I7r7LcP0N+XW1+lpSV694MPKjuE5PizJ/Ydl/scStf1okicmeXiSbUl+aOmnn19K8r4k35DkmCS/NUWN98qkpf30JR7zmZns1/sNNb9kme1n2N5yr5O9fj6T4/SwJA9N8qh9Hnu5/fe6JD82vBe+LcnZC7a9K8lNSValWwocDIRoWD/eObQo7q6q3ZmE26XclOT+VXVkd9/Q3R9ZZtlnJnlZd1/a3TckeVGSp9fka+YfSvLH3f1X3X1jJoGj91n/b7r7nd19S3fv6e7zu/sj3X1zd+/MJLx99z7r/Fp3X9fdFyb5TJL3Ddv/p0yC5lInkS1X6zRemklA+GiSXUlevcLyX0tyeJITqupO3b2zu/drFV3Gb3X3l4aA8ZdJzu3uT3T3vyR5R5Z+npuTXD/lNn6zu6/o7n9M8seZhKNkEs5e093ndvfXhv67X80kQO3rpiR3T/LgJNXdF+3tiz2leR3PxdyU5H92903d/adJbkjyoKqqTPbJz3b3P3b39Un+V5KnL/E4T0vyiu7+wrBv/98VtvktSb65u/+lu/9qmWWT5JYkvzh8WNizxDKvWrDtX8nSH0TGemYm++eq7r46yf+T5FkL7l90/y2474Squkd3X9vdH9/nsa/P5LULRIiG9eSp3b1570/2b91d6LRMWq8+V1Ufq6onLrPsNye5bMH0ZUnumOTo4b4v7L2ju7+S/Vsov7BwoqoeWFV/UlVfrEkXj/+VSav0Ql9acHvPItNHZHHL1bqi7r4pyRsyaUX7je7e9wPBvstfkuT5mYTvq6rqzVX1zdNsa3Cgz/PaTELtNL644PZXFjzmtyR5wT4fvI7NZB/eRnefncm3D6/O5HmeUVX3mHL7yZyO5xK+vLcFf7B3nxyV5K5Jzl+wP947zF+qtoWv7cuWWC6ZtL5Xko/W5EoYP7JCjVcPH6SWs++2x7zulrPYPl/42Evtv2TyzcoPJLls6L7ymH0e++5Jdq9SnbDhCdGwAXX3xd39jCT3TvKrSd5WVXfL/q3ISXJFJoFrr+OS3JxJELoyk6+nkyRDN4hv3Hdz+0z/nySfS/KAoTvJizMJGKthuVpXVFVbkvxikt9N8hvT9N/s7j/o7n8zbLcz2Z9J8s+ZhLK9vmmaGqb06STH79Miu2zgX8QXkvzKwg9e3X3X7v7DxRbu7t/s7u9IckImH8D++3DXLJ/n7TqeI12TSaD/1gX7454Luo3s68pMPnQsrG1R3f3F7v7R7v7mTLrP/HYtf0WOaY7lvtu+Yrh9m+NRVfsej5Uee7F9fsUSy972gbs/1t1PyeTvyjtzazexve+tOyf5/BKrwyFHiIYNqKp+uKqO6u5bcmvL0C1Jrh5+33fB4n+Y5Ger6viqOiKTluO3DK1Rb0vypJqclHfnTFpkVwrEd09yXZIbanLS3k+s1vNaodZlDV/nvyGTfp2nZRKSfmmFdR5UVScNYftfMglhtwx3fzLJD1TVvYYg8/wDfE776e7Lk1ySSX/Vvb6U2x63lbw2yY9X1XfWxN2q6pSq2q+Fu6oeOSx3p0xC2r/kts/zB6vqrkMwPO1AntMSDvh4DqbeJ8N74bVJXl5V904mwa+q/u0Sq7w1yc9U1TFDn+DtSz12Vf2Hqtr7YfPaTILs3v039rjt9VPDtu+VST/mvf2pP5XkW6vqYcPJhi/dZ72VtveHSV5SVUdV1ZGZdNF640rFVNWdq+qZVXXP4Rud63Lrc0wmXbbO7u6vTvPk4FAgRMPG9IQkF9bkCgevTPL0ob/yVzLpX3nO8JX2ozM5+ez3M7lyx99nEqCemyRDH9fnJnlzJqHzhiRXZdK3din/Lcl/yqR/5Gtz6z//1bBkrVP4mUxa0H5h6MbxnCTPqarHLrPO4ZmcjHZNJt0m7p1Jv90MdXwqkxMN35fVfZ7JpC/5wr6qr8ukP+ruqnrnSit393mZnBz3qkyC3SUZTjpcxD0yOVbXZvL1/peT/O/hvpcnuTGTcHZmJieurZbbczyTSYA8c9gnT5ti+Rdmsh8+MnQ1+kCWPhHutUn+LJNj/PEkb1/mcR+Z5Nzh/fbuJM/r7ksPsMa9/iCT19WlSf4uyS8nSXf/bZL/OdR+cZJ9+1+v9Dr55STnZfJtxwXDc5v2+uPPSrJz2Hc/nkn/6r2emeT/Tvk4cEioFboMAoeQobVwdyZdNf5+3vUczIbW708kOXnkSX6wpmpy6cTXdPe+faThkCZEwyGuJqOQfTCTbhy/keQ7kzxipZPyAOBQpjsH8JRMTjy6IskDMukactAE6Ko6riYDUyz2s+TJZACwHC3RAAAwkpZoAAAYSYgGAICRbs/Qq2vmyCOP7K1bt867DAAADmLnn3/+Nd291Eint7EhQvTWrVtz3nnnzbsMAAAOYlV12bTL6s4BAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAj3XHeBQBwqxN3nJ1du/fsN3/L5k05Z/tJc6gIgMUI0QDryK7de7Jzxyn7zd+6/aw5VAPAUnTnAACAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRjFgIcJAyhDjA7AjRAAcpQ4gDzI7uHAAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACPdcZYPXlU7k1yf5GtJbu7ubVV1ryRvSbI1yc4kT+vua2dZBwAArKa1aIl+XHc/rLu3DdPbk3ywux+Q5IPDNAAAbBjz6M7xlCRnDrfPTPLUOdQAAAAHbNYhupO8r6rOr6rTh3lHd/eVw+0vJjl6sRWr6vSqOq+qzrv66qtnXCYAAExvpn2ik/yb7t5VVfdO8v6q+tzCO7u7q6oXW7G7z0hyRpJs27Zt0WUAAGAeZtoS3d27ht9XJXlHkkcl+VJV3SdJht9XzbIGAABYbTML0VV1t6q6+97bSb4vyWeSvDvJqcNipyZ516xqAACAWZhld46jk7yjqvZu5w+6+71V9bEkb62q05JcluRpM6wBAABW3cxCdHdfmuShi8z/cpKTZ7VdAACYtVmfWAjAKtiyeVO2bj9ryfvO2X7SGlcEcGgTogE2gOVC8lLhGoDZmcdgKwAAsKEJ0QAAMJIQDQAAIwnRAAAwkhMLATa4pa7csWXzpjlUA3BoEKIBNjiXtwNYe7pzAADASEI0AACMJEQDAMBI+kQDzNCJO87Ort179ptvqG6AjU2IBpihXbv3ZOeOU/abb6hugI1Ndw4AABhJiAYAgJGEaAAAGEmIBgCAkZxYCDAHG2moblcYAdifEA0wBxspfLrCCMD+dOcAAICRhGgAABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGMWAiwCpYbGnu9WW7I8Y00kiLAPAnRAKtgqaGx16OlgrJhvAGmpzsHAACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEiG/QZgVZ244+zs2r1nv/lbNm9acshxgI1GiAZgVe3avSc7d5yy3/yt28+aQzUAs6E7BwAAjCREAwDASEI0AACMJEQDAMBIQjQAAIw08xBdVYdV1Seq6k+G6eOr6tyquqSq3lJVd551DQAAsJrWoiX6eUkuWjD9q0le3t33T3JtktPWoAYAAFg1Mw3RVXVMklOS/M4wXUlOSvK2YZEzkzx1ljUAAMBqm3VL9CuS/FySW4bpb0yyu7tvHqYvT7JlxjUAAMCqmlmIrqonJrmqu88/wPVPr6rzquq8q6++epWrAwCAAzfLlugTkzy5qnYmeXMm3ThemWRzVe0dbvyYJLsWW7m7z+jubd297aijjpphmQAAMM7MQnR3v6i7j+nurUmenuTs7n5mkg8l+aFhsVOTvGtWNQAAwCzM4zrRL0zyX6vqkkz6SL9uDjUAAMABu+PKi9x+3f3hJB8ebl+a5FFrsV0AAJgFIxYCAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMtCaXuANg/duyeVO2bj9r0fmrsTzAwUSIBiBJcs72k2a6PMDBRHcOAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRhGgAABjpjvMuAGAjOXHH2dm1e89+87ds3jSHagCYFyEaYIRdu/dk545T5l0GAHOmOwcAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjDTViIVV9e3dfcGsiwFYa8sN433O9pPmUBEAG8G0w37/dlUdnuQNSd7U3f80u5IA1s5Sw3hv3X7WHKoBYKOYqjtHdz82yTOTHJvk/Kr6g6p6/EwrAwCAdWrqPtHdfXGSlyR5YZLvTvKbVfW5qvrBWRUHAADr0VQhuqoeUlUvT3JRkpOSPKm7/9Vw++UzrA8AANadaftE/1aS30ny4u7++hk43X1FVb1kJpUBAMA6NW2IPiXJnu7+WpJU1R2S3KW7v9Ldvz+z6gDmZMvmTYueXLhl86Y5VAPAejNtiP5Aku9NcsMwfdck70vyr2dRFMC8ubwdAMuZ9sTCu3T33gCd4fZdZ1MSAACsb9OG6H+uqkfsnaiq70iy/+gEAABwCJi2O8fzk/x/VXVFkkryTUn+48yqAgCAdWyqEN3dH6uqByd50DDr89190+zKAgCA9WvalugkeWSSrcM6j6iqdPfvzaQqAA4ZJ+44O7t2799DcMvmTYue4Dl2eYBZmCpEV9XvJ7lfkk8m+dowu5MI0QDcLrt278nOHafsN3+xSwweyPIAszBtS/S2JCd0d0/7wFV1lyR/keTwYTtv6+5frKrjk7w5yTcmOT/Js7r7xnFlAwDA/Ex7dY7PZHIy4RhfTXJSdz80ycOSPKGqHp3kV5O8vLvvn+TaJKeNfFwAAJiraVuij0zy2ar6aCbhOEnS3U9eaoWh1XrvtaXvNPx0kpOS/Kdh/plJXprk/4yqGgAA5mjaEP3SA3nwqjosky4b90/y6iR/l2R3d988LHJ5ki1LrHt6ktOT5LjjjjuQzQMAwExM1Z2ju/88yc4kdxpufyzJx6dY72vd/bAkxyR5VJIHT1tYd5/R3du6e9tRRx017WoAADBzU4XoqvrRJG9L8pph1pYk75x2I929O8mHkjwmyeaq2tsCfkySXVNXCwAA68C0Jxb+VJITk1yXJN19cZJ7L7dCVR1VVZuH25uSPD7JRZmE6R8aFjs1ybvGlw0AAPMzbZ/or3b3jVWVJBlakle63N19kpw59Iu+Q5K3dvefVNVnk7y5qn45ySeSvO7ASgcAgPmYNkT/eVW9OMmmqnp8kp9M8sfLrdDdn07y8EXmX5pJ/2gAANiQpu3OsT3J1UkuSPJjSf40yUtmVRQAAKxnU7VEd/ctSV47/AAAwCFtqhBdVX+fRfpAd/d9V70iAABY56btE71twe27JPkPSe61+uUAAMD6N+1gK19e8LOru1+R5JQZ1wYAAOvStN05HrFg8g6ZtExP24oNAAAHlWmD8G8suH1zJkOAP23VqwGYkRN3nJ1du/fsN3/L5k1zqAaAjW7aq3M8btaFAMzSrt17snOHXmgArI5pu3P81+Xu7+6XrU45AACw/o25Oscjk7x7mH5Sko8muXgWRQEAwHo2bYg+Jskjuvv6JKmqlyY5q7t/eFaFAQDAejXtsN9HJ7lxwfSNwzwAADjkTNsS/XtJPlpV7ximn5rkzNmUBAAA69u0V+f4lap6T5LHDrOe092fmF1ZAACwfk3bnSNJ7prkuu5+ZZLLq+r4GdUEAADr2lQhuqp+MckLk7xomHWnJG+cVVEAALCeTdsS/e+SPDnJPydJd1+R5O6zKgoAANazaU8svLG7u6o6SarqbjOsCYCD0JbNm7J1+1mLzgfYaKYN0W+tqtck2VxVP5rkR5K8dnZlAXCwOWf7SfMuAWDVTHt1jl+vqscnuS7Jg5L8j+5+/0wrAwCAdWrFEF1VhyX5QHc/LongDADAIW/FEwu7+2tJbqmqe65BPQAAsO5N2yf6hiQXVNX7M1yhI0m6+2dmUhUAAKxj04botw8/AABwyFs2RFfVcd39D9195loVBAAA691KfaLfufdGVf3RjGsBAIANYaUQXQtu33eWhQAAwEaxUojuJW4DAMAha6UTCx9aVddl0iK9abidYbq7+x4zrQ4AANahZUN0dx+2VoUAAMBGseJgKwAAwG0J0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAw0rLDfgOsVyfuODu7du/Zb/6WzZtyzvaT5lARAIcSIRrYkHbt3pOdO07Zb/7W7WfNoRoADjW6cwAAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIMwvRVXVsVX2oqj5bVRdW1fOG+feqqvdX1cXD72+YVQ0AADALs2yJvjnJC7r7hCSPTvJTVXVCku1JPtjdD0jywWEaAAA2jJmF6O6+srs/Pty+PslFSbYkeUqSM4fFzkzy1FnVAAAAs7AmfaKramuShyc5N8nR3X3lcNcXkxy9FjUAAMBqmfmIhVV1RJI/SvL87r6uqr5+X3d3VfUS652e5PQkOe6442ZdJnCQ2LJ506KjFm7ZvGkO1XB7OJbAejbTEF1Vd8okQL+pu98+zP5SVd2nu6+sqvskuWqxdbv7jCRnJMm2bdsWDdoA+zpn+0nzLoFV4lgC69ksr85RSV6X5KLuftmCu96d5NTh9qlJ3jWrGgAAYBZm2RJ9YpJnJbmgqj45zHtxkh1J3lpVpyW5LMnTZlgDAACsupmF6O7+qyS1xN0nz2q7AAAwa0YsBACAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGGnmw34DHKgTd5ydXbv3LHqfoZ9ZDUu9xrZs3mTERGBZQjSwbu3avSc7d5wy7zI4iC31Gtu6/aw5VANsJLpzAADASEI0AACMJEQDAMBIQjQAAIzkxEIADnrLXYUD4EAI0QAc9FzpBVhtunMAAMBIQjQAAIwkRAMAwEj6RANz56QvADYaIRqYOyd9AbDR6M4BAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACPdcd4FAIeOE3ecnV279+w3f8vmTXOoBgAOnBANrJldu/dk545T5l0GANxuunMAAMBIQjQAAIwkRAMAwEhCNAAAjOTEQmDVuQoHG92WzZuydftZo5Y/Z/tJM6wIWG+EaGDVuQoHG93YQDwmcAMHB905AABgJCEaAABGEqIBAGAkIRoAAEZyYiEAB4XlrqjhyjDAahOiATgouMQcsJZm1p2jql5fVVdV1WcWzLtXVb2/qi4efn/DrLYPAACzMss+0W9I8oR95m1P8sHufkCSDw7TAACwocwsRHf3XyT5x31mPyXJmcPtM5M8dVbbBwCAWVnrq3Mc3d1XDre/mOTopRasqtOr6ryqOu/qq69em+oAAGAKc7vEXXd3kl7m/jO6e1t3bzvqqKPWsDIAAFjeWofoL1XVfZJk+H3VGm8fAABut7UO0e9Ocupw+9Qk71rj7QMAwO02y0vc/WGSv0nyoKq6vKpOS7IjyeOr6uIk3ztMAwDAhjKzwVa6+xlL3HXyrLYJAABrwYiFwIpO3HF2du3es9/8LZs3GSUOgEOSEA2saNfuPdm545T95m/dftYcqgGA+ZvbJe4AAGCjEqIBAGAkIRoAAEYSogEAYCQnFgLA7bRl86ZFT7R1BRs4eAnRAHA7LRWUXcEGDl66cwAAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEYsBIAZMRw4HLyEaACYEcOBw8FLdw4AABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGEmIBgCAkYRoAAAYSYgGAICRhGgAABhJiAYAgJGEaAAAGOmO8y4A2Li2bN6UrdvPWnQ+sLSl3jvLLX/O9pNmWBEwlhANHDD/1OHAjH3vjAncwNrQnQMAAEYSogEAYCQhGgAARhKiAQBgJCcWwkHsxB1nZ9fuPfvNX+pM/+WWBzaOpd7LiSt9wGoRouEgtmv3nuzcccp+85c603+p5YGNZbn3sit9wOrQnQMAAEYSogEAYCQhGgAARtInGuZknif+GK4bNpYDGSZ87GM54RDGEaJhTuZ54o9/lLCxrOZ7dqnHcsIhjDOX7hxV9YSq+nxVXVJV2+dRAwAAHKg1D9FVdViSVyf5/iQnJHlGVZ2w1nUAAMCBmkdL9KOSXNLdl3b3jUnenOQpc6gDAAAOyDxC9JYkX1gwffkwDwAANoR1e2JhVZ2e5PRh8oaq+vwcyjgyyTVz2C7LO2iOS/3qgd23WttYRQfNMTnIOC7rz7o+Jmv092I9WtfH5RA1r2PyLdMuOI8QvSvJsQumjxnm3UZ3n5HkjLUqajFVdV53b5tnDezPcVl/HJP1yXFZfxyT9clxWX82wjGZR3eOjyV5QFUdX1V3TvL0JO+eQx0AAHBA1rwlurtvrqqfTvJnSQ5L8vruvnCt6wAAgAM1lz7R3f2nSf50Htseaa7dSViS47L+OCbrk+Oy/jgm65Pjsv6s+2NS3T3vGgAAYEOZy4iFAACwkQnRU6iq51bV56rqwqr6tXnXw62q6gVV1VV15LxrOdRV1f8e3iefrqp3VNXmedd0qKqqJ1TV56vqkqraPu96SKrq2Kr6UFV9dvhf8rx518REVR1WVZ+oqj+Zdy1MVNXmqnrb8D/loqp6zLxrWowQvYKqelwmIyo+tLu/Ncmvz7kkBlV1bJLvS/IP866FJMn7k3xbdz8kyd8medGc6zkkVdVhSV6d5PuTnJDkGVV1wnyrIsnNSV7Q3SckeXSSn3Jc1o3nJblo3kVwG69M8t7ufnCSh2adHh8hemU/kWRHd381Sbr7qjnXw61enuTnkujYvw509/u6++Zh8iOZXAOetfeoJJd096XdfWOSN2fSEMAcdfeV3f3x4fb1mYQCo/XOWVUdk+SUJL8z71qYqKp7JvmuJK9Lku6+sbt3z7eqxQnRK3tgksdW1blV9edV9ch5F0RSVU9Jsqu7PzXvWljUjyR5z7yLOERtSfKFBdOXR1hbV6pqa5KHJzl3vpWQ5BWZNMbcMu9C+Lrjk1yd5HeHbja/U1V3m3dRi1m3w36vpar6QJJvWuSun89kH90rk6/fHpnkrVV133ZZk5lb4bi8OJOuHKyh5Y5Jd7+FNhHXAAADJklEQVRrWObnM/nq+k1rWRtsBFV1RJI/SvL87r5u3vUcyqrqiUmu6u7zq+p75l0PX3fHJI9I8tzuPreqXplke5JfmG9Z+xOik3T39y51X1X9RJK3D6H5o1V1SybjuV+9VvUdqpY6LlX17Zl8Uv1UVSWTbgMfr6pHdfcX17DEQ85y75UkqapnJ3likpN90JybXUmOXTB9zDCPOauqO2USoN/U3W+fdz3kxCRPrqofSHKXJPeoqjd29w/Pua5D3eVJLu/uvd/UvC2TEL3u6M6xsncmeVySVNUDk9w5yTVzregQ190XdPe9u3trd2/N5A33CAF6vqrqCZl8Lfrk7v7KvOs5hH0syQOq6viqunOSpyd595xrOuTV5BP/65Jc1N0vm3c9JN39ou4+Zvg/8vQkZwvQ8zf8L/9CVT1omHVyks/OsaQlaYle2euTvL6qPpPkxiSnamGDRb0qyeFJ3j98Q/CR7v7x+ZZ06Onum6vqp5P8WZLDkry+uy+cc1lMWj2fleSCqvrkMO/Fwwi+wG09N8mbhoaAS5M8Z871LMqIhQAAMJLuHAAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA2wwVTVz1fVhVX16ar6ZFV957xrAjjUuE40wAZSVY/JZFTIR3T3V6vqyEwGgQJgDWmJBthY7pPkmu7+apJ09zXdfUVV7RwCdapqW1V9eLj90qo6s6r+sqouq6ofrKpfq6oLquq9w1DUAIwkRANsLO9LcmxV/W1V/XZVffcU69wvyUlJnpzkjUk+1N3fnmRPklNmVyrAwUuIBthAuvuGJN+R5PQkVyd5S1U9e4XV3tPdNyW5IJOhwN87zL8gydbZVApwcNMnGmCD6e6vJflwkg9X1QVJTk1yc25tGLnLPqvs7fpxS1Xd1N09zL8l/g8AHBAt0QAbSFU9qKoesGDWw5JclmRnJi3USfLv17ougEONFgiAjeWIJL9VVZszaX2+JJOuHf8qyeuq6pcyaaUGYIbq1m/1AACAaejOAQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAw0v8P27GyrOJ7AW0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "xmin, xmax = -6,6\n",
    "\n",
    "fig2, ax2 = plt.subplots(figsize=(12, 6)) \n",
    "hist2 = ax2.hist(x_sum, bins=N_bins, range=(xmin, xmax), histtype='step')\n",
    "ax2.set(xlabel='Sum', ylabel='Frequency', title='Histogram of x_sum (the sum of the distributions)');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we overlay the plot with a unit Gaussian (i.e. not fitted to the data), to see if the resulting sums actually distribute themselves Gaussianly:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define your PDF / model \n",
    "def gauss_pdf(x, mu, sigma):\n",
    "    \"\"\"Normalized Gaussian\"\"\"\n",
    "    return 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-(x - mu) ** 2 / 2. / sigma ** 2)\n",
    "\n",
    "def gauss_extended(x, N, mu, sigma):\n",
    "    \"\"\"Non-normalized Gaussian\"\"\"\n",
    "    return N * gauss_pdf(x, mu, sigma)\n",
    "# Could also be written as:\n",
    "#gauss_extended = Extended(gauss_pdf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmclWX9//HXxT4IMgqoOEjjvqWS4pZSibu4oCGYGyCKS18st8LSpLJfZC64K6YiigoSKUZSJpqJW2iY+xoqi6LIiOzb9fvjPkMDzgyH4Zxzz5l5PR+Pecycc2/vswycz1zX/blDjBFJkiRJ0oZrknYASZIkSWooLLAkSZIkKUcssCRJkiQpRyywJEmSJClHLLAkSZIkKUcssCRJkiQpRyywJBWlEMLrIYTvpZ0jTSGE40MIH4cQFoQQvpV2ng0RQmgZQngjhNCplnWmhxAOKWSuQlrf1zOE8L0Qwow8ZSkPIcQQQrPM7cdCCP1ytO/uIYS3q9zO6eua738bQgi7hxCezdf+JRU/CyxJ9U51H7hCCP1DCM9U3o4x7hpjfGod+1njQ2IDdDXwfzHGNjHGf6cdZgMNAp6OMc4GCCGMDCFcmXKmQqv19cy8l7dLIRcxxiNjjPesa71sMsYY/xlj3DEXuap7n2Tzb8OGiDH+B6gIIRyTr2NIKm4WWJJUR/WgcPsG8HrKGXLlHODetEOkrCG9ntWqB78zuTIaODvtEJLqJwssSUWp6ihXCGGfEMLUEML8EMKnIYRrM6s9nflekZl2tX8IoUkI4bIQwochhDkhhFEhhHZV9nt6ZtncEMLlax1naAhhXAjhvhDCfKB/5tjPhRAqQgizQwg3hRBaVNlfDCGcF0J4N4TwVQjh1yGEbUMIz2byjq26/lqPsdqsmel0C4CmwCshhPer2fbbIYTPQwhbZW7vEUKYF0LYaR3P609DCDMzWd8OIRycuX+NkYK1p6dlnqdLQgj/CSEsDCHcGULYPDO17KsQwt9DCJvUcMwuwDbAC5nbg4BTgJ9kXrdHq6zeNXOML0MIY0IIrars5+gQwrTMa/FsCGH3Go4XQgjXZZ7T+SGEV0MI38wseyqEcGaVddcYOU3x9ax8L7+SeU76Vll2UWZ/s0MIA6rc3zKEcHUI4aOQ/F7cFkIoqSFb08y6n4cQPgB6rrV89fMSQtguhPCPzGvweQhhTE0ZK98nmffVJ8Dda793MvYOyRTReSGEuytf17Wf/yqvwXY1vU/Cmr+zLUMIw0MIszJfw0MILTPLKrPV9Pwdlcn0VeZ34uIqMZ4CDq7clyRVZYElqSG4Hrg+xrgxsC0wNnP/dzLfSzPTrp4D+me+DiL5UN8GuAkghLALcAvJh7ZOQDugbK1jHQeMA0pJ/oq9ErgA6ADsDxwMnLfWNocDewH7AT8BRgCnAlsB3wR+UMPjqjZrjHFpjLFNZp09Yozbrr1hjPFZ4HbgnsyH6vuAy2OMb9VwLEIIOwL/B+wdY2ybyT29pvWr8X3gUGAH4BjgMeBnQEeS/2/Or2G73YAPYowrMtlHkDy3V2Vet6pTsfoARwBbA7uTPD+E5Jylu0hGFdpnHvuEGj4AH0by3tiB5DXuA8xdj8eZxuv5nSrL28QYx2Rub8H/3qcDgZurFLLDMo+xK7BdZp1f1JDtLOBo4FtAN6B3zQ+fXwN/AzYBOgM3ZpFxU5IRukE17PMUkud120zmy2o5Ppnj1fY+qfRzktepK7AHsM9a+67t+bsTODvzu/BNYHKVY88ElgM5meooqWGxwJJUXz2cGYmoCCFUkBQ+NVkObBdC6BBjXBBjfL6WdU8Bro0xfhBjXABcCpwUkqlLvYFHY4zPxBiXkXwYjWtt/1yM8eEY46oY4+IY40sxxudjjCtijNNJPth/d61trooxzo8xvg68Bvwtc/wvSYqQmhoa1JY1G0NJPjy+CMwEbl7H+iuBlsAuIYTmMcbpMcavjabU4sYY46eZD5//BF6IMf47xrgE+BM1P85S4Kssj3FDjHFWjPEL4FGSD86QfHC/Pcb4QoxxZeZ8oaUkH67XthxoC+wEhBjjm5XnfmUprdezOsuBX8UYl8cY/wIsAHYMIQSS5+SCGOMXMcavgP8HnFTDfvoAw2OMH2ee29+u45jfALaMMS6JMT5Ty7oAq4ArMoXk4hrWuanKsX9DzUXq+jqF5PmZE2P8DPglcFqV5dU+f1WW7RJC2DjGOC/G+PJa+/6K5L0rSWuwwJJUX/WKMZZWfvH1UaGqBpL81futEMK/QghH17LulsCHVW5/CDQDNs8s+7hyQYxxEV8f2fi46o0Qwg4hhD+HED4JybTB/0cymlXVp1V+XlzN7TZUr7as6xRjXA6MJPnr+zUxxrWLxbXXfw/4MUlhNieE8GAIYctsjpVR18c5j6TgycYnVX5eVGWf3wAuWqso34rkOVxDjHEyyajlzSSPc0QIYeMsjw8pvZ41mFs58pdR+Zx0BFoDL1V5PiZl7q8pW9X39oc1rAfJqF0AXgxJx74z1pHxs0yRXZu1j70+77vaVPecV913Tc8fJCOyRwEfZqZE7r/WvtsCFTnKKakBscCSVPRijO/GGH8AbAb8DhgXQtiIr48+Acwi+TBeqQuwguRD8mySKU8AZKbWtV/7cGvdvhV4C9g+M0XxZyQfPnOhtqzrFEIoA64A7gauyeZ8kRjj/THGAzPHjSTPJ8BCkg/slbbIJkOW/gNsvdZITq3FYDU+Bn5TtSiPMbaOMT5Q3coxxhtijHsBu5AU55dkFuXzcW7Q67mePicp9nat8ny0qzIVcW2zSQrSqtmqFWP8JMZ4VoxxS5IpmbeE2jsHZvNarn3sWZmf13g9Qghrvx7r2nd1z/msGtZdc8cx/ivGeBzJvysP87+px5W/Wy2At2vYXFIjZoElqeiFEE4NIXSMMa7if39RXgV8lvm+TZXVHwAuCCFsHUJoQzLiNCbzV+xxwDEhaRDRgmQkZ13FUltgPrAgJA0kzs3V41pH1lplpoiNJDmPZCDJB+hfr2ObHUMIPTKF2BKSD+irMounAUeFEDbNfMj9cR0f09fEGGcA75GcH1PpU9Z83dblDuCcEMK+IbFRCKFnCOFrI2MhhL0z6zUn+QC/hDUf5wkhhNaZomFgXR5TDer8emZk/ZxkfhfuAK4LIWwGSVEQQji8hk3GAueHEDpnzkEaUtO+QwgnhhAq/xAxj6TIqXz+1vd1q/TDzLE3JTlvqvL8rVeAXUMIXTONL4autd26jvcAcFkIoWMIoQPJtN/71hUmhNAihHBKCKFdZiR4Pv97jJBMA54cY1yazYOT1LhYYElqCI4AXg9JJ7brgZMy50ctIjmfY0pmmtR+JI0Q7iXpMPhfkg/XgwEy59QMBh4kKUgWAHNIzuWpycXAySTnY9zB/z4Y5kKNWbNwPslf3i/PTA0cAAwIIXSvZZuWJI0RPieZircZyXlCZHK8QtL04m/k9nFCcu5a1XNj7iQ5/6UihPDwujaOMU4ladRwE8mH/vfINMCoxsYkr9U8kiljc4HfZ5ZdBywj+eB+D0kThVzZkNcTkuLinsxz0ieL9X9K8jw8n5m++ndqbspwB/BXktf4ZWB8LfvdG3gh8/s2AfhRjPGDOmasdD/J++oD4H3gSoAY4zvArzLZ3wXWPt9rXe+TK4GpJKOkr2YeW7bXVzsNmJ557s4hOZ+r0inAbVnuR1IjE9YxJV+SGq3MKEMFyfS//6adpyHLjJr9Gzh4PRtOSAUVkvb/t8cY1z4nS5IACyxJWkMI4RjgCZKpgdcA+wJ7rqtBhCRJEjhFUJLWdhzJSfCzgO1Jphs2mOIqhNAlJBdlre6rxsYGkiQpO45gSZIkSVKOOIIlSZIkSTmyIVePL5gOHTrE8vLytGNIkiRJaqReeumlz2OMNV2wfbWiKLDKy8uZOnVq2jEkSZIkNVIhhA+zWc8pgpIkSZKUIxZYkiRJkpQjFliSJEmSlCNFcQ6WJEmSVGjLly9nxowZLFmyJO0oKqBWrVrRuXNnmjdvXqftLbAkSZKkasyYMYO2bdtSXl5OCCHtOCqAGCNz585lxowZbL311nXah1MEJUmSpGosWbKE9u3bW1w1IiEE2rdvv0GjlhZYkiRJUg0srhqfDX3NLbAkSZIkKUcssCRJkiQpRyywJEmSpHpo+vTpfPOb31zjvqFDh3L11VfXut3UqVM5//zzAXjqqad49tlna1x30qRJ7LPPPuy000507dqVvn378tFHH214+PU0fPhwRo0alZN9XXLJJey6665ccskl3Hbbbav3O3LkSGbNmrV6vZNOOol33303J8esyi6CkiRJUgPSrVs3unXrBiQFVps2bfj2t7/9tfVee+01Bg8ezIQJE9h5550BmDBhAtOnT6dLly4Fy7tixQruuusuXn755Zzsb8SIEXzxxRc0bdp0jftHjhzJN7/5TbbccksAzj33XK666iruuOOOnBy3kgWWJEmStA4//jFMm5bbfXbtCsOH1337733ve+y77748+eSTVFRUcOedd9K9e3eeeuoprr76am666SZuu+02mjZtyn333ceNN95I9+7dV2//u9/9jp/97GeriyuAY489dvXPd9xxByNGjGDZsmVst9123HvvvbRu3Zr+/ftz9NFH07t3bwDatGnDggULmD17Nn379mX+/PmsWLGCW2+9lW9/+9sMHDiQqVOnEkLgjDPO4IILLljjcUyePJk999yTZs2a1fq4lixZwrnnnsvUqVNp1qwZ1157LQcddNAa+zr22GNZsGABe+21F5deeilvvvkmbdq0oby8nKlTp3LKKadQUlLCc889R/fu3enfvz8rVqxYfexccIqgJEmSVKRWrFjBiy++yPDhw/nlL3+5xrLy8nLOOeccLrjgAqZNm7ZGcQXw+uuvs+eee9a47xNOOIF//etfvPLKK+y8887ceeedtWa5//77Ofzww5k2bRqvvPIKXbt2Zdq0acycOZPXXnuNV199lQEDBnxtuylTprDXXnut83HdfPPNhBB49dVXeeCBB+jXr9/X2qlPmDCBkpISpk2bRt++fVff37t3b7p168bo0aOZNm0aJSUlNGnShO22245XXnml1se1vhzBkiRJktZhQ0aa6qqmduFV7z/hhBMA2GuvvZg+fXqdjzV37lwOPvhgFi1axKBBg7j44ot57bXXuOyyy6ioqGDBggUcfvjhte5j77335owzzmD58uX06tWLrl27ss022/DBBx8wePBgevbsyWGHHfa17WbPnr3GKFpNj+uZZ55h8ODBAOy000584xvf4J133mH33Xev8+PebLPNmDVr1tcKvA3hCJYkSZJUD7Vv35558+atcd8XX3xBhw4dVt9u2bIlAE2bNmXFihXrtf9dd9119XlP7du3Z9q0aQwaNIgFCxYA0L9/f2666SZeffVVrrjiitWjRc2aNWPVqlUArFq1imXLlgHwne98h6effpqysjL69+/PqFGj2GSTTXjllVf43ve+x2233caZZ575tRwlJSVfG4nakMe1PpYsWUJJSUlO92mBJUmSJNVDbdq0oVOnTkyePBlIiqtJkyZx4IEHZr2Ptm3b8tVXX1W77Cc/+Qm/+c1vePPNN1fft2jRotU/f/XVV3Tq1Inly5czevTo1feXl5fz0ksvAcmUvOXLlwPw4Ycfsvnmm3PWWWdx5pln8vLLL/P555+zatUqvv/973PllVdW28hi55135r333lvnY+nevfvqHO+88w4fffQRO+64YxbPQqK65+Kdd975WqfGDeUUQUlSvXLAsMnMrFhc7bKy0hKmDOlR4ESSlJ5Ro0bxwx/+kAsvvBCAK664gm233Tbr7Y855hh69+7NI4888rUmF7vtthvXX389p59+OvPnz6dDhw506dJl9TlPv/71r9l3333p2LEj++677+ri5KyzzuK4445jjz324IgjjmCjjTYCko6Fv//972nevDlt2rRh1KhRzJw5kwEDBqwe8frtb3/7tYxHHnkkp5122jofy3nnnce5557LbrvtRrNmzRg5cuTqka5s9O/fn3POOWd1k4v58+dTUlLCFltskfU+shFijDndYT5069YtTp06Ne0YkqQCKB8ykenDeq73MknKtTfffPNr5wYpP44//niuuuoqtt9++4Id87rrrmPjjTdm4MCBX1tW3WsfQngpxthtXft1iqAkSZKkVA0bNozZs2cX9JilpaX069cv5/t1iqAkSZJUgxhjjd38lDs77rjjep1PlQvVtYyH5DXfEI5gSZIkSdVo1aoVc+fO3eAP3CoeMUbmzp1Lq1at6rwPR7AkSZKkanTu3JkZM2bw2WefpR1FBdSqVSs6d+5c5+0tsCRJkqRqNG/enK233jrtGCoyThGUJEmSpByxwJIkSZKkHLHAkiRJkqQcscCSJEmSpByxwJIkSZKkHLGLoCSpwTtg2GRmViyucXlZaQlThvQoYCJJUkNlgSVJavBmVixm+rCeNS4vHzKxgGkkSQ2ZUwQlSZIkKUcssCRJkiQpRyywJEmSJClHLLAkSZIkKUcssCRJkiQpRyywJEmSJClHLLAkSZIkKUcssCRJkiQpRyywJEmSJClHLLAkSZIkKUcssCRJkiQpRyywJEmSJClHLLAkSZIkKUcssCRJkiQpR5rlc+chhOnAV8BKYEWMsVsIYVNgDFAOTAf6xBjn5TOHJEmSJBVCIUawDooxdo0xdsvcHgI8EWPcHngic1uSJEmSil4aUwSPA+7J/HwP0CuFDJIkSZKUc/kusCLwtxDCSyGEQZn7No8xzs78/AmweXUbhhAGhRCmhhCmfvbZZ3mOKUmSJEkbLq/nYAEHxhhnhhA2Ax4PIbxVdWGMMYYQYnUbxhhHACMAunXrVu06kiRJklSf5HUEK8Y4M/N9DvAnYB/g0xBCJ4DM9zn5zCBJkiRJhZK3AiuEsFEIoW3lz8BhwGvABKBfZrV+wCP5yiBJkiRJhZTPKYKbA38KIVQe5/4Y46QQwr+AsSGEgcCHQJ88ZpAkSZKkgslbgRVj/ADYo5r75wIH5+u4kiRJkpSWfDe5kCQpZ8pKSygfMrHGZVOG9ChwIkmS1mSBJUkqGrUVUDUVXpIkFVIaFxqWJEmSpAbJAkuSJEmScsQCS5IkSZJyxAJLkiRJknLEAkuS1GAtWgTvvAPL5rTlvfeS25Ik5ZNdBCVJDUJZaQnf+OlEls7YhIVvlrFkentWzGuTWfodtr8bQoCdd4ZDD4WTToJ9903ukyQpVyywJElFL0b40TY9+NWv4KPXYaON4IiDYJ99YOutoXVrWLAAPvgAnn8ebr0Vrr8e9t4bhg5NtpckKRcssCRJRe3992HQIJg8ORmduvtuOPHEpMiqyfz5MHo0XHMN9OwJrbbemxmDoXPnwuWWJDVMnoMlSSpad90Fu+0GU6fCbbfBq69C//61F1cAG28M554Lb7wBw4fD0hmbsvvu8NBDBYktSWrALLAkSUVnxQo4/3wYOBAOOCAplM4+G5o2Xb/9tGgBP/oRdOr/T7bbDvr0gV/8wimDkqS6s8CSJBWVJUugVy+48Ua48EJ47DEoK9uwfTbfdBFTpsAZZ8Cvfw39+sGyZbnJK0lqXDwHS5KUcwcMm8zMisXVLisrLWHKkB512u+SJXD88TBpUjIl8OyzNyTlmpo3hz/8AcrLk1GsxYvhgQegmf9TSpLWg/9tSJJybmbFYqYP61ntsvIhE+u0z+XLk+Lqr3+FO+9MRptyLQS4/PLkHK6LLoKSEhg5Epo430OSlCULLElSvRdj0pRi0iQYMSI/xVVVF16YXJT48suhfXu47rr8Hk+S1HBYYEmS6r3f/jYZtbr8cjjrrMIc8+c/h88+S7oM7rornHlmYY4rSSpuTnqQJNVrf/5zUuyccgr88peFO24IyXWyDjsMzjsPnn66cMeWJBUvCyxJUr01fTqcfjp861tJA4oQCnv8Zs1gzBjYemvo2xfmzCns8SVJxccpgpKkgiorLam10UVZaQkAS5fCiSfCqlXJBYBbtSpUwjWVlsK4cbD33kn79okT/9f0orZuibBhHRMlScXJAkuSVFDZFhxDh8LUqTB+PGy7bX4zrctuu8G118IPf5h8v/ji5P7auiVC3TsmSpKKl1MEJUn1znPPwVVXwcCBSWv2+uDcc5MsP/85vP562mkkSfWVBZYkqV5ZtCiZite5czJaVF+EALffDhtvDAMGwIoVaSeSJNVHFliSpHpl6FB49124++6kmKlPOnaEG2+Ef/3La2NJkqpngSVJqjdeey0pXM48E3rU094QfftCr17JNbmWz2uddhxJUj1jgSVJqhdiTK431a4dDBuWdpqahQA33wzNm8O8J3ZJO44kqZ6xwJIk1Qv33gv//Cf87nfQvn3aaWq35ZbJVMbF72/Oo4+mnUaSVJ9YYEmSUjdvXtL6fP/9kwYSxeD886F5+684/3xYXPOlsCRJjYwFliQpdVdeCXPnwi23/O8ivvVd8+aw6WGvMX06DB+edhpJUn1RJP+NSZIaqv/+F266KRm56to17TTrp1WXLzj22OScsc8/TzuNJKk+sMCSJKXqssugaVP45S/TTlI3w4bBggXJKJwkSRZYkqTUvPQS3H8/XHABlJWlnaZudt45aSt/yy3w/vtpp5Ekpa1Z2gEkScXpgGGTmVlRfXeHstKSdW4fI/zkJ9ChQ/I9TWWlJZQPmVjjsilDar8o19ChcN998POfw4MP5iGgJKloWGBJkupkZsVipg/rWeftH38cJk+GG25Irn2VptoKqJoKr6o6dUq6IP7qV3DJJbDXXrlMJ0kqJk4RlCQVXIxwxRWw1VZw9tlpp8mNiy6CTTYp3nPJJEm5YYElSSq4xx+H55+Hn/0MWrRIO01ubLwxXHghPPoovPxy2mkkSWmxwJIkFVSMySjPVlsVz0WFszV4sKNYktTYWWBJkgrq73+HZ5+FSy+Fli3TTpNb7dolHREnTIB//zvtNJKkNFhgSZIKpnL0qnNnOOOMtNPkx/nnQ2mpo1iS1FhZYEmSCuapp2DKFBgypOGNXlVq1w5+/GN45BFYNqdt2nEkSQVmgSVJKpirroLNNoOBA9NOkl/nnw8bbQTzX9wm7SiSpAKzwJIkFcR//gOTJsGPfgStWqWdJr822QTOOgsWvrElH36YdhpJUiFZYEmSCuL3v09Gdc49N+0khXHBBcn34cPTzSFJKiwLLElS3n34ITzwAAwalIzuNAZdusBGu8zijjtg3ry000iSCsUCS5KUd8OHQwhJ84fGZON9PmDhQrj11rSTSJIKpVnaASRJDdu8eXDHHXDSScmoTmPSYrOvOOIIuP56uPDCNc89O2DYZGZWLK52u7LSEqYM6VGglJKkXLLAkiTl1Z13wsKFcNFFaSdJx0UXwaGHwpgx0K/f/+6fWbGY6cN6VrtN+ZCJBUonSco1pwhKkvJm5Uq4+Wb4znega9e006Tj4INh553hxhuTCy1Lkho2CyxJUt5MnAjTp8PgwWknSU8IyeN/6SV4/vm000iS8s0CS5KUNzfcAJ07Q69eaSdJ12mnQbt2yfMhSWrYLLAkSXnxxhvwxBNw3nnQrJGf8dumDZxxBowbB7NmpZ1GkpRPFliSpLy46SZo2RLOOivtJPXDD3+YnJN2++1pJ5Ek5ZMFliQp5yoq4J574Ac/gA4d0k5TP2y7LRx1FNx2GyxdmnYaSVK+WGBJknLu7rth0aLG3dyiOuefD3PmwEMPpZ1EkpQveS+wQghNQwj/DiH8OXN76xDCCyGE90IIY0IILfKdQZJUOKtWwS23wAEHwJ57pp2mfjnkENhxx2T6pCSpYSrECNaPgDer3P4dcF2McTtgHjCwABkkSQXy5JPw3ntw7rlpJ6l/mjSBc86BF16AZXPaph1HkpQHeS2wQgidgZ7AHzK3A9ADGJdZ5R6gkTfvlaSGZcQI2GQT+P73005SP512WtL8Y8ErXdKOIknKg3yPYA0HfgKsytxuD1TEGFdkbs8AyqrbMIQwKIQwNYQw9bPPPstzTElSLsyZA3/6E/TrB61apZ2mfmrfHnr3hgWvlbFoUdppJEm5lrcCK4RwNDAnxvhSXbaPMY6IMXaLMXbr2LFjjtNJkvLhnntg+XJbs6/LoEEQlzVn7Ni0k0iSci2fI1gHAMeGEKYDD5JMDbweKA0hVF5ysjMwM48ZJEkFEmMyPfDAA2GXXdJOU7917w7NNl3AiBFpJ5Ek5VreCqwY46Uxxs4xxnLgJGByjPEU4Emgd2a1fsAj+cogSSqcp55KmlsMGpR2kvovBGjb9SOeew5efTXtNJKkXErjOlg/BS4MIbxHck7WnSlkkCTlWGVzi969172uYKNvzqBFCxzFkqQGpiAFVozxqRjj0ZmfP4gx7hNj3C7GeGKM0evZS1KR++wzGD8eTj8dSkrSTlMcmpYsp3dvuPdebHYhSQ1IGiNYkqQGZtQoWLbM5hbr6+yz4csv4aGH0k4iScoVCyxJ0gaJEUaOhH33hV13TTtNceneHbbbLnn+JEkNgwWWJGmDvPQSvPYaDBiQdpLiEwL07580CPngg7TTSJJywQJLkrRBRo6Eli2hb9+0kxSn005LCq1Ro9JOIknKBQssSVKdLV0K998Pxx8PpaVppylOXbrAwQcnF2letSrtNJKkDWWBJUmqs0cfhXnznB64oQYMgOnT4emn004iSdpQFliSpDq7+24oK0tGYFR3vXrBxhsnz6ckqbg1SzuAJKk4rVjQkkmT4Kc/haZN006TP2WlJZQPmVjjslxs27p1cg7b6NFw0011zypJSp8FliSpTha+XsaqVdCvX9pJ8mvKkB4F2bZ/f7jjDhg3rs6HkyTVA04RlCSttxhh4aud2X9/2HHHtNM0DPvvDzvs4DWxJKnYWWBJktbb1KmwfG5bm1vkUOU1sZ5+GpbPa512HElSHVlgSZLW2913Q2i2kj590k7SsFReE2vha53TjiJJqiMLLEnSelm6FB58EEq2/4R27dJO07B07gyHHgoLXivzmliSVKQssCRJ6+Wxx5JrX7XZdWbaURqkU0+FlfNb8+yzaSeRJNWFBZYkab2MHg0dO0KrrT9PO0qD1KtXMv0PvyjGAAAgAElEQVRy9Oi0k0iS6sICS5KUtS+/hEcfhZNOgtAkph2nQWrbNpl+OXYsLFuWdhpJ0vqywJIkZe2Pf0zOwTrllLSTNGwb7TqLL76ASZPSTiJJWl8WWJKkrI0eDdtuC/vsk3aShq2k/DM6dMBpgpJUhCywJElZmTkTnnwyGb0KIe00DVtoGunbFyZMgPnz004jSVofFliSpKw8+CDE6PTAQjnlFFiyBMaPTzuJJGl9WGBJkrJy332w996www5pJ2kc9tsPttnGaYKSVGwssCRJ6/TGGzBtmqNXhRQCnHwyTJ4Ms2alnUaSlC0LLEnSOo0eDU2bJu3ZVTinnAKrViXTMyVJxcECS5JUq1WrkgLrkENg883TTtO47LQT7LWX0wQlqZhYYEmSavXss/Dhh04PTMspp8DLL8Nbb6WdRJKUDQssSVKtRo+GkhLo1SvtJI3TSSdBkyaOYklSsWiWdgBJUv10wLDJzJi7hBkjD6HV1p+x22+mrbG8rLQkpWSNS6dOcPDBSYH1q195DTJJqu8ssCRJ1ZpZsZjbevTkqKvhgWFlHHtsWdqRGq2TT4YBA+Bf/4J99kk7jSSpNk4RlCTVaOxYaNcODj887SSNW69e0KKF3QQlqRhYYEmSqhVXNOFPf0o+3LdsmXaaxq20FI44Iil4V61KO40kqTYWWJKkai2e3oEvv4Q+fdJOIoC+fWHmTJgyJe0kkqTaWGBJkqq16K1ObLJJcv0rpe/YY5Nujk4TlKT6zQJLkvQ1S5bAonc35/jjk3N/lL42baBnTxg3DlasSDuNJKkmFliSpK/5618hLmtO375pJ1FVJ50Ec+bAU0+lnUSSVBMLLEnS14wZA01KlnHQQWknUVVHHZWMZI0Zk3YSSVJNLLAkSWtYvBgmTIDWO3xC8+Zpp1FVJSVw3HHwxz/CsmVpp5EkVccCS5K0hr/8BRYuhNY7zUo7iqrRty/Mmwd//3vaSSRJ1bHAkiStYexY6NgRWnX5Iu0oqsZhhyXXxbKboCTVT83SDiBJyq8Dhk1mZsXiapeVlZYwZUiP1bcXLoQ//xn69YO/NImFiqj10LIlnHACPPRQ0u2xVau0E0mSqrLAkqQGbmbFYqYP61ntsvIhE9e4PXEiLFqUXFz4L5MKkU510bcv3HUXPPYYHH982mkkSVVlNUUwhLBbvoNIktI3dixssQV07552EtWmRw/o0MFugpJUH2V7DtYtIYQXQwjnhRDa5TWRJCkVX32VjGD17g1Nm6adRrVp1ix5nR59NJnWKUmqP7IqsGKM3YFTgK2Al0II94cQDs1rMklSQf35z8k5PX36pJ1E2TjppGQ656OPpp1EklRV1l0EY4zvApcBPwW+C9wQQngrhHBCvsJJkgpn7FjYcks44IC0kygbBx4InTo5TVCS6pusmlyEEHYHBgA9gceBY2KML4cQtgSeA8bnL6IkKV/KSksoHzKRVUub8fGjh9C260ds87M3Vi9T/dW0aTLaeOut8OWX0M4J/JJUL2TbRfBG4A/Az2KMq3v9xhhnhRAuy0sySVLeVbZov+8+OG0l/PXGrdl//61TTqVs9e0L118PDz+ctNaXJKUv2ymCPYH7K4urEEKTEEJrgBjjvfkKJ0kqjDFjYKutYN99006i9bHfftClSzK9U5JUP2RbYP0dqDpXpHXmPklSkauogL/+NZlu1iTrM3NVH4SQvG5/+xt88UXaaSRJkH2B1SrGuKDyRubn1vmJJEkqpEcegeXL7R5YrPr2hRUrkmmCkqT0ZVtgLQwh7Fl5I4SwF7C4lvUlSUVizBgoL4e99047iepir71gm23sJihJ9UW2BdaPgYdCCP8MITwDjAH+L3+xJEmF8MUX8PjjyehVCGmnUV1UThN84gn4/PO000iSsuoiGGP8VwhhJ2DHzF1vxxiX5y+WJKkQ/vSnZHqZ0wOLywHDJjOz4n8TSZZ9ujErV3Znx1P/w07f+3x1d8h1bbe2stKSGreVJGUn2zbtAHsD5Zlt9gwhEGMclZdUkqSCGDsWtt0W9txz3euq/phZsZjpw3quvh0j7DgFtlq+O+9XTMx6u7WVD6l5W0lSdrKaIhhCuBe4GjiQpNDaG+iWx1ySpDz7/PNkWpnTA4tfCEmzi6eegpULW6QdR5IatWxHsLoBu8QYY7Y7DiG0Ap4GWmaOMy7GeEUIYWvgQaA98BJwWoxx2frFliRtqPHjYeVKpwc2FH36wJVXwqK3O6UdRZIatWybXLwGbLGe+14K9Igx7gF0BY4IIewH/A64Lsa4HTAPGLie+5Uk5cDYsbDDDrDHHmknUS5885uw886w8C0LLElKU7YFVgfgjRDCX0MIEyq/atsgJiqvndU88xWBHsC4zP33AL3qkFuStAE+/RSefNLpgQ1J5TTBpR9vyqxZaaeRpMYr2ymCQ+uy8xBCU5JpgNsBNwPvAxUxxhWZVWYAZTVsOwgYBNClS5e6HF6SVIPx42HVquQDuRqOPn1g6NDAH/8IgwennUaSGqesRrBijP8ApgPNMz//C3g5i+1Wxhi7Ap2BfYCdsg0WYxwRY+wWY+zWsWPHbDeTJGVh7NhkOtmuu6adRLm0887QvON8LzosSSnKtovgWSTT+m7P3FUGPJztQWKMFcCTwP5AaQihcuSsMzAz67SSpA02ezb84x/J6JXTAxue1jvNZsoU+PjjtJNIUuOU7TlYPwQOAOYDxBjfBTarbYMQQscQQmnm5xLgUOBNkkKrd2a1fsAj6x9bklRXf/xjct2kE09MO4nyYaOdkhOwxo1bx4qSpLzItsBaWrWVemYEal0t2zsBT4YQ/kMypfDxGOOfgZ8CF4YQ3iNp1X7n+seWJNXVmDFJx7lddkk7ifKh+aaL+Na3cJqgJKUk2yYX/wgh/AwoCSEcCpwHPFrbBjHG/wDfqub+D0jOx5IkFdjMmfDMM/DrX6edRPnUpw9ceilMnw7l5WmnkaTGJdsRrCHAZ8CrwNnAX4DL8hVKkpQfldPGvLhww1b5+j70ULo5JKkxyraL4KoY4x0xxhNjjL0zP69riqAkqZ4ZMwa6dk0uMKyGa5ttoFu3pFukJKmwsu0i+N8Qwgdrf+U7nCQpdz76CJ57ztGrxqJvX5g6Fd5/P+0kktS4ZDtFsBuwd+arO3ADcF++QkmScq9yupgFVuNQ2SXSaYKSVFjZThGcW+VrZoxxONAzz9kkSTk0dizstRdsu23aSVQI3/gG7Lef3QQlqdCynSK4Z5WvbiGEc8i+A6EkKWXTp8OLLzp61dj07QvTpsE776SdRJIaj2ynCF5T5eu3wF6A/01LUpGobHZggdW49O6dfLfZhSQVTlajUDHGg/IdRJJUNwcMm8zMisU1Li8rLWHsWNhnH6+J1Nh07gwHHphME7zMi6tIUkFkVWCFEC6sbXmM8drcxJEkra+ZFYuZPqzm02Lfew+2vxSuuaaAoVRv9OkD558Pb7yRdhJJahzWp4vguUBZ5uscYE+gbeZLklRPVXaRq5wupsald28IwWmCklQo2Taq6AzsGWP8CiCEMBSYGGM8NV/BJEm5MXYs7L8/dOmSdhKloVMn+O53k/dBPCbtNJLU8GU7grU5sKzK7WWZ+yRJ9dg77yRd5Pr2TTuJ0tSnD7z5Jiz/3EknkpRv2RZYo4AXQwhDM6NXLwD35C2VJCknxo5Npoc5PbBx+/73oUkTWPhmp7SjSFKDl20Xwd+EEB4DumfuGhBj/Hf+YkmScmHMmKSLXFlZ2kmUps02g4MOgn+83IkYk6JbkpQf2Y5gAbQG5scYrwdmhBC2zlMmSVIOvPEGvPaa175Som9fWDGvDdOmpZ1Ekhq2rAqsEMIVwE+BSzN3NQfuy1coSdKGGzMmmRbm9EABHH88EFbZTVCS8izbEazjgWOBhQAxxlnYnl2S6q0YkwLru9+FLbZIO43qgw4doFX554wZk7w/JEn5kW2BtSzGGIEIEELYKH+RJEkb6j//gbfftnug1rTRTrP573/hpZfSTiJJDVe2BdbYEMLtQGkI4Szg78Ad+YslSdoQY8ZA06ZJ9zipUskOn9C8efL+kCTlR7ZdBK8OIRwKzAd2BH4RY3w8r8kkSXVSOT3w4IOTaWEqPmWlJZQPmVjr8rpo2moFhx2WtO+/6iq7CUpSPqyzwAohNAX+HmM8CLCokqR6bupU+OAD+PnP006iupoypEfe9t2nD0ycCC+8APvtl7fDSFKjtc4pgjHGlcCqEEK7AuSRJG2gMWOgefNM1zhpLccdBy1aOE1QkvIlqymCwALg1RDC42Q6CQLEGM/PSypJUp2sWpVM/zr8cNhkk7TTqD5q1w6OPBIeegiuuSZp5S9Jyp1s/1kdD1wOPA28VOVLklSPPP88fPyx3QNVuz59YOZMePbZtJNIUsNT6whWCKFLjPGjGOM9hQokSaq7MWOgZUs49ti0k6g+O+YYaNUqGe088MC000hSw7KuEayHK38IIfwxz1kkSRtg5cpk2tdRR8HGG6edRvVZ27bJ++Shh5L3jSQpd9ZVYFVt4LpNPoNIkjbMM8/A7NlOD1R2+vaFTz5J3jeSpNxZV4EVa/hZklTPjBkDrVvD0UennUTFoGfP5P1iN0FJyq11FVh7hBDmhxC+AnbP/Dw/hPBVCGF+IQJKktZtxQoYNy4prjbaKO00KgYbbZS8X8aNS94/kqTcqLXAijE2jTFuHGNsG2Nslvm58rYz/CWpnnjqKfjsM6cHav307Zu8b/7xj7STSFLD4dUvJKkBePBBaNMmub6RlK0jj0zeN04TlKTcscCSpCIXVwbGj4fjjoOSkrTTqJiUlCQt/cePh+XL004jSQ2DBZYkFbkl0zswb57TA1U3ffrA3LkweXLaSSSpYbDAkqQit/CtLWnXDg47LO0kKkZHHJFcN23s2LSTSFLDYIElSUVsyRJY9M7mHH88tGyZdhoVo5YtoVevZJpgXBnWvYEkqVYWWJJUxCZNgrisudMDtUH69IGKClg8vWPaUSSp6FlgSVIRGz0amrReyiGHpJ1ExezQQ6G0FBa91SntKJJU9CywJKlIzZ8Pjz4KG+00m2bN0k6jYtaiBRx/fDLddMmStNNIUnGzwJKkIvWnP8HSpbDRLjPTjqIGoG/fZLrp3/6WdhJJKm4WWJJUpEaPhq23hhZbVqQdRQ1Ajx7QpGSZFx2WpA1kgSVJReiTT+CJJ+DkkyHY+E050Lw5tN7hEyZMgMWL004jScXLAkuSitDYsbBqFZxyStpJ1JC03mkWCxbAX/6SdhJJKl4WWJJUhEaPhq5dYeed006ihqRVl7lssUXy/pIk1Y0FliQVmffegxdfTKYHSrkUmiTvq4kT4Ysv0k4jScXJxr6SVE8cMGwyMyuqP/mlrLSEKUN6AHD//cl5Vz/4QSHTqbE49VS49loYNw4GDUo7jSQVHwssSaonZlYsZvqwntUuKx8yEYAYkwLrO9+Bzp0LmU6NRdeusMsucN99FliSVBdOEZSkIvLvf8Pbb9vcQvkTQvL++uc/Yfr0tNNIUvGxwJKkIjJ6dNJO+/vfTzuJGrLK8/vuvz/dHJJUjCywJKlIrFwJDz4IRx4Jm26adho1ZOXl0L073HtvMi1VkpQ9CyxJKhJPPw2zZjk9UIVx6qnw1lvJtFRJUvYssCSpSIweDW3awNFHp51EjcGJJ0KLFkmzC0lS9iywJKkIxBVNGDcOTjgBWrdOO40ag002gZ494YEHYMWKtNNIUvGwwJKkIrDovc348kunB6qwTj0VPvkEJk9OO4kkFQ8LLEkqAgtf68yWW8LBB6edRI3JUUdBaWkyPVWSlJ28FVghhK1CCE+GEN4IIbweQvhR5v5NQwiPhxDezXzfJF8ZJKkhmDMHFn/QkVNPhaZN006jxqRVq+RcrPHjYeHCtNNIUnHI5wjWCuCiGOMuwH7AD0MIuwBDgCdijNsDT2RuS5JqcP/9QGxCv35pJ1FjdOqpsGABTJiQdhJJKg55K7BijLNjjC9nfv4KeBMoA44D7smsdg/QK18ZJKkhGDUKWmxRwS67pJ1EjdGBB8JWW9lNUJKy1awQBwkhlAPfAl4ANo8xzs4s+gTYvIZtBgGDALp06ZL/kJJUD736anIdoq2PmUP5kCnVrlNWWlLgVKrPykpLKB8yscZl66tJk6S5yu9/D59+CptX+7+2JKlS3gusEEIb4I/Aj2OM80MIq5fFGGMIodprxMcYRwAjALp16+Z15CU1SvfcA82awYt37UCHDjukHUdFYMqQHjnf52mnwbBhyXTVCy7I+e4lqUHJaxfBEEJzkuJqdIxxfObuT0MInTLLOwFz8plBkorVihVJ97aePaFDh7TTqDHbZRfYZx+4+26I/slTkmqVzy6CAbgTeDPGeG2VRROAylO1+wGP5CuDJBWzxx9PrkFkcwvVBwMGJFNWX3457SSSVL/lcwTrAOA0oEcIYVrm6yhgGHBoCOFd4JDMbUnSWkaNgk03Ta5FJKWtb19o2RJGjkw7iSTVb3k7ByvG+AwQaljspTIlqRZffgkPPwwDByYfaqW0bbIJHH98ch7W1Vf7vpSkmuT1HCxJUt2MHQtLljg9UPXLgAHwxRdeE0uSamOBJUn10KhRsNNO0K1b2kmk/zn4YOjcOWl2IUmqngWWJNUz778PzzyTjF6FmiZaSylo2hROPx3++leYNSvtNJJUP1lgSVI9c9ddycVdTzst7STS1/XvD6tWwb33pp1EkuonCyxJqkdWrEi6tB15JJSVpZ1G+rrtt4cDD/SaWJJUEwssSapHJk1Kpl6deWbaSaSa9e8Pb78NL7yQdhJJqn8ssCSpHrnzTth8c+jZM+0kUs369IHWrW12IUnVydt1sCRJazpg2GRmViyucXnHpu149FG46CJo3ryAwaT11LYt9O4NDz4I110Hh95Q+3u7rLSEKUN6FDChJKXHAkuSCmRmxWKmD6t5aOqqq2DqSjjjjAKGkuqof//kcgLjx6/7vV0+ZGLhgklSypwiKEn1QIzJ9MDu3WHHHdNOI63bd78L220HI0aknUSS6hcLLEmqB555Bt55BwYOTDuJlJ0mTeCss+Cf/4Tln7dJO44k1RsWWJJUD/zhD/87r0UqFv37J+cLfvXKVmlHkaR6wwJLklL25Zfw0ENw8smw0UZpp5Gyt9lm0KsXLHytM0uWpJ1GkuoHm1xIUsoeeAAWL3Z6oOq3mrpgLm7SnlVL9mP8+OSPBJLU2FlgSVKKYoTbb4fdd4du3dJOI9Wspk6Bq1bBDlOTZhcWWJLkFEFJStULL8C0aXDuuRBC2mmk9VfZ7OIf/4C33047jSSlzwJLklJ0663Qpg2cckraSaS6698fmjWDO+5IO4kkpc8CS5JSMncujBkDp52WdBCUitXmmyfNLkaOhKVL004jSemywJKklFR+GD333LSTSBtu0KDkjwbjx6edRJLSZZMLScqhmjqtAZSVlqz+edUquO02OOAA2G23QqWT8ufgg2GbbZL39Q9+kHYaSUqPBZYk5VBNndbW9ve/w3vvwdCh+c8kFUKTJnDOOfCTn8Crr/qHA0mNl1MEJSkFt94KHTpA795pJ5FyZ+BAaNUKbrop7SSSlB4LLEkqsBkzYMIEOOMMaNky7TRS7my6adIR8777YN68tNNIUjossCSpwO64I7nA8Nlnp51Eyr0f/hAWLUqauEhSY2SBJUkFtHx5UmAdcUTSEEBqaL71raR5y803J81cJKmxscCSpAIaNw5mz07+yi81VIMHw/vvw6RJaSeRpMKzwJKkArr+eth+ezjyyLSTSPlz/PHQqZPNLiQ1ThZYklQgL7yQfA0enLS0lhqqFi2Scwwfeyy5HIEkNSb+Fy9JBXL99bDxxtC/f9pJpPwbNAiaNYNbbkk7iSQVlgWWJBXAzJnw0ENJa/a2bdNOI+Vfp07Jdd7uugtWLWuadhxJKhgLLEkqgFtvhZUrk+mBUmPxox/Bl1/Cgle3SjuKJBWMBZYk5dmSJXD77XDMMbZmV+Oy336w//7w1dRyVq5MO40kFYYFliTl2f33w+efJ3/Nlxqbiy6CFRUb8fDDaSeRpMKwwJKkPIoRrr0WdtsNDjoo7TRS4fXqBc1KF3LttWknkaTCsMCSpDx67DF4/XW45BIIIe00UuE1bQptu03n2Wfh+efTTiNJ+WeBJUl5dNVVsNVWcNJJaSeR0tNmt48pLYVrrkk7iSTlnwWWJOXJiy/CP/4BP/4xNG+edhopPU1arOScc2D8ePjvf9NOI0n5ZYElSXny+99Du3Zw1llpJ5HSN3hwMl1w+PC0k0hSfjVLO4AkFZsDhk1mZsXiapeVlZYA8N578Mc/wpAhXlhYAthySzj5ZPjDH+Dyy6FDh7QTSVJ+WGBJ0nqaWbGY6cN61rrONdck0wLPP79AoaQi8NOfwqhRcMMN8KtfpZ1GkvLDKYKSlGNz5sDIkXD66bDFFmmnkeqPnXeGE06AG2+E+fPTTiNJ+WGBJUk5Nnw4LF0KF1+cdhKp/rn0UqiogFtvTTuJJOWHBZYk5dAXX8BNN8GJJ8KOO6adRqp/9toLDj88uQD34upPZZSkomaBJUk5dP318NVXcNllaSeR6q+f/SyZSnvXXWknkaTcs8mFJK2lti6B8L9OgWv78sukwOrVC3bbLV/ppOJTVlpC+ZCJq2/HCC3L9udHP2vFVdOfIjSNNW43ZUiPQsWUpJywwJKktWTTJbA6N92UFFmXX56HUFIRq65I+st3oWdPGLrLUQwYUP12VYsySSoWThGUpBxYsACuuy75wLjnnmmnkeq/I49MfleuvBKWL087jSTljgWWJOXArbfC3LmOXknZCiG5FtYHH8A996SdRpJyxwJLkjbQwoVw9dVw6KGw775pp5GKx1FHwT77wK9/DcuWpZ1GknLDAkuSNtANNyQd0YYOTTuJVFwqR7E++siOgpIaDgssSdoAFRVw1VVw9NHw7W+nnUYqPocdlvzu/OY3sGRJ2mkkacPZRVCSNsDVVydF1pVXpp1E2nBrt1Nfe1k+VI5iHXII3HEHDB6cl8NIUsFYYElSHX36KQwfDiedBHvskXYaacOldc2pHj3gu9+F//f/YOBAaN06lRiSlBNOEZSkOvrtb5MpTb/8ZdpJpOIWQtLo4pNPkot1S1Ixy1uBFUK4K4QwJ4TwWpX7Ng0hPB5CeDfzfZN8HV+S8umjj5LW7P37ww47pJ1GKn7du8Mxx8CwYfD552mnkaS6y+cI1kjgiLXuGwI8EWPcHngic1uSik7lqNUvfpFuDqkhGTYsuWi35zRKKmZ5K7BijE8DX6x193FA5eUE7wF65ev4kpQvr7wCd98N//d/0KVL2mmkhmOXXZJzsG65JbkAsSQVo0Kfg7V5jHF25udPgM1rWjGEMCiEMDWEMPWzzz4rTDpJWocY4aKLYNNN4fLL004jNTxDh0KzZvDzn6edRJLqJrUmFzHGCMRalo+IMXaLMXbr2LFjAZNJUs3+8hd44gm44gooLU07jdTwbLll8keMBx+EpbPbpR1HktZboQusT0MInQAy3+cU+PiSVGfLl8PFFydNLc45J+00UsN1ySXQoQPMe3JnYo1/ipWk+qnQBdYEoF/m537AIwU+viTV2YgR8NZbycWFmzdPO43UcG28cdK2fenH7XnoobTTSNL6yWeb9geA54AdQwgzQggDgWHAoSGEd4FDMrclqd6rqEjODenRA44+Ou00UsN31lnQfLMvufhiWLgw7TSSlL18dhH8QYyxU4yxeYyxc4zxzhjj3BjjwTHG7WOMh8QY1+4yKEn10hVXwNy5cM01yUVRJeVX06aw6aGv8/HHSft2SSoWzdIOIEn5csCwycysWFztsrLSEqYM6ZHVfqZNg5tugnPPha5dc5lQUm1adZ7HySfD738PAwbANtuknUiS1s0CS1KDNbNiMdOH9ax2WfmQiVntY9UqOO88aN/ei59KabjqKnjkEbjwQnj44bTTSNK6pdamXZKKwciR8NxzyV/QN9kk7TRS41NWBpddlhRZkyalnUaS1s0CS5JqMHcu/OQncOCBcPrpaaeRGq8LLoAdd0xGk214Iam+s8DS/2/vzuOtrOo9jn9+HIZzFBEQDAUMBxwwVBKcvSqm11lTU9MKM4dMNLt5BRUsu2qmvvJySytz6irlkJBT4pQaWYJTioIgKYqQKUiiVwQOrPvHOsTgAQH3Oc/Ze3/er9fz2ns/e+B7eF57+K21nrUkrcT55+fZA6+5xoktpCK1awe/+AW89hpcdFHRaSRp1SywJKkRY8fmda/OOgv69i06jaS99oKTT4Yf/xiee67oNJK0ck5yIUkrWLywFSedBJtumhc7lVSM7h3rlpuQZtG6rUm1e7HLgfPp9rUn6NG5drVnA5Wk5mKBJUkreO9PWzJ3KjzyCKy7btFppOrVWPF0x+5wzDG1fLvbQfzP26s3G6gkNSeHCErSMsaNg7lPbcZpp8FAG8alFufoo+HQQ2H4cFg4Z52i40jSx1hgSVKD+fPzYqY17T/i8suLTiOpMRF54pm2bWH2vTtQX190IklangWWJDW46CKYNAk2+PcJdOhQdBpJK9OjB1x9Ncyf2cnGEEktjgWWJAGPPw6XXQYnnQR1m79TdBxJn+DLX4Z1tp7J977nrIKSWhYLLElVb84c+MpXYIstYMSIotNIWh0R0Hn/F+naNb9/P/qo6ESSlFlgSapqKcFpp8Fbb8HIkdC+fdGJJK2umrqFXH89TJyYFwaXpJbAAktSVbvpJrjjjrze1YABRaeRtKYOPBDOOAOuugruuafoNJJkgSWpik2ZAmeeCXvvDf/5n0WnkbS2rrwS+vWDQYNg2rSi00iqdhZYkqrS4gU1HHkk1NXBzTdDTU3RiSStrdra3BO9aBEceywsWFB0IknVzAJLUtVJCWbfvx2TJsGtt+YpnyWVt803hxtugPHj4dxzi04jqfmsnVwAAA+gSURBVJpZYEmqOiNGwIcvb8wll8C++xadRlKpHHUUnHVWfo/fcUfRaSRVq9ZFB5Ck5jR2LJxzDtT1foshQ7oVHUfSp9C9Yx29ht633L7UNmi38a4cd0IHttiihn79CgonqWpZYEmqGq+/Dl/6Uh5KNO+g54mwwJLK2RNDBza6/63vQo8t53HYYXU89RR0860uqRk5RFBSVXjvPTj44LwY6ejR0KpdfdGRJDWRbt1gw6Oe5t134YgjXIRYUvOywJJU8RYuhKOPhsmT4c47oU+fohNJamptPzOXm2+GcePglFPy5DaS1BwssCRVtJTg9NPh4Yfh2mud1EKqJkcemRcRv+UWGDas6DSSqoXnYEmqaJdeCtdfDxdcAF//etFpJDW3Cy6A6dPzZ0HXrnD22UUnklTpLLAkVay5z/Ri2MNwwgnwgx8UnUZSESLgmmtg1iz4zndykXXCCUWnklTJHCIoqSLdeCPMeXhbDj88X2/lp51UtWpqYORI2HtvOPFEuP/+ohNJqmT+5JBUce64A04+GWp7vcNtt0GbNkUnklS02lq46y7o2xe++EV44IGiE0mqVBZYkirK6NFw/PGw227Q9cinadeu6ESSWooOHeChh2CbbeDww+3JktQ0LLAkVYybb84LCQ8YAPfeC63aLC46kqQWZoMN4JFH8nINRxwB991XdCJJlcYCS1JF+NnP4Gtfg732ggcfhPXXLzqRpJaqc+dcZC0ZLjhqVNGJJFUSCyxJZS0luPxy+Na34NBDc2t0+/ZFp5LU0nXqlNfH698/L0R+9dVFJ5JUKSywJJWt+no44wwYMgSOPRbuvDOfyC5Jq6Njx1xkHXooDB4M55+fG20k6dOwwJJUlubOhUMOyUMDzz0Xfv1rZwuUtObWWSc3zpx2GvzwhzBoEMyfX3QqSeXMhYYllZ1p03KL88svwy9/madkl6S11bp1bqzp0QOGD4cpU3LR1b170ckklSN7sCSVld//HnbcEaZPhzFjLK4klUYEDBuWC6uXXsqfM2PHFp1KUjmywJJUFurr8/kRBx8MPXvCU0/BvvsWnUpSpTnySBg3Lq+ZNXAgjBjheVmS1owFlqQWb8YM2G+/fH7EySfDX/4CvXsXnUpSperTB8aPhwMPhLPPzud7/uMfRaeSVC48B0tSi5USjBwJZ56ZTzq/6aZ8AnopdO9YR6+hja8w2r1jXWn+EUmFWdV7fHWe+8TQgdx1V56+/Zxz8ppZN96Ye9ElaVUssCS1SG+/Dd/8JoweDbvtlourUvZaPTF0YOleTFKL82ne40sKs4g8ffs++8Dxx+eerJNOgiuuyIsVS1JjHCIoqUVZvDjPDNinT57Q4oor4I9/dEigpOJsu20+L2vIEPjVr2CbbfLSEJ6bJakxFliSWoynn4Zdd4VTT80/aJ59Ng/NqakpOpmkaldbC5ddBs88A716wQknwAEHwMSJRSeT1NJYYEkq3PTp8I1vwE47weuvwy23wGOP5V4sSWpJtt8e/vxn+MlPcq/WdtvB6afnYc2SBBZYkgo0e3buoerdOxdVZ58NkyfnluGIotNJUuNqavK5WVOn5uLquutgiy3g4ovhvfeKTiepaE5yIanZvfVWXlvmmmvggw9gg35/p+2ASYxqO49RP1z+sUtm82rM7pf9gRn/nLfSf8fZACU1pS5dck/W4MEwdCgMHw7fu3gh6/V/jfX6v0ZNbf1yj1/V55mkymGBJanZTJ0KV16ZZwRcuBCOOgq+/3046H+fZdpljc99vKpplmf8c95KnydJzWWrrfKMpxudOJZd5+7J6NFbkiZsySmnwBlnwKab5set7bTxksqLQwQlNamFC2HUKNh//zwU8Kab4MQT81DA22/3PCtJlaNdt7mMGgXPPw8HHZR76jffHI44Ah55xFkHpWphgSWpSUyeDMOGwSab5J6ql1+Giy6CadPg5z/P5ytIUiXabjv4zW/y591558Gf/gRf+ALM+Pk+XHhh7s2XVLkssCSVzGuv5WmMd9gBtt4aLr0UdtwR7rkn33fhhdCtW9EpJal5dO8Ol1ySZ0q95RZo0/n/uPji3Ju/xx75/K3p04tOKanULLAkrbVFi/J0xcOG5UJqs81ya21dHVx1Vf7hcO+9cMghrmUlqXrV1eXZUT9z7HimT88NUXPmwFln5V7+/v3zDITPP58XW5dU3pzkQtJqW7wYXnopD3d5/HF48MH8I6FVK9htN/jRj+CYY/IinJKkj+veHYYMydvkyXDXXXmCjOHD89a1KwwcmIcU7rtv/jx12QqpvFhgSfqXFac9XzSvDQv+0YEFb63P/Dc7s3BmZ+rntQFgo43gsMPyidz77QedOjVNpu4d61Y685bTsEtqCqv63Pmk563V6+4F/fdZn8Fb7cHDD+cJMW67Ld/VrRvsvDPssku+HDAA2rdf42iSmpEFliQ+/BBeeQWmPLk+g/sN5Lnn4Lnn4M03lj5myy1hZu83+Ol/bMIee+ThgM3RquqaMZKaW1N97qzqdXsNvY9Bg2DQoDzb4KRJ8Oij8OSTMG5c7umCPGJg882hb9+l2+c+lz+T27RpktiS1pAFllQF6uth5sx8TtT06fDGG3mbMiUPUXnjX4XUjlx0dy6mdtstr9/Sr1/eunSBXkMnMGjQJkX+KZJU8SLyEhZ9+uTPYYDZs2H8+Ly98AK8+CL87ndLz9mqqYHPfjYXX8tuPXrkYYkbbgit/dUnNYtC3moRcQAwAqgBrkspXVZEDqkcpZR7nObOhfffh3ffhXfeydusWR+/PnNm3lY8cXr99fNMVnvumQuqrbaC74wZyys/3ZN11y3mb5MkNW6DDeDAA/O2xLx5MHFiLrZeeQX+9re83X57/m5YVqtWebjhxhvnbaON8mt27tz4ZadO9ohJa6vZC6yIqAGuBvYD3gSeioi7U0oTmzuLVm7FxRAbWxzxkx6zNs8p8nVX9phFi/JWX7/0+oq3V3XfktsLFsD8+fDRR3lr7Pqy++bNywXU++8vLaaWbKuaZaq2Np8k3bVr7nXaZps8S1XPnstvHTp8/LlDnptrcSVJZaKuLs/guuOOH79vzhx49VWYMSM3si25nDkzr8/15JO5CKuvX/nrt2uXz/dacVt33aXX6+qgbdv82BW3xvbX1CzdWrde/vbq3hexdIPlb69q/8r2SaVWRA/WTsDUlNKrABFxK3A4UBYF1t13w/HHL7+vpRUEa/ocFa9161wY1dbmL6Al1zt0yD1NPXrk6+ut9/HLzp1zIbWkqFpnHb80JKnadeq08uJriZRyo93s2bnYmj176fU5c+CDDxrfZs3Kl++/v7RxcP785vvbmkLRxVhTvn45Zx8zJq8ZV24iNfOv7Yg4GjggpXRyw+2vAjunlAav8LhTgVMbbm4FTG7WoKvWBZhVdAg1OY9zdfA4VwePc3XwOFcHj3N1aInH+bMppa6f9KAWe7pjSula4NqiczQmIp5OKfUvOoealse5Onicq4PHuTp4nKuDx7k6lPNxblXAvzkD6LnM7R4N+yRJkiSprBVRYD0F9I6ITSOiLXAccHcBOSRJkiSppJp9iGBKqT4iBgMPkKdpvyGl9FJz5/iUWuTQRZWcx7k6eJyrg8e5Onicq4PHuTqU7XFu9kkuJEmSJKlSFTFEUJIkSZIqkgWWJEmSJJWIBdanEBFnRsTLEfFSRFxedB41nYj4bkSkiOhSdBaVXkRc0fBefiEiRkdEx6IzqXQi4oCImBwRUyNiaNF5VHoR0TMiHo2IiQ3fyd8uOpOaTkTURMRzEXFv0VnUNCKiY0T8tuG7eVJE7Fp0pjVhgbWWImIf4HBg+5TStsCVBUdSE4mInsD+wBtFZ1GTeQj4XEppO2AKcF7BeVQiEVEDXA0cCPQBvhwRfYpNpSZQD3w3pdQH2AU4w+Nc0b4NTCo6hJrUCGBMSmlrYHvK7HhbYK2904HLUkrzAVJKbxecR03nKuBcwBlhKlRK6cGUUn3DzSfJ6/OpMuwETE0pvZpSWgDcSm4cUwVJKf09pfRsw/X3yT/GuhebSk0hInoABwPXFZ1FTSMi1gf+DbgeIKW0IKX0z2JTrRkLrLW3JbBnRIyLiMcjYkDRgVR6EXE4MCOl9HzRWdRsTgLuLzqESqY7MH2Z22/iD++KFhG9gH7AuGKTqIn8N7nRc3HRQdRkNgXeAW5sGAp6XUSsW3SoNdHs62CVk4h4GOjWyF0XkP/vOpOHIgwAbo+IzZLz3pedTzjO55OHB6rMreo4p5TuanjMBeShRiObM5uk0oiI9sCdwNkppblF51FpRcQhwNsppWciYu+i86jJtAY+D5yZUhoXESOAocDwYmOtPgusVUgpfWFl90XE6cCohoJqfEQsBrqQK26VkZUd54joS25FeT4iIA8bezYidkopvdWMEVUCq3o/A0TEicAhwL42lFSUGUDPZW73aNinChMRbcjF1ciU0qii86hJ7A4cFhEHAbVAh4i4JaX0lYJzqbTeBN5MKS3phf4tucAqGw4RXHu/A/YBiIgtgbbArEITqaRSShNSShumlHqllHqR3/Cft7iqPBFxAHnIyWEppQ+LzqOSegroHRGbRkRb4Djg7oIzqcQit4JdD0xKKf246DxqGiml81JKPRq+k48D/mBxVXkafmdNj4itGnbtC0wsMNIaswdr7d0A3BARLwILgEG2ektl66dAO+Chht7KJ1NK3yw2kkohpVQfEYOBB4Aa4IaU0ksFx1Lp7Q58FZgQEX9t2Hd+Sun3BWaStPbOBEY2NIy9Cny94DxrJKwJJEmSJKk0HCIoSZIkSSVigSVJkiRJJWKBJUmSJEklYoElSZIkSSVigSVJkiRJJWKBJUkqaxFxQUS8FBEvRMRfI2LnojNJkqqX62BJkspWROwKHEJeBHx+RHQhL/wuSVIh7MGSJJWzjYBZKaX5ACmlWSmlmRExraHYIiL6R8RjDde/HxG/ioixEfF6RBwZEZdHxISIGBMRbYr7UyRJlcACS5JUzh4EekbElIi4JiL2Wo3nbA4MBA4DbgEeTSn1BeYBBzddVElSNbDAkiSVrZTSB8COwKnAO8BtEXHiJzzt/pTSQmACUAOMadg/AejVNEklSdXCc7AkSWUtpbQIeAx4LCImAIOAepY2Itau8JQlwwkXR8TClFJq2L8YvxclSZ+SPViSpLIVEVtFRO9ldu0AvA5MI/dsARzV3LkkSdXLljpJUjlrD/wkIjqSe62mkocLbgNcHxH/Re7dkiSpWcTSkRGSJEmSpE/DIYKSJEmSVCIWWJIkSZJUIhZYkiRJklQiFliSJEmSVCIWWJIkSZJUIhZYkiRJklQiFliSJEmSVCL/D47t0VpxYLjnAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "execution_count": 85,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N_scale = (xmax-xmin) / N_bins            # The scale factor between histogram and the fit. Takes e.g. bin width into account.\n",
    "x_gauss = np.linspace(xmin, xmax, 1000)   # Create the x-axis for the plot of the fitted function\n",
    "y_gauss = N_scale*gauss_extended(x_gauss, len(x_sum), 0, 1)                   # Unit Gaussian\n",
    "ax2.plot(x_gauss, y_gauss, '-', color='blue', label='Unit Gauss (no fit)') \n",
    "ax2.legend(loc='upper right')\n",
    "fig2.tight_layout()\n",
    "\n",
    "# Note, we refer to the old plot \"ax2\". Had we done another plot in between, \n",
    "# we would not have been able to plot on top of the old figure with the matlab syntax.\n",
    "\n",
    "if save_plots:\n",
    "    fig2.savefig('Histogram.pdf', dpi=600)\n",
    "    \n",
    "fig2        "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "First make sure that you understand what the Central Limit Theorem (CLT) states! Then, acquaint yourself with the program. Make sure that you read through it, as many\n",
    "of these features will be used onwards. Do you understand why the uniform distribution needs to go from $\\pm \\sqrt 3$ in order to give a distribution with a width of $1$ (i.e. unit) and why you subtract $1$ from the exponential distribution (and how this works at all)?\n",
    "\n",
    "# Questions:\n",
    "\n",
    "1. What is the mean and RMS of the input distributions?\n",
    "\n",
    "2. Why is there a $\\frac{1}{\\sqrt N}$ at the bottom of cell [7] in the line `sum_value = sum_value / np.sqrt(Ntotal)` (when summing up the various contributions to sum)?\n",
    "    Hint: Assume that I always wanted to compare the distribution of sums with a UNIT Gaussian.\n",
    "\n",
    "3. Using a sum of 10 uniform random numbers with mean 0 and width 1, what is the expected \n",
    "    width of the resulting distribution according to CLT? What is the probability of\n",
    "    obtaining a number beyond 3 sigma, i.e. how many numbers did you get beyond 3 sigma?\n",
    "    What would you expect from a true Gaussian distribution?\n",
    "    And what about the same question for 3.5 sigma? And 4.0 sigma?\n",
    "    Put additional counters into the code, and to test any effects in the tails (which have little statistics as it is)\n",
    "    increase the number of experiments run to (much) more than 1000...\n",
    "\n",
    "4. Now try to add 10 exponential (i.e. set `N_exponential=10` and rerun the program). Does that give something Gaussian? What about 1000?\n",
    "    Then try to add 10 cauchy numbers (i.e. set `N_cauchy=10` and rerun the program). Does that give something Gaussian? What about 1000?\n",
    "    If not Gaussian, why do the Cauchy distribution \"ruin\" the Gaussian distribution?\n",
    "    And is this in conflict with the Central Limit Theorem?\n",
    "\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    "5. If one used a trunkated mean of 10 Cauchy numbers (throwing away the top and bottom e.g. 10%),\n",
    "    will the truncated mean of 1000 Cauchy numbers then converge to a Gaussian (possibly not with unit width)?\n",
    "\n",
    "6. How few/many uniform random numbers needs to be added, before the probability\n",
    "    for the sum to follow a Gaussian distribution is greater than 1% (on average)\n",
    "    when using 1000 sums (i.e. `N_experiments=1000`)? Here, a $\\chi^2$ fit is needed.\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.6.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}