adversarial_histogram_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adversarial training: binned data differences\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here is some pseudo-data of a variable, pt:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt = np.array([50, 50, 100, 101, 102, 103, 104, 200, 201, 202, 203, 300, 301, 302], dtype=\"float\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Here is the bins that the measurements occur in: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = np.array([0,100,200,300,400], dtype=\"float\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  0., 100., 200., 300., 400.])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### First we digitize the values of pt in the bins to find the categories\n",
    "\n",
    "We have to offset the bin boundaries by one to get the index in the array. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "categories = np.digitize(pt, bins=bins) - 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 0, 1, 1, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "categories"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next make a histogram of the pseudodata: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n, xbins = np.histogram(pt, bins=bins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now get the \"observed\" values (called \"hepdata\"): "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "hepdata = np.array([0, 3, 2, 1], dtype=\"float\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "hepdata = hepdata / np.sum(hepdata)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = n / np.sum(n)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the two distributions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f66f2a2ed30>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvC0lEQVR4nO3deXyU1b3H8c+PsCvglaBXRFaDiiIIkVpRL24I2qIF677gxq0LLq1WUNyo+14sWpHiUhcUXIpbERUVXAkIKiCCXJAgFUQFBSQk/O4fZ8YMIQmBzJOZyXzfr9e8mGeeZ2Z+TybkN+c853eOuTsiIpK96qQ6ABERSS0lAhGRLKdEICKS5ZQIRESynBKBiEiWq5vqALZWbm6ut23bNtVhiIhklOnTp3/r7i3K25dxiaBt27YUFBSkOgwRkYxiZosr2qeuIRGRLKdEICKS5ZQIRESyXMZdIxARSYYNGzZQWFjIzz//nOpQkqphw4a0atWKevXqVfk5SgQikpUKCwtp0qQJbdu2xcxSHU5SuDsrV66ksLCQdu3aVfl56hoSkaz0888/07x581qTBADMjObNm291K0eJQESyVm1KAnHbck5KBCIiWU6JQEQkRQoLCzn22GPJy8ujQ4cOXHLJJRQVFfHII49w0UUX1VgcSgQiIing7vTv35/jjjuO+fPn88UXX/DTTz9x9dVX13gskSYCM+tjZvPMbIGZDSln/z1mNjN2+8LMfogyHhGRdPHmm2/SsGFDzjrrLABycnK45557GDNmDGvXrmXJkiX06tWLvLw8brjhBgDWrFnDMcccQ5cuXdhnn314+umnkxJLZMNHzSwHGAkcCRQC08xsgrvPiR/j7pclHD8Y2C+qeEREKnLppTBzZnJfs2tXuPfeivfPnj2b7t27b/JY06ZNad26NcXFxXz00Ud89tlnNG7cmP33359jjjmGxYsX07JlS15++WUAVq1alZRYo2wR9AAWuPtCdy8CxgLHVnL8ycBTEcYjIpIxjjzySJo3b06jRo3o378/U6dOpXPnzkyaNIkrr7ySKVOm0KxZs6S8V5QFZbsCSxK2C4FflXegmbUB2gFvVrB/EDAIoHXr1smNUkSyXmXf3KPSqVMnxo8fv8ljq1ev5quvvqJu3bqbDQM1Mzp27MiMGTN45ZVXGDZsGIcffjjXXntttWNJl4vFJwHj3b2kvJ3uPsrd8909v0WLcqfTFhHJKIcffjhr167lscceA6CkpIQ//elPDBw4kMaNGzNp0iS+++471q1bxwsvvEDPnj35+uuvady4MaeddhpXXHEFM2bMSEosUSaCpcBuCdutYo+V5yTULSQiWcTMeP755xk3bhx5eXl07NiRhg0bcvPNNwPQo0cPBgwYwL777suAAQPIz8/n008/pUePHnTt2pUbbriBYcOGJScWd0/KC232wmZ1gS+AwwkJYBpwirvPLnPcnsC/gXZehWDy8/NdC9OISHXNnTuXvfbaK9VhRKK8czOz6e6eX97xkbUI3L0YuAiYCMwFnnH32WY23Mz6JRx6EjC2KklARESSL9LZR939FeCVMo9dW2b7+ihjEBGRyqXLxWIREUkRJQIRkSynRCAikuWUCEREspwSgYhIitx0003svffe7LvvvnTt2pUPP/yQc889lzlz5mz5yUmkNYtFRFLg/fff56WXXmLGjBk0aNCAb7/9lqKiIkaPHl3jsahFICKSAsuWLSM3N5cGDRoAkJubS8uWLenVqxfxotntt9+eK664gr333psjjjiCjz76iF69etG+fXsmTJiQtFiUCEREAHr12vx2//1h39q15e9/5JGw/9tvN9+3Bb1792bJkiV07NiRCy64gLfffnuzY9asWcNhhx3G7NmzadKkCcOGDWPSpEk8//zzSZlsLk5dQyIiKbD99tszffp0pkyZwuTJkznxxBO59dZbNzmmfv369OnTB4DOnTvToEED6tWrR+fOnVm0aFHSYlEiEBEBeOutivc1blz5/tzcyvdXICcnh169etGrVy86d+7Mo48+usn+evXq/TIddZ06dX7pRqpTpw7FxcVb/X4VUdeQiEgKzJs3j/nz5/+yPXPmTNq0aZOSWNQiEBFJgZ9++onBgwfzww8/ULduXXbffXdGjRrF8ccfX+OxRDYNdVQ0DbWIJIOmoS6lriERkSynRCAikuWUCEQka2Va13hVbMs5KRGISFZq2LAhK1eurFXJwN1ZuXIlDRs23KrnadSQiGSlVq1aUVhYyIoVK1IdSlI1bNiQVq1abdVzlAhEJCvVq1ePdu3apTqMtKCuIRGRLBdpIjCzPmY2z8wWmNmQCo45wczmmNlsM3syynhERGRzkXUNmVkOMBI4EigEppnZBHefk3BMHjAU6Onu35vZTlHFIyIi5YuyRdADWODuC929CBgLHFvmmPOAke7+PYC7L48wHhERKUeUiWBXYEnCdmHssUQdgY5m9q6ZfWBmfcp7ITMbZGYFZlZQ267wi4ikWqovFtcF8oBewMnAQ2a2Q9mD3H2Uu+e7e36LFi1qNkIRkVouykSwFNgtYbtV7LFEhcAEd9/g7v8HfEFIDCIiUkOiTATTgDwza2dm9YGTgLKLbL5AaA1gZrmErqKFEcYkIiJlRJYI3L0YuAiYCMwFnnH32WY23Mz6xQ6bCKw0sznAZOAKd18ZVUwiIrI5rUcgIpIFtB6BiIhUSIlARCTLKRGIiGQ5JQIRkSynRCAikuWUCEREspwSgYhIllMiEBHJckoEIiJZTolARCTLKRGIiGQ5JQIRkSynRCAikuWUCEREspwSgYhIllMiEBHJckoEIiJZTolARCTLKRGIiGQ5JQIRkSwXaSIwsz5mNs/MFpjZkHL2DzSzFWY2M3Y7N8p4RERkc3WjemEzywFGAkcChcA0M5vg7nPKHPq0u18UVRwiIlK5KFsEPYAF7r7Q3YuAscCxEb6fZJhVq+Cgg+CMM2BO2a8HIlJjokwEuwJLErYLY4+VNcDMPjGz8Wa2W3kvZGaDzKzAzApWrFgRRaySAoMHwwcfwLPPwt57Q//+MG1aqqMSyT6pvlj8ItDW3fcFJgGPlneQu49y93x3z2/RokWNBijReOYZ+Oc/4eqr4auv4NprYfJk6NEDeveGt94C91RHKZIdokwES4HEb/itYo/9wt1Xuvv62OZooHuE8UiaWLoU/vCH8Ed/2DBo3hxuuCEkhNtvh08/hUMPhZ494aWXlBBEohZlIpgG5JlZOzOrD5wETEg8wMx2SdjsB8yNMB5JAxs3wsCBsH49PP441KtXuq9JE7jiCvi//4MHHoBly+C3v4WuXeGpp6C4OFVRi9RukSUCdy8GLgImEv7AP+Pus81suJn1ix12sZnNNrNZwMXAwKjikfRw333w+utw992Ql1f+MQ0bhhbDF1+E7qPiYjjlFNhzT3jooZBERCR5zDOs3Z2fn+8FBQWpDkO2wezZ0L07HHkkTJgAZlV73saN4fibboKCAmjZEi6/HAYNgu22izZmkdrCzKa7e355+1J9sViyRFERnHYaNG0Ko0dXPQkA1KkDxx0HH30EkybBHnvAH/8IbdrAX/4C338fWdgiWUGJQGrEtdfCzJkhCey887a9hhkccQS8+Sa89x4ceGB43dat4cor4T//SWrIIllDiUAi9/bbYTTQeedBv35bPr4qfv3r0F00a1a4oHznndC2LVx4ISxalJz3EMkWSgQSqVWrQuVwhw7hAnGy7bsvPPkkzJsX3uehh2D33eHMM1WtLFJVSgQSqcGDQ93AP/8J228f3fvsvjuMGhWGnl5yCYwfX1qtrLEFIpVTIpDIjBtXWj18wAE185677gp33QWLF5dWK++/v6qVRSqjRCCRWLoU/vd/wx/hYcNq/v1zc0O18uLF4frEJ5+oWlmkIkoEknQbN8JZZ5VfPVzTmjYN1cqLFsH9929arTx2LJSUpC42kXShRCBJd999Ybz/3XdDx46pjiZo2BDOPz9UKz/2GGzYACefHKqVR49WtbJkNyUCSarZs8OY/t/8JlT+ppt69eD00+Gzz+C552CHHcKw1g4d4N57Yc2aVEcoUvOUCCRpqlM9XNPq1IHf/S5UK7/2Wmi5XHZZqFa+8UZVK0t2USKQpElG9XBNMwtzHyVWK19zTUgIqlaWbKFEIEnxzjthdM655yaverimJVYr/+Y3qlaW7KFEINUWrx5u3x7uuSfV0VRfZdXKc7VihtRCVUoEZtbfzOab2SozW21mP5rZ6qiDk8wweDAUFoaholFWD9e0eLXywoVw8cWl1coDBqhaWWqXqrYIbgf6uXszd2/q7k3cvWmUgUlmSEX1cE1r1SoMhV28OBTHvflmKJQ76qgwoZ6K0yTTVTURfOPuahTLJlJdPVzTcnNh+PCQEG67LVxL6NULDjoIXn5ZCUEyV1UTQYGZPW1mJ8e6ifqbWf9II5O0lk7VwzWtaVP485/DBHf33x8S4m9+o2plyVxVTQRNgbVAb+C3sdtvogpK0t/f/haqh++6K32qh2tao0ahWnn+fFUrS2bTmsWy1WbPhvx8OPxwePHF9C4cq0kbN8K//hXWVp4+PcyEevnloXJZaytLqlV7zWIza2Vmz5vZ8tjtWTNrVYXn9TGzeWa2wMyGVHLcADNzMys3SEkf8erhJk3gH/9QEkgUr1aeNi1UK+flqVpZMkNVu4YeBiYALWO3F2OPVcjMcoCRQF+gE3CymXUq57gmwCXAh1UPW1IlE6uHa1q8WnnyZHj33VCoFq9WHjIEvvkm1RGKbKqqiaCFuz/s7sWx2yNAiy08pwewwN0XunsRMBY4tpzj/gLcBvxc1aAlNWpD9XBNO/DA0H02cyYccwzccUeoVr7oojD6SCQdVDURrDSz08wsJ3Y7DVi5hefsCixJ2C6MPfYLM+sG7ObuL1f2QmY2yMwKzKxgxYoVVQxZkqm2VQ/XtC5d4KmnQrXy6aeHQrXdd4eBA1WtLKlX1URwNnAC8B9gGXA8cFZ13tjM6gB3A3/a0rHuPsrd8909v0WLLTVEJAoXXwxLlkS/9nBtl1itPHhwKMiLVytPn57q6CRbVSkRuPtid+/n7i3cfSd3P87dv9rC05YCuyVst4o9FtcE2Ad4y8wWAQcAE3TBOP2MGxeGRw4bFvq7pfrKq1bOz1e1sqRGpcNHzezP7n67md0HbHagu19cyXPrAl8AhxMSwDTgFHefXcHxbwGXu3ulY0M1fLRmLV0KnTuHb7LvvptdhWM1afVq+PvfQ3L45ptwbeGqq+DoozUyS5KjOsNH472XBcD0cm4Vcvdi4CJgYux1nnH32WY23Mx0qTEDZHP1cE1LrFYeObK0Wnm//eDpp1WtLNGqUkGZmf3e3cdt6bGaoBZBzRkxAi65BB54AP7wh1RHk102bAgXl2+5BT7/PLTIhgwJF5rr1091dJKJql1QBgyt4mNSS8yZE1boOuaYMLGc1Kx69cIordmz4dlnoVmzMGy3Qwf461+1trIkV6WJwMz6xq4P7GpmIxJujwDFNRKh1LiiIjj11DA6KN3XHq7t6tSB/v1DtfLEiaFlcOmlpdXKP/yQ6gilNthSi+BrwvWBn9n02sAE4KhoQ5NUue660urh//7vVEcjEJJx796bVyu3bq1qZam+ql4jaAqscfeS2HYO0MDd10Yc32Z0jSBa77wT5tg/55ywRKOkr1mz4NZb4ZlnwnWDc86BK64IrQWRspJxjeA1oFHCdiPg9eoGJulF1cOZJV6t/PnnYSJAVSvLtqpqImjo7j/FN2L3G0cTkqRKyqqHX3opzL2wcWMNvmntkZcXWm8LF4Y5jOLVyscfr2plqZqqJoI1sXmBADCz7sC6aEKSVBg/PlQPX311DVUPxzu1160LczfvuSfssAMcemjo33j//RoIonZp1Sq05BYtCp/j66+HauU+fUKXn6qVpSJVTQSXAuPMbIqZTQWeJhSLSS3w9delaw9fc03Eb/bNN2Ew/J57wn/+Aw0awIwZMGZM6N9YsyYUMLz3XmlwvXuHv2zPPx+aLPqLVqkWLeAvf4GvvgrXEGbOhP/5Hzj4YHjlFf34pBzuXqUbUI8wN9A+QL2qPi/Zt+7du7skT0mJe+/e7o0bu8+bF+EbFRe7jxzp3qyZe/367tdc4752bfnHrl/v/tNP4f7Mme5du7rn5LiHv2HuO+3k/sYbYf/337svWxZh4Jlv7drwo2/TJvz4unRxHzs2fCSSPYACr+DvalVbBAB7EBaY6UZYZOaMZCclqXkjR4bVtCJde3jt2tDfdOGFoa/ik09g+PCw6G956tcvXduxSxf4+GP48cfQXXTffdC3b+nQmHHjYJddYLfdQhfTzTeHE/pZy1vENWoEF1wQ1lZ+9NEwZchJJ8Fee4VV5oqKUh2hpFpVh49eB/QiJIJXCKuOTXX34yONrhwaPpo8c+ZA9+4Rrj28YUPpBEWXXQY9eoS/QMl8o3nzQn9HQUGoupo/Pzy+ciXsuCM89xwsWBD6vbp1CyW6WW7jRnjhhZAzp08P1xYuvzxULmtt5dqrsuGjVe0W+pRwPWFWbHtnYFJVnpvsm7qGkmP9+tDjkpsbQc/Kxo3uTzzh3qqV+yefJPnFt+D7792nTCndPuec0i4lcN9jD/fzzivdn8X9Ixs3uk+c6P4//xN+NLm57jfeGH6EUvuQhK6hde6+ESiOFZctZ9O1BiTDRFY9PG9eWLD31FNDl01Nz0+xww5w0EGl26NHw4oV8O9/hyuoe+4ZWgtxhxwC++wTplkdORI++ihrupXi1cpvvQVTp8KvfhXWRmjdGoYOVbVyVqkoQyTegPuBHYA/APOBj4GHq/LcZN/UIqi+d95xNwtflpNq+PBwIbhZM/f778+Mb9s33eR+9NHuLVqUthr69y/d/9hj7tOnuxcVpS7GGvTxx+4nnuhep457w4buF13kvmhRqqOSZKCSFsEWrxGYmQGt3H1JbLst0NTdP4ksO1VC1wiqZ/Vq2HdfqFs3tAiSWjg2bFhYcuvOO2HnnZP4wjXAPQxNLSgI1xZ69YJvvw1jMSEMc+3SJVxrOPXUWr9U2/z5cPvt4eKyezjlIUNCg0oyU2XXCKp6sfhTd++c9Mi2gRJB9QwcGCqHp05Nwt+yJUvCggXnnRdG8rjXrqlK3UO5bkFB6cXo6dPh3nvDxD6ffx6usObnhwSRnx/KfOtszWC89FZYGEaUPfhg6DHr3z90G3XvnurIZGsl42Lxo8D+VTk26pu6hrbduHGh5+Oaa6r5QkVF7rff7r7ddu6NGrmPGpWU+DJCSUm40u7u/uGH7gceGH4G8W6lpk3d33sv7P/Pf9y//DJclc1wy5e7DxsWev3A/aij3N9+u1acWtagOl1DsUzyOZAHLALWABZyiO+bnFxVdWoRbJuvvw5rD3foUM21h997L5Qhf/ZZWEtxxAho1y6psWac4uIwy9u0aaHlcP31sNNOoW/lyitDV1N+fmnLoW/f0NWUgVavDivW3X03LF8OPXuGtZX79q1djcHaaJu7hsystbt/ZWblTmzr7ouTFGOVKRFsvY0bw3/UKVNCbdYee1TjxUaPDqNvRoyAfv30v78yCxbAG2+Udi19+mnoNvrxx5AIxowJ3WvxJJFB11XWrQvh3357mMqiS5eQEAYMgJycVEcn5dnmriFgRsL9Zys7tqZu6hraeiNGhOb8/fdvw5NLStwfesj9kUdKt+PTP8jWWbvWfdas0u0zzgjDt+LdSrvt5n7uuaX7162r+Ri3UlFR+NXYc89wCnl57v/4R2nvmaQPqlFHkPh1r/02ZKA+ZjbPzBaY2ZBy9v/BzD41s5lmNtXMOm3te0jl5syBP/8Zjj56GxagnzUrjMk/77xQoQvhG63KT7dNo0ZhyFbco4+GvpZ33gl9LQcdtOnX6a5dQ1/eiSeGkVhvvRWOTyP16sGZZ5aurdykSbiOrrWVM0xFGcI3bxHMqOzYcp6bA3xJSCD1gVlApzLHNE243w/495ZeVy2Cqlu/3n2//bahenj1avfLLgsTveXmhq98uipYszZudL/tNvfjj3dv27a01XDWWaX777vPferUtGqhqVo5fVFJi2BLf8xLgNXAj4TF6lcnbK/ewnN/DUxM2B4KDK3k+JOBVyt7TVci2CpDhoRP+IUXtvKJb7wRuiz+93/dV66MJDbZSsuXu7/6qvu0aWF7wYLS5FCnjnvnziFJvP9+auNMMHWq+zHH+C+DqYYMCQOpJDW2ORFU5wYcD4xO2D4d+Fs5x10YazksAfIqeK1BQAFQ0Lp16wh/VLXHVlcPL1hQeh3A3f2LLyKJS5Jo6VL3CRPcr73WvW/fUB397LNh3/vvu3fr5j5oUBje+/HHKauOjlcrm6laOZXSOhEk7D8FeHRLr6sWwZatWhXmnm/fPvTyVGrdOvcbbnBv0MB9xx3DkyUzbdxYOq3H1KnuRxzhvsMOpS2HBg3CX2X38Jf4s89qdBqQefPCF5N69dzr1nUfONB97twae/usV1kiiLIEcimbTkzXKvZYRcYCx0UYT9aIrz38+OPh4l2FJk0KFy+vuw6OOy4Mb2zatKbClGQzK73Y3LNn+Hy/+y4MYx07NixonJcX9j/4YJhsr1mzsHTZH/8ITz4ZaiIi0rFjGH385ZdhaYqnn4ZOneD3vw+L1EnqVKmgbJte2Kwu8AVwOCEBTANOcffZCcfkufv82P3fAtd5ReNcY1RHULlnnw2Lll9zTVj7pUJLl4ZCsDZtwqybvXvXWIySBhYtCoUl8RqHjz8OCwJ9/31IKHfdFSrG4jUObdsmvWZkxYowsuhvf4NVq+Coo0ItwiGHJPVtJKbaU0xs6w04mpAMvgSujj02HOgXu/9XYDYwE5gM7L2l11TXUMWWLg29O/vvX0F38IYN7i+9VLr9+usZMVZdasCGDe4LF5Zun3JK6MOJdys1b77pOg7ffZe0kWQ//OB+yy1hBVJw79nT/eWXNVAt2ajuFBPpRC2C8m2xevjDD+H888PODz4Ik8+LVGb9+jCVSLzVsPPOcOONITXE15qItxjy88MKdPHZWreBqpWjlbIWQRQ3tQjKV2H18HffhWGgZu4tW4aZ5/RVS6pjw4ZQw3DGGe6dOpVWR192Wdi/fr37zTe7v/baNg0/jlcr77GHq1o5iUjFqKGobkoEm5s9OwzLO/roMn/jS0pC7X9OTvhPusUhRCLb4Mcfw3jlOXPC9qxZpV1KEIavnXBCmK11KxQXu48fH0bBQlj59K9/dV+zJoJzyAKVJQJ1DWW4oiI44IAwSujTT2PLTs6fD+3bh/b0iy+GtQe7dEl1qJJNvv8+DAWKdysVFMA//gGHHQavvgqXXVa6hkN+fphOo4KpS9zhtdfg5pvDbBy5ueHpF1wQViaVqlHXUC02dGj4tvT88x6mGrjyyjBI+8EHUx2ayKbizdW333Y/9tjQVZlYHR0vYpwzx/2jj9x//nmzl5gyJbR849XKQ4e6f/NNzZ1CJkNdQ7XTJtXD//pXqCKLz0ezfHmqwxPZsnh19HXXlRa3nX9++D2uV2/T6uiEfs+PPw69TapWrrrKEoG6hjJUfO3hnByY2/sS6v99RCgQeuCBMIulSKZaujSMbEvsVmrWLNQ+AFx9dfgPkJ/Potx8bnp2Tx59PAd3OO20sBaQ1lbeXLXXLE4nSgTBOacXMfaJEl5/txG/XvN6WIn+kkuqsfSYSJpyD9VnO+0Utk8+GV56CX76KWxvtx0//e50hjV/gFGjIHfdEg4YsCtDrqpDt26pCzvdKBHUMm8Pf5vc6y5g5YG/5ZB3b011OCI1r6QEvviitMXQti1cdhkrlhXTrHVT1hXXp4DufNcun30G5rPXuT2hZctUR51SSgS1xfLlrL3wChqPf4yl9duy87iR1O13dKqjEkkf69fDE0+w/t1pfDexgOZLZ1GfDYxuPZxd/34NfQ74AbvrztIRSy1bZs1yq0oEtcGLL+JnnEHxqjXcnfNnfjftKjp2bZzqqETS2trv1/PiLZ9x7+O5fLCsDWfmvceYLw+hzsaScMB//3dICNdfD927h26oWpoYKksEUc4+KsmwcWP4t0MHvtrlV3T2T2g64kYlAZEqaPxfDTjx9u68vagNDz8MH9Q5kO02/sgJrd7j/VNGUHJEb1i4MCzBCvDEE6HuZsAAuOWW0hlcazm1CNLVqlUwbBisXAlPPsmcOeELy2GHhetktfRLi0ikSkrg+edDcdrHH0OrVnDFFXDuudC4MTB5Mjz0EEybFqbvjlu2LLQeCgrCRepu3TJuyna1CDKJOzz1VBj/NnIk7LgjRetKOO002H77UJypJCCybXJywjTt06eHAud27cJguzZtQnL4Yb9Dw7oM8+eHlsDrr8M998RK9gnTcx96aChp3msvOP10uP/+lJ5TMqhFkE6++grOPhveeCN8/f/73yE/n6uuCq3U558P68eISPJMnRr+f73ySviSf+GFcOmlpaNVN7FiRcgi06aVjljaaafQvAAYNCgs7hO/GL3vvtCgQU2eToU0xUSmWL7cvV0795Ejf6mynDIlVN+ffXaKYxOp5WbM2LRaefBg98WLq/DExMkcTzzRPTe3dOqMevXcL7ywdP+cOSlbOxpVFqexV1+Fxx4LF6nq1AnfJurWBULxZJcu4eGZM7ew7KSIJMW8eWFNhMceC9tbXa3sHlr38RbDHnvAwIHh2kLTpqGF0LVr6YR7hx4aLlBHTC2CdLRkifuAAeFbwx57uBcWbnbIwIGhNfDeeymITyTLLV7sfvHF7o0ahVbC8ce7T59ejRdcs8b9qafc//Qn90MOcd9++/D/f8SIsH/JkjBd/BNPhAn4SkqSch5xaNK5NLJhg/udd7pvt11of950U7mzLI4fHz6dYcNSEKOI/OKbb9yvuirMdgruffqECR+rrbg4dBXFp099/fXwNyHerdSsmfthh4X1HeLHV0NliUBdQzXt55+hc+fQXLzvvjBsoYyvvw6HtG8P772n6YNE0sGqVWGA0D33hGvGBx0UltLs0yeJI/mKi2HOnNClFL8g/fTT4Y9BYWEY77qNVFmcat9+C7fdFqoXt9su/Bbl5pb72+Me1h5+550K1h4WkZRauzYM477jjrAgVNeuISH07x/x2so//RTGkG+jlNURmFkfM5tnZgvMbEg5+/9oZnPM7BMze8PM2kQZT43buDH8xuy5Z/ga8fbb4fEWLSr8CjFyJEycCHfeqSQgko4aN4bBg0O92cMPh8RwwgnQqVPYLiqK6I2rkQS2JLJEYGY5wEigL9AJONnMOpU57GMg3933BcYDt0cVT4375BM4+OBQsrjXXmHYz9GVTxA3d26ocuzbF84/v2bCFJFtU79+GAw0Zw6MGxca+2efDbvvDiNGhASRKaJsEfQAFrj7QncvAsYCxyYe4O6T3T3+4/oA2PYOsHRz2WVhmtyHHw4tgX32qfTwoiJ+qR4eM0bVwyKZomy1ctu2oVq5bdtYtfIPKQ6wCqJMBLsCSxK2C2OPVeQc4NXydpjZIDMrMLOCFStWJDHEJHKH554Lc5JA6BL6/PPwlaHOln/M118f1vp+6KHSanYRyRxm4cLxO++EW35+WEytTZtwDWH58lRHWLG0mGvIzE4D8oE7ytvv7qPcPd/d81u0aFGzwVXFwoVwzDFhxsK//jU81rYtNG9epadPnRquJZ99tqaQEKkNDj44TFkxYwYcdRTcemtICBdfHGrN0k2UiWApsFvCdqvYY5swsyOAq4F+7r4+wniSb/16+MtfYO+9YcqUcEH4xhu36iVWrw7zVrVtC/feG0mUIpIi++0HzzwTrv+dfHJYUrxDh/Clb968VEdXKspEMA3IM7N2ZlYfOAmYkHiAme0HPEhIAmnccKrANdfAtddCv36hG+jSS3+ZHqKqLrkkfEP45z81hYRIbbXHHuHa35dfwgUXwNixYQzJCSeUzleXSpElAncvBi4CJgJzgWfcfbaZDTezfrHD7gC2B8aZ2Uwzm1DBy6WPZcvCFLUAl18O//53KPjYtbLLH+V77jl45JHQf3jggckNU0TST+vWofd40SIYOjQMFe/WLYwUnDIldXGpoKyqSkpCWeGwYeGTmzy5Wi+n6mERqZFq5RgtTFNd06ZBjx7hSs8BB8CoUdV6OffQR7huHTz+uJKASLZq1iy0DBYtCrUHixeHcqNu3UJtQklJzcShRLAlL70Ev/pV6BJ6+unQFZSXV62XVPWwiCRKrFYeM6YGq5VjlAjK4w5LYwOcDj8crrsuXAw+4YRqt9dUPSwiFalfH846K1QrP/NMSBDxauX77ouuWlmJoKw5c8JCEQcdFH7qjRqFRJCEharj1cPbbae1h0WkYjk58PvfhzqEV18trUEYOTKa99u6sY612Zo1oSbgrrvCOM5bb4WGDZP6FjfcED7Y556DXXZJ6kuLSC0Ur1bu0yeMKurcOZr3USKAMJfswQeHKzVnnhnWqSt35eptN3VqyC1nnw2/+11SX1pEssDBB0f32tmdCH7+OXzr33VXOOIIOOMMOOSQpL+NqodFJJ1l5zWCoqLw9bx9+zAaqE4dGD06kiQAqh4WkfSWfYng7bfDBCBDh4bagIgL6uLVw0OHqnpYRNJT9iSCkpIwJXSvXmE00IQJ8MIL0LJlZG+5bBkMGgTdu4eBRyIi6Sh7EkFOTrhddRXMng2//W2kbxevHl67VtXDIpLesuti8ejRNTZ4//77QxHyyJFhyWIRkXSVPS0CqLEkMHdumJhU1cMikgmyKxHUAFUPi0imya6uoRqg6mERyTRqESTRu++G8oSzzlL1sIhkDiWCJEmsHo6vXy8ikgnUNZQkl14apiqaMkXVwyKSWdQiSILnnguLR6h6WEQykRJBNal6WEQyXaSJwMz6mNk8M1tgZkPK2X+Imc0ws2IzOz7KWKKg6mERqQ0iSwRmlgOMBPoCnYCTzaxTmcO+AgYCT0YVR5Ti1cN33KHqYRHJXFFeLO4BLHD3hQBmNhY4FpgTP8DdF8X2bYwwjkh8/nmoHu7TBy64INXRiIhsuyi7hnYFliRsF8Ye22pmNsjMCsysYMWKFUkJrjqKiuDUU0P18Jgxqh4WkcyWEReL3X2Uu+e7e36LFi1SHc4v1cMPPaTqYRHJfFEmgqXAbgnbrWKPZTRVD4tIbRNlIpgG5JlZOzOrD5wETIjw/SIXrx5u00bVwyJSe0SWCNy9GLgImAjMBZ5x99lmNtzM+gGY2f5mVgj8HnjQzGZHFU8yxKuHtfawiNQmkU4x4e6vAK+UeezahPvTCF1GaS9ePXz11dCzZ6qjERFJnoy4WJxqqh4WkdpMiWALVD0sIrWdZh/dgnj18N/+puphEamd1CKohKqHRSQbKBFUIHHtYVUPi0htpq6hCgwfDtOnw7PPqnpYRGo3tQjK8e67cMstoXq4f/9URyMiEi0lgjJUPSwi2UZdQ2XEq4ffeUfVwyKSHdQiSBCvHh4yRNXDIpI9lAhi4tXD3bqpelhEsosSAaF6+JxzQvXwE09A/fqpjkhEpOboGgHwwAPw6quqHhaR7JT1LYLPP4c//UnVwyKSvbI6Eah6WEQky7uGVD0sIpLFLYL33gvVwwMHqnpYRLJbViaCH39U9bCISFxWdg1deiksWhSqh5s2TXU0IiKplXUtgueeCxeGVT0sIhJEmgjMrI+ZzTOzBWY2pJz9Dczs6dj+D82sbZTxqHpYRGRzkSUCM8sBRgJ9gU7AyWbWqcxh5wDfu/vuwD3AbVHFE68eXrMmrD2s6mERkSDKFkEPYIG7L3T3ImAscGyZY44FHo3dHw8cbhbNaP549fAdd8Bee0XxDiIimSnKi8W7AksStguBX1V0jLsXm9kqoDnwbeJBZjYIGBTb/MnM5m1jTLmDB/Pt4MHb+Oz0kkuZn1MGqy3nUlvOA3Qu6ao659Kmoh0ZMWrI3UcBo6r7OmZW4O75SQgp5XQu6ae2nAfoXNJVVOcSZdfQUmC3hO1WscfKPcbM6gLNgJURxiQiImVEmQimAXlm1s7M6gMnARPKHDMBODN2/3jgTXf3CGMSEZEyIusaivX5XwRMBHKAMe4+28yGAwXuPgH4B/BPM1sAfEdIFlGqdvdSGtG5pJ/ach6gc0lXkZyL6Qu4iEh2y7rKYhER2ZQSgYhIlsuaRLCl6S7SnZktMrNPzWymmRXEHtvRzCaZ2fzYv/+V6jjLMrMxZrbczD5LeKzcuC0YEfuMPjGzbqmLfHMVnMv1ZrY09rnMNLOjE/YNjZ3LPDM7KjVRl8/MdjOzyWY2x8xmm9klsccz6rOp5Dwy7nMxs4Zm9pGZzYqdyw2xx9vFpuBZEJuSp37s8eRN0ePutf5GuFj9JdAeqA/MAjqlOq6tPIdFQG6Zx24HhsTuDwFuS3Wc5cR9CNAN+GxLcQNHA68CBhwAfJjq+KtwLtcDl5dzbKfY71kDoF3s9y8n1eeQEN8uQLfY/SbAF7GYM+qzqeQ8Mu5zif1st4/drwd8GPtZPwOcFHv878D5sfsXAH+P3T8JeHpb3ztbWgRVme4iEyVO0fEocFzqQimfu79DGBGWqKK4jwUe8+ADYAczS5u14yo4l4ocC4x19/Xu/n/AAsLvYVpw92XuPiN2/0dgLqHSP6M+m0rOoyJp+7nEfrY/xTbrxW4OHEaYggc2/0ySMkVPtiSC8qa7qOyXJR058JqZTY9NuQGws7svi93/D7BzakLbahXFnamf00Wx7pIxCd1zGXMusS6F/QjfQDP2sylzHpCBn4uZ5ZjZTGA5MInQYvnB3YtjhyTGu8kUPUB8ip6tli2JoDY4yN27EWZzvdDMDknc6aF9mHFjgTM17gQPAB2ArsAy4K6URrOVzGx74FngUndfnbgvkz6bcs4jIz8Xdy9x966EmRh6AHvWxPtmSyKoynQXac3dl8b+XQ48T/gl+SbePI/9uzx1EW6ViuLOuM/J3b+J/efdCDxEaTdD2p+LmdUj/PF8wt2fiz2ccZ9NeeeRyZ8LgLv/AEwGfk3ohosX/ybGm7QperIlEVRluou0ZWbbmVmT+H2gN/AZm07RcSbwr9REuNUqinsCcEZshMoBwKqEboq0VKaf/HeEzwXCuZwUG9nRDsgDPqrp+CoS60v+BzDX3e9O2JVRn01F55GJn4uZtTCzHWL3GwFHEq55TCZMwQObfybJmaIn1VfKa+pGGPXwBaHP7epUx7OVsbcnjHSYBcyOx0/oD3wDmA+8DuyY6ljLif0pQtN8A6F/85yK4iaMmhgZ+4w+BfJTHX8VzuWfsVg/if3H3CXh+Ktj5zIP6Jvq+Mucy0GEbp9PgJmx29GZ9tlUch4Z97kA+wIfx2L+DLg29nh7QrJaAIwDGsQebxjbXhDb335b31tTTIiIZLls6RoSEZEKKBGIiGQ5JQIRkSynRCAikuWUCEREspwSgUg5zOxSM2tcyf7xZtY+dn97M3vAzL40sxmxaUDO28LrTy4782XsPR+IjSf/d3LORGTLlAhEyncpUG4iMLO9CTNWLow9NBr4HsjzMA1IH2DHLbz+U2y+NOtJwFPuvgJYZmY9tzF2ka2iRCBZy8zamtnnZvaEmc2NfctvbGYXAy2ByWY2uZynnkqsutPMOhCmLxjmYToD3H2Fu9+W8D5XmNm02ARoN8QeHg8ckzC3fNvYe06J7X8h9j4ikVMikGy3B3C/u+8FrAYucPcRwNfAoe5+aDnP6QlMj93fG5gVTwJlmVlvwjQGPQgToHU3s0Pc/TtCNWjf2KEnAc94aYVnAXBwdU9OpCqUCCTbLXH3d2P3HydMWbAluwArytthZlfHVsT6OvZQ79jtY2AGYTbJvNi+xO6hk2LbccsJLQSRyNXd8iEitVrZOVaqMufKOsI8LwBzgC5mVsfdN7r7TcBNZhZfYMSAW9z9wXJe51/APRaWfWzs7tMT9jWMvY9I5NQikGzX2sx+Hbt/CjA1dv9HwtKH5ZkL7A7g7gsI3Tg3mlkOhLVnCQkAYCJwdmy+fMxsVzPbKfbcnwgzS45h09YAQEdKZ8wUiZQSgWS7eYSFfuYC/0VY0ARgFPDvCi4Wvwz0Stg+lzBr5wIzKyCsLPVnAHd/DXgSeN/MPiVcJE5MME8BXdg8ERwaex+RyGn2UclasZE6L7n7Plv5vEaEb/I93b0kotjeAY519++jeH2RRGoRiGwld18HXEdEa92aWQvgbiUBqSlqEYiIZDm1CEREspwSgYhIllMiEBHJckoEIiJZTolARCTL/T8yZyIaTbzqmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(bins[:-1], hepdata, '-b', label=\"Obs\")\n",
    "plt.plot(bins[:-1], n, '--r', label='Sim')\n",
    "plt.ylim(0,0.75)\n",
    "plt.xlabel(\"pt (GeV)\")\n",
    "plt.ylabel(\"Fraction\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The weights are the difference between the histograms: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "weights = hepdata - n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.14285714,  0.14285714,  0.04761905, -0.04761905])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f66f28bfd30>"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAApxklEQVR4nO3deXhW1bn38e+dEAhICPOYgXmeE0gA66xQJ7SiMlm1Vazi8ajV1h57OZ3jUTv4eqzaiq2tyiACKjhXkdYiJJAwz5NkYAozAQxkWO8fzwYjTR5I8iQ7w+9zXbnYU559r+zAzdp7r3Wbcw4REZHShPkdgIiIVG9KFCIiEpQShYiIBKVEISIiQSlRiIhIUPX8DiDUWrZs6Tp27Oh3GCIiNUp6evo+51yrkvbVukTRsWNH0tLS/A5DRKRGMbOM0vbp1pOIiASlRCEiIkEpUYiISFC17hmFiAhAfn4+2dnZ5OXl+R1KtRIZGUlMTAwRERHn/D1KFCJSK2VnZxMVFUXHjh0xM7/DqRacc+zfv5/s7Gw6dep0zt+nW08iUivl5eXRokULJYlizIwWLVqUuZelRCEitZaSxL8rz89EiUJERIJSohARqSQPPPAAL7zwwun1kSNHcscdd5xe//nPf87zzz9f4vc+9thjfPHFF0E//4knnuB3v/vdv20/dOgQr7zySvmCLoEShYhIJRkxYgSLFi0CoKioiH379rF27drT+xctWsTw4cNL/N6nnnqKyy67rFznVaIQEakhhg8fzuLFiwFYu3Ytffv2JSoqioMHD3LixAnWr1+PmXHhhReSkJDAyJEj2bVrFwC33XYbs2fPBuDjjz+mZ8+eJCQkcN9993H11VefPse6deu46KKL6Ny5My+++CIAjzzyCFu3bmXgwIE8/PDDFW6HXo8VkVrvyQ/Wsm7nkZB+Zu/2TXj8mj5Bj2nfvj316tUjMzOTRYsWMWzYMHbs2MHixYuJjo6mV69ePPDAA8ydO5dWrVoxc+ZMHn30UV5//fXTn5GXl8ddd93FV199RadOnRg3btz3zrFhwwYWLFhAbm4uPXr04O677+bZZ59lzZo1rFixIiRtVaIQEalEw4cPZ9GiRSxatIgHH3yQHTt2sGjRIqKjo+nQoQN///vfufzyywEoLCykXbt23/v+DRs20Llz59PjHsaNG8eUKVNO77/qqqto0KABDRo0oHXr1uzZsyfkbVCiEJFa72z/869Mp55TrF69mr59+xIbG8vvf/97mjRpwkUXXXS6h1FeDRo0OL0cHh5OQUFBKML+Hj2jEBGpRMOHD+fDDz+kefPmhIeH07x5cw4dOsTixYsZN24ce/fuPZ0o8vPzv/ewG6BHjx5s27aN7du3AzBz5syznjMqKorc3NyQtUGJQkSkEvXr1499+/aRnJz8vW3R0dG0bt2a2bNn88tf/pIBAwYwcODA029JndKwYUNeeeUVRo0aRUJCAlFRUURHRwc9Z4sWLRgxYgR9+/YNycNsc85V+EOqk8TERKfCRSKyfv16evXq5XcYIXH06FEaN26Mc47JkyfTrVs3HnjggXJ/Xkk/GzNLd84llnS8rz0KMxtlZhvNbIuZPVLC/gfNbJ2ZrTKz+WYW70ecIiJ+eu211xg4cCB9+vTh8OHD3HXXXVV6ft8eZptZOPAycDmQDSw1s3nOuXXFDlsOJDrnjpvZ3cBvgJurPloREf888MADFepBVJSfPYqhwBbn3Dbn3EngbWB08QOccwucc8e91RQgpopjFJEarLbdWg+F8vxM/EwUHYCsYuvZ3rbS/BT4pKQdZjbJzNLMLG3v3r0hDFFEaqrIyEj279+vZFHMqXoUkZGRZfq+GjGOwswmAonAhSXtd85NAaZA4GF2FYYmItVUTEwM2dnZ6D+P33eqwl1Z+JkodgCxxdZjvG3fY2aXAY8CFzrnTlRRbCJSw0VERJSpipuUzs9bT0uBbmbWyczqA2OBecUPMLNBwKvAtc65HB9iFBGp83xLFM65AuBe4DNgPfCOc26tmT1lZtd6h/0WaAzMMrMVZjavlI8TEZFK4uszCufcx8DHZ2x7rNhy+SZjFxGRkNEUHiIiEpQShYiIBKVEISIiQSlRiIhIUEoUIiISlBKFiIgEpUQhIiJBKVGIiEhQShQiIhKUEoWIiASlRCEiIkEpUYiISFBKFCIiEpQShYiIBKVEISIiQSlRiIhIUEoUIiISlBKFiIgEpUQh1dpbi7fzs7fS+XLDHgqLnN/hiNRJvtbMFgnmi3V7eGzeWiLCw/h07W46NG3I+KQ4bkyMoXVUpN/hidQZ6lFItbR5Ty73z1xBvw7RpP/6Ml6ZMJiOLRvx2882MvyZL5k8fRmLtu7DOfUyRCqbehRS7Rw6fpI73kwjMiKcV29JICoygiv7tePKfu3YuvcoM1IzmZWezUerdtG51XlMSIpnzOAYohtF+B26SK1kte1/ZImJiS4tLc3vMKScCgqLuO2vS1nyzQFmTEomIb5Zicfl5Rfy0apdTE3NYHnmIRrUC+OaAe2ZmBzPgJhozKyKIxep2cws3TmXWNI+9SikWvnfjzewcMs+fjOmf6lJAiAyIpwbEmK4ISGGtTsPMy01k/eX72B2ejZ92jdhYnI81w5oz3kN9CsuUlHqUUi1MSsti4dnr+L2ER15/Jo+Zf7+3Lx83l+xk2kpGWzYnUtUg3pcP7gDE5Li6dE2qhIiFqk9gvUolCikWkjPOMi4KSkM7dScv90+hHrh5X/PwjnHssyDTE3J5KPVuzhZUMSQjs2YkBTPD/u1pUG98BBGLlI7KFFItbb7cB7XvLSQRvXDmTt5BE0b1Q/ZZx84dpLZ6VlMS80kY/9xmp9XnxsTYhifFEd8i/NCdh6Rmk6JQqqtvPxCbnp1MVtzjvLe5BF0b1M5t4iKihxfb93HtJRMPl8fGLz3g24tmZgcz6U9W1eoByNSG+hhtlRLzjkembOK1TsOM+WWxEpLEgBhYcYPurXiB91asftwHjOXZjFjSSZ3vZVO2yaRjB0ay9ghcbSN1kA+kTOpRyG+efWfW3nmkw08dEV37r2kW5Wfv6CwiC835DAtNZOvNu8lzIzLerVmQlI853dtSViYXrGVukM9Cql2FmzM4dlPN3BV/3ZMvrirLzHUCw/jij5tuaJPWzL3H2fakgxmpWXz2do9xLdoxPihcdyYGEvz80L3zESkJlKPQqrclpyjXP/y18Q2b8Tsu4fRqH71+f/KiYJCPl2zm2kpmSzZfoD64WFc2a8tE5PjSYhvpoF8UmupRyHVxuFv85n0Zhr164Xx2q2J1SpJADSoF87ogR0YPbADG3fnMj01g3eX7eD9FTvp2TaKCUlxXDeoA1GRmi5E6g71KKTKFBY5fvK3pSzauo9pdyQztFNzv0M6J8dOFPDByp1MTc1gzY4jNKofSCYTkuLo2yHa7/BEQiJYj8LXdwLNbJSZbTSzLWb2SAn7LzCzZWZWYGZj/IhRQue5Tzfwz017efLavjUmSQCc16AeY4fG8cG95zN38giu6teO95Znc/UfFnLdy18zKy2LvPxCv8MUqTS+9SjMLBzYBFwOZANLgXHOuXXFjukINAEeAuY552af7XPVo6ie3l2WzYPvrOSW5Hj++7q+fodTYYeP5zNnWTbTUjPYuvcYTSLrMSYhlgnJcXRp1djv8ETKrLo+oxgKbHHObQMws7eB0cDpROGc2+7tK/IjQAmNFVmHeOTd1SR3bs5j1/T2O5yQiG4UwU/O78TtIzqSsu0A01IzeCtlO69//Q3DOrdgYnI8l/duQ/16GsgnNZ+fiaIDkFVsPRtI8ikWqSQ5R/K46600Wkc14JUJCUTUshHQZsawLi0Y1qUFe3NP8E5aFtNTM5k8fRktGzdg7JBYxg6NJaZZI79DFSm36vXKSTmZ2SRgEkBcXJzP0cgpefmFTHorndy8AubcPbzWj0doFdWAyRd35WcXduGrTXuZmpLBy//Ywiv/2MLFPVozITmOC7u3JlwD+aSG8TNR7ABii63HeNvKzDk3BZgCgWcUFQ9NKso5x6PvrWFF1iH+NHEwvdo18TukKhMeZlzcszUX92xN9sHjvL0ki7eXZjH/b2mn637flBhLq6gGfocqck78vA+wFOhmZp3MrD4wFpjnYzwSQn9Z+A1zlmVz/2XdGNW3nd/h+CamWSMeGtmDxb+6hJfHDya+hVf3+9n5TJ6+jMVb96vut1R7vo6jMLMrgReAcOB159zTZvYUkOacm2dmQ4D3gGZAHrDbORe0oo3eevLfV5v2cttflzCyT1teHj9YcyadYeveo0xPzWR2ejaHv82ni1f3+wbV/RYfaZpxqTLf7DvG6JcW0r5pQ+bcPVylSIPIyy/kw1W7mJqSwYqsQ0RGhHFN//ZMUN1v8UF1fT1WapncvHzufDON8DDjtR8nKkmcRWREOGMSYhiTEMOaHYG633NX7GBWejZ9OzRhQlI8owe2r3bTnEjdox6FhERhkWPSm2n8c9Ne3vppEsO6tPA7pBopNy+f95fvYGpKJhv3qO63VB3depJK95tPN/DKP7by36P7cMuwjn6HU+M550jPOMi01Ew+WrWLk4WBut8Tk+MZ1Vd1vyX0lCikUs1buZP7Zixn3NA4/vf6vrq3HmIHjp1kVloW05cUq/udGMP4oar7LaGjRCGVZs2Ow4z50yL6d2jK1DuSNGVFJTpV93tqSgZfrM+hsMhxQfdWTEiKU91vqTAlCqkUe3NPcO1LCwkzY+69I2jZWAPIqsruw3m8vTSTt5dksftInup+S4UpUUjInSgoZMJrqazZeZg5dw+nT3vVZfBDQWER80/V/d60l/CwQN3vicnxjOiiut9y7vR6rISUc47H564lLeMgL40fpCTho3rhYYzs05aRfdqSsf8Y05dknq773bFFI8YnxTEmQXW/pWLUo5Aye2PRdh6ft5Z7L+7KQyN7+B2OnOFU3e+pKRks3X6Q+vXCuKpfOyYkxanut5RKt54kZBZt2cctry/h4h6tmXJLgm5tVHMbd+cyzav7ffREgep+S6mUKCQkMvcf59qXF9KqcQPevWe4/qGpQY6dKGDeyp1MTclg7c7v6n5PTI7TrUMBlCgkBI6eKOCGVxax+0ge8+4doff3ayjnHCuzDzM1JYMPVu7kREERA2ObMjE5nqv7tyMyQgP56iolCqmQoiLHz6amM39DDm/cPpTzu7X0OyQJgcPH85nt1f3etvcY0Q0jGJMQw/gk1f2ui5QopEKe/3wTL87fzGNX9+Yn53fyOxwJMecci7ftZ1pqJp+t2U1BkWN4l+/qfte28rVSMr0eK+X2yepdvDh/MzcmxHD7iI5+hyOVwMwY3qUlw7u0JCc3j1lp2UxPzeSeactoFdWAmxNjGZcUR4emDf0OVXyiHoWUat3OI9zwx0X0ahfFjEnJmoiuDikscvxzUw7TUjL5cmMOBlzcIzCQ74LurVT3uxbSrScps/1HT3DtS19TWOSY9x8jaB2laSHqquJ1v/cdPUFMs4aMG6q637WNEoWUSX5hERP/nMqKrEPM+tkw+sc09TskqQZOFhTx93W7mZaSyeJt+4kIN0b2acuEpHiSOzfXQL4aTs8opEye/GAtqd8c4P/GDlSSkNPq1wvj6v7tubp/e7bknKr7ncWHq3Z9V/c7IYbohhpfU9uoRyHfMzUlg1+/v4a7LuzMr37Yy+9wpJr79mQhH67aydTUTFYWq/s9MTmeAbFN/Q5PykC3nuScpGzbz8Q/p/KDbi35861D9MBSyiRQ9zuD95fv5Nv8Qvp1iGZCUhzXqu53jaBEIWeVdeA4o1/+mqaNInh/8giaaHoOKacjp+t+Z7Bpz1GiGtTjR4M7MCE5nu5tVPe7ulKikKCOnyzgR68sYsehb5k7eQSdNSpXQsA5R1rGQaalZPDx6t2cLCxiaMfmTEiOU93vaihYojinIZdm9ty5bJOaxznHQ7NWsmlPLn8YN0hJQkLGzBjSsTkvjB3E4l9dwq9+2JM9uXn859srGP7Mlzz7yQYy9x/3O0w5B+fUozCzZc65wWdsW+Wc619pkZWTehRl8+L8zTz/+Sb+68qeTLqgi9/hSC1XVORYuCVQ93v+hhyKnOOCboG635eo7revyv16rJndDdwDdDazVcV2RQFfhy5E8cNna3fz/OebuH5QB+78QWe/w5E6ICzMuKB7Ky7o3opdh7/1BvJlMumtdNpFRzJ2SBxjh8bSpokGeFYnQXsUZhYNNAOeAR4ptivXOXegkmMrF/Uozs3G3bn86JWv6dq6MTPvGqbppcU3BYVFfLE+h2mpGfxr8z7Cw4zLe7VhYnI8w7u0UHGsKhKSh9lmFg60oVgvxDmXGZIIQ0iJ4uwOHjvJtS8v5ER+EfPuPZ+20frfm1QP2/cdY8aSTN5Jy+Lg8Xw6tTyP8UPjGJMQQzPV/a5UFU4UZnYv8ASwByjyNjs9o6h58guLuPX1JaRlHGTmpGQGxTXzOySRf5OX/13d77SMQN3vq/u1Y0JyHIPjVPe7MoRiCo/7gR7Ouf0hi0p88fRH61m0dT+/v3GAkoRUW5ER4Vw3qAPXDerAht1HmJaSyXvLd/Du8h30bBvFxOR4rhvUgcYNNJCvKpxrj2IBcLlzrqDyQ6oY9ShKN3NpJr+cs5o7zu/Er6/u7Xc4ImVy7EQBc1cE6n6v23WE8+qHM3pQByYmxdO7fRO/w6vxyn3rycwe9Bb7AD2Aj4ATp/Y7554PYZwhoURRsrTtBxj3WgrJnVvw19uG6DVEqbGcc6zIOsS01MzTdb8HxTVlYlI8V6nud7lVJFE8HuyDnXNPVjC2kFOi+Hc7D33LtS8tJCoygvfvGUF0I03PIbXDoeMnmbNsx+m6300bRTBmcKDutwaPlo2m8KjDvj1ZyI2vLmL7vuO8P3k4XVtrrh2pfU7X/U7J5LO1gbrfI7q2YEKS6n6fqwo/zDazD4AzM8phIA141TmXV7EQpTI45/jFnFWs3XmEv9yaqCQhtdaZdb/fWZrFjCVZp+t+jx0Sy7ihcbRX3e9yOdc0uw04CrzmfR0BcoHu3nq5mNkoM9toZlvM7JES9jcws5ne/lQz61jec9VFf/znVj5YuZOHR/bgkp5t/A5HpEq0jork3ku68dUvLuYvtybSr0M0Ly3YwvnPfckdbyxlwcYcCotq152Uynau75YNd84NKbb+gZktdc4NMbO15TmxN4DvZeByIBtYambznHPrih32U+Cgc66rmY0FngNuLs/56pr56/fw2882cu2A9tx9oeZwkronPMy4tFcbLu3VhqwDx3l7aSYzl2bxxfocYpo1ZHxSoO53y8aq+30259qjaGxmcadWvOVTT4pOlvPcQ4EtzrltzrmTwNvA6DOOGQ284S3PBi41jbQ5qy05ufzn2yvo074Jz93QX4OTpM6Lbd6Ih0f2ZNEjl/LS+EHENGvIbz7dyLBn5vMfM5aTum0/te15bSida4/i58BCM9sKGNAJuMfMzuO7f8jLqgOQVWw9G0gq7RjnXIGZHQZaAPuKH2Rmk4BJAHFxcdRlh4/nc8cbaURGhDPllkQa1tergiKnnFn3e1pqBnPSs/lg5U66tW7MhKQ4rh+sut9nKstcTw2Ant7qxoo+wDazMcAo59wd3votQJJz7t5ix6zxjsn21rd6x+wr6TOhbr/1VFBYxO1/W0rKtv3MuDOZxI7N/Q5JpNr79mQhH6zayTSv7nfDiHCuHdCeCclx9I9p6nd4VaYi04xf4pz70sx+dMauLmaGc+7dCsS1A4gtth7jbSvpmGwzqwdEA5pGpBTPfrKBf23ex3M39FOSEDlHDeuHc1NiLDclxrI6O1D3e+6KncxMy6J/TKDu9zUD6nbd77MNuHvSOfe4mf21hN3OOfeTcp848A//JuBSAglhKTDeObe22DGTgX7OuZ95D7N/5Jy7Kdjn1tUexZz0bH4+ayW3De/IE9f28TsckRrtSF4+73kD+TbtOUpUZD1uGBzDhKQ4utXSut/VdsCdmV0JvACEA6875542s6eANOfcPDOLBN4CBgEHgLHOuW3BPrMuJorlmQe5eUoKifHNeOMnQzW4SCRETtX9npqSwSen6n53as7E5HhG9mlTq+p+h2Ka8TbA/wLtnXM/NLPewDDn3F9CG2rF1bVEsedIHtf8YSGREeHMnTxCc/aLVJL9R08wKz2b6amZZB44Tovz6nPTkFjGD40jtnkjv8OrsFAkik+AvwKPOucGeLeNljvn+oU21IqrS4kiL7+Qm6eksGVPLu/eM4IebWtnl1ikOikqcvxryz6mpWTwxfo9OODC7q2YkBTPJT1bE15DK/KFoh5FS+fcO2b2Kzj9qmphyCKUMnPO8V/vrmZl1iGm3JKgJCFSRcLCjAu7t+LCM+p+3/lmGu2jIxk3NI6bh8TSuhbV/T7XRHHMzFrgzfdkZskE5noSn/z5X9/w7vIdPHh5d67o09bvcETqpHbRDXng8u7ce0lX5nt1v3//+Sb+b/5mrujThglJgbrfNX3Q69lej70fWAT8ApgLdDazr4FWwI2VHp2U6J+b9vLMJ+u5sl9b/uOSrn6HI1LnRYSHMapvW0b1bcv2fceYviSTWWlZfLx6N51bnsf4pEDd76aNauYzxLO9Hvs7YDiBgXYbCLzG+hUwI9igNz/V9mcU2/YeZfTLXxPTrBFz7h5Wp9/tFqnO8vIL+WTNLqamZJJ+qu53/3ZMTI5nUGzTatfLCMXD7PpAIoGkMcz7OuScq3b1NGtzojiSl891L3/NoeP5zLt3BDHNav6bFiJ1wfpdR5iWmsF7y3Zw7GQhvdo1YWJyHKMHVp+636FIFNEEksMI78+mwGrn3O0hjDMkamuiKCxy/PSNpSzcvI9pdySR1LmF3yGJSBkdPVHA3BU7mJqSyfpdR2jcoB7XDWrPhKR4erXzt+53RUqhTiFQLzsXSAVSgBTn3MHKCDQUamuieOaT9bz6z238z3V9mZgc73c4IlIBzjmWZx1iWkomH64K1P0eHNeUicnxXNnPn7rfFUkUnwItgTUEHmovBta4ajwfb21MFO8v38H9M1cwISmOp6+vdkNXRKQCDh0/yWxvIN+2fYG63zcmxDA+KZ5OLc+rsjgqdOvJq//Qh8DzieFAXwLTaSx2zj0e4lgrrLYlilXZh7jxT4sZENuUqT9Non49Tc8hUhs551i8dT/TUr+r+31+15ZMTI7j0l6VX/c7JHM9mVkMgWcUw4GrgRbOuaahCjJUalOiyDmSx7UvfU14mDHv3hG0UCUukToh50ge76QF6n7vOPQtraMaMHZoHGOHxFZa3e+K3Hq6j+96EvkEbj+d+lrtnCsKfbgVU1sSxYmCQsZNSWH9rlzm3D2c3u39fdAlIlWvsMixYENgIN8/Nu3FgEt7tWFCUhwXdGtFWAinC6nIFB4dgVnAA865XSGLSIJyzvHr99awLPMQr0wYrCQhUkeFhxmX9W7DZb0Ddb9nLMnknbQsPl+3h9jmDRk/NJ6bEmMq/W6Dr9OMV4ba0KN4feE3PPXhOu67pCsPXtHD73BEpBo5WVDEZ2t3MzUlg9RvDlDfGxU+MTmeIR2blXsgXygmBZQqsnDzPp7+eD1X9G7D/Zd19zscEalm6tcL45oB7blmQHu25OQyNSWTOcuymbdyJ0M7NWfmpOSQj/pWoqhGtu87xuTpy+jaqjHP3zwwpPcfRaT26do6iieu7cMvR/Xkg1U7ycsvrJSpQZQoqoncvHzufDMNM3jtx4nVZli/iFR/p+p+Vxb9a1QNFBU5Hpi5gm37jvHWT4YS10JzOIlI9aHRW9XA859v4ov1OTx2dW+Gd23pdzgiIt+jROGzD1ft5KUFWxg7JJYfD9McTiJS/ShR+GjNjsM8NGslifHNeGp032o3P72ICChR+Gbf0RNMejON5o3q88eJCZrDSUSqLT3M9sHJgiLunprOgeMnmf2z4bSK0hxOIlJ9KVFUMeccj89by9LtB/nDuEH07RDtd0giIkHpfkcVm5qSwYwlmdxzUReuGdDe73BERM5KiaIKLd66nyc/WMelPVvzkOZwEpEaQomiimQdOM4909Lp2PI8Xhir6TlEpOZQoqgCx04UcOebaRQWOV77cSJRkRF+hyQics70MLuSFRU5fv7OSjbtyeVvtw+t0hq4IiKhoB5FJXvxy818unY3/3VlLy7o3srvcEREykyJohJ9umYXL3yxmRsGx/DT8zv5HY6ISLkoUVSSDbuP8OA7KxkY25Snr9f0HCJScylRVIIDx05yxxtpREXWY8otCURGhPsdkohIuelhdojlFxZxz7R0cnJPMOuuYbRuEul3SCIiFeJLj8LMmpvZ52a22fuzWSnHfWpmh8zsw6qOsbz++8N1pGw7wHM39GNAbFO/wxERqTC/bj09Asx3znUD5nvrJfktcEuVRVVB01MzeXNxBpMu6Mz1g2L8DkdEJCT8ShSjgTe85TeA60o6yDk3H8itopgqZMk3B3hs7hou7N6KX47q6Xc4IiIh41eiaOOc2+Ut7wba+BRHSGQfPM7dU9OJa96IF8cNIlzTc4hILVJpD7PN7AugbQm7Hi2+4pxzZuYqeK5JwCSAuLi4inxUmR0/WcCkN9M5WVDEa7cmEt1Q03OISO1SaYnCOXdZafvMbI+ZtXPO7TKzdkBOBc81BZgCkJiYWKGkU8bz8vCsVazffYTXbx1Cl1aNq+rUIiJVxq9bT/OAW73lW4G5PsVRIS8v2MJHq3fxy1E9ubhna7/DERGpFH4limeBy81sM3CZt46ZJZrZn08dZGb/AmYBl5pZtpmN9CXaEny+bg+/+/smrhvYnrsu6Ox3OCIilcaXAXfOuf3ApSVsTwPuKLb+g6qM61xt2pPL/W8vp39MNM/e0F/Tc4hIraYpPMro0PGT3PlmGg3r1+NVTc8hInWAEkUZFBQWce/05ew6lMertyTQLrqh3yGJiFQ6zfVUBk9/vJ6FW/bxmzH9SYgvcdYREZFaRz2Kc/ROWhZ//Xo7t4/oyE2JsX6HIyJSZZQozkF6xkF+/d4azu/akkev7OV3OCIiVUqJ4ix2Hf6Wu95Kp13TSF4aP4h64fqRiUjdomcUQeTlFzLpzXTy8guZcWcSTRvV9zskEZEqp0RRCuccv5yzijU7D/PaLYl0axPld0giIr7QfZRSvPrVNuau2MlDV/Tgst41enJbEZEKUaIowYINOTz36Qau7t+Oey7q4nc4IiK+UqI4w5aco9w3Yzm92zXht2MGaHoOEanzlCiKOfxtPpPeTKNBRBhTfpxIw/qankNERA+zPYVFjvtmLCfr4HGm35lMh6aankNEBNSjOC3zwHFW7zjMU6P7MqRjc7/DERGpNtSj8HRqeR5f/vxCjZUQETmDehTFKEmIiPw7JQoREQlKiUJERIJSohARkaCUKEREJCglChERCUqJQkREglKiEBGRoJQoREQkKCUKEREJSolCRESCUqIQEZGglChERCQoJQoREQlKiUJERIJSohARkaCUKEREJCglChERCUqJQkREglKiEBGRoHxJFGbW3Mw+N7PN3p/NSjhmoJktNrO1ZrbKzG72I1YRkbrOrx7FI8B851w3YL63fqbjwI+dc32AUcALZta06kIUERHwL1GMBt7wlt8ArjvzAOfcJufcZm95J5ADtKqqAEVEJMCvRNHGObfLW94NtAl2sJkNBeoDW0vZP8nM0swsbe/evaGNVESkjqtXWR9sZl8AbUvY9WjxFeecMzMX5HPaAW8Btzrniko6xjk3BZgCkJiYWOpniYhI2VVaonDOXVbaPjPbY2btnHO7vESQU8pxTYCPgEedcymVFKqIiATh162necCt3vKtwNwzDzCz+sB7wJvOudlVGJuIiBTjV6J4FrjczDYDl3nrmFmimf3ZO+Ym4ALgNjNb4X0N9CVaEZE6zJyrXbf0ExMTXVpamt9hiIjUKGaW7pxLLGmfRmaLiEhQShQiIhKUEoWIiASlRCEiIkEpUYiISFBKFCIiEpQShYiIBKVEISIiQSlRiIhIUEoUIiISlBKFiIgEpUQhIiJBKVGIiEhQShQiIhKUEoWIiASlRCEiIkHVusJFZrYXyKjAR7QE9oUoHD/VlnaA2lJd1Za21JZ2QMXaEu+ca1XSjlqXKCrKzNJKq/JUk9SWdoDaUl3VlrbUlnZA5bVFt55ERCQoJQoREQlKieLfTfE7gBCpLe0AtaW6qi1tqS3tgEpqi55RiIhIUOpRiIhIUEoUIiISlBKFx8xGmdlGM9tiZo/4HU9Zmdl2M1ttZivMLM3b1tzMPjezzd6fzfyOsyRm9rqZ5ZjZmmLbSozdAl70rtMqMxvsX+TfV0o7njCzHd51WWFmVxbb9yuvHRvNbKQ/UZfMzGLNbIGZrTOztWb2n972mnhdSmtLjbo2ZhZpZkvMbKXXjie97Z3MLNWLd6aZ1fe2N/DWt3j7O5b75M65Ov8FhANbgc5AfWAl0NvvuMrYhu1AyzO2/QZ4xFt+BHjO7zhLif0CYDCw5myxA1cCnwAGJAOpfsd/lnY8ATxUwrG9vd+zBkAn7/cv3O82FIuvHTDYW44CNnkx18TrUlpbatS18X62jb3lCCDV+1m/A4z1tv8JuNtbvgf4k7c8FphZ3nOrRxEwFNjinNvmnDsJvA2M9jmmUBgNvOEtvwFc518opXPOfQUcOGNzabGPBt50ASlAUzNrVyWBnkUp7SjNaOBt59wJ59w3wBYCv4fVgnNul3NumbecC6wHOlAzr0tpbSlNtbw23s/2qLca4X054BJgtrf9zGty6lrNBi41MyvPuZUoAjoAWcXWswn+i1QdOeDvZpZuZpO8bW2cc7u85d1AG39CK5fSYq+J1+pe73bM68Vu/9WYdni3LAYR+B9sjb4uZ7QFati1MbNwM1sB5ACfE+jtHHLOFXiHFI/1dDu8/YeBFuU5rxJF7XG+c24w8ENgspldUHynC/Q/a+S70DU5duCPQBdgILAL+L2v0ZSRmTUG5gD3O+eOFN9X065LCW2pcdfGOVfonBsIxBDo5fSsivMqUQTsAGKLrcd422oM59wO788c4D0Cv0R7TnX/vT9z/IuwzEqLvUZdK+fcHu8vdxHwGt/dwqj27TCzCAL/sE5zzr3rba6R16WkttTka+OcOwQsAIYRuM1Xz9tVPNbT7fD2RwP7y3M+JYqApUA37+2B+gQe/MzzOaZzZmbnmVnUqWXgCmANgTbc6h12KzDXnwjLpbTY5wE/9t6ySQYOF7sVUu2ccZ/+egLXBQLtGOu9mdIJ6AYsqer4SuPdy/4LsN4593yxXTXuupTWlpp2bcyslZk19ZYbApcTeN6yABjjHXbmNTl1rcYAX3q9wLLz+0l+dfki8NbGJgL3/B71O54yxt6ZwFsaK4G1p+IncD9yPrAZ+AJo7nespcQ/g0DXP5/APdaflhY7gTc/Xvau02og0e/4z9KOt7w4V3l/cdsVO/5Rrx0bgR/6Hf8ZbTmfwG2lVcAK7+vKGnpdSmtLjbo2QH9guRfvGuAxb3tnAolsCzALaOBtj/TWt3j7O5f33JrCQ0REgtKtJxERCUqJQkREglKiEBGRoJQoREQkKCUKEREJSolCpBzM7H4zaxRk/2wz6+wtNzazP5rZVjNb5k2zcudZPn/BmbOWeuf8o/c+/aehaYnI2SlRiJTP/UCJicLM+hCYbXSbt+nPwEGgmwtMszIKaH6Wz59BYOBncWOBGc65vcAuMxtRzthFykSJQqQUZtbRzDaY2TQzW+/1EhqZ2X1Ae2CBmS0o4Vsn4I2ONbMuBKaG+LULTBWBc26vc+65Yud52MyWepPTPeltng1cVay2QEfvnP/y9r/vnUek0ilRiATXA3jFOdcLOALc45x7EdgJXOycu7iE7xkBpHvLfYCVp5LEmczsCgJTRAwlMDldgpld4Jw7QGA07Q+9Q8cC77jvRsimAT+oaONEzoUShUhwWc65r73lqQSmgzibdsDeknaY2aNeNbWd3qYrvK/lwDICs4F28/YVv/001ls/JYdAD0Ok0tU7+yEiddqZc9ycy5w33xKYZwdgHTDAzMKcc0XOuaeBp83sVAEaA55xzr1awufMBf6fBcqKNnLOpRfbF+mdR6TSqUchElycmQ3zlscDC73lXAJlNUuyHugK4JzbQuA20f+YWTgEah8TSBAAnwE/8WolYGYdzKy1971HCcwM+jrf700AdOe72U5FKpUShUhwGwkUgloPNCNQ7AZgCvBpKQ+zPwIuKrZ+B4FZV7eYWRqBymS/AHDO/R2YDiw2s9UEHmIXT0AzgAH8e6K42DuPSKXT7LEipfDeNPrQOde3jN/XkEBPYIRzrrCSYvsKGO2cO1gZny9SnHoUIiHmnPsWeJxKqrNsZq2A55UkpKqoRyEiIkGpRyEiIkEpUYiISFBKFCIiEpQShYiIBKVEISIiQf1/81FbsmCkIk8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(bins[:-1], weights, '-', label=\"Weight\")\n",
    "plt.ylim(-0.25, 0.25)\n",
    "plt.xlabel(\"pt (GeV)\")\n",
    "plt.ylabel(\"Weight\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Then you broadcast the weights to each event based on the category"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.14285714, -0.14285714,  0.14285714,  0.14285714,  0.14285714,\n",
       "        0.14285714,  0.14285714,  0.04761905,  0.04761905,  0.04761905,\n",
       "        0.04761905, -0.04761905, -0.04761905, -0.04761905])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "weights[categories]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}