{
 "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 2018"
   ]
  },
  {
   "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",
    "from probfit import BinnedLH, Chi2Regression, Extended # Helper tool for fitting\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 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 = 1       # 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 (a simple line) and fit these both **analytically** and **numerically**.\n",
    "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/Stat2018/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!\n",
    "\n",
    "The whole thing is contained in a loop, allowing us to repeat the generation and fitting, so that we can see the result on more than just a case-by-case. For this reason, we save the result of the fit (in pre-defined arrays)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "lines_to_next_cell": 2
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  Fit: a0= 4.124+-0.726  a1=0.283+-0.129  p=0.6994\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\"  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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzs3XtclFX+B/DPA8gdGxQQkauXSMEQQfMOmpuWa1nZxcsmaeulNHPxVrox9sq01JbaLiqudNFcN3M3t35amSJ4BxVvpKYIrMhFBTRFUJjz++PEDCOowzjDzMDn/XrNK+ecZ575ztOIH4/nOUcRQoCIiIiIiCQ7SxdARERERGRNGJCJiIiIiGphQCYiIiIiqoUBmYiIiIioFgZkIiIiIqJaGJCJiIiIiGphQCYiIiIiqoUBmYiIiIioFgZkIiIiIqJaHMxxUi8vLxEcHGyOUxMRERERGeXAgQMXhRDedzvOLAE5ODgYGRkZ5jg1EREREZFRFEXJNeQ4TrEgIiIiIqqFAZmIiIiIqBYGZCIiIiKiWhiQiYiIiIhqYUAmIiIiIqqFAZmIiIiIqBYGZCIiIiKiWhiQiYiIiIhqYUAmIiIiIqqFAZmIiIiIqBYGZCIiIiKiWhiQiYiIiIhqYUAmIiIiIqqFAZmIiIiIqBYGZCIiIiKiWhiQiYiIiIhqYUAmIiJqxtRqNRRF0T7UarWlSyKyOEUIYfKTRkdHi4yMDJOfl4iIiMxDURSYIxMQWRNFUQ4IIaLvdpxNjSDPnz9f+zfc8PBwo86hVquNfi0RERERNX02FZDnzp2LgoICxMfHW7SOmzdvYurUqWjVqhVatmyJF198EdeuXTPb+02fPh2KomDDhg0Net1//vMfDBkyBJ6enlAUBTk5OXr9SUlJ6NGjB1q2bAlvb28888wzdY6pkZubi/vuu69R/nIxd+5c9OvXz6TnvHz5MsaMGQMPDw94eXkhPj4e1dXVBr/+1n+CrHmEhYUBAKqqqhAfH4/Q0FC4uroiMDAQM2fORHl5ud55srKy8Nhjj8HDwwMtW7bEww8/bNLPSURERPfOpgKyu7s7fH194e7ubtE63nzzTXz99dfYsGEDtmzZgpSUFLz22mtmea+ffvoJhw8fNuq1ly9fRp8+ffDXv/613v6dO3fipZdewu7du7F161aUlJRg6NChqKqq0jtOo9Fg3LhxiIiIMKqOhtb86aefYu7cuSY97+TJk3HgwAH8/PPPWLNmDZKTk7F48WKDXz9z5kwUFBToPe6//348/fTTAIDKykocPXoU77zzDg4fPozk5GT861//wrRp07TnKC4uRkxMDFq2bInt27fj4MGDmDx5skk/JxEREZmAEMLkj6ioKGGMzZs3i/79+wuVSiVcXFzEoEGDRGZmZp3jEhISRFhYWL3nACASEhJERESEcHJyEkOHDhVFRUV1XpuYmCh8fHyEt7e3SExM1PZXVVWJCRMmiJCQEOHo6CgCAgJEQkKCqK6u1va3atVKfPDBB9rXrFmzRjg5OYmrV68a9blv59KlS6Jjx47i1KlTAoD4+uuvjTpPenq6ACDOnj17x+MOHjwoANS55osXLxYvvfTSHa/73Zw9e1bIr9udLV68WISHhwuNRmPU+9Tn4sWLws7OTnz77bfatrffflv4+/sbfc59+/YJRVFEdnb2bY95//33hUql0j5PSEgQnTp10n6XiIisiSE/o4lsHYAMYUCWtaoR5MLCQowbNw47d+7EoUOH4Ofnh2HDhuHmzZsNOs9HH32EhQsXYt++fSgqKtIbxQOA7Oxs/Prrr0hLS8OUKVMQHx+Pc+fOAQCqq6thb2+P5ORknDhxAp9++ikSExOxYsUK7WtLSkrQv39/7fkGDBigHUE0pcmTJ2PSpEno1KmTSc97O2VlZQAAlUqlbcvMzMSKFSuwdOlSs79/ZWUlEhMTMXv2bCiKYrLzHjx4EBqNps7/s3PnzqGwsNCocyYlJSE2NhYhISG3PaasrAyenp7a5ykpKYiNjUVcXBx8fHzQvXt3/Otf/zLq/YmIiMh8rCogx8XFYcKECQgLC0NoaCjeeust5OfnNzh4vvjiixg2bBgiIiLw9ttv45tvvtGGPwCws7PDsmXLcP/99+ONN96ARqPBoUOHAACOjo5YsWIFYmJiEBISgmHDhuGJJ57A5s2bAQAXLlwAAHh5eeGZZ57BQw89BC8vL70+U/jyyy+Rm5uLGTNmmOycd1JdXY0333wTI0eORFBQEACgoqICY8eOxYcffoj77rvP7DV89tlncHJywqhRo0x63gsXLsDe3h4qlQo9evTAc889d0//z65evYp//vOfmDBhwm2POX/+PD766CO9+fIFBQVYv349vLy88MMPP+D555/H888/j/T09IZ/KCIiIjIbqwrIZ86cwejRo9G+fXt4eHiga9euAGQgaYjaN5KFhYWhurpa7+Yzf39/ODk5AQCcnJzg6uqKkpISbf/y5csRHR0Nb29vuLu7Y926dfXW0LZtWwQGBjaoNkP873//Q3x8PFavXg17e3uTn78+06ZNw6VLl5CUlKRte+ONN9CtWzf88Y9/NOqceXl5cHd3h7u7u/Zmtprn7u7uyMvL0x6r0WiwdOlSxMfHw8HBoc65Jk+erPfahqq5qS4gIAB+fn5GfZ4a69atg729PZ566ql6+69evYoRI0Zg6NCheOWVV7TtGo0G3t7eWLZsGSIjIzF79mz06NEDa9euvad6iIiIyLTqJhELGj58ONq0aYPVq1ejXbt2KC4uRr9+/aDRaO753KLW2o71BbCa/vXr12P69OlYtmwZYmJi4OLignnz5qGoqAgA4O3tDQC4ePEiPvzwQwAy0Nbuu1cHDhzAxYsXERUVpdc+evRofPXVV9i4caNJ3qfGnDlzsGXLFqSlpelNr9i2bRuysrK0q2dUVVVBo9HA2dkZu3fvRvfu3e94Xj8/P2RmZgIA8vPzERsbq31e019jw4YNKCsru+2o7FtvvYWZM2ca9fm8vb1RVVWF0tJS7bVLS0vT9jVUUlISRo0aBRcXlzp9169fx/Dhw+Hj44PPPvtMr8/LywsqlUpv+kj79u2Rn5/f4BqIiIjIfKwmIF+6dAm//PILPvnkE8TGxgIATp48We+x7u7uuH79+m3Pdfz4cb1f29nZITg42KA6du7cid69e2Pq1KnatuzsbLi5uQGQgcbT0xNpaWmIjIwEAKSmpsLJyUk74n2vBg8ejKysLL22zp07Y+nSpXVGLSsqKlBYWAiVSqUXbg2lVquxbt06pKamol27dnp9//nPf1BRUaF9/tFHH+HHH3/Epk2bDLqeDg4O6Nixo/bXALTPb/Xuu+9i2rRpcHV1rbffx8cHPj4+hnykOiIjI2FnZ4e0tDQ8/vjjAKD9vL6+vnrHlpWVoaysDL6+vnB2dq5zrsOHDyM9PR0ff/xxnb7Kyko8+eSTcHBwwIYNG9CiRYs6dfz88896bXl5eejRo4dRn4uIiIjMw2qmWHh6esLLywvJyck4c+YMfv75Z8yfP7/eY6OionD27Fl89913KCwsrLPW7OrVq7F582YcPnwY8+fPx1NPPaV3s9Sd3H///Th06BC2bt2KU6dOYe7cuThx4oS2397eHhMnTsTChQuxbds27N69G/Pnz8ef/vQnbYi+V+7u7njggQf0HoAccfX399c7du/evQgJCUFiYmKd85SUlCAzMxOnTp0CINfgzczM1E4XWbx4MT744AN89dVXcHZ2RmFhIQoLC7V/+QgODtarwcvLC46OjnjggQfqDY/G+umnn3Dy5Em9v5SYUs188dmzZ2P//v344YcfsGzZMkyZMqXOsYmJiQgJCcHevXvrPVdSUhLCw8PrhNqbN29i5MiRKCoqwsqVK1FWVqa9njXrLY8fPx6nT59GQkICTp8+jeXLl2PPnj0YPXq06T80ERERGc+QpS4a+jB2mbdt27aJ8PBw4eTkJLp16yY2bdokAIjt27fXOXbWrFnC09NTABB/+9vfai/fIf7617+KsLAw4eTkJIYMGSIKCwu1/fUtVebm5iaSk5OFEEJUVlaKCRMmCJVKJVQqlZg2bZqYNGmSiImJ0R5fWVkpXnnlFaFSqYSHh4cYN26cyZd4uxVus8zb9u3btUvb3So5OVkAqPOouZ5BQUH19tdci1vdyzJvdzJo0CDx2muvmfy8tZWVlYlRo0YJNzc30apVK/GXv/xFVFVV1TkuISHhtt+58vJyoVKpxPvvv1+nr2YZu/oetZfXW79+vQgNDRVOTk6iS5cuYt26dab8mERERgOXeaNmAAYu86YIM+y7Hh0dLTIyMkx+XkMoioKvv/4aI0eOtMj7U8Okp6ejb9++yM7OrjM6TkREjUdRFJgjExBZE0VRDgghou92nNVMsaDmqaSkBKtWrWI4JiIiIqth0E16iqLMAPAS5D8ZHwXwohCi4s6vIrq7IUOGWLoEIiIiIj13DciKorQD8CqALkKI64qi/AvA8wA+M3NtRuE/DxERERHRvTB0ioUDABdFURwAuAI4b76SiIiIiIgs564BWQiRD2ApgDwABQAuCyF+NHdh9Zk/f752R7Tau+U1hFqtNvq1RERERNT03TUgK4riCeAJACEA/AC4KYoytp7jJiqKkqEoSsaFCxdMXymAuXPnoqCgAPHx8WY5v6FqNpxo06YNFEVBSkqKWd5nwYIF8PX1hZubG0aMGIHi4mKDX1taWorJkycjJCQELi4u6NChAxYuXKi3K+GNGzcQHx8Pf39/uLq6olu3bvj222+1/SdPnsTAgQPh4+MDV1dXPPjgg1izZo1JP2N95s6di379+pn0nJcvX8aYMWPg4eEBLy8vxMfHa9cnNtSiRYvQu3dvODs719kopaqqCvHx8QgNDYWrqysCAwMxc+ZMvTW6b9y4gTfeeANBQUFwdnZGVFQUUlNTTfHxiIiIyIQMmWIxGMBZIcQFIcRNABsB9Ln1ICHESiFEtBAi2lRbLt/K3d0dvr6+cHd3N8v5DXXlyhWEh4fjvffeM9t7rFixAu+99x4++eQTpKWlIT8/H2PGjDH49UVFRbh06RI+/vhjHD9+HO+//z6WLFmChQsXao9ZvHgxvvzyS3z++ec4fvw4RowYgZEjR+LXX38FIHe/Gzt2LLZu3YqsrCxMnz4d48ePx9atW03+eWtcvnwZn376KebOnWvS806ePBkHDhzAzz//jDVr1iA5ORmLFy9u0DnKy8vx7LPP4oUXXqjTV1lZiaNHj+Kdd97B4cOHkZycjH/961+YNm2a9pj33nsPSUlJWLFiBY4dO4aYmBg89thjOH+eM5aIiIisyt0WSgbwEIDjkHOPFQCfA5h2p9cYu1HI5s2bRf/+/YVKpRIuLi5i0KBBIjMzs85xd9qwAr9vmhERESGcnJzE0KFDRVFRUZ3XJiYmCh8fH+Ht7S0SExO1/VVVVWLChAkiJCREODo6ioCAAJGQkCCqq6vrvNeFCxduu6nEvYqIiBAzZszQPt+5c6cAIH799Vejz/nqq6+Kbt26aZ8PGzZMxMXFaZ9XVVUJe3t7sWHDhtueo3v37uL1119v8HvXbKRxN4sXLxbh4eFCo9E0+D1u5+LFi8LOzk58++232ra3335b+Pv7G3W+JUuWiKCgoLse9/777wuVSqV93rNnTzFr1izt8+rqauHl5SWWLFliVB1ERKZkyM9oIlsHAzcKMWQO8j4AGwAchFzizQ7AStNHdaCwsBDjxo3Dzp07cejQIfj5+WHYsGG4efNmg87z0UcfYeHChdi3bx+Kior0RvEAIDs7G7/++ivS0tIwZcoUxMfH49y5cwCA6upq2NvbIzk5GSdOnMCnn36KxMRErFixwmSf825qRiP79++vbevduzdatGiBe9mApaysTG/L7QEDBiAtLQ15eXkQQuCf//wnPDw80KdPnX8ggBACP/30E7KyshAVFWV0DXdSWVmJxMREzJ49G4qimOy8Bw8ehEaj0bueAwYMwLlz51BYWGiy97nVrdf7xo0bcHJy0j63s7ODo6MjDh06ZLYaiIiIyAiGpOiGPowdQb5Vdna2ACAOHDig1363EeSZM2dqn3///ffC3t5elJaWal/r5uYmKioqhBBCVFRUCEVRxKZNm25bxwsvvCCGDx9ep91cI8j5+fkCgEhNTRXx8fEiJCREXLt2Tfj4+IgPP/zQqHMeOXJEODo6iu+++07bptFoxPz584WiKMLBwUF4e3uL3bt313lt7969RYsWLYSjo6NYvny5Ue9vyAjy8uXLRVBQkLh586ZR73E7a9euFfb29kKj0Yjo6Gjx7LPPiqysLAFAHDlypMHnM2QEOT8/X7Rq1Up89NFH2rbp06eL9u3bizNnzoiqqirx8ccfCwcHB/GHP/yhwTUQEZna3X5GEzUFMNUIcmM6c+YMRo8ejfbt28PDwwNdu3YFAFy9erVB56m9SkVYWBiqq6uRk5OjbfP399eO5Dk5OcHV1RUlJSXa/uXLlyM6Ohre3t5wd3fHunXrGlyDqXh7eyMwMBAODgbt6VKvwsJCjBgxAvHx8Rg2bJi2fcOGDVi7di3+/e9/IyMjA5MnT8aIESP0rhUArF+/HgcOHMCiRYvwxhtvYPfu3Qa9b15eHtzd3eHu7o6wsDAA0D53d3dHXl6e9liNRoOlS5ciPj6+3s86efJkvdc2VM3qJwEBAfDz82vw6xvi6tWrGDFiBIYOHYpXXnlF265WqxEWFoaOHTvCyckJ33zzDYYOHWrS0XIiIiK6d8anLjMYPnw42rRpg9WrV6Ndu3YoLi5Gv3799FZeMJaotYFIfQGspn/9+vWYPn06li1bhpiYGLi4uGDevHkoKiq65xoM1bp1a9jZ2eHixYuYM2cO5syZA41Gg9LSUjT0BsiLFy9i8ODBGDx4MN555x29vlmzZmHmzJl44oknAAARERH4/vvvkZSUpHczX0BAAAICAtC1a1ccP34cb7/9Nv7v//7vru/t5+eHzMxMAEB+fj5iY2O1z2v6a2zYsAFlZWWYMGFCved66623MHPmTMM/eC3e3t6oqqpCaWkpNm7cCECuRFLTZ0rXr1/H8OHD4ePjg88++0yvT6VSYdOmTSgvL8eVK1fg6+uLvn37IjQ01KQ1EBER0b2xmoB86dIl/PLLL/jkk08QGxsLQC4zVh93d3dcv379tuc6fvy43q/t7OzqLMt1Ozt37kTv3r0xdepUbVt2djbc3NwMer0pODk5ITw8HGlpaXjyyScBAHv27MHNmzcRHR2td2xFRQUKCwuhUqmgUqn0+kpLS/HII48gKioKy5cvr/M+ZWVlsLPT/0cEBweHO15be3t7XLt2zaDP4eDggI4dO2p/DUD7/Fbvvvsupk2bBldX13r7fXx84OPjY9D73ioyMhJ2dnba5fkAIDU1Fe3atYOvr6/esWVlZSgrK4Ovry+cnZ0b9D6VlZV48skn4eDggA0bNqBFixb1Hufq6gpXV1fk5uZi7969mDRpklGfi4iIiMzDagKyp6cnvLy8kJycjICAAOTk5GD+/Pn1HhsVFYWzZ8/iu+++Q3R0NFq2bKkXrFavXo2BAwfCz88P8+fPx1NPPaV3s9Sd3H///fjiiy+wdetWBAYGYvXq1Thx4oTejWlXr17F6dOnUVZWBgA4ffo0VCoVAgMD0apVq3u4Cjo1Nw/269cPwcHBeO211zB48OA6AXPv3r0YOHAgEhISoFarte1XrlzBkCFD4O3tjUWLFumNgNeEwqFDh+K9995DaGgoQkJC8O233yI9PR2LFi0CACQnJ0Oj0aBnz55wdXXFtm3b8MUXX+Ddd981yWes8dNPP+HkyZN6fykxJS8vLzzzzDOYPXs2fH19UVpaimXLltW7nnZiYiIWLFiA7du3a/+iViMvLw8lJSU4f/48bty4oR0N79atG27evImRI0eiqKgIGzdu1H43ADlKbW9vj+LiYmzZsgV9+/bFhQsX8OqrryI0NBTPPfecWT43ERERGcmQicoNfRh7k962bdtEeHi4cHJyEt26dRObNm267U1ws2bNEp6engKA+Nvf/lZ78rX461//KsLCwoSTk5MYMmSIKCws1PbXd4Ofm5ubSE5OFkIIUVlZKSZMmCBUKpVQqVRi2rRpYtKkSSImJkZ7/Pbt2wWAOo+ac5hKQkKC8PHxES4uLuKJJ57Q+xy31pKQkFBve32PGqWlpWLy5MnCz89PuLi4iK5du4q1a9dq+9evXy+ioqJEy5YthbOzs+jcubN4//33TboEmxBCDBo0SLz22msmPeetysrKxKhRo4Sbm5to1aqV+Mtf/iKqqqrqHJeQkHDb79y4ceNuez1rbkKs73H27FkhhBBFRUWie/fuwtnZWahUKvH888+L/Px8c35sIiKDgTfpUTMAA2/SU0StubmmEh0dLe5lObJ7oSgKvv76a4wcOdIi708Nk56ejr59+yI7Oxv+/v6WLoeIqNlSFAXmyARE1kRRlANCiOi7HWdVq1hQ81NSUoJVq1YxHBMREZHVsJo5yNQ8DRkyxNIlEBEREelpcgGZ/zxERERERPeCUyyIiIiIiGphQCYiIiJqILVard2lVVEUvaVWyfbZVEBOSUmBoii4ePHiXY9Vq9V6W04TERERmYpardZO6xRCMCA3MVYXkE+fPo2nnnoKKpUK7u7u6Nu3Lw4ePNjg88ycORM7duyo056UlIQePXqgZcuW8Pb2xjPPPIOcnBwTVE62oGY3QkP/olVbWVkZJk6cCB8fH7i4uKBbt244deoUACAnJ0dvJKH2Iz09HQAQGxtbb/8rr7xi8s9JRERExrOqm/QKCgrQp08fhIWF4bvvvkPr1q2xdetWFBQUNPhc7u7ucHd3r9O+c+dOvPTSS+jbty+qq6vxl7/8BUOHDsWxY8e02yFT07VgwQKjtg0XQmDEiBEoLS3FunXr0L59e5w+fVp7roCAgDrf008++QRffPGFdnvwjRs34saNG9r+4uJiREZG4umnn76HT0RERESmZlWJcNGiRaisrMR3332nDR6dO3euc9y2bdvw5ptv4vz583jqqaewatUqbbhdvHgxXn/9dQBAWFgYjh07pvfazz//XO/50qVL0b17dxw/fhwRERHm+FhkJXbv3o3//ve/WLJkCVJTUxv02h07dmDXrl04deoUQkJCAED7XwCwt7fXbuFd45tvvkFcXBwURQGAOtuQf/HFFwgKCsLAgQON+ThERERkJlY1xWLLli149NFH7zrCt2rVKmzYsAFr1qzBmjVrsGHDBm3ftGnTUFBQgPj4eIPes6ysDACgUqmML5ys3tWrVxEXF4eVK1fC0dGxwa9PSUlBWFgY1q9fD39/f4SGhiIhIQFVVVX1Hr9r1y788ssvePHFF297zlWrVukFaCIiIrIOVhWQ8/LyDNpR7c0330R4eDgef/xx9OzZUzvHEwDc3Nzg6+tb7/SKW1VXV+PNN9/EyJEjERQUdE+1k3WbPn06Hn/8cTz00ENGvb6goABnz57F1q1bsXHjRrzzzjv48MMPkZiYWO/xSUlJePjhh2/7vUpJScGZM2cQFxdnVD1ERERkPlY1xcJQHTt21P66VatWKCkpMeo806ZNw6VLl/Df//7XVKWRFdq0aRPS0tJw+PBho8+h0Whw5coVJCcnIyAgAD179sT+/fuxZs0azJw5U+/Yy5cv4+uvv8aqVatue76kpCQMHjwYgYGBRtdERERE5mFVI8iBgYHIz8+/63G33kxnzO55c+bMwZYtW/DTTz9xekUTt23bNmRnZ8PT0xPOzs545JFHAAD+/v63HQG+lZeXF5ydnREQEKBta9++fb3f17Vr18LJyQlPPvlkvecqKSnBN998g/HjxxvxaYiIiMjcrCogDxkyBJs3b0Z5eblZ30etVmPdunXYtm0b2rVrZ9b3Ist74403cOzYMWRmZiIzM1M7spuSkoIXXnhB79iysjLk5OSgoqJCrz0yMhIVFRU4f/68ti0vL6/e709SUhJGjRoFZ2fneuv58ssv4ebmhhEjRtzrRyMiIiIzsKqAPHfuXLRo0QLDhw/H7t27ceLECXz88cf4/vvvDT5HYWEhCgsLcfXqVVRVVWmf1yyvtXjxYnzwwQf46quv4OzsrO2/fv26uT4WWZiPjw8eeOAB7aNmWkPHjh3rrCyRmJiIkJAQ7N27V6992LBhaNu2LSZPnoxffvkFP/74I1auXIkxY8boHZeRkYHMzMw7jg4nJSVh9OjRcHJyMtEnJCIiIlOyqjnI7dq1w86dOzFnzhw8+uijqKqqQteuXfHJJ58YfI62bdvW+3z79u2IjY3F8uXLUVZWhv79++sdl5yczBum6Lbc3NywefNmTJ06FZGRkfDx8cGkSZMwY8YMveOSkpLw4IMPIioqqt7z7NmzB8ePH8eXX37ZGGUTERGRERRj5u/eTXR0tMjIyDD5eYmIiMg8FEUx6p6e5o7XzbYoinJACBF9t+OsaooFEREREZGlMSATEREREdXCgExEREREVAsDMhERERFRLTYfkGNjYzF16tQ67YcOHUK3bt3QokULKIpigcqIiIiIyBbZfEC+nblz5yIgIABnzpxBQUGBpcshQmxsLBRFqfN45ZVXLF0aERER1WJV6yCb0pkzZzBp0iTtphBElrZx40bthjUAUFxcjMjISDz99NMWrIqIiIhuZVUjyHFxcfjjH/+IuXPnQqVSISAgAF9//bXeMZ9//jkCAgLg5uaGKVOmQKPR6PUHBwdDURScOXMGs2fP1o7SEVlaq1at4Ovrq31s2bIFQUFBGDhwoKVLIyIiolqsKiADcsc7Dw8PpKenY+jQoZg4caJ2G+isrCyMHz8e8fHxOHjwIABg9+7deq9PT09HQUEB/P398eabb6KgoIBTLMgqrVq1CnFxcfwLHBERkZWxuoDs5+eHefPmoVOnTpg1axbKyspw+vRpAMA//vEPREdH47XXXkNoaCj+9re/wd3dXe/13t7e8PX1hb29PTw8PLSjdUTWJCUlBWfOnOH25kRERFbI6gJyhw4dtL9u1aoVAKCkpAQAcPr0aXTp0kXb7+zsjI4dOzZugUQmkJSUhMGDB3OOPBERkRWyupv0HBxeFuBsAAAgAElEQVTqlsQ9zqkpKSkpwTfffIPPP//c0qUQERFRPaxuBPlOOnbsiKysLO3ziooK7fQLIlvx5Zdfws3NDSNGjLB0KURERFQPmwrI48ePR0ZGBj744AOcPHkS8fHxuHr1qqXLImqQpKQkjB49Gk5OTpYuhYiIiOphUwE5LCwM//jHP7B06VJ0794dGo0Gffr0sXRZRAbbs2cPjh8/jvHjx1u6FCIiIroNxRzze6Ojo0VGRobJz0tERETmoSgK7/kxAq+bbVEU5YAQIvpux9nUCDIRERERkbkxIBMRERER1cKATERmp1artdu+K4oCtVpt6ZKIqMbatUBwMKoBIDhYPidq5jgHmYgaDefqEVmZtWuBiROB8nJdm6srsHIlMGaM5eqyIfy5Zls4B/k21Go1wsPDLV0GERGR5c2bpx+OAfl83jzL1ENkJawqIMfFxWn/CdbJyQmhoaFYtGgRqqurLV0a2bjLly9jzJgx8PDwgJeXF+Lj4xv8vSorK8PEiRPh4+MDFxcXdOvWDadOndL2nzp1CsOGDYOnpydUKhXGjh2L0tJSbf+NGzfwxhtvICgoCM7OzoiKikJqaqrJPiMRUYPl5TWsnaiZsKqADACDBw9GQUEBTp48iVdeeQXz5s3D0qVLLV0W2bjJkyfjwIED+Pnnn7FmzRokJydj8eLFBr9eCIERI0Zg3759WLduHbKysrBkyRK4ubkBAKqqqvD444/Dzs4Ou3btwg8//ICjR4/qrXf83nvvISkpCStWrMCxY8cQExODxx57DOfPnzf55yUiMkhgYMPaiZoLIYTJH1FRUcIY48aNE8OGDdNrGzx4sOjVq5f2+dmzZwUA8d///lc8+uijwsXFRbRp00Zs3bpVCCFEbm6uGDZsmHB1dRUqlUrExcWJ3377Tfv6hIQEERYWJmbNmiU8PDyEn5+fSEpKMqpesg0XL14UdnZ24ttvv9W2vf3228Lf39/gc2zfvl04ODiI7OzsevuzsrIEAHH8+HFt26ZNm4SiKKKoqEgIIUTPnj3FrFmztP3V1dXCy8tLLFmypKEfyWbJHzlEZDXWrBHC1VUIQPdwdZXtZBD+XLMtADKEAVnW6kaQb+Xi4oKbN2/WaZ85cyZGjBiBI0eO4LPPPtOO5I0ePRplZWXYvXs3Nm3ahNTUVMyaNUvvtSdOnEBBQQH279+POXPmYOLEiTh27FijfB5qfAcPHoRGo0H//v21bQMGDMC5c+dQWFho0DlSUlIQFhaG9evXw9/fH6GhoUhISEBVVRUAOX0CgN720c7OzhBC4PDhw9pjavfb2dnB0dERhw4duufPSERklDFj5A15QUHQAEBQEG/QI4IVTrGoodFo8H//93/YsmULBg8eXKf/2WefxcSJE9GxY0cMHToUvXr1wrFjx7Br1y58/PHHiIiIQP/+/bFgwQKsXr1aG2AAGUz+/ve/44EHHsCrr76K6OhorF69ujE/HjWiCxcuwN7eHiqVCj169MBzzz0HLy8vbZ8hCgoKcPbsWWzduhUbN27EO++8gw8//BCJiYkAgNDQULRp0wZLly5FRUUFSkpKsGzZMiiKon2PmJgYfPXVV8jOzkZ1dTU++eQTFBcXG1wDEZFZjBkD5OTAHgBychiOiWCFAXnLli1wd3eHs7MznnzySYwbN67eNVNrjwbWOH36NOzs7PRWqYiIiMCNGzeQV+uGg3bt2kGlUmmfh4WF4cyZM6b9IGRVam7+DAgIgJ+fX4Nfr9FocOXKFSQnJ6Nnz554+umnMXHiRKxZswaAHC1et24dtmzZAjc3N7Rr1w6xsbEQQkBRFAByBZWwsDB07NgRTk5O+OabbzB06FBtPxEREVkHB0sXcKsBAwZg5cqVcHJygp+fH+zt7es9ztPT06Tvy5DSdHl7e6OqqgqlpaXYuHEjACAtLU3bZwgvLy84OzsjICBA29a+fXvk5+drnw8cOBBnz55FcXExXF1dUVBQgNdffx3t2rUDAKhUKmzatAnl5eW4cuUKfH190bdvX4SGhprqoxIREZEJWN0IsqurKzp27IiAgIDbhuPb6dChAzQajd584sOHD6NFixYIrHVHbn5+PsrKyrTPjx8/jo4dO9578WSVIiMjYWdnpw3FAJCamop27drB19dX79iysjLk5OSgoqKizjkqKir0VpzIy8vTht/afHx84O7ujnXr1sHd3R3R0frrkbu6usLX1xe5ubnYu3cvYmNjTfApiYiIyFSsLiDfi65du6J3796YOnUqDh8+jLS0NCQkJODFF1+Eo6Oj9jiNRoPp06fj5MmT+Pvf/46MjAy95bioafHy8sIzzzyD2bNnY//+/fjhhx+wbNkyTJkypc6xiYmJCAkJwd69e/Xahw0bhrZt22Ly5Mn45Zdf8OOPP2LlypUYU2uu3pYtW5Camorc3Fx89tlnePfddzFr1iy4uroCAIqLi/HFF1/gzJkz2Lt3L5555hmEhobiueeeM+8FICIiogaxuikW9+qrr77Cyy+/jN69e8PR0RFPPPEElixZonfMAw88gDZt2iA6OhotW7bEP/7xD3Tp0sVCFVNjWLFiBaZMmYJBgwbByckJL774IubOnWvw693c3LB582ZMnToVkZGR8PHxwaRJkzBjxgztMZcuXcKf//xnFBUVISAgAGq1GjNnztQ7zwcffIBJkybB2dkZQ4cOxbJly/RWtiAiIiLLU4QZ9g+Pjo4WGRkZJj8vEdk2RVFgjp85RHTv+PvTOLxutkVRlANCiOi7HdekplgQEREREd0rBmQiIiIioloYkImIiIiIamFAJiIiIiKq5a4BWVGUUEVRMms9riiK8lpjFEdERERE1NjuusybEOIkgG4AoCiKPYB8AP82c11ERERERBbR0CkWDwM4I4TINUcxRERERESW1tCA/DyAdeYohIiIiIjIGhgckBVFcQTwOICvb9M/UVGUDEVRMi5cuGCq+oiIiIiIGlVDRpAfBXBQCFFUX6cQYqUQIloIEe3t7W2a6oiIiIiIGllDAvIocHoFERERETVxBgVkRVHcAPwBwEbzlkNEREREZFl3XeYNAIQQ1wC0NnMtREREREQWx530iIiIiIhqYUAmIiIiIqqFAZmIiIiIqBYGZCIiIiKiWhiQiYiIiIhqYUAmIiIiIqqFAZmIiIiIzE6tVkNRFO1DrVZbuqTbYkAmIiIiaqi1a4HgYFQDQHCwfE53pFarIYQAAAghrDogG7RRCBERERH9bu1aYOJEoLxcjjTm5srnADBmjCUrIxPhCDIRERFRQ8ybB5SX67eVl8t2ahIYkImIiIgaIi+vYe1kcxiQiYiIiBoiMLBh7WRzGJCJiIiIGmLhQsDVVb/N1VW2U5PAgExERETUEGPGACtXAkFB0ABAUJB8zhv0mgyuYkFERETUUGPGAGPGwF5RIHJyLF0NmRhHkImIiIiIauEIMhERERGZV1kZsHMnkJKCTyxdiwEYkImIiIjItEpLgbQ0ICUF2LEDOHQIEAJwckJ7QP5aUSxc5O0xIBMRERHRvbl0CUhNlWF4xw7g8GEZgp2dgd69gYQEIDYWeOghDHVxgbDicAwwIBMRERFRQ124oAvEKSnA0aOy3cUF6NMHWLBABuKePQEnJ0tWahQGZCIiIiK6s+Ji3ehwSgpw/Lhsd3UF+vYFnntOBuIePQBHR0tWahIMyERERESkr7BQPxD/8otsd3MD+vUDxo4FYmKAqKgmEYhvxYBMRERE1NydP68Lwzt2ACdPynYPDxmI4+JkIO7eHWjRwpKVNgoGZCIiIqLm5tw5/UD866+yvWVLoH9/4KWXZCCOjAQcml9cbH6fmIiIiKi5ycvTD8Rnzsh2lUoG4smTZSDu1g2wt7doqdaAAZmIiIioqcnJ0YXhlBT5HAA8PWUQnjpV/vfBBxmI68GATERERGTLhADOntUPxHl5sq91axmEZ8yQq0yEhwN2dhYs1jYwIBMRERHZEiHkFInagfjcOdnn7S0D8axZMhB36cJAbAQGZCIiIiJrJgRw6pT+HOLz52Wfj48MwjEx8r+dO1v1Fs62ggGZiIiIyJoIAZw4oR+ICwtln6+vfiAODWUgNgMGZCIiIiJLEgLIytIPxMXFss/PDxg0SBeKO3ViIG4EDMhEREREjUmjkVs114Th1FTgwgXZ5+8PPPKILhB36MBAbAEMyERERETmpNEAR4/qB+JLl2RfYCDw6KO6QBwSwkBsBRiQiYiIiEypuho4ckQ/EJeWyr6QEGD4cF0gDg62YKF0OwzIRERERPeiqgrIzJRhuCYQX74s+zp0AJ58UheIAwMtWioZhgGZiIiIqCGqqoCDB4EdO/AdIDfjuHJF9nXqBDz7rAzDMTFyTjHZHAZkIiIioju5eRM4cEC3ysSuXcBvvwEAOgDAqFG6QOznZ8lKyUQYkImIiIhqu3EDyMjQD8TXrsm+zp2BsWO1gbhz27YQy5dbtFwyPQZkIjK/tWuBefNQDcgbUhYuBMaMsXBRRES/q6wE0tN1gXj3bqC8XPaFhwNxcTIQDxgAtGljyUqpkTAgE5F5rV0LTJwIlJfDDgByc+VzgCGZiCyjshLYt0+3ysTu3UBFhex78EFgwgR5U13//oC3tyUrJQtRhBAmP2l0dLTIyMgw+XmJyAYFB8tQfKugICAnp7GrIaLbUBQF5sgEVqGiAti7VxeI9+yRIVlRgIgI3bbN/fvLG+4aoElfNzOx5DVTFOWAECL6bsdxBJmIzCsvr2HtRET3qrxcPxDv3SvnFdvZAd26AS+/rAvEnp6WrpasEAMyEZlXYGD9I8hcC5SITOXaNTkqXBOI9+2TK0/Y2QHduwOvvipHifv1A1QqS1dLNsDO0gUQURO3cCHg6qrf5uoq24nI4tRqNZTftzZWFAVqtdqyBRni6lXgxx+BN94A+vaVofcPfwAWL5YjxTNmAN9/D5SUyJvvliwB/vhHhmMyGOcgE5H5/b6KhSY3F3ZBQVzFgoga5rff5FJrNSPEGRlysw4HByA6WrdLXd++gIdHo5bGOcgNZwtzkBmQiajR8A8SIjLI5cvAzp26ZdcOHgSqq4EWLYCePXU31fXuDbi7W7RU/lxrOFsIyJyDTERERJZVVgakpekC8aFDgEYDODoCDz0EvP66LhDfOmWLyAwYkImIiKhxlZToB+LMTEAIwMkJ6NULmD9fBuJevQAXF0tXS80QAzIRERGZ16VLQGqqLhAfOSIDsbOzHBVOSJCB+KGHZBuRhTEgExERkWlduKAfiI8ele0uLkCfPsBbb8l5xD17ylFjIivDgExERET3prhYhuGaQHz8uGx3dZUrSzz/vAzEPXrIecVEVo4BmYiImgy1Wo0FCxZonyckJNjGur62prBQF4Z37AB++UW2u7vLQDx2rAzE0dFy5QkiG8Nl3ogagH/43hsuh0SNhd81Ezt/Xj8Qnzwp2z085HbNMTHy0b17swvE/K41nC0s88aATGQE/kA0Dq8bNRZ+1+7RuXP6gfjXX2V7y5bAgAG6dYi7dZObdTRj/K41nC0E5Ob9rSYiIiIgL08/EJ85I9tVKhmIJ0+WgTgiArC3t2SlRI2CAZmIiKi5ycnRheGUFPkcAFq1koF46lQZiLt2ZSCmZokBmYiIqCkTAjh7Vj8Q5+XJvtat5XSJGTPkf7t2BezsLFktkVVgQCYiImpKhJBTJGoH4nPnZJ+3twzCs2fL/3bpwkBMVA8GZCIiIlsmhLyJrnYgPn9e9rVpo7uhLiYG6NwZUBQLFktkGwwKyIqiqACsAhAOQAAYL4TYY87CiIiIqB5CACdO6N9UV1go+9q21Q/EoaEMxERGMHQE+QMAW4QQIxVFcQTgasaaiIiIqIYQQFaWfiAuLpZ97doBDz+sC8UdOzIQE5nAXQOyoij3ARgAIA4AhBA3ANwwb1lERETNlEYjt2quCcOpqcCFC7IvIAAYMkQXiNu3ZyAmMgNDRpBDAFwAkKwoSgSAAwCmCyGu1T5IUZSJACYCQGBgoKnrJCIiapo0GuDoUf1AfOmS7AsKAh57TBeIg4MZiIkawV130lMUJRrAXgB9hRD7FEX5AMAVIcRfb/ca7qRHTR13TjIOrxs1Fqv+rlVXA0eO6Afi0lLZFxKimz8cEyMDMVk1q/6uWammspPeOQDnhBD7fn++AcDceymOiIio2aiqAjIzZRiuCcSXL8u+Dh2Ap57SBWL+CyyRVbhrQBZCFCqK8j9FUUKFECcBPAwgy/ylERER2aCqKuDgQd1NdTt3AleuyL5OnYBnn9UFYn9/i5ZKRPUzdBWLaQDW/r6CRTaAF81XEhERkQ25eRM4cEAXiHftAn77TfY98AAwapScNjFgAODnZ8lKichABgVkIUQmgLvO1yAiImrybtwAMjL0A/G13+9b79IFGDtWF4h9fS1ZKREZiTvpERER3UllJZCerrupbvduoLxc9oWHA3FxukDs42PBQonIVBiQiYiIaquoAPbv1w/EFRWy78EHgQkTZCDu3x/w9rZkpURkJgzIRETUvF2/DuzbpwvEe/bIUWNFASIigEmTdIG4dWtLV0tEjYABmYiImpfychmCa+YQ79sn5xXb2QHdugEvv6wLxJ6elq6WqOlYuxaYNw/VgFzje+FCYMwYCxdVPwZkIiIrpFarsWDBAu3zhIQEqNVqyxVky65dk9MkagLx/v1y5Qk7O6B7d+DVV+WSa/36ASqVpaslaprWrgUmTgTKy2EHALm58jlglSH5rjvpGYM76VFTx52TjMPr1nC8Zka4ehVDPDzww+uvy0Ccni7XJra3B6KidDvV9e0L3HefpaslG8ffowYKDpah+FZBQUBOTqOVYehOegzIREbgD0Tj8Lo1HK+ZAa5ckUut1YwQZ2TI7ZwdHIDoaP1A7OFh6WqpieHvUQPZ2QH1XSdFATSaRivDlFtNExERWY/Ll+XudDU31R08KANxixZAjx7A7Nl4ZNEi/FhaCri7W7paIgLkNur1jSBb6fbqDMhERGTdysqAtDRdID50SI44tWgB9OoFvP66HCXu3RtwdQUA/LRoEcMxkTVZuFA7B1nL1VW2WyEGZCIisi4lJfqBODNT/tOso6MMwfPny0Dcqxfg4mLpaonIEDU34s2bB01uLuyCgqx6FQvOQSYyAuecGYfXreGaxTW7dAlITdUF4iNHZCB2dpaBOCZGBuKHHpJtBmgW140siivNGM+Svz95kx6RGfEPX+PwujVck7xmFy7IIFzzOHpUtru4AH366AJxz56Ak5NRb9EkrxtRE2ELAZlTLIiIyLyKinRhOCUFyMqS7a6ucmWJ556TgbhHDzmNgojIwhiQiYjItAoK9APxiROy3c1NbsYxdqwMxFFRDMREZJUYkImI6N7k5+sH4lOnZLuHhwzEL74op0107y5XniAisnIMyERE1DD/+59+ID59Wra3bAn07w/8+c8yEEdGys06iIhsDH9yERHRneXm6gfi7GzZrlLJQDxligzE3brJ7ZyJiGwcAzIREenLyZFBuGbZtZwc2e7pCQwYAEybJgPxgw8yEBNRk8SATETUnAkhR4RrRod37ADy8mRf69YyEM+YIQNx166AnZ1FyyUiagwMyEREzYkQcs5w7UB87pzs8/KSQXjWLPnfsDAGYiJqlhiQiYiaMiHkqhK1A/H587LPx0cut1azMUfnzoCiWLBYIiLrwIBMRNSUCCHXHa4diAsLZZ+vr34gDg1lICYiqgcDMhGRLRNC7kxXOxAXF8s+Pz9g0CBdIO7UiYGYiMgADMhERLZEowGOH9eF4dRU4MIF2efvDzzyiC4Qd+jAQExEZAQGZCIiK6YAwOHD+oH40iXZGRgIPPqoLhCHhDAQExGZAAMyEZE1qa4GjhzRBuKLgNyAAwCCg4Hhw3WBODjYUlUSETVpDMhERJZUXQ1kZuqPEF++LPs6dMC/AUz4/HMZioOCLFkpEVGzwYBMRNSYqqqAgwd1WzenpQFXrsi+Tp2AZ57RrTTh74+XFAUTXnjBoiUTETU3DMhEROZ08yZw4IBulYldu4DffpN9oaHA88/rArGfnyUrJSKi3zEgExGZ0o0bQEaGfiC+dk32de4MjB0rw/CAAUDbthYtlYiI6seATER0LyorgfR0XSDevRsoL5d9YWFAXJwuELdpY8lKiYjIQAzIREQNUVEB7N+vu6lu927ZBgBduwITJugCsbe3RUslIiLjMCATEd1JRQWwd68uEO/ZI0eNFQWIiAAmTdIF4tatLV0tERGZAAMyEVFt5eX6gXjvXjmvWFGAyEjg5ZflTXX9+wOenpauloiIzIABmYiat2vX5DSJmjnE+/fLlSfs7IDu3YFp02Qg7tcPUKksXS0RETWCJhOQ1Wo1FixYoH2ekJAAtVptuYKIyDpdvSpXlqgJxOnpcm1ie3sgKgp47TUZiPv2Be67z9LVEhGRBShCCJOfNDo6WmRkZJj8vIZQFAXm+ExEtfF7ZhyLXLcrV/QDcUaG3L3OwQGIjtatQdy3L+Dh0bi1GYDfNePwuhFZL0v+/lQU5YAQIvpuxzWZEWQiIgBym+adO3VziA8elIG4RQugRw9gzhwZiPv0AdzdLV0tERFZIQZkIrJtZWVyu+aaQHzoEKDRyEDcqxfw+utylLh3b8DV1dLVEhGRDWBAJiLbUlKiH4gzMwEhACcnGYjnz5eBuFcvwMXF0tUSEZENsrN0AUTU9KnVaiiKAkDOPWvQDbQXLwIbNwLTpwPdugFeXsCIEcDy5XJViYQEGZbLyuR/FywABg5kOCYiIqPxJj0iI/B7ZkYXLsiR4Zqb6o4dk+0uLnLecEyMHCHu2VOOGjdx/K4Zh9eNyHrxJj0iorspKtIPxFlZst3VVa4sMWqUDMU9egCOjhYtlYiImgcGZCJqXAUF+oH4xAnZ7u4uA/Gf/iQDcXS0vNGOiIiokTEgE5F55efrB+JTp2S7h4fcne7FF2Ug7t6dgZiIiKwCAzIRmdb//qcLwzt2AKdPy/aWLYH+/YE//1kG4shIuVkHERGRleGfTkR0b3Jz9QNxdrZsV6mAAQOAKVPkTXUREXI7ZyIiIivHgExEhhMCyMnRD8Q5ObLP01OODE+bJgNx164MxEREZJMYkIno9oSQI8K1A3Fenuxr3VoG4hkzZCAODwfsuLQ6WdjatcC8eagGgOBgYOFCYMwYCxdFRLaGAZmoIZr6H75CyDnDtQPxuXOyz9tbBuJZs2Qg7tKFgZisy9q1wMSJQHm53AUrN1c+B5rW71MiMjtuFEJkqFp/+Gq5ugIrV9ruH75CyFUlasJwSopchg0A2rTRbcoREwN07gz8vhseNR7+TGuA4GAZim8VFKSbCkREFmcLG4UwIBMZqin84SuEXHe4JhDv2AEUFsq+tm31A3FoKAOxFeDPtAaws5Pf8VspCqDRNH49RFQvWwjInGJBZKiaubeGtlsDIeTOdLUDcXGx7GvXDhg0SBeIO3ViICbbFhhY/19iAwMbvxYismkMyESGsoU/fDUa4Phx/UB88aLsCwgAHnlEF4g7dGAgpqZl4cL6p0EtXGi5mojIJjEgExnKGv/w1WiAo0f1A3FJiewLCgKGDdNNmwgOZiCmpq3mXoB586DJzYVdUFDTu5GWiBoF5yATNcTvq1hY7A/f6mrg8GFdGE5NBUpLZV9IiP4c4uDgxquLzIY/04zD60ZkvWxhDjIDMpERGu17VlUFZGbqVphISwMuX5Z9HTrownBMjHVN9SCT4c804/C6EVkvWwjInGJBZE2qqoCDB3WBeOdO4MoV2depE/Dss7pA7O9v0VKJiIiaKgbkZkytVmPBggXa5wkJCVCr1ZYrqDm6eRM4cEA/EF+9KvtCQ4FRo3SB2M/PoqUSERE1FwZNsVAUJQfAbwCqAVTdbWiaUyxsC69Zwxl9zW7cADIydIF41y7g2jXZ17mz/pQJX19Tlkw2ir8/jcPrRmS9mtoUi4FCiIv3UBNR81NZCaSn61aZ2L1btwpGeDgQFyfD8IABcuc6IiIisjhOsSAypcpKYN8+/UBcUSH7HnwQmDBBjhL37w94e1uyUiIiIroNQwOyAPCjoigCwAohxMpbD1AUZSKAiQAQyLvpqbmoqAD27tUF4j17ZEhWFCAiApg0SReIW7e2dLVERERkAEMDcj8hRL6iKD4AflIU5YQQIrX2Ab+H5pWAnINs4jqJrEN5ObB3LxYAcmrE3r1yXrGiAJGRwMsv6wKxp6eFiyUiIiJjGBSQhRD5v/+3WFGUfwPoCSD1zq8iagKuXZOjwjUjxPv2ATdvYh4gw/K0aTIQ9+sHqFSWrZWIiIhM4q4BWVEUNwB2Qojffv/1IwDeMntlRJZw9apcWaJmlYn0dLk2sb09EBUFvPYaEBuLVsOG4XJ6uqWrJSIiIjMwZAS5DYB/K4pSc/xXQogtZq2KqLH89ptce7gmEGdkyO2cHRyA6Ghg5kw5laJvX8DDQ/uyK5armIiIiMzsrgFZCJENIKIRaiEyv8uX9QPxwYMyELdoAfToAcyZIwNxnz6Au7ulqyUiIiIL4DJv1LSVlQFpabo5xIcOARqNDMQPPQS8/roMxL17A25ulq6WiIiIrAADMjUtJSX6gTgzExACcHICevUC5s+XN9X16gW4uFi6WqLbW7sWmDcP1QAQHAwsXAiMGWPhooiImgcGZLJtly4Bqam6QHzkiAzEzs5yVDghQQbihx6SbUS2YO1aYOJEoLwcdgCQmyufAwzJRESNQDHHXtjR0dEiIyPD5Oc1hCX397ZVNnXNLlyQQbjmcfSobHdxkfOGY2JkIO7ZU44am4lNXTOyPcHBMhTfKigIyMlp7GpsEn+PElkvS/7+VBTlgBAi+m7HcQSZrFtRkS4Mp6QAWVmy3dVVrizx/PMyFPfoATg6WrRUIpPJy2tYOxERmRQDMlmXggL9QHzihGx3c5ObcfzpTzIQR0fLG+2ImqLAwPpHkNA0CeUAAA7JSURBVAMDG78WIqJmyM7SBVAzl58PfPUVMGkSEBoK+PkBo0bJOZghIcC778rtnEtLgS1bgLlz5dxihmNqyhYulP9KUpurq2wnIrJRarUav++rAUVRoFarLVvQHXAOMjXuNfvf//RHiE+flu0tWwIDBsjR4ZgYIDJSbtZhpfg9I7P7fRULTW4u7IKCuIpFA/H3KBHVx9A5yAzIZN5rlpurC8M7dgDZ2bJdpdIF4thYICJCbudsI/g9o8bC75pxeN2IqD68SY8anxDyDvvagbjmjntPTxmGp02TgbhrV5sKxERERNR8MCCT8YSQI8K1A3HNXfatW8tAPGOGDMTh4YAdp7wTERGR9WNAJsMJIecM14ThlBR5kx0AeHvLQDxrlgzEXbowEBMREZFNYkCm2xMCOHVKPxAXFMi+Nm1084djYoDOnYHf70wlIiIismUMyKQjhFx3uCYQ79gBFBbKvrZt9QNxaCgDMRERETVJDMjNmUYDZGXhZQB49lkZiIuLZV+7dsCgQbpA3KkTAzERERE1CwzIzYlGAxw7ppsukZoKXLyIjwG5GceQIbpR4vbtGYiJiIioWWJAbso0GuDIEf1AXFIi+4KCgGHDgJgYtB8/Htm5uQzERERERGBAblqqq4HDh3WBOC1NbtEMyG2bn3hCt1NdcLD2ZWfHj2c4JiIiIvodA7Itq6oCMjN1N9WlpQGXL8u+Dh2Ap57SBeLAQIuWSkRERGQrGJBtSVUVcPCgfiD+7TfZd//98ka7mpvq2rWzZKVERERENosB2ZrdvAkcOKALxDt3Alevyr4HHgBGj5aBeMAAwM/PkpUSERERNRkMyNbkxg0gI0MXiHftAq5dk31dugB/+pMuEPv6WrJSIiIioiaLAdmSKiuB/ft1m3Ls2gVcvy77wsOBuDhdIPbxsWSlRERERM0GA3JjqqgA9u3TrTKxZ49sA4AHHwReekkXiL28LFkpERERUbPFgGxO16/LDThqAvHevXLUWFGAiAhg0iQZiPv3B1q3tnS1RERERAQGZNMqL5ejwjWBeN8+Oa/Yzg7o1g145RW5wkT//oCnp6WrJSIiIqJ6MCDfi2vXgN27dYF4/3658oSdHdC9O/DqqzIQ9+sHqFSWrpaIiIiIDPD/7d1fjKVlfQfw749dEAZBG902VGS3V6aNF6LHBaUlDdSqldCmVzYjF73ZNjEN0hitYMLxgosmTeNdkw20pXGVUNQbtAYTsdQLqbNA/Ydp+mfZYm2hobDdblMVfr04h/ra7jIzy8y+58x8PsnJznnz7rvfeTKz53ue87zvqyBvxsmTsxPpXizEX/va7NrEe/Ykb3lLcsstPyrEl146dloAAM6CgvxSTpyYFeIXL7u2tja7nfPevclkknzwg7NCfM01ySWXjJ0WAIAtoCAPPffc7GYcLxbio0eTF15Izj8/OXgw+fCHZyfVve1tyStfOXZaAAC2we4uyM8+O7td84uF+NFHZ4X4gguSq65Kbr31R4V4ZWXstAAAnAO7qyA/88yPF+LHHku6k1e8Irn66uSjH50V4quvTi66aOy0AACMYOcU5CNHkttuy/NJcuBAcscdyTvfmTz00I/uVPf1r88K8YUXzmaFb799Voivumq2DQCAXW9nFOQjR5JDh5JTp3JekjzxRHLTTbMynMxmg9/+9uRjH5sV4oMHZ7PGAADwf+yMgnzbbbObdAx1J696VfK5zyVvfetsXTEAAKxjZxTk48dPv/3Eidkl2AAAYIPOGzvAlrjiis1tBwCAM9gZBfmOO/7/ZdhWVmbbAQBgE3ZGQV5dTQ4fTvbvzwtJsn//7Pnq6tjJAABYMjtjDXIyK8Orq9lTlT52bOw0AAAsqZ0xgwwAAFtEQQYAgAEFGQAABhRkAAAYUJABAGBAQQYAgAEFGQAABhRkAAAYUJABAGBAQQYAgAEFGQAABhRkAHaM6XSaqkqSVFWm0+m4gYClpCDDJnjxhcU2nU7T3f/78DsKnI3q7i0/6GQy6bW1tS0/7kZUVbbje9rJjBksLr+fAFunqo5292S9/cwgAwDAgIK8mx05khw4kOeT5MCB2XMAgF1u79gBGMmRI8mhQ8mpU7N3SU88MXueJKurYyYDABiVGeTd6rbbklOnfnzbqVOz7QAAu5iCvFsdP7657QAAu8SGC3JV7amqR6vq/u0MxDlyxRWb2w4AsEtsZgb55iSPb1cQzrE77khWVn5828rKbDsAwC62oYJcVZcneU+SO7c3DufM6mpy+HCyf39eSJL9+2fPnaAHAOxyG72KxceTfCjJJWfaoaoOJTmUJFf4mH45rK4mq6vZU5U+dmzsNAAAC2HdGeSquiHJU9199KX26+7D3T3p7sm+ffu2LCAAAJxLG1licU2SG6vqWJJ7klxXVZ/Y1lQAADCSdQtyd3+kuy/v7gNJ3pvkS939vm1PBgAAI3AdZAAAGNjUraa7+8tJvrwtSQAAYAGYQQYAgAEFGQAABhRkAAAYUJABFtB0Ok1VJUmqKtPpdNxAALtIdfeWH3QymfTa2tqWH3cjqirb8T3tZMYMANgNqupod0/W288MMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAwsGMK8nQ6TVUlSaoq0+l03EAAACyl6u4tP+hkMum1tbUtPy7bo6qyHT8HAACLpKqOdvdkvf12zAwyAABsBQUZAAAGFGQAABhQkAEAYEBBBgCAAQUZAAAGFGQAABhQkAEAYEBBBgCAAQUZAAAGFGQAABhQkAEAYEBBBgCAAQUZAAAGFGQAABhQkAEAYEBBBgCAAQUZAAAGFGQAABhQkAEAYEBBBgCAAQV5F5tOp6mqJElVZTqdjhsIAGABVHe/9A5VFyZ5KMkrkuxNcl933/5Sf2cymfTa2tqWhQQAgJerqo5292S9/fZu4Fj/neS67j5ZVecn+UpV/UV3f/VlpwQAgAWzbkHu2RTzyfnT8+ePl552BgCAJbWhNchVtaeqHkvyVJIvdvfDp9nnUFWtVdXa008/vdU5AQDgnNhQQe7u57v7TUkuT3Kwqt54mn0Od/ekuyf79u3b6pwAAHBObOoqFt39bJIHk7xre+IAAMC41i3IVbWvql49//qiJO9I8p3tDgYAAGPYyFUsLktyd1XtyaxQ39vd929vLAAAGMdGrmLx9SRXnoMsAAAwOnfSAwCAAQUZAAAGFGQAABhQkAEAYEBBBgCAAQUZAAAGFGQAABhQkAEAYKC6e+sPWvV0kie2/MAb89ok/zbSv72sjNnmGbOzY9w2z5idHeO2ecbs7Bi3zRtzzPZ39771dtqWgjymqlrr7snYOZaJMds8Y3Z2jNvmGbOzY9w2z5idHeO2ecswZpZYAADAgIIMAAADO7EgHx47wBIyZptnzM6Ocds8Y3Z2jNvmGbOzY9w2b+HHbMetQQYAgJdjJ84gAwDAWVOQAQBgYEcU5Kr646p6qqq+OXaWZVJVr6+qB6vq21X1raq6eexMi66qLqyqv66qv5mP2cfGzrQsqmpPVT1aVfePnWVZVNWxqvpGVT1WVWtj51kGVfXqqrqvqr5TVY9X1dvGzrToquoN85+xFx8nquoDY+dadFV1y/x14JtV9amqunDsTMugqm6ej9m3FvnnbEesQa6qa5OcTPJn3f3GsfMsi6q6LMll3f1IVV2S5GiSX+vub48cbWFVVSW5uLtPVtX5Sb6S5Obu/urI0RZeVf1ukkmSS7v7hrHzLIOqOpZk0t1uQrBBVXV3kr/q7jur6oIkK9397Ni5lkVV7Uny3SRXdfdYN/xaeFX1usz+//+57v6vqro3yee7+0/HTbbYquqNSe5JcjDJ95N8Iclvd/ffjRrsNHbEDHJ3P5TkmbFzLJvu/l53PzL/+j+SPJ7kdeOmWmw9c3L+9Pz5Y/nfZW6zqro8yXuS3Dl2FnauqnpVkmuT3JUk3f195XjTrk/y98rxhuxNclFV7U2ykuSfR86zDH42ycPdfaq7f5jkL5P8+siZTmtHFGRevqo6kOTKJA+Pm2TxzZcKPJbkqSRf7G5jtr6PJ/lQkhfGDrJkOskDVXW0qg6NHWYJ/EySp5P8yXw5z51VdfHYoZbMe5N8auwQi667v5vkD5IcT/K9JM919wPjploK30zyC1X1mqpaSfIrSV4/cqbTUpBJVb0yyaeTfKC7T4ydZ9F19/Pd/aYklyc5OP/IiDOoqhuSPNXdR8fOsoR+vrvfnOTdSd4/X07Gme1N8uYkf9TdVyb5zyS/N26k5TFfknJjkj8fO8uiq6qfSPKrmb0p++kkF1fV+8ZNtfi6+/Ekv5/kgcyWVzyW5PlRQ52BgrzLzdfRfjrJke7+zNh5lsn8o9sHk7xr7CwL7pokN87X096T5Lqq+sS4kZbDfJYq3f1Uks9mtm6PM3syyZODT3Xuy6wwszHvTvJId//r2EGWwC8l+cfufrq7f5DkM0nePnKmpdDdd3X3W7r72iT/nuRvx850OgryLjY/4eyuJI939x+OnWcZVNW+qnr1/OuLkrwjyXfGTbXYuvsj3X15dx/I7OPbL3W3mZZ1VNXF85NnM18m8MuZfTzJGXT3vyT5p6p6w3zT9UmcdLxxvxHLKzbqeJKrq2pl/lp6fWbn8bCOqvrJ+Z9XZLb++JPjJjq9vWMH2ApV9akkv5jktVX1ZJLbu/uucVMthWuS3JTkG/M1tUlya3d/fsRMi+6yJHfPz/Q+L8m93e2yZWyHn0ry2dlrb/Ym+WR3f2HcSEvhd5IcmS8X+IckvzlynqUwfxP2jiS/NXaWZdDdD1fVfUkeSfLDJI9mCW6fvCA+XVWvSfKDJO9f1BNpd8Rl3gAAYKtYYgEAAAMKMgAADCjIAAAwoCADAMCAggwAAAMKMgAADCjIAAAw8D/qW/TguOlJDAAAAABJRU5ErkJggg==\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.6.6"
  },
  "main_language": "python"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}