CP5_PDFs_fits_z-scores.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c2decc5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Making the necessary imports\n",
    "import os\n",
    "import numpy as np\n",
    "from numpy import matrix\n",
    "import scipy.optimize as theMinimiser\n",
    "from iminuit import Minuit\n",
    "import math\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "import scipy.integrate as integrate\n",
    "from scipy.stats import norm\n",
    "from scipy.special import erfinv\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dd0ad2fc",
   "metadata": {},
   "source": [
    "## Problem 1 (3 points):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b3460c51",
   "metadata": {},
   "outputs": [],
   "source": [
    "# background parameters\n",
    "XMIN = 0.\n",
    "XMAX = 20.\n",
    "intercept = 20.\n",
    "slope = -1.\n",
    "# signal parameters\n",
    "mean = 10.\n",
    "sigma = 0.5\n",
    "# chose number of bins that is appropriate for the size of the statistics sample\n",
    "NBINS = 100"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19111187",
   "metadata": {},
   "source": [
    "### Creating the background and signal classes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ea9bae9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Linear:\n",
    "    def __init__(self, intercept, slope, XMIN, XMAX):\n",
    "        self.mass = []\n",
    "\n",
    "        self.intercept = intercept\n",
    "        self.slope = slope\n",
    "        self.lolimit = XMIN\n",
    "        self.hilimit = XMAX\n",
    "        \n",
    "    def maxval(self):\n",
    "        return self.intercept \n",
    "        \n",
    "    # Evaluate method (un - normalised )\n",
    "    def evaluate(self, x):\n",
    "        linear_dist = self.intercept + (self.slope * x)\n",
    "        return linear_dist\n",
    "    \n",
    "\n",
    "    def next(self):\n",
    "        doLoop = True\n",
    "        while (doLoop):\n",
    "        # start with uniform random number in [ lolimit , hilimit )\n",
    "            x = np.random.uniform(self.lolimit, self.hilimit)\n",
    "            y1 = self.evaluate(x)\n",
    "            y2 = np.random.uniform(0, self.maxval())\n",
    "            #print(self.maxval())\n",
    "            if (y2 < y1):\n",
    "                filtered_x = x\n",
    "                self.mass.append(filtered_x)\n",
    "                return filtered_x   \n",
    "            \n",
    "            "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "31ba025f",
   "metadata": {},
   "outputs": [],
   "source": [
    "class Gaussian:\n",
    "    def __init__(self, mean, sigma):\n",
    "        self.mass = []\n",
    "        \n",
    "        self.mean = mean\n",
    "        self.sigma = sigma\n",
    "        \n",
    "    def next(self):\n",
    "        gaussian_dist = np.random.normal(self.mean,self.sigma)\n",
    "        self.mass.append(gaussian_dist)\n",
    "        return gaussian_dist\n",
    "    \n",
    "    def integrate(self,lower_lim,upper_lim) :\n",
    "        gaussian = lambda x: (1/(self.sigma*np.sqrt(2*math.pi)))*np.exp((-1/2)*((x-self.mean)/self.sigma)**2)\n",
    "        integral, _ = integrate.quad(gaussian,lower_lim,upper_lim)\n",
    "        return integral"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a80025f",
   "metadata": {},
   "source": [
    "### Creating class to sum the background and signal distributions,  weighted by the fractions of background and signal events"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bf953c36",
   "metadata": {},
   "outputs": [],
   "source": [
    "class SignalWithBackground:\n",
    "    def __init__(self, mean, sigma, sig_fraction, intercept, slope, XMIN, XMAX):\n",
    "        self.mean = mean\n",
    "        self.sigma = sigma\n",
    "        self.sig_fraction = sig_fraction\n",
    "        self.intercept = intercept\n",
    "        self.slope = slope\n",
    "        self.XMIN = XMIN\n",
    "        self.XMAX = XMAX\n",
    "        self.signal = Gaussian(mean,sigma)\n",
    "        self.background = Linear(intercept, slope, XMIN, XMAX)\n",
    "        self.mass_sig = []\n",
    "        self.mass_bgd = []\n",
    "        self.mass = [] \n",
    "    \n",
    "    # Draw random number from distribution\n",
    "    def next(self):\n",
    "        q = np.random.uniform()\n",
    "        if(q < self.sig_fraction):\n",
    "            # if here , we will draw x from signal distribution\n",
    "            filtered_x = self.signal.next()\n",
    "            self.mass_sig.append(filtered_x)\n",
    "        else:\n",
    "            # if here , we will draw x from background distribuion\n",
    "            filtered_x = self.background.next()\n",
    "            self.mass_bgd.append(filtered_x)\n",
    "        self.mass.append(filtered_x)\n",
    "        return filtered_x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af5cb92d",
   "metadata": {},
   "source": [
    "### Creating function to plot the distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "e083168f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Main code to generate and plot a single experiment\n",
    "def singleToy(nevents_sig, nevents_bgd):\n",
    "    sig_fraction = nevents_sig/(nevents_bgd+nevents_sig)\n",
    "    \n",
    "    # Create the pdf\n",
    "    pdf = SignalWithBackground(mean, sigma, sig_fraction, intercept, slope, XMIN, XMAX)\n",
    "    for i in range (nevents_sig + nevents_bgd): pdf.next()\n",
    "        \n",
    "    # retrieve the mass values for signal, background and their sum\n",
    "    data = pdf.mass\n",
    "    sig_data = pdf.mass_sig\n",
    "    bgd_data = pdf.mass_bgd\n",
    "    \n",
    "    # plot things on same page\n",
    "    myRange = (XMIN, XMAX)\n",
    "    fig, axs = plt.subplots(3, 1, sharex =\"col\" ,figsize=(5,10))\n",
    "    axs[0].set_title(\"Signal distribution (\" + str(len(sig_data)) + \" entries)\")\n",
    "    axs[1].set_title(\"Background distribution (\" + str(len(bgd_data)) + \" entries)\")\n",
    "    axs[2].set_title(\"Total distribution (\" + str(len(data)) + \" entries)\")\n",
    "    axs[2].set_xlabel(\"X\")\n",
    "    axs[0].hist(sig_data, bins=NBINS, range=myRange)\n",
    "    axs[1].hist(bgd_data, bins=NBINS)\n",
    "    axs[2].hist(data, bins=NBINS)\n",
    "    fig.tight_layout()\n",
    "    plt.savefig(\"Signal, background and combined distributions.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "99b81616",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAPdCAYAAABFjSsGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAACTDklEQVR4nOzdeVxUVeM/8M8AMuyDoGzJJi64m6hImoKSiKYZqLlUqChZoAmVSlku2YNZj/poLi2GZfI1NdMnTczdTDTFTM0kMVxKQdNgBHXYzu+PftzH6wyrM3KBz/v1mlfdc88999xh5MO599w7KiGEABERESmSWW13gIiIiMrHoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoKYa8/HxwdixY2u7GwCACxcuQKVSYfXq1TXafvXq1VCpVLhw4YJUFhwcjODgYKP0rzIqlQqzZ8+WlmfPng2VSoW//vrroey/tn+WCxYsgL+/P0pLS2utD3WJKT+bZ86cgYWFBU6fPm2S9qn6GNSk59SpUxg2bBi8vb1hZWWFRx55BE888QSWLl1a211TvEOHDmH27NnIzc2t7a7oUWrftFot3n33XUyfPh1mZv/7lRQfH48uXbrAyckJNjY2aNOmDWbPno38/HzZ9kePHkVcXBzatWsHW1tbeHl5YcSIEfjtt9/09vXxxx+jT58+cHV1hVqthq+vL8aNGyf7A+1hunLlCmbPno0TJ07Uyv4Nadu2LQYNGoS33nqrtrtC/59FbXeAlOXQoUMICQmBl5cXJk6cCDc3N1y+fBmHDx/Gf/7zH0yePFmqm5GRIfvFWt9899131d7m0KFDmDNnDsaOHQtHR8cqb3fnzh1YWJj2n2NFfavNn+Wnn36K4uJijBo1SlZ+9OhRPP744xg3bhysrKzw008/Yf78+di1axcOHDgg9ffdd9/FDz/8gOHDh6Njx47Izs7GBx98gC5duuDw4cNo37691OZPP/0EX19fDBkyBI0bN0ZWVhY+/vhjbN26FT///DM8PDwe6rFfuXIFc+bMgY+PDzp37lzl7Wry2ayOSZMmYeDAgTh//jz8/PxMui+qHIOaZN555x1oNBocPXpU75f5tWvXZMtqtfoh9uzhs7S0NGn7paWlKCwshJWVFaysrEy6r8rU5s8yOTkZQ4YM0XsPDh48qFfXz88Pr776Kn788Uf06NEDAJCQkICUlBTZz+uZZ55Bhw4dMH/+fHzxxRdS+fLly/XaHDp0KLp27YrPP/8cM2bMMNZhmcTt27dhY2Nj8s9maGgoGjdujM8++wxz58416b6ocvV3OEQ1cv78ebRr187gaNDFxUW2bOi65smTJ9GnTx9YW1ujWbNmmDdvHpKTk/Wu//r4+ODJJ5/EwYMH0b17d1hZWaF58+b4/PPPZe3dvHkTr776Kjp06AA7Ozs4ODggPDwcP//8c42P8ZdffkHfvn1lfTR0bdTQdcClS5eiXbt2sLGxQePGjdG1a1ekpKQA+Oe68muvvQYA8PX1hUqlkh23SqVCXFwc1q5di3bt2kGtViM1NVVad+816jJ//fUXRowYAQcHBzg7O+Pll1/G3bt3pfUVXZu/t83K+mboZ/n7779j+PDh0qnnHj16YNu2bbI6+/btg0qlwvr16/HOO++gWbNmsLKyQr9+/ZCZmanXp/tlZWXh5MmTCA0NrbRuWT8ByE7fP/bYY3rB1bJlS7Rr1w6//vprjdqsyBdffIGAgABYW1vDyckJI0eOxOXLl2V1goOD0b59e5w5cwYhISGwsbHBI488ggULFkh19u3bh27dugEAxo0bJ/1Myn6WZW2kp6ejd+/esLGxweuvvy6tu/+zqdPpMGvWLLRo0QJqtRqenp6YNm0adDqdrN7OnTvRq1cvODo6ws7ODq1bt5baLdOoUSMEBwdjy5YtVXpPyLQ4oiYZb29vpKWl4fTp07JThlXx559/IiQkBCqVComJibC1tcUnn3xS7mgtMzMTw4YNQ3R0NKKiovDpp59i7NixCAgIQLt27QD8ExabN2/G8OHD4evri5ycHHz44Yfo06cPzpw5U+1TldnZ2QgJCUFxcTFmzJgBW1tbfPTRR7C2tq50248//hhTpkzBsGHDpMA8efIkjhw5gtGjRyMiIgK//fYb/u///g+LFi1CkyZNAABNmzaV2tizZw/Wr1+PuLg4NGnSRAqJ8owYMQI+Pj5ISkrC4cOHsWTJEvz99996f9BUpip9u1dOTg4ee+wx3L59G1OmTIGzszM+++wzDBkyBBs3bsTTTz8tqz9//nyYmZnh1VdfRV5eHhYsWIAxY8bgyJEjFfbr0KFDAIAuXboYXF9cXIzc3FwUFhbi9OnTmDlzJuzt7dG9e/cK2xVCICcnR/oc3e/GjRsoKSnBpUuXpBFjv379KmwT+OeM05tvvokRI0ZgwoQJuH79OpYuXYrevXvjp59+kv2B+/fff2PAgAGIiIjAiBEjsHHjRkyfPh0dOnRAeHg42rRpg7lz5+Ktt95CTEwMHn/8cQD//OFxbz/Dw8MxcuRIPPvss3B1dTXYr9LSUgwZMgQHDx5ETEwM2rRpg1OnTmHRokX47bffsHnzZgD//JH65JNPomPHjpg7dy7UajUyMzPxww8/6LUZEBCALVu2QKvVwsHBodL3hkxIEN3ju+++E+bm5sLc3FwEBQWJadOmiR07dojCwkK9ut7e3iIqKkpanjx5slCpVOKnn36Sym7cuCGcnJwEAJGVlSXbFoA4cOCAVHbt2jWhVqvFK6+8IpXdvXtXlJSUyPablZUl1Gq1mDt3rqwMgEhOTq7w+KZOnSoAiCNHjsj2q9Fo9PrYp08f0adPH2n5qaeeEu3atauw/ffee0+vnTIAhJmZmfjll18Mrps1a5a0PGvWLAFADBkyRFbvpZdeEgDEzz//LISo+Ljvb7Oivt3/syx7n77//nup7NatW8LX11f4+PhIP5O9e/cKAKJNmzZCp9NJdf/zn/8IAOLUqVN6+7rXzJkzBQBx69Ytg+vT0tIEAOnVunVrsXfv3grbFEKINWvWCABi1apVBter1WqpTWdnZ7FkyZJK27xw4YIwNzcX77zzjqz81KlTwsLCQlbep08fAUB8/vnnUplOpxNubm4iMjJSKjt69Gi5P7+yNlauXGlw3b2fzTVr1ggzMzPZz0sIIVauXCkAiB9++EEIIcSiRYsEAHH9+vVKjzclJUXv3wrVDp76JpknnngCaWlpGDJkCH7++WcsWLAAYWFheOSRR/Df//63wm1TU1MRFBQkmxTj5OSEMWPGGKzftm1baRQB/DO6a926NX7//XepTK1WS5OGSkpKcOPGDel03fHjx6t9fN9++y169OghG5E1bdq03D7ey9HREX/88QeOHj1a7f2W6dOnD9q2bVvl+rGxsbLlssl83377bY37UBXffvstunfvjl69eklldnZ2iImJwYULF3DmzBlZ/XHjxslOP5f9XO/9WRpy48YNWFhYwM7OzuD6tm3bYufOndi8eTOmTZsGW1tbvVnf9zt79ixiY2MRFBSEqKgog3W2b9+Ob7/9Fv/+97/h5eWFgoKCCtsEgE2bNqG0tBQjRozAX3/9Jb3c3NzQsmVL7N27V1bfzs4Ozz77rLRsaWmJ7t27V/qe3EutVmPcuHGV1tuwYQPatGkDf39/Wd/69u0LAFLfykb8W7ZsqfRWuMaNGwPAQ7tFkMrHoCY93bp1w6ZNm/D333/jxx9/RGJiIm7duoVhw4bp/YK+18WLF9GiRQu9ckNlAODl5aVX1rhxY/z999/ScmlpKRYtWoSWLVtCrVajSZMmaNq0KU6ePIm8vLxqH9vFixfRsmVLvfLWrVtXuu306dNhZ2eH7t27o2XLloiNjTV4yrAivr6+1ap/f1/9/PxgZmZm8tuJLl68aPA9adOmjbT+Xvf/LMt+yd/7s6wJBwcHhIaG4qmnnsK7776LV155BU899VS5cxSys7MxaNAgaDQabNy4Eebm5gbrhYSEIDw8HAkJCdiwYQPmzJmDDz74oMK+nDt3DkIItGzZEk2bNpW9fv31V73Jls2aNYNKpZKV3f/5rswjjzxSpYlj586dwy+//KLXr1atWgH430TQZ555Bj179sSECRPg6uqKkSNHYv369QZDWwgBAHrHQA8fr1FTuSwtLdGtWzd069YNrVq1wrhx47BhwwbMmjXLKO2X90u07BcEAPzrX//Cm2++ifHjx+Ptt9+Gk5MTzMzMMHXq1If+cIw2bdogIyMDW7duRWpqKr766issX74cb731FubMmVOlNqpyLbwi9//SLO+XaElJyQPtp7qq8rM0xNnZGcXFxbh16xbs7e0r3U9ERASee+45rFu3Dp06dZKty8vLQ3h4OHJzc/H9999Xef6Cn58fHn30UaxduxZxcXHl1istLYVKpcL27dsNHu/9ZwVq+p7cq6qfl9LSUnTo0AELFy40uN7T01Nq78CBA9i7dy+2bduG1NRUfPnll+jbty++++47WZ/L/qAom89AtYdBTVXStWtXAMDVq1fLrePt7W1wpm9VZv+WZ+PGjQgJCcGqVatk5bm5uTX6BeLt7Y1z587plWdkZFRpe1tbWzzzzDN45plnUFhYiIiICLzzzjtITEyElZWV0Ucf586dk43CMzMzUVpaKk1CKxu53j9j+f4RL1C9kZG3t7fB9+Ts2bPSemPw9/cH8M/s744dO1ZaX6fTobS0VO9syt27dzF48GD89ttv2LVrV7UuLwD/3Md+/+zo+/n5+UEIAV9fX2mk+qCM9Xnx8/PDzz//jH79+lXappmZGfr164d+/fph4cKF+Ne//oU33ngDe/fulc2+z8rKgpmZmdGOlWqOp75JZu/evQb/4i+7JlrRKeKwsDCkpaXJnrJ08+ZNrF27tsb9MTc31+vPhg0b8Oeff9aovYEDB+Lw4cP48ccfpbLr169XqY83btyQLVtaWqJt27YQQqCoqAjAP0EOVP1Wn8osW7ZMtlz2dLjw8HAA/5wabtKkCQ4cOCCrZ+h+4er0beDAgfjxxx+RlpYmlRUUFOCjjz6Cj49PtYOwPEFBQQCAY8eOycpzc3Ol9/Ren3zyCYD//eEI/HP24JlnnkFaWho2bNggtXm/4uJig6edf/zxR5w6dUrWpiEREREwNzfHnDlz9D6TQgi9z0dVGOvzMmLECPz555/4+OOP9dbduXNHugZ/8+ZNvfVlc0ru/0MlPT0d7dq1g0ajeaC+0YPjiJpkJk+ejNu3b+Ppp5+Gv78/CgsLcejQIXz55Zfw8fGpcGLLtGnT8MUXX+CJJ57A5MmTpduzvLy8cPPmzRqNHp588knMnTsX48aNw2OPPYZTp05h7dq1aN68eY2Ob9q0aVizZg0GDBiAl19+Wbo9y9vbGydPnqxw2/79+8PNzQ09e/aEq6srfv31V3zwwQcYNGiQdNo2ICAAAPDGG29g5MiRaNSoEQYPHiz9Qq6urKwsDBkyBAMGDEBaWhq++OILjB49Wnbad8KECZg/fz4mTJiArl274sCBAwYfn1mdvs2YMQP/93//h/DwcEyZMgVOTk747LPPkJWVha+++spoTzFr3rw52rdvj127dmH8+PFS+b59+6Rb4Vq2bInCwkJ8//332LRpE7p27SqbpPXKK6/gv//9LwYPHoybN2/KHnACQKqbn58PT09PPPPMM9LjRk+dOoXk5GRoNBq8+eabFfbVz88P8+bNQ2JiIi5cuIChQ4fC3t4eWVlZ+PrrrxETE4NXX321Wsfv5+cHR0dHrFy5Evb29rC1tUVgYGC15zI899xzWL9+PSZNmoS9e/eiZ8+eKCkpwdmzZ7F+/Xrs2LEDXbt2xdy5c3HgwAEMGjQI3t7euHbtGpYvX45mzZrJJg4WFRVh//79eOmll6rVDzKRWpptTgq1fft2MX78eOHv7y/s7OyEpaWlaNGihZg8ebLIycmR1b3/lh4hhPjpp5/E448/LtRqtWjWrJlISkoSS5YsEQBEdna2bNtBgwbp7f/+207u3r0rXnnlFeHu7i6sra1Fz549RVpaml69qt6eJYQQJ0+eFH369BFWVlbikUceEW+//bZYtWpVpbdnffjhh6J3797C2dlZqNVq4efnJ1577TWRl5cna//tt98WjzzyiDAzM5O1CUDExsYa7BPKuT3rzJkzYtiwYcLe3l40btxYxMXFiTt37si2vX37toiOjhYajUbY29uLESNGiGvXrum1WVHfDP0sz58/L4YNGyYcHR2FlZWV6N69u9i6dausTtntWRs2bJCVV+fnsXDhQmFnZydu374tlWVmZornn39eNG/eXFhbWwsrKyvRrl07MWvWLJGfny/bvuw2pvJeZXQ6nXj55ZdFx44dhYODg2jUqJHw9vYW0dHRBm9ZK89XX30levXqJWxtbYWtra3w9/cXsbGxIiMjQ9YnQ7fyRUVFCW9vb1nZli1bRNu2bYWFhYXsPSuvjbJ19342hRCisLBQvPvuu6Jdu3ZCrVaLxo0bi4CAADFnzhzpM7p7927x1FNPCQ8PD2FpaSk8PDzEqFGjxG+//SZra/v27QKAOHfuXJXfFzIdlRDVmNlAVANTp07Fhx9+iPz8/HIn2FDDlZeXh+bNm2PBggWIjo6u7e4Q/nmsqkqlwtdff13bXSEADGoyqjt37shmqt64cQOtWrVCly5dsHPnzlrsGSnZu+++i+TkZJw5c6Zef9FLXfDrr7+iQ4cOOHHiRLWfTkimwaAmo+rcuTOCg4PRpk0b5OTkYNWqVbhy5Qp2796N3r1713b3iIjqHE4mI6MaOHAgNm7ciI8++ggqlQpdunTBqlWrGNJERDXEETUREZGC8WIQERGRgtXJU9+lpaW4cuUK7O3t+RxaIiKqc4QQuHXrFjw8PCqdQFkng/rKlSvSs2uJiIjqqsuXL6NZs2YV1qmTQV32FKjLly/zC82JiKjO0Wq18PT0rNKX0dTJoC473e3g4MCgJiKiOqsql285mYyIiEjBGNREREQKxqAmIiJSMAY1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKVicfIUpED5/PjG2y5QvzB9VST4gaFo6oiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKZvSgnj17NlQqlezl7+8vrb979y5iY2Ph7OwMOzs7REZGIicnx9jdICIiqhdMMqJu164drl69Kr0OHjworYuPj8c333yDDRs2YP/+/bhy5QoiIiJM0Q0iIqI6zyQPPLGwsICbm5teeV5eHlatWoWUlBT07dsXAJCcnIw2bdrg8OHD6NGjh8H2dDoddDqdtKzVak3RbSIiIsUxyYj63Llz8PDwQPPmzTFmzBhcunQJAJCeno6ioiKEhoZKdf39/eHl5YW0tLRy20tKSoJGo5Fenp6epug2ERGR4hg9qAMDA7F69WqkpqZixYoVyMrKwuOPP45bt24hOzsblpaWcHR0lG3j6uqK7OzscttMTExEXl6e9Lp8+bKxu01ERKRIRj/1HR4eLv1/x44dERgYCG9vb6xfvx7W1tY1alOtVkOtVhuri0RERHWGyW/PcnR0RKtWrZCZmQk3NzcUFhYiNzdXVicnJ8fgNW0iIqKGzuRBnZ+fj/Pnz8Pd3R0BAQFo1KgRdu/eLa3PyMjApUuXEBQUZOquEBER1TlGP/X96quvYvDgwfD29saVK1cwa9YsmJubY9SoUdBoNIiOjkZCQgKcnJzg4OCAyZMnIygoqNwZ30RERA2Z0YP6jz/+wKhRo3Djxg00bdoUvXr1wuHDh9G0aVMAwKJFi2BmZobIyEjodDqEhYVh+fLlxu4GERFRvaASQoja7kR1abVaaDQa5OXlwcHBoba7Q9Qg+MzYJlu+MH9QLfWEqO6rTo7xWd9EREQKxqAmIiJSMAY1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKxqAmIiJSMAY1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKZlHbHSAiZfKZsa22u0BE4IiaiIhI0RjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRg/JpLogbq/q+xvDB/UC31hIgqwhE1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESmYyYN6/vz5UKlUmDp1qlR29+5dxMbGwtnZGXZ2doiMjEROTo6pu0JERFTnmDSojx49ig8//BAdO3aUlcfHx+Obb77Bhg0bsH//fly5cgURERGm7AoREVGdZLKgzs/Px5gxY/Dxxx+jcePGUnleXh5WrVqFhQsXom/fvggICEBycjIOHTqEw4cPm6o7REREdZLJgjo2NhaDBg1CaGiorDw9PR1FRUWycn9/f3h5eSEtLc1gWzqdDlqtVvYiIiJqCEzyZLJ169bh+PHjOHr0qN667OxsWFpawtHRUVbu6uqK7Oxsg+0lJSVhzpw5pugqERGRohl9RH358mW8/PLLWLt2LaysrIzSZmJiIvLy8qTX5cuXjdIuERGR0hk9qNPT03Ht2jV06dIFFhYWsLCwwP79+7FkyRJYWFjA1dUVhYWFyM3NlW2Xk5MDNzc3g22q1Wo4ODjIXkRERA2B0U999+vXD6dOnZKVjRs3Dv7+/pg+fTo8PT3RqFEj7N69G5GRkQCAjIwMXLp0CUFBQcbuDhERUZ1m9KC2t7dH+/btZWW2trZwdnaWyqOjo5GQkAAnJyc4ODhg8uTJCAoKQo8ePYzdHSIiojqtVr7mctGiRTAzM0NkZCR0Oh3CwsKwfPny2ugKERGRoj2UoN63b59s2crKCsuWLcOyZcsexu6JiIjqLD7rm4iISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUrFa+5pKIlMdnxrba7gIRGcARNRERkYIxqImIiBSMQU1ERKRgDGoiIiIF42QyogaCk8WI6iaOqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGCeTEVGN3D857cL8QbXUE6L6jSNqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUzOhBvWLFCnTs2BEODg5wcHBAUFAQtm/fLq2/e/cuYmNj4ezsDDs7O0RGRiInJ8fY3SAiIqoXjB7UzZo1w/z585Geno5jx46hb9++eOqpp/DLL78AAOLj4/HNN99gw4YN2L9/P65cuYKIiAhjd4OIiKheUAkhhKl34uTkhPfeew/Dhg1D06ZNkZKSgmHDhgEAzp49izZt2iAtLQ09evQwuL1Op4NOp5OWtVotPD09kZeXBwcHB1N3n6he8JmxzaTtX5g/yKTtE9UnWq0WGo2mSjlm0mvUJSUlWLduHQoKChAUFIT09HQUFRUhNDRUquPv7w8vLy+kpaWV205SUhI0Go308vT0NGW3iYiIFMMkQX3q1CnY2dlBrVZj0qRJ+Prrr9G2bVtkZ2fD0tISjo6Osvqurq7Izs4ut73ExETk5eVJr8uXL5ui20RERIpjYYpGW7dujRMnTiAvLw8bN25EVFQU9u/fX+P21Go11Gq1EXtIRERUN5gkqC0tLdGiRQsAQEBAAI4ePYr//Oc/eOaZZ1BYWIjc3FzZqDonJwdubm6m6AoREVGd9lDuoy4tLYVOp0NAQAAaNWqE3bt3S+syMjJw6dIlBAUFPYyuEBER1SlGH1EnJiYiPDwcXl5euHXrFlJSUrBv3z7s2LEDGo0G0dHRSEhIgJOTExwcHDB58mQEBQWVO+ObiIioITN6UF+7dg3PP/88rl69Co1Gg44dO2LHjh144oknAACLFi2CmZkZIiMjodPpEBYWhuXLlxu7G0RERPXCQ7mP2tiqc/8ZEf2D91ETKYdi7qMmIiKiB8OgJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYJZ1HYHiKh+8JmxTbZ8Yf6gWuoJUf3CETUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGC8clkRPXU/U8KI6K6iSNqIiIiBTN6UCclJaFbt26wt7eHi4sLhg4dioyMDFmdu3fvIjY2Fs7OzrCzs0NkZCRycnKM3RUiIqI6z+hBvX//fsTGxuLw4cPYuXMnioqK0L9/fxQUFEh14uPj8c0332DDhg3Yv38/rly5goiICGN3hYiIqM4z+jXq1NRU2fLq1avh4uKC9PR09O7dG3l5eVi1ahVSUlLQt29fAEBycjLatGmDw4cPo0ePHnpt6nQ66HQ6aVmr1Rq720RERIpk8mvUeXl5AAAnJycAQHp6OoqKihAaGirV8ff3h5eXF9LS0gy2kZSUBI1GI708PT1N3W0iIiJFMGlQl5aWYurUqejZsyfat28PAMjOzoalpSUcHR1ldV1dXZGdnW2wncTEROTl5Umvy5cvm7LbREREimHS27NiY2Nx+vRpHDx48IHaUavVUKvVRuoVERFR3WGyEXVcXBy2bt2KvXv3olmzZlK5m5sbCgsLkZubK6ufk5MDNzc3U3WHiIioTjJ6UAshEBcXh6+//hp79uyBr6+vbH1AQAAaNWqE3bt3S2UZGRm4dOkSgoKCjN0dIiKiOs3op75jY2ORkpKCLVu2wN7eXrrurNFoYG1tDY1Gg+joaCQkJMDJyQkODg6YPHkygoKCDM74JiIiasiMHtQrVqwAAAQHB8vKk5OTMXbsWADAokWLYGZmhsjISOh0OoSFhWH58uXG7goREVGdZ/SgFkJUWsfKygrLli3DsmXLjL17IiKieoXP+iYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERApmUdsdICLj8Jmxrba7QEQmwBE1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESmY0YP6wIEDGDx4MDw8PKBSqbB582bZeiEE3nrrLbi7u8Pa2hqhoaE4d+6csbtBRERULxg9qAsKCtCpUycsW7bM4PoFCxZgyZIlWLlyJY4cOQJbW1uEhYXh7t27xu4KERFRnWf0+6jDw8MRHh5ucJ0QAosXL8bMmTPx1FNPAQA+//xzuLq6YvPmzRg5cqSxu0NERFSnPdRr1FlZWcjOzkZoaKhUptFoEBgYiLS0tHK30+l00Gq1shcREVFD8FCDOjs7GwDg6uoqK3d1dZXWGZKUlASNRiO9PD09TdpPIiIipagTs74TExORl5cnvS5fvlzbXSIiInooHmpQu7m5AQBycnJk5Tk5OdI6Q9RqNRwcHGQvIiKihuChBrWvry/c3Nywe/duqUyr1eLIkSMICgp6mF0hIiKqE4w+6zs/Px+ZmZnSclZWFk6cOAEnJyd4eXlh6tSpmDdvHlq2bAlfX1+8+eab8PDwwNChQ43dFSIiojrP6EF97NgxhISESMsJCQkAgKioKKxevRrTpk1DQUEBYmJikJubi169eiE1NRVWVlbG7goREVGdpxJCiNruRHVptVpoNBrk5eXxejXR/6e076O+MH9QbXeBSLGqk2N1YtY3ERFRQ8WgJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBTM6N+eRUQE6H9JCL+kg6hmOKImIiJSMAY1ERGRgjGoiYiIFIxBTUREpGCcTEZEDwUnlxHVDEfURERECsagJiIiUjAGNRERkYIxqImIiBSMk8mI6qj7J2cRUf3EETUREZGCMaiJiIgUjEFNRESkYAxqIiIiBeNkMqI6or5NHuOTyoiqhiNqIiIiBWNQExERKRiDmoiISMEY1ERERArGyWREpAicXEZkGEfURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFq7WgXrZsGXx8fGBlZYXAwED8+OOPtdUVIiIixaqVoP7yyy+RkJCAWbNm4fjx4+jUqRPCwsJw7dq12ugOERGRYtXKA08WLlyIiRMnYty4cQCAlStXYtu2bfj0008xY8YMvfo6nQ46nU5azsvLAwBotdqH02EiBSjV3a7tLjxU/PdN9VnZ51sIUXll8ZDpdDphbm4uvv76a1n5888/L4YMGWJwm1mzZgkAfPHFF1988VWvXpcvX640Nx/6iPqvv/5CSUkJXF1dZeWurq44e/aswW0SExORkJAgLZeWluLmzZtwdnaGSqV64D5ptVp4enri8uXLcHBweOD2agOPQTnqw3HwGJSjPhwHj0GfEAK3bt2Ch4dHpXXrxLO+1Wo11Gq1rMzR0dHo+3FwcKizH6IyPAblqA/HwWNQjvpwHDwGOY1GU6V6D30yWZMmTWBubo6cnBxZeU5ODtzc3B52d4iIiBTtoQe1paUlAgICsHv3bqmstLQUu3fvRlBQ0MPuDhERkaLVyqnvhIQEREVFoWvXrujevTsWL16MgoICaRb4w6ZWqzFr1iy90+t1CY9BOerDcfAYlKM+HAeP4cGohKjK3HDj++CDD/Dee+8hOzsbnTt3xpIlSxAYGFgbXSEiIlKsWgtqIiIiqhyf9U1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoyyMfHB08++WRtd6PWqFQqzJ49u0bb7tu3DyqVCvv27ZPKxo4dCx8fH6P0rTI+Pj4YO3astLx69WqoVCocO3bsoew/ODgYwcHBD2Vfhqxfvx5OTk7Iz8+vtT7UR7Nnz4ZKpTJJ20VFRfD09MTy5ctN0n5dx6BWsLJfsPe+XFxcEBISgu3bt9d298jEzpw5g9mzZ+PChQu13RU9Su1bSUkJZs2ahcmTJ8POzk4qLyoqwpw5c9C8eXOo1Wo0b94c8+bNQ3FxsV4b6enpGDBgABwcHGBvb4/+/fvjxIkTsjoXLlzQ+7d572vixIlS3bFjx1ZY988//zTZ+1Gef/3rX9i8efND3295GjVqhISEBLzzzju4e/dubXdHcSxquwNUublz58LX1xdCCOTk5GD16tUYOHAgvvnmmwY96q1LPv74Y5SWllZrmzNnzmDOnDkIDg6u1mg8IyMDZmam/Ru8or599913Jt13Rb755htkZGQgJiZGVv7ss89iw4YNGD9+PLp27YrDhw/jzTffxKVLl/DRRx9J9Y4fP45evXrB09MTs2bNQmlpKZYvX44+ffrgxx9/ROvWrQEATZs2xZo1a/T2n5qairVr16J///5S2QsvvIDQ0FBZPSEEJk2aBB8fHzzyyCPGfAuq5F//+heGDRuGoUOHVnmbmTNnYsaMGSbr07hx4zBjxgykpKRg/PjxJttPXcSgrgPCw8PRtWtXaTk6Ohqurq74v//7vzod1AUFBbC1ta3tbjwUjRo1Mmn7QgjcvXsX1tbWUKvVJt1XZSwtLWtt38nJyejZs6cs/I4ePYr169fjzTffxNy5cwEAkyZNQpMmTbBw4ULExcWhY8eOAIA333wT1tbWSEtLg7OzM4B/Qr5Vq1Z4/fXX8dVXXwEAbG1t8eyzz+rtf/Xq1XBwcMDgwYOlsqCgIAQFBcnqHTx4ELdv38aYMWOM+waYQNm/UwsLC1hYmC4yHB0d0b9/f6xevZpBfR+e+q6DHB0dYW1trfeP5v3338djjz0GZ2dnWFtbIyAgABs3bjTYxhdffIHu3bvDxsYGjRs3Ru/evSsdCX322WewsLDAa6+9JpXduHEDzz33HBwcHODo6IioqCj8/PPPUKlUWL16tVRv7NixsLOzw/nz5zFw4EDY29tLv6QKCgrwyiuvwNPTE2q1Gq1bt8b7778PIYS0fdmpxnvbLHP/9eSya2mZmZkYO3YsHB0dodFoMG7cONy+fVu2rU6nQ3x8PJo2bQp7e3sMGTIEf/zxR4Xvw73++OMPDB06FLa2tnBxcUF8fDx0Op1ePUPXqNetW4eAgADY29vDwcEBHTp0wH/+8x8A//zCHz58OAAgJCREOk1adt27bA7Bjh070LVrV1hbW+PDDz+U1t17jbrM7du38cILL8DZ2RkODg54/vnn8ffff1f4Xpa5t83K+mboGvW1a9ekPzCtrKzQqVMnfPbZZ7I6ZT/j999/Hx999BH8/PygVqvRrVs3HD16VK9P97t79y5SU1P1Rq/ff/89AGDkyJGy8pEjR0IIgS+//FJWNzQ0VAppAHB3d0efPn2wdevWCq97X716FXv37kVERASsrKwq7GtKSgpUKhVGjx5d6XEBwJEjRzBgwABoNBrY2NigT58++OGHH2R1qvq5V6lUKCgowGeffSb97Mp+tmVtnDlzBqNHj0bjxo3Rq1cv2br7ffHFFwgICIC1tTWcnJwwcuRIXL58WVbn3LlziIyMhJubG6ysrNCsWTOMHDkSeXl5snpPPPEEDh48iJs3b1bpfWkoOKKuA/Ly8vDXX39BCIFr165h6dKlyM/P1/uL/j//+Q+GDBmCMWPGoLCwEOvWrcPw4cOxdetWDBo0SKo3Z84czJ49G4899hjmzp0LS0tLHDlyBHv27JGdsrvXRx99hEmTJuH111/HvHnzAAClpaUYPHgwfvzxR7z44ovw9/fHli1bEBUVZbCN4uJihIWFoVevXnj//fdhY2MDIQSGDBmCvXv3Ijo6Gp07d8aOHTvw2muv4c8//8SiRYtq/L6NGDECvr6+SEpKwvHjx/HJJ5/AxcUF7777rlRnwoQJ+OKLLzB69Gg89thj2LNnj+y9qsidO3fQr18/XLp0CVOmTIGHhwfWrFmDPXv2VLrtzp07MWrUKPTr10/qz6+//ooffvgBL7/8Mnr37o0pU6ZgyZIleP3119GmTRsAkP4L/HOKe9SoUXjhhRcwceJE6bRseeLi4uDo6IjZs2cjIyMDK1aswMWLF6XJb1VVlb7d686dOwgODkZmZibi4uLg6+uLDRs2YOzYscjNzcXLL78sq5+SkoJbt27hhRdegEqlwoIFCxAREYHff/+9wjMT6enpKCwsRJcuXWTlZX84WVtby8ptbGyk7e6te3+9srqFhYU4ffo0evToYXD/69atQ2lpaaWj5KKiIqxfvx6PPfZYlS5p7NmzB+Hh4QgICMCsWbNgZmaG5ORk9O3bF99//z26d+8uq1/Z537NmjWYMGECunfvLl0i8PPzk7UxfPhwtGzZEv/6179kfzDf75133sGbb76JESNGYMKECbh+/TqWLl2K3r1746effoKjoyMKCwsRFhYGnU6HyZMnw83NDX/++Se2bt2K3NxcaDQaqb2AgAAIIXDo0KE6fbbQ6AQpVnJysgCg91Kr1WL16tV69W/fvi1bLiwsFO3btxd9+/aVys6dOyfMzMzE008/LUpKSmT1S0tLpf/39vYWgwYNEkII8Z///EeoVCrx9ttvy+p/9dVXAoBYvHixVFZSUiL69u0rAIjk5GSpPCoqSgAQM2bMkLWxefNmAUDMmzdPVj5s2DChUqlEZmamEEKIrKwsvTbLABCzZs2SlmfNmiUAiPHjx8vqPf3008LZ2VlaPnHihAAgXnrpJVm90aNH67VpyOLFiwUAsX79eqmsoKBAtGjRQgAQe/fulR2/t7e3tPzyyy8LBwcHUVxcXG77GzZs0GunjLe3twAgUlNTDa6LioqSlss+RwEBAaKwsFAqX7BggQAgtmzZIpWVd9z3t1lR3/r06SP69OkjLZe9T1988YVUVlhYKIKCgoSdnZ3QarVCiP/9jJ2dncXNmzelulu2bBEAxDfffKO3r3t98sknAoA4deqUrLzsc7pmzRpZ+cqVKwUA0b59e6msQ4cOolWrVrKfi06nE15eXgKA2LhxY7n7DwgIEO7u7nr/ru73zTffCABi+fLlFdYT4p9/ky1bthRhYWGyf5+3b98Wvr6+4oknnpDKqvq5F0IIW1tb2c/z/jZGjRpV7royFy5cEObm5uKdd96R1Tt16pSwsLCQyn/66ScBQGzYsKHS471y5YoAIN59991K6zYkPPVdByxbtgw7d+7Ezp078cUXXyAkJAQTJkzApk2bZPXuHQn8/fffyMvLw+OPP47jx49L5Zs3b0ZpaSneeustvQlHhkZVCxYswMsvv4x3330XM2fOlK1LTU1Fo0aNZDNczczMEBsbW+6xvPjii7Llb7/9Fubm5pgyZYqs/JVXXoEQ4oFmt0+aNEm2/Pjjj+PGjRvQarXSvgHo7Xvq1KlVav/bb7+Fu7s7hg0bJpXZ2NjoTWQyxNHREQUFBdi5c2eV9mWIr68vwsLCqlw/JiZGNiJ98cUXYWFhIb0PpvLtt9/Czc0No0aNksoaNWqEKVOmID8/H/v375fVf+aZZ9C4cWNp+fHHHwcA/P777xXu58aNGwAg2xYABg4cCG9vb7z66qvYtGkTLl68iPXr1+ONN96AhYUF7ty5I9V96aWX8NtvvyE6OhpnzpzB6dOn8fzzz+Pq1asAIKt7r99++w3p6ekYOXJkpRP5UlJS0KhRI4wYMaLCegBw4sQJnDt3DqNHj8aNGzfw119/4a+//kJBQQH69euHAwcO6E1SrOxzXxX3t2HIpk2bUFpaihEjRkj9+uuvv+Dm5oaWLVti7969ACCNmHfs2KF36el+ZT+7v/76q8p9bQh46rsO6N69u2wy2ahRo/Doo48iLi4OTz75pDR5Z+vWrZg3bx5OnDghu056bwCfP38eZmZmaNu2baX73b9/P7Zt24bp06fLrkuXuXjxItzd3aVTiGVatGhhsD0LCws0a9ZMrw0PDw/Y29vLystOo168eLHSfpbHy8tLtlz2S+Dvv/+Gg4MDLl68CDMzM73TfpWdQi5z8eJFtGjRQu8PnKps/9JLL2H9+vUIDw/HI488gv79+2PEiBEYMGBAlfYN/BPU1dGyZUvZsp2dHdzd3U1+i9XFixfRsmVLvQAr72dc0c+tKsR9p2qtrKywbds2jBgxApGRkQAAtVqNBQsW4J133pHdxjVp0iRcvnwZ7733nnQNvWvXrpg2bZpe3XutXbsWACo97Z2fn48tW7YgLCxMdh28POfOnQOAci8nAf9cGrv3j5PKPvdVUZXP1rlz5yCE0PtclSn7o9DX1xcJCQlYuHAh1q5di8cffxxDhgzBs88+KzvtDfzvZ2eq+7XrKgZ1HWRmZoaQkBD85z//wblz59CuXTt8//33GDJkCHr37o3ly5fD3d0djRo1QnJyMlJSUmq0n3bt2iE3Nxdr1qzBCy+8UO1guJ9ara7xbUPl/cMtKSkpdxtzc3OD5ff/Iq8NLi4uOHHiBHbs2IHt27dj+/btSE5OxvPPP683yao8hq6lmkpF77Ox1fTnVhZ8f//9t94fhO3atcPp06dx5swZ/P3332jbti2sra0RHx+PPn36yOq+8847ePXVV/HLL79Ao9GgQ4cOeP311wEArVq1MrjvlJQUtG7dGgEBARX2cfPmzdWa7V02Wn7vvffQuXNng3Xu/+PBGJ/7qny2SktLoVKpsH37doP7vLdf//73vzF27Fhs2bIF3333HaZMmYKkpCQcPnxY9rMq+2OsSZMmVe5rQ8CgrqPKHtRQNgv1q6++gpWVFXbs2CG7PSc5OVm2nZ+fH0pLS3HmzJly/+GXadKkCTZu3IhevXqhX79+OHjwIDw8PKT13t7e2Lt3L27fvi0bVWdmZlb5OLy9vbFr1y7cunVLNqo+e/astB7436ggNzdXtv2DjLi9vb1RWlqK8+fPy0bBGRkZVd7+9OnTEELI/pCo6vaWlpYYPHgwBg8ejNLSUrz00kv48MMP8eabbxocqT+oc+fOISQkRFrOz8/H1atXMXDgQKmscePGeu9xYWGhdOq3THX65u3tjZMnT6K0tFT2h9r9P+MH5e/vDwDIyspChw4d9NarVCq0a9dOWv72229RWlqqN0scgGy2MwDs2rULzZo1k/ZxryNHjiAzM1O69asia9euhZ2dHYYMGVKlYyo72+Pg4GCwnzVljM+Wn58fhBDw9fUt9w+Ye3Xo0AEdOnTAzJkzcejQIfTs2RMrV66UJqcC//zsgPInJjZUvEZdBxUVFeG7776DpaWl9IE2NzeHSqWSjXwuXLig9/ShoUOHwszMDHPnztW7tmXoL+5mzZph165duHPnDp544gnpOiAAhIWFoaioCB9//LFUVlpaimXLllX5WAYOHIiSkhJ88MEHsvJFixZBpVIhPDwcwD+/qJo0aYIDBw7I6j3IIwfL2l6yZImsfPHixVXafuDAgbhy5YrsFrjbt2/LHqBRnnvfR+CfsyRl9/KWXbYou8f8/uCsqY8++ghFRUXS8ooVK1BcXCy9D8A/v3zvf48/+ugjvRF1dfo2cOBAZGdny26DKi4uxtKlS2FnZ6c3oq2pgIAAWFpaVulRqXfu3MGbb74Jd3d32bVzQ7788kscPXoUU6dONXhGqOyMVWW3Wl2/fh27du3C008/rXe5qDwBAQHw8/PD+++/b/DWsOvXr1epnfvZ2to+8OcqIiIC5ubmmDNnjt7vDiGE9BnXarV6T4Dr0KEDzMzM9G5lTE9Ph0ql0rvvvKHjiLoO2L59uzT6uHbtGlJSUnDu3DnMmDFDuuY0aNAgLFy4EAMGDMDo0aNx7do1LFu2DC1atMDJkyeltlq0aIE33ngDb7/9Nh5//HFERERArVbj6NGj8PDwQFJSkt7+W7Roge+++w7BwcEICwvDnj174ODggKFDh6J79+545ZVXkJmZCX9/f/z3v/+V7oGsyl/tgwcPRkhICN544w1cuHABnTp1wnfffYctW7Zg6tSpsuvHEyZMwPz58zFhwgR07doVBw4cwG+//Vbj97Vz584YNWoUli9fjry8PDz22GPYvXt3lc8ITJw4ER988AGef/55pKenw93dHWvWrKnSL+EJEybg5s2b6Nu3L5o1a4aLFy9i6dKl6Ny5s/THV+fOnWFubo53330XeXl5UKvV6Nu3L1xcXGp0vIWFhejXrx9GjBiBjIwMLF++HL169ZKN7iZMmIBJkyYhMjISTzzxBH7++Wfs2LFD71RkdfoWExODDz/8EGPHjkV6ejp8fHywceNG/PDDD1i8eLHe/ISasrKyQv/+/bFr1y690e2IESPg4eGBtm3bQqvV4tNPP8Xvv/+Obdu2yfZ/4MABzJ07F/3794ezszMOHz6M5ORkDBgwQO82MuCfSwJffvklevTooTfX4X5ffvkliouLq/WQEzMzM3zyyScIDw9Hu3btMG7cODzyyCP4888/sXfvXjg4OOCbb76pcntlAgICsGvXLixcuBAeHh7w9fVFYGBgtdrw8/PDvHnzkJiYiAsXLmDo0KGwt7dHVlYWvv76a8TExODVV1/Fnj17EBcXh+HDh6NVq1YoLi7GmjVrYG5uLs0ZKLNz50707NmzStfvG5TamWxOVWHo9iwrKyvRuXNnsWLFCtntGkIIsWrVKtGyZUuhVquFv7+/SE5O1rulosynn34qHn30UaFWq0Xjxo1Fnz59xM6dO6X1996eVebIkSPC3t5e9O7dW7oV7Pr162L06NHC3t5eaDQaMXbsWPHDDz8IAGLdunXStlFRUcLW1tbgcd66dUvEx8cLDw8P0ahRI9GyZUvx3nvv6R3f7du3RXR0tNBoNMLe3l6MGDFCXLt2rdzbs65fv27w/czKypLK7ty5I6ZMmSKcnZ2Fra2tGDx4sLh8+XKVbs8SQoiLFy+KIUOGCBsbG9GkSRPx8ssvi9TU1Epvz9q4caPo37+/cHFxEZaWlsLLy0u88MIL4urVq7L2P/74Y9G8eXNhbm4ua9PQz6dMebdn7d+/X8TExIjGjRsLOzs7MWbMGHHjxg3ZtiUlJWL69OmiSZMmwsbGRoSFhYnMzEy9Nivq2/23ZwkhRE5Ojhg3bpxo0qSJsLS0FB06dNC71a7s9qz33ntP75iq+vPYtGmTUKlU4tKlS7Lyd999V/j7+wsrKyvRuHFjMWTIEPHTTz/pbZ+ZmSn69+8vmjRpIv07SkpKEjqdzuD+yn7WS5YsqbRvPXr0EC4uLhXekleen376SURERAhnZ2ehVquFt7e3GDFihNi9e7dUpzqf+7Nnz4revXsLa2trAUD62ZbXxr3r7vfVV1+JXr16CVtbW2Frayv8/f1FbGysyMjIEEII8fvvv4vx48cLPz8/YWVlJZycnERISIjYtWuXrJ3c3FxhaWkpPvnkk2q/P/WdSggFzKyhemXz5s14+umncfDgQfTs2bO2u0MNSElJCdq2bYsRI0bg7bffru3uUDUsXrwYCxYswPnz5x/qRMm6gEFND+TOnTuyf1QlJSXo378/jh07huzsbP6Do4fuyy+/xIsvvohLly6VezsVKUtRURH8/PwwY8YMvPTSS7XdHcVhUNMDmTBhAu7cuYOgoCDodDps2rQJhw4dwr/+9S8kJibWdveIiOo8BjU9kJSUFPz73/9GZmYm7t69ixYtWuDFF19EXFxcbXeNiKheYFATEREpGO+jJiIiUrA6eR91aWkprly5Ant7ez4TloiI6hwhBG7dugUPD49KH61cJ4P6ypUr8PT0rO1uEBERPZDLly/rPZv+fnUyqMueJHT58uUqfxsMERGRUmi1Wnh6elbpyXx1MqjLTnc7ODgwqImIqM6qyuVbTiYjIiJSMAY1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKxqAmIiJSMAY1ERGRgtXJR4gqgc+MbdL/X5g/qBZ7QkRE9RlH1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUrFpBnZSUhG7dusHe3h4uLi4YOnQoMjIyZHXu3r2L2NhYODs7w87ODpGRkcjJyZHVuXTpEgYNGgQbGxu4uLjgtddeQ3Fx8YMfjQn5zNgmexERET0M1Qrq/fv3IzY2FocPH8bOnTtRVFSE/v37o6CgQKoTHx+Pb775Bhs2bMD+/ftx5coVRERESOtLSkowaNAgFBYW4tChQ/jss8+wevVqvPXWW8Y7KiIionpCJYQQNd34+vXrcHFxwf79+9G7d2/k5eWhadOmSElJwbBhwwAAZ8+eRZs2bZCWloYePXpg+/btePLJJ3HlyhW4uroCAFauXInp06fj+vXrsLS0rHS/Wq0WGo0GeXl5cHBwqGn3q6WiUTS/5pKIiKqjOjn2QNeo8/LyAABOTk4AgPT0dBQVFSE0NFSq4+/vDy8vL6SlpQEA0tLS0KFDBymkASAsLAxarRa//PKLwf3odDpotVrZi4iIqCGocVCXlpZi6tSp6NmzJ9q3bw8AyM7OhqWlJRwdHWV1XV1dkZ2dLdW5N6TL1petMyQpKQkajUZ6eXp61rTbREREdYpFTTeMjY3F6dOncfDgQWP2x6DExEQkJCRIy1qt1uRhzQljRESkBDUK6ri4OGzduhUHDhxAs2bNpHI3NzcUFhYiNzdXNqrOycmBm5ubVOfHH3+UtVc2K7yszv3UajXUanVNukpERFSnVevUtxACcXFx+Prrr7Fnzx74+vrK1gcEBKBRo0bYvXu3VJaRkYFLly4hKCgIABAUFIRTp07h2rVrUp2dO3fCwcEBbdu2fZBjISIiqneqNaKOjY1FSkoKtmzZAnt7e+maskajgbW1NTQaDaKjo5GQkAAnJyc4ODhg8uTJCAoKQo8ePQAA/fv3R9u2bfHcc89hwYIFyM7OxsyZMxEbG8tRMxER0X2qFdQrVqwAAAQHB8vKk5OTMXbsWADAokWLYGZmhsjISOh0OoSFhWH58uVSXXNzc2zduhUvvvgigoKCYGtri6ioKMydO/fBjoSIiKgeeqD7qGvLw7iPujqTyXgfNRERVcdDu4+aiIiITItBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBavwIUfqf+2eIcxY4EREZC0fURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIF4+1ZdRBvByMiajgY1ApQWfBW55u8iIiofuGpbyIiIgXjiLoeqGjEzdPiRER1G0fURERECsYR9f/H68BERKREHFETEREpWLVH1AcOHMB7772H9PR0XL16FV9//TWGDh0qrVepVAa3W7BgAV577TUAgI+PDy5evChbn5SUhBkzZlS3O3UCR+tERFRT1R5RFxQUoFOnTli2bJnB9VevXpW9Pv30U6hUKkRGRsrqzZ07V1Zv8uTJNTsCIiKieqzaI+rw8HCEh4eXu97NzU22vGXLFoSEhKB58+aycnt7e7269PDx4SlERMpm0mvUOTk52LZtG6Kjo/XWzZ8/H87Oznj00Ufx3nvvobi4uNx2dDodtFqt7EVERNQQmHTW92effQZ7e3tERETIyqdMmYIuXbrAyckJhw4dQmJiIq5evYqFCxcabCcpKQlz5swxZVfp/+MIm4hIWUwa1J9++inGjBkDKysrWXlCQoL0/x07doSlpSVeeOEFJCUlQa1W67WTmJgo20ar1cLT09N0HSciIlIIkwX1999/j4yMDHz55ZeV1g0MDERxcTEuXLiA1q1b661Xq9UGA5yIiKi+M9k16lWrViEgIACdOnWqtO6JEydgZmYGFxcXU3WHiIioTqr2iDo/Px+ZmZnSclZWFk6cOAEnJyd4eXkB+OfU9IYNG/Dvf/9bb/u0tDQcOXIEISEhsLe3R1paGuLj4/Hss8+icePGD3AoRERE9U+1g/rYsWMICQmRlsuuHUdFRWH16tUAgHXr1kEIgVGjRultr1arsW7dOsyePRs6nQ6+vr6Ij4+XXYMmIiKif1Q7qIODgyGEqLBOTEwMYmJiDK7r0qULDh8+XN3dNih8khkREZXhs76JiIgUjN+eRSZ179kB3pNNRFR9HFETEREpGIOaiIhIwXjq2wQ4GYyIiIyFI2oiIiIF44iaag2/AISIqHIcURMRESkYg5qIiEjBGNREREQKxmvUZFSc8U5EZFwcURMRESkYR9RULZypTUT0cDGo6zkGKxFR3cZT30RERArGoCYiIlIwnvqmB8JZ3kREpsURNRERkYIxqImIiBSMQU1ERKRg1Q7qAwcOYPDgwfDw8IBKpcLmzZtl68eOHQuVSiV7DRgwQFbn5s2bGDNmDBwcHODo6Ijo6Gjk5+c/0IEQERHVR9WeTFZQUIBOnTph/PjxiIiIMFhnwIABSE5OlpbVarVs/ZgxY3D16lXs3LkTRUVFGDduHGJiYpCSklLd7lADwnvCiaghqnZQh4eHIzw8vMI6arUabm5uBtf9+uuvSE1NxdGjR9G1a1cAwNKlSzFw4EC8//778PDw0NtGp9NBp9NJy1qttrrdpv+Ps7SJiOoWk1yj3rdvH1xcXNC6dWu8+OKLuHHjhrQuLS0Njo6OUkgDQGhoKMzMzHDkyBGD7SUlJUGj0UgvT09PU3SbiIhIcYx+H/WAAQMQEREBX19fnD9/Hq+//jrCw8ORlpYGc3NzZGdnw8XFRd4JCws4OTkhOzvbYJuJiYlISEiQlrVaLcP6IeEInIiodhk9qEeOHCn9f4cOHdCxY0f4+flh37596NevX43aVKvVete5iYiIGgKT357VvHlzNGnSBJmZmQAANzc3XLt2TVanuLgYN2/eLPe6NhERUUNl8keI/vHHH7hx4wbc3d0BAEFBQcjNzUV6ejoCAgIAAHv27EFpaSkCAwNN3R2iGuGMcyKqLdUO6vz8fGl0DABZWVk4ceIEnJyc4OTkhDlz5iAyMhJubm44f/48pk2bhhYtWiAsLAwA0KZNGwwYMAATJ07EypUrUVRUhLi4OIwcOdLgjG+qP3i9m4io+qp96vvYsWN49NFH8eijjwIAEhIS8Oijj+Ktt96Cubk5Tp48iSFDhqBVq1aIjo5GQEAAvv/+e9k15rVr18Lf3x/9+vXDwIED0atXL3z00UfGOyoiIqJ6otoj6uDgYAghyl2/Y8eOSttwcnLiw02IiIiqgM/6JiIiUjB+HzUpFq9pExFxRE1ERKRoHFFTvcVbqoioPmBQU53FICaihoBBTfUGr2kTUX3Ea9REREQKxqAmIiJSMJ76pgaD17SJqC7iiJqIiEjBGNREREQKxlPfpBictU1EpI8jaiIiIgXjiJoInGhGRMrFoCYyAQY/ERkLT30TEREpGIOaiIhIwRjURERECsagJiIiUrBqB/WBAwcwePBgeHh4QKVSYfPmzdK6oqIiTJ8+HR06dICtrS08PDzw/PPP48qVK7I2fHx8oFKpZK/58+c/8MEQ1RafGdtkLyIiY6n2rO+CggJ06tQJ48ePR0REhGzd7du3cfz4cbz55pvo1KkT/v77b7z88ssYMmQIjh07Jqs7d+5cTJw4UVq2t7ev4SEQNWycYU5Uv1U7qMPDwxEeHm5wnUajwc6dO2VlH3zwAbp3745Lly7By8tLKre3t4ebm1uV9qnT6aDT6aRlrVZb3W4TERHVSSa/jzovLw8qlQqOjo6y8vnz5+Ptt9+Gl5cXRo8ejfj4eFhYGO5OUlIS5syZY+quEkl4+pqIlMKkQX337l1Mnz4do0aNgoODg1Q+ZcoUdOnSBU5OTjh06BASExNx9epVLFy40GA7iYmJSEhIkJa1Wi08PT1N2XUiIiJFMFlQFxUVYcSIERBCYMWKFbJ194Zux44dYWlpiRdeeAFJSUlQq9V6banVaoPlRERE9Z1JgrospC9evIg9e/bIRtOGBAYGori4GBcuXEDr1q1N0SUiReOEMCIqj9GDuiykz507h71798LZ2bnSbU6cOAEzMzO4uLgYuztERER1WrWDOj8/H5mZmdJyVlYWTpw4AScnJ7i7u2PYsGE4fvw4tm7dipKSEmRnZwMAnJycYGlpibS0NBw5cgQhISGwt7dHWloa4uPj8eyzz6Jx48bGOzIiBeNkNSKqqmoH9bFjxxASEiItl11vjoqKwuzZs/Hf//4XANC5c2fZdnv37kVwcDDUajXWrVuH2bNnQ6fTwdfXF/Hx8bLr1kRERPSPagd1cHAwhBDlrq9oHQB06dIFhw8fru5uiYiIGiR+HzU1WDz9TER1AYOaiCp07x80nI1O9PDx27OIiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISME465voIeCtYERUUxxRExERKRiDmoiISMEY1ERERArGoCYiIlIwTiYjqmfun7jGx34S1W0MaqIaeJizuBm8RA0bg5qogWHwE9UtvEZNRESkYBxRE9FDw9E8UfVxRE1ERKRgDGoiIiIFq3ZQHzhwAIMHD4aHhwdUKhU2b94sWy+EwFtvvQV3d3dYW1sjNDQU586dk9W5efMmxowZAwcHBzg6OiI6Ohr5+fkPdCBEDYXPjG2yFxHVb9UO6oKCAnTq1AnLli0zuH7BggVYsmQJVq5ciSNHjsDW1hZhYWG4e/euVGfMmDH45ZdfsHPnTmzduhUHDhxATExMzY+CiIionqr2ZLLw8HCEh4cbXCeEwOLFizFz5kw89dRTAIDPP/8crq6u2Lx5M0aOHIlff/0VqampOHr0KLp27QoAWLp0KQYOHIj3338fHh4eeu3qdDrodDppWavVVrfbREREdZJRr1FnZWUhOzsboaGhUplGo0FgYCDS0tIAAGlpaXB0dJRCGgBCQ0NhZmaGI0eOGGw3KSkJGo1Genl6ehqz20RERIpl1KDOzs4GALi6usrKXV1dpXXZ2dlwcXGRrbewsICTk5NU536JiYnIy8uTXpcvXzZmt4mIiBSrTtxHrVaroVara7sbRA8NJ4kRURmjBrWbmxsAICcnB+7u7lJ5Tk4OOnfuLNW5du2abLvi4mLcvHlT2p6IjIehT1S3GTWofX194ebmht27d0vBrNVqceTIEbz44osAgKCgIOTm5iI9PR0BAQEAgD179qC0tBSBgYHG7A4R1QCDnUhZqh3U+fn5yMzMlJazsrJw4sQJODk5wcvLC1OnTsW8efPQsmVL+Pr64s0334SHhweGDh0KAGjTpg0GDBiAiRMnYuXKlSgqKkJcXBxGjhxpcMY3ERFRQ1btoD527BhCQkKk5YSEBABAVFQUVq9ejWnTpqGgoAAxMTHIzc1Fr169kJqaCisrK2mbtWvXIi4uDv369YOZmRkiIyOxZMkSIxwOEdUmPsubyPiqHdTBwcEQQpS7XqVSYe7cuZg7d265dZycnJCSklLdXRMRETU4dWLWNxEpA0fMRA8fv5SDiIhIwRjURERECsagJiIiUjBeoyaiGuM910SmxxE1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgXjZDIiqjXVnYzGB6xQQ8QRNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGGd9E5HJ8BGjRA+OI2oiIiIF44iaqIHjqJdI2RjURFQv3P8Hx/0PR6lsPZFSGf3Ut4+PD1Qqld4rNjYWABAcHKy3btKkScbuBhERUb1g9BH10aNHUVJSIi2fPn0aTzzxBIYPHy6VTZw4EXPnzpWWbWxsjN0NIiKiesHoQd20aVPZ8vz58+Hn54c+ffpIZTY2NnBzczP2romIiOodk876LiwsxBdffIHx48dDpVJJ5WvXrkWTJk3Qvn17JCYm4vbt2xW2o9PpoNVqZS8iIqKGwKSTyTZv3ozc3FyMHTtWKhs9ejS8vb3h4eGBkydPYvr06cjIyMCmTZvKbScpKQlz5swxZVeJiB4IJ6uRqZg0qFetWoXw8HB4eHhIZTExMdL/d+jQAe7u7ujXrx/Onz8PPz8/g+0kJiYiISFBWtZqtfD09DRdx4mIiBTCZEF98eJF7Nq1q8KRMgAEBgYCADIzM8sNarVaDbVabfQ+EhERKZ3Jgjo5ORkuLi4YNKji0z8nTpwAALi7u5uqK0TUAPFBLlRfmCSoS0tLkZycjKioKFhY/G8X58+fR0pKCgYOHAhnZ2ecPHkS8fHx6N27Nzp27GiKrhAREdVpJgnqXbt24dKlSxg/frys3NLSErt27cLixYtRUFAAT09PREZGYubMmaboBhERUZ1nkqDu378/hBB65Z6enti/f78pdklERFQv8duziIiIFIxBTUREpGAMaiIiIgVjUBMRESkYg5qIiEjBTPoIUSKi+uphPlCFzxFv2DiiJiIiUjCOqImozuJjQqkhYFATUYNU3dPJ/KOAaguDmojIAAYzKQWvURMRESkYR9RERDDtCJqztulBcERNRESkYBxRE1GdwevG1BBxRE1ERKRgHFETEZkAR/9kLBxRExERKRhH1ERECsPRON2LQU1EVMsYzFQRo5/6nj17NlQqlezl7+8vrb979y5iY2Ph7OwMOzs7REZGIicnx9jdICIiqhdMco26Xbt2uHr1qvQ6ePCgtC4+Ph7ffPMNNmzYgP379+PKlSuIiIgwRTeIiIjqPJOc+rawsICbm5teeV5eHlatWoWUlBT07dsXAJCcnIw2bdrg8OHD6NGjhym6Q0REVGeZZER97tw5eHh4oHnz5hgzZgwuXboEAEhPT0dRURFCQ0Oluv7+/vDy8kJaWlq57el0Omi1WtmLiIioITB6UAcGBmL16tVITU3FihUrkJWVhccffxy3bt1CdnY2LC0t4ejoKNvG1dUV2dnZ5baZlJQEjUYjvTw9PY3dbSIiIkUy+qnv8PBw6f87duyIwMBAeHt7Y/369bC2tq5Rm4mJiUhISJCWtVotw5qIiBoEkz/wxNHREa1atUJmZibc3NxQWFiI3NxcWZ2cnByD17TLqNVqODg4yF5EREQNgcnvo87Pz8f58+fx3HPPISAgAI0aNcLu3bsRGRkJAMjIyMClS5cQFBRk6q4QEdULld13za/RrF+MHtSvvvoqBg8eDG9vb1y5cgWzZs2Cubk5Ro0aBY1Gg+joaCQkJMDJyQkODg6YPHkygoKCOOObiIjIAKMH9R9//IFRo0bhxo0baNq0KXr16oXDhw+jadOmAIBFixbBzMwMkZGR0Ol0CAsLw/Lly43dDSIionrB6EG9bt26CtdbWVlh2bJlWLZsmbF3TUREVO/w27OIiIgUjF/KQUT0kPFLOKg6OKImIiJSMAY1ERGRgjGoiYiIFIxBTUREpGCcTEZEVM/dP3mNTy6rWziiJiIiUjAGNRERkYIxqImIiBSM16iJiBowXr9WPo6oiYiIFIxBTUREpGA89U1E1MCY8lnjPJVufAxqIqJ6hl/6Ub/w1DcREZGCcURNRESSykbjPJX98HFETUREpGAMaiIiIgUzelAnJSWhW7dusLe3h4uLC4YOHYqMjAxZneDgYKhUKtlr0qRJxu4KEREpjM+MbbIXVc7o16j379+P2NhYdOvWDcXFxXj99dfRv39/nDlzBra2tlK9iRMnYu7cudKyjY2NsbtCRES1jGH84Iwe1KmpqbLl1atXw8XFBenp6ejdu7dUbmNjAzc3N2PvnoiIqF4x+TXqvLw8AICTk5OsfO3atWjSpAnat2+PxMRE3L59u9w2dDodtFqt7EVERNQQmPT2rNLSUkydOhU9e/ZE+/btpfLRo0fD29sbHh4eOHnyJKZPn46MjAxs2rTJYDtJSUmYM2eOKbtKRESkSCYN6tjYWJw+fRoHDx6UlcfExEj/36FDB7i7u6Nfv344f/48/Pz89NpJTExEQkKCtKzVauHp6Wm6jhMRUa2495o279n+h8mCOi4uDlu3bsWBAwfQrFmzCusGBgYCADIzMw0GtVqthlqtNkk/iYiIlMzoQS2EwOTJk/H1119j37598PX1rXSbEydOAADc3d2N3R0iIjIizuJ++Iwe1LGxsUhJScGWLVtgb2+P7OxsAIBGo4G1tTXOnz+PlJQUDBw4EM7Ozjh58iTi4+PRu3dvdOzY0djdISIiqtOMHtQrVqwA8M9DTe6VnJyMsWPHwtLSErt27cLixYtRUFAAT09PREZGYubMmcbuChERUZ1nklPfFfH09MT+/fuNvVsiIqJ6ic/6JiIiUjB+zSURESnS/RPXGurtWhxRExERKRiDmoiISMEY1ERERArGoCYiIlIwTiYjIqI6oaFOLuOImoiISMEY1ERERArGU99ERFQvVPcLQ+rKqXMGNRER1ZoH+TauhvJNXjz1TUREpGAMaiIiIgVjUBMRESkYg5qIiEjBGNREREQKxqAmIiJSMAY1ERGRgvE+aiIiIij3WeK1NqJetmwZfHx8YGVlhcDAQPz444+11RUiIiLFqpWg/vLLL5GQkIBZs2bh+PHj6NSpE8LCwnDt2rXa6A4REZFiqYQQ4mHvNDAwEN26dcMHH3wAACgtLYWnpycmT56MGTNmVLq9VquFRqNBXl4eHBwcjNKnhvIoOiIiqhpTnvquTo499GvUhYWFSE9PR2JiolRmZmaG0NBQpKWlGdxGp9NBp9NJy3l5eQD+OVBjKdXdNlpbRERU9xkzY8pruypj5Yce1H/99RdKSkrg6uoqK3d1dcXZs2cNbpOUlIQ5c+bolXt6epqkj0RERJrFpt/HrVu3oNFoKqxTJ2Z9JyYmIiEhQVouLS3FzZs34ezsDJVK9cDta7VaeHp64vLly0Y7lf6w8RiUoz4cB49BOerDcfAY9AkhcOvWLXh4eFRa96EHdZMmTWBubo6cnBxZeU5ODtzc3Axuo1aroVarZWWOjo5G75uDg0Od/RCV4TEoR304Dh6DctSH4+AxyFU2ki7z0Gd9W1paIiAgALt375bKSktLsXv3bgQFBT3s7hARESlarZz6TkhIQFRUFLp27Yru3btj8eLFKCgowLhx42qjO0RERIpVK0H9zDPP4Pr163jrrbeQnZ2Nzp07IzU1VW+C2cOiVqsxa9YsvdPrdQmPQTnqw3HwGJSjPhwHj+HB1Mp91ERERFQ1/FIOIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQkyLs27cPKpUK+/btq9H2s2fPhkqlkpX5+Phg7NixD965Sly4cAEqlQqrV6+WysaOHQs7OzuT77uMSqXC7NmzH9r+7vfSSy/hiSeeqLX9k2k/AytXroSXlxd0Op1J2qeKMagbMJVKVaVXVcLzX//6FzZv3mzyPpvat99+W6uBVxGl9i0rKwuffPIJXn/9dVn5ihUrMHz4cHh5eUGlUlX4R1Nubi5iYmLQtGlT2NraIiQkBMePH9erFx8fjy5dusDJyQk2NjZo06YNZs+ejfz8fL26Op0O06dPh4eHB6ytrREYGIidO3ca3P+hQ4fQq1cv2NjYwM3NDVOmTDHY5sNw6NAhzJ49G7m5ubWyf0PGjh2LwsJCfPjhh7XdlYZJUIO1Zs0a2euJJ54QAPTKs7OzK23L1tZWREVF1bgve/fuFQDE3r17a7T9rFmzxP0f57t374rCwsJqtRMbG6vXTmVKS0vFnTt3RHFxsVQWFRUlbG1tq9XOg/Ttzp07oqioyKj7q6qXX35ZtGrVSq/c29tbODk5iQEDBggLC4tyPx8lJSXiscceE7a2tmL27Nnigw8+EG3bthX29vbit99+k9Xt2bOnmDJliliyZIn46KOPxIsvvijUarXo2bOnKCkpkdUdOXKksLCwEK+++qr48MMPRVBQkLCwsBDff/+9rN5PP/0krKysxKOPPipWrFgh3njjDaFWq8WAAQMe7I2poffee08AEFlZWdXaztSfgWnTpglvb29RWlpqsn2QYQxqktQkpMooMahrojrvQVFRkdDpdAbXPeygri2FhYWiSZMmYubMmXrrLly4IP1Sr+jz8eWXXwoAYsOGDVLZtWvXhKOjoxg1alSlfXj//fcFAJGWliaVHTlyRAAQ7733nlR2584d4efnJ4KCgmTbh4eHC3d3d5GXlyeVffzxxwKA2LFjR6X7N7bqBHVJSYm4c+eO6TslhDh27JgAIHbv3v1Q9kf/w1PfVKGCggK88sor8PT0hFqtRuvWrfH+++9DCCHVUalUKCgowGeffSadLi87zXnx4kW89NJLaN26NaytreHs7Izhw4fjwoULNe7TwYMH0a1bN1hZWcHPz6/c03H3X6MuKirCnDlz0LJlS1hZWcHZ2Rm9evWSToeOHTsWy5Ytk46p7AX87zr0+++/j8WLF8PPzw9qtRpnzpwxeI26zO+//46wsDDY2trCw8MDc+fOlb135V2bv7/NivpWVnb/afGffvoJ4eHhcHBwgJ2dHfr164fDhw/L6qxevRoqlQo//PADEhISpFPPTz/9NK5fv274B3CPgwcP4q+//kJoaKjeOm9vb715A4Zs3LgRrq6uiIiIkMqaNm2KESNGYMuWLZVeF/Xx8QEA2anijRs3wtzcHDExMVKZlZUVoqOjkZaWhsuXLwMAtFotdu7ciWeffRYODg5S3eeffx52dnZYv359pf3X6XSYNWsWWrRoAbVaDU9PT0ybNk2v3yqVCnFxcdi8eTPat28PtVqNdu3aITU1Vaoze/ZsvPbaawAAX19f6edc9u+lrI21a9eiXbt2UKvV0vaGPgN//vknxo8fD1dXV2l/n376qd4xLF26FO3atYONjQ0aN26Mrl27IiUlRVYnICAATk5O2LJlS6XvCRmXRW13gJRLCIEhQ4Zg7969iI6ORufOnbFjxw689tpr+PPPP7Fo0SIAwJo1azBhwgR0795d+sXo5+cHADh69CgOHTqEkSNHolmzZrhw4QJWrFiB4OBgnDlzBjY2NtXq06lTp9C/f380bdoUs2fPRnFxMWbNmgVXV9dKt509ezaSkpKkvmq1Whw7dgzHjx/HE088gRdeeAFXrlzBzp07sWbNGoNtJCcn4+7du4iJiYFarYaTkxNKS0sN1i0pKcGAAQPQo0cPLFiwAKmpqZg1axaKi4sxd+7cah13Vfp2r19++QWPP/44HBwcMG3aNDRq1AgffvghgoODsX//fgQGBsrqT548GY0bN8asWbNw4cIFLF68GHFxcfjyyy8r3M+hQ4egUqnw6KOPVut47vXTTz+hS5cuMDOTjxu6d++Ojz76CL/99hs6dOgglRcXFyM3NxeFhYU4ffo0Zs6cCXt7e3Tv3l3WZqtWrWThW9YmAJw4cQKenp44deoUiouL0bVrV1k9S0tLdO7cGT/99FOFfS8tLcWQIUNw8OBBxMTEoE2bNjh16hQWLVqE3377TW/exsGDB7Fp0ya89NJLsLe3x5IlSxAZGYlLly7B2dkZERER+O233/B///d/WLRoEZo0aQLgnz9cyuzZswfr169HXFwcmjRpIv2hcr+cnBz06NFDCvemTZti+/btiI6OhlarxdSpUwEAH3/8MaZMmYJhw4bh5Zdfxt27d3Hy5EkcOXIEo0ePlrXZpUsX/PDDDxW+J2QCtTyiJwW5/9Tq5s2bBQAxb948Wb1hw4YJlUolMjMzpbLyTm3evn1brywtLU0AEJ9//rlUVtVT30OHDhVWVlbi4sWLUtmZM2eEubm53mlhb29vWZ86deokBg0aVGH75Z1ezsrKEgCEg4ODuHbtmsF1ycnJUllUVJQAICZPniyVlZaWikGDBglLS0tx/fr1Co/bUJsVnfoGIGbNmiUtDx06VFhaWorz589LZVeuXBH29vaid+/eUllycrIAIEJDQ2XXHuPj44W5ubnIzc01uL8yzz77rHB2dq6wjhAVn/q2tbUV48eP1yvftm2bACBSU1Nl5WWfn7JX69at9d6/du3aib59++q1+csvvwgAYuXKlUIIITZs2CAAiAMHDujVHT58uHBzc6vwuNasWSPMzMz0rnuvXLlSABA//PCDVAZAWFpayv7d/PzzzwKAWLp0qVRW0alvAMLMzEz88ssvBtfd+xmIjo4W7u7u4q+//pLVGzlypNBoNNK/zaeeekq0a9euwuMsExMTI6ytratUl4yHp76pXN9++y3Mzc0xZcoUWfkrr7wCIQS2b99eaRvW1tbS/xcVFeHGjRto0aIFHB0dDc7qrUhJSQl27NiBoUOHwsvLSypv06YNwsLCKt3e0dERv/zyC86dO1et/d4rMjJSNrqpTFxcnPT/ZSObwsJC7Nq1q8Z9qExJSQm+++47DB06FM2bN5fK3d3dMXr0aBw8eBBarVa2TUxMjOw09eOPP46SkhJcvHixwn3duHEDjRs3fqD+3rlzB2q1Wq/cyspKWn+vtm3bYufOndi8eTOmTZsGW1tbvRnaVW2z7L/l1b1/3/fbsGED2rRpA39/f/z111/Sq2/fvgCAvXv3yuqHhoZKZ5sAoGPHjnBwcMDvv/9e4X7u1adPH7Rt27bCOkIIfPXVVxg8eDCEELK+hYWFIS8vT/r35+joiD/++ANHjx6tdN+NGzfGnTt3cPv27Sr3lx4cg5rKdfHiRXh4eMDe3l5W3qZNG2l9Ze7cuYO33npLusbdpEkTNG3aFLm5ucjLy6tWf65fv447d+6gZcuWeutat25d6fZz585Fbm4uWrVqhQ4dOuC1117DyZMnq9UHX1/fKtc1MzOTBSUAtGrVCgAe6Bp9Za5fv47bt28bfE/atGmD0tJS6RptmXv/8AEghe/ff/9d6f7EPdfca8La2trgdei7d+9K6+/l4OCA0NBQPPXUU3j33Xfxyiuv4KmnnsLPP/9c7TbL/lte3fv3fb9z587hl19+QdOmTWWvsp/ztWvXZPXvf5+Bf97rqrzPZaryGbx+/Tpyc3Px0Ucf6fVt3Lhxsr5Nnz4ddnZ26N69O1q2bInY2NhyT2+X/ayrMveAjIfXqMmkJk+ejOTkZEydOhVBQUHQaDRQqVQYOXJkudd2TaV37944f/48tmzZgu+++w6ffPIJFi1ahJUrV2LChAlVaqOyX9zVVd4vvJKSEqPupzLm5uYGyysLYWdn52qFjCHu7u64evWqXnlZmYeHR4XbR0RE4LnnnsO6devQqVMnqc0///yz0jbd3d1l5ffXrWzfpaWl6NChAxYuXGhwvaenp2y5pu/zvaryGSz7t/Xss88iKirKYJ2OHTsC+OePt4yMDGzduhWpqan46quvsHz5crz11luYM2eObJu///4bNjY2Rv93QBVjUFO5vL29sWvXLty6dUs2qj579qy0vkx5gbNx40ZERUXh3//+t1R29+7dGj3MoWnTprC2tjZ46jojI6NKbTg5OWHcuHEYN24c8vPz0bt3b8yePVsKamOOFEpLS/H7779LoysA+O233wD8b6Zy2cj1/vfD0NmKqvatadOmsLGxMfienD17FmZmZnoBUlP+/v5Yu3Yt8vLyoNFoatRG586d8f3336O0tFQ2oezIkSOwsbGRvX+G6HQ6lJaWys7QdO7cGXv37oVWq5VNKDty5Ii0HgDat28PCwsLHDt2DCNGjJDqFRYW4sSJE7IyQ/z8/PDzzz+jX79+RvvsGKOdpk2bwt7eHiUlJQZn5N/P1tYWzzzzDJ555hkUFhYiIiIC77zzDhITE6XLBcA/D7cpO6NGDw9PfVO5Bg4ciJKSEnzwwQey8kWLFkGlUiE8PFwqs7W1NRi+5ubmeqOFpUuX1mjEaG5ujrCwMGzevBmXLl2Syn/99Vfs2LGj0u1v3LghW7azs0OLFi1kpz1tbW0B6AdnTd373gkh8MEHH6BRo0bo168fgH/+2DE3N8eBAwdk2y1fvlyvrar2zdzcHP3798eWLVtkp9hzcnKQkpKCXr166c2GrqmgoCAIIZCenl7jNoYNG4acnBxs2rRJKvvrr7+wYcMGDB48WLp+nJubi6KiIr3tP/nkEwCQzdweNmwYSkpK8NFHH0llOp0OycnJCAwMlP5Q0Wg0CA0NxRdffIFbt25JddesWYP8/HwMHz68wr6PGDECf/75Jz7++GO9dXfu3EFBQUFV3gIZY3wGzc3NERkZia+++gqnT5/WW3/vrXf3/7uwtLRE27ZtIYTQe7+PHz+Oxx57rMb9oprhiJrKNXjwYISEhOCNN97AhQsX0KlTJ3z33XfYsmULpk6dKpsUExAQgF27dmHhwoXw8PCAr68vAgMD8eSTT2LNmjXQaDRo27Yt0tLSsGvXLjg7O9eoT3PmzEFqaioef/xxvPTSSyguLpbuAa3senPbtm0RHBws3Q967NgxbNy4UTbhKyAgAAAwZcoUhIWFwdzcHCNHjqxRX62srJCamoqoqCgEBgZi+/bt2LZtG15//XVpQppGo8Hw4cOxdOlSqFQq+Pn5YevWrXrXNqvbt3nz5mHnzp3o1asXXnrpJVhYWODDDz+ETqfDggULanQ8hvTq1QvOzs7YtWuXNIGqzDfffCNdNy4qKsLJkycxb948AMCQIUOkU6/Dhg1Djx49MG7cOJw5cwZNmjTB8uXLUVJSIjv1um/fPuk2opYtW6KwsBDff/89Nm3ahK5du+LZZ5+V6gYGBmL48OFITEzEtWvX0KJFC3z22We4cOECVq1aJevnO++8g8ceewx9+vRBTEwM/vjjD/z73/9G//79MWDAgAqP/7nnnsP69esxadIk7N27Fz179kRJSQnOnj2L9evXY8eOHXq3flWm7Of8xhtvYOTIkWjUqBEGDx4sBXhVzZ8/H3v37kVgYCAmTpyItm3b4ubNmzh+/Dh27dqFmzdvAgD69+8PNzc39OzZE66urvj111/xwQcfYNCgQbIzaenp6bh58yaeeuqpavWDjKCWZpuTAhm6/efWrVsiPj5eeHh4iEaNGomWLVuK9957T+8xgmfPnhW9e/cW1tbWAoB0K87ff/8txo0bJ5o0aSLs7OxEWFiYOHv2rN6tU9V5Mtn+/ftFQECAsLS0FM2bNxcrV640+GSy+/cxb9480b17d+Ho6Cisra2Fv7+/eOedd2SPGS0uLhaTJ08WTZs2FSqVSmqz7Hape590Vaa827NsbW3F+fPnRf/+/YWNjY1wdXUVs2bN0nvU5fXr10VkZKSwsbERjRs3Fi+88II4ffq0Xpvl9U0I/VtzhBDi+PHjIiwsTNjZ2QkbGxsREhIiDh06JKtTdnvW0aNHZeXV+XlMmTJFtGjRQq+87BY1Q697j0sIIW7evCmio6OFs7OzsLGxEX369NHrU2Zmpnj++edF8+bNhbW1tbCyshLt2rUTs2bNEvn5+Xr7v3Pnjnj11VeFm5ubUKvVolu3bnq3epX5/vvvxWOPPSasrKxE06ZNRWxsrNBqtZUeuxD/PJ3t3XffFe3atRNqtVo0btxYBAQEiDlz5siedgZAxMbG6m1//+dUCCHefvtt8cgjjwgzMzPZrVrltVG27v7PQE5OjoiNjRWenp6iUaNGws3NTfTr10989NFHUp0PP/xQ9O7dWzg7Owu1Wi38/PzEa6+9Juu7EEJMnz5deHl58RGitUAlxANO2SSiBu3333+Hv78/tm/fLp3Sp/pFp9PBx8cHM2bMwMsvv1zb3WlweI2aiB5I8+bNER0djfnz59d2V8hEkpOT0ahRI0yaNKm2u9IgcURNRESkYBxRExERKRiDmoiISMEY1ERERArGoCYiIlKwOvnAk9LSUly5cgX29vZ8ODwREdU5QgjcunULHh4eet/Ffr86GdRXrlwx2rOKiYiIasvly5fRrFmzCuvUyaAue6zd5cuXjfbMYiIioodFq9XC09NT72uEDamTQV12utvBwYFBTUREdVZVLt9yMhkREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUrE4+QlQJfGZsk/7/wvxBtdgTIiKqzziiJiIiUjAGNRERkYIxqImIiBSMQU1ERKRgDGoiIiIFY1ATEREpGIOaiIhIwRjURERECsagJiIiUjAGNRERkYIxqImIiBSMQU1ERKRg1QrqpKQkdOvWDfb29nBxccHQoUORkZEhq3P37l3ExsbC2dkZdnZ2iIyMRE5OjqzOpUuXMGjQINjY2MDFxQWvvfYaiouLH/xoTMhnxjbZi4iI6GGoVlDv378fsbGxOHz4MHbu3ImioiL0798fBQUFUp34+Hh888032LBhA/bv348rV64gIiJCWl9SUoJBgwahsLAQhw4dwmeffYbVq1fjrbfeMt5RERER1RMqIYSo6cbXr1+Hi4sL9u/fj969eyMvLw9NmzZFSkoKhg0bBgA4e/Ys2rRpg7S0NPTo0QPbt2/Hk08+iStXrsDV1RUAsHLlSkyfPh3Xr1+HpaVlpfvVarXQaDTIy8uDg4NDTbtfLRWNovk1l0REVB3VybEHukadl5cHAHBycgIApKeno6ioCKGhoVIdf39/eHl5IS0tDQCQlpaGDh06SCENAGFhYdBqtfjll18M7ken00Gr1cpeREREDUGNg7q0tBRTp05Fz5490b59ewBAdnY2LC0t4ejoKKvr6uqK7Oxsqc69IV22vmydIUlJSdBoNNLL09Ozpt0mIiKqUyxqumFsbCxOnz6NgwcPGrM/BiUmJiIhIUFa1mq1Jg9rThgjIiIlqFFQx8XFYevWrThw4ACaNWsmlbu5uaGwsBC5ubmyUXVOTg7c3NykOj/++KOsvbJZ4WV17qdWq6FWq2vSVSIiojqtWqe+hRCIi4vD119/jT179sDX11e2PiAgAI0aNcLu3bulsoyMDFy6dAlBQUEAgKCgIJw6dQrXrl2T6uzcuRMODg5o27btgxwLERFRvVOtEXVsbCxSUlKwZcsW2NvbS9eUNRoNrK2todFoEB0djYSEBDg5OcHBwQGTJ09GUFAQevToAQDo378/2rZti+eeew4LFixAdnY2Zs6cidjYWI6aiYiI7lOtoF6xYgUAIDg4WFaenJyMsWPHAgAWLVoEMzMzREZGQqfTISwsDMuXL5fqmpubY+vWrXjxxRcRFBQEW1tbREVFYe7cuQ92JERERPXQA91HXVsexn3U1ZlMxvuoiYioOh7afdRERERkWgxqIiIiBWNQExERKRiDmoiISMEY1ERERApW40eI0v/cP0Ocs8CJiMhYOKImIiJSMAY1ERGRgjGoiYiIFIxBTUREpGAMaiIiIgVjUBMRESkYb8+qg3g7GBFRw8GgVoDKgrc63+RFRET1C4O6HqgoyDnaJiKq23iNmoiISMEY1ERERArGU9//H68DExGREnFETUREpGDVHlEfOHAA7733HtLT03H16lV8/fXXGDp0qLRepVIZ3G7BggV47bXXAAA+Pj64ePGibH1SUhJmzJhR3e7UCRytExFRTVU7qAsKCtCpUyeMHz8eEREReuuvXr0qW96+fTuio6MRGRkpK587dy4mTpwoLdvb21e3K2QEvCebiEjZqh3U4eHhCA8PL3e9m5ubbHnLli0ICQlB8+bNZeX29vZ6dYmIiEjOpNeoc3JysG3bNkRHR+utmz9/PpydnfHoo4/ivffeQ3Fxcbnt6HQ6aLVa2YuIiKghMOms788++wz29vZ6p8inTJmCLl26wMnJCYcOHUJiYiKuXr2KhQsXGmwnKSkJc+bMMWVX6f/jqXAiImUxaVB/+umnGDNmDKysrGTlCQkJ0v937NgRlpaWeOGFF5CUlAS1Wq3XTmJiomwbrVYLT09P03WciIhIIUwW1N9//z0yMjLw5ZdfVlo3MDAQxcXFuHDhAlq3bq23Xq1WGwxwIiKi+s5k16hXrVqFgIAAdOrUqdK6J06cgJmZGVxcXEzVHSIiojqp2iPq/Px8ZGZmSstZWVk4ceIEnJyc4OXlBeCfU9MbNmzAv//9b73t09LScOTIEYSEhMDe3h5paWmIj4/Hs88+i8aNGz/AoRAREdU/1Q7qY8eOISQkRFouu3YcFRWF1atXAwDWrVsHIQRGjRqlt71arca6deswe/Zs6HQ6+Pr6Ij4+XnYNmoiIiP5R7aAODg6GEKLCOjExMYiJiTG4rkuXLjh8+HB1d9ug8ElmRERUhs/6JiIiUjB+exaZ1L1nB3hPNhFR9XFETUREpGAMaiIiIgXjqW8T4GQwIiIyFo6oiYiIFIwjaqo1/AIQIqLKcURNRESkYAxqIiIiBWNQExERKRivUZNRccZ7/cOH1hDVLo6oiYiIFIwjaqoWztQmInq4GNT1HIOViKhu46lvIiIiBWNQExERKRhPfdMD4SxvIiLT4oiaiIhIwRjUREREClbtoD5w4AAGDx4MDw8PqFQqbN68WbZ+7NixUKlUsteAAQNkdW7evIkxY8bAwcEBjo6OiI6ORn5+/gMdCBERUX1U7WvUBQUF6NSpE8aPH4+IiAiDdQYMGIDk5GRpWa1Wy9aPGTMGV69exc6dO1FUVIRx48YhJiYGKSkp1e0ONSC81YyIGqJqB3V4eDjCw8MrrKNWq+Hm5mZw3a+//orU1FQcPXoUXbt2BQAsXboUAwcOxPvvvw8PD4/qdomIiKjeMsms73379sHFxQWNGzdG3759MW/ePDg7OwMA0tLS4OjoKIU0AISGhsLMzAxHjhzB008/rdeeTqeDTqeTlrVarSm63SBwljYRUd1i9MlkAwYMwOeff47du3fj3Xffxf79+xEeHo6SkhIAQHZ2NlxcXGTbWFhYwMnJCdnZ2QbbTEpKgkajkV6enp7G7jYREZEiGX1EPXLkSOn/O3TogI4dO8LPzw/79u1Dv379atRmYmIiEhISpGWtVsuwfkg4Aiciql0mvz2refPmaNKkCTIzMwEAbm5uuHbtmqxOcXExbt68We51bbVaDQcHB9mLiIioITD5k8n++OMP3LhxA+7u7gCAoKAg5ObmIj09HQEBAQCAPXv2oLS0FIGBgabuDlGNcMY5EdWWagd1fn6+NDoGgKysLJw4cQJOTk5wcnLCnDlzEBkZCTc3N5w/fx7Tpk1DixYtEBYWBgBo06YNBgwYgIkTJ2LlypUoKipCXFwcRo4cyRnfRERE96l2UB87dgwhISHSctm146ioKKxYsQInT57EZ599htzcXHh4eKB///54++23ZfdSr127FnFxcejXrx/MzMwQGRmJJUuWGOFwSMl4vZuIqPqqHdTBwcEQQpS7fseOHZW24eTkxIebENVBvARA9PDxWd9EREQKxq+5JMXiqXIiIo6oiYiIFI1BTUREpGA89U31Fic+EVF9wKCmOotBTEQNAYOa6g1OPiOi+ohBTdTA8EwEUd3CyWREREQKxhE1NRgcSVYNLyEQKQtH1ERERArGoCYiIlIwnvomxeAp15rhKX2i+o1BTQSGHREpF099ExERKRhH1EQmwBE6ERkLR9REREQKxqAmIiJSMAY1ERGRgvEaNZERKPnWMiX3jYgqV+0R9YEDBzB48GB4eHhApVJh8+bN0rqioiJMnz4dHTp0gK2tLTw8PPD888/jypUrsjZ8fHygUqlkr/nz5z/wwRAREdU31R5RFxQUoFOnThg/fjwiIiJk627fvo3jx4/jzTffRKdOnfD333/j5ZdfxpAhQ3Ds2DFZ3blz52LixInSsr29fQ0PgYjuxRE0Uf1S7aAODw9HeHi4wXUajQY7d+6UlX3wwQfo3r07Ll26BC8vL6nc3t4ebm5uVdqnTqeDTqeTlrVabXW7TUREVCeZ/Bp1Xl4eVCoVHB0dZeXz58/H22+/DS8vL4wePRrx8fGwsDDcnaSkJMyZM8fUXSWScFRKREph0qC+e/cupk+fjlGjRsHBwUEqnzJlCrp06QInJyccOnQIiYmJuHr1KhYuXGiwncTERCQkJEjLWq0Wnp6epuw6ERGRIpgsqIuKijBixAgIIbBixQrZuntDt2PHjrC0tMQLL7yApKQkqNVqvbbUarXBcqL6gk8yI6LymOQ+6rKQvnjxInbu3CkbTRsSGBiI4uJiXLhwwRTdISIiqrOMPqIuC+lz585h7969cHZ2rnSbEydOwMzMDC4uLsbuDhFVgtfjiZSt2kGdn5+PzMxMaTkrKwsnTpyAk5MT3N3dMWzYMBw/fhxbt25FSUkJsrOzAQBOTk6wtLREWloajhw5gpCQENjb2yMtLQ3x8fF49tln0bhxY+MdGZGCMRyJqKqqHdTHjh1DSEiItFx2vTkqKgqzZ8/Gf//7XwBA586dZdvt3bsXwcHBUKvVWLduHWbPng2dTgdfX1/Ex8fLrlsTUd3Aa+tEplftoA4ODoYQotz1Fa0DgC5duuDw4cPV3S2R0XFUS0R1Ab+Ug4iISMEY1ERERArGoCYiIlIwBjUREZGCMaiJiIgUjEFNRESkYCb/9iwi4q1gRFRzHFETEREpGIOaiIhIwRjURERECsZr1ET00PDZ4ETVx6AmIqNhEBMZH099ExERKRhH1EQ18DBvt+IolahhY1AT1TG8J5uoYeGpbyIiIgVjUBMRESkYg5qIiEjBGNREREQKVu2gPnDgAAYPHgwPDw+oVCps3rxZtl4Igbfeegvu7u6wtrZGaGgozp07J6tz8+ZNjBkzBg4ODnB0dER0dDTy8/Mf6ECIiIjqo2oHdUFBATp16oRly5YZXL9gwQIsWbIEK1euxJEjR2Bra4uwsDDcvXtXqjNmzBj88ssv2LlzJ7Zu3YoDBw4gJiam5kdBRERUT1X79qzw8HCEh4cbXCeEwOLFizFz5kw89dRTAIDPP/8crq6u2Lx5M0aOHIlff/0VqampOHr0KLp27QoAWLp0KQYOHIj3338fHh4eeu3qdDrodDppWavVVrfbREREdZJRr1FnZWUhOzsboaGhUplGo0FgYCDS0tIAAGlpaXB0dJRCGgBCQ0NhZmaGI0eOGGw3KSkJGo1Genl6ehqz20RERIpl1AeeZGdnAwBcXV1l5a6urtK67OxsuLi4yDthYQEnJyepzv0SExORkJAgLWu1WoY11Wt8qAkRlakTTyZTq9VQq9W13Q0iIqKHzqinvt3c3AAAOTk5svKcnBxpnZubG65duyZbX1xcjJs3b0p1iIiI6B9GDWpfX1+4ublh9+7dUplWq8WRI0cQFBQEAAgKCkJubi7S09OlOnv27EFpaSkCAwON2R0iIqI6r9qnvvPz85GZmSktZ2Vl4cSJE3BycoKXlxemTp2KefPmoWXLlvD19cWbb74JDw8PDB06FADQpk0bDBgwABMnTsTKlStRVFSEuLg4jBw50uCMbyIiooas2kF97NgxhISESMtlk7yioqKwevVqTJs2DQUFBYiJiUFubi569eqF1NRUWFlZSdusXbsWcXFx6NevH8zMzBAZGYklS5YY4XCIiIjqF5UQQtR2J6pLq9VCo9EgLy8PDg4ORmmTs2yJjO/+787md2sT/aM6OcZnfRMRESkYg5qIiEjBGNREREQKxqAmIiJSMAY1ERGRgjGoiYiIFIxBTUREpGB14ks5iKh+4n3VRJVjUBORYjC4ifQxqInIZB70iX8MbiJeoyYiIlI0BjUREZGCMaiJiIgUjEFNRESkYAxqIiIiBWNQExERKRiDmoiISMF4HzUR1QuV3XPNe7KpruKImoiISMGMHtQ+Pj5QqVR6r9jYWABAcHCw3rpJkyYZuxtERET1gtFPfR89ehQlJSXS8unTp/HEE09g+PDhUtnEiRMxd+5cadnGxsbY3SAiIqoXjB7UTZs2lS3Pnz8ffn5+6NOnj1RmY2MDNze3Krep0+mg0+mkZa1W++AdJSIiqgNMeo26sLAQX3zxBcaPHw+VSiWVr127Fk2aNEH79u2RmJiI27dvV9hOUlISNBqN9PL09DRlt4mIqs1nxjbZi8hYTDrre/PmzcjNzcXYsWOlstGjR8Pb2xseHh44efIkpk+fjoyMDGzatKncdhITE5GQkCAta7VahjURETUIJg3qVatWITw8HB4eHlJZTEyM9P8dOnSAu7s7+vXrh/Pnz8PPz89gO2q1Gmq12pRdJSIiUiSTnfq+ePEidu3ahQkTJlRYLzAwEACQmZlpqq4QERHVWSYbUScnJ8PFxQWDBlX8UIETJ04AANzd3U3VFSJqgHidmOoLkwR1aWkpkpOTERUVBQuL/+3i/PnzSElJwcCBA+Hs7IyTJ08iPj4evXv3RseOHU3RFSIiojrNJEG9a9cuXLp0CePHj5eVW1paYteuXVi8eDEKCgrg6emJyMhIzJw50xTdICIiqvNMEtT9+/eHEEKv3NPTE/v37zfFLomIiOolPuubiIhIwRjURERECsagJiIiUjAGNRERkYKZ9MlkRET11cO8T/v+fV2YX/HzKah+4YiaiIhIwRjURERECsZT30RUZ/ExodQQMKiJqEGq7nVf/lFAtYVBTURkAIOZlIJBTUQE0wYzZ23Tg+BkMiIiIgVjUBMRESkYT30TUZ3B68bUEHFETUREpGAcURMRmQBH/2QsDGoiIoVhyNO9eOqbiIhIwTiiJiKqZRxBU0WMPqKePXs2VCqV7OXv7y+tv3v3LmJjY+Hs7Aw7OztERkYiJyfH2N0gIiKqF0xy6rtdu3a4evWq9Dp48KC0Lj4+Ht988w02bNiA/fv348qVK4iIiDBFN4iIiOo8k5z6trCwgJubm155Xl4eVq1ahZSUFPTt2xcAkJycjDZt2uDw4cPo0aOHwfZ0Oh10Op20rNVqTdFtIiIixTHJiPrcuXPw8PBA8+bNMWbMGFy6dAkAkJ6ejqKiIoSGhkp1/f394eXlhbS0tHLbS0pKgkajkV6enp6m6DYREZHiGD2oAwMDsXr1aqSmpmLFihXIysrC448/jlu3biE7OxuWlpZwdHSUbePq6ors7Oxy20xMTEReXp70unz5srG7TUREpEhGP/UdHh4u/X/Hjh0RGBgIb29vrF+/HtbW1jVqU61WQ61WG6uLREREdYbJ76N2dHREq1atkJmZCTc3NxQWFiI3N1dWJycnx+A1bSIioobO5PdR5+fn4/z583juuecQEBCARo0aYffu3YiMjAQAZGRk4NKlSwgKCjJ1V4iI6oXK7rvm913XL0YP6ldffRWDBw+Gt7c3rly5glmzZsHc3ByjRo2CRqNBdHQ0EhIS4OTkBAcHB0yePBlBQUHlzvgmIiJqyIwe1H/88QdGjRqFGzduoGnTpujVqxcOHz6Mpk2bAgAWLVoEMzMzREZGQqfTISwsDMuXLzd2N4iIiOoFowf1unXrKlxvZWWFZcuWYdmyZcbeNRERUb3DL+UgIiJSMH4pBxHRQ8Yv4aDq4IiaiIhIwRjURERECsagJiIiUjAGNRERkYJxMhkRUT13/+Q1PrmsbuGImoiISMEY1ERERArGU99ERA0YT4srH0fURERECsagJiIiUjAGNRERkYLxGjURUQNjymeN85q38TGoiYjqGX7pR/3CoCYiIkllIc8R8sPHa9REREQKxqAmIiJSMAY1ERGRghk9qJOSktCtWzfY29vDxcUFQ4cORUZGhqxOcHAwVCqV7DVp0iRjd4WIiBTGZ8Y22YsqZ/TJZPv370dsbCy6deuG4uJivP766+jfvz/OnDkDW1tbqd7EiRMxd+5cadnGxsbYXSEiolrGMH5wRg/q1NRU2fLq1avh4uKC9PR09O7dWyq3sbGBm5ubsXdPRERUr5j8GnVeXh4AwMnJSVa+du1aNGnSBO3bt0diYiJu375dbhs6nQ5arVb2IiIiaghMeh91aWkppk6dip49e6J9+/ZS+ejRo+Ht7Q0PDw+cPHkS06dPR0ZGBjZt2mSwnaSkJMyZM8eUXSUiIlIkkwZ1bGwsTp8+jYMHD8rKY2JipP/v0KED3N3d0a9fP5w/fx5+fn567SQmJiIhIUFa1mq18PT0NF3HiYioVtx7TZsPV/mHyYI6Li4OW7duxYEDB9CsWbMK6wYGBgIAMjMzDQa1Wq2GWq02ST+JiIiUzOhBLYTA5MmT8fXXX2Pfvn3w9fWtdJsTJ04AANz/X3v3HlNl/ccB/H1SDogpFMYtAS/LSwZYlmfYxaVMvCyxnLe5FDMrw+YlN3Op6Gqzsjmmc+CaYM1pyaay0mmgcipFMSFFUkLGMCcHygYiyCXO5/dH4/w8njvCc77A+7WxeZ7zfb5+vs/n+T4fn8Pj+YaFdXY4RETUifgUt/Y6vVAnJydj//79yM7OxoABA2AymQAAAQEB6NevH8rLy7F//35Mnz4dQUFBuHz5MlavXo1XXnkFMTExnR0OERFRt9bphTotLQ3Af19qcr/MzEwkJSVBr9cjNzcXqampaGhoQEREBGbPno0NGzZ0dihERETdXpd89O1MREQEjEZjZ/+1REREPRK/65uIiEhhXI+aiIiU9OCDa731v2vxjpqIiEhhLNREREQKY6EmIiJSGAs1ERGRwvgwGRERdQu99eEy3lETEREpjHfURETUI3j6PeTd5Y6chZqIiLzmYRb56C0LhPCjbyIiIoWxUBMRESmMhZqIiEhhLNREREQKY6EmIiJSGAs1ERGRwlioiYiIFMZCTUREpDB+4QkRERHU/S5xr91R79q1C0OGDIGfnx8MBgMKCgq8FQoREZGyvHJH/d1332HNmjVIT0+HwWBAamoqEhISUFpaiuDgYG+EREREvUx3+QpSr9xRb9++HcuWLcOSJUvw9NNPIz09Hf7+/sjIyPBGOERERMrS/I66paUFFy9exPr16y3bHnnkEcTHxyM/P9/uPs3NzWhubra8rqurAwDcuXOn0+IyNzd2Wl9ERNT9dWaNcdS3iLhsq3mh/vvvv9HW1oaQkBCr7SEhIbh27ZrdfbZu3YotW7bYbI+IiOiSGImIiAJSu/7vqK+vR0BAgNM23eKp7/Xr12PNmjWW12azGf/88w+CgoKg0+keuv87d+4gIiICf/75JwYOHPjQ/XkDx6COnjAOjkEdPWEcHIMtEUF9fT3Cw8NdttW8UA8aNAh9+vRBdXW11fbq6mqEhoba3cfX1xe+vr5W2wIDAzs9toEDB3bbk6gdx6COnjAOjkEdPWEcHIM1V3fS7TR/mEyv12PcuHE4efKkZZvZbMbJkycRFxendThERERK88pH32vWrMHixYvx/PPPY/z48UhNTUVDQwOWLFnijXCIiIiU5ZVCPW/ePPz111/YtGkTTCYTxo4di+PHj9s8YKYVX19fpKSk2Hy83p1wDOroCePgGNTRE8bBMTwcnbjzbDgRERF5BRflICIiUhgLNRERkcJYqImIiBTGQk1ERKSwXlOoPV1WMysrC6NGjYKfnx+io6Nx7NgxjSK1tXXrVrzwwgsYMGAAgoODMWvWLJSWljrdZ+/evdDpdFY/fn5+GkVs3+bNm21iGjVqlNN9VMoDAAwZMsRmDDqdDsnJyXbbq5CHn376Ca+99hrCw8Oh0+lw5MgRq/dFBJs2bUJYWBj69euH+Ph4lJWVuexX66VqnY2jtbUV69atQ3R0NPr374/w8HAsWrQIt27dctpnR87JrhoDACQlJdnEM3XqVJf9apkLV2OwNz90Oh22bdvmsE+t8+DONbWpqQnJyckICgrCo48+itmzZ9t8UdeDOjqXXOkVhbp9Wc2UlBQUFhYiNjYWCQkJqKmpsdv+7NmzWLBgAZYuXYqioiLMmjULs2bNwpUrVzSO/D9GoxHJyck4d+4ccnJy0NraiilTpqChocHpfgMHDkRVVZXlp7KyUqOIHRszZoxVTL/88ovDtqrlAQAuXLhgFX9OTg4AYM6cOQ738XYeGhoaEBsbi127dtl9/4svvsCOHTuQnp6O8+fPo3///khISEBTU5PDPj2dU53B2TgaGxtRWFiIjRs3orCwEIcOHUJpaSlmzpzpsl9PzsmH5SoXADB16lSreA4cOOC0T61z4WoM98deVVWFjIwM6HQ6zJ4922m/WubBnWvq6tWr8f333yMrKwtGoxG3bt3CG2+84bTfjswlt0gvMH78eElOTra8bmtrk/DwcNm6davd9nPnzpUZM2ZYbTMYDPLuu+92aZzuqqmpEQBiNBodtsnMzJSAgADtgnJDSkqKxMbGut1e9TyIiKxcuVKGDx8uZrPZ7vuq5QGAHD582PLabDZLaGiobNu2zbKttrZWfH195cCBAw778XROdbYHx2FPQUGBAJDKykqHbTw9JzuTvTEsXrxYEhMTPerHm7lwJw+JiYkyadIkp228mQcR22tqbW2t+Pj4SFZWlqXN1atXBYDk5+fb7aOjc8kdPf6Oun1Zzfj4eMs2V8tq5ufnW7UHgISEBIfttda+zOfjjz/utN3du3cRFRWFiIgIJCYmoqSkRIvwnCorK0N4eDiGDRuGhQsX4saNGw7bqp6HlpYW7Nu3D2+99ZbTxWFUzEO7iooKmEwmq+McEBAAg8Hg8Dh3ZE55Q11dHXQ6nct1ATw5J7WQl5eH4OBgjBw5EsuXL8ft27cdtlU9F9XV1Th69CiWLl3qsq038/DgNfXixYtobW21Oq6jRo1CZGSkw+Pakbnkrh5fqJ0tq2kymezuYzKZPGqvJbPZjFWrVuHFF1/EM88847DdyJEjkZGRgezsbOzbtw9msxkTJkzAzZs3NYzWmsFgwN69e3H8+HGkpaWhoqICL7/8Murr6+22VzkPAHDkyBHU1tYiKSnJYRsV83C/9mPpyXHuyJzSWlNTE9atW4cFCxY4XUDB03Oyq02dOhXffPMNTp48ic8//xxGoxHTpk1DW1ub3faq5+Lrr7/GgAEDXH5k7M082Lummkwm6PV6m3/kuaob7W3c3cdd3WKZS/q/5ORkXLlyxeXvb+Li4qwWOZkwYQJGjx6N3bt345NPPunqMO2aNm2a5c8xMTEwGAyIiorCwYMH3foXt2r27NmDadOmOV2mTsU89HStra2YO3cuRARpaWlO26p2Ts6fP9/y5+joaMTExGD48OHIy8vD5MmTNY/nYWVkZGDhwoUuH6D0Zh7cvaZ6U4+/o+7IspqhoaEetdfKihUr8MMPP+D06dMYPHiwR/v6+Pjg2WefxfXr17soOs8FBgZixIgRDmNSNQ8AUFlZidzcXLz99tse7adaHtqPpSfHuSNzSivtRbqyshI5OTkeL0fo6pzU2rBhwzBo0CCH8aici59//hmlpaUezxFAuzw4uqaGhoaipaUFtbW1Vu1d1Y32Nu7u464eX6g7sqxmXFycVXsAyMnJ8doynCKCFStW4PDhwzh16hSGDh3qcR9tbW0oLi5GWFhYF0TYMXfv3kV5ebnDmFTLw/0yMzMRHByMGTNmeLSfankYOnQoQkNDrY7znTt3cP78eYfHWdWlatuLdFlZGXJzcxEUFORxH67OSa3dvHkTt2/fdhiPqrkA/vvEady4cYiNjfV4367Og6tr6rhx4+Dj42N1XEtLS3Hjxg2Hx7Ujc8mTgHu8b7/9Vnx9fWXv3r3y+++/yzvvvCOBgYFiMplEROTNN9+Ujz76yNL+zJkz0rdvX/nyyy/l6tWrkpKSIj4+PlJcXOyV+JcvXy4BAQGSl5cnVVVVlp/GxkZLmwfHsGXLFjlx4oSUl5fLxYsXZf78+eLn5yclJSXeGIKIiHz44YeSl5cnFRUVcubMGYmPj5dBgwZJTU2NiKifh3ZtbW0SGRkp69ats3lPxTzU19dLUVGRFBUVCQDZvn27FBUVWZ6G/uyzzyQwMFCys7Pl8uXLkpiYKEOHDpV79+5Z+pg0aZLs3LnT8trVnNJ6HC0tLTJz5kwZPHiw/Pbbb1bzpLm52eE4XJ2TWo6hvr5e1q5dK/n5+VJRUSG5ubny3HPPyVNPPSVNTU0Ox6B1LlydTyIidXV14u/vL2lpaXb78HYe3LmmvvfeexIZGSmnTp2SX3/9VeLi4iQuLs6qn5EjR8qhQ4csr92ZSx3RKwq1iMjOnTslMjJS9Hq9jB8/Xs6dO2d5b+LEibJ48WKr9gcPHpQRI0aIXq+XMWPGyNGjRzWO+P8A2P3JzMy0tHlwDKtWrbKMNyQkRKZPny6FhYXaB3+fefPmSVhYmOj1ennyySdl3rx5cv36dcv7queh3YkTJwSAlJaW2rynYh5Onz5t9/xpj9NsNsvGjRslJCREfH19ZfLkyTZji4qKkpSUFKttzuaU1uOoqKhwOE9Onz7tcByuzkktx9DY2ChTpkyRJ554Qnx8fCQqKkqWLVtmU3C9nQtX55OIyO7du6Vfv35SW1trtw9v58Gda+q9e/fk/fffl8cee0z8/f3l9ddfl6qqKpt+7t/HnbnUEVzmkoiISGE9/nfURERE3RkLNRERkcJYqImIiBTGQk1ERKQwFmoiIiKFsVATEREpjIWaiIhIYSzURERECmOhJiIiUhgLNREB+G/BkAkTJtisHVxXV4eIiAh8/PHHXoqMqHfjV4gSkcUff/yBsWPH4quvvsLChQsBAIsWLcKlS5dw4cIF6PV6L0dI1PuwUBORlR07dmDz5s0oKSlBQUEB5syZgwsXLnRouUIiengs1ERkRUQwadIk9OnTB8XFxfjggw+wYcMGb4dF1GuxUBORjWvXrmH06NGIjo5GYWEh+vbt6+2QiHotPkxGRDYyMjLg7++PiooK3Lx509vhEPVqvKMmIitnz57FxIkT8eOPP+LTTz8FAOTm5kKn03k5MqLeiXfURGTR2NiIpKQkLF++HK+++ir27NmDgoICpKenezs0ol6Ld9REZLFy5UocO3YMly5dgr+/PwBg9+7dWLt2LYqLizFkyBDvBkjUC7FQExEAwGg0YvLkycjLy8NLL71k9V5CQgL+/fdffgRO5AUs1ERERArj76iJiIgUxkJNRESkMBZqIiIihbFQExERKYyFmoiISGEs1ERERApjoSYiIlIYCzUREZHCWKiJiIgUxkJNRESkMBZqIiIihf0PNuJt3hXFoScAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 500x1000 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "singleToy(nevents_sig=300, nevents_bgd=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8098feff",
   "metadata": {},
   "source": [
    "## Problem 2 (1 point):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bcb70b40",
   "metadata": {},
   "source": [
    "The p-value corresponds to the probability of obtaining test results that are at least as extreme as the results that have been observed. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "fc20bda2",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = 10000\n",
    "sigma = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c659f074",
   "metadata": {},
   "outputs": [],
   "source": [
    "area_of_observation = Gaussian(mean,sigma).integrate(10300,(10000+(10*sigma))) # this is normalised to is technically \n",
    "                                                                               # already the p-value\n",
    "total_area = Gaussian(mean,sigma).integrate((10000-(10*sigma)),(10000+(10*sigma)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "6040e37b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p-value is:  0.0013498980316301115\n"
     ]
    }
   ],
   "source": [
    "pvalue = area_of_observation/total_area\n",
    "print(\"The p-value is: \",pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "44d1f6fe",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Z-score is:  3.2051549205989076\n"
     ]
    }
   ],
   "source": [
    "# Calcualting the Z-score\n",
    "n_sigmas = erfinv(1-pvalue)*np.sqrt(2)\n",
    "print(\"The Z-score is: \",n_sigmas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f9cc095",
   "metadata": {},
   "source": [
    "## Problem 3 (1 point):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bf6559d7",
   "metadata": {},
   "source": [
    "Confirming the total number of background events between 5<x<15 is 50% of the total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "3d48da6b",
   "metadata": {},
   "outputs": [],
   "source": [
    "pdf = SignalWithBackground(mean, sigma, 0, intercept, slope, XMIN, XMAX)\n",
    "for i in range (10000): pdf.next()\n",
    "bgd_data = np.array(pdf.mass_bgd)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "53fed5b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "number_in_subregion = np.where((bgd_data<15) & (bgd_data>5))[0].shape[0]\n",
    "total_number_background = bgd_data.size\n",
    "fraction = number_in_subregion/total_number_background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "ebca87ac",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fraction of background events in the sub-region 5<x<15 is:  0.5017\n"
     ]
    }
   ],
   "source": [
    "print(\"Fraction of background events in the sub-region 5<x<15 is: \",fraction)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5eb5a615",
   "metadata": {},
   "source": [
    "We expect this fraction to be 50% as the distribution is linear between 0 and 20, but the spread fluctuates and the fraction is slightly different at each generation of the pdf, with it averaging to 50% over many initialisations of the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "ea4539c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = 5000\n",
    "sigma = np.sqrt(5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "1c7ca7b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "area_of_observation = Gaussian(mean,sigma).integrate(5300,(5000+(10*sigma))) # this is normalised to is technically \n",
    "                                                                               # already the p-value\n",
    "total_area = Gaussian(mean,sigma).integrate((5000-(10*sigma)),(5000+(10*sigma)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "de50f431",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p-value is:  1.1045248499292971e-05\n"
     ]
    }
   ],
   "source": [
    "pvalue = area_of_observation/total_area\n",
    "print(\"The p-value is: \",pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "0011fc32",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Z-score is:  4.3956289210855894\n"
     ]
    }
   ],
   "source": [
    "# Calcualting the Z-score\n",
    "n_sigmas = erfinv(1-pvalue)*np.sqrt(2)\n",
    "print(\"The Z-score is: \",n_sigmas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "43f99860",
   "metadata": {},
   "source": [
    "Why did the number of standard deviations (Z-score) increase after we filtered the background?:"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "69d81d1d",
   "metadata": {},
   "source": [
    "The deviation is a lot more significant now as the fraction of events that are observed above the expected number is now a larger portion of the expected number of events. This is because we filtered down the number of background events while the number of signal events was kept the same.\n",
    "\n",
    "This results in the number of observed events to deviate more extremely from the expected number of background events (in our case the deviation is more than 4 standard deviations away from the average results that would be expected if there was no signal). This significantly increases the likelihood of the alternate hypothesis of a signal being present which does not conform to the background distribution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e247b39e",
   "metadata": {},
   "source": [
    "## Problem 4 (5 points):"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55dbba25",
   "metadata": {},
   "source": [
    "### Resetting the parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "b3afc826",
   "metadata": {},
   "outputs": [],
   "source": [
    "# background parameters\n",
    "XMIN = 0.\n",
    "XMAX = 20.\n",
    "intercept = 20.\n",
    "slope = -1.\n",
    "# signal parameters\n",
    "mean = 10.\n",
    "sigma = 0.5\n",
    "# chose number of bins that is appropriate for the size of the statistics sample\n",
    "NBINS = 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b341627a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the new signal count\n",
    "nevents_sig = 150\n",
    "nevents_bgd = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bbe27d8f",
   "metadata": {},
   "source": [
    "### Plotting the new distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "06873d2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAPdCAYAAABFjSsGAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/av/WaAAAACXBIWXMAAA9hAAAPYQGoP6dpAACPkUlEQVR4nOzdeVxVdeL/8fcF5QKyCYpA4m7uS1Ea5YJpKppmYqbZhEtqhZoyk8ZM5VIzmDZli0vOlLboVLZoY4m5W4mWNmarqbmVgqYBinplOb8/+nG+Xtn1Xjng6/l4nIeez/mcz/mcey+8+Zzt2gzDMAQAACzJo6I7AAAAikdQAwBgYQQ1AAAWRlADAGBhBDUAABZGUAMAYGEENQAAFkZQAwBgYQQ1AAAWRlDjkjVo0EDDhw+v6G5Ikg4cOCCbzabFixdf0vqLFy+WzWbTgQMHzLKYmBjFxMS4pH+lsdlsmjZtmjk/bdo02Ww2/fbbb1dk+xX9Xs6aNUvNmzdXfn5+hfWhMnHn+5WSkiI/Pz8dP37cLe2j/AhqFPLNN99o0KBBql+/vry9vXXNNdfotttu04svvljRXbO8LVu2aNq0acrIyKjorhRi1b5lZWXp6aef1pQpU+Th8X+/kt5++23de++9atq0qWw2W5n/aPr73/8um82m1q1bF1r2ySefaNSoUWrdurU8PT3VoEEDF+3Fpfn+++81bdo0pz8QK1rv3r3VpEkTJScnV3RX8P9Vq+gOwFq2bNmibt26qV69eho9erTCwsJ0+PBhbd26Vc8//7zGjx9v1t29e7fTL9aq5pNPPin3Olu2bNH06dM1fPhwBQUFlXm9s2fPqlo19/44ltS3inwvX331VeXm5mro0KFO5fPnz9eOHTt044036sSJE2Vq65dfftE//vEP1ahRo8jlS5cu1dtvv63rr79eERERl933y/X9999r+vTpiomJKdcfDe5+v8aOHau//OUvmj59uvz9/d22HZQNQQ0nf//73xUYGKgvv/yy0C/zY8eOOc3b7fYr2LMrz8vLy63t5+fn6/z58/L29pa3t7dbt1WainwvFy1apP79+xd6Dd544w1dc8018vDwKHJ0XJS//OUvuummm5SXl1fkaYN//OMf+te//qXq1avr9ttv17fffuuSfbgSDMPQuXPn5OPj4/b3Ky4uTuPHj9eyZcs0cuRIt24Lpau6wyFckn379qlVq1ZFjgZDQ0Od5os6T7Zr1y517dpVPj4+qlu3rp566iktWrSo0PnfBg0a6Pbbb9dnn32mDh06yNvbW40aNdLrr7/u1N7Jkyf1l7/8RW3atJGfn58CAgIUGxurr7/++pL38bvvvtOtt97q1Meizo0WdY76xRdfVKtWreTr66uaNWvqhhtu0NKlSyX9cV75kUcekSQ1bNhQNpvNab9tNpvGjRunJUuWqFWrVrLb7UpJSTGXXXiOusBvv/2mwYMHKyAgQCEhIXr44Yd17tw5c3lJ5+YvbLO0vhX1Xv7888+66667FBwcLF9fX91000366KOPnOps3LhRNptN77zzjv7+97+rbt268vb2Vvfu3bV3795CfbrY/v37tWvXLvXo0aPQssjIyHKNGjdv3qx3331Xc+bMKbZORESEqlevXuY2L5afn685c+aoVatW8vb2Vp06dTR27Fj9/vvvTvXK8vlevHix7rrrLklSt27dzPdk48aNTm2sXr1aN9xwg3x8fPTyyy+byy5+vzIyMjRx4kRFRkbKbrerSZMmevrppwt9tt966y1FRUXJ399fAQEBatOmjZ5//nmnOqGhoWrbtq1WrFhxya8VXIcRNZzUr19fqamp+vbbb8s8iinw66+/mr9wkpKSVKNGDf373/8u9q//vXv3atCgQRo1apTi4+P16quvavjw4YqKilKrVq0k/REWy5cv11133aWGDRsqPT1dL7/8srp27arvv/++3Icv09LS1K1bN+Xm5urRRx9VjRo1tHDhQvn4+JS67r/+9S9NmDBBgwYNMgNz165d2rZtm+655x4NHDhQP/30k/7zn//oueeeU61atSRJtWvXNttYv3693nnnHY0bN061atUq9XDn4MGD1aBBAyUnJ2vr1q164YUX9Pvvvxf6g6Y0ZenbhdLT03XzzTfrzJkzmjBhgkJCQvTaa6+pf//+evfdd3XnnXc61Z85c6Y8PDz0l7/8RZmZmZo1a5aGDRumbdu2ldivLVu2SJKuv/76cu3PxfLy8jR+/Hjdf//9atOmzWW1VZKxY8dq8eLFGjFihCZMmKD9+/frpZde0v/+9z99/vnnTn8ElPb57tKliyZMmKAXXnhBf/3rX9WiRQtJMv+V/jjEPXToUI0dO1ajR49Ws2bNiuzXmTNn1LVrV/36668aO3as6tWrpy1btigpKUlHjx41/3hZs2aNhg4dqu7du+vpp5+WJP3www/6/PPP9fDDDzu1GRUVpeXLl7vw1cMlM4ALfPLJJ4anp6fh6elpREdHG5MnTzZWr15tnD9/vlDd+vXrG/Hx8eb8+PHjDZvNZvzvf/8zy06cOGEEBwcbkoz9+/c7rSvJ2Lx5s1l27Ngxw263G3/+85/NsnPnzhl5eXlO292/f79ht9uNGTNmOJVJMhYtWlTi/k2cONGQZGzbts1pu4GBgYX62LVrV6Nr167m/B133GG0atWqxPZnz55dqJ0CkgwPDw/ju+++K3LZ1KlTzfmpU6cakoz+/fs71XvooYcMScbXX39tGEbJ+31xmyX17eL3suB1+vTTT82yU6dOGQ0bNjQaNGhgvicbNmwwJBktWrQwHA6HWff55583JBnffPNNoW1d6LHHHjMkGadOnSqxXqtWrZzei4u99NJLRmBgoHHs2DHDMP5470p7r/r27WvUr1+/xDoX+vTTTw1JxpIlS5zKU1JSCpWX9fO9bNkyQ5KxYcOGQtsraCMlJaXIZRe+X08++aRRo0YN46effnKq9+ijjxqenp7GoUOHDMMwjIcfftgICAgwcnNzS93ff/zjH4YkIz09vdS6cC8OfcPJbbfdptTUVPXv319ff/21Zs2apV69eumaa67Rhx9+WOK6KSkpio6OVvv27c2y4OBgDRs2rMj6LVu2VOfOnc352rVrq1mzZvr555/NMrvdbh7+zMvL04kTJ+Tn56dmzZrpq6++Kvf+ffzxx7rpppvUoUMHp+0W18cLBQUF6ZdfftGXX35Z7u0W6Nq1q1q2bFnm+gkJCU7zBRfzffzxx5fch7L4+OOP1aFDB3Xq1Mks8/Pz05gxY3TgwAF9//33TvVHjBjhdE6/4H298L0syokTJ1StWjX5+fldcl9PnDihJ554Qo8//nixRwhcYdmyZQoMDNRtt92m3377zZyioqLk5+enDRs2ONUvy+e7NA0bNlSvXr3K1LfOnTurZs2aTn3r0aOH8vLytHnzZkl/fIazs7O1Zs2aUtusWbOmJF2xWwRRPIIahdx44416//339fvvv+uLL75QUlKSTp06pUGDBhX6BX2hgwcPqkmTJoXKiyqTpHr16hUqq1mzptP5vvz8fD333HNq2rSp7Ha7atWqpdq1a2vXrl3KzMws974dPHhQTZs2LVRe3CHFC02ZMkV+fn7q0KGDmjZtqoSEBH3++efl2n7Dhg3LVf/ivjZu3FgeHh5uv53n4MGDRb4mBYdlDx486FR+8XtZ8Ev+4nO37vDYY48pODjY6Y4Ed9izZ48yMzMVGhqq2rVrO02nT58udLFlWT7fpSnr52XPnj1KSUkp1K+Cc/8FfXvooYd07bXXKjY2VnXr1tXIkSPN6yQuZhiGpD+udUDF4hw1iuXl5aUbb7xRN954o6699lqNGDFCy5Yt09SpU13SvqenZ5HlBb8gpD+u0n388cc1cuRIPfnkkwoODpaHh4cmTpx4xR+O0aJFC+3evVsrV65USkqK3nvvPc2bN09PPPGEpk+fXqY2ynIuvCQX/9Is7pdoXl7eZW2nvMryXhYlJCREubm5OnXq1CXdBrRnzx4tXLhQc+bM0ZEjR8zyc+fOKScnRwcOHFBAQICCg4PL3fbF8vPzFRoaqiVLlhS5/OLR/KW+Jhcq6+clPz9ft912myZPnlzk8muvvVbSHxeJ7dy5U6tXr9aqVau0atUqLVq0SPfdd59ee+01p3UK/qAouJ4BFYegRpnccMMNkqSjR48WW6d+/fpFXulblqt/i/Puu++qW7dueuWVV5zKMzIyLukXSP369bVnz55C5bt37y7T+jVq1NDdd9+tu+++W+fPn9fAgQP197//XUlJSfL29nb56GPPnj1Oo6q9e/cqPz/fvAitYOR68UNMLh7xSuUbGdWvX7/I1+THH380l7tC8+bNJf1x9Xfbtm3Lvf6vv/6q/Px8TZgwQRMmTCi0vGHDhnr44YdLvBK8rBo3bqy1a9fqlltuuew/uAq46vPSuHFjnT59usir5y/m5eWlfv36qV+/fsrPz9dDDz2kl19+WY8//rjT0a/9+/ebR7BQsTj0DScbNmwo8i/+gnOiJR0i7tWrl1JTU7Vz506z7OTJk8WOQMrC09OzUH+WLVumX3/99ZLa69Onj7Zu3aovvvjCLDt+/HiZ+njxQze8vLzUsmVLGYahnJwcSTIftOGqp3/NnTvXab7g6XCxsbGSpICAANWqVcs8B1lg3rx5hdoqT9/69OmjL774QqmpqWZZdna2Fi5cqAYNGpTrPHtJoqOjJUnbt2+/pPVbt26tDz74oNDUqlUr1atXTx988IFGjRrlkr4OHjxYeXl5evLJJwsty83NvaT33FWfl8GDBys1NVWrV68utCwjI0O5ubmSCn+GPTw8zD+QHA6H07IdO3aY7w8qFiNqOBk/frzOnDmjO++8U82bN9f58+e1ZcsWvf3222rQoIFGjBhR7LqTJ0/Wm2++qdtuu03jx483b8+qV6+eTp48eUmjh9tvv10zZszQiBEjdPPNN+ubb77RkiVL1KhRo0vav8mTJ+uNN95Q79699fDDD5u3Z9WvX1+7du0qcd2ePXsqLCxMt9xyi+rUqaMffvhBL730kvr27Wseto2KipIk/e1vf9OQIUNUvXp19evXr9gnZZVm//796t+/v3r37q3U1FS9+eabuueee9SuXTuzzv3336+ZM2fq/vvv1w033KDNmzfrp59+KtRWefr26KOP6j//+Y9iY2M1YcIEBQcH67XXXtP+/fv13nvvueypWI0aNVLr1q21du3aQg/W2Lx5s/kHyPHjx5Wdna2nnnpKktSlSxd16dJFtWrV0oABAwq1WzCCvnjZrl27zIsi9+7dq8zMTLPNdu3aqV+/fsX2tWvXrho7dqySk5O1c+dO9ezZU9WrV9eePXu0bNkyPf/88xo0aFC59r99+/by9PTU008/rczMTNntdt16662FnllQmkceeUQffvihbr/9dvMWsOzsbH3zzTd69913deDAAdWqVUv333+/Tp48qVtvvVV169bVwYMH9eKLL6p9+/ZOt4UdO3ZMu3btKnQxIypIBV5xDgtatWqVMXLkSKN58+aGn5+f4eXlZTRp0sQYP358ods0Lr5FxDAM43//+5/RuXNnw263G3Xr1jWSk5ONF154wZBkpKWlOa3bt2/fQtu/+Jaoc+fOGX/+85+N8PBww8fHx7jllluM1NTUQvXKenuWYRjGrl27jK5duxre3t7GNddcYzz55JPGK6+8UurtWS+//LLRpUsXIyQkxLDb7Ubjxo2NRx55xMjMzHRq/8knnzSuueYaw8PDw6lNSUZCQkKRfVIxt2d9//33xqBBgwx/f3+jZs2axrhx44yzZ886rXvmzBlj1KhRRmBgoOHv728MHjzYOHbsWKE2S+pbUe/lvn37jEGDBhlBQUGGt7e30aFDB2PlypVOdQpuz1q2bJlTeXnej2effdbw8/Mzzpw541Re8BoUNV28Xxcr7vasRYsWFdvmxftfnIULFxpRUVGGj4+P4e/vb7Rp08aYPHmyceTIEbNOWT/fhmEY//rXv4xGjRoZnp6eTrdqFddGwbKL+3vq1CkjKSnJaNKkieHl5WXUqlXLuPnmm41nnnnGvL3y3XffNXr27GmEhoYaXl5eRr169YyxY8caR48edWpr/vz5hq+vr5GVlVWm1wTuZTOMclzZAFyCiRMn6uWXX9bp06eLvcAGV6/MzEw1atRIs2bNctlhalye6667TjExMXruuecquisQ56jhYmfPnnWaP3HihN544w116tSJkEaRAgMDNXnyZM2ePZuvubSAlJQU7dmzR0lJSRXdFfx/jKjhUu3bt1dMTIxatGih9PR0vfLKKzpy5IjWrVunLl26VHT3AKDS4WIyuFSfPn307rvvauHChbLZbLr++uv1yiuvENIAcIkYUQMAYGGcowYAwMIq5aHv/Px8HTlyRP7+/jyHFgBQ6RiGoVOnTikiIqLU5xJUyqA+cuSIIiMjK7obAABclsOHD6tu3bol1qmUQV3wFKjDhw8rICCggnsDAED5ZGVlKTIyskxfRlMpg7rgcHdAQABBDQCotMpy+paLyQAAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALq5SPEAXgfg0e/chp/sDMvhXUE+DqxogaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALc3lQz58/X23btlVAQIACAgIUHR2tVatWmcvPnTunhIQEhYSEyM/PT3FxcUpPT3d1NwAAqBJcHtR169bVzJkztWPHDm3fvl233nqr7rjjDn333XeSpEmTJum///2vli1bpk2bNunIkSMaOHCgq7sBAECVYDMMw3D3RoKDgzV79mwNGjRItWvX1tKlSzVo0CBJ0o8//qgWLVooNTVVN910U5nay8rKUmBgoDIzMxUQEODOrgNXrQaPfuQ0f2Bm3wrqCVD1lCfH3HqOOi8vT2+99Zays7MVHR2tHTt2KCcnRz169DDrNG/eXPXq1VNqamqx7TgcDmVlZTlNAABcDdwS1N988438/Pxkt9v1wAMP6IMPPlDLli2VlpYmLy8vBQUFOdWvU6eO0tLSim0vOTlZgYGB5hQZGemObgMAYDluCepmzZpp586d2rZtmx588EHFx8fr+++/v+T2kpKSlJmZaU6HDx92YW8BALCuau5o1MvLS02aNJEkRUVF6csvv9Tzzz+vu+++W+fPn1dGRobTqDo9PV1hYWHFtme322W3293RVQAALO2K3Eedn58vh8OhqKgoVa9eXevWrTOX7d69W4cOHVJ0dPSV6AoAAJWKy0fUSUlJio2NVb169XTq1CktXbpUGzdu1OrVqxUYGKhRo0YpMTFRwcHBCggI0Pjx4xUdHV3mK74BALiauDyojx07pvvuu09Hjx5VYGCg2rZtq9WrV+u2226TJD333HPy8PBQXFycHA6HevXqpXnz5rm6GwAAVAlX5D5qV+M+asD9uI8acB/L3EcNAAAuD0ENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIW5PKiTk5N14403yt/fX6GhoRowYIB2797tVCcmJkY2m81peuCBB1zdFQAAKj2XB/WmTZuUkJCgrVu3as2aNcrJyVHPnj2VnZ3tVG/06NE6evSoOc2aNcvVXQEAoNKr5uoGU1JSnOYXL16s0NBQ7dixQ126dDHLfX19FRYW5urNAwBQpbj9HHVmZqYkKTg42Kl8yZIlqlWrllq3bq2kpCSdOXOm2DYcDoeysrKcJgAArgYuH1FfKD8/XxMnTtQtt9yi1q1bm+X33HOP6tevr4iICO3atUtTpkzR7t279f777xfZTnJysqZPn+7OrgIAYEk2wzAMdzX+4IMPatWqVfrss89Ut27dYuutX79e3bt31969e9W4ceNCyx0OhxwOhzmflZWlyMhIZWZmKiAgwC19B652DR79yGn+wMy+FdQToOrJyspSYGBgmXLMbSPqcePGaeXKldq8eXOJIS1JHTt2lKRig9put8tut7ulnwAAWJnLg9owDI0fP14ffPCBNm7cqIYNG5a6zs6dOyVJ4eHhru4OAACVmsuDOiEhQUuXLtWKFSvk7++vtLQ0SVJgYKB8fHy0b98+LV26VH369FFISIh27dqlSZMmqUuXLmrbtq2ruwMAQKXm8qCeP3++pD8eanKhRYsWafjw4fLy8tLatWs1Z84cZWdnKzIyUnFxcXrsscdc3RUAACo9txz6LklkZKQ2bdrk6s0CAFAl8axvAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMJcHtTJycm68cYb5e/vr9DQUA0YMEC7d+92qnPu3DklJCQoJCREfn5+iouLU3p6uqu7AgBApefyoN60aZMSEhK0detWrVmzRjk5OerZs6eys7PNOpMmTdJ///tfLVu2TJs2bdKRI0c0cOBAV3cFAIBKr5qrG0xJSXGaX7x4sUJDQ7Vjxw516dJFmZmZeuWVV7R06VLdeuutkqRFixapRYsW2rp1q2666SZXdwkAgErL7eeoMzMzJUnBwcGSpB07dignJ0c9evQw6zRv3lz16tVTampqkW04HA5lZWU5TQAAXA3cGtT5+fmaOHGibrnlFrVu3VqSlJaWJi8vLwUFBTnVrVOnjtLS0opsJzk5WYGBgeYUGRnpzm4DAGAZbg3qhIQEffvtt3rrrbcuq52kpCRlZmaa0+HDh13UQwAArM3l56gLjBs3TitXrtTmzZtVt25dszwsLEznz59XRkaG06g6PT1dYWFhRbZlt9tlt9vd1VUAACzL5SNqwzA0btw4ffDBB1q/fr0aNmzotDwqKkrVq1fXunXrzLLdu3fr0KFDio6OdnV3AACo1Fw+ok5ISNDSpUu1YsUK+fv7m+edAwMD5ePjo8DAQI0aNUqJiYkKDg5WQECAxo8fr+joaK74BgDgIi4P6vnz50uSYmJinMoXLVqk4cOHS5Kee+45eXh4KC4uTg6HQ7169dK8efNc3RUAACo9lwe1YRil1vH29tbcuXM1d+5cV28eAIAqhWd9AwBgYW676htA1dLg0Y+c5g/M7FtBPQGuLoyoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAAL42sugasUX1sJVA6MqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsDDuowYgqfB91QCsgRE1AAAWRlADAGBhBDUAABbm8qDevHmz+vXrp4iICNlsNi1fvtxp+fDhw2Wz2Zym3r17u7obAABUCS4P6uzsbLVr105z584ttk7v3r119OhRc/rPf/7j6m4AAFAluPyq79jYWMXGxpZYx263KywszNWbBgCgyqmQc9QbN25UaGiomjVrpgcffFAnTpwosb7D4VBWVpbTBADA1eCKB3Xv3r31+uuva926dXr66ae1adMmxcbGKi8vr9h1kpOTFRgYaE6RkZFXsMcAAFScK/7AkyFDhpj/b9Omjdq2bavGjRtr48aN6t69e5HrJCUlKTEx0ZzPysoirAEAV4UKvz2rUaNGqlWrlvbu3VtsHbvdroCAAKcJAICrQYUH9S+//KITJ04oPDy8orsCAIDluPzQ9+nTp51Gx/v379fOnTsVHBys4OBgTZ8+XXFxcQoLC9O+ffs0efJkNWnSRL169XJ1VwAAqPRcHtTbt29Xt27dzPmCc8vx8fGaP3++du3apddee00ZGRmKiIhQz5499eSTT8put7u6KwAAVHouD+qYmBgZhlHs8tWrV7t6kwAAVFkVfo4aAAAUj6AGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACzsin/NJYCK0eDRj9za3oGZfV3aPoA/MKIGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCXB7UmzdvVr9+/RQRESGbzably5c7LTcMQ0888YTCw8Pl4+OjHj16aM+ePa7uBgAAVYLLgzo7O1vt2rXT3Llzi1w+a9YsvfDCC1qwYIG2bdumGjVqqFevXjp37pyruwIAQKVXzdUNxsbGKjY2tshlhmFozpw5euyxx3THHXdIkl5//XXVqVNHy5cv15AhQ4pcz+FwyOFwmPNZWVmu7jYAAJZ0Rc9R79+/X2lpaerRo4dZFhgYqI4dOyo1NbXY9ZKTkxUYGGhOkZGRV6K7AABUuCsa1GlpaZKkOnXqOJXXqVPHXFaUpKQkZWZmmtPhw4fd2k8AAKzC5Ye+3cFut8tut1d0NwAAuOKu6Ig6LCxMkpSenu5Unp6ebi4DAAD/54oGdcOGDRUWFqZ169aZZVlZWdq2bZuio6OvZFcAAKgUXH7o+/Tp09q7d685v3//fu3cuVPBwcGqV6+eJk6cqKeeekpNmzZVw4YN9fjjjysiIkIDBgxwdVcAAKj0XB7U27dvV7du3cz5xMRESVJ8fLwWL16syZMnKzs7W2PGjFFGRoY6deqklJQUeXt7u7orAABUejbDMIyK7kR5ZWVlKTAwUJmZmQoICKjo7gCVQoNHP3Jr+wdm9nVr+0BVUp4c41nfAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIVVq+gOAHCPBo9+VNFdAOACjKgBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAAL4/YsoIqo6NuxLt7+gZl9K6gnQNXCiBoAAAsjqAEAsDCCGgAACyOoAQCwsAoJ6mnTpslmszlNzZs3r4iuAABgaRV21XerVq20du3a/+tINS5ABwDgYhWWjtWqVVNYWFiZ6jocDjkcDnM+KyvLXd0CAMBSKiyo9+zZo4iICHl7eys6OlrJycmqV69ekXWTk5M1ffr0K9xDAJeD+6oB16iQc9QdO3bU4sWLlZKSovnz52v//v3q3LmzTp06VWT9pKQkZWZmmtPhw4evcI8BAKgYFTKijo2NNf/ftm1bdezYUfXr19c777yjUaNGFapvt9tlt9uvZBcBALAES9yeFRQUpGuvvVZ79+6t6K4AAGAplgjq06dPa9++fQoPD6/orgAAYCkVEtR/+ctftGnTJh04cEBbtmzRnXfeKU9PTw0dOrQiugMAgGVVyDnqX375RUOHDtWJEydUu3ZtderUSVu3blXt2rUrojsAAFhWhQT1W2+9VRGbBQCg0rHEOWoAAFA0ghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAursO+jBnB14fupgUvDiBoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAvjPmqgkrr4vmQAVRMjagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMK4PQtAheBrL4GyYUQNAICFEdQAAFgYQQ0AgIUR1AAAWFiFBfXcuXPVoEEDeXt7q2PHjvriiy8qqisAAFhWhQT122+/rcTERE2dOlVfffWV2rVrp169eunYsWMV0R0AACyrQm7PevbZZzV69GiNGDFCkrRgwQJ99NFHevXVV/Xoo48Wqu9wOORwOMz5zMxMSVJWVtaV6TBgQfmOMxXdBZfi5xlXk4LPu2EYpVc2rjCHw2F4enoaH3zwgVP5fffdZ/Tv37/IdaZOnWpIYmJiYmJiqlLT4cOHS83NKz6i/u2335SXl6c6deo4ldepU0c//vhjkeskJSUpMTHRnM/Pz9fJkycVEhIim8122X3KyspSZGSkDh8+rICAgMturyKwD9ZRFfaDfbCOqrAf7ENhhmHo1KlTioiIKLVupXgymd1ul91udyoLCgpy+XYCAgIq7YeoAPtgHVVhP9gH66gK+8E+OAsMDCxTvSt+MVmtWrXk6emp9PR0p/L09HSFhYVd6e4AAGBpVzyovby8FBUVpXXr1pll+fn5WrdunaKjo690dwAAsLQKOfSdmJio+Ph43XDDDerQoYPmzJmj7Oxs8yrwK81ut2vq1KmFDq9XJuyDdVSF/WAfrKMq7Af7cHlshlGWa8Nd76WXXtLs2bOVlpam9u3b64UXXlDHjh0roisAAFhWhQU1AAAoHc/6BgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsjKAGAMDCCGoAACyMoAYAwMIIaqhBgwa6/fbbK7obFcZms2natGmXtO7GjRtls9m0ceNGs2z48OFq0KCBS/pWmgYNGmj48OHm/OLFi2Wz2bR9+/Yrsv2YmBjFxMRckW0V5Z133lFwcLBOnz5dYX242l38GXSllJQU+fn56fjx425pv7IgqC2i4BfshVNoaKi6deumVatWVXT34Gbff/+9pk2bpgMHDlR0Vwqxat/y8vI0depUjR8/Xn5+fmb5J598olGjRql169by9PQs8Y+m/Px8zZo1Sw0bNpS3t7fatm2r//znP0XW/eGHH9S7d2/5+fkpODhYf/rTn4oMkPK0+dJLL6lFixay2+265pprlJiYqOzs7PK9EC5ixfe5d+/eatKkiZKTkyu6KxXLgCUsWrTIkGTMmDHDeOONN4zXX3/dmD17ttGqVStDkvHf//7XbduuX7++0bdvX7e1b3WSjKlTp17Suhs2bDAkGRs2bDDLzp8/b5w7d65c7SxbtqxQO2Vx7tw54/z58+Z8wefoyy+/LFc7l9o3h8NhOBwOl22rPD744APDZrMZv/zyi1N5fHy84e3tbdx8881G3bp1jfr16xfbxqOPPmpIMkaPHm0sXLjQ6Nu3ryHJ+M9//uNU7/Dhw0atWrWMxo0bG88//7zx97//3ahZs6bRrl27Qvtf1jYnT55sSDIGDRpkzJ8/3xg/frxRrVo1o2fPnpf3wlwiV30GXW3evHmGr6+vkZWV5bZtWB1BbRHF/YI9efKkUb16deOee+5x27avRFCfPn3are1fDlcH9aUozy/J/Px848yZM0Uuu9JBXZH69+9vdOrUqVD5r7/+agZH3759iw3qX375xahevbqRkJBgluXn5xudO3c26tata+Tm5prlDz74oOHj42McPHjQLFuzZo0hyXj55ZfL3eaRI0eMatWqGX/605+c+vTiiy8akowPP/ywHK+Ea7jqM+hq6enphqenp/HKK69cke1ZEYe+LS4oKEg+Pj6qVq2aU/kzzzyjm2++WSEhIfLx8VFUVJTefffdItt488031aFDB/n6+qpmzZrq0qWLPvnkkxK3+9prr6latWp65JFHzLITJ07oT3/6kwICAhQUFKT4+Hh9/fXXstlsWrx4sVlv+PDh8vPz0759+9SnTx/5+/tr2LBhkqTs7Gz9+c9/VmRkpOx2u5o1a6ZnnnlGhmGY6x84cKBQmwUuPp88bdo02Ww27d27V8OHD1dQUJACAwM1YsQInTlzxmldh8OhSZMmqXbt2vL391f//v31yy+/lPg6XOiXX37RgAEDVKNGDYWGhmrSpElyOByF6hV1jvqtt95SVFSU/P39FRAQoDZt2uj555+X9Mdpj7vuukuS1K1bN/PUR8F574JrCFavXq0bbrhBPj4+evnll81lRZ0fPHPmjMaOHauQkBAFBATovvvu0++//17ia1ngwjZL61tR56iPHTumUaNGqU6dOvL29la7du302muvOdUpeI+feeYZLVy4UI0bN5bdbteNN96oL7/8slCfLnbu3DmlpKSoR48ehZZFRESoevXqpbaxYsUK5eTk6KGHHjLLbDabHnzwQf3yyy9KTU01y9977z3dfvvtqlevnlnWo0cPXXvttXrnnXfK3WZqaqpyc3M1ZMgQpz4VzL/11lul9j8/P19z5sxRq1at5O3trTp16mjs2LGF3ueCz89nn32mDh06yNvbW40aNdLrr79u1nH1ZzAjI0MTJ040f86bNGmip59+Wvn5+U71Svq5KBAaGqq2bdtqxYoVpb4mVVW10qvgSsrMzNRvv/0mwzB07Ngxvfjiizp9+rTuvfdep3rPP/+8+vfvr2HDhun8+fN66623dNddd2nlypXq27evWW/69OmaNm2abr75Zs2YMUNeXl7atm2b1q9fr549exbZh4ULF+qBBx7QX//6Vz311FOS/vil0K9fP33xxRd68MEH1bx5c61YsULx8fFFtpGbm6tevXqpU6dOeuaZZ+Tr6yvDMNS/f39t2LBBo0aNUvv27bV69Wo98sgj+vXXX/Xcc89d8us2ePBgNWzYUMnJyfrqq6/073//W6GhoXr66afNOvfff7/efPNN3XPPPbr55pu1fv16p9eqJGfPnlX37t116NAhTZgwQREREXrjjTe0fv36Utdds2aNhg4dqu7du5v9+eGHH/T555/r4YcfVpcuXTRhwgS98MIL+utf/6oWLVpIkvmvJO3evVtDhw7V2LFjNXr0aDVr1qzEbY4bN05BQUGaNm2adu/erfnz5+vgwYPmxW9lVZa+Xejs2bOKiYnR3r17NW7cODVs2FDLli3T8OHDlZGRoYcfftip/tKlS3Xq1CmNHTtWNptNs2bN0sCBA/Xzzz+XGLY7duzQ+fPndf3115d5Xy72v//9TzVq1Ci0Lx06dDCXd+rUSb/++quOHTumG264oVAbHTp00Mcff1zuNgv+wPPx8XGq5+vra+5facaOHavFixdrxIgRmjBhgvbv36+XXnpJ//vf//T55587vX579+7VoEGDNGrUKMXHx+vVV1/V8OHDFRUVpVatWrn0M3jmzBl17dpVv/76q8aOHat69eppy5YtSkpK0tGjRzVnzhxJpf9cXCgqKkrLly8v9TWpsip4RI//r+CQ5cWT3W43Fi9eXKj+xYedzp8/b7Ru3dq49dZbzbI9e/YYHh4exp133mnk5eU51c/Pzzf/f+Gh7+eff96w2WzGk08+6VT/vffeMyQZc+bMMcvy8vKMW2+91ZBkLFq0yCyPj483JBmPPvqoUxvLly83JBlPPfWUU/mgQYMMm81m7N271zAMw9i/f3+hNgvoosPUU6dONSQZI0eOdKp35513GiEhIeb8zp07DUnGQw895FTvnnvuKdOh7zlz5hiSjHfeeccsy87ONpo0aVLocGF8fLzT4daHH37YCAgIcDqUerGSDjvWr1/fkGSkpKQUuSw+Pt6cL/gcRUVFOZ03nDVrliHJWLFihVlW3H5f3GZJfevatavRtWtXc77gdXrzzTfNsvPnzxvR0dGGn5+feZ6x4D0OCQkxTp48adZdsWJFma7J+Pe//21IMr755psS65V06Ltv375Go0aNCpVnZ2c7fX6//PJLQ5Lx+uuvF6r7yCOPGJLMaxLK2uaOHTsMSYV+zlJSUgxJhp+fX4n79emnnxqSjCVLlhS5/oXlBZ+fzZs3m2XHjh0z7Ha78ec//9ksc9Vn8MknnzRq1Khh/PTTT071Hn30UcPT09M4dOiQYRhl+7ko8I9//MOQZKSnp5datyri0LfFzJ07V2vWrNGaNWv05ptvqlu3brr//vv1/vvvO9W78C/x33//XZmZmercubO++uors3z58uXKz8/XE088IQ8P57e6qFHVrFmz9PDDD+vpp5/WY4895rQsJSVF1atX1+jRo80yDw8PJSQkFLsvDz74oNP8xx9/LE9PT02YMMGp/M9//rMMw7isq9sfeOABp/nOnTvrxIkTysrKMrctqdC2J06cWKb2P/74Y4WHh2vQoEFmma+vr8aMGVPqukFBQcrOztaaNWvKtK2iNGzYUL169Spz/TFjxjiNqB588EFVq1bNafTnDh9//LHCwsI0dOhQs6x69eqaMGGCTp8+rU2bNjnVv/vuu1WzZk1zvnPnzpKkn3/+ucTtnDhxQpKc1i2vs2fPym63Fyr39vY2l1/4b1nrlqXe9ddfr44dO+rpp5/WokWLdODAAa1atUpjx45V9erVzXrFWbZsmQIDA3Xbbbfpt99+M6eoqCj5+flpw4YNTvVbtmxpvraSVLt2bTVr1qzU1/lCZf0MLlu2TJ07d1bNmjWd+tajRw/l5eVp8+bNksr3c1HwPv/2229l7m9VwqFvi+nQoYPTIbahQ4fquuuu07hx43T77bfLy8tLkrRy5Uo99dRT2rlzp9N50gsDeN++ffLw8FDLli1L3e6mTZv00UcfacqUKU7npQscPHhQ4eHh5qG5Ak2aNCmyvWrVqqlu3bqF2oiIiJC/v79TecHhtYMHD5baz+JceO5Q+r8f7N9//10BAQE6ePCgPDw81LhxY6d6pR1CLnDw4EE1adKk0B84ZVn/oYce0jvvvKPY2Fhdc8016tmzpwYPHqzevXuXadvSH78ky6Np06ZO835+fgoPD3f7rTcHDx5U06ZNC/1hWNx7XNL7VhbGBdc2lJePj0+R1xicO3fOXH7hv2WtW5Z60h/nve+++26NHDlSkuTp6anExERt2rRJu3fvLrHve/bsUWZmpkJDQ4tcfuzYMaf5i19n6Y/Xuqyvs1T2z+CePXu0a9cu1a5du8S+lefnouB9Ls9pm6qEoLY4Dw8PdevWTc8//7z27NmjVq1a6dNPP1X//v3VpUsXzZs3T+Hh4apevboWLVqkpUuXXtJ2WrVqpYyMDL3xxhsaO3ZsuYPhYna7vdAv67Iq7ocxLy+v2HU8PT2LLL+cX+SuEhoaqp07d2r16tVatWqVVq1apUWLFum+++4rdJFVcS4+l+lOJb3Ornap71tISIikPwL94j8Iyyo8PFwbNmyQYRhOn7mjR49K+uOitIJ6F5Zf6OjRowoODjZH0WVtU5KuueYaffbZZ9qzZ4/S0tLUtGlThYWFKSIiQtdee22Jfc/Pz1doaKiWLFlS5PKLQ9IVPx9l/Qzm5+frtttu0+TJk4tcXrBv5fm5KPiDolatWmXub1XCoe9KIDc3V5LMpy+999578vb21urVqzVy5EjFxsYWefVr48aNlZ+fr++//77UbdSqVUtr165V9erV1b17dx05csRpef369XX06NFCV1Lv3bu3zPtRv359HTlyRKdOnXIq//HHH83l0v+NqjIyMpzqXc6Iu379+srPz9e+ffucyksbuVy4/r59+wr9Yivr+l5eXurXr5/mzZunffv2aezYsXr99dfN18/VI4U9e/Y4zZ8+fVpHjx51uhq9Zs2ahV7j8+fPFwqk8vStfv362rNnT6Grey9+jy9X8+bNJUn79++/5Dbat2+vM2fO6IcffnAq37Ztm7lc+iNQa9euXeTT3r744guzXnnavFDTpk3VuXNnhYWF6fvvv9fRo0eL/Hm+UOPGjXXixAndcsst6tGjR6GpXbt2pe1+Ia76DDZu3FinT58usl89evRwGt2X9nNRYP/+/apVq1axo/SqjqC2uJycHH3yySfy8vIyDx96enrKZrM5jXwOHDhQ6KrIAQMGyMPDQzNmzCj0i7Oov6Tr1q2rtWvX6uzZs7rtttvM84CS1KtXL+Xk5Ohf//qXWZafn6+5c+eWeV/69OmjvLw8vfTSS07lzz33nGw2m2JjYyVJAQEBqlWrlnkuq8C8efPKvK2LFbT9wgsvOJUXXIFamj59+ujIkSNOt8CdOXNGCxcuLHXdC19H6Y+jJG3btpX0f4dTa9SoIanwHyeXauHChcrJyTHn58+fr9zcXPN1kP74hXrxa7xw4cJCI+ry9K1Pnz5KS0vT22+/bZbl5ubqxRdflJ+fn7p27Xopu1NIVFSUvLy8LutRqXfccYeqV6/u9LkyDEMLFizQNddco5tvvtksj4uL08qVK3X48GGzbN26dfrpp5/M25rK2+bF8vPzNXnyZPn6+ha65uJigwcPVl5enp588slCy3Jzcy/pc+Sqz+DgwYOVmpqq1atXF1qWkZFhDjzK8nNRYMeOHYqOjr6sflVmHPq2mFWrVpmjj2PHjmnp0qXas2ePHn30UQUEBEiS+vbtq2effVa9e/fWPffco2PHjmnu3Llq0qSJdu3aZbbVpEkT/e1vf9OTTz6pzp07a+DAgbLb7fryyy8VERFR5GP5mjRpok8++UQxMTHq1auX1q9fr4CAAA0YMEAdOnTQn//8Z+3du1fNmzfXhx9+qJMnT0oq21/j/fr1U7du3fS3v/1NBw4cULt27fTJJ59oxYoVmjhxotP54/vvv18zZ87U/fffrxtuuEGbN2/WTz/9dMmva/v27TV06FDNmzdPmZmZuvnmm7Vu3boyHxEYPXq0XnrpJd13333asWOHwsPD9cYbbxQ6Z1+U+++/XydPntStt96qunXr6uDBg3rxxRfVvn1784+v9u3by9PTU08//bQyMzNlt9t16623FnsOsjTnz59X9+7dNXjwYO3evVvz5s1Tp06d1L9/f6d+PfDAA4qLi9Ntt92mr7/+WqtXry50eLE8fRszZoxefvllDR8+XDt27FCDBg307rvv6vPPP9ecOXMKXZ9wqby9vdWzZ0+tXbtWM2bMcFq2a9cuffjhh5L+OOKTmZlp3mbYrl079evXT9Iff5hOnDhRs2fPVk5Ojm688UYtX75cn376qZYsWeJ0uPivf/2rli1bpm7duunhhx/W6dOnNXv2bLVp00YjRoww65WnzYcffljnzp1T+/btlZOTo6VLl+qLL77Qa6+9VuQ55Qt17dpVY8eOVXJysnbu3KmePXuqevXq2rNnj5YtW6bnn3/e6cLHsnDVZ/CRRx7Rhx9+qNtvv928BSw7O1vffPON3n33XR04cEC1atUq08+F9MfvwV27dpV44WqVVzEXm+NiRd2e5e3tbbRv396YP3++0+1UhmEYr7zyitG0aVPDbrcbzZs3NxYtWmTeqnSxV1991bjuuusMu91u1KxZ0+jatauxZs0ac3lRTybbtm2b4e/vb3Tp0sW8Fez48ePGPffcY/j7+xuBgYHG8OHDjc8//9yQZLz11lvmuvHx8UaNGjWK3M9Tp04ZkyZNMiIiIozq1asbTZs2NWbPnl1o/86cOWOMGjXKCAwMNPz9/Y3Bgwcbx44dK/b2rOPHjxf5eu7fv98sO3v2rDFhwgQjJCTEqFGjhtGvXz/j8OHDZX4y2cGDB43+/fsbvr6+Rq1atYyHH37YvB2mpNuz3n33XaNnz55GaGio4eXlZdSrV88YO3ascfToUaf2//WvfxmNGjUyPD09ndos6clxxd2etWnTJmPMmDFGzZo1DT8/P2PYsGHGiRMnnNbNy8szpkyZYtSqVcvw9fU1evXqZezdu7dQmyX17eLbswzjjydJjRgxwqhVq5bh5eVltGnTptCtdgW3Z82ePbvQPpX1/Xj//fcNm81m3u5z8WtQ1HTxfuXl5Rn/+Mc/jPr16xteXl5Gq1atnG4tu9C3335r9OzZ0/D19TWCgoKMYcOGGWlpaYXqlbXNRYsWGe3atTNq1Khh+Pv7G927dzfWr19f6n5faOHChUZUVJTh4+Nj+Pv7G23atDEmT55sHDlyxKxT3OenqPfOFZ9Bw/jj5zwpKclo0qSJ4eXlZdSqVcu4+eabjWeeeca8bbCsPxfz58+/6h8hajMMC1xtg0pr+fLluvPOO/XZZ5/plltuqeju4CqSl5enli1bavDgwUUeAkbVcN111ykmJuayHohU2RHUKLOzZ886XfmZl5ennj17avv27UpLS7uiVyYDkvT222/rwQcf1KFDh5y+QQtVQ0pKigYNGqSff/75kk8DVQUENcrs/vvv19mzZxUdHS2Hw6H3339fW7Zs0T/+8Q8lJSVVdPcAoEoiqFFmS5cu1T//+U/t3btX586dU5MmTfTggw9q3LhxFd01AKiyCGoAACyM+6gBALCwSnkfdX5+vo4cOSJ/f/+r9tmvAIDKyzAMnTp1ShEREaU+brlSBvWRI0cUGRlZ0d0AAOCyHD58uNTn1VfKoC54utHhw4fNp3UBAFBZZGVlKTIyskxP66uUQV1wuDsgIICgBgBUWmU5fcvFZAAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFVcpHiF4JDR79yGn+wMy+FdQTAMDVjBE1AAAWVq6gTk5O1o033ih/f3+FhoZqwIAB2r17t1Odc+fOKSEhQSEhIfLz81NcXJzS09Od6hw6dEh9+/aVr6+vQkND9cgjjyg3N/fy9wYAgCqmXEG9adMmJSQkaOvWrVqzZo1ycnLUs2dPZWdnm3UmTZqk//73v1q2bJk2bdqkI0eOaODAgebyvLw89e3bV+fPn9eWLVv02muvafHixXriiSdct1cAAFQRNsMwjEtd+fjx4woNDdWmTZvUpUsXZWZmqnbt2lq6dKkGDRokSfrxxx/VokULpaam6qabbtKqVat0++2368iRI6pTp44kacGCBZoyZYqOHz8uLy+vQttxOBxyOBzmfMH3eGZmZrrtay45Rw0AcJesrCwFBgaWKccu6xx1ZmamJCk4OFiStGPHDuXk5KhHjx5mnebNm6tevXpKTU2VJKWmpqpNmzZmSEtSr169lJWVpe+++67I7SQnJyswMNCcIiMjL6fbAABUGpd81Xd+fr4mTpyoW265Ra1bt5YkpaWlycvLS0FBQU5169Spo7S0NLPOhSFdsLxgWVGSkpKUmJhozheMqF3p4hE0AABWcMlBnZCQoG+//VafffaZK/tTJLvdLrvd7vbtAABgNZd06HvcuHFauXKlNmzYoLp165rlYWFhOn/+vDIyMpzqp6enKywszKxz8VXgBfMFdQAAwB/KFdSGYWjcuHH64IMPtH79ejVs2NBpeVRUlKpXr65169aZZbt379ahQ4cUHR0tSYqOjtY333yjY8eOmXXWrFmjgIAAtWzZ8nL2BQCAKqdch74TEhK0dOlSrVixQv7+/uY55cDAQPn4+CgwMFCjRo1SYmKigoODFRAQoPHjxys6Olo33XSTJKlnz55q2bKl/vSnP2nWrFlKS0vTY489poSEBA5vAwBwkXIF9fz58yVJMTExTuWLFi3S8OHDJUnPPfecPDw8FBcXJ4fDoV69emnevHlmXU9PT61cuVIPPvigoqOjVaNGDcXHx2vGjBmXtydXELduAQCulHIFdVluufb29tbcuXM1d+7cYuvUr19fH3/8cXk2DQDAVYlnfQMAYGEENQAAFkZQAwBgYXwftRtwsRkAwFUYUQMAYGEENQAAFsah7zKqrF/awWF4AKjcGFEDAGBhBDUAABZGUAMAYGEENQAAFsbFZJUAF4QBwNWLETUAABZGUAMAYGEENQAAFkZQAwBgYVxMVgG4OAwAUFYEtQUQ3ACA4nDoGwAACyOoAQCwsHIH9ebNm9WvXz9FRETIZrNp+fLlTsttNluR0+zZs806DRo0KLR85syZl70zVUWDRz9ymgAAV69yn6POzs5Wu3btNHLkSA0cOLDQ8qNHjzrNr1q1SqNGjVJcXJxT+YwZMzR69Ghz3t/fv7xdgQtwfhwArK3cQR0bG6vY2Nhil4eFhTnNr1ixQt26dVOjRo2cyv39/QvVBQAAztx6jjo9PV0fffSRRo0aVWjZzJkzFRISouuuu06zZ89Wbm5use04HA5lZWU5TQAAXA3cenvWa6+9Jn9//0KHyCdMmKDrr79ewcHB2rJli5KSknT06FE9++yzRbaTnJys6dOnu7OrlRaHrgGganNrUL/66qsaNmyYvL29ncoTExPN/7dt21ZeXl4aO3askpOTZbfbC7WTlJTktE5WVpYiIyPd1/FyKu2CLy4IAwBcKrcF9aeffqrdu3fr7bffLrVux44dlZubqwMHDqhZs2aFltvt9iIDHACAqs5t56hfeeUVRUVFqV27dqXW3blzpzw8PBQaGuqu7gAAUCmVe0R9+vRp7d2715zfv3+/du7cqeDgYNWrV0/SH4emly1bpn/+85+F1k9NTdW2bdvUrVs3+fv7KzU1VZMmTdK9996rmjVrXsauAABQ9ZQ7qLdv365u3bqZ8wXnjuPj47V48WJJ0ltvvSXDMDR06NBC69vtdr311luaNm2aHA6HGjZsqEmTJjmdg0bFKe18OherAcCVVe6gjomJkWEYJdYZM2aMxowZU+Sy66+/Xlu3bi3vZgEAuCrxrG8AACyMoAYAwMIIagAALMytDzzB1YWnpAGA6zGiBgDAwhhRVzE8rhQAqhZG1AAAWBgjapTLhSN2zkEDgPsxogYAwMIIagAALIygBgDAwghqAAAsjIvJKiFuwQKAqwdBDcvgyWYAUBiHvgEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsLByB/XmzZvVr18/RUREyGazafny5U7Lhw8fLpvN5jT17t3bqc7Jkyc1bNgwBQQEKCgoSKNGjdLp06cva0cAAKiKyv3Ak+zsbLVr104jR47UwIEDi6zTu3dvLVq0yJy32+1Oy4cNG6ajR49qzZo1ysnJ0YgRIzRmzBgtXbq0vN1BJcIDTQCg/Mod1LGxsYqNjS2xjt1uV1hYWJHLfvjhB6WkpOjLL7/UDTfcIEl68cUX1adPHz3zzDOKiIgotI7D4ZDD4TDns7KyytttAAAqJbeco964caNCQ0PVrFkzPfjggzpx4oS5LDU1VUFBQWZIS1KPHj3k4eGhbdu2FdlecnKyAgMDzSkyMtId3QYAwHJcHtS9e/fW66+/rnXr1unpp5/Wpk2bFBsbq7y8PElSWlqaQkNDndapVq2agoODlZaWVmSbSUlJyszMNKfDhw+7utsAAFiSy7+UY8iQIeb/27Rpo7Zt26px48bauHGjunfvfklt2u32Que5AQC4Grj99qxGjRqpVq1a2rt3ryQpLCxMx44dc6qTm5urkydPFnteGwCAq5Xbv+byl19+0YkTJxQeHi5Jio6OVkZGhnbs2KGoqChJ0vr165Wfn6+OHTu6uzu4gvjebAC4fOUO6tOnT5ujY0nav3+/du7cqeDgYAUHB2v69OmKi4tTWFiY9u3bp8mTJ6tJkybq1auXJKlFixbq3bu3Ro8erQULFignJ0fjxo3TkCFDirziGwCAq1m5D31v375d1113na677jpJUmJioq677jo98cQT8vT01K5du9S/f39de+21GjVqlKKiovTpp586nWNesmSJmjdvru7du6tPnz7q1KmTFi5c6Lq9AgCgiij3iDomJkaGYRS7fPXq1aW2ERwczMNNKgiHowGgcnH7OWpUXYQ+ALgfQY1K48I/DHj8KICrBd+eBQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGA88wVXp4qeqlfYAFR62AqCiENSwrJIeUVreoAWAyopD3wAAWBgjalwV+AIRAJUVQQ2UE4fdAVxJBDWqBEbMAKoqghoQo2QA1sXFZAAAWBgjasDFGJ0DcCVG1AAAWFi5g3rz5s3q16+fIiIiZLPZtHz5cnNZTk6OpkyZojZt2qhGjRqKiIjQfffdpyNHjji10aBBA9lsNqdp5syZl70zAErW4NGPnCYA1lfuQ9/Z2dlq166dRo4cqYEDBzotO3PmjL766is9/vjjateunX7//Xc9/PDD6t+/v7Zv3+5Ud8aMGRo9erQ57+/vf4m7ALgeIQbAKsod1LGxsYqNjS1yWWBgoNasWeNU9tJLL6lDhw46dOiQ6tWrZ5b7+/srLCysvJsHAOCq4vZz1JmZmbLZbAoKCnIqnzlzpkJCQnTddddp9uzZys3NLbYNh8OhrKwspwkAgKuBW6/6PnfunKZMmaKhQ4cqICDALJ8wYYKuv/56BQcHa8uWLUpKStLRo0f17LPPFtlOcnKypk+f7s6uAgBgSW4L6pycHA0ePFiGYWj+/PlOyxITE83/t23bVl5eXho7dqySk5Nlt9sLtZWUlOS0TlZWliIjI93VdQAALMMtQV0Q0gcPHtT69eudRtNF6dixo3Jzc3XgwAE1a9as0HK73V5kgAMAUNW5PKgLQnrPnj3asGGDQkJCSl1n586d8vDwUGhoqKu7AwBApVbuoD59+rT27t1rzu/fv187d+5UcHCwwsPDNWjQIH311VdauXKl8vLylJaWJkkKDg6Wl5eXUlNTtW3bNnXr1k3+/v5KTU3VpEmTdO+996pmzZqu2zPgCuFWLgDuVO6g3r59u7p162bOF5w7jo+P17Rp0/Thhx9Kktq3b++03oYNGxQTEyO73a633npL06ZNk8PhUMOGDTVp0iSnc9AAAOAP5Q7qmJgYGYZR7PKSlknS9ddfr61bt5Z3swCKwHPFgaqPZ30DAGBhBDUAABZGUAMAYGF8HzUqDFdLA0DpCGqggnFBGICScOgbAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMG7PAtyM268AXA5G1AAAWBgjauAKK+8T2S6sz2gcuPoQ1EAlUlrIc5gdqHo49A0AgIUxogZgSRwdAP7AiBoAAAsjqAEAsDCCGgAACyOoAQCwsHIH9ebNm9WvXz9FRETIZrNp+fLlTssNw9ATTzyh8PBw+fj4qEePHtqzZ49TnZMnT2rYsGEKCAhQUFCQRo0apdOnT1/WjgAorMGjHzlNVmkLQNmVO6izs7PVrl07zZ07t8jls2bN0gsvvKAFCxZo27ZtqlGjhnr16qVz586ZdYYNG6bvvvtOa9as0cqVK7V582aNGTPm0vcCwCUhfAHrK/ftWbGxsYqNjS1ymWEYmjNnjh577DHdcccdkqTXX39dderU0fLlyzVkyBD98MMPSklJ0ZdffqkbbrhBkvTiiy+qT58+euaZZxQREXEZuwMAQNXi0nPU+/fvV1pamnr06GGWBQYGqmPHjkpNTZUkpaamKigoyAxpSerRo4c8PDy0bdu2Itt1OBzKyspymgAAuBq49IEnaWlpkqQ6deo4ldepU8dclpaWptDQUOdOVKum4OBgs87FkpOTNX36dFd2FbAsDkEDuFCleDJZUlKSEhMTzfmsrCxFRkZWYI8A8OQw4Mpw6aHvsLAwSVJ6erpTeXp6urksLCxMx44dc1qem5urkydPmnUuZrfbFRAQ4DQBAHA1cOmIumHDhgoLC9O6devUvn17SX+Mfrdt26YHH3xQkhQdHa2MjAzt2LFDUVFRkqT169crPz9fHTt2dGV3AFwmDsMDFa/cQX369Gnt3bvXnN+/f7927typ4OBg1atXTxMnTtRTTz2lpk2bqmHDhnr88ccVERGhAQMGSJJatGih3r17a/To0VqwYIFycnI0btw4DRkyhCu+AQC4SLmDevv27erWrZs5X3DuOD4+XosXL9bkyZOVnZ2tMWPGKCMjQ506dVJKSoq8vb3NdZYsWaJx48ape/fu8vDwUFxcnF544QUX7A4AAFVLuYM6JiZGhmEUu9xms2nGjBmaMWNGsXWCg4O1dOnS8m4aAICrDs/6BgDAwghqAAAsjKAGAMDCKsUDTwBYHw9AAdyDoAZQIQh2oGw49A0AgIUxogZgcueTyHjKGXBpGFEDAGBhBDUAABZGUAMAYGGcowZgCZzDBorGiBoAAAsjqAEAsDAOfQOodHhYCq4mjKgBALAwghoAAAsjqAEAsDDOUQOoFEq6fYtz1qjKGFEDAGBhBDUAABbGoW8AbsGTxgDXcPmIukGDBrLZbIWmhIQESVJMTEyhZQ888ICruwEAQJXg8hH1l19+qby8PHP+22+/1W233aa77rrLLBs9erRmzJhhzvv6+rq6GwAAVAkuD+ratWs7zc+cOVONGzdW165dzTJfX1+FhYWVuU2HwyGHw2HOZ2VlXX5HAaAYXEUOK3HrxWTnz5/Xm2++qZEjR8pms5nlS5YsUa1atdS6dWslJSXpzJkzJbaTnJyswMBAc4qMjHRntwEAsAy3Xky2fPlyZWRkaPjw4WbZPffco/r16ysiIkK7du3SlClTtHv3br3//vvFtpOUlKTExERzPisri7AGAFwV3BrUr7zyimJjYxUREWGWjRkzxvx/mzZtFB4eru7du2vfvn1q3Lhxke3Y7XbZ7XZ3dhUAAEty26HvgwcPau3atbr//vtLrNexY0dJ0t69e93VFQAAKi23BfWiRYsUGhqqvn1Lvghj586dkqTw8HB3dQUAgErLLYe+8/PztWjRIsXHx6tatf/bxL59+7R06VL16dNHISEh2rVrlyZNmqQuXbqobdu27ugKAACVmluCeu3atTp06JBGjhzpVO7l5aW1a9dqzpw5ys7OVmRkpOLi4vTYY4+5oxsAAFR6bgnqnj17yjCMQuWRkZHatGmTOzYJAECVxJdyAABgYQQ1AAAWxrdnAbjq8U1fsDJG1AAAWBhBDQCAhRHUAABYGOeoAVR5fG0lKjNG1AAAWBhBDQCAhXHoGwDKgcPouNIIagBwI4Idl4ugBnDVKe8DTnggCioSQQ2gyiFYUZVwMRkAABZGUAMAYGEc+gaASoIL065OjKgBALAwRtQAcBlcPcpl1IyLMaIGAMDCCGoAACzM5UE9bdo02Ww2p6l58+bm8nPnzikhIUEhISHy8/NTXFyc0tPTXd0NAKgSGjz6kTnh6uSWc9StWrXS2rVr/28j1f5vM5MmTdJHH32kZcuWKTAwUOPGjdPAgQP1+eefu6MrAHBFEahwNbcEdbVq1RQWFlaoPDMzU6+88oqWLl2qW2+9VZK0aNEitWjRQlu3btVNN93kju4AgGVwsRjKyy3nqPfs2aOIiAg1atRIw4YN06FDhyRJO3bsUE5Ojnr06GHWbd68uerVq6fU1NRi23M4HMrKynKaAAC4Grh8RN2xY0ctXrxYzZo109GjRzV9+nR17txZ3377rdLS0uTl5aWgoCCnderUqaO0tLRi20xOTtb06dNd3VUAqNTKOzpnNF85uTyoY2Njzf+3bdtWHTt2VP369fXOO+/Ix8fnktpMSkpSYmKiOZ+VlaXIyMjL7isAAFbn9tuzgoKCdO2112rv3r0KCwvT+fPnlZGR4VQnPT29yHPaBex2uwICApwmAACuBm4P6tOnT2vfvn0KDw9XVFSUqlevrnXr1pnLd+/erUOHDik6OtrdXQEAoNJx+aHvv/zlL+rXr5/q16+vI0eOaOrUqfL09NTQoUMVGBioUaNGKTExUcHBwQoICND48eMVHR3NFd8AABTB5UH9yy+/aOjQoTpx4oRq166tTp06aevWrapdu7Yk6bnnnpOHh4fi4uLkcDjUq1cvzZs3z9XdAACgSnB5UL/11lslLvf29tbcuXM1d+5cV28aACodHpCC0vCsbwAALIygBgDAwghqAAAszC3P+gYAXHk8eaxqYkQNAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhXHVNwBcpbhKvHJgRA0AgIUR1AAAWBiHvgEALsdhddchqAEAl41vAXMfghoAqijCs2rgHDUAABbGiBoAIMm955U5Z33pCGoAQKkI2opDUAMAyo3z31cO56gBALAwghoAAAtzeVAnJyfrxhtvlL+/v0JDQzVgwADt3r3bqU5MTIxsNpvT9MADD7i6KwCAy9Dg0Y/MCRXH5UG9adMmJSQkaOvWrVqzZo1ycnLUs2dPZWdnO9UbPXq0jh49ak6zZs1ydVcAAKj0XH4xWUpKitP84sWLFRoaqh07dqhLly5mua+vr8LCwsrUpsPhkMPhMOezsrJc01kAACzO7eeoMzMzJUnBwcFO5UuWLFGtWrXUunVrJSUl6cyZM8W2kZycrMDAQHOKjIx0a58BALAKm2EYhrsaz8/PV//+/ZWRkaHPPvvMLF+4cKHq16+viIgI7dq1S1OmTFGHDh30/vvvF9lOUSPqyMhIZWZmKiAgwCV95RwMAFw5V/t92FlZWQoMDCxTjrn1PuqEhAR9++23TiEtSWPGjDH/36ZNG4WHh6t79+7at2+fGjduXKgdu90uu93uzq4CAGBJbjv0PW7cOK1cuVIbNmxQ3bp1S6zbsWNHSdLevXvd1R0AAColl4+oDcPQ+PHj9cEHH2jjxo1q2LBhqevs3LlTkhQeHu7q7gAAKgEeUVo8lwd1QkKCli5dqhUrVsjf319paWmSpMDAQPn4+Gjfvn1aunSp+vTpo5CQEO3atUuTJk1Sly5d1LZtW1d3BwCASs3lQT1//nxJfzzU5EKLFi3S8OHD5eXlpbVr12rOnDnKzs5WZGSk4uLi9Nhjj7m6KwAAVHpuOfRdksjISG3atMnVmwUAoEri27MAAJZ2tZ+/5ks5AACwMIIaAAALI6gBALAwghoAAAsjqAEAsDCu+gYAVClV7SpxRtQAAFgYI2oAwBVX2lcLl7T8avtaYkbUAABYGEENAICFEdQAAFgYQQ0AgIVxMRkA4KpS2W7fYkQNAICFMaIGAFRplf12LkbUAABYGEENAICFcegbAHBVs/rFZYyoAQCwsAoL6rlz56pBgwby9vZWx44d9cUXX1RUVwAAsKwKCeq3335biYmJmjp1qr766iu1a9dOvXr10rFjxyqiOwAAWFaFBPWzzz6r0aNHa8SIEWrZsqUWLFggX19fvfrqqxXRHQAALOuKX0x2/vx57dixQ0lJSWaZh4eHevToodTU1CLXcTgccjgc5nxmZqYkKSsry2X9yneccVlbAIDKy5XZUto2DMMote4VD+rffvtNeXl5qlOnjlN5nTp19OOPPxa5TnJysqZPn16oPDIy0i19BABcvQLnXLltnTp1SoGBgSXWqRS3ZyUlJSkxMdGcz8/P18mTJxUSEiKbzXbZ7WdlZSkyMlKHDx9WQEDAZbdXEdgH66gK+8E+WEdV2A/2oTDDMHTq1ClFRESUWveKB3WtWrXk6emp9PR0p/L09HSFhYUVuY7dbpfdbncqCwoKcnnfAgICKu2HqAD7YB1VYT/YB+uoCvvBPjgrbSRd4IpfTObl5aWoqCitW7fOLMvPz9e6desUHR19pbsDAIClVcih78TERMXHx+uGG25Qhw4dNGfOHGVnZ2vEiBEV0R0AACyrQoL67rvv1vHjx/XEE08oLS1N7du3V0pKSqELzK4Uu92uqVOnFjq8XpmwD9ZRFfaDfbCOqrAf7MPlsRlluTYcAABUCJ71DQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIUR1LCEjRs3ymazaePGjZe0/rRp02Sz2ZzKGjRooOHDh19+50px4MAB2Ww2LV682CwbPny4/Pz83L7tAjabTdOmTbti27vYQw89pNtuu63Ctg/3fgYWLFigevXqyeFwuKV9lIygvorZbLYyTWUJz3/84x9avny52/vsbh9//HGFBl5JrNq3/fv369///rf++te/OpXPnz9fd911l+rVqyebzVbiH00ZGRkaM2aMateurRo1aqhbt2766quvCtV7++23de+996pp06ay2WyKiYkpsr2CP/yKmrZu3Vqo/pYtW9SpUyf5+voqLCxMEyZM0OnTp8v1OrjKli1bNG3aNGVkZFTI9osyfPhwnT9/Xi+//HJFd+WqVK2iO4CK88YbbzjNv/7661qzZk2h8hYtWpTa1j/+8Q8NGjRIAwYMcGUXL8vu3bvl4VG+v0U//vhjzZ07t1yBWL9+fZ09e1bVq1cvZw/Lp6S+nT17VtWqVcyP8/PPP6+GDRuqW7duTuVPP/20Tp06pQ4dOujo0aPFrp+fn6++ffvq66+/1iOPPKJatWpp3rx5iomJ0Y4dO9S0aVOz7vz587Vjxw7deOONOnHiRKl9mzBhgm688UansiZNmjjN79y5U927d1eLFi307LPP6pdfftEzzzyjPXv2aNWqVWV5CVxqy5Ytmj59uoYPH66goKAyr+fOz4C3t7fi4+P17LPPavz48YWOXsG9COqr2L333us0v3XrVq1Zs6ZQeWVlt9vd2n5ubq7y8/Pl5eUlb29vt26rNBW1/ZycHC1ZskQPPPBAoWWbNm0yR9MlnQZ49913tWXLFi1btkyDBg2SJA0ePFjXXnutpk6dqqVLl5p133jjDV1zzTXy8PBQ69atS+1f586dzTaL89e//lU1a9bUxo0bFRAQIOmP0yajR4/WJ598op49e5a6nYqSn5+v8+fPy9vb2+2fgcGDB2vWrFnasGGDbr31VrduC8449I0SZWdn689//rMiIyNlt9vVrFkzPfPMMzIMw6xjs9mUnZ2t1157zTy8WHCY8+DBg3rooYfUrFkz+fj4KCQkRHfddZcOHDhwyX367LPPdOONN8rb21uNGzcu9nDcxeeoc3JyNH36dDVt2lTe3t4KCQlRp06dtGbNGkl/HN6bO3euuU8Fk/R/56GfeeYZzZkzR40bN5bdbtf3339f5DnqAj///LN69eqlGjVqKCIiQjNmzHB67Yo7N39xmyX1raDs4pH2//73P8XGxiogIEB+fn7q3r17ocO+ixcvls1m0+eff67ExETz0POdd96p48ePF/0GXOCzzz7Tb7/9ph49ehRaVr9+/TKNvN59913VqVNHAwcONMtq166twYMHa8WKFU7nRSMjI8t9lOTUqVPKzc0tcllWVpb5x2lBSEvSfffdJz8/P73zzjultu9wODR16lQ1adJEdrtdkZGRmjx5cqHzuTabTePGjdPy5cvVunVr2e12tWrVSikpKWadadOm6ZFHHpEkNWzY0HyfC35eCtpYsmSJWrVqJbvdbq5f1Gfg119/1ciRI1WnTh1ze6+++mqhfXjxxRfVqlUr+fr6qmbNmrrhhhuc/kCSpKioKAUHB2vFihWlviZwLUbUKJZhGOrfv782bNigUaNGqX379lq9erUeeeQR/frrr3ruueck/THKuf/++9WhQweNGTNGktS4cWNJ0pdffqktW7ZoyJAhqlu3rg4cOKD58+crJiZG33//vXx9fcvVp2+++UY9e/ZU7dq1NW3aNOXm5mrq1KmqU6dOqetOmzZNycnJZl+zsrK0fft2ffXVV7rttts0duxYHTlypMjD/wUWLVqkc+fOacyYMbLb7QoODlZ+fn6RdfPy8tS7d2/ddNNNmjVrllJSUjR16lTl5uZqxowZ5drvsvTtQt999506d+6sgIAATZ48WdWrV9fLL7+smJgYbdq0SR07dnSqP378eNWsWVNTp07VgQMHNGfOHI0bN05vv/12idvZsmWLbDabrrvuunLtz4X+97//6frrry8UwB06dNDChQv1008/qU2bNpfU9ogRI3T69Gl5enqqc+fOmj17tm644QZz+TfffKPc3FynMkny8vJS+/bt9b///a/E9vPz89W/f3999tlnGjNmjFq0aKFvvvlGzz33nH766adC12189tlnev/99/XQQw/J399fL7zwguLi4nTo0CGFhIRo4MCB+umnn/Sf//xHzz33nGrVqiXpjz9cCqxfv17vvPOOxo0bp1q1aqlBgwZF9i09PV033XSTGe61a9fWqlWrNGrUKGVlZWnixImSpH/961+aMGGCBg0apIcffljnzp3Trl27tG3bNt1zzz1ObV5//fX6/PPPS3xN4AYG8P8lJCQYF34kli9fbkgynnrqKad6gwYNMmw2m7F3716zrEaNGkZ8fHyhNs+cOVOoLDU11ZBkvP7662bZhg0bDEnGhg0bSuzjgAEDDG9vb+PgwYNm2ffff294enoaF3+c69ev79Sndu3aGX379i2x/YtfgwL79+83JBkBAQHGsWPHily2aNEisyw+Pt6QZIwfP94sy8/PN/r27Wt4eXkZx48fL3G/i2qzuL4ZhmFIMqZOnWrODxgwwPDy8jL27dtnlh05csTw9/c3unTpYpYtWrTIkGT06NHDyM/PN8snTZpkeHp6GhkZGUVur8C9995rhISElFjHMIr/fBQsGzlyZKHyjz76yJBkpKSkFLleq1atjK5duxa57PPPPzfi4uKMV155xVixYoWRnJxshISEGN7e3sZXX31l1lu2bJkhydi8eXOhNu666y4jLCysxP164403DA8PD+PTTz91Kl+wYIEhyfj888/NMkmGl5eX08/N119/bUgyXnzxRbNs9uzZhiRj//79hbYnyfDw8DC+++67Ipdd+BkYNWqUER4ebvz2229O9YYMGWIEBgaaP5t33HGH0apVqxL3s8CYMWMMHx+fMtWF63DoG8X6+OOP5enpqQkTJjiV//nPf5ZhGGW60MbHx8f8f05Ojk6cOKEmTZooKCioyKt6S5KXl6fVq1drwIABqlevnlneokUL9erVq9T1g4KC9N1332nPnj3l2u6F4uLinEY3pRk3bpz5/4KRzfnz57V27dpL7kNp8vLy9Mknn2jAgAFq1KiRWR4eHq577rlHn332mbKyspzWGTNmjNNh6s6dOysvL08HDx4scVsnTpxQzZo1L6u/Z8+eLfJ6goJzrmfPni13mzfffLPeffddjRw5Uv3799ejjz6qrVu3ymazKSkpyWnbUtHXM3h7e5e67WXLlqlFixZq3ry5fvvtN3MqOIe7YcMGp/o9evQwjzZJUtu2bRUQEKCff/65zPvWtWtXtWzZssQ6hmHovffeU79+/WQYhlPfevXqpczMTPPnLygoSL/88ou+/PLLUrdds2ZNnT17VmfOnClzf3H5CGoU6+DBg4qIiJC/v79TecFV4KX9Epf++EX4xBNPmOe4a9Wqpdq1aysjI0OZmZnl6s/x48d19uxZp6uACzRr1qzU9WfMmKGMjAxde+21atOmjR555BHt2rWrXH1o2LBhmet6eHg4BaUkXXvttZJ0WefoS3P8+HGdOXOmyNekRYsWys/P1+HDh53KL/zDR5IZvr///nup2zMuOOd+KXx8fIq8P/fcuXPmcldo0qSJ7rjjDm3YsEF5eXlObRe3/dK2vWfPHn333XeqXbu201TwPh87dsyp/sWvs/THa12W17lAWT6Dx48fV0ZGhhYuXFiobyNGjHDq25QpU+Tn56cOHTqoadOmSkhIKPbwdsF7zVXfVxbnqOFW48eP16JFizRx4kRFR0crMDBQNptNQ4YMKfbcrrt06dJF+/bt04oVK/TJJ5/o3//+t5577jktWLBA999/f5nacFVoFCjuF15BkFwpnp6eRZaXFsIhISHlCpmihIeHF3n7VkFZRETEZbV/ocjISJ0/f17Z2dkKCAhQeHi407Yu3n5p287Pz1ebNm307LPPFru9C13q63yhsnwGC3627r33XsXHxxdZp23btpL++ONt9+7dWrlypVJSUvTee+9p3rx5euKJJzR9+nSndX7//Xf5+vq6/OcAJSOoUaz69etr7dq1OnXqlNOo+scffzSXFygucN59913Fx8frn//8p1l27ty5S3qYQ+3ateXj41Pkoevdu3eXqY3g4GCNGDHCvMioS5cumjZtmhnUrhwp5Ofn6+effzZHV5L0008/SZJ5AVDByPXi16OooxVl7Vvt2rXl6+tb5Gvy448/ysPDo1CAXKrmzZtryZIlyszMVGBg4CW10b59e3366afKz893uqBs27Zt8vX1dXr9LtfPP/8sb29v83ax1q1bq1q1atq+fbsGDx5s1jt//rx27tzpVFaUxo0b6+uvv1b37t1d9tlxRTu1a9eWv7+/8vLyirwi/2I1atTQ3Xffrbvvvlvnz5/XwIED9fe//11JSUlOt33t37+/TM9VgGtx6BvF6tOnj/Ly8vTSSy85lT/33HOy2WyKjY01y2rUqFFk+Hp6ehYaLbz44ouXNGL09PRUr169tHz5ch06dMgs/+GHH7R69epS17/4ARl+fn5q0qSJ02HPGjVqSCocnJfqwtfOMAy99NJLql69urp37y7pjz92PD09tXnzZqf15s2bV6itsvbN09NTPXv21IoVK5wOsaenp2vp0qXq1KmT061IlyM6OlqGYWjHjh2X3MagQYOUnp6u999/3yz77bfftGzZMvXr1++S7ocv6tayr7/+Wh9++KF69uxp/kEQGBioHj166M0339SpU6fMum+88YZOnz6tu+66q8TtDB48WL/++qv+9a9/FVp29uxZZWdnl7vvrvgMenp6Ki4uTu+9956+/fbbQssvfH0u/rnw8vJSy5YtZRiGcnJynJZ99dVXuvnmmy+5X7g0jKhRrH79+qlbt27629/+pgMHDqhdu3b65JNPtGLFCk2cONHpopioqCitXbtWzz77rCIiItSwYUN17NhRt99+u9544w0FBgaqZcuWSk1N1dq1axUSEnJJfZo+fbpSUlLUuXNnPfTQQ8rNzTXvAS3tfHPLli0VExNj3g+6fft2vfvuu04XfEVFRUn644lWvXr1kqenp4YMGXJJffX29lZKSori4+PVsWNHrVq1Sh999JH++te/mhekBQYG6q677tKLL74om82mxo0ba+XKlYXObZa3b0899ZTWrFmjTp066aGHHlK1atX08ssvy+FwaNasWZe0P0Xp1KmTQkJCtHbt2kIPwfjvf/+rr7/+WtIfFxLu2rVLTz31lCSpf//+5qHXQYMG6aabbtKIESP0/fffm08my8vLK3TodfPmzeYfNcePH1d2drbZZpcuXdSlSxdJ0t133y0fHx/dfPPNCg0N1ffff6+FCxfK19dXM2fOdGrz73//u26++WZ17dpVY8aM0S+//KJ//vOf6tmzp3r37l3i/v/pT3/SO++8owceeEAbNmzQLbfcory8PP3444965513tHr16kK3fpWm4H3+29/+piFDhqh69erq16+fGeBlNXPmTG3YsEEdO3bU6NGj1bJlS508eVJfffWV1q5dq5MnT0qSevbsqbCwMN1yyy2qU6eOfvjhB7300kvq27ev05G0HTt26OTJk7rjjjvK1Q+4QAVdbQ4LKur2n1OnThmTJk0yIiIijOrVqxtNmzY1Zs+e7XQrj2EYxo8//mh06dLF8PHxMSSZt+L8/vvvxogRI4xatWoZfn5+Rq9evYwff/yx0K1TZb09yzAMY9OmTUZUVJTh5eVlNGrUyFiwYIExderUUm/Peuqpp4wOHToYQUFBho+Pj9G8eXPj73//u3H+/HmzTm5urjF+/Hijdu3ahs1mM9ssuF1q9uzZhfpT3O1ZNWrUMPbt22f07NnT8PX1NerUqWNMnTrVyMvLc1r/+PHjRlxcnOHr62vUrFnTGDt2rPHtt98WarO4vhlG4VtzDMMwvvrqK6NXr16Gn5+f4evra3Tr1s3YsmWLU52C27O+/PJLp/LyvB8TJkwwmjRpUqi84Ba1oqYL98swDOPkyZPGqFGjjJCQEMPX19fo2rVroT4ZhmG+z0VNF+7/888/b3To0MEIDg42qlWrZoSHhxv33nuvsWfPniL34dNPPzVuvvlmw9vb26hdu7aRkJBgZGVllbrvhmEY58+fN55++mmjVatWht1uN2rWrGlERUUZ06dPNzIzM816koyEhIRC61/8OTUMw3jyySeNa665xvDw8HC6Vau4NgqWXfwZSE9PNxISEozIyEijevXqRlhYmNG9e3dj4cKFZp2XX37Z6NKlixESEmLY7XajcePGxiOPPOLUd8MwjClTphj16tUr9LMP97MZxmVesgngqvbzzz+refPmWrVqlXlIH1WLw+FQgwYN9Oijj+rhhx+u6O5cdThHDeCyNGrUSKNGjSp0SBlVx6JFi1S9evUin+kO92NEDQCAhTGiBgDAwghqAAAsjKAGAMDCCGoAACysUj7wJD8/X0eOHJG/vz8PhwcAVDqGYejUqVOKiIgo9F3sF6uUQX3kyBGXPasYAICKcvjwYdWtW7fEOpUyqAsea3f48GGXPbMYAIArJSsrS5GRkYW+RrgolTKoCw53BwQEENQAgEqrLKdvuZgMAAALI6gBALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAsjqAEAsLBK+QjRK6HBox85zR+Y2beCegIAuJqVa0SdnJysG2+8Uf7+/goNDdWAAQO0e/dupzrnzp1TQkKCQkJC5Ofnp7i4OKWnpzvVOXTokPr27StfX1+FhobqkUceUW5u7uXvDQAAVUy5gnrTpk1KSEjQ1q1btWbNGuXk5Khnz57Kzs4260yaNEn//e9/tWzZMm3atElHjhzRwIEDzeV5eXnq27evzp8/ry1btui1117T4sWL9cQTT7hurwAAqCJshmEYl7ry8ePHFRoaqk2bNqlLly7KzMxU7dq1tXTpUg0aNEiS9OOPP6pFixZKTU3VTTfdpFWrVun222/XkSNHVKdOHUnSggULNGXKFB0/flxeXl6lbjcrK0uBgYHKzMx027dncegbAOAu5cmxy7qYLDMzU5IUHBwsSdqxY4dycnLUo0cPs07z5s1Vr149paamSpJSU1PVpk0bM6QlqVevXsrKytJ3331X5HYcDoeysrKcJgAArgaXHNT5+fmaOHGibrnlFrVu3VqSlJaWJi8vLwUFBTnVrVOnjtLS0sw6F4Z0wfKCZUVJTk5WYGCgOUVGRl5qtwEAqFQu+arvhIQEffvtt/rss89c2Z8iJSUlKTEx0ZzPyspyeVhffKgbAAAruKSgHjdunFauXKnNmzerbt26ZnlYWJjOnz+vjIwMp1F1enq6wsLCzDpffPGFU3sFV4UX1LmY3W6X3W6/lK4CAFCplevQt2EYGjdunD744AOtX79eDRs2dFoeFRWl6tWra926dWbZ7t27dejQIUVHR0uSoqOj9c033+jYsWNmnTVr1iggIEAtW7a8nH0BAKDKKdeIOiEhQUuXLtWKFSvk7+9vnlMODAyUj4+PAgMDNWrUKCUmJio4OFgBAQEaP368oqOjddNNN0mSevbsqZYtW+pPf/qTZs2apbS0ND322GNKSEhg1AwAwEXKFdTz58+XJMXExDiVL1q0SMOHD5ckPffcc/Lw8FBcXJwcDod69eqlefPmmXU9PT21cuVKPfjgg4qOjlaNGjUUHx+vGTNmXN6eAABQBV3WfdQVxR33UZd2MdmF91FzjzUA4HJcsfuoAQCAexHUAABYGEENAICFEdQAAFgY30ftBlxsBgBwFUbUAABYGEENAICFcei7jCrrl3ZwGB4AKjdG1AAAWBhBDQCAhRHUAABYGEENAICFcTFZJcAFYQBw9WJEDQCAhRHUAABYGEENAICFEdQAAFgYF5NVAC4OAwCUFUFtAQQ3AKA4HPoGAMDCCGoAACys3EG9efNm9evXTxEREbLZbFq+fLnTcpvNVuQ0e/Zss06DBg0KLZ85c+Zl70xV0eDRj5wmAMDVq9znqLOzs9WuXTuNHDlSAwcOLLT86NGjTvOrVq3SqFGjFBcX51Q+Y8YMjR492pz39/cvb1fgApwfBwBrK3dQx8bGKjY2ttjlYWFhTvMrVqxQt27d1KhRI6dyf3//QnWL43A45HA4zPmsrKxy9BgAgMrLreeo09PT9dFHH2nUqFGFls2cOVMhISG67rrrNHv2bOXm5hbbTnJysgIDA80pMjLSnd2uVDhMDgBVm1tvz3rttdfk7+9f6BD5hAkTdP311ys4OFhbtmxRUlKSjh49qmeffbbIdpKSkpSYmGjOZ2VlEdYAgKuCW4P61Vdf1bBhw+Tt7e1UfmHotm3bVl5eXho7dqySk5Nlt9sLtWO324sst4rSRrKMdAEAl8pth74//fRT7d69W/fff3+pdTt27Kjc3FwdOHDAXd0BAKBScltQv/LKK4qKilK7du1Krbtz5055eHgoNDTUXd0BAKBSKveh79OnT2vv3r3m/P79+7Vz504FBwerXr16kv44h7xs2TL985//LLR+amqqtm3bpm7dusnf31+pqamaNGmS7r33XtWsWfMydgWuUNphem7fAoArq9xBvX37dnXr1s2cLzjfHB8fr8WLF0uS3nrrLRmGoaFDhxZa326366233tK0adPkcDjUsGFDTZo0yem8NQAA+EO5gzomJkaGYZRYZ8yYMRozZkyRy66//npt3bq1vJsFAOCqxLO+AQCwMIIaAAALI6gBALAwtz7wBFcXvuADAFyPETUAABbGiLqK4XGlAFC1ENQolwv/EODQNgC4H4e+AQCwMIIaAAALI6gBALAwghoAAAvjYrJKiCu7AeDqwYgaAAALY0QNy+DJZgBQGCNqAAAsjKAGAMDCCGoAACyMoAYAwMIIagAALIygBgDAwghqAAAsrNxBvXnzZvXr108RERGy2Wxavny50/Lhw4fLZrM5Tb1793aqc/LkSQ0bNkwBAQEKCgrSqFGjdPr06cvaEVhfg0c/cpoAAKUrd1BnZ2erXbt2mjt3brF1evfuraNHj5rTf/7zH6flw4YN03fffac1a9Zo5cqV2rx5s8aMGVP+3gMAUMWV+8lksbGxio2NLbGO3W5XWFhYkct++OEHpaSk6Msvv9QNN9wgSXrxxRfVp08fPfPMM4qIiCi0jsPhkMPhMOezsrLK220AAColt5yj3rhxo0JDQ9WsWTM9+OCDOnHihLksNTVVQUFBZkhLUo8ePeTh4aFt27YV2V5ycrICAwPNKTIy0h3dBgDAclwe1L1799brr7+udevW6emnn9amTZsUGxurvLw8SVJaWppCQ0Od1qlWrZqCg4OVlpZWZJtJSUnKzMw0p8OHD7u62wAAWJLLv5RjyJAh5v/btGmjtm3bqnHjxtq4caO6d+9+SW3a7XbZ7XZXdREAgErD7bdnNWrUSLVq1dLevXslSWFhYTp27JhTndzcXJ08ebLY89oAAFyt3P41l7/88otOnDih8PBwSVJ0dLQyMjK0Y8cORUVFSZLWr1+v/Px8dezY0d3dwRXELVgAcPnKHdSnT582R8eStH//fu3cuVPBwcEKDg7W9OnTFRcXp7CwMO3bt0+TJ09WkyZN1KtXL0lSixYt1Lt3b40ePVoLFixQTk6Oxo0bpyFDhhR5xTcAAFezcgf19u3b1a1bN3M+MTFRkhQfH6/58+dr165deu2115SRkaGIiAj17NlTTz75pNM55iVLlmjcuHHq3r27PDw8FBcXpxdeeMEFu4PSMMoFgMql3EEdExMjwzCKXb569epS2wgODtbSpUvLu2kAAK46bj9HjaqL0TkAuB9BjUrjwj8MDszsW4E9AYArh2/PAgDAwhhRA1eRi09XcGQCsD5G1AAAWBhBDQCAhRHUAABYGEENAICFEdQAAFgYQQ0AgIVxexauSuW9TYmHrQCoKAQ1LKukR5RyPzCAqwVBjasCzyUHUFlxjhoAAAtjRA2UE4fdAVxJBDWqBA5tA6iqCGpAjJIBWBdBDbgYoQ/AlbiYDAAACyOoAQCwsHIH9ebNm9WvXz9FRETIZrNp+fLl5rKcnBxNmTJFbdq0UY0aNRQREaH77rtPR44ccWqjQYMGstlsTtPMmTMve2cAuFaDRz8yJwAVo9znqLOzs9WuXTuNHDlSAwcOdFp25swZffXVV3r88cfVrl07/f7773r44YfVv39/bd++3anujBkzNHr0aHPe39//EncBcD2CCYBVlDuoY2NjFRsbW+SywMBArVmzxqnspZdeUocOHXTo0CHVq1fPLPf391dYWFh5Nw8AwFXF7eeoMzMzZbPZFBQU5FQ+c+ZMhYSE6LrrrtPs2bOVm5tbbBsOh0NZWVlOEwAAVwO33p517tw5TZkyRUOHDlVAQIBZPmHCBF1//fUKDg7Wli1blJSUpKNHj+rZZ58tsp3k5GRNnz7dnV0FAMCS3BbUOTk5Gjx4sAzD0Pz5852WJSYmmv9v27atvLy8NHbsWCUnJ8tutxdqKykpyWmdrKwsRUZGuqvrAABYhluCuiCkDx48qPXr1zuNpovSsWNH5ebm6sCBA2rWrFmh5Xa7vcgABwCgqnN5UBeE9J49e7RhwwaFhISUus7OnTvl4eGh0NBQV3cHcDuuEAfgTuUO6tOnT2vv3r3m/P79+7Vz504FBwcrPDxcgwYN0ldffaWVK1cqLy9PaWlpkqTg4GB5eXkpNTVV27ZtU7du3eTv76/U1FRNmjRJ9957r2rWrOm6PQMAoAood1Bv375d3bp1M+cLzh3Hx8dr2rRp+vDDDyVJ7du3d1pvw4YNiomJkd1u11tvvaVp06bJ4XCoYcOGmjRpktM5aABlw3PFgaqv3EEdExMjwzCKXV7SMkm6/vrrtXXr1vJuFgCAqxLP+gYAwML4mkvgKsahc8D6GFEDAGBhjKhRYbitCQBKR1ADFYzDzwBKQlADVYg7Q58/KICKwTlqAAAsjKAGAMDCCGoAACyMc9RAFVbeK+vLU59z1sCVQVADbkagAbgcHPoGAMDCCGoAACyMQ9/AFebO88YAqh5G1AAAWBgjagCWxEV4wB8YUQMAYGEENQAAFkZQAwBgYQQ1AAAWRlADAGBh5Q7qzZs3q1+/foqIiJDNZtPy5cudlhuGoSeeeELh4eHy8fFRjx49tGfPHqc6J0+e1LBhwxQQEKCgoCCNGjVKp0+fvqwdAQCgKip3UGdnZ6tdu3aaO3dukctnzZqlF154QQsWLNC2bdtUo0YN9erVS+fOnTPrDBs2TN99953WrFmjlStXavPmzRozZsyl7wUAAFVUue+jjo2NVWxsbJHLDMPQnDlz9Nhjj+mOO+6QJL3++uuqU6eOli9friFDhuiHH35QSkqKvvzyS91www2SpBdffFF9+vTRM888o4iIiELtOhwOORwOcz4rK6u83QYAoFJy6QNP9u/fr7S0NPXo0cMsCwwMVMeOHZWamqohQ4YoNTVVQUFBZkhLUo8ePeTh4aFt27bpzjvvLNRucnKypk+f7squApbFI0MBXMilF5OlpaVJkurUqeNUXqdOHXNZWlqaQkNDnZZXq1ZNwcHBZp2LJSUlKTMz05wOHz7sym4DAGBZleIRona7XXa7vaK7AQDAFefSoA4LC5MkpaenKzw83CxPT09X+/btzTrHjh1zWi83N1cnT5401wdQ+fBsbsA9XHrou2HDhgoLC9O6devMsqysLG3btk3R0dGSpOjoaGVkZGjHjh1mnfXr1ys/P18dO3Z0ZXcAAKj0yj2iPn36tPbu3WvO79+/Xzt37lRwcLDq1auniRMn6qmnnlLTpk3VsGFDPf7444qIiNCAAQMkSS1atFDv3r01evRoLViwQDk5ORo3bpyGDBlS5BXfAABczcod1Nu3b1e3bt3M+cTERElSfHy8Fi9erMmTJys7O1tjxoxRRkaGOnXqpJSUFHl7e5vrLFmyROPGjVP37t3l4eGhuLg4vfDCCy7YHQAAqhabYRhGRXeivLKyshQYGKjMzEwFBAS4pE1uiQFc63LPUXPOG1VZeXKMZ30DAGBhBDUAABZWKe6jBlD5lXZ6iUPbQNEYUQMAYGEENQAAFkZQAwBgYQQ1AAAWRlADAGBhBDUAABZGUAMAYGHcRw3AEniML1A0RtQAAFgYI2oAlQ5f2IGrCSNqAAAsjBE1ALfgnDPgGoyoAQCwMEbUACqFkkbonLNGVcaIGgAACyOoAQCwMIIaAAALc3lQN2jQQDabrdCUkJAgSYqJiSm07IEHHnB1NwAAqBJcfjHZl19+qby8PHP+22+/1W233aa77rrLLBs9erRmzJhhzvv6+rq6GwAAVAkuD+ratWs7zc+cOVONGzdW165dzTJfX1+FhYW5etMAAFQ5br096/z583rzzTeVmJgom81mli9ZskRvvvmmwsLC1K9fPz3++OMljqodDoccDoc5n5WV5c5uA7jKcbsXrMStQb18+XJlZGRo+PDhZtk999yj+vXrKyIiQrt27dKUKVO0e/duvf/++8W2k5ycrOnTp7uzqwAAWJJbg/qVV15RbGysIiIizLIxY8aY/2/Tpo3Cw8PVvXt37du3T40bNy6ynaSkJCUmJprzWVlZioyMdF/HAQCwCLcF9cGDB7V27doSR8qS1LFjR0nS3r17iw1qu90uu93u8j4CAGB1bruPetGiRQoNDVXfviWf29m5c6ckKTw83F1dAQCg0nLLiDo/P1+LFi1SfHy8qlX7v03s27dPS5cuVZ8+fRQSEqJdu3Zp0qRJ6tKli9q2beuOrgAAUKm5JajXrl2rQ4cOaeTIkU7lXl5eWrt2rebMmaPs7GxFRkYqLi5Ojz32mDu6AQBApeeWoO7Zs6cMwyhUHhkZqU2bNrljkwAAVEk86xsAAAsjqAEAsDCCGgAAC3PrA08AoDK4+JGhgJUwogYAwMIIagAALIygBgDAwjhHDaDK42srUZkxogYAwMIYUQNAOTA6x5XGiBoAAAtjRA0AbsQIHJeLoAZw1SnvA054IAoqEkENoMohWFGVcI4aAAALI6gBALAwDn0DQCXBhWlXJ0bUAABYGCNqALgMrh7lMmrGxRhRAwBgYQQ1AAAW5vKgnjZtmmw2m9PUvHlzc/m5c+eUkJCgkJAQ+fn5KS4uTunp6a7uBgBUCQ0e/ciccHVyyznqVq1aae3atf+3kWr/t5lJkybpo48+0rJlyxQYGKhx48Zp4MCB+vzzz93RFQC4oghUuJpbgrpatWoKCwsrVJ6ZmalXXnlFS5cu1a233ipJWrRokVq0aKGtW7fqpptuKrI9h8Mhh8NhzmdlZbmj2wDgdlwshvJyS1Dv2bNHERER8vb2VnR0tJKTk1WvXj3t2LFDOTk56tGjh1m3efPmqlevnlJTU4sN6uTkZE2fPt0dXQWASqu8oc8fCZWTy89Rd+zYUYsXL1ZKSormz5+v/fv3q3Pnzjp16pTS0tLk5eWloKAgp3Xq1KmjtLS0YttMSkpSZmamOR0+fNjV3QYAwJJcPqKOjY01/9+2bVt17NhR9evX1zvvvCMfH59LatNut8tut7uqiwAAVBpuvz0rKChI1157rfbu3auwsDCdP39eGRkZTnXS09OLPKcNAMDVzu1Bffr0ae3bt0/h4eGKiopS9erVtW7dOnP57t27dejQIUVHR7u7KwAAVDouP/T9l7/8Rf369VP9+vV15MgRTZ06VZ6enho6dKgCAwM1atQoJSYmKjg4WAEBARo/fryio6OLvZAMAICrmcuD+pdfftHQoUN14sQJ1a5dW506ddLWrVtVu3ZtSdJzzz0nDw8PxcXFyeFwqFevXpo3b56ruwEAQJXg8qB+6623Slzu7e2tuXPnau7cua7eNABUOjwgBaXhWd8AAFgYQQ0AgIUR1AAAWJhbHiEKALjyeERo1cSIGgAACyOoAQCwMIIaAAALI6gBALAwghoAAAvjqm8AuEpxlXjlwIgaAAALY0QNAHA5Ruuuw4gaAAALY0QNALhsfAuY+xDUAFBFEZ5VA0ENAJDk3vPKnLO+dJyjBgDAwhhRAwBKxYi44hDUAIBy4/z3lcOhbwAALIwRNQCgSIyarcHlI+rk5GTdeOON8vf3V2hoqAYMGKDdu3c71YmJiZHNZnOaHnjgAVd3BQCASs/lQb1p0yYlJCRo69atWrNmjXJyctSzZ09lZ2c71Rs9erSOHj1qTrNmzXJ1VwAAqPRcfug7JSXFaX7x4sUKDQ3Vjh071KVLF7Pc19dXYWFhZWrT4XDI4XCY81lZWa7pLAAAFuf2i8kyMzMlScHBwU7lS5YsUa1atdS6dWslJSXpzJkzxbaRnJyswMBAc4qMjHRrnwEAsAqbYRiGuxrPz89X//79lZGRoc8++8wsX7hwoerXr6+IiAjt2rVLU6ZMUYcOHfT+++8X2U5RI+rIyEhlZmYqICDAJX3logkAuHKu9vuws7KyFBgYWKYcc+tV3wkJCfr222+dQlqSxowZY/6/TZs2Cg8PV/fu3bVv3z41bty4UDt2u112u92dXQUAwJLcduh73LhxWrlypTZs2KC6deuWWLdjx46SpL1797qrOwAAVEouH1EbhqHx48frgw8+0MaNG9WwYcNS19m5c6ckKTw83NXdAQBUAjyitHguD+qEhAQtXbpUK1askL+/v9LS0iRJgYGB8vHx0b59+7R06VL16dNHISEh2rVrlyZNmqQuXbqobdu2ru4OAACVmsuDev78+ZL+eKjJhRYtWqThw4fLy8tLa9eu1Zw5c5Sdna3IyEjFxcXpsccec3VXAACo9Nxy6LskkZGR2rRpk6s3CwBAlcSzvgEAlna1n7/m27MAALAwghoAAAsjqAEAsDCCGgAACyOoAQCwMK76BgBUKVXtKnGCGgBwxZX2jYUlLb/avu2QQ98AAFgYQQ0AgIUR1AAAWBhBDQCAhXExGQDgqlLZrgpnRA0AgIUR1AAAWBiHvgEAVVplv++aETUAABbGiBoAcFWz+sVljKgBALAwghoAAAursKCeO3euGjRoIG9vb3Xs2FFffPFFRXUFAADLqpCgfvvtt5WYmKipU6fqq6++Urt27dSrVy8dO3asIroDAIBlVUhQP/vssxo9erRGjBihli1basGCBfL19dWrr75aEd0BAMCyrvhV3+fPn9eOHTuUlJRklnl4eKhHjx5KTU0tch2HwyGHw2HOZ2ZmSpKysrJc1q98xxmXtQUAqLxcmS2lbcMwjFLrXvGg/u2335SXl6c6deo4ldepU0c//vhjkeskJydr+vTphcojIyPd0kcAwNUrcM6V29apU6cUGBhYYp1KcR91UlKSEhMTzfn8/HydPHlSISEhstlsl91+VlaWIiMjdfjwYQUEBFx2exWBfbCOqrAf7IN1VIX9YB8KMwxDp06dUkRERKl1r3hQ16pVS56enkpPT3cqT09PV1hYWJHr2O122e12p7KgoCCX9y0gIKDSfogKsA/WURX2g32wjqqwH+yDs9JG0gWu+MVkXl5eioqK0rp168yy/Px8rVu3TtHR0Ve6OwAAWFqFHPpOTExUfHy8brjhBnXo0EFz5sxRdna2RowYURHdAQDAsiokqO+++24dP35cTzzxhNLS0tS+fXulpKQUusDsSrHb7Zo6dWqhw+uVCftgHVVhP9gH66gK+8E+XB6bUZZrwwEAQIXgWd8AAFgYQQ0AgIUR1AAAWBhBDQCAhV01QV3er9VctmyZmjdvLm9vb7Vp00Yff/zxFeppYcnJybrxxhvl7++v0NBQDRgwQLt37y5xncWLF8tmszlN3t7eV6jHRZs2bVqhPjVv3rzEdaz0PkhSgwYNCu2DzWZTQkJCkfWt8D5s3rxZ/fr1U0REhGw2m5YvX+603DAMPfHEEwoPD5ePj4969OihPXv2lNrulf6q2pL2IycnR1OmTFGbNm1Uo0YNRURE6L77/l979x8S5R3HAfwtdWezNW22sh9qKstJK9eCRNuINPxRkLYoFVnJ2i9no6KgwQqLgm0Z/VGESagVQVFQi2XUdHnHlj+XzqlrruIwJG9BcM66nM777I/hrct77lfzuWf6foHg3X2eL5/v83m+z8c7xe8G3L9/3+WYvlyTozUHAMjPzx+RT3p6uttx1ayFuzk4Wx8BAQEoLi5WHFPtOnhyT+3v70dhYSFCQ0Px4osvYu3atSP+UdezfF1L7oyLRu3ttpq1tbXIzc3Fpk2b0NLSgqysLGRlZaG9vV3lzP9hNBpRWFiI+vp6VFVVYXBwEKmpqXj8+LHL41566SX09PTYv7q6ulTKWNn8+fMdcvrhhx8UY7VWBwBoampyyL+qqgoAsG7dOsVj/F2Hx48fIz4+HkePHnX6+oEDB3D48GEcO3YMDQ0NmDx5MtLS0tDf3684pj+2qnU1D6vViubmZuzevRvNzc24cOECOjs7sXr1arfjenNNPi93tQCA9PR0h3zOnDnjcky1a+FuDk/n3tPTg/LycgQEBGDt2rUux1WzDp7cU7dt24ZvvvkG58+fh9FoxP379/HOO++4HNeXteQRGQeWLFkihYWF9sdDQ0Mya9Ys+eKLL5zGr1+/XlatWuXwXEJCgnz00UejmqenHjx4IADEaDQqxlRUVEhwcLB6SXmgqKhI4uPjPY7Xeh1ERLZs2SIxMTFis9mcvq61OgCQixcv2h/bbDYJCwuT4uJi+3MWi0UCAwPlzJkziuN4u6b+a8/Ow5nGxkYBIF1dXYox3l6T/yVnc9i4caNkZmZ6NY4/a+FJHTIzMyU5OdlljD/rIDLynmqxWESn08n58+ftMbdu3RIAUldX53QMX9eSJ8b8O+rhbTVXrFhhf87dtpp1dXUO8QCQlpamGK+24W0+X375ZZdxjx49QmRkJMLDw5GZmYmOjg410nPp9u3bmDVrFqKjo5GXl4d79+4pxmq9DgMDAzh9+jTee+89l5vDaLEOw0wmE8xms8N5Dg4ORkJCguJ59mVN+UNvby8CAgLc7gvgzTWpBoPBgOnTpyM2NhYFBQV4+PChYqzWa/H777+jsrISmzZtchvrzzo8e0+9efMmBgcHHc7ra6+9hoiICMXz6sta8tSYb9SuttU0m81OjzGbzV7Fq8lms2Hr1q1YunQpXn/9dcW42NhYlJeX49KlSzh9+jRsNhuSkpLQ3d2tYraOEhIScOLECVy9ehUlJSUwmUx4++230dfX5zRey3UAgK+//hoWiwX5+fmKMVqsw9OGz6U359mXNaW2/v5+7Ny5E7m5uS43UPD2mhxt6enpOHXqFL777jt89dVXMBqNyMjIwNDQkNN4rdfi5MmTmDJlituPjP1ZB2f3VLPZDL1eP+KHPHd9YzjG02M89b/Y5pL+VVhYiPb2dre/v0lMTHTY5CQpKQlxcXEoLS3Fvn37RjtNpzIyMuzfL1y4EAkJCYiMjMS5c+c8+olba8rKypCRkeFymzot1mGsGxwcxPr16yEiKCkpcRmrtWsyJyfH/v2CBQuwcOFCxMTEwGAwICUlRfV8nld5eTny8vLc/gGlP+vg6T3Vn8b8O2pfttUMCwvzKl4tmzdvxuXLl1FTU4M5c+Z4daxOp8OiRYtw586dUcrOeyEhIZg3b55iTlqtAwB0dXWhuroa77//vlfHaa0Ow+fSm/Psy5pSy3CT7urqQlVVldfbEbq7JtUWHR2NadOmKeaj5Vp8//336Ozs9HqNAOrVQemeGhYWhoGBAVgsFod4d31jOMbTYzw15hu1L9tqJiYmOsQDQFVVld+24RQRbN68GRcvXsT169cRFRXl9RhDQ0Noa2vDzJkzRyFD3zx69Ah3795VzElrdXhaRUUFpk+fjlWrVnl1nNbqEBUVhbCwMIfz/Mcff6ChoUHxPGt1q9rhJn379m1UV1cjNDTU6zHcXZNq6+7uxsOHDxXz0WotgH8+cVq8eDHi4+O9Pna06+Dunrp48WLodDqH89rZ2Yl79+4pnldf1pI3CY95Z8+elcDAQDlx4oT88ssv8uGHH0pISIiYzWYREXn33Xfls88+s8ffuHFDJk6cKAcPHpRbt25JUVGR6HQ6aWtr80v+BQUFEhwcLAaDQXp6euxfVqvVHvPsHPbu3SvXrl2Tu3fvys2bNyUnJ0cmTZokHR0d/piCiIhs375dDAaDmEwmuXHjhqxYsUKmTZsmDx48EBHt12HY0NCQREREyM6dO0e8psU69PX1SUtLi7S0tAgAOXTokLS0tNj/GvrLL7+UkJAQuXTpkvz888+SmZkpUVFR8uTJE/sYycnJcuTIEftjd2tK7XkMDAzI6tWrZc6cOfLTTz85rJM///xTcR7urkk159DX1yc7duyQuro6MZlMUl1dLW+++aa8+uqr0t/frzgHtWvh7noSEent7ZWgoCApKSlxOoa/6+DJPfXjjz+WiIgIuX79uvz444+SmJgoiYmJDuPExsbKhQsX7I89WUu+GBeNWkTkyJEjEhERIXq9XpYsWSL19fX215YtWyYbN250iD937pzMmzdP9Hq9zJ8/XyorK1XO+F8AnH5VVFTYY56dw9atW+3znTFjhqxcuVKam5vVT/4p2dnZMnPmTNHr9TJ79mzJzs6WO3fu2F/Xeh2GXbt2TQBIZ2fniNe0WIeamhqn189wnjabTXbv3i0zZsyQwMBASUlJGTG3yMhIKSoqcnjO1ZpSex4mk0lxndTU1CjOw901qeYcrFarpKamyiuvvCI6nU4iIyPlgw8+GNFw/V0Ld9eTiEhpaam88MILYrFYnI7h7zp4ck998uSJfPLJJzJ16lQJCgqSNWvWSE9Pz4hxnj7Gk7XkC25zSUREpGFj/nfURERE/2ds1ERERBrGRk1ERKRhbNREREQaxkZNRESkYWzUREREGsZGTUREpGFs1ERERBrGRk1ERKRhbNREBOCfDUOSkpJG7B3c29uL8PBwfP75537KjGh8478QJSK73377DW+88QaOHz+OvLw8AMCGDRvQ2tqKpqYm6PV6P2dINP6wURORg8OHD2PPnj3o6OhAY2Mj1q1bh6amJp+2KySi58dGTUQORATJycmYMGEC2tra8Omnn2LXrl3+Toto3GKjJqIRfv31V8TFxWHBggVobm7GxIkT/Z0S0bjFPyYjohHKy8sRFBQEk8mE7u5uf6dDNK7xHTUROaitrcWyZcvw7bffYv/+/QCA6upqBAQE+DkzovGJ76iJyM5qtSI/Px8FBQVYvnw5ysrK0NjYiGPHjvk7NaJxi++oichuy5YtuHLlClpbWxEUFAQAKC0txY4dO9DW1oa5c+f6N0GicYiNmogAAEajESkpKTAYDHjrrbccXktLS8Nff/3Fj8CJ/ICNmoiISMP4O2oiIiINY6MmIiLSMDZqIiIiDWOjJiIi0jA2aiIiIg1joyYiItIwNmoiIiINY6MmIiLSMDZqIiIiDWOjJiIi0jA2aiIiIg37G4hUH/370iAIAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 500x1000 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plotting the distributions again\n",
    "singleToy(nevents_sig=nevents_sig, nevents_bgd=nevents_bgd)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "96d01bd2",
   "metadata": {},
   "source": [
    "### Generating the updated datasets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c4b95e9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "nevents_sig = 150\n",
    "nevents_bgd = 10000\n",
    "\n",
    "sig_fraction = nevents_sig/(nevents_bgd+nevents_sig)\n",
    "    \n",
    "# Create the pdf\n",
    "pdf = SignalWithBackground(mean, sigma, sig_fraction, intercept, slope, XMIN, XMAX)\n",
    "for i in range (nevents_sig + nevents_bgd): pdf.next()\n",
    "data = pdf.mass\n",
    "sig_data = pdf.mass_sig\n",
    "bgd_data = pdf.mass_bgd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d542bc71",
   "metadata": {},
   "source": [
    "### Defining the probability density functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "8cb70683",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Generate the normalised exponential PDF\n",
    "def Linear_PDF(x,intercept,slope,x_low,x_high):\n",
    "    linear_dist = lambda x: intercept + slope*(x)\n",
    "    area, _ = integrate.quad(linear_dist,x_low,x_high) # Calculate the normalisation constant\n",
    "    return (1/area)*(intercept + slope*(x))  # Return the normalised exponential PDF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "18260402",
   "metadata": {},
   "outputs": [],
   "source": [
    "def gauss_PDF(x):\n",
    "    return scipy.stats.norm.pdf(x,mean,sigma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "5ee751ea",
   "metadata": {},
   "outputs": [],
   "source": [
    "def combined_PDF(x,intercept,slope,F,x_low,x_high):\n",
    "    PDF1 = Linear_PDF(x,intercept,slope,x_low,x_high)\n",
    "    PDF2 = gauss_PDF(x)\n",
    "    \n",
    "    return ((1-F) * PDF1) + (F * PDF2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "3c5b85db",
   "metadata": {},
   "outputs": [],
   "source": [
    "bin_width = XMAX/NBINS"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cba1fab",
   "metadata": {},
   "source": [
    "### Fitting assuming no signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "586a6f5d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10150"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(np.histogram(data,bins = np.linspace(0,20,101))[0]).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "b50dad39",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "10149.999999999998"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(nevents_bgd+nevents_sig)*bin_width*Linear_PDF(np.linspace(0.1,19.9,100),intercept,slope,XMIN,XMAX).sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "2e5b952f",
   "metadata": {},
   "outputs": [],
   "source": [
    "#create a function for the chi squared linear fit\n",
    "def chisq_linear(slope,intercept):\n",
    "    n_observed = np.histogram(data,bins = np.linspace(0,20,101))[0]\n",
    "    n_expected = (nevents_bgd+nevents_sig)*bin_width*Linear_PDF(np.linspace(0.1,19.9,100),intercept,slope,XMIN,XMAX)\n",
    "    chi_sq_sum = 0\n",
    "    for i in range(0,NBINS):\n",
    "        if n_expected[i] < 0:\n",
    "            n_expected[i] = 1e-3\n",
    "        chi_sq_sum += (n_expected[i]-n_observed[i])+(n_observed[i]*np.log(n_observed[i]/n_expected[i]))\n",
    "    chi_sq = 2*chi_sq_sum\n",
    "    return chi_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "01104132",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting the range of our mass data for use in normalisation\n",
    "x_low = 0 \n",
    "x_high = 20.\n",
    "\n",
    "# Set the minimisation parameters\n",
    "m = Minuit(chisq_linear, intercept=19, slope=-1.1)\n",
    "m.limits[\"intercept\"] = (15,25) \n",
    "m.limits[\"slope\"] = (-1.5,-0.5)\n",
    "\n",
    "# Set the errors on the initial start values\n",
    "m.errors=[0.1,0.1]\n",
    "\n",
    "# Set the increase in Chi squared for the 1-SD error calculation\n",
    "m.errordef=1.0\n",
    "\n",
    "# Run the minimisation\n",
    "m = m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "4fd5ff86",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept =   20.3489  +/-  0.3937\n",
      "\n",
      "\n",
      "slope =   -1.0197  +/-  0.0197\n"
     ]
    }
   ],
   "source": [
    "# Printing the final results and their error\n",
    "fstr = \"{0:2.4f}\"\n",
    "val = m.values[\"intercept\",\"slope\"]\n",
    "err = m.errors[\"intercept\",\"slope\"]\n",
    "print(\"{:2s}\".format(\"intercept\"), \"=  \", fstr.format(val[0]), \" +/- \", fstr.format(err[0]))\n",
    "print(\"\\n\")\n",
    "print(\"{:2s}\".format(\"slope\"), \"=  \", fstr.format(val[1]), \" +/- \", fstr.format(err[1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "1db75c3c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual intercept: 20.0\n",
      "Actual slope: -1.0\n"
     ]
    }
   ],
   "source": [
    "print(\"Actual intercept:\", intercept)\n",
    "print(\"Actual slope:\", slope)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2374f6d",
   "metadata": {},
   "source": [
    "All of the fitted values are correct to within the error margin, showing that a good fit was made"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "d5bed42b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collecting the final fit values\n",
    "h0_fit_values = m.values[\"intercept\",\"slope\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d1a02973",
   "metadata": {},
   "source": [
    "### Fitting assuming a gaussian signal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "b14405d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create function for chi squared combined fit\n",
    "def chisq_combined(slope,intercept,F):\n",
    "    n_observed = np.histogram(data,bins = np.linspace(0,20,101))[0]\n",
    "    n_expected = (nevents_bgd+nevents_sig)*bin_width*combined_PDF(np.linspace(0.1,19.9,100),intercept,slope,F,XMIN,XMAX)\n",
    "    chi_sq_sum = 0\n",
    "    for i in range(0,NBINS):\n",
    "        if n_expected[i] <= 0:\n",
    "            n_expected[i] = 1e-3\n",
    "        chi_sq_sum += (n_expected[i]-n_observed[i])+(n_observed[i]*np.log(n_observed[i]/n_expected[i]))\n",
    "    chi_sq = 2*chi_sq_sum\n",
    "    return chi_sq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "dd3b4ac5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Setting the range of our mass data for use in normalisation\n",
    "x_low = 0 \n",
    "x_high = 20.\n",
    "\n",
    "# Set the minimisation parameters\n",
    "m = Minuit(chisq_combined, intercept=19, slope=-1.1,F=0.01)\n",
    "m.limits[\"intercept\"] = (15,25) \n",
    "m.limits[\"slope\"] = (-1.5,-0.5)\n",
    "m.limits[\"F\"] = (0.0,0.1)\n",
    "\n",
    "# Set the errors on the initial start values\n",
    "m.errors=[0.1,0.1,0.1]\n",
    "\n",
    "# Set the increase in Chi squared for the 1-SD error calculation\n",
    "m.errordef=1.0\n",
    "\n",
    "# Run the minimisation\n",
    "m = m.migrad()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "4a2654e1",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "intercept =   20.4529  +/-  1.4304\n",
      "\n",
      "\n",
      "slope =   -1.0160  +/-  0.0718\n",
      "\n",
      "\n",
      "signal fraction =   0.0159  +/-  0.0033\n"
     ]
    }
   ],
   "source": [
    "# Printing the final results and their error\n",
    "fstr = \"{0:2.4f}\"\n",
    "val = m.values[\"intercept\",\"slope\",\"F\"]\n",
    "err = m.errors[\"intercept\",\"slope\",\"F\"]\n",
    "print(\"{:2s}\".format(\"intercept\"), \"=  \", fstr.format(val[0]), \" +/- \", fstr.format(err[0]))\n",
    "print(\"\\n\")\n",
    "print(\"{:2s}\".format(\"slope\"), \"=  \", fstr.format(val[1]), \" +/- \", fstr.format(err[1]))\n",
    "print(\"\\n\")\n",
    "print(\"{:2s}\".format(\"signal fraction\"), \"=  \", fstr.format(val[2]), \" +/- \", fstr.format(err[2]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "8d75c581",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Actual intercept: 20.0\n",
      "Actual slope: -1.0\n",
      "Actual signal fraction: 0.013694581280788177\n"
     ]
    }
   ],
   "source": [
    "print(\"Actual intercept:\", intercept)\n",
    "print(\"Actual slope:\", slope)\n",
    "print(\"Actual signal fraction:\", (len(sig_data)/(len(sig_data)+len(bgd_data))))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7f7b01b",
   "metadata": {},
   "source": [
    "All of the fitted values are correct to within the error margin, showing that a good fit was made"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "b3dfb35a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collecting the final fit values\n",
    "h1_fit_values = m.values[\"intercept\",\"slope\",\"F\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce48c7ce",
   "metadata": {},
   "source": [
    "### Calculating the chi-squared difference between the two hypotheses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "36d5799b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Using the combined \"nested\" hypothesis to determine chi-squared for h0 and h1\n",
    "h0_chi_squared = chisq_combined(intercept=h0_fit_values[0], slope=h0_fit_values[1], F=0)\n",
    "h1_chi_squared = chisq_combined(intercept=h1_fit_values[0], slope=h1_fit_values[1], F=h1_fit_values[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "bf35ba92",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared for h0:  114.91716261389219\n"
     ]
    }
   ],
   "source": [
    "print(\"Chi-squared for h0: \",h0_chi_squared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "59693436",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared for h1:  89.57851539337628\n"
     ]
    }
   ],
   "source": [
    "print(\"Chi-squared for h1: \",h1_chi_squared)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "9499c323",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Chi-squared difference between the two hypotheses:  25.33864722051591\n"
     ]
    }
   ],
   "source": [
    "chi_squared_diff = h0_chi_squared - h1_chi_squared\n",
    "print(\"Chi-squared difference between the two hypotheses: \",chi_squared_diff)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "957f619d",
   "metadata": {},
   "source": [
    "### Finding the p-value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "d1c2836b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The p-value is:  4.809748138656644e-07\n"
     ]
    }
   ],
   "source": [
    "pvalue = 1 - scipy.stats.chi2.cdf(chi_squared_diff,1) \n",
    "print(\"The p-value is: \",pvalue)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dde68ce",
   "metadata": {},
   "source": [
    "### Finding the Z-score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "32febf63",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The Z-score is:  5.033750810370506\n"
     ]
    }
   ],
   "source": [
    "n_sigmas = erfinv(1 - pvalue )* np.sqrt(2)\n",
    "print(\"The Z-score is: \",n_sigmas)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a115d77",
   "metadata": {},
   "source": [
    "These results show that for this signal to be a fluctuation, this fluction would occur over 4 standard deviations away from the average results that would be expected if there was no signal. This means the hypothesis that there is a signal is strongly favoured as such a fluctuation is highly improbable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "641f5646",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}