{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Machine Learning for $t\\bar{t}Z$ Opposite-sign dilepton analysis \n", "This notebook uses ATLAS Open Data http://opendata.atlas.cern to show you the steps to implement Machine Learning in the $t\\bar{t}Z$ Opposite-sign dilepton analysis, following the ATLAS published paper [Measurement of the $t\\bar{t}Z$ and $t\\bar{t}W$ cross sections in proton-proton collisions at $\\sqrt{s}$ = 13 TeV with the ATLAS detector](https://journals.aps.org/prd/pdf/10.1103/PhysRevD.99.072009).\n", "\n", "The whole notebook takes less than an hour to follow through.\n", "\n", "Notebooks are web applications that allow you to create and share documents that can contain for example:\n", "1. live code\n", "2. visualisations\n", "3. narrative text\n", "\n", "Notebooks are a perfect platform to develop Machine Learning for your work, since you'll need exactly those 3 things: code, visualisations and narrative text!\n", "\n", "We're interested in Machine Learning because we can design an algorithm to figure out for itself how to do various analyses, potentially saving us countless human-hours of design and analysis work.\n", "\n", "Machine Learning use within ATLAS includes: \n", "* particle tracking\n", "* particle identification\n", "* signal/background classification\n", "* and more!\n", "\n", "This notebook will focus on ROC curves.\n", "\n", "By the end of this notebook you will be able to:\n", "1. run machine learning algorithms to classify signal and background\n", "2. know some things you can change to improve your machine learning algorithms\n", "\n", "Feynman diagram pictures are borrowed from our friends at https://www.particlezoo.net" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Introduction (from Section 1)\n", "\n", "Properties of the top quark have been explored by the\n", "Large Hadron Collider (LHC) and previous collider experiments in great detail. \n", "\n", "Other properties of the top quark are\n", "now becoming accessible, owing to the large center-of-mass energy and luminosity at the LHC.\n", "\n", "Measurements of top-quark pairs in association with a Z boson ($t\\bar{t}Z$) provide a direct probe of the\n", "weak couplings of the top quark. These couplings\n", "may be modified in the presence of physics beyond the\n", "Standard Model (BSM). Measurements of the $t\\bar{t}Z$ production cross sections, $\\sigma_{t\\bar{t}Z}$, can be used to\n", "set constraints on the weak couplings of the top quark. \n", "\n", "The production of $t\\bar{t}Z$ is often an important\n", "background in searches involving final states with multiple\n", "leptons and b-quarks. These processes also constitute an\n", "important background in measurements of the associated\n", "production of the Higgs boson with top quarks.\n", "\n", "This paper presents measurements of the $t\\bar{t}Z$ cross section using proton–proton (pp) collision data\n", "at a center-of-mass energy $\\sqrt{s} = 13 TeV.\n", "\n", "The final states of top-quark pairs produced in association with a\n", "Z boson contain up to four isolated, prompt leptons. In this analysis, events with two opposite-sign\n", "(OS) leptons are considered. The dominant backgrounds\n", "in this channel are Z+jets and $t\\bar{t}$, \n", "\n", "(In this paper, lepton is used to denote electron or muon, and prompt lepton is used to denote a lepton produced in a Z or W\n", "boson decay, or in the decay of a τ-lepton which arises from a Z or W boson decay.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Data and simulated samples (from Section 3)\n", "\n", "The data were collected with the ATLAS detector at a proton–proton (pp) collision\n", "energy of 13 TeV. \n", "\n", "Monte Carlo (MC) simulation samples are used to model the expected signal and background distributions\n", "in the different control, validation and signal regions described below. All samples were processed through the\n", "same reconstruction software as used for the data. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Opposite-sign dilepton analysis (from Section 5A)\n", "\n", "The OS dilepton analysis targets the $t\\bar{t}Z$ process, where both top quarks decay hadronically and the Z boson\n", "decays to a pair of leptons (electrons or muons). Events are required to have exactly two opposite-sign leptons.\n", "Events with additional isolated leptons are rejected. The leading (subleading) lepton is required to have a\n", "transverse momentum of at least 30 (15) GeV.\n", "\n", "The OS dilepton analysis is affected by large backgrounds from Z+jets or $t\\bar{t}$ production, both characterized\n", "by the presence of two leptons. \n", "\n", "The signal region\n", "requirements are summarized in Table 1 below." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "| Variable | Selection |\n", "|------|------|\n", "| Leptons | = 2, opposite sign |\n", "| $p_T$ (leading lepton) | > 30 GeV |\n", "| $p_T$ (subleading lepton) | > 15 GeV |\n", "\n", "Table 1: Summary of the event selection requirements in the OS dilepton signal regions.\n", "\n", "This is a subset of Table 2 of the ATLAS published paper [Measurement of the $t\\bar{t}Z$ and $t\\bar{t}W$ cross sections in proton-proton collisions at $\\sqrt{s}$ = 13 TeV with the ATLAS detector](https://journals.aps.org/prd/pdf/10.1103/PhysRevD.99.072009)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents: \n", "\n", "[Running a Jupyter notebook](#running)
\n", "[To setup first time](#setupfirsttime)
\n", "[To setup everytime](#setupeverytime)
\n", "  [File path](#fraction)
\n", "  [Get data from files](#get_data_from_files)
\n", "\n", "[Machine learning](#MVA)
\n", "  [Training and Testing split](#train_test_split)
\n", "  [Training](#MVA_training)
\n", " \n", "[Going further](#going_further)
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running a Jupyter notebook\n", "\n", "To run the whole Jupyter notebook, in the top menu click Cell -> Run All.\n", "\n", "To propagate a change you've made to a piece of code, click Cell -> Run All Below.\n", "\n", "You can also run a single code cell, by using the keyboard shortcut Shift+Enter." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## First time setup on your computer (no need on mybinder)\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 open on mybinder, you don't need to run this cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import sys\n", "!{sys.executable} -m pip install --upgrade --user pip # update the pip package installer\n", "!{sys.executable} -m pip install -U pandas sklearn --user # install required packages" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "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.\n", "\n", "We're going to be using a number of tools to help us:\n", "* pandas: lets us store data as dataframes, a format widely used in Machine Learning\n", "* numpy: provides numerical calculations such as histogramming\n", "* matplotlib: common tool for making plots, figures, images, visualisations" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [], "source": [ "import pandas as pd # to store data as dataframes\n", "import numpy as np # for numerical calculations such as histogramming\n", "import matplotlib.pyplot as plt # for plotting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## File path" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#csv_path = \"~/13tevatlasopendataguide2020/docs/visualization/CrossFilter/13TeV_ttZ.csv\" # local \n", "csv_path = \"http://opendata.atlas.cern/release/2020/documentation/visualization/CrossFilter/13TeV_ttZ.csv\" # web address" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get data from files" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "define function to get data from files\n", "\n", "The datasets used in this notebook have already been filtered to include exactly 2 leptons per event, so that processing is quicker." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
ChannelNJetsMETMllLepDeltaPhiLepDeltaRSumLepPtNBJets
145370055.7690.781.321.40139.591
1453710497.1691.320.991.37175.511
1453720525.6688.161.671.7091.101
1453730564.0691.002.692.7122.762
1453740552.3375.130.101.12152.971
...........................
1542401414.7865.870.660.67196.762
1542411833.9791.720.600.98177.023
1542421556.0762.151.391.4988.680
154243110143.1783.462.782.7960.341
1542441742.1777.311.381.5593.172
\n", "

8875 rows × 8 columns

\n", "
" ], "text/plain": [ " Channel NJets MET Mll LepDeltaPhi LepDeltaR SumLepPt \\\n", "145370 0 5 5.76 90.78 1.32 1.40 139.59 \n", "145371 0 4 97.16 91.32 0.99 1.37 175.51 \n", "145372 0 5 25.66 88.16 1.67 1.70 91.10 \n", "145373 0 5 64.06 91.00 2.69 2.71 22.76 \n", "145374 0 5 52.33 75.13 0.10 1.12 152.97 \n", "... ... ... ... ... ... ... ... \n", "154240 1 4 14.78 65.87 0.66 0.67 196.76 \n", "154241 1 8 33.97 91.72 0.60 0.98 177.02 \n", "154242 1 5 56.07 62.15 1.39 1.49 88.68 \n", "154243 1 10 143.17 83.46 2.78 2.79 60.34 \n", "154244 1 7 42.17 77.31 1.38 1.55 93.17 \n", "\n", " NBJets \n", "145370 1 \n", "145371 1 \n", "145372 1 \n", "145373 2 \n", "145374 1 \n", "... ... \n", "154240 2 \n", "154241 3 \n", "154242 0 \n", "154243 1 \n", "154244 2 \n", "\n", "[8875 rows x 8 columns]" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "all_data = pd.read_csv(csv_path) # read all data\n", "signal_df = all_data[all_data['type']==1].drop(['type','weight'], axis=1) # get signal dataframe\n", "background_df = all_data[(all_data['type']!=0) & (all_data['type']!=1)].drop(['type','weight'], axis=1) # background dataframe\n", "signal_df # print the dataframe to take a look" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Machine learning\n", "\n", "Organise data ready for BDT" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# for sklearn data is usually organised \n", "# into one 2D array of shape (n_samples x n_features) \n", "# containing all the data and one array of categories \n", "# of length n_samples \n", "\n", "X = np.concatenate([signal_df.values, background_df.values]) # concatenate the list of MC dataframes into a single 2D array of features, called X\n", "y = np.concatenate([np.ones(signal_df.shape[0]), np.zeros(background_df.shape[0])]) # concatenate the list of lables into a single 1D array of labels, called y" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The Training and Testing split\n", "One of the first things to do is split your data into a training and testing set. This will split your data into train-test sets: 75%-25%. It will also shuffle entries so you will not get the first 75% of X for training and the last 25% for testing. This is particularly important in cases where you load all signal events first and then the background events.\n", "\n", "Here we split our data into two independent samples. The split is to create a training and testing set. The first will be used for training the classifier and the second to evaluate its performance.\n", "\n", "We don't want to test on events that we used to train on, this prevents overfitting to some subset of data so the network would be good for the test data but much worse at any *new* data it sees." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from sklearn.model_selection import train_test_split\n", "\n", "# make train and test sets\n", "X_train, X_test, y_train, y_test = train_test_split(X, y, \n", " random_state=492 )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Training your machine learning algorithm\n", "\n", "We'll use SciKit Learn (sklearn) in this tutorial. Other possible tools include keras and pytorch. \n", "\n", "After instantiating our GradientBoostingClassifier, call the fit() method with the training sample as an argument. This will train the tree, now we are ready to evaluate the performance on the held out testing set.\n", "\n", "A useful plot to judge the performance of a classifier is to look at the Receiver Operarting Characteristic (ROC) curve directly." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEGCAYAAABo25JHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+4UlEQVR4nO3deXwV1fn48c9zl+RmDwmEBJKw78iiyOJCtdqKFtFWW9Ra67etuFRrV+v3p69Wrd20+q1W/VqsftXW1qVWRIv7UhdUNkH2VZYAIQuQ/SZ3eX5/zE0IEOACublJ7vP2lVfuzJw78wzIPDPnzDlHVBVjjDGJyxXvAIwxxsSXJQJjjElwlgiMMSbBWSIwxpgEZ4nAGGMSnCfeARytnj17av/+/eMdhjHGdCmLFy+uUNVebW3rcomgf//+LFq0KN5hGGNMlyIiWw61zaqGjDEmwVkiMMaYBGeJwBhjEpwlAmOMSXCWCIwxJsHFLBGIyGMiUiYiKw6xXUTkfhHZICKficiJsYrFGGPMocXyieBxYNphtp8LDIn8zAL+N4axGGOMOYSY9SNQ1fdEpP9hilwAPKnOONgfi0i2iBSo6s5YxWSMMTGnCoF68FdDYzX4q1B/FcGGaoKhEIFgiGAoTG2gkbqgn0AwRCAUpCkUIhgOEQgGCYbDBMLO+ppQA8tDO5BQI18YciHnnXZFu4cczw5lfYFtrZZLIusOSgQiMgvnqYHi4uIOCc4Y04ZwGBp2Q/V20LBz0UNBcZbRVuvCh/is+8qGAs4FMxyAcAjCQWdbOAQaWfZXQagp8n3234+z4rDrwqqEwko4HCYUDhMOhQlpmFAYZ10oRFDDThlVAuEQocY9NIUbadIQwXCQXRJgmzvIXglRK2EaCaEoSni/32ERwkBQwC9Ck0jL7yaBgAhBEQJAo+soK2RckLP1Pc6jeyWCqKnqbGA2wIQJE2wmHWOOhSpUlTh3q8FGqNkJQb9zMQ42gn8v1O921tWUQqDBuUCHmqC2DBproGrbEQ/TbJPXwzMZGYQEIpd9NHKhVJrXCdpqe1CE3S4XTSKUe9yUud34VHErqDhlWk6n+fvsW9/yW1qXaeN7IgcH3HxddgPets/JrZAaduHVJFwqOP+5cEnzJ8GFC7d4cEsSbknC5/aR5krB4/bhcXnxurwkuTwku5PJ9GbgdXvwuD143G6S3G48Lg9Jbhdej5cktwef18uQ7AGM6DUC8WVF/ed/NOKZCLYDRa2WCyPrjDGHEg47F/DqHbB3C9SVO3fNB9xNh0NBwnu3QOVGtMm58HurPo/qEAGXj0ZPOvWebIJ4CIqHJnxUufvgzzyBWlJZ7x1KXchLIKQ0hZVACJpCSlMImsJhGkOwOXsBlZnrcIeSIXLJdK7Isu9HIYwbEFQj60I+NOSDgAd1NVET7IkbH26X4Ha5cIvgcrnwugS3uPatd4Hb5cLjcjm/RXC7nTIed/N6iZSRlmWPa992j8tFkttDWpKPNG8qqd4kUrzJeNxuhmQPYWDWQKStJNLFxTMRzAWuF5GngUlAlbUPmG6vYQ9UrIc9W6CpxrkbDzVBU51zca/dBY21EGp07tLrKiDYgIaCEAogIX9Uh1EVqkljg/alUjNpIodKHUYTXlaG+9GEl0Y8+EmmTLNpwkOjJlFOFs5FGjwuIcnjItnjItnjJkmdz0keF8luF8nJ7n3bvW6S3C4yvPvKlFatpGd4IN8s+p993/O49//sdZHkdrX63Xq7s647Xng7m5glAhH5B3AG0FNESoBfEnngUtWHgXnAecAGoB74r1jFYkxMqToXcf9e5+K9d6tzdx5sdKpT9myGncugvgIqN0Tqxw8W8OVSn5xHvaRQH/ZQH0qjPNSb3QEv1U3QpG6CuAjiZl24kK1SSDgzn54Z6fTKSiUjJZmkJA9ej5dkr2f/i7jHRa/I8vADLsI+r4skt7PcfBFOcjt3ycfj/X9VMCF3FN87feBx7cfEXizfGrr0CNsV+H6sjm9Mu/NXOXfoez6HTf9xLv7la2BXm11l9hN2edmbOZzNRd9mg6s/JY2pfObvTUlNiF11IfyaRMC/759jitdNQZaP/Fwf+Vk+53NWCgWZzvLlWT5yUpNwuTrn3XIgHGBH7Q6m9T/cG+Sms+gSjcXGdIhwyLm471jiNI6WrYLqnVBdAvV7nN/N3EmEM/oScCUjyTl83mc6G5NHUu6H3fVhNjf42FCbRHmDUE4Wist57i2FTJ+HgqwU8rN8TCj00Tuz+ULva1mf6fN06SqRnbU7CWmI4kx7y68rsERgurdwyLmL91dBzQ7nTZj63fuqcUIB502Y3Ztg+xIINuz39fr0flQn51OZXEx57mms0WI2+dN5p24A5aXufQWrnF+5aUnOBb2nj/GDIhf2yF18fpaP/Ewfacnd/5/d1pqtABRnWCLoCrr//5GmewuHnYt3qMm5sFdudOrjy1Y6F/aqEqdTz2H43Wlsdfdji+d0lrkGsaXBxyeBwVSRRqM/CQCXQF6Gj95ZPgoKfHwla/+7+IIsH3mZySR73Ic9VqLYUu3MgWJPBF2DJQLTNfiroHQFbF/kVNeUrXIaYesqIFB3cPmkDCg6GQpPhvwTIDkD0nsTTMpk2a5GXt0svLp2D9trQrhcbkb1zSIrxRuptvFxdeTi3lw/3ys9+bgbTxPJtpptpHpSyfXlxjsUEwVLBKZzqSmFnZ85DbK7Vjp3+Fs+ZL8uQUnpkDMAiiaCLxsyC8CTAp4kSMmB3qMhdzBEem7WNwV5b105ry3axVurd1HtD5LidfOFob358ejefHFYb7JSD9GDyByTrdVbKc4s7tLtHInEEoHpOIEGpzG2Yr3zts3uTc7v+t37qnZqS/eV92VBz6Ew9hJAYMT5UDAGsgqPeKi99U28ubqM11aW8v76cvyBMNmpXr40Mp9zRvXm9CG9SEmyapxY2VazjSE9hsQ7DBMlSwQmNkJB2PQOrHsVVr/s9IDV0P5lPD7nQp8/GtzJ4PZCai7kjYTCCZAzEI7ijnLH3gbeWLWL11aW8snnuwmFlYIsHzMnFHHOqHwmDsix6p0OEAwHKakt4azis+IdiomSJQLTfj5/z3m/ftWLzoXfvxe8qVAwDnqPgn6nQGYfJwHkjYRew47qQn+ghqYQ68tqeH99Ba+tLOWzEufVncF56VzzhYGcMyqfE/pmWfVEB9tZt5NgOEi/zH7xDsVEyRKBOT4VG2D967DyX1Cy0FmXM9Cpvjn39zDsPPBlHtchgqEwmyvrWVtaw9pdNawtrWZtaQ1bdte3DDY5tiibm6YN45xR+QzqlX6cJ2WOx7ZqZ2C6ooyiI5Q0nYUlAnP0QgHnor/4CfjsaWddciYMnw7TfgvZx/bKoKpSWu1nTWkN60prWFtaw5rSGjaU19IUdIZlcAn075nGiIJMLhzfl+H5GYwr6kF+lq+9zs4cp5Y+BPbqaJdhicAcmr/a6WXrr3J+tn4CVVudKiAAccPEWTDum1Aw9qiqearqA/vu7nc5F/21pTVU+4MtZfIzfQzLz+C0IT0Z1juDYfkZDM5Lx+e1Rt7ObGvNVlI8KfRK6RXvUEyULBEYR8Mep1G3cj3s3QbbPnEmH2lLr+Ew4Tsw+iJI63nY3foDITaU1bK2tIZ1u5w7/LWlNZRW7xtFM8PnYXh+BueP7cPw/AyG5WcytHc62alJ7XmGpoNsrd5KUUaRtc10IZYIEtmGt2Dj27BjKWz5YN/6jAIomgT5/wV9ToT03uBNcer9Pclt7ioUVrburmdtabVTtRO56G+uqCMcqcdP8rgY3CudUwblMjTfucMfnp9BfqbPLhrdyNaarQzKGhTvMMxRsESQaFRhxfPw6d+c1zvB6YA19jKnemfirJaOWG1/XSmvaWy5s2+u1llfVoM/4NTji0C/nFSG5WcwfUyflmqd/rmp9vpmNxcKhyipKeGMojPiHYo5CpYIEkFjDcz/E2z+ECrWQV2ZMwTDiBlw9m2Q2/bdW7U/wPrInf26SMPt2l017K0PtJTplZHM8PwMLp/Uj6GRO/wheRnWWStB7arfRSAcsMHmuhhLBN1Zw1747Bl45abICnHq9QdMhXGXOR24cHrhri+rZf2uWtaX1bChrJYNZbXsrNpXj5+e7GFo73TOHV3A8PwMhkbu8nPSrB6/u1BVqpuqaQg2ENIQYQ23+VNWX+ZM1K7hlnIhDREOh9lYtRHA+hB0MZYIuqv5D8Drt0QWBB3/LSq/eDfry+rYUFbD+pfXRi78tVTUNrZ8LcXrZkjvdKYMymVwXnpLtU7f7BSrx+8i6gP1lNaX0hhsZI9/DzWBGoLhIMFwkB11O6hurKY2UMu2mm2s27OuZVsgHDjyzqPgdXkZmGWzknUllgi6k9Uvwefvw5InW8bV39j3Qu5Kvp6Fn+1l90dvtRTNSPYwuHc6XxzeiyF5GQzunc6QvHT6ZKV02lmvEpWqtlyo64P1VDRU4A/6qWqsYrd/NysrV7KzbieldaWs27Muqn3mp+WT7k3n5N4n0ye9D8nuZLxuLw2BBvqk9yHFk4Lb5UYQ3OLGJa6Wn5CG6JHcgzRvWss6t7hxuZzfWUlZZPuyY/uHYtqVJYKurn43rHwB3rwdGp0hFvamD2ZBaBg37rmIho0++mbXcPaIPIbnZzKkdzpD8jLonZlsd/hxpqrUBeooayijurGanXU72VG7g6XlS9lWvY3GUCP+kJ/d/t2EDzHPcWuTCiZx7oBzCYaDnFl0JqneVFI8KfRM6YnX5cXj8uB1eclOzsbnsQ54Zh9LBF2Rv9qp+y9dDsuehpBTtfMf1yRm1V9Loz+JopwUfvKV/nxxeB4DeqbZRb8DqSpbqrdQWl/KXv9eNlZtxB/0s2b3GjwuD3v8e9jj30Olv5LGUONB3++R3IPxeeNJ8aaQ5EoizZtGji8Hr8uL1+0lGA4yKHsQPreP/LR8MpIyyErOisOZmu7CEkFXEg7Bh/fBh390evoCVZlDuT94IY/tHkNeZgp/uHQkJ/TNojgn1ap4YkRV2bh3I29ufZNl5csIhAIEwgGC4SANoQa212ynPli/33dSPCl4xENtoJZT+p7CwKyB5KbkkuvLJcWTQqo3lb7pfclKzqJ/Zn/cLnvrynQcSwRdxYa34P17nY5ffU6k6ew7uW2hm79/WsmIgkzu/voAzh9bYFMltpOwhnlp40ssr1iOP+inIdjA3sa9LChdgFvchDSEIAzpMYQ0bxpelxefx0eOK4eJ+RMpyigiNyWX4oxi+mf2J9WbGu9TMuaQLBF0dqGA8xTw9q+c5S//mrLR3+Pqvy3m062V/OCsIfzwrCF293+U6gJ17PHvYbd/N2X1ZVQ3VbN+z3o+r/qc3f7dlNSUUBOoIdWTSmZyJj63j+zkbE7teyrp3nQm5k/kzKIz6ZVq4+mYrs8SQWdVsR4WPw4fPbBv3U/W8emeJK554ANq/EEevvxEpo0uiFuIXUFYwywsXciHOz7k1c9fJcmdRH2gnvKG8oPKCsKg7EHk+nL5Uv8vMaH3BL4y8Cu4xHpDm+7NEkFn07AX3vylkwQAsvvBpGvg5O/y3NIybnlhCb2zkvnXdacwPP/4xvnv6krrStnbuJfNVZtZvXs1gXCAZWXLCGmI6qZqKhoqaIi8RttsSsEU8lLz6JvRl/zUfHr4etA7tTfZydlkJWdZFY5JSJYIOotwCD57Fl6/FeorwJsGZ/8SJl1NMBTm1/NW838fbubUwbk8cOmJ9EiAHr0NwQb2+PewsnIlO2p3sLxieUvVTUVDRZvfSfOmMSh7EMN6DOP0vqeTmZxJdnI2ZxWfRV5qnt3dG9MGSwSdQTgEfzoR9mx2ls+/D066EoA9dU1c/48lfLihku+cOoD/d97wbjVwW1VjVcsFfln5MsrrywmEAzQEG9hcvZlgeN/8BHmpeYzMGcnQHkPp4etBZlImg7MH43F5GJ4znPy0/DieiTFdlyWCeNu9CV68wUkCfSfAFS9CsjPV4prSaq56chG7qhv5w9fHcvFJhfGN9ThsqtpEWX0ZgZBzkf+07FPe3/4+W6q3tJTJ9eUyMHsgKZ4UclNymVwwmaE9hpLkTmJc3jj6pveN4xkY031ZIoiX2jL4+0xnBjCA038CZ94CkffHX1m+k588t4z0ZA/PzJrM+OIecQz26L277V3m75jP1uqtbKvZ1jJ9YWsucTGt/zSm9Z/G4B6DKc4oto5vxsSBJYJ4qC2DR74IVdtg6Lkw8SoYfBYA4bDyxzfXcf/bGxhfnM3Dl59E78zONxxAVWMVC0sXOnf54QC7/bsprStlU9Um1uxe01KuX2Y/CjMKGd1zNF8f+nWS3El4XV4ykjIozOi6TzjGdCeWCDrajqXwzLecJHDOb2HKdS2bavwBfvTMMt5cvYtvTCjkVxeO7lQdxFSVeZ/P4+2tb/P6ltf32yYIvdN60y+zHzOHzcTn9nHVmKts6ANjuoCYJgIRmQbcB7iBv6jq7w7YXgw8AWRHytysqvNiGVNcvX8PvHWH8/n8++Gkb7ds2lxRx1VPLmJTRR23zxjFFVP6dZpqkqZQE0vLlvL02qd5Y8sbAFw05CKmFk6lOKOY/LR8fB4fHpfdVxjTFcXsX66IuIEHgS8BJcBCEZmrqqtaFbsVeFZV/1dERgLzgP6xiiluAg3w6s1O34Cew+BbL0DWvobP/6wr54a/L8HtEv76nYmcMvjwE8J3hIqGChbsXMDrW15n/o75NAQbEITLhl/GtWOvtWGGjelGYnkLNxHYoKqbAETkaeACoHUiUKC5V1QWsCOG8cRHOAwPTd73aug3n2tJAqrKX97/nN++spqhvTN45IoJFOXEr0NTSU0Jj698nF11u3i35F0AeqX0YsagGZze93TG9x5PZlJid2IzpjuKZSLoC2xrtVwCTDqgzG3A6yJyA5AGnN3WjkRkFjALoLi4i82FumC2kwQGnQXf+hfgJICPNlXyxzfXs+Dz3Zw9Io/7Lx1PalLHVq2U1JSwbs865n0+j493fkxVZD6DIT2GcEbhGZxRdAYXDL7AqnyM6ebi/S/8UuBxVb1HRKYAfxWR0ar7z8KhqrOB2QATJkzQOMR5bNa9Dq/+HPJGwTf/yZbKOl5dUcq/lmxn7a4aeqR6uf7Mwfzw7CEd0kns7a1vs6JiBXsa97C8fDlr96wFICMpg5PyTmJkz5GcN+A8m2/WmAQTy0SwHShqtVwYWdfad4FpAKr6kYj4gJ5AWQzj6hhzroOlT0FqLlz2NB9s3M01f1tMbWOQkQWZ3HXRGGaM64PPG9u3gppCTby2+TX+tvpvrKp0auUykzLpm96XWWNmMblgMmN6jSHZnRzTOIwxnVcsE8FCYIiIDMBJAJcAlx1QZitwFvC4iIwAfMDBw0J2NQsfdZKAxwdXv88Lm5SfPbeAwXnpzP7WBIpzY9sOUNFQwX1L7mPj3o1s2LuhZeC1oT2GMvtLs8lNyY3p8Y0xXUvMEoGqBkXkeuA1nFdDH1PVlSJyB7BIVecCPwEeEZEf4TQcX6mqXafqpy31u503hAB+voW3N1bxo2cWMa4omye/O5FMnzcmh20MNfLSxpf4+5q/s37P+pb13xzxTab2ncq4vHE2sqYxpk0xbSOI9AmYd8C6X7T6vAo4NZYxdKiaXfDa/4NQE+Fpv+et9VX88OlPGZ6fwRPfaf8kEAwH+bTsU+5bch8rKlYQ0hDJ7mTG9BrDzGEzmTFoRrsezxjTPcW7sbj7WPY0vHA1ABXDL+f8d4aws2oRRTkpPP5fE8lKab8koKq8uPFF7l9yf8sEK6NzR3P12Kv5QuEXOk1HNGNM12CJoD1serclCWyY8nu++tEA6psaue+ScZx3QgHednwj6NOyT7l30b0sLV/KyNyR3DzxZsb0GmNDMBtjjpklguMVDsG/nCRwz5i5PPhuLQN6JvPAZScyoqB9Ol8tLVvK/y77X9bsXsNu/24Afjrhp1w+4nLcrs4zFpExpmuyRHC8/vlfUFvKR8Wz+NOCWi4c14c7v3oC6cnH/0dbVl/G7xb8jje3vInitKGfO+Bcvj/u+/auvzGm3VgiOB6rXoRVLxLqMYiflU1jQr8U/njJ+OPe7dbqrTy24jGeX/88AKf1PY27pt5FRlLGce/bGGMOZIngWIUC8NotqDuZ6U2/oay2ibsuHntcu6xqrOKuhXcx7/N5BMNB8tPyuXXSrUwtnGoNwMaYmLFEcKze+Q1UbePdAT9h9eoQ/7hqMlMGHXtHrYZgA5f9+7KWmbzmXDCHQdmD2itaY4w5pO4zC3pHqquAD+6FwpN5LHgOIwoyjysJLCtfxqUvX8rWmq386KQfsfzbyy0JGGM6jCWCYzH/TwBUT72dJVv2cGJx9jHv6sPtH3LFK1ewo24HD3zxAb4z+jvtFKQxxkTHqoaOlr8KPvwjDJjKb1dk4g9WccWU/se0q2A4yE3v3URYwzwz/RkGZA1o11CNMSYaUT8RiIgNVAPw/FUAbC84m6cXbuXKU/ozLP/o3uapaqzihfUvcPrTp1PdVM0tk26xJGCMiZsjPhGIyCnAX4B0oFhExgJXq+p1h/9mN1RTChveQAecwffXnUjP9EZ+ePaQo9rF4l2LufLVKwFwiYtz+5/LzGEz2z9WY4yJUjRVQ/8DnAPMBVDVZSIyNaZRdVbzfgoaZl7Rj1j6ejX/M3MsGUcxkNz7Je9z3VtO/vz2yG9z40k34nXFZjRSY4yJVlRtBKq67YD32EOxCacT274EVr9EOLOIWz9o4uT+PbhwXN8jfy+iPlDPdW9dR6onlbkXzqV3Wu8YBmuMMdGLpo1gW6R6SEXEKyI/BVbHOK7OpWEPvHQjAHcXPUBVQ4DbZ4yOupOXqvKNl78BwO2n3G5JwBjTqUSTCK4Bvo8zGf12YByQWO0DL1wLpZ9RccIs/rykjium9Gdkn+gGlPt458ecP+d8tlRv4Ucn/YhpA6bFOFhjjDk60VQNDVPVb7ZeISKnAh/GJqROJtAA619DvWnM2nUhPVLr+dGXhh7xayU1Jdz8/s0sK18GwHdHf5crRl4R62iNMeaoRZMI/gScGMW67mn+A6Bhlo39JUs+2MtdF4054iQzf1n+F+5bch8Akwsm84spv6Aoo6gjojXGmKN2yEQgIlOAU4BeIvLjVpsyceYg7v7CIVj8OAA3rx/OwF7KRScVHrK4qvLg0gf582d/pl9mP353+u8Y3XN0BwVrjDHH5nBPBEk4fQc8QOseU9XAxbEMqtMoXQ7VJWwe+h3WfFbLPV8fi9vVdgOxqnLWc2dR3lBO3/S+PD/jeZLdyR0csDHGHL1DJgJV/Q/wHxF5XFW3dGBMncfWjwG4r2wsRTkpzBjX55BF71p4F+UN5fTL7MfcC+fiEhvGyRjTNUTTRlAvIncDowBf80pV/WLMouosVr9EoMcQXtjZk1u/0r/NuYdVlXsW3cPfVv+NYT2G8ez5z1oSMMZ0KdFcsZ4C1gADgNuBzcDCGMbUOdSWw9b5rOpxJiCcOrhnm8V+t+B3PLHqCQZlDeKJc5+wJGCM6XKiuWrlquqjQEBV/6Oq3wG6/9PAwkdAwzxcPpoT+ma1ORH9Hv8e5myYA8A/Z/yTNG9aBwdpjDHHL5pEEIj83ikiXxGR8UBODGPqHJb+HRU3r5TncsnEg1/9rGqs4vJ5l1MfrOe3p/8Wj8tG9DbGdE3RXL3uFJEs4Cc4/QcygR/GMqi427MZqrbxQc+ZpOzyMGPs/o3ENU01nPb0aQCcN+A8pg+cHocgjTGmfRwxEajqy5GPVcCZ0NKzuPta/RIAvy47heljCvYbYbS6qZpv/tvpaH39uOu5euzVcQnRGGPayyGrhkTELSKXishPRWR0ZN10EZkPPNBhEcbDov9jb8YQ1jT14pKJxS2rVZUb3rqBzdWb+WLRFy0JGGO6hcM9ETwKFAELgPtFZAcwAbhZVed0QGzxEWiA3RtZ5PsSQ/LS95uP+P3t77OkbAkn55/MfV+8L34xGmNMOzpcIpgAjFHVsIj4gFJgkKpWdkxocbL2FQD+r3oCl5xX3DLUdFjD/Ow/PwPgN6f9Jm7hGWNMezvcW0NNqhoGUFU/sOlok4CITBORtSKyQURuPkSZb4jIKhFZKSJ/P5r9x8Tm92l0p/GJjuKr4/dNPLOgdAH1wXquHXst+Wn5cQzQGGPa1+GeCIaLyGeRzwIMiiwLoKo65nA7FhE38CDwJaAEWCgic1V1VasyQ4D/Bk5V1T0ikncc59I+1r7Kes8QhuZnk5OWBDhPAze/dzOCcNnwy+IcoDHGtK/DJYIRx7nvicAGVd0EICJPAxcAq1qVuQp4UFX3AKhq2XEe8/jU7IKaHazXIUwY2aNl9erdq6n0V3Lx0IvJ9mXHLz5jjImBww06d7wDzfUFtrVaLgEmHVBmKICIfIgztPVtqvrqgTsSkVnALIDi4uIDN7efz54B4J+BU/hGv32J4NXPnZCuHmNvCRljup94D4zjAYYAZwCXAo+ISPaBhVR1tqpOUNUJvXr1il00u1YCsCg8jJMiiaCmqYbHVz5OQVqBtQ0YY7qlWCaC7TivnzYrjKxrrQSYq6oBVf0cWIeTGOJj8wdUevPpkZlJ3+wUAC6fdzkAP5/487iFZYwxsRRVIhCRFBEZdpT7XggMEZEBIpIEXALMPaDMHJynAUSkJ05V0aajPE77CDRA7S6WhQZwUv8eiAgLSxeyqWoT43qN46zis+ISljHGxNoRE4GInA8sBV6NLI8TkQMv6AdR1SBwPfAasBp4VlVXisgdIjIjUuw1oFJEVgHvAD+LWz+FFc9DOMD/+b/AhEi1UHPbwG9Ot34DxpjuK5pB527DeQPoXQBVXSoiA6LZuarOA+YdsO4XrT4r8OPIT3xtfBuAD8KjualfDmEN8/Kmlzmp90k28bwxpluLahhqVa06YJ3GIpi4WjWXNVmn4fN6GV6QweJdi6kP1nPBoAviHZkxxsRUNE8EK0XkMsAd6QD2A2B+bMPqYDWlEA7wWVMfxhVl43W7+GTnJwBMLpgc5+CMMSa2onkiuAFnvuJG4O84w1H/MIYxdbwyp4/bnOphjI8MMlfeUA5A77Te8YrKGGM6RDRPBMNV9RbgllgHEzcV6wFYF+rDlUXZAOys3cmwHsNsDmJjTLcXzVXuHhFZLSK/ap6XoNupWEeTJ4MKMlvmJm4MNZKZfPA8xcYY090cMRGo6pk4M5OVA38WkeUicmvMI+tIFeuo8PXDJUJ+lg9wEkGyOznOgRljTOxFVe+hqqWqej9wDU6fgl8c/htdTMV6troKKchKwet2/kgaQ4343L44B2aMMbEXTYeyESJym4gsx5m8fj7OcBHdg78aanby4d4cJg3I2bc66CfJnRTHwIwxpmNE01j8GPAMcI6q7ohxPB0v0lC8JpjPj6cObFndFGrC57EnAmNM93fERKCqUzoikLipWAdAqbeY4fkZLav9Ib+1ERhjEsIhE4GIPKuq34hUCbXuSRzVDGVdRsU6gnhw5fZvmZ8YrLHYGJM4DvdEcGPk9/SOCCRutsxnh+RR2DOrZZWqWiIwxiSMQzYWq+rOyMfrVHVL6x/guo4JL/a0fjf+kIv+uakt65rCTQDWRmCMSQjRvD76pTbWndvegcRFOIRUrmOjFtAvN61ltT/oByDJZW8NGWO6v8O1EVyLc+c/UEQ+a7UpA/gw1oF1iGrnJahV4X6cltPqiSBkTwTGmMRxuDaCvwOvAL8Fbm61vkZVd8c0qo5SWwbANs2jf89WTwQh54nA2giMMYngcIlAVXWziHz/wA0iktMtksGezwHY5O5PXsa+i35jsBGwRGCMSQxHeiKYDizGeX1UWm1TYGBbX+pStjg1XJIzYP9XR8OWCIwxieOQiUBVp0d+RzUtZZfUsAeAnj167Le65YnAY4nAGNP9RTPW0Kkikhb5fLmI3CsixbEPrQNsmc8SGUlu+v5vB1kbgTEmkUTz+uj/AvUiMhb4CbAR+GtMo+oIqmiggb0hHz3S9k8ELW8N2eijxpgEEE0iCKqqAhcAD6jqgzivkHZtjdVIYzULQkPJSbUnAmNM4oomEdSIyH8D3wL+LSIuwBvbsDpAXQUApZpz0BOBtREYYxJJNIlgJs7E9d9R1VKcuQjujmlUHSHSh6CSzIOeCBpD9taQMSZxRDNVZSnwFJAlItMBv6o+GfPIYq2uHIBKzTz4icASgTEmgUTz1tA3gAXA14FvAJ+IyMWxDizm6p2qoUrNJOcQicCGmDDGJIJoZii7BThZVcsARKQX8Cbwz1gGFnPlawHYS/ohE4ENOmeMSQTRtBG4mpNARGWU3+vcIlVDIVcymb7982Fj0JmLoHVvY2OM6a6ieSJ4VUReA/4RWZ4JzItdSB1FCIiXopzUgy74Nk2lMSaRRDNn8c9E5GvAaZFVs1X1hdiG1QFqSlnnGrzfhDTNmkJNlgiMMQnjcPMRDAH+AAwClgM/VdXtHRVYrGnJQtYHJ+03/HQzeyIwxiSSw9X1Pwa8DFyEMwLpn4525yIyTUTWisgGEbn5MOUuEhEVkQlHe4xj0lSHhBqpCKUxsI1E0BhstDeGjDEJ43BVQxmq+kjk81oRWXI0OxYRN/AgzlSXJcBCEZmrqqsOKJcB3Ah8cjT7Py41pQCs00LObysRhBpJctsbQ8aYxHC4ROATkfHsm4cgpfWyqh4pMUwENqjqJgAReRpnvKJVB5T7FfB74GdHGfuxi/Qq3qm5DDhEIrAB54wxieJwiWAncG+r5dJWywp88Qj77gtsa7VcAkxqXUBETgSKVPXfInLIRCAis4BZAMXF7TACdu0uAPa6etAnK+Wgzf6QnzTPwQnCGGO6o8NNTHNmLA8cGbzuXuDKI5VV1dnAbIAJEybocR888kSQ0qMAl+vgvgJNoSZyknOO+zDGGNMVxLJj2HagqNVyYWRdswxgNPCuiGwGJgNzO6TBuL6CMEKPnvltbvYH/TbyqDEmYcQyESwEhojIABFJAi4B5jZvVNUqVe2pqv1VtT/wMTBDVRfFMCbn2FXb2avp9OuV2eb2xlCjvT5qjEkYMUsEqhoErgdeA1YDz6rqShG5Q0RmxOq40Qju3soW7U3vzLYbhC0RGGMSyRF7Fosz/sI3gYGqekdkvuJ8VV1wpO+q6jwOGI5CVX9xiLJnRBVxO/CUfES5jiEvo+2LvSUCY0wiieaJ4CFgCnBpZLkGp39A16SKhAM0kHzoRGAdyowxCSSaRDBJVb8P+AFUdQ/QdXtb+asAKNcs8tqoGgqGgwQ1aB3KjDEJI5pEEIj0ElZomY8gHNOoYqmxGoC1WtTmE0FTqAnAOpQZYxJGNIngfuAFIE9Efg18APwmplHFUlMdAEF3KmnJBzeR+EN+wKapNMYkjmiGoX5KRBYDZ+EML3Ghqq6OeWSx0lgLQFJq26+ONj8RWCIwxiSKaN4aKgbqgZdar1PVrbEMLGZqdgLgS81oc7M/GHkisA5lxpgEEc0MZf/GaR8QwAcMANYCo2IYV8xoOIgAja5D9yEAayMwxiSOaKqGTmi9HBko7rqYRRRjTXt2kAwMGTysze0tE9fbW0PGmARx1D2LI8NPTzpiwU4qUFNGSIXUHr3b3G5PBMaYRBNNG8GPWy26gBOBHTGLKMZC1bsoJ5vMlLbbAKyNwBiTaKJpI2jdqhrEaTN4PjbhxF64fg97NZ0MX9unbm8NGWMSzWETQaQjWYaq/rSD4ok5bayhDh+ZKd42t5c3lAOWCIwxieOQbQQi4lHVEHBqB8YTe0211KmPzEM8ESwodcbSy/HZxDTGmMRwuCeCBTjtAUtFZC7wHFDXvFFV/xXj2GLC1VRLHbltPhGEwiEWlC7g/IHnk5WcFYfojDGm40XTRuADKnHmKG7uT6BAl0wE3kA1VVrcZhvBmt1rqGqs4tS+3eshyBhjDudwiSAv8sbQCvYlgGbHP29wnLhDfppcPpI97oO2zd8xH4DJBZM7OixjjImbwyUCN5DO/gmgWddMBOEw3lADIU9qm5s/3PEhI3JGkJuS28GBGWNM/BwuEexU1Ts6LJKOEKjDTYhG78EDztUF6lhWtoxvj/p2HAIzxpj4OVzP4raeBLq2yMijeNMO2rSwdCFBDXJKn1M6OChjjImvwyWCszosio7SWANAuufgeXU+3P4hKZ4UxuWN6+CgjDEmvg6ZCFR1d0cG0iEiw0c0pPY5aNNHOz/i5PyTbbA5Y0zCOepB57q0Bie3eXz7Vw2V1JSwpXqLVQsZYxJSQiUC9TvzFXuT939r6KOdHwEwpc+UDo/JGGPiLaESQaDJqRpyZfTab/387fMpSCtgQOaAeIRljDFxlVCJwF/nNBb7UtJb1gXDQT7Z+Qmn9DkFke73opQxxhxJQiWCRn89AKlp+9oIVlSsoCZQY9VCxpiElVCJwO9vACAtdV8bwfwd83GJy4aVMMYkrIRKBLJ3CwDpafuqhubvmM/o3NE22qgxJmElVCKoE6dKKCPNeSKobqpmecVyqxYyxiS0hEoEoaYG6jSZTJ8zF8GCnQsIa9iGnTbGJLRo5iM4ZiIyDbgPZyTTv6jq7w7Y/mPgezhzIZcD31HVLbGKJ7m2hCDulklpPtzxIenedEb3HB2rQxoTV4FAgJKSEvx+f7xDMR3E5/NRWFiI19v2dLxtiVkiiMx3/CDwJaAEWCgic1V1VatinwITVLVeRK4F7gJmxiqmxrCLLKlHPS5Ulfnb5zMxfyJeV/R/YMZ0JSUlJWRkZNC/f397PToBqCqVlZWUlJQwYED0/aJiWTU0EdigqptUtQl4GrigdQFVfUdV6yOLHwOFMYyHcCjIRgoREbbWbGVH3Q4bVsJ0a36/n9zcXEsCCUJEyM3NPeonwFgmgr7AtlbLJZF1h/Jd4JW2NojILBFZJCKLysvLjz2iYAMBlw9wRhsFOKWvJQLTvVkSSCzH8vfdKRqLReRyYAJwd1vbVXW2qk5Q1Qm9evVqq0hU3MEGAq4UAD7a8RFFGUUUZRQd8/6MMaY7iGUi2A60vsoWRtbtR0TOBm4BZqhqYwzjISncQJM7BVVl0a5F1onMmA6Qnp5+5ELH6eGHH+bJJ5+M+XFamzNnDqtWrTpywVZUlR/84AcMHjyYMWPGsGTJkjbLLV68mBNOOIHBgwfzgx/8AFVnduDnnnuOUaNG4XK5WLRo0XGfQ7NYJoKFwBARGSAiScAlwNzWBURkPPBnnCRQFsNYAEgOO08E9cF6agO1FGcUx/qQxph2EgqFDrntmmuu4YorrujQYx5LInjllVdYv34969evZ/bs2Vx77bVtlrv22mt55JFHWsq++uqrAIwePZp//etfTJ069aiOeyQxe2tIVYMicj3wGs7ro4+p6koRuQNYpKpzcaqC0oHnIvVaW1V1Rqxi8moTIXcyFQ0VADZJvUkot7+0klU7qtt1nyP7ZPLL80dFXf7uu+/m2WefpbGxka9+9avcfvvtAFx44YVs27YNv9/PjTfeyKxZswDnaeLqq6/mzTff5MEHH2TatGnceOONvPzyy6SkpPDiiy/Su3dvbrvtNtLT0/npT3/KGWecwaRJk3jnnXfYu3cvjz76KKeffjr19fVceeWVrFixgmHDhrFjxw4efPBBJkyYsF+M/fv3Z+bMmbzxxhvcdNNN1NTUMHv2bJqamhg8eDB//etfWbp0KXPnzuU///kPd955J88//zwA3//+9ykvLyc1NZVHHnmE4cOH77fvF198kSuuuAIRYfLkyezdu5edO3dSUFDQUmbnzp1UV1czebJTY3HFFVcwZ84czj33XEaMGHH0f0lRiGk/AlWdB8w7YN0vWn0+O5bHP5Bbg+DytiSCnik9O/LwxiS0119/nfXr17NgwQJUlRkzZvDee+8xdepUHnvsMXJycmhoaODkk0/moosuIjc3l7q6OiZNmsQ999wDQF1dHZMnT+bXv/41N910E4888gi33nrrQccKBoMsWLCAefPmcfvtt/Pmm2/y0EMP0aNHD1atWsWKFSsYN27cIWPNzc1tqbaprKzkqquuAuDWW2/l0Ucf5YYbbmDGjBlMnz6diy++GICzzjqLhx9+mCFDhvDJJ59w3XXX8fbbb++33+3bt1NUtK/GvLCwkO3bt++XCLZv305hYeFBZWIppomgs3FrENweSwQmIR3NnXssvP7667z++uuMHz8egNraWtavX8/UqVO5//77eeGFFwDYtm0b69evJzc3F7fbzUUXXdSyj6SkJKZPnw7ASSedxBtvvNHmsb72ta+1lNm8eTMAH3zwATfeeCPgVLGMGTPmkLHOnLmvO9OKFSu49dZb2bt3L7W1tZxzzjkHla+trWX+/Pl8/etfb1nX2BjTJs92lVCJwIM9ERgTL6rKf//3f3P11Vfvt/7dd9/lzTff5KOPPiI1NZUzzjij5T14n8+H2+1uKev1eltej3S73QSDwTaPlZycfMQyh5PWaqj6K6+8kjlz5jB27Fgef/xx3n333YPKh8NhsrOzWbp06WH327dvX7Zt2/dWfUlJCX379j2oTElJyWHLtLdO8fpoR8mgHtxJVDZU4hGPjThqTAc655xzeOyxx6itrQWcKpCysjKqqqro0aMHqamprFmzho8//jgmxz/11FN59tlnAVi1ahXLly+P6ns1NTUUFBQQCAR46qmnWtZnZGRQU+NMdpWZmcmAAQN47rnnACfpLVu27KB9zZgxgyeffBJV5eOPPyYrK2u/aiGAgoICMjMz+fjjj1FVnnzySS644IKD9tWeEicRRF6/StdqKhoqyEnJwSWJc/rGxNuXv/xlLrvsMqZMmcIJJ5zAxRdfTE1NDdOmTSMYDDJixAhuvvnmlkbS9nbddddRXl7OyJEjufXWWxk1ahRZWUe+GfzVr37FpEmTOPXUU/dr/L3kkku4++67GT9+PBs3buSpp57i0UcfZezYsYwaNYoXX3zxoH2dd955DBw4kMGDB3PVVVfx0EMPtWxr3Wbx0EMP8b3vfY/BgwczaNAgzj33XABeeOEFCgsL+eijj/jKV77SZjXVsZDm91O7igkTJuixvD+rwUbkzjw+KL6OvxdXUemv5Jnpz8QgQmM6j9WrV8fsTZOuJhQKEQgE8Pl8bNy4kbPPPpu1a9eSlJQU79DaXVt/7yKyWFUntFU+YdoIAo0NJAF4nNdHe6Ueew9lY0zXU19fz5lnnkkgEEBVeeihh7plEjgWCZMIQo3O2Hbi9lDZUMmIXLtLMiaRZGRktGtv3O4kYSrJQwHnVS5XOEClv5Jcn3UmM8YYSKBEEA42AbAnOZWQhuzVUWOMiUiYRBCKJIJqCQDWh8AYY5olTCIItyQC57clAmOMcSRMIggFnSeBGnHaCiwRGNMxbBjqfaIdhvqWW26hqKioQ/7sIIESgTY1AFCllgiM6YoSaRjq888/nwULFhzVvo9H4rw+Gnb+QmuoJ8WTQqo3Nc4RGdPBXrkZSqMbViFq+SfAub+LurgNQ33kYaiBmPWuPpTEeSIIOwNPVYfr7GnAmDhoPQz10qVLWbx4Me+99x4Ajz32GIsXL2bRokXcf//9VFZWArQMQ71s2TJOO+20lmGoly1bxtSpU3nkkUfaPFbzMNR//OMfW5JN62Gof/WrX7F48eJDxto8DPUll1zC1772NRYuXMiyZcsYMWIEjz76KKeccgozZszg7rvvZunSpQwaNIhZs2bxpz/9icWLF/OHP/yB66677qD9HmoY6nhLnCeCyCNedbjWEoFJTEdx5x4LNgx155UwiSAcCgPOE0FRyqA4R2NM4rFhqKMbhjoeEqZqKBxy/meoCtVYr2Jj4sCGoY5uGOp4SJxEEA7RBNSF622uYmPiwIahjn4Y6ptuuonCwkLq6+spLCzktttuO65zP5KEGYZ63X+eIfO9a/lScV9+OeWXXDz04hhEZ0znYsNQ72PDUNsw1ITDQSoidY3WWGxM4rFhqA8tYRKBhkOWCIxJYDYM9aElTBuBq6GSCo9zupYIjDFmn4RJBAF3WssTQY4vJ87RGGNM55EwiSAcDlPhdpPuSSfJbfWCxhjTLGESgYaCVLrd9PBmxzsUY4zpVBImEYTDYSrdLnokZcc7FGMSig1Dvc+aNWuYMmUKycnJ/OEPf4hRZEcvYd4aQp3XRwclW/uAMV1RKBTab7iJ1q655poOP+acOXOYPn06I0eOjHp/OTk53H///cyZM6edImwfCZMImtsITrZEYBLU7xf8njW717TrPofnDOfnE38edflEH4Y6Ly+PvLw8/v3vfx/PH3u7S5iqoYagnwaXi1x7Y8iYuLBhqDuvhHki2BuuB7BxhkzCOpo791iwYag7r5gmAhGZBtwHuIG/qOrvDtieDDwJnARUAjNVdXMsYqkK1wHQyzqTGRMXNgx15xWzqiERcQMPAucCI4FLReTAVpXvAntUdTDwP8DvYxXP3rAzZ3HP1F6xOoQx5jBsGOrOK5ZtBBOBDaq6SVWbgKeBCw4ocwHwROTzP4GzpDndt7MqdZ4I8lLticCYeLBhqKG0tJTCwkLuvfde7rzzTgoLC6murm7X8zwWMRuGWkQuBqap6vciy98CJqnq9a3KrIiUKYksb4yUqThgX7OAWQDFxcUnbdmy5ajj+cvcX/J+6as8/K23SEmJ/XvNxnQGNgz1PjYMdRcfhlpVZwOzwZmP4Fj28b0Zt/M9bm/XuIwxXYcNQ31osUwE24GiVsuFkXVtlSkREQ+QhdNobIwx7cqGoT60WLYRLASGiMgAEUkCLgHmHlBmLvDtyOeLgbe1q02ZZkwnZ/+kEsux/H3HLBGoahC4HngNWA08q6orReQOEZkRKfYokCsiG4AfAzfHKh5jEpHP56OystKSQYJQVSorK/H5fEf1vYSZs9iYRBQIBCgpKWl5L990fz6fj8LCQrxe737ru3xjsTHm2Hi9XgYMGBDvMEwnlzBjDRljjGmbJQJjjElwlgiMMSbBdbnGYhEpB46+a7GjJ1BxxFLdi51zYrBzTgzHc879VLXNwda6XCI4HiKy6FCt5t2VnXNisHNODLE6Z6saMsaYBGeJwBhjElyiJYLZ8Q4gDuycE4Odc2KIyTknVBuBMcaYgyXaE4ExxpgDWCIwxpgE1y0TgYhME5G1IrJBRA4a0VREkkXkmcj2T0SkfxzCbFdRnPOPRWSViHwmIm+JSL94xNmejnTOrcpdJCIqIl3+VcNozllEvhH5u14pIn/v6BjbWxT/bxeLyDsi8mnk/+/z4hFnexGRx0SkLDKDY1vbRUTuj/x5fCYiJx73QVW1W/0AbmAjMBBIApYBIw8ocx3wcOTzJcAz8Y67A875TCA18vnaRDjnSLkM4D3gY2BCvOPugL/nIcCnQI/Icl684+6Ac54NXBv5PBLYHO+4j/OcpwInAisOsf084BVAgMnAJ8d7zO74RDAR2KCqm1S1CXgauOCAMhcAT0Q+/xM4S0SkA2Nsb0c8Z1V9R1XrI4sf48wY15VF8/cM8Cvg90B3GIc5mnO+CnhQVfcAqGpZB8fY3qI5ZwUyI5+zgB0dGF+7U9X3gN2HKXIB8KQ6PgayRaTgeI7ZHRNBX2Bbq+WSyLo2y6gzgU4VkNsh0cVGNOfc2ndx7ii6siOec+SRuUhV/92RgcVQNH/PQ4GhIvKhiHwsItM6LLrYiOacbwMuF5ESYB5wQ8eEFjdH++/9iGw+ggQjIpcDE4AvxDuWWBIRF3AvcGWcQ+loHpzqoTNwnvreE5ETVHVvPIOKsUuBx1X1HhGZAvxVREarajjegXUV3fGJYDtQ1Gq5MLKuzTIi4sF5nKzskOhiI5pzRkTOBm4BZqhqYwfFFitHOucMYDTwrohsxqlLndvFG4yj+XsuAeaqakBVPwfW4SSGriqac/4u8CyAqn4E+HAGZ+uuovr3fjS6YyJYCAwRkQEikoTTGDz3gDJzgW9HPl8MvK2RVpgu6ojnLCLjgT/jJIGuXm8MRzhnVa1S1Z6q2l9V++O0i8xQ1a48z2k0/2/PwXkaQER64lQVberAGNtbNOe8FTgLQERG4CSC8g6NsmPNBa6IvD00GahS1Z3Hs8NuVzWkqkERuR54DeeNg8dUdaWI3AEsUtW5wKM4j48bcBplLolfxMcvynO+G0gHnou0i29V1RlxC/o4RXnO3UqU5/wa8GURWQWEgJ+papd92o3ynH8CPCIiP8JpOL6yK9/Yicg/cJJ5z0i7xy8BL4CqPozTDnIesAGoB/7ruI/Zhf+8jDHGtIPuWDVkjDHmKFgiMMaYBGeJwBhjEpwlAmOMSXCWCIwxJsFZIjCdkoiERGRpq5/+hylb2w7He1xEPo8ca0mkh+rR7uMvIjIy8vn/HbBt/vHGGNlP85/LChF5SUSyj1B+XFcfjdPEnr0+ajolEalV1fT2LnuYfTwOvKyq/xSRLwN/UNUxx7G/447pSPsVkSeAdar668OUvxJn1NXr2zsW033YE4HpEkQkPTKPwhIRWS4iB400KiIFIvJeqzvm0yPrvywiH0W++5yIHOkC/R4wOPLdH0f2tUJEfhhZlyYi/xaRZZH1MyPr3xWRCSLyOyAlEsdTkW21kd9Pi8hXWsX8uIhcLCJuEblbRBZGxpi/Ooo/lo+IDDYmIhMj5/ipiMwXkWGRnrh3ADMjscyMxP6YiCyIlG1rxFaTaOI99rb92E9bPzi9YpdGfl7A6QWfGdnWE6dXZfMTbW3k90+AWyKf3TjjDfXEubCnRdb/HPhFG8d7HLg48vnrwCfAScByIA2nV/ZKYDxwEfBIq+9mRX6/S2TOg+aYWpVpjvGrwBORz0k4o0imALOAWyPrk4FFwIA24qxtdX7PAdMiy5mAJ/L5bOD5yOcrgQdaff83wOWRz9k4YxGlxfvv237i+9Pthpgw3UaDqo5rXhARL/AbEZkKhHHuhHsDpa2+sxB4LFJ2jqouFZEv4ExW8mFkaI0knDvpttwtIrfijFPzXZzxa15Q1bpIDP8CTgdeBe4Rkd/jVCe9fxTn9Qpwn4gkA9OA91S1IVIdNUZELo6Uy8IZLO7zA76fIiJLI+e/GnijVfknRGQIzjAL3kMc/8vADBH5aWTZBxRH9mUSlCUC01V8E+gFnKSqAXFGFPW1LqCq70USxVeAx0XkXmAP8IaqXhrFMX6mqv9sXhCRs9oqpKrrxJnr4DzgThF5S1XviOYkVNUvIu8C5wAzcSZaAWe2qRtU9bUj7KJBVceJSCrO+DvfB+7HmYDnHVX9aqRh/d1DfF+Ai1R1bTTxmsRgbQSmq8gCyiJJ4EzgoDmXxZmHeZeqPgL8BWe6v4+BU0Wkuc4/TUSGRnnM94ELRSRVRNJwqnXeF5E+QL2q/g1nML+25owNRJ5M2vIMzkBhzU8X4FzUr23+jogMjRyzTerMNvcD4Ceybyj15qGIr2xVtAaniqzZa8ANEnk8EmdUWpPgLBGYruIpYIKILAeuANa0UeYMYJmIfIpzt32fqpbjXBj/ISKf4VQLDY/mgKq6BKftYAFOm8FfVPVT4ARgQaSK5pfAnW18fTbwWXNj8QFex5kY6E11pl8EJ3GtApaIM2n5nznCE3skls9wJma5C/ht5Nxbf+8dYGRzYzHOk4M3EtvKyLJJcPb6qDHGJDh7IjDGmARnicAYYxKcJQJjjElwlgiMMSbBWSIwxpgEZ4nAGGMSnCUCY4xJcP8fNcAD+XbnK2QAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from sklearn.ensemble import GradientBoostingClassifier # BoostType\n", "from sklearn.metrics import roc_curve, auc\n", "\n", "for learning_rate_i in [0.01,0.1,1]:\n", " bdt = GradientBoostingClassifier(learning_rate=learning_rate_i)\n", " bdt.fit(X_train, y_train) # fit BDT to training set\n", "\n", " decisions = bdt.decision_function(X_test).ravel() # get probabilities on test set\n", "\n", " # Compute ROC curve and area under the curve\n", " fpr, tpr, _ = roc_curve(y_test, # actual\n", " decisions ) # predicted\n", "\n", " plt.plot(fpr, tpr, label='learning rate '+str(learning_rate_i)) # plot test ROC curve\n", "\n", " plt.xlabel('False Positive Rate') # x-axis label\n", " plt.ylabel('True Positive Rate') # y-axis label\n", " plt.legend() # add legend" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Notice how the ROC curve changes with each different learning rate." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Putting everything into a machine learning algorithm means we only have 1 variable to optimise. The signal and background distributions are separated much better when looking at machine learning output, compared to individual variables. Using machine learning algorithms also achieves much higher S/B values than on individual variables.\n", "\n", "machine learning algorithm can achieve better S/B ratios because they find correlations in many dimensions that will give better signal/background classification.\n", "\n", "Hopefully you've enjoyed this discussion on using machine learning algorithms to select for signal to background." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Going further\n", "\n", "If you want to go further, there are a number of things you could try: \n", "\n", "* **Modify some BDT hyper-parameters** in '[Training your machine learning algorithm](#MVA_training)'. Cell -> Run All Below. You may find the [sklearn documentation on GradientBoostingClassifier](https://scikit-learn.org/stable/modules/generated/sklearn.ensemble.GradientBoostingClassifier.html) helpful.\n", "* **Try other machine learning algorithms** in '[Training your machine learning algorithm](#MVA_training)'. Cell -> Run All Below. You may find [sklearn documentation on supervised learning](https://scikit-learn.org/stable/supervised_learning.html#supervised-learning) helpful.\n", "\n", "With each change, keep an eye on the ROC curve." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] } ], "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.3" } }, "nbformat": 4, "nbformat_minor": 4 }