{ "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": 2, "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": 3, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAGmZJREFUeJzt3Xm0JlV57/HvTwZxAFqg0WaywaDR6wDeFiEmDhiNCIIa\nNc5IiKjRG4zeBDS5DGqWeuOY6EokwYgiCKIXEPEqImCMyowgQa5gWkXQbpFJELDluX9UHXhpz1B9\n6Dqnznu+n7Xe9Vbtmp59OH0e9q5du1JVSJI0NPeb7wAkSZqMCUqSNEgmKEnSIJmgJEmDZIKSJA2S\nCUqSNEgmKKknST6R5F1TbHtNkm/0cM2zk/zZ+j5v3/r6eWhhM0FpMJKsTPKrJLckuTHJN5O8Pkmn\n39Mky5NUkg37jlVS/0xQGprnVdWmwMOB9wCHAEfPb0iLj0leQ2CC0iBV1U1VdSrwJ8D+SR4LkGTv\nJBcnuTnJj5McMXLY19vvG5P8MskeSR6R5GtJrk/y8ySfTrJkqusm+XB73puTXJjkD0a2HZHkxCSf\nbFt5lydZMbJ91yQXtdtOADaZoZpJ8o9JbkryvSTPHNlwQJIr2nP9IMnr1jpwvySXtHFeneQ5k5x8\nWZJLk/zPdn3HJF9vz/nVJB9Ncmy7baL1eWCSHwFfa8v3bet5Y9t9+OiR81eS3xlZv7tLM8nTk1yT\n5K1JViW5LskBI/tumeTUNv7zgEfM8LPSImSC0qBV1XnANcBEorgVeDWwBNgbeEOS57fbntp+L6mq\nB1fVt4AA7wa2AR4NbA8cMc0lzwd2AbYAjgM+m2Q00ewLfKa9/qnARwCSbAycDHyqPfazwB/PUL0n\nAz8AtgIOBz6fZIt22ypgH2Az4ADgg0me2F5rN+CTwF+1cTwVWDl64iTLgXOAj1TV+9ri44DzgC3b\nn8GrJonpaTQ/pz9K8kjgeODNwFLgdOALbV27eBiwObAtcCDw0SQPabd9FLgdWAb8afuR7sUEpYXg\nWpo/+lTV2VV1WVXdVVWX0vwBfdpUB1bVVVV1RlXdUVWrgQ/MsP+xVXV9Va2pqvcD9wceNbLLN6rq\n9Kr6DU0yekJbvjuwEfChqvp1VZ1Ek+yms2pk/xOAK2mSLlX1xaq6uhrnAF/hniR9IPDxtl53VdVP\nqup7I+d9DHA2cHhVHQWQZAfgScBhVXVnVX2DJsGu7YiqurWqfkXTev1ie51fA+8DHgD83gz1mvBr\n4B1t/U4Hfgk8KskGNMn7sPZa3wWO6XhOLSImKC0E2wK/AEjy5CRnJVmd5Cbg9TQtkEkl2TrJZ5L8\nJMnNwLEz7P/WtmvtpiQ30rQARvf/6cjybcAm7f2abYCf1L1nX/7hDPWabP9t2jj2SvLtJL9o43ju\nSBzbA1dPc95XAD8BThop2wb4RVXdNlL240mOHS3bZrQOVXVXu33baWt1j+uras3I+m3Ag2laYxuu\nda2ZflZahExQGrQkT6L5gzgxBPk4mv/z376qNgf+maYbD2Cyqfnf3ZY/vqo2A145sv/a1/oDmkEZ\nLwEeUlVLgJum2n8t1wHbJhndd4cZjpls/2uT3B/4HE2L5aFtHKePxPFjpr9ncwTwc+C4trUyEd8W\nSR44st/2kxw7+jO8lmawCtDcMGuP+UlbdBswer6HTRPTqNXAmrWuP9PPSouQCUqDlGSzJPvQ3O85\ntqouazdtStMSuL29F/PykcNWA3cBO42UbUrTtXRjkm1p7ttMZVOaP5yrgQ2THEZzD6iLb7XH/kWS\nDZO8ENhthmO2bvffKMmLae79nA5sTNO1uBpYk2Qv4Nkjxx0NHJDkmUnul2TbJL87sv3XwIuBBwGf\nSnK/qvohcAFwRJKNk+wBPG+G+E4E9m6vsxHwVuAO4Jvt9kuAlyfZoB2kMWXX6ai2e/TzbSwPTPIY\nYP8ux2pxMUFpaL6Q5BaaVsLf0NwzOmBk+58D72j3OYzmjygAbffV3wH/0Y462x04EngiTUvoizR/\nGKfyZeBLwP+j6XK6ncm7wX5LVd0JvBB4DXADzf2b6a4FcC6wM01r5++AF7X3v24B/qKt2w00Sfju\n+0XtwJEDgA+29TqHkZbOWvFsDXw8zbNkrwD2AK4H3gWcQJNwpqrTlTQtzn9sY3wezWMAd7a7HNyW\n3die++QZ6jvqTTTdfT8FPgH82zocq0UivrBQWpzaofDfq6rD5zsWaTK2oKRFIsmT0jwXdr+2S24/\n1q3VI80pnxaXFo+H0XQ7bknzbNkbquri+Q1JmppdfJKkQbKLT5I0SAuii2+rrbaq5cuXz3cYkqT1\n4MILL/x5VS2dab8FkaCWL1/OBRdcMN9hSJLWgySdZg6xi0+SNEgmKEnSIJmgJEmDZIKSJA2SCUqS\nNEgmKEnSIJmgJEmDZIKSJA2SCUqSNEgLYiYJadE7YvNZHnfT+o1DmkO2oCRJg2SCkiQNkglKkjRI\nJihJ0iCZoCRJg2SCkiQNkglKkjRIJihJ0iCZoCRJg2SCkiQNUu8JKskGSS5Oclq7vmOSc5N8P8kJ\nSTbuOwZJ0sIzFy2og4ErRtbfC3ywqnYGbgAOnIMYJEkLTK8JKsl2wN7Av7brAfYETmp3OQZ4fp8x\nSJIWpr5nM/8Q8NfApu36lsCNVbWmXb8G2HayA5McBBwEsMMOO/QcpjSmZjMLujOgayB6a0El2QdY\nVVUXjhZPsmtNdnxVHVVVK6pqxdKlS3uJUZI0XH22oJ4C7JvkucAmwGY0LaolSTZsW1HbAdf2GIMk\naYHqrQVVVW+rqu2qajnwUuBrVfUK4CzgRe1u+wOn9BWDJGnhmo/noA4B3pLkKpp7UkfPQwySpIGb\nk1e+V9XZwNnt8g+A3ebiupKkhcuZJCRJg2SCkiQNkglKkjRIJihJ0iCZoCRJg2SCkiQNkglKkjRI\nJihJ0iCZoCRJg2SCkiQNkglKkjRIJihJ0iCZoCRJg2SCkiQNkglKkjRIJihJ0iCZoCRJg2SCkiQN\n0oyvfE/yFOAI4OHt/gGqqnbqNzRJ0mI2Y4ICjgb+ErgQ+E2/4UiS1OiSoG6qqi/1HokkSSO6JKiz\nkvw98HngjonCqrqot6gkSYtelwT15PZ7xUhZAXuu/3AkSWrMmKCq6hlzEYgkSaOmTFBJXllVxyZ5\ny2Tbq+oD/YUlSVrspmtBPaj93nQuApEkadSUCaqqPtZ+Hzl34UiS1JhxJokkOyX5QpLVSVYlOSWJ\nD+lKknrVZaqj44ATgWXANsBngeP7DEqSpC4JKlX1qapa036OpRlmLklSb6YbxbdFu3hWkkOBz9Ak\npj8BvjgHsUmSFrHpRvFdSJOQ0q6/bmRbAe/sKyhJkqYbxbfjXAYiSdIo3wclSRokE5QkaZCmTVBp\nbD9XwUiSNGHaBFVVBZw8R7FIknS3Ll18307ypN4jkSRpRJf3QT0DeH2SlcCtNMPOq6oe32dgkqTF\nrUuC2qv3KCRJWsuMXXxV9UNge2DPdvm2Lscl2STJeUm+k+TyJEe25TsmOTfJ95OckGTj+1oJSdL4\n6ZJoDgcOAd7WFm0EHNvh3HfQJLUnALsAz0myO/Be4INVtTNwA3DgbAKXJI23LoMkXgDsS3P/iaq6\nlg4vMazGL9vVjdpPAXsCJ7XlxwDPX8eYJUmLQJcEdWc73LwAkjxohv3vlmSDJJcAq4AzgKuBG6tq\nTbvLNcC2Uxx7UJILklywevXqrpeUJI2JLgnqxCQfA5YkeS3wVeBfupy8qn5TVbsA2wG7AY+ebLcp\njj2qqlZU1YqlS5d2uZwkaYzMOIqvqt6X5FnAzcAjgcOq6ox1uUhV3ZjkbGB3mkS3YduK2g64dt3D\nliSNu65z8V0G/Dvw9XZ5RkmWJlnSLj8A+EPgCuAs4EXtbvsDp6xLwJKkxaHLKL4/A84DXkiTWL6d\n5E87nHsZzcsOLwXOB86oqtNoRgS+JclVwJbA0bMNXpI0vro8qPtXwK5VdT1Aki2BbwIfn+6gqroU\n2HWS8h/Q3I+SJGlKXbr4rgFuGVm/BfhxP+FIktTo0oL6CXBuklNoRtztB5yX5C0AVfWBHuOTJC1S\nXRLU1e1nwsSghhkf1pUkaba6DDM/ci4CkSRplK98lyQNkglKkjRIJihJ0iB1eVD3fyfZLMlGSc5M\n8vMkr5yL4CRJi1eXFtSzq+pmYB+aZ6IeSfPwriRJvemSoDZqv58LHF9Vv+gxHkmSgG7PQX0hyfeA\nXwF/nmQpcHu/YUmSFrsuLajDgT2AFVX1a+A2mjfsSpLUmy4J6ltVdUNV/Qagqm4FvtRvWJKkxW7K\nLr4kD6N5HfsDkuwKpN20GfDAOYhNkrSITXcP6o+A19C89XZ0QthbgLf3GJMkSVMnqKo6BjgmyR9X\n1efmMCZJkjqN4jstycuB5aP7V9U7+gpKkqQuCeoU4CbgQuCOfsORJKnRJUFtV1XP6T0SabE4YvP5\njkBaELoMM/9mksf1HokkSSO6tKB+H3hNkv+i6eILUFX1+F4jkyQtal0S1F69RyFJ0lpm7OKrqh8C\n2wN7tsu3dTlOkqT7osv7oA4HDgHe1hZtBBzbZ1CSJHVpCb2AZnLYWwGq6lpg0z6DkiSpS4K6s6oK\nKIAkD+o3JEmSuiWoE5N8DFiS5LXAV4F/6TcsSdJiN+Movqp6X5JnATcDjwIOq6ozeo9MkrSozZig\nkvwl8FmTkiRpLnXp4tsM+HKSf0/yxiQP7TsoSZK6PAd1ZFX9N+CNwDbAOUm+2ntkkqRFbV0euF0F\n/BS4Hti6n3AkSWp0eVD3DUnOBs4EtgJe6zx8kqS+dZmL7+HAm6vqkr6DkSRpQpd7UIcCD05yAECS\npUl27D0ySdKi5lx8kqRBci4+SdIgORefJGmQnItPkjRIvc3Fl2R74JPAw4C7gKOq6sNJtgBOAJYD\nK4GXVNUNs66BJGksdRlmTpuQ1nUuvjXAW6vqoiSbAhcmOQN4DXBmVb0nyaHAoTSDMCRJultvr26v\nquuq6qJ2+RbgCmBbYD/gmHa3Y4Dn9xWDJGnh6i1BjUqyHNgVOBd4aFVdB00Sw2mTJEmTmLKLL8mZ\nVfXMJO+tqll3wSV5MPA5mtkobk7S9biDgIMAdthhh9leXtK6OmLzWR530/qNQ4vedPegliV5GrBv\nks8A98osE91300myEU1y+nRVfb4t/lmSZVV1XZJlNJPQ/paqOgo4CmDFihU1c1UkSeNkugR1GM0A\nhu2AD6y1rYA9pztxmqbS0cAVVTV6/KnA/sB72u9T1jFmSdIiMGWCqqqTgJOS/K+qeucszv0U4FXA\nZUkmJpp9O01iOjHJgcCPgBfP4tySpDHX5TmodybZF3hqW3R2VZ3W4bhvsFa34Ihndg9RkrQYdZks\n9t3AwcB/tp+D2zJJknrT5UHdvYFdquougCTHABdzz+zmkiStd12fg1oysjzLMaiSJHXXpQX1buDi\nJGfR3FN6KraeJEk96zJI4vgkZwNPoklQh1TVT/sOTJK0uHWdLPY6mueXJEmaE3MyF58kSevKBCVJ\nGqRpE1SS+yX57lwFI0nShGnvQVXVXUm+k2SHqvrRXAWl/i0/9IvrfMzK9+zdQySSNLkugySWAZcn\nOQ+4daKwqvbtLSpJ0qLXJUEd2XsUkiStpctzUOckeTiwc1V9NckDgQ36D02StJh1mSz2tcBJwMfa\nom2Bk/sMSpKkLsPM30jzbqebAarq+8DWfQYlSVKXBHVHVd05sZJkQ5o36kqS1JsuCeqcJG8HHpDk\nWcBngS/0G5YkabHrkqAOBVYDlwGvA04H/rbPoCRJ6jKK7672JYXn0nTtXVlVdvFJkno1Y4JKsjfw\nz8DVNK/b2DHJ66rqS30HJw3ZbGbjAFi5yXoOZIFzVhNNpcuDuu8HnlFVVwEkeQTwRcAEJUnqTZd7\nUKsmklPrB8CqnuKRJAmYpgWV5IXt4uVJTgdOpLkH9WLg/DmITZK0iE3Xxfe8keWfAU9rl1cDD+kt\nIkmSmCZBVdUBcxmIJEmjuozi2xH4H8Dy0f193YYkqU9dRvGdDBxNM3vEXf2GI0lSo0uCur2q/qH3\nSCRJGtElQX04yeHAV4A7Jgqr6qLeopIkLXpdEtTjgFcBe3JPF1+165rErGcY8On4+8yfvdTd0Gfx\n6JKgXgDsNPrKDUmS+tZlJonvAEv6DkSSpFFdWlAPBb6X5HzufQ/KYeaSpN50SVCH9x6Fxtps7wvN\npbmMcfntx63zMSs3eXkPkWhoFsK/lbnU5X1Q58xFIJIkjeoyk8QtNKP2ADYGNgJurarN+gxMkrS4\ndWlBbTq6nuT5wG69RSRJEt1G8d1LVZ2Mz0BJknrWpYvvhSOr9wNWcE+X31jzhqU0Xob+YKrurcso\nvtH3Qq0BVgL79RKNJEmtLvegZvVeqCQfB/aheWX8Y9uyLYATaF7dsRJ4SVXdMJvzS5LG23SvfD9s\nmuOqqt45w7k/AXwE+ORI2aHAmVX1niSHtuuHdIxVkrSITDdI4tZJPgAH0iGpVNXXgV+sVbwfcEy7\nfAzw/HUJVpK0eEz3yvf3Tywn2RQ4GDgA+Azw/qmOm8FDq+q69vzXJdl6qh2THAQcBLDDDjvM8nKS\n5swRm8/ywHWfWUOLw7TDzJNskeRdwKU0yeyJVXVIVa3qO7CqOqqqVlTViqVLl/Z9OUnSwEyZoJL8\nPXA+cAvwuKo6Yj0MaPhZkmXt+ZcBvSc6SdLCNF0L6q3ANsDfAtcmubn93JLk5lle71Rg/3Z5f+CU\nWZ5HkjTmprsHtc6zTIxKcjzwdGCrJNfQzIr+HuDEJAcCPwJefF+uIUkaX10e1J2VqnrZFJue2dc1\n1S9n1tBiNNvfe2eguO/uUytJkqS+mKAkSYNkgpIkDZIJSpI0SL0NkpC0/iy/fXazLazc5OXrOZJh\ncMDO4mALSpI0SCYoSdIgmaAkSYNkgpIkDdKiGSSxEG6qLoQYpanMdiDHuPLf831nC0qSNEgmKEnS\nIJmgJEmDZIKSJA3SohkkIS1Gsxm4MK6zT2jhsQUlSRokE5QkaZBMUJKkQTJBSZIGyUESku7FGSE0\nFLagJEmDZIKSJA2SCUqSNEgmKEnSIJmgJEmDZIKSJA2SCUqSNEgmKEnSIJmgJEmDZIKSJA2SCUqS\nNEgmKEnSIJmgJEmDZIKSJA2SCUqSNEgmKEnSIJmgJEmDZIKSJA2SCUqSNEgmKEnSIM1LgkrynCRX\nJrkqyaHzEYMkadjmPEEl2QD4KLAX8BjgZUkeM9dxSJKGbT5aULsBV1XVD6rqTuAzwH7zEIckacA2\nnIdrbgv8eGT9GuDJa++U5CDgoHb1l0muvI/X3Qr4+X08x5CNe/1g/Os47vUD67jg5b3rpX4P77LT\nfCSoTFJWv1VQdRRw1Hq7aHJBVa1YX+cbmnGvH4x/Hce9fmAdx8Fc1m8+uviuAbYfWd8OuHYe4pAk\nDdh8JKjzgZ2T7JhkY+ClwKnzEIckacDmvIuvqtYkeRPwZWAD4ONVdfkcXHq9dRcO1LjXD8a/juNe\nP7CO42DO6peq37r9I0nSvHMmCUnSIJmgJEmDNBYJKsnHk6xK8t2Rsi2SnJHk++33Q9ryJPmHdpql\nS5M8cf4i7y7J9knOSnJFksuTHNyWj0U9k2yS5Lwk32nrd2RbvmOSc9v6ndAOrCHJ/dv1q9rty+cz\n/nWRZIMkFyc5rV0fmzomWZnksiSXJLmgLRuL39EJSZYkOSnJ99p/j3uMUx2TPKr97zfxuTnJm+ej\njmORoIBPAM9Zq+xQ4Myq2hk4s12HZoqlndvPQcA/zVGM99Ua4K1V9Whgd+CNaaaIGpd63gHsWVVP\nAHYBnpNkd+C9wAfb+t0AHNjufyBwQ1X9DvDBdr+F4mDgipH1cavjM6pql5FnZcbld3TCh4H/W1W/\nCzyB5r/l2NSxqq5s//vtAvx34Dbg/zAfdayqsfgAy4HvjqxfCSxrl5cBV7bLHwNeNtl+C+kDnAI8\naxzrCTwQuIhmhpGfAxu25XsAX26Xvwzs0S5v2O6X+Y69Q922o/nHvSdwGs2D62NTR2AlsNVaZWPz\nOwpsBvzX2v8dxqmOa9Xr2cB/zFcdx6UFNZmHVtV1AO331m35ZFMtbTvHsd0nbVfPrsC5jFE9266v\nS4BVwBnA1cCNVbWm3WW0DnfXr91+E7Dl3EY8Kx8C/hq4q13fkvGqYwFfSXJhmunKYIx+R4GdgNXA\nv7XdtP+a5EGMVx1HvRQ4vl2e8zqOc4KaSqeploYqyYOBzwFvrqqbp9t1krJB17OqflNNt8J2NJMK\nP3qy3drvBVe/JPsAq6rqwtHiSXZdsHUEnlJVT6Tp9nljkqdOs+9CrN+GwBOBf6qqXYFbuaerazIL\nsY4AtPdC9wU+O9Ouk5StlzqOc4L6WZJlAO33qrZ8wU61lGQjmuT06ar6fFs8dvWsqhuBs2nutS1J\nMvFA+Wgd7q5fu31z4BdzG+k6ewqwb5KVNLP470nTohqbOlbVte33Kpr7FrsxXr+j1wDXVNW57fpJ\nNAlrnOo4YS/goqr6Wbs+53Uc5wR1KrB/u7w/zT2bifJXtyNPdgdummi2DlmSAEcDV1TVB0Y2jUU9\nkyxNsqRdfgDwhzQ3n88CXtTutnb9Jur9IuBr1XaAD1VVva2qtquq5TRdJ1+rqlcwJnVM8qAkm04s\n09y/+C5j8jsKUFU/BX6c5FFt0TOB/2SM6jjiZdzTvQfzUcf5vgm3nm7kHQ9cB/yaJpsfSNNXfybw\n/fZ7i3bf0Lww8WrgMmDFfMffsY6/T9NsvhS4pP08d1zqCTweuLit33eBw9rynYDzgKtouhru35Zv\n0q5f1W7fab7rsI71fTpw2jjVsa3Hd9rP5cDftOVj8Ts6Us9dgAva39WTgYeMYR0fCFwPbD5SNud1\ndKojSdIgjXMXnyRpATNBSZIGyQQlSRokE5QkaZBMUJKkQTJBSZIGyQQlSRokE5Q0z5I8qX2Pzibt\nbAyXJ3nsfMclzTcf1JUGIMm7aGaOeADNXG/vnueQpHlngpIGoJ05+nzgduD3quo38xySNO/s4pOG\nYQvgwcCmNC0padGzBSUNQJJTaV7BsSPN20jfNM8hSfNuw5l3kdSnJK8G1lTVcUk2AL6ZZM+q+tp8\nxybNJ1tQkqRB8h6UJGmQTFCSpEEyQUmSBskEJUkaJBOUJGmQTFCSpEEyQUmSBun/A07j8qmrLlT7\nAAAAAElFTkSuQmCC\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": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNXd+PHPNyQmLAlrAgSIAYuAbAECBbUqUq3W3Wqt\ntYo+Vuzz+Fj7s1WpWpdKq89Tq631acXWVlqta+tSiysFxZUGQVRAkX1PWAJB1pDv749zJw6Q5WaS\nO/fO5Pt+veZ1Z+7c5XuGmXw55557jqgqxhhjTNRkhB2AMcYYUxdLUMYYYyLJEpQxxphIsgRljDEm\nkixBGWOMiSRLUMYYYyLJEpQxARGRh0VkSj3vXSoibwZwzlki8t2WPm7Qgvo8TGqzBGUiQ0RWiMgu\nEakSkUoReVtEvicivr6nIlIsIioimUHHaowJniUoEzVnqGoucDhwF3AD8FC4IbU+luRNFFiCMpGk\nqttU9XngAmCiiAwBEJHTRGSeiGwXkdUiclvcbm94y0oR2SEi40TkCBH5l4hsFpFNIvKoiHSq77wi\n8mvvuNtFZK6IfCXuvdtE5EkR+bNXy/tYRErj3h8hIu977z0B5DRSTBGR34jINhFZLCIT4t64TEQW\necdaJiJXHrTjWSIy34tzqYicUsfBe4rIAhH5kfe6r4i84R3zNRH5PxF5xHsvVvu8XERWAf/y1p/p\nlbPSaz4cFHd8FZEvxb2ubdIUkRNEZI2I/FBEykVkvYhcFrdtVxF53ot/DnBEI5+VaYUsQZlIU9U5\nwBoglig+By4BOgGnAf8pImd77x3nLTupagdVfQcQ4E6gEBgE9AFua+CU/wZKgC7AX4GnRCQ+0ZwJ\nPO6d/3ngfgAROQx4FviLt+9TwDcaKd6XgWVAN+BW4O8i0sV7rxw4HcgDLgPuFZGR3rnGAH8GrvPi\nOA5YEX9gESkGXgfuV9W7vdV/BeYAXb3P4OI6Yjoe9zl9TUSOBB4DfgDkA9OBf3hl9aMH0BHoBVwO\n/J+IdPbe+z9gN9AT+A/vYcwBLEGZVLAO90cfVZ2lqh+qao2qLsD9AT2+vh1V9TNVfVVV96hqBXBP\nI9s/oqqbVbVaVX8JZAMD4jZ5U1Wnq+p+XDIa7q0fC2QBv1LVfar6NC7ZNaQ8bvsngE9wSRdV/aeq\nLlXndeAVvkjSlwN/9MpVo6prVXVx3HGPAmYBt6rqgwAiUgSMBm5R1b2q+iYuwR7sNlX9XFV34Wqv\n//TOsw+4G2gLHN1IuWL2AT/1yjcd2AEMEJE2uOR9i3euj4BpPo9pWhFLUCYV9AK2AIjIl0VkpohU\niMg24Hu4GkidRKRARB4XkbUish14pJHtf+g1rW0TkUpcDSB++w1xz3cCOd71mkJgrR44+vLKRspV\n1/aFXhynisi7IrLFi+PrcXH0AZY2cNyLgLXA03HrCoEtqrozbt3qOvaNX1cYXwZVrfHe79Vgqb6w\nWVWr417vBDrgamOZB52rsc/KtEKWoEykicho3B/EWBfkv+L+599HVTsCD+Ca8QDqGpr/Tm/9MFXN\nA74Tt/3B5/oKrlPGN4HOqtoJ2Fbf9gdZD/QSkfhtixrZp67t14lINvA3XI2luxfH9Lg4VtPwNZvb\ngE3AX73aSiy+LiLSLm67PnXsG/8ZrsN1VgHcBTNvn7Xeqp1A/PF6NBBTvAqg+qDzN/ZZmVbIEpSJ\nJBHJE5HTcdd7HlHVD723cnE1gd3etZhvx+1WAdQA/eLW5eKalipFpBfuuk19cnF/OCuATBG5BXcN\nyI93vH2/LyKZInIuMKaRfQq87bNE5HzctZ/pwGG4psUKoFpETgVOjtvvIeAyEZkgIhki0ktEBsa9\nvw84H2gP/EVEMlR1JVAG3CYih4nIOOCMRuJ7EjjNO08W8ENgD/C29/584Nsi0sbrpFFv02k8r3n0\n714s7UTkKGCin31N62IJykTNP0SkCldLuAl3zeiyuPf/C/ipt80tuD+iAHjNVz8D3vJ6nY0FbgdG\n4mpC/8T9YazPy8CLwKe4Jqfd1N0MdghV3QucC1wKbMVdv2noXADvAf1xtZ2fAed517+qgO97ZduK\nS8K114u8jiOXAfd65XqduJrOQfEUAH8Udy/ZRcA4YDMwBXgCl3DqK9MnuBrnb7wYz8DdBrDX2+Qa\nb12ld+xnGylvvP/GNfdtAB4G/tSEfU0rITZhoTGtk9cVfrGq3hp2LMbUxWpQxrQSIjJa3H1hGV6T\n3Fk0rdZjTFLZ3eLGtB49cM2OXXH3lv2nqs4LNyRj6mdNfMYYYyLJmviMMcZEUko08XXr1k2Li4vD\nDsMYY0wLmDt37iZVzW9su5RIUMXFxZSVlYUdhjHGmBYgIr5GDrEmPmOMMZFkCcoYY0wkWYIyxhgT\nSSlxDcoYk/r27dvHmjVr2L17d9ihmCTJycmhd+/eZGVlJbS/JShjTFKsWbOG3NxciouLOXAQd5OO\nVJXNmzezZs0a+vbtm9AxrInPGJMUu3fvpmvXrpacWgkRoWvXrs2qMVuCMsYkjSWn1qW5/97WxGdM\nFCx/AypXQbcjoU9j00gZ0zpYDcqYMO2vhunXw7Qz4Lmr4KGT4F8/g5qasCNLSx06dGh0m+9+97ss\nXLgQgJ///OcHvHf00UcnfI42bdpQUlJS+1ixYgVlZWV8//vfB2DWrFm8/fbbde778MMPk5+fX7vv\nJZdc0mgcTXHwuR944AH+/Oc/t+g5EmE1KGPCNPuXMGcqjL0KvjwJXv8FvPG/0KEAxlwRdnSt0h/+\n8Ifa5z//+c+58cYba1/Xl0D8aNu2LfPnzz9gXXFxMaWlpYBLEh06dKg3CV5wwQXcf//9CZ+/IQef\n+3vf+14g52kqq0EZE5aKT2H23TDkPDjl59C5GM66H444EV67HbavCzvCtDVr1ixOOOEEzjvvPAYO\nHMhFF11EbGaHE044gbKyMiZPnsyuXbsoKSnhoosuAr6oHe3YsYMJEyYwcuRIhg4dynPPPZdwHKef\nfjorVqzggQce4N5776WkpITZs2f72j8WK8CmTZuIjVn68MMPc+6553LKKafQv39/rr/++tp9Xnrp\nJUaOHMnw4cOZMGFCnee+7bbbuPvuuwGYP38+Y8eOZdiwYZxzzjls3bq19tw33HADY8aM4cgjj/Qd\nc1NYDcqYsLxyM2S1hVPu/GKdCJx2D/x2HMy4A875XXjxBenFybDhw5Y9Zo+hcOpdvjefN28eH3/8\nMYWFhRxzzDG89dZbHHvssbXv33XXXdx///2H1HrA3d/zzDPPkJeXx6ZNmxg7dixnnnlmg50CYskO\noG/fvjzzzDO17xUXF/O9732PDh068KMf/ajO/Z944gnefPNNAK655houu+yyBss3f/585s2bR3Z2\nNgMGDODqq68mJyeHK664gjfeeIO+ffuyZcsWunTpcsi5Z8yYUXucSy65hN/85jccf/zx3HLLLdx+\n++386le/AqC6upo5c+Ywffp0br/9dl577bUGY2oqS1DGhGHzUljyMpxwo2vOi9elL4y8BOb+CU66\n/dD3TYsYM2YMvXv3Bqi9JhSfoBqiqtx444288cYbZGRksHbtWjZu3EiPHj3q3aeuJr6maGoT34QJ\nE+jYsSMARx11FCtXrmTr1q0cd9xxtfcldenSpcFjbNu2jcrKSo4//ngAJk6cyPnnn1/7/rnnngvA\nqFGjWLFiRVOK44slKGPCMOdByMiCUZfW/f6YSe7a1NyH4fjr694mlTWhphOU7Ozs2udt2rShurra\n976PPvooFRUVzJ07l6ysLIqLi0MZISMzM5Mar0PNweevq3yq2qJd/WPnaOrn55ddgzIm2fbuhHmP\nwpBzIbd73dt0+xJ86atQ9keo2Z/c+EytrKws9u3bd8j6bdu2UVBQQFZWFjNnzmTlSl+zRzQoNzeX\nqqqqJu1TXFzM3LlzAXj66acb3X7cuHG8/vrrLF++HIAtW7Y0eO6OHTvSuXPn2utLf/nLX2prU8lg\nCcqYZFvyCuytghHfaXi7Ed+BqvWw6p3kxGUOMWnSJIYNG1bbSSLmoosuoqysjNLSUh599FEGDhzY\n7HOdccYZPPPMM03qJPGjH/2I3/3udxx99NFs2rSp0e3z8/N58MEHOffccxk+fDgXXHBBo+eeNm0a\n1113HcOGDWP+/PnccsstTS9cgiTWcyXKSktL1SYsNGnjqcvcjbk/+hQy2tS/3Z4d8IsvwciL4eu/\nSF58AVm0aBGDBg0KOwyTZHX9u4vIXFUtbWxfq0EZk0z7dsGnL8OgMxpOTgDZHaD/V2Hh83bjrmmV\nLEEZk0yfzYB9n8NRZ/nb/qizYccGWP1esHEZE0GBJSgRGSAi8+Me20XkByLSRUReFZEl3rJzUDEY\nEzlLXoHsPCj2152Z/ieBtIHPXg02LmMiKLAEpaqfqGqJqpYAo4CdwDPAZGCGqvYHZnivjUl/qrBs\nJhR/Bdr4nMAtpyP0LoWlM4ONzZgISlYT3wRgqaquBM4CpnnrpwFnJykGY8K1ZZkbsfyI8U3br994\nWDcPdm4JJi5jIipZCepbwGPe8+6quh7AW9Z5m7yITBKRMhEpq6ioSFKYxgRomVcL6tfEBHXEiYC6\nnn+tzAVT3+GCqdbNvrUKPEGJyGHAmcBTTdlPVR9U1VJVLc3Pzw8mOGOSaelM6NgHuh7RtP16jXLX\nrZZZM19zxaa8GD58OCNHjqwdnXzdunWcd955gZ+/uLi4zvuViouLGTp0KEOHDuWoo47i5ptvZs+e\nPQ0eq7Kykt/+9rdBhRoJyahBnQq8r6obvdcbRaQngLcsT0IMxoRLFVa+DX2PcwPCNkWbTDj8aFjx\nVjCxRdje6hoWrttOeVXLDCMUGw/vgw8+4M477+THP/4xAIWFhb5GYgjSzJkz+fDDD5kzZw7Lli1j\n0qRJDW5vCaplXMgXzXsAzwMTvecTgcTGqTcmlWxaAru2QNHYxPYvGgubl8DnjY8WkE7WVu6iak81\n9722pMWPvX37djp3dp2IV6xYwZAhQ4CGp6ro0KEDN910E8OHD2fs2LFs3Oj+311RUcE3vvENRo8e\nzejRo3nrLfefic2bN3PyySczYsQIrrzySvwMjNChQwceeOABnn32WbZs2VLv1B6TJ09m6dKllJSU\ncN1117XYFCCRoqqBPYB2wGagY9y6rrjee0u8ZZfGjjNq1Cg1JqWVPax6a55qxaeJ7b/yHbf/wn+0\nbFxJtHDhQt/bHnnTdD38hhcOeRx50/RmxZCRkaHDhw/XAQMGaF5enpaVlamq6vLly3Xw4MGqqvqn\nP/1J+/btq5WVlbpr1y4tKirSVatWqaoqoM8//7yqql533XV6xx13qKrqhRdeqLNnz1ZV1ZUrV+rA\ngQNVVfXqq6/W22+/XVVVX3jhBQW0oqLikLgOP/zwQ9YPHz5c3333Xd23b59u27ZNVVUrKir0iCOO\n0JqamgNiVtV6twtbXf/uQJn6yCGBjmauqju9hBS/bjOuV58xrceqd6FdV+j6pcT2LxwBbbJh9bsw\n6PSWjS2CZl8/ninTF/HCB+uoUcjJyuBrg3tw02nNGyopfsqLd955h0suuYSPPvrokO3qmqqiT58+\nHHbYYZx+uvv8R40axauvuvvTXnvttdpp4sHVzqqqqnjjjTf4+9//DsBpp51WW2PzQ73altYztUdd\n2zd1CpCos+k2jEmG1e9C0bimX3+KycyGXiNdomsFCvJyyM3OpEbdR7anuobc7EwKcnNa7Bzjxo1j\n06ZN1NVLuL6pOLKysmqnq4hfX1NTwzvvvEPbtm0POVYi01tUVVWxYsUKjjzySN9Te0RlCpCWZEMd\nGRO0HeXuHqg+X27ecYrGuvuh9u1qmbgibtOOPRTkZjO4Zx4XfflwKnY03KutqRYvXsz+/fvp2rVr\n4xs34uSTTz5gMsFYLe24447j0UcfBeDFF1+snS69ITt27OC//uu/OPvss+ncuXO9U3scPEVGEFOA\nhM1qUMYEbe37btl7dPOO03s01FTD+gVQ1MxklwKmXlxaew/UlLOHtMgx46ddV1WmTZtGmzaNDNrr\nw3333cdVV13FsGHDqK6u5rjjjuOBBx7g1ltv5cILL2TkyJEcf/zxFBUV1XuM8ePHo6rU1NRwzjnn\n8JOf/ARwU3ucccYZlJaWUlJSUju1R9euXTnmmGMYMmQIp556KjfccEOd26Uym27DmKDNuss9frzG\njVCeqKoN8MsBcMpdMPY/Wy6+JElkuo1YgnriynFBhGSSoDnTbVgNypigrZsP3Y5sXnICyO0BuYWw\ndm7LxJUCLDG1bnYNypigrZsHhSUtc6xeI79oMjQmzVmCMiZIVRvcfE49WzBBbVkKuxq/2B5FqXBJ\nwbSc5v57W4IyJkjrXG+uFqtBFY70jjuvZY6XRDk5OWzevNmSVCuhqmzevJmcnMRvDbBrUMYEaf18\nQKDHsJY5XuEIt1w3zxvlPHX07t2bNWvW1HnfkUlPOTk59O7dO+H9LUEZE6SW6iAR07YTdCqCDYeO\nfhB1WVlZ9O3bN+wwTAqxJj5jgrR+fss178X0GAYbUy9BGdNUlqCMCUrVBqha33IdJGK6D4HNn8He\nnS17XGMixhKUMUFp6Q4SMT2GgtZA+aKWPa4xEWMJypigtHQHiZge3rA/Gxa07HGNiRhLUMYEZd18\n6Na/5TpIxHQ63E0Bb9ehTJqzBGVMUDYsgJ7DW/64Iu461IYPW/7YxkSIJShjgrCrEravhe6Dgzl+\nj6Gw8WOoqQnm+MZEgCUoY4JQsdgt85s3A2y9egyBvTtg6/Jgjm9MBFiCMiYIsR52BQHNydNjqFta\nM59JY40mKBE5RkReFZFPRWSZiCwXkWV+Di4inUTkaRFZLCKLRGSciHTxjrfEW3ZufjGMiZiKxZDV\nHjrWP0Fds+QPAmljHSVMWvNTg3oIuAc4FhgNlHpLP34NvKSqA4HhwCJgMjBDVfsDM7zXxqSX8oWQ\nPwAyAmqkyMpxQyhZDcqkMT+/nm2q+qKqlqvq5tijsZ1EJA84DpfgUNW9qloJnAVM8zabBpydYOzG\nRFf5Yig4Kthz9BiSkmPyGeOXnwQ1U0R+4TXPjYw9fOzXD6gA/iQi80TkDyLSHuiuqusBvGVB4uEb\nE0Gfb4bPy4O7/hTTYyhsXwM7twR7HmNC4mc08y97y/j54xVobKz/TGAkcLWqviciv6YJzXkiMgmY\nBFBUFFA7vjFBqPA6SATVgy8m1oW9fCEUHxvsuYwJQaMJSlXHJ3jsNcAaVX3Pe/00LkFtFJGeqrpe\nRHoC5fWc90HgQYDS0lKb4cykjtoefAEnqFgCLF9kCcqkpXoTlIh8R1UfEZFr63pfVe9p6MCqukFE\nVovIAFX9BJgALPQeE4G7vOVzCUdvTBSVL3JDEeUVBnuevEJ3ntg9V8akmYZqUO29ZW4zjn818KiI\nHAYsAy7DXfd6UkQuB1YB5zfj+MZET8ViV3sSCfY8IpA/0HXIMCYN1ZugVHWqt7w90YOr6nwOvHYV\nMyHRYxoTaarumtCgM5NzvoKBsPifyTmXMUnm50bdfiLyDxGpEJFyEXlORPolIzhjUs6Octi1Nfjr\nTzH5g2DnZthRkZzzGZNEfrqZ/xV4EugJFAJPAY8FGZQxKasiSR0kYmJd2Sts8kKTfvwkKFHVv6hq\ntfd4BNfN3BhzsPIkdTGPqe3JZ9ehTPppqBdfF+/pTBGZDDyOS0wXANbobUxdyhdB2y7QIUn3n+f2\ngJyO7rqXMWmmoV58c3EJKdYV6cq49xS4I6igjElZ5YuS04MvRsTVoqyruUlDDfXi65vMQIxJeaou\nUQxN8p0TBQPh42fd+ZOVGI1JApsPypiWsn0d7NmevA4SMfmDYHcl7NiY3PMaEzBLUMa0lGQNcXSw\nWE++cuvJZ9JLgwlKnD7JCsaYlJasQWIPFjufXYcyaabBBKWqCjybpFiMSW3li6B9AbTvmtzzdihw\nPQetBmXSjJ8mvndFxO8Musa0XuWLgp8Dqi4irlnRalAmzfhJUONxSWqpiCwQkQ9FZEHQgRmTUmpq\noOKT4GfRrU9s0Fi1e+hN+vAzYeGpgUdhTKrbtgr2fe4SRRgKBsGebVC1PvhpPoxJkkZrUKq6EugD\nnOg93+lnP2NaldhQQ2HWoMCuQ5m04mc081uBG4Afe6uygEeCDMqYlBMbaih/QDjnL7CefCb9+KkJ\nnQOcCXwOoKrraN4khsakn4rFkFsIbTuFc/723aBdNxuTz6QVPwlqr9fdXAFEpH0j2xvT+sTG4AtT\nwSAb1dykFT8J6kkRmQp0EpErgNeA3wcbljEppGY/bPo0/ASVP9D1JLSefCZNNNqLT1XvFpGTgO3A\nkcAtqvpq4JEZkyq2roDq3eEnqIKBsLcKtq2BTjYAjEl9frqZA3wItMU1830YXDjGpKDaDhJh16Di\nOkpYgjJpwE8vvu8Cc4BzgfNwN+3+h5+Di8gK78be+SJS5q3rIiKvisgSb9m5OQUwJnSx6z5h9eCL\nidXgrKu5SRN+rkFdB4xQ1UtVdSIwCtft3K/xqlqiqqXe68nADFXtD8zwXhuTuioWQaciyO4Qbhzt\nurixAK2ruUkTfhLUGqAq7nUVsLoZ5zwLmOY9nwac3YxjGRO+8kXhN+/FFAy0GpRJG34S1FrgPRG5\nzbtp913gMxG5VkSubWRfBV4RkbkiMslb111V1wN4y4K6dhSRSSJSJiJlFRUV/kpjTLLt3webloTf\nQSImf5DryVdTE3YkxjSbn04SS71HzHPe0s/Nuseo6joRKQBeFRHfbQ+q+iDwIEBpaan1mzXRtHkp\n1OwLLUFdMPUdAJ64cpxbUTDIjQm4bTV0PjyUmIxpKX66md+e6MG9USdQ1XIReQYYA2wUkZ6qul5E\negLliR7fmNDVTlIY0iCxB4sf8sgSlElxgQ36KiLtRSQ39hw4GfgIeB6Y6G02kS9qZMaknvLFgITf\ngy/GBo01acTvfVCJ6A48IyKx8/xVVV8SkX/jRqe4HFgFnB9gDMYEq3whdOkLWW3DjsRp2wlye1pP\nPpMWAktQqroMGF7H+s3AhKDOa0xSVSwOb4qN+uRbTz6THvzcqPu/IpInIlkiMkNENonId5IRnDGR\nVr3HdZKIyvWnmIJBbmxA68lnUpyfa1Anq+p24HTcPVFH4m7eNaZ127QEdH90upjH5A+EfTuhcmXY\nkRjTLH4SVJa3/DrwmKpuCTAeY1JH7DpP1BKUTV5o0oSfBPUP7/6lUmCGiOQDu4MNy5gUUL4QpA10\n/VLYkRwo1qPQrkOZFOcnQd0KjANKVXUfsBM3w64xrVv5Yuh6BGRmhx3JgXI6Ql6vBhPUBVPfqb3J\n15io8pOg3lHVraq6H0BVPwdeDDYsY1JA+cLoNe/F5A/84iZiY1JUvd3MRaQH0AtoKyIjAPHeygPa\nJSE2Y6Jr7043UeGwC8KOpG4Fg+Dfb7nZfjPahB2NMQlp6D6orwGXAr2Be+LWVwE3BhiTMdFXsRjQ\naNegqne7JNr1iLCjMSYh9SYoVZ0GTBORb6jq35IYkzHRF7u+031wuHHUJ37yQktQJkX5GUniBRH5\nNlAcv72q/jSooIyJvPKFkJkDXfqFHUndYj35KhbBoNPDjcWYBPlJUM8B24C5wJ5gwzEmRWz82CWB\nqF7fyc6FjkVfTEdvTAryk6B6q+opgUdiTCopXwRHjA87CvZW1/BZ+Q7Kq3ZTkJtz4JsFA+1mXZPS\n/HQzf1tEhgYeiTGpYucW2LEhEoPErq3cRdWeau57bcmhb+YPdGPy7a9OfmDGtAA/NahjgUtFZDmu\niU8AVdVhgUZmTFRt/Ngtu4eXoAbc/CJ7qr8YDPaR91bxyHuryM7M4JMpp7qVBYNg/17Yuhy69Q8p\nUmMS5ydBnRp4FMakkvKFblkQXg++2dePZ8r0RbzwwTpqFHKyMvja4B7cdFpct/fayQsXWoIyKanR\nJj5VXQn0AU70nu/0s58xaat8IeR0gtweoYVQkJdDbnYmNQoisKe6htzszAOvQ+UPBMmAjQtDi9OY\n5mi0BiUit+IGih0A/Ak3uvkjwDHBhmZMRG1c6O5/Eml82wBt2rGHgtxsCnKzKSnqTEXVQWM4H9YO\nuhwBGz8KJ0BjmslPTegc3OCwnwOo6jogN8igjIksVdeDLwIdJKZeXErfbu1pn53JlLOHMPXi0kM3\n6jHEEpRJWX4S1F5VVUABRKR9sCEZE2HbVsPequgOcXSw7oPdcEe7t4cdiTFN5idBPSkiU4FOInIF\n8Brwe78nEJE2IjJPRF7wXvcVkfdEZImIPCEihyUWujEhiF3PieoQRwfr7t0hUn7gdai91TUsXLed\n8oObBY2JED+dJO4Gngb+hrsOdYuq/qYJ57gGiB/3/3+Ae1W1P7AVuLwJxzImXOVeF/NUqkHBIc18\nDd4/ZUxE+Okk8f+Ap1T11aYeXER6A6cBPwOuFREBTgS+7W0yDbgN+F1Tj21MKMoXQcc+blLAVNCx\nt4t1g0tQvu6fMiYi/DTx5QEvi8hsEblKRLo34fi/Aq4HYr+IrkClqsZubV+Dm3PqECIySUTKRKSs\noqKiCac0JkAbPopEBwnfRFwzn3dz8ezrx3NmSSEZXgfEnKwMziopZPYN4Q/bZMzB/DTx3a6qg4Gr\ngELgdRF5rbH9ROR0oFxV58avrusU9Zz3QVUtVdXS/Pz8xk5nTPD27XJDB/VMsUFUug92Caqmxt/9\nU8ZEhJ+RJGLKgQ3AZqDAx/bHAGeKyNeBHFxN7Fe4zhaZXi2qN7CuaSEbE5KNC0H3Q48US1A9hsC+\nz6FyBXTp1/j9U8ZERKM1KBH5TxGZBcwAugFX+BmHT1V/rKq9VbUY+BbwL1W9CJgJnOdtNhE3nYcx\n0bfhA7dMxRoU1F6H8nX/lDER4Oca1OHAD1R1sKreqqrNHTflBlyHic9w16QeaubxjEmO9Qtch4NO\nh4cdSdPkD/KGPLIbdk1q8XMNajLQQUQuAxCRfBHp25STqOosVT3de75MVceo6pdU9XxVtUkQTWrY\nsMA174U8xFGT1Q559HHYkRjTJH6a+G7F1Xp+7K2KjcVnTOuxv9r9gU+1608xPYbAhg/DjsKYJrGx\n+IzxY/OdLOpaAAAYk0lEQVQSqN6detefYroPhsqVsHtb2JEY45uNxWeMH+sXuGWq1qB6jnDLWDmM\nSQGBj8VnTFrYsAAyc6DbkWFHkpjCErdcNy/cOIxpgkbvg1LVu0XkJGA7X4zF1+Rhj4xJaes/cCNI\ntGnKrYPBe+LKcf42bN8NOhZZgjIpxdevzUtIlpRM66TqalCDzwk7kubpNcISlEkpNnW7MY2pXOU6\nF6Tq9aeYwhGwdTns3BJ2JMb4YgnKmMZs8DoW9BwebhzNVRjrKDE/3DiM8aneBCUiM7zl/yQvHGMi\naP0HbiSGVBrFvC6xBGvNfCZFNHQNqqeIHI8b8PVxDhqJXFXfDzQyY6JiTZm7j+iwdmFH0jxtO0OX\nfl6COjrsaIxpVEMJ6hZgMm7E8XsOek9xEw8ak95qamDtXBh6XuPbpoLCEbB6DqR4rjWtQ70JSlWf\nBp4WkZ+o6h1JjMmY6Nj0KezZDr1Hhx1JyygcAR/9jbzsSra36RR2NMY0yM99UHeIyJnAcd6qWar6\nQrBhGRMRa/7tlumUoIDfn5QJ/X3eQ2VMSPwMFnsncA2w0Htc460zJv2tLXNTbHQ5IuxIWkaPYYBY\nRwmTEvzcqHsaUKKqNQAiMg2YxxejmxuTvtaUQa9SyEiTOzJy8qBbf1hrfZxM9Pn91cU3VncMIhBj\nImdPFZQvTJ/mvZhepa7pUjXsSIxpkJ8EdScwT0Qe9mpPc4GfBxuWMRGwbh5oTfolqKKxsHMTbP4s\n7EiMaZCfThKPicgsYDTuXqgbVHVD0IEZE7pYB4leI8ONo6UVeZ0jVr7tmvuMiShfTXyqul5Vn1fV\n5yw5mVZjTRl0/RK06xJ2JC2rW39o1xVWvRt2JMY0KLArvyKSIyJzROQDEflYRG731vcVkfdEZImI\nPCEihwUVgzEJU3UJKt2a9wBEXC1q1TthR2JMg4LsmrQHOFFVhwMlwCkiMhb4H+BeVe0PbAUuDzAG\nYxJTuQo+L4fepWFHEoyicW5k8yprEDHR1WCCEpEMEfkokQOrs8N7meU9YkMkPe2tnwacncjxjQlU\nrHbRe0y4cQQldh3KalEmwhpMUN69Tx+ISFEiBxeRNiIyHyjHTXi4FKhU1WpvkzVAr3r2nSQiZSJS\nVlFRkcjpjUncijfdDbrdByfldBdMfYcLpiYxWfQcBlntYKUlKBNdfm7U7Ql8LCJzgM9jK1X1zMZ2\nVNX9QImIdAKeAQbVtVk9+z4IPAhQWlpqN2yY5Fr5NhQdDRltwo4kGG2yXPOl1aBMhPlJULc39ySq\nWul1VR8LdBKRTK8W1RtY19zjG9OiqjbAlqVQ+h9hRxKsonHwxi9g93Y3woQxEdNoJwlVfR1YAWR5\nz/8NNDpOiojkezUnRKQt8FVgETATiM1dMBF4LqHIjQnKijfdsviYcOMIWtFYdyPymjlhR2JMnfwM\nFnsFrlPDVG9VL+BZH8fuCcwUkQW4pPaqNwr6DcC1IvIZ0BV4KJHAjQnMijchO88bWDWN9R4D0sY1\nZxoTQX6a+K4CxgDvAajqEhEpaGwnVV0AjKhj/TLveMZE07KZUHxs+l5/isnuAL1GwbJZMOGWsKMx\n5hB+7oPao6p7Yy9EJJN6OjYYk/K2LIetK6Df+LAjSY4jTnQjm+/cEnYkxhzCT4J6XURuBNqKyEnA\nU8A/gg3LmJAsm+mWR7SWBDUeUFj+RtiRGHMIPwlqMlABfAhcCUwHbg4yKGNCs3Qm5PV2Y/C1Br1G\nwWG5XyRmYyLEz2jmNd40G+/hmvY+UbWJZEwa2l8Ny1+HQWe48epagzZZ0PcrsPRfbvzB1lJukxL8\n9OI7DTcCxH3A/cBnInJq0IEZk3Sr34Pd26D/yWFHklz9T3JjD1Z8EnYkxhzATxPfL4HxqnqCqh4P\njAfuDTYsY0Lw6UuQkeU6DiTZ3uoaFq7bTnnV7qSfmyNPcctPX0r+uY1pgJ8EVa6q8VNvLsONrWdM\nevn0Jde9PDs36adeW7mLqj3V3PfakqSfm7xCd8+XJSgTMfVegxKRc72nH4vIdOBJ3DWo83E33hqT\nPrYsg02fQmlyZ38ZcPOL7KmuqX39yHureOS9VWRnZvDJlCS2pB95Csy+23U3T7cJGk3KaqgGdYb3\nyAE2AscDJ+B69HUOPDJjkmnRC2454JSknnb29eM5s6SQDK9vQk5WBmeVFDL7hiR3cx9wqhv26JMX\nk3teYxpQbw1KVS9LZiDGhGrhs9CzBDoXJ/W0BXk55GZnUuN1oNtTXUNudiYFuTlJjYPCEdCpyH0O\nIy5K7rmNqUej3cxFpC9wNVAcv72f6TaMSQlbV8LaufDV20I5/aYdeyjIzaYgN5uSos5UhNFRQgSO\nOgvefQB2bYW2hzaSxOareuLKccmOzrRSfsbiexY3oOs/gJpGtjUm9Sz0BtQ/KpzJnadeXFr7x3/K\n2UNCiQGAwefA27+BxdOtFmUiwU+C2q2q9wUeiTFh+fAp18TVpW/YkYSrcCR0Otx9HpagTAT46Wb+\naxG5VUTGicjI2CPwyIxJhg0fwYYFMPzbYUcSPhEY/i03uvm2tWFHY4yvBDUUuAK4C3fT7i+Bu4MM\nypik+eAxd3Pu0PMa37Y1GP4tQGHB42FHYoyvJr5zgH7xU24Ykxb274MFT7gu1nbvj9OlHxx+DMx7\nFI691sbmM6HyU4P6AOgUdCDGJN3C5+DzChh5SdiRRMuIi2HLUtfUZ0yI/CSo7sBiEXlZRJ6PPYIO\nzJjAzXnQ1RiOmBB2JNEy+Bxo19V9PsaEyE8T362BR2FMsq2b70Yv/9qdkOHn/2mtSFYOjLoUZt/j\nZhdO8s3LxsQ0+stU1dfrejS2n4j0EZGZIrJIRD4WkWu89V1E5FURWeItbdikiLhg6ju19+Okvbd+\nDYd1gBLrvVen0sshIxPevj/sSEwr5mc+qCoR2e49dovIfhHZ7uPY1cAPVXUQMBa4SkSOws3QO0NV\n+wMzvNfGJM+mJfDxMzDmCmhrl1fr1LEXlFwI7/8ZqjYAIU8JYlolPzWoXFXN8x45wDdwExc2tt96\nVX3fe14FLAJ6AWcB07zNpgHh3L5vWq/Zv4TMHBh7VdiRRNux/w9q9sFb7j79UKcEMa2Sn2tQB1DV\nZ0WkSbUeESkGRuCmje+uquu9Y60XkYKmxmBMwjZ8CB88Dkf/N3TIDzuaaOvSD4ZfyIBZo9kz65+1\nq0ObEsS0On4Giz037mUGUIqbF8oXEekA/A34gapuF5/3VYjIJGASQFFRkd/TGVM/VXjlZtes95Uf\nhh3NASI7AOuJNzN7wXim5PyAF7YWUaNuSpCvDe7BTacNCjs6k+b8dF86I+7xNaAK10zXKBHJwiWn\nR1X1797qjSLS03u/J/XMzquqD6pqqaqW5ufb/3RNC1j4rLu35/jJdY7WbeqQV0jBsZeQu+0TalTD\nnRLEtDqN1qASnRdKXFXpIWCRqt4T99bzwETc0EkTgecSOb4xTbJzC0y/zs35NPq7YUeTWr5yLZve\n+A3fkLdYVvBVBhd1D2dKENPqNDTl+y0N7Keqekcjxz4GuBj4UETme+tuxCWmJ0XkcmAVbgp5Y4Kj\nCv+4xs1zdPEz0KbJl16bJO3mTcpqy9RLj4aHT+O13cv56tnBj9OXdp+hSUhDv9TP61jXHrgc6Ao0\nmKBU9U2gvgtOduu+SZ73psKi5+GkO6DH0LCjSU3Fx/Js+29y9udPwoInYdg3w47ItAINTfn+y9hz\nEckFrgEuAx7HjWhu0sze6ho+K99BedXu9Lm+8OnL8PKNMODrcPTVYUeT0p7InciAfQsZ9Nx/Q8c+\ncLjVbkywGuwk4Y36MAVYgEtmI1X1BlWts2ODSW1pd5/LijfhqUtdrenc39vI3M1UI234ZeefQKc+\n8Ni3YN28sEMyaa7eBCUivwD+jeu1N1RVb1PVrUmLzCTNgJtfpHjyPymv2gO4+1yKJ/+TATe/GHJk\nzbD4n/DIN6BTEVz0FGR3CDuitFCV0RG+8zfIzoNpZ8KyRkc9MyZhDdWgfggUAjcD6+KGO6ryOdSR\nSRGzrx/PmSWFZHgVjJysDM4qKWT2DePDDSwR+/fBjDvg8W9DwVFw6XToYPeCt6jOxfAfL0FuT/jL\nOW5cw5r9YUdl0lC9CUpVM1S17UFDHeXFXiczyHQR1cFYC/JyyM3OpEZJ7ftcVrwFD46H2XdDyXfg\nsunQvisQ3c8+ZXXsBd99DQZ+HV69Bf54Cqx9P+yoTAuI0m8l2P62JmVs2rGHgtxsCnKzKSnqnDr3\nuVTvhcUvQNkfYcVsyC2ECx6FQaeHHVn6y8mDb/7F9ep7+cfw+/HQ/2QYdZlbBtyd36Q/+wYZAKZe\nXFr7v6YpZw8JOZoGqML2da4DxPLXXS+9nZugYxGc9FMYfQUc1i7UENOyN2R9RGD4BTDgVHj3t1D2\nJ3j8Qtf81/9k6Hc8HH4s5HZv0mFb1Wdo6tU6EpRqw6/dysb3q2u7JhwrS/e65/t2J36+OrdrmWO1\nq9mBoO6G1kaPVdfhmxbXFX8uA+D3F4+Afbugerdb7tsF1bvc6A87yuHzcthRAVuWsWPNR3TQHe4A\nbTtDvxNcc94R4yGjjb84fUr0ZtH43pBTzkmf+64a/Bxy8uCEyW6Mw09fhg8ec1OavO9NXNCuK+QP\nhM593SC97QvctcG2nSCrnRtdPqudmywxsy3lWyup2gP3vfwRU07vD5JR9wNpFb0zE/0uJrRf7Dd6\n8O85hM9Z1O8fnxCVlpZqWVlZ4gf4WU/Yt7PlAjLJlZHl/ph1OpxXNnVmbWYRl114IXQf6ns23ER+\nqE3dZ8DNL7KnuuaQ9a121O/91bD+A1j9LlQshopPoHK1+w9HTXWduwzY/TB7OOyQ9dns5ZOcSwMO\n2PiS1xuu/bhZhxCRuapa2th2raMG9ZVr6+hlVMf/Bur8H0Jd2/nZ5tB1f52zGoBvjylq4vl8bNcC\nx5r2zgoAJh7dN+C4nIfeXA7A5V/pF/e/6LZfLNt2cf/bzulUe5yHvKRxWc/h9R43LLOvH8+U6Yt4\n4YN1Nuo3uGtQvUe5R7yaGthd6WrHuysPqDXP3r6LKfOqmb46m2rakNOmhq/13MlNR22G7J+C1sQ9\n1C1bgafK3N+O80v7BL5fnfvEV2RyktdHrnUkqOOuCzsCAJ5b6P64fvsr0bwDf/oCF9/EscmJ75V5\n7nyXj47m59FUadMbMmgZGdCui3scpADI3fgh1atXuc+wJoPcXkdRcGL6NJUm4ulP3G/l/BOa9ltJ\nZL9EzxWE1pGgWlii7cF24Tc8yfrsU7Y3ZIQk+zO0gWkPFKW/U/4a8E2LSLuhhFJIsj77qReX0rdb\ne9pnZzLl7CFMvbjRZnZzEPsMwxWlv1NWg0qCgy+e25TZyWOfvTH+RPG3YjWoJEiroYRa0N7qGhau\n2055E5pwmrpPcz77ROIzh4rSyASpKtHvYlP2i+LfKUtQSZAqF8+fuHJcUtvhE2lKaOo+zfnso9TU\nYVq3RL+LTdkvin+nrIkvSezi+RcSaUpoTvNDUz/7KDZ1mNYp0e9iovtF7e+U1aASkEh1O50v/Da1\nCSeRpoTmND809bOPYlOHCU6izWfJaLpM9LuY6H5R+ztlCSoB1vTTPIk0JSSz+SGKTR0mOFH+PSf6\nXUyX77A18TWBNf20nESaEpLZ/BC1pg7T8lLl95zodzEdvsOBJSgR+SNwOlCuqkO8dV2AJ4BiYAXw\nzVSapdeGsmk5iYyenswR15tzLrvh80CJ3PiZjM8wVX7PiX4XU2aGggYE2cT3MHDKQesmAzNUtT8w\nw3udMtKl2mxMMiWzCa0p14Xs9xx9gdWgVPUNESk+aPVZwAne82nALOCGoGIIQjpUm41JhlRoQrPf\nc7Qlu5NEd1VdD+AtC+rbUEQmiUiZiJRVVFQkLcDGRK2XSxTYDa2mLqnQG7I5v2f73gcvsr34VPVB\nVS1V1dL8/PzAzhP1u9wTjS+Z5YpyLygTnnRvQkv0ex/1vzlRkuxefBtFpKeqrheRnkB5ks8fqnS7\neJ4KTTgx6fbZp4p0bEJLpe99IqL0W0l2Dep5YKL3fCLwXJLPb1pQKjThmHAlu0k8Gc1u9r1PniC7\nmT+G6xDRTUTWALcCdwFPisjlwCrg/KDOb4LX3CacRP6nlsz/3UXpf5LGn/hmtynnBDPJYRhNl4l+\nF1P9OxxkL74L63lrQlDnTESUJueqS6Lx2QR9pjVJdrNbc773Uf+bEyWR7SSRLIlc6EzmqN/JGMW4\nOaxXo4mC5jS7JfJ7bs733joV+ddqhzqK+oXOZI9ibEwqS4Ueg/bbbLpWW4OK+oXOZI9ibEyqizW7\nDe6Zx0VfPpyKHXvCDukA9ttsulZbg4r6/7ha+yjGJn0kqzk86mPP2W+z6VptgoLoX+BvzaMYG5OO\n7LfZNKKqYcfQqNLSUi0rKwvk2LH/cUW1O2ai8UW9XMYEIRW+96kQY9BEZK6qNtqzpNVegzLGGBNt\nlqCMMcZEUqu+BmWMSS+tudksHVkNyhhjTCS1+k4Sxhhjkss6SRhjjElplqCMMcZEkiUoY4wxkWQJ\nyhhjTCRZgjLGGBNJlqCMMcZEkiUoY4wxkWQJyhhjTCRZgjLGGBNJKTGShIhUACubeZhuwKYWCCeq\n0r18kP5lTPfygZUxHbRE+Q5X1fzGNkqJBNUSRKTMz9AaqSrdywfpX8Z0Lx9YGdNBMstnTXzGGGMi\nyRKUMcaYSGpNCerBsAMIWLqXD9K/jOlePrAypoOkla/VXIMyxhiTWlpTDcoYY0wKsQRljDEmktIi\nQYnIH0WkXEQ+ilvXRUReFZEl3rKzt15E5D4R+UxEFojIyPAi909E+ojITBFZJCIfi8g13vq0KKeI\n5IjIHBH5wCvf7d76viLynle+J0TkMG99tvf6M+/94jDjbwoRaSMi80TkBe912pRRRFaIyIciMl9E\nyrx1afEdjRGRTiLytIgs9n6P49KpjCIywPv3iz22i8gPwihjWiQo4GHglIPWTQZmqGp/YIb3GuBU\noL/3mAT8LkkxNlc18ENVHQSMBa4SkaNIn3LuAU5U1eFACXCKiIwF/ge41yvfVuByb/vLga2q+iXg\nXm+7VHENsCjudbqVcbyqlsTdK5Mu39GYXwMvqepAYDju3zJtyqiqn3j/fiXAKGAn8AxhlFFV0+IB\nFAMfxb3+BOjpPe8JfOI9nwpcWNd2qfQAngNOSsdyAu2A94Ev4+5Yz/TWjwNe9p6/DIzznmd620nY\nsfsoW2/cj/tE4AVA0qmMwAqg20Hr0uY7CuQByw/+d0inMh5UrpOBt8IqY7rUoOrSXVXXA3jLAm99\nL2B13HZrvHUpw2vqGQG8RxqV02v6mg+UA68CS4FKVa32NokvQ235vPe3AV2TG3FCfgVcD9R4r7uS\nXmVU4BURmSsik7x1afMdBfoBFcCfvGbaP4hIe9KrjPG+BTzmPU96GdM5QdVH6liXMn3tRaQD8Dfg\nB6q6vaFN61gX6XKq6n51zQq9gTHAoLo285YpVz4ROR0oV9W58avr2DRlywgco6ojcc0+V4nIcQ1s\nm4rlywRGAr9T1RHA53zR1FWXVCwjAN610DOBpxrbtI51LVLGdE5QG0WkJ4C3LPfWrwH6xG3XG1iX\n5NgSIiJZuOT0qKr+3VudduVU1UpgFu5aWycRyfTeii9Dbfm89zsCW5IbaZMdA5wpIiuAx3HNfL8i\njcqoquu8ZTnuusUY0us7ugZYo6rvea+fxiWsdCpjzKnA+6q60Xud9DKmc4J6HpjoPZ+Iu2YTW3+J\n1/NkLLAtVm2NMhER4CFgkareE/dWWpRTRPJFpJP3vC3wVdzF55nAed5mB5cvVu7zgH+p1wAeVar6\nY1XtrarFuKaTf6nqRaRJGUWkvYjkxp7jrl98RJp8RwFUdQOwWkQGeKsmAAtJozLGuZAvmvcgjDKG\nfRGuhS7kPQasB/bhsvnluLb6GcASb9nF21aA/8Nd3/gQKA07fp9lPBZXbV4AzPceX0+XcgLDgHle\n+T4CbvHW9wPmAJ/hmhqyvfU53uvPvPf7hV2GJpb3BOCFdCqjV44PvMfHwE3e+rT4jsaVswQo876r\nzwKd07CM7YDNQMe4dUkvow11ZIwxJpLSuYnPGGNMCrMEZYwxJpIsQRljjIkkS1DGGGMiyRKUMcaY\nSLIEZYwxJpIsQRljjIkkS1DGhExERnvz6OR4ozF8LCJDwo7LmLDZjbrGRICITMGNHNEWN9bbnSGH\nZEzoLEEZEwHeyNH/BnYDR6vq/pBDMiZ01sRnTDR0AToAubialDGtntWgjIkAEXkeNwVHX9xspP8d\nckjGhC6z8U2MMUESkUuAalX9q4i0Ad4WkRNV9V9hx2ZMmKwGZYwxJpLsGpQxxphIsgRljDEmkixB\nGWOMiSRLUMYYYyLJEpQxxphIsgRljDEmkixBGWOMiaT/D+4rshd1/oImAAAAAElFTkSuQmCC\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", "mu = 400\n", "sigma = 30\n", "B = 15\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=\"Initial Fit 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": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "A = 34.954073 +- 4.652275\n", "mu = 431.855370 +- 3.217831\n", "sigma = 25.495200 +- 3.015018\n", "B = 10.088426 +- 0.608374\n", "Covariance matrix\n", " [[ 2.16436614e+01 3.81473839e-01 -7.96170891e+00 -1.63022126e-01]\n", " [ 3.81473839e-01 1.03544389e+01 -2.13828006e-01 -9.40358913e-03]\n", " [ -7.96170891e+00 -2.13828006e-01 9.09033599e+00 -4.79433728e-01]\n", " [ -1.63022126e-01 -9.40358913e-03 -4.79433728e-01 3.70118843e-01]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl8VNX5+PHPk4WEJWEPOyQqoEAkQFAQq1XqVhGXatX6\nU2v9frXVVru5VevS0uq3tvWrrVVs7Vct7vveFhUUhaKgKDsIsi9J2JKwBJI8vz/unTBAlpnJ3Ll3\nZp736zWvO3PnLs8ZZvJwzj33HFFVjDHGmKDJ8DsAY4wxpjGWoIwxxgSSJShjjDGBZAnKGGNMIFmC\nMsYYE0iWoIwxxgSSJShjPCIij4nIpCbe+66IfOjBOaeLyH/F+7he8+rzMMnNEpQJDBFZJSK7RaRK\nRLaLyEwR+b6IRPQ9FZFCEVERyfI6VmOM9yxBmaA5S1XzgAHAPcBNwKP+hpR+LMmbILAEZQJJVXeo\n6mvAhcDlIjIMQETOFJHPRKRSRNaKyJ1hu33gLreLSLWIjBWRw0XkPRHZIiIVIvKkiHRq6rwicr97\n3EoRmSsiXwt7704ReU5EnnBreQtFpDTs/REi8qn73rNAbgvFFBH5k4jsEJElIjI+7I0rRGSxe6yV\nInL1QTueLSLz3DhXiMjpjRy8l4h8ISI/d18XicgH7jHfEZEHRWSK+16o9nmliKwB3nPXT3TLud1t\nPjwq7PgqIkeEvW5o0hSRr4vIOhH5mYiUichGEbkibNuuIvKaG//HwOEtfFYmDVmCMoGmqh8D64BQ\notgJXAZ0As4EfiAi57jvneAuO6lqB1WdBQhwN9AbOAroB9zZzCk/AUqALsBTwPMiEp5oJgLPuOd/\nDfgzgIi0AV4B/uHu+zzwrRaKdyywEugG3AG8JCJd3PfKgAlAPnAFcJ+IjHTPdQzwBHCDG8cJwKrw\nA4tIIfA+8GdV/b27+ingY6Cr+xlc2khMJ+J8TqeJyCDgaeDHQHfgLeB1t6yR6Al0BPoAVwIPikhn\n970HgT1AL+B77sOYA1iCMslgA84ffVR1uqrOV9V6Vf0C5w/oiU3tqKpfqupUVa1R1XLgjy1sP0VV\nt6hqrar+AcgBBodt8qGqvqWqdTjJaLi7fgyQDfyvqu5T1Rdwkl1zysK2fxZYipN0UdU3VXWFOt4H\n/s3+JH0l8He3XPWqul5Vl4QddwgwHbhDVR8BEJH+wGjgdlXdq6of4iTYg92pqjtVdTdO7fVN9zz7\ngN8DbYHjWihXyD7gV2753gKqgcEikomTvG93z7UAeDzCY5o0YgnKJIM+wFYAETlWRKaJSLmI7AC+\nj1MDaZSIFIjIMyKyXkQqgSktbP8zt2lth4hsx6kBhG+/Kez5LiDXvV7TG1ivB46+vLqFcjW2fW83\njjNE5D8istWN45thcfQDVjRz3EuA9cALYet6A1tVdVfYurWN7Bu+rnd4GVS13n2/T7Ol2m+LqtaG\nvd4FdMCpjWUddK6WPiuThixBmUATkdE4fxBDXZCfwvmffz9V7Qg8jNOMB9DY0Px3u+uPVtV84P+F\nbX/wub6G0ynj20BnVe0E7Ghq+4NsBPqISPi2/VvYp7HtN4hIDvAiTo2lhxvHW2FxrKX5azZ3AhXA\nU25tJRRfFxFpF7Zdv0b2Df8MN+B0VgGcC2buPuvdVbuA8OP1bCamcOVA7UHnb+mzMmnIEpQJJBHJ\nF5EJONd7pqjqfPetPJyawB73Wsx3wnYrB+qBw8LW5eE0LW0XkT44122akofzh7McyBKR23GuAUVi\nlrvvdSKSJSLnAce0sE+Bu322iFyAc+3nLaANTtNiOVArImcAp4bt9yhwhYiMF5EMEekjIkeGvb8P\nuABoD/xDRDJUdTUwB7hTRNqIyFjgrBbiew440z1PNvAzoAaY6b4/D/iOiGS6nTSabDoN5zaPvuTG\n0k5EhgCXR7KvSS+WoEzQvC4iVTi1hFtxrhldEfb+NcCv3G1ux/kjCoDbfPUb4CO319kY4C5gJE5N\n6E2cP4xN+RfwNrAMp8lpD403gx1CVfcC5wHfBbbhXL9p7lwAs4GBOLWd3wDnu9e/qoDr3LJtw0nC\nDdeL3I4jVwD3ueV6n7CazkHxFAB/F+deskuAscAWYBLwLE7CaapMS3FqnH9yYzwL5zaAve4m17vr\ntrvHfqWF8ob7IU5z3ybgMeD/otjXpAmxCQuNSU9uV/glqnqH37EY0xirQRmTJkRktDj3hWW4TXJn\nE12tx5iEsrvFjUkfPXGaHbvi3Fv2A1X9zN+QjGmaNfEZY4wJJGviM8YYE0hJ0cTXrVs3LSws9DsM\nY4wxcTB37twKVe3e0nZJkaAKCwuZM2eO32EYY4yJAxGJaOQQa+IzxhgTSJagjDHGBJIlKGOMMYHk\n6TUoEVkFVAF1QK2qlrrz3TwLFOLMYfNtVd3mZRzGmOS2b98+1q1bx549e/wOxUQhNzeXvn37kp2d\nHdP+iegkcZKqVoS9vhl4V1XvEZGb3dc3JSAOY0ySWrduHXl5eRQWFnLgAPAmqFSVLVu2sG7dOoqK\nimI6hh9NfGezf3Kyx4FzmtnWGGPYs2cPXbt2teSURESErl27tqrW63WCUuDfIjJXRK5y1/VQ1Y0A\n7rKgsR1F5CoRmSMic8rLyz0O0xgTdJackk9r/828buIbp6obRKQAmCoiS1rcw+VOVf0IQGlpqY3H\nZIwxacbTGpSqbnCXZcDLOBO4bRaRXgDusszLGIxJZRdOnsWFk2f5HUZayMzMpKSkpOGxatUq5syZ\nw3XXXQfA9OnTmTlzZsP2r7zyCosWLYr6PB06dIj4/PGyfft2/vKXvzS83rBhA+eff37cjh8rz2pQ\nItIeyFDVKvf5qcCvcCZeuxy4x12+6lUMxhgTL23btmXevHkHrCssLKS0tBRwElSHDh047rjjACdB\nTZgwgSFDhnh2/ngJJahrrrkGgN69e/PCCy94cq5oeFmD6gF8KCKfAx8Db6rqP3ES0ykishw4xX1t\njDFJZ/r06UyYMIFVq1bx8MMPc99991FSUsL777/Pa6+9xg033EBJSQkrVqxgxYoVnH766YwaNYqv\nfe1rLFniXPH46quvGDt2LKNHj+aXv/xlVOd/7LHH+OEPf9jwesKECUyfPh1wamK33norw4cPZ8yY\nMWzevBmAzZs3c+655zJ8+HCGDx/OzJkzufnmm1mxYgUlJSXccMMNrFq1imHDhgFOB5UrrriC4uJi\nRowYwbRp0xrOfd5553H66aczcOBAbrzxxtZ+nIfwrAalqiuB4Y2s3wKM9+q8xpgU9/bNsGl+fI/Z\nsxjOaP7/yrt376akpASAoqIiXn755Yb3CgsL+f73v0+HDh34+c9/DsDEiROZMGFCQ1PZ+PHjefjh\nhxk4cCCzZ8/mmmuu4b333uP666/nBz/4AZdddhkPPvhgTOdvzM6dOxkzZgy/+c1vuPHGG/nrX//K\nbbfdxnXXXceJJ57Iyy+/TF1dHdXV1dxzzz0sWLCgoYYW3nwYimn+/PksWbKEU089lWXLlgEwb948\nPvvsM3Jychg8eDA/+tGP6NevX7NxRSMpBos1xhi/taaJrbq6mpkzZ3LBBRc0rKupqQHgo48+4sUX\nXwTg0ksv5aabGr8tNNrzt2nThgkTJgAwatQopk6dCsB7773HE088ATjXtTp27Mi2bU2PlfDhhx/y\nox/9CIAjjzySAQMGNCSo8ePH07FjRwCGDBnC6tWrLUEZY9JYCzWdIKqvr6dTp05NJphYu2NnZWVR\nX1/f8Dr8nqPs7OyG42ZmZlJbWxvTOZqb1DYnJ6fheWvO0RQbi88YY+IgLy+PqqqqRl/n5+dTVFTE\n888/Dzh/9D///HMAxo0bxzPPPAPAk08+GdU5CwsLmTdvHvX19axdu5aPP/64xX3Gjx/PQw89BEBd\nXR2VlZWHxB7uhBNOaIhr2bJlrFmzhsGDB0cVZ6wsQRljTBycddZZvPzyy5SUlDBjxgwuuugi7r33\nXkaMGMGKFSt48sknefTRRxk+fDhDhw7l1VedDsz3338/Dz74IKNHj2bHjh1RnXPcuHEUFRVRXFzM\nz3/+c0aOHNniPvfffz/Tpk2juLiYUaNGsXDhQrp27cq4ceMYNmwYN9xwwwHbX3PNNdTV1VFcXMyF\nF17IY489dkDNyUvSXPUtKEpLS9UmLDTmUKF7oJ69eqzPkXhr8eLFHHXUUX6HYWLQ2L+diMxV1dKW\n9rUalDHGmECyBGWMMSaQLEEZY4wJJEtQxhhjAskSlDHGmECyBGWMSUk20nvyswRlTBqyP97RC013\nMXz4cEaOHNkwtUaipqYoLCykoqKi0fXFxcUUFxczZMgQbrvttoZhlJpy8PQaQWUJyhiTkvbW1rNo\nQyVlVbFPOR4uNBbe559/zt13380tt9wCBGNqimnTpjF//nw+/vhjVq5cyVVXXdXs9pagjDHGR+u3\n76aqppYH3lke92NXVlbSuXNngAOmpmhuCoqmpr8oLy/nW9/6FqNHj2b06NF89NFHAGzZsoVTTz2V\nESNGcPXVVzc7Jl74OR5++GFeeeUVtm7dSnV1NePHj2fkyJEUFxc3jF5x8PQaTW3nNxss1hiTUgbf\n9jY1tfsHUJ0yew1TZq8hJyuDpZPOiPm4oeku9uzZw8aNG3nvvfca3a6pKSiamv7i+uuv5yc/+QnH\nH388a9as4bTTTmPx4sXcddddHH/88dx+++28+eabPPLIIxHFGRr3b/ny5YwaNYqXX36Z/Px8Kioq\nGDNmDBMnTjxkeo3a2tpGt4t1ENt4sQRlTCrZsR7Wz4VeR0PnQr+j8cWMG09i0luLeePzDdQr5GZn\ncNrQntx6ZuuGSgqf7mLWrFlcdtllLFiw4JDtmpqCoqnpL955550DpoavrKykqqqKDz74gJdeegmA\nM888s6HGFolQbUtV+cUvfsEHH3xARkYG69evb6i5Hbx9Y9v17Nkz4nN6wRKUManii+fh5atA66FN\nBzj3YTjqLL+jSriC/FzycrKoVxCBmtp68nKyKMjLjds5xo4dS0VFBeXl5Ye819QUFE1Nf1FfX8+s\nWbNo27btIceKpQZTVVXFqlWrGDRoEE8++STl5eXMnTuX7OxsCgsLD5iSIyTS7RLNrkEZkwo2fAav\n/RD6jYHL34Dug+H570JF/K+/JIOK6hoK8nIY2iufS44dQHl1873aorVkyRLq6uro2rVrq4916qmn\n8uc//7nhdaiWFj7Nxdtvv93spIIh1dXVXHPNNZxzzjl07tyZHTt2UFBQQHZ2NtOmTWP16tXAoVOD\nNLWd36wGZUwqeOtGaNsFLvwHtO8GFz8DD4yEqbfDxU/7HV3CTb60tKEb/aRzhsXlmOFTrqsqjz/+\nOJmZma0+7gMPPMC1117L0UcfTW1tLSeccAIPP/wwd9xxBxdffDEjR47kxBNPpH///k0e46STTkJV\nqa+v59xzz+WXv/wlAJdccglnnXUWpaWllJSUcOSRRwIcML3GGWecwU033dTodn6z6TaMSWIXTp7F\n4L0L+dWWn8EZ98KxYd2LP7wP3rkTrngbBhx3yH6QPNN0xDLdRrKVMVW1ZroNq0EZk+TOqn4B2naG\nEZcc+Max34cZ98GnTxySoNKBJabkZ9egjEliHeu2MqrmP1D6PWjT/sA3s9vC0HNg0WtQU+1PgMa0\ngiUoY5LYkTvncFHNbZQVnd34BsMvhn07YfHriQ3MA8lwOcIcqLX/ZpagjEliyyoz+UQH88DnTfyU\n+49x7oda4O9QPK2Vm5vLli1bLEklEVVly5Yt5ObG3r3frkEZk4T2j5YwAmhmtAQRGHQGzP0/2Lfb\nafZLQn379mXdunWN3ndkgis3N5e+ffvGvL8lKGOS0IwbT2LSk//k36vr2UNO86MlHPENmP0QrP7I\neZ6EsrOzKSoq8jsMk2DWxGdMEirIzyVvzwZqyEZoYbSEAcdBZg582fjYccYElSUoY5JUxfYqTmmz\ngKG9WxgtoU07KBwHX76T2ACNaSVLUMYko+1rmCy/ZWyn7bTPyWLSOcOYfGkz9z0efjJULIXKjYmL\n0ZhWsgRlTDL6agYAC9sMj2z7/u6NumtnexSQMfFnCcqYZLTqQ2jXlXVZAyLbvtfRkNXWEpRJKpag\njElGa2dDvzGoRPgTzsyGPqNgzX+8jcuYOLIEZUyy2bUVtq6AvqOi26/fMbDpC9i7y5u4jIkzS1DG\nJJv1nzrLPi0OBn2g/mOgvhY2fBr/mIzxgCUoY5LN+jmAQO8R0e3Xd7SzXPdJ3EMyxgstJigRGSci\nU0VkmYisFJGvRGRlpCcQkUwR+UxE3nBfF4nIbBFZLiLPikib1hTAmLSzfi50PxJy86Pbr10XZ1y+\nDfM8CcuYeIukBvUo8EfgeGA0UOouI3U9sDjs9f8A96nqQGAbcGUUxzImvanCujnRX38K6VUCG+ex\nt7aeRRsqKavaE9/4jImjSBLUDlV9W1XLVHVL6BHJwUWkL3Am8Df3tQAnA6GhlR8HzokhbmPS0/Y1\nsHur0yMvFr2Gw7ZVbN62g6qaWh54Z3l84zMmjiIZLHaaiNwLvAQ0jKWiqpFcaf1f4EYgz33dFdiu\nqrXu63VAn8jDNSbNbfrCWfaM8Abdgwz+52Bq6p4CnGkrmhwF3ZgAiKQGdSxOs95vgT+4j9+3tJOI\nTADKVHVu+OpGNm10ghcRuUpE5ojIHBti3xjXpvkgGdBjSEy7z7huFBMzPiKbOgByszM4u6Q3M246\nKZ5RGhMXLdagVDXWb+44YKKIfBPIBfJxalSdRCTLrUX1BTY0cd5HgEcASktLbZYyYwA2fgHdBsU8\nr1NBj57k5WRSuzsDkRZGQTfGZ00mKBH5f6o6RUR+2tj7qvrH5g6sqrcAt7jH+jrwc1W9RESeB84H\nngEuB16NMXZj0s+m+TBgbKsOUdGmD+fum8mygjMo6d+ZcusoYQKquRpUe3eZ18w2sbgJeEZEJgGf\n4fQSNMa0ZNdWqFwHPY9u1WEmj9lG/fsPcVmbU5l0zrA4BWdM/DWZoFR1sru8q7UnUdXpwHT3+Urg\nmNYe05i009BBorh1xyk4igzq6VO7tvUxGeOhSG7UPUxEXheRchEpE5FXReSwRARnjAmzaYGzbGUN\nigKng0X/fatadxxjPBZJL76ngOeAXkBv4HngaS+DMsY0onwJtOsG7bu27jhdDmMv2fSrXRWXsIzx\nSiQJSlT1H6pa6z6m0ETXcGOMhyqWQ/fBrT9OZhbrs/rTr3Z1649ljIeaTFAi0kVEuuDcqHuziBSK\nyAARuRF4M3EhGmNQdaZs7zYoLodbmz2Afvu+isuxjPFKc7345uLUlEI3114d9p4Cv/YqKGPMQXZW\nwO5t8alBAWuzCjlh93uwezu07RSXYxoTb8314itKZCDGmGZULHWW8apBhaaKL1/izBNlTADZfFDG\nJINyN0HFqQa1Jtv9/2fZorgczxgvWIIyJhlULIPs9pAfn7GVt2R0Z5e0g82WoExwNZugxNEvUcEY\nY5pQvhS6DQRpbLzlGIg4zXxli1ve1hifNJugVFWBVxIUizGmKRXL4ta8F7I2q9Bp4lO7a8QEUyTz\nQf1HREar6ieeR2OMOVRNFVSub7SDxLNXxz5w7NrsAVD5NlRvhryerYnQGE9EkqBOAr4vIquAnTjd\nzlVVWzneijEmIhXLnGXca1BhHSUsQZkAiiRB2TSbxvipwp2WvVt8E9SaUFfzzYvg8JPjemxj4qHF\nXnyquhroB5zsPt8VyX7GmDgpXwoZWdAlvrcmVmV2gnZd999jZUzARDKa+R04czjd4q7KBqZ4GZQx\nJkzFMuhyOGRmx//Y3QZD+bL4H9eYOIikJnQuMBHn+hOquoH4T2JojGlK+VLoHp8RJA7RfdD+a1zG\nBEwkCWqv291cAUSkfQvbG2PipXYvbF0Z9+tPDboNgt1bnbH+jAmYSBLUcyIyGegkIv8NvAP81duw\njDGAk5y0Lm5j8B0ilPjK7TqUCZ4We/Gp6u9F5BSgEhgE3K6qUz2PzBizvwODl0184DTzFY7z5hzG\nxCiSbuYA84G2OM18870LxxhzgFAHBq9qUPl9IbudXYcygRRJL77/Aj4GzgPOxxlZ4nteB2aMwalB\ndewHbTy69JuRAV2PsCY+E0iR1KBuAEao6hYAEekKzAT+7mVgxhjcQWI9qj2FdB8Ma2Z7ew5jYhBJ\ngloHVIW9rgLWehOOMaZBfb0zikTh8XE/9AFj+HUbBPOfh707vaupGRODSBLUemC2iLyKcw3qbOBj\nEfkpgKr+0cP4jElfleugdrf3NajQ8bd8Cb2Ge3suY6IQSYJa4T5CXnWXdrOuMV4q92aQ2EOEElT5\nMktQJlAi6WZ+VyICMcYcJNTF3KubdEO6Hg6SYWPymcCxQV+NCarypdC2C7Tv6u15snKgc5F1NTeB\nYwnKmKDyYBbdJnUbZIPGmsCxBGVMUCWii3lI90FOJ4m62sScz5gIRHKj7u9EJF9EskXkXRGpEJH/\nl4jgjElbOyucQVwTVoMaDPX7YPvqxJzPmAhEUoM6VVUrgQk490QNwrl51xjjlfIEdZAIaejJZx0l\nTHBEkqBCs6R9E3haVbd6GI8xBrwfJPZg3Qa657XrUCY4IrkP6nURWQLsBq4Rke7AHm/DMibNlS9z\nBnHN75uY87XtBB16WIIygRJJDeoOYCxQqqr7gF04M+waY7xSscyp1WQksB9Tt0HWxGcCJZJv/yxV\n3aaqdQCquhN429uwjElzFcsSd/0ppPtgZ+w/1cSe15gmNNnEJyI9gT5AWxEZAYj7Vj7QrqUDi0gu\n8AGQ457nBVW9Q0SKgGeALsCnwKWqurdVpTAmldRUw4610P3yxJ632yCo2QHVmyGvZ2LPbUwjmrsG\ndRrwXaAvED4gbBXwiwiOXQOcrKrVIpINfCgibwM/Be5T1WdE5GHgSuChWII3JiVtWe4sE3UPVEh4\nTz5LUCYAmkxQqvo48LiIfEtVX4z2wKqqQLX7Mtt9KHAy8B13/ePAnViCMma/hll0fWjiA6d58bAT\nE3tuYxoRSS++N0TkO0Bh+Paq+quWdhSRTGAucATwIM6o6NtVNXS7+jqcZkRjTEjFUpBM6HJYYs+b\n1wva5FlPPhMYkSSoV4EdOImmJpqDux0rSkSkE/AycFRjmzW2r4hcBVwF0L9//2hOa0zSuXDyLMCd\nSLB8qZOcstokNggRp+egJSgTEJEkqL6qenprTqKq20VkOjAG6CQiWW4tqi+woYl9HgEeASgtLbVu\nRSZ9JHKQ2IN1GwRffeDPuY05SCTdzGeKSHG0BxaR7m7NCRFpC3wDWAxMA853N7uc/RMgGmPq9sHW\nlYnvIBHSfRBUbYA9lf6c35gwkSSo44G5IrJURL4Qkfki8kUE+/UCprnbfgJMVdU3gJuAn4rIl0BX\n4NFYgzcm5WxdCfW1Ptag3POGehIa46NImvjOiOXAqvoFMKKR9SuBY2I5pjEpL3T9x68aVPj0731G\n+RODMa4Wa1Cquhroh3NP02qcoY5sHiljvNAwirlPCapLEWRk2fTvJhAimQ/qDpxmuVvcVdnAFC+D\nMiZtVSyD/D6Q08Gf82dmQ5fDnSGPjPFZJDWhc3EGh90JoKobgDwvgzImbSVyFt2mdBtog8aaQIgk\nQe11R4VQABFp721IxqQn0Xqn5uJXB4mQ7oOdzhq1NkSm8VckCeo5EZmMc//SfwPvAH/1Nixj0k+X\n+grYtzMANahBoHWw7St/4zBpr8VefKr6exE5BagEBgO3q+pUzyMzJs30rV3jPPG7BhU+aKzfsZi0\n1mKCEpGfAM9bUjLGW31q1zpPEj1I7MFCCcqGPDI+i6SJLx/4l4jMEJFrRaSH10EZk4761K6Btp2h\nfTd/A8np4PQktARlfBbJfVB3qepQ4FqgN/C+iLzjeWTGpJk++9ZA9yOdQVv9ZtO/mwCI5obbMmAT\nsAUo8CYcY9KUKv1qVzsJKghs+ncTAJHcqPsDdyTyd4FuwH+r6tFeB2ZMOulUv5UOWg0Fjc1I44Nu\nA50ehZXr/Y7EpLFIxuIbAPxYVed5HYwx6apf7WrnSVBqUKGOGuVLoWNff2MxaSuSa1A3Ax1E5Apo\nmEajyPPIjEkjDV3MA1ODCvXksyGPjH9sLD5jAqDvvtVUSR607+53KI4OBZDb0QaNNb6ysfiMCYB2\ne7dw8d5fUFZd43coDhGnmc9qUMZHNhafMX5T5dU9JSys688D7wQoIVhXc+OzSDpJHDwW3/ewsfiM\niYvBt71NTW098HUApsxew5TZa8jJymDppJjmCo2f7oNg3hTYvc25gdiYBIukk8TvgReAF9k/Ft+f\nvA7MmHQw48aTmHh4Jrk4TXu52RmcXdKbGTed5HNk7O/JZ818xieR1KBwx+GzsfiMibOC/Fzy6rZR\nQwcygJraevJysijIy/U7NOdeKHCa+fod428sJi1FlKCMMd6pqNzFtzM/ZWGPiZT070x51R6/Q3J0\nLoTMNg09+S6cPAuAZ68e62NQJp1YgjLGZ5O7PMOi2t3clXMek84Z5nc4+2VkQlebXdf4p8lrUCLy\nrrv8n8SFY0yaUYXyJazLGuB3JI3rMRQ2L/Q7CpOmmqtB9RKRE4GJIvIMcMAQy6r6qaeRGZMOqjZC\nTSVr8/v7HUnjeg6D+c/Brq1+R2LSUHMJ6nbgZqAv8MeD3lPgZK+CMiZtlC0GCHANym1y3LwAZxAZ\nYxKnyQSlqi8AL4jIL1X11wmMyZj0sXkBAGuyAzq8Zc9iZ7lpPjDS11BM+mmxk4Sq/lpEJgInuKum\nq+ob3oZlTJrY+AXk96U6I9/vSBrXoQDaF8CmBViCMokWyWCxdwPXA4vcx/XuOmNMa22aD70CPr1a\nz2Gweb7fUZg0FMlYfGcCp6jq31X178Dp7jpjTGvs3QVblu9vRguqHsOgfCmZWut3JCbNRDrle6ew\n5x29CMSYtLN5IWg99Ax6DaoY6vbSu3at35GYNBPJjbp3A5+JyDScruYnsH9uKGNMrDZ94Sx7HQ0E\neGp1tyffgH0rWRvUzhwmJUUyWOzTwBjgJfcxVlWf8TowY1Lepi8gtxN07Od3JM3rNhAycxhQu9Lv\nSEyaiXS6r1eAAAAV3ElEQVSw2I3Aax7HYkx62TTfaT4TaXlbP2VmQ8GRDNjyld+RmDQT6TUoY0w8\n1dU616B6Dfc7ksj0KLYalEk4S1DG+GHLcqjdE/wefCE9h9Gpfjsd62zII5M4zSYoEckQkQWJCsb4\n68LJsxqmVDAe2+TeVxT0HnwhbkeJQqtFmQRqNkGpaj3wuYgEdCRLY5LUxs8hMwe6DfI7ksj0dBJU\n0b4vfQ7EpJNImvh6AQtF5F0ReS30aGknEeknItNEZLGILBSR6931XURkqogsd5edW1sIY5LOpi+g\nxxDITJIp2dp2ZoEM4omtQykLyoSKJuVF8uu4K8Zj1wI/U9VPRSQPmCsiU4HvAu+q6j0icjPOiOk3\nxXgOY5KPqtPEd9REvyOJyn11F7CkrhcPvLOcSecmybUzk9QiGSz2fREZAAxU1XdEpB2QGcF+G4GN\n7vMqEVkM9AHOBr7ubvY4MB1LUCadbF8Du7cFfww+1+Db3qamth4YCsCU2WuYMnsNOVkZLJ10hr/B\nmZQWyWCx/w28AEx2V/UBXonmJCJSCIwAZgM93OQVSmIFTexzlYjMEZE55eXl0ZzOmGBb94mz7FPa\nsOrZq8fy7NVjfQqoeTNuPImJJb3JpB6A3Ew4u6Q3M246yefITKqL5BrUtcA4oBJAVZfTRFJpjIh0\nAF4EfqyqlZHup6qPqGqpqpZ279490t2MCb71cyGrrTOdehIoyM8lLyeLOjLIYS81dUpeThYFebl+\nh2ZSXCQJqkZV94ZeiEgWzoy6LRKRbJzk9KSqvuSu3iwivdz3ewFl0YVsTJJbNwd6lzgjNCSJiuoa\nCvJy+Evbh7mk82LKq2v8DsmkgUgS1Psi8gugrYicAjwPvN7STiIiwKPAYlUNnzL+NeBy9/nlwKvR\nhWxMEqvd63Qx7zPK70iiMvnSUoq6tUdzOzNJ/8TkS2zyQuO9SBLUzUA5MB+4GngLuC2C/cYBlwIn\ni8g89/FN4B7gFBFZDpzivjYmPWyeD3U10Le05W0DaFmbo2DPDqhY6ncoJg1E0ouvXkQex+ngoMBS\nVW2xiU9VP8SZnqMx46OK0pgkEhqNo9FOD+vmOMs+yZmglrZxr5utnQ0FR3l2nmY/Q5M2IunFdyaw\nAngA+DPwpYhY31JjYrF6JuT3hU4Bn2KjCZsye0O7rrBmtt+hmDQQyY26fwBOUtUvAUTkcOBN4G0v\nAzMm5ajCmllQdILfkcROBPod69SgjPFYJNegykLJybUS63lnTPS2roTqzdA/yZut+h0LW1dAtd2f\naLzVZA1KRM5zny4UkbeA53CuQV0AfJKA2IxJLatnOssB4/yNo7UGHOcsV38IQ8/1NxaT0pqrQZ3l\nPnKBzcCJOEMUlQM2wGsK2ltbz6INlTYYqFfWzIK2XaD7YL8jaZ3eIyC7PXw1w+9ITIprsgalqlck\nMhDjv/Xbd1NVU2uDgXpBFVbNcGofQZ/ivSWZ2U45VlmCMt5qsZOEiBQBPwIKw7dX1eQaitk0af9g\noA4bDNQDW1c6g8Qed53fkcRH0ddg6lSo2gR5Pf2OxqSoSDpJvAKsAv6E06Mv9DApIjQYaIb7H/vc\n7AwbDDTeVk53loelyGda+DVnuepDf+MwKS2SbuZ7VPUBzyMxvgkNBlqvTutTTW29DQYabyunOfc/\ndT3c70jio9dwyO0EK6ZB8fl+R2NSVCQJ6n4RuQP4N9AwQqSqfupZVCkqyHfHhwYDLcjLoaR/Z8pT\nrKOEr599fR189QEcdVbyX38KyciEw0+GL6c619dSpVwmUH+nIklQxbhj6gGhCxXqvjYpYvKlpQ1f\nzEnnDPM5mhSzfq4zfl2qNO+FDDwFFr7kTF/fa7jf0ZgUFEmCOhc4LHzKDWNMFJa+BRlZcESKDUF5\nxDec5fKplqCMJyLpJPE50MnrQIxJWUv/6XTLbptitw92KIBeJbD8335HYlJUJAmqB7BERP4lIq+F\nHl4HZkyyOuCG560roXwxDP6m32F5Y/AZsPZjqNoc18PaTeMGImviu8PzKExaiuVibCIv4MZ6rgNu\neO75gbNy0OnxDi/hGv0chpwN0++GJa/D6P+K27nspvEDxfpdDFKHh1hEMh/U+4kIxJhk1+gNzxSS\nI0+wtEuRj5F5qPuR0G0QLHo1LgnKbho34SKZD6pKRCrdxx4RqRORykQEZ0wyOeSG5yzh7IwPmTF+\ntb+BeUnEqUWt+hB2VrT6cHbTuAnXYoJS1TxVzXcfucC3cCYuNMaEafSGZ9lNQek5fofmraHngtbD\nghdbfSi7adyEi6STxAFU9RXsHihjGhW64XlorzwuaT+H8twi6NTf77C81WOo08183pNxOdz+zzCf\nS44dQHl1Tcs7mZQUyWCx54W9zABKcW7UTVvJfuHReCd0w/OgvYv4dd19cNb9foeUGCWXwNs3wqYF\n0LN1N3on+qZx+z0HVyQ1qLPCHqcBVcDZXgZlTLI7Zdeb0CYPhqXJOHXFF0BGNnz6hN+RmBQSSS8+\nmxfKmCi0r69izO4PoPQyyOngdziJ0a6Lcy1q3pNw8q2Q29HviEwKaG7K99ub2U9V9dcexJPS9tbW\n82VZNWVVe+yib4Il8rM/befrtGEfjL7S0/MEzthrYP5z8NkUGHut39GYGAXp71RzTXw7G3kAXAnc\n5HFcKSn85kOTWAn77Pfu5Ixdr/Jpzmin80A66T0C+h8H/3kY6vb5HY2JUZD+TjU35XvDpIQikgdc\nD1wBPINNWBgVu/nQPwn/7Oc+Tn79Dl7pcBEj43/04Dv+J/DUBfDZP6D0e4B1QkgWQfw71WwnCRHp\nIiKTgC9wktlIVb1JVcsSEl2KSJabD5+9emzK/RFJ6Ge/ZwfM+D0L2gxnaZs0qz2FDDwF+h4D798L\n+2wcvWQSxL9TTSYoEbkX+ASn116xqt6pqtsSFlkKsZsPGxfLgKDR7tOazz7q+Gb8AXZtZUp+/Mak\nSzoiMP52qNoAM//kdzQpI9bBc6PZL4h/p5qrQf0M6A3cBmwIG+6oyoY6ip7dfHioWNq6Y9kn1s8+\nqnNtXgT/eQiGX8RX2QMjji0lFX0NhpwDH9wLW1b4HU1KiPW6ULT7Be3vlKgG/57b0tJSnTNnjt9h\nNEjXkYXj5eC27pDm2rpj2SdcNJ991Oeq2wd/+wbsWAfXzubCKcsjPlfKqtwIfx4NvYZzUc0tqGSm\n3Kj1iRDr9741v5dEfB4iMldVS1vaLuqhjow52IWTZzV8qSMRS1t3ItvHoz7XO3fCxnlw5h+gfbe4\nx5OU8nvBN38Hqz/k/Or4DIEUNNF+72MR6/c+iNeTYhHJfFDGxFUsbd2JbB+P6lzznoJZf4ZjroKh\nKT4obLRKvgNfzeC8z59mfVZ/IHg1lKCL9XsfxOtJsbAEFYMg3ciWrEJt3QV5OZT070x5BBdxY9nH\n0/gWvAivXgtFJ8Jpv/UslqQ24T6WLJ7Ptdt/D0tLnBl4Aybov+dYv/eJ/L14xZr4YhCkG9mS1eRL\nSynq1p72OVlMOmcYky9tsTk6pn08iU/V6aH24n9B/7Fw8dOQme1ZLEktO5d7O9/B6uwieOYS+ORv\nzucXIEH/Pcf6vU/k78UrVoOKQhBvZDMJtmM9vPlTWPZPZ6K+cx6CNu0P2CSIF9v9tCujA3d1+R1P\n5P0F3vwZrHwfzvidc52qCYn4DO33HHxWg4pCqlx4NDHYthr+dSv8aaTzB/a0u+GCxw9JTqZxNRlt\n4TvPwTfudJL7n0Y5nUt2rI/reaLpuGC/5+DzrAYlIn8HJgBlqjrMXdcFeBYoBFYB306mm39T5cKj\nicDeXRyxdwlH18yFR++AtR87/+jF34aTboHOhX5HmDQarvHs3EvB8T9xap7v/go+uh8+vA/6joYj\nz4TCr0HBEGjTLiFx2e85+Lxs4nsMZ2r48AlibgbeVdV7RORm97X3A8+u/cSZkvoATbSDN9o+vn9d\nRVk1g9tXM6j9Ljp270v55g2wqrLRbQ8+7pCahc7zr2pbFUN8t21i+yYvExz6RsmeRc6TZVURHHf/\nupF7ljhPlu5ocduQUXuWIABLtra4bXgMo/csdV4urnC2rdsL+3bB3l2wbyfs3g47y6F6s3Nz6fY1\n/AalHoH8EfD1m51J+Tr1a+I8pinh13gmnVsMXQ6DCx6DrSudjiaL33BqUwCI836n/pDXCzoUQE4e\ntOngTF2S3Q4yMkEyISPLfZ7RsO7IvUucf/HVkcVWUVbF4PbVDGy/i07d+1FethHW7Gx5R2DwXvf3\nvEb2xx4raX7fgXvd39jazKgOG8t+Le6Tme0MDJwAnt6oKyKFwBthNailwNdVdaOI9AKmq+rglo7T\n6ht1f9PL+WNkTFMy20CHHtC+O3Qpgu5H8ofPhMVthvG3a+x6RCyiull0x3rY8KkzI2/ZIqhcD1Wb\noLoM6m1k9EDJ7ws/XdiqQ0R6o26iO0n0UNWNAG6SKmhqQxG5CrgKoH///q0760VPNVKDopn/tTSy\nPmzbX7/h/A/jlxOGtLhtuLted/a746yhLW4bzXFbvW2T20e27a0vz0cRfntucYSHcFbe8tJ8AO4+\nrzjiGG5+6QsA7vnW8KjivfEFZ7/fnX+0sy4j22lKym7vLLNyDzn/x0u8vQkz1c248SQmvbWYNz7f\nQL0613hOG9qTW8886tCNO/ZxHkeddeh7tXthb7X72AX1taB1UF/n/K7r69zXtUx63flO3XbmkIjj\nnPTmInefRuJqdr/F+/dr1X/0W973t2855/rFN6OLMZb9WtwnK3FNoIHtxaeqjwCPgFODatXBDm/6\nomcsw3osyMlxnhwWXU+jRTluV+Qi74d8SeTwLYsyaviyrJofdxwaVfv9yjbuOF99Ip+Y4qvs3c6T\nXo0lqKatzq52nvRsJIkaT8TtGk9WG8jq4sza24L5od/Y4ZF/7+fn5Ea9D8DczBy+LKvmqoLjoipT\nLL/Nz3PczjhHRBdjLPvFei4vJLoX32a3aQ93adN2pICg30di/BO0wUfjyb733kt0Deo14HLgHnf5\naoLP76tUuz8mme4jSbXPPllMvrS0ocYw6ZxhPkcTH8n0vY9FkH4rntWgRORpYBYwWETWiciVOInp\nFBFZDpzivjZJyu4jMenIvveJ41kNSlUvbuKt8V6d0ySW3Udi0pF97xMnsJ0kTHJozYCUsTQlJLL5\nIUhNHSYyiRr4NdEDscb6XUz277ANdWRaJRUGpDSpI1EdF+x7nxhWgzLGJL1U77iQrqwGZYxJetZx\nITWlfQ0qljbrRLbrxtqmHvRJ2IyJp9Z0XEj0dRr7bUYu7WtQQb/ZLtb4gl4uY+ItWW4Ktt9m5NK2\nBhX0NutY4wt6uYzxStBvCrbfZvTStgYV9DbrWOMLermMSVf224xe2taggn6zXazxBb1cJv0k+704\n8WK/zeilbYKCxN9sF61Y4wt6uYxJV/bbjE5aJ6igt1nHGl/Qy2VMurLfZnTS9hqUMcaYYEvrGpSJ\nD7vGYNKRfe+9ZzUoY4wxgWQJyhhjTCBZE58xJmVYs1tqsRqUMcaYQLIEZYwxJpAsQRljjAkkUVW/\nY2hRaWmpzpkzx+8wjDHGxIGIzFXVFqchthqUMcaYQLIEZYwxJpAsQRljjAkkS1DGGGMCyRKUMcaY\nQLIEZYwxJpAsQRljjAkkS1DGGGMCyRKUMcaYQEqKkSREpBxY3crDdAMq4hBOUKV6+SD1y5jq5QMr\nYyqIR/kGqGr3ljZKigQVDyIyJ5KhNZJVqpcPUr+MqV4+sDKmgkSWz5r4jDHGBJIlKGOMMYGUTgnq\nEb8D8Fiqlw9Sv4ypXj6wMqaChJUvba5BGWOMSS7pVIMyxhiTRCxBGWOMCaSUSFAi8ncRKRORBWHr\nuojIVBFZ7i47u+tFRB4QkS9F5AsRGelf5JETkX4iMk1EFovIQhG53l2fEuUUkVwR+VhEPnfLd5e7\nvkhEZrvle1ZE2rjrc9zXX7rvF/oZfzREJFNEPhORN9zXKVNGEVklIvNFZJ6IzHHXpcR3NEREOonI\nCyKyxP09jk2lMorIYPffL/SoFJEf+1HGlEhQwGPA6Qetuxl4V1UHAu+6rwHOAAa6j6uAhxIUY2vV\nAj9T1aOAMcC1IjKE1ClnDXCyqg4HSoDTRWQM8D/AfW75tgFXuttfCWxT1SOA+9ztksX1wOKw16lW\nxpNUtSTsXplU+Y6G3A/8U1WPBIbj/FumTBlVdan771cCjAJ2AS/jRxlVNSUeQCGwIOz1UqCX+7wX\nsNR9Phm4uLHtkukBvAqckorlBNoBnwLH4tyxnuWuHwv8y33+L2Cs+zzL3U78jj2CsvXF+XGfDLwB\nSCqVEVgFdDtoXcp8R4F84KuD/x1SqYwHletU4CO/ypgqNajG9FDVjQDussBd3wdYG7bdOndd0nCb\nekYAs0mhcrpNX/OAMmAqsALYrqq17ibhZWgon/v+DqBrYiOOyf8CNwL17uuupFYZFfi3iMwVkavc\ndSnzHQUOA8qB/3Obaf8mIu1JrTKGuwh42n2e8DKmcoJqijSyLmn62otIB+BF4MeqWtncpo2sC3Q5\nVbVOnWaFvsAxwFGNbeYuk658IjIBKFPVueGrG9k0acsIjFPVkTjNPteKyAnNbJuM5csCRgIPqeoI\nYCf7m7oak4xlBMC9FjoReL6lTRtZF5cypnKC2iwivQDcZZm7fh3QL2y7vsCGBMcWExHJxklOT6rq\nS+7qlCunqm4HpuNca+skIlnuW+FlaCif+35HYGtiI43aOGCiiKwCnsFp5vtfUqiMqrrBXZbhXLc4\nhtT6jq4D1qnqbPf1CzgJK5XKGHIG8KmqbnZfJ7yMqZygXgMud59fjnPNJrT+MrfnyRhgR6jaGmQi\nIsCjwGJV/WPYWylRThHpLiKd3OdtgW/gXHyeBpzvbnZw+ULlPh94T90G8KBS1VtUta+qFuI0nbyn\nqpeQImUUkfYikhd6jnP9YgEp8h0FUNVNwFoRGeyuGg8sIoXKGOZi9jfvgR9l9PsiXJwu5D0NbAT2\n4WTzK3Ha6t8FlrvLLu62AjyIc31jPlDqd/wRlvF4nGrzF8A89/HNVCkncDTwmVu+BcDt7vrDgI+B\nL3GaGnLc9bnu6y/d9w/zuwxRlvfrwBupVEa3HJ+7j4XAre76lPiOhpWzBJjjfldfATqnYBnbAVuA\njmHrEl5GG+rIGGNMIKVyE58xxpgkZgnKGGNMIFmCMsYYE0iWoIwxxgSSJShjjDGBZAnKGGNMIFmC\nMsYYE0iWoIzxmYiMdufRyXVHY1goIsP8jssYv9mNusYEgIhMwhk5oi3OWG93+xySMb6zBGVMALgj\nR38C7AGOU9U6n0MyxnfWxGdMMHQBOgB5ODUpY9Ke1aCMCQAReQ1nCo4inNlIf+hzSMb4LqvlTYwx\nXhKRy4BaVX1KRDKBmSJysqq+53dsxvjJalDGGGMCya5BGWOMCSRLUMYYYwLJEpQxxphAsgRljDEm\nkCxBGWOMCSRLUMYYYwLJEpQxxphA+v8kw6ZFg8HW6gAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "popt, pcov = curve_fit(fitFunc,cbins,s_data[1],sigma=np.sqrt(s_data[1]),\n", " p0=(A,mu,sigma,B),bounds=([0,300,5,0],[100,500,50,30]))\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": 20, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAHUZJREFUeJzt3Xu8X/Od7/HXWxLiEiJsGrnYoapMW+FsKWNGibbj0lLK\nOUXb1GQa7eEMw2lF5xyJXh4449LO1MM0UyqtulN3oxHCMdqQEJE0HJdGRULikgvGJfE5f6zvrmXb\nl7V3sn6/tX+/9/Px+D1+a33X7fPd2fbH97u+67sUEZiZmVXNRvUOwMzMrDNOUGZmVklOUGZmVklO\nUGZmVklOUGZmVklOUGZmVklOUGYlkXS5pB92se0bkh4o4ZqzJP3dhj5v2cr6eVj/5gRllSFpsaT/\nlLRG0kpJD0r6lqRCv6eSWiWFpIFlx2pm5XOCsqr5YkQMAXYEzgXOAC6tb0jNx0neqsAJyiopIlZF\nxC3AfwMmSPoEgKTDJD0qabWk5yVNzR12f/peKel1SftK2lnSPZJekfSypF9LGtrVdSX9JJ13taS5\nkv46t22qpGsl/TK18hZKastt31PSI2nbNcDgHqopSf8iaZWkJyQdlNtwgqRF6VzPSjqxw4FHSJqX\n4nxG0sGdnHy4pPmS/mdaHyPp/nTOuyVdLOmKtK299TlR0p+Ae1L54ameK1P34W6584ekj+bW/9yl\nKekASUsknS5puaRlkk7I7buNpFtS/A8BO/fws7Im5ARllRYRDwFLgPZE8QbwdWAocBjwbUlfStv2\nT99DI2KLiPgdIOAcYAdgN2AUMLWbSz4MjAWGAVcC10nKJ5rDgavT9W8BfgogaWPgJuBX6djrgC/3\nUL1PA88C2wJTgBslDUvblgNfALYETgAukrRXutY44JfAd1Ic+wOL8yeW1ArcB/w0Is5PxVcCDwHb\npJ/B1zqJ6TNkP6e/kfQx4CrgVKAFuAO4NdW1iI8AWwEjgInAxZK2TtsuBt4ChgN/mz5mH+AEZf3B\nUrI/+kTErIh4PCLei4j5ZH9AP9PVgRHxdETMiIi3I2IFcGEP+18REa9ExNqIuADYBNg1t8sDEXFH\nRKwjS0Z7pPJ9gEHAjyPi3Yi4nizZdWd5bv9rgCfJki4RcXtEPBOZ+4Df8n6Snghclur1XkS8EBFP\n5M67OzALmBIR0wAkjQb2Bs6KiHci4gGyBNvR1Ih4IyL+k6z1enu6zrvA+cCmwF/2UK927wLfT/W7\nA3gd2FXSALLkfVa61gJgesFzWhNxgrL+YATwKoCkT0u6V9IKSauAb5G1QDolaTtJV0t6QdJq4Ioe\n9j89da2tkrSSrAWQ3//F3PKbwOB0v2YH4IX44OzLz/VQr8723yHFcYik30t6NcVxaC6OUcAz3Zz3\neOAF4Ppc2Q7AqxHxZq7s+U6OzZftkK9DRLyXto/otlbveyUi1ubW3wS2IGuNDexwrZ5+VtaEnKCs\n0iTtTfYHsX0I8pVk/+c/KiK2Av6VrBsPoLOp+c9J5Z+KiC2Br+b273itvyYblPFfga0jYiiwqqv9\nO1gGjJCU33d0D8d0tv9SSZsAN5C1WLZPcdyRi+N5ur9nMxV4GbgytVba4xsmabPcfqM6OTb/M1xK\nNlgFyG6YpWNeSEVvAvnzfaSbmPJWAGs7XL+nn5U1IScoqyRJW0r6Atn9nisi4vG0aQhZS+CtdC/m\nuNxhK4D3gJ1yZUPIupZWShpBdt+mK0PI/nCuAAZKOovsHlARv0vH/r2kgZKOAsb1cMx2af9Bko4h\nu/dzB7AxWdfiCmCtpEOAz+eOuxQ4QdJBkjaSNELSx3Pb3wWOATYHfiVpo4h4DpgDTJW0saR9gS/2\nEN+1wGHpOoOA04G3gQfT9nnAcZIGpEEaXXad5qXu0RtTLJtJ2h2YUORYay5OUFY1t0paQ9ZK+Eey\ne0Yn5Lb/d+D7aZ+zyP6IApC6r34E/EcadbYPcDawF1lL6HayP4xduQu4E/h/ZF1Ob9F5N9iHRMQ7\nwFHAN4DXyO7fdHctgNnALmStnR8BR6f7X2uAv091e40sCf/5flEaOHICcFGq133kWjod4tkOuEzZ\ns2THA/sCrwA/BK4hSzhd1elJshbnv6QYv0j2GMA7aZdTUtnKdO6beqhv3slk3X0vApcDv+jFsdYk\n5BcWmjWnNBT+iYiYUu9YzDrjFpRZk5C0t7LnwjZKXXJH0LtWj1lN+Wlxs+bxEbJux23Ini37dkQ8\nWt+QzLrmLj4zM6skd/GZmVkl9Ysuvm233TZaW1vrHYaZmW0Ac+fOfTkiWnrar18kqNbWVubMmVPv\nMMzMbAOQVGjmEHfxmZlZJTlBmZlZJTlBmZlZJTlBmZlZJTlBmZlZJTlBmZlZJZWeoNJU/I9Kui2t\nj5E0W9JTkq7pxeujzcysidSiBXUKsCi3fh5wUUTsQvYqgYk1iMHMzPqZUhOUpJHAYcDP07qA8bz/\nKurpwJfKjMHMzPqnsmeS+DHwXbI3lUI2i/LKiFib1peQvc77QyRNAiYBjB7tt0Gb9cnUrfpwzKoN\nH4dZH5TWgkqv614eEXPzxZ3s2ul06hExLSLaIqKtpaXHKZvMzKzBlNmC2g84XNKhwGBgS7IW1VBJ\nA1MraiSwtMQYzMysnyqtBRURZ0bEyIhoBb4C3BMRxwP3Aken3SYAN5cVg5mZ9V/1eA7qDOA0SU+T\n3ZO6tA4xmJlZxdXkdRsRMQuYlZafBcbV4rpmZtZ/eSYJMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCco\nMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOrJCcoMzOr\nJCcoMzOrpJq8D8rM1tPUreodgVnNuQVlZmaVVFqCkjRY0kOSHpO0UNLZqfxySX+UNC99xpYVg5mZ\n9V9ldvG9DYyPiNclDQIekHRn2vadiLi+xGubmVk/V1qCiogAXk+rg9InyrqemZk1llLvQUkaIGke\nsByYERGz06YfSZov6SJJm3Rx7CRJcyTNWbFiRZlhmplZBZWaoCJiXUSMBUYC4yR9AjgT+DiwNzAM\nOKOLY6dFRFtEtLW0tJQZppmZVVBNRvFFxEpgFnBwRCyLzNvAL4BxtYjBzMz6lzJH8bVIGpqWNwU+\nCzwhaXgqE/AlYEFZMZiZWf9V5ii+4cB0SQPIEuG1EXGbpHsktQAC5gHfKjEGMzPrp8ocxTcf2LOT\n8vFlXdPMzBqHZ5IwM7NKcoIyM7NKcoIyM7NKcoIyM7NKcoIyM7NKcoIyM7NK6nGYuaT9gKnAjml/\nkc0Fu1O5oZmZWTMr8hzUpcA/AHOBdeWGY2ZmlimSoFZFxJ0972ZmZrbhFElQ90r6J+BGspcQAhAR\nj5QWlZmZNb0iCerT6bstVxaApywyM7PS9JigIuLAWgRiZmaW12WCkvTViLhC0mmdbY+IC8sLy8zM\nml13LajN0/eQWgRiZmaW12WCioifpe+zaxeOmZlZpseZJCTtJOlWSSskLZd0syQ/pGtmZqUqMtXR\nlcC1ZG/I3QG4DriqzKDMzMyKJChFxK8iYm36XEE2zLz7g6TBkh6S9JikhZLOTuVjJM2W9JSkayRt\nvL6VMDOzxtNlgpI0TNIwsgd1J0tqlbSjpO8Ctxc499vA+IjYAxgLHCxpH+A84KKI2AV4DZi4/tUw\nM7NG090ovrlkLSWl9RNz2wL4QXcnjogAXk+rg9Kn/QHf41L5dLKJaC/pTdBmZtb4uhvFN2Z9Ty5p\nAFmi+yhwMfAMsDIi1qZdlgAjujh2EjAJYPTo0esbipmZ9TOlvg8qItZFxFhgJDAO2K2z3bo4dlpE\ntEVEW0tLS5lhmplZBdXkhYURsRKYBewDDJXU3nIbCSytRQxmZta/dJuglBnVlxNLapE0NC1vCnwW\nWATcCxyddpsA3NyX85uZWWPrNkGlgQ439fHcw8lGAM4HHgZmRMRtwBnAaZKeBrYheyGimZnZBxR5\n3cbvJe0dEQ/35sQRMR/Ys5PyZ8nuR5mZmXWpSII6EPiWpMXAG2TDziMiPlVmYGZm1tyKJKhDSo/C\nzMysgx5H8UXEc8AoslkhngPeLHKcmZnZ+igym/kUsoENZ6aiQcAVZQZlZmZWpCV0JHA42f0nImIp\nfomhmZmVrEiCeicNNw8ASZv3sL+Zmdl6K5KgrpX0M7IZIL4J3A38W7lhmZlZs+txFF9EnC/pc8Bq\n4GPAWRExo/TIzMysqRUZZg7wOLApWTff4+WFY2Zmlikyiu/vgIeAo8jm0Pu9pL8tOzAzM2tuRVpQ\n3wH2jIhXACRtAzwIXFZmYGZm1tyKDJJYAqzJra8Bni8nHDMzs0yRFtQLwGxJN5PdgzoCeEjSaQAR\ncWGJ8ZmZWZMqkqCeSZ927e9v8sO6ZmZWmiLDzM+uRSBmZmZ5nvTVzMwqyQnKzMwqqbQEJWmUpHsl\nLZK0UNIpqXyqpBckzUufQ8uKwczM+q8iD+r+H0lbShokaaaklyV9tcC51wKnR8RuwD7ASZJ2T9su\nioix6XPHesRvZmYNqkgL6vMRsRr4AtkzUR8je3i3WxGxLCIeSctrgEXAiPWI1czMmkiRBDUofR8K\nXBURr/b2IpJagT2B2anoZEnzJV0maesujpkkaY6kOStWrOjtJc3MrJ8rkqBulfQE0AbMlNQCvFX0\nApK2AG4ATk0tsUuAnYGxwDLggs6Oi4hpEdEWEW0tLS1FL2dmZg2iSIKaAuwLtEXEu8CbZG/Y7ZGk\nQWTJ6dcRcSNARLwUEesi4j2y90qN61PkZmbW0IokqN9FxGsRsQ4gIt4A7uzpIEkCLgUW5adDkjQ8\nt9uRwILehWxmZs2gy5kkJH2EbFDDppL2BJQ2bQlsVuDc+wFfAx6XNC+VfQ84VtJYsnn9FgMn9i10\nMzNrZN1NdfQ3wDeAkUB+Qtg1ZImmWxHxAO8ntbyGH1beOvn2Ph23+NzDNnAkZmb9V5cJKiKmA9Ml\nfTkibqhhTGZmZoVmM79N0nFAa37/iPh+WUGZmZkVSVA3A6uAucDb5YZjZmaWKZKgRkbEwaVHYmZm\nllNkmPmDkj5ZeiRmZmY5RVpQfwV8Q9Ifybr4BEREfKrUyMzMrKkVSVCHlB6FmZlZBz128UXEc8Ao\nYHxafrPIcWZmZuujyPugpgBnAGemokHAFWUGZWZmVqQldCTZ5LBvAETEUmBImUGZmZkVuQf1TkSE\npACQtHnJMZWiL9MPeeohM7P6KdKCulbSz4Chkr4J3E32mgwzM7PS9NiCiojzJX0OWA3sCpwVETNK\nj8zMzJpajwlK0j8A1zkpmZlZLRXp4tsSuEvS/5V0kqTtyw7KzMysyHNQZ0fEXwAnATsA90m6u/TI\nzMysqfXmgdvlwIvAK8B25YRjZmaWKfKg7rclzQJmAtsC3ywyD5+kUZLulbRI0kJJp6TyYZJmSHoq\nfW+9vpUwM7PGU6QFtSNwakT8RURMiYg/FDz3WuD0iNgN2Ac4SdLuwGRgZkTsQpb0JvclcDMza2xF\n7kFNBraQdAKApBZJYwoctywiHknLa4BFwAjgCGB62m068KU+xm5mZg2sJnPxSWoF9gRmA9tHxDLI\nkhhd3M+SNEnSHElzVqxY0ZvLmZlZAyh9Lj5JWwA3kHUTri56XERMi4i2iGhraWkpepiZmTWIIgnq\nnYgIoNdz8UkaRJacfh0RN6bilyQNT9uHk40ONDMz+4DS5uKTJOBSYFFEXJjbdAswIS1PAG7uXchm\nZtYMypyLbz/ga8Djkualsu8B55IlvYnAn4Bj+hS5mZk1tCKv2yAlpF7NxRcRDwDqYvNBvTmXmZk1\nH7+63czMKskJyszMKqnLBCVpZvo+r3bhmJmZZbq7BzVc0meAwyVdTYf7Se2zRJiZmZWhuwR1Ftk8\neSOBCztsC2B8WUGZmZl1maAi4nrgekn/OyJ+UMOYzMzMCj0H9QNJhwP7p6JZEXFbuWGZmVmzKzJZ\n7DnAKcAf0ueUVGZmZlaaIg/qHgaMjYj3ACRNBx7l/dnNzczMNriiz0ENzS1vVUYgZmZmeUVaUOcA\nj0q6l2yo+f40SeupdfLt9Q7BzKxpFRkkcZWkWcDeZAnqjIh4sezAzMysuRWdLHYZ2WsyzMzMasJz\n8ZmZWSU5QZmZWSV128UnaSNgfkR8okbxmDW+qR4Ia1ZEty2o9OzTY5JG1ygeMzMzoFgX33BgoaSZ\nkm5p//R0kKTLJC2XtCBXNlXSC5Lmpc+h6xO8mZk1riKj+M7u47kvB34K/LJD+UURcX4fz2lmZk2i\nyHNQ90naEdglIu6WtBkwoMBx90tqXf8QzcysGfWYoCR9E5gEDAN2BkYA/woc1Mdrnizp68Ac4PSI\neK2L605K12X06Oa4BdaXmSsWn3tYCZFYU+vrII6pqzZsHNb0ityDOgnYD1gNEBFPAdv18XqXkCW5\nscAy4IKudoyIaRHRFhFtLS0tfbycmZn1V0US1NsR8U77iqSBZG/U7bWIeCki1qXRgf8GjOvLeczM\nrPEVSVD3SfoesKmkzwHXAbf25WKShudWjwQWdLWvmZk1tyKj+CYDE4HHgROBO4Cf93SQpKuAA4Bt\nJS0BpgAHSBpL1gJbnM5nZmb2IUVG8b2XXlI4myyxPBkRPXbxRcSxnRRf2vsQzcysGRUZxXcY2ai9\nZ8hetzFG0okRcWfZwZmZWfMq0sV3AXBgRDwNIGln4HbACcrMzEpTZJDE8vbklDwLLC8pHjMzM6Cb\nFpSko9LiQkl3ANeS3YM6Bni4BrGZmVkT666L74u55ZeAz6TlFcDWpUVkZmZGNwkqIk6oZSBmZmZ5\nRUbxjQH+B9Ca3z8iDi8vLDMza3ZFRvHdRPb80q3Ae+WGY2ZmlimSoN6KiH8uPRIzM7OcIgnqJ5Km\nAL8F3m4vjIhHSovKzMyaXpEE9Unga8B43u/ii7RuddaXd0iB3yNlZtVXJEEdCeyUf+WGmZlZ2YrM\nJPEYMLTsQMzMzPKKtKC2B56Q9DAfvAflYeZmZlaaIglqSulRmJmZdVDkfVD31SIQMzOzvCIzSawh\nG7UHsDEwCHgjIrYsMzAzM2tuPQ6SiIghEbFl+gwGvgz8tKfjJF0mabmkBbmyYZJmSHoqfXvSWTMz\n61SRUXwfEBE3UewZqMuBgzuUTQZmRsQuwMy0bmZm9iFFuviOyq1uBLTxfpdflyLifkmtHYqPAA5I\ny9OBWcAZPYdpZmbNpsgovvx7odYCi8kSTV9sHxHLACJimaTt+ngeMzNrcEVG8dXlvVCSJgGTAEaP\nHl2PEBpaX6ZI8vRIZlZL3b3y/axujouI+EEfrveSpOGp9TQcWN7NBaYB0wDa2tp67FI0M7PG0t0g\niTc6+QBMpO/3jW4BJqTlCcDNfTyPmZk1uO5e+X5B+7KkIcApwAnA1cAFXR2XO+YqsgER20paQjYj\nxbnAtZImAn8Cjlmf4M3MrHF1ew9K0jDgNOB4slF3e0XEa0VOHBHHdrHpoF5FaGZmTam7e1D/BBxF\ndh/okxHxes2iskryu6c6mLpVvSMwa2jd3YM6HdgB+F/AUkmr02eNpNW1Cc/MzJpVd/egej3LhJmZ\n2YbiJGRmZpXkBGVmZpVUZKojs5rzgAwzcwvKzMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnKzMwqyQnK\nzMwqyQnKzMwqyQnKzMwqyQnKzMwqyTNJWOn6OiuE9TN9ff3I1FUbNg5rGG5BmZlZJdWlBSVpMbAG\nWAesjYi2esRhZmbVVc8uvgMj4uU6Xt/MzCrMXXxmZlZJ9UpQAfxW0lxJk+oUg5mZVVi9uvj2i4il\nkrYDZkh6IiLuz++QEtckgNGjR9cjRuuHajti8Mo+HbV48HEbOI4Nq/WtxqyX9T91aUFFxNL0vRz4\nDTCuk32mRURbRLS1tLTUOkQzM6uzmicoSZtLGtK+DHweWFDrOMzMrNrq0cW3PfAbSe3XvzIi/r0O\ncZiZWYXVPEFFxLPAHrW+rpmZ9S+e6sisxvoyCMEDEKwZ+TkoMzOrJCcoMzOrJCcoMzOrJCcoMzOr\nJA+SMOsH+jq7Q39Qy9k/Fp97WM2u1Zd61TK+/sAtKDMzqyQnKDMzqyQnKDMzqyQnKDMzqyQPkjCz\nDaKRB3JYfbgFZWZmleQEZWZmleQEZWZmleQEZWZmleQEZWZmleRRfGbWNGo5rVIt9bVeVZ9ayS0o\nMzOrpLokKEkHS3pS0tOSJtcjBjMzq7aaJyhJA4CLgUOA3YFjJe1e6zjMzKza6tGCGgc8HRHPRsQ7\nwNXAEXWIw8zMKqwegyRGAM/n1pcAn+64k6RJwKS0+rqkJ9fzutsCL6/nOaqs0esHrmMjaPT6wXrU\nUedt4EjKud6G+DfcschO9UhQ6qQsPlQQMQ2YtsEuKs2JiLYNdb6qafT6gevYCBq9ftD4daxl/erR\nxbcEGJVbHwksrUMcZmZWYfVIUA8Du0gaI2lj4CvALXWIw8zMKqzmXXwRsVbSycBdwADgsohYWINL\nb7Duwopq9PqB69gIGr1+0Ph1rFn9FPGh2z9mZmZ155kkzMyskpygzMyskhomQUm6TNJySQtyZcMk\nzZD0VPreOpVL0j+nqZbmS9qrfpEXI2mUpHslLZK0UNIpqbwh6ihpsKSHJD2W6nd2Kh8jaXaq3zVp\nYA2SNknrT6ftrfWMvzckDZD0qKTb0nrD1FHSYkmPS5onaU4qa4jf0XaShkq6XtIT6b/HfRupjpJ2\nTf9+7Z/Vkk6tRx0bJkEBlwMHdyibDMyMiF2AmWkdsmmWdkmfScAlNYpxfawFTo+I3YB9gJOUTRHV\nKHV8GxgfEXsAY4GDJe0DnAdclOr3GjAx7T8ReC0iPgpclPbrL04BFuXWG62OB0bE2NyzMo3yO9ru\nJ8C/R8THgT3I/i0bpo4R8WT69xsL/BfgTeA31KOOEdEwH6AVWJBbfxIYnpaHA0+m5Z8Bx3a2X3/5\nADcDn2vEOgKbAY+QzTDyMjAwle8L3JWW7wL2TcsD036qd+wF6jaS7D/u8cBtZA+uN0wdgcXAth3K\nGuZ3FNgS+GPHf4dGqmOHen0e+I961bGRWlCd2T4ilgGk7+1SeWfTLY2ocWx9lrp69gRm00B1TF1f\n84DlwAzgGWBlRKxNu+Tr8Of6pe2rgG1qG3Gf/Bj4LvBeWt+GxqpjAL+VNFfZdGXQQL+jwE7ACuAX\nqZv255I2p7HqmPcV4Kq0XPM6NnqC6kqh6ZaqSNIWwA3AqRGxurtdOymrdB0jYl1k3QojySYV3q2z\n3dJ3v6ufpC8AyyNibr64k137bR2B/SJiL7Jun5Mk7d/Nvv2xfgOBvYBLImJP4A3e7+rqTH+sIwDp\nXujhwHU97dpJ2QapY6MnqJckDQdI38tTeb+cbknSILLk9OuIuDEVN1QdASJiJTCL7F7bUEntD5Tn\n6/Dn+qXtWwGv1jbSXtsPOFzSYrJZ/MeTtagapo4RsTR9Lye7bzGOxvodXQIsiYjZaf16soTVSHVs\ndwjwSES8lNZrXsdGT1C3ABPS8gSy+zbt5V9Po0/2AVa1N12rSpKAS4FFEXFhblND1FFSi6ShaXlT\n4LNkN5/vBY5Ou3WsX3u9jwbuidQBXlURcWZEjIyIVrKuk3si4ngapI6SNpc0pH2Z7P7FAhrkdxQg\nIl4Enpe0ayo6CPgDDVTHnGN5v3sP6lHHet+E24A3864ClgHvkmX0iWT99TOBp9L3sLSvyF6a+Azw\nONBW7/gL1O+vyJrN84F56XNoo9QR+BTwaKrfAuCsVL4T8BDwNFlXwyapfHBafzpt36nedehlfQ8A\nbmukOqZ6PJY+C4F/TOUN8Tuaq+dYYE76Xb0J2LoB67gZ8AqwVa6s5nX0VEdmZlZJjd7FZ2Zm/ZQT\nlJmZVZITlJmZVZITlJmZVZITlJmZVZITlJmZVZITlJmZVZITlFmdSdo7vUdncJqNYaGkT9Q7LrN6\n84O6ZhUg6YdkM0dsSjbX2zl1Dsms7pygzCogzRz9MPAW8JcRsa7OIZnVnbv4zKphGLAFMISsJWXW\n9NyCMqsASbeQvYJjDNnbSE+uc0hmdTew513MrEySvg6sjYgrJQ0AHpQ0PiLuqXdsZvXkFpSZmVWS\n70GZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVklOUGZmVkl/X9XigLjX1g12gAAAABJ\nRU5ErkJggg==\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 expon # 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 = expon.rvs(scale=200,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": 21, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xlc1VX6wPHPc9lVQMQVQUHFnUTFhSxNLdM0bd+ccprW\nqZlqprKabP1Z2bROy5RNTtlkZYuapWZlWqamaZn7hiJuKaJsKus9vz++F0RkueDdgOf9et0XcO/3\ne77PReRwvuc5zxFjDEoppZSvsXk7AKWUUqoi2kEppZTySdpBKaWU8knaQSmllPJJ2kEppZTySdpB\nKaWU8knaQSnlJiLyrohMruS1P4rIj2645hIRudnV7bqbu74fqm7TDkr5DBFJFZETIpIjIpkislxE\nbhcRp35ORSRWRIyI+Ls7VqWU+2kHpXzNxcaYUKA9MAV4AJjm3ZAaHu3klS/QDkr5JGNMljFmLnA1\nMEFEegKIyGgR+VVEskVkj4g8Xua0HxwfM0UkV0SSRaSjiHwnIhkiclhEZohI08quKyL/crSbLSJr\nROTcMq89LiIfi8h7jlHeRhFJKvN6bxH5xfHaTCC4mrcpIvKqiGSJyBYRGV7mhRtFZLOjrZ0iclu5\nE8eJyFpHnCkiMrKCxtuIyDoRuc/xdZyI/OBo81sReV1E3ne8VjL6vElE0oDvHM+PdbzPTMftw25l\n2jci0qnM16W3NEXkPBHZKyL3isghETkgIjeWOTZSROY64l8FdKzme6UaIO2glE8zxqwC9gIlHcUx\n4AagKTAa+LOIXOJ4bbDjY1NjTBNjzApAgGeAKKAbEAM8XsUlfwYSgWbAB8AnIlK2oxkLfOS4/lzg\nNQARCQTmAP9znPsJcHk1b28AsBNoDjwGzBKRZo7XDgFjgDDgRuAlEenjuFZ/4D3gfkccg4HUsg2L\nSCzwPfCaMeZ5x9MfAKuASMf34PoKYhqC9X26UEQ6Ax8C9wAtgPnAF4736ozWQDjQFrgJeF1EIhyv\nvQ7kAW2APzkeSp1COyhVF+zH+qWPMWaJMWa9McZujFmH9Qt0SGUnGmN2GGO+McbkG2PSgRerOf59\nY0yGMabIGPMCEAR0KXPIj8aY+caYYqzOqJfj+YFAAPCyMabQGPMpVmdXlUNljp8JbMXqdDHGzDPG\npBjL98DXnOykbwL+63hfdmPMPmPMljLtdgeWAI8ZY94CEJF2QD/gUWNMgTHmR6wOtrzHjTHHjDEn\nsEav8xzXKQSeB0KAs6t5XyUKgScd728+kAt0ERE/rM77Uce1NgDTnWxTNSDaQam6oC1wBEBEBojI\nYhFJF5Es4HasEUiFRKSliHwkIvtEJBt4v5rj73XcWssSkUysEUDZ438v8/lxINgxXxMF7DOnVl/e\nXc37quj4KEcco0TkJxE54ojjojJxxAApVbQ7HtgHfFrmuSjgiDHmeJnn9lRwbtnnosq+B2OM3fF6\n2yrf1UkZxpiiMl8fB5pgjcb8y12ruu+VaoC0g1I+TUT6Yf1CLElB/gDrL/8YY0w48CbWbTyAikrz\nP+N4/ixjTBjwhzLHl7/WuVhJGVcBEcaYpkBWZceXcwBoKyJlj21XzTkVHb9fRIKAz7BGLK0cccwv\nE8ceqp6zeRw4DHzgGK2UxNdMRBqVOS6mgnPLfg/3YyWrANaEmeOcfY6njgNl22tdRUxlpQNF5a5f\n3fdKNUDaQSmfJCJhIjIGa77nfWPMesdLoVgjgTzHXMx1ZU5LB+xAhzLPhWLdWsoUkbZY8zaVCcX6\nxZkO+IvIo1hzQM5Y4Tj3LhHxF5HLgP7VnNPScXyAiFyJNfczHwjEurWYDhSJyChgRJnzpgE3ishw\nEbGJSFsR6Vrm9ULgSqAx8D8RsRljdgOrgcdFJFBEkoGLq4nvY2C04zoBwL1APrDc8fpa4DoR8XMk\naVR667Qsx+3RWY5YGolId2CCM+eqhkU7KOVrvhCRHKxRwsNYc0Y3lnn9DuBJxzGPYv0SBcBx++op\nYJkj62wg8ATQB2skNA/rF2NlFgILgG1Yt5zyqPg22GmMMQXAZcAfgaNY8zdVXQtgJRCPNdp5CrjC\nMf+VA9zleG9HsTrh0vkiR+LIjcBLjvf1PWVGOuXiaQn8V6y1ZOOBZCADmAzMxOpwKntPW7FGnK86\nYrwYaxlAgeOQux3PZTranlPN+y3rL1i3+34H3gXeqcG5qoEQ3bBQqYbJkQq/xRjzmLdjUaoiOoJS\nqoEQkX5irQuzOW7JjaNmox6lPEpXiyvVcLTGuu0YibW27M/GmF+9G5JSldNbfEoppXyS3uJTSinl\nk+rELb7mzZub2NhYb4ehlFLKBdasWXPYGNOiuuPqRAcVGxvL6tWrvR2GUkopFxARpyqH6C0+pZRS\nPkk7KKWUUj5JOyillFI+qU7MQSmlGrbCwkL27t1LXl6et0NRNRAcHEx0dDQBAQG1Ol87KKWUz9u7\ndy+hoaHExsZyagF45auMMWRkZLB3717i4uJq1Ybe4lNK+by8vDwiIyO1c6pDRITIyMgzGvVqB6WU\nqhO0c6p7zvTfTDsopZRSPkk7KKUaoKunruDqqSu8HUad0qRJk2qPufnmm9m0aRMATz/99CmvnX32\n2bW+hp+fH4mJiaWP1NTU6gN2UmZmJv/+979Lv96/fz9XXHGFy9o/E3WiWGxSUpLRShJKuU5J5zTz\ntmQvR+KczZs3061bN6/G0KRJE3Jzc912fFXn1KYtZ6WmpjJmzBg2bNjglvYr+rcTkTXGmKTqztUR\nlFJK1cCSJUs477zzuOKKK+jatSvjx4+n5A/98847j9WrV/Pggw9y4sQJEhMTGT9+PHBydJSbm8vw\n4cPp06cPCQkJfP7557WK49133+Uvf/lL6ddjxoxhyZIlpdd6+OGH6dWrFwMHDuTgwYMAHDx4kEsv\nvZRevXrRq1cvli9fzoMPPkhKSgqJiYncf//9pKam0rNnT8BKTrnxxhtJSEigd+/eLF68uPTal112\nGSNHjiQ+Pp6JEyfW6j1UR9PMlVJ1y4IH4ff1rm2zdQKMmuL04b/++isbN24kKiqKQYMGsWzZMs45\n55zS16dMmcJrr73G2rVrTzs3ODiY2bNnExYWxuHDhxk4cCBjx46tMqGgpLMDiIuLY/bs2VXGd+zY\nMQYOHMhTTz3FxIkT+c9//sOkSZO46667GDJkCLNnz6a4uJjc3FymTJnChg0bSmMte/vw9ddfB2D9\n+vVs2bKFESNGsG3bNgDWrl3Lr7/+SlBQEF26dOGvf/0rMTExzn0DnaQdlFJK1VD//v2Jjo4GKJ0T\nKttBVcUYwz/+8Q9++OEHbDYb+/bt4+DBg7Ru3brSc0JCQirs7CoTGBjImDFjAOjbty/ffPMNAN99\n9x3vvfceYM1rhYeHc/To0Urb+fHHH/nrX/8KQNeuXWnfvn1pBzV8+HDCw8MB6N69O7t379YOSinV\nwNVgpOMuQUFBpZ/7+flRVFTk9LkzZswgPT2dNWvWEBAQQGxsbK3WCvn7+2O320u/LttGQEBA6Yis\npvGVVVWOwpl8D5ylc1BKKeUGAQEBFBYWnvZ8VlYWLVu2JCAggMWLF7N7t1M7T5wmNjaWtWvXYrfb\n2bNnD6tWrar2nOHDh/PGG28AUFxcTHZ2NqGhoeTk5FR4/ODBg5kxYwYA27ZtIy0tjS5dutQq3trQ\nDkoppdzg1ltv5ayzzipNkigxfvx4Vq9eTVJSEjNmzKBr1661an/QoEHExcWRkJDAfffdR58+fao9\n51//+heLFy8mISGBvn37snHjRiIjIxk0aBA9e/bk/vvvP+X4O+64g+LiYhISErj66qt59913Txk5\nuZummSvVAGmaufIUTTNXSilV72gHpZRSyidpB6WUUsonaQellFLKJ2kHpZRSyidpB6WUUsonaQel\nlFJOqGjLi9WrV3PXXXcBVhHZ5cuXlx4/Z86c0q03akK33DhJSx0ppZQTKqqHFxsbS1KStZxnyZIl\nNGnSpHTfpzlz5jBmzBi6d+/utuu7SkkHdccddwAQFRXFp59+6pZr1YSOoJRSqpaWLFnCmDFjSE1N\n5c033+Sll14iMTGR77//nrlz53L//feTmJhISkoKKSkpjBw5kr59+3LuueeyZcsWAHbt2kVycjL9\n+vXjkUceqdH16/uWGzqCUkrVKc+uepYtR7a4tM2uzbryQP8Hqjymqi0vYmNjuf3222nSpAn33Xcf\nAGPHjmXMmDGlt8qGDx/Om2++SXx8PCtXruSOO+7gu+++4+677+bPf/4zN9xwQ+n2FjW9fkXqw5Yb\n2kEpVYfVtZJFddmZ3GLLzc1l+fLlXHnllaXP5efnA7Bs2TI+++wzAK6//noeeKDijrIhbrmhHZRS\nqk6pbqTji+x2O02bNq20g6lqs8Kq1PctN3QOSimlXKD8thVlvw4LCyMuLo5PPvkEsH7p//bbb4BV\nlfyjjz4CKN3awln1fcsN7aCUUsoFLr74YmbPnk1iYiJLly7lmmuu4bnnnqN3796kpKQwY8YMpk2b\nRq9evejRoweff/45YG2B8frrr9OvXz+ysrJqdM36vuWGbrehVB1W2zmoujZ3pdtt1F0+vd2GiPiJ\nyK8i8qXj6zgRWSki20VkpogEujsGpdSpCorsbNqfzaGcmm81rpSneOIW393A5jJfPwu8ZIyJB44C\nN3kgBqUaFGMMu7J2sfbQWorsp09c78s8QU5+Ea98u90L0SnlHLdm8YlINDAaeAr4u1gpJcOA6xyH\nTAceB95wZxxKNSTF9mIeWvoQC1IXAHB21Nk8N+Q5wgLD6DJpAflFJ7O+3l+Zxvsr0wjyt7F18ihv\nhaxUhdw9gnoZmAiU/I+IBDKNMSV/0u0F2lZ0oojcKiKrRWR1enq6m8NUqv54bvVzLEhdwE09b+L+\npPtZdWAV9y65F2MMSycOZWxiFDZHVnNwgI1xiVEsfWCod4NWqgJuG0GJyBjgkDFmjYicV/J0BYdW\nmKVhjHkLeAusJAm3BKlUPbPt6DZmbJ7BtV2v5Z6+9wDgZ/NjyqopLN23lMHRgwkN8sduQATyi+yE\nBvnTMjTYy5ErdTp3jqAGAWNFJBX4COvW3stAUxEp6Rijgf1ujEGpBmXqb1NpHNCYOxPvLH3uqi5X\n0T6sPS+sfoFiezGHc/NpGRpEjzZhjB/QnvTcfC9G7D5XT11Rmq2o6ia3dVDGmIeMMdHGmFjgGuA7\nY8x4YDFQUsd9AvC5u2JQqiHZlbWLb3Z/w3VdryM8KLz0+QBbAHf0uoOdWTtZ9fsqpl6fRFzzxjQO\n8mfyJT2Zen212b6Kk9td9OrViz59+pRureGprSliY2M5fPhwhc8nJCSQkJBA9+7dmTRpUmkZpcqU\n317DV3ljoe4DWAkTO7DmpKZ5IQal6p0vUr5ARLiu23WnvTa8/XBCA0OZmzLXC5F5h6tT6Utq4f32\n228888wzPPTQQ4BvbE2xePFi1q9fz6pVq9i5cye33nprlcdrB1WGMWaJMWaM4/Odxpj+xphOxpgr\njTH18/6CUh5kjOGr1K/o37o/zUOan/Z6kF8QI2NHsihtEccKj3khQs9zZyp9dnY2ERERAKdsTVHV\nFhSVbX+Rnp7O5ZdfTr9+/ejXrx/Lli0DICMjgxEjRtC7d29uu+22Kmvilb3Gm2++yZw5czhy5Ai5\nubkMHz6cPn36kJCQUFq9ovz2GpUd521aLFapemBTxib25OzhloRbKj1mbMexfLLtExalLQJaeC44\nD3NXKn3Jdhd5eXkcOHCA7777rsLjKtuCorLtL+6++27+9re/cc4555CWlsaFF17I5s2beeKJJzjn\nnHN49NFHmTdvHm+99ZZTcZbU/du+fTt9+/Zl9uzZhIWFcfjwYQYOHMjYsWNP216jqKiowuNqW8TW\nVbSDUqoe+Cr1K/xt/gxrN6zSY3q16EWLkBYs3bsUuMxzwXnY0olDmTx/M1/+th+7sVLpL+zRmodH\nn1mppLLbXaxYsYIbbriBDRs2nHZcZVtQVLb9xbfffnvK1vDZ2dnk5OTwww8/MGvWLABGjx5dOmJz\nRsloyxjDP/7xD3744QdsNhv79u0rHbmVP76i41q3bu30Nd1BOyil6oEf9/1IUqukU5IjyhMRBrUd\nxHdp39GWcQh+HozQc1qGBbs9lT45OZnDhw9T0RrNyragqGz7C7vdzooVKwgJCTmtrdqMYHJyckhN\nTaVz587MmDGD9PR01qxZQ0BAALGxsadsyVHC2eM8TauZK1XHFZLJjswdnB11drXHDmo7iOyCbE5I\nqvsD8yJ3p9Jv2bKF4uJiIiMjz7itESNG8Nprr5V+XTJKK7vNxYIFC6rcVLBEbm4ud9xxB5dccgkR\nERFkZWXRsmVLAgICWLx4Mbt37wZO3xqksuO8TUdQStVxx2zW7aHkqOorkye3ScYmNnL9NtKoqKO7\nQ/Oaqdcnla6BmnxJT5e0WXbLdWMM06dPx8/vzEehr7zyCnfeeSdnnXUWRUVFDB48mDfffJPHHnuM\na6+9lj59+jBkyBDatWtXaRtDhw7FGIPdbufSSy/lkUceAWD8+PFcfPHFJCUlkZiYSNeuXQFO2V5j\n1KhRPPDAAxUe52263YZSddjVU1ewL2Aafo23s/iqxdik+psif5j/B7YdzCau4KF6vd1GXdtSpL46\nk+02qh1BicggrIKu7R3HC2CMMR1qFa1SymUMhlzbZka0GeRU5wSQ1CqJ3w69g536vcJDO6a6z5lb\nfNOAvwFrgGL3hqOUqokCOUSxZNOvdT+nz+ndsjfItHo/D6XqPmc6qCxjzAK3R6KUqrGc4lSOp95K\n+2EJTp+T2NKaRzluS3FXWG5hjPH6uhxVM2c6heTMPYHFIvKciCSLSJ+SxxldVSnlEr8fDaD4RCyf\nry5w+pzwoHCC7G04YdvhxshcKzg4mIyMjDP+hac8xxhDRkYGwcG1T+93ZgQ1wPGx7ISWwapOrpTy\ngpPVEmIBmLFyDzNW7nG6WkKIvRPZfmuwG7vTc1feFB0dzd69eytcd6R8V3BwMNHR0bU+v9oOyhij\nO5kp5WOWThzK41+sY/6GfWACa1wtoZG9I5n+S9mVtYuOTX0/3TwgIIC4uDhvh6E8rNIOSkT+YIx5\nX0T+XtHrxpgX3ReWUqoqLcOCyTdZYPwRMTWulhBiYgHYmLGxTnRQqmGqamzf2PExtJKHUsqL9mRm\nEtB0Jd3aBNW4WkKgaY2YIDZlbKr+YKW8pNIRlDFmquPjE54LRynlrG7dfyRj9wZCC0bXuFqCYCPY\nxGgHpXxatbOjItJBRL4QkXQROSQin4uILtJVyss2HdlEsGlf6/ND7O3ZcmQLxXZd3qh8kzPpOx8A\nHwNtgCjgE+BDdwallKrakbwj/H7sd4Ltte+ggu3tOFF0gl1Zu1wYmVKu40wHJcaY/xljihyP97HS\nzJVSXlJyay7EVF5AtDplEyWU8kWVdlAi0kxEmmEt1H1QRGJFpL2ITATmeS5EpVR5JR1UsL32HVSg\naUWIfwhbjmxxVVhKuVRV66DWYI2USmqL3FbmNQP8n7uCUkpVbXPGZtqFtsMvr1Gt2xBsdAzvyPbM\n7S6MTCnXqSqLT1fFKeWjthzZQrfIbuw7w8IK8RHxfL/3e9cEpZSL+X6NE6XUKY4VHmNv7l66RHSp\ndRszb0tm5m3JdGraiSN5R8g4keHCCJVyDe2glKpjth+1bsl1juh8xm3FR8QDsCOz7hSOVQ1HlR2U\nWGI8FYxSqnrbjm4DoEuz2o+gSmgHpXxZlR2UsWrbz/FQLEopJ2w9spXQgFDaNG5zxm1FBkcSERRR\nOipTypc4c4vvJxFxfrtOpZRbbTu6jfiIeJds3icidIropB2U8knOdFBDsTqpFBFZJyLrRWSduwNT\nSp3ObuxsO7rNJbf3SsQ3jWdH5g7sxu6yNpVyBWc2LKx+9zOllEccOHaA40XHS+eOXKFTRCeOFx1n\nf+5+okNrv7mcUq5W7QjKGLMbiAGGOT4/7sx5SinXS8lMAaBT004uazO+qSZKKN/kTDXzx4AHgIcc\nTwUA77szKKVUxXZm7gSgQ7jrNhQo6ex0Hkr5GmdGQpcCY4FjAMaY/eiGhUp5xc6snUQGRxIeFO6y\nNpsENiGqcZSWPFI+x5kOqsCRbm4ARKRxNccrpdwkJSvFLVu0ayaf8kXOdFAfi8hUoKmI3AJ8C/zH\nvWEppcozxrAzc6dLb++ViG8aT2pWKoXFhS5vW6naciZJ4nngU+AzoDPwqDHm1erOE5FgEVklIr+J\nyEYRecLxfJyIrBSR7SIyU0QCz/RNKNUQHDp+iNzCXLeMoOIj4ikyRaRmp7q8baVqy9lsvPXAUuAH\nx+fOyMfK/OsFJAIjRWQg8CzwkjEmHjgK3FSzkJVqmFKyrAw+d3RQJW3uzNrp8raVqi1nsvhuBlYB\nlwFXYC3a/VN15xlLruPLAMfDAMOwRmQA04FLahG3Ug1OydbsceGu3wmnfVj7U66hlC9wZqHu/UBv\nY0wGgIhEAsuB/1Z3ooj4YW182Al4HUgBMo0xRY5D9gJtKzn3VuBWgHbtar9rqFL1RUpmCuFB4UQG\nR5Y+N/O2ZJe0HeIfQlTjKO2glE9x5hbfXiCnzNc5wB5nGjfGFBtjEoFooD/QraLDKjn3LWNMkjEm\nqUWLFs5cTql6LSUzhY7hHV1Sg68iceFx2kEpn+JMB7UPWCkijzsW7f4E7BCRv4vI3525iDEmE1gC\nDMTKBiwZuUUD+2setlINizGGlKwUOjR1fQZfibjwOFKzU7FWlSjlfc50UClYW26U/NR+DhzAWqxb\n6YJdEWkhIk0dn4cA5wObgcVYc1kAExztKaWqcCTvCFn5WXQMd32CRInYsFhOFJ3g4PGDbruGUjVR\n7RyUMeaJWrbdBpjumIeyAR8bY74UkU3ARyIyGfgVmFbL9pVqMEqy69yxBqpESfLFrqxdtG7c2m3X\nUcpZziRJ1IoxZh3Qu4Lnd2LNRzU4V09dAbhuYls1HKU1+Nx8iw+sDio5Sn9GlfdpVXKl6oCUrBQa\nBzSmVaNWbrtG85DmNAloookSymdoB6VUHbAzc6dbM/jA2l23JFFCKV/gzELdf4pImIgEiMgiETks\nIn/wRHBKKYu7M/hKaKq58iXOjKBGGGOygTFYa6I6Yy3eVUp5QFZ+FodPHHZrgkSJ2LBYDh4/yLHC\nY26/llLVcaaDCnB8vAj40BhzxI3xKKXKKRnRuKMGX3kliRJ6m0/5Amc6qC9EZAuQBCwSkRZAnnvD\nUkqVKNnm3RMjqLKZfEp5mzMd1GNAMpBkjCkEjmPtsKuU8oCUrBSC/YKJahLl9mvFhMbgJ37aQSmf\n4EwHtcIYc9QYUwxgjDkGLHBvWEqpEjszdxIXHodN3J90G+gXSHRoNKlZqW6/llLVqXShroi0xqo0\nHiIivYGS/NYwoJEHYvMIXTyrfF1KVgp9W/X12PXiwuLYla0jKOV9VVWSuBD4I1ZB1xfLPJ8D/MON\nMSmlHI4VHuP3Y7+7tQZfebHhsSzfv5xiezF+Nj+PXVep8irtoIwx07Fq6V1ujPnMgzEppRxK5oI8\nsQaqRFx4HAX2AvYf209MaIzHrqtUec7U4vtSRK4DYsseb4x50l1BKaUsJRl8nhxBlc3k0w5KeZMz\ns66fA+OAIuBYmYdSys1SslIIsAUQHRrtsWvGhTnWQmmihPIyZ0ZQ0caYkW6PRCl1mp2ZO2kf1h5/\nm9s2HjhN0+CmRARFaKKE8jpnRlDLRSTB7ZEopU6TkpnikQoS5WlNPuULnOmgzgHWiMhWEVknIutF\nZJ27A1OqocsrymNf7j6Pzj+ViA2P1Q5KeZ0z9w1GuT2KBqKgyM6OQ7kcysmjZWiwt8NRPi41OxWD\n8WgGX4m4sDhm5c0iKz+L8KBwj19fKXBiBGWM2Q3EAMMcnx935jxf8kXKF+zL3eftMNiXeYKc/CJe\n+Xa7t0NRdYA3MvhKaNFY5QuqHUGJyGNYhWK7AO9gVTd/Hxjk3tBcIzMvkymrphARHMH0kdOJDIn0\neAxdJi0gv8he+vX7K9N4f2UaQf42tk7WAaqqWEpmCn7iR/uw9h6/dkkHtTNzJ71a9PL49ZUC50ZC\nl2IVhz0GYIzZD4S6MyhXahrclNeGv8bBYwe5Y9EdXtnnZunEoYxNjMLmKBYVHGBjXGIUSx8Y6vFY\nVN2xM2snMaExBPgFVH+wi0U1iSLAFqCZfMqrnOmgCowxBjAAItLYvSG5Xu+WvXnhvBfYemQrdy++\nm4LiAo9ev2VYMKFB/tgNiEB+kZ3QIH+dh1JV8lYGH4C/zZ/2Ye3ZlakdlPIeZzqoj0VkKtBURG4B\nvgX+496wXG9w9GCeHPQkKw+s5KGlD1FsL/bo9Q/n5tMyNIgebcIYP6A96bn5Hr2+qlsKiwvZk7PH\nI3tAVSYuXIvGKu+qdg7KGPO8iFwAZGPNQz1qjPnG7ZG5wdiOYzmad5TnVz9PxKoIHh7wsMeuPfX6\npNLK6ZMv6emx66q6aXf2bopNsddGUGBtkLgobREFxQUE+gV6LQ7VcDmTJPE34JO62imVN6HHBDJO\nZPDOxndoFtwM6O3tkJQ6TUqW53bRrUxceBx2YyctO41OEZ28FodquJy5xRcGLBSRpSJyp4i0cndQ\n7va3vn9jXMdxvPHbG2T4fe3tcJQ6zc7MnQhCbHis12Io6Rx3Zu30WgyqYXNmHdQTxpgewJ1AFPC9\niHzr9sjcSER4/OzHGdF+BAdYyG8H0jiUk+ftsJQqlZKVQtsmbQnxD/FaDCXp7VpRQnlLTRbcHgJ+\nBzKAlu4Jx3P8bf5MGTwFk34VeXmNuPuzhd4OSalS3szgK9EooBFRjaN0BKW8xpk5qD8DVwMtgE+B\nW4wxm9wdmLudXDwbD8CKLTZiH5yni2eV1xXaC0nNTuXc6HO9HYoWjVVe5cwIqj1wjzGmhzHmsfrQ\nOcHpi2dttiICwn7lyWu9G5dSadlpFNmLiG8a7+1QiAuPIzU7FbuxV3+wUi7mzBzUg0ATEbkRQERa\niEic2yNzs/KLZ43xp3njUJ5Z8xBfp2rihPKe7ZlWrcZOTb2fORcXHseJohMcPHbQ26GoBqje1+Kr\nSsni2Zb7UPTKAAAgAElEQVShQSS2i+D3rOYUtFjDxB8mUmgvZHSH0d4OUTVAKZkp2MRWWg/Pm8pm\n8rVp0sbL0aiGpt7X4qvK1OuTiGvemMZB/ky+pCdvTxjAm+e/SZ9WfXho6UPM2THH2yGqBmjH0R20\nC21HsL/3S2GVdJI6D6W8oUHU4quJRgGNeH346wxsM5BHlj3CJ9s+Oe2Yq6euKK0KoZSr7cjc4RO3\n9wCaBTcjPChcM/mUV7itFp+IxIjIYhHZLCIbReRux/PNROQbEdnu+BhxZm/B9UL8Q3h1+Kuc2/Zc\nnlzxJDM2z/B2SKqByCvKIy2nisoNxsCad+HjG2D+RDiR6dZ4RIS4MM3kU97hTJLE81jp5Z9xshbf\nq060XQTca4zpBgwE7hSR7sCDwCJjTDywyPG1zwnyC+JfQ//F8HbDmbJqCu9seMfbIakGYFfWLuzG\nXvkIatET8MXdsO9X+PlteONsyHLvZpwdmnbQEZTyCqcW6hpjvjHG3G+Muc/ZmnzGmAPGmF8cn+cA\nm4G2wDhguuOw6cAlNQ/bMwL8AnhuyHOMjB3Ji2te5KU1L2Hd7VTKPXZk7gCoOMV83Sfw40uQdBPc\nsw5u+hqOH4F5f7dGVm4SFxbHkbwjZOVnue0aSlWk2iw+VxCRWKyqrCuBVsaYA2B1YiJSYVUKEbkV\nuBWgXbt2tb+4MbBxFrTqCS261Pj0AFsAU86dQmhgKP/d8F+O5h3FMALBr/YxKVVOyZzmgKTtBNgC\niAmLOfWAogJY9CRE9YaLnrPWRkQnwfBHYOE/YONs6HmZW2IrmyiR2DLRLddQqiI1KXVUKyLSBOv2\n4D3GmGxnzzPGvGWMSTLGJLVo0aL2AeRlwpd/t26L2Gu32NDP5scjAx/h9l63M3vHbPYGvIkdz256\nqBqGHUd3EBceR4Ct3C66v0yHrDQYNglsZf44GnA7tOhmjazcNIoqSTXXeSjlaZV2UCKyyPHx2do2\nLiIBWJ3TDGPMLMfTB0WkjeP1Nlg1/twnJAIufArSVsAv79a6GRHhzsQ7eaj/Q+TY1pEW+C9yCnJc\nF6dSVJLBZ7fD8lcgZiB0HH7qazY/GHAb/L4O0n5yS0xRTaIItAXqPJTyuKpGUG1EZAgwVkR6i0if\nso/qGhYRAaYBm40xL5Z5aS4wwfH5BODz2gbvtMTxEDcYvnkMsg+cUVPXdbuOtoU3cVx2cuNXN5J+\nPN1FQaqGrpgTHDh2gPiIcvNPu5dBZhr0u9m6tVfeWVdBcDismuqWuPxsfrQPb68jKOVxVc1BPYqV\nYRcNvFjuNQMMq6btQcD1wHoRWet47h/AFKzU9ZuANODKmgZdYyIw5mUr42nBRLj6f2fUXLi9P36F\njUnLeYvr5l/H68Nfp3NE52rPm3lb8hldV9Vv+bIfqKDE0W8fQmAodK2ksklgY+uPsJ/fhrwsq7Ny\nsbiwODYf2ezydpWqSqUjKGPMp8aYUcA/jTFDyz2q65wwxvxojBFjzFnGmETHY74xJsMYM9wYE+/4\neMSl76gykR1hyETYPBe2zDvj5prYe/DuyHcpthdzw4IbWLZvmQuCPJ0uCm448m1WB3XKNhsFx2Dj\nHOh5KQQ2qvzkHpdCcQFsc8+2MR2admBP9l6unPqDW9pXqiLOrIP6PxEZKyLPOx5jPBGYW5x9l5XN\nN+8+yHM6X6NS3SO788HoD2jbpC13LrqzwqoTSjkrX/YR4h9C2yZtTz65YxEUHoOEam40tE2CsLaw\nyT13zOPC4kAMBaJFY5XnVNtBicgzwN3AJsfjbsdzdY9fAFz8CuQcgG8edUmTrRu35r1R7zEwaiBP\nrniSF9e8qFsTqFrJs+0hPiIem5T5b7ntKwhuCu2quT1ss0G3sbD9G8h3ffJOh6ZWJl++/O7ytpWq\njDNp5qOBC4wx/zXG/BcY6XiuboruC8l3wpp3YLtrdq5vHNCY14a9xlWdr+KdDe9w3/f3cbzwuEva\nVg2DwU6e7KFbs24nn7QXW7fs4i+w/riqTreLoTgfUr5zeXyxYbFgpPQ2pFKe4Ow6qKZlPnf9DKyn\nDXvEWjvy+Z00trvmr01/mz+TBk7ivqT7WJS2iOsXXM/enL0uaVvVf4VyGLvkndpB7VsDxw9D55HO\nNRLTH4LC3NJBBfsHE2haki/uLaukVFnOdFDPAL+KyLsiMh1YAzzt3rDcLCAYLpsKxzO4Oeu1Gp9e\nUGRn0/5sDuXknfK8iDChxwT+PfzfHDh2gGvnXcvKAytdFbWqx/IkDYBukWU6qG0LQfyg0/nONeIX\nYC2n2PGdWxbtBpm22kEpj3ImSeJDrGKvsxyPZGPMR+4OzO3a9ILzHuTsvO85+8TiGp26L/MEOflF\nvPLt9gpfH9R2EB+N/ojI4Ehu++Y2ZmyeoTX8VJVO2NLA2E5NMd+5xCpnFNK00vNO03GYVXEiI8Xl\nMQbb21Ig6Xr7WnmMs8ViDxhj5hpjPjfG1J9Z0kH3QHQ/7s57E7Krv7feZdICYh+cx6GcfADeX5lG\n7IPz6DJpwWnHtgtrx4zRMxgcPZgpq6bwyLJHyC/Od/lbUPVDnuwh2LQl0C/Q8UQW7P8F4obUrKGO\njhUgbrjNF2SiQQwpma7v/JSqiNtr8fk0P3+4dCoUF8Ls261J6SosnTiUsYlR2ByL+YMDbIxLjGLp\nA0MrPL5xQGNeHvoyt/e6nc9TPuf6+dezJ3uPq9+FquOMMeTZ0gi2lykQm7oMjB061LCDahYHEbGw\n63uXxggQbKz09+2ZFd85UMrVGnYHBdYC3pFTrP/QS8sXzDhVy7BgQoP8sRurOEV+kZ3QIH9ahla+\nNbdNbNyZeCevDnuVvbl7ufrLq1mUtsjV70LVYYeOH6KwCNL3DDs5r7nre/APgeh+NW+w3dlW7UkX\n31YOMM0RE8j2o9pBKc+osoMSEZuIbPBUMF7T5wZIuAqWPA2pP1Z56OHcfFqGBtGjTRjjB7QnPde5\n23bnxZzHx2M+pl1YO+5ZfA/P//w8hfZCV0Sv6rjNRzZTcHgYeXmNT85r7vwe2ieDf1DNG2yfDMcz\n4LBrOxLBRpCJ0g5KeUyVHZQxxg78JiJnsCFTHSACY16EZh3g05sgt/ICsFOvTyKueWMaB/kz+ZKe\nTL0+yenLRIdG896o97i6y9VM3zSdmxbexMFjujK/IesyaQF/fP0YhZnJgJyc19zzALQ/u3aNtnOc\nl7bcZXGWCLZHs+3oNk36UR7hzC2+NsBGEVkkInNLHu4OzOOCQuHKd+HEUZh9a633jqpOoF8gkwZO\n4tlzn2XLkS1c8cUVLE6rWRahqj+WThxK29YHQKzRdHCAjXEdhKVBd1vba9RGZEdo3MIt228EmbYc\nzT9KRl6Gy9tWqjxnOqgngDHAk8ALZR71T+sEGPWslQH1o3vf4kUdLmLmmJm0adyGuxbfxeSfJpNX\nlHfacZWtuVL1Q8uwYHKL0sH4n5zXzD9IS79j0LZv7RoVgXYDYbdrR1AFRXb27+2BvagJ245uc2nb\nSlXEmXVQ3wOpQIDj85+BX9wcl/f0/SP0vAK+e8plpZAqExcex/sXvc+E7hOYuXUm13x5DVuPbD3l\nmOrWXFVGq6DXDZl5mRzLsxEWdujkvGZWDrQ+q+rq5dVplwyZu51aPuGsfZknOJ4XSEH6cJ2HUh7h\nTLHYW4BPgZLd0NoCc9wZlFeJwNhXrKrnn/4JDu9w6+UC/QK5r999TD1/Kpn5mVw37zpmbJ5RozVX\nqu7afGQzITHvE9PCWPOaYzozlaesEdCZKCkum3bmf6SU/1kszEzmsf81159F5XbO3OK7E2vzwWwA\nY8x2oKU7g/K6wMZwzQxrndRH17lka47qnN32bD4b+xkDowYyZdUU+g/4khE9Ipxec6Xqpg2HrSTZ\nEHt764nf10NRnlVX70y0PgsCGrtkHqr8+j+brYjIFin6s6jczpkOKt8YU1DyhYj4Y+2oW79FtIcr\np0PGDpjlvqSJsiJDInlt2Gs8mvwo27JXs/TAN9iNcXrNlap71qWvIy48Dj8aW0/scXQotU2QKOHn\nDzH9YPeZj6DKr/+z2/3ILc4gsrETFdaVOgPOdFDfi8g/gBARuQD4BPjCvWH5iLhzrUW82xZYa6Q8\nQES4svOVfDb2MxrbWhPQ9Cci233JZX1bOL3mStUNxhjWHV7HWc3POvlk2k/QtB2EtTnzC7Q7Gw5u\ngBOZZ9xU2fV/53QzFBc0Ynf27jOPUakqONNBPQikA+uB24D5wCR3BuVT+t8Cva+HH56DX2d47LLR\nodGsuvcW2jcPpiBkJT8V/o3RZ+/X9Sf1yN7cvRzJO8JZLRwdlDGwZ+WZj55KtBsIGNiz6oybKrv+\n7/Gx3QiJeZ+NGRvPPEalquBMFp8dmA78H1bK+XTTkH5LisDoF6HDefDFXW4pwlkZm9iILB5Oh4JH\n6RDegYd/fJhbv7lV6/nVE+vS1wHQq0UvAFoW/w65B6HdANdcIDrJ2q5j75l3UGXFhccR4h+iHZRy\nO2ey+EYDKcArwGvADhEZ5e7AfIp/IFz1HjTvAjNvoH3hTo9ePsi0Yvqo6UwaMIkNhzdw6dxLmbZ+\nmpZKquPWpa8jxD+Ejk07AtC5YJP1QoyLOqjAxtCqh0tGUGX52fzo1qxbaYKHUu7izC2+F4Chxpjz\njDFDgKHAS+4NywcFh8P4TyAolAePPEJkceXlkNzBJjau7no1c8bN4dy25/LyLy9zzZfXlP4Vruqe\ntelr6dm8J/42fwC6FG6ydsRt2d11F4npb+3MW02l/prq0bwHW45soche5NJ2lSrLmQ7qkDGm7GKg\nncAhN8Xj28Lbwh8+Jdic4KEjk+D4EY+H0KpxK14a+hIvD32ZzPxM/jD/Dzyz8hlyC3I9HouqvdyC\nXLYc2ULfVierRXQp2GjdlrP5ue5C0f2hIBcObXZdm0CPyB7kF+fr3lDKrSrtoETkMhG5DKsO33wR\n+aOITMDK4PvZYxH6mlY9eD7iUVoX7Yf3L/fIGqmKDG83nM/Hfc41Xa/hwy0fMmb2GD7f8Tl2U/t0\neK0+4Tlr09diN/bSDqqRPZeYot2uS5AoEePYrmOva//L9mzeE0DnoZRbVTWCutjxCAYOAkOA87Ay\n+iLcHpkPe/SuPxNwzXtw4Df48Boo8M4W2E0Cm/CPAf/gg9Ef0LZJWyYtm8T186/XuYE6YM3BNfiL\nf2mKeXzBFmwY1yVIlIiIg0aRLu+gYkJjCA0IZeNh7aCU+/hX9oIx5kZPBlLndL0ILnsLPrsZZv4B\nrv2wdnv3uEDP5j3530X/44uUL3hpzUtcN+86mjCM7P0XcCgnTxf3+qBfDv5C98juNAqw6u11LtyE\nHRu22haIrYyIdZvPxYkSNrHRvXl3NmToH0PKfZzJ4osTkRdFZFa93m6jNhKusOr2pSyy6vYVey+r\nziY2xnUax5eXfsmEHhNIPxpKbr6dOz+ZR6EX41KnyyvKY/3h9eXmnzax2z/O2vbF1WL6QcZ2l8+Z\n9ojswbaj2ygoLqj+YKVqodIRVBlzgGlYc0/ur/dT1/S5wbrF99UD8PEEuPIdr42kAPo+uZT8oq6l\nX/+8LZj4h7/G3w+2T74IEfFabMryy6FfKLQX0q+1Y36ouIj4wi0sCRlBnDsuGO2o67d3NXQe4bJm\ne0T2oMhexLaj20rnpJRyJWey+PKMMa8YYxYbY74vebg9srpk4O0w6p+wdZ51u6/Qe3s3lS/sGegP\nEc13ENRhMuPnj+fn312f36LJFaeq7vvx0/6fCLAFnBxBHVxPsMlj5Khx7gmobR+3LNgtTZTQeSjl\nJs50UP8SkcdEJFlE+pQ83B5ZXTPgNhjzMmz/Bj682muJE+ULexYWw0UdhzJ58EQOHT/Enxb+ib8s\n+gs7jrp3GxFVuRUHVpDYMrF0/om0ldbHM91iozJuWrDbpnEbIoIidB5KuY0zt/gSgOuBYZy8xWcc\nX6uykm60bu99fifMuNJKnAgO83gYJYU9W4YGkdgugvScPC6Nv5RRcaOYsXkG09ZP4/IvLmdcx3Hc\n1us22jZp6/EYG6qMExlsObKFu3rfdfLJPSshrC2ER7vvwjH94bePrAW7LlpnJSL0aN5Ds0aV2zgz\ngroU6GCMGWKMGep4aOdUmcTr4LL/WNsmvHsR5Pzu8RDKFvacfElPpl6fBECwfzA3JdzE/MvmM77b\neL7c+SVjZo3hiRVPcCD3gMfjbIhWHrBGS8lRySef3LPSdeWNKuOmBbu9WvQiJTOF7ALvrAdU9Zsz\nI6jfgKY01OoRtZFwBQQ3hY9vgGkXwB9mQfP4WjU187bk6g+qoabBTZnYbyI3dL+Bt9e/zazts5iz\nYw6h/oNoXnSRy6+nTvp+7/c0C25Gt2bdrCcy90D2Pvfd3itRumB3FbSuXUJDRT+LfVr2wWBYe2gt\ng6MHn0mESp3GmRFUK2CLiCysSZq5iPxXRA6JyIYyzzUTkW9EZLvjY/1d8Bt/PvzxS2suatoI2ON7\nxTdaN27NpIGTmH/ZfC7rdBkZ5jd++30nDy95loPHDno7vArV5YSMQnshS/ctZXD0YPxKbrPtccw/\nuXsEVbJg18U/hwktEvAXf345+ItL21UKnOugHsO6zfc0VuHYkkd13gVGlnvuQWCRMSYeWOT4uv5q\n2wdu+toqNDv9Ytj0ubcjqlDrxq15JPkRgg8+QPGJWD5ZeYxRs0bx+PLHSc1K9XZ49cYvB38hpyCH\n82LOO/lk6o8QGAqtE9x78ZIFuy7O5AvxD6F7ZHd+PfSrS9tVCpy4xVfblHJjzA8iElvu6XFY5ZLA\n2mNqCfBAbdqvMyI7wk3fwEfXWrf8hj4Mg++3fmH4iC6TFpBfZAesmAqODuDI0QG8u6mQWdvHcn77\n87mp5030aN7Du4HWcUv2LCHQFkhymzK3ynYvt27vubJAbGVi+lm7Qx8/Ao2auazZ3i1788GWD8gv\nzifIz3trAFX940wliRwRyXY88kSkWERqOyPayhhzAMDxsWUV171VRFaLyOr0dM9ubeFyTVrAhC+h\n17Ww+Cn49EavpaFXpPzaqeAAG+MSo/jq7/25OeFmftr/E9fMu4abv76ZFftX6K6+tWA3dhalLSI5\nKvlkenluOhzeCrGDPBNE2QW7LtSnVR8K7YWsT1/v0naVcmZH3VBjTJjjEQxcjrVxoVsZY94yxiQZ\nY5JatGjh7su5X0AwXPIGXPAkbJwD74y0Jsh9QPm1U/lFdkKD/OnaIoq7+tzF11d8zb1972Vn5k5u\n/eZWrvryKubsmEN+cb63Q68z1h5ay4FjBxgZV+au9+5l1sf253gmiLZ9wOYPaa6dw+vbqi82sbHy\n95UubVcpZ+agTmGMmUPt10AdFJE2AI6PDSszUAQG3Q3XfgQZO2HqudbCXh9QsnaqR5swxg9oT3ru\nyc6nSWAT/tjzj3x1+Vc8nvw4RfYiHln2CCM+HcHra1/nRFEmm/ZncyjHexU0fElBkf2078f8XfMJ\n9gtmWEyZ/zq7l0NAI4hK9ExggY0hqvfJjtFFwoPC6d6se2kKvVKu4swtvsvKPK4QkSlYC3VrYy4w\nwfH5BMA3swbcrctIuO17a3HmjCtg0f+5fMfTmqps7VRZgX6BXN75cmaNncVbF7xFQvMEpv42lc05\ny8nJL+TxefoLCmBf5gly8ot45dvtgJW993Xq15wXc97J23tgdRQx/cEvwHPBtR8E+35x+S3mgVED\nWZ++nmOFx1zarmrYnFkHdXGZz4uAVKxkhyqJyIdYCRHNRWQvVjbgFOBjEbkJSAOurGG89UdkR7j5\nW5h/Hyx93ko3vnwahLbydmTVEhGSo5L5478zyS8aVfr8/LW5xK6dh7+fYf0TwwnxD/FilJ53MtnE\n8v7KNN5fmUaAHwR3PsroDqNPHnz8CBzcaCXNeFLsObDsZSubr8N5Lmt2QJsBvL3+bdYcXKProZTL\nODMHdWOZxy3GmKeMMdXemjPGXGuMaWOMCTDGRBtjphljMowxw40x8Y6Pnt8z3ZcEhMC412Hcv62J\n6zeSYcs8b0fltPLJFQF+hqaR2wnq8BTDPx7O0yufZvvR7d4N0oMqSzY5J/lrWjZqyTlty8w1pf0E\nGM8lSJSIGQBig1TX3ubr3bI3QX5BrNhfN9eoKd9U6QhKRB6t4jxjjPk/N8TTMPUeD9FJ1uaHH10H\nfSbAhU9DUBNvR1al8skVRXZhTKfhXJJ8Dp9s+4RPt33Kh1s+JLFFIld2uZIL2l9wxqOqgiI7Ow7l\n+uRGjBUlm4gtj9UZ33FHrzvwt5X577Z7GfgFQZSH6y4Hh0Hrs6z5LxcK8gsiqXUSS/ct5YF6vnJE\neU5VI6hjFTwAbqK+r13yhhZd4OZFcM7f4Jf3rASKM6g+PfO2ZLeUSSqvouSKpNZJPDv4WRZduYj7\nku7jaP5RHv7xYYZ+PJRHlz3K6t9XY4ypVVWI8vM7vqb892Pd76n4iR+Xxl966oG7l0F0Pyu709Ni\nz7G2gHfxtjBDooewO3v3aYu763L1D+VdlXZQxpgXSh7AW0AIcCPwEdDBQ/E1LP6BcP7jVomk4kKr\nRNKCB6HAdyeeq0quiAiOYEKPCXxxyRf898L/cn6781mYupAbF97IqFmjOOQ/lwJxbo1bl0kLiH1w\nHodyrOzC91emEfvgPLpMWuCW91VbZb8f942MJjfyJUbGjaR149YnD8rLggO/QfuzvRNk+0FQnA/7\n1ri02ZK5p+/36nZxyjWqnINy1M6bDKzDuh3YxxjzgDNzUOoMxJ4Dd6yAfjfDyjfg3wMhZbG3o6o1\nEaFf635MPmcyi69azNPnPE10aDSH/eaxI+hhJiyYwKzts6qsiF3Z/M7SB4Z66F3U3AdbPuBE0Qlu\n6nnTqS/s+gGM3aVJCjXSPhkQl6ebt23Slk5NO/HD3h9c2q5quCrtoETkOeBnIAdIMMY8bow56rHI\nGrqgUBj9PNy4APwC4X+XwOzbIbdu/23QKKARF3e8mLdHvE18/jO0LLyUI3lHeGz5Y5w38zz+uuiv\nfLnzy9PSlStbTOxr81Alisjh/U3vMyxmGPER5SrZ71hk1d+L6e+d4EIioFVPqw6giw2JHsKag2s4\nmqe/KtSZq2oEdS8QBUwC9pcpd5RzBqWOVE21PxtuXwbn3gvrP4VXk2DlVCgu8nZkZyyAZjQvHsXc\nS+bywUUfcE3Xa9h0ZBMPLX2IITOH8LfFf+OrXV9xvNBas1PVYmJfk+4/l+NFx7m7z92nvmAMpCyC\nuMGeXf9UXuwga46zqMClzY6IHUGxKWZR2iKXtqsapkqz+IwxNa4yodwkIBiGP2rV8lsw0Xr88h6M\n+qfn05TdQERIaJFAQosE7ku6j9/Sf+OrXV/x9e6v+TbtW0L8QxgcPZjLzr2AjIWB+GHNd9VEySR9\nTRJHanMOwAlJ46jfD1zX5Vo6NC03XZuRAplpVkURb4o9F1a+aSVLuPBnqFuzbrQLbcfC1IVc0fkK\nl7WrGiZnFuoqX9E83tr8cPMX8NVD1o69XS6yEitadPF2dC5hExu9W/amd8veTOw3kV8O/cLC1IV8\ns/sbFqYuRIL8aWTvwsdb9zA0ZigtGvlWncb84nz2B0zDnzDuSLzj9AN2fGt97OjlTanjzgXxs0Zz\nLuygRIQLYy9k2oZpHMk7QrNg11VNVw2PjpLqGhHoPhb+8rM1qkr90UqimHsXZJ/Ztu2eSk0vUVHN\nurL8bH70a92PSQMn8d2V3/HuyHeJKB5KgRzi/376P4Z9Mozr5l3Hf9b9hx1Hd/hElfXnfn6OfNsB\nogonEB4UfvoBW+dB8y7QzMuJsMHh1qLdHa6/FTcybiR2Y2fBLivDsrp/Z6Uqox1UXRXYyJqXumst\n9L8N1n4Ar/aB7ybDiboxQV2TNU1+Nj/6tupL66Ir6VQwmdljZ3NX77sAeOXXV7h07qVcNOsiJv80\nmcVpi71SE27G5hnM3DqTyKIRNLFXcAvyxFGrgkPXizweW4U6DYMDa61tP1yoc0RnujXrxuztszHG\n+PzaNeW79BZfXdc4EkZNgQG3WkVnf3jOSqLofysk3+nSjelcpbKadUH+NrZOHlXFmRZB6BTRiU4R\nnbjlrFs4dPwQS/YsYenepcxNmcvMrTPxt/nTu2VvBkUNYlDbQRgMgvs2iZyxeQZTVk1hWMwwDmy/\nrOKDtn8Dphi6jK74dU/rdL71B03KIuh1jUubviz+Mh6cHkzc8vmlz9X031kp8YXbItVJSkoyq1e7\ndpO1euv3DVYntelzayuH/jdD8l+tTRN9xKHsPCbP38yXv+3Hbqw1TRf2aM3Do7tVmzZeXeJCQXEB\naw+t5cf9P7J833K2Ht0KgJ8Jo4m9G38/9yL6tupLu9B2SDW7GjuTJHGi6AQvrH6BmVtnMixmGP8c\n8s/Kd5X9eIK1F9Pft4DNB25e2O3wQhdrDurKd13adHZBNufNGEvT3D+Ruj+yxv/Oqn4TkTXGmNO3\nTChHR1D1TeuecNV0OLTFqpK+/FVY+Rb0uR4G3G5VUfcyd65pCvQLpH+b/vRv05+/9/076cfTWb5/\nOf/84XNybZt4bLm1JUjzkOb0bdWXvq36ktA8gfiI+BptV15QXMDC1IW8vvZ19uXuY0L3CdzT955T\n6+2VlZ8L2xZC4rW+0TmBFUeXUbDhMyjKB3/XbdceFhjG6M7n8smPO7GbyDqxdk35Hu2g6quWXeHy\nt2HIg/Dji7D6HVj1H+sX0sA7rGoV1Ywg3KlkTVPL0CAS20WQ7qYJ9BaNWjCu0zjeXdicHYdyeOeW\nDuw8tpY1B9ew+uBqFqYuBMBP/OjQtAPdmnWjfVh72jRuQ2bxUfYeCmflns00a+JPZn4mqdmprD20\nlqV7l3I0/yhdIrrw9oi3GdBmQNWBbJ0PRScgwcd2mOl2MfwyHXZ+D51HuLTpG7rfwIeLZxEa9jvt\nm3R267+zqp+0g6rvmneCS/5tZfz9PA1WT7N+WbZOsDqqHpda23542NTrk0pvodV0TVNtWBP1xcz5\nOW7EqyIAABKfSURBVJ/Jl17FVV2usibwc/ex+chmNmdsZsuRLfy0/yfmpswFIC97HIV5A/jDjBkE\ntzm5t2bToKYkt0nm4o4XM6jtIGzixIho/ScQHgMxA931FmsnbrBV1WLLly7voOIj4jl/wC5+2jeX\nkPynPfLvrOoXnYNqaApPwLqP4ac3IH2zlW581tXQ94/QqodHQ6ntQtiaKJ+QUaKqifrKzgnwgyUP\nJdKmcRvnOqUSxw5bcz3Jf4ELnnD+PE/55EbY9T3cu9Xl1S1+/v1n/rTwT7QqvIJvb37MpW2rusvZ\nOSgfuRmuPCYgBPpOsIrRTvgS4kfAmunwxtnwn+FWhYr8HG9H6TK1KTJb2TnLHhxO2yZta9Y5Afz6\nP7AXWZVAfFHCFXA8wy0Fifu17kfj4u4c9l9AbkGuy9tX9Zt2UA2ViFVN4PK34d4tcOEzUJALc/8K\nz8XDp3+CrV9Z237UYbVJyHBpEofdbs3/tR9kzQv6ok4XWAVk1810S/Mtiy6hWI7x1rq3nD6ntntI\n1YW9p+pCjL5C56CUtVYq+Q4Y+GerNtu6mbBhlpXdFdLMmqdKuMKqPGDz83a0NVabhAyXJXGkfAeZ\nu605QF/lH2j9G6/90Bo9B4W6tPkQE0vTokH8b9P/GNNxDJ0jOru0fVV/6QhKnSRibQEx+gW4bxtc\nOxM6DrWqVLwzyppHmXuXteC0yHcriZdX1aaKrjynQstfgSatoNvY2p3vKb2utbIM13/qluZbFl1O\naGAojy17jEJ73R6VK8/RDkpVzC8AuoyEK/4L92+3Psaea42sZlwB/+xoTa5v+MzaIVadbs8qK/ng\n7L9aoxRfFt0PWiXAz29bW4K4mD9NmDRwEhsyNvDqL6+6vH1VP+ktPlW9oFDoebn1KMq31sxs+dJK\nV984y6qKHdMfOg63qnRHJdbJW4Eu9/2z1i3SpD95O5LqiVhVR764G9J+cuy661ojYkdw1YGreGfj\nO3Rt1pWLOlRek7CgyM6OQ7kcysmr0dxfbc/zpLoQo6/QEZSqGf8ga73M2FestOQbv4Jz7rHS1xdP\nhreH/X97dx8dVXkncPz7m5kkk4QZkhAQJLyIUkh9rQKKb634AvZo1V33nFpf2uKudbeeo93VdXvo\nvp3mqK1nV+1Z1uqx2NNqsbVrV1bt+oZUQNeC+AYCQhWVJBAxQEjI20x++8fzJJlAgiEkM3euv885\n99x7n3nuzPNzhvy8z733eeDuY+Hxb8G6X7q5jwaQ7dHTs2rLC25qjbNvgcLSXLdmcE78C/fYwav/\nMaxvm/k93z7ndmYdNYtFqxfxSu0rAx4z1AFm82Fg2nxoY1DYc1Bm+LTsgvdXuCkc/rQcmne48mSV\nmxl4ypnubrbK6VkfxSKbExaS6oD757qusr95dViHEBpxL93hzvxuXO2GzRoBe9v38u1nv822vdu4\n+9y7OX/K+T2vDeW5tSM5LpvyoY3ZYs9BmewrrXR3+11xv7t1/a9fgYvvhqpZ7lrMU7fA4tlw93Hw\n62tg9U/c9BPtIXs+ZvkP4dOtbsbjfEpO4O7kLEq6JDVCRheN5uH5D1NdUc33VnyPxW8uJtWVAob2\n3NqRHJdN+dDGoLFrUGZkiLiRKY463k0FogqN78OHq+HDV9164//4uhEYWw0TT4WJp7n12OphvbFg\nKF2JQ+p+3PK8u3Nv1kKYfsHhH59rxeUuSf3hR+57GoFrUeCS1EPzH6Lm/2r46Vs/5eXtL7Po9EWc\nNPakIT2DNpIDEA+XfGhj0FiCMtkh4kZSH3MsnHqdK2vZBbXroPZ1qF3rbrx445futUjMzTw7/gSf\n6E5w4weOGpe7GD5L7Tp37e2oE2D+HbluzdCddTO88Sg8cyvc8AeIjsyfieJYMTVn1XBO1Tnc9dpd\nXP3M1ZxbdS71n1zKuEQh4xLxw3oGLVsDEB+JfGhjkNg1KBMcqrD7A/eHfucG2LnezW+1r663Tuk4\nl7DGznTXsiqnQ+UX3LNGORydnQ9Wum7LeBIWPgfJCblry3B4dxn85lo47wfw5dtG/ONaOlt45N1H\nWLppKZ+2fUpMkyTSp/HPF17G7PGzSRYmB/U+2Rjf8UjlQxtHms0HZfKPCFRMc8uJV/aW7290yWrn\nBpewdq53YwZmTutelHTJaoxPWhXToHyqW4rLRy55pdph1b3w8o+h4li4+jf5n5zATcNxwpWw4g43\noeGUM0f040oLSvnOyd9h4YkLWbl9JYteWMKe6CpueeklIhJhcmIyMypmML1sOlWJKiaOmsjRo46m\nsrjy8MdGPMBQE0Y+JJp8aOOhWIIywVdS4aaFOObc3jJVaKqDXe+5GxJ2veeWbSvh7cf6Hl+UhLIp\nUD6lN2mNngSJ8ZCY4G7uONzntvZuhw2/g9cegL0fw/F/Bpfe627TDgMRuOQeqFvnzgwXPuembhlh\nBZEC5k2ex+TOYrro5PbLSlizYw2bGjexftf6nvm7usUiMcbEx1ARr6AiXkFtQZqoJliyfiPlReUk\ni5IkC3uXRGGC0oLSPrMp2zNXwWUJyuQnERg90S3HHnAXVHuzG/9u9za/+O1Pt7pnk1IH9PtL1Ccr\nn7DiZe7ZpaJRUFACKHS2Qed+l4x2rHddkQCTz4RL74Pjzid04kn4xuOwZD784jK49gkYOyNrHx+h\ngFnjZzFrfG9PUGuqlfrmemqba6lrrqOupY7Gtka3tDayP1JPiibuef35Ad83KlEShQkShQmShUk2\nNc2mtX0G1z66lAWz9/Qkse5lVMGog7aLY8V9nmequeLEbPwnOWz50MZDsWtQ5vNFFZp3QlMtNNXD\nvnrYt8Ov/XbbXjeye3szaNodJxGXrJJHu2tek+bAjK+67sSwq38bHvGjiFz+n1B9Sa5bdEiqSmuq\nlca2RvZ17KOpo4mmjia33d7Us79k2Ul0dfVz5iydJGb+44Dvv2/TD0EPnjcrEklz+UUrehJZSUEJ\nJbESimPFbiko7rt/wFIULepzZnckgv7M1WCvQVmCMmYgqu5sSyIQLcztTRi5tucjeOxq2PG2G/h2\n3g+yejY1Ehqa2qh5ZiNPvVVHl7rnkuYfP57bFhxLcVEnLZ0ttHS20NzZ3LPd0tlC/d79/H5tnM0f\nl6AaJRpJM3ZsHZOmrqFDdvWpqwz+76sgxGPxPkmrO6FtqG0lQhEXzJzUk8z6LLG++/vbCli6qo1V\n77WiKhTFhPOqy7lt/nFMLEtQGC0c8NpdNq5bBfomCRFZANwHRIGHVPWuXLTDmEMScRM8GiibDH+1\n3N0Qsvpe2LjMXROceSlMPRvGVeddAh/ouaSqMnfH4JjiMQMe21D7Dps++ggR6NIoF0w9i5orbuxT\np0u7aEu10ZpqPWjJLN+f2t9vndZOXzfdRkp200U7q7Z/2FP2WaPCt+28HNU5ICnaU1Fe/PhpVj/9\nZM/rhZFCiqJFFEYLicfibh2Ns62wHdECbnyhsk+dgkgBhdFCKuIV3HjyjYf45OGT9QQlIlFgMXAh\nsB1YIyLLVPXdbLfFGHMYogXulvNZC2Htz9w0LL/3t6AXV7iuz4ppMLrK3dgSL3N3UMaTEC1yD15n\nrqOFvduRGESyP7DNUJ9LGsxxEYm4br6CkiNu5xWLV7O1oZmlt36552aHLu2iPd1OR7qDtlSbW6fd\nuj3dzp3LdrG1q5HRyb1UJY5md8tsrpo9s0+dniXVu/1B5z52159DQ9kqiDXQme6kI91BR1cHHekO\nKosrs5agst7FJyJzgX9R1fl+//sAqnrnQMdYF58xAbV7G2xb5UZAb3wfGj/o+9za4YrE3CLR3qR1\nyP3+yqJuydyXyIDLii27UCKcN/Oofl6XAfefeKMeFeHPT5t88OsceFx/7yGuHvRu97uGy5/s4M3W\nSq45rpOaOZ19XjvU8f/23GZUhFsvmjnge/ctg4t/9Qkb9ye4pjpGzdlFBx8Xi8Ok2UP/jgnwNSgR\nuRJYoKp/6fevBU5X1ZsGOsYSlDF5pCvtbjRp3Q1te6CtCdId7iaLdEff7VQ7pNuhqwu6Um7RtHuP\nrnTGfqq3Ts9+f2X+uD77KdAud01Ruw5aGva1ElGlsrSg39dB+z02lU4TQYlw8M0Iw2VG289p5+Ah\nv4roYHP8W7n5rGQV/O2GI/qsIF+D6q+j+qAsKSI3ADcATJ48eaTbZIwZLpGo6+Irqch1SwZlqINn\n9fnj2W/yO1RZOmNiSJ8A+1mvbO6kZvkOnt60l7QK8Zgwf/ooFn1lLJSuHfC4wbx398vdZStbUtSs\nbubpLW2kiRCPwfxjClg0NwklT/ceF83eAMi5SFDbgUkZ+1XAQX0Cqvog8CC4M6jsNM0YY4ZAxHUp\nMrwTdY4rh0SyjbQ2uRs50koiWc64KdXD+jngEnVi4zuk+ch/FiQqJjDui7l7fioXCWoNMF1EjgFq\nga8D38hBO4wxJvCyOcBs0AazzXqCUtWUiNwEPIv7340lqnpkHZrGGBNSD1w7q+fZpJrLR2YSyVx8\n1mDYg7rGGGOyymbUNcYYk9csQRljjAkkS1DGGGMCyRKUMcaYQLIEZYwxJpAsQRljjAkkS1DGGGMC\nyRKUMcaYQLIEZYwxJpDyYiQJEfkE+PAI36YS2DUMzQmqsMcH4Y8x7PGBxRgGwxHfFFUd+1mV8iJB\nDQcRWTuYoTXyVdjjg/DHGPb4wGIMg2zGZ118xhhjAskSlDHGmED6PCWoB3PdgBEW9vgg/DGGPT6w\nGMMga/F9bq5BGWOMyS+fpzMoY4wxecQSlDHGmEAKRYISkSUi0iAi6zPKKkTkeRHZ4tflvlxE5Cci\nslVE3haRU3PX8sETkUki8pKIbBSRDSJysy8PRZwiEheRP4rIWz6+f/Xlx4jIaz6+X4tIoS8v8vtb\n/etTc9n+wyEiURF5Q0Se8vuhiVFEtonIOyLypois9WWh+I12E5EyEfmtiGzy/x7nhilGEZnhv7/u\npUlEbslFjKFIUMDPgQUHlP0D8KKqTgde9PsAFwPT/XIDcH+W2nikUsDfqWo1cAbwXRH5IuGJsx2Y\np6onA6cAC0TkDOBHwD0+vt3A9b7+9cBuVT0OuMfXyxc3Axsz9sMW43mqekrGszJh+Y12uw/4X1Wd\nCZyM+y5DE6Oqbvbf3ynAacB+4HfkIkZVDcUCTAXWZ+xvBib47QnAZr/9AHBVf/XyaQGeBC4MY5xA\nCbAOOB33xHrMl88FnvXbzwJz/XbM15Nct30QsVXh/nHPA54CJEwxAtuAygPKQvMbBZLABwd+D2GK\n8YC4LgJW5yrGsJxB9ecoVa0H8Otxvnwi8HFGve2+LG/4rp4vAa8Rojh919ebQAPwPPAnYI+qpnyV\nzBh64vOv7wXGZLfFQ3Iv8PdAl98fQ7hiVOA5EXldRG7wZaH5jQLTgE+Ah3037UMiUkq4Ysz0dWCp\n3856jGFOUAORfsry5l57ERkF/Bdwi6o2HapqP2WBjlNV0+q6FaqAOUB1f9X8Ou/iE5FLgAZVfT2z\nuJ+qeRsjcJaqnorr9vmuiJx7iLr5GF8MOBW4X1W/BLTQ29XVn3yMEQB/LfRrwOOfVbWfsmGJMcwJ\naqeITADw6wZfvh2YlFGvCqjLctuGREQKcMnpUVV9wheHLk5V3QOswF1rKxORmH8pM4ae+Pzro4HG\n7Lb0sJ0FfE1EtgGP4br57iVEMapqnV834K5bzCFcv9HtwHZVfc3v/xaXsMIUY7eLgXWqutPvZz3G\nMCeoZcA3/fY3cddsusuv83eenAHs7T5tDTIREeBnwEZV/feMl0IRp4iMFZEyv10MXIC7+PwScKWv\ndmB83XFfCSxX3wEeVKr6fVWtUtWpuK6T5ap6NSGJUURKRSTRvY27frGekPxGAVR1B/CxiMzwRecD\n7xKiGDNcRW/3HuQixlxfhBumC3lLgXqgE5fNr8f11b8IbPHrCl9XgMW46xvvALNy3f5Bxng27rT5\nbeBNv3w1LHECJwFv+PjWA//ky6cBfwS24roainx53O9v9a9Py3UMhxnvV4CnwhSjj+Mtv2wAFvny\nUPxGM+I8BVjrf6v/DZSHMMYS4FNgdEZZ1mO0oY6MMcYEUpi7+IwxxuQxS1DGGGMCyRKUMcaYQLIE\nZYwxJpAsQRljjAkkS1DGGGMCyRKUMcaYQLIEZUyOichsP49O3I/GsEFETsh1u4zJNXtQ15gAEJEa\n3MgRxbix3u7McZOMyTlLUMYEgB85eg3QBpypqukcN8mYnLMuPmOCoQIYBSRwZ1LGfO7ZGZQxASAi\ny3BTcByDm430phw3yZici312FWPMSBKR64CUqv5KRKLAKyIyT1WX57ptxuSSnUEZY4wJJLsGZYwx\nJpAsQRljjAkkS1DGGGMCyRKUMcaYQLIEZYwxJpAsQRljjAkkS1DGGGMC6f8By9ad3hpqvVkAAAAA\nSUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "A = 37.186907 +- 3.394303\n", "mu = 436.793559 +- 2.024171\n", "sigma = 26.175143 +- 1.867333\n", "B = 45.288966 +- 5.719820\n", "Covariance matrix\n", " [[ 1.15212918e+01 2.22599550e-01 -3.79884374e+00 3.04784084e-01\n", " -9.42494871e-01]\n", " [ 2.22599550e-01 4.09727004e+00 -3.66278271e-01 3.16162693e-02\n", " 4.77274651e-01]\n", " [ -3.79884374e+00 -3.66278271e-01 3.48693275e+00 1.75051531e+00\n", " -6.23502950e+00]\n", " [ 3.04784084e-01 3.16162693e-02 1.75051531e+00 3.27163382e+01\n", " -6.82057356e+01]\n", " [ -9.42494871e-01 4.77274651e-01 -6.23502950e+00 -6.82057356e+01\n", " 1.77768105e+02]]\n" ] } ], "source": [ "def fitFunc(x,A,mu,sigma,B,scale):\n", " \"\"\"Function to fit is A*exp(-(x-mu)**2/2*sigma**2) + B*exp(-x/scale)\"\"\"\n", " return A*np.exp(-((x-mu)**2/(2*sigma**2))) + B*np.exp(-x/scale)\n", "\n", "# set initial values for the fit, based on obersvation of plot above\n", "A = 20\n", "mu = 420\n", "sigma = 20\n", "B = 60\n", "scale = 90\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", "# have to correct for the fact that some errors may now be zero -> inf chi2!\n", "errs = np.sqrt(s_data[1])\n", "for i in range(len(errs)):\n", " if(errs[i]==0): errs[i] = 1\n", "\n", "# plot data\n", "plt.errorbar(cbins,s_data[1],errs,fmt=\"*\",label=\"Binned Data\")\n", "# plot initial fitted values\n", "pbins = np.linspace(xMin,xMax,250) # many small steps for smooth curve\n", "plt.plot(pbins,fitFunc(pbins,A,mu,sigma,B,scale),label=\"Initial Function\")\n", "\n", "# do the fit\n", "popt, pcov = curve_fit(fitFunc,cbins,s_data[1],sigma=errs,p0=(A,mu,sigma,B,scale))\n", "# plot the fitted output\n", "plt.plot(pbins,fitFunc(pbins,popt[0],popt[1],popt[2],popt[3],popt[4]),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()\n", "\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" ] }, { "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 }