{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Table Lenght Measurements\n",
    "\n",
    "### Authors: \n",
    "- Troels C. Petersen (Niels Bohr Institute)\n",
    "- Christian Michelsen (Niels Bohr Institute)\n",
    "\n",
    "\n",
    "### Date:    \n",
    "- 22-11-2018 (latest update)\n",
    "\n",
    "***\n",
    "\n",
    "Python program for analysing measurements of the length of the lecture table in Auditorium A at NBI.  \n",
    "There are two measurements each with estimated error of the table length:\n",
    "1. Measurement with a 30cm ruler.\n",
    "2. Measurement with a 2m folding ruler.\n",
    " \n",
    "Each person was asked not only to state the measurement, but also their estimated uncertainty. None of the persons could see others measurements in order to get the largest degree of independence. Also, the 30cm ruler measurement was asked to be done first.\n",
    "\n",
    "***"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First, import the modules you want to use:\n",
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "from probfit import BinnedLH, Chi2Regression, Extended, UnbinnedLH # Helper tool for fitting\n",
    "import sys\n",
    "from scipy import stats\n",
    "from scipy.special import erfc\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And add out our custom functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../External_Functions')\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax # useful functions to print fit results on figure"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Options for the program: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "blinded = True\n",
    "apply_chauvenets_criterion = True\n",
    "save_plots = False\n",
    "\n",
    "r = np.random          # Random numbers\n",
    "r.seed(42)             # We set the numbers to be random, but the same for each run"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Functions:\n",
    "\n",
    "Function to find the bin centers and counts in a range from `xmin` to `xmax` of a histogram `hist`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "def get_bincenter_and_counts_in_range(hist, xmin=None, xmax=None):\n",
    "    \n",
    "    if xmin is None:\n",
    "        xmin = np.min(hist)\n",
    "    if xmax is None:\n",
    "        xmax = np.max(hist)\n",
    "    \n",
    "    counts, bin_edges, _ = hist\n",
    "    bin_centers = 0.5*(bin_edges[1:] + bin_edges[:-1])\n",
    "    mask1 = (xmin < bin_centers) & (bin_centers <= xmax) \n",
    "    mask2 = counts > 0\n",
    "    mask_final = mask1 & mask2\n",
    "    return bin_centers[mask_final], counts[mask_final], np.sqrt(counts[mask_final])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Function to calculate the $\\chi^2$-value:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_chi2(function, x_values, y_values, sy_values, *fitparameters):\n",
    "    # traditional loop-version\n",
    "    chi2_val = 0\n",
    "    entries = 0\n",
    "    for x, y, sy in zip(x_values, y_values, sy_values):\n",
    "        if y > 0:\n",
    "            f = function(x, *fitparameters) # calc the model value\n",
    "            residual  = ( y-f ) / sy  # find the uncertainty-weighted residual\n",
    "            chi2_val += residual**2  # the chi2-value is the squared residual\n",
    "            entries += 1 # count the bin as non-empty since sy>0 (and thus y>0)\n",
    "    \n",
    "    # numpy version\n",
    "    mask = (y_values>0)\n",
    "    yhat = function(x_values, *fitparameters)\n",
    "    chi2_val = np.sum( (y_values[mask]-yhat[mask])**2/sy_values[mask]**2)\n",
    "    entries = sum(mask)\n",
    "            \n",
    "    return chi2_val, entries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A Gaussian PDF (unit integral):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define your PDF / model \n",
    "def gauss_pdf(x, mu, sigma):\n",
    "    \"\"\"Normalized Gaussian\"\"\"\n",
    "    return 1 / np.sqrt(2 * np.pi) / sigma * np.exp(-(x - mu) ** 2 / 2. / sigma ** 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extended Gaussian PDF (including normalisation, which enables fit to histograms):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [],
   "source": [
    "def gauss_extended(x, N, mu, sigma) :\n",
    "    \"\"\"Non-normalized Gaussian\"\"\"\n",
    "    return N * gauss_pdf(x, mu, sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initial data analysis:\n",
    "\n",
    "Before we even look at the look at the data, we decide whether or not we want to blind the analysis by adding a constat to all measurements:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "if blinded:\n",
    "    blinding = r.normal(0, 0.1)      # I add a constant (Gaussian with +-10cm) to remain \"blind\"\n",
    "else:\n",
    "    blinding = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define what datafiles we want to look at. Extend it to suit your analysis: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "infiles = [\"data_TableMeasurements2018.txt\",\n",
    "           \"data_TableMeasurements2017.txt\",\n",
    "           \"data_TableMeasurements2016.txt\",\n",
    "           \"data_TableMeasurements2015.txt\",\n",
    "           \"data_TableMeasurements2014.txt\",\n",
    "           \"data_TableMeasurements2013.txt\",\n",
    "           \"data_TableMeasurements2012.txt\",\n",
    "           \"data_TableMeasurements2011.txt\",\n",
    "           \"data_TableMeasurements2010.txt\",\n",
    "           \"data_TableMeasurements2009.txt\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We read in all the data from the `infiles` files:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "Read all 10 file(s) which included 495 measurements. \n",
      "\n"
     ]
    }
   ],
   "source": [
    "L30cm = np.array([])\n",
    "eL30cm = np.array([])\n",
    "L2m = np.array([])\n",
    "eL2m = np.array([])\n",
    "\n",
    "# Loop over files and open them\n",
    "for infile in infiles:\n",
    "    \n",
    "    tmp_L30cm, tmp_eL30cm, tmp_L2m, tmp_eL2m = np.loadtxt(infile, skiprows=2, unpack=True)\n",
    "    \n",
    "    L30cm = np.append(L30cm, tmp_L30cm + blinding)\n",
    "    eL30cm = np.append(eL30cm, tmp_eL30cm)\n",
    "    L2m = np.append(L2m, tmp_L2m + blinding)\n",
    "    eL2m = np.append(eL2m, tmp_eL2m)\n",
    "    \n",
    "N_read  = len(L30cm)       # Number of measurements read in total\n",
    "\n",
    "\n",
    "print(f\"\\n\\nRead all {len(infiles)} file(s) which included {N_read} measurements. \\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAAI4CAYAAAD56sN/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XmYXnV99/H3hwQI+xr2JWzK4gICirVaRG0RF6wLapVHfbSIXZ5aF1Sq1V7Fpa0VsXYRpYqCKIIFWVyAIu5gAsi+r1lIAmSFEJLM9/njPoM3w0wyYWZyn0ner+uaa+6z/c73nPsE5nP/fufcqSokSZIkSVJvrdfrAiRJkiRJkgFdkiRJkqRWMKBLkiRJktQCBnRJkiRJklrAgC5JkiRJUgsY0CVJkiRJagEDuiRJXZLck+Tlva6jX5ITk3yt13WsjiSfSnJGr+tokySVZO9e1yFJajcDuiSpFXoRjJN8I8lJa3KfK5Pk8CTTu+dV1Weq6j1jtL9WfRgxmCRnJJmVZGGS25K8Z8DylyW5JcmjSS5PsnuvapUkaaQM6JIkqc0+C0ypqs2B1wInJTkYIMm2wPeBTwBbA1OB767pApNMGMO2J45V25Kk9jGgS5JaL8mrk1ybZH6SXyV5Tteye5J8KMl1SRYk+W6SSV3LT2h6YGcmeU//UOMkxwFvA05IsjjJBV27PHCw9pJsm+TCpo6Hk/w8yaD/L02yb5JLmvVuTXJM17KjktyUZFGSGU39mwA/BHZq6lmcZKfu4eJJpjT1vyvJ/UnmJTk+yaFNvfOTfLlrP3sl+d8kDyV5MMmZSbZsln0L2A24oNnXCc38w5pzPD/J75Ic3tXeO5Pc1dR9d5K3reRtm9Scu0VJrk7y3KaNDyc5d8C5+lKSUwZrpKpurKql/ZPNz17N9OuBG6vqe1X1GPAp4LlJ9m3a3TrJ15v3fl6S85r5hyeZ3lwbc5rr43XN+3Jb856dONSBNSMv/jPJxUkeAV6a5KfdvfvNufrFENtvmOTzSe5LMjvJfyXZaEBtH0nyAPD1lZxjSdJaxoAuSWq1JAcB/w28F9gG+ArwgyQbdq12DHAksAfwHOCdzbZHAh8AXg7sDRzev0FVnQqcCfxzVW1aVa9ZVXvAB4HpwGRge+BEOoFxYM2bAJcA3wa2A94C/EeS/ZtVTgPeW1WbAc8C/reqHgFeCcxs6tm0qmYOcVpeAOwDvBn4IvB3zTEeAByT5I/6S6HTA70TsB+wK50QS1UdC9wHvKbZ1z8n2Rm4CDiJTo/0h4Bzk0xujulLwCubuv8AuHaI+gCOBr7XtPNt4Lwk6wNnAEd2fVAwsTk/3xyqoST/keRR4BZgFnBxs+gA4Hf96zXn8M5mPsC3gI2b6e2Ak7ua3QGYBOwM/D3wVeDtwMHAi4FPJNljJcf3Z8Cngc2AQYP4SnwOeAZwIJ3rsr+G7tq2BnYHjlvNtiVJ45gBXZLUdscBX6mqK6tqRVWdDiwFDuta50tVNbOqHgYuoBN8oBO0v970wj5KE06HYaj2lgE7ArtX1bKq+nlVPSWgA68G7qmqr1fV8qq6BjgXeFNXO/sn2byq5lXV1cOsq98/VtVjVfUT4BHgrKqaU1UzgJ8DBwFU1R1VdUlVLa2qucAXgD8aulneDlxcVRdXVV9VXUJn2PhRzfI+4FlJNqqqWVV140ramlZV51TVsma/k4DDqmoW8LOuc3Ek8GBVTRuqoar6CzpB+MV0hrT396hvCiwYsPoCYLMkO9L5wOP45hwvq6orutZbBny6qe87wLbAKVW1qDmum4DnruT4zq+qXzbn6bGVrPckSULnmv7bqnq4qhYBn6HzIUW/PuCTzfu2ZLhtS5LGPwO6JKntdgc+2Ay5np9kPp2e4J261nmg6/WjdIIbzTr3dy3rfr0yQ7X3L8AdwE+aod4fXUnNLxhQ89vo9IwCvIFO6L03yRVJXjjMuvrN7nq9ZJDpTQGSbJ/kO80w+oV0eq+3XUm7uwNvGlD3HwI7Nr3TbwaOB2Yluah/KPkQnjjXVdVHZ+RB/3t2Op0PA2h+f2vlhwvNhzO/AHYB3tfMXgxsPmDVzYFFdK6Rh6tq3hBNPlRVK5rX/SF40PM4hOFeSwNNptOrP63rHP+omd9v7uqEfknS2sOALklqu/vp9HRu2fWzcVWdNYxtZ9EJdP12HbB8sN7vITW9qx+sqj3pPLDsA0leNkTNVwyoedOqel/Tzm+r6mg6w67PA85+OvUMw2eaNp/dPGTt7XSGvT9xSIPU/a0BdW9SVZ9r6v5xVb2CziiCW+gMCx/KE+c6nfv0dwH6h+yfBzwnybPojDY4czWOaSK/vwf9Rrp6uZth+Hs18+8Htu4fSj8GBp67R+gE7347MLgH6YT/A7rO8RZV1f1hwGhfB5KkccKALklqk/WTTOr6mUgnBB6f5AXp2CTJq5JsNoz2zgbelWS/JBvTedp3t9nAnsMtLp2H1e3dDFNeAKygMxx5oAuBZyQ5Nsn6zc+hTR0bJHlbki2a4dULu9qYDWyTZIvh1rQKm9HpZV7Q3F/+4QHLBx7/GcBrkvxJkgnNe3B4kl2a3vijmxC8tGl3sGPvd3CS1zfv4fubbX4D0PQOn0Pn3vSrquq+wRpIsl2StyTZtKnnT4C3Apc1q/wPnSH3b0jnQX5/D1xXVbc0Q+l/SOfe/62a9+AlwzprT8+1wOuTbJzO952/e7CVmtEEXwVOTrJdc5w7N8cmSVrHGdAlSW1yMZ3exf6fT1XVVODPgS8D8+gMMX/ncBqrqh/SebDZ5c12v2kW9d/DfBqde8Hnp3nC9yrsA1xKJ5z+GviPqrp8kP0uAv6Yzn3FM+kMmf8noP/BdscC9zTDzo+nM/ydqroFOAu4q6lpp4Ftr6Z/AJ5H58OEi+jcv93ts8DHm319qKrup/NwtxOBuXR6oT9M5++F9eg8cG8m8DCde9nfx9DOpzMkfl5zvK9vPpDodzrwbFY+vL2afUxv2vk88P6q+gFAc1/9G+g8rG0enYfndd/LfSyde81vAebQ+aBgrJwMPE7nQ4/TWfmogI/QXI/NNXAp8MwxrE2SNE5k8GfbSJK09kmyH3ADsGFVLe91PeuyJLvRCc47VNXCXtcjSVIb2IMuSVqrJfnT5nunt6LTi32B4by3mnvSPwB8x3AuSdLvGdAlSWu799IZ3nwnnXvGVzYsW2OsuYd9IfAK4JM9LkeSpFZxiLskSZIkSS1gD7okSZIkSS0wsdcFjMS2225bU6ZM6XUZkiRJkiQNadq0aQ9W1eRVrTeuA/qUKVOYOnVqr8uQJEmSJGlISe4dznoOcZckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSVJPXD99AZfdPLvXZUhSaxjQJUmS1BMn/s/1vPv0qSxeurzXpUhSKxjQJUmS1BNFdX5X9bgSSWoHA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAj0J6En+NsmNSW5IclaSSUn2SHJlkjuSfDfJBr2oTZIkSZKkXljjAT3JzsD/Aw6pqmcBE4C3AP8EnFxVewPzgHev6dokSZIkSeqVXg1xnwhslGQisDEwCzgCOKdZfjrwuh7VJkmSJEnSGrfGA3pVzQA+D9xHJ5gvAKYB86tqebPadGDnwbZPclySqUmmzp07d02ULEmSJEnSmOvFEPetgKOBPYCdgE2AI4e7fVWdWlWHVNUhkydPHqMqJUmSJElas3oxxP3lwN1VNbeqlgHfB14EbNkMeQfYBZjRg9okSZIkSeqJXgT0+4DDkmycJMDLgJuAy4E3Nuu8Azi/B7VJkiRJktQTvbgH/Uo6D4O7Gri+qeFU4CPAB5LcAWwDnLama5MkSZIkqVcmrnqV0VdVnwQ+OWD2XcDze1COJEmSJEk916uvWZMkSZIkSV0M6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZIkSS1gQJckSZIkqQUM6JIkSZIktYABXZIkSZKkFjCgS5IkSZLUAgZ0SZIkSZJawIAuSZIkSVILGNAlSZIkSWoBA7okSZJ66kPf+12vS5CkVjCgS5IkqWf+7a0HMfWeeb0uQ5JawYAuSZKkntl00sRelyBJrWFAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiRJktQCPQnoSbZMck6SW5LcnOSFSbZOckmS25vfW/WiNkmSJEmSeqFXPeinAD+qqn2B5wI3Ax8FLquqfYDLmmlJkiRJktYJazygJ9kCeAlwGkBVPV5V84GjgdOb1U4HXrema5MkSdKacfOshcxdtLTXZUhSq/SiB30PYC7w9STXJPlakk2A7atqVrPOA8D2g22c5LgkU5NMnTt37hoqWZIkSaPp3y+/g+fvsQ3P2H6zXpciSa3Ri4A+EXge8J9VdRDwCAOGs1dVATXYxlV1alUdUlWHTJ48ecyLlSRJ0th4xf7bs+FEn1ksSf168V/E6cD0qrqymT6HTmCfnWRHgOb3nB7UJkmSJElST6zxgF5VDwD3J3lmM+tlwE3AD4B3NPPeAZy/pmuTJEmSJKlXJvZov38NnJlkA+Au4F10Piw4O8m7gXuBY3pUmyRJkiRJa9yIA3qSZ1fV9auzTVVdCxwyyKKXjbQeSZIkSZLGo9EY4v4fSa5K8hfNV6hJkiRJkqTVNOKAXlUvBt4G7ApMS/LtJK8YcWWSJEmSJK1DRuUhcVV1O/Bx4CPAHwFfSnJLktePRvuSJEmSJK3tRhzQkzwnycnAzcARwGuqar/m9ckjbV+SJEmSpHXBaDzF/d+ArwEnVtWS/plVNTPJx0ehfUmSJEmS1nqjEdBfBSypqhUASdYDJlXVo1X1rVFoX5IkSZKktd5o3IN+KbBR1/TGzTxJkiRJkjRMoxHQJ1XV4v6J5vXGo9CuJEmSJEnrjNEI6I8keV7/RJKDgSUrWV+SJEmSJA0wGvegvx/4XpKZQIAdgDePQruSJEmSJK0zRhzQq+q3SfYFntnMurWqlo20XUmSJEmS1iWj0YMOcCgwpWnveUmoqm+OUtuSJEmSJK31RhzQk3wL2Au4FljRzC7AgC5JkiRJ0jCNRg/6IcD+VVWj0JYkSZIkSeuk0XiK+w10HgwnSZIkSZKeptHoQd8WuCnJVcDS/plV9dpRaFuSJEmSpHXCaAT0T41CG5IkSZIkrdNG42vWrkiyO7BPVV2aZGNgwshLkyRJkiRp3THie9CT/DlwDvCVZtbOwHkjbVeSJEmSpHXJaDwk7i+BFwELAarqdmC7UWhXkiRJkqR1xmgE9KVV9Xj/RJKJdL4HXZIkSZIkDdNoBPQrkpwIbJTkFcD3gAtGoV1JkiRJktYZoxHQPwrMBa4H3gtcDHx8FNqVJEmSJGmdMRpPce8Dvtr8SJIkSZKkp2HEAT3J3Qxyz3lV7TnStiVJkiRJWleMOKADh3S9ngS8Cdh6FNqVJEmSJGmdMeJ70Kvqoa6fGVX1ReBVo1CbJEmSJEnrjNEY4v68rsn16PSoj0bPvCRJkiRJ64zRCNL/2vV6OXAPcMwotCtJkiRJ0jpjNJ7i/tLRKESSJEmSpHXZaAxx/8DKllfVF0a6D0mSJEmS1naj9RT3Q4EfNNOvAa4Cbh+FtiVJkiRJWieMRkDfBXheVS0CSPIp4KKqevsotC1JkiRJ0jphxF+zBmwPPN41/XgzT5IkSZIkDdNo9KB/E7gqyf80068DTh+FdiVJkiRJWmeMxlPcP53kh8CLm1nvqqprRtquJEmSJEnrktEY4g6wMbCwqk4BpifZY5TalSRJkiRpnTDigJ7kk8BHgI81s9YHzhhpu5IkSZIkrUtGowf9T4HXAo8AVNVMYLNVbZRkQpJrklzYTO+R5MokdyT5bpINRqE2SZIkSZLGhdEI6I9XVQEFkGSTYW73N8DNXdP/BJxcVXsD84B3j0JtkiRJkiSNC6MR0M9O8hVgyyR/DlwKfHVlGyTZBXgV8LVmOsARwDnNKqfTeRq8JEmSJEnrhNF4ivvnk7wCWAg8E/j7qrpkFZt9ETiB3w+F3waYX1XLm+npwM6DbZjkOOA4gN12222E1UuSJEmS1A4jCuhJJgCXVtVLgVWF8v5tXg3MqappSQ5f3X1W1anAqQCHHHJIre72kiRJkiS10YgCelWtSNKXZIuqWjDMzV4EvDbJUcAkYHPgFDpD5Cc2vei7ADNGUpskSZIkSePJiIe4A4uB65NcQvMkd4Cq+n+DrVxVH6P5SramB/1DVfW2JN8D3gh8B3gHcP4o1CZJkiRJ0rgwGgH9+83PSH0E+E6Sk4BrgNNGoU1JkiRJksaFpx3Qk+xWVfdV1elPt42q+inw0+b1XcDzn25bkiRJkiSNZyP5mrXz+l8kOXcUapEkSZIkaZ01koCertd7jrQQSZIkSZLWZSMJ6DXEa0mSJEmStJpG8pC45yZZSKcnfaPmNc10VdXmI65OkiRJkqR1xNMO6FU1YTQLkSRJkiRpXTaSIe6SJEmSJGmUGNAlSZIkSWqBkdyDLkmSJK2W5Sv6+PkdD/LYshW9LkWSWscedEmSJK0xF10/i3d9/bdcevOcXpciSa1jQJckSdIas2yF384rSUMxoEuSJKlnNpy4HstW9PGx71/f61IkqecM6JIkSeqZzSatz3+9/WBum72o16VIUs8Z0CVJktRTG0z0T1JJAgO6JEmSJEmtYECXJEmSJKkFDOiSJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktMLHXBUiSJGndssVG67PTlhvxzO0363UpktQqBnRJkiStUS/bbzu+cMyBvS5DklrHIe6SJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktYECXJEmSJKkFDOiSJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqgTUe0JPsmuTyJDcluTHJ3zTzt05ySZLbm99brenaJEmSJEnqlV70oC8HPlhV+wOHAX+ZZH/go8BlVbUPcFkzLUmSJEnSOmGNB/SqmlVVVzevFwE3AzsDRwOnN6udDrxuTdcmSZIkSVKv9PQe9CRTgIOAK4Htq2pWs+gBYPshtjkuydQkU+fOnbtG6pQkSZIkaaz1LKAn2RQ4F3h/VS3sXlZVBdRg21XVqVV1SFUdMnny5DVQqSRJkiRJY68nAT3J+nTC+ZlV9f1m9uwkOzbLdwTm9KI2SZIkSZJ6oRdPcQ9wGnBzVX2ha9EPgHc0r98BnL+ma5MkSZIkqVcm9mCfLwKOBa5Pcm0z70Tgc8DZSd4N3Asc04PaJEmSJEnqiTUe0KvqF0CGWPyyNVmLJEmSJElt0dOnuEuSJEmSpA4DuiRJkiRJLWBAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiSppxK496FH+c1dD/W6FEnqKQO6JEmSeurAXbfijw/YnrOn3t/rUiSppwzokiRJ6qkJ64Xn7bZVr8uQpJ4zoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktYECXJEmSJKkFDOiSJEnjyF1zF/Po48t7XcbT9vAjS3tdgiS1lgFdkiRpnLh51kKO+Ncr+OzFt/S6lKflppkLOfmS2zl0yta9LkWSWsmALkmSNE7095w/snR89qA/+vhy9t9pc976/N16XYoktZIBXZIkSZKkFjCgS5IkSZLUAgZ0SZIkjbmq4orb5va6DElqNQO6JEmSxtxDjzzOV39+F399xN69LkWSWsuALkmSpFGzYMkyFj62bNBlm2wwkcOfud2Q2z66dAV9fTVWpUlS6xnQJUmSNGpe+cWfcfSXf7na2+293ab86s4HOfPKe8egKkkaHyb2ugBJkiStPWYueIwNJqx+H9CBu27J/3nhFOY/OnjvuyStC+xBlyRJkiSpBQzokiRJkiS1gAFdkiRJrXLNffP40Q2zel2GJK1xBnRJkiS1yifOv4Hjz7iaRUM8DV6S1lYGdEmSJLXG4yv6el2CJPWMT3GXJElSK+y13SaccM51LFvhd6FLWjfZgy5JkqRW+NODduGPD9ih12VIUs8Y0CVJknrorrmL+b/f+C3T7p037G1+eeeDq7X+eLL+eul1CZLUMwZ0SZKkHrrmvvn87y1z+OUdDw5r/efusgWH7bnNsNcfb/7uVfvzneMOY9MNvRNT0rrHgC5JktRytz6wiJ/fPpcbZixkwnpht6037nVJK9VXxQ+vn8XM+UuemHfDjAXD2nbyZhty2J7brLz9vqe23yvT7n2Ya+6bx2PLVnDeNTNY6JPnJY2AAV2SJKnlTjjndxx72lV88gc38mcv2L3X5azS8r7ifWdezRcuue2JeX991jW860VTRqX9a+6fz/vOvJqTu9rvlbeeeiX/9xu/5cLrZvH+717Lt359b69LkjSOGdAlSdJar6qoGr9PBu9rSp+4XnjjwbsAUPXk41rVMY7mORhuW319v6+tr694xx9Medr7e1K7zXT/eRmr93dl7fYve3xFH0uWrfh9TX3j9zqT1HsGdEmStNb7izOv5tBPX8qSx1f0upRRsctWG3HKZbdx6Kcv5dBPX8atDyxi30/86Ek91gN97ke3sN/f/4j7H350xPs//oxpHPrpy3hs2VPP53qBDSasx/abb8j/XDuD7189nePPmMYmG05k/Qmr/6fnL+94kD0+djHfv3o6AIuXLue935o27HpG4qyr7mfPEy9m6j0PP2XZZ394C/v//Y+fMv/kS29j2r1PXV+ShsOnb0iSpLXejPlLeHDx4zy+vI+NNpjQ63JG7M2H7sYv7niIC343E4CZC5awdHkfM+YNfU/2jHlLeGxZH/MefZxdR3gPe+d8LuXxFX1MWv/J53OTDSbyy48dweaT1ufzP76VGfOWMGP+Ek57x6FPWXc4+u8z7z+2JY+vYL3Av7zxOfzmroefVM+yQeoZiRnzH6UKZi147KnL5i1hyYAPBN548C4seXzFoOtL0nC0LqAnORI4BZgAfK2qPtfjkiStBVb0FZ+5+Gb2mrwpf/aC3VZ7+8tvmcP3r5nBNptswN+9ar+n1QskrWv+7bLbSeCvjthnxO2st174y5fu/cS8s666j9tnL+bEo/Zl4oB/j8tX9HHSRTez346bse8Om/Pfv7x7jT1MbMGSZZx04U0c9ewdeem+263WtpfdPJsZ85bw8VfvR1/BSRfexE5bbsRdDz7CfU2v9yZdTzbv/jay//zpnQBceffDnHzJbdz14CMA7LTlJB5e/DhHPXvHJ9b9xq/u4ROv2p/PXHzzU+rsr3/nrTbingcf4YQj92WnLTd6Yvm198/nG7+8m1nzfx9A73/4UT7/k1sJcMKR+z7pmNYLXHjdLGYtWEKexrenLXxsOedMm/7E9LevvI8r736IFX3FeglX3v0QJ19yG/c91Dk/J5xzHS/ddzuOOWRXfnLjA1x28xw+8Zr9h/1E+DOvvJe75j7CiUftx4SVfN3bNffN4+r7fv81d48t6+PcadPZfZvff/CxbEUfn77oZjafNJF7HnqU416yJ8/aeYtB2zvjN/dyx5zFrOgrnrPLFrzpkF2HVe9Y+PGND3D5LXP4+KsHP2+PL+/j0xfdxLN2Hps6T77kNiatP4H3Hb7XqLR370OP8K8/uY33/tGeHLDT4OdfapNWBfQkE4B/B14BTAd+m+QHVXVTbyuTNN49+vhyTvvF3Txr582fVkC/4LqZbD5pIudMm85fHbE322664RhUKa1dvnjZ7UxIRhzQT770NiZOWO8pAf266Qv4m5ftwxYbPzmgL3xsOd/41T0cuOuWvGSfbTn/2pkj2v/quGPOYr43bTor+mq1A/rvpi/gd9MX8LbDdmPZij6+1wTTo569A59/03PZdtMNWK8r5Z5w5L68+jk7sc2mG3D/w4/yD689gFMuvZ1TLrudd7xwd56xw2b83f/cAMCK5v7o/3fE3nzlZ3fxmufu1KmznlznHXMWPbFfgCP2257XdgX0y26ezXkDzueVdz/MzPlLWN5XTB3w3ezv/sM92Wu7Tdl4g4kcsNPmq3U+AG6bvYg5i5Y+cd/9t6+6lxtmLGTDievx6ufuyE9ueoBTLrud9x2+F6/Yf3uuuHUu3/3t/RxzyK6c/7uZXHTdLI594e5DBuOBvn3lfdw4cyHvf/k+bDZp/SHXu/Tm2Rw6ZWtufWARt85eBMDcRUs55S0H8Y8Xdv5snf/oMr7xq3ue2OYZ2286ZB3fvvI+bpq1EIBDdt+qpwH9/GtncPH1D/D2wwY/b/MffZzTf30vh04ZmzpPuex2Np80cdQC+m/ueogf/G4mz9xhMwO6xoVWBXTg+cAdVXUXQJLvAEcD4z6g//rOh3hk6fJelyGts/qHIS5cspxLb5q92tvPmv8Ybzh4F350wwP87La5bL6SP9wkdazoK1ZQT+vfXLe+6vTadbezcEnnq6x+etscNtngyX/OLG7+f7twybInepL7Dbb+aLpz7mKgM+R8uMd948yFT5q+6u6HWd71oLF9d9icV+y//VO223nLjdi5Cc/P220rAPbZflOJaLq4AAAgAElEQVR+dCMcPGVrXrTXNk8E9FnzH2N5Xx9/fMAObDBhPa65b/4T87vrvKOpv98NMxawcdeQ8aecz1vncuPMBey69cY8vryPG2YsYFlf3xPLt9h4fY4+cOdhnYfBXHPffLbeZAN22HwSd85dzMIlnfd24nphw4kT2Ge7zfjxjbPZf8fNed5uW9HXV1x43UwuvWk2s5th5lfe/TAPDHPIef9XpP301rlstP4E7nnw0SfOw0Zd5+GeBx9lvx03e+L9Bjp1bjEJgOunL2Dpsj663TFn8ZDXxKKlv/9qtnmPPj7ifzMjMXvhUmDo8za/+bc379FlY1bnkmUrRq3tm5p/X3eu5Pxr/Dp4963YapMNel3GqEqbnmia5I3AkVX1nmb6WOAFVfVXXescBxzXTD4TuHWNF/r0bAs82OsipGHyetV44bWq8cJrVeOJ16vGi/F0re5eVZNXtVLbetBXqapOBU7tdR2rK8nUqjqk13VIw+H1qvHCa1XjhdeqxhOvV40Xa+O12ranHM0Aum9m2aWZJ0mSJEnSWq1tAf23wD5J9kiyAfAW4Ac9rkmSJEmSpDHXqiHuVbU8yV8BP6bzNWv/XVU39ris0TLuhuVrneb1qvHCa1XjhdeqxhOvV40Xa9212qqHxEmSJEmStK5q2xB3SZIkSZLWSQZ0SZIkSZJawIC+BiQ5MsmtSe5I8tFe1yMNJcl/J5mT5IZe1yKtTJJdk1ye5KYkNyb5m17XJA0myaQkVyX5XXOt/kOva5JWJsmEJNckubDXtUgrk+SeJNcnuTbJ1F7XM1q8B32MJZkA3Aa8AphO50n1b62qm3pamDSIJC8BFgPfrKpn9boeaShJdgR2rKqrk2wGTANe539b1TZJAmxSVYuTrA/8AvibqvpNj0uTBpXkA8AhwOZV9epe1yMNJck9wCFV9WCvaxlN9qCPvecDd1TVXVX1OPAd4Oge1yQNqqp+Bjzc6zqkVamqWVV1dfN6EXAzsHNvq5KeqjoWN5PrNz/2jqiVkuwCvAr4Wq9rkdZVBvSxtzNwf9f0dPwjUpJGTZIpwEHAlb2tRBpcM2T4WmAOcElVea2qrb4InAD09boQaRgK+EmSaUmO63Uxo8WALkkat5JsCpwLvL+qFva6HmkwVbWiqg4EdgGen8RbiNQ6SV4NzKmqab2uRRqmP6yq5wGvBP6yuVVz3DOgj70ZwK5d07s08yRJI9Dcz3sucGZVfb/X9UirUlXzgcuBI3tdizSIFwGvbe7r/Q5wRJIzeluSNLSqmtH8ngP8D51bi8c9A/rY+y2wT5I9kmwAvAX4QY9rkqRxrXnw1mnAzVX1hV7XIw0lyeQkWzavN6Lz0NhbeluV9FRV9bGq2qWqptD5e/V/q+rtPS5LGlSSTZqHxJJkE+CPgbXiW4gM6GOsqpYDfwX8mM5DjM6uqht7W5U0uCRnAb8GnplkepJ397omaQgvAo6l08NzbfNzVK+LkgaxI3B5kuvofGh/SVX59VWSNDLbA79I8jvgKuCiqvpRj2saFX7NmiRJkiRJLWAPuiRJkiRJLWBAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiRJktQCBnRJktaQJNt0fS3cA0lmdE1vMMj6E5PMH6KtM5K8bjX2fVKS94+k/iHand7/Pd8D5h+R5LCVbPfGJCeu5r4uS7LF06lTkqTxYGKvC5AkaV1RVQ8BBwIk+RSwuKo+39Oixs4RwIPAb4ZY/mHgyNVs89vA8cA/jaAuSZJayx50SZJaIMkFSaYluTHJewYs+1Iz/5Ik2wyy7aFJrmi2/2GS7Vexr32S/LhZ/2dJntHMPyPJKUl+leSuJH/azJ+Q5L+S3JLkJ0l+NKD3/v1JrklyXZJnJNkLeA/w4WZ0wB8M2P/+wKKqmte1339PcmWSO5O8JMnpzf5O69r0fODPhn9WJUkaXwzokiS1wzuq6mDgUOADSbZq5m8B/LKqDgB+DXyie6MkGwKnAG9otj8D+MdV7OtU4C+a9T8GfLlr2XbAi4DXAZ9t5r0J2BnYH3gn8MIB7c2uqoOArwEfqKo7m9f/UlUHVtWvBqz/ImDagHlbVNULgBOAC+j0ku8PHJzkWQBV9SCw2WBD6iVJWhs4xF2SpHb42ySvbV7vAuwFXAssB77XzD+DzjDvbvsBBwCXJgGYAEwfaidNuD0MOLdZH57898B5VVXAdUl2bub9IXB2VfUBM5NcMaDZ7ze/pwFHreI4AXYE5g6Yd0Hz+3pgZlXd1NR7EzAFuKFZPrfZftB78yVJGs8M6JIk9ViSlwMvAQ6rqiVJfgFMGmL1Grg5cF1VvXi4uwMerKoDh1i+dMC6w9G/zQqG97fFEp56fP1t9A2ooW9Am5Oa7SVJWus4xF2SpN7bAni4CecH0Bnm3m8i8Prm9Z8Bvxiw7U3AzkmeD5Bkg6aNQTX3fc/qur98vSTPXUV9vwTemI4d6XyYsCqLgM2GWHYzsPcw2niSJBOAbYH7VndbSZLGAwO6JEm9dxGwcTOc+yTgyq5lC4AXJ7mRzlDzk7o3rKqlwBuBLyS5DrgGeMEq9vcW4PgkvwNuBF69ivXPBubQCdbfaPaxYBXbnA8c0zw87g8GLPspcMgqth/MocAvmqH2kiStddK5zUySJGloSTatqsVJJtP5AOEFVTXwPvLVae/fge9V1U9Xc5uzq2rgPfCSJK0VvAddkiQNxw+TbA6sD3xyJOG8cRJw8Gpuc43hXJK0NrMHXZIkSZKkFvAedEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktYECXJEmSJKkFDOiSJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUkapiT3JHl5r+vol+TEJF/rdR2rI8mnkpzR6zrGQpJKsnev65AkjV8GdElS6/UiGCf5RpKT1uQ+VybJ4Ummd8+rqs9U1XvGaH+t+jBioCQbJjktyb1JFiW5Nskre12XJEkjYUCXJEnj0UTgfuCPgC2AjwNnJ5kyFjtLMmEs2m3anjhWbUuSxhcDuiRpXEvy6qb3dH6SXyV5Tteye5J8KMl1SRYk+W6SSV3LT0gyK8nMJO/pH6Kc5DjgbcAJSRYnuaBrlwcO1l6SbZNc2NTxcJKfJxn0/7NJ9k1ySbPerUmO6Vp2VJKbml7hGU39mwA/BHZq6lmcZKfu4eJJpjT1vyvJ/UnmJTk+yaFNvfOTfLlrP3sl+d8kDyV5MMmZSbZsln0L2A24oNnXCc38w5pzPD/J75Ic3tXeO5Pc1dR9d5K3reRtm9Scu0VJrk7y3KaNDyc5d8C5+lKSUwY2UFWPVNWnquqequqrqguBu4GDm+0OTzK9eY/nNO/z65rze1tz7k8cqsBmBMV/Jrk4ySPAS5P8NMl7utZ5Z5JfDLH9hkk+n+S+JLOT/FeSjQbU9pEkDwBfX8m5kiStQwzokqRxK8lBwH8D7wW2Ab4C/CDJhl2rHQMcCewBPAd4Z7PtkcAHgJcDewOH929QVacCZwL/XFWbVtVrVtUe8EFgOjAZ2B44EahBat4EuAT4NrAd8BbgP5Ls36xyGvDeqtoMeBbwv1X1CPBKYGZTz6ZVNXOI0/ICYB/gzcAXgb9rjvEA4Jgkf9RfCvBZYCdgP2BX4FPN8R8L3Ae8ptnXPyfZGbgIOAnYGvgQcG6Syc0xfQl4ZVP3HwDXDlEfwNHA95p2vg2cl2R94AzgyK4PCiY25+ebK2mLZt3tgWcAN3bN3gGYBOwM/D3wVeDtdEL8i4FPJNljJc3+GfBpYDNg0CC+Ep9r6jmQzvXVX0N3bVsDuwPHrWbbkqS1lAFdkjSeHQd8paqurKoVVXU6sBQ4rGudL1XVzKp6GLiATmCCTtD+elXdWFWP0oTTYRiqvWXAjsDuVbWsqn5eVU8J6MCrgXuq6utVtbyqrgHOBd7U1c7+STavqnlVdfUw6+r3j1X1WFX9BHgEOKuq5lTVDODnwEEAVXVHVV1SVUurai7wBTrDxYfyduDiqrq46bG+BJgKHNUs7wOelWSjqppVVTcO2RJMq6pzqmpZs99JwGFVNQv4Wde5OBJ4sKqmreyAm3B/JnB6Vd3StWgZ8OlmP98BtgVOqapFTX03Ac9dSdPnV9Uvm+N9bGU1DKgndK7Nv62qh6tqEfAZOh829OsDPtmc/yXDbVuStHYzoEuSxrPdgQ82Q67nJ5lPpyd4p651Huh6/SiwafN6Jzr3MPfrfr0yQ7X3L8AdwE+aod4fXUnNLxhQ89vo9KgCvIFO6L03yRVJXjjMuvrN7nq9ZJDpTaHT45zkO80w+oV0eq+3XUm7uwNvGlD3HwI7Nj38bwaOB2YluSjJvitp64lzXVV9dEYe9L9np9P5MIDm97dWdrDNbQTfAh4H/mrA4oeqakXXscMQ52NVda6mycDGwLSuc/WjZn6/uasT+iVJ6wYDuiRpPLufTg/pll0/G1fVWcPYdhawS9f0rgOWD9b7PaSmV/aDVbUn8FrgA0leNkTNVwyoedOqel/Tzm+r6mg6w9/PA85+OvUMw2eaNp9dVZvTCcPpPqRB6v7WgLo3qarPNXX/uKpeQWcUwS10hpMP5Ylz3QTsXYD+IfvnAc9J8iw6ow3OHKqRpqf6NDq3FLyh6SkfTQPPwSN0gne/HRjcg3TC/wFd52qLqur+MGC0309J0lrAgC5JGi/WTzKp62cinRB4fJIXpGOTJK9Kstkw2jsbeFeS/ZJsDHxiwPLZwJ7DLS6dh9Xt3YTGBcAKOsOYB7oQeEaSY5Os3/wc2tSxQZK3JdmiCZsLu9qYDWyTZIvh1rQKmwGLgQXN/eUfHrB84PGfAbwmyZ8kmdC8B4cn2aXpjT+6uRd9adPuYMfe7+Akr2/ew/c32/wGoOlVPofOvelXVdV9K2nnP+ncP/+aNTRM/Frg9Uk2Tuf7zt892ErNqICvAicn2Q4gyc5J/mQN1ChJGscM6JKk8eJiOr2S/T+fqqqpwJ8DXwbm0Rli/s7hNFZVP6TzYLPLm+1+0yxa2vw+jc694POTnDeMJvcBLqUTTn8N/EdVXT7IfhcBf0znfuSZdIbM/xPQ/2C7Y4F7mmHnx9MZ/k5zb/VZwF1NTTsNbHs1/QPwPDofJlwEfH/A8s8CH2/29aGqup/Ow91OBObS6VH/MJ2/Jdaj88C9mcDDdO5lf99K9n0+nSHx85rjff2A3u/TgWezkuHtSXan83DAA4EH8vun26/s6fEjdTKdofSzmxqH7N0HPkJzXTXv5aXAM8ewNknSWiCDP79GkqR1S5L9gBuADatqea/rWZcl2Y3OMPkdqmphr+uRJGlNsQddkrTOSvKnzfdVb0WnF/sCw3lvNfekfwD4juFckrSuMaBLktZl7wXmAHfSuWd8ZcOyNcaae9gXAq8APtnjciRJWuMc4i5JkiRJUgvYgy5JkiRJUgtM7HUBI7HtttvWlClTel2GJEmSJElDmjZt2oNVNXlV643rgD5lyhSmTp3a6zIkSZIkSRpSknuHs55D3CVJkiRJagEDuiRJkiRJLWBAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiRJktQCBnRJkiRJklpgzAJ6kl2TXJ7kpiQ3JvmbZv7WSS5Jcnvze6tmfpJ8KckdSa5L8ryxqk2SJEmSpLYZyx705cAHq2p/4DDgL5PsD3wUuKyq9gEua6YBXgns0/wcB/znGNYmSZKkHvv8j2/l3d/4LctW9PW6FElqhTEL6FU1q6qubl4vAm4GdgaOBk5vVjsdeF3z+mjgm9XxG2DLJDuOVX2SJEnqrZ/eNofLbpnDY8tW9LoUSWqFNXIPepIpwEHAlcD2VTWrWfQAsH3zemfg/q7NpjfzBrZ1XJKpSabOnTt3zGqWJEmSJGlNGvOAnmRT4Fzg/VW1sHtZVRVQq9NeVZ1aVYdU1SGTJ08exUolSZIkSeqdMQ3oSdanE87PrKrvN7Nn9w9db37PaebPAHbt2nyXZp4kSZIkSWu9sXyKe4DTgJur6gtdi34AvKN5/Q7g/K75/6d5mvthwIKuofCSJEmSJK3VJo5h2y8CjgWuT3JtM+9E4HPA2UneDdwLHNMsuxg4CrgDeBR41xjWJkmSJElSq4xZQK+qXwAZYvHLBlm/gL8cq3okSZIkSWqzNfIUd0mSJEmStHIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktYECXJEmSJKkFDOiSJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AIGdEmSJEmSWsCALkmSJElSCxjQJUmSJElqAQO6JEmSJEktYECXJEmSJKkFDOiSJEmSJLWAAV2SJEmSpBYwoEuSJEmS1AJjFtCT/HeSOUlu6Jr33STXNj/3JLm2mT8lyZKuZf81VnVJkiRJktRGE8ew7W8AXwa+2T+jqt7c/zrJvwILuta/s6oOHMN6JEmSJElqrTEL6FX1syRTBluWJMAxwBFjtX9JkiRJksaTXt2D/mJgdlXd3jVvjyTXJLkiyYuH2jDJcUmmJpk6d+7csa9UkiRJkqQ1oFcB/a3AWV3Ts4Ddquog4APAt5NsPtiGVXVqVR1SVYdMnjx5DZQqSZIkSdLYW+MBPclE4PXAd/vnVdXSqnqoeT0NuBN4xpquTZIkSZKkXulFD/rLgVuqanr/jCSTk0xoXu8J7APc1YPaJEmSJEnqibH8mrWzgF8Dz0wyPcm7m0Vv4cnD2wFeAlzXfO3aOcDxVfXwWNUmSZIkSVLbjOVT3N86xPx3DjLvXODcsapFkiRJkqS269VD4iRJkiRJUhcDuiRJkiRJLWBAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiRJktQCBnRJkiRJklrAgC5JkiRJUgsY0CVJkiRJagEDuiRJkiRJLWBAlyRJkiSpBQzokiRJkiS1gAFdkiRJkqQWMKBLkiRJktQCBnRJkiRJklrAgC5JkiRJUgsY0CVJkiRJagEDuv5/e3cerVdd33v8/TEDYUoAiRgCAaRAC9QGDOCEWlsqVSvaWhtqrVotcqu9er2rrWBvoV5dtrcKxdbaolK1KC1eQKVqZajCxQFIGJIQjIIEyUASIGQQiBm+94+zYw/hhPME8py9T877tdazzt6/PX2StVfW+WQPjyRJkiSpA/pW0JNclGRlkgWDxs5NsjTJbc3nlYOWnZXkriSLkryiX7kkSZIkSeqifl5B/wxw6hDj51fVzObzNYAkRwOzgWOabf4hybg+ZpMkSZIkqVP6VtCr6nrgoR5XPw3416raUFX3AHcBJ/YrmyRJkiRJXdPGM+jvSjKvuQV+32ZsOnDfoHWWNGNPkOSMJHOSzFm1alW/s0qSJEmSNCJGuqB/AjgcmAksBz66ozuoqguralZVzZo6derOzidJkiRJUitGtKBX1Yqq2lxVW4BP8l+3sS8FDh606kHNmCRJkiRJY8KIFvQk0wbNvg7Y+ob3rwCzk+yW5DDgCOCmkcwmSZIkSVKbxvdrx0kuAV4G7J9kCXAO8LIkM4ECFgPvAKiqO5JcCiwENgHvrKrN/comSZIkSVLX9K2gV9XpQwx/+knW/xDwoX7lkSRJkiSpy9p4i7skSZIkSdqGBV2SJEmSpA6woEuSJEmS1AEWdEmSJEmSOsCCLkmSJElSB1jQJUmSJEnqAAu6JEmSJEkdYEGXJEmSJKkDLOiSJEmSJHWABV2SJEmSpA6woEuSJEmS1AEWdEmSJEmSOsCCLkmSJElSB1jQJUmSJEnqAAu6JEmSJEkdYEGXJEmSJKkDLOiSJEmSJHVA3wp6kouSrEyyYNDY3yT5fpJ5Sa5Isk8zfmiSR5Pc1nz+sV+5JEmSJEnqop4KepJffAr7/gxw6jZjVwPHVtVzgR8AZw1adndVzWw+Zz6F40mSJEmSNGr1egX9H5LclOSPkkzpZYOquh54aJuxq6pqUzP7PeCg3qNKkiRJkrTr6qmgV9XJwBuBg4G5Sb6Q5JSneew/AL4+aP6wJLcmuS7JyU9z35IkSZIkjSrje12xqn6Y5M+BOcDHgOOSBDi7qi7fkYMmeT+wCfh8M7QcmFFVDyZ5HvClJMdU1dohtj0DOANgxowZO3JYSZIkSZI6q9dn0J+b5HzgTuDlwG9U1S800+fvyAGTvAV4NfDGqiqAqtpQVQ8203OBu4Ejh9q+qi6sqllVNWvq1Kk7cmhJkiRJkjqr1yvofwd8ioGr5Y9uHayqZc1V9Z4kORX4U+ClVfXIoPGpwENVtTnJc4AjgB/1ul9JkiRJkka7Xgv6q4BHq2ozQJJnAJOq6pGq+pehNkhyCfAyYP8kS4BzGHhr+27A1QN3x/O95o3tLwE+kGQjsAU4s6oeGmq/kiRJkiTtinot6NcAvwqsb+b3AK4CXri9Darq9CGGP72ddS8DLusxiyRJkiRJu5xev2ZtUlVtLec003v0J5IkSZIkSWNPrwX9J0mO3zrTvGn90SdZX5IkSZIk7YBeb3F/D/DFJMuAAM8GfqdvqSRJkiRJGmN6KuhVdXOSnweOaoYWVdXG/sWSJEmSJGls6fUKOsAJwKHNNscnoao+15dUkiRJkiSNMT0V9CT/AhwO3AZsboYLsKBLkiRJkrQT9HoFfRZwdFVVP8NIkiRJkjRW9foW9wUMvBhOkiRJkiT1Qa9X0PcHFia5CdiwdbCqXtOXVJIkSZIkjTG9FvRz+xlCkiRJkqSxrtevWbsuySHAEVV1TZI9gHH9jSZJkiRJ0tjR0zPoSf4Q+L/APzVD04Ev9SuUJEmSJEljTa8viXsn8CJgLUBV/RB4Vr9CSZIkSZI01vRa0DdU1U+3ziQZz8D3oEuSJEmSpJ2g14J+XZKzgd2TnAJ8Ebiyf7EkSZIkSRpbei3o7wNWAfOBdwBfA/68X6EkSZIkSRpren2L+xbgk81HkiRJkiTtZD0V9CT3MMQz51X1nJ2eSJIkSZKkMaingg7MGjQ9CfhtYL/hNkpyEfBqYGVVHduM7Qf8G3AosBh4Q1WtThLgAuCVwCPAW6rqlh7zSZIkSZI0qvX0DHpVPTjos7Sq/hZ4VQ+bfgY4dZux9wHXVtURwLXNPMCvA0c0nzOAT/SSTZIkSZKkXUGvt7gfP2j2GQxcUR9226q6Psmh2wyfBrysmf4s8C3gz5rxz1VVAd9Lsk+SaVW1vJeMkiRJkiSNZr3e4v7RQdObaG5Nf4rHPGBQ6b4fOKCZng7cN2i9Jc3Y4wp6kjMYuMLOjBkznmIESZIkSZK6pde3uP9yPw5eVZXkCS+fG2abC4ELAWbNmrVD20qSJEmS1FW93uL+3idbXlXn7cAxV2y9dT3JNGBlM74UOHjQegc1Y5IkSZIk7fJ6ekkcA8+c/zcGbjmfDpwJHA/s3Xx2xFeANzfTbwa+PGj89zPg+cAanz+XJEmSJI0VvT6DfhBwfFWtA0hyLvDVqvq9J9soySUMvBBu/yRLgHOAvwIuTfI24F7+61n2rzHwFWt3MfA1a2/doT+JJEmSJEmjWK8F/QDgp4Pmf8p/vdxtu6rq9O0s+pUh1i3gnT3mkSRJkiRpl9JrQf8ccFOSK5r51zLwFWmSJEmSJGkn6PUt7h9K8nXg5GborVV1a/9iSZIkSZI0tvT6kjiAPYC1VXUBsCTJYX3KJEmSpF3cV+ct594HHmk7hiR1Sk8FPck5wJ8BZzVDE4CL+xVKkiRJu7avL1jOH7z4MPacOK7tKJLUGb1eQX8d8BrgJwBVtYwd/3o1SZIk6WcOf9ZeJGk7hiR1Rq8F/afNW9YLIMme/YskSZIkSdLY02tBvzTJPwH7JPlD4Brgk/2LJUmSJEnS2NLrW9w/kuQUYC1wFPAXVXV1X5NJkiRJkjSGDFvQk4wDrqmqXwYs5ZIkSZIk9cGwt7hX1WZgS5IpI5BHkiRJkqQxqadb3IH1wPwkV9O8yR2gqv57X1JJkiRJkjTG9FrQL28+kiRJkiSpD560oCeZUVU/rqrPjlQgSZIkSZLGouGeQf/S1okkl/U5iyRJkiRJY9ZwBT2Dpp/TzyCSJEmSJI1lwxX02s60JEmSJEnaiYZ7SdwvJVnLwJX03Ztpmvmqqsl9TSdJkiRJ0hjxpAW9qsaNVBBJkiRJksayXr9mbadJchTwb4OGngP8BbAP8IfAqmb87Kr62gjHkyRJkiSpFSNe0KtqETATIMk4YClwBfBW4Pyq+shIZ5IkSZIkqW3DvSSu334FuLuq7m05hyRJklryxTlL2o4gSZ3QdkGfDVwyaP5dSeYluSjJvkNtkOSMJHOSzFm1atVQq0iSJGmU+JvXP5ePf/OutmNIUie0VtCTTAReA3yxGfoEcDgDt78vBz461HZVdWFVzaqqWVOnTh2RrJIkSeqPEw7br+0IktQZbV5B/3XglqpaAVBVK6pqc1VtAT4JnNhiNkmSJEmSRlSbBf10Bt3enmTaoGWvAxaMeCJJkiRJkloy4m9xB0iyJ3AK8I5Bw/8nyUyggMXbLJMkSZIkaZfWSkGvqp8Az9xm7E1tZJEkSZIkqQvafou7JEmSJEnCgi5JkiRJUidY0CVJkiRJ6gALuiRJkiRJHWBBlyRJkiSpAyzokiRJkiR1gAVdkiRJkqQOsKBLkiRJktQBFnRJkiRJkjrAgi5JkiRJUgdY0CVJkiRJ6gALuiRJkiRJHWBBlyRJkiSpAyzokiRJkiR1gAVdkiRJkqQOsKBLkiRJktQBFnRJkiRJkjpgfFsHTrIYWAdsBjZV1awk+wH/BhwKLAbeUFWr28ooSZIkSdJIafsK+i9X1cyqmtXMvw+4tqqOAK5t5iVJkiRJ2uW1XdC3dRrw2Wb6s8BrW8wiSZIkSdKIabOgF3BVkrlJzmjGDqiq5c30/cAB7USTJEmSJGlktfYMOvDiqlqa5FnA1Um+P3hhVVWS2najpsyfATBjxoyRSSpJkiRJUp+1dgW9qhH+W4sAAA6dSURBVJY2P1cCVwAnAiuSTANofq4cYrsLq2pWVc2aOnXqSEaWJEmSJKlvWinoSfZMsvfWaeDXgAXAV4A3N6u9GfhyG/kkSZIkSRppbd3ifgBwRZKtGb5QVf+R5Gbg0iRvA+4F3tBSPkmSJEmSRlQrBb2qfgT80hDjDwK/MvKJJEmSJElqV9e+Zk2SJEmSpDHJgi5JkiRJUgdY0CVJkiRJ6gALuiRJkiRJHWBBlyRJkiSpAyzokiRJkiR1gAVdkiRJkqQOsKBLkiRJktQBFnRJkiRJkjrAgi5JkiRJUgdY0CVJkiRJ6gALuiRJkiRJHWBBlyRJkiSpAyzokiRJkiR1gAVdkiRJkqQOsKBLkiRJktQBFnRJkiRJkjrAgi5JkiRJUgeMeEFPcnCSbyZZmOSOJO9uxs9NsjTJbc3nlSOdTZIkSZKktoxv4ZibgP9ZVbck2RuYm+TqZtn5VfWRFjJJkiRJktSqES/oVbUcWN5Mr0tyJzB9pHNIkiRJktQlrT6DnuRQ4DjgxmboXUnmJbkoyb7b2eaMJHOSzFm1atUIJZUkSZIkqb9aK+hJ9gIuA95TVWuBTwCHAzMZuML+0aG2q6oLq2pWVc2aOnXqiOWVJEmSJKmfWinoSSYwUM4/X1WXA1TViqraXFVbgE8CJ7aRTZIkSZKkNrTxFvcAnwburKrzBo1PG7Ta64AFI51NkiRJkqS2tPEW9xcBbwLmJ7mtGTsbOD3JTKCAxcA7WsgmSZIkSVIr2niL+w1Ahlj0tZHOIkmSJElSV7T6FndJkiRJkjTAgi5JkqTWTBz/DB7buJkPfXVh21EkqXUWdEmSJI2YhcvW8uq/+3/ccNcDAEyeNIG/f+Px3PLjh1tOJknta+MlcZIkSRqjFi5fy4Klax83tvdu/koqSeAVdEmSJEmSOsGCLkmSpFYM9bU+kjSWeT+RJEmSRtTLjprKq35xGr92zAFtR5GkTrGgS5IkaUTtt+dEfnvWwW3HkKTO8RZ3SZIkSZI6wIIuSZIkSVIHWNAlSZIkSeoAC7okSZIkSR1gQZckSdKIefiRn7YdQZI6y4IuSZKkEXHn8rV85KpFHD9j37ajSFInWdAlSZI0In6yYRPHHDiF33v+IW1HkaROsqBLkiRJktQBFnRJknYh9z30CMd94CouuOaHbUeRHucnGzbxx5fcyrMnT2o7iiR1lgVdkqRdyAPrN7D6kY3c88D6tqNIP3P2FfM55bzrWP/YJj52+nE7tO0tP17NCz58LVfdcX+f0klSd3SuoCc5NcmiJHcleV/beSRpq498YxGnnHcd7/zCLW1HkTRGXLNwBa84/3ruXL6WxzZuZvaF3+WU867j8zfe+7T2u2LtY/zG393AF2788U5K+kQLlq7hFedfz7cWrWTekodZtuYxpuwxgXHPyBPW3WPieBYuW8sp513HtxatBOCxjZv5nX/6Ln908S0sX/MY7//SAhbdvw6A+9f0P/9w/vnb93Da39/AA+s3PO19ffDfF/LGT32PjZu37IRkkkaz8W0HGCzJOODjwCnAEuDmJF+pqoXtJpM01l1x6xL++dv3cO5rjuEDVy583Phlc5fyt7Nnsv9eu7WYUNJo9YMV6/jzKxawpYrdJ47jgtnHsd+eEwG4fcnDLFqxjh+sWMez9t6NhcvW8o6XHs7cxat540k79qK1S+fcx5W3L+OC2cexZPUjzF+6hiPufYjfPWnGDu/n0pvvA+DoAyfzgdOOfdzysy6fz9pHN3L3qvUsWrGOsy+fz8OPbuQLbz+JI5+995D7PPrAyXz93Sfz2e8u5uzL53P4s/binN84mhvveQiA02YeyOpHNvKDFev4zHfuYd6SNdyxbC1H3rt6u/kvvfk+rpy3jI/NPo59m7/P4VxwzQ9ZtGItF8w+jgnjnvGk+5mzeDW3L1nDW/75JiaNHwfAIc/ck/tWP8Jf/eYv8pype/V0TIDv3P0gC5evZcOmLT877mhw3Q9W8fFv3sVf/9ZzmTZlEn98ya0875B9OfOlh7cdTRq1UlVtZ/iZJC8Azq2qVzTzZwFU1YeHWn/WrFk1Z86cEUz41J3z5QUsWf1o2zEkPUV3LFvLm194KG96wSEc/4GrOfmI/X82fv/ax3jeIfuyz+4TWk4pwZpHNzLn3tUcMHk3jj1wSttx1INlax7jwCmTOPNlh/P+K+YzZfcJTJ408O/J3avWs/jBRzjmwMnst+dEFt2/jve/6hf48Ne+zzEHTt6h48xfuoaV6zZwwqH7sqVg7r2refbkSU9pP7//gkM45sApvOPiuZz8c/s/bvm13x+4An7MgZP5zFtPZPGDP2H3CeM4dvrw5+NjGzczf+ka3n/FfHafOJ7b73sYgNknHMxjGzfzw5XruWPZ2p+t/2T5B/95t/59Duemex5i3YZNvOTIqUxorvRvbz/zlq5h1bqBq+cXv+0kbrrnQT72n3cBcOz0yRywd+/P2t94z0Os37CJlx45lfFD3GHQVXetWs+9Dz7CsdMns8/uE7nhrgeYsvsEZh3i1+hpZPzJqUfx88/esX/D2pJkblXNGna9jhX01wOnVtXbm/k3ASdV1bsGrXMGcEYzexSwaMSDPjX7Aw+0HULqkeerRgvPVY0WnqsaTTxfNVqMpnP1kKqaOtxKnbrFvRdVdSFwYds5dlSSOb38j4nUBZ6vGi08VzVaeK5qNPF81WixK56rXXvIZSlw8KD5g5oxSZIkSZJ2aV0r6DcDRyQ5LMlEYDbwlZYzSZIkSZLUd526xb2qNiV5F/ANYBxwUVXd0XKsnWXU3ZavMc3zVaOF56pGC89VjSaerxotdrlztVMviZMkSZIkaazq2i3ukiRJkiSNSRZ0SZIkSZI6wII+ApKcmmRRkruSvK/tPNL2JLkoycokC9rOIj2ZJAcn+WaShUnuSPLutjNJQ0kyKclNSW5vztW/bDuT9GSSjEtya5J/bzuL9GSSLE4yP8ltSea0nWdn8Rn0PksyDvgBcAqwhIE31Z9eVQtbDSYNIclLgPXA56rq2LbzSNuTZBowrapuSbI3MBd4rf+2qmuSBNizqtYnmQDcALy7qr7XcjRpSEneC8wCJlfVq9vOI21PksXArKp6oO0sO5NX0PvvROCuqvpRVf0U+FfgtJYzSUOqquuBh9rOIQ2nqpZX1S3N9DrgTmB6u6mkJ6oB65vZCc3HqyPqpCQHAa8CPtV2FmmssqD333TgvkHzS/CXSEnaaZIcChwH3NhuEmlozS3DtwErgaurynNVXfW3wJ8CW9oOIvWggKuSzE1yRtthdhYLuiRp1EqyF3AZ8J6qWtt2HmkoVbW5qmYCBwEnJvERInVOklcDK6tqbttZpB69uKqOB34deGfzqOaoZ0Hvv6XAwYPmD2rGJElPQ/M872XA56vq8rbzSMOpqoeBbwKntp1FGsKLgNc0z/X+K/DyJBe3G0navqpa2vxcCVzBwKPFo54Fvf9uBo5IcliSicBs4CstZ5KkUa158dangTur6ry280jbk2Rqkn2a6d0ZeGns99tNJT1RVZ1VVQdV1aEM/L76n1X1ey3HkoaUZM/mJbEk2RP4NWCX+BYiC3qfVdUm4F3ANxh4idGlVXVHu6mkoSW5BPgucFSSJUne1nYmaTteBLyJgSs8tzWfV7YdShrCNOCbSeYx8J/2V1eVX18lSU/PAcANSW4HbgK+WlX/0XKmncKvWZMkSZIkqQO8gi5JkiRJUgdY0CVJkiRJ6gALuiRJkiRJHWBBlyRJkiSpAyzokiRJkiR1gAVdkqQRkuSZg74W7v4kSwfNTxxi/fFJHt7Ovi5O8todOPYHk7zn6eTfzn6XbP2e723GX57k+U+y3euTnL2Dx7o2yZSnklOSpNFgfNsBJEkaK6rqQWAmQJJzgfVV9ZFWQ/XPy4EHgO9tZ/mfAKfu4D6/AJwJ/PXTyCVJUmd5BV2SpA5IcmWSuUnuSPL2bZZ9rBm/Oskzh9j2hCTXNdt/PckBwxzriCTfaNa/PsmRzfjFSS5I8p0kP0ryumZ8XJJ/TPL9JFcl+Y9trt6/J8mtSeYlOTLJ4cDbgT9p7g544TbHPxpYV1WrBx3340luTHJ3kpck+WxzvE8P2vTLwO/2/rcqSdLoYkGXJKkb3lxVzwNOAN6bZN9mfArw7ao6Bvgu8L8Gb5RkN+AC4Lea7S8G/vcwx7oQ+KNm/bOAvx+07FnAi4DXAh9uxn4bmA4cDbwFeME2+1tRVccBnwLeW1V3N9N/U1Uzq+o726z/ImDuNmNTquok4E+BKxm4Sn408LwkxwJU1QPA3kPdUi9J0q7AW9wlSeqG/5HkNc30QcDhwG3AJuCLzfjFDNzmPdgvAMcA1yQBGAcs2d5BmnL7fOCyZn14/O8DX6qqAuYlmd6MvRi4tKq2AMuSXLfNbi9vfs4FXjnMnxNgGrBqm7Erm5/zgWVVtbDJuxA4FFjQLF/VbD/ks/mSJI1mFnRJklqW5FeBlwDPr6pHk9wATNrO6rXt5sC8qjq518MBD1TVzO0s37DNur3Yus1mevvd4lGe+Ofbuo8t22TYss0+JzXbS5K0y/EWd0mS2jcFeKgp58cwcJv7VuOB32ymfxe4YZttFwLTk5wIkGRis48hNc99Lx/0fPkzkvzSMPm+Dbw+A6Yx8J8Jw1kH7L2dZXcCP9fDPh4nyThgf+DHO7qtJEmjgQVdkqT2fRXYo7md+4PAjYOWrQFOTnIHA7eaf3DwhlW1AXg9cF6SecCtwEnDHG82cGaS24E7gFcPs/6lwEoGivVnmmOsGWabLwNvaF4e98Jtln0LmDXM9kM5AbihudVekqRdTgYeM5MkSdq+JHtV1fokUxn4D4STqmrb58h3ZH8fB75YVd/awW0uraptn4GXJGmX4DPokiSpF19PMhmYAJzzdMp544PA83Zwm1st55KkXZlX0CVJkiRJ6gCfQZckSZIkqQMs6JIkSZIkdYAFXZIkSZKkDrCgS5IkSZLUARZ0SZIkSZI64P8DnBm8sdUFihgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "minL = 0.0\n",
    "maxL = 5.0\n",
    "\n",
    "## Define two histograms with all the lengths recorded:\n",
    "fig_raw, ax = plt.subplots(nrows=2, figsize=(14,8))\n",
    "ax_L30cm, ax_L2m = ax\n",
    "\n",
    "hist_L30cm = ax_L30cm.hist(L30cm, bins=1000, range=(minL, maxL), histtype='step', label='Binned Data')\n",
    "ax_L30cm.set(xlabel='Table lenght (m)', ylabel='Frequency', title='Lengths estimates by 30cm ruler')\n",
    "\n",
    "hist_L2m = ax_L2m.hist(L2m, bins=1000, range=(minL, maxL), histtype='step', label='Binned Data')\n",
    "ax_L2m.set(xlabel='Table lenght (m)', ylabel='Frequency', title='Lengths estimates by 2m ruler')\n",
    "\n",
    "fig_raw.tight_layout()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 30 cm ruler:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We consider the measurements for the 30cm ruler and focus on that for now. Below is a mean and RMS calculation along with a general Gaussian fit to all the data. Somehow, it doesn't seem optimal/right..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Initial estimate (30 cm ruler):  3.379 +- 0.010 m\n"
     ]
    }
   ],
   "source": [
    "print(f\"  Initial estimate (30 cm ruler):  {L30cm.mean():.3f} +- {L30cm.std(ddof=1)/np.sqrt(len(L30cm)):.3f} m\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We set up the fit. Can you guess why we use the binwidth `L30cm_binwidth` in the inital value for `N`in the fit?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "\n",
      "L30cm raw fit: \n",
      "\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td title=\"Minimum value of function\">FCN = 357.1519475223496</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 135</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 135</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 4.174794303695252e-06</td>\n",
       "        <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 1e-05</td>\n",
       "        <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
       "        UP = 1.0</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<table>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
       "        <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
       "        <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
       "        <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
       "        <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
       "        <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
       "        <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#MkNfrASRkt').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>N</td>\n",
       "        <td>1.40925</td>\n",
       "        <td>0.339632</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>mu</td>\n",
       "        <td>3.40276</td>\n",
       "        <td>0.0518009</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>sigma</td>\n",
       "        <td>0.270153</td>\n",
       "        <td>0.0857666</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"MkNfrASRkt\" style=\"display:none;\">\n",
       "<textarea rows=\"12\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & N & 1.40925 & 0.339632 &  &  &  &  & No\\\\\n",
       "\\hline\n",
       "1 & $\\mu$ & 3.40276 & 0.0518009 &  &  &  &  & No\\\\\n",
       "\\hline\n",
       "2 & $\\sigma$ & 0.270153 & 0.0857666 &  &  &  &  & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(\"\\n\\nL30cm raw fit: \\n\")\n",
    "\n",
    "L30cm_x, L30cm_y, L30cm_sy = get_bincenter_and_counts_in_range(hist_L30cm, 3.0, 3.7)\n",
    "L30cm_binwidth = L30cm_x[1] - L30cm_x[0]\n",
    "\n",
    "chi2_L30cm = Chi2Regression(gauss_extended, L30cm_x, L30cm_y, L30cm_sy) \n",
    "minuit_L30cm = Minuit(chi2_L30cm, pedantic=False, N=L30cm_y.sum()*L30cm_binwidth, mu=L30cm.mean(), sigma=L30cm.std(ddof=1)) \n",
    "minuit_L30cm.migrad();\n",
    "\n",
    "# the fitted values of the parameters\n",
    "L30cm_fit_N, L30cm_fit_mu, L30cm_fit_sigma = minuit_L30cm.args "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we plot the fit on the figure and extract relevant extra information:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA+gAAAI4CAYAAAD56sN/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3X98z/X+//Hb08Y2e29j5EdY8zOxSX6cImr6gUL5tVQolah0KlTUcaJDHZWT0tGhkyNCfsWU+OQc0ekHTgiz8mOYGTUz++EtG7bn94/3e+/vZsO0H++p+/VyeV/s/Xr+eD2er+1yTo/X8/l6voy1FhERERERERHxrkreDkBERERERERElKCLiIiIiIiIVAhK0EVEREREREQqACXoIiIiIiIiIhWAEnQRERERERGRCkAJuoiIiIiIiEgFoARdREQkH2NMgjHmNm/HkccY86Ix5n1vx3EpjDETjDHzvB1HRWKMscaYJt6OQ0REKjYl6CIiUiF4IzE2xnxgjJlUnue8EGNMlDEmKf8xa+2r1tqhZXS+CnUzoijGmHnGmJ+MMZnGmD3GmKHnlN9qjNlljPnFGLPOGHOVt2IVEREpKSXoIiIiUpH9FQi31gYDdwGTjDFtAYwxNYFlwJ+BUGAzsKi8AzTG+JRh375l1beIiFQ8StBFRKTCM8b0NMZsM8akG2O+Nca0yleWYIx51hizwxiTYYxZZIzxz1f+vHsG9ogxZmjeUmNjzDBgIPC8McZpjPk03ylbF9WfMaamMWalO47jxpivjDFF/n+pMaa5Mebf7nq7jTH35Cu70xjzgzHmhDHmsDv+QGA1cKU7Hqcx5sr8y8WNMeHu+B8yxhwyxqQZYx4zxrR3x5tujPl7vvM0NsZ8YYxJNcYcM8bMN8ZUc5d9CIQBn7rP9bz7+A3ua5xujNlujInK198QY8x+d9wHjDEDL/Br83dfuxPGmK3GmGvdfTxnjPn4nGs1zRjzdlGdWGvjrLXZeV/dn8bu732BOGvtEmttFjABuNYY09zdb6gxZrb7d59mjIlxH48yxiS5/zaOuv8+ert/L3vcv7MXzzcw98qLfxhjVhljTgJdjDHr88/uu6/V1+dp72eMmWKMSTTGJBtjZhhjAs6JbYwx5mdg9gWusYiI/MYoQRcRkQrNGHMd8C9gOFADmAl8Yozxy1ftHqA70BBoBQxxt+0OjAJuA5oAUXkNrLXvAfOB1621Dmttr4v1B4wGkoArgNrAi7gSxnNjDgT+DSwAagH3Au8aY1q4q8wChltrg4AI4Atr7UngDuCIOx6HtfbIeS7L9UBTYADwFvAn9xhbAvcYY27OCwXXDPSVwDVAA1xJLNbawUAi0Mt9rteNMfWAz4BJuGaknwU+NsZc4R7TNOAOd9wdgW3niQ/gbmCJu58FQIwxpjIwD+ie70aBr/v6zD1fR8aYd40xvwC7gJ+AVe6ilsD2vHrua7jPfRzgQ6Cq+3stYGq+busA/kA94CXgn8AgoC3QGfizMabhBcZ3P/AKEAQUmYhfwGSgGdAa199lXgz5YwsFrgKGXWLfIiJyGVOCLiIiFd0wYKa1dpO1NsdaOwfIBm7IV2eatfaItfY48CmuxAdcifZs9yzsL7iT02I4X39ngLrAVdbaM9bar6y1hRJ0oCeQYK2dba09a639HvgYiM7XTwtjTLC1Ns1au7WYceWZaK3NstauAU4CH1lrj1prDwNfAdcBWGvjrbX/ttZmW2tTgDeBm8/fLYOAVdbaVdbaXGvtv3EtG7/TXZ4LRBhjAqy1P1lr4y7Q1xZr7VJr7Rn3ef2BG6y1PwH/zXctugPHrLVbzteRtfYJXIlwZ1xL2vNm1B1AxjnVM4AgY0xdXDc8HnNf4zPW2i/z1TsDvOKObyFQE3jbWnvCPa4fgGsvML4V1tpv3Ncp6wL1CjDGGFx/0yOttcettSeAV3HdpMiTC4x3/95OFbdvERG5/ClBFxGRiu4qYLR7yXW6MSYd10zwlfnq/Jzv519wJW646xzKV5b/5ws5X39vAPHAGvdS77EXiPn6c2IeiGtmFKAfrqT3oDHmS2NMh2LGlSc538+nivjuADDG1DbGLHQvo8/ENXtd8wL9XgVEnxN3J6Cue3Z6APAY8JMx5rO8peTn4bnW1tpcXCsP8n5nc3DdDMD974cXHi64b858DdQHHncfdgLB51QNBk7g+hs5bq1NO0+XqdbaHPfPeUlwkdfxPIr7t3SuK3DN6m/Jd43/z308T8qlJP0iIvLboQRdREQqukO4Zjqr5ftUtdZ+VIy2P+FK6PI0OKe8qNnv83LPro621jbCtWHZKGPMreeJ+ctzYnZYax939/OdtfZuXMuuY4DFvyaeYnjV3Weke5O1QbiWvXuGVETcH54Td6C1drI77s+ttbfjWkWwC9ey8PPxXGvjek6/PpC3ZD8GaGWMicC12mD+JYzJl///DHoc+Wa53cvwG7uPHwJC85bSl4Fzr91JXIl3njoU7Riu5L9lvmscYq3NfzOgtP8ORETkMqEEXUREKpLKxhj/fB9fXEngY8aY641LoDGmhzEmqBj9LQYeMsZcY4ypimu37/ySgUbFDc64Nqtr4l6mnAHk4FqOfK6VQDNjzGBjTGX3p707jirGmIHGmBD38urMfH0kAzWMMSHFjekignDNMme4ny9/7pzyc8c/D+hljOlmjPFx/w6ijDH13bPxd7uT4Gx3v0WNPU9bY0xf9+/wGXebjQDu2eGluJ5N/5+1NrGoDowxtYwx9xpjHO54ugH3AWvdVZbjWnLfz7g28nsJ2GGt3eVeSr8a17P/1d2/g5uKddV+nW1AX2NMVeN63/kjRVVyryb4JzDVGFPLPc567rGJiMjvnBJ0ERGpSFbhml3M+0yw1m4GHgX+DqThWmI+pDidWWtX49rYbJ273UZ3Ud4zzLNwPQuebtw7fF9EU+A/uJLTDcC71tp1RZz3BNAV13PFR3AtmX8NyNvYbjCQ4F52/hiu5e9Ya3cBHwH73TFdeW7fl+hloA2umwmf4Xp+O7+/AuPc53rWWnsI1+ZuLwIpuGahn8P13wuVcG24dwQ4jutZ9sc5vxW4lsSnucfb131DIs8cIJILL2+37nMkufuZAjxjrf0EwP1cfT9cm7Wl4do8L/+z3INxPWu+CziK60ZBWZkKnMZ102MOF14VMAb336P7b+A/wNVlGJuIiFwmTNF724iIiPz2GGOuAXYCftbas96O5/fMGBOGK3GuY63N9HY8IiIiFYFm0EVE5DfNGNPH/d7p6rhmsT9Vcu5d7mfSRwELlZyLiIj8f0rQRUTkt244ruXN+3A9M36hZdlSxtzPsGcCtwPjvRyOiIhIhaIl7iIiIiIiIiIVgGbQRURERERERCoAX28HUBI1a9a04eHh3g5DRERERERE5Ly2bNlyzFp7xcXqXdYJenh4OJs3b/Z2GCIiIiIiIiLnZYw5WJx6WuIuIiIiIiIiUgEoQRcRERERERGpAJSgi4iIiIiIiFQAl/Uz6CIiIiIiIuJy5swZkpKSyMrK8nYov1v+/v7Ur1+fypUr/6r2StBFRERERER+A5KSkggKCiI8PBxjjLfD+d2x1pKamkpSUhINGzb8VX1oibuIiIiIiMhvQFZWFjVq1FBy7iXGGGrUqFGiFQxK0EVERERERH4jlJx7V0mvvxJ0ERERERERkQpACbqIiIiIeEVsUgZrf0z2dhgiUop8fHxo3bq155OQkMDmzZt56qmnAFi/fj3ffvutl6OsuLRJnIiIiIh4xYvLY4k9nMHOl7vh8NN/lor8FgQEBLBt27YCx8LDw2nXrh3gStAdDgcdO3b0RngVnmbQK4iWLVuyfv16b4chIiIiUm4s1vWvtV6ORETK0vr16+nZsycJCQnMmDGDqVOn0rp1a7766itvh1bh6FZlPuHh4SQnJ+Pj4+M5FhMTw2233Vas9sYY9u7dS5MmTS753HFxcZfcRn5/Jk6cyEsvvVTg7+zrr7/mj3/8I3v37uXqq69m1qxZtG7dGoAJEybwyiuv4OfnB0DNmjVJSEjwVvgiIiIiUl5Wj4WfY0u3zzqRcMfkC1Y5deqU579FGzZsyPLlyz1l4eHhPPbYYzgcDp599tnSje03QjPo5/j0009xOp2eT3GTc5GylpCQwNq1awscy87OJjo6mlGjRpGRkUGfPn2Ijo4uMBMxYMAAz9+zknMRERERKUt5S9y3bdtWIDmX4tEMejENGTIEHx8fdu/ezfbt24mOjuZf//oXAHfccYdneca1116LMYahQ4fy1ltvAa7EqmHDhixevJhnn32W48ePM27cOMaMGcOYMWOYPn06v/zyC2vWrClwQyAnJ4dJkyYxe/ZsTp06xYABA/jb3/5G5cqVAXj11Vd55513OHnyJM2aNWPFihXUq1evnK+MlJdnnnmGV155hU6dOnmO7d69m5MnTzJ48GAARo8ezZ///GdiY2Np1aqVt0IVEREREW+7yEy3VEyaQb8Ea9asYcGCBcTFxbFkyRL+97//AbB69WqcTicA27dvx+l0epLz/BYvXsyWLVtITk6ma9euALz22ms4nU7CwsIK1X/zzTdZsWIF33zzDfHx8fzwww9MmzYNcCVmkyZN4ttvvyU9PZ13333Xs4xZfntWrVqFn58fN954Y4HjRT2zZ4whPj7e8/3TTz+lRo0atG7dmk8//bTMYxUREREROZ+goCBOnDjh7TAqLCXo5+jduzfVqlXzfI4ePeop69WrF2FhYYSFhREZGcmePXsuqe+XXnqJmjVrUrVqVa677rqL1n///fcZP3489erVIygoiBEjRrBs2TLAlYTl5OSwa9cucnJy+MMf/kDNmjUvbbByWcjOzub555/ntddeK1TWvHlzAgMDmTt3LmfOnGHKlCn4+Pjwyy+/AK7l7fv27ePnn3/mpZde4t5772X37t3lPQQREREREcCVUy1fvlybxJ2HEvRzxMTEkJ6e7vnUqlXLUxYaGur5uUqVKmRlZV1S302bNr2k+ocOHWLw4MGemwUPPvggKSkpADRr1oyZM2cyadIkQkNDiY6OJjMz85L6l8vDlClT6NmzJ+Hh4YXK/Pz8WLJkCW+//TZ169YlOTmZRo0aERwcDMA111xDzZo1qVy5Mn379uXmm2/m888/L+cRiIiIiMjvRd7K4vyioqJYuXIl4MpjduzYwbZt2+jcuXN5h1fhKUEvRcaYC5b7+l7aI/8NGjTg888/99wsyMzMLDBrP2TIEL755hv27dvHnj17PM/Ey2/Ld999x2uvvYYxxvM31rRpU2JiYgDo1KkTW7Zs4dixY7z88sscPHjwvM+fV6pUSa+yERERERGpoJSgl6I6deqwc+fOUuvv4Ycf5qWXXuLIkSNYa9mzZw9r1qwBYN++fXzxxRecPn0aX19fcnNzPbOm8tsSExODtdbzAdi7dy+9e/cGYNeuXWRnZ5OWlsaTTz5Jly5dPLPty5cvJz09ndzcXD777DPWr19Pt27dvDUUERERERG5ACXo5+jVqxcOh8PzmTNnTrHbvvrqqzz55JPUq1ePF1544aL1c3JyPOdJTEz0nHvWrFmAa0fuzp0706lTJ0JCQujbt6/nmfjTp08zduxYatSoQbNmzejQoYNnJ2/5fVm+fDl16tQhLCyMs2fP8uGHH3rKPvroI6666iqCg4P505/+xMKFC2nevLkXoxURERERkfMxl/Ny13bt2tnNmzd7OwwRERER+RV6vvMVOw9nEjuhK0H+lb0djshl78cff+Saa67xdhi/e0X9HowxW6y17S7WVjPoIiIiIiIiIhWAEnQRERERERGRCuDSthUXERERERGRy8KNk7/gcPqpUuuvXrUAvhl7ywXr+Pj4EBkZibUWHx8f/v73v9OxY0eOHDnCU089xdKlS0stnqKEh4ezefNmatasWeh4UFAQ4NoLrG/fvowbNw5/f//z9pWens6CBQt44oknyjTm/JSgi4iIiIiI/AYdTj9FwuQepdZf+NjPLlonICCAbdu2AfD555/zwgsv8OWXX3LllVeWeXJ+MevWraNmzZo4nU6GDRvG8OHDL7gpeHp6Ou+++265Juha4i4iIiIiIiKlLjMzk+rVqwOQkJBAREQEAB988AF9+/ale/fuNG3alOeff97TxuFw8Kc//Ylrr72WG264geTkZABSUlLo168f7du3p3379nzzzTcApKam0rVrV1q2bMnQoUMpziboDoeDGTNmEBMTw/Hjx3E6ndx66620adOGyMhIVqxYAcDYsWPZt28frVu35rnnnjtvvdKkBN1twoQJGGNYvXo14Lpb4ufnR1RUlHcDk8vKoEGDqFOnDiEhIVx//fVs2LChWO3WrVtHZGQkwcHB1K1bl1GjRpGTkwPA/PnzC7z6r2rVqlSqVIljx44BsHnzZq6//nqCg4O5+uqrWb58uaffli1bFmjr6+vLH//4x9IfuIiIiIgIcOrUKVq3bk3z5s0ZOnQof/7zn4ust23bNhYtWkRsbCyLFi3i0KFDAJw8eZIbbriB7du3c9NNN/HPf/4TgKeffpqRI0fy3Xff8fHHHzN06FAAXn75ZTp16kRcXBx9+vQhMTGxWHEGBwfTsGFD9u7di7+/P8uXL2fr1q2sW7eO0aNHY61l8uTJNG7cmG3btvHGG2+ct15pUoKeT9OmTT3LLj755BOuuuoqL0ckl5vnnnuOAwcOkJGRwcSJE7n77rs9ifaFNG/enFWrVpGRkcGePXuIi4tj+vTpAAwcOBCn0+n5vP3229x0002e52oGDRpEr169SE9P5+9//zsDBw70JO9xcXGedpmZmVx55ZX069ev7C6AiIiIiPyu5S1x37VrF//3f//HAw88UGQSe+uttxISEoK/vz8tWrTg4MGDAFSpUoWePXsC0LZtWxISEgD4z3/+w5NPPknr1q256667yMzMxOl08t///pdBgwYB0KNHD8+MfXHkxWWt5cUXX6RVq1bcdtttHD582DNzf2794tQrCSXo+bRv357vv/+es2fP8vHHHxdIZHJycnj55ZcJDw+ndu3aPPXUU5w5cwaAAwcOcMstt1C9enWCg4OJjo4mPT0dcC3lMMYwdepUrrjiCho2bMimTZu8Mj4pe9deey0BAQFYazlz5gwpKSmeZPlC6tatS4MGDTDGcPbsWXJzc/nxxx+LrDt79mweeOABz/fExET69+9PpUqVuP322wkICODAgQOF2v3nP/+hUqVK3Hzzzb9+gCIiIiIixdShQweOHTtGSkpKoTI/Pz/Pzz4+Ppw9exaAypUrY4wpdDw3N5eNGzeybds2tm3bxuHDh3E4HL86thMnTpCQkECzZs2YP38+KSkpbNmyhW3btlG7dm2ysrIKtSluvZJQgp6PMYabb76ZTz75hOTkZBo3buwpe/PNN1mxYgXffPMN8fHx/PDDD0ybNg2A7OxsHn30UQ4dOkRSUhKpqam8/PLLBfrOzMzkp59+4u677y5UJr8tTzzxBP7+/vTs2ZPo6Ghq165drHaJiYmEhIQQGhrKpk2bGDJkSKE6u3fvZvv27URHR3uOPf300yxdupScnBzWrFlDUFCQ5/me/GbPns2gQYM8/4MnIiIiIlKWdu3aRU5ODjVq1ChxX127duWdd97xfM/biO6mm25iwYIFAKxevZq0tLSL9uV0OnniiSfo3bs31atXJyMjg1q1alG5cmXWrVvnmc0PCgrixIkTnnbnq1eatIv7Ofr378+9997LsGHDChx///33ef3116lXrx4AI0aMYMqUKYwePZrmzZvTvHlzT93o6GiWLFlSoP3jjz+Or68vd955JytXriz7gYjXvPvuu7z11lt8/PHHVK5cudjtwsLCyMjIYP/+/cydO5emTZsWqjN79mzuvvtuzysiwLWU58EHH2TChAn4+/uzbNkyAgICCrRLT08nJiaG7du3//qBiYiIiMhlpV61gGLtvH4p/V1M3jPo4FoSPmfOHHx8fEp87mnTpjFixAhatWrF2bNnuemmm5gxYwbjx4/nvvvuo2XLlnTs2JGwsLDz9tGlSxesteTm5tKnTx/P8/EDBw6kV69eREZG0q5dO09uV6NGDW688UYiIiK44447GDNmTJH1SpMS9HN07NiRa665hujoaL7++mvP8UOHDjF48GAqVXItOsjNzaVOnToAHD16lKeeeoqvvvqKkydPcvr0adq2bVug39DQUMD1TEVpL4OQiqdKlSrcd999tGjRgpYtW3LNNdcUu22jRo1o3rw5I0aM4KOPPvIcz8nJ4cMPP2TWrFmeY2lpadx5553Mnj2b3r17s2HDBu6++262bt1aYA+Fjz76iGuvvZZmzZqVzgBFREREpMK72DvLy8L59l8KDw9n586dAAwZMqTAatH8E5hOp9Pzc//+/enfvz8ANWvWZNGiRYX6rVGjBmvWrLloXHnPshelZs2a593cOW92Pk9xN4H+tbyyxN0YM9IYE2eM2WmM+cgY42+MaWiM2WSMiTfGLDLGVPFSbKxZs6bQ3ZAGDRrw+eefk56eTnp6OpmZmezZsweAF154AWMMP/74I+np6UyePLnUd/OTy1NOTg5xcXGX3M5aS2xsbIFjn3/+OdZabr/9ds+xvXv34uvrS79+/fDx8aFTp040atSIb7/9tkDbc59bFxERERGRiqfcE3RjTD3gKaCdtTYC8AHuBV4DplprmwBpwCPlHduFPPzww7z00kscOXIEay179uzx3Kk5ceIEwcHBBAYGcuDAAWbOnOnlaMUbfv75Z2bNmkVGRgZnz55l5syZJCQkFFpNERUVxZgxYwoc++ijj9ixYwe5ubkkJiYyZcoUOnToUKDO7Nmzuf/++wssEWrWrBmnT58mJiYGay2bN2/mxx9/LDBjHxcXx/bt27n33nvLYNQiIiIiIlJavLVJnC8QYIzxBaoCPwG3AEvd5XOA3l6KrUijR4+mc+fOdOrUiZCQEPr27cvRo0cBGD9+PJs3byYkJIR77rmHXr16eTla8QYfHx/mz59P48aNqV69Ov/4xz9YsmQJDRs2LFAvISGh0OsYUlNT6d+/P8HBwbRv355WrVoxZcoUT/nx48f55JNPCs2CV6tWjaVLlzJ+/HiCg4MZMGAAb775pue5H3Al9j169PA8ZiEiIiIiv11ayetdJb3+xhu/QGPM08ArwClgDfA0sNE9e44xpgGw2j3Dfm7bYcAwgLCwsLZlsXOeiIiIiJS9nu98xc7DmcRO6EqQf/E3VhWRoh04cICgoCBq1KihN/d4gbWW1NRUTpw4UWiSzhizxVrb7mJ9lPsmccaY6sDdQEMgHVgCdC9ue2vte8B7AO3atdPtIREREREREaB+/fokJSUV+d5xKR/+/v7Ur1//V7f3xi7utwEHrLUpAMaYZcCNQDVjjK+19ixQHzjshdhEREREREQuS5UrVy40cyuXF288g54I3GCMqWpc6y5uBX4A1gH93XUeBFZ4ITYRERERERERryj3BN1auwnXZnBbgVh3DO8BY4BRxph4oAYw67ydiIiIiIiIiPzGeGOJO9ba8cD4cw7vB/7ghXBK3YQJE4iPj2fevHneDkVEREREREQuE956zVqFFBcXx6233kq1atWoVasW/fr1K1C+fv36Ej3wL799SUlJREVFUbVqVdq0acPOnTuL1W7VqlW0b9+e4OBgGjRowKRJkwqUT5o0iZYtW1KpUiU++OADz/HExEQcDkeBjzGGjz/+GIDXX3+dpk2bEhQUREREBCtW6MkREREREZGKSgl6Pr1796ZHjx6kpqayZ88eunXr5u2Q5DIzbNgwIiIiSE1NZcCAAQwYMKBY7ZxOJ3/9619JSUlh48aNzJs3r8AKjEaNGjFt2jTatm1boF1YWBhOp9Pz+f7773E4HHTv7noxgq+vL8uWLSMjI4OZM2cyePBg9u/fX3oDFhERERGRUqME3e3YsWPEx8fzyCOP4OPjQ7Vq1Rg2bJin3OFwcMcdd3DkyBHPTGXeLGVubi7PPvssNWvWJDIyksTERG8NQ7woMzOTNWvWMHbsWAICAhg5ciQHDx5kx44dF217zz33cNttt+Hn50e9evXo3r07GzZs8JTff//93Hrrrfj5+V2wn9mzZ9OnTx8CAwMBGDVqFJGRkVSqVIkbb7yRRo0asWXLlpINVEREREREyoQSdLfQ0FDCwsIYNmwY69evJzs7u0C50+lk9erVXHnllZ7Zyrwl8IsXL2bFihXExcWxbNkyLSP+nYqPj8ff3x+Hw0Hnzp3Zv38/jRs3Zvfu3Zfc18aNG2nVqtUltcnNzWXu3Lk88MADRZanpaWxZ88eIiIiLjkeEREREREpe0rQ3SpVqsTatWvx9fWlT58+1K5dm5dffrlYbT/77DMGDRpE7dq1adq0Kb179y7jaKUiOnnyJA6HgxMnThAXF0daWhpBQUE4nc5L6mf69OlkZWUxZMiQS2q3Zs0ajDHccsstRZYPHz6cBx98kGuuueaS+hURERERkfLhlV3cK6omTZowf/58cnNzWbduHdHR0bRv354777zzgu1SUlLo2LGj53vt2rW1zP13KDAwEKfTSYMGDTh+/DgAJ06cwOFwFLuPVatW8cYbb/Df//73osvZzzV79mwGDhxIpUqF77u9+OKLpKam6s0CIiIiIiIVmGbQi1CpUiVuvfVWoqKiiIuLK3C8KLVq1eLo0aOe78nJyWUeo1Q8TZo04dSpUyQlJQFw+vRp9u3bx9VXX12s9t9++y3Dhg1j5cqVhIWFXdK509LSWLFiRZHL26dOncrnn39OTEwMVapUuaR+RURERESk/ChBd8vJyeGll17yJFfbt2/nq6++KrBrdp06dUhJSSmUgPfs2ZN58+aRnJzM3r17iYmJKdfYpWIIDg6mW7duTJ48maysLKZOnUpYWBiRkZEF6kVFRTFmzJgCx3bs2EF0dDSLFy8u8hnxM2fOkJWVhbXW83Nubq6nfMGCBURERNCiRYsC7ebMmcOMGTNYvXo1QUFBpThaEREREREpbUrQ3SpVqkR8fDzXX389DoeDvn37Mm7cuALP8zZr1oyhQ4dy7bXXUr9+fT755BMA+vfvz913303Lli3p168fd999t7eGIV4gvGSrAAAgAElEQVQ2c+ZMYmNjqV69OgsXLmTRokUYYwrUSUhIKHSTZ+rUqSQnJ9O1a1fPWwLuuOMOT/mjjz5KQECAZ5Y9ICCA//73v57y2bNnFzl7PmHCBBISEmjUqJGn31dffbWURy0iIiIiIqXBWGu9HcOv1q5dO7t582ZvhyEiIiIiv0LPd75i5+FMYid0Jci/srfDEREpM8aYLdbadherpxl0ERERERERkQpACbqIiIiIiIhIBaAEXURERERERKQCUIIuIiIiIiIiUgEoQS+mIUOGMG7cuAvWmT9/Pl27di2niEREREREROS3RAn6Od59910aN25MYGAgkZGRbNy4sdhtBw4cyJo1awoc27t3L927dyc0NJRatWrxwAMPkJGRUdphSylat24dkZGRBAcHU7duXUaNGkVOTo6n3BhDYGCg57Vl7733HgCJiYmeY3kfYwwff/zxRc+Zm5vLgAEDqF+/PsYYEhISCpSf75wAK1asoEOHDvj5+TFkyJAC7Q4fPsydd95JSEgIV111FQsWLPj1F0ZERERERMqUEvR8lixZwl/+8hfmzZtHZmYms2bNolKlkl2izMxM7rvvPg4cOEBCQgJZWVmMHj26lCKWstC8eXNWrVpFRkYGe/bsIS4ujunTpxeos337dpxOJ06nk2HDhgEQFhbmOeZ0Ovn+++9xOBx07969WOft2LEjS5cuPW95UecECAkJ4bnnnuORRx4p1ObJJ5+kQYMGHDt2jA8++ICHH36YAwcOFCseEREREREpX0rQ83n11VcZM2YMHTp0wMfHhz/84Q/84Q9/8JT/9NNPdOrUiaCgIB5++GHP8YyMDBwOB/7+/nTq1KlAn23btuXBBx8kJCSEqlWrcv/997Nhw4ZyG5Ncurp169KgQQOMMZw9e5bc3Fx+/PHHS+5n9uzZ9OnTh8DAwIvWrVSpEk8//TTt2l301YiFREVF0bdvX0JDQwuVffnllwwfPpzKlSvTpUsXWrduzaeffnrJ5xARERERkbKnBN3t7NmzxMbG0rFjx/PWWbNmDQsWLCAuLo4lS5bwv//9D3DNYDqdTmbMmHHR82zYsIFWrVqVWtxSNhITEwkJCSE0NJRNmzYVWjp+0003UbduXR566KEiH1nIzc1l7ty5PPDAA6UW08XOWRRrbaFj8fHxpRaTiIiIiIiUHiXobseOHSMnJ6fIWcg8vXr1IiwsjLCwMCIjI9mzZ88lnWPLli3Mnj2biRMnljRcKWNhYWFkZGSwb98+Ro0aRdOmTT1l3377LYmJiWzdupUjR47w1FNPFWq/Zs0ajDHccsstpRJPcc5ZlKioKGbOnEl2djZr167l+++/55dffimVmEREREREpHQpQXerWbMmPj4+HD9+/Lx18ifvVapUISsrq9j9JyQk0K9fP+bOnUuTJk1KFKuUn0aNGtG8eXNGjBjhOdahQwd8fX2pW7cuEydOLHLJ+OzZsxk4cGCJ9zC4lHMW5Z133iEpKYn69eszefJk7rrrLoKDg0slJhERERERKV2+3g6govD19SUiIoINGzZw/fXXl2rfR48epVu3bkyaNKnYG4ZJxWGtJTY2tsiySpUqFVpGnpaWxooVK9i6dWuZxFPUOc+nfv36fPbZZ57v7dq1o2fPnmUSl4iIiIiIlIxm0PMZO3Ysr732Gps2bSInJ4etW7fy3XfflajPjIwMunfvzuOPP86gQYNKKVIpSx999BE7duwgNzeXxMREpkyZQocOHQCIjY3l+++/Jycnh9TUVCZMmECvXr0KtF+wYAERERG0aNGiyP6joqIYM2ZMoePZ2dlkZ2d7fs5boXGxc+bk5JCVlUVOTo7n57NnzwJw6NAhjh8/TnZ2Nm+99RaHDh2iT58+Jb9IIiIiIiJS6pSg53Pvvffypz/9iXvvvZegoCAGDhzoSXQuZOLEiTgcDh577DE2bNiAw+GgcePGAMTExPD9998zbty4Au/HloorNTWV/v37ExwcTPv27WnVqhVTpkwBICUlhX79+hEcHEyLFi2oXbs206ZNK9B+9uzZF9wcLiEhgeTk5ELHr776as/fRvPmzQkICCjWOT/88EMCAgKYPHky8+bNIyAggEmTJgGwe/duIiIiqF69OosWLWL16tVa4i4iIiIiUkGZ4i6VrYjatWtnN2/e7O0wRERERORX6PnOV+w8nEnshK4E+Vf2djgiImXGGLPFWnvRdyprBl1ERERERESkAlCCLiIiIiIiIlIBKEEXERERERERqQCUoIuIiIiIiIhUAErQf6WoqCjef//9AseSk5Pp1KkTDoeDW265xUuRiYiIiIiIyOVICXopeu+996hVqxYnTpzgiy++8HY48jtx6tQphg8fTmhoKNWqVeOJJ54oVCctLY0rrriCQYMGeSFCEREREREpDiXopejgwYO0aNECY4y3Q5HfkZEjRxIfH09cXBypqalFJugvvvgiDRs29EJ0IiIiF/fsku3eDkFEpEJQgp5PeHg448ePp0mTJoSGhjJjxgxP2fHjx+nVqxfBwcH06dOH06dPe8rmz5+Pw+Fgzpw5vP7661riLuXm1KlTzJkzh3feeYe6devi4+NDREREgTpbtmzhwIED3HnnnV6KUkRE5Pzeue86NiekeTsMEZEKQQn6OTZt2sSOHTuYNWsWzz//PGfPngVcM5D+/v6kpKQwePBgNmzY4GkzcOBAnE4nAwcO5Pnnn8fpdGqJu5SLPXv2YIwhJiaG2rVr06JFC5YvX+4pt9by1FNPMWXKFC9GKSIicn4Of19vhyAiUmEoQT/Hww8/TNWqVenRowcnTpzg559/BuCzzz7jySefxM/Pj759+9KoUSMvRyoCmZmZnD59mv3793Pw4EGmT5/O4MGD+emnnwCYNWsWkZGRhWbVRURERESk4tEty3OEhoYCUKVKFQCysrIASElJoVatWp56tWvXLv/gRM5RtWpVcnJyGD16NP7+/nTp0oVmzZqxceNGbrnlFv76178WWO0hIiIiIiIVlxL0YqpVqxZHjx7lmmuuAVyvVBPxtkaNGmGMKbQxobWWAwcOsH///kI3k3bu3Mm2bdvKM0wRERERESkGLXEvpp49e/L3v/+d7Oxsli1bxv79+70dkgjVq1fn5ptv5s033+TMmTN89dVX7NmzhxtuuIHWrVtjrfV8xo8fz8CBA5Wci4iIiIhUUF5J0I0x1YwxS40xu4wxPxpjOhhjQo0x/zbG7HX/W90bsZ3PK6+8QlZWFldccQXz58+nQ4cO3g5JBIB//etf7N69m2rVqvHII4/w4YcfcuWVV3o7LBERERERuUTGWlv+JzVmDvCVtfZ9Y0wVoCrwInDcWjvZGDMWqG6tHXOhftq1a2c3b95cDhGLiIiISGnr+c5XjO56Nc8u3s6WP9/u7XBERMqMMWaLtbbdxeqV+wy6MSYEuAmYBWCtPW2tTQfuBua4q80Bepd3bCIiIiJSPn78KZOUE9neDkNEpELxxhL3hkAKMNsY870x5n1jTCBQ21r7k7vOz0CR26QbY4YZYzYbYzanpKSUU8giIiIiUpqmr4vnDw1r0Kx2kLdDERGpMLyRoPsCbYB/WGuvA04CY/NXsK5190WuvbfWvmetbWetbXfFFVeUebAiIiIiUjZub1EbP1/tWSwikscb/4uYBCRZaze5vy/FlbAnG2PqArj/PeqF2ERERERERES8otwTdGvtz8AhY8zV7kO3Aj8AnwAPuo89CKwo79jKQkJCAsYYzp496+1QREREREREpALz1pqiPwLzjTE7gNbAq8Bk4HZjzF7gNvf3chUeHk5AQAAOh4N69eoxatQocnJyyjsMuYwlJSURFRVF1apVadOmDTt37ix222nTplG7dm2qV6/OCy+8UKBs/fr1XH311QQGBtK7d28yMjKKdc4zZ87wyCOPEBQURFhYGIsXLy75IEVEREREpEyUOEE3xkReahtr7Tb3c+StrLW9rbVp1tpUa+2t1tqm1trbrLXHSxrbr/Hpp5/idDpZu3YtCxYs4J///Kc3wpDL1LBhw4iIiCA1NZUBAwYwYMCAYrXbtGkTEyZM4IsvvmDnzp0sXLjQk0z/8ssvREdHM378eI4ePYoxpkACf6FzTp06ldjYWA4dOsTcuXN5+OGHOXToUOkOWkRERERESkVpzKC/a4z5nzHmCfcr1H4TmjdvTufOnQvMRg4ZMoRnnnmGvn374nA4uOqqq8jMzMRay1/+8hfq169PnTp1ePrppzlz5kyB/iZPnkxISAitWrViy5Yt5T0cKQeZmZmsWbOGsWPHEhAQwMiRIzl48CA7duy4aNulS5fSt29fWrZsSb169Rg6dCgLFy4EYN26dYSEhHD//fcTGBjIs88+y6JFi4p1ziVLlvDUU09RrVo1oqKi6NChA8uXLy+7iyAiIiIiIr9aiRN0a21nYCDQANhijFlgjLm9xJF5WWxsLF9++SXXXXddgeNz587loYceIiMjg5iYGCpXrsyyZcuYPXs23377LT/88APffvst06dPL9AuNTWVlJQUHn30Ue6//35yc3PLczhSDuLj4/H398fhcNC5c2f2799P48aN2b1790Xb7t69m+bNm/P2228zevRoWrRo4WmXV/b111/TtWtXmjRpwvHjx0lJSbnoOfPaDho0iIULFxboV0REREREKpZSeQbdWrsXGAeMAW4Gphljdhlj+pZG/+Wpd+/eVKtWjd69ezNs2DAeeuihAuVdunShV69e+Pj4cN111xEQEMAnn3zC4MGDCQsLIzQ0lMcff5yYmJgC7Z555hmqVKnC448/TmJiIvHx8eU5LCkHJ0+exOFwcOLECeLi4khLSyMoKAin01nstvv27WPv3r0F2uWV/fzzz/zwww/4+fkB4HQ6L3rOvPLY2FgOHz5c7HhERERERKT8+Za0A2NMK+AhoAfwb6CXtXarMeZKYAOwrKTnKE8xMTHcdttt5y1v2rRpoWPJyclcf/31nu916tTh559/LlCnVq1aAPj6+lK9enWOHj1Ks2bNSilqqQgCAwNxOp00aNCA48ddWyicOHECh8NR7LbTpk0DYPny5Z52eWX9+/enf//+pKWlAeBwOC56zrzy7du3A/D0008XKx4RERERESl/pTGD/g6wFbjWWjvCWrsVwFp7BNes+m+Kr2/hexq1atUqkJD//PPP1K5du0Cd5ORkAM6ePUtaWlqhcrn8NWnShFOnTpGUlATA6dOn2bdvH1dfffVFWkKzZs3YtWuX5/sPP/zgaVdUWWhoKFdcccVFz3mhfkVEREREpGIpjQS9B7DAWnsKwBhTyRhTFcBa+2Ep9F/h3XXXXXz44YckJiZy/Phx/vGPf9CrV68Cdd5++23OnDnDP/7xD8LDw2nSpImXopWyEhwcTLdu3Zg8eTJZWVlMnTqVsLAwIiMLvuggKiqKMWPGFDgWHR3NsmXLiIuL4/Dhw8yaNcuzG/stt9xCRkYGCxYs4OTJk0yZMoV77rmnWOe85557mDZtGunp6Xz55Zds2LCBPn36lMPVEBERERGRS1XiJe7Af3C9tzzvwdaqwBqgYyn0fVno168fsbGxdOjQgbNnz9K/f3/++Mc/FqhTo0YNatasSXh4OAsWLMAY46VopSzNnDmTQYMGUb16dZo3b86iRYsK/a4TEhIIDw8vcOz6669n/PjxdOnShTNnzjB8+HBPgl61alWWLFnCsGHDGDp0KLfffjuTJ08u1jlHjhzJrl27aNCgAdWqVWPWrFk0aNCgbC+CiIiIiIj8KsZaW7IOjNlmrW19sWNloV27dnbz5s1lfRoRERERKWVPLthK15Z16Ni4Bt2m/pctf77sXwIkInJexpgt1tp2F6tXGkvcTxpj2uQ7cVvgVCn0KyIiIiIiIvK7URpL3J8BlhhjjgAGqAMMKIV+RURERERERH43SpygW2u/M8Y0B/K2ht5trT1T0n5FREREREREfk9KYwYdoD0Q7u6vjTEGa+3cUupbRERERERE5DevxAm6MeZDoDGwDchxH7aAEnQRERERERGRYiqNGfR2QAtb0u3gRURERERERH7HSmMX9524Noa7rCUkJGCMKfD+8kGDBjFhwgTvBeUFu3fvpnv37lSrVq3Qu7qLa+LEiRhjiI+P9xxLSkoiKiqKqlWr0qZNG3bu3OkpO3z4MHfeeSchISFcddVVLFiwwFO2atUq2rdvT3BwMA0aNGDSpEm/emwX0qlTJ+bPn1/ifi40zgtZt24dkZGRBAcHU7duXUaNGkVOjmtByrFjx7jxxhsJDQ2lRo0a9O7dm59++snT9uuvv+a6667D4XDQtm1btm3bVqDvt956i7CwMBwOB23atCE3N7fE4xQRERERkdJXGgl6TeAHY8znxphP8j6l0K9XLF26lF9++cXbYXiNr68v9957L2+88cavap+QkMDatWsLHR82bBgRERGkpqYyYMAABgz4/xv9P/nkkzRo0IBjx47xwQcf8PDDD3PgwAEAnE4nf/3rX0lJSWHjxo3MmzePefPm/brBncd3333HgQMHuOeee0rc14XGeSHNmzdn1apVZGRksGfPHuLi4pg+fToAgYGBvPfee6SkpHD06FEiIyMZMWIEANnZ2URHRzNq1CgyMjLo06cP0dHR5C1oWbhwIX/7299Yvnw5J06c4IMPPsAYU+JxioiIiIhI6SuNBH0C0Bt4Ffhbvs9l6bbbbmPRokXeDsNrGjduzJAhQ2jYsOGvav/MM8/wyiuvFDiWmZnJmjVrGDt2LAEBAYwcOZKDBw+yY8cOAL788kuGDx9O5cqV6dKlC61bt+bTTz8F4J577uG2227Dz8+PevXq0b17dzZs2HBJMYWHh5OQkHDe8rfeeosRI0ZQuXLlSxvsOS42zgupW7cuDRo0wBjD2bNnyc3N5ccffwQgICCAli1b4uPjQ05ODjk5OZ6y3bt3c/LkSQYPHoyPjw+jR48mPj6e2NhYAGbMmMELL7xA27ZtMcbQqlUrJegiIiIiIhVUiRN0a+2XQAJQ2f3zd8DWkvbrLY899hgzZ870dhiXpVWrVuHn58eNN95Y4Hh8fDz+/v44HA46d+7M/v37ady4Mbt37wagqO0L8i+Pz2/jxo20atWq1GI+fPgwK1euZPjw4SXu62LjvJjExERCQkIIDQ1l06ZNDBkypEB5q1atCAgIYPLkyZ4Z9KKuXf7HC7Zv387Ro0dp3LgxYWFhjB8/vmSDFBERERGRMlPiBN0Y8yiwFMjLausBMSXt11uuv/56srKy2L59u7dDuaxkZ2fz/PPP89prrxUqO3nyJA6HgxMnThAXF0daWhpBQUE4nU4AoqKimDlzJtnZ2axdu5bvv/++yMcMpk+fTlZWVqHEtSiTJ0+mWrVqVKtWjcTERFq1akW1atXo2bNnoT4HDBhAjRo1ft3AL2GcFxMWFkZGRgb79u1j1KhRNG3atED5jh07OH78OJMnT6Zz586Aa2l8YGAgc+fO5cyZM0yZMgUfHx/P9cvMzGTVqlVs2LCBr7/+mjlz5rB8+fISj1VEREREREpfaSxxHwHcCGQCWGv3ArVKoV+vGT58uGbRL9GUKVPo2bNnkRvLBQYG4nQ6adCgAcePH6dDhw6cOHECh8MBwDvvvENSUhL169dn8uTJ3HXXXQQHBxfoY9WqVbzxxht88skn+Pn5XTSesWPHkp6eTnp6OmFhYezYsYP09HRWrlzpqXPq1Cn++c9/8swzzxRqP3/+fBwOBw6Hg8cee6xY1+Bi4yyuRo0a0bx5c88seX4hISE88MAD9OzZk5ycHPz8/FiyZAlvv/02devWJTk5mUaNGnmuX9WqVXnooYeoVasWYWFh9O3bl/Xr119SPCIiIiIiUj5KI0HPttaezvtijPHF9R70y9bAgQOJiYkp9synuDZae+211zDGeJ5xbtq0KTExMTRp0oRTp06RlJQEwOnTp9m3bx9XX301APXr1+ezzz4jJSWFf//73xw4cIBrr73W0/e3337LsGHDWLlyJWFhYaUW89y5c2nTpg0tWrQoVDZw4ECcTidOp5MZM2YUq7+LjfNSWGs9z5EXVZaUlERaWhrg2oF+y5YtHDt2jJdffpmDBw96HgNo3LhxoWfO9UZEEREREZGKqTQS9C+NMS8CAcaY24ElwKel0K/XBAcH06NHDz7//HNvh1LurLVkZWVx5swZz8+nT58uUGfMmDFERUUVOBYTE4O11vMB2Lt3L7179yY4OJhu3boxefJksrKymDp1KmFhYURGRgJw6NAhjh8/TnZ2Nm+99RaHDh2iT58+gGtZd3R0NIsXLyYiIqJUx/n2228zcuTIUuvzYuPMExUVxZgxYwoc++ijj9ixYwe5ubkkJiYyZcoUOnToALheo7ZmzRqys7NxOp28+OKLNGnShJo1awKwa9cusrOzSUtL48knn6RLly6elQx9+vTh/fffJzU1lSNHjrB8+XK6dOlSamMWEREREZHSUxoJ+lggBYgFhgOrgHGl0K9XDR8+nKysLG+HUe4OHjxIQEAAd955J4mJiQQEBNC1a9cCdZKTky+4K3pRZs6cSWxsLNWrV2fhwoUsWrTIM7O7e/duIiIiqF69OosWLWL16tWeJdpTp04lOTmZrl27epac33HHHZd07oSEhEJL7/NuvnTr1u2S+rqYC40zfzzJyckFjqWmptK/f3+Cg4Np3749rVq1YsqUKQDk5ubywgsvUKtWLerXr8/hw4dZsWKFp+3y5cupU6cOYWFhnD17lg8//NBTNmbMGCIiImjUqBFt27blgQce8Nz8EBERERGRisVczstd27VrZzdv3uztMOQy1L17d3r37l3s58tFRESkdD25YCtdW9ahY+MadJv6X7b8+XZvhyQiUmaMMVuste0uVs+3FE50gCKeObfWNipp3yJlITs7m44dO/LAAw94OxQRERERERGPEifoQP67AP5ANBBaCv2KlAk/Pz9eeuklb4chIiIiIiJSQImfQbfWpub7HLbWvgX0KIXYRERERERERH43SmOJe5t8XyvhmlEvjZl5ERERERERkd+N0kik/5bv57NAAnBPKfQrIiIiIiIi8rtRGkvcu+T73G6tfdRau7s0gpPyN3XqVBo2bEhwcDBNmjRh9uzZl9zHxIkTMcYQHx/vOZaUlERUVBRVq1alTZs27Ny501P2+uuv07RpU4KCgoiIiCjwCjFrLc8++yy1a9emRo0a3HfffWRmZpZskEXo1KkT8+fPL3E/FxrnxUybNo3atWtTvXp1XnjhhQJlxhgCAwM9r5p77733PGUrVqygQ4cO+Pn5MWTIkALtjh07Rp8+fahevTpXXnklEyZMKMnwRERERESkDJU4QTfGjLrQpzSClPLTo0cPtmzZQmZmJmvXrmXcuHFs37692O0TEhJYu3ZtoePDhg0jIiKC1NRUBgwYwIABAzxlvr6+LFu2jIyMDGbOnMngwYPZv38/AIsWLWLx4sVs3bqVxMREjh8/zsSJE0s+0Hy+++47Dhw4wD33lHzhx4XGeSGbNm1iwoQJfPHFF+zcuZOFCxeyePHiAnW2b9+O0+nE6XQybNgwz/GQkBCee+45HnnkkUL9jhs3jjNnznDkyBE2btzIe++9x2effVayQYqIiIiISJkocYKO65nzx4F67s9jQBsgyP25LHzwwQe0a9eOBg0aMHDgQPr06UOdOnWIjY0FYMiQIYwbN85TPyoqivfff99b4ZaZZs2aERrq2oT/9OnT5ObmsmvXrmK3f+aZZ3jllVcKHMvMzGTNmjWMHTuWgIAARo4cycGDB9mxYwcAo0aNIjIykkqVKnHjjTfSqFEjtmzZAsDBgwfp2LEj9erVIzAwkJ49e/LDDz9c0pjCw8NJSEg4b/lbb73FiBEjqFy58iX1e66LjfNCli5dSt++fWnZsiX16tVj6NChLFy4sFjnjYqKom/fvp7fW34HDx6kZ8+eBAQEEBYWRocOHS75+omIiIiISPkojQS9PtDGWjvaWjsaaAuEWWtftta+XAr9l5sqVaoQGxvLxx9/zKOPPsojjzzCokWLvB1WuVuwYAGBgYE0a9aM4OBgunbtWqx2q1atws/PjxtvvLHA8fj4ePz9/XE4HHTu3Jn9+/fTuHFjdu8u/CREWloae/bsISIiAoB7772XAwcOcOjQIZxOJytXrqRHj9J7ScDhw4dZuXIlw4cPL3FflzLOc/0/9u49Tquy3v//6zMHYGCAAQREBFFQFElBUTy20VIpSS3zkIeyMjVr7525t6aVun9aaVmIu76ZbvKApzTPpzxlmnnkoCgoIsr5KDLDcYCZuX5/3PeMAwwwyAz3Pfh6Ph73Y9Z9rXVd67PW3BbvudZa95QpU9hzzz0ZNWoUF154IQMGDNig3+c//3l69OjBt7/9bSoqKhpV0/nnn89TTz3FihUrmD59OuPGjeOoo476VMcnSZIkqXk1RUDvDqyp935Ntq3F6du3L2VlZeywww7079+fXXfdlQULFuS6rG3utNNOY/ny5bz44oucffbZdOzYcbN9Vq9ezUUXXcQ111yzwboVK1ZQWlrKsmXLmDRpEkuWLKF9+/YsX758g23PPfdcvvWtb7HXXnsBsOOOO3LggQeyyy670LFjRwoLC9e5vHtjrr76asrKyigrK2PmzJnss88+lJWVMWLEiHW2+8Mf/sApp5xCly5dNjvm5mzJcW6s77Rp05g6deoG/V566SVmzpzJ+PHjmTt3Lv/xH//RqJoGDx5MRUUFHTt2ZNddd+Xss89m0KBBn/oYJUmSJDWfpgjotwGvRcQVEXEF8CpwaxOMu80VFhYCmXuiCwsLKSoqoqqqKsdV5UZEcOihhzJt2jRuuummzW5/7bXXMmLECPr06bPBunbt2rF8+XJ69erFxx9/zMEHH8yyZcsoLS1dZ7tLL72UxYsXM2rUqLq2K6+8kokTJ7Jw4UIqKiro0KFDo8LpT3xqivIAACAASURBVH7yE8rLyykvL6d3795MnDiR8vJyHn300bptVq1axU033cSPfvSjDfrfcccddQ9kO++88za7vy05zk31vf7663n44Yc36HfwwQdTVFREjx49uPLKK3nkkUcaVdOpp57Kvvvuy4oVK5gxYwZ33nnnBve2S5IkScoPTfEU918A3waWZF/fTin9cmvHzRcpJQDatGmzTlhvjieJ56OUUt19+Jvy+uuvc8011xARRAQAu+++Ow8++CD9+vVj1apVzJ49G8jc2z5t2jT69+9f13/kyJE8+eSTPPjgg7Rq1aquffz48Zx44onssMMOlJaWcuaZZzb4ELpP47bbbmO//fZjwIABG6w7/fTT6x7IdsMNNzRqvMYc58bsscce69zrP3ny5I32KygoqPtcbs748eM566yzaN26Nb179+bYY49tsvMnSZIkqWk1xQw6QFtgaUppFDA7InZtonHzxu67787rr79OSon33nuvUaG1Jbr++uuZPXs2KSVefvll7r77bg4++OB1trn44osZNmzYOm0PPvggKaW6F8DUqVM54YQT6NChA8cccwxXX301lZWVjBw5kt69e/O5z30OgFtvvZUbbriBJ554gvbt132u4ODBg3nggQdYsmQJlZWV3HPPPQ0G6i2VUmLUqFFccMEFWz1Wrc0dZ61hw4Zx8cUXr9N20kkncf/99zNp0iTmzJnD6NGj654A/9ZbbzFhwgSqq6tZvHgxV1xxBV/5ylfq+lZXV1NZWUl1dXXdcu0fkwYPHsyYMWNYu3Yt8+fP58knn2yS8ydJkiSp6TXF16xdDlwM1H5xczFw+9aOm2++/e1vs2bNGvbee2+uvPJK9t9//1yX1CzefPNNhg4dSmlpKaeddhoXX3wxp59++jrbLFiwYJNPRW/In/70J9566y06derE3XffzV/+8pe6mfYrrriC6dOns9tuu9VdVv7LX2Yuwrj00kvZdddd6d+/PzvttBMff/wx//u//7tF+54+ffoGl94/+eSTABxzzDFbNNbmbOo469ez/rMNhg4dyuWXX84RRxzBwIEDOfnkk+sC+qJFizjxxBPp0KEDAwYMoHv37lx//fV1fceMGUNJSQlXX301t99+OyUlJVx11VVA5tsJ3njjDbp27cqgQYP4t3/7N84///wmPWZJkiRJTSMae6nsRgeIeAMYDIxPKQ3Otk1MKe2zmX6FwFhgTkppRHbW/W6gCzAOODOltGZTYwwZMiSNHTt2q+rXZ9Pw4cM54YQTGn1/uSRJalo/vHM8R++9I4f07cIxI19g3M/9lhFJ26+IGJdSGrK57ZriEvc1KZPyU3bH7RrZ7z+Bd+q9vwYYmVLqR+Ze9u82QW3SBlavXs0hhxzCN7/5zVyXIkmSJEl1miKg3xMRfwLKIuJ7wDPAJh/7HRE7A8cC/5d9H8CRwF+zm9wKnNAEtUkbaN26NZdddhlt27bNdSmSJEmSVKdoawdIKV0bEUcBS4H+wGUppac30+064CKg9olgXYDylFLtY9JnAz0b6hgR5wDnAPTu3Xsrq5ckSZIkKT9sVUDP3kf+TErpCGBzoby2zwhgYUppXEQM29J9ppRuBG6EzD3oW9p/a5133nn07NmTn//859t615IkSZKk7dhWXeKeUqoGaiKi4xZ0OxQ4LiKmk3ko3JHAKDKXyNf+wWBnYM7W1NZcbrjhhu06nK9du5bvfve7tG/fnt69e3PPPfc0qt/UqVMZPnw4nTt3plu3bnzzm9+koqICgJkzZ9Y9nb32FRHcd999m93n448/zgEHHECHDh3o1atX3dPJm9phhx3GHXfcsdXjzJ49m2HDhtG2bVv2228/3n777Ub3vf766+nevTudOnXikksuWWfd7bffTt++fenQoQOHH344kyZNqls3bNgw2rRpU3duzzzzzLp1H330EV/96lfp1KkTO+20E1dcccVWH6MkSZKk5tEU96AvB96KiNERcX3ta2Mbp5QuSSntnFLqA5wK/D2ldDrwHPD17GbfAh5qgtq0hUaOHMlbb73FrFmzuO222/jOd77DrFmzNttv6dKlfOMb3+DDDz9k+vTpVFZWcuGFFwKZWxGWL19e95owYQKlpaUMHz58s/tcvnw5v/rVr1i0aBGvvPIKt99+O7ff3rTf4vf666/z4YcfcvLJJ2/1WOeccw4DBw5k8eLFnHLKKXVflbY5r776KldccQV///vfefvtt7n77rvr/lAxe/ZsvvOd7zB69GjKy8s54ogjOOuss9bp//vf/77u/I4ZM6au/Wc/+xlr165l7ty5vPLKK9x444089thjW32ckiRJkppeUwT0+4GfAy+Q+Xq02teWuhj4cUS8T+ae9NFNUNsWmTp1KocddhgdOnSgR48eXHPNNXXrHn30UUpLSykuLuZnP/vZOv1SSlx88cV06dKFgQMHcvrpp3PGGWcAcNZZZ3HkkUfSrVs3LrvsMgYMGMD+++/PypUrAfj1r3/NbrvtRtu2benXr1+jZ6yby7333st//Md/UFZWxrBhwzj44IN54IEHNttv//3351vf+hYdO3akbdu2nHbaabz88ssNbnvzzTfz1a9+lXbt2m12nyeffDJf/OIXad26NT179mT48OEbHXdj+vTps8nvbb/uuuv4wQ9+QHFx8RaNu76lS5fy1FNP8ZOf/ISSkhIuuOACZsyYwcSJEzfb969//Stf+9rX2HvvvenZsydnn302d999NwBz5syhc+fODBs2jIKCAk466SQmT57cqJpmzJjBiBEjKCkpoXfv3hx88MGN7itJkiRp2/rUAT0iegOklG5t6NWYMVJK/0gpjcguf5BSOjCl1C+ldFJKafWnre3Tuuyyy9h3330pLy9n8uTJHHrooXXrRowYwfLlyzn99NM36Hffffdx7733MnHiRB566CGeeuqpddYPHTqUMWPGcM011/Diiy9SUlLCSy+9BECnTp144oknWLFiBddddx1nnnkmixYtat4D3YQpU6aw5557csYZZ3D33XczYMAApkyZssXjvPzyy+yzzz4btNfU1HDbbbet8xVnW7LPV155pcFxP605c+bw6KOPcu655271WO+//37dpeaHH344H3zwAX379m3U+as9B6NGjeLCCy9c5xwMGTKE/v3789xzz1FdXc29997Lscceu07/Sy65hB122IGjjjqKd9755NsLzz//fJ566ilWrFjB9OnTGTduHEcd5ffMSpIkSfloa2bQH6xdiIj7mqCWnCsoKGD+/PnMmzePTp06cdhhhzWq38MPP8yZZ55Jz5496du3L8cff/w66/v160ffvn3p1q0bnTt3pk+fPixYsACA733ve/Tv35+IYMSIEZSVla0TsLa1FStWUFpayltvvcWcOXNo3749y5cv36Ixxo0bx80338yVV165wbqnnnqKiODII4/c4n3+4Q9/oLKycoPLuxty9dVXU1ZWRllZGTNnzmSfffahrKyMESNGbDDmKaecQpcuXbboGBtSexzLli1j0qRJLFmypNHnr7bvtGnTmDp16jr9CgsLOeOMM/jKV75C69atueOOO7j22mvr+l577bXMnDmTGTNmsP/++3Pcccexdu1aAAYPHkxFRQUdO3Zk11135eyzz2bQoEFbfaySJEmSmt7WBPSot7zb1haSD37961/ToUMHhgwZQr9+/bj33nsb1W/BggXsuOOOde+7deu2zvrCwkKKioooLCwEoKioiKqqzDfKjRkzhn333ZfOnTtTVlbGokWL6sJVLrRr147ly5fz5ptvcuGFF7Js2TJKS0sb3X/69OmceOKJ3HbbbfTr12+D9TfffDOnn346BQWffPQas8/HH3+c3/zmNzz88MO0bt16s3X85Cc/oby8nPLycnr37s3EiRMpLy/n0Ucfrdtm1apV3HTTTfzoRz/aoP8dd9xR99C18847r1HHXnscvXr14uOPP+bggw9u9Pmr7Xv99dfz8MMPr9PvmWee4ac//SmvvvoqlZWV/Pu//ztf/vKXSSnzJQZDhgyhpKSEdu3a8ctf/pL58+fX/ZHn1FNPZd9992XFihXMmDGDO++8M+e3UUiSJElq2NYE9LSR5RarZ8+e3HzzzcydO5eLLrqIc845p1H9unfvzvz58+ve186Ob0pKiRkzZnD22Wfzxz/+kcWLF1NeXk5ZWVld8MqFPfbYg3fffbfu/eTJk+nfv3+j+i5cuJBjjjmGq666qu4BcPUtWbKEhx56aJ3L2xuzz5deeolzzjmHRx99lN69e2/pIW3Ubbfdxn777ceAAQM2WHf66afXPXTthhtuaNR4/fr1Y9WqVcyePRuANWvWMG3atEadv02dg/Hjx3PYYYex9957U1RUxPe+9z0mTZq0zmeuVkQQEXWfofHjx3PWWWfRunVrevfuzbHHHsuzzz7bqOORJEmStG1tTUDfNyKWRsQyYJ/s8tKIWBYRS5uqwG3pgQceYO7cuXUhp3379o3qd9xxx3H77bczZ84cpk2bxiOPPNKofitWrCAi6NatG9XV1Vx77bWUl5dvzSFstZNPPpnrr7+e8vJynn/+eV5++WW++tWvrrPNxRdfzLBhw9Zpq6ioYPjw4Xz/+9+ve0De+u68804GDhy4QSDe1D4nTpzISSedxD333MPAgQOb7DhTSowaNYoLLrigycbs0KEDxxxzDFdffTWVlZWMHDmS3r1787nPfW6d7YYNG8bFF1+8TttJJ53E/fffz6RJk5gzZw6jR4+uewL84MGD+de//sWUKVNIKTFmzBi6dOlCt27dKC8v54knnmD16tWsXr2a//mf/6Fbt25153jw4MGMGTOGtWvXMn/+fJ588skG/yAhSZIkKfeKNr9Jw1JKhU1ZSD4YO3Ys559/PsuWLaNv377ceusnz7o7+uijeemll1i9ejURwXXXXcfXv/51brnlFk488URee+019tlnH3baaSe+9KUvNWoWfMCAAVx44YUceOCBFBYWcv755zfpDPGnccEFF/Duu+/Sq1cvysrKGD16NL169VpnmwULFmzwVPQHH3yQCRMm8N57763zlPv691/ffPPNG8yeb26fI0eOZMGCBRx99NF12x9++OE88cQTjT6mhp7g/uSTTwJwzDHHNHqcxvjTn/7EGWecQadOndhzzz35y1/+QkSss8306dPp06fPOm1Dhw7l8ssv54gjjmDt2rWce+65dQH9qKOO4r/+6784+uijWbJkCf379+eBBx6gsLCQtWvXcumllzJ16lSKi4s54IADeOSRR+qeSH/LLbdw/vnn07VrV9q0acNJJ53E+eef36THLEmSJKlpRC4vp95aQ4YMSWPHjs11GRv44Q9/SNu2bfn1r3+d61K0EcOHD+eEE05o9P3lkiSpaf3wzvEcvfeOHNK3C8eMfIFxP/dbRiRtvyJiXEppyOa2a4rvQf/Mq6mp4a677qKyspIZM2bw0EMP8cUvfjHXZWkjVq9ezSGHHNLgbL4kSZIk5YoBvYmMGjWKrl27csghh3Duueeuc0m28kvr1q257LLLaNu2ba5LkSTpM6equobnpiykcm11rkuRpLzzqe9B1ycKCgp45ZVXcl2GJElS3nvsrXn8591vAHDcoJ45rkaS8osz6JIkSdpm1la33OcfSVJzM6BLkiQpZ1oXFbC2uoZL7n8r16VIUs4Z0CVJkpQz7dsUc8MZ+/PegmW5LkWScs6ALkmSpJxqVeQ/SSUJDOiSJEmSJOUFA7okSZIkSXnAgC5JkiRJUh4woEuSJEmSlAcM6JIkSZIk5QEDuiRJkiRJecCALkmSJElSHjCgS5IkSZKUBwzokiRJkiTlgaJcFyBJkqTPlo4lxexUVkL/7u1zXYok5RUDuiRJkrapL+zVjd+dPCjXZUhS3vESd0mSJEmS8oABXZIkSZKkPGBAlyRJkiQpDxjQJUmSJEnKAwZ0SZIkSZLygAFdkiRJkqQ8YECXJEmSJCkPGNAlSZIkScoDBnRJkiRJkvKAAV2SJEmSpDxgQJckSZIkKQ8Y0CVJkiRJygMGdEmSJEmS8oABXZIkSZKkPLDNA3pE9IqI5yJickRMioj/zLZ3joinI2Jq9menbV2bJEmSJEm5kosZ9CrgwpTSAOAg4AcRMQD4CfBsSml34Nnse0mSJEmSPhO2eUBPKc1LKY3PLi8D3gF6AscDt2Y3uxU4YVvXJkmSJElSruT0HvSI6AMMBl4FuqeU5mVXzQe6b6TPORExNiLGLlq0aJvUKUmSJElSc8tZQI+IUuA+4EcppaX116WUEpAa6pdSujGlNCSlNKRr167boFJJkiRJkppfTgJ6RBSTCed3pJTuzzYviIge2fU9gIW5qE2SJEmSpFzIxVPcAxgNvJNS+l29VQ8D38oufwt4aFvXJkmSJElSrhTlYJ+HAmcCb0XEG9m2S4GrgXsi4rvADODkHNQmSZIkSVJObPOAnlJ6EYiNrP7CtqxFkiRJkqR8kdOnuEuSJEmSpAwDuiRJkiRJecCALkmSJElSHjCgS5IkSZKUBwzokiRJkiTlAQO6JEmScioCZixeySsfLM51KZKUUwZ0SZIk5dSgXp04eu/u3DN2Vq5LkaScMqBLkiQppwoLgv16d8p1GZKUcwZ0SZIkSZLygAFdkiRJkqQ8YECXJEmSJCkPGNAlSZIkScoDBnRJkiRJkvKAAV2SJEmSpDxgQJckSWpBPli0nJVrqnJdxqf28YrVuS5BkvKWAV2SJKmFeGfeUo787fP86vF3c13KpzJ57lJGPj2VA/p0znUpkpSXDOiSJEktRO3M+YrVLXMGfeWaKgbs1IFvHNg716VIUl4yoEuSJEmSlAcM6JIkSZIk5QEDuiRJkppdSonn31uU6zIkKa8Z0CVJktTsFq9Yw03//IB/P7JfrkuRpLxlQJckSVKTqVi1lqWVaxtc165VEcP6d9to35Wrq6mpSc1VmiTlPQO6JEmSmsyXrnuB43//ry3u169bKS9N+4g7Xp3RDFVJUstQlOsCJEmStP2YW1FJq8ItnwMa1KuMbx7ch/KVDc++S9JngTPokiRJkiTlAQO6JEmSJEl5wIAuSZKkvDJh5hL+9va8XJchSducAV2SJEl55ecPvc15t49n2UaeBi9J2ysDuiRJkvLGmuqaXJcgSTnjU9wlSZKUF/p2a8dFf53I2mq/C13SZ5Mz6JIkScoLXx28M0fvvWOuy5CknDGgS5Ik5dAHi5bznVteZ9yMJY3u869pH23R9i1JcUHkugRJyhkDuiRJUg5NmFnO399dyL/e/6hR2++7c0cO2q1Lo7dvaX567ADuPucgSlt7J6akzx4DuiRJUp6bMn8Z/5y6iLfnLKWwIOjduW2uS9qkmpR44q15zC1fVdf29pyKRvXt2r41B+3WZdPj12w4fq6Mm/ExE2YuoXJtNQ9OmMNSnzwvaSsY0CVJkvLcRX99kzNHv8blD0/itKG75LqczaqqSXz/jvH87un36tr+/a4JfPvQPk0y/oRZ5Xz/jvGMrDd+rnzjxlf5zi2v8+jEefzoL28w5uUZuS5JUgvmtUOSJGm7l1LmqeARLfP+5prsQ82LCoKv778zv31qCimte1ybO8atPgcpwdpVsHoZqWoVVK0hqtdA9Wqo+uTnMQWv04q1tKKKAxa3h7HjSSnx9ZpJnN1mD3jtH5mxSJBqMssRUFC0zutLTKbo3WXQug2poJAoKIbCYihuS+vFa9gl5lO6pjWsWkIqKoHCVkRB0849beqc1a5bU11DwdrMVQOQmd2XpE/LgC5JkrZ7598xntenf8w/LzqSklaFuS5nq+3cqYRL7n+LMa9MB4I7zh7Kcb9/kXM+vxsXHt1/ww41NYx8bByPvPYud505gB3brIXKpbB6KVRWwOpl2eVs2+pl2eWK7M/s+poqADYV8f/Uqt6b+cCjme0vLwCeavwx/iaAB2lwfwOB51sDU4FrMuurKaCguIQoLoHittCqLbTuAK3bf/Jq03Hd963bZ7ept12bDtCqFCK467VZ/PTBt7j33IMZ0qfzOjX86ol3G5wtH/nMexzSrwv779J5g3WStDkGdEmStN2bU76Kj5avYU1VTcsL6FVr6FBTQe9YQKdYBR++wCmlS1mzy9u8P2su7VlJyT8e53I+ZK+3gYUF9QJ2NnSvWcaPgR8XAnduZD9RkAmqbTpA62yQ7dATuu6VbWtfF2RHPT+TaUuqufrkIbQtaQtFraCwNRS15us3jeOW7x1Gadu2/OnFWRQXFfKPqYu45EsD2GunMiAyM+aRne2OyMyi11Rn/gBQsxZqqhgx6nn+cvb+/Ou9Bfz+mXc488CenDSoGxXLlvGLB8dx4j5dmLlgMSft04VbX3iHimVLOe+AHrRKqzMz/WuWZ89DOVTMyv6RYVmmfXOiAFq359jUlsHFxXR7pCt06Vrv/HRgwLTlfK2mimUFbamkHV0+XsMP965hSXUJiz76CHqVQRPP6Eva/uVdQI+I4cAooBD4v5TS1TkuSdJ2oLom8cvH36Fv11JOG9p7i/s/9+5C7p8why7tWvHTY/eiuNB/dEmb87/PTiUCfnjk7ls9TkFB8IMj+tW13fXaTKYuWM6lX96TovX+e6yqruGqx95hrx7t2XPHDvz5Xx9us4eJVaxay1WPTubLA3fkiH4dP5l5rg2Lda+ldcv7zprHb4tn0fPVKqaOX8Nu7ath9VJWLiunpGYFRTWruQOgdXYnt2Z+nAlQnFle824rjipsw6pl7ZhfU0ZFTQmVBTsQbXajoqiEnXvsyIQF1bw6r4oBu/bkawcN4JZxizlwzz4M3WvXTOgsbktFZRVXPTqZnp1KmP7RCi4avic7lZXUHd8bs8q55V8f8s/KPiyuWcNVex7NrJVrufapKQRw0fA9mcJianbYA9oUs6LNGp6ctIB5FSVUt+sOpR0bfS6nM5Xy9rsz+v0VTEx9md9hD+5c2JpXP1zMCwxlaK8BjJryHrOX78zNq/qytKqKd5bsyBF7duPkIb14atJ8nn1nIT8/ZcC6T4SvqV7vd1H7quDVd6aztOJjvrBrGwpWL2PWh7OYu3IBXaMQls6F1e9C5VJSZQUnpGpOKK5X8EvwhU+KJz0SrC5sx5rCUpbRlrJOXWjXoXNdwK//88VZa5ixooiV0Y7eO+3IMfvtkZntb1W6zUP+k5Pm89y7C/nZiAENPkl/TVUNv3hsMgN7duSkIb2afP8jn36PNsWFfH9Y3yYZb8biFfz2qfc49992Y++dGv/5k3IlrwJ6RBQCfwCOAmYDr0fEwymlybmtTFJLt3JNFaNf/JCBPTt8qoD+yMS5dGhTxF/HzeaHR/Zjh9LWm+8kfcZd9+xUCiO2OqCPfOY9igoLNgjoE2dX8J9f2J2ObdcNMEsrq7jlpekM6lXG53ffgYffmE0rqiilili5GKpqoHoNVK+FqtWfLNfeT13bXlUJa1fCmpWZGdm19X820LZmJSWVK7h0+VLaT6oEqjZ/cFFI76J2tCloxfLKEioqS1jRqQdVJTvzt0UrWEZbenTvxh69e1JSWkZNqw7sunMPaN2B+Wta8e7HQcdOnZlZUcUe3dsz6pmp/G3SfL518C7ssWN7fvrA2wB8bceerOlcw257tOOXL3xA78MH89t3XudrJW0ZOrRHXTnvL1zGveNm170/cq/uHFcvoD/7zgIefGPuOofw6ocfM7d8FVU1ibHrfTf7dw/bjb7dSmnbqoi9d+qw+fOxnvcWLGPhstV8ff+dAbjztRm8PWcprYsKGLFvD56aPJ9Rz07l+8P6ctSA7jw/ZRF/eX0WJw/pxUNvzuWxifM48+BdGNizXjArKISSssxrPf/f3//JpLlLeeuMo2nfppgnnnyXP8yaxv8ePpiv7LtT3XbX/u0dFixewpz5C1j40SI6sJLdO9Zw6ZE9eWzsFAZ3K2CX0mrueXESHWIl7VnJgJIa2sU8+GjKJ1c3ZG8ZOCz7yvwSgBdq38QnYb5VKdRevl/cpt5yCRSVrPu+uM2664paQ2GrzKuo1SfLhcUbLD88YRaPvb2QMw5a77xlla9cw60vz+CAPp2aJaCPenYqHdoUNVlAf+WDxTz85lz679jegK4WIa8COnAg8H5K6QOAiLgbOB5o2QH9nUeYO/HvrK2u+VTdo0meNbK1g2xl/7R1/WOr64dcH0NT2OrzsNXHkPvfQ3zKY6iqSfyqaDHtywuZc+stW9x/xILl7LpDOw4oKGf5PXezeitm0D/tMXwi1/3Jg/+mm+IYcjtAk/zvWp5/lkYVZsLa/NG3bnSbSIlI1USqgVRdt1y/7d7ipRRSw9JRVxGpGlINo5aupKZVFYV/KGYlNev0aVdTzZutK2m1qJpWi6r4cZvqT3b4+09/tDUFxVQXtqGmsITqojZUF5ZQU9iG6qISqgs7U912JyqKi3i1opIOZZ0YuOvOVBWVUl3cjqriUqqKSqkqLqW63nJNYRte+uBj/vyvD+v287P+e1FVk7j6nXcB+PGAPTjuCxv+kWNHYMc+meXB2bbdu5fyt0mwf5/OHNq3S11An1deSVVNDUfvvSOtCguYMLO8rv2ZyQvqxnx/0bqXf789p4K2xZ/cEvDBRyvWWf+PKYuYNLeCXp3bsqaqhrfnVLC25pN/73RsW8zxg3pu0Xmub8LMcjq3a8WOHdowbdFylq7KBNqigqB1USG7d2vPk5MWMKBHB/br3YmamsSjE+fyzOQFLKioBDJ/QJifXd6c2q9I+8eURZQUFzL9o5V156Gk3nmYvngVe/XoyjuLq5mWMjfdF3bsRNmBh/DS1PHMbFvCHl3b8z9Vb9b1OaHHTozY55OQT0oUVFdStHYZf3p6AisqltA+VtKvQw0nDexA0ZplFK2t96paQUF1JYUrKymorqCwqpKC6lUUVq+msHoVBVWVFKRG/FFoM/4AXNe6EEYXs7awFTUFxaSCQohCUhTSNgV/b7WWwoVFqv82xwAAIABJREFULP9dCSkKSFFIioK6bTZsW/99kHlyQJDqHjCQeT+qeAlRE8z78+3Zhw8EqfYpBOv0i7p+9Zdrb6Go7dOnopIripbTd2IpM+e0y46xrWyf+0rb8rA2c1zdO7ShzYjt64LrSHkQOmpFxNeB4Smls7PvzwSGppR+WG+bc4Bzsm/7A1O2eaGfzg7AR7kuQmokP69qKfysqqXws6qWxM+rWoqW9FndJaXUdXMb5dsM+mallG4Ebsx1HVsqIsamlIbkug6pMfy8qqXws6qWws+qWhI/r2optsfPar495WgOUP9mlp2zbZIkSZIkbdfyLaC/DuweEbtGRCvgVODhHNckSZIkSVKzy6tL3FNKVRHxQ+BJMl+z9ueU0qQcl9VUWtxl+fpM8/OqlsLPqloKP6tqSfy8qqXY7j6refWQOEmSJEmSPqvy7RJ3SZIkSZI+kwzokiRJkiTlAQP6NhARwyNiSkS8HxE/yXU90sZExJ8jYmFEvJ3rWqRNiYheEfFcREyOiEkR8Z+5rklqSES0iYjXIuLN7Gf1f3Jdk7QpEVEYERMi4tFc1yJtSkRMj4i3IuKNiBib63qaivegN7OIKATeA44CZpN5Uv03UkqTc1qY1ICI+DywHLgtpTQw1/VIGxMRPYAeKaXxEdEeGAec4P+2Kt9ERADtUkrLI6IYeBH4z5TSKzkuTWpQRPwYGAJ0SCmNyHU90sZExHRgSErpo1zX0pScQW9+BwLvp5Q+SCmtAe4Gjs9xTVKDUkovAB/nug5pc1JK81JK47PLy4B3gJ65rUraUMpYnn1bnH05O6K8FBE7A8cC/5frWqTPKgN68+sJzKr3fjb+I1KSmkxE9AEGA6/mthKpYdlLht8AFgJPp5T8rCpfXQdcBNTkuhCpERLwVESMi4hzcl1MUzGgS5JarIgoBe4DfpRSWprreqSGpJSqU0qDgJ2BAyPCW4iUdyJiBLAwpTQu17VIjXRYSmk/4EvAD7K3arZ4BvTmNwfoVe/9ztk2SdJWyN7Pex9wR0rp/lzXI21OSqkceA4YnutapAYcChyXva/3buDIiLg9tyVJG5dSmpP9uRB4gMytxS2eAb35vQ7sHhG7RkQr4FTg4RzXJEktWvbBW6OBd1JKv8t1PdLGRETXiCjLLpeQeWjsu7mtStpQSumSlNLOKaU+ZP69+veU0hk5LktqUES0yz4klohoBxwNbBffQmRAb2YppSrgh8CTZB5idE9KaVJuq5IaFhF3AS8D/SNidkR8N9c1SRtxKHAmmRmeN7KvL+e6KKkBPYDnImIimT/aP51S8uurJGnrdAdejIg3gdeAx1JKf8txTU3Cr1mTJEmSJCkPOIMuSZIkSVIeMKBLkiRJkpQHDOiSJEmSJOUBA7okSZIkSXnAgC5JkiRJUh4woEuStI1ERJd6Xws3PyLm1HvfqoHtiyKifCNj3R4RJ2zBvq+KiB9tTf0bGXd27fd8r9d+ZEQctIl+X4+IS7dwX89GRMdPU6ckSS1BUa4LkCTpsyKltBgYBBARVwDLU0rX5rSo5nMk8BHwykbW/zcwfAvHvBM4D7hmK+qSJClvOYMuSVIeiIhHImJcREyKiLPXW3d9tv3piOjSQN8DIuL5bP8nIqL7Zva1e0Q8md3+hYjYI9t+e0SMioiXIuKDiPhqtr0wIm6IiHcj4qmI+Nt6s/c/iogJETExIvaIiL7A2cB/Z68OOGS9/Q8AlqWUltTb7x8i4tWImBYRn4+IW7P7G12v60PAaY0/q5IktSwGdEmS8sO3Ukr7AwcAP46ITtn2jsC/Ukp7Ay8DP6/fKSJaA6OAE7P9bweu3My+bgTOz25/CfD7euu6AYcCJwC/yradBPQEBgBnAQevN96ClNJg4P+AH6eUpmWXf5NSGpRSemm97Q8Fxq3X1jGlNBS4CHiEzCz5AGD/iBgIkFL6CGjf0CX1kiRtD7zEXZKk/HBBRByXXd4Z6Au8AVQB92bbbydzmXd9ewF7A89EBEAhMHtjO8mG24OA+7Lbw7r/HngwpZSAiRHRM9t2GHBPSqkGmBsRz6837P3Zn+OAL2/mOAF6AIvWa3sk+/MtYG5KaXK23slAH+Dt7PpF2f4N3psvSVJLZkCXJCnHIuKLwOeBg1JKqyLiRaDNRjZP63cHJqaUDm/s7oCPUkqDNrJ+9XrbNkZtn2oa92+LVWx4fLVj1KxXQ816Y7bJ9pckabvjJe6SJOVeR+DjbDjfm8xl7rWKgK9ll08DXlyv72SgZ0QcCBARrbJjNCh73/e8eveXF0TEvpup71/A1yOjB5k/JmzOMqD9Rta9A/RrxBjriIhCYAdg5pb2lSSpJTCgS5KUe48BbbOXc18FvFpvXQVweERMInOp+VX1O6aUVgNfB34XEROBCcDQzezvVOC8iHgTmASM2Mz29wALyQTrW7L7qNhMn4eAk7MPjztkvXX/AIZspn9DDgBezF5qL0nSdicyt5lJkiRtXESUppSWR0RXMn9AGJpSWv8+8i0Z7w/AvSmlf2xhn3tSSuvfAy9J0nbBe9AlSVJjPBERHYBi4PKtCedZVwH7b2GfCYZzSdL2zBl0SZIkSZLygPegS5IkSZKUBwzokiRJkiTlAQO6JEmSJEl5wIAuSZIkSVIeMKBLkiRJkpQHDOiSJEmSJOUBA7okSZIkSXnAgC5JkiRJUh4woEuSJEmSlAcM6JIkSZIk5QEDuiRJjRQR0yPii7muo1ZEXBoR/5frOrZERFwREbfnuo7mEBEpIvrlug5JUstlQJck5b1cBOOIuCUirtqW+9yUiBgWEbPrt6WUfplSOruZ9pdXf4xYX0S0jojRETEjIpZFxBsR8aVc1yVJ0tYwoEuSpJaoCJgF/BvQEfgZcE9E9GmOnUVEYXOMmx27qLnGliS1LAZ0SVKLFhEjsrOn5RHxUkTsU2/d9Ij4r4iYGBEVEfGXiGhTb/1FETEvIuZGxNm1lyhHxDnA6cBFEbE8Ih6pt8tBDY0XETtExKPZOj6OiH9GRIP/PxsRe0bE09ntpkTEyfXWfTkiJmdnhedk628HPAHslK1neUTsVP9y8Yjok63/2xExKyKWRMR5EXFAtt7yiPh9vf30jYi/R8TiiPgoIu6IiLLsujFAb+CR7L4uyrYflD3H5RHxZkQMqzfeWRHxQbbuDyPi9E382tpkz92yiBgfEftmx/jviLhvvXN1fUSMWn+AlNKKlNIVKaXpKaWalNKjwIfA/tl+wyJidvZ3vDD7ez4he37fy577SzdWYPYKij9GxOMRsQI4IiL+ERFn19vmrIh4cSP9W0fEtRExMyIWRMQNEVGyXm0XR8R84OZNnCtJ0meIAV2S1GJFxGDgz8C5QBfgT8DDEdG63mYnA8OBXYF9gLOyfYcDPwa+CPQDhtV2SCndCNwB/DqlVJpS+srmxgMuBGYDXYHuwKVAaqDmdsDTwJ1AN+BU4P9FxIDsJqOBc1NK7YGBwN9TSiuALwFzs/WUppTmbuS0DAV2B04BrgN+mj3GvYGTI+LfaksBfgXsBOwF9AKuyB7/mcBM4CvZff06InoCjwFXAZ2B/wLui4iu2WO6HvhStu5DgDc2Uh/A8cC92XHuBB6MiGLgdmB4vT8UFGXPz22bGIvstt2BPYBJ9Zp3BNoAPYHLgJuAM8iE+MOBn0fErpsY9jTgF0B7oMEgvglXZ+sZRObzVVtD/do6A7sA52zh2JKk7ZQBXZLUkp0D/Cml9GpKqTqldCuwGjio3jbXp5TmppQ+Bh4hE5ggE7RvTilNSimtJBtOG2Fj460FegC7pJTWppT+mVLaIKADI4DpKaWbU0pVKaUJwH3ASfXGGRARHVJKS1JK4xtZV60rU0qVKaWngBXAXSmlhSmlOcA/gcEAKaX3U0pPp5RWp5QWAb8jc7n4xpwBPJ5Sejw7Y/00MBb4cnZ9DTAwIkpSSvNSSpM2OhKMSyn9NaW0NrvfNsBBKaV5wAv1zsVw4KOU0rhNHXA23N8B3JpSerfeqrXAL7L7uRvYARiVUlqWrW8ysO8mhn4opfSv7PFWbqqG9eoJMp/NC1JKH6eUlgG/JPPHhlo1wOXZ87+qsWNLkrZvBnRJUku2C3Bh9pLr8ogoJzMTvFO9bebXW14JlGaXdyJzD3Ot+subsrHxfgO8DzyVvdT7J5uoeeh6NZ9OZkYV4EQyoXdGRDwfEQc3sq5aC+otr2rgfSlkZpwj4u7sZfRLycxe77CJcXcBTlqv7sOAHtkZ/lOA84B5EfFYROy5ibHqznVKqYbMlQe1v7NbyfwxgOzPMZs62OxtBGOANcAP11u9OKVUXe/YYSPnY3N1bqGuQFtgXL1z9bdse61FWxL6JUmfDQZ0SVJLNovMDGlZvVfblNJdjeg7D9i53vte661vaPZ7o7KzshemlHYDjgN+HBFf2EjNz69Xc2lK6fvZcV5PKR1P5vL3B4F7Pk09jfDL7JifSyl1IBOGo/4hNVD3mPXqbpdSujpb95MppaPIXEXwLpnLyTem7lxnA/bOQO0l+w8C+0TEQDJXG9yxsUGyM9WjydxScGJ2prwprX8OVpAJ3rV2pGEfkQn/e9c7Vx1TSvX/GNDUv09J0nbAgC5JaimKI6JNvVcRmRB4XkQMjYx2EXFsRLRvxHj3AN+OiL0ioi3w8/XWLwB2a2xxkXlYXb9saKwAqslcxry+R4E9IuLMiCjOvg7I1tEqIk6PiI7ZsLm03hgLgC4R0bGxNW1Ge2A5UJG9v/y/11u//vHfDnwlIo6JiMLs72BYROycnY0/Pnsv+ursuA0de639I+Jr2d/hj7J9XgHIzir/lcy96a+llGZuYpw/krl//ivb6DLxN4CvRUTbyHzf+Xcb2ih7VcBNwMiI6AYQET0j4phtUKMkqQUzoEuSWorHycxK1r6uSCmNBb4H/B5YQuYS87MaM1hK6QkyDzZ7Ltvvleyq1dmfo8ncC14eEQ82YsjdgWfIhNOXgf+XUnqugf0uA44mcz/yXDKXzF8D1D7Y7kxgevay8/PIXP5O9t7qu4APsjXttP7YW+h/gP3I/DHhMeD+9db/CvhZdl//lVKaRebhbpcCi8jMqP83mX9LFJB54N5c4GMy97J/fxP7fojMJfFLssf7tfVmv28FPscmLm+PiF3IPBxwEDA/Pnm6/aaeHr+1RpK5lH5BtsaNzu4DF5P9XGV/l88A/ZuxNknSdiAafn6NJEmfLRGxF/A20DqlVJXrej7LIqI3mcvkd0wpLc11PZIkbSvOoEuSPrMi4qvZ76vuRGYW+xHDeW5l70n/MXC34VyS9FljQJckfZadCywEppG5Z3xTl2WrmWXvYV8KHAVcnuNyJEna5rzEXZIkSZKkPOAMuiRJkiRJeaAo1wVsjR122CH16dMn12VIkiRJkrRR48aN+yil1HVz27XogN6nTx/Gjh2b6zIkSZIkSdqoiJjRmO28xF2SJEmSpDxgQJckSZIkKQ8Y0CVJkiRJygMGdEmSJEmS8oABXZIkSZKkPGBAlyRJkiQpDxjQJUmSJEnKA80W0COiV0Q8FxGTI2JSRPxntr1zRDwdEVOzPztl2yMiro+I9yNiYkTs11y1SZIkSZKUb5pzBr0KuDClNAA4CPhBRAwAfgI8m1LaHXg2+x7gS8Du2dc5wB+bsTZJkiTl2LVPTuG7t7zO2uqaXJciSXmh2QJ6SmleSml8dnkZ8A7QEzgeuDW72a3ACdnl44HbUsYrQFlE9Giu+iRJkpRb/3hvIc++u5DKtdW5LkWS8sI2uQc9IvoAg4FXge4ppXnZVfOB7tnlnsCset1mZ9vWH+uciBgbEWMXLVrUbDVLkiRJkrQtNXtAj4hS4D7gRymlpfXXpZQSkLZkvJTSjSmlISmlIV27dm3CSiVJkiRJyp1mDegRUUwmnN+RUro/27yg9tL17M+F2fY5QK963XfOtkmSJEmStN1rzqe4BzAaeCel9Lt6qx4GvpVd/hbwUL32b2af5n4QUFHvUnhJkiRJkrZrRc049qHAmcBbEfFGtu1S4Grgnoj4LjADODm77nHgy8D7wErg281YmyRJkiRJeaXZAnpK6UUgNrL6Cw1sn4AfNFc9kiRJkiTls23yFHdJkiRJkrRpBnRJkiRJkvKAAV2SJEmSpDxgQJckSZIkKQ8Y0CVJkiRJygMGdEmSJEmS8oABXZIkSZKkPGBAlyRJkiQpDxjQJUmSJEnKAwZ0SZIkSZLygAFdkiRJkqQ8YECXJEmSJCkPGNAlSZIkScoDBnRJkiRJkvKAAV2SJEmSpDxgQJckSZIkKQ8Y0CVJkiRJygPNFtAj4s8RsTAi3q7X9peIeCP7mh4Rb2Tb+0TEqnrrbmiuuiRJkiRJykdFzTj2LcDvgdtqG1JKp9QuR8RvgYp6209LKQ1qxnokSZIkScpbzRbQU0ovRESfhtZFRAAnA0c21/4lSZIkSWpJcnUP+uHAgpTS1Hptu0bEhIh4PiIO31jHiDgnIsZGxNhFixY1f6WSJEmSJG0DuQro3wDuqvd+HtA7pTQY+DFwZ0R0aKhjSunGlNKQlNKQrl27boNSJUmSJElqfts8oEdEEfA14C+1bSml1SmlxdnlccA0YI9tXZskSZIkSbmSixn0LwLvppRm1zZERNeIKMwu7wbsDnyQg9okSZIkScqJ5vyatbuAl4H+ETE7Ir6bXXUq617eDvB5YGL2a9f+CpyXUvq4uWqTJEmSJCnfNOdT3L+xkfazGmi7D7ivuWqRJEmSJCnf5eohcZIkSZIkqR4DuiRJkiRJecCALkmSJElSHjCgS5IkSZKUBwzokiRJkiTlAQO6JEmSJEl5wIAuSZIkSVIeMKBLkiRJkpQHDOiSJEmSJOUBA7okSZIkSXnAgC5JkiRJUh4woEuSJEmSlAcM6JIkSZIk5QEDuiRJkiRJecCALkmSJElSHjCgS5IkSZKUBwzokiRJkiTlAQO6JEmSJEl5oNkCekT8OSIWRsTb9dquiIg5EfFG9vXleusuiYj3I2JKRBzTXHVJkiRJkpSPmnMG/RZgeAPtI1NKg7KvxwEiYgBwKrB3ts//i4jCZqxNkiRJkqS80mwBPaX0AvBxIzc/Hrg7pbQ6pfQh8D5wYHPVJkmSJElSvsnFPeg/jIiJ2UvgO2XbegKz6m0zO9u2gYg4JyLGRsTYRYsWNXetkiRJkiRtE9s6oP8R6AsMAuYBv93SAVJKN6aUhqSUhnTt2rWp65MkSZIkKSe2aUBPKS1IKVWnlGqAm/jkMvY5QK96m+6cbZMkSZIk6TNhmwb0iOhR7+1XgdonvD8MnBoRrSNiV2B34LVtWZskSZIkSblU1FwDR8RdwDBgh4iYDVwODIuIQUACpgPnAqSUJkXEPcBkoAr4QUqpurlqkyRJkiQp3zRbQE8pfaOB5tGb2P4XwC+aqx5JkiRJkvJZLp7iLkmSJEmS1mNAlyRJkiQpDxjQJUmSJEnKAwZ0SZIkSZLygAFdkiRJkqQ8YECXJEmSJCkPGNAlSZIkScoDBnRJkiRJkvKAAV2SJEmSpDxgQJckSZIkKQ8Y0CVJkiRJygMGdEmSJEmS8oABXZIkSZKkPGBAlyRJkiQpDxjQJUmSJEnKAwZ0SZIkSZLygAFdkiRJkqQ80GwBPSL+HBELI+Ltem2/iYh3I2JiRDwQEWXZ9j4RsSoi3si+bmiuuiRJkiRJykeNCugR8blPMfYtwPD12p4GBqaU9gHeAy6pt25aSmlQ9nXep9ifJEmSJEktVmNn0P9fRLwWEedHRMfGdEgpvQB8vF7bUymlquzbV4CdG1+qJEmSJEnbr0YF9JTS4cDpQC9gXETcGRFHbeW+vwM8Ue/9rhExISKej4jDt3JsSZIkSZJalKLGbphSmhoRPwPGAtcDgyMigEtTSvdvyU4j4qdAFXBHtmke0DultDgi9gcejIi9U0pLG+h7DnAOQO/evbdkt5IkSZIk5a3G3oO+T0SMBN4BjgS+klLaK7s8ckt2GBFnASOA01NKCSCltDqltDi7PA6YBuzRUP/0/7d371F61fW9x99fciEhJAFkCEm4BNJAS2gNYbiJ4K1RFArYog0qRY8KnEqPHM/SAvYUyqHL9hRRVJbHqCgWoYYDESigXMrlRLklXHKDYDAJ5EIuEDIJISGX7/nj2aHDMJN5BvLM3jPzfq31rNn7t/dv78/M2mvWfGf/9m9nTsnM5sxsbmpq6sqpJUmSJEmqrHrvoH8X+BG1u+WvbW/MzGXFXfW6RMRJwNeA92XmhlbtTcDLmbk1Ig4GxgG/r/e4kiRJkiT1dPUW6CcDr2XmVoCI2AUYlJkbMvNf2+sQETcA7wf2joglwCXUZm3fFbi7Njqeh4sZ208ELouIzcA24LzMfLm940qSJEmS1BvVW6DfA/wpsL5Y3w24C3hPRx0y88x2mn/cwb43ATfVmUWSJEmSpF6n3tesDcrM7cU5xfJujYkkSZIkSVLfU2+B/mpETNy+Usy0/toO9pckSZIkSV1Q7xD3C4AbI2IZEMC+wF82LJUkSZIkSX1MXQV6Zj4WEX8IHFo0zc/MzY2LJUmSJElS31LvHXSAo4AxRZ+JEUFm/qwhqSRJkiRJ6mPqKtAj4l+BscCTwNaiOQELdEmSJEmSdoJ676A3A4dlZjYyjCRJkiRJfVW9s7jPoTYxnCRJkiRJaoB676DvDcyLiEeBTdsbM/PUhqSSJEmSJKmPqbdAv7SRISRJkiRJ6uvqfc3aAxFxIDAuM++JiN2Afo2NJkmSJElS31HXM+gR8UXg/wI/KJpGA79sVChJkiRJkvqaeieJ+xJwPNACkJm/A/ZpVChJkiRJkvqaegv0TZn5+vaViOhP7T3okiRJkiRpJ6i3QH8gIi4GBkfEJOBG4LbGxZIkSZIkqW+pt0C/EFgFzAbOBe4A/q5RoSRJkiRJ6mvqncV9G/DD4iNJkiRJknayugr0iFhIO8+cZ+bBOz2RJEmSJEl9UF0FOtDcankQ8Algr846RcQ1wCnAysw8vGjbC/gFMAZYBHwyM9dERABXAR8DNgCfzczH68wnSZIkSVKPVtcz6Jn5UqvP0sz8NnByHV1/CpzUpu1C4N7MHAfcW6wDfBQYV3zOAb5fTzZJkiRJknqDeoe4T2y1ugu1O+qd9s3MByNiTJvm04D3F8vXAvcDf1u0/ywzE3g4IvaIiJGZubyejJIkSZIk9WT1DnH/ZqvlLRRD09/mOUe0KrpfBEYUy6OBF1rtt6Roe1OBHhHnULvDzgEHHPA2I0iSJEmSVC31zuL+gUacPDMzIt4y+VwnfaYAUwCam5u71FeSJEmSpKqqd4j7V3a0PTOv7MI5V2wfuh4RI4GVRftSYP9W++1XtEmSJEmS1OvVNUkctWfO/yu1IeejgfOAicDQ4tMVtwJnF8tnA7e0av+rqDkWWOvz55IkSZKkvqLeZ9D3AyZm5jqAiLgUuD0zP7OjThFxA7UJ4faOiCXAJcA/AVMj4vPAYv7zWfY7qL1ibQG116x9rkvfiSRJkiRJPVi9BfoI4PVW66/zn5O7dSgzz+xg04fa2TeBL9WZR5IkSZKkXqXeAv1nwKMRMa1YP53aK9IkSZIkSdJOUO8s7v8YEXcCJxRNn8vMJxoXS5IkSZKkvqXeSeIAdgNaMvMqYElEHNSgTJIkSerlbp+1nMWrN5QdQ5Iqpa4CPSIuAf4WuKhoGgBc16hQkiRJ6t3unLOc//LegxgysF/ZUSSpMuq9g/5x4FTgVYDMXEbXX68mSZIkvWHsPrsTEWXHkKTKqLdAf72YZT0BImJI4yJJkiRJktT31FugT42IHwB7RMQXgXuAHzYuliRJkiRJfUu9s7hfERGTgBbgUODvM/PuhiaTJEmSJKkP6bRAj4h+wD2Z+QHAolySJEmSpAbodIh7Zm4FtkXE8G7II0mSJElSn1TXEHdgPTA7Iu6mmMkdIDP/W0NSSZIkSZLUx9RboN9cfCRJkiRJUgPssECPiAMy8/nMvLa7AkmSJEmS1Bd19gz6L7cvRMRNDc4iSZIkSVKf1VmBHq2WD25kEEmSJEmS+rLOCvTsYFmSJEmSJO1EnU0S9+6IaKF2J31wsUyxnpk5rKHpJEmSJEnqI3ZYoGdmv+4KIkmSJElSX1bva9Z2mog4FPhFq6aDgb8H9gC+CKwq2i/OzDu6OZ4kSZIkSaXo9gI9M+cDEwAioh+wFJgGfA74VmZe0d2ZJEmSJEkqW2eTxDXah4DnMnNxyTkkSZJUkhtnLCk7giRVQtkF+mTghlbr50fErIi4JiL2bK9DRJwTETMiYsaqVava20WSJEk9xL+c8Sdcfd+CsmNIUiWUVqBHxEDgVODGoun7wFhqw9+XA99sr19mTsnM5sxsbmpq6paskiRJaoyjDtqr7AiSVBll3kH/KPB4Zq4AyMwVmbk1M7cBPwSOLjGbJEmSJEndqswC/UxaDW+PiJGttn0cmNPtiSRJkiRJKkm3z+IOEBFDgEnAua2a/3dETAASWNRmmyRJkiRJvVopBXpmvgq8q03bWWVkkSRJkiSpCsqexV2SJEmSJGGBLkmSJElSJVigS5IkSZJUARbokiRJkiRVgAW6JEmSJEkVYIEuSZIkSVIFWKBLkiRJklQBFuiSJEmSJFWABbokSZIkSRVggS5JkiRJUgVYoEuSJEmSVAEW6JIkSZIkVYAFuiRJkiRJFWCBLkmSJElSBVigS5IkSZJUARbokiRJkiRVgAW6JEmSJEkV0L+sE0fEImAdsBXYkpnNEbEX8AtgDLAI+GRmrikroyRJkiRJ3aXsO+gfyMwJmdlcrF8I3JuZ44B7i3VJkiRJknq9sgv0tk4Dri2WrwVOLzGLJEmSJEndpswCPYG7ImJmRJxTtI3IzOXF8ovAiHKiSZIkSZLUvUp7Bh14b2YujYh9gLsj4pnWGzMzIyLbdiqK+XMADjjggO5JKkmSJElSg5V2Bz0zlxZfVwLTgKOBFRG/fQPjAAAOiUlEQVQxEqD4urKdflMyszkzm5uamrozsiRJkiRJDVNKgR4RQyJi6PZl4MPAHOBW4Oxit7OBW8rIJ0mSJElSdytriPsIYFpEbM9wfWb+KiIeA6ZGxOeBxcAnS8onSZIkSVK3KqVAz8zfA+9up/0l4EPdn0iSJEmSpHJV7TVrkiRJkiT1SRbokiRJkiRVgAW6JEmSJEkVYIEuSZIkSVIFWKBLkiRJklQBFuiSJEmSJFWABbokSZIkSRVggS5JkiRJUgVYoEuSJEmSVAEW6JIkSZIkVYAFuiRJkiRJFWCBLkmSJElSBVigS5IkSZJUARbokiRJkiRVgAW6JEmSJEkVYIEuSZIkSVIFWKBLkiRJklQBFuiSJEmSJFVAtxfoEbF/RNwXEfMiYm5EfLlovzQilkbEk8XnY92dTZIkSZKksvQv4ZxbgP+RmY9HxFBgZkTcXWz7VmZeUUImSZIkSZJK1e0FemYuB5YXy+si4mlgdHfnkCRJkiSpSkp9Bj0ixgBHAI8UTedHxKyIuCYi9uygzzkRMSMiZqxataqbkkqSJEmS1FilFegRsTtwE3BBZrYA3wfGAhOo3WH/Znv9MnNKZjZnZnNTU1O35ZUkSZIkqZFKKdAjYgC14vznmXkzQGauyMytmbkN+CFwdBnZJEmSJEkqQxmzuAfwY+DpzLyyVfvIVrt9HJjT3dkkSZIkSSpLGbO4Hw+cBcyOiCeLtouBMyNiApDAIuDcErJJkiRJklSKMmZxnw5EO5vu6O4skiRJkiRVRamzuEuSJEmSpBoLdEmSJJVmYP9d2Lh5K/94+7yyo0hS6SzQJUmS1G3mLWvhlO/+P6YvWA3AsEED+N6nJ/L486+UnEySylfGJHGSJEnqo+Ytb2HO0pY3tQ3d1T9JJQm8gy5JkiRJUiVYoEuSJKkU7b3WR5L6MscTSZIkqVu9/9AmTv7jkXx4/Iiyo0hSpVigS5IkqVvtNWQgn2jev+wYklQ5DnGXJEmSJKkCLNAlSZIkSaoAC3RJkiRJkirAAl2SJEmSpAqwQJckSVK3eWXD62VHkKTKskCXJElSt3h6eQtX3DWfiQfsWXYUSaokC3RJkiR1i1c3bWH8qOF85tgDy44iSZVkgS5JkiRJUgVYoEuS1Iu88PIGjrjsLq6653dlR5He5NVNW/ibG55g32GDyo4iSZVlgS5JUi+yev0m1mzYzMLV68uOIr3h4mmzmXTlA6zfuIXvnHlEl/o+/vwajvvGvdw198UGpZOk6qhcgR4RJ0XE/IhYEBEXlp1Hkra74tfzmXTlA3zp+sfLjiKpj7hn3go+8q0HeXp5Cxs3b2XylIeYdOUD/PyRxe/ouCtaNvJn353O9Y88v5OSvtWcpWv5yLce5P75K5m15BWWrd3I8N0G0G+XeMu+uw3sz7xlLUy68gHun78SgI2bt/KXP3iIv77ucZav3cjXfzmH+S+uA+DFtY3P35mf/GYhp31vOqvXb3rHx7r83+fx6R89zOat23ZCMkk9Wf+yA7QWEf2Aq4FJwBLgsYi4NTPnlZtMUl837Ykl/OQ3C7n01PFcdtu8N7XfNHMp3548gb1337XEhJJ6qmdXrOPvps1hWyaDB/bjqslHsNeQgQA8teQV5q9Yx7Mr1rHP0F2Zt6yFc983lpmL1vDpY7o20drUGS9w21PLuGryESxZs4HZS9cybvHLfOqYA7p8nKmPvQDAYaOGcdlph79p+0U3z6bltc08t2o981es4+KbZ/PKa5u5/gvHcMi+Q9s95mGjhnHnl0/g2ocWcfHNsxm7z+5c8meH8cjClwE4bcIo1mzYzLMr1vHT3y5k1pK1zF3WwiGL13SYf+pjL3DbrGV8Z/IR7Fn8PDtz1T2/Y/6KFq6afAQD+u2yw+PMWLSGp5as5bM/eZRB/fsBcOC7hvDCmg3805//MQc37V7XOQF++9xLzFvewqYt2944b0/wwLOruPq+BfzzX/wJI4cP4m9ueIIjD9yT8943tuxoUo8VmVl2hjdExHHApZn5kWL9IoDM/EZ7+zc3N+eMGTO6MeHbd8ktc1iy5rWyY0h6m+Yua+Hs94zhrOMOZOJld3PCuL3faH+xZSNHHrgnewweUHJKCda+tpkZi9cwYtiuHD5qeNlxVIdlazcyavggznv/WL4+bTbDBw9g2KDa75PnVq1n0UsbGD9qGHsNGcj8F9fx9ZP/iG/c8QzjRw3r0nlmL13LynWbOGrMnmxLmLl4DfsOG/S2jvNXxx3I+FHDOfe6mZzwB3u/afu9z9TugI8fNYyffu5oFr30KoMH9OPw0Z1fjxs3b2X20rV8fdpsBg/sz1MvvALA5KP2Z+Pmrfxu5XrmLmt5Y/8d5W/9/W7/eXbm0YUvs27TFk48pIkBxZ3+jo4za+laVq2r3T2/7vPH8OjCl/jOfywA4PDRwxgxtP5n7R9Z+DLrN23hfYc00b+dEQZVtWDVeha/tIHDRw9jj8EDmb5gNcMHD6D5QF+jp+7x1ZMO5Q/37drvsLJExMzMbO50v4oV6GcAJ2XmF4r1s4BjMvP8VvucA5xTrB4KzO/2oG/P3sDqskNIdfJ6VU/htaqewmtVPYnXq3qKnnStHpiZTZ3tVKkh7vXIzCnAlLJzdFVEzKjnPyZSFXi9qqfwWlVP4bWqnsTrVT1Fb7xWq/aQy1Jg/1br+xVtkiRJkiT1alUr0B8DxkXEQRExEJgM3FpyJkmSJEmSGq5SQ9wzc0tEnA/8GugHXJOZc0uOtbP0uGH56tO8XtVTeK2qp/BaVU/i9aqeotddq5WaJE6SJEmSpL6qakPcJUmSJEnqkyzQJUmSJEmqAAv0bhARJ0XE/IhYEBEXlp1H6khEXBMRKyNiTtlZpB2JiP0j4r6ImBcRcyPiy2VnktoTEYMi4tGIeKq4Vv+h7EzSjkREv4h4IiL+vews0o5ExKKImB0RT0bEjLLz7Cw+g95gEdEPeBaYBCyhNlP9mZk5r9RgUjsi4kRgPfCzzDy87DxSRyJiJDAyMx+PiKHATOB0f7eqaiIigCGZuT4iBgDTgS9n5sMlR5PaFRFfAZqBYZl5Stl5pI5ExCKgOTNXl51lZ/IOeuMdDSzIzN9n5uvAvwGnlZxJaldmPgi8XHYOqTOZuTwzHy+W1wFPA6PLTSW9VdasL1YHFB/vjqiSImI/4GTgR2VnkfoqC/TGGw280Gp9Cf4RKUk7TUSMAY4AHik3idS+Ysjwk8BK4O7M9FpVVX0b+BqwrewgUh0SuCsiZkbEOWWH2Vks0CVJPVZE7A7cBFyQmS1l55Hak5lbM3MCsB9wdET4CJEqJyJOAVZm5syys0h1em9mTgQ+CnypeFSzx7NAb7ylwP6t1vcr2iRJ70DxPO9NwM8z8+ay80idycxXgPuAk8rOIrXjeODU4rnefwM+GBHXlRtJ6lhmLi2+rgSmUXu0uMezQG+8x4BxEXFQRAwEJgO3lpxJknq0YuKtHwNPZ+aVZeeROhIRTRGxR7E8mNqksc+Um0p6q8y8KDP3y8wx1P5e/Y/M/EzJsaR2RcSQYpJYImII8GGgV7yFyAK9wTJzC3A+8GtqkxhNzcy55aaS2hcRNwAPAYdGxJKI+HzZmaQOHA+cRe0Oz5PF52Nlh5LaMRK4LyJmUfun/d2Z6eurJOmdGQFMj4ingEeB2zPzVyVn2il8zZokSZIkSRXgHXRJkiRJkirAAl2SJEmSpAqwQJckSZIkqQIs0CVJkiRJqgALdEmSJEmSKsACXZKkbhIR72r1WrgXI2Jpq/WB7ezfPyJe6eBY10XE6V049+URccE7yd/BcZdsf893m/YPRsSxO+h3RkRc3MVz3RsRw99OTkmSeoL+ZQeQJKmvyMyXgAkAEXEpsD4zryg1VON8EFgNPNzB9q8CJ3XxmNcD5wH//A5ySZJUWd5BlySpAiLitoiYGRFzI+ILbbZ9p2i/OyLe1U7foyLigaL/nRExopNzjYuIXxf7PxgRhxTt10XEVRHx24j4fUR8vGjvFxH/JyKeiYi7IuJXbe7eXxART0TErIg4JCLGAl8AvlqMDnhPm/MfBqzLzDWtznt1RDwSEc9FxIkRcW1xvh+36noL8Kn6f6qSJPUsFuiSJFXD2Zl5JHAU8JWI2LNoHw78JjPHAw8B/7N1p4jYFbgK+Iui/3XA/+rkXFOAvy72vwj4Xqtt+wDHA6cD3yjaPgGMBg4DPgsc1+Z4KzLzCOBHwFcy87li+V8yc0Jm/rbN/scDM9u0Dc/MY4CvAbdRu0t+GHBkRBwOkJmrgaHtDamXJKk3cIi7JEnV8N8j4tRieT9gLPAksAW4sWi/jtow79b+CBgP3BMRAP2AJR2dpChujwVuKvaHN/898MvMTGBWRIwu2t4LTM3MbcCyiHigzWFvLr7OBD7WyfcJMBJY1abttuLrbGBZZs4r8s4DxgBziu2riv7tPpsvSVJPZoEuSVLJIuJPgROBYzPztYiYDgzqYPds2x2YlZkn1Hs6YHVmTuhg+6Y2+9Zje5+t1Pe3xWu89fvbfoxtbTJsa3PMQUV/SZJ6HYe4S5JUvuHAy0VxPp7aMPft+gN/Xix/Cpjepu88YHREHA0QEQOLY7SreO57eavny3eJiHd3ku83wBlRM5LaPxM6sw4Y2sG2p4E/qOMYbxIR/YC9gee72leSpJ7AAl2SpPLdDuxWDOe+HHik1ba1wAkRMZfaUPPLW3fMzE3AGcCVETELeAI4ppPzTQbOi4ingLnAKZ3sPxVYSa2w/mlxjrWd9LkF+GQxedx72my7H2jupH97jgKmF0PtJUnqdaL2mJkkSVLHImL3zFwfEU3U/oFwTGa2fY68K8e7GrgxM+/vYp+pmdn2GXhJknoFn0GXJEn1uDMihgEDgEveSXFeuBw4sot9nrA4lyT1Zt5BlyRJkiSpAnwGXZIkSZKkCrBAlyRJkiSpAizQJUmSJEmqAAt0SZIkSZIqwAJdkiRJkqQK+P9a8Bt8IOxsqQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create the x-axis for the plot of the fitted function\n",
    "xaxis = np.linspace(minL, maxL, 10000) \n",
    "\n",
    "# Compute the fitted function for x_fit\n",
    "yaxis = gauss_extended(xaxis, *minuit_L30cm.args) \n",
    "\n",
    "# Plot the fitted function\n",
    "ax_L30cm.plot(xaxis, yaxis, '-', label='Fit')\n",
    "\n",
    "# Compute the chi2, the number of non-empty bins, the NDOF and the chi2-probability\n",
    "L30cm_chi2, L30cm_entries = calculate_chi2(gauss_extended, L30cm_x, L30cm_y, L30cm_sy, *minuit_L30cm.args)\n",
    "L30cm_NDOF = L30cm_entries - len(minuit_L30cm.args)\n",
    "L30cm_chi2_prob =  stats.chi2.sf(L30cm_chi2, L30cm_NDOF) \n",
    "\n",
    "\n",
    "# make the results ready to be plotted\n",
    "d = {'Entries': len(L30cm), \n",
    "     'Mean': L30cm.mean(), \n",
    "     'Std': L30cm.std(ddof=1), \n",
    "     'Chi2': L30cm_chi2, \n",
    "     'ndf': L30cm_NDOF, \n",
    "     'Prob': L30cm_chi2_prob, \n",
    "    }\n",
    "for name in minuit_L30cm.parameters:\n",
    "    d[name] = [minuit_L30cm.values[name], minuit_L30cm.errors[name]]\n",
    "\n",
    "# add these results to the plot \n",
    "text = nice_string_output(d, extra_spacing=2, decimals=4)\n",
    "add_text_to_ax(0.02, 0.95, text, ax_L30cm, fontsize=12)\n",
    "\n",
    "ax_L30cm.legend()\n",
    "\n",
    "# show the actual fit and fit results\n",
    "fig_raw"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***\n",
    "\n",
    "Start by taking a close look at the data, first by inspecting the numbers in the data file (yes, open the damn thing, and look over the numbers by eye!), and then by\n",
    "considering the histograms produced by running the notebook. \n",
    "\n",
    "To begin with, only consider the 30cm ruler measurements, and disregard the estimated/guessed uncertainties. You can then expand from there, as guided below by questions.\n",
    "\n",
    "# Questions:\n",
    "\n",
    "1. Consider the mean and width. Is the result as you would expect it? And do you think that it is close to the best possible (i.e. most accurate and precise) estimate? NOTE: Make sure that you know the difference between accuracy and precision!!! See \"Common definition\" in: http://en.wikipedia.org/wiki/Accuracy_and_precision\n",
    "\n",
    "2. Do any of the measurements looks wrong/bad/suspecious? Do you see any repeated mistakes done for obvious reasons? Would you correct or exclude any of the measurements and how would you justify this? This problem requires that you discuss internally, and then each do what you think most justified/best. Apply this to the list of measurements, and perhaps produce a new list with your accepted measurements in (to save the original data). How many measurements did you throw away in the end?\n",
    "\n",
    "3. Fit your accepted length measurements with a Gaussian distribution, possibly in a (small?) range around what you believe is the true value. Should the binning be changed for these fits? And is the Gaussian distribution justified? Also, do you see any \"human\" effects? Did any of your class mates (or you?) not read to mm precision?\n",
    "\n",
    "4. Once you have selected the measurements you want to use, calculate the mean, RMS and uncertainty on the mean. How much did your result improve in precision?\n",
    "\n",
    "\n",
    "#### Now repeat the above for the 2m folding rule\n",
    " \n",
    "***\n",
    "\n",
    "\n",
    "5. How much better/worse is the single measurement uncertainty from the 30cm ruler case to the 2m folding rule?\n",
    "\n",
    "6. The \"Pull\" distribution is defined as the plot of $z_i = \\left(x_i - \\mu \\right)/\\sigma_i$ where $\\mu$ is the *sample* mean of $x$, and $x_i$ and $\\sigma_i$ are the single measurements and their corresponding uncertainties. If the measurements and uncertainties are good, then it should give a unit Gaussian. Is that the case? And thus, were the uncertainty estimates/guesses reasonable? If not, then the pull standard deviation is often used to remove overly precise measurements, and scale the errors on the remaining measurements to a reasonable level.\n",
    "\n",
    "7. Try to calculate the weighted mean. Did you get a good Chi2 probability, when doing so? Do you need to discard more dubious measurements? Did the result improve further in precision?\n",
    "\n",
    "8. Does the length of the table seems to be different when measured with a 30cm and a 2m ruler? Quantify this statement! I.e. what is the difference, and what is the uncertainty on that difference? Is this significantly away from 0?\n",
    "\n",
    "9. If you were asked for the best estimate of the length of the table, what would you do? (If posssible, read Bevington page 58 bottom!)\n",
    "\n",
    "\n",
    "### Not too advanced questions:\n",
    "10. Consider the 2018 data file with additional information in (Gender and measurement speed), and find out if gender and measurement speed have an impact on the quality of the measurements.\n",
    "\n",
    "### Advanced questions:\n",
    "11. Is there any correlation between the errors on the measurements and the distance value? I.e. do you see any effect of those measuring e.g. too long having a smaller/larger uncertainty? What would the effect of this be?"
   ]
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}