{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pyhf\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import getpass\n",
    "from matplotlib import rcParams"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(0)\n",
    "rcParams.update({'figure.autolayout': True})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Let's look at an example of a well seperated signal and background and see that pyhf makes lots of pseudoexperiments fast.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pyhf.simplemodels.hepdata_like([20], [50], [np.sqrt(50)])\n",
    "\n",
    "signal_pars = model.config.suggested_init()\n",
    "signal_pars[model.config.poi_index] = 1.0\n",
    "\n",
    "bkg_pars = model.config.suggested_init()\n",
    "bkg_pars[model.config.poi_index] = 0.0\n",
    "\n",
    "signal_pdf = model.make_pdf(signal_pars)\n",
    "bkg_pdf = model.make_pdf(bkg_pars)\n",
    "\n",
    "sample_size = 1000000\n",
    "sample_shape = (sample_size,)\n",
    "\n",
    "sample = signal_pdf.sample(sample_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de5hU1Znv8e8bRIKGy5wgXhqhlYtoGo9K03gGnVwQQREVDVGcM4YRZTRqCNGJgB7BBBrmJBNjQJMgGuIkgsaAoQEvgCFiNFxyhksD3uLh0q2HThRQx9aAvuePqmqLpqq7q3tX7apdv8/z9NPN2rvWfin07V1rr/Uuc3dERCT3PhN2ACIixUoJWEQkJErAIiIhUQIWEQmJErCISEiOCjuAbOvWrZuXlpaGHYaIFKk//elPf3X341Idi3wCLi0tZePGjWGHISJFysx2pTumIQgRkZBENgGb2Sgzm3fgwIGwQxERSSmyCdjdq9x9QpcuXcIORUQkpcgmYBGRfKcELCISksjPgmjOu+++S11dHQcPHgw7lKLSvn17unfvTufOncMORSQ0kU3AZjYKGNWnT5+057z77rvs3buXkpISOnbsiJnlLsAi5u7U19dTW1sLoCQsRSuyQxAteQhXV1dHSUkJxxxzjJJvDpkZxxxzDCUlJdTV1YUdjkhoIpuAW+LgwYN07Ngx7DCKVseOHTX0I0WtqBMwoDvfEOm9l2JX9AlYRCQsSsAiIiGJbALWUmQRaat7V74KQMXMVQyZ/Vzg/Uc2AWspsoi01X2rXwNg/Z0XULu/PvD+IzsPuC2GzH6O2v31dO/UgfV3XsC9K19t+IcAqLrlPABGzX2hoW3i0L5MGtaPipmrqHvvIwDKSjqz7NbzmbJ4CwvX72k4d93UoWytOcD1j8TKZJZ07cgfJn+lVbE++OCDzJo1i927d3Pdddcxb948APbt20f//v158cUX6d27d6v6BlizZg3jxo1j586dGb1uzJgxnHvuudx2222tvrZIPikrycJ8dXeP9NfAgQM9ne3bt6dsH/nj59O+Jht+9cddrXrdjh07vF27dr548WJ/8803/b333ms4dvvtt/u4cePaHNvvfvc779Wr12Ftv//9733UqFF+0kknOeA///nPj3jdli1b/O/+7u98//79Tfaf7t9AJB/0umNZm/sANnqa/BTZIYi2WHbr+Tm93jWDe7bqdUuXLqWsrIzRo0dz4okn8rnPfQ6ADz74gPnz5zN+/Pggw2zw/vvvU1ZWxn333Zd2HvWAAQM49dRT+eUvf5mVGERyIfFpF2DK4i2B968EnEI23uimlE5envFr+vXrxx133MHmzZsxM0aPHt1wbMWKFZgZQ4YMaWh74okn6NChA7t2fVqcf+LEifTu3Zu9e/dmdO2LL76YyspKvvrVr/KZz6T/T+jSSy9l4cKFGfUtkq+ShxGDogScQjbe6KC98MIL9OvXjxkzZvDWW2/xi1/8ouHY2rVrGThw4GELHa688koGDBjAjBkzAPjBD37AwoULefrppzn++OOzEmNFRQXr16+nvj74hxciuZD8nCcb9BCuQHXu3Jk33niDIUOGcMIJJxx2bNeuXZx00kmHtZkZlZWVjBw5kt69e1NZWcnq1avp27dv1mI86aSTOHjwIG+++WabHgSKRFVk74ALaR7w0P7dM35NdXU1hw4d4qyzzjriWH19PZ/97GePaL/wwgsZNGgQd911F4899hiDBg1qVbwtlRgf1h2wFJKKmauA2Bzgkq6fPuNYN3Vo4NeKbAL2NswDzsYb3ZSHxmWeCDdt2kSvXr3o2rXrEce6devGvn37jmh/7rnn2Lx5M+5+xLDDgw8+yDnnnENZWRlXXXVVxvGk8s477wBw3HEpd+QWyUuJaaSThvU7bHro1prgb+Yim4DbIhtvdFPGL9iQ8Ws2bdqU8u4X4Oyzz2b79u2HtW3evJnRo0czZ84cLr/8cqZMmdJwbN++fdx///1s2LCB6upqfvazn2UcTyrV1dWUlJRkbYxZJJcS8/aDpAScQjbe6KasfjnzmrhNJeDhw4ezY8cO3n77bSA2JnzRRRdx2223cd1113HPPfewcuVK1qxZA8BRRx3Fvn37+M53vsO2bdtS3lUne//999m0aRObNm3ik08+Yffu3WzatIndu3cfdt7atWsZPnx4xn83kTBlZcFFGnoIl0by1LCgV74BVI4ewDWDe1I6eflh40wt4e5s2bKF22+/PeXxAQMGUFFRwaJFixg7diwjRoxg1KhR3H333bE4y8oYM2YMU6ZM4aWXXqJTp05UV1fz5JNP8rWvfY2ZM2dy+eWXp73+xo0b+fKXv9zw52nTpjFt2jS+/vWvs2DBAgA+/PBDlixZwjPPPJPR300kbDldB5BuhUZUvlqzEi4KnnrqKe/Xr58fOnSo2XNfffXVhp9vuukmX7RoUcOfU62Ea4m5c+f6sGHDmj0vyv8GUpgm/2ZzyvbWrlhFK+GKz4gRI7j55pupqalp9twZM2Zw2mmncfbZZ2NmjBkzps3Xb9++PXPmzGlzPyK5lm4dQGtXrDZFQxAR9s1vfrNF5yUv4gjKhAkTAu9TJFtWbd/LgB5dGFy5Ou2QYOnk5eycPTLQ6yoBS5NKS0v51re+FXYYIll1/SMb2Tl7ZOAJtjkagpAmKQGLZI8SsIhIC7RmxWpzNAQhIoXv3gFwYPeR7V16wqStzb68cvSAZs9pzYrV5hRUAjazLwHfA7YBi9x9TagBiUh+OLAbpqdYwTq9ZaUIWjLDYfyCDYEn4dATsJk9DFwC1Ll7WVL7COA+oB0w391nAw68D3wWaH5+lYgUty490yfhpLvjlsxwaM2K1eaEnoCBBcBc4JFEg5m1A+4HhhFLtBvMbCmw1t1/b2bHAz8E/jH34YpIVqUbTmhKlzR3sE0NP7Tw7jibQk/A7v68mZU2aq4AXnf3NwDMbBFwmbsnKszsAzqk69PMJgATAHr2DH7ytIhkUbrhhCwZv2BDxuUAghJ6Ak6jBEhejlIDDDazK4DhQFdid80pufs8YB5AeXm5ZzFOESlwLR3XzcYc4XxNwCm5+2JgcUvONbNRwKg+ffpkNygRyVxTwwzphhNC9ui63YEvR87XBFwLnJz05x7xthZz9yqgqry8/IYgAxORAOR4mCGVGu9Gj1TjwGmmrk1dsrVoEvAGoK+ZnUIs8V4NXJOzq7fmIUBbtHCuYioPPvggs2bNYvfu3Vx33XXMmzcPiBVZ79+/Py+++GKb9mNbs2YN48aNY+fOnRm9bsyYMZx77rncdtttrb62SDb1uOfPqQ/k8OFc6AnYzBYCXwK6mVkNMM3dHzKzW4BniE1De9jdt2XYb+uHIHL927mV/+Avv/wyN910E7/+9a8599xz6dSpU8OxyspKLr744qxthvnAAw/w/e9/n7feeosvfOEL/OhHP+L88z+to3r33XfzxS9+keuvv57WbAslEZCnwwxDZj9H7f76hprcYQo9Abv72DTtK4AVbeg38kMQS5cupaysjNGjRx/W/sEHHzB//nyqqqqyct3HHnuMiRMn8sADD3DeeefxwAMPcNFFF7F9+/aGWScDBgzg1FNP5Ze//CU333xzVuKQPJcHwwyp1O6vb9UDtfnXlgcei2pBFKh+/fpxxx13sHnzZszssCS8YsUKzIwhQ4Y0tD3xxBN06NCBXbt2NbRNnDiR3r17s3fv3oyu/cMf/pBx48Zxww03cPrppzNnzhxOPPFEfvKTnxx23qWXXsrChQtb+TcUyY7WJtIBPYL/JBfZBFxI29K3xgsvvEC/fv2YMWMGb7311mE1fdeuXcvAgQMxs4a2K6+8kgEDBjBjxgwAfvCDH7Bw4UKefvrpjDbN/Nvf/saf/vQnLrzwwsPaL7zwQl588cXD2ioqKli/fr22pZe80tpEOrhydcCR5MEQRLZEfQiic+fOvPHGGwwZMoQTTjjhsGO7du3ipJNOOqzNzKisrGTkyJH07t2byspKVq9eTd++fTO67l//+lc+/vjjI5L28ccfz6pVqw5rO+mkkzh48CBvvvlm1saiJQ80VQgnDw2uXJ3zur/pRDYBR111dTWHDh1KuTNyfX19yrvaCy+8kEGDBnHXXXdRVVXFoEHBV3dK1rFjx4Z4JMLydKy3EGgIokBt2rSJXr16pdxCvlu3buzbt++I9ueee47Nmzfj7kck6AcffJBzzjmHsrIyrrrqqrTX7datG+3atTti3Hjv3r1H3Im/8847ABx33HEt/nuJ5KuxFSc3f1KGIpuA3b3K3SdEdQrUpk2bUt79Apx99tls3779sLbNmzczevRo5syZw+WXX86UKVMaju3bt4/777+fDRs2UF1dzc9+9rO01z366KMZOHAgK1euPKx95cqV/P3f//1hbdXV1ZSUlGQ0xiySba1NpLOuODPgSCKcgKOuqQQ8fPhwduzYwdtvvw3ExoQvuugibrvtNq677jruueceVq5cyZo1awA46qij2LdvH9/5znfYtm1byrvqZN/+9rdZsGAB8+fPZ8eOHUycOJE333yTG2+88bDz1q5dy/Dhw9v+lxUJUGsT6SVz1gYcicaAU2uqhmi2rpcBd2fLli3cfvvtKY8PGDCAiooKFi1axNixYxkxYgSjRo3i7rvvBqCsrIwxY8YwZcoUXnrpJTp16kR1dTVPPvkkX/va15g5cyaXX3552utfddVVvP322w0zMMrKylixYgW9evVqOOfDDz9kyZIlPPPMMxn93SRP5emiipaqmLmKuvc+AqCspDPLbj2/mVccqbr23aDDwtyjWSwsaSXcDa+99lrKc3bs2MHpp5+e28By5Omnn2bixIls376ddu3aNXnua6+91jAb4hvf+AZf/OIXG8aBW7sU+f777+e3v/0tzz77bJPnRfnfIFKmdymeB23N/bLJsGyAmf3J3VNOPo7sEETUx4CbM2LECG6++WZqaprfOGTGjBmcdtppnH322ZgZY8aMafP127dvz5w5c9rcj0gQ7l35astPnrQ19sum0VdF+98EXiNGQxAR9s1vfrNF5yUv4gjKhAkTAu9TpLXuW/0ak4b1a1Mf6++8AKYHE09CZO+AJRilpaV861vfCjsMkdBldBfdQpFNwFGfB5wrSsAiMfetTv0sqS0im4CLfQxYRD5Vdct5YYeQUmQTsIhIviv6h3DufljVMMmdqE6BlPxw78pXDxs2aGsBnqpbzoP5bY3qcEWdgNu3b099fT3HHHNM2KEUpfr6etq3bx92GJKswCqbNWXSsH5tnvmQbUWdgLt3705tbS0lJSV07NhRd8I54u7U19dTW1urOhH5JkKVzSpmropNHQvIqLkvsPOzgXUHFHkC7ty5MwBvvvkmBw8eDDma4tK+fXuOP/74hn8DyaECX1bcUomlx/kssgm4pZtydu7cWUlAikuE7nILXWRnQWgamkhxKysJ9sZq4tDMdo9picjeAYtI8UlsOT+24uRWVTxryqRh/eAPgXapBCwi0dHaLedbomLmKtYH3GdkhyBEpPismzo0a31n46GeErCIRMbWmsJ6uKgELCKRcf0jG7PWd9AP9UAJWESkRYJ+qAcRfgjX0nnAIpEVoWXF+WDK4i3MCrjPyCZgd68CqsrLy28IOxaRUBTJgovE1LOds0dSOXpA1q6zcP0eZgW8FFlDECJS0JKnnl0zuLDu7pWARaSgDe3fPewQWk0JWEQK2kPjBuXkOtmYY6wELCIFbfyCDTm5TjbmGEf2IZxIUSiS0pJNWf1yXU6uc/0jG1UPWESSFMlMh6jSEISISEiUgEWkoGWr+llj2ZhjrAQsIgXt0XVpxsADlo05xhoDFpGCs2r73obCOyVdO+ZkAUbp5OV6CGdmxwK/B6a7+7Kw4xGR3BvQo0vOhh6yKfQhCDN72MzqzKy6UfsIM3vFzF43s8lJh+4AHs9tlCKSTwZXrg47hEDkwx3wAmAu8EiiwczaAfcDw4AaYIOZLQVKgO1AwB8ERPKY5vrmhaH9u8POYPsMPQG7+/NmVtqouQJ43d3fADCzRcBlwOeAY4EzgHozW+HunzTu08wmABMAevbUf6BS4DTXNy88NG4QTA+2z9CHINIoAfYk/bkGKHH3O939W8CjwIOpki+Au89z93J3Lz/uuONyEK6I5MIlc9YCMLbi5JxfOxtLnkO/A24Nd1/Q3DkqyC4FRwXUm1Vd+y4As644M+fXXv1yXeCDn/magGuB5F9xPeJtLaaC7FJwNNRQdPJ1CGID0NfMTjGzo4GrgaUhxyQiIeveqUPYIQQq9ARsZguBl4DTzKzGzMa7+yHgFuAZYAfwuLtvy7DfUWY278AB3VGIRMX6Oy8I7drZmHccegJ297HufqK7t3f3Hu7+ULx9hbv3c/fe7j6zFf1WufuELl26BB+0iITi3pWvhnbtbCx5Dj0BZ4vugEWi577Vr4V27alLtgbeZ2QTsO6ARSTfRTYBi0g0lE5eztaaA1nZEihs+ToNTUSkwYAesU+yYRbgmX9teeBVaCJ7B6wxYBEJUuKXQJAim4A1BiwSDROH9g07BCA7Fdgim4BFJBomDesXdghZozFgkVxSacmMVcxcFeoCjGyKbAJWMR7JS6r3kLG69z4KOwQgXoFtS7B9RnYIQmPAIoVryuJYprtkzlpKunYMOZqYbFRgi2wCFpHCtXB9rBz4slvP5w+TvxJyNDGJWsRBUgIWEWmBRC3iIEU2AWsesIjku8gmYI0BixSudVOHhh3CEbJRiziysyBEQqXthdpka80Bjj8jvzY/X3/nBYFvyqkELJINmm7WJtc/sjHUug+p3LvyVSYF3GdkhyBERIKUjVrESsAiIiGJbALWLAiRwjJk9nOUTl4OQOXoASFHkxuRTcCaBSFSWGr31zeM+14zOP8eVlbdcl7gfeohnEhrqbBOoIb27x52CDmnBCzSWprpEKiHxg0KO4QmjZr7AjsDnhkX2SEIESks4xdsCDuEnFMCFpG8sPrlurBDyDklYBEJzaPrYmPopZOX503ZyXSysTWSxoBFJDRTl2zlmsE9827VWyqThvWDPwTbp+6ARURaoGLmqsD7jGwC1kIMEQlSNrZGyigBm1lPM+vczDmdzCz0SZBaiCGS/+ZfWx52CKHK9A74/wITmznnm/HzRESaNKBH4dwglZU0ee/ZKpk+hLP4l0jxUG3fQE1ZvKVhzzegIB7AQWx/ukKoB3wC8F9Z6FckHFrxFqhZV5yZlR2Gs23K4i3MCrjPZhOwmV3bqOmsFG0A7YCewP8EtgYQm4hE0CVz1sbuJgvMwvV7mBXwUuSW3AEvADz+swOXxb8aSwxNfADc0+bIRCSSsrG7cKFqSQL+5/h3Ax4GngR+m+K8j4G3gZfcfX8w4YmIRFezCdjdf5H42cy+Djzp7o9kNSoRiaxs7C6cC+umDoUfBttnRg/h3P3LwV5eRIrN+jsvCDuEVtlac4DjA+4zsivhRCQ/3bvy1bBDaJXrH9kYeJ8ZJ2Az+6KZLTOzOjM7aGYfp/g6FHiksWufbmY/NbMnzOymbFxDRLIrG7sLF6qMhiDMbCSxh3DtgN3AK0Cbkq2ZPQxcAtS5e1lS+wjgvvi15rv7bHffAdxoZp8BHgF+0pZri4iEKdOFGNOBg8BId382oBgWAHOJJVQAzKwdcD8wDKgBNpjZUnffbmaXAjcB/xHQ9UW0v1sObK0p7MUslaMHwFPB9plpAi4DFgWYfHH3582stFFzBfC6u78BYGaLiM093u7uS4GlZrYceDSoOKTIabVb1o2a+wI7Z48smKXHjV0zuGfoCfh94J1gQ0ipBNiT9OcaYLCZfQm4AugArEj3YjObAEwA6NlTdy8i0nalk5cHvilnpgl4NfA/gg2h5dx9DbCmBefNA+YBlJeXezOni4iEItNZEHcAvc3sLjPLZlW0WuDkpD/3iLe1mAqyi+SXbOypVugyvQOeBmwjVuvhOjPbBKRaduzuPr4NcW0A+prZKcQS79XANZl04O5VQFV5efkNbYhDRAIyaVi/sENok6H9u8POYPvMNAGPS/q5NP6VigMtSsBmthD4EtDNzGqAae7+kJndAjxDbBraw+6+LcNYRSRkl8xZ21B8p3unDgW7Cg7goXGDQq8HfEqwlwd3H5umfQVNPGhrjpmNAkb16dOntV2ISBsVYtnJdMYv2MBDAfeZaS2IXQFfP2s0BCFH0FzfnJuyeEtBFl9PZfXLdRDyLIiCoTtgOYLm+ubcwvV7IpOAsyHTpcgtvk1w9zS3GrmhO2ARyXeZ3gHv5NPdMZrirehbRCRv7Zw9MvSHcI+QOgF3Bc4CehFbKFEwY8USQdrFOG+smzo07BAC8+i63ZnNhW2BTB/CjUt3LF6h7H8BNwJfb1tYbacx4CKmsd68sbXmAMefEfCTq5BMXbKVawL+qwRWkN3dP3H3e4gNU8wOqt/Wcvcqd5/QpUuXsEMRKSqPrttN6eTllE5ezrSlmr7flGyM074IpNq2XkSKwDWDe8Yqh0mzsrEl0X8Djs1CvxlRLQiRcJROXh52CFkx/9rywPsMNAGb2QXAVUB1kP22hoYgRCRIA3oEn0synQf8XBP9nAwkPnd8ty1BiYjkm8GVq0OvB/ylNO0O7CNWPOcH7p4uUYtIxA3t3z3sEApGptPQtI29iDTpoXGDwg6hYER2tZrmAYvkTvKDt6H9u0cyCY+tOBm2BNtnmxKwmXUitgrugLu/G0xIwVAtCJHcKdSNNjMx64ozA0/AGQ8pmNlRZjbZzF4nthvGTmCfmb0eb4/sXbWIpPboulBrb+XEJXPWBt5nRgnYzI4GngVmEtsNYw+wPv69NN6+Kn6eiBSJqUu2hh1C1iV29ghSpner3yY2E2IZcJu7v5Y4YGa9gX8HRsXPC305skSYiqtLBGSagK8htsjicnf/JPmAu//ZzK4ANgH/iBKwZJMK7kiOde/UAQ4G22emY8B9gKcaJ9+EePtTQO+2BtZWWooskjvZWKabb7KxoWimCfhvwOeaOedYAv89kTktRRbJnWws08039658NfA+M03AW4CvmtlxqQ6aWTfgq8DmtgYmIoVjcOXqsEPIuvtWv9b8SRnKNAHPBY4D1pvZeDM71cw6mtkpZvbPwLr48blBByoi+WfK4i2UTl5OSdeOYYdSkDJdivy4mZ0FTAbmpTjFgP/t7o8HEZyI5K9L5qxl2a3na9fjNsh40YS7TzWzpcB44GygC3AA+E/gYXd/KdgQRSQfZWNebD6ruuU8mB9sn61atebufwT+GGwoIilog02JsGYTcHxV2wvAe8AId085wyF+3lPEZkGcn+48kYxovm9eGTL7Ob46sAeThvWLzYstIqPmvhB4PeCWPIT7n8BA4N+bSqru/jfg+0AFsYUYIhIxtfvrmTSsH5CdebHFpiUJ+ArgDXdf0dyJ7v408Bowpq2BtZUWYohIvmtJAj4bWJNBn88DZ7UqmgBpIYZI8KpuOS/sEEIzcWjfwPtsSQLuBuzNoM+9wOdbF46ISH5KDL0EqSUJuJ7mlx8n+xzwYevCEZF8NmruC2GHEJqKmasC77MlCXgPkEmljXIg+tWZRaSo1L33UeB9tmQe8BrgG2ZW7u4bmzrRzAYCfw/MCSA2KRaq7StFqiUJeC5wE/BrM7vY3XekOsnM+gO/Bj4GHgguRIk8zfXNexUzV7H+zguy8iCqUJSVdIa3g+2z2QTs7q+Y2XeB6cB/mtkTwHNATfyUEmAocCXQAbjb3V8JNkwRCVPi43c2HkQVimW3nh/LggFq0VJkd/+umR0CphHbFWNso1OMWA3gO919VrAhioiEb8riLQSd3FpcC8LdK83sV8B1wBDgxPiht4gtVf65u+8KOD4RyQNlJZ3DDiF0C9fvYVbAS5EzLUe5i9hdsIgUkWW3nh92CJGUaUH20JnZ5Wb2oJk9ZmYXhh2PSDGYsnhL2CFEUl4kYDN72MzqzKy6UfsIM3vFzF43s8kA7v6ku98A3AhcFUa8IsVm4fo9YYcQunVThwbeZ14kYGABMCK5wczaAfcDFwFnAGPN7IykU+6KHxcRybqtNcFPlcyLBOzuzwPvNGquAF539zfipS4XAZdZzL8BT7n7/8l1rCJSnK5/pMl1aK3Sqh0xcqSE2DLohBpgMHArcAHQxcz6uPtPG7/QzCYAEwB69tRKKpHWGDL7OWr31zP/2vKsfPyW/E7AKbn7j4EfN3POPOKbhpaXl3su4pIW0PZCBaV2fz07Z48MO4xIy+cEXAucnPTnHvG2FjGzUcCoPn36BB2XtJaWHBeUytEDwg4hr1SOHhDbdC1AeTEGnMYGoK+ZnRLfb+5qYGlLX6yC7CJtc81gfTJJlo33Iy8SsJktBF4CTjOzGjMb7+6HgFuAZ4AdwOPuvi2DPrUlkUgblE5eHnYIeSUb70deDEG4e+PaEon2FUCze9GleW0VUFVeXn5DW2ITEcmWvLgDFpH8MH7BhobvJV07hhxN9OXFHXA26CGcSOZWv1wHwEPjBoUcSf4Z2r877Ay2z8jeAeshnIgEKRu/lCKbgEVEgpQYnglSZBOwZkGIZE4LL9JLDM8EKbIJWEMQIpl7dJ02NM+lyD6Ek5Boh+OCNnXJVi3AyCElYAmWlhtLRO2cPTLwTTkjOwShMWARCVI2hmcim4A1BiySufnXlocdQt6aumRr4H1qCEKkyCXq/kJ2tt2R9JSARYrcP/Trxqwrzgw7jKKkBCzpNVVAfVLwH8ckHEq+LTP/2nJ4PNg+I5uAVQsiAOlmNEzXuHqUXDJnLctuPT/sMPLegB7B/3evh3AiRa669t2wQygIgytXB95nZBOwiEi+UwIWKXLdO3UIO4SipQQsUuTW33lB2CEUhLEVJzd/UoYi+xBORNJL3t9s4tC+TBrWL8RoCsOsK86ELcH2GdkErFkQWdSlZ/qZECq4k7cSCy6qbjmPqlvOy8pT/Si7ZM5algXcZ2QTsDblzCLNAS5ItfvrVe+3Dapr34XPBtunxoBFisTEoX3DDkEaUQIWKRIa522bbMwWUQIWKRIVM1eFHUJBy8ZsESVgkSJR995HYYdQ0O5d+WrgfSoBi4i0wH2rXwu8TyVgkSJRVtI57BCkESMndu4AABDkSURBVCVgkSKhimf5J7IJWHvCiRxuyuKAl3EVmapbzgu8z8gmYJWjFDncwvV7wg5BGolsAhYRCdKouS8E3qcSsIhISJSARYqEdjzOP5EtxiMtlG7jTVBlswIzZPZz9D+hEw+NG8T4BRtY/XJdw7Gds0eyteYAx58RcDWZIjJxaF/4Q7B9KgEXu3Qbb0rBqd1fzx8mfwWAh8YNOuL49Y9sVDW0Npg0rF/gCVhDECIR9+i63ZROXk5J145hh1LQslFLQ3fAIhGR7u72msE9uWawhpPaqu69jwKvB6wEXCzSjfVqnDcyHl23W4m2wCgBFwuN9Ube1CVblYCzqKykM7wdbJ8FNQZsZqea2UNm9kTYsYhIcclGLY3QE7CZPWxmdWZW3ah9hJm9Ymavm9lkAHd/w93HhxOpiBSzbNTSCD0BAwuAEckNZtYOuB+4CDgDGGtmZ+Q+NJHCMf/a8rBDiLRs1NIIPQG7+/PAO42aK4DX43e8fwMWAZflPDiRAqJt5gtP6Ak4jRIg+ddNDVBiZp83s58CZ5vZlHQvNrMJZrbRzDb+5S9/yXasInlhcOXqsEOQDOVrAk7J3d929xvdvbe7z2rivHnuXu7u5ccdd1wuQxSRiMpGLY18TcC1wMlJf+4Rb2sxFWQXkSBtrQk+l+RrAt4A9DWzU8zsaOBqYGkmHagguxSbsRUnN3+StNr1j2wMvM/QE7CZLQReAk4zsxozG+/uh4BbgGeAHcDj7r4tw351ByxFZdYVZ4YdgmQo9ATs7mPd/UR3b+/uPdz9oXj7CnfvFx/vndmKfnUHLEXlkjlrww5BMhR6AhaRYFTXvht2CJFWOXpA4H1GNgFrCEKKwb0rXwVipRJVbjK7slFnI7IJWEMQUgzuW/0aAOvvvKChGLtkR+nk5YH3GdkELCKS7yKbgDUEISL5LrIJWEMQUgyqbjkv7BCKxtD+3QPvM7IJWEQkSKk2Om0rJWCRAjZq7gthh1A0xi/YEHifkd2SyMxGAaP69OkTdijBS7e/W1O095tIm6x+uU6bcraUu1cBVeXl5TeEHUvgtL+bSCRoCEKkgE0c2jfsEKQNlIBFCtikYf3CDqFo7Jw9MvA+lYBFCljFzFVhh1A0Hl2X4XOXFohsAtZCDCkGde99FHYIRWPqkq2B9xnZBKyFGCKS7yKbgEWKQVlJ57BDkDZQAhYpYMtuPT/sEIrG/GvLA+9TCVikgE1ZvCXsEIrGgB7BD2cqAYsUsIXr94QdQtEYXLk68D4jm4A1C0IK1ZDZzzVML0vseAGwavteSicvZ8js58IKTQIW2QSsWRBSqGr317P+zguAT3e8gNhH4J2zR1K7vz6s0CRgkU3AIoUq3cyGxEfg5M0h100dmpOYBMZWnBx4n0rAInmmuZkNyZtDbq3REFuuzLrizMD7VAIWyTPJMxtS7XiRvDnk9Y9szElMApfMWRt4n0rAInkm3cyGbHwElparrn038D6VgEXyWPKOF9n4CCzhUgIWKRCJj8DJm0MmP5CT7OreqUPgfSoBi+SZdDMbEh+BkzeHTH4gJ9mVmBoYpMhuSVTwe8I1te+b9ncrOIkHZyVdO/KHyV8BYgsuavfXM3Fo38MKq2+tOcDxZ8Q2H0vseFExcxUlXTsCsc0hV79c13B+NgqFy5HuXfkqk4Lu1N0j/TVw4EAvSNM6hx2BZMF1P1/f8HOvO5Y1/DxoxsqU7ZI/et2xrFX/XwIbPU1+0hCESA4kdlNIHj5InmKmwurFSQlYJAcSuymMX7Ah5EgknygBi+RQ8tht8hSz5OXHmtmQn1ItimkrJWCRPJC8/FgzG4qHErBIDjS3m0Ly8uPkpcaSP5I/sQRFCVgkBxK7KSRPGUtMMQMVVi9WSsAiOZAoJZmYDQEcNvdXipMSsEgOJWZDAA27XjSWvNRY8kfyJ5agKAGLhCR57m/y8uPkucKSP7LxiaWgErCZHWtmvzCzB83sH8OOR6SlmislmVxYXXOF81O6TyxtEXoCNrOHzazOzKobtY8ws1fM7HUzmxxvvgJ4wt1vAC7NebAirZQoJZk8GyJ57m9yYfXkucKSP7KxWjH0BAwsAEYkN5hZO+B+4CLgDGCsmZ0B9AASj4s/zmGMIm2SKCWZmA0BzW89JNEXejU0d3/ezEobNVcAr7v7GwBmtgi4DKghloQ30cQvDzObAEwA6NmzFZPam6tENmnrke1NvSaddH1J5CRKSSZmQ0BsWCJxZ1zStSOPrtutRRh5rKykM7wdbJ+hJ+A0Svj0ThdiiXcw8GNgrpmNBKrSvdjd5wHzAMrLyz3jqx/YDdPTbHY4Pc029029Jp10fUlkpSsdmShR2dQ5Eq5lt54P04PtM18TcEru/l/AP4cdh0imsrGbguTWlMVbmBVwn/kwBpxKLZD82LhHvK3FzGyUmc07cEDbdkv4srGbguRWNlYr5msC3gD0NbNTzOxo4GpgaSYduHuVu0/o0kUf8yV89658NewQJA+FnoDNbCHwEnCamdWY2Xh3PwTcAjwD7AAed/dtGfarO2DJG/etfi3sECQPhT4G7O5j07SvAFa0od8qoKq8vPyG1vYhIpKwbupQ+GGwfYZ+BywiUgiSVysGJbIJWEMQkk+ysZuC5FbyasWgRDYB6yGciOS7yCZgkXySjd0UpPBFNgFrCEJEgpSNzVIjm4A1BCEiQcpGnY7IJmCRfJKN3RQkt7KxWaq5Z16rppCY2V+AXWHHkaQb8Newg2iGYgyGYgxGocfYy92PS3Ug8gk435jZRndveo/ykCnGYCjGYEQ5Rg1BiIiERAlYRCQkSsC5Ny/sAFpAMQZDMQYjsjFqDFhEJCS6AxYRCYkSsIhISJSAc8TMvmdmW8xsk5k9a2Ynxdu/ZGYH4u2bzOzuPIzRzOzHZvZ6/Pg5Icb4fTN7OR7HEjPrGm8vNbP6pPfxp/kWY/zYlPj7+IqZDQ8xxjFmts3MPjGz8qT2fHofU8YYP5YX72MyM5tuZrVJ793Fzb7I3fWVgy+gc9LP3wR+Gv/5S8CysONrJsaLgacAA84F1oUY44XAUfGf/w34t/jPpUB12O9hMzGeAWwGOgCnAH8G2oUU4+nAacAaoDypPZ/ex3Qx5s372Cje6cDtmbxGd8A54u7vJv3xWCDvnn42EeNlwCMe80egq5mdmPMAAXd/1mNbVgH8kdiGrXmliRgvAxa5+0fu/n+B14GKkGLc4e6vhHHtlmoixrx5H9tKCTiHzGymme0B/hFIHmr4H2a22cyeMrMvhBQekDbGEiB5S9iaeFvYriN2Z55wipn9p5n93szODyuoRpJjzNf3sbF8fB+T5fP7eEt86OlhM/u75k4OfU+4KDGzVcAJKQ7d6e6/dfc7gTvNbAqxTUenAf+H2Frx9+NjRk8CWavc0soYc6q5GOPn3AkcAn4VP/YW0NPd3zazgcCTZvaFRnf1YceYUy2JMYW8ex/zSVPxAj8Bvkfsk+P3gH8n9gs4LSXgALn7BS089VfENhydlvwftruvMLMHzKybu2el+EhrYgRqgZOTjvWIt2VFczGa2TjgEmCoxwff3P0j4KP4z38ysz8D/YDg95FpZYzk2fuY5jV59T6mkdP3MVlL4zWzB4FlzZ2nIYgcMbPku9rLgJfj7SeYmcV/riD2b/J27iNMHyOwFLg2PhviXOCAu7+V8wABMxsBfAe41N0/SGo/zszaxX8+ldiniDfyKUZi7+PVZtbBzE6Jx7g+jBjTyaf3sQl5+T42ei4yGqhu7jW6A86d2WZ2GvAJsfKYN8bbvwrcZGaHgHrg6qQ7pnyJcQWxmRCvAx8A/xxOeADMJfb0e2X899Yf3f1G4B+A75rZQWLx3+ju7+RTjO6+zcweB7YTG5q42d0/DiNAMxsNzAGOA5ab2SZ3H04evY/pYsyn97GR/21mZxEbgtgJ/EtzL9BSZBGRkGgIQkQkJErAIiIhUQIWEQmJErCISEiUgEVEQqIELCISEiVgEZGQKAFL3ojXonUzWxBgn50tVst4p5n9Ld7/5PixX5hZnZkdG9T1Ulx/evyapU2cc4yZ3W2xGsIfmtkeM6s0s/ZJ5wyM93N9tmKV3NNKOIm6XxGrybAC+CWxlVNLzWwQ8E/E6rf+V1jBxZevriK2nHYJ8Fti8U4BPk98NVW8LsOTwPfMbJG7vx9SyBIgJWCJLDPrTyyZPePuIxsdexZ4l1gFq1CY2dFAFdAL+LK7/yHe/j1gG3C9mU1z9/8Xf8ksYB2xYvmVIYQsAdMQhETZV+Lff5PcaGb9gAuAx929PudRfep2YCBwRyL5AsTvbpcQ+//z/KT29cQKJP2Lmen/3QjQP6LkteRx4fjPi8zsr/Gx0o1mdkmK11xpZg7cH2+aF+/Dzex0YjVaDXgszTWfjZ97ZaN2i8fhZja7jX+vjsC/Equ/Oy/FKYmKeI1rzy4CegLD2nJ9yQ9KwFIoehErOVgK/Aex5FkG/NbMvtzo3LeAe4C/EBvzvSf+NR14ldjd78fEtgtK5V+JVQL7XqI0Y9wPgK8D89x9chv/PqOBrsCj7n4wxfHPxr//rVF74k5ZCTgCNAYsheJLwHR3vyfRYGaPAk8TS5i/S7S7+4tmto5YTd4d7j496TXHAmfF21M+fHP3zWb2H8SS7T8BC8xsKvBt4HHgpgD+Pokx6RIzm57ieKLw955G7Rvi3/8hgBgkZErAUih2ATOSG9z9GTPbTeoNGc8AOhLb8ilZCdCO2F1yU/4XcBUwzcw+B8wEngH+yd0/yTz8I5wX/351M+dtT/6Dux8wsw+JDUNIgdMQhBSKTWmKbu8BUm1+eE78e+ME/Pn4931NXczd9wA/IjbkMQd4EbjC3RsPCWQsfhfeE9jm7tb4C+gMHAT2uPvOFF28A3RraxwSPiVgKRT707QfIvV/x4kE/J+N2hOzHj5L8/6S9PP4RtsLtUViB990+5hdCLQnNnc5lY58+veQAqYELFF1DrGtYTY1aq+Lf/88TTCza4g9dEvMwZ2Y4hwzs++Y2StmVh9fVfebxuelcHT8+0dpjie2fHo4xTU/Q+zhXV3jY1J4lIAlcuKbnP534DV3f6/R4beI3dme1sTrLwYWENtU8UzgFWKLIhq/5l+BccA3gP7ApcDKFoSYSOpHbG8e3/T0YuCp+Lzfxk4jNoWu8S8WKUBKwBJF/YBOHDn8QHzD0+eBbmbWp/FxMzsPeAKoAYa7+1+Au4g9sP63RqePIJYoV7v7Lnf/o7v/tLng3P2vwA5goJmdmXTtXsBC4ACxpJ7KufHvv0tzXAqIErBEUboHcAmJYYLhyY3xHW2XEUuAw9z9LQB3fwLYCFxmZucnvWQp8C0zW2Vm/2JmmTwYm0Hs/7/VZnavmc0DNhMbXhiZ5uEbxMaHPyZWM0IKnBKwRFFLEnAdcG2iIX43/DSxcePh7v7nRq+ZEv/+/USDu/+I2JDA08TuWP8cX2nXLHd/lNjwRR2xecUXE5tjXObuL6Z6jZl1AS4HlsVnaUiB07b0UpTMbAqxgjbnuPsRQxWt6O8oYtPDbnD3x5LapwPTgFOauKtt6TVuBX4MnO/uL7SlL8kPugOWYnUvsBv4bmtebGZ3mNk4MzsjXtznHmLLhtcEF+Jh1+tI7C78N0q+0aGVcFKU3P1DM/sn4MtmdmwragJ3AO4gVqPiA2J1JYa6+96AQ00oJVa0Z0GW+pcQaAhCJIuCHIKQ6NEdsEh2rYl/T7eST4qY7oBFREKih3AiIiFRAhYRCYkSsIhISJSARURCogQsIhISJWARkZD8f9Z9GN5UDk2ZAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "ax.hist(\n",
    "    signal_pdf.log_prob(sample),\n",
    "    bins=40,\n",
    "    histtype=\"step\",\n",
    "    linestyle=\"dashed\",\n",
    "    label=r\"$f\\,(x_s|1)$\",\n",
    ")\n",
    "ax.hist(\n",
    "    bkg_pdf.log_prob(sample),\n",
    "    bins=40,\n",
    "    histtype=\"step\",\n",
    "    label=r\"$f\\,(x_s|0)$\",\n",
    ")\n",
    "\n",
    "ax.semilogy()\n",
    "ax.set_xlabel(r\"$\\ln f(x_s|\\theta)$\", fontsize=20)\n",
    "ax.set_ylabel(\"Count\", fontsize=20)\n",
    "ax.legend(loc=\"best\", fontsize=14);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Now, let's recreate Figure 5b of \"Asymptotic formulae for likelihood-based tests of new physics\" ([arXiv:1007.1727](https://arxiv.org/abs/1007.1727)). The pseudoexperiment are still fast, but the test statistic evaluation is slower.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = pyhf.simplemodels.hepdata_like([6], [9], [np.sqrt(9)])\n",
    "\n",
    "signal_pars = model.config.suggested_init()\n",
    "signal_pars[model.config.poi_index] = 1.0\n",
    "\n",
    "bkg_pars = model.config.suggested_init()\n",
    "bkg_pars[model.config.poi_index] = 0.0\n",
    "\n",
    "par_bounds = model.config.suggested_bounds()\n",
    "\n",
    "signal_pdf = model.make_pdf(signal_pars)\n",
    "bkg_pdf = model.make_pdf(bkg_pars)\n",
    "\n",
    "# protect against limited resources on Binder\n",
    "sample_size = 5000 if getpass.getuser() == \"jovyan\" else 10000\n",
    "sample_shape = (sample_size,)\n",
    "\n",
    "signal_sample = signal_pdf.sample(sample_shape)\n",
    "bkg_sample = bkg_pdf.sample(sample_shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "signal_qtilde = np.asarray(\n",
    "    [\n",
    "        pyhf.utils.qmu(1.0, sample, model, signal_pars, par_bounds)\n",
    "        for sample in signal_sample\n",
    "    ]\n",
    ")\n",
    "bkg_qtilde = np.asarray(\n",
    "    [pyhf.utils.qmu(1.0, sample, model, bkg_pars, par_bounds) for sample in bkg_sample]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAFgCAYAAACFYaNMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3de3yU9Zn38c9lAIMIwYqgBCFFwEMTKiUGumh1CwgVoYvVrYcusuXwUEVT+mzbgPsI7aPgtlatWLsvC4o8CrZaKCC4Fei6graclEMECooIASRVKXhABLmeP2YSkpCQmWTmvmcm3/frlReZe+6Z3xUPX+785ve7bnN3REQkeKeFXYCISFOlABYRCYkCWEQkJApgEZGQKIBFRELSLOwCUoGZDQWGtm7dekyPHj3CLkdE0ty6devec/dz6jvPtAzthMLCQl+7dm3YZYhImjOzde5eWN95moIQEQmJAlhEJCSaA67ig48/AyCvZHHlsf4XtWfmyMsYNWsNy7eWVx7fed8Q5qzaxc19Ogdep4hkBs0BV3H6ed39yL7tMZ+fV7KYnfcNSWJFIpKONAcsIpLiNAUhkoYOHTpEeXk5R48eDbuUJqlVq1Z06tSJ005r3DWsAriKLmefEdf5M0bU+xuGSMIdOnSI/fv3k5ubS8uWLTGzsEtqUo4fP86ePXt47733aN++faPeS1MQVbRsHt/fRwWdcpJUiUjdysvLyc3N5YwzzlD4huC0006jQ4cOHDx4sPHvlYB6MsbWdw/FdX6fqcuTVIlI3Y4ePUrLli3DLqNJa968OceOHWv0+yiARdKQrnzDlah//gpgEZGQKICr+EKrFnGdf1PR+UmqRESaAgVwFblt45tXm3ZdzyRVIiJVLVq0iEsuuYRbb72Vp59+mhtvvDHskhJCAVzFm+UfxXX+tdNXJKkSkcz2m9/8hq5du9KsWTPGjh1befzAgQN06NCBt956q9r5c+bMYcGCBXTt2pWJEycyYsSImMZ56aWXyMvLi6u2G264gV/84hdxvaahtA64isNHP4/r/NI98a2aEEmWfvf9iT1/P0z71qez+q4BPLh0G79cfmJb/aLxlwMw9JGVlceK+3dnwsAeFN27jPIPjwCQn9uG5++4gonzNjJ39e7Kc1dN6s+msoOMnh1p15rbtiWvlHy9QbVu3bqV733vezz77LP07duX1q1bVz43depUrrnmGi644IJqr5k7dy4AkydPZvLkyQ0at8LLL7/M/fffz7p169i7dy9PPPEEI0eOrHz+7rvv5sorr2T06NHk5CR3qakCWCQDnNWqebVAnDCwBxMGnnxzgdp6l6y+a8BJx6Zd1/OkKbYOl2RXvn7Oql0NrnXhwoXk5+czfPjwasc/+eQTZsyYwaJFixr83rH46KOPyM/PZ8SIEbVeSRcUFNC1a1eeeuopbr/99qTWoimIKpqdFt/SkvatT09SJSLxef6OKwIdr6FdAHv06MGPf/xjNmzYgJlVC+ElS5ZgZvTr16/aa7Zt28bAgQPJzs7mggsu4IUXXiA7O5tly5Y1qIZrrrmGqVOncv3119e5lXjYsGGVV93JpACu4uLz2sR1fm1XDiJhmDhvY6DjVW3ZGo+VK1fSo0cP7rnnHvbt28eTTz5Z+dyKFSvo3bt3tTW227dvp6ioiMLCQkpLS3n44YcZPXo0R44c4dJLL230z1GXoqIiVq9ezeHDh5M2BiiAq9l/6NO4zn9w6bYkVSISn6rztamsTZs27Nixg379+nHuuefSps2Ji5533nmHjh07Vjt//PjxDBs2jGnTptGtWzeGDBnC4MGDyc3NpV27dkmrs2PHjhw9epS9e/cmbQxQAFdT8UFErKp+yCEi9SstLeXYsWO1Xr0ePnyY7Ozsyse7d+/mxRdf5Pvf/36181q0aMGXv/zlysfDhw/nrLPO4vrrr09YnRVbvXUFLCIpp/9FDesCtn79erp06ULbtm1Peq5du3YcOHCg8vHrr79OVlYW+fn51c7buHFjtQAvLi5m9uzZDaqnLh988AEA55xT742NG0UBLJIBVk3qH+h4M0de1qDXrV+/vs652169erF58+bKx2bG8ePHq/U8fuWVV3j11VervcdVV11VbSlbIpSWlpKbm0uHDh0S+r41ZWwAm1krM3vSzH5jZrfE8ppu55wZ1xgVaytFwraprPGtEeMxataaBr3uVAE8aNAgtmzZwvvvvw9A7969adGiBSUlJezYsYMFCxYwatQogGpTEPH66KOPWL9+PevXr+f48ePs2rWL9evXs2vXiaV1K1asYNCgQQ0eI1ZpFcBm9riZlZtZaY3jg83sr2b2ppmVRA9fBzzn7mOAYYEXKxKgig0SQal6g9pYuftJ0wdVFRQUUFRUxDPPPANEPgibOXMmCxYsoGfPnsyaNYvbbruNM888k27dujW49rVr19KrVy969erF4cOHmTx5Mr169eLuu+8G4NNPP2X+/PmMGTOmwWPEKt02YswCHgEqJ3zMLAv4FTAQKAPWmNlCoBOwKXpaTFvc3vxbfFuRhz6yUjfllJRRdWlYone+AUwdXsDNfTqTV7I47r4pEJlSOHTo1LtHJ0+eTHFxMePGjSMrK4tbbrmFW2458QvslClTKCgoaNStgK666ipOdTPimTNn0qdPH/r27dvgMWKVVgHs7i+bWV6Nw0XAm+6+A8DMngG+SSSMOwHrOcWVvpmNBcYCZLVJ7oS7SLLUdSGQqJ1vsYyVCIMHD+b222+nrKyMLl26nPT8xo0bGzX9EIvmzZszffr0pI5RIa2mIOqQC1RdBFkWPTYP+JaZ/Rqoc2+juz/m7oXuXph1hm4xJBK2O++8s9bwBdiwYcNJUxgDBgzghhtuYMmSJXTq1Ik///nPjRp/7NixXHjhhY16j1il1RVwPNz9Y+Bf43lNvFuLi/t3j+t8EWmcml3SgFNuSc7LyztpHXEqyYQr4D1A1c7onaLH4tahTXb9J1VRW7MTEUkdCuDkWwN0N7MvmlkL4EZgYTxvYGZDzeyx0rf3xTVw0b0NawYiIgJpFsBmNhf4M3ChmZWZ2Sh3PwaMB/4IbAF+5+5vxPO+7r7I3cd6izPiqifercsiIlWl1Rywu99Ux/ElwJKAyxERaZS0ugJOtpbNs+I6Pz83vvaVIiJVKYCr6NY+vq3IQTfBFpHMogDmxIdwb+2Jb3tl0E2wRSSzKIA58SHcJ8S3DjhdmmCLSGpSAIuIhEQBLCISEgVwFRedG9+qhqCbYIs0VYsWLeKSSy7h1ltv5emnn+bGG28Mu6SEUABz4kO4997/IK7XBd0EWyRT/OY3v6Fr1640a9aMsWPHVh4/cOAAHTp0OKnnw5w5c1iwYAFdu3Zl4sSJjBgxIqZxXnrpJfLy8uKq7YYbbuAXv/hFXK9pKDtVX8ym5vTzuvuRfbHfaDOvZLH6AUvgtmzZwsUXX1z94IMFcHBX7S9IhpzOMGFT/efVYuvWreTn5/Pss8/St29fWrduzZlnRpaA/vCHP+S9997jiSeeSEiZL730EiNHjmTnzp3Vjj/66KP8/Oc/Z9++fXzpS1/ioYce4oorIstKN23axJVXXsnbb79NTk7dHRJr/fcQZWbr3L2wvvrSaieciNTh4C6YEuBvZFMa3rp14cKF5OfnM3z48GrHP/nkE2bMmMGiRXV2j02I3/72txQXF/Poo49y+eWX8+ijj/KNb3yDzZs307lzZwoKCujatStPPfUUt99+e1Jr0RSEiASmR48e/PjHP2bDhg2YWbUQXrJkCWZGv379qr1m27ZtDBw4kOzsbC644AJeeOEFsrOzT9mG8lQeeOABRo4cyZgxY7j44ouZPn065513Hr/+9a8rzxk2bBhz585t2A8ZBwVwFfHeZmXq8IIkVSKSmVauXEmPHj2455572LdvH08++WTlcytWrKB3796YWeWx7du3U1RURGFhIaWlpTz88MOMHj2aI0eO1HlvuVP57LPPWLduHVdffXW141dffTWvvvpq5eOioiJWr17N4cOHG/BTxk4BXMUXWrWI6/yb+3ROUiUimalNmzbs2LGDfv36ce6559KmzYmVR++88w4dO3asdv748eMZNmwY06ZNo1u3bgwZMoTBgweTm5tLu3bt4h7/vffe4/PPPz/pdvMdOnTg3XffrXzcsWNHjh49yt69e+MeIx4KYE6sgtjwVnz/sKveBFFE6ldaWsqxY8dqvXo9fPgw2dknboqwe/duXnzxxZMaqrdo0aLafeGGDx/OWWedxfXXX5+wOlu2bFlZUzIpgDmxFfm07FZhlyKS0davX0+XLl1o27btSc+1a9eOAwcOVD5+/fXXycrKIj8/v9p5NW9tX1xczOzZs4lFu3btyMrKYv/+/dWO79+/n3PPPbfy8QcfRJaknnNOcm/UqwAWkcCsX7++zrnbXr16sXnz5srHZsbx48c5evRo5bFXXnmFV199tdp7XHXVVbRu3Tqm8Vu0aEHv3r1ZunRpteNLly7lH/7hHyofl5aWkpube9JURaIpgKtonR3fqrz+F7VPUiUimelUATxo0CC2bNnC+++/D0Dv3r1p0aIFJSUl7NixgwULFjBq1CiARt2a/gc/+AGzZs1ixowZbNmyheLiYvbu3cu4ceMqz1mxYgWDBg1q8BixUgBXkXd2fFMQM0delqRKRDKPu580fVBVQUEBRUVFPPPMM0Dkg7CZM2eyYMECevbsyaxZs7jttts488wz6datW4Pr+Pa3v81DDz3EPffcw6WXXsrKlStZsmQJXbp0AeDTTz9l/vz5jBkzpsFjxEobMarY+f7HcZ0/atYahbCkhpzOjdoc0aDx4mRmHDp06JTnTJ48meLiYsaNG0dWVha33HILt9xyS+XzU6ZMoaCggNNOa9y142233cZtt91W63MzZ86kT58+9O3bt1FjxEIBXMWHnx6L6/zlW+Nr4C6SNA3cFpxqBg8ezO23305ZWVnlFWlVGzdubNT0QyyaN2/O9OnTkzpGBQUwkWVowNAuZzWL6ypi5entAPWCEEmkO++8s87nNmzYwI9+9KNqxwYMGMCGDRv4+OOP6dSpE88++yxf/epXGzx+1eZAyaYAJrIMDVhU2DFrTDz76TsF+SufiJzUJQ045ZbkvLy8k9YRpxJ9CCciGUsBLCIitVIAi4iERAEsIhISBbBIGtKdbMKVqH/+CmCRNNO8efOkd+mSUzt69CjNmjV+EZkCmBPtKA9+qqsKSX3t27dnz549fPLJJ7oSDsHx48fZv3//Ke8XFyutA6bGOmCRFFfRxHzv3r3VOoVJcFq1atWghvA1KYBF0lCbNm2q3U1C0pOmIEREQqIAFhEJiQJYRCQkCmARkZAogEVEQqIAFhEJiQJYRCQkCmARkZAogEVEQqIARr0gRCQcCmAivSDcfWxOtoVdiog0IQpgEZGQKIBFREKiABYRCYkCWEQkJApgEZGQKIBFREKiABYRCYkCWEQkJApgEZGQKIBFREKiuyI3Qpm3o9OUnNhfkNMZJmxKXkEiklYUwI0wrt0snr/jithfEE9Yi0jGUwA3wvN3XMHEeRuZu3p35bFVk/qzqewgo2evrTw2dXgBN/fpHEaJIpLCzF0tGM1sKDC021k2ZvsHx5M30JQcmHIwee8vIinBzNa5e2F95+lDONSOUkTCoSmIVPdgARzcFfv5+qBPJG0ogFPdwV3xTVvogz6RtKEpCBGRkOgKOEBxrxuGyJSCiGQkBXCALj/yMDvvGxJ2GSKSIjQFISISEgVwgHT1KyJVaQoiQHNW7Ur/HXHxLosDLY0TqYMCOECT5m9K/wCOd1kcaGmcSB00BSEiEhJdAUvy5XRO7lWwpjgkTSmAAzRjRL29OTJTssPxwYL4Al6BLSlCARyggk6aC02KeMNUc9KSIjQHHKA+U5eHXYKIpBBdATd1Dem2JiIJoQBu6hqyrExEEkJTEAG6qej8sEsQkRSiAA7QtOt6hl2CiKQQTUEE6NrpK+K7i3JDxLvmVnO6IqFRAAeodM+h5A+i9a0iaUNTECIiIdEVcIDatz4dgAeXbuOXy7dXHl80/nIAhj6ysvJYcf/uTBjYg6J7l7H6rgHBFioigcjoADazrsBdQI67Xx92PRVBOmFgDyYM7HHS87X1Cy7/8EjS6xKRcKTsFISZPW5m5WZWWuP4YDP7q5m9aWYlp3oPd9/h7qOSW6mISMOk8hXwLOARYHbFATPLAn4FDATKgDVmthDIAqbVeP133b08mFKTJz+3TdgliEiSpGwAu/vLZpZX43AR8Ka77wAws2eAb7r7NODahoxjZmOBsQCdc6zB9SZL0petiUho4p6CMLOBZvYzM/uLme01s8/M7KCZbTezZ83se2aWm4xigVxgd5XHZdFjddV6tpn9J9DLzCbWdo67P+buhe5eeM4ZqRfAE+dtDLsEEUmSmALYzM4wsxIzexv4L+DfiFyNtgXKgWNAV+BbRKYI3jaz35vZV5NTdmzc/X13H+fuF0SvktPO3NW76z9JRNJSvVMQZvZd4P8C5wFbgZ8ArwBr3P1QlfMMuBDoCwwCvgn8k5k9B/zQ3eO8k2Ot9gBVGyp0ih4TSa6GdI3TphipRyxzwDOAPwDT3H1NXSe5uxMJ6K3ALDNrA9wKlAAjgZ82ulpYA3Q3sy8SCd4bgZsb+6ZmNhQY2u2s1JuCkBQRb9c4NX2XGFgkN09xgtlX3P21Bg9glg3kufvWOF83F7gKaAfsBya7+0wzuwZ4iMjKh8fd/d6G1lZTYccsX7v380S9XULsP/QpHdpkh11GZon3ahbiv6LVFXOTZmbr3L3ee5DVG8BNSSoG8LLN+xlwSYewy5Bkm5KjvswZJNYAjnkZmpm1B64AzgU+I7IaYXOC5nalDpMXvsHo2WsBmDq8gJv7dCavZHHl8/0vas/MkZcxatYalm8tJ7dtS14p+XpY5YpIHGL5EK4ZMB0YTS2rJsxsN7AYeMLd1ya8wiautjCtbcvyzJGXAVQLZxFJbbFcAf8E+F/ALmA+8D6QDQwgshTtPOB7wDgzWwLcnm5XxZn0IVxt4SwiqSmWdcD/ArwBfMndJ7j7Pe7+78ALgANnAzdEH18DvGZmfZJVcDK4+yJ3H5uTnf4BPGdVWv3dJ9KkxRLA5wDPu/vHtT3p7h+5++/d/VqgH/AhsNjMzktgnRKjSfP1SbpIuoglgHcCebG8mbv/hcjSsSzg7oYWJSLSFMQSwHOB4WbWN5Y3dPd3gN8BmowUETmFWAL458BbwItmNi7aErI+HxPZQJEWzGyomT128NP0XxM9Y0S9Sw9FJEXUG8DufhjoD2wj0mhnp5n9jEjfh5OYWT7wHSKhnRYy6UO4gk7aAiuSLmLqhubu7wJfBf4P0IpIN7RvA5jZRjN7wcz+YGZ/AV4ncvX7s+SULKfSZ+rysEsQkRjFvBPO3Y8CU83sAeA6IkvOLgfyo18VtgD3uPvcRBYqIo2k/hQpJ+47Yrj7p8Cc6BdmdgaRpugtgHfd/f2EVigiiaGObimn0bckcvdPgO31niiBuKno/PpPEpGUkLJ3RQ5SJq2CmHZdz7BLEJEYxdKMp3MCx/t71btopAp3XwQsKuyYNSbsWhrr2ukrKN0T+UfcvvXprL5rAA8u3cYvl5/4JWXR+MsBGPrIyspjxf27M2FgD4ruXUb5h0eAyB2Zn7/jCibO21jt1kjqNyGSGLE0ZD9OpOdDIvzE3RNxZ4ykSMV+wKlGDeKTJIh+wPGOoR7FDZbIfsA/JXEB/D8Jeh8Jyaayg3S4RAEskgj1BrC7TwmgDkkT8TaIr6BpC5GTNXoVhDQt8TaIF5G6aRWEJJ16FIvUTlfAZNYdMVLRpPmbuLlPIhfTZKCczvFtfNAutYyQ9AA2s/OJrH74brLHaqhMWoYmaSreMNUutYwQxBTEF4BbAxhHRCStNPoK2MxG1HOKfvds4tSjOE1pWiTpEjEFMQv4hLrXCuuDviZOPYqTIN5wrHhNPDQtknSJCOC9wJ3uPq+2J83sUmBdAsaRNNVn6nKtA040XWlmhERcna4DvnKK5x3Q8gIRkRoScQV8P3DmKZ5/E/jHBIwjIpJREtEPeEU9z39MiveA0Drg5KroURxPpzbNG0tTUGc3NDPrCJzm7mU1jr8BvBb9eh143d0zomWSuqGlhrySxZozTkfqnlYpEd3QxgB3m9kcd/+XKsdzgYuBm6sMtpPqofyau5cjIiJ1OlUA/xT4GPgPM5vi7hW3mT8PeA74BrADeJ9IKH+LyM06ATCzdzk5lNUUQEQkqs4Adnc3s68C+4ksNavwM+BK4Mqq879mlg/cDVwPbAY6AEOiXwDHTzWeSIXi/t3DLkEkEHUGopkNBgYDV7n74SpPXQ88U/PDN3cvBf7ZzO4GvgNcALQBehFZpnZpgmuXDDVhYI+wSxAJxKnWAV8EvAVsq3E8B3ivrhdFbzn0ETDJ3cvcfZG7/8Tdhze6WmkSiu5dFnYJIoE4VQA/C2QDH5hZ1Y0W26h/Xe8y4IZG1iZNVMVNQUUyXZ0B7O57gAuBnsBfqzz1OHCZmU04xft+AeiYkApFRDLUKbcie8Qb0c0UFX5F5Ar3fjN7psbVMWZ2JXAT1T+4E4lZfm6bsEsQCUTcqxLc/XMzG0IkiEcDN5jZ+8Au4Gwi7ScNmJ7IQqXpeP6OK8IuQSQQDWrG4+5H3X0s0Ad4CvicyEqHTkSmK77r7g8lrMokM7OhZvbYwU/r6qgpQZo4b2PYJYgEolHd0Nx9jbvf6u7nAS2Blu5+ibvPSkh1AYmu1Bibk61eEKlg7urdYZcgEoiEbYxwd310LSISh3qvgM2sZWMHScR7iIhkmlimIN42s2IzOz3eNzezL5vZAuDf4i9NmqpVk/qHXYJIIGKZgvgj8AAw2cx+C/wO+EuN7cmVzKwrMAgYARQBu4GfJ6ZcaQo2lR2kwyXZYZchyfZgARyMsz9Xht34s94AdvdbzewR4F5gbPTrczPbAuwDDhDZMXc2kY0b7YgsQysH7gIe1PywxGP07LXqB9wUHNwVf//gDLvxZ0wfwrn7GuBqM+sOjAL6E2muU1Dj1L8B84DfA79396MJrFVEJKPEtQrC3bcDJQBmdgaRPsBnA4eBcnffl/AKRUQyVIOXobn7J8D26JdIwkwdXvMXK5HMFMsytB+Z2UVBFCMCcHOfzmGXIBKIWJah3Qf8c8UDM7vQzLRlTJImr2Rx2CWIBCKWAD5K9amKzURWN4iISCPEMge8l8hthSoYjewhISIZKKdzfMvEcjTVFEsALwLGm9liIkvMANQ2TJKm/0Xtwy5BGiKDNkgEJZYAvgv4IpG7Gw+OHvt3M/sW1W87v75G43aRBpk58rKwSxAJRL1TCe7+obsPBS4G7oge/jvQHRgJPAy8DBw0sy1mNsfM/s3Mvp6kmhNO/YBTy6hZa8IuQSQQMc/luvtf3f3R6MNfAa2J3C9uJJG7X7xK5D5wNwI/A5YmtNIkUj/g1LJ8a3nYJYgEoiEbMa4H9rr7caA0+jW74kkz6wH0pvoHdyIiUkND7gk3r57ntxG5df3chhYlItIUaDmZpBx1QpOmImG3JBJJlDmrdjFp/oklTTNGFFLQKYc+U5dXHrup6HymXdeTa6ev4MDHR3mlJG0+8xWpZO765L9CYccsX7v387DLkDjllSzWVXNTMSUn/h7CITCzde5eWN95moKQtNe+ddx3yxJJCQpgSXur7xoQdgkiDaIAlrT34NJtYZcg0iAKYEl7v1yuewJIelIAi4iERAEsIhISBbCkvUXjLw+7BJEGUQCLiIREASxpb+gjK8MuQaRBtBVZ0l5u25bklSymuH93JgzsQdG9yyj/8AgA+blteP6OK5g4byNzV++ufI12zkkq0FbkKrQVuWnYf+hTOrTJDrsMaQhtRRZJb5vKUv9/YGkaFMDS5IyevTbsEkQABbCISGgUwCIiIdEqCGlypg4vCLsEaaiczpEP4uI5f8Km+s8LSUYHsJn9EzAEaAPMdPcXQy5JUsDNfTqHXYI0VLxhGk9YhyBlA9jMHgeuBcrdPb/K8cHAL4EsYIa731fXe7j7H4A/mNlZwP2AAlh0B42mJMWvmFM2gIFZwCNUv+V9FvArYCBQBqwxs4VEwnhajdd/193Lo9//e/R1ItKUpPgVc8oGsLu/bGZ5NQ4XAW+6+w4AM3sG+Ka7TyNytVyNmRlwH/CCu7+W3IpFROKTbqsgcoHdVR6XRY/V5Q5gAHC9mY2r7QQzG2tma81s7d8+0a7ApqD/Re3DLkEESOEr4ERw94eBh+s55zHgMYhsRQ6iLgnXzJGXMWrWGpZvLa88tvO+IcxZtYtJ80/8yjpjRCEFnXLoM3U5uW1b8krJ18MoVzJYugXwHuD8Ko87RY+JxGXmyMtOOnZzn861rpDYed8Q8koWB1GWNDHpNgWxBuhuZl80sxbAjcDCkGuSJuCmovPrP0kkTikbwGY2F/gzcKGZlZnZKHc/BowH/ghsAX7n7m8kYKyhZvbYwU81AyG1m3Zdz7BLkAyUsgHs7je5+3nu3tzdO7n7zOjxJe7ew90vcPd7EzTWIncfm5NtiXg7yUDXTl8RdgmSgVI2gEVSSemeQ2GXIBlIASwiEhIFMJoDlvq1b3162CVIBlIAozlgqd/quwaEXYJkIAWwSAweXLot7BIkAymARWLwy+Xbwy5BMpACWEQkJApgEZGQKIDRKgip36Lxl4ddgmSgdGvGkxTuvghYVNgxa0zYtUjqqtqQp7h/dyYM7EHRvcso//AIAPm5bXj+jiuYOG8jc1ef6Jq6alJ/NpUdZPTstZXHpg4v4OY+nSvfU93WmiZz11VfhcKOWb527+dhlyFN0KhZa2rt0CYBe7AADu6K/fw6bmFkZuvcvbC+l+sKWCQFKHxTRMC3MNIcsEgKGDVrTdglSAgUwOhDOAlf1btzSNOhAEZbkUUkHApgEZGQKIBFUsDO+4aEXYKEQAEskgLmrIpj6ZNkDAWwSAqYND/O5U+SEbQOWCQF5LZtSV7JYmaMKKSgUw59pi6vfO6movOZdl1Prp2+ovLWSO1bn64exRlAO+Gq0E44EYlLHTvn7CeHYtoJpwAmsg4YGNrtLBuz/YPjYRg8aNMAAAqDSURBVJcjUq8Hl25jwsAeYZchdYh1K7LmgNE6YEk/ahCfGRTAIiIhUQCLiIREASyShtQgPjMogEVEQqIAFklDQx9ZGXYJkgAKYBGRkCiAUT9gEQmHAhitA5b0U9y/e9glSAIogEXSkHbBZQYFsEgaKrp3WdglSAIogEXSUPmHR8IuQRJAASwiEhIFsEgays9tE3YJkgAKYJE09PwdV4RdgiSAAlgkDU2ctzHsEiQBFMAiaWju6t1hlyAJoHvCiTQR/e77E1/r0a7O+8s9uHRbtUbvi8ZfTkGnnLDKbRIUwFS7JVHYpYgkzZ6/H2badT2B2ueQJwzsUW2Dx6ayg4HV1lRpCgJtRZb0s2pS/6SPoY5ryacAFklDDbk6DSK0JT4KYJE0NHr22rhfoymF1KMAFmki4g1tdVxLPgWwiNRKHdeST6sgRNLQ1OEFAOSVLK481v+i9swceRmjZq1h+dbyyuM77xvCnFW7yG3bMq4xiu5dxuq7BiSmYKmVuesuEBUKO2b52r2fh12GSErIK1nMzvuGhF1GWjKzde5eWN95moIQEQmJAlhEaqWOa8mnABaRWqnjWvIpgEWkVuq4lnwKYBGplTquJZ8CWEQkJApgEZGQaCOGiNRq1aT+LNu8v9oW5qnDC7i5T+c6N4DMHHlZGKWmLQUw6gcsUpsObbLpcEl2rZsxajtWdfedxEZTEKgfsIiEQwEsIhISBbCIJIT6RsRPASwiCTFn1a6wS0g7CmARSYhJ8zeFXULaUQCLiIREASwiEhIFsIgkxIwR9fYflxoUwCKSEAWdcsIuIe0ogEUkIfpMXR52CWlHASwiEhIFsIhISBTAIpIQNxWdH3YJaUcBLCIJMe26nmGXkHYUwCKSENdOXxF2CWlHASwiCVG651DYJaQdBbCISEgUwCKSEO1bnx52CWlHASwiCbH6rgFhl5B2FMAikhAPLt0WdglpRwEsIgnxy+Xbwy4h7SiARURCkrEBbGYXm9l/mtlzZva9sOsREakpJQPYzB43s3IzK61xfLCZ/dXM3jSzklO9h7tvcfdxwD8D/ZJZr4jAovGXh11C2mkWdgF1mAU8AsyuOGBmWcCvgIFAGbDGzBYCWcC0Gq//rruXm9kw4HvA/wuiaJGmLq9kceX3xf27M2FgD4ruXUb5h0cAyM9tw/N3XMHEeRt5edt7vFLy9bBKTQnm7mHXUCszywOed/f86OOvAlPcfVD08UQAd68ZvrW912J3r/ee2YUds3zt3s8bU7aIxCivZHHG3srezNa5e723CEnVK+Da5AK7qzwuA/rUdbKZXQVcB5wOLDnFeWOBsdGHR2pOewSoHfBeExo3zLGb4s8c5th1jmv/Ed7YSXZhLCelUwDHxd1fAl6K4bzHgMcAzGxtLH9rJUNYY+tn1tiZOm6YY5vZ2ljOS8kP4eqwB6jacLRT9JiISFpKpwBeA3Q3sy+aWQvgRmBhyDWJiDRYSgawmc0F/gxcaGZlZjbK3Y8B44E/AluA37n7Gwke+rEEv186jK2fWWNn6rhhjh3TuCm7CkJEJNOl5BWwiEhToACOimeXXYLHrXXXXwDjnm9m/21mm83sDTMrDmjcbDNbbWYbouP+JIhxa9SQZWavm9nzAY+708w2mdn6WD8lT9C4baNb8rea2Zbomvogxr0w+rNWfB0ys+8HNPaE6H9fpWY218yygxg3OnZxdNw36v153b3JfxHZTfcW0BVoAWwALglo7K8BXwFKA/6ZzwO+Ev2+NbAtiJ8ZMODM6PfNgVVA34B/9h8Ac4hs9Aly3J1AuyDHjI77JDA6+n0LoG0INWQB7wJdAhgrF3gbaBl9/DtgZEA/Zz5QCpxBZJnvMqBbXefrCjiiCHjT3Xe4+2fAM8A3gxjY3V8GPghirBrj7nP316Lff0jkg83cAMZ1d/8o+rB59CuwDyLMrBMwBJgR1JhhMrMcIn/JzwRw98/c/e8hlNIfeMvd3wlovGZASzNrRiQM9wY07sXAKnf/xCMLB/6HyIawWimAI2rbZZf0MEoV0W3fvYhcjQYxXpaZrQfKgaXuHsi4UQ8BPwKOBzhmBQdeNLN10R2YQfgi8Dfgiei0ywwzaxXQ2FXdCMwNYiB33wPcD+wC9gEH3f3FIMYmcvV7hZmdbWZnANdQff9CNQrgJs7MzgR+D3zf3QO5ra27f+7ulxLZTFNkZvlBjGtm1wLl7r4uiPFqcbm7fwX4BnC7mX0tgDGbEZni+rW79wI+BgL7jAMgum5/GPBsQOOdReQ32C8CHYFWZvadIMZ29y3AfwAvAv8FrAfqbDCjAI5okrvszKw5kfB92t3nBT1+9Ffh/wYGBzRkP2CYme0kMs30dTN7KqCxK67McPdyYD6Rqa9kKwPKqvyW8RyRQA7SN4DX3H1/QOMNAN5297+5+1FgHvAPAY2Nu890997u/jXgAJHPV2qlAI5ocrvszMyIzAtucfcHAhz3HDNrG/2+JZH2oluDGNvdJ7p7J3fPI/Lv+E/uHsiVkZm1MrPWFd8DVxP5dTWp3P1dYLeZVTSH6Q9sTva4NdxEQNMPUbuAvmZ2RvS/8/5EPuMIhJm1j/7Zmcj875y6zs3YZjzxcPdjZlaxyy4LeNwTv8uuVtFdf1cB7cysDJjs7jMDGLof8C/Apuh8LMAkd6+zc1yCnAc8Ge3vfBqRHY2BLgcLSQdgfiQPaAbMcff/CmjsO4CnoxcXO4B/DWjcir9sBgL/K6gx3X2VmT0HvAYcA14n2B1xvzezs4GjwO2n+tBTO+FEREKiKQgRkZAogEVEQqIAFhEJiQJYRCQkCmARkZAogEVEQqIAFqmDmXWKtjX8RvTxaWY2JdrQR6TRtA5YpBZmdhqRHXqfARcALwCvAj8n0lJxV4jlSYbQTjiR2l0CdAfOjv75B2A4MFfhK4miK2CRWkS7xF3k7mujj1sR6az1put/GkkQBbCISEj0IZyISEgUwCI1mFkzM7szeuPQw2b2rpk9Em1v+HczC6y1oWQ2fQgnUkW0ZePzRFoorgWmA+2A7xK5aWsOsDi0AiWjKIBFqnuESPj+0N3vrzhoZk8CL0UfvhZCXZKB9CGcSJSZXQasBp5z9xtqef4tIlfB/d39T0HXJ5lHc8AiJ4yP/vmzOp5/P/rn6wBm9jUzW2hme8zMzWxksguUzKIAFjlhEPC+u6+p4/lcYKe7H4g+PpPIfd2KgcMB1CcZRgEsAphZNpH7ttW6y83M8olsxKic/3X3Je4+yd2fA44HUqhkFAWwSMTn0a+z63j+7uif+gBOEkYBLAK4+1FgO9DZzP6x4rhF3A1UfCj3ehj1SWbSMjSRE34GPA4sNrO5wAfAAKA1sJlIgx5dAUvC6ApYJMrdnwD+N7Af+A5wPbAMKCQyP/yuu78bXoWSaXQFLFKFuz8APFD1mJmdT2RueEkoRUnGUgCL1K9X9M9q0w/RlpXdog9PIzJ/fCnwgXoGSyw0BSFSv4oArvkBXGH02OtAS+An0e9/Glxpks50BSxSv1qvgN39JcACr0YyhnpBiIiERFMQIiIhUQCLiIREASwiEhIFsIhISBTAIiIhUQCLiIREASwiEhIFsIhISP4/sUWnwAW3XycAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 360x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(5, 5))\n",
    "\n",
    "ax.hist(\n",
    "    signal_qtilde,\n",
    "    bins=np.linspace(0, 10, 26),\n",
    "    histtype=\"step\",\n",
    "    density=True,\n",
    "    linestyle=\"dashed\",\n",
    "    label=r\"$f\\,(\\tilde{q}_1|1)$\",\n",
    ")\n",
    "ax.hist(\n",
    "    bkg_qtilde,\n",
    "    bins=np.linspace(0, 10, 26),\n",
    "    histtype=\"step\",\n",
    "    density=True,\n",
    "    label=r\"$f\\,(\\tilde{q}_1|0)$\",\n",
    ")\n",
    "\n",
    "ax.set_xlim(0, 9)\n",
    "ax.set_ylim(1e-3, 2)\n",
    "ax.semilogy()\n",
    "ax.set_xlabel(r\"$\\tilde{q}_1$\", fontsize=20)\n",
    "ax.set_ylabel(r\"$f\\,(\\tilde{q}_1|\\theta)$\", fontsize=20)\n",
    "ax.legend(loc=\"best\", fontsize=14);\n",
    "\n",
    "fig.savefig(\"alternative_statistic.pdf\")\n",
    "fig.savefig(\"alternative_statistic.png\")"
   ]
  }
 ],
 "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
}