lecture001-theory-and-basic-tracking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Instructions on Running this Slide Deck in Jupyter\n",
    "\n",
    "I recommend the Anaconda-based approach:\n",
    "\n",
    "* [Install Anaconda (Python 3.X) from the Continuum.io Website](https://www.anaconda.com/distribution/)\n",
    "* Enable your Anaconda environment\n",
    "* Install some other packages:\n",
    "   * Jupyter: ```conda install jupyter```\n",
    "   * [Notebook extensions](https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/install.html): ```conda install -c conda-forge jupyter_contrib_nbextensions```\n",
    "   * [RISE](https://damianavila.github.io/RISE/index.html): ```conda install -c damianavila82 rise```\n",
    "   * [PyHF](https://github.com/diana-hep/pyhf) for data modeling demonstration:\n",
    "      * ```conda install pip3```\n",
    "      * ```pip install pyhf```\n",
    "* Launch the slides:\n",
    "   * ```jupyter-notebook <SLIDES.ipynb>```\n",
    "   * Use \"Alt-r\" to go to slideshow mode and Spacebar (SHIFT+Spacebar) to go forward (backward) through the slideshow. [More info is available in the RISE documentation](https://damianavila.github.io/RISE/usage.html).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# %load units.py\n",
    "# Conversion Factors\n",
    "import math\n",
    "\n",
    "# distance\n",
    "_m=1\n",
    "_km=1000*_m\n",
    "_cm=_m/100\n",
    "_mm=_m/1000.\n",
    "_in=2.54*_cm\n",
    "_ft=12*_in\n",
    "_mile=5280*_ft\n",
    "_miles=_mile\n",
    "_ly=9460730472580800*_m\n",
    "_Angstrom=1e-10*_m\n",
    "\n",
    "\n",
    "# time\n",
    "_s=1\n",
    "_ms=_s/1000\n",
    "_mus=_s/1e6\n",
    "_min=60*_s\n",
    "_hour=60*_min\n",
    "_hours=_hour\n",
    "_day=24*_hours\n",
    "_days=_day\n",
    "_year=365.25*_days\n",
    "_years=_year\n",
    "_year_approx=365*_days\n",
    "\n",
    "# mass\n",
    "_kg=1\n",
    "_g=_kg/1000\n",
    "_lb=_kg/2.2\n",
    "_oz=(1/16)*_lb\n",
    "_amu = 1.66053886e-27*_kg\n",
    "\n",
    "# charge\n",
    "_C = 1.602e-19 # couloumb\n",
    "\n",
    "# Force\n",
    "_N = _kg*_m/_s**2\n",
    "_dyn = 1e-5*_N\n",
    "\n",
    "# Energy\n",
    "_J = _N*_m\n",
    "_W = _J/_s\n",
    "\n",
    "_eV = _J*_C\n",
    "_keV = 1e3*_eV\n",
    "_MeV = 1e6*_eV\n",
    "_GeV = 1e9*_eV\n",
    "_TeV = 1e12*_eV\n",
    "\n",
    "\n",
    "# Temperature\n",
    "\n",
    "_K = 1\n",
    "_Celsius = _K\n",
    "\n",
    "# Angles\n",
    "_rad=1\n",
    "_deg=math.pi/180*_rad\n",
    "\n",
    "# Constants of nature\n",
    "_c = 2.998e8 * _m/_s\n",
    "_h = 6.626e-34 * _kg*_m/_s\n",
    "_hbar = _h/(2*math.pi)\n",
    "_G = 6.67408e-11 * _N*_m**2/_kg**2\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0,
     2
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# %load helpers.py\n",
    "# Rounding to specific precision, from https://github.com/randlet/to-precision\n",
    "def to_precision(x,p):\n",
    "    \"\"\"\n",
    "    returns a string representation of x formatted with a precision of p\n",
    "\n",
    "    Based on the webkit javascript implementation taken from here:\n",
    "    https://code.google.com/p/webkit-mirror/source/browse/JavaScriptCore/kjs/number_object.cpp\n",
    "    \"\"\"\n",
    "\n",
    "\n",
    "    import math\n",
    "    x = float(x)\n",
    "\n",
    "    if x == 0.:\n",
    "        return \"0.\" + \"0\"*(p-1)\n",
    "\n",
    "    out = []\n",
    "\n",
    "    if x < 0:\n",
    "        out.append(\"-\")\n",
    "        x = -x\n",
    "\n",
    "    e = int(math.log10(x))\n",
    "    tens = math.pow(10, e - p + 1)\n",
    "    n = math.floor(x/tens)\n",
    "\n",
    "    if n < math.pow(10, p - 1):\n",
    "        e = e -1\n",
    "        tens = math.pow(10, e - p+1)\n",
    "        n = math.floor(x / tens)\n",
    "\n",
    "    if abs((n + 1.) * tens - x) <= abs(n * tens -x):\n",
    "        n = n + 1\n",
    "\n",
    "    if n >= math.pow(10,p):\n",
    "        n = n / 10.\n",
    "        e = e + 1\n",
    "\n",
    "\n",
    "    m = \"%.*g\" % (p, n)\n",
    "\n",
    "    if e < -2 or e >= p:\n",
    "        out.append(m[0])\n",
    "        if p > 1:\n",
    "            out.append(\".\")\n",
    "            out.extend(m[1:p])\n",
    "        out.append('e')\n",
    "        if e > 0:\n",
    "            out.append(\"+\")\n",
    "        out.append(str(e))\n",
    "    elif e == (p -1):\n",
    "        out.append(m)\n",
    "    elif e >= 0:\n",
    "        out.append(m[:e+1])\n",
    "        if e+1 < len(m):\n",
    "            out.append(\".\")\n",
    "            out.extend(m[e+1:])\n",
    "    else:\n",
    "        out.append(\"0.\")\n",
    "        out.extend([\"0\"]*-(e+1))\n",
    "        out.append(m)\n",
    "\n",
    "    return \"\".join(out)\n",
    "\n",
    "def latexify_sci(number, template=\"%.2f \\\\times 10^{%d}\"):\n",
    "    if type(number) is str:\n",
    "        number_string = number\n",
    "    else:\n",
    "        number_string = \"%.6e\" % (number)\n",
    "        pass\n",
    "    [prefix, power10] = number_string.split(\"e\")\n",
    "\n",
    "    latex_code = \"$\" + template % (float(prefix), int(power10)) + \"$\"\n",
    "    return latex_code\n",
    "\n",
    "def printy(number, decimals=2):\n",
    "    return latexify_sci(to_precision(number,decimals))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Jet Flavor Tagging - Part 1\n",
    "\n",
    "## Heavy Quark States and Charged Particle Reconstruction\n",
    "\n",
    "## Stephen Sekula (SMU)\n",
    "\n",
    "Presented at the SMU ATLAS Weekly Meeting\n",
    "\n",
    "June 19, 2019\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Overview\n",
    "\n",
    "* Jets and Heavy Flavor: Theory Aspects\n",
    "* Heavy Flavor: Experimental Knowns\n",
    "* Heavy Flavor Physics and Implications for Experiments\n",
    "* First Steps: Tracking a Charged Particle\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Jets and Heavy Flavor: Theory Aspects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# ATLAS Detector: A $\\mathbf{t \\bar{t}}$ Candidate Event\n",
    "\n",
    "<center>\n",
    "<img src=\"Images/run267638_evt193690558_ttbar_cropped.png\" style=\"\"/>\n",
    "</center>\n",
    "A $t\\bar{t}$ event should contain at least 2 bottom-quark-initiated jets. How might we spot such jets?   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Heavy Flavor: Experimental Knowns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Properties of Specific Heavy Flavor States\n",
    "\n",
    "* b-hadrons [[PDG]](https://pdg.lbl.gov)\n",
    "   * $B^0$ ($\\bar{b}d$): mass: $5279.64 \\pm 0.13$ MeV ; lifetime: $1.519 \\pm 0.004$ ps\n",
    "   * $B^{+}$ ($\\bar{b}u$): mass: $5279.33 \\pm 0.13$ MeV ; lifetime: $1.638 \\pm 0.004$ ps\n",
    "   * $\\Lambda_b^0$ ($bdu$): mass: $5619.60 \\pm 0.17$ MeV ; lifetime: $1.471 \\pm 0.009$ ps\n",
    "   * Commonalities: heavy ($m = [5000,6000]\\mathrm{MeV}$) and live about 1-2 ps.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "* c-hadrons\n",
    "   * $D^0$ ($\\bar{c}u$): mass: $1864.83 \\pm 0.05$ MeV; lifetime: $0.4101 \\pm 0.0015$ ps\n",
    "   * $D^{+}$ ($\\bar{c}d$): mass: $1869.65 \\pm 0.05$ MeV; lifetime: $1.040 \\pm 0.007$ ps\n",
    "   * Commonalities: moderately heavy ($m \\approx 2000\\mathrm{MeV}$) and live about 0.5-1 ps\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "variables": {
     "round((1-1/math.e)*100)": "63"
    }
   },
   "source": [
    "   \n",
    "Reminder: particle decay is *exponential* in nature, e.g. $P(h \\to X|t) = 1- e^{-t/\\tau}$; $t$ is the time after production; the lifetime, $\\tau$, tells you the time by which there is a $\\approx$ {{round((1-1/math.e)*100)}}% probability the particle (in its rest frame) has decayed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5665000000000002e-12\n",
      "8.300333333333334e-13\n",
      "5.682205604927609\n"
     ]
    }
   ],
   "source": [
    "# Particle lifetimes, flight lengths, etc.\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "m_B0 = 5.27964 # GeV\n",
    "m_Dp = 1.86965 # GeV\n",
    "tau_B0 = 1.519e-12 # s\n",
    "tau_Bp = 1.638e-12 # s\n",
    "tau_Lb0 = 1.471e-12 # s\n",
    "tau_D0 = 0.4101e-12 # s\n",
    "tau_Dp = 1.040e-12 # s\n",
    "\n",
    "# Average heavy flavor lifetime\n",
    "tau_b_avg = np.mean([tau_B0, tau_Bp, tau_Bp, tau_Lb0])\n",
    "tau_c_avg = np.mean([tau_D0, tau_Dp, tau_Dp])\n",
    "\n",
    "print(f\"{tau_b_avg}\")\n",
    "print(f\"{tau_c_avg}\")\n",
    "\n",
    "ctau_b_avg = _c * tau_b_avg\n",
    "ctau_c_avg = _c * tau_c_avg\n",
    "\n",
    "# example \n",
    "p_B0 = 30 # GeV\n",
    "betagamma_B0 = p_B0/m_B0\n",
    "print(f\"{betagamma_B0}\")\n",
    "bgct_B0 = betagamma_B0 * _c * tau_B0\n",
    "\n",
    "p_Dp = 30 # GeV\n",
    "betagamma_Dp = p_Dp/m_Dp\n",
    "bgct_Dp = betagamma_Dp * _c * tau_Dp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Properties of Specific Heavy Flavor States (cont.)\n",
    "\n",
    "* b-hadrons\n",
    "   * Commonalities: heavy ($m = [5000,6000]\\mathrm{MeV}$) and live about 1-2 ps.\n",
    "* c-hadrons\n",
    "   * Commonalities: moderately heavy ($m \\approx 2000\\mathrm{MeV}$) and live about 0.5-1 ps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "variables": {
     "to_precision(ctau_b_avg/_mm,3)": "0.470",
     "to_precision(ctau_c_avg/_mm,3)": "0.249"
    }
   },
   "source": [
    "* How far will a typical $b$ or $c$ hadron travel?\n",
    "   * Estimate using $c\\tau$\n",
    "   * For b-hadrons, the $c\\tau_b \\approx$ {{to_precision(ctau_b_avg/_mm,3)}} mm\n",
    "   * For c-hadrons, the $c\\tau_c \\approx$ {{to_precision(ctau_c_avg/_mm,3)}} mm\n",
    "   * Given the actual energy (E) and momentum (p) of a hadron, scale decay length using $\\beta \\gamma = (p/E)(E/m) = p/m$: $c\\tau \\to \\beta \\gamma c\\tau$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "variables": {
     "to_precision(bgct_B0/_mm,2)": "2.6"
    }
   },
   "source": [
    "   * Example: a $B^0$ meson with $p=30\\mathrm{GeV}$ travels {{to_precision(bgct_B0/_mm,2)}} mm, on average, before decaying.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    },
    "variables": {
     "to_precision(bgct_Dp/_mm,2)": "5.0"
    }
   },
   "source": [
    "   * Example: a $D^{+}$ meson with $p=30\\mathrm{GeV}$ travels {{to_precision(bgct_Dp/_mm,2)}} mm, on average, before decaying.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# More computations of flight lengths, etc.\n",
    "m_pi0 = 134.9770*_MeV\n",
    "m_rho0 = 775.26*_MeV\n",
    "m_K0s = 497.611*_MeV\n",
    "tau_pi0 = 8.6e-17*_s\n",
    "gamma_rho0 = 147.8*_MeV\n",
    "tau_rho0 = _hbar/gamma_rho0\n",
    "tau_K0s = 8.954e-11*_s\n",
    "ctau_pi0 = _c*tau_pi0\n",
    "ctau_rho0 = _c*tau_rho0\n",
    "ctau_K0s = _c*tau_K0s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    },
    "variables": {
     "latexify_sci(ctau_pi0*1000)": "$2.58 \\times 10^{-5}$",
     "latexify_sci(ctau_rho0*1000)": "$1.34 \\times 10^{-12}$",
     "latexify_sci(tau_K0s/1e-12)": "$8.95 \\times 10^{1}$",
     "latexify_sci(tau_rho0/1e-12)": "$4.45 \\times 10^{-12}$",
     "to_precision(ctau_K0s*1000,2)": "27"
    }
   },
   "source": [
    "# Compare to Some Light(er) Flavor States\n",
    "\n",
    "* b-hadrons\n",
    "   * Commonalities: heavy ($m = [5000,6000]\\mathrm{MeV}$) and live about 1-2 ps.\n",
    "* c-hadrons\n",
    "   * Commonalities: moderately heavy ($m \\approx 2000\\mathrm{MeV}$) and live about 0.5-1 ps\n",
    "* light hadrons\n",
    "   * $\\pi^0$: mass: $135$ MeV; lifetime: $8.6 \\times 10^{-5}$ ps $\\to$ $c\\tau_{\\pi^0}=$ {{latexify_sci(ctau_pi0*1000)}} mm\n",
    "   * $\\rho^0$: mass $775$ MeV; lifetime: {{latexify_sci(tau_rho0/1e-12)}} ps $\\to$ $c\\tau_{\\rho^0}=$ {{latexify_sci(ctau_rho0*1000)}} mm\n",
    "   * $K^0_s$: mass $498$ MeV; lifetime: {{latexify_sci(tau_K0s/1e-12)}} ps $\\to$ $c\\tau_{K^0_s}=$ {{to_precision(ctau_K0s*1000,2)}} mm\n",
    "\n",
    "Some light-flavour mesons (e.g. Neutral Kaons) could create the appearance of a long-lived particle inside a jet (because they are long-lived), albeit with less mass per hadron.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Summary: What to Expect from Hadrons in Jets\n",
    "\n",
    "* Many hadrons are \"long-lived\" (e.g. $K_s^0$, $B^0$)\n",
    "* Jets initiated by light quarks (u,d,s) will *typically* contain hadrons that decay promptly, and are typically low in mass (e.g. $m_{LF} < 1500$ MeV)\n",
    "* Jets initiated by heavy quarks (c,b) will *typically* contain hadrons that decay after a longer period of time, and are typically heavy in mass (e.g. $m_{HF} > 1500$ MeV)\n",
    "\n",
    "A strategy for identifying heavy-flavor-initiated jets will, in a brief sketch...\n",
    "* recognize that jet content will typically be heavier in mass\n",
    "* will contain \"substructure\" inside the hadronization/shower evolution that indicates a particle was produced, flew, and decayed some time later\n",
    "   * as a result of this feature, a particle detector for tracking charged particles will be crucial, as one will want to look for (a) the place where the jet itself originated and (b) the place where longer-lived particles later decay (ideally into charged particles)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Identifying the trajectory of charged particles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A Tracking Example: A Cosmic Ray\n",
    "\n",
    "<img src=\"Images/JiveXML_246892_4472609.png\" style=\"width: 50%; height: auto; float: right;\"/>\n",
    "\n",
    "A single cosmic-ray muon passes near the center of ATLAS. It's a great way, in a clean environment, to understand how we hope for tracks to look.\n",
    "\n",
    "We can also see how ATLAS originally deployed an inner tracking system to do this work. The Pixel Detector provides (originally) the inner 4 layers (the Insertable B-Layer, or IBL, was added for Run 2 and is here seen as the closest layer to the beamspot). The Silicon Tracker (SCT) provides the next 4 layers.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A Tracking Example: A Higher \"Pileup\" Event\n",
    "\n",
    "<img src=\"Images/JiveXML_266904_25884805-YX-RZ-EventInfo-2015-06-04-09-21-33_lowres.jpg\" style=\"width: 50%; height: auto; float: right;\"/>\n",
    "\n",
    "A more realistic representation of recent running conditions. Multiple proton-proton interactions have happened, resulting in many places where particles appear to originate. This is in-time \"pileup\". There may also be long-lived states that provide track origins displaced from the z-axis.\n",
    "\n",
    "Can you spot them? ;-)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# A Tracking Example: A Jet with a Long-Lived Subcomponent\n",
    "\n",
    "<img src=\"Images/bjet_fig_7.png\" style=\"width: 50%; height: auto; float: right;\"/>\n",
    "\n",
    "Finally, here is an example of an event with a jet candidate. The jet has additionally be determined to contain a long-lived component, consistent with the presence of a b-hadron. Only the Pixel Detector is shown.\n",
    "\n",
    "Note the displaced secondary \"vertex\" (blue). This is an old event from Run 1 - the IBL has not yet been installed!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# An Over-Idealized Particle Detector in an Over-Idealized Universe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjUAAAGeCAYAAABsJvAoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAGcFJREFUeJzt3X2spmddJ/Dvj+lmYFugWWRg8JlSGmoNdrEhZbUbKLjRNRCCKyDdGteGwUQCC7sZG1eR1V0xWcLKrCFisNsOxZW1FHBJILPlJQGxaC3aFgZ5KaDDzKkub4biTC3S4bd/nOfA4XBm5szpec5zzjWfT3Jnnuu6r/t+ftMrp+c792t1dwAAtruHzLsAAICNINQAAEMQagCAIQg1AMAQhBoAYAhCDQAwBKEGABjCTENNVT20qt5ZVXdX1V1VdUtVXThdt2va/kxVfbyqnnaK/bx4Ou5zVXVdVZ0zy7oBgO1nM47UXJfkku6+LMm7p+0keU2S27r74iQvSvKW1cJKVT0hyauTPC3JE5M8NsmLN6FuAGAbmWmo6e77u/tgf/uxxbcluWj6+YVJ3jAd95EkX8hicFnpBUn+T3d/YbqfNya5epZ1AwDbz2afxnlFkndV1aOSPKS7v7Rs3eEkF6yyzQVJPr+GcamqfUn2LbV37NjxvY997GMfZMkAwGa65557/rG7d57pdpsWaqrqlUkuTvKSJA9LsvKlU3WKzZePPem47t6fZP9SezKZ9MLCwpkXCwDMTVV96fSjvtum3P1UVdcmeV6SZ3X3fd39lWn/o5cNe3ySI6tsfiTJhWsYBwCcxWYeaqanhK5O8mPd/dVlq96W5GXTMU/N4gXAt66yi3ck+cmqekxVVRaP9Nw026oBgO1m1rd0T5K8Lsn5ST4wva37z6ar/1OSf1lVn0lyY5J/190PTLe7vqqemyTd/VdJfi3Jh5N8LskXk9wwy7oBgO2nvn1j0nhcUwMA209V3dPdkzPdzhOFAYAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADCEmYeaqnp9VR2uqq6qS6d951fVXcuWu6vqgar6Z6ts/8yqum/F+IfNum4AYHs5ZxO+4+1JXpvk1qWO7v5qksuW2lV1bZJndPffnWQfn+juy2daJQCwrc081HT3h5Kkqk417EVJfmXWtQAA45r7NTVVdUWSRyV59ymGXVJVd1TVR6rqpZtUGgCwjWzG6afT2Zvk97r7gZOsvyPJpLvvrapJkoNV9eXuvnnlwKral2TfUvuRj3zkTAoGALaeuR6pqapzk1yV5MDJxnT317r73unnhSR/kOTpJxm7v7snS8t55503i7IBgC1o3qeffirJx7r7UycbUFW7q+oh088PT/KcJHduUn0AwDaxGbd0v6GqFpJMkry/qj67bPWLk9ywyjbXV9Vzp83nJzlUVR9NcluS9yV504zLBgC2meruedcwM5PJpBcWFuZdBgBwBqrqnu6enOl28z79BACwIYQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxh5qGmql5fVYerqqvq0mX9h6vqU1V113S56hT7eFVVfW66vHrWNQMA2885m/Adb0/y2iS3rrLuBd398VNtXFVXJrk6yZOTPJDkw1V1a3e/Z8MrBQC2rZkfqenuD3X3woPYxVVJbuzu49399SQHshhyAAC+Zd7X1Lylqg5V1fVV9eiTjLkgyeeXtQ9P+75LVe2rqoWl5dixYxtcLgCwVc0z1FzZ3T+Y5ClJvpLkzacY28s+10kHde/v7snSct55521QqQDAVrcZ19SsqruPTP/8RlX9VpK7TzL0SJILl7UfP+0DAPiWuRypqapzq+r8ZV1XJ7nzJMPfluSa6TY7k+xNctOsawQAtpfNuKX7DVW1kGSS5P1V9dkkj0nygar6WFUdSvKMJD+7bJuDVXV5knT3B5PcnORQkk8meW933zLrugGA7aW6+/SjtqnJZNILCw/mxisAYLNV1T3dPTnT7eZ99xMAwIYQagCAIQg1AMAQhBoAYAhCDQAwBKEGABiCUAMADEGoAQCGINQAAEMQagCAIQg1AMAQhBoAYAhCDQAwBKEGABiCUAMADEGoAQCGINQAAEMQagCAIQg1AMAQhBoAYAhCDQAwBKEGABiCUAMADEGoAQCGINQAAEMQagCAIQg1AMAQhBoAYAhCDQAwBKEGABiCUAMADEGoAQCGINQAAEMQagCAIQg1AMAQhBoAYAgzDzVV9fqqOlxVXVWXTvseWlXvrKq7q+quqrqlqi48yfbPrKr7puOWlofNum4AYHvZjCM1b0/ytCSfX9F/XZJLuvuyJO+etk/mE9192bLlH2ZUKwCwTc081HT3h7p7YUXf/d19sLt72nVbkotmXQsAMK6tck3NK5K86xTrL6mqO6rqI1X10pMNqqp9VbWwtBw7dmzjKwUAtqRz5l1AVb0yycVJXnKSIXckmXT3vVU1SXKwqr7c3TevHNjd+5PsX2pPJpNeOQYAGNNcQ01VXZvkeUl+tLvvW21Md39t2eeFqvqDJE9P8l2hZqUTJ07k6NGjG1UuAGeBXbt2ZefOnfMug3WYW6ipqn1Jrs5ioPnqKcbtTvKF7v5mVT08yXOS3LCW7zh+/HgOHDiwIfUCcHbYu3dv9uzZM+8yWIfNuKX7DVW1kGSS5P1V9dnpaaTXJTk/yQemt2n/2bJtrq+q506bz09yqKo+msULit+X5E2zrhsA2F5mfqSmu1+W5GWrrKpTbPNzyz7/dpLfnkFpAMBAtsrdTwAAD0p9+1Ex49m9e3fffvvt8y4DgG3EhcLzV1X3dPfkTLeb+y3ds7Rjxw4XewHAWcLpJwBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBDOOdXKqrp5Dfv4u+5+yQbVAwCwLqcMNUl+KMmvnmbML21QLQAA63a6UPPW7n7zqQZU1Q9sYD0AAOtS3T3vGmZmMpn0wsLCvMsAAM5AVd3T3ZMz3e50R2qWf8ETk1y0fJvuPnimXwgAMAtrCjVV9bokP5Pk00lOTLs7iVADAGwJaz1S8xNJntDd982yGACA9Vrrc2o+n+QfZ1kIAMCDsdYjNb+Q5F1V9b4k9y91dvfvzKQqAIAztNYjNb+cZHeSy5I8dbpcvpYNq+r1VXW4qrqqLl3Wf3FV/UlV3V1Vt1fVk06xj1dV1eemy6vXWDMAcBZZ65GapyT5vl7f/d9vT/LaJLeu6P/dJNd1941V9YIkNyS5YuXGVXVlkquTPDnJA0k+XFW3dvd71lELADCotR6p+VSSc9fzBd39oe7+jofFVNWuLAal3592vSPJE6rqwlV2cVWSG7v7eHd/PcmBLIYcAIBvWeuRmr9P8hdV9Z585zU1v7jO792T5G+6+4HpfrqqjiS5IMnhFWMvSPJHy9qHk7xgnd8LAAxqraHm09NlI608lVVrHHvScVW1L8m+pfYjH/nI9VUGAGw7awo13f1fN/h7jyaZVNU53f1AVVUWj94cWWXskSQXLms//iTj0t37k+xfak8mk3HfAQEAfIdTXlNTVXtPt4O1jFmpu7+Y5M4sPqU4SZ6f5HB3H15l+NuSXFNV51bVziR7k9x0pt8JAIztdEdqrq2qP82pTw39xyxevLuqqnpDFp9I/Ngk76+qY939xCQ/n+TGqnplkq8luWbZNgeT/Gp3/3l3f7Cqbk5yaLr6pu6+5XR/MQDg7HLKt3RX1eF897UvK32pu//FRha1UbylGwC2n5m8pbu7L1x3RQAAm2itz6kBANjShBoAYAhCDQAwhDWFmqo6f9aFAAA8GGs9UvOZqvqfVfWDM60GAGCd1hpqnpjkk0neXlV/XFUvrKodM6wLAOCMrCnUdPe93b2/uy9O8pokv5nkSFX9SlWt6+3dAAAbac0XClfVw6vq5VkMNH+Z5OVJHpfE030BgLlb0wstq+qNWXzVwTuS/JvuXnpj9x9W1SdnVRwAwFqtKdQk+WyS7+/ue1dZ9682sB4AgHVZU6jp7t88xbq/3bhyAADWZ61HaralEydO5OjRo/MuA4BtZNeuXdm5c+e8y2Adhg41x48fz4EDB+ZdBgDbyN69e7Nnz555l8E6eE0CADAEoQYAGIJQAwAMobp73jXMzO7du/v222+fdxkAbCMuFJ6/qrqnuydnut3QFwrv2LHDxV4AcJZw+gkAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAM4Zx5fXFVnZ/kg8u6/mmSi5Ls6u6/WzbumUkOJrl72dgruvsfNqFMAGCbmFuo6e6vJrlsqV1V1yZ5xvJAs8wnuvvyTSsOANh2ttLppxcluWHeRQAA29OWCDVVdUWSRyV590mGXFJVd1TVR6rqpafYz76qWlhajh07NpN6AYCtZ26nn1bYm+T3uvuBVdbdkWTS3fdW1STJwar6cnffvHJgd+9Psn+pPZlMemYVAwBbytyP1FTVuUmuSnJgtfXd/bXuvnf6eSHJHyR5+uZVCABsB3MPNUl+KsnHuvtTq62sqt1V9ZDp54cneU6SOzexPgBgG9gKoebFWXGBcFVdX1XPnTafn+RQVX00yW1J3pfkTZtbIgCw1VX3uJedTCaTXlhYmHcZAMAZqKp7untypttthSM1AAAPmlADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhnDOvAuYpRMnTuTo0aPzLgOAbWTXrl3ZuXPnvMtgHeYaaqrqcJL7p0uS/Lfufusq416V5EXT5v/u7v+8lv0fP348Bw4c2IhSAThL7N27N3v27Jl3GazDVjhS84Lu/vjJVlbVlUmuTvLkJA8k+XBV3drd79msAgGArW87XFNzVZIbu/t4d389yYEshhwAgG/ZCqHmLVV1qKqur6pHr7L+giSfX9Y+PO37LlW1r6oWlpYZ1AoAbFHzPv10ZXcfqap/kuQ3krw5ybNXGdfLPtfJdtbd+5PsX2rv3r279+7du1G1AnAW2LVr17xLYJ3mGmq6+8j0z29U1W8luXuVYUeSXLis/fhp32nt2LHDxV4AcJaY2+mnqjq3qs5f1nV1kjtXGfq2JNdMx+9MsjfJTZtRIwCwfczzmprHJPlAVX2sqg4leUaSn02SqjpYVZcnSXd/MMnNSQ4l+WSS93b3LfMpGQDYqqq7Tz9qm5pMJr2w4HphANhOquqe7p6c6XZb4e4nAIAHTagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGgBgCEINADAEoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhnDOvL66qhya5KcmTktyX5P8leUl3H14x7plJDia5e1n3Fd39D6f7jhMnTuTo0aMbVTIAZ4Fdu3Zl586d8y6DdZhbqJm6Lsn/7e6uqn8/bf/rVcZ9orsvP9OdHz9+PAcOHHiwNQJwFtm7d2/27Nkz7zJYh7mdfuru+7v7YHf3tOu2JBfNqx4AYHvbStfUvCLJu06y7pKquqOqPlJVLz3ZDqpqX1UtLC2zKRMA2IrmffopSVJVr0xycZKXrLL6jiST7r63qiZJDlbVl7v75pUDu3t/kv1L7Uc84hG9cgwAMKb69tmfORVQdW2Sf5vkR7v7q2sY/8tJHtfdLz/d2N27d/ftt9++AVUCcLZwofD8VdU93T050+3meqSmqvYluTqnCDRVtTvJF7r7m1X18CTPSXLDWva/Y8cOF3sBwFlibtfUTE8lvS7J+Uk+UFV3VdWfTdddX1XPnQ59fpJDVfXRLF5M/L4kb5pHzQDA1jX300+zNJlMemHB9cIAsJ2s9/TTVrr7CQBg3YQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYglADAAxBqAEAhiDUAABDEGoAgCEINQDAEIQaAGAIQg0AMAShBgAYwlxDTVVdXFV/UlV3V9XtVfWkk4x7VVV9brq8erPrBAC2vnkfqfndJNd19/cleW2SG1YOqKork1yd5MlJnpTkWVX145taJQCw5c0t1FTVriRPSfL70653JHlCVV24YuhVSW7s7uPd/fUkB7IYcgAAvmWeR2r2JPmb7n4gSbq7kxxJcsGKcRck+fyy9uFVxgAAZ7lz5vz9vaJdaxh3sjGpqn1J9i3r+mZV/e06a2NjnZfk2LyLIIm52ErMxdZiPraOx65no3mGmqNJJlV1Tnc/UFWVxaM3R1aMO5LkwmXtx68yJknS3fuT7F9qV9VCd082tGrWxVxsHeZi6zAXW4v52DqqamE9283t9FN3fzHJnUl+Ztr1/CSHu/vwiqFvS3JNVZ1bVTuT7E1y06YVCgBsC/O+++nnk/x8Vd2d5JeSvDhJqupgVV2eJN39wSQ3JzmU5JNJ3tvdt8ynXABgq5rrNTXd/ekkV6zS/+wV7V9P8uvr+Ir9px/CJjEXW4e52DrMxdZiPraOdc1FLd50BACwvc379BMAwIYQagCAIQwZatb6Tik2XlW9vqoOV1VX1aXL+s3JJquqh1bVO6f/ze+qqluWnthdVbum7c9U1cer6mnzrXZ8VfXeqvrYdC7+uKoum/b72ZiTqvq15f+vMhfzMf2d8anpz8ZdVXXVtP+M52PIUJM1vFOKmXl7kqflO58CnZiTebkuySXdfVmSd0/bSfKaJLd198VJXpTkLVU174dxju6F3f3k6Vy8LouvfEn8bMxFVT0lyQ/nO597Zi7m5wXdfdl0eeu074znY7gLhafvlLo7yfcse6jf3yb54VWegcOMVNXhJM/p7o+bk61h+piEm7r7iVV1LMkTuvtL03W3J/nF6SMUmLGquibJy5M8O342Nt30mWcfTPLTST6Q5DlJvhhzMRfLf18s61vX740Rj9Ss9Z1SbB5zsjW8Ism7qupRSR6yFGimDsd8zFxV/V5VHU3yG0muiZ+Nefn1JL/f3X+9rM9czNdbqupQVV1fVY/OOudjxFCTrP2dUmweczJHVfXKJBcn+ZVpl/mYg+7+2e7ek+RVSf77UveKYeZihqrqiiRPTfI7q6w2F/NxZXf/YJKnJPlKkjdP+894PkYMNd96p1SSnOKdUmweczJHVXVtkucleVZ339fdX5n2P3rZsJO+U42N191vTvIjSRbiZ2OzPSPJ9yf56+lpj0mS9yS5NOZiLrr7yPTPbyT5rSRPzzp/bwwXas7gnVJsEnMyP9M311+d5Me6+6vLVr0tycumY56axTfi3rr5FZ4dquoRVfW4Ze2fzOK/SP1sbLLufk13P667L+zuC7MYLH98GjTNxSabvtfx/GVdVye5c72/N4a7UDhJquqSJDcmeVSSryW5prv/cq5FnSWq6g1JfiKLvyS/nOTY9MJUc7LJqmqSxX/t/FWSv592f727f6iqHpPkfyV5QpJ/TPLS7v6j+VQ6vqrak+QdSR6W5JtJvpTk2u6+y8/GfK24qcFcbLKquiiLPxs7snh66a+S/IfuPrye+Rgy1AAAZ5/hTj8BAGcnoQYAGIJQAwAMQagBAIYg1AAAQxBqAIAhCDUAwBCEGmDLq6obq2qhqt64jm2vqqpPVNWXZ1EbsHV4+B6w5VXVjUn+vLt/e53bXzjd/ns2sCxgi3GkBpibqnp2Vd21bLm/qq5Zw3bPnI5/Y1Udqqo7qurSqnrr9KjM+6rqvM34OwBbh1ADzE13H+zuy7r7siS/meSTSf5wjZv/QJI3dvc/T/KnSW5J8gvd/aQk30jy07OoGdi6zpl3AQBV9SNJ/kuSK7v7708zfMmnu/uu6ec7kjy+uxem7b9IctHGVglsdY7UAHNVVZcmeVOSn+juvzmDTe9f9vnEKm3/aIOzjFADzE1VfW+SdyZ5UXf/5bzrAbY3oQaYp59L8ugk/2PZxcLPnXdRwPbklm5gy3NLN7AWjtQA28G9SV6x3ofvJXlXki9seFXAluJIDQAwBEdqAIAhCDUAwBCEGgBgCEINADAEoQYAGML/B+PVx0RmNtxRAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Drawing plane\n",
    "fig, ax = plt.subplots(figsize=(8, 6), dpi=80)\n",
    "\n",
    "ax.set_xlim(( 0, 50))\n",
    "ax.set_ylim((0, 20))\n",
    "ax.set_xlabel(\"z [m]\")\n",
    "ax.set_ylabel(\"y [m]\")\n",
    "\n",
    "# \"Detector Planes\"\n",
    "plt.axhline(y=2.5, xmin=0, xmax=0.8, linewidth=4, color = 'gray')\n",
    "plt.axhline(y=5.05, xmin=0, xmax=0.8, linewidth=4, color = 'gray')\n",
    "plt.axhline(y=8.85, xmin=0, xmax=0.8, linewidth=4, color = 'gray')\n",
    "plt.axhline(y=12.25, xmin=0, xmax=0.8, linewidth=4, color = 'gray')\n",
    "\n",
    "line, = ax.plot([], [], lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "To think basically about tracking, let's boil down the problem. The above illustrates a side-view (y-z plane) of a \"ATLAS Pixel-like\" detector, consisting of 4 planes (gray lines) that are part of a cylindrical barrel enclosing the z-axis. Each plane is made from sensitive units (e.g. \"pixels\") that respond to energy deposited by a passing particle. I will ignore the use of a magnetic field to bend particles (for now) merely to isolate the basic tracking questions. All action happens in the plane shown.\n",
    "\n",
    "Assume each pixel in each plane is 100% efficient with 0% noise."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# An Over-Idealized Particle Detector in an Over-Idealized Universe"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": [
    "# interactive particle simulation\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import Button, Layout\n",
    "from ipywidgets import interact, interactive, fixed, interact_manual\n",
    "import numpy as np\n",
    "import scipy as sp\n",
    "\n",
    "class DetectorPlane:\n",
    "    def __init__(self,x=-1,y=-1, width=-1):\n",
    "        self.y = y\n",
    "        self.x = x\n",
    "        self.width = width\n",
    "\n",
    "class Pixel:\n",
    "    def __init__(self, x=-1, y=-1, xwidth=-1):\n",
    "        self.x = x\n",
    "        self.y = y\n",
    "        self.xwidth = xwidth\n",
    "    \n",
    "    \n",
    "# Draw the particle passing through the planes\n",
    "def particle_line_y(x, slope, intercept):\n",
    "    return slope*x + intercept\n",
    "\n",
    "def particle_line_x(y, slope, intercept):\n",
    "    if slope == 0.0:\n",
    "        slope = 1e-6\n",
    "    return (y - intercept)/slope\n",
    "\n",
    "def GenerateParticle(angle=math.pi/4, intercept=0.0):\n",
    "    slope = math.tan(angle*math.pi/180)\n",
    "    particle_x = [-1, 50]\n",
    "    particle_y = [particle_line_y(-1,slope,intercept),particle_line_y(50,slope,intercept)]\n",
    "    return [particle_x, particle_y]\n",
    "\n",
    "\n",
    "def HitChiSquareFitter(detector_planes=[], pixelhits=[]):\n",
    "    # do a basic chi-square minimization and return the \"best fit\" particle trajectory hypothesis\n",
    "    angle_hypo = np.arange(start=5, stop=90, step=1)\n",
    "    global_chi_square = 1e9\n",
    "    best_angle = 1e9\n",
    "    \n",
    "    for angle in angle_hypo:\n",
    "        #print(f\"Angle Hypothesis: {angle} degrees\")\n",
    "        # get the x coordinates in each plane\n",
    "        hits_hypo = []\n",
    "        for plane in detector_planes:\n",
    "            x_hypo = particle_line_x(plane.y, math.tan(angle*math.pi/180), 0.0)\n",
    "            y_hypo = plane.y\n",
    "            hits_hypo.append([x_hypo/plane.width, y_hypo])\n",
    "        #print(f\"  {hits_hypo}\")\n",
    "        # Compute the chi-square of this hypothesis\n",
    "        chi_square = 0\n",
    "        #print(pixelhits)\n",
    "        for hit_hypo in hits_hypo:\n",
    "            for pixelhit in pixelhits:\n",
    "                #print(f\"      {hit_hypo[1]} vs. {pixelhit.y}\")\n",
    "                if hit_hypo[1] == pixelhit.y:\n",
    "                    this_chi_square = (hit_hypo[0] - pixelhit.x)**2/pixelhit.xwidth**2\n",
    "                    #print(f\"      {this_chi_square}\")\n",
    "                    chi_square += (hit_hypo[0] - pixelhit.x)**2/pixelhit.xwidth**2\n",
    "        #print(f\"   X^2 = {chi_square}\")\n",
    "        if chi_square < global_chi_square:\n",
    "            global_chi_square = chi_square\n",
    "            best_angle = angle\n",
    "    \n",
    "    #print(best_angle)\n",
    "    \n",
    "    return best_angle\n",
    "    \n",
    "    \n",
    "\n",
    "def DrawPerfectTracker(planes):\n",
    "    # \"Detector Planes\"\n",
    "    for y in planes:\n",
    "        plt.axhline(y=y, xmin=0, xmax=0.8, linewidth=4, color = 'gray')\n",
    "    return\n",
    "\n",
    "def DrawGranularTracker(planes, granularity=10):\n",
    "    UnitsPerPlane = granularity\n",
    "    PhysicalPlaneLength = 50 # cm\n",
    "    RelativePlaneLength = 0.8 # relative to x-axis shown, so 0.8 = 80%\n",
    "    PixelWidth = RelativePlaneLength/(2*UnitsPerPlane)\n",
    "    xcenters = np.arange(start=PixelWidth, stop=2*PixelWidth*UnitsPerPlane, step=PixelWidth*2)\n",
    "    PixelList = [] # [y, center, width] for each pixel\n",
    "    for y in planes:\n",
    "        for x in xcenters:\n",
    "            this_hit = Pixel(x=x, y=y, xwidth=PixelWidth)\n",
    "            PixelList.append(this_hit)\n",
    "            #PixelList.append([y, x*PhysicalPlaneLength, PixelWidth*PhysicalPlaneLength])\n",
    "            plt.axhline(y=y, xmin=x-PixelWidth, xmax=x+PixelWidth, linewidth=4, color = 'gray')\n",
    "        \n",
    "    return PixelList\n",
    "\n",
    "def DrawPerfectEvent(angle=45,showhits=False):\n",
    "    fig, ax = plt.subplots(figsize=(8, 6), dpi=80)\n",
    "\n",
    "    ax.set_xlim(( 0, 50))\n",
    "    ax.set_ylim((0, 20))\n",
    "    ax.set_xlabel(\"z [m]\")\n",
    "    ax.set_ylabel(\"y [m]\")\n",
    "\n",
    "    plane_y_coordinates = [2.5, 5.05, 8.85, 12.25]\n",
    "\n",
    "    DrawPerfectTracker(plane_y_coordinates)\n",
    "        \n",
    "    intercept = 0.0\n",
    "    the_particle = GenerateParticle(angle=angle, intercept=intercept)\n",
    "    \n",
    "    plt.plot(the_particle[0], the_particle[1] ,marker = 'o')\n",
    "    \n",
    "    \n",
    "    # If we are to show hits, draw them too\n",
    "    if showhits == True:\n",
    "        for y in plane_y_coordinates:\n",
    "            slope = math.tan(angle*math.pi/180)\n",
    "            x_hit = particle_line_x(y, slope, intercept)\n",
    "            if x_hit > 40:\n",
    "                continue\n",
    "            #print([particle_line_x(y, slope, intercept), y])\n",
    "            plt.plot(x_hit, y, marker='s', markeredgecolor='k', markersize=10, color='r')\n",
    "\n",
    "           \n",
    "    line, = ax.plot([], [], lw=2)\n",
    "\n",
    "def DrawLossyEvent(angle=45, hitefficiency=1.0):\n",
    "    fig, ax = plt.subplots(figsize=(8, 6), dpi=80)\n",
    "\n",
    "    ax.set_xlim(( 0, 50))\n",
    "    ax.set_ylim((0, 20))\n",
    "    ax.set_xlabel(\"z [m]\")\n",
    "    ax.set_ylabel(\"y [m]\")\n",
    "\n",
    "    # \"Detector Planes\"\n",
    "    plane_y_coordinates = [2.5, 5.05, 8.85, 12.25]\n",
    "    DrawPerfectTracker(plane_y_coordinates)\n",
    "        \n",
    "    intercept = 0.0\n",
    "    the_particle = GenerateParticle(angle=angle, intercept=intercept)\n",
    "    \n",
    "    plt.plot(the_particle[0], the_particle[1] ,marker = 'o')\n",
    "    \n",
    "    # Draw hits, folding in hit efficiency\n",
    "    for y in plane_y_coordinates:\n",
    "        slope = math.tan(angle*math.pi/180)\n",
    "        x_hit = particle_line_x(y, slope, intercept)\n",
    "        if x_hit > 40:\n",
    "            continue\n",
    "\n",
    "        showhit = (np.random.uniform(0,1) < hitefficiency)\n",
    "        if showhit:\n",
    "            plt.plot(x_hit, y, marker='s', markeredgecolor='k', markersize=10, color='r')\n",
    "\n",
    "           \n",
    "    line, = ax.plot([], [], lw=2)\n",
    "\n",
    "\n",
    "def DrawGranularEvent(angle=45,showhits=False):\n",
    "    fig, ax = plt.subplots(figsize=(8, 6), dpi=80)\n",
    "\n",
    "    ax.set_xlim(( 0, 50))\n",
    "    ax.set_ylim((0, 20))\n",
    "    ax.set_xlabel(\"z [m]\")\n",
    "    ax.set_ylabel(\"y [m]\")\n",
    "\n",
    "    plane_y_coordinates = [2.5, 5.05, 8.85, 12.25]\n",
    "\n",
    "    pixel_list = DrawGranularTracker(plane_y_coordinates)\n",
    "        \n",
    "    intercept = 0.0\n",
    "    the_particle = GenerateParticle(angle=angle, intercept=intercept)\n",
    "    \n",
    "    plt.plot(the_particle[0], the_particle[1] ,marker = 'o')\n",
    "    \n",
    "    \n",
    "    # If we are to show hits, draw them too\n",
    "    HitsList = []\n",
    "    if showhits == True:\n",
    "        for y in plane_y_coordinates:\n",
    "            slope = math.tan(angle*math.pi/180)\n",
    "            x_hit = particle_line_x(y, slope, intercept)\n",
    "            for pixel in pixel_list:\n",
    "                if y == pixel.y and ((pixel.x-pixel.xwidth) < x_hit/50 and x_hit/50 < (pixel.x+pixel.xwidth)): \n",
    "                    HitsList.append(pixel)\n",
    "                    plt.axhline(y=pixel.y, \n",
    "                                xmin=(pixel.x-pixel.xwidth), \n",
    "                                xmax=(pixel.x+pixel.xwidth), \n",
    "                                linewidth=4, color = 'r')\n",
    "\n",
    "           \n",
    "    line, = ax.plot([], [], lw=2)\n",
    "\n",
    "    \n",
    "def DrawGranularLossyEvent(angle=45, \n",
    "                           granularity=10, \n",
    "                           showhits=False, \n",
    "                           hitefficiency=1.0, \n",
    "                           showtruth=True, \n",
    "                           showfit=False):\n",
    "    \n",
    "    fig, ax = plt.subplots(figsize=(8, 6), dpi=80)\n",
    "\n",
    "    ax.set_xlim(( 0, 50))\n",
    "    ax.set_ylim((0, 20))\n",
    "    ax.set_xlabel(\"z [m]\")\n",
    "    ax.set_ylabel(\"y [m]\")\n",
    "\n",
    "    plane_y_coordinates = [2.5, 5.05, 8.85, 12.25]\n",
    "    DetectorPlanes = []\n",
    "    for plane_y in plane_y_coordinates:\n",
    "        this_plane = DetectorPlane(x=-1, y=plane_y, width=50)\n",
    "        DetectorPlanes.append(this_plane)\n",
    "        \n",
    "\n",
    "    pixel_list = DrawGranularTracker(plane_y_coordinates, granularity)\n",
    "        \n",
    "    intercept = 0.0\n",
    "    the_particle = GenerateParticle(angle=angle, intercept=intercept)\n",
    "    \n",
    "    if showtruth == True:\n",
    "        plt.plot(the_particle[0], the_particle[1] ,marker = 'o')\n",
    "    \n",
    "    \n",
    "    # If we are to show hits, draw them too\n",
    "    HitsList = []\n",
    "    \n",
    "    for y in plane_y_coordinates:\n",
    "        slope = math.tan(angle*math.pi/180)\n",
    "        x_hit = particle_line_x(y, slope, intercept)\n",
    "        x_hit = particle_line_x(y, slope, intercept)\n",
    "        for pixel in pixel_list:\n",
    "            if y == pixel.y and ((pixel.x-pixel.xwidth) < x_hit/50 and x_hit/50 < (pixel.x+pixel.xwidth)): \n",
    "                detected_hit = (np.random.uniform(0,1) < hitefficiency)\n",
    "                if detected_hit:\n",
    "                    HitsList.append(pixel)\n",
    "\n",
    "                if showhits and detected_hit:\n",
    "                    plt.axhline(y=pixel.y, \n",
    "                                xmin=(pixel.x-pixel.xwidth), \n",
    "                                xmax=(pixel.x+pixel.xwidth), \n",
    "                                linewidth=4, color = 'r')\n",
    "    \n",
    "    # Best fit to the hits\n",
    "    angle_hypo = HitChiSquareFitter(detector_planes=DetectorPlanes, pixelhits=HitsList)\n",
    "    if showfit:\n",
    "        fitted_track = GenerateParticle(angle=angle_hypo, intercept=intercept)\n",
    "        plt.plot(fitted_track[0], fitted_track[1], color='g' ,marker = 'o')\n",
    "  \n",
    "    \n",
    "    line, = ax.plot([], [], lw=2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "code_folding": [
     0
    ],
    "scrolled": false,
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a6a23f430c6e4bfabbc09c78059857d7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatSlider(value=0.0, continuous_update=False, max=90.0, step=2.0), ToggleButton(value=False, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0befab8d226a43b6a8d7cf5025e8d7b7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This example excludes showing \"hits\" in the detector\n",
    "angle_slider = widgets.FloatSlider(min=0.0, max=90, step=2, continuous_update=False)\n",
    "intercept_slider = widgets.FloatSlider(min=0.0, max=2.5, step=0.01, continuous_update=False)\n",
    "showhits_button = widgets.ToggleButton(value=False,description='Show Hits?')\n",
    "\n",
    "input_widgets = widgets.HBox([angle_slider, showhits_button])\n",
    "#interactive_plot = interactive(DrawPerfectEvent, angle=angle_slider, showhits=)\n",
    "#output = interactive_plot.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "#interactive_plot\n",
    "interactive_output = widgets.interactive_output(DrawPerfectEvent, \n",
    "                                                {'angle': angle_slider, 'showhits': showhits_button})\n",
    "#output = interactive_output.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "\n",
    "display(input_widgets, interactive_output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The above program will let you generate a \"particle\" (blue line) that passes through the layers of the detector. The particle alway arises from the origin of the coordinate system. \n",
    "\n",
    "Toggle the \"hits\" option to see \"strikes\" on the detector planes at 100% efficiency. This detector is modeled as if infinitely pixellated - real detectors have finite sizes to the pixels, but this is an idealization (we'll complicate this more in a moment)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# A Less-Ideal Tracker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d36b526b90a64ccb9d025ca2b594f1ef",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatSlider(value=45.0, continuous_update=False, max=90.0, step=2.0), FloatSlider(value=1.0, co…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d5d1f4b9fd7e4f10bb92d754056ab939",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This example excludes showing \"hits\" in the detector, with loss of efficiency\n",
    "angle_slider = widgets.FloatSlider(value=45, min=0.0, max=90, step=2, continuous_update=False)\n",
    "hitefficiency_slider = widgets.FloatSlider(value=1.0, min=0.0, max=1.0, step=0.01, continuous_update=False)\n",
    "\n",
    "input_widgets = widgets.HBox([angle_slider, hitefficiency_slider])\n",
    "interactive_output = widgets.interactive_output(DrawLossyEvent, \n",
    "                                                {'angle': angle_slider, 'hitefficiency': hitefficiency_slider},\n",
    "                                               )\n",
    "#output = interactive_output.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "\n",
    "display(input_widgets, interactive_output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "In this example, the tracker becomes \"lossy\" - even if a particle leaves energy in a pixel, that pixel is not guaranteed to yield exactly 1 \"hit\" as a result. The *efficiency per hit* can be controlled with the above. You can already see at 90-95% per-hit efficiency, you can start noticeably losing hits on tracks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# A Granular Tracker"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "code_folding": [
     0
    ],
    "slideshow": {
     "slide_type": "-"
    }
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "179a8907b6ae4f8bb2db68d2d2d69e1d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatSlider(value=0.0, continuous_update=False, max=90.0, step=2.0), ToggleButton(value=False, …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ebfc720dcbae43f18faa10bd5feb0a05",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This example is a granular detector, but perfect in efficiency\n",
    "angle_slider = widgets.FloatSlider(min=0.0, max=90, step=2, continuous_update=False)\n",
    "intercept_slider = widgets.FloatSlider(min=0.0, max=2.5, step=0.01, continuous_update=False)\n",
    "showhits_button = widgets.ToggleButton(value=False,description='Show Hits?')\n",
    "\n",
    "input_widgets = widgets.HBox([angle_slider, showhits_button])\n",
    "#interactive_plot = interactive(DrawPerfectEvent, angle=angle_slider, showhits=)\n",
    "#output = interactive_plot.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "#interactive_plot\n",
    "interactive_output = widgets.interactive_output(DrawGranularEvent, \n",
    "                                                {'angle': angle_slider, 'showhits': showhits_button})\n",
    "#output = interactive_output.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "\n",
    "display(input_widgets, interactive_output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "The above illustrates a slightly more realistic tracker (albeit in this case with 100% per-hit efficiency). Now the tracker planes are build from 10 sub-units. A \"hit\" is now defined by the x-coordinate of the center of a pixel unit and the y-coordinate of its vertical plane. The width of the pixel unit defines the *resolution* on the hit... now we suddenly have a sense of uncertainty in hit tracking and from where it arises."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# A Granular, Lossy Tracker - Some Basic Theory\n",
    "\n",
    "For the idealized world I have created - a single ionizing particle always emitted from the origin of the coordinate system - we can now consider a *granular, lossy* tracking detector and how one uses the information about which pixels are hit to form an hypothesis about the trajectory of the particle that left the hits.\n",
    "\n",
    "* Starting all hypothetical trajectories at the origin and assuming a straight line, scan over angle (slope) hypotheses, e.g. $\\theta = [5,89]^{\\circ}$ in 1-degree steps;\n",
    "* For each angular (slope) hypothesis, compute the x-coordinate of where that trajector would have intersected each y-plane of the detector system, building up hypotheses for hits: $x_{hypo,i}$ where $i$ runs over layer planes.\n",
    "* Using all hypothetical hits in layer $i$, compute the $\\chi^2_{i}$ with the actual hit in that layer, $$\\chi^2_i = \\left( \\frac{x_{hypo,i} - x_{i}}{\\sigma_{i}} \\right)^2$$ where $x_{i}$ is the actual center of the pixel hit in layer $i$ and $\\sigma_{i}$ is the width of the pixel that was hit.\n",
    "* Compute the total $\\chi^2 = \\sum_{i} \\chi^2_{i}$ - **the \"best fit\" will be defined as the hypothesis that achieves the lowest $\\chi^2$**\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# A Granular, Lossy Tracker - Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "76ec7d1e3724440cbe46ca9f20968957",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(HBox(children=(FloatSlider(value=45.0, continuous_update=False, max=90.0, step=2.0), FloatSlide…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1a8e49cd6626492dbd3c682f8dcd65e7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This example is a granular detector, with per-hit efficiency losses\n",
    "angle_slider = widgets.FloatSlider(value=45, min=0.0, max=90, step=2, continuous_update=False)\n",
    "granularity_slider = widgets.FloatSlider(value=10, min=10, max=100, step=5, continuous_update=False)\n",
    "hitefficiency_slider = widgets.FloatSlider(value=1.0, min=0.0, max=1.0, step=0.01, continuous_update=False)\n",
    "showhits_button = widgets.ToggleButton(value=False,description='Show Hits?')\n",
    "showtruth_button = widgets.ToggleButton(value=True,description='Show Truth?')\n",
    "showfit_button = widgets.ToggleButton(value=False,description='Show Best Fit?')\n",
    "\n",
    "slider_row = widgets.HBox([angle_slider, hitefficiency_slider, granularity_slider])\n",
    "button_row = widgets.HBox([showhits_button, showtruth_button, showfit_button])\n",
    "\n",
    "input_widgets = widgets.VBox([slider_row, button_row])\n",
    "#interactive_plot = interactive(DrawPerfectEvent, angle=angle_slider, showhits=)\n",
    "#output = interactive_plot.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "#interactive_plot\n",
    "interactive_output = widgets.interactive_output(DrawGranularLossyEvent, \n",
    "                                                {'angle': angle_slider, \n",
    "                                                 'granularity': granularity_slider,\n",
    "                                                 'hitefficiency': hitefficiency_slider, \n",
    "                                                 'showhits': showhits_button,\n",
    "                                                'showtruth': showtruth_button,\n",
    "                                                'showfit': showfit_button})\n",
    "#output = interactive_output.children[-1]\n",
    "#output.layout.height = '450px'\n",
    "\n",
    "display(input_widgets, interactive_output)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "-"
    }
   },
   "source": [
    "This example allows control of granularity, per-hit efficiency, and true particle angle. Play around and see the effect of granularity and per-hit efficiency on the ability to resolve the \"truth\" about the particle's trajectory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Closing Thoughts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "# Reality - the next steps\n",
    "\n",
    "* In reality\n",
    "    * ...each proton-proton collision produces something like dozens, hundreds or thousands of particles somewhere along the $z$ axis\n",
    "    * ... due to immersion in a solenoidal magnetic field, trajectories are not straight - they are helical, having 5 parameters to describe them, not 2.\n",
    "    * ... due to the fact that long-lived particles decay *away* from the interaction point, you can't assume that a trajectory begins at the interaction point\n",
    "    * ... detectors have noise and dead material; particle energy losses are not perfect and uniform; etc.\n",
    "* This talk highlights the basic ideas needed to even begin to think about the problem:\n",
    "    * Expect jets to create prompt charged particles and displaced charged particles;\n",
    "    * Tracking detectors need layering and granularity in the layering to attain sufficient resolution to see these features\n",
    "    * Tracking algorithms need to be complex and adapt to the varying conditions of the detector and under which the charged particles are produced in the first place"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Learn More\n",
    "\n",
    "* Andy Salzburg. [\"Tracking Challenge\"](https://erez.weizmann.ac.il/pls/htmldb/f?p=101:58:::NO:RP:P58_CODE,P58_FILE:5446,Y). Hammers and Nails Workshop. 2017. Weizmann Institute of Science."
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "rise": {
   "footer": "<h3>Stephen Sekula - SMU - June 19, 2019</h3>"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}