make_STS.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "f52780f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Defaulting to user installation because normal site-packages is not writeable\n",
      "Requirement already satisfied: pyyaml in /global/home/users/rrmastandrea/.local/lib/python3.8/site-packages (6.0.1)\n",
      "\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.2.1\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.3.1\u001b[0m\n",
      "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n"
     ]
    }
   ],
   "source": [
    "!pip install pyyaml\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib\n",
    "import numpy as np\n",
    "import pickle\n",
    "\n",
    "import os\n",
    "from helpers.datasets import *\n",
    "from helpers.plotting import *\n",
    "\n",
    "plt.style.use(\"science.mplstyle\")\n",
    "cm = matplotlib.cm.get_cmap('Greens')\n",
    "from matplotlib.backends.backend_pdf import PdfPages as pp\n",
    "\n",
    "np.random.seed(8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "121ed9da-d6af-458f-9a46-260fcb79fecf",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "with open(\"workflow.yaml\", \"r\") as file:\n",
    "    workflow = yaml.safe_load(file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b09d0c06",
   "metadata": {},
   "outputs": [],
   "source": [
    "# directories\n",
    "feta_dir = \"/global/home/users/rrmastandrea/FETA/\"\n",
    "STS_dir = \"/global/home/users/rrmastandrea/FETA/LHCO_STS_{project_id}/data/\".format(project_id = workflow[\"project_id\"])\n",
    "os.makedirs(STS_dir, exist_ok=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "71038a0c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Generating datasets...\n",
      "events_anomalydetection_herwig_features.h5\n",
      "events_anomalydetection_v2.features.h5\n",
      "events_anomalydetection_Z_XY_qqq.features.h5\n",
      "SIM shape: (999986, 7)\n",
      "DAT shape: (999881, 7)\n",
      "SIGNAL shape: (99969, 7)\n"
     ]
    }
   ],
   "source": [
    "# dataset generation parameters\n",
    "\n",
    "bands_dict = workflow[\"bands_dict\"]\n",
    "              \n",
    "print(\"Generating datasets...\")\n",
    "\n",
    "bands_to_analyze = [\"sb1\", \"sr\", \"sb2\"]\n",
    "\n",
    "df_sim = get_lhco(sim_type = \"herwig\")\n",
    "data_sim = select_lhco_features(df_sim, feature_type = 0).to_numpy()\n",
    "\n",
    "df_dat = get_lhco(sim_type = \"pythia\")\n",
    "data_dat = select_lhco_features(df_dat, feature_type = 0).to_numpy()\n",
    "\n",
    "#df_signal = get_lhco(sim_type = \"pythia\", sm = \"signal\")\n",
    "#data_signal = select_lhco_features(df_signal, feature_type = 0).to_numpy()\n",
    "df_signal = get_lhco(sim_type = \"3_prong\", sm = \"signal\")\n",
    "data_signal = select_lhco_features(df_signal, feature_type = 0).to_numpy()\n",
    "\n",
    "# add an extra column to keep track of sig, bkg information\n",
    "data_sim = np.hstack((data_sim, np.zeros((len(data_sim), 1))))\n",
    "data_dat = np.hstack((data_dat, np.zeros((len(data_dat), 1))))\n",
    "data_signal = np.hstack((data_signal, np.ones((len(data_signal), 1))))\n",
    "\n",
    "print(\"SIM shape:\", data_sim.shape)\n",
    "print(\"DAT shape:\", data_dat.shape)\n",
    "print(\"SIGNAL shape:\", data_signal.shape)\n",
    "\n",
    "\n",
    "\n",
    "n_features = 6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0f1ad115",
   "metadata": {},
   "outputs": [],
   "source": [
    "def pull_from_mass_range(dataset, mass_ranges):\n",
    "\n",
    "        \"\"\"\n",
    "        mass_ranges is a list of lists [low_bound, high_bound]\n",
    "        \"\"\"\n",
    "\n",
    "        selected_data = []\n",
    "\n",
    "        for mass_range in mass_ranges:\n",
    "            loc = np.where((dataset[:, 5] >= mass_range[0]) & (dataset[:, 5] < mass_range[1]))[0]\n",
    "            selected_data.append( dataset[loc, :] )\n",
    "        selected_data = np.concatenate(selected_data)\n",
    "        np.random.shuffle(selected_data)\n",
    "\n",
    "        return selected_data\n",
    "\n",
    "def pull_n_signal_events(n, sig_dataset):\n",
    "    \n",
    "    selected_sig_indices = np.random.choice(sig_dataset.shape[0], size=n, replace=False)\n",
    "    selected_sig_events = sig_dataset[selected_sig_indices]\n",
    "    \n",
    "    return(selected_sig_events)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7996643f",
   "metadata": {},
   "source": [
    "## Count number of events ber band"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "4470a897",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of events in band sb1:\n",
      "Herwig bkg: 728210\n",
      "Pythia bkg: 734634\n",
      "Pythia sig: 9308\n",
      "\n",
      "Number of events in band sr:\n",
      "Herwig bkg: 121978\n",
      "Pythia bkg: 121339\n",
      "Pythia sig: 79781\n",
      "\n",
      "Number of events in band sb2:\n",
      "Herwig bkg: 141038\n",
      "Pythia bkg: 135990\n",
      "Pythia sig: 10863\n",
      "\n"
     ]
    }
   ],
   "source": [
    "herwig_bkg_dict = {}\n",
    "pythia_bkg_dict = {}\n",
    "pythia_sig_dict = {}\n",
    "\n",
    "num_sig_pythia_events = {}\n",
    "num_bkg_pythia_events = {}\n",
    "\n",
    "for b in bands_to_analyze:\n",
    "    print(f\"Number of events in band {b}:\")\n",
    "    loc = pull_from_mass_range(data_sim, [bands_dict[b]])\n",
    "    herwig_bkg_dict[b] = loc\n",
    "    print(f\"Herwig bkg: {len(loc)}\")\n",
    "    \n",
    "    loc = pull_from_mass_range(data_dat, [bands_dict[b]])\n",
    "    pythia_bkg_dict[b] = loc\n",
    "    num_bkg_pythia_events[b] = len(loc)\n",
    "    print(f\"Pythia bkg: {len(loc)}\")\n",
    "    \n",
    "    loc = pull_from_mass_range(data_signal, [bands_dict[b]])\n",
    "    pythia_sig_dict[b] = loc\n",
    "    num_sig_pythia_events[b] = len(loc)\n",
    "    print(f\"Pythia sig: {len(loc)}\")\n",
    "          \n",
    "    print()\n",
    "    \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "76f5c1f7",
   "metadata": {},
   "source": [
    "## Pull from sb1 -> sb2 to visualize and get column rescales"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cb31e03e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SIM shape: (991226, 7)\n",
      "DAT shape: (991963, 7)\n",
      "SIGNAL shape: (99952, 7)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzQAAAIxCAYAAACbyP1wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABdZklEQVR4nO3dT3PjVprn+99TUYtZlZnK3Q07rovq2VdL6dl4U9Mm/QIqJGv2Y4vVd12TtBYOZ4YjWqZq1jNFZfULSIvVsx6L2V2b2nRKcu+nxfRE1XKspHtewLkLABRIASQAAsQffj8RjEyRAHgIECCe8+c55pwTAAAAANTRT8ouAAAAAABkRUADAAAAoLZ+WnYBAp9++ql79913Ey37/fff6/3338+9DEVsl7JS1m0va9Q2nz9//vfOuU/TbCfNNSLufTexLu9dv/ded33eO//3LvoaUdR1mW2zbba9mW0/uEY45yrx+PLLL11SaZZNo4jtUlbKuu1ljdqmpGeuwGtE3PtuYl3eu37vve76vHf+6xd9jSjqusy22Tbb3sy2F68RdDkL+eUvf1mLbRa13aLKWgT2a32OV13KGbZumddZv8z3XledPzf7bfPv3VRF7he2vVl13Sd13fZaXMoakKIeaaKyf/qnf0q8LMrH8aqPTR0rbaCFhu9deuyzbNhv2Szbb0VfIzhmzcGxbI40x3LxGlHLFprKRoeIxPGqjyYdqyZ9lk1hn2XDfsumrq1iqBaOZXOscyxrGdAAAAAAgERAAwAAAKDGCGgAAAAA1BYBDQAAAIDaqszEmmWy5xb7mvvSbbAkAJAM1y0AWI7r5PYgoAEAAMBWiQp2CHLqi4AmJPxFXhbVA0BVcN0CgOW4TjYfAQ0AAAC2QlQrDEFO/ZEUAAAAAEBtEdAAAAAAqC26nAEAAGDrkSigvmihAQAAAFBbtNCs4dyqPYjs2FGrAAAAmm+dgf0kCqg/WmgAAAAA1BYtNDmoWktI1VuOsppMJhoOh5pOpxoOh2UXJxeTyUSTyUSdTif1ujc3N3r58qUkaTAY5F00jUYjHRwc5L5dAACKwpiX7URAg1hnZ2eSpFarJUlqt9va2dmRJO3t7c2Wuby81OXlZeHlabfb6na7Ojw8zBTQTCYTtdvtB89v8jOEBYFZOBjp9/sajUa6u7vTycnJ7PkffvhBk8lEvV5vFvzs7e1pMpmo3+8XEtB0Oh2dnZ3p6dOnuW8bAAAgLwQ0iNTr9XR4eDjXcnBzc/MgmDg4OJgFN5vw5MmTzOuORqPIm/NNf4bA6enpXNAi3be03NzcPCjrdDrV/v6++v2+jo+PJanQcrdaLbXbbd3c3JSyfwAAAJJgDA0ijcfjB92g9vb2HnRBarfbmbpLlSGuBaaMzzCZTCTdt34l0Wq11Ov1CmmNiXNwcKB+v7+x9wMAoErsuUU+UC0ENIgU3HAv6vV6hb3ndDp98Pfic1m32+v1Yj9TGUajkY6OjlKvd3t7myoIysPe3p5ubm42+p4AAABJ0eUMkQ4ODtTtdjUcDufGnbTb7dnfwZiOyWSi29tbSV5XqX6/r8lkosvLy1kQcXFxMavpn0wmmk6nury81GAw0GQy0enpqcbjsd6+fTvbzmeffabpdDrb9jKj0UiSdHd3p+l0Otdd65tvvlGr1dLd3d1sXFCn05mNQVn8DIHz8/PZmKG7uzs9efJk1vVq1eeMGqsTdnl5mXpsyng81ng81sXFRewyNzc32t/f197engaDwazl6ezsTO12W3d3d5KknZ2d2ZiiVQP/j46O9PLlS7qdAQC2RlxyAVpnqomABpEuLi7U6/W0u7urVqulTqejo6OjuZvfdrut4XCo/f392XN7e3uz58LZu4JWkvBN9t3dnfr9vobDoV68eKGf//znc9t58eKFDg8PV5b1/Px8NqZE8m7ew4PZj4+PNZlMIsfQRH0GyRuc3+1257qiHR4e6uTkRHt7e0s/Z7/fXxp0TKfTWWAR5+7uTuPxeLb869evZ8HTqmBpOBzO7Y9+v6/d3d3Zset2uxoMBokzmLVaLVpoAABAZdHlDLGGw6Gcc7q4uFC73Z7d5IcFLRiLz02n07lgoN1uazwez9Xyt9ttXV1dxb5/0q5Vl5eXc+M8Op1Oqoxli58hCH4Wx9X0ej199tlnc+tFfc5VN/93d3crg5KdnR11Oh11Oh0dHBxoMBjMArzz8/PIdYIAKBzMSF7rVTiZwt7e3izdcxLtdrtS3fUAAADCCGiwUqfT0WAwmHXJCrptLbMYjAQZs8KigqEsXrx4oZOTE41GI43HY11dXa1sAVlmPB5HBhxPnjx5EKxEfc5VJpNJps/earX04sWLyPFAo9EodnxT0N0uMJ1O9fjx49TvDwAAUEUENHggaKGI0u/3Nz5fyypXV1ezFNOdTidRaudVLQ7rBERFClq4gtYYySvr3t6eLi4uIrvo9Xq9uS5wV1dXD1pxAACoA7KOIQoBDSLFdUna2dlZ2V0qL0m7OXW7XV1cXES2jsR1/1rWLazT6US+99XVVS6fPTw4P62gXOFyBMckSKu92FKzs7Mz66p2fn6uV69ebTxTGgAAQFFICpCDc2tezcB4PH4w2F7yJoMsYh6UNMFIWFTgEWRRC/6/t7f3IIhYFlC022198sknGo1GcwPnh8Ph3KSiWQUZxrLo9Xqzlqgog8FAu7u7c/MIvX79euUxm06nGo/HkYkCptPpxoJYAACSiMtChu1ECw0e2NnZ0YsXL/TkyRP1+/1Z1rB+v6+Tk5O5tM1BauVwSubwc9PpdC7FcbDceDxWv9+fvSZ5Y2HOzs40Ho9ng/Ink4kODw81nU4jt91utzUYDGZpn4Mb+U6no/Pz87kb8RcvXqjf788N+I/6DJIXvNzd3en8/Fyj0UhnZ2c6OTmJXS/uc0ZZ1joSlO/q6mq234N93+121e12Z13+lr1ft9udPXd0dKRHjx5pf39f+/v7c68Frq6u5hIehE0mkwfJIAAAAKrCnKtGhPvs2TP37NmzUt476HsZjvajngPy0u/3dXR0VPjcLpPJRMPhUCcnJ2q1WrOU0ePxWMPhUNfX14nK2uv1cm2lMbPnzrlnadYp8xpRRVy30GRcIxCn7Otc2e8Pz+I1ghYaoAS9Xi9V6uSshsOhjo6OZq1CQba54+PjyKxtUYIJOAEAAKqIgAYoQbvd1uPHj2djfYrS7XbnMqKFBeOLllmWDhoAAKAKSAoAlOTp06c6OzvT06dPC3uPYMzP+fm5dnZ25rqdhVM5RwnGLRVZPgAAgHUR0AAlOjg4mMtIVoSs2x6PxwQzAACg8ghogBK12+3Kjk+JSuEMAABQNYyhAQAAAFBbtNAAAAAAKQTpm8NI5VweAhoAAABUTlTQAEQhoAEAAAASiGqFIfAqHwENAAAAKouuXFiFgAYAAABY03sffx35/J+//XzDJdk+ZDkDAAAAUFu00KwhLhKvCmoEAAAANmPxvqvq94lNQkADAAAAFCQqsKHSOV8ENDmo2pdym2sEJpOJJpOJOp1O6nVvbm708uVLSdJgMMi7aBqNRjo4OMh9uwAAoHix91cfbrYceIiABpH6/b5Go5Hu7u50cnIy99rt7a16vZ729vYybXsymajdbudRzDnT6VTD4XAuGIn7HD/88IMmk4l6vd4s+Nnb29NkMlG/3y8koOl0Ojo7O9PTp09z3zYAAKiWqArvba50LhIBDSIFN/Q3NzcPbsCn06n6/b52d3cz3ZyPRqNCbupPT08fBF+rPsf+/r76/b6Oj48lKXOQlkSr1VK73dbNzU2h7wMAAIqzGKjY85OYJbEpZDlDaq1WS8PhUC9fvtR4PE69/uXlZe5lmkwmkryyJdVqtdTr9QppjYlzcHCgfr+/sfcDAABoOlpokNnJyYl6vZ5ub28TLR+07ATBR55Go5GOjo5Sr3d7e5sqCMrD3t4erTQAAGwx5qzJFwENMut0OrNB+MGYmNFoJEm6u7vTdDqd6+b1zTffqNVq6e7uTmdnZ7NthG/sl62/zOXlZepubOPxWOPxWBcXF7HL3NzcaH9/X3t7exoMBrPxNmdnZ2q327q7u5Mk7ezszPbDqoH/R0dHevnyJQENAAA+e25lFwE1RkCDzIKWjZubG7XbbZ2fn8/GokjeTX94EPzx8bEmk0nsGJpV68eZTqezwCLO3d3drHvcdDrV69evNZlMdHl5uTJBwXA4nCtXMH4oCFy63a4Gg0HiDGatVks3NzeJlgUAANVybgvB17Pk68a1wJAsYD0ENFhbEExcXl7q9vZ2Nial0+mo3++namXJsv7d3d3KoGRnZ2culfPBwYGm06k++ugj9Xq9uYAlMB6PtbOz8+C10Wg016qzt7eXqsWl3W4X0u0OAIC6c1+6sotAcFFDBDTIbDqdSpKePHkiSXrx4oUk74a/1WppMpmsbDkJy7r+ZDLRzs5OytJ7LSUvXrzQ/v6+Op3OXFA0Go3U7/cju6MF3eYC0+lUu7u7qd8fAADU07G7D7x6dJcrHVnOkNnV1ZVardasZeLq6kqHh4fqdDrqdDqzQGeZcEtFlvXXFZQ9nK3t7u5Oe3t7uri40OHh4YN1er3eXKBzdXUV2cIDAACq7dzswSPwxeWJ/vzt5w8e4WAG1UALDTLr9/uzVhXJG0vy9u3byKxhcVm9gvE3WdeXNDc4P60goAq3zuzs7Mz+Pjg4UK/X03A4nHv98PBQ5+fnkqRXr15tPFMaAAColgdja3wEQMWjhQaZ9Ho9dTqd2UD4qDEhk8lk1i0tHDiEg4/g/0nXjxJkGFvnc4TH14QNBoNZNrTA69ev1el0dHx8rOPj48hgZjqdzjK2Rb22aswPAADYnGPnZo/wc6gHWmhy0MTBY/1+X6PRaC7FcuD29nbWNSzQbrc1GAx0enqqbrcrSbNA4fz8fK772IsXL9Tv9/XBBx/MtpFm/UXLWkfiPscPP/ygm5sbdbvdWdKBm5sbnZ6eajKZqN/vz024GSw3GAx0dHSkR48ezYKSnZ2dWVrnwNXVlT777LPIzGeTyWT2GQEAQDMsBkBxLTbIHwENIg0Gg7kb9CSispGFu2oFDg4OIm/0k64fpdPpRHZLS/M5gnEzi8ITh04mE718+VJv3rxRq9WapYwej8fa39/X9fX1rDxv376NfJ+XL1+q1+slKhMAAACWI6BZA7O5VkcwzqXoySqHw6GOjo5mrUKtVkutVkvHx8e6vr5eOtYnEJ6IFAAAFI/WkmZLFdCY2fvOue8LKguQWbvd1uPHjzWdTgsdoN/tdjUejyODlslksjKYGY1GtM4AAFARX3VPvX9zGD5gi+mbn0nDZ2tvFgmkbaG5lvS4iIIA63r69KnOzs4ST+SZRTDm5/z8XDs7O3PdzqK6q4VNp1NNJpNCywcAAOItjnPJI5DJU9y4bHoFLZc2oHlkZv/ZOff3hZQGWNPBwYHG43Fs1rI8ZN32eDwmmAEAoILWCRjclw+zoT1orUGhsoyh+a2ZfSxp6Jz7x7wLBKyj3W5XdnxKVCIEAADQbEnmp4kLqJqYSbcIaeehmTjndpxzR5LMzH5nZr8xs58VUTgAAAAAWCZVC41z7q9C/38l6ZWZvSOpZ2Y7ki5ptQEAAACiJ+ck41r+1k7b7Jz7UdJvJcnMPjKz30n6V0nnzrl/W3f7AAAAABAnbZezpfxWm6GkjyW9NbOXZvY3eb4HAAAAAARym1jTzD6V1JfUlmSSRpK+ltQ2s28k/bNotQEAAECB6NK1fdJOrPmpc+73ob/flxfEHMsLYiaSPpcXuPzoL/adpD+Y2V9LOjOzq/A2AAAAsD3KSGmc5wSaqJ60Xc4GZvb/mtmvzOy1pFtJPUl/kNR1zv2Vc+63oWBmxjn3nXPu15LemNmv1i86AAAAEO3YudkDzZZ6Yk15rTBxrTErOedemdmppH9I+d4AAABoiKgJKYu2zgSaqK4sY2i+k9T3EwCk4ndRu5P0Q4b3rZyqzwJbxoUCAAAA2KS0Xc6mzrknWYIZ30TStaTHGdcHlppMJhqPx6WWYTQalfr+AAAA2yRtC803a77f55Kcc+63a26nUqrWEpJHy1G/39doNNLd3Z1OTk7mXru9vVWv19Pe3l6mbU8mE7Xb7bXLuGg6nWo4HGowGMw9f3Z2JklqtVqSpHa7rZ2dHUmafYazszNdXl7q8vJy7XJ0Oh2dnZ3p6dOna28LAAA0U1Q2Nsb7ZJMqoPEH9Ucys59J2pF0F5ea2Tl3lq54KEsQFNzc3Dy4MZ9Op+r3+9rd3c100z4ajQq52T89PX0QfPV6PR0eHqrT6cyeu7m50eHhoYbD4ey5g4ODzAHaolarpXa7rZubm9y2CQAAgGipupz5g/nj7ErqSjoys69XLIsaa7VaGg6HevnyZabuXXm0giyaTCaS7lthAuPxeC6YkbxWmYODg7nn2u32g+XWcXBwoH6/n9v2AABAM4Szr5GFLR9pu5x1JJ1EveCc+05ewgBJkpm9Y2anzrnI5VF/Jycn6vV6ur29TbR80LITBB95Go1GOjo6evB83Hv1er1CyhG2t7dHKw0AAAVhAk0E0gY0ib85zrkfjS9ao3U6HU0mk7kxMcGA+Lu7O02n07muZd98841arZbu7u5m41o6nc7cDf+y9Ze5vLyMXPbg4EDdblfD4XBu3E673Z79PZlMZgHOYnB2dnamdrutu7s7SdLOzs7s8y628iw6OjrSy5cvCWgAANhSUeOas4y9fi9iQlBSUN9LG9AkPgL+mJrYFh3UX9C96+bmRu12W+fn5zo+Pp69fnZ2Njc4/vj4WJPJJHYMzar140yn01nAseji4kK9Xk+7u7tqtVrqdDo6OjqaC0ba7baGw6H29/fn1g3GCQXLdrtdDQaDlYFMoNVq6ebmJtGyAAAgG7psITKgMbOP5I2HWfyGtM3s77S6paYt6UASgwi2QBBMXF5e6vb2dpZQoNPpqN/vp2plybL+3d3d0qxpw+FQw+FQ4/FYl5eX6vf7Gg6Hc2N5gqxnYaPRSBcXF7O/9/b2UrW4tNvtwru1AQBQVWXM1xfVklGGqFaYLPsjqhWmKp+xSiIDGufcKzObyAtMDiUdywtuTF7q5VUmkn7tnHuRV0FRPdPpVJL05MkTSdKLF97hHo1GarVamkwmsS0nUbKuP5lMIgOSRZ1OR51OR4PBQN1ud2XrT9A9LjCdTrW7u5vgkwAAAGBTYrucOefeSHoj6ZWZnUm6knQrrxtZLOfcj7mWEJV1dXWlVqs1a7G4urrSYDDQxcXFrLtVODVylPD4myzrr9r2zc1NZBexfr+vwWCwNKDp9Xq6uLiYZT8LygcAAJIrY74+xpdsl0RjaJxzEzM7ltQnYEGg3+/PWlUkb4zJ27dvH6ROlhSb7SsYf5N1fUlzg/YXvXz5MjKg2dnZWTm5587Ojg4PD3V+fi5JevXqVWTZAAAAUJ4089BcSnpVVEFQL71eT51OZxYsRI0VmUwms25pweuLwUfw/6TrRwkyj0UZj8ezgCTs9PR05Twxr1+/VqfT0fHxsY6PjyODmel0OsvMFvXaqqAJAAAA60mc5cxvmUnVfmdmf+Oc+8fUpaqZMga9Fa3f72s0Gs2lWA7c3t7q8PBwbiLKdrutwWCg09NTdbtdSfdjVs7Pz2fjbCRvrEy/39cHH3ww20aa9RfFtZrs7OzoxYsXarfb6vf7evz4sSTphx9+0MnJyVza5n6/P5snJ+hWdnR0pEePHs2W29nZ0d7e3ly3s6urK3322WeRrUCTyWT2WQAAAFCMtGmb0xpK+vcFvwcKMBgMUo8XiRqPEjUG5uDgIDIASLp+lE6n86BbWqvVmr3Pssxk7XZ7LpuZ5AUjL1++1Js3b9RqtWapocfjsfb393V9fT1737dv30Zu9+XLl+r1eonKDwAAojGBJlaJS9v8G0k3i60rZnYqqZVi+43ub1PGIDdE6/V6Gg6HuU1iORwOdXR0NGv9abVaarVaOj4+1vX19dIxPYFwwgMAAAAU40FA489BcybpraTHCy8fKl2Qwh0/NqLdbuvx48eaTqe5DNzvdrsaj8eRQctkMlkZzIxGI1pnAADIERNoIk5UC82VpJG8JACLJpIGSeeXMbMf1igbkMrTp09Xzi2TVDC25/z8XDs7O3Pdzha7py2aTqeaTCa5lAMAAADLPQho/MH/n8QsP1F0oBPnTZZCAVkdHBxoPB7PJSzIKus2xuMxwQwAAAV77+Ovyy5C7uLGC9E6tVyqpADOuV+nXD4+NRVQgHa7Xfq4laiEB8AmxWVeZNwfAKCJis5yBgAAABTmz9+mmlWkkuJaYMjwlkyqgMbMfrHw1MQ592/+a7+R1JO0I2ksqe+c+z6HMgJAKsu6ITThhy9OXAtME+fKAgAg8JOUy3cl3Ug6kTTrTmZmX0saSHolqSPphaShmb2fTzEBAAAA4KG0Xc6mkrrOuVfBE2b215KeSho65/429PxrSZ/LC34AYOPCrTFNHDwKAE1DF6vk4n7XmtwTIU7agOadcDDjO5I338zctPLOuR/N7G6dwgHAOuZ+GLun5RUEAAAUJm1A82PEcx1J05jxMqTUAQAAaLi8x+qRpjheXAvMNvdESDuGZu7bZWbvSNqTNxlnlMdZCgUAeTh2bvYAAADNlLaFprXwdzDz4IOp0/1ghy5nAAAAW4L5rlCGtC00b8zsV5JkZj+TN27mrXPu9+GF/PTOXzvnfptLKQEAAAAgQqoWGufcH8zsd2b2e9231uxJsxaZY3lJAvYkOTO7Xgx2AAAAgKS2eWwIkknb5UzOuV+bWV9S2zn33cLLN/4jQJczAAAAAIVJHdBIXkpmSd9FPDeX0tnM/iZ70QAAANBEWeab2cb5VZBM2jE0aQ1WLwIAAAAA2WRqofEH/bcl7SxZrCV/fA0AAACwiLT6yEPqgMbMriT9taQkbYV8SwEAAAAUJlVAY2a/kzfQ/xNJkxWLP5b0PzOWCwAAAGgcex7dJsAcPtllyXL2cdJlzWwxCxoAlCoq/ScDTQEAVRaXRIEue560Ac1tyuUPUy4PAAAANE5cC0xciw2Sy5QUICnn3Jsitw8ASX1xeSJpvjaLydoAID1uwDcnrgUmS9rrJksb0EzM7BfOuX9JsrCZ/cY591/TFwsAAADbhEomZJUqoHHO/cHM/ouZtZ1z/5BglZ4kAhoAAICGSTKInZYEbELaLGen/n87ZvZC0pXis53tyJurBgAqY+7HtXs6e46BlQBQDSRqQVppu5z1JL2j+zlouiuW5w4BAABgy1FphCKlDWjuJD11zv0+ycJm9kP6IgFA/qJ+TL+ivzYAALX3k5TLTySNUyxPljMAAAAAhUmbFCDxpJr+8k/SFQcAAABAVnHZ4po8NiltC82Mmf3CzD41s18tPP+ruHUAAAAAIE+pJ9Y0s19IeiFpT15ygEtJ4RTO35nZ7yT9Lul8NQAAAACyi2uB2Yb5fdKmbX5H0j/KS9f8ibwxNZ3wMs65N5J+7c9XM3XOfZ9TWQEAAFBBaeeb2YabbGxO2i5nX0v6zDn3sXPuD8657xSTmtk591t5aZ4BAAAAoBBpu5z96Jz7Q4rlmR4WAABgS6Sdb6bJA9WxOWlbaP5PyuXfSbk8AAAAACSWtoXmccHLAwAAoCLsOZ1tUH1pW2jMzP7j4nMxC/53Sf+cqVQAAAAAkEDaiTU/N7N/NbNvJJ065/6vFpIC+GmdB5J2nHN/m1tJAaAgX3VP9dVCxh36dQPAPfdlurExwCalnodG0seSvpXUN7OppDszO/Jfa0tqSbrRQjpnAAAAAPmJS5edNjlD3aUOaJxzE0l/ZWZP5aVl3g29PJH0uXPuRU7lA4DC/Pnbz2c/BsHFv2pzI9B/HUDVpJlzpmrXVDRTlhYaSZJz7kzSmSSZ2c8l3TnnfsyrYAAAAMC2iKrAiuvqF9cCk3aC06bIHNCEOefe5LEdACjL7Eegezr7u0pN9vRfB1A1aa6RjEtEkRIHNGb2M3njYj6QN1ZGku4kXUsaO+e+z710AAAAQINFVVjR3TidlQFNKGtZeJB/eC87f7lreZnP/keeBQSAIi3WMC5mOwMAANW2NKAxs1NJT3UfwIwlTeUN/pe8jGZtSU/8x8jMLpxz/6mIwgIAAABAWGxAY2YvJR1KGslreflu2Yb8xACfS/rMzNrOuf+Qa0kBAAAAYMFPop40s68ldSV1nXOfrApmJC8xgHOuJ6+l5q/M7L/lW1QAAAAAmPeghcZvaXkqac859y9pN+icuzGzJ5L+1cwGzrn/vX4xAQAAsEnnZtKz0P+XYL4ZlCmqheZY0lmWYCbgT775W3ld0AAAAACgEFFjaA7kdTdb199JusphO5mtqk2o0hwTAAAAZYlME/zs/r9J75mYb6a6olrRmnK8olpoHucxp4xz7kfNp3cGAAAAgFxFtdBMIp7LaprjtjJbrFVY1XIDAF91Tx/MSdOUmiwAiBOe5DG4X6JHS71F/XY1bcxTVECT57e29DMg6qZE3VPvtY+/5gYFAAAAqLGlE2sCwLb587efP6iVbFpNFgCEpclmBlTR1gQ04ZaYczN95bfSAAAAIBkqeFBFUQFNO8ft57ktAAAAFIjxMqijqIDmkZn9IOluzW3vSGqtuY3SRaUxDA+YA9Bcs64XfosuA2QBwMMY5GqL6jrY5N+uuC5nj/zHupq75wAAAACULi6gaUt6u+a2H0v6X2tuozRRrTCRk04BaJzFWqwHmRIBAKigqFaYbUj0EDkPTR4Ta0r60cze5LAdAAAAAIj0k4jn+jluP89tAQAAAMCcBy00zrk/5LXxPLcFAACA9aTpfkSKZtTF1sxDE+e9j7+WPgz930f2DgAAAJSJbLvJbH1AAwBJBRPyfkXlB4CaWxw83luS+IjrHKpuKwOaY+d0HPrbnp9Ikr64PJndsAAAADRR79nyAAblI9tuOlsZ0ABAGkHtZHhiTfqWAwBQDY0JaLYhxzYAAEBeGIuBpmhMQAMAm3JuJvndU8OtNgBQR7Q4b6+4Y1+3cVONC2iY4RsAAGw7eq5gmzQuoAGAooQrTKgsAdA0dauVR3Zxx7qurXUENAAAAA0V1R2WDGdoGgIaAAAAoOHiuiE2YQxoowKar7qndAMBAAAAtkijAhoA2LRgMt5wZQr90AEAVRHXAtOkxBGNDGi4mQAAAFg+u3xdB4ADixoZ0ABA0YKKk/A8NNwcANi0JtWyA1kR0AAAADSc+/K+0oWeLGgaApoYkTWtH26+HADq4dxM8sfThGtMm5A9BkD1xV1rSNGMbfCTsgsAAAAAAFnRQhNy7JyOQ38HtaxBFiMAWBSuFQ0ynR07R792AAA2hIAGAACgYWZd5z9c+BtoIAIaAMjZex9/PRtP8xWDcAEAOYpKxe2+3O7xmgQ0AAAAFbduN1YqVZBGXIteVb9HBDQAkJPwhZ4xeACAPEW1wiybOHWbENAAwAaQyhlAHrh+oEhxLTBVH4NFQAMAAABsqbjujHUKnrcuoKFpDsAmBD8ESVM5c20CkEXVa86BTdi6gAYAAKCJ/vLhSdlFQI3EtcDUcR61rQ1o1klvt1gbUtWMDwCqIyqVs/Tw+rHtqTcBZPPnbz+XPSegwXba2oAGAACgitatIadiBNuGgGaJoCnu2P87qPn44tL7l3SsAFaJSuUs3V8/Zs8922SpAABoDgIaAACACqpTlimgTD8puwAAsC2OnZs9ws8BAIDsaKEBgBLdJws4mf1NohEAi0jPDMQjoAEAAABqLG4us3USRMwlp6j4uHECGgAowWKygN6Hxb9nET94ADaLFlyUKa6lsOzvJQFNDsIRLP3hAaxrMWUr1xWgmbKmZ46rnMD2iauQWuc7EvWb81XFuzwS0OQgnL552YR5ACBF/NA8m/9zMSX8V2uOqyniBw8AsD2CKUsWg52qjO0ioMkgOJhVj1YBAED1xd0krrrPoLso4CGgyWBWqxnq8z585v177FxlolUA1bbqZuTcjAl8AQBYgYAGAGqAcTXA9qHrOpAMAU0KUbWp9EEHECfP60Pe42oAAGgKApqcnZvNcnUHNarUpALI4tg5xuoBDZA1m5lExSmQBAFNgYIa1cUbEmpVgWZZdcOxzsDdxetFeFxN1E0SFShAfcTdJwB5ifp9amIyCQKanIRvIrgwAc1UtZrSqIQBxyWUA0AyixUOae4XmngTivp5UJFWkcQ1BDQFCGpUF7uckf0MaLZN3HDQDQ1oHnpuIG/bNu6bgGYDZtEsY2uA2lh24S+7pjTq5ieoMKEbGgAgb3G/I0EFW1Sl/SYDdQKanETd/Axjlo3qM0vtDFCOptVY0Q0NKNc6CQAAZENAU6B1+soCqIayW2MA1B+D/9FUX1yeSJq/5y1jiAUBzZrS9FFcHFsj6UG2IrqGAOv7y4feBdaenyRep+6By7JuaFE/Ll9cnkjPii4VsL2yJAtqWosxqqmJmc8IaCqCFM8AyhJ3E1X3HzigavhNB4pRy4CmLv1T4yJgUjwDm7HtN+TLWm28SpTkLVgANmfbr10oRpGZz+buzUO9jzbV86iWAU2TRHVDk+iKBqB47/5pPoHAF5cnOnaObi9ARu99/PXsZo4KS2BzahnQRGXxqZJlEXDafot0RQOQp1WtNl99/LX0ofd8E/tZA3mJ7C2S4v6EigNUSdbrfVSFexnBfC0Dmqax5xYxOPfkQe1pWNn5vgFsp+DaE7TmAHiI8wPYrFoHNHW6gc/Sb/HP334e2xUtSlyavDrtJwCb9/Aa4f0dvqYEmeOCf3sfSr3QNYyWG2yrqHGxaYIZzh2UqchxNbNW/wVF3JfWOqCpu7iL2Fz3tGeLrz4cxLusJUci0AGQTfgasSoF9nsffx05HwFQd3GJiNLMLUP3MtRRURkwi0gW0OiAZlsuIEGNaVr2PLpbW5ZAJ8u+plaqev74xz/ql7/8ZdnFQEmWHf/gfF2sIAm33PT8sTe9FdeD4bP7/zch+OG8yabM/VbUe2eZAwvr4fyrnqgeRtLqMfDrHMtGBzR1lfRGP4+ALSoY2tSFmAHH1cMPw3ZLcvwXKzyyXC96z0L/z7niqYxrCOdNNlULaFa1xCzKUpnIb1z+OP/yl6gH0SrPHm5nVWsmAY1v3aaxIk6Kok60P/7xj2tdGItqvVrV/U1a/iOwLa1qVRB9rP6d/j5DdpJlrXp1/LFZt8zrrF/me2cRvg6teu9NnN+lXENMev78+ebft0aqeCP/oDv2itrjNAFMUZ+3yPO7rtsuUl33SVX29+L1+F353TQjxtasOwzCXEWa+83s95L+knDx9yV9X0AxithuEdssartFbLOo7RaxzaK2W8Q2i9ruprb5rnPu0zQbSXmNiHvfTazLe9fvvdddn/fOf/2irxHL3ntdbJtts+3itz13jahMQAMAAAAAaf2k7AIAAAAAQFYENAAAAABqq1ZJAcxsT9ITSRNJbUkT59y43FJtFzNrSTqWNJJ0J2lHUk/SZfhYJDlWHM9imNmBpGnUvszruFT12FW1XGXhfE1um8+bdcTtt7K+e1n3cROPTd1lOSZJv3fYvGXX2Jjl0x1/51wtHv6HuVx47kJSu+yybdPDPw4u9Hgr6SDtseJ4FnZ8Ov4x6cQcu7WPS1WPXVXLVYF9wvm6ej9t7XlT8H7b6Hcv6z5u4rGp+2PNY7n0e8ejlOMZe63I6/jXqctZT9Jw4bmhpEEJZdl2XUmPJO065x4550YLryc5VhzPHJlZ28yG8i4CdzGL5XVcqnrsqlqusnG+xuC8ySbhfpM2/93Luo8bc2waZJ1jsup7hw1Jca1YlP74lx21pYjubrUQmUlqeR+h/PJty8P/Ui6NsJMcK45nocfoNuoY5XVcqnrsqlqukvcJ52vyfbWV502B+23j372s+7ipx6bOjzWO5crvHY9Sj2nSFprUx78WLTR+n8gH0Z1zbuq/3t58qRAlybHieG5eXselqseuquWqujof801o+nlTprz3W9Z9zLGpHo7Jdst6/OuSFGBHuv8wEdryBg1hM9r+4C7JOzZ37r5JN8mxUoJlOJ75yuu4JFmmjGPHNSIe52t2TT9virbJ717WawDXjupZ95gs+96h+jId/7oENK2yC4CZO0kKXxzM7MLMgudaCbaRZBnkq7XBZcrQKrsAFcX5up7WBpdpmk1/99Ism8d6KE5rjXVXfe9Qfa0sK9Wiyxmqwzk3dc6dLzzN4EmggjhfURa+eygD37vtVauAxu9Xh+qZyGvibQVPJDlWHM/Ny+u4VPXYVbVcFcP5mlLTz5sNKvy7l3Ufc2yqJ8dj8uB7h+pLe7zqEtDM9ZENhD4s/Vs3xMyeRjwdDNx60J85tF7L/+8k4TLIV17HparHrqrlKhXn69qaft4UpoTvXtZ9vHXHpgYyH5ME3ztUX6bjX4uAxh8YNNHDfnU78mYd5YKzAX5miUFEhongSzdJcqw4npuX13Gp6rGrarnKxPm6vqafN0Up47uXdR9v27Gpg6zHJMn3LsdioiBZj38tAhrfWNKThef2/OexAf6XqBfxZepIugllpEhyrDiem5fXcanqsatquUrB+Zqbpp83uSvxu5d1H2/NsamR1MckxfcO1Zf+nEwywU0VHvIitcuF5y61MPEOj8KPw0F4n/vH5VrSXppjxfEs9Bi9lXQQ8Xwux6Wqx66q5Sp5n3C+Jt9XW3neFLjfNv7dS7G9W0nHTT82dX6scSxXfu94lHZMl11j1z4nzV+oFsxsT9KRpNfy+kLeOOeoQdkwP797W9JjeV+6gVuoEUlyrDie+fH7lp7I248H8pprx/IuCOH0lbkcl6oeu6qWq0ycr/E4b7JJsd82/t1btaxf9jeS+i6UDaspx6ZJ1jiWK7932Iwk14q8zslaBTQAAAAAEFanMTQAAAAAMIeABgAAAEBtEdAAAAAAqC0CGgAAAAC19dOyCwAAAABkYWYdeXOUPPb/vXHO9Zcs35aXPvhsQ0XEBhDQAAAAoHb8YEbh4MTMLszs2jm3H7F8W958JqPF11BvdDnDxpjZtZm9DR5ll6dIZjYIf1Y/nzoAAMhPN2Juks8k7fnz0Ujy5joxs6Hu50JBwxDQYJNakg4l/dx/NNmp7j/nRNJOucUBmo0KE2ArdfxWlxnn3FTe7243/Jxzrue35Ew3WkJsBAENNu3Ov7BMly1kZntmNjSzW/9xbWaX/g95K7Rc28wG6xTI3+atmbnQ49bMjhOse7Cw3lsz6wSfcdXnBKIk+e7hgZaoMAG2TVveuBlsOQIaVIofoFxKeiXpWtK+c27XObfvnOtKGkp6EWpKXiuYkSTnXN85tyspGEQ48d/zPMG6I0mP/D97zrlHEc3fQGL+d3tIUJMJFSbAFvF/c+fGw/jncDBWBluCgAaV4d/I3Uqa+hep88UfbOfcxDl3KCm40TiI2FQmoUGF7fBNTUKjJAEQkMAH/r+xWXqQDRUmwFY4lneeMfB/ixDQoBL8TCUX8gKDw1XL+8FHe9VyGQQ3GWlqxz+R1wUEWEvQOiPve9gOD2rFeqgwAZrPH09zotD4GWwHAhqUzv9xv5A3UO+zFKumWTapof9vL8U6XefcTQFlwfZpO+cmuq/NX7uFAFSYAFvkQtJH/nUUW4SABlVwIm9A72maPuH+srk2KfuByY28WtTOquX97EKv8ywDtpPfgjCSZt9tWmlyQIUJsB38tMyfcb5sJwIaVMFT/98swcnLPAviS3PT0dN9rSuwjvZCreJg4V9kQ4UJ0HBm9lTSMBzMZOjaiRojoKkYMzv2s+9ch7PgmFnHf37oz4J7GZ57IJS1Zxhaf+UPZtlCZZxmbCIey+sXv+w9Bv7jqf/vcNmFLtRX/WDZcsFrZBrCuvxWmLnB3v75MBKtNOuiwgRoMP/6eBPRMkOmyC3y07ILgHt+gLLrnOv5J+iFmV1J6kiSc64XWvZC0rWZPZLXx3p/4fWBpEsz2614X9IgKMtUxlDXnAf8wYGX8rIDjUPPH0h6Y2YfLWmaPpd3MTyWdBazzCfyurIA6/ogJiNPX97A9IES3JD7lSD7kp7Iq60895/vyJujRfLmLGlJ6gfff//aE77BfuK/XuuMWjlVmCyd4yWU1vkHSY91v2+nUcs75879rjEHZtaKW44KE2wT/3d5oGRj126cc/v+eh1/nfHCBLNPxASaW4WAplp6uh8MHPyIXkgaRGS4uZR3o/NC0utwMBN6/am/TNwNeRXs+v8WEXRdyhsEvFjzPTKznrx9tx+z7lBeMNNT/P7rJhlgDCzj/yBHtgQ45yZmNpY3G/bBsjSkW1ohsgoVJkDF+d3FBvK6YwbnW9t/RFWqXPrrtbR8rpm5TGf+da0t775o6p/DE+ccKfIbgICmWtqh2rjgRnsck64zCARaoXSgc9vKu3AFa+W5Mf8C2VZ8dqChvBu+vaibDufcjZndSNrzJ7Kbu6jStx056q74Qe3LmzNlVSvNNlaIrEKFCVBhfoXOkaRH4dZI/zd8uixdub+8JX0vApdmYwxNRfg1DeHBv0/8f+NOwKDmMW7AcPBDmviH3J+V+iI0q/XturNkJxCMf1narSODnryL4TTm9WC/PIl5XVre152+7Vib/2O+dDZrP+Aea/VYmm2uEFmllefGElaY7C10gZlZlRyAChNskb68NMvTheePJF1tvjioKwKainDOTRdq+vbkNYVOY1Z54q8X18c9+JFM3AfeORfM0RDcyK+qOc5D0DqS6AbKT4Zw7Qdbi4/wAMC2pDs/EcCDh7z9M5Z0F/decckB6NuOHB0mHKcSnIcnUS/WtUJkA5UoVJgA1Xa4eB7517PI3hNAHLqcVVCoxi7yRsfv99la8Xpb3sC5aYYidJV9EG0qzrmxmUlSK67718Lysz6x/sDaYy1Mlud/fsn7DMu6yyTpShPV152+7Vib3U/2uNKqLpD+eV56hYikUag/fHfVNSTLOimlrjDRfdKEReGue21JE7/ccVZWmEQlB6DCBNsk5nv+ie7PXSARWmiqKbg5j+uK0lnxetAtZRjz+iodSd9kXDeLIFA4Srnetf/v4g1Q7E1EBlG1qN26Z39CJSRtnQkErS2rWj3qWCFSSCVKaP+24rp/LSzfdc7tO+d25Xfzk7cfdkMZ4+YqTJY8usuSOPiCACncukyFCbbdoehuhpQIaKppVe3o4YrXg5vvWVCSdB6L0I/+0n79OTuVl17xaehmITP/JmyiHMYBLPZ1p2878pBk7Mwi/+Z81kqzZNE6VogUWYlChQlQLx3dn39AIgQ0FeN3NwhSCU5jFuvIqx180CQbql0dh7owdJS8D3nwo7+xH1S/nMFNWJqayd0lrw20olbWzJJOWBi+6aBvO/JwmKD2PkowCH1ZK02tKkQ2UIlChQlQE6HKmtQtNCUlNkJFENBUT3AyR97shH78V/V9DwcGvWWpDyPWfxBMmdlgRa3wWvwayUN5tc/XCW88YoMV//PeaPmNXz/JTWU4OYD/9zRB2YBIWVpnAv73daKYVpqaVogUWolChQlQK4fSLNhPpaTERqgIAprqCQa9Z+0uMgn/62f+ipy0L8aeQsGUX2sYzEdRKP9mLbiJuDWzoZnthTOMhcp0Ke/Hf1lA8pGkHX/g7Rz/M6XpUhPcZNC3HevqZWydCSwbS1PJCpEC1kmFChOgNvJICLCxxEaoDrKcVdPNij7UU8X0N/ezhp1J6pvZoaTrpDdPoRrfWbDkl2NsZt3otfLlX4D2/bIcyruxapvZVF6QdifpIsh2ZmZ93d+ELW5r6m/rqXkzo9/J23eSNEx5sRtKepBdCkjD/17v+d/HrFr+v1EZz6pQIZJ2Is4s66TmT3i5K++acmtm5/LO67lgyj9GfXkVJssClo8kvTKz4eLEpBkrTI5FhQm2mF/h0tL6CQE6oqVz6xDQVEzEjN2Lr59pxY//Gk2sXX/90m/ag0AqwXITrbhwrUjdnLQ8N1reBQVIoi+vS1deE1f29fA8qUyFSBHrrIMKE6DSduSdQ1kTkpSV2AgVQECzhfwuXK/k9ZMPBz8HIvc7UJjwPEoFbb9uFSKlVKJQYQJUj39ePlpzMxtPbIRqYAzNdnoir5vHLJuOX6vR1n0WJQB4wMxa/jiUxe5YsRUiWdYBgAxKSWyE8hHQbKeJHg5YvZB0vuZgZQDNl6VChEoUAJtQWmIjlIsuZ1vIOTcxs4l/kk/ldftIlJEHwNbLUiFCJQqAQlUhsRHKQ0CzpRZuIsgGAiCRLBUiVKIA2IDKJDbC5hHQYNOu/YxCcs6tO/ivsvyxAsf+n60SiwLkLkuFCJUoAPJAYiNEIaDBUn4/9yN5A+12zOzlGll99vMrWeWdKjQ2gMnygI2gwgRovmBM3uw3NjQm77CsQqFc5pwruwwAAKzFr7WdaXIlwjZ9VmCRmbUlHYQrV83sVl6LTWTqejO7lDSgO1pz0UIDAKi9bbqp36bPCixiTB6iENAAAACgNhiTh0XMQwMAAACgtmihAQAAQOPknNgIFUZSAAAAAAC1RZczAAAAALVFQAMAAACgtghoAAAAANQWAQ0AAACA2iKgAQAAAFBbBDQAAAAAaouABgAAAEBtEdAAAAAAqC0CGgAAAAC1RUADAAAAoLYIaAAAAADU1k/LLkDg008/de+++26iZb///nu9//77hZSjyG0XvX22zbbL3nbS7T9//vzvnXOfptlummtE0nKwDuuwTjXXqeo1Yh2bfr8y3pPPyHtu6v0eXCOcc5V4fPnlly6pNMumVeS2i94+22bbZW876fYlPXMFXiOSloN1WId1qrlOVa8R69j0+5XxnnxG3nNT77d4jaDL2YJf/vKXtd1+XbddpLruk7puu0qyfM5NrZNFlT8P+4B9kHWdMtWtvFls+jOWsU85js14z7Xfz6WsASnqkSYy+6d/+qdUURyqjePZLEmOpzZQ+1plfOfZB86xD5yL3wdNvEZwvJuB41gNi9eIWrbQbEM0vk04ns3C8VyNfcQ+kNgH0nbtg236rE3GcaymWgY0AAAAACAR0AAAAACoMQIaAAAAALVFQAMAAACgtiozsSbqxZ5b7GvuS7fBkgBomqjrC9cVAGHL7kMkrhnbhhYaAAAAALVFCw3WEq4BWVVbAgBpuC8d1xUASy22xHDN2E4ENJjJ6yJAdxEAAABsCl3OAAAAANQWLTR4IGtrStR6NP0CAACgSLTQAAAAAKgtWmiwEXEtNYytAQAAwDpooQEAAABQW7TQbKFNjmuJa4FhbA0AAADyQAsNAAAAgNqihWaLMX4FAAAAdUcLDQAAAIDaIqABAAAAUFt0OUOpopID0BUOAAAASdFCAwAAAKC2aKFBKaJaYUjlDAAAgLRooQEAAABQW7TQNBytHgAAAGgyWmgAAAAA1BYtNFuCzGEAAABoIlpoAAAAANQWAQ0AAACA2iKgAQAAAFBbBDQAAAAAaouABgAAAEBtEdAAAAAAqC3SNgMAAKCymCQcq9BCAwAAAKC2aKEBAABA5TFJOOIQ0KBy4pqWuZABAABgEV3OAAAAANQWLTSojLgWGAYDAgAAIA4BTYNw4w8AAIBtQ5czAAAAALVFC00DMXgeAAAA24IWGgAAAAC1RUADAAAAoLYIaAAAAADUFmNoaohsZgAAAICHFhoAAAAAtUULTY2RzQwAAADbjhYaAAAAALVFQAMAAACgtghoAAAAANQWAQ0AAACA2iKgAQAAAFBbBDQAAAAAaouABgAAAEBtEdAAAAAAqC0CGgAAAAC1RUADAAAAoLZSBTRm9n5B5QAAAACA1NK20FwXUgoAAAAAyCBtQPPIzP5zISUBAAAAgJR+mmGd35rZx5KGzrl/zLtAwLnZ/BPPYp4POXauuAIBAACgstIGNBPn3F9Jkpl9ZGa/k/Svks6dc/+We+mAFb7qnnr/fvx15Ot//vbzTRYHAAAAG5YqoAmCGf//ryS9MrN3JPXMbEfSJa02yEvQ6tJ7bnN/h8UFMgAAANgOWbqczXHO/SjptxKtNkWw5/HdrLbNe0uCly8uT+b+DlpuAAAA0GxrBzRhzrlXZnYnaSBpYGYjMdYGMZaNiYnylw9PHjz37p+WBy5xQRBd0QAAAJoht4DGzD6V1JfUlmSSRpK+ltQ2s28k/bNotcnMfcmg9yj3gcl8gEJXNAAAgO2QKqAxs0+dc78P/f2+vCDmWF4QM5F3Z3nud0WTpO8k/cHM/lrSmZldhbcBLI6NCVpVgqDkXd23wgQBzKqueEEXtLhtAwAAoBnSttAMzOxS0r6kE0l7um+NGfqJAiI5576T9Gt/nM2vnHP/kLXQqJ+03csAAACAJNIGNI/ktcLEtcas5I+zOZVEQIOl8hjn8iCQ8pMFRLXUMK4GAACgfrKMoflOUn9Za0wcv4vanaQfMrwvGoAJMAEAAJCntAHN1Dn3ZI33m0i6lddFDVuuyPEsUYHTuZm+uDxhXA0AAECDpA1ovlnz/T6X5Jxzv11zO6goxsoAAABgk1IFNM65X8e9ZmY/k7Qj6S4uNbNz7ixd8bANGLsCAACArNKmbT51zj2c3dCzK+mJv9yuvJaYuGXRcOFuXeEuXcwPAwAAgDyl7XLWkZeu+QE/LfN3wd9m9s6KAAg1R/cyAAAAlC1tQJP4DtY596Nxw4sQupYBAAAgb2kDmsQ5d/0xNbEtOphnz+sb/NUtFXOauWkkAjEAAIAqiwxozOwjSV09DGDaZvZ3Wt1S05Z0IKm/dglROrqWAQAAoKoiAxrn3Cszm8gLTA4lHcsLbkxe6uVVJpJ+7Zx7kVdBt4X7sl6tHWFVn88ltiXJD9iYnwYAAKB+YrucOefeSHoj6ZWZnUm6kjcpZmfZBp1zP+ZaQlRG3bqWAQAAoPkSjaFxzk3M7FhSn4AFqzDmBAAAAJuSJinApfx5ZgAAAACgChIHNH7LTKqqdzP7G+fcP6YuFUqRZvA/40sAAABQBT8pePvDgrcPAAAAYIvFpW3+jaSbxdYVMzuV1Eqx/Xb2oqEsaQb/M14GAABUTdT8fnXOJIvlHgQ0/hw0Z5LeSnq88PKh0gUpfHNQG3ETbgIAAKC6olporiSN5CUBWDSRNEg6v4yZ/bBG2VAgJstMLm68EK1TAABUS1QrTFRrDZrlQUDjD/7/JGb5iaIDnThvshQK2KSoLnYEfAAAAPWQJm2znHO/Trk8aZ4rbtV4mW3OZvbF5Unk/tnmfQIAAFA1RWc5AwAAAIDCpGqhMbNfLDw1cc79m//abyT1JO1IGkvqO+e+z6GMqADGiwAAAKCKUgU0krqSBpIu5I2lmUiSmX0t6b9IeiFv7pnHkoZm1iOoKVfSsSB0owIAAEAdpQ1oppK6zrlXwRNm9teSnkoaOuf+NvT8a0mfSzrJoZyNQrYNAAAAIB9pA5p3wsGM70jefDOD8JPOuR/N7G6dwiE/SSfLpGsZAAAA6iRtQPNjxHMdSdOYrmVMrLlE3jPWkmoYAAAA2yZtQDN3B25m70jaU/zcNI+zFApIKqr7Xt6BIgAAAKorbUDTWvi74/97sbigH+zQ5awEzC0DoGnClRdUWgAAwtIGNG/M7FfOuX8ws5/JGzfz1jn3+/BCfnrnXjhJAJCnqBsaki0AAABsn1QBjXPuD2b2OzP7ve5ba/akWYvMsbwkAXuSnJldLwY7qA4SAACounDlBZUWAIAoaVto5Jz7tZn1JbWdc98tvHzjPwJ0OUOtRSZa6J5Kiu66R5AIAACwWakDGslLySzpu4jn5lI6m9nfZC8a4pDNDAAAAPBkCmhSGEj6oOD3QAwG/2cXl1jh3ExfXJ48eJ19DeTvwXn1YTnlAABUW6aAxh/035a0s2SxlvzxNShG0skyAQAAgKZKHdCY2ZWkv5aUpN8Td9wVwLgOAHUUvnbRCgoAiJMqoDGz38kb6P+JpMmKxR9L+p8ZywUAAAAAK2XJcvZx0mXNbDELGlJKkgCAmksAAABsq5+kXP425fKHKZcHAAAAgMQKzXLmnHtT5Pa3SZIEAIyVAQAAwLZJG9BMzOwXzrl/SbKwmf3GOfdf0xerOZjZGgAAAChOqi5nzrk/SOqa2a8SrtJLXyQAAAAASCZtlrNT/78dM3sh6Urx2c525M1VA0nuSzJYAwAAAHlL2+WsJ+kd3c9B012xPHfxCZHNDAAAAEgvbUBzJ+mpc+73SRY2sx/SFwmovgcBaPc0ekEAAAAUKnVSAEnjFMuT5SwlspnVW1wrGscMAACgGKkCmjSTavrLP0lXHKDaogLOJN0FAQAAUIzM89CY2S8kPZF055z7h9Dzvwr/DWyDLy5PIoMdxj0Bqz2oFHg2/3ySlmsAwPZKlbZZ8gIZM3st6VrSuR6mZv7OzH7nBzwAAAAAUJi0aZvfkfSP8tI1fyJvTE0nvIxz7o2kX5vZfzGzqXPu+5zKulWo2QewbYKWmB4TEgMAUkjb5exrSZ/5E2xKkszso6gFnXO/9eetOVmjfI3DeAsAAAAgP2kDmh/DwUwC3L2viexYAAAAQLy0Ac3/Sbn8OymX3xoMcgUAAADWlzYpwOOClwcAAACAxNIGNGZm/3HxuZgF/7ukf85UKgAAAABIIO3Emp+b2b+a2TeSTp1z/1fSXN8pP13zQNKOc+5vcytpQ5HNLH8WkyHJfUk3PwAAgKbJMrHmx5K+ldQ3s6mkOzM78l9rS2pJutFCOmcAAAAAyFvqgMY5N5H0V2b2VN6kmruhlyeSPnfOvcipfFuDbGbri2uBiWuxAQAAQP1laaGRJDnnziSdSZKZ/VzSnXPux7wKBgAAAACrpE0KEMk594Zg5qHwJJrnZkyqCQAAAOQscQuNmf1M3riYD+SNlZGkO0nXksbOue9zLx0AAAC2Bt3EkcXKgCaUtSw8yD/8bXP+ctfyMp/9jzwL2BRfdU/n/ya7GQAAALC2pQGNmZ1Keqr7AGYsaSpv8L/kZTRrS3riP0ZmduGc+09FFBaossguhX4gG5eem2QQAAA8xFQLSCM2oDGzl5IOJY3ktbx8t2xDfmKAzyV9ZmZt59x/yLWkDcENLAAAAJCfyIDGzL6W1JXUdc69SrIh59wbST0zG0oam9l/c879f/kVFaimYxddi3Rupi8uTyJfZ0JVAACAfDwIaPyWlqeS9pxz/5J2g865GzN7IulfzWzgnPvf6xez+iIHsT3beDEAoHHOzWbdN2d/K74yAQCwXaLSNh9LOssSzAT8yTd/K68LGgAAAAAUIqrL2YG87mbr+jtJVzlsp1bCg9jCtYh0MQKAdMItMOHMkL1n/usbLg8AoJqiAprHecwp45z70Ww7ZpI8N5t1L2PyTAAAAGBzorqcTSKey2qa47YAAFvs3T+dksoVAPBAVAtNnr8WW/fLwyBVAAAAYHOWTqyJbBgvAwAAAGwGAQ0AoHS9Z1IvKv09AAArRAU07Ry3n+e2aufP35K1GgAAAChSVEDzyMx+kHS35rZ3JLXW3AYAYIu4L0lzDwBIJ67L2SP/sS5GyAMAAAAoTFxA05b0ds1tP5b0v9bcBlB7kXMTdU8lxSeQoLsiAABAMlEBzSSPiTUl/Whmb3LYDpALixhwzJwWAIAmytJ1k8o01FVUQNPPcft5bqsSImvbgQhRcxIF35+4Hw3GDgAAAKTzIKBxzv0hr43nua0q6T0ruwRII6oVJqq1BgCApvni8mTlMl/53aCBumIempxRww4Aycyul91TSSfzzwHIRbi3AL1M0FQENGt490/UaAAAgHqI6gotSV9RkYCaI6DJAYPoACC7Ly5P1PvQ+z/XU6A8ZN5EXf2k7AIAAAAAQFa00AAAaimc3IMU7EB2cYkDSBaAuiCgAQAA2GKMrUHdEdAAAGrjvY+/1ru6rzX+y4erU9ICTUZmQICABgAAAEuQLABVR0ADAKi88I1TMJfGsXOy57TQABLBBbYbAQ0AAAAeaFqygHAikQAJRZqBgCYGs+kCAIBtRrIA1AUBTQIPayLo4oD1xAbM/nctqr8y3QkAAEgvqhUmqrUG9cXEmivE1U4AAAAAKB8tNCkENeTBIFRqzOsvroamqD61cQFy0GIT9Z0iJScAoIroTYCqoIUGAFBL4a6b52aMfQSALUULDbZSXAsMfWoBAFguKvtZXTOfoRkIaAAAtRLuutnzKyF6z+7/Jg0rUKyo7tNkPkOZ6HIGAAAAoLZooQEA1Na7f/K6uXxxeTJrpQFQnrhENiQLQJEIaJboPfO7M3zo/R1kNwMAAABQDQQ0AAAAWEtUogCJZAHYDAKaBIIuDTSXAgAAPBQ3zxrJArAJJAUAAAAAUFu00AAVEjkxIM31wEpetxavy0swKJlWdaA6SBaAIhHQ6P4mcq6fZ+jHEagCfgyAZP7yoXfttucnzEkDAFuAgAaogLi+x5EtNgBmgoD+3Ey9D0suDFCg2N+Dirfix1W6xVXSAVkQ0EQITr4gTTM14CjTF5cnkQEPPwbAvOEzr3LgvY+/nrXSAE1BtjAgHgENEMGeP6wJo+sKUH3nZnM11kGtdlwrKACg/ghoAAAAaqIpvUbCY92AdRHQACFRrTBRrTUAqiXcAsO8FwCwXQhoAAAb8WBQM2MCgK0XrkgMxoZ+cXkifVlWiVBHTKwJAGis3jPvQUsrADQXLTQAgI0KuocFXcOOnVNvAwEHE24C9fBV9zSy6yjnLuIQ0AA1EDn/AN11gFjv/un+/CCFMwA0GwENUHNx89FQk4WqKyrLUfi7z4SbQH0wCSeyIqABKixq7ozY2aKBCnvv469nrYqzriQEGgBSoAIPcQhofL1nkkROdNQHNVloik1NWhtuERo+855jwk0AqL+tDWjmbvq6pwqCGQBAcb64PJkFEVQeAUiiyAq8uAyIm6poQT62KqCZ66oTM6D63T+d0nSJSFEXPS54QLUdO6dj///vffw1CQJQC7S0p0NXNGxVQBPli0vvxy0YNMqXHwCa76uF8Txc+4HtElchyZxV9bSVAc2xc3PzH0jayBwIqKeoix4XPKDewi014RTPQJUQaC/HWFIEtjKgAZogNtsZ89MAkf787eeM2wG2SFRgQ5DYTAQ0QENxIQceClpcz8387Jaa/R1G9jMAqA8CGqBm4m60mJ8GVfLg+0jLIYANiaq8Cyr5SCDQTAQ0QMMsu5AD8Bw7Nzd28quFgItkAdgErs1APghogDWQvx6YN7tBi2mRqWJXLhIEAMWqQiKdVQkEZteuD+efp1KjHghoAAAbU4Ubm2WCVP6LLTZAkbhpBtazlQGNPbdZBE7GG2RR5fz1kWNpuDnDhgWBwWKLTJWuueHzOO7cDZ9PVWxdAuqkir0XFoPJ4BoVtNyGr1nhFlyC0GppfEAz1z91dlNXnR9UYJMYDIm8ZR38X8UbmyjhlhrG1SArxso0D5lEq6XxAc0y7/7plC8fChFV21v0DVxU7TGZz4B07lM5nzCeBsigCj0V1rFsMu0/f/v50uA07jeX1t3iNTKgmftC+bVrQfcHSer53c0IZrAtmE0ZRXjv469rNfg/reC8CX5TgtYa5qxBVtx31N+yTKJxY++OCy0RpIYGNEBZltXslCm2pYaxNdiAKpwDqyQZTxMgxTPiUEl0ry7dSjeB7t7Fa3RAc+zc7IcmXIPWq8GPK7BJXGyxyrJkE19cnjS6lSL4bF+tuFml5QbbpA4VFXmL606etRdE1OvhHkVhXE+Wa3RAA1RJGXPWxF0AGVuDNJZ1LYtTxjiyvIU/Q9RNy2JXtKiWGyoFtlPTjvs2Bi95yBLoxHVbi6tQadp3LSsCGmCLRdUEHTtHt4ktlzZz2bbWHCZpuaH1s7maep1MG7zUraIii2XdybNU3iwbh5NWmvWafN2pZUDzxz/+Ub/85S/nnotKz9xjvplaiDqeTbJqzpqmZURr+vHMQ1X2UdQx/6p7GhvARHUts+e29Poa912uyj6IkmY8jfTwJuHcbOXEnKR89VT5e5C3Mj9rXi0s2xC8rJLkOGbpkZH2/M8SAKVdJ67Ss4pqH9DMfpAZ3Fxb2/SDllSZ3XXOzWbnU5baZY7napveR2kDlzjBD1keN0d1+54sdj9b5ti5B1mNFrunRdnG1pw6fA9W3QQmPT6b+qwEL8UKjmOWJEDrHpvwe6a5LmRt/Ym6Xq0aS7iOda51tQxopNDBCe3sYFbXRWlOyqIvOEVuv67bLlKV90nai2HhfZifhf94eC6F5+RYvDhW8YYryQ1kef6d/j7hj0LcdS2VZ1FPpt9ulpaYPGU554pYJ9O5+Cz4z8lsXz179kz/z/PnmVpzqqpOA5qz7Neo8zFxDxCTnj9/nvo981T0eVrGvUKV7k9W9chYV+btfJjL2+cuz7m+zFXkImNmv5f0l4SLvy/p+4KKUuS2i94+22bbZW876fbfdc59mmajKa8RScvBOqzDOtVcp6rXiHVs+v3KeM9Nv18Z77np99uW90z7fnPXiMoENAAAAACQ1k/KLgAAAAAAZEVAAwAAAKC2KpUUwMz2JD2RNJHUljRxzo2LWg/FynJczKwl6VjSSNKdpB1JPUmXHNPymdmBpGnSY8G5Oc/ff21Ju/6/Q+fcqNxSbZaZdSR1Jf0gbz9cO+fOyy1VOfzr3cA51yu7LJtS52t8kvM3yTUvr2WwvrhzkONYQ865SjzkHejLhecuJLWLWI9HpY+nCz3eSjoo+/PwcJLU8Y9Hp8jvQFMfkg4k7YX+bkm6lXRcdtk2uA86i98fSdeSnpZdtpL2x1DSRdnl2PBnruU1Psn5m+Sal9cyPHI7rg/OQY5jPR9V6nLWk/fFChtKGhS0Hoq1znHpSnokadc598htWQ121ZhZ28yG8i7OdylW5dyc13bO3QR/OOem8vbF4j5qsqiWiHHM841mZm15rRPbqI7X+CTnb5JrXl7LYE1LzkGOYw1VKaA5kHSz8NyV/3wR66FYax0X59zUOTfJvVRIzTk3cc71XPpuQZybPr9bw5H/b9jYf7296TKVqBvx3HTThaiAjqTLsgtRljpd41Ocv0mueXktg/XFnYMcxxqqREDjXyQe1P76NSCxP/ZZ10OxOC7gOzDP/9xt/7G1nHOHzrn+wtMHkl6WUZ6y+OOIvim7HEgmyfmb5JqX1zKZPwhm4s5BjmN9VSUpwI50f6AjtOUNqMprPRRr3ePS9gdfBtu6q0mXBNzj3FzgnHsU8XRHXpKFrdoXATM7lnTjnDsruywb1nLOTc3ymT28hmp3jV91/gY3qCuuecppma28XuQs7hxM8tulnJbhOOaoKgFNa8ProVitNda9k6Twj5uZXZiZqv6DhzmtsgtQEz1Jp2UXYtP8m9mu5LXalFycjTKzgy2/ljXpGh8+f1sJls9rGaxhxTnYSrCJvJZBjirR5QwI+P2qF8dqMIgOjeO3TtxtYeuEnHMj56VJ7ZvZtZ/atPH8bijTkotRqqZc47f5/K0zzsHmqlRAEzHgrtD1UKwcj8tEXheFvLaHDeGYRfO7p/Scc1ED5LeG3x1jKOlVyUXZlE8cc1BEqdU1ftn5m+Qz5LUMMkl0DnIc66cqAU3Qj3AufV7oixDXzzDreihW5uNiZk8jng4G1TGIrj44N5cbSPqo7EJUxFhSyx+k21h+K9RV2eUoW0Ou8VHnb5JrXl7LIIOE5yDHsaYqMYbGH5g10cM+hztaMmA263ooVtbj4td6DcxstLBMcEHgeNYE52Y8f06f/pLBoo3kn9/Xkj4Kz+cR0tpsiTbuiaRdMzsKPbcnr2ViIOl1DceQpNKEa3zc+Zv0mpfXMsgk0TnIcaynSgQ0vrG8L1v4h27Pf76I9VCs1MfFzxTTizjRO/IyIU1zLyWKxLm5wO93Pwh/x/2WickW/MC1NF9zGQhq5aOCnMaImsfJb634ICKVdSPV/Rqf4PxNcs3LaxmklOIc5DjWUFW6nElSX9Jitpue/7wkr6nOzG79i0ri9VCKrMfzLpyf3W+e7Un6rMCyIrkdRdSkc26uFkpT2zKzPf/RkXS4BcGM/FaZqPlm+pLOtmEfRHhcdgFKUMtrfMLzN8k1L69lkI+oc5DjWEPmnCu7DDN+/8YjSa/l1drdhAdv+Re+N/Kae8+TrodyrHE8D/zlH8u7eR5s6c1OJfjH6UTeMTmQV8M+lnQZdJHh3FzO3z9vY16eOOd2N1icUvlB766kH/x/r6NqTpvMv6HvS/pE3jXuXNIwpite49TtGp/m/E1yzctrGWS36hzkONZPpQIaAAAAAEijSl3OAAAAACAVAhoAAAAAtUVAAwAAAKC2CGgAAAAA1BYBDQAAAIDaIqABAAAAUFs/LbsA2G7+xGR78uYj2JOXoz120ik/d/yBc+5sQ0UEAABAhdFCg9L4wYycc2fOub5zriupbWbXMcu3JV1qO2fXBirHnyQTAEphZtdm9jZ4lF2eIpnZIPxZ/Uk74SOgQZm6ETPmfiZpz59JWpI3S7OZDXU/Sz2Akvnn6DCPoCb0Q+0WHm8XHrdmdkkgBcDXknQo6ef+o8lOdf85J5J2yi1OtRDQoEwdv9Vlxjk3lXeidsPPOed6fjez6UZLCCDOB/6/sV1Ek/JbaB9J2g89ve+cexR++K9fSAoCoM667w2g9u78+4TpsoXMbM/Mhn7FyK3funPpV6i0Qsu1zWywToH8bd4uVNDcJqmMMbODiIqdTvAZV33ObUVAgzK15Y2bAVAjQeuMpHN53UQPVqyS1BP/34lz7mbxRf/H/FzSR/JqZi8XK0UAIMwPUC4lvZJ0La+yZNc5t+93dR9KehG6jq0VzEizSppd3Vf4TPz3PE+w7kjSI//Pnl+hs9ibBQsIaFAa/yQdhZ/za0mCsTIAqqntnJvo/sd67RsAX9Ayu/TH2w92pv6fvZzeG0DD+EHKraSpf89xvtjC4ZybOOcO5VXODOR1b89FKIFRO9wKlNAoSQAEDwENquZYXk3GaOWSADbOv0EYSbMuonm20gRdyC5SrEMrL4AH/C6pF/ICg8NVy/vBRxEtvkFQkmbs3yfyxswgIQIaVIbfdeREofEzAConaJ0JDBb+zcQ//1uStKp7hZ/dp+X/SWsugDl+a8iFvJbcz1KsmmbZpIb+v2lak7tR3W4Rj4AGVXIh6aOFmyUAFeG3wswFG/75OtL6rTRB60ySH/EgeJrqvvYTAAIn8io9TtMMoveXzbWHiB+Y3Mi7Rq5MZOJX2LzOswzbgIAGleCnZf6MGgmg0j6IOUfzGEsTtMy+XLaQmT2VF/xM5Q3una7xngCa6an/b5bgZOk1KKM0rTQ9UVGT2k/LLgDg36AMwzdKZtbiRgWoDr9mMfKH3jk3MbOxvFTsBxnHwAU1l5Hdzfz37/vLjeRVgEwzvA/QKH4q4H15WQKHwUBy/5wJxo7syGux6Ae/tX5LQPgG+4n/eq0zaoVaQaYZe3yMtWKOl1Ba5x/kTfbdkrfvplHLO+fOg/n0lt3fBIkDuLalR0CDUvldVG4ian2PJZ1FrAKgHF3n3LI5Z/ryUqIOlLJWNDx+RlLPzIKXduQN+m/Lm59qpNANGbDt/KBk1znX839PL8zsSn4FgXOuF1r2QtK1mT2SN+h8f+H1gbxU6Ls17/odJArJ9BlCyU4e8K9Vl/LSKY9Dzx9IemNmHy25Pp3Lu7dZdn/zidIlRYGPLmdYmz8J1OIEUnGP69B6HXk3Knf+hFfB41hMoAlUhn+uLh187/+Ij5VtLE2w/NifRDd4HPpzORzKuya8JpgB5vR0nw0raFW4kNc6sXjTHJzDLyS1wsHMwuu5pS0uya7/bxFB2aW8rGmLYwlHkq7k7ds4SbqddeveQlYWWmiwFr+72EDegLegRqPtP6JOykt/vZaW3yDNZTrza47a8i60U7+WZLKixhhAPg4jbn6iBK00J0rXShOc75HXBOfcyD/nL8ysl2ZuBj+4OtJCS8+qa0fW9YANa4e6J+37/45jzpHgRr8VEexIxaQsLlMrz4359zttxadTHsq7Ru3FTAx8Y2Y3kvbMrLMYuJAMYD0ENMjMr7U9kvQo3N/TP+mny246/OUt7vWI5bmJAEoQmsthpVU/2EssHT+z8FpfKQbM+jWno1DlSzdJd5qs6wGb4lcMhhNxPPH/jfu9DLpixSXvCAKixN/zLIH/BioLbv1/l46DyaAn795mGvN6sN+eKD5b49B/9PTwetdT/LHDCuacK7sMqCkzu5RXcztdeP5aZCwDGsHMhglbZ4Llg+5pN865/QTL78lr1ZFzLraSIxgfsGq5JetfSnrinHu0ifWATTMzJ6/nwm7M62/ltc5Enj9mdisvwHi05KY97r2DwD/x+Jss60Rs41befUg4qVBwDZomOW/9czxImrBoEEqy4OQFLcOI5QJdeYkZYluo/e1Iof0cBKdJr7X+fVbtEzjkiRYarCMqmGlJimxuBVAvScbOLHLOjVO20iRpnZHuu6VN05Rn4X2ypELNuh6wMaHMXnFZAoPEG8teb8uriJhmKEJX6bOKZVlnJf8aJEmtuO5fC8vPurj7mciO5bUYHYaeD7rjRY1NCkuSzCgqOQDJANZEUgBkFnPR+0TJJsYDUH2HGVMwB33Mk8xLE9xMrLpuBN1p7tIWxm8FklIGZ1nXA0oQ3HzHfVc7K14PEgEsa31YpiPpmw2sk1QQKBylXC9IXLQYZKW+7iwRlRyAZABrIqBB3g7lZfoAUGNZWmcCfhA0kd9Ks2LxVTdagaCGNEttbnBTk/aGIet6wKatauk8XPF6cHM9CzCSZivMEvhvoLLgVF5r7tNQ60pmfgXuRDkkTvBbjG7kZYTskAwgHwQ0yFtH9zUcAOqrl7F1JhAMbo1tpQnd1ChB7WTL/3cuoPHn1lilI29swTTBsnmsB2yM39U7yPw5jVmsI6+71IOW0FB3s3FoTEdHyQfVZwn8C60s8D9HEMSl6coVOf7IN5DfjS1uATNLmrY+3ErTE91a18YYGuQmVBObuoWGFKlAdfjn8l7CYCFOy/932ViaoFZ4mmB7D2pH/RuLJF1B9pRtot6s6wGbFPz2RlZAhG7A44KHYP3w+d4LjyFJ8P5pA//CKwv8sTSH8lIpX8vrQruqhTc2WHHOnZtZT37Ww5jF+kkG9vvbGsrr6ndOpcn6CGiQp0Np1pyaCilSgUrp634+qby2N5Ykf+LcgeYzCrVCGYQOY64hh5JeSQq6aOwowU1XqKIl7fiZTOsBJVg6j5NWd+uchP/1z9GXKd4/S+C/kcoCfw6rXXnB2q2ZnctrHZkLpvzzvS+vkmXZ2L+PJL2Kyv7o37+kGYMUJAcgGUAOCGiQpzwSAhSS9QRAcuGsPwVs+1wZulf4Qc4jvzW3Iy8bU5Ia5K6/ftquLVnXA8pws+K7OlXMAHy/JeNMUt9v0bhO2t00S+C/6coC/35i33/fQ3kBRNvMpvKCuDtJF8F1z8z6ug8CF7c19bf11G/BvtN9C/Mw5b3LUFKa+bqwBPPQIBehuSTO08xZEbEdt+42AGwffxzBK3njAPqh52/lVZJEzomTdT0AkpkNJD1NMzdUlnWWbOvBPDTbgHloHiIpAPKyI6+WImvKR1KkAljHE3ndWGbZgvxrSlv3aaTzXA/YGmbWMrNrPxgJO1BMz4ws6wBZEdAgF865sXPu0Zq1JKRIBZDVRF6NZbirzIW8Ft9l3WeyrgdskyyBP5UF2BjG0KBK5rKe+F1Bjv3XPpB0um3NygCScc5NzGziD8ydyhsDsxio5LYesGWyBP5UFmBjGEODyvDHz5wF/djDWUT8PPnXkvZJGAAAwGb5CTnaug/8X64KTLKsk7JMjKGBJFpoUBGLWU/8AOY2eD2oRZXX95Z5IQAA2KCFQCRRpsIs6wBZMIYGVbGYIrWl6FzwjzdVIAAAUHnXZvbWzN6WXZAimdkg9DljJwDdVnQ5w0alSZFqZnvhZmS/S1qXJlYAAODfU8yEJ8tsmm36rFnQ5QybFmQ9mWU4CWU9mZskbyGYOZYXBBHMAACArbqp36bPmgUtNNgof2zMgXPuLPTcrbxgJXIyTb9WYjaLLwAAABAgoMHGpc16YmZDedk8phspIAAAAGqDgAaV5s8NMQpSNS+OqwEAAMB2I8sZKstvybmRdGdmLX+szZOSiwUAAIAKISkAKskfa3MR8RLjaAAAADBDlzMAAAAAtUWXMwAAAAC1RUADAAAAoLb+f0ZkjhJ0/c/3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x576 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pulls = [bands_dict[b] for b in bands_to_analyze]\n",
    "\n",
    "in_range_herwig_bkg = pull_from_mass_range(data_sim, pulls)\n",
    "in_range_pythia_bkg = pull_from_mass_range(data_dat, pulls)\n",
    "in_range_pythia_sig = pull_from_mass_range(data_signal, pulls)\n",
    "\n",
    "\n",
    "# list of (min, max) to preprocess every column\n",
    "col_minmax = []\n",
    "for col in range(n_features):\n",
    "    col_minmax.append([np.min(in_range_herwig_bkg[:, col]), np.max(in_range_herwig_bkg[:, col])])\n",
    "    \n",
    "\n",
    "print(\"SIM shape:\", in_range_herwig_bkg.shape)\n",
    "print(\"DAT shape:\", in_range_pythia_bkg.shape)\n",
    "print(\"SIGNAL shape:\", in_range_pythia_sig.shape)\n",
    "\n",
    "feature_labels = [r\"$f_0$\", r\"$f_1$\", r\"$f_2$\", r\"$f_3$\", r\"$f_4$\", r\"$f_5$\"]\n",
    "feature_labels = [r\"$m_{J_1}$ [GeV]\", r\"$\\Delta m_{JJ}$ [GeV]\", r\"$\\tau^{21}_{J_1}$\", r\"$\\tau^{21}_{J_2}$\", \n",
    "                  r\"$\\Delta R_{JJ}$\", r\"$m_{JJ}$ [GeV]\"] \n",
    "\n",
    "\n",
    "# Plot some features\n",
    "#fig = plot_feature_histograms(in_range_herwig_bkg, in_range_pythia_bkg, 40, n_features, feature_labels, \n",
    "#                              third_plot = in_range_pythia_sig)\n",
    "#fig.show()\n",
    "\n",
    "dpi = 800\n",
    "\n",
    "bins = 40\n",
    "sim_color = \"#990000\"\n",
    "dat_color = \"#1c4587\"\n",
    "linewidth = 2\n",
    "\n",
    "fig, ax = plt.subplots(2, 3, figsize = (12, 8))\n",
    "fig.subplots_adjust(right=0.8, left = 0.2)\n",
    "\n",
    "for ax_x in range(2):\n",
    "    for ax_y in range(3):\n",
    "        i = ax_x*3 + ax_y\n",
    "        ax[ax_x, ax_y].hist(in_range_herwig_bkg[:,i], bins = bins, density = True, label = \"Simulation (Bkg.)\", histtype = \"step\", color = sim_color, linewidth = linewidth)\n",
    "        ax[ax_x, ax_y].hist(in_range_pythia_bkg[:,i], bins = bins, density = True, label = \"Data (Bkg.)\", histtype = \"step\", color = dat_color, linewidth = linewidth)\n",
    "        ax[ax_x, ax_y].hist(in_range_pythia_sig[:,i], bins = bins, density = True, label = \"Data (Sig.)\", histtype = \"step\", color = \"green\", linewidth = linewidth)\n",
    "    \n",
    "        ax[ax_x, ax_y].set_xlabel(feature_labels[i], fontsize = 30)\n",
    "        ax[ax_x, ax_y].tick_params(axis='both', which='major', labelsize=20)\n",
    "        ax[ax_x, ax_y].set_yticks([])\n",
    "\n",
    "ax[0,0].set_ylabel(\"Density\", fontsize = 30)  \n",
    "ax[1,0].set_ylabel(\"Density\", fontsize = 30) \n",
    "ax[0,0].legend(fontsize = 19, loc = \"upper right\")  \n",
    "plt.tight_layout()\n",
    "\n",
    "dpi = 800\n",
    "\n",
    "\n",
    "#fig.savefig(f\"paper_plots/features.pdf\", dpi = dpi)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d53df80",
   "metadata": {},
   "source": [
    "## Generate a standard test set and write out\n",
    "\n",
    "This is the first 20k signal, 20k background events IN THE SR\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3391372a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20000, 6)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW0AAADsCAYAAAC/mvfrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQLElEQVR4nO3dP3Pb5pbH8d/Zud02vHKbQgP1KSC5dwFXt8vQcZkmkd+BNU5jsYlHegdSZme2deRJt5WZMXs5LDJ3S2H1BqSLZuuzBR/IsARSgAiSeFbfz4xGwMEfHpHm4eMHR5C5uwAAcfi3TScAAGiOog0AEfnbOh/sxx9/9G+++Wbu9svLS21vb7fetspjH2Neqzw3efXn3OQVT16j0eg/3P1HSZK7r+3rhx9+8Ldv3/qnT5+8ztu3b2vj921b5bGrPHdf81rlucmrP+cmr/Udu8y5P3365JL+00MdXetIe3t7W4eHh3O3P3v27MHnXnTsfedddvtDj+1rXvdtJ69228mr3Xbyqt12eRPwNY607/skWmTe6HzTyKsd8mqvr7mRVzvL5CXp0EMdNV9jy9/h4aEvGmkDAO4ys5G7H0p0jwBAVCjaABARijYARGStRfvy8lKHh4eaTCbrfFgAiFaol9vleq9a/gA8br/8/len5/v5u2/v3ef4+FhJkkiS8jzX69evdXx8rKurKz19+lQvXrzQycmJiqLQxcWFTk5Oao9ZJM9zHRwc6OzsrPXPcLvlb61FG+irrotFG00KC1ZjPB4rTVNlWSZJevXqlSRpMBjcFOIkSbS/vy9JOj09nXvMIkmS6Ndff+0k59rpETPLwtdRJTYMsddtYwDQR0VR6OPHjzfrZQHe29ur3X9vb2/uMYvkea48z5fMdubOSNvMMkkv3P2VmR2YWVpuc/exmSVtYu4+7SRTAOjYcDjUu3fvtLu7q5cvX96MrtM0rd0/TVOlaVp7TNV4PP5qPUkSHRwc6OPHjyqKQqenp0qSRHmeazAY3Izkm7gz0nb3sbuXHx1JKLovJRUhlkvKWsQAoLf+/PNPHR0d6fz8vNGouckx5Ug8yzIlSaIkSTQYDCTNCvpgMNBwONT5+Xmrgi0t6B4J0xtlNgNJ15XNT1rEAKCXyimLLMt0dnbWaAqjyTFv3rzRx48ftbu7q6Iovto2HA4lzYr30dHRnWPvM/dCpLsfm9mZmX1ufdY5ypY/aXZFdJkbrADAssbj8c1oWNLNaHjZY6oF+fj4+Kvplul0+tXx80wmk2p79Ha5UDennUpSmBbJJe1rNuWxFXYZSLoKy01js0el5Q/32GQXBzZv3Z00W1tbmk6nmk6nur6+1ps3b262FUWh3377TXme6/j4WMPhUEmSLDymdH5+frM8HA5v9s/zXGmaand392ba5M2bN7WFvzqwHY1Gl2W8bqSdSSovHg4knUsaSyovpyZhXS1iANA75VRFnfIC4e0550XHlOqmPS4uLiTNRt5//PGHBoOB8jzXu3fvWk2T1BXtU0nfm9m+JLn7B0kys73QWVKUHSFNYwCAmTRNby5GFkWh58+ftzr+TtF290Kzwn07/uAYAGCm/KWch+KGUQAQEW4YBQA9xg2jACAit28YxfQIAESEog0AEaFoA0BEKNoAEBGKNgBEhJY/AOgxWv4AICK0/AFAxCjaABARijYARISiDQARWeuFSAB3beqv9az7r8SgG7T8AUCP0fIHABGh5Q8AIkbRBoCIULQBICIUbQCICEUbACJCyx8A9BgtfwAQEVr+ACBitUXbzPbD11EldlRuq8SGZpaZ2etFMQBAN+5Mj5hZJmns7rmZnZlZ5u5jSftmNpT0KuyXSpK7j80sKddvx9x9uqafBR3a1P0wACxWN9JOJGVhOQ/rkvSTu++EAi5JLyUVlf2yOTEAQEfujLTd/bSymkp6H5aTMApP3f1Y0kDSdWXfJ3NiN8ruEWk2uR4m2AEAt0wmk2qn3Xa5MLd7JEx3TMvpjVCoZWbPQ/Fuje4RAGimOrAdjUaXZXxR90jm7gfSzYXJYYhfaTZlUkjaCrFBiNfFAAAdqR1pm9l+ZWSdSfqs2Ry1JO1IOgmxvRBLJJVz3XUxAEAH7oy0Q5E+MrMLM/uXJIUpku/DaPvC3W+mTcL+xbzY2n4SAHgE6i5EjiX9vSZ++tAYAKAb/EYkAESEG0YBQI9xwygAiAg3jAKAiFG0ASAiFG0AiAhFGwAiQtEGgIjQ8gcAPUbLHwBEhJY/AIgYRRsAIkLRBoCIULQBICIUbQCIyFq7R8qWP/6oL7B5v/z+18Ye++fvvt3YY8eGlj8AiAgtfwAQMYo2AESEog0AEaFoA0BEuGEUAPQY3SMAEBG6RwAgYhRtAIhI7fSIme2HxR13PwixoaRCUurux21iAIBu3CnaZpZJGrt7bmZnYf1aktx9bGaJmaXl/vfF3H26jh8EAB6DuumRRFIWlvOw/lKz0XMZy1rEAAAduTPSdvfTymoq6b2kXYXRdvBE0qBh7EbZ8ieJm0YBwAKTyaTaHr1dLsxt+QvTHVN3n5pZJ0nQ8gcAzVQHtqPR6LKML+rTzsqLkJpNeWyF5YGkq7DcNAYA6MDc7pFKN0im2RTJXticSBqH5aYxAEAH7lyIDEX6yMwuzOxfklR2gIRthbtPm8bW9YMAwGNQdyFyLOnvNfHTh8YAAN3gNyIBICLc5Q8Aeoy7/AFARLjLHwBEjKINABGhaANARCjaABARijYARISWPwDoMVr+IvLL739tOgUAG0bLHwBEjKINABGhaANARCjaABARukcAoMfoHgGAiNA9AgARo2gDQEQo2gAQEYo2AESEog0AEVlr90jZ8vfs2bPyiiiAR2hT99X5+btvN/K4y6DlDwAi0rjlz8zSW+tH4ft+JTY0s8zMXi+KAQC6UVu0zSyTdHYrvG9mF5LysE8qSe4+llSYWVoXW1nmAPAI1RbtUHTzW+Gf3H0nbJOkl5KKsJxLyubEAAAdadM9ktya9hhIuq5sfzInBgDoSOMLke5+LElm9jxMnwAA1qxR0Q4XH6/d/YOkK0mJZtMgW2GXQYhrTkzSl5Y/SbT9AcACk8mkekfU7XKh6Uj7s77Mce9IOgmxvRBLJJVz3XWx2aPS8gcAjVQHtqPR6LKMz+seGUraC9/l7lNJ34f1C3efhljZaVLMi63uRwKAx6d2pB2mQT7cip3W7NcoBgDoBvceAYCIULQBICIUbQCICH/YFwB6jLv8AUBE+MO+ABAxijYARISiDQARoWgDQEToHgGAHqN7BAAiQvcIAERsrSPtWP3y+1+bTgEAJDHSBoCoULQBICIUbQCICC1/ANBjtPwBQERo+QOAiFG0ASAiFG0AiAhFGwAiQtEGgIjQ8gcAPUbLHwBEpHHLn5mlt9aHZpaZ2eu2MQBAN2qLtpllks4q66kkuftYUmFmadPYqn8AAHhMaot2KLp5JfRSUhGWc0lZixgAoCNN57QHkq4r609axACgFzZ5b/yfv/u2k/PQ8gcAEWk60i4kbYXlgaSrsNw0JulLy580uyIarooCAG6ZTCbV9ujtcqFp0X4vaS8sJ5LGYblpbPaotPwBQCPVge1oNLos4/O6R4aS9sJ3ufs0xDNJhbtPm8ZW9hMBwCNUO9J29w+SPtyKndbs1ygGAOgGFyIBICIUbQCICDeMAoAe44ZRABAR/kYkAESMog0AEaFoA0BE1jqnvYxN3ugFAPqCkTYARISWPwDoMVr+ACAitPwBQMQo2gAQEYo2AESEog0AEaFoA0BEaPkDgB6j5Q8AIkLLHwBEjKINABGhaANARCjaABARijYARISWPwDoMVr+ACAitPwBQMQaF20zOwrf9yuxoZllZvZ6UQwA0I02I+19M7uQlEuSmaWS5O5jSYWZpXWxrhMGgMesTdH+yd13QkGWpJeSirCcS8rmxAAAHWlzITIxs0xS6u7HkgaSrivbn8yJ3Si7R6TZ5HqYYAcA3DKZTKqddtvlQuOiHQq1zOx5KN6t0T0CAM1UB7aj0eiyjDeaHjGzfTMbhtUrSYlm0yBbITYI8boYAKAjTUfanxUuQErakXQSYnshlkgq57rrYgCADjQq2u4+DaPta0kX7j6VJDPbC1MlxaIYAKAbbea0Tx8aAwB0g9+IBICIcMMoAOgxbhgFABHhhlEAEDGKNgBEhKINABGhaANARCjaABARWv4AoMdo+QOAiNDyBwARo2gDQEQo2gAQEYo2AESEog0AEaHlDwB6jJY/AIgILX8AEDGKNgBEhKINABGhaANAROgeAYAeo3sEACJC9wgARKzzkbaZDSUVklJ3P+76/ADwmHU60jazVJLcfSypKNe7kP/zvKtTdYq82iGv9vqaG3m109W1vK6nR15qNsqWpFxS1tWJ/+e/P3d1qk6RVzvk1V5fcyOvdvpatAeSrivrT9ocvMwPtejT9b5P3mW3P/TYvuZ133byaredvNptf4x5tal95u4PfqA7JzM7kXTi7lMzyyQ9d/eDyvb/kvTvYfVSlSuiwXZNrMm2VR67ynMvc2ys517m2FWee5ljYz33Mseu8tzLHLvKcy9z7EPOva0vrX7/6+7/kLq/EFlI2grLA0lX1Y3lgwIAHqbr6ZH3kpKwnEgad3x+SDKzLHwd1Wx7XVkehv0WxjaQ11H4vr+OvID/Tzot2u4+lWZvXklFuX6fujf7nNjGi9Cc2NqKUHhuX4QOnbTaoVNOSYXlO508q+zuaZpXsG9mF5pdrF5p11H5+A99zVb9YdIit7V+0M3JKw2POVyUQ0/y2ujzFXJyM7sIXyfzcmibV+e/XOPup+4+dvfTJvvXvdnnxDZehBYUprUVofDcvgqryYIPxrpOnpV197TIS5J+cved8PzMy7UTy7xm6/gwaZJb03zXkNcbd/8gKenLe7IurxDb9PO15e7m7juSXkg66ur52vhvRNa92ecUgI0XoQWFaW1FqBQ+lV9V1tPK40v1nTx1sXXnJc3eXNWRxcryWvI1W+nr2CK3pvmuLK8wij0P24/78p6ck5e04efr1mu35+75nBxa57Xxol26/WaviQ3UgyI0J7a2IlQKv236yswGIbS1YPe1aZJXeHONJT0Jo5SVe+BrVhfbRG51sZXndiuvp5q9Xummn7MGeUmbf77KWCbptwU5tM6rN0W75s1eG+tjXussQtX/Umn2ybw/ZzRb6G4nT11srXmZ2X5l7vFKswvWK8urtMnXrIvcNpFvzb/9q8p1q+HcA3uQV0+eL2nW9lx0+ThrvctfncqczlRf3uzj2zHNf2OvrAg1ycvMCknXYU5tHUUok1T+F3Cg2X8NEzNLwmNuhdzfS9oL+1U7eepi68zrs77M1e5IOgmxleQ153Us1Pw1W9mHSdPcbHYxbW3/xuryCucvX7dCsxHuvBw2mpeZbWnzz1d536XqHPW8HFrl1YeRdqavk87nxOraCVfZYtg0r8+Vx90J66vM61Rf3shy9w/lVyWv2k6eutiG8vo+jIgu1pDXMq/ZqltYm+a27n9jdXl9qDzeQLMP5T68J+vy6sPzpTBgqerk+er0NyIfIvxX4vuwuuvur+piYd99zZ6QxEN3Sl1sQ3ldhxyOV5kX2ln2NVvl6/iA3Nbyb6xBXk89/Kbzpt+T9+S16ecrkXRQuUjZyfO18aINAGiuD9MjAICGKNoAEBGKNgBEhKINABH5P6uvFBICSPX3AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1.08956001e+02 3.95438019e+02 1.83133811e-01 2.34884039e-01\n",
      " 3.13677979e+00 3.52986694e+03]\n",
      "[8.78901978e+01 4.54921814e+02 4.88065094e-01 4.65882450e-01\n",
      " 3.11865163e+00 3.46049561e+03]\n"
     ]
    }
   ],
   "source": [
    "\"\"\"n_sts = 20000\n",
    "\n",
    "sts_bkg = pythia_bkg_dict[\"sr\"][:n_sts,:-1]\n",
    "sts_sig = pythia_sig_dict[\"sr\"][:n_sts,:-1]\n",
    "\n",
    "print(sts_bkg.shape, sts_sig.shape)\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(sts_bkg[:,-1], alpha = 0.5, label = \"STS bkg\")\n",
    "plt.hist(sts_sig[:,-1], alpha = 0.5, label = \"STS sig\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "print(sts_bkg[0])\n",
    "print(sts_bkg[1])\n",
    "print(sts_sig[0])\n",
    "print(sts_sig[1])\"\"\"\n",
    "\n",
    "\n",
    "\n",
    "n_sts = 20000\n",
    "\n",
    "sts_sig = pythia_sig_dict[\"sr\"][:n_sts,:-1]\n",
    "\n",
    "print(sts_sig.shape)\n",
    "\n",
    "\n",
    "plt.figure()\n",
    "plt.hist(sts_sig[:,-1], alpha = 0.5, label = \"STS sig\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "print(sts_sig[0])\n",
    "print(sts_sig[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "44563eee",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(79781, 7)\n"
     ]
    }
   ],
   "source": [
    "#np.save(os.path.join(STS_dir, \"STS_bkg.npy\"), sts_bkg)\n",
    "np.save(os.path.join(STS_dir, \"STS_sig.npy\"), sts_sig)\n",
    "np.save(os.path.join(STS_dir, \"col_minmax.npy\"), col_minmax)\n",
    "\n",
    "\n",
    "print(pythia_sig_dict[\"sr\"].shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "647a979d",
   "metadata": {},
   "source": [
    "## Generate the supervised dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "4f257abf",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABrAAAAFZCAYAAADQPm/jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6OElEQVR4nO3dTY9c15kg6PcIXsyqOi3PqlECXEnUvpSUBrMrt51UL2pXlbT+gJXZP6CaKQ5glAwCQyXd+xnS1atZUUzXbjZiqqqWg6ZJ97oBpQVotqKjvB+cWWREKjLy3ogbEfc7ngcQqLzxcc6NuPeNe897PlLOOQAAAAAAAKAv3um6AgAAAAAAADBPAgsAAAAAAIBe+cG2b/CLX/wi//mf//mt7d988038+Mc/tr1kex/rZHv3Ze/S9l/96lf/Nef8i1tPXuFv/uZv8ocfflj5+cu+57E9v0912bXn96kuu/b8qs/dJOaUXeNsW5dtXtNGGX2t1y7vu3oNq15NxZuxfNZ9fk1f69Xn1/S1Xm29put6bXpPtRhzyt5/6Nv7WCf7MMztfaxT2/tQV7ypoy5VtrdRxhjqNKS6qtP46rrssRsxJ+e81X//8A//kIvYvnx7H+tke/dl79L2iPgsbxBz/vqv/7qwjDLLvuexPb9Pddm15/epLrv2/KrP3STmNL2fm7ymjTI2ec1YytjkNeq13mt2oV5NxZuxfNZ9fk1f69Xn1/S1Xm29put6bXpPtfh+Ze8/9O19rJN9GOb2Ptap7X2oK97UUZcq29soYwx1GlJd1Wm77UOr03zMaWwKwZ/85Ce1PL+u7X0rt8uylbv6fXZtn7v8rNdVlrGvS137XNfzm3zvPtV9k/fvU/2HXPdN3r9P9W/6s2laG59lG59RG/Vq67v2nTT/mjbK6Ovx1aYxHQObltPX42BM382Y9mWTcvq8L+sY8/3m0O+xuyy7zvcfyudd57E0lH2u+722LXPd5w/p/O+y7DHEni7LHtJxNqTvuvJr8gbZ81whk17mX/7lX9Z6fl12rdwuy961crsse8jlRs29d5qsa5eGXP8h1z3nYdd/yHXPuZn6bxJzdi3ezBvLvoxlP3Iez76MZT9yLt+XpuLNLnx2QzSmfcl5XPuzC/vS9D3VkO83h1b2Lu5zl2Xv4j5vW/YutOHsYtm7uM/KHka58zEnXf29uc8++yx/9tlnW70HsHtSSr/KOX+27uvEHGATm8Qc8QbYhHgDtMU9FdAW8QZo03zMaWwKQQAAAAAAANiEBBYAAAAAAAC9IoEFAAAAAABAr0hgAQAAAAAA0CsSWAAAAAAAAPSKBBYAAAAAAAC9IoEFAAAAAABAr0hgAQAAAAAA0Cs/6LoCTXvvo88Lt3/75act1wRgXMria4QYCzCzLFZGiJcwNs55oCvav4ChW3UdFSGm7SIjsAAAAAAAAOiV0Y/AmpllZ6tkcgGobr73ixgLUGyxp6B4CePmnAe6ov0L6LtNRqyLabtrZxJYAADQN0U3YqbFAAAA2JxpncdDAguA2pl/HQAAAIAy2oioQgILAABaZloMAABgzPpwf2Na5+GTwAKgNmW9Z1wgAABDZQoaAADohgQWAAAAAABQO5192IYEFgAAAKxgChoAgG5Zc333SGABAACw8zZNSElkAQBAMySwAACgR/QqBACWkTgHdo0113eXBBYAAABMVU0WSyoDAPSDRNZ4SWABAEAP6FUI41d0PkuEAZsSPwAYOwksAAAAAACg1KqOdV0k1SXyx08CCwAAtmSUFLBMUePKLG5Y9w4AAIpJYAEAAAAAACstdrJZ1SkHtiGBBQAANTFiAqjKuncAALCcBBYAAAAAALA2nfhokgQWAAAMgHVyAAAA2CUSWAAAANAzktYAAOw6CSwAAOgx6+QAAABtcZ9Bn0hgAQAAQE9IWgMAwBUJLAAAAHaKZBAAwHKmLaYP3um6AgAAAAAAADDPCCwAAAB2kp7FAADQXxJYAAAAAADATiibTlrnpv4xhSAAAAAAAAC9YgQWAAAAAAAwamUjrMpGZNE9I7AAAAAAAADoFSOwAACgIj3zAAAAoB1GYAEAAAAAANArRmABAMCayuZOBwAAYJjKZtxw/9cdI7AAAAAAAADoFSOwAAAAAACAnVQ2wsoayN2TwAIAAGCUNDoAAMBwmUIQAAAAAACAXjECCwAAgFGz8DYAwPeMUmcoJLAAAABgIIoanCToAAAYIwksAAAAAADYMTrB0HcSWAAAANBzRQ1Mpv8BAGDM3um6AgAAAAAAADDPCCwAAAAAAHrnWUpLHz/OuaWaAF2QwAIAgAWm5QIAAIBuSWABAAAAANBbiyOtVo3MAsZBAgsAAEp8++WnXVcBAAAAdtI7XVcAAAAAAAAA5hmBBQAAAAAAI2NtX4bOCCwAAAAAAAB6xQgsAAAAAAA68yyl1t7rOOfayhoKa/tuZ9lINp9tsySwAABgwMpuptxIAQAAMGQSWAAAAAAAdK7O0VGL71XnKC92w7JOgdYXa4cEFgAADNDsZmrxRvzRvce3tu/iNCkAAAAMmwQWAAAAAAA7oWwklk5f0D8SWAAAMAKzG+5H06ksjnM2TQoAAKPmehfGTQILAACAQbMGAQCwStkIK0kw6C8JLAAAAAAABsN0f7AbJLAAAAAYhW+//LTrKgAAsEPKZgJwXVqPd7quAAAAAAAAAMwzAgsAAEaubF5/U68AAACsr2yElbVZ62UEFgAAAAAAAL1iBBYAAAxA2SiqZcpGWG3yXgAAANAmCSwAAHaaKR4AAIAhc0/DWElgAQDAgNS5bpW1sQAAAOgrCSwAAIjyRXgBAACGwD0NYyOBBQAAO8baWAAAsJlV18xmM4D6vNN1BQAAAAAAAGCeEVgAAAAAALCGxZFWZjOA+klgAQAAAACw08oSUE1OCWg6QlhOAgsAAAAAAGrQRSIMxkoCCwAAAACAnVSWWGpzSkDTEUIxCSwAAAAYsPc++rxw+7dfftpyTQBgd/UhEQZjI4EFAAA9UtcNblGDtsZsAABYT13X5xJZsD4JLAAAABigsqR02YgsAAAYEgksAADooU0XeS5q0NaYDQAA69n0eryp94Fd9E7XFQAAAAAAAIB5RmABAAA3lM3Pr/coAAB0x8wK7BojsAAAAAAAAOgVI7AAAKADZaOculQ2wqqPdWV36XkMAOy6onVvYYwksAAAAAAAAGpS1ulK8nE9ElgAANAh60rBZtz8AwDAuElgAQAAAAAAbKmsk5VpsDcjgQUAAAAAAAO1as1asz4wVO90XQEAAAAAAACYZwQWAAAAAAAM3OJIq1Ujs6DvJLAAAABgQFY2Rt173E5FAIBG3frN9xs/eGVrYZWtnbXrTCEIAAAAAABArxiBBQAAAANkmiAAGKfF3/iZRyWjd+i/shFWZSOyuCKBBQAAACO2KrFV1kgGAPSTTivsCgksAAB2gp5tAADALipLeOnEQt9JYAEAAMAOMOUgAIxDk4knI7fpEwksAAB2Stnc4wAAAGNSlmzSiYWh2NkE1uIUMhoyAAAAGBONUwDApozcpg/eKdqYUnredkUAAAAAAAAgonwE1hcppfcj4ocRcZlz/qa9KjVrcaSVxbwBAGiSnorAptaNH7Oe0o+m97lV16hYVo51LgCgfX1usy67bnDNQBMKR2BFxMuc8+/jKoH1LKX09ymln7ZYLwAAAAAAAHZU2QisNyml1xHxPOf80WxjSulvc87/1E7VAABgPPRIBDbVVPxY9r5GjwJA9xZnE+tS2XWDawaaVJbAOs05/zYiIqX0VxFxGRF3IuJ/iQgJLAAAABrV56lz6qbhBwAAbitLYF2nU3PO/31u5NXv26kWAAAA0DXrXADAeOlEQ9/dSGCllP4uIu5FxAcppZPp5klcjcAy8goAAHaYhmy60Kepc5rmXAIAxmJVcsx1D1XcSGDlnH+bUrqIiA9yzl91VCcAAACgQ9a5AIDxkjxiKG5NIZhz/reIuJG8Sin9OOf8TVuVAgAA+kNDNgAA1G8X1vxcvJdwD8E6rhNYKaXnEfFJRNyJiLOI+OPsoYh4PyL+svXaAQAAAAAAvSYxRRPmR2B9mnP+U0ppEhEnOec/zB5IKb3fes0AAAAAAGDkdmnNT4qVjcbb9WPjOoE1S1jlnP+QUvppSilHxF5EHEbEeTfVAwAAAAAA+sh6WjTp1hpYMznnb1JKr3LOH6aU/i4ivmmvWgAAADAeY5tWp2x/NGIBAFRXNsJqF9ZHq6IsgZVSSj+NiK+mf7sCBQCAgTMtBQAAAENROgIrIj6KiMfT0VcfRMQ/tVMlAAAAGKehj1Aqq//YRpgBAO2qci0x9Oso1leYwMo5fxXfj776bUrpZXtVAgCAzZlq4TbTUsBuMuoSAIAhK0xgpZT+KiI+nv0ZEe9HxH9sqU4AAAAAAMBIrbOeppHeu6tsCsF7EfFs7u+jFuoCAAC1McIA2FVGXQIAMAZlCazXOec/zP4whSAAAAAAALCNPq5jtWqEVx/rvCvKElifppTOIuJtXE0h+BcR8Zet1QoAAAbGtBZAhFgAAAB1KUtgneWcv5r9kVL6WUv1AQAAAAAAWKrukVOLz9cxqXtlCaz9lNJhRHwXV2thGSMHAAAVmF4CiBALAAC6VpaAcp02HGUJrK9zzr9JKb2fc/5TkmkEAAAAAAB6xsip8SpLYN2dJq32Uko5Iu5GxD+3VisAAAB2wnsffd51FQAAGICqiamyEVaz10twDUdZAutZRDyMiP2I+G8551+3VyUAAAAAAAB2WWECK+f8bxHxact1AQAAYEd9+6VbUAAAbqtrzSprXw3PO/N/pJR+llL6XUrpbUrp/0spvUop/W1XlQMAAAAAAGD3XI/ASil9EldTBt7POf9huu3fRcRxSukXOed/7KiOAAAA0CvWTgAAgGbNTyG4l3N+OP/gdCrBX6eU/nO71QIAAIairCHfFB2we8QDALjtvY8+77oKMEjzCayvlzxv2WMAAACwkyRmAADGTQed7swnsD5MKV2WPO/DiPinFuoDAAAMRNkNm6nVYPeIBwCw2rdfftp1FWBQ5hNY9+NqDayiq8v3I+JhwXYAAAAAAIBR0UGnezcSWDnn3xc9KaX0fkv1AQAAAAAA6LVliay6phcsWz9tV0bzXSewypJXqx4DAACAsdLDFhiyoobPXWn0BGD4frD6KQAAQISGbAAAgF23bHRVXfeMZZ0NykZkjZUEFgAAAKxQ1zQwAG0oavjctUZPAIZPAguAldzoANykIRsAAIAyZSOx3Euup1ICK6X0ZznnPzVdGQAAAAAAgDFad4rBXU94XSewUkp/X/KcFBGHEfEfW6kRAL1lsV9gV1jrCgAAgHWVJZzcY25mfgTW/xwRz6f/fxgRF+1XBwAAqjG9KQAAAEOw7kiqVQmvsvvhsXU+v05g5Zyv9yyl9MOc8+/n/267YgAA0LVdn64BAAAAulK2BtbddJXhexsR70bEQUT8c1uVAgCAqsbWwwzohmldAADou7L737HOUFKYwMo5/zql9ElE3I+I1znn/9JutQAAAGAzklEAAIxB5evae4+3f48FfZiRpDCBNU1e3YmIryPiRUrppzlnI7AAAAAYtT7cqAMAAOVTCH6dc/5NSun9nPO/Jb3XAAAAGBjJKAAAhmjd69hHFaYQrPqefZrNYNUaWHsppRwRd8MaWAAAAGxhrHPzAwAA9StLYD2LiIcRsR8R/y3n/Ov2qgQAAAAAAMA6+jR6qg5lCawvIuIk5/xNSunvUkrPc84ft1kxAAAAxunbLz/tugoAAEDPLRuB9SSl9HXO+aE1sAAAAAAAAPprbGvAliWwcs755ymln6WUvoyI1xHx2xbrBQAAAIxE2XQ2Y2tkAQCgPmUJrDsppZ/mnL+KiK9SSlbaBQAAAAAA2EGr1tdqomNSYQIr5/zrhU3/Z+0lAwAAwBbGtkj1GJU1ZPjuABib9z4yBgTqdp3ASik9j4hPIuJORJxFxB9nD0XE+xHxl63XDmALXfQKAIAhKrvZ/vbLT1uuCQAAAH222KbaZMek+RFYn+ac/5RSmkTESc75D7MHUkrvN1YDAAAA2IKOSQBAX+gExlj0YcT8dQJrlrCaT1zNcTcADFabvQIAYEjKbq5NfwIAAEDX5qcQ/DK+nzYw4mrqwBymEASgJkUNonomAVVIqAAAAMByy+6dq7bB9Wl2g/kpBM9yzl8VPckUggAAwLqWjXru000RAAAA/TM/hWBZ8urHEfEXEfH7luoEwMgU9fAwmgLYRN2jNk0rCwAAwNAtu1cechvcD4o2ppQ+iYiTiPgurqYQfB0R/9RivQAAgIFaNrpK0hAAAGB8yu71tpl9ozCBFRGRc/4gpfSznPNXKaWfbVwCAAD0nOnsAAAAoF/KElhvU0p/HxF/SCn9IiL2IqJwikEAAABgOIqmkal7itaqmuipCwBAe8qu2+qYfaMwgZVz/m1K6f2c8++n0wm+2bokAAAA2ICpJwEAYPeUrYH1ZxHxw5TS30bE24i4HxH/3GbFAAAAgPoUjbLqalHvJnvqAgBwW9l1X1cj8asom0LwSUS8nvv7Ry3UBaBVpisBABgW12kAALA7yhJYL3LO12tepZRetlQfAAAAAADopa5GL8OmykZYDeFYLktg7aWUnkfEZUSkiPhZRHzYWq0AGmS6EoDdI8YDAADAsJQlsPYjYj4t910LdQEAAGCHSTYDAEPR53WDYCzKElivc85/mP1hCkEAAMbA+jkAAAAwDGUJrE9TSmcR8TauphD8i4j4y9ZqBQAAwM6SbAYAAMoSWGc5569mf6SUftZSfQAAAAAAANhxZQms91NKr3LOf4qImE9mAQAAAADdW7V2oBGtsLn3Pvq86yrAzitLYF3OklcRESmlv8o5//d2qgQAAMCQafABAAC2VZbA+k/TNbDexNUaWO+HNbAAAACowaoRAwCsZ3GklTgL9fn2y0+7rgLsrLIE1tOc829nf1gDCwAAgHWVNfg8Sw9brgkAAFCkbPaEPiRvCxNYOeffppT+c0R8EBGvcs7/pd1qAQAAMHbWZgEAAMoUJrBSSr+Iq+kDzyNiP6X095JYAAA0zbo5AAAA0LyyEVZ9ui8vm0LwDznnr2b/n8ybCwAAAAAAQEvKElj7KaUcEZcRsR8R70fEVyXPBQCAWvVhrm0AAACgO+8Ubcw5/yYi7kbEk4i4Z/pAAAAAAAAA2nI9Aiul9H9ExOncY8+m/+WU0p/lnP/UduUAAGAdz0x9Xauyuc+NkAOA4Sr6fffbDkAfzU8h+KOIeDciUkTk6bZPI+J+XCW2/rHdqgEAAAC7qKxDwnHOhdsBABif+QTWJznnf4uISCm9HxG/iYhXEbE/2w6wC5b13nfDDDAM4vV2ynphl43IAgD6r+j33W87AH12ncDKOf9bSunHcbXu1b+LiPs55z90VTEAAABgt5R1QDBFLADA7plfA+vziPi7iDjJOf/z/JNSSj9d3AYwNst667thBgBYn2soAABgU/NTCB5GxH+KiJRS+ml8vxZWiojPI+LD9qsHAACMkfVtAAAAWGZxDazfFz0ppVQ8CT4AAACsIDEJAACsa34NrMLk1fSxr9qpDgAAMGbWtwEAoE/e++jzrqsAlHin6woAAAAAAADAvB+sfgoAAAAAAIzXt19aRQf6xggsAAAAAAAAesUILAAABsd6SQAAADBuElgAAAAAAABce++jzwu3tzndpgQWAACDdZxz11UAGIU+NFAAAMA8CSwAAFpX1lAKAADQBPcgUE1ZB6YuziEJLAAAANhRfWqgAACAeRJYAAB0xtRUAABAm9yDwHBIYAGs4VlKhdutwQIAAAAAUB8JLAAAAACgUWUdQmd0DAUYhrKpppsY3SiBBVBB2YX0qgtwALYjzsIw3DpX7z3upiIAAMBoSGABAAAAAFur0vlosYOoDksAw1A2wqpsRFYdJLAAAOg9U8rAMMzO1UfTm1jnLsAwFDU+NjEVFACsQwILAAAAAKjMelYAtEECCwAAgI1c99ifrnn1qMHpQwCoX9EoqyanglqmLCkmGQawuySwAAAAAGCEmk4KSS4B0CQJLAAAALbyy5cPI0JDJgDrK/vtWDVNIQDjJ4EFAAAAACMiKQTAGEhgAQAAAIOwrPHdCECon4QXAFXdWkNxuk7u8RbvKYEFAAD0hgXcAQAAiJDAAgCgY3r2ArDKsiS23xHY3o3zaNpjfn6bjiQAlPn2y08Lt98akbUBCSyAGugtDgDbsVYHAAAA8ySwphazgWVZQwAAmiHpDwDQjmWdMB9N28hcmwHQtcEmsFb1xPQjC7RBb3EAdlHRVBA6gAEAAFCnwSaw6rJ4o13HvIwAAAAA0Dc6fLMLtO/CeAw+gbX4w2vUAwAANKNolJUGAgAA2uYaFHbD4BNYANTLRSAAAAAwdKa4huGTwAIAAGClwtku7j1uvyIAAFOSVDBuElgAFHIRCAAAsLuKZudwnwhAmySwAABojGlJYXzm1yF+ND3HF9cmBqAerqUA2GW9T2AVTlMBAMAouNYD6LeyxnOjMGC8is5viTQAutD7BBYAAMP37ZefSlYBAGxol5PGZdeQRv8CjN9gElh+lIAhcqENcJsYCNB/ZY3lRmEAANCWwSSwAAAA6M6je4+v/pXAoKd0HoNxKTt3jeoH2B0SWAANcKEN7LIbsW7a4C3+AQAA6zLyF3abBBYAAACV7fI6LPSTzmMAAOM02gTW9YVqQa9f0wcAADTvOOfrqcZcfwEAAJvSgQZ20ztdVwAAAAAAAADmjW4E1mLv3vlev0ObPmAo9a2jR/XFxUW8fPkyIiJ+9KMfxd7eXhwfH8eTJ0/iwYMHERFxcnISFxcXcXl5GUdHR/Hxxx/H0dFRREScnp7G+fl5XF5exuvXr+Pg4GBpWffu3YuDg4M4OTmJ4+PjresPAABjUXofMp3dAgAAoA2jS2AxPOfn5/Hq1as4Ozu73nZ5eRknJydxeXl5ncB6+vRpPHv2LE5PT+PFixc33uPs7Cw+/PDDeP78eTx9+jSePn1aWt5kMomIiN/85jdLE10AAAAAAEA3RpfAem864mpM+rpmRF0jxE5PT+Prr7++sW1/fz/u379/I6lVxcnJSdy/f780gTWZTGJvb2/TqgIAS1xfh01HaTwa4XUZ3Sm79uzrtfIYlM1uAQAA0AZrYNG5t2/fFm4/PDxcO9l0eHgYEVfTBBa5uLi4fg506VlKhf8BAAAAADDCEVgz3375addVYA1liaWTk5O13+vnP/95PH36VKIKADryy5cPI8LIGOpRdhzp+AEAMB6PzOIAFBhtAovhePjwYdy7dy/Ozs7i+Pj4xqirTZJQJycncffu3VvbLy8vY39/f5uqwtY0wgEAAAAArCaBRecePHgQEVdrYZ2ensb+/n4cHR3FycnJRgmng4OD2N/fj/Pz8zg6Orre/ubNmxt/AwAAAAxBXzo9Fq09bxYk6uR4AuZZA4teePDgQeSc4+XLl3F0dBTn5+dx586dePPmzUbvd3R0FE+fPq25lgDAPOv4AQAAAE0xAoteOTw8jMPDwzg7O4uTk5O4f/9+fP3112u/z8nJSdy5cycmk0ns7e3Fmzdv4uDgoIEaAwAAADSvy/VFi0bFFI3GAoA67eQIrMXewnoNd6tslNXTp0/j8vIyJpPJ2u+5v78fBwcH8cUXX0SE9a8AoEnHOd9oUFn8GwAAAGBdO5nAoj8mk0k8f/689PH9/f14+/btRu99cnJyPY3g3t7eRu8BbZNgBwAAAADYsSkEy3oCaxzu1vn5eZydnZU+vunIqZ///OdxcnIS5+fncXh4uGn1AACooGwaIQtxAwAAsAkjsOiF09PTW1MFnp6exunpaeX3uLy8jMvLy+u/9/b24vDwMJ4/f24EFr03m25r8T8AAAAAgF20UyOw6Kezs7M4PDyMx48fX2+bTCZx7969ODo6ut52cnISFxcXMZlM4v79+/Hxxx9fP356ehrPnj2L58+fx8OHD6+3n5yc3Cjr4uLielrBTz75JE5OTuL4+LjpXQQAGK2yEVYWdodxKzrHjbgEAKBOElgDMOYpDvf29q6TTcumEYyI68RTkbOzs8LXzyfAIiIODw9NJwgAaxrztQgAANtZTGhLZgNQFwksAAAAoJKihmkjLgEAaIIEVo9Z/wYA6BPXJgAAzCwmtCWzAaibBBYAAAAwWmVT4eqYAQDQbxJYAAMw5ptuC4ADAACwrjHfJwNwRQILYIeZ4gEAgLEqa8Qua/QGAKBfJLAAemzMN90WAAcAAGBdY75PBuAmCaw5z1KKuPf4+/+fMvQYGDtT9gHrkGyGcShr6Hs0vSd65FwHAAA69E7XFQAAAAAAAIB5RmDFzRFWs16Gxzkbegz03rI4ZfQosK51r32M3oRxWLxmmN0TOccB2ETRaH2/KQBsQgILAAAYrLLEq44cAAAAwyaBBTBAyxrljB4FtqXhHwCAdRWNsrJ2KgDbkMAq8d5Hn0fMLV5sqDMAAPRHWaJVRw4AAIBxeKfrCgAAAAAAAMA8I7AWzI+0epZSPJqOwurCUIZZbzM67fT0NM7Pz+Py8jKOj49jb28vIiImk0lcXl7GvXv34sGDB5Xe68mTJ3F4eBgHBwc3tl9eXsadO3fi4OAgPvjgg9jb24vJZBLPnj27fv6svIuLi3jx4kUcHR1tvE/QF9YEAQCAcq6XAQD6TQKroqILWxe12zs7O4s7d+7E6elpPH369NbjT548ibt378br169XvterV6/i66+/vvU+k8kkjo+Pb21/9uxZnJ6exuHh4Y1tl5eXG+4NAAAAAIuKOmlbrgOAVSSwBqCvP+htjBB78OBBfPfdd3H//v148eJF6fPevHkTDx8+jLt3795KVL19+zZOT08rlXd8fBxPnjzZqs7QNWuCANAnZdeMfb3GBTY3lPPd9TIAwDD0JoHVxwvF45zj0fQCfP4Ct491HbOzs7NIKcWbN29uTQ847+DgIPb39+P8/PzGFICTyST29/crlzebxhDG5L2PPo+YTon6aCDTkwLdG8p0xgAAYzDGa6+iBPYY9xOAZvQmgQXLHB4exuPHj5eOwoqIODk5iefPn99IYK2TvIqI+OCDDzaqIwzZs5RMiwpArcpGXGi0gvFxvgMA0ITeJbA0oFLk4OAgzs/PCx+bH5l1dHQUp6enMZlMrkdSLRu1VVYWjNUvXz7sugrAAPVt6icAgDFbvPYyExAAu6p3CSwocufOnbi8vFz5vP39/Tg4OIgvvvgijo+PW6gZDEtRJwE3Q7BbnPMAAADAELzTdQWgbicnJyunGgQAAAAAAPrLCKwtlPVgNg1i/SaTSeFaVufn5/Hq1at4/vz5jedeXFzcmEYQWE1Mg93i3AYAoEtF6+SZuhqAeRJYDMKrV68KE1h7e3txdnZ2a/vvfve7ePbsWTx48KCN6kGvWCwbAAAAABg6CawNlPVYtqZEcy4uLm5NCziZTOLdd98tfP7HH38cz58/l8CCCsQ0AAAA2lI0ykpHTACKSGDRe0+ePIkPPvggDg8Pb2y/uLiIo6OjwtccHR3F6elpXF5eFo7cgl1Q19QLphYEYIj8fq326N7jq381GgL0is6E23ENADAeElj02pMnT+Lly5e3Rl9FRLx8+bI0gbW/vx/7+/vx9OnTwikGJ5PJjX9hqPRSA8po+KDvyn7DrH1Rv9J4ME1gAQAA9JEEVkUWlmzG6elpnJ+fx2QyiZOTk9jb24uIq8TS27dv48MPP4yXL1/eeM3FxUWcnp7Gmzdv4u7du3F8fFz4vpeXl/HkyZOIiDg5OYn9/f24vLyMp0+fxps3b66f9/Lly7h///6tEV6w60wtCLtDMpwx8fu1Pvc1cJPRG3TNsbYZ1wAA4yOBNQBjblQ6OzsrHCG1zOHhYbx+/Xqj993f31+7PBgCDU9AGQ0g9E3Zb9aYr3n7YjEemDqQthhxCQDAJiSwVthkYUm9tYBdINbB+GhIBGCXGb0B3TMDEgDzJLB6zA80ME/PdAAAhsSISwAAtiGBVSO9tYBdsCrWLYt5RmcBAABI5C7aZAYkAMZPAqtFGnWBqpZdqA95dGbVaQdNGwHV6SgDAAAAjJEE1hYWG1g1rgK7bFkiXgM7AABDZw1YmjDflvT9MaZ9CQAiJLBaoVEXWGaso63mVZ128JcFrzNtBFSj8QwAAAAYEwmsDSw2KGtcBQAAgPUMZdpo611Dt4YSKwConwRWT5iKAMZlk8T2Ll6AawyA9ul4Q9813UjluhuArj1LKeLe4+//HwAoJIFVI2tiwW7SGNycshs7jYwAAMNVdK/smhpYVHesWJYsdI8J0E8SWB2rui5M1dcB/SKRDfSV+ETfNN2gbdQvAF259Xs27aQXoX0HAJaRwKqBNbGACI3BdZq/iXs0janHOUvus9M0sgOriBMAcNuy+0S/nQD9JoHVoG2mFDQyC7onGd0v7330+Y2eihERv3z5sKPaAMCwPLr3+LpTCPRZ2TW4zmKMwew4nrXpaMMBgOUksICdJ1E1XLMbPsl9donjGVilKE5IXkFzrKsDANAMCawGlE0puM2IrBkjs2A7mySr9PbsVh1romhUYBdJzjMmRmTUx2dGX5Udm37PGKKy49Z0dfWpo41tnvY0gH4adAJr6Beydf/YLuOHmF30aDrdXNUexxp0hmtVcn8Z8REAgE1YVwcAoFmDTmANRdmIrEVVEnJlDezbNN7CLpKsGodVPfK3aVSQ2KIrdf92i3cMmREZwIyRmPTZreu3aWfK2Zq9399DOF63VbWNrSrtaQD9NooE1tAuWOv+sV3GlIMwvBjB9ori6uJx4EYFAICmuecev/c++vw6YcX4OIcBujWKBNbQVWlcr3MdrSo2acD14w20pcke+Zsm/td9P1il6rFjJAq7rM4RGdZLhH4yEpM+Wfe48/vRniqdGOsksQXQjtYTWF31bP/Xf/3X+MlPfjK6cje+aJ/rHTT/g77N91PXd7vpj31X33GXZe9aubtoyJ91F3UvukEp6wBQ9vyZuurfxY3NkI+biOHWf6j1LjKWfRnLfkSMZ1/Gsh8R7e+Lz66fxrQvEcPcn7JrvK73pc7R/rt239f1vf3/+A//ISK+X1O5ql++fKg9YyBlVyl303N41f3fLn7eQ7GL380u7rOyuyl7G4MagbVND6uhB4Impx288V4VLtBmdSmrw2yO501tkwj7Hwt/t9XzZejH11DK3UVD/qyHUPeVsfR//38i4irurRtPNr2x2cRi3Ybw2S8z1PqvW+91j4FVx2udPUyH+h0sGst+RIxnX5q6Np5Z1mmh7HWbrpfYZQJrWb0e3Xscj3o+MmUsx3PEuPYlYjz7cx0DptdyM32aWnydTk67dt9Xd7lrt5+s0S5y0+bH1642aNZR9rJOjE2UW0dia77tqs0RW2OJ8U0Z+rkwpHKVvXtlb6OxBNaqD6SssW3bxMy6X0TZ8+va3lTZ8z/QZc//7LPP4rPPPit9n00/61WvW7eXUpPqvHlflpj79xHx7Fe/qq2sqvpU7nHOtZ4fZb755pta3qdMXTGkruevY8h1b7I+s3g5//x1499G8bLFWHg71v1P8V9riH/bdkjY1KrYNn8N0acLsH//q19tFJPLfjfr+g1r4zOq8zqozjI20Ua9xvKdbPqaddVZRtV43qdG7jb19RjYtJyxxJk2X7OuLvdlk0brNqcdHFIcWfWdtD9dYz3XsnVZNqKv721H25Rd5/t3VXb5sbv+MfZ//W//69Kyy9o52+zYuMr8fVaVe6ohtF9sctzt4rkwhtjTZdlDOs6G9F1XfU3KW2b6U0r/GBH/b8FDP46Ib2wv3d7HOtnefdm7tP3Pc86/KHjuUiml/zsiXq3xkrI6jfH5farLrj2/T3XZtedXfe7aMWfJNc62ddnmNW2UsclrxlLGJq9po4xNXtNGGZu8po0yNnlNnWU0FW/Kyqvr+V7T33r1+TVtlNHn17RRxrLXbHpPtRhzyt5/6Nv7WKd1t/exTutu72Od1t3exzqtu33b96or3tRRlyrb2yhjDHUq265Ow61T2fah1ek65mydwAIAAAAAAIA6vdN1BZqQUjpY+PsopXSYUnqwbBubWfW5+qxhN4gFVLXp77RjiGWqHh+OHzaxGLcWHhObdoQ4Q1vEHBi+snO0jfucJWUfT/87m9t2NnusjrK5bfrdHs5/7nOPNXYPvEa5t44BvzX0SWsJrLYO/JTSYUS8mPv7ICIi53wREZOU0kHRthrKLfoRaONH6VYwarPRb/p535v+f1ufdaXA2sBnfTB9z6O2yp2WmVNKX0//e9pGuQXv2epnXZe+1mtRV/GjTl3Egjp0cV7XaYjn6Ka/0307hlZ9pn36zJepsB+34lMfVT0+5mNVX1X4Tm7Frb5a4zzpdYPJYtxaeKyJa91RxJcZcaafxJr+ajPmrDr/Su5Ttm7srlBuY42py8pO5ffftTTwpyWNydPHG7mOr1BuY0mNCmU38l0vK7fp73lWfhT8FrRxn7Ok7MOIuMg5P4uI/enfERHHKaWvI+Jym3Kr1m3xeynZVvu5sEbZtR6T08/5/vT7PZj/fue/q7qPjarlTt04Buo4Jrf5bBv6rovKrj3+lJTdSrvTGmXXfozPl7skxm68z60ksOoOxstMy5gPuh9HxGT6/5cRcViybWNFPwIt/ijdCEZtlLtE45/11MrA2tA+P8w5n8fVd9zWZ/1uzjnlnO9ExP2IOGtrf6fvcTk7p3pwfK2lr/Va1FX8aFhbsaAOXZzXtRjqObrF73RvjqFVn2nfPvMyFfaj7Ca3j3pzfGyj4rFzI261WsE1VDxPbsSwDqpZSUHcmlf3fcUo4suMONNPYk2/tRVzVp1/Sx7fqrG74nlfe2NqxbJv3X8X1WeLsgsbk6ePN3IdX6HcxpIaq8ouKqeNfY4Gv+cKurzP2Z9778vp3xERn+Sc70w/r8aUtGG20q5Ztezp02s9JnPOFznnk+mf+znnNyVPrfXYWKPciNvHwFbH5DafbRPfdVHZVetTU9mNtztVLbvu/S4pt/a27FYSWNHtxf5eRLyd+/tHJdu2UfQj0PiPUkkwau3HMKV0sPADtxfNf9YR1QJr3Y0JRxHxKiIi5/ykrc964fP9IOd82Ua5c2YXcq0fXzXoa70WdRI/6tRhLNhKV+d1zYZ8js7sRbXjpWhbV1Z9pn3/zGdW1bPsJreP9mLF8VEQq/po6XdSErf6qsp5sBjDhmgv6o1NY4kvM+JMP4k1w7UX9cWcVedfU43dVc77WhtTq5Zdcv9dVJ+1VWhMbuQ6vkK5jSU1Kjag1/5dryq3ye85YuVvwV40eJ+zrOyc87NpojIi4iAifjf9/1lH2kZHdhd9L221a65RdkRD8Wf6+Z7M/d1KO0qFciNuHwNblbvlZ1v7d11SdtX6bFV2W+1Oa5QdUeN+l51bc0+ppS27rQTWXvSn0al2JT8Ce9FS49tCMGqt3Ih4t6b3WVeVwFq0bRsfRsSPphniNsuNiOuM9hdtljsNbJcppT/OvXcrZddkL/pZrxu6jh816SoWbKvT83pbIzhHh2wvln+mqx7vi71YUs8lN7lDNYRYtRfLj52iuNVXe7H8+CqKYYwnvszshTjTR3sh1uy8VedfU43dFc/7WhtT1yx78f67qD4bW2xMnrMXzSY2CsttI6mxZJ+LytmLhvd57vGmvucufwtWlj0d5fBm1pg9bdi+iKu433jHmKLvpa12zQplRzQXf55ExElKaW+6qZXjpEq5TR0DG362RduaKLtqfbYtu9V2pwplRzSw3yXnVm1t2W0lsLo0ie9Pzr2I+K5k29YWfwTaUhCMGleSsZ9EC5912z+uc76bfbep/bng7+WcJ20WOD2eJhHxOCJ+k1Lqc4/Ywesqfmyry1hQky7P662M6BydRLXjpWgbLRhIfJrEkuNjQKMiqhhs3Jo38hjGmsSZXhJr+mkSNcecVedfU43dy8pt+p6/Qsy5cf9dZ326aL+pUm6TSY1lZTf5XVf4rGv/niv8FkyiofucNX6HDnPOp9PXHM/F+O+ihZHQRd9LW+dFlbLrPibT3FRpcTXK5LiNdpSq5ZYcAxuXO6/pz3bbspuqT8Fx1to1VZWym9jvknO4trbsH9TxJhVMorsbu+cR8cH0//cjYnaiFm3b1vWPQJTvc22fwywQzXqVRcRxG+VO7U8v/t+NiHendWn8s05XC8y9zVfzd64KrHXu83fx/Typk7jKYrdR7sz8PKBtlXscEY9zzpOU0mVEHLVYdh0m0c96lWk1ftSok1hQk67P620N/RydWed46csxNInln+mqx/tiEtXqOR+f+qrwOEop7U0vmm/Fqp42lE9i+XdSFLfOW6jXJiaxfF+KYtiT1mq3pbljqyyGbWoS44gvM5MQZ/poEmLNoDQYcyJWn383Grvj9v14reWuec9fa9lzru+/69rnkvab+WNxEg1cx1cod6b273lV2U1912vsc+3fc5T8FlQ4h+s4r1eVHSml42lD82x0xO/i+3h/JyKeblj2SkXfS0rpYnFbNHAuVC07pTSJ+o/Jw4iYXQ/sxdWUbm20o1Qtt+gY+N0W5ZZ93pNooT23atlNxJ+S2NNKu1PVslNK70aN+70i3tbWlv1O1Qpt6Xl8H/wbbXSaZhM/mMsqzrKMhxExyVdzMd7aVkO5iz8CRftc9+dwGDe/6MuWyo2c8/n0YJ+V3dZn/bv4vv53pn+3sc/nc++3F1eBv5XPuqCXYCvlzpt+15Muyt5CX+t1S0fxoxYdxoI6dHZe121I5+imv9M9O4YKP9O53ka9+syXWLUfRfGpl5YcH19NH78Vq3pq1XdSFLf6auXxNTMXw3ppMW5NzY6tumPTWOLLjDjTT2JNj7UZc4rOvxXnZ9H9eN3lVr3n30iFfV68/65ln6O4/WZVfK9jv1eV29j3XKHspr7rKvvcyPe85Leg9Byu67xeVfb0/c9SSl+nq2lVZ/X5+TTefN3wPVbR99JWu2bVsps4Jp/F9wmT6++phXaUdcq9cQzUcExu89k28V0Xld3Ed11UdlvtTlXLrnu/y+JtrW3ZKee8Rp02Nz1hLiNiP38/v+4oTE/oF3E1b+O7EXE/53xRtM91fg7TH9+fT/+8m6eLpjVdbtdmWfK42pcnc9sa3ee5cj9c6JnUdLn7EXGav18Ur7XvOF3NYXoZEe8uK6evx1df6zWvq/jBla7O67oM/RwdspLP+XXO+W7Z4320bD/K4lN3td0NFY+tG3Grryrsy60Yxnjiy4w4009iDUvuQ5aen0X343WWO31OpXv+hsouu//eeJ+n77EXxe03S+P7tvu9qtymvuc197nW77piuY18z5Qr+l7aatfcoOza488u2fazrfu7XlF2bd91hbIba3faoOxa9ntJubW2ZbeWwAIAAAAAAIAq2ppCEAAAAAAAACqRwAIAAAAAAKBXJLAAAAAAAADoFQksAAAAAAAAekUCCwAAAAAAgF75QdcVYLxSSgcRcRgRb3LOF13XBxgv8QZoi3gDtEnMAdoi3gBtEW9YhxFYNCKl9CAifhMRFxGxn1L6OqW0122tgDESb4C2iDdAm8QcoC3iDdAW8YZ1pZxz13VgZFJKRxHxIiLu5pzfpJRmB9ndnPObDqsGjIx4A7RFvAHaJOYAbRFvgLaIN2zCFII04WFExFzguYirIaECEVA38QZoi3gDtEnMAdoi3gBtEW9YmxFY1CqltB8RX0fERc75Xtf1AcZLvAHaIt4AbRJzgLaIN0BbxBs2ZQ0s6nYw/VfmHGiaeAO0RbwB2iTmAG0Rb4C2iDdsRAKLus0y6C87rQWwC8QboC3iDdAmMQdoi3gDtEW8YSMSWNTtg+m/v+u0FsAuEG+Atog3QJvEHKAt4g3QFvGGjVgDi62llM4i4igi3o2Ivenmy9njOec7HVQLGCHxBmiLeAO0ScwB2iLeAG0Rb6iDBBa1SSkdRMTriDjPOd/vuj7AeIk3QFvEG6BNYg7QFvEGaIt4wzZMIUidZkNBX3VaC2AXiDdAW8QboE1iDtAW8QZoi3jDxn7QdQUYlbvTf98UPZhS2o+rYaNv4mrY6Ls552ftVA0YmaXxZmYWd3LOT5qvEjBSq65vDiPiICJ+NP33Tc75tKW6AeNTJebci4jv4iru7EXEac550kblgFGpdE81k1J6YeQEsKFV1zdnEfE85/wmpbQXEYcR8aH7KiJMIUiNUkqv46rh5oeLN1DT4PMi53xvbttxROxpWAbWtSzezD1nPyJextUQdRc9wEZWXN8cRkTknC/mtr2IiP2c890AWNOKmHMQER/MdwKcNvgcijnAuqrcU8099yiu2nRSG3UDxqVCm/EfF17yLOd80k7t6DtTCFKng4iYlFz4HMdVQ/K16Y3Xw2mgAlhHabxJKe2llJ7G1YjPy8XHAda07Prm3nzyauqTiDiYNvQArGtZzDmJiLcL257HVcw5aLpiwOgsizeLPmy4LsC4rYo3JxHxw4i4k3NOklfMk8CiFtORDhERvyt5ysdRPEz0Mq6GhQJUsire5JwnOeeT6ejOSWsVA0anwvXN4dxzIuIqBsXV9c29wlcAlKgQcyKu7quKTOqtDTBmFePN7LlHEfG02RoBY1U13kzbcnRC5hYJLOoyS0KVzZ18EMUjITTwAOtaFW8A6rIq3uzH1TUOQB2WxpxpB53F9Wc+jqu19zT4AOuodE81bXieNF4bYMy04bCVH3RdAUZjNuf6q8UHKkwR+G7ttQHGrDTeANRsabzJOf9wcdv0ume2Bh/AOta6xplOG3gYEYtJLYBVqsabw5zzs8UR5wBrqBJv9hamYN+fzqoDEljU5oPpv0XZ9GUJqrdx1cgDUNWyeANQp03izXFEXOaczxuoDzBulWJOSukwrkZ/3ouI+0ZfARtYGW+mseaLdqoDjNjSeJNznqSUfjSfsEopHaWUXuaczdpFpJxz13VgBFJKOa4W4yvqiXwQEa/jaiG+y4XHXkTEnoAEVLUs3hQ890VcNSSfNl8zYGzWiTfT5+/H1TXPXQ3KwLrWjTnT17yIiFd6KQPrqBJvUkpHsw4502ucr3POqa06AuOwyfXN3Ovu5ZwvmqkZQ2ENLLZWYTG+yfTfvYLH9sJ8ykBF6yw2DLCNDePNi4j4meQVsK4trnE+iYizhWl3AEpViTcppWOjyYFtbdmGcxlXo83ZcRJY1GG2eHnZUNBZI07RVILvxtU0ggBVLI03ADVaK96klJ5GxCc5Z/EJ2MRG1zg550lcNfCc1F0hYLSWxptpg7MOg0AdVl7fpJQeTNcRXnQZlp0hrIFFPT6c/rtssfI3UT4CyyLnQFVV4g1AHSrHm5TSg4h4Op+8SintTRuWAapYGnPmpii9XzCVziSWrzsMMG/VNc5hRNxJKX08t20vIiKldBYRYYp2oKJV1zd7EXEWV+3Gi9c375a9jt1iDSy2llJ6HRH7K+ZOLrzISSn9MSL+QgMPUEWVeLPwfGtgARupGm+m03ZNFhuUU0oPrEkDVLUq5kwTWC/iKoG1uK5wjognrneAKta9p5q+5iAiXlsDC1hHxTbjW/dN08TWHyPijunZMYUga0spHS/MsX4QEc9WvOxxXPXiufE+EfFY8goos2G8AVjbJvEmpXQYV9NavE0pHcz9dxzW+ASWWDfmTBtvnsfC9OvTEaA66wCl3FMBbdkw3lxOk+TzzuKqc47kFUZgsZ5pEHoREZFzTtOGm5cR8cNViahpr8GTiHgV0ykucs4umoBCW8abs7hqVD6Kq0bki9C4A5TYJN7M9Qosc69gmi+Aba9xjuJqOp7vIuJHcTWt16lOgUCRbeLN3HucxVWH5IO4uq96aZQ5sKiG65vZeld34irOnDdYXQbEGlhsYhIRs0bg04g4qXLhM82aazwG1jGJzeKNWAOsaxJrxJvpY6bRATY1ic2ucc4jQoMOsI5JbBBvZtxbAWuYxObXN1DICCzWNu19sxdXmXE9b4DGiDdAW8QboE1iDtAW8QZoi3hDEySwAAAAAAAA6JV3uq4AAAAAAAAAzJPAAgAAAAAAoFf+f8kjqRUuwU2fAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1728x360 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\"\"\"\n",
    "\n",
    "true_sup_bkg = pythia_bkg_dict[\"sr\"][n_sts:,:-1]\n",
    "true_sup_sig = pythia_sig_dict[\"sr\"][n_sts:,:-1]\n",
    "\n",
    "feature_labels = [r\"$f_0$\", r\"$f_1$\", r\"$f_2$\", r\"$f_3$\", r\"$f_4$\", r\"$f_5$\"]\n",
    "\n",
    "\n",
    "fig = plot_feature_histograms(true_sup_bkg, true_sup_sig, 40, n_features, feature_labels)\n",
    "fig.show()\n",
    "\n",
    "np.save(os.path.join(STS_dir, \"true_sup_bkg.npy\"), true_sup_bkg)\n",
    "np.save(os.path.join(STS_dir, \"true_sup_sig.npy\"), true_sup_sig)\n",
    "\"\"\"\n",
    "\n",
    "#n_left = 55299\n",
    "n_left = 59781 # for 3-prong signal\n",
    "\n",
    "true_sup_bkg = pythia_bkg_dict[\"sr\"][:n_left,:-1]\n",
    "true_sup_sig = pythia_sig_dict[\"sr\"][n_sts:,:-1]\n",
    "\n",
    "feature_labels = [r\"$f_0$\", r\"$f_1$\", r\"$f_2$\", r\"$f_3$\", r\"$f_4$\", r\"$f_5$\"]\n",
    "\n",
    "\n",
    "fig = plot_feature_histograms(true_sup_bkg, true_sup_sig, 40, n_features, feature_labels)\n",
    "fig.show()\n",
    "\n",
    "np.save(os.path.join(STS_dir, \"true_sup_bkg.npy\"), true_sup_bkg)\n",
    "np.save(os.path.join(STS_dir, \"true_sup_sig.npy\"), true_sup_sig)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f7ec0b86",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "9af5df63",
   "metadata": {},
   "source": [
    "## Save out the nonSTS\n",
    "\n",
    "This is a little subtle because although we only pulled the STS from the SR, we need to save out only a portion of the SB to get the proper ratios of events in SB1/SR/SB2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3525fa95",
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"\n",
    "pythia_bkg_nonSTS = {}\n",
    "pythia_sig_nonSTS = {}\n",
    "\n",
    "frac_pythia_sig_dropped = n_sts / num_sig_pythia_events[\"sr\"]\n",
    "frac_pythia_bkg_dropped = n_sts / num_bkg_pythia_events[\"sr\"]\n",
    "\n",
    "\n",
    "for b in bands_to_analyze:\n",
    "    num_drop_bkg = int(frac_pythia_bkg_dropped*num_bkg_pythia_events[b])\n",
    "    pythia_bkg_nonSTS[b] = pythia_bkg_dict[b][num_drop_bkg:]\n",
    "    print(f\"Dropping {num_drop_bkg} bkg events from band {b}\")\n",
    "    \n",
    "    \n",
    "    num_drop_sig = int(frac_pythia_sig_dropped*num_sig_pythia_events[b])\n",
    "    pythia_sig_nonSTS[b] = pythia_sig_dict[b][num_drop_sig:]\n",
    "    print(f\"Dropping {num_drop_sig} sig events from band {b}\")\n",
    "    \n",
    "pickle.dump( herwig_bkg_dict, open( \"post_STS_wide/herwig_bkg_dict.p\", \"wb\" ) )\n",
    "#pickle.dump( pythia_bkg_nonSTS, open( \"post_STS_wide/pythia_bkg_nonSTS.p\", \"wb\" ) )\n",
    "#pickle.dump( pythia_sig_nonSTS, open( \"post_STS_wide/pythia_sig_nonSTS.p\", \"wb\" ) )\n",
    "\"\"\"\n",
    "\n",
    "pythia_sig_nonSTS = {}\n",
    "\n",
    "frac_pythia_sig_dropped = n_sts / num_sig_pythia_events[\"sr\"]\n",
    "\n",
    "for b in bands_to_analyze:\n",
    "    \n",
    "    num_drop_sig = int(frac_pythia_sig_dropped*num_sig_pythia_events[b])\n",
    "    pythia_sig_nonSTS[b] = pythia_sig_dict[b][num_drop_sig:]\n",
    "    print(f\"Dropping {num_drop_sig} sig events from band {b}\")\n",
    "    \n",
    "  \n",
    "pickle.dump( herwig_bkg_dict, open( \"post_STS_wide_3prong/herwig_bkg_dict.p\", \"wb\" ) )\n",
    "pickle.dump( pythia_bkg_dict, open( \"post_STS_wide_3prong/pythia_bkg.p\", \"wb\" ) )\n",
    "pickle.dump( pythia_sig_nonSTS, open( \"post_STS_wide_3prong/pythia_sig_nonSTS.p\", \"wb\" ) )\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0777452",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efd04207",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "aa934c55",
   "metadata": {},
   "outputs": [],
   "source": [
    "for b in bands_to_analyze:\n",
    "\n",
    "\n",
    "    loc_pulls = [bands_dict[b]]\n",
    "\n",
    "    loc_in_range_herwig_bkg = pull_from_mass_range(data_sim, loc_pulls)\n",
    "    loc_in_range_pythia_bkg = pull_from_mass_range(data_dat, loc_pulls)\n",
    "    loc_in_range_pythia_sig = pull_from_mass_range(data_signal, loc_pulls)\n",
    "\n",
    "\n",
    "   \n",
    "    # Plot some features\n",
    "    fig = plot_feature_histograms(loc_in_range_herwig_bkg, loc_in_range_pythia_bkg, 40, n_features, feature_labels, \n",
    "                                  third_plot = loc_in_range_pythia_sig)\n",
    "    fig.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9182b31b",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8a2447d9",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1cafbac1",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5a31dfad",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f130901f",
   "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.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}