{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import tensorflow as tf\n",
"import matplotlib.pyplot as plt\n",
"import os\n",
"import pickle\n",
"import math\n",
"\n",
"from sklearn.metrics import accuracy_score, roc_auc_score\n",
"trunc_normal= tf.truncated_normal_initializer(stddev=1)\n",
"normal = tf.random_normal_initializer(stddev=1)\n",
"\n",
"from xgboost import XGBClassifier\n",
"from architectures.data_processing import *\n",
"from architectures.utils.toolbox import *\n",
"from architectures.DNN import *"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# IMPORTING THE DATASET"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"l_index=0\n",
"mother_ID=[\"Ds\",\"Dplus\"]"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Bkg data amounts to 14066 while signal MC amounts to 1466 Ds and 1974 Dplus samples\n"
]
}
],
"source": [
"MC_Dplus_sig_dict, MC_Ds_sig_dict, data_bkg_dict = load_datasets(l_index)\n",
"\n",
"m_plus=MC_Dplus_sig_dict[\"Dplus_ConsD_M\"].shape[0]\n",
"m_s=MC_Ds_sig_dict[\"Ds_ConsD_M\"].shape[0]\n",
"n=data_bkg_dict[\"Ds_ConsD_M\"].shape[0]\n",
"\n",
"print('Bkg data amounts to {0} while signal MC amounts to {1} Ds and {2} Dplus samples'.format(n,m_s,m_plus))"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"#Normalising the Chi2 vertex fits to the NDoF\n",
"\n",
"MC_Ds_sig_dict[\"Ds_ENDVERTEX_CHI2\"]=MC_Ds_sig_dict[\"Ds_ENDVERTEX_CHI2\"]/MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"MC_Ds_sig_dict[\"Ds_OWNPV_CHI2\"]=MC_Ds_sig_dict[\"Ds_OWNPV_CHI2\"]/MC_Ds_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
"MC_Ds_sig_dict[\"Ds_IPCHI2_OWNPV\"]=MC_Ds_sig_dict[\"Ds_IPCHI2_OWNPV\"]/MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"\n",
"del MC_Ds_sig_dict[\"Ds_OWNPV_NDOF\"]\n",
"del MC_Ds_sig_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"#del MC_sig_dict[\"Ds_M\"]\n",
"\n",
"MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_CHI2\"]=MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_CHI2\"]/MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
"MC_Dplus_sig_dict[\"Dplus_OWNPV_CHI2\"]=MC_Dplus_sig_dict[\"Dplus_OWNPV_CHI2\"]/MC_Dplus_sig_dict[\"Dplus_OWNPV_NDOF\"]\n",
"MC_Dplus_sig_dict[\"Dplus_IPCHI2_OWNPV\"]=MC_Dplus_sig_dict[\"Dplus_IPCHI2_OWNPV\"]/MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
"\n",
"del MC_Dplus_sig_dict[\"Dplus_OWNPV_NDOF\"]\n",
"del MC_Dplus_sig_dict[\"Dplus_ENDVERTEX_NDOF\"]\n",
"\n",
"data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"]=data_bkg_dict[\"Ds_ENDVERTEX_CHI2\"]/data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"data_bkg_dict[\"Ds_OWNPV_CHI2\"]=data_bkg_dict[\"Ds_OWNPV_CHI2\"]/data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
"data_bkg_dict[\"Ds_IPCHI2_OWNPV\"]=data_bkg_dict[\"Ds_IPCHI2_OWNPV\"]/data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"\n",
"del data_bkg_dict[\"Ds_OWNPV_NDOF\"]\n",
"del data_bkg_dict[\"Ds_ENDVERTEX_NDOF\"]\n",
"\n",
"data_bkg_dict[\"phi_ENDVERTEX_CHI2\"]=data_bkg_dict[\"phi_ENDVERTEX_CHI2\"]/data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
"#data_bkg_dict[\"phi_OWNPV_CHI2\"]=data_bkg_dict[\"phi_OWNPV_CHI2\"]/data_bkg_dict[\"phi_OWNPV_NDOF\"]\n",
"data_bkg_dict[\"phi_IPCHI2_OWNPV\"]=data_bkg_dict[\"phi_IPCHI2_OWNPV\"]/data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
"\n",
"#del data_bkg_dict[\"phi_OWNPV_NDOF\"]\n",
"del data_bkg_dict[\"phi_ENDVERTEX_NDOF\"]\n",
"del MC_Ds_sig_dict[\"phi_M\"]\n",
"del MC_Dplus_sig_dict[\"phi_M\"]"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"#Number of input features\n",
"\n",
"dim=len(return_branches(mother_index=1, l_index=l_index))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"15"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"dim"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAtYAAAHoCAYAAABtiKiGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xuc1VW9//HXB7nI/ZJcFDEUvGWmKd7CYyCixsk75lSaph28hWUqapqQlqQ/EktNhVS0Y4x4qdDQxAI7dTQF9YSF4jUU5aaiJGBc1u+PvZlmcAb2DN/Zew/zej4e+8F813ft/f3MYmb7Zrn2+kZKCUmSJEmbp0WpC5AkSZK2BAZrSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQMGa0mSJCkDLUtdQENts802qW/fvqUuQ5IkSVu42bNnL00pdd9UvyYbrPv27cusWbNKXYYkSZK2cBHxj0L6uRREkiRJykBJgnVEbBURz0bEQ/njHSPiLxHxUkTcExGtS1GXJEmS1FClmrH+FjC32vE1wPiU0s7Ae8AZJalKkiRJaqCiB+uI2B74T+Dn+eMADgXuy3e5Ezi22HVJkiRJm6MUM9bXA6OAdfnjTwDLUkpr8sdvAr1re2JEjIiIWRExa8mSJY1fqSRJklSgou4KEhFfBBanlGZHxKD1zbV0TbU9P6U0AZgAMGDAgFr7SJKkpu2DDz5g8eLFrF69utSlaAvXqlUrevToQadOnTJ5vWJvtzcQODoihgFbA53IzWB3iYiW+Vnr7YG3ilyXJEkqAx988AGLFi2id+/etG3bltyKUSl7KSVWrlzJggULADIJ10VdCpJSujSltH1KqS9QAfwhpfRVYAYwPN/tVOA3xaxLkiSVh8WLF9O7d2/atWtnqFajigjatWtH7969Wbx4cSavWS77WF8MfCciXia35vq2EtcjSZJKYPXq1bRt27bUZagZadu2bWbLjkp258WU0kxgZv7rV4H9S1WLJEkqH85Uq5iy/HkrlxlrSZIkqUkzWEuSJEkZMFhLkiRlaMyYMUQEEUGLFi3o2rUr++23H5dddhkLFy7M5BoPPvggAwcOpEuXLnTq1Ik99tiDs846i3/+859VfSKCG2+8MZPrNcTMmTOJCJ5//vk6+7z++utVY/XnP//5Y+evuuoqIoK+fft+7Nwrr7zCGWecQZ8+fWjdujXdu3dn+PDhPPnkk1l+G/VSsjXWkiRJhRo/fV5Jrnv+0F0a9LzOnTvzyCOPAPD+++/zzDPPcPPNNzNhwgQeeeQR9t133wbXNHnyZL7yla9w5plncvnllxMR/PWvf+XOO+9k2bJldOjQAYAnnniCHXfcscHXKaYOHTowefJkBg4cWKP9nnvuqfp+qvvzn//MsGHD2Hnnnbnyyivp168fS5Ys4YEHHmDgwIG8++67dO7cuVjlVzFYS5IkZaxly5YceOCBVcdHHHEEZ599NocccggnnXQSL774IltttVWDXvvGG29k2LBh3HLLLVVtRx55JKNGjSKlf98/r/r1y91RRx3Ffffdx09+8pOqcZkzZw5z587lS1/6Ek888URV35UrV3LSSSex3377MW3aNFq3bl117oQTTuAb3/gGrVq1Kvr3AC4FkSRJKoouXbpw7bXX8sorrzB9+vSq9rFjx9K/f3+23nprevbsyZFHHrnRJSPLli2jV69etZ6rvsPFhktBUkp873vfq7rT4Omnn05lZSURweuvvw78e2nGlClTOPPMM+ncuTPbb789o0ePZt26dVWv9cILL1BRUUGfPn1o164de+yxB9dff32NPvVxzDHHsHz5cmbMmFHVVllZycEHH0zv3r1r9L333ntZsGAB48ePrxGq1xs8eDDt2rVrUB2by2AtSZJUJIMHD6Zly5ZV64Dvuusurr76ar7zne/wu9/9jptvvpn+/fvz4Ycf1vka++yzD5MnT+bGG2/krbcKv1n19ddfz9VXX81ZZ53FfffdR9u2bRk1alStfUeNGkWHDh247777OPnkk7nyyiu57777qs4vWLCAXXfdlZ/97GdMmzaN//qv/2L06NFcc801BddTXfv27fniF7/I5MmTq9oqKyv58pe//LG+jz/+ONtttx177rlng67VmFwKIkmSVCRt2rRhm222YdGiRQA89dRTHH744ZxzzjlVfY4//viNvsbVV1/NnDlzGDlyJCNHjmTHHXfk2GOPZdSoUXXOZK9du5Zrr72Ws846iyuvvBKAww8/nNdee4033njjY/0POeQQfvzjHwMwdOhQHnnkER544AG+9KUvATBkyBCGDBkC5GbCDz74YFasWMHEiRO59NJL6zkqORUVFZxxxhncfPPNPPfcc8yfP5/hw4fzox/9qEa/BQsWsMMOOzToGo3NGWtJkqQiqr4Oeu+992batGmMHj2ap556irVr127y+X369GH27Nk89thjXHDBBXTr1o3x48fzmc98hjfffLPW57zxxhssXLiQo48+ukb7hsfrHX744TWOP/WpT9V47VWrVjF69Gj69+9PmzZtaNWqFZdddhmvvfYaa9as2eT3UJthw4axdu1afve731FZWcmQIUPYZpttau1brjcRMlhLUmOYMbbUFUgqQ6tWreKdd96hZ8+eAJx++ulcffXVTJkyhQMOOICePXvyve99b5MBe6uttmLIkCGMGzeOWbNm8bvf/Y533323apZ5Q+vXbHfv3r1G+4bH63Xp0qXGcevWrVm1alXV8cUXX8y4ceMYMWIE06ZN4+mnn+byyy+v+h4bok2bNhx77LH88pe/ZMqUKVRUVNTar3fv3syfP79B12hsBmtJkqQimTFjBmvWrOGggw4CoEWLFpx//vnMnTuX+fPnc+GFF3L11VczceLEer3u4Ycfzl577cULL7xQ6/n1S0SWLFlSo33D40Lde++9jBw5klGjRnHYYYcxYMAAWrbc/BXGFRUVTJkyhaVLl3LcccfV2mfQoEEsWLCAv/3tb5t9vawZrCVJkopg2bJlXHzxxfTv35/DDjvsY+f79OnDJZdcQv/+/fn73/9e5+ssXrz4Y22rVq3izTffrJoJr+21e/XqxW9+85sa7VOnTq3nd5GzcuVK2rRpU3W8du1aKisrG/Ra1Q0dOpQTTjiBUaNG1bkP9fDhw+nduzfnn38+q1ev/tj5mTNnsmLFis2upSH88KIkSVLG1qxZU7Xzx/Lly5k9ezY333wzK1as4JFHHqnaq/nMM8+kW7duHHjggXTu3JkZM2bw0ksvbXR3jSOOOILddtuNo446ij59+rBw4UJuvPFG3nvvPc4888xan7PVVltx0UUXcdFFF9G9e3cGDhzI1KlTmTNnDpCbOa+PoUOHctNNN9G/f3+6devGTTfdxEcffVSv16hNy5YtmTJlykb7tG3blnvuuYcvfOELDBw4kHPPPZeddtqJpUuX8utf/5q7776bd955Z7NraQiDtSRJKnsNvQNiqbz//vscdNBBRASdOnWif//+nHzyyYwcObLGzh0HHXQQEydO5NZbb2XVqlX079+fiRMncuyxx9b52qNGjaKyspKLL76YxYsX0717d/bZZx/+9Kc/sf/++9f5vPPPP5/33nuPn/3sZ1x33XUcffTRfPe73+Wcc86hU6dO9fr+brjhBs466yzOPfdc2rZty6mnnspxxx3HiBEj6vU6DTVw4ECeeeYZrr76ai677DIWLVpEly5dOPjgg5k+fXpJ7roIENU/mdqUDBgwIM2aNavUZUhS7dZ/eHFww7adkpqruXPnsvvuu5e6jGbjG9/4BtOnT+cf//hHqUspqU393EXE7JTSgE29jjPWkiRJzcDzzz/PPffcw+c+9zlatGjBww8/zB133NHgm7ro4wzWkiRJzUD79u3505/+xI033siHH37IJz/5Sa655houuOCCUpe2xTBYS1IzM376vCa3XlXS5ttxxx2ZMWNGqcvYorndniRJkpQBg7UkSZKUAYO1JEmSlAGDtSRJkpQBg7UkSZKUAYO1JEmSlAGDtSRJkpQB97GWJEnlb8bY0lx38KX1fsqYMWP4/ve/D0BE0LlzZ/r378/hhx/OyJEj6dWrV9ZVsnr1am644QZuv/12Xn31Vdq1a0e/fv047rjjuOSSSwCYOXMmgwcPZs6cOXz605/OvIZCjBkzhhtvvJGlS5fW2WfSpEl8/etfp2PHjixatIi2bdvWOD9kyBD+8Ic/cOqppzJp0qQa52bOnMm4ceN48skn+eCDD9huu+04+uijufDCC9lhhx0a41uqwRlrSZKkjHXu3JknnniC//3f/6WyspLjjz+eX/ziF+y5557Mnj078+t985vf5IorruCrX/0qDz30EBMmTODzn/88Dz74YFWfffbZhyeeeIJ+/fplfv3GkFLit7/9bY22RYsW8fjjj9OhQ4eP9f/pT3/KoYceStu2bbn11lt57LHHGD16NM8++yzHHHNMUWp2xlqSJCljLVu25MADD6w6PuKIIzj77LM55JBDOOmkk3jxxRfZaqutan3uoEGDGDRoEGPGjCnoWitWrOCOO+7ghz/8IRdddFFV+/HHH09Kqeq4U6dONWoqd0cddRSVlZUMHz68qm3KlCn069ePjh071uj77LPP8p3vfIfLL7+cK6+8sqr9kEMO4etf/zoPPfRQUWp2xlqSJKkIunTpwrXXXssrr7zC9OnTM3vdDz/8kNWrV9e6xCQiqr6eOXMmEcHzzz9f1fbee+9RUVFB+/bt2W677bjmmmu48MIL6du3b1WfSZMmERHMmTOHoUOH0r59e3bbbTceeOCBGtf67W9/y9ChQ+nRo0dViH/00Ucb/H1VVFTw29/+luXLl1e1VVZWUlFR8bG+N9xwA9tssw3f+973an2tL37xiw2uoz4M1pIkSUUyePBgWrZsyZNPPpnZa3bv3p0+ffowZswYHnjggRpBdFNOO+00pk+fzk9+8hMmTJjAo48+yj333FNr36985SscffTR/OpXv2LnnXemoqKCN998s+r8a6+9xlFHHcUvfvEL7r//fj73uc/xhS98gT//+c8N+r4GDRpE165d+fWvfw3A/PnzeeKJJ2oN1o8//jhDhgyhVatWDbpWVlwKIkmSVCRt2rRhm222YdGiRVVta9eurbFkI6XEunXrWLNmTVVbixYtaNGi7vnQSZMmUVFRwQknnECLFi347Gc/S0VFBeeddx6tW7eu9TnPP/88U6dOZcqUKZx44olA7oOBffr0qXUN8/nnn8/pp58OwL777kvPnj156KGHOOuss4DcOu/11q1bx+DBg/nb3/7GbbfdxsCBAwsZnhpatGjBiSeeSGVlJaeccgqVlZV85jOfYffdd/9Y3wULFhTlw4mb4oy1JElSEVUP0QD9+vWjVatWVY8//vGPXHXVVTXaqq8brs2hhx7KK6+8wuTJkzn99NN55513uOiiizj00ENZt25drc+ZNWsWkFvLvF7btm057LDDau1/+OGHV339iU98gh49etSYsX7zzTc59dRT6d27Ny1btqRVq1Y8+uijzJs3b+MDshEVFRVMnz6dd999t85lIOtVX/ZSKs5YS5IkFcmqVat455136NmzZ1Xbgw8+yEcffVR1fOaZZ7LvvvsyYsSIqrbttttuk6/dsWNHKioqqKioIKXE6NGjueqqq3jwwQdr3RVj4cKFdOzYka233rpGe/fu3Wt9/S5dutQ4bt26NatWrQJyM9RHH300y5cv58orr6R///60b9+eK664gsWLF2+y9rocdNBBbLfddlx99dU8++yz3H///bX26927N/Pnz2/wdbJisJYkSSqSGTNmsGbNGg466KCqtj333LNGn44dO7LddtsxYMCABl8nIrjooou46qqreOGFF2oN1r169WL58uWsWrWqRrhesmRJva/38ssv8+yzz/Lwww9z5JFHVrWvXLmyYd9ANSeddBLjxo3jgAMOYMcdd6y1z6BBg5g2bRpr1qyhZcvSxVuXgkiSJBXBsmXLuPjii+nfv3+dyy0aYvXq1Sxbtuxj7S+99BJAjdnx6tYH96lTp1a1rVy5skE7lqwP0G3atKlq+8c//tHgDy5Wd+qpp3LUUUfxne98p84+I0eOZMmSJfzwhz+s9fy0adM2u45COGMtSZKUsTVr1lTt/LF8+XJmz57NzTffzIoVK3jkkUfq3MO6Id5//3122WUXTj31VAYPHkznzp158cUXGTt2LL179+a4446r9Xmf/vSnOeqoozj77LNZvnw5vXr14rrrrqNdu3Yb/aBkbXbbbTe23357LrjgAq666iqWL1/O6NGj6d2792Z/f5/61Keqdgapy9577811113Ht7/9bf7+979TUVHBNttsw2uvvcbtt9/O+++/z7Bhwza7lk0xWEuSpPLXgFuLl9L777/PQQcdRETQqVMn+vfvz8knn9wotzTv1KkTo0aNYtq0afzyl7/kgw8+oHfv3hxxxBFcfvnldO7cuc7nTpo0ibPPPpvzzjuPDh06cO6557LTTjvx9NNP16uGNm3a8MADD3DuuecyfPhwtt9+ey677DJmzpxZY9/sxnTeeeex5557Mm7cOL7xjW/UGIfqN85pTLHhJ1ObigEDBqT1n2aVpLIzY2zuzzIMA+Onz+P8obuUugypVnPnzq11OzUVx5o1a/j0pz/NAQccwJ133lnqcopmUz93ETE7pbTJRe/OWEuSJDVT9957L2+99RZ77rknH3zwARMnTuSll17irrvuKnVpTZLBWpIkqZlq3749d9xxBy+//DJr165lzz335MEHH2T//fcvdWlNksFakiSpmRo2bFhRPtTXXLjdniRJkpQBg7UkSZKUAYO1JEmSlAGDtSRJkpSBogbriNg6Ip6KiP+LiL9FxPfz7ZMi4rWIeC7/2LuYdUmSJEmbq9i7gnwEHJpS+mdEtAL+FBEP589dlFK6r8j1SJIkSZko6ox1yvln/rBV/tE0b/0oSZJUizFjxhARVY927dqx5557MmHChBr9Zs6cSUQU7ZbfAMOHD2fQoEH1es68efMYM2YMy5Yta5yitiBF38c6IrYCZgP9gZtSSn+JiLOBH0bEFcDvgUtSSh/V8twRwAiAHXbYoYhVS5KkUvrZcz8ryXXP2fucBj2vc+fOPPLIIwB8+OGHPPjgg5x55pl06NCBr3zlK1mW2OjmzZvH97//fU477TS6dOlS6nLKWtE/vJhSWptS2hvYHtg/Ij4NXArsBuwHdAMuruO5E1JKA1JKA7p37160miVJkuqjZcuWHHjggRx44IEMGTKE66+/nv32249f//rXpS5Njahku4KklJYBM4EjU0pv55eJfATcAXgfTUmStEXp2LEjq1ev3mifyspKWrduzS233FLVdu+997LzzjvTtm1bBg8ezLPPPktEMGnSpI2+1htvvMGwYcNo27Ytffv25ec///nH+rzwwgtUVFTQp08f2rVrxx577MH111/PunXrgNxylaOOOgqAHXfckYigb9++ALz99tucfvrp7LTTTrRt25ZddtmFyy+/nH/961/1GJUtS1GXgkREd2B1SmlZRLQFDgOuiYhtU0pvR0QAxwLFW2wkSZLUCNasWQPAihUrmDp1Ko8//ji33357nf0nTZrEiBEjmDBhAqeddhoAs2bNoqKiguHDh3PDDTcwd+5cTjrppE1eO6XEMcccw9KlS7ntttvYeuutGT16NO+++y4777xzVb8FCxaw66678tWvfpWOHTvy3HPPMXr0aFauXMmll17KPvvsw7hx47jwwgt54IEH2HbbbWnTpg0AS5cupVu3blx33XV07dq1ai32kiVLuPXWWzdj5JquYq+x3ha4M7/OugUwJaX0UET8IR+6A3gOOKvIdUmSJGXmnXfeoVWrVjXazjvvPL72ta/V2v+WW27hW9/6FnfddRcVFRVV7ddccw277747lZWVRARHHnkkq1ev5uKLa101W+Xhhx/m2Wef5cknn+SAAw4AYN9996Vfv341gvWQIUMYMmQIkAvjBx98MCtWrGDixIlceumldOrUiV133RWAz372s1Wz1QB77rkn48aNqzoeOHAg7du35/TTT+eGG26gdevWBYzUlqWowTql9Ffgs7W0H1rMOiRJkhpT586deeyxxwD46KOPmD17NldccQXdunVj9OjRNfr+9Kc/5a677qKyspLjjjuuxrmnn36aL3/5y+T+p37O0Ucfvclg/dRTT9GzZ8+qUA3wyU9+kn333bdGv1WrVjF27Fjuvvtu5s+fX2Opypo1a2jZsu6omFLiJz/5CRMmTOC1115j1apVVefmz59P//79N1rjlqjou4JIkiRt6Vq2bMmAAQOqjgcOHMjq1av57ne/y8iRI+nWrVvVufvvv5/+/ftz2GGHfex1Fi5cyIYbNhSygcPChQvp0aPHx9p79OjB8uXLq44vvvhifv7znzN69Gj22WcfunTpwm9+8xt+8IMfsGrVKjp06FDnNa6//nouvPBCLrnkEj7/+c/TtWtXnn76ac4999waIbs58ZbmkiRJRfCpT32Kf/3rX7zyyis12u+++24+/PBDjjrqKFauXFnjXK9evViyZEmNtg2Pa9OrVy8WL178sfYN2+69915GjhzJqFGjOOywwxgwYMBGZ6k3fO6JJ57ID3/4Qw4//HD2228/2rdvX9Bzt1QGa0mSpCJYfyOYPn361Gjffvvt+f3vf89LL73E8OHDayzH2G+//XjwwQdJ6d/305s6deomr7XffvuxaNEi/vKXv1S1zZ8/n2eeeaZGv5UrV1Z9GBFg7dq1VFZW1uizfq30hrPQGz4Xcv9IaM5cCiJJkpSxNWvW8OSTTwLwr3/9i9mzZ/ODH/yAY445hl69en2s/0477cRjjz3GIYccwsknn8zkyZNp0aIFF198MQcccAAVFRV8/etfZ+7cuUycOBGAFi3qnh8dNmwYe+21FyeeeCLXXHMNW2+9NVdcccXHlocMHTqUm266if79+9OtWzduuukmPvqo5j361n948dZbb6WioqLqTpJDhw7lpz/9KQcccAD9+vXj7rvv5uWXX96scWvyUkpN8rHvvvsmSSpbf7g69yhD1z36YqlLkOr097//vdQlbLbRo0cnoOrRqlWr1L9//zRq1Kj0wQcfVPWbMWNGAtKcOXOq2p555pnUuXPndPrpp6d169allFK65557Ur9+/VKbNm3SwIED0/Tp0xOQfvWrX220jn/84x/piCOOSFtvvXXaYYcd0i233JJOOOGE9PnPf76qz8KFC9Oxxx6bOnbsmHr06JEuuuiiNGHChASk5cuXV/UbN25c2mGHHdJWW22VPvnJT6aUUlq+fHk67bTTUteuXVPXrl3TGWeckR588MGPfU9NwaZ+7oBZqYB8Gqna/1poSgYMGJBmzZpV6jIkqXYzxub+HHxpaeuoxfjp8zh/6C6lLkOq1dy5c9l9991LXUZZ++///m9OOeUUXn31VXbcccdSl7NF2NTPXUTMTikNqLNDnktBJEmSytjZZ5/N0KFD6dq1K8888ww/+MEP+M///E9DdRkyWEuSJJWxd955h3POOYd33nmHT3ziE5x00klce+21pS5LtTBYS5IklbEpU6aUugQVyO32JEmSpAwYrCVJUllpqhsrqGnK8ufNYC1JkspGq1atPnb3QakxrVy5klatWmXyWgZrSZJUNnr06MGCBQtYsWKFM9dqVCklVqxYwYIFCz5245yG8sOLkiSpbHTq1AmAt956q8atvaXG0KpVK3r27Fn1c7e5DNaSJKmsdOrUKbOgIxWTS0EkqTGtvwOjJGmLZ7CWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMlDUYB0RW0fEUxHxfxHxt4j4fr59x4j4S0S8FBH3RETrYtYlSZIkba5iz1h/BByaUtoL2Bs4MiIOBK4BxqeUdgbeA84ocl2SJEnSZilqsE45/8wftso/EnAocF++/U7g2GLWJUmSJG2uoq+xjoitIuI5YDEwHXgFWJZSWpPv8ibQu47njoiIWRExa8mSJcUpWJIkSSpA0YN1SmltSmlvYHtgf2D32rrV8dwJKaUBKaUB3bt3b8wyJUmSpHop2a4gKaVlwEzgQKBLRLTMn9oeeKtUdUmSJEkNUexdQbpHRJf8122Bw4C5wAxgeL7bqcBvilmXJEmStLlabrpLprYF7oyIrciF+ikppYci4u9AZUT8AHgWuK3IdUmSJEmbpajBOqX0V+CztbS/Sm69tSRJktQkeedFSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQMGa0mSJCkDBmtJakbGT59X6hIkaYtlsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWpGbCPawlqXEZrCVJkqQMGKwlSZKkDBisJUmSpAwYrCUpazPGlroCSVIJGKwlqRkaP32eH2aUpIwZrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBQ1WEdEn4iYERFzI+JvEfGtfPuYiFgQEc/lH8OKWZckSZK0uVoW+XprgAtSSs9EREdgdkRMz58bn1IaV+R6JEmSpEwUNVinlN4G3s5/vTwi5gK9i1mDJEmS1BhKtsY6IvoCnwX+km/6ZkT8NSJuj4iudTxnRETMiohZS5YtNQLtAAAgAElEQVQsKVKlkiRJ0qaVJFhHRAfgfuDbKaUPgJuBfsDe5Ga0f1zb81JKE1JKA1JKA7p37160eiVJkqRNKXqwjohW5EL13SmlBwBSSotSSmtTSuuAicD+xa5LkiRJ2hzF3hUkgNuAuSml66q1b1ut23HA88WsS5IkSdpcxd4VZCBwCjAnIp7Lt30X+HJE7A0k4HXgzCLXJUmSJG2WYu8K8icgajk1rZh1SJIkSVnzzouSJElSBgzWkiRJUgYM1pLU2GaMLXUFkqQiMFhLkiRJGTBYS5IkSRkwWEuSJEkZMFhLkiRJGTBYS5IkSRkwWEuSJEkZMFhLkiRJGTBYS5IkSRkwWEuSJEkZaHCwjoiuEbF3RLTJsiBJkiSpKSooWEfE9yPiR9WODwXmA7OBVyJij0aqT5IkSWoSCp2x/irwQrXjHwN/AgYCLwJjM65LkiRJalIKDdbbAa8CREQfYC9gdErpSeA64MDGKU+SJElqGloW2G850Dn/9aHAeymlp/LHq4B2WRcmScrG+OnzSl2CJDULhQbrx4FLImIdcCHwm2rndgHeyLowSZIkqSkpdCnI+cBHQCWwDLis2rmvAX/MuC5JkiSpSSloxjqltIDcEpDaHAGszKwiSZIkqQkqdLu9P0TEbnWc7gX8LruSJEmSpKan0KUgg4BOdZzrBBySSTWSJElSE1WfOy+mDRsiojW5JSILM6tIkiRJaoLqXGMdEaOBK/KHCXgyIurq/v8yrkuSJElqUjb24cVpwFIggJ+Su9vi6xv0+RfwQkrpfxqlOkmSJKmJqDNYp5SeBp4GiIjlwG9TSkuLVZgkSZLUlBS63d6djV2IJG0RZowtdQWSpBIpKFhHRCvgW8DxwPbA1hv2SSn1yLY0SZIkqeko9Jbm44EzgYeAGeTWVkuSJEnKKzRYnwhcklL6cWMWI0mSJDVVhe5jHcBfG7MQSZIkqSkrNFhPBL7cmIVIkiRJTVmhS0EWAV+NiBnAdGDZBudTSunmTCuTJEmSmpBCg/X1+T93AD5fy/kEGKwlSZLUbBW6j3WhS0YkSZKkZsnALEmSJGWg4GAdET0i4pqI+H1EzIuIPfLt34qIgxqvREmSJKn8FRSsI2J/4CXgBOB1oB/QJn96W+CCxihOkiRJaioKnbEeT+6Oi7uQuwNjVDv3FLB/xnVJkiRJTUqhu4LsAxyTUloXEbHBuXeAHtmWJUmSJDUthc5Yvw90r+PcTuT2uZYkNTHjp88rdQmStMUoNFj/Bvh+ROxUrS1FxDbAhcADmVcmSZIkNSGFButLgA+AvwN/zLfdArwIrASuyL40SZIkqeko9AYx70XEgcApwBDgQ+Bd4OfAXSmljxqvREmSJKn8FfrhRVJK/wJuyz8aJCL6AHcBvYB1wISU0k8iohtwD9CX3HZ+X0opvdfQ60iSJEnFVug+1o9HxNkRUdcHGAu1BrggpbQ7cCBwbkR8itxSk9+nlHYGfp8/liRJkpqMQtdYLwHGAW9FxPSIOD0iutb3Yimlt1NKz+S/Xg7MBXoDxwB35rvdCRxb39eWJEmSSqmgYJ1SGk5ur+qvAf8EbgIWRsRDEXFKRHSs74Ujoi/wWeAvQM+U0tv5a72N+2JLampmjM09JEnNVqEz1qSUPkwpTU4pHUcu+P5X/tREYGF9LhoRHYD7gW+nlD6ox/NGRMSsiJi1ZMmS+lxSkiRJalQFB+vq8ss4XgFeI7cNX9tCnxsRrciF6rtTSuv3v14UEdvmz28LLK7juhNSSgNSSgO6d9/c5d6SJElSduoVrCNi/4j4cUTMJ7ef9eeBnwA7F/j8ILeryNyU0nXVTk0FTs1/fSq5G9JIkiRJTUZB2+1FxI+ALwGfBF4C7gAqU0pz63m9geT2wp4TEc/l274L/AiYEhFnAPOBE+v5upIkSVJJFbqP9ZeAKeTC9HOb6lyXlNKfgKjj9JCGvq4kSZJUaoXeeXGnxi5EkiRJasoKXmMdEW3yN4m5LSIejYid8+0nRcTujVeiJEmSVP4KXWO9CzAd6AzMBgYB6/eu/g/gP8ntcS1JkiQ1S4XOWP+U3IcK+wJHUHOd9OPAwdmWJUmSJDUthX548T+AE1NKyyJiqw3OLQK2zbYsSZIkqWkpdMZ6FXXfBKY3sCybciRJkqSmqdBgPR34bkR0rtaWIqINMBKYlnllkiRJUhNS6FKQi4A/Ay+TC9kJuALYA2gNHN8o1UmSJElNREEz1imlN4C9gFvIfYDxFXLrqu8F9k0pLWysAiVJkqSmoNAZa1JK7wHfyz8kSZIkVVPwDWIkSZIk1c1gLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUnFMGNs7iFJ2mIVHKwj4msR0aUxi5EkSZKaqvrMWN8B7AAQOVdERK/GKUuSJElqWurcxzoifgv8X/7xVyDI3XERcoF8NPAQ4M1hJEmS1Oxt7AYx04HPAl8EdiMXqm+MiBnA09QM2pIkSVKzVmewTildv/7riGgDrASeAXYFTiEXqn8REY8Aj6WUHmnkWiVJkqSyVeca64gYGREHR0THlNJH+eY7UkpfJheuA5gMdABubPxSJUmSpPK1saUgRwOXA9tExOvkZqgrIqItMCff5+GU0jONW6IkSZJU/uqcsU4pDU0p9QS2B75Jbob6MOBh4F1yQfvsiBiSXyoiSZIkNVub3G4vpfR2Sunh/OE3UkrdgAHkgnYfYBLwXqNVKEmSJDUBDb3z4tz8n99NKfUB9s2oHkmSJKlJ2tga6xpSStVDeAL+AXyUPze31idJkiRJzUTBwbq6lNI6YMeMa5EkSZKarIYuBZEkSZJUjcFakiRJyoDBWpIkScqAwVqSJEnKgMFakpq58dPnlboESdoiGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlSZKkDBQ1WEfE7RGxOCKer9Y2JiIWRMRz+cewYtYkSZIkZaHYM9aTgCNraR+fUto7/5hW5JokSZKkzVbUYJ1S+iPwbjGvKUnN2fjp80pdgiQ1G+WyxvqbEfHX/FKRrnV1iogRETErImYtWbKkmPVJkiRJG1UOwfpmoB+wN/A28OO6OqaUJqSUBqSUBnTv3r1Y9UmSJEmbVPJgnVJalFJam1JaB0wE9i91TZIkSVJ9lTxYR8S21Q6PA56vq68kSZJUrloW82IRMRkYBGwTEW8Co4FBEbE3kIDXgTOLWZMkSZKUhaIG65TSl2tpvq2YNUhSSc0YC4MvLXUVkqRGUPKlIJIkSdKWwGAtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUmSJGXAYC1Jm2vG2FJXIEkqAwZrSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQMGa0naQo2fPq/UJUhSs2KwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJYkSZIyYLCWJEmSMmCwliRJkjJgsJakhpoxttQVSJLKiMFakiRJyoDBWpIkScqAwVqSJEnKgMFakiRJyoDBWpIkScqAwVqSJEnKQFGDdUTcHhGLI+L5am3dImJ6RLyU/7NrMWuSJEmSslDsGetJwJEbtF0C/D6ltDPw+/yxJEmS1KQUNVinlP4IvLtB8zHAnfmv7wSOLWZNkiRJUhbKYY11z5TS2wD5P3uUuB5JkiSp3sohWBcsIkZExKyImLVkyZJSlyNJkiRVKYdgvSgitgXI/7m4ro4ppQkppQEppQHdu3cvWoGSJEnSppRDsJ4KnJr/+lTgNyWsRZIkSWqQYm+3Nxl4Atg1It6MiDOAHwFDI+IlYGj+WJIkSWpSWhbzYimlL9dxakgx65AkSZKyVg5LQSRJkqQmz2AtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSVug8dPnlboESWp2DNaSJElSBgzWkiRJUgYM1pIkSVIGDNaSJElSBgzWkiRJUgYM1pIkSVIGDNaSJLfnk6QMGKwlSZKkDBisJUlAbtbamWtJajiDtSRJkpQBg7UkSZKUAYO1JG2OGWNLXYEkqUwYrCVJkqQMGKwlqSGcqZYkbcBgLUmSJGXAYC1JkiRlwGAtSZIkZcBgLUmSJGXAYC1JkiRlwGAtSVsIb0cuSaVlsJYkSZIyYLCWJEmSMmCwlqRimzHWG8xI0hbIYC1JkiRloGWpC5AkNQ0/e+5nVV+fs/c5JaxEksqTM9aSJElSBpyxlqQtiFvuSVLpOGMtSZIkZcBgLUmSJGXApSCSpHrzg4yS9HHOWEvSFsC11ZJUegZrSZIkKQMGa0mSJCkDBmtJkiQpAwZrSZIkKQPuCiJJ9TFjLAy+tNRVlBV3CJGkHGesJUmSpAwYrCWplGaMbfBT3WJPksqLwVqSJEnKgMFakiRJykDZfHgxIl4HlgNrgTUppQGlrUiS6rAZyze2dH6QUVJzVjbBOm9wSmlpqYuQpKbCddaSVD5cCiJJkiRloJxmrBPwaEQk4NaU0oQNO0TECGAEwA477FDk8iRJ9eGyEEnNTTnNWA9MKe0DfAE4NyIO2bBDSmlCSmlASmlA9+7di1+hJGXJtdqStEUpm2CdUnor/+di4FfA/qWtSJIkSSpcWQTriGgfER3Xfw0cDjxf2qokqXnyA5GS1DDlssa6J/CriIBcTb9MKT1S2pIkSY3BtdeStlRlEaxTSq8Ce5W6DklSzvpZ6/OH7lLiSiSp6SiLYC1J2rJVn6WWpC2VwVrSFq+uUOcyBElSlgzWklQO1m+9N/jS0tZRZK63lrQlKYtdQSRJkqSmzmAtSZIkZcBgLUmSJGXANdaSmi3X90qSsmSwlrTFcEu3ps1/6Ehq6gzWkrQpM8aWxW4dzemmLYZsSU2RwVqSVNYM2ZKaCj+8KEmltn4Pa0lSk+aMtSQVwvArSdoEZ6wlSXUaP31e1dpuSdLGOWMtqUlzJxBJUrkwWEuSmgw/yCipnBmsJTU5zX2Wevz0ec1iy71N2fDnwKAtqdRcYy1JkiRlwGAtSZIkZcBgLUmSJGXANdaSmoTmvq5am+YHGyWVmjPWkiRJUgYM1pLUBHnTFkkqPwZrSZIkKQOusZZUtlxXnVPKfauf+eAeAH723CdKcv2Gcr21pFJwxlqSJEnKgDPWksqGM9SSpKbMYC1J2qK5LERSsRisJUnNhiFbUmNyjbUklZnqW+m5rZ4kNR3OWEuSmqW61vQ7ky2poQzWkrQxM8aW5LIbzlQ7cy1J5c9gLamkynYnkBljYfClpa5CktSEuMZakiRJyoAz1pKKrmxnqTdUomUgKi13DpHUUAZrSZLqYMiWVB8uBZGk6pylliQ1kDPWkoqiySz/kCSpgQzWkiQVwGUhkjbFYC1JUj15cxlJtXGNtSRtqITrrL0RjCQ1Xc5YS2o0rqtWc+NyEal5M1hL0mv/A++9X+oqtIVxuYjU/LgURJIkScqAM9aSMuXyD2nj6vs74gy31HQYrCVJKmOu25aajrIJ1hFxJPATYCvg5ymlH5W4JEkFcpZaKo66QrbhWyoPZRGsI2Ir4CZgKPAm8HRETE0p/b20lUmqzgAtlY+6fh8L/T2tbwA3vEubVhbBGtgfeDml9CpARFQCxwAGa6kEDNDSlm9zfs/d8USqXaSUSl0DETEcODKl9I388SnAASmlb27QbwQwIn+4K/BiUQttuG2ApaUuoglxvOrPMasfx6v+HLP6c8zqx/GqP8esfjZnvD6ZUuq+qU7lMmMdtbR9LPGnlCYAExq/nGxFxKyU0oBS19FUOF7155jVj+NVf45Z/Tlm9eN41Z9jVj/FGK9y2cf6TaBPtePtgbdKVIskSZJUb+USrJ8Gdo6IHSOiNVABTC1xTZIkSVLBymIpSEppTUR8E/gdue32bk8p/a3EZWWpyS1fKTHHq/4cs/pxvOrPMas/x6x+HK/6c8zqp9HHqyw+vChJkiQ1deWyFESSJElq0gzWkiRJUgYM1hmJiNsjYnFEPF/H+a9GxF/zj/+NiL2KXWM52dR4Veu3X0Ssze913qwVMmYRMSginouIv0XE48WsrxwV8HvZOSIejIj/y4/Z14tdYzmJiD4RMSMi5ubH41u19ImI+GlEvJx/P9unFLWWgwLHy/f+agoZs2p9m/37f6Hj5Xv/vxX4e9l47/0pJR8ZPIBDgH2A5+s4/zmga/7rLwB/KXXN5Txe+T5bAX8ApgHDS11zqR8F/Ix1IXe30h3yxz1KXXOpHwWM2XeBa/JfdwfeBVqXuu4Sjte2wD75rzsC84BPbdBnGPAwufsPHNic38sKHC/f++s5Zvlzvv8XOF6+9zdozBrtvd8Z64yklP5I7i+mrvP/m1J6L3/4JLm9uputTY1X3kjgfmBx41dU/goYs68AD6SU5uf7N/txK2DMEtAxIgLokO+7phi1laOU0tsppWfyXy8H5gK9N+h2DHBXynkS6BIR2xa51LJQyHj53l9TgT9j4Ps/UPB4+d5fTYFj1mjv/Qbr0jiD3IyP6hARvYHjgFtKXUsTsgvQNSJmRsTsiPhaqQtqAm4Edid3Q6o5wLdSSutKW1J5iIi+wGeBv2xwqjfwRrXjN6k9GDUrGxmv6nzvr6auMfP9v3Yb+Rnzvb8OGxmzRnvvL4t9rJuTiBhM7s314FLXUuauBy5OKa3N/YNSBWgJ7AsMAdoCT0TEkymleaUtq6wdATwHHAr0A6ZHxP+klD4obVmlFREdyM0WfruWsajtF7JZ79u6ifFa38f3/mo2MWa+/29gE+Ple38tNjFmjfbeb7Auooj4DPBz4AsppXdKXU+ZGwBU5t9UtwGGRcSalNKvS1tWWXsTWJpS+hD4MCL+COxFbn2Zavd14Ecpt9Du5Yh4DdgNeKq0ZZVORLQi9x+ju1NKD9TS5U2gT7Xj7cnN+jRLBYyX7/0bKGDMfP+vpsDfSd/7qylgzBrtvd+lIEUSETsADwCnNPd/RRYipbRjSqlvSqkvcB9wTnN9U62H3wD/EREtI6IdcAC5tWWq23xyszxERE9gV+DVklZUQvn1hrcBc1NK19XRbSrwtfzuIAcC76eU3i5akWWkkPHyvb+mQsbM9/9/K/B30vf+agocs0Z773fGOiMRMRkYBGwTEW8Co4FWACmlW4ArgE8AP8v/K3xNSmlAaaotvQLGSxvY1JillOZGxCPAX4F1wM9TShvdznBLV8DP2VXApIiYQ26Jw8UppaUlKrccDAROAeZExHP5tu8CO0DVmE0jtzPIy8AKcjM/zVUh4+V7f02FjJn+bZPj5Xv/xxTyM9Zo7/3e0lySJEnKgEtBJEmSpAwYrCVJkqQMGKwlSZKkDBisJUmSpAwYrCVJkqQMGKwlZSIixkREyj/WRcR7EfF0RPwwInqVur5iiohRETGoyNdsnf872DvD1/xmRGyRW0dFxAkR8XJEbLVB+14R8cuIWBAR/4qIdyPisYg4KSIK3qI2Im6MiHfyN6qo7fyFEbE2IraNiAH5vp039/uSVFoGa0lZeh84CPgcUEH+xhjk9hPdt5SFFdkocvtnF1Nrcvt0Zxast1QR0QL4PvD/Ukprq7UPB54GtgUuBQ4DTid344i7gFPrcZnJQDfg8DrOVwAzU0pvp5Rmkbu98vn1/FYklRmDtaQsrUkpPZl//C6lNBb4DPA2cM+Gs4ObEhFbN0qVZSIi2pa6hmZqCNAP+OX6hojoDUwC7gYOTSndlVL6Y0rp1ymlEeR+jl+pxzX+F/gHuQBdQ0T0B/YlF77XuwM4qz6z4pLKj8FaUqNKKS0jN4PbDxhaV7+IOC2/jGT/iJgZESuBi/Lnto6IayPijYj4KCL+LyKG1fIa/xURcyJiVUQsioj7qv/v9Yj4Uv78R/nX+mH1IFOthj0jYnpEfBgRL0TE8Rtc5+CI+J+I+CD/eC4iTsyfe53cnfZGV1saMyh/LkXEdyLi+ohYAsxZ/5yIGFfHeHSo1vaJiLg1It7Of48vRsS386eX5/+8o9p1+xY6fhHRJr98YVl++cN48nep3Jj88pOlEXFARMyKiJUR8aeI2DEiekTEryPinxExNyIO3eC5X8v3fTe/dGhGRAzYoM8eEfFIvs+H+dc5t5C/i404FXg0pbS8Wts3yN2N+IJUy53TUkovppRmblDbMfnveVVELMyP8fo7eybgHuCYWv6BWAGsBu6v1jaV3Az3EZuoXVIZM1hLKoYZwBrgwAL6TgYeInfb7IfybfcBpwFXA0eR+9/1U6PaeuKIuBy4FXgcOBY4m9zSlA7584eTCzrPAMcANwAXAjfWUsMvyQWd44CXgMqI2D7/Op3ydb0KnAAMB34BdMk/97j8dW8jtyzmoPw117uI3FKDU4DzChiP9d9fW2Bm/nu7Kj8+Pwa2y3dZH1p/UO26b+fbNjl+wI/IhcurgK8CnwQuKLC8dsAEYDzwZXK3Dv4Fub/LPwHHAwuAeyOiXbXn9SW3xOJE4CvAm8AfI2Knan2mAmuBk4Gjyf29dcyPyab+LupyKLkZ5eoOAWallN4t5BuOiC+RW+r0VL6u7wMjgLHVuk3O1/qfGzy9AngkpfTe+oaU0gfA38gtP5HUVKWUfPjw4WOzH8AYYOlGzr8N3LyR86cBCfjWBu1D8u2f36D9j8C9+a+7ACuA6zby+k8CMzZoG0UutG2/QQ2nV+vzCXL/KDgrfzwg36fjRq61FBhTS3sCnq2l/XVgXB3j0SF/fCawDti7jmt2yPc/rQHj9wlgJXBxtfMtgBfIT75u4u/9/7d3biFaVVEc/y2Z0ujqBbsaFmVoGln5UBRUWBhCRkka9NCFshs9VNDlIewiYUUQUUiGRlGWBZVhaaZ0QQxSShMvZQqVNZVOUlnOTLh6WPs459t+38xxGNPJ/w82Z87a17POzLDOOmuvUzM+cFuSPViSjUiyyxqM04fwGK8r+gGDUp9RDfp0eS/q9Dku9RmfydcCc+q0byqVPklmRJjH7KztDUmPA0uyNYWe0/npaf5r6sz1IrC06rWoqKjsf0UeayHEf4VVbDc/Ox8LNANLzaypKMBiwrCC8M4eQsSp7j5xxHafBbyRVb1OGHTnZvIPih/cfSvwC3BCEn0L/Am8mkIBuvKO5uTXV5WLCaP8yz3sV0V/o4B+wDtFJ3ffWT7vgjbg09L5hnRcUkd2fCEws+Fm9paZ/Uw84LQDpwHDUpMW4HtghkVWjsHZvN25F0WGmi2Z3AiDt0MQYSntpTI3VQ0jvPJzM50uIfQ4sjTMa8D4UkjPZOIhcF6dtW0prU8I0QuRYS2E2OukGNOBwM8VmudtBhHGRntWpgJDUpuB6fgT9RlExAvnYxfnAzL5tuy8jTCY8Hh9f2kaby7wq5nNz8IXOqOKDuoxkMbX1xlV9FcYc79kffPzRvyRDPGCtnTcpUd3L2T9AMzscOIBZghwF3ABMAZYSYeudxK6bgZmAc0pnnp0qu/OvSjinVsz+WY6Hp4K1qQ1jaE2nGdQOr5HrU43JfmQUts5xEPf5el8EjDP3bfXWVtraX1CiF6Idh8LIf4LLiL+3yyr0DbfONZCGD1XdNJnazoey+6eSJKsHcg9nkeX5qiMuy8DxqW457HAU0RcdpUY8np5oXcQ6fLK5Mb+VuCUPVlnoor+mtNxMLW6yPXVk5xLGLKXuPu6QmhZLudUd1XaFHgBMB2Yb2YnuPvObtyL4vpy7/YnwH1m1j8Z7Lj7X8DytK7yRsdijJuBL+rMURjYuPs3ZrYCmGxm64FTidj+ehzFHv4uCiH2L+SxFkLsVdLr+elEKMCH3RhiMeFR/dPdl+cltVlGxLbWzTPskat4BbFJrszVRNxyFYO/3rh/u/u7hDd1RKlql4e7Ij8AwzNZnkFlMTDazM5oMEaNRzjr15X+viKM+wlFJ4tczxPYexSpBnd5js3sPGJD4264e7u7LyEM52PJDONO7kXOJkJXJ2XyF4hwlCcqrH098bAytJ5OU/hQmTlEto9bCC/+ggbjDgW+rjC/EGI/RR5rIURP0mRmhafwcCJX761E1ohxXvoYxx6wCFgILDKz6UTmhCOID6H0c/f73X2bmT0CTDOzg4lX9H2JbAwPuftm4uMpC81sNhH3OorIgDHT3X+ouhgzG09sUnsb+I6IGZ5CbTzxOiKudgERA7zea1O75bwFPGNmDxAZO64kNrmVeQm4HfjAzKYSxt1JwDB3v8/d28xsE3C1ma0mDOVVFfW31cyeBx4ys39Sm5tIGVX2Ep8RuplpZo8T3uuphMEKQHqIeJKIhd8I9AfuBVa6e0vFe1GDu7cmD/LZlGLy3X2zmV0PvJJCSWYTm0oPI2LRzyDFRbv7TjO7G3g5ZSZ5nzDWTybeDExM3u6C1wmD/UZgViksJucc4iFUCNFb2de7J1VUVP4fhY7sEE54gbcRr9GnAcdU6H8dpSwYWV1fIp3ZBsKAaSa8fnlmhylEXGxrajMXOKJUP4nwzrYRXuJpQFNXa6CUtYPYXPcmsamuNY0zAxhQan82YThuT+NdmOQO3FHn+g4iPLHNwG/A00SYQc1aiDjrmUTs8w7CgL+zVH8pYUzvSH2HVtVfavMckSrwNyKt3V1UywqyJZNdmOYfmclrrh8YB6wm3jasIlIIfgS8meoHE+nzNqZraia8vydWvRcN1nwPsKFB3Zlpjh+J8KEWwlCfAhyUtb2M2LS5Hfid+Hrio+XfqVLbj9P1j20w72ji72bovv5bVlFR6X4x93rhfkIIIcT/EzM7mvBwn+/un+/r9QCY2WPAGHdXHmsheiySXkIAAAB6SURBVDEyrIUQQhxwmNmzwJHufu1+sJZDibzYEz37uqMQonehzYtCCCEORB4B1qYc5/uaE4GHZVQL0fuRx1oIIYQQQogeQB5rIYQQQgghegAZ1kIIIYQQQvQAMqyFEEIIIYToAWRYCyGEEEII0QPIsBZCCCGEEKIH+BcxG7ncCxibAAAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 864x576 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"#Convert data dictionaries to arrays for NN\n",
"\n",
"MC_Ds_sig = extract_array(MC_Ds_sig_dict, return_branches(mother_index=0, l_index=l_index), dim, m_s)\n",
"MC_Dplus_sig = extract_array(MC_Dplus_sig_dict, return_branches(mother_index=1, l_index=l_index), dim, m_plus)\n",
"data_bkg = extract_array(data_bkg_dict, return_branches(mother_index=0, l_index=l_index), dim, n)\n",
"\n",
"plt.hist(MC_Ds_sig[:,dim-1]/1000,alpha=0.5, label='Ds Signal MC',density=True, bins=30);\n",
"plt.hist(MC_Dplus_sig[:,dim-1]/1000,alpha=0.5, label='D+ Signal MC', density=True, bins=30);\n",
"plt.hist(data_bkg[:,dim-1]/1000,alpha=0.5, label='Bkg data',density=True, bins=200);\n",
"plt.ylabel('# events', fontsize=15)\n",
"plt.xlabel('D reconstructed mass (GeV)', fontsize=15)\n",
"plt.legend(fontsize=15)\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(12,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/mu/bdt_train.png', format='png', dpi=100)\n",
"#Cut on Ds Mass\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"#Cuts on mass\n",
"\n",
"#data_bkg_cut=data_bkg[np.where(data_bkg[:,dim-1]<1930)]\n"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"#Add 0/1 label for bkg/sig\n",
"\n",
"MC_Dplus_sig_labelled=add_labels(MC_Dplus_sig,signal=True)\n",
"MC_Ds_sig_labelled=add_labels(MC_Ds_sig,signal=True)\n",
"data_bkg_labelled=add_labels(data_bkg,signal=False)\n",
"\n",
"#Merge MC sig and data bkg, shuffle it\n",
"\n",
"data=np.concatenate((data_bkg_labelled,MC_Dplus_sig_labelled), axis =0)\n",
"data=np.concatenate((data,MC_Ds_sig_labelled), axis =0)\n",
"np.random.seed(1)\n",
"np.random.shuffle(data)\n",
"\n",
"#get train size\n",
"train_size=data.shape[0]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"#Strip away the label column and convert it to a one-hot encoding\n",
"\n",
"X=data[:,0:dim]\n",
"Y_labels=data[:,dim].astype(int)\n",
"Y_labels=Y_labels.reshape(train_size,1)\n",
"Y_labels_hot = to_one_hot(Y_labels)\n",
"Y_labels=Y_labels_hot\n"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"#Divide the dataset k \"equi populated\" sets\n",
"test=2\n",
"k=6 #number of subsets\n",
"i=0 #number of subset that is train set\n",
"\n",
"X_train, Y_train, X_test, Y_test, X_dict, Y_dict = k_subsets(i, k, X, Y_labels)"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"X_mean=X_train.mean(axis=0)\n",
"X_train_1=X_train-X_mean\n",
"X_std=X_train_1.std(axis=0)\n",
"X_train_2=X_train_1/X_std\n",
"\n",
"\n",
"\n",
"X_mean=X_test.mean(axis=0)\n",
"X_test_1=X_test-X_mean\n",
"X_std=X_test_1.std(axis=0)\n",
"X_test_2=X_test_1/X_std\n"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(14585, 14)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"X_train_2.shape"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# SETTING UP THE NETWORK"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"task='TRAIN'\n",
"\n",
"PATH=l_flv[l_index]+'/BDTs/test_'+str(test)\n",
"\n",
"if not os.path.exists(PATH):\n",
" os.mkdir(PATH)\n",
"\n",
"if not os.path.exists(PATH+'/variables_used.pickle'):\n",
" with open(PATH+'/variables_used.pickle', 'wb') as handle: \n",
" pickle.dump(return_branches(mother_index=0, l_index=l_index), handle, protocol=2)\n",
"\n",
"PATH=l_flv[l_index]+'/BDTs/test_'+str(test)+'/NN_'+str(i)\n",
"\n",
"if not os.path.exists(PATH):\n",
" os.mkdir(PATH)"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"\n",
"LEARNING_RATE = 0.001\n",
"BETA1 = 0.5\n",
"BATCH_SIZE = 64\n",
"EPOCHS = 50\n",
"VAL_PERIOD = 10\n",
"SEED=1\n",
"LAMBD=0.1\n",
"\n",
"sizes = {\n",
"'dense_layers': [\n",
" #(16, 'bn', 0.8, lrelu, tf.glorot_uniform_initializer()),\n",
" #(8, 'bn', 0.5, lrelu, tf.glorot_uniform_initializer()),\n",
" #(16, 'bn', 1, tf.nn.relu, tf.random_normal_initializer()),\n",
" (200, 'bn', 0.95, tf.nn.relu, tf.random_normal_initializer()),\n",
" (200, 'bn', 0.95, tf.nn.relu, tf.random_normal_initializer()),\n",
" #(200, 'bn', 0.90, tf.nn.relu, tf.glorot_uniform_initializer()),\n",
" (100, 'bn', 0.95, tf.nn.relu, tf.random_normal_initializer()),\n",
" ],\n",
"'n_classes':2,\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Input for propagation (?, 14)\n",
"Logits shape (?, 2)\n",
"Input for propagation (?, 14)\n",
"Logits shape (?, 2)\n"
]
}
],
"source": [
"tf.reset_default_graph()\n",
"model_NN = DNN(dim-1, sizes,\n",
" lr=LEARNING_RATE, beta1=BETA1, lambd=LAMBD,\n",
" batch_size=BATCH_SIZE, epochs=EPOCHS,\n",
" save_sample=VAL_PERIOD, path=PATH, seed=SEED)\n",
"\n",
"init_op = tf.global_variables_initializer()\n",
"saver = tf.train.Saver()\n",
"gpu_options = tf.GPUOptions(per_process_gpu_memory_fraction=0.13)\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" Training...\n",
"\n",
" ****** \n",
"\n",
"Training DNN for 50 epochs with a total of 14585 samples\n",
"distributed in 227 batches of size 64\n",
"\n",
"The learning rate set is 0.001\n",
"\n",
" ****** \n",
"\n",
"Evaluating performance on validation/train sets\n",
"At iteration 0, train cost: 2630, train accuracy 0.5969, test accuracy 0.7789\n",
"Average of signal as predicted by NN 0.2796\n"
]
},
{
"ename": "KeyboardInterrupt",
"evalue": "",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mKeyboardInterrupt\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-17-873a3d0d91d5>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 4\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'\\n Training...'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mmodel_NN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_session\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msess\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 6\u001b[0;31m \u001b[0mmodel_NN\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfit\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX_train_2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY_train\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mX_test_2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mY_test\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 7\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 8\u001b[0m \u001b[0msave_path\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msaver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msave\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msess\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mPATH\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m'/NN_model.ckpt'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/Rphipi/architectures/DNN.py\u001b[0m in \u001b[0;36mfit\u001b[0;34m(self, X_train, Y_train, X_test, Y_test)\u001b[0m\n\u001b[1;32m 288\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtrain_accuracy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 289\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 290\u001b[0;31m \u001b[0mfeed_dict\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfeed_dict\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 291\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 292\u001b[0m )\n",
"\u001b[0;32m~/miniconda3/envs/deepnet/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m 887\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 888\u001b[0m result = self._run(None, fetches, feed_dict, options_ptr,\n\u001b[0;32m--> 889\u001b[0;31m run_metadata_ptr)\n\u001b[0m\u001b[1;32m 890\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mrun_metadata\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 891\u001b[0m \u001b[0mproto_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtf_session\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTF_GetBuffer\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrun_metadata_ptr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/deepnet/lib/python3.6/site-packages/tensorflow/python/client/session.py\u001b[0m in \u001b[0;36m_run\u001b[0;34m(self, handle, fetches, feed_dict, options, run_metadata)\u001b[0m\n\u001b[1;32m 1099\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1100\u001b[0m \u001b[0mfeed_dict_tensor\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msubfeed_t\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp_val\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1101\u001b[0;31m \u001b[0mfeed_map\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcompat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_bytes\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msubfeed_t\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0msubfeed_t\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msubfeed_val\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 1102\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 1103\u001b[0m \u001b[0;31m# Create a fetch handler to take care of the structure of fetches.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m~/miniconda3/envs/deepnet/lib/python3.6/site-packages/tensorflow/python/util/compat.py\u001b[0m in \u001b[0;36mas_bytes\u001b[0;34m(bytes_or_text, encoding)\u001b[0m\n\u001b[1;32m 57\u001b[0m \u001b[0mTypeError\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mIf\u001b[0m\u001b[0;31m \u001b[0m\u001b[0;31m`\u001b[0m\u001b[0mbytes_or_text\u001b[0m\u001b[0;31m`\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0ma\u001b[0m \u001b[0mbinary\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0municode\u001b[0m \u001b[0mstring\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 58\u001b[0m \"\"\"\n\u001b[0;32m---> 59\u001b[0;31m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytes_or_text\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_six\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtext_type\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 60\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0mbytes_or_text\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mencode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mencoding\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 61\u001b[0m \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytes_or_text\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbytes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mKeyboardInterrupt\u001b[0m: "
]
}
],
"source": [
"with tf.Session(config=tf.ConfigProto(gpu_options=gpu_options)) as sess:\n",
" sess.run(init_op)\n",
" \n",
" print('\\n Training...')\n",
" model_NN.set_session(sess)\n",
" model_NN.fit(X_train_2, Y_train, X_test_2, Y_test)\n",
" \n",
" save_path = saver.save(sess, PATH+'/NN_model.ckpt')\n",
" hyper_dict={\n",
" 'k':k,\n",
" 'LEARNING_RATE':LEARNING_RATE,\n",
" 'BETA1':BETA1,\n",
" 'BATCH_SIZE':BATCH_SIZE,\n",
" 'EPOCHS':EPOCHS,\n",
" 'VAL_PERIOD':VAL_PERIOD,\n",
" 'SEED':SEED,\n",
" 'sizes':sizes,\n",
" 'LAMBD':LAMBD,\n",
" 'PATH':PATH,\n",
" }\n",
" with open(PATH+'/hyper_parameters.pkl', 'wb') as f: \n",
" pickle.dump(hyper_dict, f)\n",
" print(\"Model and hyperparameters saved in path: %s\" % save_path)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"task='TEST'\n",
"\n",
"with open(PATH+'/hyper_parameters.pkl', 'rb') as f: \n",
" hyper_dict = pickle.load(f)\n",
" #for key, item in hyper_dict.items():\n",
" # print(key+':'+str(item))\n",
"\n",
"k=hyper_dict[\"k\"]\n",
"LEARNING_RATE=hyper_dict[\"LEARNING_RATE\"]\n",
"BETA1=hyper_dict[\"BETA1\"]\n",
"BATCH_SIZE=hyper_dict[\"BATCH_SIZE\"]\n",
"EPOCHS=hyper_dict[\"EPOCHS\"]\n",
"VAL_PERIOD=hyper_dict[\"VAL_PERIOD\"]\n",
"SEED=hyper_dict[\"SEED\"]\n",
"sizes=hyper_dict[\"sizes\"]\n",
"LAMBD=hyper_dict[\"LAMBD\"]\n",
"PATH=hyper_dict[\"PATH\"]\n",
" \n",
"if not os.path.exists(PATH+'/hyper_parameters.pkl'):\n",
" print(\"No saved sizes dict\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"with tf.Session(config=tf.ConfigProto(gpu_options=gpu_options)) as sess:\n",
" \n",
" print('\\n Evaluate model on test set...')\n",
" saver.restore(sess,PATH+'/NN_model.ckpt')\n",
" print('Model restored.')\n",
" model_NN.set_session(sess)\n",
" model_NN.test(X_test_2,Y_test)\n",
" output_NN=model_NN.predict(X_test_2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" \n",
"plt.suptitle('S/B selection from NN', fontsize='15')\n",
"Ds_mass_MC =[MC_Ds_sig_dict[\"Ds_ConsD_M\"][i][0]/1000 for i in range(m_s)]\n",
"Dplus_mass_MC =[MC_Dplus_sig_dict[\"Dplus_ConsD_M\"][i][0]/1000 for i in range(m_plus)]\n",
"NN_selected = X_dict[i][np.argmax(output_NN,1).astype(np.bool)]\n",
"Ds_mass_sel_NN = NN_selected[:,dim-1]/1000\n",
"Ds_mass_train_NN =X_dict[i][:,dim-1:dim]/1000\n",
"\n",
"plt.subplot(1,2,1)\n",
"plt.hist(Ds_mass_MC+Dplus_mass_MC,bins=70, label='MC signal events');\n",
"plt.legend(fontsize='15')\n",
"plt.ylabel('# events', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"plt.subplot(1,2,2)\n",
"\n",
"plt.hist(Ds_mass_sel_NN,alpha=0.6,bins=70, label='S selected from test set');\n",
"plt.hist(Ds_mass_train_NN,alpha=0.2,bins=70, label='Test set (S+B)');\n",
"plt.legend(fontsize='15')\n",
"plt.ylabel('# events', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"fig=plt.gcf();\n",
"fig.set_size_inches(16,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+PATH+'/D_s_NN.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
" \n",
"true_positives_NN=output_NN[:,1][np.where(Y_test==1)[0]]\n",
"false_positives_NN=output_NN[:,1][np.where(Y_test==0)[0]]\n",
"\n",
"true_positives_NN=output_NN[:,1][np.where(Y_test[:,1]==1)]\n",
"false_positives_NN=output_NN[:,1][np.where(Y_test[:,0]==1)]\n",
"plt.hist(true_positives_NN,alpha=0.5,bins=80,density=True,label=\"True positives\");\n",
"plt.hist(false_positives_NN,alpha=0.5,bins=80,density=True, label=\"False positives\");\n",
"plt.xlabel(\"NN BDT output\", fontsize='15')\n",
"plt.ylabel(\"Events (a.u.)\", fontsize='15')\n",
"plt.legend()\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(16,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+PATH+'/fp_vs_tp_NN.png', format='png', dpi=100)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# SETTING UP XGBOOST"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"XGBClassifier(base_score=0.5, booster='gbtree', colsample_bylevel=1,\n",
" colsample_bytree=1, gamma=0, learning_rate=0.1, max_delta_step=0,\n",
" max_depth=3, min_child_weight=1, missing=None, n_estimators=100,\n",
" n_jobs=1, nthread=None, objective='binary:logistic', random_state=0,\n",
" reg_alpha=0, reg_lambda=1, scale_pos_weight=1, seed=None,\n",
" silent=True, subsample=1)"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"PATH=l_flv[l_index]+'/BDTs/test_'+str(test)+'/XG_'+str(i)\n",
"if not os.path.exists(PATH):\n",
" os.mkdir(PATH)\n",
"# fit model to training data\n",
"model = XGBClassifier()\n",
"\n",
"model.fit(X_train_2, Y_train[:,1])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [],
"source": [
"output_XG=model.predict_proba(X_test_2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7sAAAIgCAYAAABETH7CAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3XmczXX///HnC8Mw2YZZmDDKGkmZLGmRLVkSEulXRFlSfVHX0EWhK1rkqqsrFBWSy5KyFCFd0qbQ1UZos+9LURnbeP/++Jw5nZk5M4bGjDnzuN9un5tz3p/35/N5fY5zPmde5/3+vN/mnBMAAAAAAKGkQG4HAAAAAABAdiPZBQAAAACEHJJdAAAAAEDIIdkFAAAAAIQckl0AAAAAQMgh2QUAAAAAhBySXQDIh8ysh5l9YWa/mdkvZvalmf0zg7rfmVk/32MXsJwys51mNsvMKudQ3JvN7Jls3me0mY0ws/g05U1851k7O493mljam9l6MztuZptz6rhZYWbX+V6PzmnKC/rePx+bmQWUFzGzgWa22vc+O2pmP5jZs4HvFzOLT/O+SjazrWY2ycyicvIcA2JqaWYDcuPYAIDsQ7ILAPmMmT0s6WVJSyR1lHSnpPmSbgpSt7KkmpIWBRSPldRIUmNJD0m6QtJCMyt0biM/Z6IlDZcUn6b8f/LO86ecCMLMCkp6TdLXkppK6pATx80q59wKSa9LetbMigesul9SbUn9nHNOksysmKT3JY2Q9z67RdKNkp6X1EzSh0EO8ZC81/taSY/Jez9OPxfnkgUtJZHsAkAel1f/MAEAnL37JL3knPt7QNnbZjYySN02ktY657YElG12zn3me7zSzH6VtFBSNUnfnZOIc4Fz7rCkz05bMfuUk1RC0n+ccx9nVMnMijrnknIurFQekrRB0khJg8ysvLzE9Hnn3LcB9UZJqiupgXNuXUD5cjMbJ6lXkH1vDHhffWJmhSW9YGYXOOd+z/YzAQCEPFp2ASD/KSVpd9rClFa5NNrIS2Qz85vv37DMKplZLzNbZ2ZJZrbfzFaYWa2A9eFm9rSZbTOzY2b2tZm1Ps2xZWZX+/Z1xMwO+Lq/Fk9Tp5KZzfAd94iZfWNm3Xxdl1OStOUpXWl926TrxmxmxczseTPb7euWu9rMWqY51gdmNse3/x/N7LCZvWtmF2ZyDj0kbfM9ne877gjfOmdmg8zsOTPbFxCvzOw+X9fgY75jDUyz3xG+c25gZmt8r/3HZlbZ1317npn97us63fR0r7Vzbo+kYZLuN7M6kp6TdFhey7j/NZLUW9L4NIluyj5OOecmne5Y8t5XJqlgwL4L+s5pq++c15lZt7QbmtmtZvatr842MxsV2PPAzEqZ2cvmdcM/mtJtOuU1k/SgpEoBXaunZCFeAMB5hpZdAMh//icvWdkq6R3n3IFglXxJSxNJo9OsKuBLHExSZXmtfD9IWpvRAc3sWkkvSnpU0kp5LZiNJJUMqDZHUn15idNPkm6VtMDMEpxzX2Ww38byusvOk9dVtoykJyWV9j2XmUX7jnlEXsvkNnndbitI2iXpdnndZfv7XpvMTJLXvfbvkn6UdI+8LtzXp2mNbSCpvLykqaikf0maKCmj5H2hvC7lb/li/ETS9oD1f5PX9fcO+X6oNrN7JP1b0j/ldRW+XtJYMyvinHsyYNtivmM/LekPeV2Jp0k6JuldSeMlJUp6w8wqOOeOnOY1mCDpLnld3+Ml3Zqm5bWe75hLT7OftFLeVwUl1fCd83Ln3KGAOo/5Yh0pabWkTpKmm5lzzs2QvPttJc2S1yX8b5LqSPqHvPdGX99+/inpKkkD5f3wU0Fe92nJ6+JfVam7ku87w3MBAJwPnHMsLCwsLPlokffH/8+SnKRTktbJSyJKpKnXTtJBSQUDylyQZZukS09zzIckfZHJ+ma+fV2XpvxDSW8EPN8s6ZmA5x/JS4gCt2nq21dt3/Mn5CV55TI4dm1f/SZpypuk2U9N3+vVPaBOAXlJ/pKAsg8kHZJUOqBsgG9fRTN5DeJ9ddqmKXeSvkxTVkDSDkmT05SP9x073Pd8RNrXVdK9vrJHA8ou8ZXdmMX3UFtf/RVB1nXxraseJOZCKUuQ8067rJMUF1Av0vf/ODzNfhfJ6wKd8vyzIO+JREnJki70PV8r6f5Mzu8Zed31c/3zysLCwsJy9gvdmAEgn3HOfSMvcbtJXnJkkh6RtMbMLgio2kZeEpecZhdjJF3pW9pI+kbSIjOLy+SwX0m63LyReK/13Y8ZqLm8FrZPzKxQyiKv1TYh2A59Lc+NJM1Os83Hkk7Ia2GUvOR3sXNuVybxZcWV8l6rN1IKnHOnfM+vTlN3tXPul4DnKfcyZ/YaZSZtV/IL5bUcv5GmfJa8VvNLA8qOy/tRIMWPvn//G6Qsq/H1lu+HADMrk2ZdyojMabvFL5D3/3JC0glLP8r1QHmvcX15LaqHJb0b8J6sLa/FONg5V/N1yy4ob8C0YHUKyHu/SN778W9mdq+ZVTvdyQIA8iaSXQDIh5xzx5xzbzvn7nPOXSLpbnldNwMHDmqt4PfrbnXOrfEti+R1vw2Xl6xkdLxl8rq+Xiuv5XO/mY03swhflbKSYhWQDPmWEfK6mAZTWl6X1/Fptjkm7/7hlO3KyOuu/FeVk/S7S9/Nd4+kYmZWJKDs1zR1jvv+DT/LY+8JEkuw8pTnkQFlv/mS8rSx+GN0zmU5PjNrJ6/Vv5u8lu4n01TZ4fs37T3KA+Qls30V3I++99Rq59w8eT/G1JLUw7f+dOdcWt77KCyTOimvy33yur4/Kmmj777nrhnEBQDIo7hnFwAg59wrZva0vHsl5Rt8KE7S4ixse8zMfpbXWpxZvamSppo3d2pHSc/Ka70bIq+79A5JN59B2L/Kaz0codRTI6XY6fv3gP5MlP6KXZIuMLNiaRLeGElHnHPHsuEYGUnbSpqSvEenKY/x/XvwXARhZkXl3fP7hnNupq/VdaKZTXLOrfJV+0Le/dEtFdB67Jz70bePC5QFzrl9ZrZff76vAs858D7zwHM+KO8Hj0xfF+fcr5IekPSA772eKO/e32+ccyEzojgA5He07AJAPuMbsCltWZS8waJSWsDaSPrcObc/C/sLl3Sx/hxNOFPOuX3OuZfkda29xFf8vryW3d8DWo39Swb7+UPe/ZnVg23jnNsZsO8bzCwm2H6U9VbX1fKSzltSCszMfM8znCroHNkuL5nvnKb8Vnk/IHybbovs8Yi81tOUVvxXJK2SNN7MCkiS74eAiZL6m1mmP4Bkxvf/VVZ/vq/Wykuig53z9773VbK8ZDtYnVPyBipLxdet/2/y/iaq4Ss+rrNvhQcAnCdo2QWA/OdbM5svb7TcvZIqyRtA6oikqb46mU05FG9mDX2Po+QNeFRSXuITlHlz+EbK14VZ0uWSrpPXqitJ78kbUfg9M3tK3uBEJeTN1RrunHs4g10nSnrfzE7JG835N0kVffEPdc59L68F+U5JH5nZKHnJU01JEc65pyVtlZQkqbuZHZJ0IliC7Zxbb2Yz5M39WkJ/jsZcQ1K/jM79XHDOnfJNkfOSmR2Q9/pd54vj7865o9l9TDOrIW906aHOuR2+OJyZ3Svvh4A+8kZqlqSh8u69XWlmL8j7YeOovN4C3eUNFpU2xuq+llzz1fubpN8lzfAd66CZPSdpmJmdlLRGXg+B1pJuC9jPcElLzGyypJny7l/+h6RJzrntvnP5WNJceQm0k/f/+Ie8xF3y5hKOMW9KqLWS9jvnNp/VCwcAyDUkuwCQ/zwmqb287qiR8gaG+lRSF+fcJjOLlNRQ3n2NwTzoWySvO+m3klo651ZncszV8loDu0oqLmmLvO7H/5L8SVNHeVP6DJCXsB6UN5DQvzPaqXPuY9+0RiPlTadT0LfvxfK1Uvu6wzaWN/XOc5KKyJsq6Qnf+qO+aXyGS1oh755PU3D3SHpKXgtnKd+5t3Wppx3KEc65Sb77hAdI+j95rb0POueePUeHHCfvdXsuTRz/M7MXJY0yszecc/udc0d88/b2l3dv7wB5f3Nsk9fSfllKt+YAzwQ83iMvme3jnNsSUP6opJPykvoYeT84/D/n3MyAeJb67r8dJm9aqb2SxipgLmB5Lbw95I0EnSzpS3kjUadM9zRb3lROT8v7QWeq/rx3GACQR5hzaW8DAgDkZ2bWTdLTzrm0AwwBAADkGSS7AAAAAICQwwBVAAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5hXI7gOxWtmxZFx8fn9thAABCxBdffLHfOReV23HkZXw3AwCyU1a/m0Mu2Y2Pj9eaNWtyOwwAQIgwsy25HUNex3czACA7ZfW7mW7MAAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkBNyozEDOLcOHz6svXv36sSJE7kdCvCXhYWFKTo6WiVKlMjtUPI9ri3A+YdrJPI6kl0AWXb48GHt2bNHcXFxKlq0qMwst0MCzppzTklJSdqxY4ck8cdcLuLaApx/uEYiFNCNGUCW7d27V3FxcSpWrBh/jCLPMzMVK1ZMcXFx2rt3b26Hk69xbQHOP1wjEQpIdgFk2YkTJ1S0aNHcDgPIVkWLFqXrbC7j2gKcv7hGIi8j2QVwRmh1QajhPX1+4P8BOD/x2UReRrILAAAAAAg5JLsA8pURI0bIzFS1atWg66tUqSIz04gRI9Kte+utt9S0aVOVKlVKRYoUUbVq1TRs2DDt37//L8VTtmzZs94+OzRp0kS33HJLrsZwNiZOnKh58+bldhiA35QpU1SvXj0VL15cpUuX1uWXX65Bgwad8+N+8MEHMjOtXbs2W/e7dOlSPffcc9m6z7Jlywa9vgbatWuXWrdurZIlS8rM9MEHH2RrDH/Vubz25OR17fjx4xoxYoS++uqrHDkekBtIdgHkO+Hh4dq0aZPWrFmTqnz16tXasmWLwsPD023z4IMPqnPnzrrooos0bdo0LV26VAMHDtTbb7+te+6556xjufvuu7VkyZKz3j4/I9nF+eSJJ57Q3XffrRtuuEFvvfWWXnvtNbVv314LFizI7dDO2rlIdrNi1KhR+vrrrzVjxgytXLlSV1xxRY7HkJlQSnZHjhxJsouQxtRDAP6y+CELc+W4m59sc1bbRURE6IorrtDMmTOVkJDgL585c6aaNm2qL774IlX9t99+W//85z/1yiuvqGfPnv7y6667Tr1799bSpUvP7gQkXXjhhbrwwgvPensglD381re5ctwnOl56xtu88MIL6tOnj0aPHu0va9eunYYPH56doeULGzZsUIMGDdS6desM6yQnJys5OVmFCxfOwcgA5DW07ALIl7p27arZs2fLOSfJm09w9uzZ6tq1a7q6zz77rK644opUiW6KggUL6sYbb8zwOL/++qvuvvtulS9fXuHh4apYsWKqluBg3Zi/+eYbXXXVVQoPD1etWrW0aNEiJSQkqEePHv46PXr0UEJCgt577z3VqVNHERERuvrqq7Vu3bpU+xo7dqyuvPJKlSxZUjExMWrXrp1+/PHHLL1Gac2fP18JCQkKDw9XbGysEhMT/SN0Ll++XGaW7vi//PKLChcurFdeecVf9vHHH+u6665TsWLFVKZMGd1zzz367bff/OunTJkiM9O3336rFi1aKCIiQjVq1NBbb73lr9OkSRN98cUXmjp1qsxMZqYpU6ZIkhYsWKB69eopIiJCpUuXVoMGDbRixYqzOmcgq3799VfFxsamKz/d4D6nu0ZI0tq1a9WmTRsVL15cxYsXV+fOnbV79+5M93vq1Ck9+eSTqlKliv+2i6lTp6arN3fuXNWvX19FixZVmTJl1Lp1a23ZskUjRozQ2LFjtWXLFv9nLPAadLrPsSR9+OGHuuyyyxQeHq569erp008/zTRmyXu93n//fc2dO1dmpvj4eEl/XvPmzZunWrVqKTw8XJ9//rkk6auvvlKzZs1UrFgxlS5dWrfffrv27Nnj3+fmzZtlZpo5c6buuusulShRQhdeeKFef/11SdLTTz+t8uXLKyoqSoMHD9apU6cyjC+za48kvfzyy6pVq5aKFCmiSpUq6emnn061/bp169SqVStFRkYqIiJCNWvW1Lhx47K077SeeOIJValSReHh4YqJiVGrVq1SvS8OHjyoPn36KCYmRuHh4brqqqv8r5kkFS9eXJJ01113+Y+3efPmjP9zgDyIZBdAvtSxY0ft2bNHH3/8sSTpo48+0r59+9ShQ4dU9U6cOKFPP/1UrVq1OqvjDBo0SB9//LGeffZZLVmyRKNHj870j98jR47ohhtuUFJSkmbMmKFhw4Zp4MCB2rp1a7q6W7du1d/+9jcNHTpUM2bM0N69e3Xrrbf6E3hJ2r59u+677z7Nnz9fkyZNUnJysho3bqxDhw6d0XnMnj1bHTt2VP369bVgwQINHz5cEydO1MMPPyzJa+UuV66cZs+enWq7uXPnSpL/df3kk0/UrFkzxcbGas6cOXruuee0aNEi3XXXXemO2a1bN910002aO3euqlatqq5du2r79u2SpPHjx6tGjRpq3bq1Vq5cqZUrV6pNmzb66aefdMstt6hp06Z6++23NX36dLVt21YHDx48o/MFztQVV1yhf//735o6daoOHDiQ5e1Od4348ccf1bhxYx09elTTpk3TlClTtG7dOrVr1y7VZz2t+++/X48//rh69+6thQsXqkOHDurZs6feeecdf51p06apY8eOuvjiizV79mxNnjxZ1apV0759+3T33XerW7duio2N9X/GHnnkEUlZ+xzv3LlTN954oyIjIzVnzhz16dNHt99+u44cOZLp67Fy5Updfvnluv7667Vy5Ur/NUTyktbExEQ9/PDDWrRokSpXrqx9+/apSZMmOnLkiP7zn//o3//+t1asWKEWLVro+PHjqfY9ePBglStXTm+++aauueYade/eXQ8++KBWrVqlV199VQMGDNDTTz+d7joWKKNrjySNGTNG/fr1080336x33nlH/fr10yOPPKIXXnjBv/1NN92kggUL6vXXX9eCBQt0//33+38kyGzfab322msaPXq0Bg0apCVLlmjChAmqUqWK/vjjD0nSsWPH1Lx5c7333nsaM2aM5s2bp6ioKDVv3tyfEP/3v/+VJA0bNsx/vHLlymX6/wPkNXRjBpAvlSpVSq1atdLMmTN1zTXXaObMmWrVqpVKlSqVqt6BAwd07NgxVaxY8ayOs2rVKvXv319dunTxl/2///f/Mqw/efJkHThwQGvWrFFcXJwk6eKLL1aDBg3S1T148KA++eQT/2Bbp06dUocOHbRx40bVqFFDktcqnSI5OVktWrRQdHS05s+frzvvvDNL5+Cc09/+9jfdeeedGj9+vL+8SJEi6t+/vx5++GGVKVNGnTt31qxZszRy5Eh/nVmzZqlly5aKjIyUJA0ZMkRXXXWVZs2a5a8TFxenZs2aae3atapdu7a/fODAgf7W9Hr16ikmJkbvvPOO+vbtq0suuUQRERGKiopSw4YN/dusWLFCxYsX15gxY/xlmXWFBLLLuHHjdPPNN6tHjx4yM9WsWVOdOnXSQw89pBIlSmS43emuESNHjlRsbKzeffddf5fdOnXqqEaNGlq0aFHQZOjHH3/UhAkTNHnyZHXv3l2S1Lx5c+3atUsjR45U27ZtderUKQ0ZMkQdOnTQjBkz/NvedNNN/sflypVTkSJFUn3GpKx9jp977jmFh4dr4cKFKlasmCTvFpLMrn+S1LBhQ5UoUUKRkZHpjnvgwAEtW7ZMdevWTRWLJC1ZssT/OlerVk0NGjTQm2++qdtuu81ft2nTpv5u5g0aNNCcOXO0YMECbdiwQQULFlSrVq00f/58zZ07N2gvH0kZXnsOHz6skSNHatiwYf6u6y1atNCRI0f0+OOPq1+/fvrll1/0888/a968ebr0Uq+rfLNmzU6772BWrVqlli1b6t577/WXdezY0f/49ddf19q1a7Vu3Tr/d0Tz5s1VvXp1jR07VmPGjNGVV14pyfuOOd3xgLyKll0A+VbXrl01Z84cHTt2THPmzMnwjxvp7OcZrFu3rsaMGaPx48fr+++/P2391atXq169ev5EV5Lq16+vmJiYdHXj4+NTjSp9ySWXSJK/9VOSPvvsM7Vo0UJlypRRoUKFVKxYMf3+++9ZiiXF999/r61bt+rWW2/VyZMn/UvTpk119OhR/wiwXbp00caNG/X1119Lkvbv36///ve//j/ijxw5opUrV6bbz9VXX62wsLB090q3bNnS/7hMmTKKjo5OdW7BXHrppTp06JC6d++upUuX+ls5gHOtTp06Wr9+vRYsWKB7771Xzjn94x//UEJCgn7//fcMtzvdNWLZsmXq0KGDChQo4P/MVK5cWfHx8ekG2Uvx/vvvq0CBAurQoUOqz1qzZs301VdfKTk5WRs3btTOnTuD9qrITFY/x6tWrVKLFi38ia6UOhk7G3FxcakS3ZTjtGzZMtUPCvXr11d8fLy/506KwMSyRIkSioqK0nXXXaeCBQv6y6tUqaIdO3accWwrV67UH3/8oc6dO6e7Tu7Zs0fbt29XZGSkKlSooL59+2rWrFnau3fvGR8nRd26dbVo0SINHz5cq1atUnJycqr1y5YtU7169VS5cmV/LJLXCyej9w0Qikh2EXLihywMugBp3XTTTfr99981dOhQ/fHHH2rXrl26OmXKlFGRIkWCdiPOihdeeEE333yzHnvsMVWvXl1Vq1bVzJkzM6y/e/duRUVFpSsPVpa2FTql1efo0aOSvG7OLVu2lHNOL730kj755BOtXr1a0dHR/jpZkTK1UuvWrRUWFuZfKleuLEnatm2bJKlRo0aqWLGiv7XnzTffVKFChXTzzTdL8u7fTU5O1r333ptqP0WKFNGJEyf8+8ns/E4Xd/Xq1TV//nz9/PPPat26tcqWLatu3bpp3759WT5f4GwVKVJE7dq10wsvvKDvvvtOL7/8sn744YdU96yndbprxP79+/XUU0+l+syEhYXp559/TveZCdwmOTlZJUuWTLVNjx49dPLkSe3atcvf1fpMu61m9XO8e/duRUdHp9q2aNGiuuCCC87oeIGC/ei3a9euoOUxMTHpbl8Idk05m+tMMCnXyVq1aqV6Xa6//npJ3nWyQIECWrp0qWJjY9WzZ0/Fxsbqmmuu0ZdffnnGx+vZs6dGjx6t2bNnq0GDBoqJidEjjzziT3r379+vzz77LN37ZvLkyRm+b3Ldzi+DL8BfQDdmAPlWRESE2rZtq2effVadO3dWREREujphYWFq3LixlixZoscff/yMj1GqVCk9//zzev755/XNN9/o6aef1u233646der4W2IDxcbGauPGjenKzyZZW7x4sY4cOaL58+f7z+3kyZNnfP9qShfkiRMn6vLLL0+3PiXpNTPdeuutmjVrlkaPHq1Zs2bpxhtv9A+CUqpUKf8cxsG6FpcvX/6M4spImzZt1KZNGx06dEgLFy7UgAEDdP/992f6IwNwLvTq1UuJiYnasGFDhnVOd42IjIxUhw4ddPfdd6fbNqM5uiMjI1WoUCF98sknKlAgfbtGdHS0/z7RXbt2ndE5ZfVzHBsbm67lMikpKdNW7tMJ1sOmXLlyQVtI9+zZo3r16p31sc5UynXynXfeCZp8V69eXZJUo0YNvfnmmzpx4oQ++ugjDR48WG3atNH27duD/l9lpECBAho4cKAGDhyobdu2afr06Ro6dKji4uLUt29fRUZGKiEhQRMmTEi3bZEiRc7yLIG8h2QXQL7Wr18/HTt2TH379s2wzoABA3TTTTdp6tSp/vvfUpw6dUpLly7N0gBWderU0ZgxYzR9+nRt2LAhaLJ75ZVXavr06dqxY4e/K/OqVatSjSyaVUlJSSpQoIAKFfrzUj979mx/d7asql69uuLi4rR58+bTzinctWtXPfPMM3rnnXe0YsWKVPcCRkREqGHDhtq4caMeffTRMzuZIE7XAlOyZEl169ZNK1as0MqVK//y8YDM7N27N11L5r59+3To0KGgyU8wwa4RKffB1qtXL8u3UzRt2lTJyck6dOiQWrRoEbROyud66tSpQXu1SME/Y1n9HF955ZV69dVXdeTIEX9X5sAR1bNLgwYNNGHCBP3222/+H9ZWr16tzZs36+qrr87240nBX5dGjRqpaNGi2rlzZ4aDSgUKCwtT06ZNNWjQIHXr1k2//vqrIiMjz6pluUKFChoyZIgmT56s7777TpLXZXvp0qWqWLFiuvdl4HlIOquWbCCvINkFkK81adJETZo0ybROu3btNGjQIPXq1UuffPKJ2rdvrwsuuEAbNmzQiy++qPj4+AyT3auvvlodOnRQ7dq1ZWaaNGmSIiIiVL9+/aD177rrLj3++ONq27athg8frqSkJA0fPlxRUVFn9Ku/9OcfvHfddZd69eqldevW6ZlnnknXbe90ChQooLFjx+qOO+7Q4cOHdeONN6pw4cL+gVbmzJnj/2O2Xr16qlKlinr37q2iRYuqbdu2qfb19NNPq1mzZipQoIBuueUWFS9eXFu3btXChQs1atQoVatWLctx1ahRQ0uWLNGSJUtUpkwZVa5cWXPmzNHKlSvVqlUrlS9fXj/88IMWwxieAAAgAElEQVTeeOONLA/GBZytSy+9VO3bt1fLli0VHR2tLVu26JlnnlGxYsXS/UgW6HTXiBEjRqh+/fpq06aNevbsqbJly2rHjh1677331KNHj6DXr+rVq6tv377q2rWrEhMTlZCQoKNHj2rdunX6/vvv9fLLL6tAgQL+VuTbb79dt912m8xM//3vf3XbbbcpISFBNWrU0J49ezRlyhTVrl1bZcuWVXx8fJY+xwMGDNC4cePUtm1bDRo0SDt37tQTTzyhokWLZuvrPmjQIE2YMEE33HCDBg8erN9//11DhgzRpZdeqk6dOmXrsVIEu/aUKVNGI0aM0P/93/9py5Ytuvbaa3Xq1Cl9//33Wr58uebOnatvvvlGDz30kLp06aKLLrpIv/zyi5566ilddtll/pbhjPadVp8+ffyDeJUsWVLLly/XDz/8oKeeekqSdOedd+rFF19UkyZN9NBDD+miiy7SgQMHtGrVKsXGxmrgwIEqXLiwKleurNmzZ6t27doKDw9XnTp1mLsYIYVkFwCyYOzYsbrqqqv0wgsvqFu3bkpKSlJ8fLxuuukmPfTQQxlu16hRI02ZMkWbN29WwYIFdfnll+vdd9/VhRdeGLR+sWLFtHjxYvXr109dunTx/2GZmJiY6YiuwVx66aWaPHmyRo4cqblz5+qyyy7TG2+8kWrU16zq0qWLSpQoodGjR+vVV19VwYIFddFFF6lt27bp/jDq0qWLRo0apa5du6YanEby/rD/8MMPNXz4cN1xxx1KTk5WpUqV1KpVqyy3fqUYNmyYf+Csw4cPa/LkyapTp44WLFigQYMG6eDBgypXrpzuuecePfbYY2d8zsCZePTRRzV//nw98MADOnjwoGJjY/0jFqd09Q/mdNeIatWq6bPPPtOwYcPUu3dvJSUl+Uc+rlKlSob7HTdunKpVq6ZJkybp0UcfVYkSJXTJJZeoV69e/jrdunVTeHi4Ro0apVtuucXfapsyRsCtt96q5cuXKzExUfv27VP37t01ZcqULH2O4+LitGjRIj3wwAPq1KmTatasqddff13t27fPjpfbLyoqSsuXL9eDDz6o2267TYULF1br1q317LPPnrOkLdi1p0ePHkpMTFT58uX17LPPauzYsQoPD1e1atX819zY2FjFxMRo1KhR2rlzp0qVKqXrr7/en6Bmtu+0GjVqpEmTJumll17S0aNHVaVKFU2aNMk/RkJ4eLiWL1+uRx99VMOHD9eePXsUHR2t+vXrpxpx+8UXX9RDDz2k5s2b69ixY9q0aZN/bmMgFFhmc7TlRQkJCY5R5vK3jAaj2vzk6bsVIXPr169XzZo1czuMfGfTpk2qVq2aJk6ceMYjpyJrMntvm9kXzrmEHA4ppJzuu5lrC3B+y5HPaEaDUZVPP1YEkNXvZlp2AeA888QTT6h8+fKqVKmStm7dqieeeEJRUVHnrEseAABAKMrRZNfMKkh6TVKspFOSJjrn/mVmkZJmSYqXtFnSrc65X8wbieFfklpLOiKph3PufzkZMwDkNDPTyJEjtXPnThUpUkTXXHONnnnmmTPuxgwAAJCf5fQ8uyclPeicqympoaT+ZnaJpCGS3nfOVZX0vu+5JN0oqapv6S0p/fjpABBihgwZop9//llHjx7VoUOH9M4776hGjRq5HRYAAECekqPJrnNuV0rLrHPuN0nrJcVJai9pqq/aVEk3+x63l/Sa83wmqZSZndns5wAAAACAfCenW3b9zCxe0uWSPpcU45zbJXkJsaSUCcHiJG0L2Gy7rwwAAAAAgAzlSrJrZhdIelPSAOfc4cyqBilLN3y0mfU2szVmtmbfvn3ZFSaAIEJtBHeA9zQAAKEpx5NdMwuTl+hOd8695Svek9I92ffvXl/5dkkVAja/UNLOtPt0zk10ziU45xJS5oYDkP3CwsKUlJSU22EA2SopKUlhYWG5HQYAAMhmOT0as0l6RdJ659w/A1YtkNRd0pO+f+cHlN9nZjMlNZB0KKW7MyBlPKcuzo3o6Gjt2LFDcXFxKlq0qLyPNJA3OeeUlJSkHTt2KCYmJrfDAQAA2Syn59ltLOkOSd+a2Ve+sr/LS3Jnm1kvSVsldfatWyRv2qEf5U09dFfOhgsgUMrUNzt37tSJEydyORrgrwsLC1NMTEy+ntYpk2kBR0i6R1LK/UF/d84t8m3zsKRekpIlPeCcW5LjgQMAcBo5muw65z5W8PtwJalZkPpOUv9zGhSAM1KiRIl8nRgAIShlWsD/mVlxSV+Y2Xu+dc86554JrOybMrCrpFqSyktaZmbVnHPJORo1AACnkWujMQMAgNyXybSAGWkvaaZz7phzbpO83lf1z32k5y8zO+3ywQcfZMuxvvvuO40YMUK///57tuzvdBYtWqQXXnghy/U3b96sCy64QNu2/TmZxs8//6xu3bqpQoUKCg8PV8WKFdWhQwetXLnyL8UWGxub6jWOjo5Wu3bttG7dulT1evXqpf79aTsB8qOc7sYMAADOU2mmBWwsb9yMOyWtkdf6+4u8RPizgM3O3bSAO788J7s9rfKXn1H1wKQtKSlJTZs21bBhw9SmTRt/+SWXXJItoX333XcaOXKk+vbtqwsuuCBb9pmZRYsWadmyZbrvvvuyVH/kyJHq3LmzKlTwxhfdt2+fGjRooMqVK2vMmDGKiYnRpk2bNG/ePH3++edq1KjRX4qvR48e6tOnj5xz2rFjhx5//HHdcMMNWr9+vYoXLy5JSkxMVN26dTV48GBVrFjxLx0PQN5CsgsAANJNC2hmEyT9Q96Uf/+QNFZST53BtICSeksK+QSjYcOG/scpLa4XX3xxqvL84MCBA/rPf/6jZcuW+ctmzpypX3/9VUuWLFHp0qUlSddff7169uyZ6bRfixcv1s0336yjR49mesy4uLhUr/NFF12kevXqafXq1WratKkkqXr16qpXr55eeukljRo16q+cIoA8hm7MAADkc8GmBXTO7XHOJTvnTkmapD+7KjMt4F+wadMmde7cWaVKlVJERITatGmjn376yb/eOafHHntMF110kcLDwxUbG6vWrVvrwIEDWrx4sTp39sbwLFeunMxMNWrUyPBYX3/9tVq0aKHSpUvrggsuUK1atTRp0qRUdebMmaMrrrhC4eHhKl++vIYOHarkZO/26yFDhmjcuHHauHGjv6tw3759MzzejBkzFBkZqauvvtpf9uuvv6po0aIqWbJkuvrnYkT/lNbctIModurUSdOmTcv24wE4v5HsAgCQj2U0LaBv3vsUHSSt9T1eIKmrmRUxs8qSqkpalVPx5mV79+5V48aNtXnzZr388suaMWOG9u/fr5YtW+r48eOSpEmTJmns2LEaPHiwli5dqnHjxqlSpUpKSkpSo0aNNHr0aEnSwoULtXLlSs2aNSvosU6dOqU2bdooIiJC//nPfzR//nz169dPhw4d8td57bXX1KVLF11zzTVasGCBHn74YT3//PMaPny4JKl///665ZZbVKlSJa1cuVIrV67U4MGDMzy/999/Xw0bNkyVxF5xxRX67bff1KNHD3311VeZtuaeDeecTp48qZMnT2rLli0aMmSIoqOj1bhx41T1rrrqKm3btk0bN27M1uMDOL/RjRkAgPwto2kBbzOzuvK6KG+W1EeSnHPrzGy2pO/kjeTcn5GYs2bMmDE6deqUli1b5m/pbNSokSpXrqxp06apV69eWrVqldq2bas+ffr4t+vUqZP/cdWqVSV5SWRsbGyGx9q5c6d27Nih5cuX+7dp1uzPiS+Sk5M1ePBg9e7dW//6178kSS1btlTBggWVmJioxMREVahQQTExMQoPD89Sl+wvvvhC3bt3T1XWpk0b3XvvvRo/frymTZumEiVKqGXLlurfv7+aNGnir+ec87coS16yLkknT570l5mZChYsmGr/o0eP9v8AIEllypTRvHnz0t3PfNlll0mSVq1aperVq5/2XACEBlp2AQDIx5xzHzvnzDlXxzlX17cscs7d4Zy71Fd+k3NuV8A2o5xzFzvnqjvn3s3N+POSZcuWqVWrVoqIiPC3RpYuXVqXXXaZ1qxZI0mqW7eu5s2bp8cee0xr1qzxJ31nKiYmRrGxsbrnnnv0xhtvaN++fanWr127Vrt371bnzp39sZw8eVJNmzbVH3/8ofXr15/xMffs2aOyZcumK0/pCv3UU0/pmmuu0cKFC9W0aVNNmTLFX+ell15SWFiYf2nTpo2OHTuWqiwiIiLdvnv27KnVq1dr9erVWrx4sW644Qa1b98+Xfzh4eGKiIjQ7t27z/i8AORdJLsAAAA5YP/+/Zo6dWqqBC4sLEyffvqpf6qefv36afjw4Zo+fbquvPJKxcbGauTIkWec9IaFhem9995TqVKl1L17d8XGxqpJkyb69ttv/bFIXmtvYCw1a9aUpFRTB2VFcnKyTpw4oSJFigRdX61aNSUmJuqdd97Rpk2bVKtWLT388MP+9Z06dfInratXr9bzzz+vwoULpyr79NNP0+23XLlySkhIUEJCgm644QZNmzZNUVFRQQeiKlKkyGkHvAIQWujGDAAAkAMiIyPVsGHDoPe9pnRrDuxGvGXLFr322msaPny4KlWqpB49epzR8WrXrq158+bp+PHjWrFihRITE9WuXTtt3rxZkZGRkqSpU6cGnRbp4osvPqNjFSxYUCVKlNCvv/562roxMTG68847lZiYqEOHDqlkyZKKiopS4EBm+/fvl5kpISHhjOIoUKCAqlevnq5l1zmnQ4cO+c8bQP5AsgsAAJADmjVrpnfffVd16tRR4cKFT1u/UqVKeuSRR/Tyyy/ru+++kyT/dmfSQlm4cGG1aNFCDzzwgHr27Kk//vhDl156qaKiorRlyxbdeeedmW6b1WNVr15dmzZtSlW2d+9eRUdHp6v7ww8/KCIiImjX5L/i1KlTWr9+fboEfvv27UpOTla1atWy9XgAzm8kuwAAADkgMTFRM2fOVLNmzdS/f3+VK1dOu3fv1gcffKDmzZurU6dOuuuuuxQXF6f69eurRIkSWrp0qbZt26brr79ekvxTDY0fP16dOnXyTymU1qpVqzR8+HDdeuutqly5svbv36+xY8eqQYMG/gRzzJgxuueee3Tw4EG1bNlShQoV0k8//aS5c+dq0aJFKliwoGrUqKFt27Zp+vTpql69uqKjozOcN7lx48b68MMPU5VNnDhR8+bN0x133KE6dero+PHjWrx4sV5++WU9+OCDKlTor/0pumPHDn322WeSpIMHD2rq1Kn64Ycf/INupVizZo0KFiyY7+Y+BvI7kl0AAIAcEBsbq88//1xDhw7VAw88oMOHD6tcuXK69tprVbt2bUneFDmvvvqqxo0bp+PHj6tq1aqaMmWKbrzxRkneva+jR4/WhAkTNHbsWFWtWlUbNmxId6y4uDiVLl1ajz32mHbt2qXSpUurefPmevLJJ/11unfvrsjISD3xxBN66aWXVKhQIVWpUkXt2rVTgQLesC633367PvroIw0YMED79+9Xnz599OKLLwY9v44dO+pf//qXdu/e7R8pun379tqxY4defPFFbdu2TWFhYapSpYpeeukl9erV6y+/plOmTPEPdFWqVCnVrFlT8+fP979eKRYvXqwWLVr45+EFkD9Yds93ltsSEhJcyoiGCH3xQxZmue7mJ9ucw0gAhCoz+8I5d2Y3DiKV0303r1+/3j8wEvK2GjVqqH///rr//vtzOxS/EydOKC4uTuPHj9ctt9yS2+HkSTnyGd35ZfDy8pef2+MiT8rqdzOjMQMAACBbDB06VP/+97/Pesqkc2H69OmKiopSx44dczsUADmMbswAAADIFrfffru2b9+uXbt2KS4uLrfDkeSNFD1p0iR/12wA+QfJLgAAALJFgQIFUs2fez644447cjsEALmEn7gAAAAAACGHZBcAAAAAEHJIdgEAQK4LtdkhgFDBZxN5GckuAADIVWFhYUpKSsrtMAAEkZSUpLCwsNwOAzgrJLsAACBXRUdHa8eOHTpy5AitSMB5wjmnI0eOaMeOHYqOjs7tcICzwmjMAAAgV5UoUUKStHPnTp04cSKXowGQIiwsTDExMf7PKJDXkOwCAIBcV6JECf6gBgBkK7oxAwAAAABCDskuAAAAACDkkOwCAAAAAEIOyS4AAAAAIOSQ7AIAAAAAQg7JLgAAAAAg5JDsAgAAAABCDskuAAAAACDkkOwCAAAAAEIOyS4AAAAAIOQUyu0AgJwSP2RhurLNT7bJhUgAAAAAnGu07AIAAAAAQg7JLgAAAAAg5JDsAgAAAABCDskuAAAAACDkkOwCAAAAAEIOyS4AAAAAIOSQ7AIAAAAAQk6OJrtm9qqZ7TWztQFls8zsK9+y2cy+8pXHm1lSwLoXczJWAAAAAEDeVSiHjzdF0guSXkspcM51SXlsZmMlHQqo/5Nzrm6ORQcAAAAACAk5muw65z40s/hg68zMJN0qqWlOxgQAAAAACD3n0z2710ja45z7IaCsspl9aWYrzOya3AoMAAAAAJC35HQ35szcJmlGwPNdkio65w6YWT1J88yslnPucNoNzay3pN6SVLFixRwJFgAAAABw/jovWnbNrJCkjpJmpZQ554455w74Hn8h6SdJ1YJt75yb6JxLcM4lREVF5UTIAAAAAIDz2HmR7EpqLmmDc257SoGZRZlZQd/jiyRVlfRzLsUHAAAAAMhDcnrqoRmSVkqqbmbbzayXb1VXpe7CLEnXSvrGzL6WNEdSX+fcwZyLFgAAAACQV+X0aMy3ZVDeI0jZm5LePNcxAQAAAABCz/nSjRkAAAAAgGxDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQUyi3AwAAAACAoHZ+mb6s/OU5HwfyJFp2AQAAAAAhh2QXAAAAABBySHYBAAAAACGHZBcAgHzMzCqY2XIzW29m68zs/3zlkWb2npn94Pu3tK/czOx5M/vRzL4xsyty9wwAAAiOZBcAgPztpKQHnXM1JTWU1N/MLpE0RNL7zrmqkt73PZekGyVV9S29JU3I+ZABADg9kl0AAPIx59wu59z/fI9/k7ReUpyk9pKm+qpNlXSz73F7Sa85z2eSSplZuRwOGwCA0yLZBQAAkiQzi5d0uaTPJcU453ZJXkIsKdpXLU7StoDNtvvKAAA4r5DsAgAAmdkFkt6UNMA5dzizqkHKXJD99TazNWa2Zt++fdkVJgAAWUayCwBAPmdmYfIS3enOubd8xXtSuif7/t3rK98uqULA5hdK2pl2n865ic65BOdcQlRU1LkLHgCADJDsAgCQj5mZSXpF0nrn3D8DVi2Q1N33uLuk+QHld/pGZW4o6VBKd2cAAM4nhXI7AAAAkKsaS7pD0rdm9pWv7O+SnpQ028x6SdoqqbNv3SJJrSX9KOmIpLtyNlwAALKGZBd5RvyQhbkdAgCEHOfcxwp+H64kNQtS30nqf06DAgAgG9CNGQAAAAAQckh2AQAAAAAhh2QXAAAAABBySHYBAAAAACGHZBcAAAAAEHIYjRkAAAA4H+z8Mn1Z+ctzPg4gRORoy66ZvWpme81sbUDZCDPbYWZf+ZbWAeseNrMfzWyjmd2Qk7ECAAAAAPKunO7GPEVSqyDlzzrn6vqWRZJkZpdI6iqplm+b8WZWMMciBQAAAADkWTma7DrnPpR0MIvV20ua6Zw75pzbJOlHSfXPWXAAAAAAgJBxvgxQdZ+ZfePr5lzaVxYnaVtAne2+MgAAAAAAMnU+JLsTJF0sqa6kXZLG+sotSF0XbAdm1tvM1pjZmn379p2bKAEAAAAAeUauJ7vOuT3OuWTn3ClJk/RnV+XtkioEVL1Q0s4M9jHROZfgnEuIioo6twEDAAAAAM57uZ7smlm5gKcdJKWM1LxAUlczK2JmlSVVlbQqp+MDAAAAAOQ9OTrPrpnNkNREUlkz2y5puKQmZlZXXhflzZL6SJJzbp2ZzZb0naSTkvo755JzMl4AAAAAQN6Uo8muc+62IMWvZFJ/lKRR5y4iAAAAAEAoyvVuzAAAAAAAZLccbdkFAAAAkIN2fhm8vPzlORsHkAto2QUAAAAAhBySXQAAAABAyCHZBQAAAACEHJJdAAAAAEDIIdkFAAAAAIQcRmNGvhY/ZGG6ss1PtsmFSIDsF+z9LfEeBwAA+QMtuwAAAACAkEPLLgAAAIC8jfmEEQQtuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJBDsgsAAAAACDkkuwAAAACAkEOyCwAAAAAIOSS7AAAAAICQQ7ILAAAAAAg5hXI7AAAAACDf2fllbkcAhDxadgEAAAAAIYdkFwAAAAAQckh2AQAAAAAhh2QXAAAAABBySHYBAAAAACGHZBcAAAAAEHJIdgEAAAAAIYdkFwAAAAAQckh2AQAAAAAhh2QXAAAAABByCuV2AACAvy5+yMLcDgEAgJyx88vcjgB5BC27AAAAAICQQ7ILAAAAAAg5JLsAAAAAgJCTo8mumb1qZnvNbG1A2Rgz22Bm35jZXDMr5SuPN7MkM/vKt7yYk7ECAAAAAPKunG7ZnSKpVZqy9yTVds7VkfS9pIcD1v3knKvrW/rmUIwAAAAAgDwuR5Nd59yHkg6mKVvqnDvpe/qZpAtzMiYAAAAAQOg53+7Z7Snp3YDnlc3sSzNbYWbX5FZQAAAAAIC85byZZ9fMhko6KWm6r2iXpIrOuQNmVk/SPDOr5Zw7HGTb3pJ6S1LFihVzKmQAAAAAwHnqvGjZNbPuktpKut055yTJOXfMOXfA9/gLST9JqhZse+fcROdcgnMuISoqKqfCBgAAAACcp866ZdfMSkuqJGm9c+7YX9hPK0mDJV3nnDsSUB4l6aBzLtnMLpJUVdLPZ3sc5C3xQxbmdggAkOdk13czAAChIEstu2Y20syeDHjeVNJWSV9I+snMamVxPzMkrZRU3cy2m1kvSS9IKi7pvTRTDF0r6Rsz+1rSHEl9nXMHg+4YAIB8Jru+mwEACFVZ7cZ8u6QNAc/HSvpYUmNJGyU9kZWdOOduc86Vc86FOecudM694pyr4pyrkHaKIefcm865Ws65y5xzVzjn3j6D8wIAINRly3czAAChKqvJbnn5uhCbWQVJl0ka7pz7TNI/JTU8N+EBAIAM8N0MAEAmsprs/iappO9xU0m/OOdW+Z4flVQsuwMDAACZ4rsZAIBMZHWAqhWShpjZKUkPSZofsK6apG3ZHRgAAMgU380AAGQiq8nuQEnTJM2U9JWkoQHr7pT0YTbHBQDIAKOVw4fvZgAAMpGlZNc5t0NeF6lgbpCUlG0RAQCA0+K7GQCAzGV16qH/mlmNDFbHSlqSfSEBAIDT4bsZAIDMZXWAqiaSSmSwroS8OXEBAEDOaaJs+G42s1fNbK+ZrQ0oG2FmO8zsK9/SOmDdw2b2o5ltNLMb/soJAABwLmU12ZUkl7bAzArL60K1O9siAgAAWZUd381TJLUKUv6sc66ub1nk2/clkrpKquXbZryZFTybwAEAONcyTHbNbLiZJZtZsrwv089SngeUJ8mbtP71HIoXAIB861x8NzvnPpR0MIshtJc00zl3zDm3SdKPkuqf+ZkAAHDuZTZA1SJJ+yWZpOcljZW0OU2d45I2OOc+OifRAQCAQDn53Xyfmd0paY2kB51zv0iKk/RZQJ3tvjIAAM47GSa7zrnVklZLkpn9Jmmhc25/TgUGAABSy8Hv5gmS/iGv9fgf8pLqnvKS7HRhBduBmfWW1FuSKlaseA5CBAAgc1mdemjquQ4EAPD/27vzuMmq8sDjv0cQFEQBabDZbFE00SQCtgRjTFQSBMkIqCBohDEqJiNMzGiSHiZRszhpE9foDAaVLRNQjBKRxqXTYMwiCoqgSFTAFjvdLMoiAZXtmT/OKaiurqr3vt1vbff9fT+f+tRb595b9Zy6972nzj3LHY9+9+ldu/LwCUSiLTHKsjkzb+r8HREfBC6sL9cBe3WtuiewfsB7nAacBrB8+fK+FWJJkkapUWU3Ih4O/C7wYkrB9ojedTJz14UNTZIkDTLKsjkilmbmhvryKKAzU/MFwDkR8S5gd2Bf4Mub8xmSJI1ao8ou8G7gdZQru5dQxgNJkqTJWZCyOSLOpdzGaJeIWAe8BXhuROxH6aK8tn4OmXl1RJwHfBO4D3h9Zt6/ZdmQJGk0mlZ2jwZWZOY7RxmMJElqbEHK5sw8rk/yh4es/zbgbVvymZIkjUPT++wGcNUoA5EkSfNi2SxJ0hBNK7sfBPpd+ZUkSZNh2SxJ0hBNuzHfBLwiIi4BVgO39yzPzDx1QSOTJEnDWDZLkjRE08rue+rz3sCv9lmelHvySZKk8bBsliRpiKb32W3a3VmSJI2BZbMkScNZUEqSJEmSWqdxZTcido2It0fEmoj4dkQ8rab/bkQ8a3QhSpKkfiybJUkarFFlNyIOBL4DvIRyc/knAtvWxUuBN44iOEmS1J9lsyRJwzWdoOrdwCXAiykV5Fd1Lfsy8PIFjkuamGUrVvVNX7vy8DFHIklDWTZLkjRE08ruAcARmflARETPsh8Cuy5sWJIkaQ6WzZIkDdF0zO4dwJIBy/ah3OtPkiSNj2WzJElDNK3sfhL4k4jYpystI2IX4E3AJxY8MkmSNIxlsyRJQzTtxrwCWAN8E7n7uGsAACAASURBVPhKTfsA8CTgu8CbFz40SZI0hGWzJM1l/RWbpu2+//jj0EQ0quxm5m0RcRDwSuBg4C7gVuBDwNmZ+dPRhShJknpZNkuSNFzTll0y8x7gw/UhLYhBMx9LkuZm2SxJ0mBN77P7TxHxOxExaCIMSZI0RpbNkiQN13SCqluAdwDrI2J1RPxWROw0wrgkSdJwls2SJA3RqLKbmS+l3K/veOA/gf8D3BgRF0bEKyNihxHGKEmSelg2S5I0XNOWXTLzrsw8NzOPohSur62LPgjcOIrgJEnSYJbNkiQN1riy2y0z7wSuo9za4EfAIxcyKEmSND+WzZIkbazxbMwAEXEg8DLgaGAP4GrgvcBHFj40SZI0F8tmaQb0u9erpJFrVNmNiJXAMcDjge8AZwAfycxrRhibJEkawLJZkqThmrbsHgOcRylEvzbCeCRJUjOWzZIkDdGospuZ+4w6EEmS1JxlsyRJwzWeoCoitq03r/9wRHwuIvat6S+LiJ8dXYiSJKkfy2ZJkgZrOmb3ycBq4DHAV4DnAp379z0HOJxynz9JkjQGls2SJA3XtGX3r4EbgGXAC4DoWvZPwC8vbFiSJGkOls2SJA3RtLL7HOAvMvN2IHuW3QQsbfqBEXF6RNwcEd/oSts5IlZHxHfq8041PSLiryPi2oi4KiIOaPo5kiS13IKVzZIktVHTyu5PGHxz+j2A2+fxmWcCh/akrQDWZOa+wJr6GuAwYN/6OBE4dR6fI0lSmy1k2SxJUus0reyuBk6JiMd0pWVEbAucDFzU9AMz8wvArT3JRwBn1b/PAo7sSj87i0uBHSPCK9WSJC1g2SxJUhs1vc/u7wP/ClxLKVwTeDPwNGAb4MVbGMdumbkBIDM3RMSuNX0P4Ptd662raRu28PMkSZp1oy6bJUmaaY1adjPz+8DTgQ9QJsK4jjIW6GPAMzLzxhHFF33SesclEREnRsTlEXH5LbfcMqJQJEmaHhMsmyVJmglNW3bJzNuAP66PhXZTRCytrbpLgZtr+jpgr6719gTW94ntNOA0gOXLl29SGZYkqY1GXDZLkjTTmo7ZHbULgBPq3ycAn+xKP77OynwQcEenu7MkSZIkSYM0btldKBFxLuXG97tExDrgLcBK4LyIeDXlnoFH19UvAl5IGY90N/CqcccrSZIkSZo9Y6/sZuZxAxYd3GfdBF4/2ogkSZIkSW0zLd2YJUmSJElaMFZ2JUmSJEmt07iyGxHHR8SOowxGkiQ1Z9ksSdJg8xmzewawP3B7RATlNgeneR8/NbVsxapJhyBJbWPZLEnSAAMruxGxCriyPq4CAujcw/ZhlFmULwQsUCVJGgPLZkmSmhvWsruacrX4N4CfoRSm74+IS4DL2LiAlSRJo2fZLElSQwMru5n5ns7fEbEt8GPgq8BTgFdSCtO/jYjPAP+YmZ8ZcaySJC1qls2SJDU3rBvzycAVwJWZeWcZCsQZmXlVRGwN3AOcC+wFvB940hjilaRWGjSmfe3Kw8cciaaZZbMkSc0N68b8IuCPgF0iYi3lavGxEfFI4Ot1nU9n5ldHG6IkSaosmyVNr/VXbJq2+/7jj0OqhnVj/nWAiFgK7AesAn4N+G1gO0oB+zsR8RHgXzLzp6MPV5I0Cv1alm1Vnj6WzZIkNTfnfXYzc0Nmfrq+fE1m7gwsp0yCsRdwJnDbyCKUJEkbsWyWJGluc1Z2B7imPp+SmXsBz1igeCRJ0uaxbJYkqcuwMbsbyczuinEC3wN+Wpdd03cjSZI0MpbNkiQN1riy2y0zHwCesMCxSJKkzWTZLEkNOZHWorG53ZglSZIkSZpaVnYlSZIkSa2zWd2YJUmSJPXo1z1W0sTYsitJkiRJah0ru5IkSZKk1rEbsyRNsWUrVk06BEmSpJlky64kSZIkqXWs7EqSJEmSWsfKriRJkiSpdazsSpIkSZJax8quJEmSJKl1rOxKkiRJklrHyq4kSZIkqXWs7EqSJEmSWsfKriRJkiSpdazsSpIkSZJax8quJEmSJKl1rOxKkiRJklrHyq4kSZIkqXWs7EqSJEmSWsfKriRJi1xEnB4RN0fEN7rSdo6I1RHxnfq8U02PiPjriLg2Iq6KiAMmF7kkSYNZ2ZUkSWcCh/akrQDWZOa+wJr6GuAwYN/6OBE4dUwxStLorL9i04dmnpVdSZIWucz8AnBrT/IRwFn177OAI7vSz87iUmDHiFg6nkglSWrOyq4kSepnt8zcAFCfd63pewDf71pvXU2TJGmqWNmVJEnzEX3ScpOVIk6MiMsj4vJbbrllDGFJkrQxK7uSJKmfmzrdk+vzzTV9HbBX13p7Aut7N87M0zJzeWYuX7JkyciDlSSpl5VdSZLUzwXACfXvE4BPdqUfX2dlPgi4o9PdWZKkabL1pAMAiIinAB/tStoHeDOwI/BaoNP/6ZTMvGjM4UkALFuxapO0tSsPn0AkkrSwIuJc4LnALhGxDngLsBI4LyJeDdwAHF1Xvwh4IXAtcDfwqrEHLElSA1NR2c3MbwH7AUTEVsB/AOdTCtB3Z+Y7JhieJEmtlpnHDVh0cJ91E3j9aCOSJGnLTUVlt8fBwHWZ+b2IfnNgSJIkSdqE94aVNjKNY3aPBc7ten1SRFwVEadHxE6TCkqSJEmSNDumqrIbEdsALwI+VpNOBZ5I6eK8AXjngO28vYEkSZIk6UFTVdkFDgO+mpk3AWTmTZl5f2Y+AHwQOLDfRt7eQJIkSZLUbdoqu8fR1YW5c3+/6ijgG2OPSJIkSZI0c6ZmgqqI2A74deB1Xcl/GRH7AQms7VkmSZIkSVJfU1PZzcy7gcf2pL1yQuFIkiRJkmbYtHVjliRJkiRpi1nZlSRJkiS1jpVdSZIkSVLrWNmVJEmSJLWOlV1JkiRJUutMzWzMkqTpsmzFqr7pa1cePuZIJEmS5s+WXUmSJElS61jZlSRJkiS1jpVdSZIkSVLrWNmVJEmSJLWOlV1JkiRJUutY2ZUkSZIktY6VXUmSJElS63ifXUmSJEnjs/6KSUegRcKWXUmSJElS61jZlSRJkiS1jpVdSZIkSVLrWNmVJEmSJLWOlV1JkiRJUutY2ZUkSZIktY6VXUmSJElS61jZlSRJkiS1jpVdSZIkSVLrbD3pAKRZtmzFqk3S1q48fAKRSJIkaUGtv6J/+u77jzcObTZbdiVJkiRJrWNlV5IkSZLUOlZ2JUmSJEmt45hdNdZvfCo4RlWSJEnS9LFlV5IkSZLUOlZ2JUmSJEmtY2VXkiRJktQ6VnYlSZIkSa1jZVeSJEmS1DrOxixJYzZoZnNJkiQtHCu7kiRJ0nytv2LSEUiag5VdjYQtV5IkSZImyTG7kiRJkqTWsbIrSZIkSWodK7uSJEmSpNaxsitJkiRJah0ru5IkSZKk1rGyK0mSJElqnam69VBErAXuBO4H7svM5RGxM/BRYBmwFjgmM2+bVIySNB/ehkuStKh5P2JN0DS27D4vM/fLzOX19QpgTWbuC6ypryVJkiRJGmgaK7u9jgDOqn+fBRw5wVgkSZIkSTNg2iq7CXwuIr4SESfWtN0ycwNAfd51YtFJkiRJkmbCVI3ZBZ6dmesjYldgdUT8e5ONasX4RIC99957lPFJ0qLXbxzy2pWHTyASSZKkwaaqZTcz19fnm4HzgQOBmyJiKUB9vrnPdqdl5vLMXL5kyZJxhixJkiRJmkJTU9mNiO0jYofO38AhwDeAC4AT6monAJ+cTISSJEmSpFkxTd2YdwPOjwgocZ2TmZ+JiMuA8yLi1cANwNETjFGSJEmSNAOmprKbmdcDT++T/kPg4PFHJEmSJEk9+t07ePf9xx+H5jQ13ZglSZIkSVooVnYlSZIkSa1jZVeSJEmS1DpWdiVJkiRJrWNlV5IkSZLUOlMzG7MkzbplK1ZNOgRJkiRVtuxKkiRJklrHll1tMVuzJEmSJE0bW3YlSZIkSa1jZVeSJEmS1DpWdiVJkiRJrWNlV5IkSZLUOk5QJUmS+oqItcCdwP3AfZm5PCJ2Bj4KLAPWAsdk5m2TilGSpEFs2ZUkScM8LzP3y8zl9fUKYE1m7gusqa8lSZo6VnYlSdJ8HAGcVf8+CzhygrFIkjSQlV1JkjRIAp+LiK9ExIk1bbfM3ABQn3edWHSSJA3hmF1JkjTIszNzfUTsCqyOiH9vumGtHJ8IsPfee48qPkmSBrJlV5Ik9ZWZ6+vzzcD5wIHATRGxFKA+3zxg29Myc3lmLl+yZMm4QpYk6UG27EqSpE1ExPbAwzLzzvr3IcCfAhcAJwAr6/MnJxelJE2J9Vdsmrb7/uOPQxuxsitJkvrZDTg/IqD8XjgnMz8TEZcB50XEq4EbgKMnGKMkSQNZ2ZUkSZvIzOuBp/dJ/yFw8PgjkiRpfqzsSmOwbMWqvulrVx4+5kgkSZKkxcEJqiRJkiRJrWNlV5IkSZLUOlZ2JUmSJEmtY2VXkiRJktQ6VnYlSZIkSa3jbMzqa9DswZLUVL/ziDOQS5KkcbFlV5IkSZLUOrbsSpK2mL1BJEnStLFlV5IkSZLUOrbsLiKDWl4cQ7ewbOGSJKll1l8x6Qi22PsuvhaAk5//pAlHIo2PLbuSJEmSpNaxZVeaIGerlSRJkkbDll1JkiRJUuvYsitJQ9j6LkmaZZ2xutJiZMuuJEmSJKl1rOxKkiRJklrHyq4kSZIkqXWs7EqSJEmSWsfKriRJkiSpdaZiNuaI2As4G3gc8ABwWma+NyLeCrwWuKWuekpmXjSZKNur32yzUpt5zEuSJLXfVFR2gfuAN2bmVyNiB+ArEbG6Lnt3Zr5jgrFJkiRJ0vysv2LTtN33H38ci9hUVHYzcwOwof59Z0RcA+wx2agkSZKkKdWvIqXpN2i/WQkeiakbsxsRy4D9gS/VpJMi4qqIOD0idppYYJIkSZKkmTEVLbsdEfEo4OPAGzLzRxFxKvBnQNbndwK/1We7E4ETAfbee+/xBSyNQL/xpGtXHj6BSCRJUmvZxVaLwNS07EbEwykV3b/LzE8AZOZNmXl/Zj4AfBA4sN+2mXlaZi7PzOVLliwZX9CSJEmSpKk0FS27ERHAh4FrMvNdXelL63hegKOAb0wiPkmzwVbxdnF/SpKkLTEVlV3g2cArga9HxNdq2inAcRGxH6Ub81rgdZMJT5IkSZI0S6aispuZ/wJEn0XeU7cBWz8kSZIkaWNTM2ZXkiRJkqSFYmVXkiRJktQ6VnYlSZIkSa0zFWN2JQ3Xb1z2II7XHr357A9tbNB353ErSZIWmi27kiRJkqTWsbIrSZIkSWodK7uSJEmSpNZxzK6kVnN8rSRJ0uJkZVeSJElqmfddfO3o3nz9FaN7b2kBWdmVNBH9WlydkVeSJEkLxTG7kiRJkqTWsbIrSZIkSWodK7uSJEnSIvG+i68d7XheaYo4ZrelnIFWbecxvjg51luSJDVly64kSZIkqXVs2ZUkSZKkadPvFk+77z/+OGaYLbuSJEmSpNaxZXfGOE5RUht5bpOkAfq17s2QzmRYJz//SROORIuRLbuSJEmSpNaxZVfS1HCmXUmS2sHbG83TjLfgTysru5IkSZIWhJVcTRMru5KkmTZovK+9AiRJWtwcsytJkiRJah0ru5IkSZKk1rGyK0mSJElqHcfsSlowjp2UJKllnCVYM8yWXUmSJElS69iyO2beR1SSJEnSZhnU0r77/uONY0bYsitJkiRJah1bdqfUoLGP0kKZz/jaSR6P/i9IkiRpc1jZlSRJklrifRdfO+kQ+urEdfLznzThSLSYWNldAM5AK0mS1BIzOvvwtFZyNSb9jlvH8TpmV5IkSZLUPrbsTgHHJGoheTxJkrR4TEuL7qjjsBu0Noctu5IkSZKk1rGyK0mSJElqHbsxS5IkafGZ0YmoFtrmdg+eVLdiuzNrPmzZlSRJkiS1ji27Q/Sb6MfbCantRjHBlZNmSZIkadys7EqSJEmLTN/Zkzeja/dCz8Lc2015Wmab1myysjtPtlBJkiRJ0vSzsitJkqR2m9LJqHpbLQe1ZvabjGlWWzwHxb25E2Rt6fu0WtPjfvf9RxvHBM1EZTciDgXeC2wFfCgzV044JEmSFjXLZqm5pl1zZ7UCO03m+q4HVYZHNcvzTMwePahS3IJK8NTPxhwRWwH/BzgMeCpwXEQ8dbJRSZK0eFk2S5JmwSy07B4IXJuZ1wNExEeAI4BvTjSqBhzfK0lqqcmVzf1aIGa59WE+LSqTzvskP3/KWp42t/vxLLXcjnuiqN73H9XnbWm+NrfleCZtSTfoKfmfnfqWXWAP4Ptdr9fVNEmSNBmWzZKkqReZOekYhoqIo4EXZOZr6utXAgdm5sld65wInFhfPgX41tgDXTi7AD+YdBBbYNbjB/MwLWY9D7MeP5iHjsdn5pKFCKYtZqBsbsOx29Riyissrvya13YyrwujUdk8C92Y1wF7db3eE1jfvUJmngacNs6gRiUiLs/M5ZOOY3PNevxgHqbFrOdh1uMH86ChprpsXkz7fTHlFRZXfs1rO5nX8ZqFbsyXAftGxBMiYhvgWOCCCcckSdJiZtksSZp6U9+ym5n3RcRJwGcptzc4PTOvnnBYkiQtWpbNkqRZMPWVXYDMvAi4aNJxjMmsd8ee9fjBPEyLWc/DrMcP5kFDTHnZvJj2+2LKKyyu/JrXdjKvYzT1E1RJkiRJkjRfszBmV5IkSZKkebGyOyIRcXpE3BwR3xiw/DER8amIuDIiro6IV9X0/SLiizXtqoh4Wdc2Z0bEdyPia/Wx3zTmoS67vyvOC7rSnxARX4qI70TER+vEJlOXh4h4Xlf8X4uIn0TEkXXZtO2HnSLi/Hq8fDkifq5r2aER8a2IuDYiVnSlj20/bG78EbFXRFwSEdfUffO7Xdu8NSL+o2sfvHBU8W9JHuqytRHx9Rrn5V3pO0fE6roPVkfETtOYh4h4Ss//wo8i4g112dj2w7DjoWudiIi/rsf7VRFxQNeyE+p3/Z2IOKEr/Rl1/1xbt41R5UELq8Ex/Yp6HFwVEf8WEU8fd4wLaa78dq33zChl8EvHFdtCa5LXiHhuPe9cHRH/NM74FlKD43jgb61Zs6Xn8VnSMK+tOEc1yWvXuuM/P2WmjxE8gF8BDgC+MWD5KcDb699LgFuBbYAnA/vW9N2BDcCO9fWZwEunPQ/19X8O2OY84Nj69weA35nWPHSts3NN325K98NfAW+pf/8MsKb+vRVwHbBPPbauBJ467v2wBfEvBQ6of+8AfLsr/rcCb5r2fVBfrwV26bPNXwIr6t8rOsfhNOaha52tgBsp97Yb634Ydjx0rfNC4NNAAAcBX6rpOwPX1+ed6t871WVfBp5Vt/k0cNi4jisfW3xMzHVM/1LXfj6sczzM6mOu/NZ1tgIupoylHls5NYF9uyPwTWDv+nrXScc8wrzO+TtlVh5bch6ftUfDvLbiHNUkr3XZRM5PtuyOSGZ+gXJCGrgKsENtRXhUXfe+zPx2Zn6nvsd64GbKyW3sNjcPg1au6z0f+PuadBZw5MJEOyDAhcnDS4FPZ+bdo4lyuAZ5eCqwpq7778CyiNgNOBC4NjOvz8x7gI8AR4x7P2xu/Jm5ITO/WtPvBK4B9hhVnMNswT4Y5gjKdw/T8b/QJA8HA9dl5vdGE+VgDY+HI4Czs7gU2DEilgIvAFZn5q2ZeRuwGji0Lnt0Zn4xS0l8NiPeD1o4cx3TmflvdX8DXEq5F/DMavA/DHAy8HHKb4eZ1SCvLwc+kZk31PVnNr8L/Vtrmm3heXymNMlrW85R8/i9NpHzk5XdyXk/8LPAeuDrwO9m5gPdK0TEgZQWueu6kt9Wuzu8OyK2HVu0/Q3LwyMi4vKIuDRq91/gscDtmdk5Sa9jQpWXLnPuB8r9I8/tSZum/XAl8GJ48Jh5POWEuQfw/a71Ot/3tO2HQfE/KCKWAfsDX+pKPqnug9NjxF2AGxiWhwQ+FxFfiYgTu7bZLTM3QCkogF3HGG8/c+4H+v8vjH0/DDgeYPAxPyx9XZ90tc+rKa1FrRURewBHUXrrtN2TgZ0i4vP13Hr8pAMaoSa/U2bOZpzHZ9aQvHZrxTlqUF4neX6ysjs5LwC+RumqvB/w/oh4dGdhvYr1t8Cruk5q/5PSvfCZlO54fzjWiDc1LA97Z+ZyytXX90TEEyldUnpNejrwJvvh5yn3kuyYtv2wklLof41y1ewKylXfQd/3tO2HQfEDEBGPolwJfENm/qgmnwo8kbLPNgDvHGvEmxqWh2dn5gGULkqvj4hfmVCMc5lrP2wDvAj4WNc2Y98PA46HBxf32WTYMT9t/wsagYh4HuWH5KTP1aP2HuAPM/P+SQcyBlsDzwAOp5TjfxwRT55sSCMz9HfKLNrM8/hMmiOvnXVacY6aI68TOz/NxH12W+pVwMrade7aiPgupQL15XoSWwX8Ue3CATzY+gPw04g4A3jTuIPuMTAPtQs2mXl9RHyecpXn45TuKFvXVsU9KVcqJ2lgHuryY4DzM/PezgbTth/qCaUzsVYA362P7YC9ulbtfN8/YIr2w5D4iYiHU46bv8vMT3Rtc1Pn74j4IHDhOGPuNSwPXf8LN0fE+ZTu5V8AboqIpZm5oV5UmWg3vGF5qA4Dvtr93Y97Pww6Hrqso/8xvw54bk/652v6nn3WV0tExC8AH6KMxf7hpOMZseXAR8q/L7sAL4yI+zLzHyYb1kisA36QmXcBd0XEF4CnU8YKts1cv1Nmyhacx2dOg7y25hzVIK8TOz/Zsjs5N1DGv1HHxT0FuL62npxPGa/Q3YLSaWXs/BA9Ehg6I+MYDMrDTp2uvRGxC/Bs4Jv1RH0JZQwswAnAJ8ce9cb65qFr+XH0dNuctv0QETvGQ7Mpvwb4Qq24XAbsG2Xm5W0oXVAvmLb9MCj++v1+GLgmM9/Vs033+J2jmNJ9EBHbR8QOdZ3tgUN4KNYLKN89TMH/wpDjqGPg/0I10v0w7HjocgFwfBQHAXfUi1OfBQ6p56adKPvhs3XZnRFxUH3/45n8OUkLJCL2Bj4BvDIz21gJ2khmPiEzl2XmMsqcDP+tpRVdKP+nz4mIrSNiO+AXKeME22iu3ykzYwvP4zOlSV7bco5qktdJnp+i/O7VQouIcyktCbsANwFvAR4OkJkfiIjdKbP6LqV02ViZmf8vIn4TOAO4uuvt/mtmfi0iLqZMVhWULi2/nZn/OYV5+CXgb4AHKBdU3pOZH67vuQ9loqSdKd0kfzMzfzpteajbLgP+Fdire3zMFO6HZ1Em1rmfMjvlqzsTHkS5Fcx7KDPgnZ6Zb6vpY9sPmxt/RPwy8M+UMUqd7/+UzLwoIv6W0p0rKbMdv26UheEW5GEfysUrKD1pzunaB4+lzIq9N+XHzNGZOdfkM2PPQ912O8oYqn0y846u9xzbfhh0PFC+v04egjK+7VDgbsowkMvr9r9V1wd4W2aeUdOXU84Bj6SMlzo5LRhnQoNj+kPAS4DOhGr31eE1M2mu/PaseyZwYWb+PTOoSV4j4vcprZ4PAB/KzPdMJNgttCW/U2bNlp7HZ0nDvLbiHNUkrz3rn8kYz09WdiVJkiRJrWM3ZkmSJElS61jZlSRJkiS1jpVdSZIkSVLrWNmVJEmSJLWOlV1JkiRJUutY2dW8RcRbIyK7HjdGxIX1xtgzKyK2qXnbb8yfu2v93GUL+J7viIi1c6zTvR8fiIjbIuKyiHhbRDxuoWIZ8vlr62f/rz7LntMV27JRx9JERLwxIi7pk/6rEfHJiLg5Iu6tz5+KiMPqLRSavv+FEfH1IcvfX/fRthFxdER8KyK22tz8SGoXy+YF/1zL5k2XWTZvutyyecpZ2dXmugN4Vn28AXgysDoidp5oVFtmG8r97cZaoAK71s9dNubPhYf24y8Bx1Jvbg58PSKeMYbP/0/guD7px9ZlUyEiHgX8IbCyJ/0NwCWUe9OeDBwMnESJ/ULg+fP4mHOBn4uIp/X5/K2AlwKfqPdD/jjlfouvnHdmJLWZZfPCsWzelGXzxp9j2TwDrOxqc92XmZfWx0eA4ykFw6ETjmssIuKRk45hgXTvx89m5l8AvwBsAD46hquTFwJPjYif6yR0FR4XjPiz5+M44KfA5zoJEXEA8A7gTzPzxZn50cz8Qmael5nHAb8M/GAen/FJ4G7Kj4lezwN2oxS6ZOYDwNmUQlySOiyb28GyuRnLZs3Jyq4WypX1ea/uxIjYOSL+JiJuioifRMS/RcQv9qyzVUT8z4j4dkT8NCLWRcSZPeucFBHfqcuvjYjf61n+1oj4QUTsHxGXRsTdEXFFRDynZ70XRcRXIuKu2u3kSxHxq3XxnfX5jO5uOvWREfGKiDg7Im4HPlXfLyPipH6x9KQ9PiLOrTHeHRFXRcTLazegTveYSzqfO8/vb8eIOKfmaUO/rkfzkZm3A38APBH49SbbRMTTIuIzEXFrjeOaiHh9g03/A/gXNi5Eng88ij4Fau2udFlE3FG/k09FxJN61vnliPjniPhRfXwtIo7uWj7sGBjkBMqV2+xKOxm4Gfjzfhtk5hcz88rutIh4TURcXY/j70XEH3St37ni/LI+b3cscBPlSnXHx4EDun+MSFIPy+aeWHrSLJv7s2x+aH3L5hlnZVcLZe/6/N1OQkRsC/wj5aT8+8CRwC3AP8bG407+BvgT4DzgN4A3Att3vc9rgfdRTrD/BfgY8M6IWNETw3bAWfX9XkK52nd+RGxX3+eJwN8DF9f3eQXlBNbp3tXp1vLnPNQNbEPX+7+DUugeDfzvRt9K+dxdgS8CzwTeVD/7w5QfHxtqHACv7/rc+Xx/ZwCHUbqsnQgcQv8rkPNxCXAfcFDD9S+gdBf6TeBFlP21Q8Ntz2XjeI+j/GC5q8+6ewLvB44AXgtsBfxrRDwGICIeTdmn11OOgZcCfwvsWJfPdQxsIiK2B34R+LeeRb8CXJyZ9zXJZET8PnAq8A+UNAcLpAAACGNJREFU4/xU4M96fpCdC+wbXd3UIuLhwFHAeZl5fyc9M68BbqN0z5KkfiybB7BsnpNl80Msm2dZZvrwMa8H8FZKF5Ct6+OJwGrgCmDbrvVeDdwD7NuVtjVwHfBX9fXPAAn89wGf9TDKFcYzetL/L2VMyyO6Ykrg+V3r7FfTDq2vXwr8cEi+HlXX/6896ctq+vl9tkngpH7fT9frv6AUDksHfO7P1fd5bk96k+/vaXXbl/Xk41ZgbZP9OGT5BuDUBsfDLjWGn5/ncbSW8iNlCXAv5QfHNpRC4khKoZPAsgHbbwU8kvIj5/iatrxus8OAbYYeAwO2+aX6nk/rSf8x8Bc9adH1f7E18LCa/mjKWKG39Kz/p8CNwFb19bY1/3/VtU7ne3hWn9g+D/zdfPLjw4ePdj6wbO5eZtls2dydZtm8iB+27GpzPZZyErwXuBbYH3hxlgH6Hb8GfAX4bkRsHRFb1/R/opz4oIx3ADhzwOfsCexOuWLc7aOUk9TPd6XdSznBdHyz6z2gdEl6TEScFRGH1KuC87Fqnut3PB/4TGZumHPNjTX5/p5Znx/sVpSly83qzYy1W9PZCm8Fvg98ICJeVq+WN5aZt1Cu5h5LGVcWwKf7BhRxUESsjogfUq5u3035AfHkusp1lILrnIg4IiJ27HmLzTkGOlfq+43xyZ7XL+Gh/4t7gb+s6c+itIh8rLMv6/68mDLeZ0+A+v9zPnBMxIOzRb4M+B5waZ/P/0FXfJJk2dycZfMQls2WzW1hZVeb6w7Kyfwg4HWUq37nRET3MbVLXX5vz+NVPDR+6LHAXZn5owGfs7Q+39ST3nnd3cXlR1kmBwAgM++pfz6ivv4WpYvNPsBFwA/qeJolc+a2fwxNPZaNu1w11eT7exxwZ2b+uGfbmzcv1CIiHkGJe8481+/8EMpV0NOBG+u4nP3n8ZEfAY4BXg78Q88Ps05Me1MmoQjKMfdsyjF4Mw/t49tqLA+ndL27JSJWRcQ+dfnmHAOPqM+9Ma3noR9rHWtqTM9k432+S32+mo33ZWecT/d4unMpXQ+fVffDEcC5WS8X9/hpV3ySZNncnGXz3CybH2LZPKO2nnsVqa/7MvPy+veXIuLHlBnojqZc2YVyVfFy4Hf6bN85Of0Q2D4iHj2gUO2clHqvSO7W9RmNZeYqYFUdR3I48B7KGJYm42gGndC26UnrHWPyQx76YTAfTb6/G4EdIuKRPYXqvK7g9vE8yvnhi01Wzsx/B15Sx7A8B3g75Xves/tHzhCfAD5AOX4OH7DOoZSxX0dk5l0A9QrsRt93Zn4RODTKrJy/BrwLOIc6xmkzjoHOMbYjcHtX+heAQyJiq6zjdWqBfnmN7Z4+7/Eb9P+R8q2uvy+u6xxLOW52oM702MeOzPN/QFKrWTYXls1YNls2C2zZ1cL5f5QrY3/YlbYGeBJwQ2Ze3vPozHJ4cX0+fsD7rqNcpTu6J/0Y4Ec8NFvivGTmHZl5DqVbylNr8kZXmxtaB/xs50W9et57/7Y1wAsiYjf6G/S5Tb6/y+rzi7pieBQNZ2rsp3YvejulC9w/zmfbzLw3My+mFGJLqZNPNNjujvqZHx/ymY8EHqB0keo4hgEX7TLzx5n5KcoV7af2Wd7vGOinU9g9oSf9fZQfdqcM2bbji5RxRLv32ZeXZ2ZntlFq4fwxyjH/cuCazLxqwPsuA77d4PMlLU6WzVg2WzYPZNm8CNiyqwWRmRkR/xv4u4g4ODPXUK4m/zbw+Yh4B2UWvscCBwI3Zua7M/NbEXEaZQbHXSlX5HYEXpqZx2bmAxHxVuBv6liQ1cCvUq6onpKZP2kaY0S8jjI+4zOUQnpfyknr7JqHeyLiu5QxGd8AfgIMOpF1nA+8PiKuqPl7DWW8Urd3U34w/HNEvI0yhuZnge0z8y+BGygn2xMi4g7g3nplvsn3d3VEXACcWmc73ECZHfLuhl/L1hHRmdVxB+AZlO92O8rkIfcP3LKKiF+gTGbx0RrjTpQfVldmZuMrm5n55jlWuZgy8cUZEfFhygQgb6Lrim5EHA78FmVWxRuAPSjdqi6uy4ceAwPi+m5EbKB8N5d0pX81It4EvCsi9qv53wA8hnIF/XGUMUpk5u31OH5vRDyecpw/jDKe6XmZeVTPx54LnESZ6bHv91LHNP0M8MeDYpe0uFk2WzZj2WzZvNjlFMyS5WO2HgyYKZBysvs28NmutMcA76UUIvdQrrZ+Anh2z3anUE7GnXXO6HnvkyhXM++p6/1ew5genJGRciJdRTmR/oRyK4a3s/EslYdQCtGf1G2X8dCMj7/R5/0fRbmlwq2Ubkt/1C8W4PGUE+5tlMLuSuDYruWvqN/dPeXfcl7f306UcTV3UbrYvJlSwK1tsB+zPh6gFEyXA28DHjeP42FXyi0Erq/f243UsS1zbLcWeMeQ5ZvM+Ej5YXId5QfIpZTbDjz4PsBTKLcv+D6lO9k6ShesnZseAwNieT+wZsCy51ImIbmFMtbnZsqYo2OB6Fn3NykTm/y4HgtfAv5Hn/eMGlsCTxrwuUdRZrvcftLnBB8+fEz+0a/sqemWzZbNls2WzYv2EXXHSJIGqBN6XAbsmZk3TjoegIg4lzKBzGsmHYskSeNm2awmrOxKUgMRsQq4IjP/aApi2YsyXukXMvPaSccjSdIkWDZrLk5QJWmoiHhY9/3neh+Tjm+M3kjpDjUN9gR+28JUkhYny+YHWTZrKFt2JQ0VEWcCJwxZ5QmZuXY80UiSJMtmqRkru5KGiohlPHTj9X6uysx7hiyXJEkLyLJZasbKriRJkiSpdRyzK0mSJElqHSu7kiRJkqTWsbIrSZIkSWodK7uSJEmSpNaxsitJkiRJap3/D37jPZPkwGPJAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1152x576 with 2 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"plt.suptitle('S/B selection from XGBoost', fontsize='15')\n",
"Ds_mass_MC =[MC_Ds_sig_dict[\"Ds_ConsD_M\"][i][0]/1000 for i in range(m_s)]\n",
"Dplus_mass_MC =[MC_Dplus_sig_dict[\"Dplus_ConsD_M\"][i][0]/1000 for i in range(m_plus)]\n",
"XG_selected = X_dict[i][np.argmax(output_XG,1).astype(np.bool)]\n",
"Ds_mass_sel_XG = XG_selected[:,dim-1]/1000\n",
"Ds_mass_train_XG =X_dict[i][:,dim-1:dim]/1000\n",
"\n",
"plt.subplot(1,2,1)\n",
"plt.hist(Ds_mass_MC+Dplus_mass_MC,bins=70, label='MC signal events');\n",
"plt.legend(fontsize='15')\n",
"plt.ylabel('# events', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"plt.subplot(1,2,2)\n",
"\n",
"plt.hist(Ds_mass_sel_XG,alpha=0.6,bins=70, label='S selected from test set');\n",
"plt.hist(Ds_mass_train_XG,alpha=0.2,bins=70, label='Test set (S+B)');\n",
"plt.legend(fontsize='15')\n",
"plt.ylabel('# events', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"fig=plt.gcf();\n",
"fig.set_size_inches(16,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+PATH+'/D_s_XG.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"task='TEST'"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"scrolled": true
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAHoCAYAAACW1sg1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3XuYXWV9N/zvjwCCohUhKC1q1Ad9oBiiDQEF5aRAKyoH5SCiqC2oD0XxGNq3FQ99QTzQR18tDwXxhBClangFq4jgEcWgEQMpohghiIqoFOUYcj9/7J04TCaZPclkZlb4fK5rX3vvte691m/tWQz5zn2ve1VrLQAAANBFG012AQAAALC2hFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgszae7ALW1tZbb91mzJgx2WUAAACwHlx11VW/aa1NH61dZ0PtjBkzsmDBgskuAwAAgPWgqn4+SDvDjwEAAOgsoRYAAIDOEmoBAADorM5eUwsAADDcfffdl6VLl+buu++e7FIY0GabbZbtttsum2yyyVp9XqgFAAA2GEuXLs3DH/7wzJgxI1U12eUwitZabrvttixdujRPeMIT1mobhh8DAAAbjLvvvjtbbbWVQNsRVZWtttpqnXrWhVoAAGCDItB2y7r+vIRaAAAAOss1tQAAwAbr9Et+PK7bO/G5T17j+ttuuy377rtvkuSXv/xlpk2blunTpydJrrzyymy66abjWs+6esUrXpG5c+dm++23z2mnnZa5c+cmSe6///7stdde+cY3vjHJFY5OqAUAABgnW221VRYuXJgkOfnkk7PFFlvkTW960wPatNbSWstGG03+wNlzzjknSbJs2bKceuqpK0PttGnTOhFoE8OPAQAA1ruf/OQn2WmnnfLqV786T3/603PTTTflkY985Mr1559/fv72b/82SfKrX/0qhxxySGbPnp05c+bkO9/5zirbO+uss3LwwQdn//33z1Oe8pS8613vWrnutNNOy0477ZSddtopH/zgB5Mkd9xxR/76r/86O++8c3baaadccMEFSZI99tgjCxcuzNy5c3PHHXdk1qxZednLXpZly5atrO/QQw/Nl7/85ZXbf+lLX5r58+dn2bJlecMb3pA5c+Zk5syZOeuss5IkN998c/bYY4/MmjUrO+20U7797W+P87f5QHpqAQAAJsC1116bc845J2eccUaWLVu22nYnnHBC3vKWt2S33XbLkiVLcuCBB2bRokWrtLvyyiuzaNGibLrpptlll11y4IEH5t577825556bK6+8Mvfff3/mzJmTPffcM4sXL86MGTPyxS9+MUly++23P2Bbp556as4666yVvcxD6zviiCMyb9687Lfffrn77rvzta99LWeffXbOPPPMbLPNNrnyyitzzz33ZLfddst+++2X8847L89//vPz1re+Nffff3/uuuuu8fj6VkuoBQAAmABPetKTsssuu4za7itf+Uquu+66le9/97vf5a677srmm2/+gHb7779/ttxyyyTJQQcdlG9+85u55557cuihh+ahD33oA5bvvffemTt3bubOnZvnP//52X333Qeu+3nPe17e+MY35r777stFF12UffbZJw95yEPy5S9/OYsXL87555+fpBeUr7/++uyyyy457rjjcvfdd+eggw7KzjvvPPC+1oZQCwAAMAEe9rCHrXy90UYbpbW28v3Q+7S21gaaVGr4rXCq6gHbHGqHHXbIggULcvHFF+fNb35zDjzwwPzDP/zDQHU/9KEPze67755LLrkk8+bNyyte8YqVdX74wx9eOTHWUJdffnkuuuiiHHXUUTnppJNy1FFHDbSvteGaWgAAgAm20UYbZcstt8z111+f5cuX53Of+9zKdc95znPyoQ99aOX7FUOCh/vyl7+c3//+97nzzjszf/787L777nn2s5+dz33uc7nrrrvyhz/8IfPnz8+znvWs3Hzzzdliiy1y9NFH5w1veEO+//3vP2BbG2/c6+9c3bDoI444ImeffXauuOKKPOc5z0nS6yn+8Ic/vPIz1113Xe666678/Oc/z2Me85gce+yxOeaYY/KDH/xg7b+oAeipBQAANlij3YJnMr373e/OAQcckMc97nHZcccdc8899yRJPvShD+U1r3lNzjnnnCxbtix77733A0LuCnvssUde8pKX5Kc//WmOPvrozJo1K0ly5JFHrhzm/JrXvCZPfepTc/HFF2fu3LnZaKONsummm+aMM85YZXuvetWrMnPmzMyePTsf+chHHrDugAMOyMtf/vK8+MUvziabbJIkOe6443LjjTeu3O8222yT+fPn59JLL8373//+bLLJJtliiy3yyU9+cvy+tBHU6rqnp7rZs2e3BQsWTHYZAADAFLJ48eLssMMOk13GenfWWWdl0aJF+dd//dfJLmVcjPRzq6qrWmuzR/us4ccAAAB0luHHAAAAHbPinrYItevXZaeM3mbvk9Z/HQAAABsow48BAADoLKEWAACAzhJqAQAA6CzX1AIAABuuQea5GYsB5sSZNm1anvrUp658//nPfz4zZswYse2SJUty4IEHZtGiReNV4VpbsGBBPv7xj+cDH/hALr/88my66aZ55jOfmSQ544wz8tCHPjQve9nLJrnKVQm1AAAA42jzzTfPwoULJ7uMMZs9e3Zmz+7dFvbyyy/PFltssTLUvvrVr57M0tbI8GMAAID1bMmSJXnWs56Vpz/96Xn605+eb3/726u0ueaaazJnzpzMmjUrM2fOzPXXX58k+eQnP7ly+XHHHZf7779/lc/OmDEjb33rWzNnzpzMmTMnP/nJT5IkP//5z7Pvvvtm5syZ2XfffXPjjTcmST7zmc9kp512ys4775xnP/vZSXpB9sADD8ySJUtyxhln5PTTT8+sWbPyjW98IyeffHLe+973ZvHixZkzZ84DjmvmzJlJkquuuip77rln/uqv/ir7779/brnlliTJBz7wgey4446ZOXNmjjjiiHH8VnuEWgAAgHF01113ZdasWZk1a1YOPvjgJMk222yTSy65JN///vczb968nHDCCat87owzzsjrXve6LFy4MAsWLMh2222XxYsXZ968efnWt76VhQsXZtq0aTn33HNH3O8jHvGIXHnllTn++OPz+te/Pkly/PHH52Uve1muvvrqHHXUUSv3+453vCNf+tKX8sMf/jAXXnjhA7YzY8aMvPrVr86JJ56YhQsX5lnPetbKdTvssEPuvffe3HDDDUmSefPm5bDDDst9992Xv//7v88FF1yQq666Kq985Svzj//4j0mSU089NT/4wQ9y9dVX54wzzljHb3dVhh8DAACMo5GGH9933305/vjjVwbTH//4x6t87hnPeEb+5V/+JUuXLs0hhxyS7bffPpdeemmuuuqq7LLLLkl6gXmbbbYZcb9HHnnkyucTTzwxSXLFFVfks5/9bJLk6KOPzlve8pYkye67755jjjkmhx12WA455JAxHd9hhx2WT3/605k7d27mzZuXefPm5brrrsuiRYvy3Oc+N0ly//33Z9ttt02SzJw5M0cddVQOOuigHHTQQWPa1yCEWgAAgPXs9NNPz6Mf/ej88Ic/zPLly7PZZput0uYlL3lJdt1111x00UXZf//9c9ZZZ6W1lpe//OU55ZTRJ7yqqhFfj9TmjDPOyHe/+91cdNFFmTVr1piuAT788MPz4he/OIccckiqKttvv31+9KMf5S//8i9zxRVXrNL+oosuyte//vVceOGFeec735lrrrkmG288flHU8GMAAID17Pbbb8+2226bjTbaKJ/4xCdGvC72hhtuyBOf+MSccMIJecELXpCrr746++67by644IL8+te/TpL89re/zc9//vMR9zFv3ryVz894xjOSJM985jNz/vnnJ0nOPffc7LHHHkmSn/70p9l1113zjne8I1tvvXVuuummB2zr4Q9/eO64444R9/OkJz0p06ZNyzvf+c4cfvjhSZKnPOUpufXWW1eG2vvuuy/XXHNNli9fnptuuil77713TjvttPz+97/PH/7whzF9d6PRUwsAAGy4BrgFz0R47Wtfm0MPPTSf+cxnsvfee+dhD3vYKm3mzZuXT37yk9lkk03ymMc8Jv/8z/+cRz3qUXnXu96V/fbbL8uXL88mm2ySD33oQ3n84x+/yufvueee7Lrrrlm+fHnOO++8JL1Jml75ylfmPe95T6ZPn55zzjknSfLmN785119/fVpr2XfffbPzzjvna1/72sptPf/5z8+LXvSizJ8/Px/84AdX2dfhhx+eN7/5zfnZz36WJNl0001zwQUX5IQTTsjtt9+eZcuW5fWvf32e/OQn56UvfWluv/32tNZy4okn5pGPfOS4fKcrVGttXDc4UWbPnt0WLFgw2WWs2SD3xJoi/5EBAMCGYPHixdlhhx0mu4wJN2PGjCxYsCBbb731ZJeyVkb6uVXVVa212aN91vBjAAAAOsvwYwAAgI5bsmTJZJcwafTUAgAAG5SuXmL5YLWuPy+hFgAA2GBsttlmue222wTbjmit5bbbbhvxFkeDMvwYAADYYGy33XZZunRpbr311skuhQFtttlm2W677db68xMaaqtqsyRfT/KQ/r4vaK29raqekOT8JI9K8v0kR7fW7p3I2gAAgO7bZJNN8oQnPGGyy2ACTfTw43uS7NNa2znJrCQHVNVuSd6d5PTW2vZJfpfkVRNcFwAAAB00oaG29fyh/3aT/qMl2SfJBf3lH0ty0ETWBQAAQDdN+ERRVTWtqhYm+XWSS5L8NMnvW2vL+k2WJvmL1Xz22KpaUFULjJEHAABgwkNta+3+1tqsJNslmZNkh5GareazZ7bWZrfWZk+fPn19lgkAAEAHTNotfVprv09yeZLdkjyyqlZMWrVdkl9MVl0AAAB0x4SG2qqaXlWP7L/ePMlzkixOclmSF/WbvTzJ/ImsCwAAgG6a6PvUbpvkY1U1Lb1A/enW2heq6tok51fVu5L8IMnZE1wXAAAAHTShoba1dnWSp42w/Ib0rq8FAACAgU3aNbUAAACwroRaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOmtCQ21VPbaqLquqxVV1TVW9rr/85Kq6uaoW9h9/M5F1AQAA0E0bT/D+liV5Y2vt+1X18CRXVdUl/XWnt9beO8H1AAAA0GETGmpba7ckuaX/+o6qWpzkLyayBgAAADYck3ZNbVXNSPK0JN/tLzq+qq6uqo9U1ZaTVRcAAADdMSmhtqq2SPIfSV7fWvvvJP+W5ElJZqXXk/u+1Xzu2KpaUFULbr311gmrFwAAgKlpwkNtVW2SXqA9t7X22SRprf2qtXZ/a215kn9PMmekz7bWzmytzW6tzZ4+ffrEFQ0AAMCUNNGzH1eSs5Msbq29f8jybYc0OzjJoomsCwAAgG6a6NmPd09ydJIfVdXC/rJ/SHJkVc1K0pIsSXLcBNcFAABAB0307MffTFIjrLp4IusAAABgwzBpsx8DAADAuhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOEmoBAADoLKEWAACAzhJqAQAA6CyhFgAAgM6a0FBbVY+tqsuqanFVXVNVr+svf1RVXVJV1/eft5zIugAAAOimie6pXZbkja21HZLsluR/VdWOSeYmubS1tn2SS/vvAQAAYI02HusHqupJSR6TZLMkv03yk9baHYN8trV2S5Jb+q/vqKrFSf4iyQuT7NVv9rEklyd561hrAwAA4MFl1FBbVRsl+eskL0uyb5Itk1R/dUuyvKquSXJBko+31m4cZMdVNSPJ05J8N8mj+4E3rbVbqmqb1Xzm2CTHJsnjHve4QXYDAADABmyNw4+r6qgk1yU5L70A+44k+yR5apInJ9k1yZFJ/jPJi5L8pKrOrqrtRtnuFkn+I8nrW2v/PWixrbUzW2uzW2uzp0+fPujHAAAA2ECN1lP71iRvT/KZ1to9q2mzIL1e2rlV9ZQkr08v6L5npMZVtUl6gfbc1tpn+4t/VVXb9ntpt03y6zEeBwAAAA9Cawy1rbWZY9lYa+26JK9Z3fqqqiRnJ1ncWnv/kFUXJnl5klP7z/PHsl8AAAAenMY8UdQ62j3J0Ul+VFUL+8v+Ib0w++mqelWSG5O8eILrAgAAoIPGJdRW1WOT1GiTRLXWvpk/TTI13L7jUQsAAAAPHuPVU3tDepNOTRun7QEAAMCoxivUviqr74EFAACA9WJcQm1r7ePjsR0AAAAYizXepxYAAACmsoF7aqvq06O1aa0dtm7lAAAAwODGMvx4+gjLHpXkKUluS3LduFQEAAAAAxo41LbW9h5pef92Pp9Lcvp4FQUAAACDWOdraltrNyU5Jclp614OAAAADG68Joq6P8l247QtAAAAGMhYJoracYTFmybZIck7k3xvvIoCAACAQYxloqhFSdoIyyu9QPu341IRAAAADGgsoXakiaLuTrK0tXbzONUDAAAAAxvL7MdfW5+FAAAAwFiN10RRAAAAMOHGJdRW1fVV9dPx2BYAAAAMaizX1K7J16PXFwAAgAk2LqG2tfaq8dgOAAAAjIXeVQAAADprTD21VVVJdk/y5CSbDV/fWvvwONUFAAAAoxo41FbVo5NcmmTHJC1J9Ve1Ic2EWgAAACbMWIYfvy/J7Ukem16g3TXJjCT/lOT69HpvAQAAYMKMZfjxnklel+SW/vtqrd2Y5P+tqo3S66Xdf5zrAwAAgNUaS0/tI5Pc2lpbnuS/k2wzZN23kzxzPAsDAACA0Ywl1P4sybb919ckOWrIuucn+e14FQUAAACDGMvw44uS7Jfk00nelWR+VS1Ncl+SxyV56/iXBwAAAKs3cKhtrZ005PUXq+qZSQ5OsnmSS1prX1wP9QEAAMBqjek+tUO11hYkWTCOtQAAAMCYrPGa2qqaMdYNVtVGVfXYtS0IAAAABjXaRFE/qqpzq2qfqqo1Nayq7arqjUl+kuSIcasQAAAAVmO04cf/M8n/k+T/T3JnVV2ZZFGS3yS5J73b/DwhyV8l2THJ1Une0lq7YL1VDAAAAH1rDLWttZuTvKaq3pLk8CT7JDkoyWOSbJbebXyuSy/0vqK1dtX6LRcAAAD+ZKCJolprdyQ5q/8AAACAKWG0a2oBAABgyhJqAQAA6CyhFgAAgM4SagEAAOgsoRYAAIDOWqdQW1VbVtWsqnrIeBUEAAAAgxo41FbV26vq1CHv90lyY5Krkvy0qv5yPdQHAAAAqzWWntqjkvzXkPfvS/LNJLsnuS7JKeNYFwAAAIxqLKH2z5PckCRV9dgkOyd5W2vtO0nen2S38S8PAAAAVm8sofaOJH/Wf71Pkt+11q7sv787yUPHszAAAAAYzcZjaPu1JHOranmSNyWZP2Tdk5PcNJ6FAQAAwGjG0lN7YpJ7kpyf5PdJ/nHIupcl+fo41gUAAACjGrintrV2c3rDjkeyf5K7xqUiAAAAGNBYbunz1ar6n6tZ/ZgkXxqfkgAAAGAwYxl+vFeSR6xm3SOSPHudqwEAAIAxGEuoTZI2fEFVbZresORfjktFAAAAMKA1XlNbVW9L8s/9ty3Jd6pqdc3fM451AQAAwKhGmyjq4iS/SVJJPpDkfUmWDGtzb5L/aq19Y9yrAwAAgDVYY6htrX0vyfeSpKruSHJRa+03E1EYAAAAjGYst/T52PosBAAAAMZq4FBbVZskeV2SQ5Jsl2Sz4W1aa9uMX2kAAACwZgOH2iSnJzkuyReSXJbetbSsq8tOGazd3iet3zoAAAA6aCyh9sVJ5rbW3re+igEAAICxGMt9aivJ1eurEAAAABirsYTaf09y5LrsrKo+UlW/rqpFQ5adXFU3V9XC/uNv1mUfAAAAPHiMZfjxr5IcVVWXJbkkye+HrW+ttX8bZRsfTfL/Jfn4sOWnt9beO4ZaAAAAYEyh9l/7z49LsucI61uSNYba1trXq2rGGPYJAAAAqzXw8OPW2kajPKatQx3HV9XV/eHJW67DdgAAAHgQGcs1tevLvyV5UpJZSW5JstrZlavq2KpaUFULbr311omqDwAAgClqTKG2qrapqndX1aVV9eOq+sv+8tdV1TPWpoDW2q9aa/e31panNxnVnDW0PbO1Nru1Nnv69OlrszsAAAA2IAOH2qqak+T6JIcmWZJe7+pD+qu3TfLGtSmgqrYd8vbgJItW1xYAAACGGstEUacnuSzJIemF4VcMWXdlkpeMtoGqOi/JXkm2rqqlSd6WZK+qmpXeRFNLkhw3hpoAAAB4EBtLqH16khe21pZXVQ1bd1uSbUbbQGttpPvcnj2GGgAAAGClsVxTe3uS1V3I+sT07mMLAAAAE2YsoXZ+krdX1ROHLGtVtXWSNyX57LhWBgAAAKMYS6idm+S/k1yb5Ov9ZWckuS7JXUn+eXxLAwAAgDUb+Jra1trvqmq3JEcn2TfJH5P8NslZST7eWrtn/ZQIAAAAIxvLRFFprd2b3sROJncCAABg0o3lPrVfq6rXVNXqJosCAACACTWWa2pvTfLeJL+oqkuq6pVVteV6qgsAAABGNXCoba29KL170b4syR+SfCjJL6vqC1V1dFU9fD3VCAAAACMaS09tWmt/bK2d11o7OL2A+3f9Vf+e5JfjXRwAAACsyZhC7VCttTuS/DTJz9K71c/m41UUAAAADGLMobaq5lTV+6rqxvTuV7tnkv+dZPvxLg4AAADWZOBb+lTVqUkOS/L4JNcnOSfJ+a21xeupNgAAAFijsdyn9rAkn04vyC5cT/UAAADAwAYOta21J67PQgAAAGCs1nhNbVW9pKoeNWzZ46pq42HL/ryq/mF9FAgAAACrM9pEUZ9I8j9WvKmqaenNdjxzWLvHJnnn+JYGAAAAazZaqK0BlwEAAMCEW+v71AIAAMBkE2oBAADorEFCbRtwGQAAAEyoQW7p86WqWjZs2aXDlo3lfrcAAAAwLkYLo2+fkCoAAABgLawx1LbWhFoAAACmLBNFAQAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdNaGhtqo+UlW/rqpFQ5Y9qqouqarr+89bTmRNAAAAdNdE99R+NMkBw5bNTXJpa237JJf23wMAAMCoJjTUtta+nuS3wxa/MMnH+q8/luSgiawJAACA7poK19Q+urV2S5L0n7dZXcOqOraqFlTVgltvvXXCCgQAAGBqmgqhdmCttTNba7Nba7OnT58+2eUAAAAwyaZCqP1VVW2bJP3nX09yPQAAAHTEVAi1FyZ5ef/1y5PMn8RaAAAA6JCJvqXPeUmuSPKUqlpaVa9KcmqS51bV9Ume238PAAAAo9p4InfWWjtyNav2ncg6AAAA2DBMheHHAAAAsFaEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6a+PJLoABXXbK6G32Pmn91wEAADCF6KkFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLOEWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOmvjyS6AcXTZKYO12/uk9VsHAADABJkyobaqliS5I8n9SZa11mZPbkUAAABMdVMm1Pbt3Vr7zWQXAQAAQDe4phYAAIDOmkqhtiX5clVdVVXHjtSgqo6tqgVVteDWW2+d4PIAAACYaqZSqN2zR20lAAARNklEQVS9tfb0JH+d5H9V1bOHN2itndlam91amz19+vSJrxAAAIApZcqE2tbaL/rPv07yuSRzJrciAAAApropEWqr6mFV9fAVr5Psl2TR5FYFAADAVDdVZj9+dJLPVVXSq+lTrbX/nNySAAAAmOqmRKhtrd2QZOfJrgMAAIBumRLDjwEAAGBtTImeWibYZacM1m7vk9ZvHQAAAOtITy0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWUItAAAAnSXUAgAA0FkbT3YBTGGXnTJ6m71PWv91AAAArIaeWgAAADpLqAUAAKCzhFoAAAA6S6gFAACgs4RaAAAAOkuoBQAAoLPc0ocHvdMv+fFA7U587pPXcyUAAMBY6akFAACgs/TUsm4uO2WwdnuftH7rAKBzjJQBYDzoqQUAAKCz9NQCAEwxerEBBqenFgAAgM7SUzvBrrjhtoHaPeOJW63nSkY2aH2DWnkcg1x767pbNjAPtp6WB9vxwlTgvzsAPbUAAAB0mJ5aOmUy/yLtr+HrxvfHChvSubAhHQsbtgfjufpgO+YH2/HCUHpqAQAA6Cw9tR033tfATpYrbrgt31k22F8Yx9Vlp2S3G0f/Dr/zuGMnoBhgqC6MzIAHK72C3eVnx4ZITy0AAACdpaeWKWO3G88cqJ1e05H5yyvA4Ma7N97vVoDJo6cWAACAztJTS+cM1KN72Vbjet/bQfa5ofUgb0jXFOrFZoUN6byeDP5bWj3n1si6cM50oUZgzYRa1qtJncjqslMmb99T2FT/h9eD8R8XU30Y5FQ/Z7pgqp/X4/0znurbg7Ux1f87Hm8PtuOl2ww/BgAAoLP01LJBmswe4vH+y6YeCmBt+f0BE08PJ0w8PbUAAAB0lp5aGCddvyXRZNS/LvtcHz1QG0qv1oZyHKyenzHry1Q/t6Z6fayeHuzV892sO6EWJtjKIHfZVmtuOI6zNwMAwIZKqIVJMtp1v99ZNvpf7Qbt6QSAqezB2AO7ofTObUiznU/175rVE2oBpoCuD18HAJgsQu0UNan3d2VK0As7NkLhAw3yfWwo38WG1MOzIR0LMD78Xpg4U/2+8ayeUAuMyVSeUGo8TdVQOFW/i2TDCckAQLcItcB6oaf5Tx4s38V4HudUmWUbAKYKPcmrN2VCbVUdkOR/J5mW5KzW2qmTXBLAGj1Ywirrx3ifP4OEcgF//ZmqozsAHgymRKitqmlJPpTkuUmWJvleVV3YWrt2cisDNjSC6J9M5e9iKg+zHsSDIbxM5fNnEIP+jKbqcU7lS0Em4w8s/qgw9T0Y/6jmeuiJU621ya4hVfWMJCe31vbvvz8pSVprp6zuM7Nnz24LFiyYoArX0mWrlm8CKAAAxmq8w95U/YPNeJqqf9SZKrow/LiqrmqtzR613RQJtS9KckBr7W/7749Osmtr7fhh7Y5NsuJMeUqS6ya00LHZOslvJrsIWAvOXbrIeUtXOXfpKucuE+HxrbXpozWaEsOPk9QIy1ZJ2621M5N04s9KVbVgkL8qwFTj3KWLnLd0lXOXrnLuMpVsNNkF9C1N8tgh77dL8otJqgUAAICOmCqh9ntJtq+qJ1TVpkmOSHLhJNcEAADAFDclhh+31pZV1fFJvpTeLX0+0lq7ZpLLWledGCYNI3Du0kXOW7rKuUtXOXeZMqbERFEAAACwNqbK8GMAAAAYM6EWAACAzhJq11FVHVBV11XVT6pq7gjrH1JV8/rrv1tVMya+SnigAc7bN1TVtVV1dVVdWlWPn4w6YbjRzt0h7V5UVa2q3G6CKWGQc7eqDuv/7r2mqj410TXCcAP8e+FxVXVZVf2g/2+Gv5mMOsE1teugqqYl+XGS56Z3W6LvJTmytXbtkDavTTKztfbqqjoiycGttcMnpWDIwOft3km+21q7s6pek2Qv5y2TbZBzt9/u4UkuSrJpkuNbawsmulYYasDfu9sn+XSSfVprv6uqbVprv56UgiEDn7dnJvlBa+3fqmrHJBe31mZMRr08uOmpXTdzkvyktXZDa+3eJOcneeGwNi9M8rH+6wuS7FtVNYE1wnCjnrettctaa3f2334nvXtHw2Qb5HdukrwzyWlJ7p7I4mANBjl3/y7Jh1prv0sSgZYpYJDztiV5RP/1nyX5xQTWBysJtevmL5LcNOT90v6yEdu01pYluT3JVhNSHYxskPN2qFcl+eJ6rQgGM+q5W1VPS/LY1toXJrIwGMUgv3efnOTJVfWtqvpOVR0wYdXByAY5b09O8tKqWprk4iR/PzGlwQNNifvUdthIPa7Dx3MP0gYm0sDnZFW9NMnsJHuu14pgMGs8d6tqoySnJzlmogqCAQ3ye3fjJNsn2Su90THfqKqdWmu/X8+1weoMct4emeSjrbX3VdUzknyif94uX//lwZ/oqV03S5M8dsj77bLqsIuVbapq4/SGZvx2QqqDkQ1y3qaqnpPkH5O8oLV2zwTVBmsy2rn78CQ7Jbm8qpYk2S3JhSaLYgoY9N8L81tr97XWfpbkuvRCLkyWQc7bV6V3LXhaa1ck2SzJ1hNSHQwh1K6b7yXZvqqeUFWbJjkiyYXD2lyY5OX91y9K8tVmdi4m16jnbX8I5/9JL9C6roupYo3nbmvt9tba1q21Gf2JSr6T3jlsoigm2yD/Xvh8kr2TpKq2Tm848g0TWiU80CDn7Y1J9k2SqtohvVB764RWCRFq10n/Gtnjk3wpyeIkn26tXVNV76iqF/SbnZ1kq6r6SZI3JFntLShgIgx43r4nyRZJPlNVC6tq+P/EYMINeO7ClDPgufulJLdV1bVJLkvy5tbabZNTMQx83r4xyd9V1Q+TnJfkGJ03TAa39AEAAKCz9NQCAADQWUItAAAAnSXUAgAA0FlCLQAAAJ0l1AIAANBZQi0Ak6KqvlpVP6yqjYctP7SqWlU9d9jynavqU1V1c1XdW1W/raqvVNXhQ7dRVR/tf37F446q+l5VHTJRxzas7k2r6uSqmjVA272G1X5fVd1QVadV1cOGtV0ypN09VfWLqrq4qo6uqo1W0251j2PWw6EPP7Zjq+qgru8DgKln49GbAMB68dokVyc5Icn7k6Sqtkjyr+ndD/GSFQ2r6kVJPpXkW0lOSrIkyaOS/E2Sj6d3X+Wzh2z7v5K8ov/6EUmOSe++y3u21r653o5oZJsmeVt6NS8c8DNHJbkhySZJnp7kX5I8Msmxw9p9KskHk0xLsm2S/dP7Ho6qqhe01u5NcnCShwz5zH8muSDJWUOW/XTww1lrxyZZlOTzHd8HAFOMUAvApGit/VdVvS/J26tqXmvt5iRvT/JnSU5c0a6q/iLJR5Ocm+SV7YE3WP98fxvbDtv8H1tr3xmyja8k2TvJC5JMdKhdG1e31hb1X3+j/x0ck1VD7S1DjzPJBVX16SRfTC/8v7219oOhH6iqZUmWDvscAHSW4ccATKZ3JrktyelVNTO9Xtu3tdZ+MaTN36b3R9g3Dgu0SZLW2nWttcvXtJPW2vIkd6bX87lSVc2qqkur6s6q+l1VnVtVjx7WZuuq+lhV3dZvd3lVzR7W5gVVdVVV/bG/ne9W1Z791Xf0n88ZMtx3xpq/llXcMbz2NRzrJen1xL5mjPsYUVXt0z+eu6vqV1X14X6P+or1x/SPaYthn1tSVe/tv748yV8lefnwIc8r2lXVP1XVL6vqD/2fw5+N1z4A2LAJtQBMmtbanUlel+TFSeYnuTa94bRDPTvJgtbab8ey7arauP94VFW9KcmM/j5WrJ+e5PIkD03ykiR/n2TPJJdU1aZDNvX59Ib1vinJ4en9v/Oyqvof/e08Kb0Q+dUkz09v6PAX0hsenST79J/fleQZ/ccto5Q/rV/75lW1e5LjknxuDId/SZJHr0V4foCq2jG94cq/SXJoesOoX5Le8Y7Fa9MbEn5x/vQdXDRk/ZFJnpPk75K8Icnz8sDh0eOxDwA2UIYfAzCpWmvzq+qq9HrZ9mmtLRvW5M8zwrWo9cAJppb3e2NX+Ksk9w1dn+Qtw3p039h/3r+19t/9bf44yXfTC3DnVdUBSXZPsldr7Wv9Nl9N7/rYN6cXNp+W5I7W2puHbPviIa+/13/+6RiG/A4/3m+mF/4HtbT//Oh+rWvrn5P8PMkLWmv3J0lV/TbJvKp6RmvtikE20lq7tqr+mOTW1XwHmyd5XmvtD/19/DHJJ6pqh9ba4nHaBwAbKD21AEyq/lDepyVpSfYaqUl/3fDP3Dfk8elhn1mcZJf+Y8/0wtm/DBuOOifJl1cE2iRprV2ZXgjcY0ibW1cE2n6bP6bXE7uizY+S/Fl/iPJ+w2cpXktH9GvfNb1ezK3Tu3540P9v1zjUkPSO/3MrAm3ffyRZlj8d/3i4ZEWg7ftsesewyzjuA4ANlJ5aACZNP6T9W5Ir0hsyO7eqPtZau2FIs5uTbDfso9fmT4Hn/4yw6TtbawuGvP96VT0myWn97bf0Jpe6ZoTP/ip/Gjq8bf/9atu01q6rqhcmmZteD+19VfW5JK9rrd060nEP4JohE0VdWVXXJ1mQ3mzPXxjg838xpM51scrxt9bur6rb8qfvaDz8etg+7qqqP2TVCcAAYBV6agGYTK9Or5f2tUlOTS/AfmBYm68n2aWqtlyxoLV2Z2ttQT+43pHBXJtkenq9nknvutZtRmj36CS/HUObtNYuaq09K8lWSV6V3vWhw68NXhfX9p93GLD9fkl+2Vpbso77XeX4q2paese54vjv7j8PvQ45SbbM4IbvY/P0btO04trj8dgHABsooRaASVFV26R3/9UPttaubq3dk97sx8/r93yucFaS+5O8Zx13uVOSu9KbbTnpXTu7f1U9fEhNu6Q3odQ3h7TZpqqePaTNQ9ObyGiVWwO11m5vrX0qvUmdduwvvrf/vNk61p4kN43WsKqem+RF6fWAr6vvJjm4H2RXOCS9kV4rjn/F9bsrA3dV7Zre/YGHujer/w6eO2xm40PSG3K+ord9PPYBwAbK8GMAJst70wuZb1uxoLV2cVXNT/KvVfXl1tpdrbWbq+oVSc6tqicmOSe96163SDI7ycwkFw7b9sOqarf+682TPCu9mXU/PGRCqfend9ubL1XVu/vbOzW9a2T/o1/Pl6rqW+lNjDQ3vUD8pv4235MkVXVcejPt/meSXyTZPr3ZnD/e38a9VfWzJIdV1aL0eh2vbq2tCLsjmdkPeRsleWKSf0pyY3r3nx1q2/5xTkvymPRmaT4mvaHcp6xh+4N6V5IfpHc977+lNwz83Um+NGSSqCvT72Gvqn9Kb1jyW5L897Bt/Vd6f0TYP73v8WettRV/YLgryUVV9Z70hhy/J71reVf0UI/HPgDYULXWPDw8PDw8JvSR3m16WpIjR1j3+CR/TPLOYctnJTkvveB4X3rDX7+a3gzEmwxp99H+tlc87kpv+O7cJJsO2+bT+tu4M8nvk3wqyaOHtZmeXkD9XX9bX0uyy5D1K24d84v0AuvP0gt+DxnSZr8kV/fXtyQzVvO97DWs9uXp9VKel+RJw9ouGdLu3vSG6n4xydFJNlrDd/+bJCeP4We1b3o9tnend+3rh5NsMazNLunN8nxneiF493597x3S5olJvpLk9n7Nxww5jvclOTm963f/2D/eR47XPjw8PDw8NuxHtbbKfewBACZEVS1JckFr7U2TXQsA3eSaWgAAADpLqAUAAKCzDD8GAACgs/TUAgAA0FlCLQAAAJ0l1AIAANBZQi0AAACdJdQCAADQWf8XWiSnBd2Wc98AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1152x576 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"true_positives_XG=output_XG[:,1][np.where(Y_test[:,1]==1)]\n",
"false_positives_XG=output_XG[:,1][np.where(Y_test[:,0]==1)]\n",
"plt.hist(true_positives_XG,alpha=0.5,bins=80,density=True,label=\"True positives\");\n",
"plt.hist(false_positives_XG,alpha=0.5,bins=80,density=True, label=\"False positives\");\n",
"plt.legend()\n",
"plt.xlabel(\"XGBoost BDT output\", fontsize='15')\n",
"plt.ylabel(\"Events (a.u.)\", fontsize='15')\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(16,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+PATH+'/tp_vs_fp_XG.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'true_positives_NN' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-23-690618031886>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0msig_eps_vals_NN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msig_eff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrue_positives_NN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0mbkg_eps_vals_NN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbkg_eff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfalse_positives_NN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m<ipython-input-23-690618031886>\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 6\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0msig_eps_vals_NN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msig_eff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtrue_positives_NN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 8\u001b[0m \u001b[0mbkg_eps_vals_NN\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mbkg_eff\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfalse_positives_NN\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mthreshold_range\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 9\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'true_positives_NN' is not defined"
]
}
],
"source": [
"threshold_range=np.linspace(0.0,1.,num=30)\n",
"sig_eps_vals_XG=[sig_eff(true_positives_XG,threshold_range[i]) for i in range(len(threshold_range))]\n",
"bkg_eps_vals_XG=[bkg_eff(false_positives_XG,threshold_range[i]) for i in range(len(threshold_range))]\n",
"if task=='TEST': \n",
"\n",
" \n",
" sig_eps_vals_NN=[sig_eff(true_positives_NN,threshold_range[i]) for i in range(len(threshold_range))]\n",
" bkg_eps_vals_NN=[bkg_eff(false_positives_NN,threshold_range[i]) for i in range(len(threshold_range))]\n",
" \n",
" plt.plot(threshold_range,threshold_range, 'black', linestyle='dashed')\n",
" plt.plot(bkg_eps_vals_NN,sig_eps_vals_NN, 'r', label=\"NN ROC Curve\")\n",
" pAUC_NN=roc_auc_score(Y_test,output_NN)\n",
" print(\"pAUC from NN {0}\".format(pAUC_NN))\n",
"\n",
"plt.plot(threshold_range,threshold_range, 'black', linestyle='dashed')\n",
"plt.plot(bkg_eps_vals_XG,sig_eps_vals_XG,'b',label=\"XG Boost ROC Curve\")\n",
"plt.xlabel(\"Background selection efficiency\", fontsize='15')\n",
"plt.ylabel(\"Signal selection efficiency\", fontsize='15')\n",
"pAUC_XG=roc_auc_score(Y_test,output_XG)\n",
"plt.text(0.69,0.1,\"NN AUC {0:.4g}\\n \\n XGBoost AUC {1:.4g}\\n\".format(pAUC_NN,pAUC_XG), bbox=dict(boxstyle=\"round\", facecolor='blue', alpha=0.10), horizontalalignment='center', verticalalignment='center',fontsize='15')\n",
"plt.legend()\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(8,8)\n",
"\n",
"print(\"pAUC from XG Boost {0}\".format(pAUC_XG))\n",
"plt.savefig('/home/hep/davide/Rphipi/'+l_flv[l_index]+'/BDTs/test_'+str(test)+'/roc_comparison_'+str(i)+'.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [],
"source": [
"# save XGBOOST model to file\n",
"pickle.dump(model, open(PATH+\"/XG_\"+str(i)+\"_.pickle.dat\", \"wb\"), protocol=2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini