{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rJtApcSrSh36"
   },
   "source": [
    "<CENTER><img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/images/ATLASOD.gif?raw=1\" style=\"width:50%\"></CENTER>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KJu173g1Sh38"
   },
   "source": [
    "# How to rediscover the Higgs boson yourself!\n",
    "This notebook uses ATLAS Open Data http://opendata.atlas.cern to show you the steps to rediscover the Higgs boson yourself!\n",
    "\n",
    "ATLAS Open Data provides open access to proton-proton collision data at the LHC for educational purposes. ATLAS Open Data resources are ideal for high-school, undergraduate and postgraduate students.\n",
    "\n",
    "Notebooks are web applications that allow you to create and share documents that can contain for example:\n",
    "1. live code\n",
    "2. visualisations\n",
    "3. narrative text"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "HiEwRXmPSh39"
   },
   "source": [
    "### What is the Higgs boson?\n",
    "The Higgs boson is a fundamental particle predicted by the Standard Model. It is a manifestation of the Higgs field, which gives mass to fundamental particles. However, it is incredibly difficult to produce. At the LHC, a Higgs particle is produced about once every 10 billion collisions, making it very challenging to detect.\n",
    "\n",
    "Despite this tiny fraction, years of data collection led to the discovery of the Higgs boson in 2012 by the CMS and ATLAS experiments at CERN. In this tutorial, we will follow their example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fImwDqm8Sh39"
   },
   "source": [
    "### Detecting the Higgs\n",
    "This analysis loosely follows the [discovery of the Higgs boson by ATLAS](https://www.sciencedirect.com/science/article/pii/S037026931200857X) (Section 5) and one of the subsequent ATLAS measurements making use of more data and allowing for more detailed studies: [Measurements of Higgs boson properties in the diphoton decay channel](https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.052005) (Section 5).\n",
    "\n",
    "The Higgs boson can be generated through various mechanisms. In particle physics, we use Feynman diagrams to illustrate these production modes. These diagrams help us visualize particle interactions and serve as essential tools for computations. For additional details on Feynman diagrams, see this [link](https://cds.cern.ch/record/2759490/files/Feynman%20Diagrams%20-%20ATLAS%20Cheat%20Sheet.pdf).\n",
    "\n",
    "There are four primary production modes for the Higgs boson, each represented by its own Feynman diagram:\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",
    "<CENTER><img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/ImagesHiggs/ggH.png?raw=1\" style=\"width:40%\"> <img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/ImagesHiggs/VBFH.png?raw=1\" style=\"width:35%\"></CENTER>\n",
    "<CENTER><img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/ImagesHiggs/WH.png?raw=1\" style=\"width:40%\"> <img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/ImagesHiggs/ttbarfusion.png?raw=1\" style=\"width:35%\"></CENTER>\n",
    "\n",
    "\n",
    "The Higgs boson has an extremely brief lifetime, approximately $10^{-22} \\,\\text{s}$. It decays almost immediately after it is produced, making direct detection of the particle impossible. However, by using the Standard Model, we can predict the various decay products of the Higgs, such as photons, Z bosons, quarks, and others, each occurring with different probabilities. These **decay channels** help us identify the presence of the Higgs boson. In this notebook, we will focus on analyzing one specific decay channel:\n",
    "$$H \\rightarrow \\gamma\\gamma$$\n",
    "\n",
    "<CENTER><img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/feynman_diagrams/Hyy_feynman.png?raw=1\" style=\"width:40%\"></CENTER>\n",
    "\n",
    "We refer to this as our desired **signal**. Ideally, we aim to identify collisions that produce two photons, which would indicate the presence of a Higgs boson. However, along with our signal, many photons detected do not originate from Higgs boson decay but rather from other processes, forming the **background**.\n",
    "\n",
    "Backgrounds are classified into two categories: reducible and irreducible. **Reducible backgrounds** can be significantly minimized using experimental techniques such as data cuts, particle identification, and isolation criteria. For instance, in our case, a reducible background might involve events where a jet is misidentified as a photon. By applying stricter criteria to ensure that the detected particles are indeed photons (and not misidentified jets), this background can be reduced.\n",
    "\n",
    "On the other hand, irreducible backgrounds cannot be easily distinguished from the signal because they involve the same final states or processes that the signal would produce. In the scenario of Higgs decay into two photons, an **irreducible background** would be the direct production of two photons from other Standard Model processes, such as quark-antiquark annihilation. These events are fundamentally indistinguishable from the signal events based on final state particles alone.\n",
    "\n",
    "    \n",
    "To address this, we can consider the total invariant mass of the photon products. By conservation of energy and momentum, the invariant mass of the products must equal the Higgs mass, whereas other background processes will have different invariant masses. The final step is to plot the invariant mass of each event and identify the peak around 125 GeV, which corresponds to the mass of the Higgs 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": "TUWvQTHFSh3_"
   },
   "source": [
    "### Running a Jupyter notebook\n",
    "A Jupyter 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": "y7CtkAXTSh3_"
   },
   "source": [
    "## ATLAS Open Data Initialisation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "tvmgp1E4Sh4A"
   },
   "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",
    "If this is opened on mybinder, you don't need to run this cell."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ROgtJAUuSh4A"
   },
   "outputs": [],
   "source": [
    "#install required packages\n",
    "import sys\n",
    "%pip install atlasopenmagic\n",
    "from atlasopenmagic import install_from_environment\n",
    "install_from_environment()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FGqUO7iQSh4B"
   },
   "source": [
    "### To setup everytime\n",
    "We're going to be using a number of tools to help us:\n",
    "* uproot: lets us read .root files typically used in particle physics into data formats used in python\n",
    "* awkward: lets us handle complex and nested data structures efficiently\n",
    "* numpy: provides numerical calculations such as histogramming\n",
    "* matplotlib: common tool for making plots, figures, images, visualisations\n",
    "* lmfit: tool for statistical fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "id": "PXnTKJOiSh4B",
    "tags": []
   },
   "outputs": [],
   "source": [
    "import uproot # for reading .root files\n",
    "import time # to measure time to analyse\n",
    "import math # for mathematical functions such as square root\n",
    "import awkward as ak # for handling complex and nested data structures efficiently\n",
    "import numpy as np # # for numerical calculations such as histogramming\n",
    "import matplotlib.pyplot as plt # for plotting\n",
    "from matplotlib.ticker import MaxNLocator,AutoMinorLocator # for minor ticks\n",
    "from lmfit.models import PolynomialModel, GaussianModel # for the signal and background fits\n",
    "import vector #to use vectors\n",
    "import requests # for HTTP access\n",
    "import aiohttp # HTTP client support"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-9MszouWSh4B"
   },
   "source": [
    "## Example 1: Reading data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "x0fmpWqmSh4B"
   },
   "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 import the module and load the Open Data release."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "A-b2rSCLSh4B",
    "outputId": "f1970dad-1fa5-4dad-9b3c-6b04a5825dbf",
    "tags": []
   },
   "outputs": [],
   "source": [
    "import atlasopenmagic as atom\n",
    "atom.available_releases()\n",
    "atom.set_release('2025e-13tev-beta')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_z-yxvIZSh4D"
   },
   "source": [
    "We would like to load the Open Data ntuples which already have been selected (skimmed) to contain at least two photons."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "id": "bgDOIEgZSh4D",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Select the skim to use for the analysis\n",
    "skim = \"GamGam\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "l9CcEXQUSh4E"
   },
   "source": [
    "The data is organized by the collection periods throughout the year. In this notebook, we will use the 2015 data from periods D, E, F, G, H, and J, as well as the 2016 data from periods A, B, C, D, E, F, G, K, and L."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "fuB8PSq3Sh4E",
    "outputId": "f9190538-efe5-49a5-f836-686f2adfa7b3",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Let's get the list of files to go through\n",
    "# Notice that we use \"cache\" so that the files are downloaded locally and not streamed\n",
    "files_list = atom.get_urls('data', skim, protocol='https', cache=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4KYHQIh2Sh4E"
   },
   "source": [
    "Let's try accessing `data15_periodD` in the ATLAS Open Data directory as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "ReYkBgwISh4F",
    "outputId": "ec0c89e1-76ab-494f-88f4-f92fa38c5235",
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data15_periodD = 'simplecache::https://opendata.cern.ch/eos/opendata/atlas/rucio/opendata/ODEO_FEB2025_v0_GamGam_data15_periodD.GamGam.root'\n"
     ]
    }
   ],
   "source": [
    "# We will use the first entry in the list for testing\n",
    "data15_periodD = files_list[0]\n",
    "print(f\"{data15_periodD = }\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "q-5C7M_OSh4F"
   },
   "source": [
    "Next, let's open the `data15_periodD` file to examine its contents. The file, known as a `tree`, contains multiple entries, each representing an event. For each event, a dictionary stores all relevant information with keys, such as the event number (`eventNumber`), the photon transverse momentum (`photon_pt`), and more.\n",
    "\n",
    "Details on the variables in the dictionary can be viewed [here](https://cds.cern.ch/record/2707171/files/ANA-OTRC-2019-01-PUB-updated.pdf) in Appendix A.\n",
    "\n",
    "More information on trees can be viewed [here](https://uproot.readthedocs.io/en/latest/uproot.behaviors.TTree.TTree.html) and [here](https://hsf-training.github.io/hsf-training-uproot-webpage/03-trees/index.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "rDPsBl7-Sh4G",
    "tags": []
   },
   "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",
    "# The number of entries in the tree can be viewed\n",
    "print(\"The number of entries in the tree are:\", tree.num_entries)\n",
    "\n",
    "# All the information stored in the tree can be viewed using the .keys() method.\n",
    "print(\"The information stored in the tree is:\", tree.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "f90oP_edSh4G"
   },
   "source": [
    "Perhaps we'd like to see the energies of the photons.\n",
    "We can access this from our tree using the key `photon_e`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "waywCz3GSh4G",
    "tags": []
   },
   "outputs": [],
   "source": [
    "tree[\"photon_e\"].arrays()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ncI7E3HkSh4G"
   },
   "source": [
    "Depending on the analysis, some variables are more revelant than others, providing directly or indirectly better discrimination of the signal with respect to the backgrounds than others.\n",
    "The most important variables for the $H \\rightarrow \\gamma\\gamma$ analysis can be stored in a list and retrieved later from the tree, and correspond to the following:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "id": "agJyXBCCSh4G",
    "tags": []
   },
   "outputs": [],
   "source": [
    "variables = [\"photon_pt\",\"photon_eta\",\"photon_phi\",\"photon_e\",\n",
    "                            \"photon_isTightID\",\"photon_ptcone20\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "mwAuz6P1Sh4G"
   },
   "source": [
    "Now that we understand how to access the information in the `data15_periodD` tree, we can begin our analysis. As mentioned in the introduction, there are two key steps to complete for each event entry:\n",
    "#### Cuts\n",
    "We need to account for photons selection rules in the event.\n",
    "Based on the [Higgs boson discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X) and the [Higgs boson decay to photons measurement paper](https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.052005), one can define some main selection criteria:\n",
    "- Photon candidates are required to pass **identification criteria** because there is a risk of misidentifying jets and other particles that can mimic photon signals as actual photons.\n",
    "\n",
    "- The leading (sub-leading) photon candidate is required to have **$E_t$ (transverse energy) > 50 GeV (30 GeV)** because background processes frequently generate photons with lower $E_t$. By setting these constraints, we can reduce the impact of background processes and thereby improve the signal-to-background ratio. Given that photons are massless particles, enforcing these requirements on transverse energy is effectively equivalent to enforcing them on transverse momentum.\n",
    "\n",
    "- Photon candidates are required to have a **calorimeter isolation**, consisting on the sum of the transverse energies of energy clusters in the calorimeter around a spatial cone centered around the photon, in order to make sure the photons detected are not originating from jets. Additional photon transverse energy relative to the diphoton mass isolation is also required.\n",
    "\n",
    "- Since the transition between the the barrel and end-cap of the calorimeter can introduce uncertainties in the energy measurements of particles this issue is resolved by **excluding the calorimeter barrel/end-cap transition** region 1.37 < |η| < 1.52.\n",
    "\n",
    "- The leading and sub-leading photon are required to satisfy a **$p_T$ relative to the invariant mass threshold**, $p_T/m_{\\gamma\\gamma}>0.35$. This keeps both photons sufficiently energetic relative to the system and symmetric. This reduces backgrounds events since they tend to be more asymmetric.\n",
    "\n",
    "We need to filter the data such that in each event, the criteria mentioned above are satisfied.\n",
    "\n",
    "#### Mass calculation\n",
    "The data to be plotted is the di-photon invariant mass, which can be calculated using the equation: $$m_{\\gamma\\gamma} = \\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 the photon 4-momenta: (`photon_pt`,`photon_eta`,`photon_phi`,`photon_e`).\n",
    "\n",
    "From this,\n",
    "    we can see why we chose those six important variables earlier.\n",
    "Let's try to perform this two-step analysis for one event in `data15_periodD`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "evPnMKSvSh4H",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# This selects the first entry of the tree\n",
    "entry = tree.arrays(library=\"ak\")[:1]\n",
    "\n",
    "# Cut on the photon reconstruction quality\n",
    "photon_reconstruction = entry['photon_isTightID']\n",
    "# isTightID==True means a photon identified as being well reconstructed, only the events which have True for both photons are kept\n",
    "photon_reconstruction_cut_bool = (photon_reconstruction[:, 0] == True) & (photon_reconstruction[:, 1]==True) # [:, i] selects the i-th lepton in each event\n",
    "print(\"The reconstruction quality of the two photons:\", photon_reconstruction[:,0], photon_reconstruction[:,1])\n",
    "print(f\"Keep the event based on reconstruction quality requirements?{photon_reconstruction_cut_bool}\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "# Cut on the transverse momentum\n",
    "photon_pt = entry['photon_pt']\n",
    "# Only the events where the leading photon has transverse momentum (pt) > 50 GeV and the sub-leading photon has pt > 30 GeV are kept\n",
    "# Since the two photons for each entry are ordered, the first photon is the leading one and the second one is the sub-leading one\n",
    "photon_pt_cut_bool = (photon_pt[:,0] > 50) & (photon_pt[:,1] > 30)\n",
    "print(\"The transverse momentum of the two photons:\", photon_pt[:,0], photon_pt[:,1])\n",
    "print(f\"Keep the event based on the transverse momentum requirements?{photon_pt_cut_bool}\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "# Cut on the calorimeter isolation\n",
    "photon_ptcone20 = entry['photon_ptcone20']\n",
    "# Only the events where the invidivual photon calorimeter isolation is less than 5.5% are kept\n",
    "photon_caloiso_cut_bool = (((photon_ptcone20[:,0]/photon_pt[:,0]) < 0.055) & ((photon_ptcone20[:,1]/photon_pt[:,1]) < 0.055))\n",
    "print(\"The calorimeter isolation of the two photons:\", (photon_ptcone20[:,0]/photon_pt[:,0]), (photon_ptcone20[:,1]/photon_pt[:,1]))\n",
    "print(f\"Keep the event based on the calorimeter isolation requirements?{photon_caloiso_cut_bool}\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "# Cut on the pseudorapidity in the barrel/end-cap transition region\n",
    "photon_eta = entry['photon_eta']\n",
    "# Only the events where modulus of photon_eta is outside the range 1.37 to 1.52 are kept\n",
    "condition_0 = (np.abs(photon_eta[:, 0]) < 1.52) | (np.abs(photon_eta[:, 0]) > 1.37)\n",
    "condition_1 = (np.abs(photon_eta[:, 1]) < 1.52) | (np.abs(photon_eta[:, 1]) > 1.37)\n",
    "photon_eta_cut_bool = (condition_0 & condition_1)\n",
    "print(\"The eta of the two photons:\", photon_eta[:,0], photon_eta[:,1])\n",
    "print(f\"Keep the event based on the eta requirements?{photon_eta_cut_bool}\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "# This calculates the invariant mass of the 2-photon state\n",
    "p4 = vector.zip({\"pt\": entry['photon_pt'], \"eta\": entry['photon_eta'], \"phi\": entry['photon_phi'], \"e\": entry['photon_e']})\n",
    "invariant_mass = (p4[:, 0] + p4[:, 1]).M # .M calculates the invariant mass\n",
    "print(f\"The invariant mass of the 2-photon state is: {invariant_mass} GeV\")\n",
    "print(\"-\"*80)\n",
    "\n",
    "# Cut on the pT relative to the invariant mass\n",
    "# Only the events where the invididual photon relative pT is larger than 35% of the invariant mass are kept\n",
    "photon_massiso_cut_bool = ((photon_pt[:,0]/invariant_mass) > 0.35) & ((photon_pt[:,1]/invariant_mass) > 0.35)\n",
    "print(\"The invariant mass based isolation of the two photons:\", (photon_pt[:,0]/invariant_mass), (photon_pt[:,1]/invariant_mass))\n",
    "print(f\"Keep the event based on the invariant mass based isolation requirements?{photon_massiso_cut_bool}\")\n",
    "print(\"-\"*80)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VEgKHwAkSh4H"
   },
   "source": [
    "Based on our analysis, this entry should be removed because the photons do not match all our requirements.\n",
    "We can turn these checks and calculations into a set of functions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "ZpEt22z5Sh4H",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Cut on the photon reconstruction quality\n",
    "def cut_photon_reconstruction(photon_isTightID):\n",
    "    # Only the events which have True for both photons are kept\n",
    "    return (photon_isTightID[:,0]==True) & (photon_isTightID[:,1]==True)\n",
    "\n",
    "# Cut on the transverse momentum\n",
    "def cut_photon_pt(photon_pt):\n",
    "# Only the events where photon_pt[0] > 50 GeV and photon_pt[1] > 30 GeV are kept\n",
    "    return (photon_pt[:,0] > 50) & (photon_pt[:,1] > 30)\n",
    "\n",
    "# Cut on the energy isolation\n",
    "def cut_isolation_pt(photon_ptcone20, photon_pt):\n",
    "# Only the events where the calorimeter isolation is less than 5.5% are kept\n",
    "    return ((photon_ptcone20[:,0]/photon_pt[:,0]) < 0.055) & ((photon_ptcone20[:,1]/photon_pt[:,1]) < 0.055)\n",
    "\n",
    "# Cut on the pseudorapidity in barrel/end-cap transition region\n",
    "def cut_photon_eta_transition(photon_eta):\n",
    "# Only the events where modulus of photon_eta is outside the range 1.37 to 1.52 are kept\n",
    "    condition_0 = (np.abs(photon_eta[:, 0]) < 1.52) | (np.abs(photon_eta[:, 0]) > 1.37)\n",
    "    condition_1 = (np.abs(photon_eta[:, 1]) < 1.52) | (np.abs(photon_eta[:, 1]) > 1.37)\n",
    "    return condition_0 & condition_1\n",
    "\n",
    "# This function calculates the invariant mass of the 2-photon state\n",
    "def calc_mass(photon_pt, photon_eta, photon_phi, photon_e):\n",
    "    p4 = vector.zip({\"pt\": photon_pt, \"eta\": photon_eta, \"phi\": photon_phi, \"e\": photon_e})\n",
    "    invariant_mass = (p4[:, 0] + p4[:, 1]).M # .M calculates the invariant mass\n",
    "    return invariant_mass\n",
    "\n",
    "# Cut on null diphoton invariant mass\n",
    "def cut_mass(invariant_mass):\n",
    "    return (invariant_mass != 0)\n",
    "\n",
    "# Cut on the pT relative to the invariant mass\n",
    "# Only the events where the invididual photon pT is larger than 35% of the invariant mass are kept\n",
    "def cut_iso_mass(photon_pt, invariant_mass):\n",
    "    return ((photon_pt[:,0]/invariant_mass) > 0.35) & ((photon_pt[:,1]/invariant_mass) > 0.35)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c2oZITTaSh4H"
   },
   "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": 13,
   "metadata": {
    "id": "yYc1zwjVSh4H",
    "tags": []
   },
   "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\"):\n",
    "\n",
    "    photon_isTightID = data['photon_isTightID']\n",
    "    data = data[cut_photon_reconstruction(photon_isTightID)]\n",
    "\n",
    "    photon_pt = data['photon_pt']\n",
    "    data = data[cut_photon_pt(photon_pt)]\n",
    "\n",
    "    data = data[cut_isolation_pt(data['photon_ptcone20'],data['photon_pt'])]\n",
    "\n",
    "    photon_eta = data['photon_eta']\n",
    "    data = data[cut_photon_eta_transition(photon_eta)]\n",
    "\n",
    "    data['mass'] = calc_mass(data['photon_pt'], data['photon_eta'], data['photon_phi'], data['photon_e'])\n",
    "\n",
    "    data = data[cut_mass(data['mass'])]\n",
    "\n",
    "    data = data[cut_iso_mass(data['photon_pt'], data['mass'])]\n",
    "\n",
    "    # Append data to the whole sample data list\n",
    "    sample_data.append(data['mass'])\n",
    "\n",
    "# turn sample_data back into an awkward array\n",
    "data15_periodD = ak.concatenate(sample_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5lGA1m1MSh4H"
   },
   "source": [
    "We can now plot the data using Matplotlib.\n",
    "The data will be turned into a histogram,\n",
    "    with bins of width 3 GeV.\n",
    "Note that much of the code written here is meant for the aesthetics of the plot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 455
    },
    "id": "lQtxQqkpSh4H",
    "outputId": "95b1be55-2482-46e9-db64-47dcefe2d7a0",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjsAAAG2CAYAAACZEEfAAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAATwtJREFUeJzt3XlcVPX+P/DXcRBUFBRFARnDJbNJs43UCpckXK5pEWWZqaXesrE0Tc3uNe22mFam3dBu5ZJmK6KpkX3dwTTjYtwWl9Q0EcEthRgVYfj8/vDH5DgMs50zc86Z1/PxmEfOWT7nfd7nDOfd2T6SEEKAiIiISKfqBDoAIiIiIiWx2CEiIiJdY7FDREREusZih4iIiHSNxQ4RERHpGosdIiIi0jUWO0RERKRrLHaIiIhI10ICHYCcqqqqcOzYMTRq1AiSJAU6HCIiInKDEAJ//vkn4uLiUKeO/OdhPCp2rFYrDAaD7EHI5fDhw2jbtm2gwyAiIiIvHDx4EG3atJG9XY+KnZYtW2LEiBF47LHH0L59e9mD8VXdunUBALt370bLli29bicxMRG5ubkBm18NMZSWlsJoNKKgoAAREREBiUGONgIdgxx59DUGOeZXQxvMpXzz8/ctTxvcJ+Vro7CwECaTyXYcl5tHxY7ZbMaHH36I119/HbfddhtGjhyJBx54AA0aNFAkOE9VX7pq1KiRTzuewWAI6PxqiQEAIiIivG5HL3mQow1f8ihHDGrJA3OpnhgA/r7laoP7pO9tlJaWAoBit6B4dGFs2rRpOHDgADZu3Ig2bdpg7NixiI2NxejRo7Fz505FAgwEs9kc0PnVEoOv9JIHPeRSLXlgLtUTg6/0kgc95FIteVBDLp0SPvjzzz/F+++/L26//XYhSZIwmUzizTffdHv+V199Vdxyyy2iYcOGIjo6WgwaNEjs3bvXbpoePXoIAHafxx9/vMb2CgoKBABRUFDgy2qREKKkpEQAECUlJYEORdOYR/kwl/JhLuXBPMpH6eO3T7c8N2zYEKNGjcK2bduwZs0aFBcXY9KkSW7Pv3XrVpjNZnz33XdYv349KioqkJKSAovFYjfd6NGjUVRUZPvMnj27xvbCwsLs/kveCwsLw/Tp05lLHzGP8mEu5cNcyoN5lI/Sx29JCCG8nfncuXP4/PPPsXjxYmzbtg1t27bFY489hueee86r9k6ePInmzZtj69at6N69OwCgZ8+euOGGGzB37lyX85eWliIyMhIlJSWyXM8mIiIi5Sl9/PbqzM727dsxatQoxMbGwmw2IyEhAZs3b8avv/7qdaEDACUlJQCAqKgou+HLly9Hs2bN0LFjR0ydOhXnzp2rtZ3S0lK7T3l5udcxERERkbzKy8sdjtVK8ujMzuzZs7F48WL8+uuvuOWWWzBy5Eg89NBDaNSokc+BVFVVYeDAgTh79iy2bdtmG/7ee+/hqquuQlxcHH788UdMmTIFt956KzIzMx3aqK4MrzR9+nTMmDHD5xiJiIjIdzNmzMCLL77oMFypMzseFTvR0dEYOnQoRo4ciY4dO8oayJgxY/D1119j27ZtiI+Pdzrdpk2b0Lt3bxw4cMDhBYLVxc6V7zwICwvjNVUiIiKVKC8vt7vqUv3OIqWKHY/es3Ps2DGHF/5cuHAB9erV8ymIsWPHYu3atcjOzq610AGALl26AECNxU41X995QERERMrx90kIj+7ZqS50qqqq8NJLL6Fly5Zo2LAhfvvtNwCX3sOzcOFCt9sTQmDs2LFYuXIlNm3ahNatW7ucJz8/HwAQGxvrSehEREQBY7FYIEkSJElyeOKYlOfVDcovv/wylixZgtmzZyM0NNQ2vGPHjvjggw/cbsdsNuOjjz7Cxx9/jEaNGqG4uBjFxcU4f/48gEt9ZLz00kvIy8vD4cOHsXr1agwbNgzdu3fH9ddf703oREREFGS8KnaWLl2K9957Dw8//LBdx6CdO3fG3r173W5nwYIFKCkpQc+ePREbG2v7fPbZZwCA0NBQbNiwASkpKejQoQMmTpyI++67D2vWrPEmbCIiooCwWq22f2dnZ9t9J+V5VewUFhaiXbt2DsOrqqpQUVHhdjtCiBo/I0aMAAAYjUZs3boVp0+fxoULF7B//37Mnj2b9+MQEZFmZGZmwmQy2b73798fCQkJNT5VLJcRI0bYLpvVrVsXLVq0wF133YVFixahqqrK7XaWLFmCxo0bKxanv3hV7JhMJuTk5DgMz8jIwI033uhzUERERHqQmZmJtLQ0FBYW2g0vLCxEWlqaogVP3759UVRUhMOHD+Prr79Gr169MG7cOAwYMACVlZWKLVeNvCp2XnjhBYwdOxazZs1CVVUVMjMzMXr0aLzyyit44YUX5I6RiIhIc6xWK8aNG4ea3vBSPWz8+PGKXdIKCwtDTEwMWrZsiZtuugnPP/88vvzyS3z99ddYsmQJAGDOnDno1KkTwsPDYTQa8eSTT6KsrAwAsGXLFjz66KMoKSmxnSWqfmfdsmXLcMstt6BRo0aIiYnBkCFDcOLECUXWQw5eFTuDBg3CmjVrsGHDBoSHh+OFF17Anj17sGbNGtx1111yx0hERKQ5OTk5OHr0qNPxQggUFBTUeKVEKXfeeSc6d+5sO6NUp04dvP322/jll1/w4YcfYtOmTZg8eTIA4LbbbsPcuXMRERFh65vy2WefBQBUVFTgpZdewv/+9z+sWrUKhw8ftt2CokYevWfncklJSVi/fr2csRAREelGUVGRrNPJpUOHDvjxxx8BXDqzVC0hIQEvv/wynnjiCcyfPx+hoaGIjIyEJEmIiYmxa+Oxxx6z/btNmzZ4++23kZiYiLKyMjRs2NAv6+EJr4sdAPjzzz/tTs/VqVNHlStJRETkb+6+D87f740TQkCSJADAhg0bMHPmTOzduxelpaWorKzEhQsXcO7cOTRo0MBpG3l5eZgxYwb+97//4cyZM7abno8cOWJ3M7ZaeHQZKz8/H/3797d9j4uLQ5MmTWyfxo0bIzc3V/YgiYiItCYpKQnx8fG2wuJKkiTBaDQiKSnJr3Ht2bMHrVu3xuHDhzFgwABcf/31WLFiBfLy8pCeng4AuHjxotP5LRYL+vTpg4iICCxfvhy5ublYuXKly/kCyaNi59///jfuuOMOu2HLli3Dpk2bsHHjRgwZMgRvv/22rAESERFpkcFgwLx58wDAoeCp/j537ly799UpbdOmTfjpp59w3333IS8vD1VVVXjzzTfRtWtXtG/fHseOHbObPjQ01OEG6r179+L06dN47bXXkJSUhA4dOqj65mTAw2Jn+/bt6Nevn92wrl27okePHujZsyfMZjOys7NlDZCIiEirUlNTkZGRgbi4OLvh8fHxyMjIQGpqqmLLLi8vR3FxMQoLC7Fr1y68+uqrGDRoEAYMGIBhw4ahXbt2qKiowL///W/89ttvWLZsGd599127NhISElBWVoaNGzfi1KlTOHfuHFq1aoXQ0FDbfKtXr8ZLL72k2HrIQnigfv36oqCgwPZ9zpw5oqSkxPb9999/F2FhYZ40KauSkhIBQLRv315ce+214p133glYLERERNWqj08ARFZWlqisrFR0ecOHD7ctLyQkRERHR4vk5GSxaNEiYbVabdPNmTNHxMbGivr164s+ffqIpUuXCgDizJkztmmeeOIJ0bRpUwFATJ8+XQghxMcffywSEhJEWFiY6Natm1i9erUAIH744QeP4nznnXfEtddeK9q3by8A2NUUcpKEqOEFAE5ERUVhzZo1uP3222sc/+233+Luu+/GH3/84WsN5pXS0lJERkYq1kU8ERGRNywWi+0BnrKyMoSHhwc4InVR+vjt0WWsG2+8EatWrXI6PjMzk29QJiIiukJ4eLitSyQWOv7n0aPnTz75JB588EEkJCRgzJgxqFPnUq1ktVoxf/58/Pvf/8bHH3+sSKBERERE3vDoMhYATJkyBa+//joaNWqENm3aAAB+++03lJWVYcKECXj99dcVCdQdvIxFRESkPUofvz0udgDgu+++wyeffIL9+/cDAK6++mo89NBD6Nq1q+wBeoLFDhERkfYoffz26g3KXbt2DXhhQ0REROQOrzoCJSIiItIKFjtERESkayx2iIiISNdY7BAREZGusdghIiIiXZO12NmzZ4/t3TtEREREaiBrsXPx4kX8/vvvcjZJRERE5BOP3rMzYcKEWsefPHnSp2DkkpiYCIPBALPZDLPZHOhwiIiIqAbp6elIT0+H1WpVdDkevUHZYDDghhtucPp2w7KyMuzatUvxoJ3hG5SJiIi0R1VvUG7Xrh2eeeYZDB06tMbx+fn5uPnmm2UJjIiIiEgOHt2zc8sttyAvL8/peEmS4EVXW0RERESK8ejMzptvvony8nKn4zt37oyqqiqfgyIiIiKSi0fFTkxMjFJxEBERESmCLxUkIiIiXWOxQ0RERLrGYoeIiIh0jcUOERER6ZpHxU5SUhLeeOMN/Prrr0rFQ0RERCQrj4qd0aNHY8eOHbj55ptx7bXXYsqUKfj222/5bh0iIiJSLY+6i6hWXl6OjRs34ssvv8SaNWtgtVrxt7/9DQMHDkSfPn1Qv359JWJ1id1FEBERaY/Sx2+v7tkJCwtD//798Z///AfHjh3D6tWrERsbi2nTpqFp06YYMGAAvv32W7ljJSIiIvKYV2d2anPw4EGsXr0aRqMRaWlpcjbtEs/sEBERaY/Sx2/Zi51Aqk5W+/btYTAYYDabYTabAx0WERER1SA9PR3p6emwWq349ddfWey4g2d2iIiItEeV9+wQERERaQWLHSIiItI1WYodq9WK/Px8nDlzRo7miIiIiGTjVbEzfvx4LFy4EMClQqdHjx646aabYDQasWXLFjnjIyIiIvKJV8VORkYGOnfuDABYs2YNDh06hL179+KZZ57BP/7xD1kDJCIiIvKFV8XOqVOnEBMTAwDIysrC/fffj/bt2+Oxxx7DTz/9JGuARERERL7wqthp0aIFdu/eDavVinXr1uGuu+4CAJw7dw4Gg0HWAImIiIh8EeLNTI8++igeeOABxMbGQpIkJCcnAwB27tyJDh06yBogERERkS+8KnZmzJiBjh07oqCgAPfffz/CwsIAAAaDAc8995ysARIRERH5wqtiZ+nSpRg8eLCtyKn20EMP4dNPP5UlMCIiIiI5eNVdhMFgQFFREZo3b243/PTp02jevDmsVqtsAXqC3UUQERFpjyq7ixBCQJIkh+FHjx5FZGSkz0ERERERycWjy1g33ngjJEmCJEno3bs3QkL+mt1qteLQoUPo27ev7EF6KjExkb2eExERqdzlvZ4ryaPLWC+++KLtvxMnTkTDhg1t40JDQ5GQkID77rsPoaGh8kfqBl7GIiIi0h6lj98endmZPn06ACAhIQGDBw9GvXr1ZA+IiIiISE5ePY01fPhwAMDFixdx4sQJVFVV2Y1v1aqV75ERERERycCrYmf//v147LHHsH37drvh1TcuB+ppLCIiIqIreVXsjBgxAiEhIVi7dq3tLcpEREREauRVsZOfn4+8vDx2DUFERESq59V7dkwmE06dOiV3LERERESy86rYmTVrFiZPnowtW7bg9OnTKC0ttfsQERERqYVX3UXUqXOpRrryXp1A36DM9+wQERFpj6res1Nt8+bNsix85syZyMzMxN69e1G/fn3cdtttmDVrFq655hrbNBcuXMDEiRPx6aefory8HH369MH8+fPRokULWWIgIiIiffPqzI5c+vbtiwcffBCJiYmorKzE888/j59//hm7d+9GeHg4AGDMmDH46quvsGTJEkRGRmLs2LGoU6cOvv32W4f2eGaHiIhIe5Q+fntd7OTk5OA///kPfvvtN3zxxRdo2bIlli1bhtatW+OOO+7wKpiTJ0+iefPm2Lp1K7p3746SkhJER0fj448/RlpaGgBg7969uPbaa7Fjxw507drVbn4WO0RERNqjyl7PV6xYgT59+qB+/frYtWsXysvLAQAlJSV49dVXvQ6mpKQEABAVFQUAyMvLQ0VFBZKTk23TdOjQAa1atcKOHTu8Xg4REREFD6+KnZdffhnvvvsu3n//fdStW9c2/Pbbb8euXbu8CqSqqgrjx4/H7bffjo4dOwIAiouLERoaisaNG9tN26JFCxQXFztt68qnw6qLMSIiIgq88vJyvz7J7VWxs2/fPnTv3t1heGRkJM6ePetVIGazGT///DM+/fRTr+a/nNFoRGRkpO0zc+ZMn9tUisVigSRJkCQJFovF4/FERERaM3PmTLvjtNFoVHR5Xj2NFRMTgwMHDiAhIcFu+LZt29CmTRuP2xs7dizWrl2L7OxsxMfH2y3n4sWLOHv2rN3ZnePHjyMmJsZpewUFBXbX/MLCwjyOiYiIiJQxdepUTJgwwfa9tLRU0YLHqzM7o0ePxrhx47Bz505IkoRjx45h+fLlePbZZzFmzBi32xFCYOzYsVi5ciU2bdqE1q1b242/+eabUbduXWzcuNE2bN++fThy5Ai6devmtN2IiAi7D4sdIiIi9QgLC3M4VivJqzM7zz33HKqqqtC7d2+cO3cO3bt3R1hYGJ599lk89dRTbrdjNpvx8ccf48svv0SjRo1s9+FERkaifv36iIyMxMiRIzFhwgRERUUhIiICTz31FLp16+bwJBYRERFRTXx6z87Fixdx4MABlJWVwWQyoWHDhp4t3Elv6YsXL8aIESMA/PVSwU8++cTupYI1XcbS4qPnFovFlreysjLb+4XcHU9ERKR1qnzPzkcffYTU1FQ0aNBA9oB8wWKHiIhIe1T5np1nnnkGzZs3x5AhQ5CVlRWwvrCIiIiIXPGq2CkqKsKnn34KSZLwwAMPIDY2FmazGdu3b5c7PiIiIiKfeFXshISEYMCAAVi+fDlOnDiBt956C4cPH0avXr3Qtm1buWMkIiIi8ppXT2NdrkGDBujTpw/OnDmD33//HXv27JEjLiIiIiJZeHVmBwDOnTuH5cuXo3///mjZsiXmzp2Le++9F7/88ouc8RERERH5xKszOw8++CDWrl2LBg0a4IEHHsC0adNqfckfERERUaB4VewYDAZ8/vnn6NOnDwwGg9wxEREREcnGq2Jn+fLlcsdBREREpAiP7tnp378/SkpKbN9fe+01u17OT58+DZPJJFtw3kpMTITJZEJ6enqgQyEiIiIn0tPTYTKZkJiYqOhyPHqDssFgQFFREZo3bw7gUoeb+fn5tp7Ojx8/jri4uIC9ZJBvUCYiItIeVb1B+cq6yIdutYiIiIj8wutHz4mIiIi0wKNiR5Ikh57KnfVcTkRERKQGHj2NJYTAiBEjEBYWBgC4cOECnnjiCdt9JOXl5fJHSEREROQDj4qd4cOH230fOnSowzTDhg3zLSIiIiIiGXlU7CxevFipOIiIiIgUwRuUiYiISNdY7BAREZGusdghIiIiXWOxQ0RERLrGYoeIiIh0jcUOERER6RqLHSIiItI1XRY7iYmJMJlMSE9PD3QoRERE5ER6ejpMJhMSExMVXY4kdNR1udJdxCvBYrGgYcOGAICysjJb1xvujiciItI6pY/fujyzQ0RERFQtaIsdi8Vi68XdYrEEOhwiIiJSSNAWO0RERBQcWOwQERGRrrHYISIiIl1jsUNERES6xmKHiIiIdI3FDhEREekaix0iIiLSNRY7REREpGssdoiIiEjXWOwQERGRrumy2GGv50REROrHXs+94EmvqWrpTZy9nhMRUbBjr+dEREREPmCxQ0RERLrGYoeIiIh0jcUOERER6RqLHSIiItI1FjtERESkayx2iIiISNdY7ASY1Wq1/Ts7O9vuuzvjiYiIqHYsdgIoMzMTJpPJ9r1///5ISEhAZmamW+OJiIjINRY7AZKZmYm0tDQUFhbaDS8sLERaWhomT55c63gWPERERO5hdxHwfzcMVqsVCQkJOHr0qNNpDAaD00tWkiQhPj4ehw4dgsFgUCpMIiIiv2B3ETqUk5NTa6EDoNZ7c4QQKCgoQE5OjtyhERER6Y4uix2193peVFSkqnaIiIgCwV+9noco2nqA5ObmKnIaTC6xsbGqaoeIiCgQzGYzzGaz7TKWUnR5ZkftkpKSEB8fD0mSnE5jMBicjpckCUajEUlJSUqFSEREpBssdgLAYDBg3rx5AOBQ0EiSBEmSMGHCBKfjAWDu3Lm8OZmIiMgNLHYCJDU1FRkZGYiLi7MbHh8fj4yMDMyePbvW8ampqf4Ml4iISLP46Dn8/+j55S6/TpmVlYWUlBS7MzauxhMREWkdHz3XucsLl+7duzsUMq7GExERUe1Y7BAREZGusdghIiIiXWOxQ0RERLrGYoeIiIh0LaDFTnZ2Nu6++27ExcVBkiSsWrXKbvyIESNs752p/vTt2zcwwRIREZEmBbTYsVgs6Ny5c619WPXt2xdFRUW2zyeffOLHCImIiEjrAto3Vr9+/dCvX79apwkLC0NMTIyfIiIiIiK9UX1HoFu2bEHz5s3RpEkT3HnnnXj55ZfRtGnTWucpLS21+x4WFoawsDAlwyQiIiI3lZeXo7y83Pb9yuO23FR9g3Lfvn2xdOlSbNy4EbNmzcLWrVvRr18/WK3WWuczGo2IjIy0fWbOnOmniImIiMiVmTNn2h2njUajostT9ZmdBx980PbvTp064frrr0fbtm2xZcsW9O7d2+l8BQUFdq+b5lkdIiIi9Zg6daqtw2vg0pkdJQseVRc7V2rTpg2aNWuGAwcO1FrsREREKNK3BhEREfnO37eXqPoy1pWOHj2K06dPIzY2NtChEBERkUYE9MxOWVkZDhw4YPt+6NAh5OfnIyoqClFRUXjxxRdx3333ISYmBgcPHsTkyZPRrl079OnTJ4BRExERkZYEtNj573//i169etm+V1+/Gz58OBYsWIAff/wRH374Ic6ePYu4uDikpKTgpZde4j04RERE5LaAFjs9e/aEEMLp+G+++caP0RAREZEeaeqeHSIiIiJPsdghIiIiXQvaYufyFxNmZ2e7fFEhERERaVNQFjuZmZkwmUy27/3790dCQgIyMzM9asdisdh6Y7dYLB6Pp78wV0REpJSgK3YyMzORlpaGwsJCu+GFhYVIS0vzuOAhIiIidQuqYsdqtWLcuHE1PgFWPWz8+PG8pEVERKQjQVXs5OTk4OjRo07HCyFQUFCAnJwcP0ZFREREStJlsZOYmAiTyYT09HS74UVFRW7N7+50RERE5L309HSYTCYkJiYquhxNdQTqrtzc3Bo7AnW3Ty32vUVERKQ8s9kMs9mM0tJSREZGKrYcXZ7ZcSYpKQnx8fGQJKnG8ZIkwWg0Iikpyc+RERERkVKCqtgxGAyYN28eADgUPNXf586dC4PB4PfYiIiISBlBVewAQGpqKjIyMhAXF2c3PD4+HhkZGUhNTQ1QZERERKQEXd6z40pqaiqSk5Nt1wezsrKQkpLCMzpEREQ6FHRndqpdXth0796dhQ4REZFOBW2xQ0RERMGBxQ4RERHpGosdIiIi0jUWO0RERKRrLHaIiIhI11js1MJisUCSJEiSBIvFEuhwNEsvefTHeuglV0REasJih4iIiHSNxQ4RERHpmi6LncTERJhMJqSnpwc6FCIiInIiPT0dJpMJiYmJii5Hl91F5ObmIiIiItBhEBERUS3MZjPMZjNKS0ttXTgpQZdndoiIiIiqsdghIiIiXWOxQ0RERLrGYoeIiIh0jcUOERER6RqLHSIiItI1FjtERESkayx2iIiISNdY7BAREZGusdghIiIiXWOxQ0RERLrGYoeIiIh0TZfFTmJiIjp06ABJkiBJEiwWS6BDcio8PBxCCAghEB4e7vF4f7BYLJrIJflnW3F/ICK5+KvXc10WO7m5ucjLywt0GERERFQLs9mM3bt3Izc3V9Hl6LLYISIiIqrGYoeIiIh0jcUOERER6RqLHSIiItI1FjtERESkayx2iIiISNdY7BAREZGusdghIiIiXWOx4wOr1Wr7d3Z2tt13d8bTX5grIiJSCosdL2VmZsJkMtm+9+/fHwkJCcjMzHRrPP2FuSIiIiWx2PFCZmYm0tLSUFhYaDe8sLAQaWlpmDx5cq3jeRD/i6tcMldEROQrFjseslqtGDduHIQQDuOqO+ycM2eO0/EAMH78eF6mgetcAswVERH5jsWOh3JycnD06NFap6nt4CyEQEFBAXJycuQOTXNc5ZK5IiIiOYQEOgAlJCYmQpIkRdouKipSVTta5m4OmCsiIn1KT09Henq64mfwdXlmJzc3F3l5eYq0HRsbq6p2tMzdHDBXRET6ZDabsXv3buTm5iq6HF0WO0pKSkpCfHx8rWeODAaD0/GSJMFoNCIpKUmWeCwWCyRJgiRJsFgssrTprxhc5VLuXClNDdvCHVqJszZ6WAci8h8WOx4yGAyYN28eADgcpKv/+E6YMMHpeACYO3cuDAaDH6JVN1e5BJgrIiLyHYsdL6SmpiIjIwNxcXF2w+Pj45GRkYHZs2fXOj41NdWf4aqaq1wyV0RE5Ctd3qDsD6mpqUhOTkZkZCQAICsrCykpKbazEK7G01+YKyIiUhLP7Pjg8oNx9+7dHQ7OrsbTX5grIiJSCosdIiIi0jUWO0RERKRrLHZqwZ645aGXPPpjPfSSKyIiNQlosZOdnY27774bcXFxkCQJq1atshsvhMALL7yA2NhY1K9fH8nJydi/f78syw4PD7f1ZRUeHu4wnj1xy0MvefTHeuglV0REahPQYsdisaBz585IT0+vcfzs2bPx9ttv491338XOnTsRHh6OPn364MKFC4rGxZ645aGXPPpjPfSSKyIiVRIqAUCsXLnS9r2qqkrExMSI119/3Tbs7NmzIiwsTHzyySc1tlFSUiIAiJKSElFWViYACACirKzM7TgqKytFfHy8bd4rP5IkCaPRKCorK10uw9sYPKGGZdQ03pM8+ms9vCHn/iDHMlxxJwZfc62GbaWGGIhIPpcfv5Wg2nt2Dh06hOLiYiQnJ9uGRUZGokuXLtixY0et85aWlqK0tNT2vby83O3lsidueeglj/5YD73kiojIXeXl5bZj9ZXHbCWottgpLi4GALRo0cJueIsWLWzjnDEajXZv5H3jjTfcXi574paHXvLoj/XQS66IiNw1c+ZMREZG2j5Go1HR5am22PFFQUEBjh07Zvv+7LPPuj0ve+KWh17y6I/10EuuiIjcNXXqVJSUlNg+BQUFii5PtcVOTEwMAOD48eN2w48fP24b50xERAQiIiJs38PCwtxert564g4UveTRH+uhl1wREbkrLCzMdqy+8pitBNUWO61bt0ZMTAw2btxoG1ZaWoqdO3eiW7duii2XPXHLQy959Md66CVXRERqFdBip6ysDPn5+cjPzwdw6abk/Px8HDlyBJIkYfz48Xj55ZexevVq/PTTTxg2bBji4uJwzz33KBoXe+KWh17y6I/10EuuiIhUSZFnvNy0efPmGh+1HT58uBDi0uPn06ZNEy1atBBhYWGid+/eYt++fU7bk+PR85raAyCysrIcHv3lo+fujXeVR3+th6983R/8sQw5Hj1Xw37tihpiIH3hPhVYSj96HqJIBeWmnj17QgjhdLwkSfjXv/6Ff/3rX36M6i/siVseesmjP9ZDL7kiIlIT1d6zQ0RERCQHFjtERESka7otdth7tHyYS+1gz+xERI50WeysXr2avUfLhD1xawd7Ziciqpkui51HHnmEvUfLgD1xawd7Zicick6XxU5Nqp/6Gj9+PE+7u8FqtWLcuHE1Pi3HXKqLP7YV9wci0rKgKXYA9h7tCfbErR3smZ2IqHZBVexUY+/RrrEnbu1gz+xEpFXp6ekwmUxITExUdDlBWeyw92jX2BO3drBndiLSKrPZjN27dyM3N1fR5QRVscPeo93Hnri1gz2zExHVTrfFDnuP9g174tYO9sxORFQ7XRY7y5YtY+/RMmBP3NrBntmJiJzTZbEzcOBA7N692/Y9KysLhw4d4h9jL6SmpjKXGuGPbcX9gYi0SJfFDsDeo+XEXGoHe2YnInKk22KHiIiICGCxQzphsVggSRIkSYLFYgl0OEREpCIsdoiIiEjXWOxo3OV9EWVnZwekbyI1xKAGWsmDVuKsjR7WgYj8h8WOhmVmZsJkMtm+9+/fHwkJCX7tfVoNMaiBVvKglThro4d1ICI/EzpSUlIiAIiSkhJRVlYmAAgAoqyszKv2XLXh63hfrFixQkiSZGu/+iNJkpAkSaxYsUK2ZTlbD3djcCcPvuZKyVy7Woac20LJfU6u7eWPXDvjz/2egksg92uyP34rgWd2NMhqtWLcuHEQQjiMqx42fvx4RU/tqyEGNdBKHrQSZ230sA5EFBgsdjQoJycHR48edTpeCIGCggLk5OToOgY10EoetBJnbfSwDkQUGCGBDkAJiYmJTjss1IOioiJZp9NqDGqglTxoJc7a6GEdiMheeno60tPTFT8jq8szO7m5ucjLywt0GIqJjY2VdTqtxqAGWsmDVuKsjR7WgYjsmc1m7N69G7m5uYouR5fFjt4lJSUhPj7e6dkrSZJgNBqRlJSk6xjUQCt50EqctdHDOhBRYLDY0SCDwYB58+YBgMMf/urvc+fOVbTPIjXEoAZayYNW4qyNHtaBiAKDxY4PwsPDIYSAEALh4eEej/dFamoqMjIyEBcXZzc8Pj4eGRkZfumFWs4YlMyV0tSwLdyhlThro4d1IKIAUOSB9gDx93t21KB6nQGIrKwsUVlZKfsyXOXBVQyBfAeOP5chx7bwx7udfN1eavhd+GO/p+Cihv06mPE9O1Sry0/Zd+/ePSCn8NUQgxpoJQ9aibM2elgHIvIfFjtERESkayx2iIiISNdY7BBpiD96+2aP4hSMuN/rG4sdIo3wR2/f7FGcghH3e/1jsUPkRxaLBZIkQZIkWCwWt+fLzMxEWloaCgsL7YYXFhYiLS1Nlj/K/lhGNW/zIGf7SsegFlpYz0DG6M/93h+0sL0DgcUOkcr5o7dv9ihOwYj7ffBgsUOkcv7o7Zs9ilMw4n4fPNjrOZHK+aO3b/YoTsGI+33gsddzH+i913MKLv7o7Zs9ilMw4n4feOz1nIgA+Ke3b/YoTsGI+33wYLFTCy13TqkmesmjP9ajpmX4o7dv9ihOwYj7ffBgsUOkAf7o7Zs9ilMw4n4fHHR5gzKRHqWmpiI5ORmRkZEAgKysLKSkpMj6f53+WAaR2nC/1z+e2SHSEH/09s0exSkYcb/XNxY7REREpGssdoiIiEjXWOwQERGRrrHYIV24/O2b2dnZ7MuGiIhsWOyQ5mVmZsJkMtm+9+/fHwkJCZrrrZiIiJTBYoc0LTMzE2lpaSgsLLQbXlhYiLS0NBY8RETEYoe0y2q1Yty4cRBCOIyrHjZ+/Hhe0iIiCnIsdkizcnJycPToUafjhRAoKChATk6OH6MiIiK10WWxk5iYiMTERLzzzjua74/JHywWCyRJgiRJsFgsDuPV2rdVUVGRrNO5yoMWBKr/rsvJkUctbAs1xKiGGOgvrraHr+O1wpP1nDNnDkwmExITExWNSZfdReTm5iIiIiLQYZDCYmNjZZ2OiIj86/HHH8eECRNQWlpq665DCbo8s0PBISkpCfHx8Q69FVeTJAlGoxFJSUl+joyIiNSExQ5plsFgwLx58wDAoeCp/j537lz2cUNEFORY7JCmpaamIiMjA3FxcXbD4+PjkZGRgdTU1ABFRkREaqHLe3YouKSmpiI5Odl2vTcrKwspKSk8o0NERAB4Zod04vLCpnv37ix0iIjIhsUOERER6RqLHSIiItI1FjtEfsTe2S9ROg/utB8s20IL66mFGLWCuayZqoudGTNm2N6yWP3p0KFDoMMi8gp7Z79E6Ty4036wbAstrKcWYtQK5rIWQsWmT58urrvuOlFUVGT7nDx50un0JSUlAoAoKSnxY5SBVVZWJgAIAKKsrCwgbcgRg6/UnocVK1YISZJs46s/kiQJSZLEihUr/BKjO22oIQ/exuhO+2raFr7yxz6nJDXFqIbfli8xaCmXNY1X+vit+mKnc+fObk8fjMWOHNTwR9tXaih2nKmsrBTx8fEOf4Qu/2NkNBpFZWVlwGKUcxnO5vckD97E4E778fHxmtoW3sYg5z6nFLXFqOViR2u5DESxo+rLWACwf/9+xMXFoU2bNnj44Ydx5MgRl/OUlpbafcrLy/0QKVHN2Dv7JUrnwZ32jx49GhTbQgv7nBZi1Aot5/LyY7WSVF3sdOnSBUuWLMG6deuwYMECHDp0CElJSfjzzz9rnc9oNCIyMtL2mTlzpp8iJnIkd+/sWqV0HuTMn9a3hRb2OS3EqBVazmVcXBwiIyNhNBoVXY6q36Dcr18/27+vv/56dOnSBVdddRU+//xzjBw50ul8BQUFdr2eh4WFKRonUW3YO/slSudBzvxpfVtoYZ/TQoxaoeVcHjt2DOHh4SgtLVW04FH1mZ0rNW7cGO3bt8eBAwdqnS4iIsLuw2KHAom9s1+idB7caT8+Pj4otoUW9jktxKgVWs7l5cdqJWmq2CkrK8PBgwdVWZ0SOcPe2S9ROg/utD9v3ryg2BZa2Oe0EKNWMJduUOS2Z5lMnDhRbNmyRRw6dEh8++23Ijk5WTRr1kycOHGixun5NJYjd+7uV8NTJa744ykHpfOwYsUK0bJlS7unJIxGo0ePhGr5aaxq7uTB10dwXbUvRwxqeAJQjlz7g6tH5JXcH+SIUa42lN6ntLC9nY0P6kfPBw8eLGJjY0VoaKho2bKlGDx4sDhw4IDT6VnsOFLDQV4Oeih2hPhrHwUgsrKyPH4UVAt/9N2Z31UefI3BnTz7GoMWih0hfN/n5OAqTqX3BzlilKMNf+xTWtjegSh2VH2D8qeffhroEIhkxd7ZL1E6D+60HyzbQgvrqYUYtYK5rJmm7tkhIiIi8hSLHSIiItI1FjtERESkayx2CFar1fbv7Oxsu+/BRAt5CA8Ph7j0YAHCw8NVuQw58qiFbaGGGNUQg17I8dtytT18Ha8ValxPFjtBLjMzEyaTyfa9f//+SEhIQGZmZgCj8j/mQR5y5FEL20INMaohBvqLq+3h63it8HY9V69erWxgijzjFSB89NyRq/dbSJJk904G4FIPuZIk+f3dDLVR8pFOf+ZBC4/5e8uTPDrLg1zbQo5XDSgdo7tx1kSOXPuTGh7JVpKr7TFp0iSfxqvtXVzO+JKH6n8H5Xt2PMVix5GzHb+yslLEx8c77HSX73xGozEg72ioiVJ/LP2dB7X/0faWp3msKQ9ybgulih259xdv9gc5cu1vei52XG0PAMJgMHg93h/7lBzkyIOSx29exgpSOTk5OHr0qNPxQggUFBQgJyfHj1H5H/MgDznyqIVtoYYY1RAD/cXV9gDg8p6U2sZrZXvKkQcl6bLYSUxMhMlkQnp6eqBDUa2ioiJZp9Mq5kEecuRRC9tCDTGqIQb6i7/yrPbtqfb4VP0GZW/l5uYq3oOq1rnbmareO11lHuQhRx61sC3UEKMaYqC/+CvPat+eao9Pl2d2yLWkpCTEx8c79JBbTZIkGI1GJCUl+Tky/2Ie5CFHHrWwLdQQoxpioL+42h7ApS4cvB2vle0pRx6UxGInSBkMBsybNw8AHHa+6u9z587Vfb8qzIM85MijFraFGmJUQwz0F1fbQ5IkTJgwwevxgDa2p695UJwitz0HCJ/GcuTqzvwVK1aIli1b2t0NbzQaVfXYuRDKP83hrzyo+akSObibR1evAvB1Wyj56LlcMbobpzNy5Npf9Pw0VjVX28PX8e4KdK68Xc9ly5YpevxmsaNz7uz41XkDILKyslTzuPnl/PHH0h95CPQfIn9wJ4+u8uDrtlC62JEjRnfjrI0cufaHYCh2hHC9PXwd7w415Mqb9VT6+K3LG5TJM5efHu3evbvqT5fWpPpV777QQx7UQI48amFbqCFGNcRAf3G1PXwdrxVqXE/es0NERES6xmKHiEilLBaL7eZOi8Wiyhi0EKNWlsEYlMNiR+eCpRddpeen4OOPfSZY9kv+finQWOzoWLD0oqv0/BR8/LHPBMt+yd8vqYIitz0HCJ/G+ovWekV2xtceptWWBzXnWi7+eBJKyRjc3Wd8fWxcjmXIkWtXfH1E35f19Me2kCsP7kzjj6fSlM6FUnlQ+vjNMzs6ZLVaMW7cuBqfTqoeNn78eNWfCvZ1PfSSB/Iff+wzwbJf8vdLasJiR4f00iuyr+uhlzyQ//hjnwmW/ZK/X1ITXRY7wd7ruV56RfZ1PfSSB/Iff+wzwbJf8vdL7vjPf/4Dk8mExMRERZejy5cKBnuv53rpFdnX9dBLHsh//LHPBMt+yd8vuePxxx/HhAkTUFpaisjISMWWo8szO8FOL70i+7oeeskD+Y8/9plg2S/5+yU1YbGjQ3rpFdnX9VBjHqq7tRBCIDw83G/LVRul8+BO+zVN4499xt/7ZaD2OU/WM1DbQo/4N8YJRZ7xChA+em5PS70i18bX3oD1kgetkCOPgd4W7uwzvsYoxzLU8Mi1q/H++P2qIQ/uTKOHDlG1+ug5ix2d00qvyK742huwXvKgBXoodoRwvc/IEaOvy1DDQd6dGJT+/aohD+5Mw2LH+Xi+Z4d8Eiy96Co9PwUff+wzwbJf8vdLgcZih4iIiHSNxQ4RERHpGosdIiIi0jUWO0REKnV5v0/Z2dkB6QfKVQxaiFEry2AMymGxQ0SkQpmZmTCZTLbv/fv3R0JCAjIzM1UTgxZi1MoyGIPCFHnGK0D46LkjOR4TVANfY9RLHrRAL4+e++MxYGdtrFixQkiSZPduGQBCkiQhSZLtHTNK5slVDJMmTXIrRjni9DVPSuZBzvciBXI9XVE6D3zPjgdY7DjSy0GexY52sNjxbRmVlZUiPj7e4aBy+cHFaDSKyspKxfLkKgYAwmAwuBWjHLnyNU9K5cHT9fRmvD/WUw154Ht2iIiCSE5ODo4ePep0vBACBQUFyMnJCVgMAGq9V0MNMcoRgxa2RbDE4CtdFjuJiYkwmUxIT08PdCiawL5ULmEe5CFHHrWwLZSKsaioyO3p5IjBYrFAkiRIkgSLxeJRDO7EqBRP8gTUvJ5yL8PV9vBmvKcxuOIqD77sD97kYcmSJTCZTEhMTHRrGd4KUbT1AMnNzUVERESgwyAi8lhsbKys0ykZg7/a8aVtX2LQ0rbQagxmsxlmsxmlpaWIjIz0eH536fLMDhGRViUlJSE+Pt6hp+9qkiTBaDQiKSkpYDEAl7pwUHOMcsSghW0RLDH4isUOEZGKGAwGzJs3DwAcDi7V3+fOnato/1CuYpAkCRMmTFB1jHLEoIVtESwx+EyR254DhE9jOVLDUy1y8MfTWESXC3RP2itWrBAtW7a0e+rFaDTK/pixLzG4G6OSv19/xOCvbeGPGHx5YkzJPPDRcw+w2HGkl4M8ix3yt0AXO0L89TcNgMjKylLk8WJfY3AnRqV/v/6IwR/bwh8x+Pp4vFJ54KPnRERB6vLLAt27dw/IZQJXMWghRq0sgzEoh8UOERER6RqLHSIiItI1FjtEREFMq71YeypY1tMVLfRirwQWOzqnhTfRusPX9dBLHsh//LHPBHq/9Fcv1oH+/Wq6t24ZaaEXe6Ww2CEiCkKZmZlIS0tDYWGh3fDCwkKkpaXp4gAHBM96uuIqD5MnT9Z1nljsEBEFGavVinHjxkEI4TCuetj48eM1fwkjWNbTFVd5EEJgzpw5us4Tix0ioiCjh16s3REs6+mKFnqxV5ouix32ek5E5JzcPWmrVbCspytq7sU+PT2dvZ57i72eExE5p4aetP0hWNbTFTX3Ys9ez4mISBF66MXaHcGynq5ooRd7pbHYISIKMrroxdoNwbKermihF3ulsdghIgpCqampyMjIQFxcnN3w+Ph4ZGRkIDU1NUCRyStY1tMVV3mYPXu2rvMkiZqeNdOo6mt+JSUlvGeHiDTPYrGgYcOGAICysjJFXjx4+b0SWVlZSElJsfs/eFcx+CNGV9yJwdV6yrEMpckRg6s8+Lo/eEvp4zeLHSKiIKaGg7g/BMt6uqLW4lXp4zcvYxEREZGusdghIiIiXdNVsVNeXm73X/JeeXk5ZsyYwVz6iHmUjxZzabVasWXLFnzyySfYsmWLal63r8VcqhHzKB+lj9+aKHbS09ORkJCAevXqoUuXLvj+++9rnI7FjnzKy8vx4osvMpc+Yh7lo7VcZmZmIiEhAb169cKQIUPQq1cv1fQgrbVcqhXzKJ+gL3Y+++wzTJgwAdOnT8euXbvQuXNn9OnTBydOnFBsmb52MyFHNxVqiMFXesmDHnKpljwESy6re5i+sj+i6h6kR48erXgMSlPD/qCGGOSglzyoIZdOCZW79dZbhdlstn23Wq0iLi5OzJw502HagoICAUAUFBT4tMxrr702oPOrIYaSkhIBQJSUlAQsBjnaCHQMcuTR1xjkmF8NbWgll5WVlSI+Pl4AqPEjSZIICQkRlZWVisXgyuW5LCsrs8VWVlbmtxjkaMOT+Z2tZzDsk5dzlQdv9wdP47iSXMdvZ1TdN9bFixeRl5eHqVOn2obVqVMHycnJ2LFjh8P04v8/RX9lZ2VhYWEICwtze7lWqxWlpaVeRu37/GqIoXreYM+Dr23IkUdfY5BjfjW0oZVcutPTdmVlJdatW+f16/fl/H1f+Y4Vd+8rCvT+4On8FovF9u/L1zMY9snLucqDs/Fyx1FeXm53yar6uC0UehuOqt+zc+zYMbRs2RLbt29Ht27dbMMnT56MrVu3YufOnXbT//bbb2jbtq2/wyQiIiIZHDx4EG3atJG9XVWf2fFUQkICDh48iLp169r17+HpmR0iIiJSzpVndoQQqKioQEJCgiLLU3Wx06xZMxgMBhw/ftxu+PHjxxETE+MwfZ06dRSpCImIiEi7VP00VmhoKG6++WZs3LjRNqyqqgobN260u6xFRERE5Iyqz+wAwIQJEzB8+HDccsstuPXWWzF37lxYLBY8+uijgQ6NiIiINED1xc7gwYNx8uRJvPDCCyguLsYNN9yAdevWoUWLFoEOjYiIiDRA1ZexACA7OxvffPMNKioqcPHiRUydOhVdunSxjRdC4IUXXkBsbCzq16+P5ORk7N+/366NP/74Aw8//DAiIiLQuHFjjBw5EmVlZf5elYDLzs7G3Xffjbi4OEiShFWrVtmNz8zMREpKCpo2bQpJkpCfn+/QxoULF2A2m9G0aVM0bNgQ9913n8M9VcGgtlxWVFRgypQp6NSpE8LDwxEXF4dhw4bh2LFjdm1wv3S9T86YMQMdOnRAeHg4mjRpguTkZIenMJnHS1zl8nJPPPEEJEnC3Llz7YYzl5e4yuWIESMgSZLdp2/fvnbTMJfu7ZN79uzBwIEDERkZifDwcCQmJuLIkSO28XIdc1Rf7FgsFnTu3Nnpmxlnz56Nt99+G++++y527tyJ8PBw9OnTBxcuXLBN8/DDD+OXX37B+vXrsXbtWmRnZ+Pvf/+7v1ZBNVzl0mKx4I477sCsWbOctvHMM89gzZo1+OKLL7B161YcO3YMqampSoWsWrXl8ty5c9i1axemTZuGXbt2ITMzE/v27cPAgQPtpuN+6XqfbN++Pd555x389NNP2LZtGxISEpCSkoKTJ0/apmEeL3GVy2orV67Ed999h7i4OIdxzOUl7uSyb9++KCoqsn0++eQTu/HMpes8Hjx4EHfccQc6dOiALVu24Mcff8S0adNQr1492zSyHXMUeVWhQgCIlStX2r5XVVWJmJgY8frrr9uGnT17VoSFhYlPPvlECCHE7t27BQCRm5trm+brr78WkiSJwsJCv8WuNlfm8nKHDh0SAMQPP/xgN/zs2bOibt264osvvrAN27NnjwAgduzYoWC06lZbLqt9//33AoD4/fffhRDcL2viTh6r31i7YcMGIQTz6IyzXB49elS0bNlS/Pzzz+Kqq64Sb731lm0cc1mzmnI5fPhwMWjQIKfzMJeOasrj4MGDxdChQ53OI+cxR/Vndmpz6NAhFBcXIzk52TYsMjISXbp0sb1heceOHWjcuDFuueUW2zTJycmoU6eOw+lwql1eXh4qKirs8t2hQwe0atWqxjda019KSkogSRIaN24MgPulNy5evIj33nsPkZGR6Ny5MwDm0RNVVVV45JFHMGnSJFx33XUO45lLz2zZsgXNmzfHNddcgzFjxuD06dO2ccyla1VVVfjqq6/Qvn179OnTB82bN0eXLl3sLnXJeczRdLFTXFwMAA43K7do0cI2rri4GM2bN7cbHxISgqioKNs05J7i4mKEhobaDtjVLs83Obpw4QKmTJmChx56CBEREQC4X3pi7dq1aNiwIerVq4e33noL69evR7NmzQAwj56YNWsWQkJC8PTTT9c4nrl0X9++fbF06VJs3LgRs2bNwtatW9GvXz9b1wrMpWsnTpxAWVkZXnvtNfTt2xf/93//h3vvvRepqanYunUrAHmPOap/GotIyyoqKvDAAw9ACIEFCxYEOhxN6tWrF/Lz83Hq1Cm8//77eOCBB7Bz506Hgwk5l5eXh3nz5mHXrl12b5cn7zz44IO2f3fq1AnXX3892rZtiy1btqB3794BjEw7qqqqAACDBg3CM888AwC44YYbsH37drz77rvo0aOHrMvT9Jmd6rco1/aG5ZiYGJw4ccJufGVlJf74448a38JMzsXExODixYs4e/as3XBnb7QOdtWFzu+//47169fbzuoA3C89ER4ejnbt2qFr165YuHAhQkJCsHDhQgDMo7tycnJw4sQJtGrVCiEhIQgJCcHvv/+OiRMn2l7Pz1x6r02bNmjWrBkOHDgAgLl0R7NmzRASEgKTyWQ3/Nprr7U9jSXnMUfTxU7r1q0RExNj94bl0tJS7Ny50/aG5W7duuHs2bPIy8uzTbNp0yZUVVXZPcJOrt18882oW7euXb737duHI0eO8I3WV6gudPbv348NGzagadOmduO5X3qvqqrK1qcO8+ieRx55BD/++CPy8/Ntn7i4OEyaNAnffPMNAObSF0ePHsXp06cRGxsLgLl0R2hoKBITE7Fv3z674b/++iuuuuoqAPIec1R/GausrMxWLQOXbkrOz89HVFQUWrVqhfHjx+Pll1/G1VdfjdatW2PatGmIi4vDPffcA+BSldi3b1+MHj0a7777LioqKjB27Fg8+OCDNT56qWeucvnHH3/gyJEjtvfBVO+EMTExiImJQWRkJEaOHIkJEyYgKioKEREReOqpp9CtWzd07do1IOsUKLXlMjY2Fmlpadi1axfWrl0Lq9Vqu74cFRWF0NBQ7pf/X215bNq0KV555RUMHDgQsbGxOHXqFNLT01FYWIj7778fAH/fl3P1+76y4K5bty5iYmJwzTXXAGAuL1dbLqOiovDiiy/ivvvuQ0xMDA4ePIjJkyejXbt26NOnDwDmspqrfXLSpEkYPHgwunfvjl69emHdunVYs2YNtmzZAgDyHnM8enYrADZv3iwAOHyGDx8uhLj0+Pm0adNEixYtRFhYmOjdu7fYt2+fXRunT58WDz30kGjYsKGIiIgQjz76qPjzzz8DsDaB5SqXixcvrnH89OnTbW2cP39ePPnkk6JJkyaiQYMG4t577xVFRUWBWaEAqi2X1Y/u1/TZvHmzrQ3ul7Xn8fz58+Lee+8VcXFxIjQ0VMTGxoqBAweK77//3q4N5vESV7/vK1356LkQzGW12nJ57tw5kZKSIqKjo0XdunXFVVddJUaPHi2Ki4vt2mAu3dsnFy5cKNq1ayfq1asnOnfuLFatWmXXhlzHHEkIITwrj4iIiIi0Q9P37BARERG5wmKHiIiIdI3FDhEREekaix0iIiLSNRY7REREpGssdoiIiEjXWOwQERGRrrHYISIiIl1jsUNERES6xmKHSEMSEhLQs2dP2/ctW7ZAkiQsWbJE9mUtWbIEkiTZ+qnRIyXzpyStxk3KmDFjBiRJsn22bdsW6JBq9d1339nFO2PGDMWXyWKHiPwqPz8fM2bMwOHDhwMdCrmB20s73nrrLSxbtszWuevlLly4gPnz5+POO+9EdHQ06tati8aNGyMxMRFTpkzB3r17vVrm/fffD0mSkJ+f73QaIQRat26Nxo0b4/z582jXrh2WLVuGt956y6tlekP1vZ4TkXPdu3fH+fPnUbdu3UCH4rb8/Hy8+OKL6NmzJxISEgIaixbzB/g3bjVtL6rdPffcU+M2+u233zBgwADs2bMHPXr0wDPPPIPY2FiUlZUhPz8fixYtwhtvvIEjR46gZcuWHi1z5MiRyMjIwOLFizFv3rwap9m8eTMOHz6Mxx9/HPXr10f9+vUxdOhQHD58GM8884w3q+oxFjtEGlanTh3Uq1cv0GFoltbyZ7VaUV5ejgYNGmgqbgqc8+fP429/+xsOHjyIzMxM3HvvvQ7TXLhwAW+99RYkSfK4/ZSUFBiNRixfvhyvv/46QkNDHaZZvHgxgEuFUaDwMhaRChUUFOCBBx5AZGQkIiIicPfdd+PgwYMO03ly70b1PTgbNmzAjBkzcNVVVyEsLAzXX389Pv30U6fzVVVV4Y033kDbtm0RFhaG9u3b48MPP3SY7tSpUzCbzTAajQgNDYXRaITZbMbp06dt08yYMQOPPvooAKBXr162a/YjRozwqJ3L12fTpk1uxVeTmvLnSbtff/01JEnC22+/XWP73bp1Q3R0NCoqKgAAf/75J/75z3+iS5cuaNasGcLCwtCuXTs899xzOHfuXI3rt2HDBrz00kto27Yt6tWrh88//7zGuL1p29U6urO9alLd/saNG/Gvf/0LV111FerXr48uXbrgu+++AwBs3boVd9xxB8LDwxEbG4uXXnrJoR1P1unChQuYMWMGrrnmGjRo0ACNGzdGp06dMGnSJI+mUXqd/O2DDz7A3r17MWnSpBoLHQCoV68epk6diri4OLvh5eXlePXVV3HdddehXr16aNy4Me6++2788MMPtmnq1KmDESNG4PTp01i9erVD26WlpVixYgU6duyIxMREeVfOAzyzQ6QyZ8+eRffu3VFQUIAnnngCJpMJW7duRa9evXD+/Hmf258yZQosFguefPJJAJf+r+uhhx7ChQsXajyIPf/88zh//jwef/xxhIWFYcGCBRgxYgTatWuH22+/HQBQUlKC2267DQcOHMBjjz2Gm266CT/88AMWLFiATZs24fvvv0ejRo2QmpqKoqIivPfee3j++edx7bXXAgDatm3rUTuexucNd9pNSUlBTEwMli5diqefftpu/v379+O7777D008/bbvcVFhYiA8++AD33XcfhgwZgpCQEGzduhWzZ8/GDz/8gG+++cYhjmeffRYVFRUYPXo0IiIicM0116C8vNxhOm/adrWOrraXK8899xysVivGjRuHixcv4s0330RKSgqWLl2KkSNH4u9//zsefvhhfP7553jhhRfQunVrDB061Kt1MpvNWLRoEYYNG4YJEyagsrIS+/fvx6ZNmzyaRul18reMjAwAwKhRozyar6KiAn379sX27dvxyCOPYOzYsSgpKcH777+P22+/HdnZ2bjlllsAAI8++ihefvllLF68GGlpaXbtfPrppzh//nxAz+oAAAQRqcrUqVMFALFo0SK74ePGjRMARI8ePWzDNm/eLACIxYsXu2x38eLFAoBo1aqVOHv2rG342bNnRatWrUSTJk3EuXPnHKa/4YYbRHl5uW340aNHRWhoqHjwwQdtw55//nkBQKSnp9st85133hEAxD//+U+Hdjdv3uwQozftuBOfMzXlz9N2n332WQFA/PLLL3bD//nPfwoAIi8vzzasvLxcXLx40SGO6ml37tzpEEf79u2FxWJxGbc3bbuzjrVtL2eq57nxxhvt2v/yyy8FABESEiJyc3PtYo+JiRFdu3a1a8eTdWrSpIno169frXG5M43S66SE6dOnCwDi0KFDDuOioqJERESEw/DKykpx8uRJu8/lv/85c+YIAGLdunV285WUlAij0Wj3d0gIIe68805hMBjEsWPH7IZ37dpVhIaGipMnTzrEcOjQIQFATJ8+3f2V9RIvYxGpzKpVq9CiRQsMGzbMbviUKVNkaX/MmDGIjIy0fY+MjMQTTzyBM2fO1PiY+ZNPPml3Hb5ly5Zo37499u/fbxu2cuVKREdH4+9//7vdvI8//jiio6OxcuVKt2Lzph134vOGu+0OHz4cALB06VLbMCEEPvroI3Ts2BE33XSTbXhoaKjtLE9lZSXOnDmDU6dOITk5GQCwc+dOhzjGjBmDBg0auIzXm7aVyt3lsV/eflJSEgCgS5cutrMC1bHfeuutDsv1ZJ0iIyPxyy+/4Oeff3YajzvTKL1O/lZaWoqIiAiH4Xv27EF0dLTdJz093Tb+o48+QocOHXDzzTfj1KlTts/Fixdx1113Ydu2bXZnmkeOHAmr1Wr3O9i7dy++++47DBw4EM2aNVN2RV1gsUOkMr/99huuvvpqGAwGu+GxsbFo3LhxrfNarVYUFxfbfUpKSuymqb4UcTmTyWRb9pXatGnjMKxp06Z299AcOnQI11xzDUJC7K+Mh4SEoH379jW2WxNv2nEnPm+42251QbN8+XJUVVUBALKzs3H48GGHghUA5s+fj+uvvx5hYWGIiopCdHS07d1JZ86ccZi+ffv2bsfsadtK5c5Z+02aNAEAtG7d2mHaJk2a1Lhcd9dp7ty5OHPmDDp16oS2bdti1KhR+PLLL23bxN1p/LFO/hQREYHS0lKH4a1bt8b69euxfv16vPHGGw7j9+zZg7179zoURNHR0Vi0aBGsVitOnTplmz41NRWNGze23YwMAIsWLQIAPPbYYwqsmWd4zw6RjhQUFDj80R0+fLhPL5+7suiqJoTwuk05KRWfJ+0OGzYM48ePx6ZNm5CcnIylS5fCYDA43KsxZ84cTJw4ESkpKXj66acRFxeH0NBQFBYWYsSIETUedN05q+Nt20pvW2ftOxt+JU/WadCgQTh8+DCysrKwdetWbNiwAQsXLkRSUhI2bNiA0NBQt6ZRep38rWPHjsjOzsahQ4fs/jaEh4fbzpBd+T8XwKV9oFOnTpgzZ47TtqOjo23/rlevHoYMGYL58+dj+/bt6NKlC5YtW4b4+Hj06dNHxjXyDosdIpVp06YN9u/fD6vVavcHtKioCGfPnq113piYGKxfv95u2JVPWOzZsweDBg2yG7Z7927bsr2Ned++faisrLT7w1lZWYlff/3Vrt3aHm/1pB01GTJkCCZNmoSlS5fi9ttvR0ZGBu666y7ExsbaTbds2TIkJCTg66+/Rp06f51YX7dunc8xKNW2N48jy8XTdYqKisLQoUMxdOhQCCHw3HPPYfbs2fjyyy9x//33uz2NUm644QbExMRg7dq1tv07KysLo0aNQtOmTdGyZcsaxxUUFHhdTKWlpSE7OxsffPABXnnlFbfnu/rqq3Hy5EnceeeddrmvzciRIzF//nwsXrwYf/zxB4qLi/GPf/zD7fmVFPgIiMjOoEGDcPz4cbtr3wAwa9Ysl/PWq1cPycnJdp/qS1TVFixYYHdpq6SkBO+++y4aN26MHj16eBXzPffcg5MnT+KDDz6wG/7+++/j5MmTdo+8NmzYEADwxx9/+NSOmkRHR6Nfv37IzMzE8uXLUVpaaruX53IGgwGSJNmdOamsrMRrr73mcwxKtV3b9lKau+tktVod/kdAkiTceOONAC7F7s40StuxYwcKCgqwatUq27D3338fw4cPx/fff+90nC9njUaNGoUOHTrg9ddfd3rvnLOzlcXFxU7P7Bw/ftxh2E033YQbbrgBn332GdLT0yFJkiouYQE8s0OkOpMnT8bHH3+M0aNHIy8vD9dddx22bNmCHTt2yHKTX7NmzdClSxfb+1MWL16MI0eO4IMPPnD7kklNMX/xxRcwm83YtWsXbrzxRvzwww9YuHAhrrnmGkyePNk2bWJiIurUqYNXXnkFZ86cQXh4OFq3bo0uXbp41I7aDB8+HKtXr8bEiRMRGRmJe+65x2GatLQ0TJ06Ff369UNqaipKS0vx8ccfy/ImZKXarm17Kc3ddfrzzz8RGxuLgQMH4sYbb0Tz5s1x6NAhLFiwAE2aNMHdd9/t1jRKq1+/PkaOHInPPvsMaWlpKC4uxldffYXZs2fXOs7XZX711VcYMGAAUlNT0bNnT9srE0pLS7F371589tlnMBgMMBqNtvnGjRuH9evXY9KkSdi0aRPuvPNORERE4MiRI9i4cSPq1auHzZs3Oyxv5MiReOqpp7Bu3Tr07NlTNWdjWewQqUyTJk2Qk5ODCRMm2M7u9OjRA5s3b0bv3r19bn/WrFnIyclBeno6jh8/jvbt22P58uUYMmSI121GRkbi22+/xfTp07F69WosXrwYLVq0wBNPPIEXX3zR7t04rVq1wqJFizBr1iyMGTMGFRUVGD58OLp06eJRO2ozYMAAREVF4Y8//sCoUaNqfMPxpEmTIITAwoULMW7cOMTExGDw4MF49NFHHc7AeUqptmvbXkpzd50aNGiA8ePHY+PGjdiwYQPKyspshU31y/IuXrzochp/uP/++zFt2jSUlZVhyZIl6Nq1K66++mqX43zRpk0b5OXlYdGiRcjIyMCbb76JkpIShIeHo127dhg1ahRGjhxp16dW3bp18dVXX2H+/PlYtmwZpk+fDuDSZfFbb721xjOXAPDwww9j0qRJuHDhgmrO6gDge3aIgoU370shIvl17dpVLF++XFx99dViyZIlbo9zpvo9O7t27RInT54UFRUVSoQtm4qKCnHy5Emxa9cuvmeHiIhIj+6//348//zzKC4udrgpurZxrtx0002Ijo62dV+hVv/9738RHR1t9w4qpfEyFhERkR+lpaVh4sSJGD16tMN9crWNc2bYsGG44447bN87deoka7xyM5lMdk+N+uO+HhY7REREftSqVSsYjUbbe27cHedMmzZtVHMjsDsiIiI8Wj858DIWUZAYMWIEhBC2t88SUeDUqVPH6ftnahtH3mE2iYiISNckIVTyznciIiIiBfDMDhEREekaix0iIiLSNRY7REREpGssdoiIiEjXWOwQERGRrrHYISIiIl1jsUNERES6xmKHiIiIdI3FDhEREekaix0iIiLSNRY7REREpGssdoiIiEjX/h+oRbkzBws2dgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# x-axis range of the plot\n",
    "xmin = 100 #GeV\n",
    "xmax = 160 #GeV\n",
    "\n",
    "# Histogram bin setup\n",
    "step_size = 1 #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(data15_periodD),\n",
    "                        bins=bin_edges ) # histogram the data\n",
    "data_x_errors = np.sqrt( data_x ) # statistical error on the data\n",
    "\n",
    "\n",
    "# *************\n",
    "# Main plot\n",
    "# *************\n",
    "main_axes = plt.gca() # get current axes\n",
    "\n",
    "# plot the data points\n",
    "main_axes.errorbar(x=bin_centres, y=data_x, yerr=data_x_errors,\n",
    "                    fmt='ko', # 'k' means black and 'o' is for circles\n",
    "                    label='Data')\n",
    "\n",
    "# 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'di-photon invariant mass $\\mathrm{m_{\\gamma\\gamma}}$ [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": "_exurt_nSh4I"
   },
   "source": [
    "Great,\n",
    "    we managed to plot `data15_periodD`!\n",
    "Now that we have understood how to manage the data, apply the cuts and calculate the mass, we can proceed to analyze the data!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wWuQV5WGSh4I"
   },
   "source": [
    "## Final Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "AHO7QyURSh4I"
   },
   "source": [
    "For the final analysis, we'll begin by applying the cuts and calculating the invariant masses across all the data. Once that's done, we'll fit the data to uncover the Higgs boson peak. Let's kick things off by applying the cuts and calculating those invariant masses!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "E85MyMWQSh4I",
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Controls the fraction of all the events analysed. All of the data is used by\n",
    "# default to run this analysis (implemented in the loop over the tree). Reduce\n",
    "# this if you want the code to run quicker. This can take several minutes to\n",
    "# process the entire dataset\n",
    "fraction = 1\n",
    "# Holder for the masses as we process the files\n",
    "sample_data = []\n",
    "\n",
    "# Loop over all the files in our list\n",
    "for afile in files_list:\n",
    "\n",
    "    # Print which sample is being processed\n",
    "    print(f'Processing file {afile} ({files_list.index(afile)+1}/{len(files_list)})')\n",
    "\n",
    "    # Open file\n",
    "    tree = uproot.open(afile + \":analysis\")\n",
    "\n",
    "    numevents = tree.num_entries\n",
    "\n",
    "    # Perform the cuts for each data entry in the tree and calculate the invariant mass\n",
    "    for data in tree.iterate(variables, library=\"ak\", entry_stop=numevents*fraction):\n",
    "\n",
    "        photon_isTightID = data['photon_isTightID']\n",
    "        data = data[cut_photon_reconstruction(photon_isTightID)]\n",
    "\n",
    "        photon_pt = data['photon_pt']\n",
    "        data = data[cut_photon_pt(photon_pt)]\n",
    "\n",
    "        data = data[cut_isolation_pt(data['photon_ptcone20'],data['photon_pt'])]\n",
    "\n",
    "        photon_eta = data['photon_eta']\n",
    "        data = data[cut_photon_eta_transition(photon_eta)]\n",
    "\n",
    "        data['mass'] = calc_mass(data['photon_pt'], data['photon_eta'], data['photon_phi'], data['photon_e'])\n",
    "\n",
    "        data = data[cut_mass(data['mass'])]\n",
    "\n",
    "        data = data[cut_iso_mass(data['photon_pt'], data['mass'])]\n",
    "\n",
    "        # Append data to the whole sample data list\n",
    "        sample_data.append(data['mass'])\n",
    "\n",
    "print('Done processing data')\n",
    "\n",
    "# turns sample_data back into an awkward array\n",
    "all_data = ak.concatenate(sample_data)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ra4OaN9KSh4I"
   },
   "source": [
    "We are now ready to fit our data to effectively detect the Higgs boson! We will use a combination of a 4th order polynomial and a Gaussian function. The polynomial function represents the background, while the Gaussian function represents our signal. The Gaussian model is used to fit the signal due to the nature of the detector's resolution. The fourth-order polynomial is chosen for the background because it offers enough flexibility to capture the overall shape without overfitting, thereby reducing the influence of spurious data—random, irrelevant fluctuations or noise that do not correspond to the true signal or background."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 493
    },
    "id": "eByzAy0dSh4I",
    "outputId": "0e5a27bb-1a64-49cd-a74c-262a74e147fb",
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAoAAAAHcCAYAAABGaRsyAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjcsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvTLEjVAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAu3ZJREFUeJzs3XdYU9cbB/BvmGHIlKk4MUhEcKHiHijOUsDWWgfOaouto2prhxGtu9qqpW7FWX9VorWuaq2rah1QEBe4xYHWAchGeH9/RK7GhGkCAd7P8+SBe+7JuefeIL6cKSIiAmOMMcYYqzL0yrsCjDHGGGOsbHEAyBhjjDFWxXAAyBhjjDFWxXAAyBhjjDFWxXAAyBhjjDFWxXAAyBhjjDFWxXAAyBhjjDFWxRiUdwV0XV5eHu7fv49q1apBJBKVd3UYY4wxVgUQEZ4/fw5nZ2fo6Wm+vY4DwCLcunUL9evXL+9qMMYYY6wKun79OurVq6fxcosdAObm5kJfX1/jFdB1hoaGAIB69erBwMAAo0aNwkcffVSqsry9vXH27Nm3rlNlLSclJQUuLi5ISEiAhYWFTtSpMpaji89Zk2XpUjm6+Kwrazn8rMumHF18zposSxfKWblyJVatWoUXL17gxo0bQhyiacUOAGvUqIGhQ4di+PDhkEgkWqmMLsrv9j169Chq1qz5VmXp6+u/9T+YylxOPgsLi7cuT9fuTdfKAXTrOWuyLF0rB9CtZ11Zy8nHz1r75QC69Zw1WZYulDNp0iRMmjQJd+/ehYuLi9aGnxW7UzkkJATbt2+Hu7s72rdvj/DwcKSnp2ulUpVVSEgIl1NGdO3edK0cTdFkfXTtGVXWZ11Zy9EkXbs3XStHU/j3R/kSERGV5A1HjhzBunXrEBERAX19fbz//vsYOXIkWrVqpa06lqv8CDwhIeGtWwBZ4VJSUmBpaYnk5GSNtggwZfycyw4/67LDz7ps8HMuO9qOP0o8raRTp05Yv349EhMTsXDhQly+fBk+Pj5o1KgRFi1apPEKljdjY2Olr0x7jI2NIZPJ+FlrGT/nssPPuuzwsy4b/JzLjrbjjxK3AKqzZ88eDBkyBElJScjNzdVEvXQG/7XDGGOMsbKm7fij1AvLpKenIzw8HB07dsQ777wDW1tbzJo1S5N1Y4wxxhhjWlDidQBPnjyJtWvXYtu2bXjx4gX69euHmTNnokOHDtqoH2OMMcYY07BiB4Dz58/HunXrEB8fjxYtWmDBggUYMGAAqlWrps36McYYY4wxDSt2ALhgwQIMGjQI27Ztg4eHhzbrxBhjjDHGtKjYAeD9+/dVVqPOzMyEWCzWeKUYY4wxxpj2FHsSSH7wl5eXh5kzZ6JGjRowNzfHjRs3AADffvst1qxZo51aMsYYY4wxjSnxLODvvvsO4eHhmD9/PoyMjIR0Dw8PrF69WqOVY0xTnj17BhMTE4hEImzcuFHpnEgkKvbr1q1buHXrFkQiEcaOHVuiOnzxxRcQiURo0KBBofliYmIwYMAAuLq6QiwWo3r16vD09MTo0aPx77//Fvt6KSkpmDlzJpo1a4Zq1arB1NQUUqkUkydPxsOHD0tUd10QHh6u9FkYGhrC1tYWzZo1w5gxY3DixIm3vsaRI0cwffp0JCUlvX2FGWNMh5V4FvCGDRuwcuVKdO3aFWPGjBHSvby8cOXKFY1WjjFN2bx5M7KyslC3bl2sXbsWgwcPFs69GRAeP34cK1euxEcffYT27dsrnbOzs8N///1X4uu/ePECGzZsQP369XHt2jUcPXoUHTt2VMm3e/duvPvuu7Czs8OQIUPg6uqKpKQkxMfHY8+ePWjQoAGaNm1a5PXi4+Ph5+eH27dvIzAwECNGjIChoSH++ecfLF68GOvWrcPvv/8OHx+fEt9Lefvss8/g7e2NvLw8JCcn48KFC5DL5VixYgU+/PBDrFu3TumP05I4cuQIQkNDMXToUFhZWWm24owxpkuohMRiMd26dYuIiMzNzen69etERHTx4kUyMzMraXE6Lzk5mQCQRCIhd3d3+umnn8q7SqwUmjRpQl26dKHFixeTSCQSfm7VWbduHQGgdevWqT1/8+ZNAkAhISHFvv7OnTsJAB06dIjs7e1pyJAhavN5eHiQubk5JSQkqJzLzc2lR48eFXmttLQ0kkgkZGhoSLt371Y5f/bsWbK0tCQ7OztKTEws9j2Ut/zPZdu2bSrn0tPTacCAAQSAxowZU+pryGQyAkA3b958i5oyxljp/fTTT+Tu7k4SiYQAUHJyslauU+IuYKlUiuPHj6ukb9++vVgtExXV2bNncenSpQqxwTNTFhUVhejoaAQHB+PDDz+EgYEB1q5dW6Z1WLNmDerVq4fOnTtj4MCB2L59O1JSUlTyXb16FW5ubmr3fdTT04OdnV2xrhUfH4/x48ejd+/eKudbtGiB2bNn47///sOCBQuE9CNHjkAkEiE8PBxLly6FRCKBWCyGRCLB0qVL1V7r6tWrGDx4MJycnGBkZIQ6depg8uTJSEtLU8o3dOhQiEQiJCcn4+OPP4a9vT3EYjHatm2L06dPF3lPRTExMUF4eDjq1auHVatW4datW8K5K1eu4JNPPkGjRo2ErvDmzZurDFkZOnQoQkNDAQB169YVupqnT58OQDER7vPPP0eTJk1gbW0NsVgMqVSKefPmVbodkBhj5SckJASXLl3C2bNntXqdEgeA06ZNw9ixYzFv3jzk5eVBLpdj1KhRmDVrFqZNm6aNOjL2VtasWQNzc3MEBQWhevXq6NOnD9avX4+8vLwyuX5iYiL27duHIUOGQCQSYejQoUhPT8fWrVtV8tavXx8XL17EyZMnS3297du3AwA++uijAvMMHToUhoaGiIiIUDm3dOlSzJ07F4MGDcKcOXNgZWWFzz77TAiO8kVGRqJFixY4duwYRo8ejbCwMPTp0wdLlixBt27dkJOTo1K2n58f7t69i2nTpmHq1Km4cOECevfujefPn5f6fvMZGRlh8ODByM3NxR9//CGkHzlyBMeOHUOfPn2wYMECzJw5E4aGhhg1ahTmzJkj5Bs9ejQCAgIAAD/88AM2btyIjRs3IjAwEABw/vx5yOVydOnSBd999x3mzp2LWrVq4csvv8Qnn3zy1vVnjLEyVZpmw2PHjpGvry/Z2dmRiYkJtW3blv744w9Nt07qhPwuYG01wTLtysjIICsrKwoODhbS8rtj9+7dq/Y9mu4Cnjt3LolEIrpx44aQ1qRJE2rZsqVK3m3btpFIJCIA1LhxYxo9ejStWbOmRF2SNjY2VK1atSLzNW7cmADQ8+fPiYjo8OHDBEClCzorK4u8vb3JwMBAKd3T05Pc3NwoJSVFqVy5XK7y/IKDgwkAffzxx0p5f/31VwJAy5cvL7K+hXUB54uIiCAANHHiRCEtNTVVJV9ubi517NiRLCwsKDs7W0gvrAs4PT2d8vLyVNIHDRpEenp6dP/+/SLvgTHGikvb8Uep9gJu3749Dh48iEePHiE9PR1///03unfv/vbRKGMaJpfLkZSUhODgYCGtV69esLOzK7Nu4LVr16J9+/aoW7eukDZ06FCcOXMGFy9eVMrbr18/HDt2DP369UNCQgJWrFiBESNGoG7duvD39y/WBJT8DcSLkr+5eHJyslL6wIEDlbqgjYyMMGHCBLx48QK///47ACA2Nhbnz5/Hhx9+iKysLDx+/Fh4tWvXDmZmZjhw4IDKNSdMmKB03KVLFwCKrmRNyL+n17vXzczMhO8zMzPx5MkTPH36FN27d0dKSkqxJ6/lzyIHgOzsbDx9+hSPHz+Gn58f8vLycO7cOY3cA2OMlYUSzwLO9/z5cxCRcKynpwdzc3ONVIoxTVmzZg3s7OxQs2ZNXLt2TUjv3r07tm3bhsePH6N69epau/7x48cRHx+PQYMGKV2/VatW0NPTw5o1a7Bo0SKl97Rr1w7t2rUDEeHq1as4fPgwfv75Z+zatQuDBg1S6t5Ux8LCQu34wjfl53kzWHR3d1fJK5VKAUBY9/Py5csAAJlMBplMprZ8dUvN1KtXT+nY1tYWAPDkyZMi61sc+feUHwgCQGpqKqZPn45ff/0VCQkJKu959uxZscp+8eIF5s6diw0bNuDatWtKv/9KUg5jjOmCYgeA0dHR+Oqrr7B3714AgLOzM9LT04XzIpEIp06dgre3t+ZryVgp3Lx5E4cPHwYRQSKRqM2zadMmjB8/Xmt1yF8cfdq0aWrHyG7atAnz5s1T2WUHUPybkkgkkEgkCA4ORqNGjXDgwAHcvXtX7SSRfB4eHjh27BiuXbsGV1dXtXnS09Nx5coV1KlTp1R/uOUHP59//jl69OihNo+1tbVKmr6+fqHlva3z588DANzc3IS0Dz/8ELt378ZHH32EDh06wNbWFvr6+ti7dy9++OGHYo8FnThxIpYuXYr+/fvj66+/hr29PQwNDREVFYUvvviizMaUMsaYJhQ7AFy6dCnatWunlLZx40bUqFEDRIS1a9diyZIlKmuqMVZe1q1bByLCqlWr1K7p9s0332Dt2rVaCwCfP3+O7du3o1u3bmonZJw/fx4zZ87Erl27EBQUVGhZYrEYTZo0wY0bN3Dv3r1CA8DAwEAcO3YMq1evxty5c9Xm2bBhA3JycoQJDq/Lb9173aVLlwC8asHLX8xaX18fvr6+hda9rGRnZ2Pjxo3Q19eHn58fACApKQm7d+/G4MGDsXz5cqX8f/75p0oZ+V286mzcuBEdOnRQmbzzessuY4xVGMUdLNiwYUOKiooSjl9fA5CI6J9//qFatWppZmSiDuFJIBVTbm4uubi4UOPGjQvMM336dAJAZ86cUUrX1CSQlStXEgCKiIhQez4tLY1MTU2pV69eQtq+ffvUTjR49OgROTg4kIGBAT158qTQ66amppKrqysZGhrSvn37VM5HRkaSlZUV2dnZ0YMHD4T0oiaB6Ovr0507d4iIKC8vjzw8PMjCwkLtmoo5OTlK9cyfBKIOAKVJOgUp7jqAn3zyiZCekpJCIpGIBg8erJT//v37ZGtrSwDo8OHDQvr3339PACgyMlLlGra2ttS+fXultNTUVGGtroJ+XhhjrDS0HX8UuwXw9u3bSmuQzZgxQ2nslJOTU4XcXopVTgcOHEBCQgJGjBhRYJ6goCBMnz4da9asKdXQhXPnzuG7775TSTcwMMCXX36JNWvWwNTUtMAuUlNTU/Ts2RM7d+7EvXv3UKNGDfTr1w/29vbo06cPpFIpDAwMcOPGDWzcuBEPHz7EtGnTYGNjU2i9zMzMsGvXLvTo0QO9e/dGUFAQOnXqBAMDA5w5cwYbN26Eubk5du7cCUdHR5X3SyQStGrVCmPGjEG1atWwZcsWnD17Ft9++y1cXFwAQNhSr0uXLvD09MTw4cPRqFEjpKen49q1a5DL5ZgzZw6GDh1a4udalOPHjyMzMxNEpLQTyH///YdBgwbhxx9/FPJWq1YN3bt3x6ZNm2BiYgJvb2/cvn0bK1asQN26dVXGHrZu3RqAYtu+gQMHQiwWw8PDAx4eHujXrx9WrFiB/v37w9fXFw8fPsTatWuFcYyMMVahFDdStLa2pr///rvA83///TdZW1trIijVKdwCWDH169ePAND58+cLzSeRSMjS0pLS09OFtOK2ABb0MjY2pgsXLhAACgwMLPT6W7ZsIQA0a9YsIlIsizJs2DCSSqVkZWVFBgYGZG9vTz169KDt27eX6BkkJSVRaGgoeXl5kZmZGYnFYnJzc6PPP/9cqeUvX34L4Lp162jx4sXk6upKRkZG5OrqSj/++KPaa9y6dYtGjx5NtWvXJkNDQ7KxsaFmzZrRl19+KbQWEmm2BTD/pa+vT1ZWVtSkSRMaPXo0nThxQu37/vvvPxoxYgQ5OTmRsbExeXh40MqVK4XyXm8BJCKaN28e1a1blwwMDAgAyWQyIlK02E6aNIlq1apFxsbG5OrqSnPmzKE///yTWwAZYxqn7fhDRFS80dddu3ZFs2bNlHYOeN3nn3+O6OhoHDp0qJShqG7KX1IjOTlZaWYhY5XNkSNH0LlzZ6xbt04rLXeMMcaKT9vxR7G7gD/55BN88MEHqFOnDj7++GPo6SmWEMzNzcXPP/+MpUuXYsuWLRqvIGOMMcYY06xiLwQdFBSEiRMn4tNPP4W1tTWaNm2Kpk2bwsbGBuPHj8e4cePQr18/bda1XHl7e0MqlSIsLKy8q8IYY4yxSiosLAxSqVTry+oVuws43z///INffvlFWLm/QYMGGDBggDB4urLhLmBWVXAXMGOM6Q5txx8lDgCrGg4AGWOMMVbWtB1/lGovYMYYY4wxVnFxAMgYY4wxVsVwAMgYY4wxVsVwAMgYY4wxVsVwAMgYY4wxVsVoLAC8fPky6tWrp6niGGOMMcaYlmgsAMzOzsbt27c1VRxjGrF582Y0btwYzZo1Q3R0dJlcc/369QgODtZK2devX8fGjRuF47y8PHTo0AHp6elauR5jjLHKqdgB4MSJEwt9ff/994W+/9ixY+jbty+cnZ0hEomwc+dOpfNEhGnTpsHJyQkmJibw9fUVFpvO9/TpUwwcOBAWFhawsrLCiBEjkJqaqpTn/PnzaN++PcRiMVxcXDB//nyVumzbtg0NGzaEWCxG48aNsXfv3uI+BlaB5OXl4fPPP8fevXsRFRWFJk2aFOt9ubm5b3Xd06dPo0WLFm9VRkF+//13REZGCsd6eno4duwYTE1NtXI9xhhjlVOxA8DFixfj6NGj+Pfff9W+rly5Uuj709LS4OXlVeBWavPnz8eSJUuwfPlynD59GmZmZvDz80NmZqaQZ+DAgbh48SIOHjyI3bt349ixY/joo4+E8ykpKejevTtq166NyMhILFiwANOnT8fKlSuFPCdPnsSAAQMwYsQI/Pvvv3j33Xfx7rvv4sKFC8V9FKwCePbsGRo2bIjnz5+jb9++WLZsGW7fvo133nkHXl5ekEql2Lp1q5C/Ro0amDlzJnx8fDBmzBiV8q5cuQJ/f3+0atUK7u7uWL9+vXDu2LFjaNmyJTw9PfHVV1/h9OnTwhY+zs7OuHPnjpC3efPmQgB369YtvPfee2jevDkkEglmzJgBANizZw98fHzQrFkzuLm5YcOGDQCADRs2QCaTYdeuXWjSpAl27dqFVatWCf8G8vLyMH36dHh5eaFRo0YIDg5GTk4OAODrr7/GqFGjEBAQADc3N3Tq1IlbDRljrCqjYpJIJLRx48YCz//777+kp6dXrLIA0I4dO4TjvLw8cnR0pAULFghpSUlJZGxsTL/88gsREV26dIkA0NmzZ4U8+/btI5FIRPfu3SMiop9//pmsra0pKytLyPPFF1+Qm5ubcPz+++9T7969lerTqlUrGj16tNq6JicnEwBKTk4u1r0x3fHLL79Q//79iYgoIyODJBIJ/fHHH0RElJCQQBYWFpSSkkIPHjwgABQWFqa2nEePHlHr1q3p2rVrRET09OlTsre3p8zMTLp48SLVr1+frl69SkRECxcuJENDQ8rIyKD79++TnZ2dUE5WVhaZmZlRVlYWPX78mGrVqqX07+D+/ftERPTff/9Rbm4uERHdvXuXnJ2dhTyNGjWi2NhY4Xj06NG0fPlyIiKaOnUqDR8+nHJycigvL4/69OlDK1asICIiPz8/8vf3p7S0NCIiatGiBR06dKiUT5Yxxpi2aTv+KHYLYIsWLZS6nt4kEolApdxV7ubNm0hMTISvr6+QZmlpiVatWuHUqVMAgFOnTsHKykqpa83X1xd6eno4ffq0kKdDhw4wMjIS8vj5+SEuLg7Pnj0T8rx+nfw8+dcpSEpKitIrKyurVPfKys65c+fQvHlzAIqxgM2aNUP37t0BADVr1oSpqSmSkpIQGRkJb29vfPLJJ2rLWbp0KRISEhAUFIQmTZqgc+fOMDQ0RG5uLubNm4dx48bB1dUVANC4cWM0atQIYrEYkZGRwvUB4MKFC2jQoAGMjIywZMkS+Pv749133xXOOzk5gYiwefNmtGvXDl5eXvD19YVYLAagaEVPSEiAu7u78J7IyEg0a9YMaWlpWL58ORYtWgQDAwOIRCJ4eHjg4cOHAICoqCgsWLBA6CrOzs6GjY2Nhp40Y4yxt5WVlaUSa2hTsQPAhQsXYvz48QWe9/LyQl5eXqkqkZiYCABwcHBQSndwcBDOJSYmwt7eXum8gYEBbGxslPKoK+P1axSUJ/98QVxcXGBpaSm85syZU5JbZOXg9QAsNjYWzZo1E87dv38fubm5qFmzJqKiotCtW7cCy/n333+xcuVKREdHC6+7d+/C1NQUUVFRaNmypZD31KlTwh8p+cFZvhMnTgjHUVFR6Nixo8q1wsLC8Ndff2H37t2IiYnB+PHj0bRpU+E9jRs3hr6+PgAgJycHV65cgaenJ65evQonJydYWloKZUVFRaFp06ZISEiAvr4+GjRoAADIzMzE9evX0ahRo5I9UMYYY1ozZ84cpTjDxcVFq9crdgDo6OiI2rVra7MuOi0hIQHJycnCa+rUqeVdJVYIIsK///4rBFwuLi6IjY0FAGRkZGDMmDGYMmUKRCIRIiMjC520UbNmTezYsUNo4X748CHu3r0LALCzs0NUVBQAICYmBj/88IMw/u+///6Dra0tAODx48f48ccfhfo4OzsLLdcA8ODBAwCKoK1t27awsbHBtWvXMHPmTKFu169fR61atYT3XLx4Ea6urjA2NoaTkxMePnyIp0+fAlC0eCYmJqJXr15CC2e+mJgYSKVSGBoalurZMsYY07ypU6cqxRkJCQlavZ5OLATt6OgIAEJ3Vb6HDx8K5xwdHfHo0SOl8y9evMDTp0+V8qgr4/VrFJQn/3xBLCwslF7GxsYluUVWxq5evQpbW1tYWVkBAMaMGYO0tDR4eHigWbNmaNu2LSZNmgQARQaA06ZNw7179+Dm5oYmTZpgwIABwuSK7777DgsXLkTz5s0xffp0WFhYCGUNHjwYa9aswcCBAzFv3jyYmpoKAeC3336LqKgoSKVSNGnSRJioFBISgtWrV6NZs2aYNWsWnJ2dhVbMjh07IjY2Fl5eXpDL5UotjA4ODpg1axbatm0Ld3d37NixA/v374eenp7K/RV1v4wxxsqesbGxSqyhTSIq7cC9t7moSIQdO3YI45+ICM7Ozpg0aRI+//xzAIoxd/b29ggPD8cHH3yAy5cvQyqVKo3rOnDgAHr06IG7d+/C2dkZy5Ytw9dff42HDx8KrRtfffUV5HK5MEu5f//+SE9Px++//y7Up02bNvD09MTy5ctV6pqSkgJLS0skJydr/cNgjDHGGAO0H3+UWQtgamqqMH4KUEz8iI6Oxp07dyASiTB+/Hh899132LVrF2JjYzFkyBA4OzsLQaK7uzt69OiBUaNG4cyZMzhx4gTGjh2LDz74AM7OzgCADz/8EEZGRhgxYgQuXryI//3vf1i8eDEmTpwo1GPcuHHYv38/Fi5ciCtXrmD69Ok4d+4cxo4dW1aPgjHGGGOsfGllbrEahw8fJgAqr+DgYCJSLAXz7bffkoODAxkbG1PXrl0pLi5OqYwnT57QgAEDyNzcnCwsLGjYsGH0/PlzpTwxMTHUrl07MjY2pho1atDcuXNV6vLrr7+SRCIhIyMjatSoEe3Zs6fAevMyMIwxxhgra9qOP4rdBdy+fXv4+/vjnXfegUQi0VpAqmu4C5gxxhhjZU1nuoBHjRqFU6dOoXnz5nB3d8cXX3yBEydOlHrtP8ZYxbd161bUqFEDDRs2RExMTHlXhzHGWDEVOwAcMmQIIiIi8PjxYyxcuBBJSUl477334OjoiOHDh2Pnzp3IyMjQZl0ZYzokISEBQ4cOhaurK4yMjNRuocd0U0BAAKytrdGvX7/yrgpjrJyUeBKIsbExevXqhRUrVuD+/fvYtWsXnJyc8O2338LW1hZ9+vTBiRMntFFXxjRKJBKV+MVe+eeff5CVlYVvv/0W8+bNQ6NGjZCbm1ve1WLFMG7cOGGPacZY1fTWs4BbtWqFWbNmITY2FrGxsejatauwqC1jumrJkiUYPHgwiKhEL/bKvXv3ACi2sOvZsydWr14t7FLCdFunTp1QrVq18q4GY6wcaXQZmPr162PChAncrcB0GhFh+/btWL9+fXlXpUJLTU0FAJibm2ul/E6dOhW6/WRVQET46KOPYGNjA5FIhOjoaH4ujDGN0ImdQBgrS0eOHEHHjh213qV77Ngx9O3bF87OzhCJRNi5c6fS+WXLlsHT01NY8d3Hxwf79u0rtMyiuqinT59eZL369u2LHj16qD13/PhxiEQinD9/vshyCgsAhw4dKqzhqSllFfgU93O5d+8eBg0aBFtbW5iYmKBx48Y4d+5cgeUW9fOgzv79+xEeHo7du3fjwYMH8PDwKNY9NGnSBB4eHiqv+/fvF+v9jLHKz6C8K1BReHt7Q19fHyEhIQgJCSnv6rC3sGnTJnz55Zdav05aWhq8vLwwfPhwBAYGqpyvWbMm5s6diwYNGoCIsH79evj7++Pff/9Fo0aN1Jb5+vCK//3vf5g2bRri4uKEtOK0xo0YMQJBQUG4e/cuatasqXRu3bp1aNGiBTw9PYssR9stgOWlOJ/Ls2fP0LZtW3Tu3Bn79u2DnZ0drl69Cmtr6wLLLernQZ3r16/DyckJbdq0KdE95C+4zxireMLCwhAWFqb9MdVaWV2wEuGFoCuX9PR08vPzU0nPzc2lWbNmkaurKxkbG5O9vb2wSLkmAKAdO3YUmc/a2ppWr15drDLXrVtHlpaWKum5ubk0e/ZsqlOnDonFYvL09KRt27YJ53NycsjBwYFmzpyp9L7nz5+Tubk5LVu2rFjXDw4OJkNDwwLP+fv7F6scIqLU1FQaPHgwmZmZkaOjI33//ffUsWNHGjdunFAe3lhE/ubNm8Uu/229+bl88cUX1K5du1KXV5yfhzfvuXbt2kRE1LFjRwoJCaGQkBCysLAgW1tb+uabbygvL69EdTh8+DAFBQWV8g4YY9qm7fjjrbuAc3NzER0djWfPnr1tUYxp1OnTp/Hee+9h3rx5Qtpvv/0Gf39/lbxz5szB1q1bsXLlSsTFxWHHjh3o0KGDUp7Zs2fD3Ny80NedO3dKVdfc3Fxs3boVaWlp8PHxKVUZr9/Lhg0bsHz5cly8eBETJkzAoEGDcPToUQCAgYEBhgwZgvDwcKWJLdu2bUNubi4GDBhQrOs8f/5cY61/kydPxtGjR/Hbb7/hwIEDOHLkCKKiooTzixcvho+PD0aNGoUHDx7gwYMHcHFxUVuWJj+ngj6XXbt2oUWLFnjvvfdgb2+Ppk2bYtWqVW/3EN6wePFizJgxAzVr1sSDBw9w9uxZ4dz69ethYGCAM2fOYPHixVi0aBFWr15d7LJ9fX3x3nvvYe/evahZsyZOnTql0bozxiqAkkaM48aNE/4SfvHiBbVt25ZEIhGZmZnR4cOHNRyflj9uAay4Hj16ROvWrSOJRCKkBQQE0NOnT1Xytm/fnr766qtCy3vy5AldvXq10FdOTo7a96KAFp/z58+TmZkZ6evrk6WlZaHbEr5JXQtgZmYmmZqa0smTJ5XSR4wYQQMGDBCOL1++TACU/s22b9+eBg0aVOzrd+/enVxcXNSeK0kL4PPnz8nIyIh+/fVXIe3JkydkYmIitAASkVKLYGHe5nPKV9TnYmxsTMbGxjR16lSKioqiFStWkFgspvDw8GLdc0E/D2/64YcfhJa/fB07diR3d3elFr8vvviC3N3di3VtxljFoO34o8RjALdv345BgwYBAH7//XfcvHkTV65cwcaNG/H111/zGoBMZ9jZ2aFfv34ICQnB2bNnUbt2bRgZGakdp/XOO+/giy++wLlz5/Dee+8hKChIJZ+NjQ1sbGw0Wkc3NzdER0cjOTkZ27dvR3BwMI4ePQqpVFqq8q5du4b09HR069ZNKT07OxtNmzYVjhs2bIg2bdpg7dq16NSpE65du4bjx49jxowZxb5WamqqRloAr1+/juzsbLRq1UpIs7GxgZubW6nK08TnVNTnkpeXhxYtWmD27NkAgKZNm+LChQtYvnw5goOD3+raxdG6dWulSUw+Pj5YuHAhcnNzeSkexlixlLgL+PHjx3B0dAQA7N27F++99x4kEgmGDx+O2NhYjVeQsbdhbm4Of39/bN68Gb/88kuB3ZuTJk3C5cuX0bVrV/zwww9wdXXFzZs3lfJoowvYyMgIrq6uaN68OebMmQMvLy8sXry41PebPzFjz549iI6OFl6XLl3C9u3blfKOGDECEREReP78OdatW4f69eujY8eOJbqWLk4A0cTnVNTn4uTkpBKku7u7l3oIAGOMlbUStwA6ODjg0qVLcHJywv79+7Fs2TIAQHp6Ov/lyXTSwIEDMWLECHh6emLPnj0F5pNIJJgyZQo+++wzWFhY4NKlS6hbt65wfsyYMXj//fcLvZazs/Nb1TUvLw9ZWVmlfr9UKoWxsTHu3LlTZDD3/vvvY9y4cdiyZQs2bNiAjz/+uERL46SmpsLW1rbUdc1Xv359GBoa4vTp06hVqxYAxSzb+Ph4pXswMjIq1qw4bXxOb34ubdu2VZp9DQDx8fGoXbt2icotrdOnTysd//PPP2jQoAH/DmaMFVuJA8Bhw4bh/fffh5OTE0QiEXx9fQEofiE1bNhQ4xVk7G35+fkhNzcXbm5uMDQ0VDk/f/58ODo6wtvbG3p6elixYgVsbW1Vlt4oaddiamoqrl27JhzfvHkT0dHRsLGxQa1atTB16lT07NkTtWrVwvPnz7FlyxYcOXIEf/zxR6nvtVq1apg0aRImTJiAvLw8tGvXDsnJyThx4gQsLCyUuifNzc3Rv39/TJ06FSkpKRg6dGiJrpWamlrobhLJyckqy5HY2tqqTN4wNzfHiBEjMHnyZNja2sLe3h5ff/019PSUOyjq1KmD06dP49atWzA3N4eNjY1KHuDtu4CL87lMmDABbdq0wezZs/H+++/jzJkzWLlyJVauXAkA+Omnn7Bjxw4cOnRIeE9RPw8lcefOHUycOBGjR49GVFQUli5dioULF5b6nhljVVBpBg5u27aNFi1aRAkJCUJaeHg47dy5U2ODE3UFTwKpHD799FM6c+aM2nOhoaEkkUhILBZT9erVyd/fny5duvTW1zx8+LDK0iUAhOVlhg8fTrVr1yYjIyOys7Ojrl270oEDB4pdfkHLwOTl5dGPP/5Ibm5uZGhoSHZ2duTn50dHjx5VyXvy5EkCQL169VJbfmG/IkxMTOjDDz9Ue07dsi0AaMSIEWrzP3/+nAYNGkSmpqbk4OBA8+fPV5n0ERcXR61btyYTExOtLgNT3M/l999/Jw8PDzI2NqaGDRvSypUrhXMymUxl8kZRPw/qFDQJ5JNPPqExY8aQhYUFWVtb01dffVXiZWAYY7pN2/GHiKhkG5xu2LAB/fv3h7GxsVJ6dnY2tm7diiFDhpQ6GNVFKSkpsLS0RHJyMiwsLMq7OoyVGZlMhqNHj+LIkSMq5/Ly8qCvr4+PPvoIK1asKPvKMcZYJaft+KPEk0CGDRuG5ORklfTnz59j2LBhGqkUY6z87du3D/Pnz1dK27BhAw4cOIAnT54AUMy0ZowxVvGUeAwgEakdKH737l1YWlpqpFKMsfJ35swZlTS5XI74+Hhhn19vb+8yrhVjjDFNKHYA2LRpU2HD+a5du8LA4NVbc3NzcfPmzQI3mGeMVQ7Tp09H165dMWfOHAwcOBDvvPNOeVeJMcZYKRQ7AMz/iz86Ohp+fn5K638ZGRmhTp06CAoK0ngFGWO6o0mTJrhz5w6ePn1a4FZsjDHGdF+JJ4GsX78e/fv3h1gs1laddEr+IEyJRAJ9fX2EhIQgJCSkvKvFGGOMsUooLCwMYWFhyM3NRXx8vNYmgZQ4AMyXnZ2NR48eIS8vTym9pOtZ6TqeBcwYY4yxsqbt+KPEk0CuXr2K4cOH4+TJk0rp+ZNDirNSP2OMMcYYKz8lDgCHDh0KAwMD7N69W9gNhDFd1rFjR9SvXx9r165VSv/5558xZcoUpKSkqN1RgjHGGKusShwARkdHIzIykrd9YxUCEeHff/9VuzfsuXPn0KRJEw7+GGOMVTkl/p9PKpXi8ePH2qgLYxp39epVPH/+HM2bN1c5d+7cObXpBcnNzcXly5cLPC+RSODj44OMjAwhjYjQunVrTJ06tcjzjDHGWFkpcQA4b948TJkyBUeOHMGTJ0+QkpKi9GJMl0RGRkJfXx9eXl5K6RkZGbh06RKaNWtW7LL27NmD1q1bq10gGQD+97//ISoqCidOnBDSNm/ejNu3b+Orr74q8jxjjDFWVkocAPr6+uKff/5B165dYW9vD2tra1hbW8PKygrW1tbaqCNjpRYVFYXc3FyYmpoKC5mLRCKYmpoiNze3RAHgO++8gylTpqB79+44d+6cyvmmTZuiSZMmuHLlCgAgPT0dU6dOxXfffYdq1aoVeZ4xxhgrKyUeA3j48GFt1IMxrYiKikJAQACmTZumlL5161YsWbIEUqlUKf3IkSPo3LlzkeX2798f169fV0mXSCSIi4sDAMyfPx/Vq1dX2iO7qPOMMcZYWShxANixY0dt1IMxrYiKikJoaCiaNGmilP7zzz/D09MT+vr6SuktW7YsdJzfr7/+ipkzZ2LBggVqz7u5ueHYsWO4e/cuFixYgD179ihNMinqPGOMMVYWSvU/z/HjxzFo0CC0adMG9+7dAwBs3LgRf//9t0Yrx9jbuHHjBpKSktR280ZFRamdAGJqaoqGDRuqfSUkJGDOnDn45ZdfEBgYqPaa+S18X375Jbp3745OnTqV6DxjjDFWFkocAEZERMDPzw8mJiaIiopCVlYWACA5ORmzZ8/WeAUZK63IyEjo6emptP7l5OTgwoULJRr/BwCtW7fGzp070a9fvwLzSCQSJCQkYPv27WpbCYs6zxhjjJWFEgeA3333HZYvX45Vq1bB0NBQSG/bti2ioqI0WjnG3kZUVBQaNGgAc3NzpfRLly4hKyurxAFgtWrV4OfnV2geiUQCABg7dixcXV1LfJ4xxhgrCyUOAOPi4tChQweVdEtLSyQlJWmiTjrJ29sbUqkUYWFh5V0VVkxz5swRZty+zsvLC0SEpk2bavyamZmZICIMGTKkVOcZY4xVbWFhYZBKpfD29tbqdUo8CcTR0RHXrl1DnTp1lNL//vtv1KtXT1P10jlnz57VymbMrHKJiYmBkZER3N3dS3WeMcZY1RYSEoKQkBCkpKTA0tJSa9cpcQvgqFGjMG7cOJw+fRoikQj379/H5s2bMWnSJHz88cfaqCNjFUZMTAykUqnS8IiSnGeMMcbKgoiIqCRvICLMnj0bc+bMQXp6OgDA2NgYkyZNwsyZM7VSyfKUH4EnJydzCyBjjDHGyoS2448SB4D5srOzce3aNaSmpkIqlaoMtK8sOABkjDHGWFnTdvxR4i7gTZs2IT09HUZGRpBKpWjZsmWlDf4YY4wxxiqjEgeAEyZMgL29PT788EPs3bsXubm52qgXY4wxxhjTkhIHgA8ePMDWrVshEonw/vvvw8nJCSEhITh58qQ26scYY4wxxjSs1GMAASA9PR07duzAli1b8Oeff6JmzZq4fv26JutX7ngMIGOMMcbKmrbjjxKvA/g6U1NT+Pn54dmzZ7h9+zYuX76sqXoxxhhjjDEtKXEXMKBo+du8eTN69eqFGjVq4Mcff0RAQAAuXrxY6ork5ubi22+/Rd26dWFiYoL69etj5syZeL2Bkogwbdo0ODk5wcTEBL6+vrh69apSOU+fPsXAgQNhYWEBKysrjBgxAqmpqUp5zp8/j/bt20MsFsPFxQXz588vdb0ZY4wxxiqaEgeAH3zwAezt7TFhwgTUq1cPR44cwbVr1zBz5kw0bNiw1BWZN28eli1bhp9++gmXL1/GvHnzMH/+fCxdulTIM3/+fCxZsgTLly/H6dOnYWZmBj8/P2RmZgp5Bg4ciIsXL+LgwYPYvXs3jh07ho8++kg4n5KSgu7du6N27dqIjIzEggULMH36dKxcubLUdWeMMcYYq0hKPAZw4MCBGDhwIPz8/KCvr6+xivTp0wcODg5Ys2aNkBYUFAQTExNs2rQJRARnZ2d8/vnnmDRpEgAgOTkZDg4OCA8PxwcffIDLly9DKpXi7NmzaNGiBQBg//796NWrF+7evQtnZ2csW7YMX3/9NRITE2FkZAQA+PLLL7Fz5061+8byGEDGGGOMlTWdWwcwv+tXk8EfALRp0waHDh1CfHw8AMWWWX///Td69uwJALh58yYSExPh6+srvMfS0hKtWrXCqVOnAACnTp2ClZWVEPwBgK+vL/T09HD69GkhT4cOHYTgDwD8/PwQFxeHZ8+eFVi/lJQUpVdWVpbmbp4xxhhjVVpWVpZKrKFNxQ4Ae/XqheTkZOF47ty5SEpKEo6fPHkCqVRa6op8+eWX+OCDD9CwYUMYGhqiadOmGD9+PAYOHAgASExMBAA4ODgovc/BwUE4l5iYCHt7e6XzBgYGsLGxUcqjrozXr6GOi4sLLC0thdecOXNKfa+MMcYYY6+bM2eOUpzh4uKi1esVexbwH3/8odTqNXv2bLz//vuwsrICALx48QJxcXGlrsivv/6KzZs3Y8uWLWjUqBGio6Mxfvx4ODs7Izg4uNTlakpCQoJSE6yxsXE51oYxxhhjlcnUqVMxceJE4TglJUWrQWCxA8A3hwq+xfKBak2ePFloBQSAxo0b4/bt25gzZw6Cg4Ph6OgIAHj48CGcnJyE9z18+BBNmjQBADg6OuLRo0dK5b548QJPnz4V3u/o6IiHDx8q5ck/zs+jjoWFBY8BZIwxxphWGBsbl2njUqmWgdGG9PR06OkpV0dfXx95eXkAgLp168LR0RGHDh0SzqekpOD06dPw8fEBAPj4+CApKQmRkZFCnr/++gt5eXlo1aqVkOfYsWPIyckR8hw8eBBubm6wtrbW2v0xxhhjjOmKYgeAIpEIIpFIJU1T+vbti1mzZmHPnj24desWduzYgUWLFiEgIEC41vjx4/Hdd99h165diI2NxZAhQ+Ds7Ix3330XAODu7o4ePXpg1KhROHPmDE6cOIGxY8figw8+gLOzMwDgww8/hJGREUaMGIGLFy/if//7HxYvXqzU7FqovDxAw62fjDHGGGNlqURdwEOHDhWaJzMzMzFmzBiYmZkBwFvPil26dCm+/fZbfPLJJ3j06BGcnZ0xevRoTJs2TcgzZcoUpKWl4aOPPkJSUhLatWuH/fv3QywWC3k2b96MsWPHomvXrtDT00NQUBCWLFkinLe0tMSBAwcQEhKC5s2bo3r16pg2bZrSWoEFkcvlCJ0wAfF37kBSrRpk06cjsLiBI2OMMcaYjij2OoDDhg0rVoHr1q17qwrpmvx1eDZu3IjBgwdDBIAA4WuEpSUCP/4Y+PZbwNS0fCvLGGOMsUpB2+sAlngh6Kom/wNo1KgRLl26pDT5RQTAE0A0AIhEQPPmgEwG9OlTPpVljDHGWKWgcwtBV1XXrl1TnQkNIE4kAvT0FOMCz50D+vbFdldX2L4cM9m4cWPI5XKV8rZs2SKMqywoD2OMMcaYNnAAWEyurq5qJ8G4eXoCz58DX30FODpCDuC969eRv6fIxQsXEBQUBPn27cL75HK5sMA1AFy8eFGRh4NAxhhjjJUBDgCL6csvv1Tu/hWJQESQyWSKsX+zZgEPHiDU3V0YHwi8Gi844/33gXffBeLiEBoaqhRMEhFEIhFmzJhRhnfEGGOMsaqKxwAW4fU++D///BMzZsxAXFwc3NzcIJPJhGVq8pmYmCAzM1OlHDGAjPw8AFRzAGKxGBkZGWrOMMYYY6wq0fYYwGIvA8OAwMBABAYGFppHIpEgNjZWpbXQzdISyM4G0tMhARCLV62EQh43N63UmzHGGGPsddwFXEze3t6QSqUICwsrNJ9MJhO6dIHXuorXrgXS0oAtWyCrVUvoGsbLr0QEmZubIkhkjDHGWJUUFhYGqVQKb29vrV6Hu4CLUJomWLlcXmRXsfyXXzDjiy8Qd/cu3IggAxAAAPr6QLt2wLx52HL9ujBZxMPDA6GhoUW2QDLGGGOs4uN1AMuZtj8AAMCJE8BnnwH//itsMycHEPRalvyWxIiICA4CGWOMsUqO1wGsCtq2BSIjgadPgbFjASsrhOJVFzHAM4UZY4wxpjkcAOoSKytg6VLg2TPEGxrizaZZIkLc+fPAtGk8VpAxxhhjpcYBoI6SuLurLjwNwI0ImDlTsfZg166Qf/89vLy8YGJiAi8vL7WLScvl8iLzMMYYY6zq4DGARSiTMYBqyOVyBAUFCWP/8r/K69VDwM2bAJEwTjB/4Wl14wTzy8nHYwkZY4wx3cdjAKuowMBAREREwNPTE2KxGJ6enpDL5Qi4fl0xVnDMGITq6SnvOpI/TnDaNKEc3nWEMcYYY2/iFsAilFcLYHEUuutIUBDQsSNMJk9GZlaWah7edYQxxhjTWdwCyAokkUjUjxMEgIgI4LPPIMnOhuiN9/GuI4wxxljVxgFgBZa/60g+kUgEAiBbsABo3lyRh0h515H8nUlksjKvL2OMMcZ0AweAFVj+OEEvLy/lcYKTJgHnzgG3byPwvfcQYWAATyi6hj2JIBeLEfD0KcBdwIwxxliVxGMAi6DLYwCLLS8PWL4cmD8fuH37Vbq1NTBsGPDuu5D/9x9CQ0MRHx8PiUQCmUzGs4QZY4yxcsJbwZWzShEAvu7iRWD9emDbNuDWLQC87RxjjDGma3gSiI7w9vaGVCpFWFhYeVfl7TRqpGgJvHYN2L0baNyYt51jjDHGdERYWBikUim8vb21eh1uASxCpWsBVMPE2BiZaraWEwPIiIgAAgIA0ZtzidWTy+XclcwYY4y9JW4BZFonadiw4OVkgoIADw8gLAxITi60nPxdR86fP4/MzEzExsYiKCiIt55jjDHGdAwHgKzg5WQ8PRV7Dl+6BIwdC9jYKLqQIyLUlsO7jjDGGGMVAweArODlZGJigPv3gSVLgFq1FLOJL10C+vUDLC2BceOA1FShnPj4eLw5ooCIEBcXV9a3xBhjjLFC8BjAIlSFMYDFkpcH+fDhCN28GfEvXkACQAYgUE8PaNsW+PFHeA0bhtjYWJXWRE9PT0RHR5dXzRljjLEKh8cAMp0g37kTQevX4/yLF8gEEAvF0jHyvDzg+HHA2xsyU1Oh2xfgXUcYY4wxXcUBICsWlfF9UAR4M6pXB+ztgbw8BP7zDyIAeOrrQ6ynB093d0VXckBAudWbMcYYY6o4AGTFUuD4vtRU4OFD4Px5YMwYBJqaIvrFC2Tk5SH60iUETJsG7NpVTrVmjDHGmDocALJikUgkqkvFiERwc3NTHDRuDCxbBly/rlg6xshIkX7hAuDvr9h2bvJkIDOzjGvOGGOMsTdxAMiKRe1SMerG9zk6Atu3AxkZih1HnJwU6UlJwPffK7qLZ84EEhPLrvKMMcYYU8IBICuWApeKKWh8n56eosXv/n3g5EmgTRvFbiLPn0M+bRq8nJxgoq8Pr5o1Id++vWxvhjHGGKviOABkxRYYGIjo6GhkZGQgOjq6+JM7fHyAEycUwd+4cQiCYhZxZl4eYu/dQ9B770HeqRPw+LHwFrlcDi8vL5iYmMDLy4t3E2GMMcY0iNcBLAKvA6hZXl5eqmsFAvAEEC0SAc2aQd6nD4JCQ1+df9ndHBERwfsKM8YYqxK0HX9wAFiE/A9AIpFAX18fISEhCAkJKe9qVVgmJibIVDMRRAwg4+X3XlC0EL7+g8kLSjPGGKsKwsLCEBYWhtzcXMTHx3MAWF64BVCz1LYAikTwbNwY0QEBQFgYTB4/hrq5wmKxGBkZGWrOMMYYY5UL7wTCKpX82cQqu4VMnw5Mnw789x8k9epB9Mb7RADciIBffinjGjPGGGOVDweArEzlzyb29PQscDaxbMECYacRQBH8EQBZVhbw4YeAhQXw6adAamq53ANjjDFW0XEXcBG4C7h8yOVyzJgxA3FxcXCTSCBr3RoB+/cDd+68yqSnB7RtC2zeDLi4lF9lGWOMMQ3jSSDljANAHRMdDUycCBw9CnleHkIBxAOQWFpC9vnnCPzmG8V6g4wxxlgFxmMAGXtdkybAX39BvmnTq/UEAcQmJyNo2jTIjY0V3cQvdxrh9QQZY4wxVRwAsgopdO5cxQSSl8cExVjBGTk5iokiTk6Q16mDoKAgnD9/HpmZmYiNjUVQUJBKEMhBImOMsaqGu4CLwF3AuqnA9QT19JChrw/k5BRrPUG5XI6goCCl87zoNGOMsfJWpbqA7927h0GDBsHW1hYmJiZo3Lgxzp07J5wnIkybNg1OTk4wMTGBr68vrl69qlTG06dPMXDgQFhYWMDKygojRoxA6huzRc+fP4/27dtDLBbDxcUF8+fPL5P7Y5ojkUiEWcL5RCIR3Bo3BtLTgdmzEQ/l4A9Q/AzFXbkiHIeGhiqVk79EzYwZM7RYe8YYY6x86UwA+OzZM7Rt2xaGhobYt28fLl26hIULF8La2lrIM3/+fCxZsgTLly/H6dOnYWZmBj8/P6WWoIEDB+LixYs4ePAgdu/ejWPHjuGjjz4SzqekpKB79+6oXbs2IiMjsWDBAkyfPh0rV64s0/tlb6fA9QRlMsDAAJg6FRJPT9UgEYBbdrZiLOG6dYiPj8ebjeBEhLi4uDK6E8YYY6wckI744osvqF27dgWez8vLI0dHR1qwYIGQlpSURMbGxvTLL78QEdGlS5cIAJ09e1bIs2/fPhKJRHTv3j0iIvr555/J2tqasrKylK7t5uam9rrJyckEgJKTk9/q/pjmRUREkJeXF4nFYvLy8iK5XK5yHgCJRCKlr3KA6OXLUyQikaKhUHiJRCLy8vIqn5tijDHGSPvxh860AO7atQstWrTAe++9B3t7ezRt2hSrVq0Szt+8eROJiYnw9fUV0iwtLdGqVSucOnUKAHDq1ClYWVmhRYsWQh5fX1/o6enh9OnTQp4OHTrAyMhIyOPn54e4uDg8e/aswPqlpKQovbKysjR276x0AgMDER0djYyMDERHRystJp1/Xu2i09u2AfXqAQBkRMIEEuCNlkTGGGOsjGRlZanEGtqkMwHgjRs3sGzZMjRo0AB//PEHPv74Y3z22WdYv349ACDx5bIeDg4OSu9zcHAQziUmJsLe3l7pvIGBAWxsbJTyqCvj9Wuo4+LiAktLS+E1Z86ct7hbVlbUBon9+gHXrwOXLiHQzw8RenrwBCAG4EkEuZsbAkxMgLy88q4+Y4yxKmLOnDlKcYaLljc4MNBq6SWQl5eHFi1aYPbs2QCApk2b4sKFC1i+fDmCg4PLuXZAQkKC0iwcY2PjcqwN0wh3d2D/fgRmZyNw+nRgxQrg6VMgLg7o2RNwdQV69QLGjwfq1i3v2jLGGKvEpk6diokTJwrHKSkpWg0CdaYF0MnJCVKpVCnN3d0dd15u/eXo6AgAePjwoVKehw8fCuccHR3x6NEjpfMvXrzA06dPlfKoK+P1a6hjYWGh9OIAsBIxMgJmzwaePAHOnQPGjQMsLYFr14AlS4B69SB3doZX7dq8ViBjjDGtMDY2Vok1tElnAsC2bduqzLyMj49H7dq1AQB169aFo6MjDh06JJxPSUnB6dOn4ePjAwDw8fFBUlISIiMjhTx//fUX8vLy0KpVKyHPsWPHkJOTI+Q5ePAg3NzclGYcsyqqeXPgxx+Be/eAWbMAsRhyAEEPHiD2zh3FgtLnz6tdUJoxxhirMLQytaQUzpw5QwYGBjRr1iy6evUqbd68mUxNTWnTpk1Cnrlz55KVlRX99ttvdP78efL396e6detSRkaGkKdHjx7UtGlTOn36NP3999/UoEEDGjBggHA+KSmJHBwcaPDgwXThwgXaunUrmZqa0ooVK9TWi2cBV3G5ueRZo4bqTGGAvCwtiW7dKu8aMsYYq4S0HX/oTABIRPT777+Th4cHGRsbU8OGDWnlypVK5/Py8ujbb78lBwcHMjY2pq5du1JcXJxSnidPntCAAQPI3NycLCwsaNiwYfT8+XOlPDExMdSuXTsyNjamGjVq0Ny5cwusEweATCwWKwV/+S8xQKSnR/TOO0QREUQ5OeVdVcYYY5WEtuMP3gquCLwVHPPy8kJsbKzSgtEiAJ7m5oh+fZcZQ0OgRw/ghx+A+vXLvqKMMcYqjSq1FRxjukjtriMAZBs2AJcvA2PHAiIRkJMD/P67YvawRAJs2FC+FWeMMcYKwAEgY0UocEHpgACgYUNg6VIgMRHyjh3hJRLBBIDX1auQBwcD5uaKADEjQyhPLpfDy8uLZxQzxhgrN9wFXATuAmbFIZfLERQUJByLoBgoGAEgEACsrYFhwyCvWxdBn376Kt/LnUciIiIQGBhYxrVmjDGmq7Qdf3AAWAQOAFlxqB0nKBIpxglaWgJ37yryAYiFIjhUyufpiejo6LKsMmOMMR3GYwAZqwDi4+Px5t9SRIS4nBzg1i1gzx6gd2/EQzn4E/K9sQYmY4wxpk0cADKmARKJRJgkkk8kEsHNzQ3Q11dsKbd7NyT160P0xntFANxevAC+/17Yf5jHCTLGGNMmDgAZ0wC1M4WJIJPJlPPNnw96eR54NVZQ9uIFMHky4O4O+bBhCAoKwvnz5xU7j8TG8s4jjDHGNIoDwGLy9vaGVCpFWFhYeVeF6aBCZwoXls/LC/LJkxHQqBFgbAzExyM0PFyplTA/sJwxY0bZ3hRjjLEyFxYWBqlUCm9vb61ehyeBFIEngbAy8/w5sHkzTD75BJlq/lmKjY2RkZkpHMvlcoSGhiI+Ph4SiQQymYxnEjPGWCXBs4DLGQeArKwVuPMIgGgPD2D6dMhFIuVlZ3g5GcYYq1Q4ACxnHACyspa/pmB+UJc/TlAOIL9D2cvEBLGZmarLzvByMowxVinwMjCMVTFqxwlu24aA+fOBmjUBkQjxGRnql53h5WQYY4wVA7cAFoFbAJnOuX8fXs2bIzYxUXlBaQCe1asj+vp1gH9WGWOsQuMWQMaYMmdnyMLC1C8n8/gxYGUFtGwJHDpUjpVkjDGmyzgAZKwCUukmbtAAcqkUAXp6ABFw9izg6ws0bQps2wZkZ5d3lRljjOkQ7gIuAncBswolMxOYMQNYvRr4779X6fb2ioBw+HCga9fyqx9jjLFi4S5gxljxicXA7NnAo0fA+fOQBwXBy8AAJo8ewWvLFsh9fQEbG2DcOMW6g+Bt5xhjrCriFsAicAsgq6gKWk4mAkAgAIhEkNerh6Dr11/l4fUEGWNMJ3ALIGOsVEJDQ4WADoAwaWSGmRnwcqxg6PXrisAwP08B285xKyFjjFUuHAAyVknFx8erXyswNxdISwO++QbxAN7sAiAixF28COTlAXjVknj+/HlkZmYiNjYWQUFBHAQyxlgFxgFgMXl7e0MqlSIsLKy8q8JYsUgkEmGZmHwikQhubm6KsYIzZ0Li6amaB4DbixeAqSnwzjsI/fJLpTwFtRIyxhh7e2FhYZBKpfD29tbqdXgMYBF4DCCrqFTGAL78KpfLERAQUHgevNp2zgRAppryxWIxMjIyyuhuGGOsauExgIyxUlFZK9DTUyn4KzTPjRvABx8ApqaQQNEq+DoRADcnJyA1VUjjcYKMMVZxcAtgEbgFkFVpeXmQT5mCoIULhVnE+V/lAAL09YGWLSF3dkZQRITwtoJmE8vlcoSGhiI+Ph4SiQQymUxltnFx8jDGWGWn7fiDA8AicADImCIomzFtGuKuXoWboyNkDRog4OJFIDERAOAFIBZQ3ZvY1BTR3t6AoSHkT58iKCpKJZCMCAxEYJMmgKkp5HFxCFq16lUZvCwNY6yK4gCwnHEAyFgB1qwBxo8HUlMLHicIIH+UYIFBIoDowvKIRPD09ER0dDQYY6yq4DGAjDHdNGKEogUwPBwSU1PVcYIiEdzq1AF+/RXYtAnxBgaqS84AiNPXB0aNAgYNQryenvplaWJjgYsXtXYrjDFW1XAAyBgrPTMzIDgYso0bhW5dvPxKRJAtWgS89x7g4gKJm5v6ZWk8PICVK4GNGyHx8FC/LE1eHtC4MTB4MHDtWlncGWOMVWocADLG3powm9jLC2JjY3g2bPhqxnF6OtC3L2QvF6bOD/Dyx/fJZDKhHJlMppoHgKxVK4AI2LQJaNgQ8u7d4SWV8oxjxhgrJR4DWAQeA8jYW7pyBejXD7h4EXIAMwDEiURws7eHbPp0BIwZo5RdLpdjxowZiIuLg5ubG2QymSKQjIwEvvkG8v37EYRXk0h4oghjrDLiSSDljANAxjSACDh3Dli7FtiyBUhJeXVu1Spg5MhiF+Xl6orY69d5oghjrFLjSSCMsYpPJAK8vYFly4AHD4DNmwFfX0BPD+jc+VW+qCjg2DFFwFiA+Hv31E8UuXRJO3VnjLFKiANAxljZMjUFPvwQOHgQuH8fqF//1bmZM4GOHYEGDRTf376t8na1exwDcMvJAaZPB/LytFt/xhirBDgAZIyVHweHV98TAY6OgLk5cP06MG0aUKeOooUwPBx4/hxAIRNFACA0FOjdG3j6tIxvhDHGKhYOAIvJ29sbUqkUYWFh5V0VxionkUjRRZyYCKxfD3Tpokg7cgQYNgzo1QtAIfsXr18PiMXA/v1A8+aQf/89703MGKtwwsLCIJVK4e3trdXr8CSQIvAkEMbK0Z07wMaNioBw7Fjgs88U6cnJwPz5inUBGzZ8lT8mBggMhPzGDQS9VgzPFGaMVTQ8C7iccQDImA4gAl68AAwNFcdr1ryaOeztrQgEP/gAsLMDnj2DV+3aiH3+nGcKM8YqLJ4FzBhjItGr4A9QTBzp0wcwMADOnlW0DDo7A337AgcPIj47W/1M4bi4Mq02Y4zpKg4AGWMVT6dOwO+/K2YRL1kCtGihaCHcvRsYMACSevXUzxTW01OsR8gYY1UcB4CMsYrLzg749FNFK+Dly8BXXyn2Jv7uO+WZwlDsGiJLTwdatlR0Hz96VK5VZ4yx8sQBIGOscmjYEJg1C1i79tVMYXd3iAF4ApADCLC2VownXLMGkEiAxYuBnJxyrjhjjJU9DgAZY5VSYGAgoiMjkfHrr4j290eAoSHw7NmrDMnJkI8fD69q1WBiZMRLxTDGqhSeBVwEngXMWCXx5AmwfTuwaRPw99+QAwjCq+7h/K8RGzcicNCg8qwpY4xV3VnAc+fOhUgkwvjx44W0zMxMhISEwNbWFubm5ggKCsLDhw+V3nfnzh307t0bpqamsLe3x+TJk/HixQulPEeOHEGzZs1gbGwMV1dXhIeHl8EdMcbKla0tMHo0cPw4cPMmQh0dhV1EgFdB4IzBgxWLTp8+DQCQy+W8oDRjrNLRyQDw7NmzWLFiBTw9PZXSJ0yYgN9//x3btm3D0aNHcf/+faVFXXNzc9G7d29kZ2fj5MmTWL9+PcLDwzFt2jQhz82bN9G7d2907twZ0dHRGD9+PEaOHIk//vijzO6PMVbO6tRBfFIS3uwAIQBxALBvH9C6NeSWlggKCsL58+eRmZmJ2NhYBAUFqQSBHCQyxioc0jHPnz+nBg0a0MGDB6ljx440btw4IiJKSkoiQ0ND2rZtm5D38uXLBIBOnTpFRER79+4lPT09SkxMFPIsW7aMLCwsKCsri4iIpkyZQo0aNVK6Zv/+/cnPz09tfZKTkwkAJScna/I2GWPlzNPTk0QiEUER9xEAEolE5GVhQWRkRASQJ0Ci184Leby8hHIiIiKE9Ne/RkRElN/NMcYqPG3HHzrXAhgSEoLevXvD19dXKT0yMhI5OTlK6Q0bNkStWrVw6tQpAMCpU6fQuHFjOLy2wbyfnx9SUlJw8eJFIc+bZfv5+QllFCQlJUXplZWV9Vb3yRgrXzKZTHmpmJfbxcnCw4HUVGDLFsS/1kWcj4gQFxsL7NwJECE0NFR4b/55kUiEGTNmlOn9MMYqtqysLJVYQ5t0KgDcunUroqKiMGfOHJVziYmJMDIygpWVlVK6g4MDEhMThTyvB3/55/PPFZYnJSUFGRkZBdbNxcUFlpaWwktdHRljFYewVIynJ8RiMTw9PSGXyxEQEKDYdWTAAEgaN1a/oHReHhAQAJiYIP7CBdWuZN51hDFWQnPmzFGKM1xcXLR6PQOtll4CCQkJGDduHA4ePAixWFze1VGRkJCgNAvH2Ni4HGvDGNOEwMBApXHEb5LJZAgKChJa+IRWQktLIDkZyMqCBEAsoLLvsJubm7arzxirRKZOnYqJEycKxykpKVoNAnWmBTAyMhKPHj1Cs2bNYGBgAAMDAxw9ehRLliyBgYEBHBwckJ2djaSkJKX3PXz4EI6OjgAAR0dHlVnB+cdF5bGwsICJiUmB9bOwsFB6cQDIWOVXYCthUhJw4wYwdChk5ubCDGK8/EpEkDVrBiQklF/lGWMVirGxsUqsoU06EwB27doVsbGxiI6OFl4tWrTAwIEDhe8NDQ1x6NAh4T1xcXG4c+cOfHx8AAA+Pj6IjY3Fo9e2eDp48CAsLCwglUqFPK+XkZ8nvwzGGHtdYGAgoqOjkZGRgejoaEUXMQDUrQusW4fA588R8cMP8LS0VN51ZN06oFYtwNISCAoCIiPL8S4YY0yZTi8E3alTJzRp0gQ//vgjAODjjz/G3r17ER4eDgsLC3z66acAgJMnTwJQLAPTpEkTODs7Y/78+UhMTMTgwYMxcuRIzJ49G4BiGRgPDw+EhIRg+PDh+Ouvv/DZZ59hz5498PPzU6kDLwTNGCuRe/eAiAhg2zbg77+FZDmAUD09xAOQ1K8P2dy5hXY/M8aqtiq7ELQ6P/zwA/r06YOgoCB06NABjo6OSutt6evrY/fu3dDX14ePjw8GDRqEIUOGKM3Gq1u3Lvbs2YODBw/Cy8sLCxcuxOrVq9UGf4wxVmI1agCffaZYcPrMGeDddyEXixEEIDYvD5l5eYi9elWxnuCwYcCdO+VdY8ZYFaTTLYC6gFsAGWNvy8vLC7GxsUqzhUVQdBdHA4C3tyJwHDcO6NSpXOrIGNMt3ALIGGMVXHx8vPpdR0QiQCQCzp5VrCvYuTPkxsbwqlYNJkZGvKsIY0xrOABkjDEtk0gkqusJikRw8/QE7t8HJk8GrK0hBxCUnY3Y1FRk5uQg9vx5RVfx1KkAd9YwxjSIA0DGGNOyAncdkckAR0dg/nzg6VOESiSKZWRevi9/eZkZc+cqZh1PmAD88Qfkv/7Kew8zxt4KjwEsAo8BZIxpglwux4wZMxAXFwc3NzfIZLJXS8q8ZGJigszMTJX3igHk71MkBxAECIFifjAZERHBs4oZq0S0HX9wAFgEDgAZY2VF7WQRkQieHh6InjkTkMvhtWGD6s4jADxtbREdGQnUrl3W1WaMaQFPAmGMsSqiwK7i0FDA3x9Yvx7xxsZ48692AhD35AlQpw7g6wuEhUG+alWR3cRyuZy7khmrorgFsAjcAsgYK0tFdRUXuKSMSITol2nF6SaWy+UICgp6VQZ3JTOmU7gLuJzlfwASiQT6+voICQlBSEhIeVeLMVZF5Qdu+QFb/le5XI4AiQTYtw9eMhli09NVu4nFYkRPmwaEhMCrfXv13c2enoiOji7r22KMvRQWFoawsDDk5uYiPj6eA8Dywi2AjDFdU1QrYXEmk5iIRMhU8+tfLBYjIyNDOJbL5QgNDUV8fDwkEglkMhm3EDJWBrgFsJxxAMgYq2gK7CY2NkZ0Tg6QlwcvQHUyyRstgNxNzFj54UkgjDHGSkTtZBIAsl9+AbKygPBwyIKChHUG8fIrEUEWGws0bgzMmIHQadOUFrDOL/P1/dUZYxUTB4CMMVbJBAYGIiIiAp6enhCLxfD09FSMEQwIAAwMgOBgBG7frsjj7g6xvj489fUhBxCQlwdcuADIZIi/eFF1CzsixMXFKaXxbGLGKh7uAi4CdwEzxqqMx4+BJUuAiAggLg5eubnq1xx0dkb01auAqSl3EzOmJTwGsJxxAMgYq5Ly8iCfPBlBixa9Wk7m5Vc5gAAAcHaGV0YGYpOSeDYxYxrGYwAZY4yVPT09BC5cqOgm9vKC2NgYnjVqQG5vj4D8cYH37yP+2bNidRMzxnQLB4CMMcYKFBgYiOjoaGRkZiL67l0EPHwIZGcD4eFAp06QiEQQvfEeEQA3Q0Pggw+AI0fKvtKMsSJxAMgYY6xkXk4kweHDkG3frjqbGIDs+XPgf/8DOncGjIwAT09g2jQgMbH86s0YE3AAyBhjrNSEGcdeXq9mHH//PQIGDwYsLRWZcnKA2Fhg5kzAyQmwsgI++ww4cwbIzS3X+jNWVfEkkCLwJBDGGHsLd+8CP/8M+aZNCE1IQDwACQAZgEAAsLYG3N2BunWBjz4COnQo1+oypit4EghjjLGKq2ZNyFu0QFBCAmJFImRCsQNJEAC5qSnw7Blw8iSweTPkHTvCSySCiZ4evOzsIP/pJ5XieM1BxjSDA8Bi8vb2hlQqRVhYWHlXhTHGKpTQ0FBhfUDg5ZIyIhFmuLoCp04B3bpBbmKCICiCw0wixD5+jKBPP4VcLAYmTQL++APyLVsQFBSE8+fPIzMzE7GxsQgKClIJAjlIZBVZWFgYpFIpvL29tXod7gIuAncBM8bY2zExMUFmZqZKulgsRkZGBoBC9i8GEP3ymPcvZlUJdwEzxhir0CQSidKewoAiMHNzcxOO4+PjVdcTBBBnYAAMGwbUqoV4KAd/wMs1B2NjgbVrgcxMobXx9fO8fzFjqjgAZIwxplUymUwIxIBXrXIymUzIU2CQ2KiRIri7dQuSBg3UrzmYlweMGAFYWCD+wgVemJqxYuAAkDHGmFYJS8V4er5aKkYuR0BAgJCnyCBRJIJs7lxh/KCQB4DMygoQiYCcHEjy8tQHiXZ2gJpuaMaqKg4AGWOMaZ2wo0hGBqKjo5WCv/zzRQWJBeZ59kwR3P3zD2TBweoXpk5IAExMADs7oEcPYMcOxfqEjFVRPAmkCDwJhDHGKha5XI4Z06cj7soVuInFkOnpISA5WTkPgFCRCPEiESTW1pD16YPA778HqldXKic0NBTx8fGQSCSQyWQ8kYSVGW3HHxwAFoEDQMYYqwSys4GICODXXyE/ehRBz54JrYP5XyMABJqbA+7ukNeqhaCICKErmmcTs7LGAWA54wCQMcYqF00tOcOYNvEyMIwxxpgGFbjkjJERMHIk4OZW8JIzMTGAiwsQEACsXw/5L7/wotOsQuIWwCJwCyBjjFUualsA32jdK04roRyKLe1UupK5m5hpALcAMsYYYxpUnHUJ1eYBIBs+HOjeHbC3RyheBX3AqyBwxoABwNChwPDhwPbtkP/6a5GthLx9HStzxAqVnJxMACg5Obm8q1KmOnbsSOPGjSvvapRLPd7mmo8fPyY7Ozu6efOmRuvE3k7//v3p+++/L+9qMB0SERFBXl5eJBaLycvLi+RyeYnziI2NCYq4T+klBoheviJepone+BoREaF0HQAkEomUvr6eh1U92o4/OAAsQmUNAB89ekRjxowhFxcXMjIyIgcHB+revTv9/fffRET05MkTSklJ0WodgoODlX5p2tjYkJ+fH8XExAh5tBEAvnnd/NfVq1eJSPXeS1KHCRMm0MiRIzVaX205evQo9enTh5ycnAgA7dixQyXP7NmzqUWLFmRubk52dnbk7+9PV65cEc7LZDKV5+jm5las6z948IA+++wzql+/PhkbG5O9vT21adOGfv75Z0pLSytWGX369CE/Pz+1544dO0YAKCYmhmJjY8na2pqSkpKKVS5jxeHp6SkEa/kvkUhEXnXrEn38MZG1NXm+FvThtSDQS0+P6IMPiH74gTzr11dfjpdXed8iK0fajj+4C7iYvL29IZVKERYWVt5V0YigoCD8+++/WL9+PeLj47Fr1y506tQJT548AQDY2NigWrVqb3WNTp06ITw8vNA8PXr0wIMHD/DgwQMcOnQIBgYG6NOnz1tdtzhev27+q27dugBKf+/p6elYs2YNRowYoenqakVaWhq8vLwK/Zk+evQoQkJC8M8//+DgwYPIyclB9+7dkZaWJuRp1KiR0nP8+++/i7z2jRs30LRpUxw4cACzZ8/Gv//+i1OnTmHKlCnYvXs3/vzzz2Ldw4gRI3Dw4EHcvXtX5dy6devQokULeHp6wsPDA/Xr18emTZuKVS5jxVFgV/LChcDPPwNPnyLe2Fh1MgmAuLw8YOtWYMIExF+/rn77uitXVK7JXcWVX1hYGKRSKby9vbV7Ia2ElZVIiSLwvDyi1NTyeeXlFfuenj17RgDoyJEjBeZ5s9UrJSWFPvzwQzI1NSVHR0datGhRkS1jHTt2pHXr1hV4Pjg4mPz9/ZXSjh8/TgDo0aNHauuxe/dusrCwoE2bNpW6Xuqu+2a989+vrrWwoO7dbdu2kZ2dnVLazZs3CQBt376d2rdvT2KxmFq0aEG3b9+mY8eOUatWrcjExIS6dOlCz549K7BO2oYCWgDf9OjRIwJAR48eJSJFC2BpWin8/PyoZs2alJqaqvZ83suf59zcXJo9ezbVqVOHxGIxeXp60rZt24R8OTk55ODgQDNnzlR6//Pnz8nc3JyWLVsmpIWGhlK7du1KXFfGClNUN3GhrYTffUfUty956umpbyUEiKysiAYMINq2jSKWL1cpB9xVXGlxC2BFkp4OmJuXzys9vdjVNDc3h7m5OXbu3ImsrKxivWfixIk4ceIEdu3ahYMHD+L48eOIiooq7ZNSKzU1FZs2bYKrqytsbW1Vzm/ZsgUDBgzA5s2bMXDgwDKp1+LFi+Hj44NRo0YJLVwuLi5q8x4/fhzNmzdXSouJiQEALFu2DLNnz8bJkyfx8OFDDBo0CHPnzsVPP/2Ew4cPIyYmBuvWrStx/WbPni18ngW97ty5U/IbL0Dyy90UbGxshLSrV6/C2dkZ9erVw8CBA4u83pMnT3DgwAGEhITAzMxMbZ78FpU5c+Zgw4YNWL58OS5evIgJEyZg0KBBOHr0KADAwMAAQ4YMQXh4uFILyrZt25Cbm4sBAwYIaS1btsSZM2eK/TPPWHEUtcVdoa2EX38N7NoF2f/+p377OgBISgJ++QV47z2EjhmjtM9xfrkzZsxQuia3ErJi0UpYWYmUKAJPTRUG/pb5q4CWlIJs376drK2tSSwWU5s2bWjq1KkFjr1LSUkhQ0NDpZaXpKQkMjU1fesWQH19fTIzMyMzMzMCQE5OThQZGalSj59++oksLS2VWi1LW683r2tmZkb9+vVTe+/qjgvi7+9Pw4cPV0qbPn062djY0OPHj4W0QYMGUZ06dZTGufXo0YOmTJkiHC9cuJC+/vpr4TgnJ4csLCwoJydHqfwnT57Q1atXC329+R51UIwWwNzcXOrduze1bdtWSNu7dy/9+uuvFBMTQ/v37ycfHx+qVatWoeNH//nnHwKg0lJia2srfB5TpkyhzMxMMjU1pZMnTyrlGzFiBA0YMEA4vnz5MgGgw4cPC2nt27enQYMGKb0vJiaGANCtW7cKvU/GNK1UE062biXavp1o6FCiESOImjcnsZqxy8ifdGJvT9ShA0V8+ilPKKkktN0CaFAOMWflZWoKpKaW37VLICgoCL1798bx48fxzz//YN++fZg/fz5Wr16NoUOHKuW9ceMGcnJy0LJlSyHN0tISbm5uSvlmz56N2bNnC8cZGRn4559/MHbsWCHt0qVLqFWrlnDcuXNnLFu2DADw7Nkz/Pzzz+jZsyfOnDmD2rVrAwC2b9+OR48e4cSJE0pjIopbL3Vevy6AAluiSiIjIwNisVgpLSYmBgEBAUotmnfu3EH//v1h+tpndufOHfj7+wvHFy5cQK9evYTjy5cvo06dOjAwUP4na2Njo9Qap00hISG4cOGC0hi/nj17Ct97enqiVatWqF27Nn799dcSj4U8c+YM8vLyMHDgQGRlZeHatWtIT09Ht27dlPJlZ2ejadOmwnHDhg3Rpk0brF27Fp06dcK1a9dw/PhxlVYRExMTAIqxmoyVpcDAwCLXBSwwT1CQ8K2kcWPEXryosjahGwA8egQ8eoTQY8cULYgv8xCRYmmaL75A4LvvAnqKjj/e55hxAKhJIhGggUCirIjFYnTr1g3dunXDt99+i5EjR0Imk6kEgMU1ZswYvP/++8LxwIEDERQUpPRLxdnZWek9ZmZmcHV1FY5Xr14NS0tLrFq1Ct999x0AoGnTpoiKisLatWvRokULoSvlbbx5XU2oXr06nj17ppQWHR2NqVOnKqXFxMRgwoQJwnFmZibi4uLg5eUlpF24cEHpfdHR0WjcuLHKNd8MutV5M+gujbFjx2L37t04duwYatasWWA+KysrSCQSXLt2rcA8rq6uEIlEiIuLU0qvV68egFeBWurLP6b27NmDGjVqKOU1NjZWOh4xYgQ+/fRThIWFYd26dahfvz46duyolOfp06cAADs7u8JulTGdJQsNRVBQkMr+xLIxY4DHj4HY2IJ3Obl2DTAwAKysILe2RtCNG8L7Y2NjERQUxAtYVzE8BpAJpFKp0uzOfPXq1YOhoSHOnj0rpCUnJyM+Pl4pn42NDVxdXYWXiYkJ7O3tldLebMF6k0gkgp6eHjIyMoS0+vXr4/Dhw/jtt9/w6aeflrheb8vIyAi5ublF5mvatCkuXbokHKekpODWrVtKrVU3b95EcnKyUlr+bgP5AR4R4cqVK+jWrRvq1KmDOnXqYOzYsfDw8FC55pgxYxAdHV3o682guySICGPHjsWOHTvw119/CTOlC5Kamorr16/DycmpwDy2trbo1q0bfvrpJ7U/b/mkUimMjY1x584dpZ8hV1dXlXGY77//PvT09LBlyxZs2LABw4cPV/lD4cKFC6hZsyaqV69ejDtnTPcEBgYiIiICnp6eEIvF8PT0hFwuR8CyZcC2bcCVK5A0bqzysy+0EhIBz54h9MYN9a2EI0YAK1cCiYkAeCxhZcctgFXQkydP8N5772H48OHw9PREtWrVcO7cOcyfP1+pGzJftWrVEBwcjMmTJ8PGxgb29vaQyWTQ09N769a4rKwsJL78ZfPs2TP89NNPSE1NRd++fZXySSQSHD58GJ06dYKBgQF+/PFHrdbrdXXq1MHp06dx69YtmJubw8bGBnp6qn87+fn5YerUqXj27Bmsra0RExMDfX19pcAtOjoaNjY2Qvd2flr9+vVhbm4OQBEkSiQSnDt3TsjTs2dPtS2Ab9MFnJqaqtRSd/PmTaF++S2GISEh2LJlC3777TdUq1ZN+KwsLS1hYmKCSZMmoW/fvqhduzbu378PmUwGfX19pckX6vz8889o27YtWrRogenTp8PT0xN6eno4e/Ysrly5gubNm6NatWqYNGkSJkyYgLy8PLRr1w7Jyck4ceIELCwsEBwcLJRnbm6O/v37Y+rUqUhJSVHbin38+HF07969VM+KMV1RVHeyTCZT30q4ebOil+rAAcSvX6++lTApCRg9Ghg9GnJ9fQTl5goTUt6mlZC7m3WUVkYWViKVcSHozMxM+vLLL6lZs2ZkaWlJpqam5ObmRt988w2lp6cTUfGWgWnZsiV9+eWXBV6nOJNA8NpA5mrVqpG3tzdt375dqYzX63Hp0iWyt7eniRMnlrpeJVkGhogoLi6OWrduTSYmJoUuA0NE1LJlS1q+fDkRES1dupQaNWqkdF4mk1HXrl2V0kJCQpQmofz2228qExhq1aql8ckLhw8fVjugPDg4WMij7jwA4XPt378/OTk5kZGREdWoUYP69+9P165dK9b179+/T2PHjqW6deuSoaEhmZubU8uWLWnBggXCBJm8vDz68ccfyc3NjQwNDcnOzo78/PyEZWhed/LkSQJAvXr1UjmXkZFBlpaWdOrUqZI/KMYqmFItTQOQl74+kZ4eEVDwAtZmZkRt2hANGUK0di3Rf/8VWRfwpJRSqTI7gRS14wCR4pf4J598QjY2NmRmZkaBgYGUmJiolOf27dvUq1cvMjExITs7O5o0aZLKLMjDhw9T06ZNycjIiOrXr19okFIZA0BNSE1NJUtLS1q9enV5V0VJeddr9+7d5O7uTrm5uaUuY9asWTR79mzhOP+eWOn9/PPP1K1bt/KuBmM6oaCgTAgUY2NJbGBQ8Izj114RL4NFMUCeZmYU0a8fUVwc0YsXRFTIOoi8y0mRqsw6gMXZcWDChAn4/fffsW3bNhw9ehT3799XakbOzc1F7969kZ2djZMnT2L9+vUIDw/HtGnThDw3b95E79690blzZ0RHR2P8+PEYOXIk/vjjjzK934rm33//xS+//ILr168jKipKWIdPXZdxVa5X79698dFHH+HevXulLuPChQtwd3cXjq9cuQKpVKqJ6lVZhoaGWLp0aXlXgzGdUOBYwvw1DD08IJFKVccSikRwq1kTaN4csLGBHEAQgFgAmQBi09IQtH075G5uivVpmzVD/Msxzq8jIpVJYDzesBxoJazUgDd3HEhKSlJZ8y1//a/8bp29e/eSnp6eUqvgsmXLyMLCgrKysoiIaMqUKSrdcv379y9wP1FuAVSIioqiZs2akZmZGVlbW5Ovry+dP3++vKuls/VijLGKrMhWQiqkK1kkEloIC+xKNjYmGjmS6M8/KWL1au4mVqPKdAG/6erVqwSAYmNjiYjo0KFDBEBlu6xatWrRokWLiIjo22+/VWlWvnHjBgGgqKgoIlIsEPvmor5r164lCwsLtfXI/wASEhIoOTlZeGVmZr79TTLGGGM6qqixhGKxWH03sbExUXw80Zo1FOHhIQR9r3+Vv9aNrDZIFInISyIhem04TUREBHl6egrbQla2ADEzM1MpzkhISKgaXcCvy8vLw/jx49G2bVthBmViYiKMjIxgZWWllNfBwUGYmZiYmAgHBweV8/nnCsuTkpKitPTIm1xcXGBpaSm85syZ81b3yBhjjOmyora5k0gk6ruJGzYEGjQAhg9HYGysorvZywtiY2N41q4Nec+eCGjbFmjbFmjQAPFQRH2vIyLExccD+vpAtWqQe3ggKCgIsbGxyMzMFGYlV6au4jlz5ijFGQVtO6opOrkMjLodB8pbQkICLCwshOM3F6JljDHGqpICl5yRyZTyFbV0TaE7nABAaipCL15Uv3bhBx8gsH17wMsLaN8e8rQ0hC5YUCGXnJk6dSomTpwoHKekpGg1CNS5FsD8HQcOHz6stOOAo6MjsrOzkZSUpJT/4cOHcHR0FPI8fPhQ5Xz+ucLyWFhYCDsQqGNhYaH04gCQMcZYVVbkZJJikoWGCgEkoGhFJACyTZuA7duBzz5DvL6+aishgLicHOCvv4AffoA8MBBBgwcj9vx5RSvh+fOKVsKRI4HLl4G8PADFm3BSHpNSjI2NVWINrdJKx3Ip5OXlUUhICDk7O1N8fLzK+fxJIK+vEXflyhW1k0AePnwo5FmxYgVZWFgIY/amTJlCHh4eSmUPGDCAJ4Ewxhhj5aTUaxdWq0bk6kpkZlbwhJP88YYiEUW8HLcoem2sId6YcKIraxdWmUkgH3/8MVlaWtKRI0fowYMHwit/YWIiojFjxlCtWrXor7/+onPnzpGPjw/5+PgI51+8eEEeHh7UvXt3io6Opv3795OdnR1NnTpVyHPjxg0yNTWlyZMn0+XLlyksLIz09fVp//79auvFASBjjDFWvoozK1lsbFzk2oUFBolGRkTe3kReXuRpaal+UkoZr11YZQJAdR8aAKVFmvMXgra2tiZTU1MKCAigBw8eKJVz69Yt6tmzJ5mYmFD16tXp888/V7sQdJMmTcjIyIjq1avHC0EzxhhjOq5UrYT5gduTJ0TbtpFYX7/IIFFcQDwi1tcnWr2a6OhRotu3KWLbtiJnJb/NzOUqEwDqKg4AGWOMMd1X6rULRSLycnUlmjiRyNubPA0MCu9KhmIHFKWu5JdfI0aNIrpwgSg3t9hdyQUFidqOP0REbyzRzZSkpKTA0tISycnJ2h+QyRhjjLFSk8vlmDFjBuLi4uDm5gaZTKY0KUUul6udufz65JUC8/TogQAiID4eXjdvIhbKy9eIAHgCiAYAMzN45eYiNjNTOY9IBM/GjREdE1PotSIiIuDr66vd+EMrYWUlUlVbAG/evEkA6N9//9V42Y8fPyY7Ozu6efPmW5cFgHbs2PHW5ZRW7dq16Ycffii36xelvJ9PcZXm5y04OJj8/f2F4/79+9P333+v+coxxiqVorqSi5OnOOMNC+xKBojEYqIaNcjTxKTA8YZVZi9gVnaGDh0KkUgkvGxtbdGjRw+cP3++TK4/a9Ys+Pv7o06dOirnnjx5gpo1a0IkEikt+TN9+nQ0adKkTOrHyp6LiwsePHggLPxeGt988w1mzZqF5ORkDdaMMVbZFLXAdXHySNzc1C+C7eUFJCUBly5BUqcOlHO8tr5hZiZw7x7iMzLUL4L9xl7J2sABYBXVo0cPPHjwAA8ePMChQ4dgYGCAPn36aP266enpWLNmDUaMGKH2/IgRI+Dp6an1epREbm4u8l6uH6Vp2dnZWim3otHX14ejoyMMDEq/Nr2Hhwfq16+PTZs2abBmjDGmSiaTqa5dmL8ItqUl4O4O2cKFoJfnhDwAZN98A8hkQEAAJCYmqkGiSAQ3NzdoGweA2pCWVvArM7P4ed/cmq6gfKVgbGwMR0dHODo6okmTJvjyyy+RkJCA//77T23+3NxcDB8+HA0bNsSdO3cAAFeuXEG7du0gFoshlUrx559/QiQSYefOnQVed+/evTA2Nkbr1q1Vzi1btgxJSUmYNGmSUnp4eDhCQ0MRExMjtFqGh4cL5x8/foyAgACYmpqiQYMG2LVrV6H3/uzZMwwZMgTW1tYwNTVFz549cfXqVaXrWVlZYdeuXZBKpTA2NsadO3fw6NEj9O3bFyYmJqhbty42b96sUnZSUhJGjhwJOzs7WFhYoEuXLoh5OdYDeNWSuXr1atStWxdisVhtHfPrsHPnTjRo0ABisRh+fn5ISEhQeWb169eHkZER3NzcsHHjxgLvu0uXLhg7dqxS2n///QcjIyMcOnQIAFCnTh3Mnj0bw4cPR7Vq1VCrVi2sXLlS6T2xsbHo0qULTExMYGtri48++gipqanC+aFDh+Ldd9/F7Nmz4eDgACsrK8yYMQMvXrzA5MmTYWNjg5o1a2LdunXCe27dugWRSITo6GgAip+3ESNGoG7dujAxMYGbmxsWL15c4L3l69u3L7Zu3VpkPsYYexvFWQS7wDwzZwLTpwNyOWSbNqkGiWp2U9EKrXQsVyKl6oN/bQyAyqtXL+W8pqYF5+3YUTlv9erq85XQm2Onnj9/TqNHjyZXV1fKfbnx9utjsjIzMykgIICaNm1Kjx49IiLFmotubm7UrVs3io6OpuPHj1PLli2LHHP22WefUY8ePVTSL168SI6OjnT79m06fPgwAaBnz54REVF6ejp9/vnn1KhRI5X1IQFQzZo1acuWLXT16lX67LPPyNzcnJ48eVJgHd555x1yd3enY8eOUXR0NPn5+ZGrqytlZ2cTEdG6devI0NCQ2rRpQydOnKArV65QWloa9ezZk7y8vOjUqVN07tw5atOmDZmYmCiNAfT19aW+ffvS2bNnKT4+nj7//HOytbUV6iOTycjMzIx69OhBUVFRFBMTo7aO+XVo0aIFnTx5ks6dO0ctW7akNm3aCHnkcjkZGhpSWFgYxcXF0cKFC0lfX5/++usvIc/rn8fmzZvJ2tpaWBSdiGjRokVUp04dysvLIyLFmEYbGxsKCwujq1ev0pw5c0hPT4+uXLlCRESpqank5OREgYGBFBsbS4cOHaK6detScHCwUGZwcDBVq1aNQkJC6MqVK7RmzRoCQH5+fjRr1iyKj4+nmTNnkqGhISUkJBCR6hjA7OxsmjZtGp09e5Zu3LhBmzZtIlNTU/rf//6ndJ3Xf46JiPbt20dGRkZK98gYY7qsoPGGvAxMOausAaC+vj6ZmZmRmZkZASAnJyeKjIwU8uT/h3z8+HHq2rUrtWvXjpKSkoTz+/btIwMDA6V1GA8ePFhkAOjv70/Dhw9XSsvMzCRPT0/auHEjEZFKAEikCJzULcIJgL755hvhODU1lQDQvn371F4/Pj6eANCJEyeEtMePH5OJiQn9+uuvRKQIvgBQdHS0kCcuLo4A0JkzZ4S0y5cvEwAhADx+/LjSrjP56tevTytWrBDuw9DQUAikC5Jfh3/++UfleqdPnyYiojZt2tCoUaOU3vfee+9Rr9d+xl7/PDIyMsja2lopiPL09KTp06cLx7Vr16ZBgwYJx3l5eWRvb0/Lli0jIqKVK1eStbU1paamCnn27NlDenp6lJiYSESKn6/atWsLf0wQEbm5uVH79u2F4xcvXpCZmRn98ssvRFS8SSAhISEUFBQkHKsLAGNiYggA3bp1q8ByGGOsIuBJIDrC29sbUqkUYWFhRWdOTS34FRGhnPfRo4Lz7tunnPfWLfX5SqFz586Ijo5GdHQ0zpw5Az8/P/Ts2RO3b99WyjdgwACkpaXhwIEDsLS0FNLj4uLg4uIi7LEMAC1btizyuhkZGSrdnlOnToW7uzsGDRpUqnt5fcygmZkZLCws8OjRI7V5L1++DAMDA7Rq1UpIs7W1hZubGy5fviykGRkZKZWb/77mzZsLaQ0bNoSVlZVwHBMTg9TUVNja2sLc3Fx43bx5E9evXxfy1a5dG3Z2dkXel4GBAby9vVWul1/Py5cvo23btkrvadu2rdJ9vE4sFmPw4MFYu3YtACAqKgoXLlzA0KFDlfK9ft8ikQiOjo7C87x8+TK8vLxgZmamdM28vDylQcuNGjWCnt6rXy8ODg5o3LixcKyvrw9bW9sCPycACAsLQ/PmzWFnZwdzc3OsXLlSGH5QkPz9vNPT0wvNxxhjuiosLAxSqVTp9782lH7EdRVz9uzZ4q/D89p/juWWt8iizODq6iocr169GpaWlli1ahW+++47Ib1Xr17YtGkTTp06hS5durz1datXr45nz54ppf3111+IjY3F9u3bAShmQOXn/frrrxEaGlpomYaGhkrHIpHorSdtmJiYqMzwKkpqaiqcnJxw5MgRlXOvB4pmGvwcS2rkyJFo0qQJ7t69i3Xr1qFLly6oXbu2Uh5NPE91ZZSk3K1bt2LSpElYuHAhfHx8UK1aNSxYsACnT58u9LpPnz4FgGIF2IwxpotCQkIQEhIirEOsLdwCyAAo/jPW09NDxhsTTz7++GPMnTsX77zzDo4ePSqku7m5ISEhAQ8fPhTSzp49W+R1mjZtikuXLimlRUREICYmRmiRXL16NQDg+PHjCAkJAaBokcvNzS31/eVzd3fHixcvlAKJJ0+eIC4uDlKptMD3NWzYEC9evEBkZKSQFhcXp7RUTbNmzZCYmAgDAwO4uroqvapXr17iur548QLnzp1TuZ67u7twLydOnFB6z4kTJwq9j8aNG6NFixZYtWoVtmzZguHDh5eoTu7u7oiJiUHaa5OPTpw4AT09PY3OWjtx4gTatGmDTz75BE2bNoWrq6tSK2pBLly4gJo1a5bqeTPGWFXCAWAVlZWVhcTERCQmJuLy5cv49NNPkZqair59+6rk/fTTT/Hdd9+hT58++PvvvwEA3bp1Q/369REcHIzz58/jxIkT+OabbwCg0JYzPz8/XLx4UakVsH79+vDw8BBedevWBaAINuzt7QEoZqfevHkT0dHRePz4MbKyskp13w0aNIC/vz9GjRqFv//+GzExMRg0aBBq1KgBf3//At/n5uaGHj16YPTo0Th9+jQiIyMxcuRIocsRAHx9feHj44N3330XBw4cwK1bt3Dy5El8/fXXSoFccRkaGuLTTz8Vrjd06FC0bt1a6GqfPHkywsPDsWzZMly9ehWLFi2CXC5XmUX9ppEjR2Lu3LkgIrXrXxVm4MCBEIvFCA4OxoULF3D48GF8+umnGDx4MBwcHEp8jwVp0KABzp07hz/++APx8fH49ttvi/UHxvHjx9G9e3eN1YMxxiorDgCrqP3798PJyQlOTk5o1aoVzp49i23btqFTp05q848fPx6hoaHo1asXTp48CX19fezcuROpqanw9vbGyJEj8fXXXwNAgUubAIoWqGbNmuHXX38tUX2DgoLQo0cPdO7cGXZ2dvjll19K9P7XrVu3Ds2bN0efPn3g4+MDIsLevXtVuijVvc/Z2RkdO3ZEYGAgPvroIyFABRSB7969e9GhQwcMGzYMEokEH3zwAW7fvl2q4MjU1BRffPEFPvzwQ7Rt2xbm5ub43//+J5x/9913sXjxYnz//fdo1KgRVqxYgXXr1hX4GeYbMGAADAwMMGDAgEI/q4Lq9Mcff+Dp06fw9vZGv3790LVrV/z0008lvr/CjB49GoGBgejfvz9atWqFJ0+e4JNPPin0PZmZmdi5cydGjRql0bowxlhlxHsBF4H3Ai6+EydOoF27drh27Rrq169fYL49e/Zg8uTJuHDhgtJEAfZKeHg4xo8fr9TFrCm3bt1C/fr1cfbsWTRr1kzj5ZeXZcuWYceOHThw4EB5V4Uxxt6atuMPngTCSm3Hjh0wNzdHgwYNcO3aNYwbNw5t27YtNPgDgN69e+Pq1au4d+8eXFxcyqi2LCcnB0+ePME333yD1q1bV6rgD1B0mS9durS8q8EYYxUCB4Cs1J4/f44vvvgCd+7cQfXq1eHr64uFCxcW673jx4/XbuWYihMnTqBz586QSCTCjOvKZOTIkeVdBcYYqzC4C7gI3AXMGGOMsbKm7fiDB2AxxhhjjFUxHAAyxhhjjFUxHAAyxhhjjFUxHAAyxhhjjFUxHAAWk7e3N6RSKcLCwjRedlpaGkQiEUQikdIWW4wxxhirWsLCwiCVSuHt7a3V6/As4CKUxSzgtLQ0mJubAwBSU1NhZmamleswxhhjrGLgWcBMK4YOHSq0OhoaGsLBwQHdunXD2rVrkZeXV+xywsPDYWVlpb2KMsYYY0zjOACswnr06IEHDx7g1q1b2LdvHzp37oxx48ahT58+ePHiRXlXjzHGGGNawgGgDvjtt9+E71u3bg25XF4m1zU2NoajoyNq1KiBZs2a4auvvsJvv/2Gffv2ITw8HACwaNEiNG7cGGZmZnBxccEnn3yC1NRUAMCRI0cwbNgwJCcnC62J06dPBwBs3LgRLVq0QLVq1eDo6IgPP/wQjx49KpP7YowxxljhOAAsZ3K5HAMHDhSOL168iKCgoDILAt/UpUsXeHl5CdfX09PDkiVLcPHiRaxfvx5//fUXpkyZAgBo06YNfvzxR1hYWODBgwd48OABJk2aBECx7+zMmTMRExODnTt34tatWxg6dGi53BNjjDHGlPEkkCJoexCml5cXYmNj8frHIBKJ4OnpiejoaI1fL9/QoUORlJSEnTt3qpz74IMPcP78eVy6dEnl3Pbt2zFmzBg8fvwYgGIM4Pjx45GUlFTo9c6dOwdvb288f/5cmPDCGGOMMfV4EkglFx8fjzdjcCJCXFxcOdVIcX2RSAQA+PPPP9G1a1fUqFED1apVw+DBg/HkyROkp6cXWkZkZCT69u2LWrVqoVq1aujYsSMA4M6dO1qvP2OMMcYKxwFgOZNIJEKwlU8kEsHNza2cagRcvnwZdevWxa1bt9CnTx94enoiIiICkZGRwjqI2dnZBb4/LS0Nfn5+sLCwwObNm3H27Fns2LGjyPcxxhhjrGxwAFjOZDKZSvcvEUEmk5VLff766y/ExsYiKCgIkZGRyMvLw8KFC9G6dWtIJBLcv39fKb+RkRFyc3OV0q5cuYInT55g7ty5aN++PRo2bMgTQBhjjDEdwgFgOQsMDMTmzZuFYw8PD8jlcgQEBGj92llZWUhMTMS9e/cQFRWF2bNnw9/fH3369MGQIUPg6uqKnJwcLF26FDdu3MDGjRuxfPlypTLq1KmD1NRUHDp0CI8fP0Z6ejpq1aoFIyMj4X27du3CzJkztX4/jDHGGCseDgB1gL+/v/D9qVOnyiT4A4D9+/fDyckJderUQY8ePXD48GEsWbIEv/32G/T19eHl5YVFixZh3rx58PDwwObNmzFnzhylMtq0aYMxY8agf//+sLOzw/z582FnZ4fw8HBs27YNUqkUc+fOxffff18m98QYY4yxovEs4CLwVnCMMcYYK2s8C5gxxhhjjGkUB4DF5O3tDalUKsyC1SQzMzMQEYiIW/8YY4yxKiwsLAxSqRTe3t5avQ53ARehLLqAGWOMMcZex13AjDHGGGNMozgAZIwxxhirYjgAZIwxxhirYjgAZIwxxhirYjgAZIwxxhirYjgAZIwxxhirYjgALEJWVpbSV6Y9WVlZmD59Oj9rLePnXHb4WZcdftZlg59z2dF2/MHrABbh7t27cHFxQUJCAmrWrFne1anUeM3FssHPuezwsy47/KzLBj/nsqPt+INbAMuQpnYRqazlaJKu3ZuulaMpmqyPrj2jyvqsK2s5mqRr96Zr5WgK//4oZ8QKlZCQQAAoISHhrctyd3fXQI0qbznJyckEgJKTk9+6LF27N10qRxefsybL0qVydPFZV9Zy+FmXTTm6+Jw1WZYulaPJ+EMdg3KMPSuE3NxcAMD9+/ffuqzs7GzcvXuXyynA8+fPAQD37t1DSkqKTtSpMpaji89Zk2XpUjm6+Kwrazn8rMumHF18zposS5fKyY878uMQTeMxgEU4fPgwunTpUt7VYIwxxlgV9Ndff6Fz584aL5dbAItQr149AMClS5dQo0aNtyrL29sbZ8+efes6VdZyUlJShAGvbzu4WNfuTZfK0cXnrMmydKkcXXzWlbUcftZlU44uPmdNlqVL5dy7dw9SqVSIQzSNA8Ai6OvrAwDeeecdGBoaIiQkBCEhIaUuSxOzpiprOfksLCzeujxduzddKwfQreesybJ0rRxAt551ZS0nHz9r7ZcD6NZz1mRZulBOWFgYwsLCkJOTI5SlDRwAFtPhw4ffehp2aQPHqlKOJunavelaOZqiyfro2jOqrM+6spajSbp2b7pWjqbw74+C3xsSEiIsA6MtPAawCLwOYNnh9aXKhq4/57S0NJibmwMAUlNTYWZmVs41Kj1df9aVCT/rssHPuezwOoDlzNjYWOkr0x5jY2PIZDJ+1lrGz7ns8LMuO/ysywY/57Kj7fiDWwCLwH/tMFa2KlMLIGOMlZa24w9uAWSMMcYYq2IqVAC4bNkyeHp6CrOPfHx8sG/fPuF8ZmYmQkJCYGtrC3NzcwQFBeHhw4dKZdy5cwe9e/eGqakp7O3tMXnyZLx48aKsb4UxxhhjrNxUqACwZs2amDt3LiIjI3Hu3Dl06dIF/v7+uHjxIgBgwoQJ+P3337Ft2zYcPXoU9+/fR2BgoPD+3Nxc9O7dG9nZ2Th58iTWr1+P8PBwTJs2rbxuiTH2ht9++034vnXr1pDL5eVYG8YYq5wq/BhAGxsbLFiwAP369YOdnR22bNmCfv36AQCuXLkCd3d3nDp1Cq1bt8a+ffvQp08f3L9/Hw4ODgCA5cuX44svvsB///0HIyMjlfJ5DCBjZUculyMoKEg4FolEICJEREQo/THHGGOVHY8BLEBubi62bt2KtLQ0+Pj4IDIyEjk5OfD19RXyNGzYELVq1cKpU6cAAKdOnULjxo2F4A8A/Pz8kJKSIrQiFiQlJUXplZWVpZ0bY6wKCw0NhUgkEo6JCCKRCDNmzFDKl5aWBpFIBJFIhLS0tLKuJmOMaVxWVpZKrKFNFS4AjI2Nhbm5OYyNjTFmzBjs2LEDUqkUiYmJMDIygpWVlVJ+BwcHJCYmAgASExOVgr/88/nnCuPi4gJLS0vhNWfOHM3dFGMMABAfH483OyWICHFxceVUI8YYKxtz5sxRijO0uQg0UAF3AnFzc0N0dDSSk5Oxfft2BAcH4+jRo1q/7pv7HvIaSIxpnkQiQWxsrFIQKBKJ4ObmVo61Yowx7Zs6dSomTpwoHOfvu6wtpQoAX6/g60QiEcRiMVxdXeHv7w8bG5u3qpw6RkZGcHV1BQA0b94cZ8+exeLFi9G/f39kZ2cjKSlJqRXw4cOHcHR0BAA4OjrizJkzSuXlzxLOz1MQTex7yBgrnGziRAQNHQoRAAIUX4kgA4CEBEDLfxEzxlh5MTY2LtPGpVJ1Af/7779Ys2YNVq5ciaNHj+Lo0aNYtWoV1qxZg0OHDmHixIlwdXXFpUuXNF1fFXl5ecjKykLz5s1haGiIQ4cOCefi4uJw584d+Pj4AAB8fHwQGxuLR48eCXkOHjwICwsLSKVSrdeVMVa4wD//RAQAyctjDycnyA0NERATA3h4AOHhQMWet8YYYzqhVAGgv78/fH19cf/+fURGRiIyMhJ3795Ft27dMGDAANy7dw8dOnTAhAkTNFrZqVOn4tixY7h16xZiY2MxdepUHDlyBAMHDoSlpSVGjBiBiRMn4vDhw4iMjMSwYcPg4+OD1q1bAwC6d+8OqVSKwYMHIyYmBn/88Qe++eYbhISEcJcuY+Vt715g0yYE6Okhv6391NWrCIiNBVq3BlJSgGHDAH9/iN5Y35MxxlgJUSk4OzvTxYsXVdIvXLhAzs7OREQUGRlJtra2pSm+QMOHD6fatWuTkZER2dnZUdeuXenAgQPC+YyMDPrkk0/I2tqaTE1NKSAggB48eKBUxq1bt6hnz55kYmJC1atXp88//5xycnIKvGZycjIBoOTkZI3eC2PsNcnJRC4uRABlf/opQdEDTKmpqYrzOTlEc+YQGRkRAbTd3JwsX+bx8PCgiIiI8q0/Y4xpmLbjj1KtA2hubo7du3ejU6dOSulHjhxB37598fz5c9y4cQNNmjTR+jRmbeN1AEuG93FlpRISAvz8M1CvHtJOnYL5y9n5Kj9DsbGQ+/sj6ObNV+MEea1AxlglpJPrAPr7+2P48OHYsWMH7t69i7t372LHjh0YMWIE3n33XQDAmTNnIJFICi+IMcaOH1cEfwCwahVQ2B8NjRsj1NxcCP6AgtcKZIwxVrBSzQJesWIFJkyYgA8++EDYR9fAwADBwcFYtGgRAMUizKtXr9ZcTRljFZ5KC7GeHjBihOLkqFFAly5AEQs7x1+9ije7LYjXCmSMsRJ5q63gUlNTcePGDQBAvXr1hF/slUl+E6xEIoG+vj5CQkIQEhJS3tXSWdwFzAqj8vMxcyYwbx7g7AxcvAi8sZC7Ol5eXmrXCvT09ER0dLSWas4YY2UjLCwMYWFhyM3NRXx8vNa6gEsVAB4+fBidO3dWey4sLKxSBUg8BrBkOABkhXn95yP9+HGYdOoE5OYCv/0GvPNOscpQ2S8Yiu5g+YYNCBg8WPOVZoyxcqCTYwADAwMRGRmpkr548WJMnTr1rSvFGKvcDAAYffKJIvjr37/YwR+g+P2zefNm4bixsTHkAALKYN1RxhirLEoVAC5YsAA9e/bElStXhLSFCxdi2rRp2LNnj8YqxzQjLS0NIpEIIpEIaUWMr2KsLEwGoB8bC9jaAkuWlPj9/v7+wvenN2xAAAAsXgw8eKCxOjLGWGVWqgBw5MiRmDRpEnx9fXHr1i3MmzcPM2bMwN69e9G+fXtN15FVIL/99pvwfevWrSGXy8uxNsXDAXLZaggotnYDFEGbvf1blZfbqxfg4wNkZAAzZ75t9RhjrEooVQAIAFOmTMHAgQPRokULzJ07F3/88Qfatm2rybpVapUx6JDL5Rg4cKBwfPHiRQQFBVWIIJCVjd927kQOACsAnubmkIvFb1+oSATMnav4ftUq4Nq1ty+TMcYquWJPAllSQDfN999/jw4dOqBly5ZC2meffaaZ2ukAbQ3CLMvJEmV1rYo6O5MnrpSN/MkbmljAWe1n1rMnsH8/MGAAsGULf66MsQpN25NAih0A1q1bt3gFikTC0jCVAQeAxWdiYoLMzEyVdLFYjIyMDK1cUxM4UCgbXo0aIfbSJaU1/DT6B8K//wLNmgnfpzVowJ8rY6zC0nYAWOyFoG/evKnxi7PKRSKRqG0BdHNzK8daMV0RHxen3QWcmzYFPvgA2LoV+Ppr4NdfNVMuY4xVQqUeA1jVeHt7QyqVIiwsrLyrorNkMplK8EdEkMlkhbyLVQknTkCSmwvRG8ka/wNhxgxAXx/Yuxd6J05ortwiVMYxvYyx8hEWFgapVApvb2+tXqfEAeDVq1cREREhtAju2bMHHTp0gLe3N2bNmoW32FhEp509exaXLl2qkItcl9XM3DfXZ/Pw8IBcLkdAQIBWrqcpFXHmcoWSkwN8/DFkgEr3r8b/QGjQABg5EgBgNH265splFQoH5KwiCwkJwaVLl3D27FntXohKQC6Xk4GBARkZGZGxsTGtX7+exGIx9ejRg3r37k0GBgY0d+7ckhSp85KTkwkAJScna7Tc1NRUguL/Q0pNTdVo2a+LiIgQrgOARCIRAaCIiAitXK+s7ktTyvr5VFRv9bkuXEgEENna0ubly4VyGjduTHK5XPOVvXePSCymCICqvbyWh4eHVj/TivZzn6+i1rsolfW+WNWirfgjX4kCwObNm9NXX31FeXl5tHbtWjIxMaEffvhBOL9ixQpq2LChputYrrT1AWzevFn4BaXN/5w8PT2FoOb1IMfLy0sr19PUL96y+gVe1s+noirO56E2T0ICkbm5IgBcs6bs/vB5913FZ8l/+BSq1J+rjquIdWbsTdoOAEvUBRwXF4fhw4dDJBIhODgY2dnZ8PX1Fc53794dt2/fLkmRVVJZrpcXHx+v0i1Pmhx4X8FV5OdTIbq5JkwAUlOBNm2AoUPL7LKhV68Ky80Ais9UJBJhxowZZVYHxhjTZSUKANPS0lCtWjXFG/X0YGJiAlNTU+G8iYkJsrKyNFvDSig0NBQi0avh8Nr8z0kikShdC+CZua/j56NF+/cD27crJmUsWwbold2cs/jr17U74/gNVX0cqa79MVLVPw9do2s/H0yhRL+R8z/Ago5Z8RS31UkT/2h4Zm7h+PloSWYmMHas4vtx4wBPzzK9fFkG9hV5B5zKGChV1M+DgyRW5krSXywSicjKyoqsra3J2tqaRCIRWVpaCsdWVlb0//bOOy6Kq+vjv13aCoam6IKCNSBErEFssSQENQUVS2JFgxqNSYxRrDFgnsTEEs3zvFjeWDAqxkTBEltiQY2vikbgMRrFElSigthAFFZg7/sH7LjLFmaX2d1ZON/PZz6wM2fOnHvu3Dtn7r1zr1QqFbCH2vqYow+e77gzocaxqI83NNvA+wrENgaQjx5L+oevTWIbm2WsPYrZs8vH/TVqxFhBgVlt04XWxz0Vf5O3bRP8WrY6jpTvB1C2di/aan7U5HGLNTlt5sTcYwB5TwQNAAkJCcaIE3qIjY3FoEGDuN/mbnXq378/9//JkyfNuiKCi4uLzU0FZEn/2CqVW4rmz5+vd/m2lgAcli4t//Hdd0DFsBFLopqSaMSIEZACCAYQC2Dgr78CkZHl6wcLhK2OI1UNRVHZztSGoqjnrTF5LwZsNT9qMrZ2D9UaqhtBZmdns7KysuqqES2qCNzf358FBgay+Ph4QfTyaXUSW2uaJbFk2i3tH7G1AFalh29LkeqedgBYG4AltWvHmFJpsl3VRT1dRevXMyaVlrdKTplilF01tcVJJpNp2KzaZDIZJ2Ns3gPWn3LHVvODrw9trT6nqbaMJz4+ngUGBjJ/f3+ztgBWOwB84YUX2LVr14SwRZRYcx7AmhwEVYVQDxQxpl0omyzlIz4PVH1drtas5LXSlZBQHgACjM2ZY7qeSuh7wJl7KEF14ZOvJuW9lafcscX8MMaHYqurq8JWA3IxDGsQ1TyAuqhbty4FgCZgyVYgMQZBhhDygSJG/whxPUv6iE9LkRgreZ3pWr78eRD41Vem66mEpceRCgGfQEmMec/nxUeo/BDjfKS2Vp/zuYfECAWAPKAA0DTEFtyJqcIQ8oEiRv8IcT1L+qjKa126xGRSqegqeb3pWrSICwKLFy4U3f1hSaoKlPjcZ5Z8wAv54QofLJWvxviQTwAspvSL8eWQD5bshdJ3LdEHgAsWLGAPHz4UwBRxYkoGiC1ws7UHmJAPFDH6x9gWDF0ylvSR3paixETGZs5kzMGBtcHzbl+xVPIG0xUbywWB40R2f1gSIbpTLfmAt/QMCpbKe77psvSX20KkjW+XvJjKkCV7WAxdS/QBIGOMlZaWsvT0dPbgwQMh1IkKawaAtvamJxRCPlAs+RbHBz4VCx8ZS7eSarQUtW7NkqdNY6xxYy6ISurQgVclLxqUSsamT2cMYGUAGyGScsgXS5b7qloJLTnmju+Lj60FgHx9yLfcC3XPCuVHS370KASWfAYZupYoA8ApU6awNWvWMMbKg79u3boxiUTCXFxcWEpKipD2mY34+HjWpEkT5uTkxDp16sRSU1N1ylkrAKypXR18EOqBYumxhHww98B7Ux66xjwsAgFW2rPn83F0zZsz9ssvWnpsYhycUsmejR/PpaWsVSvGJk5kbPPm8jWMmXi/YBRb74Cl8l6M4w0t6UNLf7kthm5QayBUDwufvDB0LVEGgI0aNWJnzpxhjDG2fft25u3tzTIzM9lnn33GunbtKqiB5mDLli3M0dGRrVu3jl24cIGNHz+eubu7s9zcXC1ZnRlQWKh/KyrSvNlzc3XKFebmMlnlAqF2vE3r1oYruidPOD3OAHNWv9aTJ5qJqJDVuVXIqmyWGbCZVS64T58a9oU6RUVGySauW8d1KQa/9BJL3rz5uaz6FB7FxXp16vVhcLDuc9SnM1IotPJLw88GZLW20lJOVObkpLuwOzlxsnorBJVMxZa4YcPzB0br1po+qryVlDy3t6SE25+kVulqVFBr1zKWnl4e3K1YwRSzZrGNAHumCvycnBibO1djkufC/Hzte1F9Uyie21Baathn6rJlZcLJFhc/l1UqWeKaNaxRxX3fBmBJqvQBjDVtytpUTHhf5T1UVKR5D1dRR/CWffpUp6zOcl9JtvDuXf35YUQdUXj3rubDW0e517BHvTwLXEfoffGpfO/zrCMMyeotG6qHd0W5V6+rWgcFsaTERK36pPDBA8Nlo7SUqepgBwN1MJ86Ta9Mmzb801ZSol8mMZGxZ8+e+0ytPtG5qcvqKfcGr6XvPlPvcVQr9zrLhpF1hN4XjcrlvlJ9wju/2rblJWPuAFDCmPGz9spkMly9ehWNGzfGhAkT4OzsjO+++w5ZWVlo27YtCgoKjFVpUUJDQxESEoL4+HgAgFKphK+vLz766CPMmjVLQ7agoABubm7Izs6GnZ0dZvj4YKMh5d7eKO3SBUkVyw4NtbODpKxMt6yXF9C79/PfO3cCFWsp1wFQrOMUmVSKosGDgd27gadPdet1dQX69n3+e/9+QF+eODsDb72F0tJSJCUn43UAnvrS5uQEqE2ajJQUIC9Pt6ydHaA22TV+/x24c0efZmDo0Of/nzgB/POPftnISMDeHqWlpchOTkYzPWJ6fQigSNcJb74JqCaB/u9/AbWJY5MBzAdwGYA/gNh27RDp719+8Px54K+/9NsbFgZ4lnu17a5d+LO4GOqFTgKgDYAMAOjVC23T0/Fnfr5+mQrKunTBtpMnAQCDO3aE3dmz+m3o3Bnw8QGUynLf/vFHuT0A/gSqvJZe4uOByZMBAEX79qHOG2/ol120CIiJKf//zBmgUyf9srGxQFxc+f8XLgCtW+uXnT4dWLy4/P/r14Fm+u4IAB98ACxfDgBITkjAoPfegwTl6Vf9TQKgmqKW9z3UoEF5ehwdyyeZTkrSa8JtAA0iI2FvXzEPf1ISYEIdoYWHB/D669xPtns3JALUEczZGT9X6BkUGQn7lBTg4UOdssUA7NXTZoY6YuuNG4hNTUUWgABUTPBdWbaijgAAnD5dfl/oIyICkMnK/z97Frh2DYCBsuHmhow+fYD//hfJmZkYBOi+h/r0AdzcAADKc+cgvXRJvw1hYSh1dUVScjICALTTI5YM6LxespoP9N6zEgmKKh73VaYtOxttT57UXzeEhDwvZ7dvA8eP609b+/bAiy+W/3/3LnDkiJaIKfWQslUrSFVLTObnA7/+qt+GgACgbdvy/588Afbs0S/bogWSGzbEoBMnDPoZANC06fN6rLQUqLTcoMHnuFIJQEeeVkzOvmnTJvTs2RO+vr7Iz8+Hq6urfptNxZSo0c/Pj/3666+stLSU+fr6st27dzPGGDt//jxzd3cXLjw1AwqFgtnZ2bHt27dr7B89ejSLiIjQkldF4KptsXoLgRk3nYPqAdbWTNdLqrimzpYQG92E8mGS2rnqf03xkT5dyUbKCLWpWqErbzKAsYYNGWvWzLAOtYnRn+7da1h20aLnBev0acOysbHPZc+fNyw7ffpz2awsw7IffMCJtgkKMnx/hIWxNg0aWLQcGnMf1bTyKrbNYNmo2PjWMULmV1KFflnF38r1Ah+b+KSNj4wlfc3Xj0L5uio/89mMuT9U12qgww+iWApOxdixYzF06FB4e3tDIpEgLCwMAJCamopWrVqZotJi3Lt3D2VlZWjYsKHG/oYNG+KSgTc0VQvgEB8f3AHwxfz5cHF21haUSqFQKjG9oqVjyVdfwcnRUbdSiaS8xUCF2pt97J9/YtCGDWqiFcvFRUeXv8k8K++M+yk9HR9u3oyHAILkcsSFhyOyTRtNvRWy+mxIvngRg9au5Xb9ifI3kqTRoxEZHKwp7+TET29l2ZKS8tYnAIpnzzBn7lwAwAKVfyrJKoqLtWVUVLSyKBQKzIqJgVSXDAz4cPRooHK6AMDBAZBKy/8vLeVaZeYvXQrJnTvc2ymr0PWFtzciZ83SkNWJmt7I0lJs+eMPfLR5Mx6gPM/mh4djoMoeBwdESqXY8scf+HDDBi5fNWQqUJSVYfqsWWAAvl24EE52dvptcHR8bodSWb5JpfD/8kv8efs2mFo+SiQSBAQHl7eClpbqb3FS6a1A2a0bVIvo3c3N1V5Sz8Hh+f8dOgCFhfr1qssGBvKX9fMzLGv/vMq7/PffqHz3MgCZTk7A/fuAvT1i9+zRvWzjnDnAyy+X+0ahKPcTY+VlAijfp142GINCocDsuXPBAHyzeDGcVPd8Jf9qlI8FC+BUt67aQQWSK93XXHmNikJkx47c/p9OndKuG1T3UOW6p4o6grcsUK06wmC5r5CtqtwDeN4SW6FXVffoRI+s/9Kl+FOt3AMVZcPHB5g5EygtxeUZM8BKSzXUMQCZ9vbAkiWAVIrk//636vpVT92jizeVSqTMnIlXACz55hutcq+33hs7FggK4pe2sjL4L1qkW0YuB2bMKG/BBcptreQDDeztn8sqleU+roRee+Ry4NNPAUD/fR8dXe7HkhL9MmPGILJDB86Gn06f1l82pFLAwQFvKhRIiYnBa9B/nylKSzF99mwAwJJFi+BUaYlJg89xtVjpzWfP8H9z5+IVALH79sGxa1cA5T2Qvr6++jxbfUyNHLdu3cqWLl3KsisGSzPG2Pr169mOHTuqGZOal1u3bjEA7MSJExr7Y2JiWKdOnbTk1fvg+X4IINRgVlO/vjP2epYeWG1rg9OFnufMFqdoqImY8kWlJSYVrkpOjKtzCIUlv0y11BQ4Qtevlvpy25LTt1jS10J+YGnuSclF+RHIDz/8wIrVBz9WoFAo2A8//FBto8yJqV3AfANAS355KlTFYumZ2sX2BWNVWKMCF9sUDTURS89PJtRar2JcnUMoLDkbgaUCKWushCHUy7Elp2+p6lp8/GjJsmGJmTpEGQBKpVKdX8zeu3ePSaXSahtlbjp16sQ+/PBD7ndZWRlr1KgR+/rrr7VkjQ0ALTk/m1AVixhbAIXSI8ZWMksGd0LZU1Ox1APOmBfD6swbpsJWl9+y5FQgluplsEYwbq1eFmuvSW7JsmGJfBVlACiRSNjdu3e19mdkZDAPD49qG2VutmzZwpycnNj69evZX3/9xSZMmMDc3d1ZTk6OlqyxAaAlV2gw95uMuVqCLNmiIsZWMgoAxYPYWpGrM1G4tVbnEApLd1tbao4/sQ6zMGfaLLlihrHd1rbUcyaqALBdu3asffv2TCqVsuDgYNa+fXtua9OmDXvhhRfYkCFDzGKo0PzP//wP8/PzY46OjqxTp07s1KlTOuXE3AJorsmAzd0NKMQboxgnyhZqjJeQdvOBAkDz5wXfh4VQ4xLFGnQYwpJBq6Unia+pwyysOTG3qd3WQpWNWtcCGBcXx+Li4phEImHTp0/nfsfFxbEFCxawzZs3M4X6hIs1APUM4BO4CFnxWrJisXZ3qrk+XBFjACg2bNVuIbBU2vner0IudWZrQYclu61tddiL2LDV8eNClA1LvGSJKgBUsX79elZUeTb7GooqAzZu3Mg7cLG1oMyS17Kl5nd1bDUAtFW7LYWl0i70Wq+2mK+WGtLCBzHWH7aIrQbSYhwapAtRBoAqFAoFy87OZjdu3NDYahKqDHjppZd43+i2FpRZ8lq2NABXHVt84DImTptqK3weFkJ+mSy2vBfTWDlbHCMpRmxxqIGQmLuMiTIAvHz5MuvevTuTSqUam0QisYmvgI2h8kogfAIXWxy/ZWstgGL8cEVsD1zGxGlTbYVvXlh63kFLIaZu69oeuAiJrQ01sAXi4+NZYGAg8/f3F18A2LVrV9ajRw+2d+9elp6ezjIyMjS2moQ1WwD5YGsBoK1+uMIHsT1wGROnTbUVIT8SskXE9hIltvrDVqmp96sYMHcLYMXaM8aRkZGB//3f/0W/fv3Qrl07tG3bVmOricyaNUtruSzGGGJjY61ole0RGRmJxMRE7nfr1q2RnJyMgQO1lnOvkv79+3P/nzx50iQdQrJz507u/86dOyO50sLgBFGbEVv5EFv9QRCWxqQAMCgoCPfu3RPaFlETEREhWOAiRlxcXMDKW4S113AVmJpY8SYnJ2PEiBHc7wsXLmDQoEFWf8gRhBig8lFzseSzgxAWkwLAhQsXYsaMGThy5Aju37+PgoICja2mIsbAxZKFT2xv8GJi/vz5kKgtBM4Yg0QiwRdffGFFqyjPxATfsloTH6h8y0dNTDtBiBUJU+/X5IlUWh43qhdo4HmhLisrE8Y6EVBQUAA3Nzfk5+fDzs4OdevWBQAUFhbqraCePHnCS86WSE5OxqBBg7jfqi7wpKQkREZGGqVLKP+Iyc916tRBcXGx1n6ZTIaioiIrWCRsnhFEdRBj+SAIsaMef7i6ugqu396Uk1JSUoS2o0aheoutSaje4FXpUn+Dp2AC8Pf3x59//qk1TjQgIMBqNlGeEWJBjOWDIGo7JgWAPXv2FNoOQuRcvnxZK6hljCEzM9NKFomL2NhYna1t1vxIiPKMEAtiLB8EUdsxaQwgAPz+++8YOXIkunbtilu3bgEANm7ciOPHjwtmHCEe/P39tbr8rf0GL6bxQkJ+3SwUYswzonYixvJBELUdkwLApKQk9OnTB3Xq1EFaWhoUCgUAID8/HwsWLBDUQLEQEhKCjh07WtsMqxEbGyvYNDhiCtyERGwfCQmZZwRRXcRWPghCrCxfvhxBQUEICQkx63VMCgC//PJLrFq1CqtXr4aDgwO3v1u3bkhLSxPMODFx5swZnD171tpmWA16g7c9KM8IgiBsj8mTJ+Ovv/7CmTNnzHodk8YAZmZmokePHlr73dzc8OjRo+raJFpq4scdxlD5Db4mtd7VVCjPCLFQ2+tPghAbJrUAyuVyXL16VWv/8ePH0bx582obRRAEQRAEQZgPkwLA8ePHY8qUKUhNTYVEIsHt27eRmJiI6dOnY9KkSULbSBAEQRAEQQiISV3As2bNglKpxGuvvYanT5+iR48ecHJywvTp0/HRRx8JbSNB2ATUxUUQBEHYCiatBKLi2bNnuHr1KgoLCxEUFMStylCTMPdM3LaEmFbeIAiCIIiajLnjD5O6gDdt2oSnT5/C0dERQUFB6NSpU40M/giCIAiCIGoiJgWAU6dORYMGDTB8+HDs3bu3Rq39S+inps7fRxAEQRC1DZMCwDt37mDLli2QSCQYOnQovL29MXnyZJw4cUJo+0RDSEgIgoKCsHz5cmubQhAEQRBEDcVSE0FXawwgADx9+hTbt2/H5s2bcfDgQTRu3BjXrl0Tyj6rQ2MACYIgCIKwNOaOP0z6ClgdZ2dn9OnTBw8fPsSNGzdw8eJFIewiCIIgCIIgzIRJXcBAectfYmIi3njjDTRq1AjfffcdBg4ciAsXLghpH0EQBEEQBCEwJrUAvvvuu9i9ezecnZ0xdOhQzJs3D126dBHaNoIgCIIgCMIMmNQCaGdnh59//hl37txBfHy8xYK/iIgI+Pn5QSaTwdvbG6NGjcLt27c1ZM6dO4dXXnkFMpkMvr6+WLRokZaerVu3olWrVpDJZAgODsbevXstYj9BEARBEIQYMCoAfOONN5Cfn891/S5evBiPHj3ijt+/fx9BQUFC28jRu3dv/Pzzz8jMzERSUhKuXbuGwYMHc8cLCgoQHh6OJk2a4OzZs1i8eDHi4uLw/fffczInTpzAsGHDEB0djfT0dAwYMAADBgzA+fPnzWY3QRAEQRCEmDDqK2A7OzvcuXMHDRo0AAC4uroiIyMDzZs3BwDk5ubCx8fHYvMC7tq1CwMGDIBCoYCDgwNWrlyJuXPnIicnB46OjgDKl63bsWMHLl26BAB455138OTJE+zevZvT07lzZ7Rr1w6rVq3SugZ9BUwQBEEQhKUR1UoglWNFa657+uDBAyQmJqJr165wcHAAAJw8eRI9evTggj8A6NOnDzIzM/Hw4UNOJiwsTENXnz59cPLkSYPXKygo0NgUCoXAKSIIgiAIoraiUCi0Yg1zYvJXwNZi5syZcHFxQb169XDz5k3s3LmTO5aTk4OGDRtqyKt+5+TkGJRRHdeHr68v3NzcuO3rr78WIjkEQRAEQRD4+uuvNeIMX19fs17PqABQIpFAIpFo7asOs2bN4vTq21TdtwAQExOD9PR0/Pbbb7Czs8Po0aMt0hKZnZ2N/Px8bps9e7bZr0kQBEEQRO1g9uzZGnFGdna2Wa9n1DQwjDGMGTMGTk5OAIDi4mJMnDiRWxfWlG7RadOmYcyYMQZlVGMMAaB+/fqoX78+/P39ERgYCF9fX5w6dQpdunSBXC5Hbm6uxrmq33K5nPurS0Z1XB+urq40BpAgCIIgCLPg5OTExVeWwKgAMCoqSuP3yJEjtWRGjx5tlAFeXl7w8vIy6hwVSqUSwPPAs0uXLpg7dy5KSkq4cYEHDhxAQEAAPDw8OJlDhw7hk08+4fQcOHCA5jEkCIIgCKLWUO21gC1Famoqzpw5g+7du8PDwwPXrl3DvHnzkJubiwsXLsDJyQn5+fkICAhAeHg4Zs6cifPnz+O9997DsmXLMGHCBADl08D07NkT33zzDd58801s2bIFCxYsQFpaGlq3bq11XfoKmCAIgiAISyOqr4CtibOzM5KTk/Haa68hICAA0dHRaNOmDY4ePco1mbq5ueG3335DVlYWOnbsiGnTpuHzzz/ngj8A6Nq1KzZv3ozvv/8ebdu2xbZt27Bjxw6dwR9BEARBEERNxGZaAK1FXl4eGjRogLt375rcVU3wQ6FQ4Ouvv8bs2bMtOg6itkF+thzka8tBvrYM5GfLYe74gwLAKvjnn3/g6+uL7OxsNG7c2Nrm1Giou90ykJ8tB/nacpCvLQP52XKYO/6wmS7gmsDy5ctJj4UQW9rEpkcohLRHbD6qqb6uqXqERGxpE5seoaD6w8owwiDZ2dkMAGvZsiULDAxk8fHxJusKDAwUxKaaqic/P58BYPn5+dXWJba0iUmPGP0spC4x6RGjr2uqHvK1ZfSI0c9C6hKDnvj4eBYYGMhatmzJALDs7GxBbKqMUdPA1EZYRQ/5rl270KhRIwAweXmWsrIyQZZ2qal6VDrEZFNN1CNGPwupS0x6xOjrmqqHfG0ZPWL0s5C6xKBn1KhRGDVqFG7duoWgoCCzLXZBYwCr4O+//0aLFi2sbQZBEARBELWQa9euaSyIIRQUAFaBUqnE7du38cILL1R72TuCIAiCIAg+MMbw+PFj+Pj4QCoV/pMNCgAJgiAIgiBqGfQVMEEQBEEQRC2DAkCCIAiCIIhaBgWABEEQBEEQtQwKAAmCIAiCIGoZtTIAPHbsGN5++234+PhAIpFgx44dGscZY/j888/h7e2NOnXqICwsDFeuXNGQefDgAUaMGAFXV1e4u7sjOjoahYWFFkyFbVCVr5OTkxEeHo569epBIpEgIyNDS0dxcTEmT56MevXqoW7duhg0aBByc3MtkwAbwpCvS0pKMHPmTAQHB8PFxQU+Pj4YPXo0bt++raGD7uuqqeqejouLQ6tWreDi4gIPDw+EhYUhNTVVQ4b8zI+qfK3OxIkTIZFI8N1332nsJ1/zoypfjxkzBhKJRGPr27evhgz5umr43NMXL15EREQE3Nzc4OLigpCQENy8eZM7LtQzsVYGgE+ePEHbtm31LtWyaNEi/Oc//8GqVauQmpoKFxcX9OnTB8XFxZzMiBEjcOHCBRw4cAC7d+/GsWPHMGHCBEslwWaoytdPnjxB9+7dsXDhQr06pk6dil9++QVbt27F0aNHcfv2bURGRprLZJvFkK+fPn2KtLQ0zJs3D2lpaUhOTkZmZiYiIiI05Oi+rpqq7ml/f3/Ex8fjzz//xPHjx9G0aVOEh4cjLy+PkyE/86MqX6vYvn07Tp06BR8fH61j5Gt+8PF13759cefOHW778ccfNY6Tr6umKj9fu3YN3bt3R6tWrXDkyBGcO3cO8+bNg0wm42QEeyaaZX0RGwIA2759O/dbqVQyuVzOFi9ezO179OgRc3JyYj/++CNjjLG//vqLAWBnzpzhZPbt28ckEgm7deuWxWy3NSr7Wp2srCwGgKWnp2vsf/ToEXNwcGBbt27l9l28eJEBYCdPnjSjtbaNIV+rOH36NAPAbty4wRij+9oU+PhZtXTWwYMHGWPkZ1PR5+t//vmHNWrUiJ0/f541adKELVu2jDtGvjYNXb6Oiopi/fv313sO+dp4dPn5nXfeYSNHjtR7jpDPxFrZAmiIrKws5OTkICwsjNvn5uaG0NBQnDx5EgBw8uRJuLu74+WXX+ZkwsLCIJVKtbp6iOpx9uxZlJSUaORHq1at4Ofnx+UHYRr5+fmQSCRwd3cHQPe1OXj27Bm+//57uLm5oW3btgDIz0KiVCoxatQoxMTE4KWXXtI6Tr4WliNHjqBBgwYICAjApEmTcP/+fe4Y+br6KJVK7NmzB/7+/ujTpw8aNGiA0NBQjW5iIZ+JFABWIicnBwDQsGFDjf0NGzbkjuXk5KBBgwYax+3t7eHp6cnJEMKQk5MDR0dHLkhRoZ4fhPEUFxdj5syZGDZsGFxdXQHQfS0ku3fvRt26dSGTybBs2TIcOHAA9evXB0B+FpKFCxfC3t4eH3/8sc7j5Gvh6Nu3LzZs2IBDhw5h4cKFOHr0KPr164eysjIA5GshuHv3LgoLC/HNN9+gb9+++O233zBw4EBERkbi6NGjAIR9JtoLZThBELZBSUkJhg4dCsYYVq5caW1zaiS9e/dGRkYG7t27h9WrV2Po0KFITU3VekASpnP27Fn8+9//RlpaGi3TaQHeffdd7v/g4GC0adMGLVq0wJEjR/Daa69Z0bKag1KpBAD0798fU6dOBQC0a9cOJ06cwKpVq9CzZ09Br0ctgJWQy+UAoPVFTW5uLndMLpfj7t27GsdLS0vx4MEDToYQBrlcjmfPnuHRo0ca+9Xzg+CPKvi7ceMGDhw4wLX+AXRfC4mLiwtatmyJzp07Y+3atbC3t8fatWsBkJ+F4vfff8fdu3fh5+cHe3t72Nvb48aNG5g2bRqaNm0KgHxtTpo3b4769evj6tWrAMjXQlC/fn3Y29sjKChIY39gYCD3FbCQz0QKACvRrFkzyOVyHDp0iNtXUFCA1NRUdOnSBQDQpUsXPHr0CGfPnuVkDh8+DKVSidDQUIvbXJPp2LEjHBwcNPIjMzMTN2/e5PKD4Icq+Lty5QoOHjyIevXqaRyn+9p8KJVKKBQKAORnoRg1ahTOnTuHjIwMbvPx8UFMTAx+/fVXAORrc/LPP//g/v378Pb2BkC+FgJHR0eEhIQgMzNTY//ly5fRpEkTAMI+E2tlF3BhYSH31gKUf/iRkZEBT09P+Pn54ZNPPsGXX36JF198Ec2aNcO8efPg4+ODAQMGACiPxvv27Yvx48dj1apVKCkpwYcffoh3331X5zQEtZmqfP3gwQPcvHmTm49OdePL5XLI5XK4ubkhOjoan376KTw9PeHq6oqPPvoIXbp0QefOna2SJrFiyNfe3t4YPHgw0tLSsHv3bpSVlXHjRTw9PeHo6Ej3NU8M+blevXr46quvEBERAW9vb9y7dw/Lly/HrVu3MGTIEABUfxhDVfVH5ZcYBwcHyOVyBAQEACBfG4MhX3t6emL+/PkYNGgQ5HI5rl27hhkzZqBly5bo06cPAPI1X6q6p2NiYvDOO++gR48e6N27N/bv349ffvkFR44cAQBhn4lGfTNcQ0hJSWEAtLaoqCjGWPlUMPPmzWMNGzZkTk5O7LXXXmOZmZkaOu7fv8+GDRvG6taty1xdXdnYsWPZ48ePrZAacVOVrxMSEnQej42N5XQUFRWxDz74gHl4eDBnZ2c2cOBAdufOHeskSMQY8rVqmh1dW0pKCqeD7uuqMeTnoqIiNnDgQObj48McHR2Zt7c3i4iIYKdPn9bQQX7mR1X1R2UqTwPDGPmaL4Z8/fTpUxYeHs68vLyYg4MDa9KkCRs/fjzLycnR0EG+rho+9/TatWtZy5YtmUwmY23btmU7duzQ0CHUM1HCGGPGhYwEQRAEQRCELUNjAAmCIAiCIGoZFAASBEEQBEHUMigAJAiCIAiCqGVQAEgQBEEQBFHLoACQIAiCIAiilkEBIEEQBEEQRC2DAkCCIAiCIIhaBgWABEEQBEEQtQwKAAmCIAiCIGoZFAASBFFjOHLkCCQSCdavXy+47vXr10MikXBrctZUzOlDc2GLNhPmJS4uDhKJhNuOHz9ubZMMcurUKQ174+LizH5NCgAJghAlqoe6arOzs4OHhwdat26NqKgo7N+/HzVlJcuMjAzExcXh+vXr1jaF4AHll+2wbNkybNy4EQEBAVrHiouLsWLFCrz66qvw8vKCg4MD3N3dERISgpkzZ+LSpUsmXXPIkCGQSCTIyMjQK8MYQ7NmzeDu7o6ioiK0bNkSGzduxLJly0y6pinYW+xKBEEQJjBs2DC88cYbYIzh8ePHyMzMxI4dO7BhwwaEhYVh69atcHd3BwD06NEDRUVFcHBwsK7RRpKRkYH58+ejV69eaNq0qVVtsUUfWtpmMeUXYZgBAwbozKO///4bb731Fi5evIiePXti6tSp8Pb2RmFhITIyMrBu3TosWbIEN2/eRKNGjYy6ZnR0NLZt24aEhAT8+9//1imTkpKC69ev4/3330edOnVQp04djBw5EtevX8fUqVNNSarRUABIEISo6dChA0aOHKmxb+nSpZgxYwaWLl2KYcOGYd++fQAAqVQKmUxmDTNrDLbkw7KyMigUCjg7O9uMzYT1KSoqwptvvolr164hOTkZAwcO1JIpLi7GsmXLIJFIjNYfHh4OX19fJCYmYvHixXB0dNSSSUhIAFAeLFoL6gImCMLmsLOzw7fffovu3btj//793PgeY8eCqcb1HTx4EHFxcWjSpAmcnJzQpk0bbNmyRec5SqUSS5YsQYsWLeDk5AR/f3/88MMPOmXv3buHyZMnw9fXF46OjvD19cXkyZNx//59TiYuLg5jx44FAPTu3Zvr8h4zZoxRetTTc/jwYd42VkaXD/nq3bdvHyQSCf7zn//o1N2lSxd4eXmhpKQEAPD48WN89tlnCA0NRf369eHk5ISWLVti1qxZePr0qc60HTx4EP/617/QokULyGQy/PzzzzptNkU3H7/xyS9dqK5x6NAhfPHFF2jSpAnq1KmD0NBQnDp1CgBw9OhRdO/eHS4uLvD29sa//vUvDR3GpKm4uBhxcXEICAiAs7Mz3N3dERwcjJiYGJPkzJUma7BmzRpcunQJMTExOoM/AJDJZJg9ezZ8fHy4fQqFAgsWLMBLL70EmUwGd3d3vP3220hPT9c4VyqVYsyYMbh//z527dqlpbugoABJSUlo3bo1QkJChE2cEVALIEEQNkt0dDSOHz+OPXv2oHv37ibrmTlzJp48eYIPPvgAQPnb+bBhw1BcXKz1YJ8zZw6Kiorw/vvvw8nJCStXrsSYMWPQsmVLdOvWjZPLz89H165dcfXqVbz33nvo0KED0tPTsXLlShw+fBinT5/GCy+8gMjISNy5cwfff/895syZg8DAQABAixYtjNJjio3GUpXe8PBwyOVybNiwAR9//LHGuVeuXMGpU6fw8ccfc121t27dwpo1azBo0CAMHz4c9vb2OHr0KBYtWoT09HT8+uuvWjZMnz4dJSUlGD9+PFxdXREQEACFQqElZ4puPn6rKr+qYtasWSgrK8OUKVPw7NkzfPvttwgPD8eGDRsQHR2NCRMmYMSIEfj555/x+eefo1mzZlwLuDFpmjx5MtatW4fRo0fj008/RWlpKa5cuYLDhw9r2MNXzlxpsgbbtm0DAIwbN473OSUlJejbty9OnDiBUaNG4cMPP0R+fj5Wr16Nbt264dixY3j55Zc5+bFjx+LLL79EQkICBg8erKFry5YtKCoqsmrrHwCAEQRBiJCUlBQGgC1evFivzNmzZxkAFhkZqXFOQkICr2skJCQwAMzPz489evSI2//o0SPm5+fHPDw82NOnTzVk27VrxxQKBSf7zz//MEdHR/buu+9q6J4zZw4DwJYvX66xPz4+ngFgn332mZYdKSkpWjaaooevjbrQ5UNj9E6fPp0BYBcuXNDQ+9lnnzEA7OzZs9w+hULBnj17pmWDSjY1NVXLBn9/f/bkyZMqbTZFN1+/GcovfajOad++vcY1du7cyQAwe3t7dubMGQ375XI569y5s0lp8vDwYP369avSLr5y5kqTuYiNjWUAWFZWltYxT09P5urqqrW/tLSU5eXlaWyq8r906VIGgO3fv1/jnPz8fObr68t69uyppe/VV19ldnZ27Pbt2xr7O3fuzBwdHVleXp7WOVlZWQwAi42N5Z9YE6EuYIIgbBZXV1cA5V0q1WHSpElwc3Pjfru5uWHixIl4+PCh1rQvH3zwgcaYnkaNGsHf3x9XrlzRkNu+fTu8vLwwYcIEjf3vv/8+vLy8sH37dl62maKHr43GwkdvVFQUAGDDhg3cPsYYNm3ahNatW6NDhw7cfkdHR641sLS0FA8fPsS9e/cQFhYGAEhNTdWyYdKkSXB2dq7SVlN0m8tv6kyaNEnjGq+88goAIDQ0VKMFydHREZ06ddK4tjFpcnNzw4ULF3D+/HmD9vCVM1earEFBQQFXd6hz8eJFeHl5aWzLly8HAGzatAmtWrVCx44dce/ePW579uwZXn/9dRw/fhxFRUUa+qKjo1FWVqZRFi5duoRTp04hIiIC9evXN29Cq4ACQIIgbBZV4KerMgfKPxLIycnR2PLz87XkVN146gQFBQEo/1pQnebNm2vJ1qtXT2s8XlZWFgICAmBvrznSxt7eHv7+/lp69WGKHr42GgsfvaogLzExEUqlEgBw7NgxXL9+HaNHj9Y6f8WKFWjTpg2cnJzg6ekJLy8v9OrVCwDw8OFDLXl/f3/e9hqr21x+M3QNDw8PAECzZs20ZD08PLSuzTdN3333HR4+fIjg4GC0aNEC48aNw86dO7k8MVbOnGmyNK6urjpfGps1a4YDBw7gwIEDWLJkicaxixcv4tKlS1oBopeXF9atW4eysjLcu3dP45zIyEi4u7tzH3wAwLp16wAA7733nhlSZhw0BpAgCJvl3LlzAKBzji8AyM7O1noIRUVFVWvCYDs7O537mYjmJDSXjXz1jh49Gp988gkOHz6MsLAwbNiwAXZ2djq/5p42bRrCw8Px8ccfw8fHB46Ojrh16xbGjBmjMwjh0/pnqm5L5K2+a+jbr44xaerfvz+uX7+OvXv34ujRozh48CDWrl2LV155BQcPHuRa7PjKmStN1qB169Y4duwYsrKyNOoHFxcXrjW18gsXYwzBwcFYunSpXr1eXl4av2UyGYYPH44VK1bgxIkTCA0NxcaNG9G4cWP06dNHwBSZBgWABEHYLGvXrgUAvPnmmzqPy+VyHDhwQGOf+ld9Ki5evIj+/ftr7Pvrr78A6G4V4kPz5s2RmZmJ0tJSjYdJaWkpLl++rKHX0FQTxugRC8OHD0dMTAw2bNiAbt26Ydu2bXj99dfh7e2tIbdx40Y0bdoU+/btg1T6vENq//791bbBnLpNmRpECIxNk6enJ0aOHImRI0eCMYZZs2Zh0aJF2LlzJ4YMGWK0nNC0a9cOcrkcu3fv5u7tvXv3Yty4ccjOzkbHjh0NHjc1wBw8eDCOHTuGNWvW4KuvvuJ1zosvvoi8vDy8+uqrGr6viujoaKxYsQIJCQl48OABcnJyMHfuXKN0mAvrW0AQBGEkZWVlmD59Oo4fP4433nhD75etMpkMYWFhGpuqa1edlStXanQN5+fnY9WqVXB3d0fPnj1NsnHAgAHIy8vDmjVrNPavXr0aeXl5GtNP1K1bFwDw4MGDaukRC15eXujXrx+Sk5ORmJiIgoICbmygOnZ2dpBIJBotbKWlpfjmm2+qbYM5dRvKL3PCN01lZWV49OiRxj6JRIL27dsDeG43XzlzcfLkSWRnZ2PHjh3cvtWrVyMqKgp2dnZVHjeVcePGoVWrVli8eLHesbi6WrVzcnL0tgDm5ubq3N+hQwe0a9cOP/30E5YvXw6JRCKK7l+AWgAJghA5aWlp2LRpEwBorARy48YNhIeHY/PmzdW+Rv369REaGsrN75aQkICbN29izZo1vLscKzNjxgxs3boVkydPRlpaGtq3b4/09HSsXbsWAQEBmDFjBicbEhICqVSKr776Cg8fPoSLiwuaNWuG0NBQo/SIiaioKOzatQvTpk2Dm5sbBgwYoCUzePBgzJ49G/369UNkZCQKCgqwefNmQVb0MKduQ/llTvim6fHjx/D29kZERATat2+PBg0aICsrCytXroSHhwfefvtto+TMRZ06dRAdHY2ffvoJgwcPRk5ODvbs2YNFixbxOl6d6+7ZswdvvfUWIiMj0atXL24Ko4KCAly6dAk//fQT7Ozs4OvrCwCYMmUKDhw4gJiYGBw+fBivvvoqXF1dcfPmTRw6dAgymQwpKSk6rxcdHY2PPvoI+/fvR69evUTTak8BIEEQoubHH3/Ejz/+CKlUirp166Jx48bo2bMnhg0bhr59+wpyjYULF+L333/H8uXLkZubC39/fyQmJmL48OEm63Rzc8P//d//ITY2Frt27UJCQgIaNmyIiRMnYv78+Rpz9/n5+WHdunVYuHAhJk2ahJKSEkRFRSE0NNQoPWLirbfegqenJx48eIBx48bpXKkjJiYGjDGsXbsWU6ZMgVwuxzvvvIOxY8fqbKk1BnPqNpRf5oRvmpydnfHJJ5/g0KFDOHjwIAoLC7lAT31yY75y5mTIkCGYN28eCgsLsX79enTu3Bkvvvgi7+Om0rx5c5w9exbr1q3Dtm3b8O233yI/Px8uLi5o2bIlxo0bh+joaG58sYODA/bs2YMVK1Zg48aNiI2NBVA+pKRTp046W7hVjBgxAjExMSguLhZN6x8AmgeQIIjaiynzuREEISydO3dmiYmJ7MUXX2Tr1683+rguVPMApqWlsby8PFZSUiK02YJSUlLC8vLyWFpaGs0DSBAEQRBEzWfIkCGYM2cOcnJydH50UtVxQ3To0AFeXl7c0nRi5Y8//oCXl5fGPJnmhrqACYIgCIKwGoMHD8a0adMwfvx4nWNuqzqui9GjR2ssDxkcHCyYveYgKChIY8YCS4wTpACQIAiCIAir4efnB19fX24OPmOP66J58+ai+diCD66urkalTwioC5ggiFrLmDFjwBjjVlIgCMI6SKVSg3PjVXWcMB7yJkEQBEEQRC1DwpiI1i8iCIIgCIIgzA61ABIEQRAEQdQyKAAkCIIgCIKoZVAASBAEQRAEUcugAJAgCIIgCKKWQQEgQRAEQRBELYMCQIIgCIIgiFoGBYAEQRAEQRC1DAoACYIgCIIgahkUABIEQRAEQdQyKAAkCIIgCIKoZVAASBAEQRAEUcv4fyL6GcKnffVyAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "data_x,_ = np.histogram(ak.to_numpy(all_data),\n",
    "                            bins=bin_edges ) # histogram the data\n",
    "data_x_errors = np.sqrt( data_x ) # statistical error on the data\n",
    "\n",
    "# data fit\n",
    "polynomial_mod = PolynomialModel( 4 ) # 4th order polynomial\n",
    "gaussian_mod = GaussianModel() # Gaussian\n",
    "\n",
    "# set initial guesses for the parameters of the polynomial model\n",
    "# c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "pars = polynomial_mod.guess(data_x, # data to use to guess parameter values\n",
    "                            x=bin_centres, c0=data_x.max(), c1=0,\n",
    "                            c2=0, c3=0, c4=0 )\n",
    "\n",
    "# set initial guesses for the parameters of the Gaussian model\n",
    "pars += gaussian_mod.guess(data_x, # data to use to guess parameter values\n",
    "                        x=bin_centres, amplitude=100,\n",
    "                        center=125, sigma=2 )\n",
    "\n",
    "model = polynomial_mod + gaussian_mod # combined model\n",
    "\n",
    "# fit the model to the data\n",
    "out = model.fit(data_x, # data to be fit\n",
    "                pars, # guesses for the parameters\n",
    "                x=bin_centres, weights=1/data_x_errors ) #ASK\n",
    "\n",
    "# background part of fit\n",
    "params_dict = out.params.valuesdict() # get the parameters from the fit to data\n",
    "c0 = params_dict['c0'] # c0 of c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "c1 = params_dict['c1'] # c1 of c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "c2 = params_dict['c2'] # c2 of c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "c3 = params_dict['c3'] # c3 of c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "c4 = params_dict['c4'] # c4 of c0 + c1*x + c2*x^2 + c3*x^3 + c4*x^4\n",
    "\n",
    "# get the background only part of the fit to data\n",
    "background = c0 + c1*bin_centres + c2*bin_centres**2 + c3*bin_centres**3 + c4*bin_centres**4\n",
    "\n",
    "# data fit - background fit = signal fit\n",
    "signal_x = data_x - background\n",
    "\n",
    "# *************\n",
    "# Main plot\n",
    "# *************\n",
    "plt.axes([0.1,0.3,0.85,0.65]) # left, bottom, width, height\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' means circles\n",
    "                label='Data', markersize=4 )\n",
    "\n",
    "# plot the signal + background fit\n",
    "main_axes.plot(bin_centres, # x\n",
    "            out.best_fit, # y\n",
    "            '-r', # single red line\n",
    "            label='Sig+Bkg Fit ($m_H=125$ GeV)' )\n",
    "\n",
    "# plot the background only fit\n",
    "main_axes.plot(bin_centres, # x\n",
    "            background, # y\n",
    "            '--r', # dashed red line\n",
    "            label='Bkg (4th order polynomial)' )\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",
    "                    labelbottom=False, # don't draw tick labels on bottom axis\n",
    "                    right=True ) # draw ticks on right axis\n",
    "\n",
    "# write y-axis label for main\n",
    "main_axes.set_ylabel('Events / '+str(step_size)+' GeV',\n",
    "                    horizontalalignment='right')\n",
    "\n",
    "# set the y-axis limit for the main axes\n",
    "main_axes.set_ylim( bottom=0, top=np.amax(data_x)*1.5 )\n",
    "\n",
    "# set minor ticks on the y-axis of the main axes\n",
    "main_axes.yaxis.set_minor_locator( AutoMinorLocator() )\n",
    "\n",
    "# avoid displaying y=0 on the main axes\n",
    "main_axes.yaxis.get_major_ticks()[0].set_visible(False)\n",
    "\n",
    "# Add text 'ATLAS Open Data' on plot\n",
    "plt.text(0.2, # x\n",
    "        0.92, # y\n",
    "        'ATLAS Open Data', # text\n",
    "        transform=main_axes.transAxes, # coordinate system used is that of main_axes\n",
    "        fontsize=13 )\n",
    "\n",
    "# Add text 'for education' on plot\n",
    "plt.text(0.2, # x\n",
    "        0.86, # y\n",
    "        'for education', # text\n",
    "        transform=main_axes.transAxes, # coordinate system used is that of main_axes\n",
    "        style='italic',\n",
    "        fontsize=8 )\n",
    "\n",
    "lumi = 36.1\n",
    "lumi_used = str(lumi*fraction) # luminosity to write on the plot\n",
    "plt.text(0.2, # x\n",
    "        0.8, # y\n",
    "        r'$\\sqrt{s}$=13 TeV,$\\int$L dt = '+lumi_used+' fb$^{-1}$', # text\n",
    "        transform=main_axes.transAxes ) # coordinate system used is that of main_axes\n",
    "\n",
    "# Add a label for the analysis carried out\n",
    "plt.text(0.2, # x\n",
    "        0.74, # y\n",
    "        r'$H \\rightarrow \\gamma\\gamma$', # text\n",
    "        transform=main_axes.transAxes ) # coordinate system used is that of main_axes\n",
    "\n",
    "# draw the legend\n",
    "main_axes.legend(frameon=False, # no box around the legend\n",
    "                loc='lower left' ) # legend location\n",
    "\n",
    "\n",
    "# *************\n",
    "# Data-Bkg plot\n",
    "# *************\n",
    "plt.axes([0.1,0.1,0.85,0.2]) # left, bottom, width, height\n",
    "sub_axes = plt.gca() # get the current axes\n",
    "\n",
    "# set the y axis to be symmetric about Data-Background=0\n",
    "sub_axes.yaxis.set_major_locator( MaxNLocator(nbins='auto',\n",
    "                                            symmetric=True) )\n",
    "\n",
    "# plot Data-Background\n",
    "sub_axes.errorbar(x=bin_centres, y=signal_x, yerr=data_x_errors,\n",
    "                fmt='ko',markersize=4 ) # 'k' means black and 'o' means circles\n",
    "\n",
    "# draw the fit to data\n",
    "sub_axes.plot(bin_centres, # x\n",
    "            out.best_fit-background, # y\n",
    "            '-r' ) # single red line\n",
    "\n",
    "# draw the background only fit\n",
    "sub_axes.plot(bin_centres, # x\n",
    "            background-background, # y\n",
    "            '--r' )  # dashed red line\n",
    "\n",
    "# set the x-axis limits on the sub axes\n",
    "sub_axes.set_xlim( left=xmin, right=xmax )\n",
    "\n",
    "# separation of x-axis minor ticks\n",
    "sub_axes.xaxis.set_minor_locator( AutoMinorLocator() )\n",
    "\n",
    "# x-axis label\n",
    "sub_axes.set_xlabel(r'Di-photon invariant mass $\\mathrm{m_{\\gamma\\gamma}}$ [GeV]',\n",
    "                    x=1, horizontalalignment='right',\n",
    "                    fontsize=13 )\n",
    "\n",
    "# set the tick parameters for the sub axes\n",
    "sub_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",
    "# separation of y-axis minor ticks\n",
    "sub_axes.yaxis.set_minor_locator( AutoMinorLocator() )\n",
    "\n",
    "# y-axis label on the sub axes\n",
    "sub_axes.set_ylabel( 'Events-Bkg' )\n",
    "\n",
    "\n",
    "# Generic features for both plots\n",
    "main_axes.yaxis.set_label_coords( -0.09, 1 ) # x,y coordinates of the y-axis label on the main axes\n",
    "sub_axes.yaxis.set_label_coords( -0.09, 0.5 ) # x,y coordinates of the y-axis label on the sub axes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Hjy3wlBoSh4J"
   },
   "source": [
    "And there it is: a clear peak in the invariant mass spectrum around 125 GeV, signaling the presence of the Higgs boson!\n",
    "\n",
    "This represents a very close result to the one obtained by ATLAS in [Higgs decay to photons measurement paper](https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.052005), which can be seen below:\n",
    "<CENTER><img src=\"https://github.com/atlas-outreach-data-tools/notebooks-collection-opendata/blob/master/13-TeV-examples/uproot_python/images/ImagesHiggs/HyyPaper.png?raw=1\" style=\"width:50%\"></CENTER>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ii29QaMPSh4K"
   },
   "source": [
    "While our main task may be done, there's still more to explore. Here are some additional tasks you can try with this notebook:\n",
    "* Check how many events are being thrown away by each cut in '[Applying a cut](#applying_cut)'\n",
    "* Add more cuts from the [Higgs boson discovery paper](https://www.sciencedirect.com/science/article/pii/S037026931200857X#se0090) or the [Higgs decay to photons measurement paper](https://journals.aps.org/prd/abstract/10.1103/PhysRevD.98.052005) in '[Changing a cut](#changing_cut)' and '[Applying a cut](#applying_cut)'\n",
    "* Find the reduced chi-squared for the fit in '[Plotting](#plotting)'\n",
    "* Find the mean of the fitted Gaussian in '[Plotting](#plotting)'\n",
    "* Find the width of the fitted Gaussian in '[Plotting](#plotting)'\n",
    "* Try different initial guesses for the parameters of the fit in '[Plotting](#plotting)'\n",
    "* Try different functions for the fit in '[Plotting](#plotting)'\n",
    "* Your idea!",
    "\n",
    "\n",
    "<div class=\"alert alert-block alert-info\">\n",
    "We welcome your feedback on this notebook or any of our other materials! Please <a href=\"https://forms.gle/zKBqS1opAHHemv9U7\">fill out this survey</a> to let us know how we're doing, and you can enter a raffle to win some <a href=\"https://atlas-secretariat.web.cern.ch/merchandise\">ATLAS merchandise</a>!\n",
    "</div>"
   ]
  }
 ],
 "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.11.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}