{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Example of fitting a peak in data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The data is binned, there is assumed to be a single peak, fit the peak and plot the fitted result to compare with data" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "# Load some of the libraries we'll need\n", "%matplotlib inline\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from scipy.optimize import curve_fit\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG0xJREFUeJzt3XmUZVV59/HvT0BBZLChwWaywaDRNyqYFiEkDhgNigIx\nmsQRCRGTmDcYfRMxyQvtkKUmjjGuRBKMKAIiGkDFKCLgMirYDDIEeAWCyiDdMjQIAWx43j/OKbi0\n1VWnqvveOnXr+1nrrnvOPtOzq6v76b33OfukqpAkqW8eMdcBSJI0GROUJKmXTFCSpF4yQUmSeskE\nJUnqJROUJKmXTFDSkCT5ZJJ3r2Pb65N8awjXPCfJH27o8w7bsH4emt9MUOqNJNcl+Z8kdya5Pcm3\nk/xRkk6/p0mWJqkkGw87VknDZ4JS37y0qrYAHg+8F3gbcOzchrTwmOTVByYo9VJVra6q04HfAw5J\n8isASQ5IclGSO5L8OMnygcO+2X7fnuRnSfZJ8oQk30hyS5KfJvlMkq3Xdd0kH2nPe0eSC5L8xsC2\n5UlOTvKptpV3eZJlA9v3THJhu+2zwKbTVDNJPppkdZIrkzx/YMOhSa5oz3VtkjeudeBBSS5u47wm\nyf6TnHxJkkuS/J92fdck32zP+fUkH0tyfLttovV5WJIfAd9oyw9s63l723345IHzV5JfGlh/sEsz\nyXOTXJ/krUlWJrkpyaED+26T5PQ2/vOBJ0zzs9ICZIJSr1XV+cD1wESiuAt4HbA1cADwx0kObrc9\nu/3euqoeU1XfAQK8B9gBeDKwM7B8ikt+D9gDWAScAHwuyWCiORA4qb3+6cA/AiR5JHAq8On22M8B\nvzNN9Z4FXAtsCxwNfCHJonbbSuAlwJbAocCHkjyjvdZewKeAv2jjeDZw3eCJkywFzgX+sare3xaf\nAJwPbNP+DF47SUzPofk5/VaSJwInAm8GFgNnAF9s69rF44CtgB2Bw4CPJXlsu+1jwD3AEuAP2o/0\nMCYozQc30vyjT1WdU1WXVtUDVXUJzT+gz1nXgVV1dVWdWVX3VtUq4IPT7H98Vd1SVWuq6gPAo4An\nDezyrao6o6rup0lGT2/L9wY2AT5cVT+vqlNokt1UVg7s/1ngKpqkS1V9uaquqca5wNd4KEkfBnyi\nrdcDVXVDVV05cN6nAOcAR1fVMQBJdgGeCRxVVfdV1bdoEuzallfVXVX1PzSt1y+31/k58H5gM+DX\npqnXhJ8D72zrdwbwM+BJSTaiSd5Htde6DDiu4zm1gJigNB/sCNwKkORZSc5OsirJauCPaFogk0qy\nXZKTktyQ5A7g+Gn2f2vbtbY6ye00LYDB/X8ysHw3sGk7XrMDcEM9fPblH05Tr8n236GN40VJvpvk\n1jaOFw/EsTNwzRTnfTVwA3DKQNkOwK1VdfdA2Y8nOXawbIfBOlTVA+32Haes1UNuqao1A+t3A4+h\naY1tvNa1pvtZaQEyQanXkjyT5h/EiVuQT6D5n//OVbUV8M803XgAk03N/562/GlVtSXwmoH9177W\nb9DclPG7wGOramtg9br2X8tNwI5JBvfdZZpjJtv/xiSPAj5P02LZvo3jjIE4fszUYzbLgZ8CJ7St\nlYn4FiV59MB+O09y7ODP8Eaam1WAZsCsPeaGtuhuYPB8j5sipkGrgDVrXX+6n5UWIBOUeinJlkle\nQjPec3xVXdpu2oKmJXBPOxbzqoHDVgEPALsNlG1B07V0e5IdacZt1mULmn84VwEbJzmKZgyoi++0\nx/5Zko2TvAzYa5pjtmv33yTJK2jGfs4AHknTtbgKWJPkRcALB447Fjg0yfOTPCLJjkl+eWD7z4FX\nAJsDn07yiKr6IbACWJ7kkUn2AV46TXwnAwe019kEeCtwL/DtdvvFwKuSbNTepLHOrtNBbffoF9pY\nHp3kKcAhXY7VwmKCUt98McmdNK2Ev6YZMzp0YPufAO9s9zmK5h9RANruq78F/rO962xv4B3AM2ha\nQl+m+YdxXb4KfAX4fzRdTvcweTfYL6iq+4CXAa8HbqMZv5nqWgDnAbvTtHb+Fnh5O/51J/Bnbd1u\no0nCD44XtTeOHAp8qK3XuQy0dNaKZzvgE2meJXs1sA9wC/Bu4LM0CWdddbqKpsX50TbGl9I8BnBf\nu8sRbdnt7blPnaa+g/6UprvvJ8AngX+bwbFaIOILC6WFqb0V/sqqOnquY5EmYwtKWiCSPDPNc2GP\naLvkDmJmrR5ppHxaXFo4HkfT7bgNzbNlf1xVF81tSNK62cUnSeolu/gkSb00L7r4tt1221q6dOlc\nhyFJ2gAuuOCCn1bV4un2mxcJaunSpaxYsWKuw5AkbQBJOs0cYhefJKmXTFCSpF4yQUmSeskEJUnq\nJROUJKmXTFCSpF4yQUmSeskEJUnqJROUJKmX5sVMEpJmaflWszhm9YaPQ5oFW1CSpF4yQUmSeskE\nJUnqJROUJKmXTFCSpF4yQUmSeskEJUnqJROUJKmXTFCSpF4yQUmSeskEJUnqpaHOxZfkOuBO4H5g\nTVUtS7II+CywFLgO+N2qum2YcUiS5p9RtKCeV1V7VNWydv1I4Kyq2h04q12XJOlh5qKL7yDguHb5\nOODgOYhBktRzw05QBXwtyQVJDm/Ltq+qmwDa7+0mOzDJ4UlWJFmxatWqIYcpSeqbYb8Pat+qujHJ\ndsCZSa7semBVHQMcA7Bs2bIaVoCSpH4aaguqqm5sv1cC/w7sBdycZAlA+71ymDFIkuanoSWoJJsn\n2WJiGXghcBlwOnBIu9shwGnDikGSNH8Ns4tve+Dfk0xc54Sq+o8k3wNOTnIY8CPgFUOMQZI0Tw0t\nQVXVtcDTJym/BXj+sK4rSRoPziQhSeolE5QkqZdMUJKkXjJBSZJ6yQQlSeolE5QkqZdMUJKkXjJB\nSZJ6yQQlSeolE5QkqZdMUJKkXjJBSZJ6yQQlSeolE5QkqZdMUJKkXjJBSZJ6yQQlSeqlad+om2Rf\nYDnw+Hb/AFVVuw03NEnSQtblle/HAn8OXADcP9xwJElqdElQq6vqK0OPRJKkAV0S1NlJ/h74AnDv\nRGFVXTi0qCRJC16XBPWs9nvZQFkB+234cCRJakyboKrqeaMIRJKkQetMUEleU1XHJ3nLZNur6oPD\nC0uStNBN1YLavP3eYhSBSJI0aJ0Jqqo+3n6/Y3ThSJLUmHYmiSS7JfliklVJViY5LYkP6UqShqrL\nVEcnACcDS4AdgM8BJw4zKEmSuiSoVNWnq2pN+zme5jZzSZKGZqq7+Ba1i2cnORI4iSYx/R7w5RHE\nJklawKa6i+8CmoSUdv2NA9sKeFeXCyTZCFgB3FBVL0myK02yWwRcCLy2qu6baeCSpPG2zi6+qtq1\nqnZrv9f+zOQmiSOAKwbW3wd8qKp2B24DDptd6JKkcTbU90El2Qk4APjXdj00UySd0u5yHHDwMGOQ\nJM1PXebiWx8fBv6Shx723Qa4varWtOvXAztOdmCSw4HDAXbZZZchhyn13PKt5joCaeSmbEGlsfNs\nTpzkJcDKqrpgsHiSXSe9I7CqjqmqZVW1bPHixbMJQZI0j03ZgqqqSnIq8KuzOPe+wIFJXgxsCmxJ\n06LaOsnGbStqJ+DGWZxbkjTmuoxBfTfJM2d64qp6e1XtVFVLgd8HvlFVrwbOBl7e7nYIcNpMzy1J\nGn9dEtTzaJLUNUkuSXJpkkvW45pvA96S5GqaMalj1+NckqQx1eUmiRet70Wq6hzgnHb5WmCv9T2n\nJGm8TduCqqofAjsD+7XLd3c5TpKk9dFlNvOjabrl3t4WbQIcP8ygJEnq0hL6beBA4C6AqroRX2Io\nSRqyLgnqvqoq2ueVkmw+zf6SJK23Lgnq5CQfp3l+6Q3A14F/GW5YkqSFbtq7+Krq/UleANwBPBE4\nqqrOHHpkkqQFretcfJcCm9F08106vHAkSWp0uYvvD4HzgZfRzADx3SR/MOzAJEkLW5cW1F8Ae1bV\nLQBJtgG+DXximIFJkha2LjdJXA/cObB+J/Dj4YQjSVKjSwvqBuC8JKfRjEEdBJyf5C0AVfXBIcYn\nSVqguiSoa9rPhInZx31YV5I0NF1uM3/HKAKRJGmQk75KknrJBCVJ6iUTlCSpl7o8qPt3SbZMskmS\ns5L8NMlrRhGcJGnh6tKCemFV3QG8hOaZqCfSPLwrSdLQdElQm7TfLwZOrKpbhxiPJElAt+egvpjk\nSuB/gD9Jshi4Z7hhSZIWui4tqKOBfYBlVfVz4G6aN+xKkjQ0XRLUd6rqtqq6H6Cq7gK+MtywJEkL\n3Tq7+JI8DtgR2CzJnkDaTVsCjx5BbJKkBWyqMajfAl4P7AQMTgh7J/BXQ4xJkqR1J6iqOg44Lsnv\nVNXnRxiTJEmd7uL7UpJXAUsH96+qdw4rKEmSuiSo04DVwAXAvcMNR5KkRpcEtVNV7T/0SCRJGtDl\nNvNvJ3nq0CORJGlAlxbUrwOvT/LfNF18AaqqnjbUyCRJC1qXBPWi2Zw4yabAN4FHtdc5paqOTrIr\ncBKwCLgQeG1V3Teba0iSxte0XXxV9UNgZ2C/dvnuLsfRtLb2q6qnA3sA+yfZG3gf8KGq2h24DThs\ntsFLksZXl/dBHQ28DXh7W7QJcPx0x1XjZwPHbAIUsB9wSlt+HHDwDGOWJC0AXVpCv00zOexdAFV1\nI7BFl5Mn2SjJxcBK4EzgGuD2qlrT7nI9zXRKkx17eJIVSVasWrWqy+UkSWOkS4K6r6qKpvVDks27\nnryq7q+qPWimS9oLePJku63j2GOqallVLVu8eHHXS0qSxkSXBHVyko8DWyd5A/B14F9mcpGquh04\nB9i7Pc/EzRk7ATfO5FySpIWhy00S76cZM/o88CTgqKr66HTHJVmcZOt2eTPgN4ErgLOBl7e7HUIz\nU4UkSQ8z7W3mSf4c+FxVnTnDcy+hmWx2I5pEeHJVfSnJfwEnJXk3cBFw7EyDliSNvy7PQW0JfDXJ\nrTTPL51SVTdPd1BVXQLsOUn5tTTjUZIkrVOXLr53VNX/At4E7ACcm+TrQ49MkrSgdblJYsJK4CfA\nLcB2wwlHkqRGlwd1/zjJOcBZwLbAG5yHT5I0bF3GoB4PvLmqLh52MJIkTegyBnUk8Jgkh8KDt4/v\nOvTIJEkL2tDm4pMkaX0MdS4+SZJma6hz8UmSNFsjmYtPkqSZmvYuvqp6f5IXAHfw0Fx8M532SJKk\nGelymzltQjIpSZJGZiYzSUiSNDImKElSL60zQSU5q/1+3+jCkSSpMdUY1JIkzwEOTHISkMGNVXXh\nUCOTJC1oUyWoo4AjaV7L/sG1thWw37CCkiRpnQmqqk4BTknyf6vqXSOMSZKkTs9BvSvJgcCz26Jz\nqupLww1LkrTQdZks9j3AEcB/tZ8j2jJJkoamy4O6BwB7VNUDAEmOAy7iodnNJUna4Lo+B7X1wPJW\nwwhEkqRBXVpQ7wEuSnI2za3mz8bWkzS+ls/y/6DLV2/YOLTgdblJ4sQk5wDPpElQb6uqnww7MEnS\nwtZ1stibgNOHHIskSQ9yLj5JUi+ZoCRJvTRlgkryiCSXjSoYSZImTJmg2mefvp9klxHFI0kS0O0m\niSXA5UnOB+6aKKyqA4cWlSRpweuSoN4xmxMn2Rn4FPA44AHgmKr6SJJFwGeBpcB1wO9W1W2zuYY0\nL832OSNpgZn2JomqOpcmkWzSLn8P6PIuqDXAW6vqycDewJuSPIXmFR5nVdXuwFntuiRJD9Nlstg3\nAKcAH2+LdgROne64qrpp4qWGVXUncEV77EHAce1uxwEHzzxsSdK463Kb+ZuAfYE7AKrqB8B2M7lI\nkqXAnsB5wPbtg78TDwBPeq4khydZkWTFqlWrZnI5SdIY6JKg7q2q+yZWkmxM80bdTpI8Bvg88Oaq\nuqPrcVV1TFUtq6plixcv7nqYJGlMdElQ5yb5K2CzJC8APgd8scvJk2xCk5w+U1VfaItvTrKk3b4E\nWDnzsCVJ465LgjoSWAVcCrwROAP4m+kOShLgWOCKqvrgwKbTgUPa5UOA02YSsCRpYegym/kD7UsK\nz6Pp2ruqqrp08e0LvBa4NMnFbdlfAe8FTk5yGPAj4BWzilySNNamTVBJDgD+GbiG5nUbuyZ5Y1V9\nZarjqupb7f6Tef5MA5UkLSxdHtT9APC8qroaIMkTgC8DUyYoSZLWR5cxqJUTyal1Ld7YIEkasnW2\noJK8rF28PMkZwMk0Y1CvoJlNQpKkoZmqi++lA8s3A89pl1cBjx1aRJIkMUWCqqpDRxmIJEmDutzF\ntyvwv2lmH39wf1+3IUkapi538Z1K88DtF2lemyFJ0tB1SVD3VNU/DD0SSZIGdElQH0lyNPA14N6J\nwolXaUiSNAxdEtRTaaYs2o+HuviqXZckaSi6JKjfBnYbfOWGNHZm8xr25as3fBySHtRlJonvA1sP\nOxBJkgZ1aUFtD1yZ5Hs8fAzK28wlSUPTJUEdPfQoJElaS5f3QZ07ikAkSRrUZSaJO2nu2gN4JLAJ\ncFdVbTnMwCRJC1uXFtQWg+tJDgb2GlpEkiTR7S6+h6mqU/EZKEnSkHXp4nvZwOojgGU81OWnObb0\nyC/P6rjr3nvABo5kflt6zwkzP2i2P/tNZ3WYtOB0uYtv8L1Qa4DrgIOGEo0kSa0uY1C+F0qSNHJT\nvfL9qCmOq6p61xDi0YjMtmtwNuxOnDuz6bq8btNXDSESaeamakHdNUnZ5sBhwDaACUqSNDRTvfL9\nAxPLSbYAjgAOBU4CPrCu4yRJ2hCmHINKsgh4C/Bq4DjgGVV12ygCkyQtbFONQf098DLgGOCpVfWz\nkUUlSUMwm7HXUY6hjnJsGPo/PjzVg7pvBXYA/ga4Mckd7efOJHeMJjxJ0kI11RjUjGeZkCRpQzEJ\nSZJ6qctMEmOh733P0ry3fKtZHrd6w8ahsTG0FlSSTyRZmeSygbJFSc5M8oP2+7HDur4kaX4bZhff\nJ4H91yo7EjirqnYHzmrXJUn6BUNLUFX1TeDWtYoPonmeivb74GFdX5I0v416DGr7qroJoKpuSrLd\nunZMcjhwOMAuu+wyovA2DF+Boan0fX68Wb16ZD1cN9KrzZx/n+dOb+/iq6pjqmpZVS1bvHjxXIcj\nSRqxUSeom5MsAWi/V474+pKkeWLUCep04JB2+RDgtBFfX5I0TwxtDCrJicBzgW2TXA8cDbwXODnJ\nYcCPgFcM6/obwqjnxRr19frMfn9Nxb8rG0bfnw8dWoKqqleuY9Pzh3VNSdL46O1NEpKkhc0EJUnq\npQUzF580n4362SSpD2xBSZJ6yQQlSeolE5QkqZccg5I0p8b1maZxrdco2YKSJPWSCUqS1EsmKElS\nLzkGpaEbZV+8/f7S+LAFJUnqJROUJKmXTFCSpF4yQUmSeskEJUnqJROUJKmXTFCSpF4yQUmSeskE\nJUnqJROUJKmXTFCSpF4yQUmSeskEJUnqJROUJKmXTFCSpF4yQUmSeskEJUnqJROUJKmXTFCSpF6a\nkwSVZP8kVyW5OsmRcxGDJKnfRp6gkmwEfAx4EfAU4JVJnjLqOCRJ/TYXLai9gKur6tqqug84CTho\nDuKQJPXYxnNwzR2BHw+sXw88a+2dkhwOHN6u/izJVet53W2Bn67nOfps3OsH41/Hca8fWMd5L+/b\nIPV7fJed5iJBZZKy+oWCqmOAYzbYRZMVVbVsQ52vb8a9fjD+dRz3+oF1HAejrN9cdPFdD+w8sL4T\ncOMcxCFJ6rG5SFDfA3ZPsmuSRwK/D5w+B3FIknps5F18VbUmyZ8CXwU2Aj5RVZeP4NIbrLuwp8a9\nfjD+dRz3+oF1HAcjq1+qfmH4R5KkOedMEpKkXjJBSZJ6aSwSVJJPJFmZ5LKBskVJzkzyg/b7sW15\nkvxDO83SJUmeMXeRd5dk5yRnJ7kiyeVJjmjLx6KeSTZNcn6S77f1e0dbvmuS89r6fba9sYYkj2rX\nr263L53L+GciyUZJLkrypXZ9bOqY5Loklya5OMmKtmwsfkcnJNk6ySlJrmz/Pu4zTnVM8qT2z2/i\nc0eSN89FHcciQQGfBPZfq+xI4Kyq2h04q12HZoql3dvP4cA/jSjG9bUGeGtVPRnYG3hTmimixqWe\n9wL7VdXTgT2A/ZPsDbwP+FBbv9uAw9r9DwNuq6pfAj7U7jdfHAFcMbA+bnV8XlXtMfCszLj8jk74\nCPAfVfXLwNNp/izHpo5VdVX757cH8KvA3cC/Mxd1rKqx+ABLgcsG1q8ClrTLS4Cr2uWPA6+cbL/5\n9AFOA14wjvUEHg1cSDPDyE+BjdvyfYCvtstfBfZplzdu98tcx96hbjvR/OXeD/gSzYPrY1NH4Dpg\n27XKxuZ3FNgS+O+1/xzGqY5r1euFwH/OVR3HpQU1me2r6iaA9nu7tnyyqZZ2HHFs66Xt6tkTOI8x\nqmfb9XUxsBI4E7gGuL2q1rS7DNbhwfq121cD24w24ln5MPCXwAPt+jaMVx0L+FqSC9JMVwZj9DsK\n7AasAv6t7ab91ySbM151HPT7wInt8sjrOM4Jal06TbXUV0keA3weeHNV3THVrpOU9bqeVXV/Nd0K\nO9FMKvzkyXZrv+dd/ZK8BFhZVRcMFk+y67ytI7BvVT2DptvnTUmePcW+87F+GwPPAP6pqvYE7uKh\nrq7JzMc6AtCOhR4IfG66XScp2yB1HOcEdXOSJQDt98q2fN5OtZRkE5rk9Jmq+kJbPHb1rKrbgXNo\nxtq2TjLxQPlgHR6sX7t9K+DW0UY6Y/sCBya5jmYW//1oWlRjU8equrH9XkkzbrEX4/U7ej1wfVWd\n166fQpOwxqmOE14EXFhVN7frI6/jOCeo04FD2uVDaMZsJspf1955sjeweqLZ2mdJAhwLXFFVHxzY\nNBb1TLI4ydbt8mbAb9IMPp8NvLzdbe36TdT75cA3qu0A76uqentV7VRVS2m6Tr5RVa9mTOqYZPMk\nW0ws04xfXMaY/I4CVNVPgB8neVJb9HzgvxijOg54JQ9178Fc1HGuB+E20EDeicBNwM9psvlhNH31\nZwE/aL8XtfuG5oWJ1wCXAsvmOv6Odfx1mmbzJcDF7efF41JP4GnARW39LgOOast3A84HrqbpanhU\nW75pu351u323ua7DDOv7XOBL41THth7fbz+XA3/dlo/F7+hAPfcAVrS/q6cCjx3DOj4auAXYaqBs\n5HV0qiNJUi+NcxefJGkeM0FJknrJBCVJ6iUTlCSpl0xQkqReMkFJknrJBCVJ6iUTlDTHkjyzfY/O\npu1sDJcn+ZW5jkuaaz6oK/VAknfTzByxGc1cb++Z45CkOWeCknqgnTn6e8A9wK9V1f1zHJI05+zi\nk/phEfAYYAualpS04NmCknogyek0r+DYleZtpH86xyFJc27j6XeRNExJXgesqaoTkmwEfDvJflX1\njbmOTZpLtqAkSb3kGJQkqZdMUJKkXjJBSZJ6yQQlSeolE5QkqZdMUJKkXjJBSZJ66f8D6+VR80DF\nxpoAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "nDataGen = 100 # Number of signal (peak) events to simulate\n", "nDataBkg = 300 # Number of background events (flat distribution)\n", "\n", "xMin = 100. # minimum of data \n", "xMax = 700. # maximum of data\n", "\n", "genMean = 435.67 # peak poisition\n", "genSigma = 23.647 # peak width\n", "\n", "from scipy.stats import norm # normal distribution i.e. Gaussian\n", "from scipy.stats import uniform # uniform distribution i.e. flat\n", "# generate the data for signal and background as \"events\"\n", "sigArray = norm.rvs(genMean,genSigma,size=nDataGen)\n", "bkgArray = uniform.rvs(xMin,xMax,size=nDataBkg)\n", "\n", "bins = np.linspace(xMin,xMax,26) # 600 is the range, 25 is a nice factor of that + 1 as more bin edges than bins\n", "s_data,s_bins,patches = plt.hist([bkgArray,sigArray],bins,histtype=\"barstacked\")\n", "\n", "plt.xlabel('x')\n", "plt.ylabel('Number of events per bin')\n", "plt.title('Data and background')\n", "plt.tight_layout() # make space for labels and titles\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VOW5wPHfExIShIQ1YYcEyk4gQEBQqyIWN0S0eq1S\npdZWvba29rZu1WpVemuvXq1ebcUWt2oVtXUp1VZlESgWBNk3WQRkSwIIISyBJM/945yJA2SZTHLm\nnJk8389nPmfmzFmed5jhyfue97yvqCrGGGNM0CT5HYAxxhhTFUtQxhhjAskSlDHGmECyBGWMMSaQ\nLEEZY4wJJEtQxhhjAskSlDEeEZHnRWRyNe99R0TmeXDO2SLyvYY+rte8+jxMfLMEZQJDRDaLyGER\nOSAi+0RkvojcJCIRfU9FJFtEVESSvY7VGOM9S1AmaC5W1XSgO/AQcAcw1d+QGh9L8iYILEGZQFLV\n/ar6DnAlMElEBgKIyEUiskREikXkCxH5Zdhuc9zlPhEpEZFRItJTRGaKyB4R2S0iL4tIq+rOKyKP\nu8ctFpHFIvL1sPd+KSKviciLbi1vlYjkh70/REQ+dd+bBqTVUkwRkf8Tkf0islZExoS9cZ2IrHGP\ntUlEbjxhx0tEZKkb50YROb+Kg3cUkeUi8jP3dY6IzHGP+aGIPCUiL7nvhWqf14vIVmCmu368W859\nbvNhv7Djq4h8Lex1ZZOmiJwtIttE5KciUigiO0XkurBt24rIO278C4GetXxWphGyBGUCTVUXAtuA\nUKI4CFwLtAIuAv5TRCa4753pLlupagtV/RgQ4NdAJ6Af0BX4ZQ2n/ATIA9oAfwZeF5HwRDMeeNU9\n/zvAkwAi0hR4C/iTu+/rwDdrKd6pwCagHXAf8FcRaeO+VwiMAzKA64DHRGSoe64RwIvAbW4cZwKb\nww8sItnAR8CTqvqIu/rPwEKgrfsZXFNFTGfhfE7niUhv4BXgViATeBf4m1vWSHQAWgKdgeuBp0Sk\ntfveU8ARoCPwXfdhzHEsQZl4sAPnP31UdbaqrlDVClVdjvMf6FnV7aiqG1T1A1UtVdUi4NFatn9J\nVfeoapmq/i+QCvQJ22Seqr6rquU4yWiwu34kkAL8VlWPqeobOMmuJoVh208D1uEkXVT176q6UR0f\nAe/zVZK+HnjWLVeFqm5X1bVhx+0PzAbuU9VnAESkGzAcuFdVj6rqPJwEe6JfqupBVT2MU3v9u3ue\nY8AjQDPgtFrKFXIMeMAt37tACdBHRJrgJO973XOtBF6I8JimEbEEZeJBZ2AvgIicKiKzRKRIRPYD\nN+HUQKokIlki8qqIbBeRYuClWrb/qdu0tl9E9uHUAMK33xX2/BCQ5l6v6QRs1+NHX95SS7mq2r6T\nG8cFIvJvEdnrxnFhWBxdgY01HHcisB14I2xdJ2Cvqh4KW/dFFfuGr+sUXgZVrXDf71xjqb6yR1XL\nwl4fAlrg1MaSTzhXbZ+VaYQsQZlAE5HhOP8hhrog/xnnL/+uqtoSeBqnGQ+gqqH5f+2uH6SqGcC3\nw7Y/8Vxfx+mU8R9Aa1VtBeyvbvsT7AQ6i0j4tt1q2aeq7XeISCrwF5waS3s3jnfD4viCmq/Z/BLY\nDfzZra2E4msjIqeEbde1in3DP8MdOJ1VAOeCmbvPdnfVISD8eB1qiClcEVB2wvlr+6xMI2QJygSS\niGSIyDic6z0vqeoK9610nJrAEfdazNVhuxUBFUCPsHXpOE1L+0SkM851m+qk4/zHWQQki8i9ONeA\nIvGxu++PRCRZRC4DRtSyT5a7fYqIXIFz7eddoClO02IRUCYiFwBjw/abClwnImNEJElEOotI37D3\njwFXAM2BP4lIkqpuARYBvxSRpiIyCri4lvheAy5yz5MC/BQoBea77y8FrhaRJm4njWqbTsO5zaN/\ndWM5RUT6A5Mi2dc0LpagTND8TUQO4NQS7sa5ZnRd2Ps3Aw+429yL858oAG7z1a+Af7m9zkYC9wND\ncWpCf8f5j7E6/wTeAz7DaXI6QtXNYCdR1aPAZcB3gC9xrt/UdC6ABUAvnNrOr4DL3etfB4AfuWX7\nEicJV14vcjuOXAc85pbrI8JqOifEkwU8K869ZBOBUcAeYDIwDSfhVFemdTg1zv9zY7wY5zaAo+4m\nP3bX7XOP/VYt5Q33Q5zmvl3A88BzddjXNBJiExYa0zi5XeHXqup9fsdiTFWsBmVMIyEiw8W5LyzJ\nbZK7hLrVeoyJKbtb3JjGowNOs2NbnHvL/lNVl/gbkjHVsyY+Y4wxgWRNfMYYYwIpLpr42rVrp9nZ\n2X6HYYwxpgEsXrx4t6pm1rZdXCSo7OxsFi1a5HcYxhhjGoCIRDRyiDXxGWOMCSRLUMYYYwLJEpQx\nxphAiotrUMaY2Dp27Bjbtm3jyJEjfodi4lhaWhpdunQhJSUlqv0tQRljTrJt2zbS09PJzs7m+AHX\njYmMqrJnzx62bdtGTk5OVMewJj5jzEmOHDlC27ZtLTmZqIkIbdu2rVct3BKUMaZKlpxMfdX3O2RN\nfMYEwedzYN9WaNcbutY2jZQxjYPVoIzxU/kxePc2eOFiePsHMPUbMONBqKjwOzLfNWnShLy8vMrH\n5s2bG+zY+/bt43e/+13l6x07dnD55Zc3yLGzs7PZvXt3xNuXlpZy7rnnkpeXx7Rp0xokhpDZs2cz\nf/78ytff+c53eOONN6I+3uuvv06/fv0YPXo0ixYt4kc/+lGV52koVoMyxk8f/Q8sfAZG/RBGfB/m\nPAxzH4HmmTDyJr+j81WzZs1YunSpJ8cOJaibb74ZgE6dOtXrP+76WLJkCceOHatTWcvLy2nSpEmt\n282ePZsWLVpw2mmn1SfESlOnTuV3v/sdo0ePBiA/P9+T84RYDcoYvxSshnmPwaAr4bxfQetsGP8k\nfO1cmPEA7ItoMt9G5fnnn+eHP/xh5etx48Yxe/ZsAFq0aMHdd9/N4MGDGTlyJAUFBQAUFBRw6aWX\nMnjwYAYPHsz8+fO588472bhxI3l5edx2221s3ryZgQMHAk4Hkeuuu47c3FyGDBnCrFmzKs992WWX\ncf7559OrVy9uv/32auN8+OGHGTFiBCNGjGDDhg0AFBUV8c1vfpPhw4czfPhw/vWvf1FYWMi3v/1t\nli5dSl5eHhs3bmTGjBkMGTKE3Nxcvvvd71Ja6kx6nJ2dzQMPPMAZZ5zB66+/zsaNGzn//PMZNmwY\nX//611m7du1xMWzevJmnn36axx57jLy8PObOnQvAnDlzOO200+jRo8dxSfnhhx9m+PDhDBo0iPvu\nO3kOywceeIB58+Zx0003cdtttzF79mzGjRtX7XkagtWgjPHL+3dDajqc999frROBcY/BUyPhw1/C\n5VN9C6/Se3fCrhUNe8wOuXDBQzVucvjwYfLy8gDIycnhzTffrHH7gwcPMnLkSH71q19x++2384c/\n/IF77rmHH/3oR5x11lm8+eablJeXU1JSwkMPPcTKlSsray3hzYdPPfUUACtWrGDt2rWMHTuWzz77\nDIClS5eyZMkSUlNT6dOnD7fccgtdu3Y9KZaMjAwWLlzIiy++yK233sr06dP58Y9/zE9+8hPOOOMM\ntm7dynnnnceaNWv44x//yCOPPML06dM5cuQIZ599NjNmzKB3795ce+21/P73v+fWW28FnPuK5s2b\nB8CYMWN4+umn6dWrFwsWLODmm29m5syZlTFkZ2dz00030aJFC372s58BTg1o586dzJs3j7Vr1zJ+\n/Hguv/xy3n//fdavX8/ChQtRVcaPH8+cOXM488wzK4937733MnPmTB555BHy8/Mr/zCo6jwNxRKU\nMX4oXAsbZ8I5v4Dm7Y5/r1U3GPYdWDgFih+EjE6+hOi3ujbxNW3alHHjxgEwbNgwPvjgAwBmzpzJ\niy++CDjXtVq2bMmXX35Z7XHmzZvHLbfcAkDfvn3p3r17ZYIaM2YMLVu2BKB///5s2bKlygR11VVX\nVS5/8pOfAPDhhx+yevXqym2Ki4s5cODAcfutW7eOnJwcevfuDcCkSZN46qmnKhPUlVdeCUBJSQnz\n58/niiuuqNw3VNOqzYQJE0hKSqJ///6Vtcz333+f999/nyFDhlQef/369cclKD9YgjLGDwuehuQ0\nGHZd1e+P+D78+3fwyVQY84vYxnaiWmo6sZScnExFWAeS8HtsUlJSKrs1N2nShLKysqjOUdMkrqmp\nqZXPazpHePfq0POKigo+/vhjmjVrFtW5AZo3b155rFatWkV1jS68DKHzqSp33XUXN954Y52P5yW7\nBmVMrB0phmWvQu4V0Lxt1du0yYE+F8Li55yefgZwmpOWLl1KRUUFX3zxBQsXLqx1nzFjxvD73/8e\ncDoXFBcXk56eflLtJeTMM8/k5ZdfBuCzzz5j69at9OnTp05xhnrjTZs2jVGjRgEwduxYnnzyycpt\nqkouffv2ZfPmzZXXrf70pz9x1llnnbRdRkYGOTk5vP7664CTYJYtW3bSdjWVM9x5553Hs88+S0lJ\nCQDbt2+nsLCw1v3qep66sgRlTKx99g8oOwxDrql5uyHfhkN74POPYhNXHDj99NPJyckhNzeXn/3s\nZwwdOrTWfR5//HFmzZpFbm4uw4YNY9WqVbRt25bTTz+dgQMHcttttx23/c0330x5eTm5ublceeWV\nPP/888fVOiJRWlrKqaeeyuOPP85jjz0GwBNPPMGiRYsYNGgQ/fv35+mnnz5pv7S0NJ577jmuuOIK\ncnNzSUpK4qabqu7N+fLLLzN16lQGDx7MgAEDePvtt0/a5uKLL+bNN9+stfPC2LFjufrqqxk1ahS5\nublcfvnldUo4kZ6nrqS2KmUQ5Ofnq01YaBLGK1fDjiXwk1WQVMPfiMeOwMNfgwET4JInq9/OA2vW\nrKFfv34xPadJTFV9l0Rksarm17av1aCMiaUjxbDhQ+h/Sc3JCSAlDfpcAGunWzOfaZQsQRkTS+vf\nh/JSp1YUiQET4PCXsLnhmk2MiReeJSgR6SMiS8MexSJyq4i0EZEPRGS9u2ztVQzGBM76D+CUttAl\nwvH2ep4DTVJhwwxv4zImgDxLUKq6TlXzVDUPGAYcAt4E7gRmqGovYIb72pjEp+rc+9Tj7Nqb90JS\nmkH3Uc5+xjQysWriGwNsVNUtwCXAC+76F4AI2zqMiXMFq+BgoVMrqoue50Dhaije6U1cxgRUrBLU\nt4BX3OftVXUngLvMqmoHEblBRBaJyKKioqIYhWmMh0K1oB6j67ZfKKFtmtWw8TSwK6d8zJVTPvY7\nDJNAPE9QItIUGA+8Xpf9VPUZVc1X1fzMzExvgjMmljbOhMy+0LJz3fbLGgDNsxpdM19ouo3Bgwcz\ndOjQyukcGnJqjJpUN21GdnY2ubm55Obm0r9/f+65555ahxk6cXoPE5lY1KAuAD5V1QL3dYGIdARw\nl5HfrmxMvCo/Bl8shOyv133fpCTIPgM2/8u5jhVQR8sqWL2jmMID0U/xHS40Ft+yZcv49a9/zV13\n3QX4OzVGyKxZs1ixYgULFy5k06ZN3HDDDTVubwkqOrFIUFfxVfMewDvAJPf5JODk25+NSTS7VsCx\ng06Hh2h0Pw0O7HBm3Q2o7fsOc6C0jCc+XN/gxy4uLqZ1a6fDb/jUGDVNgVHd9BtVTXsBsGfPHsaO\nHcuQIUO48cYbax0XL3SOp59+mrfeeou9e/dSUlLCmDFjGDp0KLm5uZWjO5w4vUd125kTqKpnD+AU\nYA/QMmxdW5zee+vdZZvajjNs2DA1Jq7Nf1L1vgzV/duj23/ncmf/pa82bFzVWL16dcTb9r77Xe1+\nx/STHr3vfrdeMSQlJengwYO1T58+mpGRoYsWLVJV1c8//1wHDBigqqrPPfec5uTk6L59+/Tw4cPa\nrVs33bp1q6qqAvrOO++oquptt92mDz74oKqqXnXVVTp37lxVVd2yZYv27dtXVVVvueUWvf/++1VV\ndfr06QpoUVHRSXF17979pPWDBw/Wf//733rs2DHdv3+/qqoWFRVpz549taKi4riYVbXa7RJRVd8l\nYJFGkEM8Hc1cVQ+5CSl83R6cXn3GNB5bP4ZW3aOfOiOrP6S2hK3zYfCVDRtbPc29fTST313D9GU7\nqFBIS0nivAEduPui+g2VFD7dxscff8y1117LypUrT9quuikwqpt+o7ppL+bMmcNf//pXAC666KLK\nGlskNGxU8J///OfMmTOHpKQktm/fXllzO3H7qrbr0KFDxOdsDGy6DWO8pgpb/w096/F3WVIT6DrC\nOU7AZGWkkZ6aTIU68y2WllWQnppMVnpag51j1KhR7N69m6p69FY3BUZ102/UNO1F+DQZkTpw4ACb\nN2+md+/evPzyyxQVFbF48WJSUlLIzs4+bkqQkEi3a+xsqCNjvLZ3Exwsgm4j63ec7qOgaC0c2tsw\ncTWg3SWlZKWnMqBjBhNP7U5RSWST50Vq7dq1lJeX07ZtNdOT1EF1016ET7Px3nvv1TipYUhJSQk3\n33wzEyZMoHXr1uzfv5+srCxSUlKYNWsWW7ZsAU6ejqK67czxrAZljNd2LHGWXWodvLlmXYa7x/sU\nvnZu/Y7VwKZck195D9TkCQMb5JjhU76rKi+88AJNmjSp93GfeOIJfvCDHzBo0CDKyso488wzefrp\np7nvvvu46qqrGDp0KGeddRbdunWr9hijR49GVamoqODSSy/lF79wJpWcOHEiF198Mfn5+eTl5dG3\nb1+A46b3uOCCC7jjjjuq3M4cz6bbMMZr/7wbPvkj3LUNmqREf5wjxfBQNxj9czjr9tq3r4doptsI\nJahpN0bZU9EkpPpMt2E1KGO8tmMptB9Yv+QEkJYBmX1g++KGiauBWWIyDc2uQRnjpYoK2LkUOg1p\nmON1HuYkqDho+TCmvixBGeOlPRvgaEkDJqihToeL/V80zPFqEA/N/ybY6vsdsgRljJd2Oj3E6JTX\nMMfrPMxZbvP2mmxaWhp79uyxJGWipqrs2bOHtLTobzewa1DGeGnHEkhuBu36NMzxsgZAUgrsWg4D\nL2uYY1ahS5cubNu2rcr7joyJVFpaGl26dIl6f0tQxnhpxxLoOAiaNNBPLbkpZPV1xvbzUEpKCjk5\nOZ6ew5jaWBOfMV6pKIedy6BjAzXvhXQY5HmCMiYILEEZ45Xd6+HYoYbrIBHSIRdKCuDAyWO8GZNI\nLEEZ45XQCBJeJCiwWpRJeJagjPHKjiWQ0hza9WrY47Z3hxLatbxhj2tMwFiCMsYroQ4SSfUfP+44\nzVpBq25WgzIJzxKUMV6oqICCVU6HBi9YRwnTCFiCMsYL+7c6U7y37+/N8TvkuqNUHPTm+MYEgCUo\nY7xQuMZZZnmVoAYBCgWra93UmHhlCcoYLxS6iSOzgUaQOFFlT77oOkpcOeXjyukxjAmqWhOUiJwu\nIh+IyGcisklEPheRTZEcXERaicgbIrJWRNaIyCgRaeMeb727bF3/YhgTMIVrIaMLpLX05vgtu0Ba\nK7sOZRJaJDWoqcCjwBnAcCDfXUbiceAfqtoXGAysAe4EZqhqL2CG+9qYxFK4BrLqNuFfnYg4tShL\nUCaBRZKg9qvqe6paqKp7Qo/adhKRDOBMnASHqh5V1X3AJcAL7mYvABOijN2YYCovg93rvE1Q4FyH\nKljlDKlkTAKKJEHNEpGH3ea5oaFHBPv1AIqA50RkiYj8UUSaA+1VdSeAu8yqamcRuUFEFonIIhtR\n2cSVvZug/Kh3HSRCOgyEssOwZ6O35zHGJ5EMsXyquwyfP16BcyI49lDgFlVdICKPU4fmPFV9BngG\nID8/3yalMfEj1EHC6xpUKAEWrobM3t6eyxgf1JqgVHV0lMfeBmxT1QXu6zdwElSBiHRU1Z0i0hEo\njPL4xgRT4RpAvOvBF5LZxzlP4RoYYC3lJvFUm6BE5Nuq+pKI/FdV76vqozUdWFV3icgXItJHVdcB\nY4DV7mMS8JC7fDvq6I0JosLV0KYHpDTz9jwpzaBNDhSt8fY8xvikphpUc3eZXo/j3wK8LCJNgU3A\ndTjXvV4TkeuBrcAV9Ti+McHjdQ++cFn9v7op2JgEU22CUtUp7vL+aA+uqks5/tpVyJhoj2lMoB07\nAns3xq7JLbMvrHsPykohOTU25zQmRiK5UbeHiPxNRIpEpFBE3haRHrEIzpi4s/sz0IoY1qD6gZY7\nkyMak2Ai6Wb+Z+A1oCPQCXgdeMXLoIyJW0VrnaXXXcxDQokwdF5jEkgkCUpU9U+qWuY+XsLpZm6M\nOVHhakhKgTY9Y3O+tr0gKfmrru3GJJCaevG1cZ/OEpE7gVdxEtOVwN9jEJsx8adwjTODbnLT2Jwv\nuamTDK2jhElANfXiW4yTkMR9fWPYewo86FVQxsStwtXQJdKhKhtIVj/YuSy25zQmBmrqxZcTy0CM\niXulB2DfVhh6bWzPm9UfVr8NRw9B01Nie25jPGTzQRnTUIrWOctYdZAIyeoLqDNArTEJxBKUMQ0l\nVmPwnahyTD67DmUSS40JShxdYxWMMXGtcA0kN4NW2bE9b+scaNK0TgnqaFkFq3cUU3jgiIeBGVM/\nNSYoVVXgrRjFYkx8K1ztDOCaFOOGiSbJ0K5PnRLU9n2HOVBaxhMf2g2+JrgimW7j3yIyXFU/8Twa\nY+JZ4VroWdssNB7J6gtbPq51sz73vEdpWUXl65cWbOWlBVtJTU5i3eQLvIzQmDqL5E+90ThJaqOI\nLBeRFSKy3OvAjIkrh/ZCya7YX38KyeoHxdvgSHGNm829fTTj8zqR5N48kpaSxCV5nZh7R7Sz6hjj\nnUhqUPZnlTG1CTWvxboHX0hm2JBHXUdUu1lWRhrpqclUKIhAaVkF6anJZKWnxShQYyJXaw1KVbcA\nXYFz3OeHItnPmEbFrx58IaHzRnAdandJKVnpqQzomMHEU7tTVFLqcXDGRKfWGpSI3IczZUYf4Dkg\nBXgJON3b0IyJI4VrILUlZHTy5/ytukPKKRElqCnX5HPlFOd61eQJA72OzJioRVITuhQYDxwEUNUd\n1G8SQ2MST2iSQpHat/VCUpLTg9Bm1zUJJJIEddTtbq4AItK8lu2NaVxUnSa+rL7+xmGz65oEE0mC\nek1EpgCtROT7wIfAH7wNy5g4cmAXHNnnXweJkMy+UFLg9Cg0JgHUeg1KVR8RkW8AxUBv4F5V/cDz\nyIyJF353kAipHPJoNWSf4W8sxjSASLqZA6wAmuE0862I9OAishk4AJQDZaqa784zNQ3IBjYD/6Gq\nX0YesjEBE+tZdKsTamIsXGMJyiSEWpv4ROR7wELgMuBynJt2v1uHc4xW1TxVzXdf3wnMUNVewAz3\ntTHxq3A1NM+E5u38jSOjM6Rm2PTvJmFEUoO6DRiiqnsARKQtMB94NspzXgKc7T5/AZgN3BHlsYzx\nX6gHn99EnOtQ1lHCJIhIOklsw2mmCzkAfBHh8RV4X0QWi8gN7rr2qroTwF1mVbWjiNwgIotEZFFR\nUVGEpzMmxioqnDH4/G7eC8nq5yQoVb8jMabeIqlBbQcWiMjbOAnnEmChiPwXgKo+WsO+p6vqDhHJ\nAj4QkYjbHlT1GeAZgPz8fPu1mWDavxWOHQxGDQqcRPnpC1BSCOnt/Y7GmHqJJEFtdB8hb7vLWm/W\ndW/qRVULReRNYARQICIdVXWniHQECusYszHB4fcYfCcKdZQoWmMJysS9SLqZ3x/Ngd0bepNU9YD7\nfCzwAPAOMAl4yF2+Xf1RjAm4UBfzzD7+xhESPrtuj7P9jMSYeou0m3k02gNvijP0SzLwZ1X9h4h8\ngnPz7/XAVuAKD2MwxluFayCjC6S19DsSR/NMaNbGOkqYhOBZglLVTcDgKtbvAcZ4dV5jYiooPfhC\nRGzII5MwbNoMY6JVXga7PwtWggLnOlTRWuvJZ+JeJDfq/o+IZIhIiojMEJHdIvLtWARnTKDt3QTl\nR4PTQSIkqx+UFkPxdr8jMaZeIqlBjVXVYmAczj1RvXFu3jWmcQvKGHwnCs2uW2gjSpj4Fsk1qBR3\neSHwiqruFb/mvDEmSArXABKcHnwhlbPrroZe51a5ybQbR8UwIGOiE0mC+pt7g+1h4GYRyQSOeBuW\nMXGgcDW06QEpzXwLITQz7nEJ55Q20KKDdZQwcS+SJr77gFFAvqoeAw7hzLBrTOMWtB584bL62uy6\nJu5FkqA+VtUvVbUcQFUPAu95G5YxAXfsCOzdGOAE1R+K1jljBRoTp6pt4hORDkBnoJmIDAFCF54y\ngFNiEJsxwbV7HWhF8HrwhWT2hWOHYN8WaJPjdzTGRKWma1DnAd8BugDhA8IeAH7uYUzGBF+B24Ov\n/QB/46hO+JBHlqBMnKo2QanqC8ALIvJNVf1LDGMyJvgKVkKTVGjT0+9IqhbqWVi0Bvpe6G8sxkQp\nkl5800Xkapwp2iu3V9UHvArKmMArXO0kgSZeDmdZD2kZ0LKr9eQzcS2SX9fbwH5gMVDqbTjGxImC\n1dBztN9R1Cyzr92sa+JaJAmqi6qe73kkxsSLg3ugZFdwO0iEZPWDzz9yxgwMak3PmBpE0s18vojk\neh6JMfGicJWzbB8HCar8KHz5ud+RGBOVSBLUGcBiEVknIstFZIWILPc6MGMCq7IH30B/46hN6B6t\nglX+xmFMlCKp91/geRTGxJPCVc6kgC0CPqV6Zl+QJCdBDZjgdzTG1FmtNShV3QJ0Bc5xnx+KZD9j\nElbBauf+p6APmpzSDNr2crrEGxOHIpkP6j7gDuAud1UK8JKXQRkTWBUV7hh8Ab/+FNJhIOyyBGXi\nUyQ1oUtxBoc9CKCqO4B0L4MyJrD2bYZjB4M7gsSJ2g+E/Vvh8Jd+R2JMnUWSoI6qqgIKICLN63IC\nEWkiIktEZLr7OkdEFojIehGZJiJN6x62MT4J2BBHR8sqWL2jmMID1cyA08HtgGsdJUwciiRBvSYi\nU4BWIvJ94EPgD3U4x4+B8NvZfwM8pqq9gC+B6+twLGP8FZpFN7Ovv3G4tu87zIHSMp74cH3VG4R6\nGlozn4lDtfbiU9VHROQbQDHQB7hXVT+I5OAi0gW4CPgV8F/iTMV7DnC1u8kLwC+B39c9dGN8ULAK\nWmdDagtykZefAAAapElEQVRfw+hzz3uUln01lcZLC7by0oKtpCYnsW5yWMfb9A5wSjsoWOFDlMbU\nTySdJH4CrFHV21T1Z5EmJ9dvgduB0C+pLbBPVcvc19twpvSo6rw3iMgiEVlUVFRUh1Ma46GClYG4\n/2nu7aMZn9eJJLcjYVpKEpfkdWLuHScMvyRiHSVM3IqkiS8D+KeIzBWRH4hIRDd/iMg4oFBVF4ev\nrmJTrWp/VX1GVfNVNT8zMzOSUxrjrSPFsGcDdBzsdyRkZaSRnppMhTo5qLSsgvTUZLLS007euP1A\np+dhednJ7xkTYJE08d0P3C8ig4ArgY9EZJuqnlvLrqcD40XkQiANJ9H9FudaVrJbi+oC7KhXCYyJ\nldD9RAFIUAC7S0rJSk8lKz2VvG6tKaqpo0R5qZNcs4Jx7cyYSNTlhttCYBewB8iqbWNVvUtVu6hq\nNvAtYKaqTgRmAZe7m03CGS3dmODb6Y7w1WGQv3G4plyTT0675jRPTWbyhIFMuSa/6g1DTZJ2w66J\nM5Fcg/pPEZkNzADaAd9X1fr8Qu/A6TCxAeea1NR6HMuY2Nm5DJpnOR0P4km73pCUArtsCE0TXyIZ\ni687cKuqLo32JKo6G5jtPt8EjIj2WMb4Ztdy6Dgo+EMcnSi5qdMt3jpKmDgTyVh8dwItROQ6ABHJ\nFJEczyMzJkiOHXE6GgTk+lOddRhoTXwm7thYfMZEonA1aHkcJ6hcKCmAAwV+R2JMxGwsPmMisXOZ\nswxIB4k66zTEWe5Y4m8cxtSB52PxGZMQdi2H1JbOKBLxqMMgZ24oS1AmjsRiLD5j4t/OZfHZQSIk\ntQW06wM7PvU7EmMi5ulYfMYkhPIyZwy+4d/zO5L66TwUPvsnqMZvojWNSiTdzHETkiUl0zjt/gzK\njsRvB4mQTkNg6cuwfxu06up3NMbUyqZuN6Y28d5BIqTTUGdpzXwmTliCMqY2O5ZAyinQrpffkdRP\nh4HOiBLbLUGZ+FBtghKRGe7yN7ELx5gA2vYJdB4GSU38jqR+klOhfX/ryWfiRk01qI4ichbOiORD\nRGRo+CNWARrjq2OHnS7mXaoZiDXedBoKO5ZCRUXt2xrjs5o6SdwL3IkzJcajJ7ynODPjGpPYdi6D\nijLokiDDR3YeCoufg72boN3X/I7GmBpVm6BU9Q3gDRH5hao+GMOYjAmOLxY6yy7D/Y2jGtNuHFW3\nHcJHlLAEZQIuksFiHxSR8SLyiPsYF4vAjAmEbZ84o0e0SJBZnTP7QXIz68ln4kIkg8X+GvgxsNp9\n/NhdZ0xiU3USVEBrT1FpkuyMiLF9sd+RGFOrSG7UvQjIU9UKABF5AVjCV6ObG5OYirfDgZ2Jc/0p\npOsIWDDFmUIkJc3vaIypVqT3QbUKe97Si0CMCZzK608J0oMvpNtpUH7UalEm8CKpQf0aWCIiswAB\nzsRqT6Yx2LYIktOg/UC/I2lY3UY6y63zIft0f2MxpgaRDBb7iojMBobjJKg7VHWX14EZ47ttC51e\nb8lN/Y6kYZ3SxuksseVjvyMxpkYRNfGp6k5VfUdV3440OYlImogsFJFlIrJKRO531+eIyAIRWS8i\n00QkwX79JiGUlTr3QCVa815I91FOE2ZFud+RGFMtL8fiKwXOUdXBQB5wvoiMBH4DPKaqvYAvges9\njMGY6Oxc5lynSbQOEiHdToOjB2DXCr8jMaZaniUodZS4L1PcR2gEijfc9S8AE7yKwZiobZ7nLLvV\n8UbYeNHdLddWa+YzwVVjghKRJBFZGe3BRaSJiCwFCnHmk9oI7FPVMneTbUDnava9QUQWiciioqKi\naEMwJjpb/gWZfRPnBt0TtewCLbvBlvl+R2JMtWpMUO69T8tEpFs0B1fVclXNwxnPbwTQr6rNqtn3\nGVXNV9X8zMwE/U/CBFN5GWz9N2Sf4Xck3uo+yqlBaZU/QWN8F0kTX0dglYjMEJF3Qo+6nERV9wGz\ngZFAKxEJ9R7sAuyoy7GM8dzOZXC0JPETVLdRcLAI9mz0OxJjqhTJfVD3R3NgEckEjqnqPhFpBpyL\n00FiFnA58CowCXg7muMb45nNc51l9wS/R6j7ac5y63wbONYEUiSDxX4EbAZS3OefAJGMNNkRmCUi\ny919PlDV6cAdwH+JyAagLTA1ytiN8cbnH7nXn7L8jsRb7XpD8yz4fI7fkRhTpVprUCLyfeAGoA3Q\nE6dTw9PAmJr2U9XlwJAq1m/CuR5lTPAcO+x0HMj/rt+ReE8Eeo6GDTOcCQyTvLzrxJi6i+Qb+QPg\ndKAYQFXXAwn+p6VptLZ+DGVHoGcjmY+z5zlwaDcU1H4/1JVTPubKKdYt3cROJAmqVFWPhl64HRys\n249JTBtnQpOmX12fSXQ9znaWG2f6GYUxVYokQX0kIj8HmonIN4DXgb95G5YxPtk4yxlMtWlzvyOJ\njfQOzmC4G2b4HYkxJ4kkQd0JFAErgBuBd4F7vAzKGF/s3w4FK6FnjZdXE0/Pc5z7vo4U+x2JMceJ\npBdfBc6QRA/idDl/QdXu7DMJ6LP3nGWfC/yNI9b6XAAVx2Cj1aJMsEQy5ftFOEMUPQE8CWwQkUb2\nCzaNwrp/QOscp/t1Y9JlBDRr7ZTfmACJ5Ebd/wVGq+oGABHpCfwdeM/LwIyJqaMHnfuBhl/vdL+O\nsVDvuGk3+jA4bZNk6DUW1r/vTL+R1CT2MRhThUiuQRWGkpNrE87gr8Ykjg0zoLwUep/vdyT+6HMB\nHN7rXIsyJiCqrUGJyGXu01Ui8i7wGk738itwRoYwJnGsfgtOaZv4wxtV52vfcKa3X/2WTQNvAqOm\nGtTF7iMNKADOAs7G6dHX2vPIjImVo4ec6y/9xjvNXY1RagunmW/12zbLrgmMan+NqnpdLAMxxjfr\n34djB2HApX5H4q8Bl8Kad5yhnnK+ftLbR8sq2FBYQuGBI2Slp/kQoGlsIunFlyMij4rIX6OdbsOY\nQFv5F2ie2Xib90J6nwcppzifRxW27zvMgdIynvhwfYwDM41VJO0Zb+GMOP43oMLbcIyJsUN74bN/\nQP71jbd5L6Rpc+g7Dlb+Fc7/NaQ0A6DPPe9RWvbVT/+lBVt5acFWUpOTWDfZ7jgx3omkF98RVX1C\nVWep6kehh+eRGRMLK16H8qMwZKLfkQTDkIlQuh/W/r1y1dzbRzM+rxNJbu/7tJQkLsnrxNw7RvsU\npGksIklQj4vIfSIySkSGhh6eR2ZMLCx5CToMgg65fkcSDNlnQstuzufiyspIIz01mQp1bhErLasg\nPTXZrkMZz0XSppELXAOcw1dNfOq+NiZ+bV8Mu5bDBQ/7HUlwJCU5tajZDzlTwbftCcDuklKy0lPJ\nSk8lr1trig4c8TlQ0xhEUoO6FOihqmep6mj3YcnJxL8Fz0DTdBj8Lb8j4WhZBat3FFMYhP/4h33H\nGU3ikz9WrppyTT457ZrTPDWZyRMGMuWafP/iM41GJAlqGdDK60CMiakDBU5vtbyrIS3D72iC1UMu\nvQP0n+A085WW+B2NacQiaeJrD6wVkU+A0tBKVR3vWVTGeG3hM1BRBiNu8DWMwPaQO/UmWPkGfPoi\njLrZvzhMoxZJgrovmgOLSFfgRaADzrWrZ1T1cRFpA0wDsoHNwH+o6pfRnMOYqBz+0klQ/cdDu6/5\nGsrc20cz+d01TF+2gwp1esidN6ADd1/Uz9e46Docsr8O/3oc8r8LKdYhwsReJPNBfVTVI4JjlwE/\nVdV+wEjgByLSH2cCxBmq2guY4b42JnYWTIHSYjjzNr8jCXYPubNuh5JdTi3KGB9EMpLEAREpdh9H\nRKRcRGqdelNVd6rqp+7zA8AaoDNwCc4EiLjLCdGHb0wdlRTC/CedG1ID0rU81ENuQMcMJp7anaKS\n0tp3ioXsr0O302DOw1B6wO9oTCMUSQ0qXVUz3Eca8E2ciQsjJiLZwBBgAdBeVXe6x94JZFWzzw0i\nskhEFhUVFdXldMZUb9avoOwwnHt/rZteOeXjynmavBTYHnIiMHYyHCyEeY/5HY1phCLpxXccVX2L\nOtwDJSItgL8At6pqrTWvsPM8o6r5qpqfmZlZ1zCNOdmOJU5z1fDv+37tKW50GQa5/wHzn6RD2Xa/\nozGNTCRNfJeFPS4XkYdwbtStlYik4CSnl1X1r+7qAhHp6L7fEZv80MRC2VF4+4fQoj2cbZc96+Qb\nD0ByGjfu/y2iNhyniZ1IalAXhz3OAw7gXEeqkYgIziCza1T10bC33gEmuc8nAW/XJWBjovLRQ1Cw\nEi56FJrZbX11ktERzv9v+h9dwQUH3/I7GtOI1NrNvB7zQp2OM0TSChFZ6q77OfAQ8JqIXA9sxZmh\n1xjvrP8A5v4vDLkG+l7odzTxKW8iC//xMhMPTIUvroCuI/yOyDQCNU35fm8N+6mqPljTgVV1HiDV\nvD0mgtiMqb+CVfDG9dA+Fy60MfeiJsKIW1+BZ86CVyfC9z6A1tl+R2USXE1NfAereABcD9zhcVym\nnmLVAy3Q9myEl74JTU+Bq16pnN/Iawn72TdrBVe/5kxP8qfLoHiH3xGZBFdtglLV/w09gGeAZsB1\nwKtAjxjFZ0x0di6HZ893/jP99l+gVVe/I0oMmX1g4uvO/WTPnge7N/gdkUlgNXaSEJE2IjIZWI7T\nHDhUVe9QVet5Z4JJFZb+GaaOhSYpcN0/oP2AqA4VqBHGg6TrCJj0Dhw9CH8YDWum+x2RSVDVJigR\neRj4BKfXXq6q/tLGzEts0TZNBaZJa+/n8MpV8NZ/QudhcMNsyOwd9eECNcJ40HQe6ny+bXrAtInw\nl+9Zk59pcDX14vspzujl9wB3O73GAafjg6qq/3MUGKPqTDy46DlY/iokpTijH4y82ZnTKAp+jDA+\n7cZRnhzXU626wfUfwLxHYc4jsOZvkDfRmU+q4yC/ozMJoNoEpap1HmXCmJg4dgS2fQKfz4G1f4fC\nVZDSHIZOcgaAzehYr8MHdoTxIEpu6tz4PPhb8NH/OHNILZoKnYZAnwsh50zoNNTZzpg6imS6DWO8\npQplR+DYYecRen54r3Mx/mCRs/zycyhaB7vXQ8UxkCTonA/jHoOBlzfYxIOBHmE8qFpnw4TfObXX\n5a/Bsldg1n87Yx82SYV2vZ0OFq2zoUUWNG8HzbOgWWund2VKM0hO+2op1d2hYhqTxpGgfpvr/IcX\nolWN1FTFuki2i/hYDXS+Krc7eZsXj5U7TyYn1bhd+LH+VF7h3Lj2oFS7TVXHernCfX7S+Kv1KOOJ\npInTpJTZF3qNhW4joftpkNay9n2jEBphPCs9lbxurSmyjhKROaUNjLzJeRzaC5vnwRcLYPdnsG0h\nrHoTtLz240hSDQ/56jlS72S27/AxAFo1S6nXcbwUbYwNXzZxWihunNNAx6tZ40hQfS50uhsfp4ov\ndZVf9Ei2i/JYUZ+vqt2O3+a1JXt5/sAIXhm2lazUYxGd892lzkXuCUM61ymudz51BhG9bGjnk7eL\n9HMI/fWc0gySmzkT5DVr7fyV3SIL0lpBUuxanadck1/Z8WPyhIF12vdoWQUbCksoPHCkcde6Tmnj\nTArZP2zy7YqKsJpxIRze91WNOXypFVU8tOr19bRwdQEAY/u1r/exvBJtjA1attAflzEcKqxxJKgL\nfuN3BDH3f598SGFZKU+UDWPyxZHNe/TKJuc/5Ann1u2C/bQNzn6XjYnDC/0eCO/9N/nSYMw5FRhJ\nSW7zXjugv9/RADB1u/P9HTsuuN/faGOMh7LVpHEkqBqE/kqOy15UVahPD7Ro//K3GoPDj95/8SDR\nfmMmdqynXoKZe/toxud1IsltRUtLSeKSvE7MvWN0rftGe9+P3S/kqM9nb4w5WaOvQSWaaHqgRfuX\nv9UYjme9/4xpWFaDiqFYjbgQ6oE2oGMGE0/tTlFJaY3bR/uXv9UYTlbXz974L9ohrWI5gkq0Mcb7\ncF2WoBLQlGvyyWnXnOapyUyeMJAp1+TXuH20f/lbjeFkdf3sjf/ioYm6sTa/WxOfAaK/78fuFzLx\nKh6aqBt783ujr0FFUwWOZdU+VueK9i//WNcYYvnZT7txlPU8awCxbGaqy/cjHpqo/Wh+D8zgz1iC\nivsqsDFBF9TfWDw0UTf25vdG28SXKFVgY4IqHn5j8dBE3Zib3z2rQYnIsyJSKCIrw9a1EZEPRGS9\nu2zt1flr40f1Pt571BhTF/HQhFafJupY/Z7jpfndC1428T0PnH/CujuBGaraC5jhvvaFH1XgoDZ1\nGOOFRGlmqo79nr3nWROfqs4RkewTVl8CnO0+fwGYDdzhVQy1iVUV2CbAa1zss/9KIjQznSgemi4T\nRaw7SbRX1Z0A7jKrug1F5AYRWSQii4qKijwJJtoqcF2r9vVp6rBmwePZ5xFfYt3MFIvvRzw0XdZH\nkH5jge3Fp6rPqGq+quZnZmb6Hc5x6lq1r09ThzUjHM8+j8Yh2q7Osfh+1LfpMkjduKsSpN9YrHvx\nFYhIR1XdKSIdgcIYn79e6lO1r2tTRzw1C8aiScuaVUxNYv39sKbL2Ih1DeodYJL7fBLwdozPXy/1\nqdrXtakj0ZsR6so+D1OTWH8/EqGH3ImC+BvzrAYlIq/gdIhoJyLbgPuAh4DXROR6YCtwhVfn90Is\neyUleg8oqNs8QQ3RrBLpuUz8qc/3Ix6+E7Fo3Qji/zle9uK7qpq3xnh1zliIZdU+EZsR6sM+j8Yj\nmkkw4+X7EeQJPoP2GTbakSRC6vqXyZRr8iv/Gp88YaAXIflyrnhgn0d8iuav//AL9ZMvzY1on3j5\nfkRTtlgJ2mfY6BOUMSY4gnihvqEkctm8Ethu5ibxxfJ+iyDd22GqF8QL9Q0lkcvmFUtQMWTTNxwv\nlvdbBOneDlO9IF6or05df8/xVLagsCY+E3OxbOqwZpX4E7QL9Q0pkcvmBVFVv2OoVX5+vi5atMjv\nMEwDKSw+wuR31zB92Q4q1GnqOG9AB+6+qF+D/zUZy3OZhpPItwUkctkiJSKLVbXWm8esic/EnN1P\nZoyJhDXxGV/Y/WTGmNpYgjK+sPvJjDG1sSY+Y4wxgWQ1KGNM4CRyB4JELltDsxqUMcaYQLJu5sYY\nY2LKupkbY4yJa5agjDHGBJIlKGOMMYFkCcoYY0wgWYIyxhgTSJagjDHGBJIvCUpEzheRdSKyQUTu\n9CMGY4wxwRbzBCUiTYCngAuA/sBVItI/1nEYY4wJNj9qUCOADaq6SVWPAq8Cl/gQhzHGmADzYyy+\nzsAXYa+3AaeeuJGI3ADc4L4sEZF19TxvO2B3PY8RZIlePkj8MiZ6+cDKmAgaonzdI9nIjwQlVaw7\nabwlVX0GeKbBTiqyKJKhNeJVopcPEr+MiV4+sDImgliWz48mvm1A17DXXYAdPsRhjDEmwPxIUJ8A\nvUQkR0SaAt8C3vEhDmOMMQEW8yY+VS0TkR8C/wSaAM+q6qoYnLrBmgsDKtHLB4lfxkQvH1gZE0HM\nyhcX020YY4xpfGwkCWOMMYFkCcoYY0wgJUSCEpFnRaRQRFaGrWsjIh+IyHp32dpdLyLyhDvM0nIR\nGepf5JETka4iMktE1ojIKhH5sbs+IcopImkislBElrnlu99dnyMiC9zyTXM71iAiqe7rDe772X7G\nXxci0kRElojIdPd1wpRRRDaLyAoRWSoii9x1CfEdDRGRViLyhoisdX+PoxKpjCLSx/33Cz2KReRW\nP8qYEAkKeB44/4R1dwIzVLUXMMN9Dc4QS73cxw3A72MUY32VAT9V1X7ASOAH4gwRlSjlLAXOUdXB\nQB5wvoiMBH4DPOaW70vgenf764EvVfVrwGPudvHix8CasNeJVsbRqpoXdq9MonxHQx4H/qGqfYHB\nOP+WCVNGVV3n/vvlAcOAQ8Cb+FFGVU2IB5ANrAx7vQ7o6D7vCKxzn08Brqpqu3h6AG8D30jEcgKn\nAJ/ijDCyG0h2148C/uk+/ycwyn2e7G4nfsceQdm64Py4zwGm49y4njBlBDYD7U5YlzDfUSAD+PzE\nf4dEKuMJ5RoL/MuvMiZKDaoq7VV1J4C7zHLXVzXUUucYx1YvblPPEGABCVROt+lrKVAIfABsBPap\napm7SXgZKsvnvr8faBvbiKPyW+B2oMJ93ZbEKqMC74vIYnGGK4ME+o4CPYAi4Dm3mfaPItKcxCpj\nuG8Br7jPY17GRE5Q1YloqKWgEpEWwF+AW1W1uKZNq1gX6HKqark6zQpdcAYV7lfVZu4y7sonIuOA\nQlVdHL66ik3jtozA6ao6FKfZ5wcicmYN28Zj+ZKBocDvVXUIcJCvmrqqEo9lBMC9FjoeeL22TatY\n1yBlTOQEVSAiHQHcZaG7Pm6HWhKRFJzk9LKq/tVdnXDlVNV9wGyca22tRCR0Q3l4GSrL577fEtgb\n20jr7HRgvIhsxhnF/xycGlXClFFVd7jLQpzrFiNIrO/oNmCbqi5wX7+Bk7ASqYwhFwCfqmqB+zrm\nZUzkBPUOMMl9Pgnnmk1o/bVuz5ORwP5QtTXIRESAqcAaVX007K2EKKeIZIpIK/d5M+BcnIvPs4DL\n3c1OLF+o3JcDM9VtAA8qVb1LVbuoajZO08lMVZ1IgpRRRJqLSHroOc71i5UkyHcUQFV3AV+ISB93\n1RhgNQlUxjBX8VXzHvhRRr8vwjXQhbxXgJ3AMZxsfj1OW/0MYL27bONuKzgTJm4EVgD5fscfYRnP\nwKk2LweWuo8LE6WcwCBgiVu+lcC97voewEJgA05TQ6q7Ps19vcF9v4ffZahjec8GpidSGd1yLHMf\nq4C73fUJ8R0NK2cesMj9rr4FtE7AMp4C7AFahq2LeRltqCNjjDGBlMhNfMYYY+KYJShjjDGBZAnK\nGGNMIFmCMsYYE0iWoIwxxgSSJShjjDGBZAnKGGNMIFmCMsZnIjLcnUcnzR2NYZWIDPQ7LmP8Zjfq\nGhMAIjIZZ+SIZjhjvf3a55CM8Z0lKGMCwB05+hPgCHCaqpb7HJIxvrMmPmOCoQ3QAkjHqUkZ0+hZ\nDcqYABCRd3Cm4MjBmY30hz6HZIzvkmvfxBjjJRG5FihT1T+LSBNgvoico6oz/Y7NGD9ZDcoYY0wg\n2TUoY4wxgWQJyhhjTCBZgjLGGBNIlqCMMcYEkiUoY4wxgWQJyhhjTCBZgjLGGBNI/w+gPMymgD8g\ntQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "def fitFunc(x,A,mu,sigma,B):\n", " \"\"\"Function to fit is A*exp(-(x-mu)**2/2*sigma**2) + B\"\"\"\n", " return A*np.exp(-((x-mu)**2/(2*sigma**2))) + B\n", "\n", "# set initial values for the fit, based on obersvation of plot above\n", "A = 60\n", "B = 15\n", "mu = 400\n", "sigma = 30\n", "\n", "# correct from bin edges to the middle of the bins\n", "cbins = np.array([(s_bins[i]+s_bins[i+1])/2. for i in range(len(s_bins)-1)])\n", "\n", "#s_data[1] is the \"stacked\" data + bkg\n", "plt.errorbar(cbins,s_data[1],np.sqrt(s_data[1]),fmt=\"*\",label=\"Binned Data\")\n", "pbins = np.linspace(xMin,xMax,250) # many small steps for smooth curve\n", "plt.plot(pbins,fitFunc(pbins,A,mu,sigma,B),label=\"Function before the fit\")\n", "plt.legend()\n", "plt.xlabel('x')\n", "plt.ylabel('Number of events per bin')\n", "plt.title('Data and background')\n", "plt.tight_layout() # make space for labels and titles\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Now do the fit then print and plot the results" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A = 42.626871 +- 5.481169\n", "mu = 438.965524 +- 2.707316\n", "sigma = 22.870076 +- 2.521536\n", "B = 9.575476 +- 0.604269\n", "Covariance matrix\n", " [[ 3.00432152e+01 8.88937936e-01 -8.11948504e+00 -1.42304384e-01]\n", " [ 8.88937936e-01 7.32956217e+00 -4.73054726e-01 -6.72937724e-04]\n", " [ -8.11948504e+00 -4.73054726e-01 6.35814304e+00 -3.56470507e-01]\n", " [ -1.42304384e-01 -6.72937724e-04 -3.56470507e-01 3.65141177e-01]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcVNWd///Xpxe6EZqtoRHZGgVxAWmgiRCNRklcIuL+\nM4mjxjGjiU5iMolbNBpnyE/zTaLRiRMxXxPNYOIal7hMBlmUuEAgoiiLgLK0oN00W7N008vn+8e9\n3RbQS/Vyq6qr38/Hox5VddfPKar6wzn33HPM3REREUk1GckOQEREpDFKUCIikpKUoEREJCUpQYmI\nSEpSghIRkZSkBCUiIilJCUokImb2sJnNaGLdN8zsbxGcc76ZfbOjjxu1qD4P6dyUoCRlmNk6M9tr\nZhVmtt3M3jCzb5lZXN9TMys0MzezrKhjFZHoKUFJqjnb3fOA4cBdwI3AQ8kNqetRkpdUoAQlKcnd\nd7j788DFwOVmNgbAzM4ys7fNbKeZbTSzn8Ts9lr4vN3MdpnZFDM7wszmmlm5mW0xs0fNrE9T5zWz\ne8Pj7jSzJWb2hZh1PzGzJ8zsD2Et730zK45ZP97M/hGuexzIbaGYZmb/aWY7zGylmU2NWXGFma0I\nj/WhmV19wI7nmNnSMM61ZnZGIwcfZGbvmtkPw/cjzOy18JivmNn9ZjYrXFdf+7zSzDYAc8Pl08Ny\nbg+bD4+OOb6b2ciY9w1Nmmb2RTMrMbMfmFmpmW02sytits03s+fD+BcBR7TwWUkXpAQlKc3dFwEl\nQH2i2A1cBvQBzgK+bWbnhutOCp/7uHtPd38TMOBO4DDgaGAo8JNmTvl3oAjoB/wReNLMYhPNdOCx\n8PzPA78GMLNuwLPAf4f7Pglc0ELxjgc+BPoDtwN/NrN+4bpSYBrQC7gCuMfMJoTn+hzwB+D6MI6T\ngHWxBzazQuBV4Nfu/otw8R+BRUB++Blc2khMJxN8Tqeb2ZHAn4DvAQOAl4C/hGWNx6FAb2AwcCVw\nv5n1DdfdD1QCg4B/Dh8i+1GCks5gE8Effdx9vrsvc/c6d3+X4A/oyU3t6O5r3H22u1e5exlwdwvb\nz3L3cnevcfdfAjnA6JhN/ubuL7l7LUEyGhcunwxkA79y92p3f4og2TWnNGb7x4FVBEkXd3/R3dd6\n4FXgf/ksSV8J/C4sV527f+zuK2OOewwwH7jd3R8EMLNhwCTgNnff5+5/I0iwB/qJu+92970EtdcX\nw/NUA78AugOfb6Fc9aqBfw/L9xKwCxhtZpkEyfu28FzvAY/EeUzpQpSgpDMYDGwFMLPjzWyemZWZ\n2Q7gWwQ1kEaZWYGZPWZmH5vZTmBWC9v/IGxa22Fm2wlqALHbfxLzeg+QG16vOQz42PcffXl9C+Vq\nbPvDwjjONLO3zGxrGMdXYuIYCqxt5riXAB8DT8UsOwzY6u57YpZtbGTf2GWHxZbB3evC9YObLdVn\nyt29Jub9HqAnQW0s64BztfRZSRekBCUpzcwmEfxBrO+C/EeC//kPdffewAMEzXgAjQ3Nf2e4/Dh3\n7wX8U8z2B57rCwSdMv4/oK+79wF2NLX9ATYDg80sdtthLezT2PabzCwHeJqgxjIwjOOlmDg20vw1\nm58AW4A/hrWV+vj6mdkhMdsNbWTf2M9wE0FnFSC4YBbu83G4aA8Qe7xDm4kpVhlQc8D5W/qspAtS\ngpKUZGa9zGwawfWeWe6+LFyVR1ATqAyvxXw9ZrcyoA44PGZZHkHT0nYzG0xw3aYpeQR/OMuALDO7\njeAaUDzeDPf9rpllmdn5wOda2Kcg3D7bzC4iuPbzEtCNoGmxDKgxszOB02L2ewi4wsymmlmGmQ02\ns6Ni1lcDFwE9gP82swx3Xw8sBn5iZt3MbApwdgvxPQGcFZ4nG/gBUAW8Ea5fCnzdzDLDThpNNp3G\nCptH/xzGcoiZHQNcHs++0rUoQUmq+YuZVRDUEm4huGZ0Rcz6a4B/D7e5jeCPKABh89VPgdfDXmeT\ngTuACQQ1oRcJ/jA25a/Ay8AHBE1OlTTeDHYQd98HnA98A9hGcP2muXMBLARGEdR2fgpcGF7/qgC+\nG5ZtG0ESbrheFHYcuQK4JyzXq8TUdA6IpwD4nQX3kl0CTAHKgRnA4wQJp6kyrSKocf5nGOPZBLcB\n7As3uS5ctj089rMtlDfWvxI0930CPAz8vhX7ShdhmrBQpGsKu8KvdPfbkx2LSGNUgxLpIsxskgX3\nhWWETXLn0Lpaj0hC6W5xka7jUIJmx3yCe8u+7e5vJzckkaapiU9ERFKSmvhERCQldYomvv79+3th\nYWGywxARkQ6wZMmSLe4+oKXtOkWCKiwsZPHixckOQ0REOoCZxTVyiJr4REQkJSlBiYhISlKCEhGR\nlNQprkGJSNdWXV1NSUkJlZWVyQ5FWiE3N5chQ4aQnZ3dpv2VoEQk5ZWUlJCXl0dhYSH7DwAvqcrd\nKS8vp6SkhBEjRrTpGGriE5GUV1lZSX5+vpJTJ2Jm5Ofnt6vWqwQlIp2CklPn095/MyUoERFJSUpQ\nIl3QxTPf5OKZbyY7jE4lMzOToqKihse6detYvHgx3/3udwGYP38+b7zxRsP2zz77LMuXL2/1eXr2\n7Bn3+TvK9u3b+a//+q+G95s2beLCCy/ssOO3lTpJiIjEoXv37ixdunS/ZYWFhRQXFwNBgurZsyef\n//zngSBBTZs2jWOOOSay83eU+gR1zTXXAHDYYYfx1FNPRXKu1oi0BmVmfczsKTNbaWYrzGyKmfUz\ns9lmtjp87htlDCIiUZk/fz7Tpk1j3bp1PPDAA9xzzz0UFRXx6quv8vzzz3P99ddTVFTE2rVrWbt2\nLWeccQYTJ07kC1/4AitXrgTgo48+YsqUKUyaNIkf//jHrTr/ww8/zL/+6782vJ82bRrz588HgprY\nLbfcwrhx45g8eTKffvopAJ9++innnXce48aNY9y4cbzxxhvcdNNNrF27lqKiIq6//nrWrVvHmDFj\ngKCDyhVXXMHYsWMZP3488+bNazj3+eefzxlnnMGoUaO44YYb2vtxHiTqGtS9wP+4+4Vm1g04BPgR\nMMfd7zKzm4CbgBsjjkNE0sXLN8Enyzr2mIeOhTPvanaTvXv3UlRUBMCIESN45plnGtYVFhbyrW99\ni549e/LDH/4QgOnTpzNt2rSGprKpU6fywAMPMGrUKBYuXMg111zD3Llzue666/j2t7/NZZddxv33\n39+m8zdm9+7dTJ48mZ/+9KfccMMN/Pa3v+XWW2/lu9/9LieffDLPPPMMtbW17Nq1i7vuuov33nuv\noYYW23xYH9OyZctYuXIlp512Gh988AEAS5cu5e233yYnJ4fRo0fzne98h6FDhzYbV2tElqDMrBdw\nEvANAHffB+wzs3OAL4abPQLMRwlKRFJce5rYdu3axRtvvMFFF13UsKyqqgqA119/naeffhqASy+9\nlBtvbPzPYWvP361bN6ZNmwbAxIkTmT17NgBz587lD3/4AxBc1+rduzfbtm1r8jh/+9vf+M53vgPA\nUUcdxfDhwxsS1NSpU+nduzcAxxxzDOvXr+8cCQo4HCgDfm9m44AlwHXAQHffDODum82soLGdzewq\n4CqAYcOGRRimiHQqLdR0UlFdXR19+vRpMsG0tTt2VlYWdXV1De9j7znKzs5uOG5mZiY1NTVtOkdz\nk9rm5OQ0vG7POZoS5TWoLGAC8Bt3Hw/sJmjOi4u7P+juxe5ePGBAi9OGiIgkVV5eHhUVFY2+79Wr\nFyNGjODJJ58Egj/677zzDgAnnHACjz32GACPPvpoq85ZWFjI0qVLqaurY+PGjSxatKjFfaZOncpv\nfvMbAGpra9m5c+dBscc66aSTGuL64IMP2LBhA6NHj25VnG0VZYIqAUrcfWH4/imChPWpmQ0CCJ9L\nI4xBRCQhzj77bJ555hmKiopYsGABX/3qV/n5z3/O+PHjWbt2LY8++igPPfQQ48aN49hjj+W5554D\n4N577+X+++9n0qRJ7Nixo1XnPOGEExgxYgRjx47lhz/8IRMmTGhxn3vvvZd58+YxduxYJk6cyPvv\nv09+fj4nnHACY8aM4frrr99v+2uuuYba2lrGjh3LxRdfzMMPP7xfzSlK1lz1rd0HN1sAfNPdV5nZ\nT4Ae4arymE4S/dy92e4fxcXFrgkLRTpO/T1Qj189JcmRxGfFihUcffTRyQ5D2qCxfzszW+LuxS3t\nG3Uvvu8Aj4Y9+D4EriCotT1hZlcCG4CLmtlfRES6qEgTlLsvBRrLklOjPK+IiHR+GupIRERSkhKU\niIikJCUoERFJSUpQIpKWNGJ756cEJSISh/rpLsaNG8eECRMaptZI1NQUhYWFbNmypdHlY8eOZezY\nsRxzzDHceuutDcMoNeXA6TVSlRKUiKSlfTV1LN+0k9KKtk85Hqt+LLx33nmHO++8k5tvvhlIjakp\n5s2bx7Jly1i0aBEffvghV111VbPbK0GJiCTRx9v3UlFVw32vrO7wY+/cuZO+fYOZgmKnpmhuCoqm\npr8oKyvjggsuYNKkSUyaNInXX38dgPLyck477TTGjx/P1Vdf3eyYeLHneOCBB3j22WfZunUru3bt\nYurUqUyYMIGxY8c2jF5x4PQaTW2XbJqwUETSyuhbX6aq5rMBVGct3MCshRvIycpg1Ywz23zc+uku\nKisr2bx5M3Pnzm10u6amoGhq+ovrrruO73//+5x44ols2LCB008/nRUrVnDHHXdw4oknctttt/Hi\niy/y4IMPxhVn/bh/q1evZuLEiTzzzDP06tWLLVu2MHnyZKZPn37Q9Bo1NTWNbtfWQWw7ihKUiKSV\nBTecwoyXVvDCO5uoc8jNzuD0Yw/llrPaN1RS7HQXb775JpdddhnvvffeQds1NQVFU9NfvPLKK/tN\nDb9z504qKip47bXX+POf/wzAWWed1VBji0d9bcvd+dGPfsRrr71GRkYGH3/8cUPN7cDtG9vu0EMP\njfucUVCCEpG0UtArl7ycLOoczKCqpo68nCwK8nI77BxTpkxhy5YtlJWVHbSuqSkompr+oq6ujjff\nfJPu3bsfdKy21GAqKipYt24dRx55JI8++ihlZWUsWbKE7OxsCgsL95uSo1682yWarkGJSNrZsquK\ngrwcjh3Ui0uOH07ZruZ7tbXWypUrqa2tJT8/v93HOu200/j1r3/d8L6+lhY7zcXLL7/c7KSC9Xbt\n2sU111zDueeeS9++fdmxYwcFBQVkZ2czb9481q9fDxw8NUhT2yWbalAiknZmXlrccA/UjHPHdMgx\nY6dcd3ceeeQRMjMz233c++67j2uvvZbjjjuOmpoaTjrpJB544AFuv/12vva1rzFhwgROPvnkZidu\nPeWUU3B36urqOO+88/jxj38MwCWXXMLZZ59NcXExRUVFHHXUUQD7Ta9x5plncuONNza6XbJFOt1G\nR9F0GyId67z7X2dN6S7m/PDkDm36ikpbptvobFOKpKv2TLehJj6RLijKLtip4vGrpyg5dXJq4hPp\nQqLqgi0SBdWgRLqQBTecwvSiw8gIO4flZmdwTtFhLLjxlOQGFofOcDlC9tfefzMlKJEuJBFdsKOQ\nm5tLeXm5klQn4u6Ul5eTm9v275aa+ES6mPou2AV5ORQN60tZB41VF6UhQ4ZQUlLS6H1Hkrpyc3MZ\nMmRIm/dXghLpYqLogh217OxsRowYkewwJMHUxCciIilJNSiRdFVXC0sehg1vwrDJMOmbyY5IpFWU\noETS1YK7Yd4M6N4Xlj0Ju8rglJuTHZVI3NTEJ5KOShbD/Dth7EVw/VoY9zV49S74+B/JjkwkbkpQ\nIuno1f8DPfrDWb+EjEw482eQ0wtevzfZkYnETQlKJN3sKIE1s2H8pZAbzEtEbm+YdCWseB7K1yY3\nPpE4RZqgzGydmS0zs6Vmtjhc1s/MZpvZ6vA5/lm4RKRlb88Cr4MJl+6//PhvAwZLH01KWCKtlYga\n1CnuXhQzcu1NwBx3HwXMCd+LSEdwDxLQ4adA38L91+UNhOGfh5UvJSU0kdZKRhPfOcAj4etHgHOT\nEINIetryAWzfAMdMb3z9UWdB2QoG1mxKbFwibdBigjKzE8KmuA/M7EMz+8jMPozz+A78r5ktMbOr\nwmUD3X0zQPhc0LbQReQga+YEz0dMbXz96K8AMKnyzQQFJNJ28dwH9RDwfWAJUNvK45/g7pvMrACY\nbWYr490xTGhXAc3OJCkiMda8AvmjoO/wxtf3HQ4DxzBh21u80POCxMYm0krxNPHtcPeX3b3U3cvr\nH/Ec3N03hc+lwDPA54BPzWwQQPhc2sS+D7p7sbsXDxgwIK7CiHRp1Xth/eswsonaU73Dv8iofSvJ\n8n0JCUukreJJUPPM7OdmNsXMJtQ/WtrJzHqYWV79a+A04D3geeDycLPLgefaGLuIxNq4EGoqm27e\nqzdsCt2o5ojq9J1NV9JDPE18x4fPsfPHO3BqC/sNBJ4xs/rz/NHd/8fM/g48YWZXAhuAi1oXsog0\nauMiwGDY8c1vN2wyAEftey/6mETaocUE5e5tmmrT3T8ExjWyvBxo4b94ItJqJX+HAaM/uzm3KT36\n83HmUEbvez8xcYm0UZMJysz+yd1nmdm/Nbbe3e+OLiwRaRX3IEEdNS2uzVd0G8PkygVQVwcZGlBG\nUlNz38we4XNeEw8RSRXla2HvNhgyKa7NV3U7lp6+C7asijgwkbZrsgbl7jPD5zsSF46ItEnJouA5\nzgS1NntU8GLTUig4OqKgRNonnht1Dzezv5hZmZmVmtlzZnZ4IoITkTiVLIZueTDgqLg235Q1hErL\ngc3vRByYSNvF0/j8R+AJYBBwGPAk8KcogxKRVvrkXRg0Lu7rSW6ZrM86AjYvjTgwkbaL59ts7v7f\n7l4TPmYRdDMXkVRQVwufvg+Hjm3Vbh9mj4TN7wYdJURSUJMJKpwWox/Bjbo3mVmhmQ03sxuAFxMX\noog0a+uHUL2nVQnq8auncOaXz4Dq3bBV80NJamruPqglBDUlC99fHbPOgf+IKigRic/FM99kyt5X\n+R60ugbFoPA2xU1Lof+ojg5NpN2a68U3IpGBiEjbFFavhYzsuDtINBhwFGTmwCfvwHEa0EVSj+7Q\nE+nkhtd8GCSbrG6t2zEzCwYcCaUroglMpJ2UoEQ6ucLqta1v3qtXcKwSlKSsZhOUBYYmKhgRaZ28\nuh30rdsGA49t2wEKjoadHwejUIikmGYTlLs78GyCYhGRVhpSvT54UdDK60/1Co4JnkvjnktUJGHi\naeJ7y8ziGz9FRBJqaE2YoFrbQaLewPoEpZHNJfXEMx/UKcC3zGwdsJug27m7+3FRBiYiLRtcs4E9\ndgiH9BrctgP0Ggw5vXQdSlJSPAnqzMijEJE2GVqznpKsYRxp1vLGjTELrkMpQUkKarGJz93XA0OB\nU8PXe+LZT0SiN6RmPSVZw9t3ECUoSVHxjGZ+O3AjcHO4KBuYFWVQIhKH3VvoXbeDje1NUPmjYO9W\n2F3eMXGJdJB4akLnAdMJrj/h7pvQhIUiyRfWekqy25mg+h8ZPJevbmdAIh0rngS1L+xu7gBm1qOF\n7UUkEbZ8AMDHWcPad5z+I8PjKUFJaoknQT1hZjOBPmb2L8ArwG+jDUtEWlS+lkrLYWtGfvuO02c4\nZHZrSHgiqaLFXnzu/gsz+zKwEzgSuM3dZ0cemYg0r3wNn2QOxq2dfZYyMqHfEVC+pmPiEukg8XQz\nB1gGdCdo5lsWXTgiErfyNWzOauP9TwfqP1I9+STlxNOL75vAIuB84EKCkSX+OerARKQZtdWwbV0H\nJqgjYdu64LgiKSKeGtT1wHh3Lwcws3zgDeB3UQYmIs3Yth68ls2ZQzrmePmjoK4Gtn4UTMEhkgLi\nabwuASpi3lcAG6MJR0TiEl4v6tAaFKiruaSUeGpQHwMLzew5gmtQ5wCLzOzfANz97uZ2NrNMYDHw\nsbtPM7MRwGNAP+AfwKXuvq8dZRDpesIE9RGDWb5pJ6UVlRTk5bb9eOpqLikonhrUWoIpNzx8/xyw\nmeBm3Xhu2L0OiL36+jPgHncfBWwDrow7WhEJlK+B7v1YvTOTiqoa7nulnYkltzf0KFCCkpQSTzfz\nO9p6cDMbApwF/BT4NzMz4FTg6+EmjwA/AX7T1nOIdEWj3ziFKv8yUAXArIUbmLVwAzlZGaya0cbx\nnfsfqSY+SSlRD/r6K+AGoC58nw9sd/ea8H0J0GgjupldZWaLzWxxWVlZxGGKdC4L+v6U6fklZISD\nmOdmZ3BO0WEsuPGUth+0/0jdrCspJbIEZWbTgFJ3XxK7uJFNvZFluPuD7l7s7sUDBgyIJEaRTqlq\nFwV7VpHX4xDqPJgxo6qmjrycrHZehzoymPpdg8ZKioj3Rt22OAGYbmZfAXKBXgQ1qj5mlhXWooYA\nmyKMQST9bP0QgC3eh4K8HAryciga1peyisr2HTd/VPC85QPoMaWdQYq0Xzw36v4fM+tlZtlmNsfM\ntpjZP7W0n7vf7O5D3L0Q+Cow190vAeYR3PALcDlBpwsRiVfYg2/mBSMY0b8HPXKymHHuGGZeWty+\n4/YPE5SuQ0mKiKeJ7zR33wlMI7hmdCTBzbttdSNBh4k1BNekHmrHsUS6nvK1wXO/wzv2uH2GQWaO\nrkNJyoiniS87fP4K8Cd332qtnF7a3ecD88PXHwKfa9UBROQz5Wug1xDodkjHHjcjM0h69QlQJMni\nSVB/MbOVwF7gGjMbALSzsVtE2qx8DeQfEc2x849QDUpSRjxNfLcDU4Bid68G9hDMsCsiieYeXCPK\nHxnN8fNHBuPx1dVGc3yRVognQb3p7tvcvRbA3XcDL0cblog0as9WqNwRbYKqq4btG6I5vkgrNNnE\nZ2aHEtxE293MxvPZPUy9gA5u/BaRuNRPKhhlgoLgOlS/EdGcQyROzV2DOh34BsG9SrEDwlYAP4ow\nJhFpSkOCiuoaVH2CWgOjvhTNOUTi1GSCcvdHgEfM7AJ3fzqBMYlIU8rXQEYW9BkezfF79IecXpr+\nXVJCPL34XjCzrwOFsdu7+79HFZSINKF8DfQdAZkRDQJjFtTOlKAkBcTzLX8O2AEsoX7oZBFJjvK1\n0V1/qpc/EjYsjPYcInGIJ0ENcfczIo9ERJpXVwdb18IR7RixPB75I2HZU1C9F7K7R3sukWbE0838\nDTMbG3kkItK8nR9DTWVialB4cD+USBLFk6BOBJaY2Soze9fMlpnZu1EHJiIHiLqLeb36HoK6DiVJ\nFk8TXxun5xSRDpWoBNUvTFBbNSafJFeLNSh3Xw8MBU4NX++JZz8R6WDlayG7B+QdGu15cntBz4Gq\nQUnSxTMf1O0EU2TcHC7KBmZFGZSINKJ+kNhWzibQJvkjNaq5JF08NaHzCAaH3Q3g7puAvCiDEpFG\nlK+Jvnmvnu6FkhQQT4La5+4OOICZ9Yg2JBE5SM0+2L4+gQlqJOwug73bE3M+kUbE00niCTObCfQx\ns38B/hn4bbRhich+tq0DrzsoQT1+9ZRozld/nq1rYfDEaM4h0oIWE5S7/8LMvgzsBEYDt7n77Mgj\nE5HPJKoHX736nnzlSlCSPC0mKDP7PvCkkpJIEjUkqMMTc75+IwDTdShJqniuQfUC/mpmC8zsWjMb\nGHVQInKA8jVwSH/o3jcx58vKgT7DlKAkqeK5D+oOdz8WuBY4DHjVzF6JPDIR+Ux9F/NEUldzSbLW\n3HBbCnwClAMF0YQjIo3ashryRyX2nPUJyj2x5xUJxXOj7rfNbD4wB+gP/Iu7Hxd1YCISqtwBu0uh\nf4I6SNTLHwn7KmBXaWLPKxKKp5v5cOB77r406mBEpBFb6jtIJLoGFTNobJ4uPUvixXMN6iagp5ld\nAWBmA8xsROSRiUigfHXw3D8JTXygjhKSNJGNxWdmuWa2yMzeMbP3zeyOcPkIM1toZqvN7HEz69ae\nAoikvS2rwTKDqd4TqfcQyMxRgpKkiXIsviqCEdDHAUXAGWY2GfgZcI+7jwK2AVe2JXCRLqN8NfQt\nhKwE/18uIxP6Ha6efJI0kY3F54Fd4dvs8OHAqcBT4fJHgHNbFbFIV7NlTeKb9+pp0FhJongS1IFj\n8b1CnGPxmVmmmS0l6KI+G1gLbHf3mnCTEmBw68MW6SLq6oLx8BI1xNGB8o+AbR9BXW1yzi9dWqRj\n8bl7LVBkZn2AZ4CjG9ussX3N7CrgKoBhw4bFczqR9LNjI9RUJrEGNRJq9wVx9C1MTgzSZcXTzZww\nIbV5LD533x7eSzWZoCaWFdaihgCbmtjnQeBBgOLiYt0pKF1TfQ++RHcxrxfbk08JShIssqnbw+7o\nfcLX3YEvASuAecCF4WaXA89FFYNIp7clSV3M6zUkKHWUkMSLqwbVRoOAR8wskyARPuHuL5jZcuAx\nM5sBvA08FGEMIp3bltWQ0xt6DEjO+XsMgJxe6ighSdFkgjKzOe4+1cx+5u43tvbA7v4uML6R5R8C\nn2vt8US6pPLVwRBHZsk5v5l68knSNFeDGmRmJwPTzewxYL9fiLv/I9LIRCToYj7ipOTGkD8SNi5M\nbgzSJTWXoG4DbiLoyHD3Aevq72cSkahU7YKKTYkfJPZA+SNh2VNQXQnZucmNRbqUJhOUuz8FPGVm\nP3b3/0hgTCICMbPoJqmDRL3+owDnhplPsz77cB6/ekpy45EuI577oP7DzKYD9e0M8939hWjDEpGG\nBNX/yOTGMSC4fXFozTrWZydoynkR4hss9k7gOmB5+LguXCYiUdqyGrBgPLxkyh8JGVkMqV6f3Dik\ny4mnm/lZQJG71wGY2SME3cNvbnYvEWmfshVBckr2dZ+sbtDvCIZUKEFJYsV7o26fmNe9owhERA7w\n6XIoaGx0sCQoOJqhNUpQkljxJKg7gbfN7OGw9rQE+P+jDUuki6uuDAaJLTgm2ZEECo6moPYTunll\nsiORLiSeThJ/CsfRm0RwL9SN7v5J1IGJdGlbPgCvS50a1ICjyMAZXLMx2ZFIFxLvYLGbgecjjkVE\n6pWuCJ5TqAYFqKOEJFRkg8WKSPwunvkmF89887MFpcshIzsYZigV9DucGrIYUrMh2ZFIF6IEJZKK\nSlcE9z81m5jQAAAVnklEQVRlZic7kkBmNpuyhqijhCRUswnKzDLM7L1EBSMiodLlMDBFmvdCJVnD\nGaIEJQnUbIIK7316x8w0pa1IolTuDGawTZUOEqGNWcMYWPsJ7Nud7FCki4ink8Qg4H0zWwQ0fDPd\nfXpkUYl0ZWUrg+dU6SARKskaHrwoWwWDJyQ3GOkS4klQd0QehYh8pnR58JxqNajs+gS1UglKEqLF\nThLu/iqwDsgOX/8d0FxQIlEpXQHZPaB3arWsb/RDuajqNko3avJCSYx4Bov9F+ApYGa4aDDwbJRB\niXRppcuh4CjISK1Otht37GOxH8l9K3okOxTpIuJp4ruWYIr2hQDuvtrMCiKNSqQr+3Q5jD4j2VE0\nGH3ry1TV1IXvMpi19Shm3fQiOVkZrJpxZlJjk/QWz3/Rqtx9X/0bM8simFFXRDrarjLYswUKjk12\nJA0W3HAK04sOI8OC97lUcc6YfBbceEpyA5O0F0+CetXMfgR0N7MvA08Cf4k2LJEuqvT94DmFOkgU\n9MolLyeLOodM6qgim7ya7RTkafp3iVY8CeomoAxYBlwNvATcGmVQIl3WpqXB86BxyY3jAFt2VVGQ\nl8OEQd24JPMVyrZvT3ZI0gXEM5p5XTjNxkKCpr1V7q4mPpEobH4n6L13SL9kR7KfmZcWN4wVOKP/\nXBhSAXw1uUFJ2msxQZnZWcADwFqC6TZGmNnV7v5y1MGJdBX7aupYU7qL0j1rKBh0XLLDad6g42Dz\nu8mOQrqAeJr4fgmc4u5fdPeTgVOAe6INS6Rr+Xj7Xiqqarhvy0Q4rCjZ4TTv0OOgfA1U7Up2JJLm\n4klQpe4ee2feh0BpRPGIdCmjb32ZwptepLSiCoBZtV+m8KVjGH1rCjdQDBoHOHyyLNmRSJprMkGZ\n2flmdj7BOHwvmdk3zOxygh58f2/pwGY21MzmmdkKM3vfzK4Ll/czs9lmtjp87tthpRHpZDplF+7B\nE4Pnj5ckNw5Je83VoM4OH7nAp8DJwBcJevTFk1RqgB+4+9HAZOBaMzuGoFfgHHcfBcwJ34t0SbFd\nuLOpCbpw9+iR2l24ew6APsPg48XJjkTSXJOdJNz9ivYcOJwmfnP4usLMVhAMk3QOQaIDeASYD9zY\nnnOJdGb1Xbjvrv45/9PzXMp2DUp2SC0bPFE1KIlcPGPxjTCzu83sz2b2fP2jNScxs0JgPEFX9YFh\n8qpPYo0Om2RmV5nZYjNbXFZW1prTiXSIg6Zhj8jMS4sp6lvFibzNjM9nMvPS4sjP2W6DJ8L2DcHI\nFyIRiWcsvmeBhwiuPdW1sO1BzKwn8DTwPXffaWZx7efuDwIPAhQXF+u+K0lro6rDOaCGfC65gcQr\n9jpUCo0bKOklngRV6e73teXgZpZNkJwedfc/h4s/NbNB7r7ZzAahHoEijNq3khqyyEqxESSaNKgI\nLBNK/q4EJZGJp5v5vWZ2u5lNMbMJ9Y+WdrKgqvQQsMLd745Z9Txwefj6cuC5VkctkmZGVa/go+wj\nIDuFO0fE6nZIcMPuhreSHYmksXhqUGOBS4FT+ayJz8P3zTkh3G+ZmYUDjPEj4C7gCTO7EtgAXNTa\noEXSSm01R1R/wNzuZzAq2bG0xrDPw+KHoKYKsnKSHY2koXgS1HnA4bFTbsTD3f9GMDRSY6a25lgi\naW3T2+R6FSu6jeUryY6lNYZPgbfuh01vw7DJyY5G0lA8CeodoA+6ViQSjXULAFjRbUySA2ne41dP\n2X/BsPD9+jeUoCQS8SSogcBKM/s7UFW/0N2nRxaVSFey7m9syCqkIrNPsiNpnR79of+RsCH6rvjS\nNcWToG6PPArpcPX37xz0v16JXKs++9pq2LCQ5d06aav38BNg2VNBOTKzkx2NpJl45oN6NRGBiHRJ\nm96G6t283yPFp9hoyhGnwJLfB/dDqZlPOlg8I0lUmNnO8FFpZrVmtjMRwYmkvbVzAWNFzthkR9I2\nI04CywjLIdKxWkxQ7p7n7r3CRy5wAfDr6EMTSa59NXUs37ST0orK6E6yejYMnkhFRu/ozhGl7n2D\nUSWUoCQC8dyoux93f5aW74GSTqitY88lasy6RGuYRPCV1dGcYPeWoGls1Gk8fvWUznu98IhTg3Ls\n3ZbsSCTNxDPl+/kxbzOAYoIbdUXS0uhbX6aq5rNhJ2ct3MCshRvIycpg1YwzO+5Ea+cCDqO+1HHH\nTIaRX4JXfwZr5sDYC5MdjaSReGpQZ8c8TgcqCKbMEElLB00imJ3BOUWHdfwkgh/8FQ7pD4PGd+xx\nE21wMfQcCCv+kuxIJM3E04uvXfNCiXQ2sZMImkFVTR15OVkdO4lgdWWQoI49FzJa3dKeWjIyYPRX\nYNmTQbk6y3iCkvKaTFBmdlsz+7m7/0cE8YikhPpJBAvyciga1peyju4osXYO7KuAY8/r2OMmy9HT\ngu7mH87vdKObd4Z7BtsaY2coW3Oa+6/b7kYeAFeiGXAlzc28tJgR/XvQIyeLGeeO6fhJBN9/Brr3\nC7ppp4PCkyCnNyx/NtmRSBppMkG5+y/rHwQTB3YHrgAeAw5PUHzSRgnpIi2NavGz37cHVr0c1DrS\nZfSFrG5w7Dmw/Hmo2pXsaCRNNNv4bWb9zGwG8C5Bc+AEd7/R3TVwbIqLvIu0NKnFz375s7BvFxx3\ncWIDi1rRJVC9G1Y8n+xIJE00dw3q58D5BLWnse6elv8t6uxttAdKWBdpOUjcn/2SRyB/ZDCOXToZ\nejz0OxyW/hGKvt6wON1+Y5I4zdWgfgAcBtwKbIoZ7qhCQx2lrvZ0kW5rs6CaEwNxffalK2HjWzDh\nsqCLYDoxg/H/FEwfUroy2dFIGmjuGlSGu3c/YKijXvXvExmkxK89XaTb2iyo5sRAXJ/9wt9AZjcY\n9/WmD9SZTfgGZOXCW/+V7EgkDcQz3YZ0Mq3tIt3WZkE1Jx6s2c9+5+aw+esS6DkgeUFGqUc+jPta\nUM5Tf5y+5ZSE6OR3CHYuiRqzrrVdpNvaLJiwERc6kWY/+7fuh7oaOOG7yQswEaZcC3XV8Pqvkh1J\npNJ1DMpUogQlbW4WTMiIC+li+wZY9FsYe1HQkSCd9R8V1KIWPQjbNyY7mhZ1hmuoXfX6sBKUAJ81\nTR07qBeXHD+csl1Vke7X5cy+DTCY2twALWnkizcDBnPuSHYkLeoM11C76vVhXYMSIGiaqm+umHHu\nmMj361JWzw5Gjjj5Jug9JNnRJEafoXDi9+DVnzG+7xjezj0+2REdpDNcQ+3q14e7fA2qLVXgRLY9\nq517f4n8PDpkjqZdZfDst6HgWDjx+x0TWGfxhR9CwTGct/1hVm3ampBmptZ8PzrDNdRkXB9Opb85\nXT5BdfYqsKSw6kp44jKo3AkX/LbrjfKd1Q3Om8nD+05lR1Ud981eleyI9tMZrqF29evDXbaJL12q\nwJKasnwf/PmbsOENuOAhGHhsskNKuM9+Y1MBmLWohFmLSlLqNxb5qPUdoK0xdoaytSSyGpSZ/c7M\nSs3svZhl/cxstpmtDp/7RnX+lnSG6r10Tj3rdnLT1tuCCfxOv7PLzjJ70G+MKs7ptZoF3ylKbmAx\nIh+1vgO0NcbOULaWRNnE9zBw4MQwNwFz3H0UMCd8nxTJqAJ39i6f0gJ3WPkiPy/7Nkfvew/O/Q1M\nuSbZUSXNQb8xupG3ZwMFj34JVrwQfF6dmH7P0YssQbn7a8DWAxafAzwSvn4EODeq88cj0V2kdb0r\nTe3bA8uegodOg8e+zu6MHtzS/1f7DZjaVe33G5tcSNmwr0BuL3j8Evi/Xwo+t317kh1mm+j3HL1E\nX4Ma6O6bAdx9s5kVJPj8+0lUF+lkXO/SyNERqKuDPVtg5ybYtg42L4VNS2HDW1CzF/oWwll3M3TC\nZfwsXeZ5aqeDf2NjoPaL8I+H4Y3/hKevhKzuMKQYhkyCw4qgzzDoPQwO6ZeSA+rq+nXimEdYzTaz\nQuAFdx8Tvt/u7n1i1m9z90avQ5nZVcBVAMOGDZu4fv36tgfyp69DTWPVcOedjdsBGDe0zwGrmv5c\n5m+o5ud7z+b3R8ynIGtvi/uVVh/CjNLj+evOEVSRTa5Vc3reR9wy8M1g/2bOtaCkhp/tPYffjZgb\n17liy9b44qb3WbF5BwBHD2pkLOBm9lv1STC4/ehD81qxXzOxN3Outz5xflp5AQ8NeZmCrAP+593B\n58LroHov7NsdPnaB1362PiMr6Pww9Hg46qxgVtmMLt8x9iBNTrdRVwvrX4eVL8LGhfDJsmAoqAYG\n3XpCTs/gOSs3+HwtEzIyY54zGh7zNtTwi73T+P0Rr1GQfcDvpQlN/g1oQml1d2Z8ejx/3TE8/D3X\ncHqvddwycFHc52ytpSVBjEVDDoixhQT+dli28XGWDWDO+mp+uecrPDzyb42X55D+cP7MuI/XGDNb\n4u4tXhRLdA3qUzMbFNaeBgFNTnzo7g8SzEVFcXFx+7Jo5Y4mEhR09/CPXFVmI2sb/8d/qupEltcO\n5r5PjmPGwPkt7ldgVeRRyT4y6UY1VZ5FnlVSYDuhtvF96j1W9Xnerx3CfWUTmXHoa3HHCNb4l7eZ\n73OtZeFYM7O8Nr7zPusW7JfdvVX7Nf/janzdw9UTea92GPdVfJEZhzV2r0bHnQszyD4EuvUI/kB2\n6wE9B0KvQcENtwOO7npdxztSRmYw5X39tPf79sCWD2BHCb9/eQF5dRVcOKY3VFUEj9p9QVLz2s+e\n3cPX+8DreLrhtzmGGQVz4wqj+b8BByuggjzfHfN7ziTPd1FQ+0nM77lpq0uDqfVGFfSM63wAPep2\nBy8qY7+rLf9Z7FkXTuO3N/4/oc9WnsSK2sHct/loZgyYc/AGGYlrHUh0DernQLm732VmNwH93P2G\nlo5TXFzsixcvjiTG1kymdmDVvl48Vfur/3sxb2/Yvl+Xz+Z61bTnXOlIn0fn1ZYJC1u7T3u+H22J\nr7W/5/aer62TPibq71trxVuDirKb+Z+AN4HRZlZiZlcCdwFfNrPVwJfD951Ge7qmJ2qE8XSlz0Oa\nk+jvRzp04T5QKv7GImvic/evNbFqalTnjFoiu6any53gHUWfhzRH34/2S8XPsMuOJNFWibw7Ox3u\nBG9Oa5st2vN5tLWJRDqPtn4/OsN3oq0xtna/VPubowTVSokcvVsjhe9Pn4c0R9+P9ku1z1B9YkUk\n5aTzKA3pXLaO1uUTVIdMqSAiTWrLbyydR2lI57J1NDXxiUjKSOdRGtK5bFFRgkog1dT2t6+mjjWl\nuyitqIy8p1AizyVtt+CGU5jx0gpeeGcTdR50dT792EO55ayjkx3aQVr7e+5MZUsVXb6JT5InkU0d\nalbpHFKxq3NHSeeyRUU1KEm4RDZ1qFml80m1rs4dKZ3LFoVIhzrqKFEOdSSJV7qzssmmjo7+32Qi\nzyUdJ53vW0vnssUr6UMdiTRFI3KISDzUxCdJoRE5RKQlSlCSFBqRQ0RaoiY+ERFJSUpQIiKSktTE\nJyIpJ517uKVz2TqaalAiIpKSlKBERCQl6UZdERFJKN2oKyIinZoSlIiIpCQlKBERSUlKUCIikpKU\noEREJCUpQYmISEpSghIRkZSkBCUiIilJCUpERFJSpxhJwszKgPXtPEx/YEsHhJOq0r18kP5lTPfy\ngcqYDjqifMPdfUBLG3WKBNURzGxxPENrdFbpXj5I/zKme/lAZUwHiSyfmvhERCQlKUGJiEhK6koJ\n6sFkBxCxdC8fpH8Z0718oDKmg4SVr8tcgxIRkc6lK9WgRESkE1GCEhGRlJQWCcrMfmdmpWb2Xsyy\nfmY228xWh899w+VmZveZ2Roze9fMJiQv8viZ2VAzm2dmK8zsfTO7LlyeFuU0s1wzW2Rm74TluyNc\nPsLMFoble9zMuoXLc8L3a8L1hcmMvzXMLNPM3jazF8L3aVNGM1tnZsvMbKmZLQ6XpcV3tJ6Z9TGz\np8xsZfh7nJJOZTSz0eG/X/1jp5l9LxllTIsEBTwMnHHAspuAOe4+CpgTvgc4ExgVPq4CfpOgGNur\nBviBux8NTAauNbNjSJ9yVgGnuvs4oAg4w8wmAz8D7gnLtw24Mtz+SmCbu48E7gm36yyuA1bEvE+3\nMp7i7kUx98qky3e03r3A/7j7UcA4gn/LtCmju68K//2KgInAHuAZklFGd0+LB1AIvBfzfhUwKHw9\nCFgVvp4JfK2x7TrTA3gO+HI6lhM4BPgHcDzBHetZ4fIpwF/D138FpoSvs8LtLNmxx1G2IQQ/7lOB\nFwBLpzIC64D+ByxLm+8o0Av46MB/h3Qq4wHlOg14PVllTJcaVGMGuvtmgPC5IFw+GNgYs11JuKzT\nCJt6xgMLSaNyhk1fS4FSYDawFtju7jXhJrFlaChfuH4HkJ/YiNvkV8ANQF34Pp/0KqMD/2tmS8zs\nqnBZ2nxHgcOBMuD3YTPt/zWzHqRXGWN9FfhT+DrhZUznBNUUa2RZp+lrb2Y9gaeB77n7zuY2bWRZ\nSpfT3Ws9aFYYAnwOOLqxzcLnTlc+M5sGlLr7ktjFjWzaacsInODuEwiafa41s5Oa2bYzli8LmAD8\nxt3HA7v5rKmrMZ2xjACE10KnA0+2tGkjyzqkjOmcoD41s0EA4XNpuLwEGBqz3RBgU4JjaxMzyyZI\nTo+6+5/DxWlXTnffDswnuNbWx8yywlWxZWgoX7i+N7A1sZG22gnAdDNbBzxG0Mz3K9KojO6+KXwu\nJbhu8TnS6ztaApS4+8Lw/VMECSudyljvTOAf7v5p+D7hZUznBPU8cHn4+nKCazb1yy8Le55MBnbU\nV1tTmZkZ8BCwwt3vjlmVFuU0swFm1id83R34EsHF53nAheFmB5avvtwXAnM9bABPVe5+s7sPcfdC\ngqaTue5+CWlSRjPrYWZ59a8Jrl+8R5p8RwHc/RNgo5mNDhdNBZaTRmWM8TU+a96DZJQx2RfhOuhC\n3p+AzUA1QTa/kqCtfg6wOnzuF25rwP0E1zeWAcXJjj/OMp5IUG1+F1gaPr6SLuUEjgPeDsv3HnBb\nuPxwYBGwhqCpISdcnhu+XxOuPzzZZWhleb8IvJBOZQzL8U74eB+4JVyeFt/RmHIWAYvD7+qzQN80\nLOMhQDnQO2ZZwsuooY5ERCQlpXMTn4iIdGJKUCIikpKUoEREJCUpQYmISEpSghIRkZSkBCUiIilJ\nCUpERFKSEpRIkpnZpHAendxwNIb3zWxMsuMSSTbdqCuSAsxsBsHIEd0Jxnq7M8khiSSdEpRICghH\njv47UAl83t1rkxySSNKpiU8kNfQDegJ5BDUpkS5PNSiRFGBmzxNMwTGCYDbSf01ySCJJl9XyJiIS\nJTO7DKhx9z+aWSbwhpmd6u5zkx2bSDKpBiUiIilJ16BERCQlKUGJiEhKUoISEZGUpAQlIiIpSQlK\nRERSkhKUiIikJCUoERFJSf8PqKPkS9aQrIMAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "popt, pcov = curve_fit(fitFunc,cbins,s_data[1],sigma=np.sqrt(s_data[1]),p0=(A,mu,sigma,B))\n", "# print the results\n", "for (name,i) in [('A',0),('mu',1),('sigma',2),('B',3)]:\n", " print(\"{:7} = {:f} +- {:f}\".format(name,popt[i],np.sqrt(pcov[i][i])))\n", "print(\"Covariance matrix\\n\",pcov)\n", "# replot data\n", "plt.errorbar(cbins,s_data[1],np.sqrt(s_data[1]),fmt=\"*\",label=\"Binned Data\")\n", "# plot fitted values\n", "plt.plot(pbins,fitFunc(pbins,popt[0],popt[1],popt[2],popt[3]),label=\"Fitted Function\")\n", "plt.legend()\n", "plt.xlabel('x')\n", "plt.ylabel('Number of events per bin')\n", "plt.title('Data and background')\n", "plt.tight_layout() # make space for labels and titles\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The background is now known to be a decaying exponential\n", "Change the code to simulate and fit the new distributions with \"scale = 200\" for the background, i.e.\n", "\n", "$bkg = B\\cdot \\exp\\left( -\\frac{x}{200}\\right)$\n", "\n", "See\n", "https://docs.scipy.org/doc/scipy-0.19.1/reference/tutorial/stats/continuous_expon.html\n", "for details of the exponetial distribution" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "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.1" } }, "nbformat": 4, "nbformat_minor": 2 }