Newer
Older
Master_thesis / archive / .ipynb_checkpoints / discovery_zfit_freq-checkpoint.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from lauztat.parameters import POI\n",
    "from lauztat.hypotests import Discovery\n",
    "from lauztat.calculators import FrequentistCalculator\n",
    "from lauztat.config import Config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
      "For more information, please see:\n",
      "  * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
      "  * https://github.com/tensorflow/addons\n",
      "If you depend on functionality not listed there, please file an issue.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "import zfit\n",
    "from zfit import ztf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Signal + background fit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAO0ElEQVR4nO3dcYxlZ1nH8e/Ptkgr4KK7CG47DhjSBAmEOkGwCZJWTG1Ja2L/aBOwRcgkGqAYDS4SJfGvTTQoigkZpYLaFJKCWikolUIaCWzc1hZaFqFihV2qLTQFlI248fGPubsOd2fmnrnn3Dvzznw/yWTvPfe9533mnTu/PXPmnmdSVUiS2vM9212AJGk6BrgkNcoAl6RGGeCS1CgDXJIade48J9u/f38tLi7Oc0pJat4999zztao6ML59rgG+uLjI0aNH5zmlJDUvyb+tt91TKJLUKANckhplgEtSowxwSWqUAS5JjTLAJalREwM8yc1JHk3ywDqP/VqSSrJ/NuVJkjbS5Qj8PcAV4xuTXAS8AvjywDVJkjqYGOBVdTfw+DoP/R7wZsCG4pK0Daa6EjPJ1cCJqro/yaSxy8AywMLCwjTTAXDp4bs48cRJAA7uO59PHrps6n1J0m6w5QBPcgHwVuBnuoyvqhVgBWBpaWnqo/UTT5zk4cNXAbB46I5pdyNJu8Y070L5UeDZwP1JHgYuBO5N8swhC5MkbW7LR+BV9VngGafvj0J8qaq+NmBdkqQJuryN8FbgU8DFSY4nee3sy5IkTTLxCLyqrp/w+OJg1UiSOvNKTElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNarLX6W/OcmjSR5Ys+13knw+yWeS/GWSfbMtU5I0rssR+HuAK8a23Qk8v6peAHwBeMvAdUmSJpgY4FV1N/D42LaPVtWp0d1PAxfOoDZJ0ibOHWAfvwi8f6MHkywDywALCwsDTLfzXXr4Lk48cRKAg/vO55OHLtvmiiTtRr1+iZnkrcAp4JaNxlTVSlUtVdXSgQMH+kzXjBNPnOThw1fx8OGrzgS5JA1t6iPwJDcArwQur6oariRJUhdTBXiSK4BfB36qqr49bEmSpC66vI3wVuBTwMVJjid5LfBO4KnAnUnuS/KuGdcpSRoz8Qi8qq5fZ/O7Z1CLJGkLvBJTkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaNUQ72W21tnXrWrZxlbTbNR/gp1u3jls8dMc2VCNJ8+MpFElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNmhjgSW5O8miSB9Zs+4Ekdyb54ujfp8+2TEnSuC5H4O8Brhjbdgj4WFU9F/jY6L4kaY4mBnhV3Q08Prb5GuC9o9vvBX5u4LokSRNMew78h6rqEYDRv8/YaGCS5SRHkxx97LHHppxOkjRu5r/ErKqVqlqqqqUDBw7MejpJ2jOmDfD/SPIsgNG/jw5XkiSpi2kD/HbghtHtG4C/HqYcSVJXXd5GeCvwKeDiJMeTvBY4DLwiyReBV4zuS5LmaOKfVKuq6zd46PKBa5EkbYFXYkpSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBLUqMmXom5Ex3cdz6Lh+44c3snG6/1k4cu2+aKJO0WTQZ4SyG4ttbTQS5JQ/AUiiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWqUAS5JjTLAJalRBrgkNapXgCf5lSQPJnkgya1JnjxUYZKkzU0d4EkOAm8Elqrq+cA5wHVDFSZJ2lzfUyjnAucnORe4APhq/5IkSV1M3cyqqk4k+V3gy8BJ4KNV9dHxcUmWgWWAhYWFaafr5dLDd3HiiZPA5h0Bu47bbq3UKWm2+pxCeTpwDfBs4IeB70vyqvFxVbVSVUtVtXTgwIHpK+3hxBMnefjwVTx8+Kozwddn3HZrpU5Js9XnFMpPA/9aVY9V1f8AHwR+cpiyJEmT9AnwLwMvSXJBkgCXA8eGKUuSNMnUAV5VR4DbgHuBz472tTJQXZKkCXr9RZ6qehvwtoFqkSRtgVdiSlKjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSo3pdidmig/vOZ/HQHWdub9SKteu4adkSVlJfey7A1wbl6YDuM25ap1vCzmr/knY/T6FIUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJalSvAE+yL8ltST6f5FiSlw5VmCRpc317obwD+NuqujbJk4ALBqhJktTB1AGe5GnAy4AbAarqO8B3hilLkjRJnyPw5wCPAX+a5IXAPcBNVfVfawclWQaWARYWFnpMN7y1LWNP39+KrbaEHW9RO8Q+Je1dfQL8XOAS4A1VdSTJO4BDwG+uHVRVK8AKwNLSUvWYb3B9w3GrLWG7zGebWUld9fkl5nHgeFUdGd2/jdVAlyTNwdQBXlX/DnwlycWjTZcDnxukKknSRH3fhfIG4JbRO1C+BLymf0mSpC56BXhV3QcsDVSLJGkLvBJTkhplgEtSowxwSWqUAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIa1fdS+h2rS+vWafc1rxav2zXvbma7Xu0muzbAh/zGXLuvebZ43a55dzPb9Wo38RSKJDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1qneAJzknyT8l+dAQBUmSuhniCPwm4NgA+5EkbUGvAE9yIXAV8CfDlCNJ6qrvEfjvA28G/neAWiRJWzB1N8IkrwQerap7krx8k3HLwDLAwsLCtNPtGEO2qR1in2ufu9mYPt0Z59mCde1ca211XtvGai/o0072UuDqJFcCTwaeluQvqupVawdV1QqwArC0tFQ95tsRZhEEffbZ5bl926bOswXr2rnW2uq8to3VXjD1KZSqektVXVhVi8B1wF3j4S1Jmh3fBy5JjRrkL/JU1SeATwyxL0lSNx6BS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhplgEtSowa5ElOz6VI4z/0PaRadADfqurh2/+PzrvfcLmvXtf6hPk87J2paBvhAZv1N19I39Sw6AW70+a/d/0adDLe6dl3rH+rztHOipuUpFElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNmjrAk1yU5ONJjiV5MMlNQxYmSdpcn14op4Bfrap7kzwVuCfJnVX1uYFqkyRtYuoj8Kp6pKruHd3+FnAMODhUYZKkzQ3SjTDJIvAi4Mg6jy0DywALCwtDTKeBrG1jupnNWrNOauW60Xw7rW1ql1a047qM2+jz7LJeXWudtI4bfZ3n/TXYKXXsJr0DPMlTgA8Ab6qqb44/XlUrwArA0tJS9Z1Pw9mo/epm1n6jdWnlutF8O61tat9WtF3a3XaZr4utruNGc837a7BT6thNer0LJcl5rIb3LVX1wWFKkiR10eddKAHeDRyrqrcPV5IkqYs+R+CXAq8GLkty3+jjyoHqkiRNMPU58Kr6ByAD1iJJ2gKvxJSkRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYN0k5WO9tmrU6H2u9G+9pozFA1damhJVttS7tRe9+tztWnpet4m9ih9tW1Fe+8WtH2beM7i1oN8D1gVi/wLvvdaMxQNe22PtJbbUu7UXvfrc7Vp6XreJvYofbVpRXvPFvR9m3jO4taPYUiSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVG9AjzJFUn+OclDSQ4NVZQkabKpAzzJOcAfAT8LPA+4PsnzhipMkrS5PkfgLwYeqqovVdV3gPcB1wxTliRpklTVdE9MrgWuqKrXje6/GviJqnr92LhlYHl092Lg68DXpq54d9qPa7KW63E21+Rse2lNfqSqDoxv7NNONutsO+t/g6paAVbOPCk5WlVLPebddVyT7+Z6nM01OZtr0u8UynHgojX3LwS+2q8cSVJXfQL8H4HnJnl2kicB1wG3D1OWJGmSqU+hVNWpJK8H/g44B7i5qh7s8NSVyUP2HNfku7keZ3NNzrbn12TqX2JKkraXV2JKUqMMcElq1MwCfNJl9km+N8n7R48fSbI4q1p2ig5rcmOSx5LcN/p43XbUOS9Jbk7yaJIHNng8Sf5gtF6fSXLJvGucpw7r8fIk31jz+vitedc4T0kuSvLxJMeSPJjkpnXG7KnXyFmqavAPVn+p+S/Ac4AnAfcDzxsb88vAu0a3rwPeP4tadspHxzW5EXjndtc6xzV5GXAJ8MAGj18JfITVaw5eAhzZ7pq3eT1eDnxou+uc43o8C7hkdPupwBfW+Z7ZU6+R8Y9ZHYF3ucz+GuC9o9u3AZcnWe/ioN3C1gNjqupu4PFNhlwD/Fmt+jSwL8mz5lPd/HVYjz2lqh6pqntHt78FHAMOjg3bU6+RcbMK8IPAV9bcP87ZC39mTFWdAr4B/OCM6tkJuqwJwM+PfhS8LclF6zy+l3Rds73kpUnuT/KRJD+23cXMy+gU64uAI2MP7enXyKwCvMtl9p0uxd9Funy+fwMsVtULgL/n/39C2av22mtkkntZ7YnxQuAPgb/a5nrmIslTgA8Ab6qqb44/vM5T9sxrZFYB3uUy+zNjkpwLfD+7+8fHiWtSVV+vqv8e3f1j4MfnVNtOZbuGNarqm1X1n6PbHwbOS7J/m8uaqSTnsRret1TVB9cZsqdfI7MK8C6X2d8O3DC6fS1wV41+K7FLTVyTsXN3V7N6zm8vux34hdE7DV4CfKOqHtnuorZLkmee/j1Rkhez+v379e2tanZGn+u7gWNV9fYNhu3p10ifboQbqg0us0/y28DRqrqd1S/Mnyd5iNUj7+tmUctO0XFN3pjkauAUq2ty47YVPAdJbmX1nRX7kxwH3gacB1BV7wI+zOq7DB4Cvg28ZnsqnY8O63Et8EtJTgEnget2+UHPpcCrgc8muW+07TeABdibr5FxXkovSY3ySkxJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKANckhr1fyzUmmQlfg4vAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bounds = (0.1, 3.0)\n",
    "\n",
    "# Data and signal\n",
    "\n",
    "np.random.seed(0)\n",
    "tau = -2.0\n",
    "beta = -1/tau\n",
    "data = np.random.exponential(beta, 300)\n",
    "peak = np.random.normal(1.2, 0.1, 25)\n",
    "data = np.concatenate((data,peak))\n",
    "data = data[(data > bounds[0]) & (data < bounds[1])]\n",
    "\n",
    "plt.hist(data, bins=100, histtype='step');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = zfit.Space('x', limits=bounds)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From C:\\Users\\sa_li\\.conda\\envs\\rmd\\lib\\site-packages\\tensorflow\\python\\ops\\resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Colocations handled automatically by placer.\n"
     ]
    }
   ],
   "source": [
    "mean = zfit.Parameter(\"mean\", 1.2, 0.1, 2., floating=False)\n",
    "sigma = zfit.Parameter(\"sigma\", 0.1, floating=False)\n",
    "lambda_ = zfit.Parameter(\"lambda\",-2.0, -4.0, -1.0)\n",
    "Nsig = zfit.Parameter(\"Nsig\", 20., -20., len(data))\n",
    "Nbkg = zfit.Parameter(\"Nbkg\", len(data), 0., len(data)*1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal = Nsig * zfit.pdf.Gauss(obs=obs, mu=mean, sigma=sigma)\n",
    "background =  Nbkg * zfit.pdf.Exponential(obs=obs, lambda_=lambda_)\n",
    "tot_model = signal + background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the negative log likelihood\n",
    "from zfit.core.loss import ExtendedUnbinnedNLL, UnbinnedNLL\n",
    "data_ = zfit.data.Data.from_numpy(obs=obs, array=data)\n",
    "nll = ExtendedUnbinnedNLL(model=[tot_model], data=[data_], fit_range=[obs]) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load and instantiate a tensorflow minimizer\n",
    "from zfit.minimizers.minimizer_minuit import MinuitMinimizer\n",
    "minimizer = MinuitMinimizer()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td title=\"Minimum value of function\">FCN = -1145.2067313908976</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 70</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 70</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 2.0741373583835403e-06</td>\n",
       "        <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 5e-06</td>\n",
       "        <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
       "        UP = 0.5</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<table>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
       "        <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
       "        <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
       "        <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
       "        <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
       "        <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
       "        <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#UxIJHoHYmM').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>lambda</td>\n",
       "        <td>-1.93315</td>\n",
       "        <td>0.140803</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-4</td>\n",
       "        <td>-1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>19.4765</td>\n",
       "        <td>7.12684</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-20</td>\n",
       "        <td>271</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>Nbkg</td>\n",
       "        <td>251.519</td>\n",
       "        <td>16.7673</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>298.1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"UxIJHoHYmM\" style=\"display:none;\">\n",
       "<textarea rows=\"12\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & $\\lambda$ & -1.93315 & 0.140803 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & Nsig & 19.4765 & 7.12684 &  &  & -20.0 & 271 & No\\\\\n",
       "\\hline\n",
       "2 & Nbkg & 251.519 & 16.7673 &  &  & 0.0 & 298.1 & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Create the minimization graph to minimize mu and sigma and run it (minimize does it directly)\n",
    "minimum = minimizer.minimize(loss=nll)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plotfitresult(pdf, bounds, nbins, data):\n",
    "    x = np.linspace(*bounds, num=1000)\n",
    "    pdf = zfit.run(tot_model.pdf(x, norm_range=bounds) * tot_model.get_yield())\n",
    "    _ = plt.plot(x, ((bounds[1] - bounds[0])/nbins)*(pdf), \"-r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAD8CAYAAABw1c+bAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XmcFNW5//HPM+wiCMiAOIAgKtFoWJyLGK6KbAoqBKMCGhU37nWJGr1eR2LUQESuxiQq9xeCgmiuIGBUUEBE1OCGMhg2g4RFUBYBRUQBQfD5/VE9MAw9Mz3d01O9fN+vV726ltNVT9H6VM2pOueYuyMiItkjJ+wARESkainxi4hkGSV+EZEso8QvIpJllPhFRLKMEr+ISJZR4hcRyTJK/CIiWUaJX0Qky1QPO4BoGjdu7K1atQo7DBGRtLFgwYIv3D03lrIpmfhbtWpFYWFh2GGIiKQNM1sba1lV9YiIZBklfhGRLKPELyKSZZT4RUSyjBK/iEiWKTfxm1kLM3vDzJaZ2UdmdktkfSMzm21mKyKfDUv5/pWRMivM7MrKPgEREamYWO749wK3u/uJQGfgRjM7CSgA5rj78cCcyPJBzKwRcC9wGtAJuLe0C4SIiFSNchO/u2909w8j898Ay4A8oB/wVKTYU8DPonz9HGC2u29196+A2cC5lRG4iIjEp0J1/GbWCugAvA80dfeNEFwcgCZRvpIHfFZseV1knYiIhCTmlrtmdjjwN+BWd99uZjF9Lcq6qKO7m9kQYAhAy5YtYw3rEF1Gvs76bbv2L+c1qMM7Bd3i3p+ISKaJKfGbWQ2CpP+Muz8fWb3JzJq5+0YzawZsjvLVdUDXYsvNgTejHcPdxwBjAPLz86NeHGKxftsu1ow8b/9yq4Lp8e5KRCQjxfJWjwFjgWXu/odim6YBRW/pXAlMjfL1WUAvM2sYeajbK7JORERCEksdfxfgcqCbmS2MTH2AkUBPM1sB9IwsY2b5ZvYEgLtvBYYD8yPTsMg6EREJSblVPe7+NtHr6gG6RylfCFxbbHkcMC7eAEVEpHKp5a6ISJZR4hcRyTJK/CIiWUaJX0Qkyyjxi4hkGSV+EZEso8QvIpJllPhFRLKMEr+ISJZR4hcRyTJK/CIiWUaJX0Qkyyjxi4hkGSV+EZEso8QvIpJllPhFRLKMEr+ISJYpdwQuMxsHnA9sdveTI+smAW0jRRoA29y9fZTvrgG+AfYBe909v5Lijm7zZnK//SqphxARSXflJn5gPDAKeLpohbsPKJo3s4eBr8v4/tnu/kW8AcZs50445hiu+Ukf4BdJP5yISLoqt6rH3ecCUQdINzMDLgEmVnJcFXfYYXDGGfRc+UHYkYiIpLRE6/jPADa5+4pStjvwqpktMLMhCR6rfH370mbrOli+POmHEhFJV4km/kGUfbffxd07Ar2BG83szNIKmtkQMys0s8ItW7bEF03fvsHntGnxfV9EJAvEnfjNrDpwITCptDLuviHyuRl4AehURtkx7p7v7vm5ubnxBdWyJR81OVaJX0SkDInc8fcAPnb3ddE2mlldM6tXNA/0ApYmcLyYzD7+NHj3XYj3rwYRkQxXbuI3s4nAe0BbM1tnZtdENg2kRDWPmR1tZjMii02Bt81sEfABMN3dX6m80KObfdxp8MMPMH16sg8lIpKWyn2d090HlbJ+cJR1G4A+kfnVQLsE46uwj5q2gebNYepUGDy4qg8vIpLyMq/lrlnwkPfVV2HXrrCjERFJOZmX+AH69QsadM2ZE3YkIiIpJzMT/1lnQb16ertHRCSKzEz8tWpB797w0kuY/xB2NCIiKSUzEz8E9fyff06HDWrFKyJSXOYm/vPPhxo1OHf5u2FHIiKSUjI38R9xBPTqRZ/l74B72NGIiKSMzE38ABddRPPtm2HBgrAjERFJGZmd+Pv25fucavDcc2FHIiKSMjI78TdqxHstfxIkflX3iIgAmZ74gXc7nA2rVtHn6lF0Gfl62OGIiIQu4xN/wV8KICeHGXmbWL9NXTiIiGR84ic3F7p2hSlTVN0jIkI2JH6Aiy6Cf/2LE75YG3YkIiKhy47E378/mAXv9IuIZLnsSPxHHQVnnEGfj5X4RUSyI/EDDBjACV9+CkuWhB2JiEioYhl6cZyZbTazpcXW3Wdm681sYWTqU8p3zzWz5Wa20swKKjPwCrv4YvZaDkyYEGoYIiJhi+WOfzxwbpT1f3T39pFpRsmNZlYN+F+gN3ASMMjMTkok2ITk5vJ2qw4wcWIwJq+ISJYqN/G7+1xgaxz77gSsdPfV7r4HeBboF8d+Ks3Uk86CtWvhvffCDENEJFSJ1PHfZGaLI1VBDaNszwM+K7a8LrIuKjMbYmaFZla4ZcuWBMIq3avHd4batYO7fhGRLBVv4v8z0AZoD2wEHo5SxqKsK7UFlbuPcfd8d8/Pzc2NM6yy7ah1WDBAy+TJ8P33STmGiEiqiyvxu/smd9/n7j8AjxNU65S0DmhRbLk5sCGe41WqQYNgyxYNxC4iWSuuxG9mzYot9geWRik2HzjezFqbWU1gIBD+6Oe9eweDtKi6R0SyVCyvc04E3gPamtk6M7sGeNDMlpjZYuBs4FeRskeb2QwAd98L3ATMApYBk939oySdR+xq1Qq6cHj+edilTttEJPtUL6+Auw+KsnpsKWU3AH2KLc8ADnnVM3SXXgpjx8LLL8PFF4cdjYhIlcqelrvFnXUWHH00/PWvYUciIlLlsjPxV6sGl18OM2bApk1hRyMiUqWyM/EDDB4M+/bBM8+EHYmISJXK3sT/ox9B587w5JMaoEVEskr2Jn4I7vqXLoUPPww7EhGRKlPuWz0ZbcAAuOUWGD8eTj31kM1dRr5+0Di9eQ3q8E5Bt5h2nch3RUSSKbsTf4MGwehcEybA738fvONfzPptu1gz8rz9y60Kpse860S+KyKSTNld1QNBdc/WrcE7/SIiWUCJv0eP4J3+8ePDjkREpEoo8VerBldcATNnwsaNYUcjIpJ0SvwAV10VvNOvu34RyQJK/AAnnABnnw2PP65hGUUk4ynxFxkyBD75BGbPDjsSEZGkUuIv0r8/NG4MY8aEHYmISFIp8RepVSt4tXPqVD3kFZGMllUNuPIa1CmzIdVpfgqT9u2DcePg17+uwshERKpOViX+8rpMaFUwHbp1Cx7yFhRUUVQiIlUrlqEXx5nZZjNbWmzdQ2b2sZktNrMXzKxBKd9dExmicaGZFVZm4EkzZAisXauHvCKSsWKp4x8PnFti3WzgZHf/CfAv4K4yvn+2u7d39/z4Qqxi/ftDbi6MHh12JCIiSVFu4nf3ucDWEutejQymDjAPaJ6E2MJRsyZcfTW89BJ5X28OOxoRkUpXGW/1XA3MLGWbA6+a2QIzG1LWTsxsiJkVmlnhli1bKiGsBNxwAwCX/0M9aopI5kko8ZvZr4G9QGnjF3Zx945Ab+BGMzuztH25+xh3z3f3/Nzc3ETCSlzLltC/PwMXzYKdO8ONRUSkksWd+M3sSuB84DL36GMXuvuGyOdm4AWgU7zHq3I330yD777VmLwiknHiSvxmdi5wJ9DX3aPeEptZXTOrVzQP9AKWRiubks44g382aQ2PPaYxeUUko8TyOudE4D2grZmtM7NrgFFAPWB25FXN0ZGyR5vZjMhXmwJvm9ki4ANguru/kpSzSAYznjz1AliyBP7+97CjERGpNOU24HL3QVFWjy2l7AagT2R+NdAuoehCNu3Es3jog2fg0Ueha9ewwxERqRRZ1XK3onbXqBU06Pqf/4E1axLaV8nuIjT4uoiERYm/PNdfDw8+GNT114g/UZdM8hp8XUTCot45y9OiBQwYAGPGUP+7b8OORkQkYUr8sbjjDvj2Wy5dmD7PpkVESqPEH4v27aFnT64unAq7d4cdjYhIQpT4Y3XHHTTZ8ZUadIlI2lPij1WPHnzU5Fh46CENyC4iaU2JP1Zm/OW0C+Hjj2G63sgRkfSlxF8BM9r+OxxzTPB6p4hImlLir4C91arDbbfB228Hk4hIGlIDroq69lq4/34YPhxmzdq/usvI11m/bdf+5bwGdcKITkSkXEr8FXXYYXD77XDnnTBvHnTuDMD6bbtYM/K8kIMTESmfqnriccMNcOSRwV2/iEiaUeKPx+GHB3X9M2bAggVhRyMiUiFK/PG66SZo0EB3/SKSdpT441W/Ptx6K0ydCgsXhh2NiEjMlPgTcfPNwQVAd/0ikkZiSvxmNs7MNpvZ0mLrGpnZbDNbEflsWMp3r4yUWREZoD1zNGwY3PU//zwnf74y7GhERGIS6x3/eODcEusKgDnufjwwJ7J8EDNrBNwLnAZ0Au4t7QKRtm67DRo14o65T4cdiYhITGJK/O4+F9haYnU/4KnI/FPAz6J89RxgtrtvdfevgNkcegFJb0ccAXfdxVmffKhB2UUkLSTSgKupu28EcPeNZtYkSpk84LNiy+si6w5hZkOAIQAtW7ZMIKz4RRsXNyY33sjnv32Ao4YODbpyMEtShCIiiUt2y91oGdCjFXT3McAYgPz8/Khlki3uwc/r1OHRLoMYMet/g3f7z1MLXhFJXYm81bPJzJoBRD43RymzDmhRbLk5sCGBY6asyaf0hDZt4Ne/Vn/9IpLSEkn804Cit3SuBKZGKTML6GVmDSMPdXtF1mWcvdWqw7BhsGgRTJoUdjgiIqWK9XXOicB7QFszW2dm1wAjgZ5mtgLoGVnGzPLN7AkAd98KDAfmR6ZhkXWZaeDAYHzeu+6C774LOxoRkahiquN390GlbOoepWwhcG2x5XHAuLiiSzc5OfDww9C9OzzySNCDp4hIilHL3crWrRtccEHQZ//maI89RETCpcSfDA89BLt2wX33hR2JiMghlPiToW1buP56+Mtf4KOPwo5GROQgSvzJcs89UK8e3HFH2JGIiBxEiT9ZGjeG3/wGZs6El18OOxoRkf005m6cYhpc/Ze/hCeegFtugR49oHbtKoxQRCQ6Jf44xTS4es2aMGpUkPQffDCo/hERCZmqepKte3cYMAAeeABWrw47GhERJf4q8fDDUK1aMGiLiEjIlPirQl5e8E7/Sy8Fk4hIiJT4q8ott8BJJwUPfL/9NuxoRCSLKfFXlRo1YMwY+PRTuPvusKMRkSymxF+VunSBG26ARx+lw/qPw45GRLKUEn9Ve+ABaN6cka88Cnv2hB2NiGQhJf6qVq8ejB5N2y8+DS4CIiJVTA24wtCnD6+2607XYb/jvE9zWZF7TJnF8xrUiX88YBGREuJO/GbWFig+xuCxwD3u/qdiZboSDMn4SWTV8+4+LN5jZpJesyfCiScye8l4eOed4OFvKVoVTK+6wEQk48Vd1ePuy929vbu3B04FdgIvRCn6VlE5Jf1icnPhz3+G+fNhxIiwoxGRLFJZdfzdgVXuvraS9pcdLr4YfvELGD48uACIiFSBykr8A4GJpWw73cwWmdlMM/txJR0vczz2GDRrFlwAdu4MOxoRyQIJJ34zqwn0BaZE2fwhcIy7twMeA14sYz9DzKzQzAq3bNmSaFjpo0EDGD8e/vUvDc4uIlWiMu74ewMfuvumkhvcfbu7fxuZnwHUMLPG0Xbi7mPcPd/d83NzcyshrDTSvXvQgduoUTBrVtjRiEiGq4zEP4hSqnnM7Cgzs8h8p8jxvqyEY2aeESOCvnyuvBI+/zzsaEQkgyWU+M3sMKAn8Hyxdf9pZv8ZWbwIWGpmi4BHgYHu7okcM2PVqQPPPgtffx3U9+/bF3ZEIpKhEkr87r7T3Y9096+LrRvt7qMj86Pc/cfu3s7dO7v7u4kGnNFOOSV42Dtnjl7xFJGkUcvdVHPNNfDGG0H//WeeCWedRV6DOmU24krHlr3RxixOt3MQSVdK/KnGDEaPhsJCGDQIFi4sNyGmY8vekmMWp+M5iKQrddKWiurVg8mTYetWuPRS2Ls37IhEJIMo8aeqdu2CLh3mzIGCgrCjEZEMoqqeVHbVVbBgQTBYe8eOwd2/iEiCdMef6v74RzjjjOCh7z/+EXY0IpIBlPhTXY0aMGUKNG4M/fvDF1+EHZGIpDkl/nTQtCm88ELQovfCC2H37rAjEpE0psSfLvLz4ckn4a234OqrQQ2gRSROeribTgYNgtWr4e67oU0bGKZxbUSk4pT4083QoUHyHz4cjj0WBg8OOyIRSTNK/BVQvOuEvAZ1wgmiqGXv2rVw3XXQokWFvp7MrhIS2XfJbikq8l11/yBSMUr8FZAyyaRGDXjuOfj3f4ef/YxTLhwe81eT2VVCIvsu+W9bke+q+weRitHD3XTVoEEwaEvjxjw15V5YtizsiEQkTSjxp7O8PJg9m305OdCzZ1D9IyJSDiX+dHfccVx+yXDYsSNI/psOGQFTROQgSvwZ4OMmrWH6dFi3Dnr0gM2bww5JRFJYwonfzNaY2RIzW2hmhVG2m5k9amYrzWyxmXVM9JgSxU9/Ci+/DKtWBYO3K/mLSCkq647/bHdv7+75Ubb1Bo6PTEOAP1fSMaWkbt0OJP9u3ZT8RSSqqqjq6Qc87YF5QAMza1YFx81ORcl/9WolfxGJqjISvwOvmtkCMxsSZXse8Fmx5XWRdZIs3boFdf6rV0PXrkHdf7r49FN4801+tPkT+OGHsKMRyUiVkfi7uHtHgiqdG83szBLbLcp3DulhzMyGmFmhmRVu2bKlEsLKcmefDa+8AuvXQ5cusHx5lYdQ1Bq3aOoy8vXSC8+dC6efDsccA2efzStP/hJatQo6plOHdCKVKuHE7+4bIp+bgReATiWKrAOK9yvQHNgQZT9j3D3f3fNzc3MTDUsAzjwT3nwTvvsuaOW7YEGVHv6dgm6sGXne/ql4twr7/fAD3HUXnHUWbNgAv/89vPYat533K2jePOiJ9IorNO6wSCVKKPGbWV0zq1c0D/QClpYoNg24IvJ2T2fga3ffmMhxpQI6dIC334a6daFrV05fuzjsiA7YuxcGDICRI2HIkKD18e23Q/fuPH9y9yDu4cPh//4vSP6q+hGpFIn21dMUeMHMivY1wd1fMbP/BHD30cAMoA+wEtgJXJXgMaWijj8e3nkHzjmH8VPugfNbwmWXhRuTO9x0U9Dn0EMPBQnfStQK5uQEXVBXrx78VXDSScGyiCQkocTv7quBdlHWjy4278CNiRxHKkFeHsydy4cdunL6L34BK1bAvfeGFs4N86bA3Kfhzjvhv/6r7MJ33glLlsA99wTVV2eWfIwkIhWhlrvZpFEjrrhkGFx1Ffz2t3DZZdTau6fq43jmGf577tPBwDIjRpRf3gzGjAke9l57bfDMQkTipsSfZb6vVgPGjoUHHoCJE5kwcWjwULWqvP46XHUV77U8JXhjJyfG/wTr1g2S/4oVGnlMJEFK/NnIDAoK4Lnn+NGWNdCxI/z970k/7Alb1kD//nDCCfxH/19DrVoV20GPHsFD3ocfDtooiEhclPiz2c9/zs8ufzjo27979yChJuud+fXrGT/lPjj8cJgxg+21D49vPyNGBA97hw6t1PBEsokSf5ZbkXsMfPAB9OsXPGTt3x8quwHd1q3Qqxf1du8IWhS3bBn/vvLygjgnTYJ58yovRpEsoqEXM1y08WgPUb9+8Frln/7Enjvu5OtWJ3Bn75tZnn/WQUMixjW27Y4dcN55sHIlQy68j4nt2yd8HoftOYW3Dm/IkQUFQQO1FFI8To39K6lKiT/DlRyPtlRm8KtfccHSmsxa8Djj/jaMZ1f0giHtoVGjqPsqd2zbbdugb9/gL4opU3jvgwrW6ZdxHr9d/Cr3znkc3ngj7n0mQ/E4NfavpCpV9chBlue2gvffhzvv5OIlr8EJJ8ATT1S81WxRB3Hz5sGECXDhhZUa54R250KzZkFbBPXlI1IhSvxyqFq1YORIzhv8CJx4Ilx3HZx8Mv0+eqP8PnP27YNx44KuItauhZdeCrplqGS7a9QKWvO+9RY/Xbuo0vcvksmU+KVUHzdpHfSaOWkS5OTwyMsPQ4sWQfcKr71G/e++De629+2DVav4xT9mBAn/mmvg5JPhH/+Ac85JXoDXXQd5edz6zgTd9YtUgBK/lM0MLrkEFi/mugvvhtNOg8ceg549WfzIQKhdO3i98rjj+N2r/y/4zoQJQQdrrVolN7batWHoUDqt+yfMmZPcY4lkED3cldjk5DD7+M4wcjhs3w7z5vG7kZO4+98aQ506cPTR9Hx/L7OfuP7QztaS6ZprWF9wH3m/+U3QFqEqjy2SppT4peLq14devXji9e+5u9ibNitWT6/6xFurFqN+OoAHZo2CmTOhT5+qPb5IGlJVj6S9Kaf0gNat4Te/UV2/SAyU+CXt7a1WPXit88MP4cUXww5HJOUp8UtmuOyyoM3BPfdgrpG6RMqiOn6JW9Fg6sWXYy1fXtkKq149GGNg0CDOb/0WcEGpRUt2PVEyRnWzIJku7sRvZi2Ap4GjgB+AMe7+SIkyXYGpwCeRVc+7uzpTzxAVTZBJT6iXXAL33x+81//9/VCjRtRiZXVjoW4WJBskUtWzF7jd3U8EOgM3mtlJUcq95e7tI5OSviRPTg7cfz9ttq6HP/857GhEUlbcid/dN7r7h5H5b4BlQF5lBSYSlwsuYG6rDsHD3sruXlokQ1TKw10zawV0AN6Psvl0M1tkZjPN7MeVcTyRUpkxrPt18M03weudInKIhBO/mR0O/A241d23l9j8IXCMu7cDHgNKfdfOzIaYWaGZFW7RnZokYGXjlnDjjcEYvfPnhx2OSMpJKPGbWQ2CpP+Muz9fcru7b3f3byPzM4AaZtY42r7cfYy757t7fm5ubiJhiQQDsh99NAweDN99F3Y0Iikl7sRvZgaMBZa5+x9KKXNUpBxm1ilyvC/jPaZIzI44Ah5/HP75z+A1TxHZL5H3+LsAlwNLzGxhZN1QoCWAu48GLgKuN7O9wC5goLva1EsV6d076CL6wQehWzfo2TPsiERSQtyJ393fBsrskcvdRwGj4j2GSMIeeSQYUWzQICgsTH5X0SJpQC13M0DJFrQlt1Xku5XeojZO5cVVVivgkt897Zw7mPTEzdCvX4UHZ49rgPkkKKu1MajFsVSMEn8GSOR/+FRNFuXFVdb2kttaFUyHKVPg/PPhvPM4rPNtMcdR4QHmk6Ss1sagFsdSMeqkTbJDz57w7LPwwQc8PfkeNe6SrKbEL9mjf3+YPJmTN62Czp1h6dKwIxIJhRK/ZJcLL2TgoAdgxw7Iz4c//CEYLF4kiyjxS9ZZeHRbWLwYzj0Xbr8dOnSAV17R6F2SNZT4JTs1aQIvvACTJ8O33wbv/HfsyIVL58DOnWFHJ5JUSvySvczg4oth2bKgX589e/jD9D8GF4VLL4UXX+Tw3boISObR65witWrBddfBtdcy8LKRPFt/LTz3HEycyELLgfmPBG8FdemiC4FkBN3xixQxY17Ln8Do0bBxI7zxBqM7XwR79gT9/fTsyeI/DYCTT4Zrrw36Apo/H3aV3rBKJBXpjl9KlaqtepPpkHPu+x/cVNANtm2D+fMZ++AE2qxeSodnJtNw7FgA9lkOaxoezbImrbmr5fHwskO7dgk/LC7eWjeeFtjFG7KV3FdlNdyraMvmVGkJne2U+KVU2fg/ZKnn3KAB9OzJdUUdvbnD6tWweDHVFi2izaJFtFm8GF59C14dB8DSmnVgzknwox8dPB13HNSuXW4s5bXWLSvuki15i++rMlv5VrRlc6q0hM52Svwi8TCDNm2CqX//A+u3b4clS2DJEiY/MZOrG++Gt9+GZ545UCYnB1q3huOPh2OPPTC1bh181q9f9ecjWUWJX6Qy1a8PXbpAly4MW9OCq4vubnfsgBUr4OOPg2nZMli5Mug59KuvDt7HkUfCscfy2De1wd6GFi2gefMDn40bBxcekTgp8YtUhbp1oX37YCrpq6/gk0+CqqOi6ZNPOHnlEvj9e7B378Hla9WCvLwDF4KiqVkzaNoUmjYN3j5y1wVColLiFwlbw4bB1LHjQavPLpjOmhG9YdMmWLfu4Omzz4LPd98NPr///qDvLgUYXRuOOgqaNuXxL4GtU6FpU65YsAUm74BGjYK/Loo+69bVhSJLKPGLpLKcnOBOvlkz+Ld/i17mhx+C3kY//zyYNm1ixJNvMvTURsFFY9Mmmq9eAdOmwZYtDPvhB3jtL4fup2bN4CJQdCEoflFo1Ch4wH3EEUF1VuQz7+vNwV8s9etDtWrJ/beQSpNQ4jezc4FHgGrAE+4+ssT2WsDTwKkEY+0OcPc1iRxTRErIydlfxUO7dgCM+eeRDC329kzvgunB2zT79pF/60QK/6M9bN0aTF9+Gf1z1aqgncKXX5Y6YP07AKOvDhbq1mWe1YKpTQ5cIIqmunWhbl1ufmcdPPzx/uVzln8Ms6rvX6ZuXTjssAPzNWsm+R8vO8Wd+M2sGvC/QE9gHTDfzKa5+z+LFbsG+MrdjzOzgcD/AAMSCVhEElCtGl/UbRg0QquInTvh66+Dafv2/Z93jH2Lh3q12r/ujdeWMqjtEQfKrF8fzO/cCTt2cNvu3fD2gTec/gLw4ojSj1s9clGoUyd4BbZWreCz+BRtXVlla9YMpho1gqmi89Wrp32VWCJ3/J2Ale6+GsDMngX6AcUTfz/gvsj8c8AoMzMNuC6SZg47LJiaNTto9ZT5tXno1gN/WdxVMJ1BZbQ9aPPf01h1d9fgLacdO+g9YiYzrz11/4Wh1Gn37uCvjpLTtm2lbyvlr5RKUdbFoXr1oNqrevWKzzdsCI8+mry4IxJJ/HnAZ8WW1wGnlVbG3fea2dfAkcAXCRxXRNLUvpxqB6p/gGVNjoWf/jQ5B3MPHnoXvxDs2hWs27Mn+EzG/L59wZtYRZ/F53fvDi5yRetLljvyyOT8W5Rg8d58m9nFwDnufm1k+XKgk7v/sliZjyJl1kWWV0XKfBllf0OAIZHFtsDyYpsbk3kXi0w7p0w7H8i8c8q084HMO6dEzucYd8+NpWAid/zrgBbFlpsDG0ops87MqgNHAFuj7czdxwBjom0zs0J3z08g1pSTaeeUaecDmXdOmXY+kHnnVFXnk0jvnPOB482stZnVBAYC00qUmQZcGZm/CHhzu7RLAAADrElEQVRd9fsiIuGK+44/Umd/EzCL4HXOce7+kZkNAwrdfRowFvirma0kuNMfWBlBi4hI/BJ6j9/dZwAzSqy7p9j8d8DFiRwjImoVUJrLtHPKtPOBzDunTDsfyLxzqpLzifvhroiIpCeNwCUikmVSKvGb2blmttzMVppZQZTttcxsUmT7+2bWquqjjF0M5zPYzLaY2cLIdG0YccbKzMaZ2WYzW1rKdjOzRyPnu9jMOkYrl0piOKeuZvZ1sd/onmjlUoWZtTCzN8xsmZl9ZGa3RCmTNr9TjOeTbr9RbTP7wMwWRc7pt1HKJDfXuXtKTAQPiFcBxwI1gUXASSXK3ACMjswPBCaFHXeC5zMYGBV2rBU4pzOBjsDSUrb3AWYCBnQG3g875ko4p67Ay2HHWYHzaQZ0jMzXA/4V5b+7tPmdYjyfdPuNDDg8Ml8DeB/oXKJMUnNdKt3x7+8Cwt33AEVdQBTXD3gqMv8c0N0sZTvNiOV80oq7z6WUdhgR/YCnPTAPaGBmzcooH7oYzimtuPtGd/8wMv8NsIygBX1xafM7xXg+aSXy7/5tZLFGZCr5sDWpuS6VEn+0LiBK/sAHdQEBFHUBkYpiOR+An0f+3H7OzFpE2Z5OYj3ndHN65M/ymWb247CDiVWkeqADwR1lcWn5O5VxPpBmv5GZVTOzhcBmYLa7l/obJSPXpVLij3Y1K3kVjKVMqogl1peAVu7+E+A1Dlzh01U6/T6x+pCgKXw74DHgxZDjiYmZHQ78DbjV3beX3BzlKyn9O5VzPmn3G7n7PndvT9DjQSczK9ldalJ/o1RK/BXpAoLyuoBIAeWej7t/6e67I4uPE4xbkM5i+Q3TirtvL/qz3IN2KzXMrHHIYZXJzGoQJMln3P35KEXS6ncq73zS8Tcq4u7bgDeBc0tsSmquS6XEn2ldQJR7PiXqVfsS1F+ms2nAFZG3RjoDX7v7xrCDSoSZHVVUt2pmnQj+nzmkk8FUEYl1LLDM3f9QSrG0+Z1iOZ80/I1yzaxBZL4O0AP4uESxpOa6lBl60TOsC4gYz+dmM+sL7CU4n8GhBRwDM5tI8AZFYzNbB9xL8GAKdx9N0Iq7D7AS2AlcFU6ksYvhnC4CrjezvcAuYGAK32wAdAEuB5ZE6pABhgItIS1/p1jOJ91+o2bAUxYMZpUDTHb3l6sy16nlrohIlkmlqh4REakCSvwiIllGiV9EJMso8YuIZBklfhGRLKPELyKSZZT4RUSyjBK/iEiW+f/20BvYUUtQiAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "nbins = 80\n",
    "plt.hist(data, bins=nbins, histtype='step', range=bounds);\n",
    "plotfitresult(tot_model, bounds, nbins, data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Discovery significance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lossbuilder(model, data, weights=None):\n",
    "    loss = ExtendedUnbinnedNLL(model=model, data=data, fit_range=[obs]) \n",
    "    return loss\n",
    "    \n",
    "config = Config(tot_model, data_, lossbuilder, MinuitMinimizer())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Get fit best values!\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td title=\"Minimum value of function\">FCN = -1145.2067314770634</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 36</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 36</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 1.9878782071289407e-06</td>\n",
       "        <td title=\"Maximum EDM definition of convergence\">GOAL EDM = 5e-06</td>\n",
       "        <td title=\"Error def. Amount of increase in FCN to be defined as 1 standard deviation\">\n",
       "        UP = 0.5</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<table>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Validity of the migrad call\">Valid</td>\n",
       "        <td align=\"center\" title=\"Validity of parameters\">Valid Param</td>\n",
       "        <td align=\"center\" title=\"Is Covariance matrix accurate?\">Accurate Covar</td>\n",
       "        <td align=\"center\" title=\"Positive definiteness of covariance matrix\">PosDef</td>\n",
       "        <td align=\"center\" title=\"Was covariance matrix made posdef by adding diagonal element\">Made PosDef</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" title=\"Was last hesse call fail?\">Hesse Fail</td>\n",
       "        <td align=\"center\" title=\"Validity of covariance\">HasCov</td>\n",
       "        <td align=\"center\" title=\"Is EDM above goal EDM?\">Above EDM</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" title=\"Did last migrad call reach max call limit?\">Reach calllim</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">True</td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "        <td align=\"center\"></td>\n",
       "        <td align=\"center\" style=\"background-color:#92CCA6\">False</td>\n",
       "    </tr>\n",
       "</table>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<table>\n",
       "    <tr>\n",
       "        <td><a href=\"#\" onclick=\"$('#LsvqnqDxbE').toggle()\">+</a></td>\n",
       "        <td title=\"Variable name\">Name</td>\n",
       "        <td title=\"Value of parameter\">Value</td>\n",
       "        <td title=\"Hesse error\">Hesse Error</td>\n",
       "        <td title=\"Minos lower error\">Minos Error-</td>\n",
       "        <td title=\"Minos upper error\">Minos Error+</td>\n",
       "        <td title=\"Lower limit of the parameter\">Limit-</td>\n",
       "        <td title=\"Upper limit of the parameter\">Limit+</td>\n",
       "        <td title=\"Is the parameter fixed in the fit\">Fixed?</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>0</td>\n",
       "        <td>lambda</td>\n",
       "        <td>-1.93314</td>\n",
       "        <td>0.140801</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-4</td>\n",
       "        <td>-1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>1</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>19.4765</td>\n",
       "        <td>7.12641</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-20</td>\n",
       "        <td>271</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>Nbkg</td>\n",
       "        <td>251.519</td>\n",
       "        <td>16.7609</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>298.1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"LsvqnqDxbE\" style=\"display:none;\">\n",
       "<textarea rows=\"12\" cols=\"50\" onclick=\"this.select()\" readonly>\n",
       "\\begin{tabular}{|c|r|r|r|r|r|r|r|c|}\n",
       "\\hline\n",
       " & Name & Value & Hesse Error & Minos Error- & Minos Error+ & Limit- & Limit+ & Fixed?\\\\\n",
       "\\hline\n",
       "0 & $\\lambda$ & -1.93314 & 0.140801 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & Nsig & 19.4765 & 7.12641 &  &  & -20.0 & 271 & No\\\\\n",
       "\\hline\n",
       "2 & Nbkg & 251.519 & 16.7609 &  &  & 0.0 & 298.1 & No\\\\\n",
       "\\hline\n",
       "\\end{tabular}\n",
       "</textarea>\n",
       "</pre>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<hr>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<zfit.minimizers.fitresult.FitResult at 0x1a44fc99b0>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config.bestfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = FrequentistCalculator(config, ntoysnull=5000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Toys successfully read from 'toys_Disco_Nsig.hdf5' !\n"
     ]
    }
   ],
   "source": [
    "calc.readtoys_from_hdf5(Nsig, \"toys_Disco_Nsig.hdf5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "poinull = POI(Nsig, value=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "discovery_test = Discovery(poinull, calc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute qobs for the null hypothesis!\n",
      "\n",
      "p_value for the Null hypothesis = 0.0008\n",
      "Significance = 3.155906757921808\n"
     ]
    }
   ],
   "source": [
    "discovery_test.result();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "#calc.toys_to_hdf5(\"toys_Disco_Nsig.hdf5\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute qobs for the null hypothesis!\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEKCAYAAAAb7IIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAFSpJREFUeJzt3X+sXOWd3/H3t7aDCTj4YhtKbcM1yCFlg9axLYfGBdxACYYsBNQ0WFWxQiQXFSugtBJOkfiRFRXbslsXRKnYxsEkNJDtQuM03hKHkJgkcsAmDjELrE3qhAsu9poYY3504/DtH3OuM7meuT/nztx7n/dLGs3Mc54z853j4/O55zlnzkRmIkkqz9/rdAGSpM4wACSpUAaAJBXKAJCkQhkAklQoA0CSCmUASFKhDABJKpQBIEmFmtzpAvozc+bM7O7uHt7ML75Yuz/zzJbVI0njwbZt2/42M2cN1G9MB0B3dzdbt24d3szLltXuv//9VpUjSeNCRPxyMP0cApKkQhkAklQoA0CSCjWmjwFI0lD85je/oaenh3fffbfTpbTF1KlTmTNnDlOmTBnW/AaApAmjp6eHadOm0d3dTUR0upxRlZns37+fnp4e5s2bN6zXcAhI0oTx7rvvMmPGjAm/8QeICGbMmDGivR0DQNKEUsLGv9dIP6sBIEmFMgAkqYOWLVs2/C+8jtCEPgj8zK8OcOWabx/VPnv6sfxozcc7UJEkjR0TOgD+7vBv2X3HpUe1dzcIBUlqldtvv50HHniAuXPnMmvWLBYtWsSFF17Itddey9tvv80ZZ5zBunXr6OrqAuBrX/san//85zl48CDr1q1jyZIl/OAHP+D6668HamP9mzdvZtq0aS2tc0IHgKSC3XADbN/e2tdcsADWru23y7Zt23jooYf46U9/yuHDh1m4cCGLFi3i6quv5u677+b888/n5ptv5rbbbmNt9VpvvfUWP/7xj9m8eTPXXHMNO3bs4M477+See+5h6dKlHDp0iKlTp7b2s+AxAElqqSeffJIrrriC97///XzgAx/gsssu46233uLAgQOcf/75AKxcuZLNmzcfmWfFihUAnHfeeRw8eJADBw6wdOlSvvCFL3DXXXdx4MABJk9u/d/r7gFImpgG+Et9NA319My+/SOCNWvWcOmll7Jx40bOOeccvvvd7/KhD32olWW6ByBJrXTeeefx6KOP8s477/Dmm2/yrW99i+OOO46uri6efPJJAL761a8e2RsAePjhhwH44Q9/yAknnMAJJ5zASy+9xNlnn82NN97I4sWLeeGFF1peq3sAktRCCxcu5DOf+QwLFizgtNNO49xzzwVg/fr1Rw4Cn3766XzlK185Mk9XVxcf+9jHjhwEBli7di1PPPEEkyZN4qyzzmL58uUtr9UAkKQWu+mmm7jpppsAuPXWWwFYsGABW7ZsOarv95v8aNXdd989WuUd4RCQJBXKPQBJGkW9ewBjkXsAklQoA0CSCmUASFKhDABJKpQHgSVNWEvv+B6vHHinZa/X7isJd3d3s3XrVmbOnMnxxx/PoUOHWvr6BoCkCeuVA+80vCLwcE20Kwk7BCRJLXb77bdz5plncuGFF7JixQruvPNOli1bxo033siSJUv44Ac/eOSyEPfffz+rV68+Mu8nP/nJpl8OazUDQJJaqP5y0I888ghPP/30kWmHDx/mqaeeYu3atdx2220drLJmwACIiLkR8UREPB8Rz0XE9VX7iRGxKSJ2VvddVXtExF0RsSsino2IhXWvtbLqvzMiVo7ex5Kkzmh0OeheV155JQCLFi1i9+7dHarwdwazB3AY+DeZ+Q+Bc4DrIuIsYA3weGbOBx6vngMsB+ZXt1XAvVALDOAW4KPAEuCW3tCQpImk2eWgjznmGAAmTZrE4cOHAZg8eTLvvffekT7vvvvu6BdYGTAAMnNPZj5TPX4TeB6YDVwOrK+6rQc+VT2+HHgga7YA0yPiFOATwKbMfD0zfw1sAi5u6aeRpA5rdDno/nR3d7N9+3bee+89Xn75ZZ566qk2VTrEs4Aiohv4CPAT4OTM3AO1kIiIk6pus4GX62brqdqatUvSqJg9/diWnrkze/qxA/ZpdjnoZpYuXcq8efM4++yz+fCHP8zChQv77d9Kgw6AiDge+Evghsw82M8v3jSakP20932fVdSGjjj11FMHW54kHaWd5+zXa3Q56Poze2bOnHnkGEBE8OCDDzZ8nfrjBK3+DgAM8iygiJhCbeP/YGY+UjW/Vg3tUN3vrdp7gLl1s88BXu2n/fdk5n2ZuTgzF8+aNWson0WSNAQD7gFE7U/9LwPPZ+af1U3aAKwE7qjuv1nXvjoiHqJ2wPeNaojoMeDf1x34vQj4Yms+hiSNTWP5ctCDGQJaCvxL4OcRsb1q+3fUNvzfiIjPAb8CPl1N2whcAuwC3gY+C5CZr0fEHwO9J8V+KTNfb8mnkKRKZg75R9nHq8yjRtGHZMAAyMwf0nj8HuCCBv0TuK7Ja60D1g2lQEkarKlTp7J//35mzJgx4UMgM9m/fz9Tp04d9mt4LSBJE8acOXPo6elh3759nS6lLaZOncqcOXOGPb8BIGnCmDJlCvPmzet0GeOG1wKSpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKZQBIUqEMAEkqlAEgSYUyACSpUAaAJBXKAJCkQhkAklQoA0CSCmUASFKhDABJKpQBIEmFMgAkqVAGgCQVygCQpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKNWAARMS6iNgbETvq2m6NiFciYnt1u6Ru2hcjYldEvBgRn6hrv7hq2xURa1r/USRJQzGYPYD7gYsbtP+nzFxQ3TYCRMRZwFXAH1Tz/JeImBQRk4B7gOXAWcCKqq8kqUMmD9QhMzdHRPcgX+9y4KHM/H/A/4mIXcCSatquzPwFQEQ8VPX96yFXLElqiZEcA1gdEc9WQ0RdVdts4OW6Pj1VW7N2SVKHDDcA7gXOABYAe4A/rdqjQd/sp/0oEbEqIrZGxNZ9+/YNszxJ0kCGFQCZ+Vpm/jYz3wP+nN8N8/QAc+u6zgFe7ae90Wvfl5mLM3PxrFmzhlOeJGkQhhUAEXFK3dMrgN4zhDYAV0XEMRExD5gPPAU8DcyPiHkR8T5qB4o3DL9sSdJIDXgQOCK+DiwDZkZED3ALsCwiFlAbxtkN/CuAzHwuIr5B7eDuYeC6zPxt9TqrgceAScC6zHyu5Z9GkjRogzkLaEWD5i/30/924PYG7RuBjUOqTpI0avwmsCQVygCQpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKZQBIUqEMAEkqlAEgSYUyACSpUAaAJBXKAJCkQhkAklQoA0CSCmUASFKhDABJKpQBIEmFMgAkqVAGgCQVygCQpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKZQBIUqEMAEkq1IABEBHrImJvROyoazsxIjZFxM7qvqtqj4i4KyJ2RcSzEbGwbp6VVf+dEbFydD6OJGmwBrMHcD9wcZ+2NcDjmTkfeLx6DrAcmF/dVgH3Qi0wgFuAjwJLgFt6Q0OS1BkDBkBmbgZe79N8ObC+erwe+FRd+wNZswWYHhGnAJ8ANmXm65n5a2ATR4eKJKmNhnsM4OTM3ANQ3Z9Utc8GXq7r11O1NWs/SkSsioitEbF13759wyxPkjSQVh8EjgZt2U/70Y2Z92Xm4sxcPGvWrJYWJ0n6neEGwGvV0A7V/d6qvQeYW9dvDvBqP+2SpA4ZbgBsAHrP5FkJfLOu/erqbKBzgDeqIaLHgIsioqs6+HtR1SZJ6pDJA3WIiK8Dy4CZEdFD7WyeO4BvRMTngF8Bn666bwQuAXYBbwOfBcjM1yPij4Gnq35fysy+B5YlSW00YABk5oomky5o0DeB65q8zjpg3ZCqkySNGr8JLEmFMgAkqVAGgCQVygCQpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKZQBIUqEG/D2AiWj29GPpXvPthu0/WvPxDlQkSe1XZAA028g3CgVJmqgcApKkQhkAklQoA0CSCmUASFKhDABJKpQBIEmFMgAkqVAGgCQVygCQpEIZAJJUKANAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFcoAkKRCGQCSVCgDQJIKZQBIUqFGFAARsTsifh4R2yNia9V2YkRsioid1X1X1R4RcVdE7IqIZyNiYSs+gCRpeFqxB/BPMnNBZi6unq8BHs/M+cDj1XOA5cD86rYKuLcF7y1JGqbRGAK6HFhfPV4PfKqu/YGs2QJMj4hTRuH9JUmDMNIASOA7EbEtIlZVbSdn5h6A6v6kqn028HLdvD1VmySpAyaPcP6lmflqRJwEbIqIF/rpGw3a8qhOtSBZBXDqqaeOsDxJUjMj2gPIzFer+73Ao8AS4LXeoZ3qfm/VvQeYWzf7HODVBq95X2YuzszFs2bNGkl5kqR+DDsAIuK4iJjW+xi4CNgBbABWVt1WAt+sHm8Arq7OBjoHeKN3qEiS1H4jGQI6GXg0Inpf579n5v+OiKeBb0TE54BfAZ+u+m8ELgF2AW8Dnx3Be4+K2dOPpXvNtxu2/2jNxztQkSSNnmEHQGb+AvjDBu37gQsatCdw3XDfrx2abeQbhYIkjXd+E1iSCmUASFKhDABJKpQBIEmFMgAkqVAGgCQVygCQpEIZAJJUKANAkgplAEhSoQwASSrUSH8PoAheJE7SRGQADIIXiZM0ETkEJEmFMgAkqVAGgCQVygCQpEJ5EHgEPDtI0nhmAIyAZwdJGs8cApKkQhkAklQoA0CSCmUASFKhDABJKpRnAY0CTw+VNB4YAKPA00MljQcGQBu5ZyBpLDEA2sg9A0ljiQeBJalQBoAkFcohoDHAYwOSOsEAGAM8NiCpEwyAMcw9A0mjyQAYw9wzkDSaDIBxqNmeQe809w4kDYYBMA71t4F370DSYBkAE4zHDSQNVtsDICIuBv4zMAn4b5l5R7trmMiabeSX3vG9Ie0dNAuMpXd8j1cOvDPo/pLGrrYGQERMAu4B/inQAzwdERsy86/bWUeJhrpxbhYYs6cfy+47Lj2q3aEnafxp9x7AEmBXZv4CICIeAi4HDIAxZqiBMdShJ/ckpM5rdwDMBl6ue94DfLTNNWgUDHXoqdmexFCHqlrJsFJp2h0A0aAtf69DxCpgVfX0UES8OIL3m0nE345g/naaCYyXWmGE9f4SiC+2rphBGLDeodY0yp9hPK0P46lWKKPe0wbTqd0B0APMrXs+B3i1vkNm3gfc14o3i4itmbm4Fa812sZTrWC9o2081TueagXrrdfuq4E+DcyPiHkR8T7gKmBDm2uQJNHmPYDMPBwRq4HHqJ0Gui4zn2tnDZKkmrZ/DyAzNwIb2/R2LRlKapPxVCtY72gbT/WOp1rBeo+IzBy4lyRpwvEXwSSpUOM+ACLi4oh4MSJ2RcSaBtOPiYiHq+k/iYju9ld5pJa5EfFERDwfEc9FxPUN+iyLiDciYnt1u7kTtdbVszsifl7VsrXB9IiIu6rl+2xELOxEnVUtZ9Ytt+0RcTAibujTp6PLNyLWRcTeiNhR13ZiRGyKiJ3VfVeTeVdWfXZGxMoO1fofI+KF6t/60YiY3mTeftebNtZ7a0S8UvfvfUmTefvdjrSx3ofrat0dEdubzNua5ZuZ4/ZG7UDyS8DpwPuAnwFn9enzr4H/Wj2+Cni4g/WeAiysHk8D/qZBvcuA/9XpZVtXz25gZj/TLwH+itp3PM4BftLpmuvWjf8LnDaWli9wHrAQ2FHX9h+ANdXjNcCfNJjvROAX1X1X9birA7VeBEyuHv9Jo1oHs960sd5bgX87iHWl3+1Iu+rtM/1PgZtHc/mO9z2AI5eWyMy/A3ovLVHvcmB99fh/ABdERKMvpI26zNyTmc9Uj98Enqf27ejx7HLggazZAkyPiFM6XRRwAfBSZv6y04XUy8zNwOt9muvX0fXApxrM+glgU2a+npm/BjYBF49aoTSuNTO/k5mHq6dbqH2XZ0xosmwHYzDbkZbrr95qG/XPga+PZg3jPQAaXVqi7wb1SJ9qxX0DmNGW6vpRDUV9BPhJg8n/KCJ+FhF/FRF/0NbCjpbAdyJiW/Ut7b4G82/QCVfR/D/PWFq+ACdn5h6o/ZEAnNSgz1hcztdQ2/trZKD1pp1WV0NW65oMr43FZXsu8Fpm7mwyvSXLd7wHwICXlhhkn7aKiOOBvwRuyMyDfSY/Q23Y4g+Bu4H/2e76+liamQuB5cB1EXFen+ljcfm+D7gM+IsGk8fa8h2sMbWcI+Im4DDwYJMuA6037XIvcAawANhDbVilrzG1bCsr6P+v/5Ys3/EeAANeWqK+T0RMBk5geLuJLRERU6ht/B/MzEf6Ts/Mg5l5qHq8EZgSETPbXGZ9Pa9W93uBR6ntLtcbzL9Buy0HnsnM1/pOGGvLt/Ja77BZdb+3QZ8xs5yrA9CfBP5FVgPSfQ1ivWmLzHwtM3+bme8Bf96kjjGzbOHIdupK4OFmfVq1fMd7AAzm0hIbgN4zJv4Z8L1mK+1oq8b1vgw8n5l/1qTP3+89RhERS6j9G+1vX5W/V8txETGt9zG1A4A7+nTbAFxdnQ10DvBG73BGBzX962ksLd869evoSuCbDfo8BlwUEV3VMMZFVVtbRe0HnW4ELsvMt5v0Gcx60xZ9jkdd0aSOsXaJmguBFzKzp9HEli7f0T7SPdo3ameh/A21o/g3VW1foraCAkylNhSwC3gKOL2Dtf5jaruWzwLbq9slwLXAtVWf1cBz1M5E2AJ8rIP1nl7V8bOqpt7lW19vUPuRn5eAnwOLO7w+vJ/aBv2EurYxs3ypBdMe4DfU/vL8HLVjUo8DO6v7E6u+i6n9al7vvNdU6/Eu4LMdqnUXtfHy3vW39wy7fwBs7G+96VC9X63Wy2epbdRP6Vtv9fyo7Ugn6q3a7+9dX+v6jsry9ZvAklSo8T4EJEkaJgNAkgplAEhSoQwASSqUASBJhTIAJKlQBoAkFartPwkpjXfVdXCupvaFqH3Atsy8s7NVSUNnAEhDEBGLqF0q4CPU/v88A2zraFHSMBkA0tCcCzya1XVwIqKT14yRRsRjANLQef0UTQgGgDQ0m4ErIuLY6oqMf9TpgqThcghIGoLMfCYiHqZ2JcxfAk92uCRp2LwaqDQCEXErcMizgDQeOQQkSYVyD0CSCuUegCQVygCQpEIZAJJUKANAkgplAEhSoQwASSrU/wcHFg5kpyxMKAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "discovery_test.plot_qdist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}