wiso/efficiency tpr fpr roc curve.ipynb

## efficiency tpr fpr roc curve.ipynb
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sample</th>\n",
       "      <th>ph_cl_E</th>\n",
       "      <th>el_cl_E</th>\n",
       "      <th>y_pred_pt_eta</th>\n",
       "      <th>el_amb_type</th>\n",
       "      <th>ph_amb_type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>ph</td>\n",
       "      <td>9.359830e+04</td>\n",
       "      <td>9.427870e+04</td>\n",
       "      <td>0.961212</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>el</td>\n",
       "      <td>8.457673e+04</td>\n",
       "      <td>8.554388e+04</td>\n",
       "      <td>0.104820</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>el</td>\n",
       "      <td>4.567817e+04</td>\n",
       "      <td>4.595964e+04</td>\n",
       "      <td>0.002681</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>el</td>\n",
       "      <td>4.142598e+05</td>\n",
       "      <td>4.122794e+05</td>\n",
       "      <td>0.005727</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>ph</td>\n",
       "      <td>1.646737e+04</td>\n",
       "      <td>1.570206e+04</td>\n",
       "      <td>0.920399</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>el</td>\n",
       "      <td>1.125313e+06</td>\n",
       "      <td>1.124633e+06</td>\n",
       "      <td>0.016469</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>ph</td>\n",
       "      <td>8.486954e+04</td>\n",
       "      <td>8.528577e+04</td>\n",
       "      <td>0.972340</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>ph</td>\n",
       "      <td>6.920021e+05</td>\n",
       "      <td>6.890341e+05</td>\n",
       "      <td>0.966844</td>\n",
       "      <td>5</td>\n",
       "      <td>5</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>ph</td>\n",
       "      <td>2.373500e+05</td>\n",
       "      <td>2.344050e+05</td>\n",
       "      <td>0.962757</td>\n",
       "      <td>2</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>el</td>\n",
       "      <td>2.644520e+05</td>\n",
       "      <td>2.683211e+05</td>\n",
       "      <td>0.133515</td>\n",
       "      <td>3</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  sample       ph_cl_E       el_cl_E  y_pred_pt_eta  el_amb_type  ph_amb_type\n",
       "0     ph  9.359830e+04  9.427870e+04       0.961212            2            2\n",
       "1     el  8.457673e+04  8.554388e+04       0.104820            0            0\n",
       "2     el  4.567817e+04  4.595964e+04       0.002681            0            0\n",
       "3     el  4.142598e+05  4.122794e+05       0.005727            0            0\n",
       "4     ph  1.646737e+04  1.570206e+04       0.920399            5            5\n",
       "5     el  1.125313e+06  1.124633e+06       0.016469            0            0\n",
       "6     ph  8.486954e+04  8.528577e+04       0.972340            2            2\n",
       "7     ph  6.920021e+05  6.890341e+05       0.966844            5            5\n",
       "8     ph  2.373500e+05  2.344050e+05       0.962757            2            2\n",
       "9     el  2.644520e+05  2.683211e+05       0.133515            3            3"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = pd.read_pickle(\"Downloads/df_diff.pkl\")\n",
    "df.head(10)[['sample', 'ph_cl_E', 'el_cl_E', 'y_pred_pt_eta', 'el_amb_type', 'ph_amb_type']]\n",
    "\n",
    "# as expected el_amb_type == ph_amb_type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAQ+0lEQVR4nO3df5BdZX3H8ffXJLCt5YcmSybDghtHBDIyAq4pDoxVUxlGOoaxmYAFuy1po6E67diZNq1/lLbMiDNWS5nYNgOWWEBD01oy2traSIaBAXQDUSMBi3SDSwNZA6S2GCXw7R/3gNfNbu7J3h+7z+b9mtnZ83Of77N388mzzz3nbGQmkqTyvGqmC5AkTY8BLkmFMsAlqVAGuCQVygCXpELN72VjixYtysHBwV42KUnF27Fjxw8ys3/i9p4G+ODgICMjI71sUpKKFxF7JtvuFIokFcoAl6RCGeCSVKiezoFLUje98MILjI2NcfDgwZkuZVr6+voYGBhgwYIFtY43wCXNGWNjY5xwwgkMDg4SETNdzlHJTPbv38/Y2BhLly6tdY5TKJLmjIMHD7Jw4cLiwhsgIli4cOFR/fZggEuaU0oM75cdbe0GuCQVyjlwSXPW4Povd/TrjV5/6fTqqG5iXLRoUUfrKSbAm1+I6X4TJWkucQpFkjro1ltvZfny5Zx77rl88IMf5MUXX+xaWwa4JHXI7t272bx5M/feey87d+5k3rx53HbbbV1rr5gpFEma7bZt28aOHTt461vfCsCPfvQjTjnllK61Z4BLUodkJsPDw3z84x//me233HJLV9pzCkWSOmTFihVs2bKFffv2AfDMM8+wZ8+kT4LtCEfgkuasXl+xtmzZMq677jouvvhiXnrpJRYsWMCGDRu61l6tAI+Ik4GbgDcBCVwNPApsBgaBUWB1Zj7blSolqRCXX345l19++c9sGx0d7UpbdadQbgC+kplnAW8GdgPrgW2ZeQawrVqXJPVIywCPiJOAtwM3A2TmTzLzOWAlsKk6bBNwWbeKlCQdrs4IfCkwDvxdRDwUETdFxKuBxZm5tzrmKWDxZCdHxNqIGImIkfHx8c5ULUmqFeDzgfOBv87M84D/Y8J0SWYmjbnxw2Tmxswcysyh/v7D/qiyJGma6gT4GDCWmQ9U61toBPrTEbEEoPq8rzslSpIm0zLAM/Mp4PsRcWa1aQXwMLAVGK62DQN3dqVCSdKk6l4H/hHgtog4Dngc+E0a4X9HRKwB9gCru1OiJE3TtSd1+OsdmPap3XikbK0Az8ydwNAku1Z0rBJJ0lHxVnpJ6qDR0VHOOussrrzySs4++2xWrVrF888/D8CNN97I+eefzznnnMMjjzzSdlsGuCR12KOPPso111zD7t27OfHEE/nMZz4DwKJFi3jwwQdZt24dn/zkJ9tuxwCXpA477bTTuPDCCwG46qqruOeeewB43/veB8Bb3vKWjtxeb4BLUodN/OvyL68ff/zxAMybN49Dhw613Y4BLkkd9sQTT3DfffcBcPvtt3PRRRd1pR0fJytp7mrjsr92nHnmmWzYsIGrr76aZcuWsW7dOm688caOt2OAS1KHzZ8/n1tvvfVntjXPeQ8NDbF9+/a223EKRZIKZYBLUgcNDg6ya9eunrRlgEuaUxoPRy3T0dZugEuaM/r6+ti/f3+RIZ6Z7N+/n76+vtrn+CampDljYGCAsbExSv3jMX19fQwMDNQ+3gCXNGcsWLCApUuXznQZPeMUiiQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RC1bqVPiJGgR8CLwKHMnMoIl4LbAYGgVFgdWY+250yJUkTHc0I/J2ZeW5mDlXr64FtmXkGsK1alyT1SDtTKCuBTdXyJuCy9suRJNVVN8AT+PeI2BERa6ttizNzb7X8FLC449VJkqZU93GyF2XmkxFxCvDViHikeWdmZkRM+gT1KvDXApx++ultFStJ+qlaI/DMfLL6vA/4IrAceDoilgBUn/dNce7GzBzKzKH+/v7OVC1Jah3gEfHqiDjh5WXgYmAXsBUYrg4bBu7sVpGSpMPVmUJZDHwxIl4+/vbM/EpEfAO4IyLWAHuA1d0rU5I0UcsAz8zHgTdPsn0/sKIbRUmSWvNOTEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqFqB3hEzIuIhyLiS9X60oh4ICIei4jNEXFc98qUJE10NCPw3wV2N61/Avh0Zr4BeBZY08nCJElHVivAI2IAuBS4qVoP4F3AluqQTcBl3ShQkjS5uiPwvwT+AHipWl8IPJeZh6r1MeDUyU6MiLURMRIRI+Pj420VK0n6qZYBHhG/AuzLzB3TaSAzN2bmUGYO9ff3T+dLSJImMb/GMRcC742I9wB9wInADcDJETG/GoUPAE92r0xJ0kQtR+CZ+UeZOZCZg8AVwNcy80rgLmBVddgwcGfXqpQkHaad68D/EPhoRDxGY0785s6UJEmqo84UyisyczuwvVp+HFje+ZIkSXV4J6YkFcoAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQLQM8Ivoi4usR8c2I+E5E/Gm1fWlEPBARj0XE5og4rvvlSpJeVmcE/mPgXZn5ZuBc4JKIuAD4BPDpzHwD8CywpntlSpImahng2fC/1eqC6iOBdwFbqu2bgMu6UqEkaVK15sAjYl5E7AT2AV8Fvgc8l5mHqkPGgFOnOHdtRIxExMj4+HgnapYkUTPAM/PFzDwXGACWA2fVbSAzN2bmUGYO9ff3T7NMSdJER3UVSmY+B9wFvA04OSLmV7sGgCc7XJsk6QjqXIXSHxEnV8s/B7wb2E0jyFdVhw0Dd3arSEnS4ea3PoQlwKaImEcj8O/IzC9FxMPAFyLiOuAh4OYu1ilJmqBlgGfmt4DzJtn+OI35cEnSDPBOTEkqlAEuSYUywCWpUAa4JBXKAJekQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFcoAl6RCGeCSVCgDXJIKZYBLUqHmtzogIk4DPgcsBhLYmJk3RMRrgc3AIDAKrM7MZ7tXqiTNctee1LR8oOvN1RmBHwJ+PzOXARcAvxMRy4D1wLbMPAPYVq1LknqkZYBn5t7MfLBa/iGwGzgVWAlsqg7bBFzWrSIlSYc7qjnwiBgEzgMeABZn5t5q11M0plgmO2dtRIxExMj4+HgbpUqSmtUO8Ij4BeAfgd/LzP9p3peZSWN+/DCZuTEzhzJzqL+/v61iJUk/VSvAI2IBjfC+LTP/qdr8dEQsqfYvAfZ1p0RJ0mRaBnhEBHAzsDszP9W0ayswXC0PA3d2vjxJ0lRaXkYIXAh8APh2ROystv0xcD1wR0SsAfYAq7tToiRpMi0DPDPvAWKK3Ss6W44kqS7vxJSkQhngklQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVAGuCQVygCXpEIZ4JJUKANckgplgEtSoQxwSSqUAS5JhTLAJalQBrgkFarOX+SZdQbXf/mV5dHrL53BSiRp5jgCl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYUq8jJCSZo1rj1pxppuOQKPiM9GxL6I2NW07bUR8dWI+M/q82u6W6YkaaI6Uyi3AJdM2LYe2JaZZwDbqnVJUg+1DPDMvBt4ZsLmlcCmankTcFmH65IktTDdNzEXZ+beavkpYPFUB0bE2ogYiYiR8fHxaTYnSZqo7atQMjOBPML+jZk5lJlD/f397TYnSapMN8CfjoglANXnfZ0rSZJUx3QDfCswXC0PA3d2phxJUl0trwOPiM8D7wAWRcQY8CfA9cAdEbEG2AOs7maRkjSrzOC1381aBnhmvn+KXSs6XIskzR3NIX/tga404a30klQoA1ySCmWAS1KhDHBJKpQBLkmFMsAlqVA+D1yS6pgl1343cwQuSYUywCWpUMVPoQyu//Iry6PXXzqDlUhSbxUf4JLUNbNw3ruZUyiSVCgDXJIKZYBLUqGcA5ekZrN83ruZI3BJKpQjcEkqaNTdzBG4JBVqTo3AvalHUm2FjrqbOQKXpELNqRG4JB3RHBh1N5uzAe50iqS5bs4GuKRj2BwbaU/lmAhwR+PSHHeMBPZEbQV4RFwC3ADMA27KzOs7UlUXNYc5GOjSrNEcwtceaH2Mph/gETEP2AC8GxgDvhERWzPz4U4V12y079cm3T548Pa2vu7EQG9Zh4E/N0wVBFMFx2w0W/pQp46jDV6DupZ2RuDLgccy83GAiPgCsBLoSoBPZapgb9Yc8nWOP9L5dX6w6rRX9z+e5vOnOmeqY472XDE3gmO29GG21DGHRWZO78SIVcAlmflb1foHgF/MzA9POG4tsLZaPRN4dJq1LgJ+MM1zS2Wfjw32ee5rt7+vy8z+iRu7/iZmZm4ENrb7dSJiJDOHOlBSMezzscE+z33d6m87d2I+CZzWtD5QbZMk9UA7Af4N4IyIWBoRxwFXAFs7U5YkqZVpT6Fk5qGI+DDwbzQuI/xsZn6nY5Udru1pmALZ52ODfZ77utLfab+JKUmaWT6NUJIKZYBLUqFmXYBHxCUR8WhEPBYR6yfZf3xEbK72PxARg72vsrNq9PmjEfFwRHwrIrZFxOtmos5OatXnpuN+NSIyIoq+5KxOfyNidfU6fyci2rvFeBao8XN9ekTcFREPVT/b75mJOjspIj4bEfsiYtcU+yMi/qr6nnwrIs5vq8HMnDUfNN4M/R7weuA44JvAsgnHXAP8TbV8BbB5puvuQZ/fCfx8tbzuWOhzddwJwN3A/cDQTNfd5df4DOAh4DXV+ikzXXcP+rwRWFctLwNGZ7ruDvT77cD5wK4p9r8H+FcggAuAB9ppb7aNwF+5PT8zfwK8fHt+s5XApmp5C7AiIqKHNXZayz5n5l2Z+Xy1ej+Na+5LVud1Bvhz4BPAwV4W1wV1+vvbwIbMfBYgM/f1uMZOq9PnBE6slk8C/ruH9XVFZt4NPHOEQ1YCn8uG+4GTI2LJdNubbQF+KvD9pvWxatukx2TmIeAAsLAn1XVHnT43W0Pjf/CStexz9avlaZl5dE8bm53qvMZvBN4YEfdGxP3Vkz5LVqfP1wJXRcQY8C/AR3pT2ow62n/vR3RMPA98roiIq4Ah4JdmupZuiohXAZ8CfmOGS+ml+TSmUd5B4zesuyPinMx8bkar6q73A7dk5l9ExNuAv4+IN2XmSzNdWClm2wi8zu35rxwTEfNp/Oq1vyfVdUetRxJExC8DHwPem5k/7lFt3dKqzycAbwK2R8QojbnCrQW/kVnnNR4DtmbmC5n5X8B3aQR6qer0eQ1wB0Bm3gf00Xjo01zW0UeQzLYAr3N7/lZguFpeBXwtq3cHCtWyzxFxHvC3NMK79LlRaNHnzDyQmYsyczAzB2nM+783M0dmpty21fm5/mcao28iYhGNKZXHe1lkh9Xp8xPACoCIOJtGgI/3tMre2wr8enU1ygXAgczcO+2vNtPv2k7xLu13abyD/bFq25/R+AcMjRf5H4DHgK8Dr5/pmnvQ5/8AngZ2Vh9bZ7rmbvd5wrHbKfgqlJqvcdCYNnoY+DZwxUzX3IM+LwPupXGFyk7g4pmuuQN9/jywF3iBxm9Va4APAR9qep03VN+Tb7f7c+2t9JJUqNk2hSJJqskAl6RCGeCSVCgDXJIKZYBLUqEMcEkqlAEuSYX6f1YpsubBj9V2AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "for k, v in df.groupby('sample')['y_pred_pt_eta']:\n",
    "    v.hist(bins=np.linspace(0, 1, 100), density=True, label=k, grid=False)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider as signal \"photon\". We define the photon selection: `y_pred_pt_eta > threshold`. This means that the efficiency is $\\int_{cut value}^\\infty \\frac{dP}{d\\texttt{y_pred_pt_eta}}\\,d\\texttt{y_pred_pt_eta}$ where $\\frac{dP}{d\\texttt{y_pred_pt_eta}}$ if the pdf of `y_pred_pt_eta`. We have the histogram of it from the test-sample."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "thresholds = np.linspace(0, 1, 100 + 1)\n",
    "df_true_ph = df.query('sample == \"ph\"')\n",
    "df_true_el = df.query('sample == \"el\"')\n",
    "\n",
    "# signal = ph\n",
    "# the efficiency is just the integral of the prediction\n",
    "# eff = integral_cutvalue^inf\n",
    "h_ph, _ = np.histogram(df_true_ph['y_pred_pt_eta'], bins=thresholds)\n",
    "h_el, _ = np.histogram(df_true_el['y_pred_pt_eta'], bins=thresholds)\n",
    "\n",
    "cdf_ph, cdf_el = np.ones(len(thresholds)), np.ones(len(thresholds))\n",
    "\n",
    "cdf_ph[1:] = np.cumsum(h_ph) / len(df_true_ph)\n",
    "cdf_el[1:] = np.cumsum(h_el) / len(df_true_el)\n",
    "\n",
    "effs_ph = 1 - cdf_ph\n",
    "effs_el = 1 - cdf_el"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.         0.99459581 0.9894168  0.98581401 0.98378744 0.98176086\n",
      " 0.98086017 0.97995947 0.9788336  0.97680703]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.9894167980184643"
      ]
     },
     "execution_count": 127,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(effs_ph[:10])\n",
    "\n",
    "# check manually one point, for example if the threshold is at 0.02 (the tested thresholds are 0, 0.01, 0.02, ...)\n",
    "(df_true_ph['y_pred_pt_eta'] > 0.02).sum() / len(df_true_ph)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import roc_curve\n",
    "fpr, tpr, _ = roc_curve((df['sample'] == 'ph').astype(int), df['y_pred_pt_eta'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'P[reco-ph BDT | true-el] = fpr')"
      ]
     },
     "execution_count": 90,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAF0CAYAAAAuMT6NAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZhcVZ3/8fenqjtLZycJEBKyYVgiWyBiGB1wQQSeERRcQBk3FHVcBlFn9KcPo/jMz1Fn/KmjMyPKosgiKONkJIgKifooQRICSALRbJ2EYBKSzh6S7q7v7497u1Pp9FJJ6lZVpz6v56mn6t46dc/3ZrnfOufUPUcRgZmZ1a9ctQMwM7PqciIwM6tzTgRmZnXOicDMrM45EZiZ1TknAjOzOtdQ7QAOxZgxY2Ly5MnVDsPMrN9YuHDhCxExtrv3+mUimDx5MgsWLKh2GGZm/Yak5p7ec9eQmVmdcyIwM6tzTgRmZnXOicDMrM45EZiZ1TknAjOzOudEYGZW55wIzMzqnBOBmVmdyzQRSLpF0gZJT/fwviR9U9IySU9JOivLeMzM7EBZtwhuAy7q5f2LgWnp41rgPzOOx8zMusg0EUTEb4DNvRS5DPhBJOYDIyWNyzImMzPbX7XHCMYDa4q216b7zMys2G+/Bj+8IpNDVzsRlEzStZIWSFqwcePGaodjZlZZW1bD809lcuhqJ4LngOOLtiek+w4QETdFxMyImDl2bLdTapuZ2SGodiKYDbwz/fXQLGBrRDxf5ZjMzOpKpgvTSLoLeBUwRtJa4J+ARoCI+C9gDnAJsAzYBbwny3jMzOxAmSaCiLiqj/cD+HCWMZiZWe+q3TVkZmZV5kRgZlbnnAjMzOqcE4GZWZ1zIjAzq3NOBGZmdc6JwMyszjkRmJnVOScCM7M650RgZlbnnAjMzOqcE4GZWZ1zIjAzq3NOBGZmdc6JwMyszjkRmJn1B1EAKZNDOxGYmfUHhXbINWZyaCcCM7P+oNAKuXwmh3YiMDPrDwptkHeLwMysfrW3Qi6bZeadCMzM+gOPEZiZ1TmPEZiZ1TmPEZiZ1blCOyibS7YTgZlZvxCAbygzM6tfEb6z2MzMnAjMzOqXWwRmZvXOYwRmZvXNLQIzs3oXmR3ZicDMrD9wi8DMzDxGYGZW19wiMDOrb+FfDZmZ1Tm3CMzM6ptbBGZm9c4tAjOz+uYWgZlZvXOLwMysvrlFYGZmWXEiMDPrDwptkGvI5NBOBGZm/UHrbmgclMmhnQjMzPqDQivkGjM5tBOBmVl/4V8NmZnVseyWI3AiMDPrP9wiMDOrY16hzMzM+uMYgaSLJC2VtEzSp7t5f6KkuZIWSXpK0iVZxmNm1m9FP2wRSMoD3wYuBqYDV0ma3qXY54B7ImIGcCXwH1nFY2bW//W/FsE5wLKIWBERe4G7gcu6lAlgePp6BLAuw3jMzPqx7FoE2dyvnBgPrCnaXgu8vEuZzwO/kPRRYAhwQYbxmJn1b9k0CKo+WHwVcFtETAAuAW6X1G1Mkq6VtEDSgo0bN1Y0SDOzquuPYwTAc8DxRdsT0n3FrgHuAYiIR4BBwJjuDhYRN0XEzIiYOXbs2AzCNTOrdf1vjOAxYJqkKZIGkAwGz+5SZjXwWgBJp5AkAn/dNzM7QD9sEUREG/AR4EHgGZJfBy2WdKOkS9NinwDeL+lJ4C7g3REZtn/MzPqzjO4jyHKwmIiYA8zpsu+GotdLgFdkGYOZ2RGhn44RmJlZWfW/MQIzMysbtwjMzKw/zjVkZmZl4jECMzPzGIGZWV1zi8DMzDxGYGZWxzxGYGZmHiMwM6trbhGYmZnHCMzM6lg1xwgk5SXNzSwCMzMrUZVaBBHRDhQkjcgkAjMzK0H11yzeAfxR0i+BnR07I+JjmURlZmb7i6j6egT3pQ8zMzvClJQIIuL76XKTJ5O0T5ZGxN5MIzMzsyJBVmMEJSUCSZcA3wGWp5FMkfSBiHggk6jMzOxAVe4a+hrw6ohYlsSiE4D7AScCM7NKyHA191LvI9jekQRSK4DtGcRjZmY9qm6LYIGkOcA9JHnpLcBjki4HiAgPJJuZZanQmtmhS00Eg4D1wPnp9kZgMPAGksTgRGBmlpVCAVp3wYAhmRy+10Qg6csR8Y/AnIi4N5MIzMysNPkBmRy2rzGCSyQJ+EwmtZuZWQkyHCmm766hnwMtwFBJ24r2C4iIGJ5ZZGZmtr9qzD4aEZ+KiJHA/RExvOgxzEnAzKxCMpx5FEr8+WhEXJZpFGZmVgKvR2BmVqdqoEVgZmY1IJsGgROBmVm9KzkRSPp68bOZmVVILQwWp85Ln8/vtZSZmWXEg8VmZnWqdloEZmZWTdW4oczMzGpADY0RmJlZVVW/RXBn+nxHFoGYmVl1lJwIIuJfi5/NzKxS3DVkZmbgwWIzs7rlwWIzM0tUYfF6SdeXcIydEfGdMsVjZmYHqG6L4FPAUGBYL49PZBmgmZmlMhoj6Gupytsj4sbeCkgaUsZ4zMyswvpaqvIf+jpAKWXMzOwwZDxYfFhjBBHxtfKGY2ZmPatO19CwTGo1M7ODUMUWQUR8IdPazcysdNW8oUzSiZIekvR0un26pM9lEpGZme2vRm4o+y7wGaAVICKeAq7MKigzM+tOdaeYaIqIP3TZ11buYMzMrDu10SJ4QdIJpNFIejPwfF8fknSRpKWSlkn6dA9l3ippiaTFku7sroyZmVG1G8o6fBi4CThZ0nPASuAdvX1AUh74NvA6YC3wmKTZEbGkqMw0ki6nV0REi6SjD+EczMzsMJSUCCJiBXBBehdxLiK2l/Cxc4Bl6WeRdDdwGbCkqMz7gW9HREtaz4aDCd7MrC7UyGAxABGxE7irxOLjgTVF22vTfcVOBE6U9DtJ8yVddDDxmJnVl+p2DRXrejE/3PqnAa8CJgC/kXRaRGzpWlDStcC1ABMnTixjCGZmta6GWgSpRSWWew44vmh7Qrqv2FpgdkS0RsRK4E8kieEAEXFTRMyMiJljx4492JjNzPq/aq9QJmmwpJMi4r0lfuQxYJqkKZIGkNx3MLtLmZ+StAaQNIakq2hFqTGZmdWFWhgjkPQG4Ang5+n2mZK6XtT3ExFtwEeAB4FngHsiYrGkGyVdmhZ7ENgkaQkwF/hURGw6tFMxMzvSVXeM4PMkvwKaBxART0ia0teHImIOMKfLvhuKXgdwffowM7MqKLVrqDUitnbZl21bxczMKqLUFsFiSW8H8ulNYB8Dfp9dWGZmdoAqDxZ/FHgpsIfkPoJtwHWZRGRmZvt7Me2Q2X3AL+vLotQ7i3cBn00fZmZWSe2tyfOoyZkcvqREIGku3YwJRMRryh6RmZntLwrJc/5Q7gHuW6lH/WTR60HAFXgaajOzyoj25Fn5TA5fatfQwi67fiep6/oEZmaWhY4WQa6KiUDSUUWbOeBsYEQmEZmZ2f4KHS2CQ5kVqG+ldg0tJBkjEEmX0ErgmkwiMjOz/bXtSZ6r1TUkKQdcHRG/yyQCMzPr3a505p1CayaH77OdEREF4FuZ1G5mZn3LNybPw8ZlcvhSO5weknSFlNFtbWZm1rOOweIq31n8AeBeYI+kbZK2S9qWSURmZra/jkRQzdlHI2JYJrWbmVnfOtYjyOhXQ6WuR/BQKfvMzCwDGXcN9doikDQIaALGSBrFvnbJcMq7drGZmfUo2xZBX11DHyCZZfQ4knsJOhLBNvxLIjOzyqjmGEFEfAP4hqSPRsS/ZxKBmZn1rhbGCJwEzMyqqEZ+PmpmZlVTAy0CMzOros6VydwiMDOrTx0rlDUMzOTwff18dCXdrExWXCR9/+sR8c1yBmZmZqmOhWkGZTP7f1+/GpqSSa1mZla6jhZBx+RzZeauITOzWrfpz8lzzonAzKw+DUine3OLwMysTkUhGR+oxn0Ekv5vJrWamVnpoj2zZSqh7xbBRZnVbGZmpSm0Qy67RNDXpHP5LrOO7iciNpc/JDMz20/GLYK+EsHJ7D/raLEAppY9IjMz21+hUNUWwZKImJFZ7WZm1re926GxKbPD+1dDZma1bu8uGJjdisF9JYJvZFazmZmVZncL5AdkdvheE0FE3CbpXZIel7QzfSyQ9M7MIjIzs/1tWQ2tuzI7fF+Tzr2LZKnK64HHSQaNzwK+Kiki4vbMIjMzs8TgkZlNOAd9dw19CHhTRMyNiK0RsSUiHgauAD6cWVRmZra/ptGZHbqvRDA8IlZ13ZnuG55FQGZm1kVEZquTQd+JYPchvmdmZuUShczmGYK+7yM4RdJT3ewXvpnMzKxCgqyWqYQSEkFmNZuZWWkiqtciiIjmrvskjQE2RURvS1iamVnZZNsi6Gsa6lmS5km6T9IMSU8DTwPrJXlmUjOzSsh4sLivrqFvAf8HGAE8DFwcEfMlnQzcBfw8s8jMzCyRcddQXymmISJ+ERH3An+JiPlJTPFsZhGZmVkXVewaAgpFr7v+XNRjBGZmlVDNwWLgDEnbSFLR4PQ16fagzKIyM7MiVfz5aERktxKCmZmVptBe1TECMzOrth1/gfbWzA7vRGBmVvPUf1sEki6StFTSMkmf7qXcFZJC0sws4zEz65dyeRh+XHaHz+rAkvLAt4GLgenAVZKmd1NuGPD3wKNZxWJm1q9FgWr+fBQASZdL+rOkrZK2Sdpe9AuinpwDLIuIFRGxF7gbuKybcl8Evgy8eFCRm5nViypPQ93hK8ClETEiIoZHxLCI6Gs9gvHAmqLttem+TpLOAo6PiPtLjtjMrJ5EALWRCNZHxDPlrFhSDvga8IkSy1+brpe8YOPGjeUMxcysdnXM71mtuYYkXZ6+XCDpR8BPgT0d70fEfb18/Dng+KLtCem+DsOAU4F5SkbDjwVmS7o0IhZ0PVhE3ATcBDBz5kzf1Wxm9SHSCR6qeGfxG4pe7wIuLNoOoLdE8BgwTdIUkgRwJfD2zg9HbAXGdGxLmgd8srskYGZWvzpaBNW7s/g9h3rgiGiT9BHgQSAP3BIRiyXdCCyIiNmHemwzs7rR2SKo3jTUSf3SVOAbwCyS9PQIcF1ErOztcxExB5jTZd8NPZR9VSmxmJnVlRe3Js+7WzKrotQUcydwDzAOOA64l+TnoGZmlqW29Jf1Y07KrIpSE0FTRNweEW3p44d49lEzs+x1zDGUb8ysipK6hoAH0iki7ibpGnobMEfSUQARsTmj+MzM6ltHIsiVerk+eKUe+a3p8we67L+SJDFMLVtEZma2z46/APDQ06sZObyFsyeNKnsVJSWCiJhS9prNzKxPSzfs5CTge39sZ9GS+dzxvlllTwYH/XskSTeVNQIzM+vRs+u2ANAWOVrbCsxfsansdRzKD1M9VbSZWYWccuxQACTR2JBj1tTRZa/jUEYfNpQ9CjMz69aJY5sAeOvLp/CPZ5a/WwgOIRFExEVlj8LMzLpXaAfgzTMnwvjyJwEofT2CCZL+W9JGSRsk/UTShEwiMjOzfbY913eZw1TqGMGtwGz23Vn8v+k+MzPLUsf9A01jei93OFWUWG5sRNxadGfxbcDYzKIyM7NEJF1D5AdkVkWpiWCTpKsl5dPH1UD5f8NkZmb7S8cIyOUzq6LURPBekruL/wI8D7wZOOQpqs3MrEQd01BXe4qJiGgGLs0sCjMz696mZclztdcslvR9SSOLtkdJuiWzqMzMLNGyKnkeOCyzKkpNMadHxJaOjYhoAWZkE5KZmXVqGASDRtTEGEFOUuedDOn009l1WJmZWaLQBiMmZlpFqRfzfwMekXRvuv0W4J+zCcnMzDqtWwTDjs20ilIHi38gaQHwmnTX5RGxJLuwzMwMgB3robEp0yoOZhj6KGBnRHwL2CjJaxSYmWUt1wjjz862ilIKSfon4B+Bz6S7GoEfZhWUmZkB7W3QvgdGZDu1W6ktgjeR3EewEyAi1gHZ/ZbJzMxgz7bkudCWaTWlJoK9EREk6xMjaUh2IZmZGbBv4fpRkzOtptREcI+k7wAjJb0f+BXw3ezCMjMzCmkiyHB6CSjhV0OSBPwIOBnYBpwE3BARv8w0MjOzerdpefLcMd9QRvpMBBERkuZExGmAL/5mZpXSvjd5HnNiptWU2jX0uKSXZRqJmZntr2OQuHFwptWU2vH0cuAdkppJfjkkksbC6ZlFZmZW7zoSQb4x02pKTQSvzzQKMzM7UMfMoxlOQQ0Htx6BmZlVUkcCGHpMptX0mmYkPd7XAUopY2Zmh2B9OqVbhmsRQN8tglMkPdXL+wJGlDEeMzPrMDhdD6xhYKbV9JUITi7hGO3lCMTMzLqIAgwcnnk1fSWC9cAHgZcAfwRujohsJ70wM7NEFEDKvJq+hqK/D8wkSQIXkyxQY2ZmlRCFzH8xBH23CKandxQj6WbgD5lHZGZmiQolgr5qaO144S4hM7MKq5EWwRmS0gmxETA43e64szj7UQwzs3pVC4kgIvKZR2BmZt3b8CxEZF5N9qnGzMwOTdNo2L0582qcCMzMalW0w9HTM6/GicDMrFYV2iGXfQ+9E4GZWa2KdpATgZlZ3dq2aw/Pb29lYXNLpvU4EZiZ1aCFzS0U1j3Jc1t2847vzc80GTgRmJnVoPkrNrExRjCQvbS2FZi/YlNmdTkRmJnVoFlTR9OodlbEcTQ25Jg1dXRmdZW6VKWZmVXQ2ZNGgf7C7uNmcscls5LtjDgRmJnVqlwDp4wMyDAJQMZdQ5IukrRU0jJJn+7m/eslLZH0lKSHJE3KMh4zs/5FMPbEzGvJLBFIygPfJlnHYDpwlaSut8gtAmZGxOnAj4GvZBWPmVm/cwTcR3AOsCwiVkTEXuBu4LLiAhExNyJ2pZvzgQkZxmNm1n9EJLOP9vM7i8cDa4q216b7enIN8ECG8ZiZ9R9RSJ4r0CKoicFiSVeTLIl5fi9lrgWuBZg4cWKFIjMzq5I925PnQmvv5cogyxbBc8DxRdsT0n37kXQB8Fng0ojY09PBIuKmiJgZETPHjh1b9mDNzGrK5uXJ89BjMq8qy0TwGDBN0hRJA4ArgdnFBSTNAL5DkgQ2ZBiLmVn/sn5J8jxmWuZVZZYI0jWOPwI8CDwD3BMRiyXdKOnStNhXgaHAvZKekDS7h8OZmdWXF/6UPE84J/OqMh0jiIg5wJwu+24oen1BlvWbmfVbHSuTNQ7OvCrPNWRmVova9sCgkSBlXpUTgZlZDdqz8hF2t0XmaxGAE4GZWc1Z2NxC2/YXWL93UOZrEYATgZlZzZm/YhOt5FkYJ2a+FgE4EZiZ1ZxZU0fTSBstMSzztQigRu4sNjOzfc6eOBK0hxlTxnLHBdmuRQBOBGZmtWdX8tPRs8cPyXwtAnDXkJlZ7dm6OnkeWZl51ZwIzMxqTcuq5HnsyRWpzonAzKzGrPvzEwAs3t5UkfqcCMzMasjC5hZ+9PjzAFz5kw2+oczMrN7MX7GJXLoGwa42ZX4PATgRmJnVlFlTR3NqvhmAxoZ85vcQgH8+amZWU86eNIpdR+2EFrjjfdnfQwBuEZiZ1Zx86w62DK7ckrxOBGZmNWTRn1czcMdaVu4cUJEJ58CJwMyspjyz9FkA7m5/dUUmnAMnAjOzmnLm+CEAbGVoRSacAw8Wm5nVlMbdGwGYdeJxvP9VHiw2M6srC5tb+MXP/weAO5Y1VqxeJwIzsxoxf8Umhhe2AbCqdXRFxgfAXUNmZjVjVNMAXpZ7hl0xkFbyjGoaUJF63SIwM6sRLbv2Mi33HM1xDEq3K8GJwMysRkxgAwCLCi8hwC0CM7N6M+KZOwD4beE0cnKLwMysrixsbmHc+l8D8GDhZTTkVJF7CMCJwMysJvzk8bUMYTcvRiMFcrzqpKMrcg8BOBGYmdWEwubVTNAL3NX+GgDGDhtYsbqdCMzMqmxhcws7Vz4KwOKYTENeXH7WhIrV70RgZlZlP3l8Ldfl7gHg9+0v5TUV7BYCJwIzs6pb+ZfNnJBL1ilex5iK1+9EYGZWRQubWxiz9pcAfLftEqCy4wPgRGBmVlU/eXwtb8s9DMBtba8nJyo6PgBOBGZmVbV9+WO8Mr8YgHWMZuakURUdHwAnAjOzqvna/y7g33d8HIAL93yZIMe0Y4ZVPA4nAjOzKljY3MLYR78EwK1tr+dPcTyqQrcQeBpqM7Oq+Nx9T/E/+bkAfKnt7QB84K+nVrxbCNwiMDOruOvuXsSbN/0nA9TOz9pfzl4aGTt0AJ++5JSqxONEYGZWQXc+uponnlzINQ0PAPDZ1msAmDGx8i2BDk4EZmYVsrC5hVvvn8e8gZ8A4MbWv2UrQxHwgfNPqFpcHiMwM6uAj9+1kGOevon7G+4F4Pa2C7il/WIA/vlNp1VlbKCDE4GZWYYWNrfwiR8tYt6uN0Fjsu+KPf/EwjgJgA+eN5W3v3xiFSN0IjAzy8TC5ha+8+vlbH7m19w14FugZP/0F29hF4MAeOOZx1VtgLiYE4GZWRksbG7hvsfXsnH7HrbufJGj1vySjzf8hFMGrgFgQeFErt77GV5kIGOGDuD6151U9ZZABycCM7PD0PHN/1dLnuet+Xn8XX4uZ+aWQ7ru/NZo4rrWDzO3MAOA86aN4QfXvLyKER/IicDM7CB0/eY/cM1vOT/3JDcNeqCzTEsM5Za2i7i1/SJ20NS5/4PnTa2JrqCunAjMzPpw56OrueV3K9m9azt/s/tnvC//MOO0mUFq7fzmv7QwgfmFU/iXtqvYnY4BQDI0cMH0Y/jg+SdU9ZdBvXEiMDMr0tHVs2LjDkY3NTB8VzNTNv+Wm/MPMSm3ofOXP2sKY3mqMIVnCpP4eeFlLItkjiAB50wexcimAYwZNpArzppQswmggxOBmdWF4i4dgC279rKnrcCUMUNYuXEHx+U2cULDJtauWc55PMuHcys5Y/uK5MON0Bp5HmqfwUOFs3iofQbrOarz2DklF/+XHDOsX1z4u3IiMLN+p+OiHsAV6Wyd81dsYlTTABav29p5sR/f1MoZY/I8tvx5li3/MxO1ntO0gXHazDStZZLW07ixjSHas+/g6VWxLXI8WjiZPxamcH/7LJ6KqbST7ywm4GX9+OJfzInAzDK3sLmF+Ss2MWvqaM6eNKr07SmjyO/dzqIVzzNzwhDyrdv5zeLVLH72GY6OTUzR86xc1MZJrGYGTTTSymu1gVYamKAXOut/I3T25RdbUxjL8jiOp2IKBXI0F45hTYzlT3E82xhyQPkxQwcwdcyQI+LiXyzTRCDpIuAbQB74XkT8S5f3BwI/AM4GNgFvi4hVWcZkVsu6XhArcYxSyi9sbmH+8hc4d9JQzho/BAptyaO9FdpehPZWljy3mT+u2czp44dxytFNEO1QaGd582pu+9WfGFLYzs78LhpOGsWKpYsZEXmWzy0wfuxOdrywh1OjQOu8dnYO3swZu9dzKnkG/roVgDOLYpkOB1y5VhaOYSQ7WB1HsyQmMVy7mNt+Jg208+eYQEsMJa8CzYVjWBnjeIHhdN7h1UVDXpx1/Eg279zLUUMG9Ku+/kOliMjmwFIe+BPwOmAt8BhwVUQsKSrzd8DpEfFBSVcCb4qIt/V17JkzZ8aCBQsOOqY7H13NA08/z8WnjsvsRo5y/Efuz2rl/MsSR9ve5GJHQBTSR+z/TPDkmhYWNW/irIkjOf24Ycn+1t3Jc6E1vWC2QevOdH/Ai1th92bID4BCO0Q7azdt597HVhGFdiblXuAVp5/EscMai+osdNa533YUkgvyljW0MIwnmjeSj1YmaQNjxhzNkAG5NJZC58WZaIco0LprG3pxCzsYDIihA5VcY9P3iQJRaCcKBXIq/7ViQ4xEDQMZ3LaVZTGevTTSNHgw7bu3siFGsjgm08QemuMY9tLA3mgkgHUxmh00sTqOZm9+MDmgrRAU+ggxJ5g5KRnI3W+M4IWdDGzIHXHf9ItJWhgRM7t7L8sWwTnAsohYkQZxN3AZsKSozGXA59PXPwa+JUmRQXa689HV3P8/d5KnwIPLYdzGybz6pKOBtKr9qjy0fcs27uDWXy6lrRA8mxMjLziRE8YOOcw6KLFcL/s69wdsfx4Gj9pXrmP/IT/Tub1uyy4WP/4sQwrB43MHMP70cRw7LG2Pd1w8Oz7beUHr6TX7Xre3phemogtg5zELXfYH23bvZe9zLZwVwch5W9h11HCaBg1KL4KF5MIc7clzyypoGp0eo33fxbJ114F/Bz04I30wv+SPdGsC8PEc++YEfnoeNA4B5dIHRa+7PADa9jCgLRgTw2ilgfUxihG7tsBR0yGX31c2lwcl28s37KR53Xp2MJidDOaMsUdxxvGjQer8zMLVW/n9ihYKiKO0nROmvZRXnDgO8g2Qa4BcI79ZvIq7n9lLW4iC8lw2YyJvOHMC5BpYumEnN9y/jPVtQ9neMIJPXjKDL8xZSmtbgcaGHDdc+FJu/Nnifduv2bedz+cggvZCJK/Z9/pVJ47lrPSbOnQ/RjBm2EBOPW4Ei9dt7RxPOBIv8ocry0QwHlhTtL0W6Ho7XWeZiGiTtBUYDbzQpRySrgWuBZg48eC/zT/w9PN8t/FrNHUMCi1IH2X0EuBbeegcT5pb3uPXuuOAd6YXstbIw5IGyKcXIJRcXDqeO1/n+n6db+i8cHXuP6Dsvv07d+ylIdooIDbGCEa/uJOm0RMg15hc3DouhLkGmHgu7FgPR01N9+X3XSx3bYLRLym62BbXldQ3788v8MslG0lSU47XvXQcrznlWNi7E5qOgoZBkO+otzFJMk1jkpZAwwAYOAJyOVCex5/bwbtuW8CeNsg1NHDH+8496IvWs80tvON78zsvqndcOavXY+xsbuHvi8tfOAu6lFdzC/+xsqjM+QeWGXJUCw8vTcvkc3zoZfvKnDQV/mHc/i20E48btd/2SccO63Eb6PZ11/PyBf7QZdk19Gbgooh4X7r9t8DLI+IjRWWeTsusTbeXp2UOSATFDqVr6M5HV3PvT+9D6bfkD5z/El7/0mM7AumIqOgEDnjRQ7l9+5Y8v43/899/pLU9aMjn+NLlpzN93PCSPnvo+ziIcgICGgbvf1HufKaH/aU9L1y9hXfc/CitbZFcMN7X+0UoKwu7XgwzjKPcddX0GEEZylj19NY1lGUiOBf4fES8Pt3+DEBEfKmozINpmUckNQB/Acb21fKaRvcAAAleSURBVDXkMYLaVSvnX8k4auWczXpTrUTQQDJY/FrgOZLB4rdHxOKiMh8GTisaLL48It7a17EPNRGYmdWrqgwWp33+HwEeJOk1vyUiFku6EVgQEbOBm4HbJS0DNgNXZhWPmZl1L9P7CCJiDjCny74bil6/CLwlyxjMzKx3XrzezKzOORGYmdU5JwIzszrnRGBmVuecCMzM6pwTgZlZnXMiMDOrc04EZmZ1zonAzKzOZTbXUJYkbQSaD/HjY+hmmusjnM/5yFdv5ws+54M1KSLGdvdGv0wEh0PSgp4mXjpS+ZyPfPV2vuBzLid3DZmZ1TknAjOzOlePieCmagdQBT7nI1+9nS/4nMum7sYIzMxsf/XYIjAzsyJHbCKQdJGkpZKWSfp0N+8PlPSj9P1HJU2ufJTlU8L5Xi9piaSnJD0kaVI14iynvs65qNwVkkJSv/+FSSnnLOmt6d/1Ykl3VjrGcivh3/ZESXMlLUr/fV9SjTjLRdItkjZIerqH9yXpm+mfx1OSzjrsSiPiiHuQLI25HJgKDACeBKZ3KfN3wH+lr68EflTtuDM+31cDTenrD/Xn8y31nNNyw4DfAPOBmdWOuwJ/z9OARcCodPvoasddgXO+CfhQ+no6sKracR/mOZ8HnAU83cP7lwAPAAJmAY8ebp1HaovgHGBZRKyIiL3A3cBlXcpcBnw/ff1j4LWSVMEYy6nP842IuRGxK92cD0yocIzlVsrfMcAXgS8DL1YyuIyUcs7vB74dES0AEbGhwjGWWynnHMDw9PUIYF0F4yu7iPgNyRruPbkM+EEk5gMjJY07nDqP1EQwHlhTtL023ddtmYhoA7YCoysSXfmVcr7FriH5RtGf9XnOaZP5+Ii4v5KBZaiUv+cTgRMl/U7SfEkXVSy6bJRyzp8Hrpa0lmSN9I9WJrSqOdj/733KdPF6qz2SrgZmAudXO5YsScoBXwPeXeVQKq2BpHvoVSStvt9IOi0itlQ1qmxdBdwWEf8m6VzgdkmnRkSh2oH1F0dqi+A54Pii7Qnpvm7LSGogaVJuqkh05VfK+SLpAuCzwKURsadCsWWlr3MeBpwKzJO0iqQvdXY/HzAu5e95LTA7IlojYiXwJ5LE0F+Vcs7XAPcARMQjwCCSOXmOVCX9fz8YR2oieAyYJmmKpAEkg8Gzu5SZDbwrff1m4OFIR2L6oT7PV9IM4DskSaC/9xtDH+ccEVsjYkxETI6IySTjIpdGxILqhFsWpfy7/ilJawBJY0i6ilZUMsgyK+WcVwOvBZB0Ckki2FjRKCtrNvDO9NdDs4CtEfH84RzwiOwaiog2SR8BHiT51cEtEbFY0o3AgoiYDdxM0oRcRjIwc2X1Ij48JZ7vV4GhwL3pmPjqiLi0akEfphLP+YhS4jk/CFwoaQnQDnwqIvprS7fUc/4E8F1JHycZOH53P/5Sh6S7SJL5mHTc45+ARoCI+C+ScZBLgGXALuA9h11nP/7zMjOzMjhSu4bMzKxETgRmZnXOicDMrM45EZiZ1TknAjOzOudEYGZW55wIrCSS2iU9IelpSfdKapI0WdJuSU9UOz4ASbdJenMJ5TrO5UlJj0v6q3R/x/kskvSMpD9Ienf63nvSzzwhaa+kP6av/6Wb46/qod7rJDUd3lkePEmfl/TJbvYPLjqfw7oTV9IbJU0/nGNY9RyRN5RZJnZHxJkAku4APgjcByzv2F8snclVNTrfS/G5vB74EvvmXloeETPS96YC90lSRNwK3JruXwW8OiJeOMh6rwN+SHIT0H4k5SOi/VBO5lBFxG7gzJ4S10F6I/AzYEmpH5DUkE74aFXmFoEdit8CL+m6M/1GvVTSD4CngeMlfUrSY+kCGl8oKvvOdN+Tkm4v+vzD2rd4zsTuKpe0StJX0m/lf5BUHMt5kn4vaUUprQOS6YtbunsjIlYA1wMfK+E4vZL0MeA4YK6kuem+HZL+TdKTwLnpeY1J35spaV76eoiSxUr+kLZWuptuG0nzJH2jqOV2TtHb09P3V6SxlE3aoroU+Gpa9wk9xZK2Tm6X9Dvg9nLGYYfOLQI7KEom6LsY+HkPRaYB74qI+ZIuTLfPIVlEY7ak80gm9/sc8FcR8YKko9LP/jvw/Yj4vqT3At8k+abZna0RcZqkdwJfB/4m3T8OeCVwMsmcLD/u5rOD0+6sQWn51/Ryyo+nxzosEfFNSdezf0tiCMmiIp8AUM/LYXyWZC6s90oaCfxB0q8iYmc3ZZsi4sz0z/kWkon3SM/h1SST8S2V9J8R0dpbzJJ+m5bv6pMR8auic/u9pNnAzyLix0Xn0lMs04FXpi0SqwFOBFaqjosnJC2Cm0m+4XbVnC6WAXBh+liUbg8lSQxnAPd2XBAjomMRjnOBy9PXtwNf6SWeu4qe/1/R/p+m3VFLJB3Tw2eLu4bOBX4g6dQeyma5WFE78JMSyl0IXFrUzz8ImAg8003ZuyBZ3ETS8DRxANyfzji7R9IG4BiSmUp7FBF/XUJsvekpltlOArXFicBK1Xnx7NDDN9jib6kCvhQR3+nyuYNaOETSgyQXrgUR8b50d/EkWcWvi6fX7vMiHhGPpN0xY3soMoPuL7jl8GKXcYE29nXXDiraL+CKiFha/GFJt6bxrYuIjnV6u04e1rFd/OfSTgn/90ttEfSip1i6a8lYFXmMwLL0IPBeSUMBJI2XdDTwMPAWSaPT/R1dQ79n3yyw7yBpeRARr4+IM4uSAMDbip4fOdQAJZ1MMqvlATN0SpoM/CtJl1U5bKf7C2uHVcDZ6esrivY/CHw0HYDvmFKciHhP+udSvFj729IyryTpPtt6qMFGxF+nx+/66C4JdHduZYvFsuUWgWUmIn6hZH74R9Jr2A7g6nQa4X8Gfi2pnaTr6N0kSwzeKulTJPPJ9za97ihJT5F8073qIEMr7uYSyZhGexrjCZIWkXwj3w58MyJuO8jj9+Qm4OeS1kXEq7t5/wvAzZK+CMwr2v9FknGQp5SsvLaSfWMiXb2Yxt8IvLdMcZfibpKpoD9Gsr5HNWOxg+RpqO2Qpd+YfxYRPfWvZ1XvKmDmIfx8syIkrUoXw6l0vfNIum0OavGdLP48DzUWqw53DdnhaAdGqEZuKLODo/SGMpJv7LV4v4dViFsEZmUm6bqI+Hq14zArlROBmVmdc9eQmVmdcyIwM6tzTgRmZnXOicDMrM45EZiZ1bn/D7MIcUNr4X+XAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# the blue points are computed by us above\n",
    "# the orange is from sklearn\n",
    "\n",
    "# true positive rate = tpr = P[reco-ph BDT | true-ph] = probability to pass photon selection given true-photon\n",
    "# false positive rate = fpr = P[reco-ph BDT | true-el] = probability to pass photon selection given true-electron\n",
    "\n",
    "# the photon selection is as described above BDT > threshold\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.plot(effs_ph, effs_el, '.')\n",
    "ax.plot(tpr, fpr)\n",
    "ax.set_xlabel('P[reco-ph BDT | true-ph] = tpr')\n",
    "ax.set_ylabel('P[reco-ph BDT | true-el] = fpr')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.014, 0.025)"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAFlCAYAAAAarsR5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAe90lEQVR4nO3df5Bd5X3f8fdX6IdxIoyyyAOWrBUK0FREtZPd4o2LGw/UFOOx5cikBjsxIih00tCZ2k1bMp44KeGPOOOMk9RMMoyw4mAV7Ch1UMdumDRAE4jXwy42lgUhFho2SJYckNfetZHFrvbbP+5ZdHXZ1bnae+/eu/e+XzM7uvfcc4+ecxD66Ps853lOZCaSJJ3JsnY3QJLU+QwLSVIpw0KSVMqwkCSVMiwkSaUMC0lSqeXtbsDZuOCCC3Ljxo3tboYkLRkXXHABDz744IOZeW0jx1lSYbFx40ZGRkba3QxJWlIi4oJGj2E3lCSplGEhSSplWEiSShkWkqRShoUkqZRhIUkqZVhIkkoZFpKkUoaFJKmUYSFJKmVYSJJKGRaSpFKGhSSplGEhSSplWEiSShkWkqRShoUkqZRhIUkqZVhIkkoZFpKkUoaFJKmUYSFJKmVYSJJKGRaSpFJ1hUVEXBsRz0TEgYi4fY7PV0XE54rPvxIRG4vt74iI0YjYV/x61Rzf3RsR32j0RCRJrVMaFhFxDnAX8E5gM3BjRGyu2e0WYDwzLwE+CXy82P4i8O7M3ALcBNxbc+xtwPcbOgNJUsvVU1lcARzIzIOZ+TJwP7C1Zp+twGeK13uAqyMiMvOrmfmtYvt+4NyIWAUQET8KfAS4s9GTkCS1Vj1hsQ54vur9oWLbnPtk5jTwPaCvZp/3AU9k5oni/W8Dvwe8dKbfPCJujYiRiBh54YUX6miuJKnZFmWAOyIup9I19e+L928Gfjwzv1D23cy8OzMHM3Nw7dq1LW6pJGku9YTFYeCNVe/XF9vm3CcilgOvA44V79cDXwA+lJnPFvv/DDAYEc8BjwKXRcQjCzsFSVKr1RMWjwOXRsTFEbESuAHYW7PPXioD2ADXAw9lZkbE+cAXgdsz87HZnTPzjzLzDZm5EbgS+IfMfHtjpyJJapXSsCjGIG4DHgSeBj6fmfsj4o6IeE+x2z1AX0QcoDJoPXt77W3AJcDHIuJrxc/rm34WkqSWisxsdxvqNjg4mCMjI+1uhiQtKRExmpmDjRzDGdySpFKGhSSplGEhSSplWEiSShkWkqRShoUkqZRhIUkqZVhIkkoZFpKkUoaFJKmUYSFJKmVYSFKXGB0b566HDzA6Nt70Yy9v+hElSYtudGycB+65k+t4jJOPwORF57F61YqmHd/KQpK6wPDBY1zHY2yOMWYSJo5PN/X4VhaS1AWGNvVx8hF4KvvZnr/J7m1DrOtfU/nwl6Lh4xsWktQFBl74C4inOXz+ALu3DTEwGxRNYlhIUjfYtweAdW/7xVMVRRMZFpLUqUZ2vRICpY7ug/4rYfDmljTFAW5J6lT79lRCoB4XboEt17esKVYWktTJLtwCN3+x3a2wspAklTMsJEml7IaSpE5SPah9dF+lG6oDWFlIUiepHtRu8aD12bCykKRO0yGD2tWsLCSpU4zsgrFH292KOVlZSFKbjY6NM3zwGB965j5WQ8d0PVUzLCSpjaqXFs8YY/Kit7C6RbOwG2E3lCS1UfXS4k9lPyOr/027mzQnw0KS2mhoUx/L4tTS4udd+cvtbtKc7IaSpDYa6F/D5EXnMXF8uiVLizeLYSFJbbZ61QpWr1rRkqXFm8VuKElSKcNCklTKsJAklTIsJKmdOnjWdjXDQpLaaXaF2Q6ctV3Nu6EkqZXKnqPd4mdnN4uVhSS1UtlztDtoGfIzsbKQpFbrwCXHz5aVhSSplGEhSSplWEiSShkWkqRShoUkNcno2Dh3PXyA0bHxdjel6bwbSpLO1hxzJyZPTHHyyAQDCScfgcmLzmP1qhWV22Yv3NKedjaRlYUkna055k5MHJ9mJiuvZ7LyHlgy8yjKWFlI0kLUzJ04OjbO9p3DTE3PsGL5MnZvG+ro51OcLcNCks7G7MJ//Veetnmgfw27dwwxfPAYQ5v6OvaJdwtlWEjS2TjDwn8D/Wu6LiRm1TVmERHXRsQzEXEgIm6f4/NVEfG54vOvRMTGYvs7ImI0IvYVv15VbH9tRHwxIv4+IvZHxO8086QkqaWWwMJ/zVYaFhFxDnAX8E5gM3BjRGyu2e0WYDwzLwE+CXy82P4i8O7M3ALcBNxb9Z1PZOZPAD8F/KuIeGdDZyJJapl6KosrgAOZeTAzXwbuB7bW7LMV+Ezxeg9wdUREZn41M79VbN8PnBsRqzLzpcx8GKA45hPA+kZPRpLUGvWExTrg+ar3h4ptc+6TmdPA94C+mn3eBzyRmSeqN0bE+cC7gb+uv9mSpMW0KAPcEXE5la6pa2q2LwfuA/4wMw/O891bgVsBNmzY0OKWSpLmUk9lcRh4Y9X79cW2OfcpAuB1wLHi/XrgC8CHMvPZmu/dDXwzM39/vt88M+/OzMHMHFy7dm0dzZUkNVs9YfE4cGlEXBwRK4EbgL01++ylMoANcD3wUGZm0cX0ReD2zHys+gsRcSeVUPlPjZyAJKn1SsOiGIO4DXgQeBr4fGbuj4g7IuI9xW73AH0RcQD4CDB7e+1twCXAxyLia8XP64tq46NU7q56oti+o7mnJklqlrrGLDLzS8CXarZ9rOr1D4Gfn+N7dwJ3znPYqL+ZkqR2ciFBSVIpl/uQpDLVS5J3yZLjZ8vKQpLKVC9J3iVLjp8tKwtJKjF5YoqJ11zC0as+27ULBZYxLCSpWs1T8CZPTJFH9vF89rN95zC7dwz1ZGDYDSVJ1WqegjdxfJqnsp8HTr6VqekZhg8ea2Pj2sfKQpJmVT/YqHgK3itPwJupPAFvaFPtsne9wbCQpFlzPNio25+AVy/DQpKqzfFgo25+Al69HLOQJJUyLCRJpQwLSRrZBbveddpdUDqdYSFJs7fL9ujs7Ho4wC1JUAmK4nZZvZqVhaTeNju3QmdkWEjqbXPMrdCr2Q0lqTfNrgF1dN+ccyt0OisLSb3JQe2zYmUhqXc5qF03KwtJUinDQpJUyrCQJJUyLCRJpQwLSVIpw0KSVMqwkNSTJk9Mcfi7xxkdG293U5YEw0JSzxkdG+fpIxM8P/4SH9w5bGDUwbCQ1HOGDx5jJiuvp6ZnGD54rL0NWgIMC0k9Z2hTH8ui8nrF8mUMbeprb4OWAJf7kNQbZhcOBAaA6VWH+Pa5l7F72xAD/Wva27YlwMpCUm+YXTiwsPyiN7Hubb9oUNTJykJS73DhwAWzspAklTIsJEml7IaS1L2qBrVfedCRFsTKQlJXGB0b566HD5w+wa56UNsn4jXEykLS0lVUDpMnpjh5ZIKBhJOPwORF57F61YpT1YSD2g2zspC0dBWVw8Tx6VdmZM8kTByfrryxmmgaKwtJS9uFWzh61WfZvnOYqekZVixfxu5tQ6xz/kRTGRaSlryB/jXs3jHE8MFjDG3qc6JdCxgWkrrCQP8aQ6KFHLOQtDSN7IKxR9vdip5hWEhammbnTziAvSgMC0lLV/+VMHhzu1vREwwLSVIpw0KSVMqwkCSVMiwkSaUMC0lSKcNCklSqrrCIiGsj4pmIOBARt8/x+aqI+Fzx+VciYmOx/R0RMRoR+4pfr6r6zkCx/UBE/GFERLNOSpLUXKVhERHnAHcB7wQ2AzdGxOaa3W4BxjPzEuCTwMeL7S8C787MLcBNwL1V3/kj4JeBS4ufaxs4D0m9YmQX7HrXqedUaFHUU1lcARzIzIOZ+TJwP7C1Zp+twGeK13uAqyMiMvOrmfmtYvt+4NyiCrkIOC8zhzMzgT8F3tvw2UjqfrMPNHL58UVVz0KC64Dnq94fAt4y3z6ZOR0R3wP6qFQWs94HPJGZJyJiXXGc6mOuO8u2S+pVPtBo0S3KqrMRcTmVrqlrFvDdW4FbATZs2NDklknqWNXPz67ms7Tbop5uqMPAG6very+2zblPRCwHXgccK96vB74AfCgzn63af33JMQHIzLszczAzB9euXVtHcyV1hernZ1ez+6kt6qksHgcujYiLqfyFfgPwgZp99lIZwP4ycD3wUGZmRJwPfBG4PTMfm905M49ExEREDAFfAT4E/I+Gz0bS0lVbSfj87I5SWllk5jRwG/Ag8DTw+czcHxF3RMR7it3uAfoi4gDwEWD29trbgEuAj0XE14qf1xef/QdgJ3AAeBb4P806KUlLUG0lYQXRUaJyM9LSMDg4mCMjI+1uhqQWmPzja5g4Ps3RbX/uE++aLCJGM3OwkWP4WFVJi2uOgevJE1PkkX08n/1s3znM7h1DBkaHcbkPSYtrjoHriePTPJX9PHDyrUxNzzB88FibGqf5WFlIao2yW1+rBq6Pjo2zfecwUzMzrFi+jKFNfYvYUNXDsJDUGtUzravNMXA90L+G3TuGGD54jKFNfXZBdSDDQlLrnMWtrwP9awyJDuaYhSSplGEhSSplWEiSShkWkqRShoWklpg8McXh7x5ndGy83U1RE3g3lKTG1cypcEZ297GykNS4mlnZzsjuPlYWkhozsgvGHoX+K1+ZU+GM7O5jWEhqzGz3U9WsbGdkdx/DQlLj+q+EwZtP2+SM7O5iWEia33yLAVbzmdg9wQFuSfOb7znY1XyiXU+wspB0Zj4HW1hZSJLqYFhIkkrZDSXpdNWD2g5eq2BlIel01YPaDl6rYGUh6dUc1FYNKwtJUinDQpJUyrCQJJVyzEISo2Pjpxb9a3dj1JEMC6mXjexicuQ+Th6ZYCDh5CMwveoQyy96U7tbpg5jN5TUy/btYeUL+5nJytuZhG+fe5m3y+pVrCykHvfy2svZfvjDTE1XHlS0e9sQ61xaXDUMC6lXFU+4W91/pQ8qUinDQupVVU+480FFKuOYhdTL5njCnTQXw0KSVMqwkCSVcsxC6iUuP64FsrKQeonLj2uBrCykXuPy41oAKwtJUinDQpJUym4oqZtVD2iDg9paMCsLqZtVD2iDg9paMCsLqVsVaz/Rf6UD2mqYlYXUrarWfpIaZVhI3cy1n9QkdkNJ3aB2IBsczFZTWVlI3aB2IBsczFZTWVlIS50D2VoEVhbSUudAthaBYSF1Awey1WKGhbRUjeyCXe969ViF1AJ1hUVEXBsRz0TEgYi4fY7PV0XE54rPvxIRG4vtfRHxcER8PyI+VfOdGyNiX0R8PSL+MiIuaMYJSb1gdGycw397L9NHnnQgW4uidIA7Is4B7gLeARwCHo+IvZn5VNVutwDjmXlJRNwAfBx4P/BD4DeAnyx+Zo+5HPgDYHNmvhgRvwvcBvxWU85K6lYju5gcuY+TRyZYzRhP0M85V32Wgf417W6Zulw9lcUVwIHMPJiZLwP3A1tr9tkKfKZ4vQe4OiIiM3+QmY9SCY1qUfz8SEQEcB7wrYWehNQz9u1h5Qv7mUl4Kvt5YPqtDB881u5WqQfUc+vsOuD5qveHgLfMt09mTkfE94A+4MW5DpiZUxHxK8A+4AfAN4FfnWvfiLgVuBVgw4YNdTRX6m4vr72c7Yc/zNT0DCuWL2P3pr52N0k9oC0D3BGxAvgV4KeANwBfB359rn0z8+7MHMzMwbVr1y5iK6UOU8ynWL1qBbt3DPGRa/4Zu3cM2QWlRVFPZXEYeGPV+/XFtrn2OVSMR7wOOFNt/GaAzHwWICI+D7xq4FxSlar5FAP9awwJLap6wuJx4NKIuJhKKNwAfKBmn73ATcCXgeuBhzIzz3DMw8DmiFibmS9QGTx/+mwbL3W96jWfju5zPoXapjQsijGI24AHgXOAT2fm/oi4AxjJzL3APcC9EXEA+A6VQAEgIp6jMoC9MiLeC1yTmU9FxH8H/iYipoAxYHtzT03qArNrPl24xVtk1VZx5gKgswwODubIyEi7myEtnl3vqvzqmk9qQESMZuZgI8dwBrfUiZydrQ7jqrNShxgdG2f44DGGNvUxUN39ZNeTOoBhIbVL1eD15IkpTh6ZYCDh5CMwveoQyy96k91P6hh2Q0ntUvXAoonj08wUw4czCd8+9zIrCnUUKwupnS7cAjd/kaNj42zfOXxqVva2IdY5j0IdxLCQOsBA/xp27xg6NWZhUKjDGBZSh3BWtjqZYxaSpFJWFlIrVC/TMZ/ZW2OlJcDKQmqFqjud5uUcCi0hVhZSM81WFLNVg/Mk1CWsLKRmcua1upSVhdRsVhTqQlYWUrMUT7KTupFhITVL1ZPspG5jWEjN5JPs1KUMC0lSKcNCklTKsJAa5VPt1AMMC6lRzq1QD3CehdQMzq1Ql7OykCSVMiwkSaUMC0lSKcNCqtPo2Dh3PXyA0bHxdjdFWnQOcKt31fOAosLkiSlOHplgIOHkIzB50XmsXrWi8qEPMVIPsLJQ76rnAUWFiePTzGTl9UxW3r/CW2bVA6ws1NvqvOX16Ng423cOMzU9w4rly9i9bYh1/WsWoYFSZzAs1Fuqu57OovtooH8Nu3cMMXzwGEOb+hgwKNRjDAv1lurZ1mfZfTTQv8aQUM8yLNR9zjRw7bOxpQVxgFvd50wD1w5GSwtiZaElZXRsvL5xA6sHqakMCy0NI7uYHLlv/rkO1Zz3IDWd3VBaGvbtYeUL++ef61DNriap6aws1LlqbnN9ee3lbD/8Yec6SG1gWKhz1dzmunrL9ex+p3MdpHYwLNTZagaqB8CQkNrAsFD7zTcvwoFqqWM4wK32m29ehAPVUsewstDicFa1tKRZWWhxOKtaWtKsLLR4rB6kJcuwUGvUdjs5WC0taXZDqTVqu53sapKWNCsLtY7dTlLXsLKQJJUyLNQSkyemOPzd44yOjbe7KZKawLBQ84zsgl3vYvKPryGP7OP58Zf44M5hA0PqAoaFmqcY1J44Ps1T2c8DJ9/K1PQMwwePtbtlkhpU1wB3RFwL/AFwDrAzM3+n5vNVwJ9SWeftGPD+zHwuIvqAPcC/BP4kM2+r+s5K4FPA24EZ4KOZ+ecNn5EWzzy3xx696rNs3znM1ExlKfGhTX3ta6OkpigNi4g4B7gLeAdwCHg8IvZm5lNVu90CjGfmJRFxA/Bx4P3AD4HfAH6y+Kn2UeCfMvOyiFgG/FjDZ6PFVb2EOLxye+xA/xp273Apcamb1FNZXAEcyMyDABFxP7AVqA6LrcBvFa/3AJ+KiMjMHwCPRsQlcxz3l4CfAMjMGeDFBZ2B2mue22MH+tcYElIXqScs1gHPV70/BLxlvn0yczoivgf0MU8ARMT5xcvfjoi3A88Ct2Xmt+fY91bgVoANGzbU0Vy1xFwLATorW+oZ7RrgXg6sB/4uM38a+DLwibl2zMy7M3MwMwfXrl27mG1UtbkWAnRWttQz6qksDgNvrHq/vtg21z6HImI58DoqA93zOQa8BPyv4v2fURn3UCdzRrbUs+qpLB4HLo2Ii4s7mG4A9tbssxe4qXh9PfBQZuZ8Byw++99U7oQCuJrTx0AkSR2ktLIoxiBuAx6kcuvspzNzf0TcAYxk5l7gHuDeiDgAfIdKoAAQEc8B5wErI+K9wDXFnVT/rfjO7wMvADc399QkSc1S1zyLzPwS8KWabR+rev1D4Ofn+e7GebaPAf+63oaqhc70FLtZDmZLPc0Z3D1udGycw397L9NHnjzzjg5mSz3NJcp71cguJkfu4+SRCVYzxhP0c85Vn3VuhKQ5WVn0qn17WPnCfmaSyjpO0291DSdJ87KyWMJGx8YbWlLj5bWXs/3wh5marqzhtNs1nCTNw7DoRHUMOE+emOLkkQkGEk4+ApMXncfqVSvq/z2O7mP1hVtcw0lSXQyLTlS7QN8cJo5PM1PMZJnJyvuzCouqRf8MCUllDItOVTJb+ujYeGUZ8NkupG1DrPMvfUktYlgsUS4DLmkxGRZLmF1IkhaLYdEpqge1nS0tqcM4z6LJRsfGuevhA4yOjZ/dF6uXAHe2tKQOY2XRRKNj4zxwz51cx2NnfzvrbDXhEuCSOpCVRRMNHzzGdTzG5hh75XbWullNSOpgVhZNNLSpj5OPVJbP2J6/6e2skrqGYdFEAy/8BcTTHD5/gN3bhrxTSVLXMCyaqbibad3bftGKQlJXMSyaYfa216P7oP9KGPShf5K6iwPczVC9lpOD1JK6kJVFs3jbq6QuZlgslDOuJfUQu6EWyhnXknqIlUUj7HqS1COsLCRJpQwLSVIpw0KSVMqwkCSVcoD7bHi7rKQeZWVxNrxdVlKPsrI4W94uK6kHGRZl7HqSJLuhar3qGdp2PUlSj1cW1VUDMHliipNHJhhITj1De/zv7XqS1PN6u7KorhqoPDN7JiuvX3mGttWEJPV4ZQGnVQ1Hx8bZvnOYqekZVixf5jO0JanQW2FR0+1UO2A90L+G3TuGGD54jKFNfT5DW5IKvRUW1U+0gzm7mAb61xgSklSjt8ICHKyWpAXo7QFuSVJdDAtJUqneCIuRXbDrXafdJitJql9vhEX1wLZzJiTprHXfAHft7bFwKigc2JakBem+yqJmVjZgRSFJDeq6ymLyxBQTr7mEo1d91vkSktQkSyssXvxmZaB6HpMnpsgj+3g++9m+c5jdO4YMDElqgq7qhpo4Ps1T2c8DJ9/K1PQMwwePtbtJktQVllZlccGlZxykfmUhwJnKQoBDm/oWsXGS1L2WVliUcCFASWqNrgoLcCFASWqFrhqzkCS1Rl1hERHXRsQzEXEgIm6f4/NVEfG54vOvRMTGYntfRDwcEd+PiE/Nc+y9EfGNRk5CktRapWEREecAdwHvBDYDN0bE5prdbgHGM/MS4JPAx4vtPwR+A/i1eY69Dfj+wpouSVos9VQWVwAHMvNgZr4M3A9srdlnK/CZ4vUe4OqIiMz8QWY+SiU0ThMRPwp8BLhzwa2XJC2KesJiHfB81ftDxbY598nMaeB7QNl9q78N/B7wUl0tlSS1TVvuhoqINwM/npkfnh3fOMO+twK3Fm9POL7xiguAF9vdiA7htTjFa3GK16LiNTThOtQTFoeBN1a9X19sm2ufQxGxHHgdcKbp0z8DDEbEc0UbXh8Rj2Tm22t3zMy7gbsBImIkMwfraHPX81qc4rU4xWtxiteiorgO1zZ6nHq6oR4HLo2IiyNiJXADsLdmn73ATcXr64GHMjPnO2Bm/lFmviEzNwJXAv8wV1BIkjpDaWWRmdMRcRvwIHAO8OnM3B8RdwAjmbkXuAe4NyIOAN+hEigAFNXDecDKiHgvcE1mPtX8U5EktUpdYxaZ+SXgSzXbPlb1+ofAz8/z3Y0lx34O+Ml62kHRHSXAa1HNa3GK1+IUr0VFU65DnKG3SJIkwOU+JEl16IiwqGM5kQ3FsiFfjYivR8R1xfaNEXE8Ir5W/Pzx4re+uRZ6LYrP/kVEfDki9kfEvoh4zeK2vrka+HPxwao/E1+LiJnidu0lq4FrsSIiPlP8eXg6In598VvfXA1ci5URsau4Fk9GxNsXvfFNVse16I+Ivy6uwyMRsb7qs5si4pvFz021332VzGzrD5VB82eBTcBK4Elgc80+dwO/UrzeDDxXvN4IfKPd59Ah12I58HXgTcX7PuCcdp9TO65FzT5bgGfbfT5t/HPxAeD+4vVrgeeAje0+pzZdi18FdhWvXw+MAsvafU4tvhZ/BtxUvL4KuLd4/WPAweLXNcXrNWf6/TqhsqhnOZGkckcVVOZwfGsR27eYGrkW1wBfz8wnATLzWGaeXIQ2t0qz/lzcWHx3KWvkWiTwI8X8p3OBl4GJ1je5ZRq5FpuBhwAy85+A7wJLeR5GPdfilXMGHq76/N8Cf5WZ38nMceCvgDPOxeiEsKhnOZHfAn4hIg5RuSvrP1Z9dnFRbv6/iHhbS1vaeo1ci8uAjIgHI+KJiPivrW5sizX652LW+4H7WtHARdTItdgD/AA4Avwj8InM/E5LW9tajVyLJ4H3RMTyiLgYGOD0CcdLTT3X4klgW/H654DVEdFX53dP0wlhUY8bgT/JzPXAdVTmdCyj8j/Ahsz8KSqLEv7PiDjvDMfpBvNdi+VUJjh+sPj15yLi6vY1c1HMdy0AiIi3AC9lZi8sETPftbgCOAm8AbgY+M8Rsal9zVwU812LT1P5S3EE+H3g76hcm272a8DPRsRXgZ+lstrGgs65E8KinuVEbgE+D5CZX6ay1skFmXkiM48V20ep9N9d1vIWt86CrwWV/wn+JjNfzMyXqPyL6qdb3uLWaeRazLqBpV9VQGPX4gPAX2bmVNH18hhLu+ulkb8vpjPzw5n55szcCpwP/MMitLlVSq9FZn4rM7cV/6D+aLHtu/V8t1YnhEU9y4n8I3A1QET8cyr/8V+IiLVRed4Gxb+WLqUyULNULfhaUJlhvyUiXlv0T/8ssJRnyjdyLSj+JfnvWPrjFdDYtfhHKgObRMSPAEPA3y9Su1uhkb8vXltcAyLiHcB0Lu3VJEqvRURcUFVt/zqV6goqf19cExFrImINlTHPB8/4u7V7RL8Ymb+OSsI/C3y02HYH8J48dUfDY1T6375GZckQgPcB+4ttTwDvbve5tOtaFJ/9QnE9vgH8brvPpc3X4u3AcLvPod3XAvhRKnfE7Kfyj4f/0u5zaeO12Ag8AzwN/F+gv93nsgjX4nrgm8U+O4FVVd/9JeBA8XNz2e/lDG5JUqlO6IaSJHU4w0KSVMqwkCSVMiwkSaUMC0lSKcNCklTKsJAklTIsJEml/j/Ph/TIxj9AXQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# same plot, zoomed: sklearn and our computation agree\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.plot(effs_ph, effs_el, '.')\n",
    "ax.plot(tpr, fpr)\n",
    "ax.set_xlim(0.85, 0.9)\n",
    "ax.set_ylim(0.014, 0.025)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '1 - fpr = electron efficiency')"
      ]
     },
     "execution_count": 89,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAFzCAYAAAAzNA41AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debydVX3v8c937zNknkgYzEACDWBERBpxFq6gAleJl2oLVasUwaGU3qvlFltfDlivtRRv6xUrWKmIA4JamwqUVgSpCkookDKIxDAkgGQgZD7DPvt3/3iefbJzcoad5Kw9ZH/fr9d+7Wfaz/o9GfZvr2c9ay1FBGZm1r4KjQ7AzMway4nAzKzNORGYmbU5JwIzszbnRGBm1uacCMzM2lxHowPYW7Nnz46FCxc2Ogwzs5Zyzz33bIiIOcPta7lEsHDhQlasWNHoMMzMWoqkJ0ba51tDZmZtzonAzKzNORGYmbU5JwIzszbnRGBm1uacCMzM2pwTgZlZm3MiMDNrc04EZmZtLlkikHS1pHWSHhhhvyR9XtIqSSslnZAqFjMzG1nKGsFXgdNG2X86sDh/XQD8fcJYzMxsBMnGGoqIOyQtHOWQZcDXIps0+S5JMyQdFhHPpIjnyY07+PWGbYPrqrxLw2yjapuG2TZ0YfTjditDe3y06jNjHacRP6thPsswx0XAhM4iUnYeSfl7do7B82n4fcoDqF6XRCHfhqC7I/t9UdCuYwaXhwvOzBqqkYPOzQXWVK2vzbftkQgkXUBWa2DBggX7VNjNDzzDZ27+5T591sbXHsmhklTy5UIlCQmKBdFR0G7Jpnp/Ych79XkLedLZuL2Pw2dNoljQ4Kug7Lw9pQEKErOndFOQKBayMkFs7y0xb+ZECnnZ5O+FPEEWClVls3uy29rTzwumT8ziL4qOQoFiAUDMnNRJR7HAxM4ik7uLWSxFUdSu+DqLBTqLhV3lOYlaQi0x+mhEXAVcBbB06dLYl3O89aVzOXHRrOx8g+fdrZQ9tg13XOQr1R+t7A/2/PBYx8Wwx+1ZBqN9dpiYhzvP8zv6mNjVMbgvIjtfBPlyvj3/0OA+2P3YvNDB7flyOeCZ53cya0rXYEzl8q7jyoOf2VV2OYacO99XrhwfQam863Pl8q7zlGPXZ8pDzl8uZ+d8bMN2jjl0KgPlYKAc9JXKDERQLgc9/WWezuN9fON2ymWy4yJYv7U3S0TSbmU1UiXpVBJhYXBdRATb+waYO2MinUXRUSwQkV3jC2ZMGDyuWBCFQvbZorLl9Vt7WXjQpDzhZMlw6HIlqW7tyZJjRyEro6MgdvQNMHNyJxM6i7sls0rCLZXLHDZ9ItMndjBtQifdHUW6OwtM6Cw29g/UBjUyETwFzK9an5dvS+KQaRM4ZNqEVKe3NrFb0hn6TvY+MBD05gmnNFCmVA5KA8GWnn4A+ktlnt7cQ2dRDJSzRFeuvEewtafEjr4Sk7o6KJdjt6RYWS4PJsBsec1zO5g5qSsrq1ymr1Tm1+u3sWj2ZAby85bKZXpLu84xUA4e37CdOVO7uXdHX3aucuX8wUB5VxkD5WBLT2nc/zyndHdQEOzsH+DwgybTVSzQURQ7+wY4Ys5k5s2cRGmgzAtmTKS7o0BPqczhsybRUSwwb+ZE5s+axJTulvg929Qa+Se4HLhQ0nXAy4HNqdoHzMaLJIqCIu17mybyxFDKXzv7BiiVy4O1rl2JJ0+AO/vpGyizcVsfPaUBevvLPL5xO1O6O+gtlXl03TZmTeqkb6BMXyl4+JktPLe9j/XbernlwWdriqlYEIsPnkJ3Z5GuougrlTlk2gRmT+3myDlT6OooIOCIOZM5YvYUDprSRWfRT89XJEsEkr4FnAzMlrQW+DjQCRARXwJuAs4AVgE7gHNTxWJm40d5m0ZHfmcn5S/yyq3BvlKZ3lKZLTv72dE3wG+27OSxDTt4bMM21m/tpVgQ/XlN7NFnt9JXKnP/2s1jnn/R7MksPGgSC2ZNorNY4ITDZzJrchfHzZvOpK72qWkoGn3jcy8tXbo0PEOZmdViZ98AfaUy2/pKPP38Th5bv51tvSUeeGozj67bxhMbt1MoiOd39O/xWQlef/TBnHD4TE4/9lCOmDOlAVcwfiTdExFLh93nRGBm7a5cDp7d2sNTm3by9OYeblz5NA8/s5Unn9ux23EfPPlI5s2cxNkvm0+h0Fq3B50IzMz2QW9pgEef3cZVd6zm3jWbWPPczsF9x8+fwfxZkzjm0Km8+1ULm77R2onAzGwcbNjWy8f/+UEeXbeVJzbuoLdUHtx39CFTueiUxZx+7KFNWVtwIjAzS6Cnf4AbVz7Dt37xJCue2DS4fe6MifzV77yYly86iK6O5ng6yYnAzCyxNc/t4OYHnuFvf/goO/sHBjsgLjlsGt+64BVMn9jZ0PicCMzM6mhbb4mrf/IYn/v3Xw1u+8jpx/C+k45sWEyjJYLmqLOYmR1ApnR3cNEpi1n9f87gI6cfA8Bnbv4lH/j6PQ2ObHhOBGZmiRQK4n0nHclNF70WgJsf+A3nffVutvWO/3Ad+8OJwMwssSUvmMYv/vwU3vSiQ7j1l+s49uO3sGrd1kaHNciJwMysDg6eNoEr37WUj/73FwJw6ufuYGffQIOjyjgRmJnV0XtfewSnH3soAMuu+EmDo8k4EZiZ1dkX35FN0f6rZ7exbmtPg6NxIjAzqztJXP72lwDw4evvb3A0TgRmZg1x1glzAfiPRzewvcFPETkRmJk1gCT+4JWHA3D+1xrbSdaJwMysQS5ddiyvP+ZgfvbrjVx71xMNi8OJwMysgf7u7OMB+O49axsWgxOBmVkDTZ3QycFTu7lvzfM0auw3JwIzswb77cNnAvD1Bt0eciIwM2uwT731WABuf2R9Q8p3IjAza7DZU7qZN3Mit/5yXUPKdyIwM2sCB03uAmBLT3/dy3YiMDNrAu959UIAbn342bqX7URgZtYEXrZwFgArHt80xpHjz4nAzKwJzJs5CYAbVtS/P4ETgZlZk3jJvOn0DZQpl+vbn8CJwMysSZx09MEA3Lf2+bqW60RgZtYkTjpqDgB3P/ZcXct1IjAzaxLHzZsOwCO/qe98xk4EZmZNorOYfSX3DpTrWq4TgZlZEzlhwQzWPLejrmU6EZiZNZHZU7rpK7lGYGbWtjo7CjzpGoGZWfua3FVkR98Aq9dvq1uZTgRmZk3kDUsOBeCJOtYKnAjMzJrI7CnZKKQPPrW5bmU6EZiZNZGjD50KQD1nrXQiMDNrIhM7iwD013G8IScCM7MmIgkJHtuwvW5lOhGYmTWZCOr6CKkTgZlZk5k9pZve/oG6ledEYGbWZH7r4Mk8sdE1AjOztnXotAns7B8g6vTokBOBmVmTmdiVPTnUW6cxh5wIzMyazBGzpwBQqtMjpE4EZmZNpqMoAHrq1GDsRGBm1qTqNVOZE4GZWZM5+pBsmInHN9anU5kTgZlZkzny4KyNoOw2AjOz9jS5uwOAVevqMydB0kQg6TRJj0haJemSYfYvkHSbpHslrZR0Rsp4zMxawZQ8ERQKqkt5yRKBpCJwBXA6sAQ4R9KSIYd9FLg+Il4KnA18MVU8ZmatZP6sidz2y3V1KStljeBEYFVErI6IPuA6YNmQYwKYli9PB55OGI+ZWcuYNqGTtZt21qWslIlgLrCman1tvq3aJ4B3SloL3AT8ccJ4zMxaxmHTJ7ZNh7JzgK9GxDzgDOBaSXvEJOkCSSskrVi/fn3dgzQzq7cXHja1bmWlTARPAfOr1ufl26qdB1wPEBF3AhOA2UNPFBFXRcTSiFg6Z86cROGamTWPYt5QXI9HSFMmgruBxZIWSeoiawxePuSYJ4FTACS9kCwR+Ce/mbW9orJEMFCHEUiTJYKIKAEXArcAD5M9HfSgpEslnZkf9mHgfEn3A98C3hP1GnfVzKyJVR4dHahDjaAj5ckj4iayRuDqbR+rWn4IeHXKGMzMWlHl1lDfQJkJ+YT2qTS6sdjMzIZRGXl0847+5GU5EZiZNaFFsycD0D+QfnIaJwIzsybUUci+nrf0lJKX5URgZtaEgqyReJsTgZlZezpk2gRgV0JIyYnAzKwJFVS/x0edCMzMmtBgz+JW7lBmZmb7rmOwQ1n6spwIzMyaUOXW0FObdqQvK3kJZma21yqzlHUn7lUMTgRmZk1pQmf29ew2AjOzNiVVGovTl+VEYGbWhCrz1tdjQGYnAjOzJlRpLG71iWnMzGwfFXxryMysvVVmb3djsZlZm6pMVblhW1/yssZMBJLukfRHkmYmj8bMzADo6si+nisT1KRUS43g94AXAHdLuk7Sm1R5rsnMzJLoLBaY1FXk8Y3bk5c1ZiKIiFUR8RfAUcA3gauBJyR9UtKs1AGambWrHX0DTGyWnsWSjgMuBy4Dvgu8HdgC/ChdaGZm7e2YQ6dSqsNjQx1jHSDpHuB54CvAJRHRm+/6uaRXpwzOzKyddRRVl/kIxkwEwNsjYvVwOyLirHGOx8zMch2FAk8/vzN5ObXcGnqvpBmVFUkzJf1lwpjMzAzY0Vdi044meHwUOD0inq+sRMQm4Ix0IZmZGcCL585ApH9Is5ZEUJTUXVmRNBHoHuV4MzMbB92dBUrl9FOU1dJG8A3gVkn/mK+fC1yTLiQzMwPoLKguPYvHTAQR8VlJK4FT8k2fiohb0oZlZmZbe0oA9A+U6SymGxGolhoBEXEzcHOyKMzMbA/HHDYV7s2GmUiZCGoZa+gsSY9K2ixpi6StkrYki8jMzAAGG4ofejrtV24tNYK/Bt4SEQ8njcTMzHZz9KFTAdjeV0paTi11jWedBMzM6m/W5C4A+gfS9i6upUawQtK3ge8DleEliIjvJYvKzMwGh6J+bnvaJ4dqqRFMA3YAbwTekr/enDIoMzODGZM6AbjvyefHOHL/1PL46LlJIzAzs2EdPHUCAJO60w5FXctTQ0dJulXSA/n6cZI+mjQqMzMDslpB6hFIa7k19GXgI0A/QESsBM5OGZSZmWWKSj8UdS2JYFJE/GLItrTPMpmZGQCFgihH4xPBBklHAgEg6W3AM0mjMjMzADoKSv7UUC2Pj/4RcBVwjKSngMeAdyaNyszMAOgtlekrpR2BtJanhlYDp0qaDBQiYmvSiMzMbNDcGROTlzFiIpD0zoj4uqQPDdkOQER8LnFsZmZtr1AQiTsWj1ojmJy/T00bgpmZjaQoKCd+amjERBARV+bvn0wagZmZjajYDE8NSbpmmMnrr04alZmZAdnt+M07+5OWUcvjo8cNM3n9S9OFZGZmFdt7Szy7pSdpGbUkgoKkmZUVSbOocWYzMzPbPwdN6WbqhM6kZdTyhX45cKekGwABbwM+nTQqMzMD4KDJXTy+YXvSMmrpR/A1SSuA1+ebzoqIh5JGZWZmQNZYnLqNYLR+BNMiYkt+K+g3wDer9s2KiOeSRmZmZuzoKzUuEZB98b8ZuId8nKGc8vUjEsZlZmbAQZO7yfvxJjNaY/Ff5e8vjIgjql6LIqKmJCDpNEmPSFol6ZIRjvldSQ9JelDSN4c7xsysXU2b2EExcSYYLRH8Xf7+s305saQicAVwOrAEOEfSkiHHLCab6+DVEfEi4H/uS1lmZgeqgtJ3KBvt1lC/pKuAeZI+P3RnRFw0xrlPBFblg9Yh6TpgGVDd0Hw+cEXeN4GIWLc3wZuZHegEJB5hYtRE8GbgVOBNZO0Ee2susKZqfS3w8iHHHAUg6adAEfhERPzr0BNJugC4AGDBggX7EIqZWWuqGuhzcHm8jZYILo6IP5O0ICKuSVJ6Vv5i4GRgHnCHpBdX92QGiIiryOZEYOnSpYlzo5lZ8ygMJgKSNRqP1kZwhrL0s6/zEz8FzK9an5dvq7YWWB4R/RHxGPArssRgZmZAIf/yT9lOMFoi+FdgE3CcpC2Stla/13Duu4HFkhZJ6iJLKMuHHPN9stoAkmaT3SpavbcXYWZ2oNJgIkhXxoiJICIujogZwI0RMS0ipla/j3XiiCgBFwK3AA8D10fEg5IulXRmftgtwEZJDwG3kd2O2rjfV2VmdoCotAukrBHUMsTEMkmHA4sj4oeSJgIdtUxZGRE3ATcN2faxquUAPpS/zMxsiC09Wa/izTv7mdBZTFJGLfMRnA98B7gy3zSP7JaOmZkltuigbLLIRrURVPwR8GpgC0BEPAocnCwiMzMb1NA2giq9EdFXWZHUwe5jD5mZWSJiVz+CVGpJBD+W9OfARElvAG4A/iVZRGZmtkteI0g5ykQtieASYD3wX8D7yBp/P5ouJDMzq0g88ChQ21NDZeDL+cvMzOqoumdxsjLSndrMzPaXGtyz2MzMGqySCFI+oeNEYGbWxOrx1NCYbQSSjgIuBg6vPj4iXj/ih8zMbFzUo0YwZiIge1z0S2SNxQMJYzEzsyFUh8biWhJBKSL+Pl0IZmY2ksrjo43uUPYvkj4o6TBJsyqvZBGZmdmgyq2hgUa2EQDvzt8vrtoWwBHjH46ZmVXr7S8DsH5rL8ccmqaMWjqULUpTtJmZjWXBQZOAXU8PpVDLU0OdwAeA1+WbbgeujIj+ZFGZmRkAncXsDn7/QDlZGbXcGvp7oBP4Yr7+rnzbe1MFZWZmmcld2WQ0lQlqUqglEbwsIl5Stf4jSfenCsjMzHaZ1J19Tff0p3t6v5anhgYkHVlZkXQE7k9gZlYXxfyxoYR3hmqqEfwpcJuk1WSPtB4OnJsuJDMzqyjkP9cb9viopCLwEmAxcHS++ZGI6E0WkZmZDarUCMoJ56oc9dZQRAwA50REb0SszF9OAmZmdVIsVG4NNbZD2U8lfQH4NrC9sjEi/jNZVGZmBuxKBKs3bEtWRi2J4Pj8/dKqbQF49FEzs8SmTujc7T2FWhLBeRGxunpD/uSQmZnVwcTOYtJbQ7U8PvqdYbbdMN6BmJnZ8DoKakzPYknHAC8Cpks6q2rXNGBCsojMzGw3Eqx5bmey8492a+ho4M3ADOAtVdu3Aucni8jMzHbT1VFkQme6mYVHTAQR8c/AP0t6ZUTcmSwCMzMb1fSJHQ2fvP79kmZUViTNlHR1wpjMzKyKpIbPUHZcRDxfWYmITcBLk0VkZma7KSjtnMW1JIKCpJmVlXyayloeOzUzs3FQkCg3eKrKy4E7JVUeGX078OlkEZmZ2R4SdiOoaarKr0lawa6exGdFxEPpQjIzs2oFqeG3hgBmAdsj4gvAekmex9jMrE4kGttYLOnjwJ8BH8k3dQJfTxaRmZntpiA1/PHR/wGcST7yaEQ8DUxNGJOZmVWRSNpYXEsi6IusThJZQJqcLBozM9uDaPycxddLuhKYIel84IfAl5NFZGZmu9nZP8ATG3ckO38tTw39jaQ3AFvIxh/6WET8e7KIzMxsNxM6i3R1NGCsoWr5F7+//M3MGuDgqd08/XxPsvOPNgz1Vhi2oVpARMS0ZFGZmdmghvUsjgg/GWRm1gRSJ4KabjpJeo2kc/Pl2e5QZmZWP8WCGjtV5TAdyrpwhzIzs7opFMRz2/vSnb+GY9yhzMysgTZt72NLTynZ+d2hzMysyc2dMZFJncVk53eHMjOzJjehs0ChoGTnd4cyM7MmJ6VLAuAOZWZmbS9dn2VA0mmSHpG0StIloxz3O5JC0tKU8ZiZtapGT16/TyQVgSuA04ElwDmSlgxz3FTgT4Cfp4rFzMxGlrJGcCKwKiJWR0QfcB2wbJjjPgV8Fkg3kIaZWYtr9MQ0SPrb6vcazQXWVK2vzbdVn/cEYH5E3DhG+RdIWiFpxfr16/ciBDOz1pe4rbjmGsHr8veTxqtgSQXgc8CHxzo2Iq6KiKURsXTOnDnjFYKZWetogsnr98VTwPyq9Xn5toqpwLHA7ZIeB14BLHeDsZnZ7kTaKkHKRHA3sFjSIkldwNnA8srOiNgcEbMjYmFELATuAs6MiBUJYzIzsyGSJYKIKAEXArcADwPXR8SDki6VdGaqcs3MDkQpG4tr6lC2ryLiJuCmIds+NsKxJ6eMxcysVTVLY/E38/dvpArEzMxG1vAOZRHxN9XvZmZWP4krBGmHmDAzs/HR8A5lZmbWOA1tI5BUkPSqtCGYmdlYEjYRjJ4IIqJMNnCcmZk1SOr5CGq5NXRrPkx06vYKMzNrgFoSwfuAG4A+SVskbZW0JXFcZmZWJRI2F9cyVeXUZKWbmdmYUt+OqalnsaSzgNeQPcH0HxHx/aRRmZnZbhrWWAwg6YvA+4H/Ah4A3i/JDchmZvWSuEpQS43g9cALI+/fLOka4MGkUZmZ2W4a3aFsFbCgan1+vs3MzOog9XwEtdQIpgIPS/oFWVI6EVghaTlARHhIaTOzFjZiIpDUHRG9wLDDRpuZWR0lvDc0Wo3gTuAE4L0R8a50IZiZ2WhSd+cdLRF0Sfp94FX546O7iYjvpQvLzMyqNapD2fuBdwAzgLfsERM4EZiZ1UHDOpRFxE+An0haERFfSRyHmZmNoqEdypwEzMwaq1nmLDYzswZqdIcyMzNroNQdyvYqEUj6RKI4zMysQfa2RuBexGZmDRAJW4v3NhF4ljIzszprtsbi304ShZmZjappGovzyezNzKyOUt+K8VNDZmbNTiICyuU09QInAjOzJje5qwjA1t5SkvPvUyKQdO54B2JmZsPrKOZf1YkaCva1RvDJcY3CzMxG1LBB5yStHGkXcEiacMzMbCSphqIebRjqQ4A3AZuGbBfwsyTRmJnZHho5Mc0PgCkRcd/QHZJuTxaRmZkNK1Xn4tHmIzhvlH2/nyYcMzMbyv0IzMwMSNe72InAzKzJKXEjgROBmVmLSDUCqROBmVmTa7bRR83MrEHcRmBm1qb81JCZmQHp+hE4EZiZNTs/NWRmZpBurCEnAjOzJuc2AjMzy7iNwMysPbkfgZmZAe5HYGbWtpS4lcCJwMysRbRkPwJJp0l6RNIqSZcMs/9Dkh6StFLSrZIOTxmPmVkratk2AklF4ArgdGAJcI6kJUMOuxdYGhHHAd8B/jpVPGZmra4V+xGcCKyKiNUR0QdcByyrPiAibouIHfnqXcC8hPGYmbWkVu5HMBdYU7W+Nt82kvOAmxPGY2bW0uo+Z3E9SXonsBQ4aYT9FwAXACxYsKCOkZmZNV7LthEATwHzq9bn5dt2I+lU4C+AMyOid7gTRcRVEbE0IpbOmTMnSbBmZs2uFfsR3A0slrRIUhdwNrC8+gBJLwWuJEsC6xLGYmbWslq2H0FElIALgVuAh4HrI+JBSZdKOjM/7DJgCnCDpPskLR/hdGZmbS/VnMVJ2wgi4ibgpiHbPla1fGrK8s3MDggt3EZgZmbjqCV7FpuZ2f5r5X4EZmbWApwIzMyanDxnsZmZgdsIzMzaVmcxqxH0l8tJzu9EYGbW5CZ0FgHo6R9Icn4nAjOzJtdRyGoEA+XWG4bazMzGkdsIzMzaVCuPPmpmZuOgMuhcK44+amZm4yGvEaQadM6JwMysyVXuDLlGYGbWpio9i91YbGbWpna1FfvWkJlZW9JgG0Ga8zsRmJk1OT81ZGbW5lwjMDMzwI+Pmpm1LT8+ambW7nxryMysve1qLPatITOztqTE94acCMzMmpzbCMzM2pyHmDAza3OD/QjcRmBm1p4KeSZINFOlE4GZWbOrzFlcGignOb8TgZlZk+so5onAk9ebmbWnjkL2VV0acCIwM2tLu2oEvjVkZtaWdrURuEZgZtaWOorZV/WA2wjMzNpTpUbQ71tDZmbtyRPTmJlZUk4EZmZtzonAzKzNORGYmbU5JwIzszbnRGBm1uacCMzM2pwTgZlZm3MiMDNrc04EZmZtzonAzKzNORGYmbU5JwIzszbnRGBm1uacCMzM2pwTgZlZm0uaCCSdJukRSaskXTLM/m5J3873/1zSwpTxmJnZnpIlAklF4ArgdGAJcI6kJUMOOw/YFBG/Bfxf4LOp4jEzs+GlrBGcCKyKiNUR0QdcBywbcswy4Jp8+TvAKVJlUjYzM6uHlIlgLrCman1tvm3YYyKiBGwGDhp6IkkXSFohacX69esThWtm1py6igVeu3g2h02fkOT8HUnOOs4i4irgKoClS5cmmr7ZzKw5zZjUxbXnvTzZ+VPWCJ4C5letz8u3DXuMpA5gOrAxYUxmZjZEykRwN7BY0iJJXcDZwPIhxywH3p0vvw34UUT4F7+ZWR0luzUUESVJFwK3AEXg6oh4UNKlwIqIWA58BbhW0irgObJkYWZmdZS0jSAibgJuGrLtY1XLPcDbU8ZgZmajc89iM7M250RgZtbmnAjMzNqcE4GZWZtzIjAza3NOBGZmbc6JwMyszTkRmJm1OScCM7M2p1Yb2kfSeuCJffz4bGDDOIbTCnzN7cHX3B7255oPj4g5w+1ouUSwPyStiIiljY6jnnzN7cHX3B5SXbNvDZmZtTknAjOzNtduieCqRgfQAL7m9uBrbg9Jrrmt2gjMzGxP7VYjMDOzIQ7IRCDpNEmPSFol6ZJh9ndL+na+/+eSFtY/yvFVwzV/SNJDklZKulXS4Y2IczyNdc1Vx/2OpJDU8k+Y1HLNkn43/7t+UNI36x3jeKvh3/YCSbdJujf/931GI+IcL5KulrRO0gMj7Jekz+d/HislnbDfhUbEAfUimxbz18ARQBdwP7BkyDEfBL6UL58NfLvRcdfhmv8bMClf/kA7XHN+3FTgDuAuYGmj467D3/Ni4F5gZr5+cKPjrsM1XwV8IF9eAjze6Lj385pfB5wAPDDC/jOAmwEBrwB+vr9lHog1ghOBVRGxOiL6gOuAZUOOWQZcky9/BzhFkuoY43gb85oj4raI2JGv3gXMq3OM462Wv2eATwGfBXrqGVwitVzz+cAVEbEJICLW1TnG8VbLNQcwLV+eDjxdx/jGXUTcQTaH+0iWAV+LzF3ADEmH7U+ZB2IimAusqVpfm28b9piIKAGbgYPqEl0atVxztfPIflG0sjGvOa8yz4+IG+sZWEK1/D0fBRwl6aeS7pJ0Wt2iS6OWa/4E8E5Ja8nmSP/j+oTWMHv7/31MSSevt+Yj6Z3AUuCkRseSkqQC8DngPQ0Opd46yG4PnUxW67tD0osj4vmGRpXWOcBXI+JySa8ErpV0bESUGx1YqzgQawRPAfOr1ufl28CSgMwAAAV8SURBVIY9RlIHWXVyY12iS6OWa0bSqcBfAGdGRG+dYktlrGueChwL3C7pcbJ7qctbvMG4lr/ntcDyiOiPiMeAX5ElhlZVyzWfB1wPEBF3AhPIxuQ5UNX0/31vHIiJ4G5gsaRFkrrIGoOXDzlmOfDufPltwI8ib4VpUWNes6SXAleSJYFWv28MY1xzRGyOiNkRsTAiFpK1i5wZESsaE+64qOXf9vfJagNImk12q2h1PYMcZ7Vc85PAKQCSXkiWCNbXNcr6Wg78Qf700CuAzRHxzP6c8IC7NRQRJUkXAreQPXFwdUQ8KOlSYEVELAe+QlZ9XEXWKHN24yLefzVe82XAFOCGvF38yYg4s2FB76car/mAUuM13wK8UdJDwABwcUS0bG23xmv+MPBlSf+LrOH4Pa38w07St8iS+ey83ePjQCdARHyJrB3kDGAVsAM4d7/LbOE/LzMzGwcH4q0hMzPbC04EZmZtzonAzKzNORGYmbU5JwIzszbnRGB1IWmGpA82Oo4KSdv28vi3SlqSKp6qci6S9LCkb+Sj5P5Q0n2Sfk/SP4wWg6QzRxuF1WwkfnzU6iIf6vsHEXHsXn6uGBEDCeLZFhFT9uL4r5LF/53xjmVIOb8ETo2ItXlnob+MiFNTlmnmGoHVy18BR+a/bi+TdLKkOyTdmI81/6V8fCAkbZN0uaT7gVfua4GjlZHv/7Sk+/PB2Q7Jty2U9CPtmrdhgaRXAWcCl+XxHynp+PxzKyX9k6SZ+edvl/RZSb+Q9CtJrx0htosl3Z1//pP5ti+RDbd8s6Q/A74OvKyqzNsrQ2QoG6P/P/P4b823vUfSF/LlOZK+m5dxt6RX59s/oWy8+9slrZZ0UVVMf5DHc7+kayVNlfSYpM58/7TqdTuANHrsbb/a4wUspGp8dbKekz1kX3xF4N+Bt+X7AvjdEc5zMXDfMK/PD3PsWGW8JV/+a+Cj+fK/AO/Ol/8Q+H6+/NXKZ/P1lcBJ+fKlwN/my7cDl+fLZwA/HCauN5KNoS+yH2M/AF6X73scmF0V/w+qPnc72YCBc8hGn1yUb5+Vv78H+EK+/E3gNfnyAuDhfPkTwM+AbrLxeDaS9Vp9Edm4RLOHnPMfgbfmyxdUrs2vA+t1wA0xYS3lFxGxGga71b+GbH6IAeC7w30gIi4jGy5jf8voI/sCBrgHeEO+/ErgrHz5WrIksRtJ04EZEfHjfNM1wA1Vh3yv6rwLh4npjfnr3nx9CtnAcHfUeE2vAO6IbFA5ImK4setPBZZo1zQb0yRVboXdGNmgg72S1gGHAK8HboiIDUPO+Q/A/yYbw+hcsvkO7ADjRGCNNLSBqrLeEyO0C0i6GHjHMLvuiIiLhtk+Uhn9EVFZHmB8/y9URnYd6bwCPhMRV45jmUMVgFdExG4T8uSJoXrk2VGvPSJ+mt8uOxkoRsSw0ydaa3MbgdXLVrKhoaudqGxUyQLwe8BPxjpJRFwWEccP8xouCexLGT9j1yCE7wD+Y2j8EbEZ2FR1//9dwI+p3S3AH1Z+oUuaK+ngvfj8XcDrJC3KPz9rmGP+jaoJWiQdP8Y5fwS8XdJBw5zza2S3mv5xL2K0FuJEYHUR2QiYP5X0gKTKrZ27gS8ADwOPAf+UoOi9LeOPgXMlrST7gv+TfPt1wMXKJkg/kmwY88vy444nayeoSUT8G9kX652S/ovsVtXQJDna59eT3a//Xt6g/u1hDrsIWJo3/j4EvH+Mcz4IfBr4cX7Oz1Xt/gYwE/hWrTFaa/Hjo9YQ+a2GP42IN7dyGe1A0tuAZRHxrkbHYmm4jcDMRiTp/wGnkz0BZQco1wjMzNqc2wjMzNqcE4GZWZtzIjAza3NOBGZmbc6JwMyszTkRmJm1uf8PkQ3YzEcnR94AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we usually call \n",
    "#  eff_photon = P[reco-photon | true-photon]\n",
    "#  eff_electron = P[reco-electron | true-electron] = 1 - P[reco_photon | true-electron] = 1 - fpr\n",
    "# we have used the fact that P[reco-electron] + P[reco-photon] = 1, since the problem is binary\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.plot(tpr, 1 - fpr)\n",
    "ax.set_xlabel('tpr = photon efficiency')\n",
    "ax.set_ylabel('1 - fpr = electron efficiency')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the ambiguity tool.\n",
    "\n",
    "If we consider the ambiguos as photon then reco-photon := amb_type != 0\n",
    "   * tpr = photon efficiency = rate of amb_type != 0 on the true-photon sample\n",
    "   * 1 - fpr = electron efficiency = 1 - (rate of amb_type != 0) on true-electron sample = (rate of amb_type == 0) on true-electron sample\n",
    "\n",
    "    \n",
    "If we consider the ambiguos as electron then reco-photon := amb_type == 6 \n",
    "   * tpr = photon efficiency = rate of amb_type == 6 on the true-photon sample\n",
    "   * 1 - fpr = electron efficiency = 1 - (rate of amb_type == 6) on true-electron sample = (rate of amb_type != 6) on true-electron sample\n",
    "   \n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [],
   "source": [
    "eff_ph_amb_as_ph = (df_true_ph['ph_amb_type'] != 0).sum() / len(df_true_ph)\n",
    "eff_el_amb_as_ph = (df_true_el['el_amb_type'] == 0).sum() / len(df_true_el)\n",
    "\n",
    "eff_ph_amb_as_el = (df_true_ph['ph_amb_type'] == 6).sum() / len(df_true_ph)\n",
    "eff_el_amb_as_el = (df_true_el['el_amb_type'] != 6).sum() / len(df_true_el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9101553704120694 0.9204535838199205\n",
      "0.003377617653681603 0.9996615046119997\n"
     ]
    }
   ],
   "source": [
    "print(eff_ph_amb_as_ph, eff_el_amb_as_ph)\n",
    "print(eff_ph_amb_as_el, eff_el_amb_as_el)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'only when both particles are reconstructed')"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAGDCAYAAAAmphcsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZwcdZ3/8de7ZyaZ3CcYICEJGCABQhIjC+sBcskNP0BFBQVR1F0FdWGFXZHorrueKC64iKuAqCDiFa5FOQIiIAQhLKdACCQQIAkh92QyM5/fH1U96Qw9Mz3JVHfP9Pv5ePRj6vx+P9U9XZ+u+lZ9SxGBmZnVrlylAzAzs8pyIjAzq3FOBGZmNc6JwMysxjkRmJnVOCcCM7Ma50RQJSSFpLdmWP48SR/PqvyCeg6UtKQXy7tS0r/3VnnbStK7JD1dwnKnSbqnHDFZdSjXdywLTgRWNapx59kxQUfEnyJi90rGVKuy/LEkaVJafn0W5Vc7JwKzIvr6DmFb4++L298XY64WTgS9SNLU9PDwDUmPSzq2YN6Vki6VdJOkNZL+ImnXImW8XdKrkuoKpp0gaUGRZSendeXS8R9Jeq1g/tWSPlewykRJf07r/4OksQXL7ifp3rS8BZIOLJg3T9K/dbZuJ+/Fv0haLmmRpA8XTB8h6aeSlkl6QdKXJOUkTQUuA/aXtFbSGwXFjerufUvLzv+qO1PSy5KWSjqnYP6+ku5Lt3GppEskDSiYH5L+UdIzwDOS7k5nLUhj+kDHU1+SJkj6Tbo9KyRd0klse0j6o6TXJT0t6f0F846U9ES6fS8VxtyhjF0l3ZHWs1zSzyWNLJi/SNIXJT0KrJNU39XnWqT8Hq0vabSkK9L3eqWk3xXM+4SkZ9PtnStpxw7v86ckPZOWe6kkpfPeKukuSavSbfxlOr3TzyKN+RXgChU5qlTBkYSkQZK+k/7vrZJ0j6RBQL78N9Ly90+X/5ikJ9Ptu1XSxIJyD5X0VFrOJYA6e2+rXkT41QsvoAF4FvgXYABwELAG2D2dfyWwAtgXqAd+DlxbsH4Ab02HnwCOKJj3W+CfOqn3ReBt6fDTwEJgasG8menwPOA5YDdgUDr+9XTeTmlsR5L8ODg0Hd+uu3WLxHMg0AJcBAwEDgDWFbwPPwV+DwwDJgF/A85I550G3NOhvC7ftw7LTkrfx2uAIcDewDLgkHT+24D90nImAU8Cn+vwGfwRGA0M6vi5FGzfknS4DlgAfDetrxF4Z8dtSectBk5P654JLAempfOXAu9Kh0cBszrZvremn81AYDuSndf3CuYvAh4BJqSfU5efa5Hye7Q+cBPwyzTmBuCAdPpB6fbNSmP9L+DuDu/zjcBIYOf0Mzo8nXcN8K9pfe3vZxefRQvwjbSeQRT/Hyr8bl1K8v+7U/r5/X267qR0ufqC9Y4j+U5PTT+3LwH3pvPGkny/T0q3/fNpLB+v9L5oq/ZflQ6gv7yAdwGvALmCadcAc9LhK4H/KZh3JPBUwXjhP+sXgZ+nw6OB9cAOndR7NfAFYBxJIvgm8ClgMvBGPp70n/9LBev9A/C/BfVd3aHcW4GPdrdukXjyX84hBdOuAy5Iv3jNpDvAdN4ngXnpcLEvcZfvW4dl81/mPQqmfRP4cSfLfw74bYfP4KAOy3SVCPYn2YnVFym7fVuADwB/6jD/h8CF6fCL6fswvIf/c8cDDxeMLwI+VjDe5edapLyS1wd2ANqAUUXK+THwzYLxocAmYFLBe1q4g78OOC8d/ilwOTC+SLnFPotmoLHY+95xPZLksgHYp4v/ncJEcAvpj5R0PEfyXZwIfAS4v2CegCX00UTgU0O9Z0dgcUS0FUx7geSXR94rBcPrSb4gxfwMOEbSEOD9JDuRpZ0sexfJF+LdJL8Q55H8Cj8gXa8wns7qnwi8Lz1MfyM9LfNOki97T2MHWBkR6wrGXyB5f8aS/Hp6ocO8wveomJ7UDcmv7451I2k3STdKekXSauA/0pg6W7c7E4AXIqKlm+UmAn/X4f39MEnyBjiRJMG9kJ4W2b9YIZLeIuna9PTRapL/k67iL+Vz7ajU9ScAr0fEyiJl7EjBZxwRa0mOJEr5LvwzyU71ASWnVz/WRawAyyKiqZtl8saSHGU8V+LyE4GLC7b99TS2nUi/7/kFI8kGPfnfqSpOBL3nZWCC0vP1qZ2Bl3paUES8BNwHnACcSvKrvzN3kRyNHJgO3wO8gyQR3FVilYtJfvmNLHgNiYiv9zT21Kg0ieXtTPL+LCf5ZTixw7z8e9RbXeFOKFI3wH8DTwFTImI4yWm8jud1exLDYmBndd9IuRi4q8P7OzQiPg0QEQ9GxHHA9sDvSH4hF/MfaXx7p/Gf0k38W/O5lrr+YmB0YRtFgZcp+IzT/4UxlPBdiIhXIuITEbEjyVHSD9T1lUIdP691wOCCuscVzFsONAHF2piKfe6LgU922P5BEXEvyem89v+ztI1jQpEy+gQngt7zF5JfNv8sqSFtVDsGuHYry/spya+jvYHfdLZQRDxDcrh7CsnOZjXwKsmvzFITQf4I5L2S6iQ1pg1x47cydoCvSBog6V3A0cCvIqKVZCf3NUnD0oa3L6T1k8Y9XgUNuFvpAkmDJe1Jcl7+l+n0YcBqYK2kPYBPl1DWq8Auncx7gGSH8HVJQ9L37R1FlrsR2E3Sqen/RoOSiwKmpu/RhyWNiIhNaXxtRcrIx78WWCVpJ+DcbmLf1s+10/XTI9RbSHbUo9Jtene63jXA6ZJmSBpIksD+EhGLuqtQ0vsK4ltJsoPOvx9dfRZ5C4A907obgTn5GenR8U+AiyTtmG7T/mmMy9J6Csu/DDg//T/KX+jwvnTeTWk9J6Q/BM5i8xFen+NE0Esioplkx38EyS+PHwAfiYintrLI35L8qvptRKzvZtm7gBURsbhgXMBfS6koXe84kl/Iy0h+CZ3L1v9/vELyJX6ZpHH3UwXvw2dJfrUtJDl6+QXJlxPgDuBx4BVJy7eybki2/1ngduDbEfGHdPo5wIdIGvl+xOYE0ZU5wFXp6YH3F85IE9sxJOefXyQ5R/yBjgVExBrgMOBkkvfkFTY3cEJy1LcoPd3zKZLTRsV8haQBdhXJjqjTHwhpvdv0uZaw/qkkR3hPAa+RtLkQEbeRtAn9miRR7ppueyneDvxF0lpgLnB2RCxM582hk8+iIOa/AV8FbgOeIfkfK3QO8H/AgySner5B0o62Hvga8Oe0/P0i4rfp/GvTz+Yxku83EbEceB/wdZLTXlOAP5e4jVVHaUOHVSFJz5Ecmt5W6Vj6AkmTgOeBhhLO25tZykcEVUrSiSSHxXdUOhYz6998J14VkjQPmAac2uGqHzOzXudTQ2ZmNc6nhszMapwTgZlZjetzbQRjx46NSZMmVToMM7M+5aGHHloeEdsVm9fnEsGkSZOYP39+pcMwM+tTJL3Q2TyfGjIzq3FOBGZmNc6JwMysxjkRmJnVOCcCM7Ma50RgZlbjnAjMzGqcE4GZWY1zIjAzq3GZJQJJP5H0mqTHOpkvSd+X9KykRyXNyioWMzPrXJZHBFcCh3cx/wiSx7tNAc4kebB4ZlbdcAPPHHQwT06dxjMHHcyqG27Isjozsz4js76GIuLu9NGBnTkO+GkkD0S4X9JISTukD8XuVatuuIG//sd3WdwwArYbDi2Qu/hqRq+EIfvt176c8n9FwTQVmdZxoOvlVDAxP1iwWME63S2nTtdVkXUpslwENDbUISXlSEr/JmW0l6fi85QGUDguiVw6DcHA+uT3RU6bl2kfLhacmVVUJTud24nkYdh5S9Jpb0oEks4kOWpg55137nFFr333e9wzZnd+vNfRW854AnjigR6XZ9vmTckhn1TS4Vw+CQnqcqI+py2STeH8XIe/heXm0qSzYl0zE0cPpi6n9ldOSblNLa3kJMYOHUhOoi6X1Ali3cYWxo8aRC6tm/RvLk2QuVxB3WyZ7NY0bWLHEYOS+OtEfS5HXQ5AjBrcQH1djkENdQwZWJfEUifqtDm+hrocDXW5zfU5iVqG+kTvoxFxOXA5wOzZs3v8SLWWpUt5z8C17LViYVJefoZyTLzmF/la0roK6m2vf4tYtiyjYH7w5pW7Wy6KLvfmOuhq3SIxFyvnjfXNDBpQ3z4vIikvgnQ4nZ6u1D4Ptlw2rbR9ejrcFrD0jQ2MHjqgPaa2ts3LtbWvs7nutuhQdjqvLb98BC1tm9dra9tcTltsXqetQ/ltbUmZzy9fxx7jhtHaFrS2Bc0tbbRG0NYWNG1q4+U03kUr1tHWRrJcBMvWbEwSkbRFXZWUTzr5RJhrHxcRwbrmVnYaOYiGOlFflyMi2cYdRza2L1eXE7lcsm6dkuFlazYyaczgNOEkybDjcD6prmlKkmN9LqmjPifWN7cyakgDjQ11WySzfMJtaWtjhxGDGDGonuGNDQysr2NgQ47GhrrKvqHWrpKJ4CVgQsH4+HRar6vfYQfGvPwyY5pWbzl9xx2ZMnFUFlVaP7VF0un4l+Rva2uwMU04La1ttLQFLa3B6qZNAGxqaePlVU001InWtiTRteX/RrCmqYX1zS0MHlBPW1tskRTzw23tCTAZXvz6ekYNHpDU1dZGc0sbzy1by+SxQ2hNy21pa2Njy+YyWtuCRcvXsd2wgTy8vjkpqy1fftDatrmO1rZgdVNLr7+fQwfWkxNs2NTKxDFDGFCXo75ObGhuZZfthjB+1GBaWtvYceQgBtbnaGppY+LowdTX5Rg/ahATRg9m6MA+8Xu2qlXyHZwLfEbStcDfAauyaB8A2P7zn2PpBV8mmprap6mxke0//7ksqrN+TBJ1gjpq9zRNpImhJX1taG6lpa2t/ahrc+JJE+CGTTS3trFibTNNLa1s3NTGohXrGDqwno0tbTzz2lpGD26gubWN5pbgyaWreX1dM8vWbuTWx18tKaa6nJiy/VAGNtQxoE40t7TxluGNjB02kF23G8qA+hwCdtluCLuMHcqYoQNoqPPV83mZJQJJ1wAHAmMlLQEuBBoAIuIy4GbgSOBZYD1welaxjDjmGCBpK2hZupT6HXZg+89/rn26mZVOaZtGfXpmJ8tf5PlTg80tbWxsaWP1hk2sb27lldUbeH75ep5fvpZlazZSlxOb0iOxZ15dQ3NLGwuWrOq2/MljhzBpzGB2Hj2YhrocsyaOYvSQAUwfP4LBA2rnSENR6ROfPTR79uzwE8rMrBQbmltpbmljbXMLL7+xgeeXrWPtxhYee2kVz7y2lhdWrCOXE2+s3/SmdSU4aPftmTVxFEfsNY5dthtagS3oPZIeiojZRec5EZhZrWtrC15d08RLKzfw8qombnr0ZZ5cuoYXX1+/xXL/cOCujB81mJPfPoFcrm+dHnQiMDPbChtbWnnm1bVcfvdCHl68ksWvb2ifN2PCSCaMHswe44bx0b+fVPWN1k4EZma9YPnajVz4+8d55rU1vLBiPRtb2trn7f6WYZx18BSO2GtcVR4tOBGYmWWgaVMrNz26lGseeJH5L6xsn77TyEF8/cS9+bvJYxhQXx1XJzkRmJllbPHr67nlsaV877Zn2LCptf0GxGk7DOeaM/djxKCGisbnRGBmVkZrN7bwk3ue56I//q192vlH7MEnD9i1YjF1lQiq45jFzKwfGTqwnrMOnsLC/ziS84/YA4D/vOUpPv2zhyocWXFOBGZmGcnlxCcP2JWbz3oXALc89gpnXPkgazf2fncd28KJwMwsY9N2HM4D/3Iw793zLdz+1GvsdeGtPPvamkqH1c6JwMysDLYf3sgPT53Nl46aCsAhF93NhubWCkeVcCIwMyujj79rF47YaxwAx116T4WjSTgRmJmV2Q8+nDyi/W+vruW1NU3dLJ09JwIzszKTxHfetw8A/3TdggpH40RgZlYRJ8zaCYA/PbOcdRW+isiJwMysAiTxkf0nAvCJn1b2JlknAjOzHrpp4U0cdv1hTL9qOoddfxg3Lbxpq8r56nF7cdAe23Pvcyu4+v4XejnK0jkRmJn1wE0Lb2LOvXNYum4pQbB03VLm3Dtnq5PBxSfPAODXDy3pzTB7xInAzKwHLv7rxTS1bnmlT1NrExf/9eKtKm9YYwPbDxvII4vfoFJ9vzkRmJn1wCvrXunR9FK8beIoAH5WodNDTgRmZj0wbsi4Hk0vxb8dvxcA855ettVlbAsnAjOzHjh71tk01jVuMa2xrpGzZ5291WWOHTqQ8aMGcftTr21reFuluh+yaWZWZY7a5SggaSt4Zd0rjBsyjrNnnd0+fWuNGTKAJSs3sLppE8Mby/sQGycCM7MeOmqXo7Z5x9/Rae+YxOd/uYDbn3yV/zdzfK+W3R2fGjIzqwJvnzQagPmLVnazZO9zIjAzqwLjRw0G4Ffzy38/gROBmVmV2Gf8CJpb22hrK+/9BE4EZmZV4oDdtwfgkSVvlLVeJwIzsypxwG7bAfDg86+XtV4nAjOzKjF9/AgAnn6lvM8zdiIwM6sSDXXJLnlja1tZ63UiMDOrIrN2Hsni19eXtU4nAjOzKjJ26ECaW3xEYGZWsxrqc7zoIwIzs9o1ZEAd65tbWbhsbdnqdCIwM6sih05LurN+oYxHBU4EZmZVZOzQAQA8/tKqstXpRGBmVkV2HzcMgHI+tdKJwMysigxqqANgUxn7G3IiMDOrIpKQ4Pnl68pWpxOBmVmViaCsl5A6EZiZVZmxQweycVNr2epzIjAzqzJv3X4IL6zwEYGZWc0aN7yRDZtaiTJdOuREYGZWZQYNSK4c2limPoecCMzMqswuY4cC0FKmS0idCMzMqkx9nQBoKlODsROBmVmVKteTypwIzMyqzO5vSbqZWLSiPDeVORGYmVWZXbdP2gja3EZgZlabhgysB+DZ18rzTIJME4GkwyU9LelZSecVmb+zpDslPSzpUUlHZhmPmVlfMDRNBLmcylJfZolAUh1wKXAEMA34oKRpHRb7EnBdRMwETgZ+kFU8ZmZ9yYTRg7jzqdfKUleWRwT7As9GxMKIaAauBY7rsEwAw9PhEcDLGcZjZtZnDG9sYMnKDWWpK8tEsBOwuGB8STqt0BzgFElLgJuBz2YYj5lZn7HDiEE1c0PZB4ErI2I8cCRwtaQ3xSTpTEnzJc1ftmxZ2YM0Myu3qTsMK1tdWSaCl4AJBePj02mFzgCuA4iI+4BGYGzHgiLi8oiYHRGzt9tuu4zCNTOrHnVpQ3E5LiHNMhE8CEyRNFnSAJLG4LkdlnkROBhA0lSSROCf/GZW8+qUJILWMvRAmlkiiIgW4DPArcCTJFcHPS7pq5KOTRf7J+ATkhYA1wCnRbn6XTUzq2L5S0dby3BEUJ9l4RFxM0kjcOG0LxcMPwG8I8sYzMz6ovypoebWNhrTB9pnpdKNxWZmVkS+59FV6zdlXpcTgZlZFZo8dggAm1qzfziNE4GZWRWqzyW759VNLZnX5URgZlaFgqSReK0TgZlZbXrL8EZgc0LIkhOBmVkVyql8l486EZiZVaH2O4v78g1lZma29erbbyjLvi4nAjOzKpQ/NfTSyvXZ15V5DWZm1mP5p5QNzPiuYnAiMDOrSo0Nye7ZbQRmZjVKyjcWZ1+XE4GZWRXKP7e+HB0yOxGYmVWhfGNxX38wjZmZbaWcTw2ZmdW2/NPb3VhsZlaj8o+qXL62OfO6uk0Ekh6S9I+SRmUejZmZATCgPtk95x9Qk6VSjgg+AOwIPCjpWknvVf66JjMzy0RDXY7BA+pYtGJd5nV1mwgi4tmI+FdgN+AXwE+AFyR9RdLorAM0M6tV65tbGVQtdxZLmg58B/gW8GvgfcBq4I7sQjMzq217jBtGSxkuG6rvbgFJDwFvAD8GzouIjemsv0h6R5bBmZnVsvo6leV5BN0mAuB9EbGw2IyIOKGX4zEzs1R9LsfLb2zIvJ5STg19XNLI/IikUZL+PcOYzMwMWN/cwsr1VXD5KHBERLyRH4mIlcCR2YVkZmYAe+80EpH9RZqlJII6SQPzI5IGAQO7WN7MzHrBwIYcLW3ZP6KslDaCnwO3S7oiHT8duCq7kMzMDKAhp7LcWdxtIoiIb0h6FDg4nfRvEXFrtmGZmdmaphYANrW20VCXXY9ApRwREBG3ALdkFoWZmb3JHjsMg4eTbiayTASl9DV0gqRnJK2StFrSGkmrM4vIzMwA2huKn3g5211uKUcE3wSOiYgnM43EzMy2sPu4YQCsa27JtJ5SjjVedRIwMyu/0UMGALCpNdu7i0s5Ipgv6ZfA74B89xJExG8yi8rMzNq7on59XbZXDpVyRDAcWA8cBhyTvo7OMigzM4ORgxsAeOTFN7pZctuUcvno6ZlGYGZmRW0/rBGAwQOz7Yq6lKuGdpN0u6TH0vHpkr6UaVRmZgYkRwVZ90BayqmhHwHnA5sAIuJR4OQsgzIzs0Sdsu+KupREMDgiHugwLdtrmczMDIBcTrRF5RPBckm7AgEg6SRgaaZRmZkZAPU5ZX7VUCmXj/4jcDmwh6SXgOeBUzKNyszMANjY0kZzS7Y9kJZy1dBC4BBJQ4BcRKzJNCIzM2u308hBmdfRaSKQdEpE/EzSFzpMByAiLso4NjOzmpfLiYxvLO7yiGBI+ndYtiGYmVln6gRtGV811GkiiIgfpn+/kmkEZmbWqbpquGpI0lVFHl7/k0yjMjMzIDkdv2rDpkzrKOXy0elFHl4/M7uQzMwsb93GFl5d3ZRpHaUkgpykUfkRSaMp8clmZma2bcYMHciwxoZM6yhlh/4d4D5JvwIEnAR8LdOozMwMgDFDBrBo+bpM6yjlPoKfSpoPHJROOiEinsg0KjMzA5LG4qzbCLq6j2B4RKxOTwW9AvyiYN7oiHg908jMzIz1zS2VSwQkO/6jgYdI+xlKKR3fJcO4zMwMGDNkIOl9vJnpqrH46+nfqRGxS8FrckSUlAQkHS7paUnPSjqvk2XeL+kJSY9L+kWxZczMatXwQfXUZZwJukoEF6d/792agiXVAZcCRwDTgA9KmtZhmSkkzzp4R0TsCXxua+oyM+uvcsr+hrKuTg1tknQ5MF7S9zvOjIizuil7X+DZtNM6JF0LHAcUNjR/Arg0vTeBiHitJ8GbmfV3AjLuYaLLRHA0cAjwXpJ2gp7aCVhcML4E+LsOy+wGIOnPQB0wJyL+t2NBks4EzgTYeeedtyIUM7O+qaCjz/bh3tZVIjg3Ir4oaeeIuCqT2pP6pwAHAuOBuyXtXXgnM0BEXE7yTARmz56dcW40M6seufZEQGaNxl21ERypJP1s7fOJXwImFIyPT6cVWgLMjYhNEfE88DeSxGBmZkAu3fln2U7QVSL4X2AlMF3SaklrCv+WUPaDwBRJkyUNIEkoczss8zuSowEkjSU5VbSwpxthZtZfqT0RZFdHp4kgIs6NiJHATRExPCKGFf7truCIaAE+A9wKPAlcFxGPS/qqpGPTxW4FVkh6AriT5HTUim3eKjOzfiLfLpDlEUEpXUwcJ2kiMCUibpM0CKgv5ZGVEXEzcHOHaV8uGA7gC+nLzMw6WN2U3FW8asMmGhvqMqmjlOcRfAK4HvhhOmk8ySkdMzPL2OQxycMiK9VGkPePwDuA1QAR8QywfWYRmZlZu4q2ERTYGBHN+RFJ9WzZ95CZmWVEbL6PICulJIK7JP0LMEjSocCvgBsyi8jMzDZLjwiy7GWilERwHrAM+D/gkySNv1/KLiQzM8vLuONRoLSrhtqAH6UvMzMro8I7izOrI7uizcxsW6nCdxabmVmF5RNBllfoOBGYmVWxclw11G0bgaTdgHOBiYXLR8RBna5kZma9ohxHBN0mApLLRS8jaSxuzTAWMzPrQGVoLC4lEbRExH9nF4KZmXUmf/lopW8ou0HSP0jaQdLo/CuziMzMrF3+1FBrJdsIgI+mf88tmBbALr0fjpmZFdq4qQ2AZWs2sse4bOoo5YayydlUbWZm3dl5zGBg89VDWSjlqqEG4NPAu9NJ84AfRsSmzKIyMzMAGuqSM/ibWtsyq6OUU0P/DTQAP0jHT02nfTyroMzMLDFkQPIwmvwDarJQSiJ4e0TsUzB+h6QFWQVkZmabDR6Y7KabNmV39X4pVw21Sto1PyJpF3w/gZlZWdSllw1leGaopCOCc4A7JS0kuaR1InB6diGZmVleLv25XrHLRyXVAfsAU4Dd08lPR8TGzCIyM7N2+SOCtgyfVdnlqaGIaAU+GBEbI+LR9OUkYGZWJnW5/Kmhyt5Q9mdJlwC/BNblJ0bEXzOLyszMgM2JYOHytZnVUUoimJH+/WrBtADc+6iZWcaGNTZs8TcLpSSCMyJiYeGE9MohMzMrg0ENdZmeGirl8tHri0z7VW8HYmZmxdXnVJk7iyXtAewJjJB0QsGs4UBjZhGZmdkWJFj8+obMyu/q1NDuwNHASOCYgulrgE9kFpGZmW1hQH0djQ3ZPVm400QQEb8Hfi9p/4i4L7MIzMysSyMG1Vf84fWfkjQyPyJplKSfZBiTmZkVkFTxJ5RNj4g38iMRsRKYmVlEZma2hZyyfWZxKYkgJ2lUfiR9TGUpl52amVkvyEm0VfhRld8B7pOUv2T0fcDXMovIzMzeJMPbCEp6VOVPJc1n853EJ0TEE9mFZGZmhXJSxU8NAYwG1kXEJcAySX6OsZlZmUhUtrFY0oXAF4Hz00kNwM8yi8jMzLaQkyp++ej/A44l7Xk0Il4GhmUYk5mZFZDItLG4lETQHMkxSSQBaUhm0ZiZ2ZuIyj+z+DpJPwRGSvoEcBvwo8wiMjOzLWzY1MoLK9ZnVn4pVw19W9KhwGqS/oe+HBF/zCwiMzPbQmNDHQPqK9DXUKF0x++dv5lZBWw/bCAvv9GUWflddUO9Boo2VAuIiBieWVRmZtauYncWR4SvDDIzqwJZJ4KSTjpJeqek09Phsb6hzMysfOpyquyjKovcUDYA31BmZlY2uZx4fV1zduWXsIxvKDMzq6CV65pZ3dSSWfm+oczMrMrtNHIQgxvqMivfN5SZmVW5xoYcuZwyK983lJmZVTkpuyQAvqHMzKzmZXfPMiDpcElPS3pW0nldLHeipJA0O8t4zMz6qko/vH6rSKoDLgWOAKYBH5Q0rchyw4Czgb9kFYuZmXUuyyOCfeAkXWEAABckSURBVIFnI2JhRDQD1wLHFVnu34BvANl1pGFm1sdV+sE0SPpe4d8S7QQsLhhfkk4rLHcWMCEibuqm/jMlzZc0f9myZT0Iwcys78u4rbjkI4J3p38P6K2KJeWAi4B/6m7ZiLg8ImZHxOztttuut0IwM+s7quDh9VvjJWBCwfj4dFreMGAvYJ6kRcB+wFw3GJuZbUlke0iQZSJ4EJgiabKkAcDJwNz8zIhYFRFjI2JSREwC7geOjYj5GcZkZmYdZJYIIqIF+AxwK/AkcF1EPC7pq5KOzapeM7P+KMvG4pJuKNtaEXEzcHOHaV/uZNkDs4zFzKyvqpbG4l+kf3+eVSBmZta5it9QFhHfLvxrZmblk/EBQbZdTJiZWe+o+A1lZmZWORVtI5CUk/T32YZgZmbdybCJoOtEEBFtJB3HmZlZhWT9PIJSTg3dnnYTnXV7hZmZVUApieCTwK+AZkmrJa2RtDrjuMzMrEBk2FxcyqMqh2VWu5mZdSvr0zEl3Vks6QTgnSRXMP0pIn6XaVRmZraFijUWA0j6AfAp4P+Ax4BPSXIDsplZuWR8SFDKEcFBwNRI72+WdBXweKZRmZnZFip9Q9mzwM4F4xPSaWZmVgZZP4+glCOCYcCTkh4gSUr7AvMlzQWICHcpbWbWh3WaCCQNjIiNQNFuo83MrIwyPDfU1RHBfcAs4OMRcWp2IZiZWVeyvp23q0QwQNKHgL9PLx/dQkT8JruwzMysUKVuKPsU8GFgJHDMm2ICJwIzszKo2A1lEXEPcI+k+RHx44zjMDOzLlT0hjInATOzyqqWZxabmVkFVfqGMjMzq6CsbyjrUSKQNCejOMzMrEJ6ekTgu4jNzCogMmwt7mki8FPKzMzKrNoai9+WSRRmZtalqmksTh9mb2ZmZZT1qRhfNWRmVu0kIqCtLZvjAicCM7MqN2RAHQBrNrZkUv5WJQJJp/d2IGZmVlx9XbqrzqihYGuPCL7Sq1GYmVmnKtbpnKRHO5sFvCWbcMzMrDNZdUXdVTfUbwHeC6zsMF3AvZlEY2Zmb1LJB9PcCAyNiEc6zpA0L7OIzMysqKxuLu7qeQRndDHvQ9mEY2ZmHfk+AjMzA7K7u9iJwMysyinjRgInAjOzPiKrHkidCMzMqly19T5qZmYV4jYCM7Ma5auGzMwMyO4+AicCM7Nq56uGzMwMsutryInAzKzKuY3AzMwSbiMwM6tNvo/AzMwA30dgZlazlHErgROBmVkf0SfvI5B0uKSnJT0r6bwi878g6QlJj0q6XdLELOMxM+uL+mwbgaQ64FLgCGAa8EFJ0zos9jAwOyKmA9cD38wqHjOzvq4v3kewL/BsRCyMiGbgWuC4wgUi4s6IWJ+O3g+MzzAeM7M+qS/fR7ATsLhgfEk6rTNnALdkGI+ZWZ9W9mcWl5OkU4DZwAGdzD8TOBNg5513LmNkZmaV12fbCICXgAkF4+PTaVuQdAjwr8CxEbGxWEERcXlEzI6I2dttt10mwZqZVbu+eB/Bg8AUSZMlDQBOBuYWLiBpJvBDkiTwWoaxmJn1WX32PoKIaAE+A9wKPAlcFxGPS/qqpGPTxb4FDAV+JekRSXM7Kc7MrOZl9cziTNsIIuJm4OYO075cMHxIlvWbmfULfbiNwMzMelGfvLPYzMy2XV++j8DMzPoAJwIzsyonP7PYzMzAbQRmZjWroS45ItjU1pZJ+U4EZmZVrrGhDoCmTa2ZlO9EYGZW5epzyRFBa1vf64bazMx6kdsIzMxqVF/ufdTMzHpBvtO5vtj7qJmZ9Yb0iCCrTuecCMzMqlz+zJCPCMzMalT+zmI3FpuZ1ajNbcU+NWRmVpPU3kaQTflV8fD6bbVp0yaWLFlCU1NTpUPpdxobGxk/fjwNDQ2VDsWsZmV91VC/SARLlixh2LBhTJo0KfNe+mpJRLBixQqWLFnC5MmTKx2OWc3K+oigX5waampqYsyYMU4CvUwSY8aM8ZGWWZXw5aPdcBLIht9Xs8rz5aN9RF1dHTNmzGCfffZh1qxZ3HvvvQAsWrSIQYMGMXPmTKZOncq+++7LlVdeCcAVV1zBjBkzmDFjBgMGDGDvvfdmxowZnHfeeRXcEjOrOm4s7hsGDRrEI488AsCtt97K+eefz1133QXArrvuysMPPwzAwoULOeGEE4gITj/9dE4//XQAJk2axJ133snYsWMrswFmVrU2Nxb71FCfsXr1akaNGlV03i677MJFF13E97///TJHZWZ9lTI+N9Tvjgi+csPjPPHy6l4tc9qOw7nwmD27XGbDhg3MmDGDpqYmli5dyh133NHpsrNmzeKpp57q1RjNrP/Kuo2g3yWCSik8NXTffffxkY98hMcee6zoslm1/JtZ/5R1FxP9LhF098u9HPbff3+WL1/OsmXLis5/+OGHmTp1apmjMrO+qv0+ArcR9B1PPfUUra2tjBkz5k3zFi1axDnnnMNnP/vZCkRmZn1RLs0EGT2psv8dEVRKvo0AklM/V111FXV1yQOnn3vuOWbOnElTUxPDhg3jrLPO4rTTTqtgtGbWl+SfWdzS2pZN+ZmUWoNaW1uLTp80aRIbNmzodv1Fixb1ckRm1l/U16WJwA+vNzOrTfW5ZFfd0upEYGZWkzYfEWRzasiJwMysym1uI/ARgZlZTaqvS3bVrW4jMDOrTfkjgk0+NWRmVpv8YJoacOWVV/KZz3ym7PXOmzePo48+uuz1mll1cSIwM6txNZkIblp4E4ddfxjTr5rOYdcfxk0Lb9rmMo8//nje9ra3seeee3L55Ze3Tx86dCjnnnsue+65J4cccggPPPAABx54ILvssgtz585tX27x4sUceOCBTJkyha985StF6/j0pz/N7Nmz2XPPPbnwwgvbp5933nlMmzaN6dOnc84557xpvTlz5nDqqaey//77M2XKFH70ox+1z1u7di0nnXQSe+yxBx/+8IfdIZ5ZLYqIPvV629veFh098cQTb5rWmRufuzFmXz079rpyr/bX7Ktnx43P3VhyGcWsWLEiIiLWr18fe+65ZyxfvjwiIoC4+eabIyLi+OOPj0MPPTSam5vjkUceiX322SciIq644ooYN25cLF++vH39Bx98sNM6Wlpa4oADDogFCxbE8uXLY7fddou2traIiFi5cuWb1rvwwgtj+vTpsX79+li2bFmMHz8+Xnrppbjzzjtj+PDhsXjx4mhtbY399tsv/vSnP71p/Z68v2bW+15dvSEmfvHGuPq+RVtdBjA/Otmv1twRwcV/vZim1i0fxt7U2sTFf714m8r9/ve/zz777MN+++3H4sWLeeaZZwAYMGAAhx9+OAB77703BxxwAA0NDey9995bdCtx6KGHMmbMGAYNGsQJJ5zAPffc86Y6rrvuOmbNmsXMmTN5/PHHeeKJJxgxYgSNjY2cccYZ/OY3v2Hw4MFF4zvuuOMYNGgQY8eO5T3veQ8PPPAAAPvuuy/jx48nl8sxY8YMd3VhVoNqLhG8su6VHk0vxbx587jtttu47777WLBgQXsHcwANDQ3tfYnncjkGDhzYPtzS0tJeRseHxHccf/755/n2t7/N7bffzqOPPspRRx1FU1MT9fX1PPDAA5x00knceOON7Umno87Kz8cDyXOXC2Mys9pQc4lg3JBxPZpeilWrVjFq1CgGDx7MU089xf3339/jMv74xz/y+uuvs2HDBn73u9/xjne8Y4v5q1evZsiQIYwYMYJXX32VW265BUjO8a9atYojjzyS7373uyxYsKBo+b///e9pampixYoVzJs3j7e//e0931Az65dqrvfRs2edzZx752xxeqixrpGzZ5291WUefvjhXHbZZUydOpXdd9+d/fbbr8dl7Lvvvpx44oksWbKEU045hdmzZ28xf5999mHmzJnsscceTJgwoT1RrFmzhuOOO46mpiYigosuuqho+dOnT+c973kPy5cv54ILLmDHHXfkb3/7W8831sz6HUUfu0pk9uzZMX/+/C2mPfnkkz164tdNC2/i4r9ezCvrXmHckHGcPetsjtrlqN4OtWrMmTOHoUOHFr2iqBQ9fX/NrHe9tqaJfb92O/9+/F6cst/ErSpD0kMRMbvYvJo7IgA4apej+vWO38ysJ2oyEdSaOXPmVDoEM6tiNddYbGZmW3IiMDOrcU4EZmY1zonAzKzGORFUAXdDbWaVlGkikHS4pKclPSvpvCLzB0r6ZTr/L5ImZRmPmZm9WWaJQFIdcClwBDAN+KCkaR0WOwNYGRFvBb4LfCOreAqtuuEGnjnoYJ6cOo1nDjqYVTfcsM1lVnM31OvWreNjH/sY++67LzNnzuT3v//9Nm+vmfUjnXVLuq0vYH/g1oLx84HzOyxzK7B/OlwPLCe927mz17Z2Q/3G3Lnx5D4z4ond92h/PbnPjHhj7tySyyimmruhPv/88+Pqq69unz9lypRYu3Zt3HnnnXHUUUd1u23uhtqssvpyN9Q7AYsLxpek04ouExEtwCpgTMeCJJ0pab6k+cuWLdumoF777veIpi27oY6mJl777ve2qdxq7ob6D3/4A1//+teZMWMGBx54IE1NTbz44ovbtL1mVj4D6nK8a8pYdhjRmEn5feLO4oi4HLgckr6GtqWslqVLezS9FIXdUA8ePLh9Zwu93w31gw8+yKhRozjttNO26Ib69ttv5/rrr+eSSy7hjjvu2GLdiODXv/41u++++xbTX3311a3eZjMrn5GDB3D1GX+XWflZHhG8BEwoGB+fTiu6jKR6YASwIsOYqN9hhx5NL0W1d0P93ve+l//6r/9qfwzlww8/vBVbaWb9VZaJ4EFgiqTJkgYAJwNzOywzF/hoOnwScEfk91YZ2f7zn0ONWx5eqbGR7T//ua0u8/DDD6elpYWpU6dy3nnnbVM31NOnT+fEE0/sshvqD33oQ1t0Q3300Uczffp03vnOdxbthvqCCy5g06ZNTJ8+nT333JMLLrhg6zbUzPqlTLuhlnQk8D2gDvhJRHxN0ldJGi3mSmoErgZmAq8DJ0fEwq7K7I1uqFfdcAOvffd7tCxdSv0OO7D95z/HiGOO6dnG1RB3Q23W91WsG+qIuBm4ucO0LxcMNwHvyzKGYkYcc4x3/GZmKd9ZbGZW45wIzMxqXL9JBBm3Mdcsv69m/V+/SASNjY2sWLHCO61eFhGsWLGCxsZsbmIxs+rQJ24o68748eNZsmQJ23rXsb1ZY2Mj48ePr3QYZpahfpEIGhoamDx5cqXDMDPrk/rFqSEzM9t6TgRmZjXOicDMrMZl2sVEFiQtA17YytXHkjzzoJZ4m2uDt7k2bMs2T4yI7YrN6HOJYFtImt9ZXxv9lbe5Nniba0NW2+xTQ2ZmNc6JwMysxtVaIri8+0X6HW9zbfA214ZMtrmm2gjMzOzNau2IwMzMOuiXiUDS4ZKelvSspPOKzB8o6Zfp/L9ImlT+KHtXCdv8BUlPSHpU0u2SJlYizt7U3TYXLHeipJDU568wKWWbJb0//awfl/SLcsfY20r4395Z0p2SHk7/v4+sRJy9RdJPJL0m6bFO5kvS99P341FJs7a50ojoVy+Sx2I+B+wCDAAWANM6LPMPwGXp8MnALysddxm2+T3A4HT407Wwzelyw4C7gfuB2ZWOuwyf8xTgYWBUOr59peMuwzZfDnw6HZ4GLKp03Nu4ze8GZgGPdTL/SOAWQMB+wF+2tc7+eESwL/BsRCyMiGbgWuC4DsscB1yVDl8PHCxJZYyxt3W7zRFxZ0SsT0fvB/p6l6KlfM4A/wZ8A2gqZ3AZKWWbPwFcGhErASLitTLH2NtK2eYAhqfDI4CXyxhfr4uIu0me4d6Z44CfRuJ+YKSkHbalzv6YCHYCFheML0mnFV0mIlqAVcCYskSXjVK2udAZJL8o+rJutzk9ZJ4QETeVM7AMlfI57wbsJunPku6XdHjZostGKds8BzhF0hKSZ6R/tjyhVUxPv+/d6hfdUFvpJJ0CzAYOqHQsWZKUAy4CTqtwKOVWT3J66ECSo767Je0dEW9UNKpsfRC4MiK+I2l/4GpJe0VEW6UD6yv64xHBS8CEgvHx6bSiy0iqJzmcXFGW6LJRyjYj6RDgX4FjI2JjmWLLSnfbPAzYC5gnaRHJudS5fbzBuJTPeQkwNyI2RcTzwN9IEkNfVco2nwFcBxAR9wGNJH3y9Fclfd97oj8mggeBKZImSxpA0hg8t8Myc4GPpsMnAXdE2grTR3W7zZJmAj8kSQJ9/bwxdLPNEbEqIsZGxKSImETSLnJsRMyvTLi9opT/7d+RHA0gaSzJqaKF5Qyyl5WyzS8CBwNImkqSCPrz4wrnAh9Jrx7aD1gVEUu3pcB+d2ooIlokfQa4leSKg59ExOOSvgrMj4i5wI9JDh+fJWmUOblyEW+7Erf5W8BQ4Fdpu/iLEXFsxYLeRiVuc79S4jbfChwm6QmgFTg3Ivrs0W6J2/xPwI8kfZ6k4fi0vvzDTtI1JMl8bNrucSHQABARl5G0gxwJPAusB07f5jr78PtlZma9oD+eGjIzsx5wIjAzq3FOBGZmNc6JwMysxjkRmJnVOCcCKwtJIyX9Q6XjyJO0tofLHy9pWlbxFNRzlqQnJf087SX3NkmPSPqApP/pKgZJx3bVC6tZZ3z5qJVF2tX3jRGxVw/Xq4uI1gziWRsRQ3uw/JUk8V/f27F0qOcp4JCIWJLeLPTvEXFIlnWa+YjAyuXrwK7pr9tvSTpQ0t2Sbkr7mr8s7R8ISWslfUfSAmD/ra2wqzrS+V+TtCDtnO0t6bRJku7Q5uc27Czp74FjgW+l8e8qaUa63qOSfitpVLr+PEnfkPSApL9JelcnsZ0r6cF0/a+k0y4j6W75FklfBH4GvL2gznn5LjKU9NH/1zT+29Npp0m6JB3eTtKv0zoelPSOdPocJf3dz5O0UNJZBTF9JI1ngaSrJQ2T9LykhnT+8MJx60cq3fe2X7XxAiZR0L86yZ2TTSQ7vjrgj8BJ6bwA3t9JOecCjxR5fb/Ist3VcUw6/E3gS+nwDcBH0+GPAb9Lh6/Mr5uOPwockA5/FfheOjwP+E46fCRwW5G4DiPpQ18kP8ZuBN6dzlsEjC2I/8aC9eaRdBi4HUnvk5PT6aPTv6cBl6TDvwDemQ7vDDyZDs8B7gUGkvTHs4LkrtU9SfolGtuhzCuA49PhM/Pb5lf/evW7LiasT3kgIhZC+2317yR5PkQr8OtiK0TEt0i6y9jWOppJdsAADwGHpsP7Ayekw1eTJIktSBoBjIyIu9JJVwG/KljkNwXlTioS02Hp6+F0fChJx3B3l7hN+wF3R9KpHBFRrO/6Q4Bp2vyYjeGS8qfCboqk08GNkl4D3gIcBPwqIpZ3KPN/gH8m6cPodJLnHVg/40RgldSxgSo/3hSdtAtIOhf4cJFZd0fEWUWmd1bHpojID7fSu9+FfM+unZUr4D8j4oe9WGdHOWC/iNjigTxpYijsebbLbY+IP6enyw4E6iKi6OMTrW9zG4GVyxqSrqEL7aukV8kc8AHgnu4KiYhvRcSMIq9iSWBr6riXzZ0Qfhj4U8f4I2IVsLLg/P+pwF2U7lbgY/lf6JJ2krR9D9a/H3i3pMnp+qOLLPMHCh7QImlGN2XeAbxP0pgiZf6U5FTTFT2I0foQJwIri0h6wPyzpMck5U/tPAhcAjwJPA/8NoOqe1rHZ4HTJT1KsoM/O51+LXCukgek70rSjfm30uVmkLQTlCQi/kCyY71P0v+RnKrqmCS7Wn8Zyfn636QN6r8ssthZwOy08fcJ4FPdlPk48DXgrrTMiwpm/xwYBVxTaozWt/jyUauI9FTDORFxdF+uoxZIOgk4LiJOrXQslg23EZhZpyT9F3AEyRVQ1k/5iMDMrMa5jcDMrMY5EZiZ1TgnAjOzGudEYGZW45wIzMxqnBOBmVmN+/+GBzac0GeY6AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.plot(tpr, 1 - fpr, label='BDT')\n",
    "ax.set_xlabel('tpr = photon efficiency')\n",
    "ax.set_ylabel('1 - fpr = electron efficiency')\n",
    "ax.scatter(eff_ph_amb_as_ph, eff_el_amb_as_ph, color='C2', label='amb as ph', marker='o')\n",
    "ax.scatter(eff_ph_amb_as_el, eff_el_amb_as_el, color='C3', label='amb as el', marker='o')\n",
    "ax.legend()\n",
    "ax.set_title('only when both particles are reconstructed')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv3",
   "language": "python",
   "name": "venv3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}