dimu/dimuon-Julia.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1beb0e7a-ad3c-4f5d-acba-701fe0b2feb8",
   "metadata": {},
   "source": [
    "# Analysis of the di-muon spectrum using data from the CMS detector\n",
    "\n",
    "This analysis extracts the mass spectrum of di-muons produced in proton-proton collisions at sqrt(s)=7TeV data from the CMS experiment data recorded in 2012 during Run B and C. It is implemented using the [Julia](https://julialang.org/) language and can run in a Jupyter notebook. Python and C++ version of this analysis can be found [here](https://opendata.web.cern.ch/record/12342).\n",
    "\n",
    "The spectrum is computed from the data by calculating the invariant mass of muon pairs with opposite charge. In the resulting plot, you are able to rediscover particle resonances in a wide energy range from the [eta meson](https://en.wikipedia.org/wiki/Eta_meson) at about 548 MeV up to the [Z boson](https://en.wikipedia.org/wiki/W_and_Z_bosons) at about 91 GeV.\n",
    "\n",
    "The analysis code produces a plot with the dimuon spectrum. Note that the bump at 30 GeV is not a resonance but an effect of the data taking due to the used trigger. The technical description of the dataset can be found in the respective record linked below.\n",
    "\n",
    "The result of this analysis can be compared with [an official result of the CMS collaboration using data taken in 2010](https://cds.cern.ch/record/1456510), see the plot below:\n",
    "\n",
    "![](http://cds.cern.ch/record/1456510/files/pictures_samples_dimuonSpectrum_40pb-1_mod-combined.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77ca9f8c-191b-40b2-a304-cd4ff56d0398",
   "metadata": {},
   "source": [
    "# Dataset description\n",
    "\n",
    "The dataset consists of the following columns.\n",
    "\n",
    "| Column name | Data type | Description |\n",
    "|-------------|-----------|-------------|\n",
    "| `nMuon` | `unsigned int` | Number of muons in this event |\n",
    "| `Muon_pt` | `float[nMuon]` | Transverse momentum of the muons (stored as an array of size `nMuon`) |\n",
    "| `Muon_eta` | `float[nMuon]` | Pseudorapidity of the muons |\n",
    "| `Muon_phi` | `float[nMuon]` | Azimuth of the muons |\n",
    "| `Muon_mass` | `float[nMuon]` | Mass of the muons |\n",
    "| `Muon_charge` | `int[nMuon]` | Charge of the muons (either 1 or -1) |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e8887f9",
   "metadata": {},
   "source": [
    "## Software prerequeries\n",
    "\n",
    "This notebook needs a [Julia kernel](https://julialang.github.io/IJulia.jl/stable/) to be executed in Jupyter.\n",
    "\n",
    "In addition it uses the Julia packages NBInclude (used to include an auxiliary notebook), UnROOT (used to read the input data file), Gafly (used for plotting). The following commands will install these two packages on the system if they are not installed yet."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "743d447c",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\u001b[32m\u001b[1m    Updating\u001b[22m\u001b[39m registry at `~/.julia/registries/General.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "┌ Warning: The active manifest file at `/home/pgras/.julia/environments/v1.7/Manifest.toml` has an old format that is being maintained.\n",
      "│ To update to the new format run `Pkg.upgrade_manifest()` which will upgrade the format without re-resolving.\n",
      "└ @ Pkg.Types /home/pgras/git.d/julia/usr/share/julia/stdlib/v1.7/Pkg/src/manifest.jl:287\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Manifest.toml`\n",
      "┌ Warning: The active manifest file is an older format with no julia version entry. Dependencies may have been resolved with a different julia version.\n",
      "└ @ nothing /home/pgras/.julia/environments/v1.7/Manifest.toml:0\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Manifest.toml`\n",
      "\u001b[32m\u001b[1m   Resolving\u001b[22m\u001b[39m package versions...\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Project.toml`\n",
      "\u001b[32m\u001b[1m  No Changes\u001b[22m\u001b[39m to `~/.julia/environments/v1.7/Manifest.toml`\n"
     ]
    }
   ],
   "source": [
    "import Pkg\n",
    "Pkg.add(\"NBInclude\")\n",
    "Pkg.add(\"UnROOT\")\n",
    "Pkg.add(\"Gadfly\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4d532d88",
   "metadata": {},
   "source": [
    "Import the packages used in this notebook:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "681332c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "using NBInclude\n",
    "using UnROOT"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2e2264ef",
   "metadata": {},
   "source": [
    "## Some utilities we will need\n",
    "\n",
    "Define four-momentum and histogram tools"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "1586f179-4745-40a9-9b0a-39181b5dc3f0",
   "metadata": {},
   "outputs": [],
   "source": [
    "@nbinclude(\"JuliaUtils.ipynb\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33185ba2-9f9e-4f0a-8d5d-b8942602163c",
   "metadata": {},
   "source": [
    "## The data file\n",
    "\n",
    "To run this example, we need to download the CMS data file from [here](https://opendata.cern.ch/record/12341/files/Run2012BC_DoubleMuParked_Muons.root). If you don't have the command `wget` on your system you can download the file using your preferred web browser and save it with the name `Run2012BC_DoubleMuParked_Muons.root` in the directory containing this notebook.\n",
    "\n",
    "The file to download is 1.2 GByte big."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "3368a5f0-2f28-4355-b42c-2dee3aa507fc",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "--2021-11-16 09:14:40--  https://opendata.cern.ch/record/12341/files/Run2012BC_DoubleMuParked_Muons.root\n",
      "Resolving opendata.cern.ch (opendata.cern.ch)... 188.184.28.138, 188.184.93.89, 188.185.82.144, ...\n",
      "Connecting to opendata.cern.ch (opendata.cern.ch)|188.184.28.138|:443... connected.\n",
      "HTTP request sent, awaiting response... 416 Requested Range Not Satisfiable\n",
      "\n",
      "    The file is already fully retrieved; nothing to do.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "run(`wget -c --progress dot:giga -O Run2012BC_DoubleMuParked_Muons.root \"https://opendata.cern.ch/record/12341/files/Run2012BC_DoubleMuParked_Muons.root\"`);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "31e43a01",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"Run2012BC_DoubleMuParked_Muons.root\""
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fname = \"Run2012BC_DoubleMuParked_Muons.root\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de45313d",
   "metadata": {},
   "source": [
    "## The analysis\n",
    "\n",
    "The function that will analyse the events. We select events with exactly two muons with opposite charges, compute the mass of the \"dimuon\" object formed by the two objects, and fill an histogram that will represent the dimuon mass spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d2897de7",
   "metadata": {},
   "outputs": [],
   "source": [
    "function analyze_tree(t, maxevents = -1)\n",
    "    bins = 30_000 # Number of bins in the histogram\n",
    "    low = 0.25 # Lower edge of the histogram\n",
    "    up = 300.0 # Upper edge of the histogram\n",
    "    h = H1{Float64}(Axis(bins, low, up))\n",
    "    for (ievt, evt) in enumerate(t)\n",
    "        maxevents >= 0 && ievt > maxevents && break\n",
    "        evt.nMuon ==2 || continue\n",
    "        evt.Muon_charge[1] != evt.Muon_charge[2] || continue\n",
    "        dimuon_mass = m(ptetaphim(evt.Muon_pt[1], evt.Muon_eta[1], evt.Muon_phi[1], evt.Muon_mass[1])\n",
    "                        + ptetaphim(evt.Muon_pt[2], evt.Muon_eta[2], evt.Muon_phi[2], evt.Muon_mass[2]))\n",
    "        hfill!(h, dimuon_mass)\n",
    "    end\n",
    "    h\n",
    "end;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "98f94ddc",
   "metadata": {},
   "source": [
    "We will now run the function on the data stored in the input `fname` `ROOT` file. It will take about 40 seconds to run on the 61.5 million of events contained in the file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "adf3b121",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " 32.669432 seconds (165.32 M allocations: 22.806 GiB, 12.62% gc time, 2.41% compilation time)\n"
     ]
    }
   ],
   "source": [
    "t = LazyTree(ROOTFile(fname),\"Events\")\n",
    "@time h = analyze_tree(t);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c9632ca-be11-4583-a3c5-a3ef31549938",
   "metadata": {},
   "source": [
    "# Let's plot the result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "small-faith",
   "metadata": {
    "code_folding": [],
    "run_control": {
     "marked": false
    }
   },
   "outputs": [
    {
     "ename": "LoadError",
     "evalue": "UndefVarError: xedges not defined",
     "output_type": "error",
     "traceback": [
      "UndefVarError: xedges not defined",
      "",
      "Stacktrace:",
      " [1] top-level scope",
      "   @ In[1]:2",
      " [2] eval",
      "   @ ./boot.jl:373 [inlined]",
      " [3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)",
      "   @ Base ./loading.jl:1196"
     ]
    }
   ],
   "source": [
    "import Gadfly as gf\n",
    "p = gf.plot(x=xedges(h), y=vcat(h.sumw[2:end-1], [h.sumw[end-1]]), gf.Geom.step, \n",
    "        gf.Scale.x_log10(minvalue=0.25, maxvalue=300.), gf.Scale.y_log10,\n",
    "        gf.Guide.xlabel(\"Dimuon mass [GeV/c²]\"), gf.Guide.ylabel(\"Event count\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6e8c251",
   "metadata": {},
   "source": [
    "[Back to the presentation](../01-Julia-dream-JIII2021.ipynb)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.7.0-rc1",
   "language": "julia",
   "name": "julia-1.7"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}