How-to-rediscover-the-Higgs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to rediscover the Higgs boson yourself!\n",
    "This notebook uses ATLAS Open Data http://opendata.atlas.cern to show you the steps to rediscover the Higgs boson yourself!\n",
    "\n",
    "The idea is that you add extra cuts to increase the ratio of signal ($H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$) to background ($Z, t\\bar{t}, ZZ \\rightarrow \\ell\\ell\\ell\\ell$)\n",
    "\n",
    "First, try to reduce the amount of $Z$ and $t\\bar{t}$ background, since these are quite different to the signal.\n",
    "\n",
    "Then, try to reduce the amount of $ZZ \\rightarrow \\ell\\ell\\ell\\ell$, whilst keeping $H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$ signal\n",
    "\n",
    "The datasets used in this notebook have already been filtered to include at least 4 leptons per event, so that processing is quicker."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<CENTER><img src=\"HZZ_feynman.pdf\" style=\"width:40%\"></CENTER>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## First time setup\n",
    "This first cell only needs to be run the first time you open this notebook on your computer. \n",
    "\n",
    "If you close jupyter and re-open on the same computer, you won't need to run this first cell again.\n",
    "\n",
    "If you re-open on binder, you will need to run this cell again.\n",
    "\n",
    "If you run into a problem of \"uproot not being available\", Kernel -> Restart & Run All"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x10ecc4390>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pip/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x10ecc42b0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pip/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x10ecc4278>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pip/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x10ecc4048>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pip/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x10ecc4470>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pip/\u001b[0m\n",
      "Requirement already up-to-date: pip in /Users/meirinevans/.local/lib/python3.7/site-packages (19.2.1)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44400>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/numpy/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44710>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/numpy/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44908>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/numpy/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44a20>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/numpy/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44438>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/numpy/\u001b[0m\n",
      "Requirement already up-to-date: numpy in /Users/meirinevans/.local/lib/python3.7/site-packages (1.16.4)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a6bc50>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pandas/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a5ffd0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pandas/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a5fe80>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pandas/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a5fef0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pandas/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a44390>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/pandas/\u001b[0m\n",
      "Requirement already up-to-date: pandas in /Users/meirinevans/.local/lib/python3.7/site-packages (0.25.0)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78630>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/uproot/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78860>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/uproot/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78550>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/uproot/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a784a8>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/uproot/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78198>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/uproot/\u001b[0m\n",
      "Requirement already up-to-date: uproot in /Users/meirinevans/.local/lib/python3.7/site-packages (3.8.0)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a72748>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/matplotlib/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78438>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/matplotlib/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78160>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/matplotlib/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78d68>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/matplotlib/\u001b[0m\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a78f60>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/matplotlib/\u001b[0m\n",
      "Requirement already up-to-date: matplotlib in /Users/meirinevans/.local/lib/python3.7/site-packages (3.1.1)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a6b0f0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/keras/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a72978>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/keras/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a72c18>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/keras/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a72b00>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/keras/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a729b0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/keras/\u001b[0m\n",
      "Requirement already up-to-date: keras in /Users/meirinevans/.local/lib/python3.7/site-packages (2.2.4)\n",
      "\u001b[33mWARNING: Retrying (Retry(total=4, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a76b38>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/scikit-learn/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=3, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a769b0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/scikit-learn/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=2, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a76940>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/scikit-learn/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=1, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a76550>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/scikit-learn/\u001b[0m\n",
      "\u001b[33mWARNING: Retrying (Retry(total=0, connect=None, read=None, redirect=None, status=None)) after connection broken by 'NewConnectionError('<pip._vendor.urllib3.connection.VerifiedHTTPSConnection object at 0x104a767f0>: Failed to establish a new connection: [Errno 8] nodename nor servname provided, or not known')': /simple/scikit-learn/\u001b[0m\n",
      "Requirement already up-to-date: scikit-learn in /Users/meirinevans/.local/lib/python3.7/site-packages (0.21.2)\n",
      "Requirement already satisfied, skipping upgrade: pytz>=2017.2 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from pandas) (2018.9)\n",
      "Requirement already satisfied, skipping upgrade: python-dateutil>=2.6.1 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from pandas) (2.8.0)\n",
      "Requirement already satisfied, skipping upgrade: awkward>=0.12.0 in /Users/meirinevans/.local/lib/python3.7/site-packages (from uproot) (0.12.2)\n",
      "Requirement already satisfied, skipping upgrade: uproot-methods>=0.7.0 in /Users/meirinevans/.local/lib/python3.7/site-packages (from uproot) (0.7.0)\n",
      "Requirement already satisfied, skipping upgrade: cachetools in /Users/meirinevans/.local/lib/python3.7/site-packages (from uproot) (3.1.1)\n",
      "Requirement already satisfied, skipping upgrade: pyparsing!=2.0.4,!=2.1.2,!=2.1.6,>=2.0.1 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from matplotlib) (2.3.1)\n",
      "Requirement already satisfied, skipping upgrade: kiwisolver>=1.0.1 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from matplotlib) (1.0.1)\n",
      "Requirement already satisfied, skipping upgrade: cycler>=0.10 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from matplotlib) (0.10.0)\n",
      "Requirement already satisfied, skipping upgrade: h5py in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from keras) (2.9.0)\n",
      "Requirement already satisfied, skipping upgrade: pyyaml in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from keras) (5.1)\n",
      "Requirement already satisfied, skipping upgrade: keras-applications>=1.0.6 in /Users/meirinevans/.local/lib/python3.7/site-packages (from keras) (1.0.8)\n",
      "Requirement already satisfied, skipping upgrade: six>=1.9.0 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from keras) (1.12.0)\n",
      "Requirement already satisfied, skipping upgrade: scipy>=0.14 in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from keras) (1.2.1)\n",
      "Requirement already satisfied, skipping upgrade: keras-preprocessing>=1.0.5 in /Users/meirinevans/.local/lib/python3.7/site-packages (from keras) (1.1.0)\n",
      "Requirement already satisfied, skipping upgrade: joblib>=0.11 in /Users/meirinevans/.local/lib/python3.7/site-packages (from scikit-learn) (0.13.2)\n",
      "Requirement already satisfied, skipping upgrade: setuptools in /Users/meirinevans/anaconda3/lib/python3.7/site-packages (from kiwisolver>=1.0.1->matplotlib) (40.8.0)\n"
     ]
    }
   ],
   "source": [
    "import sys\n",
    "!{sys.executable} -m pip install --upgrade --user pip\n",
    "!{sys.executable} -m pip install -U numpy pandas uproot matplotlib keras scikit-learn --user"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## To setup everytime\n",
    "Cell -> Run All Below\n",
    "\n",
    "to be done every time you re-open this notebook"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot\n",
    "import pandas as pd\n",
    "import time\n",
    "import math\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import glob\n",
    "import random\n",
    "\n",
    "import infofile"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "lumi = 1000\n",
    "                                                                                                                                  \n",
    "tuple_path = \"Input/\"\n",
    "\n",
    "stack_order = ['data',r'$Z,t\\bar{t}$','ZZ',r'$H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "samples = {\n",
    "\n",
    "    'data': {\n",
    "        'list' : ['DataEgamma','DataMuons']\n",
    "    },\n",
    "\n",
    "    r'$Z,t\\bar{t}$' : {\n",
    "        'list' : ['Zee','Zmumu','ttbar_lep'],\n",
    "        'color' : \"#8700da\"\n",
    "    },\n",
    "\n",
    "    'ZZ' : {\n",
    "        'list' : ['ZZ'],\n",
    "        'color' : \"#f90000\"\n",
    "    },\n",
    "\n",
    "    r'$H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$' : {\n",
    "        'list' : ['ggH125_ZZ4lep','VBFH125_ZZ4lep'],\n",
    "        'color' : \"#4faeff\"\n",
    "    }\n",
    "\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_data_from_files(cutval=0):\n",
    "\n",
    "    data = {}\n",
    "    \n",
    "    for s in samples:\n",
    "        print(s+':')\n",
    "        frames = []\n",
    "        for val in samples[s]['list']:\n",
    "            prefix = \"MC/skim.mc_\"\n",
    "            if s == 'data':\n",
    "                prefix = \"Data/skim.\"\n",
    "            else: prefix += str(infofile.infos[val][\"DSID\"])+\".\"\n",
    "            fileString = tuple_path+prefix+val+\".root\"\n",
    "            print(fileString)\n",
    "            f = glob.glob(fileString,recursive=False)[0]\n",
    "            if f != \"\":\n",
    "                temp = read_file(f,val,cutval)\n",
    "                frames.append(temp)\n",
    "            else:\n",
    "                print(\"Error: \"+val+\" not found!\")\n",
    "        data[s] = pd.concat(frames)\n",
    "\n",
    "    return data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mllll_window(mllll):\n",
    "    return 120 < mllll < 130"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_weight(mcWeight,scaleFactor_PILEUP,scaleFactor_ELE,\n",
    "                scaleFactor_MUON, scaleFactor_TRIGGER):\n",
    "    return mcWeight*scaleFactor_PILEUP*scaleFactor_ELE*scaleFactor_MUON*scaleFactor_TRIGGER"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_xsec_weight(totalWeight,sample):\n",
    "    info = infofile.infos[sample]\n",
    "    weight = (lumi*info[\"xsec\"])/(info[\"sumw\"]*info[\"red_eff\"])\n",
    "    weight *= totalWeight\n",
    "    return weight"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def correlations(data, title, **kwds):\n",
    "    \"\"\"Calculate pairwise correlation between features.\n",
    "    \n",
    "    Extra arguments are passed on to DataFrame.corr()\n",
    "    \"\"\"\n",
    "    # simply call df.corr() to get a table of\n",
    "    # correlation values if you do not need\n",
    "    # the fancy plotting\n",
    "    corrmat = data.corr(**kwds)\n",
    "\n",
    "    fig, ax1 = plt.subplots(ncols=1, figsize=(6,5))\n",
    "    \n",
    "    opts = {'cmap': plt.get_cmap(\"RdBu\"),\n",
    "            'vmin': -1, 'vmax': +1}\n",
    "    heatmap1 = ax1.pcolor(corrmat, **opts)\n",
    "    plt.colorbar(heatmap1, ax=ax1)\n",
    "\n",
    "    ax1.set_title(title+\" Correlations\")\n",
    "\n",
    "    labels = corrmat.columns.values\n",
    "    for ax in (ax1,):\n",
    "        # shift location of ticks to center of the bins\n",
    "        ax.set_xticks(np.arange(len(labels))+0.5, minor=False)\n",
    "        ax.set_yticks(np.arange(len(labels))+0.5, minor=False)\n",
    "        ax.set_xticklabels(labels, minor=False, ha='right', rotation=70)\n",
    "        ax.set_yticklabels(labels, minor=False)\n",
    "        \n",
    "    plt.tight_layout()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_data(data):\n",
    "\n",
    "    bins = [80 + x*5 for x in range(35) ]\n",
    "    data_x = [82.5 + x*5 for x in range(34) ]\n",
    "\n",
    "    data_mllll = []\n",
    "    data_mllll_errors = []\n",
    "\n",
    "    mc_mllll = []\n",
    "    mc_weights = []\n",
    "    mc_colors = []\n",
    "    mc_labels = []\n",
    "    mc_in_mllll_window = [] # list for numbers of MC events with 120 < mllll < 130 GeV\n",
    "\n",
    "    for s in stack_order:\n",
    "        if s == \"data\":\n",
    "            data_mllll,_ = np.histogram(data[s].mllll.values, bins=bins)\n",
    "            data_mllll_errors = np.sqrt(data_mllll)\n",
    "        else:\n",
    "            mc_labels.append(s)\n",
    "            mc_mllll.append(data[s].mllll.values)\n",
    "            mc_colors.append(samples[s]['color'])\n",
    "            mc_weights.append(data[s].totalWeight.values)\n",
    "            mc_in_mllll_window.append([data[s].totalWeight.values[mllll_iter] for mllll_iter in range(len(data[s].mllll.values)) if 120 < data[s].mllll.values[mllll_iter] < 130])\n",
    "    \n",
    "    HZZ_in_mllll_window = sum(mc_in_mllll_window[2]) # number signal MC events with 120 < mllll < 130 GeV\n",
    "    background_in_mllll_window = sum(mc_in_mllll_window[0]+mc_in_mllll_window[1]) # number background MC events with 120 < mllll < 130 GeV\n",
    "    SoversqrtB = HZZ_in_mllll_window/math.sqrt(background_in_mllll_window) # calculate significance\n",
    "    print('Signal/sqrt(Background) for 120<mllll<130 '+str(SoversqrtB))\n",
    "    \n",
    "    top = np.amax(data_mllll)+math.sqrt(np.amax(data_mllll))\n",
    "\n",
    "    plt.figure()\n",
    "    plt.hist(mc_mllll,bins=bins,weights=mc_weights,stacked=True,color=mc_colors, label=mc_labels)\n",
    "    plt.errorbar( x=data_x, y=data_mllll, yerr=data_mllll_errors, fmt='ko', label='Data')\n",
    "\n",
    "    #X = np.arange(115,135,5) # gives list [115,120,125,130]\n",
    "    #data_in_window = data_mllll[7:11] # gives list of data y value for [115,120,125,130]\n",
    "    #x = np.sum(X*data_in_window)/np.sum(data_in_window) # Gaussian mean\n",
    "    #width = np.sqrt(np.abs(np.sum((X-x)**2*data_in_window)/np.sum(data_in_window))) # Gaussian width\n",
    "    #fit = lambda t : np.amax(data_in_window)*np.exp(-(t-x)**2/(2*width**2)) # Gaussian fit\n",
    "    #plt.plot(X, fit(X), '-') # plot Gaussian\n",
    "    \n",
    "    plt.xlabel(r'$M_{\\ell\\ell\\ell\\ell}$ [GeV]',fontname='sans-serif',horizontalalignment='right',x=1.0,fontsize=11)\n",
    "\n",
    "    plt.ylabel(r'Events',fontname='sans-serif',horizontalalignment='right',y=1.0,fontsize=11)\n",
    "    #plt.yscale('log')                                                                                                                                                                        \n",
    "    plt.ylim(bottom=0,top=top)\n",
    "\n",
    "    ax = plt.gca()\n",
    "    plt.text(0.05,0.97,r'$\\mathbf{{ATLAS}}$ Open Data',ha=\"left\",va=\"top\",family='sans-serif',transform=ax.transAxes,fontsize=13)\n",
    "    plt.text(0.05,0.92,'for education only',ha=\"left\",va=\"top\",family='sans-serif',transform=ax.transAxes,style='italic',fontsize=8)\n",
    "    plt.text(0.05,0.9,r'$\\sqrt{s}=8\\,\\mathrm{TeV},\\;\\int L\\,dt=1\\,\\mathrm{fb}^{-1}$',ha=\"left\",va=\"top\",family='sans-serif',transform=ax.transAxes)\n",
    "\n",
    "    plt.legend()\n",
    "\n",
    "    plt.savefig(\"plot.pdf\")\n",
    "    \n",
    "    # scatter plot of signal and background lep_n vs mllll\n",
    "    #plt.figure()\n",
    "    #for s in stack_order:\n",
    "    #    if s != 'data': plt.scatter(data[s].lep_n,data[s].mllll,color=samples[s]['color'],label=s)\n",
    "    #plt.xlabel(r'Leptons',fontname='sans-serif',horizontalalignment='right',y=1.0,fontsize=11)\n",
    "    #plt.ylabel(r'$M_{\\ell\\ell\\ell\\ell}$ [GeV]',fontname='sans-serif',horizontalalignment='right',x=1.0,fontsize=11)\n",
    "    #plt.legend()\n",
    "    \n",
    "    # scatter plot of signal and background lep_pt[0] vs mllll\n",
    "    #plt.figure()\n",
    "    #for s in stack_order:\n",
    "    #    if s != 'data': plt.scatter(data[s].lep_pt.apply(lambda x: x[0]),data[s].mllll,color=samples[s]['color'],label=s)\n",
    "    #plt.legend()\n",
    "    \n",
    "    # boxplot of signal and background mllll\n",
    "    #plt.figure()\n",
    "    #for s in stack_order:\n",
    "    #    if s != 'data': data[s]['Process'] = s\n",
    "    #data_all = pd.concat([data[r'$Z,t\\bar{t}$'],data['ZZ'],data[r'$H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$']])\n",
    "    #data_all.boxplot(by='Process',column=[\"mllll\"])\n",
    "    \n",
    "    # plot correlation matrices for signal and background\n",
    "    #correlations(pd.concat([data[r'$Z,t\\bar{t}$'],data['ZZ']]),\"background\")\n",
    "    #correlations(data[r'$H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$'], \"signal\")\n",
    "\n",
    "    return"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_mllll(lep_pts,lep_etas,lep_phis):\n",
    "    theta_0 = 2*math.atan(math.exp(-lep_etas[0]))\n",
    "    theta_1 = 2*math.atan(math.exp(-lep_etas[1]))\n",
    "    theta_2 = 2*math.atan(math.exp(-lep_etas[2]))\n",
    "    theta_3 = 2*math.atan(math.exp(-lep_etas[3]))\n",
    "    p_0 = lep_pts[0]/math.sin(theta_0)\n",
    "    p_1 = lep_pts[1]/math.sin(theta_1)\n",
    "    p_2 = lep_pts[2]/math.sin(theta_2)\n",
    "    p_3 = lep_pts[3]/math.sin(theta_3)\n",
    "    pz_0 = p_0*math.cos(theta_0)\n",
    "    pz_1 = p_1*math.cos(theta_1)\n",
    "    pz_2 = p_2*math.cos(theta_2)\n",
    "    pz_3 = p_3*math.cos(theta_3)\n",
    "    px_0 = p_0*math.sin(theta_0)*math.cos(lep_phis[0])\n",
    "    px_1 = p_1*math.sin(theta_1)*math.cos(lep_phis[1])\n",
    "    px_2 = p_2*math.sin(theta_2)*math.cos(lep_phis[2])\n",
    "    px_3 = p_3*math.sin(theta_3)*math.cos(lep_phis[3])\n",
    "    py_0 = p_0*math.sin(theta_0)*math.sin(lep_phis[0])\n",
    "    py_1 = p_1*math.sin(theta_1)*math.sin(lep_phis[1])\n",
    "    py_2 = p_2*math.sin(theta_2)*math.sin(lep_phis[2])\n",
    "    py_3 = p_3*math.sin(theta_3)*math.sin(lep_phis[3])\n",
    "    sumpz = pz_0 + pz_1 + pz_2 + pz_3\n",
    "    sumpx = px_0 + px_1 + px_2 + px_3\n",
    "    sumpy = py_0 + py_1 + py_2 + py_3\n",
    "    sumE = p_0 + p_1 + p_2 + p_3\n",
    "    mllll = sumE**2 - sumpz**2 - sumpx**2 - sumpy**2\n",
    "    return math.sqrt(mllll)/1000."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_mll(lep_pts,lep_etas,lep_phis,lep_charges,lep_types):\n",
    "    # this is only pseudo-code to tell you what to do!\n",
    "    # you need to decide how to find i & j yourself\n",
    "    # same-type-opposite-flavour\n",
    "    mll = 2*lep_pts[i]*lep_pts[j!=i]\n",
    "    cosh = math.cosh(lep_etas[i]-lep_etas[j!=i])\n",
    "    cos = math.cos(lep_phis[i]-lep_phis[j!=i])\n",
    "    mll *= ( cosh - cos )\n",
    "    # Z mass 90 GeV\n",
    "    return math.sqrt(mll)/1000."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncommenting a new cut\n",
    "If you add a cut: Cell -> Run All Below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def read_file(path,sample,cutval):\n",
    "    start = time.time()\n",
    "    print(\"\\tProcessing: \"+sample)\n",
    "    mc = uproot.open(path)[\"mini\"]\n",
    "    data = mc.pandas.df([\"lep_n\",\"lep_pt\",\"lep_eta\",\"lep_phi\",\"lep_charge\",\"lep_type\",\"lep_etcone20\",\"lep_trackd0pvunbiased\",\"lep_tracksigd0pvunbiased\",\n",
    "                         \"mcWeight\",\"scaleFactor_PILEUP\",\"scaleFactor_ELE\",\"scaleFactor_MUON\", # add more variables here if you make cuts on them\n",
    "                         \"scaleFactor_TRIGGER\"], flatten=False)\n",
    "\n",
    "    nIn = len(data.index)\n",
    "\n",
    "    if 'Data' not in sample:\n",
    "        data['totalWeight'] = np.vectorize(calc_weight)(data.mcWeight,data.scaleFactor_PILEUP,data.scaleFactor_ELE,data.scaleFactor_MUON,data.scaleFactor_TRIGGER)\n",
    "        data['totalWeight'] = np.vectorize(get_xsec_weight)(data.totalWeight,sample)\n",
    "\n",
    "    data.drop([\"mcWeight\",\"scaleFactor_PILEUP\",\"scaleFactor_ELE\",\"scaleFactor_MUON\",\"scaleFactor_TRIGGER\"], axis=1, inplace=True)\n",
    "    \n",
    "    # cut on minimum lepton pt\n",
    "    \n",
    "    # cut on lepton etcone20\n",
    "    \n",
    "    # cut on lepton d0\n",
    "    \n",
    "    # example of adding column that takes the return of the function cut_lep_pt_min\n",
    "    #data['lep_pt_min'] = data.apply(cut_lep_pt_min,axis=1)\n",
    "    \n",
    "    # example of cut on minimum number of leptons passing baseline requirements\n",
    "    #fail = data[ np.vectorize(cut_n_lep_min)(data.lep_pt_min) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # cut on number of leptons\n",
    "    fail = data[ np.vectorize(cut_n_lep)(data.lep_n) ].index\n",
    "    data.drop(fail, inplace=True)\n",
    "    \n",
    "    # dataframe contents can be printed at any stage like this\n",
    "    #print(data)\n",
    "\n",
    "    # cut on lepton charge\n",
    "    #fail = data[ np.vectorize(cut_lep_charge)(data.lep_charge) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # cut on lepton type\n",
    "    #fail = data[ np.vectorize(cut_lep_type)(data.lep_type) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # cut on lepton pt\n",
    "    #fail = data[ np.vectorize(cut_lep_pt)(data.lep_pt) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # cut on deltaR\n",
    "    #fail = data[ np.vectorize(cut_deltaR)(data.lep_eta,data.lep_phi...\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # cut on minimum opposite-charge-same-type lepton pair invariant mass\n",
    "    #fail = data[ np.vectorize(cut_OCST)(data....\n",
    "\n",
    "    # calculation of Z boson candidate 1 invariant mass\n",
    "    #data['mZ1'] = np.vectorize(calc_mZ1)(data.lep_pt,data.lep_eta,data.lep_phi)\n",
    "    \n",
    "    # cut on mZ1\n",
    "    #fail = data[ np.vectorize(cut_mZ1)(data.mZ1) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # calculation of Z boson candidate 2 invariant mass\n",
    "    #data['mZ2'] = np.vectorize(calc_mZ2)(data....\n",
    "    \n",
    "    # cut on mZ2\n",
    "    #fail = data[ np.vectorize(cut_mZ2)(data.mZ2) ].index\n",
    "    #data.drop(fail, inplace=True)\n",
    "    \n",
    "    # calculation of 4-lepton invariant mass\n",
    "    data['mllll'] = np.vectorize(calc_mllll)(data.lep_pt,data.lep_eta,data.lep_phi)\n",
    "\n",
    "    mllll_window_list = data[ np.vectorize(mllll_window)(data.mllll) ] # return events with 120 < mllll < 130 GeV\n",
    "    \n",
    "    # example of expanding lep_pt list column into individual columns whilst requiring exactly 4 leptons\n",
    "    # need to change cut_lep_n to require exactly 4 leptons\n",
    "    #data[['lep1_pt','lep2_pt','lep3_pt','lep4_pt']] = pd.DataFrame(data.lep_pt.values.tolist(), index= data.index)\n",
    "\n",
    "    # example of expanding lep_pt list column into individual columns without requiring exactly 4 leptons\n",
    "    # need to do this for columns that you wish to use for fit_BDT\n",
    "    #if max(data.lep_n) < 5: \n",
    "    #    df_split = pd.DataFrame(data['lep_pt'].values.tolist(), columns=['lep1_pt','lep2_pt','lep3_pt','lep4_pt'], index=data.index)\n",
    "    #    df_split['lep5_pt'] = 0\n",
    "    #    df_split['lep6_pt'] = 0\n",
    "    #elif max(data.lep_n) < 6: \n",
    "    #    df_split = pd.DataFrame(data['lep_pt'].values.tolist(), columns=['lep1_pt','lep2_pt','lep3_pt','lep4_pt','lep5_pt'], index=data.index)\n",
    "    #    df_split['lep6_pt'] = 0\n",
    "    #else: df_split = pd.DataFrame(data['lep_pt'].values.tolist(), columns=['lep1_pt','lep2_pt','lep3_pt','lep4_pt','lep5_pt','lep6_pt'], index=data.index)\n",
    "    #df_split.fillna(0, inplace=True)\n",
    "    #data = pd.concat([data, df_split], axis=1)\n",
    "    \n",
    "    #print(data[['lep_eta','lep_phi']])\n",
    "\n",
    "    nOut = len(data.index)\n",
    "\n",
    "    elapsed = time.time() - start\n",
    "    print(\"\\t\\tTime taken: \"+str(elapsed)+\", nIn: \"+str(nIn)+\", nOut: \"+str(nOut))\n",
    "    \n",
    "    return data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Changing an already uncommented cut"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you change a cut: Cell -> Run All Below\n",
    "\n",
    "If you uncomment a cut here, you also need to uncomment the corresponding cut in the cell above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# cut on number of leptons\n",
    "def cut_n_lep(lep_n):\n",
    "    # return when number of leptons is less than 4\n",
    "    return lep_n < 4\n",
    "\n",
    "# cut on lepton charge\n",
    "def cut_lep_charge(lep_charge):\n",
    "    # return when sum of lepton charges is not equal to 0\n",
    "    # exclamation mark (!) means \"not\"\n",
    "    # so != means \"not equal to\"\n",
    "    # first lepton is [0], 2nd lepton is [1] etc\n",
    "    return lep_charge[0] + lep_charge[1] + lep_charge[2] + lep_charge[3] != 0\n",
    "\n",
    "# cut on lepton type\n",
    "def cut_lep_type(lep_type):\n",
    "# for an electron lep_type is 11\n",
    "# for a muon lep_type is 13\n",
    "    sum_lep_type = lep_type[0] + lep_type[1] + lep_type[2] + lep_type[3]\n",
    "    return (sum_lep_type != 44) and (sum_lep_type != 48) and (sum_lep_type != 52)\n",
    "\n",
    "# cut on lepton pt\n",
    "#def cut_lep_pt(lep_pt):\n",
    "# want to throw away events where the 2nd highest pt lepton used has lep_pt[1] < 15000 \n",
    "# or\n",
    "# 3rd highest pt lepton used has lep_pt[2] < 10000\n",
    "\n",
    "# cut on minimum opposite-charge-same-type lepton pair invariant mass\n",
    "#def cut_mOCST():\n",
    "# want to throw away events if the invariant mass of any opposite-charge-same-type lepton pair is < 5\n",
    "\n",
    "# cut on invariant mass of Z boson candidate 1\n",
    "#def cut_mZ1(mZ1):\n",
    "# want invariant mass of same-type-opposite-charge lepton pair that's closest to Z mass (91 GeV) to be in range 50 < m < 106 GeV\n",
    "\n",
    "# cut on invariant mass of Z boson candidate 2\n",
    "#def cut_mZ2(mZ2):\n",
    "# want invariant mass of remaining lepton pair that's closest to Z mass (91 GeV) to be in range 17.5 < m < 115 GeV\n",
    "# advanced: vary the lower range monotically from 17.5 at mllll=120 to 50 at mllll=190, and constant above mllll=190\n",
    "\n",
    "# cut on deltaR\n",
    "# want to throw away leptons that are separated from all other leptons by deltaR = math.sqrt(delta(lep_eta)**2 + delta(lep_phi)**2) < 0.2\n",
    "# want to throw away leptons that are separated from other leptons of the same type by deltaR = math.sqrt(delta(lep_eta)**2 + delta(lep_phi)**2) < 0.1\n",
    "\n",
    "# example of returning list where every element passes minimum lep_pt requirement\n",
    "#def cut_lep_pt_min(data):\n",
    "#    return [data.lep_pt[i] for i in range(len(data.lep_pt)) if data.lep_pt[i] > 6000]\n",
    "\n",
    "# cut on minimum lepton pt\n",
    "# want to throw away muons with lep_pt < 6000\n",
    "# want to throw away electrons with lep_pt < 7000\n",
    "\n",
    "# cut on maximum lepton etcone20\n",
    "# want to throw away muons with lep_etcone20/lep_pt < 0.3\n",
    "# want to throw away electrons with lep_etcone20/lep_pt < 0.2\n",
    "\n",
    "# cut on maximum lepton d0\n",
    "# want to throw away muons with lep_trackd0pvunbiased/lep_tracksigd0pvunbiased < 3.5\n",
    "# want to throw away electrons with lep_trackd0pvunbiased/lep_tracksigd0pvunbiased < 6.5\n",
    "\n",
    "# example of cutting on length of list passing minimum requirements\n",
    "#def cut_n_lep_min(lep_pt_min):\n",
    "#    return len(lep_pt_min) < 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data:\n",
      "Input/Data/skim.DataEgamma.root\n",
      "\tProcessing: DataEgamma\n",
      "\t\tTime taken: 0.19355392456054688, nIn: 89, nOut: 89\n",
      "Input/Data/skim.DataMuons.root\n",
      "\tProcessing: DataMuons\n",
      "\t\tTime taken: 0.2074899673461914, nIn: 112, nOut: 112\n",
      "$Z,t\\bar{t}$:\n",
      "Input/MC/skim.mc_147770.Zee.root\n",
      "\tProcessing: Zee\n",
      "\t\tTime taken: 0.41375064849853516, nIn: 720, nOut: 720\n",
      "Input/MC/skim.mc_147771.Zmumu.root\n",
      "\tProcessing: Zmumu\n",
      "\t\tTime taken: 0.2219090461730957, nIn: 1057, nOut: 1057\n",
      "Input/MC/skim.mc_117050.ttbar_lep.root\n",
      "\tProcessing: ttbar_lep\n",
      "\t\tTime taken: 0.22835612297058105, nIn: 897, nOut: 897\n",
      "ZZ:\n",
      "Input/MC/skim.mc_105986.ZZ.root\n",
      "\tProcessing: ZZ\n",
      "\t\tTime taken: 0.307448148727417, nIn: 1973, nOut: 1973\n",
      "$H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$:\n",
      "Input/MC/skim.mc_160155.ggH125_ZZ4lep.root\n",
      "\tProcessing: ggH125_ZZ4lep\n",
      "\t\tTime taken: 2.3683791160583496, nIn: 21783, nOut: 21783\n",
      "Input/MC/skim.mc_160205.VBFH125_ZZ4lep.root\n",
      "\tProcessing: VBFH125_ZZ4lep\n",
      "\t\tTime taken: 2.1469781398773193, nIn: 20954, nOut: 20954\n",
      "Signal/sqrt(Background) for 120<mllll<130 0.14425894026493355\n",
      "Time taken: 8.500607967376709\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhTZdr48e9DZQepIshmW5iRAgVEYARZbFHZBEEQESwIiFZGeQdEfy6DCAh11FFQXhy1KMpoBRkRQR0X9AXUQWRYKpStaGnLLmWzWCilvX9/JI1Nm7RJm7W5P9eVqzlPznLnEO6cPOdZjIiglFIqdFTzdwBKKaV8SxO/UkqFGE38SikVYjTxK6VUiNHEr5RSIUYTv1JKhZhL/B2AK6644gqJiory2fH27t0LQHR0tM+OqZRSnrZly5ZsEWlUsjwoEn9UVBSbN2/22fHi4uIAWLdunc+OqZRSnmaMyXRUrlU9SikVYjTxK6VUiHEp8RtjphljGlifv2OM2WOM6efd0JRSSnmDq3X840VknjGmD9AYuAdYAHzptciUUgEtPz+fgwcPcv78eX+HEvJq1apFixYtqF69ukvru5r4C6x/+wDJIrLBGKPVREqFsIMHD1K/fn2ioqIwxvg7nJAlIpw4cYKDBw/SsmVLl7ZxNXmfM8Y8BowGvjSWf+UaFYxTKVUFnD9/noYNG2rS9zNjDA0bNnTrl5eriX880BR4TESOAq2AZLcjVEpVKZr0A4O7/w6uJv7rRGSqiHwIICI/AwfcjC2g3HPPPRhjMMYwfPhwu9c2btzI+vXrba+XfERFRREXF4cxhvHjx1fqWGD5yTxu3DhatGhBjRo1uOKKK+jSpQuJiYll7vf7779nwIABNGjQgLp169KlSxeWLFni1nnwhri4OGrWrEn9+vVp0KABrVq1YuzYsWzZssXt/cydO9dLUSpl0bt3bzp16kTXrl39HYrPuJr4p7lYZmOMucoYs9YYs8sYs9MYM8VaPssYc8gYk2J93OJu0JWVk5PD8uXLbcuffPIJx44dsy03adKE5s2bM2XKFJo3bw5A27ZtmTJlClOmTOGee+7x2LEAhg0bxj//+U/q1avHxIkT6devH3l5eXz11VdO9/vll1/Sp08frr/+etLT0zl+/DiPPfYYU6dOZebMmS7H5y0zZswgJyeHM2fOsHbtWiIjI+nevTsrV670d2hK2fn2229JSUnxaSdRvxMRpw+gK/AgkAU8UOzxBJBazrZNgc7W5/WBNKAdMAt4pKxtSz66dOkinvT6668LINahIASQ5557zvZ6bGysxMbG2p4DMm7cOLt9OCt391gnT54UQADZsmWL3bbZ2dlO9/vHP/5Rxo8fX6r8rbfekrCwMNm/f7+tLDIyUmbPni09e/aUunXrSpcuXWTTpk2213/77Td5+OGHJSoqSi677DLp37+/7Nu3z277xMREufHGG6Vu3boSExMj//nPf5zGFhsbK3PmzClVfs8990jz5s2lsLBQREReeukliY6Olnr16slVV10ljz/+uFy8eFFERB588EGpVq2a1KhRQ+rWrSutW7cudxvlW7t27fJ3CKoYR/8ewGZxkFPLu+Jvbk3+dYE/FXs0wVLvX9YXyhER2Wp9ngPstu7P79544w0ARowYwZ133gnAm2++6Zdj1a9fn0svvRSAwYMHM27cOF555RXS0tJo2LChw32mpaXx008/MWbMmFKv3XXXXYgIa9assSt/7bXXePnllzl58iQjRozglltu4ddffwXgvvvuY8+ePWzcuJGjR4/SrVs3Bg8eTH5+vm37xYsXs2DBAs6cOUPfvn0ZN26c2+di1KhRHDp0yDYWUosWLfjss8/49ddfWbVqFYsXL7adr4ULF9K7d29mzJjB2bNnXdpGhaYJEybQqVMn26NJkyZcfvnlTtc/ePAg77//vtPlkODo26DkA+jnynplbB+F5VfDpViu+DOA7cBi4DIn2yQAm4HNERERHvtW3L59u+0Ke926dbJt2zbb8vr160XEc1f8rhxLRGT16tXSuHFj22tFj3Hjxtmujov77rvvBJDdu3c7PG7jxo1l7ty5tuXIyEh58sknbcuFhYVy1VVXSXJyshw/flwAyczMtL1eUFAgl156qXz77be27Z9//nnb66mpqQLI6dOnHR7f2RX/rl27BJDvvvvO4XYPP/yw3HHHHeXup6xtlO8E4hV/enq6NG/eXFatWuV0nbffflseffRRp8vByp0rfpfa8YvIl8aYm4A/UKztv4j8o7xtjTH1gBXAVBH51RjzKjDHmtzmAC9i6RBW8phJQBJA165dPTYj/KJFiwBo1qwZvXv3plq1arRp04Y9e/bwxhtvcMMNN3jqUC4f69Zbb+Xw4cP88MMPfP/99yxfvpxNmzaxZMkSxo0bR58+fez226iRZbC9Q4cO0aZNG7vXLly4QHZ2tm2dIsVHNzXGEBERwcGDB9m/fz8AHTt2tFs/Pz+fAwd+v3/ftGlT2/O6desClvsXDRo0cPl8HDx4EMD2S2bp0qXMmzeP9PR0Ll68yIULF+jevXuZ+6jINso3njX7PLq/x+Vqt9bPzs5mwIABzJgxgyFDhjhc57vvvmPatGmEh4fzxRdfsGDBArvlDz/8kFatWnki/IDm6pANbwMvA734vbqn3FvgxpjqWJJ+svzeIuiYiBSISCGwCLiuYqG7Ly8vj+RkSyvUw4cPExYWhjGGPXv2APDBBx9w5swZnx4rPz+ftWvXUq1aNXr06MHDDz/M559/bttPTk5OqX1fffXVtGrVivfee6/Ua8uWLcMYQ9++fe3KMzIybM9FhKysLFq0aEFkZCQA+/bt4/Tp07ZHbm4uo0ePrvR5KO7999+nefPmREdHc+DAAcaMGcOTTz7JkSNHOHPmDA8++GDRrz0AqlWz/3i6so0KTbm5udx6662MHDmS+++/3+l6vXr14k9/+hOrVq0iJSWFG264wW45FJI+uN5ztwcQIyL55a5pZe3k9SawW0TmFStvKiJHrIvDgFRX91lZK1as4OTJkwAMGjSISy75/e1//PHHnDt3zpasXfXpp5+WuuLcuHGjy8e6++67ufHGG4mKiqJLly40btyYjRs3AnDZZZdx/fXXlzqmMYaFCxdy22230bJlSx544AFq167Np59+ytSpU3nsscdK9eBbvHgxw4YNo0OHDsyfP5/c3FwGDRpEgwYNuOuuu3jggQd46aWXaN68OadPn2bt2rX07duXevXquXU+HDlw4ABvvPEGb7/9Nu+//z7GGM6ePUthYSGNGjWievXqbNy4kXfeeYe2bdvatmvSpAk//fSTbdmVbVToKSgoYNSoUbRp04Y5c+aUu/7evXvtfimXXA4Frib+irTZ7wmMBXYYY1KsZX8FRhtjOmGp6skAnH89e1jRTcC4uDg++eQTu9eGDx/OypUreeONN2w3W12RnZ1NdnZ2hY+VkJDAI488wjfffMP69es5c+YMl19+OYMGDWLmzJmlqmyKDBw4kK+//pqnn36a559/noKCAlq3bs0LL7zgsLlpQkICf/nLX0hJSSE6OppPP/3UVk2zaNEinnnmGeLi4jh69Cjh4eH07t2bfv0qPg7fnDlzeO6552y9Cnv06MGGDRu47jrLD7y2bdsye/Zshg4dyoULF+jTpw+jR48mJSXFto+HHnqICRMmEB4eTvPmzdm5c2e526jQ88ADD5Cfn2+rWi1LdnY2DRo0sF2IlVwOFcaVn8nWevkY4CPA1i/YlTp+T+jatat4oo3tTz/9xNixY+nSpQsLFy50up47E7F89NFH/PLLLyQkJFQ4rqVLlzJy5EjCwsJ4/vnnmThxotMWPRURFRXF3Llzba2ABgwYYFedpFRF7N692+7Xlj/q+GfPns3HH3/MunXrHP46vemmm/jnP/9p64+zefNmnnrqKf797387XA5mJf89AIwxW0SkVLW8qx24agE/Ax1wo44/0LzzzjvMnz+/zKQPuFVnvHv3bmJiYioV1+LFiwkLCwPg0Ucf9WjSL+m3337zSPWNUv6WkZHBrFmzOHHiBL169bI15yxqNl1YWMhPP/1k17SzTZs2ZGdn0759ezZs2FBqOVS4lPhFZIKDh+vdVwPA//3f/7Fw4UIef/xxfvzxR6ZOnUqPHj0YMGAAv/32G1u3bqV///4MGzaMEydO2LY7dOgQQ4cOpXfv3syaNQuAbdu20atXL/r3789//vMfYmJiePTRR9m+fTtg+emZnp5Obm4ud999N3369GHEiBEATJs2zdZF/PDhwyxcuJBt27YRFxfHnj17GDhwIACff/453bt3p1u3bqxduxaA2NhYpkyZQrt27Urdizh//jxjxoyhZ8+ejBo1isLCQpYvX058fDzHjh1j5syZXLhwgV27dtG2bVtWrFjByy+/DGA7H0oFk6ioKESE/fv3k5KSYnsUtcnftWsXt99+O7Vr17ZtU69ePTZt2kRqaio9evQotRwyHLXxLPkA6mBpeplsXW4D3ObKtp54eKrnbv/+/UXE0m732WefFRGRl19+WZYuXSpLliyRe++9V0Ts2/GPHDlS0tLSRETk1ltvlby8POnfv79kZ2dLTk6OtGzZUkREBg8eLOfPnxcRkQEDBkhhYaFMmzZNPvnkExERuXDhgoiInD17VkREFi9eLO+99578/PPP8uCDD4qISGZmptx3331y5swZiYuLk9zcXDl9+rQMGDBA8vLyJDIyUk6dOiXp6ely99132723WbNmybJly0RE5KGHHpLvv/9ennrqKVv7+7vvvlvS09PlrbfekmXLlsn+/ftlwoQJIiJy5513ypEjRzxyjlXoCMR2/KHMkz13i7wKVAc6WZcPAv4fEMYNeXl51KpVC7C0qvnzn/8MWJoM1q5dm9TU1FL19BcuXGD9+vXcd999xMXFsX//frKysmjdujUNGzbkkksusbUGyM3NpWbNmogIhYWFGGPYuXMngwYNAqB69eocOnSIiRMn0qdPH5555hmaNWvG9u3bbW3oU1NTad++Pd988w3Dhw+ndu3ahIWFUatWLdLS0hg0aBDh4eFkZ2eXarWzYcMG20/cove0c+dO2yByZ86coUWLFqSmphITE0NUVBQHDhzghx9+ICIigiZNmnjnxCulAo6rib+jiDwOXAAQkbNubBsQit/4KCgoID8/n7y8PFasWEFsbCx79+4t1Ynp7NmzXH/99axbt45169axdetWTp8+TUGBZV6apKQkoqOjuXDhgq2Ofv369bYbSadPn7bdL7h48SJz5szhL3/5C2vWrKFJkyZ06NCB3bt3Ex0dDfye+Is6JhUdY9CgQezcuZNOnSzfuz/++GOp+wr5+flcvHiRU6dOsWnTJtq3b8/p06dtrYLOnTtH9erVSUtLsx2vVq1aPPPMMzz66KOePdlKqYDmahumvOILxphaBFni37FjBx06dABgypQpDBw4kLCwMJ544gnCw8PJy8ujZs2adttcfvnltGrVip49e1K7dm0mT57M4MGD2b17NzfeeCMXL15k7Nix1KhRgxo1ajB16lRycnJo3749ALfffjvXXXcdtWvXZsmSJfTo0YMHHniAnj17kp+fz+WXX07Hjh35n//5H6ZNm0Zqairjx4+nfv36DBs2jA8//JCOHTuycOFCZs+ebfv1sGPHDiZNmmQXa0JCAj169KB69eq88MILFBQUUKdOHQBOnTpFeHg4WVlZrF27lo4dO7Ju3TratWtnGwZaKRU6XG3O+TxwGhiDZXTOacB2EXnSu+FZeKo5p6vcac4ZTB555BGOHDnCpZdeyrXXXssXX3zBm2++SXh4uL9D86j09HQSExM5c+YMH3zwgb/DqbIcNR9U/uON5pzTAQPkAM8Dm7AMtlYllJxoZf369Q4nYgl2GRkZtG/fnpkzZ7J69Wpuu+22Kpf0AVq1auW10VaVqgpcTfyXiEiiiHQTketEZK6IXPRqZD4yduxY1q1bZ3fHOzY2ltjYWEetm1w2f/58YmJiaN++PaNHj7bNh3nixAm74WObN29uWy6q1y+pT58+fPHFF3ZlL730ku0GtatycnKoV68eTZo04ZNPPmHs2LG8/vrrPPjggy5tf++999p6IVd2KNt77rmHxo0b26rFSlqwYAFt27YlPj7e6To7duxg8ODBdo9ffvmlwjEpFSpcTfxZxpg3jDE9vRqNjx05coTo6GhiY2M9ut9Dhw6xYMECNm/eTGpqKgUFBSxbtgywjExZ1N540qRJPPTQQ7blGjUcz18/evRo2/ZFli1b5vYgamfPni3Veav4vY/ybNu2zXaD+euvv2br1q1uHb+48ePHl9l7+B//+Adr1qwpc/rJDh068Mknn9g9GjduXOGYlAoVrib+aCAFeNkYk2aM+asxpoUX4/KJ999/n/j4eNvykiVL6NKlC5s3b2bbtm2V2vfFixc5d+4cFy9eJDc3l2bNmrm03bvvvst1111Hp06duP/++ykoKGDEiBF8+umntl8EGRkZHD58mN69e7sVk6PEv337dqeJPy0tjV69etGhQwcSExM5evQoLVq0sA1t+8EHH9CpUyfS09PdigPghhtucDpZxqRJk0hPT2fgwIGsXLmSixcvEh8fT9u2bRkxYgS5ubll7vvEiRNMmjSJbdu28be//c3t2FTwWLlypd0kLJ06daJatWq88sordmXt27fHGMPu3bv9HXJgcNS4v6wH0B54C8h3d9uKPjzVgSs7O1s+//xz2/LUqVNtz3/99Vdp27at5OXlSWxsrPTs2dPhPnr16iXXXHNNqceaNWvs1nvppZekbt26csUVV8hdd93lcF8zZ86Uv//977blXbt2yeDBg22dvf785z/LkiVLRERk0KBB8tFHH4mIyN/+9jd5+OGH3X7/rVq1kn//+992ZZdddpnDCVXOnz8v7dq1kx9++MEWy4033mh7vX///rJjx45S27l6fkRE9u/fLzExMQ5jjYyMlOPHj8v+/fvtJm+ZMGGC3TlT/lOyw1AeePRREa+//rrccMMNUlBQYFf+xBNPSHx8fIXfazDw+EQsAMaYasAtWKZcvAF429NfQt6WkZHBI488Qrdu3Thy5IjdlW5YWBjnzp3j4YcfJicnh/r16zvcx7ffflvucU6dOsWqVavYv38/4eHh3HHHHbz77rsOp0os7uuvv2bLli386U9/Aixt74uqLoqqe4YOHcqyZcsqdPPy7Nmzdu/rwIED1K9f3+FkKh999BFdu3a1jaYZExNj6wAHzoeydeX8uOuqq66iZ09LLeOYMWNYsGABjzzyiMePo4JbWloaTz/9NBs2bLCby+Gbb75h+fLllaqarGpcSvzGmHnAKCxj5y8BxorIOW8G5g1dunThjjvu4F//+hfHjh1jypQpttfq1KlDamoqH3/8MW+//bbTnqy9e/d2ODnKCy+8wM033wzAV199RcuWLW2dp4YPH86GDRvKTfwiwrhx4xxWTwwdOpSHHnqIrVu3kpubS5cuXVx+30VKVvWUVb+/Y8cOu2Ns2bLF1sy1rKFsXTk/7irZoqoqtLBSnpWfn89dd93Fiy++SEREhK389OnTjB8/nnfeecet4darujITvzEmQkSygBPAdSJysNhrncU6mXowGTNmDOPHj6dXr152V7/79u3j6quvZtSoUTz22GMUFhY63N6VK9qIiAg2btxIbm4utWvX5uuvv6Zr1/IHM73ppptsCb5x48acPHmSnJwcIiMjqVevHn369OGee+4pdVO35NCzjhQWFpKbm2uX+Muq32/YsCGpqZY5crZs2cLSpUttX5QZGRlO71l444o/KyuL77//nuuvv5733nuPXr16efwYKrjNmDGDmJgY27AlRSZNmsTYsWNtvxiVRXk3dz8CEJFE4MMSr73hlYi8rFWrVhQUFNiuXoskJiYSHR1N586dOX/+fJlJtDzdunVjxIgRdO7cmQ4dOlBYWOjSeP3t2rVj7ty59OvXj44dO9K3b1+OHDlie3306NH8+OOPdonf0dCzjpw9exag1BV/UlISUVFRREVF2c32NXbsWFJSUujUqRPPP/884eHhtGvXDig9tG1FjB49muuvv569e/fSokWLMquuoqOjeeWVV2jbti2nTp1yuxmrqtrWrVvHihUrSg23vmTJEjIzM5kxY4afIgtgjir+ix7ANkfPHS178+Gpm7tFvv3221I3f4q8++67UrNmTQEkMjJS3n33XY8e29N27NghDz30ULnrHTp0SADJycnxQVQqFATCzd2TJ09KVFSUbNiwwa78559/lmbNmsnPP//ssfcb6Dx5c1ecPHe0HDScVRUkJyeTkJBAXp5laKLMzEzblXrxZp+BpH379sybN8/p6ydOnKBbt24899xzXHHFFToJi6pSXnvtNX755ZdSvwLPnDlDbm4uw4cPtyv/3//9X7ebQVdF5SX+WsaYtliGayj+HCyzclUp06dPL9VGPDc3l+nTpwds4i9Pw4YN6du3L/feey/z58/3dziqCqvhZu92T3jiiSd44oknfH7cYFde4q8DFJ+MsvjzoL3idyYrK8ut8mDx6quv8uqrr/o7DKVUgCgz8YtIlI/iCAgRERFkZmY6LFdKqaoiqMbU97bExETbGPZF6tSpU+Z4MUopFWw08RcTHx9PUlKSbUKWyMhIkpKSgrZ+XymlHHF5yIZQER8fz6JFi4CqNxGLUkqBXvErpVTI0cSvlFIhRhO/UkqFGE38Sqmg9vrrr5fqudu+ffsKTboyYcIEuwlcmjRpYhsHq6zX3PXbb78xefJkNm7c6HDZ1XUqSm/uKqU84raVnt3fR8NcW2/Hjh107tzZtnz+/HkyMjJo3bq128d86623bM/3799P7969SUpKKvc1d7322mucP3+e7777ju7du5dadnWditIrfqVUUNu+fbtd4t+xYwetW7cmLCyszG3KmjI0OzubAQMGMGPGDIYMGeLya64e5/PPPyc6Oto2h3XJZVfXqShN/EqpoLZz506GDx9uG1584MCBdOzYscxtzp8/z2233eYwKefm5nLrrbcycuRI7r//fpdfc/U458+fp6CggK1btxIbG1tq2dV1KkOrepRSQevAgQM0atSIPXv22MomT55My5Ytbcvvvvsuzz77bKltjxw5wqhRo9i0aZOtrKCggFGjRtGmTRvmzJljt35Zr7lznH379lFQUECbNm2oXr06e/bssVt2dZ3K8FriN8ZcBfwTuBLLgG5JIvKyMeZy4H0gCsgARorIKW/FoZSqunbs2EFMTIxd2a5duxg6dKhtecyYMaWmPc3KymLIkCGlRqx94IEHyM/Pt3XidPU1d45z/Phx0tLSWLVqlcNlV9epDG9W9VwEHhaRdkB34EFjTDvgceBrEbka+Nq6rJRSbtu+fbttZrgiO3fudDqlaJG9e/fyj3/8w25KxtmzZ7Nlyxb+9a9/lZpPuqzX3D3O4cOHuf322yksLOTUqVOlll1dpzK8lvhF5IhY5+QVkRxgN9AcGIplwnasf2/zVgxKqaptx44ddon/5MmTiAhNmjQpc7u+ffvSo0cP23JGRgazZs3ixIkT9OrVy9Zk88477yzztfKUPM7FixfZunUrR48eZdKkSYSFhZVadmWdyjLig8kTjDFRwDdAeyBLRMKt5QY4VbRcYpsEIAEgIiKii6Phkr2laD5eHatHKed2795N27Zt/R2GsnL072GM2SIiXUuu6/VWPcaYesAKYKqI/Fr8NeuckA6/eUQkSUS6ikjXRo0aeTtMpZQKGV5N/MaY6liSfrKIfGgtPmaMaWp9vSnwizdjUEopZc9rid9ajfMmsFtEis8GvhoYZ30+DvDMbWqllFIu8WY7/p7AWGCHMSbFWvZX4FlguTFmIpAJjPRiDEoppUrwWuIXke8A4+Tlm7x1XKWUUmXTIRuUUj4TFxdnazWn/EcTv1IqaIWFhdGpUydiYmK45pprePHFFyksLCxzm4yMDN577z0fRRiYNPErpXwiOTmZjRs3sn79eqKiokhOTq70PmvXrk1KSgo7d+5kzZo1fPbZZ8yePbvMbTTxa+JXSvlAcnIyCQkJ5OXlAZCZmUlCQoJHkn+Rxo0bk5SUxMKFCxERMjIy6N27N507d6Zz585s2LABgMcff5xvv/2WTp06MX/+fKfrVWkiEvCPLl26iC/FxsZKbGysT4+pVLDZtWuXy+tGRkYWdda0e0RGRlYqhrp165Yqa9CggRw9elR+++03OXfunIiIpKWlSVEeWbt2rQwaNMi2vrP1go2jfw9gszjIqToss1LK67Kystwq94T8/HwmT55MSkoKYWFhpKWlVWq9qkQTv1LK6yIiInA03lZERIRHj5Oenk5YWBiNGzdm9uzZXHnllfz4448UFhZSq1Yth9vMnz/fpfWqEq3jV0p5XWJiInXq1LErq1OnDomJiR47xvHjx5k0aRKTJ0/GGMOZM2do2rQp1apV45133qGgoACA+vXrk5OTY9vO2XpVmSZ+pZTXxcfHk5SURM2aNQGIjIwkKSmJ+Pj4Su333LlztuacN998M/369WPmzJmAZeKUJUuWcM0117Bnzx7q1q0LQMeOHQkLC+Oaa65h/vz5TterynwyLHNlde3aVTZv3uyz4+mwzEqVryLDMuv/Le9xZ1hmreNXSvmMJvzAoFU9SikVYjTxV5COOaKUClaa+JVSFRYM9whDgbv/Dpr4lVIVUqtWLU6cOKHJ389EhBMnTrjV/0Bv7iqlKqRFixYcPHiQ48eP+zuUkFerVi1atGjh8vqa+B3QlgdKla969eq0bNnS32GoCtCqHqWUCjGa+JVSKsRo4ldKqRCjiV8ppUKMJn6llAoxmvi9xFM9e7WHsFLK0zTxK6VUiNHEr5RSIUYTv1JKhRhN/EopFWI08SulVIjRxK+UUiFGE79SSoUYTfxKKRVivJb4jTGLjTG/GGNSi5XNMsYcMsakWB+3eOv4oUQ7eSml3OHNK/63gQEOyueLSCfr499ePL5SSikHvJb4ReQb4KS39q+UUqpi/FHHP9kYs91aFXSZH46vlFIhzdeJ/1XgD0An4AjworMVjTEJxpjNxpjNOqenUkp5jk8Tv4gcE5ECEe5lgEUAABN/SURBVCkEFgHXlbFukoh0FZGujRo18l2QSilVxfk08RtjmhZbHAakOltXKaWUd1zirR0bY5YCccAVxpiDwEwgzhjTCRAgA7jfW8dXSinlmNcSv4iMdlD8preOp5RSyjXac7cCkpOT2bhxI+vXrycqKork5OQK7Uc7Ximl/EETv5uSk5NJSEggLy8PgMzMTBISEiqc/JVSytc08btp+vTp5Obm2pXl5uYyffp0P0WklFLu0cTvpqysLLfKlVIq0Gjid1NERIRb5UopFWg08bspMTGROnXq2JXVqVOHxMREP0WklFLu0cTvpvj4eJKSkqhZsyYAkZGRJCUlER8f7+fIlFLKNV5rx1+VxcfHs2jRIgDWrVvn32CUUspNesXvBZ5q5++p/SilVHGa+D3MU+38tb+AUspbNPF7mKfa+Wt/AaWUt2ji9zBPtfPX/gJKKW/RxO9hnmrnr/0FlFLeoonfwzzVzl/7CyilvEUTv4d5qp2/9hdQSnmLtuP3Ak+189f+Akopb9Ar/iCnbf2VUu7SxB/EtK2/UqoiNPEHMW3rr5SqCE38Qcwfbf11ukilgp8m/iCmbf2VUhWhiT+IaVt/pVRFaOIPYtrWXylVEdqOP8hpW3+llLv0it9PtP29UspfNPH7gba/V0r5kyZ+P9D290opf9LE7wc61r5Syp808fuBtr9XSvmTJn4/0Pb3Sil/0sTvB9r+XinlT9qO30+0/b1Syl+8dsVvjFlsjPnFGJNarOxyY8waY8w+69/LvHV8pZRSjnmzqudtYECJsseBr0XkauBr63JQWrdunV6pK6WCktcSv4h8A5wsUTwUWGJ9vgS4zVvHV0op5Zivb+5eKSJHrM+PAlf6+PhKKRXy/NaqR0QEEGevG2MSjDGbjTGbjx8/7sPIlFKqavN1q55jxpimInLEGNMU+MXZiiKSBCQBdO3a1ekXhHLutpXlr/PRMO/HoZQKLL5O/KuBccCz1r+rfHz8KuOCMbbn4qAMoDtp5e9IrvZgVEqpYOC1xG+MWQrEAVcYYw4CM7Ek/OXGmIlAJjDSW8ev6uYVS+oHGGMte9df4SilgojXEr+IjHby0k3eOmYg8VRTT20yqpTyNB2yQblMJ49RqmrQxK9copPHKFV1aOJXLtHJY5SqOnSQtgD0rNnn7xBK0cljlKo69IpfuUQnj1Gq6tDEr1yik8coVXVoVU8VkOCD9vtFk8RMnDiRvLw8IiMjSUxM1MljlApCmvhDRJK1k1fJL4niwzp8NyMOgF5z1tmtUzSsg6uTx8TFxZW7jlLKf7SqRymlQowmfqWUCjFa1eNjxZtqZnGuVJmvLR/++8BufR2UASA6OKpSVYle8SulVIjRK/4Q58oon0E7MbJSyiFN/H7ki2aYSilVklb1KKVUiNErflWu4m39U7NLl4FO4ahUMNHEr8qlLX+UqlpCLvFrr1L36Q1gpaoWreNXSqkQo4lfKaVCjCZ+pZQKMZr4lVIqxGjiV0qpEKOJXymlQowm/hCwjdVkkcJ+NvEccWxjdYXWcUVycjIbN25k/fr1REVFkZycXNnwlVIepom/itvGalbyJAVcAOA0h1nJk3aJ3ZV1XJGcnExCQgJ5eXkAZGZmkpCQoMlfqQCjib+K+5J55HPeriyf83zJPLfWccX06dPJzc21K8vNzWX69OluRq2U8qaQ67kbak5zpNxyV9YpUtaIollZWW6VBxvt9V05ev4Ch17xV3HhNC233JV1XBEREeFWuVLKPzTxV3H9mEZ1atmVVacW/Zjm1jquSExMpE6dOnZlderUITEx0c2olVLepIm/iruWIQxjLmHUACCcZgxjLtcyxK11XBEfH09SUhI1a9YEIDIykqSkJOLj4z30bpRSnqB1/CHgWobwX5YDzuvoXVnHFfHx8SxatAjQulylApVfrviNMRnGmB3GmBRjzGZfHVfbmCullH+v+PuISLavDuasjTmgVRFKqZASMnX82sZcKaUs/JX4BfjSGLPFGJPgaAVjTIIxZrMxZvPx48crfcCq3sZcKaVc5a/E30tEOgMDgQeNMTeUXEFEkkSkq4h0bdSoUaUPqG3MlVLKwi+JX0QOWf/+AqwErvP2MbWNuXddMMb2kPXrkfXr7couGFP+TvwoLi7O1rNUuU/PX3Dx+c1dY0xdoJqI5Fif9wOe9vZxi27gTpw4kby8PCIjI0lMTNQbux6iE7IrFTz80arnSmClsVwBXgK8JyKf++LA2sbcv1y56q8h4oNIlAptPk/8IpIOXOPr4yr/K/6rwJlA/VVQ1AckLy+PqKgo/bXoJj1/gSVkmnMqVVE6z0Dl6PkLPJr4lSqH9gGpHD1/gUfH6vGgZ80+f4egvED7gFSOnr/Ao1f8KqCUbALq6OFr2gekcvT8BR694neRawmn/JuXqmyBeAM4MTGRhIQEu+oK7QPiOj1/gUcTv4tcSUiqatI+IJWj5y/waOIPEa6MsV+Zcfi9sZ9Aon1AKkfPX2DROn4VsJIYQ5K1F7AzwThUQKDFHGjxBJKqem408augU964QO7SCXoqR89f8NGqHhV0yhsXyJ2bvzpBT+Xo+QtOmvjR9vehrKzORZq4yqfnLzhp4ldVjitf5I/L1YB2LqosPX/BKeQSv7YoUPB7v4yrAEcp6iodJdQlERERZGZmOixXgUtv7qqQNI805pFGd16gOrXsXqtOLbrzgp8iCy46wVFw0sSvQtq1DGEYcwmjBgDhNGMYc7mWIX6OLDjEx8eTlJREzZo1AYiMjCQpKUnr9wNcyFX1KFXStQzhvywHqmbnM2/TzlnBR6/4VUDaxmqySGE/m3iOOLaxukLrBJpAa/MeaPEEkqp8bjTxq4CzjdWs5EkKuADAaQ6zkiftErsr6wSaQJuQJNDiCSRV/dxo4lcB50vmkc95u7J8zvMl89xap7JKDgdd2V7CvpiQxJ1hrXWCFOeq+rmp8nX82jkrMBWNweOoTv00RxxuU7zclXVcOVZZSo7IOtBW/jtXegkXJVtnLduzMjNLfYn0tf5dU6ys5ET0RWPIVLRe3ddt8F2Js7LvyZ39lLWOO+fGUzH7UpVP/Cr4hNOU0xx2WO7OOoGi6AukAXEOY25AM+axzq6sskNRuELb4DtX1c+NJn4VcPoxjZU8aVeVU51a9GOaW+sEmsrGXPJXgTgod2ciG50gxbmqfm60jl8FHFfa1gdK+3t36tQDJeYi2gbfuap+bvSKXwUkV9rWB0L7e3dnZqtMzCWP5ag6yF3O2uC7cuO65D2HqqYq90/QK37lc75sfx+Mbf293YfB062V3DlWWb+IPNVu3pX9+PJYgUiv+JVPOWt/D3i8ysOdYwVKj11XYvblOSxL8aT9pYMyd5Q3rr+r+10KPADkWZczMzNJGDOGi2PGMBrLrxRPzSEQzHMR6BW/8ilftL/3x7E8JVD6MBQNYlfWw1MuGMNfx4xx2G7+r2PGuPVl8hSQW6Is11pexFNt9IO5rb9e8Sufcqf9fTAdy1M83YfBEW/cK6isA26WV3Q/F4wptz+Fq/cugnkuAk38yqd82f4+mNr6F6mqfRjKMo3WzudF8PB+ip6XtU7xXxiOmsyCpcoomNv6a1WP8ql+THM4/r032t/78lie4krMwfi+yvM0UKdEWR1ruavKm1+h6MvDE8eC4J6LwEgQNMnq2rWrbN68uULb6pANgWcbq1nBXyngAuE0ox/TvHZT0pfH8hRXYg7G91UeT72nsvYzjdaA5Sbw/VhuAkdgSfqjKxBzWfsJhOauxpgtItK1VLkmfuUPFR0/J9CP5SmuxByM76s8nnpPlTl/RV8OrnI0rpInVeYLxFni90sdvzFmAPAyEAa8ISLP+iMOpTxpG6v5knmc5gjhNPX7Vbin4gmk9+WLWFy9n+DqfsqKufgvkKew3IS+ior/AnGVzxO/MSYMeAXLF+VB4L/GmNUissvXsSjlKUVt64uaWfqrbb2n4wmk9xVIsRT/ciirVZQrMRf1PShqGJplXQbvJX9/3Ny9DvhJRNJF5AKwDBjqhziU8phA6zPgqXgC6X0FUiyuKi/meaQxlWYO+x5MpZlH+0sU5/M6fmPMCGCAiNxrXR4LdBORySXWSwASrIvRwF4Hu7sCyPZiuN6gMXufP+LtUsZrW1zY3tMxVzYeV/aTiW/Psyfek68/G67EXN46lYk5UkQalSwM2Hb8IpIEJJW1jjFms6MbF4FMY/a+YIsXNGZf0Zgt/FHVcwj7/hQtrGVKKaV8wB+J/7/A1caYlsaYGsAoCIIhE5VSqorweVWPiFw0xkwGvsDSnHOxiOys4O7KrAoKUBqz9wVbvKAx+4rGTJB04FJKKeU5OlaPUkqFGE38SikVYoIm8Rtjoo0xKcUevxpjphpjZhljDhUrv8XPcS42xvxijEktVna5MWaNMWaf9e9l1nJjjFlgjPnJGLPdGNM5QOL9uzFmjzWmlcaYcGt5lDHmXLFz/Zqv4y0jZqefA2PME9ZzvNcY0z+AYn6/WLwZxpgUa7nfz7Mx5ipjzFpjzC5jzE5jzBRreSB/lp3FHLCf5zJi9u7nWUSC7oHlpvBRIBKYBTzi75iKxXYD0BlILVb2PPC49fnjwHPW57cAnwEG6A78ECDx9gMusT5/rli8UcXXC7Bz7PBzALQDfgRqAi2Bn4GwQIi5xOsvAk8FynkGmgKdrc/rA2nWcxnIn2VnMQfs57mMmL36eQ6aK/4SbgJ+FpHSsyD4mYh8A5wsUTwUWGJ9vgS4rVj5P8ViIxBujPHpbBqO4hWRL0XkonVxI5a+FgHDyTl2ZiiwTETyRGQ/8BOWYUN8qqyYjTEGGIll2JaAICJHRGSr9XkOsBtoTmB/lh3GHMif5zLOszMe+TwHa+Ifhf1/ksnWn3GLi356BpgrRaRoXryjwJXW582xny3uIGX/o/vDPViu5Iq0NMZsM8asN8b09ldQTjj6HATDOe4NHBOR4mOIB8x5NsZEAdcCPxAkn+USMRcXsJ9nBzF77fMcdInfWDp9DQH+ZS16FfgD0Ak4guUnc8ASy++1oGhDa4yZDlwEkq1FR4AIEbkWmAa8Z4y51F/xlRBUn4MSRmN/IRMw59kYUw9YAUwVkV+Lvxaon2VnMQfy59lBzF79PAdd4gcGAltF5BiAiBwTkQIRKQQW4Yef8S44VvSz1/r3F2t5wA5fYYwZDwwG4q3/wbH+vDxhfb4FS/2ie7NWeEkZn4OAPccAxphLgOHA+0VlgXKejTHVsSSjZBH50Foc0J9lJzEH9OfZUcze/jwHY+K3uzoqUY84DEgttYX/rQbGWZ+PA1YVK7/b2iKiO3Cm2M9ovzGWiXIeBYaISG6x8kbGMp8CxphWwNVAun+itFfG52A1MMoYU9MY0xJLzJt8HV8Zbgb2iMjBooJAOM/W+w5vArtFpPi4xwH7WXYWcyB/nsuI2bufZ3/e0Xb3AdQFTgANipW9A+wAtltPSlM/x7gUy0+zfCz1bxOBhsDXwD7gK+By67oGy6Q0P1vfQ9cAifcnLPWIKdbHa9Z1bwd2Wsu2ArcG0Dl2+jkAplvP8V5gYKDEbC1/G5hUYl2/n2egF5ZqnO3FPge3BPhn2VnMAft5LiNmr36edcgGpZQKMcFY1aOUUqoSNPErpVSI0cSvlFIhRhO/UkqFGE38SikVYjTxK6VUiNHEr5QKScaY3sYYMcYsKFG+3lrukbGGrPvaboy52bpc3Trs8l7rUMw/GmM+MMa0K2c/nxljJpUoM8aYdGNMrHX5K2PMSWOZ3tYpTfxKqVDVGdgCtC8qMMaMxDI88i8i4skhJ3qIyFfW528BHYFuIhKDZTyet4DocvaxGBhfoiwOKAS+ARCRm7F0+CqTJn6lVKjqjGXAttYAxphawEws4yZt9cYBjTFXYxmCYaKInAbLYHci8qmIrLSu081YJmfZYn0Msm6+CvijMaZtsV1OAN4SN3viauJXSoWqzliGQM42xlwJ/D8sSb8x1sRvjHnKGPN6Wc8dLZfhWmCfiJxy9KKxzA72GnCXiHTBMrDc68aYcBG5gOWLaoJ13fpY5kNY4mhfZdHEr5QKOdar+9ZYZrP6ERgA3AX8HcsXwlbrAG4tgcPOnlv3ZbfsZhztjGVqxTRjzMtAD+u+PjOWqTg/wzKWzx+tmywGxlqPeSfwHyk2wJ+rLnF3A6WUqgKuAdJF5DdjzI/AAiyD5Z0zxnTCcsV/DZZB9b4v4zkOlsuyDbjaegV/WkR2AZ2sN2O7YhnsbruI3OBoYxH50RhzGMvw9BOAlyrw3vWKXykVkjrzez3+x8AMEVlqHZ5ZxDKtYRssF8dryniOg2WnxDLL2ipgkTGmQbGX6lr/bsDyxdCn6AVjzJ+swzcXWYxlTt7W/D4stls08SulQlFnLFffiMg+ESlq0nltUTnQAPgIy1W4s+cl13PFeGAP8F9rc87vgC7AAmvd/xBgprWZ524sSb74vt8DYoD3rPX+btOqHqVUyBGR+5yUrwBWWK+wG2MZGz/M2XNjzEclXlvpwrEvADOsD0ev/xdLM01n258Capd3nLLoePxKKeVFxpijwDHg4WJt+b11rK+wzNX7tIi85XQ9TfxKKRVatI5fKaVCjCZ+pZQKMZr4lVIqxGjiV0qpEKOJXymlQowmfqWUCjGa+JVSKsRo4ldKqRDz/wE1C7gkZCbwTAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "if __name__==\"__main__\":\n",
    "    start = time.time()\n",
    "    \n",
    "    cutmin = 0\n",
    "    cutmax = 1\n",
    "    \n",
    "    for c in range(cutmin,cutmax):\n",
    "        data = get_data_from_files(c)\n",
    "    \n",
    "    plot_data(data)\n",
    "    elapsed = time.time() - start\n",
    "    print(\"Time taken: \"+str(elapsed))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}