{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Moderne Methoden der Datenanalyse SS2022\n",
    "# Practical Exercise 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4: Combination of correlated measurements\n",
    "\n",
    "A common problem in science is the combination of several measurements to one single result, e.g., the average value. Not only the uncertainties of the individual measurements have to be taken into account, but also the correlations between them. A wrong treatment of correlations or common systematic effects can lead to biased results."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.1: Combination of *W* mass measurements (voluntary)\n",
    "\n",
    "At the LEP accelerator at CERN the mass of the *W* boson $m_W$ was\n",
    "  measured in two different channels:\n",
    "  \n",
    "  \\begin{eqnarray*}\n",
    "    e^+e^- &\\rightarrow \\ \\  W^+W^- \\  \\rightarrow& q_1 q_2 q_3 q_4 \\\\\n",
    "    e^+e^- &\\rightarrow \\ \\  W^+W^- \\  \\rightarrow& l \\nu q_1 q_2 \n",
    "  \\end{eqnarray*}\n",
    "  \n",
    "  The experimental signature in the detector for the first channel with four quarks are four reconstructed jets. The second channel is identified by a lepton (electron or muon) and two jets.\n",
    "  The neutrino is not detected. The measured *W* masses and the uncertainties are:\n",
    "  \n",
    "  \\begin{eqnarray*}\n",
    "    \\mbox{4 jets channel:} & \\,m_W = & \n",
    "    (80457 \\pm 30 \\pm 11 \\pm 47 \\pm 17 \\pm 17) \\mbox{ MeV} \\\\\n",
    "    \\mbox{lepton + 2 jets channel:} & \\,m_W = & \n",
    "    (80448 \\pm 33 \\pm 12 \\pm \\ 0 \\pm 19 \\pm 17) \\mbox{ MeV}\n",
    "  \\end{eqnarray*}\n",
    "  \n",
    "  To facilitate the interpretation of the results, different uncertainties are given, originating from different sources: The first two uncertainties are the statistical and systematic experimental\n",
    "  uncertainties, which are uncorrelated. The third uncertainty originates from the theory applied for the analysis and is only present in the first channel. The fourth uncertainty comes from a common theoretical model applied for both channels, and thus is 100\\% correlated. Also the last uncertainty is 100\\% correlated between both measurements, since it represents the uncertainty on the LEP accelerator beam energy.\n",
    "\n",
    "\n",
    "  - Construct a covariance matrix of the two *W* mass measurements taking into account all uncertainties and their correlations. Use this covariance matrix to define a $\\chi^2$ expression containing the average *W* mass $\\bar{m}_W$ as a free parameter. Determine $\\bar{m}_W$ and its uncertainty by minimizing the $\\chi^2$ expression using the IMinuit package.\n",
    "\n",
    "  For this exercise, you have to write your own $\\chi^2$-function to be minimized; see the previous exercise to learn how this can be done. There we have already used `iminuit`, however this time we'll be using the object-oriented interface. You can find some hints on how to use this in the documentation [here](https://iminuit.readthedocs.io/en/stable/notebooks/basic_tutorial.html#Initialize-the-Minuit-object).\n",
    "  \n",
    "*Hint*: Make sure, that you are using a current version of `iminuit`, e.g. 2.11.2! Check this with the command in the next cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip list | grep iminuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from iminuit import Minuit\n",
    "from scipy.linalg import inv\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi_square_creator(measurement_vector, inv_cov):\n",
    "    \n",
    "    # Using this outer function to pass the measured values (measurement_vector) and the inverated covariance matrix (inv_cov) to the function so we don't need to define them globally\n",
    "    \n",
    "    def chi_square_function(par):\n",
    "\n",
    "        chi2_value = 0\n",
    "\n",
    "        # TODO: Add code to calculate chi2\n",
    "  \n",
    "        return chi2_value # return the chi2 value\n",
    "\n",
    "    return chi_square_function # return the function which calculates the value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "jp-MarkdownHeadingCollapsed": true,
    "tags": []
   },
   "source": [
    "As this is using a different interface to IMinuit than in Exercise 3.2, here is a hint: These are the base two lines you'll need:\n",
    "\n",
    "```\n",
    "minuit_instance = Minuit(function, parametername=initial_parametervalue)\n",
    "res = minuit_instance.migrad()\n",
    "```\n",
    "\n",
    "The first line initializes the minuit object, gives the cost function and passes the initial parameter names and values.\n",
    "The second line performs the minimization using the `MIGRAD` algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Add code here to calculate the (inverse) covariance matrix\n",
    "# You can use scipy.linalg.inv to invert the matrix\n",
    "\n",
    "# Perform the minimization of the chi2 function\n",
    "\n",
    "# Get results of the minimization and plot or print them"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  - Because the minimization of the $\\chi^2$ expression in exercise 4.1 is a linear problem it can be solved analytically.  Determine\n",
    "  $\\bar{m}_W$ and its error analytically and compare them to the result from above.\n",
    "  \n",
    "**TODO:** Add your calculations here using the Latex syntax\n",
    "\n",
    "  - Estimate the contributions from statistical, systematic, theoretical, and accelerator based uncertainties to the error of the combined *W* mass measurement.  Use the quadratic difference between the total error and the error calculated with a covariance matrix where one component is removed.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Add your code here to calculate the contribution of each uncertainty component"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 4.2: Normalisation uncertainty (obligatory)\n",
    "\n",
    "Two measurements $y_1=8.0$ and $y_2=8.5$ of the same physical\n",
    "  quantity with an uncorrelated relative statistical error of 2 %\n",
    "  and a common normalisation error of 10 % should be combined.\n",
    "  \n",
    "  -  Construct a covariance matrix and a $\\chi^2$ expression and\n",
    "  determine its minimum with `iminuit` or analytically. If you need hints on how to use `iminuit`, consult the exercises 3.2 and 4.1. \n",
    "\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ROOT # Only need for plotting... Feel free to replace the plot method with your own pure python implementation!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "   \n",
    "def chi2_creator_4_2_1(measurement_vector, inv_cov):\n",
    "    \n",
    "    # We are using this outer function to pass the measured values (measurement_vector) and the inverated covariance matrix (inv_cov)\n",
    "    # to the function so we don't need to define them globally...\n",
    "\n",
    "    def chi2_function(par): \n",
    "        '''\n",
    "        calculate chi2 using 1 parameter for the mean\n",
    "        '''\n",
    "        chi2_value = 0\n",
    "        \n",
    "        # TODO: Calculate chi2 here\n",
    "        \n",
    "        return chi2_value\n",
    "    \n",
    "    return chi2_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " - Is the result reasonable?  What could be the cause\n",
    "  for the unexpected value?  Make a plot of the covariance ellipse in\n",
    "  the $y_1^\\prime y_2^\\prime$ plane defined by\n",
    "  $$\n",
    "  \\Delta y^T V^{-1} \\Delta y = c^2, \\quad \\Delta y =\n",
    "  \\left( \\begin{array}{c} y_1-y_1^\\prime \\\\\n",
    "      y_2-y_2^\\prime \\end{array}\\right)\n",
    "  $$\n",
    "  for $c=1$ and $c=2$ together with the line $y_1^\\prime=y_2^\\prime$.\n",
    "  $V$ is the covariance matrix.  To draw the ellipse a TGraph\n",
    "  object can be used.  The points on the ellipse can be calculated as\n",
    "  a function of the angle $\\phi$ if $\\Delta y$ is expressed by $\\phi$\n",
    "  and the radius $r$.  You can use the function below to draw the ellipse.\n",
    "  Pay attention to where the bisector intersects with the ellipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drawCovEllipse(CI: list, vec: list, mean: float):\n",
    "    '''\n",
    "    CI : inverse covariance matrix\n",
    "    vec: measured values i.e. [y1, y2]\n",
    "    mean: calculated mean value from chi2 minimization\n",
    "        \n",
    "    Draws and returns the covariance ellipses as well as the bisect line.\n",
    "    The outputs have to be stored in some variables or they will not be displayed.\n",
    "    '''\n",
    "    \n",
    "    npoints = 200\n",
    "    ellipse1 = ROOT.TGraph(npoints + 1)\n",
    "    ellipse2 = ROOT.TGraph(npoints + 1)\n",
    "    for i in range (npoints + 1):\n",
    "        phi = 2 * i * np.pi / npoints\n",
    "        V = [np.cos(phi), np.sin(phi)]\n",
    "        v = np.matrix(V)\n",
    "        sp = np.dot(v, np.dot(CI, v.getT()))\n",
    "        R = np.sqrt(1.0 / sp)\n",
    "        ellipse1.SetPoint(i, vec[0] + R * V[0], vec[1] + R * V[1])\n",
    "        ellipse2.SetPoint(i, vec[0] + 2 * R * V[0], vec[1] + 2 * R * V[1])\n",
    "        \n",
    "    ellipse2.SetLineWidth(2)\n",
    "    ellipse1.SetLineWidth(2)\n",
    "    \n",
    "    ellipse2.Draw(\"AL\")\n",
    "    ellipse1.Draw(\"Ls\")\n",
    "    \n",
    "    x0 = ellipse2.GetXaxis().GetXmin()\n",
    "    y0 = ellipse2.GetYaxis().GetXmin()\n",
    "    x1 = ellipse2.GetXaxis().GetXmax()\n",
    "    y1 = ellipse2.GetYaxis().GetXmax()\n",
    "    \n",
    "\n",
    "    line = ROOT.TLine(max(x0, y0), max(x0, y0), min(x1, y1), min(x1, y1))\n",
    "    line.SetLineWidth(2)\n",
    "    line.SetLineColor(4)\n",
    "    line.Draw('same')\n",
    "    marker = ROOT.TMarker(mean, mean, 2)\n",
    "    marker.SetMarkerColor(2)\n",
    "    marker.SetMarkerSize(3)\n",
    "    marker.Draw('same')\n",
    "    \n",
    "    ROOT.gPad.Update()\n",
    "    \n",
    "    return ellipse1, ellipse2, line, marker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: Draw the error ellipse\n",
    "# Hint: you have to assign the return values of the function to some variables, otherwise they will not be shown in the plot. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "  - Use an additional normalisation parameter $N$ for the treatment of\n",
    "  the common normalisation uncertainty $\\sigma_{norm}$ instead of\n",
    "  taking it into account in the covariance matrix of $y_1$ and $y_2$.\n",
    "  Add a term to the $\\chi^2$ expression for the normalisation with an\n",
    "  expected value of 1 and an error of 10 %.  The normalisation\n",
    "  factor $N$ can be applied either to the measured values $y_i$ or to\n",
    "  the fit parameter $\\bar{y}$.  Try out both ways (using\n",
    "  `iminuit` for the $\\chi^2$ minimisation) and compare the\n",
    "  results.  Which one is the more meaningful result and why?\n",
    "  \n",
    "  Determine $\\bar{y}$ from the correct $\\chi^2$ expression in an analytical way.  How does the normalisation error affect\n",
    "  the averaged value and its error?\n",
    "  \n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2_creator_4_2_2(measurement_vector, inv_cov, sigma_norm):\n",
    "   \n",
    "    # Now we additionally have to set the normalisation uncertainty sigma_norm in the cost function.\n",
    "\n",
    "    def chi2_function(par):\n",
    "        '''\n",
    "        calculate the chi2 including an additional parameter for the normalization of the measured values (version 1)\n",
    "        '''\n",
    "        # TODO: Calculate chi2 here. Apply the normalization factor to the measured values.\n",
    "        return chi2_value\n",
    "    \n",
    "    return chi2_function\n",
    "\n",
    "def chi2_creator_4_2_3(measurement_vector, inv_cov, sigma_norm):\n",
    "    \n",
    "    def chi2_function(par):\n",
    "        '''\n",
    "        calculate the chi2 including an additional parameter for the normalization of the mean value (version 2)\n",
    "        '''\n",
    "        # TODO: Calculate chi2 here. Apply the normalization factor to the fit parameters.\n",
    "        return chi2_value\n",
    "    \n",
    "    return chi2_function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "  - Construct a covariance matrix of $y_1$ and $y_2$ containing the\n",
    "  normalisation uncertainty of 10\\,\\% relative to the average value\n",
    "  $\\bar{y}$.  Solve the corresponding $\\chi^2$ minimisation with\n",
    "  `iminuit` and plot the covariance ellipse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def chi2_creator_4_2_4(measurement_vector, cov):\n",
    "\n",
    "    # Here we pass only the non-inverted covariance matrix as the inverted matrix should depend on the normalisation. This means we have to do the matrix inversion inside chi2_function(par).\n",
    "    \n",
    "    def chi2_function(par):\n",
    "        '''\n",
    "        calculate chi2 with normalization error relative to fitted mean\n",
    "        '''\n",
    "        # TODO: Calculate chi2 here\n",
    "        return chi2_value\n",
    "\n",
    "    return chi2_function"
   ]
  }
 ],
 "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}