examples/notebooks/confidence_interval_zfit_asy.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 ConfidenceInterval\n",
    "from lauztat.calculators import AsymptoticCalculator\n",
    "from lauztat.config import Config"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import zfit\n",
    "from zfit import ztf\n",
    "from zfit.core.loss import ExtendedUnbinnedNLL, UnbinnedNLL\n",
    "from zfit.minimizers.minimizer_minuit import MinuitMinimizer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAADnZJREFUeJzt3W+MXXWdx/HPZ6WaRsgWbcFm7Oy4hhgxUcAJYdPEIMRdxMRqVhN4gEUhY1QUEp9MeKDGfdIHismuRjOExu4GWTeC2qX1DwJJI9HGKalQMipoqrY0tEj4Y2w0NV8fzClcL3N7zz1/7p/vfb+Sydx77rn3fH/3dD65/c33/MYRIQDA5PuHURcAAGgGgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJDEWcM82MaNG2Nubm6YhwSAiXfgwIGnI2JTv/2GGuhzc3NaXl4e5iEBYOLZ/m2Z/ZhyAYAkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkCHQASIJAB4AkhnqlaFO27nhAR589KUma2bBeDy1eMeKKAGD0JjLQjz57Uod3vEeSNLe4Z8TVAMB4YMoFAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgCQIdAJIg0AEgib6BbnuL7Qdtr9h+zPbNxfbP2T5q+2DxdXX75QIAeinzF4tOSfp0RDxs+xxJB2zfVzz2pYj4QnvlAQDK6hvoEXFM0rHi9gu2VyTNtF0YAGAwA82h256TdLGk/cWmm2w/Ynun7XMbrg0AMIDSgW77bEl3S7olIp6X9FVJb5R0kVY/wX+xx/MWbC/bXj5x4kQDJf+9mQ3rNbe4R3OLe7R1xwONvz4ATIoyc+iyvU6rYX5nRNwjSRHxVMfjt0u6d63nRsSSpCVJmp+fj7oFd3to8YoXb88t7mn65QFgYpTpcrGkOyStRMRtHds3d+z2fkmHmi8PAFBWmU/oWyVdJ+lR2weLbbdKutb2RZJC0mFJH22lQgBAKWW6XH4syWs8tLf5cgAAVXGlKAAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBIEOgAkQaADQBJ9A932FtsP2l6x/Zjtm4vtr7F9n+3Hi+/ntl8uAKCXMp/QT0n6dES8WdJlkj5h+0JJi5Luj4gLJN1f3AcAjEjfQI+IYxHxcHH7BUkrkmYkbZO0q9htl6T3tVUkAKC/gebQbc9JuljSfknnR8QxaTX0JZ3XdHEAgPJKB7rtsyXdLemWiHh+gOct2F62vXzixIkqNQIASigV6LbXaTXM74yIe4rNT9neXDy+WdLxtZ4bEUsRMR8R85s2bWqiZgDAGsp0uVjSHZJWIuK2jod2S9pe3N4u6bvNlwcAKOusEvtslXSdpEdtHyy23Spph6T/s32DpN9J+mA7JQIAyugb6BHxY0nu8fCVzZYDAKiKK0UBIAkCHQCSINABIAkCHQCSINABIAkCHQCSINABIAkCHQCSINABIAkCHQCSKLOWy1jYuuMBHX32pCRpZsP6EVfzct31PbR4xYgrAjBtJibQjz57Uod3vGfUZfTUWd/c4p4RVwNgGjHlAgBJEOgAkASBDgBJEOgAkASBDgBJTEyXS1PO1F5I6yGASTZ1gX6m9kJaDwFMMqZcACAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASAJAh0AkiDQASCJvoFue6ft47YPdWz7nO2jtg8WX1e3WyYAoJ8yn9C/LumqNbZ/KSIuKr72NlsWAGBQfQM9IvZJemYItQAAaqgzh36T7UeKKZlzG6sIAFBJ1dUWvyrpPyRF8f2Lkj6y1o62FyQtSNLs7GzFw9XTvSzuaTMb1v/dqoqdjwHApKkU6BHx1Onbtm+XdO8Z9l2StCRJ8/PzUeV4dXUui9uJ9c4BZFJpysX25o6775d0qNe+AIDh6PsJ3fZdki6XtNH2EUmflXS57Yu0OuVyWNJHW6wRAFBC30CPiGvX2HxHC7UAAGrgSlEASIJAB4AkUv2R6M42xJkN6+liATBVUgV6Z4B39pcDwDRgygUAkiDQASAJAh0AkiDQASAJAh0AkiDQASCJVG2LwFq6l0/m+gRkRaAjvc7lk7k+AZkx5QIASRDoAJAEgQ4ASRDoAJAEgQ4ASdDlUkN3O1yZfWiZq4b3EeiPQK+hsx2uzD60zFXH+wj0x5QLACRBoANAEgQ6ACRBoANAEgQ6ACSRtstlZsP6F7sherUUjiPa8wbD+wW8JG2gT+oPNu15g+H9Al7ClAsAJEGgA0ASfQPd9k7bx20f6tj2Gtv32X68+H5uu2UCAPop8wn965Ku6tq2KOn+iLhA0v3FfQDACPUN9IjYJ+mZrs3bJO0qbu+S9L6G6wIADKjqHPr5EXFMkorv5zVXEgCgitbbFm0vSFqQpNnZ2bYP1wp6ncdLr2sMurdznjBtqgb6U7Y3R8Qx25slHe+1Y0QsSVqSpPn5+ah4vJGi13m89Arqzu2cJ0yjqlMuuyVtL25vl/TdZsoBAFRVpm3xLkk/kfQm20ds3yBph6R32X5c0ruK+wCAEeo75RIR1/Z46MqGawEA1MCVogCQBIEOAEmkXW1x3A3aCknrJIB+CPQRGbQVktZJAP0w5QIASRDoAJAEgQ4ASRDoAJAEgQ4ASRDoAJAEbYsD6rV0a53n1nnNbOi3B6oj0AdUJ2DKLPs67ei3B6pjygUAkiDQASAJAh0AkiDQASAJAh0AkqDLpYdxbiXM8Nfty7QnjvM5AMYRgd7DOIdkhr9uX6Y9cZzPATCOmHIBgCQIdABIgkAHgCQIdABIgkAHgCTochljTbXt9WoRbGtlQ1ZMBEaDQB9jTQVhrxbBtlY2ZMVEYDSYcgGAJAh0AEii1pSL7cOSXpD0V0mnImK+iaIAAINrYg79nRHxdAOvAwCogSkXAEiibqCHpB/aPmB7oYmCAADV1J1y2RoRT9o+T9J9tn8REfs6dyiCfkGSZmdnax5uMkzLsq+DLoHbvR1As2oFekQ8WXw/bvvbki6VtK9rnyVJS5I0Pz8fdY43KablQhqWwAXGS+UpF9uvtn3O6duS/lXSoaYKAwAMps4n9PMlfdv26df5RkR8v5GqAAADqxzoEfEbSW9rsBYAQA20LQJAEgQ6ACTBaotJdbcUnlalpXISl8MtM84y4+rcZ61jNPFeTOL7i/FEoCfV2VLYqUpYTOJyuGXGWWZcvd7HMz1nUJP4/mI8MeUCAEkQ6ACQBIEOAEkQ6ACQBIEOAEkQ6ACQBG2LE66737pqD3OvZW5PP9amYS433OtYTb2PTddHTzoGQaBPuM4f+Do9zKMMjmEeu9exmnof6xqXOjCZmHIBgCQIdABIgkAHgCQIdABIgkAHgCTocklkGO1/w2wxHJWyYzxTq2eZY1Tt7mljud3MS/hmHls3Aj2RYfxDzfzDcFrZMdZ5L+q0JLax3G7mJXwzj60bUy4AkASBDgBJEOgAkASBDgBJEOgAkARdLmjENLQzDluZlRfL7NPZtnemY631mt37lHn9Mq2BdVoJez23jdeso8r7UheBjkZMQzvjsJVZebHMPp1te4Met1PZ1y/TGlinlbDXc9t4zTqqvC91MeUCAEkQ6ACQRK1At32V7V/afsL2YlNFAQAGVznQbb9C0lckvVvShZKutX1hU4UBAAZT5xP6pZKeiIjfRMRfJP2vpG3NlAUAGFSdQJ+R9PuO+0eKbQCAEXBEVHui/UFJ/xYRNxb3r5N0aUR8smu/BUkLxd03SfplcXujpKcrHXxyMebpwJinwzDH/E8RsanfTnX60I9I2tJx//WSnuzeKSKWJC11b7e9HBHzNY4/cRjzdGDM02Ecx1xnyuVnki6w/Qbbr5R0jaTdzZQFABhU5U/oEXHK9k2SfiDpFZJ2RsRjjVUGABhIrUv/I2KvpL0Vn/6yaZgpwJinA2OeDmM35sq/FAUAjBcu/QeAJFoP9H7LA9h+le1vFo/vtz3Xdk1tKzHm622fsH2w+LpxFHU2xfZO28dtH+rxuG3/Z/F+PGL7kmHX2LQSY77c9nMd5/gzw66xaba32H7Q9ortx2zfvMY+qc51yTGPz7mOiNa+tPrL0l9L+mdJr5T0c0kXdu3zcUlfK25fI+mbbdbU9lfJMV8v6cujrrXBMb9D0iWSDvV4/GpJ35NkSZdJ2j/qmocw5ssl3TvqOhse82ZJlxS3z5H0qzX+bac61yXHPDbnuu1P6GWWB9gmaVdx+1uSrrTtlutq09QtiRAR+yQ9c4Zdtkn671j1U0kbbG8eTnXtKDHmdCLiWEQ8XNx+QdKKXn51eKpzXXLMY6PtQC+zPMCL+0TEKUnPSXpty3W1qeySCP9e/Jf0W7a3rPF4JtO6TMS/2P657e/Zfsuoi2lSMTV6saT9XQ+lPddnGLM0Jue67UBf65N2d1tNmX0mSZnx/L+kuYh4q6Qf6aX/oWSV7RyX8bBWL9d+m6T/kvSdEdfTGNtnS7pb0i0R8Xz3w2s8ZeLPdZ8xj825bjvQyywP8OI+ts+S9I+a7P/K9h1zRPwhIv5c3L1d0tuHVNuolFomIpOIeD4i/ljc3itpne2NIy6rNtvrtBpsd0bEPWvsku5c9xvzOJ3rtgO9zPIAuyVtL25/QNIDUfymYUL1HXPXnOJ7tTovl9luSR8qOiAuk/RcRBwbdVFtsv26078Lsn2pVn/W/jDaquopxnOHpJWIuK3HbqnOdZkxj9O5bvWPREeP5QFsf17SckTs1uqb9T+2n9DqJ/Nr2qypbSXH/Cnb75V0Sqtjvn5kBTfA9l1a/U3/RttHJH1W0jpJioivafVq4qslPSHpT5I+PJpKm1NizB+Q9DHbpySdlHTNhH9QkaStkq6T9Kjtg8W2WyXNSmnPdZkxj8255kpRAEiCK0UBIAkCHQCSINABIAkCHQCSINABIAkCHQCSINABIAkCHQCS+Bt3yeaUArX8mAAAAABJRU5ErkJggg==\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, 500)\n",
    "peak = np.random.normal(1.2, 0.1, 80)\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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean = zfit.Parameter(\"mean\", 1.2, 0.5, 2.0)\n",
    "sigma = zfit.Parameter(\"sigma\", 0.1, 0.02, 0.2)\n",
    "lambda_ = zfit.Parameter(\"lambda\",-2.0, -4.0, -1.0)\n",
    "Nsig = zfit.Parameter(\"Nsig\", 20., 0., len(data))\n",
    "Nbkg = zfit.Parameter(\"Nbkg\", len(data), 0., len(data)*1.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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": null,
   "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": null,
   "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 = -2349.8670049241327</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 143</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 143</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 1.004524152061599e-05</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=\"$('#heYmoEIjeV').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.94354</td>\n",
       "        <td>0.119214</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>mean</td>\n",
       "        <td>1.20665</td>\n",
       "        <td>0.0178914</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.5</td>\n",
       "        <td>2</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>sigma</td>\n",
       "        <td>0.0995561</td>\n",
       "        <td>0.0208026</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.02</td>\n",
       "        <td>0.2</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>3</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>73.9469</td>\n",
       "        <td>13.173</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>490</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>4</td>\n",
       "        <td>Nbkg</td>\n",
       "        <td>416.003</td>\n",
       "        <td>22.6755</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>539</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"heYmoEIjeV\" style=\"display:none;\">\n",
       "<textarea rows=\"16\" 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.94354 & 0.119214 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & mean & 1.20665 & 0.0178914 &  &  & 0.5 & 2 & No\\\\\n",
       "\\hline\n",
       "2 & $\\sigma$ & 0.0995561 & 0.0208026 &  &  & 0.019999999552965164 & 0.2 & No\\\\\n",
       "\\hline\n",
       "3 & Nsig & 73.9469 & 13.173 &  &  & 0.0 & 490 & No\\\\\n",
       "\\hline\n",
       "4 & Nbkg & 416.003 & 22.6755 &  &  & 0.0 & 539 & 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": null,
   "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": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEPCAYAAABP1MOPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3XecFdX9//HXh2WBxQbKCgTFtSD2ioafWAiyBiWKBXtBLGiMsUUNxm7kG9RE7AUrUWIBjBLBgohdQUBAQI2KGEWiaERRmsDn98eZhWXdZWeXnTu3vJ+PxzzunblTPoer97NzzplzzN0REZHC1SjtAEREJF1KBCIiBU6JQESkwCkRiIgUOCUCEZECp0QgIlLglAhERAqcEoGISIFTIhARKXCN0w4gjlatWnlZWVnaYYiI5JRJkyZ97e6lte2XE4mgrKyMiRMnph2GiEhOMbNP4+ynqiERkQKnRCAiUuCUCERECpwSgYhIgVMiEBEpcEoEIiIFTolARKTAKRGIiBS4nHigrCF1Gfgic+YvWrnerkUJr/fvlmJEIiLpKrhEMGf+ImYP7Llyvaz/qBSjERFJn6qGREQKXGKJwMyamdkEM5tqZjPM7Opo+4Nm9omZTYmWXZKKQUREapdk1dASoJu7/2BmxcBrZvZM9NlF7j48wWuLiEhMiSUCd3fgh2i1OFo8qeuJiEj9JNpGYGZFZjYF+AoY4+7jo48GmNk0MxtkZk1rOLafmU00s4nz5s1LMkwRkYKWaCJw9+XuvguwCbCnme0AXAJsA+wBbAj8sYZjB7t7J3fvVFpa67wKIiJSTxnpNeTu84GXgB7uPteDJcADwJ6ZiEFERKqXZK+hUjNrEb0vAboD75tZ22ibAYcC05OKQUREapdkr6G2wBAzKyIknMfd/Wkze9HMSgEDpgBnJhiDiIjUIsleQ9OAXavZrvEcRESyiJ4sFhEpcEoEIiIFTolARKTAKRGIiBQ4JQIRkQKnRCAiUuCUCERECpwSgYhIgVMiEBEpcEoEIiIFTolARKTAKRGIiBQ4JQIRkQKnRCAiUuCUCERECpwSgYhIgVMiEBEpcEoEIiIFTolARKTAKRGIiBS4xBKBmTUzswlmNtXMZpjZ1dH2zc1svJl9aGaPmVmTpGIQEZHaJXlHsATo5u47A7sAPcysM3AdMMjdOwDfAqcmGIOIiNQisUTgwQ/RanG0ONANGB5tHwIcmlQMIiJSu8ZJntzMioBJwFbA7cDHwHx3Xxbt8jnQroZj+wH9ANq3b1/vGLoMfJE58xetXG/XoqTe5xIRyUeJJgJ3Xw7sYmYtgH8C21a3Ww3HDgYGA3Tq1KnafeKYM38Rswf2rO/hIiJ5LyO9htx9PvAS0BloYWYVCWgT4IvELjxoEPcOvzqx04uI5IMkew2VRncCmFkJ0B14DxgH9I526wM8lVQMLF9O94/fhk8+SewSIiK5Lsk7grbAODObBrwNjHH3p4E/AheY2UfARsB9iUVwxBHhdcSIxC4hIpLrEmsjcPdpwK7VbJ8F7JnUdVez+eZMbdOBnYcNgwsvzMglRURyTd4/WfxMxy4wYQJ8+mnaoYiIZKW8TwSjO3YJb1Q9JCJSrbxPBP9p2RZ23RWGDUs7FBGRrJT3iQCA3r3hrbfgs8/SjkREJOsUTiIAVQ+JiFSjMBLB1lvDTjvB8OG17ysiUmAKIxEAHHkkvP46zJmTdiQiIlmlcBKBqodERKpVOIlgm21gxx3h0UfTjkREJKsUTiIAOO44ePNNmDUr7UhERLJGYSWCY48Nr7orEBFZqbASwWabwd57w9Ch4PWe4kBEJK8UViIAOP54mDkTpk1LOxIRkaxQeImgd29o3DjcFYiISAEmglatoEcPeOQRWLEi7WhERFJXeIkAQu+hzz+H115LOxIRkdTVmgjM7FwzW9+C+8xsspkdkIngEnPIIbDOOqoeEhEh3h3BKe7+PXAAUAr0BQYmGlXS1lkHDj0Uhg2jePlPaUcjIpKqOInAoteDgAfcfWqlbbnr+OPh22/51ccT045ERCRVcRLBJDN7npAInjOz9YBaW1nNbFMzG2dm75nZDDM7N9p+lZnNMbMp0XLQ2hWhnsrLoW1bjnz3hVQuLyKSLeJMXn8qsAswy90XmtlGhOqh2iwD/uDuk6PkMcnMxkSfDXL3v9Yv5AbSuDGceCK/uuGv8OWX0Lp1quGIiKQlzh3BGHef7O7zAdz9G2BQbQe5+1x3nxy9XwC8B7Rbm2Ab3Mkn09hXqNFYRApajXcEZtYMaA60MrOWrGoXWB/4RV0uYmZlwK7AeKALcLaZnQRMJNw1fFvNMf2AfgDt27evy+Xi23Zb3mnbkV0feADOPx9szU0fXQa+yJz5i1aut2tRwuv9uyUTm4hIhqzpjuAMYBKwTfRasTwF3B73Ama2LjACOC/qfXQnsCWhumku8LfqjnP3we7eyd07lZaWxr1cnQ3bqTtMnw6TJtW675z5i5g9sOfKpXJSEBHJVTUmAne/2d03By509y3cffNo2dndb4tzcjMrJiSBoe7+RHTeL919ubuvAO4B9myActTb09vsA82awQMPpBmGiEhqam0jcPdbzWwvMzvOzE6qWGo7zswMuA94z91vrLS9baXdDgOm1yfwhvJ9s3XhsMPgH/+AxYvTDEVEJBVxnix+CPgrsDewR7R0inHuLsCJQLcqXUWvN7N3zWwa8Cvg/HpH31D69oX582HkyLQjERHJuDjdRzsB27nXbQB/d3+N6h88G12X82REt26w6aZw//1w1FFpRyMiklFxEsF0oA2hYTc/FRXBySfDtdfC7NlQVgZU30tIRCTfxEkErYCZZjYBWFKx0d0PSSyqNJx2GgwYAPfeGxICq3oJiYjksziJ4Kqkg8gK7dvDQQfBfffBlVdCcXHaEYmIZEScXkMvA7OB4uj928DkhONKxxlnwH//q0ZjESkocXoNnQ4MB+6ONrUDnkwyqNQceGBoNL7rrrQjERHJmDhjDf2O0BX0ewB3/xDYOMmgUlNUBKefDi+8AB99lHY0IiIZEScRLHH3pRUrZtYYqFNX0lzSa1FHllkj7jr+j+olJCIFIU4ieNnM/gSUmFk5MAz4V7JhpWfqinVofNihnDnrFV4/v0va4YiIJC5OIugPzAPeJQxENxq4LMmgUnfGGfD11/DEE2lHIiKSuDjdR3sBf3f3e5IOJmt07w5bbQW33QbHHpt2NCIiiYpzR3AI8G8ze8jMekZtBPmtUSM4+2x44w2YqDmNRSS/xXmOoC+wFaFt4DjgYzO7N+nAUte3L6y3Htx8c9qRiIgkKs4dAe7+E/AM8ChhcppeSQaVFdZfPySDxx4LD5mJiOSpOA+U9TCzB4GPgN7AvUDbNR6UL84+G5Yt0wNmIpLX4twRnEx4knhrd+/j7qPdfVmyYWWJDh3C+EN33glLltS+v4hIDorTRnAM8A6wD4CZlZjZekkHljXOPRe++ipUEYmI5KH6jDW0Cfk61lB1uneH7bYLjcZ1m5tHRCQnxOkK+jvCBPPjIYw1ZGZ5M9ZQuxYllPUftdr6aszgnHPgzDPh1Vdh330zHKGISLLiJIIl7r40zEWff2MNvd6/W+07nXgiXHYZXH+9EoGI5B2NNRRH8+bhrmDUKJg+Pe1oREQaVGJjDZnZpmY2zszeM7MZZnZutH1DMxtjZh9Gry3XpgAZc9ZZISHccEPakYiINKg4vYZWuPs97n6ku/eO3sepGloG/MHdtwU6A78zs+0IiWWsu3cAxkbr2W+jjcJcBf/4B3z2WdrRiIg0mFhPFteHu89198nR+wXAe4TZzXoBQ6LdhgCHJhVDgzv//NBzaNCgtCMREWkwiSWCysysDNiV0POotbvPhZAsqGG2MzPrZ2YTzWzivHnzMhFm7TbbLIxGOngwfPtt2tGIiDSIGhOBmT0UvZ67Nhcws3WBEcB57v593OPcfbC7d3L3TqWlpWsTQsO6+GL48Ue44460IxERaRBruiPY3cw2A04xs5ZRI+/KJc7JzayYkASGunvFLC9fmlnb6PO2wFdrU4CM23HHMMn9zTdTsnRx2tGIiKy1NSWCu4BngW0II45WXmodpN/Cgwf3Ae+5+42VPhoJ9Ine9wGeqnvYKbvsMpg3j+OnjE47EhGRtVZjInD3W6IeP/e7+xbuvnmlZYsY5+4CnAh0M7Mp0XIQMBAoN7MPgfJoPbfstRd0784Z45+AhQvTjkZEZK3U+mSxu//WzHYmGnQOeMXdp8U47jXAavh4//ghZqkrr6R0n33g7rtDbyIRkRxVayIws3OAfkBFHf9QMxvs7rcmGlm223tvXt9sJ7pcf30Yh6ikZI27dxn4InPmL1q53q5FSbzhLUREEhan++hpwC/d/Qp3v4LwcNjpyYaVG27Z69gwe9ngwbXuO2f+ImYP7LlyqZwURETSFCcRGLC80vpyaq7yKSjj2+8I++0H110Hi9WDSERyU5xE8AAw3syuMrOrgLcIvYEE4MorYe7cWHcFIiLZKM5YQzcCfYH/Ad8Cfd39pqQDyxldu4ZlwAD44Ye0oxERqbNYQ0y4++SoO+nN7v5O0kHlFDP4y1/CdJY3KT+KSO7JyFhDea9zZ+jVKwxR/c03aUcjIlInSgQNZcAAWLAABube83EiUtjWmAjMrMjMXshUMDlt++3DlJa33gqff552NCIisa0xEbj7cmChmW2QoXhy29VXw4oVcM01aUciIhJbnMnrFwPvmtkY4MeKje5+TmJR5aqyMvjtb+H228OwE9tum3ZEIiK1itNGMAq4HHiF1Ucglepcdhmssw5ceGHakYiIxBJn0LkhZlYCtHf3DzIQU24rLQ3J4OKL4fnn4YAD0o5IRGSNar0jMLODgSmEuQkws13MbGTSgeW0c86BLbaACy6AZcvSjkZEZI3iVA1dBewJzAdw9ynA5gnGlPuaNg3PFMyYAffem3Y0IiJrFCcRLHP376ps8ySCySuHHRYGpLv8cviu6j+fiEj2iJMIppvZcUCRmXUws1uBNxKOK/eZwY03hieNr7027WhERGoUp/vo74FLgSXAI8BzwJ+TDCpfdHl+Pufs0J3DbxzEXudut9pn7VqUUNZ/1Grra5qopvLENoUyqU0hllkkDXF6DS0ELjWz68KqL0g+rPwwZ/4ijh47FDp25B9THgI/OdwpwM9+1ConhZrONXtgz1j75otCLLNIGuL0GtrDzN4FphEeLJtqZrsnH1qeKC2F//s/GDcOHnkk7WhERH4mThvBfcBZ7l7m7mXA7wiT1ayRmd1vZl+Z2fRK264yszlmNiVaDqp35Lnk9NOhUyf4wx/UcCwiWSdOIljg7q9WrLj7a0Cc6qEHgR7VbB/k7rtEy+h4Yea4oiK480748ku44oq0oxERWU2NicDMdjOz3YAJZna3mXU1s/3M7A7gpdpO7O6vEGY1Ewh3BGeeCbfdBu9obh8RyR5raiz+W5X1Kyu9X5vnCM42s5OAicAf3P3b6nYys35AP4D27duvxeWyyIABMGJEqCp66y1oHKfTlkDde1mJSHw1/hK5+68SuN6dhK6nHr3+DTilhusPBgYDdOrUKT8eYGvZMtwRHHVUeMbg4ovTjihn1LWXlYjEV+ufpGbWAjgJKKu8f32GoXb3Lyud9x7g6bqeI+f17h2eOr7ySjj0UNh667QjEpECF6exeDQhCbzLWg5DbWZtK60eBkyvad+8ZRbmK2jWDE47LUxkI7Vbvhw++QQWLkw7EpG8E6eSupm7X1DXE5vZI0BXoJWZfU5oY+hqZrsQqoZmA2fU9bx5oW3bUDV0yilw111w1llpR5TdJkyA446Djz+G5s3DkB3eIe2oRPJGnETwkJmdTqjGWVKx0d3X2CPI3Y+tZvN9dQsvj518cnjA7I9/hJ49044me02aBN26wcYbwx13wOjRcMEFnLx/P+A3aUcnkhfiJIKlwA2E8YYqGm0d2CKpoAqCGQweDDvsAH37Ynucn3ZEWadk6WI45pjQyP7GG9CmDZxxBhx+OJeNvBfePRt23BFYfVwiUK8ikbqIkwguALZy96+TDqbglJXBzTfDaadxKpsDB6cdUVbpN+EJ+OijMDxHmzZhY6NGcP/9LGhXRstzz4WxY8FstXGJQL2KROoiTmPxDEAtdEk55RQ49FAueuXvMG1a2tFkj3nzOGPCiNDLqmvX1T/bcENu2vu4kCBeeSWV8ETySZxEsByYEj1dfEvFknRgBSOqIvq+2bpw/PGweHHaEWWHO+6g+U9L4Jprqv340Z0OCAP63XBDhgMTyT9xEsGTwADCZDRr1X1UalBaykUHngfTp8Oll6YdTfoWLYLbb2fslnvAtttWu8uS4qaht9Xo0fDppxkOUCS/xJmPYEgmAskH1TVYxvXSlp3CD9uNN0KPHlBenkSIueHxx2HePO7d/3z2r2GXdi1K6PLpZrzqcNMJl9KuZ7UPqItIDHGeLP6EasYWcnf1GqqiaoNlnd1wA7z0EpxwAkyZEp43KEQPPAAdOvBm+x1r3GVlj6APh3LBrDe44OK/Zyg4kfwTp2qoE7BHtOwD3AI8nGRQBat5cxg2DH74AY49FpYtSzuizJs1C15+OTxnEc3mtkZ9+oQnjt98M/HQRPJVrYnA3b+ptMxx95sAddBOynbbhbkLXn4Zrroq7Wgy78EHQxfRk06Kt/8hh0BxMfzzn4mGJZLP4kxVuVulpZOZnQmsl4HYCtdJJ8Gpp4YpLp97Lu1oMscdHn4YuneHTTaJd8wGG8D++4dE4PkxSK1IpsV5oKzyvATLCGMEHZVINLLKrbeGMXZOOAEmT047msyYPDlU81x+ed2OO+yw8MTx9OkrnzQWkfji9BpKYl4CqaLaiVeGDQszmx1+OE33vSTF6DJk2LAwWU+vXnU7rlevMPvbP/+pRCBSD3F6DTUFjuDn8xFU/6SP1Eu1E6907AhDh0KvXvxl0W3w18PiNaDmIncYPjwMMLfhhnU7tnVr6NwZ/vUvzQktUg9xeg09BfQiVAv9WGmRTDjkEPjznzl8xjgYNCjtaJIzdWoYZvrII+t3/IEHhpFKv9aQWCJ1FScRbOLuR7v79e7+t4ol8chklUsvZfTWe8FFF8Hzz6cdTTKGD4eiojBrW338+tfhrmLMmIaNS6QAxEkEb5iZKl7TZMaFPc+H7beHo49m8//NSTuihuUe2ge6doVWrep3jt13D1VKhdTLSqSBxEkEewOTzOwDM5tmZu+amYbJzLCFTUrgqaegcWMeGHYVzJuXdkgNZ/p0+Pe/618tBOFuorw8JAJ1IxWpkzjdRw9MPAr5mep6EbH55vCvf9Fmn/1C28GLL0JJSe5PyjJ8ODRqxG/+sxHTq5a5Lnr0gMce03DeInUUp/uohnZMQY0/5J07c+7BF3L3k38Jzxg8/njuT8oyfDjssw/Tl5es3VhNBxwQXp97Dti+QUITKQRxqoYkyzy39V6hB9ETT4QG5Fz2/vswcyYcccTan+sXvwhTf6rBWKROEksEZna/mX1lZtMrbdvQzMaY2YfRa8ukrp/3zj0XzjsPBg2i3/gRaUdTfyOi2A8/vGHOV14Or75K05+WNMz5RApAkncEDwI9qmzrD4x19w7A2Ghd6uuvf4Vjj+VPLz0Ad9+ddjT1M2JEeBisXbuGOV/37rBkCXt8PrNhzidSABJLBO7+CvC/Kpt7ARUT3QwB6tlpXIDQU2bIkDCT129/C488knZEdTNrFrzzTsNUC1XYbz8oLmbvT6c03DlF8lycXkMNqbW7zwVw97lmtnFNO5pZP6AfQPv27TMUXg4qLuasXv35YPKtYdTS9dYDcmQYiieeCK8NmQjWWQf22ou931ciEIkraxuL3X2wu3dy906lpaVph5PVlhQ3hZEjYdddoXdv9vkkR0YrHTECdtstdIttSOXl7PDlx/n1rIVIgjKdCL40s7YA0etXGb5+/lpvPXjmGdh2W+4d8Wd49tm0I1qzOXPgrbca9m6gQsV8z2PHNvy5RfJQphPBSKBP9L4PYUA7aSgbbQRjx/Jhq/ZhaOZRWfw8weOPh9ckEsHuuzO/2brqRioSU5LdRx8B3gQ6mtnnZnYqMBAoN7MPgfJoXRrShhty3DEDYKedwoQtI0emHVH1Hn44zLXQsWPDn7uoiDfa7xQSgYabEKlVYo3F7n5sDR/tn9Q1M63aYSCywPcVfw3/+tfhL+5HHoHevTMeR41DX7z3HkyezM09z2JQQv9+r5ftwkHPvwEffghbb91g5xXJR5nuNZRXsno8nxYtwpDVBx0ERx0Fd9wRZvHKoBqHvhg6FBo14uGy/7d2Q0qswatlu4Y3Y8YoEYjUImt7DUkD2GCD8EN40EHhOYOrr06/qmTFipAIundn3rrJPVj+n5ZtQ28ktROI1EqJIN81bx7m8u3TB666Cs4+G5YvTy+el16C2bPDgHlJKy+HceNg2bLkryWSw5QICkFxMTzwAFx8cagiOvpoWLgwnVjuugtatsxMm0V5OXz/PUyYkPy1RHKYEkGhMIPrrls1aul++8EXX2Q0hNIf/hfuTvr2hZIMNKx36xbKreohkTVSY3EOWqveSuedB1tuCccdB3vsQd/DL2fcOpuudq6kGsFP/3eopum6oCOz+49KvpfVhhuGKSzHjIErr0z2WiI5TIkgB631D/XBB8Prr8PBB3PHXedR8ujQlQ92JTapzcKF9Js6Cnr04KXB/ZK5RnXKy+H660MV0frrZ+66IjlEVUOFaqedYMIE3tu4LNTXX3JJso2q994LX30Fl16a3DWqU14eGsdffjmz1xXJIUoEhax1a4499i/Qrx8MHAjl5bT68duGv86SJeGv8n33hb33bvjzr8lee4WeU2onEKmRqoYK3JLGTcKkNl26wJlnMmriNDhqM9hnn4a7yE03hUHmhgypfd+G1rRpSEBKBCI10h2BBCedBOPH82NxM+jaFS67DJYuXfvzfvEF/PnPcMghsH9Ko4uUl4e5kT//PJ3ri2Q53REUoMpjAK3Wc2fHHTmkz01Mnz8aBgwIw1o/9BBst129rmO+Ak49NbQ93HhjQ4ReP927h9cxY0LXVRFZje4IClDFGECzB/b8WQ+kH5o2Dw+fPfEE/Oc/YeKYa6+FxYvrfJ3TJjwZ5kX4299Cl9W07LgjtG4NL7yQXgwiWUyJQKp32GEwfXroanr55eHHdPTo+GMVPfwwl750f+iWetZZycZaG7NwV/DCC2GsIxFZjRKB1Kx1axg2LIxiWlQEPXuGXjhPP13zeEVLl4aHt048kTfb7xjmHbAsmEO5vDx0X506Ne1IRLKO2gikduXlMG0a3H9/6GZ68MHwi1+Ev/Y7d4ZNNgljF739dqhW+uQTOOkkTi49lA+aNUs7+uDAA6FRozDExa67ph2NSFbRHYHE06RJmM/gww/hscfC7GL33APHHx/GLTrwQLjiCmjXLrQLDBkSuqZmi403DnEOH552JCJZR3cEUjfFxWGim6OOgp9+CrONzZsHTZvS89mvmLGsKYxbBuMyMJbQGlQ3HtPrvXvD734HM2eu1hOqxpnURAqEEoHUX3FxGKoiMuPpUYnNOFZXVX/Iy/qPgj6HhfkYhg8Pdy+RGmdSEykQqVQNmdlsM3vXzKaY2cQ0YpAC1LZteIJ6xIi0IxHJKmm2EfzK3Xdx904pxiCF5sgjQ8P3jBlpRyKSNdRYLIXl2GOhcePQu0lEgPTaCBx43swcuNvdB1fdwcz6Af0A2rdvn+HwclvVhtLqPl+bYyvXv9c4XEW2Ki0N3V8fegj+8pfQziFS4NJKBF3c/Qsz2xgYY2bvu/srlXeIksNggE6dOsV8nFVg7Sauqe3YqkmiakNrTujbNzxP8MwzYTA8kQKXStWQu38RvX4F/BPYM404pEAdeCC0aQN33pl2JCJZIeOJwMzWMbP1Kt4DBwDTMx2HFLDGjcPzBM8+G54pEClwadwRtAZeM7OpwARglLs/m0IcUsjOPBOaNQuT5ogUuIy3Ebj7LGDnTF9XZDWtWkGfPvDgg5Se1oCzsYnkID1ZLAWhcm+olT2fLroI7ruPc19/BDix2n1X2z+GqsNVVBeHhq+QbKNEIAWh8o/vyh/5LbeEfv045s67wmB6HTr8bN/V9o+htl5UGr5CspEeKJPCdvnlLG1cDBdeGH/SHZE8o0Qgha1NG27qchyMHAmPP552NCKpUCKQgnffHofCHnuEkUnnzk07HJGMUyKQgre8URE8+GCYZe2II2DJkrRDEskoNRZLnVTXoybXVFuG7baDv/8deveGU04J74uKqj2+6vhK6gUkuU6JQOokH370aizDEUfA//0f/OlPoeH4wQfDFJ1VVO4ZpF5Akg+UCEQqu+QSMAuvn34a5mcWyXNqIxCpqn//kACmToXtt+fkiSNh2bK0oxJJjBKBSHWOOgreeQc6d+aqsYNh663httvgu+/SjkykwSkRiNSkQwd49llOPeLyMGz1738PrVtz1z8HwNCh6moqeUNtBCJrYsb7nbpSNv+X7LzVB/Sa+TKH/Ps1OOEEAMZt3J6hz23HtDYdmN5mKzbbsmOdTl+XHki1jWNUWa70ZqpaplyJO98oEYjUYtUPU0/gAli+HKZMgRdfZPOxY9n8rTc5fko0knqTJvDkjrDDDqFLasVSVgaNfn4DXpceSHWZDS5XejNVLVOuxJ1vlAhE6qqoCHbfPSwXXQQrVsDHH8OkSWF55x14/nkYMmTVMSUlsM023Lx4fWgyPgx4t9VWlP7wbeiqapZeeaTgKRGIrK1GjUJ7QocOcMwxq7Z/+y28916YBW3mTJgxg93GT4UBr4bkAbwN8OCZsOWW3Ll4XVjxMmyxBbRvv2pZf/1UiiWFQ4lAJCktW8Jee4Ulsk//Ucy+phxmz4aPP+bKG5/i6u2bwUcf0eHNqXDzRFi6dPXzbLABbLop9y9sBt89DZtuGhLEJptA27ZhWW893VVIvSkRiGRakyahO+rWWzPk5RVcHdWR9xn4InP/9wMb//Atv1gwjx2Wf8c1u64P//kPfPYZG781HYYNg2+++fk5S0pCQmjTBtq25W9zlnHDvo8xb52WfLVuSxq1acP9Fx4YZmZr3pwu142L3Uhb1wY/eZZTAAAJw0lEQVTdNTVq17UxuCEbk9UwXTMlApEsUd2EONdctKoh9Tf9R4WG1YUL4bPPwvLf/4Zl7txVrzNncsTcuTB//uoXuC16bdqU4cXr0nbLTWCjjaBVK/7+4Y+w6OWV62y0EbRoAS1asPTzL5j9195hjmfWrlG7ro3BDdmYrIbpmikRiOSa5s2hY8ewrMnixSsTRb/rn2Zwz83D3cTXX/Pqs+9w1GbN4euvYcoUfvPZf2HK6Gon53kb4PYTw51Mixa8+FNjeLFdSBQbbLD66/rrc+S0j+HxH2HddWGddcJrtKy35MfwlHZj/fRkk1S+DTPrAdwMFAH3uvvANOIQyWvNmoVuq2VlPL/1N3Dqqr+GL240iqMq/XW8W/9RzB7QIzRwR8mC776D+fO5bMjrXNt1k3CH8d13zHh5JltstE5Y/+yzlfuxKFS73ADwzC3VhvQuwE1A06aw7rq89lMRPF26Mmnc++kCmDUkVHU1a8YVk/4Ll7wWylJSQp9Js+CeL8Ln0T6rvVZ+36TJqkWJZ40y/q9jZkXA7UA58DnwtpmNdPeZmY5FRCopKgrVQq1arXa38fC0Dbj2klVJ4/f9R3FwdVU/S5fC99/T5fJ/8frZe8IPP4Tlxx9Xvv/zoxO4/FebrVwf/8r7bLJ1i7DPggW0XfANvLsgJJXFi+n97QKY9gz89BMAVwO8UL/ifVDUGG5vtjI5vLZwOYzYIKwXF6+eOGpaiotDUqm8FBUlt62oCDbbLCTKBKWRJvcEPnL3WQBm9ijQC1AiEMllTZpAq1bM2WBj2H77ane5b1Y7Lr98VRL5Q/9RHFEpqfSsaAeJ7FSxvnw5LF7Mrpc8xTsX7xsSRZQsfva+4vWnn0Jyipb7np/JWXu1X7n+5huzOHKnjVfbZ+Xy44+r3lc+z5IlIZbly0MVV8WSpGeegR49Er2EeYYn7Daz3kAPdz8tWj8R+KW7n11lv35Av2i1I/BBlVO1Ar5OONxMyrfyQP6VKd/KA/lXpnwrD6xdmTZz99LadkrjjqC6zs4/y0buPhgYXONJzCa6e6eGDCxN+VYeyL8y5Vt5IP/KlG/lgcyUKY3RRz8HNq20vgnwRQpxiIgI6SSCt4EOZra5mTUBjgFGphCHiIiQQtWQuy8zs7OB5wjdR+939xn1OFWN1UY5Kt/KA/lXpnwrD+RfmfKtPJCBMmW8sVhERLKLZigTESlwSgQiIgUuqxOBmfUwsw/M7CMz61/N503N7LHo8/FmVpb5KOsmRplONrN5ZjYlWk5LI864zOx+M/vKzKbX8LmZ2S1ReaeZ2W6ZjrEuYpSnq5l9V+n7uSLTMdaFmW1qZuPM7D0zm2Fm51azT659R3HKlDPfk5k1M7MJZjY1Ks/V1eyT7G+du2flQmhI/hjYAmgCTAW2q7LPWcBd0ftjgMfSjrsBynQycFvasdahTPsCuwHTa/j8IOAZwvMjnYHxace8luXpCjyddpx1KE9bYLfo/XrAv6v5by7XvqM4ZcqZ7yn6d183el8MjAc6V9kn0d+6bL4jWDkUhbsvBSqGoqisF1AxH+BwYH+zrJ6dI06Zcoq7vwL8bw279AL+7sFbQAsza5uZ6OouRnlyirvPdffJ0fsFwHtAuyq75dp3FKdMOSP6d/8hWi2Olqq9eBL9rcvmRNAO+KzS+uf8/MteuY+7LwO+AzbKSHT1E6dMAEdEt+jDzWzTaj7PJXHLnEv+X3Qb/4yZVT+oThaKqhN2JfzFWVnOfkdrKBPk0PdkZkVmNgX4Chjj7jV+R0n81mVzIogzFEWs4SqySJx4/wWUuftOhHEWh/z8kJySa99RbSYTxm/ZGbgVeDLleGIxs3WBEcB57v591Y+rOSTrv6NaypRT35O7L3f3XQgjLexpZjtU2SXR7yibE0GcoShW7mNmjYENyO7b+lrL5O7fuPuSaPUeYPcMxZaUvBpSxN2/r7iNd/fRQLGZtUo5rDUys2LCD+ZQd3+iml1y7juqrUy5+D0BuPt84CWg6nCjif7WZXMiiDMUxUigT/S+N/CiR60pWarWMlWpmz2EUP+Zy0YCJ0U9UzoD37n73LSDqi8za1NRN2tmexL+H6pmEuHsEMV6H/Ceu99Yw2459R3FKVMufU9mVmpmLaL3JUB34P0quyX6W5e10/Z4DUNRmNk1wER3H0n4j+EhM/uIkB2PSS/i2sUs0zlmdgiwjFCmk1MLOAYze4TQQ6OVmX0OXElo7MLd7wJGE3qlfAQsBPqmE2k8McrTG/itmS0DFgHHZPkfH12AE4F3ozpogD8B7SE3vyPilSmXvqe2wBALk3Y1Ah5396cz+VunISZERApcNlcNiYhIBigRiIgUOCUCEZECp0QgIlLglAhERAqcEoGISIFTIhARKXBKBCJZwswONbN7zOwpMzsg7XikcCgRSMEzszIzW1TpKVXMrLWZ/cPMZpnZJDN708wOq+U8L5nZr6tsO8/M7qi0freZdanueHd/0t1PJzxNfnS0f0k0scrSXBgrR3KTEoFI8HE0+mPFWDZPAq+4+xbuvjvhkf5NajnHI/z80f9jou0Vfgm8Vct5LgNuB3D3RVFcWT0InOQ2JQLJO9Ff+O+b2b1mNt3MhppZdzN73cw+jAYhW5NuwNJozBoA3P1Td7+10jVOiKYXnBL9lV9EmDDkN2bWtCIO4BfAa9H6tsC/3X25mZ0UzTkx1cweij43M7sOeKZi4hWRTFAikHy1FXAzsBOwDXAcsDdwIWGAsjXZnjCefbWiH/SjgS7RX+vLgePd/RtgAquGEK6YUrBiQK8DgWejSVIuBbpF4+VXzLn7e8LIk73N7Mw6lFVkrWTt6KMia+kTd38XwMxmAGPd3c3sXaCsLicys9sJSWSpu+8B7E+YJ+LtaKTjEsLMUrCqeuip6PWUSqf6NWFkzyOA4e7+NYC7/y96vQW4pc4lFVlLSgSSr5ZUer+i0voKav/vfgbhxxoAd/9d1FA7MdpkwBB3v6SaY58EbjSz3YCSiioeM2sOtHD3L6I2CA37K1lDVUMiP/ci0MzMfltpW/NK78cSqm82BjCzDc1sM4BoVqyXgPtZvZH4V8C4SscfZWYbVRyfRCFE4lIiEKkiqtM/FNjPzD4xswmEuaP/GH0+k9Cz53kzmwaMIUwuUuERYGfg0UrbDgSejY6fAQwAXjazqUBNM4eJZIQmppGCF/Xuedrdq04Y3pDXmAz80t1/qufxs4FOFe0KIg1JdwQiodfPBpUfKGto7r5bfZJAxQNlhOkyVzR8ZCK6IxARKXi6IxARKXBKBCIiBU6JQESkwCkRiIgUOCUCEZECp0QgIlLglAhERAqcEoGISIH7/4kNrLxq42RKAAAAAElFTkSuQmCC\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_ci_ex.png\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Confidence interval on mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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": null,
   "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 = -2349.8670095885573</td>\n",
       "        <td title=\"Total number of call to FCN so far\">TOTAL NCALL = 62</td>\n",
       "        <td title=\"Number of call in last migrad\">NCALLS = 62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td title=\"Estimated distance to minimum\">EDM = 5.373753935591247e-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=\"$('#lkrQCVKxnd').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.94355</td>\n",
       "        <td>0.119218</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>mean</td>\n",
       "        <td>1.20658</td>\n",
       "        <td>0.0178857</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.5</td>\n",
       "        <td>2</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>2</td>\n",
       "        <td>sigma</td>\n",
       "        <td>0.0995381</td>\n",
       "        <td>0.0207983</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0.02</td>\n",
       "        <td>0.2</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>3</td>\n",
       "        <td>Nsig</td>\n",
       "        <td>73.9469</td>\n",
       "        <td>13.1775</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>490</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "        <td>4</td>\n",
       "        <td>Nbkg</td>\n",
       "        <td>416.003</td>\n",
       "        <td>22.6921</td>\n",
       "        <td></td>\n",
       "        <td></td>\n",
       "        <td>0</td>\n",
       "        <td>539</td>\n",
       "        <td>No</td>\n",
       "    </tr>\n",
       "</table>\n",
       "<pre id=\"lkrQCVKxnd\" style=\"display:none;\">\n",
       "<textarea rows=\"16\" 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.94355 & 0.119218 &  &  & -4.0 & -1 & No\\\\\n",
       "\\hline\n",
       "1 & mean & 1.20658 & 0.0178857 &  &  & 0.5 & 2 & No\\\\\n",
       "\\hline\n",
       "2 & $\\sigma$ & 0.0995381 & 0.0207983 &  &  & 0.019999999552965164 & 0.2 & No\\\\\n",
       "\\hline\n",
       "3 & Nsig & 73.9469 & 13.1775 &  &  & 0.0 & 490 & No\\\\\n",
       "\\hline\n",
       "4 & Nbkg & 416.003 & 22.6921 &  &  & 0.0 & 539 & 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 0x1a414d2438>"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "config.bestfit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "calc = AsymptoticCalculator(config)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "poinull = POI(mean, value=np.linspace(1.15, 1.26, 100))\n",
    "ci = ConfidenceInterval(poinull, calc, qtilde=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Compute qobs for the null hypothesis!\n"
     ]
    }
   ],
   "source": [
    "ci.interval();"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f, ax = plt.subplots(figsize=(10, 8))\n",
    "ci.plot(ax=ax)\n",
    "f.savefig(\"ci_1_cl_plot.png\")"
   ]
  },
  {
   "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
}