DS_Tutorial_How_much_Higgs?.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table><tr><td><img src='hzz.png' width=\"570\"></td><td><img src='ATLAS-figaux_01.png' width=\"800\"></td></tr></table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  How much Higgs? - Intro\n",
    "- **Overview**\n",
    "    - Welcome to the tutorial that covers the first half of the Data Science in High Energy Physics module. The notebook is based on statistical analysis of open data and simulation from the ATLAS experiment at CERN.  Even if the physics context of this data is not entirely clear to you, it should not be a problem as understanding anything beyond the next paragraph is not neccessary.\n",
    "\n",
    "- **The physics in a nutshell**\n",
    "    - The Higgs boson was theorised in the 1960s as the fundamental particle that is responsible for giving mass to other fundamental particles. The problem is, the Higgs boson is extremely difficult to produce in lab conditions. Hence the Large Hadron Collider had to be built at CERN in order to provide enough high energy particle collisions to produce a significant Higgs bosons such that we could confirm the existence of this elusive boson. In 2012, the ATLAS and CMS experiments at the LHC announced that they had finally discovered the Higgs boson and the Nobel Prize in physics was subsequently awarded to Peter Higgs and Francois Englert for the boson's prediction. One of the easiest ways to find a Higgs boson in the ATLAS data is to look for collisions (often called 'events') containing 4 *leptons* which can be either electrons or muons. Muons are the slightly more massive cousin of the electron. This is because the Higgs boson can decay into two Z bosons, which in turn each decay into two leptons, leading to the four lepton signal. This process is depicted in the left hand diagram above where (reading left to right) two gluons collide to form a Higgs boson, which subsequently decays to two Z bosons and then four leptons. The right had image above shows a computerised illustration of **real** collision from the ATLAS data that contained four muons. In this notebook we will search through the ATLAS data for these events and compare what we find to predictions. \n",
    "\n",
    "- **Technical details**\n",
    "    - The notebook contains code that performs a simple analysis of the data and simulation. All code is Python 3 and utilised quite a few python package which must be installed before running the notebook. All required packages are listed in the 'requirements.txt' file.  The code is interspersed with exercises where you are either asked to discuss the analysis in the previous cell and/or extend the existing analysis. Specifically, we will look at collisions contains 4 leptons (electrons or muons) and use a histogram of the masses of the 4 lepton system ($m_{llll}$) to play with some statistical ideas. In the last few excercises you are required to perform some analyses with less guidance.\n",
    "    - This notebook has been tested with python 3.10. I recommend setting up a python virtual environment or simailr and NOT using your system python (this will quickly become impossible to manage). \n",
    "\n",
    "- **Report**\n",
    "    - You need to submit a 3-4 page report on this tutorial containg:\n",
    "        - an introduction explaining, in your own words, the overall concept of the tutorial and the exercises therein. \n",
    "        - your solution for each exercise (including any plots you make or need to explain your answer)\n",
    "            - code you developed to perform the last two exercises may be includes as an appendix to the 3-4 page report\n",
    "        - a brief conclusion."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we import some crucial python modules that will allow us to analyse the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import uproot3 as uproot\n",
    "import uproot3_methods.classes.TLorentzVector as LVepm\n",
    "import matplotlib.pyplot as plt\n",
    "import time\n",
    "import infofile\n",
    "import numpy as np\n",
    "import mplhep as hep\n",
    "from scipy import stats"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need the following helper function to help use to calibrate the simulation, we dont need to worry about what it actually does."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_xsec_weight(sample):\n",
    "    info = infofile.infos[sample] # open infofile\n",
    "    xsec_weight = (lumi*1000*info[\"xsec\"])/(info[\"sumw\"]*info[\"red_eff\"]) #*1000 to go from fb-1 to pb-1\n",
    "    return xsec_weight # return cross-section weight"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we are just defining the names of the files from where we will read the data and simulation. All data and MC samples should be downloaded from the Data Science folder on the Physics Honours Amathuba site.\n",
    "The last file is the real collision data from ATLAS, the others are simulations of the physics processes we expect to have aoccured in the collisions. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "files = [\n",
    "\"mc_361106.Zee.4lep.root\",\n",
    "\"mc_361107.Zmumu.4lep.root\",\n",
    "\"mc_410000.ttbar_lep.4lep.root\",\n",
    "\"mc_363490.llll.4lep.root\",\n",
    "\"mc_363492.llvv.4lep.root\",\n",
    "\"mc_363356.ZqqZll.4lep.root\",    \n",
    "\"mc_345060.ggH125_ZZ4lep.4lep.root\",\n",
    "\"mc_341964.WH125_ZZ4lep.4lep.root\",\n",
    "\"mc_344235.VBFH125_ZZ4lep.4lep.root\",\n",
    "\"mc_341947.ZH125_ZZ4lep.4lep.root\",\n",
    "\"data.4lep.root\"  \n",
    "]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to select out the most interesting collisions from the data and simulation so that we can as clear as possible signal of Higgs boson production. The exact details of the cuts are not very important. In summary, we select collision in which four leptons (electrons or muons) have been detected. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n",
      "File has been successfully opened!\n",
      "<class 'awkward0.array.jagged.JaggedArray'>\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 432x288 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#------------DEFINING SOME VARIABLES AND OBJECTS FOR THE ANALYSIS----------------------------------------------------\n",
    "\n",
    "lumi = 10#fb^-1\n",
    "nBins = 34\n",
    "\n",
    "minMass = 80\n",
    "maxMass = 250\n",
    "\n",
    "bins_ar = np.linspace(minMass, maxMass, num=(nBins+1))\n",
    "mc_hist_list = []\n",
    "sample_names = []\n",
    "\n",
    "f = plt.figure()\n",
    "\n",
    "H_125 = np.zeros([nBins])\n",
    "H_bkg = np.zeros([nBins])\n",
    "\n",
    "btagWP77 = 0.6459\n",
    "\n",
    "#--------------------------------------------------------------------------------------------------------------------\n",
    "\n",
    "for file in files: #looping over data and simulation files\n",
    "    sample_name = file.split(\".\")[1] \n",
    "    sample_names.append(sample_name)\n",
    "    tree = uproot.open(file)[\"mini\"]\n",
    "\n",
    "    mcWeight, SumWeights, XSection, trigM, trigE, scaleFactor_PILEUP, scaleFactor_ELE, scaleFactor_MUON,scaleFactor_LepTRIGGER, lep_type, lep_pt, lep_eta, lep_phi, lep_E, lep_charge, lep_etcone20, lep_ptcone30, jet_n, jet_pt, jet_eta, jet_phi, jet_E, jet_MV2c10 = tree.arrays([\"mcWeight\", \"SumWeights\", \"XSection\",\"trigM\", \"trigE\",\"scaleFactor_PILEUP\", \"scaleFactor_ELE\", \"scaleFactor_MUON\",\"scaleFactor_LepTRIGGER\",\"lep_type\",\"lep_pt\", \"lep_eta\",\"lep_phi\", \"lep_E\", \"lep_charge\", \"lep_etcone20\", \"lep_ptcone30\", \"jet_n\", \"jet_pt\", \"jet_eta\", \"jet_phi\",\"jet_E\", \"jet_MV2c10\"], outputtype=tuple)\n",
    "    print(\"File has been successfully opened!\")\n",
    "    print (type(lep_pt))\n",
    "    \n",
    "    leplv = LVepm.TLorentzVectorArray.from_ptetaphi(lep_pt, lep_eta, lep_phi, lep_E)\n",
    "\n",
    "    lep_reliso_pt = (lep_ptcone30 / lep_pt)\n",
    "    lep_reliso_et = (lep_etcone20 / lep_pt)\n",
    "    sum_lep_type = lep_type.sum()\n",
    "    \n",
    "    jetlv = LVepm.TLorentzVectorArray.from_ptetaphi(jet_pt, jet_eta, jet_phi, jet_E)\n",
    "    jetlv = jetlv[jet_MV2c10.argsort()]    \n",
    "    tags = jet_pt[jet_MV2c10 > btagWP77]\n",
    "\n",
    "    trig_cut = ( (trigM==1) | (trigE==1))\n",
    "    lep_kinematics_cut  = ( (lep_pt.max() > 20000) & (lep_pt.min() > 7000) & (lep_eta.min() >-2.5) & (lep_eta.max() < 2.5))\n",
    "    lep_type_cut  = ((sum_lep_type == 44) | (sum_lep_type == 48) | (sum_lep_type == 52))\n",
    "    lep_iso_cut =  ((lep_reliso_pt.max() < 0.3) & (lep_reliso_pt.max() < 0.3))\n",
    "    lept_count_cut = (leplv.counts ==4)\n",
    "    lept_charge_cut = (lep_charge.sum()==0)\n",
    "    \n",
    "    # filtering events according to cuts above that select out interesting events\n",
    "    event_cut = ( lep_kinematics_cut & lep_type_cut & lep_iso_cut  & lept_count_cut & lept_charge_cut)\n",
    "\n",
    "    first_lep_p4 =  leplv[event_cut,0]\n",
    "    second_lep_p4 = leplv[event_cut,1]\n",
    "    third_lep_p4 =  leplv[event_cut,2]\n",
    "    fourth_lep_p4 = leplv[event_cut,3]\n",
    "    mcWeight = mcWeight[event_cut] \n",
    "    \n",
    "    scaleFactor_PILEUP = scaleFactor_PILEUP[event_cut] \n",
    "    scaleFactor_ELE = scaleFactor_ELE[event_cut] \n",
    "    scaleFactor_MUON = scaleFactor_MUON[event_cut] \n",
    "    scaleFactor_LepTRIGGER = scaleFactor_LepTRIGGER[event_cut] \n",
    "    \n",
    "    #construct 4-lepton system by adding 4-vectors of the leptons vectorially\n",
    "    llll_p4 = first_lep_p4 + second_lep_p4 + third_lep_p4 + fourth_lep_p4  \n",
    "    \n",
    "    \n",
    "    # make histograms of the m_llll distribution for simulation and data\n",
    "    if(file.split(\"_\")[0] == \"mc\"):\n",
    "        finalWeight = get_xsec_weight(sample_name)*(mcWeight)*(scaleFactor_PILEUP)*(scaleFactor_ELE) *(scaleFactor_MUON)*(scaleFactor_LepTRIGGER)\n",
    "        H, b = np.histogram(llll_p4.mass/1000.0, weights=finalWeight, bins=bins_ar)                        \n",
    "        mc_hist_list.append(H)\n",
    "        if(\"H125\" in file):\n",
    "            H_125 = np.add(H, H_125)\n",
    "        else:\n",
    "            #print(\"Sample Name = \" + str(file) + \" exp. num. events = \" + str(np.sum(finalWeight)) )\n",
    "            H_bkg = np.add(H, H_bkg)\n",
    "    else:\n",
    "        finalWeight = np.ones(len(mcWeight)) \n",
    "        sample_name = \"Data (10 fb^-1)\"\n",
    "        H_data, b = np.histogram(llll_p4.mass/1000.0, weights=finalWeight, bins=bins_ar)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With our events selected, we can compare our data with our expectations from theory. Our best theory of particle physics is known as the Standard Model, so the predictions here represent the expectations from the Standard Model after a simulation of how the ATLAS detector detects particles has been applied. \n",
    "\n",
    "Specifically, we plot a histogram of the number of collisions (\"events\") as a function of the mass of the four-lepton system. We use this mass because it distinguishes between the the Higgs signal, which has a peak at the higgs mass (125 GeV) and is zero elesehere, and the backgrounds which have a peak around 90 GeV and are flatter elsewhere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAELCAYAAADeNe2OAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnDUlEQVR4nO3de3wV9bnv8c9DAFFAkatUykVR1HoKSAoiKhfRiArY1qIcL0hRbI/d9VLrhuou4Lan2N3avS1279qiYMWox8uW4iuAF6JIMQoCKqJC3Ygg91YxFBGS5/yxJulKyEpWwpp1m+/79cqLtSazZj0z8+PJzDO/+Y25OyIiEh3NMh2AiIiklxK/iEjEKPGLiESMEr+ISMQo8YuIRIwSv4hIxISa+M2snZk9aWbvmdk6MxtsZu3N7HkzWx/8e2yYMYiISE1hH/H/B7DQ3U8B+gLrgCnAi+5+EvBi8F5ERNLEwrqBy8yOAVYDJ3jcl5jZ+8Awd99qZl2BUnfvE0oQIiJyiOYhLrsXsBN4yMz6AiuBm4Au7r41mGcb0KWuD5vZZGAyQOvWrQeccsopIYYqIpJ/Vq5cucvdO9WeHuYRfyHwGjDE3cvM7D+APcA/uXu7uPn+5u711vkLCwt9xYoVocQpIpKvzGyluxfWnh5mjX8zsNndy4L3TwJnANuDEg/BvztCjEFERGoJLfG7+zbgYzOrqt+fB7wLzAcmBNMmAM+GFYOIiBwqzBo/wD8B88ysJfAhMJHYH5snzGwS8BEwLuQYREQkTqiJ391XA4fUl4gd/UfGgQMH2Lx5M1988UWmQxGRPNSqVSu6detGixYtkpo/7CN+ATZv3kzbtm3p2bMnZpbpcEQkj7g7u3fvZvPmzfTq1Supz2jIhjT44osv6NChg5K+iKScmdGhQ4dGVRSU+NNESV/y0fTp0zMdgtD4/KJSTwb0nPJcSpazcebF9f9+40b69+9P37592b9/P4WFhfzmN79JevnTp0+nd+/eXHXVVYcb6mHr3bs3GzZsCO8Lph+TouV8Vu+vN27cyHXXXccLL7xQPa13797MmjWLnTt3cvXVV6cmjpBVVFRQUlLCjBkzKCwsZNSoURQUFDR5eZMnT+aDDz6gtLSUqVOnsnz5crZt28b+/fvp0aMHvXv35uyzz2bz5s3ceeed9X6+yrRp01i0aBEtW7akoKCAJUuWHPK5hx9+mN/97nc0b96czp07M2vWLLp0qfOe0sMyZ84cZsyYQY8ePdi3bx/9+/fnnnvu4ZhjEre7++67jx/+8IcpjwV0xJ/3BgwYQGlpKcuXL+fdd99l7dq1oX9nRUVF6N+Rby688MKcSvpFRUWMHz8egPHjx1NUVNTk/f7ll1+yZs0a2rZty6ZNm/j5z39OaWkpU6ZMYdKkSZSWlvKHP/wh6c8DrFu3jmXLlvHaa6/xyiuv8NRTTx3yuaVLl/Lggw/ywgsv8PLLLzNx4sRQD3Kq1qWsrIw+ffpw00031Tv/fffdF1osSvwRcfDgQfbt20fbtm0pKipi2LBhDBw4kOXLlwOwadMmLrzwQoYOHcrIkSOprKys/uyePXu49NJLWbhwIZ9//jkXXXQRI0eO5NZbb2XYsGFA7Ozg2muvZcyYMTzxxBPMmDGDwYMHM2jQIJ57LnaGc+211/Lqq68C8Mgjj1SXCYYNG8bNN9/MBRdcwHnnncf+/fsB+PGPf8zgwYOZOHEiX375JQClpaUMHDiQ4cOHM3HixHRsurSYM2cOd999NwCPP/44ffv25dvf/jZFRUXVR7E/+tGPGDx4MN/73vfo0aMHkJntUVJSQllZGeXl5QCUl5dTVlZGSUlJk5b33HPPMWbMGCZMmMCjjz6aks+3atWKbdu2sXLlSiorK2nfvv0hn3v44Yf5yU9+wpFHHgnARRddRHl5OVu2bGHOnDmMHj2a0aNH079/f5YuXQrA22+/zciRIxkxYgTjxo1j3759AHTv3p0bbriBM888k9tuu63BmG+55RaWLl1KZWUljz76KMOHD2fw4MFcd911uDv33nsvW7ZsYdiwYcyePZslS5YwfPhwzjnnHMaOHXvYPQSV+PPcypUrGTZsGKeddhrdunWje/fuPP3005SWljJ37lzuuOMOAG677TZuueUWXn75ZRYvXkyzZrGmsXXrVi677DLuuOMOLrzwQn7/+99z7rnn8sILLzBgwIAa33XEEUcwf/58Tj31VJYuXcqf//xnFi1axC233FLjD0ldhg0bxuLFiznxxBN5/vnnWbVqFW+//TbLly9n2rRpbN0aG97p6aef5u6772bJkiXMnj07hC0Wvqp9UvUTr6Kign/5l3/h1Vdf5bHHHmPz5s0AvPnmm6xdu5bly5czdepUPvnkEyAz22PVqlXs3bu3xrS9e/eyevXqJi2vuLiYq6++mtGjRzfpj0ddn+/Vqxc/+9nPuPPOO+nZsydTphw6CPDmzZvp3r17jWk9evRgy5YtQKwb9p/+9CeeeeYZbrnlFgBuvPFGHnzwQV566SWGDBlSvc137NjBjBkzWL58OQsWLGDPnj0Nxt2pUyd27drF2LFjWbJkCcuXL+fzzz9n6dKl3HrrrRx//PGUlpYyadIkBg4cyJIlS1i6dCmnnHIKTzzxRKO3Uzwl/jxXVer54IMP6Nq1K8XFxdx0002cc845fP/73+fjjz8GYO3atYwYMQKgOulD7HRzxIgRfOMb3wBg/fr11a8HDRpU47vOOussAN5//33OPPNMzIx27drRuXNndu3aVeMCVO0xoqr+iHTv3p3du3fzwQcfVH9Pz549q+uuP/7xj5k/fz5XXnklDz30UGo2UppV7ZOqn3i7du2iS5cutG3blhYtWtC/f3+g5nbv0aNHRrdH//79ad26dY1prVu3pl+/fo1e1meffcayZcuYPHkyY8eOZePGjaxZsyYlnx87diwlJSWsX7+e1atXH1LjP/7446tLQ1U2bdrE8ccfD1Cj/X32Wezazdq1a7nmmmsYNmwYxcXFbNu2rXpZxx13HGZGt27d+Nvf/tZg7Dt37qRjx4688sorjBgxgqFDh1JWVlb9fzLe2rVrueCCCxg6dCjPPvtsnfM0hhJ/hBx77LFs2LCBgoICli5dym9/+9vqBPy1r32tOgnFH53/67/+K2vWrGHOnDlA7EJk1YB5b7zxRo3lV13cO/nkk3nttddwdz799FN27NhBx44dad++ffUR7MqVK2t8tvYfhZNOOql6nk2bNrF9+3YAOnTowKxZs3jkkUeYOXNmUkdWuaRjx45s376d8vJyDh48WH0U3bt376zZHqNGjWLQoEG0adMGgDZt2jBo0CBGjRrV6GU9+eSTTJ06lYULF7Jw4UIefPBB5s2bd9if/+tf/8qOHbFhwI444gjat29/yFnn1VdfzcyZM6vLNQsXLuTII4+sTvzx2/voo48G4PTTT6e4uJjS0lJee+01fvrTnwKH9qppaPDL++67jyFDhtCsWTOmTJnCvHnzePnllxk0aFD1Z+MPwH72s58xY8YMXn75ZcaMGdPg8huiXj15rqqs4O4cffTR3H///Vx22WWMHDmSIUOGVM/3y1/+kuuvv567776bFi1asHjxYgCaN2/OI488wsSJEzlw4ADXX38948aNY/HixZxyyim0bNnykO/s378/Z511FoMHD6ayspJf/epXNGvWjOuuu47x48fz6KOP0rFjR9q1a5cw7jPOOINTTz2VwYMHc/rpp/OVr3wFgHvvvZfFixdTWVnJ+eefX/0fMl8UFBQwffp0zj77bHr16kXnzp1p2bIlAwYM4OSTT67eHlXJKRPbo6CggEWLFlFSUsLo0aMpLi5ucq+eefPm8cADD1S/P/vss7nxxhuZOXNmnfPPmTOn+gDl3HPPpbS0tM7P33DDDXz3u9+loqICd2fAgAHVZ7RVhg4dytVXX83IkSNp3rw5HTt25JFHHqn+/VFHHcXFF1/MJ598wq9//WsA7r//fq699loOHDgAwNSpUzn//POTWtfZs2fzwgsvsG/fPr7+9a9XX7y95pprOP/886k99PzgwYP55je/yeWXX84VV1zBpEmT6NOnD8ccc8xh7+fQhmVOpVwflnndunWceuqpmQ4jZQ4ePEjz5s2ZN28ey5cvZ9asWZkOKa8cOHCAFi1acODAAQYMGMDixYs57rjjqqd/9NFHjB07tsk19VQys8M++sxGc+bMSdh1NFvVlWcSDcusI35plMrKSoYPH46ZYWb88Y9/zHRIeWfOnDnMmzePPXv2cM0113DccccBcPPNN/POO+9QXl7OL3/5ywxHGTNt2rRMhyBNoCP+NMi3I34RyT6NOeLXxd00yYU/sCKSmxqbX5T406BVq1bs3r1byV9EUq5qdM5WrVol/RnV+NOgW7dubN68mZ07d2Y6FBHJQ1Xj8SdLiT8NWrRokfQ42SIiYVOpR0QkYpT4RUQiRolfRCRilPhFRCJGiV9EJGKU+EVEIkaJX0QkYpT4RUQiRolfRCRilPhFRCIm1CEbzGwj8DlQARx090Izaw88DvQENgLj3L3hB1SKiEhKpOOIf7i794sbE3oK8KK7nwS8GLwXEZE0yUSpZywwN3g9F7g0AzGIiERW2InfgcVmttLMJgfTurj71uD1NqBLyDGIiEicsIdlPtvdt5hZZ+B5M3sv/pfu7mZW59NJgj8UkwG6d+8ecpgiItER6hG/u28J/t0BPAMMBLabWVeA4N8dCT77gLsXunthp06dwgxTRCRSQkv8ZtbazNpWvQYuAN4B5gMTgtkmAM+GFYOIiBwqzFJPF+AZM6v6nkfdfaGZvQE8YWaTgI+AcSHGICIitYSW+N39Q6BvHdN3A+eF9b0iIlI/3bkrIhIxSvwiIhGjxC8iEjFK/CIiEaPELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiEaPELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiEaPELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiEaPELyISMUr8IiIRo8QvIhIxSvwiIhGjxC8iEjFK/CIiERN64jezAjNbZWYLgve9zKzMzDaY2eNm1jLsGERE5B/SccR/E7Au7v09wK/dvTfwN2BSGmIQEZFAqInfzLoBFwN/CN4bMAJ4MphlLnBpmDGIiEhNYR/x/ztwO1AZvO8AfOruB4P3m4Hj6/qgmU02sxVmtmLnzp0hhykiEh2hJX4zuwTY4e4rm/J5d3/A3QvdvbBTp04pjk5EJLqah7jsIcAYM7sIaAUcDfwH0M7MmgdH/d2ALSHGICIitYR2xO/uU929m7v3BK4AXnL3K4ElwGXBbBOAZ8OKQUREDpWJfvz/DNxqZhuI1fxnZyAGEZHICrPUU83dS4HS4PWHwMB0fK+IiBxKd+6KiESMEr+ISMQo8YuIRIwSv4hIxCjxi4hEjBK/iEjEKPGLiESMEr+ISMQo8YuIREyDid/MhphZ6+D1VWZ2r5n1CD80EREJQzJH/P8J/N3M+gI/Av4CPBxqVCIiEppkEv9Bd3dgLDDL3e8H2oYbloiIhCWZQdo+N7OpwFXAuWbWDGgRblgiIhKWZI74Lwf2A5PcfRuxh6f8W6hRiYhIaJI54r/F3f+56o27bzKzr4UYk4iIhCiZI/7z65g2KtWBiIhIeiQ84jez7wP/BzjBzN6K+1Vb4M9hByYiIuGor9TzKFAC/ByYEjf9c3f/a6hRiYhIaBImfnf/DPgMGG9mBUCXYP42ZtbG3TelKUYREUmhBi/umtkPgOnAdqAymOzA18MLS0REwpJMr56bgT7uvjvkWEREJA2S6dXzMbGSj4iI5IFkEv+HQKmZTTWzW6t+wg4s202fPj3TIYiINEkyiX8T8DzQklhXzqqfSJsxY0amQxARaZIGa/zuPgPAzI5y97+HH5KIiIQpmfH4B5vZu8B7wfu+Zvbb0CMTEZFQJFPq+XegCNgN4O5rgHNDjElEREKU1KMX3f3jWpMqGvqMmbUys9fNbI2ZrTWzqpJRLzMrM7MNZva4mbVsQtwiItJESXXnNLOzADezFmZ2G7Auic/tB0a4e1+gH3ChmZ0J3AP82t17A38DJjUtdBERaYpkEv/3gBuB44EtxJL4jQ19yGPKg7ctgh8HRgBPBtPnApc2KmIRETksydy5a+5+ZVMWHozxsxLoDdxP7Hm9n7r7wWCWzcT+oNT12cnAZIDu3bs35etFRKQOyRzxLzOzxWY2yczaNWbh7l7h7v2IPbVrIHBKIz77gLsXunthp06dGvO1IiJSjwYTv7ufDNwJfA1408wWmNlVjfkSd/8UWAIMBtqZWdWZRjdi5SMREUmTZHv1vO7utxI7av8rsdp8vcysU9UZgpkdSexJXuuI/QG4LJhtAvBs48MWEZGmSmZY5qOBbwJXACcCzxD7A9CQrsDcoM7fDHjC3RcEN4M9ZmZ3A6uA2U0NXkREGi+Zi7trgP8G7nL35cku2N3fAvrXMf1DkvvDISIiIUgm8Z/g7m5mR4UejYiIhC6ZGv+ZGqtHRCR/aKyeOmisfRHJZ6GN1ZPLNNa+iOSzZGr8NcbqAW4iubF6REQkC4U2Vo+IiGSnZO7c3eXuV7p7F3fv7O5XufvudASXjSoqKliwYAEACxYsoKIiv6peur4hkv/M3TMdQ4MKCwt9xYoVafs+M6Ou7VJRUUFRURFlZWWUl5fTpk0bBg0axKJFiygoKEhbfGFKtO4iknvMbKW7F9aentTFXYkpKSmpTvoA5eXllJWVUVJSkuHIRESSp8TfCKtWrWLv3r01pu3du5fVq1dnJiARkSZI5mHrd8a9PiLccLJb//79ad26dY1prVu3pl+/fpkJSESkCRImfjP7ZzMbzD9G0gRIeqyefDRq1CgGDRpEmzZtAKpr/KNGjcpwZCIiyavviP894DvACWa21Mx+D3Qwsz7pCS37FBQUsGjRIoqLiwEoLi7Oqwu7IhINCXv1mNlQoAz4M/AN4FTgOeAloI+7n5WuILOlV09j58lF+bpeIlGUqFdPfXfuFgE/JTYG/73AW8Bed58YTogiIpIOCUs97v4Tdz8P2Aj8ESgAOpnZq2b2pzTFJyIiKZbMWD2L3H0FsMLMvu/uZ5tZx7ADExGRcCQzZMPtcW+vDabtCisgEREJV6Nu4ArG4hcRkRymO3dFRCJGiT9Ovo+8KSICSvzVqkbeHD9+PADjx4+nqKhIyV9E8o4Sf0Ajb4pIVCjxBzTypohEhRJ/IOojb+r6hkh0KPEHojzypq5viESLEn8gyiNv6vqGSLSElvjN7KtmtsTM3jWztWZ2UzC9vZk9b2brg3+PDSuGxiooKOCSSy4B4JJLLqk36U+bNi1dYYVO1zdEoiXMI/6DwI/c/TTgTOBGMzsNmAK86O4nAS8G73PO9OnTMx1CykT9+oZI1ISW+N19q7u/Gbz+HFgHHA+MBeYGs80FLg0rBklOlK9viERRWmr8ZtYT6E/swS5d3H1r8KttQJcEn5lsZivMbMXOnTvTEWZkRfn6hkgUhZ74zawN8BRws7vvif+dxx71VOfjntz9AXcvdPfCTp06hR1m5DXm+oZkj2RKjvlUlpTUSPjoxZQs3KwFsIDYmP73BtPeB4a5+1Yz6wqUunu9z/HNxkcv5qsor3suivJjQqVhiR69GGavHgNmA+uqkn5gPjAheD0BeDasGERE5FDJPIGrqYYAVwNvm9nqYNpPgJnAE2Y2CfgIGBdiDCIiUktoid/dXwUswa/PC+t7w9JzynMNzrNx5sVpiERE5PDozl0RkYhR4hcRiRglfhGRiFHiF5HIiuo9Dkr8dcinAdhEJLEZM2ZkOoSMUOKvQ1SPAkQkGpT4RUQiRolfalCZS+qis+D8osQvNeg/uNQlqrXwfKXELyISMUr8IjmooqKCBQsWALBgwQIqKiqaNI9EkxK/SI6pqKigqKiI8ePHAzB+/HiKiopqJPZk5pHoUuIXyTElJSWUlZVRXl4OQHl5OWVlZZSUlDRqHomuMIdlzhnJjLwpki1WrVrF3r17a0zbu3cvq1evrn6KWjLzSHTpiF8kx/Tv35/WrVvXmNa6dWv69evXqHkkupT4pdHU5TOzRo0axaBBg2jTpg0Abdq0YdCgQYwaNapR82SC2k52UOKXRlOf7swqKChg0aJFFBcXA1BcXMyiRYsoKCho1DyZoLaTHVTjj7rpxyQxz2fhxyGNUlBQUF2rT1SzT2aeqKqoqKi+0L1gwQJGjRqV8T+K6aQjfhGJFHV1VeIXkRRJVf0+7OsAmerqmk3XN5T4RSShxtz9m6r6fdjXAerr6hqmbLq+ocQvInXK15KIuroq8YtIAvl692+2dnVNJyV+EalTpkoiYcvWrq7ppMQfomy6mCOpEaV9ms0lkcPdD7W7utaX9JP5rlxrF+bumY6hQYWFhb5ixYrQlp+qsXo2zry4xnszI+u3bxP68efEeoUk29Y9mXiaGnNVjb+q3FNVEkl0dJyqWMJcp0zGk4m2Y2Yr3b2w9nQd8YvksGQeldnUx2mqJJK/Qkv8Zvagme0ws3fiprU3s+fNbH3w77Fhfb9IFIRdhmhMSURyR5hH/HOAC2tNmwK86O4nAS8G7yUs049p+KcRGtOnO1du5smUfF2vRPQ0sOSkq12EWuM3s57AAnc/PXj/PjDM3beaWVeg1N37NLQc1fibqJGJPfFyPgul3puMjG/DWvJ1vRpyODXsTF0rSEa21fhT3S6ypcbfxd23Bq+3AV3S/P3SRPnap1vCp7aTfTJ2cddjf9YS/mkzs8lmtsLMVuzcuTONkUld8rVPtxy+hso4ybadfC0HZeN6pTvxbw9KPAT/7kg0o7s/4O6F7l7YqVOntAWYq+qqDfb84tEGf5IVRp/udF4HiFpNPR16TnmOHrfPp/WJAxjzre8AMOZb36H1iQPocfv86hJqMm0nX4eHyNb1Snfinw9MCF5PAJ5N8/fnrbAHgArjNvd0DuqVTQNk5ZN9H67ky0/exw98AYAf+IIvP3mffR+urJ4nmbaTqXJQU7u6Jitby1xhducsBpYDfcxss5lNAmYC55vZemBk8F5ygPp0S12+3P4X/MD+GtP8wH6+3PFh9ftk2k6mSolhnwlma4k0tCdwufv4BL86L6zvzBa59HSfZMo9G4PeQQVA1XOcLllxJcR3tMrzp3Tl0j5Np5ZdTsRaHFF9xA9gLY6gZecTgNo95gyAH7xq8OrC6qkbZ15cXQ6qOjKGukuJ6dwPyXxXQ/Mku17ppkcvplh81zWI1fTq67om2S+K+7R2F+djhoyvs9vzkScMoOVX+lSXe6xFK1p+pQ9HnjCgUd/llZUc6HAitv8fyznQ4URufKWSH7z6HBtnXpzW/ZDMdyUzT1WZq3ZX1kyPBKrEn2L11fQy/dzTT1+dR7uzr8xoDNls+vTpdZ76Z/M+ra0p96Qk0y4S/d6aFdBl3F3s+3AlO5+6i45jbufIEwZgzRqXiJNZTkP7ofa617Vete+1SSSZfZ7MPFVlrpKSEkaPHk1xcXFWnC1qrJ4Uy9aaHsBny4ozHUJWS3QBOJv3aSocbruwZgUc1XsgAEf1HtjopJ/schq7Hw5nvZL5rmTjycZhL5T4Uyybh7JNh1T1Wc6mvs+p3KfZtF65pOeU5/jNmoPQ/Iiav2h+BPetPpCyu++rJLPPc7ldqNSTYtla00uHVNVgs62mnqp9erjrlerklmtScT0hWcns82xpF00RySP+MLtwhdHtMVduPkpVn+Vs6/vc2H2aaH81dr1yZb+nS9V1gI5jbgeg45jb6TLuriaXluqTzD5P1f/1TLT3SCb+sG/mSXVNr754M1E6mDa0ZZ3TU3Vrfipv8U/V9mnMPk3VtYJ8vOnsmCGJenkn53CvJ/Sc8lyDP1Xq2+dV8554x8JY91Ri3VRPvGPhIctpSCauIUUy8eeLTN0OPn1Yqzqnp+rW/FQtJ9tul4/69R9I3DsoyjLRLpT4c1i2lUQaujW/55Tn6Hr5Xbz0yrIaMb/0yjK6Xn5X0suB5NY917aPZIe6zgASnRWkQibahS7uZqG6GlZd0z5d9hjl5TVPEcvL93LVLx6nXXD6mZJ4krm7l+T6LDfmFv/6llPf6XHV6Xky86RTtvbpTpXDLeOkejnZJtG4QJloF0r8Oayh2+WreGVF9aBZf9/wepNusElW7bpoU2NuaDnJ3AqfjcMANLRe6ZTqdpGqMs7hLied7b0xal+sP/Rg7tAhLZK94ayxlPhzWDLd27yygu1P/JQvP3kfgF3zf0HLr/QJrTdEKmJORqq622Vb19F0ybZ2kSrZtF7Z3P1WiT9EYQ/5msxt7vUNm1vVOyKd6os50X+URNP9jB9yZLuVlD91F0cW3cr6EwZw4h0Lq4+S6juFrlrm3ze8zq5XllVvn/hrDjuePLRXTdj7NF2yrV00VqJyUBjrlUzpKdfKU5FK/OkeYfFw+2Enc8pau3tbbfXV1DP1H7yhmFO5nIZKK8lcc4hX3z49nPYV/8ctfr93vmxaKKWKbGwXjZGoHBTGeiVTesq13kqRSfz1ndKnSjKndsnW7FJ1yppsTf1wJXtBOh0ac+ZwuNunapm199eYb32nen999IsxSceerlJFutpFuuXreqVaZLpzZlvXvoYk82SjZFTV1K1FrO99mLe5x8uVHh6p2j6N3V+JSkap2u8NyVS7CFuur1e6SkaRSfy5NsJiY0sQiaTzNvd42dLDoyGp2j6pKhmlar83JFPtImy5vl7pKhlFptRTb9e+VzMXVyKpPGVNVU09XyXaPo0pVdW3v1K1nFTL13aRr+uVSpE54s+1uyYbe8qaa70KIL3loFwpGeV6qUJyQ94f8dfoLZGg+18mJXrqU2OfbJRrvQogveWgdJWM0vEkKpHDlfeJP142nALWPu3/6J4ZzPniG3XOmw3xSvLS2U1V5HBEptQjuVkOSqeobp98Xe98Xa9UUOKPkFwsB6VTVLdPvq53vq5XKkQy8efakUCuxRt1uXIPg0RXJBN/NhwJeGUFf9/wOhAbjsErEz8cJBvileTlyj0MEl2RuribLbJpBEERiZ5IHvFnWrpuyxcRqYsSfwak67Z8EZG6ZCTxm9mFZva+mW0wsymZiCGTqm7Lj6cRBEUkXdKe+M2sALgfGAWcBow3s9PSHUcm6bZ8EcmkTFzcHQhscPcPAczsMWAs8G4GYsmI+Nvy97zxDEd/45u6LV9E0iYTif944OO495uBQbVnMrPJwOTgbbmZvV/HsjoCu1IeYbgOiXnnprczFErScm0751q8oJjTJaditnsOO94edU3M2u6c7v4A8EB985jZCncvTFNIKaGYw5dr8YJiTpdcizmseDNxcXcL8NW4992CaSIikgaZSPxvACeZWS8zawlcAczPQBwiIpGU9lKPux80sx8Ai4AC4EF3X9vExdVbCspSijl8uRYvKOZ0ybWYQ4nX3D2M5YqISJbSnbsiIhGjxC8iEjE5k/jNrI+ZrY772WNmN5vZdDPbEjf9ogzH+aCZ7TCzd+KmtTez581sffDvscF0M7P7gqEr3jKzM7Ik3n8zs/eCmJ4xs3bB9J5mti9uW/9XuuOtJ+aE7cDMpgbb+H0zK8qimB+Pi3ejma0Opmd8O5vZV81siZm9a2ZrzeymYHo2t+VEMWdte64n5nDbs7vn3A+xi8LbiN2cMB24LdMxxcV2LnAG8E7ctF8AU4LXU4B7gtcXASWAAWcCZVkS7wVA8+D1PXHx9oyfL8u2cZ3tgNiwIGuAI4BewF+AgmyIudbvfwX8NFu2M9AVOCN43Rb4INiW2dyWE8Wcte25nphDbc85c8Rfy3nAX9z9o0wHUpu7vwL8tdbkscDc4PVc4NK46Q97zGtAOzPrmpZAA3XF6+6L3f1g8PY1YvdaZI0E2ziRscBj7r7f3f8H2EBs2JC0qi9mMzNgHFCc1qDq4e5b3f3N4PXnwDpid91nc1uuM+Zsbs/1bOdEUtKeczXxX0HN/yQ/CE7jHqw69cwyXdx9a/B6G9AleF3X8BX17fRM+C6xI7kqvcxslZm9bGbnZCqoBOpqB7mwjc8Btrv7+rhpWbOdzawn0B8oI0facq2Y42Vte64j5tDac84lfovd9DUG+H/BpP8ETgT6AVuJnTJnLY+dr+VEH1ozuwM4CMwLJm0Furt7f+BW4FEzOzpT8dWSU+2glvHUPJDJmu1sZm2Ap4Cb3X1P/O+ytS0nijmb23MdMYfannMu8RMbzvlNd98O4O7b3b3C3SuB35OB0/gkbK867Q3+3RFMz9rhK8zsWuAS4MrgPzjB6eXu4PVKYvXFkzMWZJx62kHWbmMAM2sOfAt4vGpatmxnM2tBLBnNc/eng8lZ3ZYTxJzV7bmumMNuz7mY+GscHdWqI34TeOeQT2TefGBC8HoC8Gzc9GuCHhFnAp/FnUZnjJldCNwOjHH3v8dN72Sx5ylgZicAJwFZ8diwetrBfOAKMzvCzHoRi/n1dMdXj5HAe+6+uWpCNmzn4LrDbGCdu98b96usbcuJYs7m9lxPzOG250xe0W7sD9Aa2A0cEzftj8DbwFvBRuma4RiLiZ2aHSBWf5sEdABeBNYDLwDtg3mN2ENp/hKsQ2GWxLuBWB1xdfDzX8G83wbWBtPeBEZn0TZO2A6AO4Jt/D4wKltiDqbPAb5Xa96Mb2fgbGJlnLfi2sFFWd6WE8Wcte25nphDbc8askFEJGJysdQjIiKHQYlfRCRilPhFRCJGiV9EJGKU+EVEIkaJX0QkYpT4RYjd2Wlms4LX083stuD1HDO7LHhdamaFdXx2o5m9XfU7M2tuZv/XYkMXVw2re0cD3/+Qmd1Qa9qlZlZiZkcGy/jSzDqmap0lupT4RVJjuLuvCF7fDXwF+F/u3o/YIGwtGvh8MbHBB+NdARS7+75gOZ+kLlyJMiV+yRvBgzXeC47SPzCzeWY20syWBUffoY/jZGZHAdcD/+TuX0BsuF13nx43z1Vm9npwFP+7YNiAF4FT4sbBaU1sOIf/DjtmiR4lfsk3vYmNZHhK8PO/id0WfxvwkzR9/yaPja1+CDM7FbgcGBIcxVcQGzisgthAXeOCWUcDpV5rREyRVFDil3zzP+7+tsdGNVwLvOixcUneJvbEpbQys4nBkf3HZvZVYg8RGgC8YbFHLZ4HnBDMHl/uqf3MCZGUaZ7pAERSbH/c68q495Wkp71vALqbWdugxPMQ8JDFnrVbQGwws7nuPrWOz/4Z6GpmfYGzOLTmL5ISOuIXSSGPDfs7G5hlZq0Aghp+y2CWF4HLzKxz8Lv2ZtYj+KwTG5d/LlBSdY1AJNWU+EVS7w5iQzC/Y2argKXEkvkn7v4ucCew2MzeAp4n9sDtKsVAX1TmkRBpWGaRw2RmG4mNP78rH75H8p+O+EUO307gxbpu7kqFqhu4iN0LUBnGd0i06IhfIsXMJgI31Zq8zN1vzEQ8IpmgxC8iEjEq9YiIRIwSv4hIxCjxi4hEjBK/iEjE/H+0EiHUMk+pfQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "hep.histplot([H_bkg, H_125], bins=bins_ar, stack=True, label=[\"Backgrounds\", \"Higgs\"], histtype='fill')\n",
    "hep.histplot([H_data], bins=bins_ar, stack=False, yerr=True, histtype=\"errorbar\", color=\"black\",label=\"ATLAS Open Data\")\n",
    "\n",
    "plt.legend(loc=1, ncol=3, fontsize=9)\n",
    "plt.xlabel(\"m_llll [GeV]\")\n",
    "plt.ylabel(\"# events\")\n",
    "plt.ylim([0.0, 60])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1 \n",
    " Discuss if the model(coloured histograms) by the is a good model for the data (black points) and why.\n",
    " Calculation of statistics to make your answer more quantitative is encouraged."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You first must define a chi-squared function that will quantify how closely the model and data agree. It should take the observed bin heights of the data and  prediect bin heights of the model as inputs. The data can be assumed to be poisson-distributed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calcChiSq(obs, preds):\n",
    "    #REPLACE THIS FUNCTION WITH A VALID CHI-SQUARED CALCULATION\n",
    "    chiSq=1\n",
    "    ndf=1\n",
    "    return chiSq, ndf"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting s_s\n",
    "  We can make the crude assumption that the only thing wrong with our model is that the predicted total number of events from the Higgs signal is wrong by some single multiplicative factor which we will call s_s. We can investigate then how the agreement between the model and the data could be improve by applying various values of s_s and observing the change in the chi-squared. \n",
    "  \n",
    "  This code will not give any sensibel results without a valid chi-squared calculation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_s_ar = np.linspace(1.0, 3.0, 100) # array of s_s values we will investigate\n",
    "\n",
    "chi2_ar = np.empty( len(s_s_ar) ) # empty array of that will hold the chi2 values we will calculate\n",
    "\n",
    "minChi2 = 1000000\n",
    "bestFit_s_s = 0.0 #starting values for min chi2 and best value of s_s\n",
    "\n",
    "for s_s in range(0, len(s_s_ar)): #looping over out s_s values \n",
    "    pred = (s_s_ar[s_s]*H_125) + (H_bkg) # generating a prediction according to this s_s\n",
    "    chi2 = calcChiSq(H_data, pred)[0] # calculating chi2\n",
    "    chi2_ar[s_s] = chi2 # adding chi2 value to chi2 array \n",
    "    if(chi2 < minChi2): # check if this is the lowest chi2 we have seen so far\n",
    "        minChi2 = chi2 # update lowest chi2 value seen \n",
    "        bestFit_s_s = s_s_ar[s_s] # update value for best fit s_s\n",
    "        \n",
    "deltaChi2_ar = chi2_ar - minChi2 # make array of delta chi2 values"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming you have implmented a correct chi-squared calcualtion in the earlier cell. You should now have an array of delta chi2 values. The next cell will give an \"out of bounds error\" if you have not yet implmented your proper chi-"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[]\n"
     ]
    },
    {
     "ename": "IndexError",
     "evalue": "index 1 is out of bounds for axis 0 with size 0",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mIndexError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-15-36af91d7835e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m     14\u001b[0m \u001b[0;31m#shading in the confience interval band\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     15\u001b[0m \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcrit\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 16\u001b[0;31m \u001b[0mpx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcrit\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mcrit\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m0.001\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     17\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m \u001b[0mfig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0max\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubplots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIndexError\u001b[0m: index 1 is out of bounds for axis 0 with size 0"
     ]
    }
   ],
   "source": [
    "\n",
    "\n",
    "\n",
    "#we expect the chi2 vs. mZ curve to be quadratic, so let's fit that function to it.\n",
    "z = np.polyfit(s_s_ar, deltaChi2_ar, 2) #\"2\" for a second-order polynomial\n",
    "\n",
    "\n",
    "p = np.poly1d(z)\n",
    "\n",
    "# we can display the estimated uncertianty on mZ via critical values of the delta chi-squared curve\n",
    "y0 = 1.0\n",
    "crit = (p - y0).roots # roots of the polynominal -1, i.e., the mz values where p = 1 \n",
    "\n",
    "#shading in the confience interval band \n",
    "print(crit)\n",
    "px=np.arange(crit[1],crit[0],0.001)\n",
    "        \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(s_s_ar, deltaChi2_ar, 'k', linewidth=2)\n",
    "ax.fill_between(px,p(px),alpha=0.5, color='g', label=\"uncertainty\")\n",
    "\n",
    "plt.xlabel(\"s_s\")\n",
    "plt.ylabel(\"delta chi-squared\")\n",
    "plt.show()\n",
    "\n",
    "print(\"Best fit value of s_s = \" + str(round(bestFit_s_s, 3)) + \" +/- \" + str(round(np.abs(crit[0] - bestFit_s_s),2 )) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 2\n",
    "  Discuss the results of the fitting of s_s. Is the fitted model now a successful mode for the data? How can you quantify this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting s_b\n",
    "  Alternatively, we can make the crude assumption that the only thing wrong with our model is that the predicted total number of events from the *backgrounds* is wrong by some single multiplicative factor which we will call s_b.\n",
    " We can investigate then how the agreement between the model and the data could be improve by applying various values of s_s and observing the change in the chi-squared. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_b_ar = np.linspace(1.15, 1.4, 100)\n",
    "\n",
    "chi2_ar = np.empty( len(s_b_ar) )\n",
    "\n",
    "minChi2 = 1000000\n",
    "bestFit_s_b = 0.0\n",
    "\n",
    "for s_b in range(0, len(s_b_ar)):\n",
    "    pred = (H_125) + (s_b_ar[s_b]*H_bkg)\n",
    "    chi2 = calcChiSq(H_data, pred)[0]\n",
    "    chi2_ar[s_b] = chi2\n",
    "    if(chi2 < minChi2):\n",
    "        minChi2 = chi2\n",
    "        bestFit_s_b = s_b_ar[s_b]\n",
    "        \n",
    "deltaChi2_ar = chi2_ar - minChi2\n",
    "\n",
    "z = np.polyfit(s_b_ar, deltaChi2_ar, 2)\n",
    "p = np.poly1d(z)\n",
    "\n",
    "y0 = 1.0 # this is the value of the delta chi-squared function that defines the 68% CI for a one parameter fit.\n",
    "         # we'll invetigate if this parmaeter makes sense in the final exercise.\n",
    "crit = (p - y0).roots \n",
    "\n",
    "px=np.arange(crit[1],crit[0],0.001)\n",
    "        \n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(s_b_ar, deltaChi2_ar, 'k', linewidth=2)\n",
    "ax.fill_between(px,p(px),alpha=0.5, color='g', label=\"uncertainty\")\n",
    "plt.xlabel(\"s_b\")\n",
    "plt.ylabel(\"delta chi-squared\")\n",
    "plt.show()\n",
    "\n",
    "print(\"Best fit value of s_b = \" + str(round(bestFit_s_b, 3)) + \" +/- \" + str(round(np.abs(crit[0] - bestFit_s_b),2 )) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 3\n",
    "Discuss the results of the fitting of s_b. Is the fitted model now a successful mode for the data? Is it better than the previous s_s model? How can you quantify this?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fitting s_s and s_b\n",
    "More realistically, we can assume that the the predicted total numbers of events from both the Higgs signal AND the backgrounds are wrong by some two seperate multiplicative factors. We can investigate then how the agreement between the model and the data could be improve by applying various values of (s_s, s_b) and observing the change in the chi-squared. We now need to plot the delta chi-squared as a 2-D function of (s_s, s_b) as it depends on both parmaters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "s_s_ar = np.linspace(0.1, 2.7, 100)\n",
    "s_b_ar = np.linspace(1.1, 1.5, 100)\n",
    "\n",
    "chi2_ar = np.empty( (len(s_s_ar), len(s_b_ar) ))\n",
    "\n",
    "bestFit_s_s = 0.0\n",
    "bestFit_s_b = 0.0\n",
    "\n",
    "minChi2 = 1000000\n",
    "\n",
    "for s_s in range(0, len(s_s_ar)):\n",
    "    for s_b in range(0, len(s_b_ar)):\n",
    "        pred = (s_s_ar[s_s]*H_125) + (s_b_ar[s_b]*H_bkg)\n",
    "        chi2 = calcChiSq(H_data, pred)[0]\n",
    "        chi2_ar[s_s, s_b] = chi2 \n",
    "        if(chi2 < minChi2):\n",
    "            minChi2 = chi2\n",
    "            bestFit_s_s = s_s_ar[s_s]\n",
    "            bestFit_s_b = s_b_ar[s_b]\n",
    "\n",
    "deltaChi2_ar = chi2_ar - minChi2\n",
    "\n",
    "levels = [2.3] # this is the value of the delta chi-squared function that defines the 68% CI for a two parameter fit.\n",
    "         # we'll invetigate if this parmaeter makes sense in the final exercise.\n",
    "fig = plt.figure(figsize=(6,6))\n",
    "ax = fig.gca()\n",
    "\n",
    "cfset = plt.contourf(s_b_ar, s_s_ar, deltaChi2_ar,  cmap='coolwarm')\n",
    "cset = plt.contour(s_b_ar, s_s_ar, deltaChi2_ar, levels=levels, colors=['white'])\n",
    "\n",
    "bf_point = plt.scatter(bestFit_s_b, bestFit_s_s, color='black')\n",
    "\n",
    "\n",
    "cbar = plt.colorbar(cfset)\n",
    "cbar.set_label('Delta Chi-squared', fontsize=15, rotation=270)\n",
    "ax.clabel(cset, inline=1, fontsize=10)\n",
    "ax.set_xlabel('s_b',fontsize=18)\n",
    "ax.set_ylabel('s_s',fontsize=18)\n",
    "plt.show()\n",
    "\n",
    "#extract 2D result\n",
    "print(\"best-fit s_s = \" + str(bestFit_s_s))\n",
    "print(\"best-fit s_b = \" + str(bestFit_s_b))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 4\n",
    "Discuss the results of the simulatenous fit of s_s and s_b. How do the results of the three fitting approaches compare? If we interpret the Confidence Intervals on s_s and s_b as uncertainties on measurements of these parameters, then how do the uncertainties compare? Have we measured one parameter better than the other? If so, why? Does it make sense given the histograms we used to perform the fit of these parameters? How do the best-fit vales of the s_s and s_b comapre to the values when we fit one parameter at a time? Discuss what can be learned from the shape of the contours of the 2-D delta-chi-squared function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a 3D version of the 2-D contour Confidence Interval plot we made earlier, for no reason other than it looks cool and making cool plots is one of the main reasons data analysis is fun!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpl_toolkits import mplot3d\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig = plt.figure()\n",
    "ax = plt.axes(projection='3d')\n",
    "ax.contour3D(s_b_ar, s_s_ar, deltaChi2_ar, 80, cmap='coolwarm')\n",
    "ax.set_xlabel('s_b')\n",
    "ax.set_ylabel('s_s')\n",
    "ax.set_zlabel('delta chi-squared')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have estimates of the 'best-fit' values of s_s and s_b. Let's apply these factors to the predictions for signal and background and see how the updated predctions compare to data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = plt.figure()\n",
    "\n",
    "H_125_bf = bestFit_s_s*H_125\n",
    "H_bkg_bf = bestFit_s_b*H_bkg\n",
    "\n",
    "hep.histplot([H_bkg_bf, H_125_bf], bins=bins_ar, stack=True, label=[\"Backgrounds (best-fit)\", \"Higgs (best-fit)\"], histtype='fill')\n",
    "hep.histplot([H_data], bins=bins_ar, stack=False, yerr=True, histtype=\"errorbar\", color=\"black\",label=\"ATLAS Open Data\")\n",
    "plt.legend(loc=1, ncol=2, fontsize=9)\n",
    "plt.xlabel(\"m_llll [GeV]\")\n",
    "plt.ylabel(\"# events\")\n",
    "plt.ylim([0.0, 60])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data and our best-fit predictions agree very well... right?.......right??"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Goodness of fit\n",
    "We see that our fit succeeded it produced predictions for our signal and background that seem to agree much more closely with data than the original predictions. However, are we sure that our fitted model is a 'good' model for the data? We could as the question: \"If this is the correct model, then how probable is our observed data?\". This is answered by working our the p-value of our observed data acording to our fitted model. \n",
    "\n",
    "In the cell below, we take out fitted model and assume each bin is distrubuted as a gaussian with mean equal to the prediction in the bin and standard deviation equal to the square root of the mean. We generate N toy from the this fitted model and calculate the chi-squared for each toy with respect to the fitted model.\n",
    "\n",
    "Overlaying the chi-squared value of our ATLAS data with respect to our model \"If this is the correct model, then how probable is our observed data?\". This is answered by working our the p-value of our observed data acording to our fitted model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# generate the chi-2 distribution of N toy experiments based on the fitted model \n",
    "# assume a gaussian-distributed bin height with mean = model pred. and sigma equal root model pred\n",
    "# what then is the p-value of the data with respect to the distribution?\n",
    "    \n",
    "#function to generate a toy histogram given the mean for each bin and assuming the data is\n",
    "# gaussian-distrubtu\n",
    "def generate_toy(means): \n",
    "        toy = np.empty(len(means))\n",
    "        for i in range(0, len(means)):\n",
    "            toy_bin = np.random.normal(means[i], np.sqrt(means[i]), 1)\n",
    "            toy[i] = toy_bin\n",
    "        #print(\"toy = \" + str(toy))\n",
    "        return toy\n",
    "\n",
    "ntoys = 10000 # start with a small number, increase when you undertsand your results\n",
    "\n",
    "means = H_bkg_bf + H_125_bf # take the results of the fit as the means of the fitted model\n",
    "\n",
    "#print(\"means\" + str(means))\n",
    "\n",
    "chi2_toys = np.empty(ntoys)\n",
    "\n",
    "for t in range(0, ntoys):\n",
    "    toy = generate_toy(means) \n",
    "    chi2_toys[t] = calcChiSq(toy, means)[0]\n",
    "    #print(\"chi2 \" + str(chi2_toys[t]))\n",
    "    \n",
    "plt.figure()\n",
    "fig, ax1 = plt.subplots()\n",
    "\n",
    "#plot distribtion of chi-squared values from toys\n",
    "bins_ar = np.linspace(0.0, 200, num=(nBins+1))\n",
    "chi2Hist, chi2bins = np.histogram(chi2_toys, bins=bins_ar, density=True)                        \n",
    "hep.histplot([chi2Hist], bins=chi2bins, histtype='fill', label=\"toys\")\n",
    "\n",
    "#plot expected distrutuon of chi-squared values from theory - chi-squared distribtuions with ndof = nbins -2\n",
    "df = nBins - 2\n",
    "x = np.linspace(stats.chi2.ppf(0.000001, df), stats.chi2.ppf(0.999999, df), 200)\n",
    "label = '$\\chi^{2}$ (ndof = ' + str(df) + ') pdf'\n",
    "ax1.plot(x, stats.chi2.pdf(x, df), 'r-', lw=5, alpha=0.6, label=label)\n",
    "\n",
    "# overlay the chi2 squared value from the ATLAS OpenData\n",
    "chi2_data = calcChiSq(H_data, means)[0]\n",
    "plt.plot(chi2_data,0.002 , 'y*', ms=14, label='ATLAS OpenData')\n",
    "plt.legend()\n",
    "plt.xlabel(\"chi^2\")\n",
    "plt.ylabel(\"# toys\")\n",
    "plt.xlim(0,100)\n",
    "\n",
    "print(\"chi2 value of ATLAS Open data is = \" + str(chi2_data))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 5 - Goodness of fit and p-value\n",
    "Calculate the p-value of the ATLAS data with respect to the chi-squared distribution shown above. Comment if this value indicates if the fitted model is a good model for the ATLAS data or not. Repeat the procedure to calculate the p-value of the data with respect to the model **before** any fitting occured. How do the two p-values compare? What do you conclude about the fitted model? \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Exercise 6 - How do we know the critical value for a Confidence Interval?\n",
    "Earlier we used critical value of 2.3 of the delta chi-squared function to define the 68% Confidence Interval on our fit parameters. This value of 2.3 is based on the propertires of the Gaussian distribution. In this exercise we aim to check that this value of 2.3 indeed gives us the correct 68 % CL. Recall the definition of a X % CL: *an interval constructed such that if we constructed this interval in the same way for a large number of repeated experiments, it would contain the true value of the parameter in X % of those experiments*.\n",
    "\n",
    "Use random numbers to work out the critical value of the delta chi-squared corresponding to the boundaries of the 68 % Confidence Interval for the fit performed above. Some hints are given towards the steps needed in the cell below. The steps are not neccessarily in the correct order, they rather they repreent the crucial pieces of information you need to build a solution.\n",
    "\n",
    "N.B. your setup will likely not initialy give you exaclty the correct critical values with respect to the theoretical expectation? What needs to be done such that the exact critical values can be determined?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# suggested steps\n",
    "\n",
    "    # generate N toys\n",
    "\n",
    "        # for each toy, run 2-parameter fit and estimate CL based on critical value d\n",
    "        # (note that you will fit the toy to the 'true' distributions, not real data)\n",
    "\n",
    "        # check if true value is within the contour of the delta chiquared function defined by d\n",
    "\n",
    "    # calculate fraction of toys (f) in which the CL estimation for critical value\n",
    "    # contains true values of s_s and s_b\n",
    "\n",
    "    # if you can check the fraction of toys with CI containing the true value \n",
    "    # or a given d, you just need to repeat this  for a range of d values.\n",
    "    # this will take some time, so start with a small number of toys and a few \n",
    "    # d values, when you understand your results you can increse both.\n",
    "\n",
    "    # plot f as a function of d\n",
    "\n",
    "    # read form this plot the critical value that gives you the 68% CI.\n",
    "    \n",
    "# Best of luck!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}