dimuonSpectrum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of the di-muon spectrum using data from the CMS detector\n",
    "\n",
    "This analysis takes data from the CMS experiment recorded in 2012 during Run B and C and extracts the di-muon spectrum. The di-muon 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 opens an interactive plot, which allows to zoom and navigate in the 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",
   "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",
   "metadata": {},
   "source": [
    "# Data analysis\n",
    "\n",
    "The following part of the notebook performs the data analysis of the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ROOT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enable multi-threading\n",
    "\n",
    "The default here is set to a single thread. You can choose the number of threads based on your system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ROOT.ROOT.EnableImplicitMT(1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create dataframe from NanoAOD files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = ROOT.RDataFrame(\"Events\", \"root://eospublic.cern.ch//eos/opendata/cms/derived-data/AOD2NanoAODOutreachTool/Run2012BC_DoubleMuParked_Muons.root\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select events with exactly two muons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_2mu = df.Filter(\"nMuon == 2\", \"Events with exactly two muons\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Select events with two muons of opposite charge"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_os = df_2mu.Filter(\"Muon_charge[0] != Muon_charge[1]\", \"Muons with opposite charge\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute invariant mass of the dimuon system\n",
    "\n",
    "The following code just-in-time compiles the C++ function to compute the invariant mass, so that the function can be called in the Define node of the ROOT dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ROOT.gInterpreter.Declare(\n",
    "\"\"\"\n",
    "using namespace ROOT::VecOps;\n",
    "float computeInvariantMass(RVec<float>& pt, RVec<float>& eta, RVec<float>& phi, RVec<float>& mass) {\n",
    "    ROOT::Math::PtEtaPhiMVector m1(pt[0], eta[0], phi[0], mass[0]);\n",
    "    ROOT::Math::PtEtaPhiMVector m2(pt[1], eta[1], phi[1], mass[1]);\n",
    "    return (m1 + m2).mass();\n",
    "}\n",
    "\"\"\")\n",
    "df_mass = df_os.Define(\"Dimuon_mass\", \"computeInvariantMass(Muon_pt, Muon_eta, Muon_phi, Muon_mass)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Book histogram of dimuon mass spectrum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bins = 30000 # 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",
    "hist = df_mass.Histo1D(ROOT.RDF.TH1DModel(\"\", \"\", bins, low, up), \"Dimuon_mass\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Request cut-flow report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report = df_mass.Report()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create canvas for plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ROOT.gStyle.SetOptStat(0)\n",
    "ROOT.gStyle.SetTextFont(42)\n",
    "c = ROOT.TCanvas(\"c\", \"\", 800, 700)\n",
    "c.SetLogx()\n",
    "c.SetLogy();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw histogram"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hist.GetXaxis().SetTitle(\"m_{#mu#mu} (GeV)\")\n",
    "hist.GetXaxis().SetTitleSize(0.04)\n",
    "hist.GetYaxis().SetTitle(\"N_{Events}\")\n",
    "hist.GetYaxis().SetTitleSize(0.04)\n",
    "hist.SetStats(False)\n",
    "hist.Draw();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw labels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "label = ROOT.TLatex()\n",
    "label.SetTextAlign(22)\n",
    "label.DrawLatex(0.55, 3.0e4, \"#eta\")\n",
    "label.DrawLatex(0.77, 7.0e4, \"#rho,#omega\")\n",
    "label.DrawLatex(1.20, 4.0e4, \"#phi\")\n",
    "label.DrawLatex(4.40, 1.0e5, \"J/#psi\")\n",
    "label.DrawLatex(4.60, 1.0e4, \"#psi'\")\n",
    "label.DrawLatex(12.0, 2.0e4, \"Y(1,2,3S)\")\n",
    "label.DrawLatex(91.0, 1.5e4, \"Z\")\n",
    "label.SetNDC(True)\n",
    "label.SetTextAlign(11)\n",
    "label.SetTextSize(0.04)\n",
    "label.DrawLatex(0.10, 0.92, \"#bf{CMS Open Data}\")\n",
    "label.SetTextAlign(31)\n",
    "label.DrawLatex(0.90, 0.92, \"#sqrt{s} = 8 TeV, L_{int} = 11.6 fb^{-1}\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Draw interactive canvas"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "c.Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Print cut-flow report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "report.Print()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}