timing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Timing Distributions of HNLs at SBN/DUNE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from hnl_flux import *\n",
    "from hnl_decays import *\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm, ticker\n",
    "from matplotlib.pylab import rc\n",
    "rc('font', **{'family': 'serif', 'serif': ['Computer Modern']})\n",
    "rc('text', usetex=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_timing_distribution(mN, mZp, Ualpha=1.0, flavor=1):\n",
    "\n",
    "    flux = FluxHNLFromProtonBrem(proton_energy=120000.0, mixing_flavor=flavor,\n",
    "                                mixing_angle=1.0, hnl_mass=mN, coupling_BL=1e-4,\n",
    "                                zprime_mass=mZp, n_samples=100000)\n",
    "    \n",
    "    flux.simulate()\n",
    "    flux.decay_to_hnl()\n",
    "    flux.propagate(Ualpha)\n",
    "\n",
    "    return flux.hnl_timing, DUNE_EFFICIENCY * DUNE_POT_TOTAL * flux.decay_axion_weight\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "test_mN = 200.0  # MeV\n",
    "test_mZp = 3*test_mN\n",
    "test_Ualpha = 1e-4\n",
    "\n",
    "\n",
    "timing, weights = get_timing_distribution(test_mN, test_mZp, test_Ualpha, flavor=2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4913\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAG5CAYAAAAu8Fa/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAh9klEQVR4nO3dz27iSMPv8V8etZTVJG5mN3pz1HLugMAVDNFZzBa6ryCwPrOAZvWoVwyszhb6Crrtba/iuYIkXr1b6hkp0uzGMYx0JDbjs2jhNwRIzJ9UMPl+JDSxKdtlTA2/LpftgyRJEgEAAFjyr5euAAAAeF0IHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsevPSFXjon3/+0Z9//qkffvhBBwcHL10dAACQQZIk+vvvv/XTTz/pX/96vG9j58LHn3/+qZOTk5euBgAAWMPt7a3+67/+69EyOxc+fvjhB0nfK390dPTCtQEAAFmMx2OdnJykv+OP2bnwMT3VcnR0RPgAACBnsgyZYMApAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsevPSFbDt3cdvT5b547dfLNQEAIDXiZ4PAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYNVKdzg1xsj3fbmuK2OM6vW6HMdZWjYIAhUKBRljVK1W5bruNuoMAABybKXwUavVdHNzI+l7uLi4uJDneQvL+r6vZrOZTjcaDfX7/Q2qCgAA9kHm0y7GmJlp13UVBMHS8l++fFm/VgAAYG9lDh/TUyj3FQoFhWG4sHyhUNDZ2Vl6+uX8/HyzmgIAgL2QOXzEcbxwfhRFC+dPT8ecnp7K8zxVq9WF5SaTicbj8cwLAADsr5XGfCyyLJQEQaButytjjBqNhiQtHPPR6XT06dOnTasBAAByInPPh+M4c70cURQtvNrFGKOrqytVKhXV63UNh0N9/fp1btyIJLXbbY1Go/R1e3u7+l4AAIDcyBw+KpXKwvmlUmluXhiGKpfL6bTrumq32wt7SQ4PD3V0dDTzAgAA+ytz+Hh4jw5jjEqlUtrzEYZh2rNRLBZ1dXU1U/6vv/5SsVjcsLoAACDvVhrz4XmeWq2WyuWyrq6uZu7x0el0VC6X1Ww25bquzs/P1ev10nAyHfcBAABet4MkSZKXrsR94/FYx8fHGo1Gz3IK5t3Hb0+W+eO3X7a+XQAA9tkqv9882wUAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVW9WKWyMke/7cl1XxhjV63U5jrOwrO/7qlQqkrS0DAAAeH1WCh+1Wk03NzeSvgeRi4sLeZ63tOxD3W5XzWZzjWoCAIB9kfm0izFmZtp1XQVBsLBsHMfyPE9JkqQvggcAAJBWCB9BEKhQKMzMKxQKCsNwYflqtZr+7fv+zDQAAHi9Mp92ieN44fwoiubm3R/jEcexoiiS67oLl59MJppMJun0eDzOWiUAAJBDG1/tsiyUTLVaLb1//37p+51OR8fHx+nr5ORk0yoBAIAdljl8OI4z18sRRdGjV7LEcawgCB4t0263NRqN0tft7W3WKgEAgBzKHD6ml80+VCqVli5zfX395GW2h4eHOjo6mnkBAID9lTl8PByzYYxRqVRKw0UYhnNXxIRhODdIFQAAvG4rjfnwPE+tVku+76vf78/c46PT6cj3/blllg00BQAAr9NBkiTJS1fivvF4rOPjY41Go2c5BfPu47cny/zx2y9b3y4AAPtsld9vnu0CAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwKo3qxQ2xsj3fbmuK2OM6vW6HMdZWj4IAhlj5LquJKlSqWxUWQAAkH8rhY9araabmxtJ34PIxcWFPM9bWDYIAnmep36/L2OMzs/PNRwON68xAADItczhwxgzM+26roIgWFq+0WikQcV1XV1eXq5ZRQAAsE8yj/kIgkCFQmFmXqFQUBiGc2WNMYqiSI7jKAxDxXGcnnoBAACvW+bwEcfxwvlRFM3NC8NQhUIhHR8yGAzk+/7C5SeTicbj8cwLAADsr5XGfCyyKJREUSRjjCqVihzHUb1e19u3b5UkyVzZTqejT58+bVoNAACQE5l7PhzHmevlmJ5aech1XTmOk743/e+iUzTtdluj0Sh93d7eZq89AADInczhY9llsqVSaW7eKuM7Dg8PdXR0NPMCAAD7K3P4eBgojDEqlUozvRrTK2Jc11WpVEpPyUzv9VEsFrdTawAAkFsrjfnwPE+tVkvlcllXV1cz9/jodDoql8tqNpszZc/OznRzc8OltgAAQJJ0kCwaBfqCxuOxjo+PNRqNnuUUzLuP354s88dvv2x9uwAA7LNVfr95tgsAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAqjerFDbGyPd9ua4rY4zq9bocx1lYNgxDSVKxWJQxRnEcq1gsblxhAACQbyuFj1qtppubG0nfg8jFxYU8z1tYtt/vazAYSJIqlcrScgAA4HXJHD6MMTPTrusqCIKl5c/OznR3dydJS3tHAADA65M5fARBoEKhMDOvUCgoDMOlp1OyhI7JZKLJZJJOj8fjrFUCAAA5lHnAaRzHC+dHUbS0vO/78n1frVZrrudkqtPp6Pj4OH2dnJxkrRIAAMihlcZ8LLIslNwfjOq6rs7PzzUcDufKtdtt/frrr+n0eDwmgAAAsMcy93w4jjPXyxFF0dJTK/d7OqZXxyzq/Tg8PNTR0dHMCwAA7K/M4aNSqSycXyqV5uaFYaiff/55bv7DMSMAAOD1yRw+XNedmTbGqFQqpT0fYRimPRuu66rb7aZlgyBQtVrlqhcAALDamA/P89RqtVQul3V1dTVz745Op6NyuaxmsynHcVQqldTr9eQ4jobDIff5AAAAkqSDJEmSl67EfePxWMfHxxqNRs8y/uPdx29Plvnjt1+2vl0AAPbZKr/fPNsFAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFUrhQ9jjHq9nnzfV6/XUxzHmZZrtVqZywIAgP22Uvio1WpqNpuqVquqVqu6uLh4cpkwDNXr9dauIAAA2C+Zw4cxZmbadV0FQZBpOdd1V68ZAADYS5nDRxAEKhQKM/MKhYLCMFy6jO/7qlarj653MploPB7PvAAAwP7KHD6WjdmIomhpecdxnlxvp9PR8fFx+jo5OclaJQAAkEMbX+2yLJR8/fpVlUrlyeXb7bZGo1H6ur293bRKAABgh73JWtBxnLlejiiKFvZuBEGg9+/fZ1rv4eGhDg8Ps1YDAADkXObwUalU1O/35+aXSqWF5b9+/Zr+bYxRp9PRhw8fVCwW16gmAADYF5nDx8MrVowxKpVKac9HGIZyHEeu686dbmk0Gmo0Glz1AgAAVhvz4XmeWq2WfN9Xv9+X53npe51OR77vz5SP4zi9x0e32330yhgAAPA6HCRJkrx0Je4bj8c6Pj7WaDTS0dHR1tf/7uO3J8v88dsvW98uAAD7bJXfb57tAgAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwKo3L12BXfTu47cny/zx2y8WagIAwP6h5wMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVSvd58MYI9/35bqujDGq1+tyHGdh2SAIJElxHOvq6kofPnxQsVjcuMIAACDfVgoftVpNNzc3kr4HkYuLC3met7Ts77//rkqloiiKVKvVNBwON68xAADItcynXYwxM9Ou66a9G4t4njfT07GshwQAALwumcNHEAQqFAoz8wqFgsIwXFi+Uqmkf3uep0ajsWYVAQDAPsl82iWO44XzoyhaukwYhvry5YvOz89Vr9cXlplMJppMJun0eDzOWiUAAJBDG1/tsiyUSFKxWFS73dZwOJTv+wvLdDodHR8fp6+Tk5NNqwQAAHZY5vDhOM5cL0cURU+O5XAcR7VaTbVabWFQabfbGo1G6ev29jZrlQAAQA5lDh/3x3DcVyqV5uYFQaC3b9+m067rSpoftCpJh4eHOjo6mnkBAID9lTl8TAPElDFGpVIp7fkIwzANF4VCYSashGEox3G4zwcAAFjtPh+e56nVaqlcLuvq6mrmHh+dTkflclnNZlPFYlEfPnzQYDCQJF1eXqb3BwEAAK/bQZIkyUtX4r7xeKzj42ONRqNnOQXz7uO3raznj99+2cp6AADYB6v8fvNsFwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFVvXroCefXu47dM5f747ZdnrgkAAPlCzwcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArFrpUltjjHzfl+u6MsaoXq/LcZyFZcMwVBAEkqSrqyt9/vx5aVkAAPB6rBQ+arWabm5uJH0PIhcXF/I8b2HZIAjUbDYlSb1eTz///HO6LAAAeL0yn3YxxsxMu66b9mw8FIahOp1OOl2tVhWG4dw6AADA65M5fARBoEKhMDOvUCgoDMO5ssViUZ8/f06n4zhOywMAgNct82mXaYB4KIqihfOr1Wr695cvX1SpVBaO+ZhMJppMJun0eDzOWiUAAJBDGz/bZVkouf++7/tLx3t0Oh19+vRp02rkWpbnxPCMGADAvsh82sVxnLlejiiKnryCpdVq6fLycmm5drut0WiUvm5vb7NWCQAA5FDm8FGpVBbOL5VKS5fp9XpqtVpyXVdxHC/sJTk8PNTR0dHMCwAA7K/M4cN13ZlpY4xKpVLao/Hwahbf91UsFtPg8fXrV+7zAQAAVhvz4XmeWq2WyuWyrq6uZu7x0el0VC6X1Ww2ZYxRrVabWdZxHNXr9e3UGgAA5NZK4cN1XXW7XUmzV7NImgkirusqSZItVA8AAOwbnu0CAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsIrwAQAArCJ8AAAAq968dAX23buP3166CgAA7BR6PgAAgFWEDwAAYBXhAwAAWEX4AAAAVhE+AACAVYQPAABgFeEDAABYRfgAAABWrRQ+jDHq9XryfV+9Xk9xHD9aPgxDnZ2dbVI/AACwZ1a6w2mtVtPNzY2k70Hk4uJCnuctLOv7vlzXVRiGm9cSAADsjczhwxgzM+26roIgWFq+Wq2uXysAALC3Mp92CYJAhUJhZl6hUKBnAwAArCRzz8ey8R1RFG1Ugclkoslkkk6Px+ON1gcAAHbbxle7PDXo9CmdTkfHx8fp6+TkZNMqAQCAHZY5fDiOM9fLEUWRHMfZqALtdluj0Sh93d7ebrQ+AACw2zKHj0qlsnB+qVTaqAKHh4c6OjqaeQEAgP2VOXy4rjszbYxRqVRKez7CMJy7ImZq01MzAABgf6w05sPzPLVaLfm+r36/P3OPj06nI9/30+kgCNRqtRa+BwAAXq+DJEmSl67EfePxWMfHxxqNRs9yCubdx29bX6cNf/z2y0tXAQCApVb5/ebZLgAAwCrCBwAAsIrwAQAArCJ8AAAAqwgfAADAKsIHAACwivABAACsInwAAACrCB8AAMAqwgcAALDqzUtXANlkuS08t2AHAOQBPR8AAMAqwgcAALCK8AEAAKxizMceYVwIACAP6PkAAABWET4AAIBVhA8AAGAV4QMAAFhF+AAAAFYRPgAAgFWEDwAAYBXhAwAAWMVNxrAQNywDADwXej4AAIBVhA8AAGAV4QMAAFhF+AAAAFYx4PSVyTKQdJvrYlAqAOAhej4AAIBV9HwgF7L22NjsaaHnBwDWs1L4MMbI9325ritjjOr1uhzH2bgsAAB4PVYKH7VaTTc3N5K+h4uLiwt5nrdxWbxu2xyHgtdtF3vIAMzLHD6MMTPTrusqCIKNywLIhtM8dvF5I6/y8N3NHD6CIFChUJiZVygUFIahisXi2mWx32z3auSh0S1C708+bev7Ro8NXpvM4SOO44XzoyjaqOxkMtFkMkmnR6ORJGk8Hmet2kr+mfy/Z1kv8uN//Z+nT//996f//WSZLN+lrN/jbX0vn6vd5EXWzzHL52TzmGyz3sA2/9+0iuk6kyR5suzGV7ssCxpZy3Y6HX369Glu/snJyQa1AjZz/H93az27ur28svk5bXNbHF9sy3N+l/7++28dHx8/WiZz+HAcZ67nIoqihVewrFK23W7r119/Taf/+ecfRVGkH3/8UQcHB1mrl8l4PNbJyYlub291dHS01XXvgn3fP2n/95H9y79930f2L/+eax+TJNHff/+tn3766cmymcNHpVJRv9+fm18qlTYqe3h4qMPDw5l5z31J7tHR0d5+qaT93z9p//eR/cu/fd9H9i//nmMfn+rxmMp8h1PXdWemjTEqlUppUAjDML3K5amyAADg9VppzIfneWq1WiqXy7q6upq5b0en01G5XFaz2XyyLAAAeL1WCh+u66rb7UqSqtXqzHsPw8VjZV/K4eGh/v3vf8+d5tkX+75/0v7vI/uXf/u+j+xf/u3CPh4kWa6JAQAA2BKeagsAAKwifAAAAKsIHwAAwCrCBwAAsGrj26u/BGOMut2uTk9PJX2/KVm9Xt/Kcuuue5vWrUMQBLq8vFQcxzLGqFarzS13dnamdrutSqUiSRoMBpKUXiJty7r7mKX+eT2GvV5Pw+FQZ2dncw9mlP7nqrFdOYbTbQ+Hw/TKtqfkpQ1Orbp/eWqD022vsn95aX/3rbKPeWuDWb5vi+xEO0xyZjgcJo7jJHd3d+m8ZrOZdLvdjZdbd93btG4dLi8vZ8rc3d0lrusm9Xp9ppykmdfD923Y5HN+qv55Pob1en1u/6avYrGYlnvpYzgcDpN6vZ7U6/XEcZyk2WxmXi4vbXCd/ctLG1x3/5IkH+1vWo919jEvbTBJsn/fHtqVdpi78FGv1+e+SHd3d8lTOSrLcuuue5vWrUO1Wp2b1+/3E0nJcDicWX+/30/6/f7MfJs2+Zyfqn+ej+Gy/2l0u92dO4ZTxWJxpf+x56EN3rfK/uWpDU6tsn9Jko/299Cq39FFdrENZv2+PbQr7TB3Yz6+fv2adgNNTW/bHgTBRsutu+5tWrcOvu+r1WrNzJs+S+f+cqenp6rX66rX63O3wbdlk8/5qfrn+Rien5/PzQvDUI7jzOzrLhzDdeSlDa4rT21wXXlof5vIUxvM+n17aFfaYa7CRxzHiuN44cF2HEdhGK693Lrr3qZN6lCtVue+LI9tJwgCK/u0aNubfs7L6r8Px/Chfr+/8DzrSx7DdeSlDW4iL21wU7vc/jaVpza4yvdtapfaYa7Cx/TBdYsUCgX99ddfay+37rq3aZM6eJ4310Cur68lKR0UJUmXl5cKgiBNyOfn51Ybz6af82P1z/sxfKjVai0cJPfSx3AdeWmDm8hLG9zErre/bdvlNpj1+3bfLrXDXF7tskwcx8+23Lrr3qZV69DtdtXtdmcSbL/fT6eLxaIajYZqtZqGw+E2q7q2p/Zxk/rn6RgaYxTH8cInQe/6MVxHXtrgqvLYBh+T9/a3ijy2wUXft1XYbIe56vlY9CWYiqJoo+XWXfc2bbMOtVpNlUpl7tKvh1/KYrEoY4y187Gb7uNj9d+nY9jtdheef5Ze/hiuIy9tcJt2tQ1uYtfb3zblrQ0u+77dt0vtMFfhY3rd9aLktSyhZl1u3XVv07bqMBgMVCgU1O/3Z+a3Wq25rsHpNh/ratumTfbxqfrv0zH8+vWrisXi3PxdOIbryEsb3JZdboPrykP726Y8tcFl37eHdqkd5ip8OI4jx3GWpq9lKTXLcuuue5u2UQff9xXH8cyXcPol6vV66TnBqem2bI3Y3mQfn6r/vhzDIAiWDvjahWO4jry0wW3Y9Ta4rjy0v23JUxt87Pv20C61w1yFD0l6//793Hm1adpcNsgm63LrrnubNqlDGIaKomim2206Ilv63o34cIDStLvU1v5J6+9jlvrn/RjeL7vIrhzDdeSlDW4iL21wHXlpf9uQlzb41PdtkZ1ph1u7Y4glw+EwcV13Zl6z2Uz6/X46fXd3l1QqleTm5mal5bKUeW6b7F+1Wk08z5t51ev1tNzl5WXied7MelzXnZlnw7r7mKX+eT6G98sua5q7cgynlt1RMc9t8L5V9y8vbXBqlf3LS/t7aJV9nMpDG8zyfdvldniQJEmynRhjTxiG+vLli8rlcprGHj5b4OzsTJ7nzaS0p5bLWua5rbN/b9++XdrVdv8QT58FMF1Po9F4kX+RrHsMs9Q/r8dwajAYqNvtLh05/9LHMI5jdTodxXGswWAgx3H0/v17nZ6epvuY5za47v7lpQ1ucvzy0v422Udp99uglO37tsvtMJfhAwAA5FfuxnwAAIB8I3wAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwBs0XM8On7XHzoHrIrwAWyo0WjI9/2XrgZ2QBzHarVaW19vGIYv+rh2YNsIH8AGjDEaDAYr/eD4vq9arfaMtdrt7e+K5/gcLi4u1O12t7pOSapWq/I8jx4Q7A3CB7CBfr+vbrcrY0zm3o8vX76oUCg8c812d/u7Ytufw2AwSB9J/hy63a4ajcazrBuwjfABrMkYkz6oynGczL0fQRDo/Pz8mWu3u9vfFdv+HPr9/tyj1rfJcRy5rsvpF+wFwgewpvs/Nu12W8aYJ38Y4jhWHMcv8iThXdj+rtj25xAEgUql0lbW9ZharaZ+v//s2wGeG0+1BdYwfVT3/UdMv337Vq7r6ubmZq789BHcYRjq+vp6JrQ8Vzf9Lm1/VzzX59BoNHR+fq5qtZrOC8NQFxcXMsbo999/T8drXF1dSdLc2BDf9xVFkQqFgqIo0s3NjWq12lxAevv2re7u7tauK7ATEgArazabc/O63W4iKbm8vFy6XL1eT6rV6nNW7VEvvf1dse3PwXXdZDgcLnxPUlKv15O7u7t0nuM4ied56fRwOEzq9frMct1ud+F36bFtAXnBaRdgRXEc68cff5ybn2XsxyrjDKb/ms76yjLg9ant93o9HRwcqNFoKAzDdJmzszOdnp5qMBjMlB8MBjo4OJjb516vp7dv36rX6y3dRq1Wm1ufLdse7xFFkVzXXfie4zjpa6pUKqU9INL3XpKHV7Lc70W5z3VdrnpB/r10+gHyptlszvwr9uF7kpKbm5u59+7u7hJJL/av1izbHw6HiaS5/atUKkm3211Yvl6vJ47jzL23qPxj27DlOY7DY/8rfdjLkSRJUq1WZ3rP7u7uEsdxEtd1k2az+WjvWbVaTfr9/uaVBl4QPR/ACqZ3r1w2PqDdbkvSwt6P6+vr9IqFl5Bl+2EYynXduf0LgmDh4MwwDNXtdlUoFGZ6XowxSwdzLtvGY3U6OzvL/Jr22Cyz7eOwjTuaOo6j//znP6pUKvJ9X+fn5zo9PX2Wu6UCu+DNS1cAyJPBYJAGjEUcx1Gz2VSv11MYhioWi+l7l5eXM1dExHH86A9wo9FYqXu90Wgs7arPuv3Ly8u50DD9Mb+/L1NRFMlxHDUaDXU6nXT7QRAsvex00TYeUywWFw7iXdeqx+Ep2xiwO/2uTK9kMcakn+nDgalxHL9YgAW2hfABZBTHsf76668nf2za7bZ6vZ5arZYuLy/T+WEYzowzeHi1zEPbvqQyy/aDIFCxWJwZi3F5ebkweNxXr9fVarVkjHnyhzEIgme5BXlWT30Ovu/P7MdjgW7KcZyNQowxRlEUpaHMdV31+/2FNxWbXhED5BnhA8io0+noxx9/zDSws1gsKgiCud6P6Q/aYDB41htSLfPY9uM4ljFGl5eXMwEiS0+F4ziqVqvqdruPhqbpNl76PiPLPofBYKDLy0t5npdeKpslfEwHgT4V0qam9xm5r9vtznwuxpiFg2KzBDxg1xE+gIwGg8HK5+A7nY48z5P0fRzI9O9SqWT9/hpPbX/ZWIggCPT58+e59T0cBzK9OqfVai39cXxsG7YCybLPIY5jNRqN9B4axWJRv//+e6Z1fvjwQdfX1zPhIwxD9ft9xXGsTqejKIr0/v17dTodXV9f6/r6WtL/9HDVajX1er2Z+jzsmZqecnlN92bBfuImYwAkff9RDsNw5lRRHMfpTa0e/uAt6j05PT2V67ryPG/hD+SibUjfL7997BSUDb7vq9PprDW+xBgzE2qey/RGZC/RawZsE1e7AJA0PxZC+t4jsexf2ot6gabjPpb9y3zRNlqt1oufhpG+nzp6OJYi63NUpqddnvvqlE6no/fv3z/rNgAbCB/AKxfHsXq9noIg0HA4TK+wmfYETP+eml6J0e/3536c379/v3CMxMNtDAYD9Xq99OZoWcdKPKdKpZIOtvV9X4PBYKWxFZ8/f04/r+cQhqEqlQqnXLAXOO0CAFvS6/XSELNNcRzr4uLi2U/rALbQ8wEAW9JsNhUEwdZPvwwGg4WDfoG8oucDAABYRc8HAACwivABAACsInwAAACrCB8AAMAqwgcAALCK8AEAAKwifAAAAKsIHwAAwCrCBwAAsOr/Aw9kHTijESwQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "tof_light = (DUNE_DIST / (1e-2 * C_LIGHT))\n",
    "\n",
    "timing_bins = np.logspace(-3, 3, 100)*1e-9\n",
    "timing_bins = np.linspace(0.0, 2.0, 50)\n",
    "\n",
    "print(timing.shape[0])\n",
    "\n",
    "plt.hist(1e9*(timing - tof_light), weights=weights, bins=timing_bins)\n",
    "#plt.xscale('log')\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r\"$\\Delta t = t_{HNL} - t_c$ (ns)\", fontsize=14)\n",
    "plt.xticks(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAh8AAAG5CAYAAAAu8Fa/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAySklEQVR4nO3db2wbaYLn959aI1OSNXaJsh1Eak62Sxrvi+kBDEpq4IIDjMFSiwVyeTNHWjvI3b6TBGRysdHeFaFXi35xYMi5AXpy6CCk3s0esCOT6EGAvReJqmMMcq/GFmFkb4BLt1RerNbuhS1TZbWsP1bbzAstOaZEUiWJLP77fgABJuth1VOiyfrpqedPVz6fzwsAAMAj7zW6AgAAoLMQPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BThAwAAeIrwAQAAPEX4AAAAnvpOoytw1Nu3b/X06VN997vfVVdXV6OrAwAAXMjn8/rmm280PDys996r3rbRdOHj6dOnCgQCja4GAAA4g/X1db3//vtVyzRd+Pjud78r6bDyly5danBtAACAG1tbWwoEAsXreDVNFz4Kt1ouXbpE+AAAoMW46TJBh1MAAOApwgcAAPBU0912AQB0nrdv3+r169eNrgZOcOHChRNHsrhB+AAANNTr16/1+PFjvX37ttFVwQnee+89ffDBB7pw4cK59kP4AAA0TD6f19dff63u7m4FAoGa/FWN+ijMw/X111/re9/73rnm4iJ8AAAa5ttvv9XOzo6Gh4fV39/f6OrgBFevXtXTp0/17bffqqen58z7IWICABrmzZs3knTuZnx4o/A+Fd63syJ8AAAajuU0WkOt3iduuwAAms4TZ1ebr7wb/TJ48YJGjD7PjtfpCB8AgKbyxNlV6Oe/0e7B+Zr2T6Ovp1vW3ZuuAkg2m1UymVQqldL8/Lymp6dlmqbu3bunaDSqiYkJRSIRTUxMaGlpSYlEQslkUrOzs5Ik27YVj8dlWVaxfGF/s7OzmpqaUjgcPnZc27aVyWRkGIYkyTRN2bZd3G8r6crn8/lGV+JdW1tbunz5sl6+fMn06gDQ5vb29vT48WN98MEH6u3tlST95ycv9S/+/X/Sp9M3NHZtoO51WH22rTtLj/Q3/+af68ORy65eY9u2RkdHtbm5WQwDkjQ+Pq7p6WnNz89LkhzHUTQaVSqVOlY2kUgUy1Xa37vGx8e1srJS8voXL14oHo+f7oTPodz7VXCa6zctHwCApjR2bcB1GPCa3+93XTYSiSiXy2lmZkbpdLr4/Lsh46T92bZ97Ln5+XklEgnX9WgmdDgFAKDOFhcXZVmWMpnMmV5fuMWSSqVKnm/FWy4SLR8AmoCbzoV0CEQzSqVSJS0YjuOULWcYhuLxuGZmZsr253BjcXFRkUhEc3NzCoVCikajCoVCZ9pXoxE+ADSU286Fp+kQCHhldna2JHwkk8mqZdPptObm5qqWqyQcDmttbU2WZWl5eVlTU1NKp9NnDjONRPgA0FCbr15r9+BN1c6FhQ6Bm69eEz7Q0pLJpEZHRzU3N+eqfCaTUTgcluM4MgxDpmlqdnZWs7OzSqVSisViLRk+6PMBoCkUOheW+/FixAPgBdM0FY/HFYlEXJUvdDS1bVvZbLZk261btyre5ml2hA8AAE4pl8ududz8/LwMwygJDuXKOY6j5eXl4uNoNFqy3bKslmz1kLjtAgBoUqvPtpvyOIVJxiQpFouVTDJm27aWlpZkmqZM01Q0GlUul9PCwkJJUFhcXNTDhw/L7m90dFRra2vKZDIyTbP4mkgkokQiUexjsra25ukcH7XEJGMAGqowoVS1CZ7clEFrKjdpVbPPcNrJmGQMQN24XVeD4a+ohxGjT9bdm6zt0sYIHwBKnOavTv5aRL2MGH38v2pjhA8AJdwMfZUY/grg7AgfAMpq5nU1ALQ2htoCAABPET4AAICnCB8AAMBT9PkAADQfZ13aeeHd8fqHJCPg3fE6HOEDANBcnHXps4+kgx3vjtnTL/30twQQjxA+AADNZefFYfD48aJ05Xr9j7fxpfT5zOFxXYSPbDarpaUlJRIJJZNJzc7OSjpc/C0ej8uyLEWj0eLz9TQ+Pn5s6vZWQPgAADSnK9el4RuNrsUxwWBQpmnKcRzNzc3p1q1bxeXuk8mkEolETYKH4zjFdVwqicfjmpiYOPexvEaHUwAAziASiSgcDmtmZqbk+ZMCgxu2bevevXsnlguFQjU5ntcIHwAAnNHi4qIsy1Imk6npflt1tVq3CB8AAJyRYRiKx+PHWj8KLMtSIpFQJpNRNBotPjc+Pq5EIiFJymQyGh0dlWVZxe0PHz7U8vKyUqmUbNuWZVkaHR1VKpVSKpXS+Pi4stmsxsfHlUqlSvabSqVkWZZSqZQikUhJfbLZbLE+iURCtm1LOrzFU3jdu3WtF/p8AABwDrOzs0qn05qbm1MymSw+b9u2otGoVlZWJEm5XE6JRELz8/Oanp4ulguHw1paWio+DoVCCoVCGh0dLfYdMU1ToVBIKysrSiaT8vv9CgaDJfspvG55eVnpdFqSlE6nlc1mFQwGi/VZXl4uvmZ8fFxffPGFUqmUgsGgQqFQsa71RPgAAOCcksmkRkdHNTc3V/Kc3+8vtmhI0oMHD858DMMwNDQ0JEkVR7cMDQ0VyxReUwgSyWRSwWCwpLxpmrp3757C4bDGx8dlmqamp6frPlKH2y4AAJyTaZqKx+PHbnMUWhNCoVCxheS0HMcpOU49+P1+bW5uanFxUS9evDh2HrVGywcAoDltfNm0xyl3W2J+fl5LS0vFsDA9PX2sL4hlWcURKmtrayXPv3sL5ehrCi0d5Y77bjipplx9stmsFhcXFYvFNDc3p2AwqGAwSPgAAHSY/qHDGUc/L9+Jsy56+g+P60I2m1U0GlUulzs2wdfi4qIePnwo6bDVIx6PKxqNanJyUpKKfSpu3bqlaDRavCUTCoWKt0VM09Tc3Jzi8bhSqZRCoZAsy5JlWcpms8X+H4XJzvx+v8LhsBzHKfYdCYVCsm1b2WxWyWRSpmkW65NIJGSaph48eKB0Ol28nWNZlvx+v3K5XMUgVCtd+Xw+X9cjnNLW1pYuX76sly9f6tKlS42uDtB2nji72nz1uuL21WfburP0SH/zb/65Phy5XLHcf37yUv/i3/+nE8udxM1+anUsNJ+9vT09fvxYH3zwgXp7e3+/gbVdmlLF90unu37T8gF0kCfOrkI//412D95ULdfX063Bixc8qhVQhhEgDLQxwgfQQTZfvdbuwRt9On1DY9cGKpYbvHhBI0afhzUD0EkIH0AHGrs2wO0LAA3DUFsAAOApwgcAoOGabOwDKqjV+8RtFwBAw/T09Kirq0vPnz/X1atX1dXV1egqoYJ8Pq/nz5+rq6tLPT0959oX4QMA0DDd3d16//339Q//8A/6u7/7u0ZXByfo6urS+++/r+7u7nPth/ABAGiogYEBff/739fBwUGjq4IT9PT0nDt4SIQPAEAT6O7urslFDa2BDqcAAMBThA8AAOApwgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BThAwAAeIrwAQAAPMWqtgCOc9alnRdVi/RubGtYGx5VCEA7IXwAKOWsS599JB3sVC02Jsny+bS+PSnpsidVA9AeCB8ASu28OAweP16UrlyvWGz9q0cK3L+t7r2ch5UD0A5OFT5s21Ymk5FpmrJtW7OzszIM48TXRaNRLSwsuCoLoElcuS4N36i4ef/5tnd1AdBWTtXhNBKJaH5+XuFwWOFwWDMzMye+JpvNKpFInLmCAACgvbgOH7Ztlzw2TVOWZbl6nWmap68ZAABoS67Dh2VZ8vv9Jc/5/X5ls9mKr8lkMgqHw2evHQAAaDuu+3w4jlP2+VyufGczx3Fc9fHY39/X/v5+8fHW1pbbKgEAgBZ07tEulULJvXv3NDs7e+LrY7GYPvnkk/NWA4CkJ86uNl+9rrh99RmdRAE0nuvwYRjGsVaOXC5XtnXDsizdunXL1X4XFhb08ccfFx9vbW0pEAi4rRaAf/LE2VXo57/R7sGbquX6ero1ePGCR7UCgONch49QKKRkMnns+YmJibLl7927V/y3bduKxWKanp5WMBgsKefz+eTz+dxWA0AFm69ea/fgjT6dvqGxawMVyw1evKARo69mx/U5q9LTysdT/5Bk8AcFgN9zHT6OjlixbVsTExPFlo9sNivDMGSapkKhUEnZubk5zc3NMeoF8MDYtQF9OFL/GUff9Pq1k/cpcP+2dL9KwZ5+6ae/JYAAKDpVn490Oq1oNKrJyUk9ePBA6XS6uC0Wi2lyclLz8/PF5xzHUSqVkiTF43HNzc0da/kA0JoOBkYU2v+ZfvmTUY1drdDysfGl9PnM4ayphA8A/+RU4cM0TcXjcUk6NoT23SBSYBiG5ufnSwIJgPbxVFe0d+WH0jBruwBw71QznAIAAJwX4QMAAHiKVW0BtA9n/bB/yUkYgQM0FOEDQHtw1qXPPpIOdk4uywgcoKEIHwDaw86Lw+Dx40XpyvXK5RiBAzQc4QNAww1rQ70bfyt1lR+y27uxrWFtuNvZlevS8I3aVa6Kk6azl2o/qRvQDggfABqqZ/uJLN9fqP/X+xXLjEmyfD6tb09Kao5hvaeZzt66e5MAAryD8AGgobr3curv2tf6j36hwPdvlC2z/tUjBe7fVvde+VW0G8HNdParz7Z1Z+mRNl+9JnwA7yB8AGgK+8ZYxdsl+8+bdzVer6azB9oJ83wAAABP0fIB4FxWn1Vulejd2NaYdDjCpAKfs1r7SgFoaoQPoNOcNBFXlaDwrsGLF9TX0607S48qlhnWhiyfT/2fz1QsE5C0k/fpTa/f1XEBtD7CB9BJ3E7E1dN/OAtoFSNGn6y7N6sONV19tq3QUvWVb1efb+vP/npNqYGRE6sPoD0QPoBO4nYiLpfTj48YfSeO4jhp5du9/Es91csTj1Vrq8+3tZevfFzm5wDqh/ABdCIPJ+JqNs+293VN0u1fPdLvqoQP5ucA6ofwAaCjbO0e6JqkP//jP9TV6x+VLcP8HEB9ET4AdKSAv09jzM8BNATzfAAAAE8RPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BSjXQC0DJ+zKj0tP1Oq22nhATQe4QNA03vT69dO3qfA/dvS/SoFXUwLD6DxCB8Amt7BwIhC+9XXiJHkelp4AI1F+ADQEk5aI0aSnji72nxSfZ2Y57ldjdW6cgBOhfABoC08cXYV+vlvtHvwpmq5H3Q91o980qW+Ho9qBuAowgeAtrD56rV2D97o0+kbGrtW+dZM78Zl6dfStQGfh7UD8C7CB4C2MnZtQB9WW7Olq0qfkTpZfbZddfvgxQssYIeOQvgAgDoZvHhBfT3durP0qGq5vp5uWXdvEkDQMQgfAFAnI0afrLs3tfnqdcUyq8+2dWfpkTZfvSZ8oGMQPgCgjkaMPkIFcATTqwMAAE/R8gG0Cmdd2nlRcXPvxraGteFhhQDgbAgfQCtw1qXPPpIOdioWGZNk+Xxa356UVGW0BwA0GOEDaAU7Lw6Dx48XpSvXyxZZ/+qRAvdvq3sv53HlAOB0CB9AK7lyXRq+UXbT/vPqc0kAQLOgwykAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BThAwAAeIrwAQAAPEX4AAAAniJ8AAAATxE+AACApwgfAADAU4QPAADgKcIHAADwFOEDAAB46juNrgAA91afb2sv/7Lstue5XY15XJ+WtvFlxU29G9sa1oaHlQE6C+EDaAHPtvd1TdLtXz3S7yqEjx90PdaPfNKlvh5vK9dq+oeknn7p85mKRcYkWT6f1rcnJV32rGpApyB8AC1ga/dA1yT9+R//oa5e/6hsmd6Ny9KvpWsDPm8r12qMgPTT30o7LyoWWf/qkQL3b6t7L+dhxYDOQfgAWkjA36exkQp/iXcNeFuZVmYEDn8q2H++7WFlgM5Dh1MAAOApwgcAAPAU4QMAAHiK8AEAADx1qg6ntm0rk8nINE3Ztq3Z2VkZhlG2rGVZkiTHcfTgwQNNT08rGAyeu8IAAKC1nSp8RCIRraysSDoMIjMzM0qn0xXLfvHFFwqFQsrlcopEIlpbWzt/jQEAQEtzfdvFtu2Sx6ZpFls3ykmn0yUtHZVaSAAAQGdxHT4sy5Lf7y95zu/3K5vNli0fCoWK/06n05qbmztjFQEAQDtxfdvFcZyyz+dylWcAzGazWlpa0tTUlGZnZ8uW2d/f1/7+fvHx1taW2yoBAIAWdO7RLpVCiSQFg0EtLCxobW1NmUymbJlYLKbLly8XfwKByrMOAgCA1uc6fBiGcayVI5fLndiXwzAMRSIRRSKRskFlYWFBL1++LP6sr6+7rRIAAGhBrsPHu3043jUxMXHsOcuyNDg4WHxsmqak451WJcnn8+nSpUslPwAAoH257vNRCBAFtm1rYmKi2PKRzWZlGIZM05Tf7y8JK4VtzPMBeGDjy7NtAwCPnGqej3Q6rWg0qsnJST148KBkjo9YLKbJyUnNz88rGAxqenpaqVRKkrS8vFycHwRAnfQPST390ucz1cv19B+WBYAGOVX4ME1T8XhckhQOh0u2HZ1s7N3tlUa6AKghIyD99LfSzovq5fqHqi4nDwD1dqrwAaDJGQGCBYCmx8JyAADAU4QPAADgKcIHAADwFOEDAAB4ivABAAA8RfgAAACeInwAAABPET4AAICnCB8AAMBThA8AAOApwgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BThAwAAeIrwAQAAPPWdRlcAANxafbZ9pm0AmgvhA0DTG7x4QX093bqz9Khqub6ebg1evFCz4/qcVenpQNltvRvbGtZGzY4FdBLCB4CmN2L0ybp7U5uvXlctN3jxgkaMvnMf702vXzt5nwL3b0v3y5cZk2T5fFrfnpR0+dzHBDoJ4QNASxgx+moSLNw4GBhRaP9n+uVPRjV2tXzLx/pXjxS4f1vdezlP6gS0E8IHAJTxVFe0d+WH0nD5Vo395/QxAc6K0S4AAMBThA8AAOApwgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKcIHwAAwFPMcAo0A2dd2nlRcbPPWfWwMgBQX4QPoNGcdemzj6SDnYpFApJ28j696fV7Vy8AqBPCB9BoOy8Og8ePF6Ur18sWWX2+rT/76zWlBkY8rhwA1B7hA2gWV65LwzfKbtrLv9RTvfS2PgBQJ3Q4BQAAniJ8AAAATxE+AACAp+jzAaDuVp9tn2kbgPZE+ABQN4MXL6ivp1t3lh5VLdfX063Bixe8qRSAhiN8AKibEaNP1t2b2nz1umq5wYsXNGL0eVQrAI1G+ABQVyNGH8ECQAk6nAIAAE8RPgAAgKcIHwAAwFP0+QCAc/A5q9LTgcoF+ockI+BdhYAWQPgAgDN40+vXTt6nwP3b0v0qBXv6pZ/+lgACvIPwAQBncDAwotD+z/TLn4xq7GqFlo+NL6XPZw5XLiZ8AEWEDwA4o6e6or0rP5SGLze6KkBLocMpAADwFOEDAAB4itsuQD0564f3+6vZ+NKbuqBxqrzHvRvbGtaGh5UBGo/wAdSLsy599pF0sHNy2Z7+wyGZaC/9Q4fv7eczFYuMSbJ8Pq1vT0qi7wg6A+EDqJedF4fB48eL0pXr1csyF0RTWn22faZtRUbgcJhtldav9a8eKXD/trr3cmepItCSCB9AvV25Lg3faHQtcAqDFy+or6dbd5YeVS3X19OtwYsXqu/MCFQNlvvPXYQYoM0QPgDgiBGjT9bdm9p89bpqucGLF1ixFzgDwgcAlDFi9BEsgDphqC0AAPDUqVo+bNtWJpORaZqybVuzs7MyDKNs2Ww2K8uyJEkPHjzQ4uJixbIAAKBznCp8RCIRraysSDoMIjMzM0qn02XLWpal+fl5SVIikdAf/dEfFV8LAAA6l+vbLrZtlzw2TbPYsnFUNptVLBYrPg6Hw8pms8f2AQAAOo/r8GFZlvx+f8lzfr9f2Wz2WNlgMKjFxcXiY8dxiuUBAEBnc33bpRAgjsrlyk+MEw6Hi/9eWlpSKBQq2+djf39f+/v7xcdbW1tuqwQAAFrQuUe7VAol727PZDIV+4bEYjFdvny5+BMIMMsjAADtzHXLh2EYx1o5crnciSNYotGolpeXK5ZbWFjQxx9/XHy8tbVFAAGAMp44u0x8hrbgOnyEQiElk8ljz09MTFR8TSKRUDQalWmaxRaSoyHE5/PJ5/O5rQYAdKQnzq5CP/+Ndg/eVC3X19Mt6+5NAgiamuvwYZpmyWPbtjUxMVEME9lsVoZhFMtlMhkFg8Fi8Lh3755mZ2drV3MA6CCbr15r9+CNPp2+obFrA2XLrD7b1p2lR9p89ZrwgaZ2qnk+0um0otGoJicn9eDBg5J+HLFYTJOTk5qfn5dt24pEIiWvNQyD8AEA5zR2bUAfjlxudDWAczlV+DBNU/F4XFLpaBZJJUHENE3l8/kaVA8AALQb1nYBAACeInwAAABPET4AAICnCB8AAMBThA8AAOApwgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKdOtbYLAKA+fM6q9LT8arWS1LuxrWFteFgjoH4IHwDQQG96/drJ+xS4f1u6X7ncmCTL59P69qQkVrVFayN8AEADHQyMKLT/M/3yJ6Mau1q55WP9q0cK3L+t7r2ch7UD6oPwAQAN9lRXtHflh9Jw5RaN/efbHtYIqC86nAIAAE8RPgAAgKcIHwAAwFOEDwAA4CnCBwAA8BThAwAAeIqhtgDQBFafVR9K+zy3qzGP6gLUG+EDABpo8OIF9fV0687So6rlftD1WD/ySZf6erypGFBHhA8AaKARo0/W3ZvafPW6arnejcvSr6VrAz6PagbUD+EDABpsxOjTiNFXvVBX5anXgVZDh1MAAOApwgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8x1BYox1mXdl5UL9M/JBkBb+oDAG2E8AEc5axLn30kHexUL9fTL03/ldR/pfz2jS9rXzcAaAOED+ConReHwePHi9KV6xXKbEhL/1r6D/+y+r56+g9bSAAARYQPoJIr16XhG5W3//S33JqB96q0qPVubGtYGx5WBjgbwgdwVkaAYAHv9A8dtqR9PlOxyJgky+fT+vakpMueVQ04LcIHALQCI3Bia9v6V48UuH9b3Xs5DysGnB7hAwBaxQmtbfvPtz2sDHB2zPMBAAA8RfgAAACe4rYL0ASeOLvafPW64vbVZzSnA2gfhA+gwZ44uwr9/DfaPXhTtVxfT7cGL17wqFYAUD+ED6DBNl+91u7BG306fUNj1wYqlhu8eEEjRp+HNQOA+iB8AE1i7NqAPhxhbgYA7Y8OpwAAwFO0fABndFInUYlbJQBQDuEDOIPTdBK17t4kgADAOwgfwBm46SS6+mxbd5YeafPVa8IHALyD8AGcA51EAeD06HAKAAA8RfgAAACeInwAAABPET4AAICn6HAKVLD6fFt7+Zflt7HQGwCcGeEDOOLZ9r6uSbr9q0f6XYXwIbHQGwCcFeEDOGJr90DXJP35H/+hrl7/qGI5Zi8FgLMhfAAVBPx9GmMODwCoOTqcAgAAT9HygY5z0oJwz3O7GvOwPgDQaQgf6ChuFoT7Qddj/cgnXerr8bBmANA5CB/oKG4WhOvduCz9Wro24PO4dgDQGQgf6EhVF4TrKh9KzuqkOUGYMwRApyF8AHUyePGC+nq6dWfp0YllmTMEQCc5VfiwbVuZTEamacq2bc3OzsowjIrls9msZmZmtLKyct56Ai1nxOiTdfdm1c6tBcwZAqCTnCp8RCKRYpCwbVszMzNKp9NlyxZCSjabPX8tgRY1YvQRKtB0ThrxJRGIUV+uw4dt2yWPTdOUZVkVy4fD4bPXCgBwZj5nVXpavu/Ss+19/eu/+kr2wWDVffT1dMu6e5MAgrpwHT4sy5Lf7y95zu/3K5vNKhgM1rxiAIDTedPr107ep8D929L98mWuSfqb93z6f/77/6iRP7hetszqs23dWXqkzVevCR+oC9fhw3Gcss/ncrlzVWB/f1/7+/vFx1tbW+faHwB0qoOBEYX2f6Zf/mRUY1fLt3ysf/VIgfu3NTqwz/IBaJhzj3apFErcisVi+uSTT85bDQCApKe6or0rP5SGyweL/ecM7UbjuV7bxTCMY60cuVyu6mgXNxYWFvTy5cviz/r6+rn2BwAAmpvr8BEKhco+PzExca4K+Hw+Xbp0qeQHAAC0L9fhwzTNkse2bWtiYqLY8pHNZo+NiCk4760ZAADQPk7V5yOdTisajWpyclIPHjwomeMjFotpcnJS8/Pzkg5HxywvL5dsY/gtALQPN/OFSMwZguNOFT5M01Q8Hpd0fB6Po5ONhUIhhUKhYnkAQPtws0J0AXOG4CjWdgEAnJqbFaIl5gxBeYQPdJxhbah3428rr1678aW3FQJaWNUVooEKCB/oKD3bT2T5/kL9v94/oWC/1D/kTaUAoMMQPtBRuvdy6u/a1/qPfqHA929ULtg/JBkBz+oFAJ2E8IGOtG+MScM3Gl0NAOhIruf5AAAAqAXCBwAA8BThAwAAeIo+HwCAult9Vn01XWZB7SyEDwBA3QxevKC+nm7dWXpUtRyzoHYWwgcAoG5GjD5Zd29WXQOGWVA7D+EDADqQz1mVnpaf5bd3Y1vD2qjZsUaMPkIFShA+AKCDvOn1ayfvU+D+bel++TJjkiyfT+vbk5KYOh21R/hAe3HWpZ0XFTf7nFUPKwM0n4OBEYX2f6Zf/mRUY1fLt3ysf/VIgfu31b2X87h26BSED7QPZ1367CPpYKdikYCknbxPb3r93tULaDJPdUV7V34oDZdv1dh/Xn1kCnBehA+0j50Xh8Hjx4vSletli6w+39af/fWaUgMjHlcOAFBA+ED7uXK94rote/mXeqqX3tYHAFCCGU4BAICnCB8AAMBThA8AAOApwgcAAPAUHU7RGk6Yv0OStPGlN3UBAJwL4QPNz8X8HUU9/VL/UP3rBAA4M8IHmp+L+TuK+ockI+BNvQAAZ0L4QOuoMn8HgNrzcvE5dBbCBwCgBIvPod4IHwCAEiw+h3ojfAAAjmHxOdQT83wAAABP0fIBADizap1SixiFhiMIHwCAU3PTKbWop1/66W8JICgifABAm1l9Vrk/RrVtp+GmU6qkw5mHP585nK+H8IF/QvgAgDYxePGC+nq6dWfpUdVyfT3dGrx44dzHO6lTKlAJ4QONd9K6LazZArgyYvTJuntTm69eVy03ePGCRow+j2oFHEf4QGO5XbeFNVsAV0aMPoIFmh7hA43ldt0WessDQNsgfKA5sG4L0N6q3D5lnZjOQ/gAANRP/9DhbdPPZyoWYZ2YzkP4QFNYfb6tvfzLitvpIAe0KCNwOMdHlU7lrBPTeQgfaKhn2/u6Jun2rx7pd1XCR19Pt6y7NwkgQCsyAlX7bLFOTOchfKChtnYPdE3Sn//xH+rq9Y/Klll9tq07S4+0+eo14QMA2gDhA00h4O/T2Aj3egGgE7CqLQAA8BThAwAAeIrwAQAAPEWfD7SVJ85u1XUtarWiJ4DaW8/tau8JQ+47AeEDbeOJs6vQz3+j3YM3VcvVakVPALVxqa9HkvTv/q//T7/7Pyv/8cCQ+/ZB+EDb2Hz1WrsHb/Tp9A2NXRuoWI6/noDmcm3AJ0n6xZ/e0N6VH5Ytw5D79kL4QNsZuzagDxm2C7ScsasD0jCf3U5A+GgjJ/V3kPirHwDQeISPNnGa/g6e3jN11quu6eBzVr2pBwCgaRA+2oSb/g6e3zN11qXPPpIOdioWCUjayfv0ptdf//oAaHknjVijdbc1ED7aTFP1d9h5cRg8frwoXbletsjq82392V+vKTUw4nHlAJyk2oXe62HrgxcvqK+nW3eWHlUtx4iY1kD4QP1duS4N3yi7aS//Uk9VeVw/AO+d5kLv1bD1EaNP1t2bJ87jw4iY1kD4AACUcHOhl7y/xTGiDY10Ve5D1vvetoa14Vl9cHaED7SMk5p5mb0UqJ0Ro6+5Wg9c9CEbk2T5fFrfnpTUJLefURbhA03PbROwxOylQEvb+LL6thP6kK1/9UiB+7fVvZerUwVRK4QPnN0Jw2irfpGcgtsmYIme7kBL6h+Sevqlz2eql+vpl773zyQjUHbz/nNaP1sF4QNn46IJVNLhl0X/0LkP13RNwABqxwhIP/1t9T9mpMPvkgrBA62F8IHy3LRqnNAEKokvCwDuGAG+KzoI4QPHnaZVo0oTKAAA5RA+cJyLycEk6R+/vaiNV5ekV5Xn6aAPBgDgKMIHKqsyOdjv15Kp3qmU2QYBAEedKnzYtq1MJiPTNGXbtmZnZ2UYxrnLwjvD2lDvxt9KXeXXf5HkapRKU64lAwBunNSnrYA+a3VzqvARiUS0srIi6TBczMzMKJ1On7tsx3P7Qaiid2NbP+h6rN6NyxWDxcD638vy/YX6f71/8g5djlJpqrVkAOAkbvu0SYffg9N/JfVfqVyGgHImrsOHbdslj03TlGVZ5y7b8U7zQahiTNJ/9En6deUyfyBpRz793Z/8lf7ge9+rvkM+UABalM9ZlZ5WaN11O1JvZ0Nvf/Wv9N5/+JdVj/X2O3167396wPflKbkOH5Zlye8vXfbc7/crm80qGAyeuWzLqlWzncvOnSdZfb6t2796pF/86Q2NXa1wG6SwgmzgpjRMawWA9vKm16+dvE+B+7el+5XLvf1On77yfaiDfOXVtF+8HdEn+/9Ofd86FcuMdT3RL/S/KfdffiP/9z48R83VcX/wuQ4fjuOUfT6XOz6N7WnK7u/va3//97cBXr48HDmxtbXltmqn8uLrv5ez8eRc+/jO3qYC9/9nvfdm78Syb7t7tf6j/1Xf9g6W3e57aev9/bzWvrms/Qtnn4zL3urV3+79V9rw/Te6NlA+WGxsvdQ/7H+t/9f+WtvfnO/3az9/pbf7O9r+ZktbW11ly2x/s6W3+ztVj+dmPwDgRr5vSH+y+2/VWyUwSJKT/67+8X//naTfVS3X29OnT//0v5W/v6fs9qd/v6p//L8X1f9//I867xXrpGtFrRlXRjT0X5/QAn5Khet2Pp8/sey5R7tUChpuy8ZiMX3yySfHng8E2iUBfiP925+cXOx/+ZOaHO2ffXpymf/BRRmvj+dmPwDgtf/uZ14dyeW1ogV88803uny5euu66/BhGMaxlotcLld2BMtpyi4sLOjjjz8uPn779q1yuZyGhobU1VXbv4S3trYUCAS0vr6uS5cu1XTfzaDdz09q/3Pk/Fpfu58j59f66nWO+Xxe33zzjYaHh08s6zp8hEIhJZPJY89PTEycq6zP55PP5yt5rt5Dci9dutS2/6mk9j8/qf3PkfNrfe1+jpxf66vHOZ7U4lHwntsdmqZZ8ti2bU1MTBSDQjabLY5yOaksAADoXKfq85FOpxWNRjU5OakHDx6UzNsRi8U0OTmp+fn5E8sCAIDOdarwYZqm4vG4JCkcDpdsOxouqpVtFJ/Pp7/8y788dpunXbT7+Untf46cX+tr93Pk/FpfM5xjV97NmBgAAIAacd3nAwAAoBYIHwAAwFOEDwAA4CnCBwAA8NS5p1dvN1NTU1peXm50NWrGsixls1lJ0tramkZHR4vDoduBZVlaXl6W4ziybVuRSESzs7ONrlZdpFIpra2tFUeRtRrbthWPxzU6OirpcDLBdnuvWv09qqYTPmvt/n15VCOvd4SPdyQSCVmW1ehq1Ew2m5XjOCUfntHRUa2trZWdgbbVFL4oCl/0juNofHxcKysrbXF+0u8v2JJ07969lv2yt21b4+Pjevz4cXGywWg0qkQi0fJf7u3yHlXTCZ+1dv++PKrh17s88vl8Pr+2tpYPh8P5dvqVzM7O5k3TPPZcu5xjOBw+9lwymcxLyq+trTWgRvUVDAbz8/Pzja7GmczOzh6r++bmZtv8Xyxo5feomk74rLX79+W7muF6R5+Pf5LJZDQ9Pd3oatTU1NSUQqFQo6tRN5lMRtFotOS5wvpB7dSC1Q7u3btXvN1SUGgB4b1qfp3wWWv378t3NcP1jvChwzeiWWZhraVwOHysufDevXttcz86HA4fu6Ch+TiOI8dxjq35JB0GkMI9djSvTvistfv3ZUGzXO86vs+H4zjK5XIyTbNtvwQzmYwePHigbDardDrdNum+3HpBDx8+lKS2Ocd2UFhwshy/368XL154WBucRSd91tr1+1Jqrutdx7d8pFKptuwg9q5wOKx4PK65uTlFo9GqF4NWF4/HFY/Hy/6VjebkOE6jq4AzaNfPWjt/XzbT9a6jw4dlWW2Vak8SDocVCoU0Pj7ell/4kUhEoVCo5UdPtJtC345ycrmcdxVBzXTCZ63dvi+b7XrX8rddstmsZmZmXJdfXFxUMBgsvrbZPzznOb9ypqamlEgklEqlmuLca3V+qVRKfr+/KYfE1fo9bDV+v19S+RYOx3GqhhM0n2b+rNVas31fnkezXe9aPnwEg0GtrKyc+nWFyYDe7cFduAcWjUY1NDTUFG/UWc9PkgYHB7WwsFByHoULwdraWk3qd17nOb+CTCYjx3FKvgyb6aJWi3NsZYZhyDCMiq0cU1NTHtcIZ9Xsn7XzaIXvy7Nqxutdy4ePsyp33yuVSsmyrLbo3Vz4K/Po/djC/cvx8XGvq1QX2WxWuVyu5IPjOI4sy2qKHt04dOvWrWNf4IX/i83UFIzK2vmz1u7fl814vevY8FFOO9zXKyhMXX20+T6ZTCoYDDZNp6PzsG1bsVhM09PTymQyxeeXl5c1NzfXwJrVR2HIaiuKRqOampoq+aJLJpNt13Tfyu9RNe3+WeuE78ujGv3/tCufz+cbWoMmYNu2ksmkMpmMbNtWOBzW1NRUW/yHS6VSWllZkWEYsm1bpmm2RcuOdNhMWukD1C7/rR3HUSwWk+M4SqVSMgxDt27dask1J7LZrJaWljQ5OVn8i7LVzqGcdnqPKumEz5rU3t+XBc1yvSN8AAAAT3X0UFsAAOA9wgcAAPAU4QMAAHiK8AEAADxF+AAAAJ4ifAAAAE8RPgAAgKcIHwBQQ/WYObKdlnUHJMIHcG5zc3MlU06jczmOU7J4V61ks1lZllXz/QKNQvgAzsG2baVSqVNdcDKZjCKRSB1r1dzHbxb1+D3MzMzUZTrucDisdDpNCwjaBuEDOIdkMql4PC7btl23fiwtLRWX6m6ERh+/WdT695BKpTQ1NVW35eXj8XhbLOIGSIQP4Mxs2y4uHmYYhuvWD8uyNDU1VefaNe/xm0Wtfw/JZLKui3MZhiHTNLn9grZA+ADO6N2LzcLCgmzbPvHCUFhyPRQKeVHFpjt+s6j178GyLE1MTNRkX9VEIhElk8m6HweoN1a1Bc6gsHz6u0umDw4OyjRNraysHCtvWZaWl5eVzWb18OHDktBSr2b6Zjp+s6jX72Fubk5TU1MKh8PF57LZrGZmZmTbtr744otif40HDx5I0rG+IZlMRrlcTn6/X7lcTisrK4pEIscC0uDgoDY3N89cV6Ap5AGc2vz8/LHn4vF4XlJ+eXm54utmZ2fz4XC4nlWrqtHHbxa1/j2YpplfW1sru01SfnZ2Nr+5uVl8zjCMfDqdLj5eW1vLz87OlrwuHo+X/b9U7VhAq+C2C3BKjuNoaGjo2PNu+n6cpp9B4a9ptz9uOryedPxEIqGuri7Nzc0pm80WXzM+Pq7R0VGlUqmS8qlUSl1dXcfOOZFIaHBwUIlEouIxIpHIsf15pdb9PXK5nEzTLLvNMIziT8HExESxBUQ6bCU5OpLl3VaUd5mmyagXtL5Gpx+g1czPz5f8FXt0m6T8ysrKsW2bm5t5SQ37q9XN8dfW1vKSjp1fKBTKx+PxsuVnZ2fzhmEc21aufLVjeKUe70O1r9KjrRz5fD4fDodLWs82NzfzhmHkTdPMz8/PV209C4fD+WQyef5KAw1EywdwCoXZKyv1D1hYWJCksq0fDx8+LI5YaAQ3x89mszJN89j5WZZVtnNmNptVPB6X3+8vaXmxbbtiZ85Kx6hWp/Hxcdc/hRabSmr9PtRiRlPDMPT48WOFQiFlMhlNTU1pdHS0LrOlAs3gO42uANBKUqlUMWCUYxiG5ufnlUgklM1mFQwGi9uWl5dLRkQ4jlP1Ajw3N3eq5vW5ubmKTfVuj7+8vHwsNBQu5u+eS0Eul5NhGJqbm1MsFise37KsisNOyx2jmmAwWLYT71md9n04SS067Bb+rxRGsti2XfydHu2Y6jhOwwIsUCuED8Alx3H04sWLEy82CwsLSiQSikajWl5eLj6fzWZL+hkcHS1zVK2HVLo5vmVZCgaDJX0xlpeXywaPd83Ozioajcq27RMvjJZl1WUKcrdO+j1kMpmS86gW6AoMwzhXiLFtW7lcrhjKTNNUMpksO6lYYUQM0MoIH4BLsVhMQ0NDrjp2BoNBWZZ1rPWjcEFLpVJ1nZCqkmrHdxxHtm1reXm5JEC4aakwDEPhcFjxeLxqaCoco9HzjFT6PaRSKS0vLyudTheHyroJH4VOoCeFtILCPCPvisfjJb8X27bLdop1E/CAZkf4AFxKpVKnvgcfi8WUTqclHfYDKfx7YmLC8/k1Tjp+pb4QlmVpcXHx2P6O9gMpjM6JRqMVL47VjuFVIKn0e3AcR3Nzc8U5NILBoL744gtX+5yentbDhw9Lwkc2m1UymZTjOIrFYsrlcrp165ZisZgePnyohw8fSvp9C1ckElEikSipz9GWqcItl06amwXtiUnGAEg6vChns9mSW0WO4xQntTp6wSvXejI6OirTNJVOp8teIMsdQzocflvtFpQXMpmMYrHYmfqX2LZdEmrqpTARWSNazYBaYrQLAEnH+0JIhy0Slf7SLtcKVOj3Uekv83LHiEajDb8NIx3eOjral8LtOiqF2y71Hp0Si8V069atuh4D8ALhA+hwjuMokUjIsiytra0VR9gUWgIK/y4ojMRIJpPHLs63bt0q20fi6DFSqZQSiURxcjS3fSXqKRQKFTvbZjIZpVKpU/WtWFxcLP6+6iGbzSoUCnHLBW2B2y4AUCOJRKIYYmrJcRzNzMzU/bYO4BVaPgCgRubn52VZVs1vv6RSqbKdfoFWRcsHAADwFC0fAADAU4QPAADgKcIHAADwFOEDAAB4ivABAAA8RfgAAACeInwAAABPET4AAICnCB8AAMBT/z+RoU+dcxLpfAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "\n",
    "beam_spill_times = norm.rvs(loc=0.0, scale=1.0, size=timing.shape[0])  # nanoseconds\n",
    "\n",
    "timing_with_spill = 1e-9*beam_spill_times + timing\n",
    "nu_timing = 1e-9*beam_spill_times + tof_light\n",
    "\n",
    "timing_bins = np.linspace(-4, 4, 50)\n",
    "\n",
    "\n",
    "plt.hist(1e9*(timing_with_spill - tof_light), weights=weights,\n",
    "         bins=timing_bins, histtype='step', label=\"HNLS\", density=True)\n",
    "plt.hist(1e9*(nu_timing - tof_light), bins=timing_bins, histtype='step', density=True, label=\"Neutrinos\")\n",
    "#plt.xscale('log')\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r\"$\\Delta t = t_{HNL} - t_c$ (ns)\", fontsize=14)\n",
    "plt.legend()\n",
    "plt.xticks(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sbn_constants import *\n",
    "def get_timing_distribution_sbnd(mN, mZp, Ualpha=1.0, flavor=1):\n",
    "\n",
    "    flux = FluxHNLFromProtonBrem(proton_energy=8000.0, mixing_flavor=flavor,\n",
    "                                mixing_angle=1.0, hnl_mass=mN, coupling_BL=1e-4,\n",
    "                                zprime_mass=mZp, det_dist=SBND_DIST, det_length=SBND_LENGTH,\n",
    "                                det_area=SBND_AREA, target=Material(\"Be\"), n_samples=100000)\n",
    "    \n",
    "    flux.simulate()\n",
    "    flux.decay_to_hnl()\n",
    "    flux.propagate(Ualpha)\n",
    "\n",
    "    return flux.hnl_timing, DUNE_EFFICIENCY * DUNE_POT_TOTAL * flux.decay_axion_weight\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_flux.py:394: RuntimeWarning: invalid value encountered in sqrt\n",
      "  zp_p = sqrt(zp_energy**2 - self.ma**2)\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_flux.py:397: RuntimeWarning: invalid value encountered in arcsin\n",
      "  theta = arcsin(sinTheta)\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/zprime_decay_widths.py:22: RuntimeWarning: invalid value encountered in sqrt\n",
      "  return np.nan_to_num((gBL**2 * mZp / 12 / pi) * sqrt(1 - 4*power(ml/mZp, 2)) * (1 + 2*power(ml/mZp, 2)))\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_flux.py:417: RuntimeWarning: invalid value encountered in sqrt\n",
      "  pzp = sqrt(self.axion_energy[i]**2 - self.ma**2)\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:155: RuntimeWarning: invalid value encountered in sqrt\n",
      "  * sqrt((1-power((M_PI - ml)/mN, 2))*(1-power((M_PI + ml)/mN, 2))))\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:169: RuntimeWarning: invalid value encountered in sqrt\n",
      "  kallen = sqrt(kallen_alplib(1.0, xP**2, xl**2))\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:126: RuntimeWarning: invalid value encountered in sqrt\n",
      "  log_arg_num = np.clip((1 - 3*xl2 - (1-xl2)*sqrt(1 - 4*xl2)),a_min=1.1e-16,a_max=np.inf)\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:127: RuntimeWarning: invalid value encountered in sqrt\n",
      "  log_arg_denom = ((1+sqrt(1-4*xl2))*xl2)\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:141: RuntimeWarning: invalid value encountered in sqrt\n",
      "  * ((1 - 14*xl2 - 2*xl4 - 12*xl6)*sqrt(1-4*xl2) + 12*xl4*(xl4-1)*logterm) \\\n",
      "/home/thompson/Dropbox/coding/pyALPs/HNL/hnl_decays.py:142: RuntimeWarning: invalid value encountered in sqrt\n",
      "  + 4*(c2*(1-delta) + c4*delta) * (xl2*(2 + 10*xl2 - 12*xl4)*sqrt(1-4*xl2) + 6*xl4*(1-2*xl2+2*xl4)*logterm)))\n"
     ]
    }
   ],
   "source": [
    "test_mN = 200.0  # MeV\n",
    "test_mZp = 2.1*test_mN\n",
    "test_Ualpha = 1e-4\n",
    "tof_light_sbn = (SBND_DIST / (1e-2 * C_LIGHT))\n",
    "\n",
    "\n",
    "timing_sbn, weights_sbn = get_timing_distribution_sbnd(test_mN, test_mZp, test_Ualpha, flavor=2)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi0AAAG5CAYAAABC9gS6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAl70lEQVR4nO3dv28bZ+Ln8Y/8Syt7I41obLAbrQLseMsUezRZHBDAhanqUpL2/QOiqhyMGCChauGKIA8LLJCKVHedTJbpOCkWCK6RNHCR0pwURrILByZHjL2x7G/CK7ScE8UfmtEPUg/5fgEEzJmH8zwzNDkfPc8zw7lOp9MRAADAJXdl0g0AAAAIg9ACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADDCtUk34Lz8+uuv+uGHH/TBBx9obm5u0s0BAAAhdDod/fTTT/roo4905crovpSpCS0//PCDVldXJ90MAABwCi9evNAf//jHkWWmJrR88MEHkg53enFxccKtAQAAYbTbba2urgbn8VGmJrR0h4QWFxcJLQAAGCbM1A4m4gIAACMQWgAAgBGmZngIAIDz8ssvv+j9+/eTbsbUuHHjxolXBoVBaAEA4D86nY7+9a9/yff9STdlqly5ckV/+tOfdOPGjTNth9ACAMB/dAPLhx9+qJs3b3Lfr3PQvY/aP//5T3388cdnOqaEFgAAdDgk1A0st2/fnnRzpsrvfvc7/fDDD/qv//ovXb9+/dTbYSIuAABSMIfl5s2bE27J9OkOC/3yyy9n2g6hBQCAIxgSOn/ndUwjDQ95nqdarSbbtuV5nrLZrCzLilzWcRxJku/72tnZ0cOHDxWPxyPXAQDAOHzv/6zWm3djq2/51g2tWAtjq88UkUJLJpPR3t6epMNwsb6+rmq1GrlsJpPR119/rVQqpWazqUwmo0ajEbkOAAAu2vf+z0r97R/6+f3ZhjaiWLh+Vc7je6GCi+u6KpfLqlQqyuVyevjwoWzb1tOnT5XP55VIJJTJZJRIJLS9va1SqaRyuaxsNivp8FxbLBblOE5Qvru9bDartbU1pdPpvnq7nQzdjoWjnQ0XJXRo8Tyv57lt20GPSdSy1Wo16FmRFOxwlDoAABiH1pt3+vn9L/r7w7/ozx/+9sLre/7ytR5tP1PrzbtQoSUejyufz6tSqWhzczM4p2azWZXLZa2trQVBwrZt+b6vjY0NPXjwQJZlybZtlctllUqloFx3e8Vicehox9FOBkkqlUp69erV2Xb+BKFDi+M4isViPctisZhc1+0JIGHKplKpYHm1WtXGxkbkOgAAGKc/f/hbfbKyNOlmDHT83DlKJpNRs9nsG8k4Gk5O2t7xTgZJyuVyKpVKodtxGqEn4g670U6z2TxVWdd1lc/nexJglDoODg7Ubrd7HgAA4GRbW1tyHEe1Wu1Ur+8OBVUqlZ7lFzk0JJ3DfVqi3DXwaNl4PC7btpXP51Wr1QaOl42qo1Ao6MmTJxFaClwOYSb0MQkPwGlUKpWeHpNh52jLslQsFrW+vj7y/DvK1taWMpmMNjY2lEqllM/ne0ZSLkLo0GJZVl+PR7PZHDjWFbasZVnKZDJaW1tTq9WKVMfm5qa++OKL4Hm73dbq6mrY3QEmIuyEviiT8ACg6/gVt+VyeWTZ7hSNUeWGSafTajQachxH9Xpda2trqlarpw5BYYQeHhqWnhKJRKSyjuNoeXk5WGbbtqTD8bEodczPz2txcbHnAVx2Ryf0ffX5pwMff3/4F/38/pexXl4JYDZ1rxJyXTdU+e5wUrcHx7btIPyUy2UVCoWLaqqkCKGlGy66PM9TIpEIEp3rusHEnFFlY7FYTzhxXVeWZQXDRaPqAKZFd0LfoMc4rk4AAOnwfF0sFpXJZEKV757nPc/rCzoPHjy48B+ajDSnpVqtKp/PK5lMamdnp2fWcaFQUDKZVC6XG1k2Ho/r4cOHweSder3ec8nUqDoAAEC/QReshC2Xy+W0vb3dEziGXWRTr9eD83w+n1e9Xg/WO45zoUNDUsTQ0k1kkvoadjxcjCp79PnxmcajXgcAwKQ8f/n6UtbTvbmcdNiBcPTmcp7naXt7W7ZtBxe/NJtNbW5u9pxjt7a2tLu7O3B7d+7cUaPRCO5W35XJZFQqlYLRkEajEZy/L8pcp9PpXGgNY9Jut7W0tKT9/X3mt+DS+vb7fX325Tf66vNPh97vIUwZAOfv7du3+u677/SnP/1Jv/nNb4Lll/2OuCYYdmylaOfvM1/yDADANFuxFuQ8vsdvD10ChBYAAE6wYi0QIi6B0FcPAQAATBKhBQAAGIHQAgDAEVNyfcqlcl7HlNACAICk69evS5L+/e9/T7gl0+fdu8NJzFevXj3TdpiICwCADk+olmXp5cuXkqSbN29qbm5uwq0y36+//qoff/xRN2/e1LVrZ4sdhBYAAP7j97//vSQFwQXn48qVK/r444/PHAIJLQAA/Mfc3Jz+8Ic/6MMPP9T79+8n3ZypcePGDV25cvYZKYQWAACOuXr16pnnX+D8MREXAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGCEa5NuAIDBnr98PXL98q0bWrEWxtQaAJg8QgtwySzfuqGF61f1aPvZyHIL16/KeXyP4AJgZhBagEtmxVqQ8/ieWm/eDS3z/OVrPdp+ptabd4QWADOD0AJcQivWAmEEAI5hIi4AADACoQUAABgh0vCQ53mq1WqybVue5ymbzcqyrMhlXdeV4ziSpJ2dHW1tbfWsk6R4PC7P8+T7vuLx+On2DgAATI1IoSWTyWhvb0/SYShZX19XtVqNXNZxHOVyOUlSqVTS/fv3g7LlclmVSkWSlEqlhm4fAADMltDDQ57n9Ty3bTvoLYlS1nVdFQqFYF06nZbrusFr7t69q1arpVarpXq9PrQnBwAAzJbQocVxHMVisZ5lsVgsGM4JWzYej2traytY7vt+sL7LsizCCgAA6BF6eKgbLo5rNpuRy6bT6WDZ9va2UqlUEFJ831etVpN0ON9lY2NDtm33bevg4EAHBwfB83a7HWY3AACAoc58n5ZhASVM2W5A6c5nkdQzYde2ba2tranRaPRtq1Ao6MmTJ6dpMgAAMFDo4SHLsvp6VZrN5sBhnLBl8/l837yVo/NhulceHZ8jI0mbm5va398PHi9evAi7KwAAwEChQ0sqlRq4PJFInKpsqVRSPp+XbdvyfV++78t1Xd2/f7/vdcfnx0jS/Py8FhcXex4AAGB6hQ4tx+eVeJ6nRCLRc3+Vbo/ISWVrtZri8XgQWJ4+fSrLsmTbtorFYvA6x3GUTqeZlAsAAKLNaalWq8rn80omk9rZ2em5h0qhUFAymQzuvzKsrOd5ymQyPdu1LCuYy5JIJFQqlWRZlhqNBvdpAQAAkqS5TqfTmXQjzkO73dbS0pL29/cZKsKl9e33+/rsy2/01eef6pOVpYlvBwAmLcr5m98eAgAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBGuTboBAE7v+cvXJ5ZZvnVDK9bCGFoDABeL0AIYaPnWDS1cv6pH289OLLtw/aqcx/cILgCMR2gBDLRiLch5fE+tN+9Glnv+8rUebT9T6807QgsA4xFaAEOtWAsEEQAzhYm4AADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACNeiFPY8T7VaTbZty/M8ZbNZWZYVuazrunIcR5K0s7Ojra2tYF2UOgAAwOyIFFoymYz29vYkHYaL9fV1VavVyGUdx1Eul5MklUol3b9/PygbpQ4AADA7Qg8PeZ7X89y27aC3JEpZ13VVKBSCdel0Wq7ryvO8SHUAAIDZEjq0OI6jWCzWsywWi8l13Uhl4/G4tra2guW+7wfro9RxcHCgdrvd8wAAANMrdGjphovjms1m5LLpdDpYtr29rVQqJcuyItVRKBS0tLQUPFZXV0fvAAAAMNqZrx4aFjTClPV9X7Va7cQ5K4Pq2Nzc1P7+fvB48eJF6HYAAADzhJ6Ia1lWX49Hs9kceGVP2LL5fF71ej1YHqWO+fl5zc/Ph20+AAAwXOiellQqNXB5IpE4VdlSqaR8Pi/btuX7vnzfj1QHAACYLaFDi23bPc89z1Mikei590r36p+TytZqNcXj8SCwPH36VJZlnfg6AAAwuyLdp6VarSqfzyuZTGpnZ6dnLkqhUFAymQzuvzKsrOd5ymQyPdu1LEvZbPbEOgAAwOya63Q6nUk34jy0220tLS1pf39fi4uLk24OMNC33+/rsy+/0Veff6pPVpamrj4AiCrK+ZvfHgIAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGOHapBsATJPv/Z/VevNu6PrnL1+PsTUAMF0ILcA5+d7/Wam//UM/v/9lZLmF61e1fOvGmFoFANOD0AKck9abd/r5/S/6+8O/6M8f/nZoueVbN7RiLYyxZQAwHQgtwDn784e/1ScrS5NuBgBMHSbiAgAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADBCpF959jxPtVpNtm3L8zxls1lZlnWqsq7ran19XXt7ez2vc11XkhSPx+V5nnzfVzwej7ZXAABg6kQKLZlMJggZnudpfX1d1Wo1ctlumOkGlKPK5bIqlYokKZVKDd0+AACYLaFDi+d5Pc9t25bjOKcqm06nh9Zz9+5dtVotSRraiwMAAGZP6DktjuMoFov1LIvFYgN7S6KUHcSyrBMDy8HBgdrtds8DAABMr9Chxff9gcubzeaZyg56ba1WU61WUz6f7+u16SoUClpaWgoeq6urJ24bAACYK9KclkGGBZTTlj06Yde2ba2tranRaPSV29zc1BdffBE8b7fbBBcAAKZY6J4Wy7L6ekqazebAYZwoZY872rPSvfJoUG/L/Py8FhcXex4AAGB6hQ4tqVRq4PJEInGmske5rqv79+/3LT8+PwYAAMye0KHFtu2e557nKZFIBL0nrusGPSInlT3q6JCRbdsqFovBc8dxlE6nuYoIAABEm9NSrVaVz+eVTCa1s7PTcw+VQqGgZDKpXC53YlnHcVSv13te1w0niURCpVJJlmWp0WhwnxYAACBJmut0Op1JN+I8tNttLS0taX9/n/ktmIhvv9/XZ19+o68+/1SfrCxNujmSLmebAOCoKOdvfnsIAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwwrUohT3PU61Wk23b8jxP2WxWlmWdqqzrulpfX9fe3t6p6wAAALMjUmjJZDJByPA8T+vr66pWq5HLdkOJ67pnqgMAAMyO0KHF87ye57Zty3GcU5VNp9NnrgMAAMyW0HNaHMdRLBbrWRaLxQb2lkQpe9rXHRwcqN1u9zwAAMD0Ch1afN8fuLzZbJ6p7GlfVygUtLS0FDxWV1dHbhsAAJjtzFcPDQsaZy170us2Nze1v78fPF68eHGqbQMAADOEntNiWVZfj0ez2Rx4ZU+Usqd93fz8vObn50O1HZh1z1++Hrl++dYNrVgLY2oNAJxO6NCSSqVULpf7licSiTOVPY/XARhs+dYNLVy/qkfbz0aWW7h+Vc7jewQXAJda6NBi23bPc8/zlEgkgl4Q13VlWZZs2z6x7FG+7wfLo7wOwMlWrAU5j++p9ebd0DLPX77Wo+1nar15R2gBcKlFuk9LtVpVPp9XMpnUzs5Oz/1TCoWCksmkcrnciWUdx1G9Xu95Xfcy6FGvAxDdirVAGAEwFeY6nU5n0o04D+12W0tLS9rf39fi4uKkm4MZ9O33+/rsy2/01eef6pOVpUk3JzRT2w1gOkQ5f/PbQwAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACMQWgAAgBEILQAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACNci1LY8zzVajXZti3P85TNZmVZVuSyo9a5ritJisfj8jxPvu8rHo+fegcBAMB0iBRaMpmM9vb2JB0Gj/X1dVWr1chlR60rl8uqVCqSpFQqNXT7AABgtoQOLZ7n9Ty3bVuO40Que9J27t69q1arJUlDe3EAAMDsCT2nxXEcxWKxnmWxWCwYzglbNsx2LMsisAAAgB6he1p83x+4vNlsRip70nZ831etVpMk7ezsaGNjQ7Zt95U/ODjQwcFB8Lzdbo9oPQAAMF2kOS2DDAshUct21x2dlGvbttbW1tRoNPrKFwoFPXnyJEJLAQCAyUIPD1mW1der0mw2Bw7jjCp70naOznnpXl10fB6MJG1ubmp/fz94vHjxIuyuAAAAA4UOLalUauDyRCIRqeyoda7r6v79+33rjs+BkaT5+XktLi72PAAAwPQKPTx0fF6J53lKJBI991exLEu2bY8se7xn5ug627ZVLBaDdY7jKJ1OMykXAABEm9NSrVaVz+eVTCa1s7PTcw+VQqGgZDKpXC53Ytlh6yzLUiKRUKlUkmVZajQa3KcFGJPnL1+PXL9864ZWrIUxtQYA+s11Op3OpBtxHtrttpaWlrS/v89QESbi2+/39dmX3+irzz/VJytLk25OaN/7Pyv1t3/o5/e/jCy3cP2qnMf3CC4AzlWU8/eZrx4CYLYVa0HO43tqvXk3tMzzl6/1aPuZWm/eEVoATAyhBYBWrAXCCIBLj195BgAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxwbdINAGCO5y9fj1y/fOuGVqyFMbUGwKwhtAA40fKtG1q4flWPtp+NLLdw/aqcx/cILgAuBKEFwIlWrAU5j++p9ebd0DLPX77Wo+1nar15R2gBcCEILQBCWbEWCCMAJoqJuAAAwAiEFgAAYARCCwAAMAKhBQAAGIHQAgAAjEBoAQAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAACMQGgBAABGILQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACNcm3QDAEyX5y9fj1y/fOuGVqyFMbUGwDQhtAA4F8u3bmjh+lU92n42stzC9atyHt8juACIjNAC4FysWAtyHt9T6827oWWev3ytR9vP1HrzjtACILJIocXzPNVqNdm2Lc/zlM1mZVlW5LKnXQfgcluxFggjAC5MpNCSyWS0t7cn6TBcrK+vq1qtRi572nUAAGB2hb56yPO8nue2bctxnMhlT7sOAADMttA9LY7jKBaL9SyLxWJyXVfxeDx02d3d3VOtO14HAHOddIWRxFVGAPqFDi2+7w9c3mw2I5U97brjDg4OdHBwEDzf39+XJLXb7YHbOKsf22/14+uDkwtiZnk/vtGvB//W65/aarfnJt2cS+naL29149e3+l//5/+eWPY316/o7//zvyl28/oYWgYgjN/9dl6/W/zNuW6ze97udDonlj3z1UPDgkbUslHXFQoFPXnypG/56upq6PYAF+G//33SLZge/+N/T7oFAMblp59+0tLS0sgyoUOLZVl9PR7NZnPglT2jyp523XGbm5v64osvgue//vqrms2mbt++rbm58/0rt91ua3V1VS9evNDi4uK5btsEs77/EseA/Z/t/Zc4BrO+/9LFHYNOp6OffvpJH3300YllQ4eWVCqlcrnctzyRSEQqa9v2qdYdNz8/r/n5+Z5lF31p9OLi4sz+Z5XYf4ljwP7P9v5LHINZ33/pYo7BST0sXaFDi23bPc89z1MikQiCguu6sixLtm2PLHs8WIRdBwAAZlukOS3ValX5fF7JZFI7Ozs9908pFApKJpPK5XInlj3tOgAAMMM6ONHbt287f/3rXztv376ddFMmYtb3v9PhGLD/s73/nQ7HYNb3v9O5HMdgrtMJcY0RAADAhIW+Iy4AAMAkEVoAAIARCC0AAMAIhBYAAGCEM9/GH4fW1tZUr9cn3YyxcRxH9Xpdvu/L8zxlMhlls9lJN+vceZ6nYrGoO3fuSDq8geE07ucws/I+RzFrn3Xp8HNQLpd1+/ZtvXr1SslkUul0etLNGhvXdeU4jiTp1atXun37dnB7j2lUqVTUaDRULBb71k38O3Fi1y1NkWKx2JmlQ1mv1zvFYjF43mq1OrZtd7LZ7ARbdf4ajUbHsqxOq9UKluVyuZ59n2az8j5HMWuf9U7n8P9BOp0OnrdarU48Hp9gi8ar0Wj0feb39vZ6jsk0aDQanWw228lmsx3Lsjq5XG5gmUl/J87Wp+8CNBqNTjqdnqkvskEf1nK53JHUaTQaE2jRxchms30f3FarNTPv9ay8z2HN4me91Wr1naTq9XrHsqzJNWrMstnswP/vqVRqAq0Zj3g8PjC0XIbvROa0nFGtVtPDhw8n3YyxqtVqyufzPcu6vw/V7UKdBk+fPg26QLu6PykxTfs5zKy8z2HN4me9UCj0/ZRKKpVSq9WaXKPGrNlsDhwmOf7jvrPgMnwnElrOoFarzdS4blc6ne77jzttfN+X7/t9v6MlHX5IXdedQKvGaxbe57Bm9bNeq9W0trYm6fCkNAv/74/b2NhQpVJRJpOR7/uSpFKppI2Njck2bMwuy3cioeWUfN9Xs9kc+AZOu2q12jfxand3V9LhX2HTwPO8oetisZhevXo1xtZMxiy8z2HM8me9+zmoVCpBL9va2tpMhZdUKqVisaharabl5WVlMhmlUqmZm5B+Wb4TuXrolCqVylTPHo+qWCyqWCzOzBd79y+uWTNr77M0u5/17kmqXq8HV0vF43Hl83ndv39/poaI0um0dnZ25HmearWaJGlra6tn2GzWjes7kZ6WU3AcZ6b+0jxJ9y+PafpiH/VlNItj2dJ0vs8n4bN+GFSOSqVS8n1flUplQi0aL9d1lc/nVa1Wtbe3F/S63L17d9JNG6vL8p04sz0trutqfX09dPmtra3gw+u6rvFf3GfZ/6MqlYpisZjK5fJ5Nm/iYrGYpMF/Pfi+P3N/YU3r+3ySafisn1b3MzBsXtPe3t44mzMx6+vrPfuay+WUTqd19+5dVSqVmRkmuizfiTMbWuLx+Kk+dN2b7hy9qqI7vpvP54256dBp9/+oWq0m3/d7TmTTckK3LEuWZQ39C6I7OXEWTPP7PMq0fNZPq/sZGNbtPwuTtD3PC07WR9m2rc3NzZkJbtLl+U6c2dByWoNSdaVSkeM4Ay+Lm1au66rZbPZ8afu+L8dxpuYqiwcPHqjRaPQs647zz8qQwSy8z8PwWT/8DOzs7PQs64aYWfgM2LY9dAKqZVkzN0R0Gb4TmdNyDmZtUqbneSoUCorFYqrVasEjn89P1QTNfD4fTLrrKpfLMzNEMivvcxSz9lkvFotyXbfnxJ3P55VOpwcOF0+jdDqtUqnUs8z3fdXr9akdGupe3nzcZfhOnOt0Op2x1TZlur/HUavV5Hme0um01tbWpvY/ctfy8vLQL+9p++/kuq62t7eVTCaDL+5pHhI4apbe55PM6mddOjyB5fP5YDhw2ofFBqlUKtrb25vqY+D7vgqFQjDJ2rIsPXjwQHfu3OnZ10l/JxJaAACAERgeAgAARiC0AAAAIxBaAACAEQgtAADACIQWAABgBEILAAAwAqEFAAAYgdACAJfARdxtd9gt6AFTEVqACdnY2Oi7JTZmU/eus+fNdV05jnPu2wUmhdACTIDneapUKpFOVLVaTZlM5gJbdbnrvywu4jisr69fyI8wptNpVatVelwwNQgtwASUy2UVi0V5nhe6t2V7e1uxWOyCW3Z5678szvs4VCoVra2tBb9rc96KxaI2NjYuZNvAuBFagDHzPC/4ETLLskL3tjiOo7W1tQtu3eWt/7I47+NQLpcv9IcXLcuSbdsME2EqEFqAMTt6ktrc3JTneSeeULo/FZ9KpcbRxEtX/2Vx3sfBcRwlEolz2dYomUxG5XL5wusBLhq/8gyMUfdn34/+lPvy8rJs29be3l5fecdxVK/X5bqudnd3e8LORQ0nXKb6L4uLOg4bGxtaW1tTOp0Olrmuq/X1dXmep6+//jqYj7KzsyNJfXNfarWams2mYrGYms2m9vb2lMlk+oLV8vKyWq3WqdsKXAodAGOTy+X6lhWLxY6kTr1eH/q6bDbbSafTF9m0kSZd/2Vx3sfBtu1Oo9EYuE5SJ5vNdlqtVrDMsqxOtVoNnjcajU42m+15XbFYHPh/aVRdgCkYHgLGxPd93b59u295mLktUeZRdP96D/sIMxH4pPpLpZLm5ua0sbEh13WD19y9e1d37txRpVLpKV+pVDQ3N9e3z6VSScvLyyqVSkPryGQyfdsbl/Oez9JsNmXb9sB1lmUFj65EIhH0uEiHvTLHrww62mtzlG3bXEUE8006NQGzIpfL9fzVfHydpM7e3l7fular1ZE0sb+Sw9TfaDQ6kvr2L5VKdYrF4sDy2Wy2Y1lW37pB5UfVMS4X8T6M+go+3qvS6XQ66XS6p7eu1Wp1LMvq2LbdyeVyI3vr0ul0p1wun73RwATR0wKMQfdup8PmP2xubkrSwN6W3d3d4AqQSQhTv+u6sm27b/8cxxk4adV1XRWLRcVisZ6eHs/zhk5yHVbHqDbdvXs39KPbQzTMeb8P53EHXMuy9N133ymVSqlWq2ltbU137ty5kLvrApfBtUk3AJgFlUolCCaDWJalXC6nUqkk13UVj8eDdfV6vecKE9/3R564NzY2Ig0DbGxsDB1SCFt/vV7vCxvdEHB0X7qazaYsy9LGxoYKhUJQv+M4Qy//HVTHKPF4fODk5tOK+j6c5DwmMnf/r3SvDPI8Lzimxyfs+r4/seALnBdCC3DBfN/Xq1evTjxJbW5uqlQqKZ/Pq16vB8td1+2ZR3H86qPjzvvS1jD1O46jeDzeM9ekXq8PDCxHZbNZ5fN5eZ534gnVcZwLudV9WCcdh1qt1rMfo4Jgl2VZZwo/nuep2WwGYc62bZXL5YE3k+teYQSYjNACXLBCoaDbt2+HmvAaj8flOE5fb0v3RFipVC70RmTDjKrf9315nqd6vd4TPML0jFiWpXQ6rWKxODJsdeuY9H1ihh2HSqWier2uarUaXLIcJrR0J8eeFO66uveJOapYLPYcF8/zBk4WDhMMgcuO0AJcsEqlEnmOQaFQULValXQ4z6X770QiMfb7o5xU/7C5Ho7jaGtrq297x+e5dK92yufzQ0+qo+oYV5AZdhx839fGxkZwD5R4PK6vv/461DYfPnyo3d3dntDiuq7K5bJ831ehUFCz2dSDBw9UKBS0u7ur3d1dSf+/Ry2TyahUKvW053hPWHdoaJburYPpxM3lAJxJPp+X67o9Q1q+7wc3Mzt+ohzUW3Pnzh3Ztq1qtTrwxDqoDunwMuhRQ2XjUKvVVCgUTjV/xvO8njB0Ubo3oJtELx1wnrh6CMCZHJ/rIR32gAz7y35Qr1N3XsuwnoBBdeTz+YkPF0mHQ1zH54qE/Z2f7vDQRV/tUygU9ODBgwutAxgHQguAU/F9X6VSSY7jqNFoBFcsdXseuv/u6l7ZUi6X+07qDx48GDgH5HgdlUpFpVIpuCle2LkgFymVSgWTkGu1miqVSqS5I1tbW8Hxugiu6yqVSjE0hKnA8BAATFipVArCz3nyfV/r6+sXPvwEjAs9LQAwYblcTo7jnPswUaVSGTgZGjAVPS0AAMAI9LQAAAAjEFoAAIARCC0AAMAIhBYAAGAEQgsAADACoQUAABiB0AIAAIxAaAEAAEYgtAAAACP8PyURbTp8bkaGAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiYAAAG5CAYAAAC6Fv9NAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAycklEQVR4nO3dT2wbZ+Lm+cdxS2rJbrtEuQ1s3GogZcOHRRYwKCrAAgMEC1O3OfWSFuYwR5KHxuwECUBCpx98WBDkINgAgz6Quu2eZHJzmsuCFTSCnVPLqgmwgz1kwurfrtbu3rRFlpXYiux2uAcuq0XxX1F/yCL5/QACxKqXb70kRfLR+7711pVGo9EQAABAALw37gYAAAC0EEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBg/GLcDRjGzz//rOfPn+tXv/qVrly5Mu7mAAAAHxqNhn744Qe9//77eu+9/n0iExVMnj9/rtXV1XE3AwAAnMH+/r5+85vf9C0zUcHkV7/6laTmA7tx48aYWwMAAPw4PDzU6uqq9z3ez0QFk9bwzY0bNwgmAABMGD/TMJj8CgAAAoNgAgAAAmOihnIAALPn3bt3evv27bibgQHm5+cHnnHjB8EEABBIjUZDf/3rX+W67ribAh/ee+89ffDBB5qfnz9XPQQTAEAgtULJ7du3tbS0xPpVAdZaZ+wvf/mLfvvb357rtSKYAAAC5927d14oWVlZGXdz4MOvf/1rPX/+XH//+981Nzd35nqY/AoACJzWnJKlpaUxtwR+tYZw3r17d656CCYAgMBi+GZyXNRrxVAOAGCiPHOPVH/1ZmTHW742rzvG4siON+sIJgCAifHMPVL086919PZ8wwXDWJy7Kuuzj32FE9u2VSgUVCwWlU6ntbm5KdM09eTJE2UyGUUiEcXjcUUiEe3s7Cifz6tQKCiZTEqSHMdRLpeTZVle+VZ9yWRSGxsbisViHcd1HEflclmGYUiSTNOU4zhevZPkSqPRaIy7EX4dHh7q5s2bevnyJUvSA8AU++mnn/TnP/9ZH3zwgX75y1962//zs5f6l//+P+qLzQe6d/v6pbfju+9/1Cc73+g//Jt/oQ/v3PR1H8dxdPfuXdXrdS8oSNLa2po2NzeVTqclSa7rKpPJqFgsdpTN5/NeuV71nbS2tqa9vb22+x8cHCiXyw33gM+h12smDff9TY8JAGDi3Lt93XdQGLVQKOS7bDweV61WUyKRUKlU8rafDCCD6nMcp2NbOp1WPp/33Y4gYfIrAABjtL29LcuyVC6Xz3T/1rBNsVhs2z6JwzgSPSZA4PmZ6MfkPCB4isViW89HrxVsDcNQLpdTIpHoOn/Ej+3tbcXjcaVSKUWjUWUyGUWj0TPVNW4EEyDA/E70G2ZyHoDRSCaTbcGkUCj0LVsqlZRKpfqW6yUWi6larcqyLFUqFW1sbKhUKp056IwTwQQIsPqrNzp6+67vRL/W5Lz6qzcEE2CCFQoF3b17V6lUylf5crmsWCwm13VlGIZM01QymVQymVSxWFQ2m53IYMIcE2ACtCb6dfsZxZkJAC6faZrK5XKKx+O+yrcmvTqOI9u22/Y9evRoYi9+SDABAOAC1Wq1M5dLp9MyDKMtVHQr57quKpWKdzuTybTttyxrIntLJIZyAAAT6LvvfwzkcVoLrElSNpttW2DNcRzt7OzINE2ZpqlMJqNaraatra22ELG9va2nT592re/u3buqVqsql8syTdO7TzweVz6f9+a0VKvVka5hcpFYYA0IsNZiUv0Wd/JTBpg0vRbrCvrKr7OMBdYAADPnjrEo67OPuVbOFCOYAAAmyh1jkaAwxZj8CgAAAoMeE8wGd196fTC43NKKZKxefnsAAF0RTDD93H3pDx9Jb18PLju3JP3+T4QTABgTggmm3+uDZij53bZ0637vci++lb5MNMsTTABgLAgmmB237kvvPxh3KwAAfRBMAACTxe+csYvC3LORIpgAACbHMHPGLgpzz0aKYAIAmBx+54xdlCHnntm2rZ2dHeXzeRUKBSWTSUnNC+3lcjlZlqVMJuNtv0xra2sdy91PAoIJAGDyBHTOWDgclmmacl1XqVRKjx49kmEYMk1ThUJB+Xz+QkKJ67redXF6yeVyikQi5z7WqBFMgCkx6GJjLKsNjE48HletVlMikVCpVPK2DwoTfjiOI8uyBgacaDR67mONA8EEmHDL1+a1OHdVn+x807ccFyIDRmt7e1sffPCByuXyhQ6n5HI5ra2tXVh9QTNUMHEcx7vUsuM4SiaTvtJfJpPR1taWV/as9QDo5OeiZt99/6M+2flG9VdvCCbAiBiGoVwup0Qi0TWYWJYl27ZlmqZ2d3fb5qBsbm4qnU6rXC4rk8moUCgoGo3Ksiw9ffpUtVpNUrNXxHEcpVIpZTIZSVKhUND29rYSiYRSqZSSyaRXbyqV8r57K5VKW2+ObduyLMvbH4vFvGGpJ0+eeL+32npZhgom8Xhce3t7kprh4nQXVTe2bSufz2tra+tc9QDojYuaAcGUTCZVKpWUSqVUKBS87Y7jKJPJeN+FtVpN+Xxe6XRam5ubXrlYLKadnR3vdjQaVTQa1d27d72hHNM0FY1Gtbe3p0KhoFAopHA43FZP634nw0ipVJJt2wqHw157KpWKd5+1tTV99dVXKhaLCofD3tBQKxRdFt/BxHGcttumacqyLF/3M03z3PUAADCJCoWC7t69q1Qq1bYtFAq1ff/t7u6e+RiGYWhlZUWSeg4brayseGVa92mFjEKhoHA43FbeNE09efJEsVhMa2trMk1Tm5ubl35Gke+rC1uWpVAo1LYtFArJtu2e9+k2rnaWegAAmFSmaSqXyykej7dtb/VCRKNRr2dlWK7rth3nMoRCIdXrdW1vb+vg4KDjcVw03z0mJx/8Sb26dHqdyjRMPcfHxzo+PvZuHx4eDmwnAGAGvPg2sMfp9n2WTqe1s7PjfQdubm4qkUi0lbEsS9FoVIZhqFqttm0/OSxz+j6tDoBux+31nXtat/bYtq3t7W1ls1mlUimFw2GFw+HgBJNeej3oJ0+eDNXd062ebDarx48fn7FlAICps7TSXIn1y8Tgshdlbql5XB9s21Ymk1GtVutY3Gx7e1tPnz6V1OwtyeVyymQyWl9fl/SP03sfPXqkTCbjDfNEo1FvqMU0TaVSKeVyORWLRW9C7MmJtNFo1FvoLRQKKRaLyXVdb65Ka8KsbdsqFAoyTdNrTz6f9ybjlkolb4ioNdpRq9V6hqSL4juYnByLaqnVal17RSzL0qNHj85dz9bWlj799FPv9uHhoVZXWRIYAGaWsdpcHj6g18oJh8NtE0hP7zs5j6M1jHOaYRhtE2VPl2kt1nbydmsS7cljnd528nY4HO6YanGyPSf3pdPpro/nsvgOJq3EdlqvVeWePHni/e44jrLZrDY3N4eqZ2FhQQsLC36bCACYBcYq162ZYr6DyelJNY7jKBKJeD0dtm17y+6eTnepVMo7d/q00/UAAIDZ5fusHKl5znMmk1G5XFahUGibQZzNZlUul9vKu66rfD4vqblSXevMm371AACA2TXU5NfWKU9S53nS3cKFYRhKp9Md41P96gEAALNrqB4TAABG6eeffx53E+BTo9G4kHq4iB8AIHDm5+f13nvv6fnz5/r1r3+t+fl5XblyZdzNQg+NRkN/+9vfdOXKFc3NzZ2rLoIJACBw3nvvPX3wwQf6y1/+oufPn4+7OfDhypUr+s1vfqOrV6+eqx6CCQAgkObn5/Xb3/5Wf//73/Xu3btxNwcDzM3NnTuUSAQTAECAtYYGzjs8gMnB5FcAABAYBBMAABAYDOUApw26mugQ180AAAyHYAK0+L1q6dxS8yJihBMAuHAEE6DFz1VLX3zbDC6vDwgmAHAJCCbASVy1FADGismvAAAgMAgmAAAgMAgmAAAgMJhjgsnn7g+esAoAmAgEE0w2d1/6w0fS29f9y80tNU8HBgAEGsEEk+31QTOU/G5bunW/dzkWRQOAiUAwwXS4dV96/8G4WwEAOCcmvwIAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMBgSXpgTJ65R6q/etO3zHff/zii1gBAMBBMgDF45h4p+vnXOnr7bmDZxbmrWr42P4JWAcD4DRVMHMdRuVyWaZpyHEfJZFKGYXQta1mWJMl1Xe3u7mpzc1PhcFiSZNu2JCkcDstxHLmu6+0DZkH91RsdvX2nLzYf6N7t633LLl+b1x1jcUQtA4DxGiqYxONx7e3tSWqGlEQioVKp1LPsV199pWg0qlqtpng8rmq1KkkqFAoqFouSpGg02rMOYNrdu31dH965Oe5mAEBg+A4mjuO03TZN0+sV6aZUKrX1gpzsWVlbW1O9Xu/YDgAAZpvvs3Isy1IoFGrbFgqFvGGZ06LRqPd7qVRSKpVq228YBqEEAAC08d1j4rpu1+21Wq3nfWzb1s7OjjY2NpRMJtvqKpfLkqTd3V2lUimZptlx/+PjYx0fH3u3Dw8P/TYXAABMoHOfldMrsEjNya2maSqTyahcLisWi0lS26RZ0zS1sbHhzT85KZvN6vHjx+dtIgAAmBC+h3IMw+joHanVagOHYwzDUDweVzwe90LMyfkqrTN8Ts9hkaStrS29fPnS+9nf3/fbXAAAMIF8B5OTc0ZOikQiHdssy9Ly8rJ3uzVM4ziObNvWw4cPO+5zev6KJC0sLOjGjRttPwAAYHr5Hso5PQfEcRxFIhGvx8S2bRmGIdM0FQqF2oJMa184HJbrusrlct4+y7IUi8WYCIvJ8uLbvrv/+vdrenH1ds/9rOgKAN0NNcekVCopk8lofX1du7u7beuPZLNZra+vK51OKxwOa3Nz01urpFKpeOufGIahSCSifD4vwzBUrVZZxwSTY2lFmluSvkz0LXajsaDfHf87PdetnmVY0RUAOl1pNBqNcTfCr8PDQ928eVMvX75kWAdNz7+Rih9Lya+l9x+M5pjuvvT6oOfu/f/yjVb/+G/1x4/L+vX9j3qWG+WKrv/52Uv9y3//H/Uf/s2/YEE3ACM3zPc318oBhmWsNn96OP5bc5hmNbSoe4QAABiK78mvAAAAl41gAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAuMX424A0Je7L70+6L3/xbejawsA4NIRTBBc7r70h4+kt6/7l5tbkpZWRtMmAMClIpgguF4fNEPJ77alW/d7l1takYzV0bULAHBpCCYIvlv3pfcfjORQz9wj1V+96Vtm+dq87hiLI2kPAMwaggnw/3vmHin6+dc6evuub7nFuauyPvuYcAIAl2CoYOI4jsrlskzTlOM4SiaTMgyja1nLsiRJrutqd3dXm5ubCofDQ9cDjEr91RsdvX2nLzYf6N7t613LfPf9j/pk5xvVX70hmADAJRgqmMTjce3t7UlqhotEIqFSqdSz7FdffaVoNKparaZ4PK5qtTp0PcCo3bt9XR/euTnuZgDATPIdTBzHabttmqbXK9JNqVTyekgkeT0iw9YDTKoF9zvpefeeF0lM2gWALnwHE8uyFAqF2raFQiHZtt0WQFqi0aj3e6lUUiqVOlM9wKR598uQXjcWtPrHfyv9sU/BuSXp938inADACb6Dieu6XbfXarWe97FtWzs7O9rY2FAymRy6nuPjYx0fH3u3Dw8P/TYXGJu31+8oevzv9D//q7u69+sePSYvvpW+TDRPiSaYAIDn3Gfl9AoakhQOh2WapjKZjMrlsmKx2FD1ZLNZPX78+LxNBEbuuW7pp1v/jfQ+c1UAYBi+r5VjGEZHr0atVht4No1hGIrH44rH43Jdd6h6tra29PLlS+9nf3/fb3MBAMAE8h1MTs4ZOSkSiXRssyxLy8vL3m3TNCU1J74OU8/CwoJu3LjR9gMAAKaX76GcVrhocRxHkUjE6+mwbVuGYcg0TYVCobYA0trXbXLr6XoAAMDsGmqOSalUUiaT0fr6unZ3d9vWHslms1pfX1c6nVY4HNbm5qaKxaIkqVKpeOuWDKoHAADMrqGCiWmayuVyktQxkfV0uDi5v3VGjp96AADA7PI9xwQAAOCyEUwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBg/GKYwo7jqFwuyzRNOY6jZDIpwzC6lrVtW5ZlSZJ2d3e1vb3tlbVtW5IUDoflOI5c11U4HD77owAAAFNhqGASj8e1t7cnqRlSEomESqVS17KWZSmdTkuS8vm8Hj586N23UCioWCxKkqLRaM86AADAbPE9lOM4Tttt0zS9HpHTbNtWNpv1bsdiMdm27dWxtramer2uer2uSqXSs9cFAADMFt/BxLIshUKhtm2hUMgbljkpHA5re3vbu+26rle+xTAMAgkAAGjjeyinFS5Oq9VqXbfHYjHv952dHUWjUS+IuK6rcrksqTn/JJVKyTTNjjqOj491fHzs3T48PPTbXAAAMIGGmmPSTa/AcnJ/uVz25pdIaps0a5qmNjY2VK1WO+6bzWb1+PHj8zYRAABMCN/BxDCMjt6RWq02cDgmk8l0zCNxHMc7C6d1ho/jOB29JltbW/r000+924eHh1pdXfXbZCD4Xnw7uMzSimTwdw9gNvgOJtFoVIVCoWN7JBLpeZ98Pq9MJiPTNL2eFcdx9PDhQ9Xr9bayp+evSNLCwoIWFhb8NhEYme++//FM+zxLK9LckvRlYnDZuSXp938inACYCb6DyeneDMdxFIlE2tYmMQzDK1culxUOh71Q8uTJEyWTSZmmqVwu59VjWZZisRgTYTERlq/Na3Huqj7Z+aZvucW5q1q+Nt+7gLHaDBuvD/of8MW3zfDy+oBgAmAmDDXHpFQqKZPJaH19Xbu7u23rj2SzWa2vryudTstxHMXj8bb7GobhzS2JRCLK5/MyDEPVapV1TDAx7hiLsj77WPVXb/qWW742rzvGYv/KjFXCBgCcMlQwOdnbcfKsG0lt4cI0TTUajZ71hMNhVnrFxLpjLA4OHQCAM+FaOQAAIDDOfbowMAmeuUcDh198TVoFAFwqggmm3jP3SNHPv9bR23cDyw6ctAoAuFQEE0y9+qs3Onr7Tl9sPtC929f7lvU1aRUAcGkIJpgZ925f14d3bo67GQCAPpj8CgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoPThTE+7n7/q+u++HZ0bQEABALBBOPh7kt/+Eh6+7p/ubklaWllNG0CAIwdwQTj8fqgGUp+ty3dut+73NKKZKyOrl0AgLEimGC8bt2X3n8w7lYAAAKCya8AACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwCCYAACAwuLowJt4z90j1V2967v/u+x9H2BoAwHkQTDDRnrlHin7+tY7evutbbnHuqpavzY+oVQCAsxoqmDiOo3K5LNM05TiOksmkDMPoWta2bVmWJUna3d3V9va2V3aYeoB+6q/e6OjtO32x+UD3bl/vWW752rzuGIsjbBkA4CyGCibxeFx7e3uSmuEikUioVCp1LWtZltLptCQpn8/r4cOH3n2HqQfw497t6/rwzs1xNwMAcE6+J786jtN22zRNr0fkNNu2lc1mvduxWEy2bctxnKHqAQAAs8V3MLEsS6FQqG1bKBSSbdsdZcPhsLa3t73brut65YepBwAAzBbfQzmtcHFarVbruj0Wi3m/7+zsKBqNyjCMoeo5Pj7W8fGxd/vw8NBvcwEAwAQ69zomvYLGyf3lcnngHJJu9WSzWd28edP7WV1dPUdLAQBA0PkOJoZhdPRq1Gq1gWfTZDIZVSoVr9ww9Wxtbenly5fez/7+vt/mAgCACeQ7mESj0a7bI5FIz/vk83llMhmZpinXdeW67lD1LCws6MaNG20/AABgevkOJqZptt12HEeRSMTr6WidddNSLpcVDoe9UPLkyRMZhjGwHgAAMLuGWsekVCopk8lofX1du7u7bfNGstms1tfXlU6n5TiO4vF4230Nw1AymRxYDwAAmF1DBRPTNJXL5SS1n3UjqS1cmKapRqNxpnoAAMDs4urCAAAgMAgmAAAgMAgmAAAgMAgmAAAgMAgmAAAgMIY6KwcAnrlHqr96M7Dc8rV53TEWR9AiANOEYALAt2fukaKff62jt+8Gll2cuyrrs48JJwCGQjAB4Fv91RsdvX2nLzYf6N7t6z3Lfff9j/pk5xvVX70hmAAYCsEEl8Pdl14f9N7/4tvRtQUX7t7t6/rwzs1xNwPAFCKY4OK5+9IfPpLevu5fbm5JWloZTZsAABOBYIKL9/qgGUp+ty3dut+73NKKZKyOrl2TbFAPE88lgClBMMHluXVfev/BuFsx2ZZWmj1LXyb6l5tbkn7/J8IJgIlHMAGCzFhtBo5B83W+TDTLEEwATDiCCRB0xiqBA8DMYOVXAAAQGAQTAAAQGAzlADPku+9/7LufZeQBjBvBBJgBy9fmtTh3VZ/sfNO3HMvIAxg3ggkwA+4Yi7I++7jvxfdYRh5AEBBMgBlxx1gkcAAIPCa/AgCAwKDHBECbfhNkB02eBYDzIpgAkDTcBNnla/OjaRSAmUMwASDJ3wRZiVOKAVwuggkADxNkAYwbk18BAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgDBVMHMdRPp9XuVxWPp+X67p9y9u2rbW1ta7bbdv26mz9DgAAZttQpwvH43Ht7e1JagaKRCKhUqnUtWy5XJZpml1DR6FQULFYlCRFo9GedQAAgNniO5g4jtN22zRNWZbVs3wsFuu5b21tTfV6XZJkGIbfJgAAgCnneyjHsiyFQqG2baFQ6MzDMIZhEEoAAEAb3z0mveaT1Gq1oQ/quq7K5bIkaXd3V6lUSqZpdpQ7Pj7W8fGxd/vw8HDoYwEAgMlx7iXpB02A7SaZTHq9JaZpamNjQ9VqtaNcNpvV48ePz9lCAAAwKXwP5RiG0dE7UqvVzjQcc3K+immachynYw6LJG1tbenly5fez/7+/tDHAgAAk8N3MIlGo123RyKRoQ5o27YePnzYsf30/BVJWlhY0I0bN9p+AADA9PIdTE7PAXEcR5FIxOsxsW27a6+H1D7cY5qmcrmcd9uyLMViMSbCAgCA4eaYlEolZTIZra+va3d3t239kWw2q/X1daXTaUnNwFGpVNr2tQJIJBJRPp+XYRiqVqusYwIAACQNGUxO9nacXqfkdLiIRqOKRqNtvSMt4XBY4XB42LYCAIApx7VyAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYBBMAABAYPxi3A0AcEFefNt//9KKZKyOpi0AcEYEE2DSLa1Ic0vSl4n+5eaWpN//iXACINAIJsCkM1abgeP1Qe8yL75tBpfXBwQTAIFGMAGmgbFK4AAwFZj8CgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoOzcjAcd7//aanS4IW+AADogWAC/9x96Q8fSW9fDy47t9Rc+AsAgCEQTODf64NmKPndtnTrfv+yLH8OADgDggmGd+u+9P6DcbcCADCFmPwKAAACg2ACAAACg2ACAAACY6g5Jo7jqFwuyzRNOY6jZDIpwzB6lrdtW4lEQnt7e+eqBwAAzIahgkk8HvdChuM4SiQSKpVKXcu2godt2+eqBwAAzA7fwcRxnLbbpmnKsqye5WOx2IXUAwAAZofvYGJZlkKhUNu2UCgk27YVDod9H/Ci6gEQfN99/2Pf/cvX5nXHWBxRawBMAt/BxHXdrttrtdpQBxymnuPjYx0fH3u3Dw8PhzoWgPFYvjavxbmr+mTnm77lFueuyvrsY8IJAM+5F1jrFTQuop5sNqvHjx9fSP0ARueOsSjrs49Vf/WmZ5nvvv9Rn+x8o/qrNwQTAB7fwcQwjI5ejVqtNvTZNMPUs7W1pU8//dS7fXh4qNVVljkHJsEdY5HAAWBovtcxiUajXbdHIpGhDjhMPQsLC7px40bbDwAAmF6+g4lpmm23HcdRJBLxejps2+4446bl5DDNoHoAAMDsGmqOSalUUiaT0fr6unZ3d9vWHslms1pfX1c6nZbUPPumUqm07WudQtyvHgAAMLuGCiamaSqXy0nqXKfkdLiIRqOKRqNeeb/1AACA2cW1cgAAQGAQTAAAQGAQTAAAQGAQTAAAQGCce+VXADgPrqcD4CSCCYCx4Ho6ALohmCDQnrlHA6+3gsnE9XQAdEMwQWA9c48U/fxrHb1917fc4txVLV+bH1GrcJG4ng6A0wgmCKz6qzc6evtOX2w+0L3b13uWYw4CAEwPggkC797t6/rwzs1xNwMAMAKcLgwAAAKDYAIAAAKDYAIAAAKDYAIAAAKDYAIAAAKDYAIAAAKDYAIAAAKDYAIAAAKDYAIAAAKDlV+BWfLi2/77l1YkY3U0bQGALggmwCxYWpHmlqQvE/3LzS1Jv/8T4QTA2BBMgFlgrDYDx+uD3mVefNsMLq8PCCYAxoZgAswKY5XAASDwmPwKAAACg2ACAAACg2ACAAACgzkm+Ad3f/DkSAAALhHBBE3uvvSHj6S3r/uXm1tqnnoKAMAlIJig6fVBM5T8blu6db93ORbgAgBcIoIJ2t26L73/YNytAADMqKGCieM4KpfLMk1TjuMomUzKMIyhy9q2LUkKh8NyHEeu6yocDp/rgQAAgMk3VDCJx+Pa29uT1AweiURCpVJp6LKFQkHFYlGSFI1Ge9YBAABmi+9g4jhO223TNGVZ1pnKrq2tqV6vS1LPHhcAADB7fAcTy7IUCoXatoVCIdm23TEM46csgQTP3CPVX73puf+7738cYWsAAEHgO5i4rtt1e61WG7qs67oql8uSpN3dXaVSKZmm2VH++PhYx8fH3u3Dw0O/zUXAPXOPFP38ax29fde33OLcVS1fmx9RqwAA43bus3J6hZB+ZU9OhDVNUxsbG6pWqx3ls9msHj9+fN4mQgrc4mn1V2909Padvth8oHu3r/cst3xtXneMxRG2DAAwTr6DiWEYHb0jtVqt65DMoLKO43hDOq2zdhzH6eg12dra0qeffurdPjw81Ooqa2gMLcCLp927fV0f3rk50mMCAILLdzCJRqMqFAod2yORyFBlbdvWw4cPvcmvLafnpEjSwsKCFhYW/DYRvbB4GgBgQvgOJqd7MxzHUSQSaVubxDAMmabZt6xpmsrlct4+y7IUi8WYDDsKLJ4GAAi4oeaYlEolZTIZra+va3d3t239kWw2q/X1daXT6b5lDcNQJBJRPp+XYRiqVqusYwIEyaD5RvSsAbhEQwWTk70dsVisbd/pcNGvbDgcZqVXIGiWVprzjL5M9C83tyT9/k+EEwCXgmvlAGgyVpuBY9DZW18mmmUIJgAuAcEEwD8YqwQOAGP13rgbAAAA0EIwAQAAgcFQDi4F18EBAJwFwQQXjuvgAADOimCCC8d1cAAAZ0UwwVAGDdFI/xim4To4AIBhEUzgm98hGolhGgDA2RBM4JvfIRqJYRoAwNkQTDA0hmgAAJeFYAIg8AadXk4PHTA9CCYAAmv52rwW567qk51v+pZbnLsq67OPCSfAFCCYAAisO8airM8+HrhY3yc736j+6g3BBJgCBBMAgXbHWCRwADOEa+UAAIDAIJgAAIDAYChnGrj70uuD3vtffDu6tmA2+PmbWlqRjNXLbwuAqUIwmXTuvvSHj6S3r/uXm1tqflEA57G00vxb+jIxuOzckvT7PxFOAAyFYDLpXh80Q8nvtqVb93uX479XXARjtRk2+vXQSc0elS8TzXL83QEYAsFkWty6L73/YNytwCwwVgkbAC4Nk18BAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgEEwAAEBgcLpw0LGqKybZoL/PAK6v88w96ns1Y0lavjbPhQWBS0IwCTJWdcWk8rtCbMBWh33mHin6+dc6evuub7nFuauyPvuYcAJcAoJJkLGqKyaVnxViA7g6bP3VGx29facvNh/o3u3rXct89/2P+mTnG+3+uaZ6jzISvSrAWQ0VTBzHUblclmmachxHyWRShmEMXXaYeiBWdcVkmuAVYu/dvq4P79zsum/52rwW567qk51v+tZBrwpwNkMFk3g8rr29PUnNcJFIJFQqlYYuO0w9AODHd9//2Hf/RfVg3DEWZX32cd95KK1elfqrNwQTYEi+g4njOG23TdOUZVlDlx2mnok1aMKqX0xsxSw459/57VfHWpv7v7T95M99yx39wtD/8tl/f2HhhMABXA7fwcSyLIVCobZtoVBItm0rHA77Lvv06VPf9QTSoNDx+oW0868HT1j1i4mtmFZ+J8gOcFvS/3pV0tX+5V43FrT//4Yl478+1/EAXC7fwcR13a7ba7XaUGWHqef4+FjHx8fe7ZcvX0qSDg8P+zf2jA7+8n/LffGs5/5f/FTX6h//B7337qe+9fx89Zfa/+8K+vsvl8/dpp9/uay3f5X01/1z13Vezt9e6efj1/rxh0MdHl4Zd3Mw6d67Kf3rr6Sjzvf+Rft/qv+HfvO/Z/SX//S/6fWL3u+lZ/Uj3X/zX/T2n3+hwx+vnfl4b1+80v03/6f++T+91dt/pmcFk8W4dUcr/9VvL7TO1vd2o9EYWPbcZ+X0ChrDlu22L5vN6vHjxx3bV1eDPqHuB+l//FfjbsSl+W+/GHcLgLP61Fepyv90MUe7qHqAafHDDz/o5s3uE8tbfAcTwzA6ejVqtVrXs2n6lR2mnq2tLX366T8+SH7++WfVajWtrKzoypWL/Y/98PBQq6ur2t/f140bNy607kkw649f4jmY9ccv8Rzw+Gf78UuX9xw0Gg398MMPev/99weW9R1MotGoCoVCx/ZIJDJUWdM0fdezsLCghYWFtm2XfVrxjRs3ZvYPUuLxSzwHs/74JZ4DHv9sP37pcp6DQT0lLb6DiWmabbcdx1EkEvGCgm3bMgxDpmn2LXs6WJyuBwAAzK6h5piUSiVlMhmtr69rd3e3be2RbDar9fV1pdPpgWX77QMAADOsgUaj0Wj89NNPjX/6p39q/PTTT+NuyljM+uNvNHgOZv3xNxo8Bzz+2X78jUYwnoMrjYaPc3cAAABG4L1xNwAAAKCFYAIAAAKDYAIAAAKDYAIAAALj3EvSz5KNjQ1VKpVxN2NkLMtSpVKR67pyHEfxeFzJZHLczbpwjuMol8vp7t27kpqL+E3j4+xlVl7nYczae11qvg8KhYJWVlZ0cHCg9fV1xWKxcTdrZGzb9q50f3BwoJWVFW/5i2lULBZVrVaVy+U69o39M3Fs5wNNmFwu15ilp6tSqTRyuZx3u16vN0zTbCSTyTG26uJVq9WGYRiNer3ubUun022PfZrNyus8jFl7rzcazb+DWCzm3a7X641wODzGFo1WtVrteM/v7e21PSfToFqtNpLJZCOZTDYMw2ik0+muZcb9mThb774zqlarjVgsNlMfVt3ekIVCoSGpUa1Wx9Ciy5FMJjvenPV6fWZe61l5nf2axfd6vV7v+CKqVCoNwzDG16gRSyaTXf/eo9HoGFozGuFwuGswCcJnInNMfCiXy9rc3Bx3M0aqXC4rk8m0bWtdz6jV3TkNnjx54nVXtrQujzBNj7OXWXmd/ZrF93o2m+24LEg0GlW9Xh9fo0asVqt1HdI4fcHZWRCEz0SCyQDlcnmmxllbYrFYxx/ntHFdV67rdlzbSWq+EW3bHkOrRmsWXme/ZvW9Xi6XtbGxIan5xTMLf/enpVIpFYtFxeNxua4rScrn80qlUuNt2IgF5TORYNKH67qq1WpdX6RpVyqVOiY7PX36VFLzv6lp4DhOz32hUEgHBwcjbM14zMLr7Mcsv9db74Nisej1lm1sbMxUQIlGo8rlciqXy1peXlY8Hlc0Gp25SeBB+UzkrJw+isXiVM/KHlYul1Mul5uZD+/Wf06zZtZeZ2l23+utL6JKpeKdhRQOh5XJZPTw4cOZGs6JxWLa3d2V4zgql8uSpO3t7bYhrlk3qs9Eekx6sCxrpv5jHKT1H8Q0fXj3+8CZxbFlaTpf50F4rzfDyEnRaFSu66pYLI6pRaNl27YymYxKpZL29va83pO1tbVxN22kgvKZONU9JrZtK5FI+C6/vb3tvUFt2574D+fzPP6TisWiQqGQCoXCRTZv7EKhkKTu/wW4rjtz/ylN6+s8yDS818+q9R7oNc9ob29vlM0Zm0Qi0fZY0+m0YrGY1tbWVCwWZ2ZIJyifiVMdTMLh8JneWK2FZ06erdAab81kMhOz8M5ZH/9J5XJZruu2fVlNy5e2YRgyDKPnfwKtCYGzYJpf536m5b1+Vq33QK8u+lmYGO04jveFfJJpmtra2pqZcCYF5zNxqoPJWXVLx8ViUZZldT2lbFrZtq1ardb2wey6rizLmpqzFx49eqRqtdq2rTXuPivd+7PwOvfCe735Htjd3W3b1goqs/AeME2z56RPwzBmbjgnCJ+JzDHxadYmQjqOo2w2q1AopHK57P1kMpmpmhSZyWS8iW4thUJhZoYzZuV1HsasvddzuZxs2277cs5kMorFYl2HdqdRLBZTPp9v2+a6riqVytQO47RODT4tCJ+JVxqNRmNkR5tAretHlMtlOY6jWCymjY2Nqf1jbVleXu75AT1tfzK2bWtnZ0fr6+veh/M0d9+fNEuv8yCz+l6Xml9SmUzGG7qb9iGsborFovb29qb6OXBdV9ls1pvYbBiGHj16pLt377Y91nF/JhJMAABAYDCUAwAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAAAoNgAgAjchmryvZaTh2YVAQT4BKlUqmO5Z0xm1qrq14027ZlWdaF1wuMC8EEuCSO46hYLA71ZVQulxWPxy+xVcE+flBcxvOQSCQu5cKAsVhMpVKJnhNMDYIJcEkKhYJyuZwcx/Hda7Kzs9P1EuyjMu7jB8VFPw/FYlEbGxvedVguWi6XUyqVupS6gVEjmACXwHEc78JYhmH47jWxLEsbGxuX3LrgHj8oLvp5KBQKl3oxQMMwZJomQzqYCgQT4BKc/CLa2tqS4zgDvzRalyGPRqOjaGLgjh8UF/08WJalSCRyIXX1E4/HR3ppeuCycHVh4IK1Lil+8jLhy8vLMk1Te3t7HeUty1KlUpFt23r69GlboLmsrv8gHT8oLut5SKVS2tjYUCwW87bZtq1EIiHHcfTVV19580N2d3clqWMuSrlcVq1WUygUUq1W097enuLxeEd4Wl5eVr1eP3NbgUBoALhQ6XS6Y1sul2tIalQqlZ73SyaTjVgsdplN62vcxw+Ki34eTNNsVKvVrvskNZLJZKNer3vbDMNolEol73a1Wm0kk8m2++Vyua5/S/2OBUwKhnKAC+S6rlZWVjq2+5lrMsy8htZ/4X5//Ey+HXT8fD6vK1euKJVKybZt7z5ra2u6e/euisViW/lisagrV650POZ8Pq/l5WXl8/mex4jH4x31jcpFzy+p1WoyTbPrPsMwvJ+WSCTi9ZxIzd6V02fcnOx9Ock0Tc7OweQbdzICpkk6nW777/f0PkmNvb29jn31er0haWz/7fo5frVabUjqeHzRaLSRy+W6lk8mkw3DMDr2dSvf7xijchmvQ7+P2dO9I41GoxGLxdp63er1esMwjIZpmo10Ot231y0WizUKhcL5Gw2MET0mwAVprerZaz7C1taWJHXtNXn69Kl3ZsU4+Dm+bdsyTbPj8VmW1XWiqG3byuVyCoVCbT02juP0nFja6xj92rS2tub7p9XT08tFvw4XsdKrYRj685//rGg0qnK5rI2NDd29e/dSVpEFguAX424AMC2KxaIXProxDEPpdFr5fF62bSscDnv7KpVK25kbruv2/XJOpVJDddmnUqme3f9+j1+pVDoCReuL/uRjaanVajIMQ6lUStls1ju+ZVk9T53tdox+wuFw1wnFZzXs6zDIRUwebv2ttM64cRzHe05PT5J1XXds4Ra4KAQT4AK4rquDg4OBX0RbW1vK5/PKZDKqVCredtu22+Y1nD6r57SLPi3Uz/Ety1I4HG6b+1GpVLqGkpOSyaQymYwcxxn4pWlZ1qUs2+7XoOehXC63PY5+Ya/FMIxzBRzHcVSr1bzAZpqmCoVC1wXVWmfuAJOMYAJcgGw2q5WVFV+TTMPhsCzL6ug1aX3ZFYvFS12Mq5d+x3ddV47jqFKptIULPz0chmEoFospl8v1DVStY4x7HZVez0OxWFSlUlGpVPJO9/UTTFoTUgcFuJbWOion5XK5tufFcZyuE3T9hD8g6AgmwAUoFotDj/lns1mVSiVJzXknrd8jkcjI1w8ZdPxecy8sy9L29nZHfafnnbTOIspkMj2/OPsdY1Rhpdfz4LquUqmUt0ZIOBzWV1995avOzc1NPX36tC2Y2LatQqEg13WVzWZVq9X06NEjZbNZPX36VE+fPpX0j56xeDyufD7f1p7TPVqtYZxZWnsG04kF1gAMlMlkZNt22/CT67regl6nvwy79brcvXtXpmmqVCp1/fLsdgypeQpxv2GtUSiXy8pms2eaz+I4TlvguSytRdjG0dsGXCTOygEw0Om5F1KzJ6PXf+jdeo9a80x6/Uff7RiZTGbsQztSczjq9NwNv9elaQ3lXPZZNNlsVo8ePbrUYwCjQDAB0JPrusrn87IsS9Vq1TsTqNWD0Pq9pXXGSKFQ6PjifvToUdc5GaePUSwWlc/nvYXh/M7NuEzRaNSb+Fsul1UsFoeay7G9ve09X5fBtm1Fo1GGcTAVGMoBgBHI5/NewLlIrusqkUhc+lARMCr0mADACKTTaVmWdeFDOsVisesEZGBS0WMCAAACgx4TAAAQGAQTAAAQGAQTAAAQGAQTAAAQGAQTAAAQGAQTAAAQGAQTAAAQGAQTAAAQGAQTAAAQGP8fcN++miVlSAUAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "beam_spill_times = norm.rvs(loc=0.0, scale=1.0, size=timing_sbn.shape[0])  # nanoseconds\n",
    "\n",
    "timing_with_spill = 1e-9*beam_spill_times + timing_sbn\n",
    "nu_timing = 1e-9*beam_spill_times + tof_light_sbn\n",
    "\n",
    "timing_bins = np.linspace(-4, 10, 50)\n",
    "\n",
    "\n",
    "plt.hist(1e9*(timing_sbn - tof_light_sbn), weights=weights_sbn,\n",
    "         bins=timing_bins, histtype='step', label=\"HNLS\")\n",
    "#plt.hist(1e9*(nu_timing - tof_light_sbn), bins=timing_bins, histtype='step', density=True, label=\"Neutrinos\")\n",
    "#plt.xscale('log')\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r\"$\\Delta t = t_{HNL} - t_c$ (ns)\", fontsize=14)\n",
    "plt.legend()\n",
    "plt.xticks(fontsize=14)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.hist(1e9*(timing_with_spill - tof_light_sbn), weights=weights_sbn,\n",
    "         bins=timing_bins, histtype='step', label=\"HNLS\", density=True)\n",
    "plt.hist(1e9*(nu_timing-tof_light_sbn), bins=timing_bins, histtype='step', density=True, label=\"Neutrinos\")\n",
    "#plt.xscale('log')\n",
    "#plt.yscale('log')\n",
    "plt.xlabel(r\"$\\Delta t = t_{HNL} - t_c$ (ns)\", fontsize=14)\n",
    "plt.legend()\n",
    "plt.xticks(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}