{ "cells": [ { "cell_type": "markdown", "metadata": { "id": "cJjTKA3n6oe4" }, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": { "id": "uW9nIM636oe5" }, "source": [ "# Find the Z boson yourself!\n", "This notebook uses ATLAS Open Data to guide you through the steps to find the Z boson. It takes less than an hour to complete, and by the end, you’ll be able to identify the Z boson in ATLAS detector data from the 2015 and 2016 runs.\n", "\n", "## The Z boson\n", "In the [Standard Model of particle physics](https://opendata.atlas.cern/docs/documentation/introduction/SM_and_beyond) the Z boson is one of the two bosons responsible for mediating the weak interaction. It has a very short lifetime, of about $10^{25}$ seconds, and decays via the weak force into leptons and quarks. Below you can see a [Feynman diagram](https://en.wikipedia.org/wiki/Feynman_diagram) of the decay of the Z boson into two leptons.\n" ] }, { "cell_type": "markdown", "metadata": { "id": "AzwFKITg6oe6" }, "source": [ "
" ] }, { "cell_type": "markdown", "metadata": { "id": "fw4O3JeN6oe6" }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "id": "Nej9yxUr6oe6" }, "source": [ "Due to its short lifetime, we don't measure the Z boson directly in the ATLAS detector. What we measure is its decays products (in this case the two leptons) and reconstruct its mass. Ready to learn how to do this? Let's go!.\n", "\n", "**Contents:**\n", "\n", "[Running a Jupyter notebook](#running)
\n", "[To setup everytime](#setupeverytime)
\n", "[Where's my data](#fraction)
\n", "[Calculate that invariant mass!](#good_leptons)
\n", "[Can we process the data yet?!](#load_data)
\n", "[Plot Data](#plot_data)
\n", "[Your tasks](#tasks)
\n", "[Going further](#going_further)
" ] }, { "cell_type": "markdown", "metadata": { "id": "63dlPOvU6oe6" }, "source": [ "## Running a Jupyter notebook\n", "\n", "To run the whole Jupyter notebook, in the top menu click Cell $\\rightarrow$ Run All.\n", "\n", "To propagate a change you've made to a piece of code, click Cell $\\rightarrow$ Run All Below.\n", "\n", "You can also run a single code cell, by using the keyboard shortcut Shift+Enter." ] }, { "cell_type": "markdown", "metadata": { "id": "CGrKbBKW6oe6" }, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": { "id": "VSEhvige6oe6" }, "source": [ "## First time setup on your computer (no need on mybinder)\n", "This first cell only needs to be run the first time you open this notebook on your computer.\n", "\n", "If you close Jupyter and re-open it on the same computer, you won't need to run this first cell again.\n", "\n", "If you open the notebook on mybinder, you don't need to run this cell." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "bl53Y7166oe7", "outputId": "2cf12c15-91c3-464d-d323-66f34463e900" }, "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": "1_-X1tru6oe7" }, "source": [ "## To setup everytime\n", "We're going to be using a number of tools to help us:\n", "* urllib: let us download files\n", "* uproot: lets us read .root files typically used in particle physics into data formats used in Python\n", "* pandas: lets us store data as dataframes, a format widely used in Python\n", "* numpy: provides numerical calculations such as histogramming\n", "* matplotlib: common tool for making plots, figures, images, visualisations\n", "* awkward: a python package for dealing conveniently with lists that are different lengths\n", "* vector: a python package for making and manipulating vectors\n", "\n", "In case you see an error when importing one of these packages, try restarting the kernel (the \"session\" in Colab) to make sure the packages that were installed in the last cell have been properly loaded." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "scrolled": true, "id": "f-FtOWXM6oe8" }, "outputs": [], "source": [ "import urllib.request # for downloading files\n", "import pandas as pd # to store data as dataframes\n", "import numpy as np # for numerical calculations such as histogramming\n", "import uproot # to read .root files as dataframes\n", "import matplotlib.pyplot as plt # for plotting\n", "from matplotlib.ticker import MaxNLocator,AutoMinorLocator # for minor ticks\n", "import awkward as ak # for handling complex and nested data structures efficiently\n", "import vector # For convenient 4-vector manipulation" ] }, { "cell_type": "markdown", "metadata": { "id": "dofdCbKC6oe8" }, "source": [ "## Where's my data?" ] }, { "cell_type": "markdown", "source": [ "We're going to use a package called `atlasopenmagic` to find the data we want to use." ], "metadata": { "id": "xKtXjFoL7lP3" } }, { "cell_type": "code", "execution_count": 3, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "NiQJJMpe6oe8", "outputId": "c00e7192-767c-47d0-bbc0-1b386e8426eb" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Active release set to: 2025e-13tev-beta. Metadata cache cleared.\n", "Fetching and caching all metadata for release: 2025e-13tev-beta...\n", "Successfully cached 374 datasets.\n" ] } ], "source": [ "# Import the atlasopenmagic package so that we can use it\n", "import atlasopenmagic as atom\n", "# Make sure that we are using the newest release of Education and Outreach open data\n", "atom.set_release('2025e-13tev-beta')\n", "# We will use a skim (pre-selection) of the data so that we start from events\n", "# with exactly two muons\n", "skim = '2muons'\n", "# And now we can get the files we want; note that we use \"cache=True\" to copy\n", "# the files locally rather than streaming them\n", "files_list = atom.get_urls('data', skim, protocol='https', cache=True)" ] }, { "cell_type": "markdown", "metadata": { "id": "6UzGr5SY6oe8" }, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": { "id": "4xIG8vMg6oe8" }, "source": [ "## Calculate that invariant mass!" ] }, { "cell_type": "markdown", "metadata": { "id": "reNmdNRi6oe8" }, "source": [ "In particle physics, the momentum of a particle is often written as a 4-vector: `(E, px, py, pz)`, where `E` is the energy of the particle, `px` is the momentum of the particle in the `x` direction, and so on. When we add them together, we add each component one-by-one. If a particle decays into two particles (for example, a Z boson decays into two leptons), then the sum of the 4-vectors of the two leptons is equal to the 4-vector of the original Z boson.\n", "\n", "One other useful trick is that the mass of an object can easily be calculated from its 4-vector:\n", "\n", "$$M^2=E^2-p^2,$$\n", "\n", "where $p$ is the momentum $(px, py, pz)$. We are using \"natural units\" here, where the speed of light is 1 ($c=1$); if we weren't, that formula would be written:\n", "\n", "$$(M \\times c^2)^2=E^2-p^2,$$\n", "\n", "That might start to look familiar. If a particle is at rest (it has no momentum, so $p=0$), that means:\n", "\n", "$$(M \\times c^2)^2 = E^2,$$\n", "\n", "or simply\n", "\n", "$$M \\times c^2 = E$$\n", "\n", "which is Einstein's famous formula!\n", "\n", "The `vector` module has a handy way to put 4-vectors together. We can give it the components of the 4-momentum that we have, and it will conveniently put it into the `(E, px, py, pz)` form for us. In particle physics, and in our Open Data files, we store:\n", "\n", "* `pt`, also written $p_\\text{T}$, which is the part of the momentum transverse (perpendicular) to the beam's original direction ($p_\\text{T}=\\sqrt{px^2+py^2}$).\n", "* `eta`, [pseudorapidity](https://en.wikipedia.org/wiki/Pseudorapidity), which is a convenient way to describe the angle of the particle's momentum in the detector.\n", "* `phi`, the angle of the particle's momentum in the $x, y$ plane.\n", "* `e`, the energy of the particle\n", "\n", "We can stick these four components together and calculate all the others from them.\n", "\n", "After that, all we have to do is add together the 4-momenta of the first two leptons, and then ask for the mass (also called the \"invariant mass\" because it doesn't matter how the particle is moving or if it is at rest, the mass is always the same)." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "id": "rPretFv76oe9" }, "outputs": [], "source": [ "# calculate dilepton invariant mass\n", "def calc_mll(lep_pt,lep_eta,lep_phi,lep_e):\n", " p4 = vector.zip({\"pt\": lep_pt, \"eta\": lep_eta, \"phi\": lep_phi, \"e\": lep_e})\n", " invariant_mass = (p4[:, 0] + p4[:, 1]).M # .M calculates the invariant mass\n", " return invariant_mass" ] }, { "cell_type": "markdown", "metadata": { "id": "2O7fWFCS6oe9" }, "source": [ "## Can we process the data yet?!" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "colab": { "base_uri": "https://localhost:8080/" }, "id": "yfEnCbtP6oe9", "outputId": "e531cfa2-c19c-4431-bbb4-c664ccadb889" }, "outputs": [ { "output_type": "stream", "name": "stdout", "text": [ "Working on file simplecache::https://opendata.cern.ch/eos/opendata/atlas/rucio/opendata/ODEO_FEB2025_v0_2muons_data15_periodD.2muons.root (0/16)\n" ] } ], "source": [ "mass_list = [] # This list will hold all our output masses\n", "\n", "# Loop through all the files that we have\n", "for afile in files_list:\n", " print(f'Working on file {afile} ({files_list.index(afile)}/{len(files_list)})')\n", "\n", " tree = uproot.open(afile+\":analysis\") # open the file, and the tree called mini\n", " numevents = tree.num_entries # number of events\n", "\n", " for data in tree.iterate(['lep_pt','lep_eta','lep_phi','lep_e'],\n", " entry_stop=numevents*0.5, # stop after fraction of events we want to process\n", " library=\"ak\"):\n", " # calculation of 2-lepton invariant mass\n", " data['mll'] = calc_mll(data.lep_pt, data.lep_eta, data.lep_phi, data.lep_e)\n", " # Keep the mll as a flat array\n", " mass_list.append( data['mll'] )\n", " # For simplicity and speed, we'll just work with the first file to start\n", " # Stop here!\n", " break" ] }, { "cell_type": "markdown", "metadata": { "id": "3WulUUJM6oe9" }, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": { "id": "iriGq3y46oe9" }, "source": [ "## Function to plot the data" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "scrolled": true, "colab": { "base_uri": "https://localhost:8080/", "height": 449 }, "id": "9M8n6gRp6oe9", "outputId": "3191fe2c-fa54-4e7b-b9b1-11b175182bee" }, "outputs": [ { "output_type": "display_data", "data": { "text/plain": [ "
" ], "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAGwCAYAAABFFQqPAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAARkBJREFUeJzt3X10VOW59/HfZjSAAYK8B2ZIBItlVEA0UqzBUBAQRTQnj4oKEaytmtYgrRXLarGtFVtbTaw5tfW0WKioNYzoqmgVTCRYagGharUqNtoQE+RFCAQNMrmfPzwzJ8NMSCbzPvv7WStLZu979lw3e8i+vF8tY4wRAACADXVLdAAAAACJQiIEAABsi0QIAADYFokQAACwLRIhAABgWyRCAADAtkiEAACAbZEIAQAA2zoh0QEkUmtrqz766CP17t1blmUlOhwAANAJxhgdPHhQQ4cOVbdukbXp2DoR+uijj+RyuRIdBgAA6IK6ujo5nc6IrmHrRKh3796SvviL7NOnjyQpLy9Pmzdv7tT7m5qa5HK5At7fkXCun2zlw61vMsUebvlku7d2qmusy/M9Tlw8dqprrMvb/XvsO+Z7jkfC1omQrzusT58+/r9ch8PR6X9EPm3f35Fwr59s5aXO1zfZYk/le2unusajvMT3OBHx2Kmu8Sgv8T2OxrAWBksfo6SkJKmun2zlY3ntZCsfrljGY6e6xqN8LK+dbOXDxfc4dcrH8trJVj6WrEh2n9+xY4d27NihGTNm+I+9+uqruuuuu7Rv3z4VFxfrG9/4RlQCjYWmpiZlZWXpwIEDYWem0Xh/qrFTfalr+rJTfalr+rJTfUPVNZr1j6hr7Pbbb9e+ffv8idCePXt00UUX6dChQ+rZs6duuukmDRo0SJdddllEQSar7t27a+nSperevXuiQ4kLO9WXuqYvO9WXuqYvO9U31nWNqEUoJydH3/jGN7RkyRJJ0q9+9Svdeuut2r59u0aNGqWCggJ1795dVVVVUQs4muyUUQMAkC6i+fyOaIzQ7t27NXToUP/r559/Xl/96ld1xhlnKCMjQ1dddZXeeuutiAIEAACIlYgSoczMTO3fv1+S5PV6tXHjRk2aNMl/vmfPnmpqaoooQAAAgFiJKBE6/fTTtWLFCu3du1cPP/ywDh06pAsvvNB//sMPP9TAgQMjDhIAACAWIkqEbrvtNr3xxhsaNGiQSkpKdNZZZyk/P99//oUXXtD48eMjDjLW8vLy5Ha7VVFRkehQAABAOyoqKuR2u5WXlxe1a0Y0WFqSXn75ZT3zzDPKysrSt771LfXr10+StHfvXt1www2aO3euLr/88qgEG20MlgYAIPVE8/kdcSKUykiEAADJqLm5Wb169ZIkHTp0SJmZmQmOKLkkzawxh8OhVatWtXv+iSeekMPhiOQjACCmvF6vqqur9dhjj6m6ulperzfRIQGIo4gSoY4ak2zc2AQgBXg8HuXm5mry5Mm6+uqrNXnyZOXm5srj8SQ6NNhc24R8w4YNMU/Qr7vuOlmWJcuydOKJJ2rw4MG68MIL9fvf/16tra2dvs4jjzyivn37xi7QGIjpXmP/+c9/orIzLABEm8fjUVFRkXbu3BlwvL6+XkVFRe0mQ7QgIdY8Ho/cbrf/9cyZM+OSoM+YMUMNDQ364IMP9Nxzz2ny5MkqLS3VJZdcoqNHj8b0sxPKhGnNmjVm/vz5Zv78+cayLHPBBRf4X7f9mT17tsnMzDTTp08P9yPi5sCBA0aSOXDgQKJDARBHR48eNU6n00gK+WNZlnG5XObo0aMB71u9enXQ+5xOp1m9enWCaoJ0s3r1amNZVsjvpGVZMfuuFRcXm9mzZwcdX79+vZFkHn74YWOMMb/85S/NGWecYU466STjdDrNTTfdZA4ePGiMMaaqqioo7qVLlxpjjFmxYoU5++yzTa9evczgwYPNnDlzzK5du7ocbzSf32EnQnfeeaf/hnTr1s3/52N/evfubS688ELz3nvvRRxkrJAIAfYU6hd2qJ+qqir/exL1gIJ9dDVBj4b2EiFjjBk7dqy56KKLjDHG3H///eall14ytbW1Zv369ea0004zN910kzHGmJaWFlNWVmb69OljGhoaTENDgz9J+t3vfmfWrl1r3n//fbNp0yYzceJE/zW7IqGJUFuWZZlHH3004iAShUQIsKdVq1Z1KhFatWqVMSaxDyjYR1cS9Gg5XiJ05ZVXmtGjR4c89+STT5r+/fv7Xy9fvtxkZWV1+HmbN282kvyJUrii+fyOaIxQbW1t2u4sDyB9ZWdnh1WupqYmaCxRW8YY1dXVqaamJirxwZ4aGhqiWi5ajDGyLEuStG7dOk2ZMkXDhg1T7969NXfuXO3du1eHDx8+7jW2bt2qWbNmafjw4erdu7cuuOACSV+MJU60iBKhnJwcnXTSSdGKBQCioqMBzfn5+XI6nf5f7seyLEsul8u/Un6yPqCQXsJN0OPl7bff1imnnKIPPvhAl1xyicaMGaPVq1dr69at/h0Zjhw50u77m5ubNX36dPXp00ePPvqoNm/erKeeeqrD98XLCZFeYNOmTXrwwQf13nvvae/evUFT5i3L0vvvvx/pxwBAp3g8HpWWlga04DidTpWXl6uwsFDSF2uglZeXq6ioSJZlBfze8iVHZWVl/nXQkvUBhfTiS9Dr6+tDLj9jWZacTmfAVlax9tJLL+mNN97Qrbfeqq1bt6q1tVW//OUv1a3bF+0of/rTnwLKZ2RkBP2Px7/+9S/t3btX99xzj1wulyRpy5Yt8alAJ0TUIrRixQqdf/75Wr16tT777DMNHz5cOTk5AT/Dhw+PVqwAcFzhTIkvLCxUZWWlhg0bFlDW6XSqsrLSnzRJ4bcgAV3hS9AlBX3XQiXo0dbS0qLGxkbV19frtdde0913363Zs2frkksu0bx583Tqqafq888/169+9Sv9+9//1sqVK/XQQw8FXCM3N1eHDh3S+vXrtWfPHh0+fFjDhw9XRkaG/33PPPOMfvKTn8SkDl0SyQCjUaNGmdGjR5v6+vqIByslAoOlgfTR1QHNR48eNVVVVWbVqlWmqqqq3QHPvlljx84cY9YYom316tVm2LBhAd8zl8sV0+9YcXGx/7NOOOEEM3DgQDN16lTz+9//3ni9Xn+5++67z2RnZ5uePXua6dOnmxUrVhhJ5pNPPvGXufHGG03//v0Dps+vWrXK5Obmmu7du5uJEyeaZ555xkgy27Zt61K80Xx+R7TXWI8ePXTvvffq29/+dpcTsURirzEgfVRXV2vy5MkdlquqqlJBQUGXPiNUt5vL5VJZWVlACxIQKd/zSZLWrl2radOmsWVVG9F8fkc0RsjpdKqlpSWiAAAgGuIxoLmwsFCzZ89WTU2NGhoalJ2drfz8fB5QiLq236lJkybxHYuhiBKhG2+8UY8++qhuvfVWbhKAhIrXgGaHw9HlFiWgszIzM9mvM04iSoTOPvtsrV69Wueee65KSkp0yimnhEyIJk2aFMnHAECHknHGDYDkF9EYId/0Of/Fjhnlbv53EaZk3ZSQMUJAevHNGpMUckr8sbPBAKSmpBkjtHz58og+HACiyTclPtQ6QgxoBhBKRC1CqY4WISA9eb1eBjQDaSxpWoQAIBkl04BmkjIguUW0srQk1dXVacGCBXI6ncrIyNBLL70kSdq9e7cWLFigzZs3RxwkAKQij8ej3NxcTZ48WVdffbUmT56s3NzcgBWuASRWxLvPn3POOVq9erVOP/30gEHRAwcO1JYtW/Q///M/EQcJAKkmnO0+ACRORInQkiVL1K1bN7355pt69NFHg6aszpw5Uxs3bowoQABINV6vV6WlpSGn8fuOLVy4MGln1AJ2ElEitG7dOt18881yuVwhNyPMyckJ+r8hAEh3NTU1x/3dZ4xRXV2dampq4hgVgFAiSoSampqOu0rrkSNHdPTo0Ug+AgBSTjy2+wAQHRElQi6XS//85z/bPf+3v/1Np556aiQfAQApJ17bfQCIXESJUGFhoX7/+9/rzTff9B/zdZGtXr1aTz75pK644orIIoyDvLw8ud1uVVRUJDoUAGnAt91HqCED0he/J10uF9t9AGGqqKiQ2+1WXl5e1K4Z0YKKTU1Nmjhxoj744ANNmjRJL7zwgqZOnaqmpib9/e9/17hx4/TKK6+oR48eUQs4mlhQEUCssN0HEDvRfH5H1CLUp08fbdq0SV//+te1ZcsWGWP04osv6p133tHNN9+sqqqqpE2CAKQOr9er6upqPfbYY6qurk6J2Va+7T6GDRsWcNzpdJIEAUkkqlts7N69W8YYDRw4sN0m4WRCixCQ/DweT8i9w8rLy1MimWBlaSD6ovn8jigRev311zVmzJiIAkgkEiEgufm6l479NUX3EmBvSdM1Nm7cOI0fP17l5eXavXt3RIEAQFssSgggHiJKhG6//Xbt3btXt956q5xOpy699FKtXr1aR44ciVZ8AGyKRQkBxENEidCyZcv0wQcf6MUXX9RVV12l6upqXXHFFcrOzlZJSYleffXVaMUJwGZYlBBAPES8+7xlWZoyZYr+8Ic/qLGxUY888ojGjx+v3/zmNzrvvPM0evToaMQJwGZYlBBAPER11lhbq1at0s0336yDBw8mbR8+g6WB5OX1epWbm6v6+vqQ44Qsy5LT6VRtbS2zsACbSZrB0sfasWOHfvjDH2rEiBGaO3euDh8+rEsuuSSaHwHAJhwOh8rLyyUpaDkO3+uysjKSIAARiTgR2r9/vx566CGdd955Ou2003TXXXcpKytLv/zlL1VfX6+nn346GnECsCEWJQQQaxF1jRUVFenZZ59VS0uLBg8erGuuuUbz5s1LmbWF6BoDUgOLEgJoK2kWVOzZs6cuvfRSFRcXa/r06Sn3i4lECACA1BPN5/cJkby5sbFRWVlZ7Z4/fPiwGhsbNWLEiEg+BgAAICbCHiOUkZGhxx9/XJKUlZWlgwcP6tJLL9Ubb7wRVPapp57Sl770pcijBAAAiIGwE6GjR4+qtbXV//rIkSP685//zBYbAAAg5UR1+jwAAEAqiWiMEAAgupghB8QXiRAAJAmPx6PS0tKAzWadTqfKy8tZMwmIEbrGACAJeDweFRUVBSRBklRfX6+ioiJ5PJ4ERQakty61CK1du1aNjY2Svpgib1mWnnzySW3fvj2g3NatWyMOEADSndfrVWlpacg91YwxsixLCxcu1OzZs+kmA6Is7AUVu3ULrxHJsiw2XQWA46iurtbkyZM7LFdVVaWCgoLYBwQkuYQuqFhVVRXRBwIAAjU0NES1HIDOCzsRuuCCC2IRBwDYVnZ2dlTLAeg8BksDQILl5+fL6XTKsqyQ5y3LksvlUn5+fpwjA9IfiRAAJJjD4VB5ebkkBSVDvtdlZWUMlAZigEQIAJJAYWGhKisrNWzYsIDjTqdTlZWVrCMExEjYs8bSCbPGACQbVpYGOpbQWWMAgNhxOBxMkQfiiK4xAABgWyRCAADAtkiEAACAbZEIAQAA2yIRAgAAtsWsMQBxxxRxAMmCRAhAXHk8HpWWlmrnzp3+Y06nU+Xl5SwaCCDu6BoDEDcej0dFRUUBSZAk1dfXq6ioSB6PJ0GRAbArEiEAceH1elVaWqpQi9n7ji1cuFBerzfeoQGwMRIhAHFRU1MT1BLUljFGdXV1qqmpiWNUAOyORAhAXDQ0NES1HABEA4kQgLjIzs6OajkAiAYSIQBxkZ+fL6fTKcuyQp63LEsul0v5+flxjgyAnaVdIrR//36dc845GjdunM444ww9/PDDiQ4JgL7YVb28vFySgpIh3+uysjLWEwIQV2mXCPXu3VsbNmzQ9u3b9eqrr+ruu+/W3r17Ex0WAEmFhYWqrKzUsGHDAo47nU5VVlayjhCAuEu7BRUdDodOOukkSVJLS4uMMSGn6wJIjMLCQs2ePZuVpQEkhaRrEdqwYYNmzZqloUOHyrIsrVmzJqhMRUWFcnNz1aNHD02YMEF///vfA87v379fY8eOldPp1G233aYBAwbEKXoAneFwOFRQUKA5c+aooKCAJAhAwiRdItTc3KyxY8eqoqIi5PknnnhCixYt0tKlS/Xaa69p7Nixmj59uj7++GN/mb59++of//iHamtrtWrVKu3atSte4QMAgBRimSTuN7IsS0899ZQuu+wy/7EJEyYoLy9PDz74oCSptbVVLpdL3/72t7V48eKga9x888362te+pqKioqBzTU1NysrKUl1dnfr06eM/3r17d3Xv3j36FQIAAGFraWlRS0uL/3VTU5NcLpcOHDgQ8PzuiqRrETqeI0eOaOvWrZo6dar/WLdu3TR16lRt2rRJkrRr1y4dPHhQknTgwAFt2LBBp5122nGv63K5lJWV5f9ZtmxZ7CoBAADCsmzZsoDntMvlitq1U2qw9J49e+T1ejV48OCA44MHD9a//vUvSdKHH36ob3zjG/5B0t/+9rd15plnHve6oVqEAABAcrjjjju0aNEi/2tfi1A0pFQi1Bnnnnuutm/fHtZ7+vTpE3HTGgAAiI1YDllJqa6xAQMGyOFwBA1+3rVrl4YMGZKgqAAAQKpKqUQoIyNDZ599ttavX+8/1traqvXr12vixIkJjAwAAKSipOsaO3TokHbs2OF/XVtbq+3bt6tfv34aPny4Fi1apOLiYp1zzjk699xzVVZWpubmZs2fPz+BUQNAYni9XhanBCKQdInQli1bNHnyZP9r3+Co4uJiPfLII7ryyiu1e/du/fCHP1RjY6PGjRun559/PmgANQCkO4/Ho9LSUu3cudN/zOl0qry8nO1KgE5K6nWEYs23jtCoUaPkcDhUUlKikpKSRIcFpBxaJeLP4/GoqKgoaAsh3wa27N2GdFRRUaGKigp5vV69++67UVlHiEQoKysqf5GAXdEqEX9er1e5ubkBf+dtWZYlp9Op2tpaElKkpWg+v1NqsDSA5OJrlTj2gVxfX6+ioiJ5PJ4ERZbeampq2k2CJMkYo7q6OtXU1MQxKiA1kQgB6BKv16vS0tKgrhlJ/mMLFy6U1+uNd2hpr6GhIarlADsjEQLQJbRKJE52dnZUywF2RiIEoEtolUic/Px8OZ1O/8DoY1mWJZfLpfz8/DhHBqQeEiEAXUKrROI4HA6Vl5dLUlAy5HtdVlbGQGmgE0iEAHQJrRKJVVhYqMrKSg0bNizguNPpZOo8EAamz7OOENBlvlljkgIGTbOWTfywhhPshHWEoox1hIDIhVpHyOVyqaysjCQIQExE8/lNIkQiBESMVgkA8RTN53fS7TUGIPU4HA4VFBQkOgwACBuDpQEAgG2RCAEAANsiEQIAALZFIgQAAGyLRAgAANgWiRAAALAtEiFJeXl5crvdqqioSHQoAACgHRUVFXK73crLy4vaNVlQkQUVAQBIKdF8ftMiBAAAbItECAAA2BaJEAAAsC0SIQAAYFtsugogCLvJA7ALEiEAATwej0pLS7Vz507/MafTqfLychUWFiYwMkQDSS4QiK4xAH4ej0dFRUUBSZAk1dfXq6ioSB6PJ0GRIRo8Ho9yc3M1efJkXX311Zo8ebJyc3O5r7A11hFiHSFA0hctBbm5uUFJkI9lWXI6naqtraUFIQX5ktxjf+VbliVJqqyspMUPKYN1hABEXU1NTbtJkCQZY1RXV6eampo4RoVo8Hq9Ki0tDUqCJPmPLVy4UF6vN96hAQlHIiS22AAkqaGhIarlkDxIcpEuYrHFBoOlJW3evJmuMdhednZ2VMsheZDkIl2UlJSopKTE3zUWDbQIAZAk5efny+l0+seMHMuyLLlcLuXn58c5MkSKJBdoH4kQAEmSw+FQeXm5JAUlQ77XZWVlDJROQSS5QPtIhAD4FRYWqrKyUsOGDQs47nQ6mVWUwkhygfYxfZ7p80AQFt1LT6EWy3S5XCorKyPJRUqJ5vObRIhECICNkOQiHUTz+c2sMQCwEYfDoYKCgkSHASQNxggBAADbIhECAAC2RSIEAABsi0QIAADYFokQAACwLRIhsekqAACpIBabrrKOEOsIwQZYOwZAOmEdIQCdFmo1YafTqfLyclYTBmB7dI0Baczj8aioqCggCZKk+vp6FRUVyePxJCgyAEgOJEJAmvJ6vSotLVWo3m/fsYULF8rr9cY7NABIGiRCQJqqqakJaglqyxijuro61dTUxDEqAEguJEJAmmpoaIhqOQBIRyRCQJrKzs6OajkASEckQkCays/Pl9PplGVZIc9bliWXy6X8/Pw4RwYAyYNECEhTDodD5eXlkhSUDPlel5WVsZ4QAFsjEQLSWGFhoSorKzVs2LCA406nU5WVlawjBMD2WFmalaVhA6wsDSCdsLI0YHPhJjYOh0MFBQXxCxAAUgSJEJBi2DIDAKKHMUJACmHLDACILhIhSXl5eXK73aqoqEh0KEC72DIDieD1elVdXa3HHntM1dXVfL+QUBUVFXK73crLy4vaNRkszWBppIjq6mpNnjy5w3JVVVWMB0JU0A2LZBXN5zctQkCKYMsMxBPdsLALEiEgRbBlBuKFbljYCYkQkCLYMgPxUlNTE9QS1JYxRnV1daqpqYljVEBskAgBKYItMxAvdMPCTkiEgBTClhmIB7phYSfMGmPWGFIQW2Yglrxer3Jzc1VfXx9ynJBlWXI6naqtreV7h4Rgiw3A5tgyA7Hk64YtKiqSZVkByRDdsEg3dI0BAILQDQu7oGuMrjEAaBfdsEhGdI0BAOKCblikO7rGAACAbZEIAQAA2yIRAgAAtkUiBAAAbIvB0kCSYHYOAMQfiRCQBDwej2655RbV19f7jzmdTpWXl7NeCwDEEF1jQIJ5PB4VFRUFJEGSVF9fr6KiInk8ngRFBgDpj0QISCCv16vS0tKQ+zn5ji1cuFBerzfeoQGALZAIScrLy5Pb7VZFRUWiQ4HN1NTUaOfOne2eN8aorq5ONTU1cYwKAJJTRUWF3G638vLyonZNxghJ2rx5M1tsICEaGhqiWg4A0llJSYlKSkr8W2xEAy1CQAJlZ2dHtRwAIDwkQkAC5efny+l0yrKskOcty5LL5VJ+fn6cIwMAeyARAhLI4XCovLxckoKSId/rsrIy1hMCgBghEQISrLCwUJWVlRo2bFjAcafTqcrKStYRAoAYskyoebs24RtsdeDAAQZLI+FYWRoAOieaz29mjQFJwuFwqKCgINFhABEhoUeqIRECAEQFW8UgFTFGCAAQMbaKQaoiEQIARIStYpDKSIQAABFhqxikMsYIATHCoFHYBVvFIJWRCAExwKBR2AlbxSCV0TUGRBmDRmE3bBWDVEYiBEQRg0ZhR2wVg1RGIgR0UlNTkyzLkmVZeu6550ImMwwahV35tooZOnRowHG2ikGyIxECOsHj8cjtdvtfz5w5U7m5uUHdXAwahZ0VFhbqww8/VFVVlVatWqWqqirV1taSBCGpMVga6IBvzM+x3V2+MT9t/2+XQaOwO7aKQaph01U2XcVxeL1e5ebmttvdZVmWnE6namtr5XA4/OXr6+tDjhM6tjwAIHzRfH7TNQYcR7hjfhg0CgCphUQIOI6ujPlh0CgApA7GCAHH0dUxP4WFhZo9ezYrSwNAkmOMEGOEcByM+QGA5MMYISBOGPMDAOmNRAjogG/Mz7BhwwKOM+YHAFIfXWN0jaGT2E0eAJJDNJ/fDJYGOomF4oDo8j3MJGnt2rWaNm0a/3OBuKNrDAAQd53dtgaINRIhSXl5eXK73aqoqEh0KACQ9nzb1tTX1wcc921bQzKE9lRUVMjtdisvLy9q12SMEGOEACBuwt22BgiF6fNAFDQ1NcmyLFmWpeeee05erzfRIQFpL9xta4BYIxGCLTE+AUiMrmxbA8QSiRBsh/EJQOJ0ddsaIFYYI8QYIVthfAKQWGxbg2hgjBDQRYxPABKLbWuQbEiEYCuMTwASz7dtzdChQwOOs20NEoGVpWErjE8AkkNhYaFmz57NtjVIOMYIMUbIVhifAACpjzFCQBcxPgEA0BaJEGyH8QkAAB+6xugaS1qx3pna6/UyPgEAUhBdY0h78Vj52eFwqKCgQHPmzFFBQQFJEADYEIkQkg4rPwMA4oVECEnF6/WqtLQ05Iwu37GFCxeyQSoAICpIhJBUWPkZQHuam5tlWZYsy1Jzc3Oiw0GaIBFCUolk5eempib/L8nnnnuOViMgzbT9N71hwwb+jSMqSISQVLq68nM8BlcDSBz+jSNWmD7P9Pmk0pWVn32Dq48t71sgkbWBgNTGv3Eci+nzSFvhrvzM4GogvfFvHLFGIoSkE87KzwyuBtIb/8YRayRCaS5VZ1kUFhbqrbfe8r9eu3atamtrg5q/IxlcDSD58W8csUYilOZiOcsi1klWnz59ZIyRMUYXXXRRyJWfuzq4GkBq4N84Yo1EKI3FepZFMkxlzc/Pl9PpDBpP5GNZllwul/Lz8+McGYBo4N84Yo1EKE3FepuKZJnKGu7gagCphX/jiDUSoTQU61kWybYXWDiDqwGknq7+G0/VMZKIL9YRSsN1hKqrqzV58uQOy1VVVamgoCDgWHNzs3r16iVJOnTokDIzMwPO+9b5aW8WR6h1fuLFdz+lLwZXT5s2jf9LBNJIuP/G+Z2QvlhHCMcVySyLjsb9JPNU1ra/4CZNmsQvPCDNdGYChU+ydN8j+ZEIpZjONPXGcpuKSJKsWDdTZ2Zm+n9JHtuSBcA+kq37HsmNRCjFdGamVldmWXT2F0ckU1mTYZYZgPTGStQIF4lQCulsU28st6no6lRWmqkBxEMyd98jOZEIpYhwm3pjtU1FV6ay0kwNIF5YiRrhIhFKAV1t6o3VNhXhJFk0UwOIJ1aiRrhIhBKsMwOII2nqjdU2FZ1NsmimBhBPrESNcJEIRVksZkbFuqm3q784OpNk0UwNIJ5YiRrhIhGKsnBnRnWmfKybemP5i4NmagDxxmrzCIuxsQMHDhhJ5sCBA1G53urVq82wYcOMJP+P0+k0q1evjqj80aNHjdPpNJZlBZT1/ViWZVwulzl69GjU43e5XO3G3xnxih0AjuX7HS/JrF27lt8zaSSaz28SoSj9Ra5evTrkw96yLGNZVlAy0dXyx76nvfJdFYtfHPGKHQDaOnTokP/3zaFDhxIdDqIoms9v9hqLwl4l4e6/1dX9ujwej2655ZaAaegul0tlZWVRa+rtaK+xropH7AAAe4jmXmMkQlH4iwx3k9NINkVN5U0EUzl2AEDyiGYidEKUYrK1cGdGRTKTyjdTKxWlcuwAgPSUdrPG6urqVFBQILfbrTFjxujJJ5+M+WeGOzOKmVQAACSHtOsaa2ho0K5duzRu3Dg1Njbq7LPP1rvvvhtyrEu0xwjV19eHbPFob4xQZ8sDAID/E82usbRrEcrOzta4ceMkSUOGDNGAAQO0b9++mH5muOvwsOAXAADJIekSoQ0bNmjWrFkaOnSoLMvSmjVrgspUVFQoNzdXPXr00IQJE/T3v/895LW2bt0qr9crl8sV46jDX8CLBb8AILnEYmcAJL+kS4Sam5s1duxYVVRUhDz/xBNPaNGiRVq6dKlee+01jR07VtOnT9fHH38cUG7fvn2aN2+efvvb38YjbEmd33+rq+UBAEB0JfUYIcuy9NRTT+myyy7zH5swYYLy8vL04IMPSpJaW1vlcrn07W9/W4sXL5YktbS06MILL9QNN9yguXPntnv9aPYx+oS7Dk+s1u0BAISHJT5Sh23HCB05ckRbt27V1KlT/ce6deumqVOnatOmTZK+2M38uuuu09e+9rXjJkFtNTU1Bfy0tLR0OcbMzEz/RqSdSWrCLQ8AiD6PxyO32+1/PXPmTOXm5srj8SQwKvi0tLQEPaujJaUSoT179sjr9Wrw4MEBxwcPHqzGxkZJ0iuvvKInnnhCa9as0bhx4zRu3Di98cYbx72uy+VSVlaW/2fZsmUxqwMAILl4PB4VFRUFrHwvSfX19SoqKiIZSgLLli0LeE5Hc+xv2i2oeP7556u1tTWs99TV1QU0rXXv3j3aYQEAkpDX61VpaWnIpUyMMbIsSwsXLtTs2bPpJkugO+64Q4sWLfK/bmpqiloylFKJ0IABA+RwOLRr166A47t27dKQIUO6fN0+ffpEbYwQACB11NTUtLvvo/RFMlRXV6eampqgLY8QP927d49ZI0VKdY1lZGTo7LPP1vr16/3HWltbtX79ek2cODGBkQEAUlEkWx4hPSRdi9ChQ4e0Y8cO/+va2lpt375d/fr10/Dhw7Vo0SIVFxfrnHPO0bnnnquysjI1Nzdr/vz5CYwaAJCK2PIISTd9vr2d2YuLi/XII49Ikh588EHde++9amxs1Lhx4/TAAw9owoQJYX9WLKbPAwBSB1sepaZoPr+TLhGKJxIhAIBv1pikgGTIt+URq/0nH9uuIxQreXl5crvd7a5mDQBIX2x5lDoqKirkdruVl5cXtWvSIkSLEABArCydSmgRAgAgytomPZMmTSIJsgkSIQAAYFtJN30eAIBE8O39CHuhRQgAgC5obm6WZVmyLEvNzc2JDgddRCIEAABsi0QIAADYFomQWEcIABA+r9fr//OGDRsCXodCV1rkWEcoylhHCADQFR6PR7fccovq6+v9x5xOp8rLy9tdgLG5uVm9evWS9MW+mpmZmXGJNR2xjhAAAAni25KjbRIkSfX19SoqKpLH40lQZOgKEiEAADrJ6/WqtLQ05DR737GFCxeG7CYLtysN8UEiBABAJ9XU1Gjnzp3tnjfGqK6uTjU1NQHHPR6P3G63//XMmTOVm5tL61ESIBECAKCTGhoawi5HV1pyIxECAKCTsrOzwyoXSVca4oNECACATsrPz5fT6ZRlWSHPW5Yll8ul/Px8SV3vSkP8kAgBANBJDodD5eXlkhSUDPlel5WV+Xeu70pXGuKLRAgAgDAUFhaqsrJSQ4cODTjudDpVWVkZsI5QuF1piD8WVMzK0qhRo+RwOFRSUqKSkpJEhwUASAG+Z4gkrV27VtOmTfO3BPl4vV7l5uaqvr4+5Dghy7LkdDpVW1sb9F4Eq6ioUEVFhbxer959992oLKhIIsTK0gCALujsStG+WWOSApIhX1fasa1I6BgrSwMAkCLC6UpD/NEiRIsQACAOOtOVhs6hRQgAgBTTNumZNGkSSVCSIBECAAC2RSIEAABs64REBwAAgB1kZmaGnEKPxKJFCAAA2BaJEAAASai5uVmWZcmyLDU3Nyc6nLRFIgQAAGyLREhSXl6e3G63KioqEh0KAABhs0vrUUVFhdxut/Ly8qJ2TRZUZEFFAEAS6uwWHuGWTQcsqAgAQJrzer3+P2/YsCHgNaKHRAgAgCTj8Xjkdrv9r2fOnKnc3Fx5PJ4ERpWeSIQAAEgivt3q6+vrA47X19erqKiIZCjKSIQAAEgSXq9XpaWlIRde9B1buHBhUDcZ3WhdRyIEAECSqKmp0c6dO9s9b4xRXV2dampq/MfoRosMiRAAAEmioaEhrHJ0o0WORAgAgCSRnZ3d6XJd7UZDIBIhAACSRH5+vpxOpyzLCnnesiy5XC7l5+d3qRsNwUiEAABIEg6HQ+Xl5ZIUlAz5XpeVlcnhcITdjdaWXVai7gwSIQAAkkhhYaEqKys1dOjQgONOp1OVlZUqLCyUFF43GtpHIiT2GgMAJJfCwkK99dZb/tdr165VbW2tPwmSwutGSxfsNRZl7DUGAEhWndk/zDdrTFLAoGlfctS2BSncaycz9hoDAACd7kZD+0iEAABIYZ3pRjsWK1H/HxIhAABSnMPh8P950qRJAa+PxUrUgUiEAABIQpmZmTLGyBgTtTE8rEQdjEQIAAAbYCXq0EiEAACwAVaiDu2ERAcAAAAi4+tGO55IVqJOZ7QIAQBgA6xEHRqJEAAANmDHlag7g0QIAAAbCGdDVzshEQIAwCZYiToYe42x1xgAwGZ8zz/pi5Wop02bllItQew1BgAAuiyclailLzZptSxLlmWpubk51uHFFYkQAACwLRIhSXl5eXK73aqoqEh0KAAAoB0VFRVyu93Ky8uL2jUZI8QYIQCAzTQ3N6tXr16SpEOHDnW4l1m45WONMUIAAABRwBYbAADYTGe25LALWoQAAMBxtd2RfsOGDWm1Qz2JEAAAaJfH45Hb7fa/njlzpnJzc+XxeBIYVfSQCAEAgJA8Ho+KiopUX18fcLy+vl5FRUVpkQyRCAEAgCBer1elpaUhxxL5ji1cuDDlu8lIhAAAQJCamhrt3Lmz3fPGGNXV1ammpiaOUUUfiRAAAAjS0NAQ1XLJikQIAAAEyc7Ojmq5ZEUiBAAAguTn58vpdMqyrJDnLcuSy+VSfn5+nCOLLhIhAAAQxOFwqLy8XJKCkiHf67Kysg53rk92JEIAACCkwsJCVVZWaujQoQHHnU6nKisrVVhYGPJ9zc3NsixLlmWpubk5HqF2GVtsAACAdhUWFmrq1KnKysqSJK1du1bTpk1L+ZYgH1qEAADAcbVNeiZNmpQ2SZBEIgQAAGyMRAgAANgWiRAAALAtEiEAAGBbJEKS8vLy5Ha7VVFRkehQAABIeW03Yt2wYUPUNmatqKiQ2+1WXl5eVK4nSZYJta2sTTQ1NSkrK0sHDhxQnz59Eh0OAABJqbm5Wb169ZIkHTp0SJmZme2W9Xg8uuWWW1RfX+8/5nQ6VV5e3u66Q+GK5vObFiEAAHBcmZmZMsbIGNNhElRUVBSQBElSfX29ioqK5PF4Yh1q2EiEAABAxLxer0pLSxWqo8l3bOHChVHrJosWEiEAABCxmpoa7dy5s93zxhjV1dWppqYmjlF1jEQIAABErKGhIarl4oVECAAARCw7Ozuq5eKFRAgAAEQsPz9fTqdTlmWFPG9Zllwul/Lz8+Mc2fGRCAEAgIg5HA6Vl5dLUlAy5HtdVlaWdBu2kggBAICoKCwsVGVlpYYOHRpw3Ol0qrKyMmrrCEUTCyqyoCIAAFHle75K0tq1azVt2rSotgSxoCIAAEhabZOeSZMmJV13WFskQgAAwLZIhAAAgG2RCAEAANsiEQIAALZFIgQAAGyLRAgAANgWiRAAALAtEiEAAGBbJyQ6AAAAkF4yMzOVKhtX0CIEAABsi0QIAADYFokQAACwLRIhAABgWyRCAADAtkiEAACAbZEIAQAA2yIRAgAAtkUiJCkvL09ut1sVFRWJDgUAALSjoqJCbrdbeXl5UbumZVJl6ccYaGpqUlZWlg4cOKA+ffokOhwAANAJ0Xx+0yIEAABsi0QoAi0tLbrzzjvV0tKS6FDiwk71pa7py071pa7py071jXVd6RqLoGnNbl1rdqovdU1fdqovdU1fdqpvqLrSNRZDsR4wHe71k618LK+dbOXDFct47FTXeJSP5bWTrXy4+B6nTvlYXjvZyseUsbEDBw4YSebAgQP+Y6NHj47o/R0J5/rJVj7c+iZT7OGWT7Z7a6e6xro83+PExWOnusa6vN2/x1253+05IUH5V1Iw/9sr2NTU5D/m9XoDXh+Pr1xny4d7/WQrH259kyn2cMsn2721U11jXZ7vceLisVNdY13e7t9j359NFEb32HqM0M6dO+VyuRIdBgAA6IK6ujo5nc6IrmHrRKi1tVUfffSRevfuLcuyEh0OAADoBGOMDh48qKFDh6pbt8iGO9s6EQIAAPbGrDEAAGBbJEIAAMC2SIQ6ITc3V5ZlBf2UlJRIkj777DOVlJSof//+6tWrl/7rv/5Lu3btSnDUXdNRXQsKCoLO3XjjjQmOumu8Xq9+8IMf6JRTTlHPnj01cuRI/eQnPwmYhWCM0Q9/+ENlZ2erZ8+emjp1qt57770ERt01nanrddddF3RvZ8yYkcCoI3Pw4EEtXLhQOTk56tmzp8477zxt3rzZfz5d7q3UcV1T+d5u2LBBs2bN0tChQ2VZltasWRNwvjP3cd++fbrmmmvUp08f9e3bV9dff70OHToUx1p0TjTqGup3+D333BPHWnROR3X1eDyaNm2a+vfvL8uytH379qBrRO3ZG/EEfBv4+OOPTUNDg//nxRdfNJJMVVWVMcaYG2+80bhcLrN+/XqzZcsW85WvfMWcd955iQ26izqq6wUXXGBuuOGGgDLRWMchEX7605+a/v37mz//+c+mtrbWPPnkk6ZXr16mvLzcX+aee+4xWVlZZs2aNeYf//iHufTSS80pp5xiPv300wRGHr7O1LW4uNjMmDEj4N7u27cvgVFH5oorrjBut9u8/PLL5r333jNLly41ffr0MTt37jTGpM+9NabjuqbyvV27dq1ZsmSJ8Xg8RpJ56qmnAs535j7OmDHDjB071vztb38zNTU15tRTTzVz5syJc006Fo265uTkmB//+McB9/rQoUNxrknHOqrrihUrzI9+9CPz8MMPG0lm27ZtQdeI1rOXRKgLSktLzciRI01ra6vZv3+/OfHEE82TTz7pP//2228bSWbTpk0JjDI62tbVmC8SodLS0sQGFSUXX3yxWbBgQcCxwsJCc8011xhjjGltbTVDhgwx9957r//8/v37Tffu3c1jjz0W11gj1VFdjfniYTl79uw4RxYbhw8fNg6Hw/z5z38OOD5+/HizZMmStLq3HdXVmPS5t8c+MDtzH9966y0jyWzevNlf5rnnnjOWZZn6+vq4xR6urtTVmC8Sofvvvz+OkUYuVCLkU1tbGzIRiuazl66xMB05ckR//OMftWDBAlmWpa1bt+rzzz/X1KlT/WW+/OUva/jw4dq0aVMCI43csXX1efTRRzVgwACdccYZuuOOO3T48OEERtl15513ntavX693331XkvSPf/xDGzdu1EUXXSRJqq2tVWNjY8C9zcrK0oQJE1Lu3nZUV5/q6moNGjRIp512mm666Sbt3bs3EeFG7OjRo/J6verRo0fA8Z49e2rjxo1pdW87qqtPutzbtjpzHzdt2qS+ffvqnHPO8ZeZOnWqunXrpldffTXuMXdVON/Ze+65R/3799dZZ52le++9V0ePHo13uDEXzWevrVeW7oo1a9Zo//79uu666yRJjY2NysjIUN++fQPKDR48WI2NjfEPMIqOraskXX311crJydHQoUP1+uuv6/bbb9c777wjj8eTuEC7aPHixWpqatKXv/xlORwOeb1e/fSnP9U111wjSf77N3jw4ID3peK97aiukjRjxgwVFhbqlFNO0fvvv6/vf//7uuiii7Rp0yY5HI4ERh++3r17a+LEifrJT36i0aNHa/DgwXrssce0adMmnXrqqWl1bzuqq5Re97atztzHxsZGDRo0KOD8CSecoH79+qXUve7sd/aWW27R+PHj1a9fP/31r3/VHXfcoYaGBt13331xjTfWovnsJREK0+9+9ztddNFFGjp0aKJDiblQdf3GN77h//OZZ56p7OxsTZkyRe+//75GjhyZiDC77E9/+pMeffRRrVq1Sqeffrq2b9+uhQsXaujQoSouLk50eFHVmbpeddVV/vJnnnmmxowZo5EjR6q6ulpTpkxJVOhdtnLlSi1YsEDDhg2Tw+HQ+PHjNWfOHG3dujXRoUVdR3VNt3uL9i1atMj/5zFjxigjI0Pf/OY3tWzZMnXv3j2BkSUvusbC8OGHH2rdunX6+te/7j82ZMgQHTlyRPv37w8ou2vXLg0ZMiTOEUZPqLqGMmHCBEnSjh074hFWVN12221avHixrrrqKp155pmaO3eubr31Vi1btkyS/Pfv2FkIqXhvO6prKCNGjNCAAQNS8t5K0siRI/Xyyy/r0KFDqqur09///nd9/vnnGjFiRFrdW+n4dQ0l1e+tT2fu45AhQ/Txxx8HnD969Kj27duXUve6q9/ZCRMm6OjRo/rggw9iGV7cRfPZSyIUhuXLl2vQoEG6+OKL/cfOPvtsnXjiiVq/fr3/2DvvvKP//Oc/mjhxYiLCjIpQdQ3FN6UxOzs7DlFF1+HDh4OWZnc4HGptbZUknXLKKRoyZEjAvW1qatKrr76acve2o7qGsnPnTu3duzcl721bmZmZys7O1ieffKK//OUvmj17dlrd27ZC1TWUdLm3nbmPEydO1P79+wNaAl966SW1trb6/0cuFXT1O7t9+3Z169YtqHsw1UX12Rvu6G678nq9Zvjw4eb2228POnfjjTea4cOHm5deesls2bLFTJw40UycODEBUUZHe3XdsWOH+fGPf2y2bNliamtrzdNPP21GjBhhJk2alKBII1NcXGyGDRvmn1Lu8XjMgAEDzPe+9z1/mXvuucf07dvXPP300+b11183s2fPTskp1h3V9eDBg+a73/2u2bRpk6mtrTXr1q0z48ePN1/60pfMZ599luDou+b55583zz33nPn3v/9tXnjhBTN27FgzYcIEc+TIEWNM+txbY45f11S/twcPHjTbtm0z27ZtM5LMfffdZ7Zt22Y+/PBDY0zn7uOMGTPMWWedZV599VWzceNG86UvfSkpp89HWte//vWv5v777zfbt28377//vvnjH/9oBg4caObNm5fIaoXUUV337t1rtm3bZp599lkjyTz++ONm27ZtpqGhwX+NaD17SYQ66S9/+YuRZN55552gc59++qm5+eabzcknn2xOOukkc/nllwfcrFTTXl3/85//mEmTJpl+/fqZ7t27m1NPPdXcdtttKbuOUFNTkyktLTXDhw83PXr0MCNGjDBLliwxLS0t/jKtra3mBz/4gRk8eLDp3r27mTJlSsjvQLLrqK6HDx8206ZNMwMHDjQnnniiycnJMTfccINpbGxMcORd98QTT5gRI0aYjIwMM2TIEFNSUmL279/vP58u99aY49c11e9tVVWVkRT0U1xcbIzp3H3cu3evmTNnjunVq5fp06ePmT9/vjl48GACanN8kdZ169atZsKECSYrK8v06NHDjB492tx9991JmfB2VNfly5eHPL906VL/NaL17GXTVQAAYFuMEQIAALZFIgQAAGyLRAgAANgWiRAAALAtEiEAAGBbJEIAAMC2SIQAAIBtkQgBAADbIhECAJt65JFHZFmWqqurO1W+oKBAubm5MY0JsWdZlv9n6tSpiQ6nQ48//nhAzI888khUr08iBAAprrq6OuBB4XA4dPLJJ+uMM85QcXGxnn/+eaXaJgL79+/XnXfe2ekkLdZyc3NlWZb69++vlpaWkGVmz57tvwfJvtt7fn6+Vq5cqTvuuCPk+X/84x+64YYbNGrUKGVmZqpHjx7KyclRYWGhVq5cqc8//zzsz/znP/8py7J02WWXHbfc8uXLZVmW7r77bklfbJy7cuVKff/73w/7MzvjhJhcFQAQd3PmzNHMmTNljNHBgwf1zjvvaM2aNVqxYoWmTp2qJ598Un379vWXnzt3rq666iplZGQkLuh27N+/Xz/60Y8kfdESlQx69Oihffv26ZlnntH/+3//L+Dcrl27tHbtWvXo0UOfffZZgiLsvBEjRujaa68Nee7uu+/WD37wA5188sm68sordcYZZygjI0M7d+7UunXrNG/ePG3cuFG/+c1vwvrM008/XRMmTNCzzz6rjz/+WIMGDQpZbvny5XI4HCouLpYk5eTkKCcnR9XV1f7kKJpIhAAgTYwfPz7o4Xbffffpe9/7nu677z7NmTNHzz33nP+cw+GQw+GId5gpa+TIkerWrZuWL18elAitWLFCkjRr1iw9+eSTiQgvKh555BEtWbJEU6ZM0erVq5WVlRVwfunSpdq+fbteeeWVLl3/+uuv16uvvqo//vGPWrRoUdD5HTt2qKamRjNnztSwYcO69BnhomsMANKYw+HQL3/5S51//vl6/vnntXHjRv+5cMcItee9997T3LlzlZ2drYyMDOXm5uq2225Tc3NzQLnrrrtOlmVp9+7dmjdvnvr376/MzExNmTJFr732mr9cdXW1TjnlFEnSj370I393U9vxSUePHtXPfvYzud1u9ejRQ/3799fll1+uN954I+AzP/jgA1mWpTvvvFN//vOflZeXpx49eig7O1u33Xabjh49GlZd58+frxdeeEEfffRRwPHly5fr4osvDtnK8dFHH+k73/mOxo0bp5NPPlk9evSQ2+3Wz372M3m93oCyn332me68806ddtppOumkk9S3b1+deeaZuu222wLKPfvss7rgggs0YMAA9ezZU8OHD1dhYaHefffdsOrT1pEjR/T9739fvXv31p/+9KegJMhn3LhxKikpCTq+bt06TZs2TX379lWPHj00ZswYPfTQQwFlrrrqKmVmZmr58uUhr+07vmDBgi7XI1wkQgBgA9dff72kLx6g0bR161adc8452rBhg775zW+qoqJCl1xyiR544AFdeOGFIceSzJgxQw0NDbrzzju1cOFCbdmyRRdccIHefPNNSdLo0aN1//33S5Iuv/xyrVy5UitXrlRZWZn/Gtdcc40WL14sp9Ope++9VzfeeKOqqqo0ceJEbdu2Legz165dqwULFuiiiy7S/fffr7Fjx+oXv/iFfv7zn4dV32uvvVbdunXTH/7wB/+xv/3tb3r77bfbfXi//vrr8ng8+trXvqa77rpL99xzj4YPH67Fixfr5ptvDihbUlKiH/3oR/rKV76i+++/Xz/96U81ZcoUvfTSS/4yL7/8si699FLt379fd9xxhx588EHdcMMN2rt3r3bs2BFWfdp65ZVX1NDQoMsvv1z9+vUL672//e1vNW3aNB06dEhLlizRfffdp5EjR+qmm24KSOJ69+6toqIivfnmm9q8eXPANVpbW7VixQoNHDhQl156aZfrETYDAEhpVVVVRpK599572y2zdetWI8kUFhb6jy1fvtxIMlVVVZ36nAsuuMDk5OQEHBszZow57bTTTFNTU8Bxj8djJJnly5f7jxUXFxtJ5vLLLzetra3+41u2bDGWZZnp06f7j9XW1hpJZunSpUFxvPDCC0aSueKKKwKus337duNwOMz5558fdJ2TTjrJ1NbW+o+3traa008/3QwZMqRTdc/JyTGnn366McaYwsJCM2rUKP+5G264wQwZMsR8/vnnpqSkxEgK+KzDhw8HxOlz7bXXmm7dupmPPvrIf+zkk082F1100XFjufXWW40ks2vXrk7FfixJpri4OOj4Aw88YCSZ++67L+hcU1OT2b17t/9nz549/nMfffSR6d69u5kzZ07Q+2655RbTrVs38/777/uPbdiwwUgyN910U0DZ559/3kgyt956a8i4fd/ztt+paKBFCABsoE+fPpKkpqamqF3zjTfe0Ouvv66rr75aLS0t2rNnj//n/PPPV2Zmpl544YWg933ve9+TZVn+12effbYuvPBCrVu3TocOHerwc5966ilJ0pIlSwKuM3bsWM2aNUsbN27U7t27A95z2WWXBXStWZalyZMnq7GxsVOf2daCBQv07rvv6pVXXtGnn36qJ554QnPnztUJJ4QedtuzZ09/nEeOHNG+ffu0Z88eTZ8+Xa2trdqyZYu/bFZWlv75z3/6W8dC8XVZrV69OuyuvePxfTd835W25s+fr4EDB/p/cnJy/OcqKyvV0tKi66+/PuA7sGfPHs2aNUutra1at26dv3x+fr5GjRqlxx57LGBgua9bzNd6GS8kQgBgA8d7yPl4vV41NjYG/Bw4cKDd8m+//bakLwbQtn1IDhw4UIMGDVJzc7N27doV9L7Ro0cHHXO73fJ6vfrwww87rEttba26desW8jqnn366v0xbI0aMCCrbv39/SdLevXs7/My2ZsyYoezsbC1fvlyVlZVqamrS/Pnz2y1/9OhR3XXXXRo1apR/PNPAgQM1d+5cSdInn3ziL1tWVqZPPvlEZ555pkaOHKmvf/3revrpp9Xa2uov861vfUtnnXWWbr75ZvXr108zZ87UAw88EJT8het4yfLSpUv14osv6sUXX9SYMWMCzvm+B1OnTg36Hlx44YWSFPQ9WLBggfbv3+9Paj/55BOtWbNG5557rv8exguzxgDABl5//XVJ0mmnndZumbq6Ov8gZZ/i4uJ2F7Az/7s20Xe+8x3NmDEjZJmTTz65C9FG3/Fmx5kw11hyOByaN2+e/vu//1v//Oc/9ZWvfCVkUuazaNEi/epXv9KVV16pJUuWaNCgQTrxxBP12muv6fbbbw9IcmbPnq0PPvhAa9eu1csvv6x169bpd7/7nfLz87Vu3TplZGSof//+2rx5s2pqavTiiy9qw4YNuvXWW7V06VKtXbtWEydODKs+PmeccYYkafv27UHnzjzzTJ155pmSgu+p7+9vxYoVys7ODnntYxPRefPmacmSJVq+fLnmzJmjVatW+VuV4o1ECABs4He/+50k6eKLL263zJAhQ/Tiiy8GHBs6dGi75b/0pS9J+iIxCGeF4rfffltf+cpXAo699dZbcjgc/i6Xtl1exxoxYoRaW1v19ttvB7VOvPXWW5IUlNBF24IFC/Szn/1Mf/vb3/Tb3/72uGVXrlypSZMm6fHHHw843t7A5n79+unaa6/VtddeK2OMFi9erJ///Od6+umn/dP2HQ6HCgoK/Gssvf766zr77LN11113dXlA/Fe/+lUNGTJETz31lO677z5/i1lHfN+DAQMGdPp7kJ2drZkzZ+rZZ59VXV2dli9frpNOOklXXXVVl2KPBF1jAJDGvF6vvvvd72rjxo2aOXOmvvrVr7ZbtkePHpo6dWrAj9vtbrf8WWedpTPOOEMPPfSQ/v3vfwedP3r0qPbt2xd0/Oc//3lAK8xrr72mdevWacqUKerVq5ck+f8b6v2+lYmXLVsWcJ0333xTzzzzjM4//3wNHDiw3bijYdSoUSovL9fSpUt15ZVXHresw+EIanVqbm72z4zz8Xq92r9/f8Axy7J01llnSfq/v4s9e/YEfcaXv/xl9ezZM+TfV2dlZGTo7rvv1sGDB3XllVe22y16bF2uuOIKde/eXUuXLtWnn34aVP7AgQMhV+O+/vrr1draqu985zvaunWrioqKjtt1Gyu0CAFAmnjttdf0xz/+UZICVpb+8MMPNW3aNK1atSqqn2dZllauXKmvfe1rGjNmjBYsWKDTTz9dhw8f1o4dO+TxeLRs2TJdd911Ae/78MMPNX36dF166aVqaGjQgw8+qJ49e+ree+/1l+nfv79OPfVUPf744xo5cqQGDx6szMxMzZo1SxdeeKGuuOIKPf744/rkk090ySWXqLGxURUVFerRo4ceeOCBqNazPbfcckunyhUVFek3v/mNrrzySk2dOlW7du3S73//+6AWl4MHDyo7O1uXXnqpzjrrLA0aNEi1tbX69a9/rZNPPlmzZs2SJN1www3auXOnpk2bppycHP+A7YMHD2revHkR1Wn+/Pn66KOP9IMf/EAjR470ryx94oknqqGhQS+88II2btzo70aTJKfTqV//+tf6+te/rtGjR2vu3LnKycnR7t279cYbb2jNmjV66623gvapu/jiizV48GD/ApTxXDsoQFTnoAEA4s43rdj3061bN9OnTx/jdrvNvHnzzHPPPRfyfdGYPm+MMR988IH55je/aXJycsyJJ55o+vXrZ8aPH28WL15s/vOf//jL+abPf/zxx+baa681/fr1Mz179jSTJ082W7ZsCbruq6++as477zxz0kknGUkBn/3555+be+65x3z5y182GRkZ5uSTTzazZ882r7/+esA1jjcNf+nSpUFT3dvTdvr88YSaPt/c3Gy++93vmuHDh5vu3bubU0891SxbtsysW7cuYDp4S0uLWbx4scnLyzP9+vUzGRkZJicnx8yfP9+8++67/uutXr3azJo1ywwbNsxkZGSYAQMGmEmTJpnKysoO4zOm/enzbW3bts1cf/315tRTTzU9e/Y03bt3Ny6Xy1x22WVm5cqV5siRI0Hv2bhxo7nsssvMwIEDzYknnmiys7NNQUGB+cUvfmE+/fTTkJ/zve99z0gyI0eODLnEQFuxmj5vGZNiO/EBAFLSddddpz/84Q8ptwFsurEsS1dddZV+9atfKSMjIyHdUeE4cuSImpqa9Morr+iyyy7T8uXLg1oZI8EYIQAAbObxxx/XwIEDVVhYmOhQOuTxeDRw4MAOd63vKsYIAQBgI21nBnZ2ZlgiTZ48OSDmaK8zRCIEAICNhLPUQTIYPHiwBg8eHLPrM0YIAADYFmOEAACAbZEIAQAA2yIRAgAAtkUiBAAAbItECAAA2BaJEAAAsC0SIQAAYFskQgAAwLZIhAAAgG39fxGGkSmGfZmCAAAAAElFTkSuQmCC\n" }, "metadata": {} } ], "source": [ "# In case we ran over multiple files / batches, merge them together\n", "full_mass_list = ak.concatenate(mass_list)\n", "\n", "bin_edges = np.arange(start=70, # The interval includes this value\n", " stop=110, # The interval doesn't include this value\n", " step=1 ) # Spacing between values\n", "bin_centres = (bin_edges[:-1] + bin_edges[1:]) / 2 # central values of each bin\n", "\n", "# histogram the data\n", "data_x,_ = np.histogram(full_mass_list, bins=bin_edges )\n", "\n", "# statistical error on the data\n", "data_x_errors = np.sqrt(data_x)\n", "\n", "# get current axes\n", "main_axes = plt.gca()\n", "\n", "# plot the data points\n", "main_axes.errorbar(x=bin_centres,\n", " y=data_x,\n", " yerr=data_x_errors,\n", " fmt='ko', # 'k' means black and 'o' is for circles\n", " label = 'Data')\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('Di-lepton Mass [GeV]', fontsize=13, x=1, horizontalalignment='right')\n", "\n", "# y-axis label\n", "main_axes.set_ylabel('Events', fontsize=13, y=1, horizontalalignment='right')\n", "\n", "# make the y-axis log scale\n", "main_axes.set_yscale('log')\n", "\n", "# add minor ticks on x-axis for main axes\n", "main_axes.xaxis.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": "YhX4PJq56oe9" }, "source": [ "Compare this with \"Total Cross Section\" as a function of \"Center of Mass Energy\".\n", "\n", "\"Di-muon\n", "\n", "\n", "\n", "## Some calculations by hand\n", "\n", "1. Calculate the invariant mass for this event, in GeV.\n", " * lep_px[0] = 32387 MeV,\n", " * lep_px[1] = -18772 MeV,\n", " * lep_py[0] = 7047 MeV,\n", " * lep_py[1] = -9563 MeV,\n", " * lep_pz[0] = 144152 MeV,\n", " * lep_pz[1] = -2942 MeV,\n", " * lep_E[0] = 147913 MeV,\n", " * lep_E[1] = 21272 MeV\n", " * (can you check these values yourself?)\n", " * (can you check your calculated value?)\n", "2. If lep[0] in the first event is an electron, what's its Lorentz factor, $\\gamma$? Use an electron mass of 0.5 MeV in the formula $E = \\gamma m$\n", "3. Calculate lep_p[0], the magnitude of the total momentum of lep[0].\n", "4. Using lep_p[0], calculate the speed of this electron. Use the formula $p = \\gamma mv/c$ (if you would like, you can compare the momentum to the formula from Newton, $p=mv$, and compare the kinetic energy to $m v^2 / 2$, to see how important relativity is at these energies).\n", "\n", "## Concepts\n", "\n", "1. Write the possible Z decays into charged leptons. Give an argument as to whether or not each decay happens at the same rate.\n", "2. Besides charged leptons, what are the other possible decays of the Z?" ] }, { "cell_type": "markdown", "metadata": { "id": "GhPBvo166oe9" }, "source": [ "[Back to contents](#contents)" ] }, { "cell_type": "markdown", "metadata": { "id": "HU-Si7JM6oe9" }, "source": [ "We hope you enjoyed finding the Z boson yourself. For more analysis notebooks visit the [ATLAS Open Data website](https://opendata.atlas.cern) or directly our [educational notebooks section](https://opendata.atlas.cern/docs/13TeV25Doc/13tutorial). Happy analyzing!" ] }, { "cell_type": "markdown", "metadata": { "id": "sGOIpz-Z6oe-" }, "source": [ "## Going further\n", "\n", "If you want to go further, there are a number of things you could try:\n", "* **Estimate the mass of the Z boson** from your graph.\n", "* **Estimate the width of the Z boson** from your graph. Look up the real width, and see if you can explain why yours is different.\n", "* **Increase the fraction** of events processed. Remove the line `entry_stop=numevents*0.5` in '[Can we process the data yet?!](#load_data)' to process all the events in the file.\n", "* **Use even more data**. Remove the line `break` from the function in '[Can we process the data yet?!](#load_data)' to run through all the available files. We can use all the data ATLAS collected during 2015 and 2016. The analysis will take a bit longer to run!\n", "* **Make a fit to your graph**. See if you can automatically fit the peak position and width of the Z boson.\n", "* **Extract the Z boson mass from your fit**.\n", "* **Extract the Z boson width from your fit**.\n", "* **Extend your graph**. Try a larger range in mass. Can you see the Higgs boson? Why or why not? Can you see any particles with much lower mass (around 5-10 GeV)? Why or why not? Can you figure out what some of those particles might be?\n", "\n", "With each change, keep an eye on your graph.", "\n", "\n", "
\n", "We welcome your feedback on this notebook or any of our other materials! Please fill out this survey to let us know how we're doing, and you can enter a raffle to win some ATLAS merchandise!\n", "
" ] }, { "cell_type": "markdown", "metadata": { "id": "j8nKESDq6oe-" }, "source": [ "[Back to contents](#contents)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.12" }, "colab": { "provenance": [] } }, "nbformat": 4, "nbformat_minor": 0 }