{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "
" ] }, { "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", "ATLAS Open Data provides open access to proton-proton collision data at the LHC for educational purposes. ATLAS Open Data resources are ideal for high-school, undergraduate and postgraduate students.\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", "The idea is that cuts 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, the amount of $Z$ and $t\\bar{t}$ background is reduced, since these are quite different to the signal.\n", "\n", "Then, the amount of $ZZ \\rightarrow \\ell\\ell\\ell\\ell$ is reduced, whilst keeping as much $H \\rightarrow ZZ \\rightarrow \\ell\\ell\\ell\\ell$ signal as possible.\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.\n", "\n", "This analysis loosely follows the [discovery of the Higgs boson by ATLAS](https://www.sciencedirect.com/science/article/pii/S037026931200857X) (mostly Section 4 and 4.1)\n", "\n", "By the end of this notebook you will be able to:\n", "1. rediscover the Higgs boson yourself!\n", "2. know some general principles of a particle physics analysis\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": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Contents: \n", "\n", "[Running a Jupyter notebook](#running)
\n", "[First time setup on your computer (no need on mybinder)](#setup_computer)
\n", "[To setup everytime](#setup_everytime)
\n", "[Lumi, fraction, file path](#fraction)
\n", "[Samples](#samples)
\n", "[Changing a cut](#changing_cut)
\n", "[Applying a cut](#applying_cut)
\n", "[Plotting](#plotting)
\n", "[What can you do to explore this analysis?](#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 clicking Cell -> Run Cells, or using the keyboard shortcut Shift+Enter." ] }, { "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", "# update the pip package installer\n", "!{sys.executable} -m pip install --upgrade --user pip\n", "# install required packages\n", "!{sys.executable} -m pip install --upgrade --user uproot awkward vector numpy matplotlib" ] }, { "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", "* uproot: lets us read .root files typically used in particle physics into data formats used in python\n", "* awkward: lets us store data as awkward arrays, a format that generalizes numpy to nested data with possibly variable length lists\n", "* vector: to allow vectorized 4-momentum calculations\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": {}, "outputs": [], "source": [ "import uproot # for reading .root files\n", "import awkward as ak # to represent nested data in columnar format\n", "import vector # for 4-momentum calculations\n", "import time # to measure time to analyse\n", "import math # for mathematical functions such as square root\n", "import numpy as np # for numerical calculations such as histogramming\n", "import matplotlib.pyplot as plt # for plotting\n", "from matplotlib.ticker import AutoMinorLocator # for minor ticks\n", "\n", "import infofile # local file containing cross-sections, sums of weights, dataset IDs" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Lumi, fraction, file path\n", "\n", "General definitions of fraction of data used, where to access the input files" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "#lumi = 0.5 # fb-1 # data_A only\n", "#lumi = 1.9 # fb-1 # data_B only\n", "#lumi = 2.9 # fb-1 # data_C only\n", "#lumi = 4.7 # fb-1 # data_D only\n", "lumi = 10 # fb-1 # data_A,data_B,data_C,data_D\n", "\n", "fraction = 1.0 # reduce this is if you want the code to run quicker\n", " \n", "#tuple_path = \"Input/4lep/\" # local \n", "tuple_path = \"https://atlas-opendata.web.cern.ch/atlas-opendata/samples/2020/4lep/\" # web address" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Samples\n", "\n", "samples to process" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "samples = {\n", "\n", " 'data': {\n", " 'list' : ['data_A','data_B','data_C','data_D'],\n", " },\n", "\n", " r'Background $Z,t\\bar{t}$' : { # Z + ttbar\n", " 'list' : ['Zee','Zmumu','ttbar_lep'],\n", " 'color' : \"#6b59d3\" # purple\n", " },\n", "\n", " r'Background $ZZ^*$' : { # ZZ\n", " 'list' : ['llll'],\n", " 'color' : \"#ff0000\" # red\n", " },\n", "\n", " r'Signal ($m_H$ = 125 GeV)' : { # H -> ZZ -> llll\n", " 'list' : ['ggH125_ZZ4lep','VBFH125_ZZ4lep','WH125_ZZ4lep','ZH125_ZZ4lep'],\n", " 'color' : \"#00cdff\" # light blue\n", " },\n", "\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Units, as stored in the data files" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "MeV = 0.001\n", "GeV = 1.0" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "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 at least 4 leptons per event, so that processing is quicker." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def get_data_from_files():\n", "\n", " data = {} # define empty dictionary to hold awkward arrays\n", " for s in samples: # loop over samples\n", " print('Processing '+s+' samples') # print which sample\n", " frames = [] # define empty list to hold data\n", " for val in samples[s]['list']: # loop over each file\n", " if s == 'data': prefix = \"Data/\" # Data prefix\n", " else: # MC prefix\n", " prefix = \"MC/mc_\"+str(infofile.infos[val][\"DSID\"])+\".\"\n", " fileString = tuple_path+prefix+val+\".4lep.root\" # file name to open\n", " temp = read_file(fileString,val) # call the function read_file defined below\n", " frames.append(temp) # append array returned from read_file to list of awkward arrays\n", " data[s] = ak.concatenate(frames) # dictionary entry is concatenated awkward arrays\n", " \n", " return data # return dictionary of awkward arrays" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "define function to calculate weight of MC event" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def calc_weight(xsec_weight, events):\n", " return (\n", " xsec_weight\n", " * events.mcWeight\n", " * events.scaleFactor_PILEUP\n", " * events.scaleFactor_ELE\n", " * events.scaleFactor_MUON \n", " * events.scaleFactor_LepTRIGGER\n", " )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "define function to get cross-section weight" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "def get_xsec_weight(sample):\n", " info = infofile.infos[sample] # open infofile\n", " xsec_weight = (lumi*1000*info[\"xsec\"])/(info[\"sumw\"]*info[\"red_eff\"]) #*1000 to go from fb-1 to pb-1\n", " return xsec_weight # return cross-section weight" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "define function to calculate 4-lepton invariant mass.\n", "\n", "Note: `lep_(pt|eta|phi|E)` are variable length lists of lepton momentum components for each event, represented by awkward arrays." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "def calc_mllll(lep_pt, lep_eta, lep_phi, lep_E):\n", " # construct awkward 4-vector array\n", " p4 = vector.zip({\"pt\": lep_pt, \"eta\": lep_eta, \"phi\": lep_phi, \"E\": lep_E})\n", " # calculate invariant mass of first 4 leptons\n", " # [:, i] selects the i-th lepton in each event\n", " # .M calculates the invariant mass\n", " return (p4[:, 0] + p4[:, 1] + p4[:, 2] + p4[:, 3]).M * MeV" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Changing a cut\n", "\n", "If you change a cut: Cell -> Run All Below\n", "\n", "If you change a cut here, you also need to make sure the cut is applied in the \"[Applying a cut](#applying_cut)\" cell." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# cut on lepton charge\n", "# paper: \"selecting two pairs of isolated leptons, each of which is comprised of two leptons with the same flavour and opposite charge\"\n", "def cut_lep_charge(lep_charge):\n", "# throw away when sum of lepton charges is not equal to 0\n", "# first lepton in each event 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", "# paper: \"selecting two pairs of isolated leptons, each of which is comprised of two leptons with the same flavour and opposite charge\"\n", "def cut_lep_type(lep_type):\n", "# for an electron lep_type is 11\n", "# for a muon lep_type is 13\n", "# throw away when none of eeee, mumumumu, eemumu\n", " sum_lep_type = lep_type[:, 0] + lep_type[:, 1] + lep_type[:, 2] + lep_type[:, 3]\n", " return (sum_lep_type != 44) & (sum_lep_type != 48) & (sum_lep_type != 52)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Applying a cut\n", "If you add a cut: Cell -> Run All Below" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def read_file(path,sample):\n", " start = time.time() # start the clock\n", " print(\"\\tProcessing: \"+sample) # print which sample is being processed\n", " data_all = [] # define empty list to hold all data for this sample\n", " \n", " # open the tree called mini using a context manager (will automatically close files/resources)\n", " with uproot.open(path + \":mini\") as tree:\n", " numevents = tree.num_entries # number of events\n", " if 'data' not in sample: xsec_weight = get_xsec_weight(sample) # get cross-section weight\n", " for data in tree.iterate(['lep_pt','lep_eta','lep_phi',\n", " 'lep_E','lep_charge','lep_type', \n", " # add more variables here if you make cuts on them \n", " 'mcWeight','scaleFactor_PILEUP',\n", " 'scaleFactor_ELE','scaleFactor_MUON',\n", " 'scaleFactor_LepTRIGGER'], # variables to calculate Monte Carlo weight\n", " library=\"ak\", # choose output type as awkward array\n", " entry_stop=numevents*fraction): # process up to numevents*fraction\n", "\n", " nIn = len(data) # number of events in this batch\n", "\n", " if 'data' not in sample: # only do this for Monte Carlo simulation files\n", " # multiply all Monte Carlo weights and scale factors together to give total weight\n", " data['totalWeight'] = calc_weight(xsec_weight, data)\n", "\n", " # cut on lepton charge using the function cut_lep_charge defined above\n", " data = data[~cut_lep_charge(data.lep_charge)]\n", "\n", " # cut on lepton type using the function cut_lep_type defined above\n", " data = data[~cut_lep_type(data.lep_type)]\n", "\n", " # calculation of 4-lepton invariant mass using the function calc_mllll defined above\n", " data['mllll'] = calc_mllll(data.lep_pt, data.lep_eta, data.lep_phi, data.lep_E)\n", "\n", " # array contents can be printed at any stage like this\n", " #print(data)\n", "\n", " # array column can be printed at any stage like this\n", " #print(data['lep_pt'])\n", "\n", " # multiple array columns can be printed at any stage like this\n", " #print(data[['lep_pt','lep_eta']])\n", "\n", " nOut = len(data) # number of events passing cuts in this batch\n", " data_all.append(data) # append array from this batch\n", " elapsed = time.time() - start # time taken to process\n", " print(\"\\t\\t nIn: \"+str(nIn)+\",\\t nOut: \\t\"+str(nOut)+\"\\t in \"+str(round(elapsed,1))+\"s\") # events before and after\n", " \n", " return ak.concatenate(data_all) # return array containing events passing all cuts" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is where the processing happens" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Processing data samples\n", "\tProcessing: data_A\n", "\t\t nIn: 39,\t nOut: \t27\t in 1.7s\n", "\tProcessing: data_B\n", "\t\t nIn: 156,\t nOut: \t86\t in 1.3s\n", "\tProcessing: data_C\n", "\t\t nIn: 237,\t nOut: \t146\t in 1.3s\n", "\tProcessing: data_D\n", "\t\t nIn: 400,\t nOut: \t248\t in 1.6s\n", "Processing Background $Z,t\\bar{t}$ samples\n", "\tProcessing: Zee\n", "\t\t nIn: 898,\t nOut: \t243\t in 1.2s\n", "\tProcessing: Zmumu\n", "\t\t nIn: 684,\t nOut: \t257\t in 1.0s\n", "\tProcessing: ttbar_lep\n", "\t\t nIn: 1031,\t nOut: \t334\t in 1.1s\n", "Processing Background $ZZ^*$ samples\n", "\tProcessing: llll\n", "\t\t nIn: 554279,\t nOut: \t523957\t in 6.9s\n", "Processing Signal ($m_H$ = 125 GeV) samples\n", "\tProcessing: ggH125_ZZ4lep\n", "\t\t nIn: 164716,\t nOut: \t161451\t in 4.9s\n", "\tProcessing: VBFH125_ZZ4lep\n", "\t\t nIn: 191126,\t nOut: \t186870\t in 6.1s\n", "\tProcessing: WH125_ZZ4lep\n", "\t\t nIn: 15379,\t nOut: \t9772\t in 2.7s\n", "\tProcessing: ZH125_ZZ4lep\n", "\t\t nIn: 14485,\t nOut: \t11947\t in 1.4s\n", "Time taken: 31.3s\n" ] } ], "source": [ "start = time.time() # time at start of whole processing\n", "data = get_data_from_files() # process all files\n", "elapsed = time.time() - start # time after whole processing\n", "print(\"Time taken: \"+str(round(elapsed,1))+\"s\") # print total time taken to process every file" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting\n", "If you only want to make a change in plotting: Cell -> Run All Below\n", "\n", "Define function to plot the data" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def plot_data(data):\n", "\n", " xmin = 80 * GeV\n", " xmax = 250 * GeV\n", " step_size = 5 * GeV\n", "\n", " bin_edges = np.arange(start=xmin, # The interval includes this value\n", " stop=xmax+step_size, # The interval doesn't include this value\n", " step=step_size ) # Spacing between values\n", " bin_centres = np.arange(start=xmin+step_size/2, # The interval includes this value\n", " stop=xmax+step_size/2, # The interval doesn't include this value\n", " step=step_size ) # Spacing between values\n", "\n", " data_x,_ = np.histogram(ak.to_numpy(data['data']['mllll']), \n", " bins=bin_edges ) # histogram the data\n", " data_x_errors = np.sqrt( data_x ) # statistical error on the data\n", "\n", " signal_x = ak.to_numpy(data[r'Signal ($m_H$ = 125 GeV)']['mllll']) # histogram the signal\n", " signal_weights = ak.to_numpy(data[r'Signal ($m_H$ = 125 GeV)'].totalWeight) # get the weights of the signal events\n", " signal_color = samples[r'Signal ($m_H$ = 125 GeV)']['color'] # get the colour for the signal bar\n", "\n", " mc_x = [] # define list to hold the Monte Carlo histogram entries\n", " mc_weights = [] # define list to hold the Monte Carlo weights\n", " mc_colors = [] # define list to hold the colors of the Monte Carlo bars\n", " mc_labels = [] # define list to hold the legend labels of the Monte Carlo bars\n", "\n", " for s in samples: # loop over samples\n", " if s not in ['data', r'Signal ($m_H$ = 125 GeV)']: # if not data nor signal\n", " mc_x.append( ak.to_numpy(data[s]['mllll']) ) # append to the list of Monte Carlo histogram entries\n", " mc_weights.append( ak.to_numpy(data[s].totalWeight) ) # append to the list of Monte Carlo weights\n", " mc_colors.append( samples[s]['color'] ) # append to the list of Monte Carlo bar colors\n", " mc_labels.append( s ) # append to the list of Monte Carlo legend labels\n", " \n", "\n", "\n", " # *************\n", " # Main plot \n", " # *************\n", " main_axes = plt.gca() # get current axes\n", " \n", " # plot the data points\n", " main_axes.errorbar(x=bin_centres, y=data_x, yerr=data_x_errors,\n", " fmt='ko', # 'k' means black and 'o' is for circles \n", " label='Data') \n", " \n", " # plot the Monte Carlo bars\n", " mc_heights = main_axes.hist(mc_x, bins=bin_edges, \n", " weights=mc_weights, stacked=True, \n", " color=mc_colors, label=mc_labels )\n", " \n", " mc_x_tot = mc_heights[0][-1] # stacked background MC y-axis value\n", " \n", " # calculate MC statistical uncertainty: sqrt(sum w^2)\n", " mc_x_err = np.sqrt(np.histogram(np.hstack(mc_x), bins=bin_edges, weights=np.hstack(mc_weights)**2)[0])\n", " \n", " # plot the signal bar\n", " main_axes.hist(signal_x, bins=bin_edges, bottom=mc_x_tot, \n", " weights=signal_weights, color=signal_color,\n", " label=r'Signal ($m_H$ = 125 GeV)')\n", " \n", " # plot the statistical uncertainty\n", " main_axes.bar(bin_centres, # x\n", " 2*mc_x_err, # heights\n", " alpha=0.5, # half transparency\n", " bottom=mc_x_tot-mc_x_err, color='none', \n", " hatch=\"////\", width=step_size, label='Stat. Unc.' )\n", "\n", " # set the x-limit of the main axes\n", " main_axes.set_xlim( left=xmin, right=xmax ) \n", " \n", " # separation of x axis minor ticks\n", " main_axes.xaxis.set_minor_locator( AutoMinorLocator() ) \n", " \n", " # set the axis tick parameters for the main axes\n", " main_axes.tick_params(which='both', # ticks on both x and y axes\n", " direction='in', # Put ticks inside and outside the axes\n", " top=True, # draw ticks on the top axis\n", " right=True ) # draw ticks on right axis\n", " \n", " # x-axis label\n", " main_axes.set_xlabel(r'4-lepton invariant mass $\\mathrm{m_{4l}}$ [GeV]',\n", " fontsize=13, x=1, horizontalalignment='right' )\n", " \n", " # write y-axis label for main axes\n", " main_axes.set_ylabel('Events / '+str(step_size)+' GeV',\n", " y=1, horizontalalignment='right') \n", " \n", " # set y-axis limits for main axes\n", " main_axes.set_ylim( bottom=0, top=np.amax(data_x)*1.6 )\n", " \n", " # add minor ticks on y-axis for main axes\n", " main_axes.yaxis.set_minor_locator( AutoMinorLocator() ) \n", "\n", " # Add text 'ATLAS Open Data' on plot\n", " plt.text(0.05, # x\n", " 0.93, # y\n", " 'ATLAS Open Data', # text\n", " transform=main_axes.transAxes, # coordinate system used is that of main_axes\n", " fontsize=13 ) \n", " \n", " # Add text 'for education' on plot\n", " plt.text(0.05, # x\n", " 0.88, # y\n", " 'for education', # text\n", " transform=main_axes.transAxes, # coordinate system used is that of main_axes\n", " style='italic',\n", " fontsize=8 ) \n", " \n", " # Add energy and luminosity\n", " lumi_used = str(lumi*fraction) # luminosity to write on the plot\n", " plt.text(0.05, # x\n", " 0.82, # y\n", " '$\\sqrt{s}$=13 TeV,$\\int$L dt = '+lumi_used+' fb$^{-1}$', # text\n", " transform=main_axes.transAxes ) # coordinate system used is that of main_axes\n", " \n", " # Add a label for the analysis carried out\n", " plt.text(0.05, # x\n", " 0.76, # y\n", " r'$H \\rightarrow ZZ^* \\rightarrow 4\\ell$', # text \n", " transform=main_axes.transAxes ) # coordinate system used is that of main_axes\n", "\n", " # draw the legend\n", " main_axes.legend( frameon=False ) # no box around the legend\n", " \n", " return" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Call the function to plot the data" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEKCAYAAADzQPVvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABUMklEQVR4nO2deXgUxda43yKRJQQJEHYkkYtsAQKEVZAkIggCAgIqBgUFERXF5S568bviT+N2UdSrfhLcUKIgIipwXTHDByqCQQIERATCjrIFErJAkvP7o2eGmWTWZCYzSep9nn6mu7q6+nRP9+nq06fOUSKCRqPRaKo3tQItgEaj0Wj8j1b2Go1GUwPQyl6j0WhqAFrZazQaTQ1AK3uNRqOpAWhlr9FoNDWA0EAL4Ijw8HDp1KmT23pnzpyhYcOGlV7v+PHjNG3a1CftBbNsvq4XzLIFSr5gli1Q8gWzbL6u52vZ0tPTz4qI44oi4nACdgCPAX9xVsdfU506deTOO++Uzz//XFxx5513ulzvr3pxcXE+ay+YZfN1vWCWTSQw8gWzbN7U0/dE+er5SrbPP/9c7rzzTgHOiBO96sqMMwmoD3ytlNqolHpQKdXK7aPFB4SEhJCSksLo0aNd1nO33l/1PMWT9oJZNn/UC8Q+g/ncBbNs3tTzFH1PlA9PdGFKSgrAn04rOXsK2E5Af2A+cABIA+70ZLvyTmFhYR497QKFp0/jQKBlKz/BLF8wyyYS3PLVJNmAn6UcPXvbB8IGEXkQuA2IAF71ZLvyEhkZ6c/mK8yMGTMCLYJTtGzlJ5jlC2bZILjl07IZKHETG0cp1QfDpDMe2AcsAZaJyEm3jSsVAbwJdAUEuAPYBSwFooEs4EYROW27Xe/eveXnn3/27kg0Go2mhqOUSheR3o7WOe3ZK6WeVkrtAV4HDgMDRSRBRN7wRNGbeRn4UkQ6AbHATuARYI2IXAGsMS/7leTkZJRSLFq0CIADBw4QHh5unUJDQ6ldu7Z1OSYmBoDo6GgWL17ssu3vv/8epRS33357mXXHjx9n2rRptG7dmvDwcFq2bMmIESM4evSo0/ZKSkp48cUXiYmJISwsjCZNmjB+/Hh27NhRgTNQcbKyslBKUb9+fRo0aECjRo3o3bs3c+fO5cyZMx63YzKZCA0NSicwjaZa48qMUwAMF5E+IvKCiBzypmGlVENgMPAWgIicF5FsYAywyFxtETDWW6G9oaSkhIULF9K4cWPLBwzatm1Lbm6udUpISOCf//yndTkzM9Pj9hcsWEDjxo356KOPyii9yZMnk5OTwy+//EJubi4ZGRlMmjQJpZTT9m6//XZefPFF5s+fT3Z2Ntu3b6d58+b069ePrVu3lu8k+JBdu3aRk5PDn3/+ySuvvMKaNWvo3bs3J096+vzXaDQBwZkx3zIBYcD/AAvNy1cAozzYrgewEXgX+AXDnFMfyLapo2yXLVPbtm0lLi7OOi1YsKDcHyz++9//SmhoqKxatUoA2bZtW5k6Q4YMkccff7xMeVRUlLz//vtO2z516pTUrVtXUlNTpUmTJvKf//zHbn14eLisXr3aY1nXrVsngJhMpjLrEhISZMiQIdZlQObPny+xsbESHh4uCQkJsnv3buv6CxcuSHJyslxxxRXSsGFDufLKK2XTpk3W9VOmTJHJkyfL9OnTpWHDhtKqVSt54403nMq2b98+AeTgwYN25SdPnpSmTZvKP//5TxEROXfunIwbN06aN28uDRo0kJ49e8rXX38tIiKHDx+WunXrCiD169eX+vXry7vvvisiIlOnTpU2bdpIeHi4dO7cWVJTUz0+bxpNTWbBggVWXQlkiTOd7GyFXFTIS4G/A9vlovLf4sF2vYEioJ95+WXgydLKHThdeltffqEeO3asjB49WkREunfvLrNmzSpTp7zK/qWXXpLIyEgpLCyU+++/X7p162a3/rrrrpMuXbrIggULZPPmzVJUVORS1kcffVTatGnjcN2bb74pISEhkpeXJyKGsu/cubPs3r1b8vLy5N5775XOnTtb9/HPf/5T+vbtK3v27JGioiJ58803pUmTJnLq1CkRMZR93bp15bPPPpPi4mJZvny5hIaGSlZWlsP9O1P2IiK33HKL9OvXT0REcnJy5P3335ezZ8/K+fPn5fnnn5cGDRrIn3/+KSIiaWlpEhIS4vD4Tpw4IUVFRfLhhx/KJZdcIpmZmS7Pl0bjTwYNGiSxsbFB7c1TGlx443ii7H82//5iU5bhwXYtbJ8ywFXAaowPtC3NZS2BXaW39dXJPXz4sISGhsqKFStEROTll1+WiIgIq8K0UF5l36VLF5k9e7aIiGRkZAggP/zwg3V9Tk6OPP3009K/f3+pW7euNGzYUGbPni35+fkO25s+fbpVaZbmv//9rwBy6NAhETGU/Ztvvmldf+7cOaldu7Z8//33UlJSIuHh4bJ27Vq7Nrp27Wo9nilTpsh1111ntz4yMlI+/fRTh/t3pez//ve/S/v27R1uJyLSpEkT6xuOM2Vfmri4OHnttdfc1tNoNBdxpew9cb08r5Sqh+FNg1LqL0Chu41E5BhwUCnV0Vw0BGNU7ufAFHPZFOAzD2QoF2+99RaNGzdm1KhRgGFDz8/PZ+nSpRVue926dezYsYM77rgDgO7du9O7d28WLFhgrRMeHs6jjz7Kjz/+yJkzZ3jvvfd45513ePrppx222bRpUw4fPuxw3ZEjRwgJCaFx48bWsujoaOt8WFgYTZs25dChQ5w4cYLc3FxGjx5NRESEddq7dy+HDl389NKyZUu7fdSvX5+cnByvz8WhQ4do0qQJAPn5+cyaNYt27dpx6aWXEhERwenTpzl+/LjT7UtKSvjXv/5Fx44dadiwIREREWRkZLjcRqPReIcnyv5x4EvgMqVUKoYHzd89bP8+IFUptRXDhv808CwwVCm1G7jGvOxzSkpKeOutt8jOzqZNmza0aNGCLl26UFxcbKeQy4vlY++wYcNo0aIFLVq0YMeOHXz00UdkZ2eXqV+7dm2uv/56rrnmGrZs2eKwzeHDh3Po0CHWrVtXZt0HH3xAfHw89erVs5ZlZWVZ5/Py8jh+/Dht2rQhMjKS+vXr8+2335KdnW2dzp07xyOP+Nb56fTp03zzzTdcffXVALz44ov83//9H2vWrOHMmTNkZ2fTqFEjy9sdtWqVveQ+/PBD3nzzTZYvX87p06fJzs4mNjbWuo2mYoSEhNCjRw9iY2Pp1asXP/zwg9dtZGVl0bVrVz9IV3Hmzp3LvHnzHK67/fbb6dGjh3Vq0aKFXYepNIcOHbJ2Bm3nqwNulb2IfAPcAEwFPgR6i4jJk8ZFZIuI9BaR7iIyVkROi8hJERkiIleIyDUicqoiB+CML7/8koMHD/LDDz+wZcsW67Rq1So2bNjAtm3bPGrnwoULFBQUWKfCwkJOnTrFxx9/zGuvvWbX9o4dO6hbty7vv/8+AA899BCbNm2ioKCAkpISTCYTaWlpXHXVVQ73NXjwYG655RaSkpL49ttvOX/+PMeOHeO+++7jp59+KnNBz58/nz179lBQUMAjjzxCu3bt6NevH0opZs+ezV//+ld2794NQG5uLl999RVHjhypwFm9SFFRERs2bGDcuHE0aNCAhx56CICzZ89Sp04dmjRpwvnz5/l//+//2T38WrRoQXFxMfv27bOWnT17ltDQUJo2bUpJSQlvv/02GRkZPpFTA/Xq1WPLli1kZGTwzDPP8Oijj1bq/kWEkpKSSt2nhXfeecd6f65YsYLQ0FDeffddp/XXrFnD5s2by8xXB5w6PCulQoB6IpIrIieVUieB2kAXpdQvIuL9+34lsmDBAsaOHUtcXJxdeYsWLRgwYAALFizg1VfdDwS+4447rKYagDp16vDMM8/QqFEjpk+fTu3ate3qz5w5kwULFnDfffdRUlLC7bffzoEDB1BK0bp1a/7617/y8MMPO93fe++9x/z587n//vvJysqibt26xMfHs2HDhjI9q+nTp3PDDTewd+9eevXqxWeffUZISAgATzzxBK+88gpjxozh0KFD1K9fn/79+/Of//zH7TG7omPHjtSqVYuQkBDatWvHyJEjefjhh4mIiACMB9zmzZtp1aoVERERPPDAA3bmpg4dOnD33XfTt29fLly4wH/+8x+mTJnCd999R/v27QkLC+PWW291+kCsysy6ZbtP23v1A+972mfPnqVRo0YAjB07loMHD1JQUMDs2bOtoznfe+895s2bh1KK7t27WzsvFvbu3cv48eNJSUmhT58+PPnkkyxevJimTZty2WWXERcXx4QJE7j22mvp168f6enp/Pe//2X58uW8/fbbgHHtPvDAA2RlZTFq1Ci2bzfOzbx588jNzWXu3LlkZWUxYsQIBg0axA8//EDr1q357LPPqFevHsnJySxatIhmzZpZ9+mKEydOMHz4cP7nf/6H66+/3mGd9evX89BDDxEREcELL7xAgwYNaNy4MV999RWffPIJ7dq18/p8BxNOR9AqpeYBf4rI8+blvcB2oB6wWUT+4S+hrrjiCklMTGT06NE+D3ZUXVBKsW7dOgYNGhRoUTQeEihlHxISQrdu3SgoKODo0aN89913xMXFcerUKRo3bkx+fj59+vRh7dq1HDt2jHHjxvHDDz8QGRlprWNRysuXL+fmm2/m3XffJTY2lk2bNnHnnXeyYcMGLly4QK9evbjrrruYMGEC7dq144cffqB///6kp6czdepUNmzYgIjQr18/Fi9eTKNGjVwq+/bt2/Pzzz/To0cPbrzxRq6//no6d+7M1KlT+emnnygqKqJXr17MnDmTv/71rw6PPy8vjyFDhnDNNdfw5JNPujxXw4cPZ968eXTt2tVuPthZuXIlK1euZOHChb+LMWC1DK6GMg4B+tgsnxGR65UxIqisUdmHNGzY0GoT12g0FcNixgH48ccfue2229i+fTuvvPIKK1asAODgwYPs3r2bTZs2MXHiRGt8Klv79vHjxxkzZgyffPIJXbp0AYwR5GPGjKFu3brUrVvXrnMWFRVF//79AaPXPG7cOOrXrw/ADTfcwLp165z2si1cfvnl9OjRA4C4uDiysrI4ceIE48aNIywsDMBlG8XFxdx888106tTJraIHY9CgJZeG7XywY+kYL1y40Olwdlc2+1oiUmSz/A8As3tPuI9k1Gg0lciAAQM4ceIEy5Yt49tvv+XHH38kIyODnj17UlBQ4HLbhg0b0rZtW9avX+/RviyK3RWhoaF29vzSMtSpU8c6HxISQlFREd5wzz33cOHCBRYuXOi27okTJ2jYsCGhoaF289UFV8q+tlKqgWVBRL4GaxiEuv4WTOMaEdEmHI3X/PrrrxQXF1OnTh0aNWpEWFgYv/76Kxs2bADg6quvZtmyZdbwF6dOXfSfqF27NitWrOC9997jgw8+AGDgwIGsXLmSgoICcnNzWbVqlcP9XnXVVXz66afk5eVx7tw5VqxYwVVXXUXz5s35888/OXnyJIWFhU63t2Xw4MF8+umn5Ofnk5OTw8qVKx3We+KJJ0hPT2fZsmVllPaQIUPKuDlnZWXRqlWrMvPVBVfKfiGwVCnV1lKglIrC8Mh509+C+YLff/+dAQMGMGvWLJ+1+emnn1bYxPThhx9SXFwMwPPPP6/jymj8Sn5+vtX18KabbmLRokUMHz6coqIiOnfuzCOPPGI1t8TExDBnzhzi4+OJjY21ellZqF+/PqtWrWL+/Pl8/vnn9OnTh+uvv57u3bszYsQIunXr5jB9Xq9evZg6dSp9+/alX79+TJ8+nZ49e3LJJZfwr3/9i759+zJ06FCPzCa9evXipptuIjY2lhEjRtCnT58ydbKyspg7dy4nT55k0KBBdsdfUlLC77//XsYFs1OnTpw4cYKuXbuSl5dnnS+Pq2ow4jLEsVJqJvBPjJg2ALnAsyLyv/4Uylchjh9//HFGjBhhvZCdUVJS4tD/2xHPPPMMgwcPZuDAgeWWa+jQoXzzzTfl3l5TNQkGbxx/kJubS3h4OHl5eQwePJiUlBR69eoVaLGcsn37dt5++21efPHFQIvic1yFOPYocxTQAGjgSV1fTL4Il7BmzRpp3LixxMfHy5YtW2T27NkyYMAAufbaayU3N1fS09Nl2LBhMnbsWLsQAYcOHZLrr79eBg0aZA2hsHnzZhk4cKAMGzZMRo4cKadPn5a//e1vkpGRISIid999t+zZs0fOnTsnt956qyQkJMj48eNFROTBBx+0xtg4fPiw/Oc//5EmTZpIfHy87Ny5U4YPHy4iIl988YX069dP+vbtK999952IiAwePFjuv/9+6dy5syxevLjC50Sj8QeTJk2S2NhY6dixozz99NOBFqdGQ0Vi4wRiat++vUcJx91x7bXXiojIu+++K88++6yIGPFxPvzwQ1m0aJFMnz69zDY33nij/PbbbyIiMnr0aCksLJRrr71WTpw4ITk5OXL55ZeLiMioUaOkoKBARESGDx8uJSUl8tBDD8mqVatEROT8+fMiIpKbmysiIm+//bZ88MEHsmfPHrn33ntFRGT//v1y5513ypkzZyQhIUHy8vIkOztbhg8fLoWFhRIVFSWnT5+WvXv3ym233Vahc6HRaKovNgnHd4sTvRqUn5p94XpZWFhI3brGd+SVK1daB3PUqlWLevXqsXnz5jIpwc6fP8/atWu58847ATh58iQHDhygQ4cONGnShIKCAqtNMS8vjzp16lhHByqlyMzM5IUXXgDgkksu4fDhwzz88MP88ccfHDp0iDfffJOtW7fSvXt3wHid7Nq1K//3f//HDTfcQL169SguLqZu3br89ttvjBw5koiICHbv3s3ll19eofOh0WiqLxV1vazS7Ny5k86dOwOGr+2FCxcoLCxk+fLlxMfHs2vXLqvStZCbm8uAAQMwmUyYTCY2b95Mdna29WNqSkoKHTt25Pz589aRqmvXrqV169YAZGdnW+O5FBUV8eSTT3L//ffzzTff0KJFC7p168bOnTvp2NGIDWdR9kVFRZw/f966j5EjR5KZmWn1L87IyLBmz9JoNJry4LJnr5TqhJFZqrW56DDwuYjs9LdgFWXbtm1069YNgNmzZzNixAhCQkJ49NFHiYiIoLCw0M6HF4wBJO3atWPgwIHUq1ePWbNmMWrUKHbu3MnVV19NUVERt956K7Vr16Z27do88MAD5OTkWEfYjR8/nr59+1KvXj0WLVrElVdeyT333MPAgQO5cOECjRs3pnv37tx333089NBDbN++nalTp9KgQQPGjRvHJ598Qvfu3Xn11Vd54oknGDlypPVYZs6cWbknUKPRVCtchUv4B0ai8SWAJS5uG+BmYImI+CVaJeiE4xqNRlMeXHnjuOrZTwNiRORCqcZeBDLxU2hijUaj0fgeVzb7EsDRELKW5nWaIObAgQN07dqVzp0788cffwRaHI1GE2BcKfsHgDVKqS+UUinm6UuM5CWz/SnUmTNnmDFjhtNh0Br3vPLKK8TGxpKQkMBnn/ktGVi1Ze/evUybNo0JEyYEWhSNxi0rV660eBeWHb5sxt0I2lpAX+w/0G4SkWIfylkGf9vsjcCd7nF1boKdCRMmEBcXV+mJKqobEyZM4OOPPw60GBqNR7iy2bt0vRSREhHZICLLzdMGfyt6f3PrrbdiMpk8HTnslDvuuINmzZrZxbouKCigb9++xMbGEhMTw+OPP15mu5MnT9qlSGvdurV12eJ+WZrExES++uoru7KXXnqJu+++26l8OTk5hIeXDU7qqMwZtunesrOzef311z3e1hWOzp2FL7/8ko4dO9K+fXuefdbxZyFP6oDxdtO5c2eSkpKcxiTftm0bo0aNspv+/PPP8h2YRhPMeKL0HCjBVeXZztPJF+ESHHHkyBF58sknfdLW2rVrJT09XWJiYqxlJSUlkpOTIyLGCNq+ffvKjz/+6LSNxx9/XP7973+73deCBQtk6tSpdmX9+vWTtWvXOt3myiuvlLfffrtMef369d3uz5F8+/btszvWiuDo3ImIFBUVSbt27WTPnj1SWFgo3bt3l8zMTK/rWOjYsaMcPHiwQrJbwl5oahaZmZnyzjvvyIEDB+Ts2bOBFsdjcBEuobyDqu703eOm8li6dClJSUnW5UWLFhEXF0f37t29Dhc8ePDgMlHzlFLWnvOFCxe4cOGCxyYjC4sXL6Zv37706NGDu+66i+LiYiZMmMDq1autPf+srCyOHDniMnWfJTiVtyQnJ9OhQwcGDRrErl27rOWPPPIIe/bsoUePHvztb3/zul1bHJ07gI0bN9K+fXvatWtH7dq1ufnmm8t8b/CkDhjpIffu3cuIESNYsWIFRUVFJCUl0blzZyZMmEBeXp5LGU+ePMnMmTP55ZdfeOaZZyp0vIFGJxwvm3DcXSJyS8rMFStWlOs+CkbKFS5BRI76WhB/cPLkSX7++WeuvfZaAPbv328NO5CTk8Nzzz3Hli1bqF27tl1SbDDib+fklE2zO2/ePK655hqn+ywuLiYuLo7ff/+de++9l379+nks786dO1m6dCnff/89l1xyCffccw+pqancdttt9O3bly+++IIxY8awZMkSbrzxRpcPkvIo+/T0dJYsWcKWLVus6d4suT2fffZZtm/fbs14VJryni9bDh8+zGWXXWZdbtOmDT/99JPXdQDeeOMNvvzyS9LS0sjNzeWhhx7irbfeYuDAgdxxxx28/vrrTtPYATRp0oQ33njDI7mDHdtMVV999RWPPvooa9eurbT9W3uWHkaW9SXvvPOOdX7fvn1cddVVpKSk2GW3si23cPDgQW6//XbatWtHTk4Ol156aaXK7Q+cnn2l1HCb+YZKqbeUUluVUh8opZpXjngVIysri7/+9a9kZ2ezc+dO64haMHo7+fn5PPzww/z888/WhNkW1q1bZ81Kbzu5U1whISFs2bKFQ4cOsXHjRmtuTU9Ys2YN6enp9OnThx49erBmzRr27t0LwKRJk1iyZAkAS5YsYdKkSS7bys3NpUGDBi7rlGbdunXWdG+XXnqp25Rxpbctz/mqLC677DJrWOrJkyd7nG3Jpyjl26kclE44HhcXR0xMjJ2ie++99+jevTuxsbHceuutZdrYu3cvPXv2ZNOmTQA8+eSTdOzYkUGDBjFp0iTmzZtHVlYWHTt25LbbbqNr164cPHiQF198ka5du9K1a1deeukloOwbw7x585g7d651XefOnbnzzjuJiYlh2LBh5OfnA87fQJ3hLOG4s/JRo0YxYcIErrvuumqh6MF1z/5p4Evz/AvAUWA0cAOwABjrL6EsrpcVTTgeFxfHxIkTWbZsGX/88QezZ1/0GA0LC2P79u1Wl6Xp06dzzz33WNdXtKcaERFBYmIiX375pcevvyLClClTHJoNxowZw4MPPsjmzZvJy8uz9ridUV4zTnnxRc++devWHDx40Lp86NAha9whb+o4ovRbkLfmtaqMJXmJbcJxgLffftsu4fj48eM5duwYTz31lF3CcVt27dpVJuH48uXLycjIsCYct1ybu3fvZtGiRdaE4++88w4//fQTIkbC8fj4eOuDxxm7d+/mww8/ZOHChdx4440sX76czp07O30DdUReXh6jR4/mxhtv5K677nJbbqFFixYen+NAY0k4jgvXS0/fq3qLyGMisl9E5gPRPpDPKZaolxVR9BYmT57M+++/T15enl1Pd/fu3dSvX5+bb76ZUaNGlcl9WZ6e6vHjx63moPz8fL755huvEhYPGTKEjz/+2OoNcurUKfbv3w8YXjSJiYnccccdZXr1pVOslZSUkJeX57Wyd5XurUGDBg6VuQVf9Oz79OnD7t272bdvH+fPn2fJkiVl3i48qeOIAwcO8OOPPwLwwQcf1KiUjhYzzq+//sqXX37JbbfdhohYx2L079/fmnD8u+++c5twPDU1ldjYWMA+4XiDBg08SjgeHh5uTTjuDkcJx715A3WWcNzbROTBzujRoy1vZ+WKetlMKfWQUuph4FJl3xWqMtEy27VrR3FxMQkJCXblycnJdOzYkV69erFv3z67Xr0nTJo0iQEDBrBr1y7atGnDW2+9xdGjR0lMTKR79+706dOHoUOHMmrUKI/b7NKlC0899RTDhg2je/fuDB06lKNHL34emTRpEhkZGXbK3lGKtdzcXMCxm2VeXh5t2rSxTrbZelyle2vSpAkDBw6ka9euFf5A6+jcgZF8+tVXX+Xaa6+lc+fO3HjjjdZon9dddx1HjhxxWccVHTt25LXXXqNz586cPn3apdtqdUYnHHddXq1x5qYDPF5qamoubwG852y7Um1kAduALZhdgoDGwDfAbvNvo9Lb+dr1ct26dVJcXOzTNoOFbdu2yYMPPmhXdvjwYQGsbqCaIAF8O3mIrbvtzp07pUmTJvLpp5/KqFGjrGV16tSRtLQ02b59u1xxxRVy4sQJERE5efKkiFx0vc3NzZWBAwdKamqqiIhs3LhRevbsKfn5+ZKTkyNXXHGF/Pvf/y7j7pqeni7dunWTc+fOSW5ursTExMjmzZvl/Pnz0qRJEzlx4oQUFBRIv379rBniSrfx73//Wx5//HFrW3l5eXL27Flp3769QxfmuXPnSlxcXJn7wFl5dQAXrpdObfYi8oST8mPAbV48TxJF5ITN8iPAGhF5Vin1iHn5H1605zXV+ZW9a9eu1t75yZMn6devH8899xyRkZHVxmVMUzEsNnswOneLFi3immuu4Y033qBz58507NjRYcLxkJAQevbsybvvvmtty5JwfOjQoYSHh3P99ddbE443b97co4TjgDXhOGBNON66dWuvE443a9bMZcLx6Ohou/u/Y8eOfPTRRw7Lly5d6v5kVmFchkuocONKZWHY+0/YlO0CEkTkqFKqJWASkY622+kQx+Xn7rvvZsmSJcyfP5+pU6cGWhyNLb7+KBwk4TyqWsLx6oyrcAn+Vvb7gNOAAAtEJEUplS0iEeb1CjhtWbZQHZX9jh072LhxI0OGDCEiIsJrt0iNJli55ZZb2LFjBwUFBUyZMkXHYwog5Ypnr5RqJSJHKrjvQSJyWCnVDPhGKfWr7UoREaVUmafN8ePH6d37orwzZswoky/WUxYsWMCWLVv43//9X2tZ165dWbZsmTVtoafcfvvt/PLLL9blY8eOcf78eU6dOuVyHVwckXf27Fnuu+++ch2LLcXFxfTu3ZvWrVuzatUqAM6dO8c//vEPJk+ebH0t12j8zQcffBBoEWo0KSkptuMkIp3Vc+Vn/6ZSqjFgwvC3Xy8iXn0KF5HD5t8/lVIrMCJo/qGUamljxikTdapp06b4qme/bds2u1fKgoICsrKy6NChg9dtORuN524d+H5E3ssvv0znzp05e/asteyNN96goKCA9evXa2Wv0dQQbDvDSqkTzuo5daEUkeuABAxlPw7YoJT6RCk1QynV1p0ASqn6SqkGlnlgGLAd+ByYYq42BfBrsPWtW7faKftt27bRoUMHa8JwZ9tYRq46wtmoO1frvB2R50qGQ4cOsXr1aqZPn25XbokGafkYp9FoNBbchTguEJEvRWS22Q70MMbbwKtKqY1u2m4OrFdKZQAbgdUi8iVGOsOhSqndwDX4Ob1hZmYmN9xwA9HR0URHRzNixAi6d+/ucpuCggLGjh3rUNm6GnXnyxF5rmR44IEHeP755+1ijRQUFFBcXMzmzZuJj4/3eD8ajaZm4FUgNBHZB7wOvK6Uqu2m7l4g1kH5SWCIN/stLwcPHqRp06b8+uvFTwWzZs2yBkMDI8qko5joR48e5eabb2bjxovPNFej7ioyIs8bGVatWkWzZs2Ii4vDZDJZy3fv3k1xcTGdOnXikksu8Wr/Go2mBuDMAT+Qk68GVa1evVpuuOEGu7LExET5+uuvXW63f/9+iY2NlfXr19uVz5gxQ4YPHy4XLlwos42rdeXBmQyPPPKItG7dWqKioqR58+ZSr149SUpKkjVr1kiLFi3k9OnTPtm/RqOpeuCHePZVgq1bt9KlSxe7sszMTLvol47YtWsXr7/+ujVKIsATTzxBeno6y5YtIzTU/oXI1bry4kgGgGeeeYZDhw6RlZXFkiVLuPrqq1m8eDFHjhxh/PjxlJSUcPr0aZ/IoNFoqg9eaSalVCPgMhHZ6id5AN9FvbSknLNw6tQpRMSt7Xzo0KF2y65G4z333HNO11VkRF5pGVxRVFTE5s2bOXbsGDNnzuTNN98s9341Gk3Vw5Ool24HVSmlTMD1GA+GdAxXye9F5CGfSVqK6jioSqPRaPxNuROOm2koImcx4ti/JyL9MLxoNBqNRlNF8ETZh5oHP90IrPKzPBqNxsckJycTExND9+7d6dGjhzWN45VXXunzfTkLvpefn098fDzFxcU+36c33HHHHTRr1qxMQqGDBw+SmJhIly5diImJ4eWXX7aui46Oplu3bvTo0cNuZH9p/vjjD2655RbatWtHXFwcAwYMYMWKFS7lSUxM5KuvvrIre+mll7j77rsZPHiw1yGdXeGJzf4J4CuMEbSblFLtMMITazQaL1CbfdueeBBr7Mcff2TVqlVs3ryZOnXqcOLECWvi+vIkHi8vb7/9NjfccIPLwYyVwdSpU5k1axa33WYfuDc0NJQXXniBXr16kZOTQ1xcHEOHDrU6eKSlpVkTujhCRBg7dixTpkyxho/Yv38/n3/+uUt5LOlGLXmywUg7+vzzz5OWlsbSpUtJSkoq7+Ha4UnP/qiIdBeRe8DqP/+im200Gk0QcPToUSIjI61JQCIjI2nVqhVg3wt3lkfWWQ5YZ/lrnZGamsqYMWOsyxMnTmTWrFkMGjSIqKgo1q9fz6233kqHDh2YNm2aL0+BHYMHD7ZL9GOhZcuW1pH2DRo0oHPnznbZ39zx3XffUbt2bWbOnGkti4qKssbBWrx4MX379qVHjx7cdddd1jecCRMmsHr1ausDOCsriyNHjnDVVVcxduxYUlNTy32spfFE2f/HwzKNRhNkDBs2jIMHD9KhQwfuuece1q5dW6aObR7ZL774wi4u1e7du7n33nvJzMwkIiKC5cuXA0ZPPT09nZ9//plXXnmFkydPOpXh/Pnz7N27l+joaGvZtm3baNeuHevXr+euu+5i2rRpPP/88+zYsYPVq1dTWFjo8TFeddVV9OjRo8z07bffetyGLVlZWfzyyy/069cPMPIVDxs2jLi4OKcPtszMTKdhnXfu3MnSpUv5/vvv2bJlCyEhIVYl3rhxY/r27csXX3wBGL36G2+8EaUUXbt2tSZ19wWuol4OAK4EmiqlbD1vLgX8+i7mK9dLjaamEx4eTnp6OuvWrSMtLY2bbrqJZ5991i7XgW0e2bp169rdc45ywAK88sorVnu0JX9tkyZNHMpw4sQJIiIirMsFBQVkZ2fzwAMPAIYynTZtGi1btgSM9IO1axsD9MeMGcNnnxnhsyZOnMiSJUvKmII8yWXrKbm5uYwfP56XXnrJGsNq/fr1tG7dmj///JOhQ4fSqVMnBg8e7LKde++9l/Xr11O7dm2mTJlCenq6NclKfn4+zZo1s9a1mHLGjBnDkiVLrGk6LechJyfHbUh0T1wvXdnsawPh5jq2ezoLTHC55wpiSTiu0WgqTkhICAkJCSQkJNCtWzcWLVrkcWKb0jlg8/PzMZlM1vy1YWFhJCQkuMxfW69ePbv1ll6wJbZTRkaGNSfwoUOHaNWqFUopDh48aH0AgJFv2ZHN/6qrriInJ6dM+bx587xKeH/hwgXGjx9PUlISN9xwg7W8devWADRr1oxx48axcePGMso+JibG+tYD8Nprr3HixAl69+6NiDBlyhSeeeYZh/sdM2YMDz74IJs3byYvL4+4uDjrusLCQurWretWdkvHeOHChd4nHBeRtWKkJuwvIk/YTC+KiP5Aq9FUAXbt2sXu3Rdv1y1bthAVFWVXZ+DAgaxcuZKCggJyc3Ot+RGccebMGRo1akRYWBi//vorGzZscFm/UaNGFBcXWxX+tm3biI29GDZr69at1uCEGRkZ1vn09HR27NjBzJkzmTJlivVbQ2nWrVvHli1bykzeKHoRYdq0aXTu3JmHHrpoyDh37pz1QXLu3Dm+/vrrMp48AFdffTUFBQV2eTPy8vIAGDJkCB9//DF//mlEcz916hT79++31gsPDycxMZE77riDSZMmWctPnjxJZGSkz2JdeeKNU0cplQJE29YXkat9IkEASEhIALALJKbRVEdyc3O57777yM7OJjQ0lPbt25d5a+7Tp49HeWQtDB8+3GH+WlcMGzaM9evXc80117Bt2zZrLtqCggLy8/Np1KgRYK/409PTeeGFF+jTpw+rV6/m+PHj5T0NViZNmoTJZOLEiRO0adOGJ554gmnTpvH999/z/vvvW10sAZ5++mk6derEuHHjAGOk+i233MLw4cPLtKuU4tNPP+XBBx/k+eefp2nTptSvX5/nnnuOLl268NRTTzFs2DBKSkq45JJLeO211+weupMmTWLcuHEsWbLEWpaWlsbIkSMrfMxWnAXNsUxABnA3RuKROMvkbruKTL4KhOaM+Ph4iY+P9+s+NJqqRE5OjoiInDt3TuLi4iQ9Pd2n7aenp8vkyZO92mbs2LFSWFgoIiL/+te/ZOvWrT6VKdgZN26c7Nq1y6ttcBEIzZOefZGI/K/7ahqNpqoyY8YMuzyyvk4Y3qtXLxITEykuLvbY1952QNITTzzhU3mCnfPnzzN27NhyZdRzhiexceZixMNZAVj9oUTklM+kKIW/Y+NoM45Go6mOlCvhuA2WFIJ/sykToF1FBXOGdr3UaDQaz/FJ1MtAoHv2Go1G4z0VinqplApTSj1m9shBKXWFUmqUu+00Go1GEzx4Ei7hHeA8xmhagMPAU36TSKPRaDQ+xxNl/xcReR64ACAieYDyq1QajUaj8SmeKPvzSql6GB9lUUr9BRuvHI1G4zssYQ00Gl/jibKfC3wJXKaUSgXWAH/3p1AajcY3hISE0KNHD2JiYoiNjeWFF16gpKTE5TZZWVnWmOya6oNbZS8iX2OkJJwKfAj0FhGTP4WyuF6aXYk0mhpBamoqGzZsYO3atURHR/sklnm9evXYsmULmZmZfPPNN3zxxRduByhpZV/1WLlyJTNmzAAXrpeehEtYCUwC6rur66tJh0vQ1DQWL14sYWFhgmEuFUDCwsJk8eLFFWq3fv36dst79uyRxo0bS0lJiezbt08GDRokPXv2lJ49e8r3338vIiL9+vWTSy+9VGJjY+XFF190Wk8TfOAiXIInyj4eeB3YD3yMEd64rrvtKjJpZa+paURFRdkpessUFRVVoXZLK3sRkYYNG8qxY8fk3Llzkp+fLyIiv/32m1juu7S0NBk5cqS1vrN6muDDlbJ3O4JWRNYCa5VSIcDVwJ3A2xhJTNxi3u5n4LCIjFJKXQ4sAZoA6cCtInLek7Y0murKgQMHvCr3BRcuXGDWrFnW7Em//fZbheppghtPPtBi9sYZD8wE+gCLvNjHbGCnzfJzwHwRaQ+cBvyXcFKjqSK0bdvWq/LysnfvXkJCQmjWrBnz58+nefPmZGRk8PPPP1vzoJbG03qa4MaTEbQfYSjrq4FXMfzu7/OkcaVUG2Ak8KZ5WZnb+dhcZREw1mupNZpqRnJyMmFhYXZlYWFhJCcn+2wfx48fZ+bMmcyaNQulFGfOnKFly5bUqlWL999/35oEu0GDBnaZn5zV01QtPOnZv4Wh4GeKSJqIuPbbsuclDDdNyzZNgGwRKTIvHwJae9GeRlMtSUpKIiUlxZoGMCoqipSUFJKSkirUbn5+vtX18pprrmHYsGE8/vjjANxzzz0sWrSI2NhYfv31V+rXrw9A9+7dCQkJITY2lvnz5zutp6liODPmA3+3mZ9Yat3TzrazqTMKeN08nwCsAiKB323qXAZsL71t27ZtJS4uzjotWLDAZx8wFi9eLHXq1LF+/Kqot4NG40u084DGWxYsWGDVlUCWONHJTqNeKqU2i0iv0vOOlp1s/wxwK1AE1MX4oLsCuBZoISJFSqkBwFwRudZ2W39FvUxNTWXGjBnW3JBgvCr7ogel0Wg0gaa8US+Vk3lHy2UQkUdFpI2IRAM3A9+JSBKQhuG+CUas/M/cteUr5syZY6fowUgKPGfOnMoSQaPRaAKCK2UvTuYdLXvDP4CHlFK/Y9jw36pAW14RCPc2jUajCQZc+dnHKqXOYvTi65nnMS/X9WYnYoRXMJnn92IkL6902rZty/79+x2WazQaTXXGac9eREJE5FIRaSAioeZ5y/IllSmkr6gM9zaNRqMJRjwaVFVd8Jd7m0aj0QQ7Qans/Rn1Mikpif79+xMfH09WVpZW9JpqTXJyMjExMXTv3p0ePXrw008/AfDSSy+VcVZwhKf1pk6dyscff2xXFh4eXj6hNV7jSdTLoFT2DRs2JCUlhdGjRwdaFI3GZ5ROcG8ymWjatGmFy53x448/smrVKjZv3szWrVv59ttvueyyywDfK3tNYBk9ejQpKSkAZ5zVCUplr9FUR2wzUJlMJiZOnMiyZcsqVO6Ko0ePEhkZaTVbRkZG0qpVK1555RWOHDlCYmIiiYmJANx999307t2bmJgY6whbR/XKg8lkIiEhgQkTJtCpUyeSkpIsgyrZtGkTV155JbGxsfTt29cuTIPGxzgbbRXISYc41lRn0tLSJDIyUtLS0ipU7o6cnByJjY2VK664Qu6++24xmUzWdVFRUXL8+HHr8smTJ0VEpKioSOLj4yUjI8NhPWdMmTJFli1bZldmCa+clpYml156qRw8eFCKi4ulf//+sm7dOiksLJTLL79cNm7cKCIiZ86ckQsXLnh1jBp7cBHiWPfsNZpKxFc9ek9MOeHh4aSnp5OSkkLTpk256aabePfddx3W/eijj+jVqxc9e/YkMzOTHTt2eHVcRoxD52V9+/alTZs21KpVix49epCVlcWuXbto2bIlffr0AeDSSy8lNNRt1HVNOdFnVqOpJHyp6CdOnMjx48fd7jMkJMSaxLxbt24sWrSIqVOn2tXZt28f8+bNY9OmTTRq1IipU6dSUFDg1bE1adKE06dPW5dPnTpFZGSkddliSrLIVFRUhKZy0T17jaaS8KWiX7Zsmdv97dq1i927d1uXt2zZQlRUFGAfxvjs2bPUr1+fhg0b8scff/DFF19Ytykd7tgZCQkJLF261Brr/t1333Vr5+/YsSNHjx5l06ZNAOTk5OiHgB8Jyp69xfVy9OjR2iNHU23wpaL35CNtbm4u9913H9nZ2YSGhtK+fXuLxwYzZsxg+PDhtGrVirS0NHr27EmnTp247LLLGDhwoLWN0vWmT5/OzJkz6d3bPtbWqFGjSE9PJy4ujpCQEP7yl7/wxhtvuJSvdu3aLF26lPvuu4/8/Hzq1avHt99+y9mzZ5k+fTr//e9/3R6jxmDlypUWV3WnrpdOo14GEn9FvbRguVE8dWHTaDSaqkB5o15qNBqNppqglb1Go9HUALSy12g0mhqAVvYajUZTA9DKXqPRaGoAQans/Rn1UqPRaKobnkS91K6XGo1GU02oMa6XlmHhGo1Go7GnWil7jUaj0TgmKMMl+BttvtFoNDUN3bPXaDSaGoBW9hqNRlMDCEplr10vNRqNxnNqnOuldqnUaDQ1mRrjeqnRaDQax/hN2Sul6iqlNiqlMpRSmUqpJ8zllyulflJK/a6UWqqUqu0vGTQajUZj4M+efSFwtYjEAj2A4Uqp/sBzwHwRaQ+cBqb5UQaNRqPR4EdlLwa55sVLzJMAVwMfm8sXAWP9JYNGo9FoDPxqs1dKhSiltgB/At8Ae4BsEbFkFT4EtPanDBqNRqPxs7IXkWIR6QG0AfoCnTzZ7vjx4/Tu3ds6WZIka7xHxwvSaKo3KSkpVl0JRDqrVynhEkQkWymVBgwAIpRSoebefRvgcOn6TZs2xZ9RLzUajaa6MGPGDIuPPUqpE87q+dMbp6lSKsI8Xw8YCuwE0oAJ5mpTgM/8JYNGo9FoDPzZs28JLFJKhWA8VD4SkVVKqR3AEqXUU8AvwFt+lEGj0Wg0+FHZi8hWoKeD8r0Y9nufkpqayoYNGygsLCQ6Oprk5GSSkpJ8vRuNRqOpklSLEbSpqanMmDGDwsJCAPbv38+MGTNITU0NsGQajUYTHFQLZT9nzhzy8vLsyvLy8pgzZ06AJNJoNJrgIiiVvbdRLw8cOOBVuUaj0VQnakzUy+joaPbv31+mPCoqiqysLB9KVvXQkUA1mppDtY96mZycTFhYmF1ZWFgYycnJAZJIo9FogotqoeyTkpJISUmhTp06gNGjT0lJ0d44Go1GY6baJBxPSkpi4cKFgDZZaDQaTWmqRc9eo9FoNK7Ryl6j0WhqAEGp7HXCcd9gGVW8du1aoqOj9SAzjaaaUmNcLy1oN8OLWEYV2w42CwsL0x+uNZpqTLV3vdSURY8q1mg0tmhlX03Ro4o1Go0tWtlXU9q2betVuUajqd5oZV9N0aOKNRqNLVrZV1P0qGKNRmNLUI6gtbhejh49mtGjRwdanCqLHlWs0dQMVq5caXFVd+p6GZTKvmHDhqSkpARaDI1Go6kSWDrGCxcuPOOsjjbjaDQaTQ1AK3uNRqOpAQSlGae8aLu0RqPROEb37DUajaYGUCOUfekev8lkomnTpvpNQKPR1BiCUtn7OuqlJUAaGIp+4sSJLFu2rEy5RqPRVEVqXNRLpygFgAmYCCwDEmxWW8qPB+G5qCg6EqhGU3PQUS9xr+iXVbpEGo1GU3n4zRtHKXUZ8B7QHBAgRUReVko1BpYC0UAWcKOInPaXHBYSgONelGs0Gk11wp89+yLgYRHpAvQH7lVKdQEeAdaIyBXAGvOyRqPRaPyI35S9iBwVkc3m+RxgJ9AaGAMsMldbBIz1lwwajUajMagUm71SKhroCfwENBeRo+ZVxzDMPBqNRqPxI34fQauUCgeWAw+IyFll9owBEBFRSpVxgTl+/Di9e1/8oDxjxgyLW5HGS7QXTs1Ge2NVf1JSUmwDR0Y6q+dXZa+UugRD0aeKyCfm4j+UUi1F5KhSqiXwZ+ntmjZtik9dLzUajaaaYtsZVkqdcFbPb2YcZXTh3wJ2isiLNqs+B6aY56cAn/lLBo1Go9EY+LNnPxC4FdimlNpiLvsn8CzwkVJqGrAfuNGPMmg0Go0GPyp7EVkPKCerh/hrvxpNTSBQtnj9DaDqUuVH0DoLaqaDnWk0Gs1FqrSydxXUzFG5RqPR1FSCMnmJpwnHJyYmGrFuEhOtZSZsYuDYlGs0Gk11pdonHHcV1CyhbHWNRqOplniScDwolb2nJDhY1kHNNBqNpixVWtlryqI2e1ZPevlXDo1GE1xU6Q+0Go3GOampqWzYsIG1a9cSHR1NampqoEWqUiQkJFQrBw+t7DWaakhqaiozZsygsLAQgP379zNjxgyt8GswWtlrNNWQOXPmkJeXZ1eWl5fHnDlzAiSRJtAEpbL3dcJxjaY64Yl55sCBAw63dVbuq/1qAkP1TziunEVjKCdBeC68RX+grd5YzDO2vfawsDBSUlJISkqylkVHR7N///4y20dFRZGVleW3/XpDsIdeCHb5HKETjms01QRPzTPJycmEhYXZlYWFhZGcnOzX/WqCF63sNZoqhKfmmaSkJFJSUqhTpw5g9Ogr0gv3h1lIU7loP3uNpgrRtm1bh+aZtm3blilLSkpi4cKFQMVNEd7sVxOc6J59FcHf0T2D2ac4ULIF4znxtXnGX/sNxnNX09HKvoqgo3tqwPfmmWDfr8Z3BKUZx9OolzWJE8nLSLw0ASzeNj+b4B8T4blS5Zpqjy/NM1Vhv4HA4mZaWFhIdHQ0ycnJQf1g8yTqZVD27C1RL7Wit6F3wsV5G0VfpnxIU+PXFmflGo2mDFVx9PHo0aMtkYKdRr0MSmWvcYErRe9NeTnRtlhNZVPZ11wg3Uz9eaxBacbRlEXilH28/rscJ2xJ9PQB0MtmWaOpZILZTFJd3Ux1z76KYMJxYhZn5W57+hpNgAh2M4kzd9Kq7maqlX0VwaeK/rll/hRVo3FJsI/GDZR7q78JSjOO9sYpi08Vfe+EsjF0co2f0uX+jKFTFWOPVDWC8dz6y0ziq+vJYk6aNm0ahYWFREVFBZWZCcoeq/bGqUYk2MybcN7Tj7wrkbTsE8hdiUicQuIUaXGqTLkdX6TCtg2QvhZGRhvLQUKgIi3qCI/lx925qwpmkqSkJPr37098fDxZWVlBpegdob1xqiEmvDPpOCu38kUqPDUDzhv2U47tN5aDQOEHyrYb7DblYOaxxx5j8uTJZc7dY489RtOmTYHqayYJdrSyr2IkYCRVT/Cg3ITzB4CVV+dAgb39lII8ozzABMq2G+w25WDFZDLxzDPPlCnPy8vjmWeeYdky41uRHo0bIETELxPwNvAnsN2mrDHwDbDb/NvI0bZxcXHiEUYEet9NwYyXx5IGEmn+dVROuhiTUgKUnZQS0u1FWLx4sdSpU0cAiYqKksWLF1fokOLj4yU+Pt7peuVENqVUhfbrjkDttypjua6Uo2vJXF76HnP3/4t4d8150p43+Lo9d/jiWIGfxYlO9mfP/l1geKmyR4A1InIFsMa8HLQE0wCiWZO22U2tm/WmdbPeZcrBM5OOleZO7KTN29qNuPXWtOGLcxco225VsCk7o7zB8ir6f1muK2dnqDxnrqqY0zw9d67qVcax+k3Zi8j/AadKFY8BFpnnFwFj/bX/mooJL233100u20jdMKPcxh8/EKaNqhLhMZgIVLA8y3WVDISVWhdmLvcWV9dcMHoZVYTKuL8q22bfXESOmuePAc0dVTp+/Di9e/e2TuavzBoP8Nof/5MFcMccqG3YT2kRBbc8aJTb+OMHYlShjvDoPZ4oen8oSkvrSUAKYL6aiDIvW86cN28erq65iROr18DAitxfKSkpVl0JRDqt6My+44sJiMbeZp9dav1pR9sFi82+sm12rrh30ja7qVWz3tKqWe8y5Z4eKwvShIhI4zddhLh4YypVbiEqKsqhLTYqKkpERNLS0uzkjY2NldDQ0DLlFgJpi/WUYPr/PcbDbzae/l9paWkSGRnpvrzUfuLNU5lvSTbtxMfHS2xsrNP2a9Wq5fCaq1WrVpn6vv6eZJHPk//fF/Xc3V+2uDpWAmSzd8QfSqmWAObfPyt5/0FPpX0ncDTwKifbadA0d6aN0j3HzMxMYmJiHB6Lv+yTwfSNpTQVlc2bHrEJD97wPPi/XL0ZlNc0lJCYCEqBUmSvXUtmRgbLTpywKzcpxcTERB4tKSljEgJ4tKTEqG/G3fXkr4Q/vsRT02FF7p3KHkH7OTAFeNb8+1kl77/ScXfDHD9+3Kf7s3ykdUfah93sgqn1ADIxXKRsyzHewDweVWg5rpiYGCIiIuzKLcfvyj5ZFUwlgcBThWvCQ1OeUnblMUBERkaZ8mVgp1idlXuLCeN6i3EjZ2dgGlCIYXN+FHjKtr7JxG233UZJSYld+3l5edx22220bt3a4wdY6QdXZecKAPf3V4XuHWdd/opOwIfAUeACcAjjP2uC4YWzG/gWaOxo22Aw4/jitdDjV2Cb8tDQUImNjS1TbmuqGTrgWQmpVVsACQ9rKUMHPFvGnONuKv1qHwoS6+CVvzSuXkUtxzVnzhy7czdnzhyJjIy01vPWtdHXr9Oe4sv2fNmWq+vHm+vdYtKZA1LHYjYwL7s0ATm6Tjw043hbx51JKDIy0uPryd39WB481ROe1quoOzIuzDh+69mLyCQnq4b4a5++wtmrEuBVz9ObV2BXPeKJEydy09A0AH7LWk3axrkUl5wHIDfvKGkb5wLQIXqk18dqwqZn56Dc0/cOk1JMBO4C5icnYx6Py/79+0lOTsbWp0Anry4HHva4vXlPTABeAmbAxf8LmI/xUTXBpq7dfh2U+/b91DMmJiayDJiKIXdp2op49abiLZ7qCV/pE6jYvVMlRtBW1eQFtgo9ISGBHj16uH0AOFL0lpGHAD9mvExRcYHdfoqKC/gx42WvZAP7Cz3CSbmnWOovBvIcrF9sMx9o10ZfX0+VdX2a8GL8hIfMoez/lWcutyUB1yO3y4OJUqO5vcRTd08T9ufNMpUut8rl7FuIudyCp3rCl26VFbl3dNRLB/jKzdBWAWRnZ5OZmck333xjV65STPa5ZHMdl9/LdgBy84453Jezcpfy4fhGdVbuCssN4+wM2ZZXhaiCwYaJcsY+coMn/1ewkmD+tVw1Ftt+FIaiT7Kp5+l17slHaQue6glfui07u3cuvfRSy9tCQ2fbBmXPPtBRL302gtLGuyAzI4OYoqIyXge20SgXxynqpK+F9LXUuiuROQ6iVIaHtXC4K2flvsBdT8dkMllvPE9HUHoSVdDSfnZ2tsNyW1xFWvREfn/iiwiaCXgeE8kbfDniNZAkAf2BeCCLi4reGyxeQM68g6zlZtzpCct15esR2Y7uHR31spz40sxgwrVN3NIjS8XedlqCYTstrRYGxM4mNKSuXVloSF0GxM72WrbScpoclXvpfudgPC5h5nJHijU7O9upIra0v2XLFut6Rz2s1NRUa08H7N3RnMmZnZ3t1n3QZDK5fRCYTCa+//57hw8ki2zBPOTflyNebTFRMRNNIPAqmqzJxMmTJ60D7yzUqVOHkydP2nn3JCcnl6kXkBHZzr7cBnIq7Y3j9At1sHvj2Hgv2HoUOPJqiHLwhR1zeRr41RvH1VTae8HpQBibbRaX8u5YbOM9YUuFB/KYad68ucNz17x5c4f13Q348hSLPLGxsXbXp613hzeDZUQq8Vr34P8Klsn23vFnvTQHy+68kmzPXXOQcCf1w0EuKXWO08Dz/9/DerjwxnFYGOjJVtm7VLq+vrDKeeKd4c0F4SpaYCSUUdbORtD6evL0BvDohrLBlZupt+6qjs6bZSpdf86cOWKrcC3Xk7eK352bqaU9d65ypY/Lcs25G6Hq68lTRRnMkyfH4IsHW+nrPx7Dbdndg8FWPku5t9ebu+vElbIPyg+0FnzpshQIjlPWPPMHxmvy4VJ12+LEfQzDBjnLPyK6JcFm3oTzV13bZaf4eCCPpbw5xnktTfNS9R/D3jyxf/9+ZkyezM7Jk1kAHBfx5CispqG77rqL+fPn212fycnJzJkzx/oK78pVriImJo33lL4X95uXwTsbv6OPvRHAFgflCU7qOyt3hi+uk6C22Vd2pEVvbMqeMGvSNmaGtXTo2jYzrKXdaFd/2U59hQnXNs3SoZZ3N+vN7lIhmEu3E+FF+67KX6DsuatjLretXzathvFfPIOXbqbmj3WLk5PLXJ8Ai5OTrQ+wyZPLfsUICwtj8uTJTm/gzMxMrej9gKdupoHAnRNB6Y/GzkJNuCIoe/YW18vKjrTo7MYr7S7pDZ66SrpzHws0CbjuoXjy5mHCXnGbnJQ7q++u/HaM4drNMRR9Uqn6zvrtgnfeLJ66mZpMJhYsWMCcOXOYN2+e1VVu8uTJLFiwwLjeSr2pWEMIlCr3Rr6aisnN+qB0M1XKszdam01MGJaBb2zKVwILjNmGznYVlMre4nr59ddfV+pIy8RLE2CzeeFnE2zPhL/ElCmXGQketxke1oLcvKMOy6FsLJsma24HYNSQd/gR+NGbAwhyvPUTT8DzV+DWGEo7FvvXadv2p+LcVOYNCTbbOWvPhP2r9/r16wGYO3eu05g2owEFZGCEi002H1egRqhWN1z9X4HChG/ui9FYrx+nrpdBqewtJCcnM2PGDLtXZX+6LFn82U0YJ7o1EPFbBqZS5czwzLYLhqtk2sa5dqNefeEqWRXx1QAuRyQAA920n4xho7V9la+IqWyyg20tbqYToYxfdjYXh/iX7sGNxngruWAu24/xlncJRq9NU3F8/f/7ggT8e1/YEtQ2+0AkkTBx8Ym6BcemBm/oED2SxL5zCalVG4DwsJYk9p1brjg2towb8g7jhrxToTb8STDKl4TvEmuAEX1xsYP2nqLsQKdsDBONs55aEy5+OLRQaC4vXd8fmKh6fvHe4u7/Ly8mPDt3ntbzF0HdswdD4S9cuBConJCjCfj+SdsheiSZez4GCDoFWNNIAhaa502l1pUnbrur9rApdxbO17IclPbkaogn/1d1JeiVfbBiqwjmzp1rLbMf4ek8Q1hVwdP4+NUBT909E7xtF+cmJgvBaE/WVC+CUtnn5OTYLVu8YhwlGAgEnoYQ+Dhlu9dt655/+R8wXS3bu6m32/wRfFapc/3qh90A/4TzNblZH4z2ZE31IiiVfclvvzFDKUYDDQg+dzQ9EKZ8VIW3BBPuvSNKH4e7h4cnBLvrrSa4WYn1Q34Vc73E+HBiwd0rcGVj62WRCvzARf/uw7ap2qqActPYk4D7bzYf+2nfNdmerKkYFtfLhVXV9dKCKdAClCLB/JuK0ROzuMv9wcXh161LbaPNM8GDr/+LyvhvPX0r8uZtoqZiCrQAAaJKKPtgZQ5l3eUsw6+z8F8PUFM18ERBa+WsqSyC2s8+2NHuchqNpqqge/YVQLvLaSpK6d6/px5FnrbnDF+/UWgzU/CjlX05sFzYHbJWc9hBKIQOfecyq4IjZDWVT1XwFqoJ6AeHfwhKZX8G40On5QtzsGIJefDdT/+iuOQ84WEtGRA7u8KhEDQaje+oCQ+PauN6WVmUp2enQyFoNBpHVOYH+mrjeumMmvDE1mj8QXW4d3xtdgvmc+LxsbqQrUore0/xty1W9+g1Gs/Q30UCR0CUvVJqOPAyEAK8KSLP2q4P9kQNmb8vI6b9xECL4RAtW/kJZvkCJZunyrkmnTtfPrBSuDgQ099Uup+9UioEeA0YAXQBJimlutjWOeFhW/sOmwJSz2Kj90V7wSybr+sFs2wQGPmCWTZv6ul7onz1HmvcpUz+ZkeTp/vExQfaQAyq6gv8LiJ7ReQ8sAQYY1vhfEhdj05AlocnwNf1PMWT9oJZNn/UC8Q+g/ncBbNs3tTzFH1PlA8v2opwtkKJeJ5izxcopSYAw0Vkunn5VqCfiMyyqVMAFNtsdhzHHf6GuPj67Md6kU7kKU97wSybr+sFs2wQGPmCWTZv6ul7onz1fCFbJNDUPF9LROo5qhSUH2hFpG6gZdBoNJrqRCDMOIeBy2yW25jLNBqNRuMnAqHsNwFXKKUuV0rVBm4GPg+AHBqNRlNjqHRlLyJFGHGevgJ2Ah8Bw5RSmUqp7UqpD5VSdc0Pg5+UUr8rpZaaHwx+Ryn1tlLqT6XUdpuyxkqpb5RSu82/jczlSin1ilnGrUqpXgGS799KqV/NMqxQSkXYrHvULN8updS1lS2bzbqHlVKilIo0LwfFuTOX32c+f5lKqedtygN67pRSPZRSG5RSW5RSPyul+prLK/XcKaUuU0qlKaV2mM/RbHN5wO8LF7IFyz3hUD6b9ZV3X4hIQCeMPB/7gHrm5Y+Aqebfm81lbwB3V5I8g4FewHabsueBR8zzjwDPmeevA74AFNAf+ClA8g0DQs3zz9nI1wXIAOoAlwN7gJDKlM1cfhnGw30/EBlk5y4R+BaoY15uFiznDvgaGGFzvkyBOHdAS6CXeb4B8Jv5/AT8vnAhW7DcEw7lMy9X6n0RLPHsQ4F6SqlQjDzLR4GruZj/YxEwtjIEEZH/A06VKh5jlqG0LGOA98RgAxChlGpZ2fKJyNdivDEBbMD4DmKRb4mIFIrIPuB3DNfXSpPNzHzg74Ct61dQnDvgbuBZESk01/nTRr5AnzsBLjXPNwSO2MhWaedORI6KyGbzfA7GG3lrguC+cCZbEN0Tzs4dVPJ9EXBlLyKHgXkYOT+OYrgXpQPZNn/WIcpm+qtMmovIUfP8MYx0s2DIdNCmXqDlBLgDo2cAQSCfUmoMcFhEMkqtCrhsZjoAV5lNhmuVUn3M5cEg3wPAv5VSBzHukUfN5QGTTSkVDfQEfiLI7otSstkSFPeErXyBuC8CruzNdr4xGK9UrYD6wPCACuUCMd61KndwgocopeYARRjpcQOOUioM+Cfwr0DL4oJQoDHGK/PfgI+UUiqwIlm5G3hQRC4DHgTeCqQwSqlwYDnwgIictV0X6PvCmWzBck/YymeWp9Lvi4Are+AaYJ+IHBeRC8AnwECM1xfLOIBAu2f+YXmVMv9aXvWDxo1UKTUVGAUkmW88CLx8f8F4iGcopbLM+9+slGoRBLJZOAR8Yn5t3giUYAxSCQb5pmDcDwDLuGhuqHTZlFKXYCirVBGxyBQU94UT2YLmnnAgX0Dui2BQ9geA/kqpMHOPagiwA0gDJpjrTAE+C5B8YLiGTnEgy+fAbeYv6P2BMzavtZWGMgLL/R24XkTybFZ9DtyslKqjlLocuALYWFlyicg2EWkmItEiEo2hWHuJyDGC5NwBn2J8pEUp1QGojTGiMaDnzswRIN48fzWw2zxfqefOfF++BewUkRdtVgX8vnAmW7DcE47kC9h94asvvRWZgCeAX4HtwPsYX8rbYfwJv2P0aupUkiwfYnw7uGD+E6YBTYA1GDfbt0Bjc12FEdRtD7AN6B0g+X7HsPNtMU9v2NSfY5ZvF2bPjsqUrdT6LC56HQTLuasNLDZfe5uBq4Pl3AGDML5fZWDYoeMCce7Mcgiw1eYauy4Y7gsXsgXLPeFQvkDcF5UeG0ej0Wg0lU8wmHE0Go1G42e0stdoNJoagFb2Go1GUwPQyl6j0WhqAFrZazQaTQ1AK3uNRqOpAWhlr9FoNDUArew1QYVSyqSUeizQcniLUuoLpdTfa7oMNQ2lVJZSqkAplRmg/X+tlMpXShW5q6uVvaZSUErVUkr9YE7U0Mb9Fj7bb6U8PERkhIg8775m1ZChqj50A8R0EYmxLVBKxSmllisjIU2u+aGwXCl1tScNKqU+U0q952RdmlLqVQARGQaM8KRNrew1lcWDQJ7bWhqvMQfa0gQJSqmhwPcYIQ96YyQt6QZ8AIzzsJkFwARlk2HL3PYVGPGSFngtmD/jQuhJTyICRsz4PUAPjDghbVzUNQGP2SyHYcRy34eR3ONLoL3N+iyMULHrgVzgZ6CPed2rQDFQaF63y1zeBHgPIwb7MYzEG41LtflPjLgvuRhxc650c4xWud1tD9wLbCm1/eVmWaOB2RixonIwAgU+g002JZtjTjO3f7MDGTxpw6GMzs6bk+POAh6zkWUb0B2YhBGf5gzwJhezRrmT637zf52DEe3xaU/WVVSuClzbWcDkUmW/A296sK3TaxujI74fuK/UNv8GfixVlgAUud1foBWBnqr3ZL5o12P0aKLxXtmnAqswEmPU5mLQvEvM67MwokPGmdc/AhwHLnXUnrnsS2Al0Mg8rQZW26zPMt+wMUAIRkah3W6O07ofd9ub91kA9LApewJYY54fj6H8FUayiz+Au0rJd9C8TnExpaetDJ604UrGMufNyXFnYQRC6wxcghFUbg+QgpGboi1G6OMkd3JhdArygBjzcgTQ3906X8hVges7Cxtlb5ZTgCEebOvu2v4fYKtN/dpmmaeWaicBrez1FOgJw3zzsXk+Gi+UPUZceQHa2qyvhdErG2RezgKetFmvMHqMt5Ruz7zcytzmFTZlHc1lLW3a/JvN+hjz+oYeyu12e2Ap8LKNzFnOFA9G7+8jm+Us4F+uZPCwDacyumqrVLul27nO3E5Tm7KPgPnu5MKIdJsP3AiEl6rndJ0v5cJ46zuNWYFjXLPfutmPrbIfaN5PJ5uy64FsjOu2wItruxVGFNR+5uWbzLLVKyVDAh4oe22z1/gNpVR74GFglpP1SeaPV7lKqVwHVS43/25VSmUrpbIxXncvwT7BQ5ZlRoyr/wAXc46WxrLdPpuyPaXWgRFu2MI5828DJ206wt327wC3mO3tV2P0VD8BUEpNUkptUkqdVEqdwTD7NC3VfparnXvYRkWP0VE7eUCxiBwvVdbAnVwishdIAu4Ejiil1iulhrlb5wu5bHgM4020vJww/1qvPxH5XEQigJEY4dvBg2tbRI5g9PxnmOvOABaLSH55BNPKXuNPBmHcyNuVUicw4sWDcYHfIyKpIhJumRxsv9/8e4WIRNhMYSLyoU29aMuMOVlEW4yY8GBknrLlYOltMHqNtusqg28wbOKjgakYSbDzlVKXYZgcnsJ402iIEd+8dKrE0sdlxYs2XOG0/fLiiVwi8omIDMXo+X4EfKaM9JYu1/lIvvYYPfv0CjTzG7AXuNlNPU+v7RTgJqVUT4wkO95/mDWjlb3Gn3yEkYKth3m6zlw+DOMDqUtE5E8MD4bXlVKtAZRSEUqpccrI6WnhDqVUL3Mv+W8YH75Wm9cdA9rbtHkE+Bp4wdxWI+AF4AupxExZIlKMcQ7uB24A3javCse4L48DF8zZim71snlftGF33nyES7mUUh2VUsPNCvwChklDgBJX63wo3//DsJuXG/Ob5b3ArUqp55RSl5mzToUB/WzqeXptf4XxtrAc48Ps9vLKppW9xm+ISJ6IHLJMGAoE4JiIODLbOOJOjIxCJqVUDoZXxUTsk1unAK9g2DNvAkaKyBnzuvlAb/OrsmXgy2QMj45dGB/EsoHbynOMFeQdDDe6fWLkv0VEdgKPY6T4y8b44PyhswYc4Ys2cHzeKoQHctXG8DI6al5/PzBeRArcrKswSqkrgZMissdtZTeIyJcYb7UdMN5mc4FMDHu+rZ+922tbREqAhRhmn5SKyKUzVWmqNMpI2PyYiCwOtCyaqotS6n4MT6F8jDeac8BMjIfLmyJyjZPtdgEtgSwR6V5J4tru/wuMh0gtJ6ZQK6GVI5JGo9EELyLyCsbbIUqpucDvIvKjUirazXYd/S+dy/17NHoWtLLXaDQaO0Rkrs18FuCwV1/V0GYcjUajqQHoD7QajUZTA9DKXqPRaGoAWtlrNBpNDUAre41Go6kBaGWv0Wg0NQCt7DUajaYGoJW9RqPR1AC0stdoNJoawP8HuhpjqYRiNfQAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "plot_data(data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## What can you do to explore this analysis?\n", "\n", "* Increase the fraction of data used in '[Lumi, fraction, file path](#fraction)'\n", "* Check how many events are being thrown away by each cut in '[Applying a cut](#applying_cut)'\n", "* Add more cuts from the [Higgs discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X#se0040) in '[Changing a cut](#changing_cut)' and '[Applying a cut](#applying_cut)'\n", "* Add a plot to show the ratio between Data and MC other than Higgs\n", "* Add a plot to show the invariant mass distribution of the sub-leading lepton pair, like [Figure 1 of the Higgs discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X#fg0010)\n", "* Get the estimated numbers of events, like [Table 3 of the Higgs discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X#tl0030)\n", "* Add a plot of m12 against m34, like [Figure 3 of the Higgs discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X#fg0030)\n", "* Your idea!" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Back to contents](#contents)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.9.0" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 2 }