{ "cells": [ { "cell_type": "markdown", "id": "8d1db2ab", "metadata": {}, "source": [ "### Calculate 95% limit with Bayesian method with a flat prior\n", "Marginalize the likelihood via MC integration.\n", "Then interpret the marginalized likelihood as a posterior pdf" ] }, { "cell_type": "code", "execution_count": 1, "id": "a49afec5", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "from iminuit import Minuit\n", "from scipy.stats import norm\n", "import math\n", "\n", "# Using the lognormal from numba_stats\n", "# unstead of scipy_stats speeds things\n", "# up by a big factor.\n", "# Be careful however because the call sequences\n", "# are not exactly the sameand not all\n", "# functions are implemented.\n", "# from scipy.stats import lognorm\n", "from numba_stats import lognorm" ] }, { "cell_type": "markdown", "id": "aa49c280", "metadata": {}, "source": [ "### This is fake data with 3 point" ] }, { "cell_type": "code", "execution_count": 2, "id": "e1da7c51", "metadata": {}, "outputs": [], "source": [ "# This is some fake data in 3 bins\n", "data = np.array([ 7, 6, 1]) # observed\n", "sig = np.array([ 2, 2, 2]) # signal predicted with mu=1\n", "bg = np.array([10, 5, 2]) # bg predicted\n", "err = 0.2*bg # bg uncertainty" ] }, { "cell_type": "code", "execution_count": 3, "id": "e2f96841", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD8CAYAAAB6paOMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkW0lEQVR4nO3dfXRU9b3v8feXgAKSoiIqARXqMTwJARIRgdKggmhRr0hbEY/Y1juW6rV1cYq2eH2iWcej2NOeo62Nzz0FS02RigWtVwFRfICwEgR5EisWgoCoCCIo+L1/zEBDmJDJzCR7z8zntVZWZvb3t/d8fw7+vtlPv23ujoiI5K4WQScgIiLBUiEQEclxKgQiIjlOhUBEJMepEIiI5DgVAhGRHNdgITCzU8xsvpm9bWYrzezHseXHm9kLZrYu9vu4etafEGuzzswmpLsDIiKSGmvoPgIz6wR0cvdlZpYPVAL/C7gG+Mjd7zazW4Dj3P3mOuseDywFSgCPrVvs7h+nuyMiIpKcBvcI3H2zuy+Lvd4JrAI6A5cCT8SaPUG0ONR1AfCCu38UG/xfAEalIW8REUmTlo1pbGZdgf7AG8BJ7r45FvoAOCnOKp2Bf9R6vzG2LN62I0AE4Jhjjinu0aNHY1ITEclplZWVH7p7x2TWTbgQmFk74M/AT9z9UzM7GHN3N7OU5qpw93KgHKCkpMSXLl2ayuZERHKKmW1Idt2Erhoys1ZEi8B0d58VW7wldv7gwHmErXFW3QScUut9l9gyEREJiUSuGjLgEWCVu/+yVugZ4MBVQBOAv8RZ/XlgpJkdF7uqaGRsmYiIhEQiewRDgH8FzjWzqtjPRcDdwAgzWwecH3uPmZWY2cMA7v4RMBVYEvu5K7ZMRERCosHLR4OgcwQizefLL79k48aN7NmzJ+hUJAGtW7emS5cutGrV6pDlZlbp7iXJbLNRVw2JSPbZuHEj+fn5dO3aldoXgUj4uDvbt29n48aNdOvWLW3b1RQTIjluz549dOjQQUUgA5gZHTp0SPvemwqBiKgIZJCm+K5UCEREcpzOEYjIIWpq0ru9goKG2+Tl5dGnTx/cnby8PO6//34GDx4MwJtvvsnkyZPZtGkT+fn5dOrUibvvvps+ffqkN9EcpkIgIoFr06YNVVVVADz//PP87Gc/Y+HChWzZsoXvfOc7zJgx42BheOWVV1i/fr0KQRqpEIhIqHz66accd1x0Vvv777+fCRMmHCwCAEOHDg0qtaylQiAigfv888/p168fe/bsYfPmzbz00ksArFy5kgkT9BiTpqaTxSISuAOHhlavXs1zzz3H1VdfTbybXc8++2x69uzJj3/84wCyzF4qBCISKueccw4ffvgh27Zto3fv3ixbtuxg7I033mDq1Kns2LEjwAyzjwqBiITK6tWr2b9/Px06dOD666/n8ccfZ/HixQfju3fvDjC77KRzBCJyiEQu90y3A+cIIDqNwhNPPEFeXh4nn3wyM2fO5Oabb2bTpk2ceOKJnHDCCdx2223Nn2QWUyEQkcDt37+/3tigQYNYuHBhM2aTe3RoSEQkx6kQiIjkOBUCEZEcp0IgIpLjVAhERHJcg1cNmdmjwGhgq7ufGVs2E+gea3Is8Im794uz7nvATmA/sC/Zx6iJiEjTSeTy0ceB+4HfH1jg7t898NrM7gOOdJvfcHf/MNkERUSkaTVYCNz9ZTPrGi9m0UflfAc4N815iYhIM0n1HME3gC3uvq6euAN/M7NKM4uk+Fki0lwmTYL77ou+LiiIPq1mwQIoLY0ui0SgvDz6Oj8fdu6EOXPg4oujy668EmbMiL5O4dGKtaefbmrf//73OfHEEznzzDOb7TPDwuLN8HdYo+gewbMHzhHUWv5b4B13v6+e9Tq7+yYzOxF4Afg/7v5yPW0jQATg1FNPLd6wYUOjOiIiyVm1ahU9e/b854LiYqisDC6hgLz88su0a9eOq6++mhUrVgSdzhEd9p0BZlaZ7HnYpPcIzKwlMAaYWV8bd98U+70VeBoYeIS25e5e4u4lHTt2TDYtEUnVgb/00+XAnsERfPbZZ3zrW9+iqKiIM888k5kzo8NKu3btAJg6dSrdu3dn6NChjBs3jmnTpsXdTmlpKatXrwZg+/btjfrrftiwYRx//PEJt88mqcw1dD6w2t03xgua2TFAC3ffGXs9Ergrhc8TkeaQn5/e7T37bPRQ0RE899xzFBQU8Ne//hXgkGmmlyxZwp///Geqq6v58ssvGTBgAMXFxXG3884771BYWAjA8uXLDz7O8hvf+AY7d+48rP20adM4//zzk+pWNknk8tEngVLgBDPbCNzu7o8AVwBP1mlbADzs7hcBJwFPR88n0xKY4e7PpTd9EUm70aNh7dr0bS+BPYI+ffowadIkbr75ZkaPHs03vvGNg7FXX32VSy+9lNatW9O6dWsuPnAeoo4NGzbQuXNnWrSIHuhYvnw5ffv2BWDRokVp6Ej2SuSqoXH1LL8mzrIa4KLY63eBohTzE5Hmls4iANETyHPmHLFJYWEhy5YtY+7cudx6662cd955jZ5qurq6+uDAD1BZWcl3vxu90l17BEemO4tF5FB33JHe7UUavmCwpqaGtm3bctVVV/HTn/70kKeSDRkyhDlz5rBnzx527drFs88+G3cbVVVV7NmzB4B169bxl7/85eChoUWLFlFVVXXYj4pAlAqBiBzqzjv/WQwKC6N7CJWV0auJoPGXlh5YfgRvvfUWAwcOpF+/ftx5553ceuutB2NnnXUWl1xyCX379uXCCy+kT58+tG/f/rBtVFdX89VXX1FUVMRdd91Fr169eOKJJxLu9rhx4zjnnHNYs2YNXbp04ZFHHkl43UyX0OWjza2kpMSXLl0adBoiOSHepYhhs2vXLtq1a8fu3bsZNmwY5eXlDBgw4JA2Z5xxBsuWLSM/3Se7Qyjdl4/qCWUiEnqRSIS3336bPXv2MGHChMOKwM6dOzGznCgCTUGFQERCb0YDVx7l5+ezNt0nuXOIzhGIiOQ4FQIRkRynQiAikuNUCEREcpwKgYhIjlMhEBHJcSoEInKIBQugY8fo76aI1/Xee+9l3MNgfvWrX7F79+4G21177bW8/fbbR2wze/bsBts0NRUCETlEaSk89RR8+9vxB/NU45lu//79CReChx9+mF69eh2xjQqBiIRSOotBIvbt28f48ePp2bMnY8eOPTjIzp07lx49elBcXMyNN97I6NGjD1v38ccf54Ybbjj4fvTo0SyIJdWuXTumTJlCUVERgwYNYsuWLQBs2bKFyy67jKKiIoqKili8eDEAf/jDHw7OeXTdddexf//+g9uZNGkSRUVFlJWVUVNTw/Dhwxk+fDgAEydOpKSkhN69e3P77bfX+u9QyoHpcuLlsnjxYp555hl++tOf0q9fP9avX3/IXdPr1q077C7qpqBCICJxpasYJGLNmjX86Ec/YtWqVXzta1/jN7/5DXv27OG6665j3rx5VFZWsm3btkb34bPPPmPQoEFUV1czbNgwHnroIQBuvPFGvvnNb1JdXc2yZcvo3bs3q1atYubMmbz66qtUVVWRl5fH9OnTD27n7LPPprq6mttuu42CggLmz5/P/PnzASgrK2Pp0qUsX76chQsXsnz58oRyGTx4MJdccgn33nsvVVVVnH766bRv356qqioAHnvsMb73ve81ut+NpUIgIvVKRzFIxCmnnMKQIUMAuOqqq3jllVdYvXo1X//61+nWrRsQnR20sY466qiDexHFxcW89957ALz00ktMnDgRgLy8PNq3b8+LL75IZWUlZ511Fv369ePFF1/k3XffPdjm8ssvr/dz/vSnPzFgwAD69+/PypUr4x7qqS+Xuq699loee+wx9u/fz8yZM7mygae7pYPmGhKRI6o92D/11OGDe0PxRMSeZFjv+yNp2bIlX3311cH3B55JANCqVauD28rLy2Pfvn31bsfdmTBhAv/+7/9+WKx169bk5eXFXe/vf/8706ZNY8mSJRx33HFcc801h+TQ2Fwuv/xy7rzzTs4991yKi4vp0KFDvTmni/YIROQQQZwgfv/993nttdeA6ARzQ4cOpXv37rz77rsH/3I+8ED7urp27UpVVRVfffUV//jHP3jzzTcb/LzzzjuP3/72t0D05O+OHTs477zzqKioYOvWrQB89NFHbNiwIe76+fn5B5949umnn3LMMcfQvn17tmzZwrx58xrV99rbgmjRueCCC5g4cWKzHBYCFQIRqaO+v+hLS2HbtuTjR9K9e3ceeOABevbsyccff8zEiRNp06YNv/nNbxg1ahTFxcXk5+fHfSDNkCFD6NatG7169eLGG29M6OTqr3/9a+bPn0+fPn0oLi7m7bffplevXvziF79g5MiR9O3blxEjRrB58+a460ciEUaNGsXw4cMpKiqif//+9OjRgyuvvPLgIa5EXXHFFdx7773079+f9evXAzB+/HhatGjByJEjG7WtZDX4YBozexQYDWx19zNjy+4A/jdw4OzNz919bpx1RwG/BvKIPtT+7kSS0oNpRJpPmB9Mc+CBNO7O9ddfzxlnnMFNN90UdFpNbtq0aezYsYOpU6fGjQfxYJrHgfuB39dZ/p/uPq2+lcwsD3gAGAFsBJaY2TPuHuwFsyKSMR566CGeeOIJvvjiC/r37891110XdEpN7rLLLmP9+vW89NJLzfaZDRYCd3/ZzLomse2BwDvu/i6Amf0RuBRQIRCRhNx00005sQdQ29NPP93sn5nKOYIbzGy5mT1qZsfFiXcG/lHr/cbYsrjMLGJmS81saTLXC4uISHKSLQS/BU4H+gGbgftSTcTdy929xN1LOnbsmOrmREQkQUkVAnff4u773f0r4CGih4Hq2gScUut9l9gyEREJkaQKgZl1qvX2MmBFnGZLgDPMrJuZHQVcATyTzOeJSLiUlpZSmsx1ohJKDRYCM3sSeA3obmYbzewHwD1m9paZLQeGAzfF2haY2VwAd98H3AA8D6wC/uTuK5uoHyKSRe644w6mTav3osRQzNiZTRK5aijeBB+P1NO2Brio1vu5wGH3F4hI5po+fTqvv/46e/fupWvXrpSVlTF+/PhmzWH27NmMHj26wSmeJTG6s1hEEjZ9+nQikQh79+4FYMOGDUQikYOzdKairKyMwsJChg4dypo1a4DofQRnnXUWRUVFXH755ezevTvu1M3x2kniVAhEJGFTpkw5bJDdvXs3U6ZMSWm7lZWV/PGPf6Sqqoq5c+eyZMkSAMaMGcOSJUuorq6mZ8+ePPLII3Gnbo7XThKn2UdFJGHvv/9+o5YnatGiRVx22WW0bdsWgEsuuQSAFStWcOutt/LJJ5+wa9cuLrjggrjrJ9pO4tMegYgk7NRTT23U8lRdc8013H///bz11lvcfvvtcad3bkw7iU+FQEQSVlZWdvCv9gPatm1LWVlZStsdNmwYs2fP5vPPP2fnzp3MmTMHgJ07d9KpUye+/PLLQ85D1J26ub52khgVAhFJ2Pjx4ykvL+foo48G4LTTTqO8vDzlq4YGDBjAd7/7XYqKirjwwgs566yzAJg6dSpnn302Q4YMoUePHgfb1526ub52kpgGp6EOgqahFmk+yUxDfeBmsgXJPIVGUhbENNQiIodQAcguOjQkIpLjVAhEhDAeIpb4muK7UiEQyXGtW7dm+/btKgYZwN3Zvn07rVu3Tut2dY5AJMd16dKFjRs3ogdCZYbWrVvTpUuXtG5ThUAkx7Vq1Ypu3boFnYYESIeGRERynAqBiEiOUyEQEclxKgQiIjlOhUBEJMcl8sziR81sq5mtqLXsXjNbbWbLzexpMzu2nnXfiz3buMrMNHmQiEgIJbJH8Dgwqs6yF4Az3b0vsBb42RHWH+7u/ZKdDElERJpWg4XA3V8GPqqz7G/uvi/29nUgvXc3iIhIs0nHOYLvA/PqiTnwNzOrNLPIkTZiZhEzW2pmS3WHo4hI80mpEJjZFGAfUN8jgYa6+wDgQuB6MxtW37bcvdzdS9y9pGPHjqmkJSIijZB0ITCza4DRwHivZ7Yqd98U+70VeBoYmOzniYhI00iqEJjZKGAycIm7766nzTFmln/gNTASWBGvrYiIBCeRy0efBF4DupvZRjP7AXA/kA+8ELs09MFY2wIzmxtb9STgFTOrBt4E/uruzzVJL0REJGmJXDU0zt07uXsrd+/i7o+4+7+4+ymxy0L7ufsPY21r3P2i2Ot33b0o9tPb3cuSSXDBAujYMfpb8cyLi0gGcPfQ/RQXF3tt8+e7n3BC9Hc8ioc7LiJND1jqSY65gQ/68X7qFgL34AczxROLi0gwsq4QFBYeXgjcwzPYKV5/XESCkXWFoGXL4lAPdopr0BcJm6wrBIWFxaEe7BRXMRAJm6wrBMXFxYEPZoqrGIhkkqwsBO7BD2aKqxiIZIqsLQTuwQ9miqsYiGSCrC4E7sEPZoqrGIiEXdYXAvfgBzPFdR+BSJhlXSHQfQSZGxeRYGRdIdB9BJkbF5FgZF0h0H0EmR0XkeaXdYVA9xFkflxEmldWFgL34AczxVUMRDJF1hYC9+AHM8VVDEQyQVYXAvfgBzPFVQxEwi7rC4F78IOZ4rqPQCTMmrwQAI8CW4EVtZYdD7wArIv9Pq6edSfE2qwDJiTyebqPIHPjIhKM5igEw4ABdQrBPcAtsde3AP8RZ73jgXdjv4+LvY5bMGr/6D6CzI2LSDCa5dAQ0LVOIVgDdIq97gSsibPOOOB3td7/DhjX0GfpPoLMjotI80ulEFh0/YaZWVfgWXc/M/b+E3c/NvbagI8PvK+1zr8Brd39F7H3/xf43N2nxdl+BIgAdO58avGvfrWBSATKy2Hw4MPzWbwYxUMYLyg4vK2IND0zq3T3kmTWbZGOBGLVKLGKUv82yt29xN1LOnToyODB0UEmEokOOnUpHu64iGSOVArBFjPrBBD7vTVOm03AKbXed4ktS0jQg5niKgYiuSCVQvAM0SuCiP3+S5w2zwMjzew4MzsOGBlblrCgBzPFVQxEsl1ChcDMngReA7qb2UYz+wFwNzDCzNYB58feY2YlZvYwgLt/BEwFlsR+7oota5SgB7N0xceNK2XEiNLQ5peOuIhknoRPFjen008v8UWLlh62PKwnSBONjxhRytq18OSTC0KZXzriY8cevlxEml7gJ4vTbcOG4P+ybYp4+/ZQWBje/NIRF5HME8pCcNpp4R7sUi0GYc5P5wREck8oC0G7duEe7BRXMRDJJqEsBBD8YKa4ioFIrghtIYDgBzPFVQxEckGoCwEEP5ilM75jR7jzUzEQyU2hLwQQ/GCWrvjatdlfDEQk84SyEOzadfiyMA12ycYLC6PFIKz5pSMuIpknlIVA9xFkblxEMk8oC4HuI8jcuIhknlAWAt1HkNlxEcksoSwEEPxgpriKQV2lpaWUlpYGnYZI2oW2EEDwg5niKgYiuSDUhQCCH8x0H4GKgUi2C30hgOAHM91HoPsIRLJZKAuB7iPI3LiIZJ5QPpimqKjE5807/ME0mW7s2FIAKioWBJpHUyooCDqDpnPgRPGCBQsCzUMknqx7MI2IiDSfpAuBmXU3s6paP5+a2U/qtCk1sx212tyWcsYiIpJWLZNd0d3XAP0AzCwP2AQ8HafpIncfnezniIhI00rXoaHzgPXuviFN2xMRkWaSrkJwBfBkPbFzzKzazOaZWe/6NmBmETNbamZLt2/flqa0RESkISkXAjM7CrgEeCpOeBlwmrsXAf8NzK5vO+5e7u4l7l7SoUPHVNMSEZEEpWOP4EJgmbtvqRtw90/dfVfs9VyglZmdkIbPFBGRNElHIRhHPYeFzOxkM7PY64Gxz9uehs/MOLNmTaey8nVee20hAwd2Zdas6UGnJCICpHDVEICZHQOMAK6rteyHAO7+IDAWmGhm+4DPgSs8jHewNbFZs6YzeXKEL77YC8CmTRuYPDk6H8OYMeODTE1ERHcWN4eBA7uyadPhF1R17nwab775XvMn1IR0Z7FIMHRnccjV1LzfqOUiIs1JhaAZFBSc2qjlIiLNSYWgGdxySxlt2rQ9ZFmbNm255ZaygDISEfknFYJmMGbMeO65p5yjjjoaiJ4buOeecp0oFpFQSOmqIUncmDHjmTHjISC7p6EWkcyjPQIRkRynQiAikuNUCEREcpwKgUgCpk+fzuuvv87ChQvp2rUr06drihDJHioEIg2YPn06kUiEvXujU4Rs2LCBSCSiYiBZQ4VApAFTpkxh9+7dhyzbvXs3U6ZMCSgjkfRSIZC0qm8angULoGPHzIy//378qUBqLw9z/o2JS25SIZC0+va34w82paXw1FOZGT/11PhTgdReHub8GxOX3KRCIGkVhsEs3fGysjLatj10ipCjj25LWVlZQutnUlxykwqBpFUYBrN0x8ePH095eTlHHx2dIuSkk06jVatyOnc+fIqQMObfmLjkJhUCSbugB7OmKgaDBg3im9/8Jh988B5z5owPVX7pjEvuUSGQJhH0YKa4ioEkToVAmkzQg5niKgaSmJQLgZm9Z2ZvmVmVmR32fEmL+i8ze8fMlpvZgFQ/UzJH0INZOuOffBLu/FQMJFnp2iMY7u796nle5oXAGbGfCPDbNH2mhFAYB7N0xVeuzP5iILmpOQ4NXQr83qNeB441s07N8LkSgDAMZk0V7907WgzCmp/uI5BkmbuntgGzvwMfAw78zt3L68SfBe5291di718Ebnb3pXXaRYjuMdC586nFb765IaW8wmjs2FIgux9Ms3gxRCJQXg6DB2dXfOzYUnbsgM2bF4Qyv3TECwoOXyaZwcwq6zkq06B07BEMdfcBRA8BXW9mw5LZiLuXu3uJu5d06NAxDWlJEAYPjg4ykUh00Mm2ePv24c4v1bjkppQLgbtviv3eCjwNDKzTZBNwSq33XWLLJEsFPZgprmIgjZNSITCzY8ws/8BrYCSwok6zZ4CrY1cPDQJ2uPvmVD5Xwi/owUxxFQNJXKp7BCcBr5hZNfAm8Fd3f87MfmhmP4y1mQu8C7wDPAT8KMXPlAwR9GCmuIqBJCalQuDu77p7Ueynt7uXxZY/6O4Pxl67u1/v7qe7e5+6J4kluwU9mKUzvmNHuPNTMZBk6c5iSaswDmbpiq9dm/3FQHKTCkEzqqhYkNWXjkI4BrOmihcWRotBWPNLR1xykwqBpFUYBrOmvHS0sDC8+aUjLrlJhUDSKgyDme4j0DkBaRwVAkm7oAczxVUMpHFSnmKiKRQVlfi8ebq4SMIlF6YI0RQTmSvoKSZERCSDqRCIiOQ4FQIRkRynQiAikuNaBp2ASKbI5pPEktu0RyAikuNUCEREcpwKgYhIjlMhEBHJcSoEIiI5ToVARCTHqRCIiOS4pAuBmZ1iZvPN7G0zW2lmP47TptTMdphZVeznttTSFRGRdEvlhrJ9wCR3X2Zm+UClmb3g7m/XabfI3Uen8DkiItKEkt4jcPfN7r4s9nonsAronK7ERESkeaTlHIGZdQX6A2/ECZ9jZtVmNs/Meqfj80REJH1SnmvIzNoBfwZ+4u6f1gkvA05z911mdhEwGzijnu1EgAhA586nppqWiIgkKKU9AjNrRbQITHf3WXXj7v6pu++KvZ4LtDKzE+Jty93L3b3E3Us6dOiYSloikkYLFkDHjtHfioc3nopUrhoy4BFglbv/sp42J8faYWYDY5+3PdnPFJGmFW+wKS2Fp56Cb39b8TDHU5HKHsEQ4F+Bc2tdHnqRmf3QzH4YazMWWGFm1cB/AVd4GB+SLCJA8IOZ4snHU6GH14vIQWvXRgebp56KP7gsWKB4WONZ9/D6ln9/B4Bjr7+SNk/PAKCgswHQ5ukZHHv9lQAcP+Fijv7bHGzXTk4uzAeg7R/KaT85AkCHsaUctXgBLT6o4aQBBQAc8+B9fO3OSQCcMKqYVssryVu/lhOHFgKQf98d5N93BwAnDi0kb/1aWi2v5IRRxQB87c5JHPPgfQCcNKCAFh/UcNTiBXQYWwpA+8kR2v6hHICTC/OxXTs5+m9zOH7CxeqT+hT6PvV5bw4P/XInxcPzqaiAT+4p57OrItTUwN7BpZz54QIeK6uh+7kFVFTAjtvuY9d1k6ipgS/6FtN7TyX/83/X0uW8QioqYOekO9g56Q5qamDf1wvpmbeWJ/+tkmPPL6aiAnZdN4kdt91HTQ3sP7mAHl+roeKGBbQcUUpFBXx2VYRP7imnpga+apdP94KdzP7BHPaMvJiKCth92ZV8/MAMamoAMwoL4dkrZ/DhyCupqIA9Iy5m++Nz2Lx2J5/lteTnPy/luTHlbLggQkVFtE8fVizgg2U17D+5gMJC+H8X3cfKUZOoqIj2adtzlWxZtJZ9Xy+ksBAWDr+DNy68g4qKaJ+2LFrLtucq+aJvMYWF8OqgSbx40X1UVET79MGyGj6sWMDewaUUFsIbRRH+8q1yKiqifdq8difbH5/DnhEXU1gIld2vZProGVRURPtUUwMfPzCD3ZddSWEhVJ92Mb8bPYfZ/7OTr9rlU1MT/Z5SEco9gn7dzvC5r64LOg2RnLV4MUQiUF4OgwdnR3xsrGBWVCwIZX6pxjt3zrI9Aj8mP+gURHLa4MHRQSYSiQ46imdWvLFCWQharaoOOgWRnBf0YKZ48xWDUB4a0sliEUm32oeGslHWHRpqsX1b0CmIiOSMUBYC+3x30CmIiOSMUBaC/V1OCzoFEZGcEcpC0HL9mqBTEJEsMmvWdCorX+e11xYycGBXZs2aHnRKoRLKQrD/pIKgUxCRLDFr1nQmT47wxRd7Adi0aQOTJ0dUDGoJZSHwo1sHnYKIZIm7757C53XOO37++W7uvntKQBmFTygLQat1q4JOQUSyRE3N+41anotCWQi+7NU36BREJEsUFMR/0FV9y3NRKAtBi21bgk5BRLLELbeU0aZN20OWtWnTlltuKQsoo/AJZSGwL78MOgURyRJjxoznnnvKOeqoowHo3Pk07rmnnDFjxgecWXik/MziprC/oEvQKYhIFhkzZjwzZjwEZO8UE6kI5R5BS50sFhFpNqEsBLqzWESk+YSyEHiLUKYlIpKVUhpxzWyUma0xs3fM7JY48aPNbGYs/oaZdU1kuwceVSkiIk0v6UJgZnnAA8CFQC9gnJn1qtPsB8DH7v4vwH8C/5HItvf1ODPZtEREpJFS2SMYCLzj7u+6+xfAH4FL67S5FHgi9roCOM/MrKEN522pSSEtERFpjFQuH+0M/KPW+43A2fW1cfd9ZrYD6AB8WHdjZhYBIrG3ezt3thUp5BZmJxCn/1lE/ctsWd+/zp0tW/vXPdkVQ3MfgbuXA+UAZrY02UeuhV029w3Uv0yn/mUuM0v6+b6pHBraBJxS632X2LK4bcysJdAe2J7CZ4qISJqlUgiWAGeYWTczOwq4AnimTptngAmx12OBl9zdU/hMERFJs6QPDcWO+d8APA/kAY+6+0ozuwtY6u7PAI8A/2Nm7wAfES0WiShPNq8MkM19A/Uv06l/mSvpvpn+QBcRyW26hVdEJMepEIiI5LjACkEC01NcY2bbzKwq9nNtEHkmy8weNbOtZvHvh7Co/4r1f7mZDWjuHJOVQN9KzWxHre/utubOMRVmdoqZzTezt81spZn9OE6bTP7+EulfRn6HZtbazN40s+pY3+6M0yapqW/CIMH+NX7sdPdm/yF6cnk98HXgKKAa6FWnzTXA/UHkl6Y+DgMGACvqiV8EzAMMGAS8EXTOaexbKfBs0Hmm0L9OwIDY63xgbZx/n5n8/SXSv4z8DmPfR7vY61bAG8CgOm1+BDwYe30FMDPovNPcv0aPnUHtESQyPUVGc/eXiV4pVZ9Lgd971OvAsWbWqXmyS00Cfcto7r7Z3ZfFXu8EVhG9S762TP7+EulfRop9H7tib1vFfupeEZPU1DdhkGD/Gi2oQhBveop4/xAvj+12V5jZKXHimSzR/waZ6pzY7us8M+sddDLJih026E/0L6/asuL7O0L/IEO/QzPLM7MqYCvwgrvX+925+z7gwNQ3GSGB/kEjx84wnyyeA3R1977AC/yzgkv4LQNOc/ci4L+B2cGmkxwzawf8GfiJu38adD7p1kD/MvY7dPf97t6P6GwHA80sq6YzTqB/jR47gyoEDU5P4e7b3X1v7O3DQHEz5dZcEpmiIyO5+6cHdl/dfS7QysxOCDitRjGzVkQHyenuPitOk4z+/hrqXzZ8h+7+CTAfGFUnlBVT39TXv2TGzqAKQYPTU9Q53noJ0eOY2eQZ4OrY1SeDgB3uvjnopNLBzE4+cMzVzAYS/XeWMf+jxXJ/BFjl7r+sp1nGfn+J9C9Tv0Mz62hmx8ZetwFGAKvrNMvYqW8S6V8yY2cgs496YtNT3GhmlwD7iJ6YvCaIXJNlZk8SvfLiBDPbCNxO9MQO7v4gMJfolSfvALuB7wWTaeMl0LexwEQz2wd8DlyRKf+jxQwB/hV4K3YsFuDnwKmQ+d8fifUvU7/DTsATFn1wVgvgT+7+rKVn6pswSKR/jR47NcWEiEiOC/PJYhERaQYqBCIiOU6FQEQkx6kQiIjkOBUCEZEcp0IgIpLjVAhERHLc/wd4e3t3r11UFQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# plot it\n", "x = np.array([1.,2.,3.])\n", "b = np.array([0.5, 1.5, 2.5, 3.5, 4.5])\n", "ax = plt.subplot(111)\n", "ax.hist(x, b, weights=bg, alpha=0.1, color='blue', label='BG', histtype='stepfilled')\n", "\n", "n1,b1,p1 = ax.hist(x, b, weights=bg+err, histtype='stepfilled',\n", " facecolor='none', edgecolor='none', linestyle='--')\n", "n2,b2,p2 = ax.hist(x, b, weights=bg-err, histtype='stepfilled', \n", " facecolor='none', edgecolor='none', linestyle='--')\n", "ax.bar(x=b1[:-1], height=n2-n1, bottom=n1, width=np.diff(b1),\n", " align='edge', linewidth=0, edgecolor='blue', color='none', zorder=-1, \n", " label='bg uncertainty', hatch='\\\\\\\\\\\\')\n", "\n", "ax.hist(x, b, weights=sig, color='red', label='sig $\\mu$=1', histtype='step',\n", " linestyle='dotted')\n", "ax.errorbar(x, data, yerr=np.sqrt(data), fmt='o', color='black', label='data')\n", "# ax.bar(x, sig, align='center',width=1, color='None', edgecolor='red', label='sig ($\\mu$=1)')\n", "\n", "ax.set_xlim((0.5, 3.5))\n", "ax.legend()\n", "_ = ax.set_ylim((0,20))" ] }, { "cell_type": "code", "execution_count": 4, "id": "ef904097", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAD4CAYAAAAEhuazAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABNn0lEQVR4nO2deXgURfrHP5WbkHAkhCuAhCsEAoRwqQGMCwKiggIqCCoiorj+1FV21fVYD8SbBVbURUDxRF0OkVNRooDcR7iPiAECCIKY+5pM/f6omRCSTDJJemYyM/V5nnl6uru6+judzjs11VXvV0gp0Wg0Go134eNqARqNRqNxPjr4azQajReig79Go9F4ITr4azQajReig79Go9F4IX6uOnGDBg1ku3btXHV6u8nOzqZu3bqullEpWqexuINOd9AIWqfR7Nix47yUMqKm9bgs+Ddp0oTt27e76vR2k5SURGJioqtlVIrWaSzuoNMdNILWaTRCiONG1KO7fTQajcYL0cFfo9FovBAd/DUajcYLcVmfv0ZjBIWFhaSlpZGXl2dovfXr1+fgwYOG1mk07qARtM7qEhQURIsWLfD393dI/Tr4a9yatLQ0QkNDad26NUIIw+rNzMwkNDTUsPocgTtoBK2zOkgpuXDhAmlpaURFRTnkHLrbR+PW5OXlER4ebmjg12hcjRCC8PBww3/RlkQHf43bowO/xhNx9H2tg79Go9F4ITr4azQajReig79Go9F4IZUGfyHEfCHEOSHEvkrK9RJCmIQQo4yTp9F4L6mpqcTGxtaojpCQEIPUuI4JEybQuHHjMtdi9erVREdH065dO1599VUXqbvEyZMnufbaa+nUqROdO3dm5syZ5ZarLbrtafl/CAypqIAQwhd4DfjWAE0ajUZTzPjx41m9evVl24qKivjrX//KqlWrOHDgAJ9//jkHDhxwkUKFn58fb731FgcOHGDz5s3Mnj27jKbapLvS4C+l/An4o5Ji/wcsAs4ZIUqjcSdSU1OJiYnhvvvuo3PnzgwaNIjc3NwyLfc333yT559/ntTUVDp27Mj48ePp0KEDY8eOZe3atSQkJNC+fXu2bt1afIzJZGLs2LHExMQwatQocnJyALj55pvp378/nTt3Zs6cOZVqnD59OrGxscTGxjJjxozi7S+99BLR0dH07duXMWPG8Oabb1b4mSojMTGRQ4cOAXDhwoUa/3IB6N+/P2FhYZdt27p1K+3ataNNmzYEBAQwevRovv766yrVm5ycTP/+/enUqRM+Pj7Uq1eP5557rto6mzVrRnx8PAChoaHExMRw6tQpw3UbRY0neQkhIoFbgGuBXpWUnQRMAoiIiCApKammp3c4WVlZWqeBGK2zfv36ZGZmAhD4xBP47N1rSL11pMQkBOYuXch/7bUKy2ZlZXH06FHmzp3L9OnTufvuu/nkk0+46qqrMJvNxfry8/PJz88nKyuLlJQUPvzwQ2bOnEliYiILFixg1apVrFy5khdffJHPP/+crKwsDh8+zH/+8x/ee+89HnzwQf7973/z8MMPM3PmTOrXr09BQQGJiYkMGjSI8PDwcvX99NNPzJs3j++//x4pJX/5y1/o2bMnJpOJr776ig0bNlBYWEi/fv2IjY0lMzPT5mcaPXp0hdfi6NGjNGvWjMzMTDZv3kzHjh0pKioqvgYAgwcPJisrq8yxU6dO5dprr7V5jUtey5SUFJo2bVq8Hh4ezvbt2y87T0Xk5eVx66238t///peePXvy0ksvkZuby5QpU2qsFeD48ePs3LmTTp06XVZfVXXn5eU57P/aiBm+M4AnpJTmysalSinnAHMAoqOjpTukT3WXNK/eqvPgwYOXZmUGBICvryH1moqK8PP1hYAAAiqZ9RkSEkJUVBQJCQkA9OnTh7NnzxISEoKPj0+xvsDAQAoLC4vLX3nllQB06dKFwYMHU69ePXr37s2rr75KaGgoISEhtGzZkuuuuw6Ae+65h1mzZhEaGspbb73FokWL8PHx4dSpU/z222+0bt26XH27du1i5MiRNG3aFIBRo0axc+dOzGYzt9xyCxERKjX88OHDCQwMLD53eZ+pohmwx48fp0WLFtSvXx9Qga5Hjx74+vpedtzPP/9c4fW0dY1LXss6derg7+9/2XpAQIDdM3R//PFHevbsWRzAe/XqxTfffEO9evUuK1cdrVlZWdx9993MnDmTyMjIy/ZVVXdQUBDdu3evsgZ7MCL49wQWWgJ/I2CoEMIkpVxqQN0ajf2U6M6oKblVnOofGBhY/N7X15fc3Fz8/Pwwm83F20vO1ixZ3sfHp3jdx8cHk8lUvK90g0oIQVJSEmvXrmXt2rU0adKExMREh8wELe8zVURycjJdu3YtXt+xYwe33357mXL9+vUrt6X75ptvMnDgQLu0RUZGcvLkyeL1tLS0MoG2Ivbt20eXLl2K13fu3Em3bt1qrLWwsJCRI0cyduxYRowYYbhuI6lx8JdSFieeEEJ8CCzXgV+jUYZF586d48KFC4SEhLB8+XKGDKlw7EQZTpw4waZNm7jqqqv47LPP6Nu3L+np6TRs2JDg4GAOHTrE5s2bK6yjX79+jB8/nieffBIpJUuWLOHjjz/GZDJx//3389RTT2EymVi+fDmTJk2q9ufdvXt38ZfQ0aNH+frrr5k6dWqZcuvXr6/2Oaz06tWLo0eP8uuvvxIZGcnChQv57LPP7D4+PDycH374AYAjR46wePFi1qxZUyOtUkruvfdeYmJieOyxxxyi20jsGer5ObAJiBZCpAkh7hVCPCCEeMDx8jQa98Xf35/nnnuO3r17c91119GxY8cq1xEdHc3s2bOJiYnh4sWLTJ48mSFDhmAymejZsydPPvlkcfeRLeLj4xk/fjy9e/emT58+TJw4ke7du9OrVy+GDRtG165duf766+nSpUtxl011SE5Oxmw2061bN1588UU6derEggULql2flTFjxnDVVVdx+PBhWrRowbx58/Dz8+Ptt99m8ODBxMTEcNttt9G5c+cq1ZmVlUVsbCyTJk3i888/t/nMxF42btzIxx9/zA8//EBcXBxxcXGsXLkSgKFDh3L69Oka6zYUKaVLXh06dJDuwLp161wtwS68VeeBAwcMrc9KRkaGQ+o1EiM0ZmZmSimlzM7Olj169JA7duyodl3t2rUrV5M7XEspa6fO8u5vYLs0IAbrlM4ajRczadIkDhw4QF5eHnfffXfxUMWqkpmZiRCi1qRE1lSODv4ajRdjVH9zaGgoR44cMaQujXPQuX00Go3GC9HBX6PRaLwQHfw1Go3GC9HBX6PRaLwQHfw1Go3GC9HBX6PRaLwQHfw1Go3GC9HBX6PRaLwQHfw1Gs1leJptItQ+7bUBHfw1Gs1leJJtItRO7bUBHfw1GgOwZYd4880306NHj8vsFm3ZO2ZnZ3PDDTfQrVs3YmNj+eKLLyo9rzvbJgohHG6baJR2T0Tn9tF4FuW5hN12Gzz4IOTkwNChZfePH69e58/DqFEA1CkqUq5gdljobdu2jUWLFpGcnExhYSHx8fH06NEDgPnz5xMWFkZubi69evVi5MiRNutZvXo1zZs3Z8WKFQCkp6dXeu6UlBQ6dOgAwJ49ey4zKLFihHnKqVOnaNmyZfF6ixYt2LJli13HgjKyuf322/noo4/o3bs3zz77LHl5ebzwwguG6E1NTWXXrl306dPHcO2eig7+Gk0N2bhxI8OHDycoKIigoCBuuumm4n2zZs1iyZIlgOqjPnr0aLGdYmm6dOnC448/zhNPPMGNN95Iv379KjzviRMniIyMxMdH/YDfs2fPZU5aVowwT6kpa9euJT4+nt69ewPQtWtXVq9eXcapDKquNysri5EjRzJjxowyNowa2+jgr/EsKmqpBwdXvL9Ro+L9VbVxLF+KslvctGkTwcHBxXaLtuwdO3TowM6dO1m5ciXPPPMMAwYMqLBbZN++fW5tm2grfXRV9FZmm2iEdk9FB3+NpoYkJCSUa4doy27Rlr3j6dOnCQsLY9y4cTRo0IC5c+dWeN49e/a4tW2iLXN0e/VKO2wTjdDuqegHvhpNDbFlh2i1W4yJibnMbtGWvePevXvp3bs3cXFxvPDCCzzzzDMVnnffvn3aNtGGbSLUUuvE2oQRdmDVeWkbR2PxVp21xcbRSDtEe2nTpk2ttB4sjTtolLJ26tQ2jhpNLccoO0R70baJmppSafAXQswHbgTOSSnLDCIWQowFngAEkAlMllImGy1Uo6nNOLsPOTQ0lF27djn1nBrPwp4+/w+BIRXs/xW4RkrZBXgJmGOALo1Go9E4kEpb/lLKn4QQrSvYX/KR/WaghQG6NBqNRuNAjO7zvxdYZWunEGISMAkgIiKCJDtmT7qarKwsrdNAjNZZv379cseE15SioiKH1Gsk7qARtM6akJeX57D/a8OCvxDiWlTw72urjJRyDpZuoejoaJlY3lT8WkZSUhJap3EYrfPgwYMOeeiZacAkL0fjDhpB66wJQUFBdO/e3SF1GxL8hRBdgbnA9VLKC0bUqdFoNBrHUeNJXkKIVsBi4E4p5ZGaS9JoNBqNo7FnqOfnQCLQSAiRBvwL8AeQUr4HPAeEA+9YkjSZpJQ9HSVYo9FoNDXHntE+YyrZPxGYaJgijUaj0TgcndtHo6mllDZ9qQ4hISEGqXEdrVu3pkuXLsTFxdGzZ+3oVLBldVmS2m4dqYO/RqOp9axbt47du3ezfft2V0sByre6LIk7WEfq4K/R1JDU1FRiYmK477776Ny5M4MGDSI3N9emXWNqaiodO3Zk/PjxdOjQgbFjx7J27VoSEhJo3749W7duLT7GZDIxduxYYmJiGDVqFDk5OYCyh+zfv/9l9pAVMX36dGJjY4mNjWXGjBnF223ZT9r6TJXhCFvJmmK0fSSUb3VZEnewjtSJ3TQew6OrH2X3b7sNqauoqAhfX1/imsYxY8iMSssfPXqUzz//nPfff5/bbruNRYsW0bevzSkvpKSk8NVXXzF//nx69erFZ599xoYNG1i2bBnTpk1j6dKlABw+fJh58+aRkJDAhAkTeOedd5gyZQrz58/H398fPz+/YntIWymSd+zYwQcffMCWLVuQUtKnTx+uueYaTCaTTftJW59p3LhxFV4HR9hKCiEYNGgQQgjuv/9+Jk2aVKGGkthrH9mvXz/S09OLXdEq01QZ7mAdqYO/RmMAUVFRxMXFAdCjRw9SU1MrDP5RUVHFgbFz584MGDAAIQRdunQhNTW1uFzLli1JSEgAYNy4ccyaNYspU6Ywa9YsFi1ahI+PT7E9pK3gv2HDBm655Rbq1q0LwIgRI1i/fj1ms9mm/aStz1QRx48fd4it5IYNG4iMjOTcuXPF/gf9+/e361h77SPXr19fKyd5ORId/DUegz0tdHupaiAIDAwsfu/r60tubq5Nu8bS5X18fIrXfXx8MJlMxftKBykhRLE95Nq1a2nSpEmxPaTRlPeZKiI5OdkhtpJWy8XGjRtzyy23sHXrVruDv732kUa3/N3BOtIr+/y/+AI6d4b5812tROPJlLRrzM/PZ/ny5VWu48SJE2zatAlQaaP79u1r0x7SFv369WPp0qXk5OSQnZ3NkiVL6NevHwkJCXzzzTfk5eWRlZVVLX0l2b17dxlbyfK6fdavX8/u3bvLvMoLstnZ2cVfFNnZ2Xz77bdVeo4QHh7Onj17gEv2kaNHjy5X08aNG+3SZA8lrSMLCgpYuHAhw4YNq1ZdjsIrg//8+XDgAPztb1CikaXRGIotu8aqEB0dzezZs4mJieHixYtMnjy52B6yZ8+el9lD2iI+Pp7x48fTu3dv+vTpw8SJE+nevbtN+8nqkpycbLit5NmzZ+nbty/dunWjd+/e3HDDDQwZUlGG+ctxhH2ktd7SVpfgZtaRRtiBVeflKhvHnBwpg4KkbNpUSpBy27aKy3urPaKj8FQbR1dghEYj7SfbtWtXriZ3uJZS1k6djrRx9LqW/4YNkJcH06apdTfIgqzROIxJkyYRFxdHfHw8I0eOrLb9pLaVdD+87oHv/v1qedNN0K4d/PxzxeU1Gk/GKPvJ0NBQjhzReR3dCa9r+R87BqGhEB4OMTGQkuJqRRqNRuN8vC74//orREWBEKrl/8svIKWrVWk0Go1z8drgDyr45+TAb7+5VpNGo9E4G68K/lKq4N+mjVpv104tddePRqPxNrwq+J87p1r61pZ/27ZqqYO/RqPxNrwq+J84oZZXXHFp6eurHgJrNBqNN+FVwf/sWbVs2lQt/fzU+1OnXKdJo9FoXIFXBf9z59SyceNL25o3h9OnXaNHo9FoXIVXBv+IiEvbmjfXLX+NpiTaNtE7qDT4CyHmCyHOCSH22dgvhBCzhBApQog9QojqzQ93AufOQd266mVFt/w1mrJo20TPx56W/4dARWn0rgfaW16TgHdrLssxnDt3eZcPQGQk/PGHyvej0VQXW3aIN998Mz169LjMbtGWvWN2djY33HAD3bp1IzY2li+++KLS87qDbWK9evW0bWItpNLcPlLKn4QQrSsoMhz4yJJtbrMQooEQopmU8oxRIo3i3Dlo0uTS+o4daqYvqNa/dfy/xn1J/DCxzLbbOt/Gg70eJKcwh6GfDi2zf3zceMbHjed8znlGfTkKuGTjmDQ+qdJzbtu2zaYd4vz58wkLCyM3N7fYbtEWq1evpnnz5qxYsQKA9PT0Ss/tDraJGRkZ5domVkVTZbiDbWJtw4jEbpHAyRLraZZtZYK/EGIS6tcBERERJDk5peaxYz1p0iSPpCTVg/Xdd014//0oIIgVK3bRpUvZf7asrCyn66wO3qqzfv36lwWRoqKiMmXy8vLIzMwkpzCnwv1ZuVnF+6WUFBUVlRugSvP9998zZMgQCgsLARg8eDD5+flkZmbyxhtvFJuknDx5kt27d9OkSRPMZnNx3fn5+eTn5xMVFcW3337L3/72N4YMGcLVV19d4fl//fVXmjZtSnZ2NqBav9HR0WWOWblypc06yqvf+iX0+++/M3z4cFq1alVsJVkZq1atokuXLsTExJCZmUn79u357rvvyMrKqpEmUPdOyetWktzcXAoLC4v35ebmUlBQYNffz4q9f29nkpeX57D/a6dm9ZRSzgHmAERHR8vExERnnp6cHOjUKQTreU+fht9/V/uaNOlOeXKSkpJwts7q4K06Dx48eFka4fX32vaHDSW04v2hl/ZXxcYxKCiIwMDA4vIBAQEEBgayY8cO1q9fz5YtWwgODiYxMRFfX18aNGhQfD5QXzSBgYHEx8eza9cuVq5cybRp0xgwYECF3SUHDx6ke/fuxfXs37+f22+/vYzuqrayo6Oji/WNHDmSffv22W2gcuzYMeLj44s1HDx4kLi4uBprAggJCcHHx6fcv0u7du345JNPivdduHCBqKioKqWYro0evkFBQXTv3t0hdRsR/E8BLUust7Bsq1WYzSrQN24MW7bAK6/Am29CWJjq87eOBNJoqkpCQgL3338/Tz31FCaTieXLlzNp0iSbdosl7R1DQkJYvnw5Q4YM4fTp04SFhTFu3DgaNGjA3LlzKzzvnj17ytgmTp06tUy5qhimZ2dnYzabCQ0NLbZNrEp/fXh4OD/88ANwyTZxzZo1NdJkDyVtEyMjI1m4cKFh6ao9FSOGei4D7rKM+rkSSK+N/f0ZGcqyMTwcli9Xr4gI6N1b7T9/3rX6NO6LLTtEq91iTEzMZXaLtuwd9+7dS+/evYmLi+OFF17gmWeeqfC8+/bt07aJ7mSbWNuozOoL+BzVf1+I6s+/F3gAeMCyXwCzgV+AvUBPeyzEnG3j+Msvyrbxgw+kHDZMyk6d1Pbnn1fbJ04s/zhvtUd0FJ5q42ikHaK9tGnTplZaD5bGHTRKWTt1OtLG0Z7RPmMq2S+Bv1b/68c5XLyolg0bwt690KuXWr/lFnjnHd3y19SMSZMmceDAAfLy8rj77rurbYdoL9o2UVNTvMbG0Rr8AwNVWucJE9R6167Qvj3YMapOo7GJs/uXQ0ND2bVrl1PPqfEsvCa9gzX4A/TpAyUbZoGBcPJk2WM0Go3GU/G6ln9sLFgGXRSzezeUGoas0Wg0Ho3XtfwbNiy7LyICCgrUcFCNRqPxBrwq+Pv7w1NPwQ03XL4vMlIt//jD+bo0Go3GFXhV8LeO9Cn9cNdq65ic7HxdGo1G4wq8LvinpkLr1pfvsxq5HzzobFUajUbjGrwq+DdooEb1WFv6Vvr0Ucv69Z0uS6PRaFyC1wT/P/6AOnWgqKhsy79tW7XMz3e6LI1Go3EJXhP8//wT6tWDUaOgdMpzq63j1q1Ol6XR2KS06Ut1CAkJMUiNazh58iTXXnstnTp1onPnzsycOdPVkgD7LSOLioro3r07N954oxPV2YfXjPNPT1ejet55p+y+OnWUqYslGaFGo6kl+Pn58dZbbxEfH09mZiY9evTguuuuo1OnTi7TZLWM/O6772jRokVxYr/yNM2cOZOYmBgyMjJcoLRivKLlL6UK/vXq2S4TGHj5LGCNxl5SU1OJiYnhvvvuo3PnzgwaNIjc3Fybdo2pqal07NiR8ePH06FDB8aOHcvatWtJSEigffv2bC3xE9RkMjF27FhiYmIYNWoUOTk5gLKH7N+//2X2kBUxffp0YmNjiY2NZcaMGcXbbdlP2vpMlWG0rWSzZs2K8ySFhoYSExPDqVNVyxhf2lZSCFEjW0l7LSPT0tJYsWIFEydOrPa5HIlXtPzz8lQ6502boEUL9dDXat9oJSREpX3WuC+PPqpmaxtBUVEdfH0hLg5KxEqbHD16lM8//5z333+f2267jUWLFtG3b1+b5VNSUvjqq6+YP38+vXr14rPPPmPDhg0sW7aMadOmsXTpUgAOHz7MvHnzSEhIYMKECbzzzjtMmTKF+fPn4+/vj5+fX7E9pK3UyTt27OCDDz5gy5YtSCnp06cP11xzDSaTyab9pK3PNG7cuAqvgyNsJa2kpqaya9cu+lhHaNhBebaSeXl55dpKpqen4+NzeXu4PE32WkY++uijvP7667XOHcyKVwR/67j+3FwICCgb+EGNBDp/Xs30DQhwqjyNBxAVFUVcXBwAPXr0IDU1tcLgHxUVVRwYO3fuzIABAxBC0KVLF1JTU4vLtWzZsthCcdy4ccyaNYspU6Ywa9YsFi1ahI+PDydPnuTo0aM2g/+GDRu45ZZbqFu3LgAjRoxg/fr1mM1mhg8fTlBQEEFBQdx0002VfqaKOH78OJGRkcUBdM+ePXTt2rVMueoYuWRlZTFy5EhmzJhBvYp+wpdi7dq1xMfH09ti3NG1a1dWr16NKBUE1q9fb6iT1/Lly2ncuDE9evSotfaqXhH8rS367Gxo3rz8MuHhkJKivgBsldHUbuxpodtLZmZulQJBYGBg8XtfX19yc3Px8/PDXCJniNV1q3R5Hx+f4nUfHx9MJlPxvtJBSghBUlISa9euZe3atTRp0oTExMTL6jaK8j5TRSQnJ18W7Hfs2MHtt99eplxVW/6FhYWMHDmSsWPHMmLEiKp8BPbt23fZr4+dO3eWm267Ki3/yMhITpbIBJmWlkakNU2AhY0bN7Js2TJWrlxJXl4eGRkZjBs3jk8++aRK+h2JV/T5W1v+6enQrFn5Zbp3V/3+1pE/Gk1NKWnXmJ+fX2zkXhVOnDjBpk2bAJU2um/fvjbtIW3Rr18/li5dSk5ODtnZ2SxZsoR+/fqRkJDAN998Q15eHllZWdXSV5Ldu3eXsZUsr9tn/fr17N69u8yrvMAvpeTee+8lJiaGxx57rMqawsPD2bNnD3DJVnL06NHlatq4caNdmkpaRhYUFLBw4UKGDRt2WZlXXnmFtLQ0UlNTWbhwIX/5y19qVeAHL2v5X7xou1V/xRVqnH9hocoBpNHUlJJ2jZGRkcV2jVUhOjqa2bNnM2HCBDp16sTkyZPx9fXlvffeo2fPnsTExBTbQ9oiPj6e8ePHF3d9TJw4sdgU3Go/2aRJk2L7yeqSnJxMUFAQ3bp1o2vXrsW2ko8++mi169y4cSMff/wxXbp0Ke6CmjZtGkOHDrXr+DFjxrBs2TJiY2Np1KiRIbaSJS0ji4qKmDBhQrFl5NChQ5k7dy7N3aH7wAg7sOq8nGnjuGiRsmq84w4pv/66/DL/+Y8qs3r15du91R7RUXiqjaMrMEKjkfaT7dq1K1eTO1xLKWunTpfaOHoC1pb/1KllUztYadBALdesgcGDnSJLo3E5RtlPaltJ98Mrgr+1z7+i7hzr85oqDiHWaNwao+wnQ0NDOXLkiCF1aZyDXQ98hRBDhBCHhRApQogny9nfSgixTgixSwixRwhhX4eck7C2/Fu2VGP9y6NRI7X87TfnaNJoNBpXUmnwF0L4ArOB64FOwBghROl5zM8AX0opuwOjgXKSKLiO9PRLrX5bz3qs28+fd44mjUajcSX2tPx7AylSymNSygJgITC8VBkJWGde1AdOGyex5qSnQ1CQem9rKKc1+OsUDxqNxhuwp88/EjhZYj0NKD2/+nngWyHE/wF1gXLnaAshJgGTACIiIpw28y0lpRPQEB8fP3bv/rHcGb4AgYH96NXrDElJKcXbsrKyau0MvZJ4q8769es7ZPp8UVFRrZ2Wb8UdNILWWRPy8vIc9n9t1APfMcCHUsq3hBBXAR8LIWKllJdZoksp5wBzAKKjo2ViYqJBp6+Y115TE7iCg+Haa22fs3FjaNCgBYmJLYq3JSUl4SydNcFbdR48eNAhI0yMnOrvKNxBI2idNSEoKKh4TobR2BP8TwEtS6y3sGwryb3AEAAp5SYhRBDQCDhnhMiakpGhZvbefXfF5erUgZ9/hsxMqGX3gEaj0RiKPX3+24D2QogoIUQA6oHuslJlTgADAIQQMUAQ8LuRQmtCZia0aQOPP15xOX9/ld9HD/fUaDSeTqXBX0ppAh4C1gAHUaN69gshXhRCWBNaPA7cJ4RIBj4HxltmotUKMjPBxweysiouZx3uea5W/F7RaDQax2HXOH8p5UopZQcpZVsp5cuWbc9JKZdZ3h+QUiZIKbtJKeOklN86UnRVycyEVavgH/+ouFyTJmp59qzjNWk0tRFtm+g9eEVWz4wMZehSWcZO6yxfHfw13orVNvHAgQNs3ryZ2bNnc+DAAZdqstomrlq1igMHDvD555/b1GS1TdRUjscH/4IClakTLnXr2KKFZZBPWppjNWk8D1t2iDfffDM9evS4zG7Rlr1jdnY2N9xwA926dSM2NpYvvvii0vO6g21ivXr1tG1iLcTjc/uUHLZbWcvf+uUwYYLj9GgcS3mjSG+7DR58EHJyoLxMwOPHq9f58zBqlNpmtXG0Z4j1tm3bbNohzp8/n7CwMHJzc4vtFm2xevVqmjdvzooVKwBItyalqgB3sE3MyMgo1zbRXk2eYptY2/Cq4F9Zy1/P8tVUh40bN9q0Q5w1axZLliwBKLZbbNq0abn1dOnShccff5wnnniCG2+8kX79+lV43hMnTriFbeI333xTrm2ikbiDbWJtw2uC//jx0Kl0RqJSWIP/7NlQhcaOphZR0f99cHDF+xs1urS/qjaO5WtRdoubNm0iODi42G7Rlr1jhw4d2LlzJytXruSZZ55hwIABFXaX7Nu3zy1sE7t161YjTZ5im1jb8Jrgf/vtlXvzWoP/2rWO1aTxLBISErj//vt56qmnMJlMLF++nEmTJtm0Wyxp7xgSEsLy5csZMmQIp0+fJiwsjHHjxtGgQQPmzp1b4Xn37NlTxjZx6tSpZcpVpZUtDbBN/OGHH4BLtolr1qypkaaStomRkZEsXLiwTCrqV155hVdeeQVQX7pvvvmmDvyV4DXB355sndbgb00BrdHYQ69evcq1QxwyZAjvvfceMTExREdHF9st2rJ33Lt3L3//+9/x8fHB39+fd999t8Lz7tu3j5CQkDK2ic8++2y1P4u2TfQijLADq87LWTaOX32l7Bnr1q28bFGRKgtSmkxqm7faIzoKT7VxNNIO0V7atGlTK60HS+MOGqWsnTq1jWMNsLb8w8IqL+vjo/qFc3LUQ9/KHhBrNFaMskO0F22bqKkpXhP8KxvmaaVePTUv4M8/dfDX2I9Rdoj2Ehoayq5du5x6To1n4TXBv3Fj+8pfcQXExkK7do7TpNFoNK7G42f4WoO/NW9PZTRqBBcuOE6PRqPR1Aa8IviHhMCkSfaVDw2Fw4fBMslS4wbI2pNAVqMxDEff114R/MPD4eqr7SsfEaEe+CYnO1aXxhiCgoK4cOGC/gLQeBRSSi5cuECQ1XzcAXh8n//FiyCEGudvzwNca/fQb785VpfGGFq0aEFaWhq//26sd1BeXp5D//GMwB00gtZZXYKCgmjRokXlBauJxwf/s2chNRU2boThwysvbx0Sqt283AN/f3+ioqIMrzcpKclh3qlG4Q4aQeusrXh8t481MaK9Qz2tkxF1Tn+NRuPJeHzwt1o32jtm3xr8/Tz+N5FGo/FmPD74Z2erZVVb/g8/7Bg9Go1GUxvw+OCfl6ce+Navb195a/DXY/01Go0nY1fwF0IMEUIcFkKkCCGetFHmNiHEASHEfiGEc+e620BKZeN4660qb489WIP/K69Abq7jtGk0Go0rqbRnWwjhC8wGrgPSgG1CiGVSygMlyrQHngISpJQXhRB2JlNwLLm56gvAkpnWLoKDVX//r7/alwZao9Fo3BF72sO9gRQp5TEpZQGwECg9aPI+YLaU8iKAlPKcsTKrh/Vhb1WDuDVRog7+Go3GU7FnTEskcLLEehpQ2uSwA4AQYiPgCzwvpVxduiIhxCRgEkBERITDvTZPnQoCrmTFigxuummn3cfVqdObixeD+f77ZDp2zHILT9CsLK3TSNxBpztoBK2ztmLUgEY/oD2QCLQAfhJCdJFS/lmykJRyDjAHIDo6WiYmJhp0+vLZvVstr7iiHlU5V2QknD4NzZt3IyTkYpWOdRVJSUlap4G4g0530AhaZ23Fnm6fU0DLEustLNtKkgYsk1IWSil/BY6gvgxcyp9/qqW9wzytNGkCgYFQi2Z6azQajaHYE/y3Ae2FEFFCiABgNLCsVJmlqFY/QohGqG6gY8bJrB5nzqilvbn8rTRvroaGjhhhvCaNRqOpDVQa/KWUJuAhYA1wEPhSSrlfCPGiEGKYpdga4IIQ4gCwDvi7lNLlI+Wtwb9Zs6odFx4Of/yhRgppNBqNJ2JXn7+UciWwstS250q8l8BjlletoU4dtbzuuqodFx4OJhM88ACMGWO8Lo1Go3E1Hj3D1zpJq3Xrqh1nzeypc/prNBpPxaOD/5EjaunvX7XjrLN8DU4Rr9FoNLUGjw7++/appa9v1Y6zBv+LF43Vo9FoNLUFjw7+6ekqqVtVh2xag396un7oq9FoPBOPDv6ZmVVv9cOleQEdO0JBgUdfIo1G46V4dGTLzoaAgKof17Ch+tK4+WYIDDQbrkuj0WhcjUcH/7y86s3S9fFRrX/9wFej0XgqHh38W7WCnj2rd2xoKHzyCWzeHGasKI1Go6kFeHTwz8mpemoHKxERap7AH39Uo99Io9FoajkeG/zz8uDsWbWsDs2bq2V6ehUnCWg0Go0b4LHB/9w51XK3GrhXFR38NRqNJ+Oxwf/0abUMq2aXvbW7SHf7aDQaT8Rjg//x42pZ3T5/63GRkdrFXaPReB4eG/xPWownmzSp3vHW4H/llS7PTK3RaDSG47HB39rtExlZveOtwf/PP3W3j0aj8Tw8NvgPs9jMNG1aveOtwf+FFzoZI8hATmWcYsHuBbyz7R22pG1B6gREGo2mihhl4F7rsI7yqVevesdb8/vk5vpRWFj1tNCOoKCogGd/eJbpm6djMpuKt1/V4ioeaPaAC5VpNBp3w2Nb/l9+qZbVDf6hoeBn+Wo8f94YTTWhoKiAUV+O4vWfX+eurnexb/I+0v6Wxrs3vMuRC0eYvHMyP6b+6GqZGo3GTfDY4P/TT2pZ3eAvhDJxh9oR/B9d/SjfHPmG2UNnM2/4PDo37kxkvUge6PkAux/YTURgBDd9fhMHfj/gaqkajcYN8Njgn56ulqGh1a/Dmtff1cF/6aGlvLv9XaZcNYUHez1YZn+Lei14vcvrBPsHM3zhcC7mahcajUZTMXYFfyHEECHEYSFEihDiyQrKjRRCSCFENdOpGYOUKpc/QN261a+neXNo0KCg2nMFjCC7IJsHVzxIXNM4pg2YZrNc46DGLL59Mcf/PM4DK3T/v0ajqZhKg78QwheYDVwPdALGCCHKDIERQoQCjwBbjBZZVbKzwWRSufx9avDbplUrlc+/c2fjtFWVmVtmcibrDG9f/zb+vhU/db665dU8d81zfLn/S5YdXuYkhRqNxh2xZ7RPbyBFSnkMQAixEBgOlO5cfgl4Dfi7oQqrwR9/KDOW4OCa1dO4MVy86E9WFoSEGKOtKlzIucBrG19jWPQwElollNl/Nusszyc9z8HzB8lMz2Rg4UDGdBnDl/u/ZPKKySS2TqReYDUfemg0Go/GnuAfCZwssZ4G9ClZQAgRD7SUUq4QQtgM/kKIScAkgIiICJKSkqos2F4SEjqRmlqXpKRt1a4jM7MlBQVtGTPmNI8/fsRAdfYx59gcsvKzuDn05nKvVX5RPov2LaJJYBMKTAVM3zSd139+nbGtxvLpuU958LMHmRg10em6KyIrK8uhf3ejcAed7qARtM5ai5SywhcwCphbYv1O4O0S6z5AEtDasp4E9Kys3g4dOkhHMniwlL1716yOTz6REqQcOtQYTVUhKz9LNni1gbz1y1sv277t1DY55JMhMis/S0opZZG5SEop5bp16+Tv2b/LZ75/Rp7PPi/vWHSHDJoaJE+mn3S69opYt26dqyXYhTvodAeNUmqdRgNsl5XEV3te9vSInwJallhvYdlmJRSIBZKEEKnAlcAyVz70XbMGdu6EOnVqVo81rfNvv9VcU1X5KPkj/sz7k0f6PFK8bd+5fQz6eBCHzh/itywlykdc+hM2Cm7ES395ifDgcJ7r/xwFRQX8/TuX98JpNJpaiD3BfxvQXggRJYQIAEYDxU8TpZTpUspGUsrWUsrWwGZgmJRyu0MU28HOncp/t7pj/K1Y8wI5e6inWZqZtXUWPZr14OqWVwOQlpHGoI8HUce/Dj/c9QNtw9pWWEeeKQ+BYOG+hXrsv0ajKUOlwV9KaQIeAtYAB4EvpZT7hRAvCiGGOVpgdTh3Tk3Sqm4ufyvWlv+ff9ZYUpX44dcfOHT+EI/0eQQhBCaziTsW3UFGfgZrxq0hqmFUpXV0a9qN/1z/HwBGfDHC0ZI1Go2bYddASCnlSillByllWynly5Ztz0kpy4wnlFImurLVD5eCf01b/iEhEBhYRJcuxuiylw93f0j9wPrc2vlWAM5knuFM1hneu/E9YhvH2l3P5F6TiW8Wz+ELh5m/a76j5Go0GjfEI2f4/v47mM01m91rpWnTvGpnBq0OGfkZLD64mDGxYwjyCwKgZf2W7J28l3Fdx1W5vv/d+j8EginfTsEszUbL1Wg0bopHBn9rQraatvwBwsIKOHas+l7AVeWr/V+Ra8plfNx4isxFTN80nayCrOIvgqoS1TCK0bGjyczPJC0jzWC1Go3GXfHI4D/f0sNhRPD39zeza9elRHGOZkHyAqLDo+kd2Zt3tr3D498+zuqU1TWq89WBr4KANze+ybnscwYp1Wg07oxHBv+LlrxmDRvWvK6mTZWH76lTlRQ0gFMZp1h/Yj13dr2T9Px0nv/xeQa2GcjImJE1qrdV/VaM6zqO2dtnM+jjQbr7R6PReF7w/+03GGfpGm/QoOb1tWiRB0BKSs3rqoylh5YCMLLTSF5Z/woXcy/yxnVvIISocd3/uPofmKWZ5LPJfLHvixrXp9Fo3BuPC/4nT6px/mBM8G/aVAX/Y8dqXldlLD60mE4RnQj2D2bmlpnc2e1O4prGGVJ3TEQMN3W4CV/hy9M/PE1BUYEh9Wo0GvfE44J/ye4ZI4J/eHg+AGkOflZ6Puc8P6b+yIiOIygsKmRQ20G8dO1Lhp7jyb5PUiSL+PXPX5m3c56hdWs0GvfCo4O/EX3+ERGqhdyrV83rqohlh5dRJIsYETOCtmFtWTZmGa3qtzL0HFe3vJq+LfsS6BvIguQF2vhdo/FiPC74nz59KYe/1YaxJoSFFSCEMb8iKmLxwcW0btCarae28ssfvzjsPE/2fZL8onwmxU8y5FmCRqNxTzwu+Nerp3LyBAWpV03x9ZVERKjnCI5qKGfkZ/Ddse9IvCKRySsm8/7O9x1zImBo+6HENo5l+ubp5BbmkmfKc9i5NBpN7cXjgv8TT8CQIca21AMCYPlyZRLjCFYcWUFBUQFnss5Qx78Oj131mGNOBAgheCLhCfb/vp+W/27JO9vecdi5NBpN7cXjgj+oRGxGBn9rgrczZ4yrsySLDy2mUXAjvjv2HZN7TqZxXceaBt/e+XZa1W+FyWzijZ/f0K1/jcYL8bjg36sX7NljbPBv00YtHTHRK7cwl5VHV9KoTiMCfAOYcvUU409SCn9ffx6/6nHS89P5Les3nfRNo/FCPCr45+TA9u2QlWVs8I+JUcsDDkiL/+0v35JTmEPHRh15qNdDNA1xTha5e7vfS1hQGA2DGvLqhlf1uH+NxsvwqOB/+rRamkzGBv+4OLU8eNC4Oq0sPrSYhkEN+fLWL3lj0BvGn8AGdQPq8nCfh7mYd5GTGSf5+tDXTju3RqNxPR4Z/HNzjRnjb6VzZ7W0dv8YRUFRAUsPLeXKFlfi7+tvbOV28FDvh6jjV4fBbQczqtMop59fo9G4Do8K/tY++cxMCA83rt5WrcDXV3UnGUlSahIZ+Rl8+8u3xZ68ziQ8OJz74u/j+1+/52TGST3pS6PxIjwq+DdoAP37q/H4jRoZV6+/PzRtChs2GFcnwCd7PgFgXNdxTuvrL411WOn4JePpM7cPReYil+jQaDTOxc/VAozk+uuhXTvo0MHYlr+VjRuNq6vIXMSig4sQCJ7p/4xxFVeRKxpcwZjYMXy5/0vyi/JZdHARt3W+zWV6HI6Uaizwn38qhx6TSfl1hoaqFoOvr6sVajROwaOCP8D582ppdPBv3lx1K2VnQ926Na9v5dGV5BTm0LdlX9qFtat5hTXgHwn/4OM9H9MouBFTf5rKqE6j8BEe8qPwxAnlxLN+vZqmnZKiAn95BARAVBRER9MqIkJ9MVx9NQQHO1WyRuMM7Ar+QoghwEzAF5grpXy11P7HgImACfgdmCClPG6w1kq5+upLvr1GdvsAtG0L27ap4Z5GJHlbkLwAgOmDp9e8shoS2ziWGzvcSFJqEnvP7WXZ4WXc3PFmV8uqPqdOwaefwldfqbG/oBI99eoFd9yhntyHh6tvcT8/9TAnI0N9UfzyC+zfT5tly2DePJUjZOBAuOkmGDYMpxo6azQOpNLgL4TwBWYD1wFpwDYhxDIpZclR77uAnlLKHCHEZOB14HZHCK6IgwchPl69N7rlbx3xs21bzYO/lJJtp7dxfbvr6RXp4HShdvJEwhMsP7KciOAIpv40leHRw90v8dvPP8PMmbBoERQVQc+e8NprMHgwxMZWqUtnw/Ll9PXzg1WrYNkyld9j8mQYOhQmTIAbblC/FDQaN8Weln9vIEVKeQxACLEQGA4UB38p5boS5TcD44wUaQ/Wbtw6ddS60S1/65fKnj01r2vxwcWcSD/BC4kv1Lwyg+jbqi9Xt7yaX/74hXdveNe9Av/mzfDss7B2rXrq/+ij8MAD6gFQeZjNKvVrUZE67sQJlbjpwgVIT4fRozElJsI118C//w0JCer96dOwaZP6IoiIgAcfhL/+Vb3XaNwMe4J/JHCyxHoa0KeC8vcCq8rbIYSYBEwCiIiIICkpyT6VdpCSEgL0JD39d3x9w9mx4yeMiF9ZWVkkJSVRVOQL9CM7O5WkpNRq15dblMvYn8cC0PD3hoZdA6vOmnBDvRt4+uTTLF6/mOym2YboKo0ROq0Enj1Lu3feIeKnnyho0IATkydz+qabMNepo9x3LA48fhkZNNyxg/p799Jgzx6yo6I4+PTTAFz1/vuY/fww1atHYb16mJo25WJ2NllZWfy0Zg3dTp0iYM8eAi9cwKewEIBz11yDb34+4S+8gPnllzl77bWcuPtuciMjDflc9mLktXQkWmctRUpZ4QsYhernt67fCbxto+w4VMs/sLJ6O3ToII1kyRIpQcqbb5ayaVPj6l23bl3x+xYtpBw7tmb1vbbhNcnzyF5zetWsolKU1FldisxFsvt73WXL6S3lmP+Nkd+mfFtzYaUwQqfMz5fy1VelDA6Wsk4dKV98UcrMzPLLjh8vpa+vujnq1pVy4EApZ826tL+oyD6dZrOUaWlSrlwp5f79atuCBape66trV7XfSRhyLZ2A1mkswHZZSXy152VPy/8U0LLEegvLtssQQgwEngaukVLm1+D7qFo0awZ33aV+vTtimCeoET8//VT94zPyM5i2fhoAk3pMMkiVcfgIH14b+BqDPhnEqqOrOJF+goFtBtauLqAdO9Qf+sABuPlm1S3TurXal52t+vuXLYOFC9XD3G7d1B/uppugRw81aaMkPnaOahJCGUWUbN2PG6dyfyxerB4O79mjngncdBO8/baaHajR1FLsufO3Ae2FEFFCiABgNLCsZAEhRHfgv8AwKeU542VWTp8+sGCB6rI1ur/fSmGhMog3map3/PRN00nPT8dX+DIiZoSx4gziurbXcV2b6ygwF7Dx5EZWp6x2tSSFyQRTp8KVV6qHO998A0uWqMB/9Cg88ohqAdx9N+zaBamp6rhHH4WXX1bHlQ78NcXHB7p2heefVzfGrl3wl7/A6tXQvr0aGTBypGox6NnTmlpGpcFfSmkCHgLWAAeBL6WU+4UQLwohhlmKvQGEAF8JIXYLIZbZqM5hXLyo/r9++81xo/Gs2T23bq36sVJKVhxdQbB/MIPaDiKsTpix4gzktYGvkVOYQ4OgBjz1/VOYpdm1glJToV8/9XB21CjYuxduvFHt274doqPh3Xdh+HAVaFNSbD/sdSRxcfD992q46D33qF8pixerh8Xdu6tfI9VtOWg0BmPXb14p5UopZQcpZVsp5cuWbc9JKZdZ3g+UUjaRUsZZXsMqrtF4unSB++9XhivNmjnmHH0sj7m//77qxwohmDF4BjmFObV+Bm33Zt0Z13UcWflZJJ9NZuG+ha4Ts3q1Gmp18CB89pl6bdmiullA7XvrLTVi5+OP1ZeEq7upWraE995Tv0hGWRLm7d0LY8bA//2fa7VpNBY8YhpnVpaa19O8uXrvqOA/eLBabtpUteN+z/6dzPxMFh1cRIBvgFtMoHpt4GsE+QXRPqw917a+1vkCzGZ46SXVh96ypWrhN2qkZvINHarG81uHbP7tb7Vz8lXbtmqi2Y4dqjsI4Ouv4aOPYP9+ePNN9ZxCo3EBHhH8jxxRS2tfv6OCf/v2KtZU1dTl4dUP0+29bnyx7wsGtx1Mg6AGDtFnJM1Dm/PSX17i6B9H2XjSwKRG9vDnn6oL57nnYOxY1aUzYQIMGqS+5f/7X/VlYO/DWlcTHw/ffadezZur5xIDBsDf/65mG7/1lnIi0miciJv891SMNRhbUzs4Kvj7+Khf8VlZ9j+/23hiIwv3LaRfq36czjrN7Z2dPvG52jzU+yHimsbx1xV/ZcBHAziX7YRn+Xv2qJm5q1erkTwffQSBgXDsGMyerbpSJk1yz9m1AweqB0ZffHHpZi0ogClTVE6h2bNdq0/jVXhE8N+zR8UCazxwZA9AYqKaCHrcjsxFZmnmkdWPEBkaSaG5kHqB9bgl5hbHiTMYPx8/3r3hXc7lnCMpNYl/fv9Px57wk08ujeaJi4PkZNV/36MH/PqrmlEbGOhYDY7Gxwduu021WN5559LnkVLlDrG+L9C2mhrH4hHBf+hQlcLlnKVh6qiWP1yyh/z228rLfrj7Q3ac2cHzic+z9NBSRnceTbC/e2WIvLLFlTx25WOYpZl5u+ax9VQ1hjpVRkGBehB6550qkZr129XqnwnGD9N0Nf7+KldQSop6tpGXp37lTJyofhm0baseGusvAY2D8Ijgn5iohnOfOaNa/2EOHEVZv75aLl9eednvf/2efq36YTabyTXlMqH7BMcJcyDTBkyjU6NO+Agf7v/mfkxmA4crnj4N116rJkWBahm/9ppq6T/yiHHnqa2EhMAzz6jhoY88okYs3X23av1PnqyGsc6fr4eIagzH7YN/ZqYafZOfr+bZREY6dqRf795qac9Y/09u+YRvxnzDvN3z6BTRid6RvR0nzIEE+gXyxa1f4IMPu8/u5v2d7xtT8QcfqDG6yckq+E+dqoL+P/5hjGmCOxERoZ5xHD4Mt9+uHmzXrauC/r33qhaOniimMRC3D/4//qhG/23Zop4JGm2yXppGjVTXz9mzlwzjS7Pu13Wk/pmKEIKD5w+y9dRW7u9xf+1Kk1BFYhvH8tagtwA4lV4mu0fV2LNHzYydMEF9a2/ZorJjPv30pQeh3krr1qr7JzlZTQ5LS1M/ZTt3VllITSY1XNTs4ol3GrfH7YP/5s0qTXuPHs4J/qDOBeX3+5/OPM1t/7uNCV+rLp6ZW2ZSL7Ae98Td43hhDub/+vwf4+PG8/KGl1mwewEFRVXsj05OJvaf/1T5dvbuhY4d1Xh3q1mC5hJdu8KKFZCUpMYYz5mjpphPnqxyGnXvTqMNG/SvAU218Yjg362b+h84d845wX/wYPWFU7rf32Q2cceiO8gpzGH20NmcTD/JV/u/4r74+wgNdP8WrRCC9254j96Rvbnn63u45+sqfqHNm0f45s2qX+7ll9WIlyuucIxYT+Gaa1S/5pIl6ifn3LlQrx6cOkXss8+q/EErV+ovAU2VcevgX1ioegyuvFJ1FYNzgv9jj6nULWvWQG7upe3/Wvcvfjz+I+/e8C4xETHM2DwDieSh3g85XpSTCPQLZNnoZdQNqMtnez/j9Y2vl19QSli3TjleffaZmpg1Zw4FDRsqP91//tP1aRjcBSFUa3/rVvjhB9XPeeECZj8/Ne/h0UcvldVfAho7cevgv2mTmnA1cKDq8gE1V8bR+PqqodpZWWouEsA3h79h2oZp3Bd/H3d1u4tTGaeYvW02d3a9k9YNWjtelBNpEtKEbfdtI8A3gCfWPsGX+7+8tLOwUPnn9uypUhps2aJG7zzwAFxzDdvnzlXOWJqqI4QaGbVqFSQnc+7aa9XM4KNH1bVesEB1F/3732okhEZTAW4d/Pv0UY3LgQPV/Q9qeLQz2L5dpYv/6CO1PqDNAF669iXeHqqGLE79aSpmaeZf1/zLOYKcTMdGHfnhrh/wFb6M/t9ovtr/ldpx440qz31OjnqI6+OjRrDMnAmrVlHYsKFrhXsKXbty6J//VKOCXnlFZT4dP15d68cegxYt4PHH1bpGUw5uHfwDA9UIuNBQlUq9ZUvHjvEvSXi4GnixbJlk39E/CfYP5pn+zxDgG8Dh84eZu2su98XfR1RDJ/wUcREJLa7im44vEmz2ZezisXyc/LFKsrZggXqIO3u2ymWzYwc8/LD75OJxJxo3hiefVJPFli9X+Y+EgIwMmD5dPVTf6OTcTBq3wG3/G48eVQ0b63DLnTtVynRnMWwYCCExmyWjpvxYvN0szdz3zX2EBITw7DXPOk+QM0lNhWnToH17rr/9aU68H0Lf+l24a+ldjN/9Arl/e0gFopdfVv3UejSP4/H1Vc9Xli9Xsx1nzrw0Q7pfP/XguH9/NYN45079bEDjvsH/009V1yaorLiHDzs3+H939mO44kd8AnM5sfam4tQSc3bMYf2J9UwfNJ2mIbUwzXB1sY4rP35cPVh5+mnVtfDpp4T9cpo1rZ+hbaY/C/I3Ez2pgK1Jn6qHuu6YgM3dadJE/dLatUulvP3Xv5S/6fr1ygehRw81YWXiROWToPFK3DL45+erwSODBqlehd27VUPGGcE/z5TH5OWTuWvpXXQcsh5zfl3y83x48knYdWYXj3/7OAPbDGR83HjHi3EkUqrJWNOmwVVXqdTKoIZmvv++Gl6VlKQSKQ0ejP/wEaxY15xIvzBOBuXTZ80o7ll6D6cyajghTFMz2rdXwX/vXjUqYto01RX0xx/qi6BrVzV2eepUNZcg3+n22xoX4ZbBf8ECZdc4ZYpaX7VKdSf37ev4cxeZi/j22Lf84+p/sPWtJ3n4YeUg9sEHMOiZtwmvE87Ht3zs1rN5efllFeS7dVMtfLP50sw2ULlnfv750oieI0fgnXeI3nyUPY8dZWj7oQAsSF5A21lteWzNY/x68VcXfRhNMVFR8NRTqrV/8aIK9nfdpR4aP/usupGDgtSvggEDlM9AaqruIvJUpJQueXXo0EFWhwsXpGzUSMqEBCnNZrWta1cp+/WrVnWVsm7dOrklbYscu2iszCnIkVJK+Wfun5eVOXQmVQa23CcJ/FN+uvqgY4RUwrp166p2QFaWlBs3Svn221LeeaeUHTtKmZ+v9r30kpQjRkg5b56UZ86obWazlLt2Sfn441I2ayYlqGPmzJEyJ+eyqs1ms3x/x/uy9YzW8tYvb5V+L/pJ8byQtyy8Rb741YsytzC3xp/X0VT5eroAQzWmpEj5yCNSduokZWCg+vtaX40aSdmzp5TDhkn55ptS/vST+kd0hU4H4i46ge3SgBjs5+ovn6ri768aJU8+qQY17NuneifeeMPY81htF2fsnkHyj8mE1wln/+/76dm8J/WD6heXW5OyhjvnPY/J/3Ua1W/H/aM6IuYou9ZaQV6eGgly+LAaGhUerlp8kydf6sdv3Fg9DLx4UfUXP/OM2p6fr0aK/Pvf6kHigQPqDzB0qOovHjq03BE8Qggmxk/knrh78PXx5UT6CXrO6cnKoytZUrSEN46+wZB2QxgQNYABbQbQtmFb9/6l5Am0bQszZlxaP3JEPViLiFD9qitWqPHNy5ZdKhMcrBLztWmjchK1bavet2ypugO9LTmfm2FX8BdCDAFmAr7AXCnlq6X2BwIfAT2AC8DtUspUI4WeO6ce7EZFwUKLn7iU6tdqvXoqR1hNyMjPIKsgi+ahzTmdeZq2s9qSZ8ojsk4krw18jck9J1+WomHbqW288fMbfHXgK6ICBlKUcRU+Pn60bg133KH+b558Us1nclhcy8lRue9//53wn39W/3itWqnRHI8/rn6yHz9+6Wf7ypVw/fUqNemzzyp7wfh4lQq1sFD14//4o/on375dTdDKyVEBPyEBHnpIzW4LD7dLnq+PLwARwRHcE3cP7+98n/yifPJN+Sw/spyvDqi5Ac1DmxPXNI64JnF0bdKVdmHtaN2gNWF1wvSXgqvo0AFeeOHS+sWLasju99/Dhg1w6JCa3h4SonKsfP552TqCguCKK+hWp476x42MVA/pwsPVmOyGDdXS+goJ0bO+nYiQlfTnCSF8gSPAdUAasA0YI6U8UKLMg0BXKeUDQojRwC1Sygr9CqOjo+XhEhNQpFQNUV9f9f6331TsOnIE1q6VLFmiHuiuWiWQEs6cLWT6W4L/vufH8y8V8Njf8xFCEOwfjJSSzPxMMgoyyMrPJLsgm+yCLIL8gujevAcSyVsb3+DohRTOZJ3h8IXDHM84wW0xo5h38weYzWZeTnqBga2uJScli+i4DpzL+Z2U9F/ZdWEf3x1by+ELhwkNCOWx7g/yVI9HSDkazM13hpJyzIfWV0jOn5dkZfvQrLGJhLhsYtvmEtvVl8bt6xNWv4iGv2ynrn8BfmbLqzAX3+h2alhkerpKcZydrV45OWo68ejRahLVoUOqpV4ytwTArFmqP37fPjXlv2lTNSKnWTPVom/USM38PH8efv9dfaOmpqpc8idPXvqSCAhQ/f19+qin6tbJFDWkoKiA6Uunc8j3EN//+j2PXfkYQX5BLDuyjB9+/YHCokIkl+7HOn51aNOwDWF1wgj2Dya8TjhN6jYhom4E9QLr0SCoASEBIQT5BRHsH1z88vf1x1f44ufjh5+PH74+6r11W8l1W18uSUlJJCYm1vgzOxKXaywqUv+woEZgfPedeqh8+rRqlISFQf/+pB86RP0jR+wzpgkMVK+gINWqa9ZMfSnk5KhtwcFqWaeOynXUqJEqn5t7+bGBgeqeDQ1Vv05NJjUr09//0svPT718fcHXl83btnFlQkLxevE+669bIS5/ld5mT5mS26qJEGKHlLJntSuw1mNH8L8KeF5KOdiy/hSAlPKVEmXWWMpsEkL4Ab8BEbKCyoXoIWF7yS2AtLzsfQ5thj7/gSGPqsOdQHAB9D8ONx2BcXugXonBEdkE8zYP8Qnj6Mx+rmcVK7iB/zEKaddnunS5fCnCDxMSKKCsdaEfJvwpRCLII6jMfn8K8MeEGUEedcrsDyQPf1FEEb7kyrLHBwUJ/PzU/0xeXlmlQUFUuL9OHfV/U1hYegCJBAR16kjLfkF+fjm3iX82fr6+mE2+mE3luHj5Z4OQUBSgXmX2Z6l7whQA5nL2B2SppSkQzKXrlxCQbdkfBObSP5DNEGAxXC8MAllqvzCDv3V/HZC+pfYXgb/li7swGKRP1fb7mMDPctEL6lLm5vcpBL/8yvdLoDCEMvgUgF+B7f2+BeolBRSW07Xjmw++hUp3YTnOdb554GsCsw+Yytnvl6c+o9kXTGXvXfxywaeogv054GOGIj8oKntvX9rvD0Xl2IIW31uV7XfRvYef04L/KGCIlHKiZf1OoI+U8qESZfZZyqRZ1n+xlDlfqq5JwCTLaiywr6YfwAk0As5XWsr1aJ3G4g463UEjaJ1GEy2lrPFPcac+8JVSzgHmAAghthvx7eVotE5j0TqNwx00gtZpNEKI7ZWXqhx7+iJOAS1LrLewbCu3jKXbpz7qwa9Go9FoaiH2BP9tQHshRJQQIgAYDSwrVWYZcLfl/Sjgh4r6+zUajUbjWirt9pFSmoQQDwFrUEM950sp9wshXkRNNlgGzAM+FkKkAH+gviAqY04NdDsTrdNYtE7jcAeNoHUajSE6K33gq9FoNBrPwy1z+2g0Go2mZujgr9FoNF6Iw4O/EGKIEOKwECJFCPFkOfsDhRBfWPZvEUK0drSmcjS0FEKsE0IcEELsF0I8Uk6ZRCFEuhBit+X1nLN1WnSkCiH2WjSUGfIlFLMs13OPECLeBRqjS1yn3UKIDCHEo6XKuOR6CiHmCyHOWeamWLeFCSG+E0IctSzL9ZoUQtxtKXNUCHF3eWUcqPENIcQhy990iRCigY1jK7w/nKDzeSHEqRJ/16E2jq0wLjhB5xclNKYKIXbbONaZ17PcOOSw+9OI7HC2XqgHxL8AbYAAIBnoVKrMg8B7lvejgS8cqcmGzmZAvOV9KCqdRWmdicByZ2srR2sq0KiC/UOBVag5hlcCW1ys1xc14/uK2nA9gf5APLCvxLbXgSct758EXivnuDDgmGXZ0PK+oRM1DgL8LO9fK0+jPfeHE3Q+D0yx456oMC44Wmep/W8Bz9WC61luHHLU/enoln9vIEVKeUxKWQAsBIaXKjMcWGB5/z9ggBDOze4kpTwjpdxpeZ8JHAQinanBQIYDH0nFZqCBEKKZC/UMAH6RUh53oYZipJQ/oUaklaTkPbgAuLmcQwcD30kp/5BSXgS+A4Y4S6OU8lsppcmyuhk138al2LiW9mBPXDCMinRaYs1tQDmZ6ZxLBXHIIfeno4N/JHCyxHoaZYNqcRnLzZ0O2Jc20gFYup26A1vK2X2VECJZCLFKCOEqY1oJfCuE2CFUuozS2HPNnclobP9j1YbrCdBESnnG8v43oEk5ZWrTdZ2A+nVXHpXdH87gIUv31HwbXRS16Vr2A85KKY/a2O+S61kqDjnk/tQPfEsghAgBFgGPSikzSu3eieq66Ab8B1jqZHlW+kop44Hrgb8KIfq7SEelCDUpcBjwVTm7a8v1vAypfkPX2vHPQoinARPwqY0irr4/3gXaAnHAGVSXSm1mDBW3+p1+PSuKQ0ben44O/m6TGkII4Y+64J9KKReX3i+lzJBSZlnerwT8hRCNnCwTKeUpy/IcsAT1E7ok9lxzZ3E9sFNKebb0jtpyPS2ctXaNWZbnyinj8usqhBgP3AiMtQSBMthxfzgUKeVZKWWRlNIMvG/j/C6/llAcb0YAX9gq4+zraSMOOeT+dHTwd4vUEJZ+v3nAQSnldBtlmlqfRQgheqOunVO/pIQQdYUQodb3qIeApTOjLgPuEoorgfQSPxmdjc1WVW24niUoeQ/eDXxdTpk1wCAhRENLV8YgyzanIJSh0j+AYVLKHBtl7Lk/HEqp50u32Di/PXHBGQwEDklLNuLSOPt6VhCHHHN/OuEJ9lDUU+tfgKct215E3cQAQahugRRgK9DG0ZrK0dgX9VNqD7Db8hoKPAA8YCnzELAfNTJhM3C1C3S2sZw/2aLFej1L6hTAbMv13gv0dLZOi466qGBev8Q2l19P1JfRGaAQ1S96L+oZ0/fAUWAtEGYp2xPlXGc9doLlPk0B7nGyxhRUn671/rSOkGsOrKzo/nCyzo8t990eVNBqVlqnZb1MXHCmTsv2D633Y4myrryetuKQQ+5Pnd5Bo9FovBD9wFej0Wi8EB38NRqNxgvRwV+j0Wi8EB38NRqNxgvRwV+j0Wi8EB38NRqNxgvRwV+j0Wi8kP8HVweDyXgfE6wAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# check that lognormal makese sense, compare with gaussian\n", "color = ('red', 'green', 'blue')\n", "ax21=plt.subplot(111)\n", "for b,e,c in zip(bg, err, color):\n", " x = np.linspace(0.01, 20., 400)\n", " ax21.plot(x, lognorm.pdf(x, np.log(1+e/b), 0, b), color=c, \n", " label='numba logn $\\mu$ = '+str(b) + ' $\\sigma$ = ' + str(e))\n", " ax21.plot(x, norm.pdf(x, scale=e, loc=b), color=c, linestyle='dashed', \n", " label='gaus $\\mu$ = '+str(b) + ' $\\sigma$ = ' + str(e))\n", "\n", "ax21.grid()\n", "ax21.legend()\n", "ax21.set_xlim([0, 20.])\n", "_ = ax21.set_ylim(0,1.5)" ] }, { "cell_type": "code", "execution_count": 5, "id": "cf07671a", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAhCklEQVR4nO3df5AU9Z3/8eebH7qJIquI+0XWzGBg8cepgES5r5paJBcJsYK5oGJZipFvkYtaFXJ3nHy/Fid+62vqPBP1tO6b7zff03KxzCHRi1KW+arR3bJiCo0IGhSNq5mJSxAjwUWS76ro5/vHdC+9s/N7e2a6p1+Pqqnt6e7peW9vb7/786M/bc45REQkucY1OwAREWkuJQIRkYRTIhARSTglAhGRhFMiEBFJuAnNDgCgvb3dzZw5s9lhlPWnP/2JI444otlhlKU4wxOHGEFxhi0ucW7duvU959zUsW4nEomgo6ODF154odlhlNXX10d3d3ezwyhLcYYnDjGC4gxbXOI0s2wY21HVkIhIwikRiIgknBKBiEjCRaKNQETi5+OPP2ZgYIChoaGKPzN58mR27txZx6jCEbU429ra6OzsZOLEiXXZvhKBiNRkYGCASZMmkU6nMbOKPvPBBx8wadKkOkc2dlGK0znH3r17GRgYYMaMGXX5DlUNiUhNhoaGmDJlSsVJQGpjZkyZMqWqkle1lAhEpGZKAo1R7/2sRCAiknBKBCISjnQazEq+Jh11VNl1hl/pdLN/o1AtXryY9vZ2LrzwwhHzf/vb33L22Wczc+ZMLr30Uj766KOGx6ZE4EmnW+64E2msbBacK/n6YP/+susMv7Kh3DQbGWvWrOG+++4bNf/666/nu9/9Lv39/Rx99NHcfffdDY9NicCTzbbccSfS0jKZDCeddBJXXXUVXV1dXH755fz85z/nnHPOYdasWTz//PNAbtygq6++mrPOOou5c+fyyCOPDH/+vPPOY968ecybN49f/vKXQG54iSVLlrBs2TJOOukkLr/8csJ4kuOiRYtG9URyzvH000+zbNkyAFasWMHDDz885u+qlrqPikhs9ff385Of/IR77rmHL3zhC/z4xz/mF7/4BZs3b+Z73/seDz/8MDfffDPnn38+99xzD++//z5nnXUWX/rSlzjuuON48sknaWtr44033uCyyy4bHvPs5Zdf5pVXXuH444/nnHPO4dlnn+Xcc88d8d233nor999//6iYvvjFL3LnnXdWFP/evXtpb29nwoTcqbizs5Ndu3aNca9UT4lARGJrxowZnHbaaQCceuqpLFq0CDPjtNNOI5PJAPDEE0+wefNmvv/97wO5bq+/+93vOP7447nuuuvYvn0748eP5ze/+c3wds8880w6OzsBmDNnDplMZlQiWLNmDWvWrGnAb1l/SgQiEluHH3748PS4ceOG348bN46DBw8CueqXhx56iNmzZ4/47Pr16+no6OCll17i008/pa2tbXjZYYcdNjw9fvz44W0FhVEimDJlCu+//z4HDx5kwoQJDAwMMH369Io+Gya1EYhIS7vgggu46667huv5t23bBsDg4CDTpk1j3Lhx3HfffXzyySdVbXfNmjVs37591KvSJAC5+wMWLlzIgw8+CEBPTw9Lly6tKo4wKBGISDhSqXC7j6ZSoYS1bt06Pv74Y04//XROPfVU1q1bB8A111xDT08PZ5xxBq+99lrdH0Rz3nnncfHFF/PUU0/R2dnJ448/DsAtt9zCbbfdxsyZM9m7dy8rV66saxyFWBit4WM1e/Zs9/rrrzc1Bv/GvVK7Iy4Pq1Cc4YlDjNCcOHfu3MnJJ59c1WeiNIZPKVGMs9D+NrOtzrn5Y922SgQiIgmnRCAiknCJ7zWku4lFJOkSnwh0N7GIJJ2qhkREEi7RJYJ0emQPNZUORCSJEp0I/MESfXrGhkjt0nekyQ6GdzWVmpwiszoT2vaabfz48cPDYXzuc59j8+bNTY7okEQnAhEJT3Ywi7ux9H1J1fTPt5ta68rsM5/5DNu3b292GAWpjUBEYiluw1BHmUoEIhJbcRqGemhoiPnz5zNhwgTWrl3LRRddVJd9UgslAhGJrTgNQ53NZpk+fTpvvfUW559/Pqeddhqf//znx/Lrh0aJQERiK07DUPvDS5944ol0d3ezbdu2yCQCtRGISEuLwjDU+/bt48MPPwTgvffe49lnn+WUU04Z428WnopKBGaWAT4APgEOOufmm9kxwANAGsgAlzjn9pmZAf8CLAH+DFzlnHsx/NBFJEpSk1Oh9vRJTQ5vGOrVq1dz+umn8+mnnzJjxgweffRRrrnmGr7xjW+wYcMGFi9eXNdhqHfu3Mm3vvUtxo0bx6effsratWsjlQhwzpV9kTvRH5s375+Btd70WuAWb3oJ8DPAgAXAc+W239XV5ZoBSr/P19vbW7dYwqQ4wxOHGJ1rTpyvvvpq1Z/Zv39/HSIJXxTjLLS/gRdcBefwcq+xVA0tBXq86R7gosD8DV6cW4B2M5s2hu9pqHRaA9GJSLJUmggc8ISZbTWzVd68Dufcbm/6HaDDm54OvB347IA3LxayWQ01ISLJUmmvoXOdc7vM7DjgSTN7LbjQOefMrKo7LryEsgpg6tSp9PX1VfPxMVu+fAEdHdDXtyUwt3t4qlA8Bw4caHictVCc4YlDjNCcOCdPnswHH3xQ1Wc++eSTqj/TDFGMc2hoqH5/42rrkoD1wN8DrwPTvHnTgNe96f8NXBZYf3i9Yq9mtBEUag/IjTxUvK1A9cXhikOccYjRObURhC2KcTa1jcDMjjCzSf408GVgB7AZWOGttgJ4xJveDFxpOQuAQXeoCklERCKmkqqhDuCnuV6hTAB+7Jz7v2b2K2CTma0EssAl3vqPkes51E+u++g3Q49aRERCUzYROOfeAs4oMH8vsKjAfAdcG0p0IhIb6XQlHS0qG3kUcs8K8UaJiL3t27fz7W9/m/379zN+/HhuuOEGLr300maHNUxDTIhIKPKf71FIVcNQt9Ao1J/97GfZsGEDs2bN4ve//z1nnnkmF1xwAe3t7c0ODdAQEyISU3Eahrqrq4tZs2YBcPzxx3Pcccfxhz/8YUzbDJNKBCISW3Eahtr3/PPP89FHH0VmwDlQIhCRGIvTMNQAu3fv5oorrqCnp4dx46JTIaNEEOA/yF53FovEQ5yGod6/fz9f/epXufnmm1mwYEGVv2l9KREE+D0UWqmRSiTp/GGo77rrLsyMbdu2MXfuXAYHB+ns7GTcuHH09PTUNAx1pSWCjz76iK9//etceeWVLFu2rJZfo66iUzYRkVhLpXIXUaVeRx01qew6/isVzijUrFu3jo8//pjTTz+dU089lXXr1gFwzTXX0NPTwxlnnMFrr71W12GoN23axDPPPMO9997LnDlzmDNnTqQeZG9jbQ0Pw+zZs93rr7/e0O80K97Vrdiyvr4+uru76xpXGBRneOIQIzQnzp07d3LyySdX9Zlquo82UxTjLLS/zWyrc27+WLetEoGISMIpEYiIJJwSgYjULApVy0lQ7/2sRCAiNWlra2Pv3r1KBnXmnGPv3r0jureGTd1HRaQmnZ2dDAwMVDVUwtDQUF1PaGGJWpxtbW3DN7jVgxKBiNRk4sSJzJgxo6rP9PX1MXfu3DpFFJ64xBkWVQ2JiCScEkEBqVRubHURkSRQ1VABmYyGmRCR5EhkiSCdDu/2dRGRuEtkiaCSJymJiCRFIksEIiJyiBKBiEjCKRGIiCScEoGISMIpEYiIJJwSgYhIwikRiIgknBKBiEjCKRGIiCRcxYnAzMab2TYze9R7P8PMnjOzfjN7wMwO8+Yf7r3v95an6xS7iIiEoJoSwXeAnYH3twC3O+dmAvuAld78lcA+b/7t3noiIhJRFSUCM+sEvgr8m/fegPOBB71VeoCLvOml3nu85Yu89UVEJIIqHXTuDuAfgEne+ynA+865g977AWC6Nz0deBvAOXfQzAa99d8LbtDMVgGrAKZOnUpfX19tv0FNuiv4vtHrHDhwoMFx1kZxhicOMYLiDFtc4gyNc67kC7gQ+J/edDfwKHAs0B9Y5wRghze9A+gMLHsTOLbUd3R1dblGgtrW6e3tDT2WelCc4YlDjM4pzrDFJU7gBVfmHF7Jq5ISwTnA18xsCdAGHAX8C9BuZhNcrlTQCezy1t/lJYYBM5sATAb2ji1diYhIvZRtI3DO/VfnXKdzLg0sB552zl0O9ALLvNVWAI9405u993jLn/Yyl4iIRNBY7iO4HvhbM+sn1wZwtzf/bmCKN/9vgbVjC1FEROqpqieUOef6gD5v+i3grALrDAEXhxCbiIg0gO4sFhFJOCUCEZGEUyIQEUk4JYIiUilIp5sdhYhI/VXVWJwkmQxoYAwRSQKVCEREEk6JQEQk4ZQIREQSTolARCThlAjKSKfVe0hEWpsSQQmpVO5nNtvcOERE6kndR0vIZHI/1Y1URFpZ4koE6fShK30REUlgiSCbBT0dQUTkkMSVCEREZCQlAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUk4JQIRkYRTIhARSbiyicDM2szseTN7ycxeMbObvPkzzOw5M+s3swfM7DBv/uHe+35vebrOv4OIiIxBJSWCD4HznXNnAHOAxWa2ALgFuN05NxPYB6z01l8J7PPm3+6tJyIiEVU2EbicA97bid7LAecDD3rze4CLvOml3nu85YvM9IwvEZGoMlfBU1rMbDywFZgJ/CtwK7DFu+rHzE4Afuac+wsz2wEsds4NeMveBM52zr2Xt81VwCqAqVOnnrlp06bwfqsSFi7spre3r6bPHDhwgCOPPDL8oEKmOMMThxhBcYYtLnEuXLhwq3Nu/pg35Jyr+AW0A73AuUB/YP4JwA5vegfQGVj2JnBsqe12dXW5Rkilcq9q+eWi3t7eEKOpH8UZnjjE6JziDFtc4gRecFWcw4u9qnpUpXPufTPrBf4SaDezCc65g0AnsMtbbZeXGAbMbAIwGdg7hlwVGj2mUkRktEp6DU01s3Zv+jPAXwE7yZUMlnmrrQAe8aY3e+/xlj/tZS4REYmgSkoE04Aer51gHLDJOfeomb0KbDSz/wFsA+721r8buM/M+oE/AsvrELeIiISkbCJwzr0MzC0w/y3grALzh4CLQ4lORETqTncWi4gknBKBiEjCKRGIiCScEkEFUilIp5sdhYhIfVR1H0FSZTKgQTJEpFWpRCAiknBKBCIiCadEICKScEoEIiIJp0QgIpJwSgQiIgmnRCAiknBKBCIiCadEICKScIlJBOl0bqgIEREZKTFDTLTqYyrTd6QByKzONDUOEYmvxCSCVpUdzDY7BBGJucRUDYmISGFKBC0sfUd6uOpIRKQYVQ21MFUbiUgllAhiSlf6IhIWJYKY0tW+iIRFbQQxlL4jTWqybooQkXAoEcRQdjCr+wZEJDSqGmpBy7csp217W7PDEJGYUCKIkVINxHaT4W7M3Tq958M98GGDghKR2FPVUIxkB7MjGolTk1Mj2grUk0hEaqESQcykJqeGk4HfTuAngOxglvQdaToO78iVCjwaj0hESilbIjCzE8ys18xeNbNXzOw73vxjzOxJM3vD+3m0N9/M7E4z6zezl81sXr1/iSTwewoVOplnVmeG52cHs2xcsHFEaSG/JCEiElRJ1dBB4O+cc6cAC4BrzewUYC3wlHNuFvCU9x7gK8As77UK+GHoUVepFYagrranUDA5iIiUUjYROOd2O+de9KY/AHYC04GlQI+3Wg9wkTe9FNjgcrYA7WY2LezAq5HNQibTzAhERKKrqjYCM0sDc4HngA7n3G5v0TtAhzc9HXg78LEBb97uwDzMbBW5EgNTp06lr6+vytCr0R3C9rs5cOBAneMsLfjdpeIIxtlxeO7PsufDPU2NvZBm789KxCFGUJxhi0ucoXHOVfQCjgS2An/tvX8/b/k+7+ejwLmB+U8B80ttu6ury9UTjH0bqZRzHR3/b+wbqhHry/8SrMexHtfb21vT5xutUJxRE4cYnVOcYYtLnMALrsJzeKlXRd1HzWwi8BBwv3PuP7zZe/wqH+/nu978XcAJgY93evNiLZOBPXta5yat9B1p7CZTl1MRqajXkAF3Azudc7cFFm0GVnjTK4BHAvOv9HoPLQAG3aEqJIkAu8nIDmZxNzr1JhKRitoIzgGuAH5tZtu9ef8N+Cdgk5mtBLLAJd6yx4AlQD/wZ+CbYQYsxZUbiE5X/yJSSNlE4Jz7BWBFFi8qsL4Drh1jXJHU0TFEOt0W2R5IfnfRQo1cwXsKRESCNMREFTZu3EI2pudR3VcgIsUoEYiIJJwSgYhIwikRxICeSCYi9aREUKVUKjd2USOF+UQyJRQRyadEUKVMhtg2GIOGohaR0ZQIREQSTg+miTDdACYijaASQUQFnzomIlJPSgQRFWYDsYhIKaoairh69PJRzyERCVIiiLh6lApU0hCRIFUNiYgknBKBiEjCKRGIiCScEkEENXNsofQdad2/IJIwaiyugT/eUL0eUOM/RrIZdN+CSPIoEdQgkwEr9sw2EZGYUdWQiEjCKREkXGpyCrvJ1C4gkmCqGko4/+Yyu0l1XSJJpRKBALmSgUoF5HoBNPrJQyJN1vIlgnQ618snDtJ3pMkOZpvSdTSzOqNSAcT7qUMiNWr5RJDNgmtOT8yqNbPbaD6/dJCocYn8qwYlA0mYlk8EUpuWuJ/Ar+Lxb/gwO3RVkL8MDl01qG+wJIwSgYziV03FPhmUurKv9aq/UAIRiTk1FkdEM4eVyJdZnYlHlVAlDbv+beDlPlOoMSmdpnvhwpHr+wlEDcrSQsomAjO7x8zeNbMdgXnHmNmTZvaG9/Nob76Z2Z1m1m9mL5vZvHoG30r0RLIaZLPlr+wzmdw6/onbXz+Vyr3Mcsuy2dFX+dksfb29o7/D32YxqlqSmKmkRHAvsDhv3lrgKefcLOAp7z3AV4BZ3msV8MNwwpRGSE1ORaZUUpN0+tCJPV82e+jkn8kcejlX+KRebXczdTuVGCvbRuCce8bM0nmzlwLd3nQP0Adc783f4JxzwBYzazezac653aFFLHUT+xJJscZe/4ReTb2+v62+vsrXF4mpWhuLOwIn93eADm96OvB2YL0Bb96oRGBmq8iVGpg6dSp9lf7DVa07tG0fOHAgsK3wtrt8y3I6Du+oU5xjU7+/y9jj7PZ++tvo9qbzf3Lvvf6KFW1nQUcHbXv20NfXx4EDB4aXdS9cODy9oKODNjOGOjrYsnHjiG10U9/9VkiYf/N6UpwR5Zwr+wLSwI7A+/fzlu/zfj4KnBuY/xQwv9z2u7q6XL1AeNvq7e0dnk6lcq8wsD7EIN3IOMeC9bjU7alQtlXImOL0/wC56/aRfxD/j17pH7/QH9P7bG9v78jt5a+XH0f+94Z5oJQQ1t+83hRnuIAXXAXn8HKvWksEe/wqHzObBrzrzd8FnBBYr9Ob13KSMhR1ZLuQ5lcD5d85WE0dfzVVRvnrBg+E/B5KflwiEVdr99HNwApvegXwSGD+lV7voQXAoGvh9oFC//fVilK30XzBu5xj9+SyQr2AalXpHzr4faW6qapRWSKmbInAzP6dXDXqsWY2ANwI/BOwycxWAlngEm/1x4AlQD/wZ+CbdYg5MsIoFURpWIlSmloyqOQmrnoOKJXJlC5hBOf7B0Wwl1KwVKASgkRQJb2GLiuyaFGBdR1w7ViDEhkhePIsdjUd9p2++QNUldp+/rL8Xkr+1UIS6hIllnRnscRL8IawYpo93Kx/j0KQn8DyY1NykAjQWENSlt+O0dTqoWA9fbmr/6i15Psnf7+UUS624OB4Ig2gRNBEUW4oDvLbMZr6vIJiJ/diV/9ROpFqgDqJOCWCJopLQ3GkxfEk65du4hi7tCS1EYg0mp8AinUvFWkwlQiaJC7VQkH+c40bNiZRK58U86u68kdHFWmglk4EUX5ecVyqhYLJqqHPNU7CiTF4j0GxEVBBVUhSdy2dCKL6vOI4lQYaPiJpMAE4N7JUENWsXqtyvZtaOQlKpKiNoAni/BAav3qobvIfNhPsk1+of36rC2McE5EylAjGKGn/p5nVmegORBdH/lPSipV2Cj0NTeMVSchaumqoEaJ271JsJfXEFizh5DceB5NDsL1AVUYSMpUIpCahj0aa/yzhVmsPqETw984fPdXfP6VKDyI1UolAquY3dGcHs8PJILM6M2K6Jkmr/8+XP0hdseX+OhqKQkKiRBCSJPX080/0/onfbzNQ20FIqrni1x3KEoKWrRpq5D0EwdJ8wVji9lCXCmVWZ8Lp/RTlGz6aodreUcH2FTUkSw1atkTQyHsIypXos4PZEd0u43IPQS2quft4wfLlsGdP7o2qOIorlSTzeyuoIVlq0LKJIGqCd+XG4Y7iWpW9+zhQh9a2Z8/I5w5LYYVKB3nJYcHy5dDWNnq9JNVZSs2UCBqoVUsCVd1kVuiKVaWB6gVP7KkUDA0dmvb3cRKG6ZBQKBE0UFzvJi6nUCmgbA+idJqhjg4KXMNKtTIZtvT10d3dnXvvl7CCw3QEG5VVSpA8SgQSivzSzqgeRAXGDNqycSPd9Q0ruYIN8IXaEfJviVdSSDQlAglFqSt/4NDJJ5U6dNLp66t/YElVrreEbomXACWCkOWXuuM00mg9pJcNkPnBJ7k3uupsjEK9jPwxjYotl0RTIghZfrtcXJ47ULNRmS9Natl4ALKTPoHOTtJ/N0Dmwc6mhJdIhRJuqZ5HakxOPCUCGZtRmS9L5gfetFc1YTeZSgNRFLwBplgDshqWE6ElE0GzblTNHywycf87/o4vcIWZviM9fGOdb2hoiLbtbS3bmyo2ggduoeqj/O6oiTuwW19LJoJmPZks2C6Q/W4mOe0Dfg+UQk8V41CPInejG32z2YeNCVFKKHRi9xNC8Ia//ARfaNA7DYQXSy2ZCJot2HUy9le7lVwF+j1Qgt0Vg4sD+yA4cqnPbrLWbkeJo/widf6zEvJpfKNYa7lE0Mzxy+wmIzU5lTvZpfBKBTEvSRe6S7XYFV8Fv2hw5NKhoSH2fLhn+H3sk2YrKfa3LPbPVXTExXTp7Ukk1GX0UTNbbGavm1m/ma2tx3cUk/88j0bzn0ecyeTOl9lsiwwIGfxHT6cPjYdfY9bNrM6wccHG4cSZfwNayfGKpDn8e0D8l98l1T8O/GPBn+dXFwbmdS9cOLJxWiIh9BKBmY0H/hX4K2AA+JWZbXbOvRr2d0WF3xAKo++wLfS0QYjABVKBbp/DN31lMiPf+//M/nwYXSqosV44WEIoNExFsIE5mCxSk1MqQTRKkSq/sgexXzwPrNfnD4XhX0zkX0gEn1Tnvw+uo2Ey6qIeVUNnAf3OubcAzGwjsBSoeyJoVLVQ8MQPuZNSsTru4HGaP2x8uXXyl1VqxEiUhYrs/sk9/yYjv6E3v76/ARks/6Tuj1UUbGD2p92NrmDigJFtEIWSRf4T1YqtV8iYn8AWV7X+zcu1K0HhkkGww0F+54P8O9WDx2/+sZr/PflJJbid4GczmUP/Q6W60+ZfPBVbr5gIJTNzIbfwm9kyYLFz7r94768AznbOXZe33ipglff2L4AdoQZSH8cC7zU7iAoozvDEIUZQnGGLS5yznXOTxrqRpjUWO+d+BPwIwMxecM7Nb1YslVKc4YpDnHGIERRn2OIUZxjbqUdj8S7ghMD7Tm+eiIhEUD0Swa+AWWY2w8wOA5YDm+vwPSIiEoLQq4accwfN7DrgcWA8cI9z7pUyH/tR2HHUieIMVxzijEOMoDjDlqg4Q28sFhGReKnLDWUiIhIfSgQiIgnX0ERQbugJMzvczB7wlj9nZulGxufFcIKZ9ZrZq2b2ipl9p8A63WY2aGbbvdc/NjpOL46Mmf3ai2FUNzLLudPbny+b2bwGxzc7sI+2m9l+M1udt07T9qWZ3WNm75rZjsC8Y8zsSTN7w/t5dJHPrvDWecPMVjQ4xlvN7DXvb/pTM2sv8tmSx0cD4lxvZrsCf9slRT7bsCFpisT5QCDGjJltL/LZRu7Pguehuh2fzrmGvMg1HL8JnAgcBrwEnJK3zjXA//KmlwMPNCq+QAzTgHne9CTgNwXi7AYebXRsBWLNAMeWWL4E+BlgwALguSbGOh54B0hFZV8CXwTmATsC8/4ZWOtNrwVuKfC5Y4C3vJ9He9NHNzDGLwMTvOlbCsVYyfHRgDjXA39fwXFR8rxQ7zjzlv8A+McI7M+C56F6HZ+NLBEMDz3hnPsI8IeeCFoK9HjTDwKLzBr7hG3n3G7n3Ive9AfATmB6I2MI0VJgg8vZArSb2bQmxbIIeNM5F5nnIjrnngH+mDc7eAz2ABcV+OgFwJPOuT865/YBTwKLGxWjc+4J59xB7+0WcvfqNFWRfVmJSs4LoSkVp3euuQT493p9f6VKnIfqcnw2MhFMB94OvB9g9Al2eB3vQB8EpjQkugK8qqm5wHMFFv+lmb1kZj8zs1MbG9kwBzxhZlstN2RHvkr2eaMsp/g/WBT2pa/DObfbm34H6CiwTpT269XkSn2FlDs+GuE6rwrrniLVGFHal+cBe5xzbxRZ3pT9mXceqsvxqcbiIszsSOAhYLVzbn/e4hfJVXGcAdwFPNzg8HznOufmAV8BrjWzLzYpjpIsd2Ph14CfFFgclX05isuVsyPbv9rMbgAOAvcXWaXZx8cPgc8Dc4Dd5KpdouwySpcGGr4/S52Hwjw+G5kIKhl6YngdM5sATAb2NiS6ADObSG7n3++c+4/85c65/c65A970Y8BEMzu2wWHinNvl/XwX+Cm5YnZQVIb7+ArwonNuT/6CqOzLgD1+9Zn3890C6zR9v5rZVcCFwOXeCWGUCo6PunLO7XHOfeKc+xT4P0W+v+n7EobPN38NPFBsnUbvzyLnobocn41MBJUMPbEZ8Fu4lwFPFzvI68WrJ7wb2Omcu63IOv/Jb7sws7PI7ceGJiwzO8LMJvnT5BoQ80dw3QxcaTkLgMFAsbKRil5pRWFf5gkegyuARwqs8zjwZTM72qvu+LI3ryHMbDHwD8DXnHN/LrJOJcdHXeW1R329yPdHZUiaLwGvOecGCi1s9P4scR6qz/HZiBbwQGv2EnKt328CN3jz/ju5AxqgjVz1QT/wPHBiI+PzYjiXXHHrZWC791oC/A3wN9461wGvkOvhsAX4z02I80Tv+1/yYvH3ZzBOI/eQoDeBXwPzmxDnEeRO7JMD8yKxL8klp93Ax+TqUVeSa5N6CngD+DlwjLfufODfAp+92jtO+4FvNjjGfnJ1wP7x6fe0Ox54rNTx0eA47/OOu5fJncCm5cfpvR91XmhknN78e/1jMrBuM/dnsfNQXY5PDTEhIpJwaiwWEUk4JQIRkYRTIhARSTglAhGRhFMiEBFJOCUCEZGEUyIQEUm4/w/L4jcS07XSlgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# make sure that we know how to pick according to lognormal\n", "np.random.seed(12345)\n", "picks = []\n", "nPicks = 10000\n", "for b,e in zip(bg, err):\n", " temp = np.random.lognormal(np.log(b), np.log(1+e/b), nPicks)\n", " picks.append(temp)\n", " \n", "for b,e,y,c in zip(bg, err, picks, color):\n", " x = np.linspace(0.0, 20., 400)\n", " ax22=plt.subplot(111)\n", " _ = ax22.hist(y, x, histtype='step', color=c, label=\"mean = \"+str(b))\n", "\n", "ax22.grid()\n", "ax22.legend()\n", "_ = ax22.set_xlim([0, 20.])" ] }, { "cell_type": "code", "execution_count": 6, "id": "05506bc2", "metadata": {}, "outputs": [], "source": [ "# calculate likelihood function for given mu,b,s,counts. (ignore the constant factorials)\n", "def getlik(mu, counts, s, b):\n", " totalExpected = mu*s+b\n", " temp = 1.\n", " for n,t in zip(counts, totalExpected):\n", " temp = temp * (t**n) \n", " return temp*math.exp(-totalExpected.sum())" ] }, { "cell_type": "code", "execution_count": 7, "id": "bccb814c", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0\n", "0.2\n", "0.4\n", "0.6\n", "0.8\n", "1.0\n", "1.2\n", "1.4000000000000001\n", "1.6\n", "1.8\n", "2.0\n", "2.2\n", "2.4\n", "2.6\n", "2.8000000000000003\n", "3.0\n", "3.2\n", "3.4\n", "3.6\n", "3.8000000000000003\n", "4.0\n", "4.2\n", "4.4\n", "4.6000000000000005\n", "4.8\n", "5.0\n" ] } ], "source": [ "# scan mu, marginalize at a given mu\n", "b0 = np.array(picks[0]) # ugly but works\n", "b1 = np.array(picks[1])\n", "b2 = np.array(picks[2])\n", "nMu = 501\n", "muScan = np.linspace(0,5.,nMu)\n", "lik = np.empty(nMu)\n", "for i in range(nMu):\n", " m = muScan[i]\n", " if i%20 == 0:\n", " print(m)\n", " sumLik = 0.\n", " for bb0,bb1,bb2 in zip(b0, b1, b2):\n", " sumLik = sumLik + getlik(m, data, sig, np.array([bb0, bb1, bb2]))\n", " # print(sumLik)\n", " lik[i] = sumLik/len(b0)\n", " " ] }, { "cell_type": "code", "execution_count": 8, "id": "e0618187", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAasAAAEPCAYAAAD8nOuVAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1HklEQVR4nO3deXgV9dn/8fcnO4QlLAk7AoJhc0NFUKzgitaK1bZarftT+1RrrX1s1fb6VavdrG1trdWWqlVb614rolVRiNYqm+IGBFkVKBAQZJUtuX9/zMQcQ5ZJcpI5J7lf13WuM/Od7T5zaW6+M99FZoZzzjmXyjLiDsA555yrjycr55xzKc+TlXPOuZTnyco551zK82TlnHMu5WXFHUBrUVBQYIMHD447jJSwfft28vPz4w4jJfi9qOL3oorfiypvvPHGBjMrrG8/T1ZJ0qNHD+bOnRt3GCmhpKSE8ePHxx1GSvB7UcXvRRW/F1UkfRBlP38M6JxzLuWlZbKSdK+kMknvJZTdKqlU0juSnpRUkLDteklLJC2SdHJC+cSwbImk6xLKB0qaFZY/IimnxX6cc865faRlsgLuAyZWK5sGjDSzg4D3gesBJA0HzgFGhMfcKSlTUibwB+AUYDjw1XBfgFuA28xsMLAJuLR5f45zzrm6pGWyMrNXgI3Vyl4ws73h6kygb7g8CXjYzHaZ2XJgCTA6/Cwxs2Vmtht4GJgkScBxwOPh8fcDZzTn73HOOVe31trA4hLgkXC5D0HyqrQqLANYWa38SKAb8HFC4kvc/zMkXQZcBlBYWEhJSUkyYk9727Zt83sR8ntRxe9FFb8XDdfqkpWkHwJ7gQeb+1pmNhmYDFBcXGzeuifgLZ2q+L2o4veiit+LhmtVyUrSRcBpwPFWNZz8aqBfwm59wzJqKf8IKJCUFdauEvd3zjkXg7R8Z1UTSROB7wOnm9mOhE1TgHMk5UoaCAwBZgNzgCFhy78cgkYYU8IkNwP4Unj8hcBT9V2/3Gdacc65ZpOWyUrSQ8DrQLGkVZIuBe4AOgLTJL0l6Y8AZjYfeBRYADwHXGFm5WGt6VvA88BC4NFwX4Brge9KWkLwDuue+mJat70iqb/ROedclbR8DGhmX62huNaEYmY/BX5aQ/mzwLM1lC8jaC0Y2e4KWLVpB327tG/IYc455yJIy5pVqnppYVncITjnXKvkySpJsjPgxYXr4g7DOedaJU9WSdIuS8xc9hFbd+6JOxTnnGt1PFklSftssafc+PfiDXGH4pxzrY4nqyTJy4SC9tn+KNA555qBJ6skmlBcxIzSMsorvNOVc84lkyerJDp+WBGbduzhzQ83xR2Kc861Kp6skuhzBxSSnSl/FOicc0nmySqJOuVlc+TAbry4wJOVc84lkyerJDt+WBFL129nxYbtcYfinHOthierJDthWA/AOwg751wyebJKsn5d21Pco6MnK+ecS6Jak5WkCknlUT8tGXSqO35YEXNXbGKLj2bhnHNJUdeo698GKjsMZQP/B2wjmNupDOgBTALygV83Y4xpZ8LQIu4sWcqrizdw6oG94g7HOefSXq3JyszuqFyW9BtgFvDlhBl4kXQd8BgwsDmDTDeH9iugc7tsZpSWebJyzrkkiPrO6gLgz4mJCiBc/zPwtWQHls6yMjP43AGFzFi0ngofzcI555osarLKBIbVsm1EA87TZkwoLmTDtl3M/++WuENxzrm0F3Wm4AeBn0nKAqYQvLMqInhndRMRpn1va449oBAJppeWcWDfznGH45xzaS1qjei7wJ8IEtNC4KPw+8dh+XebJbo01q1DLgf3LWDGIp892DnnmipSsjKz3WZ2NdAXOB44FzgO6Gtm3zGz3c0YY9qaUFzE26s+5qNtu+IOxTnn0lqD3jWZ2UYzKzGzR8zsZTPb2FyBtQbHDS3CDF5+f33coTjnXFqLnKwkDZJ0l6R3Ja0Ov++UNKg5A0xnI3p3onuHXGYs8mTlnHNNEamBhaTDgBnATmAqsI6gU/BZwHmSJpjZm80WZZrKyBDjiwt5Yf5a9pZXkJXpjSadc64xov71/BUwDxhgZpeY2fVmdglBZ+B54XZXg+OGFrFl517mrfw47lCccy5tRU1Wo4FfmtmOxMJw/VfAkckOrLUYN6Q7mRliRqm3CnTOucaKmqw+AbrVsq0rweNBV4NOedkcvl8Xpnuycs65RouarJ4BfiFpXGJhuP5z4OlkB1YXSfdKKpP0XkJZV0nTJC0Ov7uE5ZJ0u6Qlkt6RNCrhmAvD/RdLujCh/LCwAcmS8Fg1Jd7jhhZRunYrazZ/0pTTOOdcm9WQTsHLgJclrZH0tqQ1wMvAcoIR2VvSfcDEamXXAS+Z2RDgpXAd4BRgSPi5DLgLguQG3EDwCHM0cENlggv3+XrCcdWv1SAThhYBMKPUWwU651xjRO0U/JGZjQM+D9wJ/Cf8PsXMxpnZR80YY03xvAJU7+M1Cbg/XL4fOCOh/AELzAQKJPUCTgamhX3HNgHTgInhtk5mNjMcqPeBhHM1ypCiDvQpaOejWTjnXCNFHRsQADN7DniumWJpqh5mtiZcXkvQtB6gD7AyYb9VYVld5atqKN+HpMsIamsUFhZSUlJSa3DFHffwyqJ1TJs+g+yMJj1VTHnbtm2r8160JX4vqvi9qOL3ouEiJ6twENuzgHEEjSo2Av8G/mFme5snvMYxM5PU7HNzmNlkYDJAcXGxjR8/vtZ9y3usY/r9c8nrN5JjhhQ2d2ixKikpoa570Zb4vaji96KK34uGi/QYUFIRMBd4iOBR4KDw+2FgjqRU+Ou7LnyER/hd+cxtNdAvYb++YVld5X1rKG+Ssft3Iycrw99bOedcI0RtYPEbgqbrY8xskJmNNbNBBI0TuoXb4zYFqGzRdyHwVEL5BWGrwDHA5vBx4fPASZK6hA0rTgKeD7dtkTQmbAV4QcK5Gq19ThZjB3WjxN9bOedcg0VNVqcC15rZ7MRCM5sDXE9Qy2oxkh4CXgeKJa2SdCnwC+BESYuBE8J1gGcJWjIuIZjV+PIw9o3AzcCc8HNTwsC8lwN3h8csBf6VjLgnFBeybMN2VmzYnozTOedcmxH1nVUusLWWbVuBnOSEE42ZfbWWTcfXsK8BV9RynnuBe2sonwuMbEqMNTluaA9ufHoB00vLuGTcwGSf3jnnWq2oNauZwLWS8hMLw/Vrw+2uHv27tWdQYb43YXfOuQaKWrP6P4JR11dKeoFg1PUigr5KAsY3S3St0ITiIv76+gds37WX/NwG9Rxwzrk2K2qn4LcIRnKYDBQCJxIkqz8CQ8zs7eYKsLU5bmgRu8sreG1pi/ajds65tBb5n/ZmtoGqIYwaRdIwgsYPo4GeQB5Bf633gVeBF8ysVQ+gd8SArnTIzWJ6aRknDu9R/wHOOecaNq19Y4RNxi+QNAeYD9xI0HfpY+ADIJNg/L4ngLWSJktqta0PcrIyGDe4OyWLygjafjjnnKtP1JmCs4GrgDMJEk1e9X3MrKiWwxeG338Fzjez0lqu0Z7gHdiXgXcl/a+Z/S1KfOlmwtBCnpu/ltK1WxnWq1Pc4TjnXMqL+hjwNuAbBFPazwB2N+AaPyQYkqnOakQ4keOTwJOS+vLZUSRalQnFQV6fXlrmyco55yKImqy+DFxnZr9u6AXM7IlGHLOKzw4m26oUdcpjZJ9OzCgt44oJg+MOxznnUl7Ud1YC3mnqxSRlScqtVnaSpO8kTorYFkwoLuLNDzexaXtDKqnOOdc2RU1WfwZqGzWiIR4hnPwQQNK3CaYc+TkwU9JpSbhGWpgwtIgKg1cW+8C2zjlXn1ofA0q6PGF1LXCepBkEkxR+XG13M7O7qN8YgoYalb4H/NrMvifpToL3W1OjBJ7uDu5bQNf8HGaUljHpkBqny3LOOReq653VHTWU9QeOraHcSKgx1aEbQeJD0oFAb4KOxQCPAedFOEerkJkhxh9QyIxFZZRXGJmtfEJG55xrilofA5pZRgM+mRGvtw4YEC5PBD4ws6XhejugotG/JA2NH1rEph17eGvlx3GH4pxzKa3ZOwVX8xhwi6RbCQbAfSBh26HA4haOJ1bHDikkM0PMKPWBbZ1zri51vbMaDiw1s13hcp3MbEGE610HbAGOIHhs+LOEbYcBj0Y4R6vRuX02h/XvwvTSMq45uTjucJxzLmXV9c7qPYIGEbPD5do69SrcVu+jQDPbC9xUy7Yz6zu+NZowtIhbnitl7ead9Oy8z8AgzjnnqDtZTQAWJCw3maRyYGz1GYfDbYcBsxvw/qtVmDC0kFueK6VkURnnjO4fdzjOOZeSak1WZvZyTctNVFeTt2xgb5KukzaKe3Skd+c8ppd6snLOudo0++x/kvpT1QIQ4FBJ1Z935QEXAsubO55UI4kJQ4t4ct5qdu0tJzerTVUsnXMukroaWKyn9vdU+6hj1PWLgRvCc9XVH+sT4H+iXq81mVBcxIOzPmTO8k2MG9I97nCccy7l1FWz+gMNSFZ1uBN4nKrxBc9j33EGdwMfmtmuJFwv7Rw1uBs5WRlMLy3zZOWcczWo653Vjcm4gJmtB9YDhJMqrjEzH701QfucLMYO6saMRWX86Av19hJwzrk2p0GdgiV1kXSMpHMldQnL8iTVep5wUsVK64EsSe1r+zTuZ6S/44YWsXzDdpZv2B53KM45l3IiJatwao9fEswx9TLBrL+VU88/QfBOqjZbJY0Ol7cBW+v5tEmVEzL6aBbOObevqK0Bfwp8HfgWwUzByxK2PQX8L7UnrEuApQnLyXgP1ur079ae/QvzmbGojEvGDaz/AOeca0OiJqsLCGYK/ouk6m2rlwKDajvQzO5PWL6vwRG2IccNLeL+1z5g+6695Oc2e68C55xLG1HfWRVQVTuqLocIQy21FElXS5ov6T1JD4Xv1AZKmiVpiaRHJOWE++aG60vC7QMSznN9WL5I0sktEfuE4iJ2l1fwnyUbWuJyzjmXNqImq/eASbVsOwV4M8pJJGVLukbSa5I+lFRW/RMxntrO3wf4NnC4mY0kSKLnALcAt5nZYGATcGl4yKXAprD8tnC/ykF8zwFGEExlcmcNNcqkO3xAVzrkZjHd31s559xnRH3W9BPgCUntCKb5MOAQSV8EvgGcHvE8t4X7TyV499UcTdizgHaS9gDtgTXAccC54fb7gRsJOidPCpch6At2hySF5Q+H/b6WS1oCjAZeb4Z4P5WTlcGxBxTyUmkZFRVGhk/I6JxzQMRkZWZPSToX+CVBIwmAu4HVwPlm9nzE632Z4N3XrxscaQRmtlrSr4APCUbEeAF4A/g4HPEdghaNlfPI9wFWhsfulbSZYDbjPsDMhFMnHtOsThzeg2feXcM7qzdzSL+Clrikc86lvEjJSpLM7FHgUUkHAN2BjcAiMzNJHc0sSrPzylEsmkXY92sSQbP6jwlqgROb8XqXAZcBFBYWUlJS0uRzZu02MgT3/Gs2Zx2Q0+TzxWHbtm1JuRetgd+LKn4vqvi9aLiojwHvIxhoFjN7H3i/coOk7sBzwOERzvNn4KvAtAZFGd0JwPJw1Awk/QM4GiiQlBXWrvoS1AgJv/sBqyRlAZ2BjxLKKyUe8ykzmwxMBiguLrbx48cn5Uf8bcXrLN6xh/HjP5eU87W0kpISknUv0p3fiyp+L6r4vWi4qA0sTpK0zwC0knoBrzTgPOuA4yTNkPQDSZdX+3wz4nlq8yEwJhwNQ8DxBHNyzQC+FO5zIUHfMIAp4Trh9ulmZmH5OWFrwYHAEIJJKFvECcN6ULp2Kys37mipSzrnXEqLWrM6EZghaaeZXQ0QNvOeTpCATol4nt+G3/2BY2vYXteo7PUys1mSHidonbgXmEdQ83kGeFjST8Kye8JD7gH+Gjag2EjQAhAzmy/pUYJEtxe4wszKGxtXQ504vAc/eWYh0xas8w7CzjlH9AYW70k6CXhJ0icEwy1NI3gc+AUzizSgnZk1aCzCxjCzG9h3NI1lBK35qu+7k6DRR03n+SnByB0tbr9u+Qwp6sCLCz1ZOeccNGAgWzObR1CDugKYQ1BDOSVqonINc8LwHsxavpHNO/bEHYpzzsWurskXL69l0zMEjwVfBC4NXg1hZlbv4ztJp9a3j5k9W98+bcGJw3twV8lSSt4vY9IhLdJq3jnnUlZdjwHvqOfY2xKWo75rmhruW723a+LgtikzdFOcDulbQPcOOby40JOVc87VNflic7xfqukFTBfgZOBi4KJmuGZaysgQxw/twbPvrWH33gpyspr9dZ9zzqWsFv0LaGYf1PB5y8xuIRgR4wctGU+qO2F4D7bu3MucFRvjDsU552JV1zur4cBSM9sVLtfJzBY0MZZ5VI3T54Bxg7uTm5XBtAXrOHpw97jDcc652NT1zuo9YAxBZ9j3qH3SRIXbGv2uKZyy4yKCQWddqF1OJscM6c6LC9dxwxeGEzZmcc65NqeuZDWBoFNs5XKTSZrDvkkvBxgAdCR4b+USnDCsBy8uLKN07VaG9eoUdzjOOReLuhpYvFzTcnXh2ID1PiYMzWffZLWTYMDZf5rZ/IjnaTOOH9YD6V2en7/Wk5Vzrs1KxtzpxwKPEuExoJldlITrtSmFHXM5Yr+uPPfeWr5zwgFxh+Occ7Hw9tBp4OSRPSldu5UVG3ywEOdc2+TJKg2cPKIHAM/PXxtzJM45Fw9PVmmgb5f2HNinM895snLOtVGerNLExJE9mffhx6zdvDPuUJxzrsXVmqwkrZdUVt8HuLcF422zTh7RE4AXFnjtyjnX9tTVGvAP1N4RuFEkvUEw4eFDZrYpmedu7QYXdWBwUQf+9e5aLhg7IO5wnHOuRdXVz+rGZrjeu8AtwK8lTSFIXNPCqeRdPSaO6MmdJUvYuH03XfNz4g7HOedaTEsPZHsR0JNgAseewHPAh5J+KmlwS8aSjiaO7EmFwYsL1sUdinPOtagWb2BhZtvN7F4zOxYYAvwFOA9YJOkVSRdJymvpuNLBiN6d6FPQzlsFOufanLhbA5ZT9V6snGBQ3DuBFZJOjC2qFCWJiSN78uriDWzd6dPdO+fajhZPVpLaS7pQ0gxgCXA2QYLqZ2bHAH2B6cCfWjq2dHDKyJ7sLq9gemlZ3KE451yLadFkJekvwFqCloYfABPMbKiZ/dLM1gGY2UbgdwQjsbtqRvXvQs9OeTz9ts+m4pxrO1q6ZjUcuAboZWYXmdm/a9lvPkmalqS1ycgQpx7Yi1feX88WfxTonGsj6popeDkN6GdlZoPq2i4pF3gamGVmW+s51zag1mlJ2rrTDu7Fvf9ZzrT56zjrsL5xh+Occ82urprVE9U+2UBngpmDp4bfnQkS3uP1XcjMdgHXAwVNithxaL8C+hS0Y+o7/407FOecaxF1dQq+pnJZ0g+ApcDnzWx7QnkHgsS1JeL1ZgOj8FpTk0jitIN6cc+ry/l4x24K2nsHYedc6xb1ndUVwK2JiQo+fVz3q3B7FN8HLpf0LUmDJOWHrQM//UQPvW37/EG92FthPm2Ic65NiJqsOgE9atnWE+gQ8TyzgP2B24HFBDWyrdU+TSKpQNLjkkolLZQ0VlJXSdMkLQ6/u4T7StLtkpZIekfSqITzXBjuv1jShU2NK9kO7NOZ/l3bM/UdbxXonGv9ok5r/zRwq6QtwBQz2y0pB5hEMNbf0xHPcwlJHhy3Br8DnjOzL4Uxtgd+ALxkZr+QdB1wHXAtcArBKBpDgCOBu4AjJXUFbgAOD+N9Q9KUVBp8t/JR4J9eWcZH23bRrUNu3CE551yziZqsvgncBzwKmKStQEeCESemhNvrZWb3NTzE6CR1Bj4HXBRebzewW9IkYHy42/1ACUGymgQ8EA6kOzOslfUK950W9vlC0jRgIvBQc8bfUKcd1Js7S5byr/fW8rUx+8UdjnPONZtIycrMNgNflDQCOILgkeBaYI6ZLWjG+BpqILAe+Iukg4E3gKuAHmZW+bxsLVWPNPsAKxOOXxWW1Vb+GZIuAy4DKCwspKSkJGk/JAozo2e++NsrC+i7c3mLXrsu27Zta/F7kar8XlTxe1HF70XDRa1ZAWBm8wk67DaapLOBrwMHAPsMWGtmRU04fRZBa8MrzWyWpN8RPPJLPL9JSsqjSDObDEwGKC4utvHjxyfjtA3ylb3v8/vpixl+2BiKOqbG+L8lJSXEcS9Skd+LKn4vqvi9aLjII1hIKpJ0i6SXJC0Ka1lIukrS2IjnOJfgMdwSgjEApxA0fc8gaGxxR0N/QDWrgFVmNitcf5wgea0LH+8RflcOrLca6JdwfN+wrLbylHP6wb0wg6k+/JJzrhWLlKwkjSZovXcWsAIYDFS+0e8F/F/E630PuJmqpu53mtklBI/vNgA7Ip6nRma2FlgpqTgsOh5YQJAUK1v0XQg8FS5PAS4IWwWOATaHjwufB06S1CVsOXhSWJZyBhd1ZETvTjw5LyVzqXPOJUXUmtVtwAyCR3ffIGhYUWk2MDrieYYA/zGzcoIpQToBhMMv3QJ8K+J56nIl8KCkd4BDgJ8BvwBOlLQYOCFcB3gWWEZQ0/szcHkYz0aCpDon/NxU2dgiFX3x0D68u3ozS8qa3PLfOedSUtR3VqOASWZWIUnVtn0ERH3PtIWqGtlqYBhByzwIEmC3iOeplZm9RdDkvLrja9jXqKVDs5ndC9zb1HhawumH9OZnzy7kyXmr+d7JQ+MOxznnki5qzWozUFjLtkFA1HnW5wAHhctTgB9J+nrY6fZWYGbE87gERR3zOGZIIf+c918qKpq7G5tzzrW8qMlqCvBjSYkjq5uk7gRTfvwj4nl+DnwYLv+I4BHiXQRT228geMToGuHMUX1Y/fEnzF6Rsk8rnXOu0aImq2sJHuEtAF4Jy/4ILAI+IUg89TKzmWb2SLj8sZlNAvKBAjM70syWNSR4V+Wk4T3Jz8nkyTe9oYVzrvWJlKzCYYbGELzf+QB4EVhO0Ifp6PrmpwKQlCdpl6Qzqp17l5lFHbXd1aJdTiYTR/bi2XfXsHNPedzhOOdcUkXuZ2Vmu83sHjM718xOMrNzzOzP4TxVUY7fSdC/aW9jg3V1O3NUH7bu2stLC8vq39k559JI1H5W0yXdJaldDduGSZoe8Xp/Ar4tKbshQbpoxgzqRo9OuTw5b1XcoTjnXFJFbbo+nuAx4FhJXzSzxIHoOgHHRjxPATASWCHpJYJWhInN18zMro14LldNZoY445A+3PPqch+J3TnXqkR+DAh8FdhGMF3GqY283lnALmA3cAzwJeDL1T6uCc4c1Ze9FcZTb/mU98651qMhyWoNQQ3rEWCKpEgtABOZ2cB6PoPqP4urS3HPjhzctzOPzl1J0OfZOefSX0OSFWa218y+SdAf6npJTwGdox4v6QJJNY5SEc7me0FD4nE1+8oR/Shdu5V3Vm2OOxTnnEuKBiWrSmZ2D0EtaxRBTSuqvxBMa1+TgeF210SnH9ybdtmZPDxnZf07O+dcGoiarO4nmNTwU+E0HIcTDKH0YU0H1aD6uIKJuhF0PHZN1DEvm1MP7MXTb/+XHbu9p4BzLv1FnSn44lrK1xFMn1GrcEr5SQlF/0/S+mq75RE0uJgTJR5Xv3NG9+OJN1fxzDtr+PLh/eo/wDnnUlityUpSezPbUblc34kq961BEXBgwvr+QM9q++wGXgB+Ut91XDSH79eFQYX5PDp3pScr51zaq6tmtVXSWDObTdBkvb6mZZk1FZrZnwnmikLSDOByM1vYmGBddJI4+/B+/PxfpSwp28bgog5xh+Scc41WV7K6BFiasNzkdtBmNqGp53DRnTmqL7c+v4hH567kB6cOizsc55xrtFqTlZndn7B8X2MvIOly4DEzWx8u18XM7K7GXst9VmHHXE4Y1oPH5q7kuyceQF52jZVf55xLeVGHW2qKO4C5BK0J76hnXyOY38olyflj9+O5+Wt59t01nDmqb9zhOOdco9TVwGIODXj0Z2ajaynPqGnZtYyj9u/GoMJ8Hnj9A09Wzrm0VVfNaj5JeE/l4iWJ88fsx4+fXsC7qzZzYN/IA44451zKqOud1UXNcUFJOcBFwGigF8GYg7OA+81sd3Ncs60767C+/PK5Rfxt5gfc8qWD4g7HOecarEUfy0kaBiwG/kAwVUh5+P0HYImk4S0ZT1vRKS+bMw7tw1Nvr2bzjj1xh+Occw0WuYGFpAHA14ADCEac+Awz+0qE00wGNgPHmNmnQzRJ6g9MBf4IfC5qTC6688fsx0OzP+SxN1byP8f44PbOufQSdabgwwjeYZ0XfoYQjAv4JYJJGbtHvN7hwI8SExVAuH4DcETE87gGGt67E4fv14W/zfyAigp/FemcSy9RHwPeCjxG8MhOwKXh3FPjCBph/DLieVZQQ60slEf0AXFdI1xw1ABWfLSDGYvK4g7FOecaJGqyOgR4CKgI1/MAzOw14MfALyKe5zrgJ5KOTCyUNAa4GfAp7ZvRKSN70rtzHnf/e3ncoTjnXINETVYG7LZg6tkyYL+EbSsJHgvWSNIcSbMlzQZ+CHQCXpO0RtLbktYA/wnLf9CYH1HDNTMlzZM0NVwfKGmWpCWSHglbJCIpN1xfEm4fkHCO68PyRZJOTkZcccvOzOCiowfw+rKPeG+1T8zonEsfUZPVAqomTXwduFrSEEn7Ad+nagzBmsyv9nkGeAB4Dngz/H4gLJ/f0B9Qi6uAxMFybwFuM7PBwCbg0rD8UmBTWH5buB9hq8RzgBHAROBOSa1irKJzRvcnPyeTe1712pVzLn1EbQ04mara1A8IpvMoDde3EzS0qFFz9deqjaS+wOeBnwLflSTgOODccJf7gRsJhnWaFC4DPA7cEe4/CXjYzHYByyUtIegX9noL/Yxm0ykvm7OP6M8Dr6/g+xOL6dW5XdwhOedcvaJOvvjXhOWFYX+powjeXc00s3rf2EvKI2i2fraZ/bNx4UbyW4LaXsdwvRvwsZlVTpm7CugTLvcheIyJme2VtDncvw8wM+Gcicd8StJlwGUAhYWFlJSUJPN3NJthmRWUVxg3P/xvvlKck/Tzb9u2LW3uRXPze1HF70UVvxcN16iBbM1sG0HtqiHH7JRUBjTbPOuSTgPKzOwNSeOb6zqVzGwyQa2T4uJiGz++2S+ZNCWb3uTfi9dz60XjyM9N7njGJSUlpNO9aE5+L6r4vaji96LhGtIpOI+gw25f9m1+HnVqjz8B35b0vJk1x1AKRwOnSzqVIMZOwO+AAklZYe2qL7A63H810A9YJSkL6Ax8lFBeKfGYVuHSYwbyzLtreHjOSi4dNzDucJxzrk6RkpWkccATQGEtu0Sd2qOAoK/WCkkvAev47GC5ZmaNbr5uZtcD14cxjweuMbPzJD1G8F7tYeBC4KnwkCnh+uvh9ulmZpKmAH+X9BugN0Frx9mNjSsVjerfhTGDujL5laV8bUx/crNaRfsR51wrFbU14O3AMuBQINfMMqp9ov6lOwvYBewGjiFIEF+u9mkO1xI0tlhC8E7qnrD8HqBbWP5dgn5gmNl84FGCVpDPAVeYWXkzxRabK48bwrotu3j8jVVxh+Kcc3WK+hiwGDjTzN5uysXMrMWeN5lZCVASLi8jaM1XfZ+d1JIgzeynBC0KW62j9u/GIf0KuKtkKV85vB/ZmT7dmHMuNUX96/QO0LM5A3EtTxJXHjeYVZs+4am3/ht3OM45V6uoNatvAvdJWmFmLzf1ouE7sNpGb7+zqed30R03tIhhvTpx54wlfPHQPmRmKO6QnHNuH1GT1TSgPTBd0m5ga/UdzKyovpNI6gG8BAwnaFhR+ZcxsZGFJ6sWVFm7uvzBN3n23TV84eDecYfknHP7iJqs/kByprj/NUHH4H4EnXGPJGgR+DXgAoKRJ1wLmziiJ4OLOnD7S4s59cBeXrtyzqWcqCNY3Jik6x1LMG7fmnBd4VxWP5OUQVCrahWDxqaTjAxx9QkHcMXf3+Spt1Zz5qi+cYfknHOf0dLNvwqA9WZWAWwBEh8dvkYwhJOLwSkjezKidydue/F9du+tqP8A55xrQbXWrCQ9ClxvZkvD5TpFnNZ+OdArXK6ceXhquP4FYGOEc7hmkJEhrjm5mIv/ModH5nzI+WMHxB2Sc859qq7HgIVAdrhcRHLeWT0DnETQ4fYnwFOSVgF7gP745IuxGn9AIaMHdOX26Uv40mH9aJfjo1o451JDrcnKzCYkLI9PxsXC4ZAql/8l6Sjgi0A7YJqZ/SsZ13GNI4nvTSzmy398nfteW8E3x+9f/0HOOdcCkjvcdgOZ2VxgbpwxuM86YkBXJhQX8seXl3Lu6P50bp9d/0HOOdfMog5k+6M6Nlc2lni7pg7DkgaaWYOmpQ1bBvYxs5UNOc4lx/cnDuXzt/+b3770Pjd8YUTc4TjnXOSa1ZUEo03kh+vbgA7h8vbwPLmS3gJOMbN1Cce+LelpgkFjZ5hZre++wll+zwauIBjF/daI8bkkGtarE2cf0Z+/vv4B5x25H4OLOtR/kHPONaOoTddPJegbdTbQzsw6EbxnOicsP4FgrqtCgo6/iYYR1LyeBsokPSPpFknfk/RtST+S9BdJ7wArCDoIf9/MPFHF6P9OOoB22Zn89JkFcYfinHORa1Z3AL8ws8cqC8xsF/CopI7A781slKSfELTyI2G/1cA3JX2fINkdB5xBMDBuHkFz9UUEyexiM3ujaT/JJUP3DrlcefxgfvZsKSWLyhhfXO9oWs4512yiJquDgLW1bFtDUHsCKAU61rSTmW0F7g4/Lg1cdNRA/j7rQ37yzEKOHtzdpxBxzsUm6l+f94GrJOUkFkrKBa4mqBlBUFtah2sVcrIy+OHnh7OkbBsPvP5B3OE459qwqDWrqwg69K6SNA1YT/B+6kSCRhenhvsdCvwj2UG6+JwwrIjxxYX85oVFnDKyJ70L2sUdknOuDYpUswpn3R0C3A/0JhhstjdwHzCkssm6mV1nZlc3S6QuFpK4edJIys24Ycr8uMNxzrVR9daswkd91wBTzex7zR+SSzX9urbn6hMO4Of/KuX5+Ws5eYRPGu2ca1n11qzCVn8/JBgx3bVRl4wbyNCeHbnhqfls3bkn7nCcc21M1AYWs4BRzRmIS23ZmRn8/MwDWbd1J796flH9BzjnXBJFTVbfBy6X9C1JgyTlS2qf+GnOIF1qOLR/Fy4cO4D7X/+A15ZsiDsc51wb0pCa1f7A7cBighEptlb7uDbg2olDGdQ9n2see5st/jjQOddCojZdv4TkzGfl0ly7nEx+c/YhnHXXa/x4ygJ+/ZWD4w7JOdcGREpWZnZfM8fh0sgh/Qq4Yvz+3D59CScO78HEkd460DnXvHz8HNcoVx4/hJF9OvGDJ99l3ZadcYfjnGvlIicrSWdLelHSh5LKqn+aM8ioJPWTNEPSAknzJV0VlneVNE3S4vC7S1guSbdLWiLpHUmjEs51Ybj/YkkXxvWbUlV2Zga/PfsQPtldzpUPzWNveUXcITnnWrFIyUrSuQSjVywB+gJTgKnh8VsIRmVPBXuB/zOz4cAY4ApJw4HrgJfMbAjwUrgOcArByBxDgMsI5tBCUlfgBuBIYDRwQ2WCc1UGF3XkZ2eOZPbyjfx62vtxh+Oca8Wi1qy+B9xMMCkiwJ1mdgkwENgA7GiG2BrMzNaY2Zvh8lZgIdAHmESQbAm/zwiXJwEPWGAmUCCpF8FwUtPMbKOZbQKmARNb7pekjy8e2pevju7PXSVLmV7qYxg755pH1GQ1BPiPmZUD5UAn+DQh3AJ8q3nCazxJAwgG1p0F9DCzNeGmtUCPcLkPsDLhsFVhWW3lrgY3fGE4w3t14upH3mbVppT4d4tzrpWJ2nR9C5AbLq8mmL+qJFwX0C25YTWNpA7AE8B3zGyLpE+3mZlJSkozfEmXETw+pLCwkJKSkmScNi1dNKSCG17bw1fvfJnvjCxv0/ci0bZt2/xehPxeVPF70XBRk9UcggkYnyd4X/UjSXuB3cCPgJnNE17DScomSFQPmlnldCXrJPUyszXhY77KBiGrgX4Jh/cNy1YD46uVl1S/lplNBiYDFBcX2/jx46vv0qb0GLyeS+6bw9+WZvH4d44lM0P1H9TKlZSU0Nb/u6jk96KK34uGi/oY8OfAh+Hyj4DZBI0R/kLwzuobyQ+t4RRUoe4BFprZbxI2TQEqW/RdCDyVUH5B2CpwDLA5fFz4PHCSpC5hw4qTwjJXh2MPKOSGLwxnXlk5tzxXGnc4zrlWJGqn4JmEtScz+xiYFE4dkmtmW5ovvAY7GjgfeFfSW2HZD4BfAI9KuhT4APhKuO1ZgokjlxA0ErkYwMw2SrqZoEYJcJOZbWyRX5DmLhg7gJfnLWLyK8sY1D2fc0b3jzsk51wrEPUx4D7CqUN2JTGWJjOzVwneodXk+Br2N6paOFbfdi9wb/KiazvOHZrDnrwu/PCf79E1P4eTfP4r51wT1ZqsJP2oAecxM7s5CfG4ViAzQ9x13ijOu3sW3/r7PO67+AiOGtw97rCcc2msrprVjcAnwHZqr61UMoJ+WM4BkJ+bxX0XH8HZf5rJ/zwwl79/fQyH9CuIOyznXJqqq4HFUiAbeINgWvtBZlZYy6eoRaJ1aaWgfQ5/vXQ03TvkctFfZvPe6s1xh+ScS1O1JqtwaKKjgPkEtaZ1kv4h6cuS2rVUgC69FXXK48H/OZL8nCy++ueZvPnhprhDcs6loTqbrpvZXDO7xsz6Eww3tJZgHMAySQ9K+lxLBOnSW7+u7Xn0f8fSLT+H8++excxlH8UdknMuzUQedd3MXjGzywk60f4ROBv4TjPF5VqZPgXtePQbY+ld0I4L753NjNKUGKjfOZcmGjJFyNGSfk/QT+mbwOPA75orMNf6FHXK4+HLxjCkRwcuvX8Of535QdwhOefSRJ3JStIoSb+U9AHB1Br9gKuBIjM7x8xebokgXevRrUMuj1w2lgnFRfy/f77HzVMXUF6RlKEanXOtWK3JStIiglErDiKY26nIzM4ws4fNzIfWdo2Wn5vF5AsO56KjBnDPq8v5xl/nsvmTPXGH5ZxLYXXVrIYQTGZ4GPBLYElNMwSn0kzBLn1kZogbTx/Bj08fQcmi9Xzh969603bnXK3q6hT84xaLwrVZFx41gJF9OnHFg/M4867XuOn0EZx9RD8Sp3Vxzrlak5WZebJyLeKw/boy9dvj+M7Db3HdP97l5ffX85MzRtKtQ279Bzvn2oTIrQGda07dO+Ry/yWjuXbiUF5aWMbJv32FF+avjTss51yK8GTlUkZmhvjm+P2ZcuXRFHbM47K/vsG3H5pH2ZadcYfmnIuZJyuXcob27MRTVxzNVccP4bn31nL8r1/m3leXs7e8Iu7QnHMx8WTlUlJOVgZXn3gAL1z9OUbt14Wbpi7gtN+/yozSMoJpyJxzbYknK5fSBnTP576Lj+CPXxvFJ3vKufi+OZw9eSZvfOAD4jrXlniycilPEhNH9uLF7x7LzZNGsGz9ds666zW+dvcs/rNkg9e0nGsDPFm5tJGdmcH5Ywfw8vfGc+3EoSxat5Xz7p7F6Xf8h2ffXePvtJxrxerqFOxcSsrPzeKb4/fn4qMH8I83V/OnV5Zy+YNv0rNTHl85oh/nHNGP3gU+5ZpzrYknK5e28rIzOffI/px9RD9eXLiOv8/6kN9PX8wd0xczvriILx7ah+OHFdE+x/8zdy7d+f/FLu1lZoiTR/Tk5BE9WblxB4/MWcljb6xkemkZ7bIzOX5YEacd1JvPHdDdE5dzacr/z3WtSr+u7bnm5GKuPvEA5qzYyNNv/5d/vbeWqe+sIScrgyMHdmVCcREThhYxsHt+3OE65yLyZOVapcwMMWZQN8YM6saNp49g1rKNzFhUxoxFZdw0dQE3TV1A/67tGT2wK6MHdOWIgV0Z0K29D6DrXIryZOVavezMDMYN6c64Id35f6cN58OPdlDyfhmvLt7ASwvX8fgbqwAo7JjLof0KGNG7MyP7dGJE78706JTrCcy5FODJyrU5/bu154KxA7hg7AAqKoyl67cxe8VG5izfyDurNjNt4Toqu25175DDsF6d2L+wA4MK8xnUPfju1TnPk5hzLciTlWvTMjLEkB4dGdKjI+cduR8A23btZeGaLcxfvZn3/ruF0rVbeHTuSnbsLv/0uHbZmQzonk+fgjx6F7SjV+d29A6Xexe0o6hjLtmZ3o3RuWTxZFULSROB3wGZwN1m9ouYQ3ItpENuFkcM6MoRA7p+WmZmrNuyi2Xrt7F0w3aWrd/Gig3bWbXpE2Yv38iWnXv3OU/ndtm0y9hLv9LX6JafS7cOOXTLz6Frfg4d87LpmJf16XfndsF3h9wssjzJObcPT1Y1kJQJ/AE4EVgFzJE0xcwWxBuZi4skenbOo2fnPI4a3H2f7dt27WXNx5/w3807+e/Hn1C2ZRcfbd9F6fJVZGSIpeu3MWfFbjbu2E19o0O1z8kkPzeLvOwM2mVn0i47k9zwu7IsL+GTk5VBdobIyswgO1NkZ2aQlSmyMzLIzhJZGUF5VkYG2Qn7ZmaIDEGGRGaGULgcrAe/OVgPy8P9M6Vw22fLMxScA0AovG9Vv2t3ubFrb/mn26v2rbrHn13/bLlr2zxZ1Ww0sMTMlgFIehiYBHiycjXqkJv16ePERCUlGxg/fuyn6+UVxsc7drN1597ws4ct4Xdi2fbde9m5p4Kde8r5ZE85O/eU8/Ene9i5uZyde8v5ZHf5p9v2lKfR2IjTnmvyKaQICY7PZsLq2+NWUVFB5ktNvxdNlSr3IwpPVjXrA6xMWF8FHFl9J0mXAZeFq7skvdcCsaWD7sCGuINIEX4vqvi9qOL3okpxlJ08WTWBmU0GJgNImmtmh8ccUkrwe1HF70UVvxdV/F5UkTQ3yn7+Jrdmq4F+Cet9wzLnnHMx8GRVsznAEEkDJeUA5wBTYo7JOefaLH8MWAMz2yvpW8DzBE3X7zWz+fUcNrn5I0sbfi+q+L2o4veiit+LKpHuhXyWVeecc6nOHwM655xLeZ6snHPOpTxPVkkgaaKkRZKWSLou7njiIuleSWXe3wwk9ZM0Q9ICSfMlXRV3THGRlCdptqS3w3vx47hjipOkTEnzJE2NO5a4SVoh6V1Jb9XXhN3fWTVRODTT+yQMzQR8tS0OzSTpc8A24AEzGxl3PHGS1AvoZWZvSuoIvAGc0Ub/uxCQb2bbJGUDrwJXmdnMmEOLhaTvAocDnczstLjjiZOkFcDhZlZvB2mvWTXdp0MzmdluoHJopjbHzF4BNsYdRyowszVm9ma4vBVYSDAySptjgW3hanb4aZP/SpbUF/g8cHfcsaQbT1ZNV9PQTG3yj5KrmaQBwKHArJhDiU346OstoAyYZmZt9V78Fvg+UBFzHKnCgBckvREOX1crT1bONSNJHYAngO+Y2Za444mLmZWb2SEEo8GMltTmHhNLOg0oM7M34o4lhYwzs1HAKcAV4auEGnmyajofmsnVKHw/8wTwoJn9I+54UoGZfQzMACbGHEocjgZOD9/TPAwcJ+lv8YYULzNbHX6XAU8SvFapkSerpvOhmdw+wkYF9wALzew3cccTJ0mFkgrC5XYEjZFKYw0qBmZ2vZn1NbMBBH8nppvZ12IOKzaS8sPGR0jKB04Cam1J7MmqicxsL1A5NNNC4NEIQzO1SpIeAl4HiiWtknRp3DHF6GjgfIJ/Pb8Vfk6NO6iY9AJmSHqH4B9308yszTfbdvQAXpX0NjAbeMbMap3ky5uuO+ecS3les3LOOZfyPFk555xLeZ6snHPOpTxPVs4551KeJyvnnHMpz5OVc865lOfJyjnnXMrzZOVcKyPpSkkm6ZIatnWWVCFpehyxOddYnqyca30OC79rGjB1FKBatjmXsjxZOdf6jAJ2AjUN+1WZyN5suXCcazpPVs61IpLygGHAO+G4ldXVVetyLmV5snKudTkIyKL2ZHQYsBVY3GIROZcEnqyca11Ghd/7JCtJnYHBwDzzEaxdmvFk5VzrUmuyIpi2RPj7KpeGPFk517pUJquaZqs+J/z291Uu7Xiycq6VkJQNHBiujqu27StA5ay077ZkXM4lg0++6FwrIelQgkd864DOwJPAR8DBwKHAdoLZWf8B3GpmM2MK1bkGy4o7AOdc0lQ+ArwJGARcBLQDZgETgGPDbfsBH8YQn3ON5jUr51oJSX8ALgeONLPZccfjXDL5OyvnWo9RQDn+Tsq1Ql6zcq4VkJQJbAGWm9nIuONxLtm8ZuVc6zAUaA+8FXMczjULr1k555xLeV6zcs45l/I8WTnnnEt5nqycc86lPE9WzjnnUp4nK+eccynPk5VzzrmU58nKOedcyvv/z7Rd7Hb3UaYAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# Plot the merginalized likelihhod\n", "ax33 = plt.subplot(111)\n", "ax33.plot(muScan, lik)\n", "ax33.set_xlabel('$\\mu$', size=20)\n", "ax33.set_ylabel('Marginalized Likelihood\\n (arbitrary units)', size=15)\n", "ax33.set_ylim(bottom=0)\n", "ax33.set_xlim(muScan[0],muScan[-1])\n", "ax33.grid()" ] }, { "cell_type": "code", "execution_count": 9, "id": "eccc4ac2", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.24 1.2464608820756984 1.25\n", "0.9490595948215217 0.9505151316467726\n" ] } ], "source": [ "# find the 95% point\n", "# i is the index of the first element past the 95% point\n", "cumlik = np.cumsum(lik)/lik.sum()\n", "i = np.nonzero(cumlik>0.95)[0][0]\n", "# do some interpolation\n", "dlik = (0.95-cumlik[i-1]) / (cumlik[i] - cumlik[i-1])\n", "dmu = muScan[i] - muScan[i-1]\n", "mu95 = muScan[i-1] + dlik*dmu\n", "print(muScan[i-1], mu95, muScan[i])\n", "print(cumlik[i-1], cumlik[i])" ] }, { "cell_type": "code", "execution_count": 10, "id": "5561b4cd", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbAAAAEPCAYAAAAj0pGKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABBy0lEQVR4nO3dd3gVZfbA8e9JJxASSoKQ0EuoUqUIalBRrIi9g2JZ27r29lu7rrt2V1CxoquIq7gCYkEliIUiCkonCAhI76EEkpzfHzMhl5AyKTeT3JzP88xz574zd+bc12sOM/MWUVWMMcaY6ibM7wCMMcaYsrAEZowxplqyBGaMMaZasgRmjDGmWrIEZowxplqK8DuAUJGQkKBt2rTxO4wqYffu3dSuXdvvMKoEq4t8Vhf5rC7yzZkzZ7OqJpbls5bAKkijRo346aef/A6jSkhPTyctLc3vMKoEq4t8Vhf5rC7yiciqsn7WbiEaY4yplqplAhORN0Rko4jMDyh7UkQWi8ivIvKxiCQEbLtHRDJEZImInBxQPtgtyxCRuwPKW4rITLd8nIhEVdqXM8YY40m1TGDAW8DgAmVTgM6qeiSwFLgHQEQ6AhcCndzPjBKRcBEJB0YCpwAdgYvcfQH+CTyrqm2AbcCI4H4dY4wxpVUtn4Gp6rci0qJA2ZcBb2cA57rrQ4D3VTULWCEiGUBvd1uGqv4OICLvA0NEZBFwPHCxu88Y4EHgpSB8FWNC3oEDB1izZg3x8fEsWrTI73CqhJpYFzExMaSkpBAZGVlhx6yWCcyDK4Fx7noyTkLLs8YtA1hdoLwP0ADYrqrZhex/CBG5BrgGIDExkfT09IqIvdrLzMy0unBZXUCdOnVo1KgRLVq0ICIiVP/klE5OTg7h4eF+h1FpVJUdO3Ywb948MjMzK+y4IfdrEpH7gGzg3WCfS1VHA6MBUlNT1VoVOayFVT6rC1i0aBEpKSlkZmYSFxfndzhVwq5du2pcXcTFxZGZmUmvXr0q7JghlcBEZDhwOnCC5g+zvxZoGrBbiltGEeVbgAQRiXCvwgL3N8aUgYj4HYLxWTB+A9W1EcdhRGQwcCdwpqruCdg0AbhQRKJFpCXQFpgFzAbaui0Oo3AaekxwE99U8p+hDQM+Ken8OTYrjTHGVKpqmcBEZCzwI5AqImtEZATwIhAHTBGRuSLyMoCqLgA+ABYCnwM3qGqOe3V1I/AFsAj4wN0X4C7gVrfBRwPg9ZJiWr87lwM5uRX6PY0xxhStWt5CVNWLCikuMsmo6mPAY4WUTwYmF1L+O/ktFT05kAtvfr+Ca45tXZqPGWOMKaNqeQVWFdWKEJ77ahnrduz1OxRjTCGef/55OnfuTKdOnXjuuecOlrdo0YIuXbrQrVu3gw0MNm3axIABA+jcuTP/+9//Du47ZMgQ/vzzzyLPUadOnYPrRx99dKniy9t/+/btjBo1qlSfLa8rr7ySpKQkOnfuXOj21atXM3DgQDp27EinTp14/vnnD9leWB1WBktgFaRBjJCTqzwyaaHfoRhjCpg/fz6vvvoqs2bNYt68eUyaNImMjIyD26dOncrcuXMPjmc6duxY/vKXvzBr1qyDyW7ixIl0796dJk2aeDrnDz/8UKoY8/b3I4ENHz6czz//vMjtERERPP300yxcuJAZM2YwcuRIFi489G9dwTqsDJbAKkhEGNw4sA2Tf1vPt0s3+R2OMSbAokWL6NOnD7GxsURERHDccccxfvz4IvePjIxkz549ZGVlER4eTnZ2Ns899xx33nmn53PWqVOHlStX0r59e4YPH067du245JJL+Oqrrxg0aBBt27Zl1qxZh+wPcPfdd7N8+XK6devGHXfcUew50tLSWLx4MQBbtmwp8gqqJMceeyz169cvcnvjxo3p0aMH4DSH79ChA2vX+t842xJYBbrmuFa0bFibByYsICs7x+9wjKma0tIOX/KuOPbsKXz7W2852zdvPnybB507d2b69Ols2bKFPXv2MHnyZFavdsYxEBFOOukkevbsyejRowG4+OKL+eSTTxg0aBD33nsvo0aN4rLLLiM2NrbUXzcjI4PbbruNxYsXs3jxYt577z2+/PJLnnrqKR5//PHD9n/iiSdo3bo1c+fO5cknnwTg1FNPLfTWZUZGBu3atQPg119/pUuXLoXGcMwxx9CtW7fDlq+++qrU32flypX88ssv9OnT52BZYXVYGaplI46qKjoinIfO7MTlb8xi9LTfuemEtn6HZIwBOnTowF133cVJJ51E7dq16dat28GRML777juSk5PZuHEjgwYNon379hx77LF8+umnAGzbto0nnniCjz/+mKuvvppt27Zx22230a9fP0/nbtmy5cHE0qlTJ0444QREhC5durBy5UpPx5g8+bC2ZqxatYrk5GTCwpzrkF9//ZUjjzyy0M9Pnz7d03lKkpmZyTnnnMNzzz1H3bp1D5YXVYfBZgmsgh3bLpHTujTmxakZnNU9mab1S/8vNmNCWnFDa8XGFr+9YcPitxdjxIgRjBjhjMt97733kpKSAkBysjNSXFJSEkOHDmXWrFmH/PF95JFHuO+++xg7diwDBgzg3HPP5eyzz+aLL77wdN7o6OiD62FhYQffh4WFkZ2dXdTHSjRv3rxDEtacOXO44IILCt33mGOOYdeuXYeVP/XUU5x44omeznfgwAHOOeccLrnkEs4+++xDtpVUh8FitxCD4P9O70B4mPDghAUl72yMqRQbN24E4I8//mD8+PFcfPHF7N69++Af9t27d/Pll18e8hxp2bJlrFmzhrS0NPbs2UNYWBgiwt69wWttHBcXV2iyKWju3Lns27fvYJyffPJJkbcQp0+fzty5cw9bvCYvVWXEiBF06NCBW2+99ZBtJdVhMFkCC4LG8bX424lt+XrxRqYs3OB3OMYY4JxzzqFjx46cccYZjBw5koSEBDZs2MCAAQPo2rUrvXv35rTTTmPw4PyZmu677z4ee8zpQnrRRRfx0ksvcdRRR3HzzTcHLc4GDRrQv39/OnfufLARR2HPwObNm0dubi5du3bl4YcfpmPHjowZM6ZM57zooovo168fS5YsISUlhddfz+9We+qppzJ9+nTeeecdvvnmm4PPz/Jua5ZUh8Ek+UMGmvJITU3VJUuWHHx/ICeX016Yzu6sHL669ThqRdWckadtANt8VhdOC8AOHTrUyAFsi1IRddG2bVt+/vnnalWneb+FQCIyR1XL1HnMrsCCJDI8jEeGdGbt9r2MnJpR8geMMcajXbt2ISLVKnkFgyWwIOrTqgFndWvC6G9/Z9WW3X6HY4wJEXFxcSxdutTvMHxnCSzI7j6lAxHhwqOf1qzZV40xJtiKTGAikisiOV6Xygy6OjkiPoYbj2/DlIUbmGYjdBhjTIUprh/YX4G8Fh6RwG1AJs7cWBuBRsAQoDbwdBBjrPZGDGjJuNmreWjiAj6/+ViiIuzC19Qs1ljMBOM3UORfUlV9UVVHqupIoBkwE+ioqner6jOqehfQwS1vWeGRhZDoiHDuP70jv2/azds/rvQ7HGMqVUxMDFu2bLEkVoOpKlu2bCEmJqZCj+t1JI7LgUu0wC9QVVVEXgXeA4LXMSIEHN8+ibTURJ77ahlndmtCUlzF/oc0pqpKSUlhzZo1bN++vcL/gFVX+/btq3F1ERMTc3D0k4riNYGF41xtFTZ2SiesMUiJRIS/n96Rwc99y5OfL+HJ87r6HZIxlSIyMpKWLVuSnp5O9+7d/Q6nSrC6qBheE8+7wOMicruItBORBPf1DpyZjt8NXoiho3ViHa7s35L/zlnD3NXb/Q7HGGOqNa8J7FbgFeBhYBGwxX19yC2/teiPmkA3Ht+GxLhoHpiwgNxceyZgjDFl5SmBqep+Vb0FSAFOAC4GjgdSVPVvqro/iDGGlLiYSO4e3J55q7fz0c9r/A7HGGOqrVI9u1LVraqarqrjVHWaqm4NVmChbGj3ZLo3S+Cfny9h574DfodjjDHVkucEJiKtROQlEflNRNa6r6NEpFUwAwxFYWHCg2d0YsvuLP799TK/wzHGmGrJUwITkZ7AXOAcYDbwtvt6DvCLiPQIVoChqmvTBM7v2ZQ3v19JxsZMv8Mxxphqx+sV2FPAL0ALVb1SVe9R1StxOjD/4m43pXTH4FRqRYXz8KSF1snTGGNKyWsC6w38S1X3BBa6758C+lR0YDVBwzrR3HxCW75duomvF230OxxjjKlWvCawvUCDIrbVB/ZVTDg1z7CjW9A6sTaPfLqQrGwbE9kYY7zymsA+BZ4QkQGBhe77fwATKzqw4ojIGyKyUUTmB5TVF5EpIrLMfa3nlouIvCAiGSLya+DzOhEZ5u6/TESGBZT3dBupZLiflWB9l8jwMO4/oxOrtuzhje9WBus0xhgTckrTkfl3YJqIrBOReSKyDpgGrMAZqb4yvQUMLlB2N/C1qrYFvnbfA5wCtHWXa4CXwEl4wAM4tz97Aw/kJT13n6sDPlfwXBXquHaJnNghiRe/WcbGnXYxa4wxXnjtyLxFVQcApwGjgO/d11NUdYCqbglijIXF8y1QsA/aEGCMuz4GOCug/G11zAASRKQxcDIwxe3btg2YAgx2t9VV1Rnu4MVvBxwraP7vtI4cyFGe+HxxsE9ljDEhwetgvgCo6ufA50GKpbwaqeo6d309znxlAMnA6oD91rhlxZWvKaT8MCJyDc5VHYmJiaSnp5frC5zYLJzxP6+lU9QWWieEl+tYfsrMzCx3XYQKq4t8Vhf5rC4qhucEJiIROP2+BuA03NgKTAfGq2p2cMIrG3eal6C3S1fV0cBogNTUVE1LSyvX8Xr1y+anp9KZsCaaj8/sT1hY0B69BVV6ejrlrYtQYXWRz+oin9VFxfDakTkJ+AkYi3MbsZX7+j4wW0QSgxahdxvc23+4r3nt0tcCTQP2S3HLiitPKaQ86OpER3DX4PbMW7PDxkk0xpgSeG3E8QxOM/q+qtpKVfupaiucBhAN3O1+mwDktSQcBnwSUH652xqxL7DDvdX4BXCSiNRzG2+cBHzhbtspIn3d1oeXBxwr6IZ2T6ZbU2ecxF02TqIxxhTJawI7FbhLVWcFFqrqbOAenKuxSiMiY4EfgVQRWSMiI4AngEEisgw40X0PMBmnBWUG8CpwvRv7VuARnCGxZgMPBwxOfD3wmvuZ5cBnlfG9wB0n8cxObM7M4sVvMirrtMYYU+14fQYWDewqYtsuIKpiwvFGVS8qYtMJheyrwA1FHOcN4I1Cyn8COpcnxvLo1jSBc3um8Mb3K7iwdzNaNqztVyjGGFNleb0CmwHcJSKH/CV139/lbjcV6M7BqURHhPPopIV+h2KMMVWS1wR2G9AJWC0i74vI8+5tvNVARyq/I3PIS4qL4abj2/D14o2kL7FxEo0xpiCvHZnn4oxIMRpIBAYBScDLQFtVnResAGuyK/q3pGXD2jw8aSH7s3P9DscYY6oUz/3AVHUz+cMzlYmIdMBpYNEbOAKIwelPthT4DvhSVfeW5xyhJCoijL+f3oEr3/qJt39cyVXH2NyhxhiTx/OMzGXlNl+/XERmAwuAB3H6Vm0HVgHhOOMVfgSsF5HRItIy2HFVF8e3b0RaaiLPf7WMTbuy/A7HGGOqDK8dmSNF5HYR+UFE/nBHgj9kKebji4B7gf8BHVW1gaoOVNXzVPVSVT1dVTsDdYHhQB3gNxG5tHxfLXT8/fSO7D2Qw1NfLPE7FGOMqTK83kJ8FrgWmARMBfaX4hz34Qw3VezQTu7kmB8DH4tICoeOhlGjtU6swxX9W/Dadyu4tG9zuqTE+x2SMcb4zmsCOw+4W1WfLu0JVPWjMnxmDYcOqFvj3XRCWz7+ZS0PTlzAh3/pRxCnKDPGmGrB6zMwAX4t78lEJEJEoguUnSQifwucaNIcrm5MJHee3J45q7bxydw//Q7HGGN85zWBvQoUNfpFaYzDnVASQET+ijM9yz+AGSJyegWcI2Sd2zOFI1Pi+cdni9idVaUmADDGmEpX5C1EEbk+4O164BIRmYoz8eP2Arurqr5EyfoCNwe8vwN4WlXvEJFROM/LJnkJvCYKCxMeOKMT57z0A6PSM7jj5PZ+h2SMMb4p7hnYi4WUNQOOK6RcCbiyKkYDnGSIiHQBmuB0hgb4L3CJh2PUaD2b12No92Renb6CC3o1o1mDWL9DMsYYXxR5C1FVw0qxeJ0+eAPQwl0fDKxS1eXu+1qADTfhwV2D2xMRJjz6qY2TaIypuYLekbmA/wL/FJEncQYBfjtgW3dgWSXHUy0dER/DDQPb8OXCDXy3bLPf4RhjjC+KTGAi0jGvxaC7Xuzi8Xx3A68A7XFuOT4esK0n8EEZv0eNM2JAS5rVj+WhiQs4kGMXrsaYmqe4Z2DzcRpdzHLXi+qILO62Em8jqmo28HAR284u6fMmX0xkOPed1oFr35nDf2as4or+NvqWMaZmKS6BDQQWBqyXm4jkAP0KzuzsbusJzCrF87Qa76SOjRjQpiHPTlnKmV2b0KBOdMkfMsaYEFFkAlPVaYWtl1Nxw0dEAta5qRREhAfO6Mjg56fz9JSlPD60i98hGWNMpfE8nUpZiUgz8lseAnQXkZgCu8UAw4AVwY4n1LRtFMdlfZsz5seVXNKnGZ2a2DiJxpiaobiOzJso+rnXYVQ1qYhNVwAPuMcqrr/YXuAqr+cz+W45sR2fzF3LQxMXMu6avjZOojGmRijuCmwkpUhgxRgFfEj+eIqXcPi4ivuBP1TVJrwqg/jYSG4/OZX7Pp7Pp7+t4/Qjm/gdkjHGBF1xz8AerIgTqOomYBOAO1HlOlUtzXQsxoMLj2rGuzP+4PFPF3FC+0bUirK2MMaY0FaqjswiUk9EjhGRi0WknlsWIyLF9ScLHOtoExAhIrFFLWX7GiY8THjwzE78uWMfL09bXvIHjDGmmvM6I3OEiPwLZ46uacA7QF7Ho49wnnEVZZeI9HbXM4FdJSymjHq3rM/pRzbm5WnLWbNtj9/hGGNMUHlthfgYcDVwI86MzL8HbPsE+AtFJ7ErgeUB6xXxXM0U4d5TO/DVog38Y/JiRl5iU6wZY0KX1wR2Oc6MzG+KSMGHK8uBVkV9UFXHBKy/VeoITak0SajFdce14dmvlnLp8i30a93A75CMMSYovD4DSyD/KqqgKDwMI1VZROQWEVkgIvNFZKz7jK6liMwUkQwRGSciUe6+0e77DHd7i4Dj3OOWLxGRk337QmVw7XGtSE6oxUMTF5Bt4yQaY0KU1wQ2HxhSxLZTgJ+9HEREIkXkdhH5QUT+EJGNBReP8RR1/GTgr0AvVe2Mk1gvBP4JPKuqbYBtwAj3IyOAbW75s+5+uIMTXwh0wpn2ZVQhV55VVt44iYvX7+K9WX/4HY4xxgSF11uIjwIfiUgtnClRFOgmIkOBa4EzPR7nWXf/STjP0oLRnD4CqCUiB4BYYB1wPHCxu30M8CBOh+oh7jo4fdVeFKcX8BDgfbdf2goRyQB6Az8GId6gOKXzEQxo05Anv1jCKZ0bkxhn4yQaY0KLpwSmqp+IyMXAv3AaYgC8BqwFLlPVLzye7zycZ2lPlzpSD1R1rYg8BfyBM7LHl8AcYLs7Ej44LSmT3fVkYLX72WwR2YEza3QyMCPg0IGfqRZEhIeHdGLwc9P5x+RFPHNBN79DMsaYCuUpgYmIqOoHwAci0g5oCGwFlqiqikicqnppAp83GkdQuH3ThuA08d+Oc7U4OIjnuwa4BiAxMZH09PRgnarMTm4ezvhf1tIucgvt61fOXdDMzMwqWRd+sLrIZ3WRz+qiYni9hfgWzmC7qOpSYGneBhFpCHwO9PJwnFeBi4AppYrSuxOBFe7oH4jIeKA/kCAiEe5VWArOlSPua1NgjYhEAPHAloDyPIGfOUhVRwOjAVJTUzUtLS0Y36lc+hydwy/PTuOjleFMHnIMkeHBn4Q7PT2dqlgXfrC6yGd1kc/qomJ4/Wt2kogcNgiviDQGvi3FcTYAx4vIVBG5V0SuL7Bc5/E4RfkD6OuO6iHACThzmk0FznX3GYbTdw1ggvsed/s3qqpu+YVuK8WWQFuciT2rnVpR4Tx0ZieWbczkje9ssH9jTOjwegU2CJgqIvtU9RYAt8n5NzhJ6RSPx3nOfW0GHFfI9uJGqy+Rqs4UkQ9xWkVmA7/gXCF9CrwvIo+6Za+7H3kdeMdtpLEVp+UhqrpARD7ASX7ZwA2qmlPWuPx2QodGnNihEc99tYwzujahSUItv0Myxphy89qIY76InAR8LSJ7cYaSmoJzK/EMVd3t8ThBv3+lqg9w+Kggv+O0Iiy47z6chiWFHecxnBFIQsIDZ3Rk0LPTeGTSQl66tKff4RhjTLl5Tiiq+gvOldYNwGycK5lTvCYv46+m9WO56fi2fDZ/PVOXlKu7nTHGVAnFTWh5fRGbPsW5pfgVMMKdPFFVtcRbfyJyakn7qOrkkvYxZXP1Ma0Y//MaHvhkAf1uaUBMZLXpm22MMYcp7hbiiyV89tmAda/Pria5+xacMjhwgF/7qxokURFhPDKkMxe/NpOX0pdzy6B2fodkjDFlVtyElsF4XtWykLJ6wMnAFcDwIJzTBDi6TUPO7NqEl6YtZ2j3ZFo0rO13SMYYUybB7xQUQFVXFbLMVdV/4ozscW9lxlNT/d9pHYgKD+P+CQtweg0YY0z1U9xMyh1FJDpgvdilAmL5BWfMQhNkSXVjuHVQO75duonP5q/3OxxjjCmT4p6BzQf64nTgnU/RE1GKu63Mz67c6U2G4wy8ayrB5f2aM/6XNTw4YQH92zQkvlak3yEZY0ypFJfABuJ05M1bLzcRmc3hiTAKaAHE4TwHM5UgIjyMfww9kiEjv+Nfny/msaFd/A7JGGNKpbhGHNMKWy/IHQvR6y3EBRyewPbhDLr7P1Vd4PE4pgJ0SYnnyv4tee27FZzVPZmjWtT3OyRjjPHM61BSxTkO+AAPtxBVdXgFnM9UoFsGteOz+eu5Z/xvfPrXAURHWC8GY0z1UKmtEE3VUzs6gkfP6kzGxkxeTv/d73CMMcYzS2CGge2TOKNrE0ZOzSBjY6bf4RhjjCeWwAwA95/ekVpR4dw7/jdyc61vmDGm6rMEZgBIjIvmvlM7MGvlVj74abXf4RhjTImKG8x3E0X3/QoUXXHhGD+d1yuF8b+s4fHJizi+QxJJcTF+h2SMMUUqrhXiSLwlMM9EZA7OJJJjVXVbRR7blJ+I8PjQLgx+fjoPTVzIyIt7+B2SMcYUqbh+YA8G4Xy/Af8EnhaRCTjJbIragHxVRqvEOtw0sA1PT1nKkK7rOanTEX6HZIwxharswXyHA0fgTIp5BPA58IeIPCYibSozFlO0a49rTYfGdbnvf/PZvme/3+EYY0yhKr0Rh6ruVtU3VPU4oC3wJnAJsEREvhWR4SJiD198FBURxpPnHsm23ft5eNLCkj9gjDE+8LsVYg75z9lycAYGHgWsFJFBvkVl6Jwcz/VprRn/81q+WbzB73CMMeYwlZ7ARCRWRIaJyFQgA7gAJ2k1VdVjgBTgG+CVyo7NHOrG49uS2iiOe8b/xo69B/wOxxhjDlGpCUxE3gTW47RwXAUMVNX2qvovVd0AoKpbgedxRqg3PoqKCOPJ845kc+Z+HrVbicaYKqayr8A6ArcDjVV1uKpOL2K/BVTQFC6mfI5MSeDaY1vx3zlrSF+y0e9wjDHmoOI6Mq+gFP3AVLVVcdvd2Z0nAjNVdVcJx8oEipzCxVSum09sy5SFG7hn/G98ccux1I2xyS+NMf4r7grsowJLJBCPM0PzJPc1HicJfljSiVQ1C7gHSChXxKbSRUeE8+R5Xdmwcx//mLzI73CMMQYoviPz7XnrInIvsBw4TVV3B5TXwUlmOz2ebxbQA7u6qna6NU3g6mNb8cq03zmtSxMGtG3od0jGmBrO6zOwG4AnA5MXHLzV95S73Ys7getF5EYRaSUitd1WiQcX76GbynbLie1olVibOz+cZ60SjTG+85rA6gKNith2BFDH43FmAq2BF4BlOFduuwos5SIiCSLyoYgsFpFFItJPROqLyBQRWea+1nP3FRF5QUQyRORXEekRcJxh7v7LRGRYeeMKBTGR4Txzfjc27MrioQkL/A7HGFPDFTeYb6CJwJMishOYoKr7RSQKGIIztuFEj8e5kgoeILgQzwOfq+q5boyxwL3A16r6hIjcDdwN3AWcgjMaSFugD/AS0EdE6gMPAL3ceOeIyAQbgNi5lXjjwDY8//UyTuzYiFO7NPY7JGNMDeU1gV0HvAV8AKiI7ALicEbOmOBuL5GqvlX6EL0TkXjgWGC4e779wH4RGQKkubuNAdJxEtgQ4G13MOEZ7tVbY3ffKW6fNERkCjAYGBvM+KuLG49vw9QlG7n349/o1bweSXVt5C9jTOXzlMBUdQcwVEQ6AUfh3E5cD8xW1arUw7UlsAl4U0S6AnOAm4FGqrrO3Wc9+bdDk4HA2RvXuGVFlR9CRK4BrgFITEwkPT29wr5IVXdRi1we+PMAV41O55ae0YjIwW2ZmZk1qi6KY3WRz+oin9VFxfB6BQaAqi7A6WRcZiJyAXA10A447J/uqppUjsNH4LRyvElVZ4rI8zi3CwOPryJSIbcxVXU0MBogNTVV09LSKuKw1cb++it5YMIC1tZqxaV9mx8sT09Pp6bVRVGsLvJZXeSzuqgYnkfiEJEkEfmniHwtIkvcqzFE5GYR6efxGBfj3MLLwBnzcAJOM/wwnAYdL5b2CxSwBlijqjPd9x/iJLQN7q1B3Ne8ISXWAk0DPp/ilhVVbgJc1rc5x7RtyGOfLmLF5t0lf8BUvkcecRZjQpCnBCYivXFaDZ4DrATaANHu5sbAbR7PdwfwCPnN7kep6pU4t/42A3s8HqdQqroeWC0iqW7RCcBCnESZ15JwGPCJuz4BuNxtjdgX2OHeavwCOElE6rktFk9yy0yAsDDhyXO7EhURxi3j5pKdk+t3SKagv//dWYwJQV6vwJ4FpuLc9rsWp/FGnllAb4/HaQt8r6o5ONOn1AVwh5b6J3Cjx+MU5ybgXRH5FegGPA48AQwSkWXAie57gMnA7zhXhK8C17vxbMVJtLPd5eG8Bh3mUEfEx/DoWZ2Zu3o7L3y9zO9wjDE1iNdnYD2AIaqaK4FP6x1bAK/PrXaSf+W2FuiA0yIQnKTYwONxiqSqc3Gavxd0QiH7KkV0wlbVN4A3yhtPTXBG1yZMW7qJf0/NoF9rG6GjSjnlFOf1s8/8jcOYIPB6BbYDSCxiWyvA64yHs4Ej3fUJwP0icrXbUfhJYIbH45gq5qEzO9GiQW1uGTeXzP3B7upnPNu711mMCUFeE9gE4CERCRxxXkWkIc70KOM9HucfwB/u+v04tx9fAt7EeQZ2rcfjmCqmdnQE/76oO1t2Z/HG/Cyci1tjjAkerwnsLpzbfwuBb92yl4ElwF6cZFQiVZ2hquPc9e2qOgSoDSSoah9V/b00wZuqpXNyPHcNbs/PG3P4z8w/Sv6AMcaUg6cE5g6h1BfnedEq4CtgBU4fq/4lze8FICIxIpIlImcVOHaWqnodzd5UcVf2b0mXhuE8OmkhS9aXe2hLY4wpkueOzO6wTK+7S6mp6j4R2Qhkl+XzpnoICxOu6hLNI7NzuGnsz0y4cQAxkeF+h1VznX663xEYEzRe+4F9IyIviUitQrZ1EJFvPJ7vFeCvImJT+oaw+GjhmfO7snRDJg9NtFHrfXX77c5iTAjyegWWhnMLsZ+IDFXVFQHb6gLHeTxOAtAZWCkiX+O0Xgx82q+qepfHY5kq7Nh2iVyX1pqX0pdzVIv6nN0jxe+QjDEhpjRjIV6EM5LGHBG5VFUnl+F85wBZ7voxhWxXnAYjJgTcNqgdP6/axn0fz6dzcjztGsX5HVLNkzfeng0ca0KQ57EQgXU4V2LjgAki4qnlYSBVbVnC0qrko5jqIiI8jH9f1J3a0RFc95857M6yx5/GmIpTmgSGqmar6nU4/bXuEZFPgHivnxeRy0Wk0NE23FmTLy9NPKbqS6obwwsXdWPF5t3cM/436x9mjKkwpUpgeVT1dZyrsR44V2RevQm0LmJbS3e7CTFHt27IrYPaMWHen9Y/zBhTYbwmsDE4E0Ue5E5Z0gtneCivf5UKjqMYqAFOZ2kTgq5Pa0NaaiKPTFzIr2u2+x2OMSYEeJ2R+YoiyjfgTDVSJBEZAgwJKPq7iGwqsFsMTqOO2V7iMdVPWJjw7PndOO2F6Vz3n5+ZeNMA6teO8jus0Hf++X5HYEzQFJnARCRWVffkrZd0oLx9C5EEdAl43xo4osA++4EvgUdLOo+pvurVjuKlS3ty3is/cuN7P/P2lb2JCC/TXWzj1fXX+x2BMUFT3BXYLhHpp6qzgEwO7a9VmEKHW1DVV3Hm2kJEpgLXq+qisgRrqr+uTRN47KzO3PHhr/zjs8X8/fSOfocU2va4/66MLfHfoMZUO8UlsCuB5QHr5W4+pqoDy3sMU/2d16spC/7cyevfraBzcl2GdrdOzkFz6qnOq/UDMyGoyASmqmMC1t8q6wlE5Hrgv6q6yV0vjqrqS2U9l6k+7jutA4vW7eTuj36jTWIcXVI898YwxhigdCNxlNWLwE84rRhfLGFfxZkfzIS4yPAwRl7SgzP//R3XvvMTE24aQMM60SV/0BhjXMU14phNKW4bqmrvIsrDCls3pmGdaF65rBfnvvwDN7z7M/+5qg+R1qjDGONRcVdgC6iA517GFKdLSjxPnNOFW8bN4/5P5vP40C6IFNdd0BhjHMU9AxsejBOKSBQwHOgNNMYZY3EmMMadc8zUMEO7p5CxMZORU5fTqmEdrj7WhsSsMMOH+x2BMUFTGc/ADhKRDsDnQBNgDrARZ3qVy3E6OA9W1YWVGZOpGm4blMqKzbt5/LNFNGsQy8mdCnYVNGViCcyEMM8JTERaAJcC7XBGzjiEqnrp8j8a2AEco6oHh58SkWbAJOBl4FivMZnQERYmPH1eN9Zu+5G/vT+X//6lH52TrWViuW3e7Lw2bOhvHMYEgdcZmXviPBO7xF3a4oyDeC7ORJde/+/oBdwfmLwA3PcPAEd5PI4JQbWiwnl1WC/q145ixJjZrN+xz++Qqr9zz3UWY0KQ1yZfTwL/xbndJ8AId+6uATgNPf7l8TgrKeTqzRWD90GBTYhKiovhtWG9yNyXzYgxs20OMWNMkbwmsG7AWCDXfR8DoKo/AA8BT3g8zt3AoyLSJ7BQRPoCj2CzMRugQ+O6vHhxDxat28lNY38hOye35A8ZY2ocrwlMgf3qzEa4EWgesG01zi3FQonIbBGZJSKzgPuAusAPIrJOROaJyDrge7f83rJ8iULOGS4iv4jIJPd9SxGZKSIZIjLObQmJiES77zPc7S0CjnGPW75ERE6uiLiMdwPbJ/HwkM58s3gj935sE2EaYw7ntRHHQpxR5KcCPwK3iMhPOKPI30n+mImFKdifbEEZ4iytm4FFOEkR4J/As6r6voi8DIzAGfFjBLBNVduIyIXufheISEfgQqATTovJr0SknarmVELsxnVp3+Zs3JXFC18vIzEumjtObu93SMaYKsRrAhtN/lXXvThTnyx23+/GacxRqGD1JyuKiKQApwGPAbeK0yv2eOBid5cxwIM4CWyIuw7wIfCiu/8Q4H1VzQJWiEgGTr+1HyvpaxjXLSe2ZdOuLEZOXU5SXAzDjm7hd0jVy3XX+R2BMUHjdULLdwLWF7n9uY7GeRY2Q1U3lnQMEYnBaUJ/gar+r2zhevIczlVhnPu+AbBdVfNaA6wBkt31ZJxboKhqtojscPdPBmYEHDPwMweJyDXANQCJiYmk24jfAGRmZlZoXZyYoCxOCufBCQvY8EcGvY+o1O6L5VLRdVFqjRo5r1Xgt+l7XVQhVhcVo0x/CVQ1E+cqrDSf2SciG4GgNSsTkdOBjao6R0TSgnWePKo6GufqlNTUVE1LC/opq4X09HQqui76H5PDpa/N5LXfdjDgqG4c3bp69GsKRl2UyurVzmvTpv7F4PK9LqoQq4uKUZqOzDE4nYxTOLwpvNdpUF4B/ioiX6jqAe9hetYfOFNETsWJsS7wPJAgIhHuVVgKsNbdfy3QFFgjIhFAPLAloDxP4GeMD2Iiw3ltWC/Oe/lHrnl7Du9e1YeuTRP8Dqvqu+wy59X+tW9CkKcEJiIDgI+AxCJ28ToNSgJOX7KVIvI1sIFDG3ioqpa5Kb2q3gPc48acBtyuqpeIyH9xntO9DwwDPnE/MsF9/6O7/RtVVRGZALwnIs/gNOJoC8wqa1ymYiTERvHOiD6c/8qPXP7GLN6/pi8dGtct+YPGmJDktRn9C8DvQHcgWlXDCizhHo9zDpCF03rxGJykcV6BJRjuwmnQkYHzjOt1t/x1oIFbfitOPzVUdQHwAU7ry8+BG6wFYtVwRHwM717Vh9iocC59bSYZGzP9DskY4xOvtxBTgbNVdV55TqaqLcvz+VKeKx1Id9d/x2lFWHCffRSRNFX1MZyWjKaKaVo/lnev6sP5r8zgktdm8MG1/WjeoLbfYRljKpnXK7BfARse3FQZrRLr8O5VfdifncvFr87kz+17/Q7JGFPJvF6BXQe8JSIrVXVaeU/qPlMralT7UeU9vqkZUo+I450Rfbjo1Rlc/OoMxl3bj0Z1ixpqs4a67Ta/IzAmaLwmsClALPCNiOwHdhXcQVWTSjqIiDQCvgY64jTeyJt6N7AhhyUw41nn5HjeuqI3l78+kwte+ZH3ru5Lk4RafodVdZxxht8RGBM0XhPYSA5NMmX1NE5n5qY4HYj74LREvBRnUsvTKuAcpobp2bweb4/ow/A3ZnHB6B9576q+NK0f63dYVcOSJc5raqq/cRgTBF5H4niwgs53HM44hevc9+LOBfa4iIThXH3ZwLmm1Ho2r8d/rurDZa/P5MLRMxh7dV+aNbAkxrXXOq/WD8yEIK+NOCpKArBJVXOBnUDgbccfcIanMqZMujZN4L2r+7J7fzbnv/IjKzbv9jskY0wQFZnAROQDEWkdsF7s4vF8K4DG7nreDM95zgC2luVLGJOnc3I8Y6/uy4GcXM5/5UeWbTjsca0xJkQUdwWWCES660nu++IWLz4FTnLXHwXOEZE1IrIC+Cvw71JFb0whOjSuy/vX9AXg/Fd+ZO7q7f4GZIwJiiKfganqwID1tIo4mTvUU976ZyJyNDAUqAVMUdXPKuI8xrRtFMeHf+nHZa/P4uJXZ/DypT05tp3Xf2cZY6oDX+elUNWfgJ/8jMGEruYNavPhdf0Y9sZsrnxrNk+f35Uh3Q6bFSe0/d//+R2BMUHjdTDf+4vZnNcgY15hnZxFpKWqrihNUG6LxGRVXV2azxlTUFJcDOOu7cvVY37i5vfnsnX3fq7oX2kjmvnvxBP9jsCYoPF6BXYTzqgZeQPOZQJ13PXd7nGiRWQucIqqbgj47DwRmYgzcO5UVS2yP5k7m/IFwA04o9s/6TE+Y4pUNyaSMVf25q9jf+GhiQvZnJnF7Sel4ky+HeLmznVeu3XzMwpjgsJrAjsVeBe4D5igqlkiEg0MwWmMcQXOqBpjcTorXxrw2Q7A/wETgT0iMguYD2zGGZk+AWgJ9MQZoeM34E5V/bBc38yYADGR4Yy6pAd//2Q+I6cuZ/XWvfzr3COJifQ6kUI19be/Oa/WD8yEIK8J7EXgCVX9b16BqmYBH4hIHPBvVe0hIo/iJDQC9lsLXCcid+JcXR0PnIUzOHAMTtP5JTgJ7gpVnVO+r2RM4SLCw3h8aBea1o/lX58vYe32vYy+rCcN6kT7HZoxpgy8JrAjgfVFbFuHc5UFsBiIK2wnVd0FvOYuxvhCRLg+rQ0tGtTmlnFzOWvU97wx7CjaNir0Z2uMqcK8jsSxFLhZRKICC93biLfgXEGBc1W1AWOquFO7NGbctf3Yuz+Xs1/6ge+WbfY7JGNMKXlNYDcDA4A1IvKuiDwnIu/iDMh7tLsdnBmbx1d8mMZUvG5NE/jkxv4kJ9Ri2JuzeOv7FRTTxsgYU8V4Hcw3XUTa4lxt9QJ64NxSfAt4TlX/dPe7O0hxGhMUyQm1+O9f+nHLuHk8OHEhv67ZwWNDu1ArKkQadzz+uN8RGBM0JSYw9zbh7cAkVb0j+CEZU7niYiIZfVlPRk7N4JmvlrJ4/S5euaxnaEzJcrSNj21CV4m3EN3WhvfhNHc3JiSFhQk3ndCWN4YdxZptezj9398xbekmv8Mqvx9+cBZjQpDXZ2AzcW4bGhPSBrZPYuJNA2gcH8PwN2fxwtfLyMmtxs/F7r3XWYwJQV4T2J3A9SJyo4i0EpHaIhIbuAQzSGMqU/MGtfn4+v6c1S2ZZ6Ys5dLXZrJh5z6/wzLGFFCaK7DWwAvAMpyxD3cVWIwJGbWiwnnm/K48ee6RzF29nVOen87UxRv9DssYE8BrR+YrgWp8H8WY0hMRzuvVlO7N6nHjez9zxVuzuWpAS+4c3J6oiMqezNwYU5DXZvRvBTkOY6qsNkl1+N8N/fnH5EW89t0KZqzYwrPnd7PRO4zxmf0z0hgPYiLDeWhIZ165rCdrt+3ltH9/x+hvl1f9Bh7PPecsxoQgzwlMRC4Qka9E5A8R2VhwCWaQXolIUxGZKiILRWSBiNzsltcXkSkissx9reeWi4i8ICIZIvKriPQIONYwd/9lIjLMr+9kqpaTOx3Bl7ccR1q7RB6fvJgLXvmRlZt3+x1W0bp1s6lUTMjylMBE5GJgDJABpAATgEnu53fijFZfFWQDt6lqR6AvcIOIdATuBr5W1bbA1+57gFOAtu5yDc4cZIhIfeABoA/QG3ggL+kZkxgXzSuX9eSZ87uyZMMuTnl+OmN+WEluVbwa++orZzEmBHm9ArsDeARnokmAUap6Jc48XpuBPUGIrdRUdZ2q/uyu7wIWAck485aNcXcbgzOdC2752+qYASSISGPgZGCKqm5V1W3AFGBw5X0TU9WJCGf3SGHKLcfRu2V9HpiwgAtG/8iyDVWsQe6jjzqLMSHIawJrC3yvqjlADlAXDiaJfwI3Bie8shORFjiDC88EGqnqOnfTeqCRu56MMyBxnjVuWVHlxhziiPgY3rriKP517pEs25jJqS9M56kvlrDvQI7foRkT8rw2o98J5M36txZn/q90970ADSo2rPIRkTrAR8DfVHVn4NTxqqoiUiH3ekTkGpxbjyQmJpJus94CkJmZWePqIgl4uG8k7y9WXpyawQczlzOsYzTNY/b6Whfdtm8HYG4V+O9RE38XRbG6qBheE9hsnEktv8B5/nW/iGQD+4H7gRnBCa/0RCQSJ3m9q6p5U7tsEJHGqrrOvUWY1+hkLdA04OMpbtlaIK1AeXrBc6nqaGA0QGpqqqalpRXcpUZKT0+nptbFmSfB9xmbue/j33jypz30bRzBs8P70Di+lj8BJSQAVIn/HjX5d1GQ1UXF8HoL8R/AH+76/cAsnAYPb+I8A7u24kMrPXEutV4HFqnqMwGbJgB5LQmHAZ8ElF/utkbsC+xwbzV+AZwkIvXcxhsnuWXGlKh/m4Z8/rdj+evxbfhpQw7HPzWNF75eZrcVjalgXjsyz8C9ylLV7cAQd5qVaFXdGbzwSq0/cBnwm4jMdcvuBZ4APhCREcAq4Hx322TgVJzWlXuAKwBUdauIPIJz5QnwsKpurZRvYEJCTGQ4t56UStPstUzdFs8zU5YybvZq7jm1Pad1aUzgbe2geuWVyjmPMT7wegvxMO40K1kVGEu5qep3OM/kCnNCIfsr+S0rC257A3ij4qIzNVFibBijTu3JjN+38NDEhdz43i+MabGSu0/pQM/mldAzIzU1+OcwxidFJjARub8Ux1FVfaQC4jEmJPVt1YBJNw3gg59W8/SXSznnpR84sUMj7jg5ldQjgjgk1cSJzusZZwTvHMb4pLgrsAeBvcBuir6qyaM4/cSMMUUIDxMu6t2MId2a8Ob3K3k5fTmDn/+Ws7olc8uJ7WjWIAizEj39tPNqCcyEoOIS2HKgOTAHeB8Y7/b7MsaUQ2xUBDcMbMMlfZrx0rTlvPX9Sib9+ifn9mzK9WmtaVrfptczxosiWyG6wy4dDSzAubraICLjReQ8EfGpTbAxoSMhNop7TunAtDsGcn6vpnw0Zw1pT6Vz2wfzWL4p0+/wjKnyim1Gr6o/qertqtoMZyil9TjjHm4UkXdF5NjKCNKYUHZEfAyPDe3CtDvTuLxfcyb9+icnPjONm8b+wqJ1VamRrzFVi+fR6FX1W1W9Hqfj78vABcDfghSXMTVO4/haPHBGJ76763iuObYV3yzawCnPT+ey12eSvmQjTqNZY0wez83oRaQ/cCFwLhAHfIg7ersxpuIkxkVzzykduO641rw78w/G/LCS4W/Opm1SHa4c0JKh3ZOJiQz3drB33glusMb4qNgE5s6PdSHO1VYj4HPgFmCCqlaJEeiNCVUJsVHcMLANVx/Tikm//slr01dwz/jfePKLJZzXK4WLezejeYPaxR+kadPitxtTjRXXD2wJznQp3+DMjTW+io26YUyNEBURxtk9UhjaPZkZv2/lze9X8Nr0Fbwy7XeOaduQS/o058QOSUSEF/JEYNw45/WCCyo3aGMqQXFXYG2BfUBPoAfwr+KGv1HVpIoNzRgTSETo17oB/Vo3YN2OvYybvZr3Z63mL/+ZQ6O60ZzTI4Wze6TQJqlO/odecu/yWwIzIai4BPZQpUVhjCmVxvG1+NuJ7bhxYBumLtnEuzNX8fK05YxKX07XlHjO7pHCGV2bUN/vQI0JoiITmKpaAjOmiosID2NQx0YM6tiIjTv3MWHen3z081oemLCARyYt5NMNu2hYJ5raB3K8N/wwppoo82C+xpiqJaluDFcd04qrjmnFonU7+fiXtWSOzWbb7v1c+cgUBqYmcUqXIxiYmkTtaPtf31R/9is2JgR1aFyXDo3ros3qsWPfAc7qnsyXC9bz6W/riI4I49h2iZzS+QhO6NCI+FqRfodrTJlYAjMmhMlHH5IAPN6wIY8M6cxPK7fy2fz1fLFgPVMWbiAiTOjZvB4D2ycxMDWJdo3qVN5cZcaUkyUwY0JZw4YHV8PDhD6tGtCnVQPuP70j89Zs56tFG5i6eBNPfLaYJz5bTOP4GNJSkxiYmsjRbRpSx241mirMfp3GhLK33nJehw8/pDgsTOjerB7dm9XjjpPbs37HPqYt3cjUxZuYOO9Pxs76g/AwoUtyPP1aN+Do1g3o1bw+taKsIYipOiyBGRPKikhgBR0RH8MFRzXjgqOasT87lzmrtvHD8s38sHwLr377Oy+lLycyXOjWNIF+rRrQq0V9ujVLoG6MPT8z/rEEZow5RFRE2MEO07cBu7Oymb1yKz/+voUZy7fw4tQMchVEoF1SHD2a16Nn83r0aJZAy4a17RmaqTSWwIwxxaodHUFaahJpqc5gO7v2HWDe6h3MWbWNOX9sY9Kvzi1HgITYSLokx9OpSTydk+vSuUk8zerHEhZmSc1UPEtgxphSiYuJZEDbhgxo6zQQyc1VMjZl8vOqbfzyx3YWrNvB69/9zoEcZ/qXuOgIOjapS0JuFlvrrqFdozhaJ9ax52mm3CyBGWPKJSxMaNcojnaN4riwdzMA9mfnsnTDLhb8uYP5a3cy/88dTF2TzRer5gHO7cem9WJp16gObRvF0TapjiU2U2qWwIwJZZMn+3LaqIgwOifH0zk5nguOcsq+/mYqzTsfxbINu1i6IZOlG3exbMMupi3ddPBqDaBxfAzNG8TSokFtmjeoTYsGsTRvUJvmDWJtBBFzCPs1GBPKYmP9juCg8DChTVId2iTV4ZQu+eUHcnJZtWU3SzdkkrExk5VbdrNqyx6+WrSBzZn7DzlGYlw0LRrEkpxQiyYHlxiaJNSicXwt6sZEWCOSGsQSmDGhbNQo5/X66/2NoxiR4WG0SYqjTVLcYdsys7JZ5Sa0lVt2s2qz8+o0HllHdq4esn+d6AiaJMTQOL4WjeNjSIyLJjEumiT3NbGOU2a3KUODJTBjQtkHHzivVTiBFadOdASdmjitGgvKyVU2Z2axdvte1m3fx5/b9zrrO/by5/Z9LPhzJ1t3Z1Egxx08rpPQnMRWv3YU9WIjSYiNol5t9zU2ioRakdSLjSIuJsJaUlZBlsCMMdVSeJjQqG4MjerGQLPC98nJVbbu3s/GXfvYtCvLWTKz8td3ZbFo3U627tnPjr0H0EKSXd654mtFkhDrJLT4WpHUiY6gTkwEcTERxEVHUCc6griYyICy/PU60RFER4TZ7c0KZgmsCCIyGHgeCAdeU9UnfA7JGFNK4WFy8DZiSXJylZ17D7Btz3627TnA9j372b7HeV/wddOuLH7flElmVjY792WzPzu3xONHhAm1osKpFRmO5Oyn3txviY0KJzYqgpjIcGLdbbWinCU2YD0mIpyoiLCDS7S7RIWHH3xfcHtUeOgnTEtghRCRcGAkMAhYA8wWkQmqutDfyIwxwRIeJtSrHUW92lGl/mxWdg67s3LYte8Au/Zls2tfNplZ2WRmHfp+7/4c9u7PYdXaP6mTEMu+Azns2Z/Nlt372bs/m70HctizP4d9B3IOaZlZVlHhAYktPIyIcCEyPIyIMCE8zF0PFyLChIiwgPXwMCLDhfCwMCLDhIi89fDD9wsXIUyc7hRhIoSHQZiIuzj1KuKcL0zytznl5ft+lsAK1xvIUNXfAUTkfWAIYAnMGHOY6IhwoiPCqe8x+aWnbyUtrVex+xzIyWXvASfh7TuQw/7sXLKyc9mfk0vWAed1f3auW+5szyvLyts3b8nJIetALjm5yoFcJSc3lwM5SnZOLtm5SnaOkp2by77svPWAbbm5ZOcoB3Kcz2XnKAdy3WNVQJItD0tghUsGVge8XwP0KbiTiFwDXOO+zRKR+ZUQW3XQENjsdxBVRNWoi6pxK6lq1EXVYHWRL7WsH7QEVg6qOhoYDSAiP6lq8f+kqiGsLvJZXeSzushndZFPRH4q62fDKjKQELIWaBrwPsUtM8YYU0VYAivcbKCtiLQUkSjgQmCCzzEZY4wJYLcQC6Gq2SJyI/AFTjP6N1R1QQkfGx38yKoNq4t8Vhf5rC7yWV3kK3NdiBbVc88YY4ypwuwWojHGmGrJEpgxxphqyRJYKYnIYBFZIiIZInJ3IdujRWScu32miLTwIcxK4aEuhovIJhGZ6y5X+RFnsInIGyKysah+gOJ4wa2nX0WkR2XHWFk81EWaiOwI+E3cX9kxVhYRaSoiU0VkoYgsEJGbC9mnRvw2PNZF6X8bqmqLxwWnQcdyoBUQBcwDOhbY53rgZXf9QmCc33H7WBfDgRf9jrUS6uJYoAcwv4jtpwKfAQL0BWb6HbOPdZEGTPI7zkqqi8ZAD3c9DlhayP8jNeK34bEuSv3bsCuw0jk4xJSq7gfyhpgKNAQY465/CJwgoTmippe6qBFU9VtgazG7DAHeVscMIEFEGldOdJXLQ13UGKq6TlV/dtd3AYtwRvkJVCN+Gx7rotQsgZVOYUNMFfyPcHAfVc0GdgANKiW6yuWlLgDOcW+NfCgiTQvZXhN4rauaop+IzBORz0Skk9/BVAb3UUJ3YGaBTTXut1FMXUApfxuWwEwwTQRaqOqRwBTyr0xNzfUz0FxVuwL/Bv7nbzjBJyJ1gI+Av6nqTr/j8VMJdVHq34YlsNLxMsTUwX1EJAKIB7ZUSnSVq8S6UNUtqprlvn0N6FlJsVU1NjSZS1V3qmqmuz4ZiBSRhj6HFTQiEonzB/tdVR1fyC415rdRUl2U5bdhCax0vAwxNQEY5q6fC3yj7hPKEFNiXRS4l38mzn3vmmgCcLnb4qwvsENV1/kdlB9E5Ii8Z8Ii0hvnb1Ao/gMP93u+DixS1WeK2K1G/Da81EVZfhs2lFQpaBFDTInIw8BPqjoB5z/SOyKSgfMw+0L/Ig4ej3XxVxE5E8jGqYvhvgUcRCIyFqcFVUMRWQM8AEQCqOrLwGSc1mYZwB7gCn8iDT4PdXEucJ2IZAN7gQtD9B94AP2By4DfRGSuW3Yv0Axq3G/DS12U+rdhQ0kZY4ypluwWojHGmGrJEpgxxphqyRKYMcaYaskSmDHGmGrJEpgxxphqyRKYMcaYaskSmDHGmGrJEpgxIUZEbhIRFZErC9kWLyK5IvKNH7EZU5EsgRkTevLGnJxTyLYeOHNPFbbNmGrFEpgxoacHsA9YUMi2vOT2c+WFY0xwWAIzJoSISAzQAfjVnY+uoOKuzoypViyBGRNajsQZpLuoBNUT2AUsq7SIjAkSS2DGhJYe7uthCUxE4oE2wC8hPAK8qUEsgRkTWopMYDhTWgj2/MuECEtgxoSWvARW2Ky+eXPT2fMvExIsgRkTItwp27u4bwcU2HY+cKn79rfKjMuYYLEJLY0JESLSHef24AYgHvgYZ0r2rkB3YDfQCBgPPKmqM3wK1ZgKEeF3AMaYCpN3+/BhoBUwHKgFzAQGAse525oDf/gQnzEVyq7AjAkRIjISuB7oo6qz/I7HmGCzZ2DGhI4eQA72jMvUEHYFZkwIEJFwYCewQlU7+x2PMZXBrsCMCQ3tgVhgrs9xGFNp7ArMGGNMtWRXYMYYY6olS2DGGGOqJUtgxhhjqiVLYMYYY6olS2DGGGOqJUtgxhhjqiVLYMYYY6ql/wdnCr9K6FFqgAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# now plot it again but zoom in a little bit\n", "ax34 = plt.subplot(111)\n", "ax34.plot(muScan, lik)\n", "ax34.set_xlabel('$\\mu$', size=20)\n", "ax34.set_ylabel('Marginalized Likelihood\\n (arbitrary units)', size=15)\n", "ax34.set_ylim(bottom=0)\n", "ax34.set_xlim(muScan[0],2.5)\n", "ax34.plot([mu95, mu95], [0, 0.5*np.amax(lik)], color='red', linestyle='dashed',\n", " label='95% limit: $\\mu$ ='+str(round(mu95,2)))\n", "ax34.grid()\n", "_ = ax34.legend()" ] }, { "cell_type": "code", "execution_count": null, "id": "9052e781", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.7" } }, "nbformat": 4, "nbformat_minor": 5 }