Exercise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Particle physics data-analysis with CMS open data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Welcome to the RAL Particle Physics masterclass computer exercise, here we will use real data from the CMS experiment at CERN for a simple particle physics data-analysis. \n",
    "\n",
    "The goal of the exercise is to understand how particles are discovered, as an example we will look at the <b> discovery of the Z boson </b>.\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 value, and the mass of the <b> Z </b> estimated.\n",
    "\n",
    "Finally, we will also look at 4-lepton events and try to identify the <b> Higgs boson </b>.\n",
    "\n",
    "The structure of the exercise is:\n",
    "- Theory background\n",
    "- Identifying events from event displays\n",
    "- Computer exercise:\n",
    "    - Introduction to computing and python \n",
    "    - Loading the data\n",
    "    - Making some plots\n",
    "    - Calculating the invariant mass\n",
    "    - Looking for Higgs to 4-lepton decays\n",
    "\n",
    "If you complete the exercise and have time left, there are two possible extension exercises:\n",
    " - The effect of pseudorapidity on the <b> Z </b> mass distribution\n",
    " - Fitting the <b> Z </b> mass distribution to determine the mass and lifetime of the <b> Z boson </b>\n",
    "  \n",
    "\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": [
    "## Part1 : Theory background"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle physics is the field of physics where structures of matter and radiation and the interactions between them are studied. In experimental particle physics,  research is performed by accelerating particles and colliding them either with other particles or with solid targets. This is done with _particle accelerators_ and  the collisions are examined with _particle detectors_.\n",
    "\n",
    "The world's largest particle accelerator, the Large Hadron Collider (LHC), is located at CERN, the European Organization for Nuclear Research. The LHC is a 27 kilometers long circle-shaped synchrotron accelerator. The LHC is located in a tunnel 100 meters underground on the border of France and Switzerland (image 1)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <center> <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> </center>\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 with a mass 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 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 the help of an electric field.\n",
    "\n",
    "The process of accelerating the protons starts before the LHC. Before the protons arrive in the LHC they are accelerated with electric fields and directed with magnetic fields in smaller accelerators(Linac 2, Proton Synchrotron Booster, Proton Synchrotron and Super Proton Synchrotron). After these the protons have an energy of 450 GeV. The protons are injected into the LHC in two different beampipes, each beam contains 2808 proton bunches located about 7.5 meters from each other. Each of these bunches include $1\\text{.}2\\cdot 10^{11}$ protons.\n",
    "\n",
    "The two beams circulate in opposite directions in two different vacuum tubes. Image 2 shows a part of the LHC accelerator opened with the two vacuum tubes visible inside. Each of the proton beams will reach the energy of about 7 TeV (7000 GeV) in the LHC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "   <center> <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> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Particle collisions are created by crossing these two beams that are heading in opposite directions. Because the bunches are travelling so fast, there will be about 40 million bunch crosses per second in the LHC. 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 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 of these collisions is 14 TeV. However in most cases the collision energy is smaller than that because when protons collide it is really their constitiuents, the quarks and gluons, which collide with each other. So not all of the energy of the protons is transmitted to the collision.\n",
    "\n",
    "When the protons collide the energy of the collision can be transformed into mass ($E=mc^2$) and new particles are produced in the collisions. These new particles are ejected from the collision area, a bit like a small explosion. By examining and measuring the particles created in collisions, researchers try to understand better the known particles which make up our universe and search for new particles which could explain puzzles such as dark matter. \n",
    "\n"
   ]
  },
  {
   "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 SHIFT + ENTER)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "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": [
    "At the LHC the proton beams are brought together to colide at four different points. In order to study the particles produced by the collisions, <b>particle detectors</b> are built around the collision points. \n",
    "The four particle detectors at the LHC are ATLAS, LHCb, ALICE and CMS (check Image 1). These detectors are like very large digital cameras and take a \"picture\" of the particles emerging from the collision.\n",
    "\n",
    "In Image 3 there is a visualisation of some particles created in one collision <b>event</b> seen at the CMS (Compact Muon Solenoid) detector."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "  <center>  <img src=\"images/eventdisplay.png\" alt=\"image missing\" style=\"height: 450px\" />\n",
    "     <figcaption>Image 3: A visualised collision event.</figcaption> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This exercise uses data recorded by the CMS detector so lets look in more detail at this detector....\n",
    "\n",
    "\n",
    "Simplified, the goal of the CMS detector is to detect particles that are created in collisions and measure different quantities about them (charge, energy, momentum, etc.). The CMS detector consists of different <b> sub-detectors</b> which form an onion-like structure around the collision point. This structure ensures that as many particles as possible from the collision are detected and measured. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <center> <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> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Different particles act differently in the different sub-detectors of CMS. Image 5 shows a cross-section of the CMS detector. The particle beams would travel in and out from the plane. Image 5 also demonstrates how different particles can be identified in the detector.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "   <center>  <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> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Let's look at the different parts of the detector: \n",
    "<dl>\n",
    "    <dt> <b>Tracker</b> </dt>\n",
    "    <dd> The 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 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. </dd>\n",
    "    <dt> <b>Calorimeter</b> </dt>\n",
    "    <dd> 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",
    " 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. The operation of the HCAL is also based on detecting light. </dd>\n",
    "    <dt> <b>Muon detector</b> </dt>\n",
    "<dd> Only muons and very weakly interacting particles like neutrinos will pass through both the ECAL and HCAL without being stopped. 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.</dd>\n",
    "</dl>\n",
    "\n",
    "Neutrinos can't be detected directly in the detector (they only interact very weakly and pass right through the detector), 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. Yet, we know that energy must be conserved. This situation indicates that something was undetected in the collision, this <b>\"missing energy\"</b> is assumed to be due to neutrinos created in the collision."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part2 : Looking at some events"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look at some more event displays by downloading the file <a href=\"Events/EventDisplays.pdf\"> here </a>\n",
    "\n",
    "If you want to look at more events they can be found at this <a href=\"http://opendata.cern.ch/visualise/events/cms#\"> link </a>\n",
    "\n",
    "Click the “folder” icon, click “Open files from the Web” and the “Education” folder"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Indirect detection of particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we have seen,  not every particle can be detected directly with the particle detectors. Interesting particles are often short-lived and <b> decay </b> essentially at the interaction point so never reach the detectors. These processes can be searched for via their long-lived decay products, this is indirect detection.\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 possible 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 in many ways (24 in fact) and in this exercise we will look at one of these: the decay of the Z to a muon ($\\mu^-$) and an antimuon ($\\mu^+$). This decay process is shown as a Feynman diagram in Image 6 below."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "   <center>   <img src=\"images/Zdecay.png\" alt=\"image missing\" style=\"height: 170px\" />\n",
    "   <figcaption>Image 6: Feynmann diagram of the process where the Z boson decays to a muon and an antimuon.</figcaption> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "As we have just seen in the event displays, the muons that are created from the decay of the Z can be detected. \n",
    "\n",
    "But just the detection of the muon and the antimuon isn't sufficient evidence for the existence of the Z as they could have originated from another process (there are many different processes which can lead to the same final state). \n",
    "Assuming that the muon, antimuon pair came from the decay of a single </b> \"mother\" </b> particle, we can use their momentum and energy to calculate the <b> invariant mass </b> of that particle.\n",
    "\n",
    "With the invariant mass it is possible to prove the existence of particles. \n",
    "\n",
    "In our example, we can take all the muon-antimuon pairs recorded by the detector and calculate the invariant mass for each pair. \n",
    "If we get a different answer each time then the muon-antimuon pair were just a random combination.\n",
    "If the answer is always the same it indicates that the muon-antimuon pair came from a single particle with a specific mass. \n",
    "\n",
    "We can make a plot showing the calculated mass value for each muon-antimuon pair. A peak in this plot (i.e. lots of pairs with the the same mass value) would prove that the the muon pairs came from a single particle with that specific mass vaule. __So the invariant mass can be used as an evidence about the existence of a particle__.\n",
    "\n",
    "In this notebook we will look at some real data from muon pairs, plot the mass of the muon pairs and look at the particles we find. Then we will find out how to calculate the mass ourselves.\n",
    "The different parts of the exercise are:\n",
    "\n",
    "1) Introduction to python, Jupyter notebooks and some simple programming\n",
    "\n",
    "2) Loading the data\n",
    "\n",
    "3) Making some simple plots\n",
    "\n",
    "4) Make a plot of the invariant mass of the muon pair \n",
    "\n",
    "5) Calculating the invariant mass yourself\n",
    "\n",
    "6) Apply the same principle to the 4-particle decay of the Higgs boson\n",
    "\n",
    "\n",
    "\n",
    "Now to get started......"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part3 : Computer exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1 : An introduction to python and programming"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a jupyter notebook, where you can have text \"cells\" (like this text here) and code \"cells\" i.e. boxes where you can write python code to be executed (like the one below). No need to install anything or find compilers, it is all done for you in the background.\n",
    "\n",
    "It is useful to save the workbook as you work through the exercise (just in case of problems),  use \"File\" -> \"Save Notebook\"\n",
    "\n",
    "We will be using python as the programming language: \n",
    "\n",
    "It is easy to get started, for example just type:  \"1 + 1\" in the cell below then click on \"Run\"->\"Run Selected Cell\"  above,  or click \"SHIFT\" & \"ENTER\" at the same time.\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "1+1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try some other maths functions for yourself: use \"-\", \"*\", \"/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now try something more advanced, for example sqrt(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sqrt(4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ooops, that failed: basic python can do some mathematical operations but not everything. For anything more complex, we need additional software packages or \"modules\".\n",
    "Here we import \"numpy\", a maths module: (run the cell below):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can try sqrt again using numpy: np.sqrt(4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Try out np.sqrt - This is a comment separated with #-symbol. \n",
    "np.sqrt(4)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can try some other values, yourself, e.g. np.sqrt(16), np.sqrt(81)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that starting a line with \"#\"  marks the line as a comment, this line doesn't affect the functionality of the code."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, you will need to be able to raise numbers to a power. This is done with \"** n\", where <b> n </b> is the power you wish to raise to. Try \"3**2\" in the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can try some other calculations as well. What is \"2\\*\\*4\", \"3\\*\\*3\" ??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 2 : Loading the data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data used in the analysis has been collected by the CMS detector in 2011. \n",
    "From the original data only those collision events with exactly two muons have been selected and the information stored on a CSV file. \n",
    "\n",
    "The CSV file used in this excercise is already saved to the same repository as this notebook file. Now let's get the file with Python and start the analysis!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the code cell below some needed Python modules _pandas_, _numpy_  are imported and named as _pd_, _np_. 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 Python's built-in functions.\n",
    "\n",
    "Run the cell below to import the data file ('DoubleMuRun2011A.csv'). Note that the file is 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"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "\n",
    "ds = pd.read_csv('DoubleMuRun2011A.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### How many events?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we want to figure out how many collision events (or in this case data rows) there are in the data file. We can do this using the `len()` function which tells us the length of an object. Inside the brackets will be written the variable whose length we want to be determine (in this case the file).\n",
    "\n",
    "Try to use the `len()` function to print 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.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "    The datafile was saved to the variable that was named as \"ds\".\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    "    \n",
    "    Use the function  len(). For example len(variablename), where 'variablename'\n",
    "    refers to the name of your variable (in this case 'ds').\n",
    "    \n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add your own code to print the number of collision events in the datafile!\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "    \n",
    "    len(ds)\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### What does the data 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`. \n",
    "\n",
    "Next we will print the five first rows of the file to look at what is inside. \n",
    "With the function _variablename_.`head(N)` you can get the first <b>N</b> elements of _variablename_. You can get the first rows of the data file by changing the _variablename_ to the name of your dataset variable.\n",
    "\n",
    "Write a code that prints the five first rows of the data file and run the code cell. First try to figure it out yourself, but if you get stuck click on the answer below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "    \n",
    "    ds.head(5)\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first row shows the information about muon pairs contained in the file. 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 muon (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 the momentum of the muon (remember momentum is a vector, $z$ is along the beamline, $x$ and $y$ are perpendicular to the beam)\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 the angle the particle makes with the beamline\n",
    "- phi = $\\phi$ = azimuth angle, also a coordinate describing an angle - this time in the x-y plane\n",
    "- Q = electrical charge of the muon"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 3  : Making some plots"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let's plot some of the values from the file in a histogram.\n",
    "\n",
    "A histogram describes how values are distributed, that is, how many values fall in each bin of the histogram. In Image 7 below 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 was 10–15 euros (12 people had this)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "  <center>   <img src=\"images/histogram.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "   <figcaption>Image 7: An example histogram from the distribution of the amount of cash.</figcaption> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Histograms can be created in python with the _matplotlib.pyplot_ module. \n",
    "\n",
    "Run the cell below to import this module as _plt_.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can plot something.... Let's try _'px1'_  (this is the x-component of the momentum vector for muon 1)\n",
    "\n",
    "The function _plt.hist()_ is used to create a histogram by giving different parameters inside the brackets. \n",
    "The full list of parameters can be seen at https://matplotlib.org/devdocs/api/_as_gen/matplotlib.pyplot.hist.html.\n",
    "\n",
    "For now, we will only use the first three:\n",
    "\n",
    "plt.hist('variable name', bins = #BINS, range=(#low end of range, # high end of range))\n",
    "\n",
    "<ul>\n",
    "    <li> 'variable name' : a variable from which values the histogram is created (here \"px1\") </li>\n",
    "    <li> 'bins' : number of bins for the histogram </li>\n",
    "    <li> 'range' : the lower and upper range of the bins  </li>\n",
    "</ul>\n",
    "\n",
    "The function _plt.show()_ is used to display the histogram\n",
    "\n",
    "Uncomment the \"plt.hist\" line in the cell below and fill in some values for bins and range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# fill variable px1 with px1 from the file\n",
    "px1 = ds['px1']\n",
    "\n",
    "# now use plt.hist to make a histogram\n",
    "# You need to add values for the number of bins and the histogram range\n",
    "plt.hist(px1, bins= #bins , range=(#low range , #high range))\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "  <summary>Answer</summary>\n",
    "   <code>\n",
    "px1 = ds['px1']\n",
    "plt.hist(px1, bins=100, range=(-20.,20.))\n",
    "plt.show()\n",
    "   </code>\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can vary the bins and range until you have something suitable.\n",
    "\n",
    "\n",
    "We can add axes labels and a title using \"plt.xlabel(' label')\", \"plt.ylabel(' label')\" and \"plt.title(' label')\".\n",
    "Try that in the cell below"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First add your plt.hist() line here\n",
    "\n",
    "\n",
    "# add labels and title\n",
    "plt.xlabel('x-component of momentum [GeV]')\n",
    "plt.ylabel('Number of events')\n",
    "plt.title('Histogram of px for muon 1. \\n')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "    <code>\n",
    "\n",
    "plt.hist(px1, bins=100, range=(-20.,20.))\n",
    "plt.xlabel('x-component of momentum [GeV]')\n",
    "plt.ylabel('Number of events')\n",
    "plt.title('Histogram of px for muon 1. \\n')\n",
    "\n",
    "plt.show()\n",
    "    </code>\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "You can also plot some of the other muon properties using the variables we printed above.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 4  : Plotting the invariant mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "Next, let's look at the invariant mass, this has already been calculated and stored in the file as \"M\".\n",
    "\n",
    "Write the code to make a plot of the invariant mass.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# First fill a variable \"invariant_mass_1\" with the invariant mass (\"M\") from the file\n",
    "\n",
    "# Use \"plt.hist\" to make a histogram of the invariant_mass_1 values. Remember to input the number of bins and the range.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "     <code>\n",
    "invariant_mass_1 = ds['M']\n",
    "\n",
    "\\# remember to input number of bins and range (0.5-150 works well)\n",
    "no_bins = 500\n",
    "         \n",
    "\\# use plt.hist to plot the invariant_mass_1 variable\n",
    "plt.hist(invariant_mass_1, no_bins, range=(0.5,120.), color=\"darkgrey\")\n",
    "\n",
    "\n",
    "plt.show()\n",
    "     </code>\n",
    "</details>\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Looking at the muon pair invariant mass spectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Below is the histogram published by the CMS experiment of the invariant mass of muon pairs. Does it look like yours??\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <center> <img src=\"images/CMShistogram.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "   <figcaption>Image 8: 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> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Not quite.... That's because the CMS plot uses log scales on the axes to make the plot clearer.\n",
    "\n",
    "We can change our plot to log axes using <b> plt.yscale('log')</b> and <b>plt.xscale('log')</b> \n",
    "\n",
    "Try that in the cell below\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#You need to add you plt.hist line here\n",
    "\n",
    "# set log scales on x and y axes\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "     <code>\n",
    "         \n",
    "plt.hist(invariant_mass_1, no_bins, range=(0.5,120.), color=\"darkgrey\")\n",
    "plt.yscale('log')\n",
    "plt.xscale('log')\n",
    "\n",
    "plt.show()\n",
    "     </code>\n",
    "</details>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now it should look more similar.\n",
    "\n",
    "The plot shows a smooth 'background' of random coincidences and on top of that some 'peaks'\n",
    "\n",
    "Each of these peaks is evidence for a particle decaying to muon pairs. \n",
    "The peaks corresponding to known particles and have been given labels in the CMS plot.    \n",
    "You can use the Particle Data Group <a href=\"https://pdg.lbl.gov/2020/tables/contents_tables.html\"> website </a> if you want to know more about these particles.\n",
    "If we saw a peak at a point where no known particle was expected, this would be evidence of a new particle discovery."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Challenge - find the mass of the J/&psi; particle\n",
    "\n",
    "Now try changing the range of the histogram to look at different parts of the mass spectrum - you can zoom in on the individual peaks (particles).\n",
    "\n",
    "For example, take the J/&psi; particle which is made up from an charm and anti-charm quark. As a challenge using the plot above, identify a range to zoom in on and determine the 'J/&psi;' particle mass to one decimal place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "invariant_mass_1 = ds['M']\n",
    "\n",
    "#Set the values num_bins, min_x, max_x below to zoom in on the J/psi mass distribution\n",
    "# Too many bins and the distribution is too sparse. Too few bins and tour mass readings will note be precise. \n",
    "# You also have to pick a good range too (min_x and max_x)\n",
    "\n",
    "num_bins = ####\n",
    "min_x = ####\n",
    "max_x = ####\n",
    "\n",
    "# You need to add a plt.hist function with your bins and range, Change the range to zoom on different regions of the plot \n",
    "plt.hist(invariant_mass_1, num_bins, range=(min_x,max_x), color=\"darkgrey\")\n",
    "\n",
    "# remember plt.show() to plot the histogram to the screen    \n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- ### Question 2 -->"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 5 :  Calculating the invariant mass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have seen that the invariant mass can be used to identify particles. \n",
    "Now let's calculate the invariant mass of the muon pairs for ourselves.\n",
    " \n",
    "#### Equation for invariant mass\n",
    "\n",
    "First, we derive loosely the equation for the invariant mass.\n",
    "\n",
    "Let's assume we have a particle with mass $M$ and energy $E$ which decays to two particles with masses $m_1$ and $m_2$, and energies $E_1$ and $E_2$. \n",
    "\n",
    "Energy $E$ and momentum $\\vec{p}$ are conserved in the decay process so\n",
    "\n",
    "$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 conservation 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",
    "$$\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) \n",
    "\n",
    "$$\n",
    "M = \\sqrt{(E)^2 - (\\vec{p})^2} = \\sqrt{(E_1+E_2)^2 - (\\vec{p_1} + \\vec{p_2})^2}, \\qquad (2)\n",
    "$$\n",
    "\n",
    "\n",
    "For those that like maths, a fuller derivation of this can be found <a href=\"images/Invariant_mass.pdf\"> here </a>\n",
    "\n",
    "\n",
    "#### How to do this in python\n",
    "\n",
    "In python, you only need to write a proper equation once - since the 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",
    "Remember that you can use 'sqrt'  from the _numpy_ module that we named as _np_. You can get a square root with the function `np.sqrt()`.  Naturally inside the brackets there will be anything that is inside the square root or brackets in the equation too.\n",
    "\n",
    "__In the cell below write code__ that will calculate the invariant mass value for muon pairs in each collision event in the data file. \n",
    "\n",
    "You need to use the muons energy and momentum and then use equation 2 to calculate the invariant mass of the parent particle:\n",
    "\n",
    "\n",
    "The energy of each particle can be calculated from:\n",
    "$$\n",
    "E_1^2 = \\vec{p_1}^2 + m_{1}^2\n",
    "$$\n",
    "\n",
    "Remember that momentum is a vector so:\n",
    "$$\n",
    "\\vec{p_1}^2 = (p_1^x)^2 +  (p_1^y)^2 +  (p_1^z)^2 \n",
    "$$\n",
    "\n",
    "where $p_1^x$ is the $x$-component of the momentum of particle 1. \n",
    "\n",
    "And\n",
    "\n",
    "$$\n",
    "(\\vec{p_1} + \\vec{p_2})^2 = (p_1^x + p_2^x)^2 +  (p_1^y + p_2^y)^2 +  (p_1^z + p_2^z)^2 \n",
    "$$\n",
    "Save the values calculated in a variable called `invariant_mass`.\n",
    "\n",
    "There are some comments in the cell below to help you with the different steps.\n",
    "There are also some hints - only use these if you are really stuck!\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    " When you write different quantities of the equation to your code, remember to refer to the variable (datafile) from where you want to get the quantities. \n",
    "For example if you want to use the quantity \"pt1\", write \"ds.pt1\" to the code.   \n",
    "  \n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 2</summary>\n",
    " Use the equations above for each step, for example to calculate the momentum squared of muon1 :\n",
    "   <code> \n",
    "  p1_squared = (ds.px1)**2 + (ds.py1)**2 + (ds.pz1)**2\n",
    "   </code>\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 3</summary>\n",
    " To calulate the energy of muon1 :    \n",
    "    <code> \n",
    "        e1 = np.sqrt(p1_squared + (muMass**2))\n",
    "    </code>\n",
    "</details>\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    " \n",
    "# You need the Mass of the Muon to calculate the energy\n",
    "muMass = 0.105658   \n",
    "\n",
    "# Momentum squared for the two individual muons\n",
    "\n",
    "\n",
    "\n",
    "# Energy of the two individual muons\n",
    "\n",
    "\n",
    "# Total Energy of the two muons \n",
    "\n",
    "\n",
    "\n",
    "# Momentum squared of the muon pair vector (p1+p2) - remember to add the vectors before squaring\n",
    "\n",
    "\n",
    "# Invariant mass of the muon pair, save this in  variable called \"invariant_mass\"\n",
    "# invariant_mass = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "   \n",
    "    # You will need the Mass of the Muon to calculate the energy\n",
    "    muMass = 0.105658   \n",
    "\n",
    "    # Momentum squared for the two individual muons\n",
    "    p1_squared = (ds.px1)**2 + (ds.py1)**2 + (ds.pz1)**2\n",
    "    p2_squared = (ds.px2)**2 + (ds.py2)**2 + (ds.pz2)**2\n",
    "\n",
    "    # Energy of the two individual muons\n",
    "    e1 = np.sqrt(p1_squared + (muMass*muMass))\n",
    "    e2 = np.sqrt(p2_squared + (muMass*muMass))\n",
    "\n",
    "\n",
    "    # Total Energy of the two muons \n",
    "    epair =  e1 + e2\n",
    "\n",
    "\n",
    "    # Momentum squared of the muon pair vector (p1+p2) - remember momentum is a vector, so add x,y,z components before squaring\n",
    "    ptpair_squared = (ds.px1 + ds.px2)**2 + (ds.py1 + ds.py2)**2 + (ds.pz1 + ds.pz2)**2\n",
    "\n",
    "    # Invariant mass of the muon pair\n",
    "    invariant_mass = np.sqrt(epair**2 - ptpair_squared)\n",
    "   \n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, if you run the cell below, the code will print the first five mass values that are calculated and will tell if the calculation is correct. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\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": [
    "#### Creating the histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "Next, let's create a histogram from the invariant mass values that you have calculated. \n",
    "\n",
    "Here we want to focus on the Z boson, so set the range wisely to get the values near the mass of the Z boson. \n",
    "\n",
    "Try different numbers of bins to make a clear histogram. You can try different values and see how they affect the histogram.\n",
    "\n",
    "Add axes labels and a title of the histogram. \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 histogram 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 we want to histogram (in our case the invariant masses). The zeroes will be replaced with the number 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>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 4</summary>\n",
    "    \n",
    "   A good mass range for the Z boson is 60-120 GeV\n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write down the code to create and plot the histogram (use plt.hist as we did earlier).\n",
    "\n",
    "# Use plt.show() to print the histogram to the screen\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "    # Write down the code to create and plot the histogram (use plt.hist as we did earlier).\n",
    "    \n",
    "    plt.hist(invariant_mass, 100, range=(60.,120.))\n",
    "\n",
    "    # Use plt.show() to print the histogram to the screen\n",
    "\n",
    "    plt.show()\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "#### Question : Can you describe the histogram. What information can you get from it?\n",
    "\n",
    "<details>\n",
    "    <summary>Answer</summary>\n",
    "\n",
    "    The position of the peak of the histogram on the x-axis tells you the mass of the Z boson\n",
    "    \n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 6 : Looking for Higgs to 4 lepton decays\n",
    "\n",
    "Now that we can reconstruct invariant masses we can look to find the mass of the Higgs via its decay to two Z bosons. As the Z boson is not stable and decays we can identify the Z boson by its decay to two leptons as above. Consequently the Higgs boson can end up decaying to 4 leptons. We can look at the final states electron-positron and electron-positron ($e^+ e^- e^+ e^-$), electron-positron and muon-antimuon ($e^+ e^- \\mu^+ \\mu^-$) as well as $\\mu^+ \\mu^- \\mu^+ \\mu^-$\n",
    "\n",
    "<figure>\n",
    "    <center> <img src=\"images/Feymann_Higgs.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    "   <figcaption>Image 9: Feymann diagrams for Higgs to 4-lepton decays</figcaption> </center>\n",
    "</figure>\n",
    "\n",
    "To calculate the invariant mass of the Higgs we need to know the mass of the particles in the final state. There are three different mass configurations here so three different calculations. We could look at the invariant mass distribution of each diagram above and then add them together to get the final distribution. But because these are high energy collisions, the masses of the electron and muon are very small compared to the momentum of the particles: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### $p >> m$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " Therefore the mass contribution to the energy is negligible, so: \n",
    "\n",
    "$ E^2 = \\vec{p}^2 + m^2 \\approx \\vec{p}^2$\n",
    " \n",
    " Now we can add all the $e^+ e^- e^+ e^-$, $e^+ e^- \\mu^+ \\mu^-$ and $\\mu^+ \\mu^- \\mu^+ \\mu^-$\n",
    "  data together as their only physical difference is their mass and we have set this to zero. The mass equation is now:\n",
    "  \n",
    "$\n",
    "M = \\sqrt{(E)^2 - (\\vec{p})^2} = \\sqrt{(E_1+E_2+E_3+E_4)^2 - (\\vec{p_1} + \\vec{p_2} + \\vec{p_3} + \\vec{p_4})^2}   \\qquad (3)  $\n",
    "\n",
    "and\n",
    "$ {E_1}^2 = \\vec{p_1}^2 $\n",
    "\n",
    "If you noticed the energy of each particle is already in the dataset. For example taking the concatenated dataset below, the energy of particle 1 is ds2.E1, the energy of particle 2 is ds2.E2 ... etc.\n",
    "  \n",
    "Have a try!  Calculate the invariant mass of all the datasets in one go!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### 4 lepton invariant mass"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we load the data for the different final sets of particles\n",
    "ds_2e2mu_2011 = pd.read_csv('2e2mu_2011.csv')\n",
    "ds_2e2mu_2012 = pd.read_csv('2e2mu_2012.csv')\n",
    "ds_4e_2011 = pd.read_csv('4e_2011.csv')\n",
    "ds_4e_2012 = pd.read_csv('4e_2012.csv')\n",
    "ds_4mu_2011 = pd.read_csv('4mu_2011.csv')\n",
    "ds_4mu_2012 = pd.read_csv('4mu_2012.csv')\n",
    "\n",
    "# Here we concatenate the 6 datasets into one called \"ds2\"\n",
    "ds2 = pd.concat([ds_2e2mu_2011, ds_2e2mu_2012, ds_4e_2011, ds_4e_2012, ds_4mu_2011, ds_4mu_2012], axis=0, ignore_index=True)\n",
    " \n",
    "#Total Energy of the four leptons\n",
    "E =\n",
    "\n",
    "#Total momentum in the x direction of the four leptons\n",
    "px =\n",
    "\n",
    "#Total momentum in the y direction of the four leptons\n",
    "py =\n",
    "\n",
    "#Total momentum in the z direction of the four leptons\n",
    "pz=\n",
    "\n",
    "# Now calculate the invariant mass using Equation (3) above and assign it to a variable called 'invariant_mass_2e2mu'\n",
    "invariant_mass_2e2mu = "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>   \n",
    "   <summary> Hint 1  </summary>\n",
    "The energy of particle 1 is 'ds2.E1'. Sum the energies for the 4 particles:\n",
    " \n",
    "    E = ds2.E1 + ds2.E2 + ds2.E3 + ds2.E4\n",
    "    \n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>   \n",
    "   <summary> Hint 2  </summary>\n",
    "The x component of the momentum for particle 1 is 'ds2.px1'. Sum the x components of mementum  for the 4 particles:\n",
    "\n",
    "     px = ds2.px1 + ds2.px2 + ds2.px3 + ds2.px4\n",
    " \n",
    "And repeat for the 'y' and 'z' components\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>   \n",
    "   <summary> Hint 3  </summary>\n",
    "Use 'np.sqrt' to calculate the invariant mass using equation (2)\n",
    "    \n",
    "    invariant_mass_2e2mu = np.sqrt(E**2 - px**2 - py**2 - pz**2)\n",
    "    \n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Full answer</summary>\n",
    "    \n",
    "    # Here we load the data for the different final sets of particles\n",
    "\n",
    "    ds_2e2mu_2011 = pd.read_csv('2e2mu_2011.csv')\n",
    "    ds_2e2mu_2012 = pd.read_csv('2e2mu_2012.csv')\n",
    "    ds_4e_2011 = pd.read_csv('4e_2011.csv')\n",
    "    ds_4e_2012 = pd.read_csv('4e_2012.csv')\n",
    "    ds_4mu_2011 = pd.read_csv('4mu_2011.csv')\n",
    "    ds_4mu_2012 = pd.read_csv('4mu_2012.csv')\n",
    "\n",
    "    # Here we concatenate the 6 datasets into one called \"ds2\"\n",
    "    ds2 = pd.concat([ds_2e2mu_2011, ds_2e2mu_2012, ds_4e_2011, ds_4e_2012, ds_4mu_2011, ds_4mu_2012],   axis=0, ignore_index=True)\n",
    "  \n",
    "    # Total Energy of the four leptons\n",
    "    E = ds2.E1 + ds2.E2 + ds2.E3 + ds2.E4\n",
    "\n",
    "    #Total momentum in the x direction of the four leptons\n",
    "    px = ds2.px1 + ds2.px2 + ds2.px3 + ds2.px4\n",
    "\n",
    "    #Total momentum in the y direction of the four leptons\n",
    "    py = ds2.py1 + ds2.py2 + ds2.py3 + ds2.py4\n",
    "\n",
    "    #Total momentum in the z direction of the four leptons\n",
    "    pz = ds2.pz1 + ds2.pz2 + ds2.pz3 + ds2.pz4\n",
    "\n",
    "    # Now calculate the invariant mass using Equation (2) above and assign it to a variable called 'invariant_mass_2e2mu'\n",
    "    invariant_mass_2e2mu = np.sqrt(E**2 - px**2 - py**2 - pz**2)\n",
    "\n",
    "  \n",
    "</details>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "<b> Next we want to plot your mass values </b>. \n",
    "\n",
    "You can run the cell below to do this. You should see peaks where the arrows are, corresponding to the Z and Higgs bosons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Here we create a histograma and fill it with the values in 'invariant_mass_2e2mu'\n",
    "plt.hist(invariant_mass_2e2mu, bins=60, range=(45,180))\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 four leptons. \\n')\n",
    "\n",
    "plt.arrow(70, 10, 18, -1.8,length_includes_head=True, width=0.2, fc='r', ec='r')\n",
    "plt.arrow(125, 9, 0, -3.5,length_includes_head=True, width=0.5, fc='r', ec='r')\n",
    "plt.text(118, 10.5, 'Higgs Boson', fontsize=12)\n",
    "plt.text(60, 10.5, 'Z Boson', fontsize=12)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Combining the 4 lepton data we start to see hints of the Higgs particle decaying to 4 leptons at a mass of around 126 $GeV/c^2$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Compare to the CMS analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compare our distribution to the CMS analysis. Bear in mind the image produced below uses more data and a more sophisticated analysis. For example looking at events with more than 4 leptons where the additional leptons can come from other particles in the event "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "    <img src=\"images/CMS-HIGGSTOZZ.png\" alt=\"image missing\" style=\"height: 350px\" />\n",
    " <center>   <figcaption>Image 10: Distribution of the reconstructed four-lepton invariant mass in the low-mass range. &copy; <a href=\"https://arxiv.org/abs/1706.09936\">CMS Collaboration</a> [6] </figcaption> </center>\n",
    "</figure>"
   ]
  },
  {
   "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 of your 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": [
    "# Further Work"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extension exercise 1 : Effect of pseudorapidity to the mass distribution\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have finished all the exercises above and would like to do more, look at the sections below on fitting the Z mass plot and the effect of pseudorapidity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this final section, we will study how the <b> pseudorapidities </b> of muons that are detected in the CMS detector alter the mass distribution.\n",
    "\n",
    "Pseudorapidity (denoted by $\\eta$) is a measure of the angle the detected particle makes with the particle beam (z-axis). \n",
    "The angle itself is called $\\theta$ (see diagram below). \n",
    "\n",
    "Pseudorapity is then determined with the equation:\n",
    "\n",
    "$$\n",
    "\\eta = -\\ln(\\tan(\\frac{\\theta}{2}))\n",
    "$$\n",
    "\n",
    "From the image one can see that, in practise, a large pseudorapidity means that the particle has continued almost among the beam-line after the collision.\n",
    " And vice versa: a small pseudorapidity means that the particle is more perpendicular to the beam-line\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<figure>\n",
    "  <center>  <img src=\"images/CMSangles.png\" alt=\"image missing\" style=\"height: 300px\" />\n",
    "     <figcaption>Image 11: Quantities $\\theta$, $\\eta$ and $\\phi$ in the CMS detector.</figcaption> </center>\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",
    "   <center>  <img src=\"images/pseudorapidities.png\" alt=\"image missing\" style=\"height: 300px\" />\n",
    "    <figcaption>Image 12: Two particles with different pseudorapidities in the CMS detector.</figcaption> </center>\n",
    "</figure>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this final section, two different histograms will be made: one using only muon pairs with small pseudorapidities and one using only those with large pseduorapidities. We can then study how the pseudorapidities of the muons that are detected in the CMS detector affect the mass distribution.  "
   ]
  },
  {
   "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 `large_etas`. To the first one we will save only collision events where pseudorapidities of  both the detected muons are small (for example under 0.38). And to the second one we save only those events there the pseudorapidities are both large (for example over 1.52). Absolute values will be used because $\\eta$ can have both positove and negative values.\n",
    "\n",
    "Complete the code cell below by determining the variables `small_etas` and `large_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: `ds[(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 are 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 absolute value function.\n",
    "</details>\n",
    "\n",
    "<details>\n",
    "    <summary>Hint 3</summary>\n",
    "    \n",
    "    The first variable with the conditions is \"large_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": null,
   "metadata": {},
   "outputs": [],
   "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",
    "# Define new variables \"large_etas\" and \"small_etas\" which contain only those events in \"ds\" which satisfy some condition\n",
    "large_etas = # \n",
    "small_etas = # \n",
    "\n",
    "# Let's print out some information about the selection\n",
    "print('Total number of events = %d' % len(ds))\n",
    "print('Number of events where the pseudorapidity of the both muons is large = %d' %len(large_etas))\n",
    "print('Number of events where the pseudorapidity of the both muons issmall = %d' %len(small_etas))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Creating the histograms\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create separate histograms from the events with small and with large values of pseudorapidities. \n",
    "You need to fill \"inv_mass_large\" and \"inv_mass_small\" with the invariant mass of events in your large and small eta datasets.\n",
    "\n",
    "The cell will get the invariant masses for both of the selections and will create the histograms out of them near to the peak of the Z boson."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<details>\n",
    "    <summary>Hint 1</summary>\n",
    "    \n",
    "  You can access the invariant mass values ('M') for the large eta slection with:  large_etas['M']\n",
    "</details>\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Let's get the invariant masses of the large and small pseudorapidity\n",
    "# events for making the histograms.\n",
    "inv_mass_large = #\n",
    "inv_mass_small = #\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_large, bins=120, range=(60,120))\n",
    "plt.ylabel('large 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 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the two histograms that were created above. In what way does the pseudorapidity of the muons affect the mass distribution?\n",
    "\n",
    "What could possibly explain your observations?\n",
    "\n",
    "First try to think of the explanation by yourself, then you can open the explanation below to see how you did."
   ]
  },
  {
   "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 is small, produces a narrower peak than the events where the muons have large  pseudorapidities. That means that the <b> resolution </b> of the invariant masses is worse with larger 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) is worse for muons with greater pseudorapidities. This can be seen for example from image 21 on page 32 of the CMS paper https://arxiv.org/pdf/1206.4071.pdf\n",
    "    \n",
    "    The explanation for the effect of the pseudorapidity on the resolution is that the particles which enter the endcap of the detector (larger pseudorapidities) will more probably interact with the material of the detector than the muons with smaller pseudorapidities (check the image 8). 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. The measurement of the transverse momentum also depends on, for example, the orientation of the muon chambers, the amount of material in the detector and the magnetic field. It can be assumed that these things are worse known for particles that have larger pseudorapidities.\n",
    "</details>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extension exercise 2 : Fitting a function to the Z mass histogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To get information about the 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": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = pd.read_csv('DoubleMuRun2011A.csv')\n",
    "\n",
    "invariant_mass = ds['M']\n",
    "\n",
    "# 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 = [#THE INITIAL GUESS FOR GAMMA, #THE INITIAL GUESS FOR M, -2, 200, 13000]\n",
    "\n",
    "\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": [
    "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 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "The width of the decay (called $\\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 and is equal to $6.58 \\times 10^{-16}$ $\\rm{eV s}$.\n",
    "\n",
    "\n",
    "\n",
    "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 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What is the physical meaning of the Z?\n",
    "\n"
   ]
  },
  {
   "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).\n",
    "\n",
    "[6] CMS Collaboration, *Measurements of properties of the Higgs boson decaying into the four-lepton final state in pp collisions at $\\sqrt{s} =$ 13 TeV*, 2017. Url: [arXiv:1706.09936](https://arxiv.org/abs/1706.09936)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}