{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simple example\n", "\n", "## Setting up" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import jax.numpy as jnp \n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from smoofit.model import Variable,Process,Model,ChannelContrib\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Building the model\n", "\n", "Let's define dummy signal and background templates:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAXgklEQVR4nO3df3DV9Z3v8eeLkBKEaoSkjiXSpCs1RfkhBlblWhnoeumCQGewihbBdcQ79WqLTrvRTkd21Bltme3u1tsfUVSca/1Rf1QHuXYdQV38uQnSqhBX6kYJRQmoCFutoO/7R75kIybk5Pzg5Hx9PWYy53x/fb7vnMm8zjef8/l+jiICMzNLl0HFLsDMzPLP4W5mlkIOdzOzFHK4m5mlkMPdzCyFBhe7AICqqqqora0tdhlmZiWlpaVlR0RU97RtQIR7bW0tzc3NxS7DzKykSHq9t23uljEzSyGHu5lZCjnczcxSaED0uZvZwLd3717a29v54IMPil3KZ05FRQU1NTWUl5dnfIzD3cwy0t7ezuc//3lqa2uRVOxyPjMigp07d9Le3k5dXV3Gx7lbxswy8sEHHzBy5EgH+yEmiZEjR/b7PyaHu5llzMFeHNm87g53M7MUcp+7mWVl6vVr2Pru+3lrb1TlUJ5qnN7nftdddx2//vWvKSsrY9CgQfzqV7/ipptu4vLLL2fs2LF5qwdg+PDh7NmzJ69tHiolH+75/gMrtEz/gM0Guq3vvk/b9bPy1l5t48N97vPMM8+watUq1q9fz5AhQ9ixYwcffvghN998c97qSIuSD/d8/4EVWiZ/wGbWs23btlFVVcWQIUMAqKqqAmDatGksX76choYGVqxYwQ033EBlZSUTJkxgyJAh3HjjjSxevJjDDz+c5uZm3nzzTX784x8zf/589uzZw9y5c3nnnXfYu3cv1157LXPnzi3mr5kX7nM3s5JxxhlnsGXLFr7yla/wne98hyeeeOIT2//0pz9xzTXX8Oyzz/LUU0/R2tr6ie3btm1j3bp1rFq1isbGRqBzDPkDDzzA+vXrWbt2LVdccQVp+PpRh7uZlYzhw4fT0tJCU1MT1dXVnH322dx2221d259//nlOP/10RowYQXl5OWedddYnjp83bx6DBg1i7NixvPXWW0DnOPKrrrqK8ePH8/Wvf52tW7d2bStlJd8tY2afLWVlZUybNo1p06Yxbtw4Vq5cmfGx+7tzgK6r8zvuuIOOjg5aWlooLy+ntrY2FXfh+srdzErGK6+8wquvvtq1vGHDBr70pS91LU+ePJknnniCd955h3379nHffff12eauXbv4whe+QHl5OWvXruX113udRbek+MrdzLIyqnJoXgcIjKoc2uc+e/bs4dJLL+Xdd99l8ODBHHvssTQ1NTF//vzONkaN4qqrrmLKlCmMGDGC+vp6jjjiiIO2ed5553HmmWcybtw4GhoaqK+vz8vvU2wOdzPLSjGG9J500kk8/fTTn1r/+OOPdz0/99xzWbJkCfv27eOb3/wm8+bNA/hE3zzQNX69qqqKZ555psfzleoYd3C3jJmlzLJly5g4cSInnHACdXV1XeH+WeMrdzNLleXLlxe7hAGhzyt3SbdI2i7ppR62XSEpJFUly5L0L5I2S/qDpEmFKNrMzA4uk26Z24CZB66UdAxwBvBGt9XfAMYkP0uAX+ReopmZ9Vef4R4RTwJv97Dpp8APgO63cs0Fbo9OzwKVko7OS6VmZpaxrD5QlTQX2BoRvz9g0yhgS7fl9mSdmZkdQv3+QFXSYcBVdHbJZE3SEjq7bhg9enQuTZlZMfx0HOx6o+/9MnXEaFj64kF3KSsrY9y4cUQEZWVl3HjjjZx66qn9PtXixYuZPXt21/j4geLxxx9n+fLlrFq1Kue2shkt81dAHfD75NtBaoD1kqYAW4Fjuu1bk6z7lIhoApoAGhoaSn+WHrPPml1vwLJd+Wtv2cFvNgIYOnQoGzZsAOB3v/sdV1555acmDyu0ffv2MXjwwB9o2O9umYh4MSK+EBG1EVFLZ9fLpIh4E3gIOD8ZNXMysCsituW3ZDMzeO+99zjyyCOBzpuNZsyYwaRJkxg3bhwPPvhg1363334748ePZ8KECSxcuPBT7fzoRz9i8eLFfPTRR6xevZr6+npOOukkLrvsMmbPng10jp1fuHAhU6dOZeHChbS1tTF9+nTGjx/PjBkzeOONzv9gFi9ezL333tvV9vDhw4HOK/Jp06Yxf/586uvrOe+887rmtnnkkUeor69n0qRJ3H///Xl7ffp8+5F0JzANqJLUDlwdESt62X018LfAZuDPwAV5qtPMjPfff5+JEyfywQcfsG3bNtasWQP897S9hx9+ODt27ODkk09mzpw5bNy4kWuvvZann36aqqoq3n77k2NDvv/977N7925uvfVW/vKXv3DxxRfz5JNPUldXx4IFCz6x78aNG1m3bh1Dhw7lzDPPZNGiRSxatIhbbrmFyy67jN/+9rcHrf2FF17g5Zdf5otf/CJTp07lqaeeoqGhgYsuuog1a9Zw7LHHcvbZZ+fttcpktMyCiDg6IsojoubAYE+u4HckzyMiLomIv4qIcRHRnLdKzewzb3+3TGtrK4888gjnn38+EdHrtL1r1qzhrLPO6vpSjxEjRnS1dc0117Br1y5++ctfIonW1la+/OUvU1dXB/CpcJ8zZw5Dh3bOf/PMM89w7rnnArBw4ULWrVvXZ+1TpkyhpqaGQYMGMXHiRNra2mhtbaWuro4xY8YgiW9/+9t5eZ3Ad6iaWYk65ZRT2LFjBx0dHaxevbrf0/ZOnjyZlpYW3n777U+Efm+GDRvW5z6DBw/m448/BuDjjz/mww8/7NrWfbrhsrIy9u3b12d7ufDcMmZWklpbW/noo48YOXJkr9P2Tp8+nd/85jfs3LkT4BPdMjNnzqSxsZFZs2axe/dujjvuOF577TXa2toAuPvuu3s996mnnspdd90FdM4Hf9pppwFQW1tLS0sLAA899BB79+496O9QX19PW1sbf/zjHwG48847s3gleuYrdzPLzhGjMxrh0q/2+rC/zx06v2xj5cqVlJWV9Tpt7/HHH88Pf/hDTj/9dMrKyjjxxBM/MTvkWWedxe7du5kzZw6rV6/m5z//OTNnzmTYsGFMnjy51zp+9rOfccEFF/CTn/yE6upqbr31VgAuuugi5s6dy4QJE7raOZiKigqampqYNWsWhx12GKeddhq7d+/u83XIhAbCdwU2NDREc3N23fO1jQ+X3Bdkl1K9Zvtt2rSJr371q8Uuo6D27NnD8OHDiQguueQSxowZw9KlS4tdFtDz6y+pJSIaetrf3TJmZombbrqJiRMncvzxx7Nr1y4uvvjiYpeUNXfLmJklli5dOmCu1HPlK3czy9hA6Mb9LMrmdXe4m1lGKioq2LlzpwP+EIsIdu7cSUVFRb+Oc7eMmWWkpqaG9vZ2Ojo6il3KZ05FRQU1NTX9OsbhbmYZKS8v77p70wY+d8uYmaWQw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlUJ/hLukWSdslvdRt3U8ktUr6g6QHJFV223alpM2SXpH0PwtUt5mZHUQmV+63ATMPWPcocEJEjAf+A7gSQNJY4Bzg+OSYn0sqy1u1ZmaWkT7DPSKeBN4+YN2/RsT+b3d9Ftg/o81c4K6I+EtE/CewGZiSx3rNzCwD+ehz/zvg/yXPRwFbum1rT9Z9iqQlkpolNXuWOTOz/Mop3CX9ENgH3NHfYyOiKSIaIqKhuro6lzLMzOwAWU/5K2kxMBuYEf89e/9W4Jhuu9Uk68zM7BDK6spd0kzgB8CciPhzt00PAedIGiKpDhgDPJ97mWZm1h99XrlLuhOYBlRJageupnN0zBDgUUkAz0bE/4qIlyXdA2yks7vmkoj4qFDFm5lZz/oM94hY0MPqFQfZ/zrgulyKMjOz3PgOVTOzFHK4m5mlkMPdzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyOFuZpZCDnczsxRyuJuZpZDD3cwshRzuZmYp5HA3M0shh7uZWQo53M3MUsjhbmaWQg53M7MUcribmaWQw93MLIX6DHdJt0jaLumlbutGSHpU0qvJ45HJekn6F0mbJf1B0qRCFm9mZj3L5Mr9NmDmAesagcciYgzwWLIM8A1gTPKzBPhFfso0M7P+6DPcI+JJ4O0DVs8FVibPVwLzuq2/PTo9C1RKOjpPtZqZWYYGZ3ncURGxLXn+JnBU8nwUsKXbfu3Jum0cQNISOq/uGT16dJZlwLohl8Gyc7M+/lBrqwCWFbsKMxswjhgNS1/Me7PZhnuXiAhJkcVxTUATQENDQ7+P369GO2DZrmwPP+RqGx+m7fpZxS7DzAaKZUcUpNlsR8u8tb+7JXncnqzfChzTbb+aZJ2ZmR1C2Yb7Q8Ci5Pki4MFu689PRs2cDOzq1n1jZmaHSJ/dMpLuBKYBVZLagauB64F7JF0IvA58K9l9NfC3wGbgz8AFBajZzMz60Ge4R8SCXjbN6GHfAC7JtSgzM8uN71A1M0shh7uZWQo53M3MUsjhbmaWQg53M7MUcribmaWQw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOzFHK4m5mlkMPdzCyFHO5mZinkcDczSyGHu5lZCuUU7pKWSnpZ0kuS7pRUIalO0nOSNku6W9Ln8lWsmZllJutwlzQKuAxoiIgTgDLgHOAG4KcRcSzwDnBhPgo1M7PM5dotMxgYKmkwcBiwDZgO3JtsXwnMy/EcZmbWT1mHe0RsBZYDb9AZ6ruAFuDdiNiX7NYOjOrpeElLJDVLau7o6Mi2DDMz60Eu3TJHAnOBOuCLwDBgZqbHR0RTRDREREN1dXW2ZZiZWQ9y6Zb5OvCfEdEREXuB+4GpQGXSTQNQA2zNsUYzM+unXML9DeBkSYdJEjAD2AisBeYn+ywCHsytRDMz669c+tyfo/OD0/XAi0lbTcDfA5dL2gyMBFbkoU4zM+uHwX3v0ruIuBq4+oDVrwFTcmnXzMxy4ztUzcxSyOFuZpZCDnczsxRyuJuZpZDD3cwshRzuZmYp5HA3M0shh7uZWQo53M3MUsjhbmaWQg53M7MUcribmaWQw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkIOdzOzFMop3CVVSrpXUqukTZJOkTRC0qOSXk0ej8xXsWZmlplcr9z/GXgkIuqBCcAmoBF4LCLGAI8ly2ZmdghlHe6SjgC+BqwAiIgPI+JdYC6wMtltJTAvtxLNzKy/crlyrwM6gFslvSDpZknDgKMiYluyz5vAUT0dLGmJpGZJzR0dHTmUYWZmB8ol3AcDk4BfRMSJwH9xQBdMRAQQPR0cEU0R0RARDdXV1TmUYWZmB8ol3NuB9oh4Llm+l86wf0vS0QDJ4/bcSjQzs/7KOtwj4k1gi6TjklUzgI3AQ8CiZN0i4MGcKjQzs34bnOPxlwJ3SPoc8BpwAZ1vGPdIuhB4HfhWjucwM7N+yincI2ID0NDDphm5tGtmZrnxHapmZinkcDczSyGHu5lZCjnczcxSyOFuZpZCDnczsxRyuJuZpZDD3cwshXK9Q9X6aVTlUGobHy52GRkZVTmUpxqnF7sMM8uCw/0QK6WwLJU3ITP7NHfLmJmlkMPdzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyOFuZpZCDnczsxTKOdwllUl6QdKqZLlO0nOSNku6O/l+VTMzO4TyceX+XWBTt+UbgJ9GxLHAO8CFeTiHmZn1Q07hLqkGmAXcnCwLmA7cm+yyEpiXyznMzKz/cr1y/yfgB8DHyfJI4N2I2JcstwOjcjyHmZn1U9bhLmk2sD0iWrI8fomkZknNHR0d2ZZhZmY9yOXKfSowR1IbcBed3TH/DFRK2j/bZA2wtaeDI6IpIhoioqG6ujqHMszM7EBZh3tEXBkRNRFRC5wDrImI84C1wPxkt0XAgzlXaWZm/VKIce5/D1wuaTOdffArCnAOMzM7iLx8WUdEPA48njx/DZiSj3bNzCw7vkPVzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyOFuZpZCDnczsxRyuJuZpVBebmKydBpVOZTaxoeLXUbGRlUO5anG6cUuw2xAcLhbr0otKEvpjcis0NwtY2aWQg53M7MUcribmaWQw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFLI4W5mlkJZh7ukYyStlbRR0suSvpusHyHpUUmvJo9H5q9cMzPLRC5X7vuAKyJiLHAycImksUAj8FhEjAEeS5bNzOwQyjrcI2JbRKxPnu8GNgGjgLnAymS3lcC8HGs0M7N+ysvEYZJqgROB54CjImJbsulN4KhejlkCLAEYPXp0Psqwz7hSmsXSM1haoeUc7pKGA/cB34uI9yR1bYuIkBQ9HRcRTUATQENDQ4/7mPVHKYVlqbwJWenKabSMpHI6g/2OiLg/Wf2WpKOT7UcD23Mr0czM+iuX0TICVgCbIuIfu216CFiUPF8EPJh9eWZmlo1cumWmAguBFyVtSNZdBVwP3CPpQuB14Fs5VWhmZv2WdbhHxDpAvWyekW27ZmaWO9+hamaWQg53M7MUcribmaWQw93MLIUc7mZmKeRwNzNLIYe7mVkKOdzNzFIoL7NCmln/lNIMluBZLEuRw92sCEotKEvpjcg6uVvGzCyFHO5mZinkcDczSyGHu5lZCjnczcxSyKNlzKxPpTZ0s5S0VRSmXYe7mfWp1IZulpRlhWnW3TJmZinkcDczS6GChbukmZJekbRZUmOhzmNmZp9WkHCXVAb8H+AbwFhggaSxhTiXmZl9WqGu3KcAmyPitYj4ELgLmFugc5mZ2QEKNVpmFLCl23I78Nfdd5C0BFiSLO6R9EqW56riH7Qjy2OLoQoolXpLqVYorXpLqVYorXpLqVbILcO+1NuGog2FjIgmoCnXdiQ1R0RDHko6JEqp3lKqFUqr3lKqFUqr3lKqFQpXb6G6ZbYCx3RbrknWmZnZIVCocP93YIykOkmfA84BHirQuczM7AAF6ZaJiH2S/jfwO6AMuCUiXi7EuchD184hVkr1llKtUFr1llKtUFr1llKtUKB6FRGFaNfMzIrId6iamaWQw93MLIVKOtxLaYoDSbdI2i7ppWLX0hdJx0haK2mjpJclfbfYNfVGUoWk5yX9Pqn1H4pdUyYklUl6QdKqYtdyMJLaJL0oaYOk5mLX0xdJlZLuldQqaZOkU4pdU08kHZe8pvt/3pP0vbyeo1T73JMpDv4D+Bs6b5L6d2BBRGwsamG9kPQ1YA9we0ScUOx6DkbS0cDREbFe0ueBFmDeQHxtJQkYFhF7JJUD64DvRsSzRS7toCRdDjQAh0fE7GLX0xtJbUBDRJTETUGSVgL/FhE3JyP1DouId4tc1kElWbYV+OuIeD1f7ZbylXtJTXEQEU8Cbxe7jkxExLaIWJ883w1sovOu4wEnOu1JFsuTnwF9xSKpBpgF3FzsWtJE0hHA14AVABHx4UAP9sQM4I/5DHYo7XDvaYqDARlApUxSLXAi8FyRS+lV0sWxAdgOPBoRA7bWxD8BPwA+LnIdmQjgXyW1JFOGDGR1QAdwa9LldbOkYcUuKgPnAHfmu9FSDncrMEnDgfuA70XEe8WupzcR8VFETKTzTugpkgZst5ek2cD2iGgpdi0Z+h8RMYnOGV4vSboXB6rBwCTgFxFxIvBfwED/LO5zwBzgN/luu5TD3VMcFFDSf30fcEdE3F/sejKR/Au+FphZ5FIOZiowJ+nLvguYLun/Frek3kXE1uRxO/AAnd2hA1U70N7tP7d76Qz7gewbwPqIeCvfDZdyuHuKgwJJPqRcAWyKiH8sdj0HI6laUmXyfCidH7C3FrWog4iIKyOiJiJq6fybXRMR3y5yWT2SNCz5QJ2ke+MMYMCO9oqIN4Etko5LVs0ABtwggAMsoABdMlDCX5B9iKc4yJmkO4FpQJWkduDqiFhR3Kp6NRVYCLyY9GUDXBURq4tXUq+OBlYmIw4GAfdExIAeXlhCjgIe6HyvZzDw64h4pLgl9elS4I7kgu814IIi19Or5A3zb4CLC9J+qQ6FNDOz3pVyt4yZmfXC4W5mlkIOdzOzFHK4m5mlkMPdzCyFHO5mZinkcDczS6H/D5pqszsrVMnKAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "sig_shape = np.array([150., 50., 30., 20., 10., 5., 5.])\n", "bkg_shape = np.array([100.]*7)\n", "x = np.arange(0, 8)\n", "_,bins,_ = plt.hist(x[:-1], weights=sig_shape, bins=x, label=\"Signal\", histtype=\"step\")\n", "plt.hist(x[:-1], weights=bkg_shape, bins=bins, label=\"Background\", histtype=\"step\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using those shapes, we create the signal and background processes:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": true }, "outputs": [], "source": [ "sig = Process(\"sig\")\n", "bkg = Process(\"bkg\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now create the signal strength variable $\\mu$, which has a default value of 1, and indicate that the signal template should be scaled linearly by the signal strength. We also specify that the signal strength should be positive:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "WARNING:absl:No GPU/TPU found, falling back to CPU. (Set TF_CPP_MIN_LOG_LEVEL=0 and rerun for more info.)\n" ] } ], "source": [ "mu = Variable(\"mu\", 1., lower_bound=0.)\n", "sig.scale_by(mu)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a next step, we'll register the templates. We do that by defining the \"contribution\" of the signal and background processes to a \"channel\". In this example there is a single channel, so we'll just give it a dummy name:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "sig_contrib = ChannelContrib(\"channel\", sig_shape)\n", "sig.add_contrib(sig_contrib)\n", "bkg_contrib = ChannelContrib(\"channel\", bkg_shape)\n", "bkg.add_contrib(bkg_contrib)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we add some systematic uncertainties. Let's start with the luminosity, as a simple log-normal uncertainty of 5%. We register those uncertainties with the `ChannelContrib` objects we created above:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true }, "outputs": [], "source": [ "lumi_nuis = Variable(\"lumi\", 0., nuisance=True)\n", "lumi_unc = 1.05\n", "sig_contrib.add_lnN(lumi_nuis, lumi_unc)\n", "bkg_contrib.add_lnN(lumi_nuis, lumi_unc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We could go on like this, but it is also possible to add several sources in one go, as an \"array\". Let's define some shape uncertainties. The shapes for different sources should be stacked vertically:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "scrolled": true }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD4CAYAAAAXUaZHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjIUlEQVR4nO3de3yU5Z338c+Pg4CAooCUgzJYqWwRiCSytmpLtbW4BJDXahUtQt0SbG2tdbdrtNsmaPusad2uu91WCeupz1LtrmhLgAd1q1bRqk2QekBs1YZKQIhQEFpRDr/njxnCBHIYZubOPfc93/frNa9MrrkP3xnbH3euue7rMndHRETipVvYAUREJP9U3EVEYkjFXUQkhlTcRURiSMVdRCSGeoQdAGDQoEGeSCTCjiEiEikNDQ3vuPvgtl4riOKeSCSor68PO4aISKSY2fr2XlO3jIhIDKm4i4jEkIq7iEgMFUSfu4hEw549e9iwYQO7d+8OO0pR6d27NyNGjKBnz54Z76PiLiIZ27BhA/379yeRSGBmYccpCu7O1q1b2bBhA6NGjcp4P3XLiEjGdu/ezcCBA1XYu5CZMXDgwCP+a0nFXUSOiAp718vmM1e3jIhkJVG5PJDjNt4yNZDjFhtduYtIpHz3u99l7NixjB8/npKSEp577jm++MUvsnbt2ryfq1+/fnk/Zpdx99AfpaWlnpOJE90h+Rg6NNlWVXWwDdzr65OP9LaqquS2Q4cebJs4Mdk2b17rbZua3Jcubd22cGFy2/S28vJkW3l563b35PYjR+b2XkVCtHbt2pbnI69f5iOvX5a3Y2dyvGeeecbPPPNM3717t7u7Nzc3e1NTU94yHKpv376BHftIpX/2BwD13k5djUe3TEPD4W3V1cnHodpaeWrjxsPbamuTj3TDhrW9f1ttdXWHt1VUwPz5h7eLSEY2bdrEoEGD6NWrFwCDBg0CYPLkydx6662UlZVx5513UlNTw4ABA5gwYQK9evXiP/7jP5g7dy7HHHMM9fX1vP3223zve9/joosuYteuXcyYMYM//elP7Nmzh+985zvMmDEjzLeZF512y5jZXWa2xcxeTmv7mZmtST0azWxNqj1hZu+lvXZHgNkPqqjoktOISLjOP/983nrrLT7ykY/w5S9/mV/96letXt+4cSM333wzzz77LE8//TTr1q1r9fqmTZtYtWoVy5Yto7KyEkiOIX/ooYdYvXo1jz/+OH//93+Px2D50Uz63O8BpqQ3uPsl7l7i7iXAEuDBtJffOPCau1+Vt6QdWbSoS06TF0uXhp1AJLL69etHQ0MDtbW1DB48mEsuuYR77rmn5fXnn3+eT37ykxx//PH07NmTiy++uNX+F154Id26deOjH/0omzdvBpJd0zfeeCPjx4/n05/+NE1NTS2vRVmn3TLu/qSZJdp6zZLjcz4HnJvnXPFVWhp2ApFI6969O5MnT2by5MmMGzeOe++9N+N9D3TnAC1X54sXL6a5uZmGhgZ69uxJIpGIxR24ufa5nwNsdvffp7WNMrMXgHeBf3L3p9ra0cwqgAqAk046KccYETJ8eNt99CIRFdSQyLa89tprdOvWjdGjRwOwZs0aRo4cycsvJ3uNzzjjDK699lr+9Kc/0b9/f5YsWcK4ceM6POaOHTs44YQT6NmzJ48//jjr17c7i26k5FrcZwH3pf2+CTjJ3beaWSnwczMb6+7vHrqju9cCtQBlZWW5Vbumppx2F5Fo2LVrF1/96lfZvn07PXr04JRTTqG2tpaLLroIgOHDh3PjjTcyadIkjj/+eMaMGcOxxx7b4TEvv/xypk2bxrhx4ygrK2PMmDFd8VYCZ5l8cZDqllnm7qeltfUAmoBSd9/Qzn5PAP/g7h2uxFFWVuY5LdZRVwfTpmW/f1cy05W7RNarr77KX/3VX4Udo0O7du2iX79+7N27l5kzZ3LllVcyc+bMsGPlrK3P3swa3L2sre1zuYnp08C69MJuZoPNrHvq+cnAaODNHM6RmenTAz9F3sybF3YCkVirrq6mpKSE0047jVGjRnHhhReGHSkUnXbLmNl9wGRgkJltAKrc/U7gUlp3yQB8ArjJzPYA+4Gr3H1bfiNH3KFj50Ukr2699dawIxSETEbLzGqnfW4bbUtIDo2U9pSWtn3TlYhIHsVibpnay3+FGS2PurrkTafpbQfucyotPdg2bFiyrbq69bYNDclHetuBm12HDTvYdmBUY0VF6203bkxmSG87cMFuq1XYRSR4GX2hGrScv1CNEH2fKlEWhS9U4+pIv1CNxdwyKpgiIajueIhh9sfdEcxxi0wsumWipLw87AQi0dW9e3dKSkqYMGECEydO5JlnnsnqOHPnzuWBBx7Ic7rcPfHEE5TnqUjE4so9SupKq4HqkFOI5FG+rrQz+EugT58+rFmzBoCHH36YG2644bDJw4K2d+9eevQo/NIZiyv38hOeDztCxqYt0NwyIvnw7rvvctxxxwHJG5fOO+88Jk6cyLhx4/jFL37Rst1PfvITxo8fz4QJE5g9e/Zhx/nWt77F3Llz2bdvHytWrGDMmDGUlpZyzTXXtFxFV1dXM3v2bM466yxmz55NY2Mj5557LuPHj+e8887jj3/8I3D4XwQHFvt44oknmDx5MhdddBFjxozh8ssvb5nbZuXKlYwZM4aJEyfy4IMPki+F/89PBuo2Two7QsaWEZE7aUUK0HvvvUdJSQm7d+9m06ZNPPbYY8DBaXuPOeYY3nnnHc4880ymT5/O2rVr+c53vsMzzzzDoEGD2Lat9W033/jGN9i5cyd3330377//PvPnz+fJJ59k1KhRzJrVehT42rVrWbVqFX369GHatGnMmTOHOXPmcNddd3HNNdfw85//vMPsL7zwAq+88grDhg3jrLPO4umnn6asrIx58+bx2GOPccopp3DJJZfk7bOKxZX7tCHRuXIXkewd6JZZt24dK1eu5IorrmhZeaitaXsfe+wxLr744pZFPY4//viWY918883s2LGDO+64AzNj3bp1nHzyyYwaNQrgsOI+ffp0+vTpA8Cvf/1rLrvsMgBmz57NqlWrOs0+adIkRowYQbdu3SgpKaGxsZF169YxatQoRo8ejZnx+c9/Pi+fE8Tkyn3ZluhcuYtIfnzsYx/jnXfeobm5mRUrVhzxtL1nnHEGDQ0NbNu2rVXRb0/fvn073aZHjx7s378fgP379/PBBx+0vJY+3XD37t3Zu3dvp8fLRSyKe5R4fQOgfneJkaCGRHZi3bp17Nu3j4EDB7Y7be+5557LzJkzue666xg4cGCrQj5lyhQ++9nPMnXqVB555BFOPfVU3nzzTRobG0kkEvzsZz9r99wf//jHuf/++5k9ezaLFy/mnHPOASCRSNDQ0MDnPvc5li5dyp49ezp8D2PGjKGxsZE33niDD3/4w9x336EzumRPxb2L1ZYtpMI1v4xINg70uUNysY17772X7t27tztt79ixY/nmN7/JJz/5Sbp3787pp5/eauWmiy++mJ07dzJ9+nRWrFjBj3/8Y6ZMmULfvn0544wz2s3xwx/+kC984Qt8//vfZ/Dgwdx9990AzJs3jxkzZjBhwoSW43Skd+/e1NbWMnXqVI4++mjOOeccdu7cmduHlKI7VLuYbriSKIv7HaoHpgt2d66++mpGjx7N17/+9bBjAV075W/BqP38k2FHEJEYWLRoESUlJYwdO5YdO3Ywf/78sCNlLRZX7lG6Go5SVpFDxf3KvZAV5ZV7lCxdGnYCkdwUwgVhscnmM1dx72KlXzw97AgiWevduzdbt25Vge9C7s7WrVvp3bv3Ee0Xi9EyS//peSAaY92Hb3kB/d9ComrEiBFs2LCB5ubmsKMUld69ezNixIgj2icWxb105klhRxApCj179my5g1MKWyy6ZYaXfijsCCIiBSUWxT1K5s0LO4GIFINOi7uZ3WVmW8zs5bS2ajNrMrM1qcffpL12g5m9bmavmdlngwoeVbVUhB1BRIpAJlfu9wBT2mj/V3cvST1WAJjZR4FLgbGpfX5sZt3zFbY988ZE5yam0kXRvSlCRKKj0y9U3f1JM0tkeLwZwP3u/j7wBzN7neQwll9nH7FjicrlMCP1MwLWMzXsCCJSBHLpc/+Kmb2Y6rY5LtU2HHgrbZsNqbbDmFmFmdWbWX2uw6qO/peROe0vIhI32Q6FvB24GfDUz38BrjySA7h7LVALyekHsswBwKt7T8NvOS2XQ3SJROVyeh79F+DosKOISMxldeXu7pvdfZ+77wcWcfAOoibgxLRNR6TaJGX2BbeEHUFEikBWxd3Mhqb9OhM4MJJmKXCpmfUys1HAaCDwNfCG2NtBnyJvTlyi0aciErxMhkLeR/IL0VPNbIOZ/R3wPTN7ycxeBD4FfB3A3V8B/htYC6wErnb3fYGlT+n9jw1BnyJvFlAddgQRKQKZjJaZ1UbznR1s/13gu7mEOlJjfrob1NshItIiFn0ED7/1t2FHEBEpKLEo7lHyoTlPhR1BRIqAinsXW3bvtWFHEJEiEIvifnb5fWFHyFgZ0fnyV0SiKxbFXUREWotFcV+1rK0BPSIixSsWxR3ADBoakg+zg4/q6uTrw4YdbCstTbZVVLTeduNGqKtr3VZbe/D4Bx7TpiXbpk1r3Q7J7dPb6uqSxzWD9TVTsV57uvRzEZHiZIWw0G1ZWZnX19dntW+icjnra6ZSAG+jU4nK5Yz98Qcsf3dm2FFEJAbMrMHdy9p6LRZX7h8+Zm3YETL2ws6/DjuCiBSBWBT3vV/6Q9gRMraJYWFHEJEiEIvi/lZNWwtFiYgUr1gU9/0EvpKfiEikxKK4R0m/XjvCjiAiRSAWxb0374UdIWOz3v9Z2BFEpAjEorgPuf6xsCNkbBEVYUcQkSIQi+L+wW3jwo4gIlJQYlHcN71/UtgRREQKSiyKO5C8v78i1eVRWnrw/v9hqXHl1dWt5wUIYa6Cxpryrvo0RKTIafqBLpSoXM7Xah7l635b2FFEJAZiP/3Ax8t/GnaEjF3HbWFHEJEi0GlxN7O7zGyLmb2c1vZ9M1tnZi+a2UNmNiDVnjCz98xsTepxR4DZW/R7oVdXnEZEJDIyuXK/Bzj0/v5HgdPcfTzwO+CGtNfecPeS1OOq/MTs2CNNWiBbRCRdp8Xd3Z8Eth3S9oi77039+iwwIoBsIiKSpXz0uV8J/L+030eZ2Qtm9iszO6e9ncyswszqzay+ubk5DzGioYrqsCOISBHIqbib2TeBvcDiVNMm4CR3Px24DvipmR3T1r7uXuvuZe5eNnjw4FxicM6xj+S0f1daoOIuIl0g6+JuZnOBcuByT42ndPf33X1r6nkD8AbwkTzk7NAfr9LSdSIi6bIq7mY2BfhHYLq7/yWtfbCZdU89PxkYDbyZj6AdWV8zNehTiIhESo/ONjCz+4DJwCAz2wBUkRwd0wt41JIrQz+bGhnzCeAmM9sD7AeucvdtbR5YREQC02lxd/dZbTTf2c62S4AluYaKs4k0AKVhxxCRmIvFHaqD2BJ2hIwNY2PYEUSkCMSiuPe9/jdhR8jYMqaFHUFEikAsivufa84IO4KISEGJRXF/hxPCjiAiUlBiUdxFRKS12BT3NtbGaLVmBiTX0khvq6tLrrmR3hbkeh9bV45rlTWL9T4K7j3lYQ0TvSe9p6J+T4kEgYjFYh0AjbcU/o1MUVpYRES6Rmlp8h+AbMR+sY6e/54IO4KISFayLeydiUVxf/29sWFHEBHJyoFunnyLRXEXEYmqRYuCOa6Kexer4RthRxCRIhCL4n52v/+NzFfk1/P9+H7tr/ek96T3dOTvKSCxGC0z6pltPP7k7Dynyj+NlhGRQ220YQzz7Oaciv1omSeeKvzCLiLSloZLbw3kuLEo7iIiUTX9/ssCOa6Kexcrpy7sCCJSBGJR3I9he9gRMraRYWFHEJEiEIviftz1T4cdIWOrtQqTiKRZSDB3McWiuP+xZkrYEUREslIxMZj5B2JR3J3uYUcQEcmKrQ6xuJvZXWa2xcxeTms73sweNbPfp34el2o3M/t3M3vdzF40s4mBJBcRkXZleuV+D3Bo30cl8Et3Hw38MvU7wAXA6NSjArg995idic5dQVVUhx1BRIpAj0w2cvcnzSxxSPMMYHLq+b3AE8D1qfafePLW12fNbICZDXX3TXlJfIjG3pdBFUShZjb2BmMH1dXHhh1FRApEOf8FTMv7cXPpcx+SVrDfBoakng8H3krbbkOqrRUzqzCzejOrb25uziEGJG76bU77i4iEpY7pgRw3oyv3zri7m9kR9Y24ey1QC8m5ZXI5/3pPQPWOXA7RJQ6sGhWFrCLSBaqPZdrJKwK5tTGXK/fNZjYUIPVzS6q9CTgxbbsRqTYRETnEsjcvCOS4uRT3pcCc1PM5wC/S2q9IjZo5E9gRVH97FP2Aa8OOICJFINOhkPcBvwZONbMNZvZ3wC3AZ8zs98CnU78DrADeBF4HFgFfznvqQ9yUqAr6FHlzHbeFHUFEikCmo2VmtfPSeW1s68DVuYQSESkWjhHEcO5Y3KH67cYFYUcQEclK7bhgOjdiUdxFRKJq/ks/CuS4Ku5dbB61YUcQkSIQi+I+0hrDjpCxZZSHHUFEikAsinvjtyeEHSFjm7RYh4ikWRrA1AMQk+JuC7Ynn1RXg9nBR0ND8pHeVl2d3HbYsINtpakFNCoqWm+7cSPU1bVuq011q6S3TUv9x5k2rXU7JLdP/d5Yk7pq37ix9XYVqcn6S0sPtg0bFon3hFnyeHpPek96T1m9p9KBLxAES45cDFdZWZnX19dnt3P1sdiCHRTA2+hUonI562umRiKriHSBHOuXmTW4e1lbr8Xiyj1KhrIx7AgiUgRiUdyNfWFHyFg5y8KOICJFIBbFfX/V8WFHyNiigBbDFZFoCmp4dCyK+4AFjWFHEBHJSi3zAzluLIr7Do4LO4KISFZKyXIwSSdiUdxFRKJqNaWBHFfFvYtpPncR6QqxKO4rPzY17AgZ03zuIpIuqOHRsSjud625IuwIIiJZ2cjwQI4bi+L+3+9dEnYEEZGsVPeuDuS4sSjuEI1pJtbXJLuPimzqDL0nvSe9p/be04IdLNgdzDKhmlumCyUqlzO35jdUe3XYUUSkEBTi3DJmdqqZrUl7vGtm15pZtZk1pbX/TbbnyFTC3gz6FHmzgOqwI4hIEchogey2uPtrQAmAmXUHmoCHgC8A/+rut+YjYCb+8O3TgR1ddToRkYKXrz7384A33H19no53RGyBCruISLp8FfdLgfvSfv+Kmb1oZneZ2XF5OoeIiGQo5+JuZkcB04H/STXdDnyYZJfNJuBf2tmvwszqzay+ubk51xiRMZGGsCOISBHIx5X7BcBqd98M4O6b3X2fu+8HFgGT2trJ3WvdvczdywYPHpxTgCjN5z5Mi3WISBfIR3GfRVqXjJkNTXttJvByHs7RoSjN574soMVwRUTS5VTczawv8BngwbTm75nZS2b2IvAp4Ou5nCMT3RZsC/oUIiKRkvVQSAB3/zMw8JC22TklyiYH3bv6lCIiBS020w+IiMhBKu5drIZvhB1BRIpALIq7Vx0bdoSMXc/3w44gIkUgFsV91E0vhB1BRKSgxKK4N/rJYUcQESkosSjuIiLSmop7F6vhG4WzUEAGix8kEl30wYhIXsVisY5z/nkZT71/Tn5DBSBRuZzGmnIisbJIilmk4opESyEu1lFIvjSuNuwIGTNUKUUkeLEo7pc3/N+wI8SWrtpFoikWxV2CUxudP4pEJI2Kexcrpy7sCEdk/vywE4hINmJR3HvyftgRMrb4qEtzG76S3g7J7dPb6uqSx83XkBwRiaRYjJZJ/iz8dVQTlcvZ8kAZf3l9SNhRMqbRMiIB0miZjtmC7WFHyNh7b0SnsAMsXRp2AhHJRiyKO6j7ICilXzw97AgikoWYFHcJyvAtmpRNJIpiUtyj0yn84pBTwo4gIkUgFsXdqwaEHSFjk0se7HwjEZEcxaK4D1iwPuwIGdv28PiwIxyRefPCTiAi2ci5uJtZo5m9ZGZrzKw+1Xa8mT1qZr9P/Twu96jt28GAIA9f1GqpCDuCiGQhX1fun3L3krTxlpXAL919NPDL1O8SQaWLdIuqSBT1COi4M4DJqef3Ak8A1wd0rkhZyjQSlVeFHSNj65kadgQRyUI+rtwdeMTMGszswN/wQ9x9U+r528Bhd+6YWYWZ1ZtZfXNzc04Brhp0R077d6Vvzbki7AgiUgTyceV+trs3mdkJwKNmti79RXd3MztsrKK71wK1kJx+IA85Cl7jLVMjdTt/onI5PY/+C3B02FFE5AjlfOXu7k2pn1uAh4BJwGYzGwqQ+rkl1/N05I53otPNETWzL7gl7AgikoWciruZ9TWz/geeA+cDLwNLgTmpzeYAv8jlPBKeE5fEYrSsSNHJ9f+5Q4BVZvZb4HlgubuvBG4BPmNmvwc+nfpdgHlEa/WLBVSHHUFEspBTn7u7vwlMaKN9K3BeLsc+EseyHSIy1r12YdgJRKQYxOJv7u1VI8OOkLHShRVHvF6GWXIbSO5zoK20NNlWUdF623ytAbK+RsMgRaIqFot12ILtuEdk2t8IDZdJVC5nfc3UqMQViR4t1tGZiBT2CPoB14YdQUSyEJPiLkG5jtvCjiAiWYhJcY9Qv0F5edgJRKQIxKK4R2k+d+rqwk4gIkUgFsX9qAWB3gCbX+0NTUlvq6tLDnlJb6tITdtTWnqwLeDhNo01+itDJKpiMlom+2+bpX2JyuWcX9NErWtOd5FAaLSMhGUZunoXiSIVd+nQJoaFHUFEshCL4r64dHbYEURECkosivvtL6lPWEQkXSyK+6oPzgk7QmwNZWNBjuRpd+IcEQGCW0NVYuL0/s/Bu218ld/W1/sNDYe3VVcfLN6d7b9x4+FttbUHi/cBw4ZFZn4ekbDE4spdgrNi58ywI4hIFmJR3BP2ZtgRYi2f0wgH2cOT6PN2l34uIoUsFjcxJX/uyF8gAQ5O+Tvy+uVhR8nI9lWj2b7qI2HHEMmcbmLqmC1QYRd44O3Phx1BpGDoC1VpV+MtU2moKaX0lja+KC0wicrlXPHGz2njK1mRohSLK3cJUH3hF/YDdDetyEFZF3czO9HMHjeztWb2ipl9LdVebWZNZrYm9fib/MVtJwv7gj5F0SorC3cd10y/gN10z9ld9pmIREHWX6ia2VBgqLuvNrP+QANwIfA5YJe735rpsfSFagGLyJqvWu9VIqkQv1B1903uvjr1fCfwKjA82+PlotuCbWGcVgpMv176B17kgLz0uZtZAjgdeC7V9BUze9HM7jKz49rZp8LM6s2svrm5OafzO91z2l86MHJk4U85YEZjTTkXf/BAdKZJSCQC+c8lckDO49zNrB/wK+C77v6gmQ0B3iG5sOnNJLturuzoGFqsQ3IVuW6ZiHR3ScAKsVsmdeCewBJgsbs/CODum919n7vvBxYBk3I5h0gsaaF0CVguo2UMuBN41d1/kNY+NG2zmcDL2cfLjFcdG/QpJCKiME2CGVSXaqF0CVYuo2XOBp4CXgL2p5pvBGYBJSS7ZRqB+e6+qaNj5dotM+qmF/jD/pOz219iIXLdMtOmJf/VkOIWYLdM1neouvsqoK0JtFdke8xsNboKu0A9pSQqbwo7RkY+WLZQd9NKoHSHqsRG+Zzbwo6QMd1NK0HT3DISC423TI3MAJREZTRm2ZRoi8WV+9lHPRV2BJEjctQQ3XAlwYpFcf/SuNrON5LYG0ljQc9/c2BkzfqaqRy1vas+FSlWsVisQzcxCRCZG4MSlcs5v6aJWq8IO4qErVBvYhIpKEuXhp0gY4tQYZdgqbhLfKTfUVRoc8kc2scDhT//TVuZJTJi0S1z1IItfOC98htKJCCRu+EKdNNVUNQt07EPqk4IO4JIvC1bFnYCOUKxKO62YHvYEUSOyA+4NuwIEnOxKO5tz4IgUrj+oe8/t3RpD5zyIonK5a26uY8+ZTOJyuUcfcrmVu2JyuUMnPJiq7YTLvoNI67+31Zt/Uv+SKJyOb0+tKOlrUf/3SQqlzPg7N+12nbo3KcYOvepVm0Dzv4dicrl9Oi/O3le/hD2RyZHKCbFXSRa9v+5D47hGFsfnsB5rz/HGV/+SUvbX974EP9n5Q854aJ6XhxySkv7cz+6gv4lb/GvZ13W0rZlySRK/vwKU+f8W0vbzt+O5NpVixk6dxVv9xuIY+zd1Ye6e77GgLN/z+IJU1q23XTvJzix7x+48m+/3dK2/elTmbVmJSOu/iWOcUw33XQVNbH4QtUWbMddV+8SHREZkg8k/1porCmPTuAo0ReqHfOqAWFHEIm1ChaGHUGOUCyK+4AF68OOIHJEqqqSPwt5moT0Ifa66Sp6YtIto+kHRIISyXH5UaFuGREJW8nc/6HOprW62h845UVe+tAph430ue3sy5hmda3ay+f+G7VWcdhIn9/2HXvYSJ+flkyh1Bpa2o7q+xf+7qJvU23Vh430abDSw0b6/O8pkxhmG1vaen1oBzdM+SoVVttq27DfE6/tCey/l67cRaRDB67cR16veejzrbH3ZYFduePuoT9KS0s9a1XH+FWDbs9+fxHpVD0Tvb7ePTlkJvmoosp96VIfOmRvS9tE6t3nzfN581pv23RCiS9d2rptIfPc6+tbtZWz1L2qysvLW2/rEyf6woWt25ZS7k0Nm1q1zWOh+8KFPnHiwbahNLmXl3tVVev9C+I9lV+TfJ4loN7bqauxuHL/0o9quL35qvyGEhEJWhT73M1sipm9Zmavm1llUOcBuOMdFXYRkXSBFHcz6w78CLgA+Cgwy8w+GsS5RETkcEFduU8CXnf3N939A+B+YEZA5xIRkUP0COi4w4G30n7fAPx1+gZmVgEtd0bsMrPXsj+dDTLjnez371KDIDJZIVp5lTU4UcobpazkWL9GtvdCUMW9U+5eC+RlZWszq2/vS4VCE6WsEK28yhqcKOWNUlYILm9Q3TJNwIlpv49ItYmISBcIqrj/BhhtZqPM7CjgUiA6qxeLiERcIN0y7r7XzL4CPAx0B+5y91eCOFdKXrp3ukiUskK08iprcKKUN0pZIaC8BXETk4iI5JcmDhMRiSEVdxGRGIp0ce/KKQ5yZWZ3mdkWM3s57CydMbMTzexxM1trZq+Y2dfCztQRM+ttZs+b2W9TeReEnakzZtbdzF4ws2VhZ+mMmTWa2UtmtsbMspwEqmuY2QAze8DM1pnZq2b2sbAztcfMTk19pgce75rZtXk7flT73FNTHPwO+AzJm6R+A8xy97WhBmuHmX0C2AX8xN1PCztPR8xsKDDU3VebWX+gAbiwgD9bA/q6+y4z6wmsAr7m7s+GHK1dZnYdUAYc4+7lYefpiJk1AmXuXvA3BpnZvcBT7v6fqZF6R7v79pBjdSpVz5qAv3b3vCwtF+Ur90hNceDuTwLbws6RCXff5O6rU893Aq+SvOu4IKVmP92V+rVn6lGwVy1mNgKYCvxn2FnixMyOBT4B3Ang7h9EobCnnAe8ka/CDtEu7m1NcVCwBSiqzCwBnA48F3KUDqW6OdYAW4BH3b2Q894G/COwP+QcmXLgETNrSE0bUqhGAc3A3akur/80s75hh8rQpcB9+TxglIu7BMzM+gFLgGvd/d2w83TE3fe5ewnJu6EnmVlBdn2ZWTmwxd0bws5yBM5294kkZ3m9OtXFWIh6ABOB2939dODPQEF/FweQ6j6aDvxPPo8b5eKuKQ4ClOq7XgIsdvcHw86TqdSf4Y8DU0KO0p6zgOmpfuz7gXPN7L/CjdQxd29K/dwCPESyS7QQbQA2pP3V9gDJYl/oLgBWu/vmfB40ysVdUxwEJPUF5Z3Aq+7+g7DzdMbMBpvZgNTzPiS/ZF8Xaqh2uPsN7j7C3RMk/zf7mLt/PuRY7TKzvqkv1Ul1cZwPFOSIL3d/G3jLzE5NNZ0HFOQggEPMIs9dMhDirJC5CmGKg5yY2X3AZGCQmW0Aqtz9znBTtessYDbwUqofG+BGd18RXqQODQXuTY046Ab8t7sX/BDDiBgCPJT8954ewE/dfWW4kTr0VWBx6oLvTeALIefpUOofzM8A8/N+7KgOhRQRkfZFuVtGRETaoeIuIhJDKu4iIjGk4i4iEkMq7iIiMaTiLiISQyruIiIx9P8BNPyt/4EE1+YAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "uncertainties = {\n", " \"sig\": {\n", " \"up\": jnp.vstack((\n", " [180., 60., 35., 25., 15., 6., 6.], # up source 0 (red in the plot below)\n", " [170., 65., 32., 22., 11., 6., 6.])), # up source 1 (blue in the plot below)\n", " \"down\": jnp.vstack((\n", " [130., 40., 25., 15., 10., 4., 4.], # down source 0\n", " [140., 45., 28., 17., 8., 4., 4.])), # down source 1\n", " },\n", " \"bkg\": {\n", " \"up\": jnp.vstack((\n", " [110.]*7, # up source 0\n", " [115]*7)), # up source 1\n", " \"down\": jnp.vstack((\n", " [90.]*7, # down source 0\n", " [80]*7)), # down source 1\n", " },\n", "}\n", "\n", "_,bins,_ = plt.hist(x[:-1], weights=sig_shape, bins=x, label=\"Signal\", histtype=\"step\", lw=2)\n", "plt.hist(x[:-1], weights=bkg_shape, bins=bins, label=\"Background\", histtype=\"step\", lw=2)\n", "colors = [\"red\", \"blue\"]\n", "for i in range(2):\n", " for s in [\"sig\", \"bkg\"]:\n", " for d in [\"up\", \"down\"]:\n", " plt.hist(x[:-1], weights=uncertainties[s][d][i,:], bins=bins, ls=\"--\", histtype=\"step\", color=colors[i])\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As before, we register those shape uncertainties with the signal and the background process. Note that the nuisance parameter is now an array:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "shape_nuis = Variable(\"shapes\", [0.,0.], nuisance=True)\n", "sig_contrib.add_shape_syst(shape_nuis, up=uncertainties[\"sig\"][\"up\"], down=uncertainties[\"sig\"][\"down\"])\n", "bkg_contrib.add_shape_syst(shape_nuis, up=uncertainties[\"bkg\"][\"up\"], down=uncertainties[\"bkg\"][\"down\"])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we define a model and add the signal and background processes to it:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "scrolled": true }, "outputs": [], "source": [ "model = Model()\n", "model.add_proc(sig)\n", "model.add_proc(bkg)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we can \"prepare\" the model, and compile the NLL and gradient functions:" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Compiling pred\n", "Compiling NLL\n", "Compiling NLL grad\n", "Tracing Hessian\n" ] } ], "source": [ "model.prepare()\n", "model.compile()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running the fit\n", "\n", "To go further we need to generate some pseudo-data. We can do that using frequentist toys, or using an \"Asimov\" toy. For that we need to use the default values of the variables, but with a twist: we'll inject a signal strength of 0.05 instead of 1:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [], "source": [ "values = model.values_from_dict({mu: 0.05})\n", "asimov = model.pred(values)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's fit those pseudo-data:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Computing covariance matrix...\n" ] } ], "source": [ "best_fit = model.fit(asimov)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can check out the best-fit value, the NLL at the minimum, and the (analytical!) covariance matrix:" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Best fit: [ 5.00000068e-02 -3.25262830e-09 -9.68462346e-09 -7.50083739e-09]\n", "Minmum NLL: 25.379123724585103\n", "Hessian at the minimum:\n", "[[ 0.00668827 -0.00310521 -0.00939339 -0.0073618 ]\n", " [-0.00310521 0.9506984 -0.10175006 -0.18456131]\n", " [-0.00939339 -0.10175006 0.78990573 -0.38105691]\n", " [-0.0073618 -0.18456131 -0.38105691 0.30866546]]\n" ] } ], "source": [ "print(f\"Best fit: {best_fit.x}\")\n", "print(f\"Minmum NLL: {best_fit.fun}\")\n", "print(f\"Hessian at the minimum:\\n{best_fit.hess_inv}\")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can now obtain the profiled (\"minos\") 1-sigma bounds on mu, and compare them with the 1-sigma intervals from the covariance matrix. For that we need to use the \"index\" of the variable mu among the returned arrays:" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Best 𝜇: 0.05 +- 0.08\n", "Profiled best 𝜇: 0.05 +0.09 -0.05\n", "Stat.-only uncertainties on 𝜇: +0.07 -0.05\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "/home/swertz/Documents/Recherche/Stats/Tests/smoofit_venv/lib/python3.8/site-packages/scipy/optimize/_hessian_update_strategy.py:182: UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.\n", " warn('delta_grad == 0.0. Check if the approximated '\n" ] } ], "source": [ "mu_idx = mu.sub_idxs[0]\n", "best_mu = best_fit.x[mu_idx]\n", "print(f\"Best 𝜇: {best_fit.x[mu_idx]:.2f} +- {jnp.sqrt(best_fit.hess_inv[mu_idx,mu_idx]):.2f}\")\n", "up,down,_,_ = model.minos_bounds(mu, best_fit, asimov)\n", "print(f\"Profiled best 𝜇: {best_mu:.2f} +{up-best_mu:.2f} -{best_mu-down:.2f}\")\n", "\n", "direction = np.zeros((len(best_fit.x,))) # using np for convenience because jax arrays are immutable\n", "direction[mu_idx] = up-best_mu\n", "up_stat = model.nll_crossings(direction, best_fit, asimov).x[mu_idx]\n", "down_stat = model.nll_crossings(-direction, best_fit, asimov).x[mu_idx]\n", "print(f\"Stat.-only uncertainties on 𝜇: +{up_stat-best_mu:.2f} -{best_mu-down_stat:.2f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, let's do a profile NLL plot and compare it with a simple NLL scan of $\\mu$ (showing the effect of the nuisance parameters), and also superimpose the Minos bounds found just above:" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABEUElEQVR4nO3dd3hUZfbA8e9Jp7eEIgECEnoJEHoRBenYUBBFxYaKurr+LKurYl/XtWBFURFFRVyxgIJSBKRD6BBCb6Em1EAS0t7fH3dgEwjJTDJzZzI5n+e5T5Jb3vcwGXLm3reJMQallFLqnABvB6CUUsq3aGJQSimVhyYGpZRSeWhiUEoplYcmBqWUUnkEeTuA4goPDzdRUVHeDkMpv3A2M4etR1KoXbkMVcuFeDuci2VnwuFNUL46VLzM29GUaKtWrUo2xkTkd6zEJ4aoqCji4uK8HYZSfuHjBTt4fWYCS57uRc1KYd4O52IL3oB5r8LfFkHVBt6OpkQTkT2XOqaPkpRS5/25+QjNalX0zaSQkw2rJ0GDnpoUPEwTg1IKgBOpGazae5xeTat7O5T87ZgHJ/dC2zu8HYnf08SglAJgwdYksnMMVzXx0cSweiKUrQZNBno7Er9X4tsYlFLu8WfCEaqVC6F1ZGVvh3KxlMOwZSZ0egCCQs/vzszMJDExkfT0dC8G59vCwsKIjIwkODjY6Ws0MSilyMrOYcHWJHo1qUFAgHg7nIutmQQ5WRc9RkpMTKRChQpERUUh4oNxe5kxhqNHj5KYmEj9+vWdvk4fJSmlWLPvBCdSM32zfSEnG1ZNhPo9IDw6z6H09HSqVaumSeESRIRq1aq5fEeliUEpxdzNRwgKELpFh3s7lIttnwMn90HsXfke1qRQsKK8PpoYlFLMSzhCh/pVqRjm/HNo26z8HMrXgCaDvB1JqaGJQalSLvF4KlsOp/hmb6Tje2DbLGh7OwT6YNICAgMDiYmJOb/t3r2bLl26uKXsqKgokpOT3VKWK7TxWalSbl7CEQDfTAyrvwQRnx67UKZMGdauXZtn35IlS7wTjJvYdscgIhNE5IiIbCzkvPYikiUiN9oVm1Kl2dyEI9QPL0eDiPLeDiWvrAxY/RVE94XKdbwdjUvKl7dey59++olevXphjOHgwYM0atSIQ4cOkZSUxJAhQ2jfvj3t27dn8eLFABw9epQ+ffrQvHlz7rnnHry1wqaddwwTgQ+Ary51gogEAv8GZtkUk1KlWmpGFkt2HOW2TvW8HcrFEn6FM0nQ/m6nTn9x+ibiD5xyawjNLqvImMHNCzwnLS2NmJgYAOrXr89PP/10/tj111/P1KlT+fDDD/n999958cUXqVmzJrfccgt///vf6datG3v37qVv375s3ryZF198kW7duvH888/z22+/8fnnn7v13+Ms2xKDMeYvEYkq5LSHgalAe89HpJRasv0oGVk5vvkYKW4CVK4Ll1/l7UgKlN+jpNzef/99WrRoQadOnRg+fDgAc+bMIT4+/vw5p06d4vTp0/z111/8+OOPAAwcOJAqVap4NPZL8Zk2BhGpDVwPXEkhiUFERgGjAOrWrev54JTyU3MTjlA+NIj2UVW9HUpeSVtg90LoNQYCAp26pLBP9t6SmJhIQEAAhw8fJicnh4CAAHJycli2bBlhYT44WSG+1StpLPCUMSansBONMeONMbHGmNiIiHynE1dKFcIYw7yEI3SPDickyJf+FABxX0BAMLS5zduRFEtWVhZ33XUXkydPpmnTprz99tsA9OnTh/fff//8eefuOHr06MG3334LwMyZMzl+/LjtMYMP3TEAscB3jsEY4cAAEckyxvzs1aiU8lPxB09x6FS67z1GykiFdd9Cs2ugfMn+4Pfaa6/RvXt3unXrRuvWrWnfvj0DBw7kvffe48EHH6RVq1ZkZWXRo0cPPv74Y8aMGcPw4cNp3rw5Xbp08doTEZ9JDMaY8xN5iMhE4FdNCkp5zp+bjyACPRv7WGLYOBXST0Ksc43O3nb69OlL7nv++efP76tQoQIJCQnnf54yZcpF11WrVo1Zs7zf98a2xCAik4GeQLiIJAJjgGAAY8zHdsWhlLLMSThCq8jKRFQILfxkO8VNgIgmUM89g8SU6+zslTTchXNHejAUpUq9w6fSWbfvBI/3aeTtUPI6sAYOrIb+b1gD25RX+FiLk1LKDnM2HwagT/OaXo7kAis/g+Cy0GqYtyMp1TQxKFUKzY4/TL1qZYmu7kOjnVOPwYYfoNVQKFPZ29GUapoYlCplTp/NYsn2o/RpVsO3pqxeMwmy0qHDKG9HUuppYlCqlFmwJYmM7ByubuZDj5Fysq3HSPW6QQ3fHKhWmmhiUKqUmRV/iKrlQmhXzzvTLeRr2yw4sRc63OvtSFx2btrtFi1acNNNN5GamurS9cOHD6dVq1a88847PP/888yZMweAnj17EhcX53Q58+fPZ9Ag96xZ4TPjGJRSnpeZncOfCUfo17wmgb60tvPyT6DCZdBkoLcjcVnuuZJuvfVWPv74Yx577LHzx7OysggKyv9P7aFDh1i5ciXbt2+3I1Sn6R2DUqXI8p3HSEnP4upmNbwdyv8kbYWd86D9XT67GI+zunfvzvbt25k/fz7du3fnmmuuoVmzZqSnp3PnnXfSsmVL2rRpw7x58wBraoz9+/cTExPDwoULGTlyJD/88MNF5c6aNYvOnTvTtm1bbrrppvMD6H7//XeaNGlC27Ztz0++5w56x6BUKTI7/hBhwQF0j/ahqSZWfgaBIdB2ZPHKmfkPOLTBLSGdV7Ml9H/dqVOzsrKYOXMm/fr1A2D16tVs3LiR+vXr89ZbbyEibNiwgYSEBPr06cPWrVuZNm0agwYNOn/Hkd8028nJybzyyivMmTOHcuXK8e9//5u3336bJ598knvvvZc///yThg0bMmyY+7r46h2DUqWEMYbZ8YfpHh1BmRDnZiz1uLMpsPZbaH59iZ0X6dx6DLGxsdStW5e777am8ujQoQP161sz/SxatIgRI0YA0KRJE+rVq8fWrVudKn/ZsmXEx8fTtWtXYmJi+PLLL9mzZw8JCQnUr1+f6OhoROR8+e6gdwxKlRKbDpziwMl0Hr3ah0Y7r/sOMlLc00XVyU/27nap9RjKlSvnlvKNMVx99dVMnjw5z/6C1oAoLr1jUKqUmBV/mACBXr4ym6oxsOJTuKwN1G7n7Wg8qnv37nzzzTcAbN26lb1799K4cWOnru3UqROLFy8+30B95swZtm7dSpMmTdi9ezc7duwAuChxFIcmBqVKiVmbDhFbryrVyvvIpHm7FkDyFuhwn9/PizR69GhycnJo2bIlw4YNY+LEiYSGOvd7iIiIYOLEiee7tXbu3JmEhATCwsIYP348AwcOpG3btlSv7r6EL95abNpdYmNjjSt9fZUqjfYdS6X7G/N4dmBT7unewNvhWL67FfYuhb/HQ3DRVjLbvHkzTZs2dXNg/ie/10lEVhljYvM7X+8YlCoFZsVbk+b5TDfVE3thywxoe0eRk4LyHE0MSpUCs+MP0ahGeepVc0+DaLHFTbC+xt7l3ThUvjQxKOXnjp/JYMWuY/TxlbmRMlJh1URoPAAq1yl2cSX9cbinFeX10cSglJ/7M+EIOcaHHiOtnwJpx6Hzg8UuKiwsjKNHj2pyuARjDEePHiUszLXHdTqOQSk/Nzv+MDUrhtGydiVvh2J1UV02Dmq1hrqdi11cZGQkiYmJJCUluSE4/xQWFkZkZKRL12hiUMqPpWdms2BrEkPa1SbAFybN2/Gn1UX1+k/c0kU1ODj4/Ohi5T76KEkpPzZ/SxJpmdn0a17L26FYlo2D8jWsKTCUz7ItMYjIBBE5IiIbL3H8VhFZLyIbRGSJiLS2Kzal/NXvGw9SpWwwHRtU9XYo1iyq22dD+3sgyEcG2al82XnHMBHoV8DxXcAVxpiWwMvAeDuCUspfnc3KZs7mI1zdrAbBgT7wcGD5xxAYql1USwDb2hiMMX+JSFQBx5fk+nEZ4FpriVIqj0Xbkjl9Nov+LX3gMVLqMVg3GVoNhXLh3o5GFcIHPkbk625g5qUOisgoEYkTkTjtjaBU/mZuPESFsCC6Xu4Df4hXfwmZqdDpAW9Hopzgc4lBRK7ESgxPXeocY8x4Y0ysMSY2IqJkzuGulCdlZucwO/4wVzetQUiQl/+bZ2das6jWvwJqNPduLMopPpUYRKQV8BlwrTHmqLfjUaqkWrrjKCfTMn3jMdLmaXBqP3Qa7e1IlJN8JjGISF3gR+A2Y4xzSxsppfI1c+NByoUE0j3aBx4jLRsHVRtAdB9vR6KcZFvjs4hMBnoC4SKSCIwBggGMMR8DzwPVgI/EGviSdakpYZVSl5aVncMfmw5zVdMahAV7eQnPfSshcSX0/w8E+MznUFUIO3slDS/k+D3APTaFo5TfWrH7GMfOZDCghQ9Mmrd8HIRWgphbvB2JcoGmcKX8zMwNhygTHEjPxl5ewvPEPtj0M7S9DULLezcW5RJNDEr5kZwcw++bDtGzcQRlQrz8GGn5x9Z8SNpFtcTRxKCUH1m19zhJKWe93xsp7YS15kKLIVBJx6qWNJoYlPIjMzYcJCQogKuaePkx0qqJkHEaOj/k3ThUkWhiUMpP5OQYft94iB7REZQP9eKM+lkZ1mOkBj2hVivvxaGKTBODUn5iXeIJDp5MZ0BLL/dG2vgDpByELg97Nw5VZJoYlPITMzceIjhQ6NXUi0t4GgNL3ofqzeHyXt6LQxWLJgal/IAxhpkbD9K1YTiVygR7L5Adc+FIvHW34IYV2pR3aGJQyg9s3H+KfcfS6O/tQW1L3ocKtazeSKrE0sSglB+Yvv4AwYFC3+ZeTAwH18HO+dDxfggK8V4cqtg0MShVwuXkGH5dd4Ae0RFULuvFP8hLPoCQ8tBupPdiUG6hiUGpEm713uMcOJnO4NaXeS+Ik4mwcaqVFMpU9l4cyi00MShVwk1bd4DQoAB6N/Nib6Rl46yvHe/3XgzKbTQxKFWCZWXnMGPDQXo3reG9QW1pJ2DVl9DiBqhcxzsxKLfSxKBUCbZs5zGST2cwuLUX50aK+xwyUnRAmx/RxKBUCTZt3X7KhwZ5b4rtzDTrMVLD3lCrtXdiUG6niUGpEupsVja/bzxEn+ZeXKltzddwJgm6Pead+pVHaGJQqoRauDWZU+lZ3uuNlJ0Ji9+DOh2hXhfvxKA8QhODUiXUtHUHqFI2mG4Nw70TwMapcHIvdPu7Tn/hZzQxKFUCpWZkMTv+MP1b1iI40Av/jXNyYNE7UL0ZRPe1v37lUS6/o0SknIi4/EBTRCaIyBER2XiJ4yIi74nIdhFZLyJtXa1DqdJi7uYjpGVmM7iVlx4jbZ0JSQnW3UKAfr70N4X+RkUkQERuEZHfROQIkAAcFJF4EfmPiDR0sq6JQL8CjvcHoh3bKGCck+UqVepMX3eA6hVC6VC/qv2VGwML34bK9aD5DfbXrzzOmREx84A5wNPARmNMDoCIVAWuBP4tIj8ZY74uqBBjzF8iElXAKdcCXxljDLBMRCqLSC1jzEFn/iFF0bNnz4v2DR06lNGjR5OamsqAAQMuOj5y5EhGjhxJcnIyN95440XHH3jgAYYNG8a+ffu47bbbLjr+f//3fwwePJgtW7Zw3333XXT82WefpXfv3qxdu5ZHH330ouOvvfYaXbp0YcmSJTzzzDMXHR87diwxMTHMmTOHV1555aLjn3zyCY0bN2b69Om89dZbFx2fNGkSderUYcqUKYwbd3Fu/uGHHwgPD2fixIlMnDjxouMzZsygbNmyfPTRR3z//fcXHZ8/fz4Ab775Jr/++mueY2XKlGHmzJkAvPzyy8ydOzfP8WrVqjF16lQAnn76aZYuXZrneGRkJF9/bb0NH330UdauXZvneKNGjRg/fjwAo0aNYuvWrXmOx8TEMHbsWABGjBhBYmJinuOdO3fmX//6FwBDhgzh6NGjeY736tWL5557DoD+/fuTlpaW5/igQYN4/PHHgeK993YlHuTr5+6iRsUwes0ue/64be+9qeN45tUFUO1y+Ln3+eP63rP/vXfu3+RuztwD9jbGvGyMWX8uKQAYY44ZY6YaY4YAU9wQS21gX66fEx37LiIio0QkTkTikpKSilRZz9ZrOb3v4jevUr7uz4QjpB/+nOQ9n3il/uYnnmHsvTlQ3otTcCiPEusDuk2VWXcMvxpjWuRz7FfgdWPMIsfPc4GnjDFxBZUZGxtr4uIKPCVfQ65YBsDUBZ1cvlYpb7p9wgpmvB3FFY0jmDrV5t5AB9bA9A5QvSlcv97eupVbicgqY0xsfsdcnlxFRG4BrgGyAQGmG2MmFy9EAPYDuSdaiXTs8whNCKokOnr6LIu3J/PEmxV5qp8XRjsvegeO1YQ7Ftpft7JNUboTXGGMudkYc6sx5hagm5timQbc7uid1Ak46cn2BaVKohkbDpKdY7zTGyl5G8RPg/b3QFgl++tXtinKdIyhIjIQqz0gEijjzEUiMhnoCYSLSCIwBggGMMZ8DMwABgDbgVTgziLE5rSn75oPwL8m9PRkNUq51Y9r9tO4RgW+eq8CIuBoj7THwrchKAwqHoW1T0OMnZUrOxUlMYwGbgBaYjUQP+TMRcaY4YUcN8CDRYinSJauqmxXVUq5xc6k06zZe4Kn+zdh8hib2xaO7YT1U6z1Fk4tsbduZTuXE4MxJhUosGuqUsr9fl6znwCB69rUxh2Nei5Z+DYEBEHXv8FyTQz+rlhDFkVksbsCUUpdWk6O4cc1++naMJwaFcPsrfz4Hlg32Vq2s0JNe+tWXlHcsexeXGRWqdJj5e5jJB5P44a2+Q7t8axF74AEQNdH7K9beUWhj5JE5H1gg2PbaIxJyXXYvkEQbhZZ84y3Q1DKaT+u3k/ZkED6Nrc+sUdG2lTxyURrzYW2t0MlR1Iqa1flylucaWPYgNXQfCvQQkROOfZtBCp4MDaP+vqPrt4OQSmnpGdmM2PDQfq3qEXZEOu/7Nd2tfItGmt97fb3/+3rok2M/q7QxGCMGZ/7ZxGJxEoUrYA/PBSXUsphdvxhUs5m2f8Y6dQBWP0lxNwClesUfr7yG073ShKRC9fuywTmiUiMMWatW6OywaO3LABg7LdXeDkSpQr24+pEalUKo1ODauf3nZvjzjHvmmcsfg9ysqH7Bf/1Vzkqb+fJypU3udJdNdaxTXf8PAhYD9wvIv81xrzh7uA8ae0mHbmpfF9Syln+2pbMqB4NCAz439iFCybvdL+Uw7DqC2g9HKpE5T123NOVK29zJTFEAm2NMacBRGQM8BvQA1gFlKjEoFRJMG3dAbJzDDe0sfkx0pL3IDvj4rsFVSq40l21OnA218+ZQA1jTNoF+5VSbvLj6kRa1q5EdA0b+3mcToK4CdByqLXmgip1XLlj+AZYLiK/YM2qOgj4VkTKAfGeCE6p0mzLoRQ2HTjFmMHN7K146fuQmQY9Hre3XuUznE4MxpiXRWQmcK6f5/251kq41e2ReVij+qe8HYJSBfpxTSKBAcLg1hePI23UyEOVphyG5eOh5U0QHp3/ORU8VbnyFa7OlZQJ5GANbMt0fzj2Gf9zD2+HoNQlZecYfl6zn56NIggvH3rR8fHj87nIHRa9Y7Ut9PzHpc/p6KnKla9wuo1BRB7BepwUjtXe8LWIPOypwJQqzZbsSObwqbPc0NbGUcYn91ttCzHDtW2hlHPljuFuoKMx5gyAiPwbWAq874nAPG3UdX8BeuegfNPUVYlUCAuiV9P8V2kbNcr66tY7h4VvgsmBHk8WfN5yR+V65+C3XEkMgrWc5znnlvYskbbuqujtEJTK18nUTGZuPMRNsZGEBQfme87WrW6u9PhuWD3JmhOpSr2Cz01xd+XK17iSGL7A6pX0k+Pn64DP3R6RUqXcL+v2czYrh5vb17Wv0gX/sWZQ1Z5ICtd6Jb0tIgv4X6+kO40xazwTllKlkzGGySv20fyyirSobdPo/OTt1noLHe+HijqTvnKxV5IxZhXWKGellAds3H+KzQdP8fK1ze2rdMHrEBSadwZVVao5sx5DCvmvuyBYSzWXyIf1Mc1PejsEpS7y3cq9hAYFcE1MwVNgxMS4qcIjm2HDD9DtUSgf4dw1VdxVufJVzky77bax+CLSD3gXCAQ+M8a8fsHxusCXQGXHOf8wxsxwV/256ayqytekZWQzbe0BBrasRaUywQWe67ZZVee9BiHlocvfnL9GZ1X1e8Vd2tNpIhIIfAj0B5oBw0XkwrH+zwLfG2PaADcDH9kVn1LeNmPDQVLOZjG0vU1rHxxcB5unQecHoWxVe+pUJYJtiQHoAGw3xuw0xmQA3wHXXnCOAc49mqoEHPBUMCP6LmZE38WeKl4pl01ZuY+oamXpWL/wP9IjRlhbsfz5KoRVhs6jXbtuyQhrU37L5cQgIoOLWFdtYF+unxMd+3J7ARghIonADCDfkdUiMkpE4kQkLikpqUjBJB4qR+KhckW6Vil325F0mhW7jzGsfV1ECh8elJhobUW2Zwls+wO6PgJhLvZ+Sk20NuW3inLH8Krbo/if4cBEY0wkMACYJCIXxWiMGW+MiTXGxEZEONlgppQP+z5uH4EBwpB2Nqy7YAzMHgMVLoNOD3i+PlXiFCUxFHW0834g98PTSMe+3O4GvgcwxiwFwrDmZlLKb2Vm5zB1VSK9mlSneoUwz1eY8BskrrAmygsu4/n6VIlTlMSQX9dVZ6wEokWkvoiEYDUuT7vgnL1ALwARaYqVGIr2rEipEmLu5iMkn87g5g42NDpnZ8HcFyG8EcSUuNnylU1cnXa7yIwxWSLyEPAHVlfUCcaYTSLyEhBnjJkG/B/wqYj8HSsBjTTGFDURFahzuxOeKFYpl01ZuZeaFcPoEe38Y9HOnYtY2dpvIHkrDPsGAov43z+8qJWrkkJc/bsrIuuNMa08FI/LYmNjTVxcXOEnKuWDDp5Mo+vrfzK6Z0Me79vYs5VlpML7baFSHbh7FjjRyK38l4isMsbE5nesKB8ZDhczHqWUww9xieQYGBprw2OkFZ9AykEY8rkmBVUgl9sYjDFXeyIQuw25YhlDrljm7TBUKZadY/hu5T66NqxG3WplXbp2yBBrc1rqMVj4DjTqB1FdCz+/IAuHWJvyW7a1Mfiaoyds6P2hVAHmJRxh/4k0nh3Y1OVrjx518YJFb8PZU9DreZfrushZVytXJY2dI5+VUrl8tWwPNSqG0rtZDc9WdGIfLB8PrYdDDRtnbVUlliYGpbxgV/IZ/tqaxC0d6hEc6OH/hvP/ZX298hnP1qP8RrHekSKikw0pVQRfL9tDUIAw3NNjFw5tgLXfQod7obJNk/OpEq+4bQwldrmnXt1PeDsEVUqlZWTz37h99GtRk+oVi9bW1auXEycZA388A2Uqu3fJzhrOVK5KMmcW6nkf2ODYNhpjUnId9sjgMzs890FPb4egSqlf1u7nVHoWt3eOKnIZzz3nxElbf4ddf0H/N6BMlSLXdZGWzlSuSjJn7hg2AC2BW4EWInLKsW8j4LZFfJQqDYwxfLV0D01qVqB9lBv/WF8oKwNmPWtNfRF7l+fqUX7JmRXcxuf+WUQisRJFK6zpLUqk/h1XAjBzeXsvR6JKk9V7jxN/8BSvXt/Cqem1L6V/f+vrzJmXOCHuczi6HW75HgILXg3OZfMclV95qcpVSed0G4OIPHbBrkxgnojEGGPWujUqG6Slu/k/i1JO+GrpHiqEBnFdIWs6FyYtrYCDqcdg/uvQoCdE9ylWPfnKLqhy5Q9c6ZUUC9yPtbhObeA+oB/WpHdPeiA2pfxKUspZZmw4yJB2kZQL9eDY0gVvWIPZ+r6mU1+oInHl3RkJtDXGnAYQkTHAb0APYBXwhvvDU8p/TFm5l8xsw4hO9TxXSfI2WPkptL1dB7OpInPljqE6cDbXz5lADWNM2gX7lVIXyMrO4Zvle+nasBoNq5f3XEWznoOgMnDlPz1Xh/J7rtwxfAMsF5FfsFZxGwR8KyLlgHhPBOdJg/qc8HYIqhSZm3CEgyfTGTPYPZ/iBw3KZ+fO+bB1JvR+AcpXd0s9+aqdX+XKn7i0HoOIxALnpmZcbIzx+kIIRV2PYcn2ZP41M4Gv7+lIpTLaEK08a8Rny9mZdJq/nrySIE9MgZGTDZ/0sNoWHlwJwTpJpCpYQesxuPoOzQRygGzH9yVWxTLBbNh/kqmrEr0divJz2w6nsGh7Mrd0rOuZpACwaiIc3gi9X9SkoIrN6XepiDyC9TgpHKu94WsRedhTgXnaQwN2EfFLDb5etgcPrR6qFACfLdxFWHAAt3R0X6Nzz57WBljdU/98GaK6Q/Pr3VbHJc3paW3Kb7ny8eVuoKMxZowx5nmgE3CvZ8KyR0hgADuTz7Bkh84vrzwjKeUsP63Zz5C2kVQtF+KZSua+COmnrKkvtHuqcgNXEoNgPUI6J9uxr8QKDhSqlgvhq6W7vR2K8lOTlu4mMyeHu7vV90wF+1fDqi+h4/1Qo5ln6lCljiuJ4QusXkkviMgLwDLgc1cqE5F+IrJFRLaLyD8ucc5QEYkXkU0i8q0r5btOGBpbh9nxhzl4UkdzKvdKy8hm0rI99GpSgwYRnuiiamDGE1AuAno+5YHyVWnldGIwxrwN3AUcc2x3GmPGOnu9iAQCHwL9gWbAcBFpdsE50cDTQFdjTHPgUWfLL6pbO9bFAJOX7/V0VaqUmbo6keOpmdzb3UN3C6cPw/446PMyhFXyTB2qVHJpXL4xZhXWKOei6ABsN8bsBBCR74BryTsG4l7gQ2PMcUd9R4pYV6GGXncSgDpVy3Jl4+pMXrmPh66KJiRIF7VTxZeTY5iwaBetIivRoX5Vt5c/9LpUWPAR1OkErYa5vfwC1R1qb33Kds6sx5BC/usuCGCMMRWdrKs2sC/Xz4lAxwvOaeSoczEQCLxgjPk9n5hGAaMA6tat62T1eY1+8Yrz39/WqR53TlzJrPhDDGpVYtceUj5kbsIRdiaf4b3hbYo1i+qljI5+Hk5+DgMW2N/g3Gi0vfUp2xX68dgYU8EYUzGfrYILScFZQUA00BMYjjVBX+V8YhpvjIk1xsRGREQUqaLUU6mknkoF4IpGEdSpWoZJS/cUNW6l8vh04U5qVy7DgBY13V/4wfWkLvmW1Bb3Q61W7i+/MFmp1qb8VqGJQZz4uOPMOcB+IPeis5GOfbklAtOMMZnGmF3AVqxE4XYDum9lQPetAAQECCM61mP5rmNsOZRSyJVKFWx94glW7DrGnV2j3D+gzVgNzgO++5EB777k3rKdNX+AtSm/5cy7dp6IPCwieZ7ZiEiIiFwlIl8CdzhRzkogWkTqi0gIcDMw7YJzfsa6W0BEwrEeLe10ouxiuym2DiFBAXy9TO8aVPF8unAXFUKDGNa+TuEnu2rdZNi3DKpEQYBO5aI8w5nE0A9rzMJkETng6Eq6E9iG9bhnrDFmYmGFGGOygIewVn3bDHxvjNkkIi+JyDWO0/4AjopIPDAPeMIYY8vos6rlQhjUqhY/rk7k9NksO6pUfmj/iTRmbDjI8I51qRDm5j/cZ47CH/+EyA5QvoZ7y1YqF2eW9kwHPgI+EpFgrCkx0owxJ1ytzBgzA5hxwb7nc31vgMccm+1u61SPH1fv56c1+7nNk3PmK7/1xaJdCDCyS5T7C5/1T2uSvMHvwtclemyp8nFOPQAVkSYi0gsINcYcPJcURKSfJ4OzW0ydyrSsXYmvl+r8Scp1p9Iz+W7lPga2qsVllcu4t/Cd863HSF0f0RHOyuOc6a76N+BBrMc/n4vII8aYXxyHXwMu6k5aEoy89fRF+0SE2zrV48mp61m+6xidGlTzQmSqpPp62R5On83i3u4N3FtwZhr8+neo2gB6PAHAyJHurcIlDbxZubKDMwPc7gXaGWNOi0gU8IOIRBlj3qUEz5U08slu+e4f3Poy/jVzM58t3KWJQTktNSOLzxbu4srGEbSo7eZRyH/9B47thNt/gWDrTkQTg/IkZx4lBZxb59kYsxur11B/EXmbEpwYkhOPkpx4cbt2mZBAbuscxZzNh9mRdPFdhVL5+Xb5Xo6dyeChq9zcu/pwPCx+F1oPhwY9z+9OTrY2r0hPtjblt5xJDIdFJObcD44kMQirEbqlh+LyuBsH7uPGgfvyPXZ753qEBgXw2UJbesqqEi49M5tP/tpJ14bVaFevivsKzsmB6Y9AaEXo82qeQzfeaG1esehGa1N+y5nEcDtwKPcOY0yWMeZ2oIdHovKy8PKhDGkXydTV+0lKOevtcJSP+z5uH0kpZ3nY3XcLqyZA4gro+xqU08eayj7OTImRaIw5nxhEZHCuY4s9FZi33du9AZnZObpWgypQRlYOH8/fQfuoKnR052R5pw7CnBeh/hXQ+mb3lauUE4oyXv/Vwk8p+eqHl6NPsxpMWraH1Awd8KbyN3V1IgdOpvPwVdHumyzPGJjxOGRnwKB3dFU2ZbuiJIZS8y4d1aMBJ1Iz+X5l/m0RqnTLzM7ho/nbaV2nMt2jw91X8MapkPAr9Hwaql3uvnKVcpJL6zE4+MXIrwfuLXx2yHb1qtKuXhU+W7SLEZ3quX9CNFWiTVt7gH3H0hgzqLn77hZSDlt3C5HtocvDlzztgQfcU12RRHuzcmWHoiQGvzDsoS5OnTeqRwPum7SK3zfpWg3qf7JzDB/O206zWhXp1bS6ewo1xhrIlpEK134EAYGXPHWYzWvz5FHPm5UrO5Taj8D7thxg35YDhZ53ddMa1A8vx/i/duo0Geq83zYcZGfyGR6+qqH77hY2/Be2/Aa9noOIRgWeum+ftXnFmX3WpvxWURLDYbdH4QW3DT3CbUMLXzk0IEC4p3t91ieeZNnOYzZEpnxdTo7hgz+3EV29PH2bu2khnpRDMOMJqNMROhW+Qtptt1mbVyy9zdqU33I5MRhjrvZEIL5sSNtIqpUL4VMd8KaAWfGH2Hr4NA9d1ZCAADfcLRgD0x+FrPRCHyEpZYdS+yjJFWHBgdzRJYo/E46w7bCu8FaaZecY3pq1lQYR5dzX5rTuO9g6E3o9D+EN3VOmUsWgicFJt3WqR1hwAOPm7/B2KMqLflydyLYjp3miT2MC3XG3cOog/P4U1OkEHe8vfnlKuYEmBidVKRfCbZ3q8fPa/ezUyfVKpfTMbMbO2UaryEr0a+GGtgVjYPrfICsDrtNHSMp3lNruqv/3aIbL19x3xeV8vWwv783dxtib23ggKuXLvlm+l/0n0njjxlbu6Ym04lPYNgv6v+HyQLb/+7/iV19kTbxZubJDqU0Mg+/s4PI14eVDub1zPT5duJOHroqmYfXyHohM+aKU9Ew+nLedbg3D6drQDaOcD8fDrGchug90GOXy5YMHF36Ox0R6s3Jlh1L7KGnLqp1sWeV6L6NRPRoQFhzIe3O3eSAq5as+W7iLY2cyeLJf4+IXlpkOU++GsIpWL6Qi3H1s2WJtXnFqi7Upv2VrYhCRfiKyRUS2i8g/CjhviIgYEYn1VCz33XWK++465fJ11cqHckeXKKavP8BW7aFUKiSfPstnC3cysGUtWkVWLn6Bc8bAkXi4bhyUjyhSEffdZ21eseI+a1N+y7bEICKBwIdAf6AZMFxELlrVXEQqAI8Ay+2KzVWjujegbHAg7+pdQ6nwwZ/bSc/K4bE+BY9Gdsq22bD8Y6sHUnSpGxKkSgg77xg6ANuNMTuNMRnAd8C1+Zz3MvBvIN2z4RjIOmutkuWiKuVCuLNrfX5bf5CEQ67fdaiSY9+xVL5ZvoehsZFcHlHMNqXTR+DnB6B6c+j9onsCVMoD7EwMtYHcE6wkOvadJyJtgTrGmN8KKkhERolInIjEJSUlFS2a7ExrpOnKT4t0+T3d61MhNIh35+hdgz97Z85WAkT4W69irs5mDPzyIKSfgiGfQXCYewJUygN8pvFZRAKAt4FC+8IZY8YbY2KNMbEREUV7RktgCAQEwewxcNT1QWuVy4ZwZ7f6zNx4iE0HThYtBuXTthxK4ac1+xnZJYpalcoUr7AV462uqX1egRoXPUFVyqfY2V11P1An18+Rjn3nVABaAPMdfcRrAtNE5BpjTJy7g3n2nzmQngkHQuDn0XDnDJcHGN3drT5fLN7Fu3O2Mf52j7WTKy95feZmyocG8UDPYi6Wc2AtzHoOovtCh3vdEtuzz7qlmKJp4c3KlR3sTAwrgWgRqY+VEG4Gbjl30BhzEjjfQVxE5gOPeyIpAPQe2tb6Zt1/4KdRsGwcdHnIpTIqlQnmnm4NeGfOVjbuP0mL2pU8EKnyhj8TDjNvSxL/HNCUymVDil5Q2nH4/nYoF26NbnbTFN29e7ulmKKp6c3KlR1se5RkjMkCHgL+ADYD3xtjNonISyJyjV1xnLP2ry2s/WsLtBoKjQfC3JcgaavL5dzZLYqKYUG8Pdv1a5VvOpuVzUvT47k8ohx3dIkqekE5Odbd6Kn9cNNEKzm4ydq11uYVx9dam/JbtrYxGGNmGGMaGWMuN8a86tj3vDFmWj7n9vTU3QLAow+n8ejDadYnuEHvQEhZ+Pl+yM5yqZyKYcE80LMhfyYcYfH2ZA9Fq+w0YdFudh9NZczg5oQEFeO/yJL3YMsMq12hjusj7Qvy6KPW5hWrHrU25bd8pvHZqyrUgAFvwv5VsPR9ly+/s2sUkVXK8PKv8WTn6CpvJdnhU+m8/+c2rm5Wgx6NitixAWD3Ipj7IjS7TmdNVSWOJoZzWgyBptfAvNeseWxcEBYcyD/6NyHhUArfx+mShyXZ6zMTyMoxPDewGD2HUg7Bf++Eqg3gmvfd1q6glF00MZwjAgPfhtAK1iCk7EyXLh/Yshax9arw1qwtpKS7dq3yDXG7j/HTmv3c16MBdauVLVoh2Vnww11wNgWGTrLmQ1KqhNHEkFv5CKu94eBaWPBvly4VEZ4b1Izk0xm6mE8JlJ1jGDNtE7UqhRWve+qfL8OexTB4rI5XUCVWqZ12+7XXLjFmodm1EDMC/noTorpBg55Ol9m6TmWub1ObzxbtYniHutSpWsRPncp2U1buY9OBU7w/vA1lQ4r43yL+F1g8FtrdCa1vdmt8F3rtNY8WX7DW3qxc2UGMKdmNpbGxsSYuzs2dlzLOwPgrrT7oDyyG8tWdvvTgyTSufHM+vZvW4INb2ro3LuURJ1Mz6fnmPKJrVGDKqE5FW4Tn4DqY0A9qNIc7ftUpL5TPE5FVxph8R+aW2kdJS37bwJLfNuR/MKQc3PQFnD0FP45yaaK9WpXKMKrH5fy6/iCr9hxzU7TKk96evYWTaZm8MLh50ZJCymGYPBzKVIVh39iSFJYssTavSFpibcpvldrE8Mwz2TzzTPalT6jRHPr/G3bOg8XvuFT2/Vc0oEbFUF76dTM52n3Vp63Ze5xJy/YwolM9ml1WhIbizHSYcqt1dzn8W6vrsw2eecbavGLdM9am/FapTQxOaXuH1Y31z1dhz1KnLysbEsQTfZuwbt8Jpq074MEAVXGczcrmyR/WU7NiGE/0LcLKbMbA9L9B4kq4/hOo1dr9QSrlBZoYCiICg8ZC5brWUoypzj8auqFNbVrWrsS/f0/gzFnXRlMre3zw53a2HTnNaze0pEJYsOsFLB4L66fAlc9CM9tndVHKYzQxFCasotXecPqINe+Nk431AQHCC9c059CpdP7zh66P62s2HTjJR/N3cEPb2vRs7HzngvMSZsCcF607yh6Puz9ApbxIE4MzLmtjzXezdSYs/cDpy9rVq8Ltnerx5dLdrN573IMBKldkZufw5A/rqVI2hOcHFWGswaGNMPUeuCwGrv1QRzYrv1NqxzGMfd/FhVc63mcNXJr9PFRvCg2dm3r4iX5NmB1/mH9MXc+vD3cv3qRsyi3G/7WTTQdO8fGItq5PqX1iL3xzI4RVgpu/heBiLuBTRGPHeqVaSztvVq7sUGr/SsX0aExMDxcaHEXgunFQvRn89y5I3u7UZeVDg3jl+hZsPXyaj+Y7d43ynG2HU3h3zjYGtqxFvxa1XLv4zFGYdANkpsKIqVDxMs8E6YSYGGvziiox1qb8VqlNDHO+X82c71e7dlFoeRg+GQKDYfIwSDvh1GVXNanBNa0v48N529l2OMX1YJVbZOcYnpy6nrKhgbxwTXPXLs44A9/eBCf3wfDvvD7dxZw51uYVh+ZYm/JbpTYxvPJqAK+8WoR/fuW6MGwSHN9jTZaWU8BYiFzGDG5G+dAgnpq6Xqfm9pIvFu9izd4TvDC4OREVQp2/MDsTvr8DDqyBGydAvS6eC9JJr7xibV6x8RVrU36r1CaGYqnXBQa+BTvmWm0OTqhWPpTnBzdj9d4TfL1sj4cDVBfafiSFN2dt4aom1bk2xoVHQMbAtIdh+2xrgsUmAz0XpFI+QhNDUbW7AzrcZ/VSWvONU5dcF1ObHo0ieOP3BPafSPNwgOqc9MxsHvxmDWVDgvjXDS1dm/ZizhhYNxmu/Ce0G+mxGJXyJZoYiqPva9bsq78+CnuXF3q6iPDa9S0wwLM/baCkT2BYUrz0azxbDqfw9tDW1KjowjxGS96Hxe9C+3ugxxOeC1ApH6OJoTgCg+DGL6BSJHw3HJK3FXpJZJWyPN6nMfO2JPHdSl3tzdOmrzvAt8v3ct8VDVwbyLZsHMx61lqas/8bOlZBlSq2JgYR6SciW0Rku4j8I5/jj4lIvIisF5G5IlLPU7F8MqEin0xww+paZavCrT+ABMBX18HJxEIvGdkliu7R4bwwbRNbDmkvJU/Zc/QMT/+4gTZ1K/N4Hxe6Ji/7GH7/BzQdDEM+g4BLrN3hRZ98Ym1e0eETa1N+y7bEICKBwIdAf6AZMFxELuzztwaINca0An4A3vBUPI3bNaBxuwbuKaza5TDiR2ua7knXW/3dCxAQILw9NIYKYcE8+O1qUjN0LiV3O5uVzUPfriFA4P3hbQgOdPKtvnw8/P4UNBlk3Q0GFmEOJRs0bmxtXlGxsbUpv2XnHUMHYLsxZqcxJgP4Drg29wnGmHnGmFTHj8uASE8FM/2LFUz/YoX7CqzVyurffm5k7NmC7wQiKoTy7s0x7Eg6zZhfNrkvDgXAv2duYcP+k/znptZEVnFyJb0Vn8LMJ3w+KQBMn25tXpE43dqU37IzMdQGcj9UT3Tsu5S7gZn5HRCRUSISJyJxSUlJRQrmrbEhvDXWxekQChPVFW6aaK3m9d2tkHW2wNO7NgznoSsb8t9Vify4uvBHUMo5s+MPM2HxLkZ2iaJv85rOXbTiU5jxODQeaCWFIDe/N9zsrbeszSsS3rI25bd8svFZREYAscB/8jtujBlvjIk1xsRGRETYG1xhGveH6z6CXQusqboLGQD3SK9oOkRV5dmfN7Ij6bRNQfqvxOOpPP7fdbSoXZGnBzRx7qKVnzmSwgArsft4UlDK0+xMDPuBOrl+jnTsy0NEegP/BK4xxhT8kdtXtb4Z+r0Om6fD9EcKXBo0KDCAd4fHEBoUwIPfrCY907mR1Opip9IzuWviSnKM4YPhbQkNKqTR2BhY8B/47f+gUX+46UtNCkphb2JYCUSLSH0RCQFuBqblPkFE2gCfYCWFIzbG5n6dHoAeT8KaSTDtIci+dANzrUpleGtoaxIOpfDKb/E2Buk/MrNzePCb1exMOsMnI9oRFV6u4Atysq27hHmvQKubrWlONCkoBdg47bYxJktEHgL+AAKBCcaYTSLyEhBnjJmG9eioPPBfx+jUvcaYkrs01pXPQEAQzH/NmoTthk8v+cfnqiY1GNWjAeP/2kn7qKpcG1NQ84vKzRjDmGmbWLgtmTeGtKJLw/CCL8hMh59GQfwv0OVv0PtFCPDJp6pKeYWt6zEYY2YAMy7Y93yu751b5MANJn1fhFW7XCUCPZ+CkLLWYKnMNBj65SXn8H+8T2PW7j3BE/+11iHu2KCa52P0A58t3MW3y/cyuuflDG1fp+CT009aHQN2L7RGrnd+0J4g3WzSJC9W3tmblSs7SEmfliE2NtbExcV5O4zCxU2AXx+D+t3h5snWFN75OJGawQ3jlpCccpYfR3ehYfUKNgdasvy+8RAPfLOKAS1q8f7wNgQEFDBCOeUQfH0jJCXA9R9DyxvtC1QpHyMiq4wxsfkdK7X3z1M+WMKUD5bYV2HsXdYfo92LrEFwl1jLoXLZEL68swMhQYGM/GIlR1LS7YuxhFm37wSPTllDTJ3KvDW0dcFJ4dBG+PxqOL4Lbv2+xCeFKVOszSv2TLE25bdKbWIY92lZxn3q5MAnd2l9s9Xz5cAa+HIwpBzO97Q6VcsyYWQsR09ncPfEOB0ZnY/E46nc81Uc4eVD+fT2WMKCC+iBtHGqlRSyM2Hkr3D5VfYF6iHjxlmbV2wbZ23Kb5XaxOA1za6xVoE7uh0+vRIOrM33tFaRlfngljZsOnCSh79dQ1b2pbu8ljb7T6Rxy6fLSc/M5ouR7Qkvf4lFd3KyYdZz1oJKNVvBqAVwWRt7g1WqBNLE4A3RV8NdfwACE/pZn2jz0atpDV66tgVzE47wwvRNOk03sO9YKsM+Wcrx1Awm3d2R6BqXaINJPQZfD4El70Hs3XDHdKhQw95glSqhbO2VpHKp1QpGzYMpI6xPtEcSoOfTF3WbHNGpHonH0/h4wQ6qVwjjb72ivRSw9+07lsrN45eRkp7JN/d0pFVk5fxPPLTB6nmUchCueR/a3m5rnEqVdJoYvKl8deuT7G+PwV9vwJF4uP6Ti3osPdm3MUdS0nl79lZSM7J5ql9j11Yh8wN7jp5h+PhlpGZm8+29nWhRu9LFJxkDq7+CmU9Bmcpw50yIzLfThVKqAKU2MfzwWyH93e0SFArXfADVm8Osf8KEvtYo3Kr/mxI8IEB488bWlA0J5OMFO0hJz+Tla1sU3AvHj+xKtpLC2axsvr2nE80uy2cdjdNHYNrfYOtMqN8DbvjMrx8d/fCDFyvv5s3KlR1KbWIIj/ShwWMi0Hk0RDSGH+6Ej7tDv39Bm9vOrxwWECC8fG0LKoQFM27+Dk6fzeLNm1o7v85ACbUj6TS3fLqMrGzD5FGdaFIzn6SQ8JuVFM6mQN9/Qcf7/X4kc3ghg7s9KsyblSs7+Pf/ngJMfGMRE99Y5O0w8mrYCx5YYvWcmfaw9Zz89P+mFRcRnurXhCf7NeaXtQe4f9Iqv550b9G2ZIaMW0J2ziWSwtkU+OVB+O4WqHgZ3LfASrB+nhQAJk60Nq/YOdHalN8qtSOfe7ZeC8D8dTHuDcgdcnJg+TiY8yKEVbQaUBv3z3PK18v28NwvG+lYvyqf3dGe8qH+c/NnjOHThTt5fWYCDauXZ/xtsRdPirdzgZU8T+6Dbn+HK/5RqibB69nT+jp/vhcqn+OovLc3KlfuoiOfS5qAAGsOn1HzoXxNmHyz41HJ/9ZrGNGpHmOHxbBy93GGfryU3clnvBevG6VlZPPId2t5bUYCfZvX5KfRXfMmhRN7Ycpt8NU11jrbd/4OvZ4vVUlBKU/TxODLajSDe+dC10es3jYftId1U86v73BtTG0+uyOW/SfSGPT+In5Ze9HyFiXKvmOpDBm3hOnrD/BE38Z8dGtbyp27E8pMg/mvW6/Bttlw1bMwehnU7ejdoJXyQ5oYfF1QKFz9kjUgrnx1a7roz6+GfSsBuLJxdWY80p0mNSvwyHdreeqH9aRllLx2h8Xbk7nmg0XsO57KhDva8+CVDa0uucZA/DT4oAPM/5e1ytrDcdDjCQgO83bYSvklTQwlRd2OcO88uG4cnEyEz3vD1HvgZCK1K5fhu1GdePDKy/l+1T6u+WARWw6leDtip5xIzeDpH9dz62fLqVY+lGkPdePKJtWthLB9LnwxAL6/DUIrwMjf4KYvoFKkt8NWyq+V2sbn1FOpAJStaPNEeu5w9jQsegeWvG89Z+882uqiWb46i7Yl8+iUtaSkZ/L84GYMb1/XJ8c7GGP4cfV+XpuxmRNpmdzVNYpHezeiXHAAbJkBC9+0JhuscBl0fwza3QmB/tPAXlyp1tuXst54+2Y5Kg8qgf931HkFNT6X2sTgF07shdljYNNPEBgCrYZC5wdJKtOAx75fy8JtyTSrVZGn+jehR3S4z4yW3n4khWd/3siyncdoW7cyr17fkqbVy8KmH2Hh25C0GarUt3obtb7ZepymlHIrTQz5+GjMAgBGv3iFu0OyX/J2WPYhrJ0MWWnQsDc5nR5ieko0b87eyr5jaXS5vBpP9WtC6zqVvRbmsTMZfLZwJ58u3EnZkCD+0b8Jw+pnELBhCqz7Dk7uhYim0P3/oPn1eodQgI8+sr6OHu2Fyrc6Km/kjcqVu2hiyIdPj2MoqjNHrZXiVoyHM0cgoilZzW9gekYsLy/P5tiZDAa2rMXjfRtT/8JxAR4Uf+AUXy7Zzc9r93M2K4dbW1fk6bqbKb/5v5C4wnoc1qAntL8HGvUvFQPUikvHMajiKigx6Ecyf1KuGlzxBHT9G2z4L6yaSND8V7keuCa8CSsu68brWxrTe9NBulwezuDWl9G3WU0qlQ12eyhZ2TnMjj/MF0t2s2LXUaKDk/hX1EGuDtlAhe1zYctZiGgCvV+0HoFVvMztMSilisbWxCAi/YB3gUDgM2PM6xccDwW+AtoBR4FhxpjddsboF4JCoc0Iazu5HxJ+JTB+Gp33fM4vAYbjFSNZcqAxS3fW5bufLyfi8rb0i6lH76Y1qBBW9CRx8GQaK3YdY8WuY2zYvJkGp9dwR5kEPq8UT4WzhyARKF8DYu+02g5qxZyfC0op5TtsSwwiEgh8CFyN9SdipYhMM8bE5zrtbuC4MaahiNwM/BsYZleMfqlSbeh4n7WdPgIJv1Jly+8M2B/HwOy5AGTsCWLzrrpM/6kBGeXrElSpOmWr1KJKRG2q16pDZGRdyoWFcib9LKmpZ0hLSyU9LZW0tFSOJh3iyJ54Mg5vo+rZfUTJYXrKISrJGQgBE1IFiepuzXha/woIj9ZkoJSPs/OOoQOw3RizE0BEvgOuBXInhmuBFxzf/wB8ICJiSnpDiK8oXx1i74LYuxBjrF5NB1YTvH8N9XeuoHHSUsJS50AqcDDvpTlGqCSGfFZBsI4jnClXE6peTrlaV1gzxdbrgtRooW0GSpUwdiaG2sC+XD8nAhfOZ3D+HGNMloicBKoByblPEpFRwCjHj6dFZEsRYwoXyVu2jwgHn4wLCoztJLAFmGFjOOf56mvm0biKcfPlhrg8cudXKn+PxVCcuOpd6kCJbHw2xowHxhe3HBGJu1SrvDf5alzgu7FpXK7RuFxT2uKy8x5/P5B72bRIx758zxGRIKASViO0Ukopm9iZGFYC0SJSX0RCgJuBaRecMw24w/H9jcCf2r6glFL2su1RkqPN4CHgD6zuqhOMMZtE5CUgzhgzDfgcmCQi24FjWMnDk4r9OMpDfDUu8N3YNC7XaFyuKVVxlfiRz0oppdxL+xEqpZTKQxODUkqpPPwqMYhIPxHZIiLbReQf+RwPFZEpjuPLRSQq17GnHfu3iEhfZ8v0ZFwicrWIrBKRDY6vV+W6Zr6jzLWOrbqNcUWJSFquuj/OdU07R7zbReQ9KcJc38WI69ZcMa0VkRwRiXEcs+P16iEiq0UkS0RuvODYHSKyzbHdkWu/Ha9XvnGJSIyILBWRTSKyXkSG5To2UUR25Xq9YuyKy3EsO1fd03Ltr+/4nW93vAdcXuy7GK/XlRe8v9JF5DrHMTter8dEJN7xu5orIvVyHXPv+8sY4xcbVoP2DqABEAKsA5pdcM5o4GPH9zcDUxzfN3OcHwrUd5QT6EyZHo6rDXCZ4/sWwP5c18wHYr30ekUBGy9R7gqgE9bop5lAf7viuuCclsAOm1+vKKAV1nxfN+baXxXY6fhaxfF9FRtfr0vF1QiIdnx/GdZ498qOnyfmPtfO18tx7PQlyv0euNnx/cfAA3bGdcHv9BhQ1sbX68pc9T3A//4/uv395U93DOen3DDGZADnptzI7VrgS8f3PwC9HBn0WuA7Y8xZY8wuYLujPGfK9Fhcxpg1xpgDjv2bgDJiTTToDsV5vfIlIrWAisaYZcZ6V34FXOeluIY7rnWXQuMyxuw2xqwHci64ti8w2xhzzBhzHJgN9LPr9bpUXMaYrcaYbY7vDwBHgAgX63d7XJfi+B1fhfU7B+s9cJ2X4roRmGmMSXWx/uLENS9XfcuwxoKBB95f/pQY8ptyo/alzjHGZGHN4VCtgGudKdOTceU2BFhtjDmba98XjtvW54rwCKK4cdUXkTUiskBEuuc6P7GQMj0d1znDgMkX7PP06+XqtXa9XoUSkQ5Yn1R35Nr9quOxxTtF+EBS3LjCRCRORJade1yD9Ts+4fidF6VMd8R1zs1c/P6y8/W6G+sOoKBri/z+8qfE4LdEpDnWTLP35dp9qzGmJdDdsd1mY0gHgbrGmDbAY8C3IlLRxvoLJCIdgVRjzMZcu735evk0xyfLScCdxphzn5KfBpoA7bEeUTxlc1j1jDXVwy3AWBG53Ob6L8nxerXEGpN1jm2vl4iMAGKB/3iqDn9KDMWZcuNS1zpTpifjQkQigZ+A240x5z/NGWP2O76mAN9i3YraEpfjkdtRR/2rsD5lNnKcH5nrettfL4eLPs3Z9Hq5eq1dr9clORL6b8A/jTHLzu03xhw0lrPAF9j7euX+fe3Eah9qg/U7ruz4nbtcpjvichgK/GSMycwVry2vl4j0Bv4JXJPr6YH7319FbSzxtQ1rFPdOrMbjc403zS8450HyNlp+7/i+OXkbn3diNQYVWqaH46rsOP+GfMoMd3wfjPXM9X4b44oAAh3fN3C82aqa/Bu7BtgVl+PnAEc8Dex+vXKdO5GLG593YTUMVnF8b9vrVUBcIcBc4NF8zq3l+CrAWOB1G+OqAoQ6vg8HtuFoiAX+S97G59F2xZVr/zLgSrtfL6zkuANHhwFPvr+cDrwkbMAAYKvjxfunY99LWNkVIMzxxtrueMFy//H4p+O6LeRquc+vTLviAp4FzgBrc23VgXLAKmA9VqP0uzj+UNsU1xBHvWuB1cDgXGXGAhsdZX6AY3S9jb/HnsCyC8qz6/Vqj/Uc9wzWp9tNua69yxHvdqxHNna+XvnGBYwAMi94f8U4jv0JbHDE9jVQ3sa4ujjqXuf4eneuMhs4fufbHe+BUJt/j1FYHzwCLijTjtdrDnA41+9qmqfeXzolhlJKqTz8qY1BKaWUG2hiUEoplYcmBqWUUnloYlBKKZWHJgallFJ5aGJQSimVhyYGpZRSeWhiUMoDxFr/oYnj+2oisrGwa5TyFZoYlPKMhlijWMGa23+DF2NRyiWaGJRyM8fKWvvN/2YqbYU1HYdSJYImBqXcrzV5E0E7NDGoEkQTg1LuF4M10R8iEo21Epc+SlIlhiYGpdyvNRAgIuuA54F44A7vhqSU83R2VaXcTES2AW2NtSiQUiWO3jEo5UYiUgEwmhRUSaZ3DEoppfLQOwallFJ5aGJQSimVhyYGpZRSeWhiUEoplYcmBqWUUnloYlBKKZWHJgallFJ5/D+wKrpWP9jmpgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "x = np.linspace(0., 0.2, 50)\n", "fixed = []\n", "prof = []\n", "for x_i in x:\n", " vals = np.copy(best_fit.x)\n", " vals[mu_idx] = x_i\n", " fixed.append(2 * (model.nll(vals, asimov) - best_fit.fun))\n", " prof.append(2 * (model.profile({mu: x_i}, best_fit, asimov).fun - best_fit.fun))\n", "\n", "plt.plot(x, fixed, label=\"Fixed\")\n", "plt.plot(x, prof, label=\"Profiled\")\n", "plt.plot([x[0], x[-1]], [1.]*2, \"--\", color=\"black\")\n", "plt.plot([up, up], [0., 1.], \"--\", color=\"orange\")\n", "plt.plot([down, down], [0., 1.], \"--\", color=\"orange\")\n", "plt.plot([up_stat, up_stat], [0., 1.], \"--\", color=\"blue\")\n", "plt.plot([down_stat, down_stat], [0., 1.], \"--\", color=\"blue\")\n", "plt.gca().set_ylim([0.,1.5])\n", "plt.gca().set_ylabel(\"$-2 (\\log L - \\log L_0)$\")\n", "plt.gca().set_xlabel(\"$\\mu$\")\n", "plt.legend()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As you can see, the 1-sigma Minos bounds obtained above are correct! Furthermore, the lower bound on $\\mu$ is properly taken into account." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Pulls and impacts\n", "\n", "Getting nuisance parameter pulls is easy:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Lumi nuisance: pull = -0.00, +0.98 -0.97\n", "Shape nuisance 0: pull = -0.00, +0.91 -0.85\n", "Shape nuisance 1: pull = -0.00, +0.65 -0.49\n" ] } ], "source": [ "pull_u_lum,pull_d_lum,_,_ = model.minos_bounds(lumi_nuis, best_fit, asimov)\n", "lumi_idx = lumi_nuis.sub_idxs[0]\n", "print(f\"Lumi nuisance: pull = {best_fit.x[lumi_idx]:4> .2f}, +{pull_u_lum:.2f} {pull_d_lum:.2f}\")\n", "\n", "pulls_u = []\n", "pulls_d = []\n", "for i in range(2):\n", " pull_u,pull_d,_,_ = model.minos_bounds(shape_nuis, best_fit, asimov, sub_idx=i)\n", " nuis_idx = shape_nuis.sub_idxs[i]\n", " print(f\"Shape nuisance {i}: pull = {best_fit.x[nuis_idx]:4> .2f}, +{pull_u:.2f} {pull_d:.2f}\")\n", " pulls_u.append(pull_u)\n", " pulls_d.append(pull_d)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And so are impacts (we use the up/down intervals from the pulls to calculate those):" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "scrolled": true }, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/home/swertz/Documents/Recherche/Stats/Tests/smoofit_venv/lib/python3.8/site-packages/scipy/optimize/_hessian_update_strategy.py:182: UserWarning: delta_grad == 0.0. Check if the approximated function is linear. If the function is linear better results can be obtained by defining the Hessian as zero instead of using quasi-Newton approximations.\n", " warn('delta_grad == 0.0. Check if the approximated '\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Lumi nuisance: impacts = -0.003, 0.003\n", "Shape nuisance 0: impacts = -0.010, 0.012\n", "Shape nuisance 1: impacts = -0.014, 0.012\n" ] } ], "source": [ "impact_up = model.profile({lumi_nuis: pull_u_lum}, best_fit, asimov).x[mu_idx] - best_mu\n", "impact_down = model.profile({lumi_nuis: pull_d_lum}, best_fit, asimov).x[mu_idx] - best_mu\n", "print(f\"Lumi nuisance: impacts = {impact_up:.3f}, {impact_down:.3f}\")\n", "\n", "for i in range(2):\n", " impact_up = model.profile({shape_nuis: (pulls_u[i], i)}, best_fit, asimov).x[mu_idx] - best_mu\n", " impact_down = model.profile({shape_nuis: (pulls_d[i], i)}, best_fit, asimov).x[mu_idx] - best_mu\n", " print(f\"Shape nuisance {i}: impacts = {impact_up:.3f}, {impact_down:.3f}\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": true }, "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.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }