Gustav1-SB.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.14/00\n"
     ]
    }
   ],
   "source": [
    "import ROOT\n",
    "import numpy as np\n",
    "#from root_numpy import random_sample\n",
    "from scipy.ndimage import uniform_filter\n",
    "import time\n",
    "import matplotlib.pyplot as plt\n",
    "import multiprocessing as mp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Example simulation</h1>\n",
    "\n",
    "Simulation of particles and how they compare with noise visually"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "61\n",
      "\n",
      "1.0000042332535612 1.0700208431186262\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Warning in <TCanvas::Constructor>: Deleting canvas with same name: aaa\n",
      "Warning in <TROOT::Append>: Replacing existing TH1: randomnumbers (Potential memory leak).\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAxwAAAI8CAIAAAD0vjrdAAAABmJLR0QAAAAAAAD5Q7t/AAAgAElEQVR4nO3dW5KjSJoGUBibdTVg1mupXkUVMD2bmNxKmTVQG2MevNPbEyRCF0eBxDlWDwpxEXiogi9/d5xynucCAIDn/Nd3HwAAwCcQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqgFuN4ziO41t/6LecAnASpRnVgRuVZVkUxYv/aIQPHYahruuj7Q0gpVIFfKy6roUn4GX++7sPAHgbbdt+9yHcYRzHaZoWbw7DMI6jpAXsQagCbtV13XcfwrPUroD9CFXAf8Rx3A+EjxC5Fhs+s8Nt6ZDzLDu//VD3Oyngvc0A8zwMw/rvQ9u26TqLPxpVVRVFMQzDulvw4g6rqlp8XPrOtU8JPw7DsH2oi3Wu/a1brznP88VuzcU6cSfhrDfWBE5LpQooxnFsmia8DgkjDEjq+774qtcvbFhVVajZhE3CTXaLvU3TlGU80+JQ44c2TTP/TD8h+oQxVesYlKrrOq6WnkLTNOubBMPKG58LnNp3pzrg+4XYsagbxSwS31n8GFdISzVp1SfdW3gnfsTDlaq4/4tbpaW1WNBarxb3FtdZVJuunftib9c2B87JlArAvys6i6pM13Vt2355x18s8MStwotFJ13eOwcv9i0+sJ9wtItTKIri2gShi7MwoApICVXAvxNJ6MyK6rruuu7LO/6uBYudAkfXdfM8r3e+nj3hFmGrjXNcLFqv+ViYAz6SUAX8JyuUZVmWZdd1x3+WS7gFr+u6uq7jEK7HXMx/0hJwL6EKKOq6HoYhxoi+75umKcvymN1bXdeVZdk0TdM0fd8/VqMqrvfxBeHcjx8ugeMQqoCiKIq6rsdxnOc5TVfTND1ZBMquruvQTVlVVTo8/IHC0peTUT16jMBJCVXAL9J0Fd45VLwIdam8T5vZOMFj1uqAYxKqgH8PpVq8Wdd1KP/sF6rWPXfb4+LTKdS/3NUtNk7w4h2RABuEKuDf1oHm9cFicQfijR5+KGGc7XORq+IpC1XA7cyoDhRt2/Z9H7JFjBEx3+wRLOI+y7Js2zb0OX6ZqC5uFaZrD++HH9cTbl17SF+4z3GapqZpFjOqF5dmwwLY8k2TjgLHcnGg92LG88UfjbDJ4vmA69WCMG3m+vF/679Ii83Dj+s50FPDMKTvL7ZN31zsbePcrz3772K7mVEdmOe5nD2yCiiK4meZJ/54rbqTV+y5u+vjLm4Vjz/tDYyvvxytlQ7Y0usHPECoAgDIwEB1AIAMhCoAgAyEKgCADIQqAIAMhCoAgAyEKgCADMyonl+YDsdUNwBwKipVmZVlGaYQbJrm4eeRAfAWSj7dfd8Hk39mNI5jeJRYeN00jeYF+GBl6TL6ye79/apU5RQe7xpej+N48YFiAMBHMqYqv1CjKpJHwwIAH0+lKr+6rud5btvWQHUAXmMcx/LnoN73VZblWw9HFqq2pN15qa7rws19i999+qD7ruumadr9EAHgg7x1LhSqrhrH8WIqKsuy7/vwuu/7tBwVBqrH13sfIQBwHELVBSEbhXFRCyEztW0bilJt207TFINUqE6VZVnXddM0wzC88KgBeAPrf3J/+Y/wjRWe2fbiCmmXy8WVHy4ZrLfd3lXepS+qdMyspO0zDMN60eKdqqrSd4ZhWGwFwEe68TIarhTxX9rhzbZtF1fkeO0Iaw7DkN5FvrjWLG4wD3tLrz6LFdY7Xyxdv3Px4BdLw6ekBxb3H89xve3i3BcNtVi6uKR+eV5hhbZt45sbe/vyF3f7yvM8nytUhd/T+v2LrZx+LdI1L36t8x4nAG/hxr//aTBKL/YhrKQhI6wfo0BcIQ0K889LT1x6MQmF9Td2Hpamm6frx4tduvP13m4JVesjiXtbnNfFI1mkru3zijsMB7C9t1t+cbevPJ8tVM2XctW13HotVMXf/bUdAnASd4WqdRkpXSd9Z1HTijsJQSemh2sfsa5apZusN18km/nXbpntg78lVKV7Duune1sHuHRpuodr57UoSi2aJf30ey/Z917fTzemquu6tm3jxPNlWQ7DcOPcB8aeA/CwxY1N86+jTdbW/YNx2+LnGN+LK4e7qdKPC6/Tu6/WS9MdrievXq9/u/X623tIl4ajCqd87bzi3WPFpUZLl3Zd92WzP+OMk3+G31DIVbcnKgDIKz7Z7K4peMImi4tXXddpelinoqqqroWqvd31WRcfRnL7kS+Wtm3b93244of5I3c98dNVqoKYx+9q3IsrK18B8KVFVui6Lk7QU9f1Ykx6cO0K9fpawKGeulb/qm3bayW9oii6roujqfq+b5rm3mck3+WMlariZ69fmH/23kqgFAXAk0KcSi9AdV3fW69a9CemS9e7etl81HtcJdNUd++U6+lk3V3X9X3fdd1O87afsVIVx1EtxlfdYlE+LYpimqZDRXgA3sLi2nF7FlkPgSp+HTkU9rze4WuuVk+GqsVFNu3rvFiRCnNDXjuSxXNv0hFaezhdqCrLMn0q3725Kvw+4uYXv9kA8KXFOKHbK0l1XYd/4V8byhJCQ9M04UXIFsU+V6vFdTDLI9riQwzHcQwTccc602Lp9nnFYewxRe191T5XqApzoC9aM+SqG3Nr6L4N06aXZTlNkwcnA3CvcP9/+VPxswxz40OR491w8WK0nmCzKIowhCjkkp1uzAqX1Hgkfd+v59u8S2iHxZHHpevz2r4Kp3t7wVX77hFFBBdvvgDgVB4YmJuKvR/hahIe5HJ7HSWsv3FH25cr5JL9g7Z3eO/HPXx49/5+hSoAeNCToYqDu/f3e667/3a9kRIAOLNzhari1/tXAeAZ/q1O6lwD1QEAdnK6ShUAZKRYRXS6ULX97dc5CMBdXDg+2L2J+XShyrcfANiDMVUAABkIVQAAGZyu+w8AXiDM4n1x0ZPPnvNIj+Oaz+Rs5wvArjYuK1VVPXzlHYZhGIbtPT92wDsJj9gLqqpaHHz68L6qqm5cNP/ahrfv8/YN27bdPq972/lYv5W9He1bCMBb+zJUPbbb7W23I9frhZgSslRMVxtL0wwUc896UWiEtm2HYQiv41lv7zPd8MuD2c5VQtUWoQqAjL4lVB3N4mhDcIlhZZF4QpS5uGFYFJPTIvGk+7m44cVP397wy3YWqra80XcUgON7OFSFa3/aFbWIGrGEE9aZfxZ15l8zxLzqelsfw7U+slyqqlrUe2Ksiad5ccP1UcUNF+c4/9qez2yYHur24c33xwZ3/wHA92iaJg0fYex513WxtyuORi/LMnRaFT8Hqgdd1/V9H7q6hmGYpinOV1mW5TRN4f22badpujZw/knjOKZD78PrcORxTH1YZz1Cfz3c/tpBLta8tuF687jmeoB/epBZnO7uPzOqA/Ay64tOVVXxKt62bYwgVVVN01QkV/308h/XXAiJKi4ahqFpmnEcw7Zt24YXr7lVMJ5s/NygrutwauGA07B4zcVs1Pd9fD/jGQlVjxObAHiZRQffwrXksb1mtFF6CS/6vi9+lo6enMfhFqHrre/7vu/Tj5umKV58y7JsmuZTr8W6/wBgL/Ulj+1n/WYIVU3TlIn4fhhO1Pd9eP8FoSp08IXAlB5wmizXPZg3igkye59dxqKXUAUAbymkgfVw6Rg4xnEMA7ZCutqpE7DrukXEiV2Z6XGmr9MjvLjPBw71Wkfnosdwp4FlgVAFAG/pYkoIo8LXLxZBJ6N0qFMwTVMYa78+wrTLcnFIYdF6nHv8lPg6ZMSL+1xvGA8mXTnI3zF6172C7+5s5wvArjYuK3EKyrX50s381yZwilMqLPacvo77uTbvQJz08qmzvWJxDItJNdOli9k41zNapWea/riYwurhDddTYW23yb2x4VwhQ6gCIKMvQ9VFcfapa6EqTj2VzlO12HN6DIudL3bygjLK4oPW01ZdO4zFQS6mjNrY52LDdNHi5oDFhovfy5fndVsD/Fu5PtsPVpbnOt+MyvJHfD3Pv33jkQAcx36Xlbuemhwe3nxxFPwecxBcO4aND7pl6R4bPrDP1L2/33OFDKHqYWX5I2Sp+AIAl5XPdu/v10B1AIAMTjf5pxnVAYA9nC5UiU0AwB50/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBQH7jOHZdt3jScNR1Xc7n+HIMQhUA5Nd1Xd/3TdNcW9T3/euP6sjqui4vSZ8qE9/cfsJMumaaay/uf3sCy7ucbp4qAHil8FS+9B1x6qJ1YW8cx2maQuvVdT1NU9u2dV2P49j3fXix3k8ISWHNruuaphmGIexk/ZDpvu83Hn19t7sev/zuzna+GRXF/y1eALBxWYmX6rZt11uFpfse3PtLW2/RkiEerTcJ7w/DEN/ZaOphGLZ/C/f+jk7X/Xet9Je3AAgARVG0bbuoS4WhVIva1TiOFzu8wvrpojgSK2xS/Hpd2/VcXiw0UTrybLvLLwhNl665sVXTNOva1TNOF6q2M+Z3Hx0AHyVkgjQk9X2/uJCP49g0TVVVwzAMw1BVVRyJFfq52raNixYRrSzLUJUJRZdbYsdbCB1/4aSCeZ7Ts7t2B8DGDhfvrEPb804XqgDgxRZX9EX0icGrrus4TiisM45jVVVd14VF64gWRg6F9auqmqZpv7N4pa7rqqpaZ8RY0ltErii2W3zn2gi2a3t4hlAFADtKewAv9v2F8dfpOzEepWO3wxwNi51/TGkqFcpUF2tIdV23bRuGSV1cIbzZNE3IoNe6RGMSzXTIP901AutlQhV08WZox6qq1oP+vlwaHPZ8j89AdYC1jctKOj66+Dl0uvg52jodZ719gV7fm5b29137xLcWrubb64TTv3bFjyE19JwWvw5dn5PfyLZ72/OIlapw2+SiWFqWZUz64UbK25cCwPeKZZWLV6iL4aBI5hFI3/x46zLVukq3faHvui40Vyz1petfrBdmcbhQFYp+izfD+bdtO47jOI5t26Ytvr0UAL5XemG6eC1fD7oKK4dEld7xt+dhHsK1Vur7/pbTv2VWsMxzU6Xuqmu9QPFz9o51MF+vdsvSxfu5j/csdP8BrG1cVhadceFSFS9tafdfeL1YFLsL08tZ2Mlnd/9dO4vFtT5dbdEVmDbaetqqeTXl1YZ72/NYlaqyLKuquhhFF6FycY/D9lIAOIKLZapwp1vf9+G+tjCHQpwBfJqmOA1ViAhN03xwyWo9bD8IySk2xca9e8MwxEZLGzMIlbC9urPuimC7WozaW1SqFqFyY+X5+kSrhzrf96JSBbCW97ISJqP68v1bRlh/quGnG9d88uPu/f0e5dl/YX6zi6lzO4/fm9YfnnB2PsfwQAC+y7Wh04v3z3wz1u3n/i2tdJRQFSaTfaAJ7t1ENgIA9nCIUBWH+qd9nOH+yTCH7HqTD+5OBgDe0SFCVbC46XGapmmaYqK6qxNQ5AIAXuwQd//FSbrSoWFh7HkIVeu7+aZpinf8bS8FAHiBQ4SqLy2mAls8WXp7KQDAC7xHqAoPUIzTToRJLNIUtbEUAOAFyve6G279BJ/blxZFUZZvdr7HUZY/5vm39AUALiuf7d7f74EGqt9iu/6kOgUAfJc3C1UA8EbGn+qf7tq2uFQvCDuMP967W3b05Azu70VrPMxjagDWti8c6V3qFy80289RKX59mnL6/trFNd9R2laL9kmfCXjL+VZVtV4tfXbLlzu5Nxi8x0D1jL5sPgB4Xtd14capeZ7HcYyX87Sq9MyjkdOLV7xbK8eBf6f4pORhGKqqStun67rwdORhGML5btfnxnFcTLcUdhKe4HLjTu52VwR7d2c734xUqgDWNi4roeKyeDPUWtLNN4pVxWalavFmSGwhw72pdYulLbBojUVLpjYKWot3QqNt/wruOoXTVaoA4Lt0XRcu5OM4hsJS0zTpDEHlT/fOthh2shhrtd5bWZbrxzMfpMS1nrg7FJM2NrlY5wsTLaXRKl05bdh1oz3rrgj27s52vhmpVAGsbVxW4kX9WvUo1pZCpSTkifBjHPdze6Vq/rXSE9YJuwpHEg5jXeDZOMIXW59vekbhdawq3ZJhihvqUtcaOS699ejD+net/e6EqocJVQBr25eVRbEkDOVZbJ6mhHUCeCxUrZNTXPSvf/2rKIp//etf6RF+dZYvcjHwpe8s6lhf7nDdgMWlPsSMoUr3HwDsIjzZNtSKwmNq086+1LWeqces+7Pi3pqmKYqi7/u45nEelRtOP3RQdl236JSMA/9j7e2BJooj+ruuq+s6tkM2d0Wwd3e2881IpQpg7d7LyroPK7wOdan1zh+rVF3LSeGz0oJQkVStDiIcfLxBLxzqehj+LQPzLzZguK8wLrrWyHEPdx28ShUA5FeW5bpiFN65NjL6yRHT6UDvi0EhlHZCQWgYhvAi1K6OI0w/EaZLjQ2ScYx52G34iHRXWQhVALCL9R1810LAOiLcGxfSwFTX9eKmuZBR0nf6vu/7/jh9f0VRhC659J0YE3Pdphc6FuOPXz4y+G531bXe3dnONyPdfwBrG5eVOAg6dvbFe/rSzePSi3e3bXT/tYn4TrpO3HY95j1ucqi+v8VxXmyQi811rStw3YDrrbaDwb2x4VwhQ8R8mFAFsLZ94VhPlVT8ekv/Ijmlq1U/XfzQ7d3Ovz6M5WI4O+ZVb9Fi6Uktzihdenuomlett3089zZRuf6AD1aW5zrfjMryxzz/lr4A4JbLSnygcnGppyn0ecX306cvP394eff2Mhu9clnO6PbRVPfGhnOFDKHqYUIVwJrLyme79/droDoAQAb//d0HwKGV5Y/vPgQAeA9CFV/Q2Qew4SBPI+YIhCoAeJABVdvONubMmCoAgAyEKgCADIQqAIAMTjemantE4am6fgGAjE4XqsQmAGAPuv8AADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMjjd5J9mVAcA9nC6UCU2AQB70P0HAJCBUAUAkIFQBQCQwYFC1TiOXdfVdV3Xddd1i6XlymKdjW0BAPZ2lFA1jmPTNH3fhx/7vk9v0xvHcXvzsizTbeu63uUoAQCuOEqoapqmKIp5nsdxHMdxGIaiKBY1p2EY5kRcGl60bRu2bdt2mib1KgDglY4SqoqiaNs2vg6lpligCi+u1Z9CjWqRsb4sbgEAZHSUeaqGYdjos0vT1TiOYexUukJVVYsfp2nKfpAAANccJVTFkBTyU6g2Lbrw4iirvu+rqkprUYuMVde1UAUAvNKBuv+CpmmappmmqW3bGJVCQmrbNoymCqOmFl2EN1rfRXijvKcJAHyYo1SqojBWveu6dKTUonOw67pxHEPSuvdGP4+pAQD2cLhKVVEUdV2P41hVVZwlYZ2cHitTAQDs5BChKow9X7x5Vwlqka6ELQDgxQ4RqoqiWM8slQaj9fzp6SQL63v9pmla3A8IALCrQ4SqkI36vo9Bquu6NBiFrsB0Jqowkj3+WCSVrfDC5J8AwCuVxxm4vbjDbjFpwmJp27ZpbIoD2y8uTXdynPN9C2X5Y55/234HAC4622X3WGd7bW7PuDS8uDbcanvi9eJ8v93nCVUAPOxsl92Tne3JfrvPE6oAeNjZLruHGFMFAPDuDjf5596250Y/VaAGADI6XagSmwCAPej+AwDIQKgCAMhAqAIAyOB0Y6p4Xln+iK9NrwAAgVDFfdIUlaYrADg53X8AABkIVQAAGQhVAAAZnG5MlRnVAYA9nC5UiU0AwB50/wEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkIFQBAGQgVAEAZCBUAQBkcLrJP82oDgDs4XShSmwCAPag+w8AIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACCD082o7jE1AMAeTheqxCYAYA+6/wAAMhCqAAAyOFD33/hTURR1XXddt1ih67qHlwIA7Ko8yBijcRybpimKoqqqoiimaSp+Hf8UBpjHpVVVhQh1y9J0tYOc77soyx/z/NtjSwE4ubNddo/S/RcS1TzPoVg1DENRFLHgFF60bRuWtm07TdONSwEAXuAoEbIsy7Zt0yRUlmUsOIVC1KJwdePSxacc5HzfhUoVAA8722X3KGOqhmGo63pjhdC1l/4YughvWQoAsLejdP/FRBW68MKPaeFqEbnu+hEAYG9HqVRFYXBVURRt24ZsdHHIebS9dG17RvUNpypgAgD3OkqlKprneRiGqqr6vr9lsPm9Ran5UY+dDgBwEocLVUVR1HU9jmPIVcWV2HRvgQoAYFeHCFVxEFVq8c5dnYAiFwDwYocIVUVRrGeWSoPR+m6+MMPnLUsBAF7gEKEqFKX6vo9Bquu6NBiFvBVrV4t7A7eXAgC8wFHu/pvnuSzLeOtfURTp7J11Xbdt2/d9vHcv3hv45VIAgBc41lSncZKqa5EoPjL5gaXF+aZ2fZ4Z1QF42Nkuuyc725P9dp8nVAHwsLNddg8xpgoA4N0dZUzVy2zPqH6qQA0AZHS6UCU2AQB70P0HAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQwelmVPeYGgBgD6cLVWITALAH3X8AABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGZxu8k8zqgMAezhdqBKbAIA96P4DAMhAqAIAyECoAgDIQKgCAMhAqAIAyECoAgDIQKgCAMhAqAIAyOB0k3+aUR0A2MPpQpXYBADsQfcfAEAGQhUAQAZCFQBABscKVV3X1T8tFpUrXddd3HbxPgDACxwoVJVl2fd9eD1NU1mW4ziGH+OLW7bt+36dyQAAdnWUUBVi0DAM4ziO4xju0WuaJl1nGIY5EStS4UXbtmHbtm2naVKvAgBe6SihapqmqqrSClPbtvF1qFRdqz+FGtUiY31Z3AIAyOgooaqqqo3aUtoP2HXdOjBVVbX4cZqmvEcIALDhKJN/rnNSHCMVxcnQ+76vqirdZFHEqutaqAIAXukolarUOI4hPw3DEN4JCalt2zCaKoyaCkHq3m6+9V2EN8p7jgDAhzlKpSqKRaZhGGL9KX1dFEXoAQyr3Xujn8fUAAB7OFClquu6siynaQoVqTQtrZPTY2UqAICdHKVSNY7jeqTUXZtv/AgAsLejVKqaptlIVOv509NJFtb3+oUJGnY5UACASw5RqYpZat3NFxZVVZVORtV1XeglDOt0Xdc0TV3XadIy+ScA8ErlEQZuj+O4mDw9ioe3uP+ubds0NnVdl07BsFgaleUhzveNlOWPef7tsaUAnNzZLrvvdLYbBa10hY37Ac/2232eUAXAw8522T1E99+Nvpw9wXOUAYDvcpSB6gAAb+2dKlVZbM+NfqoqJQCQ0elCldgEAOxB9x8AQAZCFQBABqfr/uMWZfnjuw8BAN6MUMVlJqACgLvo/gMAyECoAgDIQKgCAMhAqAIAyOB0A9XNqA4A7OF0oUpsAgD2oPsPACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgg9PNqO4xNQDAHk4XqsQmAGAPuv8AADIQqgAAMhCqAAAyEKoAADIQqgAAMhCqAAAyEKoAADIQqgAAMjjd5J9mVAcA9nC6UCU2AQB70P0HAJCBUAUAkMGxuv+6rhvHMbyOL9ZL67ruuu6upQAAuzpQpaosy77vw+tpmsqyTHNVurTv+7qur227XgoAsLejhKoQg4ZhGMdxHMcwnLxpmrA0VJ7atg1L27adpimWo7aXAgC8QHmQu+HKsqyqKi1NdV3X9304vDAPQnqo6frbSxefcpDzPbiy/DHPv+VaDYBzOttl9yiVqqqqtmtLVVUtfpym6calAAB7O0qoGsdxMRAqjpEKFkvv+hEAYG9HCVWpcRxDj94wDMWl2wAXK9+18/JRj58PAHACx5pSoSiKuq5Dz90wDLcUnO4tSp2qcxcAeJkDVaq6rivLcpqmtm3neY5p6WJsurdABQCwq6NUqsZx7Pv+4i17cYXtzW9fGQAgu6NUqpqm2UhU67v5pmmKd/xtLwUAeIFDhKqYpeqV8H6YbWHRIbiY/PPaUgCAFzhK919RFBszS9V13bZt3/fxLry2bdMUtbEUAOAF3myq0/jI5AeWFueb2vVhZlQH4Hlnu+weqFJ1i+36k+oUAPBdDjGmCgDg3b1Zpep523Ojn6pKCQBkdLpQJTYBAHvQ/QcAkIFQBQCQgVAFAJDB6cZUkVdZ/ggvTFgFwMkJVTwuBqkYrQDgtHT/AQBkIFQBAGQgVAEAZCBUAQBkcLqB6h5TAwDs4XShSmwCAPag+w8AIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgg9NN/mlGdQBgD6cLVWITALAH3X8AABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGZxu8k8zqgMAezhdqBKbAIA96P4DAMhAqAIAyECoAgDI4Iihqq7rcRwXb5YrXdelK3RdV9d1XdeL9wEAXuBwoWocx2ma1m9ub1WWZd/34XXf93Vd73BoAABXHShUjePYdV3TNNdWGIZhTsSKVHjRtu04juM4tm07TZN6FQDwSuVxphhIZ5AahiGtNnVd1/f9tUMNG6ZLy7KsqupiH+JxzvfIyvLHPP+23/oAnMHZLrsHqlSF+tMwDOtFMR6FatY6LVVVtfhx3YcIALCfd5r8M5ay+r5fFKIWg6jquhaqAIBXOlClakNISG3bhmpWGDUVgtSXY9gX1ncR3miH0wIAPsd7VKrWQ6ziTYL33uh3qs5dAOBl3qNStU5Oj5WpAAB28h6h6kuLdCVsAQAv9h6haj1/eohNoV61vtdvmqbF/YAAALt6j1BVVVXf9+lsn9M0tW0bfyySLsLwwuSfAMArvcdA9XEcw4No4rNo2raNsamu67Zt+76P9+i1betJNQDAK73TVKdxpNS1wBi+ADAAAAyOSURBVJT2CV50tqldH2ZGdQCed7bL7ntUqoIvi0+qUwDAd3mPMVUAAAf3TpWqLLbnRj9VlRIAyOh0oUpsAgD2oPsPACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIPTTf5pRnUAYA+nC1ViEwCwB91/AAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGZxuRnWPqQEA9nC6UCU2AQB70P0HAJCBUAUAkIFQBQCQgVAFAJCBUAUAkIFQBQCQgVAFAJCBUAUAkMHpJv80ozoAsIfThSqxCQDYg+4/AIAMhCoAgAyOGKrquh7Hcf1+13V1Xdd13XXdvUsBAHZ1uDFV4zhO07R+Pwwwr6qqKIq+78dxTIPX9lJuVJY/vvsQAOBdHahSNY5j13VN06wXheJT27YhLbVtO01TrEhtL+Uu8/xb+O+7DwQA3kx5nLvh0skOhmGo63qxKD3UsiyrqgrlqO2li484zvkeUFn+eCxOPbwhAB/sbJfdA1Wq5nme53kYhotLQ9de+mPaS7i9FABgbwcKVdvSwtW9PwIA7O1wA9XXtoec3zsgfXtG9Q2nKmACAPd6g1C17d6ilGwEAOzhDbr/LsYmMyYAAIfyBqEquKsTUOQCAF7sPULV+m6+aZriHX/bSwEAXuA9QlWYyTP2A4YXi8k/ry0FAHiB9whVdV2HedLLsizLcpqmtm3TFLWxFADgBd5sqtMwWOpaYNpeWpxvatd7mVEdgIzOdtl9sykVtutPqlMAwHd5j+4/AICDE6oAADJ4s+6/520/puZUXb8AQEanC1ViEwCwB91/AAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABmcbvJPM6oDAHs4XagSmwCAPej+AwDIQKgCAMhAqAIAyECoAgDI4HQD1dlJWf4IL+b5t+89EgD4FkIVGcQgFaMVAJyN7j8AgAyEKgCADE7X/WdGdQBgD6cLVWITALAH3X8AABkIVQAAGQhVAAAZCFUAABkIVQAAGQhVAAAZCFUAABkIVQAAGZxu8k8zqgMAezhdqBKbAIA96P4DAMhAqAIAyECoAgDI4G1CVbnSdV26Qtd1dV3Xdb14HwDgBd4jVI3juL1CWZZ934fXfd/Xdb33IQEApN4jVAXzr2JFKrxo23Ycx3Ec27adpkm9CgB4pfcIVduVqlCjWmSsL4tbAAAZvVmoGsex67p1YKqqavHjNE0vOTQAgKJ4r8k/42Tofd9XVZVGq8UgqrquhSoA4JXeo1IVElLbtmE0VRg1FYLUvd1867sIb7TDaQEAn+M9KlXDMKS1qNAD+FgtymNqAIA9vEelaj1FQixTXZw9wSh1AODF3iNUfUmKAgC+1xuEqnEc1/OnhxQVylTre/2maVrcDwgAsKs3CFUhOfV9H8tRXddN09S2bfyxSLoIwwuTfwIAr/QeA9XneS7Lsmma+E7btjE21XXdtm3f9/EevbZtPakGAHil8o3uhouVqmuBKe0TvKgs3+l8X68sf8zzb9+7BwA+xtkuu+9RqQq+LD6pTgEA3+UNxlQBABzfO1WqstieG/1UVUoAIKPThSqxiV2V5Y+NpQacAXyw04UqyOhihLqWnMryR7q+gAXwYYQquM8iSN2ejdI1FwHr3l0BcEBCFdwkb5FpvYc0ZklXAO9IqIItL+uwizuXrgDelFAFX3hxsrmYrl5/GADcS6iCC7Zv4nuNa2OwpCuAYxKqzu4I6eGYDpVddA4CHJ9QhWvzO0nT1fceCQALpwtVZlRnw3slFSUrgEM5XagSm9j2LgFFhyDA0ZwuVMHCe1Wn1qQrgIMQquBDIoi5GAC+l1AFn8ZcDADfQqjipN691+9GOgcBXkao4rxOlS2kK4C9CVVwLtIVwE6EKs7lJL1+tzCwHSAvoYrTERoWDGwHyEKoIjOX5LemcxDgYacLVR5Ts6vDPpnuaMdzfNIVwL1OF6rEptMSCB5zbejVYikApwtVwGPW+ckId4CUUMUn0+u3KyPcAVJCFR/OBf41jMECEKqAnLbHYF1cE+AzCFV8Gl1+B7GdmRS0gM8jVHGrsizf5d7JA16k36j1XuPemwo14DO03jO0HrcTqoDv9OVNhQDvQqjiQ7gMf4yLhcay/IdZsoCDO12oMqP6B3OJ/WyL3+92QcuXAXi904UqsQk+w0ZscuMh8C1OF6r4MHr9WLv9xsMHNge4Rqi6w/P3gDy5h4ybvz6L7Nd6N14Cv7fxn3ecL8/rN39eegAPRa5/FMX/XVz/lq/fJ7Xe222eZQ/f+Oknb72381Gh6o8//pimqSiKqqr+53/+57sP5+je+p/jClTs5Now+Svvu1ER+I/PCVVhBPrf/va3oij++c9/TtMUAhaf6q1DIZ/h5irphVsX834EcAQfEqr++OOPoij+/PPPv//97+HHf/7zn3/88Yd61TfaY75sVQHe1GP/F8RK2PPf/MUeZDXYw4f0lYYyVXouZVn+7W9/WxSr3r1zOu+Yqnv/qj726fGDnjn4n9eDfxTF/z18MThO47/jATj99918vYf7I9rVIWW3b/5d/+c+v4czb36EA3gvH1KpKn52/KU//vXXX991MGQ3z79dG9cC3OX+f1A99b9eWf6jeK7Ylr1Ql2Vzf45Y+5xQVVXVdx/CGzh+99nxjxC415OZ7Mn48nwofPJhSntEupdtnmUPJzK/vz///LMoit9//z198/fff1+f3Xc3NgCcywvjwPf7nErVLWa5CgDYx3999wFkEO74WzCfAgDwSp8QqgIpCgD4Rh8Sqtb3+v3111+L+wEBAPbzIaEqDEuPNwCGF+FNAIAX+JBQ9fe///3333//66+/yrIsy/Kvv/76/fff//d//3ccx/XKXdfVdV3Xddd19y49J23ypbqufdkeExuhruuNpRpwbRzHZ9rn5K2Xuvj/r9bbUK4s2uG8rffdtx9m9ueff/7555/zPA/DUBTFMAyLFcJZV1UVqllVVd2+9Jy0yZd82R62aIRFM2rADeFbt2i9dAWtd6PQAm3bpm9qvQ3xu5dKG/DMrfdpoWqe52EY2rZd/4Ge5zm8H3/3d/14TqERYjNqkwVftmeEP6nrFBVea8BtixQVrnOPNdcJWy+K+SA9fa23LTTataUnb70PDFVpdl5c5y7+Yy7G5O2l56RNtvmyPWN9yuEvbFyqATesL0W3t4/Wi4qfY3DXhZb1arcsPYP0/9O1k7feh4ypSoUTu1ifLFZPs6mqKp2LYXvpOWmTDb5sz6iqantEhQbcMAyD1ntSWZZVVV0cDan1NsQWC6P61g145tb7wFC1bTEY9q4fz0kjPMyXbds4jouz7vs+/VEDbojnO45jbMk0Zmm9baGtLiaqQuvdoCzLpmn6vm+a5pnm+rDWO1GouvY/zy1Lz+lim3zY/wA78WW71ziOZVkWRRHKfhrwdk3TNE0zTVPbtuF/T633pXEc+76/WGPWel8KhaV0XNQ0Tb57wZs9+2/7V/LM9V5W4GV82Rbqug5/podhuKVxNGBqnufQCxPqfF/eo671iqJomqaqqgeaQusVq/9PQw/gLV14Z2i99wtVG9Xa7V/YxaVnCM4P02IP03Q3ilGgbduNrqtAA14Tplmq67rv+zAD0HodrReFb9pikqR00q/1JlovtW6i8O+idYd+cKrWe7NQ9fxEYXcVJ0/1VbhGIzzMl21b6IK5NlK40IDXXRwgHAt+cZ3tPWz8eAaLAXzTNMU+rELrPefUrfcdtxy+wsX5GMNNB+k7RXIz5/bSc9Imt/Ble8z2KWvADYtZqYK0TbTeXRaNqfU2+O5tO1eoCm8ufrtxne2l56RNbuHL9oDYAmuLFcKPGnAh/Kt4MTGv1nvMIiVovW0Xv3uxAU/eeucKVfPPX3+0iNvbS89Jm3zJl+0B1yb3Sv8VqwG3Ldpt8c99rXe7e9tH6y2+e1ovKudV65xB6MS9NrB9e+k5aZOH+bI9SQNuiJNUPdY+J2+9L2m9DXEslO9e6qShCgAgrxNN/gkAsB+hCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIAOhCgAgA6EKACADoQoAIIP/B1RPlB+6cgSoAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "canvas = ROOT.TCanvas(\"aaa\", \"bbb\", 800, 600)\n",
    "\n",
    "\n",
    "noise_sigma = 5 # mV\n",
    "signal_mu = 100 # mV\n",
    "signal_sigma = 25 # mV\n",
    "\n",
    "threshold = 20\n",
    "\n",
    "low = -100\n",
    "high = 500\n",
    "\n",
    "fNoise = ROOT.TF1(\"ElecNoiseBg\", \"TMath::Gaus(x, 0, %f)\" % noise_sigma, low, high)\n",
    "fNoise.SetNpx(10000)\n",
    "fSignal = ROOT.TF1(\"Signal\", \"TMath::Landau(x, %f, %f)\" % (signal_mu, signal_sigma), low, high)\n",
    "fSignal.SetNpx(10000)\n",
    "\n",
    "all_signal = fSignal.Integral(low, high)\n",
    "all_noise = fNoise.Integral(low, high)\n",
    "\n",
    "threshold_signal = fSignal.Integral(threshold, high)\n",
    "threshold_noise = fNoise.Integral(threshold, high)\n",
    "\n",
    "expected_signal_above = threshold_signal/all_signal\n",
    "expected_noise_above = threshold_noise/all_noise\n",
    "\n",
    "\n",
    "nbins = 200\n",
    "hist = ROOT.TH1F(\"randomnumbers\", \"simulation\", nbins, low, high)\n",
    "prng = ROOT.TRandom3(123)\n",
    "\n",
    "event_count = 3600000\n",
    "event_ratio = 0.5\n",
    "signal_count = int(round(event_count * event_ratio))\n",
    "noise_count = event_count - signal_count\n",
    "\n",
    "noise_above_count = 0\n",
    "signal_above_count = 0\n",
    "\n",
    "for i in range(signal_count):\n",
    "    signal = prng.Landau(signal_mu, signal_sigma)\n",
    "    hist.Fill(signal)\n",
    "    if signal > threshold:\n",
    "        signal_above_count += 1\n",
    "\n",
    "for i in range(noise_count):\n",
    "    noise = prng.Gaus(0, noise_sigma)\n",
    "    hist.Fill(noise)\n",
    "    if noise > threshold:\n",
    "        noise_above_count += 1\n",
    "\n",
    "\n",
    "print(noise_above_count)\n",
    "\n",
    "\n",
    "# Här vill jag räkna antalet\n",
    "print()\n",
    "print((signal_above_count/expected_signal_above)/signal_count, (noise_above_count/expected_noise_above)/noise_count)\n",
    "\n",
    "hist.Draw()\n",
    "canvas.Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "====================================================\n",
    "\n",
    "\n",
    "<h1>Simulation of error in signal estimate  vs. S/B</h1>\n",
    "\n",
    "As argument it takes: Signal sigma, noise sigma, event ratio, (and event count)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of CPU cores to use: 8\n",
      "C\n",
      "F\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZgAAAEWCAYAAABbgYH9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl0nGd58P/vNTOa0TJabI0ky5J3y4ltxUvimOw7JIESBxogKaXQwkvblxxKKS2B9gclJfQF2lIo0B5a2kKhTUKAxKFx9hgISZx4t+RVkjft+76O5v79Mc/IY3lGmpHmmRlJ1+ecOZGe9X48ylxzb9ctxhiUUkqpRHOkugBKKaXmJw0wSimlbKEBRimllC00wCillLKFBhillFK20ACjlFLKFhpglFJK2UIDjFLznIjsEJGDItIrIu0i8pKIrLT2/bWIjIlIv/U6JiK/ndoSq/lCA4xSsyQirli2xXuNRBCRtcAPgT8D8oFVwHeBQNhhjxljvMYYL/Ap4EciUmJHedTCogFGqQhEZKmI/FRE2kTktIh8MmzfX4vIEyLyIxHpBT4SZZtHRP5RRBqt1z+KiMe6xi0iUi8inxWRZuA/Jt3fIyLdIlIZtq1IRIZEpFhEfCLyC+uYThH5tYhE+v95C3DaGPOSCeozxvzUGHMu0nMbY54D+oA1s/5HVAueBhilJrE+qJ8GDgFlwO3Ap0TkzrDDdgBPAAXAj6Ns+0vgGoIf8puB7cBfhV1jCbAYWAF8PLwMxpgR4GfAA2Gb3w/80hjTSrBGUg8UASXA54FIeZ/2A5eLyDdE5FYR8U7x3CIi7wLcwNFoxykVKw0wSl3qaqDIGPOwMWbUGFMH/Ctwf9gxrxtjnjTGBIwxQ1G2fRB42BjTaoxpA74EfCjsGgHgi8aYkbBrhPtvLg4wv2NtAxgDSoEVxpgxY8yvTYTEglbZbyEYKB8H2kXkPycFmveLSDcwAOwEvmKM6Z72X0mpaWiAUepSK4ClVvNTt/Xh+3mCNYWQ8xHOm7xtKXA27Pez1raQNmPM8BTleBnIEpG3icgKgjWhn1v7vg7UAM+LSJ2IPBTtIsaYN4wx7zfGFAE3AjcRrF2FPG6MKTDGZBNsGvs9EfnDKcqlVEw0wCh1qfME+y0Kwl65xph3hh0TqTlq8rZGgsEqZLm1baprXNhpTIBgreMBgrWXXxhj+qx9fcaYPzPGrAbeDXxaRG6f7sGMMW8RbHqrjLL/DLDLuqZSs6IBRqlLvQn0Wh3wWSLiFJFKEbk6zuv8D/BXVue8D/gC8KM4r/HfwAcINreFmscQkd8SkbUiIkAvMG69LiIiN4jI/xGRYuv3y4F7gDci3UxEyoG7gOo4y6nUJTTAKDWJMWac4Df4LcBpoB34N4LDfOPxZWAvcBg4QrDD/ctxlmUPwb6RpQRrFiEVwItAP/A68F1jzO4Il+gmGFCOiEg/8CzBZravhR3zgdA8GOAt4DcE+4uUmhXRBceUUkrZQWswSimlbKEBRimllC00wCillLKFBhillFK2sCXB3lzh8/nMypUrU10MpZSaU/bt29duTdyd0oIOMCtXrmTv3r2pLoZSSs0pInJ2+qO0iUwppZRNNMAopZSyhQYYpZRSttAAo5RSyhYaYJRSStlCA4xSSilbaIBRSillCw0wNnvyQAPdg6OpLoZSSiWdBhgb1bb186nHDvLT/Q2pLopSSiWdBhgbHTrfDUBL71TLriul1PykAcZGh+t7AGju0QCjlFp4NMDY6HC91mCUUguXBhibjI0HqG7sBaC1byTFpVFKqeTTAGOTky19jPgDFOd6aOkdxhiT6iIppVRSaYCxyRGr/+WODSUMjo7TP+JPcYmUUiq5NMDY5FB9D3mZLravXAxoP4xSauHRAGOTw/XdbCovoCQvE4CWXu2HUUotLBpgbDA8Ns6J5j42leezJD8UYLQGo5RaWDTA2OBYUy/+gGFTeQHFuR5AazBKqYVHA4wNQhMsN5Xnk+NxketxaQ1GKbXgaICxwaH6bnxeD6VW81hxnofWPg0wSqmFxdYAIyJ3icgJEakRkYci7PeIyGPW/j0isjJs3+es7SdE5M5J5zlF5ICI/CJs2yrrGqesa7rtfLapHKnvYXN5PiICQElepqaLUUotOLYFGBFxAt8B7gY2AA+IyIZJh30U6DLGrAW+AXzVOncDcD+wEbgL+K51vZA/AY5NutZXgW8YYyqALuvaSdc/4qemrZ8ryvMntpXkZWofjFJqwbGzBrMdqDHG1BljRoFHgR2TjtkB/MD6+Qngdgl+7d8BPGqMGTHGnAZqrOshIuXAu4B/C13EOuc26xpY17zXlqeaRlVDD8bA5vKCiW0leZm09ulsfqXUwmJngCkDzof9Xm9ti3iMMcYP9ACF05z7j8BfAIGw/YVAt3WNaPcCQEQ+LiJ7RWRvW1tbvM80rSNhHfwhJXkexsYNXYNjCb+fUkqlKzsDjETYNvkrfLRjIm4Xkd8CWo0x+2Zwr+BGY75njNlmjNlWVFQU6ZBZOVTfTVlBFoVez8S2C5MttR9GKbVw2Blg6oFlYb+XA43RjhERF5APdE5x7vXAPSJyhmCT220i8iOgHSiwrhHtXklxuL7notoLBGswoAFGKbWw2Blg3gIqrNFdboKd9jsnHbMT+LD1833AyybYUbETuN8aZbYKqADeNMZ8zhhTboxZaV3vZWPM71rnvGJdA+uaT9n4bBF1D45yrnOQTWH9LwDFuVqDUUotPLYFGKs/5EHgOYIjvh43xlSLyMMico912PeBQhGpAT4NPGSdWw08DhwFngU+YYwZn+aWnwU+bV2r0Lp2UoUmWG6eVIMpztPZ/Eqphcc1/SEzZ4x5Bnhm0rYvhP08DLwvyrmPAI9Mce3dwO6w3+uwRpqlSmgFy41lFwcYj8vJ4hx30mswg6N++of9FFt9QEoplUw6kz+BDtf3sNqXQ35WxiX7gguPJbcG89c7q/mtf3oV/3hg+oOVUirBNMAkUKQO/pDQXJhkGfGPs6uqmda+Efad7UrafZVSKkQDTIK09g7T3DvMFZM6+ENK8jxJTRfzWk0HfcPBaUEvHG1J2n2VUipEA0yCHIrSwR9SkpdJe/9I0pqrdlU1ketxce3qQl441qJZBJRSSacBJkGO1HfjENi4NHqACRjoGBi1vSz+8QAvHG3htvXFvGtTKWc7BjnZ0p+Qa798vIV/3l2bkGsppeY3DTAJcqi+h3UluWS5nRH3J3M2/57TnXQNjnF35RLevqEEgBeONs/6usYY/t+u43zjhZOM+nXggFJqahpgEsAYw+H67qgd/BA+m9/+kWS7qprIynBy87piSvIy2bysICH9MMea+jjZ0s/oeIBTrX0JKKlSaj7TAJMA9V1DdA2OXTKDP1yyajCBgOG56hZuuaxoojb1jg0lHKrvmfW9nzrYMPFzdUPvrK6llJr/NMAkwIUZ/NEDTGGOG4fYH2D2neuirW+EuyqXTGy70Ew281rMeMDw1MFGbr+8GK/HRVVjz6zLqpSa3zTAJMDh+m7cTgeXLcmNeozL6cDn9dgeYHYdacbtdHDb5cUT2yqKvawozJ5VgNlzuoPm3mHu3VrGhqV5VDVogFFKTU0DTAIcru9hfWkubtfU/5xL8u1d2dIYw3PVzdxY4SM380I2ARHhHRtKeL22g77hma1J89SBRrweF3esL6FyaT5Hm3oZD+jQZ6VUdBpgZikQMFQ19EzZ/xJSnJtpaw3mcH0PDd1DFzWPhbx9wxJGxwP88mT8i6wNj43zzJEm7ty4hCy3k8qyPIbHAtS1JWboc6x6hsa48Wsv86sZPINSKvk0wMxSXfsAfSN+rphiBFlISZ6H1j77ajC7qppxOWSizyXcVSsWsTjHPaNmsleOt9I34uferUsBqLSSeSa7H+ZoYy/nO4f411/XJfW+SqmZ0QAzS0caghmUp+rgDynJy6RzYJQR/3QrD8TPGMOzVU1cu6aQgmz3JfudDuG2y4t55XgrY3FmE3jyYANFuR6uW+MDYLUvh8wMB1VJHklWa9WYfn2qnfOdg0m9t1IqfhpgZunQ+R6yMpysLfZOe2xoLkyrDf0wx5v7ONMxGLF5LOTtG0roHfbz5unOmK/bMzjGK8fbuGfzUpyO4MrULqeD9aXJ7+ivaxvA7XTgEHh87/mk3lspFT8NMLN0uL6byrK8iQ/fqYTmwtiRVXlXVTMi8I4N0QPMjRU+PC5HXM1kz1Q1MToe4N4tZRdt37g0j6ONvQSS2NFf29bP2mIvN68r4vG953UZAqXSnAaYWRgbD1Dd2BtTBz+ET7ZMfA3m2aomrl65mKJcT9Rjst0ubqwo4oWjsSe/fPJAA6uLcqgsy7toe+XSfPpG/JxLYlNVXXs/q4ty+MDVy2npHWH3Ce3sVyqdaYCZhVMt/Yz4A1OmiAln12z+2rZ+Trb0c/cUzWMh79hQQkP3EEebpu8/aegeYs/pTt6zpQyRi2toye7oHx4bp75riDVFXm5fX4zP6+HRt84l5d5KqZnRADMLoSWSY63BLMrOIMMpCa/BPFsVTGQ5Vf9LyG3rixGB56unbybbebARgB2TmscAKkq8ZDglaR39ZzoGMAZWF+WQ4XTwvm3lvHy8Nalr7Cil4qMBZhYO1feQl+liZWF2TMeLiC1zYXZVNbFlWQGl+VnTHuvzerhq+aKY+mGeOtjAlcsLWB7h+TwuJ+tKcqlOUg2mrm0AgDVFwcEU91+9jICBn2hnv1JpSwPMLBxp6GZTecElzUdTKclLbLqY852DVDX0xtQ8FvL2DSUcbeqlvit6/8mxpl6ON/fxnq2X1l5CKpfmU9XQk5TFzEKTOlcX5QCwojCH69YU8tje80kdaKCUip0GmBkaHhvneFNfTBMswwXTxSQuwISax+6uLI35nNBEzBenqMU8ebABl0N416alUY+pLMuja3CMxiQ0U9W2DVCan0m22zWx7f7ty6nvGuLVmnbb76+Uip8GmBk61tSLP2CiLpEcTXFuZkLnweyqamJDaV7EZqxoVhd5WVvs5YVjkQNMIGDYebCRm9YVsTjn0kmbIRtDHf1JmA9T19Y/0TwWcufGEhZlZ2hnv1JpSgPMDB2xPlRj7eAPKcnLpG/Ez8CIf9ZlaO4ZZv+57riax0LevqGEPXWd9Axdmvxyz+lOmnqCmZOnsn5JHg6BapsDjDGG2raBieaxEI/LyXuvLOeFoy2099u/kJtSKj4aYGbo0PkefF4PpfmZcZ03MZs/ATnJnqu2mseumFmA8QcMu0+0XrLvqYMN5LidvH39pTnNwmW5gxkMqhrtHUnW1jdC/4j/khoMwAPblzE2bvjpvnpby6CUip8GmBkKLZEcTwc/XJgLk4jhtbuqmlhb7GVtcfR1aKLZUl5AUa7nkuHKw2Pj/G9Y5uTphDr67VRrjSCbXIMBWFucy7YVi3jsrfNJGWyglIqdBpgZGBjxU9PWH/MEy3AXajCzCzAd/SO8ebpzRs1jAA6HcMf6YnafaL0o+ebuE630DfunbR4L2ViWT2vfCK02LkNQOzGCLHK+t/u3L6eufYA9ceRYU0rZTwPMDASH5saWQXmyRM3mf/5oCwET2+TKaN6xYQkDo+O8Xtsxse3JA434vB6uW1MY0zUqlwZTyFTb2ExW1zZAVoaT0rzIzZHvuqKU3EwXj76pnf1KpRMNMDNwuD7YJBTvEGUAr8dFtts569n8u6qaWb44mw2ledMfHMW1awrJdjsnJl32DI3x8vFW3r25FJcztj+NDRMBxr5msrr2flb5cnBESSia5XZy75Yynqlqpntw1LZyzBWN3UOcaO5LdTGU0gAzE+e7BikryMLnjZ5YMhoRoSRvdnNhegbHeK2mnbsrl8TdBxQuM8PJzeuKePFYC4GAYdeRyJmTp5KbmcEqX46tKWNq2/oj9r+Eu3/7Mkb9AX5+oMG2cswVn/nJIT7x3/tTXQylNMDMxMM7Knnh0zfN+PziXM+s5sK8eKwFf8DMqnks5O0bSmjpHeFIQw9PHmxgtS8n7r6ljUvzbEt6GZ7kcuoy5LOpPJ9H31zYnf2dA6O8UdehOdpUWtAAM0PhM8rjVZKXSfMsajC7qpopzc+cUR/QZLddXozTIfzg9TPsOd3JjgiZk6dTWZZPfdeQLc1T4Ukup/OBq5dxoqWPA+e7E16OueJFq2+uf8TP8FjiV05VKh4aYFIglC5mJt+0B0f9/OpUG3duXBK1TyIeBdlurl65iJ/tb8AY2LElemqYaCqXBms8dnT0T05yOZV7Ni8lK8PJY28u3ASYz1pzoyA4f0ipVNIAkwLFuR5G/AF6h+KfzX/ofA+j/gA3rytKWHlCq2BuXV7ASt/0NYXJNlod/XbMh5mc5HIquZkZvHtzKU8fbqQ/AZkS5pq+4TFePdU+8W/VMaADHlRqaYBJgYmhyjOYC3PgfBcAW5bNvnks5B0bS3A7HbzvqmUzOn9RjpuygixbZvRHSnI5lfu3L2dwdHxiLZuF5OXjrYyOB/jg21YA0K41GJViGmBSYDZzYQ6c62a1L4dFUyShjFf5omxe/eytPLB9ZgEGgpmV7chJFinJ5VS2LivgspLcBZkA87nqZopyPbzDypat+dlUqmmASYHQbP54R/oYYzhwrpstyxNXewkpzsuc1ZDnyqX51LUP0Dd8afLMmTLGUBchyeVURIT7ty/jcH1P0hZDSwfDY+O8cryNOzeWUJQb/PvSJjKVahpgUiBUg4k34WV91xDt/SNsXb7IjmLNSqWVuv9YU+Im+LX1jdAXJcnlVN6ztQy3y8GjC6iz/1cn2xgaG+eujaVkZjjJ9bi0k1+lnK0BRkTuEpETIlIjIg9F2O8Rkces/XtEZGXYvs9Z20+IyJ3WtkwReVNEDolItYh8Kez4/xSR0yJy0HptsfPZZiMzw0l+VkbcTWSh4bdbE9j/kigbyxLf0T9VksupFGS7eceGEnZVNS2YOTHPVjeTn5XB21YvBqDQ69YmMpVytgUYEXEC3wHuBjYAD4jIhkmHfRToMsasBb4BfNU6dwNwP7ARuAv4rnW9EeA2Y8xmYAtwl4hcE3a9PzfGbLFeB+16tkSYydLJB851kZnh4PIl8WdPtltxbibFuZ6ETricLsnlVLavWkx7/ygN3UMJK0+6GhsP8OLRFu5YX0KGleLH5/VogFEpN/PZgtPbDtQYY+oARORRYAdwNOyYHcBfWz8/AXxbgh0BO4BHjTEjwGkRqQG2G2NeB/qt4zOs15z8ihpMFxPfB8CBc91sKi+IOU9YslWW5VOdwJQx0yW5nEpoEurh+h7KF8W+2mciGGPoGRqjZ2iMvmE/fcN++kf89I9c+D24Lfj7FWX5fOzG1TO+3xt1HfQO+y/K7ODzeiYCtFKpYmeAKQPCG8HrgbdFO8YY4xeRHqDQ2v7GpHPLYKJmtA9YC3zHGLMn7LhHROQLwEvAQ1aAuoiIfBz4OMDy5ctn/HCzVZybSW1r7GvJj/jHOdrYy+/fsNK+Qs1S5dI8dp9oZWh0PKa1ZKYzXZLLqVxemkuGUzhU3807ryiddVnCGWPoGhyjvmuQ+q6hsP8O0WD9PjA69Sx6t9NBbqYLEdh5qJGb1xVRUTKzmumuqmay3U5urPBNbCv0utlzWmswKrXsDDCRPhUm1zaiHRP1XGPMOLBFRAqAn4tIpTGmCvgc0Ay4ge8BnwUevuQixnzP2s+2bdtSVvspyfPQ2jdCIGBi+gCtbuxldDzA1mXp18EfsrEsn4CB4829CRmIUNvWP+N0OB6Xk/WleRw+n7gmuycPNPDd3TXUdw0xOCmA5GW6KF+UzYrCbK5f62NpQSaLst3kZrrwZrrI9WRc+DnThccVDMCdA6Pc+NWX+dbLNfzTA1vjLtN4wPB8dQu3Xl5MZsaFoO7zeugaHGNsPDDRbKZUstkZYOqB8IkV5cDk2W+hY+pFxAXkA52xnGuM6RaR3QT7aKqMMU3WrhER+Q/gMwl6DluU5GXiDxg6BkYnhpVO5cA5q4PfhiHKiRIaSVbVOPsAE0py+d6t5TO+xqbyfJ460BhzEJ/Of752ht4hPx+4ehnli7IpX5RlvbLJz8qY0TUX57j58HUr+edf1vLJ29bGXYvZf66L9v4R7tp4ceJTn/U31TUwSvEMmhiVSgQ7v9q8BVSIyCoRcRPstN856ZidwIetn+8DXjbBYT87gfutUWargArgTREpsmouiEgWcAdw3Pq91PqvAPcCVTY+26zFO9nywLkuygqyJs5LR0vzM1mUnZGQCZdnOwZjTnIZzabyAvpG/NS1D8y6PMNj41Q39rBj61K++O6NfPSGVdy5cQkbl+bPOLiEfOzG1WRnOPnWyzVxn/tsVTNup4NbLy++aLvPmojbph39KoVsCzDGGD/wIPAccAx43BhTLSIPi8g91mHfBwqtTvxPAw9Z51YDjxMcEPAs8AmraawUeEVEDhMMYC8YY35hXevHInIEOAL4gC/b9WyJEO/SyXZNsEwkEaGyLD8hI8lCHdTxzoEJF2peO5SA7MqH63sYGzdsW7F41teaLFSL+cXhRk61xD6PyBjDs1XN3Fjhw+u5uDEiVINp79fJlip1bG2cNcY8Y4xZZ4xZY4x5xNr2BWPMTuvnYWPM+4wxa40x20Mjzqx9j1jnXWaM2WVtO2yM2WqM2WSMqTTGPBx2/G3GmCus7b9rjEnrITQXajDTf8Ns7R2moXsoLee/TLZxaT4nmvsY9QdmdZ14klxGs7bYS7bbyeH62QeYfWeDOeCutCnIz6QWU93YS0P3EHdGWBcotBheh9ZgVApp71+KhPpdYmkim5hgmYYz+CerLMtjbNxwMo5v4pHEm+QyEqcjWKM6VD/7GtW+s52s9uVQOINVTGMxk1rMrqomnA7h7etLLtnn8wabyHQujEolDTApkuF04PO6Ywsw57rJcMpEWvx0dmFtmNl9qMeb5DKazeX5HG3qnVWNyhjDvrNdXLnC3gAfby3m2apmrlm9OGLiU6/Hhdvl0CYylVIaYFKoODe2yZYHz3exYWn+RcNQ09Xyxdl4PS6qZjHhciZJLqPZVF7AqD8wqxpVXfsAXYNjbLM5wMRTi6lp7aO2beCS0WMhIkKR16Mp+1VKaYBJodDKllPxjwc4XN8zJ/pfABwOYcPSvFl19M80yWUkEx39s+iHCfW/XGVzgIHYazHPVgVXrnxHlAADwWayds2orFJIA0wKBfORTf0N82RLP4Oj42k9/2WyyqX5HGvqxT8+s2apmSa5jGTZ4iwWZWfMasLl/rNd5GdlJCTgTSfWWsyz1c1cubxgymHrhVqDUSmmASaFinMz6RgYYWyKD+LQCpbpPIN/ssqyPIbHAjOef1LXPvMkl5OJCFeUF8yqBrP3bBdXLi9IyGTNWExXiznfOUhVQ+9Fucci8WlGZZViGmBSqCQvE2OYct2OA+e6Kcxxs2xxVhJLNjsTM/pnOOGytnXmSS4j2Vyez8mWPgZH/XGf2z04Sk1rP9tWJn7+SzTT1WKeqw42j905RfMYBIcqdw6MEgjMyXywah7QAJNCocmWU/XDHDjXxdblBbNabTLZVvtyyMxwzLijfzZJLiPZXF5AwATnjcRr/7nQ/Jfk1iCnqsU8W9XM+tI8VhRO3YRY6PXgDwQzOyuVChpgUmi6yZY9g2PUtg3Mifkv4VxOB+tLZ97RX9vWn5D+l5BNy4I1qpnM6N93tgunQ9hsXSNZotViWnuH2XeuK+rosXA6F0almgaYFLqwdHLkGszB+vRdwXI6lUvzOdrYG3fzTCjJZSI71ItzMynNz+TwDCZc7j3TxcalebOa8DlTkWoxzx9twRi4+4rpA0yRV9PFqNTSAJNChTlunA6J2kR24FwXIrBpLgaYsjz6R/yc7RyM67xEJLmMZFN5ftwpY8bGAxyq705681hIpFrMc9XNrPblUFE8fQC+kI9MazAqNTTApJDDIRTnRh+qfOBcN5eV5F6SyHAu2GjN6A/NIYlVIpJcRrKpvIAzHYN0D8b+bf5YUy/DYwG2rUxdE2V4LaZ7cJTXazu4s3JJTH1yhTnaRKZSSwNMihXnRZ5sGQgYDp7vnlPzX8KtL81jtS+H7796Oq5mskQkuYwkfAnlWO09k7wJltGE12L+5Zd1+AMmpv4XgEXZwRqyBhiVKhpgUqwk1xMxwJzuGKBnaGxOzX8J53QID962lmNNvbxwrCXm8+oSkOQykivKgzWqeJrJ9llr8JTmp3aIeKgW8y+/rKU0P5NN5bENOHA4hMU5bjq0D0aliAaYFAumi7n0G+ZcWMFyOvdsXsrKwmy+9dIpguvITa82QUkuJ8vPymC1LyfmzMrGGPadsT/BZSxCtRgIzn2JZ8h6YY5OtlSpowEmxUryMukZGmN47OI13g+c6yLX40pKehK7uJwOHrytgurGXl481jrt8YlMchlJPB39jT3DNPcO257gMlb/58bV3HZ5MR982/K4zivK9dCmNRiVIhpgUqzYGunTOqkWE1rBMlnpSexy75alLF+czTdfOjltLSaRSS4j2VReQEvvSExLJOw90wmktv8l3KIcN//+kaupKMmN6zyf16OLjqmU0QCTYhOTLcPmwgyO+jne3Dsn579M5nI6ePDWtVQ19PLy8alrMYlMchnJ5jgmXO4/20W228nlS+L7QE83oXxksTZRKpVIGmBSLBRgmnsuBJgj9T0EzNxYwTIW77myjGWLs6bti0lkkstINpTm43RITCPJ9p7tYsuyAlzOuf2/SKHXw/BYgIHR8ekPXkD6hsf4+A/38ldPHuGn++o53T6gQdgGc2+CxTwTKR9ZaInkLfOgBgPB1Ts/cctaHvrZEXafbOPWy4ojHpfoJJeTZbmdrCvJnTaz8sCIn2NNvTx461pbypFMPms2f0f/yJycT2WX/ee6ef5oC26Xgx+9cQ6AguwMti4rYOvyRVy5fBGbluWTl5mR4pLObdP+xYmIE/ikMeYbSSjPgpOflYHH5aA1LKPygXNdrPLlRFwKd65675Xl/NPLNXzzxVPcsq4o4kioRCe5jGRzeT67qpoxxkQdjXXofDcBQ1qMIJut8Hxk0yXHXEhC861e/Ytb6Roc48C5Lg6c62b/uS5eOdEGgAjHi3CdAAAgAElEQVRUFHvZumwR924t49o1haks8pw0bYAxxoyLyA5AA4wNRISSsMmWxhj2n+vmxrW+FJcssdwuB//31jX85c+r+NWpdm5eV3TJMXVtAzHP8ZipzcsKePSt85ztGGSlL/IH7t6zwRQ986GJMlSDaevTkWThTrcPkJvpoijXQ3FeJpctyeX+7cERej1DYxyu7+bAuW4OnOvif4808daZTl7+zC2pLfQcFGsD829E5NsicqOIXBl62VqyBSS4smUwwDT2DNPWNzKn579E876rlrE0P5NvvnjpiLLhsXHOdw3aPiw7FMCmaibbd7aLdcW55GfN/eaRiSayAR1JFq6ubYDVvpyItdj8rAxurCjik7dX8B+/v53fv34lZzoGGPFrP1a8Yg0w1wEbgYeBv7def2dXoRaaYLqY4AfAAWv9kfnw7Xkyt8vBH9+6lv3nunm1pv2ifXYluZxsXUkuHpcjakd/IGDYfy49JlgmwuJQPjKtwVzkdPtAzINJ1hZ7CRg40x5f4lYVY4Axxtwa4XWb3YVbKEpyg01kxhgOnOsmM8PBZXN8eGw0799WTml+Jt988eIRZXU2JbmcLMPpYOPSvKgTLk+19tM37E+bCZaz5XY5yM/K0Nn8YYZGx2noHmJVlCbSyUJ/kzWt/XYWa16KKcCISL6I/IOI7LVefy8iyV2BaR5bku9hcHSc/hE/B851samsgIw5Pjw2Go/LyR/fsoa9Z7t4vbZjYnutTUkuI9lUXkBVQy/+8cAl+/aeTa8JlokQmgujgs50xDffak2RFxE41Xrp8tVqarF+iv070Ae833r1Av9hV6EWmtBcmPOdQ1Q19s7L/pdw79+2jJI8D//40qmJbXYluYxk87J8hsbGqWm79BvpvrNd+LxuVhRm216OZAnO5tcmspA6a0JvrDWYLLeTsoIsrcHMQKwBZo0x5ovGmDrr9SVgtZ0FW0iKc4MBZvfJVkb9gXkz/yWazAwnf3TzGt483TlRi7EryWUkm6zU/ZFm9O8/28WVyxfFlVAy3flyPVqDCXPamtAba4CB4HBlDTDxizXADInIDaFfROR6YMieIi08ocmWz1U1A/Ozg3+yB7YvpyjXM5GjzM4kl5OtKswhN9N1SWbltr4RznQMpnSBMTv4cty0aYCZMJPa8tpiL3XtA4zHuQT4QhdrgPkj4DsickZEzgDfBv7QtlItMMVWE9mh+h5K8zNZkm/PTPZ0EqrFvFHXyf8eabI1yeVkDodEzKy8/1zqFxizg8/roW/Yr8NsLXXt8X+ZWVvsZdQfoL5LR5LFY9oAIyIO4DJjzGZgE7DJGLPVGHPY9tItEF6PayKNx3zvfwn3O9uX4/N6+OJT1UByOvhDNpUXcLyp76JlEvad7cLtdFBZNr/Gr/hyQ+litB8mWFvuj6t5DIIBBuBUizaTxWPaAGOMCQAPWj/3GmN6bS/VAhRqJpurK1jORJbbyR/etJqOgeAHn11JLiPZXJ6PP2A41nThz3nf2S6uKM/H43ImrRzJUJhzIV3MQtc5MErvsJ/Vvvj+1tYWBacNRBoYoqKLtYnsBRH5jIgsE5HFoZetJVtgQiPJFlINBuCD1yynMMdta5LLSEId/aEJl8Nj4xyp75l3zWOgNZhwde3WCLI4a8v52RkU5Xq0oz9OsfZy/YH130+EbTPoSLKEKcnLxOWQedc8M51st4tH3lNJbdtAUhdXK83PxOf1TKSMqW7sYXQ8MC8DTFEoH5nWYDhtDVFeE2cNBmBtkY4ki1cs2ZQdwO8aY36ThPIsWB+6dgVXrVhEZsb8ap6JxV2VpUm/p4iwZVn+xFDlvWeCHfxXzsMRfIXe5DaR9Q6PpW2a+9r2ftxOB2WLsuI+d22xlycPNEyZiVtdLNY+GM07ZrMrly/id69ZkepiLCibyguoax+gb3iMfWe7WFmYTZHVnDSfZLtdZLudSclHdqqljysffoFXTky9emmqnG4bYEVhNs4Z1JbXFnvpG/FP5A1U04u1D+Z5Eflt0bCt5pFN5fkYA0caeuZVgstIfF5PUjIq/+xAA/6A4fnqZtvvNRN17QNxjyALCY0k02ay2MUaYD4NPA6MiEiviPSJiI4mU3NaqKP/6UNNtPePsm3F/B23kox8ZMYYnj7UCMCvTran3RLE4wHD2Y7YsyhPVjERYDQnWaxiDTD5wEeALxtj8gim7n/7dCeJyF0ickJEakTkoQj7PSLymLV/j4isDNv3OWv7CRG509qWKSJvisghEakWkS+FHb/KusYp65rzZzlIZYvFOW6WLc7iZ/vrgfk3wTJcoddjexPZwfPd1HcNcdWKRTR0D1Frdaini/quQcbGDatnWIMpyvWQm+nSocpxiDXAfAe4BnjA+r2P4Gz+qKyllr8D3A1sAB4QkQ2TDvso0GWMWUtwxcyvWuduAO4nGMjuAr5rXW8EuM2a9LkFuEtErrGu9VXgG8aYCqDLurZSU9pUXsCIP0BupmviG+p8lIwmsqcPNeF2OvibHZUA/Opkm633i1doiPJMJ/SKCGs1J1lcYg0wbzPGfAIYBjDGdAHT1RC2AzVWcsxR4FFgx6RjdgA/sH5+Arjd6ufZATxqjBkxxpwGaoDtJij07mZYL2Odc5t1Daxr3hvjs6kFbLO1wuWVyxcldZh0shV53XQOjNqWS2s8YPjF4UZuuayIDUvzWO3L4den0izAxJlFORIdqhyfWAPMmFWDMAAiUgRcupjGxcqA82G/11vbIh5jjPEDPUDhVOeKiFNEDgKtwAvGmD3WOd3WNaLdC+v8j4fWtWlrS6//AVTybbb6YebLAmPRFHo9BExwJrsd3jrTSWvfCO/evBSAGyt8vFHXmVb5z06395OflTGxyudMrC320t4/SvegTlqNRawB5lvAz4FiEXkEeBX4yjTnRPo6OPnrU7Rjop5rjBk3xmwByoHtIlIZ472wzv+eMWabMWZbUVFR1MKrhWHr8kX8/vUree9V5akuiq181mRLu5rJdh5qJNvt5Pb1xQDctK6IobHxiflF6aCuLTiCbDaDYStKdCRZPGJdMvnHwF8Afws0AfcaY34yzWn1wLKw38uBxmjHiIiL4GCCzljONcZ0A7sJ9tG0AwXWNaLdS6lLuF0OvvjujZQVxD/xbi7xhSZb2tDRPzYeYNeRJu5YXzKRAv+a1YVkOIVfpVEz2ekZZFGebCInmQaYmMS8Lq8x5rgx5jvGmG8bY47FcMpbQIU1ustNsNN+56RjdgIftn6+D3jZBMc27gTut0aZrQIqgDdFpEhECgBEJAu4AzhunfOKdQ2saz4V67MpNd8VWjUYO4Yq/6amna7BsYnmMYAcj4urViziVyfbE36/mRgc9dPUMzzrJSHKFmXhcTk4pQEmJrYt/G71hzwIPAccAx43xlSLyMMico912PeBQhGpITjX5iHr3GqC826OAs8CnzDGjAOlwCsicphgAHvBGPML61qfBT5tXavQurZSigv5yOwIME8faiI308VN63wXbb9pXRHHmnpp7RtO+D3jdbp99h38AE6HsFo7+mNm6wLoxphngGcmbftC2M/DwPuinPsI8MikbYeBrVGOryM4ck0pNUlelgu300F7gjMqD4+N83x1M3dVLrlkmYObKor42rMnePVUO++9MrV9XKERZIlYc6ii2Mu+s+nTt5TObKvBKKXSh4hQaMNs/t0n2ugb8V/UPBayoTSPwhx3WsyHOd0+gAisLJx9gFlb7KWhe4jBUf/0By9wGmCUWiDsCDBPH26kMMfNdWsKL9nncAg3Vvj49al2Ailey76urZ+l+VkJyVYeyklWl2aZCtKRBhilFgif15PQRccGRvy8dKyFd15RissZ+aPkpnVFdAyMcrQptakLEzGCLGRi+WTNSTYtDTBKLRA+ryehNZgXj7UwPBaI2DwWckNFsOM/lcOVjTHUtQ3MOAfZZCsLc3A6RDv6Y6ABRqkFotDrpqN/NGFZjp8+1MSSvMwpsyAU52ayvjQvpf0w7f2j9I34Zz2CLMTtcrCiMFsDTAw0wCi1QBR5PYyOB+gdnn3ndM/gGL882cpvbSqdNofbTet87DvbxcBIajrF66zsxzNN0x+J5iSLjQYYpRYIXwLnwjxX3czYuOGeLdGbx0JurihibNzwRl3HrO87E4maAxNubbGXsx2DjPqnS8m4sGmAUWqBKJxIFzP7APP04UZWFGZzRVn+tMdetXIRWRnOlDWT1bUP4HY5EpoOaG2xF7+1gJmKTgOMUgvEhRrM7EaStfeP8Juadt69aWlMiSM9LifXrF7Mr06lJm1MXdsAqwpzErocgy6fHBsNMEotEInKqLzrSBMBw5Sjxya7aV0Rp9sHON85OKt7z0Rde3/ChiiHhHKaaYCZmgYYpRaIxTluRGbfRPb0oSbWlXi5bEluzOfcWBFcGiPZw5XHxgOc6xhMaP8LBJN5lhVk6fLJ09AAo9QC4XQIi7PdtM2iiayxe4g3z3Ty7k2x114A1hTlUFaQlfR+mPquIfwBk9ARZCFrdPnkaWmAUWoBCc7mn3kN5n8PNwHxNY9BMBfaTet8vFbTwdh48kZehYYoJ7oGA8GhyrVt/SlPg5PONMAotYD4cmeXj+zpw41sKs9n5Qw+sG+qKKJvxM+h890zvn+8QkOU1yS4DwaCHf3DYwEauocSfu35QgOMUgtIYY5nxqPIzrQPcLi+J+7msZDr1vhwCEltJqttG2BRdgYF2e6EX1uXT56eBhilFpDZNJH94nBwFfJ3bSqd0fn52RlsWVbAL5M4XPl0e78t/S8QbCIDDTBT0QCj1ALiy3UzMDrO0Oh43OfuPNTI1SsXsXQWExZvWlfE4fpuugcTu/BZNHVtA7b0vwAsynFTmOPWADMFDTBKLSC+nJmliznR3MfJln7uibNzf7IbK4owBl6tsb8W0z/ip7VvJOFzYMKtKfZq2v4paIBRagHx5Qb7ItriDDBPH2rEIXD3FTNrHgvZXJ5PXqYrKf0wp0PLJNtUg4FgR39Na3/CMlTPNxpglFpAJmbzx9nR/1x1M9esLpw4f6ZcTgc3WKtc2v2hXNee+CzKk1UUe+kd9scdsBcKDTBKLSAzyah8tmOAU6393LG+JCFluLGiiKaeYdv7LuraBhCB5YuzbbuH5iSbmgYYpRaQxTnxZ1R+8VgrQMICzE3rgmljfmlzM9np9gHKF2WRmeG07R5zMcAEAoYn9tUz4o9/oEe8NMAotYBkZjjJzXTRMRB7E9mLR1tYV+JleWFiagJlBVmsKcqxPbtyXXs/q332NY8BLMnLxOtxzakA89ShBj7zk0O8eLTV9ntpgFFqgSnyemLuM+gZHOPNM50Jq72E3LSuiD11HQyP2fMt2hjDaRuHKIeICGuKcuZMgBkeG+fvnjtJZVked1cusf1+GmCUWmAKve6Ym8h2n2xlPGC4Y0OCA0xFESP+AG+d6UzodUNa+0YYGB23JUXMZGuLc+dMgPnBa2do6B7i83evT+j6ONFogFFqgfF5PTE3kb14rBWf182W8oKEluFtqxfjdjp4tqqZlt5hugdHGR4bT9jIsrq20DLJ9jaRQbAfprVvhJ6hMdvvNRtdA6N8+5Uabr2siOvW+pJyT1dS7qKUShs+r4fX6zqmPW5sPMDuE63cXbkk4d92s90u3rZ6MT/ec44f7zl30T6Py4HH5SAzw2m9HBRku3l4x0YuX5IX0/UvDFFORg3mQkf/VSsW2X6/mfrOKzUMjPh56O71SbunBhilFhif10P34Bhj4wEynNEbMd463UnfsJ/bE9z/EvL/fnsTb9R2MOIPMDw2zrB/nOGxACP+cUbGgttC+/ac7uTB/z7A0w/eQJZ7+lFhp9sGyMxwsCQv05ayhwsFmNo0DjDnOwf54etnue+q8rgWipstDTBKLTCF3uBQ5Y7+UZbkR/8AfvFYK26Xgxsr7GlOKSvI4revKo/p2FdPtfO739/DV545xt/cWznt8XXtA6zyeZPSz7BsURZupyOtV7f8+nMncDjg02+/LKn31T4YpRaYWCZbGmN48VgL168pJNud+u+hN1T4+NgNq/ivN87y0rGWaY8/3T5ga4qYcC6ng9VpPJLscH03Ow818rEbVk/5hcIOGmCUWmCKrHxkUwWYmtZ+znUOJnz02Gz8+V2XcfmSXP7iicO0TTEKbtQf4FznYFL6X0LSdflkYwxfeeYYhTlu/vDm1Um/vwYYpRaYwomMytFHkr1g1RJuvzx9AozH5eRbD2ylf8TPnz9xKOqIs/Ndg4wHjO1zYMKtLfJyvmvQtnk9M/XKiVbeqOvkT+6oIDczI+n31wCj1ALjyw0lvIxeC3jpWCtXlOUnvUllOutKcvn8O9ez+0QbP3z9bMRjQkOU7UxyOdnaYi/GQG0a9cP4xwP87TPHWeXL4YHty1NSBg0wSi0wOe7g0N9oTWTt/SPsP9eV8Nn7ifJ7167g1suK+MozxzjZculaLKetIcrJrMGk4/LJP91fz6nWfv7izsumHC1oJw0wSi0wIkJhjidqE9nLx1sxBm5fX5zkksVGRPjafZvxelx88n8OXJK0sa5tAJ/XTX5W8pqEVvlycEhwqHI6GBz18/fPn+TK5QXclYSUMNFogFFqAfLleqLWYF461kJpfiYbl8Y2qTEVinI9fO2+TRxv7uPrz564aF9wiHLyai8Q7B9avjg7bYYqf//Xp2ntG+Hz71yPiP1DtaPRAKPUAlTkdUeswQyPjfOrk+3cvr44pR9Msbh9fQkfumYF//bqaX596kLq/7q2AduzKEeyttjLqZbUB5j2/hH+5Ze13LmxhG0rF6e0LBpglFqAfN7INZjX6zoYGhtP2/6XyT7/zvWsKcrhzx4/RNfAKL3DY7T3j7AqiUOUQ9YUeznTMYB/PJD0e4f75ounGPYH+Oxdl6e0HGBzgBGRu0TkhIjUiMhDEfZ7ROQxa/8eEVkZtu9z1vYTInKntW2ZiLwiIsdEpFpE/iTs+L8WkQYROWi93mnnsyk1lxV63XQOjBIIXDzU98WjLWS7nVyzujBFJYtPltvJN+/fStfgKA/97DCnQyPIktxEBlBRnMvYuOFs52DS7x1S29bPf795jt/Zvjypo+iisS3AiIgT+A5wN7ABeEBENkw67KNAlzFmLfAN4KvWuRuA+4GNwF3Ad63r+YE/M8asB64BPjHpmt8wxmyxXs/Y9WxKzXU+r4fxgKE7LAOwMYaXjrVyU0WRratAJlplWT5/fudlPFfdwtefC/bHJHOSZcjlVo6vZw43Jf3eIV9/9gSZLgefvL0iZWUIZ2cNZjtQY4ypM8aMAo8COyYdswP4gfXzE8DtEmz43QE8aowZMcacBmqA7caYJmPMfgBjTB9wDCiz8RmUmpcipYupbuyluXc4bUePTeVjN6zmujWFvFrTjkNg+eLkB5iNS/N41xWlfPOlU1Q19CT9/nvPdPJsdTN/ePMaiqy5TqlmZ4ApA86H/V7PpcFg4hhjjB/oAQpjOddqTtsK7Anb/KCIHBaRfxeRiGlNReTjIrJXRPa2tdm7JrhS6SqU8DJ84bEXjrYgArddPvcCjMMh/P37N5OflcHyxdm4XcnvXhYRvnxvJYtz3PzpYweTPqv/u7trKcr18LEbVyX1vlOx812INARlcm6HaMdMea6IeIGfAp8yxvRam/8ZWANsAZqAv49UKGPM94wx24wx24qKiqZ+AqXmqaJQDSZs4bGXjrdw1fJFFHrT49tvvErzs/jhH2znb9+7KWVlWJTj5mv3beJUa/9Ec10yjPoDvF7bwTsrl6RFctIQOwNMPbAs7PdyoDHaMSLiAvKBzqnOFZEMgsHlx8aYn4UOMMa0GGPGjTEB4F8JNtEppSKYaCKzajBNPUNUNfTatvZLsmxeVsC1a1I7QOGWy4r50DUr+P6rp3mtpj0p9zx4vpuhsfGkrVQZKzsDzFtAhYisEhE3wU77nZOO2Ql82Pr5PuBlE8xgtxO43xpltgqoAN60+me+DxwzxvxD+IVEpDTs1/cAVQl/IqXmifysDJwOmeiDeelYKwBv3zD3msfS0efeeTmrfTl85ieHkrKU8m+svqd0G/1nW4Cx+lQeBJ4j2Bn/uDGmWkQeFpF7rMO+DxSKSA3waeAh69xq4HHgKPAs8AljzDhwPfAh4LYIw5G/JiJHROQwcCvwp3Y9m1JzncMhFOa46bAmW754rIUVhdmsSYOhrfNBttvFP3xgCy19I3xpZ7Xt93u9toPKsvykpseJha2NddZQ4WcmbftC2M/DwPuinPsI8Mikba8SuX8GY8yHZltepRaS0GTLgRE/r9V28KFrVqT97P25ZMuyAh68dS3ffOkUt68v4V2bSqc/aQYGR/0cON/FR29I/nov09GZ/EotUKF8ZL8+1c6oPzAnhyenuwdvW8um8nz+8skjtPYO23KPt850MTZuuC7FfU+RaIBRaoHy5QTzkb10rIW8TBdXpzhv1XyU4XTwjQ9sYXhsnD9/4nDURdJm47WadtxOR1q+fxpglFqgfLke2vpHePl4K7dcVpyyNUPmuzVFXj5393p+ebKNH+05l/Dr/6a2na3LC8hyp1/2Bf2LUmqB8nndjPoDdAyMcseGuT08Od196JoV3Fjh4yv/e4y6BKb07x4cpbqxl+vWpNfw5BANMEotUIU5wbkwLodw8zqddGwnh0P4+n2bcbsc/OnjhxKWcfmNug6MgevXpl//C2iAUWrB8ln5qravWpx2w1vnoyX5mXz53koOne/mu7trE3LN39R0kO12snlZQUKul2gaYJRaoEryggFmrs/en0vevXkpO7Ys5ZsvneLQ+e5ZX++12na2r1qctv1n6VkqpZTtLivJ5VsPbOWDb1ue6qIsKA/fU0lxrof/76nZJRtp7hmmtm2A69O0/wU0wCi1YIkI92xeOqfWfpkP8rMz+KOb13C4vmdWaf1fqw3mObsuTftfQAOMUkol3Y4tS3G7HDyxr37G13ittoNF2RmsX5KXwJIllgYYpZRKsoJsN3duXMLPDzTMaN0YYwyv1bRz7ZpCHI70Te+jAUYppVLg/dvK6Rka44WjLXGfe7ZjkMaeYa5N4/4X0ACjlFIpcf0aH2UFWTy+9/z0B0/yG6v/5fo0zD8WTgOMUkqlgMMh3HdVOa/WtFPfNRjXua/VdFCan8kqX45NpUsMDTBKKZUi79tWDhBXZ38gYHi9roNr1xSm/fIKGmCUUipFyhdlc/0aHz/ZW08gEFum5ePNfXQOjKb1/JcQDTBKKZVC7796GQ3dQ7xe1xHT8XNh/kuIBhillEqhd2woIS/TxWNvxdbZ/1ptB6t9OZTmZ9lcstnTAKOUUimUmeHk3q1lPFvdTM/g2JTHjo0H2FPXMSdqL6ABRimlUu7925Yx6g/w1KGGKY87XN/NwOh42q7/MpkGGKWUSrHKsnw2lOZNOyfmtZoORODa1VqDUUopFaMPXL2MqoZeqhujJ8D8TW07G0rzWJTjTmLJZk4DjFJKpYEdW5bidjr4yd7Ic2KGRsfZf7ab69J89n44DTBKKZUGCrLdvGNjSdQEmPvOdjE6HuC6tXOj/wU0wCilVNr4wNXLoibA/E1tOy6HsH3l4hSUbGY0wCilVJqYKgHma7UdbFlWQI7HlYKSzYwGGKWUShPhCTAbuocmtvcMjXGkvntONY+BBhillEor911lJcAM6+zfU9dBwDCnOvhBA4xSSqWVZYutBJj7zk8kwHyttoPMDAdblxekuHTx0QCjlFJp5n3byqnvupAA87Xadq5euRiPy5niksVHA4xSSqWZOzcuIS/TxeN7z9PaN8zJlv45kx4m3NwZjqCUUgtEKAHmo2+dZ9uKRQBcP0cSXIbTGoxSSqWhUALMv3v+JHmZLjYuzU91keKmAUYppdJQKAFmz9AY16wuxOlI7+WRI9EAo5RSaer924JDlq+fY/NfQrQPRiml0tR925ZxrnOIezYvTXVRZkQDjFJKpSmvx8UX3r0h1cWYMW0iU0opZQtbA4yI3CUiJ0SkRkQeirDfIyKPWfv3iMjKsH2fs7afEJE7rW3LROQVETkmItUi8idhxy8WkRdE5JT130V2PptSSqmp2RZgRMQJfAe4G9gAPCAik+t6HwW6jDFrgW8AX7XO3QDcD2wE7gK+a13PD/yZMWY9cA3wibBrPgS8ZIypAF6yfldKKZUidtZgtgM1xpg6Y8wo8CiwY9IxO4AfWD8/AdwuImJtf9QYM2KMOQ3UANuNMU3GmP0Axpg+4BhQFuFaPwDutem5lFJKxcDOAFMGhC9qUM+FYHDJMcYYP9ADFMZyrtWcthXYY20qMcY0WddqAoojFUpEPi4ie0Vkb1tbW9wPpZRSKjZ2BphIs4JMjMdMea6IeIGfAp8yxvTGUyhjzPeMMduMMduKioriOVUppVQc7Aww9cCysN/LgcZox4iIC8gHOqc6V0QyCAaXHxtjfhZ2TIuIlFrHlAKtCXsSpZRScbMzwLwFVIjIKhFxE+y03znpmJ3Ah62f7wNeNsYYa/v91iizVUAF8KbVP/N94Jgx5h+muNaHgacS/kRKKaViJsHPc5suLvJO4B8BJ/DvxphHRORhYK8xZqeIZAL/RbAvpRO43xhTZ537l8AfEBw59iljzC4RuQH4NXAECFi3+bwx5hkRKQQeB5YD54D3GWM6pylfG3B2ho/nA9pneG660WdJP/PlOUCfJV3N5llWGGOm7WOwNcDMZyKy1xizLdXlSAR9lvQzX54D9FnSVTKeRWfyK6WUsoUGGKWUUrbQADNz30t1ARJInyX9zJfnAH2WdGX7s2gfjFJKKVtoDUYppZQtNMAopZSyhQaYacxmyYF0E8OzfERE2kTkoPX6WCrKOR0R+XcRaRWRqij7RUS+ZT3nYRG5MtlljEUMz3GLiPSEvR9fSHYZYzXVUhphx8yV9yWWZ0n790ZEMkXkTRE5ZD3HlyIcY+/nl9PxfXUAAAPhSURBVDFGX1FeBCeI1gKrATdwCNgw6Zj/C/yL9fP9wGOpLvcsnuUjwLdTXdYYnuUm4EqgKsr+dwK7COa0uwbYk+oyz/A5bgF+kepyxvgspcCV1s+5wMkIf19z5X2J5VnS/r2x/p291s8ZBBMDXzPpGFs/v7QGM7XZLDmQbmJ5ljnBGPMrgpkfotkB/NAEvQEUhPLUpZMYnmPOMFMvpREyV96XWJ4l7Vn/zv3WrxnWa/KoLls/vzTATG02Sw6km1ieBeC3reaLJ0RkWYT9c0GszzoXXGs1cewSkY2pLkwsIiylETLn3pcpngXmwHsjIk4ROUgw+e8Lxpio74kdn18aYKY2myUH0k0s5XwaWGmM2QS8yIVvNnPNXHlPprOfYM6nzcA/AU+muDzTmmYpjTn1vkzzLHPivTHGjBtjthDMSL9dRConHWLre6IBZmqzWXIg3Uz7LMaYDmPMiPXrvwJXJalsiRbL+5b2jDG9oSYOY8wzQIaI+FJcrKimWEojZM68L9M9y1x7b4wx3cBugkvQh7P180sDzNRms+RAupn2WSa1h99DsO15LtoJ/J41aukaoMdYq53OJSKyJNQeLiLbCf7/2pHaUkVmlTPaUhohc+J9ieVZ5sJ7IyJFIlJg/ZwF3AEcn3SYrZ9frkRdaD4yxvhF5EHgOS4sOVAtYUsOEPxD/C8RqcFaciB1JY4uxmf5pIjcQ3CJhE6Co8rSjoj8D8FRPD4RqQe+SLADE2PMvwDPEByxVAMMAr+fmpJOLYbnuA/4YxHxA0MEl7NIxy8vANcDHwKOWG3+AJ8nuHzGnHpfiO1Z5sJ7Uwr8QEScBAPg48aYXyTz80tTxSillLKFNpEppZSyhQYYpZRSttAAo5RSyhYaYJRSStlCA4xSSilbaIBRKsVE5C+tbLeHrcy8bxOR3RLMfH3Qyur78VSXU6l46TwYpVJIRK4Ffotg9t4Raza429r9QWPMXhFZDNSKyH9aiUqVmhM0wCiVWqVAeyhFjzGmHWBSQlsvMACMJ710Ss2CNpEplVrPA8tE5KSIfFdEbg7b92MROQycAP7GGKMBRs0pGmCUSiErYeJVwMeBNuAxEfmItfuDVmbr5cBnRGRFakqp1MxoE5lSKWbVTHYDu0XkCBeSD4b2t4nIfuBtwNnkl1CpmdEajFIpJCKXiUhF2KYtTAoiIpJNcNGr2mSWTanZ0hqMUqnlBf7JSqvuJ5hp+OMEl6/9sYgMAR7gP40x+1JXTKXip9mUlVJK2UKbyJRSStlCA4xSSilbaIBRSillCw0wSimlbKEBRimllC00wCillLKFBhillFK2+P8BAuxPTRs/KawAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/usr/lib64/python3.6/site-packages/matplotlib/cbook/deprecation.py:107: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n",
      "  warnings.warn(message, mplDeprecation, stacklevel=1)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAD8CAYAAABpcuN4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3Xl83Gd96PvPMzOa0TLSSNZIsjYvsmTHsuIlcUwWkpDdgRKHNkByWi63By4955LbhdKStPfSNqf0FGhLaQvti1Pa0lPaJARKHEicPYSQxIlXWbK8SLIta9/3bZbn/jG/kcfyzGj2xfq+Xy+/Iv3mtzwTy/rOs3y/j9JaI4QQQoRiSncDhBBCZDYJFEIIIcKSQCGEECIsCRRCCCHCkkAhhBAiLAkUQgghwpJAIYQQIiwJFEIIIcKSQCGEECIsS7obkAhOp1Nv2LAh3c0QQoiscvjw4WGtddlK510VgWLDhg0cOnQo3c0QQoisopS6EMl5MvQkhBAiLAkUQgghwpJAIYQQIiwJFEIIIcKSQCGEECIsCRRCCCHCkkAhhBAiLAkUEfrx0R7GZxfT3QwhhEg5CRQR6Bia5refOsYPj/SkuylCCJFyEigicPziOAADk/NpbokQQqSeBIoINHdPANA/IYFCCLH6SKCIQHO39CiEEKuXBIoVuDxeWnsnARicWkhza4QQIvUkUKzgzMAUC24v5YU2Bibn0Vqnu0lCCJFSEihWcMKYn7i7sYLZRQ/TC+40t0gIIVJLAsUKjndPUJRrYc+GNYDMUwghVh8JFCto7h5ne00xFUW5AAxMyjyFEGJ1kUARxrzLw+n+KbbXOFjr8AcK6VEIIVYXCRRhtPVN4vZqttcUU15oA6RHIYRYfSRQhOFPtNte46DAZqHQZpEehRBi1ZFAEcbx7nGcdhuVxrBTeZGNwSkJFEKI1SWiQKGU2quUOq2UaldKPRbkdZtS6inj9YNKqQ0Brz1uHD+tlLpv2XVmpdRRpdRPAo5tNO5x1rinNfa3F58T3RPsqHGglAKgoihXyngIIVadFQOFUsoMfAu4H2gEHlFKNS477TPAmNa6HvgG8FXj2kbgYWAbsBf4tnE/v98C2pbd66vAN7TWDcCYce+Um15w0z40zbU1jqVjFUW5MkchhFh1IulR7AHatdadWutF4Elg37Jz9gHfM75+BrhL+T6G7wOe1FovaK3PAe3G/VBK1QAfAf7RfxPjmjuNe2Dc88FY3li8Wnom0Bp21BQvHasoymVwSrKzhRCrSySBohq4GPB9t3Es6DlaazcwAZSucO1fA78PeANeLwXGjXuEehYASqnPKaUOKaUODQ0NRfA2onMiYCLbr6LIhsujGZt1Jfx5QgiRqSIJFCrIseUfqUOdE/S4UuqXgEGt9eEYnuU7qPV3tNa7tda7y8rKgp0Sl+Pd41QX51Fqty0du5R0J/MUQojVI5JA0Q3UBnxfA/SGOkcpZQEcwGiYa28BHlBKncc3lHWnUurfgGGg2LhHqGelRHP3xGW9CfD1KEAChRBidYkkULwPNBirkaz4Jqf3LztnP/Bp4+uHgNe0byB/P/CwsSpqI9AAvKe1flxrXaO13mDc7zWt9a8Z17xu3APjns/G8f5iMj67SNfoLNsD5icAygulRyGEWH1WDBTGfMGjwIv4Vig9rbVuVUo9oZR6wDjtu0CpUqod+ALwmHFtK/A0cBI4AHxea+1Z4ZFfAr5g3KvUuHdK+RPtdizrUZQXSXa2EGL1sax8CmitnweeX3bsywFfzwMfD3HtV4CvhLn3G8AbAd93YqyMShf/jnbbqi8PFDaLmTUF1pT3KGYX3UzPuyk35kiEECKVJDM7iObuCeqcBTjycq54zbeBUWp7FH+8v5Vf+tu3cHu8K58shBAJJoEiiGAT2X7+XIpUWXB7eKGln8GpBQ5fGEvZc4UQwk8CxTKDk/P0T85z7bKJbL+KIltKy3i83T7C1LwvreTlkwMpe64QQvhJoFjmeIiJbL+KolyGpxdSNgz0QksfhTYLN9WV8nLbgGSFCyFSTgLFMie6xzEp2FYVOlB4NYzMLCa9LW6Pl5dPDnDn1nI+sr2SCyOznBmYTsi9Xzs1wN+/0ZGQewkhrm4SKJY53j3B5opC8qzmoK+nMjv74LlRxmZd3N+0lnsaKwB4+WR/3PfVWvPnL5ziGy+fYdEtE+RCiPAkUATQWht7ZAfvTUBgdnbyVz690NJHXo6Z2zeXU1GUy47a4oTMU7T1TXFmYJpFj5ezg1MJaKkQ4momgSJA99gcY7OuKzKyA6WqR+H1al5sHeBDW8qWejf3NlZwvHsi7mc/e6xn6evWnsm47iWEuPpJoAhwKSM7dKAoLbBiUskPFIe7xhiaWmBv09qlY5eGn2LvVXi8mmeP9XLXNeXYbRZaeifibqsQ4uomgSJAc/c4VrOJLWsLQ55jMZtw2m1JDxQvnOjHajZx5zXlS8cayu2sL82PK1AcPDdC/+Q8D+6qprGqiJYeCRRCiPAkUARo7p5ga2UhVkv4/y1rHcnd6U5rzYut/dza4KQw91J2uFKKexsreKdjhKn52PbEePZoL3abhbu3VtBU5eBk3yQeryy5FUKEJoHC4PVqWnomws5P+JUX5ia1R9HcPUHP+Nxlw05+9zSuZdHj5Wdnot+sad7l4fkTfdy3bS15VjNN1UXMu7x0DiVmyW2kJuZc3Pq113gzhvcghEg9CRSGzuEZphbcl+2RHUpFkY3BqeT1KF5o6cdiUktzEoGuX1/CmgJrTMNPr58aZGrBzYO7qgBoMooepnqe4mTvJBdH5/hfP+9M6XOFELGRQGE40eOrGBtuItuvoiiX0ZlFFtwrVUyPntaaAy193LSplOJ86xWvm02KO68p5/VTg7iizA7/8bEeygpt3LzJCUCds4DcHBMtKV751GH0YH5+dpiLo7MpfbYQInoSKAzHL06Ql2Omvty+4rn+XIrBJMxTnOqf4vzIbNBhJ797GiuYnHfz3rnRiO87Mevi9VNDPLCjCrPJt+OsxWxia2XqJ7Q7h2awmk2YFDx96OLKFwgh0koChaG5e5ym6qKlX6Lh+HMpklFF9oWWfpSCextDB4pbG5zYLKaohp+eb+lj0ePlwZ3Vlx3fVlXEyd5JvCmc0O4Ymqa+3M7tm8t4+tBFKZ8uRIaTQAG4PF5aeycjmsiGwKS7xPcoDrT0ccOGNZQV2kKek2+1cGtDGS+fjLxI4I+P9lBXVkBTddFlx5uqHEwtuOlK4RBQ5/A0dWUFfPKGdQxMLvDGaZnUFiKTSaAAzg5Ms+D2hi3dEShZ2dkdQ9OcGZjm/jDDTn73NlbQMz7Hyb6V5xd6xuc4eG6Uj+2sRqnLe0ypntCed3noHptjU5mdu7aW47TbePL9rpQ8WwgRGwkUXNr6NNIeRUl+DjlmlfAexYEWX8G/cPMTfnduLUcpeKl15eGn/cd6Adi3bNgJoKHCTo5ZpWxC+/zIDFpDXVkBOWYTH99dw2unBlO6x4cQIjoSKPBVjC3KtbChND+i85VSScmleKGlj521xVQ68lY812m3cf26kojmKZ491sN164pZF+T92SxmNlcU0pqiHkXn0AwAm8p8iwYevqEWr4YfyKS2EBlLAgW+pbHba4qvGJYJp6IosWU8Lo7O0tIzGdGwk989jRWc7Jukeyz0/EJb3ySn+qf42K4rexN+TVUOWnomUrIpkj+5r66sAID1pQXcvKmUpw5dTOmEuhAicqs+UMy7PJzqm4oo0S6Qr4xH4gKFf9jp/qbKiK/xJ+S9EqZX8eNjPVhMio9srwp5TlN1EWOzLnpTMPzTMTRDpSOXfKtl6djDe9bRPTbHW+3DSX++ECJ6qz5QtPVN4vbqkFufhlJemJvQPIoXWvporCwKOjwUSl2ZnfpyOy+3BQ8UXq9m/7FebttcxpqCK5P3/Lb5J7RTkE/ROTS9NOzkd9+2Ckryc2RSW4gMteoDxQnjl2OkE9l+FUW5TC24mVlwx92G/ol5jnSNRzXs5HdPYwUHO0eZmLuySODBc6P0TfgqxYazdW0RJgWtSQ4UWms6hmaWhp38bBYzv3xdDS+fHGB4OvkbQgkhorPqA8XxixM47TYqHblRXbeUnZ2Amk8vthrDTtfGFijcXs0bpweveO3ZYz0UWM3cs/XKmlGB8qy+jPSW3uSufBqaWmB6wX1FjwLgkT21uDyaHx7uTmobhBDRW/WBwr/1aTQT2XAplyIRyzpfaOmjvtxOfXnofTBC2VlTTFmh7YplsvMuDz8NqBS7Ev+EdjJ1GCuelvcoAOrLC9m9voSn3r+Ykkl1IUTkVnWgmFlw0z40HXGiXaBLPYr4AsXI9ALvnRuNadgJwGRS3L21nDdOD15WpPCN04NMzbtXHHby21btYHBqgcEklk/vWFrxFLye1sN71tE5PMPBKGpYCSGSb1UHCt+S0Mgqxi6XqOzsl04O4NWRJdmFcm/jWmYWPbzTMbJ07MdHe3Habdy8qTSiezRV+Up7tCZx+KlzaIa8HDOVRcGH+T5ybSWFuRaefE8mtYXIJKs6UPj3yI52aSyA3WYh32qOOzv7hZZ+1q3Jp7GyaOWTQ7hpUyn5VvNS8t3EnIvXTg3y0R2VWMyR/RU3LgWK5A0/dQ5Ps9FZgClE4cU8q5kHd1bzfEs/47OLSWtHtugdn+N0/1S6myHE6g4UF8dmqS7Ow2kPXYAvFKUUFUXx5VJMzLp4u32Y+5vWRj1HEig3x8ztm8t4pW0Ar1fzwonglWLDKczNYaOzIKmlPDqGpoPOTwR6eE8ti24v/3m0J2ntyBZf/MFxPv/vR9LdDCFWd6B4Yl8TL3/htpivLy+0xZVL8UrbAG6vjmvYye+exgoGJhc40TPBj4/1UOcsiHruZVtVUdKKAwYWAwzfBgfbaxw8+d7qntQenVnk3c4RqYElMsKqDhTAZRnC0aooyqU/jh7FCy39VDpyY5ojWe7Oa8oxmxTfe+c8B8+Nsi9IpdiVNFU76B6bS8qwT2AxwJV88oZaTg9McfTieMLbkS1eMeauphfczLsSv5OiENFY9YEiHv4yHrF88p1ddPPm2SHu27Y25Jh9NIrzrdywoYQfHelBa9i3M3TJjlCaqnw9kGRMaC8vBhjOAzuqyMsx89R7q7dQ4AEjtwZ8+SdCpJMEijiUF9pYcHuZnIs+O/v4xQkW3V5u31yWsPb4d8Xbta6YDc6VP7kvt82Y0E5GPsXyYoDhFObm8NEdlTzX3Mt0AjLfs83UvIu3zg4v/b8amZGJfZFeEijisLRENoZciqMXxwDYWRv/sJPfvdsqsJpNfPz62piuLymwUl2cl5QM7WDFAMN5eM86Zhc9S3tprCavnRpk0ePlVz+wHoBh6VGINJNAEYd4cimOdo1T5yygJEyxvmjVlOTz1pfu4JE9sQUK8FWSTUbNp2DFAMPZVVvMlorCVVko8MXWfsoKbdxrVAeW+lci3SRQxMGfnR3tyhStNUe7xtm5LnG9Cb/yoty4lto2VTnoHJ5hav7KIoOx0lrTGaQYYDhKKR7eU0tz90TKNlXKBPMuD6+fGuK+bRVL+6bL0JNINwkUcfD3KKItDNg9Nsfw9AK71pUko1lx8e+h3daXuESvoakFpkIUAwznY7uqsVpMPLmKJrXfPDPEnMvD3m2V5OaYKbRZZDJbpF1EgUIptVcpdVop1a6UeizI6zal1FPG6weVUhsCXnvcOH5aKXWfcSxXKfWeUuq4UqpVKfUnAef/i1LqnFLqmPFnZ/xvMzlyc8w48nKiHnryL/vclcD5iUTZVp34Ce1wxQDDKc63cm9jBS+09K2anIoDrf048nL4QN0aAErtVhl6Emm3YqBQSpmBbwH3A43AI0qpxmWnfQYY01rXA98Avmpc2wg8DGwD9gLfNu63ANyptd4B7AT2KqVuDLjf72mtdxp/jsX1DpMsli1Rj3aNkZtj4pq10VeLTbbywlzKC20JTbxbqRhgOHs2rmF4epGe8bmEtSdTuTxeXjk5wN1bK8gxSq847TYJFCLtIlmCsgdo11p3AiilngT2AScDztkH/LHx9TPA3ynfQPk+4Emt9QJwTinVDuzRWr8DTBvn5xh/svIjo6+MR3T/kI92+fbojrQOU6o1VTtoTWApj5WKAYbjT0Zs7p6gpiTy3f8SQWvNxJyLiTkXU/NupubdTC+4mV649L3vmO/7a6sdfPbWupif927nCJPz7ssy9Z1221KgFSJdIgkU1UDgIHE38IFQ52it3UqpCaDUOP7usmurYamnchioB76ltT4YcN5XlFJfBl4FHjMCzWWUUp8DPgewbt26CN5GcpQX5tIxGPlezwtuDyd7J/n1D25IXqPi1FRVxBunB5lb9ES0l8VKVioGGM41lYXkmBXHu8f58LWR7yceCa01Y7Muusdm6R6bC/jvHD3G9zOL4bOirWYThbkWlIL9x3u5fXMZDRWx9RRfaOkn32rm1gbn0rFSu5WD56RHIdIrkkAR7F/38k//oc4Jea3W2gPsVEoVA/+plGrSWrcAjwP9gBX4DvAl4IkrbqL1d4zX2b17d9p6IxVFNganFvB6dUS/CFt7J1n0eNlVm3kT2X7bqh14NZzqn0zIhHvH0HTMZUpsFjNbK4tovpi4obAfH+3h22+00z02x+yyQFCUa6GmJJ/1pfncUu+kqjiXknwrhbkW7LkWCm05l77OtWCz+ALp6Mwit371Nf7mtXb+9pFdUbfJ49W81DrAHdeUk5tzKTg77TbGZl24PN6l4SghUi2SQNENBC7MrwGWZ0H5z+lWSlkABzAaybVa63Gl1Bv45jBatNZ9xksLSql/Br4Y2VtJj4qiXNxezcjM4tJyxnCOdhkT2UlYGpso/pVPLb3xBwp/McBf3lUT8z221zh49mhvxMF4Jf/y9nkm59x88oZaakryqSnJM/7k48jLiemeawqsfPrmDfz9zzr4zTvro+5VHOkaY3h6gb3bLi8Q6TR+psZmFimPYehOiESI5CPK+0CDUmqjUsqKb3J6/7Jz9gOfNr5+CHhN+5ap7AceNlZFbQQagPeUUmVGTwKlVB5wN3DK+L7S+K8CHgRa4nmDyRZt0t3RrjGqi/OWrstEVY5cSvJzEpJ4d2FkNuJigKFsrylmasFN5/BM3O2Zd3lo7Z1g364q/uij2/jMBzdy37a1bKtyxBwk/D57ax35OWb+5rX2qK890NKP1WzijmvKLzvuNBIyh2RCW6TRioFCa+0GHgVeBNqAp7XWrUqpJ5RSDxinfRcoNSarvwA8ZlzbCjyNb+L7APB5Y8ipEnhdKdWMLxC9rLX+iXGv7yulTgAnACfwp4l5q8kR7ZaoyUq0SySlFE3VjoSsfPJPxEabQxHIP2x1PAHVZJu7J3B5NLvXr4n7Xsv5exU/ae7l7EDkeShaaw609HNrgxO77fJOvr9HMTwtSXcifSIa9NRaP6+13qy13qS1/opx7Mta6/3G1/Na649rreu11nv8K6SM175iXLdFa/2CcaxZa71La71da92ktX4i4Pw7tdbXGsd/TWud0Us+LvUoVv7ENzg5T8/4XEbmTyy3rcrB6f4pFt3euO4TTTHAUOrL7eRbzTR3xx8oDl/w1di6LknBOpZeRWvvJD3jc9wXZF8S/6ZaI9KjEGkks2Nx8s9LRDL0tJRol4EZ2cs1VRfh8mjORPHJOJhoiwEGYzb5ejjHu+Pv4Ry+MEqds4DSGHY1jEQsvYoXWvowmxT3bK244jWn3Tf0JLkUIp0kUMQpx2zCabdGFii6xskxq6Vy3pns0t4U8f1yjrYYYCg7ahyc7JuMq4ejtebwhTGuW5/cQB1tr+JASz831q0JWiDSbrNgtZhk6EmklQSKBCgvjCzp7tjFMRqrHJctf8xU69bkY7dZ4tpDO5ZigKFsrylm0e2Nq4fTOTzD2KyL3UkOFNH0KtoHp+gYmrlitZOfUooyu01KjYu0kkCRAP6d7sJxe7w0d09kxfwEgMmkaIxzD+1YiwEGszShHcc8hX9+4vokBwqIvFdxoMW3k929IQIF+IafhqWCrEgjCRQJ4Kv3FP4T35mBaWYXPRmdP7FcU5WDtr5J3J7YhntiLQYYTO2aPEryc+JKvDtyYQxHXk5CAtdKIu1VHGjt57p1xWGXS5dKj0KkmQSKBCgvzGVkZgFXmF+o/h3tMjkje7mm6iLmXd6Y8xc6h2MvBricUopra4rj6lEcujDGdeuKE5K0F4mVehUXR2dp6Zm8rLZTME6pICvSTAJFAlQU5aI1YfcNONo1TmmBldo1eSlsWXyWMrRjTLzrGIy9GGAwO2ocnBmYYnYx+n20x2cXaR+cZveGxOdPhLJSr+LFVt+w031hhp3At0R2dGYRrzcr62aKq4AEigTwJ92Fm6c42jXGrnXFce0+l2p1zgJyc0wxT2jHUwwwmB01xXi1L+8gWke6/PkTqe3RhetVHGjpZ2tlEetLww/NldptuL2+SrZCpIMEigRYKeluYtZFx9BMVuRPBLKYTWytjH1Cu2NoOiHzE37ba309nFgytA9fGMNsUuww7pEqoXoVg5PzHO4aC7naKZDkUoh0k0CRAJe2RA3eozjWnbk72q2kqcrByd7JqIc9/MUAEzlxXF6YS6Ujl+YYEu8OnR9jW1VRXIl/sQrWq3jp5ABaw/3XrhwoyuxSxkOklwSKBCgtsGI2qZBDT0e7xlAKtmdjoKguYnrBzYXR2aiuS0QxwGC21ziiLuXh8ng53j2e8mEnv2C9ihdb+6lzFtBQvnIgvVTvSXoUIj0kUCSAyaQoLwy9RPZo1zhbKgqvKPiWDbYZGdr+HIRIJaIYYDDba4o5PzLL+Gzkn67b+iaZd3nZvSF9Q3+BvYrx2UXe6Rjhvqa1Ec1ZlRbI0JNILwkUCVJeFDzpzuvVHLs4nlX5E4G2VhZR5yzgu2+di2r4KRHFAIMJ3Bo1UofOpy7RLpTAXsU//KwTt1dHND8BUJLv67FKoBDpIoEiQSoKbUEDxbmRGSbmXFmVPxHIbFI8emc9bX2TvNw2EPF1nQkoBhjMtTW+Hk40w0+HjT1AKh3pXZrs71X8w886qHTksr0msol1k0mxpsDKiMxRiDSRQJEgvjIeV37iy4Yd7VbywI4qNpTm8zevnsW3H9XKOhJUDHA5R14Odc6CiCvJaq05fD75hQAj4e9VgC93Ipql0qUFknQn0kcCRYJUFOUyMedi3nX5HsxHu8YotFlSUjYiWSxmE4/e2UBr7ySvtA2ueH4iiwEGE82Edu/EPP2T80kvBBip/+vWOu68ppxf/cC6qK4rK7QxJD0KkSYSKBKk3FiZMrisV+Hf0S5VZSOS5cGdVaxbk883Xz2zYq8ikcUAg9leU8zA5EJEpd0PnR8F0js/EaikwMo//Z83RL2nttNuk82LRNpIoEiQpaS7gFyK2UU3p/onszJ/YjmL2cSjd9TT0jPJa6fC9yoSWQwwmB1RJN4duTBGvtXMNWuj+8Wcafz1niId+hMikSRQJIg/UPRPXAoUJ7on8Ors2NEuEh+7rpraNXkrzlUkshhgMI2VDswmFdHKp0MXxthZW4zFnN0/6qV2G/MuLzOLnpVPXkWm5l187l8P8f/++AQ/PNzNueEZCaZJkH0L+zNUsHpP/q1Pd14FPQrw7eb3+Q/V89iPTvDGmSHu2FIe9LxEFwNcLs9qZnNF4YqVZGcW3LT1TfLoHfVJaUcqBe6dnY35OMlypGucl04OYLWY+Ld3uwAozs9hV20xu9aVcN26ErbXOijKzUlzS7Ob/MQliCMvB5vFxGBABdmjXWNsdBYE3eIyW/3ydTX87WvtfPOVs3xoc1nQlTuJLgYYzI4aBy+09KO1Drl66PjFcbyajFjxFK/Aek8rFRFcTfz5Om/9/h2Mzbo42jXG0a5xjnSN8frpIQCUgoZyO7tqS3hwVzU3bSpNZ5OzkgSKBFFKURGQdKe15kjXOLfWO9PcssSyWkz833ds4g//s4U3zw5z++ayK87pHJqJOEcgVjtqi3ny/YtcGJllgzP4L85DF3ylU66GoT9/j2JoSlY+BTo3PENhroWyQhvlRblsWVvIw3t8K8om5lw0d49ztGuco11j/PREH++fH+W1L34ovY3OQtk9cJthfDvd+QJF78Q8Q1MLWZ0/EcrHr6+lypHLN1+5cgXUvMvDxbHZpC8H9geicMNPhy+Msbm8EEde9g87LA09zcjKp0CdQzPUOQuC9iodeTnc2lDGb97VwD//+h5+/ZYNnB+ZYcEt8zzRkkCRQL4yHr5/yEeN/Q+uhk+zy1ktJv77HfUc6Rrnrfbhy15LVjHA5TZXFGKzmEJOaHu9miNdmZFolwhr/PWepEdxmXPDMxEvmqgvt+PVcH44ugKXQgJFQlUU+oaetNYc7RonN8fElixflhnKJ3bXUOnI5ZuvXL4CqjNJxQCXyzGb2FZVFDLx7uzgNFPz7oxJtIuX1WLCkZcj2dkB5hY99IzPsTHE0ONy/p/J9sHpZDbrqiSBIoHWOmzMLnqYXnBztGuM7dXF5GT5ssxQbBYz//1Dmzh0YYx3OkaWjnckqRhgMNtrimnpmcQdZK/yQxcyK9EuEWTv7MudH4kuX2dTmR2l4OzgldvSivCuzt9iaeLPpbg4OkdL7+RVOT8R6BO7a6kosvHXr55dOpasYoDB7Kh1MOfy0D505SfEwxfGcNqtrC/NT3o7UsWXnS1DT36dRmJnpD2KPKuZ6uI86VHEQAJFApUX+gLFG2cGWXR7r5r8iVByc8z8t9s38d650aVeRbKKAQaz3Sg5HixD+8iFMa5bV5JVe5SvxFlokx5FgHNGYmekgQJ8y2QlUERPAkUC+ZPuXmzpB67OiezlHtmzjrJC21INqGQWA1xuY2kBhbmWKyrJDk0tcH5kNq0bFSWDs8DKkASKJbH0XuvL7XQOz+CJcmvf1U4CRQKVG0NPx7snqHTkstaRnMzkTOLvVbzbOcpPT/QltRjgciaTClpJ9khX+jcqSgan3cbUvFuWdxo6h6P/UFJfbmfR7aV7TFY+RUMCRQLZbZal8gpX+/xEoP+yZx1Ou40/erYVSM1Ett/2mmJO9U1dVt798IUxrGYTTdXJTfpLNf/e2TJP4S9Fsh4nAAAd6klEQVRlPx3VsBP4AgXA2QEZfoqGBIoE8w8/ZeuOdrHIs5r5jdvqGJnx/QJLVjHAYHbUOHB7NW19k0vHDl8Y49oaBzaLOWXtSAXZO/uS0ZlFJufd1Dmj+1mrL/MtVw+2AEKEJoEiwfwrn1ZTjwLgV29cR2mBNanFAIPZvmwP7XmXhxPdE1fdsBNIjyJQ57Cx4inK3qsjP4eyQptMaEdJaj0lWEVRLhaTuuqGPVaSb7XwlY810TE0k9JNmioduTjttqVSHq29Eyx6vFdloCjz13uSHgXnjKWxm6LsUQDUl8nKp2hJoEiwT920nuvXl5Cbc3UNe0Rib1Nlyp+plGJnrWNpieyh876J7OuuwhVnpfbUDj1Nzrsytjx3x/A0VrOJ6pK8qK+tL7fz46M9YSsPi8vJ0FOCXbeuhF+7cX26m7GqbK8ppnN4hql5F4cvjLGhNJ8yY5jmapJvtZBvNaek3tPZgSmue+JlXj+98h7p6XBuaIb1pfmYY+i91pfbmVpwL9VlEyuTQCGy3vYaB1rDiZ6Jq6oQYDBOuy0lFWR/dLQHt1fzUmt/0p8Vi87hmahXPPn5Vz7J8FPkIgoUSqm9SqnTSql2pdRjQV63KaWeMl4/qJTaEPDa48bx00qp+4xjuUqp95RSx5VSrUqpPwk4f6Nxj7PGPa+eXX9EUvgntJ873sfw9CK7169Jc4uSJxX1nrTWPHe8F4A3zwxn3NaiHq/mwkjkVWOXa1gKFFLzKVIrBgqllBn4FnA/0Ag8opRqXHbaZ4AxrXU98A3gq8a1jcDDwDZgL/Bt434LwJ1a6x3ATmCvUupG415fBb6htW4Axox7CxHSmgIrtWvy+NGRbuDqS7QLVGq3JX3o6djFcbrH5rh+fQk943N0GBPHmaJ7bBaXR1MXY4+irNBGYa5FlshGIZIexR6gXWvdqbVeBJ4E9i07Zx/wPePrZ4C7lG+WaB/wpNZ6QWt9DmgH9mgf/99SjvFHG9fcadwD454PxvjexCqyvaaYBbeXwlzL0ifGq1Eqhp6eO96H1Wzif+xrAuDNM0NJfV60/EtjY03sVEpRLzWfohJJoKgGLgZ8320cC3qO1toNTACl4a5VSpmVUseAQeBlrfVB45px4x6hniXEFXYYO95dt64kpctzU63MbmV0ZjFptYo8Xs1Pmnv50JYyGquKqHMW8POzGRYooqwaG4wskY1OJIEi2L+65T+loc4Jea3W2qO13gnUAHuUUk0RPsv3QKU+p5Q6pJQ6NDSUWT/IIvV2GPMUV8tGRaGU2m14tS8zORnePz/K4NQCH91RBcCtDU7e7RzNqPpS54anceTlLO36F4v6cjvD04uMz0ryYiQiCRTdQG3A9zVAb6hzlFIWwAGMRnKt1noceAPfHMYwUGzcI9Sz/Nd9R2u9W2u9u6ysLIK3Ia5mu9aV8Ou3bOCXr69Jd1OSKtl7Z+8/3ku+1cxdW8sBuG1zGXMuz1J+SiboHPKteIonB6KhQlY+RSOSQPE+0GCsRrLim5zev+yc/cCnja8fAl7TvqUS+4GHjVVRG4EG4D2lVJlSqhhAKZUH3A2cMq553bgHxj2fjf3tidXCajHxRx/dRnVx9AlY2cRpT97e2S6PlxdO9HH31oql0t031pWSY1a8mUHDT+diqBq73FLNJwkUEVkxUBjzBY8CLwJtwNNa61al1BNKqQeM074LlCql2oEvAI8Z17YCTwMngQPA57XWHqASeF0p1YwvEL2stf6Jca8vAV8w7lVq3FsIgW/oCZKTnf2L9mHGZl1Lw04ABTYL168v4c0zwwl/XixmF930TczHXcq+uiQPm8XEWQkUEYmohIfW+nng+WXHvhzw9Tzw8RDXfgX4yrJjzcCuEOd34ltpJYRYpiyJgeK5430U5lq4bbPzsuO3bS7jawdOMzg1v7SLY7qcG45/IhvAbFLUyYR2xCQzW4gsUpRnwWo2MZzgCrLzLg8vtfazd9vaK8qz39bgmwN862z6exX+FU+J2PNEtkWNnAQKIbKIUorSJGRnv3F6iKkF92XDTn6NlUWUFlgzIp/i3PAMSsGG0vgDRX25nZ7xOWYX3SufvMpJoBAiyyQjUDzX3EtpgZWbN5Ve8ZrJpLi1wcnPzw7jTfNe051D01Q58hJSndlf86kzwzLPM5EECiGyjNNuS+jmRTMLbl5tG+DD11ZiMQf/lXDb5jJGZhY5GbCTYDokYsWT39K2qFLzaUUSKITIMk67LaE9ilfaBph3eYMOO/l9sME3wZ3OZbK+fbJnYq7xtNyG0gLMJiXzFBGQQCFElim1WxmZXkxYVdfnjvextig3bFZ7eWEuWyuL0jpPMTy9yNSCO+4VT35Wi4n1pfkSKCIggUKILFNmt7Ho8TI5H/8k7MSsi5+dGeSXtleuWCPrts1ODl8YY2YhPZO/nUa111jLiwcjNZ8iI4FCiCzjTGAuxYut/bg8mgd2hh528ru9oQyXR/Nu50jcz41FonIoAtWX27kwMsui25uwe16NJFAIkWWW9s6eij9QPNfcy/rSfK6tdqx47vUbSsjLMadt+KlzeAarxZTQMi315XbcxkZIIjQJFEJkmUs9ivhWPg1PL/CL9mE+ur0qogJ7NouZG+vW8GaaEu86h2bYWFqQ0DLysi1qZCRQCJFlElVB9oUTfXg1YVc7LXfb5jLODc9wcXQ2rmfHonN4OmFLY/38NaMkUIQngUKILLOmwIpS8Q89PXe8j80VdrasLYz4mluNch6pXibr8njpGplN6PwE+IoeVhfnybaoK5BAIUSWMZsUa/KtDMUx9NQ7Psd750f56PbIexMAm8oKqC7OS/k8RffYHG6vTuiKJ79NUvNpRRIohMhCvuzs2HsUP23uA6IbdgJfranbNjt5u30Elyd1K4X8S2MT3aMA3xLZjqHptJcnyWQSKITIQs7C+Oo9Pdfcy/YaBxti+MV7W0MZUwtujl8cj/n50fIvjd2U4DkK8E1oz7u89IzPJfzeVwsJFEJkodICW8yrns4Pz9DcPRH1sJPfzZucmBQpHX7qGJqhJD+H4vzY98kORbZFXZkECiGyUDxDTz9p9m1D/5HtlTFd78jPYWdtMT9L4TLZc8PTSZmfAN/QE0igCEcChRBZyFloZWbRw9yiJ+pr9x/v5YYNJVTFkbh22+YymrvHGZ9N/N7dwXQOzSRlfgKgpMBKaYFVAkUYEiiEyELOgtjKeJzun+LMwDQPRDmJvdytDWVoDW+1J79XMb3gZnBqIeE5FIE2ldul3HgYEiiEyELOQt9Y/VCUgeK5472YFNx/bWzDTn47ahwU5VpSMk9xzr/9aZJ6FOCb0G4fnE5YRd6rjQQKIbLQUnZ2lBPaL7b2c2Nd6dL1sbKYTXzQ2PUu2b9cO4cTXzV2uYZyO5Pz7qgD72ohgUKILBRLBdkLIzOcHZzm7q0VCWnDrQ1l9E3MJ31sv3PIt0/2ujX5SXuG1HwKTwKFEFloTUH0FWRfaRsESFiguG2zr5zHz5I8/HRueIaaksTskx1KNgYKr1fzzOFuFtzRL2iIlgQKIbJQbo6ZwlwLIzORDz29cnKAzRV21pUm5pN5dXEem8oKkl5NtnN4mjpn8oadANYW5WK3WbIqUDx7vIcv/uA4r5wcTPqzJFAIkaXK7LaIx9QnZl28d340Yb0Jv9s2l3Gwc4R5V3I+1WqtOZfEpbF+Sik2lRVkTaCYd3n4ixfP0FRdxP1Na5P+PAkUQmSpUrs14qGnN84M4vFq7m5McKBoKGPB7eX986MJva/f4NQCM4uepJTuWK6+vDBrAsX33j5Pz/gcf3D/1oTuzxGKBAohspTTbot46OmVtkGcdis7a4oT2oYP1K3BajZxoKWfgcl5xmcXmXd5ErYSqnPIv/1pcoeewDdPMTi1wMScK+nPisfYzCJ/93o7d2wp4+Z6Z0qeaUnJU4QQCee023gngv2rXR4vb5we5P6mtQn/9JlvtfCBujV8/2AX3z/YddlrNosJm8VEbo7Z+GOiON/KE/u2cc3aoojuf2lpbCp6FJcmtK9fX5L058XqW6+3M7Pg5rH7t6bsmRIohMhSTruN8VkXLo+XHHPowYH3z40yNe/mrgTPT/j9+a9s592OERbcXuZdHubdHuZdXhbcHhZcvmP+1w6eG+XRfz/Kc49+kDzryquYzg3NkJtjYm1RblLaHsgfKDoyOFBcHJ3lX9+5wEPX10S14VS8JFAIkaVK7b4lsiPTi6x1hP5F+krbIFaLiVsbkjNMUV2cx69cXxPRuW+dHebXvnuQP3u+jf/xYNOK53cOz7DRaU/JOHxtSR5Wsymjd7v7+ounMZngC/dsSelzZY5CiCwVSdKd1ppX2ga4ZVMp+db0fy78YIOTz35wI//73Qu82jaw4vnnhmeSWrojkMVsoi6DVz41d4+z/3gvn/1gXdgPBskggUKILFVm1HsKFyjaB6fpGp1N+GqnePze3i1cs7aQ33+mmaEwq7YW3V66RmdTMj/hl6nbomqt+bPn2ygtsPIbt9el/PkSKITIUqVLFWRDr3x62fjUftc1mRMobBYzf/PILqYX3PzeM8dDrpC6ODaLx6uTnkMRqL7MzsWx2aTlhcTq9dODvNs5ym/d3UBhbk7Kny+BQogs5Sz0FwYM/an81bZBrq12pHyoYiWbKwr5gw9v5Y3TQ/zrOxeCnuNfGpvMYoDL1Zfb0Ro6Mmiewu3x8j+fP8VGZwGP7FmXljZIoBAiSxVYfUtOQw09DU8vcKRrLOHZ2Inyf9y0nju2lPFnz7dxZuDKvSDOGUtjU9mjyMRtUX94pJuzg9P8/n1bwq5uSyYJFEJkKaVU2L2zXzs1iNZw19byFLcsMkopvvbQDuw2C7/5H0evKG7XOTSD027FkZe6oZaNzgJMyrdENhPMLrr5y5fOcN26YvamoFRHKBIohMhizkJbyB7Fq20DVDpy2VYVWXJbOpQV2vjaQ9s51T/F1w+cvuw139LY1PUmwDd/sm5NfsYskf3uz88xOLXAH3x4K0olf4lwKBIohMhiZXZr0B7FvMvDm2eGuWtreVp/wUTirq0VfOrG9fzjW+f4+dlLJcs7h2aSXjU2mPpyO2cH0h8ohqcX+IefdXDftgp2b1iT1rZEFCiUUnuVUqeVUu1KqceCvG5TSj1lvH5QKbUh4LXHjeOnlVL3GcdqlVKvK6XalFKtSqnfCjj/j5VSPUqpY8afD8f/NoW4OjntwXsU73SOMOfyZOz8xHJ/8OGtbCor4HefPs7YzCKT8y6GpxfYmMKlsX6byu2cH5nB7fGm/NmBvvnKWebdXr6095q0tgMiCBRKKTPwLeB+oBF4RCnVuOy0zwBjWut64BvAV41rG4GHgW3AXuDbxv3cwO9qrbcCNwKfX3bPb2itdxp/no/rHQpxFSu1WxmdWcTrvXyJ6SsnB8i3mrmxrjRNLYtOntXMNx/exdjsIo/9qDkl+2SH0lBeiMujuTA6m/Jn+3UMTfPv73XxX/asS+mqr1Ai6VHsAdq11p1a60XgSWDfsnP2Ad8zvn4GuEv5+rv7gCe11gta63NAO7BHa92ntT4CoLWeAtqA6vjfjhCri9Nuw+PVjAdUPNVa82rbILc1lCV1V7hEa6p28Hv3beHF1gG+/qJvviKVyXZ+1xg1lJ5v7kv5s/2+fuA0uRYTv3lXQ9raECiSQFENXAz4vpsrf6kvnaO1dgMTQGkk1xrDVLuAgwGHH1VKNSul/kkplZnVuYTIAMHKeLT2TtI/OZ+xq53C+ewH67h5UylvtQ9jUrBuTeoDxbaqIj5ybSXffPUsLT0TKX/+ofOjHGjt5zdu30SZkSuTbpEEimAzYctTKUOdE/ZapZQd+CHw21rrSePw3wObgJ1AH/CXQRul1OeUUoeUUoeGhpK7Z68QmcpfGDBwA6OXTw6gFNx5TfYFCpNJ8Zef2IEjL4d1a/KxWlK/3kYpxZ8+2MSaAiu/89SxlGdpf/uNDsoKbXz21o0pfW44kfwtdAO1Ad/XAL2hzlFKWQAHMBruWqVUDr4g8X2t9Y/8J2itB7TWHq21F/hf+Ia+rqC1/o7WerfWendZWVkEb0OIq0+Zv0cRsIHRq6cGuH5dCaX2zPg0Gq1KRx7/+l/38D9/eXva2lBSYOVrD23n7OD00jBYKiy6vbzTMcKHm9ZmRBFHv0gCxftAg1Jqo1LKim9yev+yc/YDnza+fgh4TfsKuOwHHjZWRW0EGoD3jPmL7wJtWuu/CryRUqoy4NuPAS3RvikhVouloSejR9E3MUdLz2TS9p5IlR21xdy0Kb0T8R/aUs6nblzPd986x9vtwyl55rGL48y5PCnbuS5SKwYKY87hUeBFfJPOT2utW5VSTyilHjBO+y5QqpRqB74APGZc2wo8DZwEDgCf11p7gFuATwF3BlkG+zWl1AmlVDNwB/A7iXqzQlxtHHk5mE1qaY7i1bZBAO5pzL5hp0z0+Ievoc5ZwBd/cDwlW6T+wpibybTVahH1bYwlqs8vO/blgK/ngY+HuPYrwFeWHXuL4PMXaK0/FUmbhBC+Mf3SAisjRtLdK20DrC/NZ1MGLKm8GuRbLfzVJ3fyK3//Nn+yv5W/+uTOpD7vnY4RmqodKS1bEgnJzBYiy/mT7mYW3LzdMcLdWysyPhs7m+ysLebRO+r50dEefprEJbOzi26OXhzj5k2ZNewEEiiEyHr+ek8/PzvMotublctiM92jd9azvcbBH/74BIOT80l5xvvnx3B5NDeneW4mGAkUQmQ5Z4Gv3tOrbQMU5Vq4Ic11ga5GOWYT3/jkTuZdHn7vmeaQmy3F4+32YaxmU0b+/UmgECLLOQttDE0v8NqpQT60pTxtexZc7TaV2Xn8/q387MwQ/3awK+H3/0XHMLvWFZNnzbxsevmJEiLLOe1WFt1eRmYWM2pv7KvRp25cz60NTv7sp210JrAU+fjsIq29kxk5PwESKITIev69sy0mxe2bJfk0mUwmxdcf2oHVYuJ3nj6esAqz73aOoDXcUp958xMggUKIrOffO3vPxjUZt6zyarTWkcufPtjE8YvjfPuNjoTc8xftI+RbzeyoLU7I/RJNAoUQWa6iyBcosj0bO5t8dEcV+3ZW8c1Xz3L84njc93u7Y5g9G9dk7PxSZrZKCBGxLRWF/M0ju/jVD6xLd1NWlSceaKK80Mb/92x8VYb6J+bpGJrhlgydnwAJFEJkPaUUD+yoyqq9J64Gjvwc/tvtm2junoirHPnbHb46Ujdn6PwESKAQQoiY7dtZhdVi4pnD3THf4+2OEUryc9i6tiiBLUssCRRCCBGj4nwr921by38e7Ylp3wqtNW+3D3PTplJMpswtuyKBQggh4vCJ3TVMzLl4+eRA1NdeGJmld2KemzJ4fgIkUAghRFxu2eSkujiPpw9dXPnkZX5hzE/ckoH1nQJJoBBCiDiYTIqHrq/hrfZhusdmo7r27fYRKh25bHSmfm/waEigEEKIOH18dw1AVJPaXq/mnc4RbtpUmvFl4SVQCCFEnGpK8rllk5MfHOrG642ssuyp/ilGZxYzOn/CTwKFEEIkwCduqKVnfI53OkciOj8b8if8JFAIIUQC3NtYQVGuhafej2xS++2OEeqcBVQ68pLcsvhJoBBCiATIzTHz4K5qDrT2MzHrCnuuy+PlYOdIVvQmQAKFEEIkzCd217Lo9vLs8Z6w5zV3jzOz6MnY/SeWk0AhhBAJ0lTtoLGyaMWcirfbR1AKbqqTHoUQQqw6n7yhlpaeSVp7QxcK/EXHMI2VRZQUWFPYsthJoBBCiATat7MKq9nEDw4Fz6mYW/Rw5MI4N2d4NnYgCRRCCJFAxflW7t1WEbJQ4OELYyx6vNxcnx3zEyCBQgghEu6TN9SGLBT4i45hLCbFng1r0tCy2EigEEKIBAtXKPDtjhF21hZTYLOkoWWxkUAhhBAJFlgosGd8bun4xJyLE93jWTXsBBIohBAiKR663igUGDCpfbBzBK8mqyayQQKFEEIkRe0ao1Dg4YtLhQLf7hghN8fErnXFaW5ddCRQCCFEknx8dw3dY5cKBb7dMcwNG9Zgs5jT3LLoSKAQQogkuW/bWopyLTx96CKDU/OcGZjOmrIdgbJn2l0IIbKMv1Dgk+9fZPf6EgBuyZJCgIGkRyGEEEnkLxT4Fy+doSjXwrYqR7qbFDUJFEIIkUT+QoETcy5urCvFbMrsbU+DkUAhhBBJ9gljT+1bsix/wk/mKIQQIske2l1L1+gcD+yoSndTYiKBQgghksxus/Dljzamuxkxi2joSSm1Vyl1WinVrpR6LMjrNqXUU8brB5VSGwJee9w4flopdZ9xrFYp9bpSqk0p1aqU+q2A89copV5WSp01/lsS/9sUQggRqxUDhVLKDHwLuB9oBB5RSi0PjZ8BxrTW9cA3gK8a1zYCDwPbgL3At437uYHf1VpvBW4EPh9wz8eAV7XWDcCrxvdCCCHSJJIexR6gXWvdqbVeBJ4E9i07Zx/wPePrZ4C7lFLKOP6k1npBa30OaAf2aK37tNZHALTWU0AbUB3kXt8DHoztrQkhhEiESAJFNRBYK7ebS7/UrzhHa+0GJoDSSK41hql2AQeNQxVa6z7jXn1AeQRtFEIIkSSRBIpgi351hOeEvVYpZQd+CPy21noygrZceqBSn1NKHVJKHRoaGormUiGEEFGIJFB0A7UB39cAvaHOUUpZAAcwGu5apVQOviDxfa31jwLOGVBKVRrnVAKDwRqltf6O1nq31np3WVlZBG9DCCFELCIJFO8DDUqpjUopK77J6f3LztkPfNr4+iHgNa21No4/bKyK2gg0AO8Z8xffBdq01n8V5l6fBp6N9k0JIYRInBXzKLTWbqXUo8CLgBn4J611q1LqCeCQ1no/vl/6/1sp1Y6vJ/GwcW2rUupp4CS+lU6f11p7lFIfBD4FnFBKHTMe9Qda6+eBPweeVkp9BugCPp7INyyEECI6yvfBP7sppYaACzFe7gSGE9icZMj0NmZ6+yDz25jp7QNpYyJkWvvWa61XHLu/KgJFPJRSh7TWu9PdjnAyvY2Z3j7I/DZmevtA2pgImd6+UKQooBBCiLAkUAghhAhLAgV8J90NiECmtzHT2weZ38ZMbx9IGxMh09sX1KqfoxBCCBGe9CiEEEKEtaoDxUrl09MpXCn2TKOUMiuljiqlfpLutiynlCpWSj2jlDpl/L+8Kd1tWk4p9TvG33GLUuo/lFK5GdCmf1JKDSqlWgKOZcwWACHa93Xj77lZKfWfSqnidLUvVBsDXvuiUkorpbJiy7tVGygiLJ+eTuFKsWea38JXATgTfRM4oLW+BthBhrVTKVUN/CawW2vdhC+p9eH0tgqAf8G3NUCgTNoC4F+4sn0vA01a6+3AGeDxVDdqmX/hyjailKoF7sGXUJwVVm2gILLy6WmzQin2jKGUqgE+AvxjutuynFKqCLgNX+UAtNaLWuvx9LYqKAuQZ9RJy+fKWmopp7V+E1+VhUAZswVAsPZprV8yqlcDvIuvtlzahPh/CL49e36fK4urZqzVHCgiKZ+eEYKUYs8kf43vh96b7oYEUQcMAf9sDI39o1KqIN2NCqS17gH+At+nyz5gQmv9UnpbFVI2bQHwX4EX0t2I5ZRSDwA9Wuvj6W5LNFZzoIikfHraxVOKPdmUUr8EDGqtD6e7LSFYgOuAv9da7wJmyLAdE41x/n3ARqAKKFBK/Vp6W5XdlFJ/iG/o9vvpbksgpVQ+8IfAl9Pdlmit5kARSfn0tApTij1T3AI8oJQ6j2/o7k6l1L+lt0mX6Qa6tdb+ntgz+AJHJrkbOKe1HtJau4AfATenuU2hRLQFQDoppT4N/BLwqzrz1v5vwveB4Ljxb6YGOKKUWpvWVkVgNQeKSMqnp80Kpdgzgtb6ca11jdZ6A77/f69prTPm07DWuh+4qJTaYhy6C18l40zSBdyolMo3/s7vIsMm3ANk9BYASqm9wJeAB7TWs+luz3Ja6xNa63Kt9Qbj30w3cJ3xc5rRVm2gMCa9/OXT24Cntdat6W3VZW7BV4r9TqXUMePPh9PdqCz0/wDfV0o1AzuBP0tzey5j9HaeAY4AJ/D9m0x79q5S6j+Ad4AtSqluo+z/nwP3KKXO4lu18+cZ1r6/AwqBl41/L/+QrvaFaWNWksxsIYQQYa3aHoUQQojISKAQQggRlgQKIYQQYUmgEEIIEZYECiGEEGFJoBBCCBGWBAohhBBhSaAQQggR1v8P8AKpeCPDavsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "G\n",
      "37.191417999999885\n"
     ]
    }
   ],
   "source": [
    "#event_count = 100000\n",
    "#event_ratio = 0.5\n",
    "#noise_sigma = 5\n",
    "#signal_mu = 100\n",
    "#signal_sigma = 20\n",
    "#threshold = 20\n",
    "#error_of_expected_signal(event_count, event_ratio, signal_mu, signal_sigma, noise_sigma, threshold)\n",
    "#prng = ROOT.TRandom3(123)\n",
    "\n",
    "\n",
    "def get_number_above_threshold(threshold, count_, list_):\n",
    "    arr = np.fromiter(list_, dtype=np.float, count=count_)\n",
    "    above_count = (arr > threshold).sum()\n",
    "    return above_count\n",
    "\n",
    "def error_of_expected_signal(prng, event_count, event_ratio, signal_mu, signal_sigma, noise_sigma, threshold):\n",
    "    low = -100\n",
    "    high = 500\n",
    "    \n",
    "    # expected\n",
    "    fNoise = ROOT.TF1(\"ElecNoiseBg\", \"TMath::Gaus(x, 0, %f)\" % noise_sigma, low, high)\n",
    "    fNoise.SetNpx(1000)\n",
    "    fSignal = ROOT.TF1(\"Signal\", \"TMath::Landau(x, %f, %f)\" % (signal_mu, signal_sigma), low, high)\n",
    "    fSignal.SetNpx(1000)\n",
    "\n",
    "    model_expected_signal_above = 1 - fSignal.Integral(low, threshold)/fSignal.Integral(low, high)\n",
    "    model_expected_noise_above = 1 - fNoise.Integral(low, threshold)/fNoise.Integral(low, high)\n",
    "    # threshold_noise / all_noise\n",
    "    \n",
    "    \n",
    "    # simulated\n",
    "    signal_count = int(round(event_count * event_ratio))\n",
    "    noise_count = event_count - signal_count\n",
    "\n",
    "    signal_above_count = 0\n",
    "    noise_above_count = 0\n",
    "    \n",
    "\n",
    "    # like 17% faster, almost not worth it\n",
    "    signal_iter = [ prng.Landau(signal_mu, signal_sigma) for _ in range(signal_count) ]\n",
    "    signal_above_count = get_number_above_threshold(threshold, signal_count, signal_iter)\n",
    "\n",
    "    noise_iter = [ prng.Gaus(0, noise_sigma) for _ in range(noise_count) ]\n",
    "    noise_above_count = get_number_above_threshold(threshold, noise_count, noise_iter)\n",
    "    \n",
    "    \n",
    "    #det här väret är det vi får reda på från detektorn\n",
    "    total_above_count = signal_above_count+noise_above_count\n",
    "    total_count = total_above_count/(event_ratio*model_expected_signal_above + (1-event_ratio)*model_expected_noise_above)\n",
    "    computed_signal_count = total_count*event_ratio\n",
    "    deviation = abs(1 - computed_signal_count/signal_count)\n",
    "    \n",
    "    #print((signal_above_count/expected_signal_above)/signal_count, (noise_above_count/expected_noise_above)/noise_count)\n",
    "    \n",
    "    return deviation\n",
    "\n",
    "\n",
    "def find_best_threshold_and_error(prng, event_count, event_ratio, signal_mu, signal_sigma, noise_sigma):\n",
    "    low = 0 # cannot be below 0\n",
    "    high = 20\n",
    "    tests = high-low\n",
    "    \n",
    "    deviation_iter = [error_of_expected_signal(prng, event_count, event_ratio, signal_mu, signal_sigma, noise_sigma, threshold) for threshold in range(low, high)]\n",
    "    deviations = np.fromiter(deviation_iter, dtype=np.float, count=tests)\n",
    "    \n",
    "    # blur the errors, assumes errors follow a smooth curve, minimum is still minimum.\n",
    "    size = 5  # default is 3\n",
    "    deviations_blur = uniform_filter(deviations, size=size)#[1:-1] # nah no need to slice, the boundary error doesn't bother us\n",
    "    \n",
    "    # get the minimum\n",
    "    index = np.argmin(deviations_blur)\n",
    "    threshold = low + index\n",
    "    \n",
    "    # return the best found threshold and the error for it\n",
    "    #return (threshold, deviations[index])\n",
    "    return (threshold, deviations_blur[index])\n",
    "    \n",
    "    \n",
    "# Average eller expected ???\n",
    "def caluclate_SB(expected_signal_energy, noise_sigma):\n",
    "    return expected_signal_energy/noise_sigma\n",
    "\n",
    "def calculate_expected_signal_voltage(signal_mu, signal_sigma, reference_signal_mu_for_sigma):\n",
    "    #low = 0  # There are no particles hitting the detector with negative energies.\n",
    "    #high = 500\n",
    "    #fSignalEnergy = ROOT.TF1(\"Signal\", \"TMath::Landau(x, %f, %f)*x\" % (signal_mu, signal_sigma * (signal_mu/reference_signal_mu_for_sigma)), low, high)\n",
    "    # documents say minimum change if 4 fC - i dunno, just gonna say it is a scale thingy\n",
    "    #charge = 1 \n",
    "    #scale = 1 # 0.001 # mili, no sigma is also in mV so it doesn't matter\n",
    "    #return fSignalEnergy.Integral(low, high) * charge * scale\n",
    "    return signal_mu\n",
    "    \n",
    "    \n",
    "    \n",
    "def calculate_average_signal_voltage():\n",
    "    pass\n",
    "\n",
    "def linspace(start, end, steps):\n",
    "    scale = (end-start)/(steps - 1)\n",
    "    return map(lambda x : start + x*scale, range(steps))\n",
    "\n",
    "def calculate_SB_array(signal_sigma, noise_sigma, iterator, steps):\n",
    "    return np.fromiter(\n",
    "                map(lambda signal_mu :\n",
    "                        caluclate_SB(\n",
    "                                #calculate_expected_signal_voltage(signal_mu, signal_sigma), \\\n",
    "                                signal_mu,\n",
    "                                noise_sigma) ,\n",
    "                        iterator)\n",
    "             , dtype=np.float, count=steps)\n",
    "\n",
    "def calculate_error_array(prng, event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma, iterator_data, steps):\n",
    "    # signal_sigma = 0.5 + signal_sigma * (signal_mu/reference_signal_mu_for_sigma)\n",
    "    return np.fromiter(\n",
    "                 map(lambda signal_mu :\n",
    "                         find_best_threshold_and_error(prng, event_count, event_ratio, signal_mu, 0.5 + signal_sigma * (signal_mu/reference_signal_mu_for_sigma), noise_sigma)[1] ,\n",
    "                         linspace(*iterator_data))\n",
    "                 , dtype=np.float, count=steps)\n",
    "\n",
    "# Single core version\n",
    "def plot_error_vs_BS(event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma):\n",
    "    min_signal_mu = 0 #mv # cannot be below 0\n",
    "    max_signal_mu = 20  #mv\n",
    "    steps = (max_signal_mu - min_signal_mu) * 1\n",
    "    \n",
    "    signal_mu_list = list(linspace(min_signal_mu, max_signal_mu, steps))\n",
    "    SB    = calculate_SB_array(signal_sigma, noise_sigma, linspace(min_signal_mu, max_signal_mu, steps), steps)\n",
    "    \n",
    "    seed = 123\n",
    "    error = calculate_error_array(ROOT.TRandom3(seed), event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma, (min_signal_mu, max_signal_mu, steps), steps)\n",
    "    \n",
    "    fig = plt.figure()\n",
    "    plt.title(\"error vs SB\")\n",
    "    plt.xlabel('SB')\n",
    "    plt.ylabel('error')\n",
    "    plt.plot(SB, error)\n",
    "    plt.show()\n",
    "    ax = fig.add_subplot(111)\n",
    "    fig.savefig(\"SB_error_c\"+str(event_count)+\".pdf\", dpi=200)\n",
    "    plt.clf()\n",
    "    plt.plot(signal_mu_list, error)\n",
    "    plt.show()\n",
    "\n",
    "# Multi core version\n",
    "def plot_error_vs_BS_multicore(event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma):\n",
    "    # multi core, set up pool of threads/processes\n",
    "    max_nprocs = 99  # lower this to like 4 if you like\n",
    "    nprocs = min(mp.cpu_count(), max_nprocs)\n",
    "    print(f\"Number of CPU cores to use: {nprocs}\")\n",
    "    pool = mp.Pool(processes=nprocs)\n",
    "    \n",
    "    # how many steps for each mu to test\n",
    "    min_signal_mu = 0 #mv # cannot be below 0\n",
    "    max_signal_mu = 15  #mv\n",
    "    steps = (max_signal_mu - min_signal_mu) * 2\n",
    "    \n",
    "    # setup the x-axis lists\n",
    "    signal_mu_list = list(linspace(min_signal_mu, max_signal_mu, steps))\n",
    "    SB = calculate_SB_array(signal_sigma, noise_sigma, linspace(min_signal_mu, max_signal_mu, steps), steps)\n",
    "    \n",
    "    seed = 123\n",
    "    #seed = int(time.time()*1000)\n",
    "    #prng = ROOT.TRandom3(seed)\n",
    "    print(\"C\")\n",
    "    # create argument list\n",
    "    param_list = [(ROOT.TRandom3(seed + i), event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma, (min_signal_mu, max_signal_mu, steps), steps) \\\n",
    "                     for i in range(nprocs)]\n",
    "    # calculate\n",
    "    results = pool.starmap(calculate_error_array, param_list)\n",
    "    error = sum(results)*(1/nprocs)\n",
    "    pool.close()\n",
    "    \n",
    "    # plot\n",
    "    fig = plt.figure()\n",
    "    plt.title(\"error vs SB\")\n",
    "    plt.xlabel('SB')\n",
    "    plt.ylabel('error')\n",
    "    \n",
    "    plt.plot(SB, error)\n",
    "    print(\"F\")\n",
    "    plt.show()\n",
    "    \n",
    "    ax = fig.add_subplot(111)\n",
    "    fig.savefig(\"SB_error_c\"+str(event_count)+\"M\"+str(nprocs)+\"_with_dynamic_sigma.pdf\", dpi=200)\n",
    "    plt.clf()\n",
    "    plt.plot(signal_mu_list, error)\n",
    "    plt.show()\n",
    "    print(\"G\")\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "\n",
    "def main():\n",
    "    # estimate: 10^6: 5 hours\n",
    "    #           10^5: 0.5 hours\n",
    "    # 500000: 3120 s\n",
    "    event_count = 400000\n",
    "    event_ratio = 0.4/4  # = 0.1\n",
    "    noise_sigma = 5\n",
    "    #signal_mu = 100\n",
    "    signal_sigma = 25\n",
    "    #threshold = 0\n",
    "    \n",
    "    # 25 sigma at S1M-2 # s.24\n",
    "    # mu at 100 # s.15\n",
    "    reference_signal_mu_for_sigma = 100\n",
    "\n",
    "    #error_of_expected_signal(event_count, event_ratio, signal_mu, signal_sigma, noise_sigma, threshold)\n",
    "    #find_best_threshold_and_error(event_count, event_ratio, signal_mu, signal_sigma, noise_sigma)\n",
    "    t0 = time.clock()\n",
    "    plot_error_vs_BS_multicore(event_count, event_ratio, signal_sigma, reference_signal_mu_for_sigma, noise_sigma)\n",
    "    t1 = time.clock()\n",
    "    print(t1-t0)\n",
    "\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    mp.freeze_support()\n",
    "    main()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}