Week4/original/HypothesisTesting/HypothesisTesting_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hypothesis Testing\n",
    "\n",
    "Python notebook for illustrating the concept of Hypothesis Testing and specific test statistics; among them the very useful Kolmogorov-Smirnov test.\n",
    "\n",
    "The Kolmogorov-Smirnov test (KS-test) is a general test to evaluate if two distributions in 1D are the same. This program applies an unbinned KS test, and compares it to a $\\chi^2$-test and a simple comparison of means. The distributions compared are two unit Gaussians, where one is then modified by changing:\n",
    "- Mean\n",
    "- Width\n",
    "- Normalisation\n",
    "\n",
    "The sensitivity of each test is then considered for each of these changes.\n",
    "\n",
    "### References:\n",
    "- Barlow: p. 155-156\n",
    "- __[Wikipedia: Kolmogorov-Smirnov test](http://en.wikipedia.org/wiki/Kolmogorov-Smirnov_test)__\n",
    "- Though influenced by biostatistics, a good discussion of p-values and their distribution can be found here:\n",
    "  [How to interpret a p-value histogram?](http://varianceexplained.org/statistics/interpreting-pvalue-histogram/)\n",
    "\n",
    "### Authors: \n",
    "Troels C. Petersen (Niels Bohr Institute)\n",
    "\n",
    "### Date:    \n",
    "07-12-2021 (latest update)\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy.special import erfc\n",
    "from scipy import stats\n",
    "\n",
    "sys.path.append('../../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression, BinnedLH, UnbinnedLH\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the parameters of the plot:"
   ]
  },
  {
   "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)\n",
    "\n",
    "save_plots = False\n",
    "verbose = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The small function below is just a simple helper function that takes a 1D-array input along with an axis, position and color arguments an plot the number of entries, the mean and the standard deviation on the axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ax_text(x, ax, posx, posy, color='k'):\n",
    "    \n",
    "    d = {'Entries': len(x), \n",
    "         'Mean': x.mean(),\n",
    "         'STD': x.std(ddof=1),\n",
    "        }\n",
    "    \n",
    "    add_text_to_ax(posx, posy, nice_string_output(d), ax, fontsize=12, color=color)\n",
    "    return None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and finally a function that calculates the mean, standard deviation and the standard deviation (i.e. uncertainty) on mean (sdom):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mean_std_sdom(x):\n",
    "    std = np.std(x, ddof=1)\n",
    "    return np.mean(x), std, std / np.sqrt(len(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the experiment:\n",
    "\n",
    "How many experiments, and how many events in each:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_exp = 1\n",
    "N_events_A = 100\n",
    "N_events_B = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define the two Gaussians to be generated (no difference to begin with!):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "dist_mean_A  =  0.0\n",
    "dist_width_A =  1.0\n",
    "dist_mean_B  =  0.0\n",
    "dist_width_B =  1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "Define the number of bins and the range, initialize empty arrays to store the results in and make an empty figure (to be filled in later):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "N_bins = 100\n",
    "xmin, xmax = -5.0, 5.0\n",
    "\n",
    "all_p_mean = np.zeros(N_exp)\n",
    "all_p_chi2 = np.zeros(N_exp)\n",
    "all_p_ks   = np.zeros(N_exp)\n",
    "\n",
    "# Figure for the two distributions, A and B, in the first experiment:\n",
    "fig1, ax1 = plt.subplots(figsize=(10, 6))\n",
    "plt.close(fig1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Loop over how many times we want to run the experiment, and for each calculate the p-value of the two distributions coming from the same underlying PDF (put in calculations yourself):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   0:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.70206\n",
      "   1:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.70206\n",
      "   2:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.70206\n",
      "   3:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.70206\n",
      "   4:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.36819\n",
      "   5:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.28194\n",
      "   6:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.90841\n",
      "   7:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.03638\n",
      "   8:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.90841\n",
      "   9:  p_mean: 0.50000   p_chi2: 0.50000   p_ks: 0.28194\n",
      "Got to experiment number: 1000\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de5zWdZ3//8cLyBARBc/KYUhDPKXirC6mK6mZechDeUzNNNk1XZW1bx62n4fbtq62rZTr2lLSlm2e8FRtma5rmqfKwfAAgmKCjIQiiYCGILx/f7yvaT7M8RpmrrmGmcf9drtuc12f4+vzuRhuz+s97+v9jpQSkiRJkrJ+1S5AkiRJ6kkMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCX1WhExMyImVLuOnigitomIX0fE8oj4t2rXUxQRNRGRImJAtWuR1DcZkCVtkCJiXkQc2mTZmRHxeMPrlNJuKaVH2jlOXw1jE4G3gCEppYtb2ygirirdn327r7S2ld77P0fEioh4OyJ+HhEjql2XpN7DgCxJFdSDg/coYFZqY7aoiAjgdOBPwBe6q7AyHZ1SGgxsB7wB/HuV65HUixiQJfVaxVbmiNg3IuoiYllEvBER15c2+3Xp59JSi+T4iOgXEV+LiPkR8WZE3BIRmxWOe0Zp3ZKI+P+anOeqiLgrIv47IpYBZ5bO/VRELI2IP0bEjRGxUeF4KSK+HBEvl7o8/FNE7FjaZ1lE3NmwfURsGRH/UzrWnyLisYho8f/yiNg/Ip6OiHdKP/cvLf8BOfB+tXTNh7a0P3AgsD1wIXByseYWzlXONf5d6Rrfjoj/KAVwIqJ/RHwzIt6KiD8AR7Z2nqZSSiuBu4Bdy91HktpjQJbUV3wb+HZKaQiwI3BnafnflH5unlIanFJ6Cjiz9PgE8BFgMHAjQETsCtwEfJ7cerkZsEOTcx1DDm2bAz8G1gCTgC2B8cAhwJeb7HM4sA/w18BXge+WzjEC2B04pbTdxUA9sBWwDXA50KwVOCKGAT8HbgC2AK4Hfh4RW6SUzizV9Y3SNT/Uyj37AvAz4I7S66Na2Y4yr/Eo4K+APYETgU+Vlp9TWrc3UAt8ro3zrCMiBgEnAb8pdx9Jao8BWdKG7L5Si+XSiFhKDq6tWQ3sFBFbppRWpJTaClSfB65PKf0hpbQCuIzcgjqAHN5+llJ6PKW0CriC5gH1qZTSfSmltSmlP6eUpqeUfpNS+iClNA+YAhzUZJ/rUkrLUkozgReAB0vnfwe4nxweG65jO2BUSml1SumxVrpJHAm8nFL6Uem8twGzgaPbuO6/KAXPE4BbU0qryYG/1W4WZV7jtSmlpSml14BfAXuVlp8IfCultCCl9CfgX8oo8b7Se74M+CTwr+VclySVw4AsaUN2bEpp84YHzVssi84GxgCzS90N2moN3R6YX3g9HxhAbrHdHljQsCKl9B6wpMn+C4ovImJMqVvEolK3i2vILa1FbxSe/7mF14NLz/8VmAs8GBF/iIhLy7yGhuto2trdmuOAD4BflF7/GPh0RGzV0sZlXuOiwvP3aLymde5pC3W35NjSe/5h4Hzg0YjYtoz9JKldBmRJfUJK6eWU0inA1sB1wF0RsQktdE8AFpK/xNZgJDksvgH8ERjesCIiNiZ3YVjndE1ef4fcevvRUhePy4FYz+tYnlK6OKX0EXJr8D9ExCFlXEPDdbxe5qm+QA6wr0XEImAa8CEau3o01Zlr/CO5K0mxzrKklNaklO4hd/E4oNz9JKktBmRJfUJEnBYRW6WU1gJLS4vXAIuBteS+xg1uAyZFxOiIGExuDb0jpfQBuavB0aUvwG0EXE37QXBTcleAFRExFji3E9dxVETsVPqC27LSNaxpYdNfAGMi4tSIGBARJ5G/yPY/ZZxjB3If4qPI3SD2Ivcbvo7Wu1l05hrvBC6IiOERMRRorVW8pVojIo4BhgIvduCcktQqA7KkvuJwYGZErCB/Ye/klNLKUheJfwaeKPVl/mvg+8CPyCNcvAqsBP4eoNRH+O+B28ktn8uBN4H32zj3V4BTS9t+j8Yvva2PjwIPASuAp4CbWhrrOaW0hBxwLyZ3AfkqcFRK6a0yznE6MCOl9GBKaVHDg/yFv49FxO4t7NOZa/we8ADwLPAMcE8Z+/ys9F4uI79/Xyi9N5LUadHGEJiSpHaUWpiXkrsWvFrteiRJnWcLsiR1UEQcHRGDSn2Yvwk8D8yrblWSpK5iQJakjjuG/CW4heQuDye3NSOdJGnDYhcLSZIkqcAWZEmSJKlgQLULKNpyyy1TTU1NtcuQJElSHzB9+vS3UkrNJkDqUQG5pqaGurq6apchSZKkPiAiWpy50y4WkiRJUoEBWZIkSSowIEuSJEkFPaoPsiRJkjYcq1evpr6+npUrV1a7lDYNHDiQ4cOH86EPfais7Q3IkiRJWi/19fVsuumm1NTUEBHVLqdFKSWWLFlCfX09o0ePLmsfu1hIkiRpvaxcuZItttiix4ZjgIhgiy226FArtwFZkiRJ660nh+MGHa3RgCxJkiQVGJAlSZLUJWpqIKLrHuVOsHzvvfcSEcyePbtLrsOALEmSpC4xfz6k1HWP+S3Oc9fcbbfdxgEHHMDtt9/eJddhQJYkSdIGa8WKFTzxxBNMnTrVgCxJkiTdd999HH744YwZM4Zhw4bxzDPPdPqYBmRJkiRtsG677TZOPvlkAE4++WRuu+22Th/TiUIkSZK0QVqyZAkPP/wwL7zwAhHBmjVriAi+8Y1vdGr4OVuQJUmStEG66667OOOMM5g/fz7z5s1jwYIFjB49mscff7xTxzUgS5IkqUuMGtW1w7yNGtX2+W677TaOO+64dZZ99rOf5dZbb+3UddjFQpL6mpqaFsdOqu8/ihFr5jVbPmoUzGu+WJKa6e7/Kx555JFmyy644IJOH9eALEl9TcNApU0Mj2hpMRvALLKS1KXsYiFJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEldo6amawdCrqlp95T9+/dnr732Ys8992TcuHE8+eSTnb4Mh3mTJElS12hlGMn1VsY4kxtvvDEzZswA4IEHHuCyyy7j0Ucf7dRpbUGWJElSr7Bs2TKGDh3a6ePYgixJkqQN1p///Gf22msvVq5cyR//+EcefvjhTh/TgCxJkqQNVrGLxVNPPcUZZ5zBCy+8QHRiGlC7WEiSJKlXGD9+PG+99RaLFy/u1HEMyJIkSeoVZs+ezZo1a9hiiy06dRy7WEiSJKlrjBpV1sgTHTpeOxr6IAOklPjhD39I//79O3VaA7IkSZK6xrx53X7KNWvWdPkx7WIhSZIkFRiQJUmSpAIDsiRJktZb6sqZ8yqkozUakCVJkrReBg4cyJIlS3p0SE4psWTJEgYOHFj2Pn5JT5IkSetl+PDh1NfXd3rc4UobOHAgw4cPL3t7A7IkSZLWy4c+9CFGjx5d7TK6nF0sJEmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqqGhAjohJETEzIl6IiNsiYmAlzydJkiR1VsUCckTsAFwA1KaUdgf6AydX6nySJElSV6h0F4sBwMYRMQAYBCys8PkkSZKkTqlYQE4pvQ58E3gN+CPwTkrpwabbRcTEiKiLiLrFixdXqpweZbfd4JFHql2FJEmSWlLJLhZDgWOA0cD2wCYRcVrT7VJK300p1aaUarfaaqtKldNMTQ1svDEMHtz4eOih8vePgLlz1+/cM2fChAnrt29f9vWv5w8X/frBD36w7rr6+nxPBw2CcePghRca161eDWefDZtuCiNHwp13dmfVkiRpQ1PJLhaHAq+mlBanlFYD9wD7V/B8Hfazn8GKFY2PQw+tdkVqy0c+AjfcAPvs03zdxImw++6wZAmcdFJ+NJg8GZ5/HhYsgFtugbPOys8lSZJaUsmA/Brw1xExKCICOAR4sYLn6zJnnplbHA84ILc6nnVW47pPfzq3NgPsuWd+ftFFjevnzcuty9OmwahRef/rrsvrLrkkb9+vX/PW6jVr4Oqrc8v2NtvABRfkls8G11wD220HQ4ZAbS28/noFLryHO/VUOOQQ+PCH112+bBk8+CBcemn+q8CkSTB/Pjz3XF4/bVq+n5tvnluZx4+He+/t9vIlSdIGopJ9kH8L3AU8AzxfOtd3K3W+rvbgg3Drrbk7xLRp8Lvf5eX3359bmwGefTY//9a3mu9/550wfTq88QYcdlhedt11efuRI5tvf/318JOfwBNP5K4bs2bl1lKAOXNy94Inn4SlS+Gmm5qHxL5s7lwYODB/+DjwwNzdYscd832D/HPsWDjtNLjnHth118Z1kiRJTVV0FIuU0pUppbEppd1TSqenlN6v5Pk66thjc6tiw+PNNxvXHX10DrIjR8Iee8BLL3Xs2FdcAVtumfvE7r13+9vffDNceSXssENudT7vvBzmILdIr1kDs2fnn/vum4+t7N13czhevhxefhnefjvfw4YPMg3r58zJLe/FdZIkSU0NqHYB1XTffa33Ox42rPH5RhvBypUdO/ZHP9qx7RcsgNNPz90vANauhW23zc/HjIEpU3Ir8oknwuGHw9SpubuFYJNNcuAdMQIWLcrLli9v7ArTsP7pp/PrCy9sXCdJktSUU02vp4i21w/o4EePESPggQdyF4qlS3O/2mKr9Zln5u4Xr7ySl3//+x0uudfaaaf8Aaa+Pr9etSrfp513zq/HjMmt7w1mzWpcJ0mS1JQBeT1tu+26Q4l11lln5W4ZCxdCSjkEP1gaNfqVV+Dhh3PwGzAgty73xdbj1atzEE6p8XnDvTjsMLj22rxs8uTGrjGQW91vuCF/8Hj0UXjqKTjuuOpeiyRJ6rn6dEA++uh1x0H+4Q/L3/eaa+D883Of4csua3/7NWsaz/Paa43nnjo1r7/44vwFswMOyIHv+OMb+0SvWpVHaNhii9waOn587o7R15xzTh6l4skn87BuG28Mv/51XjdlSh7KbehQuP12uOOOxlb+SZPyEHAjRuQv6k2dmp9LkiS1JFJK1a7hL2pra1NdXV21y5Ck3i0i/ymmzOWtbS5JG7qImJ5Sqm26vE+3IEuSJElNGZAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVJBnwzIV12Vv5V9//359dKleba8CROqWVXfsHo1nH12nu555Ei4887y9/3612G33fJsgz/4wbrr6uvz+zdoEIwb1/IY1W+/DVttlYd6kyRJak2fDMiQp4K+6678/Kc/hZqaqpbTZ0yenMcrXrAAbrklT5CyYEF5+37kI3nCj332ab5u4sQ81vGSJXDSSfnR1OWXw+jRnatfkiT1fn02IO+7L/z+9/DBB3D33fDZzzauW7MGrr46h+ZttoELLsgtnwCvvgoHHwzDhsF22+UpoFesyOvmzcst05Mnw5Zb5jD2299284X1cNOm5fu5+ea5xXf8eLj33vL2PfVUOOQQ+PCH112+bFmedfDSS/PkIZMmwfz58NxzjdtMn57fuyOO6LJLkSRJvVSfDcgABx2UW4/feAN23LFx+fXXw09+Ak88AXPnwqxZueUS4P3384xur70Gc+bkn1dcse5xly2DRYvgmGNy0FajOXNg7NjczeGee2DXXfOyzpg7FwYOzDMTHnhg7m6x446Nx00ph/J//dfO1y9Jknq/Ph2QP/c5uPBCOPLIdZfffDNceWWeRnrTTeG883KYgxzuTjklh7EhQ+CEE2DGjHX3P/dcGDAgt1a+9FL3XMuG4t13872bMwdefz3f34YW+M4ec/lyePnl3Ne4eNypU2GPPfJDkiSpPQOqXUA17b8/7LJLDrmPP964fMECOP30/GUwgLVrYdtt8/M338ytkY89loPZqlXN+8QOG5Z/brQRrFxZ+evoiXbbLXdzgNwCP3Jkfr7JJjm4Pv10fn3hhTncdkbDMUeMyC33kMPy4MHwzjvwL/8CTz3VuXNIkqS+o08H5IjcdxXWDcgjRuRREsaPb77PZZfl/WbNgs02y10vOjISQ18xc2bLy8eMgdmzobY06/msWXD00Z0710475Q8i9fUwfHj+0PLKK7Dzzrnf8R/+kPuSF73wQvOWf0mSJOjjXSxac9ZZuV/xwoW5/+pLLzUG6eXLc9eKwYNz+Joypbq1bmhOPDF/qFi6FB59NLfsHnfcuttccknLQ+6tXp2DcEqNz9euze/HYYfBtdfmZZMn5xbrPfaAvfbK2zc8rrwSPv95w7EkSWqdAbkFF1+cv+x1wAE5fB1/fO5aATlg1dXl1uMTT+x862dfM2lSHo5txIj8Rb2pU/PzojfeyCOCNHXOOXmUiiefzMO6bbwx/PrXed2UKXn4uKFD4fbb4Y47cku/JElSR0VKqdo1/EVtbW2qq6urdhmS1LtF5D+plLm8tc0laUMXEdNTSrVNl9uCLEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUkGfDcgzZ8Ihh8Dmm8PWW8NnP5uXX3NNHuN48GDo3x8GDmx8/dpreXzegQPz7G3jx8MDD1T1MjY4q1fD2WfnqaBHjuzYJCtTpkBNTd732GPzlNIN6uvzezNoEIwblycCafDWW3ms5aFDYfvt4aqruuhiJElSr9RnA/Kxx8KRR8KSJXkikE99Ki+//PI8bfGKFXks5BtvbHzdMF3yjTfC4sVw/vk5WDsyXfkmT87jFS9YALfckidlWbCg/f1mzICvfhV+/vM8TvLq1fD//l/j+okT8/jKS5bASSflR4OvfS1vv3Ah/OY38N3v5uNIkiS1pE8G5Lfegrlzc0tm//65FXnixI4dY9CgPCPbCSfATTdVps7eaNo0uOCCfM8nTMit8Pfe2/5+v/41HHQQ7LZbvvd///dwzz153bJleabDSy/Nk4dMmgTz58Nzz+X18+fDUUfldSNH5nPOmlWxS5QkSRu4PhmQhw3LQWnixDzd8fvvr/+x9trLaYs7Ys4cGDs2z6J3zz2w6655WXuaTlKQUu5i8ac/5Q87DV1hDjwwd7fYccfG4375yzlAv/tunqFv+nT45Ce7/NIkSVIv0ScDcr9+8H//BwMG5K4W22wDV1+9fscaPDi3YKo8776b79mcOfD667k/8YoV7e930EHwyCO5e8aKFfCf/5mXv/de4zGXL4eXX87BuXjcvfeGd97J04OPHg1f+lL+YCNJktSSPhmQAXbaCX7849xn9e674dvfhl/8ouPHWbEChgzp+vp6q002yffs6adzN4nly3O4bc9ee8E3vpE/0Oy8M+y3X16+6aaNxxwxAhYtgn32Wfe4J58Me+6Zg/T8+XDrrR37cqAkSepb+mxAbtCvXx7NYsKE9euXOmOGrZEdMWYMzJ7d+HrWrBx4y/F3fwevvJJbnj/2sTyixWab5Q87K1fmrhUAq1bl7RqO+8wzcOaZ8OEP5641Rx6Z/4IgSZLUkj4ZkNesgSuuaAxUzz4Ljz2Whwcr13vv5RboadNyH1eV58QT4YYbYOnS3P/7qafyEGxFl1ySP7A09fzzsHZt7nN82WU5MENuwT/sMLj22hyUJ0/OQXiPPfL6vfeGH/0oj2SxaFEemm/XXSt6mZIkaQPWJwNyv345ZO23X/4z/PHH56HADj64vP3PPx+22gr+/d9z94za2srW25tMmpSHYxsxIn9Rb+rU/LzojTfyl+ma+sd/zGF4//3zqBTFYd6mTMkBeuhQuP12uOMOiMjrfvCD3NK/1Va5tf+gg+C88yp1heqxamryP4qmj5oaS5IkrSNS0+EBqqi2tjbVOaiwpEqIaD4cSlvLu0HVSurgiat4iySpoiJiekqpWVNnn2xBliRJklpjQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFBmRJkiSpwIAsSZIkFRiQJUmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJCrYbfd4JFHql2FJEmSWtB3A3JNDWy8MQwe3Ph46KHy94+AuXPX79wzZ8KECeu3b2+2di2cdBIMH57v77x5665/5BHYeWfYZBM49lh4553GdfX1+Z4OGgTjxsELL3Rj4ZIkqTepaECOiM0j4q6ImB0RL0bE+Eqer8N+9jNYsaLxceih1a5I++8Pd93VfPl778EJJ8CVV8Kbb+YAfdlljesnToTdd4clS3LIPumk7qtZkiT1KpVuQf428MuU0lhgT+DFCp+va5x5Jpx9NhxwAGy6KZx1VuO6T386tzYD7Llnfn7RRY3r583L4W3aNBg1Ku9/3XV53SWX5O379WveWr1mDVx9dW7Z3mYbuOACWL26cf0118B228GQIVBbC6+/XoELr7J+/eDCC/P1NfWrX8Fmm8Gpp+YW5K98Be64I69btgwefBAuvTT/VWDSJJg/H557rnvrlyRJvULFAnJEDAH+BpgKkFJalVJaWqnzdbkHH4Rbb83dIaZNg9/9Li+///7c2gzw7LP5+be+1Xz/O++E6dPhjTfgsMPysuuuy9uPHNl8++uvh5/8BJ54InfdmDULbrghr5szB77+dXjySVi6FG66CT784a6/5p5szhwYOxYefxwOPzx3tfjTn2Dx4ny/Bg7MHz4OPDB3t9hxx7yPJElSB1WyBfkjwGLgvyLi9xFxc0Rs0nSjiJgYEXURUbd48eIKltOCY4+FzTdvfLz5ZuO6o4/OQXbkSNhjD3jppY4d+4orYMstc5/Yvfduf/ubb87dB3bYIbc6n3ce3HNPXheRW5hnz84/9903H7sveffdHIAXL87Bd+DAvHzFisZ1y5fDyy/D22/ne9jwQUZSp7xKTf5/qKVHTU21y5OkLlfJgDwAGAd8J6W0N/AucGnTjVJK300p1aaUarfaaqsKltOC++7LLbINj623blw3bFjj8402gpUrO3bsj360Y9svWACnn94Y1r/whRwGAcaMgSlTcivysGG5L+6yZR07/oZuk01y4D3uOHj11cbuJ4MHN64bMQIWLYJ99slhuaErjKROqWE+pNTyY/78apcnSV2ukgG5HqhPKf229PoucmDuHSLaXj9gQMeON2IEPPBAY1hftmzdVuszz8zdL155JS///vc7XPIGbcyY3ILeYNas/GFhq61gp53yB5j6+rxu1ap8n3beuTq1SpKkDVrFAnJKaRGwICIaUsohwKxKna/bbbtt1w4ldtZZuVvGwoW5Veall3I/aMhh7+GHc/AbMCAPhzZkSNeduyd5//38aHje0HJ/8MF5WLdbb81dKr75TTjxxLxuyJDcz/vaa/P2kyc3do2RJEnqoEqPYvH3wI8j4jlgL+CaCp+vY44+et1xkH/4w/L3veYaOP/83Ge4ONxYa9asaTzPa681nnvq1Lz+4ovzF8wOOCAHvuOPb+wTvWpVHqFhiy1yS+r48bk7Rm+0886NXSPGjs2jUkDuyz1tGlx1VW41Xrs2B+IGU6bA88/D0KFw++15hIv2WvklSZJaECmlatfwF7W1tamurq7aZUjqjSLyX2fKXd4NqlZSR0/cVkFVvH+S1FkRMT2l1Gx82b47k54kSZLUAgOyJEmSVGBAliRJkgrKHossIg4A9gVeSCk9WLmSJEmSpOpptQU5In5XeH4OcCOwKXBlRDSb8EOSJEnqDdrqYvGhwvOJwCdTSlcDhwGfr2hVlXbVVfmb1/ffn18vXZpny5swoZpV9R719fleDhoE48Z1bLzoG26AbbbJk4AUh897+WU4/PC8fOut4Ywz8rjIkGfVO+MM2G472Gwz+MQnYObMLr0kSZLUd7QVkPtFxNCI2II8HNxigJTSu8AH3VJdJX30o3DXXfn5T38KNTVVLadXmTgRdt8dliyBk07Kj3L87nf5w8vDD+cxjW+/He68M69btgxOOSVPMz1vXp4Q5OKL87o1a/Jsek8/DW+/DZ/5DBx7bCWuTJIk9QFtBeTNgOlAHTAsIrYFiIjBwIY/A8O++8Lvfw8ffAB33w2f/WzjujVr4Oqrc2jeZhu44ILcSgk5oB18cG7J3G67PAX0ihV53bx5uWV68mTYcksYPRp++1v6lGXL8gyAl16aJ/mYNAnmz4fnnmt/32nT8gQpu+2WJ2D50pdySAbYZx/4whdyC/GgQXDqqfDUU3ndwIF5FsLhw6FfP/jiF2HuXFi8uHLXKUmSeq1WA3JKqSal9JGU0ujSz0WlVWuB47qnvAo76KDcevzGG7Djjo3Lr78efvITeOKJHLRmzcp/+oc8/fE55+TZ8ObMyT+vuGLd4y5bBosWwTHH5KDdl8ydmwPr4MF5ZsD6+nxv58xpf985c/Lsed/+du5eseuure/31FPwsY+1vm777fPMg5IkSR3U5jBvke0XEcdHxHERsR/w55TSq91UX2V97nNw4YVw5JHrLr/5ZrjyytyKuemmcN55cM89ed3YsflP/YMH5ymhTzgBZsxYd/9zz4UBA+CII+Cll7rnWnqKd9/N92b58txv+O238z1saGUvZ98FC/K+re03fTr813/BP/1T83XvvJPf0+uvz63JkiRJHdTqMG8RcRhwE/Ay8Hpp8XBgp4j4cq8Y6m3//WGXXXLIffzxxuULFsDppzcGrLVrYdtt8/M338xdLh57LAe6Vavyn/+Lhg3LPzfaKPeV7Us22SSH2hEjcis65LA8eHD5+37zm/n1vfc232/evNwd5pZbcr/jovffh+OO61i/Z0mSpCbaGgf528ChKaV5xYURMRr4BbBLBevqHhG5vyysG5BHjIAf/ADGj2++z2WX5f1mzcr9YW+4ofGLZMqhdeXK3LVi+PD8AeKVV2Dnndvfd8wYmD278fWsWevu9+ab8KlPwde/nke0KFqzJvdL3mmnlluWJUmSytTW36AHAPUtLH+ddSz2jxQAABaJSURBVIeA633OOiv3K164EFLK3SQagvTy5blrxeDB+Qt7U6ZUt9aeZsgQOOwwuPbaHJQnT4aRI2GPPdbdbsIEuOSSdZedcELuyjJzJrz+Okyd2tgS/M47ORSfey6cdlrz8/7t3+af3/lOl1+SJEnqW9pqQf4+8HRE3A4sKC0bAZwMTK10YVV18cX5z/UHHJBHQhg1Ko/KALlv8hln5NbjXXaBo49et/VZ+UPDaafB0KG5z/Ydd+RW96J585oPrbfffvn+fuITedSQv/3bxoB833151JGXXoKvfa1xnxUr8igZU6fmUTM226xx3f335y8KSpIkdUCklFpfGbELcAywA3lot3rgpymlWZUopra2NtXV1VXi0JL6uoj8F6Fyl3eDqpXU0RO3VVAV758kdVZETE8p1TZd3lYLMimlF4EXK1aVJEmS1MOs1zhYEXFVF9chSZIk9QjrO1Ds9C6tQpIkSeoh1isgp5R+1tWFSJIkST1BWxOFDADOJk8rvT2QgIXAT4CpKaXV3VKhJEmS1I3aakH+EbAXcBVwBHAkcDWwJ/DfFa+s0mbOhEMOgc03h623zrOzAVxzTR7jePBg6N8fBg5sfP3aa3n83oED86xv48fDAw9U9TJ6pPr6fJ8GDYJx4+CFF8rbb9WqPK331lvDFlvAV7+67rfjI/J9b3g/vvvd5sd47LG83c03d8mlSJKkvqetgDwupXRuSuk3KaX60uM3KaVzgb27q8CKOfZYOPJIWLIkj637qU/l5ZdfnsfWXbEij6F7442Nr0eOzNvceGMeH/n883Owdmi6dU2cCLvvnu9tR6Z9/s534Le/hTlz8pjHd9wBP/7xuts8+2zj+zFx4rrrPvggj1c9dmzXXIckSeqT2grIb0fECRHxl20iol9EnAS8XfnSKuitt2DuXDj77NxKvPnmzcNWewYNgs9/Ps/+dtNNlalzQ7RsWZ518NJL88QdkybliTyee679fR99NE8XPXRo/jByyil5Zr1y/fu/5w8922yz/vVLkqQ+r62AfDLwOeCNiHgpIl4CFgHHl9ZtuIYNywFs4sQcyt5/f/2PtddeMGNG19W2oZs7t7FbyoEH5u4WO+6YW4Xb03SygZTy8Yr+5m9gu+3gi1/M0083WLQI/uu/4B/+ofPXIEmS+rRWA3JKaV5K6aSU0lbAeGD/lNLWpWWvdl+JFdCvH/zf/8GAAbmrxTbbwNVXr9+xBg/OrabK3n0335Ply+Hll+Htt2HTTXOXiPZMmAC33pq7ZsybB3fdBe+917j+ySdzP/BnnoGFC+GCCxrXfeUruXvMwIFdfUWSJKmPaXMmvQYppSWVLqTb7bRT7t+6di386le5q8Rf/RUccUTHjrNiBQwZUpkaN0SbbJLvyYgRuVUXclgePLj9fc89N4fqPfbI3SyOPTa/Nw3Gj88/t9sO/umf4PDD8+snnoA//AFO3rD/sCFJknqG9Z0opPfo1y+PZjFhAsya1fH9Z8zI3SyU7bQTrFyZu1ZAHpnilVdg553b33ejjfIXIBcuzKOMRMCee7a8bb/CP92nn4annsrbR+RuM+ecAxdd1PnrkSRJfU6rAbk0DnLvtGYNXHFFY4h79tk8PNi4ceUf4733cgv0tGnw5S9Xps4N0ZAhcNhhcO21OShPnpz7e++xx7rbTZgAl1yy7rJ33819jteuhYceykO1NXx58vnn88gWa9bkLhhXXQVHH53XXXRR7q/c8DjoIPje9+Bb36r01UqSpF6orRD8m4ioB34J/DKlNK97SuoG/frlILbffvmLXttsA1/7Ghx8cHn7n38+XHhhDn133w21tZWtd0MzZQqcdlruJjF2bB6uLWLdbebNg5qadZe9914ehWLBAth++zw6SEO3isWL4UtfgjfeyN01jjoK/u3fuuNqJElSHxOp6cgBxZURo4BPA4cDOwCPA/cDj6aUOjH0Q8tqa2tTnWMKS6qEiOYjpbS1vBtUraSOnritgqp4/ySpsyJiekqpWUtnm32QU0rzU0r/mVI6Ftgf+BlwKPBYRPy8MqVKkiRJ1VN2P+OU0mrg4dKDiNihUkVJkiRJ1bLeo1iklF7vykIkSZKknsBh3iSpA2pqGkcU7Oxj1KjKnrx+QE3X34DOau0amn5pV5KqqN0uFhGxe0rphe4oRpJ6uvnzq/idtFZO3tr35IY3HT2mJ2jtBvbEWiX1WeW0IP9nRPwuIr4cEZtXvCJJkiSpitoNyCmlA4DPAyOAuoi4NSI+WfHKJEmSpCooqw9ySull4GvAJcBBwA0RMTsijq9kcZIkSVJ3azcgR8THImIy8CJwMHB0SmmX0vPJFa5PkiRJ6lbljIN8I/A94PKU0p8bFqaUFkbE1ypWmSRJklQF5QTkI4A/p5TWAEREP2BgSum9lNKPKlqdJEmS1M3K6YP8ELBx4fWg0jJJkiSp1yknIA9MKa1oeFF6PqhyJUmSJEnVU05AfjcixjW8iIh9gD+3sb0kSZK0wSqnD/JFwLSIWFh6vR1wUuVKkiRJkqqn3YCcUno6IsYCOwMBzE4pra54ZZIkSVIVlNOCDPBXQE1p+70jgpTSLRWrSpIkSaqSdgNyRPwI2BGYAawpLU6AAVmSJEm9TjktyLXArimlVOliJEmSpGorZxSLF4BtK12IJEmS1BOU04K8JTArIn4HvN+wMKX0mYpVJUmSJFVJOQH5qkoXIUmSJPUU5Qzz9mhEjAI+mlJ6KCIGAf0rX5okSZLU/drtgxwR5wB3AVNKi3YA7qtkUZIkSVK1lPMlvfOAjwPLAFJKLwNbV7IoSZIkqVrKCcjvp5RWNbyIiAHkcZAlSZKkXqecgPxoRFwObBwRnwSmAT+rbFmSJElSdZQTkC8FFgPPA38L/AL4WiWLkiRJkqqlnFEs1gLfKz0kSZKkXq3dgBwRr9JCn+OU0kcqUpEkSZJUReVMFFJbeD4QOAEYVplyJEmSpOpqtw9ySmlJ4fF6SulbwMHdUJskSZLU7crpYjGu8LIfuUV504pVJEmSJFVROV0s/q3w/ANgHnBiRaqRJEmSqqycUSw+0R2FSJIkST1BOV0s/qGt9Sml67uuHEmSJKm6yh3F4q+An5ZeHw38GlhQqaIkSZKkaiknIG8JjEspLQeIiKuAaSmlL1WyMEmSJKkayplqeiSwqvB6FVBT7gkion9E/D4i/qeDtUmSJEndrpwW5B8Bv4uIe8kz6h0H3NKBc1wIvAgM6Xh5kiRJUvcqZ6KQfwa+CLwNLAW+mFK6ppyDR8Rw4Ejg5s4UKUmSJHWXcrpYAAwClqWUvg3UR8ToMvf7FvBVYG1rG0TExIioi4i6xYsXl3lYST1dTQ1ENH/U1FTpxKNGVeMwVTWPUS1eRH3/US1eW33/lrefRzdc9KhWah1QU51/R22o2r9tSd0mUkptbxBxJXkki51TSmMiYnvyl/Q+3s5+RwFHpJS+HBETgK+klI5qa5/a2tpUV1fXoQuQ1DNFQEv/vbS2vOIn7uD2Fa+zDa2eu4O1dtV70Ob2XXqw8rfvie9PNWuStH4iYnpKqbbp8nJakI8DPgO8C5BSWkh5U01/HPhMRMwDbgcOjoj/LrtiSZIkqQrKCcirUm5mTgARsUk5B04pXZZSGp5SqgFOBh5OKZ223pVKkiRJ3aCcgHxnREwBNo+Ic4CHgO9VtixJkiSpOtod5i2l9M2I+CSwDNgZuCKl9L8dOUlK6RHgkfUpUJIkSepObQbkiOgPPJBSOhToUCiWJEmSNkRtdrFIKa0B3ouIzbqpHkmSJKmqyplJbyXwfET8L6WRLABSShdUrCpJkiSpSsoJyD8vPSRJkqRer9WAHBEjU0qvpZR+2J0FSZIkSdXUVh/k+xqeRMTd3VCLJEmSVHVtBeQoPP9IpQuRJEmSeoK2AnJq5bkkSZLUa7X1Jb09I2IZuSV549JzSq9TSmlIxauTJEmSulmrATml1L87C5EkSZJ6gjYnCpEkSZL6GgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFBmRJkiSpwIAsSZIkFRiQJUmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZKkXq6mBiOaPmppqV9aC1opt5VE/oKZrzjtqVIvHf5UuOn5bWrnmBf1rWrxsaPl2jBrVoUtrdftqau3tb01r19Yj/223ZoP6BZX6lkgpVbuGv6itrU11dXXVLkPqNSKgpV/x1pZX9dytrOjo9l12cT3xJnVHTS1o87Qdramr3ucuKmd9dNnb08PeZ6kviojpKaXapsttQZYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFBmRJkiSpwIAsSZIkFRiQJUmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVJBxQJyRIyIiF9FxIsRMTMiLqzUuSRJkqSuMqCCx/4AuDil9ExEbApMj4j/TSnNquA5JUmSpE6pWAtySumPKaVnSs+XAy8CO1TqfJIkSVJX6JY+yBFRA+wN/LaFdRMjoi4i6hYvXtwd5UhqRU0NRDR/1NR0cIcIEi0vf5XWDtYx9f1HtXj8+v6jOlpqi49WdfRAbT1GtVwro1q+tla3r6LW3ofWHvMYVdFLe5WaDv4jbuW6BrRynAgW9O/YsSRteCKlVNkTRAwGHgX+OaV0T1vb1tbWprq6uorWI/UlEdDSr3hXLW99xXrs08ryNk7RofN29Dhdd6ANX1uXXPHbUa1/AOtx0f4bkzY8ETE9pVTbdHlFW5Aj4kPA3cCP2wvHkiRJUk9QyVEsApgKvJhSur5S55EkSZK6UiVbkD8OnA4cHBEzSo8jKng+SZIkqdMqNsxbSulxoK2vukiSJEk9jjPpSZIkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFBmRJkiSpwIAsSZIkFRiQJUmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLHVRTAxHNHzU11a6suVepabHYV6npkuMwalSHa6rvP6rFY81jVJecYlTLh+94qV12oA1fa7eiW25Hayfvol+4+gE1LR6/vn8bF9ZKTR39vYKWL63Duuo/pdaO05XH6on/UUotiJRStWv4i9ra2lRXV1ftMqQ2RUBLvzatLa+qDhbb6jWsx8X1yPuh3qOjv4jd8Yvb0WN1VU2VPk41a5IqLCKmp5Rqmy63BVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFBmRJkiSpwIAsSZIkFRiQJUmSpAIDsiRJklRgQJYkSZIKDMiSJElSgQFZkiRJKjAgS5IkSQUGZEmSJKnAgCxJkiQVGJAlSZKkAgOyJEmSVGBAliRJkgoMyJIkSVKBAVmSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIkiRJUoEBWZIkSSowIEuSJEkFFQ3IEXF4RMyJiLkRcWklzyVJkiR1hYoF5IjoD/wH8GlgV+CUiNi1UueTJEmSukIlW5D3BeamlP6QUloF3A4cU8HzSZIkSZ02oILH3gFYUHhdD+zXdKOImAhMLL1cERFzKliTGm0JvFXtIjZUER1bXkVbEtHy+9xKsa1ew3pcXA+8H71V3/x97ugvYnf84nb0WB2rqfX3uauura3tK3ttatQ3f5+rZ1RLCysZkFv6DUjNFqT0XeC7FaxDLYiIupRSbbXrUGX5PvcNvs99g+9z3+D73DNUsotFPTCi8Ho4sLCC55MkSZI6rZIB+WngoxExOiI2Ak4GflrB80mSJEmdVrEuFimlDyLifOABoD/w/ZTSzEqdTx1mt5a+wfe5b/B97ht8n/sG3+ceIFJq1i1YkiRJ6rOcSU+SJEkqMCBLkiRJBQZkERFfiYgUEVtWuxZ1vYj414iYHRHPRcS9EbF5tWtS14mIwyNiTkTMjYhLq12Pul5EjIiIX0XEixExMyIurHZNqoyI6B8Rv4+I/6l2LX2dAbmPi4gRwCeB16pdiyrmf4HdU0ofA14CLqtyPeoiEdEf+A/g08CuwCkRsWt1q1IFfABcnFLaBfhr4Dzf517rQuDFahchA7JgMvBVWpjERb1DSunBlNIHpZe/IY9Jrt5hX2BuSukPKaVVwO3AMVWuSV0spfTHlNIzpefLyQFqh+pWpa4WEcOBI4Gbq12LDMh9WkR8Bng9pfRstWtRtzkLuL/aRajL7AAsKLyux+DUq0VEDbA38NvqVqIK+Ba5wWpttQtRZaeaVg8QEQ8B27aw6h+By4HDurciVUJb73NK6Selbf6R/KfaH3dnbaqoaGGZfw3qpSJiMHA3cFFKaVm161HXiYijgDdTStMjYkK165EBuddLKR3a0vKI2AMYDTwbEZD/7P5MROybUlrUjSWqC7T2PjeIiC8ARwGHJAc/703qgRGF18OBhVWqRRUUER8ih+Mfp5TuqXY96nIfBz4TEUcAA4EhEfHfKaXTqlxXn+VEIQIgIuYBtSmlt6pdi7pWRBwOXA8clFJaXO161HUiYgD5i5eHAK8DTwOnOmtp7xK5FeOHwJ9SShdVux5VVqkF+SsppaOqXUtfZh9kqfe7EdgU+N+ImBER/1ntgtQ1Sl++PB94gPzFrTsNx73Sx4HTgYNLv8MzSi2NkirEFmRJkiSpwBZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJUYECWJEmSCgzIktRDRcRxEZEiYmw7202JiI83WXZVRLxeGhJsdkR8JyL8P1+SyuB/lpLUc50CPA6c3M52+wG/aWH55JTSXsCuwB7AQV1bniT1TgZkSeqBImIweYKIs2kjIEfELsBLKaU1bRxuI/L0tW93aZGS1EsZkCWpZzoW+GVK6SXgTxExrpXtPg38spV1kyJiBvBHcoieUYE6JanXMSBLUs90CnB76fntpdct+RStB+SGLhZbA5tERHtdNSRJwIBqFyBJWldEbAEcDOweEQnoD6SI+GpKKRW2GwRsnlJa2NbxUkqrI+KXwN/QGLolSa2wBVmSep7PAbeklEallGpSSiOAV4EDmmz3CeBX7R0sIgLYH3ilyyuVpF7IgCxJPc8pwL1Nlt0NnNpkWVv9j6GxD/IL5L8Y3tRlFUpSLxaFv9ZJkjYgEfEMsF9KaXW1a5Gk3sSALEmSJBXYxUKSJEkqMCBLkiRJBQZkSZIkqcCALEmSJBUYkCVJkqQCA7IkSZJU8P8DBfyOmlS3cAQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "for iexp in range(N_exp):\n",
    "    if ((iexp+1)%1000 == 0):\n",
    "        print(f\"Got to experiment number: {iexp+1}\")\n",
    "\n",
    "    # Generate data:\n",
    "    x_A_array = r.normal(dist_mean_A, dist_width_A, N_events_A)\n",
    "    x_B_array = r.normal(dist_mean_B, dist_width_B, N_events_B)\n",
    "    \n",
    "    \n",
    "    # Test if there is a difference in the mean:\n",
    "    # ------------------------------------------\n",
    "    # Calculate mean and error on mean:\n",
    "    mean_A, width_A, sdom_A = mean_std_sdom(x_A_array) \n",
    "    mean_B, width_B, sdom_B = mean_std_sdom(x_B_array) \n",
    "\n",
    "    # Consider the difference between means in terms of the uncertainty:\n",
    "    d_mean = mean_A - mean_B\n",
    "    # ... how many sigmas is that away?\n",
    "\n",
    "    # Turn a number of sigmas into a probability (i.e. p-value):\n",
    "    p_mean    = 0.5      # Calculate yourself. HINT: \"stats.norm.cdf or stats.norm.sf may be useful!\"\n",
    "    all_p_mean[iexp] = p_mean\n",
    "    \n",
    "    \n",
    "    # Test if there is a difference with the chi2:\n",
    "    # --------------------------------------------\n",
    "    # Chi2 Test:\n",
    "    p_chi2 = 0.5         # Calculate the p-value of the Chi2 between histograms of A and B yourself.\n",
    "    all_p_chi2[iexp] = p_chi2\n",
    "\n",
    "    \n",
    "    # Test if there is a difference with the Kolmogorov-Smirnov test on arrays (i.e. unbinned):\n",
    "    # -----------------------------------------------------------------------------------------\n",
    "    p_ks = stats.ks_2samp(x_A_array, x_B_array)[1]           # Fortunately, the K-S test is implemented in stats!\n",
    "    all_p_ks[iexp] = p_ks\n",
    "\n",
    "\n",
    "    # Print the results for the first 10 experiments\n",
    "    if (verbose and iexp < 10) :\n",
    "      print(f\"{iexp:4d}:  p_mean: {p_mean:7.5f}   p_chi2: {p_chi2:7.5f}   p_ks: {p_ks:7.5f}\")\n",
    "\n",
    "    \n",
    "    # In case one wants to plot the distribution for visual inspection:\n",
    "    if (iexp == 0):\n",
    "        \n",
    "        ax1.hist(x_A_array, N_bins, (xmin, xmax), histtype='step', label='A', color='blue')\n",
    "        ax1.set(title='Histograms of A and B', xlabel='A / B', ylabel='Frequency / 0.05')        \n",
    "        ax_text(x_A_array, ax1, 0.04, 0.85, 'blue')\n",
    "\n",
    "        ax1.hist(x_B_array, N_bins, (xmin, xmax), histtype='step', label='B', color='red')\n",
    "        ax_text(x_B_array, ax1, 0.04, 0.65, 'red')\n",
    "        \n",
    "        ax1.legend()\n",
    "        fig1.tight_layout()\n",
    "\n",
    "        \n",
    "fig1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "## Show the distribution of hypothesis test p-values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAPoCAYAAADDVV/dAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5xWdb33/9cHyTgMiAiRIYgh4HFnSh5Kd6RmZplooeb5UPx+W807t+7Ebrfgru3GtlvT7a609M4sA7F2auYtlZqHnRqWqaAIJsJ4AFJOgyIKn/uPa800jMPMAGtmuIbX8/G4HnOttb5rfT9rrvUA3ny/a12RmUiSJEmSNl23zi5AkiRJkroKA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJXVxEzIyIMZ1dRzWLiPsj4ksbue8PI+KbLWyvi4gPNm0bEQdFxOyNq1iS1FkMWJJUxSJiXkQc2mTdaRHxUP1yZu6emfe3cpxhEZER0b2dStV6ZGZNZv6lmfUPZuao+uXmPmtJ0ubHgCVJanebc3DbnGuTJFUfA5YkdXGNRz4iYt+ImBERyyNiYURcWTR7oPi5tJiydkBEdIuIiyPixYhYFBE/iohtGh33lGLbaxHxz036mRQRt0XEjyNiOXBa0ffvI2JpRLwSEddGxNaNjpcRcVZEzImIFRHxjYgYXuyzPCJubdy+lXP+YUR8LyJ+XRzrdxGxY5O+zo6IOcCcYt1HI+IPEbGs+PnRJocdHhGPFdtvj4j+jY43LSJeLbY9EBG7N9l3QCu17NzMOYyJiNri/c3AUODO4vP5WkTcFRFfabLPkxExtplj1Y9Qnh4RCyJiSUT8/xHxkWKfpRFxbaP2kyLix83sbxiVpFYYsCRpy3I1cHVm9gWGA7cW6/+++NmvmLL2e+C04vUJ4INADXAtQETsBnwHOBHYHtgGGNykr6OA24B+wE+ANcB5wADgAOAQ4Kwm+xwO7APsD3wNuL7oYwiwB/DFDTjXE4FvFP09UdTQ2FhgP2C3IizdBVwDbAdcCdwVEds1an8KcAbwAeCdom29u4ERwPuAPzbTV2u1tCgzTwbmA0cWn8+3gJuAk+rbRMSHqHwGv2rhUPsVdR4HfBv438ChwO7AsRHx8Q2pS5L0bgYsSap+vyhGIJZGxFIqwWd93gZ2jogBmVmXmY+00PZE4MrM/Etm1gEXAccXoxhfAO7MzIcyczVwCZBN9v99Zv4iM9dm5puZ+XhmPpKZ72TmPOA6oOk/6C/PzOWZORN4Gphe9L+MSoj5cNt+JQDclZkPZOZbVILEARExpNH2f8vM1zPzTeAzwJzMvLmo76fAs8CRjdrfnJlPZ+ZK4J+pBJKtADLzxsxcUfQ1CfhQ49G+NtSyMW4HRkTEiGL5ZGBq8Xmszzcyc1VmTgdWAj/NzEWZ+RLwIBv2+5UkNcOAJUnVb2xm9qt/8e5RocbOBEYCzxbT4D7bQtsPAC82Wn4R6A4MKrYtqN+QmW8ArzXZf0HjhYgYGRG/LKbSLQcuozKi09jCRu/fbGa5poV6m2pcXx3welF3c/U1PVeK5cajcguabHsPlal/W0XE5Ih4vjiveUWbAc3tu55aNlgR1m4FToqIblRG925uZbcyf7+SpGYYsCRpC5KZczLzi1Smsl0O3BYRvXn36BPAy8COjZaHUpkatxB4BdihfkNE9KQytW6d7posf5fKqNCIYori14HY+LNpVcMIUUTUAP2pnFNz9TU9V6ic70vNHa/Y9jbwV+AEKtMhD6UyVXJYfbcbUEtbNPcZ3URlpPEQ4I1iamcZVgK9Gi2/v6TjSlKXZ8CSpC1IRJwUEQMzcy2wtFi9BlgMrKVyr1W9nwLnRcRORSi4jMoUtHeo3Ft1ZPFgiK2BS2k9LPUBlgN1EbEL8A+beC4ZLX+/1xERcWBR3zeARzNzwXra/goYGREnRET3iDgO2A34ZaM2J0XEbhHRC/gX4LbMXFOc11tURvB6Ufk9bUot67OQdT8fikC1FvgPWh+92hBPAH8fEUOLqY4XlXhsSerSDFiStGU5HJgZEXVUHnhxfHFPzhvAvwIPF/dy7Q/cSOUf7Q8ALwCrgK8AFPdIfQWYQmU0awWwiErQWJ8LqIz2rAC+D0zd2JOIiB2AOuCpFprdAkykMh1vHyojPc3KzNeAzwLnUwlKXwM+m5l/bdTsZuCHwKtAD+DcYv2PqEwZfAmYBTR3X1uba2nBvwEXF5/PBY3W/wjYE/hx87ttuMz8NZXP50ngcdYNmpKkFkRmczMOJElqu2KEaymV6X8vdEB/JwG7Z2azIysR8UOgNjMvbu9aOltEnAKMz8wDO7sWSVLlZmVJkjZYRBwJ/JbK1MArqIwmzeuIvjOztNGaalZMVzyLlp8cKUnqQE4RlCRtrKOoPKjhZSrfrXR8Oi2iw0TEp6jcO7eQyhRESdJmwCmCkiRJklQSR7AkSZIkqSRd8h6sAQMG5LBhwzq7DEmSJElV6vHHH/9rZg7c0P26ZMAaNmwYM2bM6OwyJEmSJFWpiHhxY/ZziqAkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVJJ2C1gRcWNELIqIpxut6x8Rv46IOcXPbYv1ERHXRMTciHgyIvZutM+pRfs5EXFqe9UrSZIkSZuqPUewfggc3mTdBOC3mTkC+G2xDPBpYETxGg98FyqBDJgI7AfsC0ysD2WSJEmStLlpt4CVmQ8ArzdZfRRwU/H+JmBso/U/yopHgH4RsT3wKeDXmfl6Zi4Bfs27Q5skSZIkbRa6d3B/gzLzFYDMfCUi3lesHwwsaNSutli3vvXvEhHjqYx+MXTo0JLLliR1JR+bfC8vLX1zo/cf3K8nD084uMSKJEldRUcHrPWJZtZlC+vfvTLzeuB6gNGjRzfbRpIkgJeWvsm8yZ/Z6P2HTbirxGokSV1JRz9FcGEx9Y/i56JifS0wpFG7HYCXW1gvSZIkSZudjg5YdwD1TwI8Fbi90fpTiqcJ7g8sK6YS3gMcFhHbFg+3OKxYJ0mSJEmbnXabIhgRPwXGAAMiopbK0wAnA7dGxJnAfGBc0fxXwBHAXOAN4HSAzHw9Ir4B/KFo9y+Z2fTBGZIkSZK0WWi3gJWZX1zPpkOaaZvA2es5zo3AjSWWJkmSJEntoqOnCEqSJElSl2XAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKkmnBKyIOC8iZkbE0xHx04joERE7RcSjETEnIqZGxNZF2/cWy3OL7cM6o2ZJkiRJak2HB6yIGAycC4zOzD2ArYDjgcuBqzJzBLAEOLPY5UxgSWbuDFxVtJMkSZKkzU5nTRHsDvSMiO5AL+AV4GDgtmL7TcDY4v1RxTLF9kMiIjqwVkmSJElqkw4PWJn5EnAFMJ9KsFoGPA4szcx3ima1wODi/WBgQbHvO0X77ZoeNyLGR8SMiJixePHi9j0JSZIkSWpGZ0wR3JbKqNROwAeA3sCnm2ma9bu0sO1vKzKvz8zRmTl64MCBZZUrSZIkSW3WGVMEDwVeyMzFmfk28HPgo0C/YsogwA7Ay8X7WmAIQLF9G+D1ji1ZkiRJklrXGQFrPrB/RPQq7qU6BJgF3Ad8oWhzKnB78f6OYpli+72Z+a4RLEmSJEnqbJ1xD9ajVB5W8UfgqaKG64ELgX+MiLlU7rG6odjlBmC7Yv0/AhM6umZJkiRJaovurTcpX2ZOBCY2Wf0XYN9m2q4CxnVEXZIkSZK0KTrrMe2SJEmS1OUYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJC0GrIjYJSIOiYiaJusPb9+yJEmSJKn6rDdgRcS5wO3AV4CnI+KoRpsva+/CJEmSJKnadG9h25eBfTKzLiKGAbdFxLDMvBqIjihOkiRJkqpJSwFrq8ysA8jMeRExhkrI2hEDliRJkiS9S0v3YL0aEXvVLxRh67PAAGDP9i5MkiRJkqpNSwHrFODVxisy853MPAX4+3atSpIkSZKq0HqnCGZmbVTsBwwGEngZeCwzH+6oAiVJkiSpWqw3YEXEYcB3gDnAS8XqHYCdI+KszJzeAfVJkiRJUtVo6SEXVwOHZua8xisjYifgV8Cu7ViXJEmSJFWdlu7B6g7UNrP+JeA97VOOJEmSJFWvlkawbgT+EBFTgAXFuiHA8cAN7V2YJEmSJFWblh5y8W8R8QvgKOAAKt99VQucmJmzOqg+SZIkSaoaLY1gkZnPAM90UC2SJEmSVNVaugdrvSJiUsl1SJIkSVLV26iABTxeahWSJEmS1AVsVMDKzDvLLkSSJEmSql1LXzTcHTgTOBr4AJDAy8DtwA2Z+XaHVChJkiRJVaKlEaybgb2AScARwGeAS4EPAT/elE4jol9E3BYRz0bEMxFxQET0j4hfR8Sc4ue2RduIiGsiYm5EPBkRe29K35IkSZLUXlp6iuDemTmqybpa4JGIeG4T+70a+L+Z+YWI2BroBXwd+G1mTo6ICcAE4ELg08CI4rUf8N3ipyRJkiRtVloawVoSEeMioqFNRHSLiOOAJRvbYUT0Bf6e4suKM3N1Zi6l8n1bNxXNbgLGFu+PAn6UFY8A/SJi+43tX5IkSZLaS0sB63jgC8DCiHiuGLV6FTim2LaxPggsBv5PRPwpIn4QEb2BQZn5CkDx831F+8HAgkb71xbr1hER4yNiRkTMWLx48SaUJ0mSJEkbZ71TBDNzHnAcQERsB0Rm/rWkPvcGvpKZj0bE1VSmA65PNFfeu1ZkXg9cDzB69Oh3bZckSZKk9tamx7Rn5mslhSuojEDVZuajxfJtVALXwvqpf8XPRY3aD2m0/w5UnmYoSZIkSZuVjf2i4Y2Wma8CCyKi/gEahwCzgDuAU4t1p1J5HDzF+lOKpwnuDyyrn0ooSZIkSZuTFr8HKzPfaad+vwL8pHiC4F+A06mEvVsj4kxgPjCuaPsrKo+Jnwu8UbSVJEmSpM1OS49pfyQiaoH/S+WR6vPK6jQznwBGN7PpkGbaJnB2WX1LkiRJUntp6SEXoyNiRyrfQ/XtiBgMPATcDfwuM9/qoBolSZIkqSq0eA9WZr6Ymd/LzLHAR4E7gUOBByPiro4oUJIkSZKqRUtTBNeRmW8D9xYvihEtSZIkSVJho58imJkvlVmIJEmSJFW7Dn9MuyRJkiR1Va0GrIjYoyMKkSRJkqRq15YRrO9FxGMRcVZE9Gv3iiRJkiSpSrUasDLzQOBEYAgwIyJuiYhPtntlkiRJklRl2nQPVmbOAS4GLgQ+DlwTEc9GxDHtWZwkSZIkVZO23IP1dxFxFfAMcDBwZGbuWry/qp3rkyRJkqSq0ZbvwboW+D7w9cx8s35lZr4cERe3W2WSJEmSVGXaErCOAN7MzDUAEdEN6JGZb2Tmze1anSRJkiRVkbbcg/UboGej5V7FOkmSJElSI20JWD0ys65+oXjfq/1KkiRJkqTq1JaAtTIi9q5fiIh9gDdbaC9JkiRJW6S23IP1VWBaRLxcLG8PHNd+JUmSJElSdWo1YGXmHyJiF2AUEMCzmfl2u1cmSZIkSVWmLSNYAB8BhhXtPxwRZOaP2q0qSZIkSapCrQasiLgZGA48AawpVidgwJIkSZKkRtoygjUa2C0zs72LkSRJkqRq1panCD4NvL+9C5EkSZKkateWEawBwKyIeAx4q35lZn6u3aqSJEmSpCrUloA1qb2LkCRJkqSuoC2Paf9dROwIjMjM30REL2Cr9i9NkiRJkqpLq/dgRcSXgduA64pVg4FftGdRkiRJklSN2vKQi7OBjwHLATJzDvC+9ixKkiRJkqpRWwLWW5m5un4hIrpT+R4sSZIkSVIjbQlYv4uIrwM9I+KTwDTgzvYtS5IkSZKqT1sC1gRgMfAU8P8BvwIubs+iJEmSJKkateUpgmuB7xcvSZIkSdJ6tBqwIuIFmrnnKjM/2C4VSZIkSVKVassXDY9u9L4HMA7o3z7lSJIkSVL1avUerMx8rdHrpcz8NnBwB9QmSZIkSVWlLVME92602I3KiFafdqtIkiRJkqpUW6YI/kej9+8A84Bj26UaSZIkSapibXmK4Cc6ohBJkiRJqnZtmSL4jy1tz8wryytHkiRJkqpXW58i+BHgjmL5SOABYEF7FSVJkiRJ1agtAWsAsHdmrgCIiEnAtMz8UnsWJkmSJEnVptXHtANDgdWNllcDw9qlGkmSJEmqYm0ZwboZeCwi/htI4GjgR+1alSRJkiRVobY8RfBfI+Ju4KBi1emZ+af2LUuSJEmSqk9bpggC9AKWZ+bVQG1E7NSONUmSJElSVWo1YEXEROBC4KJi1XuAH29qxxGxVUT8KSJ+WSzvFBGPRsSciJgaEVsX699bLM8ttg/b1L4lSZIkqT20ZQTraOBzwEqAzHwZ6FNC3/8LeKbR8uXAVZk5AlgCnFmsPxNYkpk7A1cV7SRJkiRps9OWgLU6M5PKAy6IiN6b2mlE7AB8BvhBsRzAwcBtRZObgLHF+6OKZYrthxTtJUmSJGmz0paAdWtEXAf0i4gvA78Bvr+J/X4b+BqwtljeDliame8Uy7XA4OL9YIovNS62LyvaryMixkfEjIiYsXjx4k0sT5IkSZI2XKsBKzOvoDJy9DNgFHBJZv7nxnYYEZ8FFmXm441XN9d1G7Y1rvP6zBydmaMHDhy4seVJkiRJ0kZr8THtEbEVcE9mHgr8uqQ+PwZ8LiKOAHoAfamMaPWLiO7FKNUOwMtF+1pgCJWnF3YHtgFeL6kWSZIkSSpNiyNYmbkGeCMitimrw8y8KDN3yMxhwPHAvZl5InAf8IWi2anA7cX7O4pliu33FveESZIkSdJmpdUvGgZWAU9FxK8pniQIkJnnllzLhcCUiPgm8CfghmL9DcDNETGXysjV8SX3K0mSJEmlaEvAuqt4lS4z7wfuL97/Bdi3mTargHHt0b8kSZIklWm9ASsihmbm/My8aX1tJEmSJEl/09I9WL+ofxMRP+uAWiRJkiSpqrUUsBo/Hv2D7V2IJEmSJFW7lgJWrue9JEmSJKkZLT3k4kMRsZzKSFbP4j3FcmZm33avTpIkSZKqyHoDVmZu1ZGFSJIkSVK1a/GLhiVJkiRJbWfAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLA6oJ233137r///s4uQ5IkSdriGLDaybBhw+jZsyc1NTUNr9/85jdt3j8imDt37kb1PXPmTMaMGbNR+1a72bNnc/jhh9OvXz+GDRv2ru3XXHMNgwYNon///lx00UXrbLv//vsZNWoUvXv3ZuzYsSxbtqxhW21tLWPGjKFXr17svffePP300+19KpIkSapCBqx2dOedd1JXV9fwOvTQQzu7pC6ve/fuHH/88fz7v//7u7Y99thjTJo0iXvvvZennnqKKVOmcOuttwLwxhtvMG7cOCZOnMiiRYuIiHUC2Pjx49ljjz147bXXOO644zjuuOM67JwkSZJUPQxYneC0007jzDPP5MADD6RPnz6cccYZDds+/elPU1NTA8CHPvQhampq+OpXv9qwfd68eUQE06ZNY8cdd6RPnz5cfvnlAFx44YXU1NTQrVu3d42WrVmzhksvvZRhw4YxaNAgzj33XN5+++2G7Zdddhnbb789ffv2ZfTo0bz00kvt+StoN8OHD+e0005jp512ete2adOmccwxx7D77rszePBgvvSlLzFlyhQA7rvvPrbZZhtOOOEEevfuzQUXXMDUqVMBWL58OdOnT2fChAn07NmT8847jxdffJEnn3yyQ89NkiRJmz8DVieZPn06t9xyCzNnzmTatGk89thjANx9993U1dUB8Oc//5m6ujq+/e1vv2v/W2+9lccff5yFCxdy2GGHAXD55ZdTV1fH0KFD39X+yiuv5Pbbb+fhhx9m7ty5zJo1i2uuuQaoTKv75je/yf/8z/+wdOlSvvOd7/De9763vU6908yePZtddtmFq6++mosuuojddtuN2bNnr7PtoYce4vDDD2fUqFG8/vrrLF68mLlz59KjRw9qamo46KCDqK2tZfjw4Q37SpIkSfUMWO1o7Nix9OvXr+G1aNGihm1HHnkkQ4cOZejQoey5554899xzG3TsSy65hAEDBtCrVy8+/OEPt9r+Bz/4ARMnTmTw4MH06dOHs88+m5///OdA5X6vNWvW8Oyzz7JmzRr23XdfBgwYsGEnWwVWrlxJTU0NCxYsYM6cOfTp06chzNZvW7x4MbNnz6ZHjx4A1NXVNWxbsWIFc+bMYcmSJevsK0mSJNUzYLWjX/ziFyxdurTh9b73va9hW//+/Rveb7311qxatWqDjj1ixIgNar9gwQJOPvnkhrB36qmnsnjxYgBGjhzJddddxze/+U369+/PuHHjWL58+QYdvxr07t2buro6rrjiCm677TZWrFjRMB2zftvRRx/NCy+80DB9sqampmHbkCFDePXVV9lnn33W2VeSJEmqZ8DaTEVEi9u7d+++QccbMmQI99xzT0PYW758+TqjZqeddhoPP/wwzz//PM899xw33njjRtW9ORs5ciTPPvtsw/KsWbMYNWrUerf179+fgQMHsvPOO7Nq1Spqa2sBWL16Nc8//3zDvpIkSVI9A9Zm6v3vf3+pjwI/44wzuOSSS3j55ZfJTJ577jmmT58OwPPPP8+9997L6tWr6d69O2vXrqVv376l9d2RMpNVq1bx9ttvN7xfvXo1AOPGjePnP/85M2fO5KWXXuKGG25oeBrgwQcfzLJly7jllltYuXIlV1xxBcceeywAffv25bDDDmPy5MmsWrWKq666qmFqpyRJktSYAasdHXnkket8D9ZNN93U5n0vu+wyzjnnHAYPHvyu72tqzpo1axr6mT9/fkPfN9xwAwDnn38+Bx10EAceeCB9+/blmGOOabgnbPXq1UyYMIHtttuOkSNHcsABB3DyySdv3El3shdffJGePXtyxBFHMH/+fHr27NnwEJD99tuPiRMn8olPfII99tiDY489tiFg9erVi2nTpjFp0iQGDhzI2rVrmTx5csNxr7vuOp566im23XZbpkyZwtSpU1sdZZQkSdKWJzKzs2so3ejRo3PGjBmdXYYkaTM1bMJdzJv8mU7bX5K0+YuIxzNz9Ibu5wiWJEmSJJWkwwNWRAyJiPsi4pmImBkR/6tY3z8ifh0Rc4qf2xbrIyKuiYi5EfFkROzd0TVLkiRJUlt0xgjWO8D5mbkrsD9wdkTsBkwAfpuZI4DfFssAnwZGFK/xwHc7vmRJkiRJal2HB6zMfCUz/1i8XwE8AwwGjgLqnwJxEzC2eH8U8KOseAToFxHbd3DZkiRJktSqTr0HKyKGAR8GHgUGZeYrUAlhQP238g4GFjTarbZYJ0mSJEmblU4LWBFRA/wM+GpmLm+paTPr3vXow4gYHxEzImLG4sWLyypzo0yaNImI4O677wZg6dKlbL311owZM6ZT6+oqamtrGTNmDL169WLvvfdu8/eF3X///XTr1m2dR+c/88wz62wfNWoUvXv3ZuzYsSxbtmyT+5QkSdKWpVMCVkS8h0q4+klm/rxYvbB+6l/xc1GxvhYY0mj3HYCXmx4zM6/PzNGZOXrgwIHtV3wbjRgxgttuuw2AO+64g2HDhnVuQV3I+PHj2WOPPXjttdc47rjjGr7Lqi0+8IEPUFdX1/DaddddAXjjjTcYN24cEydOZNGiRUTEOt8/til9SpIkacvRGU8RDOAG4JnMvLLRpjuAU4v3pwK3N1p/SvE0wf2BZfVTCTdn++67L3/605945513+NnPfsbnP//5hm1r1qzh0ksvZdiwYQwaNIhzzz2Xt99+G4AXXniBgw8+mP79+7P99ttz2mmnUVdXB8C8efOICK666ioGDBjATjvtxKOPPtop59dZli9fzvTp05kwYQI9e/bkvPPO48UXX+TJJ5/cpOPed999bLPNNpxwwgn07t2bCy64gKlTp7Zrn5IkSep6OmME62PAycDBEfFE8ToCmAx8MiLmAJ8slgF+BfwFmAt8HzirE2reKB//+Me54447WLhwIcOHD29Yf+WVV3L77bfz8MMPM3fuXGbNmsU111wDwFtvvcWXv/xl5s+fz+zZs5k/fz6XXHLJOsddvnw5r776KkcddRSXXnpph55TZ5s7dy49evSgpqaGgw46iNraWoYPH87s2bPbtP+iRYsYNGgQw4cP57LLLqP+i7Znz57NLrvswkMPPcThhx/OqAwHJEMAACAASURBVFGjeP3111m8ePEm9ylJkqQtR/eO7jAzH6L5+6oADmmmfQJnt2tR7eQLX/gCxx9/POPHj19n/Q9+8AO+9a1vMXhw5VkdZ599NldccQXnn38+u+yyC7vssktD23HjxjFt2rR19v+Hf/gHunfvzhFHHMEvf/nL9j+RzcjKlSupqalhxYoVzJkzhyVLltCnT5+GUb6W7LrrrsycOZPhw4fz5JNPctRRR7H99ttz+umnNxx38eLFzJ49mx49egBQV1e3SX1KkiRpy9LhAWtL8tGPfpRdd92VcePG8dBDDzWsX7BgASeffDLdulUGENeuXcv73/9+oDLCcu655/Lggw+ycuVKVq9ezT777LPOcfv37w/A1ltvzapVqzrobDYPvXv3pq6ujiFDhvDqq68CsGLFCmpqalrdd9CgQQwaNAiAvfbai7PPPps777yT008/veG4Rx99NEcffTRLliwBoKamZpP6lCRJ0palUx/T3tVFBNOnT19nRApgyJAh3HPPPSxdupSlS5eyfPlynnvuOQAuuugiIoJZs2axdOlSJk+e3DCNTbDzzjuzatUqamtrAVi9ejXPP/88o0aN2uBj1QdcgJEjR/Lss882LM+aNYv+/fszcODAUvuUJElS12bA6gRnnHEGl1xyCS+//DKZyXPPPcf06dOByshI3759qamp4YUXXuC6667r5Go3L3379uWwww5j8uTJrFq1iquuuoqhQ4ey5557rtNuzJgxXHjhheusu/fee5k/fz4AzzzzDN/97nc58sgjATj44INZtmwZt9xyCytXruSKK67g2GOP3aA+JUmSJANWJzj//PM56KCDOPDAA+nbty/HHHMMixZVnko/ceJEZsyYwTbbbMOxxx7bEAD0N9dddx1PPfUU2267LVOmTGHq1KlUHk75N/PmzWPhwoXrrPvjH//IRz7yEXr16sURRxzBl7/8ZU4//XQAevXqxbRp05g0aRIDBw5k7dq1TJ48eYP6lCRJkqIrTj8bPXp0zpgxo7PLkCRtpoZNuIt5kz/TaftLkjZ/EfF4Zo7e0P0cwZIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsNrJzJkzOeSQQ+jXrx/ve9/7+PznPw/AZZddRk1NDTU1NWy11Vb06NGjYXn+/PmMGTOGHj160Lt3bw444ADuueeeTj6TzU9tbS1jxoyhV69e7L333jz99NNt3veaa65h0KBB9O/fn4suumidbffffz+jRo2id+/ejB07lmXLlpXSpyRJkrYcBqx2MnbsWD7zmc/w2muv8dxzz/GpT30KgK9//evU1dVRV1fHQQcdxLXXXtuwPHToUACuvfZaFi9ezDnnnMPnP/95fOT8usaPH88ee+zBa6+9xnHHHcdxxx3Xpv0ee+wxJk2axL333stTTz3FlClTuPXWWwF44403GDduHBMnTmTRokVExDoBbGP7lCRJ0pbFgNUO/vrXvzJ37lzOPPNMttpqK/r168f48eM36Bi9evXixBNPZNy4cXznO99pp0qrz/Lly5k+fToTJkygZ8+enHfeebz44os8+eSTre47bdo0jjnmGHbffXcGDx7Ml770JaZMmQLAfffdxzbbbMMJJ5xA7969ueCCC5g6deom9ylJkqQtiwGrHfTv35+hQ4cyfvx4fve73/HWW29t9LH22msvnnjiiRKrq25z585tmFZ50EEHUVtby/Dhw5k9e3ar+86ePZtddtmFq6++mosuuojddtutYb/6bQ899BCHH344o0aN4vXXX2fx4sWb1KckSZK2LAasdtCtWzd++9vf0r17d8aOHcugQYO49NJLN+pYNTU1LF++vOQKq9fKlSupqalhxYoVzJkzhyVLltCnTx/q6uravO+CBQuYM2fOOvvVb1u8eDGzZ8+mR48eANTV1W1Sn5IkSdqydO/sArqqnXfemZ/85CesXbuW++67j3HjxvGRj3yEI444YoOOU1dXR9++fdupyurTu3dv6urqGDJkCK+++ioAK1asoKamps37XnHFFQD893//d8N+9duOPvpojj76aJYsWQJUAu6m9ClJkqQtiyNY7axbt24ccsghjBkzhlmzZm3w/k888QR77bVXO1RWnXbeeWdWrVpFbW0tAKtXr+b5559n1KhRre47cuRInn322YblWbNmNezX3Lb+/fszcODATepTkiRJWxYDVjtYs2YNl1xyScM/yP/85z/z4IMPsvfee7f5GG+88QY/+clPmDZtGmeddVZ7lVp1+vbty2GHHcbkyZNZtWoVV111FUOHDmXPPfdcp92YMWO48MIL11k3btw4fv7znzNz5kxeeuklbrjhhoanAR588MEsW7aMW265hZUrV3LFFVdw7LHHblCfkiRJkgGrHXTr1o25c+ey3377UVNTwzHHHMPFF1/MwQcf3Kb9zznnHAYOHMh//ud/8rOf/YzRo0e3c8XV5brrruOpp55i2223ZcqUKUydOpWIWKfNvHnzWLhw4Trr9ttvPyZOnMgnPvEJ9thjD4499tiGgNWrVy+mTZvGpEmTGDhwIGvXrmXy5Mkb1KckSZIUmdnZNZRu9OjR6XdHSZLWZ9iEu5g3+TOdtr8kafMXEY9n5gaPdDiCJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklaRqAlZEHB4RsyNibkRM6Ox6JEmSJKmpqghYEbEV8F/Ap4HdgC9GxG6dW5UkSZIkrasqAhawLzA3M/+SmauBKcBRnVyTJEmSJK2je2cX0EaDgQWNlmuB/Ro3iIjxwPhi8a2IeLqDapPqDQD+2tlFaIvjdbeR4vLO3b+Kec2pM3jdqTOM2pidqiVgRTPrcp2FzOuB6wEiYkZmju6IwqR6XnfqDF536mhec+oMXnfqDBExY2P2q5YpgrXAkEbLOwAvd1ItkiRJktSsaglYfwBGRMROEbE1cDxwRyfXJEmSJEnrqIopgpn5TkScA9wDbAXcmJkzW9jl+o6pTFqH1506g9edOprXnDqD1506w0Zdd5GZrbeSJEmSJLWqWqYISpIkSdJmz4AlSZIkSSWp6oAVEYdHxOyImBsRE5rZ/t6ImFpsfzQihnV8lepq2nDd/WNEzIqIJyPitxGxY2fUqa6jtWuuUbsvRERGhI8y1iZry3UXEccWf97NjIhbOrpGdT1t+Dt2aETcFxF/Kv6ePaIz6lTXERE3RsSi9X2HblRcU1yTT0bE3q0ds2oDVkRsBfwX8GlgN+CLEbFbk2ZnAksyc2fgKmDL/VpIlaKN192fgNGZ+XfAbcC3OrZKdSVtvOaIiD7AucCjHVuhuqK2XHcRMQK4CPhYZu4OfLXDC1WX0sY/7y4Gbs3MD1N5qvR3OrZKdUE/BA5vYfungRHFazzw3dYOWLUBC9gXmJuZf8nM1cAU4KgmbY4Cbire3wYcEhHNfWmx1FatXneZeV9mvlEsPkLle9ukjdWWP+sAvkElzK/qyOLUZbXluvsy8F+ZuQQgMxd1cI3qetpy3SXQt3i/DX4vqjZRZj4AvN5Ck6OAH2XFI0C/iNi+pWNWc8AaDCxotFxbrGu2TWa+AywDtuuQ6tRVteW6a+xM4O52rUhdXavXXER8GBiSmb/syMLUpbXlz7qRwMiIeDgiHomIlv4HWGqLtlx3k4CTIqIW+BXwlY4pTVuwDf23X3V8D9Z6NDcS1fSZ821pI22INl9TEXESMBr4eLtWpK6uxWsuIrpRmQJ9WkcVpC1CW/6s605lyswYKiP1D0bEHpm5tJ1rU9fVluvui8APM/M/IuIA4Obiulvb/uVpC7XBeaKaR7BqgSGNlnfg3cPEDW0iojuVoeSWhgCl1rTluiMiDgX+N/C5zHyrg2pT19TaNdcH2AO4PyLmAfsDd/igC22itv4de3tmvp2ZLwCzqQQuaWO15bo7E7gVIDN/D/QABnRIddpStenffo1Vc8D6AzAiInaKiK2p3Oh4R5M2dwCnFu+/ANybfrOyNk2r110xXes6KuHKexK0qVq85jJzWWYOyMxhmTmMyn1/n8vMGZ1TrrqItvwd+wvgEwARMYDKlMG/dGiV6mract3NBw4BiIhdqQSsxR1apbY0dwCnFE8T3B9YlpmvtLRD1U4RzMx3IuIc4B5gK+DGzJwZEf8CzMjMO4AbqAwdz6UycnV851WsrqCN192/AzXAtOKZKvMz83OdVrSqWhuvOalUbbzu7gEOi4hZwBrgnzLztc6rWtWujdfd+cD3I+I8KtO0TvM/z7UpIuKnVKY6Dyju7ZsIvAcgM79H5V6/I4C5wBvA6a0e02tSkiRJkspRzVMEJUmSJGmzYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRt8SJiWEQ83dl1SJKqnwFLkiRJkkpiwJIkVZVitOnZiLgpIp6MiNsioleTNlMj4ohGyz+MiM8X+z4YEX8sXh9t5vinRcS1jZZ/GRFjiveHRcTvi32nRURNO56qJKkKGbAkSdVoFHB9Zv4dsBw4q8n2KcBxABGxNXAI8CtgEfDJzNy72H5NWzuMiAHAxcChxf4zgH/cxPOQJHUxBixJUjVakJkPF+9/DBzYZPvdwMER8V7g08ADmfkm8B7g+xHxFDAN2G0D+ty/aP9wRDwBnArsuAnnIEnqgrp3dgGSJG2EbLK8TRF6AC7JzDsi4n7gU1RGqn5abDsPWAh8iMp/Mq5q5tjvsO5/QPYofgbw68z84qaXL0nqqhzBkiRVo6ERcUDx/ovALzNzr+J1R7F+CnA6cBBwT7FuG+CVzFwLnAxs1cyx5wF7RUS3iBgC7FusfwT4WETsDBARvSJiZNknJkmqbgYsSVI1egY4NSKeBPoD322mzXTg74HfZObqYt13iv0eAUYCK5vZ72HgBeAp4ArgjwCZuRg4Dfhp0e8jwC5lnZAkqWuIzKazLCRJ2nxFxDAqI1Z7dHIpkiS9iyNYkiRJklQSR7AkSZIkqSSOYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliR1URExMyLGdHYd1Swi7o+IL23kvj+MiG+2sL0uIj7YtG1EHBQRszeu4vapT5LUdgYsSapCETEvIg5tsu60iHiofjkzd8/M+1s5zrCIyIjo3k6laj0ysyYz/9LM+gczc1T9cnOfdUdoXF9E/FNEPB0RKyLihYj4p46uR5KqhQFLktRuNufgtjnXthkK4BRgW+Bw4JyIOL5zS5KkzZMBS5K6qMYjHxGxb0TMiIjlEbEwIq4smj1Q/FxaTAk7ICK6RcTFEfFiRCyKiB9FxDaNjntKse21iPjnJv1MiojbIuLHEbEcOK3o+/cRsTQiXomIayNi60bHy4g4KyLmFCMk34iI4cU+yyPi1sbtWznnH0bE9yLi18WxfhcROzbp6+yImAPMKdZ9NCL+EBHLip8fbXLY4RHxWLH99ojo3+h40yLi1WLbAxGxe5N9B7RSy87NnMOYiKgt3t8MDAXuLD6fr0XEXRHxlSb7PBkRY9fzOzkwIv6n+P0viIjTGm3etjjeioh4NCKGN1dfZn4rM/+Yme9k5mzgduBjzfUnSVs6A5YkbRmuBq7OzL7AcODWYv3fFz/7FVPCfg+cVrw+AXwQqAGuBYiI3YDvACcC2wPbAIOb9HUUcBvQD/gJsAY4DxgAHAAcApzVZJ/DgX2A/YGvAdcXfQwB9gC+uAHneiLwjaK/J4oaGhsL7AfsVoSlu4BrgO2AK4G7ImK7Ru1PAc4APgC8U7StdzcwAngf8Mdm+mqtlhZl5snAfODI4vP5FnATcFJ9m4j4EJXP4FdN94+IoUWN/wkMBPYq6qj3ReBSKiNTc4F/ba2miAjgIGDmhpyLJG0pDFiSVL1+UYxKLI2IpVSCz/q8DewcEQMysy4zH2mh7YnAlZn5l8ysAy4Cji+m1H0BuDMzH8rM1cAlQDbZ//eZ+YvMXJuZb2bm45n5SDH6MQ+4Dvh4k30uz8zlmTkTeBqYXvS/jEpA+HDbfiUA3JWZD2TmW8D/Bg6IiCGNtv9bZr6emW8CnwHmZObNRX0/BZ4FjmzU/ubMfDozVwL/DBwbEVsBZOaNmbmi6GsS8KHGo31tqGVj3A6MiIgRxfLJwNTi82jqROA3mfnTzHw7M1/LzMYB6+eZ+VhmvkMl/O3Vhv4nUfn3w//Z+FOQpK7LgCVJ1WtsZvarf/HuUaHGzgRGAs8W0+A+20LbDwAvNlp+EegODCq2LajfkJlvAK812X9B44WIGBkRvyym0i0HLqMyotPYwkbv32xmuaaFeptqXF8d8HpRd3P1NT1XiuXGo3ILmmx7D5Wpf1tFxOSIeL44r3lFmwHN7bueWjZYEdZuBU6KiG5URqFuXk/zIcDzLRzu1Ubv36CV33NEnENlRO8zRR2SpCYMWJK0BcjMOZn5RSpT2S4HbouI3rx79AngZWDHRstDqUyNWwi8AuxQvyEielKZWrdOd02Wv0tlVGhEMUXx61QemtBeGkaIIqIG6E/lnJqrr+m5QuV8X2rueMW2t4G/AidQmQ55KJWpksPqu92AWtqiuc/oJiqjU4cAbxRTO5uzgMqU0E0WEWcAE4BDMrO2jGNKUldkwJKkLUBEnBQRAzNzLbC0WL0GWAyspXKvVb2fAudFxE5FKLiMyhS0d6jcW3Vk8WCIrancv9NaWOoDLAfqImIX4B828VwyWv5+ryOKBztsTeX+p0czc8F62v4KGBkRJ0RE94g4DtgN+GWjNidFxG4R0Qv4F+C2zFxTnNdbVEbwelH5PW1KLeuzkHU/H4pAtRb4D9Y/egWVaX+HRsSxxfltFxFtmQa4jog4kcr5fbK5R8tLkv7GgCVJW4bDgZkRUUflgRfHZ+aqYorfvwIPF/dy7Q/cSOUf7Q8ALwCrgK8AFPdIfQWYQmU0awWwiErQWJ8LqIz2rAC+D0zd2JOIiB2AOuCpFprdAkykMh1vHyojPc3KzNeAzwLnUwlKXwM+m5l/bdTsZuCHVKbT9QDOLdb/iMqUwZeAWUBz97W1uZYW/BtwcfH5XNBo/Y+APYEfr2/HzJwPHEHl/F6n8oCLD21EDd+kMlL5h+JphnUR8b2NOI4kdXmR2dzMA0mSWleMcC2lMv3vhQ7o7yRg98y8aD3bfwjUZubF7V1LZ4uIU4DxmXlgZ9ciSfobv2RRkrRBIuJI4LdUpgZeQWU0aV5H9J2Z6x2t2ZIU0xXPouUnR0qSOoFTBCVJG+ooKg9qeJnKd0Adn06H6DAR8Skq984tpDIFUZK0GXGKoCRJkiSVxBEsSZIkSSpJl7wHa8CAATls2LDOLkOSJEn6f+zde7yVZZ3//9dHyThsEQkiQ2g7iHhqMiMPpRNhmYdMtFCzzFPx+02W3xz9jdj0FZxpHJqcTKep0dJvViqINamZX7HUSksNO6igCCbC9gQeOGwUUfj8/lj33m02m80C7rX3XpvX8/FYj73u677udX/WXvdD95vruq+lOvXggw++kJlDN/e4XhmwGhsbmT17dneXIUmSJKlORcRTW3KcUwQlSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJLULGBFxNURsSQiHmnTNjgi7oiI+cXPnYv2iIjLI2JBRDwUEfu3OebUov/8iDi1VvVKkiRJ0taq5QjW94Ej2rVNBn6ZmaOBXxbbAEcCo4vHJOA7UAlkwBTgQOAAYEpLKJMkSZKknqZmASszfw281K75WOCa4vk1wIQ27T/IivuAQRGxC/AR4I7MfCkzXwbuYMPQJkmSJEk9Qp8uPt+wzHwWIDOfjYi3Fu3DgcVt+jUVbRtr30BETKIy+sXIkSNLLluS1Ju8f9qdPL3s1S0+fvigftw7eXyJFUmSeouuDlgbEx20ZSftGzZmXglcCTB27NgO+0iSBPD0sldZOO3oLT6+cfKtJVYjSepNunoVweeLqX8UP5cU7U3AiDb9dgWe6aRdkiRJknqcrg5YNwMtKwGeCtzUpv0zxWqCBwHLi6mEtwOHR8TOxeIWhxdtkiRJktTj1GyKYERcD4wDhkREE5XVAKcBN0TEmcAiYGLR/efAUcAC4BXgdIDMfCki/gX4fdHvnzOz/cIZkiRJktQj1CxgZeYnN7LrsA76JnDWRl7nauDqEkuTJEmSpJro6imCkiRJktRrGbAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKkm3BKyIOCci5kTEIxFxfUT0jYjdIuL+iJgfETMiYoei75uL7QXF/sbuqFmSJEmSNqXLA1ZEDAfOBsZm5r7A9sBJwNeASzNzNPAycGZxyJnAy5m5O3Bp0U+SJEmSepzumiLYB+gXEX2A/sCzwHjgxmL/NcCE4vmxxTbF/sMiIrqwVkmSJEmqSpcHrMx8GrgEWEQlWC0HHgSWZeYbRbcmYHjxfDiwuDj2jaL/W9q/bkRMiojZETF76dKltX0TkiRJktSB7pgiuDOVUandgLcDA4AjO+iaLYd0su+vDZlXZubYzBw7dOjQssqVJEmSpKp1xxTBDwFPZubSzHwd+AnwPmBQMWUQYFfgmeJ5EzACoNi/E/BS15YsSZIkSZvWHQFrEXBQRPQv7qU6DJgL3AV8ouhzKnBT8fzmYpti/52ZucEIliRJkiR1t+64B+t+KotV/AF4uKjhSuB84B8iYgGVe6yuKg65CnhL0f4PwOSurlmSJEmSqtFn013Kl5lTgCntmv8CHNBB39XAxK6oS5IkSZK2Rnct0y5JkiRJvY4BS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJK0mnAiog9I+KwiGho135EbcuSJEmSpPqz0YAVEWcDNwFfBB6JiGPb7L641oVJkiRJUr3p08m+zwHvyczmiGgEboyIxsy8DIiuKE6SJEmS6klnAWv7zGwGyMyFETGOSsh6BwYsSZIkSdpAZ/dgPRcR+7VsFGHro8AQ4J21LkySJEmS6k1nAeszwHNtGzLzjcz8DPB3Na1KkiRJkurQRqcIZmZTVBwIDAcSeAZ4IDPv7aoCJUmSJKlebDRgRcThwLeB+cDTRfOuwO4R8fnMnNUF9UmSJElS3ehskYvLgA9l5sK2jRGxG/BzYK8a1iVJkiRJdaeze7D6AE0dtD8NvKk25UiSJElS/epsBOtq4PcRMR1YXLSNAE4Crqp1YZIkSZJUbzpb5OLfIuKnwLHAwVS++6oJ+FRmzu2i+iRJkiSpbnQ2gkVmPgo82kW1SJIkSVJd6+werI2KiKkl1yFJkiRJdW+LAhbwYKlVSJIkSVIvsEUBKzNvKbsQSZIkSap3nX3RcB/gTOA44O1AAs8ANwFXZebrXVKhJEmSJNWJzkawfgjsB0wFjgKOBi4C3gX8aGtOGhGDIuLGiHgsIh6NiIMjYnBE3BER84ufOxd9IyIuj4gFEfFQROy/NeeWJEmSpFrpbBXB/TNzTLu2JuC+iHh8K897GfB/M/MTEbED0B/4MvDLzJwWEZOBycD5wJHA6OJxIPCd4qckSZIk9SidjWC9HBETI6K1T0RsFxEnAi9v6QkjYiDwdxRfVpyZazJzGZXv27qm6HYNMKF4fizwg6y4DxgUEbts6fklSZIkqVY6C1gnAZ8Ano+Ix4tRq+eA44t9W+pvgKXA/4mIP0bE9yJiADAsM58FKH6+teg/HFjc5vimom09ETEpImZHxOylS5duRXmSJEmStGU2OkUwMxcCJwJExFuAyMwXSjrn/sAXM/P+iLiMynTAjYmOytugIfNK4EqAsWPHbrBfkiRJkmqtqmXaM/PFksIVVEagmjLz/mL7RiqB6/mWqX/FzyVt+o9oc/yuVFYzlCRJkqQeZUu/aHiLZeZzwOKIaFlA4zBgLnAzcGrRdiqV5eAp2j9TrCZ4ELC8ZSqhJEmSJPUknX4PVma+UaPzfhG4tlhB8C/A6VTC3g0RcSawCJhY9P05lWXiFwCvFH0lSZIkqcfpbJn2+yKiCfi/VJZUX1jWSTPzT8DYDnYd1kHfBM4q69ySJEmSVCudLXIxNiLeQeV7rylWEAAAIABJREFUqL4ZEcOBe4DbgF9l5mtdVKMkSZIk1YVO78HKzKcy878zcwLwPuAW4EPAbyLi1q4oUJIkSZLqRWdTBNeTma8DdxYPihEtSZIkSVJhi1cRzMynyyxEkiRJkupdly/TLkmSJEm91SYDVkTs2xWFSJIkSVK9q2YE678j4oGI+HxEDKp5RZIkSZJUpzYZsDLzEOBTwAhgdkRcFxEfrnllkiRJklRnqroHKzPnA18Bzgc+AFweEY9FxPG1LE6SJEmS6kk192D9bURcCjwKjAeOycy9iueX1rg+SZIkSaob1XwP1reA7wJfzsxXWxoz85mI+ErNKpMkSZKkOlNNwDoKeDUz1wJExHZA38x8JTN/WNPqJEmSJKmOVHMP1i+Afm22+xdtkiRJkqQ2qglYfTOzuWWjeN6/diVJkiRJUn2qJmCtioj9WzYi4j3Aq530lyRJkqRtUjX3YH0JmBkRzxTbuwAn1q4kSZIkSapPmwxYmfn7iNgTGAME8Fhmvl7zyiRJkiSpzlQzggXwXqCx6P/uiCAzf1CzqiRJkiSpDm0yYEXED4FRwJ+AtUVzAgYsSZIkSWqjmhGsscDemZm1LkaSJEmS6lk1qwg+Aryt1oVIkiRJUr2rZgRrCDA3Ih4AXmtpzMyP1awqSZIkSapD1QSsqbUuQpIkSZJ6g2qWaf9VRLwDGJ2Zv4iI/sD2tS9NkiRJkurLJu/BiojPATcCVxRNw4Gf1rIoSZIkSapH1SxycRbwfmAFQGbOB95ay6IkSZIkqR5VE7Bey8w1LRsR0YfK92BJkiRJktqoJmD9KiK+DPSLiA8DM4FbaluWJEmSJNWfagLWZGAp8DDw/wA/B75Sy6IkSZIkqR5Vs4rgOuC7xUOSJEmStBGbDFgR8SQd3HOVmX9Tk4okSZIkqU5V80XDY9s87wtMBAbXphxJkiRJql+bvAcrM19s83g6M78JjO+C2iRJkiSprlQzRXD/NpvbURnR2rFmFUmSJElSnapmiuB/tHn+BrAQOKEm1UiSJElSHatmFcEPdkUhkiRJklTvqpki+A+d7c/Mb5RXjiRJkiTVr2pXEXwvcHOxfQzwa2BxrYqSJEmSpHpUTcAaAuyfmSsBImIqMDMzP1vLwiRJkiSp3mxymXZgJLCmzfYaoLEm1UiSJElSHatmBOuHwAMR8T9AAscBP6hpVZIkSZJUh6pZRfBfI+I24NCi6fTM/GNty5IkSZKk+lPNFEGA/sCKzLwMaIqI3WpYkyRJkiTVpU0GrIiYApwPXFA0vQn4US2LkiRJkqR6VM0I1nHAx4BVAJn5DLDj1p44IraPiD9GxM+K7d0i4v6ImB8RMyJih6L9zcX2gmJ/49aeW5IkSZJqoZqAtSYzk8oCF0TEgJLO/b+AR9tsfw24NDNHAy8DZxbtZwIvZ+buwKVFP0mSJEnqcaoJWDdExBXAoIj4HPAL4Ltbc9KI2BU4GvhesR3AeODGoss1wITi+bHFNsX+w4r+kiRJktSjVLOK4CUR8WFgBTAGuDAz79jK834T+Ef+OtXwLcCyzHyj2G4ChhfPhwOLi1reiIjlRf8X2r5gREwCJgGMHDlyK8uTJEmSpM3XacCKiO2B2zPzQ8DWhqqW1/wosCQzH4yIcS3NHXTNKvb9tSHzSuBKgLFjx26wX5IkSZJqrdOAlZlrI+KViNgpM5eXdM73Ax+LiKOAvsBAKiNagyKiTzGKtSvwTNG/CRhBZXn4PsBOwEsl1SJJkiRJpanmHqzVwMMRcVVEXN7y2NITZuYFmblrZjYCJwF3ZuangLuATxTdTgVuKp7fXGxT7L+zWHRDkiRJknqUTd6DBdxaPGrtfGB6RHwV+CNwVdF+FfDDiFhAZeTqpC6oRZIkSZI220YDVkSMzMxFmXnNxvpsrcy8G7i7eP4X4IAO+qwGJtaqBkmSJEkqS2dTBH/a8iQiftwFtUiSJElSXessYLVdve9val2IJEmSJNW7zgJWbuS5JEmSJKkDnS1y8a6IWEFlJKtf8ZxiOzNzYM2rkyRJkqQ6stGAlZnbd2UhkiRJklTvqvkeLEmSJElSFQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA1YvtM8++3D33Xd3dxmSJEnSNseAVSONjY3069ePhoaG1scvfvGLqo+PCBYsWLBF554zZw7jxo3bomPr3bx58zjiiCMYNGgQjY2NG+y//PLLGTZsGIMHD+aCCy5Yb9/dd9/NmDFjGDBgABMmTGD58uWt+5qamhg3bhz9+/dn//3355FHHqn1W5EkSVIdMmDV0C233EJzc3Pr40Mf+lB3l9Tr9enTh5NOOomvf/3rG+x74IEHmDp1KnfeeScPP/ww06dP54YbbgDglVdeYeLEiUyZMoUlS5YQEesFsEmTJrHvvvvy4osvcuKJJ3LiiSd22XuSJElS/TBgdYPTTjuNM888k0MOOYQdd9yRM844o3XfkUceSUNDAwDvete7aGho4Etf+lLr/oULFxIRzJw5k3e84x3suOOOfO1rXwPg/PPPp6Ghge22226D0bK1a9dy0UUX0djYyLBhwzj77LN5/fXXW/dffPHF7LLLLgwcOJCxY8fy9NNP1/JXUDOjRo3itNNOY7fddttg38yZMzn++OPZZ599GD58OJ/97GeZPn06AHfddRc77bQTJ598MgMGDOC8885jxowZAKxYsYJZs2YxefJk+vXrxznnnMNTTz3FQw891KXvTZIkST2fAaubzJo1i+uuu445c+Ywc+ZMHnjgAQBuu+02mpubAfjzn/9Mc3Mz3/zmNzc4/oYbbuDBBx/k+eef5/DDDwfga1/7Gs3NzYwcOXKD/t/4xje46aabuPfee1mwYAFz587l8ssvByrT6r761a/y29/+lmXLlvHtb3+bN7/5zbV6691m3rx57Lnnnlx22WVccMEF7L333sybN2+9fffccw9HHHEEY8aM4aWXXmLp0qUsWLCAvn370tDQwKGHHkpTUxOjRo1qPVaSJElqYcCqoQkTJjBo0KDWx5IlS1r3HXPMMYwcOZKRI0fyzne+k8cff3yzXvvCCy9kyJAh9O/fn3e/+92b7P+9732PKVOmMHz4cHbccUfOOussfvKTnwCV+73Wrl3LY489xtq1aznggAMYMmTI5r3ZOrBq1SoaGhpYvHgx8+fPZ8cdd2wNsy37li5dyrx58+jbty8Azc3NrftWrlzJ/Pnzefnll9c7VpIkSWphwKqhn/70pyxbtqz18da3vrV13+DBg1uf77DDDqxevXqzXnv06NGb1X/x4sWccsoprWHv1FNPZenSpQDsscceXHHFFXz1q19l8ODBTJw4kRUrVmzW69eDAQMG0NzczCWXXMKNN97IypUrW6djtuw77rjjePLJJ1unTzY0NLTuGzFiBM899xzvec971jtWkiRJamHA6qEiotP9ffr02azXGzFiBLfffntr2FuxYsV6o2annXYa9957L0888QSPP/44V1999RbV3ZPtsccePPbYY63bc+fOZcyYMRvdN3jwYIYOHcruu+/O6tWraWpqAmDNmjU88cQTrcdKkiRJLQxYPdTb3va2UpcCP+OMM7jwwgt55plnyEwef/xxZs2aBcATTzzBnXfeyZo1a+jTpw/r1q1j4MCBpZ27K2Umq1ev5vXXX299vmbNGgAmTpzIT37yE+bMmcPTTz/NVVdd1boa4Pjx41m+fDnXXXcdq1at4pJLLuGEE04AYODAgRx++OFMmzaN1atXc+mll7ZO7ZQkSZLaMmDV0DHHHLPe92Bdc801VR978cUX84UvfIHhw4dv8H1NHVm7dm3reRYtWtR67quuugqAc889l0MPPZRDDjmEgQMHcvzxx7feE7ZmzRomT57MW97yFvbYYw8OPvhgTjnllC17093sqaeeol+/fhx11FEsWrSIfv36tS4CcuCBBzJlyhQ++MEPsu+++3LCCSe0Bqz+/fszc+ZMpk6dytChQ1m3bh3Tpk1rfd0rrriChx9+mJ133pnp06czY8aMTY4ySpIkadsTmdndNZRu7NixOXv27O4uQ5LUQzVOvpWF047utuMlST1fRDyYmWM39zhHsCRJkiSpJF0esCJiRETcFRGPRsSciPhfRfvgiLgjIuYXP3cu2iMiLo+IBRHxUETs39U1S5IkSVI1umME6w3g3MzcCzgIOCsi9gYmA7/MzNHAL4ttgCOB0cVjEvCdri9ZkiRJkjatywNWZj6bmX8onq8EHgWGA8cCLatAXANMKJ4fC/wgK+4DBkXELl1ctiRJkiRtUrfegxURjcC7gfuBYZn5LFRCGNDyrbzDgcVtDmsq2tq/1qSImB0Rs1u+QFeSJEmSulK3BayIaAB+DHwpM1d01rWDtg2WPszMKzNzbGaOHTp0aFllbpGpU6cSEdx2220ALFu2jB122IFx48Z1a129RVNTE+PGjaN///7sv//+VX9f2N13381222233tL5jz766Hr7x4wZw4ABA5gwYQLLly/f6nNKkiRp29ItASsi3kQlXF2bmT8pmp9vmfpX/FxStDcBI9ocvivwTFfVuqVGjx7NjTfeCMDNN99MY2Nj9xbUi0yaNIl9992XF198kRNPPLH1u6yq8fa3v53m5ubWx1577QXAK6+8wsSJE5kyZQpLliwhItb7/rGtOackSZK2Hd2ximAAVwGPZuY32uy6GTi1eH4qcFOb9s8UqwkeBCxvmUrYkx1wwAH88Y9/5I033uDHP/4xH//4x1v3rV27losuuojGxkaGDRvG2Wefzeuvvw7Ak08+yfjx4xk8eDC77LILp512Gs3NzQAsXLiQiODSSy9lyJAh7Lbbbtx///3d8v66y4oVK5g1axaTJ0+mX79+nHPOOTz11FM89NBDW/W6d911FzvttBMnn3wyAwYM4LzzzmPGjBk1PackSZJ6n+4YwXo/cAowPiL+VDyOAqYBH46I+cCHi22AnwN/ARYA3wU+3w01b5EPfOAD3HzzzTz//POMGjWqtf0b3/gGN910E/feey8LFixg7ty5XH755QC89tprfO5zn2PRokXMmzePRYsWceGFF673uitWrOC5557j2GOP5aKLLurS99TdFixYQN++fWloaODQQw+lqamJUaNGMW/evKqOX7JkCcOGDWPUqFFcfPHFtHzR9rx589hzzz255557OOKIIxgzZgwvvfQSS5cu3epzSpIkadvRp6tPmJn30PF9VQCHddA/gbNqWlSNfOITn+Ckk05i0qRJ67V/73vf49///d8ZPryyVsdZZ53FJZdcwrnnnsuee+7Jnnvu2dp34sSJzJw5c73j//7v/54+ffpw1FFH8bOf/az2b6QHWbVqFQ0NDaxcuZL58+fz8ssvs+OOO7aO8nVmr732Ys6cOYwaNYqHHnqIY489ll122YXTTz+99XWXLl3KvHnz6Nu3LwDNzc1bdU5JkiRtW7o8YG1L3ve+97HXXnsxceJE7rnnntb2xYsXc8opp7DddpUBxHXr1vG2t70NqIywnH322fzmN79h1apVrFmzhve85z3rve7gwYMB2GGHHVi9enUXvZueYcCAATQ3NzNixAiee+45AFauXElDQ8Mmjx02bBjDhg0DYL/99uOss87illtu4fTTT2993eOOO47jjjuOl19+GYCGhoatOqckSZK2Ld26THtvFxHMmjVrvREpgBEjRnD77bezbNkyli1bxooVK3j88ccBuOCCC4gI5s6dy7Jly5g2bVrrNDbB7rvvzurVq2lqagJgzZo1PPHEE4wZM2azX6sl4ALssccePPbYY63bc+fOZfDgwQwdOrTUc0qSJKl3M2B1gzPOOIMLL7yQZ555hszk8ccfZ9asWUBlZGTgwIE0NDTw5JNPcsUVV3RztT3LwIEDOfzww5k2bRqrV6/m0ksvZeTIkbzzne9cr9+4ceM4//zz12u78847WbRoEQCPPvoo3/nOdzjmmGMAGD9+PMuXL+e6665j1apVXHLJJZxwwgmbdU5JkiTJgNUNzj33XA499FAOOeQQBg4cyPHHH8+SJZVV6adMmcLs2bPZaaedOOGEE1oDgP7qiiuu4OGHH2bnnXdm+vTpzJgxg8rilH+1cOFCnn/++fXa/vCHP/De976X/v37c9RRR/G5z32O008/HYD+/fszc+ZMpk6dytChQ1m3bh3Tpk3brHNKkiRJ0Runn40dOzZnz57d3WVIknqoxsm3snDa0d12vCSp54uIBzNz7OYe5wiWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIBVI3PmzOGwww5j0KBBvPWtb+XjH/84ABdffDENDQ00NDSw/fbb07dv39btRYsWMW7cOPr27cuAAQM4+OCDuf3227v5nfQ8TU1NjBs3jv79+7P//vvzyCOPVH3s5ZdfzrBhwxg8eDAXXHDBevvuvvtuxowZw4ABA5gwYQLLly8v5ZySJEnadhiwamTChAkcffTRvPjiizz++ON85CMfAeDLX/4yzc3NNDc3c+ihh/Ktb32rdXvkyJEAfOtb32Lp0qV84Qtf4OMf/zguOb++SZMmse+++/Liiy9y4okncuKJJ1Z13AMPPMDUqVO58847efjhh5k+fTo33HADAK+88goTJ05kypQpLFmyhIhYL4Bt6TklSZK0bTFg1cALL7zAggULOPPMM9l+++0ZNGgQkyZN2qzX6N+/P5/61KeYOHEi3/72t2tUaf1ZsWIFs2bNYvLkyfTr149zzjmHp556ioceemiTx86cOZPjjz+effbZh+HDh/PZz36W6dOnA3DXXXex0047cfLJJzNgwADOO+88ZsyYsdXnlCRJ0rbFgFUDgwcPZuTIkUyaNIlf/epXvPbaa1v8Wvvttx9/+tOfSqyuvi1YsKB1WuWhhx5KU1MTo0aNYt68eZs8dt68eey5555cdtllXHDBBey9996tx7Xsu+eeezjiiCMYM2YML730EkuXLt2qc0qSJGnbYsCqge22245f/vKX9OnThwkTJjBs2DAuuuiiLXqthoYGVqxYUXKF9WvVqlU0NDSwcuVK5s+fz8svv8yOO+5Ic3Nz1ccuXryY+fPnr3dcy76lS5cyb948+vbtC0Bzc/NWnVOSJEnblj7dXUBvtfvuu3Pttdeybt067rrrLiZOnMh73/tejjrqqM16nebmZgYOHFijKuvPgAEDaG5uZsSIETz33HMArFy5koaGhqqPveSSSwD4n//5n9bjWvYdd9xxHHfccbz88stAJeBuzTklSZK0bXEEq8a22247DjvsMMaNG8fcuXM3+/g//elP7LfffjWorD7tvvvurF69mqamJgDWrFnDE088wZgxYzZ57B577MFjjz3Wuj137tzW4zraN3jwYIYOHbpV55QkSdK2xYBVA2vXruXCCy9s/YP8z3/+M7/5zW/Yf//9q36NV155hWuvvZaZM2fy+c9/vlal1p2BAwdy+OGHM23aNFavXs2ll17KyJEjeec737lev3HjxnH++eev1zZx4kR+8pOfMGfOHJ5++mmuuuqq1tUAx48fz/Lly7nuuutYtWoVl1xyCSeccMJmnVOSJEkyYNXAdtttx4IFCzjwwANpaGjg+OOP5ytf+Qrjx4+v6vgvfOELDB06lP/8z//kxz/+MWPHjq1xxfXliiuu4OGHH2bnnXdm+vTpzJgxg4hYr8/ChQt5/vnn12s78MADmTJlCh/84AfZd999OeGEE1oDVv/+/Zk5cyZTp05l6NChrFu3jmnTpm3WOSVJkqTIzO6uoXRjx45NvztKkrQxjZNvZeG0o7vteElSzxcRD2bmZo90OIIlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJIkSZJUEgOWJEmSJJXEgCVJkiRJJTFgSZIkSVJJDFiSJEmSVJK6CVgRcUREzIuIBRExubvrkSRJkqT26iJgRcT2wH8BRwJ7A5+MiL27typJkiRJWl9dBCzgAGBBZv4lM9cA04Fju7kmSZIkSVpPn+4uoErDgcVttpuAA9t2iIhJwKRi87WIeKSLapNaDAFe6O4itM3xuttC8bXuPb6Oec2pO3jdqTuM2ZKD6iVgRQdtud5G5pXAlQARMTszx3ZFYVILrzt1B687dTWvOXUHrzt1h4iYvSXH1csUwSZgRJvtXYFnuqkWSZIkSepQvQSs3wOjI2K3iNgBOAm4uZtrkiRJkqT11MUUwcx8IyK+ANwObA9cnZlzOjnkyq6pTFqP1526g9eduprXnLqD1526wxZdd5GZm+4lSZIkSdqkepkiKEmSJEk9ngFLkiRJkkpS1wErIo6IiHkRsSAiJnew/80RMaPYf39ENHZ9leptqrju/iEi5kbEQxHxy4h4R3fUqd5jU9dcm36fiIiMCJcy1lar5rqLiBOK/97NiYjrurpG9T5V/D92ZETcFRF/LP4/e1R31KneIyKujoglG/sO3ai4vLgmH4qI/Tf1mnUbsCJie+C/gCOBvYFPRsTe7bqdCbycmbsDlwLb7tdCqhRVXnd/BMZm5t8CNwL/3rVVqjep8pojInYEzgbu79oK1RtVc91FxGjgAuD9mbkP8KUuL1S9SpX/vfsKcENmvpvKqtLf7toq1Qt9Hziik/1HAqOLxyTgO5t6wboNWMABwILM/EtmrgGmA8e263MscE3x/EbgsIjo6EuLpWpt8rrLzLsy85Vi8z4q39smbalq/lsH8C9UwvzqrixOvVY1193ngP/KzJcBMnNJF9eo3qea6y6BgcXznfB7UbWVMvPXwEuddDkW+EFW3AcMiohdOnvNeg5Yw4HFbbabirYO+2TmG8By4C1dUp16q2quu7bOBG6raUXq7TZ5zUXEu4ERmfmzrixMvVo1/63bA9gjIu6NiPsiorN/AZaqUc11NxX4dEQ0AT8Hvtg1pWkbtrl/+9XH92BtREcjUe3XnK+mj7Q5qr6mIuLTwFjgAzWtSL1dp9dcRGxHZQr0aV1VkLYJ1fy3rg+VKTPjqIzU/yYi9s3MZTWuTb1XNdfdJ4HvZ+Z/RMTBwA+L625d7cvTNmqz80Q9j2A1ASPabO/KhsPErX0iog+VoeTOhgClTanmuiMiPgT8E/CxzHyti2pT77Spa25HYF/g7ohYCBwE3OxCF9pK1f4/9qbMfD0znwTmUQlc0paq5ro7E7gBIDN/B/QFhnRJddpWVfW3X1v1HLB+D4yOiN0iYgcqNzre3K7PzcCpxfNPAHem36ysrbPJ666YrnUFlXDlPQnaWp1ec5m5PDOHZGZjZjZSue/vY5k5u3vKVS9Rzf9jfwp8ECAihlCZMviXLq1SvU01190i4DCAiNiLSsBa2qVValtzM/CZYjXBg4DlmflsZwfU7RTBzHwjIr4A3A5sD1ydmXMi4p+B2Zl5M3AVlaHjBVRGrk7qvorVG1R53X0daABmFmuqLMrMj3Vb0aprVV5zUqmqvO5uBw6PiLnAWuD/y8wXu69q1bsqr7tzge9GxDlUpmmd5j+ea2tExPVUpjoPKe7tmwK8CSAz/5vKvX5HAQuAV4DTN/maXpOSJEmSVI56niIoSZIkST2KAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkrTNi4jGiHiku+uQJNU/A5YkSZIklcSAJUmqK8Vo02MRcU1EPBQRN0ZE/3Z9ZkTEUW22vx8RHy+O/U1E/KF4vK+D1z8tIr7VZvtnETGueH54RPyuOHZmRDTU8K1KkuqQAUuSVI/GAFdm5t8CK4DPt9s/HTgRICJ2AA4Dfg4sAT6cmfsX+y+v9oQRMQT4CvCh4vjZwD9s5fuQJPUyBixJUj1anJn3Fs9/BBzSbv9twPiIeDNwJPDrzHwVeBPw3Yh4GJgJ7L0Z5zyo6H9vRPwJOBV4x1a8B0lSL9SnuwuQJGkLZLvtnYrQA3BhZt4cEXcDH6EyUnV9se8c4HngXVT+kXF1B6/9Buv/A2Tf4mcAd2TmJ7e+fElSb+UIliSpHo2MiIOL558EfpaZ+xWPm4v26cDpwKHA7UXbTsCzmbkOOAXYvoPXXgjsFxHbRcQI4ICi/T7g/RGxO0BE9I+IPcp+Y5Kk+mbAkiTVo0eBUyPiIWAw8J0O+swC/g74RWauKdq+XRx3H7AHsKqD4+4FngQeBi4B/gCQmUuB04Dri/PeB+xZ1huSJPUOkdl+loUkST1XRDRSGbHat5tLkSRpA45gSZIkSVJJHMGSJEmSpJI4giVJkiRJJTFgSZIkSVJJDFiSJEmSVBIDliRJkiSVxIAlSZIkSSUxYEmSJElSSQxYkiRJklQSA5YkSZIklcSAJUmSJEklMWBJkiRJUkkMWJLUTSJiTkSM6+466llE3B0Rn93CY78fEV/tZH9zRPxN+74RcWhEzNuyissTEY0RkRHRp7trkST9lQFLkmogIhZGxIfatZ0WEfe0bGfmPpl59yZexz+iu0lmNmTmXzpo/01mjmnZ7uizrlb7zzcq/jMiHouI4VtevSSpuxiwJGkb1pODW0+urRYiIoArgHHABzLz6e6taOtta5+hJIEBS5K6TduRj4g4ICJmR8SKiHg+Ir5RdPv/j8knAAAgAElEQVR18XNZMWXt4IjYLiK+EhFPRcSSiPhBROzU5nU/U+x7MSL+d7vzTI2IGyPiRxGxAjitOPfvImJZRDwbEd+KiB3avF5GxOcjYn5ErIyIf4mIUcUxKyLihrb9N/Gevx8R/x0RdxSv9auIeEe7c50VEfOB+UXb+yLi9xGxvPj5vnYvOyoiHij23xQRg9u83syIeK7Y9+uI2KfdsUM2UcvuHbyHcRHRVDz/ITASuKX4fP4xIm6NiC+2O+ahiJjQya9me+D7wFhgXGY+XxzX6Wfd7hx3R8RXI+K3RS23RMRbIuLa4nP6fUQ0tum/0d9rROxW/L5WRsQvIuK/IuJHbfZ/LCpTXJcV592rzb6FEXF+RDwErIqIPhGxV9FvWXHcx4q+BxWfz/Ztjj+uOFaS6pIBS5J6hsuAyzJzIDAKuKFo/7vi56BiytrvgNOKxweBvwEagG8BRMTewLeBTwG7ADsB7aeaHQvcCAwCrgXWAucAQ4CDgcOAz7c75gjgPcBBwD8CVxbnGAHsC3xyM97rp4B/Kc73p6KGtiYABwJ7F2HpVuBy4C3AN4BbI+Itbfp/BjgDeDvwRtG3xW3AaOCtwB86ONemaulUZp4CLAKOKT6ffweuAT7d0ici3kXlM/h5Jy91LbAnMD4zX2zTfhob+aw34iTglOJ8o4DfAf8HGAw8CkwpatrU7/U64IFi39TiNVvezx7A9cCXgKHF+7qlXcj+JHA0lWssgFuAWVQ+hy8C10bEmMy8D1gFjG9z7MnF+SWpLhmwJKl2flr8i/2yiFhGJfhszOvA7hExJDObiz88N+ZTwDcy8y+Z2QxcAJwUlelYnwBuycx7MnMNcCGQ7Y7/XWb+NDPXZearmflgZt6XmW9k5kIq09Q+0O6Yr2XmisycAzwCzCrOv5xKiHl3db8SAG7NzF9n5mvAPwEHR8SINvv/LTNfysxXqfyRPj8zf1jUdz3wGHBMm/4/zMxHMnMV8L+BE1pGRDLz6sxcWZxrKvCudiNAm6plS9wEjI6I0cX2KcCM4vPYmMOBGzJzWbv2zj7rjvyfzHyizefyRGb+IjPfAGby189po7/XiBgJvBe4MDPXZOY9wM1tznEild/bHZn5OnAJ0A9oO7J4eWYuLj7Dg6gEw2nF690J/Iy/hvLrW55HxI7AUUWbJNUlA5Yk1c6EzBzU8mDDUaG2zgT2AB4rpmt9tJO+bweearP9FNAHGFbsW9yyIzNfAV5kfYvbbkTEHhHxs2Kq1grgYiojOm093+b5qx1sN3RSb3tt62sGXirq7qi+9u+VYrvtqNzidvveRGXq3/YRMS0inije18Kiz5COjt1ILZutCGs3AJ+OiO2ohIcfbuKwjwJTIuKMdu2dfdYdqfZz6uz3+nbgpeLaabHRzyQz1xX7N/aZvB1YXPRrfy6ojFYdHxFvBo4H/pCZ7WuTpLphwJKkHiAz52fmJ6lMofoacGNEDGDD0SeAZ4B3tNkeSWVq3PPAs8CuLTsioh+VaV7rna7d9neojF6MLqYofpnKtK5aaR0hiogGKtPXntlIfe3fK1Teb9sFIEa02/c68AKVqWbHAh+iMlWyseW0m1FLNTr6jK6hMvp0GPBKMbWzM7+lMip3WUSc3Ka9s896a3T2e30WGBwR/dvsa/s7Xu/YiIhif9vPpP1nOKIIm+3PRWbOpRK4jsTpgZJ6AQOWJPUAEfHpiBha/Ct/yzSxtcBSYB2V+29aXA+cUyxE0EBlxGlGMQ3sRirTvN5X3BNzEZsOSzsCK4DmiNgT+PutfC8ZnX+/11ERcUhR378A92fm4o30/TmwR0ScXCyWcCKwN5UpZi0+HRF7F4Hgn4EbM3Nt8b5eozKC15/K72lratmY51n/86EIVOuA/2DTo1ctx/yKygjOlRHxiaK5s896a2z091qMHs0GpkbEDhFxMOtPybwBODoiDouINwHnUvk9/3Yj57qfyn1W/xgRbyqujWOA6W36XAecTeWew5lb+d4kqVsZsCSpZzgCmBMRzVQWvDgpM1cX07T+Fbi3uJfrIOBqKn+0/xp4ElhNZeEAinukvkjlj9dngZXAEip/AG/MeVRGDlYC3wVmbOmbiIhdgWbg4U66XUdlsYWXqCyc8amNdSwWfPgolT/iX6SywMZHM/OFNt1+SGUFvueAvlT+UAf4AZWRkaeBuUBH97VVXUsn/g34SvH5nNem/QfAO4EfdXzYhjLzDir3OH0/Io6hk896a1Txe/0UlQVPXgS+SuWaeK04dh6VRTz+k8pI4TFUFvno8B6zov1jVEaoXqByL+JnMvOxNt2up7I8/Z3tPltJqjuR2dHMBklSb1CMeiyjMv3vyS4436eBfTLzgo3s/z7QlJlfqXUt3S0iPgNMysxDuruWrRURM4DHMnNKd9ciST2dXwAoSb1MMfLxSypTAy+hMpq0sCvOnZlVj9b0ZsV0xc/T+cqRPVZEvJfKqN6TVFY4PBaY1q1FSVKdcIqgJPU+x1JZWOAZKt8BdVI6XaHLRMRHqNw79zz1u2DD24C7qUz3vBz4+8z8Y7dWJEl1wimCkiRJklQSR7AkSZIkqSTdcg9WRFxNZfWiJZm5b7t95wFfB4Zm5gvF92tcRuWb3V8BTsvMP3T2+kOGDMnGxsaa1C5JkiSp93vwwQdfyMyhm3tcdy1y8X3gW1SWsG0VESOADwOL2jQfSeUegtHAgVS+EPPAzl68sbGR2bNnl1iuJEmSpG1JRDy1Jcd1yxTBzPw1ldWJ2ruUyndxtL0x7FjgB1lxHzAoInbpgjIlSZIkabP0mHuwIuJjwNOZ+ed2u4YDi9tsNxVtkiRJktSj9IjvwSq+L+SfqHzXxga7O2jbYOnDiJgETAIYOXJkqfVJkiRJUjV6ygjWKGA34M8RsRDYFfhDRLyNyojViDZ9d6Xy3S7rycwrM3NsZo4dOnSz70WTJEmSpK3WIwJWZj6cmW/NzMbMbKQSqvbPzOeAm4HPRMVBwPLMfLY765UkSZKkjnRLwIqI64HfAWMioikizuyk+8+BvwALgO8Cn++CEiVJkiRps3XLPViZ+clN7G9s8zyBs2pdkyRJkiRtrR4xRVCSJEmSegMDliRJkiSVxIAlSZIkSSUxYEmSJElSSXrEFw1LkiT1du+fdidPL3t1i48fPqgf904eX2JFkmrBgCVJktQFnl72KgunHb3FxzdOvrXEaiTVilMEJUmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkkhiwJEmSJKkkBixJkiRJKokBS5IkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJIkSZKkknR5wIqIqyNiSUQ80qbt6xHxWEQ8FBH/ExGD2uy7ICIWRMS8iPhIV9crSZIkSdXqjhGs7wNHtGu7A9g3M/8WeBy4ACAi9gZOAvYpjvl2RGzfdaVKkiRJUvW6PGBl5q+Bl9q1zcrMN4rN+4Bdi+fHAtMz87XMfBJYABzQZcVKkiRJ0mboifdgnQHcVjwfDixus6+paNtAREyKiNkRMXvp0qU1LlGSJEmSNtSjAlZE/BPwBnBtS1MH3bKjYzPzyswcm5ljhw4dWqsSJUmSJGmj+nR3AS0i4lTgo8BhmdkSopqAEW267Qo809W1SZIkSVI1esQIVkQcAZwPfCwzX2mz62bgpIh4c0TsBowGHuiOGiVJkiRpU7p8BCsirgfGAUMiogmYQmXVwDcDd0QEwH2Z+f9m5pyIuAGYS2Xq4FmZubara5YkSZKkanR5wMrMT3bQfFUn/f8V+NfaVSRJkiRJ5egRUwQlSZIkqTcwYEmSJElSSQxYkiRJklQSA5YkSZIklaTHfA+WJEmS1JH3T7uTp5e9usXHDx/Uj3snjy+xImnjDFiSJEnq0Z5e9ioLpx29xcc3Tr61xGqkzjlFUJIkSZJKYsCSJEmSpJIYsCRJkiSpJAYsSZIkSSqJAUuSJEmSSmLAkiRJkqSSGLAkSZIkqSQGLEmSJEkqiQFLkiRJkkpiwJKk/7+9e4/Tqqz3///6ACkzDAcRREJwCBUPaGZomZKE5k8xEy00y4pC6WBppm0Hdxvw+/PnHsstRaVhWpplCB4xdUeJqGFqmgYCAqMOOKiACsIgJ+H6/TE38wUZYbhnzeEeXs/H437Mva61rrU+S9dj4M211rUkSZIy0q65C5AkSU3ruPLpLFm5Nq++vboUMbNsSMYVSVLrYcCSJGk3s2TlWirLT8urb2nZAxlXI0mti7cISpIkSVJGDFiSJEmSlBEDliRJkiRlxIAlSZIkSRkxYEmSJElSRgxYkiRJkpQRA5YkSZIkZcSAJUmSJEkZMWBJkiRJUkYMWJIkSZKUEQOWJEmSJGXEgCVJkiRJGWmWgBURv42IZRHxwlZtXSPirxGxMPdzr1x7RMSEiKiIiFkRcVRz1CxJkiRJO9NcI1i3AKe8r60MeDildCDwcG4Z4FTgwNxnFHBDE9UoSZIkSbukWQJWSukx4O33NZ8B3Jr7fiswbKv236caTwJdIqJn01QqSZIkSfXXkp7B6pFSeh0g93OfXHsv4NWttqvKtW0jIkZFxDMR8czy5csbvVhJkiRJer+WFLA+SNTRlrZrSOnGlNLAlNLA7t27N0FZkiRJkrStlhSwlm659S/3c1muvQrovdV2+wGvNXFtkiRJkrRTLSlgTQW+nvv+deC+rdq/lptN8JPAO1tuJZQkSZKklqRdcxw0Iv4EDAa6RUQVMBYoByZHxEhgMTA8t/mDwFCgAngX+EaTFyxJkiRJ9dAsASuldO4HrDqxjm0TcGHjViRJkiRJDdeSbhGUJEmSpIJmwJIkSZKkjDTLLYKSJB1XPp0lK9fm3b9XlyJmlg3JsCJJkhrOgCVJahZLVq6lsvy0vPuXlj2QYTWSJGXDWwQlSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJykheASsiDo6IEyOi5H3tp2RTliRJkiQVnl0OWBFxEXAf8H3ghYg4Y6vVV2dVmCRJkiQVmnxeNHwB8PGUUnVElAJ3RkRpSunnQGRZnCRJkiQVknwCVtuUUjVASqkyIgZTE7L2x4AlSZIktRjHlU9nycq1effv1aWImWVDMqyo9csnYL0REUemlJ4HyI1kfQ74LXB4ptVJkiRJytuSlWupLD8t7/6lZQ9kWM3uIZ9JLr4GvLF1Q0rpvZTS14BPZ1KVJEmSJBWgXR7BSilVRY1PAL2ABLwGPJ1Smpl1gZIkSZJUKHY5YEXEycD1wEJgSa55P+CAiPhuSmlahvVJkiRJUsHI5xmsnwMnpZQqt26MiL7Ag8AhGdQlSZIkSQUnn2ew2gFVdbQvAT7UsHIkSZIkqXDlM4L1W+CfETEJeDXX1hv4EnBzVoVJkiRJUqHJZ5KL/46Ie4EzgGOpefdVFfCVlNLcjOuTJEmSpIKRzwgWKaV5wLyMa5EkSZKkgpbPM1gfKCLGZbk/SZIkSSokmQYs4NmM9ydJkiRJBSPTgJVSuj/L/UmSJElSIcnnRcPtgJHAmcCHgQS8BtwH3JxS2phphZIkSZJUIPKZ5OI2YCUwjv/7Pqz9gK8DfwDOyaQySZIkSSow+QSso1JK/d/XVgU8GRELMqhJkiRJkgpSPs9grYiI4RFR2zci2kTEOcCK7EqTJEmSpMKST8D6EvBFYGlELMiNWr0BnJVbl7eIuCQi5kTECxHxp4hoHxF9I+KpiFgYEXdExB4NOYYkSZIkNZZdvkUwpVRJ7jmriNgbiJTSmw0tJCJ6ARcBh6aU1kbEZGoC21BgfEppUkT8mpoJNm5o6PEkSZIkKWsNmqY9pfRWFuFqK+2AotxMhcXA68AQ4M7c+luBYRkeT5IkSZIyk/WLhvOWUloCXAsspiZYvUPNi4tXppTey21WBfSqq39EjIqIZyLimeXLlzdFyZIkSZK0jV0OWLnRpcxFxF7AGUBfat6v1QE4tY5NU139U0o3ppQGppQGdu/evTFKlCRJkqQdyicsPRkRVcD/Av+beyYrCycBr6SUlgNExN3Ap4AuEdEuN4q1HzUvNZbUghxXPp0lK9fm3b9XlyJmlg3JsCJJkqTmkc8kFwMjYn9qRpd+lpuc4u/AQ8CjKaX1edayGPhkRBQDa4ETgWeAR6iZtXASNS8zvi/P/UtqJEtWrqWy/LS8+5eWPZBhNZIkSc0nr2ewUkqLUkq/TikNo2aU6X5qRqAej4i8/qaUUnqKmsks/gXMztV2I3A58MOIqAD2Bm7OZ/+SJEmS1Nga/DxVSmkjMD332TLder77GguMfV/zy8AxeRcoSZIkSU0k8wkrcrMBStJuwefPJEnS1hplRkBJ2l34/JkkSdpa3u/BiogBWRYiSZIkSYWuIS8a/nVEPB0R342ILplVJEmSJEkFKu+AlVI6HvgK0Bt4JiJuj4jPZlaZJEmSJBWYhoxgkVJaCPyYmqnUTwAmRMSLEXFWFsVJkiRJUiFpyDNYR0TEeGAeMAQ4PaV0SO77+IzqkyRJkqSC0ZBZBH8J/Aa4IqVUO0dxSum1iPhxgyuTJEmSpALTkIA1FFibUtoEEBFtgPYppXdTSrdlUp0kSZIkFZCGPIP1N6Boq+XiXJskSZIk7ZYaErDap5Sqtyzkvhc3vCRJkiRJKkwNCVhrIuKoLQsR8XFg7Q62lyRJkqRWrSHPYP0AmBIRr+WWewLnNLwkSZIkSSpMeQeslNI/I+JgoD8QwIsppY2ZVSZJkiRJBaYhI1gARwOluf18LCJIKf2+wVVJkiRJUgHKO2BFxG1AP+B5YFOuOQEGLEmSJEm7pYaMYA0EDk0ppayKkSRJkqRC1pBZBF8A9s2qEEmSJEkqdA0ZweoGzI2Ip4H1WxpTSp9vcFWSJEmSVIAaErDGZVWEJEmSJLUGDZmm/dGI2B84MKX0t4goBtpmV5okSZIkFZa8n8GKiAuAO4GJuaZewL1ZFCVJkiRJhaghk1xcCBwHrAJIKS0E9smiKEmSJEkqRA0JWOtTShu2LEREO2regyVJkiRJu6WGBKxHI+IKoCgiPgtMAe7PpixJkiRJKjwNCVhlwHJgNvAt4EHgx1kUJUmSJEmFqCGzCG4GfpP7SJIkSdJuL++AFRGvUMczVymljzSoIqlAHVc+nSUr1+bdv1eXImaWDcmwIkmSJDW1hrxoeOBW39sDw4GuDStHKlxLVq6lsvy0vPuXlj2QYTWSJElqDnk/g5VSemurz5KU0s+ABv3ze0R0iYg7I+LFiJgXEcdGRNeI+GtELMz93Kshx5AkSZKkxtKQFw0ftdVnYER8G+jYwHp+DvxvSulg4KPAPGom03g4pXQg8HBuWZIkSZJanIbcIvg/W31/D6gEzs53ZxHRCfg0MAIg946tDRFxBjA4t9mtwAzg8nyPI0mSJEmNpSGzCH4my0KAj1Az7fvvIuKjwLPAxUCPlNLruWO+HhH7ZHxcSZIkScpEQ2YR/OGO1qeUrsujlqOA76eUnoqIn7MLtwNGxChgFECfPn128dCSJEmS1HANedHwQOA7QK/c59vAodQ8h5XPs1hVQFVK6anc8p3UBK6lEdETIPdzWV2dU0o3ppQGppQGdu/ePY/DS5IkSVLDNOQZrG7AUSml1QARMQ6YklI6P5+dpZTeiIhXI6J/Smk+cCIwN/f5OlCe+3lfA2qWJEmSpEbTkIDVB9iw1fIGoLRB1cD3gT9GxB7Ay8A3qBllmxwRI4HF1LxvS5IkSZJanIYErNuApyPiHiABZwK/b0gxKaXn2fYFxluc2JD9SpIkSVJTaMgsgv9fRDwEDMo1fSOl9Fw2ZUmSJElS4WnIJBcAxcCqlNLPgaqI6JtBTZIkSZJUkPIOWBExlpoX/o7ONX0I+EMWRUmSJElSIWrICNaZwOeBNQAppdfIb3p2SZIkSWoVGhKwNqSUEjUTXBARHbIpSZIkSZIKU0MC1uSImAh0iYgLgL8Bv8mmLEmSJEkqPA2ZRfDaiPgssAroD4xJKf01s8okSZIk7baOK5/OkpVr8+7fq0sRM8uGZFhR/eQVsCKiLfCXlNJJgKFKkiRJUqaWrFxLZflpefcvLXsgw2rqL69bBFNKm4B3I6JzxvVIkiRJUsHK+xZBYB0wOyL+Sm4mQYCU0kUNrkqSJEmSClBDAtYDuY8kSZIkiTwCVkT0SSktTind2hgFSZIkSVKhyucZrHu3fImIuzKsRZIkSZIKWj4BK7b6/pGsCpEkSZKkQpdPwEof8F2SJEmSdmv5THLx0YhYRc1IVlHuO7nllFLqlFl1kiRJklRAdjlgpZTaNkYhkiRJklTo8nrRsCRJkiRpewYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIy0uYEVE24h4LiL+nFvuGxFPRcTCiLgjIvZo7holSZIkqS4tLmABFwPztlq+BhifUjoQWAGMbJaqJEmSJGknWlTAioj9gNOAm3LLAQwB7sxtciswrHmqkyRJkqQda1EBC/gZ8B/A5tzy3sDKlNJ7ueUqoFddHSNiVEQ8ExHPLF++vPErlSRJkqT3aTEBKyI+ByxLKT27dXMdm6a6+qeUbkwpDUwpDezevXuj1ChJkiRJO9KuuQvYynHA5yNiKNAe6ETNiFaXiGiXG8XaD3itGWtUIzuufDpLVq7Nq2+vLkXMLBuScUWSJElS/bWYgJVSGg2MBoiIwcBlKaWvRMQU4IvAJODrwH3NVqQa3ZKVa6ksPy2vvqVlD2RcjSRJkrRrWswtgjtwOfDDiKig5pmsm5u5HkmSJEmqU4sZwdpaSmkGMCP3/WXgmOasR5IkSZLqoxBGsCRJkiSpIBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsBqhQ477DBmzJjR3GVIkiRJux0DViMpLS2lqKiIkpKS2s/f/va3evePCCoqKvI69pw5cxg8eHBefQvd/PnzOeWUU+jSpQulpaXbrZ8wYQI9evSga9eujB49ept1M2bMoH///nTo0IFhw4bxzjvv1K6rqqpi8ODBFBcXc9RRR/HCCy809qlIkiSpABmwGtH9999PdXV17eekk05q7pJavXbt2vGlL32Jn/70p9ute/rppxk3bhzTp09n9uzZTJo0icmTJwPw7rvvMnz4cMaOHcuyZcuIiG0C2KhRoxgwYABvvfUW55xzDuecc06TnZMkSZIKhwGrGYwYMYKRI0dy/PHH07FjR775zW/Wrjv11FMpKSkB4KMf/SglJSX84Ac/qF1fWVlJRDBlyhT2339/OnbsyDXXXAPA5ZdfTklJCW3atNlutGzTpk1ceeWVlJaW0qNHDy666CI2btxYu/7qq6+mZ8+edOrUiYEDB7JkyZLG/E/QaPr168eIESPo27fvduumTJnCWWedxWGHHUavXr04//zzmTRpEgCPPPIInTt35stf/jIdOnTgsssu44477gBg1apVTJs2jbKyMoqKirjkkktYtGgRs2bNatJzkyRJUstnwGom06ZN4/bbb2fOnDlMmTKFp59+GoCHHnqI6upqAP79739TXV3Nz372s+36T548mWeffZalS5dy8sknA3DNNddQXV1Nnz59ttv+uuuu47777mPmzJlUVFQwd+5cJkyYANTcVnfVVVfxxBNPsHLlSq6//nr23HPPxjr1ZjN//nwOPvhgfv7znzN69GgOPfRQ5s+fv826v//975xyyin079+ft99+m+XLl1NRUUH79u0pKSlh0KBBVFVV0a9fv9q+kiRJ0hYGrEY0bNgwunTpUvtZtmxZ7brTTz+dPn360KdPHw4//HAWLFiwS/seM2YM3bp1o7i4mI997GM73f6mm25i7Nix9OrVi44dO3LhhRdy9913AzXPe23atIkXX3yRTZs2ccwxx9CtW7ddO9kCsGbNGkpKSnj11VdZuHAhHTt2rA2zW9YtX76c+fPn0759ewCqq6tr161evZqFCxeyYsWKbfpKkiRJW7Rr7gJas3vvvfcDn7vq2rVr7fc99tiDdevW7dK+DzzwwF3a/tVXX+WrX/0qbdrUZOrNmzez7777AnDQQQcxceJErrrqKj437Au073sUe596MW32LN6lYwD06lLEzLIhu9yvKXTo0IHq6mquvfZaAO65557a2zG3rDvzzDM588wzWbFiBQAlJSW163r37s0bb7wBwOrVq2v7SpIkSVsYsFqoiNjh+nbtdu1/Xe/evbnllls49thj61w/YsQIRowYQe/v/5Guj/2Eb+y/ZJtnv+qrtOyBXe7TVA466CBefPHF2uW5c+fSv3//2nXXX3/9Nuu6du1K9+7d2XPPPVm3bh1VVVXst99+bNiwgZdeeqm2ryRJkrSFtwi2UPvuu2+mU4F/85vfZMyYMbz22muklFiwYAHTpk0D4KWXXmL69Ols2LAB2rRl8+bNdOrUKbNjN6WUEuvWrWPjxo213zds2ADA8OHDufvuu5kzZw5Llizh5ptvrp0NcMiQIbzzzjvcfvvtrFmzhmuvvZazzz4bgE6dOnHyySdTXl7OunXrGD9+fO2tnZIkSdLWDFiN6PTTT9/mPVi33nprvfteffXVfO9736NXr17bva+pLps2bao9zuLFi2uPffPNNwNw6aWXMmjQII4//ng6derEWWedVftM2IYNGygrK2Pvvffmtd98i2OPPZavfvWr+Z10M1u0aBFFRUUMHTqUxYsXU1RUVDsJyCc+8QnGjh3LZz7zGQYMGMDZZ59dG7CKi4uZMmUK48aNo3v37ngFrJcAABU1SURBVGzevJny8vLa/U6cOJHZs2ez1157MWnSJO64446djjJKkiRp9+Mtgo2ksrLyA9fdcsst2yzPmDFju2223LL3fqWlpaSUtmtv27btDiddaNeuHWPGjGHMmDHbrTvkkENqZzEsLXuAG8tP+8D9tHQf9N9ni4svvpiLL764znWDBw/+wMlGevfuzaOPPppJjZIkSWq9HMGSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMuIkF5IkSZJanV5diprlHa0GrEYwbtw4rrzySh588EFOPfVUVq5cyT777MOnPvWpOmcM1K6pqqrivPPO4+mnn+bggw/m97//PQMGDNhpvxkzZjBkyBCKi4tr2/75z39yyCGH8Oabb3LGGWcwb948IoJBgwZxww030LNnT6BmhsEnn3yy9gXPZ555JrfddlvjnKAkSZIabGbZkAb1j2vy6+ctgo3kwAMP5M477wRg6tSplJaWNm9BrcioUaMYMGAAb731Fuecc07tu6zq48Mf/jDV1dW1n0MOOQSADh06cOONN7J8+XKWLVvG4YcfzoUXXrhN31/+8pe1/QxXkiRJqosBq5Ecc8wxPPfcc7z33nvcddddfOELX6hdt2nTJq688kpKS0vp0aMHF110ERs3bgTglVdeYciQIXTt2pWePXsyYsSI2vdbVVZWEhGMHz+ebt260bdvX5566qlmOb/msmrVKqZNm0ZZWRlFRUVccsklLFq0iFmzZjVov0VFRRx22GG0bduWTZs2sWnTJubNm5dR1ZIkSdpdGLAa0QknnMDUqVNZunQp/fr1q22/7rrruO+++5g5cyYVFRXMnTuXCRMmALB+/XouuOACFi9ezPz581m8ePF2LwdetWoVb7zxBmeccQZXXnllk55Tc6uoqKB9+/aUlJQwaNAgqqqq6NevH/Pnz69X/2XLltGjRw/69evH1Vdfvd1LiY844giKioooLy/fbgRr9OjRdOvWjc9+9rOGL0mSJNXJgNWIvvjFL3LxxRdz2mmnbdN+0003MXbsWHr16kXHjh258MILufvuuwE4+OCDOffccykpKaFTp04MHz6c559/fpv+3/nOd2jXrh1Dhw5lwYIFTXY+LcGaNWsoKSlh9erVLFy4kBUrVtCxY8faUb4dOeSQQ5gzZw6vv/46d911FxMnTuSWW27ZZptZs2bx9ttvU15ezqBBg2rbr732WhYvXsyiRYv4+Mc/zuc///naUUdJkiRpCye5qMNx5dNZsnJt3v1X/n0B7614ncenruDNdt359eLurHtiFmtefovSsgdY/PIizjz7XIhcvk2badthL0rLHmDTmpWsfewm9nxzAWvWrGHDhg18/OMf32b/Xbt2BWCPPfZg3bp1eddZiDp06EB1dTW9e/fmjTfeAGD16tWUlJTstG+PHj3o0aMHAEceeSQXXngh999/P9/4xje22a5z58587Wtf4+ijj6ayspK2bdsycODA2vVXX301v/rVr5g3bx5HHHFEhmcnSZKkQmfAqsOSlWupLD9t5xt+gHHj/klFRRv+cM3n4JrPATWjVn9YO4sZ5afR/579ueWWWzj22GO36zty5EgmbdzMS3Pn0rlzZyZMmMDkyZPzrqW1OeCAA1i3bh1VVVXst99+bNiwgZdeeon+/fvDs6/u0r7atPngAdyUElVVVaxYsYJu3bptsy4iiIi86pckSVLr5i2CzeCb3/wmY8aM4bXXXiOlxIIFC5g2bRpQMxrTZo9iSkpKeOWVV5g4cWIzV9uydOrUiZNPPpny8nLWrVvH+PHj6dOnD4cffvg22w0ePJjLL798m7bp06ezePFiAObNm8cNN9zA6aefDsDf//53pk2bxvr166muruaKK67ggAMOoFu3bqxcuZKHHnqI9evXs379eq688kr22Wef2hkIJUmSpC0MWM3g0ksvZdCgQRx//PF06tSJs846i2XLlgEwduxYNryxkM6dO3P22WfXBgD9XxMnTmT27NnstddeTJo0iTvuuGO7EaXKykqWLl26Tdu//vUvjj76aIqLixk6dCgXXHBB7e2BmzdvZvTo0eyzzz7st99+LFmyhPvuuw+AjRs3csUVV7D33nuz77778sQTT3D//ffzoQ99qGlOWJIkSQXDWwQbwbhx47ZrO//88zn//PMBaNeuHWPGjNludkCAww47jJ4jfl7nLYqlpaXbzHo3ePBgqqqqsiu8QPTu3ZtHH310h9tUVlZu13bZZZdx2WWX1bn9pz/9aZ599tk613Xv3p3nnntul+uUJEnS7qfFjGBFRO+IeCQi5kXEnIi4ONfeNSL+GhELcz/3au5aJUmSJKkuLSZgAe8Bl6aUDgE+CVwYEYcCZcDDKaUDgYdzy5IkSZLU4rSYgJVSej2l9K/c99XAPKAXcAZwa26zW4FhzVOhJEmSJO1YiwlYW4uIUuBjwFNAj5TS61ATwoB9PqDPqIh4JiKeWb58eVOVKkmSJEm1WlzAiogS4C7gBymlVfXtl1K6MaU0MKU0sHv37o1XoCRJkiR9gBYVsCLiQ9SEqz+mlO7ONS+NiJ659T2BZc1V366YM2cOJ554Il26dGGfffbhC1/4AgBXX301JSUllJSU0LZtW9q3b1+7vHjxYgYPHszi/zmTNnu0Z89eB9Pj7P9DadkD9f4cVz69mc+88VVVVTF48GCKi4s56qijeOGFF+rV78EHH+Too4+mU6dO9O7dm6uuumqb9RMnTqS0tJSOHTsybNgwVqxYUbtu7dq1fOtb36Jr16506dKF7373u5mekyRJklqHFjNNe9S8yOhmYF5K6bqtVk0Fvg6U537e1wzl7bJhw4bxne98h2nTprF69WomT54MwBVXXMEVV1wB1Eyzft5559VO377FjTdcz5e//GXuuecevvWtb/HAjBkMHDiwXsctLXsg2xNpgUaNGsWAAQN46KGHmDBhAueccw5z5szZab/q6mr++7//m0GDBvHmm29y4oknUlpaynnnncfzzz/Pf/zHf/DEE0/Qt29fhg8fzo9+9CNuuukmAC655BIqKiqYM2cO3bt358UXX2zs05QkSVIBakkjWMcBXwWGRMTzuc9QaoLVZyNiIfDZ3HKL9uabb1JRUcHIkSNp27YtXbp0YdSoUbu0j+LiYr7yla8wfPhwrr/++kaqtPCsWrWKadOmUVZWRlFREZdccgmLFi1i1qxZO+179tlnc9JJJ7HnnnvSq1cvTjnlFP7xj38A8Nhjj3HCCSdw2GGHUVxczPe//33uvrtmEHXt2rXceuut/OIXv6Bnz560a9eOAQMGNOp5SpIkqTC1mICVUvp7SilSSkeklI7MfR5MKb2VUjoxpXRg7ufbzV3rznTt2pU+ffowatQoHn30UdavX5/3vo488kief/75DKsrbBUVFbW3VQ4aNIiqqir69evH/Pnzd3lfTz75JEcccQTANi9w3rK8YsUK3n77bRYsWEBEcO+999KjRw8OPfRQ7rnnnkzOR5IkSa1Li7lFMEsvvrG6QbfK9epS1KDjt2nThocffpixY8cybNgwUkpccskljB07dpf3VVJSwqpV9Z7ro9Vbs2YNJSUlrF69moULF7JixQo6duxIdXU1UFzv/fzqV79i3bp1jBgxAoATTjiB//qv/2L27Nn07duXX//61wC8++67rFq1ig0bNvDyyy+zaNEi/vGPf3D66aezcOFCevbs2QhnKUmSpELVKgPWxk2bqSw/rVlrOOCAA/jjH//I5s2beeSRRxg+fDhHH300Q4cO3aX9VFdX06lTp0aqsvB06NCB6upqevfuzRtvvAHA6tWrKSkpqfc+HnzwQX7605/y2GOPseeeewI1I4U/+clPGDZsGOvWrePCCy9k6tSpdOzYkeLiYjZt2sSll15K+/bt+cxnPsNBBx3Ek08+yZlnntko5ylJkqTC1GJuEWyt2rRpw4knnsjgwYOZO3fuLvd//vnnOfLIIxuhssJ0wAEHsG7dOqqqqgDYsGEDL730Ev37969X/yeeeIJRo0bx5z//mT59+myz7tvf/jYvvfQSS5Ys4YgjjqC0tJTOnTvzkY98hIigZh6WGiml7W4rlCRJkgxYjWDTpk2MGTOmNgT8+9//5vHHH+eoo46q9z7effdd/vjHPzJlyhSnBN9Kp06dOPnkkykvL2fdunWMHz+ePn36cPjhh2+z3eDBg7n88su3aZs1axbDhw9n8uTJdU5SMXv2bDZv3kxFRQWjR4/m29/+NgB77bUXJ5xwAtdddx0bN27k8ccfZ+HChXzyk59svBOVJElSQTJgNYI2bdpQUVHBJz7xCUpKSjjrrLP48Y9/zJAhQ+rV/3vf+x7du3fnF7/4BXfddVe9p2jfXUycOJHZs2ez1157MWnSJO64445tRpcAKisrWbp06TZt48ePZ+nSpZx88sm17x479dRTa9f/53/+J506deJTn/oUn/vc5/jRj35Uu+63v/0t8+fPp0uXLowcOZLbbruND3/4w417opIkSSo4rfIZrOYWEdx+++073W7GjBn1atO2evfuzaOPPrrDbSorK7dr+93vfsfvfve7D+wzderUD1zXt2/fnR5TkiRJcgRLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkiRJkjJiwJIkSZKkjBiwJEmSJCkjBixJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScqIAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKSMEErIg4JSLmR0RFRJQ1dz2SJEmS9H4FEbAioi3wK+BU4FDg3Ig4tHmrkiRJkqRtFUTAAo4BKlJKL6eUNgCTgDOauSZJkiRJ2kaklJq7hp2KiC8Cp6SUzs8tfxX4RErpe1ttMwoYlVscALzQ5IVqd9cNeLO5i9Bux+tOTc1rTs3B607NoX9KqeOudmrXGJU0gqijbZtkmFK6EbgRICKeSSkNbIrCpC287tQcvO7U1Lzm1By87tQcIuKZfPoVyi2CVUDvrZb3A15rplokSZIkqU6FErD+CRwYEX0jYg/gS8DUZq5JkiRJkrZRELcIppTei4jvAX8B2gK/TSnN2UGXG5umMmkbXndqDl53ampec2oOXndqDnlddwUxyYUkSZIkFYJCuUVQkiRJklo8A5YkSZIkZaSgA1ZEnBIR8yOiIiLK6li/Z0TckVv/VESUNn2Vam3qcd39MCLmRsSsiHg4IvZvjjrVeuzsmttquy9GRIoIpzJWg9XnuouIs3O/7+ZExO1NXaNan3r8GdsnIh6JiOdyf84ObY461XpExG8jYllE1PkO3agxIXdNzoqIo3a2z4INWBHRFvgVcCpwKHBuRBz6vs1GAitSSgcA44FrmrZKtTb1vO6eAwamlI4A7gR+0rRVqjWp5zVHRHQELgKeatoK1RrV57qLiAOB0cBxKaXDgB80eaFqVer5++7HwOSU0seomVX6+qatUq3QLcApO1h/KnBg7jMKuGFnOyzYgAUcA1SklF5OKW0AJgFnvG+bM4Bbc9/vBE6MiLpeWizV106vu5TSIymld3OLT1Lz3jYpX/X5XQfw/1IT5tc1ZXFqtepz3V0A/CqltAIgpbSsiWtU61Of6y4BnXLfO+N7UdVAKaXHgLd3sMkZwO9TjSeBLhHRc0f7LOSA1Qt4davlqlxbnduklN4D3gH2bpLq1FrV57rb2kjgoUatSK3dTq+5iPgY0Dul9OemLEytWn1+1x0EHBQRMyPiyYjY0b8AS/VRn+tuHHBeRFQBDwLfb5rStBvb1b/7FcZ7sD5AXSNR759zvj7bSLui3tdURJwHDAROaNSK1Nrt8JqLiDbU3AI9oqkK0m6hPr/r2lFzy8xgakbqH4+IASmllY1cm1qv+lx35wK3pJT+JyKOBW7LXXebG7887aZ2OU8U8ghWFdB7q+X92H6YuHabiGhHzVDyjoYApZ2pz3VHRJwE/Cfw+ZTS+iaqTa3Tzq65jsAAYEZEVAKfBKY60YUaqL5/xt6XUtqYUnoFmE9N4JLyVZ/rbiQwGSCl9A+gPdCtSarT7qpef/fbWiEHrH8CB0ZE34jYg5oHHae+b5upwNdz378ITE++WVkNs9PrLne71kRqwpXPJKihdnjNpZTeSSl1SymVppRKqXnu7/MppWeap1y1EvX5M/Ze4DMAEdGNmlsGX27SKtXa1Oe6WwycCBARh1ATsJY3aZXa3UwFvpabTfCTwDsppdd31KFgbxFMKb0XEd8D/gK0BX6bUpoTEf8HeCalNBW4mZqh4wpqRq6+1HwVqzWo53X3U6AEmJKbU2VxSunzzVa0Clo9rzkpU/W87v4CnBwRc4FNwI9SSm81X9UqdPW87i4FfhMRl1Bzm9YI//FcDRERf6LmVuduuWf7xgIfAkgp/ZqaZ/2GAhXAu8A3drpPr0lJkiRJykYh3yIoSZIkSS2KAUuSJEmSMmLAkiRJkqSMGLAkSZIkKSMGLEmSJEnKiAFLkrTbi4jSiHihueuQJBU+A5YkSZIkZcSAJUkqKLnRphcj4taImBURd0ZE8fu2uSMihm61fEtEfCHX9/GI+Ffu86k69j8iIn651fKfI2Jw7vvJEfGPXN8pEVHSiKcqSSpABixJUiHqD9yYUjoCWAV8933rJwHnAETEHsCJwIPAMuCzKaWjcusn1PeAEdEN+DFwUq7/M8APG3gekqRWxoAlSSpEr6aUZua+/wE4/n3rHwKGRMSewKnAYymltcCHgN9ExGxgCnDoLhzzk7ntZ0bE88DXgf0bcA6SpFaoXXMXIElSHtL7ljvnQg/AmJTS1IiYAfw/1IxU/Sm37hJgKfBRav6RcV0d+36Pbf8Bsn3uZwB/TSmd2/DyJUmtlSNYkqRC1Ccijs19Pxf4c0rpyNxnaq59EvANYBDwl1xbZ+D1lNJm4KtA2zr2XQkcGRFtIqI3cEyu/UnguIg4ACAiiiPioKxPTJJU2AxYkqRCNA/4ekTMAroCN9SxzTTg08DfUkobcm3X5/o9CRwErKmj30zgFWA2cC3wL4CU0nJgBPCn3HGfBA7O6oQkSa1DpPT+uywkSWq5IqKUmhGrAc1ciiRJ23EES5IkSZIy4giWJEmSJGXEESxJkiRJyogBS5IkSZIyYsCSJEmSpIwYsCRJkiQpIwYsSZIkScrI/w+KQ4Nt/OQSwQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x1008 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "N_bins = 50\n",
    "\n",
    "if (N_exp > 1):\n",
    "    fig2, ax2 = plt.subplots(nrows=3, figsize=(12, 14))\n",
    "    \n",
    "    ax2[0].hist(all_p_mean, N_bins, (0, 1), histtype='step')\n",
    "    ax2[0].set(title='Histogram, probability mu', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n",
    "    ax_text(all_p_mean, ax2[0], 0.04, 0.25)\n",
    "    \n",
    "\n",
    "    ax2[1].hist(all_p_chi2, N_bins, (0, 1), histtype='step')\n",
    "    ax2[1].set(title='Histogram, probability chi2', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n",
    "    ax_text(all_p_chi2, ax2[1], 0.04, 0.25)\n",
    "    \n",
    "    ax2[2].hist(all_p_ks, N_bins, (0, 1), histtype='step')\n",
    "    ax2[2].set(title='Histogram, probability Kolmogorov', xlabel='p-value', ylabel='Frequency / 0.02', xlim=(0, 1))\n",
    "    ax_text(all_p_ks, ax2[2], 0.04, 0.25)\n",
    "\n",
    "    fig2.tight_layout()\n",
    "\n",
    "\n",
    "    if save_plots:\n",
    "        fig2.savefig('PvalueDists.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Questions:\n",
    "\n",
    "1. First run the program with one experiment (N_exp = 1) to display the two distributions A and B, when:\n",
    "    - They are the same.\n",
    "    - The mean of A is increased (to e.g. 0.1).\n",
    "    - The width of A is enlarged (to e.g. 1.2).\n",
    "    - The normalisation of A is increased.\n",
    "\n",
    "Get a feel for how much you need to change the distribution, before you can _by eye_ see that they are not the same. I.e. can you see any difference, if `mean_A` $= 0.1$? Or how about $0.2$? How do you quantify this and when do you start to doubt? And how about `width_A` $= 1.1$? Or $1.2$? Again, can you see it by eye? Finally, try to put $1050$ events into B. Is that visible? How about $1100$?<br>\n",
    "When you have an idea about when you can see effects by eye, now check if you see an impact in the p-values?\n",
    "\n",
    "2. Could you for the test of the means have calculated how much of a change in the mean is needed for a difference to be statistically significant? Do so, and see if it somewhat matches you guess/estimate from above!\n",
    "\n",
    "\n",
    "3. Now run the tests 1000 times, where A and B are unit Gaussians and thus identical. How should the distributions of the test probabilities come out? And is this the case, approximately? If not, think of reasons for this, and what could be a remedy. HINT: Large statistics is always easier!\n",
    "\n",
    "\n",
    "4. Repeat the changes in question 1), and see which tests \"reacts\" most to these modifications. How much of a change in the mean is required for 95% of the tests (of each kind) to give a probability below 5%? How much is required for the width? And the norm?\n",
    "\n",
    "\n",
    "5. Possibly try to test different distributions than the Gaussian one (e.g. exponential, uniform, etc.), and see how the tests performs.\n",
    "\n",
    "\n",
    "NOTE: The Kolmogorov-Smirnov test has the great advantage that it can handle ANY distribution (even the Cauchy distribution - remind yourself of that one!). The reason is, that it doesn't care about any PDF, nor how far out an outlier is. It is just a matter of the difference in integrals between the two functions.\n",
    "\n",
    "\n",
    "## Advanced:\n",
    "\n",
    "6. Obviously, the test of the means is not sensitive the a change in the width. Make such a test yourself by calculating the widths and the uncertainty on the widths (or perhaps try the F-test!). Note that in a (unit) Gaussian the uncertainty on the width is of the same order as that of the means!\n",
    "\n",
    "\n",
    "## Very advanced:\n",
    "7. Implement in python the following tests:\n",
    "     - Lilliefors test\n",
    "     - Shapiro-Wilk test\n",
    "     - Anderson-Darling test\n",
    "     - Cramer-von-Mises test\n",
    "     - Jarque-Bera test\n",
    "     - Kuiper's test\n",
    "     - Mann-Whitney-Wilcoxon test\n",
    "     - Siegel-Tukey test\n",
    "     \n",
    "and quantify under various conditions and datasets the power of each and the correlation among them. Write it up, and send it to a statistics journal. :-)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}