muon_mass.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "VxPAdfoTjv0X"
   },
   "source": [
    "# Calculating Invariant Mass \n",
    "\n",
    "This activity uses data from the [CMS detector](https://cms.cern/detector) at CERN in Geneva, Switzerland. We've used this in [Quarknet's Data Camp at Fermilab](https://quarknet.org/page/data-camp-description) for several years to help teachers learn about particle physics.  \n",
    "\n",
    "To get started,\n",
    "- You won't hurt anything by experimenting. If you break it, close the tab and open the activity again to start over.\n",
    "- Is this your first time? Need a refresher? Try the 5-minute [Intro to Coding activity](./intro.ipynb) and come back here. \n",
    "\n",
    "When you're ready, run each code cell until you get down to **Part One**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "NiDzd_bwjv0Z"
   },
   "outputs": [],
   "source": [
    "# imports some software packages we'll use\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 143
    },
    "id": "WMTlKXuwjv0j",
    "outputId": "51c686ba-9bdb-449a-ce0d-ce7c28925c68"
   },
   "outputs": [],
   "source": [
    "# a hashtag tells the program \"don't read the rest of the line\"\n",
    "# That way we can write \"comments\" to humans trying to figure out what the code does\n",
    "\n",
    "data = pd.read_csv('https://github.com/QuarkNet-HEP/coding-camp/raw/main/data/muons.csv')\n",
    "# units in these files are energy, E (in GeV) and momentum, p (in GeV/c)\n",
    "\n",
    "# The .head(n) command displays the first n rows of a file.\n",
    "data.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "Bdlj4S83jv0z",
    "outputId": "1141b68c-59c6-4b93-fc7d-d73225173cef"
   },
   "outputs": [],
   "source": [
    "# .shape displays the (number of rows , number of columns)\n",
    "data.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "OezzANGxjv04"
   },
   "source": [
    "## Part One\n",
    "Get acquainted with this data set. It represents muons produced in a proton-proton collision (called an *event*). Look at the cells above to find the answers to the following questions:\n",
    "- In the table above, what do you think each of the column headings represent? After you discuss with your partner, take a look at the [key](https://github.com/QuarkNet-HEP/coding-camp/blob/main/CMS_data_headings.md).\n",
    "- How many events does this data set contain?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "KufRU5EUjv0_"
   },
   "source": [
    "## Part Two\n",
    "Each muon has momentum and since they travel in 3 dimensions the momentum has three components. You can calculate the total momentum with the Pythogorean Theorem:  \n",
    "- In 2 dimensions, it's the familiar:  \n",
    "|p|<sup>2</sup> = p<sub>x</sub><sup>2</sup> + p<sub>y</sub><sup>2</sup>  \n",
    "  \n",
    "  \n",
    "- In 3 dimensions, it looks like:  \n",
    "|p|<sup>2</sup> = p<sub>x</sub><sup>2</sup> + p<sub>y</sub><sup>2</sup> + p<sub>z</sub><sup>2</sup>  \n",
    "  \n",
    "  \n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/6/64/Coord_XYZ.svg/2000px-Coord_XYZ.svg.png\" alt=\"Drawing\" width=\"200\" height=\"200\"/>  \n",
    "\n",
    "Try editing the code in the cell below to calculate the muon's total momentum.\n",
    "\n",
    "When you're ready, scroll down to **Part 3**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 143
    },
    "id": "baDYUis1jv05",
    "outputId": "463414e4-b18d-48b2-c71d-7d1bde2d3885"
   },
   "outputs": [],
   "source": [
    "# You can specify a column by data['column name']\n",
    "# This adds a column and fills it with E+px^2 for each event\n",
    "data['oogieboogie'] = data['E'] + data['py']**2         # this calculates nothing useful at all\n",
    "data.head(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4MXlLaWNjv1F"
   },
   "source": [
    "## Part Three\n",
    "Here's where special relativity comes in. A particle's **energy** is due to both it's **mass** and **momentum**. This equation is also similar to the Pythagoream theorem:  \n",
    "\n",
    "<center>energy<sup>2</sup> = mass<sup>2</sup> + momentum<sup>2</sup></center>  \n",
    "\n",
    "*Mass* in this equation is called the particle's *invariant mass* because it is a value all observers would agree on, regardless of their reference frame. The code below creates a column for the mass of the muon, but the equation isn't quite right.  \n",
    "- Edit the code in the cell below to *really* calculate the muon's mass.  \n",
    "- How could your equation produce some imaginary numbers?  \n",
    "\n",
    "When you're ready, scroll down to **Part 4**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 143
    },
    "id": "1sgk3CLRjv1A",
    "outputId": "a9e001d2-d016-4920-ec41-483ba796bece"
   },
   "outputs": [],
   "source": [
    "data['mass'] = (data['E'] + data['px']) # this also calculates nothing useful\n",
    "data = data.fillna(0) # fills in a zero for any imaginary values your caluclation might produce\n",
    "data.head(3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "id": "-scsKIdKjv1G",
    "outputId": "2916b614-aea5-4393-82b5-ba3ab7f68847"
   },
   "outputs": [],
   "source": [
    "plt.hist(data['mass'], bins=10, range=[0,50], log=False)  # makes the histogram\n",
    "plt.title(\"Good enough Title\")\n",
    "plt.xlabel(\"x-axis label (in GeV/c^2)\")\n",
    "plt.ylabel(\"number of events\")\n",
    "plt.grid(False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zlF-KDk0jv1K"
   },
   "source": [
    "## Part Four\n",
    "The code above may take a few moments to run since it's grabbing a pretty big data set (>400,000 events!?). When it's finished, you'll see a histogram of the invariant mass values you calucalted for muon 1 in each event. This is called a *mass plot*.  \n",
    "- You'll probably need to adjust the histogram's range and number of bins to see a clear peak. The x-value of that peak occurs at the invariant mass of the particle (in this case, a muon).\n",
    "- Based on your histogram, what value does the muon's invariant mass seem to be?\n",
    "- Try replacing the histogram's title and x-axis label to something better.  \n",
    "\n",
    "## Part Five  \n",
    "Now that you've analyzed your own huge set of particle collision events, here are some follow-up questions:\n",
    "- How does the value you calculated compare to the accepted mass of the [muon](https://en.wikipedia.org/wiki/Muon)? Is it very different, if so why might that be?\n",
    "- Why doesn't energy have x, y, and z components like momentum?\n",
    "- How does E<sup>2</sup> - p<sup>2</sup> give you a unit of mass<sup>2</sup>?\n",
    "\n",
    "## More Tools  \n",
    "The cells below show some more advanced techniques.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "3Zq2yLC0jv1L"
   },
   "outputs": [],
   "source": [
    "# run this command to make your plots look like they're from xkcd.com\n",
    "plt.xkcd();\n",
    "# then re-execute your code to make a plot and see it xkcd-ified."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "SIbaFNC1jv1Q"
   },
   "outputs": [],
   "source": [
    "# run this cell to make normal-looking plots again\n",
    "mpl.rcdefaults()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Mj5aIoG4jv1U"
   },
   "source": [
    "A professional physicist would usually fit a theoretical curve to the mass plot to identify the particle's invariant mass. For reconstructing more complex particle decays they'd use a [relativistic Breit-Wigner](https://en.wikipedia.org/wiki/Relativistic_Breit%E2%80%93Wigner_distribution) curve, but we'll try it here just for practice.  \n",
    "- First, choose some parameters for the curve based on where the peak is on your histogram and how wide it looks.  \n",
    "- Then, run the code and keep adjusting the parameters until the curve tightly fits the peak."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "LOaKSIAWjv1U"
   },
   "outputs": [],
   "source": [
    "# set the Breit-Wigner parameters\n",
    "xmin = 0  # lower bound for your plot\n",
    "xmax = 1  # upper bound for your plot\n",
    "particle_mass = .5 # the x-value of the peak\n",
    "width = .02 # width of the peak halfway up; a.k.a. \"full width at half max\" or FWHM\n",
    "K = 15  # a constant that affects the height of the curve, you'll need to tinker with this some\n",
    "\n",
    "# calculate the Breit-Wigner curve\n",
    "x = np.arange(xmin, xmax, (xmax-xmin)/200) # makes a series of equally spaced x-values\n",
    "y = K / ((x**2 - particle_mass**2)**2 + (particle_mass*width)**2) # calculates the y-values for the B-W curve\n",
    "\n",
    "# make the plot\n",
    "fig, ax = plt.subplots(1,1, figsize=(6,4))\n",
    "ax.hist(data['mass'], bins=10, range=[xmin,xmax], log=False)  # plots the histogram\n",
    "ax.plot(x, y, color='r')  # plots the curve\n",
    "plt.title(\"title me!\")\n",
    "plt.xlabel(\"x-axis label (in GeV/c^2)\")\n",
    "plt.ylabel(\"number of events\")\n",
    "plt.grid(True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DcA5-9bPjv1a"
   },
   "source": [
    "If you finally have a curve that pretty closely fits the peak in your histogram, take a look at the parameters you set. We can use those to determine the particle's mass and lifetime. The invariant mass of the particle is your *particle_mass*. The particle's mean lifetime is h-bar/width.  \n",
    "- How do the values below compare with the accepted values for mass and lifetime of the muon?  \n",
    "They probably don't agree very well because we're applying the Breit-Wigner curve to a particle we observed (not valid) instead of one we're trying to reconstruct from analyzing its decay products. If you ever find yourself doing a *mass reconstruction* to discover a short-lived particle by its decay products, this will might be useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "nfxSGutbjv1a"
   },
   "outputs": [],
   "source": [
    "print(\"mass = \", np.round(particle_mass,6), \" GeV/c^2\")\n",
    "hbar = 6.6e-25      # in GeV*s\n",
    "print(\"mean lifetime = \", np.format_float_scientific(hbar/width,6), \" seconds\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "id": "IckJMDaijv1f"
   },
   "source": [
    "---  \n",
    "## Saving Your Work  \n",
    "This is running on a Google server on a distant planet and deletes what you've done when you close this tab. To save your work for later use or analysis you have a few options:  \n",
    "- File > \"Save a copy in Drive\" will save it to you Google Drive in a folder called \"Collaboratory\". You can run it later from there.  \n",
    "- File > \"Download .ipynb\" to save to your computer (and run with Jupyter software later)  \n",
    "- File > Print to ... um ... print.  \n",
    "- To save an image of a graph or chart, right-click on it and select Save Image as ...  \n",
    "\n",
    "## Credits\n",
    "This notebook was designed by [Quarknet](https://quarknet.org/) Teaching and Learning Fellow [Adam LaMee](https://adamlamee.github.io/) and UCF Physics undergrad Brooke Emison. The handy csv files were created from the CMS Run2011A primary datasets and converted from ROOT format by the masterful [Tom McCauley](https://github.com/tpmccauley). More can be found on the [CERN OpenData](http://opendata.cern.ch/?ln=en) site, like [here](http://opendata.cern.ch/record/545). The 3D vector image can be found on [WikiMedia Commons](https://commons.wikimedia.org/wiki/File:Coord_XYZ.svg). Finally, thanks to the great folks at [Binder](https://mybinder.org/) and [Google Colaboratory](https://colab.research.google.com/notebooks/intro.ipynb) for making this notebook interactive without you needing to download it or install [Jupyter](https://jupyter.org/) on your own device. Find more activities and license info at [CODINGinK12.org](http://www.codingink12.org)."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "3_Muon mass",
   "provenance": []
  },
  "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}