{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  simpleClassifier.py\n",
    "#  G. Cowan / RHUL Physics / November 2021\n",
    "#  Simple program to illustrate classification with scikit-learn\n",
    "\n",
    "import scipy as sp\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib.ticker as ticker\n",
    "\n",
    "from sklearn.discriminant_analysis import LinearDiscriminantAnalysis as LDA\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import metrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#  read the data in from files, \n",
    "#  assign target values 1 for signal, 0 for background\n",
    "sigData = np.loadtxt('signal.txt')\n",
    "nSig = sigData.shape[0]\n",
    "sigTargets = np.ones(nSig)\n",
    "\n",
    "bkgData = np.loadtxt('background.txt')\n",
    "nBkg = bkgData.shape[0]\n",
    "bkgTargets = np.zeros(nBkg)\n",
    "\n",
    "# concatenate arrays into data X and targets y\n",
    "X = np.concatenate((sigData,bkgData),0)\n",
    "y = np.concatenate((sigTargets, bkgTargets))\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.5, random_state=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearDiscriminantAnalysis()"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Create classifier object and train\n",
    "# Add code here to include other claassifiers (MLP, BDT,...)\n",
    "clf = LDA()\n",
    "clf.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "power of test with respect to signal =  0.7916666666666666\n"
     ]
    }
   ],
   "source": [
    "# Evaluate accuracy using the test data.\n",
    "# If available, use the decision function, else (e.g. for MLP) use predict_proba\n",
    "# Adjust threshold value tCut or pMin as appropriate\n",
    "\n",
    "X_bkg_test = X_test[y_test==0]\n",
    "X_sig_test = X_test[y_test==1]\n",
    "y_bkg_test = y_test[y_test==0]\n",
    "y_sig_test = y_test[y_test==1]\n",
    "if hasattr(clf, \"decision_function\"):\n",
    "    tCut = 0.\n",
    "    y_bkg_pred = (clf.decision_function(X_bkg_test) >= tCut).astype(bool)\n",
    "    y_sig_pred = (clf.decision_function(X_sig_test) >= tCut).astype(bool)\n",
    "else:\n",
    "    pCut = 0.5\n",
    "    y_bkg_pred = (clf.predict_proba(X_bkg_test)[:,1]  >= pCut).astype(bool)\n",
    "    y_sig_pred = (clf.predict_proba(X_sig_test)[:,1]  >= pCut).astype(bool)\n",
    "\n",
    "power = metrics.accuracy_score(y_sig_test, y_sig_pred)        # = Prob(t >= tCut|sig)\n",
    "print('power of test with respect to signal = ', power)\n",
    "\n",
    "#  Add code here to obtain the background efficiency\n",
    "# = size of test alpha = Prob(t >= tCut|bkg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAG5CAYAAABFtNqvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEsklEQVR4nO3de3RU9b338c+E3ENmyI0IyFXgIFfjJSoCCU+WhCJawQu3crNoaQsWe9FKpRMKSg/Hx8NDacUFj8iBItISPc/yRiyQIHBQKCBe0AINCBElMZBBQmISfs8fOZlDTDKZ7CSTmcn7tdYszd6//du/vbMz82FfvmMzxhgBAACgyULaegAAAACBiiAFAABgEUEKAADAIoIUAACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWBTa1gMIZleuXNEXX3yh2NhY2Wy2th4OAADwgjFGFy9eVNeuXRUS4vmcE0GqFX3xxRfq3r17Ww8DAABYcPr0aV177bUe2xCkWlFsbKyk6l+E3W5v49EAAABvuFwude/e3f057glBqhXVXM6z2+0EKQAAAow3t+VwszkAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWBTa1gMAAPhOgUsqLvPcJj5S6mb3zXiAQEeQAoB2osAlZWyQLld6bhcVKm2fTpgCvEGQAoB2orisOkStyJT6xtff5nixtGBbdVuCFNA4ghQAtDN946Uhndt6FEBw4GZzAAAAiwhSAAAAFhGkAAAALCJIAQAAWESQAgAAsIggBQAAYBFBCgAAwCKCFAAAgEUEKQAAAIsIUgAAABb5bZDav3+/xo0bp7i4OMXExCg1NVWbNm3yevndu3frF7/4hW666SYlJCQoMjJSAwYM0BNPPKELFy7Uu0yvXr1ks9nqfc2dO7eFtgwAAAQLv/yuvdzcXGVmZio8PFyTJ0+Ww+FQdna2pk2bppMnT2rhwoWN9nH//ferqKhII0aM0IwZM2Sz2ZSbm6vly5dr69at2rt3rzp3rvtlUw6HQwsWLKgz/eabb26JTQMAAEHE74JUZWWl5syZI5vNpl27diklJUWS5HQ6dfvtt8vpdOqBBx5Qv379PPbz2GOPacaMGerSpYt7mjFGP/3pT/X8889r8eLF+uMf/1hnuU6dOikrK6tFtwkAAAQnv7u0t2PHDp04cUJTp051hyhJio2N1aJFi1RZWal169Y12s8TTzxRK0RJks1m06JFiyRJeXl5LTtwAADQ7vjdGanc3FxJ0pgxY+rMq5nWnBAUFhYmSQoNrX/Ty8vLtX79ehUUFCguLk7Dhw/XsGHDLK8PAAAEL78LUseOHZOkei/dxcXFKTEx0d3GihdffFFS/UFNkr788kvNmjWr1rSxY8dqw4YNSkxM9Nh3eXm5ysvL3T+7XC7L4wQAAP7P7y7tlZSUSKq+6bs+drvd3aapDh8+rMWLF6tz5856/PHH68x/6KGHlJubq8LCQrlcLu3bt0/f+9739Pbbb+uee+6RMcZj/8uWLZPD4XC/unfvbmmcAAAgMPhdkGot+fn5Gj9+vKqqqrR58+Z6zy799re/VVpamhITExUbG6tbb71Vr7/+ukaMGKH/+q//0ptvvulxHU8++aRKSkrcr9OnT7fW5gAAAD/gd0Gq5kxUQ2edXC5Xg2erGnLq1CmNHj1ahYWF+utf/6rRo0d7vWxISIhmz54tSdqzZ4/HthEREbLb7bVeAAAgePldkKq5N6q++6DOnz+voqKiRksfXO3kyZNKT0/XF198oS1btmj8+PFNHlPN2avS0tImLwsAAIKX3wWptLQ0SVJOTk6deTXTato0piZEFRQU6JVXXtH3v/99S2N67733JFVXPgcAAKjhd0EqIyNDffr00aZNm3T48GH39IsXL2rJkiUKDQ2t9VRdUVGRPv30UxUVFdXq5+oQtXnzZk2YMMHjej/55JN6vzpm9+7deu655xQREaGJEyc2Z9MAAECQ8bvyB6GhoVq7dq0yMzM1cuRITZkyRXa7XdnZ2crPz9fSpUvVv39/d/tVq1Zp8eLFcjqdtSqSp6en69SpU7rtttt05MgRHTlypM66rm6/ZcsWLV++XBkZGerVq5ciIiL00UcfKScnRyEhIVq9erV69OjRmpsOAAACjN8FKUkaPXq0du/eLafTqS1btujbb7/VoEGDtGTJEk2bNs2rPk6dOiVJ2rdvn/bt21dvm6uD1OjRo3X06FEdPHhQeXl5KisrU3JysiZNmqTHHntMqampzd4uAAAQXGymseJIsKzmCcOSkhKe4APQ5j48J41/WXp9ijSk7ne2e90GCHZN+fz2u3ukAAAAAoVfXtoDEEBKSqTGSoNER0tNrP8GAIGAIAXAupISadUqqaLCc7uwMGnePMIUgKBDkAJgXWlpdYiaOFFKSqq/TWGhlJ1d3ZYgBSDIEKQANF9SktSlS1uPAgB8jpvNAQAALCJIAQAAWESQAgAAsIggBQAAYBFBCgAAwCKCFAAAgEUEKQAAAIsIUgAAABYRpAAAACwiSAEAAFhEkAIAALCIIAUAAGARQQoAAMAighQAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARQQpAAAAi0LbegAAWkFJiVRa6rlNdLTkcPhmPAAQpAhSQLApKZFWrZIqKjy3CwuT5s0jTAFAMxCkgGBTWlodoiZOlJKS6m9TWChlZ1e3JUgBgGUEKSBYJSVJXbq09SgAIKhxszkAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYxFN7QHtWWOh5PkU7AcAjghTQHkVHVxfkzM723I6inQDgEUEKaI8cjuqA5OlrZCjaCQCNIkgB7ZXDQUACgGbiZnMAAACLOCMFAEGiwCUVlzU8/3ix78YCtBcEKQAIAgUuKWODdLnSc7uoUCk+0jdjAtoDghQABIHisuoQtSJT6hvfcLv4SKmb3XfjAoIdQQoAgkjfeGlI57YeBdB+cLM5AACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWESQAgAAsMhvg9T+/fs1btw4xcXFKSYmRqmpqdq0aZPXy+/evVu/+MUvdNNNNykhIUGRkZEaMGCAnnjiCV24cKHV1gsAANoPv6wjlZubq8zMTIWHh2vy5MlyOBzKzs7WtGnTdPLkSS1cuLDRPu6//34VFRVpxIgRmjFjhmw2m3Jzc7V8+XJt3bpVe/fuVefOtYuttMR6AQBA++F3QaqyslJz5syRzWbTrl27lJKSIklyOp26/fbb5XQ69cADD6hfv34e+3nsscc0Y8YMdenSxT3NGKOf/vSnev7557V48WL98Y9/bPH1AgCA9sPvLu3t2LFDJ06c0NSpU91hRpJiY2O1aNEiVVZWat26dY3288QTT9QKUZJks9m0aNEiSVJeXl6rrBdAGyspkc6e9fwqKWnrUQIIEn53Rio3N1eSNGbMmDrzaqZ9NwQ1RVhYmCQpNLT2prf2egH4QEmJtGqVVFHhuV1YmDRvnuRw+GZcAIKW3wWpY8eOSVK9l9Di4uKUmJjobmPFiy++KKluYGqJ9ZaXl6u8vNz9s8vlsjxOABaUllaHqIkTpaSk+tsUFkrZ2dVtCVIAmsnvglTJf59ydzTwBme323XmzBlLfR8+fFiLFy9W586d9fjjj7f4epctW6bFixdbGhuAFpSUJH3n0j4AtAa/u0eqteTn52v8+PGqqqrS5s2blZiY2OLrePLJJ1VSUuJ+nT59usXXAQAA/IffnZGqOSNU0sDNoC6Xq8GzRg05deqURo8ercLCQm3dulWjR49ulfVGREQoIiKiSWMDAACBy+/OSNXco1Tf/Ujnz59XUVFRk0oQnDx5Uunp6friiy+0ZcsWjR8/3ifrBQAAwc/vglRaWpokKScnp868mmk1bRpTE6IKCgr0yiuv6Pvf/75P1gsAANoHvwtSGRkZ6tOnjzZt2qTDhw+7p1+8eFFLlixRaGioZs2a5Z5eVFSkTz/9VEVFRbX6uTpEbd68WRMmTGjR9QIAAPjdPVKhoaFau3atMjMzNXLkSE2ZMkV2u13Z2dnKz8/X0qVL1b9/f3f7VatWafHixXI6ncrKynJPT09P16lTp3TbbbfpyJEjOnLkSJ11Xd2+qesF2kxJSfWj+w0pLPTdWACgnfO7ICVJo0eP1u7du+V0OrVlyxZ9++23GjRokJYsWaJp06Z51cepU6ckSfv27dO+ffvqbXN1kGqp9QKtqikFJ6OjfTMmAGjH/DJISVJqaqreeuutRttlZWXVCURS9ffqteZ6gTbhTcFJqTpEUWwSAFqd3wYpAB5QcBIA/ILf3WwOAAAQKAhSAAAAFhGkAAAALCJIAQAAWESQAgAAsIin9gD4RmOFQinZACAAEaQAtK7o6OoCodnZntuFhUnz5hGmAAQUghSA1uVwVAekxr7WJju7ug1BCkAAIUgBaH0OBwEJQFAiSAEAWkWBSyou89wmPlLqZvfNeIDWQJACALS4ApeUsUG6XOm5XVSotH06YQqBiyAFAGhxxWXVIWpFptQ3vv42x4ulBduq2xKkEKgIUgCAVtM3XhrSua1HAbQeCnICAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWESQAgAAsIiviAH8SUmJVFra8PzCQt+NBQDQKIIU4C9KSqRVq6SKCs/twsKk6GjfjAkA4BFBCvAXpaXVIWriRCkpqeF20dGSw+G7cQEAGkSQAvxNUpLUpUtbjwLt3PFi3y3fWNv4SKmbvXnjAVoLQQoA4BYfKUWFSgu2Nb+vqNDq/pq7rqhQaft0whT8E0EKAODWzV4dWorLmt9XY2eSvFnX8eLqoFVcRpCCfyJIAQBq6Wb3XWjx5bqA1kAdKQAAAIsIUgAAABYRpAAAACwiSAEAAFhEkAIAALCIp/YAoI0VuJpfbqC5BTTbA2/3MwVA0RQEKQBoQwUuKWODdLmy+X01VgCzPWvKfqYAKJqCIAUAbai4rPrDfUWm1De+eX1xJqVh3u5nCoCiqQhSAOAH+sZLQzq39SiCH/sZLY2bzQEAACwiSAEAAFhEkAIAALCIIAUAAGARQQoAAMAighQAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgkd8Gqf3792vcuHGKi4tTTEyMUlNTtWnTJq+XP3funJYtW6b7779fvXv3ls1mk81m87hMr1693O2++5o7d25zNwkAAASZ0LYeQH1yc3OVmZmp8PBwTZ48WQ6HQ9nZ2Zo2bZpOnjyphQsXNtrHJ598ooULF8pms6lfv36Kjo5WaWlpo8s5HA4tWLCgzvSbb77ZyqYAAIAg5ndBqrKyUnPmzJHNZtOuXbuUkpIiSXI6nbr99tvldDr1wAMPqF+/fh77uf7665WXl6eUlBTFxsZqwIAB+uyzzxpdf6dOnZSVldUSmwIAAIKc313a27Fjh06cOKGpU6e6Q5QkxcbGatGiRaqsrNS6desa7Sc5OVmjRo1SbGxsaw4XAAC0Y353Rio3N1eSNGbMmDrzaqbl5eW12vrLy8u1fv16FRQUKC4uTsOHD9ewYcNabX0AACBw+V2QOnbsmCTVe+kuLi5OiYmJ7jat4csvv9SsWbNqTRs7dqw2bNigxMREj8uWl5ervLzc/bPL5WqNIQIAAD/hd5f2SkpKJFXf9F0fu93ubtPSHnroIeXm5qqwsFAul0v79u3T9773Pb399tu65557ZIzxuPyyZcvkcDjcr+7du7fKOAEAgH/wuzNSbem3v/1trZ9vvfVWvf7660pLS9Pu3bv15ptv6q677mpw+SeffFI///nP3T+7XC7CFBDMSkqkxp4Gjo6WGviHIYDA16wgtX37du3YsUN79+7VmTNnVFRUpOjoaCUlJWnIkCFKS0vT+PHjdc0113jdZ82ZqIbOOrlcrgbPVrWGkJAQzZ49W7t379aePXs8BqmIiAhFRET4bGwA2lBJibRqlVRR4bldWJg0bx5hCghSTQ5S33zzjVauXKk1a9bo888/d1/uioyMVHx8vC5fvqyPPvpIR44c0Z///GeFhobqnnvu0WOPPaY77rij0f5r7o06duyYbrrpplrzzp8/r6KiIg0fPrypw26WmnujvKlDBaCdKC2tDlETJ0pJSfW3KSyUsrOr2xKkgKDUpHukVq9erb59++qpp55Sp06dtHTpUu3YsUMul0ulpaU6c+aMvv76a1VUVOjTTz/V+vXrNWnSJOXk5GjUqFGaOHGi8vPzPa4jLS1NkpSTk1NnXs20mja+8t5770mqrnwOALUkJUldutT/aihgAQgaTQpS8+fP19ixY/Xhhx/q0KFDevLJJ5Wenq6OHTvWamez2dS/f39Nnz5dGzZs0FdffaU1a9boww8/1IYNGzyuIyMjQ3369NGmTZt0+PBh9/SLFy9qyZIlCg0NrfVUXVFRkT799FMVFRU1ZVPq+OSTT3ThwoU603fv3q3nnntOERERmjhxYrPWAaD9cZVLR4ukD8/V/zpe3NYjBNAcTbq09+mnn+q6665r8kqioqL00EMPaebMmTpz5oznAYWGau3atcrMzNTIkSM1ZcoU2e12ZWdnKz8/X0uXLlX//v3d7VetWqXFixfL6XTWqUh+deA6e/ZsnWnPPvus+7Ldli1btHz5cmVkZKhXr16KiIjQRx99pJycHIWEhGj16tXq0aNHk7cdQPt19qK09QNpfZR0zsOVvahQKT7Sd+MC0HKaFKS+G6I+//zzJoWLDh06qGfPno22Gz16tHbv3i2n06ktW7bo22+/1aBBg7RkyRJNmzbN6/WtX7/e47SsrCx3kBo9erSOHj2qgwcPKi8vT2VlZUpOTtakSZP02GOPKTU11ev1AoAkXSiXKq5Ii9Ol7gMabhcfKXWz+2pUAFpSs57a6927t373u9/pN7/5TUuNxy01NVVvvfVWo+2ysrIa/G68xuo+XS0tLc3n914BaB96x0nXd27rUQBoDc0qyGmMUVVVlcc2e/bs0datW5uzGgAAAL/U5DNSBw4c0JAhQ7yul/S3v/1Nv/vd7xoNXACgwkLP8wO0uGVoUaF01kODAN0uABaCVGpqqkJDQzVgQPUF//3792v79u0aNmxYvd9FV15ertBQCqgD8CA6urpwZXa253YBVtzySlS0KjuEyf5mtuTpu9YDbLsA/I8mJ5xf/vKXOnTokA4ePCibzaY33nhDb775piSpa9euSklJUUpKioYNG6aOHTvqr3/9K1+TAsAzh6M6SHgqehuAxS2v2B1aP3Kext9VquSGvvM8ALcLwP9ocpBavny5+/9DQkL04IMP6oYbbtChQ4d0+PBhvfHGG3r99ddls9kkVd9H9fvf/77lRgwgODkcjQYJV7lUUCRVdqh/fmiR1K1c8qcH4C5GOVSZ7JC42RwISs265jZv3jwNHz5ckydPdk+7dOmSPvjgAx05ckRff/21UlJSNG7cuGYPFED75k1Nps4l0swPpPsuVhcWB4DW1qwgtXLlyjrTYmJiNHz4cJ9/Hx6A4OZNTabTn0qf7a1uS44C4AvcBQ4goHiqyRT6lfSZb4cDoJ1rVh0pAACA9qxJQWr8+PH6+9//bmlFly9f1rPPPqvnn3/e0vIAAAD+pkmX9k6fPq3U1FSlp6dr+vTpmjhxoux2z8/HHDhwQBs3btSmTZv0zTff1Pv9dwDgc0Fa/BOAbzUpSK1bt04HDhzQM888o4ceekhz5szRgAEDdOONNyo5OVlxcXG6fPmyiouLdezYMR04cEAlJSXuMglPP/20evXq1UqbAgBeCNLinwDaRpOC1C233CKn06l//vOfeuONN7R+/Xrl5uZq48aNddqGhIRo6NChuvfeezVnzhx17dq1xQYNAJYFafFPAG2jSUHKZrPJGKOQkBDdfffdmj17tpYtW6YRI0bozJkz+vrrrxUVFaWkpCQNGjRIDt6AAPijIC3+CcD3mhSkOnfurOPHj7t/Li4u1tmzZ3X99dfr+uuvb/HBAUBboPgnAG81KUjdeeedevnll9W5c2fde++9rTQkAGhbFP8E4K0mBanly5fr448/1r//+79rxYoVstlsev755/XBBx8oJSVFN9xwg2644QZde+21rTVeAPAZin8CaEyTglRycrIOHDig3Nxcbd++XU8//bS++eYbvfbaa3r11VfdX1QcHx+vG264oVa4GjhwYKtsAAAAQFux9BUx6enpSk9P19NPP61f/epXeuyxx3T48GEdOnTI/d9du3Zp+/btkqpvUq+qqmrRgQMAALS1Zn3X3osvvqjrrrtOsbGxGjlypEaOHOmeV1FRoY8++kgHDx7U4cOHmztOwL+VlHh+nF6iwCMABKFmBalZs2Y1OC8sLEwpKSlKSUlpzioA/1dSIq1aJVVUeG5HgUcACDrNClIAVH0mqqJCmjhRSkqqvw0FHgEgKBGkgJaSlERBIT+Rf16qPNfw/PhIqRuVNAPK8WLP8/mdoq0QpAAEjU4RUliI5MyVzh1quF1UqLR9Oh+8gSA+svr3tWCb53b8TtFWCFIAgkaXWGnGMCljglSZXH+b48XVH8rFZXzoBoJu9uqAVFzWcBt+p2hLBCkAQcUeIdkTJTVQSBOBp5udgAT/FdLWAwAAAAhUBCkAAACLuLQHIKCEFhVKZxuYWVjo07EAAEEKQEC4EhWtyg5hsr+ZLeV5aBgWVl1FHgB8gCAFICBcsTu0fuQ8jb+rVMmJHhryVTwAfIggBSBgXIxyqDLZwRN5APwGN5sDAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARTy1B/iSp4KR/lpMMhDHHIga25eUdQD8EkEK8IXo6OpCkdnZntv5UzHJQBxzIGrKfp43jzAF+BmCFOALDkf1h2Bpqed2/nTWIRDHHIi82c+FhdVBq7SUfQ34GYIU4CsOR+B9CAbimAMR+xkIWNxsDgAAYBFBCgAAwCKCFAAAgEUEKQAAAIsIUgAAABbx1B4CU0mJ7x7Lb2xdFKUEgHaLIIXAU1IirVolVVR4btcSBQybsi6KUgJAu0OQQuApLa0ONhMnSklJ9bdpqQKG3qxLoiglALRTBCkErqQkqUuX4FsXACBgcLM5AACARX4bpPbv369x48YpLi5OMTExSk1N1aZNm7xe/ty5c1q2bJnuv/9+9e7dWzabTTabrdXXCwAA2g+/vLSXm5urzMxMhYeHa/LkyXI4HMrOzta0adN08uRJLVy4sNE+PvnkEy1cuFA2m039+vVTdHS0Sht5yqsl1gsAANoPvzsjVVlZqTlz5shms2nXrl1as2aNnn32WX3wwQcaNGiQnE6njh071mg/119/vfLy8lRSUqLPPvtM3bt398l6AQBA++F3QWrHjh06ceKEpk6dqpSUFPf02NhYLVq0SJWVlVq3bl2j/SQnJ2vUqFGKjY316XoBAED74XdBKjc3V5I0ZsyYOvNqpuXl5QXNegEAQODyu3ukai6f9evXr868uLg4JSYmtsoltpZYb3l5ucrLy90/u1yulh0kAADwK34XpEpKSiRJjgaKG9rtdp05c8Yv17ts2TItXry4xccGBLoCl1Rc5rlNfKTUze6b8bSk/PNS5bn65x0v9u1Y2jtP+5vfBVqL3wWpQPbkk0/q5z//uftnl8vV6E3uQLArcEkZG6TLlZ7bRYVK26cHTpjqFCGFhUjOXOncoYbbRYVWh0S0nvjI6v28YJvndvwu0Br8LkjVnBGqOUP0XS6Xq8GzRm293oiICEVERLT42IBAVlxWHaJWZEp94+tvc7y4+kOwuCxwglSXWGnGMCljglSZ3HC7QD3TFki62atDeLCe9YR/87sgVXOP0rFjx3TTTTfVmnf+/HkVFRVp+PDhQbNeoL3oGy8N6dzWo2hZ9gjJnigpyLYrEHWzE5LQNvzuqb20tDRJUk5OTp15NdNq2gTDegEAQODyuyCVkZGhPn36aNOmTTp8+LB7+sWLF7VkyRKFhoZq1qxZ7ulFRUX69NNPVVRU5NP1AgAA+N2lvdDQUK1du1aZmZkaOXKkpkyZIrvdruzsbOXn52vp0qXq37+/u/2qVau0ePFiOZ1OZWVl1err6uBz9uzZOtOeffZZJSYmWlovAACA3wUpSRo9erR2794tp9OpLVu26Ntvv9WgQYO0ZMkSTZs2zet+1q9f73FaVlaWO0i15HoBAED74JdBSpJSU1P11ltvNdouKyurzpmoGsaYVlsvAACA3wYpAGhNngo0nj7vu3HAPzVWwJNSCqhBkALQrnhTvLFziTQzpLroJtqXphT3DKQCsmg9BCkA7Yo3xRtDv5K6XZbssb4bF/yDN8dHIBaQReshSAFodxot3lglibNR7RbFPdEUfldHCgAAIFAQpAAAACwiSAEAAFjEPVJo30pKpNLShucXFvpuLACAgEOQQvtVUiKtWiVVVHhuFxYmRUf7ZkwAgIBCkEL7VVpaHaImTpSSkhpuFx0tORy+GxcAIGAQpICkJKlLl7YeBQAgAHGzOQAAgEUEKQAAAIsIUgAAABYRpAAAACwiSAEAAFjEU3sA0BBPBVkp1gpABCkAqCs6uroQa3a253YUawXaPYIUAHyXwyHNm+f564MkirUCIEgBQL0cDkISgEZxszkAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWESQAgAAsIggBQAAYFFoWw8AQGArcEnFZQ3PP17sfV+e2jalHwDwFYIUAMsKXFLGBulyped2UaFSfGTD8+Mjq9ss2Na8fgDA1whSACwrLqsOUSsypb7xDbeLj5S62Rue380ubZ/u+cyWN/0AgK8RpAA0W994aUjn5vXRzU5IAhB4uNkcAADAIoIUAACARQQpAAAAiwhSAAAAFnGzOQAEisJCz/OjoyWHwzdjASCJIAUA/i86WgoLk7KzPbcLC5PmzSNM+YnGitVKlPQIBgQpAPB3Dkd1QCotbbhNYWF10CotJUj5gaYUq90+nTAVyAhSABAIHA4CUgDxpljt8eLqav7FZQSpQEaQAgCglbREsVr4N57aAwAAsIggBQAAYBFBCgAAwCKCFAAAgEV+G6T279+vcePGKS4uTjExMUpNTdWmTZua1MeVK1e0atUqDR06VFFRUUpKStKDDz6oY8eO1du+V69estls9b7mzp3bEpsFXysslM6erf/VWHFDIBB5OubPnpVKStp6hEBQ8cun9nJzc5WZmanw8HBNnjxZDodD2dnZmjZtmk6ePKmFCxd61c/cuXO1Zs0aDRw4UPPnz9dXX32lV155RTk5Odq7d68GDhxYZxmHw6EFCxbUmX7zzTc3d7PgS00pYBgd7Zsx+RmKBQYZinb63PFia/MQXPwuSFVWVmrOnDmy2WzatWuXUlJSJElOp1O33367nE6nHnjgAfXr189jPzt37tSaNWs0cuRIvfPOO4qIiJAkzZgxQ3feead+/OMfKy8vr85ynTp1UlZWVotvF3zMmwKGUrv9Sg2KBQYhinb6THxk9d/Ggm2e20WFVrdFcPO7ILVjxw6dOHFCs2fPdocoSYqNjdWiRYs0efJkrVu3Ts8884zHftasWSNJWrp0qTtESVJGRoYyMzP19ttv6x//+If69+/fOhuCtkcBwwZRLDBIccz7RDd79T8wOKMLyQ+DVG5uriRpzJgxdebVTKvvTFJ9/cTExOiOO+6oM68mSOXl5dUJUuXl5Vq/fr0KCgoUFxen4cOHa9iwYRa2BPB/FAsErOlmJyShmt8FqZobweu7dBcXF6fExMQGbxavcenSJZ09e1aDBw9Whw4d6syv6bu+fr788kvNmjWr1rSxY8dqw4YNSkxM9Lje8vJylZeXu392uVwe2wMAgMDmd0/tlfz3EyWOBk5P2+12d5vm9HF1uxoPPfSQcnNzVVhYKJfLpX379ul73/ue3n77bd1zzz0yxnhc77Jly+RwONyv7t27e2wPAAACm98Fqbb029/+VmlpaUpMTFRsbKxuvfVWvf766xoxYoT+67/+S2+++abH5Z988kmVlJS4X6dPn/bRyAEAQFvwuyBVcxapobNOLperwTNNTenj6naehISEaPbs2ZKkPXv2eGwbEREhu91e6wUAAIKX390jdfX9SzfddFOteefPn1dRUZGGDx/usY+YmBh16dJF+fn5qqqqqnOflKf7sOpTc29UaWOP0gNAIGisGG07LQsCWOF3QSotLU3Lli1TTk6OJk+eXGteTk6Ou403/WzevFl79uzRqFGjas3btm2b1/1I0nvvvSepuvI5AAQsinYCLc7vglRGRob69OmjTZs26dFHH9UNN9wgSbp48aKWLFmi0NDQWk/VFRUVqaioSImJibWeqnvkkUe0efNmPfXUU/rb3/6m8PBwSdL27du1bds2jRo1qlbpg08++URdu3ZVp06dao1n9+7deu655xQREaGJEye22nYDQKujaCfQ4vwuSIWGhmrt2rXKzMzUyJEjNWXKFNntdmVnZys/P19Lly6tFYBWrVqlxYsXy+l01qpIPnr0aM2ZM0dr165VSkqK7rrrLvdXxNjtdj3//PO11rtlyxYtX75cGRkZ6tWrlyIiIvTRRx8pJydHISEhWr16tXr06OGr3QAArYOinUCL8rsgJVWHoN27d8vpdGrLli369ttvNWjQIC1ZskTTpk3zup8XXnhBQ4cO1QsvvKCVK1eqY8eOuvvuu/X000/XKcQ5evRoHT16VAcPHlReXp7KysqUnJysSZMm6bHHHlNqampLbyYAAAhwfhmkJCk1NVVvvfVWo+2ysrIa/G68kJAQzZ8/X/Pnz2+0n7S0NK/vmQIAAJD8sPwBAABAoCBIAQAAWESQAgAAsMhv75EC4B+OF1ubB8A7jf0dxUdK3fiiDL9FkAJQr/hIKSpUWrDNc7uo0Oq2AJqmKX9j26cTpvwVQQpAvbrZq9+8i8s8t+Nfy4A13vyNHS+uDlrFZfyd+SuCFIAGdbPz5g20Jv7GAh83mwMAAFhEkAIAALCIIAUAAGARQQoAAMAibjaH/ykpkUpLG55fWOi7sQAA4AFBCv6lpERatUqqqPDcLixMio72zZgAAGgAQQr+pbS0OkRNnCglJTXcLjpacjh8Ny4AAOpBkIJ/SkqSunRp61EAAOARN5sDAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARTy1B9+i2CYAIIgQpOA7FNtsVIFLKi7z3CY+Uupm9814AAQG3jvaDkEKvkOxTY8KXFLGBulyped2UaHS9um8IQKoxntH2yJIwfcotlmv4rLqN8IVmVLf+PrbHC+WFmyrbsubIQCJ9462RpAC/EzfeGlI57YeBYBAw3tH2+CpPQAAAIsIUgAAABYRpAAAACwiSAEAAFjEzebwTmOFNKV2W7YACEqNFcfl7x2QRJCCN5pSSHPePN5cgUAWHV39t5yd7bkdf++AJIIUvOFNIc3Cwuo33tJS3liBQOZwVAekxr7Kib93QBJBCk1BIU2gfXA4CEiAl7jZHAAAwCKCFAAAgEUEKQAAAIsIUgAAABYRpAAAACziqT20LE9F/Bor8AcAQIAhSKFlNKWIX3S0b8YUxI4Xe54fHyl1s/tmLABan6e/+cbeD9C6CFJoGd4U8ZP4Wolmio+UokKlBds8t4sKlbZPJ0wBga4pf/Pxkb4ZE2ojSKHlUMSv1XWzVwek4rKG2xwvrn7TLS4jSAGBzpu/eYmz0G2JIAUEmG523jCB9oS/ef/GU3sAAAAWEaQAAAAsIkgBAABYRJACAACwiJvNA1lJScuUG2isHz8spFngapmnWLzpx1eoBYOA0xLvDZRE8amWeJ9pqffWYHnSkCAVqEpKpFWrpIoKz+3CwqrrOzX0RtWUfvykkGaBS8rYIF2u9NyusVpK3vbjS9SCQUDwtgCvNxp7j0KL8LYelTda6r01WOrdEaQCVWlpdfiZOFFKSqq/TWFh9RtdaWnDb1Le9CP51b8ai8uq/0BXZEp94+tv400tJW/68bVg+Rcagpy3BXgb4817FFqEt/WoGtNS763BVO+OIBXokpKkLl38px8f6hsvDensP/0A7QoFeAOOr+tRtZf3Vr+92Xz//v0aN26c4uLiFBMTo9TUVG3atKlJfVy5ckWrVq3S0KFDFRUVpaSkJD344IM6duxYq64XAAC0D355Rio3N1eZmZkKDw/X5MmT5XA4lJ2drWnTpunkyZNauHChV/3MnTtXa9as0cCBAzV//nx99dVXeuWVV5STk6O9e/dq4MCBrbJeAADQPvhdkKqsrNScOXNks9m0a9cupaSkSJKcTqduv/12OZ1OPfDAA+rXr5/Hfnbu3Kk1a9Zo5MiReueddxQRESFJmjFjhu688079+Mc/Vl5eXouvFwAAtB9+d2lvx44dOnHihKZOneoOM5IUGxurRYsWqbKyUuvWrWu0nzVr1kiSli5d6g5RkpSRkaHMzEzt2rVL//jHP1p8vQAAoP3wuyCVm5srSRozZkydeTXTrj6T5KmfmJgY3XHHHXXmZWZm1umnpdYLAADaD7+7tFdzI3h9l9Di4uKUmJjo8WZxSbp06ZLOnj2rwYMHq0OHDnXm1/R9dT8tsd7y8nKVl5e7fy4pKZEkuVwuj8tZcvGivrlYroKTF1VZGlNvk9BzFxVXXK6Lf89XZcLF+tt8XaTY4nKd99CPv/nneelKmfTNRcnVQM2lby5WtznyefX/W+0nELXnbUeAuXhRKi+X8vOr/x9+r+xrKeEr6bO/S2Wd6m9z6kJ1m7ITkuuC9X68FZ8Yo6Tk2OZ18h01n9vGmMYbGz9z5513Gknm2LFj9c7v06ePCQ8P99hHQUGBkWTuuOOOeufv2rXLSDKPPPJIi67X6XQaSbx48eLFixevIHidPn3a4+e+Mcb43RmpQPbkk0/q5z//ufvnK1euqLi4WAkJCbLZbHXau1wude/eXadPn5bdHuAVyVoZ+8p77Cvvsa+8x75qGvaX9/xxXxljdPHiRXXt2rXRtn4XpBz/XeCt5rLYd7lcLneb5vRxdbuWWm9EREStG9slqVOnTh6XkSS73e43B4+/Y195j33lPfaV99hXTcP+8p6/7avGPvNr+N3N5vXdv1Tj/PnzKioqarQEQUxMjLp06aL8/HxVVVXVmV/f/VAtsV4AANC++F2QSktLkyTl5OTUmVczraZNY/1cunRJe/bsqTNv27ZtdfppqfUCAIB2pNG7qHysoqLC9OnTx0RERJhDhw65p7tcLjNo0CATGhpqPvvsM/f0wsJCc/ToUVNYWFirnx07dhhJZuTIkaa8vNw9/W9/+5ux2Wxm1KhRzVpvSygrKzNOp9OUlZW1aL/BiH3lPfaV99hX3mNfNQ37y3uBvq9sxnjzbJ9v7dy5U5mZmYqIiNCUKVNkt9uVnZ2t/Px8LV26VL/5zW/cbbOysrR48WI5nU5lZWXV6ufhhx/W2rVrNXDgQN11113ur4iJjIys9ytimrJeAAAAv7u0J0mjR4/W7t27NWLECG3ZskV/+tOflJCQoI0bNzYpzLzwwgtauXKlbDabVq5cqTfeeEN333233n///TohqiXXCwAA2ge/PCMFAAAQCPzyjBQAAEAgIEgBAABYRJBqZSdPnpTNZvP4qu/7ABviqZ/f//73rbglvjFr1qwGt2/AgAFN7m/btm1KT0+X3W5XbGys0tPT3eUvAt2lS5e0ceNGPfjgg+rfv7+ioqLUqVMnpaWl6eWXX25yf8FybO3fv1/jxo1TXFycYmJilJqaqk2bNjWpjytXrmjVqlUaOnSooqKilJSUpAcffLDR79sMFAUFBVqxYoXGjBmjHj16KDw8XNdcc43uu+8+vffee173k5ub6/G42bdvXytuhW/16tWrwe2cO3eu1/0E+7H10ksvNfqZl5GR0Wg/gXRs+V1l82DTqVMnOZ3OeucdOHBAb7zxhjIzM5vUZ8+ePTVr1qw600eMGGFliH7pZz/7WZ2q8ImJiU3q489//rN+8IMfKDExUTNnzpTNZtOWLVs0duxYbdy4UdOmTWvBEfveu+++q+nTpyshIUEZGRm67777dO7cOWVnZ2vq1Knau3ev/vCHPzSpz0A/tnJzc5WZmanw8HBNnjxZDodD2dnZmjZtmk6ePKmFCxd61c/cuXO1Zs0aDRw4UPPnz3c/8ZuTk1PvE7+B5g9/+IP+9V//Vdddd53uvPNOde7cWceOHdNrr72m1157TS+//LIefPBBr/tLS0tTenp6nenXXnttC4667TkcDi1YsKDO9JtvvtnrPoL92Lrhhhsa/Mz761//qo8//rhJn3kBcWy1bfWF9m38+PFGktm6davXy0gyaWlprTeoNjZz5kwjyeTn5zern+LiYtOpUyeTmJhoPv/8c/f0L774wlxzzTWmU6dOpri4uJmjbVuHDx82f/7zn823335ba/qXX35pevbsaSSZ999/3+v+Av3YqqioMNddd52JiIgwBw8edE+/uhbcP/7xj0b7uboG3dV1bRqqQReItm7danbt2lVn+q5du0xYWJiJj4/3qqbPzp07jSTjdDpbYZT+pWfPnqZnz57N6qM9HFsNKS8vNwkJCSY0NNR8+eWXjbYPpGOLS3tt5IsvvtBbb72lzp076+67727r4QSdv/zlL7pw4YLmz5+v7t27u6d36dJFCxYs0IULF/SXv/ylDUfYfMOGDdPUqVMVFhZWa3pycrJ+9KMfSZLy8vLaYmhtYseOHTpx4oSmTp2qlJQU9/TY2FgtWrRIlZWVWrduXaP9rFmzRpK0dOnSWt+dmZGRoczMTO3atUv/+Mc/Wn4DfGjixIkaOXJknekjR47U6NGjVVxcrA8//LANRhbc2sOx1ZBXX31VX3/9tcaPH6/k5OS2Hk6L4tJeG3nppZdUVVWlGTNm1PkgbMyFCxe0du1anTt3TklJSUpPTw+67wF84403dPHiRUVERGjo0KFKT09v0r1kubm5kqQxY8bUmZeZmalf//rXysvL0yOPPNJSQ/YrNcdUaGjT/sQD+djy9DuvmeZNsMzNzVVMTIzuuOOOOvMyMzP19ttvKy8vT/3792/egP2UlWPn2LFjWrlypUpLS9WzZ0/deeedTb4UHwjKy8u1fv16FRQUKC4uTsOHD9ewYcO8Xr49H1v/9//+X0nSnDlzmrRcQBxbbX1KrD26cuWKue6664wkc/To0SYtK6nOy2azmR/84Afm0qVLrTRi36m5tPfdV//+/c3f//53r/u5+eabjSRTVFRUZ94333xjJJlbbrmlJYfuNyorK82QIUOMzWYzH374odfLBfqxdf/99xtJ5sCBA/XOT0xMNElJSR77qDk2Bg8eXO/8119/3Ugyv/rVr5o9Xn906tQpExERYa655hpTWVnZaPuayy/ffUVFRZnly5f7YMS+U3O5/LuvsWPH1vmKsvq052Pr5MmTJiQkxHTr1s2r48qYwDq2uLTXBvLy8nTixAmNGDGiyU+i/fKXv9R7772n4uJinT9/Xjt27NCtt96qjRs36oc//GErjdh30tLStHXrVp0+fVqXL1/W0aNHtWDBAp04cUJjxozRF1984VU/JSUlkqpvDv2umJgYdejQwd0m2CxatEgffvihZs+ercGDB3u9XKAfW55+55Jkt9sb/Z1708fV7YJJRUWFpk+frvLyci1fvtyrM8BJSUn6t3/7Nx09elSXLl1SQUGBNm7cqPj4eD3++ON64YUXfDBy33jooYeUm5urwsJCuVwu7du3T9/73vf09ttv65577pFppLZ1ez621q1bpytXrmj27NleX1kIqGOrrZNcoEhISKg3HTf02rlzZ4N9/eAHPzCSzLp161pkbJcuXTJ9+/Y1ksxHH33UIn02R0vuqxoLFy40kswvf/lLr8bQr18/I8lUVFTUO79Dhw6mf//+TdmsVtOS++uFF14wkkxKSoq5ePFis8fmb8eWJ3feeaeRZI4dO1bv/D59+pjw8HCPfRQUFBhJ5o477qh3/q5du4wk88gjjzR7vP6kqqrK/b708MMPN7u/Dz/80ISHh5vk5GRTVVXVAiP0T1VVVWbEiBFGknn99dc9tm3Px1aPHj2MzWYz//znP5vdnz8eW9wj5aUpU6bo4sWLXre/5ppr6p1+4cIFbd26VXa7vUmPF3sSHR2tKVOmaMmSJdqzZ48GDRrUIv1a1VL76mo//OEP9cwzz2jPnj1e9Vnzr76SkhIlJCTUmnfp0iVVVVU1+C9DX2up/bVu3TrNnTtXQ4YM0TvvvKOOHTs2e2z+dmx5cvXvvD4ul6vR37k3fVzdLhgYY/Twww9r48aN+sEPfqDVq1c3u8/Bgwfr1ltv1bvvvqvjx48H5T0/khQSEqLZs2dr9+7d2rNnj+66664G27bHY0uS3nnnHX3++efKyMhQ7969m92fPx5bBCkvNbUeT0M2bdqky5cva8aMGYqOjm6RPqX/qbFUWlraYn1a1VL76mpN3b5+/frpwIEDOnbsWJ0gVVP4zl9uom6J/fXiiy/q4Ycf1sCBA7V9+/Y629wc/nRseVLz+zx27JhuuummWvPOnz+voqIiDR8+3GMfMTEx6tKli/Lz81VVVVXnMoS/HTvNdeXKFc2ZM0fr1q3TlClT9NJLLykkpGXu+AiU46a5vN3O9nZs1bB6k7kn/nZscY+Uj7XGQSXJXY24V69eLdqvv2jq9qWlpUmScnJy6syrqWxe0ybQvfjii5ozZ44GDBigHTt2KCkpqUX7D5Rjy9PvvGaaN7/ztLQ0Xbp0qd6zn8F07FwdoiZNmqQNGzY06clYTyorK3Xw4EHZbDb16NGjRfr0V035+2gvx1aNr7/+Wv/5n/+p+Ph4TZgwoUX69Mtjq62vLbYnhw4dMpLM0KFDPba7dOmSOXr0qDl16lSt6QcPHqz36aktW7YYm81mEhMTW+S+mLZy9uxZc/z48TrTz5w5YwYMGGAkmc2bN9ea19C+Ki4uNg6HI6gLchpjzNq1a43NZjPXX3+9V0XugvnYqqioMH369DERERHm0KFD7ulXF+T87LPP3NMLCwvN0aNH6zxxdXXRxPLycvf0YCqaWFVVZWbNmmUkmQceeKDBewlrNLSv9u7da65cuVJrWkVFhVmwYIH7ibZg8PHHH5vz58/Xmf7uu++ayMhIExERUetvqj0fW1f793//dyPJPProow22CYZjiyDlQ/PmzTOSzMqVKz22q3ns87tVpmfOnGkcDoeZOHGiWbBggfnZz35mRo4caSSZyMhI88Ybb7Ti6Fvfzp07jc1mMyNHjjQPP/yweeKJJ8ykSZNMTEyMkWRmzpxZ5w+roX1ljDEbNmwwkkxiYqKZN2+eefTRR01ycrKRZDZs2OCjrWo927dvNzabzUgyP/rRj4zT6azzevXVV2stE+zH1o4dO0xYWJjp2LGjefjhh80vfvEL07t3byPJLF26tFZbp9PZYOXkOXPmGElm4MCB5le/+pWZMWOGiYiIMA6Hw3z88cc+2prWU7PtHTt2NL/5zW/qPXauDqMN7auePXuaXr16malTp5pf/epX5uGHHzb/8i//YiSZHj16mJMnT/p2w1qJ0+k0UVFRZvz48WbevHnmF7/4hcnMzDQ2m8106NDBrFmzpk779npsXW3w4MFGkjly5EiDbYLh2CJI+cjly5dNXFyciYiIaPRMSEMfdtnZ2eb73/++6dWrl4mOjjbh4eGmd+/e5oc//GGT61H5o88//9zMmTPHDB061MTFxZnQ0FCTkJBg7rzzzjpnomp4ClLGGPPWW2+ZUaNGmY4dO5qOHTuaUaNGmbfffrsVt8J31q1b1+gTfjNnzqy1THs4tt577z0zduxY43A4TFRUlLn55pvNxo0b67Tz9GFXVVVlVq5caQYNGmQiIiJMQkKCuf/++2ud0QpkDdVru/p19VPFDe2r3//+9yY9Pd107drVhIeHm+joaDN06FDzm9/8JijO+NbIzc01Dz74oOnbt6+JjY01YWFh5tprrzWTJ0827733Xp327fnYqvHee+8ZSSY1NdVju2A4tmzGNFL8AgAAAPXiZnMAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgEUEKAADAIoIU0A6dPHlSNptNs2bN8vt1+XKs3vr222/11FNP6brrrlN4eLhsNptyc3PbelhN4o/7FQhEBCkAaKJnn31WTz/9tHr06KHHH39cTqfT777UOTc3VzabTVlZWW09FK89/vjjstlsev/999t6KIDXQtt6AACCW7du3XT06FE5HI42Wb41vPnmm+rYsaNycnIUFhbW1sOxxB/368GDB9WhQwcNGTKkrYcCeI0gBaBVhYWFacCAAW22fGv44osvlJCQELAhSvLP/Xrw4EENGDBAUVFRbT0UwGtc2gOCWFVVlf71X/9Vffv2VWRkpPr27atly5bpypUrDS6za9cu3X333UpMTFRERIT69eunp556SqWlpfW2f/fddzVhwgQlJycrIiJC3bt318SJE7V7925Jnu/F2bp1q9LS0tS5c2dFRkaqe/fuGjt2rF577TV3G0/Lr1+/Xrfddps6duyojh076rbbbtP69evrtLv6MtfBgweVmZmp2NhYORwOTZgwQSdPnvS4H2tkZWXJZrMpPz9fp06dks1mk81mU3p6uiTppZdeks1m00svveRxDM0dV2P7PCsrS6NHj5YkLV682D1Om83m7tOf9uv8+fNls9l0/vx5ffzxx/WOF/BXnJECgtgjjzyiF198Ub1799ZPf/pTlZWV6bnnntPevXvrbb969Wr95Cc/UVxcnO6++24lJSVp//79evrpp7Vz507t3LlT4eHh7vZ//OMfNX/+fEVFRWnChAnq0aOHCgoKtHv3bv31r3/ViBEjGhzb888/r5/85Cfq0qWLJkyYoISEBJ09e1bvv/++XnvtNd17770et+2xxx7TihUr1K1bN/3whz+UzWbT1q1bNWvWLH3wwQd67rnn6ixz4MAB/du//ZvS09P1ox/9SIcOHdJrr72mDz/8UB999JEiIyM9rrMmMK1YsUKStGDBAklq9v1RTRmXN/s8PT1dJ0+e1Pr165WWluYetyR16tTJ41jaYr/efvvt+uqrr/SXv/xF48aN0y233CJJCgkJUc+ePb3fkUBbMACC0s6dO40kM2zYMPPNN9+4p585c8YkJiYaSWbmzJnu6R9//LEJDQ01KSkp5uuvv67V17Jly4wk8+yzz7qnHTlyxHTo0MF07drV5Ofn12p/5coVU1BQYIwxJj8/v866jDHmxhtvNOHh4ebcuXN1xl5UVOT+//qW37Vrl5Fkrr/+enPhwgX39AsXLpgBAwYYSebdd9+tsy8kmc2bN9da1/Tp040k8/LLL9cZR0N69uxpevbsWWf6unXrjCSzbt26OvNqxuB0Oi2Py9t93tD6ruZv+/WZZ54xksy2bdu8ag/4Cy7tAUHqP/7jPyRJv/3tbxUTE+Oe3q1bN/3sZz+r0/6FF15QZWWlVq5cqfj4+FrzHn/8cSUlJenll192T1u9erWqqqq0dOnSOmdkbDabunbt2ugYw8LC6r3PKCEhweNyNZfOsrKyat0s7XA45HQ6a7W52qhRozRp0qRa0x566CFJ0v79+xsdb2vxdlwtsc89acv9evjwYUnSDTfc0LRBA22MS3tAkPrggw8kSSNHjqwzr75p+/btkyS9/fbb+tvf/lZnflhYmD799FP3zzWPqI8ZM8bS+B588EH9+te/1uDBgzV58mSlp6drxIgRjV56kqRDhw5JUq1LVjVqptV8MF/txhtvrDPt2muvlSRduHDB26G3OG/H1dx93pi23K+HDh1S165d1blzZ6/aA/6CIAUEqZKSEoWEhCgxMbHOvOTk5DrTiouLJUlPP/20V/1fuHBBNptNXbp0sTS+xx9/XAkJCVq9erWee+45/e///b8VGhqqcePGacWKFerdu3eDy7pcLoWEhCgpKanOvOTkZIWEhKikpKTOvPoe9Q8NrX4brKqqsrQdLcHbcTV3nzemrfbrN998o+PHj2vcuHEWRg20LS7tAUHK4XDoypUrKioqqjPvq6++qjPNbrdLqv4wNcY0+KrRqVMnGWN09uxZS+Oz2WyaM2eODhw4oMLCQr366quaOHGi/t//+3+66667PH4A2+12XblyRYWFhXXmnTt3TleuXHFvjy+FhFS/pVZWVtaZV18Aaarm7vPGtNV+/eCDD2SMUUpKSov3DbQ2ghQQpIYNGyap+lH576pv2q233irpfy7xNSY1NVWSlJOTY3WIbgkJCbr33nv1yiuv6H/9r/+lo0eP6vjx4w22r/nAre9rWfLy8iS1zb02cXFxkqSCgoI682oumzVHU/Z5hw4dJDXtTFtb7dcjR45I+p9jFggkBCkgSM2YMUOS9Lvf/U6XLl1yTy8oKND/+T//p077n/zkJwoNDdX8+fN1+vTpOvMvXLhQKwzMnTtXHTp00FNPPaVTp07VauvNWZNt27bVOXNTUVHhvsToqSjjzJkzJVXXSHK5XO7pLpdLixcvrtXGl2688UbZbDZt3rxZZWVl7unHjh2rd583VVP2ec0DA2fOnPG6/7bar19//bUkqWPHji3eN9DauEcKCFLp6emaPXu21q1bpyFDhmjChAkqLy/XK6+8ottuu02vv/56rfaDBw/Wn/70J/34xz/Wv/zLv2jcuHG67rrr5HK59M9//lN5eXmaNWuWVq9eLUkaMmSIVqxYoUcffVSDBg3Svffeq549e+rLL7/Url27dNddd7nrLdVn0qRJio6O1ogRI9SzZ09VVFTonXfe0SeffKJJkyapR48eDS47atQozZ8/X3/4wx80ePBg3XfffTLGKDs7W6dPn9ajjz6qUaNGtch+bIpu3bpp0qRJ2rx5s2666SaNHTtW586d06uvvqqxY8dq69atzeq/Kft8wIAB6tq1qzZv3qzo6Ghde+21stls+vGPf9zg18K01X6tORM2f/583X///YqIiFBGRka9D0UAfqctai4A8I3KykqzbNky06dPHxMeHm769OljnnnmGXP8+PF6azsZY8z7779vJk+ebLp27WrCwsJMYmKiufHGG82vf/1rc/To0Trtd+7cacaPH2/i4+NNeHi4ufbaa819991n9uzZY4xpuI7Un/70J3PPPfeYnj17msjISJOQkGBuvfVW88ILL5iKigp3u4aWN8aYF1980dxyyy0mOjraREdHm1tuucW8+OKL9Y5RDdRU8tR/QxqqI2WMMZcuXTLz5883ycnJJiIiwgwdOtT8+c9/9lhHqqnjamyf19i3b59JS0szsbGx7npPNfWn/G2/Llu2zPTu3duEhoYaSWbLli1eLQe0NZsxV909CgAAAK9xjxQAAIBFBCkAAACLCFIAAAAWEaQAAAAsIkgBAABYRJACAACwiCAFAABgEUEKAADAIoIUAACARQQpAAAAiwhSAAAAFhGkAAAALCJIAQAAWPT/AQX5JclrO4VxAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# make histogram of decision function\n",
    "plt.figure()                                     # new window\n",
    "matplotlib.rcParams.update({'font.size':14})     # set all font sizes\n",
    "tTest = clf.predict_proba(X_test)[:,1]\n",
    "if hasattr(clf, \"decision_function\"):\n",
    "    tTest = clf.decision_function(X_test)        # if available use decision_function\n",
    "else:\n",
    "    tTest = clf.predict_proba(X_test)[:,1]       # for e.g. MLP need to use predict_proba\n",
    "tBkg = tTest[y_test==0]\n",
    "tSig = tTest[y_test==1]\n",
    "nBins = 50\n",
    "tMin = np.floor(np.min(tTest))\n",
    "tMax = np.ceil(np.max(tTest))\n",
    "bins = np.linspace(tMin, tMax, nBins+1)\n",
    "plt.xlabel('decision function $t$', labelpad=3)\n",
    "plt.ylabel('$f(t)$', labelpad=3)\n",
    "n, bins, patches = plt.hist(tSig, bins=bins, density=True, histtype='step', fill=False, color='dodgerblue')\n",
    "n, bins, patches = plt.hist(tBkg, bins=bins, density=True, histtype='step', fill=False, color='red', alpha=0.5)\n",
    "plt.savefig(\"decision_function_hist.pdf\", format='pdf')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "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.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}