cls.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to JupyROOT 6.10/02\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from IPython.display import HTML, IFrame, Image, SVG, Latex\n",
    "import ROOT\n",
    "from ROOT import RooFit, RooStats\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(ws) importing dataset data_toy_signal\n",
      "[#1] INFO:ObjectHandling -- RooWorkspace::import(ws) importing dataset data_toy_nosignal\n"
     ]
    }
   ],
   "source": [
    "from create_example_ws import create_example_onoff as create_example\n",
    "\n",
    "ws = create_example()\n",
    "\n",
    "# create references to useful variables\n",
    "model_config = ws.obj('model_config')\n",
    "pdf = model_config.GetPdf()\n",
    "obs = model_config.GetObservables()\n",
    "poi = model_config.GetParametersOfInterest().first()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/svg+xml": [
       "<svg height=\"331pt\" viewBox=\"0.00 0.00 539.79 330.96\" width=\"540pt\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\n",
       "<g class=\"graph\" id=\"graph0\" transform=\"scale(1 1) rotate(0) translate(4 326.9605)\">\n",
       "<title>pdf</title>\n",
       "<polygon fill=\"#ffffff\" points=\"-4,4 -4,-326.9605 535.7939,-326.9605 535.7939,4 -4,4\" stroke=\"transparent\"/>\n",
       "<!-- pdf -->\n",
       "<g class=\"node\" id=\"node1\">\n",
       "<title>pdf</title>\n",
       "<ellipse cx=\"267.397\" cy=\"-296.0904\" fill=\"none\" rx=\"58.8803\" ry=\"26.7407\" stroke=\"#ff0000\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"267.397\" y=\"-299.8904\">RooProdPdf</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"267.397\" y=\"-284.8904\">pdf</text>\n",
       "</g>\n",
       "<!-- pdf_sr -->\n",
       "<g class=\"node\" id=\"node2\">\n",
       "<title>pdf_sr</title>\n",
       "<ellipse cx=\"198.397\" cy=\"-206.3503\" fill=\"none\" rx=\"57.0522\" ry=\"26.7407\" stroke=\"#ff0000\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-210.1503\">RooPoisson</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-195.1503\">pdf_sr</text>\n",
       "</g>\n",
       "<!-- pdf&#45;&gt;pdf_sr -->\n",
       "<g class=\"edge\" id=\"edge1\">\n",
       "<title>pdf-&gt;pdf_sr</title>\n",
       "<path d=\"M247.8089,-270.6145C240.4523,-261.0467 231.9756,-250.0221 224.1434,-239.8356\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"226.9041,-237.684 218.034,-231.8899 221.3548,-241.9508 226.9041,-237.684\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- pdf_cr -->\n",
       "<g class=\"node\" id=\"node7\">\n",
       "<title>pdf_cr</title>\n",
       "<ellipse cx=\"337.397\" cy=\"-206.3503\" fill=\"none\" rx=\"57.0522\" ry=\"26.7407\" stroke=\"#ff0000\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"337.397\" y=\"-210.1503\">RooPoisson</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"337.397\" y=\"-195.1503\">pdf_cr</text>\n",
       "</g>\n",
       "<!-- pdf&#45;&gt;pdf_cr -->\n",
       "<g class=\"edge\" id=\"edge2\">\n",
       "<title>pdf-&gt;pdf_cr</title>\n",
       "<path d=\"M287.2689,-270.6145C294.7321,-261.0467 303.3317,-250.0221 311.2774,-239.8356\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"314.0846,-241.9275 317.4753,-231.8899 308.5651,-237.6221 314.0846,-241.9275\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nobs_sr -->\n",
       "<g class=\"node\" id=\"node3\">\n",
       "<title>nobs_sr</title>\n",
       "<ellipse cx=\"59.397\" cy=\"-116.6102\" fill=\"none\" rx=\"59.2941\" ry=\"26.7407\" stroke=\"#0000ff\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"59.397\" y=\"-120.4102\">RooRealVar</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"59.397\" y=\"-105.4102\">nobs_sr</text>\n",
       "</g>\n",
       "<!-- pdf_sr&#45;&gt;nobs_sr -->\n",
       "<g class=\"edge\" id=\"edge3\">\n",
       "<title>pdf_sr-&gt;nobs_sr</title>\n",
       "<path d=\"M164.3937,-184.3974C145.6305,-172.2836 122.1489,-157.1236 102.1526,-144.2137\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"103.8306,-141.131 93.5309,-138.6475 100.0338,-147.0119 103.8306,-141.131\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nexp_sr -->\n",
       "<g class=\"node\" id=\"node4\">\n",
       "<title>nexp_sr</title>\n",
       "<ellipse cx=\"198.397\" cy=\"-116.6102\" fill=\"none\" rx=\"61.5366\" ry=\"26.7407\" stroke=\"#ff0000\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-120.4102\">RooAddition</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-105.4102\">nexp_sr</text>\n",
       "</g>\n",
       "<!-- pdf_sr&#45;&gt;nexp_sr -->\n",
       "<g class=\"edge\" id=\"edge4\">\n",
       "<title>pdf_sr-&gt;nexp_sr</title>\n",
       "<path d=\"M198.397,-179.4253C198.397,-171.3504 198.397,-162.3471 198.397,-153.7362\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"201.8971,-153.5161 198.397,-143.5162 194.8971,-153.5162 201.8971,-153.5161\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- ns -->\n",
       "<g class=\"node\" id=\"node5\">\n",
       "<title>ns</title>\n",
       "<ellipse cx=\"198.397\" cy=\"-26.8701\" fill=\"none\" rx=\"59.2941\" ry=\"26.7407\" stroke=\"#0000ff\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-30.6701\">RooRealVar</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"198.397\" y=\"-15.6701\">ns</text>\n",
       "</g>\n",
       "<!-- nexp_sr&#45;&gt;ns -->\n",
       "<g class=\"edge\" id=\"edge7\">\n",
       "<title>nexp_sr-&gt;ns</title>\n",
       "<path d=\"M198.397,-89.6852C198.397,-81.6102 198.397,-72.607 198.397,-63.9961\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"201.8971,-63.776 198.397,-53.776 194.8971,-63.7761 201.8971,-63.776\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nb_sr -->\n",
       "<g class=\"node\" id=\"node6\">\n",
       "<title>nb_sr</title>\n",
       "<ellipse cx=\"335.397\" cy=\"-26.8701\" fill=\"none\" rx=\"59.2941\" ry=\"26.7407\" stroke=\"#0000ff\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"335.397\" y=\"-30.6701\">RooRealVar</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"335.397\" y=\"-15.6701\">nb_sr</text>\n",
       "</g>\n",
       "<!-- nexp_sr&#45;&gt;nb_sr -->\n",
       "<g class=\"edge\" id=\"edge8\">\n",
       "<title>nexp_sr-&gt;nb_sr</title>\n",
       "<path d=\"M232.6145,-94.1964C250.9787,-82.1672 273.772,-67.2367 293.2244,-54.4947\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"295.1683,-57.4054 301.6156,-48.9981 291.3326,-51.5498 295.1683,-57.4054\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nobs_cr -->\n",
       "<g class=\"node\" id=\"node8\">\n",
       "<title>nobs_cr</title>\n",
       "<ellipse cx=\"337.397\" cy=\"-116.6102\" fill=\"none\" rx=\"59.2941\" ry=\"26.7407\" stroke=\"#0000ff\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"337.397\" y=\"-120.4102\">RooRealVar</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"337.397\" y=\"-105.4102\">nobs_cr</text>\n",
       "</g>\n",
       "<!-- pdf_cr&#45;&gt;nobs_cr -->\n",
       "<g class=\"edge\" id=\"edge5\">\n",
       "<title>pdf_cr-&gt;nobs_cr</title>\n",
       "<path d=\"M337.397,-179.4253C337.397,-171.3504 337.397,-162.3471 337.397,-153.7362\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"340.8971,-153.5161 337.397,-143.5162 333.8971,-153.5162 340.8971,-153.5161\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nb_cr -->\n",
       "<g class=\"node\" id=\"node9\">\n",
       "<title>nb_cr</title>\n",
       "<ellipse cx=\"472.397\" cy=\"-116.6102\" fill=\"none\" rx=\"57.0522\" ry=\"26.7407\" stroke=\"#ff0000\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"472.397\" y=\"-120.4102\">RooProduct</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"472.397\" y=\"-105.4102\">nb_cr</text>\n",
       "</g>\n",
       "<!-- pdf_cr&#45;&gt;nb_cr -->\n",
       "<g class=\"edge\" id=\"edge6\">\n",
       "<title>pdf_cr-&gt;nb_cr</title>\n",
       "<path d=\"M370.7678,-184.1673C388.7788,-172.1946 411.1893,-157.2975 430.3792,-144.5411\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"432.5519,-147.2997 438.9422,-138.849 428.6767,-141.4701 432.5519,-147.2997\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- nb_cr&#45;&gt;nb_sr -->\n",
       "<g class=\"edge\" id=\"edge10\">\n",
       "<title>nb_cr-&gt;nb_sr</title>\n",
       "<path d=\"M438.883,-94.6573C420.5278,-82.6339 397.5912,-67.6096 377.9792,-54.763\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"379.7996,-51.7714 369.5166,-49.2197 375.9639,-57.627 379.7996,-51.7714\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "<!-- tau -->\n",
       "<g class=\"node\" id=\"node10\">\n",
       "<title>tau</title>\n",
       "<ellipse cx=\"472.397\" cy=\"-26.8701\" fill=\"none\" rx=\"59.2941\" ry=\"26.7407\" stroke=\"#0000ff\"/>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"472.397\" y=\"-30.6701\">RooRealVar</text>\n",
       "<text fill=\"#000000\" font-family=\"Times,serif\" font-size=\"14.00\" text-anchor=\"middle\" x=\"472.397\" y=\"-15.6701\">tau</text>\n",
       "</g>\n",
       "<!-- nb_cr&#45;&gt;tau -->\n",
       "<g class=\"edge\" id=\"edge9\">\n",
       "<title>nb_cr-&gt;tau</title>\n",
       "<path d=\"M472.397,-89.6852C472.397,-81.6102 472.397,-72.607 472.397,-63.9961\" fill=\"none\" stroke=\"#000000\"/>\n",
       "<polygon fill=\"#000000\" points=\"475.8971,-63.776 472.397,-53.776 468.8971,-63.7761 475.8971,-63.776\" stroke=\"#000000\"/>\n",
       "</g>\n",
       "</g>\n",
       "</svg>"
      ],
      "text/plain": [
       "<IPython.core.display.SVG object>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pdf.graphVizTree(\"pdf_graph.dot\")\n",
    "!dot -Tsvg pdf_graph.dot > pdf_graph.svg; rm pdf_graph.dot\n",
    "SVG(\"pdf_graph.svg\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.Eval) \n",
    "ROOT.RooMsgService.instance().setGlobalKillBelow(RooFit.ERROR)\n",
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.NumIntegration) \n",
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.Fitting) \n",
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.Minimization) \n",
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.InputArguments) \n",
    "ROOT.RooMsgService.instance().getStream(1).removeTopic(RooFit.Eval) \n",
    "ROOT.RooMsgService.instance().setGlobalKillBelow(RooFit.ERROR)\n",
    "ROOT.RooMsgService.instance().setSilentMode(True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  1) 0x55dfee23a5f0 RooRealVar:: nobs_sr = 7  L(0 - 1500)  \"nobs_sr\"\n",
      "  2) 0x55dfee1d0430 RooRealVar:: nobs_cr = 7  L(0 - 900)  \"nobs_cr\"\n"
     ]
    }
   ],
   "source": [
    "# create a useful function to create toys specific for this workspace\n",
    "def create_toy():\n",
    "    return pdf.generate(obs, 1)\n",
    "\n",
    "# just an example\n",
    "data = create_toy()\n",
    "data.get(0).Print(\"v\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from functools import wraps\n",
    "\n",
    "def stateless(varnames):\n",
    "    def decorator(method):\n",
    "        @wraps(method)\n",
    "        def f(*args, **kwargs):\n",
    "            old_values = {v: ws.var(v).getVal() for v in varnames}\n",
    "            result = method(*args, **kwargs)\n",
    "            for k, v in old_values.iteritems():\n",
    "                ws.var(k).setVal(v)\n",
    "                ws.var(k).setConstant(False)\n",
    "            return result\n",
    "        return f\n",
    "    return decorator\n",
    "\n",
    "@stateless([poi.GetName()])\n",
    "def evaluate_test_stat(test_value, data):\n",
    "   \n",
    "    # uncond fit\n",
    "    poi.setConstant(False)\n",
    "    r = pdf.fitTo(data, RooFit.Save())\n",
    "    denominator = r.minNll()\n",
    "    \n",
    "    # special case, we are not interested in these cases\n",
    "    if poi.getVal() > test_value:\n",
    "        return 0.\n",
    "    \n",
    "    # recompute if needed\n",
    "    if poi.getVal() < 0:\n",
    "        poi.setVal(0)\n",
    "        poi.setConstant(True)\n",
    "        r = pdf.fitTo(data, RooFit.Save())\n",
    "        denominator = r.minNll()\n",
    "        \n",
    "    # conditional fit\n",
    "    poi.setVal(test_value)\n",
    "    poi.setConstant(True)\n",
    "    r = pdf.fitTo(data, RooFit.Save())\n",
    "    numerator = r.minNll()\n",
    "    \n",
    "    return numerator - denominator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\r",
      "  0%|          | 0/10000 [00:00<?, ?it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing 0.1\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [01:38<00:00, 101.44it/s]\n",
      "100%|██████████| 10000/10000 [04:01<00:00, 32.06it/s]\n",
      "  0%|          | 4/10000 [00:00<05:52, 28.35it/s]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "computing 6.0\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 10000/10000 [06:11<00:00, 22.24it/s]\n",
      "100%|██████████| 10000/10000 [10:17<00:00, 14.14it/s]\n"
     ]
    }
   ],
   "source": [
    "nsignal_toys = 10000\n",
    "nbkg_toys = 10000\n",
    "\n",
    "test_poi_values = np.linspace(0.1, 6, 2)\n",
    "\n",
    "all_qvalue_data = {}\n",
    "all_qvalues_signal_toys = {}\n",
    "all_qvalues_bkg_toys = {}\n",
    "all_clsb = {}\n",
    "all_cls = {}\n",
    "all_clb = {}\n",
    "\n",
    "for test_poi_value in test_poi_values:\n",
    "    print \"computing\", test_poi_value\n",
    "    \n",
    "    qvalue_data = evaluate_test_stat(test_poi_value, data)\n",
    "\n",
    "    qvalues_signal_toys = []\n",
    "    qvalues_bkg_toys = []\n",
    "\n",
    "    for itoy in tqdm(range(nsignal_toys)):\n",
    "        # create signal-sample toy\n",
    "        poi.setVal(test_poi_value)\n",
    "        poi.setConstant(True)\n",
    "        pdf.fitTo(data)\n",
    "        poi.setConstant(False)\n",
    "        toy_s = create_toy()\n",
    "        qvalues_signal_toys.append(evaluate_test_stat(test_poi_value, toy_s))\n",
    "\n",
    "\n",
    "    for itoy in tqdm(range(nbkg_toys)):\n",
    "        # create signal-sample toy\n",
    "        poi.setVal(0)\n",
    "        poi.setConstant(True)\n",
    "        pdf.fitTo(data)\n",
    "        poi.setConstant(False)\n",
    "        toy_b = create_toy()\n",
    "        qvalues_bkg_toys.append(evaluate_test_stat(test_poi_value, toy_b))\n",
    "    \n",
    "    \n",
    "    all_qvalue_data[test_poi_value] = qvalue_data\n",
    "    all_qvalues_signal_toys[test_poi_value] = np.array(qvalues_signal_toys)\n",
    "    all_qvalues_bkg_toys[test_poi_value] = np.array(qvalues_bkg_toys)\n",
    "    \n",
    "    pvalue_sb = (all_qvalues_signal_toys[test_poi_value] >= qvalue_data).sum() / float(nsignal_toys)\n",
    "    pvalue_b = (all_qvalues_bkg_toys[test_poi_value] >= qvalue_data).sum() / float(nbkg_toys)\n",
    "    cls = pvalue_sb / (1 - pvalue_b)\n",
    "    \n",
    "    all_clsb[test_poi_value] = pvalue_sb\n",
    "    all_clb[test_poi_value] = 1 - pvalue_b\n",
    "    all_cls[test_poi_value] = cls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA3YAAAE/CAYAAAAKZVEAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X+4nHV54P/3TfiRiHBCSUpDDvFgg15iUGmPIAatYC20\nkIJbqvzwam1ZIgL6rd+uGlddj/v9tkYvr3VlRWtElsWtgGarEpIV3QXKBlJNQBBjlhrZ03ACSkLN\nqWBAwXv/OBM8JGfmzMnMPPM8M+/XdZ3rzDzPMzP3J3Myn7k/PyMzkSRJkiRV1wHdDkCSJEmS1BoT\nO0mSJEmqOBM7SZIkSao4EztJkiRJqjgTO0mSJEmqOBM7SZIkSao4EztJkiRJqjgTO6lgEfFrEfGV\niHgiIv4pIi6sc91pEXFbRIxHxGjBYUqS1BURcX5EbKnVkz+MiNfUue5dEfGjiPiXiLgmIg4pOlap\nTEzspOJdBfwcOAq4CPhMRLx0iuueAK4B3l1gbJIkdU1EvAH4KPBnwGHAa4EHp7juDGAF8HrgBcAL\ngQ8XF6lUPpGZ3Y5B6hsRcSjwE2BJZv5j7dgXgO2ZuaLOY34XuDozhwoLVJKkLoiIu4DPZ+bnp7nu\ni8BoZv7b2v3XA3+bmb9RQJhSKdljJxXrRcDTe5K6mvuAqXrsJEnqGxExCxgG5kfE1ogYi4hPRcSc\nKS5/KRP15x73AUdFxJFFxCqVkYmdVKznA/+y17FxJoabSJLUz44CDgLOA14DvAI4EfjAFNc+n4n6\nc489t61P1bdM7KRiPQ4cvtexw4GfdiEWSZLKZHft93/KzEcycyfwH4A/mOLavevTPbetT9W3TOyk\nYv0jcGBEHDfp2MuBzV2KR5KkUsjMnwBjwOQFIOotBrGZifpzj5cDP87MxzoUnlR6JnZSgTLzCeDv\ngH8fEYdGxFLgHOALe18bEQdExGwmhqVERMyOiIOLjViSpEL9Z+AdEfHrEXEE8C7g5imuuw64OCKO\nj4i5TAzXvLa4MKXyMbGTincZMAd4FLgeeHtmbo6I10TE45Ouey0Tw1LWAYtqt79RdLCSJBXo/wM2\nMjHCZQvwHeCvImJRRDweEYsAMvPrwMeA24BtwD8BH+pOyFI5uN2BJEmSJFWcPXaSJEmSVHEmdpIk\nSZJUcSZ2kiRJklRxJnaSJEmSVHEmdpIkSZJUcQd2O4BG5s2bl0NDQ90Oo6c8uOMJAF44/9AuRyJJ\nv3L33XfvzMz53Y6jKqwf28/6UVJZNVtHljqxGxoaYtOmTd0Oo6e8+bMbALjxbad0ORJJ+pWI+Kdu\nx1Al1o/tZ/0oqayarSMdiilJUkVExLKIWDU+Pt7tUCRJJWNiJ0lSRWTmmsxcPjAw0O1QJEklU8rE\nzhZJSZIkSWpeKefYZeYaYM3w8PAl3Y5FUm/6xS9+wdjYGE8++WS3Q+krs2fPZnBwkIMOOqjboUiS\npmD92D2t1pGlTOwkqdPGxsY47LDDGBoaIiK6HU5fyEwee+wxxsbGOPbYY7sdTiVFxDJg2eLFi7sd\niqQeZf3YHe2oI0s5FFOSOu3JJ5/kyCOPtNIqUERw5JFH2grcAufYSeo068fuaEcdaWInqW9ZaRXP\nf3NJKj8/q7uj1X93EztJqoiRkRE+/vGPdzsMSZJKxzrSOXaSBMDSlbeyfdfutj3fwrlzuHPF6fv9\n+KGhIUZHR9sWj3qDc+wkFa1s9SNYR9ZjYidJwPZduxldeVbbnm9oxdqG55944gne9KY3MTY2xjPP\nPMMHP/hB3vzmN0/7vPfddx+nnHIKO3fu5D3veQ+XXOLiwf3EVaMlFa3o+hGsI/eXiV2Pq9fKMrRi\nbVtaTCTtn69//escffTRrF07UcE1u2/nd7/7Xf7hH/6BJ554ghNPPJGzzjqLo48+upOhqoc1aom3\njpDULdaR+8fErsft3cry5s9uAODGt53SVIuJpM444YQT+Mu//Eve+973cvbZZ/Oa17yGyy+/nDvv\nvBOAhx9+mFe84hUA/PEf/zHvf//7ATjnnHOYM2cOc+bM4bTTTuPb3/425557btfKoWpr1BJvHSGp\nW6wj94+JnSR1wYte9CLuuece1q1bxwc+8AFe//rXc9VVVz17fmhoiHvvvXefx+29YpYrl0mSeo11\n5P4p5aqYEbEsIlY12+0qSVXz8MMP87znPY+3vOUtvPvd7+aee+5p6nFf+9rXePLJJ3nssce4/fbb\neeUrX9nhSCVJKpZ15P4pZY+dk8PbZ/0h74SRC3914KkPTPweOZP1h8wD2jcZVlLz7r//ft797ndz\nwAEHcNBBB/GZz3ymqce97GUv47TTTmPnzp188IMf7Ku5A3JVTEn9wTpy/5QysVP7DMZOGJnU81mb\nY8fbxhkcGehOUFIJLZw7p61zihbOndPw/BlnnMEZZ5xR9/xUyziPjIy0GJWqzoZPSUUrun4E68j9\nZWInSeDqf5IkTcH6sTpKOcdOkiRJktQ8EztJkiRJqjgTO0mSJEmqOBM7SZIqwu2AJEn1mNhJklQR\nmbkmM5cPDLiqsSTpuUzsJKlLRkdHWbJkSUvPcfvtt3P22We3KaL2eutb38rq1au7HYYkqWKsH/eP\n2x1IEsAnToDxbe17voFF8K772/d8HZCZZCYHHGAbnySpDuvHyjCxkySYqLRG2jhvaaS5oXJPP/00\nF110Effccw8vfelLue666/j4xz/OmjVr2L17N69+9av57Gc/S0SwdetWLr30Unbs2MGsWbP48pe/\n/Jzn2rhxI8uXL2f16tUcfvjhXHjhhTz88MOccsopfPOb3+Tuu+/m8ccf54wzzuDkk0/m7rvvZt26\nddx111389V//NZnJWWedxUc/+lEAnv/85/P4448DsHr1am6++WauvfZa3vrWt3L44YezadMmfvSj\nH/Gxj32M8847j8zkHe94B9/85jc55phjOPjgg9v37ylJ6g7rx8rUj9VKQyWpxzzwwANcdtllbNmy\nhcMPP5xPf/rTXHHFFWzcuJHvfe977N69m5tvvhmAiy66iMsvv5z77ruPu+66iwULFjz7PHfddReX\nXnopX/va1/jN3/xNPvzhD3P66aezefNmzjvvPLZt+1Vr6w9+8AMuu+wyNm/ezEEHHcR73/tebr31\nVu699142btzIV7/61WnjfuSRR1i/fj0333wzK1asAOArX/kKDzzwAN///ve57rrruOuuu9r8ryVJ\n6hfWjzNnYidJXXTMMcewdOlSAN7ylrewfv16brvtNk4++WROOOEEbr31VjZv3sxPf/pTtm/fzhvf\n+EYAZs+ezfOe9zwAtmzZwvLly1mzZg2LFi0CYP369Zx//vkAnHnmmRxxxBHPvuYLXvACXvWqVwET\nrZive93rmD9/PgceeCAXXXQRd9xxx7Rxn3vuuRxwwAEcf/zx/PjHPwbgjjvu4IILLmDWrFkcffTR\nnH766W36V5Ik9Rvrx5kzsZOkLoqIfe5fdtllrF69mvvvv59LLrmEJ598suFzLFiwgNmzZ/Od73yn\nqdc89NBDZxzb3jEccsghz97OzKaeT5KkZlk/zpyJnSR10bZt29iwYQMAX/ziFzn11FMBmDdvHo8/\n/vizq2YddthhDA4OPjsM5KmnnuJnP/sZAHPnzmXt2rW8733v4/bbbwdg6dKlfOlLXwLgG9/4Bj/5\nyU+mfP2TTjqJv//7v2fnzp0888wzXH/99fzO7/wOAEcddRRbtmzhl7/8JV/5ylemLctrX/tabrzx\nRp555hkeeeQRbrvttv38V1E97mMnqV9YP85cYYldRJwbEZ+LiBsj4veKel1JKrMXv/jFXHXVVbzk\nJS/hJz/5CW9/+9u55JJLWLJkCWeccQavfOUrn732C1/4AldeeSUve9nLePWrX82PfvSjZ88dddRR\n3HzzzVx++eV861vf4kMf+hDf+MY3WLJkCV/+8pf5jd/4DQ477LB9Xn/BggWsXLmS0047jZe//OX8\n9m//Nueccw4AK1eu5Oyzz+bVr371c+Yr1PPGN76R4447juOPP54/+ZM/4ZRTTmnDv5Amcx87Sf3C\n+nHmopUuwoi4BjgbeDQzl0w6fibwSWAWcHVmrpx07gjg45l58XTPPzw8nJs2bdrv+MTEykOTVjJ6\n82cnWj5ufNsp+5yT+smWLVt4yUte8qsDPbac81NPPcWsWbM48MAD2bBhA29/+9u59957uxbPZPv8\n2wMRcXdmDncppMppV/04tGItoyvPmvG5XvSc+lHqY9aP3dVKHdnqdgfXAp8Crpv0wrOAq4A3AGPA\nxoi4KTO/X7vkA7XzklQeJd9TZ6a2bdvGm970Jn75y19y8MEH87nPfa7bIUmSqsj6sTJaSuwy846I\nGNrr8EnA1sx8ECAibgDOiYgtwErgv2fmPa28riSpseOOO67pyeKSJPWLXq4fOzHHbiHw0KT7Y7Vj\n7wB+FzgvIi6t9+CIWB4RmyJi044dOzoQniRJkiT1llaHYjYtM68ErmziulXAKpiYQ9DpuCT1r8zc\nZzlldZZbI0hS+Vk/dkerdWQneuy2A8dMuj9YOyZJpTF79mwee+wxE40CZSaPPfYYs2fP7nYokqQ6\nrB+7ox11ZCd67DYCx0XEsUwkdOcDF87kCSJiGbBs8eLFHQhPkmBwcJCxsTEc8l2s2bNnMzg42O0w\nJEl1WD92T6t1ZEuJXURcD7wOmBcRY8CHMvPzEXEFcAsT2x1ck5mbZ/K8mbkGWDM8PHxJK/FJUj0H\nHXQQxx57bLfDkCSpVKwfq6vVVTEvqHN8HbCuleeWJEmSJDWnE3PsJElSB0TEsohYNT4+3u1QJEkl\nU8rEzopLkqR9ZeaazFw+MDDQ7VAkSSVTysTOikuSJEmSmlfKxE6SJEmS1LxSJnYOxZQkSZKk5pUy\nsXMopiRJkiQ1r5SJnSRJkiSpeSZ2kiRJklRxJnaSJEmSVHGlTOxcPEWSJEmSmlfKxM7FUyRJkiSp\neaVM7CRJkiRJzTOxkyRJkqSKM7GTJEmSpIorZWLn4imSJEmS1LxSJnYuniJJkiRJzTuw2wFIktTP\nIuJQ4NPAz4HbM/Nvi3rt9Ye8E0YurHNuHnBWUaFIklpkYidJUptFxDXA2cCjmblk0vEzgU8Cs4Cr\nM3Ml8K+A1Zm5JiJuBApL7AZjJ4xMPe1hcMRRM5JUJaUciilJUsVdC5w5+UBEzAKuAn4fOB64ICKO\nBwaBh2qXPVNgjJKkHmJiJ0lSm2XmHcA/73X4JGBrZj6YmT8HbgDOAcaYSO7AelmStJ9KWYG4KqYk\nqQct5Fc9czCR0C0E/g74o4j4DLCm3oMjYnlEbIqITTt27OhspJKkyillYueqmJKkfpGZT2Tmn2Xm\n2xstnJKZqzJzODOH58+fX2SIkqQKKGViJ0lSD9oOHDPp/mDtmCRJLTOxkySpGBuB4yLi2Ig4GDgf\nuGkmT+BUBUlSPSZ2kiS1WURcD2wAXhwRYxFxcWY+DVwB3AJsAb6UmZtn8rxOVZAk1eM+dpIktVlm\nXlDn+DpgXcHhSJL6gD12kiRVhEMxJUn1mNhJklQRDsWUJNVjYidJkiRJFVfKxM6hJpIkSZLUvFIm\ndg41kSRpXzZ8SpLqKWViJ0mS9mXDpySpHhM7SZIkSao4EztJkiRJqjgTO0mSKsI5dpKkekzsJEmq\nCOfYSZLqMbGTJEmSpIo7sNsBqHqWrryV7bt2T3lu4dw53Lni9IIjkiRJkvqbiZ1m7MbdlzA4e+eU\n58Z2zwN+WGxAkiRJUp8zsdOMDcZOGJl64v7giPM+JKlTImIZsGzx4sXdDqU8PnECjG+b+tzAInjX\n/cXGI0ldUsrEzopLkqR9ZeYaYM3w8PAl3Y6lNMa31W1sxMZGSX2klImdFZckSd03tGLtlMedTy1J\n5VPKxE6SJHXf6MqzpjxeL+GTJHWP2x1IkiRJUsXZYydJktrHxUwkqStM7CRJqohKLC7mYiaS1BUm\ndpIkVURZFhdbOHdO3Xl2o7MLDkaSBJjYSZKkqQwsqtvDdidAnQRuLOcx2LGg+pBDWyU1ycROxbFy\nkqTq2M/P5FNXrGW03skGyWLh9cDeddJTH5j4PXJmueokh7ZKapKJnYozvo2hJ7845alRLiw4GElS\n4RolS0UnKXsnTJ/dMPH7beP7H4sNmJK6yMROhaq3JxIjhYYhSVL72bsmqYtM7FRtto5KkiRJJnYq\nj3orrC2cO4c7V5w+9YM60Dq6dOWtbN+1e+axmGRKUsMVMxt+hpZp/p0kVZCJnUqj3jDNel8Qpju/\nv0tu37j7EgZn75zy3NjuecAPp36gcwgldVgV9rGrm7gxzed5mebfSVIFmdipHBq01I7Opu4cvLGc\n1/Z5e4Oxs24v4OA0Xy6cQyipk8qyj50kqXxM7FQODVpqh1asrZswNVxWW5IkTc3pA1LPKSyxi4gX\nAu8HBjLzvKJeV9U33XwNSZI0Q67gKfWclhK7iLgGOBt4NDOXTDp+JvBJYBZwdWauzMwHgYsjYnUr\nr6n+02i+hiRJhZqup6vX2dMnlVarPXbXAp8CrttzICJmAVcBbwDGgI0RcVNmfr/F15IkSequRj1d\n/cCePqm0WkrsMvOOiBja6/BJwNZaDx0RcQNwDmBi1yfavUqlJKn3NdpqxvpDkqbXiTl2C4GHJt0f\nA06OiCOBvwJOjIj3ZeZHpnpwRCwHlgMsWtQHQxp6kCtDSpJmavuu3dYfktSCwhZPyczHgEubuG4V\nsApgeHg4Ox2X2my6DWYlSZIktV0nErvtwDGT7g/WjqkfOGlakqR9Tdfwaf0pqUWdSOw2AsdFxLFM\nJHTnAxfO5AkiYhmwbPHixR0IT5IkqWCNEjcXHZHUBq1ud3A98DpgXkSMAR/KzM9HxBXALUxsd3BN\nZm6eyfNm5hpgzfDw8CWtxCc12v/ObRQkST2t37dmkPpMq6tiXlDn+DpgXSvPLbVsYBGjdTqLx3bP\nA35YbDyS1CJHtGhGit6aweGmUlcVtnjKTFhxqS0aVCCDDnuRVEGOaFGpOdxU6qpSJnZ9WXFNN1zC\nVi5JknqTK0rX5/cjqWmlTOz6UqPhErZytd1YzqvbazeW8xgsOB5JUh8zOanP70dS00zs1JdOferK\nuhvhnrpiLaPFhiNJaqBRY1zDHi17wvqXPX3qQ6VM7Jxjp36ydOWtbN+1e8pzrt4pSY0b4xryy3v/\nsqdPfaiUiV1fzrFT37px9yUMzt455TlX75TUSxbOndNwGxpJ0v4rZWIn9ZPB2Fm3VdHVOyX1Ekcg\n9DG3QpA6zsROkiRJneVWCFLHmdipLzkcSJIkSb2klIldYYunuGJS33I4kCSpLlfTlFRBpUzsCls8\nxRWTJEnS3mzYlVRBpUzsJEmSmuGWMR1gj6VUSSZ2kiSpsrbv2l13j7ulK2+tO5+6nqEVa00I7bGU\nKsnETiq5Rl9K+v7LhyQ1MJPPxzd/dgMAN77tlBkng13TqGdtz3lJfaOUiV1hi6dIZTewiFEurHu6\nExuYO6xJKl5EvBB4PzCQmed1O56ycSXjOuxZkzRJKRO7whZPUV2NvtyPzi44mH42TaXdiQ3MGw1r\nqkwrtlSgiLgGOBt4NDOXTDp+JvBJYBZwdWaurPccmfkgcHFErO50vFVkg5IkTa+UiZ26r9GXe0YK\nDUWSyu5a4FPAdXsORMQs4CrgDcAYsDEibmIiyfvIXo//88x8tJhQJTXkVliqMBM7TWn9Ie+EkTpD\nAB2z39MavffrD5kH1En4pT6VmXdExNBeh08CttZ64oiIG4BzMvMjTPTuSSojt8JShZnYaUqDsbP+\nB5t6WqP3vhNDP6UetRB4aNL9MeDkehdHxJHAXwEnRsT7agng3tcsB5YDLFpkA5sk6blM7Io0Xfe+\n1E4OJ5EqIzMfAy6d5ppVwCqA4eHhLCIuSVJ1lDKx69lVMRt170vt5nASqZu2A8dMuj9YOyZJUkeU\nMrFzVUypOY8wnwV1krRHmM+CguOR9KyNwHERcSwTCd350GDvkib1bMOnJKllpUzspDKrt+R/N7aB\nWDCyte65U1asZbS4UKS+FRHXA68D5kXEGPChzPx8RFwB3MLESpjXZObmVl/Lhk/1pEYbrZdpqsp0\ncdab4uDUCBXExE6aIbeBkDRZZl5Q5/g6YF3B4UjVU5XEplGcjaY4ODVCBTGxk2aiKq2KknqSQzGl\nNrAuV48ysZNmoiqtipJ6kkMxpTawLlePMrGTJElq0tKVt7J91+4ZP27h3DncueL0DkQkSRNM7CRJ\nkpq0fdfu+nOtG6i38JYktcsB3Q5gKhGxLCJWjY+755skSXtYP0qS6illYpeZazJz+cCAKwVJkrSH\n9aMkqZ5SJnaSJEmSpOaZ2EmSJElSxbl4iiRJFeE+dv2n0SqcrrQpaTITO6lHLZw7p+4qbKOzCw5G\nUlu4j13/abQKpyttSprMxE7qUQ1bcUcKC0OSJEkFcI6dJEmSJFVc7/fYfeIEGN829bmBRcXGIkmS\nJEkd0PuJ3fg2GHEjV0lS9bl4itRHGnVOwEQHxbvuLy4elV7vJ3aS9jGW8xgcmXqD47Gcx2DB8Uhq\njounSH1kus6JOvW4+peJndSHTn3qyrqrrJ26Yi2jxYYjSZKkFpUysXOoiVROjXr6HBIiSZLUPaVM\n7BxqIpVTo56+/R0S4ua7ksqm0T6gC+fOKTgaqY6BRfXrXhtb+1IpEztJ1dToi1C9BM3NdyWVjQ1K\nqoRGiZvz7/qSiZ2kthjLeYzOvnDqc7vnAT8sNiBJkqQ+YmLXbu6bpx7WeHjS5+q2ctedlwesP+Sd\nMDJ1Qrj+kHlAnaGfUh9yDrokqR4Tu3Zz3zz1sE4MTxqMnXX/zzRKCKV+5Bx0tYvzm6XeY2In9SEX\nBpCk8pguyeoE5zdLvcfETupDtsRKUnk0SrIkqVkHdDsASZIkSVJr7LGTJEnqMIfAS+o0EztJkqQO\ncwi8pE4zsZPUcWM5r+4Kl2M5j8GC45EkqRQGFtXfTNxtsjRDJnaSOu7Up66suzDAqSvWMlrncY0S\nQgYWwbvub0t8UlW4j53UYzpVj02XMFp/9iQTO0ml1SghrFthST3MfewkNaVR4mb92bNcFVOSJEmS\nKs4eO0nV5DATSZKkZxWW2EXEocCngZ8Dt2fm3xb12pJ6kMNMJKlwS1feyvZdu6c8t3DuHFf/lLqo\npcQuIq4BzgYezcwlk46fCXwSmAVcnZkrgX8FrM7MNRFxI2BiJ0mSVCHbd+2uO/e53j59KhlHvPSs\nVnvsrgU+BVy350BEzAKuAt4AjAEbI+ImYBDY85fyTIuvK6kPTLehry3DkiTNkCNeelZLiV1m3hER\nQ3sdPgnYmpkPAkTEDcA5TCR5g8C9uGiLpCY0StxsGZYkSfqVTiRYC4GHJt0fqx37O+CPIuIzwJp6\nD46I5RGxKSI27dixowPhSZIkSVJvKWzxlMx8AvizJq5bBawCGB4ezk7HJUmSJElV14keu+3AMZPu\nD9aOSZIkSZI6oBOJ3UbguIg4NiIOBs4HbprJE0TEsohYNT4+3oHwJEmqJutHSVI9LSV2EXE9sAF4\ncUSMRcTFmfk0cAVwC7AF+FJmbp7J82bmmsxcPjDgyjySJO1h/ShJqqfVVTEvqHN8HbCuleeWeskj\nzGdBnSWEH2E+CwqOp2jTbVsgSZKk1hS2eMpMRMQyYNnixYu7HYrUFgtGttY/V2Ac3eJ+c5IkSZ1V\nysQuM9cAa4aHhy/pdiySJElVtHTlrWzftXvKc46WkHpPKRM7SZIktWb7rt2Mrjyr22FIKoiJnSRJ\nUgU1mr+857yk/lHKxM45dpI65hMnwPi2qc8NLIJ33V9sPJK0n5y/rLYbWAR1Fnub9nHWn11XysTO\nOXaSOmZ8G4zU2QNsfyozSZJ6xf4mZ9afpVDKxE6SJEnVMt3WNvYwSp1lYiepkhp9gRid3fix+/s4\nSVJ9jRK3RnMB1QMaDeF0mGZhTOwkVVLDlt+Rxo+tu0rcNI+TJElTaJS4OUyzMKVM7Fw8RVKr7JWT\nJEn9pJSJnYunSGqVvXKSJKmflDKxk6SyWbryVrbv2j3lORcFkCSpDrdQKIyJnSQ1Yfuu3XV7AV0U\nQK2KiHOBs4DDgc9n5je6HJIktYdbKBTGxE6SpBZExDXA2cCjmblk0vEzgU8Cs4CrM3NlvefIzK8C\nX42II4CPAyZ2UovcekH9ppSJnYunSJIq5FrgU8B1ew5ExCzgKuANwBiwMSJuYiLJ+8hej//zzHy0\ndvsDtcdJapGjLNRvSpnYlWLxFPfjkCQ1ITPviIihvQ6fBGzNzAcBIuIG4JzM/AgTvXvPEREBrAT+\ne2be09mIJbXddN8bpQKUMrErBffjkCTtv4XAQ5PujwEnN7j+HcDvAgMRsTgz/2bvCyJiObAcYNEi\nvyhKpWKDv0rAxE6SpC7LzCuBK6e5ZhWwCmB4eDiLiEv9aeHcOQ3npxWp0YrE7ksqPZeJnSRpWm73\nMGPbgWMm3R+sHZNKr0z/nxutSOy+pNJzmdhJkqbldg8zthE4LiKOZSKhOx+4sNUndXExSVI9B3Q7\ngKlExLKIWDU+Pt7tUCRJaigirgc2AC+OiLGIuDgznwauAG4BtgBfyszNrb5WZq7JzOUDA871liQ9\nVyl77EqxKqYktYFDGHtfZl5Q5/g6YF3B4UiS+lQpEztJ6hUOYZQkSUUo5VBMSZK0L6cqSJLqscdO\nkpqw/pB3wsjUa1+sP2QeUGfVNqmNnKogSarHxE6SmjAYO2Fk6l6SwREXspAkqa0GFkG9+nVgUXk2\nhf/ECTC+bepzBcdpYidJUkW43YGkvtEoISpTg+r4troNv0XHaWInqfdM18rXZo8wnwV1Xm/D7PnU\nG6bZaHgnUK4WSZWCQzElSfWUMrGzRVJSSwpOhhaMbG1wrn5rXaPhnUC5WiQlqUcsnDun7qrEbkOj\nKitlYmeLpCRJkjqhUeLmNjSqMrc7kCRJkqSKM7GTJKki3MdOklRPKYdiSpKkfTlVQVVVlXltVYmz\n71VlK4SCmdhJkiSpo6oyr60qcfa9qmyFUDCHYkqSJElSxdljJ0md1GC4yFjOY7DgcCRJUm8ysZOk\nTmowXOSpE7weAAANxklEQVTUFWsZLS4S9QD3eVUvmm5eW1k4/05lZ2InSVJFuHiKelFVEiLn36ns\nnGMnSZIkSRVnYidJkiRJFVfKoZjOIShOvaEDo7MLDkRSqa0/5J0wcmGdc/OAs4oNSJIkPUcpEzvn\nEBRndGWdL2MjhYYhqeQGYyeMjE99ro/3DJKk6SxdeSvbd+2e8pyLrqidSpnYSZKkfTmiRaqe7bt2\n121Id9EVtZNz7CRJqojMXJOZywcG7CWVJD2XiZ0kSZIkVZyJnSRJkiRVnImdJEmSJFWci6f0s4FF\nUG81u4FFxcYilcBYzqu7wuNYzmOw4HgkSZKaZWLXz951f7cjkErl1KeurLty2akr1jJabDiSpAZs\njJOey8ROkiRJlWNjnPRczrGTJKkiImJZRKwaH596s3hJUv8ysZMkqSLcx06SVI+JnSRJkiRVnImd\nJEmSJFWciZ0kSZIkVZyJnSRJkiRVXGHbHUTEC4H3AwOZeV5RrytJkiRVzdKVt7J91+665xfOncOd\nK04vMCKVXVOJXURcA5wNPJqZSyYdPxP4JDALuDozV9Z7jsx8ELg4Ila3FrIkSZLU27bv2l13nz6A\noRVrC4xGVdBsj921wKeA6/YciIhZwFXAG4AxYGNE3MREkveRvR7/55n5aMvRSpIkSZL20VRil5l3\nRMTQXodPArbWeuKIiBuAczLzI0z07kmSJEmSCtDK4ikLgYcm3R+rHZtSRBwZEX8DnBgR72tw3fKI\n2BQRm3bs2NFCeJIkSZLUHwpbPCUzHwMubeK6VcAqgOHh4ex0XJIkVUVELAOWLV68uNuhSJJKppUe\nu+3AMZPuD9aOSZKkDsjMNZm5fGBgoNuhSJJKppXEbiNwXEQcGxEHA+cDN7UjqIhYFhGrxsfH2/F0\nkiRJktTTmkrsIuJ6YAPw4ogYi4iLM/Np4ArgFmAL8KXM3NyOoGyRlCRJkqTmNbsq5gV1jq8D1rU1\nIkmSJEnSjLQyFLNjHIopSZIkSc0rbFXMmcjMNcCa4eHhS7odiyR109CKtVMeXzh3DneuOL3gaCRJ\nZbFw7hzrCD1HKRM7SdKE0ZVnTXm8XmUuSeoPjRI364j+VMqhmJIkSZKk5pUysXOOnSRJkiQ1r5SJ\nndsdSJIkSVLzSpnYSZIkSZKaZ2InSZIkSRVXysTOOXaSJEmS1LxSJnbOsZMkSZKk5pUysZMkSZIk\nNc/ETpKkLouIl0TE30TE6oh4e7fjkSRVj4mdJEktiIhrIuLRiPjeXsfPjIgHImJrRKxo9ByZuSUz\nLwXeBCztZLySpN5kYidJUmuuBc6cfCAiZgFXAb8PHA9cEBHHR8QJEXHzXj+/XnvMHwJrgXXFhi9J\n6gUHdjuAqUTEMmDZ4sWLux2KJEkNZeYdETG01+GTgK2Z+SBARNwAnJOZHwHOrvM8NwE3RcRa4Iud\ni1iS1ItK2WPnqpiSpIpbCDw06f5Y7diUIuJ1EXFlRHyWOj12EbE8IjZFxKYdO3a0N1pJUuWVssdO\nkqR+kpm3A7dPc80qYBXA8PBwdj4qSVKVlLLHTpKkitsOHDPp/mDtmCRJHWFiJ0lS+20EjouIYyPi\nYOB84KZWnzQilkXEqvHx8ZYDlCT1FhM7SZJaEBHXAxuAF0fEWERcnJlPA1cAtwBbgC9l5uZWX8s5\n6JKkeko5x85VMSV1w8K5cxhasbbuuSJfD2DD7PksGJn6C/yG2fMZqrMz2sK5c7hzxentCFFNyMwL\n6hxfh1sXSOqC6eqXTrye9U73lTKxy8w1wJrh4eFLuh2LpP5RdKU0/eudVffMgpEBRldOfb7IylzF\nsuFTUjOKrs+sd8rBoZiSJFWEQzElSfWY2EmSJElSxZnYSZIkSVLFmdhJklQRbncgSarHxE6SpIpw\njp0kqR4TO0mSJEmquFImdg41kSRJkqTmlTKxc6iJJEn7suFTklRPKRM7SZK0Lxs+JUn1mNhJkiRJ\nUsWZ2EmSJElSxUVmdjuGuiJiB/BPBb/sPGBnwa9ZpF4uXy+XDSxflfVy2aA95XtBZs5vRzD9oI31\nY6//bUJ/lBH6o5z9UEawnL2kXWVsqo4sdWLXDRGxKTOHux1Hp/Ry+Xq5bGD5qqyXywa9X75e1g/v\nXT+UEfqjnP1QRrCcvaToMjoUU5IkSZIqzsROkiRJkirOxG5fq7odQIf1cvl6uWxg+aqsl8sGvV++\nXtYP710/lBH6o5z9UEawnL2k0DI6x06SJEmSKs4eO0mSJEmquL5M7CLi1yLimxHxg9rvI+pc96e1\na34QEX9aO/a8iFgbEf87IjZHxMpio59aRJwZEQ9ExNaIWDHF+UMi4sba+W9FxNCkc++rHX8gIs4o\nMu5m7W/5IuINEXF3RNxf+3160bE3o5X3r3Z+UUQ8HhH/pqiYm9Xi3+bLImJD7f/a/RExu8jYm9HC\n3+ZBEfFfauXaEhHvKzr26TRRttdGxD0R8XREnLfXuX0+P1Ue0723vSAiromIRyPie92OpVMi4piI\nuC0ivl/7nPx/uh1TJ0TE7Ij4dkTcVyvnh7sdU6dExKyI+E5E3NztWDolIkZrdd+9EbGp2/F0SkTM\njYjVtZxhS0Sc0vEXzcy++wE+Bqyo3V4BfHSKa34NeLD2+4ja7SOA5wGn1a45GPhfwO93uTyzgB8C\nL6zFdB9w/F7XXAb8Te32+cCNtdvH164/BDi29jyzuv0etbF8JwJH124vAbZ3uzztLN+k86uBLwP/\nptvlaeN7dyDwXeDltftH9tjf5oXADbXbzwNGgaFul2mGZRsCXgZcB5w36fiUn5/dLpM/zb+3vfAD\nvBb4LeB73Y6lg2VcAPxW7fZhwD/26HsZwPNrtw8CvgW8qttxdais/y/wReDmbsfSwTKOAvO6HUcB\n5fwvwL+u3T4YmNvp1+zLHjvgHCb+san9PneKa84AvpmZ/5yZPwG+CZyZmT/LzNsAMvPnwD3AYAEx\nN3ISsDUzH6zFdAMTZZxscplXA6+PiKgdvyEzn8rM/wNsrT1fmex3+TLzO5n5cO34ZmBORBxSSNTN\na+X9IyLOBf4PE+Urm1bK9nvAdzPzPoDMfCwznyko7ma1Ur4EDo2IA4E5wM+Bfykm7KZMW7bMHM3M\n7wK/3OuxU35+FhG0mtLM323lZeYdwD93O45OysxHMvOe2u2fAluAhd2Nqv1ywuO1uwfVfnpukYiI\nGATOAq7udixqTUQMMNG49HmYyBkyc1enX7dfE7ujMvOR2u0fAUdNcc1C4KFJ98fY68MyIuYCy4D/\n2YkgZ2DaWCdfk5lPA+NM9IA089hua6V8k/0RcE9mPtWhOPfXfpcvIp4PvBco67CUVt67FwEZEbfU\nhvu9p4B4Z6qV8q0GngAeAbYBH8/MMn0JbeWzoQqfK/3M96cH1YZ5n8hEb1bPqQ1RvBd4lImGo14s\n538E3sO+jWW9JoFv1KbILO92MB1yLLAD+M+1obVXR8ShnX7Rnk3sIuJ/RMT3pvjZu8U52Y9Wn1or\n+/XAlZn5YJvCVodExEuBjwJv63YsbTYCfGJSS2YvORA4Fbio9vuNEfH67obUVicBzwBHM1EB/GVE\nvLC7IUmqoloj338D/iIzy9Tz3zaZ+UxmvoKJUVInRcSSbsfUThFxNvBoZt7d7VgKcGpm/hbw+8Dl\nEfHabgfUAQcyMRT8M5l5IhMNuR2fz9yziV1m/m5mLpni52vAjyNiAUDt96NTPMV24JhJ9wdrx/ZY\nBfwgM/9jp8owA9PF+pxraknpAPBYk4/ttlbKt2dow1eAP8nMH3Y82plrpXwnAx+LiFHgL4B/GxFX\ndDrgGWilbGPAHZm5MzN/Bqxj4kOyTFop34XA1zPzF5n5KHAnMNzxiJvXymdDFT5X+pnvTw+JiIOY\nSOr+NjP/rtvxdFptONtt9N7w7qXAH9bq8xuA0yPiv3Y3pM7IzO21348y8f2sbFOA2mEMGJvUs7ya\nAr7D9GxiN42bgD2rtP0p8LUprrkF+L2IOCImVs38vdoxIuL/Z+LL2V8UEGszNgLHRcSxEXEwEws0\n3LTXNZPLfB5wa6238ibg/NrKfccCxwHfLijuZu13+WrDZdcysVjOnYVFPDP7Xb7MfE1mDmXmEBND\nOP46Mz9VVOBNaOVv8xbghJhYifZA4HeA7xcUd7NaKd824HSA2vCMVwH/u5Com9NM2eqp+/mpUmjl\nvVWJ1Obrfh7Ykpn/odvxdEpEzK/V50TEHOANlOvzsmWZ+b7MHKzV5+czUVe8pcthtV1EHBoRh+25\nzUT90HMr12bmj4CHIuLFtUOvp4jvMJ1enaWMP0zMb/mfwA+A/wH8Wu34MHD1pOv+nInFRLYCf1Y7\nNsjE0M0twL21n39dgjL9AROrYf0QeH/t2L8H/rB2ezYTqyZuZSJxe+Gkx76/9rgH6PIKn+0uH/AB\nJrq/75308+vdLk87379JzzFCyVbFbMPf5luYWBTme8DHul2WNv9tPr92fDMTH/bv7nZZ9qNsr2Si\nVfIJJnohN0967D6fn/6U52eq97bXfpiYLvEI8Iva3+nF3Y6pA2U8lYnvJN+dVMf9Qbfj6kA5XwZ8\np1bO7wH/rtsxdbi8r6NHV8VkYjXe+2o/m3v186dW1lcAm2p/t1+lgNWho/bCkiRJkqSK6tehmJIk\nSZLUM0zsJEmSJKniTOwkSZIkqeJM7CRJkiSp4kzsJEmSJKniTOwkSZIkqeJM7CRJkiSp4kzsJEmS\nJKni/i8RfcZRjml2aAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb9d78ee210>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, len(test_poi_values), figsize=(15, 5))\n",
    "for ax, test_poi_value in zip(axs.flat, test_poi_values):\n",
    "    \n",
    "    qvalue_data = all_qvalue_data[test_poi_value]\n",
    "    qvalues_signal_toys = all_qvalues_signal_toys[test_poi_value]\n",
    "    qvalues_bkg_toys = all_qvalues_bkg_toys[test_poi_value]\n",
    "    \n",
    "    bins = np.linspace(min([np.min(qvalues_signal_toys), np.min(qvalues_bkg_toys)]),\n",
    "                       max([np.max(qvalues_signal_toys), np.max(qvalues_bkg_toys)]), 50)\n",
    "    ax.hist(qvalues_signal_toys, bins=bins, histtype='step', label='s+b', normed=True)\n",
    "    ax.hist(qvalues_bkg_toys, bins=bins, histtype='step', label='background', normed=True)\n",
    "    ax.axvline(qvalue_data)\n",
    "    ax.set_title(test_poi_value)\n",
    "    ax.set_yscale('log')\n",
    "    ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb9d7f32fd0>"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlYAAAGfCAYAAACdqpz+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl4VOXB/vHvmexASMhCIMlABIkQEiAYNiEkA6KCRcDW\nuuuLtWopxYX2qhV9RQu29ieoVKsttq9afbVVCyKiIpCNNRBEsgCyZ2EJJCRk3+b8/ojNq1b2SU4y\nuT/XxRUmc3Lmnhbhnuc853kM0zQRERERkUtnszqAiIiIiLtQsRIRERFxERUrERERERdRsRIRERFx\nERUrERERERdRsRIRERFxERUrERERERdRsRIRERFxERUrERERERfxtOqFQ0JCzKioKKteXkREROS8\nZWVlnTRNM/Rcx1lWrKKioti2bZtVLy8iIiJy3gzDOHw+x+lSoIiIiIiLqFiJiIiIuIiKlYiIiIiL\nWDbHSkRERDqmhoYGCgsLqa2ttTqKy/n6+hIZGYmXl9dF/byKlYiIiFyQwsJC/P39iYqKwjAMq+O4\njGmalJSUUFhYyGWXXXZR59ClQBEREbkgtbW1BAcHu1WpAjAMg+Dg4EsaiVOxEhERkQvmbqXq3y71\nfalYiYiIiLjIOYuVYRh/Mwyj2DCMnDM8bxiGscQwjH2GYew0DGO462OKiIiInN38+fN57rnnLM1w\nPiNWrwPXneX5ycCAr3/dB7xy6bFEREREOp5zFivTNNOB0rMcMg1402y2GQg0DKO3qwKKiIhIx5d1\n+BQvp+wj6/Apl53zzTffZMiQIQwdOpQ777zzW88tWbKEmJgYhgwZwi233OKy1zwXVyy3EAEUfONx\n4dffO+qCc4uIiEg79tRHueQdOX3WYypqG9h9rAKnCTYDBvbyx9/3zOtExYR358mpg896ztzcXBYs\nWMDGjRsJCQmhtLSUJUuWtDz/+9//noMHD+Lj40NZWdmFvalL0KaT1w3DuM8wjG2GYWw7ceJEW760\niIiIWOR0bSNOs/n3TrP58aVat24dN910EyEhIQAEBQV96/khQ4Zw++2389Zbb+Hp2XbLdrrilYoA\n+zceR379vf9gmuZfgL8AJCQkmC54bREREbHQuUaWoPky4O2vbaah0YmXp40Xb4nnyr49WjXXxx9/\nTHp6Oh999BELFy4kOzu7TQqWK0asVgB3fX134Gig3DRNyy8DbjtUyrxl2S69lisiIiIX7sq+PXj7\n3tE8cs0VvH3vaJeUqgkTJvDee+9RUlICQGnp/00HdzqdFBQU4HA4ePbZZykvL6eysvKSX/N8nLO6\nGYbxDpAMhBiGUQg8CXgBmKb5KrAKmALsA6qBma0V9nxlHT7FrUs309Bk8k5mPn+5K4GrB4VZHUtE\nRKTTurJvD5eOUg0ePJh58+aRlJSEh4cH8fHxREVFAdDU1MQdd9xBeXk5pmkyZ84cAgMDXfbaZ3PO\nYmWa5q3neN4Efu6yRC6w+UAJTV9fzHWaMOvt7TzxgxhuH9kHm809V4oVERHpbO6++27uvvvu731u\n/fr1bZymmVuuvD66XzDenjY8DPDxtHFFWDeeWJ7Dja9sJPdIudXxRERExE213TT5NvTva7mbD5Qw\nul8ww/sEsnxHEQtW7uKGlzYw86ooHp4UTVcft3z7IiIiYhG3bRbfvZY7Iz4SxxU9efbT3by2/iAf\nZx9l/g2DuXZwLwtTioiIiDtxy0uBZxLYxZvf3TiED342hu6+Xtz/9yzufWMbRWU1VkcTERERN9Cp\nitW/Xdk3iJVzxvHo5IGs33eCSYvTWJp+gMYmp9XRREREpAPrlMUKwMvDxgNJ/fn84STG9Atm4apd\nTH1pA9vzte6ViIiIXJxOW6z+zR7UhdfuTuDVO4ZzqqqeH76ykXnLsimvabA6moiIiFyA+fPn89xz\nzwGQnJzMtm3b2jxDpy9WAIZhcF1sb9bMTWLmVZfxTmY+Exel8eGOIpqX6RIRERE5NxWrb+jm48l/\nT41hxexxhAf68uC7O7jrb5kcOllldTQREZGOrSATMhY1f3WRN998kyFDhjB06FDuvPPO/3j+73//\nO8OGDSM2NpbMTNe97tm47XILlyI2IoBls8by9pbD/L9P93DNC+nMdlzO/Un98PH0sDqeiIhI+/HJ\no3As++zH1J2G4zlgOsGwQVgs+HQ/8/G94mDy7896ytzcXBYsWMDGjRsJCQmhtLSUJUuWfOuY6upq\nduzYQXp6Ovfccw85OTnn+64umkaszsDDZnDXmCjWzE3impgwFn/+FZNfzGDj/pNWRxMREelYasub\nSxU0f6299F1Q1q1bx0033URISAgAQUFB/3HMrbc278o3fvx4Tp8+TVlZ2SW/7rloxOocwrr78tJt\nw/nRlcU88WEOty3dwo3DI5g3ZRDB3XysjiciImKtc4wsAc2X/964AZrqwcMbfvga2Ee2ejTDMM76\nuDVoxOo8JV/Rk9UPJfFzR38++vIIExal8W5mPk6nJreLiIiclX0k3L0CJsxr/uqCUjVhwgTee+89\nSkpKACgtLf2PY/7xj38AzRsyBwQEEBAQcMmvey4asboAft4e/OragUwfFsG8ZTk8+q9s3s8qZOGM\nOK7o5W91PBERkfbLPtKlo1SDBw9m3rx5JCUl4eHhQXx8PFFRUd86xtfXl/j4eBoaGvjb3/7mstc+\nG8Oq5QQSEhJMK9aXcBXTNHkvq5DfrdpFRW0j9yb248GJA/Dz1uR2ERFxb7t27WLQoEFWx2g13/f+\nDMPIMk0z4Vw/q0uBF8kwDH6cYGft3GSmx0fwatp+Jj2fRsruYqujiYiIiEVUrC5RUFdvnrtpKO/e\nNxpfLw9mvr6VWW9ncay81upoIiIi0sZUrFxkdL9gVs1J5JfXRLN2VzFXL07jfzYcpEmT20VERDoN\nFSsX8va0MXvCAFY/PJ7hfXvw1Ed5TH95AzsLW3/dDBEREbGeilUr6BvclTdmjuCPt8Zz7HQt01/e\nwPwVuVTUamNnERERd6Zi1UoMw2Dq0HDWzk3ijtF9eWPTISYuSuPjnUe1sbOIiIibUrFqZd19vXh6\nWizLZo0lpJsPP//f7cx8fSsFpdVWRxMREXEr8+fP57nnnrvo511BxaqNDLMHsmL2WB6/fhCZB0uZ\n9Hwaf0rdR32j0+poIiIi4iIqVm3I08PGvYn9WPNIEknRofzh0z384I8ZbD30n8vwi4iIuJMdxTt4\nLfs1dhTvcNk533zzTYYMGcLQoUO58847v/XckiVLiImJYciQIdxyyy0t3//yyy8ZM2YMAwYMYOnS\npS7L8m/a0sYC4YF+/PnOBNbkHefJFbnc9Oombk6w8+jkgfTo6m11PBERkfP2bOaz7C7dfdZjKusr\n2XNqDyYmBgZX9LiCbt7dznj8wKCB/Hrkr896ztzcXBYsWMDGjRsJCQmhtLSUJUuWtDz/+9//noMH\nD+Lj40NZ2f/dnb9z5042b95MVVUV8fHxXH/99YSHh5/nuz03jVhZ6OqYMD5/ZDz3j+/H+9sLmbg4\njQ+yCjW5XURE3EpFQwUmzf+2mZhUNFRc8jnXrVvHTTfdREhICABBQUHfen7IkCHcfvvtvPXWW3h6\n/t840rRp0/Dz8yMkJASHw0FmZuYlZ/kmjVhZrIu3J7+ZMojp8RHMW5bN3Pe+5L2sAhZMj+Pynmdu\n8yIiIu3BuUaWoPky4E9X/5QGZwNeNi9+n/h7hvUc1qq5Pv74Y9LT0/noo49YuHAh2dnZQPNd+9/0\n3ceXSiNW7cSg3t15/4GreGZGHHlHTjP5xXQWr95DbUOT1dFEREQuybCew1h6zVJmx89m6TVLXVKq\nJkyYwHvvvUdJSQkApaX/N1/Z6XRSUFCAw+Hg2Wefpby8nMrKSgA+/PBDamtrKSkpITU1lREjRlxy\nlm/SiFU7YrMZ3DaqD5Niwnhm1S6WrNvHh18eYcH0WBIHhFodT0RE5KIN6znMpaNUgwcPZt68eSQl\nJeHh4UF8fDxRUVEANDU1cccdd1BeXo5pmsyZM4fAwECg+RKhw+Hg5MmTPPHEEy6dXwVgWDWfJyEh\nwdy2bZslr91RbNh3kseX53DwZBU3DA3n8R8Moqe/r9WxRESkk9u1axeDBg2yOkar+b73ZxhGlmma\nCef6WV0KbMfGXh7CJw8m8uDEAXyac4yJi9L4++bDOLWxs4iISLukYtXO+Xp58PCkaD59KJG4iACe\nWJ7Dja9sJPdIudXRRERE5DtUrDqIfqHdePveUTx/81AKSqu54aUNLFiZR1Vdo9XRRESkE3LXpYEu\n9X2pWHUghmEwIz6StXOT+HFCJK+tP8ikxWmszj1mdTQREelEfH19KSkpcbtyZZomJSUl+Ppe/Hxm\nTV7vwLIOl/LYv3LYc7yCSTFhzL9hMBGBflbHEhERN9fQ0EBhYSG1tbVWR3E5X19fIiMj8fLy+tb3\nz3fyuopVB9fQ5OSv6w/ywpqvsBkGD18dzcyxUXh6aDBSRETEVXRXYCfh5WHjgaT+fP5wEmP6BbNw\n1S6mvrSB7fmnrI4mIiLS6ahYuQl7UBdeuzuBV+8Yzqmqen74ykbmLcumvKbB6mgiIiKdhoqVGzEM\ng+tie7NmbhIzr7qMdzLzmbgojQ93FLndBEMREZH2SMXKDXXz8eS/p8awYvY4wgN9efDdHdz1t0wO\nnayyOpqIiIhbU7FyY7ERASybNZanpw1mR34Z17yQzpK1e6lr1MbOIiIirUHFys152AzuGhPFmrlJ\nTIoJY/HnXzH5xQw27j9pdTQRERG3o2LVSYR19+Xl24bz+swRNDQ5uW3pFh755w5KKuusjiYiIuI2\nVKw6meQrerL6oSR+7ujPR18eYcKiNN7NzNfGziIiIi6gYtUJ+Xl78KtrB7JqTiJXhPnz6L+yufkv\nm/jqeIXV0URERDo0FatObECYP/+4fzR/+NEQ9hVXMuXFDJ79dDc19ZrcLiIicjFUrDo5wzD4cYKd\ntXOTmR4fwSup+5n0fBopu4utjiYiItLhqFgJAEFdvXnupqG8e99ofL08mPn6Vma9ncWxcvfbYFNE\nRKS1qFjJt4zuF8yqOYn88ppo1u4q5urFafzPhoM0aXK7iIjIOalYyX/w9rQxe8IAVj88nuF9e/DU\nR3lMf3kDOwvLrI4mIiLSrqlYyRn1De7KGzNH8Mdb4zl2upbpL29g/opcKmq1sbOIiMj3UbGSszIM\ng6lDw1k7N4k7RvfljU2HmLgojY93HtXGziIiIt+hYiXnpbuvF09Pi2XZrLGEdPPh5/+7nZmvb6Wg\ntNrqaCIiIu2GipVckGH2QFbMHsvj1w8i82Apk55P40+p+2hoclodTURExHIqVnLBPD1s3JvYjzWP\nJJEUHcofPt3D9Usy2Hqo1OpoIiIillKxkosWHujHn+9M4LW7Eqiqa+KmVzfx6/d3cqqq3upoIiIi\nllCxkkt2dUwYnz8ynvvH9+P97YVMXJzGB1mFmtwuIiKdjoqVuEQXb09+M2UQK38xjqjgLsx970tu\nXbqZfcWVVkcTERFpMypW4lKDenfn/Qeu4pkZceQdOc3kF9NZvHoPtQ3a2FlERNyfipW4nM1mcNuo\nPqydm8wPhoSzZN0+rn0hnYy9J6yOJiIi0qpUrKTVhPr78PzNw3j73lHYDIM7/5rJnHe+oLhCGzuL\niIh7UrGSVjf28hA+eTCRBycO4NOcY0xclMbfNx/GqY2dRUTEzahYSZvw9fLg4UnRfPJQInERATyx\nPIcbX9lI3pHTVkcTERFxGRUraVP9Q7vx9r2jeP7moRSUVjP1pfUs/DiPqrpGq6OJiIhcMhUraXOG\nYTAjPpK1c5P4cUIkSzMOMmlxGqtzj1kdTURE5JKoWIllArt487sbh/DBz8bg7+vFfX/P4qdvbqOo\nrMbqaCIiIhdFxUosd2XfIFbOGcejkweSsfcEkxansTT9AI3a2FlERDqY8ypWhmFcZxjGHsMw9hmG\n8ej3PN/HMIwUwzC+MAxjp2EYU1wfVdyZl4eNB5L68/nDSYzpF8zCVbuY+tIGtuefsjqaiIjIeTtn\nsTIMwwN4GZgMxAC3GoYR853DHgf+aZpmPHAL8CdXB5XOwR7UhdfuTuDVO4ZzqqqeH76ykXnLsimv\nabA6moiIyDmdz4jVSGCfaZoHTNOsB94Fpn3nGBPo/vXvA4AjrosonY1hGFwX25s1c5OYedVlvJOZ\nz8RFaXy4o0gbO4uISLt2PsUqAij4xuPCr7/3TfOBOwzDKARWAb9wSTrp1Lr5ePLfU2NYMXsc4YG+\nPPjuDu76WyaHTlZZHU1EROR7uWry+q3A66ZpRgJTgL8bhvEf5zYM4z7DMLYZhrHtxAntGyfnJzYi\ngGWzxvL0tMF8kV/GNS+ks2TtXuoatbGziIi0L+dTrIoA+zceR379vW/6CfBPANM0NwG+QMh3T2Sa\n5l9M00wwTTMhNDT04hJLp+RhM7hrTBRr5yYxKSaMxZ9/xeQXM9i0v8TqaCIiIi3Op1htBQYYhnGZ\nYRjeNE9OX/GdY/KBiQCGYQyiuVhpSEpcLqy7Ly/fNpzXZ46gocnJrUs388g/d1BSWWd1NBERkXMX\nK9M0G4HZwGfALprv/ss1DONpwzBu+PqwucBPDcP4EngH+C9Ts4ylFSVf0ZPVDyXxc0d/PvryCBMX\np/GPrfna2FlERCxlWNV/EhISzG3btlny2uJe9h6vYN6yHDIPlTIiqgcLZ8QRHeZvdSwREXEjhmFk\nmaaZcK7jtPK6dHgDwvz5x/2j+cOPhrCvuJIpL2bw7Ke7qanX5HYREWlbKlbiFgzD4McJdtbOTWZ6\nfASvpO5n0vNppOwutjqaiIh0IipW4laCunrz3E1Defe+0fh6eTDz9a3MejuLY+W1VkcTEZFOQMVK\n3NLofsGsmpPIL6+JZu2uYq5enMb/bDhIkya3i4hIK1KxErfl7Wlj9oQBrH54PMP79uCpj/KY/vIG\ndhaWWR1NRETclIqVuL2+wV15Y+YI/nhrPMdO1zL95Q3MX5FLRa02dhYREddSsZJOwTAMpg4NZ80j\nSdw+qi9vbDrE1YvTWJV9VBs7i4iIy6hYSacS4OfFb6fHsmzWWIK7+jDr7e3c8/pWCkqrrY4mIiJu\nQMVKOqVh9kBWzB7L49cPYsvBUiY9n8afUvfR0OS0OpqIiHRgKlbSaXl62Lg3sR9rHkkiKTqUP3y6\nh+uXZLD1UKnV0UREpINSsZJOLzzQjz/fmcBrdyVQVdfETa9u4tfv7+RUVb3V0UREpINRsRL52tUx\nYXz+yHjuH9+P97cXMnFxGh9kFWpyu4iInDcVK5Fv6OLtyW+mDGLlL8YRFdyFue99ya1LN7OvuNLq\naCIi0gGoWIl8j0G9u/P+A1fxzIw48o6cZvKL6SxevYfaBm3sLCIiZ6ZiJXIGNpvBbaP6sHZuMj8Y\nEs6Sdfu49oV0MvaesDqaiIi0UypWIucQ6u/D8zcP462fjMJmGNz510zmvPMFxRXa2FlERL5NxUrk\nPI0bEMInDyby4MQBfJpzjImL0nhr82Gc2thZRES+pmIlcgF8vTx4eFI0nzyUSFxEAI8vz+HGVzaS\nd+S01dFERKQdULESuQj9Q7vx9r2jeP7moRSUVjP1pfUs/DiPqrpGq6OJiIiFVKxELpJhGMyIj2Tt\n3CR+nBDJ0oyDTFqcxurcY1ZHExERi6hYiVyiwC7e/O7GIXzwszH4+3px39+z+Omb2ygqq7E6moiI\ntDEVKxEXubJvECvnjOPRyQPJ2HuCSYvTWJp+gEZt7Cwi0mmoWIm4kJeHjQeS+vP5w0mM6RfMwlW7\nmPrSBrbnn7I6moiItAEVK5FWYA/qwmt3J/DqHcM5VVXPD1/ZyLxl2ZTXNFgdTUREWpGKlUgrMQyD\n62J7s2ZuEjOvuox3MvOZuCiND3cUaWNnERE3pWIl0sq6+Xjy31NjWDF7HOGBvjz47g7u+lsmh05W\nWR1NRERcTMVKpI3ERgSwbNZYnrphMF/kl3HNC+ksWbuXukZt7Cwi4i5UrETakIfN4O6rolg7N4lJ\nMWEs/vwrJr+Ywab9JVZHExERF1CxErFAWHdfXr5tOK/PHEFDk5Nbl27mkX/uoKSyzupoIiJyCVSs\nRCyUfEVPVj+UxM8d/fnoyyNMXJzGP7bma2NnEZEOSsVKxGJ+3h786tqBrJqTSHRPf379QTY3/2UT\nXx2vsDqaiIhcIBUrkXZiQJg//7h/NH/40RD2FVcy5cUMnv10NzX1mtwuItJRqFiJtCOGYfDjBDtr\n5yYzPT6CV1L3M+n5NFJ2F1sdTUREzoOKlUg7FNTVm+duGsq7943G18uDma9vZdbbWRwrr7U6moiI\nnIWKlUg7NrpfMKvmJPLLa6JZu6uYqxen8T8bDtKkye0iIu2SipVIO+ftaWP2hAGsfng8w/v24KmP\n8pj+8gayC8utjiYiIt+hYiXSQfQN7sobM0fwx1vjOXa6lmkvr2f+ilwqarWxs4hIe6FiJdKBGIbB\n1KHhrHkkidtH9eWNTYe4enEaq7KPamNnEZF2QMVKpAMK8PPit9NjWTZrLMFdfZj19nbueX0rBaXV\nVkcTEenUVKxEOrBh9kBWzB7L49cPYsvBUiY9n8afUvfR0OS0OpqISKekYiXSwXl62Lg3sR9rHkki\nKTqUP3y6h+uXZLD1UKnV0UREOh0VKxE3ER7ox5/vTOC1uxKoqmviplc38ev3d3Kqqt7qaCIinYaK\nlYibuTomjM8fGc/94/vx/vZCJi5O44OsQk1uFxFpAypWIm6oi7cnv5kyiJW/GEdUcBfmvvclty7d\nzL7iSqujiYi4NRUrETc2qHd33n/gKp6ZEUfekdNMfjGdxav3UNugjZ1FRFqDipWIm7PZDG4b1Ye1\nc5P5wZBwlqzbx3UvpJOx94TV0URE3I6KlUgnEervw/M3D+Otn4zCMAzu/Gsmc975guIKbewsIuIq\nKlYincy4ASF88mAiD04cwKc5x5i4KI23Nh/GqY2dRUQumYqVSCfk6+XBw5Oi+eShROIiAnh8eQ43\nvrKRvCOnrY4mItKhqViJdGL9Q7vx9r2jeP7moRSUVjP1pfUs/DiPqrpGq6OJiHRIKlYinZxhGMyI\nj2Tt3CR+nBDJ0oyDTFqcxurcY1ZHExHpcFSsRASAwC7e/O7GIXzwszH4+3px39+z+Omb2ygqq7E6\nmohIh6FiJSLfcmXfIFbOGcejkweSsfcEkxansTT9AI3a2FlE5JxUrETkP3h52HggqT+fP5zEmH7B\nLFy1i6kvbeCL/FNWRxMRaddUrETkjOxBXXjt7gRevWM4p6rqufGVjTy+PJvymgaro4mItEsqViJy\nVoZhcF1sb9bMTWLmVZfxv1vymbgojQ93FGljZxGR71CxEpHz0s3Hk/+eGsOK2eMID/TlwXd3cNff\nMjl0ssrqaCIi7YaKlYhckNiIAJbNGstTNwzmi/wyrnkhnSVr91LXqI2dRURUrETkgnnYDO6+Koq1\nc5OYFBPG4s+/YvKLGWzaX2J1NBERS6lYichFC+vuy8u3Def1mSNoaHJy69LNPPLPHZRU1lkdTUTE\nEipWInLJkq/oyeqHkvi5oz8ffXmEiYvT+MfWfG3sLCKdjoqViLiEn7cHv7p2IKvmJBLd059ff5DN\nzX/ZxFfHK6yOJiLSZlSsRMSlBoT584/7R/OHHw1hX3ElU17M4NlPd1NTr8ntIuL+VKxExOUMw+DH\nCXbWzk1menwEr6TuZ9LzaaTsLrY6mohIq1KxEpFWE9TVm+duGsq7943G18uDma9vZdbbWRw/XWt1\nNBGRVqFiJSKtbnS/YFbNSeSX10SzdlcxExel8fqGgzRpcruIuBkVKxFpE96eNmZPGMDqh8cT3yeQ\n+R/lMf3lDWQXllsdTUTEZc6rWBmGcZ1hGHsMw9hnGMajZzjmx4Zh5BmGkWsYxv+6NqaIuIu+wV15\n856R/PHWeI6drmXay+uZvyKXilpt7CwiHZ/nuQ4wDMMDeBmYBBQCWw3DWGGaZt43jhkA/AYYa5rm\nKcMwerZWYBHp+AzDYOrQcMZHh/LcZ3t4Y9MhPsk5ypNTBzM5theGYVgdUUTkopzPiNVIYJ9pmgdM\n06wH3gWmfeeYnwIvm6Z5CsA0Td36IyLnFODnxW+nx7Js1liCu/ow6+3t3PP6VgpKq62OJiJyUc6n\nWEUABd94XPj1974pGog2DGODYRibDcO47vtOZBjGfYZhbDMMY9uJEycuLrGIuJ1h9kBWzB7L49cP\nYsvBUiY9n8afUvfR0OS0OpqIyAVx1eR1T2AAkAzcCiw1DCPwuweZpvkX0zQTTNNMCA0NddFLi4g7\n8PSwcW9iP9Y8kkRSdCh/+HQP1y/JYOuhUqujiYict/MpVkWA/RuPI7/+3jcVAitM02wwTfMg8BXN\nRUtE5IKEB/rx5zsTeO2uBKrqmrjp1U38+v2dlFXXWx1NROSczqdYbQUGGIZxmWEY3sAtwIrvHLOc\n5tEqDMMIofnS4AEX5hSRTubqmDA+f2Q894/vx/vbC5mwKI0PsgoxTa19JSLt1zmLlWmajcBs4DNg\nF/BP0zRzDcN42jCMG74+7DOgxDCMPCAF+JVpmiWtFVpEOocu3p78ZsogVv5iHFHBXZj73pfctnQL\n+09UWh1NROR7GVZ9+ktISDC3bdtmyWuLSMfjdJq8u7WA33+yi9oGJw8k9WOW43J8vTysjiYinYBh\nGFmmaSac6zitvC4iHYLNZnDbqD6snZvMlLheLFm3j+teSCdjr+4wFpH2Q8VKRDqUUH8fXrglnrd+\nMgrDMLjzr5nMeecLiiu0sbOIWE/FSkQ6pHEDQvjkwUQenDiAT3OOMXFRGm9tPoxTGzuLiIVUrESk\nw/L18uDhSdF88lAicREBPL48hxtf2UjekdNWRxORTkrFSkQ6vP6h3Xj73lE8f/NQCkqrmfrSehZ+\nnEdVXaPV0USkk1GxEhG3YBgGM+IjWTs3iR8nRLI04yCTFqexOveY1dFEpBNRsRIRtxLYxZvf3TiE\nD342Bn9fL+77exY/fXMbR8pqrI4mIp2AipWIuKUr+waxcs44Hp08kIy9J7h6cRqvZRygURs7i0gr\nUrESEbf8I0EHAAAgAElEQVTl5WHjgaT+fP5wEmP6BbPg411MfWkDX+SfsjqaiLgpFSsRcXv2oC68\ndncCr94xnFNV9dz4ykYeX55NeU2D1dFExM2oWIlIp2AYBtfF9mbN3CT+66oo/ndLPhMXpfHhjiJt\n7CwiLqNiJSKdSjcfT56cOpgVs8cRHujLg+/u4K6/ZXLoZJXV0UTEDahYiUinFBsRwLJZY3nqhsF8\nkV/GNS+ks2TtXuoam6yOJiIdmIqViHRaHjaDu6+KYu3cJCbFhLH486+Y/GIGm/aXWB1NRDooFSsR\n6fTCuvvy8m3DeX3mCBqanNy6dDOP/HMHJZV1VkcTkQ5GxUpE5GvJV/Rk9UNJ/NzRn4++PMLExWn8\nY2u+NnYWkfOmYiUi8g1+3h786tqBrJqTSHRPf379QTY3/2UTXx2vsDqaiHQAKlYiIt9jQJg//7h/\nNH/40RD2FVcy5cUMnv10NzX1mtwuImemYiUicgaGYfDjBDtr5yYzPT6CV1L3c80LaaTsKbY6moi0\nUypWIiLnENTVm+duGsq7943Gx9ODmf+zlVlvZ3H8dK3V0USknVGxEhE5T6P7BbNqTiK/vCaatbuK\nmbgojdc3HKRJk9tF5GsqViIiF8Db08bsCQNY/fB44vsEMv+jPKa/vIHswnKro4lIO6BiJSJyEfoG\nd+XNe0byx1vjOXa6lmkvr2f+ilwqarWxs0hnpmIlInKRDMNg6tBw1jySxO2j+vLGpkNcvTiNVdlH\ntbGzSCelYiUicokC/Lz47fRYls0aS3BXH2a9vZ17Xt9KQWm11dFEpI2pWImIuMgweyArZo/l8esH\nseVgKZOeT+NPqftoaHJaHU1E2oiKlYiIC3l62Lg3sR9rHkkiKTqUP3y6h+uXZLDtUKnV0USkDahY\niYi0gvBAP/58ZwKv3ZVAVV0TP3p1E49+sJOy6nqro4lIK1KxEhFpRVfHhPH5I+O5f3w/3ssqZMKi\nND7IKtTkdhE3pWIlItLKunh78pspg1j5i3FEBXdh7ntfctvSLew/UWl1NBFxMRUrEZE2Mqh3d95/\n4CqemRFH7pFyJr+QweLVe6ht0MbOIu5CxUpEpA3ZbAa3jerD2rnJTInrxZJ1+7juhXQy9p6wOpqI\nuICKlYiIBUL9fXjhlnje+skoDMPgzr9mMuedLyiu0MbOIh2Z+xargkzIWNT8VUSknRo3IIRPHkzk\nwYkD+DTnGBMXpfHW5sM4tbGzSIdkWHVnSkJCgrlt27bWOXlBJrzxA2isB5sHJPwEekSBYfv6l/GN\n35/pl3Eex1/oeS7iGIwLPI/ROv+bikir23+ikieW57BxfwnD7IE8MyOOmPDuVscSEcAwjCzTNBPO\ndZxnW4Rpc4cy2OFhsK2rPwm1tQzL/LPVidrQuQrYhTz/fce6okye43mMSz/HRee4mPfsooINrjmP\nCnaH1T+0G2/fO4rlO4pYsHIXU19azz1jo3jo6mi6+rjnX9ci7sYt/0vd0aM3P+kVSr0BXgTwx9hZ\njB18G5hOMM2vv57t1xmOwTz3Medzngs+h6vO09rv+evnna7K8X3HnuPn0eWTZhdaBs/1uJWKrUs/\nCFzoMeczKt1KWc7yvg3DYEZfGxPvvYw/pR1kecY2tnyZy9xrB5F0RZgKtkg755bFaptZTYPNBpg0\nAD/LfZURJVk47A6S7clE+kdaHVFay3kV0fMsaC4prK4ok64u6md7z654nYsp2Od63PkKdnfgUeBR\nX6Ae+OjrXxektQu0K8pkW5Ta8xkFb6UsF/2+XZnFzQt2QSYcyoCoRLCPtDqNexarhLAEfDx8aHA2\n4GF4cG3UteSV5PHs1md5duuzDOgxgOTIZCb0mUBMcAw2w2Z1ZHGVlk/p+v+0U7uQgk0rF9YLKo/f\nf55GZxMb954gZfcxPDCZODCEkVE98Pje7OeRBVzwni62YF/EKHgnLtgu0+rF9VxlspWKZGUx7Pm4\n+c+Zpw/cvcLycuWek9eBHcU72HZ8GwlhCQzrOQyA/NP5pBakklKQwvbi7ThNJ6F+oSTZk3DYHYzq\nPQofD59WyyQicikKSquZvyKXtbuLGdS7O8/MiCW+Tw+rY7Vf3y3YZyzRbTE6fb4jwq46xjzHe3bx\n61zSe76E5xuqofHrJUoMD5gwDxLntsofp/OdvO62xepcymrLyCjKIKUghQ1FG6hurMbP04+x4WNJ\nticzPnI8PXz1F5aItC+mafJZ7jHmr8jjeEUtt4/qw6+uHUiAn5fV0UTaXkEmvHEDNNWDh3erjlip\nWF2A+qZ6Mo9ltoxmFVcXYzNsDAsdhsPuwNHHQd/ufa2OKSLSorKukUWr9/DGxkMEdfXhiR8M4oah\n4RjuPp9G5LvaaI6VitVFMk2TvNI8UvJTSC1IZc+pPQBcFnBZc8myO4gLicPD5mFxUhERyCkq57Fl\n2ewsLCdxQAi/nRZLVEhXq2OJuB0VKxc5UnmElILmkrXt2DYazUaCfINIikwi2Z7MmPAx+Hn6WR1T\nRDqxJqfJW5sP8/8+20N9k5PZjsu5P6kfPp76ACjiKipWreB0/Wk2FG0gJT+FjKIMKhsq8fHwYUzv\nMTj6OBgfOZ4QvxCrY4pIJ3X8dC1Pr8zj451H6RfalYXT4xjTP9jqWCJuQcWqlTU0NZBVnNVyyfBI\n1REMDOJC41ouGfYL6Kf5DiLS5lL3FPPEhzkUlNbww+GRPDZlIMHddMezyKVQsWpDpmny1amvSClI\nIaUghbySPAD6+Pch2Z6Mw+5gWM9heNrcctkwEWmHauqbeCllL39JP0BXH09+M3kgN11px2bThz2R\ni6FiZaFjVcdIL0xnXcE6Mo9m0uBsINAnkPGR40m2J3NV+FV09dLkUhFpfXuPVzBvWQ6Zh0oZEdWD\nhTPiiA7ztzqWSIejYtVOVDVUsaFoA6kFqaQXpVNeV46XzYtRvUfhsDtIikwirGuY1TFFxI2Zpsl7\nWYX8btUuKmob+en4fsyZMAA/b01uFzlfKlbtUKOzkS+Kv2hZL6ugogCAwcGDWy4ZRveI1rwsEWkV\npVX1PLNqF+9nFWIP8uPpabE4ruhpdSyRDkHFqp0zTZMD5Qda5mVln8jGxCSiWwTJ9mSS7clcGXYl\nXjatpiwirrX5QAnzlmWz/0QVU+J68eTUwYR197U6lki7pmLVwZysOUlaQRqpBalsOrqJuqY6/L38\nGRc5DofdwbiIcfh7a16EiLhGXWMTS9MP8Md1+/DysPHLa6K5c0wUHprcLvK9VKw6sOqGajYf3UxK\nQQrphemU1pbiaXiS0CsBh91Bsj2Z8G7hVscUETdwuKSKx5fnkLH3JHERATwzI464yACrY4m0OypW\nbqLJ2UT2yWzWFawjtSCVg+UHARgYNLDlkmFMUIzmZYnIRTNNk5U7j/L0yjxKKuu4a0wUc6+Jxt9X\nUxFE/k3Fyk0dKj/UMvl9x4kdOE0nPbv0bFmUdESvEXh7eFsdU0Q6oPKaBp77bA9vbTlMT38fnpw6\nmMmxvfTBTQQVq06htLaUjMIMUgpS2HhkIzWNNXTx7MLYiLE47M1b7AT4aEhfRC7MjoIyHvtXNnlH\nT+O4IpSnp8ViD+pidSwRS6lYdTJ1TXVsObqlZcPokzUn8TA8GB42nOTIZBx9HNj97VbHFJEOorHJ\nyesbD7H4869wmiYPTozm3sTL8PKwWR1NxBIqVp2Y03SSezK3ZSmHfWX7ALg88PKW9bJiQ2KxGfoL\nUkTO7khZDU99lMtnuceJDuvGMzPiSIgKsjqWSJtTsZIWBRUFpBakklqQStbxLJrMJkL8QkiKTMJh\ndzCq9yh8PbWGjYic2Zq84zy5IpeishpuGWHn0ckDCeyi+ZzSeahYyfcqrysnoyiD1IJU1hetp6qh\nCj9PP8b0HkOyPZkkexJBvvo0KiL/qbq+kRfX7OW19QcJ8PNi3pRB3Dg8QpPbpVNQsZJzqm+qZ9ux\nbS1LORyvPo6BwbCew1rWy7os4DKrY4pIO7Pr6GnmLctme34ZY/oFs2BGLP1Du1kdS6RVqVjJBTFN\nk92lu1smv+8q3QVAVPeolpI1NHQoHjZt2ioi4HSavLu1gN9/sovaBicPJPVjluNyfL30d4S4JxUr\nuSRHK4+SWphKSn4KW49vpdHZSA+fHoyPHI+jj4MxvcfQxUu3X4t0dicq6lj4cR7LdxwhKrgLv50e\nS+KAUKtjibicipW4TEV9BRuObCAlP4WMogwq6ivwtnkzOnw0DruDpMgkQrvoL1KRzmz93pM88WEO\nB09WccPQcB7/wSB6+uumGHEfKlbSKhqcDWw/vr1l9feiyiIAhoQMadli5/LAyzWZVaQTqm1o4pXU\n/bySuh8fLxu/vm4gt43sg00bO4sbULGSVmeaJnvL9jaXrPwUckpyAIjsFkmyPZkJfSYQ3zMeT5un\nxUlFpC3tP1HJE8tz2Li/hPg+gSycHkdMeHerY4lcEhUraXPF1cUt62VtObqFemc93b27Mz5yPMn2\nZMaGj6Wbt+4cEukMTNNk+Y4iFqzcRVlNA/eMjeKhq6Pp6qMPWtIxubRYGYZxHfAi4AG8Zprm789w\n3A+B94ERpmmetTWpWLm36oZqNh7ZSEpBCumF6ZTVleFl82Jkr5Etlwx7de1ldUwRaWVl1fU8++lu\n3sksIDzAl6emxTIpJszqWCIXzGXFyjAMD+ArYBJQCGwFbjVNM+87x/kDHwPewGwVK/m3RmcjX574\nkpT85i128ivyARgUNAhHHwcOu4MrelyheVkibizrcCmP/SuHPccrmBQTxlM3DCY80M/qWCLnzZXF\nagww3zTNa79+/BsA0zR/953jXgA+B34F/FLFSr6PaZocPH2QlPzm9bK+PPElJia9u/ZuGckaETYC\nLw8vq6OKiIs1NDn56/qDvLDmK2yGwSOTovmvq6Lw1MbO0gG4slj9CLjONM17v358JzDKNM3Z3zhm\nODDPNM0fGoaRyhmKlWEY9wH3AfTp0+fKw4cPX8BbEnd0suYkGYUZpBSksOnIJmqbaunm1Y1xEeNw\n2B2MixxHd29NehVxJwWl1cxfkcva3cUM6t2dZ2bEEt+nh9WxRM6qzYqVYRg2YB3wX6ZpHjpbsfom\njVjJd9U01rDl6JaW1d9La0vxNDy5MuxKHH2aV3+P6BZhdUwRcQHTNPks9xjzV+RxvKKW20f14VfX\nDiTAT6PV0j612aVAwzACgP1A5dc/0gsoBW44W7lSsZKzcZpOdp7Y2XKX4f7y/QBE94huXsrBPoFB\nwYOwGbqEINKRVdY1smj1Ht7YeIigrj488YNB3DA0XHMupd1xZbHypHny+kSgiObJ67eZppl7huNT\n0YiVuNjh04dbFiX9ovgLnKaTnn49SbIn4bA7GNl7JD4ePlbHFJGLlFNUzmPLstlZWE7igBB+Oy2W\nqJCuVscSaeHq5RamAC/QvNzC30zTXGgYxtPANtM0V3zn2FRUrKQVldWWkV6UTmpBKuuL1lPTWIOf\npx/jIsaRbE9mfMR4An0DrY4pIheoyWny1ubD/L/P9lDf5OQXjsu5L6kfPp7a2FmspwVCpVOoa6oj\n82hmyyXD4ppibIaN+J7xOOzNSzn06d7H6pgicgGOn67l6ZV5fLzzKP1Du7Jgehxj+gdbHUs6ORUr\n6XScppNdJbtIKWheL+urU18B0C+gHw578+T3IaFDNC9LpINI3VPMEx/mUFBaww+HR/LYlIEEd9Ml\nf7GGipV0ekWVRS3zsrKOZdFoNhLkG9S8XlZkMqPDR+PnqQUKRdqzmvomXkrZy1/SD9DVx5PfTB7I\nTVfatbGztDkVK5FvKK8rZ0PRBlIKUlhftJ7Khkp8PXwZHT6aCfYJJEYmEuIXYnVMETmDvccrmLcs\nh8xDpYyI6sHCGXFEh/lbHUs6ERUrkTNoaGpg2/FtLetlHa06ioHBkNAhLfOyLgu4TLd7i7Qzpmny\nXlYhv1u1i4raRn46vh9zJgzAz1uT26X1qViJnAfTNNlzak/zvKz8FHaV7gKgb/e+JEcm4+jjYGjo\nUDxtnhYnFZF/K62q55lVu3g/qxB7kB9PT4vFcUVPq2OJm1OxErkIx6qOkVaQRkpBCluObaHR2Uig\nTyDjI8fjsDu4Kvwqunh1sTqmiACbD5Qwb1k2+09UMSWuF09OHUxYd1+rY4mbUrESuUSV9ZVsOLKB\n1IJU0gvTOV1/Gm+bN6N6j2rZMLpnF31KFrFSXWMTS9MP8Md1+/DysPHLa6K5c0wUHprcLi6mYiXi\nQo3ORr4o/qLlkmFhZSEAscGxJNubLxkOCBygeVkiFjlcUsXjy3PI2HuSuIgAnpkRR1xkgNWxxI2o\nWIm0EtM02V+2v2Xy+86TOwGI6BbRXLLsDoaHDcfLps1kRdqSaZqs3HmUp1fmUVJZx11joph7TTT+\nvvpvUS6dipVIGzlRfYK0wjRSC1LZfHQzdU11+Hv7kxiRiMPuYGzEWPy9dVu4SFspr2nguc/28NaW\nw/T092H+1MFcF9tLI8pySVSsRCxQ3VDNpqObSC1IJa0gjVN1p/C0eTIibETLaFbvbr2tjinSKewo\nKOOxf2WTd/Q0jitCeXpaLPYg3XwiF0fFSsRiTc4mdp7cSUp+8xY7h04fAmBg0MCWLXYGBQ3Sp2iR\nVtTY5OT1jYdY/PlXOE2TBydGc2/iZXh5aGsruTAqViLtzMHygy2bRX9R/AUmJmFdwlpGskb0GoG3\nh7fVMUXc0pGyGp76KJfPco8THdaNZ2bEkRAVZHUs6UBUrETasdLaUtIL00nJT2HT0U3UNNbQ1asr\nY8PH4ujjIDEikQAf3dEk4mpr8o7z5IpcispquGWEnUcnDySwiz7QyLmpWIl0ELWNtWQey2Rd/jrS\nCtM4WXMSD8OD4WHDWy4Z2v3tVscUcRvV9Y28uGYvr60/SICfF/OmDOLG4RG6LC9npWIl0gE5TSc5\nJ3NILUglpSCFfWX7ALg88PKWfQwHhwzGZmh+iMil2nX0NI8ty+aL/DLG9AtmwYxY+od2szqWtFMq\nViJuoOB0AamFzSVr+/HtNJlNhPqFkmRPwmF3MLLXSHw9tYWHyMVyOk3e2ZrPs5/sprbByQNJ/Zjl\nuBxfL23sLN+mYiXiZsrrykkvTCe1IJX1ReupbqzGz9OPq8KvItmeTFJkEj18e1gdU6RDOlFRx8KP\n81i+4whRwV1YMD2OcQNCrI4l7YiKlYgbq2+qZ+uxrS2rvx+vPo7NsDEsdFjLvKyogCirY4p0OOv3\nnuSJD3M4eLKKacPCefz6GEL9fayOJe2AipVIJ2GaJrtKd7WUrN2luwGI6h6Fo0/zvKwhIUPwsOnS\nhsj5qG1o4pXU/bySuh8fLxu/vm4gt43sg00bO3dqKlYindTRyqMtJWvrsa00mo0E+QYxPnI8yfZk\nxvQeQxcvrT4tci77T1TyxPIcNu4vIb5PIAunxxET3t3qWGIRFSsRoaK+gg1FG1hXsI71heupaKjA\nx8OH0b1H47A7SLInEeKneSQiZ2KaJst3FLFg5S7Kahq4Z2wUD10dTVcfT6ujSRtTsRKRb2lwNrD9\n+HZSClJIyU/hSNURDAziQuOa52VFJtM/sL/W8hH5HmXV9Tz76W7eySwgPMCXp6bFMikmzOpY0oZU\nrETkjEzT5KtTX7Wsl5VbkguA3d/essVOfM94PG36VC7yTVmHS3nsXznsOV7BpJgwnrphMOGBflbH\nkjagYiUi5+141XHSCtNIKUhhy9EtNDgbCPAJYHxE87yssRFj6erV1eqYIu1CQ5OTv64/yAtrvsJm\nGDwyKZr/uioKT23s7NZUrETkolQ1VLHxyEZSC1JJK0yjvK4cL5sXI3uPxBHZvJRDWFddAhEpKK3m\nyRW5rNtdzKDe3XlmRizxfbSWnLtSsRKRS9bobGRH8Y6WS4b5FfkAxATHkGxPZoJ9AtE9ojUvSzot\n0zT5LPcY81fkcbyilttH9eFX1w4kwM/L6mjiYipWIuJSpmlysPwg6wrWkVqQys4TOzExCe8aTrI9\nmWR7Mgm9EvCy6R8U6Xwq6xpZtHoPb2w8RHA3H574QQxTh/TWhw43omIlIq3qZM1J0gvTSSlIYfOR\nzdQ21eLv5c+4iHEk25MZFzmO7t5a80c6l5yich5bls3OwnISB4Tw22mxRIVofqI7ULESkTZT01jD\n5iObSSlIIa0wjdLaUjwNT67sdSUOe/Pq7+Hdwq2OKdImmpwmb20+zP/7bA/1TU5+4bic+5L64eOp\n3Q86MhUrEbFEk7OJ7JPZLau/Hyg/AMAVPa5oWcohJjhGl0jE7R0/XcvTK/P4eOdR+od2ZcH0OMb0\nD7Y6llwkFSsRaRcOnz5MakEq6/LXsePEDpymk55depIcmYyjj4ORvUbi7eFtdUyRVpO6p5gnPsyh\noLSGHw6P5LEpAwnupo2dOxoVKxFpd07VniKjKIOU/BQ2HNlATWMNXTy7MDZiLA67g8SIRAJ9A62O\nKeJyNfVNvJSyl7+kH6Crjye/mTyQm660a2PnDkTFSkTatbqmOrYc3UJqQSqpBamcqDmBh+FBfM/4\nlqUc7N3tVscUcam9xyuYtyyHzEOljIjqwcIZcUSH+VsdS86DipWIdBhO00leSV7zPoYFKew9tReA\n/gH9m+dl9XEQFxKHzdDK1tLxmabJe1mF/G7VLipqG/np+H7MmTAAP29Nbm/PVKxEpMMqrChsGcna\ndnwbTWYTwb7BLetlje49Gl9PX6tjilyS0qp6nlm1i/ezCrEH+fH0tFgcV/S0OpacgYqViLiF8rpy\n1hetJ7UglfVF66lsqMTXw5cx4WNw2B2MjxxPsJ/utJKOa/OBEuYty2b/iSqmxPXiyamDCeuuDw7t\njYqViLidhqYGth7fSkp+CqmFqRyrOoaBwdDQoTj6NO9j2C+gn9UxRS5YXWMTS9MP8Md1+/DysPHL\na6K5c0wUHprc3m6oWImIWzNNkz2n9pCS3zwva1fpLgCiuke1XDIcFjoMD5vmrUjHcbikiseX55Cx\n9yRDIgN4ZkYcsREBVscSVKxEpJM5VnWsZbPozGOZNDobCfQJZHzkeCbYJzAmfAxdvLpYHVPknEzT\nZOXOozy9Mo+SyjruGhPF3Gui8ffVPpxWUrESkU6rsr6SDUc2kFKQQnphOhX1FXjbvBkdPrp5NCsy\nmdAuoVbHFDmr8poGnvtsD29tOUxPfx/mTx3MdbG9tGuBRVSsRESABmcDXxz/omUph6LKIgDiQuJa\ntti5PPBy/WMl7daOgjIe+1c2eUdP47gilKenxWIP0uhrW1OxEhH5DtM02Ve2r+WSYfbJbAAiukW0\nbBYdHxaPl02XXKR9aWxy8vrGQyz+/CucpsmDE6O5N/EyvDy0tltbUbESETmHE9UnSC1sXi9r85HN\n1Dvr6e7dncTIRJLtyYwLH0c3725WxxRpcaSshqc+yuWz3ONEh3XjmRlxJEQFWR2rU1CxEhG5ANUN\n1Ww6sqllXtapulN42jwZ2WtkyyXDXl17WR1TBIA1ecd5ckUuRWU13DLCzqOTBxLYRZuZtyYVKxGR\ni9TkbOLLE1+2zMs6fPowAIOCBuGwN6+XNTBooOZliaWq6xt5Yc1e/rr+IAF+XsybMogbh0foz2Ur\nUbESEXGRA+UHWrbY2VG8AxOTXl17kRzZPJI1otcIvDw0L0ussevoaR5bls0X+WWM6RfMghmx9A/V\nJWxXU7ESEWkFJTUlpBemk1qQysYjG6ltqqWrV1fGRYzDYXcwLmIcAT5a0FHaltNp8s7WfJ79ZDe1\nDU4eSO7PrOT++HppgVxXUbESEWlltY21bDm6hZSCFFILUimpLcHD8ODKsCtbLhlG+kdaHVM6kRMV\ndSz8OI/lO44QFdyFBdPjGDcgxOpYbkHFSkSkDTlNJ9kns5uXcshPYX/5fgAG9BhAcmQyE/pMICY4\nBpuh2+Ol9a3fe5InPszh4Mkqpg0L5/HrYwj197E6VoemYiUiYqH80/kt62VtL96O03QS6hdKkj0J\nh93BqN6j8PHQP3TSemobmngldT+vpO7Hx8vGr68byG0j+2DTxs4XRcVKRKSdKKstI6Mog5SCFDYU\nbaC6sRo/Tz/Gho8l2Z7M+Mjx9PDtYXVMcVP7T1TyxPIcNu4vIb5PIAunxxET3t3qWB2OipWISDtU\n31RP5rHMltGs4upibIaNYaHDmld/7+Ogb/e+VscUN2OaJst3FLFg5S7Kahq4Z2wUD10dTVcfT6uj\ndRgqViIi7ZxpmuSV5pGS3zz5fc+pPQBcFnBZyxY7cSFxeNh0Z5e4Rll1Pc9+upt3MgsID/DlqWmx\nTIoJszpWh6BiJSLSwRypPNJyh+G2Y9toNBsJ8g0iKTKJZHsyY8LH4OfpZ3VMcQNZh0t57F857Dle\nwaSYMJ66YTDh/7+9ew+usjzwOP59Tm6EEAi3hJCcGO6KQAG5E5L3WGulVtS6O0u7dKd2HGZ3aaeO\n6+4WrbZe6tpZt1OmHTtTtGt3cFt3LFq1LVaWcxLuhJsQ7nLLSUgCISbknpOcZ/846VnbCoKc8Oac\n/D4zjjxwPOc3j2fgx/s+7/Nk6bt1JSpWIiJx7FLXJbZWb8Vf6Wdz9WZaQi2kJaWxMHchvgIfxfnF\njErXY/Ty6YV6wry85TQ/2ngcjzE88rnJfG1RIck62PljqViJiCSIUE+IPef3RG8Znms9h8EwffT0\n6C3D8cPG6ygT+VSCDW18961DbDp6nqm5Q3nuS9OZ6c1yO1a/o2IlIpKArLUc//B49BzDwxcPA1CQ\nWYDjdXC8DrOyZ5Hs0aJkuXrWWt49VMv33jpMXXMHfzu/gH/+/M0MS9dRTX+kYiUiMgDUttZSVlXG\npuAmdtXsIhQOMSxtGMV5xfgKfCwau4iMlAy3Y0qcaOns5j/+cIxfbDvDyCFpPPHFqdwzI1dXQ1Gx\nEhEZcFpDrWyt3kogGKCsuoymziZSPCnMz52Pz+ujJL+EnAw9ASafrKK6icfeOMiBqiaWTBrFM/dO\noy8EiCYAABLPSURBVHDUwC7oKlYiIgNYd7ibfef3RffLCjYHAbh15K04Xgef18fk4ZN1JUIuqyds\nWbfjLP/+7jG6esJ80zeRlSXjSUsemNt/qFiJiAgQWT9zqulUdF3WwQsHsVjyhuRF12XdlnMbKR6t\np5G/VHepg6ffOcxvD9QwYXQGz943nYUTRrod64ZTsRIRkY9V315PabCUQDDA9prtdPZ0kpmSSVF+\nET6vj6K8IjJTM92OKf1M4Nh5nvhNBcGGdh6Ync9jX7iZkUMGznmXKlYiIvKJ2kJt7KjZQSAYoLSq\nlIaOBpJNMnPGzMHn9eF4HcYOGet2TOkn2rt6+In/BD8rO0VGWjKrl97MX9/mHRAHO6tYiYjINekJ\n93Cw/iCbgpsIBAOcbjoNwJThU/AVRErW1BFTtS5LOFHXzONvVLDrTANzC4fz/funMzknsa9yqliJ\niMh1OdN0Jrr4ff+F/YRtmOzB2dFNSeeOmUtqUqrbMcUl4bDl9T1VPPf7I7R0dLOyeDzfvH0S6amJ\nubg9psXKGHMXsAZIAl6y1j7/Z7/+CPAQ0A1cAL5urT17pfdUsRIRiR8NHQ1srtqMP+hn27lttHe3\nMzh5MIvzFuPzRo7YGZY2zO2Y4oKG1i6e+90RXt9ThXdEOk/fOw3flGy3Y8VczIqVMSYJOA58DqgC\nyoEvW2sPf+Q1PmCntbbNGPMPgGOt/Zsrva+KlYhIfOrs6WRnzc7ogdH17fUkmSRm58zGyXfwFfjw\nZnrdjik32I5TF3n8jYOcvNDK3dNzefKeqeQMHeR2rJiJZbFaCHzPWvv53vFqAGvtv13m9bOAn1hr\nF1/pfVWsRETiX9iGOVR/KLqVwweNHwAwMWtidL+saaOm4TE62Hcg6OzuYW3ZKX686QNSkjw8eudk\nvrqwkKQEWNwey2L1V8Bd1tqHesdfBeZba79xmdf/BKi11j57pfdVsRIRSTzB5iClwVL8QT976vbQ\nY3sYlT6KkvwSfF4f83PnMyg5ca5iyMc7e7GV77xZweYT9czIH8Zz909nWl583yp2pVgZY1YA3wBK\nrLWdH/PrK4GVAAUFBbedPXvFZVgiIhLHmjqb2Fy9mUAwwJbqLbSGWklPTmdh7kIcr0OJt4QRg0a4\nHVP6iLWWdw7U8PQ7h7nY0snfLSzkn+6cTOag+NyI9obfCjTG3AH8mEipOv9JH6wrViIiA0dXTxe7\na3dHt3Koa6vDYJiZPTO6X9a4YePcjil9oKk9xAvvHmPdzrNkZ6bxvXtu5a5pY+Ju245YFqtkIovX\nPwtUE1m8/hVr7aGPvGYW8DqRK1snriagipWIyMBkreVow9Ho4vcjDUcAKBxaGC1Znxn9GZI8ifnY\n/kC1P9jIY+sPcrjmEr4po3n63ml4Rwx2O9ZVi/V2C18AfkRku4WfW2u/b4x5GthtrX3LGLMRmA7U\n9P4nldbaZVd6TxUrEREBqGmpIVAVwF/pp7yunO5wN8PThlOcX4yvwMfC3IUMTomfP4Dl8rp7wryy\n7Qw/fO84YWv51mcn89CScaQk9f+HG7RBqIiIxJ3mrma2ntuKv9LP5urNNHc1k+pJZcHYBfi8Pkry\nSxg9eLTbMeU6nWts56m3D/HuoTom5wzhufunM6ewf6+3U7ESEZG4FgqH2Fu3N7r7e3VLNQAzRs3A\n8To4XoeJWRPjbq2O/L+Nh+v47luHqG5sZ/lcL99eejNZg/vnbv4qViIikjCstZxoPBEpWZV+Ki5W\nAJA/JB/H63B7we3Myp5FsifZ5aRyrVo7u1nzvyd4ectpstJTePzuW7h/Vl6/K8wqViIikrDOt50n\nEAwQCAbYWbOTrnAXQ1OHUpxfjON1WDx2MUNSh7gdU67BkZpLPPbGQfZVNrJw/EievX8aE0b3n/+H\nKlYiIjIgtIXa2HZuG/6gn7KqMho7G0nxpDBvzLzoLcMxGWPcjilXIRy2/LK8kh/8/igdoTB/70zg\nH50JDEpx/wlRFSsRERlwusPdvH/hffyVkSN2KpsrAbhlxC34Cnz4vD6mDJ/S724zyZ+60NzJ9397\nmDf3n6Nw5GCevW86RZNGuZpJxUpERAY0ay2nL53GXxnZL+v9C+9jseRm5EavZM3NmUtKUnzuBD4Q\nbDlRzxO/qeB0fSv3zhzLd+6eyujMNFeyqFiJiIh8RH17PZurNuMP+tl+bjsdPR0MSRlCUV4RPq+P\novwihqYOdTum/JmOUA8/DZzkp4GTpKV4+Ne7buYr8wrw3OCDnVWsRERELqO9u52dNTuju783dDSQ\nbJK5Lec2fAWR3d/zhuS5HVM+4uSFFp54s4JtJy8yqyCL7983naljb1wRVrESERG5CmEb5mD9wegt\nw5NNJwGYNHwSPm9kXdbUkVPxmP6/O3iis9by5v5qnn3nCI3tIb6+uBBnSjb7g40sGD+S224a3mef\nrWIlIiLyKZy9dDa6Kem+8/sI2zDZ6dmUeEvweX3My51HWpI763wkorGtix9sOMovdwUB8BhITfbw\n6kML+qxcqViJiIhcp8aORjZXR9ZlbaneQnt3O+nJ6RTlFeF4HZbkLWH4oL67SiJX9vgbB3l1Z+TJ\nzyQDj9w5hVW+iX3yWVdbrLRFrYiIyGVkDcringn3cM+Ee+js6WRXza7oxqTvnX0Pj/EwK3tW9JZh\nwdACtyMPKF+anc+v91YR6g6TkuxhwfiRbkfSFSsREZFrFbZhjlw8gj8Y2S/r+IfHARg/bDw+b2Tx\n+4zRM7Qu6wbYc/ZDdpy6qDVWKlYiIpIoqluqo+uy9tTuodt2M2LQiMh+WfkOC8YuID053e2Ych1U\nrERERFxwqesSW6q2EAgG2Fy9mZZQC4OSBrFg7AJ8Xh/F+cWMSnd3F3G5dipWIiIiLgv1hNhdtzu6\nX1ZNaw0Gw4zRM6LrssYNG6cjduKAipWIiEg/Yq3l+IfH2RTcRCAY4PDFwwDcNPQmnPzIETszs2eS\n7NFzZf2RipWIiEg/VttaS2mwFH/Qz87anXSHu8lKy6I4vxif18eisYsYnDLY7ZjSS8VKREQkTrR0\ntbD13FYCwQBlVWVc6rpEqieV+bnzowdGZw/OdjvmgKZiJSIiEoe6w93sO78vspVDpZ+qlioApo2c\nhuN18BX4mJQ1SeuybjAVKxERkThnreVk40kCVQH8lX4O1B8AIG9IXqRkeX3MzplNiifF5aSJT8VK\nREQkwVxou0BpVSmBYIAdNTvo7OkkMzWTJXlL8Hl9LM5bTGZqptsxE5KKlYiISAJrC7WxvWY7gWCA\n0mApH3Z+SLInmbk5c6NXs3KH5LodM2GoWImIiAwQPeEeDtQfwF8ZOWLnzKUzANw84uboETu3jLhF\n67Kug4qViIjIAHW66XT0sOh95/dhseQMzoleyZo7Zi6pSalux4wrKlYiIiJCQ0cDZVVl+Cv9bK/Z\nTnt3OxkpGSweuxhfgY8leUsYljbM7Zj9noqViIiI/ImO7g521e5iU+UmSqtKqW+vJ8kkMTtndvSW\noTfT63bMfknFSkRERC4rbMNU1FcQCAbwB/180PgBABOzJkbPMbx11K14jMflpP2DipWIiIhcteCl\nYGS/rKCfvXV76bE9jE4fTYm3BJ/Xx7wx8xiUPMjtmK5RsRIREZFPpamzibKqMgLBAFuqt9DW3UZ6\ncjqLxi7C8TqU5JcwfNBwt2PeUCpWIiIict26eroory3HH/QTCAaoa6vDYzzMHD0zui6rcFih2zH7\nnIqViIiIxJS1liMNR6Il62jDUQAKhxbiK4isy5oxagZJniSXk8aeipWIiIj0qZqWmmjJKq8tp9t2\nM2LQCIrzi3G8DgtzFzI4ZbDbMWNCxUpERERumOauZrZWb2VTcBNbqrbQHGomLSmNBbkL8Hl9lHhL\nGJU+yu2Yn5qKlYiIiLgiFA6xt24v/qAff6Wfc63nMBimj54eWZeV7zAha0JcHbGjYiUiIiKus9Zy\n/MPj0SN2Ki5WAODN9EaP2JmVPYtkT7LLSa9MxUpERET6nbrWOkqrSvEH/eys2UkoHGJo6lCK84vx\neX0szltMRkqG2zH/goqViIiI9GutoVa2ndtGIBigtKqUps4mUjwpzMudhy8/spVDTkaO2zEBFSsR\nERGJI93hbvaf3x89YqeyuRKAqSOn4ngdbvfezuThk11bl6ViJSIiInHJWsvpptNsCm4iEAxw4MIB\nLJaxGWNxvA6O12HOmDmkeFJuWCYVKxEREUkI9e31lFWV4Q/62XFuBx09HWSmZFKUV8S4YePosT0U\n5RUxM3tmn2VQsRIREZGE097dzo5zOwhUBXjv7Hs0dzVjMKQlpbH2zrV9Vq6utlh5+uTTRURERPpA\nenI6vgIfTy16iq9N/RoePFgsoXCI3XXuX7BRsRIREZG4NC93HqlJqSSZJFI8KczJ+cQLSn2uf+/G\nJSIiInIZM7NnsvbOteyu282cnDl9usbqaqlYiYiISNyamT2zXxSqP9KtQBEREZEYUbESERERiREV\nKxEREZEYUbESERERiREVKxEREZEYca1YHTt2jFdeeQWAUCiE4zisW7cOgLa2NhzH4bXXXgOgqakJ\nx3FYv349APX19TiOw9tvvw1AbW0tjuOwYcMGAILBII7jsHHjRgBOnTqF4ziUlpZGP9txHLZt2wZA\nRUUFjuNQXl4OwP79+3Ech/379wNQXl6O4zhUVFQAsG3bNhzH4dixYwCUlpbiOA6nTp0CYOPGjTiO\nQzAYBGDDhg04jkNtbS0Ab7/9No7jUF9fD8D69etxHIempiYAXnvtNRzHoa2tDYB169bhOA6hUAiA\nV155BcdxonO5du1a7rjjjuj4xRdfZOnSpdHxmjVrWLZsWXT8wgsv8MADD0THzz//PMuXL4+On3nm\nGVasWBEdP/nkkzz44IPR8erVq1m5cmV0/Oijj7Jq1aro+OGHH+bhhx+OjletWsWjjz4aHa9cuZLV\nq1dHxw8++CBPPvlkdLxixQqeeeaZ6Hj58uU8//zz0fEDDzzACy+8EB0vW7aMNWvWRMdLly7lxRdf\njI7vuOMO1q5dGx07jqPvnr57gL57+u7puwf67l3rd++T6IqViIiISIzorEARERGRT6CzAkVERERu\nMBUrERERkRhRsRIRERGJERUrERERkRhRsRIRERGJEdeeCjTGXADOXufbjALqYxBH/pTmtW9oXvuG\n5rVvaF77hua1b9yIeb3JWjv6k17kWrGKBWPM7qt59FGujea1b2he+4bmtW9oXvuG5rVv9Kd51a1A\nERERkRhRsRIRERGJkXgvVj9zO0CC0rz2Dc1r39C89g3Na9/QvPaNfjOvcb3GSkRERKQ/ifcrViIi\nIiL9RlwWK2PMXcaYY8aYD4wx33Y7T6IwxvzcGHPeGFPhdpZEYYzxGmP8xpjDxphDxphvuZ0pERhj\nBhljdhlj3u+d16fczpRIjDFJxph9xph33M6SKIwxZ4wxB40x+40xu93OkyiMMVnGmNeNMUeNMUeM\nMQtdzxRvtwKNMUnAceBzQBVQDnzZWnvY1WAJwBhTDLQA/2WtneZ2nkRgjMkFcq21e40xmcAe4D59\nX6+PMcYAGdbaFmNMCrAF+Ja1dofL0RKCMeYRYA4w1Fr7RbfzJAJjzBlgjrVWe1jFkDHmF8Bma+1L\nxphUYLC1ttHNTPF4xWoe8IG19pS1tgv4FXCvy5kSgrW2DGhwO0cisdbWWGv39v64GTgC5LmbKv7Z\niJbeYUrvP/H1t8R+yhiTD9wNvOR2FpErMcYMA4qBlwGstV1ulyqIz2KVBwQ/Mq5Cf1BJHDDGFAKz\ngJ3uJkkMvber9gPngfestZrX2PgR8C9A2O0gCcYCfzDG7DHGrHQ7TIIYB1wA/rP31vVLxpgMt0PF\nY7ESiTvGmCHAr4GHrbWX3M6TCKy1PdbamUA+MM8Yo9vX18kY80XgvLV2j9tZElCRtXY2sBRY1bv0\nQq5PMjAb+Km1dhbQCri+7joei1U14P3IOL/350T6pd41QL8GXrXWrnc7T6LpvfTvB+5yO0sCWAws\n610P9CvgdmPMOncjJQZrbXXvv88DbxBZ1iLXpwqo+sjV6teJFC1XxWOxKgcmGWPG9S5UWw685XIm\nkY/Vu8j6ZeCItfaHbudJFMaY0caYrN4fpxN5mOWou6nin7V2tbU231pbSOT31k3W2hUux4p7xpiM\n3odX6L1VdSegp6+vk7W2FggaY6b0/tRnAdcfDEp2O8C1stZ2G2O+AbwLJAE/t9YecjlWQjDG/BJw\ngFHGmCrgu9bal91NFfcWA18FDvauBwJ4zFr7OxczJYJc4Be9Twl7gP+x1mprAOmvcoA3In/PIhn4\nb2vtBncjJYxvAq/2Xmg5BTzocp74225BREREpL+Kx1uBIiIiIv2SipWIiIhIjKhYiYiIiMSIipWI\niIhIjKhYiYiIiMSIipWIiIhIjKhYiYiIiMSIipWIiIhIjPwfkjAIiO3fbe8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fb9e78bfe50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def get_xy(d):\n",
    "    xs = sorted(d.keys())\n",
    "    ys = [d[x] for x in xs]\n",
    "    return xs, ys\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10, 7))\n",
    "\n",
    "x, y = get_xy(all_cls)\n",
    "ax.plot(x, y, '.-', label='cls')\n",
    "\n",
    "x, y = get_xy(all_clb)\n",
    "ax.plot(x, y, '.-', label='clb')\n",
    "\n",
    "x, y = get_xy(all_clsb)\n",
    "ax.plot(x, y, '.-', label='clsb')\n",
    "\n",
    "ax.axhline(0.05, ls=\":\", color='k')\n",
    "ax.legend(loc=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy\n",
    "from sympy import symbols, diff, init_printing, simplify, lambdify\n",
    "from sympy.solvers import solve\n",
    "from sympy import stats\n",
    "from sympy.plotting import plot3d\n",
    "init_printing() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "metadata": {},
   "outputs": [],
   "source": [
    "ns, theta = symbols('ns theta', real=True)\n",
    "NBEXP, SIGMA = symbols('NBEXP, SIGMA', positive=True)\n",
    "nobs = symbols('nobs', positive=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmgAAAA3BAMAAABKhqtxAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAInarRM2ZVBDdiWbv\nuzJCz3LGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAKcElEQVR4Ae1afYxU1RX/zezM7JuvN6MmbSkk\nO64Rg7TuyIKpiuElbUMwNDsNXUoT032Crkm1MrUtCX8gT0iMUMwOioZYdR/4UUtNGWljk1p1rbVV\nU2UpRqGV7qRN+pU2LJ8uH7vT37nvvZnZnf1gQZjd7J7JnHvuueece+959+vdd4BJAoHWSdLQCdXM\n5ydUayZJY745AdvZUBSYgA3zmnTEIyZQOn8CtaWiKdG7/zYbiMy2A0evsrDzmoqi2pO+TO3bMFwL\noi35WBLbcLN/u98MGf8bTqZmvFC2VPXMEjUe4kf58Uifk+xLlPrcHCSsUBrPxFPRTPyjh85J8VIJ\nbYKv/VBBavN1A/6e48DasweAXxSbmw8Zwulvbu7MINpzKqW19dPFVxX337LcEg2BwIE1DlHC1Ta+\nWNzffOMJjGijpOoSdUkS87AUsSz+HDNCaX3RwFCZmubvwkZE6CkgKE2N9xOtkOx7QEiamsiwt91A\n3Rl60pISdj/+thAC/7USliLKqNoGNfAq/yPYKKs6lMYq8T7a9URSG0hYsez36MLagr5rX7kBWg5f\nAv4ojJcFLe6h62wSOn0UOEmiJQVoJhknoKeZh34aqD8qlMCPkcgqoj6lEqIqGwGRfpD/YW3c5+mV\n050c973Ydn+wUGeHrVfxNO4tF9aEimNJud5IHu8AiwrkbBXuCy0mAuIFP3f6qHR2d4FOS5FoK1xB\nzBKOSz8dpyDahUReUWWnVdmIyIjdIFLD2RjGaVv4wNK4Nu+b9xPoN9l4ZY8h2jWEIDqsUvWb2BFL\nOc2XE6YZ7kVEiuNdwJbtJNpk5CVJrL8+TcySXs4zNaGZCZs4nFXsstOqbARpEz6RGs7GME4LGsrk\nhEJry62Zq8hOi4PHFNKMnkFIiEQa+qMWiX7488JArN9SqSxyDSa0e4rFk2gx0GQodoXThtpIdMOp\nZ6gNpeg67dbrL78qxZNZ42tEpiqZMOiVFLSvq9YEDCY5IQMy2eLiGs7MUxwVhMPzGg9mhXOieW9S\nGIjIZCUctlH/BPDvO9buS2JRc/Neh11yWrWNlmULOdkFhthQPNdp5qGURqlv4EOuCxlVMlFQJJnG\nRku15gtsmT8lZDgtKEvEmbna+IrwmpLQ3hTOAIKS55Q866RNNz2/koV/QLQXeExtplJQclq1DS6L\nrEtgsA2H5zhNzz4G34Cs//+irZxTNEHwFfi2dvtmaYxvA9elsGpWu+CYQRTizOt+QbIf899EP3I9\nCkEXzu2LmCX0WIJDGQT6gNPQP2Eu0tr6rbdaW9OqpMoGl0UTBSkbbIOM51pbn2xtXca9RjsuT0E7\newe5eq8Ijwb1PXkseReBnhtW7QN+eKx9wRFEO7ei5wPga8fa9ycdZe07Kc/KujMIPGmPqeGJV6bP\n7OoLF08pjsYj2a+FclYz5TTOzPgRU5jH+G+yubZ1Q1M7X6hwuFsKOH8FuGf6c9A+QaiLOUJ5pFXZ\nYEWW2oCH2HAU3elZ34u6DHtcZC1jOw2v5uAvAL8BXmQPuoADQIeNZ2mTO3oo59iOWyscgtyDwHXM\njKHhiVem+/SuUnZ1ErI9Yg18SXd6LuYidtImT45pWGQ4xzTI2HsWYbUw6XJ442AqIMjd4DgSKZUv\nO63KhjqmRSyKDbHhKLpO49BN5K3PYgNrCbh9dgSGxTuOcwYAq7gjW6ijo22O/t6oQV5dGv5epoSV\n+MghOBkaDI73MTU88co0F0mVsg1m1GZGN+Gn09RGIN7pNIjkmKZzhKhjmp4mw0BUucsv05FP1eLL\nAPA7/FRlK0ZalQ11TPs7xfxDbDiartNiNrtv3wawv9GMUzQKtg8ppz0C3MmlxdatPFecoz8QjXAe\n8W6lqvXx8O1ChDyT9OganvSg9EDQKuXr+sKSmbHr+RuZREyiHP+v849gF7CUIx6d/C8tADcwLVpE\nwT4i+jjr28Zkpe8Nla1wWpUNHv0QeJ9iQ204mq7TOBCeixb2qb0pZLtGR0y0VEuSx2VtAP4VfIJG\nXEkuzErSUcAcCw/uLR7xZ8CF24ENgSNaalQNR+7l1xpNDL6L+s8Op0iwfnqNJG1FtchxSmgfn8pz\ndJNX13as+RaObG0B39sP8gWq7WwBXy0+yinQ1p8SLW3ltRaTB2cLFnDXtGobG8UWX/mrbTiKrtOW\nAFf8CpfPvYbcYMEpGhnXF8Imh1XgaON+Cq1vf0qJvs7mc0XZtSsL7ZGZR6zg8lXOIyb3fjxeb4ym\nIaqEX57BGxjlLuo2PqQybC2T50OVNoLxKrtOq1SbUZkZlo4g2mdzemTg3w4uERx1nJ5X5yR5TlDI\n5Ct/LBtVHGHYeF1ODSNriJAcKnLoQ8VdlLrYdpEIrDcFe/CyR5xfqiXPTw8PVOvtrGZ5nHWHCJas\nZ781uUiY8HGNXY4HNA6yzQE1GT8U2USSfU8k4/l6UTgArmex2eSPrOEK1qUCXZC7qHWiaIktBd49\nR3BQP+ODcp5sLVItM2atPBu9aNM1WdRzy6LfIhY0W87aKseFTU7ICcwqWbIReZeZkTVcyVghZBaG\nuYsadM9RsurrLpE1Js7h8fFc05HlK52FWTZ0+o3L88OcOwUu1cxxC7W4MAaTf/W6wpc7HipH03Al\nGxBeXBjmLmrQPYdnFZhZJmtLvTRm9aF3DISSiHS2r/qAS9rJ9gUD2HjW0to4mK4+KVNGn3c7Ubvh\nmbpnG99HOCBH1nAlvwz/zyrvorxpyRttz9akSjc3stuXGkrT0r3nuNT1X2B9ga34/AWaOA/10rTc\naJ2Hds1VZmW1e2vRiLWobztlq3uOyvPHZKB7sXDP3Fo8bE7L7/NVybvnqMVju4A65eKlBsBpuYK/\nSQpyv3zpQaZlxpe+9BV/OjXyOOH/dCyNx4pMy1vVveN4tCaM7NXQanDkuGjdv9zZSXovWgXKsH/l\nKxe3gotv/TPfLb2iYP5I1T38jlfCW6Vp0NLYbbhuGCVca5knUpx2Ga8pLcS6XEeE+Oo8Asg7soB2\nykmnNo7lUCchGgKbnGQYLF/5HJDAoSkP/FRQctpdI3ojkvOK3vSIKZ6G3WGk0TNN3D9P8pPB3mLR\nLEVlzEdw3mxO3Z1zC7xcnQbxQIvt+CGSR+Tn9828Mgn98cVPXQYnKoMfqzqSHSsYnxE0eMf8J0d2\nymPvFMEl7Z/Y4S8wNivLD4BuVEbMRthebzFwoQG6hbunvLuUA/i9x4G5knRzO5Xvti2WG5Wx2mK4\n1p3QTqPj94Z8XJ8GekC+7QQMIlns9XSY+Awwh19AVVTGewxdSC1H4BPGnPBT/XUsnoZ6E5ehFK4V\nyicsCZWSBV9FZTzEO4ndSUYI9vmyWlMS1rTH6IF/AFeWw7U6rIQBhk7ob0tkhURlcNS9y6CALflg\nFz+lTYN4wPfErlUMZODEU9cmc8RVDNGZYXBWqqiMVfJd9mbGNfKDXCO0t0RpqkOQBzMuZqvdcK2/\noM7mMa1RQoaWqKiMyM47JNQjD1/jHp7VvK12qvtN+u+Ga027YjwecMO1xqMyLeuGa007YlweGByu\nNS7VqSs8OFxr6vphXD0fHK41LtXJJ/x/r1n6sbRp0xYAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{\\sqrt{2} \\left(NBEXP \\sigma \\theta + NBEXP + ns\\right)^{nobs}}{2 \\sqrt{\\pi} \\sigma nobs!} e^{- NBEXP \\sigma \\theta - NBEXP - ns - \\frac{\\theta^{2}}{2 \\sigma^{2}}}$$"
      ],
      "text/plain": [
       "                                                                       2   \n",
       "                                                                      θ    \n",
       "                                     -NBEXP⋅SIGMA⋅θ - NBEXP - ns - ────────\n",
       "                                                                          2\n",
       "                               nobs                                2⋅SIGMA \n",
       "√2⋅(NBEXP⋅SIGMA⋅θ + NBEXP + ns)    ⋅ℯ                                      \n",
       "───────────────────────────────────────────────────────────────────────────\n",
       "                              2⋅√π⋅SIGMA⋅nobs!                             "
      ]
     },
     "execution_count": 135,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poisson = sympy.stats.density(sympy.stats.Poisson('X_ns', ns + NBEXP * (1 + SIGMA * theta)))(nobs)\n",
    "gaussian = sympy.stats.density(sympy.stats.Normal('X_G', theta, SIGMA))(0)\n",
    "likelihood = poisson * gaussian\n",
    "simplify(likelihood)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAScAAAAVBAMAAAD2sTVgAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAZpkQ3Ynvq81UMrtE\nInZCK3CHAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAESUlEQVRIDc1WTYhbVRT+kpfMJC+/bhQVbIpU\nu9AxjHZEqpOq1akgTLRFikg72BbFCo1SG7C2DfWnIwjJyo3gPHDhwsUMKhYpMm9RS0sXBhcWN02m\nmyKFTGo7M7WJE79z78t7mUxqXLTqgdx7z3fOPfe75957XgD41uL/I+ZQWsgYz/yHlIzbuhePviGI\nr96Nd+mXWjbC43u5gzta+fwJG7jUOp7ff41TZ69WjWKjBFxsHT94VO2xa7JSN+7uhQp29rGUZxp+\n4WUq6wUITkj7dzI7A4TUinsB/3W6khDu5m/wMpkqC5HIrhsGOXYDi7kLVzzTx3ikDGwToC+p6Gd/\n0ktcTXKI/s7fMpVhAa7BzCnLAhAWlPK57jrbPzqVjnEgCc/km0CwALwk9r6kgo+TxxlxjZOdj0sH\nJVujghStc9Ih3uSPxERWk4p2pEP7OO3tAPfpSCSHMB3/GSk/xoGtMjMyAYw8zbOsM3MxQaYOqUQh\nUudZkpjIalLBTcqwunkKppepRA5Rai6pHw+dv1hlBmp3ycSAXkiGIn6sSSIlI04zx9LsZzAnOjDQ\noEpJFIBKSkYuKTfmewgdvlACdsxZyu410aZKMTAmWCal7oUiNVgAUieqBo/mE3wt1hEiHXIGiU3I\nCjB9uPYto+OnYwfoLRJc1v10FuF2aXEypWPGnkemnHmf0UN2JKWd3XZwKX+gTs14XaBKFiZvwBoZ\nf2rBLG1B7DpidXwhSJdshb8ZtQWcL8N4jf1OCwXReWSLup9/58Whsh46x+fEHMgikJ1KYwEVmGnH\no90FCirFjspMC6lQCviKxcowmlKujMUH2u7EPthM2WIxjYg1giredzTPl3i9BaQFePAIVcpse73t\nk5PPTk6yBDgxH0rz0NfDWEDmVVu5djSJLKb1fAGd48P+KnDK4muuQ07xbGuGxm7Jcs2fFbjEdp5q\nA0gbVSp+a1rPuKrsunGOT8dkZatUj8oNjs5y2krJlHCfh/DGhnnRQzlCwsbPPCbTv2LUuSqeJ4tS\niZdonSBSpnDE1mVK5e4VBNQMk0XDFYeUivkwz2Nn+TL8V2Ilg6e/UjI2lj0kMgMfS4K6U/JeefJT\n6exHQN3zaY9k9YB67VKmTG5XlakvaY/b8Ck6ce7QFYeUjsl97EETI8nQBBKuizMYsGXxtrB4ciNe\nSajY2O6zfkC8Kh4rX1+AiP4+Mi5OfUOKpK6+m2/T1CJnhLhDVxxSOub94Jt7i686UkYNhnyfPPFX\nfymLpl8fnsOI1UHqJHDuTpyf+01NWFGnNu57kmCOv8HiUv7g98Cj0o8XYBYXLZxujfH0i42qmqka\nh5SOGdzBxzN8IYlY7Z4SsFm24IoxdK8eMwZlQ+1dtm7xVNjNahxSPcOFVWZ6mtrgrSH1RDt8j159\n13vgHZAipe9LB3orh6P9g6v/U/iwv+NN80j2jRR9U7mE1/b1/PcczHV8Cn8BsQAtssqYwf4AAAAA\nSUVORK5CYII=\n",
      "text/latex": [
       "$$\\left \\{ ns : - NBEXP + nobs, \\quad \\theta : 0\\right \\}$$"
      ],
      "text/plain": [
       "{ns: -NBEXP + nobs, θ: 0}"
      ]
     },
     "execution_count": 136,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "solution = solve([diff(likelihood, ns),\n",
    "                  diff(likelihood, theta)], ns, theta, dict=True, check=True)[0]\n",
    "ns_hat = solution[ns]\n",
    "theta_hat = solution[theta]\n",
    "solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAVoAAAFBCAYAAAA2QGcuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XuUXVWB5/Hvz8QEpEN4REM6AUMLokC3CGnU0bFpnlFp\ncJaoYVCDC0wzKtrLWaPQjtAdu9eI/fCxWmWVAQWVhxNEoyIYQFrHbiIB0iAEJEYlFYEAiYjaJEZ+\n88fdhYfLrapbSZ06dat+n7XuqnP22WfffRP4Zdc++5wr20RERH2e1XQHIiImugRtRETNErQRETVL\n0EZE1CxBGxFRswRtRETNGg1aSQsl3StpnaRzOhx/taTbJG2XdErbscWS7iuvxZXyIyTdWdr8pCSN\nxWeJiIlB0h6Slku6R9JaSa9oOz5T0tcl/YekuyS9fbg2GwtaSVOATwGvAQ4GTpV0cFu1+4HTgcva\nzt0LOB94GXAkcL6kPcvhzwDvAA4sr4U1fYSImJg+AVxr+0XAS4C1bcffBdxt+yXAUcA/SZo2VINN\njmiPBNbZXm97G3AFcHK1gu2f2r4DeLLt3BOAlbY3294CrAQWSpoD7G77ZrfuxLgUeH3tnyQiJgRJ\nM4FXAxcB2N5m+xdt1QzMKL8t/wGwGdg+VLtNBu1cYENlv7+U7cy5c8v2jrQZEbE/8DDwOUm3S1om\nabe2Ov8CvBj4OXAn8F7b7YPBp5laS1d7gKQlwBKAZ+3y7COes99eO9Xetm2j/0epbfVPL0/ZVvtb\njNiUJ3Jb+GTx+K82PmL7uSM9778etYu3bB4y2zq6687f3gU8USnqs91X2Z8KHA6cbXuVpE8A5wAf\nqtQ5AVgDHA28AFgp6Xu2fznY+zYZtBuBfSv780pZt+ce1XbuTaV8Xjdtlj/cPoAZB+3jl376rV2+\ndWf398/aqfM7mf6zIad9RsWM+8dfqM388damuxBj5Mab/vpnO3Lels1PctU3R/7/3Iv2e+AJ2wuG\nqNIP9NteVfaX0wraqrcDHynTk+sk/QR4EfCDwRptcurgFuBASfuXieRFwIouz70OOF7SnuUi2PHA\ndbYfAH4p6eVl/uRtwNfq6HxETDy2HwQ2SDqoFB0D3N1W7f5SjqTZwEHA+qHabWxEa3u7pHfTCs0p\nwMW275K0FFhte4WkPwWuBvYE/kLS39o+xPZmSR+mFdYAS21vLtvvBD4P7Ap8q7x6zliMZiOio7OB\nL5UB4Hrg7ZLOArB9IfBh4POS7gQEfMD2I0M12Ogcre1rgGvays6rbN/C06cCqvUuBi7uUL4aOHR0\nexoRk4XtNUD79MKFleM/p/VbdNdyZ1iMK5mfjYkoQRsRUbMEbUREzRK0ERE1S9BGRNQsQRsRUbME\nbUREzRK0ERE1S9BGRNQsQRsRUbMEbYwbuSssJqoEbUREzRK0ERE1S9BGRNQsQRsRUbMEbUREzRK0\nERE1S9BGRNQsQRsRUbMEbUREzRK049TW529rugsRMUoStBERNUvQRkTULEEbEVGzBG2MC3lyV0xk\njQatpIWS7pW0TtI5HY5Pl3RlOb5K0vxSfpqkNZXXk5IOK8duKm0OHHve2H6qiOhlkvaQtFzSPZLW\nSnrFIPX+VNJ2SacM1+bU0e9mdyRNAT4FHAf0A7dIWmH77kq1M4Attg+QtAi4AHiz7S8BXyrt/DHw\nVdtrKuedZnv1mHyQiJhoPgFca/sUSdOA57RXKPl1AfDtbhpsckR7JLDO9nrb24ArgJPb6pwMXFK2\nlwPHSFJbnVPLuRERO0XSTODVwEUAtrfZ/kWHqmcDVwGbumm3yaCdC2yo7PeXso51bG8HHgP2bqvz\nZuDytrLPlWmDD3UI5oiIwewPPEwrQ26XtEzSbtUKkuYC/w34TLeNNjZ1MBokvQz4je0fVopPs71R\n0gxa/+K8Fbi0w7lLgCUA0583Y6f7st+8R7i/f9ZOtxMRw9vyu+dw5WMLduDMr8+SVJ1W7LPdV9mf\nChwOnG17laRPAOcAH6rU+TjwAdtPdjuOazJoNwL7VvbnlbJOdfolTQVmAo9Wji+ibTRre2P5+bik\ny2hNUTwjaMsfbh/AjIP28U59kojoFY/YHiqh+4F+26vK/nJaQVu1ALiihOws4LWSttv+6mCNNjl1\ncAtwoKT9y4TzImBFW50VwOKyfQpwo20DSHoW8CYq87OSpkqaVbafDZwI/JCIiC7YfhDYIOmgUnQM\ncHdbnf1tz7c9n1YQv3OokIUGR7S2t0t6N3AdMAW42PZdkpYCq22voDUh/QVJ64DNtMJ4wKuBDbbX\nV8qmA9eVkJ0CXA98dgw+TkRMHGcDXyoDwPXA2yWdBWD7wh1psNE5WtvXANe0lZ1X2X4CeOMg594E\nvLyt7NfAEaPe0YiYNMpS0fbphY4Ba/v0btrMnWHA7lOfaLoLETGBJWgjImqWoI2IqFmCNiKiZgna\niIiaJWgjImqWoI2IqFmCNiKiZgnaiIiaJWgjImqWoI2IqFmCNiKiZgnaiIiaJWijcfmq8ZjoErQR\nETVL0EZE1CxBGxFRswRtRETNErQRETVL0EZE1CxBGxFRswRtRETNErTFsbPvbboLETFBJWgjImqW\noB1F+817pOkuRMQ41GjQSloo6V5J6ySd0+H4dElXluOrJM0v5fMl/aekNeV1YeWcIyTdWc75pCSN\n3SeKiF4naQ9JyyXdI2mtpFe0HVfJlnWS7pB0+HBtNha0kqYAnwJeAxwMnCrp4LZqZwBbbB8AfAy4\noHLsx7YPK6+zKuWfAd4BHFheC+v6DBExIX0CuNb2i4CXAGvbjr+G3+fLElqZM6QmR7RHAutsr7e9\nDbgCOLmtzsnAJWV7OXDMUCNUSXOA3W3fbNvApcDrR7/rETERSZoJvBq4CMD2Ntu/aKt2MnCpW24G\n9ijZM6gmg3YusKGy31/KOtaxvR14DNi7HNtf0u2S/lXSf63U7x+mzYiIwewPPAx8ruTLMkm7tdXp\nJrueZuro9nHMPADsZ/tRSUcAX5V0yEgakLSE1rCfmXN2raGLO2/r87cx/WfTmu5GxLjzy+27cP1D\nB+3AmV+fJWl1paDPdl9lfypwOHC27VWSPgGcA3xoJ7rbaNBuBPat7M8rZZ3q9EuaCswEHi3TAlsB\nbN8q6cfAC0v9ecO0STmvD+gDmHvIHt7pTxMRveAR2wuGON4P9NteVfaX0wraqm6y62manDq4BThQ\n0v6SpgGLgBVtdVYAi8v2KcCNti3pueViGpL+iNak9HrbDwC/lPTyMpf7NuBrY/FhIqL32X4Q2CBp\nYLh8DHB3W7UVwNvK6oOXA4+V7BlUYyNa29slvRu4DpgCXGz7LklLgdW2V9CakP6CpHXAZlphDK3J\n6qWSfgs8CZxle3M59k7g88CuwLfKKyKiW2cDXyoDwPXA2yWdBWD7QuAa4LXAOuA3wNuHa7DROVrb\n19DqdLXsvMr2E8AbO5x3FXDVIG2uBg4d3Z5GxGRhew3QPr1wYeW4gXeNpM3cGRYRUbMEbUREzRK0\nERE1S9BGRNQsQRsRUbMEbUREzRK0ERE1S9BW5OtsIqIOCdqIiJolaCMiapagjYioWYJ2lOULGiOi\nXYI2IqJmCdpJ7vH98iXBEXVL0EZE1CxBG4177AXTm+5CRK0StBERNUvQjnNbn7+t6S5ED5q2dgPT\n1m4YvmKMiV79uvGI6KA9XKv72168b3v1GCMZ0UZMEMONYDPCbU6CNmIC6DZEE7bNSNBG9LiRhmfC\nduwlaCMmoYTt2ErQRvSwnQnMhO3YSdC2ycO/I2K0NRq0khZKulfSOknndDg+XdKV5fgqSfNL+XGS\nbpV0Z/l5dOWcm0qba8rreWP3iSLGzmiMSDOqfSZJPy3ZskbS6kHqHFWO3yXpX4drs7F1tJKmAJ8C\njgP6gVskrbB9d6XaGcAW2wdIWgRcALwZeAT4C9s/l3QocB0wt3LeabY7/gFFRHThz213fOappD2A\nTwMLbd/fzWCuyRHtkcA62+ttbwOuAE5uq3MycEnZXg4cI0m2b7f981J+F7CrpNwwH5PGaI5EM6od\nsf8OfMX2/QC2Nw13QpNBOxeo/g338/RR6dPq2N4OPAbs3VbnDcBttrdWyj5XhvUfkjTmzwHMw78j\nepqBb5dpySUdjr8Q2LNMU94q6W3DNdjTt+BKOoTWdMLxleLTbG+UNAO4CngrcGmHc5cASwBmztl1\nDHobMX5NW7uhp27R3bZtKvf3z9qRU2e1zbv22e5rq/OqkiHPA1ZKusf2dyvHpwJHAMcAuwL/Lulm\n2z8a7E2bDNqNQPVvdl4p61SnX9JUYCbwKICkecDVwNts/3jgBNsby8/HJV1Ga4riGUFb/nD7AOYe\nsodH6TNFxPj2iO0FQ1WoZMgmSVfTypBq0PYDj9r+NfBrSd8FXgIMGrRNTh3cAhwoaX9J04BFwIq2\nOiuAxWX7FOBG2y6T0d8EzrH9/YHKkqZKmlW2nw2cCPyw5s8RMabqmlPNXC1I2q38Noyk3Wj9ttye\nIV8DXlXy5jnAy4C1Q7Xb2IjW9nZJ76a1YmAKcLHtuyQtBVbbXgFcBHxB0jpgM60wBng3cABwnqTz\nStnxwK+B60rITgGuBz47Zh8qInrdbODqcmlnKnCZ7WslnQVg+0LbayVdC9wBPAkssz3kgK7ROVrb\n1wDXtJWdV9l+Anhjh/P+Dvi7QZo9YjT7GDGZ9Npc7WizvZ7WNEB7+YVt+/8A/EO37ebOsB5Q98O/\n8wWNEfVK0Eb0kMyj9qYEbUQ8TcJ89CVoIyJqlqCNiKhZgraDPCoxJrtMH4yuBG2MC4+9IM8Eiokr\nQRsRUbMEbUSPGOtf5zN9MHoStDXJoxIjYkCCNiKiZgnaiIiaJWgjYlCZpx0dCdqIiJolaCMiapag\n7RF5VGJE70rQRkTULEEbEUPKBbGd19VX2Uh6HXAIsMtAme2ldXVqPDh29r1c/9BBTXcjIiaAYUe0\nki4E3gycDYjWd3g9v+Z+RURMGN1MHfwX228Dttj+W+AVwAvr7VZExMTRTdD+Z/n5G0l/CPwWmFNf\nl2KyyqMSY6LqJmi/IWkPWl+texvwU+DyOjs1UeTBMhEB3V0M+6jtrcBVkr5B64LYE/V2KyJi4uhm\nRPvvAxu2t9p+rFoWERNflnjtnEGDVtI+ko4AdpX0UkmHl9dRwHNG480lLZR0r6R1ks7pcHy6pCvL\n8VWS5leOnVvK75V0QrdtRkQMRdJPJd0paY2k1R2OnybpjlLn3yS9ZLg2h5o6OAE4HZgH/HOl/JfA\nX4+w788gaQrwKeA4oB+4RdIK23dXqp1Ba7XDAZIWARcAb5Z0MLCI1trePwSulzSwEmK4NiMihvPn\ntge7yPIT4M9sb5H0GqAPeNlQjQ0atLYvAS6R9AbbV+1wdwd3JLDO9noASVcAJwPVUDwZ+JuyvRz4\nF0kq5VeUueOfSFpX2qOLNiMidpjtf6vs3kxrMDqkbuZovy/pIknfApB0sKQzdrCPVXOB6sRPfynr\nWMf2duAxYO8hzu2mzZ6VB8tEjAkD35Z0q6Qlw9Q9A/jWcA12s+rgc+X1wbL/I+BK4KIuzh23yh/g\nEoCZc3ZtuDcRw9v24n1zUarQNjH9Z9N25NRZbfOufbb72uq8yvZGSc8DVkq6x/Z3n9EH6c9pBe2r\nhnvTbka0s2x/GXgSnhpZ/q6L84azEdi3sj+vlHWsI2kqMBN4dIhzu2kTANt9thfYXrDbnjv0FxYR\nveeRgf/vy6s9ZLG9sfzcBFzN76clnyLpT4BlwMm2Hx3uTbsJ2l9L2pvWcBpJL6f1K/zOugU4UNL+\nkqbRuri1oq3OCmBx2T4FuNG2S/misiphf+BA4Addttm1Y2ffu6OnRkQPkrSbpBkD28DxwA/b6uwH\nfAV4q+0fddNuN1MH76MVVi+Q9H3gubRCb6fY3i7p3cB1wBTgYtt3SVoKrLa9gtb0xBfKxa7NtIKT\nUu/LtC5ybQfeZft3AJ3a3Nm+RsSkMRu4unXNnanAZbavlXQWgO0LgfNoXSv6dKm33faCoRodNmht\n3ybpz4CDaD29617bv92ZT1Jp+xrgmray8yrbT9B6Wlinc/8e+Ptu2oyI6EZZsfSMdbElYAe2zwTO\nHEm73T74+8jy5ocDp0p620jeZDLL8w5GJg+WiYlo2BGtpC8ALwDW8PuLYAYurbFfERETRjdztAuA\ng8tFqIiIGKFupg5+COxTd0ciIiaqQUe0kr5Oa4pgBnC3pB8AWweO2z6p/u5FRFVuWuhNQ00d/COt\nVQYXAK+vlA+URQO2Pn/bjt4R05XH9xMz7s8sUcRoGuqhMv8KIOnZA9sDJOWe1YiILg01dfA/gHcC\nfyTpjsqhGcD36+5YRMREMdTUwWW0nkrzf4DqA7Qft7251l6NI8fOvpfrHzqo6W5ERA8baurgMVrP\nNDh17LoTEcPJBbHe0+2dYRERsYMStDHu5DbcmGgStGMgzzuImNwStBERNUvQRvSgbS/ed/hKMW4k\naCMiapag7UH5NtyI3pKgjYioWYI2okdlnrZ3JGi7kG/Djckuob5zErQxLuWmhZhIErQRPSwjzd6Q\noB0juTssYvJK0EZE1CxBG9HjMn0w/jUStJL2krRS0n3l556D1Ftc6twnaXEpe46kb0q6R9Jdkj5S\nqX+6pIclrSmvM8fqM0VMVJMtyCVNkXS7pG90OLafpO+U43dIem03bTY1oj0HuMH2gcANPP0bHIBW\nGAPnAy8DjgTOrwTyP9p+EfBS4JWSXlM59Urbh5XXslo/RYNyd1hEbd4LrB3k2P8Gvmz7pcAi4NPd\nNNhU0J4MXFK2L+Hp37I74ARgpe3NtrcAK4GFtn9j+zsAtrcBtwHzxqDPEePWZBt11kXSPOB1wGCD\nNAO7l+2ZwM+7abepoJ1t+4Gy/SAwu0OduUD1+zr6S9lTJO0B/AWtUfGAN5Qh/XJJ+a+vh2UtbTTg\n48D7gScHOf43wFsk9QPXAGd30+hQX864UyRdD+zT4dAHqzu2Lck70P5U4HLgk7bXl+KvA5fb3irp\nL2mNlo8e5PwlwBKAmXPy7enR++r4LrHxOlKesg1m3D/i2ACYJWl1Zb/Pdh+ApBOBTbZvlXTUIOef\nCnze9j9JegXwBUmH2h4smIEag9b2sYMdk/SQpDm2H5A0B9jUodpG4KjK/jzgpsp+H3Cf7Y9X3vPR\nyvFlwEeH6F9faYO5h+wx7N9Yvg03YkJ4xPaCQY69EjipXODaBdhd0hdtv6VS5wxgIYDtf5e0CzCL\nzhn2lKamDlYAi8v2YuBrHepcBxwvac9yEez4Uoakv6M1P/JX1RNKaA84icEntCMmpNEcgY7X0Wxd\nbJ9re57t+bQudN3YFrIA9wPHAEh6Ma1Afni4tpsK2o8Ax0m6Dzi27CNpgaRlALY3Ax8GbimvpbY3\nl8nqDwIHA7e1LeN6T1ny9R/Ae4DTx/JDDSd3h8VYmGwBWTdJSyWdVHb/J/COkjGXA6fbHvY34tqm\nDoZSfsU/pkP5auDMyv7FwMVtdfqBjmuPbJ8LnDuqnY2YhCZ7WNu+iTJVafu8SvndtKYYRiR3hkVM\nQJM9KMebBG0MKjctTE4J6dGXoO1hdd8dNh5kLe2OS2COHwnaiAlspGGbcK5HgjZigus2PBOy9UnQ\nRkwC216875BBmpCtVyPLu3pV7g6LXlcN1GlrNyRgx0hGtBGTVEJ27CRox1juDouYfBK0MaTxsJY2\nS7yi1yVoIyJqlqDtcZPhpoWIXpegjYioWYI2IqJmCdqIiJolaEfo2Nn3Nt2FiOgxCdroCVniFb0s\nQduAXrtpYTyspY3oZQnaiIiaJWgngKyljRjfErQRETVL0EZE1CxBGxFRswRt9Iws8YpelaDdAblp\nISJGopGglbSXpJWS7is/9xyk3uJS5z5JiyvlN0m6V9Ka8npeKZ8u6UpJ6yStkjR/bD7RyGUtbcT4\nJGmKpNslfWOIOm+QZEkLummzqRHtOcANtg8Ebij7TyNpL+B84GXAkcD5bYF8mu3DymtTKTsD2GL7\nAOBjwAV1foiImJDeC6wd7KCkGaXOqm4bbCpoTwYuKduXAK/vUOcEYKXtzba3ACuBhSNodzlwjKRJ\nMRTLWtqInSdpHvA6YNkQ1T5MaxD3RLftNhW0s20/ULYfBGZ3qDMX2FDZ7y9lAz5Xpg0+VAnTp86x\nvR14DNh7VHseERPZx4H3A092OijpcGBf298cSaO1fd24pOuBfToc+mB1x7YleYTNn2Z7YxnCXwW8\nFbh0hP1bAiwBmDln1xG+fTTlsRdMZ+aPtzbdjWjYlCe8o/8dzJK0urLfZ7sPQNKJwCbbt0o6qv1E\nSc8C/hk4faRvWlvQ2j52sGOSHpI0x/YDkuYAmzpU2wgcVdmfB9xU2t5Yfj4u6TJac7iXlnP2Bfol\nTQVmAo8O0r8+oA9g7iF7jDToI6I3PWJ7sAtYrwROkvRaYBdgd0lftP2WcnwGcChwU/kleh9ghaST\nbK/u2GLR1NTBCmBgFcFi4Gsd6lwHHC9pz3IR7HjgOklTJc0CkPRs4ETghx3aPQW40XZCNCKGZftc\n2/NszwcW0cqPt1SOP2Z7lu35pc7NwLAhC80F7UeA4yTdBxxb9pG0QNIyANubaU0631JeS0vZdFqB\newewhtYo9rOl3YuAvSWtA95Hh9UMseOyxCsmI0lLJZ20M23UNnUwFNuPAsd0KF8NnFnZvxi4uK3O\nr4EjBmn3CeCNo9rZQRw7+16uf+igsXiriBhjtm/i91OV5w1S56hu28udYQ3qtZsWImLHJGgnkMmy\nljbPPIhek6CNiKhZgjYiomYJ2oiImiVoY0SyxCti5BK00ZNyQSx6SYJ2J+QB4BHRjQRtw7KWNmLi\nS9BOMJNlLW1EL0nQRkTULEEbI5aVBxEjk6CNnpWVB9ErErQRETVL0EZE1CxBC+w55Tc7fG7W0kbE\ncBK048Bor6XNEq+I8SVBGz0tF8SiFyRoY4dkiVdE9xK0ERE1S9BGRNQsQRsRUbMEbfS8XBCL8S5B\nW7x55uodPnc01tJmiVfExJWgjR2WlQcR3WkkaCXtJWmlpPvKzz0Hqbe41LlP0uJSNkPSmsrrEUkf\nL8dOl/Rw5diZY/m5IqL3SZoi6XZJ3+hwbLqkKyWtk7RK0vxu2mxqRHsOcIPtA4Ebyv7TSNoLOB94\nGXAkcL6kPW0/bvuwgRfwM+ArlVOvrBxfVv9HifEg87Qxit4LrB3k2BnAFtsHAB8DLuimwaaC9mTg\nkrJ9CfD6DnVOAFba3mx7C7ASWFitIOmFwPOA79XY14iYJCTNA14HDDZIq2bXcuAYScPOoTUVtLNt\nP1C2HwRmd6gzF9hQ2e8vZVWLaI1gXSl7g6Q7JC2XtO+o9TgiJoOPA+8Hnhzk+FO5ZHs78Biw93CN\nTh2t3rWTdD2wT4dDH6zu2LYkd6jXjUXAWyv7Xwcut71V0l/S+pfn6EH6twRYAvCHc6fs4NtHRBP0\nxDamrd0wfMVnmiWpusSoz3YfgKQTgU22b5V01Ch08ym1Ba3tYwc7JukhSXNsPyBpDrCpQ7WNwFGV\n/XnATZU2XgJMtX1r5T0frdRfBnx0iP71AX0Ah/7JtB0N+lG137xHuL9/1qi1t/X525j+s2mj1l4n\nj+8nZtw/Lv74IrrxiO0Fgxx7JXCSpNcCuwC7S/qi7bdU6mwE9gX6JU0FZgKPPrOpp2tq6mAFsLhs\nLwa+1qHOdcDxkvYsqxKOL2UDTgUur55QQnvASQw+od1R02tpY+fkgljsDNvn2p5nez6t35ZvbAtZ\neHp2nVLqDDvSqG1EO4yPAF+WdAatVQNvApC0ADjL9pm2N0v6MHBLOWep7c2VNt4EvLat3fdIOgnY\nDmwGTq/xM0TEJCBpKbDa9grgIuALktbRyphF3bTRSNCWX/GP6VC+Gjizsn8xcPEgbfxRh7JzgXNH\nr6cRMRnZvokyVWn7vEr5E8AbR9pe7gyLiKhZgjYmlMzTxniUoJ3gxuLhMnnmQcTQErTjzGg/xSsi\nmpegjYioWYK2TdbSRsRoS9DGhJMLYjHeJGhjVOSCWMTgErSTQL7WJqJZCdqIiJolaMehLPHaeZmn\njfEkQRsRUbMEbQdZ4hURoylBG6MmKw8iOkvQxoSVedoYLxK0k0SWeEU0J0EbEVGzBO04lSVeERNH\ngjZG1Xi7IJZ52hgPErQ1yBKviKhK0A5iZ9bSRkRUJWgnkaw8iGhGgjYmvMzTRtMStONYr648GG8X\nxCKalqCNiKhZI0EraS9JKyXdV37uOUi9ayX9QtI32sr3l7RK0jpJV0qaVsqnl/115fj8+j9NREwU\nknaR9ANJ/yHpLkl/O0i9N0m6u9S5bLh2mxrRngPcYPtA4Iay38k/AG/tUH4B8DHbBwBbgDNK+RnA\nllL+sVKvEVniNb5knja6tBU42vZLgMOAhZJeXq0g6UDgXOCVtg8B/mq4RpsK2pOBS8r2JcDrO1Wy\nfQPweLVMkoCjgeUdzq+2uxw4ptTfIRNxiVdWHkQMzi2/KrvPLi+3VXsH8CnbW8o5m4Zrt6mgnW37\ngbL9IDB7BOfuDfzC9vay3w/MLdtzgQ0A5fhjpX6MsVwQi14laYqkNcAmYKXtVW1VXgi8UNL3Jd0s\naeFwbU6to6MAkq4H9ulw6IPVHduW1P4vRu0kLQGWlN2tL9rvgR92rvn1nXiXnTmXWcB4W3YwHvsE\n47Nf6VN3DtqRk365/eHrrn3oM7N24NRdJFV/Ve2z3VetYPt3wGGS9gCulnSo7Wo+TAUOBI4C5gHf\nlfTHtn8x2JvWFrS2jx3smKSHJM2x/YCkObT+5ejWo8AekqaWUes8YGM5thHYF+iXNBWYWep36l8f\n0Ff6s9r2ghH0oXbpU/fGY7/Sp+60hV7XbA87itxZtn8h6TvAQqAatP3AKtu/BX4i6Ue0gveWwdpq\naupgBbC4bC8GvtbtibYNfAc4pcP51XZPAW4s9SMihiXpuWUki6RdgeOAe9qqfZXWaBZJs2hNJawf\nqt2mgvYjwHGS7gOOLftIWiBp2UAlSd8D/i+ti1r9kk4ohz4AvE/SOlpzsBeV8ouAvUv5+xh8NUNE\nRCdzgO/FYJxzAAADv0lEQVRIuoPWCHWl7W9IWirppFLnOuBRSXfTGvT9L9sdf3MeoAz4WvO17fM0\nTUufujce+5U+dWc89qkOCdqIiJrlFtyIiJpN+qCVtFDSveW23cbndCVdLGmTpEGWm409SftK+k7l\nlsP3joM+dXWrZBPKOszb228db5Kkn0q6U9KaHb3SP9ok7SFpuaR7JK2V9Iqm+1SXST11IGkK8CNa\nVxb7aU1+n2r77gb79GrgV8Cltg9tqh9VZQneHNu3SZoB3Aq8vuE/JwG72f6VpGcD/w94r+2bm+rT\nAEnvAxYAu9s+sen+QCtogQW2x806WkmXAN+zvaw8r+Q5Q61F7WWTfUR7JLDO9nrb24AraN3G2xjb\n3wU2N9mHdrYfsH1b2X4cWMvv78Zrqk/d3Co55iTNA14HLBuu7mQmaSbwasqKIdvbJmrIQoL2qVt2\ni+rtvNFBeSLaS4H22xLHXBe3Sjbh48D7gSeb7kgbA9+WdGu5K7Jp+wMPA58r0yzLJO3WdKfqMtmD\nNkZA0h8AVwF/ZfuXTffH9u9sH0br7sAjJTU61SLpRGCT7Vub7McgXmX7cOA1wLvKFFWTpgKHA5+x\n/VLg10zgde+TPWgHbtkdUL2dNyrKPOhVwJdsf6Xp/lSVXzkHbpVs0iuBk8p86BXA0ZK+2GyXWmxv\nLD83AVfTmjZrUj/QX/ktZDmt4J2QJnvQ3gIcqNaDxKcBi2jdxhsV5cLTRcBa2//cdH+g61slx5Tt\nc23Psz2f1n9LN9p+S5N9ApC0W7mISfn1/Hiefu/+mLP9ILBB0sBDZY4BGru4WrfaHirTC2xvl/Ru\nWrfUTQEutn1Xk32SdDmt+6hnSeoHzrd90dBn1e6VtB7AfmeZEwX4a9vXNNinOcAlZeXIs4Av2x43\ny6nGmdm0nkIFrf/nL7N9bbNdAuBs4EtlkLMeeHvD/anNpF7eFRExFib71EFERO0StBERNUvQRkTU\nLEEbEVGzBG1ERM0StBERNUvQRkTULEEbPUXS/PLs0s+W59B+W9Kukt5Tnpd7h6Qrmu5nRFVuWIie\nUp4eto7Ws1XXSPoyrdumPwrsb3urpD0m8iP3ovdkRBu96Ce2B24FvhWYD9xB63bOtwDbm+pYRCcJ\n2uhFWyvbv6N1//7rgE/RegLULZIm9XM8YnxJ0MZE8CxgX9vfAT4AzAT+oNkuRfxe/tWPiWAK8MXy\n9SgCPpk52hhPcjEsIqJmmTqIiKhZgjYiomYJ2oiImiVoIyJqlqCNiKhZgjYiomYJ2oiImiVoIyJq\n9v8BWw7o8IfAwRMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1a60241690>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "example = {nobs:15, NBEXP:10, SIGMA:0.3}\n",
    "likelihood_example = likelihood.subs(example)\n",
    "example_solution = [float(sympy.N(solution[x].subs(example))) for x in (ns, theta)]\n",
    "likelihood_example = lambdify((ns, theta), likelihood_example, modules=['numpy'])\n",
    "X, Y = np.meshgrid(np.linspace(0, 6, 50), np.linspace(-.1, .1, 50))\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "pc = ax.contourf(X, Y, -2 * np.log(likelihood_example(X, Y)))\n",
    "fig.colorbar(pc)\n",
    "ax.set_xlabel(ns)\n",
    "ax.set_ylabel(repr(theta))\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAA3BAMAAAAmtoSFAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAInarRM2ZVBDdiWbv\nuzJCz3LGAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAJL0lEQVRoBdVafWxT1xX/+fv5ObZfizQNVVVM\nuqZqYY1pGOqYUJ60VYiqU1J1QVSa8IOwVP0SliahIa3jCVQ0aKe46rqhdhIeaynN/sDrWrqu03DV\ndrRsE6aTNlCb8SSkfalSkvGVAMH7nfee7Rjq2EASz0d655577jn33nPvuff63mOgCXBfXxMabWaT\nqtHM1pvQdiDdhEavbLK9KHAld27ycW1u6r2mWpddk/R1CR/8fYcB7H0MvR33Qu1o7jR75mEfeeMi\n3kVY9+XuNDuxMVa4rlGbLaVgZbxvudY6PbmyxjMVssyLaCR/BXi2ZTCBAajJj/EmRr5nliWaQeyE\nZ3DEkpY9BcB78gyw5dJx4NfF7u4RXThT3d17+hA5OZlQUlNp4LbiX1auMUXjVkE2+I9vLpHT0qWk\nfbTal/DncRjR7Ps0+1TRmiYy/+Qj2A6VVgIBmZDoFNFayf4ZCF5mGu8DIgX2+iJHwZSSsxT7UIgH\nBdnwmRk3XdJNvs/806QVqrdZQcMax5A1gcMhzZdwReYriQ0fqzSlZHAv8CdhHBS06iTNzpKI0T7/\nORK9CfbaIOMsYknmEbsAhP5LIkKuCy8hnrbJUMJOlOU9JqBKZi/QjvAqax1WYHekELaCFrnzCVGs\nrjSn5vAR0GOR85xw9/ca8EvfveO0SMw6YNHgBImUtYCYJfQHL41GmGPjQCSPeM4mXYOBLpNqefKe\nBb4O7y+x4K40Tg0j8s6XbMF5RAEMmeXmdtIQ0zbYkxGmER6DKsXRPPu6m0RKZlxM27o0ScySMfq3\nLIKFdlZQ2MBo2s5VGYxvkRfQ7YKmoi2V1hfb5B6Tk2YIaUQuIihEPInY8yaJKXhzwkDblGmnsqjb\nDZKd/JQnisVz6NXRpTNHX0/YiTPD2MeMaricJiXvJKA4m41fZxcy0g2/OGhUzKI3T0IV3ug9HSfS\nwjnbfVQTBlRxcMJoFqEXhfiE3783bjmmoae7+6hwrjT4TXIifXZBs5CqJbHdtFtfwp54E0KGk4LS\nRPTmTfo3hNelQXlPOJcRkDzd+JKTdq3Yt0ET8qec4MOIjAEv2Fu38KpnWFwplBF+02AB1inrd0nz\nnm1ch2G7I4OC23SiIL21sF+yn/Lr4hgExsiMCWd9D7OEk6adAO9Tvg/+CeACYufJVPv7H/6gvz8p\n5bJpYZRfbIxoRnjqIvw/y+KbpwdXagillg+kNLw9qUWndPhPLn9yqau8g/u9Cw1qiPQvhifCxUlb\nT+GR+1uhnNVrG0xvjo4bwjzNryvLtVyAgm3MBK3RAhPQ512gwdybvRko5xHMO8zqGW7M4MgJ4G76\nfh4qa1nCLphQLkDVWeXvgEedmnE8YLoUGtQQ8WOxfEkLmzTINozN8GiuS6/ioj2XJU+OYfTozjEM\nmfOXEeZJxRI5nG2gS7dbCBhQziCecHjVBq8m059xSmpjtV2HwRFNIkRvWUlPoewmvCEa3Rx2SVmc\niKYdigINadjSGTVR0uJWG8kyEzPgpcH2piWW7dGJ5BiOTZE4YJFIkqEjYpvqFee1gZvWkInPTHHu\n11xetcGNbVpqtCAGhxOyQnAE/4qwsvY7dKnyJWBEUmkqnHMoGtyQhi09zS/gmwibZC4c3vdVJqpB\nlOF3iB8CeeCB4yT28HvAApYzLZpEAc6CA/s4SmnPj5nZ4HnX5VUb/Bi5waxbVDPZ5h9XErTIwu30\n5TMoiKTvsGDlMry77bNPux1xXVgCdTUcMcH/+XmFjl3YLJlU0V7UdD3l00kOoo88X+p098q/scGv\n8A5xgj8qU5cs3Fd8nv1PTSVEi0CHVTbcZZLa0SlYoGTwwNGHWNU6cgKW8GeCH+AnIZ2Lp+PLSVZw\n7lHbYP4eIPgvDL8K7Fx06IvoH+wxhSVQV8MRuwp/e+101nPTM43QYe1qqZLBdon9423h1UJXcLI4\ntIusV4A7TagZmW3g5ocFq32CB3n94G3rbaFtqKtRErwi3WpMZxycnmmEjlSpOxrK9EFQE2RWDhNH\nooJDIwQuGwNtneT+lQs3zT1Q0Z6h62R7dfKiBhHGETeJPkHDGk9J1Sbs1ysXSUX2hVAIG6LT+1pi\nzpgem7G0cj2sI5aFeoQi3Ay3avataz+wC2GZ3HiaSBnDFmDMkyftQD2NklwprbovlZjwFMpkg8St\ndeRkQOoPI3/O8iBHbAIe2j3K3x4GlAJCchlnjrBbNoMHo5bQAvU1HLkyrrovlbm4pUI2RnlyM8qV\nnnhmFAKe4D6/Hnjt3OBKHeqewYGUpfR8iKcvZeFNsQA41Um7/7GxXE99DVe0cvWnizQXlL1ixFxD\n2ZXd+9JctzdD/as1T3aG4lkqKrvydnOWarzeajyHscO6XuVr0dvC28hk1r0vXYviLMv6/tjxm1mu\n8nOroyt/lz8f7fvS9DNqvmnwPWVegK681nmEnZfmajcSz9Yum8USceU+zzwN7oz9bqPBiRklZqVQ\nXPl++84/K9XdQCXePuyYh2PpBnooqu3OUr/BWhz1zsW5WalnTitZNqe1/59U/oUn/17qSe3A6o8+\nKsnwwtzaoCRxQHdNCKZr2vKQW+Ip1hRpkQKfiba829edtftceuZSJmvLtEZJWwY+uaYJPOIkn4Pl\nadoBRlJbG/hiVzJYydQ0RS0XvVdTpnUKwu70qXxzl3OJt3M7gHb/0ptvS/DVqWMZAvd0cn3vXWxh\nTevYVbOnvVmniEt4+J+LXr+pFEAzRhLKOPhKO6QNrSUV0PkA9XHNalqnoHTSLEYsrWp384nJDqDF\n0i/Acxn8jRjObjUZXWpHzMTjrWNXrZ4GDb4h6SyVdRqV/0cE7QCawhfzyBg2mbz5fEdCQUN/oNTL\nFGpxkAfUJX0MtSRInBKznQAaQmPw9clfX9oTa+A/z0AfIzN0gBaHkIGbKoHVx/l/KzeAJhMdz/2Q\nofQDGv80MOFJK10azBY3l2+FwCKUAqvKBM6UA2hcvVvNLAMiR8h8NhfII97y1jJq/uLwAL14kxNY\n9WZoHdwAWruOVyLWgAQHvsaoAd+fO6B80Oo2B3jw0mA3sMqI6gru0k4AbTWw4C2oezdKLC0HT8fr\naaC0pbe63W5gtdXNaLz/bmC1cYWWl6wOrLa8OfUNqA6s1pdveYnqwGpLmvM/ArUbFSO0YckAAAAA\nSUVORK5CYII=\n",
      "text/latex": [
       "$$\\frac{\\sqrt{2} \\left(NBEXP \\left(\\sigma \\theta + 1\\right)\\right)^{nobs}}{2 \\sqrt{\\pi} \\sigma nobs!} e^{- NBEXP \\sigma \\theta - NBEXP - \\frac{\\theta^{2}}{2 \\sigma^{2}}}$$"
      ],
      "text/plain": [
       "                                                           2   \n",
       "                                                          θ    \n",
       "                              -NBEXP⋅SIGMA⋅θ - NBEXP - ────────\n",
       "                                                              2\n",
       "                        nobs                           2⋅SIGMA \n",
       "√2⋅(NBEXP⋅(SIGMA⋅θ + 1))    ⋅ℯ                                 \n",
       "───────────────────────────────────────────────────────────────\n",
       "                        2⋅√π⋅SIGMA⋅nobs!                       "
      ]
     },
     "execution_count": 138,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "likelihood_bkg_only = likelihood.subs({ns: 0})\n",
    "simplify(likelihood_bkg_only)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABVgAAAAyBAMAAABbg615AAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMA74lUMhCrmXbNZiK7\nRN38cqFhAAAACXBIWXMAAA7EAAAOxAGVKw4bAAARHklEQVR4Ae1de3BcVRn/srs3+8puSnmoWMna\nGWwRtYuo+AdDlgFk+IdsGW3LCM2iQyjlkQDaAlK74EjDQxtHKUSULMigCDNdRwW0TLMdLA9BslDH\nGRw1axUEeSQtNH2lXb/vvO65r30krZO098zk3nO+7zu/+93f+e453717OoWO6k7wi8/AjGdgUXUH\nHHvueTPeT99BnwH48rmnwXE+Dz4Ds4OBs/1gnR0D5XsJoAXraVmfkMPMQPD4w3yBIxteBatxRa8f\nrId5rJes2nWYr3Bkw6tgBVjvB+vhHuuIH6zTodgzWGOF6cAefX2XNHLLtmD9ZSN9fBsjJzjwDNY7\nfJKaYiAx0IC5NVgjjXRpAPWIN3lN3KFXsCYO06vA2uoRUcadAXK7U+SQWIP1lUObd4WPCGarVQdr\nEMlxmVewLs44+xwCSbx0CEBmJkS0UN8vS7Aah/gL99P1rz9bLX7MHfcK1p7p3VfJo3tr2UMx+8Xx\nBl6eLMHampvSTbd59brGSzH75cN8DfII1tD+ad1h0Gvgfjct2JnduYHn2xKsXcUp3U+nR/IQH5gS\n3KzoJJ5rj2CN5qZ1Exd5BeuaacHO7M5dHlGkeW0J1g2aovGqscjjMuFy4yCzzTLE48kjWIeL07qf\n73oEqzHQFOyC59TAPLROVZuCcBqXnCJTUlNpmrnWonlXsS7UgzU+tcUr0u1BxMX6herWA7d8U9rE\nFqqqFE3x7JmhEF5NZb0L3scMPIL12/V619SHPu8RrOF0zX42ZaDYOiBEsYxRsmmn2PTMUAivprLe\nBWNz6lmAHqzhSl1zN4NzvIL1BjdrT9mpcKXUPQyXyOo0z14ZCoOtqax34QeZgRmsqzo2FFSfD6va\nVCoBwyNYlzeFliyq1LmtCDJum4JwGntmKGRaU+nEskriB6xtZyt82d6VStqSV9VmKimPYDUqzaDA\n26BS5hch2lRXT2PPDIV61FR6QkrFcJlqZrBKOZ3jH+ityAgOw1k7tgLEeid6rv4TQGRksqdnDMPn\n8eqbcHH1GZT0Tva8IJeTx7yC9ZMIa6I9Uu3peQm9kGixkYN5Y2gyDXB79dmr780my2r8wx8L5XSX\n0MXVL2WaAZO9XTKU86UOrEpjS0VpGqnMbcRI2QxnqHp3tQzBsXlZRS1Knu25di9yLdmwUBsregRr\noOCGJql1DNS7wMcfO11X/hIe9fLkLX/Aprdr7qMOLhlKqE/i2pSXnCQVDZ1bUmTmHqyxcdKp0jaJ\n1Xuo2VICWJsBaB9ANrEOvQWIXY5nJtmSpRoYeY9g5U+/iTYPIEyPhURr3QfwCMPAwWo7XluWjbHn\nGLR5WAIBmsmaAGN93TKUdySqXenxmuncO8U3rKmFVeLVPG/k2pESQJTuWVILeO/wOP4pNogeQS0s\nA49g5YuWA01S6xyoTj5YAMsn+/BqWjFysKiMbQeYcs1t1AFcMpRkSeLale6vl05quYTng+7BGq7I\na7DztpEiQIqqXQWAUby30TzGJN1j1wDw57IrBdCJUizB1at3F1gNvshP4hhhUoUWwthM7EadRAvt\nhVCObEP4L22CewCSeWpRefAguqAXfEKOwXYTYKy3maEo1yL4zPBiKnm7PS01dA7mecuxd0puWFvP\n9Q0e+WtB4kx8zwpQF0ltAu8bluKfYkOjFs5YvfZmsjbd4S24k85ONEmtc6DeFB0h1vG2rPJzaxZa\n5riBKdfcRh0jRD5HkimAW0sSWimF4H2pYGcxGCa1EkJIIuNk5h6s5gPBoE7G+EywgVtfBNiYB1iU\noQkUldEP4jk84Z6tMv6lWBXXaJmzqohgiiQ7KrQIjlOMhkahDWWeYCYRnDQjGLALWAsPgdRpONHq\nBZeR3gxAE2Cst5mhKNd+/xOJayq5ZDQrNXSW/Flek7gB37DWieQ0Xvh8HrgUn9Zt1EtSG6C1ZhlJ\nJBsWamFYOKXcIVOACh2caIpa+0AF+qgDlfnBYyz3CS0VaEUnnGDKNTsYwzEzFOVa/NESU+E429KX\nhIwQrpeDod5AFQSXBDFUvIK1Xd0Jw+pLjkOA3dAWjNGTsDZEM20RlcEJHl4wgo3eMjPHVWwizWvS\nCTCo3c+ECq1tDsDi11Gm0DqvyjGTtnFcBQ9ArI+18LAcjH/JOj8PZVmwNg7GumkZinItJYNVU/Jr\nfMZySRt/uo4Ha5e4aV3jXaclBdMgGMMnjmqS2ihyDXGSSDYs1MZGvk467dmBLLYieZI50RS19oH6\nOdlTwdxsudXx6C4WrE4w5ZodjAGZGYpiKhAtMRU9fXLa5YJARSj4SQ6GV7CG2ETsPrO25y1QfbF9\nyAOVCYALL8fzJEQK1Abo2MrPmGnFkf+lHdWqCjFQaUB8O1lVmKlCa89B6A1iWqG1TFKTJ7HDffAq\n/FkgPg3wHabRD2No3DiYcV21utslQ4kVMVhNZeITm/pX4FUW3l6GDf034ATT/zC/qBoCxahyhgdr\nd0EJGqjsYTZh2FgARpiktr3El3TMYgUbglrmowms3NmcQiFftJxoilrbQCVSiQy/63gJImkTlteS\nOPE5wUzX3EbdJUN5jIKVhQRTvnrVU7fnidFrILpiHV504Z0Zfr16wWpgxug1s1ppxwzgIM+rQhP9\nd9MMl9jb01Hkl1mUZ2fKMs9KQeg/297ZxBXsKJ346w+wGUuRzETrWtH/MtFkogX4CEJXCoLbIfTD\n07cKxGgR6BmxlAResnEwuHDNWS+iz44MZRlgsJrKcNvz0F2EzWVjv7E3i7WfwnP8sio6vILV9ohb\nnHU0OP2YAbRXcJHCDFVQC6MbLmNrHg4qZ0NQK32USNKdxKNzUPQQEzvQTGrBOlAXnH5mViDeD5ik\n2spoipITb9esYHLU7RmKkcdg1ZR9L+WN/RB/Fwnuvgdr0XJbH7+wjBNFrbw7mXMxKtxn1uGC7jxm\nAGvLfycJZZmLc0jjB+rbXC9RhZqJ+avvQmUaNrK2OEgnoLOATz/jxETDLM04AQ1NtNYdvN/6ry24\nsQgt1eq4QDSuv60oINUpmcOuyrV6YMaxEBvHvo4MJYXBqinPSQ7AcAbmAmzHnwjaU/FxOI9f0s6f\ncgTzTXZn7SlNVK8aZ3MFZgDhA4kyGktqKYEf4J0FG4Ja6aMEVu4ECel1JnagmdSCdaBwAZR3vXTh\nqRJTnd/AmgNMc80KJkbdkaEEAIPVVIbSb9FH0ZYUJFOdWdgJwxBixJkrsGew4sQkZ1ZtJ+Q4SVn2\ndfF7VH5E6wEMl04medscZBXTB0xewhAiSfgJlk1AG0+Yb8JXe3F9gFMGB/85OLiBzChph8+ymomG\nWRqsT+to3+rFJpYRhQEWROMBcumtDDNajcfGwcIDoCX1pmvBDAarpsx2l/F5C+PtTKC0O2XsWENX\nCwwOPvCRwcEc1Tmjui8iWPtIay3Kql3jmEx4FoYZQHwSnziTWkrg+yBDJoINQa3mo9UdwHfNYIo6\nONE8B4qsLYighps0dCc1XKs16hpTn6NgNQfQMPA1ZBzWZqE9dxMYO6H7Q2XywxwMQa0GIWdWluG7\nz6z0GcUsAaRyPxsIzDLZWzomLwZ/Y/1GaJIZ4opBBR/y77OKOKiZld7v72FCEw2zNPYBQaGFM10l\nZnOQHdnBjqg0EfKocbD2AkQqqjNWhGsXAQarruwE/CWypQTxPdE+GE3D36ol3i2Y52fJn2jRyTtY\nNSNrlb7bYUnhk/kkVSS1mMBD1sijRLIhqNV9pA7KHfoQE82SyImmqAXbQJG1HZFkvKyjU8rbNRuY\nfYyEa0aKglVTBscBP5jOw7kvfy8k3ofECA8eNRgmteruxFxbI1iHC+StLNswN9ydolZXgb2lsw+j\n7B02UIaxItOwYKZfvu6lpiwqWI1JCDIT+k7D0dhw9Zb5Z1aG9mlIsnQtxFxjGA5EiQxnQLzYBBgu\n7Rh8WhGufaGnZ8cKXYkzwT6jOwWB/TiYa7O/hmUihbTzp2HxYMVOjReeBmDODaPXUy9JLfuWyeZa\nyUYXd1v3kTood6C7hGRQcaKN5lFMi6JtoFCCIWOnhEnxgCO1yQVMuWYDc4yRcC3Y09N7wl1aSOBU\n3l74Cq7oi4r4xr4rnjbWF/k1ZZx4pQE4DWNxn1m78xyDH+lmx8pUJ2ya/DCvglAOBd/Dz66MS1rO\n8aV6F4SOp4os0gn83asYxc5YFBplaWxelmiRMsRYmFKmIYoDUSrwh7xIsQmw7iycnZWd6Wy6No8+\nqyjlXNzOEsA4bUm1l40TAd8Mx3k3FR2KUS7HIw9W6yOulO4VYw/JKSqTB6gmqWXfMul3FsUGp9bi\nI3VQ7kB4P9xPEhc0Sa19oJi5ftdMIA/nA9zlAqZcqzfqpmstJX0AW1LQmU3hVDsXDsDiQnQOtItL\nysFQ1CoILuGPtnuwWl8VKoi4kaHiyCXmIb9j2Holg7NsHy4mLFsdQimWZ/A1mlXEQTqBj3HqL1xW\nwRNDQ1/hla3YkmhXYL1KOPihTxU7olRsnr/gWoAKNhsDa0vjW6heTNcmQVPidoTAwMPhHH4DaC3g\ny+qL4hOmFh2KUQXHg7UrjQJjr5LWrrC5Iok2sXEylNQmsZU4DgWKDUGt5iPZa8Ea35nIMZETTVLr\nGCiytyMyDDzEt89fVWEfw9xdc4DZx0hFGmAaoilxO8IpscwqQFpXIr9tRfHl3Zw5FLUKgktqfWdt\nyUnHkf0tBwsYliS4knaepMGg88s4KBdUj4PE0MEcwE3V58kAlvbTxGkWMyLaxtnTb6K1DuGmmGcQ\nX6CFhnZk4LUqvoeGhybzCsKOKBVD1erBZsCMG2+j58AsyrVPVVdqykQF4ivTMH9dGYwbbwZ46k75\ngVfyZ9k7xfDEhrXRMrXes17GvKCtdgK2L+n4Lx5z+CeohSXEytgAmGxIajUf0V4PVhjhi5YDTVLr\nMlCEYEckGZUovgtWHGDKtfqjLpmCwMhEQRvAfwA88TMILFyDgbKuAPH+X6T5FeUyZ1IrIYQkwVYf\nMbPyfTaiJz5zA6o6vYqKCEhM5KYHdah7m641iCz58zTfyMI0WHQ34NvETN0LZnVKNc2d0Z4pIRy2\nTpprDV7DMRg2iBhba3mwGjm+z0ZCB+bI2jTPXzX7j2XM+kyoaa415o7hEYWq98dZDT9QuBaxTUzp\nuLVqNl3R3EniS8xMKpprDbrlGAwbBA9IHqytfJ+NQsbP4Ye83H3IEWcaIGYwWJZ5uHU53yamtJ1s\nHlbN6VQipen0ng19kznykgdrS4VtXVBui98CVduvNMCAIK3gYXoS3yamtN1ehsrCrygGOFk8WMU+\nG6WDZ82qX2uQAb6NzdMYX+p7M5q2Lac1/GptBvjuSx6saIn7bLRtPQ/W7utrXRjAT951ylhW25UW\nO1TvBXUuekSo8ZcaLCpYR1P6tp7u9BFxj//Xm2gp1LkcbhOTW5DQ0jB/+ajTz1fDiYwDFaxvyE04\nTBypO0v4DNoZuNUusLfxLUHflXaHXe+3vRiI8FVIBmu4z7oJh35D8UtTDPCnv0aX1fjjibYrrSVT\nw9ZX6QyIRUsG6zrbJpxHdVu/3gADdRcj2ia2T9uVFvNXrwZoZSZiERLBGsR9NsP6JpxovlEg344z\n8KtiHSZwm9gmbQsSwH11OvhqwUD8o7wigvV83GfTrW09AuPfPlVNMWD/B42OzrRN7LeWXWnLMw4j\nX+DGgPwHjTxY2T4b6yaczfUmCjfUo1jWVqhz82ybmLYFCXc3vV2ni6/mDPxREMGDle2zsW7CiWd8\nqpph4NJ6xrRNzLYr7Tf1+vh6YsAoCx5EGuCz4jMw8xnwg3Xmj5HvoWDAD1Y/FGYNA2f7/4X7rBmr\no9xR+i/cO6rsXwMd5Uz4tz/jGVhU3fE/iSX1AwD7ITQAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left [ \\frac{1}{2 \\sigma} \\left(- NBEXP \\sigma^{4} + \\sqrt{NBEXP^{2} \\sigma^{8} - 2 NBEXP \\sigma^{4} + 4 \\sigma^{4} nobs + 1} - 1\\right), \\quad - \\frac{1}{2 \\sigma} \\left(NBEXP \\sigma^{4} + \\sqrt{NBEXP^{2} \\sigma^{8} - 2 NBEXP \\sigma^{4} + 4 \\sigma^{4} nobs + 1} + 1\\right)\\right ]$$"
      ],
      "text/plain": [
       "⎡                    ____________________________________________________     \n",
       "⎢             4     ╱      2      8                4          4               \n",
       "⎢- NBEXP⋅SIGMA  + ╲╱  NBEXP ⋅SIGMA  - 2⋅NBEXP⋅SIGMA  + 4⋅SIGMA ⋅nobs + 1  - 1 \n",
       "⎢────────────────────────────────────────────────────────────────────────────,\n",
       "⎣                                  2⋅SIGMA                                    \n",
       "\n",
       "  ⎛                  ____________________________________________________    ⎞\n",
       "  ⎜           4     ╱      2      8                4          4              ⎟\n",
       " -⎝NBEXP⋅SIGMA  + ╲╱  NBEXP ⋅SIGMA  - 2⋅NBEXP⋅SIGMA  + 4⋅SIGMA ⋅nobs + 1  + 1⎠\n",
       " ─────────────────────────────────────────────────────────────────────────────\n",
       "                                    2⋅SIGMA                                   \n",
       "\n",
       " ⎤\n",
       " ⎥\n",
       " ⎥\n",
       "─⎥\n",
       " ⎦"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_hat_hat_solutions = solve(diff(likelihood_bkg_only, theta), theta)\n",
    "theta_hat_hat = theta_hat_hat_solutions[0]\n",
    "simplify(theta_hat_hat_solutions)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUYAAAFACAYAAADAjtXtAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4FNeZ7/Hv261dQruEAAkkIRazGTAgNoNtHBvsxMSx\nHeNJ5tqxE2fzzGSy3DgzdzxOZpKb5OaOM7lxZuKJs9hJvATHCbGJHeOdTeyYHQlJgMSiFS1AS+ru\nc//oEiO3BWpJra5e3s/z6FGpqrrrLVr6UXXOqSoxxqCUUuq/OewuQCmlwo0Go1JK+dFgVEopPxqM\nSinlR4NRKaX8aDAqpZQfDUallPKjwaiUUn40GJVSyk+c3QX4y83NNcXFxXaXoaLckSNHAJgyZYrN\nlahQ2blzZ5MxJi+QdcMuGIuLi9mxY4fdZagod9111wHw1ltv2VqHCh0ROR7ounoqrZRSfjQYlVLK\njwajUkr50WBUSik/GoxKKeVHg1EppfxoMCqllB8NRqWU8qPBqJRSfjQYVUwxxrC5qomGji4aO7r4\n095TXOh2212WCjMBXRIoIiuBfwecwM+MMd/1W54IPAVcAzQDdxtjakWkGDgEHLFW3WqM+VxwSlcq\ncF6v4ZUDZ3j8zSoOnGrnTGMnAH/zzG5yUhP4zLJS/nrhBFITw+4qWWWDAY8YRcQJPA6sAqYB94jI\nNL/VHgBajTFlwGPA9/osO2aMmW19aSiqkHN7vHzld3v5wm92caHbw/fvmMXc8VnMGZ/Fsw8uZPq4\nDL7758Pc+qN3OXXuot3lqjAQyKn0AqDKGFNtjOkGngVW+62zGviVNb0WWCEiErwylRqaLreHL/52\nFy/urufLH5rMhi8v5+Pzi0iIc5AY52BhaQ5P3b+A3366nObObj7+0y2cbLlgd9nKZoEE4zjgZJ+f\n66x5/a5jjHEDbUCOtaxERHaLyNsicm1/GxCRB0Vkh4jsaGxsHNQOKHU5Hq/hC7/exasHzvLIh6fx\ntysm4XT0///14rJcfvOZcjpcbu7+6RbqWjUcY9lId76cBsYbY+YAXwZ+KyLp/isZY54wxswzxszL\nywvodmlKDei/3q3m9cMNPPqRady/tGTA9WcVZvLMZxbS4XLz98/tweM1IahShaNAgrEeKOrzc6E1\nr991RCQOyACajTFdxphmAGPMTuAYMHm4RSs1kP31bfzfvxxh1YwC7l1cHPDrpo1N55urp7O9tpWf\nvnNs5ApUYS2QYNwOTBKREhFJANYA6/zWWQfca03fCbxhjDEikmd13iAipcAkoDo4pSvVP1ePhy89\nt4eslAS+c/tMBtvcffuccdwys4DHXjvK/vq2EapShbMBg9FqM3wIeBXf0JvnjTEHRORbInKbtdqT\nQI6IVOE7ZX7Ymr8MeE9E9uDrlPmcMaYl2DuhVF8/er2SqoZOfnDX1WSlJgz69SLCtz86k6yUBL7y\n/F7cHu8IVKnCWUCDtowx64H1fvMe6TPtAu7q53UvAC8Ms0alAnamzcWTG2v46OyxLJs89PbqrNQE\nvnnbdD7/m138fnc9H59XNPCLVNTQK19UVPnhhqN4jeErNw3/IVcrZxRwdVEmP3ztKK4eTxCqU5FC\ng1FFjaqGTp7fcZJPlE+gKDtl2O8nInx95RROtbl4ekvAz1FSUUCDUUWNH7x6hOR4Jw/dUBa091w8\nMZdlk/N4/K0q2l09QXtfFd40GFVUOHymnVcOnOHT15aSm5YY1Pf+nzdP4dyFHn6+sSao76vClwaj\nigq/2FhLUryDTy0pDvp7zxiXwfVT8vj11uN0ubWtMRZoMKqI19zZxYt76rljbiGZKYMfnhOI+5eW\n0NTZzZ/2nh6R91fhRYNRRbzfVpyg2+0dkaPFXkvLcpmUn8bPN9ZgjF4qGO00GFVE63Z7eWrrcZZP\nzqMsf9SIbUdEuH9pCQdPt1NRo9coRDsNRhXRXt53isaOroBuEjFct88ZR1ZKvHbCxAANRhXRfr31\nBKV5qSyblDvi20qKd3LPgvG8duis3tA2ymkwqohV3djJzuOt3D2vaNA3ihiqu+cXYQy8uNv/BlMq\nmmgwqoj1wq46HOI7xQ2VCTmpLCjJZu3OOu2EiWIajCoiebyG3++qZ/nkPPLTk0K67TuvKaSm6Ty7\nTrSGdLsqdDQYVUTacqyZ020u7rimMOTbvmXmGJLjnazdWRfybavQ0GBUEWntzpOkJ8Vx41WjQ77t\ntMQ4Vs0s4KW9p7nYrVfCRCMNRhVxOlw9vHLgDLfNHktSvNOWGu66poiOLjd/OXjGlu2rkaXBqCLO\nXw6cxdXj5WNzQ38a3au8JJtxmcn8QXuno5IGo4o46/edZmxGEnOKMm2rweEQbp01ho1VTbRd1NuR\nRRsNRhVR2l09vFvZxKqZY0I2dvFybpk5hh6P4bWDZ22tQwWfBqOKKK8fOku3x8stM8fYXQpXF2Yw\nLjOZ9fv0jjvRRoNRRZSX3zvDGJtPo3uJCLfMLODdykY9nY4yGowqYnS4eninspFVM8bgcNh7Gt2r\n93R6g55ORxUNRhUxXj/UQLfby62zCuwu5ZLZRZl6Oh2FNBhVxFi/7zQF6UnMKcqyu5RLRIRVMwp4\nt7JJH5YVRTQYVURw9Xh4p7KRm6aPDpvT6F6rZhbQ7fHy9pFGu0tRQaLBqCLC5mNNuHq8rLDhEsCB\nzC7KIjs1gQ2HtJ0xWmgwqoiw4VADqQlOFpZm213KBzgdwg1T83nzcAM9Hq/d5agg0GBUYc8YwxuH\nGrh2Uh6JcfZcGz2QG6/Kp93lZket3oosGmgwqrB34FQ7Z9pdrLgq3+5SLuvaSXkkOB16Oh0lNBhV\n2Ntw6CwicP3U8A3G1MQ4Fk3MYcOhs3pn7yigwajC3uuHGphTlEluWqLdpVzRjdNGc7z5AscaO+0u\nRQ2TBqMKa2fbXeyrbwvL3mh/K6wj2g2HGmyuRA2XBqMKa28e9oVMOLcv9hqbmcy0Mem8ru2MEU+D\nUYW1dyobKUhPYsroUXaXEpDrp+ax68Q5vQomwmkwqrDl9nh5t7KJZZNzbb/3YqCWT87H4zVsrmqy\nuxQ1DBqMKmztrTtHh8vN8snhfxrda874TEYlxvH2Ub08MJJpMKqw9faRRhwCS8ty7S4lYPFOB4vL\ncnjnaJMO24lgGowqbL1d2cTVRZlkpMTbXcqgLJ+cT/25izpsJ4JpMKqw1Hq+m/fqzrF8cp7dpQza\nssm+I9y39G47ESugYBSRlSJyRESqROThfpYnishz1vIKESn2Wz5eRDpF5KvBKVtFu3ermjAGlkVg\nMBZmpVCWn6btjBFswGAUESfwOLAKmAbcIyLT/FZ7AGg1xpQBjwHf81v+b8Cfh1+uihXvHG0kIzme\nqwvtf7bLUCyfnEdFTQsXuz12l6KGIJAjxgVAlTGm2hjTDTwLrPZbZzXwK2t6LbBCrPEVIvJRoAY4\nEJySVbQzxvDO0UaWluXiDLOb0gZq+eQ8ut1eKmqa7S5FDUEgwTgOONnn5zprXr/rGGPcQBuQIyJp\nwNeBbw6/VBUrqho6aejoYumkyOmN9regJJuEOAebdDxjRBrpzpdHgceMMVfsnhORB0Vkh4jsaGzU\ndplY1xsmkTRMx19SvJN5E7LYWKVHjJEokGCsB4r6/Fxozet3HRGJAzKAZqAc+L6I1AJfAv5BRB7y\n34Ax5gljzDxjzLy8vMhrbFfBtelYM0XZyRRlp9hdyrAsKcvl0Ol2mjq77C5FDVIgwbgdmCQiJSKS\nAKwB1vmtsw6415q+E3jD+FxrjCk2xhQDPwS+Y4z5cZBqV1HI7fGytbqZJRMj92ixV+8R7+ZjetQY\naQYMRqvN8CHgVeAQ8Lwx5oCIfEtEbrNWexJfm2IV8GXgA0N6lArEvvo2OlxulkTwaXSvGeMySE+K\nY1OltjNGmrhAVjLGrAfW+817pM+0C7hrgPd4dAj1qRjTe3S1eGKOzZUMn9MhLJqYw8Yq3+WBkXIj\nDKVXvqgws6mqiakFo8gJ87t1B2ppWS715y5yvPmC3aWoQdBgVGHD1eNhx/HWiO6N9tfbJLBRh+1E\nFA1GFTZ21LbS7fZGRftir5LcVMZmJOl4xgijwajCxuZjTcQ5hPkl2XaXEjQiwuKyXLZUN+P16m3I\nIoUGowobW6qbmVWYQVpiQH2CEWPxxBzOXejh0Jl2u0tRAdJgVGHhfJeb9+raWBQFvdH+evdpi45n\njBgajCosbK9tweM1LCyNvmAck5FMcU4KW6s1GCOFBqMKC1urW4h3CtdMyLK7lBGxaGIOFTW+8Ffh\nT4NRhYUt1c1cXZhJSkJ0tS/2WliaQ4fLzYFTbXaXogKgwahs1+HqYX99dLYv9lpUqu2MkUSDUdlu\nR20rHq+5FB7RKD89ibL8NLZoO2NE0GBUttta3UyC08HcKG1f7LWoNIdtNS30eLx2l6IGoMGobLel\nupnZ4zNJinfaXcqIWjQxhwvdHt6r03bGcKfBqGzV274YjcN0/PXuow7bCX8ajMpWO2pb8RpYGEWX\nAV5OdmoCk0enUVHTYncpagAajMpWW2uaiXcKc8ZHd/tir4WlOeyo1XbGcKfBqGy1raaFqwszSU6I\n7vbFXuUlvnbG/fXazhjONBiVbS50u9lX18aCGDiN7tW7r3o6Hd40GJVtdh5vxe01lMdAx0uvvFGJ\nTMxLpUI7YMKaBqOyTUV1C05H9F4ffTnlpTlsr23Fre2MYUuDUdlmW00LM8amR939FweysDSHzi43\nB0/r/RnDlQajsoWrx8Oek+di6jS6V+/QpIpqbWcMVxqMyha7T5yj2+NlQXHsdLz0yk9PoiQ3VQd6\nhzENRmWLbTUtiMD8GAxGgPKSbLbV6v0Zw5UGo7LFttpmphakk5ESb3cptigvzabD5eawPgcmLGkw\nqpDr8XjZdfwc5TE0ftFfeYmvbXWbjmcMSxqMKuT217dxsccTs6fRAGMzkynMStZgDFMajCrkesNg\nfklsjV/0t6Akm201LRij7YzhRoNRhdy2mhZKc1PJH5Vkdym2Ki/Jpvl8N8caO+0uRfnRYFQh5fUa\ntte2xPRpdK/edka9bjr8aDCqkDpytoN2lzumbhxxORNyUsgflajtjGFIg1GF1PZaXwhoMIKIsKAk\nm4pqbWcMNxqMKqQqaloYk5FEYVay3aWEhfKSbM60uzjZctHuUlQfGowqZIwxbK9pYUFJNiJidzlh\nofda8YoavTwwnGgwqpA53nyBho4u7XjpoywvjayUeO2ACTMajCpktmn74gc4HML84uxLba8qPGgw\nqpDZXtNCVko8ZXlpdpcSVhaUZHO8+QJn2lx2l6IsGowqZLbVtjCvOBuHQ9sX+7p03bQeNYYNDUYV\nEg3tLo43X4jJ+y8O5Koxo0hLjGObdsCEDQ1GFRK9R0PztX3xA+KcDq6ZkKUDvcNIQMEoIitF5IiI\nVInIw/0sTxSR56zlFSJSbM1fICJ7rK+9InJ7cMtXkWJ7TQvJ8U6mj023u5SwtKAkm6NnO2k93213\nKYoAglFEnMDjwCpgGnCPiEzzW+0BoNUYUwY8BnzPmr8fmGeMmQ2sBH4qIrH15CMFwLbaVq6ZkEW8\nU09S+tPbU6+90+EhkN/SBUCVMabaGNMNPAus9ltnNfAra3otsEJExBhzwRjjtuYnAXrdUwxqu9jD\n4TPtOn7xCmYVZpAQ59DT6TARSDCOA072+bnOmtfvOlYQtgE5ACJSLiIHgH3A5/oEpYoRu463Yoze\nf/FKEuOczCnK1J7pMDHi5zXGmApjzHRgPvANEfnATfhE5EER2SEiOxobG0e6JBVi22pbiHcKc4o0\nGK+kvCSb/fVtdLh67C4l5gUSjPVAUZ+fC615/a5jtSFmAO8be2CMOQR0AjP8N2CMecIYM88YMy8v\nLy/w6lVE2F7TwoxxGSQnOO0uJawtKMnBa2DXiXN2lxLzAgnG7cAkESkRkQRgDbDOb511wL3W9J3A\nG8YYY70mDkBEJgBTgdqgVK4igqvHw966c3oZYADmTsgkziE6njEMDNhDbIxxi8hDwKuAE/i5MeaA\niHwL2GGMWQc8CTwtIlVAC77wBFgKPCwiPYAX+IIxpmkkdkSFpz0nz9HjMTqwOwApCXHMGJehHTBh\nIKChM8aY9cB6v3mP9Jl2AXf187qngaeHWaOKYNtqWhCBeRM0GANRXpLNLzbV4urxkBSvTQ920UFl\nakRtr21hyuhRZKTE211KRCgvzabb42W3tjPaSoNRjRi3x8vO462Ua/tiwK6ZkI0IejptMw1GNWIO\nnGrnQrdHr48ehIzkeK4qSGdbrXbA2EmDUY2Y3qMe7XgZnAUl2ew83kq322t3KTFLg1GNmG21LRTn\npJCf/oEx/eoKykuycfV42X+qze5SYpYGoxoRXq9he22LXh89BL1NDxXV2s5oFw1GNSKqGjs5d6FH\nB3YPQW5aImX5afrkQBtpMKoR0fvUOw3GoVlQks3O2lbcHm1ntIMGoxoR22paKEhPYnx2it2lRKTy\nkmw6utwcPN1udykxSYNRBZ0xhorqZhaUZCOiD74aioWlvgdkaTujPTQYVdAdb75AQ0cX5aV6Gj1U\no9OTKM5J0XZGm2gwqqDrHb+oV7wMT3lJDttqWvB49cb3oabBqIJua00z2akJTMxLs7uUiFZemk27\ny83hM9rOGGoajCrottW0sKBY2xeHq1zbGW2jwaiCqv7cRepaL2r7YhCMy0ymMCtZ2xltoMGogqr3\n7tM6fjE4etsZvdrOGFIajCqottW0MCopjqkF6XaXEhXKS7NpvdBDZUOn3aXEFA1GFVQV1b7ro50O\nbV8MhkVWO+PWaj2dDiUNRhU0De0uqpvOs1DbF4OmMCuZcZnJGowhpsGogmarNX6x96oNNXwiwsLS\nHLZWN2s7YwhpMKqgqahuJi0xjmljtH0xmBZa7YxHGzrsLiVmaDCqoNla3cz84izinPprFUy9R+Bb\nj+npdKjob7AKioYOF8caz+tp9Agoyk6hMCuZLdrOGDIajCooLl0frcE4IhaV5lCh4xlDRoNRBcXW\n6mZSE5zMGKvtiyNhYWkO5y70cPiMtjOGggajCoqK6hbml2Rr++IIWThRxzOGkv4Wq2Fr6uyisqGT\n8hI9jR4p4zKTGZ+dou2MIaLBqIat9+4vOrB7ZC0qzaGiulnvzxgCGoxq2LZUN/naF8dl2F1KVFtc\nlkO7y80Bfd70iNNgVMO2+Zjv+S7x2r44ohZZ7YybdTzjiNPfZDUsZ9tdVDeeZ/HEXLtLiXr5o5KY\nlJ/Gpqomu0uJehqMali2WEcvvUczamQtKctle20L3W593vRI0mBUw7L5WBMZyfF6fXSILJqYg6vH\ny+4TrXaXEtU0GNWwbD7WzMLSbBx6/8WQWFiag0O0nXGkaTCqITvZcoG61ovavhhCGcnxzBiXweZj\n2s44kjQY1ZD1/nEu1vbFkFo8MZfdJ85xodttdylRS4NRDdmWY83kpiVSlq/Pjw6lxRNzcHvNpRt3\nqODTYFRDYoy51L6oz48OrfnF2SQ4HWys1NPpkaLBqIakqqGTho4ulpZp+2KoJSc4mVecxUYdzzhi\nNBjVkLxrHa0s0WC0xdJJuRw+00FDh8vuUqJSQMEoIitF5IiIVInIw/0sTxSR56zlFSJSbM3/kIjs\nFJF91vcbglu+ssumqiaKc1Ioyk6xu5SYdG1ZHoBeBTNCBgxGEXECjwOrgGnAPSIyzW+1B4BWY0wZ\n8BjwPWt+E/ARY8xM4F7g6WAVruzT4/GytbpZjxZtNH1sOlkp8ZeO3FVwBXLEuACoMsZUG2O6gWeB\n1X7rrAZ+ZU2vBVaIiBhjdhtjTlnzDwDJIpIYjMKVffacPMf5bg/XTtJgtIvDISwuy2VjZRPG6G3I\ngi2QYBwHnOzzc501r991jDFuoA3wH9x2B7DLGNM1tFJVuNhY2YRDYFGpBqOdlk3KpaHDd5NgFVwh\n6XwRken4Tq8/e5nlD4rIDhHZ0djYGIqS1DBsqmpiZmEmGSnxdpcS05ZO8rUz6ul08AUSjPVAUZ+f\nC615/a4jInFABtBs/VwIvAj8D2PMsf42YIx5whgzzxgzLy8vb3B7oEKqw9XD7pPnWFqmV7vYbVxm\nMqW5qbxbqQcTwRZIMG4HJolIiYgkAGuAdX7rrMPXuQJwJ/CGMcaISCbwMvCwMWZTsIpW9qmobsHj\nNdrxEiaunZRLRXULrh6P3aVElQGD0WozfAh4FTgEPG+MOSAi3xKR26zVngRyRKQK+DLQO6TnIaAM\neERE9lhf+UHfCxUy71Y2khTv4JoJWXaXooDlU/K42ONhR63ehiyY4gJZyRizHljvN++RPtMu4K5+\nXvevwL8Os0YVRt6pbGJRaQ6JcU67S1H4bkOW4HTw9tEGluoogaDRK19UwI43n6em6TzLJ2s7cLhI\nSYhjQUk2bx3RdsZg0mBUAXvnqO+Pb/kUbQ0JJ9dNyaOyoZP6cxftLiVqaDCqgL19tImi7GSKc/Qy\nwHDSewT/th41Bo0GowpIt9vL5mNNLJ+cp7cZCzNl+WmMzUji7aMNdpcSNTQYVUB2HG/hQreHZZO0\nfTHciAjLp+SzqaqZHo8+PTAYNBhVQN452kScdX2uCj/LJ+fR2eVm53EdthMMGowqIG8fbWRecRZp\niQGN8FIhtqQshziH8OYRPZ0OBg1GNaCz7S4OnW5nmQ7TCVujkuJZUJLNG4c0GINBg1EN6M3Dvj+2\nG6bqMJ1wtuKq0VQ2dHKi+YLdpUQ8DUY1oNcPNzAuM5kpo0fZXYq6ghXWf1yvHz5rcyWRT4NRXZGr\nx8PGyiZumJqvw3TCXHFuKhPzUnnjsJ5OD5cGo7qiipoWLvZ4uOEqPY2OBCuuGs3W6mY6XD12lxLR\nNBjVFb1x6CzJ8U4Wler9FyPBiqn59HiMPnN6mDQY1WUZY3j9cANLynJIite76USCayZkkZEcz+t6\nOj0sGozqsiobOqlrvcgNU0fbXYoKUJzTwXVT8njzcAMerz4ka6g0GNVlvX5Ih+lEohuvGk3z+W52\nndCrYIZKg1Fd1oZDZ5k+Np2CjCS7S1GDcN2UPBKcDl7df8buUiKWBqPqV0OHi10nWrl5eoHdpahB\nGpUUz5KyHF45cEafOT1EGoyqXxsONmAM3DRd2xcj0c3TC6hrvcjB0+12lxKRNBhVv149cIYJOSl6\ntUuEunHaaBwCrx7Qq2CGQoNRfUC7q4fNx5q4eXqBXu0SoXLTEplXnK3tjEOkwag+4M3DDfR4DDdN\n09PoSHbz9AKOnO2gtum83aVEHA1G9QF/OXiW3LRE5o7XZ0dHsput9uFXDuhR42BpMKr3cfV4eOtw\nAx+aNhqHQ0+jI1lhVgqzCjNYv++03aVEHA1G9T4bK5s43+25dLShItutM8fwXl0bx5v1dHowNBjV\n+7z03ikykuNZos92iQq3zhoDwEvv6VHjYGgwqktcPR5eO3iWldMLiHfqr0Y0KMxKYc74TA3GQdLf\nfnXJW0caOd/t4cNXj7G7FBVEH541lkOn2znW2Gl3KRFDg1Fd8vK+02SnJui9F6PMrTPHIAIv61Fj\nwDQYFQAXuz28fugsK2cUEKen0VGlICOJ+ROyeem9U3aXEjH0L0AB8OaRBi50e/jwTD2NjkYfvnoM\nR892cviMXjsdCA1GBfh6o3PTEijX0+iodOvMMcQ5hBd31dtdSkTQYFS0Xexhw6EGbp05BqcO6o5K\nOWmJXDclnxd31+udvQOgwahYv+803W4vH5tbaHcpagTdMXccDR1dbKrSB2UNRINR8eKuekrzUplV\nmGF3KWoE3XBVPulJcfx+V53dpYQ9DcYYd7LlAttqW7hjbqHeYizKJcY5+cjVY3nlwBk6u9x2lxPW\nNBhj3Iu7fY3xq2ePtbkSFQofmzsOV4+XP+uNJa5IgzGGGWN4cXc9C0uzKcxKsbscFQJzx2dRnJPC\n2p16On0lGowxbPfJc9Q0nedjc7TTJVaICHfNK6KipoVqvUTwsjQYY9hz206SHO9k1Ux9EmAsuWte\nIXEO4dntJ+0uJWwFFIwislJEjohIlYg83M/yRBF5zlpeISLF1vwcEXlTRDpF5MfBLV0NR4erh3V7\nT3Hb1WMZlRRvdzkqhPJHJXHjVaNZu7OOLrfH7nLC0oDBKCJO4HFgFTANuEdEpvmt9gDQaowpAx4D\nvmfNdwH/BHw1aBWroFi39xQXezzcUz7e7lKUDe4pH0/L+W5eO6hPEexPIEeMC4AqY0y1MaYbeBZY\n7bfOauBX1vRaYIWIiDHmvDFmI76AVGHkmW0nmFowiqt17GJMurYsl3GZyTyz7YTdpYSlQIJxHNC3\nMaLOmtfvOsYYN9AG6EW3YWpfXRv769v5q/LxOnYxRjkcwj0LithU1axPEexHWHS+iMiDIrJDRHY0\nNjbaXU7Ue2b7CZLiHaye7f//m4olH59XRJxDeGrLcbtLCTuBBGM9UNTn50JrXr/riEgckAE0B1qE\nMeYJY8w8Y8y8vLy8QF+mhqDd1cMfd9dz68yxZCRrp0ssy09P4tZZY3h+x0k6XD12lxNWAgnG7cAk\nESkRkQRgDbDOb511wL3W9J3AG8YYvYVHGHp++0nOd3u4b3Gx3aWoMPCpJSV0drl1wLefAYPRajN8\nCHgVOAQ8b4w5ICLfEpHbrNWeBHJEpAr4MnBpSI+I1AL/BtwnInX99GirEPF4Db/cXMv84ixmaqeL\nAmYXZTJnfCa/3FyLV29HdklcICsZY9YD6/3mPdJn2gXcdZnXFg+jPhVEGw6dpa71Iv9wy1V2l6LC\nyP1LSvibZ3bzxuEGbpymzxOHMOl8UaHxi001jMtM5ib95Vd9rJxRQEF6Ej/fVGN3KWFDgzFGHDrd\nztbqFv560QR92JV6n3ing08tKWbzsWb2nDxndzlhQf9CYsQT71STHO9kzfyigVdWMecTCyeQnhTH\nT96ssruUsKDBGANONF9g3d5TfKJ8PJkpCXaXo8JQWmIc9y0p4S8Hz3L0bIfd5dhOgzEG/Mfbx3CK\n8JllpXaXosLYpxYXk5Lg5D/eOmZ3KbbTYIxyZ9pcvLCzjrvmFTI6PcnuclQYy0pN4BPl41m39xQn\nmi/YXY6gcK3DAAAQsUlEQVStNBij3H+9W43HGD63fKLdpagI8OlrS3E6hB+9UWl3KbbSYIxijR1d\n/LbiBKuvHktRtj66QA1sdHoS9y6awO931VEZw22NGoxR7MdvVNLj8fK3KybZXYqKIJ+/royUhDh+\n8JcjdpdiGw3GKHWi+QK/3XaCu+cXUZybanc5KoJkpybwmWtLefXA2Zgd16jBGKX+7bUjOB2iR4tq\nSB64toSc1AS+/8phYvF+MBqMUejgqXb+uPcUn1pSoj3RakjSEuP4mxvK2HysOSYff6DBGGWMMXxn\n/SFGJcbxuWXaE62G7hMLJzB5dBr/8vJBXD2x9dAsDcYo88r+M2ysauIrN00hI0VvRKuGLt7p4NGP\nTOdky0WeeKfa7nJCSoMxilzs9vAvLx1kasEoPqFP/1NBsLgsl1tmFvCTt6qoP3fR7nJCRoMxivzk\nrSpOtbn41uoZegcdFTS99+/85z/uj5mOGP3riRLHGjv56dvVfHT2WBaUZNtdjooihVkpfOVDU9hw\nqIF1e0/ZXU5IaDBGAY/X8LXf7SUl0ck/3Kp351bBd//SEuaMz+Sf1x2gsaPL7nJGnAZjFHhyYzW7\nTpzjm7dNJ3+UDs9Rwed0CP/nzllc6PLwSAycUmswRriqhk5+8Jej3DRtNLddPdbuclQUK8sfxZc+\nNIk/7z8T9U8V1GCMYF1uD3//3B5SEpz86+0zEBG7S1JR7rPLJrKoNIdH/nggqm8yocEYwb798iH2\n1bfx/Ttm6Sm0CgmnQ/j3NbNJSXDy0G93R+3Abw3GCPXSe6d4astxPr20hJumF9hdjooh+elJPHb3\nbI6c7eB//SE62xs1GCNQVUMHD7+wj7njM/n6qql2l6Ni0LLJefzdikms3VkXlVfFxNldgBqcps4u\nPvXL7STFO/h/fzWXeB3IrWzydysmUdXYyXdfOUxJbmpUnbnoX1UEcfV4+MxTO2ho7+Jn985nXGay\n3SWpGOZwCP/3rquZNS6Dv3t2D7tPtNpdUtBoMEYIt8fLl57dw56T5/j3NbOZXZRpd0lKkRTv5L/+\nxzzyRiVy78+3ceBUm90lBYUGYwRwe7z8/fN7eeXAGf7p1mmsnDHG7pKUuiQ/PYnffLqctMQ4/vrJ\nbVExjEeDMcx5vIav/m4vf9p7iodXTeX+pSV2l6TUBxRlp/CbzyzE6RDWPLGV/fWRfeSowRjGXD0e\nvvibXfxhzym+dvMUfQSqCmsluak8++BCkuKdrHliK5urmuwuacg0GMNUy/luPvGzCl49eIZ/+vA0\nvnh9md0lKTWgiXlpvPD5xYzLTOa+X2zndztO2l3SkGgwhqGDp9r52E82sa++jZ/81Vwe0NNnFUEK\nMpJ4/rOLmFecxdfWvsf/+sM+ut1eu8saFA3GMGKM4bntJ7j9J5u40O3hmc+Us2qmdrSoyJOREs9T\n9y/gs8tK+fXWE9z10y1UN3baXVbANBjDRGNHFw/9djdff2Ef84uzWf9313LNBL3hrIpccU4H37jl\nKv7jE3OpbTrPLT96l19sqsHrDf9LCPXKF5t5vYbf767nX146yMVuz6VOFqdD75SjosOqmWOYOyGL\nh194j2/+6SB/3HOKR2+bHtZjcTUYbVRR3cx31h9ib10bc8dn8v07Z1GWP8ruspQKutHpSfz8vvn8\nflc9333lMB99fBMfmzuOL62YzPicFLvL+wANxhAzxrCtpoWfvHWMt482UpCexA/uuprb54zTo0QV\n1USEO64p5Kbpo/nxm1X8YlMtf9xzio/NGcfnrpvIxLw0u0u8RIMxRFw9HtbvO83TW4+z+8Q5clIT\n+NrNU7h/SQnJCU67y1MqZEYlxfONVVdx/5IS/uOtYzyz7QS/21nHtZNy+eTCCdwwNd/2m6NoMI4g\nt8dLRU0LL+87zUt7T9HuclOck8K/rJ7OXfOKSIrXQFSxa3R6Eo/eNp0vXl/Gc9tP8OutJ/js0zvJ\nSonn1lljuHXmWOYVZ9kSkhqMQXamzcWmqiberWzkncomWs53kxzv5EPTRrNmQRGLSnP0EQRK9ZE3\nKpGHbpjEZ5dP5O0jjfxhTz1rd9bx660nSE+KY9nkPJaU5bKwNIfinJSQ/P0EFIwishL4d8AJ/MwY\n812/5YnAU8A1QDNwtzGm1lr2DeABwAP8rTHm1aBVbyNjDM3nuzl8uoODp9t4r66N3SfOUX/uIgA5\nqQksm5TLyhkFLJ+cr6fLSg0g3ungxmmjuXHaaDq73GysbOT1Qw28eaSRl947DUBuWgKzCjOZVZjB\nqhljmFIwMp2VAwajiDiBx4EPAXXAdhFZZ4w52Ge1B4BWY0yZiKwBvgfcLSLTgDXAdGAssEFEJhtj\nwv5BEa4eD83nu2nq6KKho4uGDhf1rRepa73I8ZYL1DR20u5yX1p/TEYSc8dncf/SEspLspk2Jh2H\ndqYoNSRpiXGsnDGGlTPGYIzhWON5tlY3s+tEK/vq2njzSANFWSn2BSOwAKgyxlQDiMizwGqgbzCu\nBh61ptcCPxbf8e5q4FljTBdQIyJV1vttCU75sKO2hebz3Xi9Bo8xeLy+L7fH4PYa3F4v3W4vPR5D\nt9tLt8dDV4+Xiz0e31e3h84uN+e73HR2uWm/6KbtYg8X+3nIT5xDGJOZRFFWCrfNHktJbhpTC0Zx\n1Zh0slMTgrVLSqk+RISy/DTK8tP45MIJAHR2uRnJ445AgnEc0PdK8Dqg/HLrGGPcItIG5Fjzt/q9\ndtyQq+3HD/5yhK3VLQGvH+cQEuIcJMc7SYp3kpzgJDUxjrREJ/mjkkhPjiM9KZ6s1ARyUhPITk1g\ndHoSo9OTyBuVqENqlAoDaYkj2z0SFp0vIvIg8CDA+PHjB/Xa79w+kwvdHpwOufQVZ32PdzoufU+M\ncxDnEOL0GSlKqQEEEoz1QFGfnwutef2tUycicUAGvk6YQF6LMeYJ4AmAefPmDepCytIwGhSqlIoO\ngRw+bQcmiUiJiCTg60xZ57fOOuBea/pO4A3je9jsOmCNiCSKSAkwCdgWnNKVUmpkDHjEaLUZPgS8\nim+4zs+NMQdE5FvADmPMOuBJ4Gmrc6UFX3hirfc8vo4aN/DFSOiRVkrFtoDaGI0x64H1fvMe6TPt\nAu66zGu/DXx7GDUqpVRIaU+EUkr50WBUSik/GoxKKeVHg1EppfxoMCqllB8NRqWU8qPBqJRSfsR3\ngUr4EJFG4PggX5YLNI1AOZGw/Vje91jffizv+1C2P8EYkxfIimEXjEMhIjuMMfNicfuxvO+xvv1Y\n3veR3r6eSiullB8NRqWU8hMtwfhEDG8/lvc91rcfy/s+otuPijZGpZQKpmg5YlRKqaDRYFRKKT8R\nE4wicpeIHBARr4hctoteRFaKyBERqRKRh/vMLxGRCmv+c9bdyAez/WwReU1EKq3vWf2sc72I7Onz\n5RKRj1rLfikiNX2WzQ7mtq31PH3ef12f+aHY99kissX6jN4Tkbv7LBv0vl/uc+yzPNHalypr34r7\nLPuGNf+IiNw8mH0dxPa/LCIHrX19XUQm9FnW7+cQ5O3fJyKNfbbz6T7L7rU+q0oRudf/tUHY9mN9\ntntURM71WRaMff+5iDSIyP7LLBcR+ZFV33siMrfPsmHt+yXGmIj4Aq4CpgBvAfMus44TOAaUAgnA\nXmCatex5YI01/Z/A5we5/e8DD1vTDwPfG2D9bHx3M0+xfv4lcOcQ9z2gbQOdl5k/4vsOTAYmWdNj\ngdNA5lD2/UqfY591vgD8pzW9BnjOmp5mrZ8IlFjv4xzk/gay/ev7fLaf793+lT6HIG//PuDHl/m9\nq7a+Z1nTWcHctt/6f4Pvrv5B2XfrPZYBc4H9l1l+C/BnQICFQEUw9r3vV8QcMRpjDhljjgyw2qVn\nYBtjuoFngdUiIsAN+J55DfAr4KODLGG19bpAX38n8GdjzIVBbicY274kVPtujDlqjKm0pk8BDUBA\nVxn0o9/P8Qo1rQVWWPt66VnmxpgaoPdZ5kHdvjHmzT6f7VZ8D3oLlkD2/3JuBl4zxrQYY1qB14CV\nI7jte4BnBvH+AzLGvIPvoOJyVgNPGZ+tQKaIjGH4+35JxARjgPp7BvY4fM+4PmeMcfvNH4zRxpjT\n1vQZYPQA66/hg78w37YO/R8TkcQR2HaSiOwQka29p/DYsO8isgDf0caxPrMHs++X+xz7Xcfat77P\nMh/otQMZ7Hs8gO8Ipld/n8NIbP8O6990rYj0Po1zuPsf8Out5oMS4I0+s4e778OpMRifPRAmz5Xu\nJSIbgIJ+Fv2jMeaPdm6/7w/GGCMilx3nZP3vNRPfA8R6fQNfqCTgG3/1deBbQd72BGNMvYiUAm+I\nyD58gTGgIO/708C9xhivNfuK+x7JROSTwDxgeZ/ZH/gcjDHH+n+HIfsT8IwxpktEPovv6PmGIG9j\nIGuAteb9D7gLxb6PuLAKRmPMjcN8i8s9x7oZ3+F2nHV0cbnnW192+yJyVkTGGGNOW3/8DVeo4+PA\ni8aYnj7v3XvE1SUivwC+GuxtG2Pqre/VIvIWMAd4gRDtu4ikAy/j+49sa6D73o8Rf5Z5ELaPiNyI\n7z+O5caYrt75l/kcBhMOA27fGNPc58ef4WsH7n3tdX6vfSuY2+5jDfBFv7qGu+/DqXG4+/7fhttQ\nGuovrtz5EoevwbWE/244nm4t+x3v74D4wiC3+394fwfE96+w7lbger95Y6zvAvwQ+G4wt42vsTnR\nms4FKvnvjqcR33fr3/t14Ev9LBvUvl/pc+yzzhd5f+fL89b0dN7f+VLN4DtfAtl+7x/8pEA/hyBv\nf0yf6duBrdZ0NlBj1ZFlTWcHc9vWelOBWqyLRIK1733eq5jLd77cyvs7X7YFY9/ft42hvMiOL+vD\nrwO6gLPAq9b8scD6PuvdAhy1fmn/sc/8UmAbvsb43/V+gIPYfo71h18JbOj9B8d3GvUzvw+0HnD4\nvf4NYB+wH/g1kBbMbQOLrfffa31/IJT7DnwS6AH29PmaPdR97+9zxHf6fZs1nWTtS5W1b6V9XvuP\n1uuOAKuG+Ps20PY3WL+Hvfu6bqDPIcjb/9/AAWs7bwJT+7z2fuvfpQr4VLC3bf38KH7/wQVx35/B\nN6qhB9/f/APA54DPWcsFeNyqbx99DpSGu++9X3pJoFJK+Ym2XmmllBo2DUallPKjwaiUUn40GJVS\nyo8Go1JK+dFgVGFDRDJF5AvW9HUi8tIgX3+fiIwdmepULNFgVOEkE99dc4bqPnzjWpUaFh3HqMKG\niPTeyeUIvsG95/E9N3gGsBP4pDHGiMg1wL8Badby+4Al+G5vVg9cBBYBXwM+AiQDm4HPGv2FVwHQ\nYFRhQ3w3m33JGDNDRK4D/ojvEr9TwCZ8QVcBvA2sNsY0iu+GuDcbY+63rs39qjFmh/V+2caYFmv6\naXyXDf4ptHulIlFY3URCKT/bjDF1ACKyB9/llufwHUG+5rv9Ik58l4/153oR+Z9ACr7raA/guyuN\nUlekwajCWVefaQ++31cBDhhjFl3phSKSBPwE33W0J0XkUXzXVys1IO18UeGkAxg1wDpHgDwRWQQg\nIvEiMr2f1/eGYJOIpOG7o7pSAdEjRhU2jDHNIrLJegjSRXx3r/Ffp1tE7gR+JCIZ+H6Hf4jvNPmX\nwH+KSG/ny3/hu6PPGWB7aPZCRQPtfFFKKT96Kq2UUn40GJVSyo8Go1JK+dFgVEopPxqMSinlR4NR\nKaX8aDAqpZSf/w8S9tHy9GxOHQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f1aba526450>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "likelihood_bkg_only_example = likelihood_bkg_only.subs(example)\n",
    "likelihood_bkg_only_example = lambdify(theta, likelihood_bkg_only_example, modules=['numpy'])\n",
    "X = np.linspace(-1, 1, 100)\n",
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "ax.plot(X, (likelihood_bkg_only_example(X)))\n",
    "ax.axvline(simplify(theta_hat_hat).subs(example), color='k')\n",
    "ax.set_xlabel(theta)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAACKMAAAA/BAMAAADEJYTOAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAEM3dMna7RFTvIomZ\nZquuw15VAAAACXBIWXMAAA7EAAAOxAGVKw4bAAAceklEQVR4Ae1dC5BkZXU+Pf2cfky3EEytIDMO\nJIoSZlCRCMq2kqQkWsxEMmhgYdsQVkCy04WWgbLKacXgoyx3EIPoInstihQYYBrNanTQvT7QAMvu\nRAFjAtkGCh8gO7OAC/vsfOd/3Hu772Pu7Pbsdvf2qZr7v84595zz3//0///3P3eIeqAsMNGzRCdY\n4KqhTpCyW2VMmt2q2TLolRpdBqY9li23QMZsOcsew/AWuDQ86hGPeUfxiDdBRxigr9IRYnarkCvK\n3apZ6/U6rfUsexyXwQJ5YxmY9liGtUDyq2Exl4aX+ekPbYKSne3gXGboIIW/98qCxaFk5XqZFlng\nsYcnTKKND9LUxBOUmai0iG2PzVIt8KOlEoTD/x49aSGmdljZTs7MHOSELmMMDGn9u8QkWp22SH+5\nj75K0fLA3HmFVfSzWKkthDoShRisLIvWZ9CYoRkf2x0uZY3W5wDTqBHbrUm7xCRanXZI08eM0g7a\nQpnapfRzWvMJe0bYDtIdSTJkppdF29PJ/lH/m65wKakXtaEyBzZdiZbTL2kW3WESrc1hTSM1dfuB\n8eQ03USJ6pfhWM6tFw+rVEfyzSN/WCbt1+mfidgnusKl9A1pQ92rM0tNreVOl5hkqfovD/5Vim1/\nMWsWd9NgcQfdlDIGxpfnbj2ui1tgkx76i6PGJ/Eze+2uDxD9X3379jX4sY5P7ty+feUQ0aP10+io\n+mVEF9fX/vpkZmm9H8lEusKljM0pA2VqKiOS8+tlSq28sUC5+f3bf/1B1J1fX7v9qeeJcpMvjEdG\ndlYoPr9z+9/zbnV0XFF6msTfuhYn27qKUdcmlln9nzRt1tyQtMIMRa8unkjP0qm5UrSYLXatbdpe\nsRk9VEJImtgJpNcw4o1EWV4J5IcwdkrIzM9R7gmkhMGUuJ7IPm/0P93hUuZZOYbPssO0YbJE1MdV\n/chsLSIDC9Cj+BvYB9/LLWyky5GxpjfeJvG3rsVJWxdMuxwss/o/adqsH5Km+DjFv093X1ihc2cp\n9/j7utw+7ayePaF3SpmCV3DD1ZMGURX1MYyW5F5kpsaJIiYyY0N0LhKK7SFKYTH1GBcYIuNd4VIi\n+6U6FPmyysgkeR32XDOcH5uDOUzYhZeSx+Ev9jzFakhpDCZbN045kwsAH5P4W9fipK0r+XTx1TKr\n/5OmzUorKl1siE5ULW69hnBIf/wWz8XKPRgySe5AJsrx0NlU5BGCTN+L6RoSimNtFN9DyWqyyEX4\nl82b92KwdTrkoJeAhClTdc38MTzr1VxYbRBtgykyPHs7hmtGindzQqvL+KvSx+gjouhnkgDrak7K\nupJNN18ts/o/adqslJruZkt0oG4RTCvcEPd0KWZ0gTIFYCfQiytORWaEZy0GMqn9cvgkFjDjf4ne\n+ZnrGE9A2pOXbu2QlPUSYG07y2KWVhLdw/nL4VxvgNJ9sBGluWbdMzVOaBJ/8+XY6z6DXSgJniYJ\nsK7mpKyr2XRvapnV/0lTZoUNzuxeO3SmZid5ie3jUnL7KMvo+RrFzmKfsZPic5J+gxwv+SGiGZM2\n1OsYYffX6/gNf3R/RaJ08pX1EvA2laokS9vmyOQCVkZn825SvkQXcRm7KzvZRGJzJb2X+uv1BayI\nYBnT2ySmv3U1J2VdwbWrL5ZZ/Z803leCWQHrjK62hUO5AflAOWraMstrFxd4uhSseV5Q+wZPT1xS\nAVHy+e0bDEm8aVyk2DdI6V+Ns3927RWqVeJ08HWqKoXPYcnnhKspPyq3l/ZPnP/n3DJ1yu8UToaX\nhnJ76VpFHnvb1W/5hqh1XYKsqzjx9oFlXReDbqqwzDrm+6TxvpI0a9QMo3r2ZTr7A5TZM/tbbKPP\nfnuBsjvpfnTLptkthqKP19sQFmzlkt/Vkoq6thWXd1hd4OlSsObZWv4nxl1tUOQWpNg36FO089Mi\ns/rZe+8zZFXkJso57CErO/W6qSwlj5YaNbiHsi8luY0X/StqyMBFDyEBDOwSSXz/7OYLRA4YFdqm\nss1JkHUVJ3JYV5Jn23AMHIhIzcawzOr/pNlmzYVaWSe/RXH49yG6hujDRDWi7eLFwtmUGlK3/1iz\nGO5ybPYKWZm+6G/drSFqrMWvA9erztFsZ6OGnSdszi0KCVaVIUsDMrO0a+56N75XXSPWYKmxLEqe\nLgVrnpmS3jeg1RWxb5ClGFNk735ZEE4WRCKqhigZqrMtijbOjCi9xpRr0aKalN5pbS/F2QbYXjKp\nyO0PzFc4STiMgJCqKdtC3GpBkHUVJ3JYV9L9q0XeXRnLrNih8nnSbLOGO64Z75uDS4lX6T6cniJj\nHA+zGJOrbJfyAUpNztGH11Ny8vdb4DpSI7/fMmLQL14wEjvL9NH9m9caXmNUxMhet4+Sb60GUMoO\nco7HkuozZ11yWFXKBB7QCY0uJYS4fTQoH7fEFe/WjJYi7kdO1VREJZV11nmLG522qaycp0vJYHjs\nNhmH36jiHQbvG0Tk242fxPjMCmFlZEF+juKjVqnDM1qv+5v0gBEm7+Q6LPr5VRdvL1FBvATLFsdK\nKPPKyIJ9ROdZhcZMgHU1J6d1JfFvGnl0Tckyq/+T5jBr4zD0MUI8Of1qzBy//2iZ6NL0eAE/g9OM\n+hwdVZQkqXFs/I1SHMVfEX0XlX+HJ75AeIPBQRi5acqOeo1RESObu4SI5y2+lPIejvFoHaZ21DkO\nLzF+bMuoJFPXBpcSRlysDSVpZOWQ5rQEcSNztkvR4jrrfMQdwIzdBZ4uBe9Kk3urwOXDAniHIY+l\niPcdMPlKg1vEfhkygJki9Zki1/mXtHSYRLc06oINEJp6kOv4WApedcljKWLe8h2KCtOOmdwuII0X\nzCfrQlMaYF3NyWldyW+oiUmXFC2z+j9p5DArVkeLQ5z+DC/x+wuIV6ZbvlUEQUYcDLh54gJFnED6\nvZfEC4gt2PQtEP1anjnaSr9E00CN4gteY1TEyGZmymQCK4ASrWo81pAlHTvaOEZfx00MYqGSGpUF\ndW1wKaHEjbxekuZm9RNMSxA3TsKl1JiHFlfVSbbkKW6G/UMzeLqUe4C1sowL7xuIWcmmIjI11LwP\nZzJMpGLqj5QBU65z0C1dASnlKsUj7tCIfUdUPJn8WPOsRBxLORco8TLlBBWvEBXg1F/sel1oSv2t\na3FyWldQZ7k3QoPjIza5jT8MTRaMyA+2L5R8W4IbLLP6P2liOaS4zMwFsxOtt9M5VaKjif6R0jtS\nBMtVv4mGyIJF+zhy1TXCpZxGhDUqrad35JDMvB/YFJ2jBGbl7jEqYmQzaDOB5U+JRuzmiDFa4ryO\nHW0coxhaEoRhg1xKKHGPL0huA3SroRgvQdx3z64vgqrElFpcVcdVAE9xczxbbwZPlzIKLLG52DdN\ndA3vKeFMBl1ThEc1Mb/fgUIfXxQkKunTdb7T07h8eSOcqVOXKApyD/okTFAQ6kM4ukPJL6H69/ir\no4JG+KLgMrq1rPNN6SjK3ta1ODmtK6iPauIRXHR8xOZhemUwbuhW/jH3Az1X9mv3rbfM6v+kOc3a\nX/PlZDd8jh0CfYF9emqI4gZOWK1AXyCvoQQHMz5lHIPkRYq/BslLcjQN3MQog0U6r0A8Ro/bUN/t\nGKNiRnNMcjcvdwMomcoxRjl2lI9ZGFyXe+gvV+HnaFU1+YeLC7Tx/bidh0uJrf5PFkRBKYS4kQeO\nleImhKdUlOizcOISXc8kuBNZ4nLdIuImebnaDB4uJXL5C3NyOjYwggi5y6DQk4gevATHA95Z/xIl\nR16oUXZkJ8yqIHLfhTD/IYa7PO+XMzyrw1dm90lcPnnlgFdueDtKNfyxKdZUiI5n22DdGhvZVaSr\n6mfh165+KdoVHDdxj842pv7WtTk1WFeQO5+wRn5eJcdHbK6w3tR5IS6hLrIpoI/1XHkJ/ASqZVb/\nJ63BrInRxe+QW0+xcUoN8yvja98wewplzqQ7vkb0SfmzdztmEBU4mGLU/DSvXifWgmVq749LSAgH\nhgBTs7MV4jEaOe2u3YWGMYqpyafobalyICVTOcYox44efcG2r4MzxujUXL9Bp9Nz8VPjZrb8V6j0\ncCkshYSw4kbrL0hxkw9+R9OKmVRIcY/da4CshD9bXNQtIm4aPqEZsr97/rnmuhaVrRhT/9Dd5oDb\nO3/6D86bRy63nqDYk1VnC84+lZzlV96gS8/oTEPq4KTq05vXFJE99w0NeOhfsbiRsYFNTYexGGFD\nuAzqDg9XQbuOj9jcX+a1mROkkV3MrMBqK9ZcRwBL4rh6o+Bg9R4rr+fKqsLuDgujNZnMwsHyScIz\n8rwtQ7kdeKQymMacimRU7LPTbScwf14mEY/RrEn4CXWM0YyJFRNtw4QggFJQkT1GOXZ0M/UXsUux\n1/j6eZQvZGv0g8R4bihx2SvAqdGl5B2v5hcotLjgI8XljIYliCtJSkhscVFYTNxdku5QXSdLWB8V\ncLd+ZLYWkWGf9ij+fAJuIzXSp0KABOiriYQvK+ysqOszRKIvp+jMUTrTmDo4yYbj5Z6dOHTjRO1T\na8ZB01l7uPOZOZbAZdD8EH5YS2iZ1+HhXHM5bG6tUo7eaaLdAdrILmZW77iYCep3uV3KWzRb12dh\nrO7QGCL1DlFtQFmkYAVgeeFdtwZQIMeQdGfxTFXpZ6DOEv2LieGMl9gvo9KMGK/GbLw6VUbblfgT\nkDfw/TiGFLPGBsBjyJvUj6ULJXwpLaqB4eEvDv+JiB3dDUcCIsDTOC7TX6EP9Zeztdj8i/Tk8EnD\nw5hVjcpmxzXBFG0ubt0h7zJm+Z0XwIoxDQjd5YORHM5MggankPunBa265CtWKVGzsiLDnesAPBYS\nMvLmuqhTBydZ9QTRVzi3VmOotF8+QzTmzacJewlFZZUlUDhQsdfoZVBXeLgVtIsJuoTchjN0VqbK\nyAG9M6ZizS1mgrCqXYqlSPxGzdn1WRirOwTGP0s8O0TVYqEZhE1j/KwsCcS4dFIkd6RrKGfw8riK\nfcEKpfAjwo/HPZhFUHQIbZb0g00jXcTIVimzHkj+lE6qEh5thNP+14J+x4tpz+ZY3oi8mC/0Vz4H\nB9M8S0GFgl8VkGlzcTewjGr4itzyXNQDY8WY8hFcn9BdK+BW0PSP0gBevtowJQXmimaX0jhQ7JN2\nkWmb3JFzcJK1WCnNF5Ftdin5Bdm+riLTll0PeBixBBfxxW3QTUX8so6jSYeH66DdjIlKAbOpr9hG\n5BplZDczK7C6mZngkzNcLuWuPxUtuDR/FsbuDoGhXApZ+3YHbIvIXn3LsKkYlw3IJw8YKF+NXdgK\nXt0V6I4qEoNSJkVKxF8Jtb8cHC30FYFrAcfIJrHNwj7Hn9JJVWJaxI6eSidyDvkFOv1TfcWBarTw\nKL2ZfoKqBNenpvnqhOTZouPaW9yVhhD5gDvUqXBAXvG3YkwDQnetgFtBg/ccjS7lW/ZdmlxKetRu\nQi5jF1/b0KALDk6yaqQQ7FK2FTVpi9KDsrpQz21QV3i4Dtq15nDYhDm60qCBMrKbWZ8OrNax5pqZ\nID9Gn86UE0quq2qX4vosjKM7GLF1LkUcvWSWYUGNSyf61MdRyt5QpqxBmZWbt3wQjw+SkWJk/no6\naleV3it2Kpkk9vQDTkoRI3s/dnlRG0DppCoxPcJpz11lcA6zjhpdOJd++j8o9myVHn+kjCp2Kekf\n7x9HIkBNrF4dLXCxvcU91C5l26Khu/khnAUwYTg94KI7kj/4xgRPBzdeXKZTJvCmIzPxMIpwKZGN\nF5Yp+exDEyVs2puIfHjmtovRD5mJ31Df06sqILmo6Dzl58OJNKzkLsMs5dgJ7JF//N5HzucGPUvh\nmVVLQWu4VKYsZBxq8vK/2aCu8HAVtCsm6CI8HLvY8QoTOyGKxZ2bWd4KrFax5joCGO9E6+bnN299\nSPLQiuQMuBRxE/mlHO/uYJoWupSdTj0c+WMnGge/blLjUhc51btlzrply5dCcBazFIEntVATq3Hp\nUtpYXMg8UhSC62dCFJbhovhbMaYx39BdO+BWyzRVzSYupUEDpzoiuyPPF5B7rdgug0s526BLETgT\nncMAS8xhn2zNeGQ34UjMoDH4GuT6ygmT8H8YNPhw0s1J7OOwS4mdQYlxTNXfHufy4AJfERNoiKR1\nF63hEjneWgVBVBC5DOoOD1dBu+IjNjI8/PWEzZMm4GBrFzNHYLWKNVfMdEz1jOKjFTmG4FJUDDp/\nFsa7O/jWLXQpX2xSRRWTZ9C/ebaocelsi407S8ucL4Xgb7kUqYWaWMVmv/vvTNy+4rJ0k2W+WjMC\nUViGi3rmrBhT/9BdO+BWP6dn0buwQzZTxHFGOhMHqPJVLD7fASHhUs7EArb4FA2MohitUKzyRsK5\n9/4qRavrCrSHZiiGp/5aSyMfTro9WuPcWlBT7EUwfUo0DE6LhLaCVUtBa7g0psmfszz/K4hcBnWH\nh6ugXUwrSIWHH7fxe6474gwNNiP9A6vnpQ0UsxUVcTAvN/lbyUgrUoVL0THomNr7dAfTtNClYAvM\nC+6oRHhXwgV6XLoaDlWFEeJGaY0jteCJ1SfreNj7W/0I6vv4p4Z/k9Viics1i7iUo97EcCI534y7\nX8MF11BmePiELw4P1/CrtXjorgq4tWmyJhUGy3iEs/jZ258dosFqZBe/84NLiWFTLF+9X7qU/jJF\ncOQxt8BjP1/7EUd8Dd5cBuIYRU5gPd5Y9OSklSQcFWBYS2MVoj3Z6SaXMi9aNavgN5PBJlkAJ1vD\nJVkXpKl9ID9VCOMyKHY/msLDEzCbAhjPsUtqq401jwkUFzM7sFrHmitmP8J63nq4bUVSRbgUx018\nuoO+OTz81uHhU1gquT1rs8Ch4yUAc6C3iqvr8rtHLrJEbGw8DOOyUYAllKQWmFilapiXRjY9sATa\nw4TKLsXq0ObOxOPbIlA/Y1X4sDuZZb6Gp2kPMlgIFyLjyLgDbhUNv/BfhxOqfJQl/Yc+k6Yq9NF6\nCbWJWgqH0PKlmSJ+NeHBoUpqgQaG+BP+M+MnU/JlSk7ySnsMfwq8OanGuCkya2kdWO1Mv4gfXAa9\n8JkXpRZelIZL5ojf5RRMCcC10aDY/WgKD8+PWvzzvuHhbGQ3M7t3dKy5Ygafdp7FFRmlyLEEl+K8\niXd3MGULZyk+LmW/U0BnviPGpRZYaMETq76COM+l69s4XWSW0irJ5TOXrODX7UHmOTaHw23wB+KL\n8hwnRu6AW0mDF3nf4CiufZFBvP3fjQd2a+F2OmY3SNQsxbxto1hhYuHDP5H5uVfAV20y8FnLHelK\nZLUhProDdAHenFTj5yltICtnKXvTF15YFA16e7ZNFj74tERZnhckt0H5WAofp+CI/JUGiz9m8lXA\nTNEnPFwY2c3M6p0mZq6YauVSvrB9+66nnTfx7g6WpYUu5U1SueYr/2B1Pigt+gt9dHeHaDNSFoKq\nZ2LZhJb82XcsFrprB9xKmvcQXYBdlNSODLxJfzVfjnyNTiJagKjYS7kZPqr8aSk2b8/2V3FItIof\n0fX0Eq2Y65umPBrt7VkfTpJBzKS4gexa6jcp+VJKmgbzIL4ZYJMhktZdDtTqgyX6vJDCbdBNReze\n1dD4Ppz7MZE6g3ZxjOScAtc1gzCym5kVWN3EzBVTbStyY0MMund38N1b6FJ8tmfXYx7crGfnlaUW\nmFglv1/tEOlHikJQ+5mw5S7Z2YPOSf5R8MmJ8Qmn4BO6izMSCgRN+szZLaP4f8KZoYezNbznGZhL\nmHSFfIWaKNE5RXoDXfWRv2aSjIn1Tpm+mStuoYSJr+tcSQmDJtAiJkZIiXw4iTa6dfbepzi3hmJf\nwqnn5Gv/whAN2qWsLopi6y5eVg/DPbubXi/w3AZdifpripgCmnCFwpQjthPxCw8XRiY3Mx1Y7WLW\nHFNtK7KTnDfx7g4WvYUuBctcL3gvovy86jurrgO1WFkUJrafCcviBxycbnFwZAR/K8Y0IHTXEXAr\naPqwvzOaHKX0cxWaXYUN2PseIrrtoveDd2x+v5F+8AdYBNTrmK0QsPhjm3f/N2U2Yvv2uFVzlJ54\npIKWE/EnwZuTahypIywMYYP1Z+jO2d/gBXt9l8lN+kD+tjKXWggeVg/FPb0nWWNEl0Hd4eENQbt+\n4eHCyC5mVmC1FWtuMWuOqbYU+Xb9OedNvLuDRVcuxQ5RtVhw65Jgrzd2/L7HvRs6qrYDtZDrbb2/\n5rR26OB0FcDrpHXlD+SBCU1zGb3q2nHc8gzXbWVFGr7mACBRidwu4mL4GBjDuvAupSl+WtK7roEa\nBrGY7DNczA5nRaAinoLpWYrVGMgiyBbe/4vK4tvLHGoLbJA3VB3q7Lqm4HR/wVQArz8CWiJGYLNn\nY1iaVEnF3j7myQYrIvY3S4frQPIokyUQ5MEQPmwwUmuKn5YMmq9BGgaymNrezOrwloMU8ZbsVc3V\nQSwCbbH0sMHmW/fKLbVAXXKTHersOldwuu9tn1ABvL4Iy9wQeTNGfQU3wc6JJxzlWbto5dFFossY\nC/8ERkDelOniV5xO7Z9eHC0II5BFVDm5IAZd1BZoi5TqnS7St6NVifDmgQXOrnMFp1tYzZkbVABv\nc/0hKx+78cKf8M3SJe9bXuFdvVhtZHbi3iIj6Q9FRoe4FAb6m+OnwxA14gSyiJcakbu8FGiLzO4u\n177D1FP/LUNJ7ey65uB0f8XwSmG+6N98CFvu8rxXzvCsDl+Z2iNxswthaVRob1h0L7wWsPBi25F1\ngbYYCN0rHal7xwmtP/4uBXd0nSs4PVC1lYXA5g5vjKmDmOH+sZ1SVu/pHoTuLWBxEHdvL1J/W/SX\n2kvSI12a+L4mC6DrHMHpTY1+RRnA69fa8fWR56UKHFEUGji09yChBSwOUoL2Ife3xcx4+0jZkwQ7\nj837fOg6R3B6SBNFayEROxSNz48x3CKTUNezQmEFIrWARSD/Tmr0t8Vqo5P06H5ZEbjaCGc5g9Mb\nm/xLm/2buqJla1GqsTW8NlkzPK4PZgtY+HDuvOoAW1zSedp0tcT50Ub10HXZIWcEfGOzdyluetd3\nTe26OanKVDm0SqtCY/oitoCFL+9Oa/C3Rbh/s95p+nawvGNmo/DourxvBHwjql1SAbx2RbflZmpS\no6gZVjMR2hsW2RuvBSy8GXdgbYAtrO9hd6BaXSny6kqDWtx1M34R8A2YjoIO4HVUdVkWL9cFhH/l\nI0J7D84KLWBxcAK0EXWALfprbSRnTxScUWs0AnedbwR8I6pd0gG8dk235azTVDeH1EyG9oZE9kZr\nAQtvxh1YG2SLdUYHKtTNIjcGhouucwanh1JdBfCGwu1MpLQ662b/w75F9BChvYvgLNLcAhaL3KFz\nmoNscXLnqHFESCo+22VrKrrOGZxuNx3ZuROU+n3mkW2H9tP+4P8Bavvp1NESDaiw/Y5W4hAIv60o\nb8JfSOxBO1lgRfiXcO0kdvfKEj62tnttEEazfFVhXRsGu4dz6Cxw5aG7Ve9OYSygf33D4B7JOBk9\nm8PHKnvQRhbIldpImJ4osMBwzwrhLLBeoz2tM720HSxwVTsI0ZPBtkCytzVgGyMwt62gmrO9tXug\noQ5tYwSf/+pBO1kg/GnQdpL6cMgS1Zsph+PmvXv2LNApFlinf3s7ReDDJmes9+mww2b73o07yAJf\n6yBZD7OovzjM9+/dvmeBDrDAQK0DhGwTEfm/Y/WgZ4GeBQIt0PvlDTRPY+PrGou9Us8CPQs0WyD9\n1eaaXtnfAn/k39Rr6VmgZwG2QNpY1A7/D9qiXpNRuAnTAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$- \\frac{NBEXP^{2} \\sigma^{4}}{2} + NBEXP + 2 nobs \\log{\\left (nobs \\right )} - nobs - 2 \\log{\\left (\\left(- \\frac{NBEXP^{2} \\sigma^{4}}{2} + \\frac{NBEXP}{2} \\sqrt{NBEXP^{2} \\sigma^{8} - 2 NBEXP \\sigma^{4} + 4 \\sigma^{4} nobs + 1} + \\frac{NBEXP}{2}\\right)^{nobs} e^{- \\frac{NBEXP}{4} \\sqrt{NBEXP^{2} \\sigma^{8} - 2 NBEXP \\sigma^{4} + 4 \\sigma^{4} nobs + 1}} e^{\\frac{1}{4 \\sigma^{4}} \\sqrt{NBEXP^{2} \\sigma^{8} - 2 NBEXP \\sigma^{4} + 4 \\sigma^{4} nobs + 1}} \\right )} + \\frac{1}{2 \\sigma^{4}}$$"
      ],
      "text/plain": [
       "                                                         ⎛                    \n",
       "                                                         ⎜                    \n",
       "                                                         ⎜                    \n",
       "                                                         ⎜⎛                   \n",
       "       2      4                                          ⎜⎜       2      4    \n",
       "  NBEXP ⋅SIGMA                                           ⎜⎜  NBEXP ⋅SIGMA    N\n",
       "- ───────────── + NBEXP + 2⋅nobs⋅log(nobs) - nobs - 2⋅log⎜⎜- ───────────── + ─\n",
       "        2                                                ⎝⎝        2          \n",
       "                                                                              \n",
       "\n",
       "                                                                              \n",
       "                                                                              \n",
       "                                                                     nobs     \n",
       "        ____________________________________________________        ⎞      -NB\n",
       "       ╱      2      8                4          4                  ⎟      ───\n",
       "BEXP⋅╲╱  NBEXP ⋅SIGMA  - 2⋅NBEXP⋅SIGMA  + 4⋅SIGMA ⋅nobs + 1    NBEXP⎟         \n",
       "──────────────────────────────────────────────────────────── + ─────⎟    ⋅ℯ   \n",
       "                             2                                   2  ⎠         \n",
       "                                                                              \n",
       "\n",
       "                                                                 _____________\n",
       "       ____________________________________________________     ╱      2      \n",
       "      ╱      2      8                4          4             ╲╱  NBEXP ⋅SIGMA\n",
       "EXP⋅╲╱  NBEXP ⋅SIGMA  - 2⋅NBEXP⋅SIGMA  + 4⋅SIGMA ⋅nobs + 1    ────────────────\n",
       "────────────────────────────────────────────────────────────                  \n",
       "                            4                                                 \n",
       "                                                            ⋅ℯ                \n",
       "                                                                              \n",
       "                                                                              \n",
       "\n",
       "_______________________________________⎞           \n",
       "8                4          4          ⎟           \n",
       "  - 2⋅NBEXP⋅SIGMA  + 4⋅SIGMA ⋅nobs + 1 ⎟           \n",
       "───────────────────────────────────────⎟           \n",
       "               4                       ⎟           \n",
       "        4⋅SIGMA                        ⎟      1    \n",
       "                                       ⎟ + ────────\n",
       "                                       ⎠          4\n",
       "                                           2⋅SIGMA "
      ]
     },
     "execution_count": 148,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denominator = likelihood.subs({ns: ns_hat, theta: theta_hat})\n",
    "numerator = likelihood_bkg_only.subs({theta: theta_hat_hat})\n",
    "q0 = - 2 * sympy.log(numerator / denominator)\n",
    "significance = sympy.sqrt(q0)\n",
    "sympy.expand(q0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}