Week1/original/ChiSquareTest/AnalyticalLinearFit_original.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analytical linear fit:\n",
    "This Python macro illustrates how to perform a linear fit **analytically**, i.e. without any minimisation (e.g. by Minuit), and thus much faster. This can be done, since one can differentiate the Chi2, set this to zero and solve for both intercept and slope.\n",
    "\n",
    "## References:\n",
    "- Note on course webpage\n",
    "- Bevington: Chapter 6\n",
    "\n",
    "## Author(s), contact(s), and dates:\n",
    "- Author: Troels C. Petersen (NBI)\n",
    "- Email:  petersen@nbi.dk\n",
    "- Date:   8th of November 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np                                     # Matlab like syntax for linear algebra and functions\n",
    "import matplotlib.pyplot as plt                        # Plots and figures like you know them from Matlab\n",
    "import seaborn as sns                                  # Make the plots nicer to look at\n",
    "from iminuit import Minuit                             # The actual fitting tool, better than scipy's\n",
    "import sys                                             # Module to see files and folders in directories\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sys.path.append('../../../External_Functions')\n",
    "from ExternalFunctions import Chi2Regression\n",
    "from ExternalFunctions import nice_string_output, add_text_to_ax"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Program settings:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = np.random                         # Random generator\n",
    "r.seed(42)                            # Set a random seed (but a fixed one - more on that later.)\n",
    "save_plots = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nexp = 10       # Number of datasets fitted\n",
    "Npoints = 9\n",
    "alpha0 = 3.6\n",
    "alpha1 = 0.3\n",
    "sigmay = 1.0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating and fitting data:\n",
    "In the following, we generate data points following a simple line and fit these both **analytically** and **numerically**. The linear fit (and only this) can be done analytically, as discussed in Barlow chapter 6.2 and outlined here: http://www.nbi.dk/~petersen/Teaching/Stat2019/StraightLineFit.pdf.\n",
    "The numerical fit of the line (and any other function) is done iteratively by Minuit. The code is slightly shorter, but a lot slower, as the code will typically test many (50+) possible combination of fit parameters!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Ana. Fit: a0= 3.596+-0.726  a1=0.273+-0.129  p=0.8048  Num. Fit: a0= 3.596+-0.726  a1=0.273+-0.129  p=0.8048\n",
      "  Ana. Fit: a0= 2.384+-0.726  a1=0.523+-0.129  p=0.8135  Num. Fit: a0= 2.384+-0.726  a1=0.523+-0.129  p=0.8135\n",
      "  Ana. Fit: a0= 3.201+-0.726  a1=0.398+-0.129  p=0.2321  Num. Fit: a0= 3.201+-0.726  a1=0.398+-0.129  p=0.2321\n",
      "  Ana. Fit: a0= 3.303+-0.726  a1=0.435+-0.129  p=0.1219  Num. Fit: a0= 3.303+-0.726  a1=0.435+-0.129  p=0.1219\n",
      "  Ana. Fit: a0= 2.912+-0.726  a1=0.328+-0.129  p=0.8453  Num. Fit: a0= 2.912+-0.726  a1=0.328+-0.129  p=0.8453\n",
      "  Ana. Fit: a0= 4.054+-0.726  a1=0.212+-0.129  p=0.2378  Num. Fit: a0= 4.054+-0.726  a1=0.212+-0.129  p=0.2378\n",
      "  Ana. Fit: a0= 3.553+-0.726  a1=0.292+-0.129  p=0.6738  Num. Fit: a0= 3.553+-0.726  a1=0.292+-0.129  p=0.6738\n",
      "  Ana. Fit: a0= 3.917+-0.726  a1=0.268+-0.129  p=0.3798  Num. Fit: a0= 3.917+-0.726  a1=0.268+-0.129  p=0.3798\n",
      "  Ana. Fit: a0= 4.626+-0.726  a1=0.151+-0.129  p=0.4542  Num. Fit: a0= 4.626+-0.726  a1=0.151+-0.129  p=0.4542\n",
      "  Ana. Fit: a0= 2.802+-0.726  a1=0.559+-0.129  p=0.6745  Num. Fit: a0= 2.802+-0.726  a1=0.559+-0.129  p=0.6745\n"
     ]
    }
   ],
   "source": [
    "# Arrays for storing fit results:\n",
    "array_alpha0 = np.zeros(Nexp)\n",
    "array_alpha1 = np.zeros(Nexp)\n",
    "array_Chi2 = np.zeros(Nexp)\n",
    "array_Prob = np.zeros(Nexp)\n",
    "\n",
    "# Loop, repeating the data generation and fit:\n",
    "for iexp in range(Nexp) : \n",
    "\n",
    "    # Generating data by filling values into (x,y) and associated uncertainties (here chosen to be 0 for x):\n",
    "    x = np.arange(Npoints)+1\n",
    "    ex = np.zeros_like(x)\n",
    "    y = alpha0 + alpha1*x + r.normal(0, sigmay, Npoints) #+ alpha2*x**2\n",
    "    ey = sigmay*np.ones_like(x)\n",
    "\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Analytical fit:\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Calculate the relevant sums (see note):\n",
    "    s = len(x)\n",
    "    sx = np.sum(x)\n",
    "    sxx = np.sum(x**2)\n",
    "    sy = np.sum(y)\n",
    "    sxy = np.sum(x*y)\n",
    "    \n",
    "    # Use sums in calculations:\n",
    "    Delta = sxx * s - sx**2\n",
    "    alpha0_calc = (sy  * sxx - sxy * sx) / Delta\n",
    "    alpha1_calc = (sxy * s   - sy  * sx) / Delta\n",
    "    sigma_alpha0_calc = sigmay * np.sqrt(sxx / Delta)\n",
    "    sigma_alpha1_calc = sigmay * np.sqrt(s   / Delta)\n",
    "\n",
    "    # So now you have data points and a fit/theory. How to get the fit quality?\n",
    "    # The answer is to calculate the Chi2 and Ndof, and from these two values get their\n",
    "    # probability using the ChiSquare function (e.g. from scipy):\n",
    "    Chi2_calc = 0.0\n",
    "    for i in range( Npoints ) : \n",
    "        fit_value = alpha0_calc + alpha1_calc * x[i]\n",
    "        Chi2_calc += ((y[i] - fit_value) / ey[i])**2\n",
    "    \n",
    "    Nvar = 2                     # Number of variables (alpha0 and alpha1)\n",
    "    Ndof_calc = Npoints - Nvar   # Number of degrees of freedom = Number of data points - Number of variables\n",
    "    \n",
    "    # From Chi2 and Ndof, one can calculate the probability of obtaining this\n",
    "    # or something worse (i.e. higher Chi2). This is a good function to have!\n",
    "    # We will discuss in class, why/how this function works, and how it plays a central role in statistics.\n",
    "    from scipy import stats\n",
    "    Prob_calc =  stats.chi2.sf(Chi2_calc, Ndof_calc) # The chi2 probability given N degrees of freedom (Ndof)\n",
    "    \n",
    "    # Fill the arrays with fit results (to produce plots of these at the end):\n",
    "    array_alpha0[iexp] = alpha0_calc\n",
    "    array_alpha1[iexp] = alpha1_calc\n",
    "    array_Chi2[iexp] = Chi2_calc\n",
    "    array_Prob[iexp] = Prob_calc\n",
    "    \n",
    "    \n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Numerical fit:\n",
    "    # ------------------------------------------------------------------ #\n",
    "    # Define a fit function:\n",
    "    def fit_function(x, alpha0, alpha1):\n",
    "        return alpha0 + alpha1*x\n",
    "\n",
    "    # Now we define a ChiSquare to be minimised (using probfit), where we set various settings and starting parameters:\n",
    "    chi2_object = Chi2Regression(fit_function, x, y, ey) \n",
    "    minuit = Minuit(chi2_object, pedantic=False, alpha0=1, alpha1=1, print_level=0)  \n",
    "    minuit.migrad();  # perform the actual fit\n",
    "    minuit_output = [minuit.get_fmin(), minuit.get_param_states()] # save the output parameters in case needed\n",
    "    \n",
    "    alpha0_fit = minuit.values['alpha0']\n",
    "    alpha1_fit = minuit.values['alpha1']\n",
    "    sigma_alpha0_fit = minuit.errors['alpha0']\n",
    "    sigma_alpha1_fit = minuit.errors['alpha1']\n",
    "    \n",
    "    # In Minuit, you can just ask the fit function for it:\n",
    "    Chi2_fit = minuit.fval # the chi2 value\n",
    "    Prob_fit = stats.chi2.sf(Chi2_fit, Ndof_calc) # The chi2 probability given N degrees of freedom (Ndof, taken from above!)\n",
    "    \n",
    "    # Let us see what the fit gives for the first couple of datasets:\n",
    "    if (iexp < 10) :\n",
    "        print(f\"  Ana. Fit: a0={alpha0_calc:6.3f}+-{sigma_alpha0_calc:5.3f}  a1={alpha1_calc:5.3f}+-{sigma_alpha1_calc:5.3f}  p={Prob_calc:6.4f}\"+\n",
    "              f\"  Num. Fit: a0={alpha0_fit:6.3f}+-{sigma_alpha0_fit:5.3f}  a1={alpha1_fit:5.3f}+-{sigma_alpha1_fit:5.3f}  p={Prob_fit:6.4f}\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsgAAAGoCAYAAABbtxOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3XtclFX+B/DPMyh3bEBABLmZSiu4ig7eEWjd1WRLKy3N/aVlmaZWhreiYmyz3NRdatdScdV2NbOs3azXqluuBHgFFa95A4VEQETxEgrCnN8fp7kJyoADM8Dn/XrNK+acZ575Pk8Dfng4zzmKEAJERERERCSpbF0AEREREZE9YUAmIiIiIjLBgExEREREZIIBmYiIiIjIBAMyEREREZEJBmQiIiIiIhMMyEREREREJhiQiYiIiIhMMCATEREREZlo0xg79fb2FiEhIY2xayIiIiKiBtm3b99FIYRPXdtZFJAVRZkJ4DkAAsBhAM8IIW7eafuQkBBkZWVZWisRERERUaNTFCXPku3qHGKhKEoAgJcAaIQQEQAcAIy9t/KIiIiIiOyTpWOQ2wBwURSlDQBXAOcbryQiIiIiItupMyALIQoALAaQD6AQwBUhxH9v305RlMmKomQpipJVUlJi/UqJiIiIiJqAJUMsPAGMBBAKwB+Am6Iof7h9OyHECiGERgih8fGpc+wzEREREZFdsmSIxVAAZ4QQJUKIWwC+AjCwccsiIiIiIrINSwJyPoD+iqK4KoqiAPgNgB8btywiIiIiItuwZAzyHgAbAeyHnOJNBWBFI9dFRERERGQTFs2DLIRIApDUyLUQEREREdkcl5omIiIiIjLBgExEREREZIIBmYiIiIjIBAMyEREREZEJBmQiIiIiIhMMyEREREREJhiQiYiIiIhMMCATEREREZlgQCYiIiKiRqfVaqEoiuGh1WptXdIdKUIIq+9Uo9GIrKwsq++XiIiIiJo3RVHQGPnTwvfeJ4TQ1LVds7qC/MYbbxh+64iIiGjQPrRabYNfS0REREQtX7MKyPPmzUNhYSESEhJsWsetW7cwffp0eHl5oV27dnjmmWfw888/W/U9tFotHnjgAbi6uqJDhw4YP348CgsL67WPvXv3IiYmBu3atYOPjw+mT5+OiooKs22WLVuG4OBguLi4IC4uDqdPnzb0HTp0CI8++ig6duwINzc39OnTB1999ZVVju9u5s2bh+joaKvu88qVKxg/fjw8PDzg7e2NhIQEVFdXW/z62/8spH+Eh4cDAKqqqpCQkICwsDC4uroiODgYs2bNQnl5udl+jh07hhEjRsDDwwP33XcffvOb31j1OImIiOjeNauA7O7uDj8/P7i7u9u0jrfeegtffPEFNm7ciC1btiA1NRWvvPKKVd+jc+fO+Oijj3D06FFs2bIF58+fx5gxYyx+fVlZGeLj49G1a1fs27cPGzZswKZNmzBv3jzDNlu2bMH06dORmJiIvXv3wtXVFQ8//LAhOGZnZyMkJASff/45Dh8+jLFjx2LMmDH44YcfrHqspq5evYply5Zh7ty5Vt3vlClTsG/fPmzbtg1r167F6tWrsXDhQotfP2vWLBQWFpo9unXrhscffxwAUFFRgcOHD+Pdd9/FwYMHsWrVKnz++eeYMWOGYR8XLlww/MKyfft27Nu3D1OmTLHqcRIREZEVCCGs/ujTp49oiM2bN4vo6GihVquFq6urePDBB0V2dnaN7ZKSkkR4eHit+wAgkpKSRM+ePYWTk5MYPny4KC4urvHa5ORk4evrK3x8fERycrKhv6qqSkyaNEmEhoYKR0dHERgYKJKSkkR1dbWh38vLS3zwwQeG16xdu1Y4OTmJ69evN+i4LbFp0yahKIq4ceOGRdv/5z//EYqimNX04Ycfivvuu0/cunVLCCHEyJEjxaOPPmroP3funAAgvvvuuzvut0ePHuKVV16pd/1nzpwR8uN2dwsXLhQRERFCp9PV+z3u5OLFi0KlUomvv/7a0PbOO++ITp06NXife/bsEYqiiNzc3Dtus2TJEqFWqw3Pk5KSRNeuXQ2fJSIiotbIkjzQiO+dJSzIsnZ1BbmoqAgTJkxARkYG9u/fD39/f8THx+PWrVv12s/f/vY3LFiwAHv27EFxcbHZVTwAyM3NxalTp5Ceno6pU6ciISEBP/30EwCguroaDg4OWL16NY4fP46PP/4YycnJWL58ueG1ly5dMhsCMGTIEMMVxMZw6dIl/OMf/0B4eDicnZ0tek1lZSUURYGjo6OhzdnZGVeuXMGZM2cAAJmZmWbHERAQgM6dO+NON1gKIVBSUgJPT897OJo7q6ioQHJyMubMmQNFUay23/3790On09X4f3bu3DkUFRU1aJ8pKSmIjY1FaGjoHbcpLi42O1epqamIjY3FxIkT4evri969e+Pzzz9v0PsTERFR47GrgDxx4kRMmjQJ4eHhCAsLw9tvv42CgoJ6B89nn30W8fHx6NmzJ9555x18+eWXKCsrM/SrVCosWbIE3bp1w+uvvw6dTofs7GwAgKOjI5YvX46YmBiEhoYiPj4eI0eOxObNmwEAJSUlAABvb2+MGTMG/fr1g7e3t1mftXz77bdwd3dH+/bt8dNPP+G7776z+LV9+/aFk5MTFi1ahKqqKpw7d84Q8vV1lpSUwNvbGx9++CF8fX3x008/wdvb+47HsWLFCty4cQPPPffcvR9cLdasWQMnJyeMGzfOqvstKSmBg4MD1Go1oqKi8OSTT97T/7Pr16/js88+w6RJk+64zfnz57Fy5Uqz8fKFhYXYsGEDvL29sXXrVowdOxZjx45FZmZm/Q+KiIiIGo1dBeScnBw89dRT6Ny5Mzw8PNCjRw8AMpDUh/7GKf3X1dXVOHv2rKGtU6dOcHJyAgA4OTnB1dUVly5dMvQvW7YMGo0GPj4+cHd3x/r16w01iF+mJVEUBR07dkRQUFCDjtUScXFxOHDgAL7//ns4ODjUuBJ+Nx07dsSaNWuQnJwMZ2dndO/eHSNHjgQgf0HQUxQFXl5eCA4OhpOT0x2nXUlLS8OsWbPw6aefwt/f36Ia8vPz4e7uDnd3d8P/E/1zd3d35OfnG7bV6XRYvHgxEhIS0KZNmxr7mjJlitlr60MIYbipLjAwEP7+/vc0vcz69evh4OCAxx57rNb+69evY9SoURg+fDimTZtmaNfpdPDx8cGSJUsQGRmJOXPmICoqCuvWrWtwLURERGR9NZOIDT388MPo0KEDVq1ahYCAAFy4cAGDBw+GTqdr8D71Qcg0ENUWwPT9GzZswMsvv4wlS5YgJiYGLi4uSExMRHFxMQDA19cXgLzy+OGHHwKAYXiGj49Pg+usjZubG7p27YquXbuiW7duCAoKwt69e9G3b1+LXv/EE0/giSeeQGFhITw9PZGeng4AhoDr4+ODkpISzJw5E3/4wx8AAKWlpTWOY+/evXjkkUewbNkyjBgxwuL6/f39DVfmCwoKEBsba3huWgcAbNy4EWVlZXe8Kvv2229j1qxZFr+3KV9fX1RVVeHy5cuGWTj056Ih/89SUlIwbtw4uLi41Oi7ceMGHn74Yfj6+mLNmjVmfd7e3lCr1WbDRzp37oyCgoJ610BERESNx24CcmlpKX788Ud89NFHiI2NBQCcOHGi1m3d3d1x48aNO+7r6NGjhq+PHTsGlUqFkJAQi+rIyMjAgAEDMH36dENbbm4u3NzcAMhAow+bkZGRAOTVVScnJ8MV78bg4OAAADWmk7t58yaKioqgVquhVqtrfW3Hjh0ByCufoaGhhqveGo0G6enpmDlzJgAZYnNzc6HRGOfPzs7OxkMPPYTFixdj/Pjx9aq5TZs26NKli+FrAIbnt/vTn/6EGTNmwNXVtdZ+X19fwy8n9dW7d2+oVCqkp6fjkUceASD/nwUEBMDPz89s27KyMpSVlcHPz6/W8d4HDx5EZmYmli5dWqOvoqICjz76KNq0aYONGzeibdu2Zv2RkZHYtm2bWVt+fj6ioqIadFxERETUSCy5k6++j4bMYlFdXS28vb3F008/LU6fPi2+//570bNnTwFAbN++3Wzb//3vf0JRFPHNN9+IwsJC8fPPP5venSjat28v/vOf/4js7GwRGRkpRo8ebeivbQYMNzc3sXr1aiGEnOmhXbt24rvvvhMnTpwQc+fOFe7u7iImJsaw/dy5c4Wvr6/Ytm2b2LFjhwgJCRHPPfdcvY/5Tq5duyZeeuklsX37dnHmzBmxe/du8dvf/lZ06tRJXLt2zWzb7du3G2buuN0XX3wh9uzZI86cOSOWLFkiVCqV+OSTTwz9mzdvFg4ODmLFihXi0KFDYsSIEeKBBx4QVVVVQgghDh8+LNq3by/efPNNUVhYaHiUlpbW+5juNovFf//7X+Hm5tag/VrqySefFGFhYWLPnj1iy5YtwtPTU7zzzjs1tktKSqr1M6c3bdo0ERERUaO9srJS/P73vxe9evUSubm5ZudLfz4zMzOFSqUSb731ljh16pT4+OOPhaIoYs+ePVY9ViIiInt2pzzQRO9t0SwWdhOQhZDBNyIiQjg5OYlevXqJTZs23TGszJ49W3h6egoA4i9/+YvpgYs333xThIeHCycnJzFs2DBRVFRk6K8rIFdUVIhJkyYJtVot1Gq1mDFjhnjhhRfMAnJFRYWYNm2aUKvVwsPDQ0yYMMGqU7zduHFDPPbYYyIgIEA4OjoKX19fMWrUKHHs2LEa294tIC9atEj4+vqKtm3biu7duxuO0dRHH30kAgMDhZOTk4iJiREnT5409OnD4u0P03NhDQ8++GCDpo6rj7KyMjFu3Djh5uYmvLy8xKuvvmoIrqbuFpDLy8uFWq0Wf/7zn2v06X8BqO1x5swZw3YbNmwQYWFhwsnJSXTv3l2sX7/emodJRERk95pDQFZEI6yFrdFoxJ2mCmtsiqLgiy++wOjRo23y/lQ/WVlZGDhwIHJychAYGGjrcoiIiKiRKYpyTzfL3+N77xNCaOrazq5msaDWp7S0FCtXrmQ4JiIiIrthNzfpUes0bNgwW5dAREREZKbFBWRbXbInIiIiopaBQyyIiIiIiEw0q4D8xhtvGFZEi4iIaNA+tFptg19LRERERC1fswrI8+bNQ2FhIRISEmxah37BiQ4dOkBRFKSmpjbK+8yfPx9+fn5wc3PDqFGjcOHCBYtfe/nyZUyZMgWhoaFwcXFBly5dsGDBghqrEu7cuRPR0dFwdXVF+/btMXHiREOfVqs1/EJi+jBdyrsxzJs3D9HR0Vbd55UrVzB+/Hh4eHjA29sbCQkJqK6urtc+3nvvPQwYMADOzs41Fp6pqqpCQkICwsLC4OrqiuDgYMyaNQvl5eWGbSorK/H6668jODgYzs7O6NOnD9LS0qxxeERERGRFzSogu7u7w8/PD+7u7jat4+rVq4iIiMD777/faO+xfPlyvP/++/joo4+Qnp6OgoKCeq1kV1xcjNLSUixduhRHjhzBkiVLsGjRIixYsMCwzdGjRzF06FBERkZi9+7d2LlzJx599FFD/6xZs1BYWGj26NatGx5//HGrHqupq1evYtmyZZg7d65V9ztlyhTs27cP27Ztw9q1a7F69WosXLiwXvsoLy/HE088gaeffrpGX0VFBQ4fPox3330XBw8exKpVq/D5559jxowZhm3ef/99pKSkYPny5Thy5AhiYmIwYsQInD9//p6Pj4iIiKzIksmS6/to6EIhmzdvFtHR0UKtVgtXV1fx4IMPiuzs7Brb1bbYh8kE0CIpKUn07NlTODk5ieHDh4vi4uIar01OTha+vr7Cx8dHJCcnG/qrqqrEpEmTRGhoqHB0dBSBgYEiKSlJVFdX13ivkpKSu666di969uwpZs6caXiekZEhAIhTp041eJ8zZswQvXr1MjyfMGGCePDBBy1+/Z49e4SiKCI3N7fe7323lfRMLVy4UERERAidTlfv97iTixcvCpVKJb7++mtD2zvvvCM6derUoP0tWrRIBAcH17ndkiVLhFqtNjzv27evmD17tuG5fvXIRYsWNagOIiKi5siSPNCI723RQiF2dQW5qKgIEyZMQEZGBvbv3w9/f3/Ex8fj1q1b9drP3/72NyxYsAB79uxBcXGx2VU8AMjNzcWpU6eQnp6OqVOnIiEhAT/99BMAoLq6Gg4ODli9ejWOHz+Ojz/+GMnJyVi+fLnVjrMu+quRpsMMBgwYgLZt2+JeFmApLi6Gp6en4Xlqair69++Phx9+GL6+vhg4cCC2b99+x9enpKQgNjYWoaGhDa7hbioqKpCcnIw5c+ZAURSr7Xf//v3Q6XRm53PIkCE4d+4cioqKrPY+t7v9fFdWVsLJycnwXKVSwdHREQcOHGi0GoiIiKj+7CogT5w4EZMmTUJ4eDjCwsLw9ttvo6CgAIcPH67Xfp599lnEx8ejZ8+eeOedd/Dll1+irKzM0K9SqbBkyRJ069YNr7/+OnQ6HbKzswEAjo6OWL58OWJiYhAaGor4+HiMHDkSmzdvtuqx3k1paSl0Oh28vb0xa9YsdO7cGTdv3oSnpydKSkoatM/Dhw/j3//+t9n47cLCQixduhSDBw/G1q1b0b9/f4wYMcLwy4Kp69ev47PPPsOkSZMafFx1WbNmDZycnDBu3Dir7rekpAQODg5Qq9WIiorCk08+CW9vb0NfYzh//jxWrlxpdr5jYmLw6aefIjc3F9XV1fjoo49w4cKFRquBiIiIGsauAnJOTg6eeuopdO7cGR4eHujRowcAGc7qw/QmsvDwcFRXV+Ps2bOGtk6dOhmu5Dk5OcHV1RWXLl0y9C9btgwajQY+Pj5wd3fH+vXr613DvRC/zOWsKAp8fHwQFBSENm3aNHiO56KiIowaNQoJCQmIj483tOt0OvTp0wdz585FZGQkFi9eDE9PT3z55Zc19rF+/Xo4ODjgscces/h98/Pz4e7uDnd3d8P/E/1zd3d35Ofnm9WyePFiJCQkoE2bmtNzT5kyxey19SGEMNxgGBgYCH9//0adL/v69esYNWoUhg8fjmnTphnatVotwsPD0aVLFzg5OeHLL7/E8OHDoVLZ1bchERFRq2dXC4U8/PDD6NChA1atWoWAgABcuHABgwcPrjHzQn3og5BpIKotgOn7N2zYgJdffhlLlixBTEwMXFxckJiYiOLi4gbXUF/e3t5QqVQoKSnB3LlzMXfuXOh0OpSVlcHHx6de+7p48SKGDh2KoUOH4t13363xPt26dTM8V6lUCA4ORkFBQY39pKSkYNy4cXBxcbH4vf39/Q1X5gsKChAbG2t4ru/X27hxI8rKyu54hfrtt9/GrFmzLH5vU76+vqiqqsLly5fx1VdfAZAzkQCo9/msy40bNwxDVtasWWPWp1arsWnTJpSXl+Pq1avw8/PDoEGDEBYWZtUaiIiI6N7YTUAuLS3Fjz/+iI8++gixsbEAgBMnTtS6rbu7O27cuHHHfR09etTw9bFjx6BSqWpMy3UnGRkZGDBgAKZPn25oy83NhZubm0WvtwYnJydEREQgPT3dMKvErl27cOvWLWg0GrNtb968iaKiIqjVaqjVarO+y5cv43e/+x369OmDZcuW1XifyMhI5ObmGp4LIXDu3DkEBASYbXfw4EFkZmZi6dKl9TqONm3aoEuXLoavARie3+5Pf/oTZsyYAVdX11r7fX194evrW6/31+vduzdUKpVhej4ASEtLQ0BAAPz8/My2LSsrQ1lZGfz8/ODs7Fyv96moqMCjjz6KNm3aYOPGjWjbtm2t27m6usLV1RV5eXnYvXs3XnjhhQYdFxERETUOu/nbrqenJ7y9vbF69Wrk5ORg27ZteOONN2rdtk+fPjhz5gy+/fZbFBUVmc01CwCrVq3C5s2bcfDgQbzxxht47LHHzG6Wuptu3brhwIED+P7773Hy5EnMmzcPx48fN9vm+vXryM7OxpEjRwAAp0+fRnZ2ttkwjXs1depULF++HF999RX279+PV155BUOHDq0RMHfv3o3Q0FAkJyebtV+9ehXDhg2Dj48P3nvvPRQXF6OoqMjsprTnn38e3333HT7++GOcPn0aWq0WFy9erDGMIiUlBREREYiKirLa8Zn67rvvcOLECbNfSqypffv2GDNmDObMmYO9e/di69atWLJkCaZOnVpj2+TkZISGhmL37t01+vLz85GdnY3z58+jsrIS2dnZhivit27dwujRo1FcXIwVK1agrKzMcL718y1fuHAB//jHP5CTk4Pdu3djzJgx6NatG5588slGOW4iIiJqIEumuqjvo6HTvP3vf/8TERERwsnJSfTq1Uts2rTpjtOozZ49W3h6egoA4i9/+Yvp9B3izTffFOHh4cLJyUkMGzZMFBUVGfprmyLOzc1NrF69WgghREVFhZg0aZJQq9VCrVaLGTNmiBdeeEHExMQYtt++fbsAUOOh34e1JCUlCV9fX+Hi4iJGjhxpdhy315KUlFRre20PUx9++KEIDg4WLi4uonfv3uK7774z6y8vLxdqtVr8+c9/tuqxmXrwwQfFK6+80mj7F0KIsrIyMW7cOOHm5ia8vLzEq6++Kqqqqmpsl5SUdMfP3IQJE+54PvXT2NX2OHPmjBBCiOLiYtG7d2/h7Ows1Gq1GDt2rCgoKGjMwyYiIrI7t2eRJn5vi6Z5U0Qj3Kyk0WjEvUxHdi8URcEXX3yB0aNH2+T9qX6ysrIwcOBA5OTkIDAw0NblEBERUSNTFKVRb5av4733CSE0dW1nN0MsqHUqLS3FypUrGY6JiIjIbtjNTXrUOg0bNszWJRARERGZaXEB2VaX7ImIiIioZeAQCyIiIiIiEwzIREREREQmmlVATk1NhaIouHjxYp3barVaRERENEFVRERERNSS2F1APn36NB577DGo1Wp4eHhg0KBB2L9/f733M2vWLPzwww812lNSUhAVFYV27drBx8cHY8aMwdmzZ61QOTUH+tUILf1FS0+r1UJRlBqP8PDwWrcfOXIkFEWB6XSHN2/exLPPPovw8HA4ODhg4sSJ93o4RHSb279XtVqtrUsiombIrgJyYWEhBg4ciMuXL+Pbb7/F3r17MXbsWBQWFtZ7X+7u7mjfvn2N9oyMDDz33HPYsWMHvv/+e1y6dAnDhw9HVVWVNQ6B7Nz8+fMbtGz4rFmzUFhYaPbo1q0bHn/88Rrbrly5EmVlZTXaq6ur0aZNG8ycORN9+/ZtUP1EdHdardZws7YQggGZiBrGktVE6vto6Ep6M2bMEO3atRPXr1+vtV+/OtyGDRtEWFiY8PDwEBMmTBC3bt0ybPPee+8ZVjC7fcW82uzbt08AENnZ2Q2qmZqPHTt2iB49eoitW7cKAKKkpKTB+9qzZ49QFEXk5uaatZ8+fVrcf//9hs9VZmZmra+Pj48XEyZMaPD7E9HdwYYrdRHR3dny+xMWrqRnV1eQt2zZgoceeqjOK3wrV67Exo0bsXbtWqxduxYbN2409M2YMQOFhYVISEiw6D2Li4sBAGq1uuGFk927fv06Jk6ciJSUFDg6Ot7z/lJSUhAbG4vQ0FBDW3V1Nf7v//4PCxcuhJeX1z2/BxEREdmGXQXk/Px8dOrUqc7t3nrrLUREROCRRx5B3759kZmZaehzc3ODn58f3N3d69xPdXU13n33XYwePRrBwcH3VDvZt5dffhmPPPII+vXrd8/7un79Oj777DNMmjTJrP3dd9+Fv78/lzknIiJq5uxuoRBFUercpkuXLoavvby8cOnSpQa914wZM1BaWopvvvmmQa+n5mHTpk1IS0vDoUOHrLK/9evXw8HBAY899pihbf/+/Vi6dCkOHjxolfcgIiIi27GrK8hBQUE4d+5cndu1aWOe60UDVs+bO3cutmzZgu+++47DK1q4//3vfzhz5gw8PT3h7OyM3/3udwCATp06ITk5ud77S0lJwbhx4+Di4mJoS0tLQ0lJCYKDg+Hs7IywsDAAwKBBgzBz5kzrHAgRERE1Cbu6gjxs2DD885//RHl5OVxdXRvtfbRaLdavX4+0tDQEBAQ02vuQfXj99dcxZcoUw/O9e/diwoQJSE1NRbdu3cy2LSsrQ1lZGfz8/ODs7FxjXwcPHkRmZiaWLl1q1j5hwgQMHz7c8LygoABDhw7Fp59+isGDB1v5iIiIiKgx2VVAnjdvHj777DM8/PDD+OMf/wgvLy9s27YNISEhiI+Pt2gfRUVFAOQ40aqqKsNzLy8vODo6YuHChfjggw/wzTffwNnZ2dB/3333mV0RpJbD19cXvr6+huf6/+ddunSpcTNdcnIy5s+fj+3btyM2NrbGvlJSUhAREYGoqCizdk9PT3h6ehqe68N1cHAwOnToYGg/duwYKisrcfXqVahUKmRnZ8Pd3d1s2BARERHZll0F5ICAAGRkZGDu3Ll46KGHUFVVhR49euCjjz6yeB8dO3as9bk+8CxbtgxlZWWIjo4222716tVcuIHu6saNG1i3bh3eeuutBu9jxIgRyMvLMzz/5ptvEBMTg9TUVCtUSERERNagNGT8bl00Go0wXUGMiIioKSmK0qD7U4io8dny+1NRlH1CCE1d29nVTXpERERERLbGgExEREREZIIBmYiIiIjIBAMyEREREZGJZh+QY2NjMX369BrtBw4cQK9evdC2bVuLVucjIiIiIgJaQEC+k3nz5iEwMBA5OTkoLCy0dTlEiI2NhaIoNR7Tpk2zdWlERERkwq7mQbamnJwcvPDCCwgKCrJ1KUQAgK+++gqVlZWG5xcuXEBkZCQef/xxG1ZFREREt7OrK8gTJ07E73//e8ybNw9qtRqBgYH44osvzLb55JNPEBgYCDc3N0ydOhU6nc6sPyQkBIqiICcnB3PmzDFcpSOyNS8vL/j5+RkeW7ZsQXBwMOLi4mxdGhEREZmwq4AMyBXvPDw8kJmZieHDh2Py5Mm4ceMGALlM77PPPouEhATs378fALBz506z12dmZqKwsBCdOnXCW2+9hcLCQg6xILu0cuVKTJw4kb/AERER2Zk6A7KiKGGKomSbPK4qivJKYxXk7++PxMREdO3aFbNnz0ZZWRlOnz4NAPj73/8OjUaDV155BWFhYfjLX/4CNzc3s9f7+PjAz88PDg73AlduAAAgAElEQVQO8PDwMFytI7InqampyMnJ4fLmREREdqjOgCyEOCGE6CWE6AWgD4ByAP9qrILuv/9+w9deXl4AgEuXLgEATp8+je7duxv6nZ2d0bVr18YqhajRpKSkYOjQoRwjT0Q2p9VqzW4c1mq1ti6JyObqe5PebwDkCCHyGqMYAGjTpmZJput188/R1NxdunQJX375JT755BNbl0JEBK1WawjJpv/eErVm9R2DPBbA+to6FEWZrChKlqIoWSUlJfdeWS26dOmCo0ePGp5XVFQYhl8QNRf//Oc/4ebmhlGjRtm6FCIiIqqFxQFZURRHAI8A+KK2fiHECiGERgih8fHxsVZ9Zp599llkZmbigw8+wIkTJ/Dqq6/i+vXrjfJeRI0lJSUFTz31FJycnGxdChEREdWiPleQHwKwXwhR3FjF1CU8PByrVq3C4sWL0bt3b+h0OgwcONBW5RDV265du3D06FE8++yzti6FiIiI7kCxdLyRoiifAdgqhFhd17YajUZkZWXda21EREQNwvG09cdzRk3Flp81RVH2CSE0dW1n0RVkRVFcAfwWwFf3WhgRERERkT2zaBYLIUQ5gPaNXAsRERERkc3Z3Up6RERERES2xIBMRERERGSi1QVkrVaLiIgIW5dBRERERHbKrgLyxIkTDUtdOjk5ISwsDO+99x6qq6ttXRo1c1euXMH48ePh4eEBb29vJCQk1PtzVVZWhsmTJ8PX1xeurq7o1asXTp48aeg/efIk4uPj4enpCbVajT/84Q+4fPmyof9f//oXhgwZAk9PT3h6emL48OE4dOiQ1Y6RiIiIrMOuAjIADB06FIWFhThx4gSmTZuGxMRELF682NZlUTM3ZcoU7Nu3D9u2bcPatWuxevVqLFy40OLXCyEwatQo7NmzB+vXr8fRo0exaNEiuLm5AQCqqqrwyCOPQKVSYceOHdi6dSsOHz5sNt9xRkYGRo4ciW3btmHXrl3w9PTEb3/7W5SWllr9eImIiOgeCCGs/ujTp49oiAkTJoj4+HiztqFDh4r+/fsbnp85c0YAEN9884146KGHhIuLi+jQoYP4/vvvhRBC5OXlifj4eOHq6irUarWYOHGiuHbtmuH1SUlJIjw8XMyePVt4eHgIf39/kZKS0qB6qXm4ePGiUKlU4uuvvza0vfPOO6JTp04W72P79u2iTZs2Ijc3t9b+Y8eOCQDi6NGjhrZNmzYJRVFEcXFxra8pLS0VAMS///1vi+sgIsvIf96oPnjOqNGtXStEcLCoBoQIDpbPmxiALGFBlrW7K8i3c3Fxwa1bt2q0z5o1C6NGjcKhQ4ewZs0aw5W8p556CmVlZdi5cyc2bdqEtLQ0zJ492+y1x48fR2FhIfbu3Yu5c+di8uTJOHLkSJMcDzW9/fv3Q6fTITo62tA2ZMgQnDt3DkVFRRbtIzU1FeHh4diwYQM6deqEsLAwJCUloaqqCgBQWVkJAGbLRzs7O0MIgYMHD9a6z+JiuSilp6dng46LiIio2Vi3Dpg8GcjLk+EzL08+X7fO1pXVym4Dsk6nw3/+8x9s2bIFQ4cOrdH/xBNPYPLkyejSpQuGDx+O/v3748iRI9ixYweWLl2Knj17Ijo6GvPnz8eqVasMAQYAVCoV/vrXv+KBBx7ASy+9BI1Gg1WrVjXl4VETKikpgYODA9RqNaKiovDkk0/C29vb0GeJwsJCnDlzBt9//z2++uorvPvuu/jwww+RnJwMAAgLC0OHDh2wePFi3Lx5E5cuXcKSJUugKMod3+Ptt99G3759MXjwYOscKBERkb1KTATKy83bystlux2yu4C8ZcsWuLu7w9nZGY8++igmTJgArVZbYzvTq4F6p0+fhkqlMpulomfPnqisrER+fr6hLSAgAGq12vA8PDwcOTk51j0QshtCCMPNn4GBgfD396/3Epc6nQ5Xr17F6tWr0bdvXzz++OOYPHky1q5dC0BeLV6/fj22bNkCNzc3BAQEIDY2Vv6ZRlXz2+z9999HamoqPv/881r7iYiIWhSTHGZRu41ZtJJeUxoyZAhWrFgBJycn+Pv7w8HBodbtavuztD4I1eZO7XW9jpo/X19fVFVV4fLly/jqK7laenp6OgDAx8fHon14e3vD2dkZgYGBhrbOnTujoKDA8DwuLg5nzpzBhQsX4OrqisLCQrz22mvw9/c329fHH39sCMjBwcH3enhERET2LyhIDquord0O2d2lK1dXV3Tp0gWBgYF3DMd30qVLF1RXV5uNJz548CDatm1rFmwKCgpQVlZmeH7s2DF06dLl3osnu9S7d2+oVCpDKAaAtLQ0BAQEwM/Pz2zbsrIynD17Fjdv3jRrj4yMxM2bN3H+/HlDW35+PgICAmq8n6+vL9zd3bF+/Xq4u7tDo9EY+lavXo3ExERs3bqV83ETEVHrsWAB4Opq3ubqKtvtkN0F5HvRo0cPDBgwANOnT8fBgweRnp6OpKQkPPPMM3B0dDRsp9Pp8PLLL+PEiRP461//iqysLLPpuKhlad++PcaMGYM5c+Zg79692Lp1K5YsWYKpU6fW2DY5ORmhoaHYvXu3WXt8fDw6duyIKVOm4Mcff8R///tfrFixAuPHjzdss2XLFqSlpSEvLw9r1qzBn/70J8yePRuuv/xAWLduHaZNm4bVq1cjICAARUVFKCoqwvXr1xv3BBAREdna+PHAihVAcDB0ABAcLJ+b/DtqT+xuiMW9+vTTT/Hiiy9iwIABcHR0xMiRI7Fo0SKzbR544AF06NABGo0G7dq1w9///nd0797dRhVTU1i+fDmmTp2KBx98EE5OTnjmmWcwb948i1/v5uaGzZs3Y/r06YiMjISvry9eeOEFzJw507BNaWkpnn/+eRQXFyMwMBBarRazZs0y9KekpODGjRsYNWqU2b6TkpJqHWdPRETUoowfD4wfDwdFgTh71tbV3JVS35uVLKHRaERWVpbV90tERGQJRVHqfTNua8dzRk3Flp81RVH2CSE0dW3XooZYEBERERHdKwZkIiIiIiITDMhERERERCYYkImIiIiITDAgt2JardawwpyiKJxJgYiIiAicxYLAO5eJqOXhz7X64zmjpsJZLIiIiIiImhkGZCIiIiIiEwzIREREREQmGJCJiIiIiEwwIBMRERERmWBAJiIiIiIywYBMRERERGSCAZmIiIiIyAQDMhERERGRCQZkIiIiIiITDMhERERERCYYkImIiIiITDAgExERERGZYEAmIiIiqietVgtFUQwPrVZr65LIihQhhNV3qtFoRFZWltX3S41DURQ0xueAiMhW+HOt/njOGobnrf5sec4URdknhNDUtR2vIBMRERERmWBAJiIiIiIywYBMRERERGSCAZmIiIiIyAQDMhERERGRCQZkIiIiIiITDMhERERERCYYkImIqOVYtw4ICUE1AISEyOdERPXUxtYFEBERWcW6dcDkyUB5ubz6k5cnnwPA+PG2rIyImhleQSYiopYhMREoLzdvKy+X7URke+fPA+vWYbqt67AAryATEVHLkJ9fv3YialxFRUBqqnxs3w6cPAkASAAAIQBFsWFxd8eATERELUNQkBxWUVs7ETW+khLzQPzjj7K9XTtgyBA55CkuDvf36YNqOw7HAAMyERG1FAsWGMYgG7i6ynYisr7SUuCHH4yB+MgR2e7mBkRHAxMnAnFxQGQk0MYYOXU2KbZ+GJCJiKhl0N+Il5gIXV4eVMHBMhzzBj0i67h8GUhLk2E4NRU4dEgOlXB1BQYNAp56SgbiPn2Atm1tXe09UYQQVt+pRqMRWVlZVt8vNQ5FUdAYnwMiIlvhz7X64zlrmBZ93q5cAdLTjYH4wAEZiJ2dgYEDZRiOiwOiogBHR4t3a8tzpijKPiGEpq7teAWZiIiIiIBr14CMDGMg3rcP0OkAJydgwAAgKUkG4n79ZFsLxoBMRERE1Br9/DOwY4cxEGdmAtXVcnhE//5yisS4OPm1i4utq21SDMhERERErcGNG8DOncZAvHcvcOuWvIGub19g7lwZiAcOlOOKWzEGZCIiIqKW6OZNYPduYyDevRuorAQcHACNBnj1VRmIBw0C3N1tXa1dYUAmIiIiagkqKuRVYX0g3rlTtqlUcqq1l16SgXjwYDk3Md0RAzIRERFRc3Trlhw3vH27fOzcKYdRKArQsyfw4osyEEdHA2q1rattVhiQiYiIiJqDqio5s4Q+EO/YIW+0A4AePYDnn5eBeMgQwMvLtrU2cwzIRERERPaoulrOPawfMpGeLqdiA4Du3Y0r1cXEAN7etqy0xWFAJiIiIrIHOh1w8KDxCnFaGnD1quwLC5OrQuoDcYcOtq21hWNAJiIiIrIFnQ44csQ8EF++LPu6dAGefFIG4thYoGNHm5ba2jAgExERETUFIYBjx4yB+IcfgNJS2RcaCjz6qDEQd+pk01JbOwZkIiIiosYgBHDihDEQp6YCJSWyLygI+P3vjYE4ONiWldJtGJCJiIiIrEEI4PRp80BcVCT7AgKAYcOMgTg0VE7HRnaJAZmIiKg1W7cOSExENQCEhAALFsibwahuQiAUAP7+d2MgLiiQfX5+MgzrA3GXLgzEzYhFAVlRFDWAlQAiAAgAzwohdjVmYURERNTI1q0DJk8GysuhAoC8PPkcYEi+k7w8syvEuQDw3HOAr68MwvpAHBbGQNyMWXoF+QMAW4QQoxVFcQTg2og1ERERUVNITATKy83bystlOwOydO6c+ZCJM2dke/v2QGwspuXnY+mRI3JeYgbiFkMRQtx9A0VpB+AggM6iro1/odFoRFZWlhXKo6agKAos/F9LRNQs8OeahVQqOW72dooipyBrjQoLzQPx6dOy3dNTzj+sHzYRHg6oVPysNYAtz5miKPuEEJq6trPkCnJnACUAViuK0hPAPgAvCyF+vscaiYiIyJaCguSQgdraW4viYhmE9YH4xAnZft99csnmF1+UgfjXv5a/UFCrYElAbgOgN4AZQog9iqJ8AGAegDdNN1IUZTKAyQAQ1Jq+sYiIiJqrBQsMY5ANXF1le0t18aJ5ID52TLZ7eADR0XI8cVwc0KsX4OBgy0rJhiwJyOcAnBNC7Pnl+UbIgGxGCLECwApADrGwWoVERETUOPTjjBMTocvLgyo4uOXNYnHpklyQQx+IDx+W7W5uwODBwNNPy0DcuzfQhpN7kVTnJ0EIUaQoyk+KooQJIU4A+A2AY41fGhERETW68eOB8ePhoCgQZ8/aupp7V1Yml2zWB+KDB+U4axcXYNAgYOxYGYg1GqBtW1tXS3bK0l+VZgBY98sMFrkAnmm8koiIiIgsdPUqkJ5uDMQHDsgbDJ2cgIEDgfnzZSCOipJtRBawKCALIbIB1HnHHxEREVGjun4dyMgwjiPetw+orgYcHYH+/YE335SBuF8/wNnZ1tVSM8XBNkRERGS/ysuBHTuMV4gzM4GqKjk8om9f4LXXZCAeMEAOoyCyAgZkIiIish83bgC7dhkD8Z49wK1bckaJqChg9mwZiAcOlDfaETUCBmQiIiKynYoKYPduYyDetQuorJRzDvfpA8ycKQPxoEFyKjaiJsCATERERE2nshLYu9cYiHfuBG7elKv3RUYCM2bIQDx4sFysg8gGGJCJiIio8dy6BWRlGQPxjh3GhUl69gSmTJGBeMgQQK22aalEegzIREREZD1VVcD+/cZAnJEhZ54AgIgIYNIkYyBu396mpRLdCQMyERERNVx1NZCdbZx2LT1dzk0MAL/6lXGlupgYwMfHpqUSWYoBmYiIiCyn0wGHDhkDcVqaXL0OALp1A8aNk4E4Nhbo0MGWlRI1GAMyERER3ZlOBxw9agzEP/wAXLok+7p0AUaPNgZif39bVkp2TqvVYv78+QAARVGQlJQErVZr26LuQBFCWH2nGo1GZGVlWX2/1DgURUFjfA6IiGyFP9fqz3DOhAB+/NEYiFNTgYsX5UahoTII6wNxYKDtCrYT/Kw1L4qi7BNC1Lk6NK8gExERtWZCAKdOYTIAjB0rA3FxsewLDARGjDAG4pAQm5VJ1JQYkImIiFoTIYCcHPMrxOfPYzkgb7AbOlQG4rg4ecVYUWxbL5ENMCATERG1dGfOGAPx9u3AuXOy3c/PMGSi2wsv4OS5cwzERGBAJiIianny880DcV6ebPfxMY4hjosDwsIMgfjUCy8wHBP9ggGZiIiouSsoMA/EubmyvX17Of/wrFkyEHfvzhBMZAEG5NZs3TogMRHVgLzxYsECYPx4GxdFRER1KioyD8SnTsl2tVoG4pdekoE4IgJQqWxaKlFzxIDcWq1bB0yeDJSXQwXIP79Nniz7GJKJiOzLhQsyEOtD8fHjsr1dO7lk85QpMhD/+teAg4MtKyVqETgPcmsVEmIck2YqOBg4e7apqyEisqpmPzftxYtyQQ59ID56VLa7uwPR0cYxxJGRVgvEzf6c2QjPW/PCeZDp7vLz69dORESN5/Jl80B86JBsd3UFBg8G/vAHGYj79AHa8J9uosbG77LWKiio9ivIQUFNXwsRUWtz5QqQlmYMxNnZcn5iFxdg0CDgnXdkII6KAtq2tXW1RK0OA3JrtWCBYQyygaurbCciIuu6dk0uwqEPxPv3Azod4OQEDBwIaLUyEPftK9uIyKYYkFsr/Y14iYnQ5eVBFRzMWSyIiKzl55+BjAxjIM7KAqqrAUdHoH9/4I03ZCDu3x9wdrZ1tUR0G96kR7zBgIhanCb/uVZeDuzcaQzEe/cCVVVyvHC/fsbFOQYMkH+ts0P8t6BheN6aF96kR0RE1Fhu3gR27TIG4j17gMpKOaNEVJRxYY5BgwA3N1tXS0T1xIBMRGSHtFot5s+fb3ielJQErVZru4Jau4oKGYL1gXjXLtmmUsmZJV55RV4lHjwY8PCwdbVEdI84xIL45yEiO8bvz4a55/NWWQlkZsownJoqh0/cuCGXaY6MNA6ZiI4G7rvPWmXbFD9rDcPz1rxwiAUREZGlqqrkjXT6QJyRYZzlp2dPOetPXJxctc7T06alElHjY0AmIqLWp6oKOHDAGIjT04Hr12VfRAQwaZK8ShwTA7Rvb8tKicgGGJCJiKjlq64GDh40BuK0NODqVdn3q18BTz9tDMS+vraslIjsAAMyERG1OApgHoh/+AEoK5Od3boB48bJQBwbC/j52axOIrJPDMhERNT8CQEcPWoIxBcAoFcv2Xf//cDo0cZAHBBguzqJqFlgQCYiouZHCOD4ceMV4tRUoKRE9oWEYBOAZz/5RN5YFxhow0KpxVq3DkhMRDUAhIRwNdoWhgGZiIjsnxDAqVMyEOtDcXGx7AsMBB56yDj1WkgIJikKnn36aVtWTC3ZunVyZpPycqgAIC9PPgcYklsIzoNMnMORyI612u9PIYDcXPNAfP687PP3l0FYH4g7d5bzE5toteftHvCc1UNIiAzFtwsOBs6ebepqqB44DzIRETUvZ8+aB+KffpLtHTqYB+KuXWsEYqImlZ9fv3ZqdhiQiYjINn76yTwQ66+8+fjIMDxvngzEDzzAQEz2JSio9ivIQUFNXws1CgZkIiJqGufPmwfinBzZ7uUlA/Grr8pAHB7OQEz2bcECwxhkA1dX2U4tAgMyERE1jqIiGYT1gfjkSdmuVssFOaZPl4G4Rw9ApbJlpUT1o78RLzERurw8qIKDOYtFC8Ob9Ig3ZhDZsWb1/VlSYh6If/xRtrdrBwwZYhxD3LMn4ODQqKU0q/NmJ3jOGobnrXnhTXpERNS4SkvlCnX6YRNHj8p2NzcgOhqYOFEG4shIoA3/uSGi5oM/sYiIyDKXLwNpacZAfOiQbHd1BQYPln9ejosD+vQB2ra1ba1ERPeAAZmIiGp35QqQnm4MxNnZcn5iZ2dg0CDgj3+UgTgqCnB0tHW1RERWw4BMRETStWtARoYxEO/fD+h0gJMTMGAAkJQkA3G/frKNiKiFYkAmImqtfv4Z2LHDGIizsoDqajk8on9/IDFRBuL+/QEXF1tXS0TUZBiQiYhai/JyYNcuYyDeuxeoqpI30PXtC8ydKwPxwIFyXDERUSvFiSeJiFqqmzfldGtJSXKaNU9PYOhQYOFCeaU4IQHYskXefLdjh5zHdejQZh2OtVotlF8WGVEUBVqt1rYFEVGzxHmQiXM4Etmxen1/VlTIq8L6K8S7dsk2lUpOtRYXJx+DB8u5iYlM8N+ChuF5a144DzIRUUtXWSnHDesD8c6dwI0bcpnmnj2BF1+UgTg6Wq5eR0REFmFAJiJqLqqqZCDWr1aXkSHHFQNyuebnn5eBeMgQwMvLpqUSETVnLWYMsn7cmf7BcWdE1OxVV6MPACxaBIwYIccQDxgAvPYacO4c8MwzwMaNconnQ4eADz4ARo1iOCYiukctbgwyxwLVH88ZkZ2orpZBVz9kIi0NuHpV9j3wABAbK68Qx8QAHTrYtFRqefhvQcPwvDUvHINMRGTvdDrgyBHzQHz5suzr2hV48kmMS0nB+vPngY4dbVsrEVErwivIxHNG1FSEAI4dMwbiH34ASktlX+fOxivEsbFAp04A+P1JTYeftYbheWteeAWZiMjWhABOnDAG4tRUOV4YAIKCgN//3hiIg4NtWSkREZlgQCYishYhgNOnzQNxUZHsCwgAhg0zzkUcEiKnYyMiIrvDgExE1FBCAGfOmAfiggLZ5+dnDMNxccD99zMQExE1EwzIRET1kZdnHojz82W7r69xDHFcHNCtGwMxEVEzxYBMRHQ3584ZA/H27cDZs7K9fXsZiOfMkYH4V79iICYiaiEYkImITJ0/b1ypbvt2ICdHtnt6ykA8c6YMxOHhgKrFrLVEREQmGJCJqHUrLjYPxCdPyvb77pMLckybJgPxr3/NQExE1EowIBNR61JSIucf1gfiH3+U7R4ewJAhwPPPy0Dcqxfg4GDbWomIyCYYkImoZbt0yTwQHzki293cgMGDgQkTZCDu3Rtowx+JRETEgExELU1ZmVyyWR+IDx2S07G5uMhAPG6cDMQaDdC2ra2rJSIiO8SATETN29WrQHq6MRAfOCADsbMzMHAgMH++DMR9+wKOjrauloiImgEGZCJqXq5fBzIyjPMQ79sHVFfL8DtgAPDWWzIQ9+snQzIREVE9WRSQFUU5C+AagGoAVUIITWMWRURkUF4O7NhhDMSZmUBVlRwe0a8f8NprMhAPGCCHURAREd2j+lxBjhNCXGy0SoiIAODGDWDXLmMg3rMHuHVL3kAXFQXMni0D8cCB8kY7IronWq0W8+fPBwAoioKkpCRotVrbFkVkY4oQou6N5BVkjaUBWaPRiKysrHssrWEURYElx0RGPGdkUxUVwO7dxkC8axdQWSnnHNZoZBiOjZU32Lm727raJsfvTyL7xu/R5kVRlH2WjISw9AqyAPBfRVEEgOVCiBW1vOFkAJMBICgoqD61ElFrUlkJ7N1rDMQ7dwI3b8plmnv3BmbMkKE4Ohpo187W1RIRUStk6bJQg4QQvQE8BGCaoihDbt9ACLFCCKERQmh8fHysWiSRvdBqtVAUxfDgnyEtcOsWVk6ahNcUBf9VFPzs5CTDb1KSnKN4yhTg66/l11lZwOLFQHw8wzEREdmMRUMszF6gKFoA14UQi++0DYdYNC88Z/XHc3YXVVXA/v3GK8QZGXLmCQCHAfTQXyGOiQG8vGxaanPAzxqRfeP3aPNitSEWiqK4AVAJIa798vXvALxthRqJqCWorgays42BOC0NuHZN9nXvLleqi40FYmLwa19fiA8/tGW1REREdbJkDHIHAP9SFEW//adCiC2NWhUR2S+dTq5OZxqIy8pkX1gYMH68DMSxsUCHDjYslIiIqGHqDMhCiFwAPZugFiKyRzodcPSoMRD/8IMcLwwAXboAY8YYA7G/vw0LJSIisg6upEdE5oQAfvzRGIhTU4GLv8zwGBoKjBplDMSBgbark4iIqJEwIBO1dkIAJ0/KQKwPxRcuyL6gIDmjRGysvLEuONiWlRIRETUJBmSi1kYIICfHPBAXFsq+gADgd78zBuLQUDk/MRERUSvCgEzUGpw5Yx6Iz52T7X5+xpXq4uLkmGIGYiIiauUYkIlaovx880CclyfbfX2NYTg2Vs46wUBMRERkhgGZqCUoKDAPxLm5sr19exmEZ8+W/+3enYGYiIioDgzIRM1RYaEMwvpAfOqUbPf0lCvUvfyyDMQREYDK0hXliYiICGBAJmoeLlwwD8THj8v2++4DhgwBpk6VgfjXvwYcHGxYKBERUfPHgExkjy5elAty6APx0aOy3cMDiI4GJk2SgTgykoGYiIjIyhiQiezB5cvGQLx9O3D4sGx3cwMGDwb+7/9kIO7TB2jDb1siIqLGxH9piWzhyhUgLc0YiA8elPMTu7gAgwYBY8fKQBwVBbRta+tqiYiIWhUGZKKmcO0akJ5uDMQHDgA6HeDkBAwcCMyfL6dei4qSbURERGQzDMhEjeH6dWDHDuMY4qwsoLoacHQE+vcH3nxTBuJ+/QBnZ1tXS0RERCYYkImsobwc2LnTeIU4MxOoqpLjhfv1A157TQbiAQPkMAoiIiKyWwzIRA3gBBjD8PbtwJ49wK1bckaJqCi5MEdcnBw+4eZm63KJiIioHhiQiSxRUSFD8C+BuAwAHnxQLsLRpw8wc6YMxIMGyanYiIiIqNliQCaqTWWlHCahv0K8cydw86ZcpjkyEn8DMOvbb+UUbPfdZ+tqiYiIyIoUIYTVd6rRaERWVpbV92sJRVHQGMfUkvGcQQ6P2LfPGIh37JDjigGgZ095dTguTi7S4enJc9ZAPG/1x3NGZN/4Pdq8KIqyTwihqWs7XkFuxbRaLebPnw9AfoMnJSVBq9XatqimUlUlp1rTB+KMDDnzBABERMiV6uLi5DLO7dvbtlYiIiJqUryCTK1DdbVcjEMfiNPTgatXZd+vfmW8QhwTA/j41Lk7fs4ahqvCzDYAABA1SURBVOet/njOiOwbv0ebF15BptZNp5PLNesDcVoaUFYm+7p1A8aNMwZiPz/b1kpUm3XrgMREVANASAiwYAEwfryNiyIiah0YkKll0OmAY8eMgfiHH4BLl2Tf/fcDo0cbA3FAgG1rJarLunXA5MlAeTlUAJCXJ58DDMlERE2AQyyoeRICOH7cGIhTU4GLF2VfSIhxyERsLBAYaPW35+esYXjeLBQSIkPx7YKDgbNnm7oaIroL/lxrXjjEgloWIYCTJ2UQ1gfi4mLZFxgIjBghw3BcnAwXRM1Zfn792omIyKoYkMk+CQHk5JgH4vPnZZ+/PzB0qDEQd+4s5ycmaimCgmq/ghwU1PS1EBG1QgzIZD/OnjVfvvncOdneoYNxuERcHNC1KwMxtWwLFhjGIBu4usp2IiJqdAzIZDs//WQeiPVXzHx8ZBjWB+IHHmAgptZFfyNeYiJ0eXlQBQdzFgsioibEm/So6RQUGIdMbN8O5ObKdi8v80AcHm73gZifs4bheas/njMi+8bv0eaFN+mR7RUVmQfiU6dku1otp1ubMUMG4h49AJXKpqUSERER6TEgk/VcuCDnH9YH4uPHZXu7dnLJ5hdekIG4Z0/AwcG2tRIRERHdAQMyNVxpqXkgPnpUtru5AdHRwDPPyEAcGQm04UeNiIiImgemFrLc5ctyyWZ9ID50SLa7ugKDB8sbiOLigD59gLZtbVsrERERUQMxINOdXbkiA7F+HHF2tpyf2NkZGDQI+OMfZSCOigIcHW1dLREREZFV8M4oMrp2Ddi8GZgzR4ZeLy/gkUeApUuB++4DkpLkkIqyMuD774E33pBBmeGYiIhaGa1WC+WXGZcURYFWq7VtQWRVnOatNfv5Z2DHDuOQiawsoLpaDo/o39847Vr//oCLi62rtSv8nDUMz1v98ZwREVkPp3mjmsrLgV27jIF4716gqkreQNe3LzB3rgzEAwfKccVERERErRADckt286YMxPoxxHv2AJWVcoo1jQZISJCBeNAgwN3d1tUSERER2QUG5JakokKGYH0g3rVLtqlUcqq1l16SgXjwYDk3MRERERHVwIDcnFVWApmZxkC8cydw44ZcprlXL+DFF2Ugjo6Wq9cRERERUZ0YkJuTqip5I50+EGdkyHHFAP6/vbuLtawuzwD+vMxgYCiKkUGnjMyo8asYgmQHW9GRQVtrJdp6RTM0pibSNKRBe2Fq58KDCRdNTONdEwJtbZxqrB8xIdZoIuPHhZozIIqiVi2DgMpUtBRpqg5vL/bGWZgZ55zDOay9z/x+ycmcvfZizTNvzsfD2v+1di66KHnrW6eFeM+e6R0oAABYNQV5nh09mtx++7GL6j7/+eThh6fPXXjhsXeqe9WrknPPHTcrAMAmoSDPk6NHp+9O91gh/tznkocemj73ohclV189LcSXX56cd96oUQEANisFeUyPPprceefjC/FPfjJ97vnPT666alqGL7882bFjzKQAAKcMBfnJ1J184xvHCvFnP5v8+MfT55773ORNbzpWiHfuHDMpAMApS0HeSN3Jt751rBAfPJgcOTJ97oILkiuvnC6Z2Lt3+hgAgNEpyOupO/nOdx5fiH/4w+lz55+fvPa1xwrxc54zalQAAI5v8xTkAweS/ftzNEl2705uuCHZt29j/87u5HvfO3bbtYMHk/vumz63Y0dyxRXT5RJ79ybPe970/sQAAMy1zVGQDxxIrrkmeeSRnJYkhw9PHyfrX5IPHz52hvjWW5Pvf3+6/bzzjt1hYu/e5AUvUIgBABZQdfe6H3QymfTy8vK6H/eEdu+eFtdft2tXcvfdT+zY9977+EL82PHOPffYBXV79yYvfrFCfAqpqmzE985mZ26rZ2YA66eqDnX35GT7bY4zyPfcs7rtv8n99x9bMnHrrcl3vzvd/vSnT8vw298+LcQXXpicdtpaEwMAMKc2R0G+4ILjn0FeyZ0hfvSjxxfib397uv1pT5u+Q921104L8UUXKcQAAKeAzVGQb7jhV2uQf2Xbtun2X3fkyPT+w48V4rvumm4/++xkz57pcS6/PLn44mTLliclPgAA82NzFOTHLsTbvz+PHj6c03btOnYXiwcffHwhvvPO6b5nnZW88pXJm988PUN8ySXJ1s0xDgAA1m5zXKQ3cE5Vfvrxjx8rxF/96vR2bGeembziFcfuNDGZJKefPkpGFp8Lp9bG3FbPzADWz6l1kV6S3HFH8pa35MEkeeMbkzPOSF7+8uTd754W4ksvTZ7ylJFDAgAw7zZPQX7Ws5Kzz871Sa4/eDB52cumJRkAAFZh89yW4ZnPTA4ezLuT6d0nlGMAANZg8xRkAABYBwoyAAAMKMgAADCgIAMAwICCDAAAAwoyAAAMKMgAADCw4oJcVVuq6vaqumUjAwEAwJhWcwb5uiR3bVQQYBM7cCDZvTtHk2T37uljAJhTKyrIVbUzyeuT3LSxcWDOKXqrd+BAcs01yeHD0x84hw9PH5sdAHNqpWeQ35vkHUkePdEOVXVNVS1X1fKRI0fWJRzMFUVvbfbvTx555PHbHnlkuh0A5tBJC3JVXZnkge4+9Jv26+4bu3vS3ZPt27evW0CYG4re2txzz+q2A8DIVnIG+bIkb6iqu5N8MMkVVfX+DU0F80jRW5sLLljddgAY2UkLcne/s7t3dvfuJFcl+Ux3X73hyWDeKHprc8MNybZtj9+2bdt0OwDMIfdBhpVS9NZm377kxhuTXbumFzHs2jV9vG/f2MkA4Liqu9f9oJPJpJeXl9f9uCtRVdmIfxMkmV6Qt39/Hj18OKft2jUtx4reivn+XD0zA1g/VXWouycn22/rkxEGNo19+5J9+7KlKn333WOnAQA2gCUWAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMsAcWlpaSlUlSaoqS0tL4wYCOIVUd6/7QSeTSS8vL6/7cVeiqrIR/yYY8nW2NuYGwJiq6lB3T062nzPIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAObpiAvLS2lqpIkVZWlpaVxAwEAsJCqu9f9oJPJpJeXl9f9uDAvqiob8b2z2ZkbAGOqqkPdPTnZfpvmDDIAAKwHBRkAAAYUZAAAGFCQAQBgQEEGAIABBRkAAAYUZAAAGFCQAQBgQEEGAIABBRkAAAYUZAAAGFCQAQBgQEEGAIABBRkAAAZOWpCr6oyq+nJV3VFVX6+q65+MYAAAMIatK9jn/5Jc0d0PV9XpSb5QVf/e3V/c4GwAAPCkO2lB7u5O8vDs4emzj97IUAAAMJYVrUGuqi1V9ZUkDyT5dHd/6Tj7XFNVy1W1fOTIkfXOCSywpaWlVFWSpKqytLQ0biAA+A1WVJC7+2h3X5xkZ5JLq+olx9nnxu6edPdk+/bt650T5oKitzZLS0vp7l99mBsA86ymKyhW8R9UvSvJz7r7PSfaZzKZ9PLy8hPNBgAA66aqDnX35GT7reQuFtur6pzZ52cmeU2Sbz7xiAAAMH9WcheLHUneV1VbMi3UH+ruWzY2FgAAjGMld7H4apKXPglZAABgdN5JDwAABhRkAAAYUJABAGBAQQYAgAEFGQAABhRkAAAYUJABAGBAQQYAgAEFGQAABhRkAAAYUJABAGBAQQYAgAEFGQAABqq71/+gVUeSHF73A6/MuUn+a6S/e1GZ2eqZ2dqY2+qZ2dqY2+qZ2dqY2+qNObNd3b39ZDttSEEeU1Utd/dk7ByLxMxWz8zWxtxWz8zWxtxWz8zWxtxWbxFmZokFAAAMKMgAADCwGQvyjWMHWEBmtnpmtjbmtnpmtjbmtnpmtjbmtnpzP7NNtwYZAACeiM14BhkAANZMQQYAgIFNUZCr6h+r6oGqunPsLIukqp5dVbdW1V1V9fWqum7sTPOuqs6oqi9X1R2zmV0/dqZFUVVbqur2qrpl7CyLoqrurqqvVdVXqmp57DyLoKrOqaoPV9U3Zz/bfm/sTPOuql44+xp77OOhqnrb2LnmXVW9ffZ74M6q+kBVnTF2pkVQVdfNZvb1ef462xRrkKtqT5KHk/xLd79k7DyLoqp2JNnR3bdV1dlJDiX54+7+xsjR5lZVVZKzuvvhqjo9yReSXNfdXxw52tyrqr9OMkny1O6+cuw8i6Cq7k4y6W5vQrBCVfW+JJ/v7puq6ilJtnX3T8fOtSiqakuS+5K8rLvHesOvuVdV52f68/93uvt/q+pDST7R3f88brL5VlUvSfLBJJcm+XmSTyb5y+7+j1GDHcemOIPc3Z9L8uDYORZNd/+gu2+bff4/Se5Kcv64qeZbTz08e3j67GPx/y9zg1XVziSvT3LT2FnYvKrqqUn2JLk5Sbr758rxqr06yXeV4xXZmuTMqtqaZFuS+0fOswhenOSL3f1Id/8yyWeT/MnImY5rUxRknriq2p3kpUm+NG6S+TdbKvCVJA8k+XR3m9nJvTfJO5I8OnaQBdNJPlVVh6rqmrHDLIDnJjmS5J9my3luqqqzxg61YK5K8oGxQ8y77r4vyXuS3JPkB0n+u7s/NW6qhXBnkj1V9Yyq2pbkj5I8e+RMx6Ugk6r6rSQfSfK27n5o7DzzrruPdvfFSXYmuXT2khEnUFVXJnmguw+NnWUBXdbdlyR5XZJrZ8vJOLGtSS5J8g/d/dIkP0vyN+NGWhyzJSlvSPJvY2eZd1X19CRvTPKcJL+d5KyqunrcVPOvu+9K8ndJPp3p8oo7kvxy1FAnoCCf4mbraD+S5EB3f3TsPItk9tLtwSR/OHKUeXdZkjfM1tN+MMkVVfX+cSMthu6+f/bnA0k+lum6PU7s3iT3Dl7V+XCmhZmVeV2S27r7R2MHWQCvSfKf3X2ku3+R5KNJXj5ypoXQ3Td39yXdvSfT5bFzt/44UZBPabMLzm5Ocld3//3YeRZBVW2vqnNmn5+Z6Q/Jb46bar519zu7e2d378705dvPdLczLSdRVWfNLp7NbJnAH2T68iQn0N0/TPL9qnrhbNOrk7joeOX+NJZXrNQ9SX63qrbNfpe+OtPreDiJqjpv9ucFSd6UOf2a2zp2gPVQVR9IcnmSc6vq3iTv6u6bx021EC5L8mdJvjZbU5skf9vdnxgx07zbkeR9syu9T0vyoe522zI2wjOTfGz6uzdbk/xrd39y3EgL4a+SHJgtF/hekj8fOc9CmK0H/f0kfzF2lkXQ3V+qqg8nuS3TJQK3ZwHePnlOfKSqnpHkF0mu7e6fjB3oeDbFbd4AAGC9WGIBAAADCjIAAAwoyAAAMKAgAwDAgIIMAAADCjIAAAwoyAAAMPD/d4ugrmTFwm0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.errorbar(x, y, ey, fmt='ro', ecolor='k', elinewidth=1, capsize=2, capthick=1)\n",
    "ax.plot(x, fit_function(x, *minuit.args), '-r')\n",
    "plt.close()\n",
    "\n",
    "d = {'alpha0':   [alpha0_calc, sigma_alpha0_calc],\n",
    "     'alpha1':   [alpha1_calc, sigma_alpha1_calc],\n",
    "     'Chi2':     Chi2_calc,\n",
    "     'ndf':      Ndof_calc,\n",
    "     'Prob':     Prob_calc,\n",
    "    }\n",
    "\n",
    "text = nice_string_output(d, extra_spacing=2, decimals=3)\n",
    "add_text_to_ax(0.02, 0.95, text, ax, fontsize=14)\n",
    "fig.tight_layout()\n",
    "fig"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "if (save_plots) :\n",
    "    fig.savefig('AnalyticalLinearFit.pdf', dpi=600)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 2
   },
   "source": [
    "---------------------------------------------------------------------------------- \n",
    "\n",
    "\n",
    "# Questions:\n",
    "\n",
    "1) Do the analytical result(s) (Calc) agree with the numerical one(s) (Fit)?\n",
    "\n",
    "### Advanced questions:\n",
    "\n",
    "2) Can you measure/determine the difference in speed between the two methods?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "executable": "/usr/bin/env python",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}