Exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle physics data-analysis with CMS open data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welcome to the exercise where real data from CMS experiment at CERN is used for a simple particle physics data-analysis. The goal for the exercise is to discover the appearance of Z boson, determine the mass and the lifetime of Z and compare the results to the known values of these.\n",
    "\n",
    "In the exercise invariant mass values will be calculated for muon pairs that are detected in the CMS detector. A histogram will be made from the calculated invariant mass values. After that a Breit-Wigner fit will be made to the histogram. With the fitted Breit-Wigner function it will be possible to determine the mass and the lifetime of Z boson.\n",
    "\n",
    "In the end there will be also a quick look about how a pseudorapidity effects to the mass distribution of muon pairs.\n",
    "\n",
    "The structure of the exercise is following:\n",
    "- theory background\n",
    "- calculation of invariant masses\n",
    "- making the histogram\n",
    "- fitting the function to the histogram\n",
    "- analysing the histogram\n",
    "- looking the histogram of the whole range of data\n",
    "- the effect of pseudorapidity to the mass distribution\n",
    "\n",
    "Now take a relaxed position and read the theory background first. Understanding the theory is essential for reaching the goal and learning from the exercise. So take your time and enjoy the fascination of particle physics!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Theory background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle physics is the field of physics where structures of matter and radiation and interactions between them are studied. In experimental particle physics research is made by accelerating particles and colliding them to others or to solid targets. This is done with the _particle accelerators_. The collisions are examined with _particle detectors_.\n",
    "\n",
    "World's biggest particle accelerator, Large Hadron Collider (LHC), is located at CERN, the European Organization for Nuclear Research. LHC is 27 kilometers long circle-shaped synchrotron accelerator. LHC is located in the tunnel 100 meters underground on the border of France and Switzerland (image 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/LHC.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "    <figcaption>Image 1: The LHC accelerator and the four detectors around it. &copy; <a href=\"https://cds.cern.ch/record/1708847\">CERN</a> [1]</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "In 2012 the ATLAS and CMS experiments at CERN made an announcement that they had observed a new particle which mass was equal to the predicted mass of the Higgs boson. The Higgs boson and the Higgs field related to it explain the origin of the mass of particles. In 2013 Peter Higgs and François Englert, who predicted the Higgs boson theoretically, were awarded with the Nobel prize in physics."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Accelerating particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LHC mainly accelerates protons. The proton source of the LHC is a bottle of hydrogen. Protons are produced by stripping the electrons away from the hydrogen atoms with help of an electric field.\n",
    "\n",
    "Accelerating process starts already before the LHC. Before the protons arrive in the LHC they will be accelerated with electric fields and directed with magnetic fields in Linac 2, Proton Synchrotron Booster, Proton Synchrotron and Super Proton Synchrotron accelerators. After those the protons will receive energy of 450 GeV. Also the protons will be directed into constantly spreaded bunches in two different proton beams. Each beam contains 2808 proton bunches located about 7,5 meters from each others. Each of these bunches include $1\\text{,}2\\cdot 10^{11}$ protons.\n",
    "\n",
    "After the pre-accelerating the two proton beams are directed to the LHC accelerator. The beams will circulate in opposite directions in two different vacuum tubes. Image 2 shows a part of the LHC accelerator opened with the two vacuum tubes inside. Each of the proton beams will reach the energy of about 7 TeV (7000 GeV) in LHC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/LHC2.jpg\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "    <figcaption>Image 2: Part of the LHC accelerator opened. &copy; <a href=\"https://cds.cern.ch/record/905940\">CERN</a> [2]</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle collisions are created by crossing these two beams that are heading in opposite directions. When two proton bunches cross not all of the protons collide with each others. Only about 40 protons per bunch will collide and so create about 20 collisions. But because the bunches are travelling so fast, there will be about 40 million bunch crosses per one second in the LHC. That means there will be 800 million proton collisions every second in the LHC. That's a lot of action!\n",
    "\n",
    "The maximum energy in collisions is 14 TeV. However in most cases the collision energy is smaller than that because when protons collide it is really the quarks and gluons which collide with each others. So all of the energy of the protons won't be transmitted to the collision.\n",
    "\n",
    "When the protons collide the collision energy can be transformed into mass ($E=mc^2$). So it is possible that new particles are produced in the collisions. By examining and measuring the particles created in collisions, researchers try to understand better for example the dark matter, antimatter and the constitution of all matter.\n",
    "\n",
    "In image 3 there is a visualisation of some particles created in one collision event. These particles are detected with the CMS detector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/eventdisplay.png\" alt=\"image missing\" style=\"height: 450px\" />\n",
    "    <figcaption>Image 3: A visualised collision event.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Video"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The acceleration and collision processes are summarised well in the short video below. Watch the video from the start until 1:15 to get a picture about these processes. You can start the video by running the code cell below (click the cell and then press CTRL + ENTER)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/pQhbhpU9Wrg\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.display import HTML\n",
    "HTML('<iframe width=\"560\" height=\"315\" src=\"https://www.youtube.com/embed/pQhbhpU9Wrg\" frameborder=\"0\" allow=\"accelerometer; autoplay; encrypted-media; gyroscope; picture-in-picture\" allowfullscreen></iframe>')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examining particle collisions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle collisions are examined with _particle detectors_. In LHC particle beams are crossed in four different sections. These sections are the locations of four particle detectors of LHC: ATLAS, LHCb, ALICE and CMS (check the image 1). This exercise focuses on the CMS detector and on the data it collects.\n",
    "\n",
    "CMS, the Compact Muon Solenoid, is a general-purpose detector. Goals of the CMS are for example studying the standard model, searching for extra dimensions and searching for particles that could make up dark matter.\n",
    "\n",
    "Simplified, the goal of the CMS detector is to detect particles that are created in collisions and measure different quantities from them. The detector consists of different detectors that can detect and measure different particles. The structure of the CMS detector is shown in the image 4."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMS.jpg\" alt=\"image missing\" style=\"height: 360px\" />\n",
    "    <figcaption>Image 4: The CMS detector opened. &copy; <a href=\"https://cds.cern.ch/record/1433717\">CERN</a> [3]</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The detectors form an onion-like structure to the CMS. This structure ensures that as many as possible particles from the collision is detected.\n",
    "\n",
    "Different particles act differently in the detectors of the CMS. Image 5 shows the cross-section of the CMS detector. The particle beams would travel in and out from the plane. Image 5 also demonstrates how different particles act in the CMS."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMS2.gif\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "    <figcaption>Image 5: The cross-section of the CMS and different particle interactions in it. &copy; <a href=\"https://cms-docdb.cern.ch/cgi-bin/PublicDocDB/ShowDocument?docid=4172\">CERN</a> [4]</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Innermost part is the silicon tracker. The silicon tracker makes it possible to reconstruct trajectories of charged particles. Charged particles interact electromagnetically with the tracker and make the tracker to create an electric pulse. An intense magnetic field bends the trajectories of the charged particles. With the curvature of the trajectories shown by the pulses created in the tracker, it is possible to calculate the momenta of the charged particles.\n",
    "\n",
    "Particle energies can be measured with help of the calorimeters. Electrons and photons will stop to the Electromagnetic Calorimeter (ECAL). Hadrons, for example protons or neutrons, will pass through the ECAL but will be stopped in the Hadron Calorimeter (HCAL).\n",
    "\n",
    "ECAL is made from lead tungstate crystals that will produce light when electrons and photons pass through them. The amount of light produced is propotional to the energy of the particle. So it is possible to determine the energy of the particle stopped in ECAL with the photodetectors. Also the operation of the HCAL is based on detecting light.\n",
    "\n",
    "Only muons and weakly interacting particles like neutrinos will pass both the ECAL and HCAL. Energies and momenta of muons can be determined with the muon chambers. The detection of the momentum is based on electrical pulses that muons create in the different sections of the muon chambers. Energies of muons can't be measured directly, but the energies will be determined by calculating them from the other measured quantities.\n",
    "\n",
    "Neutrinos can't be detected directly with the CMS, but the existence of them can be derived with the help of missing energy. It is possible that the total energy of the particles detected in a collision is smaller than the energy before the collision. This makes a conflict with the energy conservation. The situation indicates that something has been left undetected in the collision, so there is a possibility that neutrons are created in the collision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise focuses on muons that are detected with the CMS detector. How can you describe the behaviour and detection of muons in the CMS?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Recording the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As mentioned above, there happens about billion particle collision in the CMS in one second. The detector can detect all of these but it would be impossible to record all data from these collisions. Instead right after a collision different trigger systems will decide whether the collision has been potentially interesting or not. Non-interesting collision will not be recorded. This multi-staged triggering process reduces the amount of recorded collisions from billion to about thousand collisions per second.\n",
    "\n",
    "Data collected from collisions will be saved to AOD (Analysis Object Data) files that can be opened with the ROOT program (https://root.cern.ch/). Structures of the files are very complicated so those can't be handled for example in simple data tables.\n",
    "\n",
    "In this exercise a CSV file format is used instead of the AOD format. A CSV file is just a regular text file that contains different values separated with commas (check the image 6). These files can be easily read and handled with the Python programming language."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/datafile.png\" alt=\"image missing\" style=\"height: 80px\" />\n",
    "    <figcaption>Image 6: An example of the structure of the CSV file.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indirect detection of particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not every particle can be detected directly as explained above with the CMS or other particle detectors. Interesting processes are often short-lived. These processes can be searched throughout long-lived processes so detecting is then indirect.\n",
    "\n",
    "For example the Z boson (the particle that mediates weak interaction) can't be detected directly with the CMS since the lifetime of the Z is very short. That means that the Z boson will decay before it even reaches the silicon detector of the CMS.\n",
    "\n",
    "How it is possible to detect the Z boson then? A solution to this question comes from the decay process of the Z boson. If particles that originate from the decay of the Z are prossible to detect, it is also possible to deduce the existence of the Z. So the detection is indirect.\n",
    "\n",
    "The Z boson can decay with 24 different ways. In this exercise only one of these is observed: the decay of the Z to the muon $\\mu^-$ and the antimuon $\\mu^+$. This decay process is shown as a Feynman diagram in the image 7."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/Zdecay.png\" alt=\"image missing\" style=\"height: 170px\" />\n",
    "    <figcaption>Image 7: The process where the Z boson decays to the muon and the antimuon.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Muons that are created in the decay process can be detected with the CMS. But just the detection of the muon and the antimuon isn't a sufficient evidence of the existence of the Z. The detected two muons could originate from any of processes that will happen in the collision event (there are many different processes going on the same time). Because of this the mass of the Z is also needed to be reconstructed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The invariant mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mass of the Z boson can be determined with the help of a concept called _invariant mass_. Let's next derive loosely an expression for the invariant mass.\n",
    "\n",
    "Let's observe a situation where a particle with mass $M$ and energy $E$ decays to two particles with masses $m_1$ and $m_2$, and energies $E_1$ and $E_2$. Energy $E$ and momentum $\\vec{p}$ is concerved in the decay process so $E = E_1 +E_2$ and $\\vec{p} = \\vec{p}_1+ \\vec{p}_2$.\n",
    "\n",
    "Particles will obey the relativistic dispersion relation:\n",
    "\n",
    "$$\n",
    "Mc^2 = \\sqrt{E^2 - c^2\\vec{p}^2}.\n",
    "$$\n",
    "\n",
    "And with the concervation of energy and momentum this can be shown as\n",
    "\n",
    "$$\n",
    "Mc^2 = \\sqrt{(E_1+E_2)^2 - c^2(\\vec{p_1} + \\vec{p_2})^2}\n",
    "$$\n",
    "$$\n",
    "=\\sqrt{E_1^2+2E_1E_2+E_2^2 -c^2\\vec{p_1}^2-2c^2\\vec{p_1}\\cdot\\vec{p_2}-c^2\\vec{p_2}^2}\n",
    "$$\n",
    "$$\n",
    "=\\sqrt{2E_1E_2 - 2c^2 |\\vec{p_1}||\\vec{p_2}|\\cos(\\theta)+m_1^2c^4+m_2^2c^4}. \\qquad (1)\n",
    "$$\n",
    "\n",
    "The relativistic dispersion relation can be brought to the following format\n",
    "\n",
    "$$\n",
    "M^2c^4 = E^2 - c^2\\vec{p}^2\n",
    "$$\n",
    "$$\n",
    "E = \\sqrt{c^2\\vec{p}^2 + M^2c^4},\n",
    "$$\n",
    "\n",
    "from where by setting $c = 1$ (very common in particle physics) and by assuming masses of the particles very small compared to momenta, it is possible to get the following:\n",
    "\n",
    "$$\n",
    "E = \\sqrt{\\vec{p}^2 + M^2} = |\\vec{p}|\\sqrt{1+\\frac{M^2}{\\vec{p}^2}}\n",
    "\\stackrel{M<<|\\vec{p}|}{\\longrightarrow}|\\vec{p}|.\n",
    "$$\n",
    "\n",
    "By applying the result $E = |\\vec{p}|$ derived above and the setting $c=1$ to the equation (1), it can be reduced to the format\n",
    "\n",
    "$$\n",
    "M=\\sqrt{2E_1E_2(1-\\cos(\\theta))},\n",
    "$$\n",
    "\n",
    "where $\\theta$ is the angle between the momentum vector of the particles. With this equation it is possible to calculate the invariant mass for the particle pair if energies of the particles and the angle $\\theta$ is known.\n",
    "\n",
    "In experimental particle physics the equation for the invariant mass is often in the form\n",
    "\n",
    "$$\n",
    "M = \\sqrt{2p_{T1}p_{T2}( \\cosh(\\eta_1-\\eta_2)-\\cos(\\phi_1-\\phi_2) )}, \\qquad (2)\n",
    "$$\n",
    "\n",
    "where transverse momentum $p_T$ is the component of the momentum of the particle that is perpendicular to the particle beam, $\\eta$ the pseudorapidity and $\\phi$ the azimuth angle. The pseudorapidity is defined with the $\\theta$ with the equation $\\eta = -\\ln(\\tan(\\frac{\\theta}{2}))$. So basically the pseudorapidity describes an angle. Also $\\phi$ is describing an angle.\n",
    "\n",
    "Image 8 expresses $\\theta$, $\\eta$ and $\\phi$ in the CMS detector. The particle beams will travel to the z-direction. Image 8 also shows that because of the determination of $\\eta$ it goes to 0 when $\\theta = 90^{\\circ}$ and to $\\infty$ when $\\theta = 0^{\\circ}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMSangles.png\" alt=\"image missing\" style=\"height: 300px\" />\n",
    "    <figcaption>Image 8: Quantities $\\theta$, $\\eta$ and $\\phi$ in the CMS detector.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reconstructing the Z mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With the invariant mass it is possible to prove the existence of the Z boson. In this exercise only the decay of the Z to two muons shown in the image 7 is handled.\n",
    "\n",
    "This exercise uses data that contains collisions where two muons have been detected (among with many of other particles). It is possible to calculate an invariant mass value for the muon pair in an one collision event with the equation (2). And this can be repeated for a great amount of collision events.\n",
    "\n",
    "If the invariant mass of the muon pair is equal to the mass of the Z boson it can be verified that the muon pair originates from the deacay of the Z. And if the invariant mass of the muon pair gets some other value the muons will originate from some other processes. __So the invariant mass can be used as an evidence about the existence of the Z boson__."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Identifying the Z boson"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice the identification of the Z boson goes in the following way. The invariant mass for two muons is calculaetd for the great amount of collision events. Then a histogram is made from the calcuated values. The histogram shows how many invariant mass values will be in each bin of the histogram.\n",
    "\n",
    "If a peak (many invariant mass values near the same bin compared to other bins) is formed in the histogram, it can prove that in the collision events there has been a particle which mass corresponds to the peak. After that it is possible to fit a function to the histogram and determine the mass and the lifetime of the Z from the parameters of the fitted function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's practice the calculation of the invariant mass with the following task. Let's assume that for one muon pair the following values have been measured or determined:\n",
    "\n",
    "- $p_{T1} = 58,6914$ GeV/c\n",
    "- $p_{T2} = 45,7231$ GeV/c\n",
    "- $\\eta_1 = -1,02101$\n",
    "- $\\eta_2 = -0,37030$\n",
    "- $\\phi_1 = 0,836256$ rad\n",
    "- $\\phi_2 = 2,741820$ rad\n",
    "\n",
    "Calculate the invariant mass value for this single pair of muons.\n",
    "\n",
    "Compare the calculated value to the mass of the Z boson reported by the Particle Data Group (PDG, http://pdg.lbl.gov/). What do you notice? Can you make sure conclusions from your notifications?\n",
    "\n",
    "That's the end of the theory part of this exercise. You can now move on to analysing the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Calculating the invariant mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section the data-analysis is started by calculating the invariant masses of the muon pairs that are detected in the collision events. Analysis will be done with the Python programming language.\n",
    "\n",
    "The data used in the analysis has been collected by the CMS detector in 2011. From the original data a CSV file containing only some of the collision events and information has been derived. The original data is saved in AOD format that can be read with ROOT program. Open the link http://opendata.cern.ch/record/17 and take a look how large the original datafile is from the section _Characteristics_.\n",
    "\n",
    "From the original datafile only collision events with exactly two muons detected have been selected to the CSV file. The selection is done with the code similar to the one in the link http://opendata.cern.ch/record/552. In practice the code will select wanted values from the original file and write them to the CSV file. You can get an example of a CSV file by clicking the link http://opendata.cern.ch/record/545 and downloading one of the CSV files from the bottom of the page to your computer.\n",
    "\n",
    "The CSV file used in this excercise is already saved to the same repository than this notebook file. Now let's get the file with Python and start the analysis!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initialisation and getting the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code cell below needed Python modules _pandas_, _numpy_ and _matplotlib.pyplot_ are imported and named as _pd_, _np_ and _plt_. Modules are files that contain functions and commands for Python language. Modules are imported because not all of the things needed in the exercise could be done with the Python's built-in functions.\n",
    "\n",
    "Also the data file from the repository is imported and saved to the variable named `ds`. __Don't change the name of the variable.__ The file is imported with the function `read_csv()` from the pandas module. So in the code there has to be an reference to pandas module (that we named as _pd_) in front of the function.\n",
    "\n",
    "First we want to figure out how many collision events (or in this case data rows) there are in the data file. Add to the code cell below needed code to print out the number of rows of the imported file. With Python printing is done with the `print()` function where the thing that is wanted to be printed will be written inside the brackets. The length of an object can be determined with the `len()` function. Inside the brackets will be written the variable which length is wanted to be determined.\n",
    "\n",
    "You can run the code cell by clicking it active and then pressing CTRL + ENTER. Feel free to test different solutions for printing the length of the file.\n",
    "\n",
    "After you have printed the number of the rows in the datafile, you can move on to the next section. First try to figure it out yourself, but if you get stuck click on the hints below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    The data was saved to the variable that was named as \"ds\".\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    Write the function \"len()\" inside the function \"print()\": \"print(len(variablename))\", where variablename\n",
    "    refers to the name of your variable.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "475465\n",
      "\n",
      "      Run     Event Type1        E1       px1       py1       pz1      pt1  \\\n",
      "0  165617  74601703     G   9.69873 -9.510430  0.366205   1.86329  9.51748   \n",
      "1  165617  75100943     G   6.20385 -4.266610  0.456545  -4.47930  4.29097   \n",
      "2  165617  75587682     G  19.28920 -4.212080 -0.651623  18.81210  4.26219   \n",
      "3  165617  75660978     G   7.04268 -6.326780 -0.268545   3.08025  6.33248   \n",
      "4  165617  75947690     G   7.27507  0.102966 -5.533080  -4.72122  5.53403   \n",
      "\n",
      "       eta1     phi1  ...  Type2       E2       px2       py2      pz2  \\\n",
      "0  0.194546  3.10311  ...      G  9.76332  7.327710 -1.152370  6.34728   \n",
      "1 -0.912070  3.03499  ...      G  9.66898  7.273970 -2.821120 -5.71040   \n",
      "2  2.190460 -2.98811  ...      G  9.82439  4.343940 -0.473505  8.79849   \n",
      "3  0.469033 -3.09917  ...      G  5.58571  4.474760  0.848935 -3.23185   \n",
      "4 -0.773616 -1.55219  ...      G  7.31811 -0.398831  6.940810  2.28247   \n",
      "\n",
      "       pt2      eta2      phi2  Q2         M  \n",
      "0  7.41776  0.775564 -0.155984   1  17.49220  \n",
      "1  7.80188 -0.678625 -0.369977   1  11.55340  \n",
      "2  4.36967  1.449670 -0.108575   1   9.16361  \n",
      "3  4.55458 -0.660499  0.187488   1  12.47740  \n",
      "4  6.95226  0.322677  1.628190   1  14.31590  \n",
      "\n",
      "[5 rows x 21 columns]\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "ds = pd.read_csv('DoubleMuRun2011A.csv')\n",
    "\n",
    "# This is a comment separated with #-symbol. Comments do not affect to the code.\n",
    "# Add your own code to print the number of collision events in the datafile!\n",
    "print(len(ds))\n",
    "print(\"\")\n",
    "print(ds[0:5])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What does the file look like?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file was saved as a _DataFrame_ structure (practically a table) of _pandas_ module in a variable called `ds`. Next print the five first rows of the file to look properly how does the file look. With the `print()` function it is possible to print a variable inside the brackets. With the function _variablename_`.head()` you can get the first five rows of the data file by changing the _variablename_ with the name of your variable.\n",
    "\n",
    "Write a code that prints the five first rows of the data file and run the code cell by clicking it active and pressing CTRL + ENTER. First try to figure it out yourself, but if you get stuck click on the hint below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint</summary>\n",
    "    \n",
    "    Hint: \"print(variablename.head())\"\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      Run     Event Type1        E1       px1       py1       pz1      pt1  \\\n",
      "0  165617  74601703     G   9.69873 -9.510430  0.366205   1.86329  9.51748   \n",
      "1  165617  75100943     G   6.20385 -4.266610  0.456545  -4.47930  4.29097   \n",
      "2  165617  75587682     G  19.28920 -4.212080 -0.651623  18.81210  4.26219   \n",
      "3  165617  75660978     G   7.04268 -6.326780 -0.268545   3.08025  6.33248   \n",
      "4  165617  75947690     G   7.27507  0.102966 -5.533080  -4.72122  5.53403   \n",
      "\n",
      "       eta1     phi1  ...  Type2       E2       px2       py2      pz2  \\\n",
      "0  0.194546  3.10311  ...      G  9.76332  7.327710 -1.152370  6.34728   \n",
      "1 -0.912070  3.03499  ...      G  9.66898  7.273970 -2.821120 -5.71040   \n",
      "2  2.190460 -2.98811  ...      G  9.82439  4.343940 -0.473505  8.79849   \n",
      "3  0.469033 -3.09917  ...      G  5.58571  4.474760  0.848935 -3.23185   \n",
      "4 -0.773616 -1.55219  ...      G  7.31811 -0.398831  6.940810  2.28247   \n",
      "\n",
      "       pt2      eta2      phi2  Q2         M  \n",
      "0  7.41776  0.775564 -0.155984   1  17.49220  \n",
      "1  7.80188 -0.678625 -0.369977   1  11.55340  \n",
      "2  4.36967  1.449670 -0.108575   1   9.16361  \n",
      "3  4.55458 -0.660499  0.187488   1  12.47740  \n",
      "4  6.95226  0.322677  1.628190   1  14.31590  \n",
      "\n",
      "[5 rows x 21 columns]\n"
     ]
    }
   ],
   "source": [
    "print(ds.head())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The \"\\\\\" symbols in the output tell that a row won't fit totally on a screen but continues to next rows of the output. The first row shows which information about muon pairs the file contains. For example E1 is the energy of the first muon and E2 the energy of the second etc. Here are the different values listed:\n",
    "\n",
    "- Run = number of the run where data has been collected from\n",
    "- Event = number of the collision event\n",
    "- Type = type of the muon, global muon (G) has been measured both in the silicon tracker and muon chambers, tracker myon (T) has been measured only in the silicon tracker (these classifications are hypotheses since the type cannot be known absolutely)\n",
    "- E = energy of the muon\n",
    "- px, py, pz = different coordinates of momentum of the muon\n",
    "- pt = transverse momentum, that is the component of momentum of the muon that is perpendicular to the particle beams\n",
    "- eta = $\\eta$ = pseudorapidity, a coordinate describing an angle (check the image 8)\n",
    "- phi = $\\phi$ = azimuth angle, also a coordinate describing an angle (check the image 8)\n",
    "- Q = electrical charge of the muon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Calculating the invariant mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next calculate invariant mass values for muon pairs in each event with the different values from the data file. You have to write a proper equation only once since code executes the equation automatically for each row of the file.\n",
    "\n",
    "For example if you would like to sum the electrical charges of two muons for each event and save results in a variable _charges_, it could be done with the following code:\n",
    "```\n",
    "charges = ds.Q1 + ds.Q2\n",
    "```\n",
    "\n",
    "So you have to tell in the code that Q1 and Q2 refer to values in the variable `ds`. This can be done by adding the variable name separated with a dot in front of the value that is wanted, as in the example above.\n",
    "\n",
    "There are square root, cosine and hyperbolic cosine terms in the equation of invariant mass. Those can be fetched from the _numpy_ module that we named as _np_. You can get a square root with the function `np.sqrt()`, a cosine with `np.cos()` and a hyperbolic cosine with `np.cosh()`. Naturally inside the brackets there will be anything that is inside the square root or brackets in the equation too.\n",
    "\n",
    "__Write below a code__ that will calculate the invariant mass value for muon pairs in each collision event in the data file. Save the values calculated in the variable `invariant_mass` that is already written in the code cell. Don't change the name of the variable.\n",
    "\n",
    "After running, the code will print the first five values that are calculated. Also the output will tell if the values are correct. This is done with a small piece of code at the end of the cell.\n",
    "\n",
    "You can get help from the theory part. Also use the hints below if you get stuck. But first try different solutions by yourself and try to figure it out without the hints!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    Use the equation (2) of the theory part for the calculation.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    When you write different quantities of the equation to your code, remember to refer to the variable from where you want to get the quantities. For example if you would need the quantity \"pt1\", write \"ds.pt1\" to the code.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 3</summary>\n",
    "    \n",
    "    In practice write the equation (2) to one line to the code after the text \"invariant_mass = \". Remember that you can get a cosine, a hyperbolic cosine and a square root from \"numpy\" module with the way that is described above. Also remember to tell from which variable you want to get the different quantities (hint 2).\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The first five values calculated (in units GeV):\n",
      "0    17.490875\n",
      "1    11.551382\n",
      "2     9.160909\n",
      "3    12.475632\n",
      "4    14.314303\n",
      "dtype: float64\n",
      "Invariant mass values are correct!\n"
     ]
    }
   ],
   "source": [
    "invariant_mass = np.sqrt(2*ds.pt1*ds.pt2*((np.cosh(ds.eta1-ds.eta2))-(np.cos(ds.phi1-ds.phi2))))\n",
    "print('The first five values calculated (in units GeV):')\n",
    "print(invariant_mass[0:5])\n",
    "\n",
    "# Rest of the code is for checking if the values are correct. You don't have to change that.\n",
    "if 14.31 <= invariant_mass.values[4] <= 14.32:\n",
    "    print('Invariant mass values are correct!')\n",
    "else:\n",
    "    print('Calculated values are not yet correct. Please check the calculation one more time.')\n",
    "    print('Remember: don´t change the name of the variable invariant_mass.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Making the histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's make a histogram from the calculated invariant mass values. The histogram describes how the values are distributed, that is, how many values there has been in each bin of the histogram. In the image 9 there is a histogram that represents how the amount of cash in a wallet has been distributed for some random group of people. One can see from the histogram that for example the most common amount of cash has been 10–15 euros (12 persons have had this)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/histogram.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "    <figcaption>Image 9: An example histogram from the distribution of the amount of cash.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Creating the histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histograms can be created with Python with the _matplotlib.pyplot_ module that was imported before and named as _plt_. With the function `plt.hist()` it is possible to create a histogram by giving different parameters inside the brackets. These parameters can be examined from https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.hist.html.\n",
    "\n",
    "Now only the first three of the parameters are needed: a variable from which values the histogram is created (_x)_, number of bins (_bins_) and the lower and upper range of the bins (_range_).\n",
    "\n",
    "Write down a code that will create a histogram from the invariant mass values that were calculated. Because this exercise focuses on the Z boson, set the range wisely to get the values near the mass of the Z boson. Use the Z boson mass value that you looked earlier from the Particle Data Group as a reference.\n",
    "\n",
    "Try what is the best amount of bins to make a clear histogram. You can try different values and see how they affect to the histogram.\n",
    "\n",
    "In the code there are already lines for naming the axes and the title of the histogram. Also there are comments marked with # symbols. These comments doesn't affect to the functionality of the code.\n",
    "\n",
    "If you get stuck use the hints below. But try to create the histogram without using the hints!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    The invariant mass values that you have calculated are saved in the variable \"invariant_mass\".\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    The function is in the form \"plt.hist(x, bins=0, range=(0,0))\", where x will be replaced with the name of the variable that contains the data that is wanted to be used in the histogram (in our case the invariant masses). The zeroes will be replaced with the wanted amount of bins and with the lower and upper limits of the histogram.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 3</summary>\n",
    "    \n",
    "    Try different bin values between 50 and 200.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Valor mínimo = 0.05904545300204003\n",
      "Valor máximo = 386.91839095367584\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZUAAAElCAYAAAAskX9OAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3X+8HFV9//HXmwQk/AjhR6CQgAGMPwA14BViVaSgEEAJtlBBlEiRqIUWKyqBKhGRr9AqVqpSUVLAKiFFhQhoDMiPovxIIoEQQImAEBNJMAmE3wKf7x/nrEw2u/fOvZm9e/fm/Xw89rEzZ87MnNnZnc+cM7NnFBGYmZlVYYN2F8DMzAYPBxUzM6uMg4qZmVXGQcXMzCrjoGJmZpVxUDEzs8o4qFRI0kJJ+7W7HO0k6f2SHpX0lKQ9G0x/StIubSjXMZJ+3t/r7XSSbpT00XaXo0bS6yTdKWm1pH9ud3lsbQ4qJUl6WNK769I+IumW2nhE7B4RN/awnDGSQtLQFhW13b4CnBQRm0XEnfUTc/qD/V2oiPh+RBxYxbIafRes33wWuDEiNo+I8+snDrQguD5yUBlkBkCwejWwsM1lWMMA+EysOgPu+2V1IsKvEi/gYeDddWkfAW5plAfYG5gLPAk8BpyX0x8BAngqv95GCu6fA34PLAMuBbYoLPfYPO1PwOfr1vMF4Argf/K6PprXfSuwClgKfAPYqLC8AP4ReABYDZwF7JrneRKYUcxft80Nywq8Km9PAE8Dv2syfwCvycMXA98ErsnluB3YNU/7L+ArdfNeBXwqD08Bfpfnuxd4f91++SXwNWAF8KUG++rrwKN5e+cB7yxM+0L+DC7Ny18IdOVp3wNeBp7N2/vZBtu4H7CYdFa9LO+Dw4FDgN/mMp1eyN90fwHK27EMeAK4G9gjTzskb/tq4A/ApxuU5VV5uXsU0kbm8m8LbAlcDSwHVubh0YW8NwIfLXwu/1OYNibvz6F5fAvgorwNf8if+5A87TXATXkbHgcu7+a3dlj+zFfl9b8hp/8CeAl4Ln/2r62b7+y66d8AzgT+M0/fkPTd/Lc8Pizn3bK79XbzPS71G6Luu9fgd7AF6bu2nPS7+hywQXFeUgvASuAh4OC67/qDuQwPAce0/VjZ7gJ0yoveB5VbgQ/n4c2A8Xl4jR9iTvsHYBGwS877I+B7edpu+QfyDmCj/OX6M2sGlT+TDlob5B/KW4DxwNC8vvuAT9Z9oWcCw4HdgeeB6/P6tyAdqCY1+RyalrWw7Nd08znWB5UVpIPqUOD7wPQ8bV/SQV95fEvSgXCHPH4ksEPe5g+QDhbbF/bLi8A/5eUOa7CvPgRsnaefAvwR2LjwmT5HOmgPAb4M3Nbdd6FuG/fL6z+DdCA7gXTA+AGwef7MnwN2yfmb7i/gIFLQG0EKMG8obOdScjDMn89eTcozDTi7MH4i8LM8vDXwd8AmuWz/C1xZyHsj5YPKlcC3gU1JAesO4GN52mXAv+b9tTHwjiZlfW3el+/Jn91nSd+3jerL02T+NaYD+wML8vBfk05Ebi9Mu6vMept8j0v9hug5qFxKOmHaPH+mvwWOL8z7Z9J3aAjwCWBJ/i5sSgpgr8t5twd2b/uxst0F6JQX6UDyFOkspvZ6huZB5WbSWdI2dctZ44eY064H/rEw/rr8RRpKOjBdVpi2CfACawaVm3so+yeBHxfGA3h7YXwecGph/KvAfzRZVtOyFpbdm6Dy3cK0Q4D787BItbp98/gJwC+6We58YGIe/gjwSN30tX7YddNXAm8ufKbXFabtBjzbaD83WdZ+pABYO0vfPG/3PnWf+eE97S/Sge+3pKCzQV2+R4CPAcN72P/vBh4sjP8SOLZJ3nHAysL4jZQIKsB2pAPrsML0o4Eb8vClwIUUakFN1v95YEZhfANSrWe/+vI0mX+N6bxSG9maVLs9nVSL3Iz0+zy/zHqbfI9L/YYafffy/K8hBYrngd0K0z5Gum5Um3dRYdomed6/IgWVVaSTgmHNPpP+fvmaSu8cHhEjai9S9beZ40lnP/dLmiPpvd3k3YFU7a35Pa/8UHcgnbEDEBHPkJrBih4tjkh6raSrJf1R0pPA/wO2qZvnscLwsw3GN+tDWfvij4XhZ2rrjfQLmk46MAF8kFSTAUDSsZLmS1olaRWwB2tu4xqfST1Jp0i6T9ITef4t6uavL9fGvbw286eIeCkPP5vfG37G3e2viPgFqRnnm8Bjki6UNDwv4+9Igfj3km6S9LYmZfkFMEzSPpJeTQocP87r3kTStyX9Pq/7ZmCEpCG92FZI1zo2BJYW9sm3STUWSGf+Au7Id0n+Q5PlrPH9ioiXSftyVC/LU5v/WVIz9LtItd+bgF8Bb89pN63Devv6GyrahtQCUf+bKq73L9/F/PsH2CwinibV0j9O+tyvkfT6EutsKQeVFomIByLiaNKP6lzgCkmbks4y6i0h/ShrdiI1nzxGauIYXZsgaRjprGuN1dWNXwDcD4yNiOGkszP1fWtKl7VqlwFH5APhPsAPAfL4d4CTgK1zgL+HNbex0edMnv+dwKnA35Pa00eQ2vrLfkZNl91H3e6viDg/It5CamZ5LfCZnD4nIiaSvmNXktrx1y5sOkDOIAXoDwJXR8TqPPkUUm1zn7zufXN6o8/iadKZcs1fFYYfJZ1xb1M48RoeEbvnMvwxIk6IiB1IZ+LfkvSaButY4/slScCOpFpDGY32zU2kGt+ewJw8fhCp2fXmitbbnTU+N0nFz+1xUk2//jdVar0RMSsi3kNq+rqf9LtoKweVFpH0IUkj8w96VU5+idS2/jKp7bXmMuBfJO0saTPSmerlEfEi6SL8+yT9taSNSFX2ng5+m5PaWp/KZy6fqGzDui9rpSLdkrwc+C4wKyJqn2MtOC8HkHQcqaZS1uakQLgcGCrpDFLbeFmPseb+W1dN95ekt+YaRu0i83PAS5I2yv+92SIi/pznf6nRwrMfkM5qj8nDxXU/C6yStBUwtZtlzAf2lbSTpC2A02oTImIp8HPgq5KGS9pA0q6S3pW340hJtZOjlaT916i8M4BDJR2Qt/kUUrD6VTflKmq0b24i3exyb0S8QG4iAx6KiOUVrbc7dwG7SxonaWNSMyIAuTY7Azhb0ub5hOlTpBtvuiVpO0mH5ZPV50nN8919B/qFg0rrTAAWSnqKdKfRURHxXK6+ng38MjcTjCddSP0e6azpIdKB458AImJhHp5OqrWsJt0J9Hw36/406Yx0NenM5fIKt6tpWVvkMtI1gb8cCCPiXlKb9a2kg8gbSdcJypoF/JR0reL3pG3otrmszpeBz+X99+lezNdMd/treE5bySt3AH4lT/sw8HButvo46eaDhiLidlJQ2oG07TX/Qbru8DhwG/CzbpYxO5ftbtI1hKvrshxLasq5N5f3CtIZNMBbgdvz72EmcHJEPNRgHb/J2/GfuUzvA96Xg0EZXyfVbldKqv2P5Vd5G2u1kntJ+7w2XsV6m4qI3wJfBK4j3S12S12WfyLtmwfztB+Qfmc92YAU/JaQbnh5F7lJXtI782fd72p31liHyLWDVaSmkrV+lGZm7eSaSgeQ9L58QXVT0lnqAtIdSGZmA4qDSmeYSKriLgHGkprSXMU0swHHzV9mZlYZ11TMzKwyDipmZlYZBxUzM6uMg4qZmVXGQcXMzCrjoGJmZpVxUDEzs8o4qJiZWWUcVMzMrDIOKmZmVhkHFTMzq4yDipmZVcZBxczMKuOgYmZmlXFQMTOzyjiomJlZZYa2uwD9bZtttokxY8a0uxhmZh1l3rx5j0fEyJ7yrXdBZcyYMcydO7fdxTAz6yiSfl8mn5u/zMysMg4qZmZWGQcVMzOrjIOKmZlVxkHFzMwq46BiZmaVcVAxM7PKOKiYmVllHFTMzKwyDioVGjPlmnYXwcysrRxUzMysMg4qZmZWGQcVMzOrjIOKmZlVxkHFzMwq46BiZmaVcVAxM7PKOKiYmVllHFTMzKwyLQsqkjaWdIekuyQtlHRmTr9Y0kOS5ufXuJwuSedLWiTpbkl7FZY1SdID+TWpkP4WSQvyPOdLUqu2x8zMeja0hct+Htg/Ip6StCFwi6Sf5mmfiYgr6vIfDIzNr32AC4B9JG0FTAW6gADmSZoZEStznsnAbcC1wATgp5iZWVu0rKYSyVN5dMP8im5mmQhcmue7DRghaXvgIGB2RKzIgWQ2MCFPGx4Rt0ZEAJcCh7dqe8zMrGctvaYiaYik+cAyUmC4PU86OzdxfU3Sq3LaKODRwuyLc1p36YsbpDcqx2RJcyXNXb58+Tpvl5mZNdbSoBIRL0XEOGA0sLekPYDTgNcDbwW2Ak7N2RtdD4k+pDcqx4UR0RURXSNHjuzlVpiZWVn9cvdXRKwCbgQmRMTS3MT1PPDfwN4522Jgx8Jso4ElPaSPbpBuZmZt0sq7v0ZKGpGHhwHvBu7P10LId2odDtyTZ5kJHJvvAhsPPBERS4FZwIGStpS0JXAgMCtPWy1pfF7WscBVrdoeMzPrWSvv/toeuETSEFLwmhERV0v6haSRpOar+cDHc/5rgUOARcAzwHEAEbFC0lnAnJzvixGxIg9/ArgYGEa668t3fpmZtVHLgkpE3A3s2SB9/yb5AzixybRpwLQG6XOBPdatpGZmVhX/o97MzCrjoGJmZpVxUDEzs8o4qJiZWWUcVMzMrDIOKmZmVhkHFTMzq4yDipmZVcZBxczMKuOgYmZmlXFQMTOzyjiomJlZZRxUzMysMg4qZmZWGQcVMzOrjIOKmZlVxkHFzMwq46BiZmaVaVlQkbSxpDsk3SVpoaQzc/rOkm6X9ICkyyVtlNNflccX5eljCss6Laf/RtJBhfQJOW2RpCmt2hYzMyunlTWV54H9I+LNwDhggqTxwLnA1yJiLLASOD7nPx5YGRGvAb6W8yFpN+AoYHdgAvAtSUMkDQG+CRwM7AYcnfOamVmbtCyoRPJUHt0wvwLYH7gip18CHJ6HJ+Zx8vQDJCmnT4+I5yPiIWARsHd+LYqIByPiBWB6zmtmZm3S0msquUYxH1gGzAZ+B6yKiBdzlsXAqDw8CngUIE9/Ati6mF43T7N0MzNrk5YGlYh4KSLGAaNJNYs3NMqW39VkWm/T1yJpsqS5kuYuX76854KbmVmf9MvdXxGxCrgRGA+MkDQ0TxoNLMnDi4EdAfL0LYAVxfS6eZqlN1r/hRHRFRFdI0eOrGKTzMysgVbe/TVS0og8PAx4N3AfcANwRM42CbgqD8/M4+Tpv4iIyOlH5bvDdgbGAncAc4Cx+W6yjUgX82e2anvMzKxnQ3vO0mfbA5fku7Q2AGZExNWS7gWmS/oScCdwUc5/EfA9SYtINZSjACJioaQZwL3Ai8CJEfESgKSTgFnAEGBaRCxs4faYmVkPWhZUIuJuYM8G6Q+Srq/Upz8HHNlkWWcDZzdIvxa4dp0La2ZmlfA/6s3MrDIOKmZmVhkHFTMzq4yDipmZVcZBxczMKuOgYmZmlXFQMTOzyjiomJlZZRxUzMysMg4qZmZWmR6DiqR/kzRc0oaSrpf0uKQP9UfhzMyss5SpqRwYEU8C7yV1N/9a4DMtLZWZmXWkMkFlw/x+CHBZRKxoYXnMzKyDleml+CeS7geeBf5R0kjgudYWy8zMOlGZmspU4G1AV0T8GXgGOKylpTIzs45UJqjcGhEraw/GioingZ+2tlhmZtaJmjZ/SforYBQwTNKegPKk4cAm/VA2MzPrMN1dUzkI+AgwGjivkL4aOL2FZTIzsw7VNKhExCWkZ8z/XUT8sB/LZGZmHarMNZWrJX1Q0umSzqi9eppJ0o6SbpB0n6SFkk7O6V+Q9AdJ8/PrkMI8p0laJOk3kg4qpE/IaYskTSmk7yzpdkkPSLpc0ka93H4zM6tQmaByFTAReBF4uvDqyYvAKRHxBmA8cKKk3fK0r0XEuPy6FiBPOwrYHZgAfEvSEElDgG8CBwO7AUcXlnNuXtZYYCVwfIlymZlZi5T5n8roiJjQ2wVHxFJgaR5eLek+0oX/ZiYC0yPieeAhSYuAvfO0RRHxIICk6cDEvLz9gQ/mPJcAXwAu6G1ZzcysGmVqKr+S9MZ1WYmkMcCewO056SRJd0uaJmnLnDYKeLQw2+Kc1ix9a2BVRLxYl95o/ZMlzZU0d/ny5euyKWZm1o0yQeUdwLx8TeNuSQsk3V12BZI2A34IfDL3IXYBsCswjlST+Wota4PZow/paydGXBgRXRHRNXLkyLJFNzOzXirT/HVwXxcuaUNSQPl+RPwIICIeK0z/DnB1Hl0M7FiYfTSwJA83Sn8cGCFpaK6tFPObmVkb9FhTiYjfkw7q++fhZ8rMJ0nARcB9EXFeIX37Qrb3A/fk4ZnAUZJeJWlnYCxwBzAHGJvv9NqIdDF/ZkQEcANwRJ5/EummAjMza5MeayqSpgJdwOuA/yb1Wvw/wNt7mPXtwIeBBZLm57TTSXdvjSM1VT0MfAwgIhZKmgHcS7pz7MRa1zCSTgJmAUOAaRGxMC/vVGC6pC8Bd5KCmJmZtUmZ5q/3ky6y/xogIpZI2rynmSLiFhpf97i2m3nOBs5ukH5to/nyHWF716ebmVl7lLlQ/0JuagoASZu2tkhmZtapygSVGZK+TboofgJwHfCd1hbLzMw6UY/NXxHxFUnvAZ4kXVc5IyJmt7xkZmbWccpcqP8X4H8dSMzMrCdlmr+GA7Mk/Z+kEyVt1+pCmZlZZyrzP5UzI2J34ERgB+AmSde1vGRmZtZxytRUapYBfwT+BGzbmuKYmVknK/PP+E9IuhG4HtgGOCEi3tTqgpmZWecp8+fHV5M6g5zfY04zM1uvlbmmMgXYTNJxAJJG5r65zMzM1lCm+WsqqY+t03JSre8vMzOzNZS5UP9+4DDyI4QjYgnQY99fZma2/nHfX2ZmVhn3/WVmZpVx319mZlaZMrcUk4OIA4mZmXWrN/+oNzMz65aDipmZVaZpUJF0fX4/t/+KY2Zmnay7msr2kt4FHCZpT0l7FV89LVjSjpJukHSfpIWSTs7pW0maLemB/L5lTpek8yUtknR3cR2SJuX8D0iaVEh/i6QFeZ7zJanvH4WZma2r7i7UnwFMAUYD59VNC2D/Hpb9InBKRPxa0ubAPEmzgY8A10fEOZKm5HWcChwMjM2vfYALgH0kbQVMBbryeudJmhkRK3OeycBtwLXABOCnZTbczMyq1zSoRMQVwBWSPh8RZ/V2wRGxFFiah1dLug8YBUwE9svZLgFuJAWVicCl+Y+Wt0kaIWn7nHd2RKwAyIFpQu45eXhE3JrTLwUOx0HFzKxtyvxP5SxJhwH75qQbI+Lq3qxE0hhgT+B2YLsccIiIpZJqz2YZBTxamG1xTusufXGD9Ebrn0yq0bDTTjv1puhmZtYLZTqU/DJwMnBvfp2c00qRtBnwQ1L3+U92l7VBWvQhfe3EiAsjoisiukaOHNlTkc3MrI/K3FJ8KPCeiJgWEdNI1y0OLbNwSRuSAsr3I+JHOfmx3KxFfl+W0xcDOxZmHw0s6SF9dIN0MzNrk7L/UxlRGN6izAz5TqyLgPsionihfyZQu4NrEnBVIf3YfBfYeOCJ3Ew2CzhQ0pb5TrEDgVl52mpJ4/O6ji0sy8zM2qBMNy1fBu6UdAOpyWlfXnm2SnfeDnwYWCCp9tTI04FzSJ1UHg88AhyZp10LHAIsAp4BjgOIiBWSzgLm5HxfrF20Bz4BXAwMI12g90V6M7M2KnOh/rJ8p9VbSUHl1Ij4Y4n5bqHxdQ+AAxrkD+DEJsuaBkxrkD4X2KOnspiZWf8o26HkUlLzlJmZWVPu+8vMzCrjoGJmZpXpNqhI2kDSPf1VGDMz62zdBpWIeBm4S5L/hm5mZj0qc6F+e2ChpDuAp2uJEXFYy0plZmYdqUxQObPlpTAzs0GhzP9UbpL0amBsRFwnaRNgSOuLZmZmnaZMh5InAFcA385Jo4ArW1koMzPrTGVuKT6R1OXKkwAR8QCwbbdzmJnZeqlMUHk+Il6ojUgaSpMu5s3MbP1WJqjcJOl0YJik9wD/C/yktcUyM7NOVCaoTAGWAwuAj5F6E/5cKwtlZmadqczdXy9LuoT0KOAAfpN7FDYzM1tDj0FF0qHAfwG/I3Vlv7Okj0WEn11iZmZrKPPnx68CfxMRiwAk7Qpcgx+IZWZmdcpcU1lWCyjZg7zyXHkzM7O/aFpTkfS3eXChpGuBGaRrKkfyyqN9zczM/qK75q/3FYYfA96Vh5cDW7asRGZm1rGaBpWIOG5dFixpGvBeUvPZHjntC8AJpMAEcHpEXJunnQYcD7wE/HNEzMrpE4Cvk/ob+25EnJPTdwamA1sBvwY+XPyTppmZ9b8yfX/tLOk8ST+SNLP2KrHsi4EJDdK/FhHj8qsWUHYDjgJ2z/N8S9IQSUOAbwIHA7sBR+e8AOfmZY0FVpICkpmZtVGZu7+uBC4i/Yv+5bILjoibJY0pmX0iMD0ingcekrQI2DtPWxQRDwJImg5MlHQfsD/wwZznEuALwAVly2dmZtUrE1Sei4jzK1znSZKOBeYCp0TESlLPx7cV8izOaQCP1qXvA2wNrIqIFxvkX4ukycBkgJ128kMszcxapcwtxV+XNFXS2yTtVXv1cX0XALsC44ClpP/AQPpTZb3oQ3pDEXFhRHRFRNfIkSN7V2IzMyutTE3ljcCHSc1NteavyOO9EhGP1YYlfQe4Oo8uBnYsZB0NLMnDjdIfB0ZIGpprK8X8ZmbWJmWCyvuBXaq4s0rS9hGxtLDce/LwTOAHks4DdgDGAneQaiRj851efyBdzP9gRISkG4AjSHeATQKuWtfymZnZuikTVO4CRtDLf9FLugzYD9hG0mJgKrCfpHGkms7DpF6PiYiFkmYA9wIvAidGxEt5OScBs0i3FE+LiIV5FacC0yV9CbiTdDOBmZm1UZmgsh1wv6Q5wPO1xIg4rLuZIuLoBslND/wRcTZwdoP0a0nd7denP8grd4iZmdkAUCaoTG15KczMbFAo8zyVm/qjIGZm1vnKPE9lNa/crrsRsCHwdEQMb2XBzMys85SpqWxeHJd0OL6WYWZmDZT58+MaIuJK+vAfFStnzJRr2l0EM7M+K9P89beF0Q2ALrr597qZma2/ytz9VXyuyouk/5dMbElpzMyso5W5prJOz1UxM7P1R3ePEz6jm/kiIs5qQXk6wpgp1/DwOYe2uxhmZgNOdzWVpxukbUp6GNbWwHobVMzMrLHuHidc65YeSZsDJwPHkTpw/Gqz+czMbP3V7TUVSVsBnwKOIT1dca/8UC0zM7O1dHdN5d+BvwUuBN4YEU/1W6nMzKwjdffnx1NIzzb5HLBE0pP5tVrSk/1TPDMz6yTdXVPp9b/tzcxs/ebAYWZmlXFQMTOzyjiomJlZZVoWVCRNk7RM0j2FtK0kzZb0QH7fMqdL0vmSFkm6W9JehXkm5fwPSJpUSH+LpAV5nvMlqVXbYmZm5bSypnIxMKEubQpwfUSMBa7P4wAHA2PzazJwAfzlfzJTgX1Iz3CZWgtEOc/kwnz16zIzs37WsqASETcDK+qSJ5L+REl+P7yQfmkktwEjJG0PHATMjogV+U+Xs4EJedrwiLg1IgK4tLCsQcHPVTGzTtTf11S2i4ilAPl925w+Cni0kG9xTusufXGD9IYkTZY0V9Lc5cuXr/NGgA/6ZmaNDJQL9Y2uh0Qf0huKiAsjoisiukaOHNnHIpqZWU/6O6g8lpuuyO/LcvpiYMdCvtHAkh7SRzdI71eurZiZram/g8pMoHYH1yTgqkL6sfkusPHAE7l5bBZwoKQt8wX6A4FZedpqSePzXV/HFpbVUvWBxIHFzOwVZR4n3CeSLgP2A7aRtJh0F9c5wAxJxwOPAEfm7NcChwCLgGdIXewTESsknQXMyfm+GBG1i/+fIN1hNgz4aX6ZmVkbtSyoRMTRTSYd0CBvACc2Wc40YFqD9LnAHutSRjMzq9ZAuVBvZmaDgIPKAObrNWbWaRxUzMysMg4qZmZWGQeVCriZyswsadndX9Z3jYJUfdrD5xzKmCnX8PA5h/ZXsczMeuSaSoeqBRnXksxsIHFQ6QAOHGbWKRxUBrgyAcVBx8wGCgeVQWLMlGscXMys7RxUBhkHFjNrJwcVMzOrjINKRVxDMDNzUBmUHODMrF0cVMzMrDIOKoOUaytm1g4OKoOYA4uZ9TcHlUHOgcXM+pODipmZVaYtQUXSw5IWSJovaW5O20rSbEkP5Pctc7oknS9pkaS7Je1VWM6knP8BSZPasS2dwLUVM+sv7ayp/E1EjIuIrjw+Bbg+IsYC1+dxgIOBsfk1GbgAUhACpgL7AHsDU2uByNbmwGJm/WEgNX9NBC7Jw5cAhxfSL43kNmCEpO2Bg4DZEbEiIlYCs4EJ/V3oTuLAYmat1q6gEsDPJc2TNDmnbRcRSwHy+7Y5fRTwaGHexTmtWfpaJE2WNFfS3OXLl1e4GWZmVtSuJz++PSKWSNoWmC3p/m7yqkFadJO+dmLEhcCFAF1dXQ3zmJnZumtLTSUiluT3ZcCPSddEHsvNWuT3ZTn7YmDHwuyjgSXdpLfMYGg+GgzbYGYDV78HFUmbStq8NgwcCNwDzARqd3BNAq7KwzOBY/NdYOOBJ3Lz2CzgQElb5gv0B+a0jtVfB3wHFjNrlXY0f20H/FhSbf0/iIifSZoDzJB0PPAIcGTOfy1wCLAIeAY4DiAiVkg6C5iT830xIlb032Z0tjFTruHhcw5tdzHMbJDp96ASEQ8Cb26Q/ifggAbpAZzYZFnTgGlVl9HMzPpmIN1SbP3MjyA2s6o5qJgDi5lVxkGlYp16gHatxcyq4KAyQAyUA/pAKYeZdSYHFVuLay1m1lcOKi3gA7KZra8cVKwp11jMrLccVKxHDi5mVpaDygDQKQdsBxcz64mDivWag4uZNdOuru9tECgGFvcjZmbgmkrbDZYzftdebLDx97lvlPprXH90dXXF3Llz+zTvunzJamfy69MX1bUX6zSNfp/+HieS5kVEV4/5HFTKW58CQpX8o7RO0N3v29/h8kHF11Ss5XztxQayMieLfv5QeQ4q1q/qf8D+oVqncGApx0HF2spt2NaBUHKWAAAIz0lEQVROvW3SdmDpmYOKDTgONNZqvj7aOg4q1hGaHQQePudQnz1ar6xrQKnN7+9cYx1/95ekCcDXgSHAdyPinO7y++6v9ZcPAuuvVv5215fv1Xpx95ekIcA3gfcAi4E5kmZGxL3tLZkNRFUcWNaXA0gn6++TP9/duKaODirA3sCiiHgQQNJ0YCLgoGIt4dqqdacv34/BFog6PaiMAh4tjC8G9qnPJGkyMDmPPiXpN31c3zbA432cd6AbzNsGg3v7vG2daxud2zHb9+oymTo9qKhB2loXiSLiQuDCdV6ZNLdMm2InGszbBoN7+7xtnWswbl+ndyi5GNixMD4aWNKmspiZrfc6PajMAcZK2lnSRsBRwMw2l8nMbL3V0c1fEfGipJOAWaRbiqdFxMIWrnKdm9AGsMG8bTC4t8/b1rkG3fZ1/P9UzMxs4Oj05i8zMxtAHFTMzKwyDiolSJog6TeSFkma0u7yrCtJO0q6QdJ9khZKOjmnbyVptqQH8vuW7S5rX0kaIulOSVfn8Z0l3Z637fJ8Y0fHkTRC0hWS7s/7722DbL/9S/5O3iPpMkkbd+q+kzRN0jJJ9xTSGu4rJefnY8zdkvZqX8nXjYNKDwpdwRwM7AYcLWm39pZqnb0InBIRbwDGAyfmbZoCXB8RY4Hr83inOhm4rzB+LvC1vG0rgePbUqp193XgZxHxeuDNpG0cFPtN0ijgn4GuiNiDdPPNUXTuvrsYmFCX1mxfHQyMza/JwAX9VMbKOaj07C9dwUTEC0CtK5iOFRFLI+LXeXg16cA0irRdl+RslwCHt6eE60bSaOBQ4Lt5XMD+wBU5S0dum6ThwL7ARQAR8UJErGKQ7LdsKDBM0lBgE2ApHbrvIuJmYEVdcrN9NRG4NJLbgBGStu+fklbLQaVnjbqCGdWmslRO0hhgT+B2YLuIWAop8ADbtq9k6+Q/gM8CL+fxrYFVEfFiHu/UfbgLsBz479y0911JmzJI9ltE/AH4CvAIKZg8AcxjcOy7mmb7atAcZxxUelaqK5hOJGkz4IfAJyPiyXaXpwqS3gssi4h5xeQGWTtxHw4F9gIuiIg9gafp0KauRvL1hYnAzsAOwKakZqF6nbjvejJYvqMOKiUMyq5gJG1ICijfj4gf5eTHalXu/L6sXeVbB28HDpP0MKmpcn9SzWVEblKBzt2Hi4HFEXF7Hr+CFGQGw34DeDfwUEQsj4g/Az8C/prBse9qmu2rQXOccVDp2aDrCiZfY7gIuC8izitMmglMysOTgKv6u2zrKiJOi4jRETGGtK9+ERHHADcAR+RsnbptfwQelfS6nHQA6TEPHb/fskeA8ZI2yd/R2vZ1/L4raLavZgLH5rvAxgNP1JrJOo3/UV+CpENIZ7u1rmDObnOR1omkdwD/ByzglesOp5Ouq8wAdiL9wI+MiPoLjR1D0n7ApyPivZJ2IdVctgLuBD4UEc+3s3x9IWkc6QaEjYAHgeNIJ4eDYr9JOhP4AOkOxTuBj5KuLXTcvpN0GbAfqfv+x4CpwJU02Fc5iH6DdLfYM8BxEdG3R9S2mYOKmZlVxs1fZmZWGQcVMzOrjIOKmZlVxkHFzMwq46BiZmaVcVAxM7PKOKhYx5L0VD+s44uS3t3Hecfl/zgNGJJekjRf0g55fDNJF0j6Xe5PbJ6kE3pYxo2SDqpL+6Skb0naNS+/5fvGBiYHFbMmJA2JiDMi4ro+LmIcMKCCCvBsRIyLiFoXIN8ldSc/NvcnNoH0J8PuXEbqraDoKOCyiPhdRIyrtMTWURxUrONJ2i+fPdceXvX93N3FwZJm1OX7SR6+QNLc/ECoMwt5HpZ0hqRbgCMlXSzpiDztDElz8gOkLsz/gq6duZ8r6Q5Jv5X0ztylzxeBD+Qz9w/Ulfkjkq6U9BNJD0k6SdKncm3hNklb5Xwn5HXeJemHkjbJ6Ufmctwl6eactnsuw3ylBz2N7eFz25X0aIfPRcTLALnfrXMLeT6T13934XO6AnivpFflPGNIHUDe0qsdZ4OSg4oNFnsCnyQ9SG0XUseSs0l9SW2a83wAuDwP/2tEdAFvAt4l6U2FZT0XEe+IiOl16/hGRLw1P0BqGPDewrShEbF3LsPU/OydM4DLc83gcta2B/BB0oH9bOCZXFu4FTg25/lRXmftgVy1B1SdARyU0w/LaR8Hvp5rCl2kTgq7sztwVy2g1JN0IOmhUXuTal1vkbRvRPwJuINXHkB1VN5Od89hDio2aNwREYvzAXI+MCY/g+NnwPtyL7eH8koHfn8v6dekvqR2JwWjmkYBAOBvlB5ru4DU+/HuhWm1np7nAWNKlvmGiFgdEctJzw75SU5fUFjGHpL+L6/zmMI6fwlcnK9/DMlptwKnSzoVeHVEPFuyHABI+tdcy6k1jR2YX3cCvwZeTwoysGYT2FF53MxBxQaNYgeDL5GePQIpQPw9KQjMiYjVknYGPg0cEBFvAq4BNi7M/3T9wiVtDHwLOCIi3gh8p26e2vqL6+5NmV8ujL9cWMbFwEl5nWfW1hkRHwc+R+oufb6krSPiB6Ray7PALEn797D+e4E3S9ogL/PsXMsZXtts4Mu5pjUuIl4TERflaVcCByg9S31Y7UmiZg4qNtjdSHrmyAm8UgMZTgocT0jajsYPgqpXCyCPKz3c7IjuMmergc17Vdq1bQ4sVXr+zTG1REm7RsTtEXEG8DiwY+6J+cGIOJ/UlfqbGi4xi4hFwFzgS5KG5OVuzCsPjJoF/EPeXiSNkrRtnvcp0mc7DddSrMBBxQa1iHgJuJoUOK7OaXeRmnQWkg6KvyyxnFWk2skC0ln6nBKrvwHYrdGF+l74POmRBLOB+wvp/y5pgaR7gJuBu0jXjO6RNJ/UVHVpieV/lPS45UWS5gHXAacCRMTPgR8At+bmtytYM0heBryZ1C29GeCu783WK5KeiojNBst6bOBxTcVs/fJk8c+PVav9+ZH0UCpbD7mmYmZmlXFNxczMKuOgYmZmlXFQMTOzyjiomJlZZf4/HSrtBwH9grUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Write down there a code that will create the histogram.\n",
    "plt.hist(invariant_mass, bins=300, range=(0,110))\n",
    "print(\"Valor mínimo = \" + str(invariant_mass.min()) + \"\\nValor máximo = \" + str(invariant_mass.max()))\n",
    "\n",
    "# Let's name the axes and the title. Don't change these.\n",
    "plt.xlabel('Invariant mass [GeV]')\n",
    "plt.ylabel('Number of events')\n",
    "plt.title('Histogram of invariant mass values of two muons. \\n')\n",
    "#plt.xscale(\"log\")\n",
    "#plt.yscale(\"log\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Describe the histogram. What information you can get from it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fitting the function to the histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get information about mass and lifetime of the detected resonance, a function that describes the distribution of the invariant masses must be fitted to the values of the histogram. In our case the values follow a Breit-Wigner distribution:\n",
    "\n",
    "$$\n",
    "N(E) = \\frac{K}{(E-M)^2 + \\frac{\\Gamma^2}{4}},\n",
    "$$\n",
    "\n",
    "where $E$ is the energy, $M$ the maximum of the distribution (equals to the mass of the particle that is detected in the resonance), $\\Gamma$ the full width at half maximum (FWHM) or the decay width of the distribution and $K$ a constant.\n",
    "\n",
    "The Breit-Wigner distribution can also be expressed in the following form:\n",
    "\n",
    "$$\n",
    "\\frac{ \\frac{2\\sqrt{2}M\\Gamma\\sqrt{M^2(M^2+\\Gamma^2)} }{\\pi\\sqrt{M^2+\\sqrt{M^2(M^2+\\Gamma^2)}}} }{(E^2-M^2)^2 + M^2\\Gamma^2},\n",
    "$$\n",
    "\n",
    "where the constant $K$ is written open.\n",
    "\n",
    "The decay width $\\Gamma$ and the lifetime $\\tau$ of the particle detected in the resonance are related in the following way:\n",
    "\n",
    "$$\n",
    "\\Gamma \\equiv \\frac{\\hbar}{\\tau},\n",
    "$$\n",
    "\n",
    "where $\\hbar$ is the reduced Planck's constant.\n",
    "\n",
    "With the code below it is possible to optimize a function that represents Breit-Wigner distribution to the values of the histogram. The function is already written in the code. It is now your task to figure out which the values of the maximum of the distribution $M$ and the full width at half maximum of the distribution $\\Gamma$ could approximately be. The histogram that was created earlier will help in this task.\n",
    "\n",
    "Write these initial guesses in the code in the line `initials = [#THE INITIAL GUESS FOR GAMMA, #THE INITIAL GUESS FOR M, -2, 200, 13000]`. In other words replace the two comments in that line with the values that you derived.\n",
    "\n",
    "Notice that the initial guesses for parameters _a, b_ and _A_ have been already given. Other comments in the code can be left untouched. From them you can get information about what is happening in the code.\n",
    "\n",
    "After running the code Jupyter will print the values of the different parameters as a result of the optimization. Also uncertainties of the values and a graph of the fitted function are printed. The uncertainties will be received from the covariance matrix that the fitting function `curve_fit` will return."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    Think how M and gamma could be determined with the help of the histogram. Look from the histogram that you created that which would approximately be the values of M and gamma.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    If you figured out the initial guesses to be for example gamma = 12 and M = 1300 (note that these values are just random examples!) write them to the code in the form \"initials = [12, 1300, -2, 200, 13000]\".\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The values and the uncertainties from the optimization\n",
      "\n",
      "The value of the decay width (gamma) = 3.9240962076167643 +- 0.07790605300093356\n",
      "The value of the maximum of the distribution (M) = 90.79835898805646 +- 0.03051923389810049\n",
      "a = -1.6400575548826273 +- 0.12667564341450008\n",
      "b = 167.9533195718338 +- 12.787300621152154\n",
      "A = 13097.55327566484 +- 180.36534277205286\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEWCAYAAACe8xtsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3Xl4VOX1wPHvyQJhSdiSIBAgIFsI+yJYRUQFFCtUcEOr0KpYq3WrInXtr9YqVWuta0VRVAoqiiCgiKggLiC7bLIGEoIQIBBISMhyfn/cO+MEkjCBzEwSzud57pOZd+5ykkxy5l3u+4qqYowxxpRHWKgDMMYYU/VY8jDGGFNuljyMMcaUmyUPY4wx5WbJwxhjTLlZ8jDGGFNuljxMlSAifxWRd0IdR2lE5LCItA7g+V8RkYcDdf7yEJHGIrJQRA6JyDMi8oCIvBbquExwRYQ6AGPA+efr87Q2kAcUus9vqeBrvQlcCxwFFNgI3KOqC072nKpa95jzp6nqQ6VcvwmQDpyhqrvdsgeBv5dQ1k9VL1bVP5xsbAEwBtgLxOgxN4qJSCKwDYhU1YLgh2aCxWoeplJQ1bqeDdgBXOZTNjkAl/yne616wMvAhyISXtKOIlKhH7JUdRewGTjPp/g8YEMJZQsr8trlVcrPpCWw7tjEYU4vljxMVVJDRN5ym0vWikgvzwsi0lREPhCRDBHZJiJ3+HNCVS0C/gc0BBq75xotIt+IyLMish/4q1v+exFZLyKZIjJXRFr6XF9FpI2IjAGuA8a6TVkfl3LphbiJwv0H3R147piys939EJE3ReTvPtcbKyK7RCRdRG7yXN9n3xdFZLb7s1osImf6HNtBROaJyH4R+UlErvJ57U0ReVlE5ohINjDAN2i3VjXK5/u76JgmRU+yO+C+frY/vwdT9VjyMFXJUGAqUB+YCbwAICJhwMfAKqAZcCFwl4gMPtEJ3X/SN+A0tez2eakPsBWIBx4Xkd8ADwDDgTjga2DKsedT1VeBybg1G1W9rJRLe5MHTuLYAMw/piwSWFJCzBcD9wAXAW2A/iWcfyTwf0ADnFrO4+6xdYB5OAkz3t3vJRFJ9jn2Wnf/aGDRMd/f6GO+v8+Pua4n/vru69+V8v2bKs6Sh6lKFqnqHFUtBN4GurrlvYE4Vf2bqh5V1a3ABOCaMs51r4gcALKBfwMPu+f1SFfV51W1QFWP4PS7PKGq6922/H8A3XxrH+W0AOgkIg2AfsDXqroJiPUp+15Vj5Zw7FXAG6q6VlVzcJLEsT5U1SVurJOBbm75r4EUVX3D/d6WAx8AV/gcO0NVv1HVIlXNPcnvz1RzljxMVfKzz+McIMrtj2gJNBWRA54Np5bQuIxzPa2q9YFaQC/gKRG5xOf11GP2bwk853P+/YDg1HTKJCLXuU04h0XkEwBVTQHSgHNxPq1/7e7+nU9Zaf0dTY+J79hY4fifladDvyXQ55if1XXAGSc4nzHF2GgrUx2kAttUtW15D3Q7fdeIyDfApcAnnpdKuMbjfnbeFzvWPaak477GSRJn4/Qj+Jadi9ssV4JdQILP8+Z+xOSRCixQ1YFl7HMqHeHWiX6asJqHqQ6WAFkicr+I1BKRcBHpJCK9/TlYRDrg/LNeW8ZurwB/8fQNiEg9EbmylH13A/7c87EQp78lXVWz3LJFblk9nFpISd4DficiSSJSG3jEj2t5zALaicj1IhLpbr1FJKkc5yhLBlCEf9+/qcIseZgqz+2ruAynXX8bzj0Ir+H8Ay6NZ7RQNvAZ8Abw3zKuMR0YD0wVkSxgDXBJKbu/DnR0m4U+KiOGBTid1r6d0itxmtKWuf0ZJcXyCfAf4EucznBPkskr41qeYw8Bg3D6g9JxmrfGAzVPdKw/3JgfB75xv/++FXFeU/mIDdU2pmpzaw1rgJp2Y54JFqt5GFMFicjlIlLDHZk1HvjYEocJJksexlRNt+D0L2zBmcbl1tCGY0431mxljDGm3KzmYYwxptyq7X0esbGxmpiYGOowjDGmyli2bNleVY3zZ99qmzwSExNZunRpqMMwxpgqQ0S2+7uvNVsZY4wpN0sexhhjys2ShzHGmHKrNkN1ly1bFh8REfEa0AkI27dvX8smTZqEOixjjKl0oqKiSEhIIDIysli5iCxT1V6lHFZMtekwj4iIeO2MM85IiouLywwLC9N169a1TEqqqLnejDGmelBV9u3bR1paGq1atTrp81SnZqtOcXFxWWFhYdWjKmWMMQEgIjRq1Ijc3FNb56s6JY8wSxzGGHNiInLK56hOycOYqknV2YypQix5mFK98sordO7cmW7dunHuueeybt26Evd77rnn6NSpE8nJyfz73//2lt9333106NCBLl26cPnll3PgwIFix+3YsYO6devy9NNPe8s+/fRT2rdvT5s2bXjyySe95V988QU9evSgU6dOjBo1ioKCXyaQ/eqrr+jWrRvJycn0798fgNTUVAYMGEBSUhLJyck899xz3v3379/PwIEDadu2LQMHDiQzMxOAp556im7dutGtWzc6depEeHg4+/fvB+DAgQNcccUVdOjQgaSkJL77zllC4/333yc5OZmwsLDjbkpdvXo1Z599NsnJyXTu3Jnc3FwOHTrkvUa3bt2IjY3lruuug2XLYPly3nvqKTq6MV977bXFzpeVlUWzZs24/fbbvWUXX3wxXbt2JTk5mT/84Q8UFhZSljfffBMRYf78+d6y6dOnIyJMmzatzGP9cf/999OpUyc6derEu+++6y3ftm0bffr0oW3btlx99dUcPXr80uyTJ08u9rMJCwtj5cqVALz77rt06dKF5ORkxo4d6z1mx44dDBgwgO7du9OlSxfmzJkDwJIlS7zn6dq1K9OnT/cek5iY6H1f9+r1S9/wypUr6du3r7d8yZIlAGRmZnL55ZfTpUsXzjrrLNasWeM95tlnnyU5OZlOnToxcuRIb1PQ/Pnz6dGjh/dvZ/PmzQDk5eVx9dVX06ZNG/r06UNKSgoAKSkp1KpVyxvzH/7wB+81jh49ypgxY2jXrh0dOnTggw8+KPZzmzZtGiJS7P1X0nuvwqlqtdhWrlyZoqpLPdvatWvVnJqDBw96H8+YMUMHDx583D4//vijJicna3Z2tubn5+uFF16oGzduVFXVuXPnan5+vqqqjh07VseOHVvs2OHDh+sVV1yhTz31lKqqFhQUaOvWrXXLli2al5enXbp00bVr12phYaEmJCToTz/9pKqqDz/8sL722muqqpqZmalJSUm6fft2VVXdvXu3qqqmp6frsmXLVFU1KytL27Ztq573xH333adPPPGEqqo+8cQTx8Wlqjpz5kwdMGCA9/kNN9ygEyZMUFXVvLw8zczMVFXVdevW6YYNG7R///76ww8/ePfPz8/Xzp0768qVK1VVde/evVpQUHDcdXp06aIL/vtf1e3bdeNXX2m3du10/7p1xb4XjzvuuENHjhypt912m7fM8zsqKirS4cOH65QpU467hq833nhDO3furDfeeKO37KqrrtKuXbvq+++/X+axJzJr1iy96KKLND8/Xw8fPqw9e/b0xnfllVd6Y7vlllv0pZdeKvNcq1ev1latWqmq87Nr3ry57tmzR1Wd38Xnn3+uqqo333yz91xr167Vli1bqqp634+qznshLi7O+7xly5aakZFx3DUHDhyoc+bMUVXV2bNna//+/VVV9d5779W//vWvqqq6fv16veCCC1RVNS0tTRMTEzUnJ8f7Pb7xxhuqqtq2bVtd5/4eX3zxRR01apT38S233KKqqlOmTNGrrrpKVVW3bdumycnJJf4sHnnkEX3wwQdVVbWwsLBY7FlZWdqvXz/t06eP9/3n73vPE58vYKn6+T/Xah4V6LHHHqNDhw4MHDiQkSNHej9RT5gwgd69e9O1a1dGjBhBTo6zQNzo0aO59dZbGTBgAK1bt2bBggX8/ve/JykpidGjR3vPW7duXe6//3569uzJRRddxJIlSzj//PNp3bo1M2fOBJxPLv369aNHjx706NGDb7/99pS/n5iYGO/j7OzsEttJ169fT9++falduzYRERH079/f+ylv0KBBREQ4A/r69u1LWlqa97iPPvqI1q1bk5yc7C1bsmQJbdq0oXXr1tSoUYNrrrmGGTNmsG/fPmrWrEm7du0AGDhwoPfT1//+9z+GDx9OixYtAIiPjwegSZMm9OjRA4Do6GiSkpLYuXMnADNmzGDUKGfJ8FGjRvHRR8cv9jdlyhRGjhwJOJ/4Fy5cyI033ghAjRo1qF+/PgBJSUm0b9/+uOM/++wzunTpQteuXQFo1KgR4eHhxfbZtGkTezIy6NezJzRvzoTZs7lt5Ega1KhR7HsBWLZsGbt372bQoEHFzuH5HRUUFHD06FG/2rL79evHkiVLyM/P5/Dhw2zevJlu3bqd8LgTWbduHf379yciIoI6derQtWtXPv30U1SVL774giuuuAIo/Wfuy/fnv3XrVtq1a0dcnDPl0kUXXeT9/YsIWVnOCr4HDx6kadOmAN73I0Bubq5fP5fSzrVu3TouvPBCADp06EBKSgq7d+8GnJ/7kSNHKCgoICcnx3tMaefyfe9dccUVzJ8/Hz1Bk+XEiRP5y1/+AkBYWBixsbHe1x5++GHGjh1LVFSUt8yf915FqJ7J4/e/b95y9Gg4//yK2+66q8xLLl26lA8++IAVK1bw4YcfFqtCDh8+nB9++IFVq1aRlJTE66+/7n0tMzOTL774gmeffZbLLruMu+++m7Vr1/Ljjz96q+zZ2dmcf/75LFu2jOjoaB566CHmzZvH9OnTeeQRZ/nq+Ph45s2bx/Lly3n33Xe54447SoyzX79+xZoGPNvnn39e4v4vvvgiZ555JmPHjuU///nPca936tSJhQsXsm/fPnJycpgzZw6pqanH7Tdx4kQuueQS7/czfvx4Hn300WL77Ny5k+bNm3ufJyQksHPnTmJjY8nPz/f+TKdNm+a9xsaNG8nMzOT888+nZ8+evPXWW8ddOyUlhRUrVtCnTx8Adu/ejeceoCZNmrBnz55i++fk5PDpp58yYsQIwPnnFRcXx+9+9zu6d+/OTTfdRHZ2dok/L4+NGzciIgwePJgePXrwz3/+87h9pkyZwtWDByN16oAIGzduZGN6OudccQV9+/bl008/BaCoqIg///nPPPXUUyVea/DgwcTHxxMdHe39B10WEeGiiy5i7ty5zJgxg6FDh5a6r29Tnu9W0vura9eufPLJJ+Tk5LB3716+/PJLUlNT2bdvH/Xr1/f+M/f8Xsvy7rvvepNHmzZt2LBhAykpKRQUFPDRRx95f/9//etfeeedd0hISGDIkCE8//zz3nMsXrzY22zzyiuveK8vIgwaNIiePXvy6quvevf/97//zX333Ufz5s259957eeKJJ7zf14cffgg4H3C2b99OWloazZo1495776VFixY0adKEevXqeZP7a6+9xpAhQ0hISODtt99m3LhxQPH3eEREBPXq1WPfvn2A07TXvXt3+vfvz9dffw3gbep9+OGH6dGjB1deeaU3ca1YsYLU1FR+/etfF/vZ+fPeqwjVM3mEwKJFixg2bBi1atUiOjqayy67zPvamjVr6NevH507d2by5MmsXbvW+9pll12GiNC5c2caN25M586dCQsLIzk52dseWqNGDS6++GIAOnfuTP/+/YmMjKRz587effLz87n55pvp3LkzV155Zan9E19//TUrV648brvoootK3P+2225jy5YtjB8/nr///e/HvZ6UlMT999/PwIEDve3vnj9Sj8cff5yIiAiuu+46AB599FHuvvtu6tatW2y/kj6BiQgiwtSpU7n77rs566yziI6O9l6joKCAZcuWMXv2bObOnctjjz3Gxo0bvccfPnyYESNG8O9//7tYTaosH3/8Meeccw4NGzb0XmP58uXceuutrFixgjp16hTrjylJQUEBixYtYvLkySxatIjp06cX62cAmDp1KiMvugjq1PEesyk1la9efpkp77zDTTfdxIEDB3jppZcYMmRIscTqa+7cuezatYu8vDy++OILv77Ha665hqlTpzoxuP+kS3LfffeV+H4p6YPEoEGDGDJkCL/61a8YOXIkZ599NhEREaX+XkuzePFiateuTadOnQBo0KABL7/8MldffTX9+vUjMTHR+/ufMmUKo0ePJi0tjTlz5nD99ddTVFQEQJ8+fVi7di0//PADTzzxhLfd/5tvvmH58uV88sknvPjiiyxcuBCAl19+mWeffZbU1FSeffZZb01z3LhxZGZm0q1bN55//nm6d+9OREQEmZmZzJgxg23btpGenk52djbvvPMO4PSFzJkzh7S0NH73u99xzz33AKW/x5s0acKOHTtYsWIF//rXv7j22mvJysqioKCAtLQ0zjnnHJYvX87ZZ5/NvffeS1FREXfffTfPPPPMcefz571XEarNTYLFTJyYun3duviOHTsG7ZJlVT1Hjx7NRx99RNeuXXnzzTf56quvvK/VrFkTcKqjnsee555O4cjISO8fm+9+vvs8++yzNG7cmFWrVlFUVFSsGuurX79+HDp06Ljyp59+utQEAs4/m1tvLXmxuhtvvNH7h/bAAw+QkJDgfW3SpEnMmjWL+fPne7+HxYsXM23aNMaOHcuBAwcICwsjKiqKnj17Fqu1pKWleav7Z599tvfT2GeffeZNEAkJCcTGxlKnTh3q1KnDeeedx6pVq2jXrh35+fmMGDGC6667juHDh3vP27hxY3bt2kWTJk3YtWtXseYh4Lh/qAkJCSQkJHhrLldcccUJk0dCQgL9+/f3NjEMGTKE5cuXe5s/Vq1aRUF+Pj3bt4fatb3H9O3ShciICFo1bkz79u3ZtGkT3333HV9//TUvvfQShw8f5ujRo9StW7dYDFFRUQwdOpQZM2YwcODAMmMDvB2/tWrV8jYHluSpp55i8uTJx5Wfd955JSaQBx98kAcffBCAa6+9lrZt2xIbG8uBAwcoKCggIiKi2O+1JCUltMsuu8z7gezVV1/1NsO8/vrr3hra2WefTW5uLnv37i32O01KSqJOnTqsWbOGXr16ea8dHx/P5ZdfzpIlSzjvvPOYNGmSd2DFlVdeyU033QQ4TYNvvPEG4Pydt2rVilatWjF37lxatWrlbU4bPnw43377LYMHD2bVqlXe98vVV1/t/fCXkJBAamoqCQkJFBQUcPDgQRo2bIiIeP+ue/bsyZlnnsnGjRvp2bMntWvX5vLLL/fG9frrr3Po0CHWrFnD+eefD8DPP//M0KFDmTlz5gnfexXFah4V5Nxzz+Xjjz8mNzeXw4cPM3v2bO9rhw4dokmTJuTn55f4h1gRDh48SJMmTQgLC+Ptt98uddRNeWoemzZt8j6ePXs2bdu2LfGcnmafHTt28OGHH3r/8D/99FPGjx/PzJkzqe3+g/TEkJKSQkpKCnfddRcPPPAAt99+O71792bTpk1s27aNo0ePMnXqVG+TiucaeXl5jB8/3jsaZdiwYXz99dfeNufFixeTlJSEqnLjjTeSlJTk/dTnMXToUCZNmgQ4yW3YsGHFfo4LFiwoVnbGGWfQvHlzfvrpJ8AZSXOiDyaDBw9m9erV5OTkUFBQwIIFC4odM2XKFEb+5jfOE/dn0/Wci5j+hfMpeG9qKhs3bqR169ZMnjyZHTt2kJKSwtNPP80NN9zAk08+yeHDh9m1axfgfNqcM2cOHTp0AOCFF17ghRdeKDPGJ554gn/84x9l7lOemkdhYaG3CWb16tWsXr2aQYMGISIMGDDAO5rr2J+5r6KiIt5//32uueaaYuWe339mZiYvvfSS9x97ixYtvJ+q169fT25uLnFxcWzbts37wWr79u389NNPJCYmkp2d7f3wlJ2dzWeffeat4TRt2pQFCxYAzug+z/v9wIED3tFhr732Gueddx4xMTG0aNGC77//npycHFSV+fPnk5SURIMGDTh48KD3A868efPwzHbh+96bNm0aF1xwASJCRkaG929269atbNq0idatWyMiXHbZZd4PnJ73Xr169di7d6/376hv377MnDmTXr16nfC9V1GqZ80jBHr37s3QoUPp2rUrLVu2pFevXtSrVw9wOtL79OlDy5Yt6dy5c4mf/E/VH//4R0aMGMH777/PgAEDqOM2hZyKF154gc8//5zIyEgaNGjgfdOnp6dz0003eYdFjhgxgn379hEZGcmLL75IgwYNALj99tvJy8vzfhLu27cvr7zySqnXi4iI4IUXXmDw4MEUFhby+9//3tuh/tRTTzFr1iyKioq49dZbueCCCwDnU+XFF19Mly5dCAsL46abbqJTp04sWrSIt99+2zskE+Af//gHQ4YMYdy4cVx11VW8/vrrtGjRgvfff98bw/Tp0xk0aNBxP7/nn3+e6667jqNHj9K6dWvvJ9Hp06fzpz/9iYyMDC699FK6devG3LlzadCgAffccw+9e/dGRBgyZAiXXnqp93zvvfcec157DcLCwK0l/ur8C/l24ZckXXUVEZGRPPXUUzRq1KjUn1d2djZDhw4lLy+PwsJCLrjgAm9S3bBhA+ecc06Zv19PH1RFyc/Pp1+/foDzaf2dd97xNi+NHz+ea665hoceeoju3bt7a6ozZ85k6dKl/O1vfwNg4cKFJCQk0Lp162LnvvPOO1m1ahUAjzzyiLe29Mwzz3DzzTfz7LPPIiLeociLFi3iySefJDIykrCwMF566SViY2PZunWr91N8QUEB1157rbdWMGHCBO68804KCgqIiory9oesX7+eG264gfDwcDp27Ojts+zTpw9XXHEFPXr0ICIigu7duzNmzBgiIiKYMGECI0aMICwsjAYNGjBx4kTAqaVff/31tGnThoYNGzJ16lTv9/3II48QERFBeHg4r7zyirfZdPz48Vx//fXcddddxMXFed97pTnRe6+iVJuJEVetWpXStWvXvZ7n69at6xnMZitw2tfr1q1LTk4O5513Hq+++qp3xI8xx9mwwfnq1hZWpzmdo4mZ6cSEKfiMRCuvX//613z44YfUcEduGXOs9evXc+z8f6flxIiVwZgxY1i3bh25ubmMGjXKEocpnSrk5IDPsEuP3IiaxGRnQlGRUzM5CbNmzTrVCI0pU8CSh4g0B94CzgCKgFdV9TkRaQi8CyQCKcBVqpopTm/qc8AQIAcYrarL3XONAh5yT/13VZ0UqLhPxf/+979Qh2CqitxcJzn49AV5HIl0B04cOeIdiWVMZRPIDvMC4M+qmgT0BW4TkY7AOGC+qrYF5rvPAS4B2rrbGOBlADfZPAr0Ac4CHhWRBiVcr6ioqOjUZ/syJhjcG0VLSg7e5OHZx5gKVhHdFQFLHqq6y1NzUNVDwHqgGTAM8NQcJgHukBOGAW+5d8l/D9QXkSbAYGCequ5X1UxgHnBxCZdck5GRUc8SiKkSsrOLdZb7OhoeCeHhljxMQKi7nkdpw/n9FZQ+DxFJBLoDi4HGqroLnAQjIp4B2c0A31uT09yy0sqLKSgouGnx4sXzatWq1UlEUNUKmXbYmID4+Wfnq6fTHNidecT7eH3eQcjMtARiAsKzkuCpCHjyEJG6wAfAXaqaVcY/9JJe0DLKi+nZs+eenj17dvU879Wrlx47y6kxlUJREZx1FoweDT7TaVwy7pd7gx6eP4FrV35Kp7vfY8s/S58+xJhQCehNgiISiZM4Jqvqh27xbrc5CverZ2KhNMB3/oUEIL2McmOqpo0b4fBhKGM03prGZ1KrII/W+9JK3ceYUApY8nBHT70OrFfVf/m8NBMY5T4eBczwKb9BHH2Bg27z1lxgkIg0cDvKB7llxlRNnqaqzp1L3WVtY+cmuY57tgYjImPKLZDNVucA1wM/ishKt+wB4EngPRG5EdgBXOm+NgdnmO5mnKG6vwNQ1f0i8hjwg7vf31R1fwDjNiaw3MksSUwk0aepytf2+s6sv80P7g5SUMaUT8CSh6ououT+CoDjZuhyFyK5rZRzTQQmVlx0xoRQSoozRLeMqUfyImuSUbs+TbMygheXMeVgEyMaE2zbt0PLlnCC0YDpMXEkHNxT5j7GhIolD2OCLSUFEhNPuNvOmDiaWc3DVFKWPIwJNk/N4wR21ot3mq2qyeSlpnqx5GFMMGVlOTf/+VHzSI+Jo1ZBHuzde8J9jQk2Sx7GBNP27c5Xv5qt4osfY0wlYsnDmGDyDNP1s9kKsORhKiVLHsYEk889HieS5ql57NgRsHCMOVmWPIwJpu3bnZl04+NPuOvBqLpkR0ZZzcNUSpY8jAmmlBS/7vEAQMTp97DkYSohW4bWmGDavp2FebW5oZRpSY61s14c7Sx5mErIah7GBFNKCmn1Gvu9e3pMnPV5mErJkocxwZKdDXv3klbvxP0dHjtj4mHfPudYYyoRSx7GBIvb/FSe5JFmw3VNJWXJw5hgcYfppsX432y104brmkrKkocxweLWHnbWi/P7kPSYuGLHGlNZBHIlwYkiskdE1viUvSsiK90txbNIlIgkisgRn9de8Tmmp4j8KCKbReQ/UsYi6MZUaikpEBnJnroN/T5kd92G5IeF8+Kb80tdOMqYUAjkUN03gReAtzwFqnq157GIPAMc9Nl/i6p2K+E8LwNjgO9xVhu8GPgkAPEaE1juPR4q/n9mKwoL5+foWJpl2boepnIJWM1DVRcCJS4X69YergKmlHUOEWkCxKjqd+5Kg28Bv6noWI0JCj+nYj9WekycrShoKp1Q9Xn0A3ar6iafslYiskJEFohIP7esGZDms0+aW1YiERkjIktFZGlGhv2xmUrGz0WgjpUWE0ezg/Z+NpVLqJLHSIrXOnYBLVS1O3AP8D8RiaHkNdBLXRlHVV9V1V6q2isuzv9OSWMC7sgR2L37pGoeO2PiOePwPiIKCwIQmDEnJ+jTk4hIBDAc6OkpU9U8IM99vExEtgDtcGoaCT6HJwDpwYvWmAriGWqbmAhry3doekwc4VpE48MltgIbExKhqHlcBGxQVW9zlIjEiUi4+7g10BbYqqq7gEMi0tftJ7kBmBGCmI05NeVYx+NYnnU9rNPcVCaBHKo7BfgOaC8iaSJyo/vSNRzfUX4esFpEVgHTgD/tjkGLAAAgAElEQVSoqudj1q3Aa8BmYAs20spURampztcWLcp9qOdGwWYHLXmYyiNgzVaqOrKU8tEllH0AfFDK/kuBThUanDFB9twbX/AnhHYvroLw8v3ZpcfEAtDkkK1lbioPu8PcmCBociiDjLoNKChn4gDIjYwiMyrakoepVCx5GBMETbMySI8++RGAP0c34gxLHqYSseRhTBA0ObTX2/x0MpwbBS15mMrDkocxgaZK06y97Io++eRhNQ9T2VjyMCbQ9u+nVkEeu2JOvtlqV3QsjY5kQW5uBQZmzMmz5GFMoLnDdNNPqebhHrtzZ0VEZMwps+RhTKC5yeNUax6+5zIm1Cx5GBNoFVDz8CaPtLSydzQmSCx5GBNoqankh4Wzt079kz7Fz9GNnAeWPEwlYcnDmEBLTWV33UYUhYWf9ClyatTiYM06ljxMpWHJw5hAS01l5yn0d3jsio61Pg9TaVjyMCbQUlPZdQo3CHrsiom1moepNCx5GBNIRUWQlsauU5iaxGNXdCwZG7aSOG42ieNmV0Bwxpw8Sx7GBNKePZCff0pTk3j8HB1LXM4BahTkV0BgxpwaSx7GBJK7gmDF1DycEVfxh/ed8rmMOVUnTB4iUtOfshL2mSgie0RkjU/ZX0Vkp4isdLchPq/9RUQ2i8hPIjLYp/xit2yziIzz79syppLw3iBYAX0ebgJqanNcmUrAn5rHd36WHetN4OISyp9V1W7uNgdARDrirDCY7B7zkoiEu0vTvghcAnQERrr7GlM1VMANgh6eGwVtgkRTGZS6Mo2InAE0A2qJSHdA3JdigNonOrGqLhSRRD/jGAZMVdU8YJuIbAbOcl/brKpb3Zimuvuu8/O8xoRWaipERZFZK+aUT+W5UdAWhTKVQVnLmg0GRgMJwL98yg8BD5zCNW8XkRuApcCfVTUTJ0l977NPmlsGkHpMeZ9TuLYxwZWa6qxbLnLifU8gu2ZtsmrU5oxD1udhQq/U5KGqk4BJIjLCXWO8IrwMPAao+/UZ4Pf8UqspFgIlN6tpaScXkTHAGIAWLVqcaqzGnLrUVGjevMJO93N0LE2zMirsfMacLH8WVJ4lItcCib77q+rfynsxVd3teSwiE4BZ7tM0wPcvLAFIdx+XVl7S+V8FXgXo1atXqUnGmKBJTYWBAyvsdLtiYq3mYSoFfzrMZ+D0MxQA2T5buYlIE5+nlwOekVgzgWtEpKaItALaAkuAH4C2ItJKRGrgdKrPPJlrGxN0BQWwa1eF1jx2Rcdan4epFPypeSSoakmjpsokIlOA84FYEUkDHgXOF5FuOE1PKcAtAKq6VkTew+kILwBuU9VC9zy3A3OBcGCiqq4tbyzGhER6unOHefPmsKViTvlzdCNisw8QWWg3CprQ8id5fCsinVX1x/KcWFVHllD8ehn7Pw48XkL5HGBOea5tTGUw4oF3+QAYNW8XtG5aIefcFR1LGErjw/sr5HzGnCx/kse5wGgR2Qbk4XRuq6p2CWhkxlRxno7tirjHw8Pu9TCVhT/J45KAR2FMNZSQtQeAnfXiK+ycnuTRJMuShwmtE3aYq+p2nBFPF7iPc/w5zpjTXcLB3eyrFUNOjVoVds6f3WlOrNPchJo/c1s9CtwP/MUtigTeCWRQxlQHCQf3kFavcYWe81DNOmTVqE0zt1ZjTKj4U4O4HBiKOzxXVdOB6EAGZUx14CSPimuy8thZL55mBy15mNDyJ3kcVVXFvbNbROoENiRjqgFVmmVVfM0DYGdMHM3sLnMTYv4kj/dE5L9AfRG5GfgcmBDYsIyp4vbsIargaIWsXX6snfXiLXmYkDvhaCtVfVpEBgJZQHvgEVWdF/DIjKnKUlIAAlTziCcmLxsOHoR69Sr8/Mb444TJQ0TuBt63hGFMOXiTRwD6PGLcc27fDl3sdisTGv7c5xEDzBWR/cBUYJrvBIfGmBK4ycP7j74Cee4bufFv05jfxlmxIOXJSyv8OsaUxZ/7PP5PVZOB24CmwAIR+TzgkRlTlaWksL9WDNk1T7huWrl5EpKNuDKhVJ6b/fYAPwP7gIr/OGVMdZKSEpAmK4C9deqRFx5pneYmpPy5SfBWEfkKmA/EAjfbvFbGnEBKCmkBaLICUAlzhutazcOEkD99Hi2Bu1R1ZaCDMaZaUIWUFHZ27hCwS+yMiffOnWVMKPjT5zEOqCsivwMQkTh3wSZjTEn27IHc3IAM0/Vw7vWw5GFCx+a2Mqaibd8OBGaYrsfOmDjisg9Qs+BowK5hTFkCNreViEwUkT0issan7CkR2SAiq0VkuojUd8sTReSIiKx0t1d8jukpIj+KyGYR+Y+ISHm/SWOCKoA3CHp4hus2tU5zEyKBnNvqTeDY5WvnAZ3cDveN/FKbAdiiqt3c7Q8+5S8DY3DWNW9bwjmNqVwCeI+Hhw3XNaEWsLmtVHUhsP+Yss9UtcB9+j2QUNY5RKQJEKOq37kJ7C3gN37EbEzopKRAw4YcDsA9Hh6emof1e5hQCeXcVr8H3vV53kpEVrjXeUhVvwaaAWk++6S5ZSUSkTE4tRRatGhRASEacxJSUiAxMaCX+LluIwolzGoeJmT8GaqLmywqbG4rEXkQKAAmu0W7gBaquk9EegIfiUgyznrpx4VTRpyvAq8C9OrVq9T9jAmolBToELhhugAF4RH8XLeR1TxMyAR9OVkRGQX8GrjObYpCVfNUdZ/7eBmwBWiHU9PwbdpKANKDG7Ex5eDe4xHomgfAznpxJFiHuQmRoCYPEbkYZ9jvUFXN8SmPE5Fw93FrnI7xraq6CzgkIn3dUVY3ADOCGbMx5dHzjv/BkSP8dXV2wK+1M8ZWFDShU2ryEJH57tfxJ3NiEZkCfAe0F5E0EbkReAFnmO+8Y4bkngesFpFVwDTgD6rq6Wy/FXgN2IxTI/nkZOIxJhgSDjoTTgdymK7HznrxnHFoL+FFhQG/ljHHKqvPo4mI9AeGishUjul/UNXlZZ1YVUeWUPx6Kft+AHxQymtLgU5lXcuYysIzWWEgbxD02BkTT4QW0fjwvoBfy5hjlZU8HgHG4fQz/OuY1xS4IFBBGVNVJWY6XXKpwah5uEvcWtOVCYVSk4eqTgOmicjDqvpYEGMypspqtT+d3XUbBmQdj2P9cq+HdZqb4PPnPo/HRGQoTr8EwFeqOiuwYRlTNbXK3MnWhqXeilShPDUPTz+LMcHkz8SITwB3Auvc7U63zBhzjFb7d7KtQXCSR25kFLvrNqRl5s9BuZ4xvvy5SfBSoJuqFgGIyCRgBcXnpTLG7N9PoyNZbG3YNGiXTGnQ1NvPYkww+XufR32fx/UCEYgxVd7GjQBsC1KzFcA2Sx4mRPypeTwBrBCRL3GG656H1TqMOZ4neQSp2QqcmkdczgE4eBDq2ec6Ezz+rCQ4BegLfOhuZ6vq1EAHZkyVs3EjBRJGav3AD9P12NbAbSLbtClo1zQG/J8YcRcwM8CxGFO1bdxIav3G5IdHBu2S2xr6JI9evYJ2XWOCPjGiMdXWxo1BbbIC2F6/ifPAah4myCx5GFMRiopg06agdpYD5EXWZGd0nCUPE3RlJg8RCfNdg9wYU4r0dMjJCXryAEhp2MSShwm6MpOHe2/HKhGxZfmMKYs70mprg+Dd4+GR0qCpJQ8TdP50mDcB1orIEsC7SIGqDg1YVMZUNZ7k0TDhBDtWvG0NmsLK/bB/PzRsGPTrm9OTP8nj/wIehTFV3caNULs2u6OD/887xdNJv2kT9OkT9Oub05M/93ksAFKASPfxD0CZa3l4iMhEEdnj228iIg1FZJ6IbHK/NnDLRUT+IyKbRWS1iPTwOWaUu/8mdxlbYyqXjRuhbVtUgj8GpdhwXWOCxJ+JEW/GWd3vv25RM+AjP8//JnDxMWXjgPmq2haY7z4HuARn+dm2wBjgZff6DYFHgT7AWcCjnoRjTKWxcSO0axeSS++ofwaEhXmbzowJBn8+Jt0GnANkAajqJsCvZdJUdSGw/5jiYcAk9/Ek4Dc+5W+p43ugvog0AQYD81R1v6pmAvM4PiEZEzr5+bB1a8iSR354JDui45g5bQGJ42aHJAZz+vEneeSp6lHPExGJwFlJ8GQ1du9Y99y57klEzYBUn/3S3LLSyo8jImNEZKmILM3IsAVyTHBc8MeJUFjIPStzQhaDza5rgs2f5LFARB4AaonIQOB94OMAxCIllGkZ5ccXqr6qqr1UtVdcXFyFBmdMaVpl7gSCOyHisbY1bEri/nTQU/lcZ4z//Eke44AM4EfgFmAO8NApXHO32xyF+9WzAHMa0NxnvwQgvYxyYyqFM/elAQRtBcGSpDRoSszRHBrlHAxZDOb04s9oqyKcvonHcIbtTlI9pY83MwHPiKlRwAyf8hvcUVd9gYNus9ZcYJCINHA7yge5ZcZUCkl7tpEeHcvBWtEhi8Ezu641XZlgOeF9HiJyKfAKsAWnCamViNyiqp/4cewU4HwgVkTScEZNPQm8JyI3AjuAK93d5wBDgM1ADvA7AFXdLyKP4QwRBvibqh7bCW9MyHTcs5X18a1CGkOKmzxaWfIwQeLPTYLPAANUdTOAiJwJzAZOmDxUdWQpL11Ywr6KM7KrpPNMBCb6EasxwZWby5n70vis7dkhDSOtXmPyw8JptX9nSOMwpw9/+jz2eBKHayu/9FMYc3pbt44ILQp5zaMgPIKtDZvRbu/2kMZhTh+l1jxEZLj7cK2IzAHewxnldCW/NCEZc3pbuRKAdSFOHgAb4lrRc+e6UIdhThNl1Twuc7coYDfQH6f/IgOwO7yNAVi1iuzIKLY3aBLqSNgQn0hCVgYcOBDqUMxpoNSah6r+LpiBGFMlrVzJhrjEkMxpdaz1cW7t58cfoV+/0AZjqj1/Rlu1Av4EJPrub1Oym9OeKqxaxfpWvwp1JACsj090HqxaZcnDBJw/o60+Al7Huau8KLDhGFOFbN8OBw+GvLPcY3fdRmRGRdNg9epQh2JOA/4kj1xV/U/AIzGmCkkcN5uBm75nArAuvnWow3GIsD6+Fb9atSrUkZjTgD/J4zkReRT4DMjzFKqqX2t6GFNdJe3ZRhHChrjEUIfitSEukV+t+RwKCyE8PNThmGrMn+TRGbgeuIBfmq3UfW7Maavjnq1sa9iUIzWiQh2K1/r4RFiW40wR37ZtqMMx1Zg/yeNyoLXvtOzGGKfm8WPjNqEOo5j1nia0VasseZiA8md84SqgfqADMaYqqZuXQ8sDP7OucSXp73BtatTcWVXQOs1NgPlT82gMbBCRHyje52FDdc1pq0PGNoBKM9LKIy+yJrRv79Q8jAkgf5LHowGPwpgqpuPurUDlmJbkWB9rLN0WLKbfuNmkPHlpqMMx1dQJk4eqLghGIMZUJT3SN7C7bkN2120U6lCOsz6+FZdt+JrovOxQh2KqMX/uMD/EL8u+1gAigWxVjQlkYMZUZr3S1vFDs44gJa2SHFqeprT2GSmhDcRUa/6sJBitqjHuFgWMAF4IfGjGVFKpqSRkZbA0oWOoIymR576TDntSQhqHqd7KPZubqn7EKdzjISLtRWSlz5YlIneJyF9FZKdP+RCfY/4iIptF5CcRGXyy1zamQnzzDQA/VNLksSs6lgNRdUnevSXUoZhqzJ9mq+E+T8OAXvzSjFVuqvoT0M09dziwE5iOs+zss6r69DHX7whcAyQDTYHPRaSdqhaebAzGnJJFizhcoxYbKmFnOQAirGzSnp47N4Q6ElON+VPzuMxnGwwcAoZV0PUvBLaoalnLnw0Dpqpqnqpuw1nj/KwKur4x5bdoEcubdqAwrPJO/7E0IYl2+3bA/v2hDsVUU/6Mtgrkuh7XAFN8nt8uIjcAS4E/q2om0Az43mefNLfsOCIyBhgD0KJFi4AEbE5zBw/C6tUsPefaUEdSpmXN3Ca1776DS224rql4ZS1D+0gZx6mqPnYqFxaRGsBQ4C9u0cvAYzhNYo8BzwC/B0oazlJis5mqvgq8CtCrV6+TblozplTffQeqlba/w2Nlk3YUSBgR33xjycMERFnNVtklbAA3AvdXwLUvAZar6m4AVd2tqoWqWgRM4JemqTSguc9xCUB6BVzfmPJbtAjCw1nZpH2oIynTkRpRrG18prdz35iKVmryUNVnPBvOp/laOJ3aU4GKmNBnJD5NViLiuwj05cAa9/FM4BoRqemuatgWWFIB1zemXBLHzeb7t2awMr51pZpJtzTLmiXBkiVw1OY0NRWvzA5zEWkoIn8HVuM0cfVQ1ftVdc+pXFREagMDgQ99iv8pIj+KyGpgAHA3gKquBd4D1gGfArfZSCsTCpGF+XTbtZGlzSp3k5XH0oSOkJsLK1aEOhRTDZXV5/EUMByn1tFZVQ9X1EVVNQdodEzZ9WXs/zjweEVd35iT0ennLUQVHOWHhORQh+KXpc2SnAfffAN9+oQ2GFPtlFXz+DPOfRUPAenuzXxZInJIRLKCE54xlUfvtLUALEtICnEk/tkT3QgSE63fwwREqTUPVS333efGVGf9ty3jp9gW7K3TINSh+O+cc+Dzz0G1Us7DZaouSxDG+OPgQc5KXcv8NlXs/tRzzoHdu51laY2pQJY8jPHH3LlEFhUy/8yqlTwu/sEZW3L3nS+FOBJT3VjyMMYfs2aRGRXNiqaV+/6OY22MbUFWjdre/hpjKoolD2NOpLAQ5szhyzN7UVSJ57MqSVFYON+37EL/rcudfg9jKoglD2NOZPFi2LePL87sHepITsoXrXvR7FAGrFlz4p2N8ZMlD2NOZNYsiIhgYaseoY7kpHx5Zi/nwezZoQ3EVCuWPIw5kVmzoF8/sqLqhjqSk7I7OpY1jc+05GEqlCUPY8qyfTv8+CP8+tehjuSUfNG6F3z7ra3vYSqMJQ9jyvDwzU8CMGBD1ax1eHx5Zm8oKoK5c0MdiqkmLHkYU4bL1i9kU6PmbGtY4vpjVcaqJm0hNtaarkyFseRhTGk2beKstHV80OnCUEdyyorCwvmgcWf2f/gxrcfODHU4phqw5GFMaSZNolDC+DB5QKgjqRBfntmbhkey6Ja+MdShmGrAkocxJSkshEmTWNiquzM7bTWwsFUPCiSMC7b+EOpQTDUQsuQhIinu4k8rRWSpW9ZQROaJyCb3awO3XETkPyKyWURWi0jVHHBvqo4vv4S0NKZ1uijUkVSYrKi6LG7RiSEbFtnd5uaUhbrmMUBVu6mqexcT44D5qtoWmO8+B2e987buNgZ4OeiRmtPLm29C/fp83rZ6LaI0PfkCWmemw3ffhToUU8WFOnkcaxgwyX08CfiNT/lb6vgeqH/MmufGVJyDB+HDD2HkSPIiaoQ6mgr1SbtfkRNZEyZNOvHOxpQhlMlDgc9EZJmIjHHLGqvqLgD3a7xb3gxI9Tk2zS0rRkTGiMhSEVmakZERwNBNdXb/NQ/DkSMMO9wm1KFUuOyatfmk/Tnw7rtw5EiowzFVWCiTxzmq2gOnSeo2ETmvjH1LWgLtuEZbVX1VVXupaq+4uLiKitOcTlQZvexjNsS2ZFWTdqGOJiCmdbrQqV3NmBHqUEwVFrLkoarp7tc9wHTgLGC3pznK/brH3T0NaO5zeAKQHrxozWljzhySMlL4b58R1XbZ1u9bdIYWLazpypySkCQPEakjItGex8AgYA0wExjl7jYK8Hw0mgnc4I666gsc9DRvGVMREsfNJnHcbH64ZSxpMfF8nFRWRbhqUwmD66+Hzz6DdPsMZk5OqGoejYFFIrIKWALMVtVPgSeBgSKyCRjoPgeYA2wFNgMTgD8GP2RT3fVKW0vvnet49azLKQiPCHU4gXXDDc5cV5MnhzoSU0WF5C9EVbcCXUso3wccNxeEqipwWxBCM6exP373Pntr1+O9LgNDHUrAJU7cxPvNOhL3j2dJvOceCK9aKySa0KtsQ3WNCYkOe7ZxwdalvNFzKLmRUaEOJyhe6/0bEg/sgg8+CHUopgqy5GEMcOc3UzhUoxZv97g01KEEzby2fdjSsBmMH293nJtys+RhzBdfcMnGb/lvnxFVdrXAk1EUFs5/zxoBy5fD55+HOhxTxVjyMKe3ggK44w521GvMhLOGhzqaoPsoeQA0bQpPPnninY3xYcnDnN5efhnWruXvF9xU7aYi8cfRiEi4+2744gv4wWbbNf6z5GFOXxkZ8MgjcNFFfNa2b6ijCZ0xY6BePXj88VBHYqoQSx7m9HX//XDoEDz3XLW9m9wvMTHw5z8705V8+WWoozFVhCUPc3r64AN44w1e6n05iW9tC3U0oXfvvdCyJdx1l9MPZMwJWPIwp5/UVLj5ZlY2acu/zv1tqKOpHGrVgqefhtWrYcKEUEdjqgBLHub0UljoTM1x9Ch3XnZf9Z+GpDxGjIDzz4eHHoL9+0Mdjank7C/HnF6eeAK++grefJPt62NDHU2lkDhutvdxh1Yj+HThQmcgwQsvhDAqU9lZzcOcNm4fdj88/DDTO55P4rpGoQ6nUtoQ3wpuuw1efNFuHDRlsuRhTg9ff80zs//F4oRkxl1yx+k9uupEnnwSkpKc5r29e0MdjamkLHmYasezNoe3OWb9ehg2jNR6ZzBm+EOn5c2A5ZH4ty8Z0udW8vZkwE032bxXpkSWPEz1tmEDDBwINWow+qr/42Ct6FBHVCWsa9yaf543yrn349VXQx2OqYSCnjxEpLmIfCki60VkrYjc6Zb/VUR2ishKdxvic8xfRGSziPwkIoODHbOpmpJ3b4F+/Zz7FubNI61e41CHVKVM7D0MBg2CO+6AhQtDHY6pZEJR8ygA/qyqSUBf4DYR6ei+9qyqdnO3OQDua9cAycDFwEsiYivXmDL1TFvHlCkPkHZUOH/o30icvCPUIVU5KmEwZQq0agW/+Y1TizPGFfShuu7a47vcx4dEZD3QrIxDhgFTVTUP2CYim4GzgO8CHqypMnyHm1616jMem/cSO2Pi+e01fyc9Jj6EkVVxDRvCJ59A374wZAh89x00thqcCXGfh4gkAt2BxW7R7SKyWkQmikgDt6wZkOpzWBqlJBsRGSMiS0VkaUZGRoCiNpVVjYJ8/vHpC/zz0/+wJKETI377lCWOU5Q4bjaJ/13H0IvHwc8/w6WX2g2EBghh8hCRusAHwF2qmgW8DJwJdMOpmTzj2bWEw0sc/qGqr6pqL1XtFRcXF4CoTWXVNmM7708ey7WrPuXlPlcw6qr/I7N2vVCHVW2sbtIO3n8ffvwRBgyAPXtCHZIJsZDcYS4ikTiJY7Kqfgigqrt9Xp8AzHKfpgHNfQ5PANKDFKqpxBLHzSa8qJBbFn/ArG/+R3aN2tzymweY2/5XoQ6terr0Upg1y+n/OO885ybChIRQR2VCJBSjrQR4HVivqv/yKW/is9vlwBr38UzgGhGpKSKtgLbAkmDFayqvPjt+ZMZb9zB24Vt83qYPA298yRJHoA0cCHPnwq5dcM45sGJFqCMyIRKKmsc5wPXAjyKy0i17ABgpIt1wmqRSgFsAVHWtiLwHrMMZqXWbqhYGPWoTMr6d4SlPXgqbN8PYsbw7fTrp0bH8cdg45nQ4N4QRnmbOPddZ92PYMPjVr+C11+C660IdlQmyUIy2WkTJ/RhzyjjmccCWOTvNnbk3lQ86XcCwdQvIi6jBS/2u5/Xew8iNjAp1aKefHj1g2TK48kr47W9hyRIYPx6i7HdxurBZdU2l5K1tqHL2jh8ZtfxjBm38ntzIGrzZ8zL+2+cKMuo2KPskJrDi451+j3vvhf/8x3k8aRL06hXqyEwQWPIwlVKj7AMMXb+A61Z8Qpv9aWRGRfPi2VcxsddQG0UVIr7Nh75SnnsOLrnEmQerb19ned8HH4TatYMcoQkmSx6m0uhy17tcuHkJw9YtYHHKCiK0iBVN2nPPpXczu/255EXWDHWIpjQXXwxr1jjL2P7jH/D2204z1jXX2AzG1ZQlDxM6qvDTT/DZZzBjBsu//IoILSItJp7/9hnBRx3PZ1Ncy1BHafxVvz68+SbceCPceSdcey08/zz87W9w4YWWRKoZ0Wo63XKvXr106dKloQ7jtFNq08aTlzrJYssW+PprZ6K9zz+HtDRnh6QkXqqXzLw2fVnZtJ0zr5KpclKevNR5UFjoJJJHH4WdO51RWQ8/DIMHWxKpxERkmar61WllNQ9zykpMGKrEZR+g456tdNm1EYb+1xmRs9u5FzQzKppvW3bhm8HD+DqxO6n1zwhy1CagwsOdGshvfwsTJzrL/15yCXTsCH/6E1x/PdSpE+oozSmwmoc5Kb6joeKyM2m9fydt9qU6295UkjK2EZtzEIAihK0Nm7GqSVuWJiTzQ0JHtjRKsNpFNeethQDk5cHUqU4z1rJlUK8eXH01jB7tdLJbbaRSKE/Nw5KHKZuqsxRpaqqz7dgBKSl8Ous7Whz4mZaZu6iTn+vdPTsyii2NEtgQl8j6+Fasj2/FmsZtOFzTRt6cboolDw9V+P57ePll+OADyMmBtm3hqqvgiiuga1dLJCFkyQNLHmU6csSZGXXfPmfLyHC2PXtg927mfrGK2OwDND68n4TcA3D0aPHDI2qSVi+e1HqN2d6gCdsaNGV7g6Zsim3OruhYq1GYMnmTyqFDMG2aMzJrwQIoKoIzz4Rf/9qZ/r1/f6hpI+yCyZIH1SR5qEJ+vvPPPjfX+ZqT43zNzmbUC19RKz+X2vm51M7Po/bRI9Q5msudZzV2/jCzsuDgQZav2U5MbjYxednEF+Q4TQglKELIrBXN3jr1yajTgN11G5JRtyG76zYkPTqO9Jg4dsXEsrd2fft0aE5aiTWSjAz46COYPt2Z+iQ317lPpF8/uOACZ+vWDSKsmzaQLHlwislj6VJnJJAqT32yHkEJU+WeC9s4n448W2Hhcdvkb7cSXlTkbFpIeFEREUUFXNoh1lkONT//l+3oUedrXp7z2Odr7uEcahTkE1by7POlKkLIrhFFdo1aZNeozaGatcmqWUo/Ey4AAAz5SURBVIdDNWtzMKouWVF1yapZh8xaMWTWiuZArWj21apHZu0YMmvFUBhmizSa0PAmlZwc+OorZxGqL76Adeuc8jp1oE8fZ0LGs86C3r1tYaoKZsmDU0we//433H136a+Hh3NUBRWhMCyM2lE1nNEl4eFk5BZSKGEUhIVTJGEUhoVREBZBQVi4u0XQ/cw4iIyEGjWYtyWTgrAIjoZHMqxPK6hRA6KieOX7NI6GR5IbWZO88EjyImqQF1GDI5E1yY2oyZHIX7bsyFrk1IgiJzKKI5E1rdnIVEkl1kjAWYTqq6/gm2+cbdUq58MbQPPmzjxbXbs6W+fO0Lq18/doys2SB6eWPNrcN4PwokJUwlABxUkSilhzjTFBdlxSOXzYmQr+hx+cVoKVK52bTT0JJSoKkpKcrX176NAB2rVz+lOio4P/DVQhljw4teRR2o1uxpjQ8k0kxabqf2SAMz3K2rVMeHkm7fbu4Mx9aSQcynD6Dj3i452aSatWzpaYCC1bQosWznaaz8dlyQNLHsYYN6ls2gSbN/Pk87NIzEynxYGf+VVYljP0vPCYpYEaNHCawpo14//bO/fgKqo7jn++BgMJSSCAPKVFKeOLIlqk9IVWrYK1Wjs+sO2o01bHzjBT29pRq6XK1Glp6zjaqh2rFLUVUXQsvl+Voo4oPkCCYgWpFqW85JEXCYRf/zjnlk24CblA7t4kv8/Mzu6eu3vOd09O9rfn7Nnfj6FDw3rIEBg8OKwHDQpLSUk6F9TB+BfmjuM4wIjpz8etXjDh7Ga/9WjawaCaT3jp/JHwwQewejV3P/ASQ6rXM+jN9xi8YCED6zY377lkKC8PvZiBA3l6g7GxpIJNpRV8UlLBppI+bCopZ+ZPJ0H//tCvX/D71cVminWaq5E0CbgJKALuMLPfpCzJcZxOzI6iHnzUZyAjHq8G+oXllDHNjina2UT/2s0cVLuJgbWbGFC7iYNqNzOgdjMD6jbT7+MtDK/bwtiP36VvfTXFO3fsOvnB6c3y2lpcytZeZRx86NDQw6msDEalb9/wxX1mXVER1pmloiIsPXs2e+e6W4TNPNMpjIekIuAW4GvAamCRpHlm9na6yhzH6co0HVDEuvL+rCvvz7I9HWxGeWMdlXVbqazfSmV9NX23VdNnWw1943bFthqWV9fQZ/1aKhpWUrGtlj4NNZRuz/7tVZLtBxRRU1xKbXEJNcUlPFhc8v/tuY/dSE1xCfOOPJ6H7rl8v1z7nugUxgMYD6wws/cBJN0HnEmIa+44jpM+EtU9e1PdszcfVg7J6dQeTTsoa6yjYlst5Q21lDfWUd5QR1lDHWWZ7cY6ejfWx7R6Shu3Ud5Qx+DqjZRur6essZ7FQw/roIvLojlvJe0bw4D/JPZXA59veZCkS4BL4m6NpHf3srwBwIa9PLcjcV254bpyw3XlRuHpevQGbtIN+6Kr3QF0OovxyPZxxW5vsczsduD2fS5Meq29Mw7yievKDdeVG64rN7q7rs7yKfJqYHhi/2Dg45S0OI7jdHs6i/FYBIySdIikYmAKMC9lTY7jON2WTjFsZWY7JE0FniJM1Z1pZnuc/LAP7PPQVwfhunLDdeWG68qNbq2ry35h7jiO43QcnWXYynEcxykg3Hg4juM4OdOtjYekwyQtTixbJV0mqZ+kZyS9F9eVBaLrWkkfJdJPy6euqO3HkpZJqpI0W1KvOJHhlVhfc+KkhkLQNUvSqkR9jU1B14+ipmWSLotpqbavNnSl0r4kzZS0TlJVIi1rHSlws6QVkt6SdGyB6DpB0pZE3U3Ls65z4t9yp6RxLY6/KtbXu5JO3W9CzMyX8N6nCPgv4SOZ3wJXxvQrgRkFouta4PIUtQwDVgElcf9+4KK4nhLT/gT8sEB0zQLOTrG+RgNVQClhcsqzwKi021cbulJpX8BE4FigKpGWtY6A04AnCN9+TQBeKRBdJwCPplhfRwCHAfOBcYn0I4ElQE/gEGAlULQ/dHTrnkcLTgJWmtkHBNcnd8X0u4Bvpqaqua5CoAdQIqkH4eazBjgRmBt/T6u+WuoqhO+AjgAWmlmdme0A/gmcRfrtqzVdqWBmC4BPWiS3VkdnAndbYCHQV1JuvkA6RlfeyKbLzN4xs2weNc4E7jOzBjNbBawguHvaZ9x47GIKMDtuDzKzNQBxPTA1Vc11AUyN3fWZ+R7uMLOPgN8DHxKMxhbgdWBzvAlB+KBzWNq6zOzp+PP1sb5ulNQzn7oIT/cTJfWXVEp4ah5O+u2rNV2QYvtqQWt1lM1VUT7bW1t/uy9IWiLpCUlH5VFTW3RYfbnxAOIY/RnAA2lrSZJF123ASGAs4SZ5Q571VBKeZA4BhgK9gclZDs3r/O9suiR9F7gKOBw4juBz+4p86jKzd4AZwDPAk4Thgx1tnpQH2tCVavtqJ+1yVZQCbwCfNrOjgT8AD6esJ0OH1Zcbj8Bk4A0zWxv312a6wnG9rhB0mdlaM2sys53An9lP3c8cOBlYZWbrzWw78BDwRcLQQeaD0zRcx2TVZWZr4vBGA/AX8l9fmNmdZnasmU0kDDW8RwG0r2y6CqB9JWmtjtJ2VZRVl5ltNbOauP04cKCkAXnU1RodVl9uPALn03xoaB5wYdy+EPh73hUFmulqMbZ7FmH4IZ98CEyQVCpJhPcxbwPPA5kwbWnUVzZd7yT+yUUYm853fSFpYFx/CvgW4e+ZevvKpqsA2leS1upoHnBBnHU1gTBEuSZtXZIGx3aGpPGEe+vGPOpqjXnAFEk9JR1CmBjx6n7JOR+zAwp5Ibxc3Qj0SaT1B54jPCU+B/QrEF33AEuBt2KjGJKCruuA5YQbyz2EWRyHxga5gjDE1rNAdP0j1lcV8FegLAVdLxAM7BLgpAJqX9l0pdK+CAZ1DbCd8KT8/dbqiDAMcwth1tBSEjOLUtY1FVgW63MhoeebT11nxe0GYC3wVOL4q2N9vQtM3l863D2J4ziOkzM+bOU4juPkjBsPx3EcJ2fceDiO4zg548bDcRzHyRk3Ho7jOE7OuPFwHMdxcsaNh1PwSKrJQxnTJZ28l+eOzZf78vYiqSm6Bh8a98sk3SZppaQ3Jb0u6eI95DG/pQtvhdAAt0oaGfPv8L+NU5i48XC6PZKKzGyamT27l1mMJTgXLCTqzWysmWVcUdwBbAJGmdkxwCSCv6+2mE1wzJlkCjDbzFaaWd7joziFgxsPp9MQA+7MlzRX0nJJf4tuKiZLur/FcY/E7dskvRYD5VyXOObfkqZJehE4RyFw1Nnxt2mSFikES7o94XZivqQZkl6V9C9JX4nOK6cD58Un8fNaaL5I0sOSHlEITDVV0k/i0/9CSf3icRfHMpdIejB6u80E+amK6Qti2lFRw+LoAXfUHuptJMFP1TUW/FZhwQ/YjMQxP4vlv5Wop7nA6YreiCWNIDiefDGnP5zTJXHj4XQ2jgEuIwS5ORT4EsE77ARJveMx5wFz4vbVZjYOGAMcL2lMIq9tZvZlM7uvRRl/NLPjzGw0UAKcnvith5mNjxp+aWaNwDRgTnzSn8PujAa+TbiBXw/Uxaf/l4EL4jEPxTKPBt4huJwg5n1qTD8jpl0K3BSf/McR3FK0xVHAkozhaImkUwg+j8YTelGfkzTRzDYS3M5MiodOidfpbikcNx5Op+NVM1sdb4SLgREWYok8CXwjevf9Orsc6Z0r6Q3gTcJN9MhEXtlu9ABfVQiru5QQ6CoZm+GhuH4dGNFOzc+bWbWZrSfEQHkkpi9N5DFa0guxzO8kynwJmBXfTxTFtJeBn0u6guAGvL6dOgCQdHXstWSGtE6Jy5sE1+KHE4wJNB+6ahlbxunGuPFwOhsNie0mQgRBCIbgXMLNfpGZVUcvopcTHP+NAR4DeiXOr22ZuaRewK2E8LWfJbgmT56TKT9Zdi6adyb2dybymAVMjWVelynTzC4FriG41V4sqb+Z3UvohdQDT0k6cQ/lvw0cLemAmOf1sddSkbls4Nex5zTWzD5jZnfG3x4GTlKIFV5iZm+085qdLo4bD6erMJ8Q1/lidvUoKggGYoukQWQPXNWSjKHYIKmMXa7m26IaKM9J7e6UA2skHUjoeQDhfYWZvWJm04ANwHBJhwLvm9nNBO+3Y7LmGDGzFcBrwK8kFcV8e7ErUNBTwPfi9SJpmKLLdgsxKuYDM/Feh5PAjYfTJTCzJuBRgoF4NKYtIQzFLCPc/F5qRz6bCb2NpYSn7kXtKP554MhsL8xz4BfAK4T3N8sT6b+TtFRSFbCA4PL7PKBK0mLCENPd7cj/BwR34iskvQ48S4ysaCFk773Ay3HYbC7NjeFs4Gig5bshpxvjLtkdpwsiqcbMyrpKOU7h4T0Px+mabE1+JLi/yXwkSAg85HRDvOfhOI7j5Iz3PBzHcZyccePhOI7j5IwbD8dxHCdn3Hg4juM4OfM/uQPIYHj/75QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1e79b50e6d8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's limit the fit near to the peak of the histogram.\n",
    "lowerlimit = 70\n",
    "upperlimit = 110\n",
    "bins = 100\n",
    "\n",
    "# Let's select the invariant mass values that are inside the limitations.\n",
    "limitedmasses = invariant_mass[(invariant_mass > lowerlimit) & (invariant_mass < upperlimit)]\n",
    "\n",
    "#Let's create a histogram of the selected values.\n",
    "histogram = plt.hist(limitedmasses, bins=bins, range=(lowerlimit,upperlimit))\n",
    "\n",
    "# In y-axis the number of the events per each bin (can be got from the variable histogram).\n",
    "# In x-axis the centers of the bins.\n",
    "y = histogram[0]\n",
    "x = 0.5*( histogram[1][0:-1] + histogram[1][1:] )\n",
    "\n",
    "# Let's define a function that describes Breit-Wigner distribution for the fit.\n",
    "# E is the energy, gamma is the decay width, M the maximum of the distribution\n",
    "# and a, b and A different parameters that are used for noticing the effect of\n",
    "# the background events for the fit.\n",
    "def breitwigner(E, gamma, M, a, b, A):\n",
    "    return a*E+b+A*((2*np.sqrt(2)*M*gamma*np.sqrt(M**2*(M**2+gamma**2)))/(np.pi*np.sqrt(M**2+np.sqrt(M**2*(M**2+gamma**2)))))/((E**2-M**2)**2+M**2*gamma**2)\n",
    "\n",
    "# Initial values for the optimization in the following order:\n",
    "# gamma (the full width at half maximum (FWHM) of the distribution)\n",
    "# M (the maximum of the distribution)\n",
    "# a (the slope that is used for noticing the effect of the background)\n",
    "# b (the y intercept that is used for noticing the effect of the background)\n",
    "# A (the \"height\" of the Breit-Wigner distribution)\n",
    "initials = [5, 91, -2, 200, 13000]\n",
    "\n",
    "\"\"\"def Exponencial(x, a, b, c):\n",
    "        return a * np.exp(-b * x) + c\n",
    "    \n",
    " ver o link: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html\"\"\"\n",
    "\n",
    "# Let's import the module that is used in the optimization, run the optimization\n",
    "# and calculate the uncertainties of the optimized parameters.\n",
    "from scipy.optimize import curve_fit\n",
    "best, covariance = curve_fit(breitwigner, x, y, p0=initials, sigma=np.sqrt(y))\n",
    "error = np.sqrt(np.diag(covariance))\n",
    "    \n",
    "# Let's print the values and uncertainties that are got from the optimization.\n",
    "print(\"The values and the uncertainties from the optimization\")\n",
    "print(\"\")\n",
    "first = \"The value of the decay width (gamma) = {} +- {}\".format(best[0], error[0])\n",
    "second = \"The value of the maximum of the distribution (M) = {} +- {}\".format(best[1], error[1])\n",
    "third = \"a = {} +- {}\".format(best[2], error[2])\n",
    "fourth = \"b = {} +- {}\".format(best[3], error[3])\n",
    "fifth = \"A = {} +- {}\".format(best[4], error[4])\n",
    "print(first)\n",
    "print(second)\n",
    "print(third)\n",
    "print(fourth)\n",
    "print(fifth)\n",
    "\n",
    "plt.plot(x, breitwigner(x, *best), 'r-', label='gamma = {}, M = {}'.format(best[0], best[1]))\n",
    "plt.xlabel('Invariant mass [GeV]')\n",
    "plt.ylabel('Number of event')\n",
    "plt.title('The Breit-Wigner fit')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notification 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the fitted function does not follow the histogram well, go back and check the intial guesses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notification 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In fitting the so called background of the mass distribution is taken into account. The background basically means muon pairs that come from other decay processes than from the decay of the Z boson. The background is taken into account in the code in the line that follows the command `def breitwigner`. The fit is adapted in the background with the term `a*E+b+A`, where $aE + b$ takes care of the linear part of the background and $A$ the height of the background."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notification 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Even more correct way for doing the fit and getting the values and the uncertainties from it would be to iterate the fit several times. In the iteration a next step would take initial guesses from the previous fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Analysing the histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What can you say about the appearance of the Z boson based on the histogram and the fitted function?\n",
    "\n",
    "Can you define the mass of the Z with the uncertainty? How?\n",
    "\n",
    "Explain your answers with the help from the theory part and other sources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate the lifetime $\\tau$ of the Z boson with the uncertainty by using the fit.\n",
    "\n",
    "Compare the calculated value to the known lifetime of the Z. What do you notice? What could possibly explain your observations?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When was the Z boson detected first time and what is the physical meaning of the Z?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 7"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If energy and momentum could be measured by infinite accuracy, would there be an one exact peak that differs from the other distribution, or an distribution in the histogram on the location of the mass of the Z? Justify your answer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The histogram of the whole data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an example let's also create a histogram from the all of the invariant masses in the data file without limiting near to the peak of the Z boson.\n",
    "\n",
    "Run the code cell below to make that kind of histogram. Notice that the y-axis is logarithmic and the x-axis has logarithms to base 10 of the values of the invariant masses ( $\\log_{10}(\\text{value of the mass})$ ). So for example it is possible to calculate the invariant mass value in units GeV corresponding to the x-axis value of 0.5 with the following way:\n",
    "\n",
    "$$\n",
    "\\log_{10}(\\text{mass}) = 0.5\n",
    "$$\n",
    "\n",
    "$$\n",
    "10^{\\log_{10}(\\text{mass})} = 10^{0.5}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\text{mass} = 10^{0.5} \\approx 3.1622 \\text{GeV}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAElCAYAAAD+wXUWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3Xu4HVV9//H3h0AAuQQlESEQgyRSY1XUI0rhJ2hBgxpARIFaKxRBtBSsthovVShVqxarCKKoNF6DKSImCFIFIVoROdwJAQ0BJIAQBBJA5Pr9/bHWJsNmX2afs+fsvc/5vJ7nPGfPbc2amb3nO2utmTWKCMzMzOqt1+sMmJlZf3KAMDOzhhwgzMysIQcIMzNryAHCzMwacoAwM7OGJmSAkHSspO90Ka0Fkv69xfQHJD2vG+saVJLeI+nOvC+2LDH/IZJ+2aV1f0TS17uR1gjWvUzSHr1Y93gkaUdJV0i6X9LRvc7PRLB+rzNQBUkPFAafATwMPJ6H3z2WeYmITdvNk08i34mIbavP0diStAHweeBVEXFVg+kzgZuADSLisW6vPyI+1e00O1j3C7uRjqRDgHdFxG7dSG+AfRC4MCJe2miipAtJv6OeXBCMR+OyBBERm9b+gN8D8wrjvtvr/I01SZN6uPqtgI2AZT3Mw5iSNC4vvPrAc5lA36O+EBHj+g+4GdizbtyxwCLgW8D9pC/dUGH6NsAPgNWkq9ujW6S/ADgZ+HFO6xJgh8L0AGblz28Arsvz3Qb8M7AJ8BDwBPBA/tsG2BD4AnB7/vsCsGEh3Q8Cd+Rp76pbzwLgFOAc4EFgT+CNwBXAWuBW4NhCWjPz8ofmafcCRwKvAK4G7gNOarEPGuYVeH5ef+TtuqDBsr8vTH8A2AU4BPgl8J85LzcBexeWmQJ8I2//bcC/A5Oa5O1Y0lVlcTvfmdd7N/DRwjF/CHhWYdmX5nk2AHYALgD+mMd9F9ii7nv2oby/HiaVzm8mf/eAnYGL8768AzgJmFz3PTkS+F3e5pMBAS8A/kwqAT8A3NdkOy/M++FXeb4lwJY5n2uBS4GZhfm/mI/1WuAy4P8Vpu0MDOdpdwKfz+M3Ar6T98F9Oc2t2h0TYBZwEbAm77vvt/gu7UP6Pd6Xt+kFefwFeR/8OW/f8+uW+2Td9JOA44Av5ekbkL6Ln83DG+d5n9lqvU3yGMB787G6Hzg+fz8uzvtsUe3Ykr/LDZafVdhv3yKda24BPgasV1yW5r+DQ4CVOQ83AW/v+vmz6hN0r/9oHiD+TDphTwI+Dfw6T1sv/2A+DkwGnpcPwuubpL8AuCf/qNYn/SBPb/JluIP8QwSeCbwsf94DWFWX7r8BvwaeDUwj/fCPz9PmAn8AXkiqQvs2Tw8Qa4Bd8/ZslNfxojz8YtIPf788/8y8/FfyvK/L++esvP7pwF3A7k32Qau81tJev8myT5uev/iPAofn4/MeUuBRnn4W8FVScH028Bvg3U3SP5anB4ivkU4QLyGdzIsnocMLy34O+Er+PAvYixT4pgFLgS/Ufc+uBLYDNq7/7gEvB15F+o7MBJYD76v7npwNbAHMIJ0w5hb2xy8bbV9h+QuBFaQT1RTShchvSRcH65NOQv9dmP9vSQFkfeADpO/TRnnaxcA78udNSdWDkKpnl5C+c5PyNm3e7pgAC4GPsu67uFuTbahdUOxFOqF/MG/T5MI2vqvNPnhXYfi1wDX5818BNwKXFKZdVWa9DdYTwGJgc9Jv8GHgfNK5orbv39ns2PHU3+q3gB8Bm+XvxW+Bw9r9DvJ+XgvsmOfdGnhh18+f3U6w3/5oHiB+VhieAzyUP78S+H3d/B8u/rjqpi0Avl4YfgNwfZMvw+/zj2zzujT24OkB4kbgDYXh1wM358+nAZ8uTJvF0wPEt9rsly8A/5U/z8zLTy9M/yNwYGH4BxROaB3ktZZ2pwFiRWH4GXme55CqrB4mn4Tz9IOBnzdJ/1ieHiC2LUz/DXBQ/vwucikn/whvBV7dJN39gCvqvmd/3+67V5j2PuCHdd+T3QrDi4D5hf1RJkB8tDB8AnBuYXgecGWL5e8FXpI/LyVdfU+tm+fvScH/xXXjWx4T0knw1OJ+b5KHfwUWFYbXI5VG9ihsYycBolZK2BKYD3wEWEUKescBJ5ZZb4P1BLBrYfgy4EN1+/4LzY5dXn4W6aT/MDCnMO3dpHaW2rLNfgebkEo7bynu927/jcs2iJL+UPj8J2CjXHf8XGAbSffV/khfrK06SKtZw/RbSAHkFkkXSdqlRZrbkIqcNbfkcbVptxamFT83HCfplZJ+Lmm1pDWk6oypdcvcWfj8UIPhZtvVKq8j9eQ+jYg/5Y+bko7PBsAdhePzVdJVa8dp89TjdQawi6RtgFeTfoy/AJD0bEmnS7pN0lpSVUv9/mt0HMjLP1/S2ZL+kJf/VIPly36Pmil9/CR9QNJySWvyPpxSyM9hpKvq6yVdKulNefy3gfOA0yXdLumz+SaEdsfkg6SA+5t8Z9ffN8n/U75HEfEEaZ9O72w3PLn8Q6Sqst1Jx/MiUoDbNY+7aBTrHelvpWgqqZai/rdTXG/D30FEPAgcSPod3yHpx5L+osQ6OzKRA0QztwI3RcQWhb/NIuINo004Ii6NiH1JP5yzSFeJkE5E9W4n/fBqZuRxkKqqinc8bddodXXD3yMVi7eLiCmk6iR1tAHNtcprO422vZVbSVddUwvHZ/Powh1DEXEf8L/A24C/ARZGvnQjVUMG6ep5c1IVTf3+a7UtpwDXA7Pz8h9psHzTrJWcrxRJ/4/UXvI2Uh38FqQqSQFExO8i4mDS9/QzwBmSNomIRyPiuIiYQ6qyeRPwd7Q5JhHxh4g4PCK2IV0hf1nSrAZZe8r3SJJI3+3bSm5ao/10Eak66aWkNpOLSCXcnUklpW6st5UHSVf+tbSfU5h2N6kKqf63U2q9EXFeROxFql66nlR12lUOEE/3G2CtpA9J2ljSJEl/KekVo0lU0mRJb5c0JSIeJdUf1m69vRPYUtKUwiILgY9JmiZpKqlNpPbsxiLgUEkvkPSMPK2dzYB7IuLPknYmnQC7pVVe21lNaqAv9axIRNxBOomfIGlzSetJ2kHS7iPJeAPfI5303pI/12xGbiSWNB34lw7T3Yx0zB/IV3rv6WDZO4FtJU3ucJ2t8vIYad+vL+njpPp0ACT9raRp+Ur6vjz6cUmvkfSifFfcWtLJ7fF2x0TSWyXVLmjuJZ3Ia9/9okXAGyX9dS6ZfIAUeH5Vcrvu5Onfo4tIx/O6iHiEXA1Fughc3aX1tnIV8EJJO0naiFTlCUBEPJ7X/UlJm0l6LvB+Svx2JG0laR9Jm+S8PkDjfToqDhB18kGbB+xEujPgbuDrpCL4aL0DuDlXMRxJugolIq4nnWRX5iL6NqS7QIZJd8VcA1yexxER5wInAj8nNaZdnNN/uMW63wv8m6T7SSfwRS3m7VTTvLaTi82fBP4vb/urSiz2d6Si+XWkE84ZpKuoblgMzAbujKc+t3Ec8DLSlfaPgTM7TPefSUH5ftKV3vc7WPYC0h02f5B0d4frbeQ84FxSg+gtpHr6YvXYXGBZfp7oi6Q2mj+T6r7PIAWH5aSTb+1k1uqYvAK4JKe3GDgmIm6qz1RE3ED6TXyJ9LubR7pF/ZGS2/VF4ABJ90o6MY/7FaktolZauC5vb224G+ttKiJ+S7qJ42eku57qHwD9R1IpY2We9j1SG2M765EC2e2km2R2J/3Gu6p2V4gNMEkvAK4l3Qbb9YfNzGxicgliQEl6c662eiapnniJg4OZdZMDxOB6N6kO+UZS3WMnddpmZm25isnMzBpyCcLMzBpygDAzs4YcIMzMrCEHCDMza8gBwszMGnKAMDOzhhwgzMysIQcIMzNryAHCzMwacoAwM7OGHCDMzKwhBwgzM2vIAcLMzBpygDAzs4YcIMzMrCEHCDMza2j9XmdgNKZOnRozZ87sdTbMzAbKZZdddndETGs330AHiJkzZzI8PNzrbJiZDRRJt5SZz1VMZmbWUN+UICStBxwPbA4MR8Q3e5wlM7MJrdIShKTTJN0l6dq68XMl3SBphaT5efS+wHTgUWBVlfkyM7P2qq5iWgDMLY6QNAk4GdgbmAMcLGkOsCNwcUS8H3hPxfkyM7M2Kg0QEbEUuKdu9M7AiohYGRGPAKeTSg+rgHvzPI83S1PSEZKGJQ2vXr26imybmRm9aaSeDtxaGF6Vx50JvF7Sl4ClzRaOiFOB44DLJ0+eXGU+zcwmtF40UqvBuIiIPwGHlUkgIpYAS4aGhg7vas7MzOxJvShBrAK2KwxvC9zeSQKS5kk6dc2aNV3NmJmZrdOLAHEpMFvS9pImAwcBiztJICKWRMQRU6ZMqSSD1iNqVLg0s16p+jbXhcDFwI6SVkk6LCIeA44CzgOWA4siYlmH6boEYWZWMUVEr/MwYkNDQ+GuNsaJWulhgL+PZoNC0mURMdRuvoHsasMlCDOz6g1kgHAbhJlZ9QYyQLgEYWZWvYEMEC5BmJlVbyADhJmZVc8BwszMGhrIAOE2CDOz6g1kgHAbhJlZ9QYyQJiZWfUcIMzMrKGBDBBugzAzq95ABgi3QZiZVW8gA4SZmVXPAcLMzBpygDAzs4YGMkC4kbqO38RmZhUYyAAxIRqpfdI3sx4byABhZmbVc4AwM7OGHCDMzKwhB4h+NJL2B7dZ9J6PgY0zfRMgJO0h6ReSviJpj17nx8xsoqs0QEg6TdJdkq6tGz9X0g2SVkian0cH8ACwEbCqynzVZXLMVmVmNkiqLkEsAOYWR0iaBJwM7A3MAQ6WNAf4RUTsDXwIOK7ifKXAUAsODhJmZk9TaYCIiKXAPXWjdwZWRMTKiHgEOB3YNyKeyNPvBTasMl8N1QJG8a8XiustkwcHNzOryPo9WOd04NbC8CrglZL2B14PbAGc1GxhSUcARwDMmDGjwmyy7uQbUe16zMz6UC8CRKNL3oiIM4Ez2y0cEadKugOYN3ny5Jd3PXeNSClI1P6XXQY6m78bgahMOt1aV7e4FGTWl3pxF9MqYLvC8LbA7Z0k0JOuNurbKzqphmo0X6N02s3Xbv766qlmf63maZYP64z3oY0DvQgQlwKzJW0vaTJwELC4kwR63llfmRNx/bwjCSzN0qj/3Cp/nc5Tn8/6wOFgUp73jQ24qm9zXQhcDOwoaZWkwyLiMeAo4DxgObAoIpZ1ku7AdtbX6oq/35UpBXXawD5eTeRtt3Gl0jaIiDi4yfhzgHNGmq6kecC8WbNmjTQJq0qjINFJe0e/tY90aiSlt9r2Dvq227jTN09Sd2JgSxAT1USujipz+3SjajyzPjCQAaLnbRDWuXZtHiNdvl90O3+9fh7HjAENEC5BDDCf+DrnfWU9MpABwiWIccInvvK8r6wHBjJAuARhE5JLXjbGBjJAmE1oDhI2RgYyQLiKySY8lyZsDAxkgHAVk/UVn6htnBrIAGFmmUsSViEHCOs/PuF1zvvMKjCQAcJtENY3+unE7NKEddlABgi3QUwAPtGNnAOFdUnbACFpB0kb5s97SDpa0hbVZ80s69eTXb/mq6bf82d9r0wJ4gfA45JmAd8Atge+V2muzGCwu0fvF95nNgplAsQT+R0Obwa+EBH/BGxdbbbM+twgnXgHKa/WV8oEiEclHQy8Ezg7j9uguiyZFfjk1h1ul7ARKBMgDgV2AT4ZETdJ2h74TrXZas13MZmNkIOEdaBMgNgrIo6OiIUAEXET8FC12WrNdzFZTw36SXbQ829jpkyAeGeDcYd0OR9m5biqpDu8H62Epu+kzu0OfwNsL2lxYdJmwB+rzphZS83ed131e53H20nV78G2FpoGCOBXwB3AVOCEwvj7gaurzJRZabUTXPHEXX8S7+QE2OiEOd6CQj0HCWuiaYCIiFuAW0gN1Gb9q90JvFHAaBRYyqY3HjlIWANlnqTeX9LvJK2RtFbS/ZLWVpEZSZtIukzSm6pI3wxYFwAmYiBoZbztj/G2PT1QppH6s8A+ETElIjaPiM0iYvMyiUs6TdJdkq6tGz9X0g2SVkiaX5j0IWBR+eybWVf5pGoFZQLEnRGxfITpLwDmFkdImgScDOwNzAEOljRH0p7AdcCdI1yXmXXDeAoS42lbeqBVI3XNsKTvA2cBD9dGRsSZ7RaMiKWSZtaN3hlYERErASSdDuwLbApsQgoaD0k6JyKeKLMRZtZlbpMwygWIzYE/Aa8rjAugbYBoYjpwa2F4FfDKiDgKQNIhwN3NgoOkI4AjAGbMmDHCLJhZWw4SE17bABERh3Z5nY3KfE9+CyNiQZv8nCrpDmDe5MmTX97lvJmZWVbmLqbnSzq/1tAs6cWSPjaKda4CtisMbwvc3kkC7mrDbIz4iesJrUwj9deADwOPAkTE1cBBo1jnpcBsSdtLmpzTWtxmmadwZ31mY8xBYkIqEyCeERG/qRv3WJnEJS0ELgZ2lLRK0mH53RJHAecBy4FFEbGsk0y7BGHWAw4SE06ZRuq7Je1AbieQdACpC462IuLgJuPPAc4pm8l6kuYB82bNmjXSJMzMrI0yJYh/AL4K/IWk24D3AUdWmqs2XIIw65FBKUUMSj77XJkSxC0RsaekTYD1IuL+qjPVjksQZj3k218njDIliJsknQq8Cnig4vyU4hKEWY/5Cn1CKBMgdgR+RqpquknSSZJ2qzZbrfkuJjOz6rUNEBHxUEQsioj9gZeSnqy+qPKctc6TSxBmveZnJMa9MiUIJO0u6cvA5cBGwNsqzZWZmfVc20ZqSTcBV5K64f6XiHiw8ly14UZqsz7iRutxS9HmwEraPCIqeUHQaA0NDcXw8PDIFnbR2Ky7+ilIjOa1sxOApMsiYqjdfGWqmJ7T5b6YzGw88kXXuNOLvpjMbLxykBhXKu2LqSq+zdWsjzlIjBtlAsSI+2Kqim9zNetzvQwSDlBdU6arjX8ATmVdX0w3AW+vNFdmZtZzZd4otxLoq76YzGwA+PbXgVemBAFAPzz/YGYDxkFioJV6ktrMzCaegQwQvovJbIC40XhgtQ0Qkp4h6V8lfS0Pz5b0puqz1pzvYjKzhhyMuqpMCeK/gYeBXfLwKuDfK8uRmY0/7vl1IJUJEDtExGdZ9yT1Q4CPtJnZOFcmQDwiaWPWPSi3A6lEYWbWGZckBkqZ21w/AfwE2E7Sd4FdgUOqzJSZWVf5dtsRKfOg3E8lXU56J7WAYyLi7m5nRNILgGOAqcD5EXFKt9dhZn3CJ+yBUPY2142Ae4G1wBxJry6zkKTTJN1V6yq8MH6upBskrZA0HyAilkfEkaS31bXtp9zMBly3q5tcddV1Zd4o9xngQGAZ8EQeHcDSEukvAE4CvlVIbxJwMrAX6Y6oSyUtjojrJO0DzM/LmJlZD5Vpg9gP2DEiOm6YjoilkmbWjd4ZWJH7eELS6cC+wHURsRhYLOnHwPcapSnpCOAIgBkzZnSaJTPrN92obnLpoRJlAsRKYAO6d+fSdODWwvAq4JWS9gD2BzYEzmm2cEScSupdlqGhIVdimo0HbpPoS00DhKQvkaqS/gRcKel8CkEiIo4e4TobhfqIiAuBC0slIM0D5s2aNWuEWTCzvlMLEg4WfaNVCWI4/78MWFw3bTRHbxWwXWF4W+D2ThKIiCXAkqGhocNHkQ8z6zeuKuorTQNERHwTQNIxEfHF4jRJx4xinZcCsyVtD9xGer/133SSgEsQZvYkB5XKlLnN9Z0Nxh1SJnFJC4GLgR0lrZJ0WEQ8BhwFnAcsBxZFxLKS+TWzicAn/b6gaFLXJ+lg0pX9bsAvCpM2Ax6PiD2rz15rQ0NDMTw83H7GRvwFNBsMzdojRvIbdtsGAJIui4i2z5u1aoP4FXAH6cnmEwrj7weuHl32zMxKKgYCn+DHVKs2iFuAW1jXzXffcBuE2QTlkv+YGsg3yvmFQWZm1RvIAOFXjpqZVa9pgMgPxtX6YuorLkGYmVWvVSP11pJ2B/bJ/SU9pfIvIi6vNGdmZtZTrQLEx0k9q24LfL5uWgCvrSpT7biR2sxGxN14dKTpcxBPziD9a0QcP0b56YifgzCzjjlAdOU5CAAi4vj8nobaS4IujIizR5tBMzPrb23vYpL0adKrQK/Lf8fkcWZmNo6VeR/EG4GdIuIJAEnfBK4APlxlxlpxG4SZWfXKPgexReFzz+8t9W2uZmbVK1OC+DRwhaSfk251fTU9LD2YmdnYKNNIvVDShcArSAHiQxHxh6ozZmZWCd/qWlqZEgQRcQdPf6ucmdlgcpAoZSD7YjIzs+oNZIBwZ31mZtVrGSAkrSfp2rHKTFm+i8nMrHotA0R+9uEqSTPGKD9mZtYnyjRSbw0sk/Qb4MHayIjYp7JcmZlVzQ3VbZUJEMdVngszM+s7ZZ6DuEjSc4HZEfEzSc8AJlWfNTMz66UynfUdDpwBfDWPmg6cVUVmJO0n6WuSfiTpdVWsw8zsSe72v6Uyt7n+A7ArsBYgIn4HPLvsCiSdJumu+ruhJM2VdIOkFZLm57TPiojDgUOAA8uuw8zMuq9MgHg4Ih6pDUhan/RGubIWAHOLIyRNAk4G9gbmAAdLmlOY5WN5upmZ9UiZAHGRpI8AG0vaC/gfYEnZFUTEUuCeutE7AysiYmUOPqcD+yr5DHBus3deSzpC0rCk4dWrV5fNhpmZdahMgJgPrAauAd4NnEO6wh+N6cCtheFVedw/AnsCB0g6stGCEXFqRAxFxNC0adNGmQ0zM2umzF1MT+SXBF1Cqlq6Idq9yLq9Ri1DEREnAie2XdgvDDIzq1zbACHpjcBXgBtJJ/btJb07Is4dxXpXAdsVhrcFbh9FemZmI+MH5poqU8V0AvCaiNgjInYHXgP81yjXeykwW9L2kiYDB9FBd+Lui8nMusq3uzZUJkDcFRErCsMrgbvKrkDSQuBiYEdJqyQdFhGPAUcB5wHLgUURsayDNN2bq5l1l4PE06hZc4Kk/fPHvYDnAotIbRBvJbVDfGBMctjC0NBQDA8Pj2xhfxnMrN4EqWqSdFlEDLWbr1UbxLzC5zuB3fPn1cAzR5G3UXMjtZlZ9ZqWIAaBSxBm1lW18+E4b7juRgmiltD2pOcTZhbnd3ffZmbjW5nuvs8CvkF6evqJarNTjquYzKwS47zk0KkyAeLP+QG2vhERS4AlQ0NDh/c6L2Zm41WZAPFFSZ8A/hd4uDayWV9JZmbjxgQvUZQJEC8C3gG8lnVVTJGHe8JVTGZWGd/A8qS2dzFJuh54cbHL737hu5jMrDIR47YEUfYupjJPUl8FbDH6LJmZ2SApU8W0FXC9pEt5ahuEb3M1MxvHygSIT1Seiw65DcLMrHp+ktrMrJUBPkc2080nqe9n3TuoJwMbAA9GxOajy6KZmfWzMm+U26w4LGk/0julzcxsHCtzF9NTRMRZ9PAZCDOznplgVdNlqpj2LwyuBwyxrsrJzMzGqTJ3MRXfC/EYcDOwbyW5Kcl3MZnZmGr0wNw4fYiuyHcxmZmVUf9k9QAHiFHfxSTp4y2Wi4g4fkQ5MzMbRBPworJVFdODDcZtAhwGbAk4QJjZxDRBgkXTABERJ9Q+S9oMOAY4FDgdOKHZcmZmNj60bKSW9Czg/cDbgW8CL4uIe8ciY2ZmfW+A2yHKaPochKTPAZcC9wMviohjqwwOkp4n6RuSzqhqHWZmVl6rB+U+AGwDfAy4XdLa/He/pLVlEpd0mqS7JF1bN36upBskrZA0HyAiVkbEYSPdEDOznhjH7RFNA0RErBcRG0fEZhGxeeFvsw76YVoAzC2OkDQJOBnYG5gDHCxpzgjzb2Y2tsZxQKjXcVcbnYiIpcA9daN3BlbkEsMjpEbv0g/eSTpC0rCk4dWrV3cxt2ZmVlRpgGhiOnBrYXgVMF3SlpK+ArxU0oebLRwRp0bEUEQMTZs2req8mpm1N05LFWW62ui2RnsyIuKPwJGlEnBXG2ZmletFCWIVsF1heFvg9h7kw8yse8ZhKaIXAeJSYLak7SVNBg4CFneSQEQsiYgjpkyZUkkGzcxGpBYk6v8PqEoDhKSFwMXAjpJWSTosIh4DjgLOA5YDiyJiWYfpzpN06po1a7qfaTMzA9ybq5lZdxV7fe3TJ63L9ubaiyqmUXMJwswGxgBfjA5kgHAbhJn1vQEODDUDGSDMzKx6AxkgXMVkZla9gQwQrmIyM6veQAYIMzOr3kAGCFcxmdlAqW+wHpAG7IEMEK5iMrOBNSDBAQY0QJiZWfUcIMzMrKGBDBBugzAzq95ABgi3QZiZVW8gA4SZ2cBp1zjdh43XDhBmZmOlkyDQBwHDAcLMrJf6IBA04wBhZmYNDWSA8F1MZmbVG8gA4buYzGzg9XHVUs1ABggzM6ueA4SZmTXkAGFm1k3drDrqcTWUA4SZmTW0fq8zUCNpE+DLwCPAhRHx3R5nycxsQqu0BCHpNEl3Sbq2bvxcSTdIWiFpfh69P3BGRBwO7FNlvszM+kKf38lUdRXTAmBucYSkScDJwN7AHOBgSXOAbYFb82yPV5wvMzNro9IAERFLgXvqRu8MrIiIlRHxCHA6sC+wihQkWuZL0hGShiUNr169uopsm5n1Rp+VKHrRSD2ddSUFSIFhOnAm8BZJpwBLmi0cEadGxFBEDE2bNq3anJqZ9VoPg0YvGqkbbW1ExIPAoaUSkOYB82bNmtXVjJmZ9a1aoIgYs1X2ogSxCtiuMLwtcHsP8mFmZi30IkBcCsyWtL2kycBBwOJOEnBfTGZm1av6NteFwMXAjpJWSTosIh4DjgLOA5YDiyJiWYfpujdXM7OKKcawPqvbhoaGYnh4eGQL99ndAmY2gUU8vY2heI5qNH0UJF0WEUPt5hvIrjZcgjAzq95ABgi3QZiZVW8gA4SZ2bjSp1XeAxkgXMVkZuOW1DcBYyADhKuYzGxC6VHAGMgAYWZm1RvIAOEqJjOz6g1kgHAVk5lZ9QYyQJiZWfUcIMzMrKGBDBBugzAzq95ABgi3QZiZVW8gA4SZmVXPAcLMzBoa6O6+Ja0Gbhnh4lOBu7uYnV7ytvSf8bId4G3pV6PZludGxLR2Mw10gBgNScNl+kMfBN6W/jNetgO8Lf1qLLbFVUxmZtaQA4THCsFhAAALCUlEQVSZmTU0kQPEqb3OQBd5W/rPeNkO8Lb0q8q3ZcK2QZiZWWsTuQRhZmYtOECYmVlDEyZASHqWpJ9K+l3+/8wm8z0u6cr8t3is89mKpLmSbpC0QtL8BtM3lPT9PP0SSTPHPpftldiOQyStLhyHd/Uin2VIOk3SXZKubTJdkk7M23q1pJeNdR7LKLEde0haUzgmHx/rPJYlaTtJP5e0XNIyScc0mKfvj0vJ7aj2uETEhPgDPgvMz5/nA59pMt8Dvc5rk3xNAm4EngdMBq4C5tTN817gK/nzQcD3e53vEW7HIcBJvc5rye15NfAy4Nom098AnAsIeBVwSa/zPMLt2AM4u9f5LLktWwMvy583A37b4DvW98el5HZUelwmTAkC2Bf4Zv78TWC/HuZlJHYGVkTEyoh4BDidtE1FxW08A/hrqU/efr5Ome0YGBGxFLinxSz7At+K5NfAFpK2HpvclVdiOwZGRNwREZfnz/cDy4HpdbP1/XEpuR2VmkgBYquIuAPSjgee3WS+jSQNS/q1pH4KItOBWwvDq3j6l+XJeSLiMWANsOWY5K68MtsB8JZc9D9D0nZjk7VKlN3eQbCLpKsknSvphb3OTBm5mvWlwCV1kwbquLTYDqjwuKzfzcR6TdLPgOc0mPTRDpKZERG3S3oecIGkayLixu7kcFQalQTq71EuM0+vlcnjEmBhRDws6UhSqei1leesGoNwTMq4nNR/zwOS3gCcBczucZ5akrQp8APgfRGxtn5yg0X68ri02Y5Kj8u4KkFExJ4R8ZcN/n4E3FkrQub/dzVJ4/b8fyVwISlq94NVQPFKelvg9mbzSFofmEL/VRu03Y6I+GNEPJwHvwa8fIzyVoUyx63vRcTaiHggfz4H2EDS1B5nqylJG5BOqt+NiDMbzDIQx6XddlR9XMZVgGhjMfDO/PmdwI/qZ5D0TEkb5s9TgV2B68Ysh61dCsyWtL2kyaRG6Pq7rIrbeABwQeSWrD7Sdjvq6oL3IdW9DqrFwN/lu2ZeBaypVXUOEknPqbVnSdqZdO74Y29z1VjO5zeA5RHx+Saz9f1xKbMdVR+XcVXF1MZ/AIskHQb8HngrgKQh4MiIeBfwAuCrkp4g7ej/iIi+CBAR8Ziko4DzSHcCnRYRyyT9GzAcEYtJX6ZvS1pBKjkc1LscN1ZyO46WtA/wGGk7DulZhtuQtJB0J8lUSauATwAbAETEV4BzSHfMrAD+BBzam5y2VmI7DgDeI+kx4CHgoD68+KjZFXgHcI2kK/O4jwAzYKCOS5ntqPS4uKsNMzNraCJVMZmZWQccIMzMrCEHCDMza8gBwszMGnKAMDOzhhwgJjBJD4xi2aNyT5hRfDCnVS+ZkraWdHb+PCTpxNFtQal8niNpixEuu5+kOd3OU4d5eJGkBU2m1Xq9/Xoe3qO2f0e4roY9uqpFT8h5H328MPy3+bgvy90/fL3V/s/bsLBu3NS8XRtKOl3S7MK0n0t6IN+ebhVzgLCR+j9gT+CWuvF7kx71nw0cAZxSmPZ+0pPRRMRwRBxdVeZyoFovIt4QEfeNMJn9gJ4GiIi4BthW0owms3w/P8PTDQuAuQ3GzwfOj4jZwPl5uOaDwJchdeMO/BOwd0S8kNQ77K+ArVqs80xgL0nPKIw7AFicn6Y/Ja8DgIh4DTDc2WbZSDlAWO1k+jlJ10q6RtKBefx6kr6crwbPzlfjBwBExBURcXOD5Fr1kvkW4Cc57SevdiUdm69eL5S0UtLRefxnJL23kM9jJX1A0qaSzpd0ec7vvnn6TKW+879M6qNmO0k310o4ks6SdFneniMK6T4g6ZP5ivfXkraS9Fekp7g/p9TP/g51+2yBpFPyFe1KSbvnbVhevOLP8wzndR5XGP8fkq7LV9v/mce9NR+DqyQtLaxuCR0+9Jiv+s/K6f9a0ovz+Gm5FHC5pK9KuqW2f1r06NqwJ2RJzwcejoi787SPAv8cEbfl9B6PiNMi4oY8/8slXZSPwXmSts59Cy0F5hXWdxBQK1X8AthTqesYG2tV9SPuv/7/I7/7gnTi/inpyeatSE+ab026kjuHdCHxHOBe4IC6NG4GphaGzwZ2KwyfDwwB2wOXFcbvQe7HHjiWdKW5ITCV1FXABqR+sC4qLHMd6SnS9YHN87ippKdhBcwEngBe1Sh/wLPy/42Ba4Et83AA8/LnzwIfy58X1G9vId0FpK7KRTqBrgVelPfVZcBOdeucROrb68XAs4AbWPeg6hb5/zXA9OK4/HlXYEmDPBxC4b0Zdfv0S8An8ufXAlfmzycBH86f5+ZtLx6/mdS9EwK4r2743vz/UOCEwvh7gClN9tcG+RhPy8MHkp6ih9SrwQ/z521IfSJNKiz7U+DlheELgaFe/34mwp9LEAawG6n31Mcj4k7gIuAVefz/RMQTEfEH4Ocl0mrWS+bWwOoWy/04ImpXo3eRume/Ani2pG0kvYR0Yvp9XsenJF0N/IzUTXOtGuOWSCWXRo6WdBXwa1JHbbW67UdIgQ3SyX1mie2EdNIO0on9zoi4JiKeAJYV0nibpMuBK4AXkqqs1gJ/Br4uaX9SVw+Qqu0WSDqcFFBq7iKdODuxG/BtgIi4ANhS0pQ8/vQ8/iekoD9STY+pUtvJlZJuzCXSHYG/BH6q1G3Ex0gd5EG+qJC0OfA24IyIeLyQ3Ei237rAAcKg8Um91fhWmvWS+RCwUYvlHi58fpx1/YSdQSrJHEg+sQFvB6aRrip3Au4spP1go8Ql7UFqM9klIl5COmHXlnk0n+jr191OLc9P1OX/CWB9SdsD/wz8dUS8GPgxsFGkd3XsTOqlcz9ytVtEHEk6cW4HXCmp9i6PjUj7rxPNAvVIjmmznpDrj+kyUrsDOVjuRHpr28Z5vcsiYqf896KIeF2e9yHSPngzT61eqhnJ9lsXOEAYpDrgAyVNkjSN9PrJ3wC/JL24Zz1JW5GqMNpp1kvmbyl/ZV50OumkcQApWEDqxvyuiHhU0muA55ZIZwqpBPInSX9Bes1kO/eTXvU4UpuTAtaavP/2hif7958SqXvm9wE75fE7RMQlEfFx4G7WBdrnk6rEOrGUFEhrwfHuSPX9vyRdpSPpdUDDd7PXadYT8nJgVmG+TwP/KWnbwriN8/8bgGmSdsnr3kBPfbnNQtJNDFuRSnhFzycFHxtjDhAG8EPgatL7oS8APpirlH5AKhFcC3yV9DarNQCSjlbq9XNb4GrlWy1JbRYrSe0CXyO9J5uIeBC4UVLxhNJWRCwjnaRvi3XdMX8XGJI0TDoJXl8iqZ+QruqvBo7n6SehRk4H/kXSFfWN1CXzfhWppLIMOI1UhQRpe87OebmIdOcPpAbxa5RuM11KOh4AryGVPjpxLGkfXU3qybh2gj8OeF2u9tobuIMUCGs9ul4M7ChplVLPx+Tl95L0O2CvPEzO40ul1N10DngnAufmBvhfkUpk50V6vewBwGdyNd+VwF8V8vu/pGqk7xdKc+TA+lD0WVfcE4V7c7WWJG0a6W1VW5JKFbvm4DGStN5Mqhb6WFczOY4pvZ/kIlLD/2N10w4hNdYe1WF6j0fqdn0X4JRcFTTS/H2R1Bbzs5Gm0Sb9fwLWRsQ3CuMuJN0t5dtdK+Zbx6yds5UedJoMHD/S4AAQET8s1KtbOTOA+fXBIXsI2FvS16P8sxAzSO9FWY/UOH/4KPP3KeCVo0yjlfvIje2QHpQDngc8WuE6LXMJwszMGnIbhJmZNeQAYWZmDTlAmJlZQw4QZmbWkAOEmZk19P8Be9ha3NXZeegAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "ds = pd.read_csv('DoubleMuRun2011A.csv')\n",
    "invariant_mass_1 = ds['M']\n",
    "\n",
    "no_bins = 300\n",
    "# Let's calculate the logarithms of the masses and weighs.\n",
    "inv_mass_log = np.log10(invariant_mass_1)\n",
    "weights = []\n",
    "for a in invariant_mass_1:\n",
    "    weights.append(no_bins/np.log(10)/a)\n",
    "\n",
    "# Let's plot the weighted histogram.\n",
    "plt.hist(inv_mass_log, bins=no_bins, range=(-0.5,2.5), weights=weights, color=\"red\")\n",
    "plt.yscale('log')\n",
    "\n",
    "# Naming the labels and the title.\n",
    "plt.xlabel('log10(invariant mass) [log10(GeV)]')\n",
    "plt.ylabel('Number of the events')\n",
    "plt.title('The histogram of the invariant masses of two muons \\n')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the histogram that you created to the histogram published by the CMS experiment in the image 10 below. What can you notice? Use the Particle Data Group web site if needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMShistogram.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "    <figcaption>Image 10: The histogram of the invariant masses published by the CMS experiment. &copy; <a href=\"https://arxiv.org/abs/1206.4071\">CMS Collaboration</a> [5]</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Effect of pseudorapidity to the mass distribution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this final section it will be shortly studied how does pseudorapidities of muons that are detected in the CMS detector affect to the mass distribution.\n",
    "\n",
    "As it was told in the theory part, pseudorapidity $\\eta$ describes an angle of which the detected particle has differed from the particle beam (z-axis). Pseudorapidity is determined with the angle $\\theta$ mentioned before with the equation\n",
    "\n",
    "$$\n",
    "\\eta = -\\ln(\\tan(\\frac{\\theta}{2}))\n",
    "$$\n",
    "\n",
    "For recap the image 8 is shown again below. From the image one can see that a small pseudorapidity in practice means that the particle has differed lot from the particle beam. And vice versa: greater pseudorapidity means that the particle has continued almost among the beam line after the collision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMSangles.png\" alt=\"image missing\" style=\"height: 300px\" />\n",
    "    <figcaption>Image 8: Quantities $\\theta$, $\\eta$ and $\\phi$ in the CMS detector.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The image 11 below shows a situation where two particle beams from left and right collide. The image shows two muons with different pseudorapidities. The muon with the smaller pseudorapidity hits the barrel part of the detector when the muon with the greater pseudorapidity goes to the endcap of the detector. There are also muon chambers in the both ends of the detector so these muons can also be detected."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/pseudorapidities.png\" alt=\"image missing\" style=\"height: 300px\" />\n",
    "    <figcaption>Image 11: Two particles with different pseudorapidities in the CMS detector.</figcaption>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this final section it will be studied that how does pseudorapidities of muons that are detected in the CMS detector affect to the mass distribution. For doing that, two different histograms will be made: an one with only muon pairs with small pseudorapidities and an one with great pseduorapidities. The histograms will be made with the familiar method from the earlier part of this exercise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Selecting the events"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let’s create two variables for dividing the events: `small_etas` and `great_etas`. To the first one will be saved only collision events where pseudorapidities of the both detected muons have been small (for example under 0.38). And respectively to the second those whose pseudorapidities have been great (for example over 1.52). Absolute values will be used because $\\eta$ can get also negative values.\n",
    "\n",
    "Complete the code cell below by determining the variables `small_etas` and `great_etas` in a way that the division described above will be made. You will need the following functions:\n",
    "\n",
    "- `ds[condition]` selects from the variable `ds` only events which fulfill the condition written inside the brackets. There can also be more than one condition. Then the function is in the form `ds[(condition1) & (condition2)]`\n",
    "- an example of this could be a function where from the variable `example` only rows where the values of the columns `a` and `b` have been both greater than 8 would be selected: `example[(example.a > 8) & (example.b > 8)]`\n",
    "- you can get the absolute values with the function `np.absolute()` from the _numpy_ module\n",
    "- pseudorapidity of the first muon is `ds.eta1` and the second `ds.eta2`\n",
    "- ”greater than” and ”smaller than” comparisons can be made in Python straight with the symbols > and <\n",
    "- Python uses a dot as a decimal separator (for example 0.38)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    Remember to define the small values in a way that both eta1 and eta2 have been smaller than 0.38. And same for the large values.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    Remember to tell from which variable you want to get the values of the pseudorapidities (write ds.eta1 or ds.eta2). Remember to use \"np.\" in front of the aboslute value function.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 3</summary>\n",
    "    \n",
    "    The first variable with the conditions is \"great_etas = ds[(np.absolute(ds.eta1) > 1.52) & (np.absolute(ds.eta2) > 1.52)]\" and the second \"small_etas = ds[(np.absolute(ds.eta1) < 0.38) & (np.absolute(ds.eta2) < 0.38)]\".\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Amount of all events = 475465\n",
      "Amount of the events where the pseudorapidity of the both muons have been large = 34263\n",
      "Amount of the events where the pseudorapidity of the both muons have been small = 34337\n"
     ]
    }
   ],
   "source": [
    "# Let's import the needed modules.\n",
    "import matplotlib.pyplot as plt\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "# With this line the data is imported and saved to the variable \"ds\".\n",
    "ds = pd.read_csv('DoubleMuRun2011A.csv')\n",
    "\n",
    "great_etas = ds[(np.absolute(ds.eta1) > 1.52) & (np.absolute(ds.eta2) > 1.52)]\n",
    "small_etas = ds[(np.absolute(ds.eta1) < 0.38) & (np.absolute(ds.eta2) < 0.38)]\n",
    "\n",
    "# Let's print out some information about the selection\n",
    "print('Amount of all events = %d' % len(ds))\n",
    "print('Amount of the events where the pseudorapidity of the both muons have been large = %d' %len(great_etas))\n",
    "print('Amount of the events where the pseudorapidity of the both muons have been small = %d' %len(small_etas))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the histograms\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the code cell below to create separate histograms from the events with small and with great values of pseudorapidities. The cell will get the invariant masses for both of the selections and will create the histograms out of them near to the peak that refers to the Z boson."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4cAAANpCAYAAAC4o0xkAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzs3XmYZVV9L/zvTxCcBQUNMqSJUZyjpjXGxBGNAwYI0QQ1EZWEa+LrmJuI8V6nPCYYvTHmzWB4HcBcgxLiGJwIEoka1MaBiEjkAkIDAoogBq/a+Hv/OKdD2VZ31+46Q1fV5/M856naa6+zz4/eXd18e629VnV3AAAAWNtuNu8CAAAAmD/hEAAAAOEQAAAA4RAAAIAIhwAAAEQ4BAAAIMIhAAAAEQ4BAACIcAgAAECSXeddwLTttddevW7dunmXAQAAMBdnn332N7p77+31W/XhcN26ddmwYcO8ywAAAJiLqvraUvqZVgoAAIBwCAAAgHAIAABAhEMAAAAiHAIAABDhEAAAgAiHAAAARDgEAAAgwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMmu8y4AAPhR6449davnLj7ukBlWAsBaYuQQAAAA4RAAAADhEAAAgAiHAAAARDgEAAAgwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMIhAAAAEQ4BAACIcAgAAECEQwAAACIcAgAAEOEQAACACIcAAABkzuGwqt5aVVdV1ZcWtL2uqr5SVedU1Xuqao8F515aVRdU1flV9bj5VA0AALD6LDscVtXNq+oBVXXQDrz9hCSP36LttCT36e77JfmPJC8df869khyZ5N7j9/x1Ve2yw4UDAADwX5YcDqvq16rq5Kq6w4K2uyY5N8mGJF+uqndX1a5LvWZ3n5nkmi3aPtrdm8aHZyXZb/z9YUne2d3f6+6LklyQ5MFL/SwAAAC2bsjI4bOT3KO7F4a5/5Xkp5OckeScjALcsyZXXp6d5EPj7/dNcumCcxvHbQAAACzTkHB4rySf3XxQVbdL8sQkJ3f3YzIaxftKJhQOq+plSTYlecfmpkW69Vbee0xVbaiqDVdfffUkygEAAFjVhoTDvZNcseD455PsmuSdSdLdP8joecG7LreoqjoqyZOSPL27NwfAjUn2X9BtvySXL/b+7j6+u9d39/q99957ueUAAACsekPC4fVJbr/g+BEZjdx9YkHb/01y2+UUVFWPT/KSJId29w0LTr0/yZFVtXtVHZjkbkk+s5zPAgAAYGTJi8ck+WqSJ1TV7hmFwqckOae7v7Ggz08muWqpF6yqk5I8MsleVbUxySsyWp109ySnVVWSnNXdz+nuc6vq5CRfzmi66XO7+8YB9QMAALAVQ8Lh8UnellFI/EGSdUleuEWfn8to9dIl6e6nLtL8lm30f02S1yz1+gAAACzNkqeVdveJSY5LcquMppf+5fiVJKmqR2cUGM+YbIkAAABM25CRw3T3Hyb5w62c/kSSPZP853KLAgAAYLaWPHJYVc+oqvtt7Xx3fz+j1USfNonCAAAAmJ0hq5WekOTw7fQ5LKPnEgEAAFhBhoTDpdglW9mYHgAAgJ3XpMPh3ZN8a8LXBAAAYMq2uSBNVb11i6bDq2rdIl13SXJAkoclOXUilQEAADAz21ut9JkLvu8k9x+/FtNJPp3kRcsvCwAAgFnaXjg8cPy1klyY5M+TvHGRfjcm+VZ328YCAABgBdpmOOzur23+vqpeleSMhW0AAACsDtsbOfwv3f2qaRYCAADA/Cw5HG5WVbskOSjJnhktRPNjuvvMZdYFAADADA0Kh1X1PzNacOb22+m6aGgEAABg57TkcFhVf5DkVUmuS/J3SS5NsmlKdQEAADBDQ0YOfzvJZUke2N1XT6keAAAA5uBmA/run+S9giEAAMDqMyQcXpkdWMAGAACAnd+QcHhyksdW1e7TKgYAAID5GBIOX57kiiSnVNWBU6oHAACAORgyTfTcJDdPcpckT6yq65Jcu0i/7u67TqI4AAAAZmNIOLxZRltXXLKgrRbpt1gbAAAAO7Elh8PuXjfFOgAAAJijIc8cAgAAsErtcDisqj2rav9JFgMAAMB8DAqHVXWbqvpfVfX1JN9IctGCcz9XVR+sqgdOukgAAACma8nhsKpun+TfkrwoyeVJzsuPLj7z70keluSpkywQAACA6RsycviyJPdO8szufmCSf1h4srtvSPLxJAdPrjwAAABmYUg4PCLJR7r77dvo87Uk+y6vJAAAAGZtyD6H+yX5x+30+U6S2+94OQCwMq079tStnrv4uENmWAkA7JghI4fXJ7nTdvocmNFCNQAAAKwgQ0YOP5vkSVV12+6+fsuTVbVPkicm+adJFQcAq8HWRhWNKAKwMxkycvjGJHdM8sGquufCE+Pjf0hyiyR/MbnyAAAAmIUljxx290eq6pVJXpnkS0l+kCRV9Y0ke2a0rcVLuvtTky8TAACAaRoycpjufnVGW1W8P8m3ktyYpJN8MMljuvt1E68QAACAqRvyzGGSpLvPSHLGFGoBAABgTpY8clhVe0yzEAAAAOZnyMjhFVX1/iQnJvlwd/9wSjUBwJqwrb0RAWDWhjxzeHGSpyT5QJLLqup1VXXfqVQFAADATC05HHb3PZP8XJI3Jbl5kt9L8oWqOruqnl9Ve02pRgAAAKZs6Gqln+3u5ybZJ6NRxFOT3DfJn2c0mvjeqjp88mUCAAAwTYPC4Wbd/YPu/sfuPjTJXZK8OKO9Dw9NcsoE6wMAAGAGBm9lsYhvJjk3yXlJ7jOhawLATskiMgCsVjsc5KrqHkmOSvIbGY0eVpILMlrNFAAAgBVkUDisqj2TPDWjULg+o0D47SRvSXJCd39q4hUCAAAwdUsOh1V1SpInZbRSaSf55yQnJHlPd//fqVQHAADATAwZOTwiyfkZTRv9u+6+bDolAQAAMGtDVit9aHffs7uPm1QwrKq3VtVVVfWlBW13qKrTquqr4697jturqv6iqi6oqnOq6oGTqAEAAIAB4bC7z5rC55+Q5PFbtB2b5PTuvluS08fHSfKEJHcbv45J8jdTqAcAAGBNGrzPYVX9clW9s6q+WFUXLGi/Z1X9QVXtu9RrdfeZSa7Zovmw3LTi6YlJDl/Q/vYeOSvJHlW1z9D6AQAA+HFDFqSpjEb6fmPc9N0kt1zQ5VtJ/jijFUxfu4ya7tzdVyRJd19RVXcat++b5NIF/TaO265YpNZjMhpdzAEHHLCMUgAAANaGISOHv5vkN5O8Lckdkrx+4cnu/nqSTyY5ZGLV/ahapK0X69jdx3f3+u5ev/fee0+pHAAAgNVjSDg8OskXk/x2d1+XxYPZV5McuMyartw8XXT89apx+8Yk+y/ot1+Sy5f5WQAAAGRYODwoyRndveho3dhVSZY7VPf+JEeNvz8qyfsWtD9jvGrpQ5Jct3n6KQAAAMszZJ/DTUlusZ0++yb5zlIvWFUnJXlkkr2qamOSVyQ5LsnJVXV0kkuSPGXc/YNJnpjkgiQ3JHnWgNoBAADYhiHh8MtJHllVtdjoYVXdIsmjk3x+qRfs7qdu5dTBi/TtJM9d6rUBAABYuiHTSv8uyT2SvKGqfuR9VbVLkj9LcpeMVjQFAABgBRkycvi3SQ5N8vyMpnpenyRVdUqSh2QUDN/X3e+YdJEAAABM15JHDrv7xiRPSvLqJLsluXtG20sckeRWSf4oNz0fCAAAwAoyZOQw3b0pySur6lUZhcM7JrkuyVfG4REAAIAVaFA43Gy8OMz5E64FAACAORmyIA0AAACrlHAIAACAcAgAAIBwCAAAQIRDAAAAso1wWFV/VlW/tOD4gKq63WzKAgAAYJa2tZXFC5Ncm+Sj4+OLkrwyo83uAYA5WHfsqYu2X3zcITOuBIDVZlvTSr+T5FYLjmv8AgAAYJXZ1sjhBUmOqKr3JLli3LZHVR2wvYt29yWTKA4AAIDZ2FY4fF2S/53kUwvaXjB+bUtv57oAAADsZLYa4rr7pKq6KMkhSfZN8swk5yT5wmxKAwAAYFa2OcLX3WclOStJquqZSd7T3a+eQV0AAADM0JDpn89K8vlpFQIAAMD8LDkcdveJ0ywEAACA+dnWVhaLqqojq+qfq+qbVbWpqq6pqtOq6shpFAgAAMD0LXnksKoqyduTPC2j/Q5vTHJ1kr2SHJzk0VV1aHc/bRqFAgAAMD1DRg7/W5KnJ/lcksckuUV375PkFuPjs5P8elU9Z+JVAgAAMFVDFqR5dpKLkzy8u7+7ubG7b0zysap6RJIvJTk6yZsmWSQAsGPWHXvqVs9dfNwhM6wEgJ3dkJHDe2W0lcV3Fzs5bn9vkntOojAAAABmZ0g47IyeNdyW7Z0HAABgJzQkHJ6X5IiquuViJ8fthyf58iQKAwAAYHaGhMO3JjkgyZlVdXBV7ZokVbVLVT0qyRlJfnLcDwAAgBVkyII0f5vkYUmemuSjSX5YVdckuUNGIbOSnNzdFqMBgBnb1sIzALAUSx457JGnZ7SdxceSXJdRMLxufPz07j5yKlUCAAAwVUNGDpMk3X1SkpOmUAsAAABzMuSZQwAAAFYp4RAAAADhEAAAAOEQAACACIcAAABEOAQAACADwmFVPaOqHjfNYgAAAJiPISOHb03y+GkVAgAAwPwMCYdfH9gfAACAFWJI2PtwkkdVlYAIAACwygwJei9Lctskb6mqvaZUDwAAAHOw64C+JyW5LskzkhxZVRdnNNW0t+jX3X3wZMoDAABgFoaEw0cu+H73JAeNX1vaMiwCAACwk1tyOOxuzxoCAACsUgIfAAAAO284rKoXVdW5VfWlqjqpqm5RVQdW1aer6qtV9a6q2m3edQIAAKwGg8JhVd2sqp5XVWdV1XVVtWnBuQdU1V9X1d2XW1RV7Zvk+UnWd/d9kuyS5Mgkr03yhu6+W5JvJTl6uZ8FAADAgHA4HqU7LcmfJ7lrkuuT1IIuFyV5dpKnT6i2XZPcsqp2TXKrJFckeXSSU8bnT0xy+IQ+CwAAYE0bMnL4+0keleRVSe6c5M0LT3b3tUnOTPK45RbV3ZcleX2SSzIKhdclOTvJtd29ebRyY5J9l/tZAAAADAuHT0/yye5+dXf/MItvWXFRkgOWW1RV7ZnksCQHJrlLklsnecIiXRfdNqOqjqmqDVW14eqrr15uOQAAAKvekHB4YJKzttPnmiR32PFy/stjklzU3Vd39w+SvDvJQ5PsMZ5mmiT7Jbl8sTd39/Hdvb671++9994TKAcAAGB1GxIOv5tkj+30OSDJtTtezn+5JMlDqupWVVVJDk7y5SRnJHnyuM9RSd43gc8CAABY84aEwy8k+aWtbR9RVbfP6HnDzyy3qO7+dEYLz3wuyb+P6zw+yUuSvLiqLkhyxyRvWe5nAQAAMFoRdKn+vyTvSPKOqvqRLSSqao8kb0uyZ5I3TaKw7n5Fklds0XxhkgdP4voAAADcZMnhsLtPqqrHJHlWkkMz2mcwVbUhyb2T7J7kr7r7g9MoFAAAgOkZMq003X10RnsZfjnJ3hntc/jAJBckObq7nzfxCgEAAJi6IdNKkyTdfUKSE6rqlhlNI72uu/9z0oUBAAAwO4PD4Wbd/d2MVjAFAABghRscDqvqNkl+JckDktw+yXVJPp/kPd39ncmWBwAAwCwMCodV9ZSMViPdI6PnDTfrJH9eVf+tu0+ZYH0AAADMwJLDYVU9NslJSX6Y5O1J/iXJ15P8RJJHJXlakpOq6tru/ufJlwoAAMC0DBk5fHmS7yV5WHd/botzJ1bVXyY5c9xPOAQAAFhBhmxl8YAk71okGCZJuntDkpMz2toCAACAFWRIOPxekiu20+fycT8AAABWkCHh8F+T/OJ2+vxCRlNLAQAAWEGGhMOXJLlvVR1XVbdeeKKqbl1Vf5rkPkmOnWSBAAAATN9WF6Spqrcu0nxOkt9PckxVfS7JlUnunNFzhrfPaNTwD5IcPflSAYBJWnfsqYu2X3zcITOuBICdwbZWK33mNs7tkeTRi7Q/IsnDIxwCAACsKNsKhwfOrAoAAADmaqvhsLu/NstCAAAAmJ8hC9IAAACwSm1rWumiqupmSfZNsl+Smy/Wp7ttZwEAALCCDAqHVfX7Sf57kr2203WXHa4IAACAmVtyOKyqVyZ5eZJvJjkxyWVJNk2nLAAAAGZpyMjh0UkuTPKz3X3dlOoBAABgDoYsSHPHJO8XDAEAAFafIeHwgiR7TqsQAAAA5mdIOPzrJE+qqp+YVjEAAADMx5KfOezuN1XV3ZN8sqpeneRzSRadYtrdl0yoPgAAAGZg6D6HX0zyzCRv3Uaf3oHrAgAAMEdDtrL4rSR/m9H2Ff+S5PLYygIAAGBVGDLC93tJrkry0O6+aEr1AAAAMAdDwuG6JG8WDAFYC9Yde+q8SwCAmRqyWullSW4+rUIAAACYnyHh8O1JnlhVt51WMQAAAMzHkGmlf5zkfkn+uapekuTs7r5+OmUBwPSZOgoANxkSDr83/lpJTk+SqlqsX3e3rSwAAABWkCEh7l8z2sMQAACAVWbJ4bC7HznFOgAAAJijIQvSAAAAsEoJhwAAACx9WmlVvXyJXbu7/2gH6wEAAGAOhixI88ptnNu8UE2NvxcOAQAAVpAh4fBRW2nfI8mDkjw/yalJ3rTcogAAAJitIauVfnwbp99XVe9K8pkk71x2VQAAAMzUxBak6e5/T/K+JH84qWsCAAAwG5NerfSSJPeZ8DUBAACYskmHw59L8t0JXxMAAIApG7KVxQHbuMb+SX47yS8mOXkCdQEAADBDQ1YrvTg3bVmxmEry1ST/fTkFAQAAMHtDwuHbs3g4/GGSb2W0Uun7uvt7kyisqvZI8uaMnmHsJM9Ocn6SdyVZl1FY/bXu/tYkPg8AAGAtG7KVxTOnWMdi3pjkw9395KraLcmtMloJ9fTuPq6qjk1ybJKXzLguAACAVWfSC9JMRFXdLsnDk7wlSbr7+919bZLDkpw47nZiksPnUyEAAMDqslOGwyQ/leTqJG+rqs9X1Zur6tZJ7tzdVyTJ+Oud5lkkAADAajHkmcNU1d2SvCDJg5PsmWSXRbp1d991AnU9MMnzuvvTVfXGjKaQLrXOY5IckyQHHLC1RVYBAADYbMkjh1X180m+kOR3k9w/yS0yWqF0y9ckRiM3JtnY3Z8eH5+SUVi8sqr2GdezT5KrFntzdx/f3eu7e/3ee+89gXIAAABWtyEjh3+SZPckz0ny1u7eNJ2Sku7+elVdWlUHdff5SQ5O8uXx66gkx42/vm9aNQAAAKwlQ8Lhg5Kc0t3HT6uYLTwvyTvGK5VemORZGY1KnlxVRye5JMlTZlQLAADAqjYkHH4/o0A2E939hSTrFzl18KxqAAAAWCuGPB/4qSQPmFYhAAAAzM+QcPiHSR5aVb85rWIAAACYjyHTSg9L8rEkJ1TVbyU5O8m1i/Tr7v6jSRQHAADAbAwJh69c8P3Dxq/FdBLhEAAAYAUZEg4fNbUqAAAAmKslh8Pu/vg0CwEAAGB+hixIAwAAwColHAIAACAcAgAAMGxBGgBYkdYde+q8SwCAnZ6RQwAAAIRDAAAAhEMAAAAiHAIAAJAJLkhTVS9PsinJyd19waSuCwALbWtxmYuPO2SGlQDA6jLJkcNXJvmjJOdV1UlVdb8JXhsAAIApmuRWFq/KKGzeN8ljkvxakl0meH0AAACmZGLhsLtftfC4qu4zqWsDAAAwXVNbkKa7vzStawMAADBZEwmHVbVnVd16EtcCAABg9pYcDqvq4Kr606rac0Hbnarq40m+keSaqvqzaRQJAADAdA0ZOXxekiO6+1sL2l6f5GFJLkjyzSQvqKpfm2B9AAAAzMCQcPgzST6x+aCqbpnkyUlO6+6DkhyU5NIkz5lohQAAAEzdkHB4pySXLzj+uSS3SHJCknT39Un+KaOQCAAAwAoyJBx+L8ktFxw/LEknOXNB27eT3GECdQEAADBDQ8LhRUkeveD4V5N8tbsvW9C2f0aL0wAAALCCDAmHJya5b1V9uqr+Ncl9k/z9Fn0emOT8SRUHAADAbOw6oO/fJHlIkl9PUkk+kOS1m09W1YOT3DPJSZMsEAAAgOlbcjjs7h8keVpVPWd02Ndv0eXCJA9IcvHkygMAAGAWhowcJkm6+9tbaf9GPG8IAACwIg155hAAAIBVatDIYVXdOsnvJnlckn2T7L5It+7uu06gNgAAAGZkyeGwqvZI8okk98poP8PbJbkuyW65af/Dy5P8YMI1AgAAMGVDppX+j4yC4dFJ9hy3vSHJbZI8NMnnkvyfjFYsBQAAYAUZEg4PTXJmd7+tu3tzY4+cleSJSe6R5GUTrhEAAIApGxIO989odHCzH2bBM4fdfVWSDyU5cjKlAQAAMCtDFqS5IcmNC46vS/ITW/S5MqOFagBg5tYde+q8SwCAFWvIyOGlGY0ebvblJA+vql0WtP1ikq9PojAAAABmZ0g4/HiSR1RVjY/fleSuSU6tqudW1T8keUiSD064RgAAAKZsyLTSEzPatmK/jEYR35Tk0UkOT/JL4z6fzGhVUwBghdra9NyLjztkxpUAMEtLDofd/bkkv7PgeFOSI6rqZ5P8dJKLk3y2u3846SIBAACYriEjh4vq7rOTnD2BWgAAAJiTJT9zWFUXVtXzt9PnuVV14fLLAgAAYJaGLEizLske2+mzR5Kf3OFqAAAAmIsh4XApbpPk+xO+JgAAAFO2zWcOq+qALZr2WKQtSXZJckCSJycxrRQAAGCF2d6CNBcn6QXHLxi/tqaSvHiZNd10sapdkmxIcll3P6mqDkzyziR3SPK5JL/Z3UYqAQAAlml74fDtGYXDSvKMJOck+cIi/W5M8s0kp3f3RydY3wuSnJfkduPj1yZ5Q3e/s6relOToJH8zwc8DAABYk7YZDrv7mZu/r6pnJHlPd7962kWNP2+/JIckeU2SF1dVJXl0kqeNu5yY5JURDgEAAJZtyfscdvekF6/Znj9P8gdJbjs+vmOSa7t70/h4Y5J9Z1wTAADAqrRDga+q7lFVv1JVvznpgsbXf1KSq7r77IXNi3TtRdpSVcdU1Yaq2nD11VdPo0QAAIBVZVA4rKr7V9WGJOcmOSXJCQvOPaKqbqiqX55AXb+Q5NCqujijBWgendFI4h5VtXm0c78kly/25u4+vrvXd/f6vffeewLlAAAArG5LDodVdfck/5LkoCRvTPKhLbqcmeSajLazWJbufml379fd65IcmeRj3f30JGcsuP5RSd633M8CAABg2MjhK5LsluTB3f3iJJ9deLK7O8m/JXnQ5Mr7MS/JaHGaCzJ6BvEtU/wsAACANWPJC9IkOTjJu7v7vG30uSTJY5dX0o/q7n/JaMQy3X1hkgdP8voAAAAMGzncI6MVQrd3vd12vBwAAADmYUg4vCrJT2+nz72TXLrj5QAAADAPQ8Lhx5L8clUdtNjJqnpQRlNPPzKJwgAAAJidIeHwT5JsSnJmVf1OkrskSVXde3z8gSTXJ3n9xKsEAABgqpa8IE13n19Vv5rkpCR/OW6uJOeMv16b5IjuvmTiVQIAADBVQ1YrTXd/uKoOzGiPwYdktJ3EdUnOSvK27r5m8iUCAAAwbYPCYZJ097VJ3jh+AQBk3bGnLtp+8XGHzLgSAHbUkGcOAQAAWKUGjxwCAGvT1kYHAVgdjBwCAAAgHAIAAGBaKQA7KVMYAWC2jBwCAACw9HBYVQ+vqgO202f/qnr48ssCAABgloaMHJ6R5Jnb6fOMcT8AAABWkCHhsJbYp3ewFgAAAOZk0s8cHpDk+glfEwAAgCnb5mqlVfXyLZoeWbXoAOIuGQXDI5N8YjKlAQAAMCvb28rilQu+7ySPHL+25rIkxy6rIgAAAGZue+HwUeOvleRjSU5IcuIi/W5M8s0k53f3DydWHQAAADOxzXDY3R/f/H1VnZjkvQvbAAAAWB22N3L4X7r7WdMsBAAAgPlZcjgEAJiUdceeutVzFx93yAwrAWCzQVtZVNU+VfVXVXVBVX23qm5c5LVpWsUCAAAwHUseOayqfZN8Jsmdk5ybZPckX0vyvSQ/Nb7WF5JcN/kyAQAAmKYhI4cvT/ITSR7f3T8zbntbd98jo3D4kSS3THLEZEsEAABg2oaEw8cl+XB3//OWJ7p7Y5KnZBQOXzWh2gAAAJiRIeHwJzKaTrrZjRmFwSRJd38nyWlJDptMaQAAAMzKkHD47SS7LTj+VpJ9t+hzXZK9l1sUAAAAszUkHH4tyf4Ljr+Y5NFVdaskqaqbJfmlJBsnVx4AAACzMCQcnp7kUVV18/HxiUnukuRTVfW6JJ9Mcu8k75psiQAAAEzbkreySPKWjKaS7pXkiu7+31X1s0mel+R+4z7vTPKayZYIAADAtC05HHb3V5O8dou2F1XVH2e0lcXF3X3lhOsDAABgBoaMHC6qu69OcvUEagFgDVp37KnzLgEAyA6Gw6q6R5J7JrlNd//dZEsCAABg1oYsSJOqun9Vbchov8NTkpyw4NwjquqGqvrlyZYIAADAtC155LCq7p7kX5LskuSNSe6e5AkLupyZ5JokT07ygcmVCMBKsbUpohcfd8iMKwEAhhoycviKJLsleXB3vzjJZxee7O5O8m9JHjS58gAAAJiFIeHw4CTv7u7zttHnkoz2PgQAAGAFGRIO90iycQnX223HywEAAGAehoTDq5L89Hb63DvJpTteDgAAAPMwJBx+LMkvV9VBi52sqgdlNPX0I5MoDAAAgNkZEg7/JMmmJGdW1e9k/GxhVd17fPyBJNcnef3EqwQAAGCqlryVRXefX1W/muSkJH85bq4k54y/XpvkiO6+ZOJVAgAAMFVLDodJ0t0frqoDkxyV5CFJ7pjkuiRnJXlbd18z+RIB2NlsbT9DAGDlGhQOk6S7r03yxvELAACAVWDJzxxW1Y1V9ffTLGbBZ+1fVWdU1XlVdW5VvWDcfoeqOq2qvjr+uucs6gEAAFjthixIc32Sr02rkC1sSvJ73X3PjKavPreq7pXk2CSnd/fdkpw+PgYAAGCZhoTDzye517QKWai7r+juz42/vz7JeUk0lFzHAAAgAElEQVT2TXJYkhPH3U5Mcvgs6gEAAFjthjxz+NokH6iqx3b3adMqaEtVtS7JA5J8Osmdu/uKZBQgq+pOW3nPMUmOSZIDDjhgNoUCAD9mRxYv2tp7Lj7ukOWWA8A2DAmHd0ry4SQfqqr3Jvlskq8n6S07dvfbJ1FcVd0myT8meWF3f7uqlvS+7j4+yfFJsn79+h+rDwAAgB81JByekFEQrCRHjF/Jj4bDGh8vOxxW1c0zCobv6O53j5uvrKp9xqOG+yS5armfAwAAwLBw+KypVbGFGg0RviXJed39ZwtOvT+jPRaPG39936xqAgAAWM2WHA67+8Tt95qYX0jym0n+vaq+MG77w4xC4clVdXSSS5I8ZYY1AQAArFpDRg5nprs/kdEU1cUcPMtaAAAA1oIhW1kAAACwSi155LCqLlxCtx8m+XZG+xK+u7v/cUcLAwAAYHaGTCu92bj/XcbHm5J8M8kdF1zn8oy2vLh/kiOr6oNJDu/uGydTLgAAANMwZFrp/ZJcluRfk/xiklt09z5JbpHkYeP2jUn2TXJQRnsiPjHJCyZZMAAAAJM3JBy+Jsntkxzc3Z/q7h8mSXf/sLs/meSxSfZI8pru/mpGK4leluTpE64ZAACACRsyrfRXkvx9d29a7GR3f7+qPpDkqUme3903VNXpSZ48gToBWKZ1x566aPvFxx0y40pgx/g9DDBdQ0YO75hkt+30ufm432Zfz066XQYAAAA3GRIOL0zyq1V128VOVtXtkvxqkosWNO+T5JodLw8AAIBZGDKqd3ySNyT5dFW9Jsknk1yZ5M4ZLVDzsoxWMn1xklRVJXlkki9MsF4AMpvpdVv7DABgdVpyOOzuN1bVQUmek+Tti3SpJMd39xvHx3dKclKS05ZdJQAAAFM16HnA7v7dqvr7JM/MaC/D22e06f3nk7y9u89c0PfKJC+dXKkAAABMy+DFYrr7E0k+MYVaAFimnXUq6M5aFwBwkyEL0gAAALBK2WYCYEbs0QY7v22NcvtZBVY7I4cAAAAIhwAAAJhWCqxQpn5NjsViAIDEyCEAAAARDgEAAIhppbDTWE3TJK3KCQCw8hg5BAAAQDgEAADAtNJVYTVNR2Qy/J6YnB2ZImv1T9h5mOYOsHRGDgEAADByuJIYjZi+Sf4ar4V/lV7Lvycn+d++ln8dYSUxCgmsdkYOAQAAEA4BAAAwrZRVwDSflWNHpk+6jwAAs2HkEAAAAOEQAAAA00pXvaHT+CY5hW/eUwjnvQKk6a4AO69Z/R3h7wJgJTFyCAAAgHAIAACAaaVMgM3Ah9mR/8ZJTktaaVOcZvXrtRZ+7wGr00r7cx3YeRk5BAAAwMjhvOysoxQ7a11M31r4l2e/v2F1Wok/27OY4bCzLgy3mv5egdXGyCEAAADCIQAAAKaVwoo2i6lUK3G6FsAs+XNymG39eplyCvNl5BAAAADhEAAAANNKAQAYW2n7xM5qiupKW3l1pdU7b369bmLkEAAAAOEQAAAA00oBAHYaVj5d3I78ugx9z846hXCSU2d35NdxZ/11YTpW5MhhVT2+qs6vqguq6th51wMAALDSVXfPu4ZBqmqXJP+R5LFJNib5bJKndveXF+u/fv363rBhwwwrXBr/MggAwM5uayOHK3G/ykmOJq+0RWyq6uzuXr+9fitx5PDBSS7o7gu7+/tJ3pnksDnXBAAAsKKtxHC4b5JLFxxvHLcBAACwg1bigjS1SNuPzI2tqmOSHDM+/E5VnT/1qobbK8k35l0Ec+Her13u/drl3q9d7v3atGrue712Nu/ZGe3gf/vOeu9/cimdVmI43Jhk/wXH+yW5fGGH7j4+yfGzLGqoqtqwlHm/rD7u/drl3q9d7v3a5d6vTe772rXS7/1KnFb62SR3q6oDq2q3JEcmef+cawIAAFjRVtzIYXdvqqr/J8lHkuyS5K3dfe6cywIAAFjRVlw4TJLu/mCSD867jmXaqae9MlXu/drl3q9d7v3a5d6vTe772rWi7/2K2+cQAACAyVuJzxwCAAAwYcLhjFTVHlV1SlV9parOq6qfr6o7VNVpVfXV8dc9510nk1VVB1XVFxa8vl1VL3TvV7+qelFVnVtVX6qqk6rqFuOFtD49vu/vGi+qxSpTVS8Y3/dzq+qF4zY/86tQVb21qq6qqi8taFv0XtfIX1TVBVV1TlU9cH6Vs1xbufdPGf/c/7Cq1m/R/6Xje39+VT1u9hUzKVu5968b/z/+OVX1nqraY8G5FXXvhcPZeWOSD3f3PZL8TJLzkhyb5PTuvluS08fHrCLdfX5337+775/kZ5PckOQ9ce9XtaraN8nzk6zv7vtktHjWkUlem+QN4/v+rSRHz69KpqGq7pPkt5M8OKM/659UVXeLn/nV6oQkj9+ibWv3+glJ7jZ+HZPkb2ZUI9NxQn783n8pyRFJzlzYWFX3yujvgHuP3/PXVbXLDGpkOk7Ij9/705Lcp7vvl+Q/krw0WZn3Xjicgaq6XZKHJ3lLknT397v72iSHJTlx3O3EJIfPp0Jm5OAk/6e7vxb3fi3YNcktq2rXJLdKckWSRyc5ZXzefV+d7pnkrO6+obs3Jfl4kl+Jn/lVqbvPTHLNFs1bu9eHJXl7j5yVZI+q2mc2lTJpi9377j6vu89fpPthSd7Z3d/r7ouSXJDRPyCxAm3l3n90/Gd+kpyV0T7syQq898LhbPxUkquTvK2qPl9Vb66qWye5c3dfkSTjr3eaZ5FM3ZFJThp/796vYt19WZLXJ7kko1B4XZKzk1y74C+PjUn2nU+FTNGXkjy8qu5YVbdK8sQk+8fP/FqytXu9b5JLF/TzZ8Da4d6vLc9O8qHx9yvu3guHs7Frkgcm+ZvufkCS/4wpRWvK+NmyQ5P8w7xrYfrGzxgdluTAJHdJcuuMppRtyXLRq0x3n5fR9OHTknw4yReTbNrmm1grapE2fwasDe79GlFVL8voz/x3bG5apNtOfe+Fw9nYmGRjd396fHxKRmHxys1TSsZfr5pTfUzfE5J8rruvHB+796vbY5Jc1N1Xd/cPkrw7yUMzmka2eX/Z/ZJcPq8CmZ7ufkt3P7C7H57R1KOvxs/8WrK1e70xo1HkzfwZsHa492tAVR2V5ElJnt437RW44u69cDgD3f31JJdW1UHjpoOTfDnJ+5McNW47Ksn75lAes/HU3DSlNHHvV7tLkjykqm5VVZWbfubPSPLkcR/3fZWqqjuNvx6Q0eIUJ8XP/FqytXv9/iTPGK9a+pAk122efsqq9/4kR1bV7lV1YEaLEn1mzjUxQVX1+CQvSXJod9+w4NSKu/d1U7Blmqrq/knenGS3JBcmeVZG4fzkJAdk9D+TT+nuLR9sZ4UbP3d0aZKf6u7rxm13jHu/qlXVq5L8ekbTSz6f5Lcyes7gnUnuMG77je7+3tyKZCqq6l+T3DHJD5K8uLtP9zO/OlXVSUkemWSvJFcmeUWS92aRez3+h6K/zGjFwhuSPKu7N8yjbpZvK/f+miT/b5K9k1yb5Avd/bhx/5dl9CzapiQv7O4PLXJZVoCt3PuXJtk9yTfH3c7q7ueM+6+oey8cAgAAYFopAAAAwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMIhAAAAEQ4BAACIcAgAAECEQwAAACIcAgAAEOEQAACACIcAAABEOAQAACDCIQAAABEOAQAAiHAIAABAhEMAAAAiHAIAABDhEAAAgAiHAAAARDgEAAAgwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMIhAAAAEQ4BAACIcAgAAECEQwAAACIcAgAAEOEQAACACIcAAABEOAQAACDCIQAAABEOAQAAiHAIAABAhEMAAAAiHAIAABDhEAAAgAiHAAAARDgEAAAgwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMIhAAAASXaddwHTttdee/W6devmXQYAAMBcnH322d/o7r2312/Vh8N169Zlw4YN8y4DAABgLqrqa0vpZ1opAAAAwiEAAADCIQAAAJlzOKyqt1bVVVX1pQVtr6uqr1TVOVX1nqraY8G5l1bVBVV1flU9bj5VAwAArD7zHjk8Icnjt2g7Lcl9uvt+Sf4jyUuTpKruleTIJPcev+evq2qX2ZUKAACwes01HHb3mUmu2aLto929aXx4VpL9xt8fluSd3f297r4oyQVJHjyzYgEAAFaxeY8cbs+zk3xo/P2+SS5dcG7juA0AAIBl2mnDYVW9LMmmJO/Y3LRIt97Ke4+pqg1VteHqq6+eVokAAACrxk4ZDqvqqCRPSvL07t4cADcm2X9Bt/2SXL7Y+7v7+O5e393r99577+kWCwAAsArsdOGwqh6f5CVJDu3uGxacen+SI6tq96o6MMndknxmHjUCAACsNrvO88Or6qQkj0yyV1VtTPKKjFYn3T3JaVWVJGd193O6+9yqOjnJlzOabvrc7r5xPpUDAACsLnXTrM3Vaf369b1hw4Z5lwEAADAXVXV2d6/fXr+dblopAAAAsyccAgAAIBwCAAAgHAIAAJA5r1YKAGvZumNPXbT94uMOmXElAGDkEAAAgAiHAAAARDgEAAAgwiEAAAARDgEAAIhwCAAAQIRDAAAAIhwCAAAQ4RAAAIAIhwAAAEQ4BAAAIMIhAAAAEQ4BAACIcAgAAECEQwAAACIcAgAAEOEQAACACIcAAABEOAQAACDCIQAAABEOAQAAiHAIAABABoTDqlpXVU+sqlsvaNu1ql5VVV+sqk9V1a9Mp0wAAACmadcBfV+R5NAkd17Q9j+S/M8FxydX1cO6+6xJFAcAa9G6Y08d/J6LjztkCpUAsJYMmVb680lO7+5NSVJVN0vyu0m+kuSAJA9O8p9JXjTpIgEAAJiuIeHwzkm+tuD4/kn2SvJX3b2xuzckeV+SB02wPgAAAGZgSDi8eZJecPwL4+OPLWjbmGSfCdQFAADADA0JhxuT3G/B8ROTfKO7z1vQdqck355EYQAAAMzOkAVp/inJi6rq9Un+b5LHJnnbFn3ukR+degoAAMAKMCQc/mmSw5O8eHx8WUYrmCZJquonkzw0yRsmVh0AAAAzseRw2N1XVdV9kxw8bvp4d1+/oMttMgqOH5lgfQAAAMzAkJHDdPd3M5peuti5/7+9Ow+zpCzvPv79MSwqLiAgIoiDhrgr6EhUREBMQFEhxgVXICAv5nVfIsYoKC5oYhQTNRlFB6KCiigEUIIgSvRlxwVFRGWEQUBUFhdcRu73j6oOh+Z0T1f3Waa7v5/rOledeuo5VXd19TL3PNt3ge8OIihJkiRJ0mjNODlM8ifgsKo6fJo6bwLeWlWdkk5Jkhay2SxqL0nSqHWZrTTtayb1JEmSJEnzSJfkcCY2ppnJVJIkSZI0j0zb/TPJEycVLe1TBrAE2Bp4AXDZgGKTJEmSJI3ImsYGngVU+76AfdtXPwFuBV47kMgkSZIkSSOzpuTwbTRJYYC30CSLX+1T70/AL4CvVNX3BxmgJEmSJGn4pk0Oq+qwifdJ9gW+UFUfGHZQkiRJkqTRmvGSE1W1zTADkSRJkiSNz6BnK5UkSZIkzUOdFqtPsi3wSmAHmmUrlvSpVlX1gAHEJkmSJEkakRknh0keB3wZuDOwGriu3d6h6mBCkyRJkiSNSpeWw3cBGwAHAx+rqn6JoSRJkiRpHuqSHD4GOL6qlg8rGEmSJEnSeHSZkOYPwJXDCkSSJEmSND5dksNvANsP8uJJPpbkZ0ku6Sm7Z5LTk1zebjduy5PkA0l+mOTbSR41yFgkSZIkaTHrkhz+A/D4JC8a4PVXAHtMKjsEOKOqtgXOaPcBngJs274OAj48wDgkSZIkaVHrMuZwL+BMYEWSA4ELgRv71KuqOnwmJ6yqryVZ2uc6u7TvjwbOAt7Qlh9TVQWck2SjJFtU1TUd7kGSJEmS1EeX5PCwnvc7ta9+CphRcjiFzScSvqq6Jsm92vItgat66q1qy0wOJUmSJGmOuiSHuw4tipnpt35i9a2YHETT9ZStt956mDFJkiRJ0oIw4+Swqr46zEB6XDfRXTTJFsDP2vJVwH176m0F/LTfCdrlNpYDLFu2rG8CKUmSJEm6TZcJaUblJGDf9v2+wIk95S9uZy19LHCT4w0lSZIkaTC6dCsFIMkjgOcDDwY2rKont+VLgR2A06vqhhme61iayWc2TbIKOBQ4AvhMkgNo1lV8dlv9VOCpwA+B3wL7d41dkiRJktRfp+QwydtolrSYaHHs7bK5DnAs8CrgX2dyvqp63hSHdutTt4D/O+NgJUmSJEkzNuNupUn2Af4ROB3YDnhX7/Gq+jFwAfCMQQYoSZIkSRq+LmMOX0HTpXOvqvo28Ic+dS6lWaRekiRJkjSPdEkOHw6cVlX9ksIJPwU2n1tIkiRJkqRR65IcBrh1DXU2B343+3AkSZIkSePQJTm8HHj8VAeTLAGeAHx3rkFJkiRJkkarS3L4GeBRSV47xfE3An8GfGrOUUmSJEmSRqrLUhbvp1lz8D1JnkO7jEWSfwZ2ApYB5wDLBx2kJEmSJGm4ZpwcVtUtSXYFjgReACxpD72GZiziJ4CXVdXqgUcpSZIkSRqqLi2HVNVNwH5JXgM8BtgEuAk4r6quH0J8kiRJkqQR6JQcTqiqXwKnDTgWSZIkSdKYzHhCmiSfTvKUJF0msZEkSZIkzQNdEr1nAycDVyf5pyQPH1JMkiRJkqQR65IcPg74D2B94LXAN5NckOTlSTYdSnSSJEmSpJGYcXJYVedW1d8BWwDPAU4FHkEze+nVSU5IsneSWY1jlCRJkiSNT+fxg1X1h6o6vqqeDmxJ04p4KbA38Dngp4MNUZIkSZI0bHOaXKaqrq+q9wHbA68DVtMsbyFJkiRJmkfm1AU0yQOBfYEX0rQiBrh8AHFJkiRJkkaoc3KYZGNgH5qk8DE0CeHNwFHA0VX19YFGKEmSJEkauhknh0meRpMQPo1mxtICvgwcDZxQVb8bSoSSJEmSpKHr0nJ4Urv9AU1CeExVXT34kCRJkiRJo9YlOVwOrKiqc4YVjCRJkiRpPGacHFbVwcMMRJIkSZI0PrOZkGY9YDfgwcBdq+rwtvxOwN2Bn1fVrQONUpIkSZI0VJ3WOUyyB7ASOAV4L3BYz+HtgGuA5w4oNkmSJEnSiMw4OUyyDPgCzSylrwY+1Xu8HYt4BfDXgwxQkiRJkjR8XVoO3wz8FlhWVR+g/2L35wOPHERgkiRJkqTR6ZIc7gh8oaqunabOVcAWcwtJkiRJkjRqXZLDuwI/X0Odu3Q8pyRJkiRpLdAlkbsaeOga6mwH/Hj24UiSJEmSxqFLcvhFYPckT+h3MMlTgMcDJw8iMEmSJEnS6HRJDt8F3Aj8d5J3Aw8BSLJnu/9ZmqUs/mXgUUqSJEmShmrdmVasqquT/BXwGeD1PYdOAgL8CHhmVa1pXKIkSZIkaS0z4+QQoKouSvJAYE/gccAmwE3AOcCJVbV68CFKkiRJkoatU3IIUFV/omktPGnw4UiSJEmSxsFlJyRJkiRJJoeSJEmSJJNDSZIkSRImh5IkSZIkTA4lSZIkSUyTHCZ5RJJ7jTIYSZIkSdJ4TNdyeDFw8MROkjOTvHj4IUmSJEmSRm265PBWYEnP/i7A0mEGI0mSJEkaj+mSw1XAdqMKRJIkSZI0PutOc+y/gJcluRS4pi3bL8kuazhnVdVugwhOkiRJkjQa0yWHbwLWB/YEdgaKplvp0jWcswYRmCRJkiRpdKbsVlpVv6qqg6vqvlW1BAhwWFWts4bXkqnOKUmSJElaO3VZ5/CrwMohxSFJkiRJGqPpupXeTlXtOsxAJEmSJEnjM+PkcEKSuwDPBLYHNgJuAi4CPl9VvxlseJIkSZKkUeiUHCZ5KnA0cE+aMYgTCnhfkv2r6uQBxidJkiRJGoEZJ4dJHgWcACwBPgmcSbPExRbAk4DnAccn2bGqLhxCrJIkSZKkIenScvgmmhbCnarqnEnHViT5IHAW8A/A3wwmPEmSJEnSKHSZrXQn4LN9EkMAqupc4Pi23pwleXWS7ya5JMmxSe6UZJsk5ya5PMmnk6w/iGtJkiRJ0mLXJTm8B3DVGupcCdx99uE0kmwJvAJYVlUPo+nKug/wbuB9VbUtcANwwFyvJUmSJEnqlhz+FNhhDXWW0YxDHIR1gTsnWRe4S3veJ9G0TkIzMc7eA7qWJEmSJC1qXZLDU4EnJTkkyZLeA0nWSfJa4MltvTmpqquBf6ZpibyGZrmMC4Ebq2p1W20VsGW/zyc5KMkFSS64/vrr5xqOJEmSJC14XZLDw4FrgXcAP0xyTJJ3JzkauBx4T3v87XMNKsnGwF7ANsB9gA2Bp/SpWv0+X1XLq2pZVS3bbLPN5hqOJEmSJC14M56ttKquTbIj8B/AXwL3m1TldODgqhpEt9InA1dU1fUASU4AHg9slGTdtvVwK5qurpIkSZKkOeqylAVVtRLYvZ0wZnuaSWpuAi5uu4IOypXAY5PcBbgF2A24APgK8CzgOGBf4MQBXlOSJEmSFq1OyeGENhEcZDI4+fznJjkeuAhYDVwMLAdOAY5L8va27KhhxSBJkiRJi8msksNRqKpDgUMnFf+YNc+YKkmSJEnqqMuENJIkSZKkBcrkUJIkSZJkcihJkiRJWovHHEqSpJlbesgpfctXHrHniCORJM1XthxKkiRJkmaeHCb5U5JPDTMYSZIkSdJ4dGk5/BXwk2EFIkmSJEkany5jDi8GHjKsQCRJms+mGvMnSdJ80aXl8N3AU5P85bCCkSRJkiSNR5eWw3sBXwK+mOQLwPnAtUBNrlhVxwwmPEmSJEnSKHRJDlfQJIIBntm+4PbJYdp9k0NJkiRJmke6JIf7Dy0KSZIkSdJYzTg5rKqjhxmIJEmSJGl8ukxII0mSJElaoLp0KwUgyWbA3wAPBjasqgN7yrcBvlNVtww0SkmSJEnSUHVKDpMcAHwAuBO3TT5zYHt4c+D/AQcBRw0wRkmSJEnSkM24W2m7vuFy4AfAXwMf7j1eVZcA3wX2HmSAkiRJkqTh69Jy+AbgGmDnqro5yfZ96nwbeNxAIpMkSZIkjUyXCWmWASdX1c3T1FkF3HtuIUmSJEmSRq1Ly+H6wG/WUGcj4E+zD0eSJA3S0kNOmfLYyiP2HGEkkqS1XZeWw5XAo9dQ5y+Ay2YdjSRJkiRpLLokhycCOyV5dr+DSfYHHgF8bhCBSZIkSZJGp0u30vcA+wDHJnkWcA+AJC8DdgKeCVwO/Ougg5QkSZIkDdeMk8OquiHJzsAxQG/r4Qfa7dnA86tqTeMSJUmSJElrmS4th1TVlcAuSR5Bs2TFJsBNwDlVdeEQ4pMkSZIkjUCn5HBCVX2bZk1DSZIkSdICMKvkMMl6wINpxh3eBFxaVX8cZGCSJEmSpNHpMlspSTZJ8hHgRuBi4Kx2e2OSjyTZdPAhSpIkSZKGbcYth0k2B74O3J+mtfA84Frg3sB2wAHArkl2rKrrhhCrJEljN92i8pIkzWddWg7fSZMYvh+4X1XtWlXPq6pdgfsBR7bH3zH4MCVJkiRJw9RlzOHTgLOr6jWTD1TVzcCrkywDnj6o4CRJkiRJo9Gl5fBuwP+soc7ZwF1nH44kSZIkaRy6JIffB7ZYQ50tgMtmH44kSZIkaRy6JIdHAs9N8oh+B5NsBzyHZkyiJEmSJGkemXLMYZInTiq6AjgdOC/JMcDXgOuAzYGdgRcBXwRWDiVSSZIkSdLQTDchzVlA9SkPcCDN0hW9ZQB7Ac8AlgwiOEmSJEnSaEyXHL6N/smhJEmSJGmBmTI5rKrDRhiHJEmSJGmMukxII0mSJElaoEwOJUmSJEnTjjm8gyRbAa8GtgO2AtbrU62q6gEDiE2SJEmSNCIzTg6T7AKcCtwJWE2zjMXqflUHEpkkSZIkaWS6tBy+h2aJihcDn6qqW4cTkiRJkiRp1Lokhw8Hjq2qTwwrGEmSJEnSeHSZkOYG4JfDCkSSJEmSND5dksOTgZ2HFYgkSZIkaXy6JIf/ANwjyQeTbDisgCRJkiRJozfjMYdV9fMkewDnAi9O8gPgpv5Va7dBBShJkiRJGr4uS1k8FPgKsHFbtP0UVWuuQUmSJEmSRqtLt9J/ATYB3gLcD1ivqtbp81oylEglSZIkSUPTZSmLxwEnVNXbhxWMJEmSJGk8urQc/gFYOaQ47iDJRkmOT/L9JJcmeVySeyY5Pcnl7XbjNZ9JkiRJkrQmXZLDs4AdhhRHP0cCX6qqBwGPBC4FDgHOqKptgTPafUmSJEnSHHVJDv8eeEiSQ5JkWAEBJLk78ETgKICq+kNV3QjsBRzdVjsa2HuYcUiSJEnSYtFlzOE/ApcA7wBekuSbTL2UxQFzjOv+wPXAx5M8ErgQeCWweVVd017kmiT3muN1JEmSJEl0Sw7363m/Tfvqp4C5JofrAo8CXl5V5yY5kg5dSJMcBBwEsPXWW88xFEmSJEla+Lokh1Mlg8OwClhVVee2+8fTJIfXJdmibTXcAvhZvw9X1XJgOcCyZctcd1GSJEmS1mDGyWFV/WSYgUy61rVJrkrywKq6DNgN+F772hc4ot2eOKqYJEmSJGkh69JyOGovBz6ZZH3gx8D+NBPofCbJAcCVwLPHGJ8kSZIkLRgzTg6TzHjwXlVdObtwbneObwLL+hzaba7nliRJkiTdXpeWw5U0k82sSXU8ryRJkiRpzLokccfQPzncCNgOuB9wFjCysYmSJEmSpMHoMiHNflMdS7IO8GbgYJqJYiRJkiRJ88g6gzhJVd1aVW+l6Xp6xCDOKUmSJEkanYEkhz2+AfzVgM8pSZIkSRqyQSeH9wQ2HPA5JUmSJElDNrDkMMmTgecClwzqnJIkSZKk0eiyzuGZ05zjvsDEOohvmxglhnYAAB5/SURBVGtQkiRJkqTR6rKUxS5TlBdwA3Aa8M9VNVUSKUmSJElaS3VZymLQ4xMlSZIkSWsJEz5JkiRJksmhJEmSJKnbmEMAkiwDdgA2Bpb0qVJVdfhcA5MkSZIkjU6X2UrvDpwA7ApkmqoFmBxKkiRJ0jzSpeXwn4AnAWcDHweuAlYPIyhJkiRJ0mh1SQ73Ai4Cdq2qW4cUjyRJkiRpDLpMSHMP4CsmhpIkSZK08HRJDi8HNh9WIJIkSZKk8emSHH4QeHqSLYcVjCRJkiRpPLqMOfwizYQ0X0/yVuBC4MZ+FavqygHEJkmSJEkakS7J4UqaZSoCfHSaetXxvJIkSZKkMeuSxB1Dk/hJkiRJkhaYGSeHVbXfEOOQJEmSJI1RlwlpJEmSJEkLlMmhJEmSJMnkUJIkSZJkcihJkiRJwuRQkiRJkoTJoSRJkiQJk0NJkiRJEiaHkiRJkiQGmBwm2TrJfQZ1PkmSJEnS6Ayy5XAlcFWSrybZfYDnlSRJkiQN2SCTwyuBVcCOwKlJLhzguSVJkiRJQ7TuoE5UVUsBkmwEPLF9SZIkSZLmgYElhxOq6kbgpPYlSZIkSZoHnK1UkiRJkjTzlsMkGwNbAD+qqt/3lO8P7A38Bnh/VZ038CglSRqxpYecMu4QJEkaqS7dSt8JvBC410RBkpcD7wfSFu2dZFlVfW9wIUqSJEmShq1Lt9IdgTOq6paestcBV9NMPvOctuw1A4pNkiRJkjQiXVoOtwTOmNhJ8hDgvsAbqup/2rJn4yylkiRJkjTvdGk5vDPwu579HYECvtxT9iOaJFKSJEmSNI90SQ6vBh7Us787cDPwrZ6yjYHebqeSJEmSpHmgS7fSrwD7JnkZTQviM4DPVdWtPXX+DLhqgPFJkiRJkkagS8vhu4BfA0cCy2kSxMMmDia5F7Az8I0BxidJkiRJGoEZtxxW1RVJHgo8qy06qaqu7KlyP+CDwKcGGJ8kSZIkaQS6dCulqq4F/m2KY+cD5w8iKEmSJEnSaHXpVipJkiRJWqA6tRwCJNkC2I1myYoN+lSpqjp8roFJkiRJkkanU3KY5K3AIZM+F5r1DnvfmxxKkiRJ0jwy426lSV4AvBk4m2ZSmgBHA88HPgLcChwHPGnwYUqSJEmShqnLmMOXAquAParq823Zyqo6rqoOBp4GPAe4+6CCS7IkycVJTm73t0lybpLLk3w6yfqDupYkSZIkLWZdksOHA6dW1eqesiUTb6rqNOA04PUDig3glcClPfvvBt5XVdsCNwAHDPBakiRJkrRodUkO1wN+0bN/C3CPSXUuAR4516AAkmwF7Al8tN0PTZfV49sqRwN7D+JakiRJkrTYdUkOrwG26Nm/EnjEpDpbAqsZjPcDf08zlhFgE+DGnpbLVe31JEmSJElz1CU5vJima+mEM4GdkrwoyYZJ9gT+pq03J0meBvysqi7sLe5TtfqUkeSgJBckueD666+faziSJEmStOB1SQ5PBh6aZJt2/wjgJmAFcDNwEk0C948DiGtH4BlJVnLbDKjvBzZKMrGMxlbAT/t9uKqWV9Wyqlq22WabDSAcSZIkSVrYZpwcVtWKqrpLVV3R7l8FPAb4MPDfwHLgMVV1zlyDqqo3VtVWVbUU2Ac4s6peAHyFZhkNgH2BE+d6LUmSJEnS7Rez76xNFF82oFhm4g3AcUneTtN99agRXluSpAVl6SGn9C1fecSeI45EkrQ2mHFymOQtwFlV9bVp6uwE7FpVbxtEcABVdRZwVvv+x8AOgzq3JEmSJKnRZczhYcAua6jzRODQ2QYjSZIkSRqPLsnhTKzLbUtPSJIkSZLmiUEnh48Gfj7gc0qSJEmShmzaMYdJzpxUtF+SXfpUXQLcF7gfcOxgQpMkSZIkjcqaJqTZped9AUvb12S3Ar8APg28egBxSZIkSZJGaNrksKr+t9tpkluBwwY5E6kkSZIkae3QZZ3D/WnWFpQkSZIkLTAzTg6r6uhhBiJJkiRJGp9Os5UmWSfJy5Ock+SmJKt7jm2f5ENJ/nzwYUqSJEmShmnGyWGS9YHTgfcDDwB+BaSnyhXA3wIvGGSAkiRJkqTh6zLm8PXArsBhwNuBtwBvnjhYVTcm+RqwO3DoAGOUJGlOlh5ySt/ylUfsOeJIJElae3XpVvoC4OtV9baqupVmaYvJrgC2HkhkkiRJkqSR6ZIcbgOcs4Y6vwTuOftwJEmSJEnj0KVb6S3ARmuoszVw4+zDkSRp9qbqPipJktasS8vhN4G/aiemuYMk96AZb3jeIAKTJEmSJI1Ol+TwI8B9gU8muXvvgSQbASuAjYF/H1h0kiRJkqSRmHG30qo6NsmTgf2BZwA3ACS5AHgosAHwwao6dRiBSpI0aHZD7Wa6r5czv0rS/Nel5ZCqOoBmLcPvAZvRrHP4KOCHwAFV9fKBRyhJkiRJGrouE9IAUFUrgBVJ7kzTjfSmqvrNoAOTJEmSJI1O5+RwQlXdQjODqSRJkiRpnuvUrVSSJEmStDCZHEqSJEmSTA4lSZIkSSaHkiRJkiRMDiVJkiRJmBxKkiRJkjA5lCRJkiQxzTqHSd4yy3NWVR0+y89KkiRJksZgyuQQOGyW5yzA5FCSJEmS5pHpksNdRxaFJEmSJGmspkwOq+qrowxEkiRJkjQ+TkgjSZIkSTI5lCRJkiRNP1vprTSTy3RVVTXdWEZJkiRJ0lpmuiTua8wuOZQkSZIkzTPTTUizywjjkCRJkiSNkWMOJUmSJEkmh5IkSZKk6ccc9pVkC2A3YEtggz5VqqoOn2tgkiRJkqTR6ZQcJnkrcMikz4XbJq6ZeG9yKEmSJEnzyIy7lSZ5AfBm4GzgWTSJ4NHA84GPALcCxwFPGnyYkiRJkqRh6tJy+FJgFbBHVa1OArCyqo4DjkvyeeAU4NjBhylJkiRJGqYuE9I8HDi1qlb3lC2ZeFNVpwGnAa8fUGySJEmSpBHp0nK4HvCLnv1bgHtMqnMJcPBcg5IkSeOz9JBTxh2CJGkMurQcXgNs0bN/JfCISXW2BFYjSZIkSZpXuiSHF9N0LZ1wJrBTkhcl2TDJnsDftPUkSZIkSfNIl26lJwMfSrJNVV0BHAE8F1jRvgD+CPzjIAOUJElrv6m6oq48Ys8RRyJJmq0ZJ4dVtYLbkkCq6qokjwFeCzwAWAl8qKq+M9gQJUmSJEnD1qXl8A7aFsSXDSgWSZIkSdKYzCk5lCRpJqab/XKqbod2U5QkabQ6J4dJ1qGZlXQrmuUt7qCqvjbHuCRJkiRJI9QpOUzyeuB1wKZrqLpk1hFJkiRJkkZuxslhksOAtwC/AI4GrmZIaxomuS9wDHBv4FZgeVUdmeSewKeBpTQT4Dynqm4YRgySJEmStJh0aTk8APgx8OiqumlI8UxYDby2qi5KcjfgwiSnA/sBZ1TVEUkOAQ4B3jDkWCRJkiRpwVunQ91NgJNGkBhSVddU1UXt+18Bl9KMc9yLptWSdrv3sGORJEmSpMWgS8vhD4GNhxXIVJIsBbYHzgU2r6proEkgk9xris8cBBwEsPXWW48mUEnSSEw386kkSZq9Li2HHwKeluTewwpmsiR3BT4HvKqqbp7p56pqeVUtq6plm2222fAClCRJkqQFYsYth1X170n+HPh6krcBFwF9u5hW1ZVzDSzJejSJ4Ser6oS2+LokW7SthlsAP5vrdSRJkiRJ3dc5/BbNpDAfm6ZOzeK8t5MkwFHApVX1Lz2HTgL2BY5otyfO5TqSJEmSpEaXpSwOBP6DZibRs4CfMqSlLIAdgRcB30nyzbbsH2iSws8kOQC4Enj2kK4vSZIkSYtKlxa+19J043x8VV0xpHgAqKr/ATLF4d2GeW1JkiRJWoy6JIdLgY8OOzGUJC0uzj4qSdLaoctspVcD6w0rEEmSJEnS+HRJDo8BnprkbsMKRpIkSZI0Hl2Sw3cC5wFfTrKLSaIkSZIkLRxdxhz+vt0GOAOgWXHiDqqq5rSUhSRJkiRptLokcWfTrGEoSZIkSVpgZpwcVtUuQ4xDkiRJkjRGXcYcSpIkSZIWqDmPDUyyKbAT8Fvgy1X1pzlHJUmSJEkaqRm3HCZ5aZJzk9yzp+zRwKXA8cCpwDeSbDj4MCVJkiRJw9SlW+lzaWYi/WVP2T8BGwMfp0kOHwMcPLjwJEmSJEmj0CU53Bb49sRO2510Z+Coqjqwqp4OnA88f7AhSpIkSZKGrUtyuAnws579Hdvt53vKzgbuN9egJEmSJEmj1SU5/CWwac/+zsCtwDd6ygq40wDikiRJkiSNUJfk8FLg6Uk2SbIRzRjE86vq5p46S4FrBxifJEmSJGkEuiSHRwJbAKuAq4B7Ax+aOJhkCfAE4FuDDFCSJEmSNHwzXuewqk5KcjBwUFv0yar6RE+VJ9N0KT1tgPFJkqQFaOkhp0x5bOURe44wEknShBknhwBVtRxYPsWx02iWtZAkSZIkzTNdupVKkiRJkhYok0NJkiRJksmhJEmSJMnkUJIkSZJExwlpJEmLh7NJSpK0uNhyKEmSJEkyOZQkSZIk2a1U0lpuqq6No+rWOJvrjzvmrqbrPjrOc2lh8HtCkuYPWw4lSZIkSSaHkiRJkiS7lUqSpHnOmXUlaTBsOZQkSZIkmRxKkiRJkuxWKmkt4YyGkiRJ42XLoSRJkiTJ5FCSJEmSZLdSSVqrTdXddjYzMLrYvSRJmo4th5IkSZIkk0NJkiRJkt1KF7xBdkmT5rtxd6scxfX92dZCYLfl/vy5lzRsthxKkiRJkkwOJUmSJEl2K9UCYDebO5quS9ZC+boshnucDbvjSbfX9WdiMf/+kCRbDiVJkiRJJoeSJEmSJLuVzivj7j45isW4x92dp2tc87Fr49r6tR8ku1ZKmq35+Ht9vvFrLK29bDmUJEmSJJkcSpIkSZLsVrpWWsxd4lwkfHxm8/VazN+rkqTGKP7e+jddGg1bDiVJkiRJJoeSJEmSJLuVjs0ouk+Oyriv39W4453N9WczW6okabBG8Tt3Ps5APkj+XdM4rK0/D+Ngy6EkSZIkaX4mh0n2SHJZkh8mOWTc8UiSJEnSfDfvupUmWQJ8EPhLYBVwfpKTqup7441sfrH5fH4ZZzcbu/hI0ugspN+5o7iX6a7RdUjGbP4NNIqhIoOcNXzc9zgbg/w+6hrXuO99HOZjy+EOwA+r6sdV9QfgOGCvMcckSZIkSfPafEwOtwSu6tlf1ZZJkiRJkmYpVTXuGDpJ8mxg96o6sN1/EbBDVb28p85BwEHt7gOBy0Ye6JptCvx83EFoLHz2i5fPfvHy2S9ePvvFyee+eK2tz/5+VbXZmirNuzGHNC2F9+3Z3wr4aW+FqloOLB9lUF0luaCqlo07Do2ez37x8tkvXj77xctnvzj53Bev+f7s52O30vOBbZNsk2R9YB/gpDHHJEmSJEnz2rxrOayq1UleBpwGLAE+VlXfHXNYkiRJkjSvzbvkEKCqTgVOHXccc7RWd3vVUPnsFy+f/eLls1+8fPaLk8998ZrXz37eTUgjSZIkSRq8+TjmUJIkSZI0YCaHI5JkoyTHJ/l+kkuTPC7JPZOcnuTydrvxuOPUYCV5YJJv9rxuTvIqn/3Cl+TVSb6b5JIkxya5UzuR1rntc/90O6mWFpgkr2yf+3eTvKot82d+AUrysSQ/S3JJT1nfZ53GB5L8MMm3kzxqfJFrrqZ49s9uf+5vTbJsUv03ts/+siS7jz5iDcoUz/6f2n/jfzvJ55Ns1HNsXj17k8PRORL4UlU9CHgkcClwCHBGVW0LnNHuawGpqsuqaruq2g54NPBb4PP47Be0JFsCrwCWVdXDaCbP2gd4N/C+9rnfABwwvig1DEkeBrwE2IHmd/3TkmyLP/ML1Qpgj0llUz3rpwDbtq+DgA+PKEYNxwru+OwvAZ4JfK23MMlDaP4GPLT9zIeSLBlBjBqOFdzx2Z8OPKyqHgH8AHgjzM9nb3I4AknuDjwROAqgqv5QVTcCewFHt9WOBvYeT4Qakd2AH1XVT/DZLwbrAndOsi5wF+Aa4EnA8e1xn/vC9GDgnKr6bVWtBr4K/DX+zC9IVfU14JeTiqd61nsBx1TjHGCjJFuMJlINWr9nX1WXVtVlfarvBRxXVb+vqiuAH9L8B5LmoSme/X+3v/MBzqFZhx3m4bM3ORyN+wPXAx9PcnGSjybZENi8qq4BaLf3GmeQGrp9gGPb9z77Bayqrgb+GbiSJim8CbgQuLHnj8cqYMvxRKghugR4YpJNktwFeCpwX/yZX0ymetZbAlf11PN3wOLhs19c/hb4Yvt+3j17k8PRWBd4FPDhqtoe+A12KVpU2rFlzwA+O+5YNHztGKO9gG2A+wAb0nQpm8zpoheYqrqUpvvw6cCXgG8Bq6f9kBaL9Cnzd8Di4LNfJJK8ieZ3/icnivpUW6ufvcnhaKwCVlXVue3+8TTJ4nUTXUra7c/GFJ+G7ynARVV1Xbvvs1/YngxcUVXXV9UfgROAx9N0I5tYX3Yr4KfjClDDU1VHVdWjquqJNF2PLsef+cVkqme9iqYVeYK/AxYPn/0ikGRf4GnAC+q2tQLn3bM3ORyBqroWuCrJA9ui3YDvAScB+7Zl+wInjiE8jcbzuK1LKfjsF7orgccmuUuScNvP/FeAZ7V1fO4LVJJ7tdutaSanOBZ/5heTqZ71ScCL21lLHwvcNNH9VAveScA+STZIsg3NpETnjTkmDVCSPYA3AM+oqt/2HJp3zz63JbYapiTbAR8F1gd+DOxPk5x/Btia5h+Tz66qyQPbNc+1446uAu5fVTe1ZZvgs1/QkrwVeC5N95KLgQNpxhkcB9yzLXthVf1+bEFqKJKcDWwC/BF4TVWd4c/8wpTkWGAXYFPgOuBQ4Av0edbtfxT9G82Mhb8F9q+qC8YRt+Zuimf/S+Bfgc2AG4FvVtXubf030YxFWw28qqq+2Oe0mgemePZvBDYAftFWO6eqDm7rz6tnb3IoSZIkSbJbqSRJkiTJ5FCSJEmShMmhJEmSJAmTQ0mSJEkSJoeSJEmSJEwOJUmSJEmYHErSvJdkaZJKsmLcscxVkl3aezls3LHojpLs1z6fideXxh3TbCT59aT7WDrumCRpbWByKEkSkGSFicKMnQi8FfhEv4Np7J3k00muSPKbJL9LsirJqUlekWST2V48yZ+3z+rqJEvWUHfHtu63eorf2cb/rSk+JkmL0rrjDkCSNGdXAw8Gbhp3IANwHs29/HzcgWhaX6iqFf0OJLkP8GngCcBvga8Anwd+B9wbeDxwJPCOJH9WVdd1vXhV/SDJV4GdgT2Bk6apfmC7Xd7z+Xe2sS4FHtn1+pK0UJkcStI8V1V/BL4/7jgGoap+ywK5l8UoyYbAl4CHA58FXlpVv+hT77HAu4A7z+Fyy2mSwwOZIjlMcnfgOTRJat9WTknSbexWKknz3FRjDnu7SSb5P0m+03btuy7J8iT36Kl7pyQ3JvlZkr7/cZjk39vz7dlTtneSTyT5Qdt18NdJLmy7Dd7hb0xPTPdP8vIk305yS5Kz2uN9xxwmeXSSI5N8K8kv2/u4PMl7k2zc5zoTY+P2S7JrkrOS/CrJzUlOSfLgSfUL2LfdvaJnLNrK6b/6d7jWXyY5u/06XJ/k40k2auttn+TkJDe0x0/q14V1Fve6fvv1vqg992+TrExyYpInT6q7U5L/art3/j7JtUnOSXLomu5zhl5HkxieDTyvX2IIUFXnAE8CrupzPw9qv0+uamO8LsmnkjxwUtXPAb8Antq2VvbzfOAuwGeqaiG0rEvSUJkcStLC95729S3ggzTdUF9C09UPgKr6HU1XwM2Ap0w+QZINaFpgrgNO6zl0BPAo4FzgX4H/BO5K023w6GliOhI4HPhO+/7ra7iHlwD7AJcBHwf+HbgGeA3w9SR3m+JzTwP+G7i5/czZwFOBrybZtKde7/izI9v9twLvX0NcvZ4BnAJc317rcmA/4AttS9n/0PTYOYrmfp8OnNInie56ryvamNcDjgE+AHyNJknbY6JSkj2As2i6e54BvBf4AvB74O863Od0Dmi3h1fVn6arWI3b1WljvAh4AXA+zX2dATwTOC/Jo3o+/3ua77clwP5TXOYl7fYjHe9DkhanqvLly5cvX/P4BSwFClgxqXxFW34lsHVP+bo0yUMBO/SUP64tO77PNZ7dHnvvpPIH9Km7Dk1iWMBfTBHT1cA2fT67S3v8sEnl9wOW9Kl/QFv/DZPK92vLVwO7TTr2rvbY308R29KOX//ea+086etwenvsl8ALJn3uqPbYXrO9V+AewK3ABVN8ZpOe959rP//IPvU27Xiv+/U5tnV77I/ABrP4Pt4YuIFmvOlDJh17KPBr4KJJ5Q9ur/ljIJOObd8eu2Saa87qmfvy5cvXQn3ZcihJC9/bqurKiZ2qWk3TIgWwQ0/5/wN+ADw9yT0nnWOiy+XtWgOr6keTL1ZVt9K0+ADsPkVM76mqK2Z6A1X1k+rfEvUxmlbBqa5zXFWdMalsYmKSHSZXnqNjq+qrEzvt1+E/291LquqTk+of02636y3seK8FhKb179bJH6j+3Tpv6VNvEBMA3bvd/qKaVr3babsgHzbptUtPlRcDGwGHVtX3JsX3XZrWv+2TPKSn/FKaFtltgN0mXXKi1XA5kqQZcUIaSVr4LuhTNjHWa/IYtqOBd9B0a/wQQJLNaRKSi6vq272V0yxH8Hqarpr3BzacdL4tp4jpvJkG315nPeD/tHE9hKbFrPc/OKe6Tpd7n6t+1/ppu72wz7Gr2+1WvYVd7rWqbk7yXzRdVL+Z5HM0XWfPrWZyn16fpOmeeW6ST9PMIvr1qlo1g3ubiUyENcXxvbntPxl6ndVuH9duHzl5zGnrz9vtg4He5PEjNF1lXwJ8GSDJnWnGG/6O2xJ0SdIamBxK0sJ3Y5+y1e128hpxx9CMBdyXNjmkGf+1LpNaDduJVs6nabU5r/3sL9tzbwS8Ethgipiu7XQHzXjIv6bpPnhi+/mJ1qlXTXOdO9x7Va1OAne897nqN+HJ6hkcW29Sedd7fS7wBppk6K1t2e+SHA+8rtqlIqrqhCRPA14L/C1NAkqSC4E3VtXpa7rBNbim3W6aZIPJrYdVtR9Nt1SSHMgdxwFOrHv4EqZ310n7n6EZG7p3kk3bVtDn0CTVn6iqG7rchCQtZiaHkqT/VVWrkpwJPDnJg6rq+zSJ4h+BT02qfiBNYvjWqjqs90CSx9Ekh1NeaqYxJVlGkyx9GXhqNUt3TBxbB/j7mZ5rbTebe62qW4DDgMOS3Bd4Ik0S9kKa8ag79dQ9hWYSnA2Bv6CZsOelwMlJtp/cnbOLqroyyVXARAxdk82JBPqRk1uo13Dd3yX5BPBymq6p/0KftQ0lSWvmmENJ0mQr2u2+SbYDHgF8saqun1Tvz9rt5/qcY+cBxjNxnZN6k6XWDsxtrbxeE+P8Bt2i2MWc7rWqrmrHNu5OM1vqE9quv5Pr/aaqzqyq1wDvBNanzyy1s3BUu31Tv6VM1uCcdrvTtLX6m0gCD0zyIJpupt+vqrNncS5JWrRMDiVJk51AM/HJC2m7AXJbwthrZbvdpbcwyfbAGwcYz1TXuRfN0hyDMjF5y9YDPGdXK9vtLr2FU91rks2S/EWf82wI3I2m6+of2rq7tWPxJtu83U4eozgb7wUuofnPgU/2mdhowkZ9yj5O0w340CR3mCwoyTqTJrD5X1V1CU1y+WBuSxRdvkKSOrJbqSTpdqrqliSfpVk64e9okqZT+lQ9hmYymvcn2ZWmpWpbmq6KJ9CMhRuE82nWBXxmkm/QzE65OU1L12XcNunLXJ1Bcz8facfr/Rq4sar+bUDnn4mu97olcE6SS2nWB7wKuDvNM7g38IGq+lVb973A0iRn0SShfwAeTbMY/U+A4+YafFX9OsnuwGdpJtR5RttN+bL2epvTrIu5HfAr4OKez/4iybNo1t88J8kZwHdpZmHdmmbCmk2AO01x+eXAY2laHn/P9OtsSpL6sOVQktTPina7Hs0SDX+YXKGqfkrzD/FTaLrxvYxmjb6/Aw4ZVCDtsg7PAD4M3Ad4RXu9j9J0n5zc/XK21zmNZrKWPwKvppmY53WDOHeHGLre60rgUJpJa3YFXkMzI+kVNBPUvKqn7juBL9KsGXggcDBNsvZO4DGDmril5/vimcCpNN2S/28b2x7AdW1c21TViZM+e0Zb/0M04yUPbmN9GHAmTcI5lU9z27jFE6ZYxkOSNI1UzXhOAEmStIgl2Y+m++f+VbVivNHMXZIVNBMubVNVK8cbjSSNny2HkiSpq48nqSRfGncgs5Hk10mK/usuStKi5ZhDSZI0U9/ktrUUAX44rkDmaGKG1gn91gKVpEXHbqWSJEmSJLuVSpIkSZJMDiVJkiRJmBxKkiRJkjA5lCRJkiRhcihJkiRJwuRQkiRJkgT8fzKq9V6wdkzgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x21409118ba8>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Let's differ the invariant masses of the large and small pseudorapidity\n",
    "# events for making the histograms.\n",
    "inv_mass_great = great_etas['M']\n",
    "inv_mass_small = small_etas['M']\n",
    "\n",
    "# Let's use the matplotlib.pyplot module to create a custom size\n",
    "# figure where the two histograms will be plotted.\n",
    "f = plt.figure(1)\n",
    "f.set_figheight(15)\n",
    "f.set_figwidth(15)\n",
    "plt.subplot(211)\n",
    "plt.hist(inv_mass_great, bins=120, range=(60,120))\n",
    "plt.ylabel('great etas, number of events', fontsize=20)\n",
    "plt.subplot(212)\n",
    "plt.hist(inv_mass_small, bins=120, range=(60,120))\n",
    "plt.ylabel('small etas, number of events', fontsize=20)\n",
    "plt.xlabel('invariant mass [GeV]', fontsize=20)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Question 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the two histograms that were created above. In which way the pseudorapidities of the muons affect to the mass distribution?\n",
    "\n",
    "What could possibly explain your observations?\n",
    "\n",
    "First think the explanation by yourself and after that open the explanation below to see how you managed out."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Click here to open the explanation</summary>\n",
    "    \n",
    "    From the histograms one can see that the events where the pseudorapidity of both of the muons have been small, will produce more narrow peak than the events with the muons with great pseudorapidities. That means that the resolution of the invariant masses have been worse with the great pseudorapidities.\n",
    "    \n",
    "    The worse resolution follows from the fact that the resolution of the transverse momentum (p_t, the component of momentum that is perpendicular to the particle beams) has been worse for the muons with great pseudorapidities. This can be seen for example from the image 21 on page 32 of the paper https://arxiv.org/pdf/1206.4071.pdf of the CMS Collaboration.\n",
    "    \n",
    "    The explanation for the effect of the pseudorapidity to the resolution is that the particles which head to the endcap of the detector (great pseudorapidities) will more probably interact with the material of the detector than the muons with smaller pseudorapidities (check the image 11). In these interactions muons will lose some of their energy. This messes up slightly the fitting of the trajectories of the muons and the measurement of the transverse momentum. It also effects to the situation that the measurement of the transverse momentum is easily depended for example from the orientation of the muon chambers, from the division of the material of the detector and from the magnetic field. It can be assumed that these things are worse known for the particles that have great pseudorapidities.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## In the end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you have completed the exercise. Feel free to go back and test some different values to the code and see what happens. You can also create a new code cell by clicking \"INSERT\" -> \"Insert Cell Below\" and try to write some own code too!\n",
    "\n",
    "More information about the CERN Open Data can be found from http://opendata.cern.ch/."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sources"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[1] P. Mouche, *Overall view of the LHC. Vue d'ensemble du LHC*, 2014.\n",
    "Url: [https://cds.cern.ch/record/1708847](https://cds.cern.ch/record/1708847).\n",
    "\n",
    "[2] M. Brice, *View of an open LHC interconnection. Vue d'une interconnection ouverte*, 2005.\n",
    "Url: [https://cds.cern.ch/record/905940](https://cds.cern.ch/record/905940)\n",
    "\n",
    "[3] CMS Collaboration, *Detector Drawings*, 2012.\n",
    "Url: [https://cds.cern.ch/record/1433717](https://cds.cern.ch/record/1433717).\n",
    "\n",
    "[4] M. Lapka, D. Barney, E. Quigg et al., *Interactive slice of CMS detector*, 2010.\n",
    "Url: [https://cms-docdb.cern.ch/cgi-bin/PublicDocDB/ShowDocument?docid=4172](https://cms-docdb.cern.ch/cgi-bin/PublicDocDB/ShowDocument?docid=4172).\n",
    "\n",
    "[5] CMS Collaboration, *Performance of CMS muon reconstruction in pp collision events at $\\sqrt{s} =$ 7 TeV*, 2012.\n",
    "Url: [arXiv:1206.4071](https://arxiv.org/abs/1206.4071)."
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}