Booster trainer.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import train_test_split\n",
    "import xgboost as xgb\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "import uproot\n",
    "import pickle\n",
    "import nue_booster\n",
    "import importlib\n",
    "importlib.reload(nue_booster)\n",
    "\n",
    "import awkward"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import localSettings as ls\n",
    "print(ls.ntuple_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "VARLOAD = [\"shr_tkfit_npointsvalid\",\"shr_tkfit_npoints\", # track-fitted shower\n",
    "           \"shrsubclusters0\",\"shrsubclusters1\",\"shrsubclusters2\", # sub-clusters for shower\n",
    "           \"trk_llr_pid_score_v\", # PID for proton candidate\n",
    "           \"shrmoliereavg\",\"shrmoliererms\", # moliere metrics for shower\n",
    "           \"trkshrhitdist2\", # 2D distance shower hits to track hits\n",
    "           \"trk_id\",\"trk_energy_tot\",\"shr_energy_tot_cali\",\n",
    "           \"tksh_distance\",\"shr_tkfit_dedx_Y\",\"shr_tkfit_dedx_U\",\"shr_tkfit_dedx_V\",\"tksh_angle\",\"n_showers_contained\",\"shr_score\", \n",
    "           \"nu_e\",\"slpdg\",# truth variables\n",
    "           \"category\",\"selected\",\n",
    "           #\"crtveto\",\"crthitpe\",\"_closestNuCosmicDist\", # CRT quantities\n",
    "           \"hits_ratio\",\"slclustfrac\",\"shr_energy\",\"shr_theta\",\"CosmicIP\",\"n_tracks_contained\"\n",
    "          ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fold = \"nuselection\"\n",
    "tree = \"NeutrinoSelectionFilter\"\n",
    "\n",
    "EXT3 = 'data_extbnb_mcc9.1_v08_00_00_25_reco2_G1_all_reco2.root'\n",
    "NU3  = 'prodgenie_bnb_nu_uboone_overlay_mcc9.1_v08_00_00_26_filter_run3_reco2_G_reco2.root'\n",
    "NUE3 = 'prodgenie_bnb_intrinsice_nue_uboone_overlay_mcc9.1_v08_00_00_26_run3_reco2_reco2.root'\n",
    "PI03 = 'prodgenie_nc_pi0_uboone_overlay_mcc9.1_v08_00_00_26_run3_G_reco2.root'\n",
    "\n",
    "EXT1 = 'run1_beam_off_cc0pinp.root'\n",
    "NU1  = 'run1_bnb_nu_cc0pinp.root'\n",
    "NUE1 = 'run1_bnb_nue_cc0pinp.root'\n",
    "PI01 = 'run1_bnb_ncpi0_cc0pinp.root'\n",
    "\n",
    "mc3 = uproot.open(ls.ntuple_path+NU3)[fold][tree]\n",
    "nc3 = uproot.open(ls.ntuple_path+PI03)[fold][tree]\n",
    "nue3 = uproot.open(ls.ntuple_path+NUE3)[fold][tree]\n",
    "ext3 = uproot.open(ls.ntuple_path+EXT3)[fold][tree]\n",
    "\n",
    "mc1 = uproot.open(ls.ntuple_path+NU1)[fold][tree]\n",
    "nc1 = uproot.open(ls.ntuple_path+PI01)[fold][tree]\n",
    "nue1 = uproot.open(ls.ntuple_path+NUE1)[fold][tree]\n",
    "ext1 = uproot.open(ls.ntuple_path+EXT1)[fold][tree]\n",
    "\n",
    "uproot_v = [mc3,nc3,nue3,ext3,mc1,nc1,nue1,ext1]\n",
    "\n",
    "nc3 = nc3.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "mc3 = mc3.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "nue3 = nue3.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "ext3 = ext3.pandas.df(VARLOAD, flatten=False)\n",
    "nc1 = nc1.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "mc1 = mc1.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "nue1 = nue1.pandas.df(VARLOAD + [\"weightSpline\"], flatten=False)\n",
    "ext1 = ext1.pandas.df(VARLOAD, flatten=False)\n",
    "\n",
    "ext3[\"weightSpline\"] = 1\n",
    "ext1[\"weightSpline\"] = 1\n",
    "\n",
    "df_v = [mc3,nc3,nue3,ext3,mc1,nc1,nue1,ext1]\n",
    "\n",
    "# define some energy-related variables\n",
    "for i,df in enumerate(df_v):\n",
    "    df[\"reco_e\"] = (df[\"shr_energy_tot_cali\"] + 0.03) / 0.79 + df[\"trk_energy_tot\"]\n",
    "    df[\"reco_e_qe\"] = 0.938*((df[\"shr_energy\"]+0.030)/0.79)/(0.938 - ((df[\"shr_energy\"]+0.030)/0.79)*(1-np.cos(df[\"shr_theta\"])))\n",
    "    df[\"reco_e_rqe\"] = df[\"reco_e_qe\"]/df[\"reco_e\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# how to get the LLR-PID value for the \"track candidate\" (proton for nue selection, muon for numu)\n",
    "# can be done for any variable\n",
    "# code from Giuseppe!\n",
    "\n",
    "for i,df in enumerate(df_v):\n",
    "    up = uproot_v[i]\n",
    "    trk_llr_pid_v = up.array('trk_llr_pid_score_v')\n",
    "    trk_id = up.array('trk_id')-1 # I think we need this -1 to get the right result\n",
    "    trk_llr_pid_v_sel = awkward.fromiter([pidv[tid] if tid<len(pidv) else 9999. for pidv,tid in zip(trk_llr_pid_v,trk_id)])\n",
    "    df['trkpid'] = trk_llr_pid_v_sel\n",
    "    df['subcluster'] = df['shrsubclusters0'] + df['shrsubclusters1'] + df['shrsubclusters2']\n",
    "    df['trkfit'] = df['shr_tkfit_npointsvalid'] / df['shr_tkfit_npoints']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lee_bins = [0, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.6, 0.8]\n",
    "lee_scaling = [\n",
    "    6.3744101, 6.3744101, 5.6455402, 3.7305500, 1.5091400, 1.0742800, 0.7540929,\n",
    "    0.4763070, 0.1523270\n",
    "]\n",
    "#lee_scaling = [1,1,1,1,1,1,1,1,1]\n",
    "\n",
    "nc3[\"train_weight\"] = 1\n",
    "nue3[\"train_weight\"] = 0\n",
    "mc3[\"train_weight\"] = 1\n",
    "ext3[\"train_weight\"] = 1\n",
    "nc1[\"train_weight\"] = 1\n",
    "nue1[\"train_weight\"] = 0\n",
    "mc1[\"train_weight\"] = 1\n",
    "ext1[\"train_weight\"] = 1\n",
    "\n",
    "\n",
    "for i, lee_bin in enumerate(lee_bins):\n",
    "    \n",
    "    if i == 0:\n",
    "        continue\n",
    "        \n",
    "    nue3.loc[(nue3['nu_e'] > lee_bins[i-1]) & (nue3['nu_e'] < lee_bins[i]), 'train_weight'] = lee_scaling[i-1] * nue3['weightSpline']\n",
    "    nue1.loc[(nue1['nu_e'] > lee_bins[i-1]) & (nue1['nu_e'] < lee_bins[i]), 'train_weight'] = lee_scaling[i-1] * nue1['weightSpline']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nc3[\"is_signal\"] = nc3[\"category\"] == 11\n",
    "nue3[\"is_signal\"] = nue3[\"category\"] == 11\n",
    "mc3[\"is_signal\"] = mc3[\"category\"] == 11\n",
    "ext3[\"is_signal\"] = ext3[\"category\"] == 11\n",
    "nc1[\"is_signal\"] = nc1[\"category\"] == 11\n",
    "nue1[\"is_signal\"] = nue1[\"category\"] == 11\n",
    "mc1[\"is_signal\"] = mc1[\"category\"] == 11\n",
    "ext1[\"is_signal\"] = ext1[\"category\"] == 11"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# variables to be trained on\n",
    "TRAINVAR = ['tksh_angle',\"shr_tkfit_dedx_Y\",\"shr_tkfit_dedx_U\",\"shr_tkfit_dedx_V\",\"n_showers_contained\",\"shr_score\",\"tksh_distance\",\\\n",
    "            \"trkfit\",\"trkpid\",\"subcluster\",\"shrmoliereavg\",\"shrmoliererms\",\"trkshrhitdist2\",\"hits_ratio\",\\\n",
    "            'is_signal','train_weight','nu_e']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print (ls.pickle_path)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "train_mc3, test_mc3 = train_test_split(mc3, test_size=0.5, random_state=1990)\n",
    "train_ext3, test_ext3 = train_test_split(ext3, test_size=0.5, random_state=1990)\n",
    "train_nue3, test_nue3 = train_test_split(nue3, test_size=0.5, random_state=1990)\n",
    "train_nc3, test_nc3 = train_test_split(nc3, test_size=0.5, random_state=1990)\n",
    "train_mc1, test_mc1 = train_test_split(mc1, test_size=0.5, random_state=1990)\n",
    "train_ext1, test_ext1 = train_test_split(ext1, test_size=0.5, random_state=1990)\n",
    "train_nue1, test_nue1 = train_test_split(nue1, test_size=0.5, random_state=1990)\n",
    "train_nc1, test_nc1 = train_test_split(nc1, test_size=0.5, random_state=1990)\n",
    "\n",
    "# merge run1 and run3 samples\n",
    "train_mc = pd.concat([train_mc3,train_mc1])\n",
    "train_ext = pd.concat([train_ext3,train_ext1])\n",
    "train_nue = pd.concat([train_nue3,train_nue1])\n",
    "train_nc = pd.concat([train_nc3,train_nc1])\n",
    "test_mc = pd.concat([test_mc3,test_mc1])\n",
    "test_ext = pd.concat([test_ext3,test_ext1])\n",
    "test_nue = pd.concat([test_nue3,test_nue1])\n",
    "test_nc = pd.concat([test_nc3,test_nc1])\n",
    "mc = pd.concat([mc3,mc1])\n",
    "ext = pd.concat([ext3,ext1])\n",
    "nue = pd.concat([nue3,nue1])\n",
    "nc = pd.concat([nc3,nc1])\n",
    "\n",
    "samples = {\n",
    "    \"mc\": (train_mc, test_mc),\n",
    "    \"nue\": (train_nue, test_nue),\n",
    "    \"ext\": (train_ext, test_ext),\n",
    "    \"nc\": (train_nc, test_nc)\n",
    "} \n",
    "\n",
    "fig, ax = plt.subplots(1,1)\n",
    "\n",
    "my_booster = nue_booster.NueBooster(samples, TRAINVAR, random_state=1990)\n",
    "\n",
    "print (my_booster.variables)\n",
    "\n",
    "PRESEL = \"reco_e < 1.0 and selected == 1 and n_tracks_contained > 0\"# and (crtveto == 0) and (_closestNuCosmicDist > 20)\"\n",
    "PRESEL += \" and shr_energy_tot_cali > 0.07\"\n",
    "PRESEL += ' and n_showers_contained == 1'\n",
    "PRESEL += ' and hits_ratio>0.5'\n",
    "PRESEL += ' and tksh_distance < 6.0'\n",
    "PRESEL += ' and shr_tkfit_dedx_Y < 4.0'\n",
    "PRESEL += ' and tksh_angle > -0.9'\n",
    "PRESEL += ' and trkpid < 0.1'\n",
    "PRESEL += ' and shr_score < 0.30'\n",
    "PRESEL += ' and CosmicIP > 20.'\n",
    "\n",
    "my_booster.set_preselection(PRESEL)\n",
    "\n",
    "for label, bkg_query in zip(nue_booster.labels, nue_booster.bkg_queries):\n",
    "    \n",
    "    preds = my_booster.train_booster(ax, bkg_query)\n",
    "    \n",
    "    with open(ls.pickle_path+'booster_%s.pickle' % label, 'wb') as booster_file:\n",
    "        pickle.dump(preds, booster_file)\n",
    "\n",
    "    variables = my_booster.variables.copy()\n",
    "    print ('variables are : ',variables)\n",
    "    variables.remove(\"is_signal\")\n",
    "    variables.remove(\"nu_e\")\n",
    "    variables.remove(\"train_weight\")        \n",
    "        \n",
    "    mc_prediction = preds.predict(\n",
    "        xgb.DMatrix(mc[variables]),\n",
    "        ntree_limit=preds.best_iteration)\n",
    "    nue_prediction = preds.predict(\n",
    "        xgb.DMatrix(nue[variables]),\n",
    "        ntree_limit=preds.best_iteration)\n",
    "    ext_prediction = preds.predict(\n",
    "        xgb.DMatrix(ext[variables]),\n",
    "        ntree_limit=preds.best_iteration)\n",
    "    nc_prediction = preds.predict(\n",
    "        xgb.DMatrix(nc[variables]),\n",
    "        ntree_limit=preds.best_iteration)\n",
    "\n",
    "    nc[\"%s_score\" % label] = nc_prediction\n",
    "    mc[\"%s_score\" % label] = mc_prediction\n",
    "    nue[\"%s_score\" % label] = nue_prediction\n",
    "    ext[\"%s_score\" % label] = ext_prediction\n",
    "\n",
    "\n",
    "ax.set_ylim([0, 1.05])\n",
    "ax.set_xlim([0, 1.0])\n",
    "ax.set_xlabel('False Positive Rate')\n",
    "ax.set_ylabel('True Positive Rate')\n",
    "ax.set_title('ROC')\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "fig.savefig(ls.plots_path+\"roc_single.pdf\")"
   ]
  }
 ],
 "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.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}