{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "-NFVjVZ0k1Bk" }, "source": [ "## How to rediscover the Higgs boson yourself!" ] }, { "cell_type": "markdown", "metadata": { "id": "WYFVOpVmk1Bm" }, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": { "id": "UraN3pgPk1Bm" }, "source": [ "This notebook uses ATLAS Open Data http://opendata.atlas.cern 2025 beta release 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" ] }, { "cell_type": "markdown", "metadata": { "id": "8lkh5zZrk1Bm" }, "source": [ "### What is the Higgs boson?\n", "The Higgs boson is a fundamental particle predicted by the Standard Model.\n", "It is a manifestation of the Higgs field,\n", " which gives mass to the fundamental particles.\n", "However,\n", " it is incredibly hard to produce.\n", "At the LHC,\n", " a Higgs particle is produced about once every 10 billion collisions!\n", "This tiny fraction makes it very difficult to detect.\n", "Nevertheless,\n", " after years of data collection,\n", " the Higgs boson was finally discovered in 2012 by CMS and ATLAS experiments at CERN.\n", "In this tutorial,\n", " we shall be following their example.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "e_Pa25dak1Bm" }, "source": [ "### Detecting the Higgs\n", "This analysis loosely follows the paper on the [discovery of the Higgs boson by ATLAS](https://www.sciencedirect.com/science/article/pii/S037026931200857X) (mostly Section 4 and 4.1).\n", "\n", "The Higgs boson can be produced in many different ways.\n", "In particle physics,\n", " we describe these production modes using Feynman diagrams.\n", "These diagrams allow us to visualise particle processes while also acting as powerful tools for calculations.\n", "See [here](https://cds.cern.ch/record/2759490/files/Feynman%20Diagrams%20-%20ATLAS%20Cheat%20Sheet.pdf) for more information on Feynman diagrams.\n", "\n", "There are four main production modes of the Higgs boson, and their respective Feynman diagrams:\n", "1. Gluon-gluon fusion (top left)\n", "2. Vector boson fusion (top right)\n", "3. Vector boson bremsstrahlung (bottom left)\n", "4. Top-antitop fusion (bottom right)\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "\n", "The Higgs has a very short lifetime, on the order of $10^{-22} \\,\\text{s}$. It decays extremely quickly after production, so there is no hope of directly detecting the particle. Nevertheless, we can use the Standard Model to predict its decay products: photons, Z bosons, quarks, etc., all with different probabilities. These **decay channels** can be used to identify the Higgs boson. In this notebook, we'll be looking at one particular decay channel: $$H \\rightarrow ZZ^* \\rightarrow \\ell\\ell\\ell\\ell$$\n", "\n", "
\n" ] }, { "cell_type": "markdown", "metadata": { "id": "hR8omi4uk1Bm" }, "source": [ "We refer to this as our desired **signal**.\n", "Ideally,\n", " we would search for collisions which yield four leptons as products and this would tell us that a Higgs boson is present.\n", "Unfortunately,\n", " in addition to our signal,\n", " there are many other **background** processes that lead to four reconstructed leptons in the final state.\n", "The main background is $ZZ^* \\to \\ell\\ell\\ell\\ell$,\n", " where decay products have the same properties as those in the Higgs decay.\n", "This is known as an irreducible background.\n", "
\n", "\n", "We can get around this by accounting for the total invariant mass of the lepton products.\n", "We know through conservation of energy and momentum that the invariant mass of the products must be equal to the Higgs mass, while other background processes will have different invariant masses.\n", "Our last step would be to plot the invariant mass of each event and spot the peak in mass around $125\\, \\text{GeV}$, which corresponds to the mass of the Higgs boson.\n", "\n", "We also have background contributions from $Z+$ jets and top-anti top processes, where additional charged leptons can arise either from semi-leptonic decays of heavy flavour or light flavour jets misidentified as leptons. We can also have contributions from processes with more than four leptons, but where one (or more) of the leptons are not reconstructed, and thus do not appear in our final selection. These backgrounds are difficult to remove completely.\n", "\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
\n", "\n", "For such processes,\n", " we will attempt to distinguish them from the Higgs decay using the properties of the leptons.\n", "Because the Higgs is a neutral particle with zero lepton number,\n", " the lepton products from its decay must sum to zero charge and zero lepton numbers.\n", "Thus,\n", " we can cut away all data with products that do not have these properties.\n", "These cuts increase the ratio of our signal to the reducible background.\n", "\n", "Note: $Z^*$/$W^*$ refer to a $Z$/$W$ boson that is off its mass shell.\n", "This means that its mass is not fixed to the $91/80 \\, \\text{GeV}$ of a typical $Z$/$W$ boson.\n", "\n", "By the end of this notebook you will be able to:\n", "1. Learn to process large data sets using cuts\n", "2. Understand some general principles of a particle physics analysis\n", "3. Discover the Higgs boson!\n", "\n", "See [here](https://cds.cern.ch/record/2800577/files/Signal%20and%20Background%20Physics%20Cheat%20Sheet.pdf) for more information on signals and backgrounds!" ] }, { "cell_type": "markdown", "metadata": { "id": "iTd_GrqTk1Bn" }, "source": [ "### Running a Python notebook\n", "A Python notebook consists of cell blocks,\n", " each containing lines of Python code.\n", "Each cell can be run independently of each other,\n", " yielding respective outputs below the cells.\n", "Conventionally,\n", " cells are run in order from top to bottom.\n", "\n", "\n", "- To run the whole notebook, in the top menu click Cell $\\to$ Run All.\n", "\n", "- To propagate a change you've made to a piece of code, click Cell $\\to$ Run All Below.\n", "\n", "- You can also run a single code cell, by clicking Cell $\\to$ Run Cells, or using the keyboard shortcut Shift+Enter.\n", "\n", "For more information,\n", " refer to [here](https://www.codecademy.com/article/how-to-use-jupyter-notebooks)." ] }, { "cell_type": "markdown", "metadata": { "id": "GaPCqRhpk1Bn" }, "source": [ "## ATLAS Open Data Initialisation" ] }, { "cell_type": "markdown", "metadata": { "id": "iPHL9OIbk1Bn" }, "source": [ "### First time package installation on your computer (not needed on mybinder)\n", "This first cell installs the required python packages.\n", "It only needs to be run the first time you open this notebook on your computer.\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 this is opened on mybinder, you don't need to run this cell." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "R29WD4HZk1Bn", "scrolled": true }, "outputs": [], "source": [ "import sys\n", "import os.path\n", "!pip install atlasopenmagic\n", "from atlasopenmagic import install_from_environment\n", "install_from_environment()" ] }, { "cell_type": "markdown", "metadata": { "id": "jdPF_-bek1Bn" }, "source": [ "We're going to import a number of packages to help us:\n", "* `numpy`: provides numerical calculations such as histogramming\n", "* `matplotlib`: common tool for making plots, figures, images, visualisations\n", "* `uproot`: processes `.root` files typically used in particle physics into data formats used in python\n", "* `awkward`: introduces `awkward` arrays, a format that generalizes `numpy` to nested data with possibly variable length lists\n", "* `vector`: to allow vectorized 4-momentum calculations" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "id": "RA7rxOESk1Bo" }, "outputs": [], "source": [ "import numpy as np # for numerical calculations such as histogramming\n", "import matplotlib.pyplot as plt # for plotting\n", "import matplotlib_inline # to edit the inline plot format\n", "#matplotlib_inline.backend_inline.set_matplotlib_formats('pdf', 'svg') # to make plots in pdf (vector) format\n", "from matplotlib.ticker import AutoMinorLocator # for minor ticks\n", "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 # for printing time stamps\n", "import requests # for file gathering, if needed" ] }, { "cell_type": "markdown", "metadata": { "id": "MkLLHOz_k1Bo" }, "source": [ "Unit definitions, as stored in the data files" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "id": "XRmyeETUk1Bo" }, "outputs": [], "source": [ "MeV = 0.001\n", "GeV = 1.0" ] }, { "cell_type": "markdown", "metadata": { "id": "4fsnt-oSk1Bo" }, "source": [ "We will use the [atlasopenmagic](https://opendata.atlas.cern/docs/data/atlasopenmagic) to access the open data directly from the ATLAS OpenData Portal so no need to download any samples. First we need to install the package" ] }, { "cell_type": "markdown", "metadata": { "id": "vdvTTNVtk1Bo" }, "source": [ "Import the module and load the release." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "3ecf42YPk1Bo", "outputId": "07282414-6718-4ce5-e754-f26dc5f78563" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Available releases:\n", "========================================\n", "2016e-8tev 2016 Open Data for education release of 8 TeV proton-proton collisions (https://opendata.cern.ch/record/3860).\n", "2020e-13tev 2020 Open Data for education release of 13 TeV proton-proton collisions (https://cern.ch/2r7xt).\n", "2024r-pp 2024 Open Data for research release for proton-proton collisions (https://opendata.cern.record/80020).\n", "2024r-hi 2024 Open Data for research release for heavy-ion collisions (https://opendata.cern.ch/record/80035).\n", "2025e-13tev-beta 2025 Open Data for education and outreach beta release for 13 TeV proton-proton collisions(https://opendata.cern.ch/record/93910).\n", "2025r-evgen 2025 Open Data for research release for event generation (https://opendata.cern.ch/record/160000).\n", "Fetching and caching all metadata for release: 2025e-13tev-beta...\n", "Successfully cached 374 datasets.\n", "Active release: 2025e-13tev-beta. (Datasets path: REMOTE)\n" ] } ], "source": [ "import atlasopenmagic as atom\n", "atom.available_releases()\n", "atom.set_release('2025e-13tev-beta')" ] }, { "cell_type": "markdown", "metadata": { "id": "VM4fn_Qbk1Bo" }, "source": [ "## Example 1: Reading data" ] }, { "cell_type": "markdown", "metadata": { "id": "Uq7c_dzwk1Bo" }, "source": [ "We would like to read some of the data from the open dataset." ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "id": "NAngOYOOk1Bp" }, "outputs": [], "source": [ "lumi = 36.6 # fb-1 # data size of the full release\n", "fraction = 1.0 # reduce this is if you want the code to run quicker" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "id": "pQovxMkJk1Bp" }, "outputs": [], "source": [ "# Select the skim to use for the analysis\n", "skim = \"exactly4lep\"" ] }, { "cell_type": "markdown", "metadata": { "id": "kjbH2F2jk1Bp" }, "source": [ "For convenient naming and identification purposes,\n", " we define a dictionary which stores all the important names of the samples we want to pull from the database." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "pGK-65Z2k1Bp" }, "outputs": [], "source": [ "defs = {\n", " r'Data':{'dids':['data']},\n", " r'Background $Z,t\\bar{t},t\\bar{t}+V,VVV$':{'dids': [410470,410155,410218,\n", " 410219,412043,364243,\n", " 364242,364246,364248,\n", " 700320,700321,700322,\n", " 700323,700324,700325], 'color': \"#6b59d3\" }, # purple\n", " r'Background $ZZ^{*}$': {'dids': [700600],'color': \"#ff0000\" },# red\n", " r'Signal ($m_H$ = 125 GeV)': {'dids': [345060, 346228, 346310, 346311, 346312,\n", " 346340, 346341, 346342],'color': \"#00cdff\" },# light blue\n", "}\n", "\n", "samples = atom.build_dataset(defs, skim=skim, protocol='https', cache=True)" ] }, { "cell_type": "markdown", "metadata": { "id": "HRqMu0-Jk1Bp" }, "source": [ "The key named `data` refers to the event information collected from real experiments,\n", " while the `Background` and `Signal` keys refer to Monte-Carlo (MC) simulations of the ATLAS experiments.\n", "Both real data and MC data will then be analysed and compared together to discover the Higgs!\n", "\n", "Let's try accessing `data15_periodD` in the CERN database URL as an example." ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "Lct6YiBIk1Bp", "outputId": "e7a5b2ec-4afb-47c4-9dde-f9ea38b00098" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "data15_periodD = 'simplecache::https://opendata.cern.ch/eos/opendata/atlas/rucio/opendata/ODEO_FEB2025_v0_exactly4lep_data15_periodD.exactly4lep.root'\n" ] } ], "source": [ "# We shall use the first entry in 'list', 'data15_periodD'\n", "data15_periodD = samples['Data']['list'][0]\n", "print(f\"{data15_periodD = }\")" ] }, { "cell_type": "markdown", "metadata": { "id": "xrVcQtUAk1Bp" }, "source": [ "Next, we shall try opening the `data15_periodD` file to see what is inside.\n", "In the file (called a `tree`),\n", " there are 39 entries,\n", " one for each event.\n", "In each event,\n", " a dictionary stores the all relevant information as keys, such as the event number (`eventNumber`), lepton transverse momentum (`lep_pt`), etc. \n", "Details on the variables in the dictionary can be viewed [here](https://opendata.atlas.cern/docs/data/for_education/13TeV25_details#variable-list).\n", "\n", "More information on trees can be viewed [here](https://masonproffitt.github.io/uproot-tutorial/03-trees/index.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "UfqvD5jQk1Bp" }, "outputs": [], "source": [ "# Accessing the file from the online database (\":analysis\" opens the tree in a desired manner)\n", "tree = uproot.open(data15_periodD + \":analysis\")\n", "\n", "# There are 39 entries in the tree\n", "print(tree.num_entries)\n", "\n", "# We can view all the information stored in the tree using the .keys() method.\n", "print(tree.keys())\n", "\n", "# We can also view the entire tree using the .arrays() method\n", "# This generates a 39-entry list of dictionaries\n", "print(tree.arrays())" ] }, { "cell_type": "markdown", "metadata": { "id": "bT6QTgd-k1Bq" }, "source": [ "Perhaps we'd like to see the lepton energies.\n", "We can access this from our tree using the key `lep_e`.\n", "Also,\n", " from this point on we shall be manipulating our tree arrays using the `awkward` library.\n", "We can use `library=\"ak\"` in the argument of the `.arrays()` method to use this library.\n", "If you ever see `library=\"ak\"` in the code,\n", " it means that the array is output as an `awkward` array." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "_Uxj9gnDk1Bq" }, "outputs": [], "source": [ "tree[\"lep_e\"].arrays(library=\"ak\")" ] }, { "cell_type": "markdown", "metadata": { "id": "qDqdjeSKk1Bq" }, "source": [ "In our analysis,\n", " not all the information in the tree is important.\n", "We can store the important variables in a list and retrieve them from the tree later on.\n", "As it turns out,\n", " we will need the following set of variables:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "id": "BU84o2F7k1Bq" }, "outputs": [], "source": [ "# Define what variables are important to our analysis\n", "variables = ['lep_pt','lep_eta','lep_phi','lep_e','lep_charge','lep_type','trigE','trigM','lep_isTrigMatched',\n", " 'lep_isLooseID','lep_isMediumID','lep_isLooseIso','lep_type']\n", "\n", "# To see all the data for our given variables\n", "# for data in tree.iterate(variables, library=\"ak\"):\n", "# print(data)" ] }, { "cell_type": "markdown", "metadata": { "id": "GlvwZpvik1Bq" }, "source": [ "Now that we understand how to access the information in the `data15_periodD` tree,\n", " we can begin analysis.\n", "As mentioned in the introduction,\n", " there are two key steps to be completed for each event entry:\n", "1. **Cuts** - we need to account for lepton selection rules in the event.\n", "In the [paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X),\n", " it is stated that we must\n", "\"[select] two pairs of isolated leptons, each of which is comprised of two leptons with the **same flavour** and **opposite charge**\".\n", "The datasets used in this notebook have already been filtered to include at least 4 leptons per event.\n", "We need to filter the data such that in each event, there are pairs of leptons of the **same lepton type** (`lep_type`) and summing to **zero lepton charge** (`lep_charge`).\n", "\n", "2. **Mass calculation** - the data to be plotted is the 4-lepton invariant mass, which can be found using the equation: $$m_\\text{4l} = \\sqrt{E^2_\\text{tot}-\\mathbf{p}_\\text{tot}\\cdot\\mathbf{p}_\\text{tot}}$$\n", "in units where $c=1$.\n", "$E_\\text{tot}$ is the total energy and $\\mathbf{p}_\\text{tot}$ is the total momentum.\n", "This calculation is performed using the vector array method `.M` on the sum of lepton 4-momenta (`lep_pt`,`lep_eta`,`lep_phi`,`lep_E`).\n", "\n", "From this,\n", " we can see why we chose those six important variables earlier.\n", "The physical reasoning for why we perform these steps is encapsulated in the idea of **conservation laws**.\n", "You may read more [here](https://cds.cern.ch/record/2759491/files/Conservation%20Laws%20-%20ATLAS%20Physics%20Cheat%20Sheet.pdf).\n", "\n", "Let's try to perform this two-step analysis for one event in `data15_periodD`." ] }, { "cell_type": "markdown", "metadata": { "id": "Kop2wRrmk1Bq" }, "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": 13, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "oiHAWTbyk1Bq", "outputId": "36b2cb9d-9251-41cb-b9f0-ee6a90423e72" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Cut for lepton type? [True]\n", "Cut for lepton charge? [False]\n", "Invariant mass: [5.6]\n" ] } ], "source": [ "# This selects the first entry of the tree\n", "entry = tree.arrays(library=\"ak\")[:1]\n", "\n", "# Cut lepton type (electron type is 11, muon type is 13)\n", "lep_type = entry['lep_type']\n", "sum_lep_type = lep_type[:, 0] + lep_type[:, 1] + lep_type[:, 2] + lep_type[:, 3]\n", "lep_type_cut_bool = (sum_lep_type != 44) & (sum_lep_type != 48) & (sum_lep_type != 52)\n", "print(f\"Cut for lepton type? {lep_type_cut_bool}\") # True means we should remove this entry (lepton type does not match)\n", "\n", "# Cut lepton charge\n", "# first lepton in each event is [:, 0], 2nd lepton is [:, 1] etc\n", "lep_charge = entry['lep_charge']\n", "sum_lep_charge = lep_charge[:, 0] + lep_charge[:, 1] + lep_charge[:, 2] + lep_charge[:, 3] != 0\n", "print(f\"Cut for lepton charge? {sum_lep_charge}\") # True means we should remove this entry (sum of lepton charges is not equal to 0)\n", "\n", "# Calculate invariant mass of the 4-lepton state\n", "# [:, i] selects the i-th lepton in each event\n", "p4 = vector.zip({\"pt\": entry['lep_pt'], \"eta\": entry['lep_eta'], \"phi\": entry['lep_phi'], \"E\": entry['lep_e']})\n", "invariant_mass = (p4[:, 0] + p4[:, 1] + p4[:, 2] + p4[:, 3]).M # .M calculates the invariant mass\n", "print(f\"Invariant mass: {invariant_mass}\")" ] }, { "cell_type": "markdown", "metadata": { "id": "byby5xFck1Bq" }, "source": [ "Based on our analysis, this entry should be removed because the lepton type does not match our requirements.\n", "We can turn these checks and calculations into a set of functions." ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "id": "ikYeGLvok1Br" }, "outputs": [], "source": [ "# Cut lepton type (electron type is 11, muon type is 13)\n", "def cut_lep_type(lep_type):\n", " sum_lep_type = lep_type[:, 0] + lep_type[:, 1] + lep_type[:, 2] + lep_type[:, 3]\n", " lep_type_cut_bool = (sum_lep_type != 44) & (sum_lep_type != 48) & (sum_lep_type != 52)\n", " return lep_type_cut_bool # True means we should remove this entry (lepton type does not match)\n", "\n", "# Cut lepton charge\n", "def cut_lep_charge(lep_charge):\n", " # first lepton in each event is [:, 0], 2nd lepton is [:, 1] etc\n", " sum_lep_charge = lep_charge[:, 0] + lep_charge[:, 1] + lep_charge[:, 2] + lep_charge[:, 3] != 0\n", " return sum_lep_charge # True means we should remove this entry (sum of lepton charges is not equal to 0)\n", "\n", "# Calculate invariant mass of the 4-lepton state\n", "# [:, i] selects the i-th lepton in each event\n", "def calc_mass(lep_pt, lep_eta, lep_phi, lep_e):\n", " p4 = vector.zip({\"pt\": lep_pt, \"eta\": lep_eta, \"phi\": lep_phi, \"E\": lep_e})\n", " invariant_mass = (p4[:, 0] + p4[:, 1] + p4[:, 2] + p4[:, 3]).M # .M calculates the invariant mass\n", " return invariant_mass\n", "\n", "\n", "def cut_trig_match(lep_trigmatch):\n", " trigmatch = lep_trigmatch\n", " cut1 = ak.sum(trigmatch, axis=1) >= 1\n", " return cut1\n", "\n", "def cut_trig(trigE,trigM):\n", " return trigE | trigM\n", "\n", "\n", "def ID_iso_cut(IDel,IDmu,isoel,isomu,pid):\n", " thispid = pid\n", " return (ak.sum(((thispid == 13) & IDmu & isomu) | ((thispid == 11) & IDel & isoel), axis=1) == 4)" ] }, { "cell_type": "markdown", "metadata": { "id": "9XuEr67ak1Br" }, "source": [ "You may verify on your own that these functions give the same outputs as the previous code block.\n", "Now,\n", " we shall apply these functions over the entire data tree using a `for` loop." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "id": "IqovJmoGk1Br" }, "outputs": [], "source": [ "# Define empty list to hold all data for this sample\n", "sample_data = []\n", "\n", "# Perform the cuts for each data entry in the tree\n", "for data in tree.iterate(variables, library=\"ak\"): # the data will be in the form of an awkward array\n", " # We can use data[~boolean] to remove entries from the data set\n", " lep_type = data['lep_type']\n", " data = data[~cut_lep_type(lep_type)]\n", " lep_charge = data['lep_charge']\n", " data = data[~cut_lep_charge(lep_charge)]\n", "\n", " data['mass'] = calc_mass(data['lep_pt'], data['lep_eta'], data['lep_phi'], data['lep_e'])\n", "\n", " # Append data to the whole sample data list\n", " sample_data.append(data)\n", "\n", "# turns sample_data back into an awkward array\n", "data_A = ak.concatenate(sample_data)" ] }, { "cell_type": "markdown", "metadata": { "id": "wcQ84YFzk1Bs" }, "source": [ "We can now plot the data using Matplotlib.\n", "The data will be turned into a histogram,\n", " with bins of width $2.5 \\,\\text{GeV}$.\n", "Note that much of the code written here is meant for the aesthetics of the plot." ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 453 }, "id": "6XrVtATfk1Bs", "outputId": "f459e810-879d-40db-a89d-eb523c43e8e9" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEKCAYAAADzQPVvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnkElEQVR4nO3de7xcZX3v8c+P3GBDCJdsaESSHRVjxQJmj0q9kQJyENxy2toKRqtW2Uct3vAlBxtOKec0rdfmtOWonXDxki2WikcT+wJBS7BWse5QUJA7JBECsoOHSAjkAr/zx1qzWdnMmlkz6zp7vu/Xa157Zq1nPc9vPXvmmTVrnvktc3dERGR626fsAEREJH8a7EVE+oAGexGRPqDBXkSkD2iwFxHpAxrsRUT6wMy4FWY2y913FxlMwwEHHOAveclL2pbbtm0b8+bNK7zcxMQEg4ODhbZb5diyjq/KsZUVX5Vjyzq+KsdWVnxJ69qwYcNv3L15QXdvegMeAS4BTgIsrlwetzlz5vjZZ5/ta9eu9VbOPvvsluvzKjc8PFx4u1WOrZNySeKrcmx5tJukXJVjc882virHlke5LF4Ta9eu9bPPPtuBbR4zrrY6jfPbwE+BC4Bfmtnfmdnxbd9aMjBjxgzq9TojIyMty7Vbn1e5pLJst8qxdVKujDbVd/mXSyrL57r+r3uvr9frEBykNxf3LhC9Ac8DPgz8GLgXWJlku25vAwMD7d8OS5T0SKEMVY7NvdrxKbbuVTm+Ksfmnm18wLh3cWQffUPYAlwKfAF4HHhvku26NX/+/DyrT210dLTsEGJVOTaodnyKrXtVjq/KsUFx8Zm3yI1jZvsCI8BZwKuBa4CvA9e5+9N5BVWr1Xx8fDyv6kVEpiUz2+DutWbrWs3G+RpwMnADMAa8zd2fyidEERHJU6vTONcAL3T3P3L3q4oc6Ldt28bo6Cjr1q0rqkkRkZ61bt26xumg2PmZLU/jAJjZ4cBfA0e4+6lm9lLgd9390iyDjdJpHBEpwrJlywBYv359qXFkpdVpnCRf0H4J+C6wIHx8F/CRTCITEZlGZsyYwXHHHcfRRx/Nsccey+c+9zmeeeaZltts3LiRr33ta7nHlmSwn+/uVwLPALj7HiC3L2dFRIowNjbGjTfeyA033MDQ0BBjY2Op69xvv/24+eabue2227juuuu4+uqrueiii1puU6XB/gkzOxRwgPCHVdtyjUpEJEdjY2OMjo6yc+dOADZt2sTo6GgmA37DYYcdRr1e5+KLL8bd2bhxI6973etYunQpS5cu5Uc/+hEA559/Pv/2b//Gcccdx6pVq2LLpRY3Ab9xA5YC/04wwP87wWmcY9ptl+ZW9R9BiEhvW7RokRMcwO51W7RoUap6999//+csmzdvnj/88MP+xBNP+JNPPunu7nfdddfkj6muv/56P/300yfLx5VLghY/qoqdehl5M7jJzE4AlgAG3Ok5J0hrzMYZGRnJ/CfRIiKbN2/uaHkWdu/ezTnnnMPNN9/MjBkzuOuuu1KVi1q3bl1j9mLsbJxW8+wPBA5397vdfU84C2c/YKmZfdfdf9U2gi7NmzevkedBRCRzCxcuZNOmTU2XZ+m+++5jxowZHHbYYVx00UUcfvjh3HLLLTzzzDPsu+++TbdZtWpVonJRjQPj1atXx55ib3XO/rPAayKP/waoAa8HWn/jICJSYStXrmRgYGCvZQMDA6xcuTKzNiYmJnjf+97HOeecg5mxbds2FixYwD777MNXv/pVnn46mOcyd+5cHn/88cnt4sql1WqwfwXw5cjjx939Q+7+XuBl7So2s8vM7BEzu3XK8g+a2R1mdpuZfbq7sEVEurd8+XLq9Tpz5swBYNGiRdTrdZYvX56q3ieffHJy6uXJJ5/MKaecwoUXXgjABz7wAb785S9z7LHHcscdd7D//vsDcMwxxzBjxgyOPfZYVq1aFVsurdgfVZnZz939dyKPX+but4b3b3X3lgO+mb0e2A58pVHWzH4PWAGc7u47zewwd39OSk79qEpEitBPP6pq9QXtM2b2W+7+MEBkoD+CcM59K+7+AzMbmrL4/cAn3X1nWCY+97KISM6myyCfRKvTOJ8B1pnZ681sbng7AfhWuK4bLwZeZ2Y/MbMbzOwVXdYjIiIdiD2yd/c1ZrYV+CvgaIJ5qLcBf+HuV6do7xDgeILvBK40sxf4lHNJExMT1GrPfhIZHR2tfE5qEZEy1Ov16OzF2IuBtE2ElkZ4Guc7kXP21wCfcvfrw8f3Ase7+0R0O52zFxHpXNpEaFn6FvB7AGb2YmA2sLXgGERE+k7bX9B2y8yuAJYB883sAeBC4DLgsnA65i7gnVNP4YiISPZyG+zd/ayYVW/Pq00REWku8WBvZq8FXgnc6u7X5heSiIhkLfacvZn9R+T+2cDFwFzgQjM7P8+gdFlCEZHkUl2W0Mz+091fHt7/KXCau0+Y2f7AjdFf12ZNs3FERDrX7S9o9zGzgwmO/q0xPdLdnzCzPTnEKSIiOWk12M8DNhDksHczW+DuD5nZAeEyERHpEa1+QTsUs+oZ4PdziUZERHLR8dRLd98B3J9DLCIikpOufkFrZt/JOpAozcYREUku1WycVhrn71PE1pJm44iIdC7z3Dh5DvQiIpK9Vj+qOtDM/sbMvmpmb5uy7vP5hyYiIllpdWR/OcEUy6uAM83sKjObE647PvfIREQkM60G+xe6+/nu/i13fzNwE/CvZnZoQbGJiEhGWg32c8xscr27rwRWAz8Ach3wNRtHRCS5tLlxPg1c6+7fm7L8VOAf3P2oDGPdi2bjiIh0rqvcOO5+Xszya4DcBvqiLVu2DGh9lfkkZfLYVkQkK7ldltDMLjOzR8KrUk1d9zEzczOLvTiuiIhkJ89r0H4JOHXqQjM7EjgF2Jxj2yIiEtFqnv3z0lTs7j8Aft1k1SrgPEDXnhURKUirRGiXmNkhwHrgGuCH7p4qj72ZnQE86O63mMVnSZ6YmKBWe/Y7htHR0cY3zSIiElGv16nX642HsafGW31Be5qZ7QssI0hp/Fkz20ww8F/j7h2dhjGzAeDPCU7htDQ4OIhm44iItBc9GDazrXHlWqY4dvenCAf3sKLFwBuBi83st9z9lR3E9EJgMdA4qn8+cJOZvdLdH+6gHhER6VBH+ezd/X7g88DnzWx2h9v+HDis8djMNgI1d499JxIRkWx0PRvH3Xe1Wm9mVwA/BpaY2QNm9p5u2xIRkXQ6vlJVUu5+Vpv1Q3m1LSIie+voyN7MDjazY/IKRkRE8tF2sDez9WFu+0MIMl+uNrO/zTMoJUITEUkuSSK0JEf289z9N8AfAF9x91cBJ2cTYkyD8+ZRr9cZGRnJsxkRkWlhZGSkMdd+W1yZJIP9TDNbAPwxkOuFxkVEJB9JBvuLgO8C97j7T83sBcDd+YYlIiJZSjIb5yF3n/xS1t3vy/ucfV6UblhE+lWSI/t/SLhMREQqKvbI3sx+F3g1MGhm50ZWHQjMyDOoxmyckZERfUkrItLGunXrGrMXY2fjtDqNMxs4ICwzN7L8N8BbsggwTmM2joiItNc4MF69enXsbJxWWS9vAG4wsy+5+6ZcIhQRkUIk+YJ2jpnVgaFoeXc/Ma+gREQkW0kG+38GvghcAjydbzgiIpKHJIP9Hnf/Qu6RiIhIbpJMvVxnZh8wswVmdkjjlmdQyo0jIpJcktw4SY7s3xn+/XhkmQMv6D601jQbR0QkuVSzcRrcfXE3jZvZZcCbgEfc/WXhss8AI8Au4F7g3e7+WDf1i4hIcklSHA+Y2QXhjBzM7Cgze1OCur8EnDpl2XXAy8L0C3cBn+gwXhER6UKSc/aXExyJvzp8/CDwV+02cvcfAL+esuxad98TPryR4KLjIiKSsySD/Qvd/dPAbgB33wFYBm3/KXB1BvWIiEgbSQb7XWa2H8GXspjZC4GdaRo1sxXAHmCs2fqJiQlqtRpz585l7ty5lfyydtmyZZNZNEVEylKv16nVatRqNYD5ceWSzMb5S+Aa4EgzGwNeA7yr28DM7F0EX9ye5O7erMzg4CDj4+OTg2k4pUhERKYYHR2dHCPNbGtcuSSzca41sw3A8QSnbz7s7rEVtmJmpwLnASeEp4NERKQAbQd7M1sHfA1Y6+5PJK3YzK4AlgHzzewB4EKC2TdzgOvMDOBGd39fF3GLiEgHkpzG+SzwVuCTZvZT4OvAd9z9qVYbuftZTRZf2nmIIiKSVpLTOI1UxzOAE4GzgcsILmIiIiI9IMmRPeFsnBGCI/ylwJfzDKqRG+fRRx/l0EMPzbMpEZGel+RKVUl+QXslcDvBUf3FBPPuP5hVkM00cuNooBcRaW9kZKQxRb373DgE59nPcnflshcR6VGxR/Zmdh6Au38X+IMp6/4657gKMTY2xo033sgNN9zA0NAQY2PP/Y1XkjJp6hcRKUKr0zhnRu5PTVg2NcFZzxkbG2N0dJSdO4MfA2/atInR0dG9BuQkZdLULyJSlFaDvcXcb/a456xYsYIdO/b+XdeOHTtYsWJFR2XS1C8iUpRWg73H3G/2OFPR2Th52bx5c9vlScqkqV9EJAtJrlTVarA/1sx+Y2aPA8eE9xuPfyfjWPdSxGychQsXtl2epEya+kVEspBkNk7sYO/uM9z9QHef6+4zw/uNx7PyCLhIK1euZGBgYK9lAwMDrFy5sqMyaeoXESlKkhTHlddNuuHly5dTr9eZM2cOAIsWLaJer7N8+fKOyqSpX0SkKIl+QTtdLV++nNWrVwOwfv36rsukqV9EpAjT4sheRERa02AvItIHKnkaR4nQRESSyyQRWhmUCE1EJLlUUy/TMrPLzOwRM7s1suwQM7vOzO4O/x6cV/siIvKsPI/sv8Rzc+icD3zf3Y8Cvh8+LoSSkolIP8vtnL27/8DMhqYsPoPgurQQXABlPfDf84qhIS4pmYhIvyj6nP3h7v5QeP9h4PAiGlVSMhHpd6V9QevuTkxCtYmJCWq1Ghs2bGDDhg2NLx66pqRkIjJd1et1arUatVoNYH5cuaIH+1+Z2QKA8O8jzQoNDg4yPj7O8PAww8PDqU+5KCmZiExXo6OjjI+PMz4+DrA1rlzRg/1a4J3h/XcC3y6iUSUlE5F+l+fUyyuAHwNLzOwBM3sP8EngDWZ2N3By+Dh3SkomIv0uz9k4Z8WsOimvNqMaWTAbCciqlJRsamxVUuXYRKR7lfwFrYiIZEu5cUREepxy44iI9IFSc+OIiEh1VHawr3IumyrHJiLSTCUH+1//+tdNc9lUYVCNy7NThdhEROJUcrB/8MEHK5vLRnl2RKQXVXI2zq5du5our0IuG+XZEZGq6dnZOLNnz266vAq5bJRnR0Sqpmdn4xxxxBGVzWWjPDsi0osqOdgfcsghlc1lozw7ItKLKnnOHqqVy2aqKscmItJMJY/sRUQkW5U8slduHBGR5Hp2Nk6ZuXGWLVs2mea3yG2LaKOI+KrQpki/6dnZOCIikq1SBnsz+6iZ3WZmt5rZFWa2b7d1KU+NiEh7hZ+zN7MjgA8BL3X3J83sSuBM4Eud1hWXp0ZERPZW1mmcmcB+ZjYTGAC2dFOJ8tSIiCRT+GDv7g8CnwU2Aw8B29z92miZiYkJarUaGzZsYMOGDY0vHp5DeWpEpN/V63VqtRq1Wg1gfly5wgd7MzsYOANYDDwP2N/M3h4tMzg4yPj4OMPDwwwPD8eemlGeGhHpd6Ojo4yPjzM+Pg6wNa5cGadxTgbud/cJd98NfBN4dTcVKU+NiEgyZQz2m4HjzWzAzAw4Cbi9m4qUp0ZEJJnCZ+O4+0/M7BvATcAe4D+B5iflE0ibpybJNmny3yh3johUQSnpEtz9QuDCMtoWEelH+gWtiEgfUCI0EZEep0RoIiJ9QInQREQE6MHBPknK3ConR4uLrQqpgKvcb3mpQr+LFKGS5+zTqHJytF6NTb9bEOl9PXdk306Vk6MpNhEpSyUH++hsnE5VOTmaYhORPKxbt67xSbx/ZuNUOTmaYhORPPTlbJxWydHWr19favqCKiduq3JsIpLetBvsq5wcTbGJSFmm3WwcSJ8cLU+KTUTKMO2O7EVE5LkqeWSv3DgiIskpN46ISB/oy9k4IiLyXKUM9mZ2kJl9w8zuMLPbzex3y4hjqjS5YYrIK1P1+KrQpog0V9Y5+78DrnH3t5jZbGCg3QZ5S5Mbpoi8MlWPrwptiki8wo/szWwe8HrgUgB33+XujxUdx1RpcsMUkVem6vFVoU0RiVfGaZzFwARwuZn9p5ldYmb7RwtMTExQq9XYvn0727dvb3zxkKtOc8NEf42bV16ZaPrdNG2UkfemVZtVTitc5dhEmqnX69RqNWq1GsD8uHJlDPYzgaXAF9z95cATwPnRAoODg4yPj0/eikgDnCY3TBF5ZaoeXxXaFOlHo6Ojk2MlsDWuXBmD/QPAA+7+k/DxNwgG/1KlyQ1TRF6ZqsdXhTZFJF7hg727Pwz80syWhItOAn5RdBxTpckNU0RemarHV4U2RSReWbNxPgiMhTNx7gPeXVIce0mTG6aIvDJVj68KbYpIc6UM9u5+M1Aro20RkX5UyV/QNnLjhLkeRESkhSRXqqpkIrRGbhwREWlvZGSEkZERVq9erdw4IiL9TIO9iEgf6KnBPi6xVtnXlk1LCcPKoX6XflLJc/bNTNfEWq32S/IT1+8LFy7k8MMPLzk6kez1zJH9dE2sNV33q+ri+v3+++8vKSKRfFVysG829bKMZF5FmK77VXVx/ds40hfpJUmmXlZysG9MvRwZGZlcNl0Ta03X/aq6uP5tpHcQ6SXT6rKE0zWx1nTdr6qL6/fFixeXFJFIvnpmsJ+uibXK3q8k+durlOM9aSztysX1e7MvZ5vVlXWfFNHHWbZRpeeEJNMzs3Fg+ibWmq77VXXN+r3xWGS66ZkjexER6V4lB3slQhMRSU6J0ERE+kCSRGilDfZmNgMYBx509zdlXX+Vz30rNhEpWpmncT4M3F5i+5U2XfO2VHm/qhybSFqlHNmb2fOB04GVwLllxFBl0zVfTpXzG03XPhdpKOvI/n8D5wHPlNR+pU3XfDlV3q8qxyaShcIHezN7E/CIu2+IKzMxMUGtVpu89duXtdM1X06V96vKsYm0Uq/XJ8dKYH5cuTJO47wGeLOZnQbsCxxoZmvc/e2NAoODg4yPj5cQWjUsXLiQTZs2NV3ey6q8X1WOTaSV0dHRyVOOZrY1rlzhR/bu/gl3f767DwFnAv8aHehl+ubLqfJ+VTk2kSxU8kdV/a7sfDl5qfJ+VTk2kSyU+qMqd18PrC8zhqqarvlyqrxfVY5NJK1K/oK2TGle5EUMEFWPrwptishzVfI0jnLjiIgkp9w4faKRV7zso+g0cVRlH9JIug+91k9Zt9msvrL6bjo87yBZbpxKHtmLiEi2NNj3sSS5YKqULyZpLFnG3KyurPukiD7Ou0+k+ip5GkfylyRPTZVy2SSNJcuYm9X17ne/GzNj165dz6k/z/1KI+8+qUp+I2lNR/Z9KkkumCrli0kaS5YxN6tr9+7dkwN92vrj2si6j/PuE+UQ6g2VHOw1Gyd/SXLBVClfTNJYsoy5k2267ZMi+riIPlEOoXIlmY1TycG+MRtnZGSk7FCmrbicL9HlScoUJWksWcbcyTbd9kkRfVxEnyiHULlGRkYaMxg1G0f2liQXTJXyxSSNJcuYm9U1a9YsZs+enUn9cW1k3cd594lyCPWGnhvs169f3/NzYpsper+S5IKpUr6YpLF0GnOrfm9W1+WXX85ll12WWZ8U0cdZtlGl54R0RrNx+liSXDBVyheTNJYsY46rK8s+KaKPi+gTqbaeO7IXEZHOVXKw12wcEZHklBtHRKQPKDeOiIgA5Vxw/Egzu97MfmFmt5nZh4uOQUSk35RxZL8H+Ji7vxQ4HvgzM3tpCXFMC1VJSpUmjqrsQxpFJGkro5+KSPpWVt9Nh+ddR9y91BvwbeAN0WXDw8Mu7ieccIKfcMIJsevXrFnjAwMDDkzeBgYGfM2aNZm1kUcc0frSbJt2v5Ju125Z0n1oVa4q/+us4k1a36xZs3z27Nmp+i7PdnsNMO5xY23ciiJuwBCwGTgwulyDfaDdC2rRokV7PVkbt0WLFmXWRh5xROtLs23a/Uq6XbtlSfehVbmq/K+zireT+tL2XZ7t9ppWg31ps3HM7ADgKuAj7v6b6LqJiQlqtdrk49HR0VRpZHtVux+slJmUKnqFn07jiO5Xq22zvnpRknLd/Ego6T60Kjc0NPSc5Um3zevKV2nibbas08RyaZ5jadrNUhFX1lqyZAlbtmxhyZIlAPPj6itlNo6ZzSIY6Mfc/ZtT1w8ODjI+Pj5568eBPomqJKVKE0dV9iGNIpK0ldFPWbeZJrFcHn2XtmxVLFiwgOHhYcbHxwG2xpUrYzaOAZcCt7v73xbd/nRSlaRUaeKoyj6kUUSStjL6Kes20ySWy7rvsk5o1wvKOLJ/DfAO4EQzuzm8nVZCHD2vKkmp0sRRlX1II68kbVlt262s20yTWC7rvss6oV0vKPycvbv/ELCi252uqpKUKk0cVdmHNIpI0lZGP2XdZprEcnn0Xa8/7zqhX9CKiPSBSg72SoQmIpLco48+qssSiohMd4ceeqguSygiIpSfLqHZTb+g7Uy3vxbt1po1a3zOnDmTvzhs/MQ8TRxTt41rI2ks3ZZLI+k+JOmnpNum2a+k23Ybbyf7X0bfJa2vW1n/b5IsA+7zKqZLiLtpsO9MkYN91vlSGrLONdNNubSS7kOaPDjd9FOnbbTar6R1tco9027QLaLvku5bt4rI5dNsGfC0a7Cfvooc7LPOl9KQda6ZbsqllXQf0uTB6aafOm2j1X51Ulcn+1903yXdt24Vlcun2c1jxtVKnrPXbJzqKiIfT9I2si6XpTRtFrFfWfZJ1rlnqrJfaeSxD2lVcrDXbJzqKiJHS9a5Znotr0yv5drJOvdMVfYrjaJy+XSikoO9VFcROVqyzjXTa3llei3XTta5Z6qyX2kUkcun2TLgmdiK487vlHnTOfvOFHnO3l2zcdLuQ5YzSjQbJ/l+FXnOvlUc3W6bdjaOuXvbd5qi1Wo1D9N1SgJpcmZn2WYRubuTbpumXLeS7kOaeIvo927KdBJb1suyjiVLRbwmosvMbIO7P3sxkIjSLl4i2alKEqc0cVRlH9JIug+91k9Zt9msvrL6bjo875Kq5Dl7zcYREUlOuXFERPqAcuOIiAhQ3jVoTzWzO83sHjM7f+r6iYmJMsJKLHwHraR6vc7Y2BhDQ0Pss88+DA0NMTY2VnZYk6red1XVSWzd/v/TPG/K6rskMSeNrVldaZYlbeM973lPMa/XuGk6ed2AGcC9wAuA2cAtwEujZQYGBhJNT1q7dm0p5ZJODc2y3aR1DQ0NJcrJkTa2qVPWsuy7Rl3tpsVlXS6P/2uSqX1J+rhVbNFy5557bkf//8a2cblczj333ET7Ojw8nGhqYyfP9XZ912pfk/ZdQ3QKY+MWl49m5syZbcs1+i4aR1zOGzNruw+t+iS6DLjbK5Qu4ZXAPe5+n7vvAr4OnNFNRUm/wM26XFJZtpu0ri1btrBjx469lu3YsYMVK1bkFlsn5cpos8r/107KJVGv17v6/69YsaLpdlkfsWf5XE+6r0msWLGCnTt37rVs9+7d7Nq16znL9uzZ07Zcs75r1se7d+9uHAS33IcOniMHxa0ofJ69mb0FONXd3xs+fgfwKnc/J1LmKeDpyGYTwNYm1c2jxRcSOZabHxNPnu0mrWu4xboNXdRXRt9VObY82k1SLmls3f7/k24XJ8u+y/q5niS2VnWlEY2j0zaSvl7nA4Ph/X3cfb9mhSo5z97d9y07BhGR6aSM0zgPAkdGHj8/XCYiIjkpY7D/KXCUmS02s9nAmcDaEuIQEekbhQ/27r4HOAf4LnA7cCVwipndZma3mtkVZrZv+Gbwk3B65j+FbwyFMLPLzOwRM7s1suwQM7vOzO4O/x4cLjcz+/swzp+Z2dISYvuMmd0Rtv9/zeygyLpPhLHdaWb/pejYIus+ZmZuZvPDx4X2W6v4zOyDYf/dZmafjiwvte/M7Dgzu9HMbjazcTN7Zbi86OfckWZ2vZn9IuyjD4fLS39NtIitKq+JpvFF1hf3uoibplPUDTgCuB/YL3x8JfCu8O+Z4bIvAu8vMKbXA0uBWyPLPg2cH94/H/hUeP804GrAgOOBn5QQ2ynAzPD+pyKxvZRgauscYDHBlNcZRcYWLj+S4M19EzC/jH5r0Xe/B3wPmBM+PqwqfQdcC7wx0l/rS3rOLQCWhvfnAneF/VP6a6JFbFV5TTSNL3xc6OuiKr+gnQnsZ2YzgQHgIeBE4Bvh+i8D/7WoYNz9B8Cvpyw+I4xjajxnAF/xwI3AQWa2oMjY3P1aDz4xAdxI8D1II7avu/tOd78fuIdg6mthsYVWAecRzCNuKLTfWsT3fuCT7r4zLPNIJL6y+86BA8P784AtkdiKfM495O43hfcfJ/hEfgQVeE3ExVah10Rc30HBr4vSB3t3fxD4LLCZYJDfRjDl6LHIP+sBnu2gshzu7g+F9x8GDg/vHwH8MlKu7Fj/lODIACoQm5mdATzo7rdMWVV6bKEXA68LTxneYGavCJdXIb6PAJ8xs18SvEY+ES4vLTYzGwJeDvyEir0mpsQWVYnXRDS+Ml4XpQ/24Xm+Mwg+Uj0P2B84tdSg2vDg81blLgRgZiuAPUAl8iOY2QDw58BflB1LCzOBQwg+Mn8cuNLMrNyQJr0f+Ki7Hwl8FLi0zGDM7ADgKuAj7v6b6LqyXxNxsVXlNRGNL4yn8NdF6YM9cDJwv7tPuPtu4JvAawg+vjR+B1CF6Zm/anycCv82Pu5XYiqpmb0LeBOwPHzhQfmxvZDgTfwWM9sYtn+Tmf1WBWJreAD4Zvix+T8ILus2vyLxvZPg9QDwzzx7uqHw2MxsFsFgNebujZgq8ZqIia0yr4km8ZXyuqjCYL8ZON7MBsIjqpOAXwDXA28Jy7wT+HZJ8TWsDeOAveNZC/xJ+C368cC2yEfbQpjZqQTn/t7s7tHfY68FzjSzOWa2GDgK+I+i4nL3n7v7Ye4+5O5DBAPrUnd/mAr0W+hbBF/SYmYvJsjXtJWS+y60BTghvH8icHd4v9C+C1+XlwK3u/vfRlaV/pqIi60qr4lm8ZX2usjqm940N+Ai4A7gVuCrBN+Uv4Dgn3APwVHNnALjuYLg+4Pd4T/iPcChwPcJXnDfAw4Jyxrwfwi+1f85UCshtnsIzvPdHN6+GCm/IoztTsKZHUXGNmX9Rp6ddVBov7Xou9nAmvC5dxNwYlX6DngtwfdXtxCchx4u6Tn3WoJTND+LPMdOq8JrokVsVXlNNI2vjNdFJa9BKyIi2arCaRwREcmZBnsRkT6gwV5EpA9osBcR6QMa7EVE+oAGexGRPqDBXkSkD2iwl0oxs/VmdkHZcXTKzK42s/P6PYZ+Y2YbzewpM7utpPavNbMnzWxPu7Ia7KUQZraPmf0ovFDD89tvkVm7hbx5uPsb3f3T7Uv2Rgy9+qZbkve6+9HRBWY2bGZXWXBBmu3hm8JVZnZikgrN7Ntm9pWYddeb2cUA7n4K8MYkdWqwl6J8FNjRtpR0LEy0JRVhZm8A/p0g5UGN4KIlvwN8Dfj9hNX8I/AWi1xhK6z7KIJ8Sf/YcWB55oXQTTd3hyBn/L3AcQR5Qp7foux64ILI4wGCXO73E1zc4xrgRZH1GwlSxf4Q2A6MA68I110MPA3sDNfdGS4/FPgKQQ72hwkuvHHIlDr/nCDvy3aCvDmvbrOPk3G32x74M+DmKdsvDmMdAj5MkCvqcYJEgX9D5GpKkX2+Pqz/zCYxJKmjaYxx/Raz3xuBCyKx/Bw4BjiLID/NNuASnr1qVLu4PhT+rx8nyPb410nWpY0rxXN7I/D2KcvuAS5JsG3sc5vgQHwT8MEp23wG+PGUZcuAPW3bK3sg0G1638In7Q8JjmiG6HywHwO+Q3BhjNk8mzRvVrh+I0F2yOFw/fnABHBgs/rCZdcA64CDw9u/AP8SWb8xfMEeDcwguKLQ3W32c7KddtuHbT4FHBdZdhHw/fD+HxIM/kZwsYtfAf9tSny/DNcZz17SMxpDkjpaxficfovZ740EidB+G5hFkFTuXqBOcG2KhQSpj5e3i4vgoGAHcHT4+CDg+HbrsogrxfN7I5HBPozTgZMSbNvuuf0/gJ9Fys8OY37XlHqWocFet7JvBKdvvhHeH6KDwZ4gr7wDCyPr9yE4Kntt+Hgj8L8i643giPFtU+sLHz8vrPOoyLIl4bIFkTo/Hll/dLh+XsK4224P/BPwd5GYN8YNPARHf1dGHm8E/qJVDAnriI2xVV1T6p1az2lhPYORZVcCq9rFRZDp9kngj4EDppSLXZdlXASf+v4f4QBO8Jz9Xpt2ooP9a8J2XhJZ9mbgMYLn7VMdPLefR5AF9VXh47eGse03JYZlJBjsdc5ecmNmLwI+BpwTs355+OXVdjPb3qTI4vDvz8zsMTN7jODj7iz2vsDDxsYdD579m3n2mqNTNba7P7Ls3inrIEg33PBE+HduTJ3NtNv+cuBt4fn2EwmOVL8JYGZnmdlPzexRM9tGcNpncEr9G1s1nrCOtPvYrJ4dwNPuPjFl2dx2cbn7fcBy4Gxgi5n90MxOabcui7giLiD4JNqtreHfyeefu69194OA0wnSt0OC57a7byE48h8Ny44Ca9z9yW4C02AveXotwQv5VjPbSpAvHoIn+AfcfczdD2jcmmy/Kfx7lLsfFLkNuPsVkXJDjTvhxSIWEuSEh+DKU1G/nLoNwVFjdF0RriM4Jz4CvIvgIthPmtmRBKcc/orgk8Y8gvzmUy+VOHW/JnVQRyux9XcrSVzu/k13fwPBke+VwLctuLxly3UZxfcigiP7DSmquQu4DzizTbmkz+068FYzeznBRXY6/2I2pMFe8nQlwSXYjgtvp4XLTyH4grQld3+EYAbD583sCAAzO8jMft+Ca3o2/KmZLQ2Pkj9O8MXXv4TrHgZeFKlzC3At8LmwroOBzwFXe4FXynL3pwn64EPAHwCXhasOIHhdTgC7w6sVvaPD6rOoY69+y0jLuMxsiZmdGg7guwlOaTjwTKt1Gcb3PwnOm3ct/GT5Z8A7zOxTZnZkeNWpAeBVkXJJn9vfJfi0cBXBF7O3dhubBnvJjbvvcPcHGjeCAQTgYXdvdtqmmbMJrii03sweJ5hV8UfsfXHrOvD3BOcz3wqc7u7bwnWrgFr4Ubnxw5e3E8zouJPgC7HHgD/pZh9TupxgGt39Hlz/Fne/HbiQ4BJ/jxF84XxFXAXNZFEHzfstlQRxzSaYZfRQuP5DwB+6+1Nt1qVmZq8GHnX3e9sWbsPdryH4VPtigk+z24HbCM7nR+fZt31uu/szwGqC0z71NHHpSlXS0yy4YPMF7r6m7Fikd5nZhwhmCj1J8InmCeB9BG8ul7j7yTHb3QksADa6+zEFhRtt/2qCN5F9Yk6FTppZTEgiItXl7n9P8OkQM/tL4B53/7GZDbXZbkn+0bVsP9GvZ0GDvYjIXtz9LyP3NwJNj+p7jU7jiIj0AX1BKyLSBzTYi4j0AQ32IiJ9QIO9iEgf0GAvItIHNNiLiPQBDfYiIn1Ag72ISB/4/ygNXh933QsuAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# x-axis range of the plot\n", "xmin = 80 * GeV\n", "xmax = 250 * GeV\n", "\n", "# Histogram bin setup\n", "step_size = 2.5 * GeV\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", "# Creating histogram from data\n", "data_x,_ = np.histogram(ak.to_numpy(data_A['mass']),\n", " bins=bin_edges ) # histogram the data\n", "data_x_errors = np.sqrt( data_x ) # statistical error on the data\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", "# 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", "# draw the legend\n", "main_axes.legend( frameon=False ); # no box around the legend" ] }, { "cell_type": "markdown", "metadata": { "id": "KypqMZFlk1Bs" }, "source": [ "Great,\n", " we managed to plot `data15_periodD`!\n", "Now,\n", " we have not discussed how to deal with the Monte-Carlo simulation data,\n", " or even what they are for.\n", "Let us explain." ] }, { "cell_type": "markdown", "metadata": { "id": "2pMKl0vRk1Bs" }, "source": [ "## Example 2: Reading Monte-Carlo data" ] }, { "cell_type": "markdown", "metadata": { "id": "m2LnGyIWk1Bs" }, "source": [ "Using the Standard Model, we can do a set of randomised simulations to produce a set of theoretical data points to compare to our ATLAS data. These are known as Monte-Carlo(MC) simulations. But there is one important change to be made to the MC data before we can compare them with our ATLAS data:\n", "\n", " - **Weights** - The MC data was simulated with *ideal* or *controlled* circumstances. However, the real ATLAS detector has some inefficiencies, which we can account for by attributing an appropriate weight to each data point. The weight of a data point affects how it contributes to the histogram count for its bin.\n", "\n", "Let's open an MC file." ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "id": "wfpSfVpek1Bs" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "simplecache::https://opendata.cern.ch/eos/opendata/atlas/rucio/opendata/ODEO_FEB2025_v0_exactly4lep_mc_410470.PhPy8EG_A14_ttbar_hdamp258p75_nonallhad.exactly4lep.root\n" ] } ], "source": [ "# We open an MC data file with simulated ttbar events (data set ID 410470)\n", "value = samples[r'Background $Z,t\\bar{t},t\\bar{t}+V,VVV$'][\"list\"][0]\n", "print(value)\n", "\n", "# This is now appended to our file path to retrieve the root file\n", "background_ttbar_path = value\n", "\n", "# Accessing the file from the online database\n", "tree = uproot.open(background_ttbar_path + \":analysis\")" ] }, { "cell_type": "markdown", "metadata": { "id": "xYiXBOBPk1Bs" }, "source": [ "These are the weights which are important to our analysis:\n", "- `scaleFactor_PILEUP` - scale factor for pileup reweighting\n", "- `scaleFactor_ELE` - scale factor for electron efficiency\n", "- `scaleFactor_MUON`- scale factor for muon efficiency\n", "- `scaleFactor_LepTRIGGER` - scale factor for lepton triggers\n", "\n", "Scale factors are generally related to estimates of the efficiencies and resolutions of detectors." ] }, { "cell_type": "code", "execution_count": 18, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "v0Ae8IQSk1Bs", "outputId": "75b8f479-c906-465a-d9a3-01ddef2ae4f0" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['filteff', 'kfac', 'xsec', 'mcWeight', 'ScaleFactor_PILEUP', 'ScaleFactor_ELE', 'ScaleFactor_MUON', 'ScaleFactor_LepTRIGGER']\n" ] } ], "source": [ "weight_variables = [\"filteff\",\"kfac\",\"xsec\",\"mcWeight\",\"ScaleFactor_PILEUP\", \"ScaleFactor_ELE\", \"ScaleFactor_MUON\", \"ScaleFactor_LepTRIGGER\"]\n", "\n", "# For example, see below for the weights corresponding to muon rejection\n", "tree[\"ScaleFactor_MUON\"].arrays(library = \"ak\")\n", "print(weight_variables)" ] }, { "cell_type": "markdown", "metadata": { "id": "YUGgtOInk1Bs" }, "source": [ "Additionally, there is a *cross section weight* $w_\\sigma$ associated with each MC file. This weight is meant to normalise the entire Monte-Carlo distribution based on the number of events in the data. This is its definition:\n", "$$ w_\\sigma = \\frac{\\int L \\text{d}t ~ \\sigma \\,\\, \\eta \\,\\, k \\,\\, w}{\\sum_i w_i } $$\n", "where \n", " * $\\int L \\text{d}t$ is the integrated luminosity in inverse femtobarns (`lumi`),\n", " * $\\sigma$ is the cross section in picobarn (hence the factor $1000$ in the formula below) (`event[\"xsec\"]`),\n", " * $\\eta$ is the filter efficiency of the MC generator (`event[\"filteff\"]`), \n", " * $k$ is the k-factor, a multiplicative correction factor used to account for higher-order effects in theoretical calculations of the cross section (`event[\"kfac\"]`). \n", " * $w$ are weights added to a given simulated event by the event generator. These can be set to $1$ for many sampels and depends on the process being simulated and the generator used to simulate the events (`event[\"mcWeight\"]`).\n", " * In datasets where events are assigned weights, the $\\sum_i w_i$ value represents the cumulative total of these weights across all events in the dataset (`event[\"sum_of_weights\"]`). If no weights are assigned to the events this (i.e. the weights are all set to 1) $\\sum_i w_i$ simply gives the total number of simulated events. Note that the sum of weights is calculated before any selection has been applied and thus is typically much larger than what you get if you sum up the weights of the events in the ntuple. \n", "\n", "When the integrated luminosity is multiplied by the *cross section weight*, $w_\\sigma$, of a certain process, it gives a measure of the total number of expected events during a period of data taking for the given process. \n", "\n", "For `data_D`, the integrated luminosity has a value of $0.105 \\,\\text{fb}^{-1}$.\n", "\n", "For more on cross sections and luminosities,\n", " [see this cheatsheet](https://cds.cern.ch/record/2800578/files/Cross%20Section%20and%20Luminosity%20Physics%20Cheat%20Sheet.pdf)." ] }, { "cell_type": "markdown", "metadata": { "id": "l9Wry4Zyk1Bs" }, "source": [ "Now, with all the weights we've defined, we will calculate a total weight for an event in the $t\\bar{t}$ sample scaled to the integrated luminosity of period D in data15. The total weight is the collective product of all the weights:" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "8Fx8buVGk1Bs", "outputId": "15105514-5c12-47a9-9821-e1cd95edfe76" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "total_weight = np.float64(0.0017336209613845047)\n" ] } ], "source": [ "# luminosity of data_D\n", "lumi = 0.105 #Luminosity of periodD in data 15 is 0.105 fb-1\n", "\n", "# Let's use the first event of our tree\n", "event = tree.arrays()[0]\n", "\n", "# Multiply all the important weights together\n", "total_weight = lumi * 1000 / event[\"sum_of_weights\"]\n", "for variable in weight_variables:\n", " total_weight = total_weight * event[variable]\n", "print(f\"{total_weight = }\")" ] }, { "cell_type": "markdown", "metadata": { "id": "QfR1soaJk1Bs" }, "source": [ "This calculation means that in our final histogram,\n", " this event will be represented with $0.00173$ of a single count in the bin.\n", "We can encapsulate these calculations in a single function `calc_weight`." ] }, { "cell_type": "code", "execution_count": 20, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "rgig0o5nk1Bs", "outputId": "31487cb6-d239-4796-d772-e3e3952cf8a9" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.0017336209613845047\n" ] } ], "source": [ "def calc_weight(weight_variables, events):\n", " total_weight = lumi * 1000 / events[\"sum_of_weights\"]\n", " for variable in weight_variables:\n", " total_weight = total_weight * abs(events[variable])\n", " return total_weight\n", "\n", "# Verify that we get the same answer\n", "print(calc_weight(weight_variables, event))" ] }, { "cell_type": "markdown", "metadata": { "id": "yg1zCm_wk1Bs" }, "source": [ "Now, we can apply the cuts as before to plot the MC data.\n", "The code is the same as before,\n", " but we make sure to add in `weight_variables` to our `tree.iterate()`,\n", " and we store the weights in each event using a new dictionary key." ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "id": "G2-30oDbk1Bs" }, "outputs": [], "source": [ "sample_data = []\n", "\n", "# Perform the cuts for each data entry in the tree\n", "for data in tree.iterate(variables + weight_variables+['sum_of_weights'], library=\"ak\"):\n", " # Cuts\n", " lep_type = data['lep_type']\n", " data = data[~cut_lep_type(lep_type)]\n", " lep_charge = data['lep_charge']\n", " data = data[~cut_lep_charge(lep_charge)]\n", "\n", " # Invariant Mass\n", " data['mass'] = calc_mass(data['lep_pt'], data['lep_eta'], data['lep_phi'], data['lep_e'])\n", "\n", " # Store Monte Carlo weights in the data\n", " data['totalWeight'] = calc_weight(weight_variables, data)\n", " # data['totalWeight'] = calc_weight(data)\n", "\n", " # Append data to the whole sample data list\n", " sample_data.append(data)\n", "\n", "# turns sample_data back into an awkward array\n", "background_ttbar = ak.concatenate(sample_data)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 453 }, "id": "PTFUQWExk1Bt", "outputId": "1e683d1c-7262-4d35-d1d6-721e794bbdb6" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAAsTAAALEwEAmpwYAAA6+ElEQVR4nO2de7gU1ZW33yUoiiAgoFFEwQgRUSBwABMRUUaFGCQoeE2iGZV8KpOLSfzMJDGYy2R0NI4+ON/MIYoaYwRBFIxCIsrEu4AB5KgMqGcCaBTQIMhF0fX9UdXHok9Vd127q89Z7/PUc6qrdu1a1ad7/3rvtfbaoqoYhmEYxl7VNsAwDMPIByYIhmEYBmCCYBiGYbiYIBiGYRiACYJhGIbhYoJgGIZhANA26ISI7K2qH1XSmKh06NBBjz766LLltmzZQqdOnVIpF7aujRs30r1794ralrZ9ebatWvbl2ba07cuzbdWyL8+2Ralv2bJl76tq84Kq6rsB7wC/AUYDElSumlv79u01DJdddllq5cLWNWTIkNTumUW5MPbl2bYs7humXJ5tU03XvjzblkW51vSdADaqT5taasioH7AE+DGwTkRuEZHjy0pPDhk3blxq5cLWFZY0bYtSrhr3tPcu+3JhSfOzbv/X7MuFJUJ9f/c96qcSxRtwKPBt4FngNeCXYa7LegvbQ6gGYRW9WuTZPrMtPnm2L8+2qebbvrRtA5ZqxB6CVzTeBG4H/h+wFbg0rAxlSbdu3aptQiCTJ0+utgklybN9Zlt88mxfnm2DfNtXKdtES+QyEpF9gXHA+cAXgQXAfcCfVPXjilhYgrq6Ol26dGm1zTAMw6gpRGSZqtYVHw/sIYjIvcBfgXOA3wG9VPViVV2QBzEAx6M+efJk5s+fX21TDMMwcs/8+fMLvQ3fUKTAHoKIfB2Yq6pbszMvGdZDMAzDiE7kHoKq3q2qW0XkYBG5XUQWuBUdIyKXZGmsYRiGUXnCOJXvBBYCh7iv/wf4Tkb2GIZhJOLEE09k0KBB1NU1+wGc6jUtkcCZyh66qeosEfkhgKruFpFc+BDyxuLFi5k0aRL3338/o0aNqrY5htEqefLJJytyTUskTA/hAxHpCiiAOzltS6ZW1SAmBoZh1DphBOEqYB7wWRF5Grgb+KdMraoxTAyMLGjTpg2DBg1i4MCBDB48mGeeeSZyHY2NjRx77LEZWJecqVOncuONN/qe+8Y3vsGgQYOats985jMceOCBgXWtX7+emTNnNtsvRZxrgnjttdc47rjj9ji2a9cuevfuTUNDQ+jyIsLee+8dWE+Ue8ShrCCo6ovASTjzEL4J9FfVlancvYVgYmBkwX777cfy5ctZsWIFv/rVr/jhD39Y0furKp988klF71lgxowZLF++nOXLlzN37lzatm3LnXfeGVh+0aJFvPjii832SxHlmsWLF3PxxRcHnu/duzfr16/f4/2qr69n5MiR9O/fP3T5r371q3To0MG3nn79+kW6RxxKZTs9ADhYVde4foNjgP2AwSKyUFXfTsWCFkCxGCxevNjEoQUx5YJVqdY37d7ov9jff/99unTpAsBXvvIV1q1bx86dO/n2t7/dNIv17rvv5sYbb0REGDBgAL/97W/3qOP111/n7LPPpr6+nqFDh/Lzn/+ce+65h+7du9OzZ0+GDBnCxIkTOf300xk+fDjLli3jkUceYc6cOdxxxx0AXHrppXznO9+hsbGRL3/5y6xa5bw3N954I9u2bWPq1Kk0NjYyduxYRowYwTPPPEOPHj146KGH2G+//fjlL3/JXXfdxUEHHdR0z1Js2rSJMWPG8JOf/IQzzzzTt8xTTz3FVVddRefOnbnpppvo2LEjBx54IAsXLuSBBx7gyCOPTOWaUuy1114cfvjhNDY2cuSRR7Jjxw5uuukmFi9eHLn8ypUrfY9HvUccSvUQbgRO8Lz+FVAHjASuS82CBORlYlqxGEyaNKl6xhgthh07djBo0CCOPvpoLr30Un7yk58AcMcdd7Bs2TKWLl3KrbfeyubNm2loaOAXv/gFjz/+OCtWrOCWW27Zo67Vq1dz9tlnc+eddzJ06FCWLFnCnDlzWLFiBY8++ije+Txr1qzhiiuuoKGhgU2bNjFjxgyef/55nnvuOaZPn85f/vKXsravWbOGK6+8koaGBjp37sycOXNYtmwZ9913H8uXL+eRRx5hyZIlJevYvn0748aN45xzzuGb3/xmYLkRI0YwdOhQHnroIXbv3s3w4cN56KGHWL58eWDDHueacvTr149XX30VgNtuu41x48bRq1evyOVL1RP1HsWUm5hWKspoKM4QUYGtqvotABF5KrQFGdKpUyfq6+urbUYTXl9C1vVbD6TlUxgyAnj22Wf5+te/zqpVq7j11luZO3cuAOvWrWPNmjUsWbKESZMmNeX38o63b9y4kfHjx/PAAw9wzDHHAPD0008zfvx49t13X/bdd989smQeccQRHH+8k9j4qaeeYsKECey///4AnHXWWTz55JOBv9YL9O7dm0GDBgEwZMgQGhsb2bRpExMmTKB9+/YAJev4+OOPOe+88zj66KP5+c9/Xva9Wr16NYW1Ubz7Sa8ZPnw4u3btYtu2bbz77rtNz3T99ddz+umn71G2X79+rF69mpEjRzJt2jSef/75kvcPKl+qnqj3KGbcuHGMGzeO6dOn+wYGlRKEtrrnNOavefY7R7KiFZB1Y21i0Lr5whe+wKZNm7j//vt57LHHePbZZ2nfvj2jRo1i586dJa/t1KkThx9+OE899VSTIJSi0PiXom3btnuMZRfb0K5du6b9Nm3asGPHjrJ1erniiiv46KOPmD59etmymzZtolOnTrRt23aP/TSuKTS4ixcv5s477yzpx+jXrx+LFi3illtu4cILL+Tggw8uaUNQ+VL1RL1HVEoNGX0iIp8pvFDVVQAi0gOojqcpp1SisTYxaN28+uqrfPzxx7Rr144uXbrQvn17Xn31VZ577jkATjnlFO6//342b94MwLvvvtt07T777MPcuXO5++67uffeewE44YQTmD9/Pjt37mTbtm08/PDDvvc98cQTefDBB9m+fTsffPABc+fO5cQTT+Tggw/mnXfeYfPmzezatSvwei8jR47kwQcfZMeOHWzdujVwqPe6665j2bJl3H///c0a6dGjR7Nhw4Y9jjU2NnLooYc2209yTRz69evHCy+8wB133MEPfvCDkvcvVT7oeLlzaVBKEP4NmC8iI0Wko7udBDzonjNcKtFY+zmujZZNwYcwaNAgzj33XO666y7GjBnD7t276devH9dcc03T0E7//v350Y9+xEknncTAgQO56qqr9qhr//335+GHH+bmm29m3rx5DB06lDPPPJMBAwYwduxYjjvuON+lFwcPHszFF1/MsGHDGD58OJdeeimf//zn2Xvvvbn22msZNmwYp556aqghmsGDB3PuuecycOBAxo4dy9ChQ5uVaWxsZOrUqWzevJkRI0bs8fyffPIJa9eubRZ+evTRR7Np0yaOPfZYtm/f3rT/zDPPxLomLn379uWll15i8uTJdO7cGSDw/kHlSx0vdy4NyqW/HgP8M9AfZ2JaA/Cvqvpo6pbEIC/J7YqjirKOMir0SDZu3JjZPYxPyUOUURZs27aNDh06sH37dkaOHEl9fT2DBw+utlmBrFq1ijvuuINf//rXmV6TJtW+fxBBye1KCkLeyYsgeMm6sTZfgpEWF1xwAS+//DI7d+7koosuqvg8B6N6BAlCmFxGRkgsysioJQr+BMMoEGoJTaM8FmVkGEatY4KQAhZlZBhGSyC0IIjICBG5SkROy9KgKORlprJFGRmGUQskWULzBVUd5u5fBlwJzAVOA+ar6r9mYnEE8uJUtigjwzBqichLaALeHKyTgVNV9TocQbgwZftqmkrmMsracW0YRuulVJTRXiLSBUc0RFU3AqjqByKyuyLW1RgWZWQYRi1TShA6AcsAAVREDlHVt0Skg3vM8GBRRoZh1DqBQ0aq2ktVj1TV3u7ft9xTnwATwlQuImNEZLWIrBWRa3zOjxSRF0Vkt4hMLDp3kYiscbeLojxUpbEoIyNtfvnLX9K/f38GDBjAoEGDmpKs/fu//zvbt28ve33YchdffDGzZ8/e41iHDh3iGW3UPJHDTlV1u6q+Ua6ciLQBbgPGAscA57uL7Hj5K3AxcG/RtQcCPwWGA8OAn7rDV7nEooxaNt27d2/2fi9evDiV4348++yzPPzww7z44ousXLmSxx57jJ49ewLpC4JheIk1D0FEyqc2dBrytar6uqp+CNwHjPcWUNVGdznO4uyppwN/UtV3VfU94E/AmDi2VoIRA25gdn03plywiikXrGLC6Bmp38MW4akefmLs9yMg6vEg3nrrLbp169aUQrpbt24ceuih3Hrrrbz55pucfPLJnHzyyQBcfvnl1NXV0b9/f376058C+JaLQyFabuLEiRx99NFceOGFFKISlyxZwhe/+EUGDhzIsGHD2Lp1a+z7GDlCVSNvwCEhykwEfuN5/TVgWkDZO4GJntffB37sef0T4PvF1w0ZMkTzwJXnv9S0feWUO3Tfdl0yu9cTTzyh3bp10yeeeCKzexjBBL3/UY+XYuvWrTpw4EDt06ePXn755bp48eKmc0cccYRu3Lix6fXmzZtVVXX37t160kkn6YoVK3zLBXHRRRfp/fffv8ex/fffv8n2Aw44QNetW6cff/yxHn/88frkk0/qrl27tHfv3vrCCy+oquqWLVv0o48+Cv18RvUBlqpPWxyrh6Cf+hOqysaNG6mrq2vakq6elnQYZsPbS1jw9PcYc8JNieoJwhzL1SXrnkGBDh06sGzZMurr6+nevTvnnntu4MIss2bNYvDgwXz+85+noaGBl19+OdIziTSPD/EeGzZsGIcddhh77bUXgwYNorGxkdWrV3PIIYc0pa8+4IADyi5IY1Sf+vr6prYS6OZXJlAQROQAEfmViPxWRC4oOvcfIe6/AejpeX2YeywMoa7t3r07S5cubdoKi43HpfjLHHbMF/YUgx4HN8/znhQTg+qSlhiE/Ty1adOGUaNGcd111zFt2jTmzJnTrMwbb7zBjTfeyKJFi1i5ciVnnHFG2dXTiunatSvvvfde0+t33323aRlOaL7y2e7dFnFeq0yePLmprQQ2+ZUp1UOYgRNeOgc4T0TmiEjh03F8iPsvAfqISG8R2Qc4D5gX0vaFwGki0sV1Jp/mHqsIUb/MWYsBWJRRNUlTDML4flavXs2aNWuaXi9fvpwjjjgCgI4dOzaN17///vvsv//+dOrUibfffptHH/10mRJvuVKMGjWKmTNn8uGHHwJw5513lvU7fO5zn+Ott95iyZIlAGzdutWEoqXgN47kDDGxvOj1j4Cnga7Ai0HXFV3zJeB/gNeAH7nHfgac6e4PBdYDHwCbgQbPtf8IrHW3b/jVn4UPodxYsB/7tuuiXznljj18CVnYVeq1kR1p+Ayi+BKWLl2qX/jCF7Rfv3563HHH6YQJE5r8Abfeeqv27dtXR40apaqOD6BPnz56yimn6IQJE3TGjBm+5S655BJdsmSJ7/2mTp2qxx57rA4cOFDPOussfeedd5psPuOMM5rKXXnllU31v/DCCzp8+HAdMGCADh8+XLdu3aobNmzQsWPHln0+o/oQ4EMolcvoFaC/qn7iOXYx8AOgg6oekZ4sxSPtXEZxx4InjJ6xR89gw9tLmLvoG6nZBXuu2lXokezY+W6JKwzDMPyJk8toPnCK94Cq3gl8D/gwVetyQlzHYLEYLHj6e5nZGMdxHcUXYhhG66XUTOWrVfUxn+MLVLVPtmZVh6RRIllHGcX1VYT1hUR1pBuG0bKwBXI8pCUGWTiWk9Qf1rFpjmvDaN2YIPiQdZRRnF/iaYhNuWysJgaG0bopNQ/h0EoaEocsVkwr51j2y1ETVQzi/BIvrn/D20tCXxv2vpYryTBaNuVWTCsVMvoI8Bzwr8AooG1Q2Wpt5cJOS4WKliofNZSwOOT0K6fcEan+MPilxwgbihr1vlHfN8MwagsCwk7LzSPYFyep3C3AUuABnNXTDi91XaW2UoIQp/GNG1ceNpdRkhxExfUXi06pRjyOGNg8B8NoucQShGaFoTdwBc6M4xeiXJvFFiQIcScHxZ1kVK6xLlV/WIrr905+i/NcQeVMDAyj5ZOKIOxxIewT99q0Nj9BSCvrZJR6SjXW5eoPS3H9hXukORxkYmAYrYPUBSEPm58gVFoMVIMb66j1l8IvPUZaz+U9bxhGy6fVCEKlxeCJJ56IlMso7rCMn+M6S9+ADR8ZRsslLR9CF2BAlGuy3MJGGSUdQy93vKVGGYWt3zCM2iK2IACLgQOAA4E3gOeBX5e7rhJbnCijqGPoYUSiNUQZRfV5GIaRX5IIwl/cv5cC17n7K8tdV4ntqKOO0ssuu0znzZu3x8NG/aXvPR+nniRRRmEb62pGGdnwkWG0DObNm6eXXXaZAms0piC8BBwC/BEYqjkShLZt28ZuxMsRpZ4kUUZhG+tqRRmZGBhGyyNJD2EisBL4D/f1kcCcctdVYuvbt+8eDxlVDNI6nlaUUanGuhpRRiYGhtEySSIIJ4Q5Vo3N60NIo3H3+6UfxpHrLRc3yqhc42tRRoZhpEUSQXgxzLFqbAVByPKXfhhHblhBSDKcFVacsh4uMwyj9oksCMAXcFZHWwdc5dmmAiuCrqvkNmTIkFSHg8qJQZAj108Q0m6sw4pTlmJgM5wNo2UQRxBOAn4KvOX+LWxXAX2Crqvk1rdv31R9A+XEIMiRW66HkEZjHSfKKOov/bh2Wk/CMGqLJENGR5QrU60tjSijoF/6QcNHfvWUEoS0ejBB94j7XH6k2cMwDCO/BAlCmBXT2olIvYj8UUQeL2whrsucI488MtSyl1GPB62AtuHtJZEWt0nLnlJrIMd9Lj+SrildbkU2wzDyTdsQZe4H/hP4DfBxtuZE45NPPmHy5MmMGzeOjh07Zi4GC57+Ho8ueMCnse7WzLY0xWDSpEmce+oTze6RxnN5mV3fjdn1q5rKP7Xy6sTiZxhGfpg/f35hhcloK6YVNmBZuTLV2vIcZRRlvkLUYZ+ojuuw8yQqlSvJMIzqQoIho/kicoWIHCIiBxa21CQrIWn+Ei/VMygcn3LBKqZcsIoJo2cwdsxZjBhwg69dYeuPaqeXJM+1ePFiunfv3mw4yvu8fvX7YT0Dw2gZhBkyusj9+wPPMcWZsVxVtm7d2qwxmnLBqj0atdn13Si0U+Ua33JiUO64l7D1T5rUfFgm7PDLlAs+Hd4p1YiHvW/Qc5XzDfjZb+JgGLVH2R6Cqvb22aouBgCvv/56rEYtqPFNSwy8lKs/qSM3zPOGvW+U981L2J5EGgT1bAzDSE5ZQRCR9iLyYxGpd1/3EZEvZ29aefyijPwaNe/wzuz6bk2/rOM6luOIQVD9SaKMoopf1J5E3qKMzHFtGNkSxocwA/gQ+KL7egPwizCVi8gYEVktImtF5Bqf8+1EZKZ7/nkR6eUe31tE7hKRl0TkFRH5oV/9HTt2bNqPOuyTlhhseHuJ77NHrb9cD8aPNJ6rVE8ijRDbtDAxMIzsCSMIn1XVG4CPAFR1OyDlLhKRNsBtwFjgGOB8ETmmqNglwHuqehRwM3C9e3wS0E5VjwOGAN8siIUfcYZ90hKDBU9/z9emNMUg6Bd30ucqJx5J50OkiYmBYWRPGEH4UET2w3EkIyKfBXaFuG4YsFZVX1fVD4H7gPFFZcYDd7n7s4HRIiLuvfYXkbbAfjg9lPf9bhL3l37cKCO/435UIsooC99AsUO+OKrKj0o01mHFyTCM+IQRhKnAAqCniPwOWARcHeK6HjiJ8Qqsd4/5llHV3cAWoCuOOHyAk0fpr8CNqvpu8Q3Wr1/PqaeeyoEHHsj3v/996uvrgXCNeH6ijOKJgZd4UUbxekh+VKKxrqTj2jBaIvX19dTV1VFXVwd+s2kJEXaqqn8UkWXA8ThDRd9W1U2pWtqcYTizog8FugBPishjqvq6t9DmzZv505/+tEdjMWHmjEiNuLfxnV3frWz5JI7loPqTzBCO+1x+M67TFL+NGzeWtT0OWTuuDaOlMnnyZCZPngyAiPi24WGijOYDpwGLVfXhCGKwAejpeX2Ye8y3jDs81AnYDFwALFDVj1T1HeBpoK74BmGjjLJ2LAeR5nCWH1lHGSV53iwwx7JhZEuYIaMbgROBl0VktohMFJF9Q1y3BOgjIr1FZB/gPGBeUZl5fDrxbSLwuDut+q/AKQAisj9O7+TV4hu8t6ltszHuWo0yijNDOOsoozTELy1MDAwje8JMTPtvVb0CZ2byfwHnAO+EuG43MAVYCLwCzFLVBhH5mYic6Ra7HegqImtx1lkohKbeBnQQkQYcYZmhqiuD7hWnES+el1A8w7mSUUZxZwhnHWWUVPzSxMTAMLJHnB/kZQo5UUbjgHOBwcDDqvpPGdtWloO69tcTBn0/dnRQFscLTJy8qdkv/bA9mGn3HtuskS1Mpgsiqp0TJ28qm/bDr54dO5v59unevXvmjfXionQYxa8NwwiPiCxT1WbD8GUFQURm4Th5FwAzgf9W1U8ysTIiXQ7opTs/fD+1xv2Ss/4c6r5hhlOm3Xts036hcR8x4IZEjXVUe0od93NcTxhd3iE/d9E3mt1/wugZzcr7lUuLrB3XhtHSCRKEMMntbgfOV9VcrYUAsGXbesafPD21X/phCDu27pd4LvykueRRRnHWdYjrSPcfRstGECzKyDCyI9CHICJXA6jqQuCsonP/krFdoejU4bBMh338iFI+rj3Fvo0Jo2ekWn9WUUZRxTUq5lg2jGwp5VQ+z7NfnEtoTAa2RGbvvfdv2q+EGIC/I9ePtH0VfkR9rh4HD+WSs/6cSZRR3PczLCYGhpE9pQRBAvb9XleFXR9u5YkXprK0YXpFxACChkeak7bj2o9qzJMot6hOFmIAFmVkGGkwf/78wuQ03yU0SwmCBuz7va4K7fbpSN8jzmDF6t+m0vhGIa3GOoloZTFPolBvcU8iaq6ktMkqPcZiW1/BaEWMGzeukOJni9/5UoIwUETeF5GtwAB3v/D6uAxsjcxHH32Q+bCMH3Eb6+JGFqIP+0Sxp9Rz+c3DCFN/WN+DH0ka3yxyGdkwlGHsSaAgqGobVT1AVTuqalt3v/B670oaGcSWbeszH5YpJm1Ha9ZiEOW+Ueov53vwI43GN60oIxMDw2hOqIlpeaXLAb30wi8/3PR6w9tLfBvRLHwJYUizcY9iX9z7Fs/DiDpprtx8hSSTy/wmzXnnekSlEpPpDCOvBM1DCJPLKLf4RRn5UetikPZwVlD5YtJId+ElybBP2v8vW1/BMJpT04JQIOthmThUYtgnq/sWSJorKYg4wz5B4pQnn4Rh1Do1LwhRHbyF10mjjMqR1XyFrO/rJYsoo6AeRrnGPaw4RcVmPhvGp9S0IARFGYUhybBMmvXnaZ5EuXovOevPiaOMsl4uNArmWDaMPalpQdiydR0Hdz2OD3d/ELuOrNMtVGs4K60ooyDKRRkVN85JxCDoufzq8btvOftNDIzWQpKJabmnU8eefPmk2+jdYxSQbHJZtRzLWU+ai5aQLhxhHMthGv0kYhDG0V3KN2BiYLRGkkxMyz1ho4z8yIMYFFNqhnBW943aQ4oaZRRVDIJ+0Sep3w+LMjKM5tS0IHgpNKZhyVOUURDVjjLyI2sxSGO50DA9D4syMozmtBhBiIpfGok0qcSwT5L7xq0/bJRRXDFIulxo1GEoizIyjE9ptYKQNZUc9gm6r9/rAnHrDxtlFFcMkkYZxRUD8yUYhglC5mQ97BOXrKKMwGm80xQDL3Ec1+XsNzEwDAcThIzxcxSXI+tJcwW7vPWnEWXkJYkYBE1Ssygjw8gWE4QcEnbYJw3SijKKWj7qcYsyMozsqelsp506Hq6HHTyMXj1GNc1FMMITx1dx+wMjQ5X3y4paIM7xsWPOanbfiZM3xR6G8q7/UHgfdux8N9R7YBi1yvz585k/fz7Tp09fq6p9is/XtCAc1LW/nnP6zGqbUZPEdVwXp+IOSs0dtbEud3zEgBua3feplVfH9kmETdltGC2RFpn+2ohP1rmS/BrlCaNnNFupDcKJRBZRRllPTgzyhRhGXjFBaKUknYdRzveQ1XyFtKKMKiEG5rg2ag0TBCMycRfDSTpfIc0oo6xDe9NKzW09DKOSmCAYkUiyGE7SYZ/iEN4pF6xqNgzlV48fWYf2RpkP4dfoWw/DqAaZCoKIjBGR1SKyVkSu8TnfTkRmuuefF5FennMDRORZEWkQkZdEZN8sbTXCkcZiOF5RKfgSvI17VgnvvGS9HkaS+RCVEAPreRh+ZCYIItIGuA0YCxwDnC8ixxQVuwR4T1WPAm4GrnevbQvcA/wfVe0PjAI+yspWIzxJF8MJs0JckmGfJI7lLAjqqRR6BtUQA7B5GIY/WfYQhgFrVfV1Vf0QuA8YX1RmPHCXuz8bGC0iApwGrFTVFQCqullVP87QViMkYX9ZF0cTTblgVeTlQuMM+/jdN4hqOZbjHE+bsKJrtC6yFIQewDrP6/XuMd8yqrobZ9GGrkBfQEVkoYi8KCJX+91gx873mLXw3KatYa1lrKwUaa0El/awT9bZXsOSphhk2ViXE12j5VBfX09dXR11dXUA3fzKtK2sSaFpC4wAhgLbgUXuRIpF3kL77dsFm5hWedJcCS7tYR//HkzzCWdp3LcUaYpBVo21Oa5bF5MnTy4sn4mIbPIrk2UPYQPQ0/P6MPeYbxnXb9AJ2IzTm/izqm5S1e3AI8DgDG01QpJkRbawIpGmGETpwaTZOKYpBlk01iYGhh9ZCsISoI+I9BaRfYDzgHlFZeYBF7n7E4HH1cmlsRA4TkTau0JxEvByhrYaIclaDIJCUdOo30uY+yYhz2IAlu3V8CczQXB9AlNwGvdXgFmq2iAiPxORM91itwNdRWQtcBVwjXvte8CvcURlOfCiqv4hK1uN8CQdo886yihJDyYLR24exQAsysjwJ1Mfgqo+gjPc4z12rWd/J+D7bVfVe3BCT40c4df4+q1lHVcMehw81DfxnB9pioGfL2Hjxo1l6yhFmo7ltMXB775Jn9eofWymshGLrKOM0nRcl7M/SqhrWCzKyKhFTBCMyGQdZZS249qPqDmRos7stSgjoxYxQTAiUYkoo7Tr9yNpor1yZBVlFEec/DAxMPyo6QVybMW0fJOWGAQtwhO2nmn3HtvsWr8ZzBMnb0rNweu3ItujCx7ITfSRX9oMo+VTbsW0mu4htNunIycPm2pikENK/dIvzlpaqnzUYZ8kM5CzaHzjrN9QiegjizJqnYwbN476+npwskI0o6YFwcgv1Rr2SSMdRVo5heKu31CJUNRKOq6N2iGvqSuMGiftYR8/wkcZNR8yCqJcoxw2NDOsGPiF2M6u78bs+j2P+w03pSEOFmVkeLEegpEJlR72KTUMVfzLPuiXftRon4KDN+xzBdWf5mS9sJRaXMhovZggGJlSrWGfOI1pmiGhYVNzp50SPCxZZ3s1ahMTBCMzog77FF4nXd4ybmOa5ph+FvMw0vQlmBgYftR02OlBXfurpb9u2ZRKjxGWafceGzh2X+q+T628OraD1883kEQMJk7elKpjecLoGc3uO3dR8zThRsvEXU6grvi49RCM3JLW8pZRf1mXyrqadmruLFKChyHrNaWN2sQEwcglaY5xx23Eww4fRbE/zZTgacxYznpNaaO2qGlB2PXhVp54YSpvbFhcbVOMFEnb4Zm0EY+TeC5NMUgqTkGYY7n1MX/+/MKqaZ38zte0INhM5ZZJ2o1UkkY8zSijLB3LheNhQ2xNDFonNlPZqDmSRhkFEbURr5Uoo7R6MIZhgmDkDr9Q1KTEacTTSCNRbTGI0oMxDEtdYbR4kjTixWkk/EJRw9STtMcwdkxze4JCUeP2YMDCTls71kMwWjxp/qKvVpRR0vUbgqKSLMrI8GI9BKPF49+IN8/dEzfaxy/hXSG3UtT6oxxPur6COZaNYqyHYLR4sggVjZNTqFqhqCYGRlish2C0ePzSSPiRdU6haoWi+h03MTD8qOkegk1MM8KS1voKxSmjy+VEKld/ccruQrmwK8rFjT6yKKPWiU1MM1o9UYdHkizaE9WesKQdihomtDet9BhGfrCJaUarJ0sxiPPL2q8HUI40ooyiDHOlmWo7KiZC1cMEwWjxhG3E447pV4KwIbBZikGlGulq3dcwQTBaAWEb8Shj99WK308ryiiIOGkw0iZJjiYjGZkKgoiMEZHVIrJWRK7xOd9ORGa6558XkV5F5w8XkW0i8v0s7TRaB1Eb8ax9CVEpvm9hGU6voztJVtS4aTCyIg/i1NrITBBEpA1wGzAWOAY4X0SOKSp2CfCeqh4F3AxcX3T+18CjWdlotB7iNOJ5FoNyx+M0pmn0MNIib+LUWsiyhzAMWKuqr6vqh8B9wPiiMuOBu9z92cBoEREAEfkK8AbQkKGNRisgbiOeFzGA8OJU3GOYXf/pjOxyjWktiUGlHd2thSwFoQewzvN6vXvMt4yq7sYJheoqIh2A/wtcl6F9RishbiOeVZRRHLKavOYlL42viUH1yKtTeSpws6puK1Vox873mLXw3KatYa11I43mJG3E8xpllNXM6jxGGZkYJKe+vp66ujrq6urAL5kXIKqayc1F5AvAVFU93X39QwBV/ZWnzEK3zLMi0hb4G9Ad+DPQ0y3WGfgEuFZVp3nvcVDX/nrO6TMzsd9omXhnCKdRz5gTbmLuouZpo6PMYE6DNBzg0+49tuxwjV8ivywpZY+JQ3xEZJmq1hUfz7KHsAToIyK9RWQf4DxgXlGZecBF7v5E4HF1OFFVe6lqL+DfgX8pFgPDiEpaoaJ5SwyX1mS6vDly4zjGbVJbMjITBNcnMAVYCLwCzFLVBhH5mYic6Ra7HcdnsBa4CmgWmmoYaZBF2olaFYOgYa6sx+6jNNZxxckmtSUjUx+Cqj6iqn1V9bOq+kv32LWqOs/d36mqk1T1KFUdpqqv+9QxVVVvzNJOo2WTZiPeEsQgqIeU9dh9mvMhojjGWxqFHFNZkFensmGkRpzcQUGkVU8apJE620txY1qc1TWOX8SbIK+S8yFa6nyFrJ/L1kMwjBolTTEA/3Ujmvsejg1tX5rDPrWSmC9LKvFc1kMwjBolTTEIW09Ykgz7eHskhUl2SYebStlZC47oSomc9RAMo8aJ40vwE4c0E/mlMezjve/s+m7Mrt+zB7Nj57uh7htELfUkKmVnTfcQbMU0o7WTZpRRmrmb0hSDpI5xP2pJDCC96ClbMc0wWjBppuxOM3dTkqyraTrGoziu41L8HFkMP6UVPWUrphlGKyJJY5pFIr84k8vS6KkEOa4LPoM05ytUMtQ16ygjEwTDaEFE6TEEkVYiv7hRRmmKQdieShqNeNaNtUUZGYaRmFqLMso6/UYWaTmybqwr5fMwQTCMFk6aa0qHJe0ooyzSb6TVyFaisa6UA9zCTg2jhVPcMwjqKWQdZTR2zFnNQkin3es/0S2uOBVmVntDVJ9aeXWmaTkq0VhXKtur9RAMwwCyjzJKs6cSxTGelxxNSfCbrJdFdJMJgmEYQLpRRsWN14gBN2QyHyJoprG3fJIQ2DBUMsqo3HMlFYmaFgSbmGYY6eNdRKhYDKKkech6PkRQ4x61fFAjHjWtRdZRRmk8V7mJaTXtQyhMTDMMIx3izAT2G3bJej7EhreXMGlSc99A2EbTL5FfMXlKgxFXDIrFady4cYwbN47p06fbxDTDMEpTLlV4JSeXBdlTyjeQRk+igF/9fj2GSkQZxRWDMPZ7MUEwjBSIOrxQ62Q9uawUBZEI67iOIwYQPv1GJaKMwj5XUjtresjIMPJApRqFvJDu5LLw6yuEsWd2fbdm9310wQOhehJh6vdrZP1CXecu+kbi5/Li9376PVdS0bIegmEkxMQgmLRmPke1pzh1dnHUU5ozn7N4riDiRE9F+X+ZIBhGQlrTwu6VSGEdlqi+gXLHw/oGwtw3C5I4lsN+Pk0QDCMhlYxDrzaVSGEdljTFIOkv7qjzJKKSphiU+nyaIBhGSmQdh54HKh1lVIqwUUZRRSJNMYgznOjXU0lTDEp9PmtaEGximpEXWrNj2Y+sxQD8xckvZDZLMYjaWMd9Lj/7CzPAvT6ScvZ/97vf5d5774WAiWmiqpGMzRMHde2v55w+s9pmGK2ciZM3NfsSFiY+tUSKG1O/BHXe5/fOfC6m3LVx7PE7H0YM/P6PEC7KqLh+v4R6UQmKnkrjuURkmarWFddZ0z0Ew8gDralnkGRyWbXsCT/zubmoe3MxeRvlclFGSQMN0gzttSgjw6ggrSnKKO1hnySkuRKcXyOeJIdSkkCDtEN7o3w+bWKaYSSkeFJS0PBISyCryWVxSEMMwDvz+dPnCPINTBg9Y4/5DcX/91LiEZa0Q3v97AzCegiGkRJZx6HngUpOwipH0vUVgijlyE3DQV2OrEJ7w3w+MxUEERkjIqtFZK2IXONzvp2IzHTPPy8ivdzjp4rIMhF5yf17SpZ2GkZSsoiiyTNBjUsh7r4SJFlfIYhyvoGsxQCyCe0N+/nMTBBEpA1wGzAWOAY4X0SOKSp2CfCeqh4F3Axc7x7fBIxT1eOAi4DfZmWnYSSlNYtB0uGRtO3xI2xPIsiBXK6R9Qt1TXut5rDPlXT+R5Y9hGHAWlV9XVU/BO4DxheVGQ/c5e7PBkaLiKjqX1T1Tfd4A7CfiLTL0FbDiI2JQbqTsNKwx4ufz8CPNIdfiucHTBg9I9JzpRllFOXzmaUg9ADWeV6vd4/5llHV3cAWoGtRmbOBF1V1V0Z2GkYiwv4CbQlUa3jEj7R7Zmn+4k7iaynuqZSblxHHlxBErqOMRKQ/zjDSaX7nd+x8j1kLz2163f+zE+l/VMvNI2Pkk6BhhJZI2JXLorLnJKxw16TdM8viF3ecQIM0Q2m9xxvW3s/yV+9iy7b1AM1nupGtIGwAenpeH+Ye8yuzXkTa4kyn3gwgIocBc4Gvq+prfjfYb98u2ExlIy+01igjv7z8Ae2NL8U9jEKIZDn8GsE0xCELMYhiV5Tn8pv0F3Tfzh17sfPD9xl/8nQefPwfN/nVl+WQ0RKgj4j0FpF9gPOAeUVl5uE4jQEmAo+rqopIZ+APwDWq+nSGNhpGKrQ2xzKUXrksLEmGm8L6BqKQ1vBLkuynSUN7C/+XXEUZuT6BKcBC4BVglqo2iMjPRORMt9jtQFcRWQtcBRRCU6cARwHXishydzsoK1sNIwmtUQz8yHoSVtaEmbyWJNQ1yQzkJET5fFpyO8MwEjHt3mObNXZhE9RNnLypWQ/DL5FbJbj9gZGpibpf4rmwCe+mXLAq8x8Zt/3+OEtuZxhG+kSNivGSp8WFsnZQJ52BXAlMEAzDSEQajVceFhfKInQ4TsK7ag4/miAYhpGIJPMwSqWYriZZOqiDxK/gcK7mvBYTBMMwEpEkKqYlO+T9HNFhVmSrZgLBmhYEW0LTMPJD1pOwap2oy3NC+osLvbFhMU+8MBUCltCsaUFot09HTh42ld49RlXbFMNo1aQ5CaulElUMsqB3j1GcPGwqOGmCmpHr1BWGYeSfJMM+rSntR8FH4l2sJo21l9PEBMEwjESkMexjaT8coqb9SJuaHjIyDKP6JB32acmOZT/SWns5C0wQDMNIRJKcQq1ZDPKwuFAxNmRkGEbVaO1iUJymotrzMEwQDMOoGi3VgexHGim1s8aGjAzDMCpA3sUAalwQbGKaYRi1QtLlLdOg3MQ0S39tGIaRAyo5D8PSXxuGYeSUvMzDMEEwDMOoInnyJZggGIZhVIk8iQFY2KlhGEbVSDubaVKsh2AYhmEAJgiGYRiGiwmCYRiGAdS4D6EwMa1Xj1G2SI5hGEYZ3tiwmEZnIq9NTDMMwzBsYpphGIZRBhMEwzAMAzBBMAzDMFxMEAzDMAwgY0EQkTEislpE1orINT7n24nITPf88yLSy3Puh+7x1SJyul/9O3a+l6H1yWhYW92l8MqRZ/vMtvjk2b482wb5tq9StmUmCCLSBrgNGAscA5wvIscUFbsEeE9VjwJuBq53rz0GOA/oD4wB/sOtbw927gonCGHXSwhTLmxdDa/NTu2eWZQLY1+ebcvivmHK5dk2SNe+PNuWRbnW9J0gIOw0yx7CMGCtqr6uqh8C9wHji8qMB+5y92cDo0VE3OP3qeouVX0DWOvWF4vGkG9SmHJh6wpLmrZFKVeNe9p7l325sKT5Wbf/a/blwhKhvs5+BzObhyAiE4Exqnqp+/prwHBVneIps8ots959/RowHJgKPKeq97jHbwceVdXZRffYCXzsObQR2ORjTidgSwizw5QLW1e3AFvi1pd2uTD25dm2LO4bplyebYN07cuzbVmUa+nfiW5Ad3d/L1Xdr7hATc9UVtV9q22DYRhGSyHLIaMNQE/P68PcY75lRKQtjrptDnmtYRiGkSJZCsISoI+I9BaRfXCcxPOKyswDLnL3JwKPqzOGNQ84z41C6g30AV7I0FbDMIxWT2aCoKq7gSnAQuAVYJaqNojIz0TkTLfY7UBXEVkLXAVc417bAMwCXgYWAFcC3xKRBhFZJSK/F5F9XbF53g1PnekKT0UQkTtE5B3XD1I4dqCI/ElE1rh/u7jHRURude1cKSKDq2Dbv4nIq+7954pIZ8+5siG+WdvnOfc9EVER6ea+rvp75x7/J/f9axCRGzzHK/beBfxfB4nIcyKyXESWisgw93il37eeIvKEiLzsvkffdo/n5TsRZF/VvxdBtnnOV+47oaq534AewBvAfu7rWcDF7t/z3GP/CVxeQZtGAoOBVZ5jNwDXuPvXANe7+18CHgUEOB54vgq2nQa0dfev99h2DLACaAf0Bl4D2lTaPvd4T5wfEP8LdMvRe3cy8BjQzn19UDXeuwDb/giM9bxXi6v0vh0CDHb3OwL/474/eflOBNlX9e9FkG3u64p+J2pppnJbYD9xfA3tgbeAU3DCVcEJX/1KpYxR1T8D7xYd9obReu0ZD9ytDs8BnUXkkErapqp/VKfXBvAcjl+mYFtqIb5x7XO5Gbga8Ia+Vf29Ay4H/lVVd7ll3vHYVrH3LsA2BQ5w9zsBb3psq+T79paqvujub8UZFehBfr4Tvvbl4XtR4r2DCn8nakIQVHUDcCPwVxwh2AIsA/7u+Weu59M3sVocrKpvuft/Aw5293sA6zzlqm3rP+L8woCc2CYi44ENqrqi6FQe7OsLnOgOT/63iBRWQ8+Dbd8B/k1E1uF8R37oHq+abeJkHPg88Dw5/E4U2eel6t8Lr23V+E7UhCC4447jcbpuhwL748xgzi3q9O1yt9iEiPwI2A38rtq2FBCR9sA/A9dW25YA2gIH4nTPfwDMEhGprklNXA58V1V7At/F8ctVDRHpAMwBvqOq73vP5eE7EWRfHr4XXttcWyr+nagJQQD+AXhDVTeq6kfAA8AJOF2lwlyKPISmvl3ourl/C0MLuQijFZGLgS8DF7pfTsiHbZ/FEfsVItLo2vCiiHwmJ/atBx5wu+gvAJ/gTPLJg20X4XwfAO7n02GNitsmInvjNGi/U9WCTbn5TgTYl4vvhY9tVflO1Iog/BU4XkTau7/MRuNEID2BE64KzhfjoSrZV8AbRuu1Zx7wdTc64Hhgi6cbXRFEZAzOWOSZqrrdc6rqIb6q+pKqHqSqvVS1F04DPFhV/0YO3jvgQRzHMiLSF9gHZ9Zo1d87HJ/BSe7+KcAad7+i75v7vbwdeEVVf+05lYvvRJB9efhe+NlWte9EWt7prDfgOuBVYBXwWxzv/5E4/6S1OL+O2lXQnt/j+DM+cv9ZlwBdgUU4X8rHgAPdsoKT6O814CWgrgq2rcUZd1zubv/pKf8j17bVuBErlbav6Hwjn0ZU5OG92we4x/3svQicUo33LsC2ETj+tBU4Y+JDqvS+jcAZDlrp+Yx9KUffiSD7qv69CLKtGt+Jml5T2TAMw0iPWhkyMgzDMDLGBMEwDMMATBAMwzAMFxMEwzAMAzBBMAzDMFxMEAzDMAzABMEwDMNwMUEwahIRWSwiP662HVERkUdF5OrWbkNrQ0QaRWSniDRU6f5/FJEdIrK7VDkTBCNXiMheIvKMuyDIYeWvSO2+FREYVR2rqjeUL1kbNtSqMFeJS1W1v/eAiAwRkTniLHy0zRWOOSJySpgKReQhEbk74NwTIjINQFVPA8aWq88Ewcgb3wW2ly1lRMZNoGbkBBE5FXgaJwVFHc7iOMcB9wITQlbzX8BE8az05tbdByfH1X9FMirL/CG22RZlw1l34DVgEE5ul8NKlF0M/Njzuj3OegBv4CwiswA4ynO+ESeV8FPANmApMNQ9Nw34GNjlnlvtHu8K3I2Tx/9vOAu8HFhU5z/j5OrZhpPr6ItlnrHJ7nLX4ywdu7zo+t6urb2Ab+Pk99qKkwDyV3hW9fI88xNu/ef52BCmDl8bg963gOduBH7sseUlYABwPk4+oS3Ab/h09bJydn3L/V9vxcn0+S9hziW1K8FnuxH4atGxtcBvQlwb+NnG+VH/v8A/FV3zb8CzRcdGAbtL3qvajYBttqk2fbCfwvll1IvogvA74GGcBVj24dNkiHu75xtxMoMOcc9fA2wEDvCrzz22AJgPdHG3PwB/8JxvdL/U/YE2OKtbrSnznE33KXe9e8+dwCDPseuARe7+2TgCITiLqrwNfLPIvnXuOeHTJWi9NoSpo5SNzd63gOduxElw1w/YGydZ4GtAPc76JofjpMa+sJxdOD8ctgP93dedgePLnUvDrgSf70Y8guDaqcDoENeW+2z/BFjpKb+Pa/PFRfWMwgTBtlrYcIaKZrv7vYggCDhrEyhwuOf8Xji/7ka4rxuBn3vOC84vzwuK63NfH+rW2cdz7HPusUM8df7Ac76/e75TSLvLXg/MBG7x2NwY1Djh/Iqc5XndCFxbyoaQdQTaWKquonqL6/mSW093z7FZwM3l7MLJcrwDOAfoUFQu8FyaduH0Ht/DbeRxPrOPlbmPVxBOcO9ztOfYmcDfcT63OyN8tg/FyYA73H19rmvbfkU2jKKMIJgPwag6InIU8D1gSsD5C12H2zYR2eZTpLf7d6WI/F1E/o7Ttd6bPRcSaSzsqPMN+SufrqFbTOG6NzzHXis6B0466gIfuH87BtTpR7nrZwAXuOP/p+D84n0AQETOF5ElIrJZRLbgDDF1L6q/sdTNQ9aR9Bn96tkOfKyqG4uOdSxnl6q+DlwIXAa8KSJPichp5c6lYZeHH+P0aOOyyf3b9PlT1Xmq2hk4Aye9P4T4bKvqmzg9iMlu2cnAPaq6I6pRJghGHhiB82VfJSKbcNYcAOdLcIWq/k5VOxQ2n+v/1/3bR1U7e7b2qvp7T7lehR13UZLDcdYVAGcVNC/riq/B+fXpPVcJ/oQzRj8OuBhn4fcdItITZ3jjFzg9lk44OfKLl/Ysfq4mItRRisD64xLGLlV9QFVPxfkFPQt4SJylWEueS8m+o3B6CMsSVPM/wOvAeWXKhf1s1wPnisjncRZziuZMdjFBMPLALJwlAwe525fc46fhOHVLoqrv4ERm/IeI9AAQkc4iMkGcdWoL/KOIDHZ/bf8Ax1n3B/fc34CjPHW+CfwRuMmtqwtwE/CoVnDFNlX9GOc9+BZwFnCHe6oDzvd3I/CRu3LW1yJWn0Yde7xvKVHSLhH5nIiMcRv5j3CGTxT4pNS5FO37Gc44fmzcHuqVwNdE5HoR6emugNYeGO4pF/azvRCn1zEHx5m8Ko5dJghG1VHV7aq6vrDhNDIAf1NVvyEiPy7DWdlqsYhsxYkWmcSei7rXA7fijK+eC5yhqlvcczcDdW63vDB56Ks4kSqrcZx4fwe+HucZEzIDJ4TwDXXWdEZVXwF+irMk5d9xnOS/D6rAjzTqwP99S0QIu/bBiZ56yz3/LeBsVd1Z5lxiROSLwGZVfa1s4TKo6gKc3nFfnF7xNqABx7/gnYdQ9rOtqp8A03GGmOrj2mQrphmtAnEWKv+xqt5TbVuM2kVEvoUTAbUDp2f0AfB/cAToN6r6DwHXrQYOARpVdUCFzPXe/1EcodkrYNgVgLaVM8kwDKO2UdVbcXqZiMhUYK2qPisivcpc97nsrSt5/7KzlMEEwTAMIxaqOtWz3wj49g5qCRsyMgzDMABzKhuGYRguJgiGYRgGYIJgGIZhuJggGIZhGIAJgmEYhuFigmAYhmEAJgiGYRiGiwmCYRiGAcD/B3tH2eZWU3NhAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "step_size = 2.5 * GeV\n", "\n", "mc_x = ak.to_numpy(background_ttbar[\"mass\"]) # define list to hold the Monte Carlo histogram entries\n", "mc_weights = ak.to_numpy(background_ttbar[\"totalWeight\"]) # define list to hold the Monte Carlo weights\n", "mc_colors = samples[r'Background $Z,t\\bar{t},t\\bar{t}+V,VVV$']['color'] # define list to hold the colors of the Monte Carlo bars\n", "mc_labels = r'Background $Z,t\\bar{t},t\\bar{t}+V,VVV$' # define list to hold the legend labels of the Monte Carlo bars\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] # 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 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", "# add minor ticks on y-axis for main axes\n", "main_axes.yaxis.set_minor_locator( AutoMinorLocator() )\n", "\n", "# draw the legend\n", "main_axes.legend( frameon=False ); # no box around the legend" ] }, { "cell_type": "markdown", "metadata": { "id": "V0jfjVO-k1Bt" }, "source": [ "## Final Analysis" ] }, { "cell_type": "markdown", "metadata": { "id": "KRv7JhYQk1Bt" }, "source": [ "Now that we understand all the steps of our analysis,\n", " all that's left is to import the entire ATLAS data and implement it.\n", "The `samples` dictionary will be useful for this.\n", "\n", "We will loop over all values in the `samples` dictionary.\n", "Depending on whether it is a data sample or MC sample,\n", " `fileString` will change,\n", " which opens the correct file on the open data folder.\n", "As before,\n", " the cuts,\n", " mass calculations and MC weight calculations will be performed for each sample value,\n", " and then stored in the array.\n", "The data will all be concatenated into `all_data` for plotting." ] }, { "cell_type": "code", "execution_count": null, "metadata": { "id": "chBwdGa4k1Bt" }, "outputs": [], "source": [ "# Set luminosity to 36.6 fb-1, data size of the full release\n", "lumi = 36.6\n", "\n", "# Controls the fraction of all events analysed\n", "fraction = 1.0 # reduce this is if you want quicker runtime (implemented in the loop over the tree)\n", "\n", "# Define empty dictionary to hold awkward arrays\n", "all_data = {}\n", "\n", "# Loop over samples\n", "for s in samples:\n", "\n", " # Print which sample is being processed\n", " print('Processing '+s+' samples')\n", "\n", " # Define empty list to hold data\n", " frames = []\n", "\n", " # Loop over each file\n", " for val in samples[s]['list']:\n", " if s == 'data':\n", " prefix = \"Data/\" # Data prefix\n", " else: # MC prefix\n", " prefix = \"MC/mc_\"\n", " fileString = val\n", "\n", " # start the clock\n", " start = time.time()\n", " print(\"\\t\"+val+\":\")\n", "\n", " # Open file\n", " tree = uproot.open(fileString + \":analysis\")\n", "\n", " sample_data = []\n", "\n", " # Loop over data in the tree\n", " for data in tree.iterate(variables + weight_variables + [\"sum_of_weights\", \"lep_n\"],\n", " library=\"ak\",\n", " entry_stop=tree.num_entries*fraction):#, # process up to numevents*fraction\n", " # step_size = 10000000):\n", "\n", " # Number of events in this batch\n", " nIn = len(data)\n", "\n", " data = data[cut_trig(data.trigE, data.trigM)]\n", " data = data[cut_trig_match(data.lep_isTrigMatched)]\n", "\n", " # Record transverse momenta (see bonus activity for explanation)\n", " data['leading_lep_pt'] = data['lep_pt'][:,0]\n", " data['sub_leading_lep_pt'] = data['lep_pt'][:,1]\n", " data['third_leading_lep_pt'] = data['lep_pt'][:,2]\n", " data['last_lep_pt'] = data['lep_pt'][:,3]\n", "\n", " # Cuts on transverse momentum\n", " data = data[data['leading_lep_pt'] > 20]\n", " data = data[data['sub_leading_lep_pt'] > 15]\n", " data = data[data['third_leading_lep_pt'] > 10]\n", "\n", " data = data[ID_iso_cut(data.lep_isLooseID,\n", " data.lep_isMediumID,\n", " data.lep_isLooseIso,\n", " data.lep_isLooseIso,\n", " data.lep_type)]\n", "\n", " # Number Cuts\n", " #data = data[data['lep_n'] == 4]\n", "\n", " # Lepton cuts\n", "\n", " lep_type = data['lep_type']\n", " data = data[~cut_lep_type(lep_type)]\n", " lep_charge = data['lep_charge']\n", " data = data[~cut_lep_charge(lep_charge)]\n", "\n", " # Invariant Mass\n", " data['mass'] = calc_mass(data['lep_pt'], data['lep_eta'], data['lep_phi'], data['lep_e'])\n", "\n", " # Store Monte Carlo weights in the data\n", " if 'data' not in s: # Only calculates weights if the data is MC\n", " data['totalWeight'] = calc_weight(weight_variables, data)\n", " # data['totalWeight'] = calc_weight(data)\n", "\n", " # Append data to the whole sample data list\n", " sample_data.append(data)\n", "\n", " if not 'data' in val:\n", " nOut = sum(data['totalWeight']) # sum of weights passing cuts in this batch\n", " else:\n", " nOut = len(data)\n", "\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", " frames.append(ak.concatenate(sample_data))\n", "\n", " all_data[s] = ak.concatenate(frames) # dictionary entry is concatenated awkward arrays" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "colab": { "base_uri": "https://localhost:8080/", "height": 704 }, "id": "o65VsV6pk1Bt", "outputId": "5de2102a-23c1-4058-b2b0-22a51e81ac2f" }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsoAAAHnCAYAAAChTHs1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAAsTAAALEwEAmpwYAACxYUlEQVR4nOzdeVxU5f4H8M9hEcQFVHDDAM0dTVM0tJJxQUvDpV9qiWU3bey23DK792ZY4k2u1bVsL7FcStLE1CRLw+VQqaS4ZeSWhqYp4gKlKLI8vz9gxhlmOzOc2eDzfr3m5cxznvM8z5kZme+ceZ7vkYQQICIiIiIiYz7uHgARERERkSdioExEREREZAYDZSIiIiIiMxgoExERERGZwUCZiIiIiMgMBspERERERGY4NVCWJClEkqRVkiQdkiTpoCRJ/SRJaipJUqYkSUer/m3izDEQERERETnC2WeU3wKwQQjRGUAPAAcBPA9gsxCiA4DNVY+JiIiIiDyK5KwLjkiSFAxgH4B2wqATSZIOA9AIIc5IktQKgCyE6OSUQRAREREROcjPiW23BVAAYLEkST0A7AbwNIAWQogzVXXOAmhRfcfAwEDh6+urfxwaGoqwsDCHB1JUVITg4GCH9zfXno+PD44fP4527dqhUaNG+m1//fWX2XJXja2goKBGz1V1ao7PGa8DnzvPaE/N587Tj1XN9jz5Pefp7dWl586T/78Cnn2snvzcefqx1vbnrqCgAOfPnwcAFBcXXxNC1DdbUQhh9gbA39I2JTcAMQDKANxW9fgtAC8DKKxW71L1fYOCgoSaHn30UVXbGzFihAgNDRVbt241Kt+6davZcleOrXfv3qq2p+b41D5WPnee056az52nH6ua7Xnye87T26tLz50n/38VwrOP1ZOfO08/1rr03AEoEBbiWWtzlE9LkvSRJEmDJUmSHAjWTwE4JYT4serxKgC9AORXTblA1b/nHGjbLgkJCaq299133yE9PR0ajUZfJssyxo4da1Lu6rGpTc3xqX2sfO48pz01efqx8rnznPbU5MnH6snPG+DZx+rJz52nH2tdeu4AFFraYHGOsiRJzQDcB+B+AB0AfAFguRAiW2mvkiR9D2CKEOKwJEnJABpUbboghHhFkqTnATQVQvzLcL8GDRqIK1euKO3G5WRZViVIdoaYmBjk5OS4dQzeis+d4/jcOYbPm+P43DmOz53j+Nw5zpOfO0mSdgshYsxts3hGWQhxQQixQAgxEEBfAMcBzJck6ZgkSSkK+34KQJokST8B6AngvwBeARAvSdJRAEOqHhsJDQ1V2Lx7KAmSZVl2+bgAQKvVuqXf2oDPneP43DmGz5vj+Nw5js+d4/jcOc5bnzvFWS8kSWoI4F4AzwJoJYQwWYSnlpiYGOGp3zoMWQuSx44di4KCAvcNjoiIiIhscuiMctWOgZIkjZUkaTWAXwEMQmXe49bqD9O72AqS09PT3Tc4IiIiIqoxi+nhJEn6DJVTI7IApAGYIIS45qqBeTIlQbK75yoTERERUc1Yy6O8AcBUIcRfrhqMt2CQTERERFT7WVvM94kQ4i9JklpIkvSxJEkbAECSpK6SJE123RA9D4NkIiIiotrP6hzlKksAbATQqurxEQDPOGk8XsGTs14QERERkTqUBMqhQoiVACoAQAhRBqDcqaPyErbmKhMRERGR91ISKF+puviIAABJkmIBFDl1VF6AWS+IiIiIajdri/l0ngWwDsDNkiRtAxCGyiv21VnMekFERERU+9kMlIUQeyRJigPQCYAE4LAQotTpI/NgDJKJiIiIaj+LUy8kSWosSVIHQD8vuSuAXgAekCTJaVflA4CioiJotVpkZGQ4sxuHMUgmIiIi8m4ZGRm6S2sHW6pj8RLWkiSlAtguhFhS9fhXAF8DCAJQJoR4TO0B63jLJawB69MwGDQTERHVXkIIlJcb5zfw9fWFJEke0R4p4+glrPsAWGrw+C8hxD+EEFMAdFNzgN6KWS+IiIjqrqysLPj7+xvdli5dantHF7VHNWdtjrKfMD7d/KDB/RDnDMd7MOsFERFR3da7d2/s2rXLqKxt27Ye0x7VnLVAuUKSpJZCiLMAIIT4GQAkSQpHVU7luopZL4iIiKhRo0aIiTH7i71HtEc1Z23qxf8AZEiSNECSpEZVtzgAa6u21VnuDJIfffRRSJKEadOmGfUtSZLN28MPPwyg8sqCd9xxh+I+4+PjIUkS3nrrLbPbKyoqsHjxYvTt2xdNmjRBgwYNcPPNN+P+++/Hzp07FfWxa9cu/N///R9atGiBgIAAREVF4fHHH8fp06cVj9Mdqj/39evXR5s2bTB8+HB89NFHuH79ukPt5uXlITk5GcePH1d5xETkDEuWLDH6W+Dr64vw8HCMGzcOhw8fdkqfycnJkCQJZWVlTmnfU2k0GpufsyUlJahXr57Nz8WDBw8q7vfPP/9EcnKyyT6Wyh1tryZtKjVu3DhERkZa3P7rr7+iXr16eOwxZcvR7GnPWXWdxeIZZSHEMkmSzgOYAyAalRccyQXwkhDiG6eNyAu4K0i+evUqVq5cCQD47LPP8L///Q9+fn7o1asXduzYoa935swZ3HvvvZgxYwZGjhypLw8LC7O7z1OnTmHLli0AgE8++QRPP/20SZ3nnnsOb7/9Np555hkkJyfD398fR44cwZo1a/Djjz+ib9++Vvv49NNP8be//Q133HEH3nrrLbRu3RoHDx7Ea6+9hlWrVmHTpk245ZZb7B67K7399tvo06cPSktL8ccffyAzMxNPPPEE3n33XWRmZtr93Ofl5WH27Nm444470K5dOyeNmsizPDnhZ7f2/+5nNV9+k56ejjZt2qC8vBzHjh3Dyy+/jMGDByM3NxfBwRYX1pMTfPfdd2bLly9fjrfffhsDBgxAx44dFbeXk5OD2bNn495771VU7mh7NWlTqejoaKxatQqXL19Gw4YNTbbPmDED9evXx3/+8x/V2/vggw+cUtdZrOZRFkJsALDBab17KSVBsjOyXqxduxZ//vknhg8fjq+//hobNmzAPffcg8aNGyM2NlZfLy8vDwDQrl07o3JHfPrpp6ioqND3+fPPP6NbtxsfJlevXsV7772Hp556CvPmzdOXx8fH44knnkBFhfVZOocOHcKjjz6K0aNHY+XKlfDxqfyRY8CAAbjvvvtw22234b777kNubi78/f1rdCzO1KVLF6Pnevz48Zg8eTIGDRqERx55xGNTHRKRunr27In27dsDAG6//Xa0bt0a8fHx2L59O+6++243j85+JSUlCAgIcPcw7BYQEGD28+/LL7/E+++/jz59+uCrr76Cr6+v4jb37t2LgIAAdO3aVVG5o+052mZUVBQefvhhJCcn26zbrVs3CCFw8OBB9OnTx2hbdnY2Vq1ahblz56J58+aK+ranPWfVdRYll7AmC1yd9WLp0qVo0qQJlixZgvr167tkJezSpUsRHR2NN998U//Y0JUrV3D9+nW0bNnS7P66wNeSt956C+Xl5XjnnXdM6jZr1gz//e9/cfToUaxevVpfHhUVhYkTJ2LhwoVo3749AgMD0atXL2zdutWk/aysLAwePBiNGjVCgwYNMGzYMPz8s/FZK91UlE2bNqFXr14ICgpCt27dsGbNGqtjt6Vfv3547LHH8NVXX+HYsWP68nfffRf9+vVD06ZNERISgtjYWKxfv16/XZZlDBw4EMCNaS+SJEGWZQDAihUrMGjQIISFhaFhw4a49dZbuSqayEM1btwYAFBaWnmdrl9//RUPPvgg2rZti/r166Ndu3b4+9//jkuXLpnsu3//fowZMwbNmjVD/fr10alTJ8ydO9dqfxs2bEDDhg3x5JNP6k9ULF++HJ07d0ZgYCC6d++OdevWmZ2+oJvK8fPPP2PYsGFo2LAhxo0bp2+3X79+qF+/PoKDgzF69GiTKSUPP/wwoqKiTMZUvS9dP0ePHsWIESPQsGFDREZG4j//+Y/JyZUVK1agc+fOCAgIQHR0dI3+LmdmZmL8+PHo0qULNmzYgEaNGinet0uXLnjuuedQUlICf39/SJKE//u//7NY7mh7trapRXfC65dffjHZ9s9//hORkZF45plnnNKes+o6CwNlB7k668Uff/yBTZs2Yfz48QgLC8Po0aORkZFh9o+rWn788UccPnwYDz74IDp06IB+/fohLS3NKMdjaGgo2rZti3nz5uHDDz/EyZMn7epj8+bNiImJQatWrcxuHzFiBHx8fPTTP3RkWcYbb7yBlJQUrFixAgEBAbj77ruN/nCvX78egwcPRsOGDbFs2TJ89tln+Ouvv3DnnXfi999/N2rv2LFjePrpp/Hss89i9erVaNWqFcaOHYtff/3VruOpbvjw4QCAbdu26cvy8vIwZcoUpKen4/PPP0dMTAzuuecebNhQ+eNNr1698N577wGonNKxY8cO7NixA7169QIAHD9+HPfddx/S0tKwdu1aJCQkYMqUKfjwww9rNFYiqrny8nKUlZWhpKQEBw8exAsvvIDmzZvrPyf++OMP3HTTTXjzzTexceNGvPTSS9i8ebP+b4XOzp070a9fPxw7dgzz58/H+vXr8eyzz+LUqVMW+/7kk08wcuRIPP/883j33Xfh4+ODzMxMJCYmonPnzli9ejWee+45PPPMMzhy5IjFdkaNGoW4uDisW7cO06ZNw4YNG/QB7eeff44PPvgAP//8M+64444arSMZM2YMBg0ahLVr12L06NGYNWuW0Zf+TZs2YcKECejQoQNWr16Nf/7zn3j66acdmvP9ww8/YPTo0YiIiEBmZiaaNm1q1/6ffPIJ2rVrh4SEBP3f5DfeeMNiuaPt2dqmlvbt2yMgIAC5ublG5WvXrsUPP/yAuXPnIjAw0CntOauu0wghPO7Wu3dv4cm2bt0qQkNDxdatWxWVq+HVV18VAMT27duFEEJs2LBBABAffPCBSd3ffvtNABALFy4021ZcXJy4/fbbbfb597//Xfj4+IhTp04JIYT48MMPBQDxzTffGNXbsWOHiIyMFKicxy5at24tHnnkEfHjjz/a7CMwMFDcf//9Vuu0aNFC3H333frHkZGRwt/fX5w8eVJf9ueff4omTZqIiRMn6stuvvlmMWjQIKO2ioqKRLNmzcTTTz+tL4uLixN+fn7iyJEj+rL8/Hzh4+MjUlJSrI5t69atAoDIzMw0u/3QoUMCgHjllVfMbi8vLxelpaUiPj5ejBw5UnG71fefMmWKuOWWW6zWJfJ0TzxwwK23mli8eLH+b6DhrXXr1mLnzp0W9ystLRXff/+9ACD27NmjL7/zzjtFmzZtxJUrVyzuO2vWLAFAlJaWildffVX4+fmZ/N3v16+fiI6OFhUVFfqynJwcAUDExcWZbe/NN980Ku/du7do3769KC0t1ZcdP35c+Pn5iWnTpunLJk2aJCIjI03GGRcXZ9SXrp9FixYZ1evWrZuIj4/XP+7fv7/o0qWLKC8v15ft2LHD7NitycnJEY0bNxYRERFGnxv2KCkpEfXq1RPz589XVO5oe0rbrKioEKWlpUa3yMhI8eKLLxqVlZWVWWyjR48eYsSIEfrHpaWlolOnTuK2224zer8oZU97zqrrKAA5wkJMqviMsiRJd0iS9KwkSUNVitG9ljuyXixdulR/VhcAhgwZgtatWzvtJ/eSkhL9T/zh4eEAKufdBgQEmPQZGxuLw4cP45tvvsH06dMRFRWFpUuXol+/fvjkk0+cMr7Y2FjcdNNN+seNGjXCiBEj9Isajx49imPHjiExMRFlZWX6W1BQEPr162eyyKNDhw7o0KGD/nHz5s3RvHlzu8+QVyeqUpEbXlVp9+7duOeee9CiRQv4+fnB398fmZmZis+SHD16FA888ADCw8P1Cek/+ugjp62sJyLl1qxZg127dmHnzp1Yu3YtunbtiuHDh+uzF1y/fh3//e9/0blzZ9SvXx/+/v648847AUD/f7i4uBjbtm1DYmIigoKCbPY5bdo0zJo1C6tWrcKUKVP05eXl5cjJycH//d//Gf0N6t27t9XcvGPGjNHfv3LlCvbs2YPx48fDz+/Gsqa2bdvi9ttvR1ZWlsJnxtSIESOMHnfr1k3/N7e8vBy7du3CfffdZzQtLzY21uz0Dktyc3MxbNgwBAUFYfPmzUafG/bIzc3F9evX9b/s2Sp3tD2lbZq7MMmJEyfw8ssvG5UNHjzYYhvdunUzmtKwcOFCHD58GG+88YZDVwK0pz1n1XUGi4GyJEk7De4/CuBdAI0AzJIk6Xmnj8yDuTpIzsnJwS+//IJ7770XhYWFKCwsxF9//YV7770X2dnZVn9Cc5RuWseYMWP0fQLAsGHD8OWXX+LPP/80qh8QEIC77roL8+bNw7Zt2/DLL7+gZcuWePbZZ63206ZNG/3iQ3OuXLmCgoICkz9uLVq0MKnbokUL/c+A586dAwBMnjzZ5I/JV199hQsXLhjta+5nuICAAFy7ds3q+G3RTfHQTS35/fffMXjwYFy8eBHvvPMOtm/fjl27duGuu+5S1Nfly5cRHx+P/fv345VXXsH333+PXbt24ZFHHkFJSUmNxkpENdetWzfExMSgT58+GDVqFNatWwchhH6B1YwZM5CcnIyJEydi/fr12Llzp34Nhu5vwKVLl1BRUYE2bdoo6nP58uXo1q0bhgwZYlR+/vx5lJaWml3oZO5vqI7hVLhLly5BCGF2elzLli1x8eJFRWM0p/rfXcO/ubqxW/pbr8SxY8cQHx8PoHIah26RpSP27NkDSZLQs2dPReWOtqe0Td2FSQxvrVq1wqOPPmpUtmDBAottREdHIy8vD8XFxbh8+TJmz56NsWPHon///nYdiyPtOauuM1jLemGYYkALIF4IUSBJ0jwA2QBecerIPJirs17ozuC++uqrePXVV022f/LJJ5gzZ44qfVXv84knnsATTzxhsn3lypVGZy6q69ixI8aPH4/58+fj3LlzFlekDh48GB9//DHOnDlj9g/x+vXrUVFRgUGDBhmV5+fnm9TNz8/Xn/1u1qwZAGDu3LkmHx4AUK9ePYtjV5NukZ4ub/WGDRtQVFSElStXGn0IFhcXK2pvx44dOHHiBL7//nujXNh1LY8qkbfQLdj76aefAFQuTnvooYcwc+ZMfZ3Lly8b7dOkSRP4+Pgonv+7efNmDB06FHfffTe+/vprfRqt0NBQ+Pv7608cGMrPz0dERITZ9gzP0jVp0gSSJOHs2bMm9c6ePWsU7AYGBprNHX/hwgX932SldGO39LfeWm5d4MZJiStXrmDz5s2Ijo62q//q9u7di5tvvlm/ONNWuaPtKW3T3IVJ6tWrh9atWyu+YIlhRokvv/wShYWFZmMMpexpz1l1ncHa1AsfSZKaSJLUDIAkhCgAACHEFQBO/VQuKiqCVqv1+JRarsh6cf36dSxfvhy33XYbtm7danLr2bMnPv30U/1P/Go4d+4cNmzYgFGjRpnts2XLlvpAurS01OTsrM6hQ4f0K6Qtefrpp+Hj44OnnnrKZLXzxYsX8cILL6B9+/YmuSSzs7ONFuT99ddfWL9+vX5qSqdOnRAVFYXc3FzExMSY3FyRl3nHjh348MMPMXr0aP3PnLqA2DDV3ZEjR4wW+wHQp2O6evWqUbm5/S9duoQvv/xS/QMgohorLi7GsWPH9LnUi4uLTVJdLl682OhxUFAQ7rjjDixbtszkb4A50dHRkGUZR48exd13360PvH19fRETE4MvvvjC6DNi9+7d+O233xSNv0GDBujduzfS09ONFnKfOHEC27dvN/rsi4yMRH5+PgoKCvRlx44dc2hamK+vL/r06YNVq1YZfTb8+OOPVn+FBCo/w4YMGYKCggJ89dVXqlzp7pdffjGbqs1SuaPt1aRNe+kySmzevBlvvPEGnnrqqRpdLtue9pxV114ZGRnQarUAYDFQsXZGORjAbgASACFJUishxBlJkhpWlTlNcHAwUlNTndlFjbkq68X69etx4cIFvP7662bPUE+dOhV///vfjVKKKXHhwgWsWrXKpPyWW27B+vXrUVZWhmnTpiEuLs6kzqRJk/Daa6/h+PHjaNy4MaKiojB+/HgMGTIEbdq0wYULF7BixQp88803+Ne//mU1B2eXLl2wYMECTJkyBYMHD8Zjjz2GVq1a4dChQ3jttddQWFiIzMxMkw+WFi1aYOjQoUhOTkZAQABeffVVXLlyBS+++CKAyjMi7733HkaNGoXr169j3LhxCA0NRX5+PrZv346IiAib00LscfDgQTRs2BBlZWU4c+YMvv32W3z66afo2rUrFi5cqK83ZMgQ+Pn54aGHHsL06dNx5swZzJo1CxEREUYfBh07doSfnx8WLVqEpk2bIiAgAJ06dUL//v3RuHFjPPHEE5g9ezauXLmCOXPmIDQ0FEVFRaodDxE5Zt++fTh//jyEEDhz5gzeffddXLx4EU899RQA4K677sLSpUvRvXt3tG/fHqtXr8b27dtN2pk3bx7i4uLQr18/TJ8+HW3atMHx48exb98+vPPOOyb1u3Tpov8cGDZsmD792ezZszF06FCMGTMGWq0W58+fR3JyMlq2bGkzfafOyy+/jBEjRuCee+7B448/jsuXL2PWrFkIDg7G9OnT9fXGjh2LF198ERMnTsSzzz6L8+fPY+7cuQgNDXXoudSNffTo0Zg6dSoKCgowa9Ysi+lIgcrpevHx8Thy5Ahmz54Nf39/ZGdnm9Rr166d/pfOvLw8tG3bFrNmzbKYgzgkJAR79uzBxo0bERwcjA4dOqBZs2YWy221aWk/W9vUFBUVhQYNGmDWrFlo2LCh0a8chpQ8P/a058y69kpISEBCQgIWLlxo+QPU0io/SzcAQQDa2rufPTdmvbhh1KhRolGjRhZXPhcWFor69euLSZMm6cuUZL2AmdXZAMT//vc/0aNHD3HzzTdbXE16+PBhAUDMmjVLlJSUiNdee03Ex8eL8PBw4e/vLxo1aiRiY2PFggULFK9I3bFjhxg9erQIDQ0V/v7+IiIiQkydOtXsCuXIyEiRmJgoFi5cKNq1ayfq1asnevbsKTZv3mxSd/v27WLEiBEiJCREBAQEiMjISDF+/Hh99hDd82EuC0hkZKTR82qOLjuF7hYQECBat24t7r77bvHRRx+JkpISk30+//xz0alTJxEQECC6du0qli9fbna1+Icffijatm0rfH19BQD9+2rz5s2iZ8+eIjAwULRr10689dZb+lXkRN6stmW9CAsLEwMHDhQbNmzQ1ysoKBDjx48XISEhIiQkREyYMEHs3LlTABCLFy82anPPnj3innvuEcHBwSIwMFB06tTJKIOOYdYLnSNHjojw8HARGxsrioqKhBBCpKWliY4dO4p69eqJrl27itWrV4uePXuK0aNHG/Vnrj2db775RsTGxorAwEDRuHFjMXLkSHHo0CGTemvWrBHR0dEiMDBQ3HLLLWLjxo0Ws15U78fc38HPPvvMZOzV2zOUmZlp8fPN8LZx40b9Pj///LPFLFI6Bw4cEH379hWBgYECgPj++++tlttq09J+trZZExkZKWbNmqWork6fPn0EAPHuu+9arKPk+bGnPWfXdQSsZL2QhIo/2aslJiZG5OTkuHsYFoWFhbntMtZU+e1S97MkERHZ59SpU2jfvj2SkpL0v8LVZampqUhKSsKJEycUZRlxV5vuUpuOxRJJknYLIczO0bF6CWsrDX4lhLinZsPyXgySiYjIG1y9ehXPPvsshgwZgtDQUBw/fhyvvfYagoKCrC7IrkuysrIwbdo0VYNAZ7TpLrXpWBzhUKAM4FFVR+FlXJ31goiIyBG+vr44e/YsnnzySVy4cAENGjTAnXfeifT0dItXRK1r0tLSvKJNd6lNx+IIhwJlIcQZtQfijWwt6DNc+UvqsbXimYiIKtWrVw9r1qxx9zCIvJa1C440liRpriRJn0qSNKHatvedPzTP5qqsF0RERETkHtZywyxGZRq4LwDcL0nSF5Ik6fJ8xTp9ZB5MSZDMaRdERERE3s1aoHyzEOJ5IcRaIcRIAHsAbKm6AEmdxiCZiIiIqPazNkc5QJIkHyFEBQAIIVIkSToN4DsADV0yOg/FIJmIiIio9rN2RjkDwCDDAiHEEgDTAZhezL0OUZr1goiIiIi8l8UzykKIf1ko3wCgg9NG5EWY9YKIiIio9lJ2oXcywawXRERERLUbA2UHuCPrRUFBAeLj49GkSRM88sgjqratxODBg10e/Ofm5qJjx44u7ZOIiIhIx1oe5dauHIihoqIiaLVaZGRkuGsIVrkj68XcuXPRoUMHXLp0CYsWLVK9fVtyc3Nxyy23OLWPNm3aYO/evfrH0dHROHLkiFP7JCIioropIyMDWq0WAIIt1ZGEEOY3SNLXAJoCkAFsAPCDEKJM/WGaiomJETk5Oa7oyiHVL0/tiqwXt9xyC9566y0MHDjQrv2EEBBCwMfH8R8PLl68iDZt2uDy5cs1asea8+fPo1WrVrh8+TICAgJs70BERESkAkmSdgshYsxtsxj1CCGGA9CgMlAeAyBbkqTVkiRpJUmKcMZAvYUrs15cv34dwcHBOHDgABISEtC9e3cAld+CevTogUaNGiE2Nha5ubn6fe68807MnDkTffv2Rf369XHhwgWTdg8cOIBhw4ahWbNmaN++Pb7++mv9tqKiIkyaNAmNGzdGt27dkJ6ejq5du8LHxwe//vormjRpYtTWkCFDsGLFCv3j5cuXIzo6GkFBQbj55pv1z8WXX36JXr16ITg4GK1atcKbb74JAPj1119x0003oaKiAs2aNUOzZs1QVlaGBx98UF8HAD7++GN06tQJwcHBGDp0KE6fPq3f1qpVK7z99tu49dZb0bBhQ9x7770oLy93+HknIiIi0p9xVHID0BbA4wDWAdhpz7723Hr37i28wdatW0VoaKjYunWr2XK15ObmiubNm+sfb9y4Udx0001i165doqKiQrzyyiuiV69eQgghKioqROPGjcXtt98u8vLyxLVr10zaO3TokAgLCxMrV64UZWVlIjMzUzRr1kxcuXJFVFRUCI1GI5544glx9epVcfDgQdG0aVPxt7/9TQghxMqVK8WAAQOM2gsNDRW5ublCCCHmzZsnunbtKnbv3i3Ky8vFTz/9JH777TchhBBfffWVOHr0qKioqBDff/+98Pf3F/n5+UIIId59910xbtw4o3a7d+8utmzZIoQQYuHChSI6OlocOXJElJWViccff1zce++9Qgghzp07JwCICRMmiAsXLohz586JsLAw8f3339f0qSciIqJaDkCOsBT7Wtpg6wagnqP72rp5Q6BsK0iuXl4TaWlpIj4+Xv/41ltvFenp6frHJ0+eFJIkibKyMnH8+HEhSZI4fvy4xfbGjBkjkpKSjMo6duwo9u3bJ9atWyduvvlmUVZWpt8WHx8v3njjDSGEEC+88IJ48skn9dtOnz4tAgICRFlZmTh37pxo1KiR2Ldvn6Ljuummm8Thw4eFEEI8+uij4r///a9+W0lJifD39xcXLlwQ165dEy1bthS7du3Sb//uu+9ERESEEEKIzMxM0bRpU3H16lX99q5duzJQJiIiVVVUVAghhJg1a5b47bff9I/Ju1kLlB2ecCqEqLMXHXF11ot9+/ahR48eACrn8u7fvx8jR47Ubz9//jyaN28OX19f7N+/HzExMWjbtq3ZtioqKrBhwwa8/fbbCAkJ0d9+//13BAUFYf369RgzZgx8fX31++Tn5+sX8u3du1c/FgDYv38/oqOj4evri02bNqF79+5G23WuXLmCf/7zn+jUqZO+z/z8fERFRZkcIwD88ssvaNGiBZo2bYoDBw7A19cXMTE3pg/p5jQDwE8//YQhQ4YgMDAQQOV0laNHj6Jr1652Pc9ERETWrFy5Es8//zwKCwuxc+dOTJw4Efn5+e4eFjkR08M5wNVZL/bv368PIgsKChAYGIh69erpt3/55Ze444479HV79+5tsa2rV6/i6tWrOHPmDAoLC/W34uJidOjQQR906xw8eBAHDhzQB8q//PKLfp40AGRmZuq3Xbx4ESEhIWb7ffzxx3Hx4kVs3boVly5dwvLly9GlSxfUq1cPFRUV+Pnnn9GzZ0+LxxwcbLwgtfox33rrrfptP//8M1q2bImmTZtafB6IiNS0ZMkSSJKkv/n6+iI8PBzjxo3D4cOHndJncnIyJElCWZlL1tl7DI1GY/NztqSkBPXq1TN6TczdDh48aFfd8ePH495778WiRYvw/vvv46OPPkKLFi1cc+DkFnYFypIkNZEkybk5wryAq1PDGQaN7dq1Q1BQEFavXo2ysjKsWrUK77//PmbNmqWvaxhwVtegQQN06tQJb731FkpLS1FaWoq9e/fi0KFDAIBOnTrhiy++wJ9//om8vDw89NBDaN68OcLCwgBUTtUpKioCAPzwww/48MMP9YHyrbfeih9++AH79++HEAJHjx7FwYMHAVSeMe7Tpw9atWqFH3/8EU888YR+nLrgvaKiwuwx9+jRAydOnMDOnTtRUlKCDz74AJmZmZg+fTqAyjPKhsdc/aw3EXkJSXLvTQXp6enYsWMHvvvuO8ydOxd79+7F4MGD9X83yXW+++477Nixw+T2j3/8AwAwYMAAfa5+pXXT09OxevVqPPLII3j88ceh1Wp5Rrm2szQnQ3dDZdaLxqhMFfcbgB8BvGFrv5rcvGGOso61ucpqOHPmjKhXr564fv26vuzbb78VnTp1Eg0aNBD9+vUTO3bs0G9r166dyM7Ottrm3r17xe233y4aNWokmjRpIuLi4sShQ4eEEEIUFBSIgQMHikaNGolbb71VPProo0bzoz/66CPRunVr0b9/fzF9+nQRHR0tNm3apN/+2muviTZt2ogGDRqI6OhosXv3biGEEF988YVo2bKlaNq0qXjkkUdEQkKCmD9/vn6/xx57TDRq1EiEh4cLIYQYNGiQWLFihX77smXLREREhGjYsKGIj4/Xj7e0tFQEBASIM2fO6Os+8cQTYubMmYqfYyLyEIB7bzWwePFiAUAcPXrUqDwzM1MAEF9//XWN2jdn1qxZAoAoLS1VvW0dcwvC3S0uLk7ExcU5tO/atWuFn5+f6NOnj/jzzz/trss5yrUTarKYD8Deqn+nAJhddf8nW/vV5OYtgbKrsl4QEdUJtTBQ/vHHHwUA8eWXXwohhDh69KiYOHGiiIqKEoGBgaJt27biscceExcvXjRpc9++fWL06NGiadOmIjAwUHTs2NFo0bO5QPmbb74RDRo0EE888YQoLy8XQgjx2WefiU6dOomAgADRrVs38eWXX5oNNnXtHThwQAwdOlQ0aNBAjBw5Ut9ubGysCAwMFI0bNxajRo3Sn7DQmTRpkoiMjDQ5jup96fo5cuSIGD58uGjQoIGIiIgQs2fP1o9ZZ/ny5aJTp06iXr16omvXrmL16tUOB8rffvutCAgIEN27dxcXLlxQrS55P2uBspKpF36SJLUCMA7AV+qcx/Z+Shb0ERFR3VJeXo6ysjKUlJTg4MGDeOGFF9C8eXP958Qff/yBm266CW+++SY2btyIl156CZs3b8bw4cON2tm5cyf69euHY8eOYf78+Vi/fj2effZZnDp1ymLfn3zyCUaOHInnn38e7777Lnx8fJCZmYnExER07twZq1evxnPPPYdnnnnG6lVPR40ahbi4OKxbtw7Tpk3Dhg0bMGLECDRs2BCff/45PvjgA/z888+44447jPLZ22vMmDEYNGgQ1q5di9GjR2PWrFlYunSpfvumTZswYcIEdOjQAatXr8Y///lPPP300w7N+f7hhx8wevRoREREIDMz0+oaFnvqUh1gKYLW3QDcB+AnAO9XPW4H4Atb+9Xk5ulnlF2ZGo6IqM6oBWeUq99at24tdu7caXG/0tJS8f333wsAYs+ePfryO++8U7Rp00ZcuXLF4r6GZ5RfffVV4efnJxYuXGhUp1+/fiI6OtpoikBOTo4AYPGM8ptvvmlU3rt3b9G+fXujM9fHjx8Xfn5+Ytq0afoye88oL1q0yKhet27djKb69e/fX3Tp0sXoLPOOHTvMjt2anJwc0bhxYxERESFOnjypWl2qPVDDM8pnhBC3CCEerwqsjwN4Q8VY3eu4OusFERF5hzVr1mDXrl3YuXMn1q5di65du2L48OH6hc3Xr1/Hf//7X3Tu3Bn169eHv78/7rzzTgDQnyktLi7Gtm3bkJiYiKCgIJt9Tps2DbNmzcKqVaswZcoUfXl5eTlycnLwf//3f5AMFiv27t3bYgpRoPJMr86VK1ewZ88ejB8/Hn5+fvrytm3b4vbbb0dWVpbCZ8bUiBEjjB5369YNJ0+e1I99165duO++++DjcyNUiY2N1acVVSI3NxfDhg1DUFAQNm/ejJtuukmVulR3+NmugncA9FJQVmcwSCYiInO6deuG9u3b6x8PHToUN910E5KTk/H5559jxowZeOedd/DSSy+hf//+aNSoEU6dOoV7770X165dAwBcunQJFRUVaNOmjaI+ly9fjm7dumHIkCFG5efPn0dpaalRyk8daynNdDnqdWMRQhiV6bRs2RInTpxQNEZzqk9pCAgI0D8HurGbG6fSdGzHjh1DfHw8gMppHIavS03qUt1i8YyyJEn9JEmaDiBMkqRnDW7JAHwt7aeGoqIiaLVaZGRkOLMbhykJkmVZdvm4iIjIs9SvXx/t2rXDTz/9BABYsWIFHnroIcycORODBg1Cnz59TPLPN2nSBD4+Porn/27evBknT57E3XffjcuXL+vLQ0ND4e/vj3PnzpnsYy2lmeHZ5yZNmkCSJJw9e9ak3tmzZ42C3cDAQFy/bnotsgsXLig6DkO6sZsbp5J0bL///jsGDx6MK1euYMOGDYiOjlalLtUuGRkZ0Gq1ABBsqY61qRf1ADRE5VnnRga3P1E5b9lpgoODkZqaioSEBGd2U2O2FvQREVHdVlxcjGPHjulz0RcXF8Pf39+ozuLFi40eBwUF4Y477sCyZctw9epVm31ER0dDlmUcPXrUKFjWXdH0iy++0K05AgDs3r0bv/32m6LxN2jQAL1790Z6ejrKy8v15SdOnMD27duNPvsiIyORn5+PgoICfdmxY8ccWnzn6+uLPn36YNWqVUY59n/88Ufk5eVZ3ffcuXMYMmQICgoK8NVXXxld1bUmdan2SUhIQGpqKgBYTHRuceqFECILQJYkSUuEEI7/tlJLMesFERFVt2/fPpw/fx5CCJw5cwbvvvsuLl68iKeeegoAcNddd2Hp0qXo3r072rdvj9WrV2P79u0m7cybNw9xcXHo168fpk+fjjZt2uD48ePYt28f3nnnHZP6Xbp0gSzLGDhwIIYNG4YNGzagUaNGmD17NoYOHYoxY8ZAq9Xi/PnzSE5ORsuWLY3m/lrz8ssvY8SIEbjnnnvw+OOP4/Lly5g1axaCg4P1F34CKtfvvPjii5g4cSKeffZZnD9/HnPnzkVoaKhDz6Vu7KNHj8bUqVNRUFCAWbNmoWXLlhb3uXLlCuLj43HkyBHMnj0b/v7+yM7ONqnXrl07NGjQQHFdc9NXqI6wtMpPdwPQEUAqgG8BbNHdbO1XkxuzXhAR1UG1LOtFWFiYGDhwoNiwYYO+XkFBgRg/frwICQkRISEhYsKECWLnzp0CgFi8eLFRm3v27BH33HOPCA4OFoGBgaJTp07ilVde0W83l0f5yJEjIjw8XMTGxoqioiIhhBBpaWmiY8eORrmIe/bsKUaPHm3Un7ULmFTPozxy5EiTPMpCCLFmzRoRHR0tAgMDxS233CI2btxoMetF9X7MZc347LPPTMZuLY+y7gIvtm4bN260qy7VbrCS9UISBj/HmCNJ0n4AHwLYDUD/u4sQYndNAnRrYmJiRE5OjrOar7GwsDAu6CMAwOeff44XXngB165dwzvvvIN7773X3UMiIrLq1KlTaN++PZKSkvDiiy+6ezhEbidJ0m4hhNl5N0oC5d1CiN5OGZkFnh4oy5IEjeFjAGMBpANG5bDx3JJ3y87ORv/+/eHn54d69eoBAI4cOYLWrVu7eWRERJWuXr2KZ599FkOGDEFoaCiOHz+O1157Dfn5+cjNzTWbzYKorrEWKCuZoJQhSdLjkiS1kiSpqe6m8hi9isbgvgzzQbLsstGQu6xevRpCCLz22mv47LPPcOXKFXzzzTfuHha5wMsvv4yOHTvCx8cHa9eudfdwiCzy9fXF2bNn8eSTTyI+Ph7PPvssOnTogO+++45BMpECSvIoT6r6958GZQKVV+ir02RYDpLHAigw2YNqE92lZPv06YPevXtj7969aNeuzv+3qBPi4+ORmJiIRx55xN1DIbKqXr16WLNmjbuHQeS1bAbKQgjLl++pw2RYD5JdmfPi1VdfxaJFi3DkyBGH9rc1/YbM+/PPPwFU5hkNDAxEz5493TsgcpnY2Fh3D4GIiFzA5tQLSZKCJEmaKUlSatXjDpIk3eP8oXkuGbaDZA1c5/PPP8ePP/7ocJYRe506dQpPPfUU+vXrh6CgIEiSZDGv5caNGzFo0CC0bNkSAQEBaNOmDcaNG4dffvnFah+SJNm82XMZ09GjR6NJkyYoKSkxu/2vv/5CgwYN8PDDDytuU5ertGHDhibblixZAkmS8OuvvypuzxqNRmO0SHTt2rV44w3XX0nekdfz66+/xoABA9CwYUM0btwYMTEx2LJli6L+HN130aJF6NChA+rVq6e/mENycjIkSUJZWZmivomIiJTMUV4M4DqA/lWPTwOY47QReQFPCpJzcnLQvn17kys7OdOvv/6KlStXokmTJrjzzjut1r148SJ69+6Nd999F99++y3mzp2L3NxcxMbGWr306Y4dO4xuLVu2xLBhw4zK7Pk5cdKkSSgsLMRXX31ldvuqVatQXFyMSZMmmd1uji5QbtSokeJ91OKuQNne13PBggUYNWoUevfujTVr1iA9PR1jx45FcXGxzb4c3fePP/6AVqtF//79sWXLFmzatMmuY9QtejJ327Ztm11tERGRl1NwtjGn6t+9BmX7HT17qeTm8XmUq+Xe3AqI0Kp/1crJqdRTTz0lMjIyXNKXTnl5uf7+woULBQDx22+/Kd7/0KFDAoCYN2+e4n0iIyNFYmKiPcM0UlJSIpo1ayZGjhxpdrtGoxERERGioqJCcZsdO3YUAMT169dNtulyqh49etThMRuqnjd00qRJIjw8XJW2a8rS6/nbb7+JwMBAMX/+fLvbrMm+siwLAGLz5s1G5dZyxDoqLi5OrFmzRrX2iIjI9WAlj7KSM8rXJUmqj8oFfJAk6WYA5n+/rkaSpDxJkg5IkrRPkqScqrKmkiRlSpJ0tOrfJnbE9R5BY3BfhvuyXpSWlmLLli246667zG7/8ccfMXLkSLRq1QqBgYGIjIy064ypJUqv5mRJs2bNAAB+fkrWkiqzf/9+jBw5Ek2aNEH9+vVx++234/vvv9dvr1evHh544AF88803uHDhgtG+J0+eRFZWFh588EFIkqS4z8uXLyMgIMDkcrQ1tWLFCnTu3BkBAQGIjo42OXP+8MMPY+nSpTh9+rRD01DUZun1XLRoEXx8fPDYY4/Z3aaj+z788MP6KSqDBw+GJEkm02kOHjyIgQMHIigoCK1atcJLL71kdIlcIiIiHSURTzKADQBukiQpDcBmAP+yo4+BQoie4kZ+uucBbBZCdKhq63k72vIoMqxPw1DTTz/9hI8//tio7Ouvv8bQoUPNBpy7du3CnXfeiWbNmuGjjz7C119/jRdeeAGNGzc2qieEQFlZmc1beXm5SR/2KC8vx/Xr13H06FFMnToVLVu2xAMPPFCjNnX27NmD/v374+LFi1i4cCG++OILNGvWDEOGDMHu3TeuizNp0iSUlpZixYoVRvsvW7YMQgg89NBDdvX7119/mZ2fXBObNm3ChAkT0KFDB6xevRr//Oc/8fTTT+Pw4cP6Oi+++CKGDx+OsLAwxdNQ1H6dlbyeP/zwAzp37owVK1bg5ptvhp+fH9q3b4/33nvPZvuO7vviiy/i7bffBgC899572LFjh8kFFUaPHo0hQ4Zg7dq1mDBhAl5++WX85z//UXTcOsnJyWjTpg127NiBKVOmoE2bNvosKEREVItYOtVseAPQDMAIAPcACFWyT9V+edXrAzgMoFXV/VYADlffz9OnXlibbmFUrpL169eLcePGCV9fX3HmzBl9+b333iv27t1rdp9//OMfokOHDjbb3rp1q6JLeFq6XKjSqRe9e/fWt9W+fXvxyy+/2BybIWtTLwYNGiQ6d+4sSkpK9GVlZWWic+fOYtSoUUZ1u3btKvr27WtU1rlzZxEbG2vXeIQQwsfHx+RyqzqOTr3o37+/6NKli9H0lh07dpi8BvZOvajp61ydktezU6dOolGjRiI0NFSkpqaKzZs3i8cee0wAEG+++abV9muyr+6ytNUvJa+bejF37lyj8ilTpoiGDRuKS5cuKTp2IiKqXWBl6oXN374lScoA8BmAdUKIK/bG4QC+lSRJAFgghEgF0EIIcaZq+1kALarvVFBQgJiYGxdI0Wq10Gq1dnbtPDJcu6Bv+PDhGDx4MDZu3Ijly5dj2rRpuHTpEn7//XeLKcmaN2+OX3/9Fc899xwmTZqE7t27m63Xu3dv7Nq1y+YYarpg7dNPP8Wff/6J48ePY968eYiPj8cPP/xQ4ykDV69eRVZWFl544QX4+PgYZTQYMmQI0tLSjOpPmjQJ//73v3HkyBF07NgRO3fuxKFDh/DBBx/Y1W9xcTEqKipUXchXXl6OXbt24fnnnzea3hIbG1vj50nt11nJ61lRUYG//voLS5Ys0V/ae9CgQcjLy8PcuXPxj3/8w+JUl5rsa8u4ceOMHt9///346KOP8PPPP+OOO+5wqE0iIvIuqampSE1N1T0MtVjRUgQtbpwBjgPwPoATAFYBuA9AoK39qvYNr/q3OYD9AAYAKKxW51L1/Tz9jLLNM8lOWsw3efJk0atXLyGEEO+//77VxXDXrl0TL730kmjbtq3+rN8777xjUq+iokKUlpbavJWVlZntx5HFfJcuXRLBwcFi6tSpivexdEb51KlTNs+SGp6dPX36tPDx8RFJSUlCCCGeeOIJERAQIC5evKh4LEIIcfbsWQHA4ploR84o69p89913TbbddtttNTqjXNPX2RpLr2dsbKwAIP7880+j8jfeeEMAEKdPn7bYZk32tXVG+fLly0blBw4cEADEihUrrB0mERHVUqjJYj4hRJYQ4nFUXolvAYBxAM7Z2q9q39NV/54DsAZAXwD5kiS1AoCqfxW15UnclRpu4sSJ2LNnDw4ePIjPPvsMiYmJFusGBARg9uzZOH78OHJzc9GjRw889dRT2L59u1G9rKws+Pv727wNHjxYteMICQlB+/btVckxHBISAh8fHzz11FPYtWuX2Zvh2dnWrVsjPj4ey5Ytw/Xr1/H5558jISEBTZrYt6bUGanhQkND4e/vj/z8fJNt5srs4czX2dLrGR0dbXU/a4tCa7KvLdWfS93j8PBwh9skIqLaSVHagaqsFwkAxgPoBWCpgn0aAPARQvxVdX8ogP8AWIfKy2K/UvXvl44N3X00BvdlWJ6GYfhYDXFxcbjpppvw4osvonHjxmjZsqWi/bp27YpnnnkGX3zxhcnFFlw19cJQfn4+Dh06ZDXQV6pBgwa48847sX//fvTq1UtRADVp0iRMmDABM2bMwPnz5x3KBPLXX38BMH+xEUf5+vqiT58+WLVqFZKTk/XH8uOPPyIvLw+RkZH6ugEBAbh69aritp35Olt6PceMGYOPP/4YGzduxH333acv37BhA9q0aWP1/VuTfW1ZuXIlnn/+xhriFStWoGHDhhanJxERUd2lZI7ySlSeCd4A4F0AWUIIJbmUWgBYUzWP0A/AZ0KIDZIk7QKwUpKkyaiczjHOShseTYb1ucoFKvcnSRImTJiAV1991SRzg6EnnngCV69exZAhQxAeHo5jx44hJSUFffv2xe23325Ut1GjRkbzwZVatWoVAOizSnzzzTcICwtDWFgY4uLi9PXGjBmDXr164ZZbbkHjxo1x5MgRzJ8/H35+fpg+fbrd/ZrzxhtvYMCAARg2bBgmT56MVq1a4fz589izZw/Ky8vxyiuvGNUfPXo0GjdujPnz56N58+Ym6fXy8vLQtm1bzJo1C8nJyWb7VHpGecOGDSZBXXBwMOLj483Wnz17NoYOHYrRo0dj6tSpKCgowKxZs0za6Nq1Ky5evIgPPvgAMTExCAwMtBroOfo6V2fP6zl8+HAMHDgQU6dOxfnz59GuXTukp6fj22+/xeLFi/X1srKyMHjwYCxatEifeUTpvo5YuHAhKioq0KdPH2zcuBEfffQRkpOTERwcXKN2iYioFrI0J0PcmEM8DICvrXpq3jx9jrKrs14YOnDggAgODhZXr161WOftt98W/fv3F82aNROBgYGiU6dO4sUXXzSZ71kTUJg14ZVXXhG9evUSwcHBon79+qJjx45Cq9XaNadZCNsXHPnll1/E+PHjRVhYmKhXr54IDw8XCQkJYv369WbrT5kyRQAQzzzzjMm2n3/+WQAQH3zwgcX+vv76awFAPP7442a36+Yom7tFR0dbPdbPPvtMdOzYUdSrV0907dpVrF692uSCI5cvXxb333+/CAkJEQAsZt9Qm72vZ1FRkXj88cdF8+bNhb+/v+jevbtIS0szqqPLyLF48WK79zXH1hzlAwcOCI1GIwIDA0WLFi3EzJkzjeaxk/tgt3tvNbVmzRpx5513irCwMBEYGCgiIiLEqFGjxDfffKOvo3sfupu943jqqafEiBEjnDgi1/r999/Fk08+KWJjY0X9+vUtrrVJT08X9957r4iIiBCBgYGiY8eO4vnnnzf5PLWUWSg4ONiucW3fvl2MHz9ehIeHC39/f9GoUSMRExMjZs6cKf744w+72ho1apQICQkR165dM7v9zz//FEFBQWLSpElCCCHmz58vunXrVif/HsLKHGVrAfK/DO6Prbbtv5b2U+Pm6YGyoiDZiX8If/jhB6e1TUIsWLBAhIaGiitXrliss3LlSgFA/Otf/3LhyIhqN28OlN966y0BQDzyyCNi3bp1YvPmzSI1NVWMGjVK/POf/9TX+/3338WOHTtq+EzVnD2B8q+//ir8/f3Frl27nDwq19m6dato3ry5uPvuu8XQoUMtBsq33XabGDt2rFi2bJmQZVnMnz9fBAcHi9tuu80ooNQFym+//bbYsWOH/mbPczZv3jwhSZIYNGiQWLx4scjKyhLr168XSUlJonnz5uKuu+6y6xhXr14tAIhVq1aZ3b5o0SIBQGzZskUIIURxcbFo0aKFWLRokV391AaOBsp7zN0391jtm6cHyu7KekGuMWHCBJGSkmJS/tNPP+nPaH744YcCgJgzZ46rh0dUa3lzoHzTTTeJ0aNHm93miWfo7AmUn3zySRETE+PkEbmW4WtiLXvTuXPnTMqWLl0qAIjNmzfry3SBcmZmpkPj2bJli5AkyeyvnEJU/opY/Vc3W0pKSkSzZs3EyJEjzW7XaDQiIiJCVFRU6Mv++c9/iq5du9rVT21gLVC2tvJJsnDf3OM6xV1ZL8g10tLS8MILL5iUr1y5EomJifj++++RmZkJALj55ptdPTwi8kAXL160uMjUcJFxcnKySQ7w5cuXo3Pnzvq1BuvWrYNGo9Ffjt1wv6NHj2LEiBFo2LAhIiMj8Z///MfoEuy//vorHnzwQbRt2xb169dHu3bt8Pe//x2XLl1y6LhKSkqwbNkyTJgwwWRbz5498fDDD2PhwoXo2rUr6tevj/79++PYsWMoKirCU089hRYtWqBJkyZ48skndSfaPILSzDlhYWEmZX369AEAnD59WrXxvPrqqwgNDcWrr75qdnuDBg3w8MMPG5Xt378fI0eORJMmTVC/fn3cfvvt+P777/Xb69WrhwceeADffPMNLly4YLTvyZMnkZWVhQcffNDo/Xj//ffjl19+McmQVZdZe6cIC/fNPa5TNAb3ZVhe0Ee1y4QJE+Dv748BAwbgiy++QIsWLTB8+HB3D4uIPEDfvn2xdOlS/O9//8ORI0cU75eZmYnExER07twZq1evxnPPPYdnnnnGYhtjxozBoEGDsHbtWowePRqzZs3C0qU3ElH98ccfuOmmm/Dmm29i48aNeOmll7B582aH/1ZlZ2ejsLAQd955p1H59evX8csvv2DLli1Yv349Xn31VSxYsAD79+/H3//+dwwZMgRNmjRBWloaHnzwQbz33ntYt26dQ2MwJIRAWVmZzVt5eXmN+7IkKysLANClSxeTbYmJifD19UWzZs0wYcIEnDx50mZ7ZWVlyMrKQnx8POrVq6doDHv27EH//v1x8eJFLFy4EF988QWaNWuGIUOG6BfZA5VZnkpLS00SACxbtgxCCP0Cap2ePXuiUaNG2LBhg6Jx1AmWTjUDKAfwJ4C/AJRV3dc9LrW0nxo3T596YXW6hUE51T5r1qwR/fv3F+PHjxe//vqru4dDVKt489SLw4cPi+7du+sXcTVr1kzcf//9YuPGjUb1qk956Nevn4iOjjb6+TsnJ8dkcbRuv+rzR7t16ybi4+Mtjqu0tFR8//33AoDYs2ePxXFY8sorrwhJkkRJSYlR+e7duwUAMW7cOKPy++67TwAQ6enp+rKysjLh5+cn/vvf/9rszxZLi+aq36ovLLfGngtnnTp1SoSFhYkhQ4YYle/Zs0dMnz5drFu3Tj+XOSwsTLRu3Vrk5+dbbVN3sannn3/eZFv1i0LpDBo0SHTu3NnodSkrKxOdO3cWo0aNMmqja9euom/fvkZlnTt3tnjBrDvuuMPqe6o2giOXsBZC+KobktcuMmxfxppqn9GjR2P06NHuHgYReZiOHTti79692LZtG7799ltkZ2djzZo1WLFiBV5++WXMnDnTZJ/y8nLk5ORgxowZRj9/9+7dG23btjXbz4gRI4wed+vWDXv37tU/vn79OubNm4dPPvkEJ06cwLVr1/TbDh8+jFtvvdWu4/rjjz/QuHFjkzOduj7/85//GJVfuXIFt9xyi1H+86tXr6KsrAzNmjXTl+3ZswfDhw/H2bNn9WWFhYVo1qwZLl68aDFdozty/+tcvnwZo0aNgp+fn0mayltvvdXouY2Li8OAAQPQt29fvP3225gzZ47d/Z09exatWrUyKistLUVpaSmysrLwwgsvwMfHx+j6CEOGDEFaWprRPpMmTcK///1vHDlyBB07dsTOnTtx6NAhfPDBB2b7DQsLs+tXkdpO0QVHyJgM20GyBkREVJf4+vpiwIABGDBgAIDKIPOuu+7C7Nmz8cQTT5hcAfT8+fMoLS1F8+bNTdpq0aKF2T6aNm1q9DggIMAoGJ4xYwbeeecdvPTSS+jfvz8aNWqEU6dO4d577zWqp9S1a9cQEBBgUr53715ERESgU6dOJuUPPvigUdn+/fsBAD169NCX7d69G7fddptRvZ07d6Jdu3ZWc5o3bNgQPXv2tDnu6vPAa+rq1atISEjA8ePHkZWVhTZt2tjcp1evXujYsaPNwL5Zs2YIDAw0maYRGhqq3zc1NRULFy4EUDkfvry8HC+//DJefvlls21WVFTo52FPnDgRM2bMwCeffII5c+bgk08+QUBAAMaPH2923/r169t1QavazvHrwNZhDJKJiMiW1q1bY8qUKSgrK8PRo0dNtusuW3/u3DmTbY5etn7FihV46KGHMHPmTAwaNAh9+vRBSEiIQ20BlUFcYWGhSfnevXvRq1cvo7KzZ8/i7NmzJuV79+6Fr68vbrnlFn1ZTk6OSaD8448/onfv3lbHk5WVBX9/f5u3wYMH23mklpWWluK+++5DTk4Ovv76a7uv4mkraPfz88OAAQOQmZmJ69evG5XHxMQgJiYGrVu31peHhITAx8cHTz31FHbt2mX2ZrhYsXXr1oiPj8eyZctw/fp1fP7550hISDD54qZz8eJFhIaG2nWMtZlHnlEuKiqCVqtFQkICEhIS3D0cEwySiYjI0JkzZ0x+JgeAQ4cOAYDZjBi+vr6IiYnBF198YZQNY/fu3fjtt98QERFh9ziKi4vh7+9vVFaTq1l27twZ169fx6lTp/RnUSsqKrB//378+9//Nqqrm45RfXrH3r170blzZ9SvX19ftnv3bnz55ZdITU3Vl507dw6zZ8+2Oh5XT72oqKhAYmIitmzZgq+++gqxsbGK983JycHhw4eNpqFY8q9//Qvx8fH497//jfnz51ut26BBA9x5553Yv38/evXqpSiDx6RJkzBhwgTMmDED58+fx6RJkyzW/e2339C3b1+bbdYGGRkZyMjIAACLP2N4ZKAcHBxs9J/H02gM7suwfIbZ8DEREdVe3bp1w5AhQzB8+HC0bdsWf/75J77++mt8+OGHGDdunMWgV3fZ+jFjxkCr1eL8+fNITk5Gy5YtFacwM3TXXXdh6dKl6N69O9q3b4/Vq1fXKNWXbhrJzp079YHy0aNHceXKFbNnjhs2bIgOHTqYlBvWvX79Og4cOIDTp08bnbls2bKlzTPKjRo1QkxMjMPHY2jVqlUAoM8S8c033yAsLAxhYWGIi4sDADzxxBNIT09HUlISGjRogOzsbP3+bdq00T8niYmJaNu2LXr16oWQkBDs3bsXc+fORXh4OP7xj3/YHMvgwYPxyiuv4Pnnn8dPP/2Ehx56CG3btsW1a9dw5MgRrFixAg0aNNB/mXrjjTcwYMAADBs2DJMnT0arVq1w/vx57NmzB+Xl5XjllVeM2h89ejQaN26M+fPno3nz5rjrrrvMjqOwsBBHjhzBc889Z+ez6Z10J2QXLlxYZLGSpVV+7rwx6wUREXmTDz74QCQkJIiIiAgREBAggoKCRM+ePcWrr75qlJnAXLaJtLQ0k8vW9+zZ0+gCJrr9DDMfCCHEpEmTjC5hX1BQIMaPHy9CQkJESEiImDBhgti5c6fJZeLtueBI3759xcMPP6x/vHz5cgHA5JLK9913n7j99tuNyq5fvy7q1asn5s+fry/LyckR4eHhRvXOnDkjAIhLly4pGpMaoCBjRmRkpMV6s2bN0tf773//K7p37y4aN24s/Pz8RJs2bcSjjz5q92Wnf/jhBzF27FjRunVro0tYv/TSSyZt/fLLL2L8+PEiLCxM1KtXT4SHh4uEhASxfv16s21PmTJFALB4URMhhFi2bJkICAgQ58+ft2vc3g5Wsl5Ilds9S0xMjMjJyXH3MCyTJGUL+jzwuSUiIs926tQptG/fHklJSXjxxRfdPRwsWbIETz/9NM6cOYOgoKAat7dgwQJkZGTgq6++0pdt2LABTz75JH799dcat0+Ou/vuuxEaGopPP/3U3UNxKUmSdgshzP5UwcV8DpDBBX1q0n1ZS05ORl5enkddvYmIyJmuXr2Kv//97/jiiy+QlZWFxYsXIz4+HkFBQZgyZYq7hwegMmtC69at8f7776vS3u7du00yV+zbt8/mtAtyrn379mHLli2YNWuWu4fiUTxyjrKnY5CsrpUrV2Lv3r24du0adu7ciaSkJLzxxhsW0yMREdUWvr6+OHv2LJ588klcuHBBv1ArPT3d7OJAd9DlDd6zZ48q7Zlbg/T888+r0jY57uzZs1iyZAnat2/v7qF4FktzMtx58/Q5ypbmJFcvV5vu6kE5OTkm2woKCgQAMXPmTNX7vXbtmvD397d5JaRffvnFrrqGfvzxR9GoUSMRFxcniouLVT8GW4YNGyYAiKSkJLPbZ8+eLbp16yY+//xzF4+MiIiInAmOXJmPLNMY3JfhuqwXe/fuhZ+fH7p162aybd++fQBM0/Ko5bvvvjNbvnz5crz99tsYMGAAOnbsiLKyMsV1ddLT07F792488sgj6N+/P7RaLebNm+eyM8rLly/XJ8Q35+uvv0ZqairuuecefPvttxg3bpxLxkVERETuxUC5BmRYn4ZRoHJ/+/btQ5cuXcxeJUkXKCu5YpG9AgICzOaO/PLLL/H++++jT58++Oqrr+Dr6wtfX1/FdXXuu+8+jB07FsnJyejbty/Gjh2r+lWVLLl06RKmTZuG+fPnY8KECWbrLFy4EI8++iiOHDmCm2++2SXjIiIiIvfjYj4HybA9V1lNQgj89NNPFs8Y79u3D8HBwWjbtq1d7ZaVlWH27Nm4fPmyXftlZmZi/Pjx6NKlCzZs2GA1ubuturqgODk5GVFRUTUOku05pn//+9/o1q0bHnjgAbPbS0pK8O2332LkyJHYvn07+vfvX6OxERERkfdgoOwAGa5f0Hf06FFcvnwZHTt2RGFhoclt37596Nmzp91B5unTp/Hhhx9i2LBh+PPPPxXt88MPP2D06NGIiIhAZmYmmjZtqkpdtSg9ph9++AGffPIJ3nvvPYt1srOz4efnh3r16uHixYv6JPRERERU+zFQdoA7sl7oplbMnDkTTZo0Mbnl5uY6NO0iMjISW7duxfHjxzF06FAUFVm+OA1QmdZnxIgRCA0NxebNm63OI7anrpqUHNP169cxdepUPPfcc+jUqZPFtrZt24ZevXph2bJleOCBB4wuwUpERES1G+coO8AdqeF0gfKGDRtM5ij/9NNPePrpp81Oyzh//jzCwsIU9XH27Fm88MILFs+w5ubmYtiwYQgKCsLmzZtx0003WWzLnrr2UuOYXnvtNVy9ehVJSUlW98/NzUXbtm2xaNEii4sUiYiIqHZioOwAjcF9Ga7JerF37160bdsWw4YNM9l26NAhAOYzXjRp0gQHDx602vb58+cxduxYBAQEWLy++7FjxxAfHw8A2LRpk9U8i/bUdURNj+nkyZNISUnBRx99hJKSEpSUlOi3lZSUoLCwEI0aNYKvry/++OMP5OXlYejQoVbPPBORe2k0GgCALMtuHQcR1S4MlGtAhuuyXuzbtw/9+vWzuC0gIABdunQx2ebr64vOnTtbbPfChQsYP348AgMDsXXrVkRFRZnU+f333zF48GBcuXIFmzdvRnR0tMX27KnrqJoe0/Hjx3Ht2jVMnDjRZN958+Zh3rx52Lt3L3r27InCwkIUFBRg7ty5ah8GEREReTiPDJSLioqg1WqRkJCAhIQEdw/HLBmuy3qRn5+Ps2fPWs14ER0dDX9/f7vb9vX1RceOHfH6668jIiLCZPu5c+cwZMgQFBQUYMOGDYiJMXspdLvrOpOtY+rZsye2bt1qUj5w4EBMnDgRkydP1p8Fr6iowNNPP402bdo4fdxERETkOhkZGcjIyACAYEt1PDJQDg4ONnuJS08hw7UL+vbu3QvA/NSKiooKHDhwwGJ6M1tCQkKQnm4+rL9y5Qri4+Nx5MgRzJ49G/7+/sjOzjap165dOzRo0EBx3ebNmzs0VqWsHZNuu+5n2uoiIyP125YuXYoDBw6gW7duKC8vx7Rp0/DUU0+hQ4cOThg1EXmrJUuW4G9/+5v+cVBQEMLCwnDrrbfigQcecCg3/L59+7B27Vr84x//cEm2IKK6SHdCduHChRYzGXhkoOzpXJ31wtrFRI4cOYLi4mKnXGhkx44d+OmnnwAAs2bNwqxZs8zW27hxI3x8fBTXHTp0qOpjVVtxcTHS09Px1Vdf4cUXX0S3bt3w6KOPMkgmIovS09PRpk0blJSU4OTJk1i/fj0eeOABpKamIiMjw66sOfv27cPs2bMxceJEBspEbsRA2QGuznrx/PPP4/nnnze7rXPnzqi8TLn6hgwZYlfbzhqHqxiOPygoCF999RUAYPjw4e4aEhEpkJaWhuzsbJSUlCAqKgopKSlITEx0+Th69uxptHj5wQcfxNixYzF27Fj861//wjvvvOPyMRFRzTCPsgM0BvdlWD7DTEREzpWWlgatVqvPXnPixAlotVqkpaW5eWSV/u///g+jRo3CwoULUVxcDKDyV7devXqhcePGCA0NxaBBg4ymqhlO5ejQoQMkSYIkScjLywMAvPvuu+jXrx+aNm2KkJAQxMbGYv369S4/NqK6gIFyDciwPg2DiIicKykpSR+A6hQXF9vMke5Kw4cPR0lJCXJycgBUXj102rRp+PLLL7FkyRI0b94cAwYMwIEDBwAAI0aMwMyZMwFUTufYsWMHduzYgVatWgEA8vLyMGXKFKSnp+Pzzz9HTEwM7rnnHmzYsME9B0hUi3HqhYNkuC7rBRERmXfy5Em7yt1Bl33nzJkzAICPPvpIv628vBx33XUXoqOj8dFHH+Gtt95CWFgYbr75ZgCm0zmAyjSWOhUVFRg8eDCOHDmCDz74AHfddZezD4eoTuEZZQfIcM9lrImIyJi5FJDWyt1Bt/5Bl/li06ZNGDhwIJo1awY/Pz/4+/vjyJEjOHz4sKL2du/ejXvuuQctWrTQ75+Zmal4fyJSjoGyAxgkExF5hpSUFAQFBRmVBQUFISUlxU0jMvX7778DAFq1aoU9e/Zg+PDhaNiwIT7++GNkZ2dj165d6NGjB65du6aorcGDB+PixYt45513sH37duzatQt33XWXov2JyD6ceuEABslERJ5Bl91i8uTJKCkpQWRkpNuyXliyfv16BAYGonfv3khJSYGfnx9Wr15tdJGoS5cuISQkxGZbGzZsQFFREVauXGl0IaTq87SJSB08o+wAjcF9Gcx6QUTkTomJiYiNjUVcXBzy8vI8Kkj+4osvsG7dOjz22GMICgpCcXExfH19jS5AsmXLFpM51QEBAQCAq1evGpXrAmLDIPvIkSPYtm2bsw6BqE5joFwDMpj1goiIKu3btw/Z2dn47rvvsGzZMtx///0YN24c4uPjMXfuXADAXXfdhcuXL+Phhx/G5s2b8cEHH2DixIkIDw83aqtr164AgPfeew87duxATk4Orl+/jiFDhsDPzw8PPfQQvv32WyxduhRDhw71qDnZRLUJA2UHyWDWCyIiumHs2LHo168fhg0bhqSkJJSUlGDFihXYsGEDAgMDAQDDhg3D22+/jW3btuGee+7BokWL8Mknn5hktujRoweSk5ORkZGBO+64A3369MEff/yB6OhopKWl4cSJExg5ciRee+01vPLKKxgwYIA7Dpmo1pM88WpqMTExQpdv0hPJkqRsQZ8HPrdERLWRRqMBAMiy7NZxEJH3kSRptxAixtw2j1zMV1RUBK1Wi4SEBCQkJLh7OCaY9YKIyLMwQCYie2VkZCAjIwMAgi3V4RllB8iSpCxI9sDnloiIiIhusHZGmXOUHaAxuC+DWS+IiIiIaiMGyjUgg1kviIiIiGorBsoOksGsF0RERES1GQNlB8jggj4iIiKi2o6BsgMYJBMRERHVfgyUHcAgmYiIiKj2Y6DsAI3BfRnMekFERERUGzFQrgEZzHpBREREVFsxUHaQDGa9ICJSkyzLCAsLM7nKnqvKa2Lt2rUYMGAAmjdvjvr16yMyMhKjR4/Ghg0bjPpJTk5GRUWFQ33s27cPycnJuHjxosPjjIqKwsSJE81uS05OhiRJKCsrc7h9otqGgbIDZHBBHxGR2saOHYv09HRoNBp9mSzLLit31Ntvv40xY8agQ4cO+Pjjj7F+/XrMnDkTALBlyxajvmbPnl2jQHn27Nk1CpSJyD5+7h6AN2KQTESkPncGyYbl9po3bx5Gjx6Njz/+WF82aNAgPProow4HxUTkGXhG2QEMkomI1OeNQTIAXLx4ES1btjS7zcen8mM2OTkZs2fPBgD4+/tDkiRIkqSvN2vWLPTq1QuNGzdGaGgoBg0ahOzsbP32JUuW4G9/+xsAoEOHDvr98/LyajR2a/Ly8iBJEhYsWICXXnoJrVq1QkhICBISEnDq1CmT+gsXLkSvXr1Qv359NGnSBHFxcdi+fbvTxkfkCgyUHaAxuC+DWS+IiNTkTUEyAPTt2xdLly7F//73Pxw5csRsnSlTpmDy5MkAgB9++AE7duzAjh079NtPnz6NadOm4csvv8SSJUvQvHlzDBgwAAcOHAAAjBgxQj+dIz09Xb9/q1atajx+W+bOnYtff/0VixYtwltvvYUdO3aYzHN+7rnnoNVq0atXL6xcuRLLli3DgAEDcPLkSaePj8iphBAed2vfvr149NFHxbp164RHAoQAxFZAhFb9KwxuunIiIrLP1q1bRWhoqNi6davLyx11+PBh0b17dwFAABDNmjUT999/v9i4caNRvVmzZgkAorS01Gp7ZWVlorS0VHTs2FH84x//0JcvXrxYABBHjx51eKyRkZEiMTHR7Lbq4/vtt98EABEXF2dU73//+58AIE6fPi2EEOLo0aPCx8dHTJs2zeFxEbnDunXrxKOPPioAHBUWYlKPPKMcHByM1NRUJCQkuLRfu1ZJg1kviIjU5K4zybpyR3Xs2BF79+5FVlYWkpKS0LNnT6xZswbDhg3DnDlzFLWxadMmDBw4EM2aNYOfnx/8/f1x5MgRHD582OFxqWX48OFGj7t37w4A+rPFmzZtQkVFBbRarcvHRlQTCQkJSE1NBYAiS3U8MlB2F8V/VMEFfUREavPWrBcA4OvriwEDBmDOnDnYtGkTjh8/ju7du2P27Nm4dOmS1X337NmD4cOHo2HDhvj444+RnZ2NXbt2oUePHrh27VqNxlWdn58fysvLzW4rLy+HJEnw9fU1Km/atKnR44CAAADQj+3ChQsAgDZt2qg6ViJPwEDZgOI/qmCQTESkNm/NemFO69atMWXKFJSVleHo0aNW637xxRfw8/PD6tWrMXr0aNx2222IiYmxGWA7onnz5vjjjz/Mbvvjjz8QFhZmtMhQidDQUACV86yJahsGymbY/KMKBslERGrz1iD5zJkzZssPHToEAPqMGLozsVevXjWqV1xcDF9fX6MAdcuWLSYL4Sztb4+BAwciOzvbJFi+evUqvvnmGwwcONDuNocMGQIfHx/dT9hEtYulycvuvPXu3dtJ07ZtU7QAxMzCPXML+oiIyH7uWtDnqKZNm4px48aJJUuWiKysLJGRkSH+/ve/C0mSxLhx4/T11q5dKwCIWbNmiezsbLFr1y4hhBAbNmwQAERiYqLYtGmTeP/990WrVq1EeHi40UK6ffv2CQBi6tSpYvv27WLXrl2ipKRECCHEoEGDxM0332xzrPn5+aJ169YiMjJSLFiwQGzZskV8+umnomfPnqJBgwbi559/1tfVLeZbuHChURtbt24VAIyev+nTpwtJksSjjz4qMjIyxNdffy2Sk5PFihUr9HVuvvlmMWjQILueWyJXAJAjLMSkbg+Kzd3cFSgr/qPKrBdERFTlgw8+EAkJCSIiIkIEBASIoKAg0bNnT/Hqq6/qA1khKrNZPP744yIsLExIkiRg8Dnx9ttvi6ioKBEYGChiYmJEZmamiIuLM8k4kZycLFq3bi18fHwEAPHbb78JIYSIi4sTkZGRisZ74sQJ8fDDD4tWrVoJPz8/0axZM3HvvfeKn376yaiePYGy7nno3r27qFevnmjSpImIi4sT27dv12+PjIw0OR4iT2AtUJYqt3uWmJgYkZOT4/J+w8LClP08p3T+lgc+t0RERER0gyRJu4UQMea2cY6yAVcu9CAiIiIiz8ZA2YDSBSBEREREVPsxUDbDWUnpiYiIiMh7OD1QliTJV5KkvZIkfVX1uK0kST9KkvSrJEmfS5JUz9ljsIezk9ITERERkXdwxRnlpwEcNHj8KoD5Qoj2AC4BmOyCMSjijqT0REREROSZnBooS5LUBsAIAB9VPZYADAKwqqrKUgCjnTkGezBIJiIiIiIdPye3/yaAfwFoVPW4GYBCIURZ1eNTAMKr71RQUICYmBtZOrRaLbRarXNHCma9ICIiIqoLUlNTDa8mGWqpntMCZUmS7gFwTgixW5IkjT37hoWFwR15lJVmvdCY7ElERERE3sLwJKwkSect1XPmGeXbAYyUJGk4gEAAjQG8BSBEkiS/qrPKbQCcduIYHGJrrnKB+4ZGRERERC7itDnKQogZQog2QogoAPcD2CKESASwFcB9VdUmAfjSWWNwBLNeEBERERHgnjzK/wbwrCRJv6JyzvLHbhiDWcx6QUREREQ6zl7MBwAQQsgA5Kr7xwH0dUW/9mKQTEREREQ6vDKfAQbJRERERKTDQNmA0qwXRERERFT7MVA2w9ZcZSIiIiKq/RgoV8OsF0REREQEMFA2wqwXRERERKTDQNkAg2QiIiIi0mGgbIBBMhERERHpMFA2wKwXRERERKTDQNkMZr0gIiIiIo8MlIuKiqDVapGRkeHyvpn1goiIiKj2y8jIgFarBYBgS3UkIYTrRqRQTEyMyMnJcXm/irNeSJKyBj3wuSUiIiKiGyRJ2i2EiDG3zSPPKLsLs14QERERkQ4DZQMMkomIiIhIh4GyAWa9ICIiIiIdBspmMOsFERERETFQroZZL4iIiIgIYKBsRHHWCyIiIiKq9RgoG2CQTEREREQ6DJQNMEgmIiIiIh0GygaY9YKIiIiIdBgom8GsF0RERETEQLkaZr0gIiIiIoCBshFmvSAiIiIiHY8MlIuKiqDVapGRkeHSfhkkExEREdUNGRkZ0Gq1ABBsqY4khHDdiBSKiYkROTk5Lu9XlmVlQbIkKWvQA59bIiIiIrpBkqTdQogYc9s88oyyuzDrBRERERHpMFA2g1kviIiIiIiBcjXOynqh0Wg4x5mIiIjIizBQNsCsF0RERESkw0DZAINkIiIiItJhoGyAQTIRERER6TBQNsCsF0RERESkw0DZDGa9ICIiIiIGytU4K+sFEREREXkXBsoGmPWCiIiIiHQYKBtgkExEREREOgyUDTBIJiIiIiIdBsoGmPWCiIiIiHQ8MlAuKiqCVqtFRkaGW/pn1gsiIiKi2i0jIwNarRYAgi3VkYQQrhuRQjExMSInJ8ctfSta0DdwoLLGDJ5bXVs8I01ERETkOSRJ2i2EiDG3zSPPKLsLs14QERERkY6fuwfgScYOHIh0wOiMsQxgLGBSTkRERES1G88oG0gHoDF4LMMgSK5Bu7IsY9u2bSgsLDQpJyIiIiLPxEDZgMbgvgzzQbJsZ5u6aRvR0dEICQkxKSciIiIiz8RA2QwZloNke0Jbw7nN5oJkXg6biIiIyHMxUK5GhvUg2Z7QlgsDiYiIiLwXA2UDMmwHyRooVz0YLiwsZJBMRERE5CWY9cKAmkEyAJMgOTc3F5mZmQySiYiIiLwAzygbcGbWi9zcXERHR/Ny2ERERERegoGyAY3BfRnMekFERERUlzFQNkMGs14QERER1XUMlKuRYX/WCw3MT81g1gsiIqotNBoNP7dIFd70XvLIQLmoqAharRYZGRku7VcGs14QERER1QUZGRnQarUAEGypjkdmvQgODkZqaqrL+2XWCyIiIqK6ISEhAQkJCVi4cGGRpToeeUbZXZj1goiIiIh0GCgb0Bjcl8GsF0RERER1GQNlM2Qw6wURERFRXcdAuRoZ9me9sIRZL4iIiIi8FwNlAzKY9YKIiIiIKnlk1gt3YdYLIiKyh+7vORdnE9VOPKNsgFkviIiIiEiHgbIBjcF9Gcx6QURERFSXMVA2QwazXhARERHVdQyUq5HBrBdERERE5MTFfJIkBQL4DkBAVT+rhBCzJElqC2AFgGYAdgN4UAhx3VnjsIcM52S9SEtLQ3Z2NkpKSjB48GDMmDGDQTIRERGRh3PmGeUSAIOEED0A9ARwlyRJsQBeBTBfCNEewCUAk504Brs4I+tFWloatFotSkpKAAAVFRWYP38+0tLSajxeIiIiInIepwXKotLlqof+VTcBYBCAVVXlSwGMdtYY7OWMrBdJSUkoLi42KisuLkZSUhKzXhARERF5MKfOUZYkyVeSpH0AzgHIBHAMQKEQoqyqyikA4dX3KygoQExMjP6WmprqzGHqaQzuy1An68XJkyctljPrBRERuZtGo+F0QKpzUlNT9XEmgFBL9Zx6wREhRDmAnpIkhQBYA6Czkv3CwsKQk5PjzKFZJcP6NIwCpe3IMiRJghDCZJskScx6QUREROQGWq0WWq0WACBJ0nlL9VyS9UIIUQhgK4B+AEIkSdIF6G0AnHbFGJSSoW7WixkzZiAoKMhkGxf0EREREXk2pwXKkiSFVZ1JhiRJ9QHEAziIyoD5vqpqkwB86awx2EuG+lkv5syZg9TUVAQEBAAAfHx8kJSUhDlz5qgwYiIiIiJyFmeeUW4FYKskST8B2AUgUwjxFYB/A3hWkqRfUZki7mMnjsEuzsh6AQCJiYno3Lkz/Pz8sHnzZgbJRERERF7AaXOUhRA/AbjVTPlxAH2d1W9NOCPrBVA5Vzk3NxfR0dEmFx/h9AsiIiIiz8Qr8xnQGNyXoU7WC92V+KKjo81expqIiIiIPBMDZTNkWJ+Gobgdg8tVmwuSmfWCiIiIyHMxUK5GhrpZL3SXsda3YxAkc9oFERERkedioGxAhvpZLwyD4cLCQgbJRERERF7CqRcc8TbOynoBVAbJubm5yMzMZJBMRERE5AV4RtmAO7JeEBEREZFn4hllAxqD+zIsn2E2fGyVJOnbiQYQsn8/IElG7ReYubw1EREREbkfzyibIUOlrBcG7YRYKCciIiIiz+SRgXJRURG0Wi0yMjJc3rcMFbNe2GhHY7IHEREREblCRkYGtFotAARbquORUy+Cg4ORmprq8n5lqJz1olr9QgfbISIiIiJ1JSQkICEhAQsXLiyyVMcjA2V3URokP/nAAaP9jm7+W2X54MVG5e8u766/XwggF0AmGCQTEREReQOPnHrhLk7LeoHKIDnaTPtERERE5JkYKBvQGNyXYfkMsz107UTD/II+IiIiIvJMDJTNkKE868WRvPXIP/8T/jiXg6VfDsWRvPVm2wmx0D4REZGnS0tLQ3Z2NrKyshAVFYW0tDR3D4m8lLe9lxgoVyNDedaLI3nrsXVnMsorrgMALhefwdadyfpgmVkviIjI26WlpUGr1aKkpAQAcOLECWi1Wo8PcMjzeON7iYGyARn2Bbc79r+FsvJrRm2UlV/Djv1vAWbqF1poh4iIvIssywgLC0NhYaHZcnvbqX6lVkvl7pCUlITi4mKjsuLiYiQlJblpROStvPG9xEDZgL1ngC8XnzXbjq7csH4hKhf0MUgmIvJusixj7NixSE9PR0hIiNlye9vRaDQ2y93l5MmTdpUTWeKN7yUGygbsnSbRMKil2Xaql8tg1gsiotpAzeDWG4JkAIiIiLCrnMgSb3wvMVA2oDG4LwOIBxBuplynX4+n4ecbaNSGn28g+vV42qg+s14QEdUO5oLYwsJCh4JbZwfJGo1GlXZSUlIQFBRkVBYUFISUlJQat011ize+lxgomyFDWXDbMWoEBvZNhq9PPQBAw6BWGNg3GR2jRhjVZ9YLIqLawVyQnJub61BwqyRI9oQ5yomJiUhNTUVAQAAAIDIyEqmpqUhMTHTzyMjbeON7iYFyNTLsC247Ro1Ai9Bb0Lp5DCaN+lYfJAPMekFEVNsYBrEzZ87E/v37UVZWhocffli/ct/e4NbWdA5PkJiYiNjYWMTFxSEvL8+jAxvybN72XuIlrA3IUBbcrlLYHrNeEBHVTjNnzjT6uViX5urgwYNYsGABCgoKFLWjZM4zEbkPA2UD5oLYQgvlSlRvJxdApgPtEBGR55BlGXPnzjUpLy4uxty5c7F582bF7XjLgj6iusojp14UFRVBq9UiIyPDpf2aC5LVSOkmg1kviIhqA10QK4Qwu10IUeuyXhDVVhkZGdBqtQAQbKmORwbKwcHBSE1NRUJCgkv71Rjcl6FOcCuDWS+IiGoLXRCrRporBslE7pWQkIDU1FQAKLJUxyMDZXeToU5wq6vPrBdERLWDLohVI82Vt2S9IKrLGChXI0O94JZZL4iIahddEJuYmIhp06bpyw3TXNXGrBdqUyvHM5GzMVA2IEPd4JZZL4iIaidZlrFgwQL06NHDKM2VvcGtu7NepKWlITs7G1lZWYiKitKnuCOiSgyUDbgi6wWDZCIi72YYxIaEhJgtt7cdd8xVTktLg1arRUlJCYAbKe4YLBPdwEDZQAFMg9gQC+X2kMGsF0REtYGawa27s14kJSWhuLjYqKy4uBhJSUlO7ZfImzBQdjIZzHpBRFRbmAtiCwsLHQpu3Z314uTJk3aVE9VFDJSdSAazXhAR1SbmguTc3FyHglt3Z71QI8UdUW3HQNmJmPWCiKh2qR7E5ubmIjo6ukbBrbuyXqiR4o6otmOg7ETMekFEVDvpgtjo6GizC/rsbccdWS8SExORmpqKgIAAAMYp7oioEgNlJ9IY3E8GcBoMkomIvJ0jWS/M5Q12d9YLoDJYjo2NNUpx50nUzrfsrvzN7uiXuarV4efuAdQFMixPw9CY1CYiIk8lS9KNv+cDB1ouF8JmW2MHDjRtBzBtX0FbROQcPKPsZDKsz1UmIiLvoWa+fa5hIfJ8HhkoFxUVQavVIiMjw91DqREZthf0ERGR9zAXJDt6MSnD+jIsf14QkXNkZGRAq9UCQLClOh4ZKAcHByM1NRUJCQnuHkqNMOsFEVHtojG4L0Odi0nJ4C+PRO6QkJCA1NRUACiyVMcjA+XagkEyEVHtJEOdi0np6nv6L4+yLGPbtm0oLCw0KSeqzRgoO5HG4L4MBslERLWBDHUuJmVYX6Og3F3USoVH5I0YKLuADM49IyKqDWSoF9za044sywgLCzM5g2upXC2OpMIjqk2YHs7JZFj/Y1jg8hEREdVduryytgJLS/WcmfWiJyrnPGfCNEg2zKusG1tycrJT8y07mudZ6XNM3qWuvq48o+xkGlQGwxoL5URE5D2cmfXC0sJAc0FpYWGh0y9K4q5+a4PqwWRhYSG2bdtW54LM2oCBMhERkUIag/sy1M16YWlhoLlgNTc31+nBqtJ+3TUtxJNVf35yc3MRHR3NLxdeiIEyERGRnWSon/XCXDvpUBZ0OSMYVdovzzxbZm0hJHkHBspERER2kOGerBfuyj5hq193nfH2dLYWQpJ3YKBMRESkkAznZr0otNS+m7JPKOnXXWe8PZmShZDkHRgoExERKeTMrBeFML8w0NHsEzVlb7/Mt3wDp6PUHgyULUgDkA0gC0BU1WMiIqrbvCXrRVpaGrKzs5GVlYWoqCikpdn/KWZPv8y3bIzTUWoP5lE2Iw2AFkBJ1eMTVY8BINEtIyIiIk+gMbgvw3Jwa/jYFhmOZb3IzMy0GCRrtVqUlFR+ip04cQJabeWnWGKi8k8xpf2664y3IU/L8etNWS887bnzNDyjbEYSgOJqZcVV5URERDI8L+uFLhVbUlISiouNP8WKi4uRlGTfpxizXtQcs154P48MlIuKiqDVapGRkeGW/k/aWU5ERHWHDM/MeqELSk+eNP9pZanc5jhVynpR1/ItM+uF58vIyND92hJsqY5HBsrBwcFITU1FQkKCW/qPsLN8zODFGDN4sbOGQ0REHkKG52a90AWlERHmP60iIiLsDkYN29+3b59+f0eyXrh7eoYrMeuFdbIse8QXo4SEBKSmpgJAkaU6Hhkou1sKgKBqZUFV5UREVHdpABTANBi2VG6Ns7JeTJw40aSvoKAgTJw40a4zmWpnvahL0zPq2vHWZgyUzUgEkAogoOpxZNVjLuQjIiK1aKo9DgFwu5lye7NPLFiwAElJSQgIqPwUi4yMxLRp07BgwQK7zmQ6EiQrzbdc27NA1LXjrc0YKFuQCCAWQByAPDBIJiIi5czNuS0sLHRoLq49c4B1QemcOXMQGxuLuLg4LFmyRB8k2xOkORIkKymfOXMm9u/fj7KyMjz88MP61HWe8FO8EkrmWntT1guyjunhiIiIVGZvSjdr3JV9QjNw4I32YTB3unp5aKjifmfOnImUlBsTGXWp6w4ePIgFCxagoKDArjG6mprTURg0eweeUSYiIlKZM84oqpV9wu5+YWMBox1nvOfOnWvSfnFxMebOnetRC9w0Go3Z50zN6SjkHXhGmYiIyEmsnlFU2ogkGQWryefPmy1XeubZnqDZqH1L5Xac8RZCmO1HCOEVZ1jVnI5C3oFnlImIiBwkV93MblPpjKIMdfMtu6vf9PR0q6nrHCXLMrZt24bCwkKTcrXVJEhm1gvvxECZiIhIZWqeUVQ733JN+pVr0K9Go0FKSgqCgowTsAYFBRnNW7aHWl8KHO3X4ut7/nzlXG5JAiQJhVlZyN2/37jcwX4tLSQk52CgTEREpDI1zyg6K9+yvf3KsBAkSxLGDhx4IwjMygKyskzLASQmJmLatGn6fSMjI5GamorExES7gz01vxQ42q/F59mwPipfr2jYl2fb3n7JOZwWKEuSdJMkSVslSfpFkqRcSZKeripvKklSpiRJR6v+beKsMRAREbmDYTCTlpaG7du3Y//+/fD19cXp06ftaktjcF+G+aBLhvoXuTBpH8qD50JL9avyPPfo0QNxcXHIy8vTB8nVg720tDRkZ2cjKysLUVFR+jRyunbs/VJgrT2l7J6TDCABgARgP4AoAGmwPF2nJv0aUuNYqZIzzyiXAZguhOiKypTET0iS1BXA8wA2CyE6ANhc9ZiIiKjWMAySJ0+ejNLSUgBAfn4+tFqtY0EaKoPPaFRenKR6uduyXsD+M962zgCnpaVBq9WipKQEwI00crrnzd4vBbbaU8qRILm06gYAJwBMripXyt7gXK1jpUpOC5SFEGeEEHuq7v8F4CCAcACjACytqrYUwGhnjYGIiMidpk+frg9YdIqLizF9+nS72pFxIygNsVCuNPuEo/1qFJbbc8bbUrCXlJSE4uJio7EUFxcjKSkJgP1fCmy1p5Q9QfJYAA0AGL/6lY8b2NGnvWfO1TpWquSSOcqSJEUBuBXAjwBaCCHOVG06C6BF9foFBQWIiYnR31JTU10xTCIiItXIsoz8/Hyz2yyVm20Hzs16Icuy2QDa7n5h/xlvS8HeyZMnTcZjWG7vlwJb7SlVvf34+HiEh4ebBrGofH7OWWjHXLml3M32zjlX61gdZek4PE1qaqo+zgQQaqme0wNlSZIaAvgCwDNCiD8Nt4nKhIomSRXDwsKQk5Ojv2m1WmcPk4iISDW6YKZFC5NzQQBgsdwcr8h6AfvOeNs6A6w0jZzSLwVqp6WzmQoPlc+Ppdbt6VVpSjp9205IwVcbabVafZwJ4Lylek4NlCVJ8kdlkJwmhFhdVZwvSVKrqu2tYPkLFxERkVfSBTOvv/66STq0gIAAvP7664rbUjwH2I7sE470K8Px4NmoXMEZYCVp5Oz5UqBmWjpFqfB0/QIIqLZ/UFV5Tfq1lvVC7RR8dZ3TrswnSZIE4GMAB4UQbxhsWgdgEoBXqv790lljICIicofqwczf/vY3lJaWokWLFnj99deRmJgITJyoqC2NwX0ZVuYAQ/mZZ0f6rXHWC+nGlQSjAYTs36/PJ6wrL6i6cl9iYiIAYPLkySgpKUFkZCRSUlL05fYucLPVnlL29hsOwB9ABSoX9EWiMkgOt6tX+7JeqHWsVMmZl7C+HcCDAA5IkrSvquwFVAbIKyVJmozKBaDjnDgGIiIilzMKlsLDIYRAjx49sG/fPofblGFjDjDMn3nOhP1Bsrl+7c16Ub1fw/rJFsoNJSYmYuHChZV1qs2hdiQVnrX2lLKnX91xZeDG8coG5QUK+3QkFZ4ax0qVnBYoCyF+QGXqQHMGO6vfmnjygQNGj49u/ltl+eDF7hgOERGpSJdbtqSkBFFRURbPsimtp4S1uawawz4BZKMyI0IUKs86GvYoozK4mgpgnkG9iQAWwL7sE4aPbY4f9gXJM1GZLxgAHjY4DnvbscVS1ovMzEynLiRT2q8M28erlJoXkXGUmv8nvA2vzEdERLWe0tyyauagtTWXVd8nAC1upBE7UfVY16OMG0Hy/Gr1UqrKNYb9wvqZZ8Xjh/1BsuEsWN1xzLRQv9BCuRLOSIWnZr9qHq8nBMl1OS8zA2UiIqr1lOaWVSsHrT1zSpMAFFfbv7iqHLgRXC0zUw9V5fr2DeqHWChXyt6Fe3PNtFFcVa50QaK97E2FpxalWS90CuH48dqb9UJtdT0vMwNlIiKq9ZTmllUrB61RlglJArKyUJiVZTb7hKWWdeW64MpWPRn2Z58ICwszCbJ05fZmtzDJ9VpFmKlvaVqIPdRMhaeYJBllEQnZv78yu0i1co3hOGH+eO2lNOuF2tydl9ndGCgTEVGtpzS3rFo5aO05o2gr165GQT0ZjudbtnjGslo/GlQuQKtermtfSc5g3ThrPC3EgQVuapCh3kVY7OrXjl8o1Fab8jI7cjEUBspERFTrKc0tq1YOWo3BfRnWzyimoDK3rlGfMM21a62eBnbkW1Y5yNS1b+s4ZKg4LcSBrBdqsOvLCKwfr1Lu+lKgU9fzMjNQJiKiWi8xMRGpqakICKi8/ENkZCRSU1NNVu4rraeUDNtnFBMBpOLGhSkiqx5X79FWPU219i3mW1Y56NLVTgQwzaDccHwyXJP1wtlBo+IvI/D+rBe6aTnV/0/4+Phg2rRpdSbrhTPzKBMREXkMpbll1cpBK8N63mDDPLqJABYabLc4NgX1dO1bzLfspKBLRmWquh5V/coG5apnvTCY422Ukm7gQKOUdBphaea0YzQG92W45uIv1ee0G7VveKVFO49VlmUkJCSgpKQEpaWl+rRv4eHhJu+HxMRE/O9//3NJCj5PwzPKREREKpOh3hlFR/sNsTQeJ2RRUNSvQXkh1Ml6YSslnbPIsP5lxNxc7hAL5bYY1te1b+59ZQ9dkFxaWorS0lIAlWnfJk+ejISEBLPvB1el4PM0DJSJiIhgOwuEPUGB2mdQlZAttG+xXMUsChqYDwLNlctQKesFbKekcwYZ6s21drRfjZlyxe1Uvb4NGjTQ50bWKSkpQYMGDcy+H9yRgs/S/0dX8shAuaioCFqtFhkZGe4eChER1QHOWuCmUwh1zqBao3bWC2eQoe7FUJSmpFOLrl+15lqr2a9Sutf33LlzZrcblrslBR/U/RJnTUZGBrRaLQAEW6rjkYFycHAwUlNTkZCQ4O6hEBFRHaB2FgXD2jLUyaNri7uyXiglQ/2LoShKSVcHfykICwtDYWGhcTu6/NhVr6+ttG9ue5+48EtcQkICUlNTAaDIUh2PDJSJiIhcyVlZFGSok0dXCU219l2V9UIpVbNeQGFKOjuDPZtBJpzzS4Esy9i2bZvZfqs/PxoAPWH7lwJLZ4B1x5uSkqLPZKGjS/tm6flxdgo+d6fCM4eBMhER1XnOWLgkQ508uvbSte/qrBe2qJr1oupfmynp7Aj2FAWZBn3JUPcKgzW9HLa9QWZ4eDj8/f3h7+8P4EYqRHNZLwDXpODztCAZYKBMRESkp9bCJRl1K+uFEhqD+4VQ6UwsbqSkiwOQh2p5mxUGe4bPw+nTp5GdnY2srCy0bNnSfBYIqDTXWsXg3JEvBRkZGejfvz/i4uKQl5dnMUh2VdYLTwuSAQbKVo0ZvBhjBi929zCIiMgF1Fy45I65rIBxlgkZN85uGpbruGrBVHUy1Mt6Ye+XAkvBnmGQrNVq9dkg8vPzUVpaitOnT9vVr6Lx2zstBPb9UqDkS4E95a7IeuGuL3HWMFAmIqI6T+25ke7IemEPVy6Yqk4D4HYYB3u6cqVkOJYKz+L0hqrnISkpCcXFxUZ9lZSUICkpyaF+rbHrDDDU/1Kg1nQUZ3DXlzhzGCgTEVGdVxuyXijliQum7OVoKjxbC9xOnjxptr+TJ0+aXVhnrV9bFJ8BtnRclsoVfimw2W8dyHqhBANlIiKq82pD1gulvD1IBpyXCq958+Zm+2vevLldC+uUUHwG2Ez7haj5lwKb/TLrBQAGykRERE5ZuCTDPVkvbEk/fx6agQMBSQIkCbIkYezAgSblnkxjcF+Gegvcrly5YpIyLSAgAFeuXKl1WS/UOvOsJk8LkoE6GihrNBqv+MZMRESupebCJQ2sX9bZXTQG92VY/lnfG8hQd4FbRkYGPv74Y32w3KJFC/j7+yMjI8MkbtBAhbnWbs56ocaZZzV5WpAM1NFAmYiIqDp3LVxyFxnW5746u29ZhTacscAtMTERsbGx6NGjB8rLy80GyQCQBiAbQBaAqKrHdo1f5ZR0npD1Ii0tTX8cUVFRSEuz71lROh5X8nNpb0RERB7IE+dGOpMM9+R5VosM1yxwy8zMtBgkawGUVD0+UfUYqMzjrIS1lHQBAQEWU9IlWzpeteceS5Jxv+fPmy3XPydpaUbHceLECWi1lc9KYqLSZ8X6OHXlBQWu+02GZ5SJiKjOc8fCJXeRoV6KM3dx1wI3nSQAxknkKh8n2XEMzkpJZ/NLQbW56IVZWcjdv9907rqd/Zo7juLiYv1xKMWsFwoUFRVBq9UiIyPD3UMhIqI6wB0Ll9zF24NkwPzc7xAz5c6aZmA+iZzlcnN0/dpKSaeB9bnuhuWqznmGfanwrB2HUq7+ZScjI0N31jvYUh2PDJSDg4ORmpqKhIQEdw+FiIjqgOpZIHL370d0WZlJeW3g7UGyPZx1cY0IC/1ZKrfGVko6pez+UgAbc55hXyq8iAjzR2+p3Bxz2VfMZmVRSUJCAlJTUwGgyFIdjwyUiYiI3EGG9eChNtAY3JfhWVkvZFlGWFiYyYItS+W2KJ5mYGeQmQIgqFpfQVXl9lCSkk4pk2AyKwuFWVlmg0wZChZCViu3lQovJSUFQUHGz0pQUBBSUpQ/K574JY6BMhEREZQFD7WJDPdlvbDEHVdIlGFfkAlULthLBaALbyOrHtuzZM2RlHTWKM63DDvnPEPZl8fExESkpqbqjyMyMhKpqal2LeRTOh5XYqBMREReT6OpWX58Gc6duyvDs3ITy/DMrBfVg+Hk5GScPn1atSskZgLYZ6bckSvuJQKIBRAHIA/Wg2QZpq+/YRCeOHEiYktK0ANAeX4+Mi5ftnuageE4ZQCnUXm81cvtWggJ+7486lLrxcXFIS8vz+5sF7bGqStX3I4Kv1AwPRwREdV59i5c8mYyPHdBn9IFdxqTPS2TYf8ZVEtnnu3p1xZLwXn14NZeMuw7Xkv92tuOWpT0q5Th+0eX47mkpASDBw/GjBkzFH354hllIiKq8xw5o+itPDVINmRrzrDiduCcaQZqqN6vKpfDhnpfCtzx5dHaOB3p1zBINszxXFFRgfnz5yu6IAoDZSIiqvM0BvdlqBO0eCpvD5LtWuAGx+bohlgodwZd+zUNznX11fpS4I4vj2p/idO9f6zleLY1/YKBMhERURUZzHqhK3cHtfPoqrWQzRVnUEMslCul9peC6vVd8eXRWc+/tRzPtn6hYKBMREQEZr2oXu4Oal9swrC2DPvn7nrVNAOF7Tha7o7pKGp8iZNlGZKFRZGSJNn8hYKBMhER1RmyLJv9qVWG8YeybHBz1bQEXX9q1bPVhrOyXshwfHzOuiKbDMcWuLl6mkGhhXIlDOvLFtqxVG6pX8P6+3DjdTUsV8qe7BO2xq+U7v0zY8YMkxzPABQt6GPWCyIiqvO8YYGbmjSovPyx0nJXYdaL2pv1wt4rBqqd9aJLly6YPHkySkpK4OPjgxkzZmDOnDk22+AZZSIiqvPqUpDsDZj14ka5oTQA2QCyAERVPUa1+s7MemGpXAm1gmSjcktnpKvKDdtPTExE586d4efnh82bNysKkgEGykQ1VtMLHRCRcym56IDGsByetcCtrmHWC+NynTQAWgAlVY9PVD1Oq1ZfrS8FBTANhjUWyu1h9fW1Z/wG7SQnJ+vbstZ+bm4uoqOjTcqtYaBMRES1lr1ZFGR43gK3uoRZL0zLdZIAGCc4q3ycVHXfXV8KzH0JLSwstPjl1Orrq3T89v6/riqPjo5GSEiISbk1HhkoFxUVQavVIiMjw91DISIiL6ZmkFzbsl54Ima9sFxuPsHZjXJ7gnwN1DtjXP11KSwsRG5urt1BbHp6uvLXy0w7hYWFNtuvHiSPHDkSt912GwAEWzo+j1zMFxwcjNTUVHcPg4iIvJyaQbJhOTlH+vnz0AwcqH8sw+D5Nyi3l1E7CsoLoc7COmvsDc4jUDndorqIqn8N68sW2pHNtFtTSqc32PslSAPLC04tBeeZmZl29btu3TpoNBpIklRk6fg88owyEdWMrTmZRHWFPVkUnDUnk5TTGNyXoXyuuGyh3LBdpa+vDPdcXKMQ1lPSpQConuAsqKq8Og0sH291MtQ5NlvTG1T9paAGwbmlM8+WMFAmqmXUXAhDVFuolUWBXEOGe+aK69r3xKwXiQBSAQRUPY6sepyo8rgcYW16g9pBsrl+lQbnlqaFWFPnAuW0tDRkZ2cjKysLUVFRSEurnlyFyHs5+48SkTfil0fvIsN9c8U1qDwTG2Kh3Jn93m6hX0OJAGIBxAHIg+cFyTWZbmEr+4S1fi0GyQMHApIESBJkSULu/v2ILiszKbemTgXKaWlp0Gq1KCmpTK5y4sQJaLVaBstUazBIJjLGL4/eRQbninsbRxfW1eSXHUez2TjyS0GdCpSTkpJQXGycXKW4uBhJSUkW9nAdziklNVT/49CzZ0/Ex8czGKA6i0Gyd2GQ7H3UzHqhlF3BOWqWCq9OBconT5pPrmKp3FX4syCpRekCB6K6gkGyd2GQ7H2clfXCGsXBOWr+5atOBcoRERF2lbsCfxYkZ7C1wIGornDG3EhyHo3BfRnKs16oTXZRP7WJR2a9QM3zY9epQDklJQVBQcbJVYKCgpCSYi65imswSPY+0h7jW9blylv1cndRssCBPIetaVcM4tTBrBfeRYZ7sl6QYxxZWDd24MAbebOryh3t12JwDvtS8JlTpwLlxMREpKamIiCgMrlKZGQkUlNTkZjovnWjDJJJTfyFwrvw9XKcqpfN5ZdHjyKDV0j0Jo4urKvpLwWKgvNq7TuSH9sjr8znTImJiVi4cCEAz/i5TenPgvywJCXUSK5OrsMg2XGGz09aWhq2b9+O0tJStGjRAqdPn9bX45cR7yKDC/q8jaML6zQGbejKlabgkyXJ7BUbLZbD8awXdS5Q9lS2/pgXFDgzgyPVFkov60megcGb4wyD5MmTJ6O0tBQAkJ+fD61WCwAIDw+v/HkX5j80a3pZZFIfg2Tvo/hy0lDvlwJ75h4btp9sodzaX4E6NfXCU/FnQVILs154Fy40q7np06frc+PrFBcXY/r06WbnKMpg0OXJ+Hp5H6UX9VDzS5DSuceW2renX48MlIuKiqDVapGRkeHuoThd9Q9H3Y1nlmoX3euqVj1bnJX1Qq3xkTEuNLPM2ntOlmXk5+eb3Zafn292jqK7siiQMhqD+zI8//WSod541GzLHWRUvl6ZAPaZKXckWNVY2Fa9HUtzj22def4LQOXvTwi2MATPDJSDg4ORmpqKhIQEdw/F6TinlNTGrBfehb8oOUb3/LRo0cLs9hYtWtj1sy95Hg0q56xqzJSTZ5FhXzBsqdzRfi3NPbZ15jkBQGrlpiJLfXhkoFyXKE2aTaQEFy55F75ejtM9P6+//rpJ2s+AgAC8/vrr+scymEWByFlkqBMkyzXoN0Rh+8x64YWUzinlhyUpwV8ovAuDZMdVf37+9re/6bNevP766/q0nzK4QIzImdQKku3KeuFg+8x64cVszSll1gtSglkvvAuDZMcZPj/h4eEQQqBHjx7Yt2+fvtzSlblkMEgmUotaQbKzsl4A1qfxFACQrPTFqRcegHNKvdg3acCBbGB3FjAiqvKxGzHrhXdh1ouas/fKXDIYJBOpSWNwX4Zr5iob1k8DsB3AfgC+AE5b2MdRDJTdjHMUvdg3acAcLXC9KjXV2ROVj90cLAPOy3pBzsGsF45x5Mpcnp5FgchbyXDdgj5d/TQAkwGUVj3OR2UWizSo9/+agbKbcU6pF3s3CbhWbFx2rbiyHIC058Yt63LlzbBM2uOcYfEXCu/CrBeOUfOyufwqQlQzMtyT9WI6gJJqZcVV5eb+X8uwP4DmHGU345xSL5Z/0r5yF+AvFN6Fr5fj1LxsLr+KEDlOhnoL+gwfK+nXfBb1yvKtdrRlDQNlN2PWCy/WIqJyuoW5cjfhLxTehUGy49LPnze6/HQhKlM/ZcLK5aoN9rdUTkT2sff/lwbms1tUr2eNrv0WMB8st7CzPWs49cJDcE6pF3oyBQg0zt+KwKDKcjdhXm7vwiDZcRqD+zLsuzKXpXIisp87/n/p2n8dQLVPYQRUlauFgbIHcNecUlmWERYWZrKq3lI5VXN3IjAzFagXUPm4ZWTl47sT3TYkZr3wLsx6UXMy7Lsyl6VyInKMxuC+DNcsmNWlektE5ZX1qj6FEQng46pytTgtUJYkaZEkSeckSfrZoKypJEmZkiQdrfq3ibP69xbumqPIuZEquTsR6B4L9I4D1ue5NUg2JMsyEhISIEkS9u/fj6ioKKSlpfEXCg/FrBeOkWHflbkAy/lUiahmZLhnwWwigFgAcQDyoG6QDDj3jPISAHdVK3sewGYhRAcAm6se12numlPKILn20gXJpaWlKC2tTJpz4sQJTJ48GQkJCcyi4GGY9cIxMnjGmOomGZ6Z0lAD6xf18FZOC5SFEN8BuFiteBSApVX3lwIY7az+vYW75pQySK6ddK9jgwYNUFJinDSnpKQEDRo04OvrQfjLjuPsvTIXEZEjXD1HuYUQ4kzV/bOoXJhooqCgADExMfpbamqq60boYu6aU8q5kbWT7nU8d+6c2e2Wysk9GCQ7zlyQnGumnIjInFQAMVU3AKGW6rktPZwQQkiSJMxtCwsLQ05Ojmp9hYWFmT1za65cTU8+cMBmnXer/rWV9aKgwDk/XNg6o+Wsfsk5dK9jREQETpwwTV0XEeG+1HVkikGy4zQG92VYznph+JiISEdbdQMACThvqZ6rzyjnS5LUCgCq/nXJ6S1PT5nlzqwXnBtpmUaj8Yj3hz10401JSUFAQIDRtqCgIKSkpNTaXwq8+fUCKv/fxcfHIzw83OYvO954rM4iw3rWCyKimnB1oLwOwKSq+5MAfOmKTj05ZRazXpAzhIeHw9/fH/7+/gCAyMhIpKamIjw8nFkUPBDzqDtGhu2sF0RENeHM9HDLAewA0EmSpFOSJE0G8AqAeEmSjgIYUvXYZXQfOpmZmdi3b59JuTuMHThQf4WpNElCdlYWsrKyMHjgQEzVXXlKktTvl0Gy63yTBhzIBnZnASOiKh+bkZaWhuzsbGRlZenTuTlC9zpmZGSgf//+iIuLQ15enj5IdvSXArXGR8bc9YuSN7D2npPBrBdEtVkagGwAWQCiqh67g9PmKAshHrCwabCz+rRGyRnUVW5YM6j7Y56GyrkyujwFFQDmA+gC9XMCApwb6TLfpAFztMD1qlf27InKx4BRzuW0tDRotVp9pooTJ05Aq62sl5io/B3grF8K1BofGeMvO5bZes8x6wVR7VU9JjqBG/OJXf2JUyeuzOfJH0a6XpMAFFfbVlxVLjujX2a9cI13k4Br1V7Za8WV5QaSkpJQXGxcr7i4GElJxvVscVZebrXGR8bclUfdG9h6zzHrBVHtZS0mcrU6ESh7apBs6KSVcmdOCuEVwZws38IrW6385Enz9SyVW+KshatqjY+MefpCY3ey9Z7TGJTJsJz1goi8j7WYyNXqRKDs6UGyDMDSLGQJzluQwqwXLtDCQjq2auWW0rbZm85tYGMNpD2ovKXK2P9zLsraRZuU20ut8ZExT15o7G72vOc0AG6H8YI+XTkReR9Lnyzu+MSpE4GyJ08zkFF5xngGgCAz22fAOX/sPXk6Sq3yZAoQWO2VDQyqLDeQkpKCoCDjerp0bvYQvSWI3hK29pYQOnUgepSXIe7IfpNye6k1PjKPWS9M8T1HVHelwDQmCqoqd7U6ESjreOI0A93CkzmovEqMLvOtDyrn4sxxVr8Mkl3j7kRgZipQr+qVbRlZ+fhu4+UIiYmJSE1N1ec+1qVzc2ShnAz1U2apOT4yxqwX5tn7npPBqRZEtUUijGOiyKrH7vjEcduV+VzN07NeAJVvgP+hcq5dJtQ5k6w71upnzD0hSLY0tlrn7kRgzcLK+1amPSQmJmLhwsp6tp4TS8+dDOelzLJnfGry9PeJ0vGZq+fsX3bUfu5qcqyO1AkPD0d5eTl69OhhktJTY3O0ROTNEgFUfXK69UtwnTij7MnTDAx7leG6BSmePB2FHMeUWd6FWS8s43QUIvIEdSJQ9tQg2ZAM91yG1ROno5DjmDLLuzDrhXmcjkJEnsIjA+WioiJotVpkZGSo0l664RXuJAmyJBldEU9X7i4y3HMZVma9qH00BvdlMGWWp6v+dyl3/35El5WZlNclnvwLIBHVLhnQX8gk2FIdj5yjHBwcjNRU9SYMawzuy7A8h9MdZLjnMqz8MKrdZFj/haLA5SMia2Tw9dLhdBQicpWEqttCoMhSHY88o+wsMqwHpe7grjmlDJJrLxnu+YWCHCODr5chTkchIk/ikWeUnUGG7TO3q1w9KDPjKYS6WS8s9ssgWT0OXMDDWWS45xcKcowMvl7VKb0Ii8ZkTyKqq2Qntl0nzijL8NwPI8N+ZTDrBdUMs154HlmWERYWhsLCQpNyvl6WMesFEXmCOhEoe2qQbEgGs154krS0NGRnZyMrKwtRUVFIS0tz95AUqatZLzz19bKZvQE3Xpc0ANsB7AfgC+C0hTY99VjVxKwXROQp6kSg7C1Bcm3KeqHRaFw+hUOtPtPS0qDValFSUgIAOHHiBLRarVcEJBqD+zLqRtYLe18vpe+Tmr6fFC2Y1R0DgMkASqse56NyJXYajF8vZ7031fi/oztzbnIRHAvl1trhGgoi8hR1IlDWGNyXwawXAD+MrElKSkJxcbFRWXFxMZKSktw0IvvJcM8vFO7gqa+XPdkbpgMoqbZ/cVW54evlqceq5t8TZr0gIk9SZxbzAbaD0vEuHxGzXniikydP2lXuaWTceP8kWyivTTz19bKUvSEzM9P4/x0qzyCbkw9gq8FjTz1WNYNbfX57XTswWOBsUE5E5Ap14owyoOzMrTu4a04pg2TLIiIi7Cr3JDKM3+eywc2TphupyVNfL0XZG1D5urSw0EYLGL9ennqsaqZ0M6wto25MHyIiz1UnAmUZnrugz7BfGcx64QlSUlIQFBRkVBYUFISUlBQ3jUg5T32fO5Onv15Wszeg8nV5HUBQtf0CqsoNeeqxKk3pZg8ZdWf6EBF5rjoRKHtD8CCDWS88RWJiIlJTUxEQEAAAiIyMRGpqKhITE908Mts8/X3uDJ78einNepEIIBWAf9X2FgA+rio35MnHCqiX0k0GL8JCRJ6hTgTKnh48yHDOh4KtNFLOzHpR07G5W2JiImJjYxEXF4e8vDyPCUQA68+dxqCeDM9auOoIpe8TNV8va3mPw8LCFI/P8P/R6dOn9fVatmyJhIQEo6wXABAOQADoAeAsbgTJsoPHqvb/MaV/T6ZOnYpDhw7p682cOdPk70laWhq2b9+ufz6Mnjd4x8kNIqob6kSgrDHzuMBMuTvIcM6Hgq00Uu7MeuHN6dfcTelzJ8PzLtduL3veJ7IsY9u2bWaDW3vYk7/X1vgMg2TDevn5+SgtLcXp0zcyJctQ9xcltf+PKf17MnXqVMyfP9+oXkpKCqZOnar/e5KWlobJkyejtLQyGV5+fr6+LV6EhYg8TZ0IlD2Zsz4UbKWRcmfWC09NceUNlDx3Mjxz4aq9lL5PVPu5384vj7bGp8vekDRxokm9kpISJE2cWNk+1P9FSe3/Y0r/nixbtsykHgAsW7ZMf3/69On6QNqwrenTp5tchAWoOxfNISLP5JHp4YqKiqDVapGQkICEhATF++k+xLxpEVr1P/7JUOfMia00UtVTMMmG/RqmYBLCaH81nmN3p7hyy/skVVlftsZk67mT4dyfrdV8zmy9DkreJ4ZBbHJystlypexNcWZrfLralt7VJ+G810vt/2M2/55UPT8235+yjPx888nw8vPzsXXrVpO/S6dRlRrOoK4MBs1EdYHs5PYzqm4Agi3V8cgzysHBwUhNTbUrSPZWGoP7MtSbU2orjZSz+lXCU1NceQNbz11tmttp61jVnj5kb4ozpe/j5hb6aw7n/aKk9v8xm39PBg4EJAkR1b5Y6+sJAVmSMHbsWLRoYT4ZXosWLcymzvPm6UNE5NkSULmQGkCRpToeGSjXVRqYnztd/bESStNIyXD9h5GnprjyBraeu9oSJAO2j1XtK7jZm+JMyftYBnAFlaneDAVUlTtrmoHa/8eUtjfRzL5BVeVjUfll5PXXXzdpKyAgAK+/fiMZngzbX/qIiFyBgXItpSSNlAz3fBh5eoorT2brudMY1JXh3VkvbB2rmhe5MKR0zrOt8cmofP4zUJnqTRcst0BlGrgMOCePuizLeOaZZzBt2jSjsU2bNg3PPPOMQ9NnlP49WQAgCTeONRLAtKrydFR+GdG15e9fmQyvRYsW+Pjjj43a0sDySQNPWYhNRHWDR85RJnUkJiZi4cKFAEzngcpw3s/0ujRSJSUliIqKQkpKikkQbG1szuqztlDy3Mmw/voWOG94qrJ2rErPABueYTZ8n/j6+mLGjBmYM2eOSf2pU6di3rx5+vfTxIkTsWDBApM5z9bGV/35X4jKM8anYT5Itpb1QunrVf14f/jhBwBAcnKyxTPtSv/v2PP35Ieq8mSYfx+Gh4dDCIEePXpg3759RuM3Hh0RkXsxUK6jnBkkm0sjBcBpgatafUp7qhVctlDu4WTUrZ+tbZ0BNgySDd8nFRUVmD9/Prp06YLExESbKc6SkpLsm/MM89MqzC1M070uyRbKlbJ3Oopa/3cUz7WWJOMvBfv3A5IEwPu+xBFR3cCpF3WUs+ayuiP1G9PN3TBm0CLcHdAEdwxahFUPHMCTVbfq5bWFkrzHuuDQ2vtElmWMHTgQ6efPY1lKivkUZykpN4K6qouPmJxZNSjXGJbD8rQKNb+02jsdxZ4UfNYuwqJ0rrUMXnGPiLwLA+U6SmNwX4Z6c1ndkfrNbenmcmRgcFjlv+bK3WDDtum46/bXEd6ij77sdP4us+WeqnrwWVhYiG3btpkGpZKkD241AwcCWVlAVpZpeRVr7xPD/L3W0rnpxmdPtg0Z1qdVqJn1wt4Fifam4LP4ZcRgXxmu+VJAROQKDJTrOBnqZr1wR+o3t6Sby5GBf48FXk0HYjTmyxWQ9ii7KeXtQTJgx9xjKA8yZVmGVHU2uDpJkoyCPUvvmghU/b+wY3qDDNtnUKuPX42sF0oXJKqdgk+G674UEBG5Aucou9mTCn8Gf7fa45peNOPJBw4YBVGrWvTBqqpthuX2SklJgVarNfo5t6ap32wdqzP6tEpJkKy03PCxCpQEyafzdwHopmq/zmBz7jEUzgGuqj9jxgzMnz/fZKrBjBkzjIK9iQCqv3OqpzgzN70hMzPTZHqDBubn3OrKn3zggNH/u5+2atE0+GasunuVUfkaM21YYs9FWGz933H0S8GNXm1/Kaj+ehEReRKeUa6jrAVRNTkD6Y7Uby7vU60g+d/Ou2yCrdfX06n6c39V/Tlz5hi9T3x8fJCUlGSc9QLKUpwZjsfW9AYldK9L0+CbEVCvkUm5UvaeAQ4PD4ePj48+VZvu/054eHjl3GOFc55l2DetQoY6qfCIiJyt1gTKuhRHWVlZiIqKQlpamruH5NEcCZKVPseJiYmIjY1FXFwc8vLyXJKmTWmfqrxP1AqSFU7PsJeSL0HuZu11UP3nfoP6iYmJ6Ny5M/z8/LB582aTIFkX1M0BEAsgDsASGATJhu0rnN5gi+HrYi5Ituf1susMcNU4MzIy0L9/f/3/nfDwcIv1LX0p0EB53mMZ1l8vIiJPUiumXrgjJdn/t3fuYXJUVQL/nSQQMiQkxJAhjiAYQDEfBgloXHz0gPISzaJmQSeKiiY+oohRFiTLxl2zK2o237roakLgA4ko4bEhKi/ZFK6sPAQTSORhQgISwkCAyQ4ZCCRz9o+6Panpqeqq6qmZqu45v+/rr6tOVd0+dfp231P3nntuvVOLk1zvNk59D0u88IKycpIzDruAgRspyJK476HW4f490gp5YEKfR6Ans0KeNuY5NrwhUH4Uld/LGc1XhMorieqxrpzA2EEgvCF4v6tXp3sYiXkoSJznmexS4RmGYQwGDdGjbOnB0pM8ltWnnm1cnhA36xvh9zDrGxfVliu5YE4y1EfWi3nz5oV+D/PmzRvw4f60E81CU5yl7PGuRtbfVy89iQ9HSdPzXC0MJimW9cIwjHqjkI7y9u3bmT17NqtWrUp0fm7pwRqApLGsRbVxkny2PbRH6Bolr0YBnWQoftYLz/Nob28PPdbe3p5quL8WJzlp9gmPKk5mMPWcCNx5Jx133tlXnoCk31eqek6Ch4KkDyMpU/DFkTQ8wzAMYzBYBfjjmYyNOqeQoRdjx45lyZIlic+fOHFiaOM7ceJEDjjgAFvpKYI0sawHH3wwTzzxRJ8yBjQFWwyJe/bKvcXNB8Mzfe+B5nT3sHqa7HGu5kQM389pRe5X/0BOWS/K++d85HeZfkZ/mDlzJs3NzaG/1+bmZv/7qlipLWoFt1qc5DTytD3PtWRvSDqyM3Pm+aljtqtmnxiAFHyGYRj1xofcaylsjzqnkD3KafA8jx07dvTMZC8zcuRIduzYkWpYcCiRNpZ14cKFNDU19SpjQFOwJSDtcr3MXQj79L4H9mny5Wk+lxTDx3E9zFGLlVTKG4QVK1awaNGiPnVp5MiRLFq058HMI30PcBp5R4Qc9vRwjouQV5afRfaGuN9jLU5ypT7VYq3TpuAzJ9kwjKFCXTvKwVnby5Yt63GWm5ub2WuvvVi1alXiWMGhRtrYyDzSvsWxbeEKWvcr9cQg37m1g7Xr1veR93BqG8xfAnu7h6oD3+jvn5ruHjJzklPIw4bXOzo6QuVFplQq9dSlckqy5uZmli1b1lOXPNLHspaIHtYfsEU9yCZ7Q5KH1sQT7qhtQmJ/U/AZhmE0KoUMvUhKZWOxdOlSOjo62LJlS8M5ybVOIIpyomqJZW1ra2Pp0qVVy01Dv8uodDI3rofJU6qHN5zaBjf69xCZ1SKGQGlVnboJc1qzCc/4+5msuLH397tgwYLUE8eKREtLC6rK1KlTWbNmTY88zNnz6N+Er+D5HtHOXppyo/QJysOuqSTtyE7s/0CFPgui9JRA+FAwG0aU3JVTGV5SlleGt4XdayVJzjEMw8ibuu5RTjohpREINo7BHLQnnngic+bM6dNoxk38qWyUV66ezb6jJkZmvciLRHmPy07m5CkwZlxf+QDhkU0YQJIe5qQ9ikWhVCpF/gbTDvd7ZBMTWy6nvz3AUfpEyatRjiGPii1P5SSH9ABnEY6SpBzDMIxGpa4d5YFYIauoBJ3kYA7a7u5uFi9e3ONApu15zmpFsIEgKt9uL2c56GSGOcmVi3r80YM1d0FnR195CjwGz0kOyuO+34Ei7uErcTmBbAnj1q71syiI9M6iEDyfaHumpUT12OOkZOUkp6GWEaUS/Q9H8bCUboZhDG3qOvSiTFbJ8ItMuRGsls+42opa1ZzkU45fxL3rfhwqz5Nq99rW1tYn+0TJnbOgQt4nvCGq5/mOZDXFI7vsCkmdZHmgt7x1v9KebB7B871Et5CKJE5aUrKyW9hw/2CRh9M4s7U1MhyiUp6EUmDbo0oqPMxJNgxjaFPfjrJI71RI27aFyhuJavmMa3WSW5qPg3UR8hyJy92cKptB0JlcsiBcnpAsnb3VHduqxjCXkZ+uTuZUe+sS30cSao2Nj6IRhvtLIfsD7bTnFY4Sdl+lCLlhGEYjUsjQi87Ozl77cUn1h0qPh+d5SMRCBiKSajnaMGd456udmTnJnuf1O+wlKkdzWR42cWlLiDzrRUCycpIHouc5a7J0kiF59okSA7c4hUfyyWZJzhsMSoFtj2zCUYLlrAlcH5QbhmEMdQrpKI8ZsydWNqtUSPVOOYbzwu5umkKOX9jd3XtYNqYnMMxJfmH7xkL0JJeJy91cCsg9UjqZnR01O5mVZ5cId96y7nnWOa3oNEGnCaunCRPmtPaRZ02WTnIQD0s1Vgse1etJVuU02v+nYRhGrRTSUS6TNhVSB42bDL/ceH0HWAKUl1cZBlzk5EDfiVGtrf6KZhXyyuwWL2zfyPixkwuR9aI8gtDS0tIrd3NzczPDhg2jpaWl9/nEOA9hTvLG9QO+nHQuPc8ZMxDZNjyyyT4x1PDILhylhC0nbRiGkYTCOsq1pEJq5B6qYOPYBrwFP8D8DgJOMumdqy3t9/GrO+ciwPMdj3LlypN4bPOvc816Efze29ramD59OlOnTmX37t198mN7xN9vZU/siMfWMnX3rj7y0BXx+rFSXilkv9IJSaJ/GnnmhDxkhT18BamW0i+o/7gQ/W24PxwP6wE2DMPIg0JO5uvs7EwVG+mRfkJKvVEKbHtkM0u97CSr7qJbdwHwUtdW/vueixEZwenvu7SXDp7ncddddzFlypQ+8izzVkflx7799ttTO8lh8siezDxigKvoWYv8ugHQMcnnllm+fDnnnHNOn5R+4C8yUsImiNVCCbObYRhGHhTSUX788cf7OkUxTvIK/MlcYfJGwqO6s5fGSb7lrnnsNWIUL+98vtdn7O5+lVEjx4QudhCZgu+57JrrPqnP1q2HyVP6yCeE3JdHvB0WRMhbkzjJNfQsVyPrMIysHeW0+sybN6/HSS7T1dXFvHnz2L17tzl1hmEYRl1RyNCLibt28fPWVlbFDPtC7+Fsjz2hFkF5o+ARP2xdSigvZ7d4eecLoZ8VlAcfUtasWdMTk5o2j26SFQN7EZP3OGmMpUdvO3iBV+iwdVyWjAwphez3d2JgTXqUSqGjAqnqlefR3t4eWn57e3ufeuLROKFRhmEYRv2xCvDHPBkbdU4he5Qn4U9Yg+rOXnC/GnM//lA2iuVMieTDrx7VnZxydovRTQfyUtfWPmWObjrQPz9JrHjCsIuk5ei0ivzYHX5+bG9abfmxUzmZSVLJBfIbDxZ5xaamcZJnzpxJc3NzqLPc3NzccMvKG4ZhGPXNh9xrKWyPOqeQPcplPAZ54lKD4BHvHJbDKt419VyGD9u71/Ujhu/Du6aey5b2+7JdbCJlzHlWzmHScqJSrlXK86AUsj8YIybB8j2qfC/ue1y0aFGflH4jR45k0aJ8V3k0DMMwjFooZI8yZJsKaahRIr7nuRzLuu+oiYiMYJh00627GN00iXdNPZd9R03klrvmcfPOF5Mtm6sar1fCVGNZ58euPL9EuH0sq0B1SkTXq56wKMdngNeAZmDRzp20zZoFs2YNgpaGYRiGkR2F7FHuxJyWwaA8oe/0913KgQcczesnHsvZM27rcZJPOX5R4h7FJMgD7rXEo/WMmWxbuILW/Up95GFO8mDkxx5K9S0YEx5M5zZ8+HDmz5+furxSYLsFUGAq8Ax+OkOweGTDMAyj/ihkj/LjwO0MHaelFrwMymhpPo5zPvI7f2ddhDzwedUeXpJkM6iMPS7Nieiprih/C+H1IbifBZXl9zc2vsiUe/KXL1/O7NmzezJVdHd3s3jxYo488kja2tpiSumLR+OnajQMwzCGDoXsUZ6KrRhVJDyyCYNJUk6p4poS0SuIDRRR+pTljUA5e8xFs2bR1dXV61hXVxcXzZqF5zLLJMXDFhMxDMMwGotCOspGcfBInyIsKgVcPYTTeAyt2Pgnq8jTPhSUsGWRDcMwjMaikKEXRnFI6yRXzZKRsJw8KTF0VkDzAMGPJ65EaLyHAsMwhh5J08Nees1RA6yJUa+Yo9xgZJ0zOqnTKEu8njzDfVbWc/JSRWq1sHLqgUb54y0BV+EnWw8GXzTh5zEvDb5KhmHUIY3yn2gMLEnqSRHriDnKBgBnnHhFv66fMKe1+gS9nPIPG9UpT9c7B9gJvBFYGJDXitfP6w3DMIyhRxEfunJxlEXkFODfgeHAZar63eDxYC/jcuAi/JjJg8mmEW9k1m9YwZTDBn/KWT2EVcSxhJ6lLIcUbcBSt+3VWMZQsl2Wf+R5/V6LTFL7JrVd1g1qlqN2een2tmuOSvR7zXqE0mxXDJLaLmunNcv/u8F0qAfdURaR4cCPgA8ATwH3ichNqvrn8jnb3Ptyeg8LP8GexjiNs7xpi8ehLaX+qD1g5WWt2/qN12Xa8H7wvf+RSL9gZfQId5JX4S8VmRVZlzd//Ft58ORfZlZeEv2S/thPvuaoTO+18nv9yx2f8fWpGFlI+ieTpaNc9HqS5W82r99rUpLWu6T1OA/b5aFb0vLy0i3r/7o82sRGsF1evknWtkv8cHvrmZn932WtG9ccNTbqUB49yu8ANqjq4wAi8gtgBtDjKL86fB/m/t19XLnyJLq6tva6uAv4QtMk/jDjtsSN+OaMK2OW5SUtK/GXfeuZ/VOoglrutUR47HFShyXpvb6asfOYNRdM/ii3vmNB4cqC7H8TT45/K3MTNB5JfrNZ3+vqexdka7t7F2T3m034e036m0iqW1KKXu+ypBHaibwo8r0W2XZFv9ehZDt6ZzXthWiCpYezREQ+Bpyiqp9z+58E3qmqcwPnvALsxp9XFMX9KT52LLC9BnUHo7ysdZvAnk75LCjyvZrtilNelrYr+r1mWV6R61zRyxtKtivy7xWKfa9Ftl3R77XRbTcBOMBtD1PVUWEnFXIyn6ruk7cOhmEYhmEYxtAmjwVHtgAHBfbf4GSGYRiGYRiGURjycJTvAw4XkUNFZG/gLOCmHPQwDMMwDMMwjEgGPfRCVXeJyFzgVvz0cJer6vrB1sMwDMMwDMMwqpFHjzKq+htVPUJVJ6vqQhE5T0TWi8g6EblGRPZxPc73iMgGEfml630e8ojI5SLyrIisC8jGi8jtIvIX976/k4uI/NDZ8EEROSY/zfMnwnbfF5FHnH1uFJFxgWMXOts9KiIn56J0AQizW+DYPBFREZng9q3OBYiynYh8xdW79SLyvYDc6pwj4vd6tIjcLSJrROSPIvIOJ7d65xCRg0RktYj82dWvc53c2okYqtjO2okYomwXOF6/bYWq5voCWoBNwCi3fy3wafd+lpP9BPhi3roW4QW8FzgGWBeQfQ+4wG1fAFzitk8DbgYEmA7ck7f+BbTdScAIt31JwHZvBdYCI4FDgY3A8LzvoSh2c/KD8EeGngAmOJnVuRjbAa3Ab4GRbn+ie7c6F2+724BT3fZpgBfYtnrn22IScIzbHgM85uqWtRO1287aiRpt5/bruq3IpUc5hBHAKBEZgZ8SbitwAnCdO34l8Lf5qFYsVPV3wAsV4hn4NoLetpoBXKU+dwPjRGTSoChaQMJsp6q3qeout3s3/uRS8G33C1XdqaqbgA34OcCHHBF1DmAxcD4QzDFpdS5AhO2+CHxXVXe6c551cqtzASJsp8B+bnss8LTbtnrnUNWtqvqA2+4EHsbvkLJ2IoYo21k7EU+Vegd13lbk7iir6hbgB/irVG/Fz4t3P9ARqJhPscfgRl+aVbW8MsszQLPbbgH+GjjP7Fidz+I/4YLZrioiMgPYoqprKw6Z3eI5AniPCy27U0SOc3KzXTxfA74vIn/FbzcudHKzXQgicgjwduAerJ1IRYXtglg7EUPQdo3QVuTuKLs4qRn4wxavB/YFTslVqTpG/TGNwV1FpgEQkYuAXfgrpxtVEJEm4FvAxXnrUqeMAMbjDzd+E7hWRCRfleqGLwLnqepBwHnAspz1KSwiMhq4Hviaqv5f8Ji1E9WJsp21E/EEbYdvq7pvK3J3lIH3A5tU9TlVfQ24ATgevxu+nJXDci1Xp708ZOHey0O5lrM6ASLyaeB0oM01IGC2q8Zk/AfbtSKyGd82D4jIgZjdkvAUcIMbcrwX6MZfIcpsF8/Z+G0EwAr2DHOb7QKIyF74zspyVS3by9qJBETYztqJBITYriHaiiI4yk8C00WkyfWqnAj8GVgNfMydczawMif96oGb8G0EvW11E/ApN7t0OrA9MPRmACJyCn7s1IdVtStw6CbgLBEZKSKHAocD9+ahY9FQ1YdUdaKqHqKqh+A7fseo6jNYnUvCf+FP6ENEjgD2xl/W1epcPE8D73PbJwB/cdtW7xyuHV0GPKyq/xY4ZO1EDFG2s3YinjDbNUxbkccMwsoX8G3gEWAd8DP8GaRvwq9wG/B7DkbmrWcRXsA1+LHcr+FXunOA1wF34DcavwXGu3MF+BH+TNyHgGPz1r+AttuAHye1xr1+Ejj/Ime7R3Ez7YfiK8xuFcc3s2cms9W5GNvhO8ZXu/+7B4ATAudbnatuu3fjz2FZix87Os2da/Vuj93ejR9W8WDgf+00ayf6ZTtrJ2q0XcU5ddlWiFPYMAzDMAzDMIwARQi9MAzDMAzDMIzCYY6yYRiGYRiGYYRgjrJhGIZhGIZhhGCOsmEYhmEYhmGEYI6yYRiGYRiGYYRgjrJhGIZhGIZhhGCOsmEYhmEYhmGEYI6yYRiGYdQ5IuKJyPy89UiLiNwsIucPdR2GGiKyWUReEZH1OX3+bSLysojsijvXHGXDMAzDKCgiMkxE/ldEVETeMIifOyiOt6qeqqrfG+jPGSwd6vWBJSc+p6pTggIRmSYi14vIsyLyknOorxeRE5IUKCIrReSqiGOrReRSAFU9CTg1SZnmKBuGYRhGcTkP6MpbiUZERPbKWwdjDyLyAeAu/GWtjwXGAEcBPwfOSFjMT4GPici4irIPB97njqfCHGXDMAzDKCAicgTwJeAbNVzbJCI/EJFNIvKCiNwiIocFjm8WkYtF5Peu5+6PInKcO3Yp8B7gH9yxR538dSJylYg8415Xisj4ijK/JSJ3uOvWicjfxOjZ0wMbd72IfFlE1lRcf6iI7BaRQ0TkXBF5REQ6ReRJEflXERkecs+rReQl4KMhOiQpI1THKLtF3PdmEZlf1kVEHhKRt4nIx0Vkg4hsF5HLRGREQr2+6r7rThHZIiL/kuRYf/XKmP8ErlbV81X1SfXpVNXrVfUrAR2r1e1bgOeAT1aUPRu4R1UfSquUOcqGYRiGUTBEZBhwOb6T3FFDEUuBtwDTgQOBe4BfVfSifgE4FxgPXAf8RkT2U9W5wP8A/6yqo1X1ze785cD+wJHuNQH4WcXnfhb4KjAWuB24MqXe1a7/OfAWETk6IPs04KnqZuAp/OH0/YAZrqzPVZT/eeDr+L2VK0M+P0kZoTpWsVsUZ+M/CO0PrAVuBFqBqfg9qR8GzozTyz1QfRc4XVXHAFOAm+KOZaRXJjg9JwPXJDg9sm6rajdwGf73XC57b/x7St2bDOYoG4ZhGEYRORd4RlVvTHuhiEwAPgF8SVXbVfVV4NvAJOCdgVOXqer97vglwMvA6RFlvh44Gfi6qr6oqi/iO5ynicikwKk/VdX1qrob32E5TETGplA/8nr3mSuBzzidBN8Butwdv15VN7meyD/hO/EnVpS/VFX/5M55ufLDE5bR33sss0RVH1bV1/AfAt4EXKSqO1T1ScDDD0GI02sXIMAUERmtqh2qeneCY/3Wq4wbbXhRRGa5/UNE5LcpbHGAe98SKPPDItLherFfcbIkdXsZcKSIlPfPAPYCfplCnx7MUTYMwzCMAuGGkecBcyOOt7lh8ZdcCEElh7r3B52j0QG8gO8sHBQ4b3N5Q1UVeBKImjBYvm5TQLax4hjA1sD2Dvc+JqLMMOKuvwL4hOsZPwEYB9wA4MID7hOR50VkO/Bl9jhgZTZX+/CEZfT3HsPK6QJ2q+pzFbIxcXqp6uNAG34v6tPih9OcFHcsC70CzAd+H3vH0Wxz7z31T1VvUtVxwAeBkU4cW7dV9WngV/jhFrj3q8MejJJgjrJhGIZhFIt34ztB60RkG/CAkz8oIl9S1eVuaH+0qo4Ouf4J9364qo4LvJpUNTi0fUh5w/XOHow/xA/QXVHmXyuvwe9pDB4bDG4HdgIfwg+7+IWqviwiBwFXA98BJqnqWOBH+L2pQSrvq4cUZVQjsvxaSaKXqt6gqh/AD4e5FlgpIk1xxzLS7zDgdcD9/SjmMeBx4KyY85LW7SXAmSLydvywkZrCLsAcZcMwDMMoGtfix2se7V6nOflJQGjqqyCq+iz+kPmPRaQFQETGicgZIhJ0rD8rIse43tlvAk3Ar92xZ4DDAmU+DdwGLHJl7Q8sAm5W1WAP5IDiwh2uwo8R/ggu7AIYje/TPAe8JiLT6TuhK44syuhlt4yoqpeIvFlETnHO72vAdkCB7mrHMtTvn/DDH2rGjWh8GfikiFwiIgeJTxOBcKEUdftW/F7q64E/qOq6WnUzR9kwDMMwCoSqdqnqU+UXvvMFfsxyWKhFGJ8HHgU8EekEHgJm4jtJZZYAPwRexJ+c9UFV3e6OLQaOdcPb5UUhZgGdrtxH8CcZfqqWe+wnV+Cn+tqkqvcCqOrDwD/ixzB3ABeQbGJYD1mUQbjd+kUCvfYGLsYPmejAf4j4qKq+EnOs34if8eN5Vd0Ye3IMqnoL/mjKEfijKC8B64Hj8cNsysTWbTepbyl+qMaS/uglvhNvGIZhGMZQQUQ2A/NV9eq8dTHqFxH5Kn6avZfxe9J34GdT2Qpcpqrvj7juUfwJeJtV9W2DpG7w82/Gd8CHRYQv9TAQefAMwzAMwzCMBkdVf4g/KoGILAA2qOofROSQmOviUucNKKqaaFU+MEfZMAzDMAzD6CequiCwvRkI7U2uNyz0wjAMwzAMwzBCsMl8hmEYhmEYhhGCOcqGYRiGYRiGEYI5yoZhGIZhGIYRgjnKhmEYhmEYhhGCOcqGYRiGYRiGEYI5yoZhGIZhGIYRgjnKhmEYhmEYhhHC/wNEjtGtDaKcXQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "data_x,_ = np.histogram(ak.to_numpy(all_data['Data']['mass']),\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(all_data[r'Signal ($m_H$ = 125 GeV)']['mass']) # histogram the signal\n", "signal_weights = ak.to_numpy(all_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(all_data[s]['mass']) ) # append to the list of Monte Carlo histogram entries\n", " mc_weights.append( ak.to_numpy(all_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", "# Main plot\n", "# *************\n", "fig, main_axes = plt.subplots(figsize=(12, 8))\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", "signal_heights = 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)*2.0 )\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.1, # 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=16 )\n", "\n", "# Add text 'for education' on plot\n", "plt.text(0.1, # 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=12 )\n", "\n", "# Add energy and luminosity\n", "lumi_used = str(lumi*fraction) # luminosity to write on the plot\n", "plt.text(0.1, # x\n", " 0.82, # y\n", " r'$\\sqrt{s}$=13 TeV,$\\int$L dt = '+lumi_used+' fb$^{-1}$', # text\n", " transform=main_axes.transAxes,fontsize=16 ) # coordinate system used is that of main_axes\n", "\n", "# Add a label for the analysis carried out\n", "plt.text(0.1, # x\n", " 0.76, # y\n", " r'$H \\rightarrow ZZ^* \\rightarrow 4\\ell$', # text\n", " transform=main_axes.transAxes,fontsize=16 ) # coordinate system used is that of main_axes\n", "\n", "# draw the legend\n", "my_legend = main_axes.legend( frameon=False, fontsize=16 ) # no box around the legend" ] }, { "cell_type": "markdown", "metadata": { "id": "pSP1-lqvk1Bt" }, "source": [ "### Signal Significance" ] }, { "cell_type": "markdown", "metadata": { "id": "L6bpFz1rk1Bt" }, "source": [ "We can do some analysis to study how significant the signal is compared to the background.\n", "One method is to check a quantity known as the signal significance $S$,\n", " which is defined by\n", "$$ S = \\frac{N_\\text{sig}}{\\sqrt{N_\\text{bg}}} $$\n", "where $ N_\\text{sig} $ and $N_\\text{bg}$ are the number of signal and background points respectively.\n", "A larger $S$ represents a better signal-to-background ratio,\n", " and a more significant signal peak.\n", "To calculate $N_\\text{sig}$,\n", " we can look at the plot and sum over the number of events of our Monte-Carlo signal.\n", "The signal range roughly corresponds to the bins from $115 \\,\\text{GeV}$ to $130 \\, \\text{GeV}$.\n", "$N_\\text{bg}$ then corresponds to the number of background events in those same bins." ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "FXOZx4clk1Bt", "outputId": "9de79f39-467b-4205-ca5b-3615a25c104f" }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "11.094884981667025\n", "[12.69896722 11.09488498 4.93528783]\n", "\n", "Results:\n", "N_sig = np.float64(28.729140024845144)\n", "N_bg = np.float64(10.300928352079978)\n", "signal_significance = np.float64(4.425965281174207)\n", "\n" ] } ], "source": [ "# Signal stacked height\n", "signal_tot = signal_heights[0] + mc_x_tot\n", "\n", "# Peak of signal\n", "print(signal_tot[18])\n", "\n", "# Neighbouring bins\n", "print(signal_tot[17:20])\n", "\n", "# Signal and background events\n", "N_sig = signal_tot[17:20].sum()\n", "N_bg = mc_x_tot[17:20].sum()\n", "\n", "# Signal significance calculation\n", "signal_significance = N_sig/np.sqrt(N_bg + 0.3 * N_bg**2) # EXPLAIN THE 0.3\n", "print(f\"\\nResults:\\n{N_sig = }\\n{N_bg = }\\n{signal_significance = }\\n\")" ] }, { "cell_type": "markdown", "metadata": { "id": "feedback" }, "source": [ "
\n", "We welcome your feedback on this notebook or any of our other materials! Please fill out this survey to let us know how we're doing, and you can enter a raffle to win some ATLAS merchandise!\n", "
" ] } ], "metadata": { "colab": { "provenance": [] }, "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.6" }, "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": 1 }