{ "cells": [ { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "iminuit version: 2.8.2\n" ] } ], "source": [ "# Example of maximum-likelihood fit with iminuit version 2.x\n", "# pdf is a mixture of Gaussian (signal) and exponential (background),\n", "# truncated in [xMin,xMax].\n", "# G. Cowan / RHUL Physics / August 2021\n", "\n", "import numpy as np\n", "import scipy.stats as stats\n", "from scipy.stats import truncexpon\n", "from scipy.stats import truncnorm\n", "from scipy.stats import chi2\n", "import iminuit\n", "from iminuit import Minuit\n", "import matplotlib.pyplot as plt\n", "plt.rcParams[\"font.size\"] = 14\n", "print(f\"iminuit version: {iminuit.__version__}\") # should be v 2.x" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# Select fit type, define pdf\n", "fitType = 'M' # choose least squares (LS) or multinomial ML (M)\n", "np.random.seed(seed=1234567) # fix random seed\n", "theta = 0.2 # fraction of signal\n", "mu = 10. # mean of Gaussian\n", "sigma = 2. # std. dev. of Gaussian\n", "xi = 5. # mean of exponential\n", "xMin = 0.\n", "xMax = 20." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def f(x, par): # fit function\n", " theta = par[0]\n", " mu = par[1]\n", " sigma = par[2]\n", " xi = par[3]\n", " fs = stats.truncnorm.pdf(x, a=(xMin-mu)/sigma, b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n", " fb = stats.truncexpon.pdf(x, b=(xMax-xMin)/xi, loc=xMin, scale=xi)\n", " return theta*fs + (1-theta)*fb" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "class ChiSquared: # function to be minimized\n", "\n", " def __init__(self, xHist, bin_edges, fitType):\n", " self.setData(xHist, bin_edges)\n", " self.fitType = fitType\n", " \n", " def setData(self, xHist, bin_edges):\n", " numVal = np.sum(xHist)\n", " numBins = len(xHist)\n", " binSize = bin_edges[1] - bin_edges[0]\n", " self.data = xHist, bin_edges, numVal, numBins, binSize\n", "\n", " def chi2LS(self, par): # least squares\n", " xHist, bin_edges, numVal, numBins, binSize = self.data\n", " xMid = bin_edges[:numBins] + 0.5*binSize\n", " binProb = f(xMid, par)*binSize\n", " nu = numVal*binProb\n", " sigma = np.sqrt(nu)\n", " z = (xHist - nu)/sigma\n", " return np.sum(z**2)\n", " \n", " def chi2M(self, par): # multinomial maximum likelihood\n", " xHist, bin_edges, numVal, numBins, binSize = self.data\n", " xMid = bin_edges[:numBins] + 0.5*binSize\n", " binProb = f(xMid, par)*binSize\n", " nu = numVal*binProb\n", " lnL = 0.\n", " for i in range(len(xHist)):\n", " if xHist[i] > 0.:\n", " lnL += xHist[i]*np.log(nu[i]/xHist[i])\n", " return -2.*lnL\n", "\n", " def __call__(self, par):\n", " if self.fitType == 'LS':\n", " return self.chi2LS(par)\n", " elif self.fitType == 'M':\n", " return self.chi2M(par)\n", " else:\n", " print(\"fitType not defined\")\n", " return -1" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# Generate data\n", "numVal = 200\n", "xData = np.empty([numVal])\n", "for i in range (numVal):\n", " r = np.random.uniform();\n", " if r < theta:\n", " xData[i] = stats.truncnorm.rvs(a=(xMin-mu)/sigma, \n", " b=(xMax-mu)/sigma, loc=mu, scale=sigma)\n", " else:\n", " xData[i] = stats.truncexpon.rvs(b=(xMax-xMin)/xi, loc=xMin, scale=xi)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# Put data values into a histogram\n", "numBins=40\n", "xHist, bin_edges = np.histogram(xData, bins=numBins, range=(xMin, xMax))\n", "binSize = bin_edges[1] - bin_edges[0]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# Initialize Minuit and set up fit:\n", "parin = np.array([theta, mu, sigma, xi]) # initial values (here = true)\n", "parname = ['theta', 'mu', 'sigma', 'xi']\n", "parstep = np.array([0.1, 1., 1., 1.]) # initial setp sizes\n", "parfix = [False, True, True, False] # change to fix/free param.\n", "parlim = [(0.,1), (None, None), (0., None), (0., None)]\n", "chisq = ChiSquared(xHist, bin_edges, fitType)\n", "m = Minuit(chisq, parin, name=parname)\n", "m.errors = parstep\n", "m.fixed = parfix\n", "m.limits = parlim\n", "m.errordef = 1.0 # errors from chi2 = chi2min + 1" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "par index, name, estimate, standard deviation:\n", " 0 theta = 0.201533 +/- 0.052956\n", " 3 xi = 5.148508 +/- 0.651156\n", "\n", "free par indices, covariance, correlation coeff.:\n", "0 0 0.002821 1.000000\n", "0 3 -0.018494 -0.534693\n", "3 0 -0.018494 -0.534693\n", "3 3 0.424135 1.000000\n" ] } ], "source": [ "# do the fit, get errors, extract results\n", "m.migrad() # minimize -logL\n", "parhat = m.values # max-likelihood estimates\n", "sigma_parhat = m.errors # standard deviations\n", "cov = m.covariance # covariance matrix\n", "rho = m.covariance.correlation() # correlation coeffs.\n", " \n", "print(r\"par index, name, estimate, standard deviation:\")\n", "for i in range(m.npar):\n", " if not m.fixed[i]:\n", " print(\"{:4d}\".format(i), \"{:<10s}\".format(m.parameters[i]), \" = \",\n", " \"{:.6f}\".format(parhat[i]), \" +/- \", \"{:.6f}\".format(sigma_parhat[i]))\n", " \n", "print()\n", "print(r\"free par indices, covariance, correlation coeff.:\")\n", "for i in range(m.npar):\n", " if not m.fixed[i]:\n", " for j in range(m.npar):\n", " if not m.fixed[j]:\n", " print(i, j, \"{:.6f}\".format(cov[i,j]), \"{:.6f}\".format(rho[i,j]))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "chi2min = 35.26955011049652 , ndof = 37\n", "pval = 0.5503255503742228\n" ] } ], "source": [ "# Retrieve minimized chi-squared, etc.\n", "ndof = numBins - m.nfit - 1 # for fixed numVal\n", "chi2min = chisq(parhat)\n", "print (\"chi2min = \", chi2min, \", ndof = \", ndof)\n", "pval = chi2.sf(chi2min, ndof)\n", "print (\"pval = \", pval)" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAGnCAYAAABB+tgiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAACUqElEQVR4nOzddVxV5x/A8c+hLo0gtgJ2dyc2xnR293TTtXNzM2Zv+nPTlVNnO7trdqFize4O7EARpOv8/jjCRBDhci+X+L5fr/OSe85znvM9F4Tvfc4TiqqqKkIIIYQQ4p3MTB2AEEIIIURGIYmTEEIIIUQySeIkhBBCCJFMkjgJIYQQQiSTJE5CCCGEEMkkiZMQQgghRDJJ4iSEEEIIkUwWpg4gM4uJieHBgwc4ODigKIqpwxFCCCGyBFVVefnyJXnz5sXMzLBtRJI4GdGDBw8oUKCAqcMQQgghsqS7d++SP39+g9YpiZMROTg4ANo3ztHR0cTRCCGEEFlDYGAgBQoUiPs7bEiSOBlR7OM5R0dHSZyEEEKINGaMbjLSOVwIIYQQIpkkcRJCCCGESCZJnIQQQgghkkkSJyGEEEKIZJLESQghhBAimSRxEkIIIYRIJkmchBBCCCGSKd0lTvfv3+fXX3+ladOmuLm5YWVlRe7cuWnfvj1Hjx5N9JzAwEC++uor3N3d0el0uLu789VXXxEYGJji6x87dowWLVrg7OyMnZ0d1apVY+nSpam9LSGEEEJkAoqqqqqpg3jdd999x//+9z8KFy6Mp6cnOXPm5Nq1a6xfvx5VVVm2bBmdOnWKKx8cHEydOnU4ffo0TZo0oVKlSpw5c4Zt27ZRoUIFfHx8sLOzS9a1vb298fLywsrKii5duuDk5MTatWu5desWP/zwA8OHD0/RvQQGBuLk5ERAQIBMgCmEEEKkEaP+/VXTmTVr1qj79+9PsH///v2qpaWl6uLiooaFhcXtHzVqlAqoQ4cOjVc+dv+oUaOSdd3IyEi1cOHCqk6nU0+ePBm3PzAwUC1durRqYWGhXr16NUX3EhAQoAJqQEBAis4TQgghhP6M+fc33bU4JcXLy4sdO3Zw7NgxqlSpgqqq5M+fn8DAQB49ehSvZSksLIy8efNia2vL3bt33znt+o4dO/Dy8qJv377Mmzcv3rEVK1bQpUsXhg0bxo8//pjseKXFSQghhEh7xvz7m+76OCXF0tISAAsLbYm9a9eu8eDBA2rXrp3gcZy1tTX16tXj/v37XL9+/Z11e3t7A9C0adMEx2L37du3L8k6wsPDCQwMjLcJIYQQIvPIMInTnTt32LVrF7lz56Zs2bKAljgBFC1aNNFzYvfHlktKUnU5Ozvj6ur6znomTpyIk5NT3FagQIF3XlcIIYQQGUeGSJwiIyPp2bMn4eHhTJ48GXNzcwACAgIAcHJySvS82Oa52HJJSU5d76pn2LBhBAQExG13795953WFEEIIkXFYmDqAd4mJiaFfv37s37+fAQMG0LNnT1OH9FY6nQ6dTmfqMIQQQghhJOm6xUlVVQYMGMDixYvp0aMHM2fOjHc8tnXobS1BsX2M3taKlNK6klOPEEIIITKvdJs4xcTE8MEHHzBv3jy6du3KggULMDOLH+67+jC9qw9Ucuvy9/fHz88vWfUIIYQQIvNKl4lTTEwM/fv3Z/78+XTu3JlFixbF9Wt6XdGiRcmbNy8HDx4kODg43rGwsDD2799P3rx5KVKkyDuv6enpCWjTErwpdl9sGSGEEEJkTekucYptaZo/fz4dO3Zk8eLFiSZNAIqi0L9/f4KCghg3bly8YxMnTsTf35/+/fvHm8MpMjKSy5cvc+PGjXjlGzVqRKFChVi6dCmnT5+O2//y5UvGjx+PhYUFffr0Mdh9CiGEECLjSXcTYI4ZM4axY8dib2/PF198ETdn0+vatGlDhQoVgIRLrlSuXJkzZ86wdevWRJdcuX37NgULFsTd3Z3bt2/Hq3fv3r14eXmh0+no2rUrjo6OcUuuTJgwgREjRqToXmQCTCGEECLtGfPvb7obVRebzAQFBfHDDz8kWsbDwyMucbKzs8Pb25uxY8eyevVqvL29yZ07N4MHD2b06NHJXqcOoEGDBvj4+DB69GhWrlxJREQEpUuXZvz48XTv3j21tyaEEEKIjM7gi7iIOLJWnRAiIwBUT09PU4eRbvXu3VsF1Fu3bsXt27t3rwqoo0ePTpPr3bp1SwXU3r17xyvr6empZpQ/5fPnz1cBdf78+Ua/ljH//qa7Pk5CCJGZ3L59G0VRUBSFfPnyER0dnWi5c+fOxZUrUaJEGkeZ9SxYsABFUZg0aZKpQxEZTLp7VCeEEJmRhYUFDx48YPv27bRo0SLB8blz52JhYUFUVFSax3bp0iVsbW3T/LoZxcSJE/nuu+/Ily+fyWLIly8fly5dkvkE0wFpcRJCiDRQq1YtnJycmDdvXoJjERERLFmyJNGEKi2UKFECNzc3k1w7I8iTJw8lSpSIW2jeFCwtLSlRogR58uQxWQxCI4mTEEKkARsbGzp37symTZvw8/OLd2zjxo34+fnRt2/fRM998OABo0ePpkaNGuTMmROdToeHhwcff/wxT548iVf2ypUr2Nvb4+bmhr+/f7xjsS1LHh4e8VZJUBSF+vXrxyvbp08fFEXh5s2b/PzzzxQrVgwbGxtKlSrF8uXLAW16l1GjRlGwYEGsra0pV64c27dvTxC/h4cHHh4eid5b/fr1400ZA9roakVR8Pb2Zv78+ZQtWxYbGxsKFizI77//DmgrS/z222+UKFECa2trihUrxqJFixK9RmrFvhdvjsROzIsXL6hbty7m5ubMmjUrbv/Lly8ZPXo0pUuXxsbGhmzZstGsWTN8fHySFUPsI9+3TYsTFRXF+PHjKViwIDqdjmLFijF9+vREy4aEhDBmzJi4987FxYWWLVty6NAhg5R//vw5AwcOJFeuXNja2lK1alXWrVuXrPvMCORRnRBCpJF+/foxa9YslixZwhdffBG3f968eeTMmZP33nsv0fP279/PlClTaNSoEdWrV8fS0pJTp04xY8YMtm/fzsmTJ+Me4RQvXpxff/2VAQMGMGDAAFavXg1AeHg4Xbt2jWvdSu4jn6+++oqjR4/SqlUrzM3NWb58Od26dcPZ2Zk///yT8+fP06JFC8LCwli6dCmtW7fm8uXLFCxYMJXvFvz66694e3vz/vvv07BhQ9asWcMXX3yBra0tZ86cYdWqVbz33ns0bNiQ5cuX06tXLwoWLEidOnVSfW19PHjwAC8vL65du8aqVato164doCUS9erV48KFC9StWxcvLy8CAgLYsGEDDRo0YNWqVbRp0yZV1+7atStHjx6lefPmmJubs3LlSj755BMsLS0ZMGBAXLnw8HAaNWrEkSNHqFSpEl9++SVPnjxhxYoV7NixgxUrVsTFrU/5kJAQ6tevz7lz56hZsyaenp7cvXuXzp0707Rp01TdY7ph8O7mIk5sr/6HT5+pweGRydpiYmJMHbYQwoBiR0N5eXmpqqqqpUuXVsuVKxd3/N69e6q5ubk6ZMgQVVW1EW7FixePV8fjx4/Vly9fJqh74cKFKqBOmDAhwbEOHTqogDpr1ixVVVX1yy+/fOsoMBIZVRc7sqto0aLqkydP4vYfOXJEBdRs2bKpderUUYOCguKOrVixQgXUzz//PF5d7u7uqru7eyLvTuKjwkaPHq0CqouLi3rjxo24/Xfu3FGtrKxUJycntVixYvHiOnr0qAqorVu3TvQ6b4od4TVx4sR3lk3OqLorV66o7u7uqqOjo7p3795453fr1k0F1Hnz5sXb/+jRI7VAgQJqjhw51NDQ0CSv965RddWrV483guzy5cuqhYVFgp+lcePGqYDavXv3eH9vzpw5o+p0OtXZ2VkNDAzUu3zs927AgAHxrrt9+3YVkFF1Inmq/bCbUqO2J2vrOPMwavqak1QIYUB9+/bl7NmznDhxAtBGd0VHR9OvX7+3npMzZ07s7e0T7O/ZsyeOjo7s2rUrwbHZs2dToEABvvzyS37//Xd+++03atWqxffff5+ieEeMGEGOHDniXlevXp1ChQrx4sULfvjhh3hz5bVv3x5LS0vOnDmTomu8zeeff06hQoXiXhcoUIA6deoQEBCQIK5q1apRqFAhg107JY4dO0bt2rUJCwtj37598R57+vn5sWLFCho1apTgUWyuXLn45ptvePr0aaLfw5SYOHFivIkeixcvTu3atbly5QovX76M279gwQIsLS2ZNGlSvEek5cqVo0+fPvj7+7Nhwwa9y//9999YWVklWM2jadOmNGrUKFX3mF5I4pTOHPf1JzQy8eHKQoiMr2fPnlhaWsZ1El+wYAHVq1enVKlSSZ63du1avLy8yJEjBxYWFiiKgpmZGYGBgTx48CBB+WzZsrFkyRLCw8P54osvcHR0ZMmSJW9dwuptKlasmGBfbAfl2ImIY5mbm5MzZ07u37+fomsY4tqxxwx17eQ6cOAADRs2xMnJiYMHDyaI69ixY0RHRxMWFsaYMWMSbEeOHAHg8uXLqYqjUqVKCfblz58f0PpdgTab9s2bNylSpEjcsdfFJnyxy46ltPzLly+5desWRYoUIXfu3AnK161bN4V3lT5JH6c08O+IRu+c8j0kIpoqE1L3iUMIkf7lzJmTFi1asGzZMlq3bs3169f5+uuvkzxnypQpfP311+TIkYOmTZuSP39+bGxsAK0fUHh4eKLnValShfz58+Pr60vLli3f2kE7KYn97opdCuttxyIjI1N8HUNdO62nczh16hRBQUE0b9480ff3+fPnABw8eJCDBw++tZ43F6pPqcT6rMW+V7FzhwUGBgJaS1diYpOd2IEDKS0f+2/OnDkTLf+2ejIaSZzSgK2VBbZW8lYLITT9+vVjw4YNfPDBB9jY2NC1a9e3lo0dLZU3b15Onz4d7/GUqqpMnjz5recOGTIEX19fsmfPzrJly+jdu7dJOuiamZkRERGR6LHXR/dlRJ9++in3799n3rx5WFhYsGjRoniterEJ3pAhQ/j5559NFWa8WB4/fpzo8dj9seX0Lf/mSM83y2d08qhOCCHSWIsWLcidOzf379+nffv2SbZI+/n5ERAQQI0aNeIlTQDHjx8nNDQ00fM2btzIjBkzaNCgAf/++y+Ojo707t2bp0+fGvReksPZ2ZknT54kaA0KDg7m2rVraR6PIZmZmTFnzhz69+/PsmXL6NmzZ7zZ4atWrYqiKBw+fNiEUWocHR0pVKgQ169fT/SR5r59+4D/HoPqU75gwYJcv36dR48eJSh/4MABA92JaUniJIQQaczCwoKNGzeybt26ty5mHitnzpzY2Nhw8uRJQkJC4vb7+/vz2WefJXrOw4cP+eCDD3BxcWHRokUUKlSIGTNm8OjRoyQ7oRtLlSpViIyMZMmSJXH7VFVl2LBhqX5ElR4oisKsWbMYMGAAy5Yto3v37nHJU+7cuenUqROHDh3ip59+SnTwz9GjR+N9b42pd+/eREZGMmzYsHixnD9/nvnz5+Pk5BRvaoSUlu/ZsycRERGMGjUq3nV37NjB7t27jXZfaUmeHwkhhAlUrVqVqlWrvrOcmZkZH3/8MVOmTKF8+fK0atWKwMBAtm7diru7O3nz5o1XXlVVevfujZ+fH2vWrIlbJqRr165s3bqVRYsWMW3aND799FOj3FdiPv30U+bPn0///v3ZuXMnOXLk4MCBA7x48YLy5cubZCRcrFWrVr21Y3a3bt2S/WhTURT++uuvuCRKVVWWLFmChYUF06dP58qVKwwdOpRFixZRs2ZNnJycuHv3LidOnODatWs8fPgwTZa9GTp0KJs3b2bRokVcunSJRo0a8fTpU1asWEFkZCR///03Dg4OqSq/du1aZs+ezYULF6hXrx53795l5cqVtGzZks2bNxv9Ho1NEichhEjnJk6ciIuLCwsWLGD69OnkypWLLl26MHbsWMqUKROv7JQpU9i5cyf9+/ePNzEhwJ9//snBgwf55ptvqF+/foJzjaVs2bJs27aN4cOHs3r1auzt7WnRogU//fQTnTt3TpMY3ubkyZOcPHky0WMVKlRIUZ8wRVGYOXMmZmZmzJw5E1VVWbp0KS4uLhw6dIhp06axYsUKlixZQkxMDLlz56Z8+fJ8//33uLq6GuqWkmRtbc2ePXv43//+x4oVK/jll1+wtbWlXr16DB8+PMHkoSktb2dnx759+xg2bBjr1q3j5MmTlC5dmhUrVhAQEJApEidFlUmDjCYwMBAnJycCAgKSMaouilKjtKUKLo7zks7kQgghhJ5S8vc3paSPkxBCCCFEMkniJIQQQgiRTJI4CSGEEEIkkyROQgghhBDJJImTEEIIIUQySeIkhBBCCJFMkjgJIYQQQiSTJE5CCCGEEMkkiZMQQgghRDJJ4iSEEEIIkUySOAkhhBBCJJMkTkIIIYQQyZQuE6fFixfz0UcfUaVKFXQ6HYqisGDBgkTLKoryzu3u3bvJuq6Hh8db6xg4cKAB71AIIYQQGZGFqQNIzMiRI/H19cXV1ZU8efLg6+v71rKjR49OdP/169dZsmQJJUuWpECBAsm+tpOTE19++WWC/VWqVEl2HUIIIYTInNJl4jRnzhyKFi2Ku7s7kyZNYtiwYW8tO2bMmET3f/bZZwD0798/RdfOli3bW+sUQgghRNaWLhOnxo0bp+r8sLAwlixZgpWVFT179jRQVEIIYRiLFi3izp07DB06FEtLS1OHI4RIgXSZOKXW2rVr8ff3p0OHDuTIkSNF54aHh7Nw4ULu37+Ps7MztWrVonz58sk+Nzw8PO51YGBgiq4thMj8fHx86NOnD/b29oSHhzNu3DhThySESIFMmTjNnTsXSPljOoBHjx7Rp0+fePuaNWvGokWLcHV1TfLciRMnMnbs2BRfUwiRNYSGhtKvXz8GDx5M48aNef/992nbti0VK1Y0dWhCiGRKl6PqUuPWrVvs3bsXNzc3mjRpkqJz+/Xrh7e3N0+fPiUwMJAjR47QvHlztm3bRuvWrVFVNcnzhw0bRkBAQNyW3NF8QoisYfjw4WTPnp2JEyfSrFkzBg8eTN++fYmMjDR1aEKIZMp0LU7z5s1DVVX69u2LmVnK8sJRo0bFe129enX++ecfPD098fHxYcuWLbRs2fKt5+t0OnQ6nV5xCyEyv19++SXe60mTJjFp0iQTRSOE0EemanGKiYlhwYIFmJmZ0a9fP4PUaWZmRt++fQE4ePCgQeoUQgghRMaUqRKnbdu2ce/ePZo0aYKbm5vB6o3t2xQSEmKwOoUQQgiR8WSqxCk1ncKTcvToUUCbWVwIIVLj+fPnfPnll7i7u2NtbU2ZMmVYtmyZqcMSQiRTpkmcnj59yqZNm3B1daV169ZvLRcZGcnly5e5ceNGvP0XL17kxYsXCcr7+PgwdepUdDod7dq1M3TYQogs5OrVq5QrV465c+fSuHFjBg0axOPHj+nWrRv//PNPmsZy7NgxWrRogbOzM3Z2dlSrVo2lS5cm+/z79+/z66+/0rRpU9zc3LCysiJ37ty0b98+7sOmIa6bkiW4wHRLZ6X2/dS3npTc74sXL/j888+pWbMmuXPnRqfTkS9fPho2bMiaNWveOgDKUPeWWaTLzuFz5szBx8cHgHPnzsXt8/b2BqBNmza0adMm3jl///03kZGR9OrVCysrq7fWff/+fUqWLIm7uzu3b9+O279y5UomT55Mo0aN8PDwQKfTcf78eXbs2IGZmRkzZ8406OM/IUTWEhQURMuWLYmOjubkyZMULVoUgN69e1OpUiV+/PFH3nvvvTSJxdvbGy8vL6ysrOjSpQtOTk6sXbuW7t27c/v2bYYPH/7OOv744w/+97//UbhwYZo0aULOnDm5du0a69evZ/369SxbtoxOnTql+ropWYIrVlovnWWI9zM19ST3fv38/Jg3bx41atSgTZs2uLi48OTJEzZt2kSHDh0YMGAAs2bNMsq9ZSpqOtS7d28VeOs2evToBOeULFlSBdSLFy8mWfetW7dUQHV3d4+339vbW+3UqZNapEgR1cHBQbW0tFTz58+vdunSRT169Khe9xEQEKACakBAwDvLBodHqu7f/qO6f/uPGhweqdf1hBDp1zfffKMC6saNGxMcK126tKooihoeHm70OCIjI9XChQurOp1OPXnyZNz+wMBAtXTp0qqFhYV69erVd9azZs0adf/+/Qn279+/X7W0tFRdXFzUsLCwVF93586d6u3bt1VVVdWJEyeqgDp//vy3xuXu7p7g93tKeXp6JrsOQ72f+taTkvuNiopSIyMT/n0JDAxUS5UqpQLq+fPnDX5vppCSv78plS4Tp8xCEichhKqqqr+/v2pjY6OWK1cu0eP16tVTAfXevXtGj2X79u0qoPbt2zfBseXLl6uAOmzYsFRdo2nTpiqgHjt2zKDXTY+Jk6HeT33rMcT9qqqqDh48WAXU9evXpzqm9MCYiVO6fFQnhBCZyapVqwgNDaV3796JHg8LCwNIspuBocR2eWjatGmCY7H79u3bl6prxK6/Z2Hx35+YtLhurNQsnZVShrqv1NST2vsNCwtjz549KIpCqVKlDBJTZiaJkxBCGNnWrVsBuHTpEmPGjElw/ObNm1hbW5M9e/ZEz0/snKR8+eWXZMuWLdFj165dA4jrY/U6Z2dnXF1d48ro486dO+zatYvcuXNTtmzZNLvu61KydFZi7+3t27d58eJFosfefG8NdV+pqSelS4W9ePGCX3/9lZiYGJ48ecKWLVu4e/cuo0ePjnf9tPyeZSSSOAkhhJHFTp47Z86ct5YpX778W1c7SOkamH369Hlr4hQQEABoHYoT4+joyL1791J0vViRkZH07NmT8PBwJk+ejLm5eZpc93X9+vXD09OT0qVLo9PpuHjxImPHjmXr1q20bt2agwcPoihKXPmk3tvEjr353hrqvvStJ6X3C1ri9Pq9WVpa8tNPPzFkyBCDxJTZZZrpCIQQIj16/vw5T548wdPTE1XrVxpv27x5MwC1a9d+ax2JnZfUZoo552JiYujXrx/79+9nwIAB9OzZM81jAG3pLE9PT1xdXXFwcIhbOqtOnTocPnyYLVu2xCuf2Pvn6emJu7t7unlvk5LS+wVtCgNVVYmKiuLWrVuMGzeOESNG0L59e6KiokxwFxmLJE5CCGFE9+/fByBPnjyJHt+2bRsAzZs3T5N4YlsPYlsT3hQYGPjWFoa3UVWVAQMGsHjxYnr06MHMmTPT5LrJZcylswx1X4Z8f5J7v+bm5nh4ePDdd98xYcIE1q1bx+zZs40SU2Yij+qEEMKIIiMjARJdADwiIoKVK1eSK1cuvLy83lqHIfs4xfZXuXbtGpUrV453zN/fHz8/P2rVqpXsa8XExNC/f3/mz59P165d49YLNfZ1U8pYS2cZ6r4M/f6k9H6bNm3K0KFD8fb2ZtCgQUaJKbOQxEkIIYwoV65cADx+/DjBsVmzZvH48WOmTp0aNxItMYbs4+Tp6cnEiRPZsWMHXbp0iXdsx44dcWWS4/WkqXPnzixatChevyZjXVcfxlo6y1D3Zej3J6X3++DBAyD+SEhTf8/SLYNPcCDiyDxOQghVVdXChQurNjY26p07d+L2HT58WLWzs1OrVauW6KSExhIZGakWKlRI1el06qlTp+L2vz6p4ZUrVxKcd/36dfXSpUtqRESEqqqqGh0drfbp00cF1I4dO77zHvS97uveNY/ThQsXVH9//wT7Dxw4oFpbW6s6nU719fVN8hoppc99vfle6ltPSu/31KlT6osXLxKUf/bsmVqhQgUVUBctWpSqmNILmcdJCCEysGHDhtG/f39q1apFly5dePToEStXrqRw4cJs2LAh3qd8Y7OwsGDOnDl4eXlRt25dunbtiqOjI2vXruXWrVtMmDCBYsWKJTivUaNG+Pr6cuvWLTw8PBg3bhwLFizA3t6eYsWKMWHChATntGnThgoVKqTquilZgkufpbNS+xhUn/t6873Ut56U3u+CBQuYM2cODRo0wN3dHTs7O3x9fdm8eTNBQUG0b9+ebt26peresgSDp2IijrQ4CSFiTZ8+XS1SpIhqZWWlFipUSB0xYoQaFBRksniOHj2qNmvWTHVyclJtbGzUKlWqqIsXL35reXd3dxVQb926parqu5fG4i0tQym9bkqW4NJn6ax33cObW+z9p+a+3nwv9a0npfd74MABtU+fPmqJEiVUR0dH1cLCQs2ZM6farFkzdenSpWpMTEyq7y29MGaLk6Kqb1kOWaRa7IiDgIAAHB0dkywbEhFFqVHbAbg4zgtbK2kMFEIIIfSRkr+/KSXTEQghhBBCJJMkTkIIIYQQySSJkxBCCCFEMkniJIQQQgiRTJI4CSGEEEIkkyROQgghhBDJJImTEEIIIUQySeKUBl6GhJs6BCGEEEIYgCROaeD7JcdNHYIQQgghDEASpzSw068AYRFRpg5DCCGEEKkkiVNayObGmKX/mjoKIYQQQqSSJE5pZMWt7MTEyLKAQgghREYmiVMaUMNfEuNSnMmrpa+TEEIIkZHplTgFBQVx584doqLi99tZsWIF3bt3Z8CAAZw+fdoQ8WUKxSNOAzD3nE5anYQQQogMTK/E6dtvv6VUqVKEh/83zH7GjBl069aNZcuWMXfuXOrWrcuVK1cMFmhGNqFjcdTIMCJcyzFn21lThyOEEEIIPemVOB04cIDGjRtjZ2cXt2/ixInky5eP/fv3s3LlSqKjo/npp5/0Cmrx4sV89NFHVKlSBZ1Oh6IoLFiwINGyY8aMQVGURDdra+sUX/vYsWO0aNECZ2dn7OzsqFatGkuXLtXrPmKVLpgT9yCtc/ivhyJTVZcQQgghTMdCn5Pu379P48aN416fO3eOe/fuMXnyZOrUqQPA6tWr2bdvn15BjRw5El9fX1xdXcmTJw++vr7vPKd37954eHjE22dhkbLb8/b2xsvLCysrK7p06YKTkxNr166le/fu3L59m+HDh6eovtdNbO9Bt51RBOeowqr9l+hYr6TedQkhhBDCNPRKnEJDQ7Gysop77ePjg6IoNG3aNG5foUKF2Lhxo15BzZkzh6JFi+Lu7s6kSZMYNmzYO8/p06cP9evX1+t6AFFRUfTv3x9FUdi/fz8VK1YEYPTo0dSsWZPRo0fTsWNHihYtqlf9dcq6kXPVQZ5mr82EHQF0rKd3qEIIIYQwEb0e1eXPn5+zZ//rq7N582acnZ0pW7Zs3L5nz55hb2+vV1CNGzfG3d1dr3P1tWfPHm7cuEG3bt3ikiYABwcHvv/+e6Kiopg/f36qrjH+vVyoMTG8yFGDTYevpjZkIYQQQqQxvVqcmjdvzp9//sk333yDtbU127Zto2fPniiKElfm8uXLuLm5GSzQdzlw4AD//vsv5ubmlChRgsaNG6PT6ZJ9vre3N0C8VrNYsfv0ffQYq3m1ImTfdIjnrrUYteUZrWqmqjohhBBCpDG9Eqdhw4axadMmpkyZAkDu3LkZO3Zs3PE7d+5w8OBBPv/8c8NEmQyjRo2K9zpPnjwsXLiQJk2aJOv8a9euAST6KM7Z2RlXV9e4Mm8THh4eb6RhYGBggjLjWrjy6b/wzKU6W45ep0X1IsmKTwghhBCmp9ejuty5c3PhwgU2btzIxo0bE7QuvXz5kilTpvDhhx8aLNC3qVChAgsXLuT27duEhoZy7do1xo8fz4sXL2jdujVnzpxJVj0BAQEAODk5JXrc0dExrszbTJw4EScnp7itQIECCcq0qlmMbE+PoJiZ8f0/T5IVmxBCCCHSB71anABsbGx47733Ej1WunRpSpcurXdQKdGmTZt4r4sUKcLIkSPJlSsXH374IRMmTGDVqlVpEsuwYcP46quv4l4HBgYmmjyNbubM4BPw1KUGO47fpGmVQmkSnxBCCCFSJ1VLrkRERLBlyxamTp3K+PHj4/aHhYXx5MkTYmJiUh2gvnr37o2FhQUHDx5MVvnYlqa3tSoFBga+tTUqlk6nw9HRMd6WmHZ1iuP09CiKmRnDNz5MVnxCCCGEMD29E6eNGzfi5uZGq1at+PrrrxkzZkzcsbNnz5InTx6WL19uiBj1YmVlhYODAyEhIckqH9u3KbF+TP7+/vj5+ek9FUFivm+qJVVPnGuw59Qtg9UrhBBCCOPRK3E6ePAgHTp0QKfT8dtvv9GtW7d4x6tVq0aRIkVYs2aNQYLUx7Vr1/D3908wKebbeHp6ArBjx44Ex2L3xZYxhI71SuLw9F8UM3O+W3ffYPUKIYQQwnj0SpwmTJhAtmzZOH78OJ9++mmiLTGVK1dOdsdsfb18+TLefFKx/P39+eCDDwDo2rVrvGORkZFcvnyZGzduxNvfqFEjChUqxNKlS+MtUPzy5UvGjx+PhYUFffr0MWj8Ixpq81w9ylYDn3N3DVq3EEIIIQxPr87hR44coUOHDuTIkeOtZQoUKJCqmcN9fHwAbTmX2H2xcy21adOGNm3a8OzZM8qXL0+VKlUoW7YsOXPm5P79+2zdupVnz57RpEkTBg8eHK/u+/fvU7JkSdzd3bl9+3bcfgsLC+bMmYOXlxd169ala9euODo6snbtWm7dusWECRMoVqyYXvfzNl0blmLCnmME5ajKiPX3DFq3EEIIIQxPr8QpPDz8nR2lAwICMDPTrwuVj48PCxcujLfv4MGDcR29PTw8aNOmDS4uLnzyySccOXKETZs28eLFC+zs7Chbtiw9evSgf//+mJubJ/u6DRo0wMfHh9GjR7Ny5UoiIiIoXbo048ePp3v37nrdy7uMbGjHd+fgsXM1zJ7sMso1hBBCCGEYiqqqakpPKlOmDK6urnEtQGPHjmXcuHFER0fHlSlfvjw6nY5///3XYMFmNLEj8QICAt46wg6g7MijBGavjPJ4OwAXx3lha6X3TBFCCCFElpbcv7/60KtJqH379hw4cIC///470eM///wz58+fp3PnzqkKLqsY3zybqUMQQgghRDLo1eIUFBREjRo1uHTpEo0aNSIsLIyDBw8yZMgQDh8+zKFDh6hQoQKHDh1K0XpxmU1KMt4KIw/xIsofkBYnIYQQIjXSXYuTvb09Bw4coEuXLuzduxcfHx9UVeXnn3/m0KFDdOrUiV27dmXppCmlxrX8r6P92gNXTBiJEEIIId5Grxan1z179oxjx47x/PlzHB0dqVq1Krly5TJUfBlaSjLekIgoSo3S+jjZmufg4g/V0iJEIYQQItMxZotTqp8HZc+enWbNmhkiFvFKsGslFu44R++mZU0dihBCCCFek6q16oTxTNwfaeoQhBBCCPGGZLU49evXD0VR+PHHH8mVKxf9+vVLVuWKojB37txUBZgVqVFhhGYvx9R1JxjYsnySZW0szVEUJY0iE0IIIbK2ZPVxMjMzQ1EULl26RLFixZI9saWiKPHmdspq9O3jlBJV3J1ZNbCmJE9CCCHEKybv43Tr1i0A8uXLF++1MBwbS3OquDtz3Nc/Recd9/UnNDJapi8QQggh0kCy/tq6u7sn+VqknqIorBpYk9BIrYWu488HOG9TF+XFbU5+kx/rNxKjkIhoqkyQJVqEEEKItKRX5/BChQrx6aefGjqWLE9RFGytLLC1smD2h5VRwwJQXYrw/aJjcfv/25K/Bp8QQgghDEOvxMnPzw8HBwdDxyJek9fVkfo25wHY9KwI/i9DTRyREEIIIfRKnCpUqMDVq1cNHYt4wx/9q0PgPRSHPHw866ipwxFCCCGyPL0Sp2+//ZZNmzaxd+9eQ8cjXuNkb03b3LcBOBhZnjtPAkwbkBBCCJHF6TUU69mzZzRt2pQmTZrQtm3buGVWEhsS36tXr1QHmZVN7lOTDROuEeNSlIFzvNkyvL6pQxJCCCGyLL3Wqoud1+nNU19PnFRVlXmcDDSPxA/LjjLrSXXUiBA2tg2gQpE88eZ9ujjOS6YjEEIIIV4x+TxOb5o3b55MuJiGhnWuxoLRZ4lwLcfARSc4MjaPqUMSQgghsiS9Eqc+ffoYOAyRFDMzhdGeZoy4AA+y1Wb9wSs0rVrY1GEJIYQQWY5encP//vtvzp49m2SZCxcu8Pfff+sVlEioR+MyZPc7hGJmxrCtQcTEpPgJqxBCCCFSSa/EqU+fPqxfvz7JMv/88w99+/bVp3rxFtM6F0CNCickR2X+2Hja1OEIIYQQWY5eiVNyREdHJ3sxYJE8tcoUoEToYQD+Om/Yzm5CCCGEeDejZTanTp3CxcXFWNVnWbP7V0ANeY7qLH2chBBCiLSW7M7hDRs2jPd6wYIFeHt7JygXHR3NvXv3uH37Np06dUp1gCI+99zZaGC7D++Y2qYORQghhMhykp04vZ4kKYrC7du3uX37doJyZmZmuLi40LFjR3799VcDhCje9OeHNSn9v7umDkOkkKqqhEambF4zG0tzmfpDCCHSkWQnTjExMXFfm5mZMWbMGEaNGmWUoETS7G2s6F7oEUsvaq9PXn1InTIFTBuUSJKqqnSYeZgTvv4pOq+KuzOrBtaU5EkIIdIJvfo47d27l969exs6FpECwztVifv68+W+JoxEJEdoZHSKkyaA477+KW6lEkIIYTx6TYDp6elp6DhECpmZ/dcC8Sx7DeZuO8MHzcqbMCKRXMdHNsbWyjzJMiER0VSZsCuNIhJCCJFcei9wFhERwfr16zl27BgvXrxIdE06RVGYO3duiutevHgxBw4c4MSJE5w7d46IiAjmz5+fYMbyyMhINm7cyKZNmzh69Ch37tzBzMyMUqVK0bt3bz766CPMzZP+A/U6Dw8PfH0Tb7356KOPmDlzZorvJa38eNCS3k1isDCXKSDSO1src1lbUAghMii9fnv7+vrSpEkTbty4kWCh39fpmziNHDkSX19fXF1dyZMnz1uTmRs3btChQwccHBxo2LAhrVu3JiAggE2bNvHJJ5+wbds2NmzYkKL+IU5OTnz55ZcJ9lepUiVh4XRCDQ8kyrUUg+cc4I+P6po6HCGEECLT0itxGjx4MNevX6dnz57069eP/PnzY2FhuE/Qc+bMoWjRori7uzNp0iSGDRuWaDkHBwemT59O7969sbW1jds/ZcoU6tevz6ZNm1i9ejUdO3ZM9rWzZcvGmDFjUnsLaaq21VkOUZ8NfsUZ5hdIXleZHFMIIYQwBr2ynT179tCoUSMWLlxo6HgAaNy4cbLK5cuXj0GDBiXYb2dnx1dffUW3bt3Yt29fihKnjOi3D6pT7ddbkK0gfWfuZfvIBqYOSQghhMiU9EqcYmJiqFixoqFjMShLS0uAFLeEhYeHs3DhQu7fv4+zszO1atWifPn03ena3saSQaX8mP6gIJdsarP39G0aVPAwdVhCCCFEpqNX4lSzZk0uXbpk6FgMat68eQA0bdo0Rec9evQoQSf0Zs2asWjRIlxdXZM8Nzw8nPDw8LjXgYGBKbp2agztUJW/vz9OUI4qfLb6CeclcRJCCCEMTq8hWJMmTWLv3r2sXr3a0PEYxKxZs9i6dSsNGzakRYsWyT6vX79+eHt78/TpUwIDAzly5AjNmzdn27ZttG7dOsmO8AATJ07EyckpbitQIO0mpVQU+KWNC2p0JC9zVOPH5UfT7NpCCCFEVqFXi9OmTZto0KABnTt3xtPTk4oVK+Lk5JSgnKIofP/996kOMiU2b97Mp59+iru7O4sXL07RuW/OhF69enX++ecfPD098fHxYcuWLbRs2fKt5w8bNoyvvvoq7nVgYGCaJk9NqxSi5DZvLjvU56+reRj0MhRnB5s0u74QQgiR2emVOL0+6szb2zvRxX4h7ROn7du30759e3LlysWePXvIkydPqus0MzOjb9+++Pj4cPDgwSQTJ51Oh06nS/U1U+PvT6pQ9c+HKE5u9P3Tm/Xf1TdpPEIIIURmolfitHfvXkPHkWrbtm2jbdu2uLq6snfvXgoVKmSwumP7NoWEhBisTmPJ5WxPT7dzLPbPw0mL6vicu0Odsm6mDksIIYTIFDLFkivbtm2jTZs2uLi4sHfvXooUKWLQ+o8e1foLeXh4GLReYxnfowZrR50kJEclBq54yHlJnIQQQgiDyPDrc8QmTc7Ozuzdu5eiRYsmWT4yMpLLly9z48aNePsvXrzIixcvEpT38fFh6tSp6HQ62rVrZ8jQjcbMTOG3NtledRSvzo/L/zV1SEIIIUSmoPd031FRUfzxxx8sW7aMy5cvExISQlRUFACnT59m1qxZfPnllxQrVizFdc+ZMwcfHx8Azp07F7cvti9VmzZtaNOmDZcvX6ZNmzaEh4dTv359li1blqAuDw+PeNML3L9/n5IlS+Lu7s7t27fj9q9cuZLJkyfTqFEjPDw80Ol0nD9/nh07dmBmZsbMmTNxc8s4LTdNqxSi1HZvLtnX56+ruaWjuBBCCGEAeiVOoaGhNG3alEOHDuHq6oqjoyPBwcFxxwsWLMj8+fNxcXFhwoQJKa7fx8cnwazkBw8e5ODBg4CWDLVp04ZHjx7FzZu0fPnyROvy9PRMMC9TYho0aMClS5c4efIk+/btIywsjFy5ctG5c2cGDx5MtWrVUnwfprbw4/86iveetpcN39UnNDLhYsxJsbE0T9Faf0IIIURmplfi9OOPP3Lw4EEmTZrEN998w9ixYxk/fnzccScnJzw9Pdm+fbteidOCBQtYsGDBO8vVr1//nXMrvcnDwyPRczw9PdNd363UyuVsT2/3c/z9PA+nLWvS+Oe93HgWmqI6qrg7s2pgTUmehBBCCPTs47RixQrq16/P0KFDURQl0T+qhQoV4s6dO6kOUKTO2O41sH96DMXcIsVJE8BxX/8Ut1IJIYQQmZVeLU537tyhbdu2SZZxdHQkICBAr6CE4ZiZKczpmofOm0OJTW+Pj2yMrZV5kueFRERTZcIu4wcohBBCZCB6tTg5ODjw9OnTJMvcuHGDHDly6BWUMKyapfNTw+xE3OsHT15ga2Xxji3pxEoIIYTIivRKnGrUqMGmTZve2qJ07949tmzZQr169VIVnDCc6R9Wj/u6//wrJoxECCGEyLj0Spy++eYbnj9/TuPGjTl06FDcNAQhISHs3r2bpk2bEhkZGW/dNmFattaWcV8/dqnF9I0nTRiNEEIIkTHp1cepXr16/Pnnn3z++efUrVs3br+DgwMA5ubmTJ8+ncqVKxsmSmFwk09mo0ejMBztrE0dihBCCJFh6D1z+MCBAzlz5gyffvopVatWpXDhwlSsWJGBAwdy6tQp+vfvb8g4hQGpwU9QnQvR7bfDpg5FCCGEyFD0njkcoGTJkvz222+GikWkkW4FfFn2Ii9nreuw/uAV2tQubuqQhBBCiAwhw69VJ1Lu+y5VyeZ3FMXckiHbogiLiDJ1SEIIIUSGIIlTFqQosKivB2pYAFGupen7xwFThySEEEJkCJI4ZVHlCuWijbO2gPJBtQZ7T90ycURCCCFE+ieJUxb2a//a2PqdRLG04aM1L4iKjjF1SEIIIUS6JolTFmZmpjCvS07UiGDCc1Tk4xnyyE4IIYRISrISp/3798uCvZlUzdL5aajTlmPZFlyRY5fvmzgiIYQQIv1KVuLUoEEDFixYEPe6YcOG/P3338aKSaSx2R/XwcrvHIq1I70WPSQmRjV1SEIIIUS6lKzEycLCIm5ZFQBvb29u375trJhEGrO0MOOvdg6okaGE5KzCoOn7TB2SEEIIkS4lK3EqUKAABw8eJCbmv87DiqIYLSiR9hpW9KCB5TEAtoRU5vCFuyaOSAghhEh/kjVzeJcuXfjxxx9xdnYme/bsAPzyyy/Mnz8/yfMUReHGjRupj1Kkibmf1qH0mDOEuZZn4Korpg5HCCGESHeS1eI0evRoJkyYQLly5VAUBUVRUFX1ndvrLVQi/bMwN2NuZxfUiGAiclQwdThCCCFEupOsFidLS0uGDx/O8OHDATAzM2Pw4MGMGjXKqMGJtFenTAGa7fdhe3gNU4cihBBCpDt6zeM0evRo6tevb+BQRHoxY2BtbJ6diXsdLmvZCSGEEEAqEqd69eoZOhaRTpibKczpmivudd9pB00YjRBCCJF+pGrm8EOHDvHhhx9SrVo1ihcvTtWqVfnwww/x8fExVHzCRCoVzR339UnLWizbc96E0QghhBDpQ7L6OCXm66+/5pdffkFVtckSzczMiImJ4cSJE8ydO5cvvviCqVOnGixQYTqKmSXf7bOhcaUgcmSzN3U4QgghhMno1eL0999/M3XqVIoXL86yZct4+PAhUVFRPHr0iOXLl1OiRAl+++03mV08k1CDHoNLYdpOPWbqUIQQQgiT0itxmjFjBgUKFODo0aN07tyZXLm0/jA5c+akU6dOHD58mPz58zN9+nSDBitM48syjwG4m70BPy47bOJohBBCCNPRK3E6f/487du3x8HBIdHjjo6OtGvXjgsXLqQqOJE+DGxZjqJBBwCYebMwl24/NnFEQgghhGno3Tk8tm/T26RmSZbFixfz0UcfUaVKFXQ6HYqixFtk+E2BgYF89dVXuLu7o9PpcHd356uvviIwMDDF1z527BgtWrTA2dkZOzs7qlWrxtKlS/W+l8xizZfVMPO/jmKfk46zbhIdLZObCiGEyHr0SpzKlCnDmjVrCAoKSvT4y5cvWbNmDaVLl9YrqJEjRzJr1ix8fX3JkydPkmWDg4Px9PTkl19+oXjx4gwePJhSpUrxyy+/4OnpSXBwcLKv6+3tTZ06dThw4AAdOnRg0KBB+Pn50b17d3788Ue97iWzcLLT8XszUKPCeZmrJh/8vtfUIQkhhBBpTq/EaeDAgdy7d4+aNWuyZs0a/Pz8APDz82P16tXUqlWLe/fuMWjQIL2CmjNnDrdv3+bp06cMHDgwybKTJ0/m9OnTDB06lB07djBp0iS2bt3KqFGjOH36NJMnT07WNaOioujfvz+KorB//35mz57Nzz//zJkzZyhdujSjR4/m2rVret1PZtGqRhEaWWkdxPdE12btfnkUK4QQImvRK3Hq3bs3X3zxBRcuXKBTp07kypULS0tLcuXKRefOnblw4QKffvopvXv31iuoxo0b4+7u/s5yqqoyZ84c7O3tEyz/MmzYMJydnZk7d+47HysC7Nmzhxs3btCtWzcqVqwYt9/BwYHvv/+eqKiody5qnBXM/aQ2Ts+Oo1haM3i3jsfPX5o6JCGEECLN6N3H6ZdffmH//v306dOHChUq4OHhQYUKFejbty/79u3jt99+M2Scibp27RoPHjygdu3a2NnZxTtmbW1NvXr1uH//PtevX39nXd7e3gA0bdo0wbHYffv27Ut90BmcmZnC2oGFX01RUITWU48nKzEVQgghMgO9J8AEqFOnDnXq1DFULCkW++isaNGiiR6P3X/t2rW3lklOXc7Ozri6ur7zUV14eDjh4eFxr/XpnJ4RFMnrzMhK55lwOQePcjTgmzn7+HmAp6nDEkIIIYwuVUuumFpAQAAATk5OiR53dHSMVy61db2rnokTJ+Lk5BS3FShQ4J3Xzag+bF6GSlGHAFjxvALep26aOCIhhBDC+DJ04pTeDBs2jICAgLjt7t27pg7JqFZ8WRPrZ+cxs3ai3+ogAoJCTR2SEEIIYVQZOnGKbR16W0tQ7KOyt7UipbSud9Wj0+lwdHSMt2VmOktzlvbMgRr6guic5XjvfwdNHZIQQghhVBk6cXq9D1Ni3tUHKrl1+fv74+fnl6x6sprKRXPxefHbANxxbczQ2d4mjUcIIYQwpgyfOOXNm5eDBw8mmOgyLCyM/fv3kzdvXooUKfLOujw9tc7NO3bsSHAsdl9sGRHf120rUDZca21a9rwS249eMXFEQgghhHFk6MRJURT69+9PUFAQ48aNi3ds4sSJ+Pv7x01qGSsyMpLLly9z48aNeOUbNWpEoUKFWLp0KadPn47b//LlS8aPH4+FhQV9+vQx5u1kaGsG18Tm+XnMrB356J8YnvrL/E5CCCEyH72mI2jYsCF16tRJkKwYypw5c/Dx8QHg3Llzcfti51pq06YNbdq0AWDo0KFs3LiRyZMnc+rUKSpXrsyZM2fYunUrFSpUYOjQofHqvn//PiVLlsTd3Z3bt2/H7bewsGDOnDl4eXlRt25dunbtiqOjI2vXruXWrVtMmDCBYsWKGeV+MwOdpRlr+uWh+RI/cC1Ji5928+8PDVO1ZqEQQgiR3ujV4nT06FGioqIMHUscHx8fFi5cyMKFCzl58iQABw8ejNv3eouQnZ0d3t7eDB48mMuXLzNlyhTOnz/P4MGD8fb2TjAxZlIaNGiAj48PderUYeXKlUyfPp3s2bOzePFiRowYYejbzHRKu2dnRMXHqDExPMnViEF/7DJ1SEIIIYRB6dXiVLJkyXitNYa2YMECFixYkOzyTk5OTJ06lalTp76zrIeHR5IzXVerVo2tW7cm+9oivo+al2bP5cMcoSZbIuqweMcpejSt+O4ThRBCiAxArxanzz77jI0bN3Lx4kVDxyMygSWf1yCb/ykUSxuGH87OxZsPTR2SEEIIYRB6tTgVLFiQ+vXrU6NGDT766COqVq1Krly5Eu3PUq9evVQHKTIWC3OFLZ8Vo870u8Rkc6PNrCOcHu2CrY3O1KEJIYQQqaJX4lS/fn0URUFVVaZMmZJkB+Do6Gi9gxMZV77sdvz13nP67w4hPE8N3vtxK3vGNzd1WEIIIUSq6JU4jRo1SkZLiXdqWrEAH1w9zbxHFbjh2pyvZ+7k54FNTB2WEEIIoTe9EqcxY8YYOAyRWY3uXIF/fzrKeavqrAisQU3vM7SvX97UYQkhhBB6ydATYIqMYe2XVbF/cR4znQOD99pz8eZ9U4ckhBBC6EWvFqdYp06dYtmyZVy+fJmQkBB27dLm7fH19eXo0aM0btwYFxcXgwQqMi6dpRmbB3lQf9YDcClM67+OcPJ7ZxzsbAiNNG4fOBtLc3msnEGoqprinwf5/goh0preidPQoUOZMmVK3JxIr//yUlWVbt26MWXKFL744ovURykyPI+c9sx6L4D+O4OJzFuDxuM3U8DDmRO+/ka9bhV3Z1YNrCl/XNM5VVXpMPNwin8e5PsrhEhrej2qmz9/Pj///DPvvfceZ8+eZdiwYfGOe3h4UK1aNTZu3GiQIEXm0LRCPr4o7gvA49xeRk+aAI77+hu9VUukXmhktF4/D/L9FUKkNb1anKZPn07JkiVZs2YNFhYWWFlZJShTokSJuEd3QsQa0roUp6efZH94ubh9x0c2xtbK3KDXCYmIpsoE+fnLiJLz8yDfXyGEqeiVOF28eJEBAwZgYfH203PlysWTJ0/0DkxkXgsHVaLWj6d49Oq1z4nLtPGUkXZCY2tljq1VqrpfCiGE0ej1qM7CwoKIiIgkyzx48AB7e3u9ghKZm5kC6z8tGff6s+0WXLpx14QRCSGEEMmjV+JUtmxZ9u7dS0xMTKLHY0fYVa5cOVXBiczL0ea/FgWzHMVpOcOXZ/4BJoxICCGEeDe9Eqd+/fpx5coVBg0alKDlKTAwkD59+vDo0SMGDBhgkCBF5qZGhRFdoA4Nxu0jIiLS1OEIIYQQb6V34tS1a1dmz56Nq6src+fOBaBatWrky5eP1atX07t3bzp06GDQYEXmNKTMQ9SYGAI8WtPy+zVxU1wIIYQQ6Y3eM4cvWbKEv/76i4IFC3L//n1UVeX48eO4ubkxY8YM5s2bZ8g4RSY2oFFB2ue4DMCV3J0Y+NNqE0ckhBBCJC5VS64MGDCAM2fOEBQUxL179wgMDOTChQt89NFHhopPZBFTu5eigtkFFDMzttCC/y3YYuqQhBBCiAQMsladjY0NefPmlVF0Qm+KAqs/Lk2u0CuY6ez4415lFm3cZ+qwhBBCiHhSNVlKcHAwGzZs4PTp0wQEBODk5ESFChV4//33sbOzM1SMIouwNIddXxalxtTbBDt4MOxYILmzn6RJ7UqmDk0IIYQAUpE4LVu2jE8//ZQXL17E68yrKArZsmXjzz//pEuXLgYJUmQdjtZm7BiUF8+Zj8C1KH03nGCz8zXKlypq6tCEEEII/R7Vbdq0iR49ehAWFsbHH3/M8uXL2bt3L8uXL2fQoEGEhYXRo0cPNm/ebOh4RRaQ39mKdT0dUMJfYJ6vMq3/8uXOvQemDksIIYTQL3GaMGECDg4OnDp1ij/++INOnTrh6elJp06dmDZtGidOnMDOzo7x48cbOl6RRZTLb8fsljEQFQaFGtPwh0M8f278RYGFEEKIpOiVOJ07d44uXbpQrFixRI+XKFGCLl26cPbs2VQFJ7K2JqVdGF8zADUmmsgSHaj73XqCgoJMHZYQQogsTK8+To6OjmTLli3JMtmyZcPJyUmf6kU6ExIRneyyNpbmKIpisGv3qpWLR4H3+fNKPoJK98Vz8HQO/tEPa2vrZNdhrPhVVSU0Mnl1pySGzELeHyFEZqRX4tSqVSv++ecffvjhB8zNzRMcj4qKYvPmzbRu3TrVAQrTqzJhV/LLujuzamBNgyZPQ5vl43nIXZbdLYBfmY9p/Pmv7P3zEywtLZMXkxHiV1WVDjMPc8JXHh8mRt4fIURmpdejup9++glra2uaN2/O0aNH4x07cuQIzZs3x8bGhv/9738GCVKkPRtLc6q4O6f4vOO+/sluZUiJiW0L0CzHXQDulPiU9774jejot1/H2PGHRkbrlRRUcXfGxjLhh43MRt4fIURmlawWp0KFCiXYFxERwalTp9i9ezeWlpZkz56dZ8+eERmpLdKaJ08eKleuzI0bNwwbsUgTiqKwamDNFD1qSUnLTsrjgRldC9B1/j2OvMzPpUKf0mnIVFZNHYKZWcL8Py3jPz6yMbZWyftjb+hHmRmBvD9CiMwkWS1OMTExqKoab7O0tMTNzQ03Nzfy5MmDlZUVefLkidtnaWlJTEyMseMHYMGCBSiKkuTWqFGjd9bj7e2dZB1HjhxJg7tJPxRFwdbKIpmb8VsJzBRY0ic/pa3vo1ha82+egXT7avJbf87SKn5bK/NkXycrJgXy/gghMpNktTjdvn3byGGkToUKFRg9enSix1avXs2FCxfw8vJKdn2enp7Ur18/wf78+fPrG6IwEAszWPtBPpr/9Yib1rk5mOsjeg6ZxOKpw+SPrhBCCKNL1ZIr6UWFChWoUKFCgv0RERFMmzYNCwsLevfunez66tevz5gxYwwXoDAoawv4Z0Bums16wh3bnOxzHUDvr35g4dQRkjwJIYQwKoMs8pterVu3jmfPnvHee++RK1cuU4cjDMjOCjYPyElecz/M7XOwx/kD+g0ZH2/5HyGEEMLQ9G5xCgoKYu7cuZw5c4b79+/HdQp/naIo7N69O1UBpsbcuXMB6N+/f4rOu3btGr///jshISG4u7vTpEkTXF1djRGiSAVHHWzt70rT2c957JSHHTF96P/VGGZPGZ1oh3EhhBAitfRKnE6cOEGzZs14/vx5kp/wTfnYxNfXl927d5MvXz6aNWuWonOXLl3K0qVL417b2NgwduxYvvnmmyTPCw8PJzw8PO51YGBgyoIWKZbNGrZ84ILX3Bf4ObuxLaYPPT/5jr+nTUx0jjEhhBAiNfT6WP7ZZ5/h7+/PpEmTuHPnDpGRkcTExCTYkppnx9jmz59PTEwMffv2TfYf0Bw5cvDTTz9x6dIlgoODuX//PosXL8bFxYWhQ4fy119/JXn+xIkTcXJyitsKFChgiFsR7+BqC5v7ZsPFLBCL7AXxzjGIjh8MTrQVVAghhEgNvRKnU6dO0aVLF7755hvy58+f7j7Zx8TEMH/+fBRFoV+/fsk+r3Tp0nz99deUKFECW1tb8ubNS/fu3dm2bRtWVlaMHj06ySkWhg0bRkBAQNx29+5dQ9yOSIbc9rCljyOu5i+xyF6Qo26Deb/Hx/FaAIUQQojU0itxyp49Ozly5DB0LAazc+dO7ty5Q8OGDSlYsGCq6ytTpgzVq1fn8ePHXL9+/a3ldDodjo6O8TaRdvI4wD+9HchhGYxF9oKcKT6cFp0HEBISYurQhBBCZBJ6JU7t2rVjz549aTbBZUrp2yk8KbGdw+WPcPqWxwE29bQjp1UIFtkLcrncWBq26Ym/v6yZJoQQIvX0Spx+/PFHdDod3bt35/79+4aOKVWePXvGhg0bcHFxoW3btgapMyoqipMnT6IoCm5ubgapUxhPHgfY2MOWXLpQLLIX5G7NqdRu0Y179+6ZOjQhhBAZnF6Jk729PX/99Rfbt2/Hzc2N7NmzU6hQoQRb4cKFDR3vOy1atIiIiAh69OiBTqdLtIyfnx+XL1/Gz88v3v7Dhw8nGCUYFRXFN998g6+vL15eXri4uBgtdmE4WvJkQz7bcCxc3AloPpda7/Xi0qVLpg5NCCFEBqbXdAS7d++mVatWhIWFYWlpia2tbaLTEphiMsLkPKabNm0aY8eOZfTo0fFmCO/atSuKolCrVi3y5cvHixcv2L9/P1euXMHNzY2ZM2caO3y9hES8e/RicsoYSnqJJ7c9bOimo+PKCG6Rl+gOy6nXvgOb5k2mRo0aRr9+eqeqarIXQY4li/C+nbyfQmQNeiVO3377Laqqsnz5cjp06JBuJhv8999/OX/+PNWqVaNs2bIpPn/QoEFs27YNb29v/Pz8sLCwoEiRIowYMYIhQ4bg7OxshKhTr8qEXaYOIZ70FE8OO1jXxYouq6O4TE50PdbTtGdrFk8ZSuvWrU0dnsmoqkqHmYc54Zuyvl9V3J1ZNbCm/LF/g7yfQmQdeiVOFy9epEePHnTq1MnQ8aRKtWrVktXKNWbMmETXovv222/59ttvjRCZ4dlYmlPF3ZnjevyitrE0/PQR6S2e1znbwMpOFvRaG83pJy44fbCZTl+1YuKNG3z55ZdZ8o9WaGR0iv/IAxz39Sc0Mhpbq0yxzKXByPspRNah1//WHDlyYGNjY+hYRAooisKqgTXTzaOB9BbPm5x0sLS9OX03xHD0gRM5PtrO8Dkd4pbXycqOj2yMrVXSyWtIRHS6aklMz+T9FCJz0ytx6t69O6tWrSI0NFQSKBNSFCVdfVJNb/G8yc4K/m5rxsebVXbftiHHBxtYuLQ3N1q2ZMHiZaYOz2RsrczT9fcto5H3U4jMTa/OSWPGjKFMmTJ4eXnh4+NDUFCQoeMSwiisLeCv9xTalgDF3ALXnks4FFacRo0amjo0IYQQGYBeH4tiW5lUVcXT0/Ot5RRFISoqSr/IhDASS3OY2hScrWHeaXBp9zsPdownm6kDE0IIke7plTjVrVs3S3aoFZmHmQKj6kE2a5h6BJwaD4PH2wH4+eef+X7Yt/IzLoQQIgG9Eidvb28DhyFE2lMU+KK61vI0as9/ozHHTvgfZ08eZ/78+Tg4OJgwQiGEEOlN+piASQgT6lUefmv2X+tSzoFbWbdtH9WrV+fChQsmjEwIIUR6I4mTEECT11YHsipQhXxDjnL1SRhVq1Zlzpw5JpkFXwghRPqj16O6hg2TNwJJURR2796tzyWEMJm8DvDArBBuQ09w/08vBgwYwK5du5g1axaOjo6mDk8kISYmhjt37nDlyhWePXtGQEBA3BYUFISNjQ1OTk44Ojri6OiIi4sLpUuXxsPDI92sgCCESN+M0sdJURRUVZXOtSJDWtYeBm2Fi0+dyTf4IE8X9WDFihUcO3aM5cuXU7VqVVOHKF7z66+/cuncGS5evMjly5cJCQlJcR329vaULVuWsmXLUrlyZVq0aEH+/PmNEK0QIqPTK3GKiYlJdH9gYCAnT55k+PDh5MuXj+XLl6cqOCFMIacdrOoAn22FPbctcem1ApeCFbm5Zhi1atVixIgRDP4mYyzNk5moqsrx48fZvXs3e/b7QLlBAIwcORI1MjyunKWlJUWLFiV37tw4OTnFbfb29oSGhhIYGEhAQACBgYE8evSIS5cuERQUxOHDhzl8+HBcPZUrV+b999/n/fffp2zZsvJBUAgB6Jk4vY2joyP169dn+/btlC1blh9++IFRo0YZ8hJCpAl7K5jTCsYfgPmngbrfUblYVU78rxljx45lw+at0Eh+to0tKiqKAwcOsHbtWtavX8+9e/cAUCx1uL1KnN57rxU1qlSkVKlSlCxZkkKFCmFpaZmia1y9epWzZ89y9uxZ9u3bx+HDhzlx4gQnTpxg1KhRFC1alE8//ZS+ffvKSEshsjijrAvg4OBA8+bNmT9/viROIsMyN4MxnuDuBOP2w9NcjWgw7R5nJtTizJkzuDXSykVFRYEssWEwqqri4+PDokWLWLt2Lc+ePYs7ZmdnR5MmTajj2ZA/Hmn7li9flqolTiwsLChVqhSlSpWiS5cuADx+/Jh//vmHjRs3smPHDq5du8YXX3zB999/zwcffMBnn31GwYIFU3WfQoiMyWi9Ic3MzHj48KGxqhcizfStoLU+2VnCjYhclBh3hYYdPow73rBBQ44fP266ADOJ69evM3r0aAoXLky9evWYPXs2z549I3v27PTt25eNGzfy9OlT1q1bx6CPBxk1lly5cvHBBx+wYcMGnj59yvTp0ylevDiBgYH88ssvFClShC5dunDz5k2jxiGESH+MkjjdvHmTVatW4e7ubozqhUhzjQrC2k6Q3xHuBVnwoNaUuGMnT52kWrVqfPLJJ/j7+5swyownNDSURYsWUa9ePYoWLcq4ceO4desW9vb29O3bl127dvHo0SPmzZtHq1atEl1UPCQimpCIqGRt+kwrYW9vz6BBg7h48SJbtmzBy8uLmJgYVqxYQcmSJRk6dCgBAYGGeDuEEBmAXu3b/fr1S3R/VFQU9+/fx8fHh8jISMaMGZOa2IRIV0q4wqYuMGgzHLmrENtVuFPnzqxY/DfTp09n9erV/Pzzz/To0UM6Eyfh/PnzzJo1i0WLFvHixQtAa6Vu0qQJvXr1ok2bNtja2iarrioTdiX7ulXcnVk1sKZe3xszMzOaN29O8+bNOXv2LF9//TU7d+7kp59+Yv6iJdj1npXiOoUQGY9eidOCBQuSPF6sWDG++uorPvzwwyTLCZHRuNjAorYwcjesfPxqX5e5bO3Tj8GfDuTy5cv06tWLGTNmMHnyZOrUqWPagNOR0NBQVq9ezcyZMzl06FDcfnd3d/r370/fvn3Jly9fsuqysTSnirszx31T1sJ33Nef0MjoVPWJAihXrhzbt29n27ZtDBkyhMvXb2L36tj58+epVqlCquoXQqRfev32uHXrVqL7zczMyJYtm4w6EZmalTmMqQ8r92ivN1+FWzk92XTgDGvm/sK4ceM4fPgwdevWpVWrVkycOJHSpUubNGZTunb1GgvnzWbBggVxjzItLCxo3bo1H374IY0bN8bc3DxFdSqKwqqBNQmNjE5W+ZCI6BS1TCU3hubNm9OkSRP+/GsOv9zV9tetU5cx3w/nm2++wcJCBg0Ikdno1cfJ3d090a1AgQKSNIksJ7sNXPSDdqutqN75W65du8ZHH32Eubk5mzZtoly5cvTr1++tHzgyo9DQsLivK1aqyC+//IK/vz/u7u788MMP3L17lzVr1uDl5ZXipCmWoijYWlkkc9PvGslhYWHBgAH9415HRkUyfPhw6tSpw5UrV4x2XSGEacgaA0Kk0uqOUDE3BIRDnw2w5l5eps+YyYULF2jfvj0xMTHMnz+fokWL0r17d86ePWvqkI3m3LlzfP755xQtWjRun5mZGa1atWLz5s3cuHGD4cOHkzt3bhNGaVyzZs3CycmJo0ePUqFCBf766y9Z61CITETvduSIiAjWr1/PsWPHePHiBdHRCZvMFUVh7ty5qQpQiPQutwOsaA9j9sHS8/DzYTj9CKY2Lc7q1as5cuQIo0aNYufOnSxdupSlS5fSvHlzvv32W+rVq5fhO5E/fvyYZcuW8ffff3Pq1ClAm6AydlW/S5cuUbRg1hlh261bN7waNeCDDz5g586dDBw4kCNHjjB9+vRERwUKITIWvRInX19fmjRpwo0bN5L8JCWJk8gqdBYwsRGUywWjvGHXLWixDKY3hxo1arBjxw5OnjzJ5MmTWbVqFVu3bmXr1q2UKlWKDz74gJ49e5IjRw5T30ayvXz5ks2bN7N48WK2bdsW98HJ0tKS1q1b06tffz7fr+1LbofvzKRAgQJs376dn376iWHDhrFgwQLOnj3LmjVr8PDwMHV4QohU0OtR3eDBg7l+/To9evRg7969XLt2jVu3biXYZHI4kdV0LaPN9+TmBPcCocNqWHgGVBUqVarE8uXLuXr1KoMGDcLGxoaLFy8yZMgQ8uXLR8eOHdm6dSsRERGmvo0kderYiRw5ctC1a1c2b95MdHQ01atX588//+Thw4esXr2axo0bmzpMk1MUhaFDh7Jjxw5cXV05efIklStXZseOHaYOTQiRCnq1OO3Zs4dGjRqxcOFCQ8cjRIZXNif80xW+2Qnbb2gtUMcewKRG2hp4hQsXZvr06UycOJHly5czZ84cjh8/zurVq1m9ejUODg40bdEKPLqZ9D6io6M5ceIEu3btYsceb6jyBQBbtm5BjQynaNGidOrUiZ49e1K8eHGTxpqeNWrUiBMnTtChQweOHTtGs2bNmPjTFKCYqUMTQuhBrxanmJgYKlasaOhYhMg0nHTwV0sYVQ8szGDTVWi5FM49fq2MkxMfffQRx44d48yZM3z22Wfkzp2bly9fsnbtmrhyTZt68d1337FhwwaePHlitJgDAgLw9vZm6tSptGvXDldXV6pXr86IESPYv39/XLnhw0dw7tw5rly5woQJEyRpSgY3Nzf2799P//79UVWVYcOGmTokIYSe9GpxqlmzJpcuXTJ0LEJkKooCH1SECrnh861wOwDaroShtaF/RTB7rU94uXLl+P333/n11185duwYazf8w4oY7dihQwc5uG9PXNlChQpRoUIFChcuTOHChSlUqBCFCxcmR44c2NnZYWaW+OchVVUJCQnh7t27XL15O25/z569OHPiX27cuJHgHCcnJxo2bEjdBo357b62b/jwYameQDIrsra2ZtasWZQoUYJvho2I2x8SEoqtlUzjIkRGoddvv0mTJlG3bl1Wr15Nhw4dDB2TEJlK5TywpRt8uxu2XocfDoDPHZjaFFzfWFXEzMyM6tWrU7ZiZVaM2g7A9OnTOXH0MIcPH+bixYvcvHkzyf6DdnZ22NvbY2dnR2RkJKGhoYSEhBAaGho3mEOx1OH2ldaqtW7dWtTIcAA8PDyoVKkSVapUoVGjRlSqVAkLCwtCIqL47VU8Qn+KojBkyBBy5S3AyDPavhYtWrB5w9oMNThAiKxMr8Rp06ZNNGjQgM6dO+Pp6UnFihVxcnJKUE5RFL7//vtUB/kuHh4e+Pr6Jnrso48+YubMmcmqJyYmhunTpzNr1iyuXbuGvb09DRo04Icffog3L40QKeVkDTNaaNMVjN0H+3yh2RKY0hQ83zFSv1evXgzsr60PGRAQwL///sulS5e4ceMGN2/ejPs3PFxLfoKDgwkODn5rfQ4ODrgVKkLQq9c//PADVSuWp2LFimTPnt0QtyveoV37dow8oyWix48fo2bNmuzYsYNChQqZODIhxLvolTi9vnivt7c33t7eiZZLq8QJtEcKX375ZYL9VapUSXYdAwcOZPbs2ZQqVYrPPvuMx48fs2LFCnbs2MGhQ4coVaqUASMWWY2iQPeyUCUPfLIVrj2HXuuhbwX4rjZYJ+N/o5OTE02aNKFJkybx9quqSlhYGC9fviQoKChus7S0xNbWFhsbG2xtbeNao0Ijoyn1qgXpiy++kEdvJuTm5s6NG1epW7cuu3btomTJkqYOSQiRBL1+W+7du9fQcaRatmzZ4iV0KbV3715mz55N3bp12blzJzqdDtA+7Tdp0oRBgwaxb98+A0UrsrLirrCpC0w8qE1VMP+09ujuVy8ok1O/OhVFwcbGBhsbG3Lm1LMSYRJ79uymdYtmXLhwgXr16rFjxw4ZfCNEOqZX4uTp6WnoOExu9uzZAEyYMCEuaQJtKLGXlxfbtm3j6tWrFCsmQ4hF6tlYwrj60MBDm7bg2nNoswKG1IQPK5k6OpGWcufOzb59+/Dy8uLEiRM0aNCArVu3UrNmTVOHJoRIRKZZqy48PJyFCxfy448/MmPGDM6cOZOi8729vbGzs6N27doJjnl5eQFIi5MwuAYesKMHeBWGyBiYdBA6robb/qaOLPMKiYgmJCLKwFvCJadSInv27OzevZs6deoQEBBAkyZN0mXLvhAiFWvVpTePHj2iT58+8fY1a9aMRYsW4erqmuS5wcHBPHz4kDJlyiS6Untsx/Br164lWU94eHhcB12AwMDAZEYvsjIXG23Op5UXYdx+OPEQ2qw0dVSZV5UJu0wdQqKcnJzYtm0bbdu2ZefOnbRo0YINGzbQtGlTU4cmhHhNpmhx6tevH97e3jx9+pTAwECOHDlC8+bN2bZtG61bt37nyuQBAQEAiY4MBHB0dIxX7m0mTpyIk5NT3FagQAE97kZkRYoCnUvD9u5QuwCER/137G7SP3YiGWwszani7mz061Rxd8bGMuGHr+Sys7Nj48aNtG7dmrCwMN5//3127UqfiZ4QWVWmaHEaNWpUvNfVq1fnn3/+wdPTEx8fH7Zs2ULLli2NHsewYcP46quv4l4HBgZK8iRSJL8jLG4L807ChFXavvdXwLd1oHf5+JNmiuRTFIVVA2sSGpm6R2rvYmNpjqKk7ptkbW3NqlWr6NChA5s2baJVq1b8888/NGrUyEBRCiFSI1O0OCXGzMyMvn37AnDw4MEky8a2NL2tRSn2kdvbWqRi6XQ6HB0d421CpJSZAt3K/vc6NBLG7IOOq+D6c9PFldEpioKtlYVRt9QmTbGsrKxYtWoV7733HmFhYbRq1Yo9e/a8+0QhhNFl2sQJiOvbFBISkmQ5Ozs78uTJw61bt4iOTviJNLZvk0yCKUzh+3pgZwnHH0LzpTDtXzByw4lIB3Q6HatXr6Zly5aEhoby3nvvSfIkRDqQqROno0ePAtrM4u/i6elJcHBwoq1T27dvjysjRFrrVlYbeefpDhHR8NNhaL0cTj8ydWTC2HQ6HWvWrKFFixZxydOBAwdMHZYQWVqGT5wuXrzIixcvEuz38fFh6tSp6HQ62rVrF7ffz8+Py5cv4+fnF6/8hx9+CMDIkSOJiIiI27979262b99OvXr1ZA4nYTL5HWHh+/BLU8hmDRf9tHmfRnlDYPg7TxcZWGzy5OXlRWhoKC1btuTff/81dVhCZFkZPnFauXIlefPmpVWrVnz22Wd8/fXXNGvWjHr16hEZGcm0adNwc3OLKz9t2jRKlizJtGnT4tXToEED+vfvz4EDB6hYsSJDhw6ld+/etGzZEkdHR2bMmJHWtyZEPIoC7UrC7p7QtgSoaDOPN14Em6/BOwaPigzM2tqatWvXUr9+fV6+fEmzZs1SPFedEMIwMnzi1KBBA1q1asXly5dZuHAhv//+OxcuXKBz584cOnSI/v37J7uuv/76i99//x1FUfj999/ZvHkzrVq14t9//5V16kS64WqrLc+ypC14OMHjYPh4C/TdCL4vTB2dMBZbW1s2btxIjRo18Pf3p0mTJly6dMnUYQmR5WT46Qg8PT1T1PdozJgxb13TzszMjM8++4zPPvvMQNEJYTx13GB7D5h+DKYfh7234dBdGFgZPq6avEWDRcbi4ODA1q1badiwIadOnaJx48bs37+fwoULmzo0IbKMDN/iJERWZm0BX9XUEqh6bhAeDb/9C40Wwc6b8vguM8qWLRs7duygdOnSPHjwgEaNGnH37l1ThyVEliGJkxCZQGFn+LsNzGgBee3hXiD03wR9NsjcT5mRq6sru3btokiRIvj6+tKkSROePHli6rCEyBIkcRIik1AUaFEUdveCT6qApRl4+4LXEhi3DwJk9F2mkjt3bnbt2kWBAgW4cuUKXl5eiY4wFkIYliROQmQytpYwtDbs7AFNCkFUDMw9DfUXwuKzEB1j6giFobi7u7Nr1y5y5szJ6dOnadmyJcHBwaYOS4hMTRInITKpgs4wpxUsagNFXeB5KIzYq80+vt/X1NEJQylWrBg7duwgW7ZsHDp0iLZt2xIeLs2LQhiLJE5CZHL13GFrNxjjqU2eeeUZfPSPqaMShlS+fHm2bt2KnZ0dO3fupEuXLkRFRZk6LCEyJRmwLEQ6FxJhmIXpOpeCZoVg5glYcjaa2D+rw3fD0DqQ18Fw8RgqZqOLiYF79+D6dW27dg0eP05YzsEBihT5bytYEKytDRKCod6rcpWqsGrtBtq1a8eGzVvp1a8/i+bPxdzc3CD1CyE0kjgJkc5VmbDLqPWvuwxbbkDv8vBxFXC2MW08RhUTA2fOwPbt2nb4MOjzWEtRoFQpaNoUmjWDunXB5h1v3FsY+v3M9elSAA4BVb5dyonJ3TEzk4cLQhiKJE5CpEM2luZUcXfmuK+/Ua9TIo8zDnnN+fcBzDoJy85rE2j2q6h1Mk9tPFXcnbGxNHGLR0wM7N4NixZpydKbw/YtLbUWpCJFoGhRyJsX3kw0nj//r1Xq+nV4+RIuXNC2X37RWp88PaFDB+jSBeztkwwprb6//hYuDB0+kp8n/WjU6wiRlUjiJEQ6pCgKqwbWJDTSuI+8tKRGwdsXJh/UFg/+6TAsOAOfVIVuZUBnoX88NpbmKIpinODf5elTWLAA/voLbtz4b7+dHTRsCF5e0KiRljBZpOBXoapqj/P27/+v5er+/f++/uor6NkTBg6EsmUTrcLY39+QiOi4lqypU6eS3cmBYcOGGeVaQmQ1kjgJkU4pioKtVdr8F23gAZ7usPEKTDkCdwJgzD746wR8Vg06lgIr87SLJ1UuXoQff4RVqyAiQtvn6KglMx06QK1aYGWlf/2KArlzQ6dO2qaq2jX/+QfmzNFapKZP17ZateCbb+D997Xz4lWTdu/n8OHDsbe3l+WkhDAAefAthADATIE2JWB3T/ixIeS2h4dBMHwPNPgbVl4AIzeApc7Nm9CrF5QpA0uWaElTlSpaMvPgAUybBvXrpy5pSoyiQOnS8O23cOUK7NqlJWgWFnDoELRtCzVqaPtNsAbOd99pLU2ff/45c+fOTfPrC5HZSOIkhIjHyhy6l4V9vbUpDHLYaku4fLNLS6CWnYd0NWjuwQP4+GMoXlzrx6SqWrJy7Ji2ffCB9nguLZiZaY//Vq2CO3dg2DCwtYV//4UmTbRjhw+nTSyvjBgxnCFDhgAwYMAAli5dmqbXFyKzkcRJCJEoawvoWwEO9IFhdSC7DdwNhO92Q4OFsOSciROoyEiYPFnr0D1jBkRFaaPc/v0X1q7VWptMKU8e7ZHhzZvwxRdaS9fevdrju27d4NGjNAlDURR++uknBg0ahKqq9OrVi7Vr16bJtYXIjCRxEkIkycZSG2nn0xdG1n3VAvVSe4RXbwHMPQUhkWkc1OHDULmy9ngsJERLRry9tc7ZVaumcTDvkCsX/PqrNkfUBx9orVLLlkHJkjBrljbqz8gURWHatGn06dOH6OhounTpwpYtW4x+XSEyI0mchBDJYmsJAyppCdQYT8hlp/WBGrcfas2D345CQJiRg3jxAgYNgtq14dw5yJ5dGznn46NNB5Ceublp/a3+/RcqVdLu5aOPtDmgzp0z+uXNzMyYM2cOnTt3JjIyknbt2rFrVwaek0sIE5HESQiRIq8/wpvYENydwD8Mph6BmvNgwgF48NIIF965U5t0cuZMrR9Tnz5w+TL07p1gxFq6VrkyHD2qtULZ22sdyCtVgh9+gGjjPvs0Nzdn0aJFvP/++4SHh9O6dWv27Nlj1GsKkdlI4iSE0IvOArqVhT294I9mUNIVgiNh9kmouwAGb4dLTw1wobAwGDxY67/08CEUK6b1FZo/H1xdDXABE7Cw0Po9XboEbdpo/bNGjtRG/d2+bdRLW1pasmLFClq2bEloaCjvvfce3t7eRr2mEJmJJE5CiFSxMIPWxbWFhOe3hhr5ISoG1l6GZkuh5zrY76vnSPxz56BaNa11BrTRc6dOaQlGZpA/v9aRfcECrfXJxwfKl4fFi406dYFOp2PNmjW0aNGC0NBQWrZsyb59+4x2PSEyE0mchBAGoSjQsCCsaA+bukCrYtrcUPvvQM/10GSxNpVBWNQ7q9KShj/+0Dp6nzsHOXNqE0z++ac2vD8zURTtceOZM1CzJgQGapN1du+ufW0ksclTs2bNCAkJoWXLlhw4cMBo1xMis5DESQhhcOVywbTm2lxQfSuAnSVce65NZVBjLkw+BA/f1g/q5UttvbfPP9cW4G3ZEs6e1f7NzAoV0pZxGTsWzM21kXdVq8L580a7pLW1NevWraNp06YEBwfTvHlzaXkS4h0kcRJCGI2bkzYC78gH8H1dyO+odST/8xjUng+DNsORe689lbpwQUsWVq7U+gH9+its2qQN6c8KLCxg1Cg4cEB7jHf1KlSvrs2EbiTW1tasX7+eJk2axCVPO3bsMNr1hMjoJHESQhidow76V4L9veGvllA9H0SrsOU6dF4DTZfA4nnHCKrTUFu2JF8+rfXliy8y1og5Q6lZE06ehMaNtXmqevSATz/VWuBSKCQimpCIqCQ31dyS5avX0fxVn6dWrVqxadOmJOtVVfWd9Sa4jgmWnBHC0DLAip1CiMzC3AyaFdG2y37w9xlYe1nl6jOFEVTlx++u0/aRN90G1qZ0cRdTh2taOXLAtm0wZgxMmKD17zp+XOtMnjdvsqupMiH5czVVajmcttbWrFu7lnbt2rF06VI6duyYoJyqqnSYeZgTvv7JrhugirszqwbWRMmKybDINKTFSQhhEiVc4ccKzzi6qQ2j1n1JoSdXCLZ2YLFHK1psc6HNCm1h4TSflTw9MTeH8eO1jvHOztr8T1WqaJNoJsHG0pwq7s4pvtzJOy9YsGgJ3bp1Iyoqii5durBo0aIE5UIjo1OcNAEc9/UnNF2vFC3Eu0mLkxDCNC5cgNatcbp5kw/s7Oj3uSeHqxZn6XnYdh1OPdK2sfvh/WLQubTW6TxLNla0bKktWNy6NVy8CPXqwezZ2ui7RCiKwqqBNZOdpIRERMe1TFlYWPD3339jbW3NvHnz6NWrF8+fP+eLL75I9NzjIxtja2We7PqFyOgkcRJCpL2NG7Xh9kFBULAgbNiAUrYstYBaBcAvRGttWn4BfANgyXltK+mqJVDvFwcXG1PfRBorXFhbo69nT+3969VLG204aZLWMvUGRVGwtdLvV7y5uTmzZ8/Gzs6OP/74gy+//JLHjx/zww8/JHjMZmtlrvd1hMiIMvyjuvv37/Prr7/StGlT3NzcsLKyInfu3LRv356jR48mux5vb28URXnrduTIESPehRBZhKrCxInabNlBQdCggfbYqWzZeMVcbeHjquDdG5a3hzbFQWcOl/xgzD6oNgc+/Ad23oQs9eTH0RHWrdNmGQf4+Wdo1QoCAgx+KTMzM3777TcmTJgAwMSJExkwYABRUcmZiEuIzCvDf0z4448/+N///kfhwoVp0qQJOXPm5Nq1a6xfv57169ezbNkyOnXqlOz6PD09qZ/IrMT58+c3YNRCZEHh4fDhh/D339rrjz/WphuwtHzrKWYK1MyvbePCYN0VWH0Rzj2B7Te0LbuN1gLVrgSUyZkFHuWZmWn9nsqW1dbr27pVW/R40yat9c6AFEVhxIgR5MyZk4EDBzJ37lyePn3KvL+NNz2CEOldhk+cqlWrxv79+6lbt268/QcOHKBRo0YMGjSI999/H51Ol6z66tevz5gxY4wQqRBZmJ8ftG2rLSlibg6//64lTingZA19ymvbZT9YfQnWXQK/UJh3WtsKO0PbEtCmBBRwNMqdpB+dOmmP71q10vqLVa8O69dDrVoGv9SAAQPIkSMHXbp0YePGjbz33ntQ+xuDX0eIjCDDP6pr165dgqQJoG7dujRo0IDnz59z7tw5E0QmhAC0hWyrV9eSJkdH2LIlxUnTm0q4wsi62sSa81ppy7vozOGGP/x8GOrMh/YrYeEZrb9UplW5svaos2JFePpUe/RppMky27Rpw44dO3BycuLoUem6ILKuDJ84JcXy1SMAC4vkN6xdu3aN33//nUmTJrFs2TL8/PyMFZ4Qmd+uXdpkjjdvao+RDh+Gpk0NVr2lOTQqpC3vcmIA/NwEauUHBTj+EEZ5a/2heq7TOpsHpHz+yPQvf35tpvE2bSAiQpssc9QooywSXK9ePQ4dOoS7u0fcvj179hj8OkKkZ5k2cbpz5w67du0id+7clH2j42lSli5dyhdffMGwYcPo1q0bbm5u/PTTT8k6Nzw8nMDAwHibEFnWrFnQrJnWcbl2bW0OolKljHY5Bx10LAXL2mstUSPrQvlc2gzl++/AN7ug8izos0HrJ5Wpkig7O1izBr79Vns9frw2ajEszOCXKlWqFPv2ece9btu2LTNnzjT4dYRIrzJl4hQZGUnPnj0JDw9n8uTJmCcyVPdNOXLk4KeffuLSpUsEBwdz//59Fi9ejIuLC0OHDuWvv/56Zx0TJ07EyckpbitQoIAhbkeIjCU6Gr7+Gj76SPu6Rw/YvVubCTuN5LaHAZVgYxdtoeEhNaFYdoiMgb23YcjO/5KoFRfgeWiahWY8Zmba1ARz52pr3i1bBg0bwpMnBr+Uq6tr3NfR0dEMGjSIzz//nMjIrDxbqcgqMl3iFBMTQ79+/di/fz8DBgyg51smiHtT6dKl+frrrylRogS2trbkzZuX7t27s23bNqysrBg9ejQxMTFJ1jFs2DACAgLitrt37xriloTIOIKDoX17mDJFez1unDaKLpmDM4zBIxt8Xg129oBdPeGrGvGTqKG7oPJs6LIGFpyBhy9NFqph9OsHO3ZAtmzao9Hq1bVJM41k1KjRgDbCuWHDhjx8+NBo1xIiPchUiZOqqgwYMIDFixfTo0cPgzQflylThurVq/P48WOuX7+eZFmdToejo2O8TYgs4/59bUbrDRu0RGnZMvj++3Q1P0BRF/iievwkqnQOiFHh8D0Y7Q015kGrZfD7v3DpqVG6ChlfgwZw5Ig26u72ba2f2Y4dRrnU0KHfsH79ehwdHfHx8aFixYrs37/fKNcSIj3INIlTTEwMH3zwAfPmzaNr164sWLAAMzPD3F5ss3RISGYeniNEKhw/DlWrwsmT2iO5vXuhSxdTR5Wk2CRqSzc40EfrE1U5j9ax/OwTmHIYmi2FOgu0STd97kBERppss3hxLXmqWxcCA6FFC5g+3SiXev/99zl+/Dhly5bl8ePHNGzYkKlTp6JmyKxTiKRlisQpJiaG/v37M3/+fDp37syiRYuS1a8pOaKiojh58iSKouDm5maQOoXIVNas0VqaHj7UOn8fPaq1cGQgbk5an6i1neBYf5jUCBoV1KY4uBcI809D93VQ4S/46B9thN7TYFNHnQyurrBzJ/TurfU3++QT+PxzMMLs30WLFuXw4cN0796d6OhohgwZQrt27WRkssh0MnziFNvSNH/+fDp27MjixYuTTJr8/Py4fPlygv/Mhw8fTvDpKCoqim+++QZfX1+8vLxwcXExyj0IkSGpKvzwA3ToAKGh2gi6Q4cMPnt1WsthB13LwLzWcPojmPWetj5eDlsIjoRtN7QRelXmQMtl8NMhOPYAopLuAmk6Oh3Mn68tdQPwxx/w3ntGWabFzs6ORYsWMW3aNCwtLVm/fj3lypWTKQtEppLhZw4fN24cCxYswN7enmLFisWtq/S6Nm3aUKFCBQCmTZvG2LFjGT16dLwZwrt27YqiKNSqVYt8+fLx4sUL9u/fz5UrV3Bzc5PhtllISDKexySnTKYWFgYDBsDixdrrzz/XOoSnYM60jMDWErwKa1uMCuefwO5b2nbuifb6/BOYdgwcdVDXTdvquUG+9NTFUVHgu++gWDFtkeDt27UZxjdu1PpBvSa1P/+KovDJJ59Qq1YtunXrxuXLl2ndujVuX61Jdv36srE0T7AIsUifVFUlNIULTaaX72+G/y13+/ZtAIKCgvjhhx8SLePh4RGXOL3NoEGD2LZtG97e3vj5+WFhYUGRIkUYMWIEQ4YMwdnZ2cCRi/SqyoRdpg4hfXvwQFs+5d9/teVTpk2DgQNNHZXRmSlQLpe2Da6hParbfwe8b8M+X21eqM3XtA205V/qukEdN6ieT0usTK5dO/Dw0JZpuXgRqlWDlSuhrmdcEUP9/FesWJETJ07w9ddfM3POPIPXn5gq7s6sGlgzXfxxFW+nqiodZh7mhK9/is5LL99fRZXee0YTGBiIk5MTAQEBMsIunVNVlY4zD3M8g/5HTjNHj2pJ08OH4OwMq1ZBo0amjsrkomLg9COtA/n+O9rX0a/9ZjV/lXTVLqAlUhVzg7UpP7a+kfyqU6bS0aqK0X7+N27cxCfrrmOeq1hqok6Wi+O8sLXK8G0CmVpIRBSlRm3X69zkfn+N+fdXEicjksQpY8nITcdpYuFC+PBDbVmP0qW1aQfeeMwjNAHhcPguHLgDB+/CrRfxj+vMoVIeqJlfWyKmfG6wMsx4luQLC9O+n4sWAaD260for3+AzirZVaTk5//hw4d8/PlgNmxYD0CRIkWZ9scf1KlbJ8WhvykkIjquJUsSp/Tv9cTp+MjG2L7jh1+f768x//7KT5cQryiKIr9wExMVBUOHwi+/aK/ff1/7Y+vgYNq40jEnHTQrom0A9wO1BMrnLhy6C09DtHmjDt+DqWitT5Vya4/0qudPoxYpa2stGa5QAb75BmXePGwvXYLVqyFvXoNfLk+ePKxbtZx169bxySefcO3SebwaN+DDDz9k4sSJMvgmi7K1Ms9wv3cz/Kg6IYQRPXqkPYqLTZpGjYK1ayVpSqF8jtCpNPzeTJvuYHdPmNAA3isK2W0gLAoO3YNfjmozmJeZAe1WwiQf2H0TXhh+yTmNosBXX8GWLf/NNF6pEhhxAsu2bdty8eJFPvzwQwBmzZpF0aJFmTZtmizZIjIESZyEEInz8fnvj6iDgzZf09ix2ppoQm+KAkVcoGc5+LMFnBigzWT+QwNoXQxy2WnLwZx4CDNOQL9NUP4vaLIYhu2GVRfhpr+BZzT38oJjx6BsWXj8WFvjbsoUo02bni1bNv766y+8vb0pU6YMz58/57PPPqN8+fJs365f3xch0or8BhRCxKeq8Ntv2rIdDx9q/ZmOHdNGZAmDUxRt7bwe5eCP5nD0A21h4p+baPNHFX41oPfqM1h6Hr7eCQ3+hoqzoN9GmPav1in9ZXgqAylSRGtx6tHjv4WaO3WCl8ZbvM/T05NTp04xY8YMXF1duXTpEs2aNaNly5acPn3aaNcVIjUkcRJC/CcgQFsq5csvtb5NXbtqy3YUL27qyLIMRdEWJu5YCiY3hj294OQAbSLOjypDlTxa53L/MG0+qZ8Oa7Oal52ptUp9sxOWntPml0rhWAews9MWZZ42DSwttf5OVavCmTPGuFUALCwsGDhwINeuXWPw4MFYWFiwZcsWKlasSIcOHTh//rzRri2EPiRxEkJojhzROgqvXKlNZPnbb7BkCdjbmzqyLC+7rTYJ5/A6sKYTnB8E6zvD93WhVTHI7wAqWqvUyoswbI82q3npGdB2pbbW3tpL2vHod81wrija0iz790O+fHDlClSvDn/+adQVj7Nly8bUqVO5cOFC3ITEa9asoVy5cnTp0oVLly4Z7dpCpETG6souhDC8mBiYPBlGjtQe0RQsCMuWaX8sRbpkZa6NvKuY+799T4K1+aNOP4Izj+HsYwiMgJMPtS2WrSWUcoWyuaBMDiiTU3scaPnmiPAaNeD0aejbF/75Bz79VFv3bu5cyJ7daPdWrFgxli5dyvDhwxk7diyrV69mxYoVrFixglatWjFkyBDq1auXdaYBEemOJE5CZGUPH0KvXrDr1WzOXbrAzJng5GTauESK5bSDpoW1DbQlYm6/0JKoM4/h3GO48BRCIuH4Q22LpTOH4tmhdE4tqSqZA0q6gr2rq7Ysy++/a1NSbNgAJ05oS+14eiYah6GUKVOGVatWcfr0acaOHcv69evZtGkTmzZtonLlygwZMoQOHTpgaWlp1DiEeJMkTkJkRaoKy5drj2T8/cHWVlv8tW9f7VGNyPDMFCjkrG1tS2j7omPghr/W/+ncEy2RuvAUgiLg7BNte527E5R0VShR/QtK/NOKEiMG4nZ8N+YNGsDgwTBhAtjYGPU+KlSowLp167hy5Qq//vorCxYs4MSJE3Tr1o2vv/6avn370q1XX6PGIMTrZOZwI5KZw0W69OQJfPyxNr0AaFMOLFkCJUqYNi5hEjEq3Al4lUQ9gYt+cMkPHgUlXt46Jpxi985S/OF5ikU/oXi/9yleuwS57NIm5/bz82PGjBlMmzaNJ0+0TM/MypoCg1cDcHpkQ7LZGzeZE6nz+szhyZkJPKXlQWYOF0IYypo12oK8fn5aB/Dvv4dhw7QRVCJLMns1is8jG7Qs+t/+56Fw8Slc9nu1PdM6l4eh46xbVc66VdUKXtY2RyuVItkViroQtxVxgbwO2jUMxdXVle+//55vv/2WDRs2MHv2bHZ5/zdhZ8GCBWn//nt06dKF+vXrY2Ehf+aEYclPlBBZwb178MUX2qzfoE10uHAhVKxo2rhEuuVioy1IXMftv33RMdq6e1efwdV7IVzxPsfVCAdu5ShGYIRFgo7oADYW2uPCws5aIlX41ePDgtnAJhX5upWVFR07dqRjx45cvHqDFvMuAxAQ8IK5c+cyd+5ccubMSceOHWnbti316tWT/lDCICRxEiIzi4rSOvaOGgXBwWBuDt9+q73W6UwdnchgzM205KeIC7QoagsNqsO6dYR/4sVt1YlruUtxtWlXrlVvzrVga26/gNCo//pSvSmv/ask6lUiFbvld0xklF8SPDzc0Zq+YNvWbaxdtZzVq1fz5MkT/vzzT/78808cHR1p1qwZrVq1onnz5mQ34shAkblJ4iREZnX4sPZY7uxZ7XWtWjBjBpQrZ9q4RObSti26+vUpPnw4xf/6C06v0ta9mzSJyIEDuBtkxg1/uP4cbjyHGy+0JWNehMGDIG3zuRu/SgszbW4q92yvHiM6aR3V3bJBAcekF0CuVL0mderW4X9TfsHb25u1a9ayddtW/Pz8WLVuA6vWbcDMzIyKFSrSqH49GjduRO3atbExcCd3VVUJTeEMpDaW5ulmmoWMHr8xSedwI5LO4cIkfH1hxAitwzeAi4s2T1PfvrLOnDCuo0dh0CA4dUp7XakS/PSTtvbdG56HaiP8bvpr0ybcegG3/OF2gLbocVJy24ObI7g5aYlUTtsohi9N+Rp3Yfcu8HjJt+h0OmrVqkXt2rWpXbs2NWrUIFu2bCmuL5aqqnSYeZgTvv4pOq+KuzOrBtY0efJh7Pilc7gQIn148QJ+/FF7NBf+auGyvn21pMnV1aShiSyienX491+YPl0beHDyJDRqBC1baj+HpUrFFXWx0baqeeNXEaNqI/p8A8D3hfbvrRfayL87AfAyQjv+KAj+ffDqJNUcLJ1RIlP2h946f2nyunnw4M5t9u7dy969ewFQFIXSpUtTs2ZNqlSpQuXKlSlTpgy6ZD7eDo2MTnHSAXDc15/QyOhkJQbGlNHjN7bMfXdCZAVhYfDXXzBuHDx/ru1r0AB+/ln7xC9EWrKwgM8/19Y5HDdOm1B182bYuhU++EDrX5c//1tPN1O0kXh5HaDmG8VUVVuj706AllDdCYC7gXA3QMHXsSYPA6OJecczFDMFctpG8/imNunrwIVXsY54xuPrp7h59hBnD27jxvnjnD9/nvPnzzN79mwALC0tKVu2LBUrVqRs2bJxW44cOZK83vGRjbG1SrrDVkhENFUm7Eo6cBPJ6PEbgyROQmRUwcFawvTTT/DokbavVCntk32LFjKRpTCtHDm0SVU/+wy++w7WrYPZs7XRnP36afvc3VNUpaL811JVIXeCo0RGW/AoGO4FvrG9hAcv4eFLiIyBR8EQ+79jzikFzHIDzaFkcyg5nqLmKg5KMGbBjwl7egs/3wsEPbrOpRf3OL/zLNGrtxH98hHERJMrVy5Kly5N8eLFKVGiBMWLF8e9cLG4qGytzDN0C0xGj98Y5N0QIqN5+VJ7FDJlCjx9NVTJzU1ba65vX+0TvxDpRbFi2jQYPj5a37v9+7VWqDlztOV+hg2DIkUMcilLc63PU4G3dGmJjoGnIXD9GfSYoe3rUQ6ehsL9l/DoJfiFQni0Qjj2oLOH/IWxyt8YlzcrU2OIDnxMdOADzgU+5PTzh0RvfUD0inVEBz8lZ7N+ALTv2JXiBfNRuHBhChYsiIeHB+7u7jg4OBjknkXak9+wQmQUN29qK9TPm6f1ZwIoVAiGD4eePcHKyqThCZGkOnVg3z4tcRo/Xlsfcd48WLAA3n9fa5mqX9+oLaXmZlrHcsfX/quMqAu2r70Oj4LHwfAwSGuhehT8X5+qR0HwOAiehEBUjBnmTnkwd8qT8EIxUfD4VWfm+ss4FxpA9P1HxFx+TPTLs0S/3Iku+iXZrWPIaW9OHhdrsND6fx30OUgh9/zkzZvX4CP9hGFI4iREeqaq2h+YP/7QVqiPHQRbrJj26b1bN2lhEhlLvXqwc6c2XcaECbBli/YYb906KFNGS6C6dwc7O5OEp7PQRuu5JbHOdYwKz0K0BCs2kXry2tcPA1QuPP6vvLlddsztskOe0vHqCXq13YyJQnmVaPU7Wobo3c+IDjqBefgLbAnBwSISZ10MOezNye1kRYHsdrjncqBwHmfy58lB9uzZMTdPwcRXIlXkN64Q6dHt29oK9H//Ddeu/bffy0vreNusmUwtIDK2mjW1TuMXLsC0adrP+vnz8NFHMHQodO6sPcqrVSvd9dczUyCHnbaVyZnweEiEQqnz2teH+0FQlJZYPQ0BvxC47x+Or18IDwIieR6mEBhjTnhs3dZOmNlmxzKn1k8qEnj+arsRe4EXr7YrEBMeREzwXczCA7CMDsJaDcPePAJHXQwuNgqudhbkcrIir4stBVztccvpiHtuF+zsbI32/mR2kjgJkV4EBGh9QRYu1B5pxLK31/ouffIJFC9uuviEMIbSpbWJWSdOhPnztSTq5k2YNUvbChfWEqju3bWvM5hsNpDXCorFm6hc92rTaPMUaV9v6gohUSr3noVw81Egd5+F8CQwkmchKi8iLQiO0RFuZke0lSOYWWKms8dMZw9ANBD8aotr8FL5L9G6+WpXTDRqqB9KeCAWUUFYxYRirYRjZxGFo6WKve6/oYkzNp6kgKsteZztyJ/DkbyujuissvbSNZI4CWFKDx/Chg2wfj3s2QORkdp+RdEmDezVC9q105InITKzbNlg8GBtTUVvb60FavVquHEDRo/WtnLloG1bbStXLt21RBlCERewtVIgvx2Uf/vjSlXV5rR6GhTFrUcB3H4cwIPnoTwKCMcvKBr/UJXASHOCY6wIV2yINLcnxsoRxcoWxcwcxc4V7FyJBkJfbXEzN8VEoaA9Ovz9biW4/1+qoMbEoIa/QIl4iXlUMJYxoVjGhMcdbzd5P07WCk7WZmSztcDFzhJXRx05nWzImc2G3M72ONhbG/x9S0uSOAmRliIj4dgxrd/S1q1w5Ej84yVLah29u3fXRsoJkdWYmWkfGho21AZDrFuntcLu3astH3T2LIwdCwULwnvvaRNs1q8PTkl0SsqEFAUcdeCos6Bw9uxQOnlr74VGqjx8Hozv4wDuPwvioX8ojwPCeRYcxYvQGALCzXgZqRC7tKDy4hYxVo6gc0SxtEExM0OxyQY22YgBwoHw1/poXbKrB2YWWvPXy1fbozeCiImKmw6i5I93MI8JxzImHEsi0CmRWJtHY2seg50l2OvAxuq/BHnO1jPkcrIhu4M1rk42uDrZ4upkh3UaTpkgiZMQxhQRAWfOwMGDsHu39kk6KCh+merVtU/QbdrIozghXmdnBz16aNuzZ9oAiXXrYPt2uHVLGzTxxx9aslWtmpZE1aunfZ2KJVMyMxtLhUK57CmU6+2t2PGWOBlRMG4ep+DwKB76veSe30se+Qfz6EUoz15G8ORlJLEL3rgHHiQ8xoKwGAvCsSJS0RFtbkuMhS3oHFAs3xgp6ORGtJkF0UAYWp6VQMR/LWBTbpbXErM3qJGhEBmCEhmKWXQYSqhfyt6YFMg0idOxY8cYPXo0hw8fJiIigtKlS/Pll1/SrVu3ZNcRExPD9OnTmTVrFteuXcPe3p4GDRrwww8/ULRoUSNGLzKF6GjtscLp09qyE0eOwIkT2szer8ueXfs03bix9ok5b95EqxNCvCZ7dujdW9uCg7WReTt3aq23V69q/9+OHIEfftCaY0qWhBo1tA8mlSppfalkeH+q2OksKJLPmSL5nOPtfz3R2jasdpITZgaHRXL7SQAtp2uvvyt2npchUbwIieRFSBQvw2J4GaESHKkQGqUQFmNBWLQ5sR83lRe3iLGwA0tbsLJDMdf6WymWNvAqKYsBYsIS6bVvIJkicfL29sbLywsrKyu6dOmCk5MTa9eupXv37ty+fZvhw4cnq56BAwcye/ZsSpUqxWeffcbjx49ZsWIFO3bs4NChQ5R6bZ0lkYVFRWmj3q5e1bbz57XHB+fPQ2howvIuLtov8Pr1tWSpfHkZESdEatjZaS20bdpor+/e1Vp0d+/Wpjm4cQMuXtS2efO0MmZmULSo1jeqTHmggrY/OASsZBH2tGJnbUnB3NniXvduUiZli/y+1gIWE6MSHBbOkxfBPA8MxS8wlOcvwwkIieDRU3/GG+keMnziFBUVRf/+/VEUhf3791OxYkUARo8eTc2aNRk9ejQdO3Z8Z4vR3r17mT17NnXr1mXnzp1xizn26tWLJk2aMGjQIPa9PtJJZF5BQfD4Mdy/D3fuxN9u3NBG/ES9Zfl2Gxvtk23VqlqyVLOmNityJuzEKkS6UaAA9OmjbaDNqB/bAnX0qPa43M8PrlzRtvUb4as1WtlcOSGnq/b/1N1d61sYu+XPD7lzax9+5P9wumNmpuBgq8PBVkfhNxruAwMDGf+Bca6b4ROnPXv2cOPGDfr27RuXNAE4ODjw/fff06VLF+bPn8+PP/6YZD2xCzlOmDAh3grYjRo1wsvLi23btnH16lWKFSv2tipEeqGqWt+i4GBtCwzUhvrH/hsQoC2G++yZtj1/rv1SffRI24KD330Na2vt02vRotojgfLlta1wYZCJ6IQwrRw5oFUrbQPtd8Ljx1oCdfYsnL8Yv/z9+9r2NpaWkCuXlkTlyKE9NnRx0f7Nnl3rT+XoqHVQd3QE29eWU1HVt1YrMqYMnzh5e3sD0LRp0wTHYvclp6XI29sbOzs7ateuneBYbOK0b98+/RKnjRvB1giTjSX3P2Ri5WL3vX7sbfte3//m12/bYmL++/fNLTo6/r9RUdrXsVtUlDb6LDIy/tfh4VpCFLuFhWlbaOh//4aEaIlPdHTK3ss32dpqfY9e//Tp5gYeHtqs3fnyyeM2ITIKRdGSnty5tUlkI6Lg1aOfkBu+cOcW3LwB9+7B3Xvao797d+H+A3jxapD+46falgwhltbw2RLthYODtoienZ3WIm1jo33wiv3Xygp0tlCqv3buRx+DuaIla5aWYGGurQ5gbqF9KDOPfW2m/Q4yN9f+NYv910y739e/VhRt1k4lkS32/XntdUiMGfAq+Vu1Gixea217/ZzYf6MBtP5FIWvWa/En9f7EALyakmDDRnjXZ81oiJ33KmTdRkjGr96QkJB3F9JThk+crr2aVTmxR3HOzs64urrGlXmb4OBgHj58SJkyZRKdtj627nfVEx4eTnj4f/NZBAQEABDYs2fSNyGMx8JC+8Xl5KT96+iobc7O2idGFxft6+zZIWdObcuV693zJr05Mk4IkWGEREQRE679Ya00df+rvVZAIdAVgiJoW2q8qj9QVYkKCdE+1L0tHgsrYgprA5kqZWsS/6CKNn14ZCrjSbFX8ffpTVRURNIlLayI+XQxAJX+BS3oZNbfrUfK6j+SZNE4MeHah2fVCC1+GT5xik1OnN4yh4ejoyP37t1LdR2vl3ubiRMnMnbs2AT7CyR5ljCqqCjw99c2IYRIY4ksAZxQVAT82snYoejF2PEbu/5nz5699W+7vjJ84pSeDBs2jK+++iru9YsXL3B3d+fOnTsG/8aJdwsMDKRAgQLcvXs3LvkVaUfef9OS99+05P03rYCAANzc3HBxcTF43Rk+cYpNSN7WGhQYGPjOpCU5dbxe7m10Ol28juWv1y//cUzH0dFR3n8TkvfftOT9Ny15/03LzAh9UTN879ak+h/5+/vj5+f3zqkI7OzsyJMnD7du3SI6kU7FSfWjEkIIIUTWkeETJ09PTwB27NiR4Fjsvtgy76onODiYgwcPJji2ffv2ZNcjhBBCiMwrwydOjRo1olChQixdupTTp0/H7X/58iXjx4/HwsKCPrGTogF+fn5cvnwZP7/469h8+OGHAIwcOZKIiP96+O/evZvt27dTr169FE9FoNPpGD16dKKP74TxyftvWvL+m5a8/6Yl779pGfP9V1RjjNVLY3v37sXLywudTkfXrl1xdHRk7dq13Lp1iwkTJjBixIi4smPGjGHs2LGMHj2aMWPGxKtnwIABzJkzh1KlStGyZcu4JVesra1lyRUhhBBCZPwWJ4AGDRrg4+NDnf+3d+9BUVb/H8Dfq+Iiy4IIDBcvoMiowziz3hM0AhWMBmSi8oYjeUNnbByMUCopJsW0mUpLUjMuhVaWTDURXhi5hPoli9qhEZIUU6uZMnQJCgP28/3D2B8IysPvi/vg7vs1wwxzzgHe+wD7fPY8z54zaxYOHz6MzMxMuLu7Iy8vr1PR1JN9+/Zh9+7d0Gg02L17NwoKChAdHY2vvvqKRRMRERHZxowTERERkTXYxIwTERERkTWwcCIiIiJSiIXTPXD27FlERUXBzc0NOp0O06dPx6FDh9SOZXd27twJjUYDjUaD//xH4QZH9D8REeTn5yMsLAw+Pj5wcnLCuHHjkJiYiIsXL6odzybk5eUhMTERU6dOhVarhUajQU5OTpdxLS0tOHLkCBISEjBhwgTodDro9XrMmDEDmZmZ3a5ZRz1Tevw7qqurw+rVq+Hn5wetVgsvLy+EhYXho48+sk5oG/Lzzz/j9ddfR0REBEaNGoXBgwfD29sbcXFxqKio6PZrGhoasHHjRsvx9/Pzw8aNGy2LW/eaUJ8qLi6WwYMHi7Ozs6xatUqefvppGT16tACQbdu2qR3Pbpw7d060Wq3odDoBIGfOnFE7kl3YuHGjABAfHx9Zu3atpKSkSGRkpGg0GtHr9VJVVaV2xPuen5+fABAPDw/L59nZ2V3GVVdXCwDR6/WyYMECSUlJkcTERPH19RUAEh0dLWaz2foP4D6n9Pi3O378uDg5OYmTk5MsXLhQUlNTZe3atRIcHCxr1qyxXnAbsWnTJgEgAQEBsmLFCtm8ebPExcXJwIEDZcCAAfLhhx92Gt/Y2CgGg0EAyLx582TTpk0yf/58ASAGg0EaGxt7nYGFUx9qaWmRgIAA0Wq1UllZaWlvaGiQoKAgGTRokJw/f17FhPahtbVVpk2bJtOnT5f4+HgWTlby66+/yoABA8Tf319MJlOnvtdee00AyJNPPqlSOttx4sQJuXTpkoiIbN++/Y4n7qtXr0pmZqY0NTV1am9sbJSpU6cKADl8+LA1ItsUpcdfROTy5cvi4uIigYGB8tNPP3Xpb2lpuZdRbdKRI0ekrKysS3tZWZk4ODjIsGHDpLm52dKelpYmACQlJaXT+Pb2tLS0Xmfgpbo+dPLkSVy4cAFLlizBpEmTLO16vR5btmxBa2srsrOzVUxoH3bs2AGj0YisrCwMHDhQ7Th249KlSzCbzQgJCemyN9cjjzwCAPjtt9/UiGZT5s6dCz8/vx7HDR8+HOvWrYOTk1Ondp1OZ9mMvLS09J5ktGVKjz8AZGRkoKGhAXv37sWoUaO69A8adN9vF2t1jz76KGbPnt2lffbs2QgLC0N9fT2qqqoA3Lp14MCBA3B2dkZaWlqn8ampqXBzc8M777wD6eXiAvyt9aGSkhIAQERERJe+9jY+Ud1b33//PdLT0/H8888jKChI7Th2JTAwEIMHD8apU6fw559/Qq/XW/q++OILAEB4eLha8agDBwcHADxx30sigsOHD8Pd3R3h4eH45ptvUFpaCrPZDIPBgPDw8HuyAa09u/3vura2Fr/88gsiIyOh0+k6jXV0dMSDDz6ITz/9FD/++GOv9qLlf00futtmwG5ubvDw8Oh2M2LqG62trZYbYTdv3qx2HLvj7u6Obdu24ZlnnsGECRMQExMDvV6PqqoqFBUVYc2aNXjqqafUjkkAsrKyAHT/Io/6Rl1dHerr6zFt2jSsW7cOe/fu7dQ/adIkfPbZZxgxYoRKCW3L5cuXUVRUBG9vb0ycOBHA3c/JHdtra2tZOKnFZDIBAFxdXbvtd3FxwdWrV60Zya5kZGTAaDSioqLC8sqDrCs5ORm+vr5ITEzEW2+9ZWkPDg5GfHw8fy/9wP79+1FYWIjw8HBERUWpHcdmtV+WrqysRHV1NbKzs7FgwQKYTCZkZGTg7bffxmOPPcZ3/PaBlpYWLFu2DDdv3sTOnTstt2goOSd3HKcU5wnJJhiNRmzduhXJycmYPHmy2nHs1tatW5GQkIDU1FRcuXIFjY2NKC8vR2trK8LCwpCfn692RLtWUFCA9evXw8/PD3l5eWrHsWlmsxkA0NbWhpdeegkJCQlwc3ODv78/9u/fjxkzZqCiogLl5eUqJ72/mc1mrFixAmVlZVi9ejWWLVt2z38mC6c+1F7V3ql6bWhouGPlS/+b5cuXIyAgoMvGzWQ9J0+exJYtW7B+/Xo8++yzGDFiBHQ6HUJCQvD5559jyJAhSEpKUjum3Tp27Bji4uLg5eWFkydPwsfHR+1INq3jc31MTEyX/ujoaADA119/bbVMtkZEsHr1auTl5SE+Pr7L5VAl5+SO45Ri4dSHOl4vvd3169dx7dq1Xl1HJeWMRiNqamrg6OhoWfRSo9EgNzcXADBz5kxoNBp88skn6ga1YQUFBQBubbp9O09PT0ycOBGXL1/GtWvXrB3N7h09ehSxsbHw8PBAcXExxowZo3Ykmzd27FjLJaOhQ4d26W9v+/vvv62YynaYzWasXLkSWVlZWLx4MXJycrrcbH+3c3LH9t6el3mPUx8KDQ3F9u3bcfz4cSxatKhT3/Hjxy1jqO+tXLmy2/aysjLU1tYiJiYGnp6e8Pf3t24wO/LPP/8AAH7//fdu+9vbtVqt1TLR/xVNw4YNQ3FxMcaOHat2JLug1WoRHByML7/8EufOncOsWbM69Z87dw4A+Jz0/2A2m7Fq1SpkZ2dj4cKFeO+997pdeiYwMBC+vr44deoUmpqaOr2zrrm5GWVlZfD19e39/0SvV36iO2ppaZExY8aIVquVb7/91tLecQHMH374Qb2Admj58uVcANNK3n//fQEgQUFBcuPGjU59OTk5AkCmTJmiUjrb1NMCjIWFhaLVasXb21tqamqsG84O9HT8Dx06JABkzpw5nRZlrK6uFicnJ9Hr9VJfX2+ltLahra1NEhISBIA8/vjjPS4iei8WwNSI9HLlJ7qr4uJiREZGQqvVYvHixXBxcUF+fj7q6uqwdetWPPfcc2pHtCsJCQnIzc3FmTNn8MADD6gdx6a1tbVh7ty5KCkpgaenJ2JiYuDm5gaj0YgTJ05Aq9WiqKioyytv6p0DBw5YbiiuqqpCZWUlQkJCLK+aY2NjERsbi5qaGhgMBty8eROLFi3CuHHjunwvf39/JCQkWDP+fU/p8Qdu3YPzxBNP4OOPP8a4ceMQGRkJk8mEI0eO4K+//sK7776LpUuXqvVQ7ksvvvgi0tPT4ezsjA0bNnS7FllsbCwMBgMAoKmpCbNmzcJ3332HefPmYcqUKTAajSgsLITBYEB5eXmXNZ561OtSi3pUUVEh8+fPF1dXVxkyZIhMnTpV8vLy1I5llzjjZF3Nzc2yY8cOmTx5sjg5OcmgQYNk+PDhsmTJEu5T10fa/6bv9PHCCy+IyK19M+82DoCEhoaq+ljuR0qPf7uWlhZ59dVXJSgoSLRarbi4uEhERISUlJSo8wDucz0df3QzA3jjxg1JSkqSkSNHioODg4wcOVKSkpK6zIwrxRknIiIiIoX4rjoiIiIihVg4ERERESnEwomIiIhIIRZORERERAqxcCIiIiJSiIUTERERkUIsnIiIiIgUYuFEREREpBALJyIiIiKFWDgRERERKcTCiYiIiEghFk5ERERECrFwIiIiIlKIhRMR0b9mzpwJjUaDs2fPdmq/fv06goKC4OjoiNLSUpXSEVF/wMKJiOhfL7/8MgAgLS3N0tbc3IyYmBjU1NTg4MGDCA0NVSseEfUDLJyIiP4VGhqKhx9+GEePHsXp06dhNpuxdOlSlJeX480330RcXJzaEYlIZRoREbVDEBH1F0ajEZMmTUJ4eDjGjx+PPXv2IC0tDenp6WpHI6J+gIUTEdFt4uPjcfDgQQDAmjVrsG/fPpUTEVF/wUt1RES38fDwAAC4urrijTfeUDkNEfUnLJyIiDrYtWsXdu3aBS8vL5hMJuTl5akdiYj6EV6qIyL61wcffIAlS5Zgzpw5yM3Nxfjx4zF06FCcP38ejo6Oascjon6AM05ERACKioqwfPlyGAwG5Ofnw9fXFxs2bMCVK1ewZ88eteMRUT/BGScisnuVlZV46KGH4OnpidOnT8PLywsAYDKZMHr0aAwYMAAXL16Ei4uLykmJSG2ccSIiu3bhwgVERUXB0dERx44dsxRNwK2bw5OTk/HHH3/glVdeUTElEfUXnHEiIiIiUogzTkREREQKsXAiIiIiUoiFExEREZFCLJyIiIiIFGLhRERERKQQCyciIiIihVg4ERERESnEwomIiIhIIRZORERERAqxcCIiIiJSiIUTERERkUIsnIiIiIgU+i8st/komie+4wAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Plot fit\n", "fig, ax = plt.subplots(1,1)\n", "plt.gcf().subplots_adjust(bottom=0.15)\n", "plt.gcf().subplots_adjust(left=0.15)\n", "yMin = 0.\n", "#yMax = f(0., parhat)*numVal*binSize*1.1\n", "yMax = np.max(xHist)*1.1\n", "plt.xlim(xMin, xMax)\n", "plt.ylim(yMin, yMax)\n", "plt.xticks(np.arange(xMin, xMax+1, 4.0))\n", "xCurve = np.linspace(xMin, xMax, 100)\n", "yCurve = f(xCurve, parhat)*numVal*binSize\n", "sig_parhat = np.copy(parhat)\n", "sig_parhat[0] = 1.\n", "bkg_parhat = np.copy(parhat)\n", "bkg_parhat[0] = 0.\n", "sigCurve = f(xCurve, sig_parhat)*numVal*binSize*parhat[0]\n", "bkgCurve =f(xCurve, bkg_parhat)*numVal*binSize*(1.-parhat[0])\n", "plt.plot(xCurve, yCurve, color='black')\n", "plt.plot(xCurve, sigCurve, color='red')\n", "plt.plot(xCurve, bkgCurve, color='dodgerblue')\n", "\n", "# Plot data histogram\n", "binLo, binHi = bin_edges[:-1], bin_edges[1:]\n", "xPlot = np.array([binLo, binHi]).T.flatten()\n", "yPlot = np.array([xHist, xHist]).T.flatten()\n", "plt.plot(xPlot, yPlot)\n", "plt.xlabel(r\"$x$\")\n", "plt.ylabel(r\"number of entries\")\n", "if fitType == 'LS':\n", " plt.figtext(0.55, 0.8, r\"Least Squares\")\n", "elif fitType == 'M':\n", " plt.figtext(0.55, 0.8, r\"Maximum Likelihood\")\n", "plt.figtext(0.55, 0.72, r\"$\\hat{\\theta} = $\" + f\"{parhat[0]:.4f}\" +\n", " r\"$\\pm$\" + f\"{sigma_parhat[0]:.4f}\")\n", "plt.show()\n", "plt.savefig(\"histFit.pdf\", format='pdf')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.13" } }, "nbformat": 4, "nbformat_minor": 4 }