taggerComp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Signal files\n",
    "First let's look at the available signal files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash \n",
    "export XrdSecGSISRVNAMES=cmseos.fnal.gov\n",
    "xrdfs root://cmseos.fnal.gov/ ls /eos/uscms/store/group/lpcjme/noreplica/NanoHRT/Trees/Apr08/forMCComparions/reweighted/"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Open signal files and inspect tree\n",
    "Let's open the one signal file and one background file and print the contents of the tree"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ROOT as rt\n",
    "\n",
    "filename_hh = 'root://cmseos.fnal.gov//eos/uscms/store/group/lpcjme/noreplica/NanoHRT/Trees/Apr08/forMCComparions/reweighted/BulkGravTohh_tree_noc.root'\n",
    "filename_qcd = 'root://cmseos.fnal.gov//eos/uscms/store/group/lpcjme/noreplica/NanoHRT/Trees/Apr08/forMCComparions/reweighted/qcd-pythia_tree_SLIM.root'\n",
    "\n",
    "file_hh = rt.TFile.Open(filename_hh)\n",
    "file_qcd = rt.TFile.Open(filename_qcd)\n",
    "\n",
    "tree_hh = file_hh.Get('Friends')\n",
    "tree_qcd = file_qcd.Get('Friends')\n",
    "tree_hh.Print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot tagger distributions\n",
    "Let's plot each of the different discriminatig variables distributions "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%jsroot on\n",
    "\n",
    "tagger = 'ak8_sdmass'\n",
    "#tagger = 'ak8_tau21'\n",
    "#tagger = 'ak8_ecfN2'\n",
    "#tagger = 'ak8_doubleb'\n",
    "#tagger = 'ak8_nn_HbbvsQCD'\n",
    "#tagger = 'ak8_decorr_nn_ZHbbvsQCD'\n",
    "#tagger = 'ak8_bestH'\n",
    "\n",
    "cuts = 'gen_pt > 1000 && gen_pt < 1500 && abs(gen_eta) < 2.4 && ak8_sdmass > 40 && ak8_sdmass < 200'\n",
    "\n",
    "hist_hh = rt.TH1D('hist_hh','hist_hh',85,30,200)\n",
    "hist_qcd = rt.TH1D('hist_qcd','hist_qcd',85,30,200)\n",
    "tree_hh.Project('hist_hh',tagger,'rewgt*(%s)'%cuts)\n",
    "tree_qcd.Project('hist_qcd',tagger,'genweight*(%s)'%cuts)\n",
    "\n",
    "c = rt.TCanvas('c','c',500,400)\n",
    "hist_hh.Scale(1/hist_hh.Integral())\n",
    "hist_qcd.Scale(1/hist_qcd.Integral())\n",
    "hist_qcd.SetLineColor(rt.kBlue)\n",
    "hist_hh.SetLineColor(rt.kRed)\n",
    "if hist_qcd.GetMaximum() > hist_hh.GetMaximum():\n",
    "    hist_qcd.Draw('h')\n",
    "    hist_hh.Draw('hsame')\n",
    "else:\n",
    "    hist_hh.Draw('h')\n",
    "    hist_qcd.Draw('hsame')\n",
    "c.Draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# checkpoint: \n",
    "plot the output of the ML-based algorithms: e.g. for the Higgs case for DoubleB, DeepAK8 and DeepAK8-MD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## place holder for the checkpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Convert tree to dataframe with `uproot`\n",
    "Open the file with `uproot` and convert the tree to a `pandas` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## this step may take some time \n",
    "\n",
    "import uproot\n",
    "#!pip install lz4 --user\n",
    "#!pip install xxhash --user\n",
    "\n",
    "uptree_hh = uproot.open(filename_hh)[\"Friends\"]\n",
    "uptree_qcd = uproot.open(filename_qcd)[\"Friends\"]\n",
    "\n",
    "branches = ['ak8_nn_HbbvsQCD',\n",
    "            'ak8_decorr_nn_ZHbbvsQCD',\n",
    "            'ak8_doubleb',\n",
    "            'ak8_tau21',\n",
    "            'ak8_ecfN2',\n",
    "            'ak8_bestH',\n",
    "            'ak8_sdmass',\n",
    "            'ak8_pt',\n",
    "            'ak8_eta',\n",
    "            'gen_pt',\n",
    "            'gen_eta']\n",
    "\n",
    "\n",
    "print(\"get hh\")\n",
    "df_hh = uptree_hh.pandas.df(branches=branches+['rewgt'])\n",
    "#print(df_hh)\n",
    "\n",
    "print(\"get QCD\")\n",
    "df_qcd = uptree_qcd.pandas.df(branches=branches+['genweight'])\n",
    "#print(df_qcd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Apply cuts and merge dataframes\n",
    "Let's apply the same fiducial cuts and merge the signal and background dataframes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas\n",
    "import numpy as np\n",
    "\n",
    "# mask higgs outside of fiducial cuts\n",
    "mask_hh = (df_hh['gen_pt'] > 1000) & (df_hh['gen_pt'] < 1500) & (np.abs(df_hh['gen_eta']) < 2.4) & (df_hh['ak8_sdmass'] > 40) & (df_hh['ak8_sdmass'] < 200)\n",
    "df_hh = df_hh[mask_hh]\n",
    "df_hh['sample_weight'] = df_hh['rewgt']\n",
    "df_hh['label'] = np.ones((len(df_hh)),dtype=int)\n",
    "\n",
    "# mask qcd outside of fiducial cuts\n",
    "mask_qcd = (df_qcd['gen_pt'] > 1000) & (df_qcd['gen_pt'] < 1500) & (np.abs(df_qcd['gen_eta']) < 2.4) & (df_qcd['ak8_sdmass'] > 40) & (df_qcd['ak8_sdmass'] < 200)\n",
    "df_qcd = df_qcd[mask_qcd]\n",
    "df_qcd['sample_weight'] = df_qcd['genweight']\n",
    "df_qcd['label'] = np.zeros((len(df_qcd)),dtype=int)\n",
    "\n",
    "# concatenate signal and background\n",
    "df = pandas.concat([df_hh,df_qcd])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot ROCs\n",
    "Plot the ROC for one of the tagging variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from sklearn.metrics import roc_curve, auc\n",
    "%matplotlib inline\n",
    "fpr = {}\n",
    "tpr = {}\n",
    "threshold = {}\n",
    "plt.figure()\n",
    "\n",
    "tagger = 'ak8_sdmass'\n",
    "fpr[tagger], tpr[tagger], threshold[tagger] = roc_curve(df['label'], df[tagger], sample_weight=df['sample_weight'])\n",
    "myauc = auc(fpr[tagger], tpr[tagger],reorder=True)\n",
    "plt.plot(tpr[tagger],fpr[tagger],label='%s tagger, AUC = %.3f'%(tagger, myauc))\n",
    "\n",
    "plt.legend(loc='best')\n",
    "plt.semilogy()\n",
    "plt.xlim([0, 1])\n",
    "plt.ylim([1e-4, 1])\n",
    "plt.xlabel(\"TPR\");\n",
    "plt.ylabel(\"FPR\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# checkpoint: \n",
    "Compare performance of other tagging algorithms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## place holder for checkpoint"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Beyond ROC curves: Tagger's response vs. jet mass \n",
    "Study the correlation between the tagger's response and the jet mass\n",
    "Different WP for each algorithm will be provided by the instructor\n",
    "Discussion: thoughts, ideas, concerns ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## place holder for the exercise"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Train your own simple tagger\n",
    "Train your own simple tagger with keras"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import keras\n",
    "from keras.models import Model\n",
    "from keras.layers import Input, Activation, Dense\n",
    "\n",
    "inputs = Input(shape=(2,), name = 'input')  \n",
    "x = Dense(5, name = 'dense_1', activation = 'relu')(inputs)\n",
    "outputs = Dense(1, name = 'output', activation='sigmoid')(x)\n",
    "\n",
    "# creae the model\n",
    "model = Model(inputs=inputs, outputs=outputs)\n",
    "# compile the model\n",
    "model.compile(optimizer='adam', loss='binary_crossentropy', metrics=['accuracy'])\n",
    "# print the model summary\n",
    "model.summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.vstack([df['ak8_sdmass'], df['ak8_tau21']]).T\n",
    "print(x.shape)\n",
    "y = df['label']\n",
    "print(y.shape)\n",
    "model.fit(x,y, epochs=1)\n",
    "df['ak8_sdmass_tau21'] = model.predict(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    " - Exercise 1. Plot all of the taggers on the same ROC curve\n",
    " - Exercise 2. Train your own custom \"traditional\" tagger based on mSD and tau21\n",
    " - Exercise 3. Add your tagger to the ROC curve\n",
    " - Exercise 4. Repeat with a different signal sample (W, Z, or top)\n",
    " - Exercise 5 [advanced]. Understand the impact of different sets of inputs in the jet classification problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "nanohrt-hats",
   "language": "python",
   "name": "nanohrt-hats"
  },
  "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.15+"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}