examples/notebooks/discovery_zfit_asy.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "from lauztat.parameters import POI\n",
    "from lauztat.hypotests import Discovery\n",
    "from lauztat.calculators import AsymptoticCalculator\n",
    "from lauztat.config import Config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "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": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADvFJREFUeJzt3W2MXOdZxvHrIk6o06a44C0NtodNUWUpVEUNoyjFUoliAiauYiSC5EgNSUm1Eog0RaCypYJIfFoJVCgUKVoakwCRG8kNYJoUEupWEVFqYbtO49SlDcWkdgN2EjkpqkWwuPmws+5kdsdz3uZl7/n/pNXOzHn2PPc+O3Pp7Nk59zoiBABY+75v3AUAAJpBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACSxbpSTbdy4MWZnZ0c5JQCseYcPH34xImYGjRtpoM/OzurQoUOjnBIA1jzb/1FkHKdcACAJAh0AkiDQASAJAh0AkiDQASCJgYFue4/t07aPrbLtt2yH7Y3DKQ8AUFSRI/T7Je3ofdD2Fkk3Snq+4ZoAABUMDPSIeELSy6ts+iNJH5HE/7ADgAlQ6Ry67ZslnYqIpxuuBwBQUekrRW1fLuljkn624Pg5SXOS1Gq1yk53wbaFAzp19pwkadOG9Xpy/obK+wKAjKocof+YpKskPW37hKTNko7YfttqgyNiMSLaEdGemRnYiqCvU2fP6cTCTp1Y2Hkh2AEA31P6CD0inpH01uX7nVBvR8SLDdYFACipyNsW90p6StJW2ydt3zn8sgAAZQ08Qo+IWwdsn22sGgBAZVwpCgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJDAx023tsn7Z9rOuxP7D9Ndtfsf03tjcMt0wAwCBFjtDvl7Sj57HHJb0zIt4l6euSPtpwXQCAkgYGekQ8Ienlnscei4jznbtfkrR5CLUBAEpY18A+fkXSQ/022p6TNCdJrVargekm37aFAzp19pwkadOG9Xpy/oYxVwRgGtT6o6jtj0k6L+nBfmMiYjEi2hHRnpmZqTPdmnHq7DmdWNipEws7LwQ7AAxb5SN027dLep+k7RERzZUEAKiiUqDb3iHptyX9dER8t9mSAABVFHnb4l5JT0naavuk7TslfVLSFZIet33U9r1DrhMAMMDAI/SIuHWVh+8bQi0AgBq4UhQAkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0Akmiife5Ydbeq7UbbWgDTZs0H+nKr2l6z84+MoRoAGB9OuQBAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgQ6ACRBoANAEgMD3fYe26dtH+t67AdtP277G53PbxlumQCAQYocod8vaUfPY/OSPh8R75D0+c59AMAYDQz0iHhC0ss9D++S9EDn9gOSfqHhugAAJVU9h/7DEfGCJHU+v7W5kgAAVQz9j6K252wfsn3ozJkzw54OAKZW1UD/L9tXSlLn8+l+AyNiMSLaEdGemZmpOB0AYJCqgb5f0u2d27dL+rtmygEAVFXkbYt7JT0laavtk7bvlLQg6Ubb35B0Y+c+AGCMBv7Hooi4tc+m7Q3XAgCogStFASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0Akhh46f8k2rRhvWbnH7lwe5L11vrk/A1jrghAVmsy0NdSKHbXuhzsADAMnHIBgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCRqBbrt37D9rO1jtvfafkNThQEAyqkc6LY3SfqQpHZEvFPSJZJ2N1UYAKCcuqdc1klab3udpMslfbt+SQCAKio354qIU7b/UNLzks5JeiwiHusdZ3tO0pwktVqtqtPVsm3hgE6dPSfp4h0Pi44bt7VSJ4DRqnPK5S2Sdkm6StKPSHqj7ff3jouIxYhoR0R7ZmameqU1nDp7TicWdurEws4LQVhn3LitlToBjFadUy4/I+nfI+JMRPyvpIcl/VQzZQEAyqoT6M9Lus725bYtabuk482UBQAoq3KgR8RBSfskHZH0TGdfiw3VBQAoqdZ/LIqIeyTd01AtAIAauFIUAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJKodaXoWrRpw3rNzj9y4Xa/1rNFx1VFC1wATZu6QO8OzuXArjOuquUWuMPaP4DpwykXAEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiiVqDb3mB7n+2v2T5u+z1NFQYAKKfupf+fkPQPEXGL7cskXd5ATQCACioHuu03S3qvpDskKSJek/RaM2UBAMqqc8rl7ZLOSPoL21+2/Snbb2yoLgBASXVOuayTdI2kuyLioO1PSJqX9Lvdg2zPSZqTpFarVWO65nW3yF2+X0bZFri9LXmb2CcALKsT6CclnYyIg537+7QU6K8TEYuSFiWp3W5HjfkaVzcsy7bALTIfbXUBVFX5lEtE/Kekb9ne2nlou6SvNlIVAKC0uu9yuUvSg513uHxT0gfqlwQAqKJWoEfEUUnthmoBANTAlaIAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJ1O3lMrGKtKqtuq9RtbQd17yZ0Z4YmaUN9CZfqN37GmVL23HNmxntiZEZp1wAIAkCHQCSINABIAkCHQCSINABIAkCHQCSINABIAkCHQCSINABIInagW77Ettftv3ZJgoCAFTTxBH63ZKON7AfAEANtQLd9mZJOyV9qplyAABV1T1C/2NJH5H0fw3UAgCooXK3Rdvvk3Q6Ig7bvv4i4+YkzUlSq9WqOt3EaLItbxP77P7ai42p031ylC1nu+fqVnZe2uRiGtVpn7tN0s22b5L0Bklvtv3XEfH+7kERsShpUZLa7XbUmG8iDCMY6uyzyNfWbRM7ypaz3XN1KzsvbXIxjSqfcomIj0bE5oiYlbRb0oHeMAcAjA7vQweAJBr5j0UR8UVJX2xiXwCAajhCB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkGrlSFMPpwjjK/TdpGJ0O+3WV7N5/77yrfW2RtStaf1PfJ50h0RQCvSHDfhGupRf5MDod9vv+u/ffr1Nj2bUrWn9T3yedIdEUTrkAQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkUTnQbW+x/QXbx20/a/vuJgsDAJRTp5fLeUm/GRFHbF8h6bDtxyPiqw3VBgAoofIRekS8EBFHOre/I+m4pE1NFQYAKKeRbou2ZyW9W9LBVbbNSZqTpFar1cR0aEh329aLuVgr2kGta/vNN2ltYou03u1VZFy/77PIehWtddA69vs5j/pnMCl1ZFY70G2/SdJnJH04Il7t3R4Ri5IWJandbkfd+dCcfu1mL6b7hVekdW2/+SatTWzd1rtF2vsWma+IsuvYb65R/wwmpY7Mar3LxfalWgrzByPi4WZKAgBUUeddLpZ0n6TjEfHx5koCAFRR5wh9m6TbJN1g+2jn46aG6gIAlFT5HHpE/LMkN1gLAKAGrhQFgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQIdABIgkAHgCQaaZ+LyXax1q5N7bffvvqNaaqmIjWsJWXb8PZrZ1x2rjotbHvb4ja1r6Kth0fVerdu2+JR1EqgT4FhPYmK7LffmKZqytZHu2wb3n7tjMvOVaeFbW9b3Kb2VaT18Chb79ZtWzyKWjnlAgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJEOgAkASBDgBJ1Ap02zts/6vt52zPN1UUAKC8yoFu+xJJfybp5yVdLelW21c3VRgAoJw6R+jXSnouIr4ZEa9J+rSkXc2UBQAoq06gb5L0ra77JzuPAQDGwBFR7QvtX5L0cxHxwc792yRdGxF39YybkzTXubtV0kuSXqxccU4bxZp0Yz1WYk1WmqY1+dGImBk0qE773JOStnTd3yzp272DImJR0uLyfduHIqJdY950WJPXYz1WYk1WYk1WqnPK5V8kvcP2VbYvk7Rb0v5mygIAlFX5CD0iztv+dUn/KOkSSXsi4tnGKgMAlFLrPxZFxKOSHi35ZYuDh0wd1uT1WI+VWJOVWJMelf8oCgCYLFz6DwBJDC3QB7UFsP39th/qbD9oe3ZYtUyKAmtyh+0zto92Pj44jjpHxfYe26dtH+uz3bb/pLNeX7F9zahrHKUC63G97Ve6nh+/N+oaR8n2FttfsH3c9rO2715lzFQ9RwaKiMY/tPRH0n+T9HZJl0l6WtLVPWN+TdK9ndu7JT00jFom5aPgmtwh6ZPjrnWEa/JeSddIOtZn+02SPifJkq6TdHDcNY95Pa6X9Nlx1znC9bhS0jWd21dI+voqr5mpeo4M+hjWEXqRtgC7JD3Qub1P0nbbHlI9k4BWCT0i4glJL19kyC5JfxlLviRpg+0rR1Pd6BVYj6kSES9ExJHO7e9IOq6VV6NP1XNkkGEFepG2ABfGRMR5Sa9I+qEh1TMJirZK+MXOr477bG9ZZfs0ob3ESu+x/bTtz9n+8XEXMyqdU7LvlnSwZxPPkS7DCvTVjrR7305TZEwmRb7fv5c0GxHvkvRP+t5vMNNq2p4jgxzR0iXgPyHpTyX97ZjrGQnbb5L0GUkfjohXezev8iVT+xwZVqAXaQtwYYztdZJ+QLl/3Ry4JhHxUkT8T+fun0v6yRHVNqkKtZeYFhHxakT8d+f2o5Iutb1xzGUNle1LtRTmD0bEw6sM4TnSZViBXqQtwH5Jt3du3yLpQHT+ypHUwDXpOfd3s5bOGU6z/ZJ+ufNOhuskvRIRL4y7qHGx/bblvzPZvlZLr9+XxlvV8HS+1/skHY+Ij/cZxnOkS60rRfuJPm0BbP++pEMRsV9LP6i/sv2clo7Mdw+jlklRcE0+ZPtmSee1tCZ3jK3gEbC9V0vv3Nho+6SkeyRdKkkRca+WrkK+SdJzkr4r6QPjqXQ0CqzHLZJ+1fZ5Seck7U5+ELRN0m2SnrF9tPPY70hqSdP5HBmEK0UBIAmuFAWAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEiCQAeAJAh0AEji/wHLLq/GHhJD+AAAAABJRU5ErkJggg==\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": [],
   "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.206731390898</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.073656188013206e-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=\"$('#WZxQsddGsh').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.140799</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>Nbkg</td>\n",
       "        <td>251.519</td>\n",
       "        <td>16.7531</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>298.1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>19.4765</td>\n",
       "        <td>7.12533</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-20</td>\n",
       "        <td>271</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"WZxQsddGsh\" 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.140799 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & Nbkg & 251.519 & 16.7531 &  &  & 0.0 & 298.1 & No\\\\\n",
       "\\hline\n",
       "2 & Nsig & 19.4765 & 7.12533 &  &  & -20.0 & 271 & 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\", label=\"fit result\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYwAAAEPCAYAAABRHfM8AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xu8VXP+x/HXp/tFSUrSRTVCbpGjorlRUqLMCGEomWLkfs3lF8PMiAaD0ERR7pdBodCIIYSYUo1LpNGppiKJpOvn98daJ7vTPuesc1l7nb3P+/l4rMdel+9a67PatT+t9V3f79fcHRERkZJUSzoAERHJDkoYIiISiRKGiIhEooQhIiKRKGGIiEgkShgiIhKJEoaIiESihCEiIpEoYYiISCQ1kg6gIjVp0sTbtGmTdBgiIlnj/fff/8rdm0Ypm1MJo02bNsyaNSvpMEREsoaZ/TdqWT2SEhGRSJQwREQkEiUMERGJRAlDREQiUcIQEZFIYksYZtbKzF41s4/MbL6ZXRCub2xm08xsQfi5UxH7DwzLLDCzgXHFKSIi0cR5h7EJuMTdOwBdgWFmtg8wHHjF3dsDr4TL2zCzxsC1QBegM3BtUYlFREQyI7aE4e7L3P2DcP474COgBdAPmBAWmwAcl2b3o4Bp7r7K3b8BpgG94opVRERKlpGGe2bWBjgIeAdo5u7LIEgqZrZLml1aAItTlvPDdemOPRQYCtC6desyx9ht5HSWrF73UwCN6vLm8CPKfDwRkVwTe8Iwsx2AfwAXuvsaM4u0W5p1nq6gu48FxgLk5eWlLRPFktXrWDSyz9blNsNfKOuhRERyUqxvSZlZTYJk8bC7Px2uXm5mzcPtzYEVaXbNB1qlLLcElsYZq4iIFC/Ot6QMGAd85O63pmyaDBS89TQQmJRm95eAnma2U1jZ3TNcJyIiCYnzDqMbcBpwhJnNDqejgZHAkWa2ADgyXMbM8szsPgB3XwXcALwXTteH60REJCGx1WG4+wzS10UAdE9Tfhbw+5Tl8cD4eKITEZHSUktvERGJRAlDREQiUcIQEZFIlDBERCQSJQwREYlECUNERCJRwhARkUiUMEREJBIlDBERiUQJQ0REIlHCEBGRSJQwREQkEiUMERGJRAlDREQiUcIQEZFIlDBERCQSJQwREYkkthH3zGw8cAywwt33C9c9DuwVFmkErHb3A9Psuwj4DtgMbHL3vLji5Icf4NRTOfHH1kCf2E4jIpLtYksYwAPAaGBiwQp3P6lg3sxuAb4tZv/D3f2r2KIrUK8efPQRx6z/b+ynEhHJZrE9knL314FV6baZmQEnAo/Gdf5S6duXrl/OhW+Ly18iIlVbUnUYvwCWu/uCIrY78LKZvW9mQ4s7kJkNNbNZZjZr5cqVZYumb19qbdkEL71Utv1FRKqApBLGyRR/d9HN3TsBvYFhZvbLogq6+1h3z3P3vKZNm5YtmkMP5eu6DWHSpLLtLyJSBWQ8YZhZDeC3wONFlXH3peHnCuAZoHOsQVWvzvSfdYYpU2DjxlhPJSKSrZK4w+gBfOzu+ek2mll9M2tQMA/0BObFHdS09l1g9Wp44424TyUikpViSxhm9ijwNrCXmeWb2ZnhpgEUehxlZruZ2ZRwsRkww8zmAO8CL7j7i3HFWeCNNgdBnToweXLcpxIRyUqxvVbr7icXsX5QmnVLgaPD+YVAx7jiKsq6WnWgR4+gHuO22zJ9ehGRSk8tvVP17QuLFsG82J+AiYhkHSWMVMceG3zqbSkRke0oYaTadVfo0kX1GCIiaShhFNavH7z3Hs2+i79XEhGRbKKEUdhxxwHQc8HMhAMREalclDAK69AB9tmHoz95M+lIREQqFSWMdPr3p/Pi+bB8edKRiIhUGkoY6fTvT3XfAs8+m3QkIiKVhhJGOvvtx+eNW8BTTyUdiYhIpaGEkY4ZU/fqBq++Cl/pbSkREVDCKNKsvCNg82YuP+0Guo2cnnQ4IiKJU8IowgO3D4V27biZT1myel3S4YiIJE4Joyhm0L8//POfNPzx+6SjERFJnBJGcY4/HjZt4sgF7yQdiYhI4pQwinPIIdCqFUd/MiPpSEREEqeEURwzOOEEfvHFv2HVqqSjERFJlBJGSU45hVpbNsHTTycdiYhIouIconW8ma0ws3kp664zsyVmNjucji5i315m9omZfWZmw+OKMZJOnYJGfI88kmgYIiJJi/MO4wGgV5r1t7n7geE0pfBGM6sO3AX0BvYBTjazfWKMs3hmPNfhl/Daa7BkSWJhiIgkLbaE4e6vA2V58N8Z+MzdF7r7BuAxoF+FBldKkzv8Ctzh8ceTDENEJFFJ1GGca2Yfho+sdkqzvQWwOGU5P1yXlpkNNbNZZjZr5cqVFR0rAAt3bgkHHwyPPhrL8UVEskGmE8Y9wM+AA4FlwC1pyliadV7UAd19rLvnuXte06ZNKybKdE45BWbNgk8/je8cIiKVWEYThrsvd/fN7r4FuJfg8VNh+UCrlOWWwNJMxFesk04KXrPVXYaIVFEZTRhm1jxl8TfAvDTF3gPam1lbM6sFDAAmZyK+YrVoAb/6VZAwvMgbHhGRnBXna7WPAm8De5lZvpmdCdxsZnPN7EPgcOCisOxuZjYFwN03AecCLwEfAU+4+/y44iyVU06BTz6Bf/876UhERDKuxIRhZheYWUMLjDOzD8ysZ0n7ufvJ7t7c3Wu6e0t3H+fup7n7/u5+gLv3dfdlYdml7n50yr5T3H1Pd/+Zu/+5fJdYgY4/HmrWhIceSjoSEZGMi3KHMdjd1wA9gabAGcDIWKOqrBo3hmOPhYcfho0bk45GRCSjoiSMgreWjgbud/c5pH+TqWoYNAhWrIAXX0w6EhGRjIqSMN43s5cJEsZLZtYA2BJvWJVYr16wyy7wwANJRyIiklFREsaZwHDgEHf/AahF8FiqaqpZE373O3juOY33LSJVSo0IZaa5e/eCBXf/2syeALoXs09uGzgQbr016JDw/PO329xt5PRthnVt0agubw4/ItKhy7OviEicikwYZlYHqAc0CbvwKKi3aAjsloHYKq8DDoBOnYLHUmkSxpLV61g0ss/W5TbDX4h86PLsKyISp+IeSZ0FvA/sHX4WTJMIepOt2gYNCtpjzJmTdCQiIhlRZMJw99vdvS1wqbu3c/e24dTR3UdnMMbK6ZRTgvqMCROSjkREJCNKrMNw9zvN7DCgTWp5d58YY1yV3847Q9++QSO+m24KkoeISA6L0tL7QeCvwM+BQ8IpL+a4ssOgQbByZfDGlIhIjovyllQesI+7etzbTq9e0LIljB0Lv/1t0tGIiMQqSjuMecCucQeSlWrUgDPPhJdfhi++SDoaEZFYRUkYTYD/mNlLZja5YIo7sKzx+98H42Tce2/SkYiIxCrKI6nr4g4iq7VsCX36wPjx8Mc/qvJbRHJWiXcY7v4vYBFQM5x/D/gg5riyy9ChsHw5TNaNl4jkrhLvMMxsCDAUaEwwHncLYAxVqGuQFo3qFtviulXD2rzRqhX8/e/BmBkiIjkoyiOpYQRjb78D4O4LzGyXWKOqZErqy6nN8BeCuoxrr4XPP89QVCIimRWl0nu9u28oWDCzGkCJr9ia2XgzW2Fm81LWjTKzj83sQzN7xswaFbHvonAo19lmNivKhSRu8GCoVg3uuy/pSEREYhElYfzLzK4C6prZkcCTQJSWag8AvQqtmwbs5+4HAJ8CVxaz/+HufqC7Z0cjwZYt4ZhjYNw4am/aUHJ5EZEsEyVhDAdWAnMJOiScAlxT0k7u/jqwqtC6l919U7g4E2hZqmgru/POg5Ur6fPxG0lHIiJS4aIkjH7ARHc/wd37u/u9FdTqezAwtYhtDrxsZu+b2dDiDmJmQ81slpnNWrlyZQWEVQ7du0OHDgx6/zlQw3gRyTFREkZf4FMze9DM+oR1GOViZlcDm4CHiyjSzd07Ab2BYWb2y6KO5e5j3T3P3fOaNm1a3tDKxwzOP58D/vcZvP12srGIiFSwKO0wzgD2IKi7OAX43MzKXLNrZgOBY4BTi7pTcfel4ecK4BmCt7Syw2mnsaZ2fbjjjqQjERGpUFHuMHD3jQSPjx4jGESpX1lOZma9gCuAvuH44OnK1DezBgXzQE+C/qyyQ/36PH7AkfDUU5Cfn3Q0IiIVJkr35r3M7AHgM6A/cB/QPMJ+jwJvA3uZWb6ZnQmMBhoA08JXZseEZXczsynhrs2AGWY2B3gXeMHdXyz9pSVnYqdjYMsWGDMm6VBERCpMlPqIQQR3Fme5+/qoB3b3k9OsHldE2aXA0eH8QqBj1PNURosb7QrHHht0e35NiS+UFatwK/MWjeqW2JBQRCQOUUbcG2BmuwO/AP5pZnWBGu7+XezRZbPzzw/6lnrsMaDslfGFk0NxXZSIiMQpyiOpIcBTwN/DVS2BZ+MMKicccQTsvz/89a96xVZEckKUSu9hQDdgDQR9SQFVqi+pMjGDyy6D+fP59cLs6N1ERKQ4sfUlJcCAAdCyJWe9+3TSkYiIlFucfUlJzZpw0UUc+uVceO+9pKMRESmX2PqSktCQIUFDvlGjko5ERKRcorwltQW4N5yktBo04KGDenPOP/4Bn30Ge+yRdEQiImUSqaW3lM/9B/eFGjXg1luTDkVEpMyUMDJg5Q6N4bTT4P77YcWKpMMRESmTIhOGmT0Yfl6QuXBy2GWXwfr1291ldBs5nTbDX9g6tWhUN6EARUSKV1wdxsFhC+/BZjYRsNSN7r4q/W6S1l57Ba/Zjh4Nl14KTZoAsGT1OhaN7JNwcCIiJSvukdQY4EVgb4IealMntUQri6uvhrVr4W9/SzoSEZFSKzJhuPsd7t4BGO/u7dy9bcrULoMx5o5994X+/YOxMr75JuloRERKJcoASn8ws45mdm44HZCJwHLWNdfAd99pgCURyTpROh88n2Ao1V3C6WEzOy/uwHJWx47Qr1/wWOrbb5OORkQksiiv1f4e6OLuI9x9BNAVGBJvWDnu//4PVq8OKsBFRLJElIRhwOaU5c0UemNKSungg6FPH7jlFhr++H3S0YiIRBIlYdwPvGNm15nZdcBMihg5rzAzG29mK8xsXsq6xmY2zcwWhJ87FbHvwLDMAjMbGOV8WeWGG+Cbbxjy7jNJRyIiEkmUSu9bgTOAVcA3wBnuHvW90AeAXoXWDQdecff2wCvh8jbMrDFwLdAF6AxcW1RiyVoHHQQnnsjgWZNg+fKkoxERKVGkrkHc/YPwNdvb3f3fUQ/u7q8TJJpU/YAJ4fwE4Lg0ux4FTHP3Ve7+DTCN7RNP9rvhBmpv2gA33ph0JCIiJSqxt9oYNHP3ZQDuvszM0o3e1wJYnLKcH67bjpkNBYYCtG7duoJDjaZFo7rbjLUduXuPPffkyf17cPI998BFF8Huu8cUoYhI+SWRMKJIV6medpQ/dx8LjAXIy8tLZCTAN4cfUeZ97+h2Mid/8i+4/noYF6lqSEQkEcU+kjKz6mb2zwo+53Izax4evzmQrvvWfKBVynJLYGkFx1EpLGvYFM45Bx54AD7+OOlwRESKVGzCcPfNwA9mtmMFnnMyUPDW00BgUpoyLwE9zWynsLK7Z7guN115JdSvH3yKiFRSUSq9fwTmmtk4M7ujYIpycDN7FHgb2MvM8s3sTGAkcKSZLQCODJcxszwzuw+29oR7A/BeOF2f073jNm0aJItnn4XXXks6GhGRtKLUYbwQTqXm7icXsal7mrKzCFqVFyyPB8aX5bxZ6cILYcwYuPhimDULqmlsKxGpXKKM6T3BzOoCrd39kwzEVDXVrRu8XnvqqfDggzAw99oqikh2i9L54LHAbIKxMTCzA81sctyBVUkDBkDnznDVVcG4GSIilUiU5x7XEbS2Xg3g7rOBtjHGVHVVqxYM4bp0KdxyS9LRiIhsI0rC2OTuhfvhTqS9Q5XQrRuccALcdBPk5ycdjYjIVlEqveeZ2SlAdTNrD5wPvBVvWLmr28jpLFm9buty2lbhN98Mzz0Hl1wCjz+ewehERIoWJWGcB1wNrAceJWgPcUOcQeWyJavXsWhkn+ILtWkT1GOMGAFDhkCPHhmJTUSkOFF6q/3B3a8meBX2cHe/2t1/jD+0Ku6yy+BnP4Nzz4UNG5KORkQk0ltSh5jZXOBDggZ8c8zs4PhDq+Lq1IE774RPPoHbbks6GhGRSJXe44Bz3L2Nu7cBhhEMqiRx690bjjsu6Jhw8eKSy4uIxChKwvjO3d8oWHD3GcB38YUk27jtNtiyJWgJLiKSoCIrvc2sUzj7rpn9naDC24GTgNfiD02AoAJ8xIigEvyZZ4BaSUckIlVUcW9JFW45dm3KvNphZNKll8ITT8CwYTQ8UfUZIpKMIhOGux+eyUCkGDVrwn33QefOXPnqeIKbPBGRzCqxHYaZNQJOB9qklnf38+MLS7Zz8MFwySWcPGoUvPoqHK58LiKZFaXh3hRgJjAX2BJvOFKs664jf9zDbPrtqfQafCc/1qxTbPEWjeqWa/hYEZFUURJGHXe/OPZIpGT16tHyqYfgiCP4uNpMGDmq2OJthpdpGBMRkbSivFb7oJkNMbPmZta4YCrrCc1sLzObnTKtMbMLC5X5tZl9m1JmRFnPl3MOPxzOOivozfaNN0ouLyJSQaLcYWwARhH0J1XwdpQD7cpywnAQpgMBzKw6sAR4Jk3RN9z9mLKcI+f99a/wz3/C6afDnDnQsGHSEYlIFRDlDuNiYI+wpXfbcCpTskijO/C5u/+3go5XNeywA0ycCF9+CRddlHQ0IlJFREkY84EfYjr/AIIGgekcGvZbNdXM9o3p/NnrsMNg+HAYPx4mTUo6GhGpAqI8ktoMzDazVwm6OAfK/1qtmdUC+gJXptn8AbC7u39vZkcDzwLtizjOUGAoQOvWrcsTUva59lqYOjXoAv3QQ2GXXZKOSERyWJQ7jGeBPxMMmvR+ylRevYEP3H154Q3uvsbdvw/npwA1zaxJuoO4+1h3z3P3vKZNm1ZAWFmkVi148EFYswYGDwZXA3wRiU+JdxjuPiGmc59MEY+jzGxXYLm7u5l1JkhsX8cUR3bbd9+gEvy884LxwC+5JOmIRCRHRWnp/QVp+o4qT8W3mdUDjgTOSll3dnjcMUB/4A9mtglYBwxw13+fizRsWND6e/jwYEzwrl2TjkhEclCUOoy8lPk6wAlAmdthQDCKH7BzoXVjUuZHA6PLc44qxQzGjYODDoIBA+Df/4addko6KhHJMVEeSRV+FPQ3M5sBqDFdZdKoETz+OPz853DGGfDMM7RoVLfY1t7Z2HVIt5HTWbJ63dblbLwGkWwV5ZFUp5TFagR3HA1ii0jKrnNnuOkmuPhi+NvfeHN48W00srHrkCWr17FoZJ+ty9l4DSLZKsojqdRxMTYBi4ATY4lGyu/CC+H11+Gyy+DAA9WrrYhUmCiPpPSLk03MYMIE6NIFTjwRZs2C3XdPOioRyQFRHknVBo5n+/Ewro8vLCmXhg2D1t+HHAK/+Q3MmAH16iUdlYhkuSgN9yYB/QgeR61NmaQy23NPeOQRmD0bhg5Voz4RKbcodRgt3b1X7JFIxevTB264Aa65JqjPuPTSpCMSkSwW5Q7jLTPbP/ZIJB5XXQUnnACXXw7PpOtFXkQkmigJ4+fA+2b2iZl9aGZzzezDuAOTCpJaCX7qqfDuu0lHJCJZKsojqd6xRyHxqls3qATv2hWOPRZmzoS2bZOOSkSyTJTXajW4UTmltrhu0ahuMkHssgtMmRJ0g96nD7z1Vql2j7OFdXmOXbg1e2n2VatxkdKJcoch5VRpfoT23juox+jZE/r2pfYhF5a8TyjOFtblOXbhP9vS7KtW4yKlE6UOQ3LJr38djKExYwZ3TRoJGzcmHZGIZAkljKropJPg7rvp8fl7QUeFW7YkHZGIZAEljKrq7LO5+Zenw8MPwwUXqGGfiJRICaMKu7vrCUFjvtGj4eqrlTREpFiq9K7KzODmm4MxwW+8MVj+05+CTxGRQhJLGGa2CPgO2Axscve8QtsNuB04GvgBGOTuH2Q6zpxnBvfcE8z/5S/BXcaf/6ykISLbSfoO43B3/6qIbb2B9uHUBbgn/JSKVq1akDTMgjsNUNIQke0knTCK0w+Y6O4OzDSzRmbW3N2XJR1YTqpWDe6+O5i/8UbYvBlGjlTSEJGtkkwYDrxsZg783d3HFtreAlicspwfrtsmYZjZUGAoQOvWreOLtiooSBrVqgV1G998E9x5VK+e0TBK1Xp740Z44AF48kn4/HOeXVcdar8L558PO++cmYBFqogkE0Y3d19qZrsA08zsY3d/PWV7uv/abvcaT5hoxgLk5eXpNZ/yqlYN7roLGjcOHkt98w089FBGQ4jcenv27KBDxf/8Bzp0gC5dWD9jbtCl+113BeOB9OyZgYhFqobEXqt196Xh5wrgGaBzoSL5QKuU5ZbA0sxEV8UVvC11663w1FNwzDHU27Cu5P0yadIk6NYNvv02mJ8/Hx55hJNOGQlz5sBuuwV9Zk2alHSkIjkjkYRhZvXNrEHBPNATmFeo2GTgdAt0Bb5V/UWGXXRR0DX6q6/yyGNXwfLlSUcUePll6N8f9tsvGLO8b99t61r23x/eeAMOOigY11xduotUiKTuMJoBM8xsDvAu8IK7v2hmZ5vZ2WGZKcBC4DPgXuCcZEKt4k4/HZ55hj2/+jIYU2Ne4byeWR1WLAySxT77wLRpsOuu6QvuuCNMnQrNmwflv/46s4GK5KBEEoa7L3T3juG0r7v/OVw/xt3HhPPu7sPc/Wfuvr+7z0oiVgGOPZYTT7kJNmyAww6DF19MJo7Fi7n/yeuCZDBlCjRsWHz5nXcOHqktWxbcLYlIuahrEIlk3q57BI922rWDPn047YPnM9uVyOrV0Ls39Tb8GCSLFi2i7ZeXB8OHBz30vvRSvDGK5DglDImuZUuYMQP69OGGaWNg8GBYF39leK1NG+G3v4VPP+Ws314d1FGUxtVXB2OBnHMOrF8fT5AiVYAShpTODjvAM89w+2EnB+0fDjsMFi6M73zu3DT1dnj1Vbj/ft7evWPpj1GnDtx+exBnQeNEESk1JQwpverVue0Xp8ILL8CiRXDwwTB5cjznuuIKfvOf14J+rk49tezH6dkTjjoqaKOxalWFhSdSlShhSNkdfTS8/z60bQv9+sHQofD99xV3/JtvhlGjmNCpT1APUV6jRgXtNv7yl/IfS6QKqsx9SUmCuo2czpLVP9VPtGhUN33Bdu3g7bd58KhBnHrvfSx+8jluPOlKxow5v9hjlTjO+d13wxVXwIABXNf6FAaWsU+rwue+p2N3et91F1xySZmOF5fUOCP9+YgkQAlD0lqyeh2LRvaJVrh2bf6v6+847YZz2P3007n77xfClnnB/+SbNNnuWEV29QHBm1d/+hOMGAHHHAMTJuAjplXYdfzqm6X0/nB6cPdSu0eZj1vRUuMs9s9HJEF6JCUV5xe/gDlzGHdIPxg/HvbcE0aNit6tyKpVwdtQI0bA734HTz8NtWpVaIj/3Wm3oDHimDE0/V51GSKloYQhFathQ/58xO+D/pwOOQQuv5w37xkMV14ZdBKYRoP1a+G226B9e3j++aAPq4kToWbNeGK85hrYtIlzZj4Zz/FFcpQeSUk89t03aCg3cyazTruAI0eNCsbXaN2aMXVawJoXgjYRCxfy3oy3YNMG6N4dbrkFOpbh1dnSaNcOBg7klAcmwpIl0RsBilRxusOQeHXtypDjRwQ/zKNHQ5cu7PHVYnjiiSChrF3Lwwf2DlqRT5sWf7IocM01VPMtP40wKCIlUsKQzGjWDIYNgyeeoMeQMfDVV5CfD+++yw3dhwSPrzI5ul+bNjy5/5Fw773w5ZeZO69IFlPCkCrrzsNOCmb+9KdkAxHJEkoYUmUta9gUzjoreKPr88+TDkek0lPCkKrtyiuDt7Guvz7pSEQqPb0lJRnXolHdbRqnFdmKPE35ksqWWvPmcO65cOut/GzwocUWLdxqvHCMap0tuU4JQzKutD+ssf8QX3453HMPF854BDi7yGLFtX5X62ypCjL+SMrMWpnZq2b2kZnNN7ML0pT5tZl9a2azw2lEpuOUKqRpU7jwQo79+A2N/y1SjCTqMDYBl7h7B6ArMMzM9klT7g13PzCc9IBZ4nXFFayovxOcfz5s2ZJ0NCKVUsYThrsvc/cPwvnvgI8ANbWVZDVowE2/GgTvvAMPPZR0NCKVUqJvSZlZG+Ag4J00mw81szlmNtXM9i3mGEPNbJaZzVq5cmVMkUpV8PR+h0PnzkG36mvWJB2OSKWTWMIwsx2AfwAXunvhf50fALu7e0fgTuDZoo7j7mPdPc/d85o2bRpfwJLz3KrBnXfCihVw6aVJhyNS6SSSMMysJkGyeNjdny683d3XuPv34fwUoKaZNclwmFIVde4cDK50773w8stJRyNSqSTxlpQB44CP3P3WIsrsGpbDzDoTxPl15qKUKu3662HvveHMM4M+r0QESOYOoxtwGnBEymuzR5vZ2WZW8BJ8f2Cemc0B7gAGuLsnEKtURXXqBBXfK1fCySfD5s1JRyRSKWS84Z67zwCK7ZbU3UcDozMTUdVVuMV14W2l2bfCW2CXUUlxFddqvPC+Q485j6v+8Ve46iq46aZSxVGmccxjUFzrdFALdSkdtfSuwsrzQ1FZf2RKiqu47YW3tRkOV/1hbTD+d9OmQIfIcZRqHPMYlTQ2u1qoS2koYYgU58474euv4bLLGNR9KFD0j69IrlPCEClO9erw4IOwYQPXPTsWLtkhuOOoXj3pyEQyTt2bi5SkVi146inuP/hYuPVWOPJIjdInVZIShkgU1avzxx7hYEvvvQf77x8kj/Xrk45MJGOUMERK44wz4MMP4bDDggZ+HTrAXXdRd8OPSUcmEjslDJHSatsWpk6Fl16CXXaBc89l5t0DYcgQmD5d7TYkZylhiJRVz54wcya89Rav7NEZHnsMuneH5s25ffIouP9+yM9POkqRCqOEIVJehx7KxcdcEnRa+OST0LMnh305BwYPhlatoE0bRk+6CW67Dd5+G37U4yvJTnqtVqSi1K0L/ftD//4ccsXzLDq1NbzyCszed4DHAAAK70lEQVScyYEvvgoXXxyUq1EjqPvo2BEOOAA6dqTJ2m+SjV0kAiUMqXCVtduQOG13zTvVC5LBAQcAcNLI6WzIX8JBSz+h47JP6bDiCzo8O5Xm4WBNswCeuAQ6duTaVXXgni+DDhD33ht23RWs2N50tpHaHUhZunhJbfFe+FgV1cK/tF2nVJauVqo6JQypcFXxH3KZuyT5+uvgras5c4Jp7lxO/HA+nPPcT2UaNvwpeey9N7RvD+3aBVOjRtsdsqTuQIqLq3BXIanHqshuRErbdUpl6WqlqlPCEEnSzjvD4YcHU2jfK55n0XkHwscfbztNnw4TJ267/047BYmjbdutSeTnX6yEz/aCli2DnndFKogShkhlYxb82LdsCT16bLvtu+/g88/hiy9g4cKfpg8/hMmTYcMGHgJ44v+C8k2aBBXvBcdr2XLb5RYtoF69TF+hZCklDJFs0qABHHhgMBW2ZQssXcqJVzzMEz13hcWLg9d68/ODrkzefBNWrdp+vx13hGbNoFkz7lqxBb6fGtSbNGtGjwWL4Z0m0KwZtTeqVXtVp4QhkiuqVYOWLXm31X4wsIg6jB9+gCVLfkokixfD//4XTMuXs/fK/8LD82H1agDuA3j6BgA+ARi7AzRuHDxKK+qz8LpGjaBmzUz8CUjMlDBEqpJ69YJK8/bt027uPvyFoHJ5/XpYsYJjr/kHz53YHv73P26e+DqXd2oc3KWsWhVU2OfnB5+rVgV3OEWpWzeovN9xR2jYkIdWboRPx21dvmjWCvjrR1uXt07160P9+jT77qsgidWvr+SToEQShpn1Am4HqgP3ufvIQttrAxOBgwnG8j7J3RdlOk6RKqt2bWjVirnN20Of4G7l7gW7cnlRb19t2QJr1vyUSAo+v/4avv02mNas2fpZe9l/YcGCrevO+3YNvPVYkeG8A3B3uFCzJnOq1YKHGm1NKNSvHyTD1OX69YNEVadOMNWu/dN8uqmo7TVrluq15lyW8YRhZtWBu4AjgXzgPTOb7O7/SSl2JvCNu+9hZgOAm4CTMh2riERUrVrw6KlRo+BtrRKcUHAnE2p3xfMsuuZX2yQV1qwJHqGtXcuVD87kxqPawdq1sHYtT0+bzxkdm25dZu3a4A5kyZJt1/34Y/F3PlGYBcmkVq1gqlkzmFLnCy+XZVuNGsE4KzVqlH6+bl3o0qV81xlBEncYnYHP3H0hgJk9BvQDUhNGP+C6cP4pYLSZmbt7JgMVkQwxCyr0GzQI3twq5NG5jbjxop8SzB/9Bc6I2NaETZuCxFEwrV+/7XJR61LXr1sHGzfChg3BZ3Hz69YFyS5q+YrQrFlQDxUzy/RvsJn1B3q5++/D5dOALu5+bkqZeWGZ/HD587DMV2mONxQYGi7uRVg3F2oCbLdPlsu1a8q164Hcu6Zcux7IvWsqz/Xs7u5NoxRM4g4j3cPAwlkrSplgpftYYGzaE5nNcve80oVXueXaNeXa9UDuXVOuXQ/k3jVl6nqS6K02H2iVstwSWFpUGTOrAewIpHmBXEREMiWJhPEe0N7M2ppZLWAAMLlQmcnAwHC+PzBd9RciIsnK+CMpd99kZucCLxG8Vjve3eeb2fXALHefDIwDHjSzzwjuLAaU8XRpH1VluVy7ply7Hsi9a8q164Hcu6aMXE/GK71FRCQ7acQ9ERGJRAlDREQiyYmEYWa9zOwTM/vMzIan2V7bzB4Pt79jZm0yH2V0Ea5nkJmtNLPZ4fT7JOKMyszGm9mKsH1Nuu1mZneE1/uhmXXKdIylFeGafm1m36Z8RyMyHWNpmFkrM3vVzD4ys/lmdkGaMlnzPUW8nmz7juqY2btmNie8pj+mKRPvb527Z/VEUHH+OdAOqAXMAfYpVOYcYEw4PwB4POm4y3k9g4DRScdaimv6JdAJmFfE9qOBqQTtb7oC7yQdcwVc06+B55OOsxTX0xzoFM43AD5N8/cua76niNeTbd+RATuE8zUJutjqWqhMrL91uXCHsbWrEXffABR0NZKqHzAhnH8K6G5WaXsTi3I9WcXdX6f4djT9gIkemAk0MrPmmYmubCJcU1Zx92Xu/kE4/x3wEVC4j46s+Z4iXk9WCf/cvw8Xa4ZT4beWYv2ty4WE0QJYnLKcz/Z/MbaWcfdNwLfAzhmJrvSiXA/A8eFjgafMrFWa7dkk6jVnm0PDxwdTzWzfpIOJKnyMcRBhJ7EpsvJ7KuZ6IMu+IzOrbmazgRXANHcv8juK47cuFxJGhXY1UglEifU5oI27HwD8k5/+R5Gtsun7ieoDgj56OgJ3As8mHE8kZrYD8A/gQndfU3hzml0q9fdUwvVk3Xfk7pvd/UCCHjI6m9l+hYrE+h3lQsLIta5GSrwed//a3QvGy7yXYNyQbBblO8wq7r6m4PGBu08BappZk4TDKpaZ1ST4cX3Y3Z9OUySrvqeSricbv6MC7r4aeA3oVWhTrL91uZAwcq2rkRKvp9Bz474Ez2ez2WTg9PAtnK7At+6+LOmgysPMdi14dmxmnQn+rX2dbFRFC2MdB3zk7rcWUSxrvqco15OF31FTM2sUztcFegAfFyoW629d1g/R6pntaiR2Ea/nfDPrC2wiuJ5BiQUcgZk9SvBGShMzyweuJaiww93HAFMI3sD5DPgBOCOZSKOLcE39gT+Y2SZgHTCgEv8nBaAbcBowN3xGDnAV0Bqy8nuKcj3Z9h01ByZYMAhdNeAJd38+k7916hpEREQiyYVHUiIikgFKGCIiEokShoiIRKKEISIikShhiIhIJEoYIiISiRKGiIhEooQhkmXM7Dgzu9fMJplZz6TjkapDCUMkIjNrY2brUloOY2bNzOwRM1toZu+b2dtm9psSjvOamR1VaN2FZnZ3yvLfzaxbuv3d/Vl3H0LQwv+ksHzdcBCgDdnSH5JkHyUMkdL5POwttKC/omeB1929nbsfTNAVQ8sSjvEo23fZMCBcX6ALMLOE41wD3AXg7uvCuCptZ4CS/ZQwpMoK7xg+NrP7zGyemT1sZj3M7E0zWxB2SFecI4ANYb9EALj7f939zpRz/C4cVnN2eNdQnWBgm2PMrHZBHMBuwIxwuQPwqbtvNrPTw3FP5pjZg+F2M7ObgKkFgwSJZIIShlR1ewC3AwcAewOnAD8HLiXorK44+xKMqZBW+MN/EtAt/N//ZuBUd/8aeJefuqYuGEqzoGO33sCL4YA+VwNHhGM2FIxLfR5BT6X9zezsUlyrSLlkfW+1IuX0hbvPBTCz+cAr7u5mNhdoU5oDmdldBMlmg7sfAnQnGKvkvbAX7boEI6XBT4+lJoWfg1MOdRRBT7DHA0+5+1cA7r4q/LwDuKPUVypSTkoYUtWtT5nfkrK8hZL/fcwn+FEHwN2HhRXOs8JVBkxw9yvT7PsscKuZdQLqFjxaMrN6QCN3XxrWkag7aak09EhKpOymA3XM7A8p6+qlzL9C8NhoFwAza2xmuwOEI729Boxn28ruw4FXU/Y/0cx2Ltg/josQiUoJQ6SMwjqH44BfmdkXZvYuwfjqV4Tb/0PwJtPLZvYhMI1gEJwCjwIdgcdS1vUGXgz3nw/8GfiXmc0BihoJTyQjNICSSETh20zPu/t+MZ7jA6CLu28s4/6LgLyCeg+RiqQ7DJHoNgM7pjbcq2ju3qksyaKg4R7BMLFbKj4yEd1hiIhIRLrDEBGRSJQwREQkEiUMERGJRAlDREQiUcIQEZFIlDBERCQSJQwREYlECUNERCL5f9GHmx4nTUDYAAAAAElFTkSuQmCC\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)\n",
    "plt.xlabel(\"m [GeV/c$^2$]\")\n",
    "plt.ylabel(\"number of events\")\n",
    "plt.savefig(\"fit_discovery_ex.png\")"
   ]
  },
  {
   "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.987538447270647e-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=\"$('#GScZfFUSkB').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.140798</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>Nbkg</td>\n",
       "        <td>251.519</td>\n",
       "        <td>16.7524</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>298.1</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>19.4765</td>\n",
       "        <td>7.12552</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>-20</td>\n",
       "        <td>271</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"GScZfFUSkB\" 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.140798 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & Nbkg & 251.519 & 16.7524 &  &  & 0.0 & 298.1 & No\\\\\n",
       "\\hline\n",
       "2 & Nsig & 19.4765 & 7.12552 &  &  & -20.0 & 271 & 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 0x1a3afb1908>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config.bestfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = AsymptoticCalculator(config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "poinull = POI(Nsig, value=0)\n",
    "discovery_test = Discovery(poinull, calc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute qobs for the null hypothesis!\n",
      "\n",
      "p_value for the Null hypothesis = 0.0007571045424956679\n",
      "Significance = 3.1719464825102244\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'pnull': 0.0007571045424956679, 'significance': 3.1719464825102244}"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "r = discovery_test.result();\n",
    "r"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}