Week1/original/CentralLimit/CentralLimit_original.ipynb

{
 "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",
    "- 11-11-2019 (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 similar 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": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # 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": 2,
   "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": 3,
   "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": 4,
   "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": 5,
   "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": 6,
   "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!), and use in the future. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "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 = 0.5*np.tan(pi * (r.uniform() - 0.5))              # Cauchy with mean 0... now produced in a \"complicated\" way!\n",
    "        # x = 0.5*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] = 0.5*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": 8,
   "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. Note that it is important to think about this - to many/few bins or poor choice of range may ruin the plot, and one should try to get a \"regular\" bin widths:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_bins = 100\n",
    "x_ranges = [(-2.5, 2.5), (-1.5, 8.5), (-5.0, 5.0)]"
   ]
  },
  {
   "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": 10,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xl0VuW99//3F9EChhiQQUVoHAAPQp2yHGp5DNrHqmihOBUfHCqWarWeemwV2/4cemwXnaTaqtWjVq0DOKDYqo+0Wqz1OBSPA4YaCIoSUIgIhoAUhOv3B3fyJBAkkJ35/VrrXtz72tN3J65r5eO197UjpYQkSZIkqfE6tXQBkiRJktReGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixlJiJ+FxH/X63lCyJiSURURcSuLVmbJEmS1BzC92CpWkQkYGBKqaxW29XAvimlcdt4rB2BSuDwlNLrmRYqSZIktVKOYKmp9AW6ACXbumNs5H+bkiRJanP8I1YNFhHFEVEeEZdGxNKIeD8ivlFr/Z0RcW1EDAJKc80rIuKZ3PovRsQ/IuLj3L9frLXvzIj4SUQ8D6wG9s61XRsR/527zfCPEbFrRNwbEZW5YxQ2309AkiRJ+mwGLG2r3YBdgH7AeODGiOhRe4OU0lxg/9xiQUrp6IjoCTwO3ADsClwHPL7Js1lnAhOA7sC7ubav59r7AfsALwC/B3oC/wSuyvoCJUmSpO1lwNK2Wgf8OKW0LqX0BFAFDG7AfiOBeSmlP6SUPk0p3Q+8BZxUa5s7U0olufXrcm2/TynNTyl9DDwJzE8p/SWl9CnwIHBQZlcmSZIkNZIBS7WtB3bcpG1HNoaqasty4abaaiCvAcfeg/83KlXtXTaOTFVbWM9+S2p9/6Se5YacW5IkSWoWBizV9h5QuEnbXmwejLbHYuDzm7QNABbVWnZKS0mSJLVpBizVNhX4UUTsGRGdIuLLbLyF76EMjv0EMCgizoiIzhFxOjAE+FMGx5YkSZJahc4tXYBalR/nPn8HegDzgf+TUnqzsQdOKS2LiBOB64GbgTLgxJTSh409tiRJktRa+KJhSZIkScqItwhKkiRJUkYMWJJapYi4I/dC63pvUY2NboiIsoh4IyIObu4aJbV99jWSsmbAktRa3Qkc9xnrjwcG5j4T2PhsnyRtqzuxr5GUIQOWpFYppfQ34KPP2GQUcHfa6EWgICJ2b57qJLUX9jWSstYqZhHs1atXKiwsbOkyJDXCK6+88mFKqXcznrIfdV9OXZ5re3/TDSNiAhv/zzM777zzIfvtt1+zFCgpe/Y1kprL9vY3rSJgFRYWMmvWrJYuQ1IjREQWL6TeplPW01bvtKgppVuBWwGKioqS/Y3UdtnXSGou29vfeIugpLaqHOhfa3lPYHEL1SKp/bKvkbRNDFiS2qrHgLNyM3wdDnycUtrslh1JaiT7GknbpFXcIihJm4qI+4FioFdElANXATsCpJR+BzwBnACUAauBb7RMpZLaMvsaSVkzYElqlVJKY7eyPgEXNlM5ktop+xpJWfMWQUmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsJrB/vvvz8yZM1u6DEmSJElNrMMErMLCQrp27UpeXl7N5y9/+UuD948IysrKtuvcJSUlFBcXb9e+qt/06dM54ogj+NznPsc555xTZ926desYP3483bt3Z8CAATzwwAN11t9www307duXnj17csUVV9RZN3PmTAYPHszOO+/M6NGj+fjjj2vWlZeXU1xcTLdu3Tj44IN58803m+z6JEmS1DZ1mIAF8Mc//pGqqqqaz5e//OWWLknbaZddduH73/8+48eP32zd5MmTmT17NgsXLuTuu+/m3HPPZeHChQC8/PLLXH311TzzzDPMnj2bKVOm1ASw1atXc+qpp3LVVVexdOlSIqJOAJswYQJDhw5l2bJlnH766Zx++unNc7GSJElqMzpUwNqSc845h/Hjx/OlL32J7t27c+6559asO/7448nLywPggAMOIC8vj+9+97s16xcsWEBE8OCDD/L5z3+e7t2787Of/QyAyy+/nLy8PDp16rTZaNn69eu55pprKCwspG/fvlx88cWsW7euZv1Pf/pTdt99d/Lz8ykqKmLRokVN+SNoc4qLixkzZgw9e/bcbN2DDz7IxRdfTEFBAcXFxRxxxBE88sgjNevGjBnD/vvvT79+/TjvvPOYMmUKAH/961/ZZZddOOOMM9h555353ve+x9SpUwGorKxkxowZTJw4ka5du3LJJZfw7rvv8sYbbzTfRUuSJKnVM2DlzJgxg/vuu4+SkhIefPBBXn75ZQCefPJJqqqqAHj99depqqri17/+9Wb7P/DAA7zyyissWbKEY489FoCf/exnVFVVMWDAgM22v+6665g+fTrPP/88ZWVlzJkzhxtuuAGA0tJSrr32Wv77v/+bFStWcNNNN/G5z32uqS693SktLWW//fZj3LhxTJs2jSFDhlBaWlpn3fXXX88VV1xR77q///3vHHfccQwePJiPPvqIiooKysrK6NKlC3l5eQwfPpzy8nL22Wefmn0lSZIkaEDAioj+EfHXiPhnRJRExL/n2ntGxJ8jYl7u3x659oiIGyKiLCLeiIiDm/oiGmr06NEUFBTUfJYuXVqz7qSTTmLAgAEMGDCAYcOGMXfu3G069pVXXkmvXr3o1q0bBx100Fa3v+2227jqqqvo168f3bt358ILL2TatGnAxue91q9fz1tvvcX69es59NBD6dWr17ZdbAe2atUq8vLyKC0tZdGiRXTv3r0mJFevW7hwIfPmzat3XUVFBaWlpXTp0gWAqqqqmnUrV65k3rx5LF++vM6+kiRJEjRsBOtT4NKU0r8BhwMXRsQQYCLwdEppIPB0bhngeGBg7jMBuDnzqrfTo48+yooVK2o+ffr0qVlX+1aznXbaiTVr1mzTsQcOHLhN2y9cuJAzzzyzJuydffbZVFRUADBo0CBuueUWrr32Wnr27Mmpp55KZWXlNh2/I9t5552pqqriH//4B9/5zndYuXJlzW2e1et++ctf8tBDD9W77mtf+xrvvPNOzS2beXl5Nev69+/PBx98wCGHHFJnX0mSJAkaELBSSu+nlP4n930l8E+gHzAKuCu32V3A6Nz3UcDdaaMXgYKI2D3zyptZRHzm+s6dO2/T8fr3789TTz1VE/YqKyvrjJqdc845PP/888yfP5+5c+dyxx13bFfdHdGgQYN46623apbnzJnD4MGDt2tdz5496d27N/vuuy9r1qyhvLwcgLVr1zJ//vyafSVJkiTYxmewIqIQOAh4CeibUnofNoYwoHo4qB+wsNZu5bm2TY81ISJmRcSs6pGb1my33XbLdFruc889lyuvvJLFixeTUmLu3LnMmDEDgPnz5/PMM8+wdu1aOnfuzIYNG8jPz8/s3O3B+vXrWbNmDevXr6/5/umnnwJw2mmnccMNN7BixQqeffZZXnjhBb72ta8BcOqppzJt2jRKSkpYtGgRt99+e81sgEcffTQff/wx9913H6tWreKXv/wlp512GgD5+fkce+yxTJo0iTVr1jB58uSa20klSZKkag0OWBGRBzwMfDel9Fn3q9U31JM2a0jp1pRSUUqpqHfv3g0to1FOOumkOu/Buuuuu7a+U85Pf/pTLrroIvr167fZu5Pqs379+przvPfeezXnvv322wG49NJLGT58OF/60pfIz89nzJgxNc+ErV27lokTJ7LrrrsyaNAgjjjiCM4888ztu+h26g9/+ANdu3Zl0qRJ3HPPPXTt2pVrr70WgEsuuYShQ4fSv39/xo0bx+23307//v0BOOyww7jqqqsYMWIEQ4cO5bTTTqsJWN26dePBBx/k6quvpnfv3mzYsIFJkybVnPOWW25h9uzZ9OjRgylTpjB16tStjmxKkiSpY4mUNss+m28UsSPwJ+CplNJ1ubZSoDil9H7uFsCZKaXBEXFL7vv9m263peMXFRWlWbNmZXA5klpKRLySUipq6Tq2xv5GatvsayQ1l+3tbxoyi2AAtwP/rA5XOY8BZ+e+nw1Mr9V+Vm42wcOBjz8rXEmSJElSe9GQmRmOBM4EZkfEa7m2HwCTgAciYjzwHnBqbt0TwAlAGbAa+EamFUuSJElSK7XVgJVS+jv1P1cFcEw92yfgwkbWJUmSJEltzjbNIihJkiRJ2jIDliRJkiRlpEMErKuvvpqI4MknnwRgxYoV7LTTThQXF7dsYdpu5eXlFBcX061bNw4++OBtekfZDTfcQN++fenZs2edKffnzZvHcccdR8+ePenTpw9nnXUWH3/8cSbnlCRJUsfQIQIWwMCBA3nooYcAeOyxxygsLGzZgtQoEyZMYOjQoSxbtozTTz+95l1WW/Pyyy9z9dVX88wzzzB79mymTJnCAw88AEBlZSVjx47lnXfeYcGCBaxZs4ZLL7200eeUJElSx9FhAtahhx7Kq6++yqeffsrDDz/MySefXLNu/fr1XHPNNRQWFtK3b18uvvhi1q1bB8A777zD0UcfTc+ePdl9990555xzqKqqAmDBggVEBJMnT6ZXr17stddevPTSSy1yfR1JZWUlM2bMYOLEiXTt2pVLLrmEd999lzfeeGOr+z744IOMGTOG/fffn379+nHeeecxZcoUAA455BDOPvtsdtllF7p168YZZ5zBCy+80OhzSpIkqePoMAEL4KijjuKxxx5jyZIl7LPPPjXt1113HdOnT+f555+nrKyMOXPmcMMNNwDwr3/9i29+85u89957lJaW8t5773HllVfWOW5lZSUffPABo0aN4pprrmnWa+qIysrK6NKlC3l5eQwfPpzy8nL22WcfSktLt7pvaWkp++23H9dffz1XXHEFQ4YM2eJ+L7zwAl/4whcafU5JkiR1HB0qYJ1yyin8+7//OyNHjqzTftttt3HVVVfRr18/unfvzoUXXsi0adMA2G+//Rg7dix5eXnk5+dz6qmn8tprr9XZ/4ILLqBz586ccMIJzJ07t9mup6NatWoVeXl5rFy5knnz5rF8+XK6d+9eM7LYkH0XLlzIvHnztrjfK6+8wu9//3v+8z//s9HnlCRJUsfRkBcNtxtf/OIX+bd/+zdOPfVU/v73v9e0L1y4kDPPPJNOnTbmzQ0bNrDbbrsBsHTpUi6++GKee+45Vq1axdq1aznkkEPqHLdnz54A7LTTTqxZs6aZrqb9e++99xgyZAgAn//85ykpKQFg5513pqqqiv79+/PBBx8AsHLlSvLy8rZ6zOp9f/nLXwLwyCOPbLbfggULOPnkk7n77rvZd999G31OSZIkdRwdagQrIpgxYwb77bdfnfb+/fvz1FNPsWLFClasWEFlZWXNSNQVV1xBRDBnzhxWrFjBpEmT2PguZTW1AQMGUFVVRVVVVU24Ath3331Zs2YN5eXlAKxdu5b58+czePDgrR5z0KBBvPXWWzXLc+bMqbPf0qVL+cpXvsK1117Lcccdl8k5JUmS1HF0qIC1Jeeeey5XXnklixcvJqXE3LlzmTFjBrBxlCI/P5+8vDzeeecdbrnllhauVvn5+Rx77LFMmjSJNWvWMHnyZAYMGMCwYcPqbFdcXMzll19ep+3UU09l2rRplJSUsGjRIm6//faa2QA//vhjjjvuOC644ALGjRu3XeeUJElSx2bAAi699FKGDx/Ol770JfLz8xkzZgxLly4F4KqrrmLWrFnssssunHbaaZx00kktXK0AbrnlFmbPnk2PHj2YMmUKU6dOJSLqbLNgwQKWLFlSp+2www7jqquuYsSIEQwdOpTTTjutJmA9+uijvPrqq/zoRz8iLy+v5rMt55QkSVLHFq3hdreioqI0a9asli5DUiNExCsppaKWrmNr7G+kts2+RlJz2d7+xhEsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGrG109dVXbzaFtyRJkiRBBwpYJSUlHHPMMRQUFNCnTx9OPvnkOutnzpzJnnvu2ULVaVuVl5dTXFxMt27dOPjgg3nzzTcbtN/atWu58MIL6dOnD7vuuiuXXXZZnRdH33LLLRQWFtK9e3dGjx7N8uXLa9Z98sknfOtb36Jnz54UFBTw7W9/O/PrkiRJUtvWYQLW6NGjGTlyJMuWLWPu3Ll85StfaemS1AgTJkxg6NChLFu2jNNPP73mXVZbc/PNN/PSSy9RWlrKq6++ytSpU7n33nsBeO2117jssst4/PHHWbJkCevWreP73/9+zb6XXHIJZWVllJSU8OGHHxqwJEmStJkOEbA+/PBDysrKGD9+PDvssAMFBQVMmDChZn1eXh7HH388ixcvrnm57MMPPwzAhg0b+N73vkevXr0YNmwY7733XktdhnIqKyuZMWMGEydOpGvXrlxyySW8++67vPHGG1vd99lnn+WMM86gR48eDBgwgLFjxzJt2jQA/va3v3HUUUex//77061bN77zne/UrPvkk0+46667+M1vfsPuu+9O586dGTp0aJNepyRJktqeDhGwevbsyYABA5gwYQLPPvss//rXv+qsr6qq4sknn2SPPfagqqqKqqqqmlsIH3jgAaZPn05JSQnTpk1j+vTpLXEJqqWsrIwuXbqQl5fH8OHDKS8vZ5999qG0tHSr+276Yu2UEmVlZVtct3z5cj766CPmzp1LRPDoo4/St29fhgwZwiOPPJLdRUmSJKld6BABq1OnTjz99NN07tyZ0aNH07dvX6655poG7fv4448zbtw4+vbty8CBAxk9enQTV6utWbVqFXl5eaxcuZJ58+axfPlyunfvTlVV1Vb3LS4u5r777mPZsmUsWLCAhx56iNWrVwNw1FFHMXPmTGbPnk1VVRW/+93vAFi9ejWVlZWsXbuWt99+m3fffZcbb7yRM888k/fff79Jr1WSJEltS4cIWAD77rsv9957L8uWLePhhx/m+uuv54knntjqfhUVFfTp06dmuW/fvk1Zphpg5513pqqqiv79+/PBBx9wyCGHsHLlSvLy8ra67wUXXMDhhx/OsGHDGDlyJKNHjyY/Px+AAw88kJ///OeMHj2awYMHc9hhhwHQvXt3unXrxvr167n00kvp0qULI0aMYNCgQbz44otNeq2SJElqWzpMwKrWqVMnjjnmGIqLi5kzZ06d9vr06dOHpUuX1iwvWbKkyWvUZ9t3331Zs2YN5eXlwMaZAefPn8/gwYO3uu9OO+3Eb3/7WxYvXkxJSQkRwQEHHFCz/vzzz2f+/PksWrSIL3zhCxQWFrLLLruw9957ExFERM22KaXNbiuUJElSx9YhAtb69eu58sora/4gf/3113nuuec4+OCDa7bZbbfdqKio2CxAnXjiidxzzz0sWbKEefPm8eijjzZr7dpcfn4+xx57LJMmTWLNmjVMnjyZAQMGMGzYsDrbFRcXc/nll9dpW7VqFWVlZWzYsIG//OUv3HbbbXUmPJk9ezYbNmygrKyMK664gvPPPx+AHj16cNRRR3Hdddexbt06nnvuOebNm8fhhx/e9BcsSZKkNqNzSxfQHDp16kRZWRmHHXYYH3/8MX379uVHP/oRRx99dM02gwYN4rzzzuOAAw6gc+fO3HTTTXz1q1/llFNO4eWXX2b//fdnjz32YNSoUaxdu7YFr0aw8X1V48aNo0ePHuy3335MnTq1zugSwIIFCygsLKzTtnr1akaOHMnChQvZY489uOmmmzjiiCNq1v/whz/kmWeeoVu3bowfP77ONO133HEH55xzDgUFBfTr148//OEP7LHHHk16nZIkSWpbojXc4lRUVJRmzZrV0mVIaoSIeCWlVNTSdWyN/Y3UttnXSGou29vfdIhbBCVJkiSpORiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsDaRldffTXjxo1r6TJatYigrKyspcuQJEmSml2HCVglJSUcc8wxFBQU0KdPH04++eQ662fOnMmee+6Z+XnbWiBrK+GovLyc4uJiunXrxsEHH8ybb77ZoP3Wrl3LhRdeSJ8+fdh111257LLL6rws+JZbbqGwsJDu3bszevRoli9fXrPuk08+4Vvf+hY9e/akoKCAb3/725lflyRJktq2DhOwRo8ezciRI1m2bBlz587lK1/5SkuXpEaYMGECQ4cOZdmyZZx++umcfvrpDdrv5ptv5qWXXqK0tJRXX32VqVOncu+99wLw2muvcdlll/H444+zZMkS1q1bV+c9WJdccgllZWWUlJTw4YcfGrAkSZK0mQ4RsD788EPKysoYP348O+ywAwUFBUyYMKFmfV5eHscffzyLFy8mLy+PvLw8Hn74YQA2bNjA9773PXr16sWwYcN47733GnTO5557jry8PH76058yderUmuMuXbq0ZptzzjmH7373u4wZM4a8vDw+//nPU1lZyZ133smXvvSlmu1qj4ItWLCAiGDy5Mn06tWLvfbai5deeqnOuW+66SYGDx5Mfn4+w4cP54MPPgDgySef5MADDyQ/P5999tmHSZMm1YzeHH/88eTl5QFwwAEHkJeXx3e/+92aY86ePZvi4mJ69OjBQQcdxAsvvFCz7rXXXuOggw5il1124Sc/+UmDfj6NUVlZyYwZM5g4cSJdu3blkksu4d133+WNN97Y6r7PPvssZ5xxBj169GDAgAGMHTuWadOmAfC3v/2No446iv33359u3brxne98p2bdJ598wl133cVvfvMbdt99dzp37szQoUOb9DolSZLU9nSIgNWzZ08GDBjAhAkTePbZZ/nXv/5VZ31VVRVPPvkke+yxB1VVVVRVVdXcQvjAAw8wffp0SkpKmDZtGtOnT2/QOYcPH05VVRU/+MEPOP3002uO26dPnzrb3X333XzjG9/g448/5tFHH2XHHXds0PErKyv54IMPGDVqFNdcc01N+/33388vfvELHnroIZYvX84111zDunXrgI1h8Te/+Q3Lli3jz3/+M9dddx1//OMfgY3hq6qqCoDXX3+dqqoqfv3rXwOwcuVKjj32WM444ww+/PBDfvKTnzBmzBhWr15NSokzzjiDs846i4qKCpYtW9ag+hujrKyMLl26kJeXx/DhwykvL2efffahtLR0q/tu+mLtlFLNLZH1rVu+fDkfffQRc+fOJSJ49NFH6du3L0OGDOGRRx7J7qIkSZLULmw1YEXEHRGxNCLerNU2NSJey30WRMRrufbCiPik1rrfNWXxDdWpUyeefvppOnfuzOjRo+nbt2+dUPJZHn/8ccaNG0ffvn0ZOHAgo0ePzrS2ESNGcNJJJ7HDDjtw0EEH0bVr1wbtd8EFF9C5c2dOOOEE5s6dW9P++9//nksvvZRhw4axww47cPTRR9O/f38ARo4cyfDhw9lxxx3Ze++9GTFiBK+99tpWz/WnP/2J3XbbjQkTJrDDDjtwwgkn0KdPH55//nnefvtt5s+fz0UXXcROO+3EJZdcsn0/iG2watUq8vLyWLlyJfPmzWP58uV07969JiB+luLiYu677z6WLVvGggULeOihh1i9ejUARx11FDNnzmT27NlUVVXxu99t/M939erVVFZWsnbtWt5++23effddbrzxRs4880zef//9Jr1WSZIktS0NGcG6EziudkNK6fSU0oEppQOBh4FptVbPr16XUjo/u1IbZ9999+Xee+9l2bJlPPzww1x//fU88cQTW92voqKizqhT3759M61r4MCB27Vfz549Adhpp51Ys2ZNTfvChQspLCysd5+XX36ZESNG0Lt3bwoKCnjkkUdqRrc+y8KFCykpKaGgoKDmM2/ePN5//30qKiooKCioGXnL+udTn5133pmqqir69+/PBx98wCGHHMLKlStrbnH8LBdccAGHH344w4YNY+TIkYwePZr8/HwADjzwQH7+858zevRoBg8ezGGHHQZA9+7d6datG+vXr+fSSy+lS5cujBgxgkGDBvHiiy826bVKkiSpbdlqwEop/Q34qL51ERHAacD9GdfVZDp16sQxxxxDcXExc+bMqdNenz59+tR5bmrJkiXbfL7P0rlz583aunTpwqefflqzXFlZ2eDz9e/fnwULFtS7buzYsYwaNYrFixezYsUKTjzxxM1ui9v4K938mCNGjGDFihU1n9WrV3PWWWfRp08fVqxYURPUtvXnsz323Xdf1qxZQ3l5ObBxZsD58+czePDgre6700478dvf/pbFixdTUlJCRHDAAQfUrD///POZP38+ixYt4gtf+AKFhYXssssu7L333kREnZ9PSmmzn58kSZI6tsY+gzUcWJJSmlerba+IeDUino2I4VvaMSImRMSsiJhVUVHRyDI+2/r167nyyitr/iB//fXXee655zj44INrttltt92oqKjYLCCceOKJ3HPPPSxZsoR58+bx6KOPbtO5d9ttN9566606gWlrBg4cSGlpKZWVlaxevZonn3yywfuec845/OpXv+LNN99k/fr1PPvssyxcuBDY+CxVr1696Ny5MzNnzuSpp56qt95NpzwfOXIkb775Jg8++CCffvopq1atYtq0aSxfvpy99tqLgQMHcuONN7Ju3TomT57c4Fq3V35+PsceeyyTJk1izZo1TJ48mQEDBjBs2LA62xUXF3P55ZfXaVu1ahVlZWVs2LCBv/zlL9x22211JjyZPXs2GzZsoKysjCuuuILzz984CNujRw+OOuoorrvuOtatW8dzzz3HvHnzOPzww5v8eiVJktR2NDZgjaXu6NX7wICU0kHAfwD3RUR+fTumlG5NKRWllIp69+7dyDI+W6dOnSgrK+Owww4jLy+PMWPG8KMf/Yijjz66ZptBgwZx3nnnccABB7Dnnnvy2GOPAXDKKacwatQo9t9/f04++WRGjRq1Tec+7bTTyM/Pp1+/fuy55540JEwecsghnHLKKRx88MGMGzeOI488ssHnO+OMM/iP//gPxowZQ0FBAT/84Q9rbt+76aab+OEPf0h+fj6//e1vOeGEEzbb/6c//SkXXXQR/fr144orrgA2Bponn3yS3/3ud/Tp04fCwkLuueceOnXqRERw7733cuedd9KrVy8EHCw5AAAfA0lEQVSa+ndZ7ZZbbmH27Nn06NGDKVOmMHXq1M1G3xYsWLBZYF69ejUjR44kLy+P888/n5tuuokjjjiiZn31z+eLX/wiJ554Yp1p2u+44w5KS0spKChg/Pjx/OEPf2CPPfZo2guVJElSmxINucUpIgqBP6WUhtZq6wwsAg5JKZVvYb+ZwPdSSrM+6/hFRUVp1qzP3ERSKxcRr6SUilq6jq2xv5HaNvsaSc1le/ubxoxgfRl4q3a4iojeEbFD7vvewEDg7UacQ5IkSZLajIZM034/8AIwOCLKI2J8btXX2Xxyi/8FvBERrwMPAeenlOqdIKOtq36RcH2fXr16tXR5kiRJklrA5lPYbSKlNHYL7efU0/YwG6dtb/eqXyQsSZIkSdUaO8mFJEmSJCnHgCWp1YqI4yKiNCLKImJiPesHRMRfc6+GeCMiNp8aU5K2wr5GUpYMWJJapdyEOTcCxwNDgLERMWSTzX4EPJB7NcTXgZuat0pJbZ19jaSsGbAktVaHAmUppbdTSmuBKcCmL6JLQPW79nYBFjdjfZLaB/saSZkyYElqrfoBC2stl+faarsaGBcR5cATwHfqO1BETIiIWRExqyEv+5bUodjXSMqUAUtSaxX1tG36ZvSxwJ0ppT2BE4A/RMRm/VpK6daUUlFKqah3795NUKqkNsy+RlKmtjpNuzZ35KRnWLTiE/oVdOX5iUe3dDltQvXPDPDnpoYqB/rXWt6TzW/LGQ8cB5BSeiEiugC9gKXNUqGk9sC+RlKmHMHaDotWfMKCSSNrAoO2rvpn5s9N2+AfwMCI2CsidmLjg+WPbbLNe8AxABHxb0AXwPtyJG0L+xpJmTJgSWqVUkqfAhcBTwH/ZOMMXiUR8eOI+Gpus0uBb0bE68D9wDkppU1v7ZGkLbKvkZQ1bxGU1GqllJ5g4wPltduurPV9DnBkc9clqX2xr5GUJQNWRnwuS5IkSZK3CGbE57IkSZIkOYIltRKOgkqSJLV9jmBJrYSjoJIkSW2fAUuSJEmSMuItgtI28lY+SZIkbYkjWNI28lY+SZIkbYkjWNqq6hEboE2M2jjCJEmSpJbiCJa2qnrEpq2M2jjCJEmSpJZiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIw4TXszynq6c6cjlyRJkloXR7CaUdbTnTsduSRJktS6GLAkSZIkKSMGLEmSJEnKiM9gtbDPei4r62e21Lx8Rk6SJKnjMWC1sOrnqAAKJz7e4HVZMQQ0nerfX1P97iRJktT6eItgB+dEGZIkSVJ2DFiSJEmSlBEDliRJkiRlxGew1Co5wYckSZLaIgOWWqXmmOBDkiRJytpWA1ZE3AGcCCxNKQ3NtV0NfBOoyG32g5TSE7l1VwDjgfXAxSmlp5qg7mazrbPs9SvoSuHEx9vFqEvta2/odm39mrdF9e+6+ntHunZJkiTVryEjWHcCvwXu3qR9ckrpl7UbImII8HVgf2AP4C8RMSiltD6DWlvEtk61Xf1HdnsYdak9itSQ7drDNW+L2oGqo127JEmS6rfVgJVS+ltEFDbweKOAKSmlfwHvREQZcCjwwnZXqG3SUUeT2it/n5IkSW1LY2YRvCgi3oiIOyKiR66tH7Cw1jblubbNRMSEiJgVEbMqKirq20TboS2+1+rISc9QOPFxCic+zpGTnmnpclqVtvj7lCRJ6si2N2DdDOwDHAi8D/wq1x71bJvqO0BK6daUUlFKqah3797bWYbag+oQYZCQJElSW7ddswimlJZUf4+I/wL+lFssB/rX2nRPYPF2V6d2raGTaGR9Lm+1kyRJUlPZrhGsiNi91uLXgDdz3x8Dvh4Rn4uIvYCBwMuNK1HtVfXIVXMEHm+1kyRJUnNoyDTt9wPFQK+IKAeuAooj4kA23v63APgWQEqpJCIeAOYAnwIXtuUZBCVJkiRpWzRkFsGx9TTf/hnb/wT4SWOKkiRJkqS2aLuewVLLqX6WCGiWZ5c21Z5epCxJkiRlzYDVxtT38t/q0FP9fUuymOhhW1+knMVEFoY6SZIktRUGrHagoaGjOpw1NBxlob5AuK22NdRJkiRJLcWA1QgNHTna1uM19UhN1nU3l03rdjRLkiRJrY0BqxGy/gO/uUZq2mowqV23o1mSJElqjQxY2i5tfRTMETBJkiQ1BQOWtktbDSc+zyVJkqSmZMBSHS09DbwkSZLUlhmwVEcWs/5JkiRJHZUBqxlk8S4ofbasR96yeGfY9mqrz7dJkiTJgNUstnVUqPZEDNuro0zmUDsIZTny1hLvDKvWnn9fkiRJ7Z0BqxXK4g/sjjKZg7c0SpIkqTUxYCkz7X3UzFs9JUmStDUGLGWmvY+aOVomSZKkrTFgNZHWMlFBa6mjtiyeMWuLNp2Ioz2O8kmSJHV0Bqwm0lr+eG4tddTWGmtqDrVHwNrrKJ8kSVJHZ8BSq9Jco1ubjux11NAnSZKkbBmw1Ko0V9CpfR5HkyRJkpSVTi1dgNTSqkezjpz0TEuXIkmSpDbOESwB2d6a1xon1vgs7X32Q0mSJDUfA5aAbG/N83kmSZIkdVTeIqg2qa1P9e5tiZIkSe2TI1hqk9r6KJm3JUqSJLVPBiy1W23tWTBJkiS1fQYstVttfZRLkiRJbY/PYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJWXAadclSZIETnIhZcJp1yVJkgQGrAbLesrvtv6iXEmSJEmbM2A1UNZTfjuFuCRJktT++AyW1IKqRzK3NprZEZ/xiojjIqI0IsoiYuIWtjktIuZERElE3NfcNUpqH+xvJGXJESypBTV0JLOjPeMVETsANwL/GygH/hERj6WU5tTaZiBwBXBkSml5RPRpmWoltWX2N5KyZsCSMpT1s3od2KFAWUrpbYCImAKMAubU2uabwI0ppeUAKaWlzV6lpPbA/kZSpgxYbYSTYrQNPluXmX7AwlrL5cBhm2wzCCAingd2AK5OKf3f+g4WEROACQADBgzIvFhJbVpm/Y19jSRoQMCKiDuAE4GlKaWhubZfACcBa4H5wDdSSisiohD4J1Ca2/3FlNL5TVB3h+Mf7upgop62tMlyZ2AgUAzsCTwXEUNTSis22zGlW4FbAYqKijY9jqSOLbP+xr5GEjRskos7geM2afszMDSl9AVgLhvvS642P6V0YO5juFKb0dAJJ7b3uI4+bpNyoH+t5T2BxfVsMz2ltC6l9A4b/8fOwGaqT1L7YX8jKVNbHcFKKf0tNzJVu21GrcUXgVOyLUtqfk01Sujo43b5BzAwIvYCFgFfB87YZJtHgbHAnRHRi4238LzdrFVKag/sbyRlKotp2s8Fnqy1vFdEvBoRz0bE8C3tFBETImJWRMyqqKjIoAxJ7UVK6VPgIuApNt52/EBKqSQifhwRX81t9hSwLCLmAH8Fvp9SWtYyFUtqq+xvJGWtUZNcRMQPgU+Be3NN7wMDUkrLIuIQ4NGI2D+lVLnpvt6nLOmzpJSeAJ7YpO3KWt8T8B+5jyRtN/sbSVna7hGsiDibjZNf/J9cx0NK6V/V/0cnpfQKGyfAGJRFoZIkSZLU2m1XwIqI44DLga+mlFbXau+de2EfEbE3Gx8A9R5lSZIkSR1CQ6Zpv5+N05L2iohy4Co2zhr4OeDPEQH/bzr2/wX8OCI+BdYD56eUPmqi2iVJkiSpVWnILIJj62m+fQvbPgw83NiipNaieor16u+SJEnSZ2nUJBdSe+cU65IkSdoWWUzTLkmSJEnCgCVJkiRJmTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZaVDAiog7ImJpRLxZq61nRPw5Iubl/u2Ra4+IuCEiyiLijYg4uKmKlyRJkqTWpKEjWHcCx23SNhF4OqU0EHg6twxwPDAw95kA3Nz4MiVJkiSp9WtQwEop/Q34aJPmUcBdue93AaNrtd+dNnoRKIiI3bMoVpIkSZJas8Y8g9U3pfQ+QO7fPrn2fsDCWtuV59okSZIkqV1rikkuop62tNlGERMiYlZEzKqoqGiCMiRJkiSpeTUmYC2pvvUv9+/SXHs50L/WdnsCizfdOaV0a0qpKKVU1Lt370aUIUmSJEmtQ2MC1mPA2bnvZwPTa7WflZtN8HDg4+pbCSVJkiSpPevckI0i4n6gGOgVEeXAVcAk4IGIGA+8B5ya2/wJ4ASgDFgNfCPjmiVJkiSpVWpQwEopjd3CqmPq2TYBFzamKEmSJElqi5pikgtJkiRJ6pAMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlqdWKiOMiojQiyiJi4mdsd0pEpIgoas76JLUP9jWSsmTAktQqRcQOwI3A8cAQYGxEDKlnu+7AxcBLzVuhpPbAvkZS1gxYklqrQ4GylNLbKaW1wBRgVD3b/Sfwc2BNcxYnqd2wr5GUKQOWpNaqH7Cw1nJ5rq1GRBwE9E8p/emzDhQREyJiVkTMqqioyL5SSW2ZfY2kTBmwJLVWUU9bqlkZ0QmYDFy6tQOllG5NKRWllIp69+6dYYmS2gH7GkmZMmBJaq3Kgf61lvcEFtda7g4MBWZGxALgcOAxHz6XtI3sayRlqnNLF9BaHTnpGRat+IR+BV1buhSpo/oHMDAi9gIWAV8HzqhemVL6GOhVvRwRM4HvpZRmNXOdkto2+xpJmTJgbcGiFZ+wYNLIli5D6rBSSp9GxEXAU8AOwB0ppZKI+DEwK6X0WMtWKKk9sK+RlDUDlqRWK6X0BPDEJm1XbmHb4uaoSVL7Y18jKUs+gyVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRpymXWpD+hV0pXDi4zXfn594dAtXJEmSpNoMWFIbUjtQVQctSZIktR7bHbAiYjAwtVbT3sCVQAHwTaAi1/6D3Av8JEmSJKld2+6AlVIqBQ4EiIgdgEXAI8A3gMkppV9mUqEkSZIktRFZTXJxDDA/pfRuRseTJEmSpDYnq4D1deD+WssXRcQbEXFHRPSob4eImBARsyJiVkVFRX2bSJIkSVKb0uiAFRE7AV8FHsw13Qzsw8bbB98HflXffimlW1NKRSmlot69eze2DEmSJElqcVmMYB0P/E9KaQlASmlJSml9SmkD8F/AoRmcQ5IkSZJavSwC1lhq3R4YEbvXWvc14M0MziFJkiRJrV6j3oMVEd2A/w18q1bzzyPiQCABCzZZJ0mSJEntVqMCVkppNbDrJm1nNqoiSZIkSWqjsppFUJIkSZI6PAOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZ6dzYA0TEAmAlsB74NKVUFBE9galAIbAAOC2ltLyx55IkSZKk1iyrEawRKaUDU0pFueWJwNMppYHA07llSZIkSWrXmuoWwVHAXbnvdwGjm+g8kiRJktRqZBGwEjAjIl6JiAm5tr4ppfcBcv/2yeA8kiRJktSqNfoZLODIlNLiiOgD/Dki3mrITrkwNgFgwIABGZQhSZIkSS2r0SNYKaXFuX+XAo8AhwJLImJ3gNy/S+vZ79aUUlFKqah3796NLUOSJEmSWlyjAlZE7BwR3au/A8cCbwKPAWfnNjsbmN6Y80iSJElSW9DYWwT7Ao9ERPWx7ksp/d+I+AfwQESMB94DTm3keSRJkiSp1WtUwEopvQ0cUE/7MuCYxhxbkiRJktqappqmXZIkSZI6HAOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDlqRWKSKOi4jSiCiLiIn1rP+PiJgTEW9ExNMR8fmWqFNS22d/IylLBixJrU5E7ADcCBwPDAHGRsSQTTZ7FShKKX0BeAj4efNWKak9sL+RlDUDlqTW6FCgLKX0dkppLTAFGFV7g5TSX1NKq3OLLwJ7NnONktoH+xtJmTJgSWqN+gELay2X59q2ZDzw5JZWRsSEiJgVEbMqKioyKlFSO5FZf2NfIwkMWJJap6inLdW7YcQ4oAj4xZYOllK6NaVUlFIq6t27d0YlSmonMutv7GskAXRu6QIkqR7lQP9ay3sCizfdKCK+DPwQOCql9K9mqk1S+2J/IylTjmBJao3+AQyMiL0iYifg68BjtTeIiIOAW4CvppSWtkCNktoH+xtJmTJgSWp1UkqfAhcBTwH/BB5IKZVExI8j4qu5zX4B5AEPRsRrEfHYFg4nSVtkfyMpa94iKKlVSik9ATyxSduVtb5/udmLktQu2d9IypIjWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGTFgSZIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZWS7A1ZE9I+Iv0bEPyOiJCL+Pdd+dUQsiojXcp8TsitXkiRJklqvzo3Y91Pg0pTS/0REd+CViPhzbt3klNIvG1+eJEmSJLUd2x2wUkrvA+/nvq+MiH8C/bIqTJIkSZLamkyewYqIQuAg4KVc00UR8UZE3BERPbawz4SImBURsyoqKrIoQ5IkSZJaVKMDVkTkAQ8D300pVQI3A/sAB7JxhOtX9e2XUro1pVSUUirq3bt3Y8uQJEmSpBbXqIAVETuyMVzdm1KaBpBSWpJSWp9S2gD8F3Bo48uUJEmSpNavMbMIBnA78M+U0nW12nevtdnXgDe3vzxJkiRJajsaM4vgkcCZwOyIeC3X9gNgbEQcCCRgAfCtRlUoSZIkSW1EY2YR/DsQ9ax6YvvLkSRJkqS2K5NZBCVJkiRJBixJkiRJyowBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjDRZwIqI4yKiNCLKImJiU51HUvu1tX4kIj4XEVNz61+KiMLmr1JSW2dfIylLTRKwImIH4EbgeGAIMDYihjTFuSS1Tw3sR8YDy1NK+wKTgZ81b5WS2jr7GklZa6oRrEOBspTS2ymltcAUYFQTnUtS+9SQfmQUcFfu+0PAMRERzVijpLbPvkZSpjo30XH7AQtrLZcDh9XeICImABNyi1URUdpEtQD0Aj7c1p2idf3/qe26hlakrdfP/9/e/YfaXddxHH++ck0JTcURiq6mMMOxfxwjFCSTiYwRjmDFAmnREDSSwn8S9k/oXxoaBIINFC2spN+XKAzTYUh3baRN51pMi7o4UqhWIC2H7/74fhx3d2f3fM7Z9/v9fD/3vh5w4Ht2z855vb/3c97nfL7nc8+XAdYw4RjtOv9HWr6/sX1k/m0i4oSkY8AlLKhzQb85LumVlrP2aXDjcEK154f6a6g9/0dbvj/3mtFqHydQfw3OX95U/aarCdaoozpxypWI3cDujh7/1DDS/ojY2MdjdaX2GmrPD/XXUGH+sX0k8zan9JsK98MpnL+82mtYCvnbvssR/+ZeU3l+qL8G5y9v2n7T1RLBOWD1vOtXAG909FhmtjTl9JGTt5G0ArgQ+Ecv6cxsqXCvMbNWdTXB2geslXSlpJXAdmCmo8cys6Upp4/MADvS9jbg2Yg47aiymdki3GvMrFWdLBFM65O/BDwNnAM8FhEHu3isTL0sRexY7TXUnh/qr6Gq/GfqI5LuBfZHxAzwKPAdSUdojiZvz7jrqvbDCM5fXu01OP887jVnVHt+qL8G5y9vqhrkAzBmZmZmZmbt6OxEw2ZmZmZmZsuNJ1hmZmZmZmYtWTYTLElfl/RHSQck/UTSRaUzTUrSpyUdlPSupGq+9lLSZkmHJR2RdE/pPJOQ9JikN2s9l4mk1ZKek3QojZ0vl87Uh3FjTtK5kp5KP98raU3/KReXUcPdkl5NPe3Xkto+D9lZyX3eS9omKYbW03LyS/pM+h0clPTdvjOOkzGGPpz6w4tpHG0pkfNMxvVfNb6Z6jsgaUPfGVOOqvuNe015tfcb95oRImJZXIBbgBVp+37g/tKZpqjhGpoTnu0BNpbOk5n5HOA14CpgJfAHYF3pXBPk/ziwAXildJYp818GbEjbFwB/qmn/T1nz2DEHfBF4JG1vB54qnXuKGm4CPpC27xxSDbnP+zQmnwdmh9TTMvf/WuBF4OJ0/UOlc09Rw27gzrS9DvhL6dwL8i3af4EtwC9pzlF1HbB3oPt5sP3Gvab8pfZ+414z+rJsPsGKiF9FxIl0dZbmPBdViYhDEXG4dI4JfQw4EhGvR8T/gO8DWwtnyhYRz1PxuU4i4mhE/D5t/wc4BFxeNlXncsbcVuCJtP1DYJOkUScSLWVsDRHxXES8na4OraflPu/vAx4A/ttnuAw5+W8HHo6IfwJExJs9Zxwnp4YAPpi2L2Rg56vM6L9bgW9HYxa4SNJl/aQ7qfZ+415TXu39xr1mhGUzwVrgCzQzUeve5cDf5l2fY+m/wR+ktCzlWmBv2SSdyxlzJ2+TDrwcAy7pJV2eSZ83OxlWTxubX9K1wOqI+HmfwTLl7P+rgaslvSBpVtLm3tLlyanha8BtkuaAXwB39ROtNUN4fam937jXlFd7v3GvGaGT82CVIukZ4NIRP9oVET9Lt9kFnACe7DNbrpwaKjPqKJ3PDdAzSecDPwK+EhH/Lp2nYzljbujjMjufpNuAjcCNnSaazKL5Jb0P+Abw+b4CTShn/6+gWbbzCZoj+r+RtD4i/tVxtlw5NXwWeDwiHpR0Pc15ntZHxLvdx2vFEJ7Htfcb95ryau837jUjLKkJVkTcvNjPJe0APglsirSocmjG1VChOWD1vOtXMLCPhpc6Se+nmVw9GRE/Lp2nBzlj7r3bzElaQbNkYUhLQbOeN5JuBnYBN0bE8Z6y5RiX/wJgPbAnrZS6FJiRdGtE7O8t5ZnljqHZiHgH+LOkwzRvgPb1E3GsnBp2ApsBIuK3ks4DVgFDWn60mCG8vtTeb9xryqu937jXjLBslgimj1O/Ctw6by2xdW8fsFbSlZJW0vyB70zhTMtGWuf/KHAoIh4qnacnOWNuBtiRtrcBzw7soMvYGtKyl2/R9LShvUgtmj8ijkXEqohYExFraP6uY0hveHLG0E9p/vgfSatolvC83mvKxeXU8FdgE4Cka4DzgLd6TXl2ZoDPpW/4ug44FhFHe85Qe79xrymv9n7jXjPK2XzrRk0X4AjN+smX0uWR0pmmqOFTNLPo48DfgadLZ8rMvYXm2+teo1nqWDzTBNm/BxwF3kn7fmfpTBPmv4HmY+wD88b+ltK5eqj7tDEH3EvzwgpNc/9B6gu/A64qnXmKGp5JfeC93+tM6cyT5F9w2z0M75u9xu1/AQ8BrwIvA9tLZ56ihnXACzTf+vUScEvpzAvyn9Z/gTuAO+b9Dh5O9b1cagzV3m/ca8pfau837jWnX5T+o5mZmZmZmZ2lZbNE0MzMzMzMrGueYJmZmZmZmbXEEywzMzMzM7OWeIJlZmZmZmbWEk+wzMzMzMzMWuIJlpmZmZmZWUs8wTIzMzMzM2vJ/wGWOgeBcN16ygAAAABJRU5ErkJggg==\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": 11,
   "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": 12,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtEAAAGDCAYAAADtZ0xmAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu4ZGdZJ+zfQyDQEKDBBMROQnMISEBONgiT8RTAQcJJRxEG+SIyxiOCoyNNxJFv1PkaP+Wg4GgwSOQgMFEOGkGQgIxRAuEYMGCY2EA6gSRImyANIeGZP2o12eneh1qdXbv27n3f17WvXWvVWrWeWqtq71+99a71VncHAACY3s3mXQAAAGw0QjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJETDBlFVH6+q75l3HfNUVT9QVZ+tqi9V1YPmXc9NUVW3rKp/rKpvXmaZ3VX1yLWsay2NPZ5V9T1VdemMatleVV1VNx+m31pVp67SY39nVX1ywfSqHtdZ/22oqvtX1d/P6vFhoxKiYR1Y7J9qVf1YVf3d/unuvm93v3uFx7lREDgM/XaSn+vuo7r7Q/Mu5iY6Lcl7uvtzSVJVr6yq35hzTWtt2eM5vJbvOYe60t3f391nrbTcNDV29//u7nuvRl2LvU6m+dtwU3T3R5PsrarHzWobsBEJ0cDU1kE4v2uSj8+5htXyk0leNe8i5uxwOp6LWgfvmdXymkxes8BAiIYNYmFrdVU9tKouqKqrq+rzVfXCYbH3DL/3Dl+RP7yqblZVz6uqT1fVFVX1J1V1+wWP+/8M932hqn71gO08v6rOrqpXV9XVSX5s2PY/VNXeqrq8ql5aVUcueLyuqp+pqour6pqq+vWqusewztVV9YaFyx/wHBetdej68KUkRyT5SFX9n0XW/XdVdVVVHTdMP2Co8VtX2K/Pqao9Q62frKpHDPNv1OJ3YFeCYT/916r6aFX9W1WdWVV3HroBXFNVf1NVd1him8cnuUeS84fp05I8NckvD8ftLxYs/sBhG/9aVa+vqlsteJzHVtWHh+f591V1/yW2V1X1omGf/uvwePcb7nt3Vf3nBcve6BuQOR7P/a/ljwz75EcW3PeLw+NdXlVPXzD/llX121X1mZq8L/6gqrYsUdsRw7JXVdUlSU454P5v7JequmdV/e2w766qqtcvVeP+18nwuvpckj8+8LUzeEhNuvN8sar+eP9xPXD/LzgG91zqdVI3fs/esqpeXFWXDT8vrqpbDvftr22p/feYoaZrhvfELy0o491JHrH/sQAhGjaqlyR5SXffLpMw9oZh/ncNv7cOX5H/Q5IfG36+N8ndkxyV5KVJUlUnJvn9TP4x3yXJ7ZNsO2BbT0hydpKtmbRGXZ/kF5IcneThSR6R5GcOWOfRSb49ycOS/HKSM4ZtHJfkfkmessTzWrTW7v5qdx81LPOA7r7HgSt2998n+cMkZw3B6VVJntfdn1hiW6mqeyf5uSQP6e7bJvkPSXYvtfwi/mOSRyW5V5LHJXlrktMz2Tc3S/LzS6z3bUku6e7rhtrPyGTf/tZw3BZ+bf6kTPbn3ZLcP5P9k6p6cJJXZNI6+E3Dc3/LEiHn+zJ5bdwrk+P4I0m+MOJ5zuN4fteC+4/q7tcP09+cG16nz0jysgUfVl4wPMcHJrnnsMx/W6K2n0jy2CQPSrIjyQ8t/fTz60nenuQOSY5N8ntT1HjHTFraT1viMZ+ayevtHkPNz1tm+xm2t9zrZL9fyeQ4PTDJA5I89IDHXm7/nZnkJ4f3wv2SnLtg23uSfC3JqnRLgcOBEA3rx5uGFsW9VbU3k3C7lK8luWdVHd3dX+ru9y6z7FOTvLC7L+nuLyV5bpIn1+Rr5h9K8hfd/XfdfW0mgaMPWP8fuvtN3f317t7X3R/o7vd293XdvTuT8PbdB6zzgu6+urs/nuRjSd4+bP9fMwmaS51Etlyt03h+JgHhfUkuS/KyFZa/Psktk5xYVbfo7t3dfVCr6DJ+r7s/PwSM/53k/O7+UHd/Nckbs/Tz3Jrkmim38bvdfVl3/0uSv8gkHCWTEPiH3X1+d18/9N/9aiYB6kBfS3LbJN+apLr7ou6+fMrtJ/M7nov5WpL/3t1f6+6/SvKlJPeuqspkn/xCd/9Ld1+T5H8kefISj/OkJC/u7s8O+/b/W2Gbd03yLd39le7+u2WWTZKvJ/m14cPCviWWeemCbf9mlv4gMtZTM9k/V3T3lUn+3yRPW3D/ovtvwX0nVtXtuvuL3f3BAx77mkxeu0CEaFhPntjdW/f/5ODW3YWekUnr1Seq6v1V9dhllv2WJJ9eMP3pJDdPcufhvs/uv6O7v5yDWyg/u3Ciqu5VVX9ZVZ+rSReP/5FJy+tCn19we98i00dlccvVuqLu/lqSV2bSivY73X3gB4IDl/9UkmdnEr6vqKrXVdW3TLOtwaE+zy9mEmqn8bkFt7+84DHvmuQXD/jgdVwm+/BGuvvcTL59eFmSz1fVGVV1uym3n8zpeC7hC/tb8Af798kxSW6d5AML9sfbhvlL1bbwtf3pJZZLJq3vleR9NbkSxo+vUOOV3f2VFZY5cNtjXnfLWWyfL3zspfZfMvlm5TFJPj10X3n4AY992yR7V6lO2PCEaNiAuvvi7n5Kkjtl8hX22VV1mxzcipxMWmTvumD6+CTXZRKELs/k6+kkydAN4psO3NwB0/8zySeSnDB0Jzk9k4CxGpardUVVtS3JryX54yS/M03/ze5+bXf/+2G7ncn+TJJ/yySU7bfkpegOwUeT3P2AFtllA/8iPpvkNxd+8OruW3f3ny62cHf/bnd/e5L7ZvIB7L8Od83yed6k4znSVZkE+vsu2B+3X9Bt5ECXZ/KhY2Fti+ruz3X3T3T3t2TSfeb3a/krckxzLA/c9mXD7Rsdjzr4EogrPfZi+/yyJZa98QN3v7+7n5DJ35U35YZuYhk+XB6Z5JNLrA6bjhANG1BV/WhVHdPdX88NLUPXJ7kyk6+S775g8T9N8gtVdbeqOiqTluPXD61RZyd5XE1Oyjsyk69+VwrEt01ydZIv1eSkvZ9etSe2fK3LGr7Of2Um/TqfkUlI+vUV1rl3VZ08hO2vZBLCrh/u/nCSx1TVHYcg8+xDfE4H6e5Lk1ycSX/V/T6fGx+3lbw8yU9V1XfUxG2q6pSqOqiFu6oeMix3i0xC2ldy4+f5g1V16yEYPuNQntMSDvl4DqbeJ8N74eVJXlRVd0omH6qq6j8sscobkvx8VR079AneudRjV9UPV9X+D5tfzCTI7t9/Y4/bfj87bPuOmXwQ3d+f+iNJ7ltVDxxONnz+AeuttL0/TfK8qjqmqo7OpIvWq1cqpqqOrKqnVtXth290rs4NzzFJvifJuUNXJSBCNGxUj07y8Zpc4eAlSZ489NX8cib9K88bvtJ+WCYnn70qkyt3/HMmAeqZSTL0cX1mktdlEjqvSXJFJn1rl/JLSf7TsOzLc8M//9WwZK1T+PlMugn86tCN4+lJnl5V37nMOrdMsiuTVszPZdICd/pw36syCTS7MzmpbDWfZzLpS76wr+qZmfRH3VtVb1pp5e6+IJM+wC/NJNh9KsNJh4u4XSbH6ouZfL3/hUyu0ZwkL0pybSbh7KxMTlxbLTfleCaTAHnWsE+eNMXyz8lkP7x36Gr0N1n6RLiXJ/nrTI7xB5P8+TKP+5Ak5w/vt7ckeVZ3//Mh1rjfazN5XV0y/PxGknT3PyX570PtFyc5sP/1Sq+T30hyQSbfdlw4PLdprz/+tCS7h333U0l+dMF9T03yB1M+DmwKtUKXQWATGVoL92bSVeOfV1qeQze0fn8oySNGnuQHa6qqvi3JGd19YB9p2NSEaNjkajIK2Tsz6cbxO0m+I8mDVzopDwA2M905gCdkcuLRZUlOyKRryGEToKvq+JoMTLHYz5InkwHAcrREAwDASFqiAQBgJCEaAABGuilDr66Zo48+urdv3z7vMgAAOIx94AMfuKq7lxrp9EY2RIjevn17LrjggnmXAQDAYayqPj3tsrpzAADASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASDefdwEA3OCkXedmz959B83ftnVLztt58hwqAmAxQjTAOrJn777s3nXKQfO37zxnDtUAsBTdOQAAYCQhGgAARhKiAQBgJCEaAABGEqIBAGAkIRoAAEYSogEAYCQhGgAARhKiAQBgJCMWAhymDCEOMDtCNMBhyhDiALOjOwcAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIN5/lg1fV7iTXJLk+yXXdvaOq7pjk9Um2J9md5End/cVZ1gEAAKtpLVqiv7e7H9jdO4bpnUne2d0nJHnnMA0AABvGPLpzPCHJWcPts5I8cQ41AADAIZt1iO4kb6+qD1TVacO8O3f35Uky/L7TYitW1WlVdUFVXXDllVfOuEwAAJjeTPtEJzmpuy+rqjsleUdVfWLaFbv7jCRnJMmOHTt6VgUCAMBYM22J7u7Lht9XJHljkocm+XxV3SVJht9XzLIGAABYbTML0VV1m6q67f7bSb4vyceSvCXJqcNipyZ586xqAACAWZhld447J3ljVe3fzmu7+21V9f4kb6iqZyT5TJIfnmENAACw6mYWorv7kiQPWGT+F5I8YlbbBQCAWZv1iYUArIJtW7dk+85zlrzvvJ0nr3FFAJubEA2wASwXkpcK1wDMzjwGWwEAgA1NiAYAgJGEaAAAGEmIBgCAkZxYCLDBLXXljm1bt8yhGoDNQYgG2OBc3g5g7enOAQAAIwnRAAAwkhANAAAj6RMNMEMn7To3e/buO2i+oboBNjYhGmCG9uzdl927TjlovqG6ATY23TkAAGAkIRoAAEYSogEAYCQhGgAARnJiIcAcbKShul1hBOBgQjTAHGyk8OkKIwAH050DAABGEqIBAGAkIRoAAEYSogEAYCQhGgAARhKiAQBgJCEaAABGEqIBAGAkIRoAAEYyYiHAKlhuaOz1ZrkhxzfSSIoA8yREA6yCpYbGXo+WCsqG8QaYnu4cAAAwkhANAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADCSEA0AACMZ9huAVXXSrnOzZ+++g+Zv27plySHHATYaIRqAVbVn777s3nXKQfO37zxnDtUAzIbuHAAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADDSzEN0VR1RVR+qqr8cpu9WVedX1cVV9fqqOnLWNQAAwGpai5boZyW5aMH0C5K8qLtPSPLFJM9YgxoAAGDVzDREV9WxSU5J8kfDdCU5OcnZwyJnJXniLGsAAIDVNuuW6Bcn+eUkXx+mvynJ3u6+bpi+NMm2GdcAAACramYhuqoem+SK7v7AwtmLLNpLrH9aVV1QVRdceeWVM6kRAAAOxSxbok9K8viq2p3kdZl043hxkq1VtX+48WOTXLbYyt19Rnfv6O4dxxxzzAzLBACAcWYWorv7ud19bHdvT/LkJOd291OTvCvJDw2LnZrkzbOqAQAAZmEe14l+TpL/UlWfyqSP9JlzqAEAAA7ZzVde5Kbr7ncnefdw+5IkD12L7QIAwCwYsRAAAEYSogEAYCQhGgAARhKiAQBgJCEaAABGEqIBAGCkNbnEHQDr37atW7J95zmLzl+N5QEOJ0I0AEmS83aePNPlAQ4nunMAAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEg3n3cBABvJSbvOzZ69+w6av23rljlUA8C8CNEAI+zZuy+7d50y7zIAmDPdOQAAYCQhGgAARhKiAQBgJCEaAABGEqIBAGAkIRoAAEYSogEAYCQhGgAARhKiAQBgpKlGLKyq+3X3x2ZdDMBaW24Y7/N2njyHigDYCKYd9vsPqurIJK9M8tru3ju7kgDWzlLDeG/fec4cqgFgo5iqO0d3//skT01yXJILquq1VfWomVYGAADr1NR9orv74iTPS/KcJN+d5Her6hNV9YOzKg4AANajqUJ0Vd2/ql6U5KIkJyd5XHffZ7j9ohnWBwAA6860faJfmuTlSU7v7m+cgdPdl1XV82ZSGQAArFPThujHJNnX3dcnSVXdLMmtuvvL3f2qmVUHMCfbtm5Z9OTCbVu3zKEaANabaUP03yR5ZJIvDdO3TvL2JP9uFkUBzJvL2wGwnGlPLLxVd+8P0Blu33o2JQEAwPo2bYj+t6p68P6Jqvr2JAePTgAAAJvAtN05np3kf1XVZcP0XZL8yGxKAgCA9W2qEN3d76+qb01y7ySV5BPd/bWZVgYAAOvUtC3RSfKQJNuHdR5UVenuP5lJVQBsGiftOjd79h7cQ3Db1i2LnuA5dnmAWZgqRFfVq5LcI8mHk1w/zO4kQjQAN8mevfuye9cpB81f7BKDh7I8wCxM2xK9I8mJ3d3TPnBV3SrJe5LcctjO2d39a1V1tySvS3LHJB9M8rTuvnZc2QAAMD/TXp3jY0m+eeRjfzXJyd39gCQPTPLoqnpYkhckeVF3n5Dki0meMfJxAQBgrqZtiT46yT9W1fsyCcdJku5+/FIrDK3W+68tfYvhp5OcnOQ/DfPPSvL8JP9zVNUAADBH04bo5x/Kg1fVEUk+kOSeSV6W5P8k2dvd1w2LXJpk2xLrnpbktCQ5/vjjD2XzAAAwE1N15+juv02yO8kthtvvz6Q/80rrXd/dD0xybJKHJrnPYostse4Z3b2ju3ccc8wx05QJAABrYqoQXVU/keTsJH84zNqW5E3TbqS79yZ5d5KHJdlaVftbwI9NctlS6wEAwHo07YmFP5vkpCRXJ0l3X5zkTsutUFXHVNXW4faWJI9MclGSdyX5oWGxU5O8eXzZAAAwP9P2if5qd19bVUmSoSV5pcvd3SXJWUO/6JsleUN3/2VV/WOS11XVbyT5UJIzD610AACYj2lD9N9W1elJtlTVo5L8TJK/WG6F7v5okgctMv+STPpHAwDAhjRtd46dSa5McmGSn0zyV0meN6uiAABgPZuqJbq7v57k5cMPAABsalOF6Kr65yzSB7q7777qFQEAwDo3bZ/oHQtu3yrJDye54+qXAwAA69+0g618YcHPnu5+cSbDdwMAwKYzbXeOBy+YvFkmLdO3nUlFAACwzk3bneN3Fty+LpMhwJ+06tUAzMhJu87Nnr37Dpq/beuWOVQDwEY37dU5vnfWhQDM0p69+7J71ynzLgOAw8S03Tn+y3L3d/cLV6ccAABY/8ZcneMhSd4yTD8uyXuSfHYWRQEAwHo2bYg+OsmDu/uaJKmq5yf5X939n2dVGAAArFfTDvt9fJJrF0xfm2T7qlcDAAAbwLQt0a9K8r6qemMmIxf+QJI/mVlVAACwjk17dY7frKq3JvnOYdbTu/tDsysLAADWr2m7cyTJrZNc3d0vSXJpVd1tRjUBAMC6NlWIrqpfS/KcJM8dZt0iyatnVRQAAKxn07ZE/0CSxyf5tyTp7sti2G8AADapaU8svLa7u6o6SarqNjOsCYDD0LatW7J95zmLzgfYaKYN0W+oqj9MsrWqfiLJjyd5+ezKAuBwc97Ok+ddAsCqmfbqHL9dVY9KcnWSeyf5b939jplWBgAA69SKIbqqjkjy1939yCSCMwAAm96KJxZ29/VJvlxVt1+DegAAYN2btk/0V5JcWFXvyHCFjiTp7p+fSVUAALCOTRuizxl+AABg01s2RFfV8d39me4+a60KAgCA9W6lPtFv2n+jqv5sxrUAAMCGsFKIrgW37z7LQgAAYKNYKUT3ErcBAGDTWunEwgdU1dWZtEhvGW5nmO7uvt1MqwMAgHVo2RDd3UesVSEAALBRrDjYCgAAcGNCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMtOyw3wDr1Um7zs2evfsOmr9t65act/PkOVQEwGYiRAMb0p69+7J71ykHzd++85w5VAPAZqM7BwAAjCREAwDASEI0AACMJEQDAMBIQjQAAIw0sxBdVcdV1buq6qKq+nhVPWuYf8eqekdVXTz8vsOsagAAgFmYZUv0dUl+sbvvk+RhSX62qk5MsjPJO7v7hCTvHKYBAGDDmFmI7u7Lu/uDw+1rklyUZFuSJyQ5a1jsrCRPnFUNAAAwC2vSJ7qqtid5UJLzk9y5uy9PJkE7yZ3WogYAAFgtMx+xsKqOSvJnSZ7d3VdX1bTrnZbktCQ5/vjjZ1cgcFjZtnXLoqMWbtu6ZQ7VcFM4lsB6NtMQXVW3yCRAv6a7/3yY/fmqukt3X15Vd0lyxWLrdvcZSc5Ikh07dvQs6wQOH+ftPHneJbBKHEtgPZvl1TkqyZlJLuruFy646y1JTh1un5rkzbOqAQAAZmGWLdEnJXlakgur6sPDvNOT7Eryhqp6RpLPJPnhGdYAAACrbmYhurv/LslSHaAfMavtAgDArBmxEAAARhKiAQBgJCEaAABGEqIBAGAkIRoAAEYSogEAYKSZD/sNcKhO2nVu9uzdt+h9hn5mNSz1Gtu2dYsRE4FlCdHAurVn777s3nXKvMvgMLbUa2z7znPmUA2wkejOAQAAIwnRAAAwkhANAAAjCdEAADCSEwsBOOwtdxUOgEMhRANw2HOlF2C16c4BAAAjCdEAADCSEA0AACPpEw3MnZO+ANhohGhg7pz0BcBGozsHAACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIx083kXAGweJ+06N3v27jto/ratW+ZQDQAcOiEaWDN79u7L7l2nzLsMALjJdOcAAICRhGgAABhJiAYAgJGEaAAAGMmJhcCqcxUONrptW7dk+85zRi1/3s6TZ1gRsN4I0cCqcxUONrqxgXhM4AYOD7pzAADASEI0AACMJEQDAMBIQjQAAIzkxEIADgvLXVHDlWGA1SZEA3BYcIk5YC3NrDtHVb2iqq6oqo8tmHfHqnpHVV08/L7DrLYPAACzMss+0a9M8ugD5u1M8s7uPiHJO4dpAADYUGYWorv7PUn+5YDZT0hy1nD7rCRPnNX2AQBgVtb66hx37u7Lk2T4faelFqyq06rqgqq64Morr1yzAgEAYCXr9hJ33X1Gd+/o7h3HHHPMvMsBAIBvWOsQ/fmqukuSDL+vWOPtAwDATbbWIfotSU4dbp+a5M1rvH0AALjJZnmJuz9N8g9J7l1Vl1bVM5LsSvKoqro4yaOGaQAA2FBmNthKdz9libseMattAgDAWjBiIbCik3admz179x00f9vWLUaJA2BTEqKBFe3Zuy+7d51y0PztO8+ZQzUAMH/r9hJ3AACwXgnRAAAwkhANAAAjCdEAADCSEwsB4CbatnXLoifauoINHL6EaAC4iZYKyq5gA4cv3TkAAGAkIRoAAEYSogEAYCQhGgAARhKiAQBgJCEaAABGEqIBAGAkIRoAAEYSogEAYCQjFgLAjBgOHA5fQjQAzIjhwOHwpTsHAACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIwkRAMAwEhCNAAAjCREAwDASEI0AACMJEQDAMBIQjQAAIx083kXAGxc27Zuyfad5yw6H1jaUu+d5ZY/b+fJM6wIGEuIBg6Zf+pwaMa+d8YEbmBt6M4BAAAjCdEAADCSEA0AACMJ0QAAMJITC+EwdtKuc7Nn776D5i91pv9yywMbx1Lv5cSVPmC1CNFwGNuzd1927zrloPlLnem/1PLAxrLce9mVPmB16M4BAAAjCdEAADCSEA0AACPpEw1zMs8TfwzXDRvLoQwTPvaxnHAI4wjRMCfzPPHHP0rYWFbzPbvUYznhEMaZS3eOqnp0VX2yqj5VVTvnUQMAAByqNQ/RVXVEkpcl+f4kJyZ5SlWduNZ1AADAoZpHS/RDk3yquy/p7muTvC7JE+ZQBwAAHJJ5hOhtST67YPrSYR4AAGwI8zixsBaZ1wctVHVaktOGyS9V1SdnWtXijk5y1Ry2y/IOm+NSLzi0+1ZrG6vosDkmhxnHZf1Z18dkjf5erEfr+rhsUvM6JneddsF5hOhLkxy3YPrYJJcduFB3n5HkjLUqajFVdUF375hnDRzMcVl/HJP1yXFZfxyT9clxWX82wjGZR3eO9yc5oaruVlVHJnlykrfMoQ4AADgka94S3d3XVdXPJfnrJEckeUV3f3yt6wAAgEM1l8FWuvuvkvzVPLY90ly7k7Akx2X9cUzWJ8dl/XFM1ifHZf1Z98ekug86pw8AAFjGXEYsBACAjUyInkJVPXMYpvzjVfVb866HG1TVL1VVV9XR865ls6uq/7+qPlFVH62qN1bV1nnXtFlV1aOHv1mfqqqd866HpKqOq6p3VdVFw/+SZ827Jiaq6oiq+lBV/eW8a2GiqrZW1dnD/5SLqurh865pMUL0CqrqezMZUfH+3X3fJL8955IYVNVxSR6V5DPzroUkyTuS3K+775/kn5I8d871bEpVdUSSlyX5/iQnJnlKVZ0436pIcl2SX+zu+yR5WJKfdVzWjWcluWjeRXAjL0nytu7+1iQPyDo9PkL0yn46ya7u/mqSdPcVc66HG7woyS9nkcF6WHvd/fbuvm6YfG8m14Bn7T00yae6+5LuvjbJ6zJpCGCOuvvy7v7gcPuaTEKB0XrnrKqOTXJKkj+ady1MVNXtknxXkjOTpLuv7e69861qcUL0yu6V5Dur6vyq+tuqesi8CyKpqscn2dPdH5l3LSzqx5O8dd5FbFLbknx2wfSlEdbWlaranuRBSc6fbyUkeXEmjTFfn3chfMPdk1yZ5I+HbjZ/VFW3mXdRi5nLJe7Wm6r6myTfvMhdv5LJPrpDJl+/PSTJG6rq7u2yJjO3wnE5Pcn3rW1FLHdMuvvNwzJM9bqxAAADKElEQVS/kslX169Zy9r4hlpknr9X60RVHZXkz5I8u7uvnnc9m1lVPTbJFd39gar6nnnXwzfcPMmDkzyzu8+vqpck2ZnkV+db1sGE6CTd/cil7quqn07y50Nofl9VfT2T8dyvXKv6NquljktVfVuSuyX5SFUlk24DH6yqh3b359awxE1nufdKklTVqUkem+QRPmjOzaVJjlswfWySy+ZUCwtU1S0yCdCv6e4/n3c95KQkj6+qxyS5VZLbVdWru/tH51zXZndpkku7e/83NWdnEqLXHd05VvamJCcnSVXdK8mRSa6aa0WbXHdf2N136u7t3b09kzfcgwXo+aqqRyd5TpLHd/eX513PJvb+JCdU1d2q6sgkT07yljnXtOnV5BP/mUku6u4Xzrseku5+bncfO/wfeXKScwXo+Rv+l3+2qu49zHpEkn+cY0lL0hK9slckeUVVfSzJtUlO1cIGi3ppklsmecfwDcF7u/un5lvS5tPd11XVzyX56yRHJHlFd398zmUxafV8WpILq+rDw7zThxF8gRt7ZpLXDA0BlyR5+pzrWZQRCwEAYCTdOQAAYCQhGgAARhKiAQBgJCEaAABGEqIBAGAkIRpgg6mqX6mqj1fVR6vqw1X1HfOuCWCzcZ1ogA2kqh6eyaiQD+7ur1bV0ZkMAgXAGtISDbCx3CXJVd391STp7qu6+7Kq2j0E6lTVjqp693D7+VV1VlW9fVjmB6vqt6rqwqp62zAUNQAjCdEAG8vbkxxXVf9UVb9fVd89xTr3SHJKkickeXWSd3X3tyXZN8wHYCQhGmAD6e4vJfn2JKcluTLJ66vqx1ZY7a3d/bUkF2YyFPjbhvkXJtk+m0oBDm/6RANsMN19fZJ3J3l3VV2Y5NQk1+WGhpFbHbDK/q4fX6+qr3V3D/O/Hv8HAA6JlmiADaSq7l1VJyyY9cAkn06yO5MW6iT5j2tdF8BmowUCYGM5KsnvVdXWTFqfP5VJ1477JDmzqk5Pcv4c6wPYFOqGb/UAAIBp6M4BAAAjCdEAADCSEA0AACMJ0QAAMJIQDQAAIwnRAAAwkhANAAAjCdEAADDS/wVmErG4aAXFnAAAAABJRU5ErkJggg==\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": 13,
   "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": 14,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XmclWX9//HXxSIMgoyCGg7iuItbpLilVII7LmgI5gZI4tIXy63QNCnrF5kLbrmlIooKEilGUiaaiVtomGtuobIoiIzs+/X74z5DA84Mh+Gcc8+ZeT0fj3nMOefe3mdBz2eu6/7cIcaIJEmSJGnjNUk7gCRJkiQ1FBZYkiRJkpQjFliSJEmSlCMWWJIkSZKUIxZYkiRJkpQjFliSJEmSlCMWWJKKUgjhzRDCd9LOkaYQwokhhE9CCAtDCN9IO8/GCCG0CCG8FUL4Wi3rTAshHFbIXIW0oe9nCOE7IYTpecpSHkKIIYRmmftPhBD65Wjf3UII/6lyP6fva77/2xBC2DuE8Hy+9i+p+FlgSap3qvvCFULoH0J4rvJ+jHGPGOMz69nPWl8SG6Brgf+LMbaOMf4r7TAbaRDwbIzxU4AQwogQwi9TzlRotb6fmc/yTinkIsZ4dIzxvvWtl03GGOM/Yoy75iJXdZ+TbP7bsDFijP8GKkIIx+XrGJKKmwWWJNVRPSjctgPeTDlDrpwD3J92iJQ1pPezWvXg30yujCL5zErSV1hgSSpKVUe5Qgj7hxCmhBDmhxA+CyFcn1nt2czvisy0q4NCCE1CCFeEED4KIcwOIYwMIbStst8zM8vmhhCuXOc4Q0MIY0MID4QQ5gP9M8d+IYRQEUKYFUK4JYSwSZX9xRDC+SGE90IIC0IIV4cQdsxsMz+EMKbq+us8x2qzZqbTLQSaAq+FED6oZttvhhA+DyFsm7n/9UzG3dbzuv4khDAjk/U/IYQemcfXGilYd3pa5nW6NITw7xDCohDC3SGErTNTyxaEEP4WQti8hmN2AnYEXsrcHwScBvw48749XmX1LpljfBlCGB1CaFllP8eGEKZmnufzIYS9azheCCHckHlNv8zsb8/MsmdCCN+vsu5aI6cpvp+Vn+XXMq9J3yrLLs7sb1YIYUCVx1uEEK4NIXwckn8Xt4cQSmrI1jSz7uchhA+BnussX/O6hBB2CiH8PfPafR5CGF1TxsrPSeZz9Slw77qfnYz9QjJFdF4I4d7K93Xd17/Ke7BTTZ+TsPa/2RYhhOEhhJmZn+EhhBaZZZXZanr9jslkWpD5N3FJlRjPAD0q9yVJVVlgSWoIbgRujDFuRvJFfUzm8W9lfpdmpl29APTP/BwK7AC0Bm4BCCHsDvyO5EtbB6AtULbOsU4AxgKlJH/FXgVcCLQHDgJ6AOevs81RwL7AgcCPgTszx9gW2BP4Xg3Pq9qsMcZlMcbWmXW+HmPccd0NY4zPA3cA92W+VN8PXBFjfKeGYxFC2BX4P2C/GGMb4EhgWk3rV+O7wOHALsBxwBPA5SSvTRPgghq22wv4MMa4MpP9TpLX9prM+1Z1KlYfktdze2BvkteHEMI+wD0kowrtMs99fA1fgI8g+WzsQvI+9gXmbsDzTOP9/FaV5a1jjKMz97/G/z6nA4FbqxSyv8k8xy7ATpl1flZDtrOBY4FvAF2B3jU/fa4G/gpsDnQEbs4i4xYkI3SDatjnaSSftx0zma+o5fhkjlfb56TST0nepy7A14H919l3ba/f3cA5mX8LewKTqhx7BrACyMlUR0kNiwWWpPrq0cxIREUIoYKk8KnJCmCnEEL7GOPCGOOLtax7GnB9jPHDGONC4DLglJBMXeoNPB5jfC7GuJzky2hcZ/sXYoyPxhhXxxiXxBhfiTG+GGNcGWOcRvLF/tvrbPObGOP8GOObwBvAXzPH/5KkCKmpoUFtWbMxlOTL48vATODW9ay/CmgB7B5CaB5jnBZj/MpoSi1ujjF+lvny+Q/gpRjjv2KMy4A/UvPzLAUWZHmMm2KMM2OMXwCPk3xxhqRAuCPG+FKMcVXmfKFlJF+u17UCaAPsBoQY49sxxllZHh/Sez+rswL4RYxxRYzxz8BCYNcQQiB5TS6MMX4RY1wA/D/glBr20wcYHmP8JPPa/no9x9wO2CbGuDTG+Fwt6wKsBq7KFJJLaljnlirH/hU1F6kb6jSS12d2jHEO8HPgjCrLq339qizbPYSwWYxxXozx1XX2vYDksytJa7HAklRf9Yoxllb+8NVRoaoGkvzV+50Qwj9DCMfWsu42wEdV7n8ENAO2ziz7pHJBjHExXx3Z+KTqnRDCLiGEP4UQPg3JtMH/RzJiU9VnVW4vqeZ+a6pXW9b1ijGuAEaQ/PX9uhjjusXiuuu/D/yIpDCbHUJ4OISwTTbHyqjr85xHUvBk49MqtxdX2ed2wMXrFOXbkryGa4kxTiIZtbwV+CyEcGcIYbMsjw8pvZ81mFs58pdR+ZpsCbQCXqnyekzMPF5Ttqqf7Y9qWA+SUbsAvBySjn1nrSfjnBjj0vWss+6xN+RzV5vqXvOq+67p9YNkRPYY4KPMlMiD1tl3G6AiRzklNSAWWJKKXozxvRjj94CtSKZFjQ0hbMpXR58gGcnZrsr9TsBKki/Js0imPAGQmVrXbt3DrXP/NuAdYOfMFMXLSb585kJtWdcrhFAGXAXcC1yXzfkiMcYHY4yHZI4bSV5PgEUkX9gr1dhOvQ7+DeywzkhOrcVgNT4BflW1KI8xtooxPlTdyjHGm2KM+wJ7kBTnl2YW5fN5btT7uYE+Jyn29qjyerStMhVxXbNICtKq2aoVY/w0xnh2jHEbkimZvwu1dw7M5r1c99gzM7fXej/CV9v4r2/f1b3mM2tYd+0dx/jPGOMJJP9deZT/TT0m84eHTYD/1LC5pEbMAktS0QshnB5C2DLGuJr//UV5FTCHZHrSDlVWfwi4MISwfQihNcmI0+jMX7HHAseFpEHEJiTTidZXLLUB5gMLQ9JA4rycPbHas9YqM0VsBMl5JANJvkBfvZ5tdg0hdM8UYktJvqCvyiyeChwTQtgi8yX3R3V8Tl8RY5wOvEdyfkylz1j7fVufu4BzQwgHhMSmIYSeIYSvjIyFEPbLrNec5Av8UtZ+nieFEFplioaBdXlONajz+5mR9WuS+bdwF3BDCGErSAruEMKRNWwyBrgghNAxcw7SkJr2HUI4OYRQ+YeIeSRFTuXrt6HvW6UfZI69BckfKSrP33oN2COE0CXT+GLoOtut73gPAVeEELYMIbQnmfb7wPrChBA2CSGcFkJomxkJns//niPAd4BJmemvkrQWCyxJDcFRwJsh6cR2I3BK5tyQxSTnc0zOTJM6kKQRwv0kHQb/S/LlejBA5pyawcDDJAXJAmA2ybk8NbkEODWz7l3874thLtSYNQsXkEw9uzIzNXAAMCCE0K2WbVoAw0hGPz4l+cv95Zll95N82Z1G0uAgl88TknPXqp4bczfJ+S8VIYRH17dxjHEKyTlHt5B86X+fTAOMamxG8l7NI5kyNpfkGlQANwDLSb6430fSRCFXNub9hKS4uC/zmvTJYv2fkLwOL2amr/6Nmpsy3AX8heQ9fhUYV8t+9wNeyvx7Gw/8MMb43zpmrPQgyefqw8zPLwFijO8Cv8hkfw9Y93yv9X1OfglMIRklfT3z3LK9vtoZwLTMa3cucHqVZacBt2e5H0mNTFjPlHxJarQyowwVJNP//ru+9VV3mVGzfwE9NrDhhFRQIYS9gDtjjOuekyVJgAWWJK0lhHAc8BTJ1MDrgAOAfdbXIEKSJAmcIihJ6zqB5CT4mcDOJNMNG0xxFULoFJKLslb3U2NjA0mSlB1HsCRJkiQpRxzBkiRJkqQc2ZirxxdM+/btY3l5edoxJEmSJDVSr7zyyucxxpou2L5GURRY5eXlTJkyJe0YkiRJkhqpEMJH2aznFEFJkiRJyhELLEmSJEnKEQssSZIkScqRojgHS5IkSSq0FStWMH36dJYuXZp2FBVQy5Yt6dixI82bN6/T9hZYkiRJUjWmT59OmzZtKC8vJ4SQdhwVQIyRuXPnMn36dLbffvs67cMpgpIkSVI1li5dSrt27SyuGpEQAu3atduoUUsLLEmSJKkGFleNz8a+5xZYkiRJkpQjFliSJEmSlCMWWJIkSVI9NG3aNPbcc8+1Hhs6dCjXXnttrdtNmTKFCy64AIBnnnmG559/vsZ1J06cyP77789uu+1Gly5d6Nu3Lx9//PHGh99Aw4cPZ+TIkTnZ16WXXsoee+zBpZdeyu23375mvyNGjGDmzJlr1jvllFN47733cnLMquwiKEmSJDUgXbt2pWvXrkBSYLVu3ZpvfvObX1nvjTfeYPDgwYwfP57OnTsDMH78eKZNm0anTp0KlnflypXcc889vPrqqznZ3x133MGcOXNo0aLFWo+PGDGCPffck2222QaA8847j2uuuYa77rorJ8etZIElSZIkrcePfgRTp+Z2n126wPDhdd/+O9/5DgcccABPP/00FRUV3H333XTr1o1nnnmGa6+9lltuuYXbb7+dpk2b8sADD3DzzTfTrVu3Ndv/5je/4fLLL19TXAEcf/zxa27fdddd3HnnnSxfvpyddtqJ+++/n1atWtG/f3+OPfZYevfuDUDr1q1ZuHAhs2bNom/fvsyfP5+VK1dy22238c1vfpOBAwcyZcoUQgicddZZXHjhhWs9j0mTJrHPPvvQrFmzWp/X0qVLOe+885gyZQrNmjXj+uuv59BDD11rX8cffzyLFi3igAMO4LLLLuPtt9+mdevWlJeXM2XKFE477TRKSkp44YUX6NatG/3792flypVrjp0LThGUJEmSitTKlSt5+eWXGT58OD//+c/XWlZeXs65557LhRdeyNSpU9cqrgDefPNN9tlnnxr3fdJJJ/HPf/6T1157jc6dO3P33XfXmuXBBx/kyCOPZOrUqbz22mt06dKFqVOnMmPGDN544w1ef/11BgwY8JXtJk+ezL777rve53XrrbcC8Prrr/PQQw/Rr1+/r7RTHz9+PCUlJUydOpW+ffuuebx379507dqVUaNGMXXqVEpKSmjSpAk77bQTr732Wq3Pa0M5giVJkiStx8aMNNVVTe3Cqz5+0kknAbDvvvsybdq0Oh9r7ty59OjRg8WLFzNo0CAuueQS3njjDa644goqKipYuHAhRx55ZK372G+//TjrrLNYsWIFvXr1okuXLuywww58+OGHDB48mJ49e3LEEUd8ZbtZs2atNYpW0/N67rnnGDx4MAC77bYb2223He+++y577713nZ/3VlttxcyZM79S4G0MR7AkSZKkeqhdu3bMmzdvrce++OIL2rdvv+Z+5XlGTZs2ZeXKlRu0/z322GPNeU/t2rVj6tSpDBo0iIULFwLQv39/brnlFl5//XWuuuqqNaNFzZo1Y/Xq1QDEGFm+fDkA3/rWt3j22WcpKyvjjDPOYOTIkWy++ea89tprfOc73+HWW2/l+9///ldylJSUfGUkqrrnFWPcoOeXjaVLl1JSUpLTfVpgSZIkSfVQ69at6dChA0899RSQFFcTJ07kkEMOyXofbdq0YcGCBdUu+/GPf8yvfvUr3n777TWPLV68eM3tBQsW0KFDB1asWMGoUaPWPF5eXs4rr7wCwGOPPcaKFSsA+Oijj9hqq604++yzGThwIK+++iqff/45q1ev5rvf/S5XX311tY0sOnfuzPvvv7/e5/Ktb31rTY53332Xjz/+mF133TWLVyFR3Wvx7rvvsscee2S9j2w4RVCSVK8cPGwSMyqWVLusrLSEyUO6FziRJKVn5MiR/OAHP+Diiy8G4KqrrmLHHXfMevvjjjuO3r1789hjj32lycVee+3FjTfeyJlnnsmCBQto164dnTp1WnPO09VXX80BBxzAdtttx1577bWmODn77LM54YQT2H///enRowebbropkHQs/O1vf0vz5s1p3bo1I0eOZMaMGQwYMGDNiNevf/3rr2Q8+uijOeOMM9b7XM4//3zOPfdc9tprL5o1a8aIESO+0imwNv379+fcc89d0+Ri/vz5lJSU0KFDh6z3kY2Qj6G2XOvatWucMmVK2jEkSQVQPmQC04b13OBlkpRrb7/99lfODVJ+nHjiiVxzzTXsvPPOBTvmDTfcwGabbcbAgQO/sqy69z6E8EqMsev69usUQUmSJEmpGjZsGLNmzSroMUtLS+nXr1/O9+sUQUmSJKkGMcYau/kpd3bdddcNOp8qF6prGQ8b30zDESxJkiSpGi1btmTu3Ll56V6n+inGyNy5c2nZsmWd9+EIliRJklSNjh07Mn36dObMmZN2FBVQy5Yt6dixY523t8CSJEmSqtG8eXO23377tGOoyDhFUJIkSZJyxAJLkiRJknLEAkuSJEmScsQCS5IkSZJyxAJLkiRJknLELoKSpAbv4GGTmFGxpMblZaUlTB7SvYCJJEkNlQWWJKnBm1GxhGnDeta4vHzIhAKmkSQ1ZE4RlCRJkqQcscCSJEmSpByxwJIkSZKkHLHAkiRJkqQcscCSJEmSpByxwJIkSZKkHLHAkiRJkqQcscCSJEmSpByxwJIkSZKkHLHAkiRJkqQcscCSJEmSpByxwJIkSZKkHLHAkiRJkqQcscCSJEmSpBxpls+dhxCmAQuAVcDKGGPXEMIWwGigHJgG9IkxzstnDkmSJEkqhEKMYB0aY+wSY+yauT8EeCrGuDPwVOa+JEmSJBW9NKYIngDcl7l9H9ArhQySJEmSlHP5LrAi8NcQwishhEGZx7aOMc4CyPzeqroNQwiDQghTQghT5syZk+eYkiRJkrTx8noOFnBwjHFmCGEr4MkQwjvZbhhjvBO4E6Br164xXwElSZIkKVfyOoIVY5yZ+T0b+COwP/BZCKEDQOb37HxmkCRJkqRCyVuBFULYNITQpvI2cATwBjAe6JdZrR/wWL4ySJIkSVIh5XOK4NbAH0MIlcd5MMY4MYTwT2BMCGEg8DFwch4zSJIkSVLB5K3AijF+CHy9msfnAj3ydVxJkiRJSku+m1xIkpQzZaUllA+ZUOOyyUO6FziRJElrs8CSJBWN2gqomgovSZIKKY0LDUuSJElSg2SBJUmSJEk5YoElSZIkSTligSVJkiRJOWKBJUlqsBYvhnffheWz2/D++8l9SZLyyS6CkqQGoay0hO1+MoFl0zdn0dtlLJ3WjpXzWmeWfoud74UQoHNnOPxwOOUUOOCA5DFJknLFAkuSVPRihB/u0J1f/AI+fhM23RSOOhT23x+23x5atYKFC+HDD+HFF+G22+DGG2G//WDo0GR7SZJywQJLklTUPvgABg2CSZOS0al774WTT06KrJrMnw+jRsF110HPntBy+/2YPhg6dixcbklSw+Q5WJKkonXPPbDXXjBlCtx+O7z+OvTvX3txBbDZZnDeefDWWzB8OCybvgV77w2PPFKQ2JKkBswCS5JUdFauhAsugIED4eCDk0LpnHOgadMN288mm8APfwgd+v+DnXaCPn3gZz9zyqAkqe4ssCRJRWXpUujVC26+GS66CJ54AsrKNm6fzbdYzOTJcNZZcPXV0K8fLF+em7ySpMbFc7AkSTl38LBJzKhYUu2ystISJg/pXqf9Ll0KJ54IEycmUwLPOWdjUq6teXP4/e+hvDwZxVqyBB56CJr5f0pJ0gbwfxuSpJybUbGEacN6VrusfMiEOu1zxYqkuPrLX+Duu5PRplwLAa68MjmH6+KLoaQERoyAJs73kCRlyQJLklTvxZg0pZg4Ee68Mz/FVVUXXZRclPjKK6FdO7jhhvweT5LUcFhgSZLqvV//Ohm1uvJKOPvswhzzpz+FOXOSLoN77AHf/35hjitJKm5OepAk1Wt/+lNS7Jx2Gvz854U7bgjJdbKOOALOPx+efbZwx5YkFS8LLElSvTVtGpx5JnzjG0kDihAKe/xmzWD0aNh+e+jbF2bPLuzxJUnFxymCkqSCKistqbXRRVlpCQDLlsHJJ8Pq1ckFgFu2LFTCtZWWwtixsN9+Sfv2CRP+1/Sitm6JsHEdEyVJxckCS5JUUNkWHEOHwpQpMG4c7LhjfjOtz157wfXXww9+kPy+5JLk8dq6JULdOyZKkoqXUwQlSfXOCy/ANdfAwIFJa/b64Lzzkiw//Sm8+WbaaSRJ9ZUFliSpXlm8OJmK17FjMlpUX4QAd9wBm20GAwbAypVpJ5Ik1UcWWJKkemXoUHjvPbj33qSYqU+23BJuvhn++U+vjSVJqp4FliSp3njjjaRw+f73oXs97Q3Rty/06pVck2vFvFZpx5Ek1TMWWJKkeiHG5HpTbdvCsGFpp6lZCHDrrdC8Ocx7ave040iS6hkLLElSvXD//fCPf8BvfgPt2qWdpnbbbJNMZVzywdY8/njaaSRJ9YkFliQpdfPmJa3PDzooaSBRDC64AJq3W8AFF8CSmi+FJUlqZCywJEmp++UvYe5c+N3v/ncR3/queXPY4og3mDYNhg9PO40kqb4okv+NSZIaqv/+F265JRm56tIl7TQbpmWnLzj++OScsc8/TzuNJKk+sMCSJKXqiiugaVP4+c/TTlI3w4bBwoXJKJwkSRZYkqTUvPIKPPggXHghlJWlnaZuOndO2sr/7nfwwQdpp5Ekpa1Z2gEkScXp4GGTmFFRfXeHstKS9W4fI/z4x9C+ffI7TWWlJZQPmVDjsslDar8o19Ch8MAD8NOfwsMP5yGgJKloWGBJkupkRsUSpg3rWeftn3wSJk2Cm25Krn2VptoKqJoKr6o6dEi6IP7iF3DppbDvvrlMJ0kqJk4RlCQVXIxw1VWw7bZwzjlpp8mNiy+GzTcv3nPJJEm5YYElSSq4J5+EF1+Eyy+HTTZJO01ubLYZXHQRPP44vPpq2mkkSWmxwJIkFVSMySjPttsWz0WFszV4sKNYktTYWWBJkgrqb3+D55+Hyy6DFi3STpNbbdsmHRHHj4d//SvtNJKkNFhgSZIKpnL0qmNHOOustNPkxwUXQGmpo1iS1FhZYEmSCuaZZ2DyZBgypOGNXlVq2xZ+9CN47DFYPrtN2nEkSQVmgSVJKphrroGttoKBA9NOkl8XXACbbgrzX94h7SiSpAKzwJIkFcS//w0TJ8IPfwgtW6adJr823xzOPhsWvbUNH32UdhpJUiFZYEmSCuK3v01Gdc47L+0khXHhhcnv4cPTzSFJKiwLLElS3n30ETz0EAwalIzuNAadOsGmu8/krrtg3ry000iSCsUCS5KUd8OHQwhJ84fGZLP9P2TRIrjttrSTSJIKpVnaASRJDdu8eXDXXXDKKcmoTmOyyVYLOOoouPFGuOiitc89O3jYJGZULKl2u7LSEiYP6V6glJKkXLLAkiTl1d13w6JFcPHFaSdJx8UXw+GHw+jR0K/f/x6fUbGEacN6VrtN+ZAJBUonSco1pwhKkvJm1Sq49Vb41regS5e006SjRw/o3Bluvjm50LIkqWGzwJIk5c2ECTBtGgwenHaS9ISQPP9XXoEXX0w7jSQp3yywJEl5c9NN0LEj9OqVdpJ0nXEGtG2bvB6SpIbNAkuSlBdvvQVPPQXnnw/NGvkZv61bw1lnwdixMHNm2mkkSflkgSVJyotbboEWLeDss9NOUj/84AfJOWl33JF2EklSPllgSZJyrqIC7rsPvvc9aN8+7TT1w447wjHHwO23w7JlaaeRJOWLBZYkKefuvRcWL27czS2qc8EFMHs2PPJI2kkkSfmS9wIrhNA0hPCvEMKfMve3DyG8FEJ4L4QwOoSwSb4zSJIKZ/Vq+N3v4OCDYZ990k5Tvxx2GOy6azJ9UpLUMBViBOuHwNtV7v8GuCHGuDMwDxhYgAySpAJ5+ml4/30477y0k9Q/TZrAuefCSy/B8tlt0o4jScqDvBZYIYSOQE/g95n7AegOjM2sch/QyJv3SlLDcuedsPnm8N3vpp2kfjrjjKT5x8LXOqUdRZKUB/kewRoO/BhYnbnfDqiIMa7M3J8OlFW3YQhhUAhhSghhypw5c/IcU5KUC7Nnwx//CP36QcuWaaepn9q1g969YeEbZSxenHYaSVKu5a3ACiEcC8yOMb5S9eFqVo3VbR9jvDPG2DXG2HXLLbfMS0ZJUm7ddx+sWGFr9vUZNAji8uaMGZN2EklSruVzBOtg4PgQwjTgYZKpgcOB0hBC5SUnOwJeclGSGoAYk+mBhxwCu++edpr6rVs3aLbFQu68M+0kkqRcy1uBFWO8LMbYMcZYDpwCTIoxngY8DfTOrNYPeCxfGSRJhfPMM0lzi0GD0k5S/4UAbbp8zAsvwOuvp51GkpRLaVwH6yfARSGE90nOybo7hQySpByrbG7Ru/f61xVsuud0NtkER7EkqYEpSIEVY3wmxnhs5vaHMcb9Y4w7xRhPjjF6PXtJKnJz5sC4cXDmmVBSknaa4tC0ZAW9e8P992OzC0lqQNIYwZIkNTAjR8Ly5Ta32FDnnANffgmPPJJ2EklSrlhgSZI2SowwYgQccADssUfaaYpLt26w007J6ydJahgssCRJG+WVV+CNN2DAgLSTFJ8QoH//pEHIhx+mnUaSlAsWWJKkjTJiBLRoAX37pp2kOJ1xRlJojRyZdhJJUi5YYEmS6mzZMnjwQTjxRCgtTTtNcerUCXr0SC7SvHp12mkkSRvLAkuSVGePPw7z5jk9cGMNGADTpsGzz6adRJK0sSywJEl1du+9UFaWjMCo7nr1gs02S15PSVJxa5Z2AElScVq5sAUTJ8JPfgJNm6adJn/KSksoHzKhxmW52LZVq+QctlGj4JZb6p5VkpQ+CyxJUp0serOM1auhX7+0k+TX5CHdC7Jt//5w110wdmydDydJqgecIihJ2mAxwqLXO3LQQbDrrmmnaRgOOgh22cVrYklSsbPAkiRtsClTYMXcNja3yKHKa2I9+yysmNcq7TiSpDqywJIkbbB774XQbBV9+qSdpGGpvCbWojc6ph1FklRHFliSpA2ybBk8/DCU7PwpbdumnaZh6dgRDj8cFr5R5jWxJKlIWWBJkjbIE08k175qvceMtKM0SKefDqvmt+L559NOIkmqCwssSdIGGTUKttwSWm7/edpRGqRevZLpl6NGpZ19smchAAAgAElEQVREklQXFliSpKx9+SU8/jiccgqEJjHtOA1SmzbJ9MsxY2D58rTTSJI2lAWWJClrf/hDcg7WaaelnaRh23SPmXzxBUycmHYSSdKGssCSJGVt1CjYcUfYf/+0kzRsJeVzaN8epwlKUhGywJIkZWXGDHj66WT0KoS00zRsoWmkb18YPx7mz087jSRpQ1hgSZKy8vDDEKPTAwvltNNg6VIYNy7tJJKkDWGBJUnKygMPwH77wS67pJ2kcTjwQNhhB6cJSlKxscCSJK3XW2/B1KmOXhVSCHDqqTBpEsycmXYaSVK2LLAkSes1ahQ0bZq0Z1fhnHYarF6dTM+UJBUHCyxJUq1Wr04KrMMOg623TjtN47LbbrDvvk4TlKRiYoElSarV88/DRx85PTAtp50Gr74K77yTdhJJUjYssCRJtRo1CkpKoFevtJM0TqecAk2aOIolScWiWdoBJEn108HDJjF97lKmjziMltvPYa9fTV1reVlpSUrJGpcOHaBHj6TA+sUvvAaZJNV3FliSpGrNqFjC7d17csy18NCwMo4/viztSI3WqafCgAHwz3/C/vunnUaSVBunCEqSajRmDLRtC0cemXaSxq1XL9hkE7sJSlIxsMCSJFUrrmzCH/+YfLlv0SLtNI1baSkcdVRS8K5enXYaSVJtLLAkSdVaMq09X34JffqknUQAffvCjBkweXLaSSRJtbHAkiRVa/E7Hdh88+T6V0rf8ccn3RydJihJ9ZsFliTpK5YuhcXvbc2JJybn/ih9rVtDz54wdiysXJl2GklSTSywJElf8Ze/QFzenL59006iqk45BWbPhmeeSTuJJKkmFliSpK8YPRqalCzn0EPTTqKqjjkmGckaPTrtJJKkmlhgSZLWsmQJjB8PrXb5lObN006jqkpK4IQT4A9/gOXL004jSaqOBZYkaS1//jMsWgStdpuZdhRVo29fmDcP/va3tJNIkqpjgSVJWsuYMbDlltCy0xdpR1E1jjgiuS6W3QQlqX5qlnYASVJ+HTxsEjMqllS7rKy0hMlDuq+5v2gR/OlP0K8f/LlJLFREbYAWLeCkk+CRR5Jujy1bpp1IklSVBZYkNXAzKpYwbVjPapeVD5mw1v0JE2Dx4uTiwn+eWIh0qou+feGee+CJJ+DEE9NOI0mqKqspgiGEPfMdRJKUvjFj4Gtfg27d0k6i2nTvDu3b201QkuqjbM/Buj2E8HII4fwQQmleE0mSUrFgQTKC1bs3NG2adhrVplmz5H16/PFkWqckqf7IqsCKMR4CnAZsC0wJITwYQjg8r8kkSQX1pz8l5/T06ZN2EmXjlFOS6ZyPP552EklSVVl3EYwxvgdcAfwE+DZwUwjhnRDCSfkKJ0kqnDFjYJtt4OCD006ibBxyCHTo4DRBSapvsmpyEULYGxgA9ASeBI6LMb4aQtgGeAEYl7+IkqR8KSstoXzIBFYva8Ynjx9Gmy4fs8Plb61ZpvqradNktPG22+DLL6Ft27QTSZIg+y6CtwB3AZfHGNf0+o0xzgwhXJGXZJKkvKts0f7AA3DGKvjLzdtz0EHbp5xK2erbF268ER59NGmtL0lKX7ZTBI8BHqwsrkIITUIIrQBijPfnK5wkqTBGj4Ztt4UDDkg7iTbEgQdCp07J9E5JUv2QbYH1N6DqXJFWmcckSUWuogL+8pdkulmTrM/MVX0QQvK+/fWv8MUXaaeRJEH2BVbLGOPCyjuZ263yE0mSVEiPPQYrVtg9sFj17QsrVybTBCVJ6cu2wFoUQtin8k4IYV9gSS3rS5KKxOjRUF4O++2XdhLVxb77wg472E1QkuqLbAusHwGPhBD+EUL4BzAa+L/8xZIkFcIXX8CTTyajVyGknUZ1UTlN8Kmn4PPP004jScqqi2CM8Z8hhN2AXYEAvBNjXJHXZJKkvPvjH5PpZU4PLC4HD5vEjIr/TSRZ/tlmrFrVjV1P/ze7fefzNd0h17fduspKS2rcVpKUnWzbtAPsB5RntvlGCIEY48i8pJIkFcSYMbDjjrDPPutfV/XHjIolTBvWc839GGHXybDtir35oGJC1tutq3xIzdtKkrKT1RTBEML9wLXAISSF1n5A1zzmkiTl2eefJ9PKnB5Y/EJIml088wysWrRJ2nEkqVHLdgSrK7B7jDFmu+MQQkvgWaBF5jhjY4xXhRC2Bx4GtgBeBc6IMS7fsNiSpI01bhysWuX0wIaiTx/45S9h8X86pB1Fkhq1bJtcvAF8bQP3vQzoHmP8OtAFOCqEcCDwG+CGGOPOwDxg4AbuV5KUA2PGwC67wNe/nnYS5cKee0LnzrDoHQssSUpTtgVWe+CtEMJfQgjjK39q2yAmKq+d1TzzE4HuwNjM4/cBveqQW5K0ET77DJ5+2umBDUnlNMFln2zBzJlpp5GkxivbKYJD67LzEEJT4BVgJ+BW4AOgIsa4MrPKdKCshm0HAYMAOnXqVJfDS5JqMG4crF6dfCFXw9GnDwwdGvjDH2Dw4LTTSFLjlNUIVozx78A0oHnm9j9Jzp9a33arYoxdgI7A/kDn6larYds7Y4xdY4xdt9xyy2xiSpKyNGZMMp1sjz3STqJc6twZmm8534sOS1KKsu0ieDbJtL47Mg+VAY9me5AYYwXwDHAgUBpCqBw56wg4kUGSCmjWLPj735PRK6cHNjytdpvF5MnwySdpJ5Gkxinbc7B+ABwMzAeIMb4HbFXbBiGELUMIpZnbJcBhwNvA00DvzGr9gMc2PLYkqa7+8Ifkukknn5x2EuXDprslf7ccO3Y9K0qS8iLbAmtZ1VbqmRGo9bVs7wA8HUL4N8mUwidjjH8CfgJcFEJ4H2gH3L3hsSVJdTV6dNJxbvfd006ifGi+xWK+8Q2cJihJKcm2ycXfQwiXAyUhhMOB84HHa9sgxvhv4BvVPP4hyflYkqQCmzEDnnsOrr467STKpz594LLLYNo0KC9PO40kNS7ZjmANAeYArwPnAH8GrshXKElSflROG/Piwg1b5fv7yCPp5pCkxijbLoKrY4x3xRhPjjH2ztxe3xRBSVI9M3o0dOmSXGBYDdcOO0DXrkm3SElSYWXbRfC/IYQP1/3JdzhJUu58/DG88IKjV41F374wZQp88EHaSSSpccl2imBXYL/MTzfgJuCBfIWSJOVe5XQxC6zGobJLpNMEJamwsp0iOLfKz4wY43Cge56zSZJyaMwY2Hdf2HHHtJOoELbbDg480G6CklRo2U4R3KfKT9cQwrlAmzxnkyTlyLRp8PLLjl41Nn37wtSp8O67aSeRpMYj2ymC11X5+TWwL+D/piWpSFQ2O7DAalx6905+2+xCkgonq+tgxRgPzXcQSVLdHDxsEjMqltS4vKy0hDFjYP/9vSZSY9OxIxxySDJN8AovriJJBZFVgRVCuKi25THG63MTR5K0oWZULGHasJ41Ln//fdj5MrjuugKGUr3Rpw9ccAG89VbaSSSpcdiQLoLnAWWZn3OB3UnOw/JcLEmqxyq7yFVOF1Pj0rs3hOA0QUkqlKxGsID2wD4xxgUAIYShwCMxxu/nK5gkKTfGjIGDDoJOndJOojR06ADf/nbyOYjHpZ1Gkhq+bEewOgHLq9xfDpTnPI0kKafefTfpIte3b9pJlKY+feDtt2HF5046kaR8y7bAuh94OYQwNIRwFfASMDJ/sSRJuTBmTDI9zOmBjdt3vwtNmsCitzukHUWSGrxsuwj+KoTwBNAt89CAGOO/8hdLkpQLo0cnXeTKytJOojRttRUceij8/dUOxJgU3ZKk/Mh2BAugFTA/xngjMD2EsH2eMkmScuCtt+CNN7z2lRJ9+8LKea2ZOjXtJJLUsGVVYGWmBf4EuCzzUHPggXyFkiRtvNGjk2lhTg8UwIknAmG13QQlKc+yHcE6ETgeWAQQY5yJ7dklqd6KMSmwvv1t+NrX0k6j+qB9e2hZ/jmjRyefD0lSfmRbYC2PMUYgAoQQNs1fJEnSxvr3v+E//7F7oNa26W6z+O9/4ZVX0k4iSQ1XtgXWmBDCHUBpCOFs4G/AXfmLJUnaGKNHQ9OmSfc4qVLJLp/SvHny+ZAk5Ue2XQSvDSEcDswHdgV+FmN8Mq/JJEl1Ujk9sEePZFqYik9ZaQnlQybUurwumrZcyRFHJO37r7nGboKSlA/rLbBCCE2Bv8QYDwMsqiSpnpsyBT78EH7607STqK4mD+met3336QMTJsBLL8GBB+btMJLUaK13imCMcRWwOITQtgB5JEkbafRoaN480zVOWscJJ8AmmzhNUJLyJaspgsBS4PUQwpNkOgkCxBgvyEsqSVKdrF6dTP868kjYfPO006g+atsWjj4aHnkErrsuaeUvScqdbP+zOgG4EngWeKXKjySpHnnxRfjkE7sHqnZ9+sCMGfD882knkaSGp9YRrBBCpxjjxzHG+woVSJJUd6NHQ4sWcPzxaSdRfXbccdCyZTLaecghaaeRpIZlfSNYj1beCCH8Ic9ZJEkbYdWqZNrXMcfAZpulnUb1WZs2yefkkUeSz40kKXfWV2BVbeC6Qz6DSJI2znPPwaxZTg9Udvr2hU8/TT43kqTcWV+BFWu4LUmqZ0aPhlat4Nhj006iYtCzZ/J5sZugJOXW+gqsr4cQ5ocQFgB7Z27PDyEsCCHML0RASdL6rVwJY8cmxdWmm6adRsVg002Tz8vYscnnR5KUG7UWWDHGpjHGzWKMbWKMzTK3K+87w1+S6olnnoE5c5weqA3Tt2/yufn739NOIkkNh1e/kKQG4OGHoXXr5PpGUraOPjr53DhNUJJyxwJLkopcXBUYNw5OOAFKStJOo2JSUpK09B83DlasSDuNJDUMFliSVOSWTmvPvHlOD1Td9OkDc+fCpElpJ5GkhsECS5KK3KJ3tqFtWzjiiLSTqBgddVRy3bQxY9JOIkkNgwWWJBWxpUth8btbc+KJ0KJF2mlUjFq0gF69kmmCcVVY/waSpFpZYElSEZs4EeLy5k4P1Ebp0wcqKmDJtC3TjiJJRc8CS5KK2KhR0KTVMg47LO0kKmaHHw6lpbD4nQ5pR5GkomeBJUlFav58ePxx2HS3WTRrlnYaFbNNNoETT0ymmy5dmnYaSSpuFliSVKT++EdYtgw23X1G2lHUAPTtm0w3/etf004iScXNAkuSitSoUbD99rDJNhVpR1ED0L07NClZ7kWHJWkjWWBJUhH69FN46ik49VQINn5TDjRvDq12+ZTx42HJkrTTSFLxssCSpCI0ZgysXg2nnZZ2EjUkrXabycKF8Oc/p51EkoqXBZYkFaFRo6BLF+jcOe0kakhadprL176WfL4kSXVjgSVJReb99+Hll5PpgVIuhSbJ52rCBPjii7TTSFJxsrGvJNUTBw+bxIyK6k9+KSstYfKQ7gA8+GBy3tX3vlfIdGosTj8drr8exo6FQYPSTiNJxccCS5LqiRkVS5g2rGe1y8qHTAAgxqTA+ta3oGPHQqZTY9GlC+y+OzzwgAWWJNWFUwQlqYj861/wn//Y3EL5E0Ly+frHP2DatLTTSFLxscCSpCIyalTSTvu73007iRqyyvP7Hnww3RySVIwssCSpSKxaBQ8/DEcfDVtskXYaNWTl5dCtG9x/fzItVZKUPQssSSoSzz4LM2c6PVCFcfrp8M47ybRUSVL2LLAkqUiMGgWtW8Oxx6adRI3BySfDJpskzS4kSdmzwJKkIhBXNmHsWDjpJGjVKu00agw23xx69oSHHoKVK9NOI0nFwwJLkorA4ve34ssvnR6owjr9dPj0U5g0Ke0kklQ8LLAkqQgseqMj22wDPXqknUSNyTHHQGlpMj1VkpSdvBVYIYRtQwhPhxDeDiG8GUL4YebxLUIIT4YQ3sv83jxfGSSpIZg9G5Z8uCWnnw5Nm6adRo1Jy5bJuVjjxsGiRWmnkaTikM8RrJXAxTHGzsCBwA9CCLsDQ4CnYow7A09l7kuSavDgg0BsQr9+aSdRY3T66bBwIYwfn3YSSSoOeSuwYoyzYoyvZm4vAN4GyoATgPsyq90H9MpXBklqCEaOhE2+VsHuu6edRI3RIYfAttvaTVCSstWsEAcJIZQD3wBeAraOMc6CpAgLIWxVwzaDgEEAnTp1KkRMSap3Xn89uQ7R9sfNpnzI5GrXKSstKXAq1WdlpSWUD5lQ47IN1aRJ0lzlt7+Fzz6Drbfe2ISS1LDlvcAKIbQG/gD8KMY4P4SQ1XYxxjuBOwG6du3qdeQlNUr33QfNmsHL9+xC+/a7pB1HRWDykO453+cZZ8CwYcl01QsvzPnuJalByWsXwRBCc5LialSMcVzm4c9CCB0yyzsAs/OZQZKK1cqVSfe2nj2hffu006gx23132H9/uPdeiP7JU5Jqlc8uggG4G3g7xnh9lUXjgcpTtfsBj+UrgyQVsyefTK5BZHML1QcDBiRTVl99Ne0kklS/5XME62DgDKB7CGFq5ucYYBhweAjhPeDwzH1J0jpGjoQttkiuRSSlrW9faNECRoxIO4kk1W95OwcrxvgcUNMJV14qU5Jq8eWX8OijMHBg8qVWStvmm8OJJybnYV17rZ9LSapJXs/BkiTVzZgxsHSp0wNVvwwYAF984TWxJKk2FliSVA+NHAm77QZdu6adRPqfHj2gY8ek2YUkqXoWWJJUz3zwATz3XDJ6leWVLaSCaNoUzjwT/vIXmDkz7TSSVD9ZYElSPXPPPcnFXc84I+0k0lf17w+rV8P996edRJLqJwssSapHVq5MurQdfTSUlaWdRvqqnXeGQw7xmliSVBMLLEmqRyZOTKZeff/7aSeRata/P/znP/DSS2knkaT6xwJLkuqRu++GrbeGnj3TTiLVrE8faNXKZheSVJ28XQdLkrS2g4dNYkbFkhqXb9m0LY8/DhdfDM2bFzCYtIHatIHeveHhh+GGG+Dwm2r/bJeVljB5SPcCJpSk9FhgSVKBzKhYwrRhNQ9NXXMNTFkFZ51VwFBSHfXvn1xOYNy49X+2y4dMKFwwSUqZUwQlqR6IMZke2K0b7Lpr2mmk9fv2t2GnneDOO9NOIkn1iwWWJNUDzz0H774LAwemnUTKTpMmcPbZ8I9/wIrPW6cdR5LqDQssSaoHfv/7/53XIhWL/v2T8wUXvLZt2lEkqd6wwJKklH35JTzyCJx6Kmy6adpppOxttRX06gWL3ujI0qVpp5Gk+sEmF5KUsocegiVLnB6o+q2mLphLmrRj9dIDGTcu+SOBJDV2FliSlKIY4Y47YO+9oWvXtNNINaupU+Dq1bDLlKTZhQWWJDlFUJJS9dJLMHUqnHcehJB2GmnDVTa7+Pvf4T//STuNJKXPAkuSUnTbbdC6NZx2WtpJpLrr3x+aNYO77ko7iSSlzwJLklIydy6MHg1nnJF0EJSK1dZbJ80uRoyAZcvSTiNJ6bLAkqSUVH4ZPe+8tJNIG2/QoOSPBuPGpZ1EktJlkwtJyqGaOq0BlJWWrLm9ejXcfjscfDDstVeh0kn506MH7LBD8rn+3vfSTiNJ6bHAkqQcqqnT2rr+9jd4/30YOjT/maRCaNIEzj0XfvxjeP11/3AgqfFyiqAkpeC226B9e+jdO+0kUu4MHAgtW8Itt6SdRJLSY4ElSQU2fTqMHw9nnQUtWqSdRsqdLbZIOmI+8ADMm5d2GklKhwWWJBXYXXclFxg+55y0k0i594MfwOLFSRMXSWqMLLAkqYBWrEgKrKOOShoCSA3NN76RNG+59dakmYskNTYWWJJUQGPHwqxZyV/5pYZq8GD44AOYODHtJJJUeBZYklRAN94IO+8MRx+ddhIpf048ETp0sNmFpMbJAkuSCuSll5KfwYOTltZSQ7XJJsk5hk88kVyOQJIaE/8XL0kFcuONsNlm0L9/2kmk/Bs0CJo1g9/9Lu0kklRYFliSVAAzZsAjjySt2du0STuNlH8dOiTXebvnHli9vGnacSSpYCywJKkAbrsNVq1KpgdKjcUPfwhffgkLX9827SiSVDAWWJKUZ0uXwh13wHHH2ZpdjcuBB8JBB8GCKeWsWpV2GkkqDAssScqzBx+Ezz9P/povNTYXXwwrKzbl0UfTTiJJhWGBJUl5FCNcfz3stRccemjaaaTC69ULmpUu4vrr004iSYVhgSVJefTEE/Dmm3DppRBC2mmkwmvaFNp0ncbzz8OLL6adRpLyzwJLkvLommtg223hlFPSTiKlp/Ven1BaCtddl3YSSco/CyxJypOXX4a//x1+9CNo3jztNFJ6mmyyinPPhXHj4L//TTuNJOWXBZYk5clvfwtt28LZZ6edRErf4MHJdMHhw9NOIkn51SztAJJUbA4eNokZFUuqXVZWWgLA++/DH/4AQ4Z4YWEJYJtt4NRT4fe/hyuvhPbt004kSflhgSVJG2hGxRKmDetZ6zrXXZdMC7zgggKFkorAT34CI0fCTTfBL36RdhpJyg+nCEpSjs2eDSNGwJlnwte+lnYaqf7o3BlOOgluvhnmz087jSTlhwWWJOXY8OGwbBlccknaSaT657LLoKICbrst7SSSlB8WWJKUQ198AbfcAiefDLvumnYaqf7Zd1848sjkAtxLqj+VUZKKmgWWJOXQjTfCggVwxRVpJ5Hqr8svT6bS3nNP2kkkKfdsciFJ66itSyD8r1Pgur78MimwevWCvfbKVzqp+JSVllA+ZMKa+zFCi7KD+OHlLblm2jOEprHG7SYP6V6omJKUExZYkrSObLoEVueWW5Ii68or8xBKKmLVFUl//jb07AlDdz+GAQOq365qUSZJxcIpgpKUAwsXwg03JF8Y99kn7TRS/Xf00cm/lV/+ElasSDuNJOWOBZYk5cBtt8HcuY5eSdkKIbkW1ocfwn33pZ1GknLHAkuSNtKiRXDttXD44XDAAWmnkYrHMcfA/vvD1VfD8uVpp5Gk3LDAkqSNdNNNSUe0oUPTTiIVl8pRrI8/tqOgpIbDAkuSNkJFBVxzDRx7LHzzm2mnkYrPEUck/3Z+9StYujTtNJK08ewiKEkb4dprkyLrl79MO4m08dZtp77usnyoHMU67DC46y4YPDgvh5GkgrHAkqQ6+uwzGD4cTjkFvv71tNNIGy+ta0517w7f/jb8v/8HAwdCq1apxJCknHCKoCTV0a9/nUxp+vnP004iFbcQkkYXn36aXKxbkopZ3gqsEMI9IYTZIYQ3qjy2RQjhyRDCe5nfm+fr+JKUTx9/nLRm798fdtkl7TRS8evWDY47DoYNg88/TzuNJNVdPkewRgBHrfPYEOCpGOPOwFOZ+5JUdCpHrX72s3RzSA3JsGHJRbs9p1FSMctbgRVjfBb4Yp2HTwAqLyd4H9ArX8eXpHx57TW49174v/+DTp3STiM1HLvvnpyD9bvfJRcglqRiVOhzsLaOMc4CyPzeqqYVQwiDQghTQghT5syZU7CAklSbGOHii2GLLeDKK9NOIzU8Q4dCs2bw05+mnUSS6qbeNrmIMd4ZY+waY+y65ZZbph1HkgD485/hqafgqqugtDTtNFLDs802yR8xHn4Yls1qm3YcSdpghS6wPgshdADI/J5d4ONLUp2tWAGXXJI0tTj33LTTSA3XpZdC+/Yw7+nOxJh2GknaMIUusMYD/TK3+wGPFfj4klRnd94J77yTXFy4efO000gN12abJW3bl33SjkceSTuNJG2YfLZpfwh4Adg1hDA9hDAQGAYcHkJ4Dzg8c1+S6r2KiuTckO7d4dhj004jNXxnnw3Nt/qSSy6BRYvSTiNJ2ctnF8HvxRg7xBibxxg7xhjvjjHOjTH2iDHunPm9bpdBSaqXrroK5s6F665LLooqKb+aNoUtDn+TTz5J2rdLUrFolnYAScqXg4dNYkbFkmqXlZWWMHlI96z2M3Uq3HILnHcedOmSy4SSatOy4zxOPRV++1sYMAB22CHtRJK0fhZYkhqsGRVLmDasZ7XLyodMyGofq1fD+edDu3Ze/FRKwzXXwGOPwUUXwaOPpp1Gktav3rZpl6T6YMQIeOGF5C/om2+edhqp8SkrgyuuSIqsiRPTTiNJ62eBJUk1mDsXfvxjOOQQOPPMtNNIjdeFF8KuuyajyTa8kFTfWWD9//buPV6qst7j+OfH5rK3IAKKoUChiBdMkwTvHhXSI95TU9MKMy+ZaHYyQQXLXmqmvvJwSivz1lHKS0Le8p6aWSJeUlQESVGUTEESPSKw4Tl/PJsQ3CDgzF57Zj7v12u9ZmbNmuG7Wa9ZM7/1POt5JGkFzjwzjx542WUObCEVqUMH+OUv4eWX4Zxzik4jSStngSVJzXj44Tzv1SmnwFZbFZ1G0m67wbHHwk9+Ak89VXQaSVoxB7mQpOUsXtiGY46BjTbKk51KKkbPLg3LDEizqGNbUv1u7DB0Pj2+9gi9utWv8migktRSLLAkaTnv/HlT5k6D+++Hjh2LTiPVruaKp5t2hsMOq+fbPfbhf95ctdFAJakl2UVQkj5kwgSYO3FjTjgBBntiXGp1Dj0U9t8fRo+GhXPWKjqOJH2EBZYkNZk/P09mWtfpAy68sOg0kpoTkQeead8eZt++DY2NRSeSpGVZYElSk3POgcmTYd3/nETnzkWnkbQivXrBpZfC/JldPRkiqdWxwJIk4KGH4IIL4JhjoKHvW0XHkfQxvvxlWGvzmXz/+44qKKl1scCSVPPmzIGvfAU22QTGjCk6jaRVEQHd9nqW7t3z5/eDD4pOJEmZBZakmpYSnHACvPEGjB0LnToVnUjSqqprWMiVV8Lzz+eJwSWpNbDAklTTrrkGbropz3c1aFDRaSStrqFD4aST4JJL4Lbbik4jSRZYkmrY1Klw8smw++7wve8VnUbSmrr4YhgwAIYNg+nTi04jqdZZYEmqSYsX1HHwwdDQANdeC3V1RSeStKbq63NL9KJFcPjhsGBB0Ykk1TILLEk1JyWYfefWTJ4M11+fh3yWVNn69oWrroLHHoPTTy86jaRaZoElqeaMGQPvv/rvwlIAAA+kSURBVLAh550HQ4YUnUZSqRxyCJxySv6M33RT0Wkk1aq2RQeQpJb08MNw2mnQ0O8NRozoUXQcSZ9Azy4N9Bl5xzLrUvugw4Y7csRRndlkkzoGDCgonKSaZYElqWa88gp86Uu5K9G8fZ4mwgJLqmSPjBzc7Po3vgu9Np3HAQc0MHEi9PCjLqkF2UVQUk145x3Yd988Gen48dCmQ2PRkSSVSY8esP4hj/P223DQQU5CLKllWWBJqnoLF8Khh8KUKXDzzdC/f9GJJJVb+0/N5dprYcIEOO64PLiNJLUECyxJVS0lOPFEuO8+uPxyB7WQasnBB+dJxK+7DkaNKjqNpFrhNViSqtr558OVV8JZZ8HXv150Gkkt7ayzYMaMfCzo3h1OPbXoRJKqnQWWpKo194k+jLoPjjoKfvjDotNIKkIEXHYZzJoF3/lOLrKOOqroVJKqmV0EJVWlq6+GOfdtyYEH5vttPNpJNauuDsaOhd13h6OPhjvvLDqRpGrmTw5JVeemm+DYY6G+z1vccAO0a1d0IklFq6+HW26BrbaCL34R7r676ESSqpUFlqSqMn48HHkk7LQTdD/4cTp0KDqRpNaic2e4917YYgs48EBbsiSVhwWWpKpx7bV5IuFBg+D226FNu8VFR5LUyqy7Ltx/f56u4aCD4I47ik4kqdpYYEmqCj//OXzta7DbbnDPPbDOOkUnktRadeuWi6wl3QXHjSs6kaRqYoElqaKlBBdeCN/6Fuy/fz4b3alT0akktXZdu+b58QYOzBORX3pp0YkkVQsLLEkVq7ERTjoJRoyAww+Hm2/OF7JL0qro0iUXWfvvD8OHw5ln5pM2kvRJWGBJqkhz58J+++WugaefDr/5jaMFSlp9a62VT86ccAL86EcwbBjMn190KkmVzImGJVWc6dPzGecXXoBf/SoPyS5Ja6pt23yyplcvGD0apk7NRVfPnkUnk1SJbMGSVFH+8AfYdluYMQPuusviSlJpRMCoUbmweu65fJx5+OGiU0mqRBZYkipCY2O+PmLffaF3b5g4EYYMKTqVpGpz8MEwYUKeM2vwYBgzxuuyJK0eCyxJrd7rr8Oee+brI449Fv76V+jXr+hUkqpV//7w2GMwdCicemq+3vOf/yw6laRK4TVYklqtlGDsWDj55HzR+TXX5AvQS6Fnlwb6jGx+htGeXRpK849IKszKPuOr8tpHRg7mllvy8O2nnZbnzLr66tyKLkkrY4ElqVV680345jdh/HjYaadcXJWy1eqRkYNL92aSWp1P8hlfUphF5OHb99gDjjwyt2QdcwxcdFGerFiSmmMXQUmtyuLFeWTA/v3zgBYXXQR/+pNdAiUVZ8st83VZI0bAr38NW2yRp4bw2ixJzbHAktRqPP447LgjHH98/kHz5JO5a05dXdHJJNW6+nq44AJ44gno0weOOgr23huef77oZJJaGwssSYWbMQO+8Q3Ybjt45RW47jp48MHciiVJrcnnPgd/+Qv89Ke5VWvrreHEE3O3ZkkCCyxJBZo9O7dQ9euXi6pTT4UpU/KZ4Yii00lS8+rq8rVZ06bl4uqKK2CTTeDcc+Gdd4pOJ6loDnIhqcW98UaeW+ayy+C992DdAf+g/aDJjGs/j3E/WnbbJaN5NWfnC/7I6/+at8J/x9EAJZXTeuvllqzhw2HkSBg9Gr5/7kLWHvgyaw98mbr6xmW2X9nxTFL1sMCS1GKmTYOLL84jAi5cCIccAj/4Aezzv08y/YLmxz5e2TDLr/9r3gpfJ0ktZbPN8oinGxz9MDvO3ZXx4zclTdqU446Dk06CjTbK263psPGSKotdBCWV1cKFMG4c7LVX7gp4zTVw9NG5K+CNN3qdlaTq0aHHXMaNg6efhn32yS31ffvCQQfB/fc76qBUKyywJJXFlCkwahR8+tO5peqFF+Ccc2D6dPjFL/L1CpJUjbbeGn7723y8O+MM+POf4QtfgNd/sQdnn51b8yVVLwssSSXz8st5GONttoHNN4fzz4dtt4XbbsvPnX029OhRdEpJahk9e8J55+WRUq+7Dtp1+z/OPTe35u+yS75+a8aMolNKKjULLElrbNGiPFzxqFG5kNp443y2tqEBLrkk/3C4/XbYbz/nspJUuxoa8uionzr8MWbMyCei5syBU07JrfwDB+YRCJ9+Ok+2LqmyOciFpFW2eDE891zu7vLQQ3DPPflHQps2sNNO8OMfw2GH5Uk4JUkf1bMnjBiRlylT4JZb8gAZo0fnpXt3GDw4dykcMiQfT522QqosFliS/m35Yc8XzWvHgn92ZsEb6zD/tW4snNmNxnntANhgAzjggHwh9557Qteu5cnUs0vDCkfechh2SeWwsuPOx71ujd53Nxi4xzoM32wX7rsvD4hxww35qR49YPvtYYcd8u2gQdCp02pHk9SCLLAk8f778OKLMPXRdRg+YDBPPQVPPQWvvbp0m003hZn9XuVn//VpdtkldwdsibOqzhkjqaWV67izsvftM/IOhg2DYcPyaIOTJ8MDD8Cjj8KECbmlC3KPgb59Yautli6f/Ww+JrdrV5bYklaTBZZUAxobYebMfE3UjBnw6qt5mTo1d1F59d+F1Lacc2supnbaKc/fMmBAXtZbD/qMnMSwYZ8u8k+RpKoXkaew6N8/H4cBZs+Gxx7LyzPPwLPPwu9/v/Sarbo6+MxncvH14aVXr9wtcf31oa2/+qQWUchHLSL2BsYAdcAVKaULisghVaKUcovT3Lnw7rvw9tvw1lt5mTXro/dnzszL8hdOr7NOHslq111zQbXZZvCdux7mxZ/tSseOxfxtkqTmrbsuDB2alyXmzYPnn8/F1osvwt//npcbb8zfDR/Wpk3ubrjhhnnZYIP8nt26NX/btastYtKaavECKyLqgEuBPYHXgIkRcWtK6fmWzqIVW34yxOYmR/y4bdbkNUW+74q2WbQoL42NS+8v/3hlzy15vGABzJ8PH3yQl+buf3jdvHm5gHr33aXF1JJlZaNM1dfni6S7d8+tTltskUep6t172aVz54++dsRTcy2uJKlCNDTkEVy33fajz82ZAy+9BK+/nk+yLbmdOTPPz/Xoo7kIa2xc8ft36JCv91p+6dhx6f2GBmjfPm+7/NLc+rq6pUvbtss+XtXnIpYusOzjla1f0Tqp1IpowdoOmJZSegkgIq4HDgQqosC69VY48shl17W2gmB1X6PitW2bC6P6+vwFtOR+5865palXr3x/7bU/etutWy6klhRVa63ll4Yk1bquXVdcfC2RUj5pN3t2LrZmz156f84ceO+95pdZs/Ltu+8uPTk4f37L/W3lUHQxVs73r+Tsd92V54yrNJFa+Nd2RBwK7J1SOrbp8VeB7VNKw5fb7njg+KaHmwFTWjToyq0HzCo6hMrO/Vwb3M+1wf1cG9zPtcH9XBta437+TEqp+8dtVEQLVnN17keqvJTS5cDl5Y+z+iLi8ZTSwKJzqLzcz7XB/Vwb3M+1wf1cG9zPtaGS93ObAv7N14DeH3rcC5hZQA5JkiRJKqkiCqyJQL+I2Cgi2gNHALcWkEOSJEmSSqrFuwimlBojYjhwN3mY9qtSSs+1dI5PqFV2XVTJuZ9rg/u5Nrifa4P7uTa4n2tDxe7nFh/kQpIkSZKqVRFdBCVJkiSpKllgSZIkSVKJWGB9AhFxckRMiYjnIuLCovOofCLitIhIEbFe0VlUehFxUUS8EBHPRMT4iOhSdCaVTkTs3XSsnhYRI4vOo9KLiN4R8UBETG76Tv520ZlUPhFRFxFPRcTtRWdReUREl4j4XdN38+SI2LHoTKvDAmsNRcQewIHA1imlLYGLC46kMomI3sCewKtFZ1HZ3At8NqW0NTAVOKPgPCqRiKgDLgWGAv2BL0dE/2JTqQwage+mlLYAdgBOcj9XtW8Dk4sOobIaA9yVUtoc+BwVtr8tsNbcicAFKaX5ACmlNwvOo/K5BDidZibEVnVIKd2TUmpsevgoeX4+VYftgGkppZdSSguA68knx1RFUkr/SCk92XT/XfKPsZ7FplI5REQvYF/giqKzqDwiojPwH8CVACmlBSmlfxWbavVYYK25TYFdI2JCRDwUEYOKDqTSi4gDgNdTSk8XnUUt5hjgzqJDqGR6AjM+9Pg1/OFd1SKiDzAAmFBsEpXJf5NPei4uOojKZmPgLeDqpq6gV0REx6JDrY4WnwerkkTEfUCPZp46i/x/15XcFWEQcGNEbJwc977ifMx+PhPYq2UTqRxWtp9TSrc0bXMWuavR2JbMprKKZtZ5nK5SEdEJuBk4NaU0t+g8Kq2I2A94M6X0RETsXnQelU1b4PPAySmlCRExBhgJjC421qqzwFqJlNIXVvRcRJwIjGsqqB6LiMXAeuSKWxVkRfs5IrYCNgKejgjI3caejIjtUkpvtGBElcDKPs8AETEM2A8Y4omSqvIa0PtDj3sBMwvKojKKiHbk4mpsSmlc0XlUFjsDB0TEPkA90DkirkspfaXgXCqt14DXUkpLWqF/Ry6wKoZdBNfc74HBABGxKdAemFVoIpVUSmlSSmn9lFKflFIf8gf+8xZX1Sci9gZGAAeklN4vOo9KaiLQLyI2ioj2wBHArQVnUolFPgt2JTA5pfSTovOoPFJKZ6SUejV9Jx8B/NHiqvo0/c6aERGbNa0aAjxfYKTVZgvWmrsKuCoingUWAMM86y1VrJ8BHYB7m1orH00pfbPYSCqFlFJjRAwH7gbqgKtSSs8VHEultzPwVWBSRPytad2ZKaU/FJhJ0po7GRjbdGLsJeDrBedZLWFNIEmSJEmlYRdBSZIkSSoRCyxJkiRJKhELLEmSJEkqEQssSZIkSSoRCyxJkiRJKhELLElSRYuIsyLiuYh4JiL+FhHbF51JklS7nAdLklSxImJHYD/yJODzI2I98sTvkiQVwhYsSVIl2wCYlVKaD5BSmpVSmhkR05uKLSJiYEQ82HT/BxHx64i4p2mbgyPiwoiYFBF3RUS74v4USVI1sMCSJFWye4DeETE1Ii6LiN1W4TV9gX2BA4HrgAdSSlsB85rWS5K0xiywJEkVK6X0HrAtcDzwFnBDRBz9MS+7M6W0EJgE1AF3Na2fBPQpT1JJUq3wGixJUkVLKS0CHgQejIhJwDCgkaUnEeuXe8mS7oSLI2JhSik1rV+M34uSpE/IFixJUsWKiM0iot+HVm0DvAJMJ7dsARzS0rkkSbXLM3WSpErWCfhpRHQht1pNI3cX3AK4MiLOBCYUmE+SVGNiac8ISZIkSdInYRdBSZIkSSoRCyxJkiRJKhELLEmSJEkqEQssSZIkSSoRCyxJkiRJKhELLEmSJEkqEQssSZIkSSqR/we3KNJNhOi1NAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x432 with 1 Axes>"
      ]
     },
     "execution_count": 14,
     "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, 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.7.3"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}