Open-Data-with-CMS-outreach-and-education.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open Data with CMS - outreach and education"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### What is Open Data and why is it important?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = \"https://yt3.ggpht.com/a-/ACSszfENmQ2SCS7AnguI8q_ZZlNl2Ne_CC0n2lA_xA=s900-mo-c-c0xffffffff-rj-k-no\" align = \"right\"  style = \"height:200px\" >\n",
    "\n",
    "CMS produces a huge amount of data and releases it publicly in a manner very few other experiments at the top of their relative fields do. On one hand that data is very easily available to the general public, on the other it is a bit cryptic to use without someone experienced at hand to explain a few things. Here we aim to present some tools to make it very easy and fast to use.\n",
    "\n",
    "Obviously the data is very important to the scientists or we wouldn't be doing this experiment. It is also very valuable that it can be shared with others: it can make it easier for non-scientists to understand what your work is about, it can give educators interesting tools that promote modern research, it can hook a curious student to get into sciences and so forth. Especially in education there is this problem of \"relevance\" and \"authenticity\", ie. how do things done in school relate to the students' lives or to what is actually happening in scientific research today, and using actual data from actual research (from experiments obviously unperformable in classrooms) might prove to be an important asset in closing that gap. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### So how do we do it?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src = \"https://cdn-images-1.medium.com/max/346/0*I3hkRieQ6B3qwwhy.\" align = \"right\" style = \"height:200px\">\n",
    "\n",
    "There are a few ways. You could use event visualizers, spreadsheet programs, particle physics masterclass materials and what not, all of which can be found through [the Open Data portal](http://opendata.cern.ch/), but the one we'll focus on here is Jupyter. You are currently reading a Jupyter notebook, an interactive file that can trivially incorporate text, pictures, animations and most importantly, code. It can support many languages, but we'll use Python 3 here as it is intuitive enough to be used even with people who have absolutely no programming background.\n",
    "\n",
    "Best thing is that you can add to it on the fly and let people explore the data as they want without them needing to install anything. Especially in education: people hate it when they have to install things and learn to use a bunch of programs they might use every now and then (or just once), teenagers even more so. With Jupyter you don't have to: we can use [MyBinder](www.mybinder.org) to create a virtual workspace and get to work in a minute or two without any hassle, completely browser-based. There are also others, like Google's [colab](https://colab.research.google.com), if it suits your preferences better.  \n",
    "\n",
    "CERN's own [SWAN](https://swan.web.cern.ch/) (Service for Web based ANalysis) is also be based on Jupyter, though it requires you to have CERNbox.\n",
    "\n",
    "OF NOTE, this document is meant for Jupyter and Python 3, especially with Binder, so other platforms may require a bit of juggling to get things right. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Right, let's get some examples on the table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see (and test by double clicking around), this document is made of *cells*. Those can have either text (in markdown format) or code in them. To run a cell, select it (so it gets green) and press *ctrl + Enter*. On the left side of the code cells there's a [ ]. If it is empty, the cell has not been run. If it shows an asterisk, the cell is running. If there is a number, it has been run.\n",
    "\n",
    "If the system somehow gets stuck, just select *Kernel $\\rightarrow$ Restart & Clear Output* from the upper menus. This clears any weird leftovers or hangups that might be left in the kernel. If otherwise unsure of commands, press **h** when no cells are chosen to bring up key bindings.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This notebook runs on Python 3, so we need to import some modules.\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.stats import norm\n",
    "\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When that is done, let's get some data in here. We'll use the 2012 sample and a dataset comprising of events where exactly two muons were detected without worrying about any trigger selections (using code similar to the one used in [record 552](http://opendata.cern.ch/record/552)).\n",
    "\n",
    "This is an example of a simplified derived dataset, but you could very well use any set you want to make from the openly published results. All derived sets can be found from [here](\n",
    "http://opendata.cern.ch/search?page=1&size=20&q=&experiment=CMS&type=Dataset&subtype=Derived)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "muons = pd.read_csv('http://opendata.cern.ch/record/545/files/Dimuon_DoubleMu.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Getting a picture of \"what does particle physics data actually look like\" can be a bit hard for those who don't practise it. It's easy to show them, just run the cell below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This shows top 5 rows of data, though you can put a desired number in the parentheses.\n",
    "\n",
    "muons.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Okay, seems legit. That is still a bit dense if left up to the reader, but with a bit of helpful explaining or physics background it starts to make sense. Importantly we can show that there is quite a lot of it, which means it's going to require computers and programming to analyze it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This shows the number of rows in the data.\n",
    "\n",
    "len(muons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Though that number seems small to people who are used to these experiments, it is still a huge number for others. It also doesn't tell us anything, whereas a good graphical visualization does. So let's plot a histogram."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 10))\n",
    "\n",
    "plt.hist(muons.M, 300, range = (0,150))\n",
    "\n",
    "plt.xlabel('Invariant mass (GeV/$c^2$)', fontsize = 15)\n",
    "plt.ylabel('Number of events \\n', fontsize = 15)\n",
    "plt.title('Measured distribution of events with two muons \\n', fontsize = 15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As a sidenote: if you're using a dataset that doesn't have invariant mass already computed in it, you could just as easily create a new variable from the existing data and use that for plotting."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "invMassSqrd = (muons.E1 + muons.E2)**2 - ((muons.px1 + muons.px2)**2 +(muons.py1 + muons.py2)**2  +(muons.pz1 + muons.pz2)**2)\n",
    "invMass = np.sqrt(invMassSqrd[(invMassSqrd >0)])\n",
    "\n",
    "invMass.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we see how the data plots out, we can focus on using it to show different things to the intended audience. Say you want to show them the bump at 90 GeV, for an example. Let's extract those results that fall within 80-100 GeV."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bump = muons[(muons.M >= 80) & (muons.M <= 100)]\n",
    "\n",
    "fig = plt.figure(figsize=(15, 10))\n",
    "\n",
    "plt.hist(bump.M, 300, range = (80,100))\n",
    "\n",
    "plt.xlabel('Invariant mass (GeV/$c^2$)', fontsize = 15)\n",
    "plt.ylabel('Number of events \\n', fontsize = 15)\n",
    "plt.title('Measured distribution of events with two muons \\n', fontsize = 15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's follow with a question \"how do these correlate with energies the measured particles have\"? You can set the desired threshold for high $p_t$ in the code below and see how it changes the plot. At which point does the low part stop to contribute?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "threshold = 30\n",
    "\n",
    "highE = bump[(bump.pt1 >= threshold) & (bump.pt2 >= threshold)]\n",
    "lowE = bump[(bump.pt1 < threshold) & (bump.pt2 < threshold)]\n",
    "\n",
    "fig = plt.figure(figsize=(15, 10))\n",
    "\n",
    "plt.hist(highE.M, 300, range = (80,100), alpha = 0.5 , color = 'black', label = 'High energy')\n",
    "plt.hist(lowE.M, 300, range = (80,100), alpha = 0.5 , color = 'red', label = 'Low energy')\n",
    "\n",
    "plt.xlabel('Invariant mass (GeV/$c^2$)', fontsize = 15)\n",
    "plt.ylabel('Number of events \\n', fontsize = 15)\n",
    "plt.title('Measured distribution of events with two muons \\n', fontsize = 15)\n",
    "plt.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also try to find out how the measurements are spread spatially."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "prapMuons = pd.concat([muons['eta1'],muons['eta2']])\n",
    "degMuons = prapMuons.copy()\n",
    "degMuons[:]=[2*np.arctan(np.exp(x))*360/(2*np.pi) for x in degMuons]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.hist(prapMuons, bins=100, range=(-3,3))\n",
    "\n",
    "plt.xlabel('Pseudorapidity', fontsize=15)\n",
    "plt.ylabel('Number of events', fontsize=15)\n",
    "plt.title('Spatial distribution \\n', fontsize=15)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.hist(degMuons, bins=180, range=(0,180))\n",
    "\n",
    "plt.xlabel('Degrees', fontsize=15)\n",
    "plt.ylabel('Number of events', fontsize=15)\n",
    "plt.title('Spatial distribution \\n', fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There is certainly a dip in the middle. Could it have something to do with the $p_t$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "transverse = pd.concat([muons.pt1,muons.pt2])\n",
    "plt.scatter(prapMuons, transverse, s=0.5)\n",
    "\n",
    "plt.title('$p_t$ in regard to $\\eta$ \\n', fontsize=15)\n",
    "plt.xlabel('Pseudorapidity', fontsize=15)\n",
    "plt.ylabel('Transverse momentum (GeV)', fontsize=15)\n",
    "\n",
    "plt.show()\n",
    "\n",
    "plt.scatter(prapMuons, transverse, s=0.5)\n",
    "\n",
    "plt.ylim(0,6)\n",
    "\n",
    "plt.title('$p_t$ in regard to $\\eta$ \\n', fontsize=15)\n",
    "plt.xlabel('Pseudorapidity', fontsize=15)\n",
    "plt.ylabel('Transverse momentum (GeV)', fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could also use pseudorapidity to see how accurately we can measure things dependant on the angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "high = 1.5\n",
    "low = 0.88\n",
    "\n",
    "# Check the whole data.\n",
    "\n",
    "highEta = muons[(abs(muons.eta1) >= high) & (abs(muons.eta2) >= high)]\n",
    "lowEta = muons[(abs(muons.eta1) < low) & (abs(muons.eta2) < low)]\n",
    "\n",
    "# Or uncomment and run this to check the 90 GeV bump (which was defined a bit above).\n",
    "\n",
    "# highEta = bump[(abs(bump.eta1) >= high) & (abs(bump.eta2) >= high)]\n",
    "# lowEta = bump[(abs(bump.eta1) < low) & (abs(bump.eta2) < low)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 8))\n",
    "\n",
    "h = len(highEta)\n",
    "l = len(lowEta)\n",
    "\n",
    "plt.hist(highEta.M, bins=300, range=(0,100), alpha=0.5, label='High $\\eta$, n = %i' %h, facecolor='black')\n",
    "plt.hist(lowEta.M, bins=300, range=(0,100), alpha=0.5, label='Low $\\eta$, n = %i' %l)\n",
    "\n",
    "# Adjust range and ylim accordingly, if you want to focus on the bump.\n",
    "\n",
    "plt.ylim(0,1500)\n",
    "\n",
    "plt.xlabel('Invariant mass [GeV/c²]', fontsize=15)\n",
    "plt.ylabel('Number of events', fontsize=15)\n",
    "plt.title('Detected muons in regard to $\\eta$ \\n', fontsize=15)\n",
    "plt.legend (loc='upper right', fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While more proper and complex curve fits take a bit more effort, you can also pretty easily show some rough distributions on the fly. For a more rigorous show you could use something like ROOT, but that can get a bit burdensome for some situations. Let's take a look at the J/$\\psi$ peak with nimble Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(15, 10))\n",
    "\n",
    "lower = 2.8\n",
    "upper = 3.4\n",
    "\n",
    "fit_i = 3.05\n",
    "fit_f = 3.15\n",
    "\n",
    "jpsi = muons[(muons.M >= lower) & (muons.M <= upper)]\n",
    "fit = muons[(muons.M >= fit_i) & (muons.M <= fit_f)]\n",
    "\n",
    "coefficient = len(fit)/len(jpsi)\n",
    "(mu, sigma) = norm.fit(fit.M)\n",
    "\n",
    "n, bins, patches = plt.hist(jpsi.M, 300, density = 1, facecolor = 'b', alpha=0.5, histtype='stepfilled',\n",
    "                            range=(lower, upper))\n",
    "\n",
    "y = coefficient*norm.pdf( bins, mu, sigma)\n",
    "l = plt.plot(bins, y, 'k--', linewidth = 6)\n",
    "\n",
    "plt.xlabel('Invariant mass (GeV/$c^2$)', fontsize = 15)\n",
    "plt.ylabel('Probability \\n', fontsize = 15)\n",
    "plt.title( 'J/$\\psi$ normed to one   $\\mu$ = %f , $\\sigma$ = %f \\n' %(mu,sigma), fontsize=15)\n",
    "\n",
    "cur_axes = plt.gca()\n",
    "cur_axes.axes.get_yaxis().set_ticks([])\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feedback & thanks for all the fish "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thank you for participating!\n",
    "\n",
    "If you could spare a minute or two and give us some feedback [here](https://goo.gl/forms/qE5g0dzQyLNKJvZw1), we'd appreciate it a lot."
   ]
  }
 ],
 "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.6.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}