PS191_other_channels.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "reliable-grounds",
   "metadata": {},
   "source": [
    "# Batch computation\n",
    "\n",
    "This notebook computes event rate at T2K and PS191 throughout the parameter space. Sometimes we compute the lifetime as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "likely-jefferson",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy import interpolate\n",
    "\n",
    "import pathos.multiprocessing as mp\n",
    "from pathos.multiprocessing import ProcessingPool as Pool\n",
    "\n",
    "from functools import partial\n",
    "\n",
    "from particle import *\n",
    "from particle import literals as lp\n",
    "\n",
    "from hnl_apps import *\n",
    "from hnl_apps.plot_tools import *\n",
    "from hnl_apps.const import *\n",
    "\n",
    "from Nlimits import *\n",
    "from Nlimits.constraint_dict import *\n",
    "from Nlimits.plotter import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "medical-piano",
   "metadata": {},
   "source": [
    "# Minimal version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "eastern-depression",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set up grid\n",
    "gridx = 100\n",
    "gridy = 50\n",
    "\n",
    "M4 = np.linspace(0.02,(m_neutral_kaon - m_mu), gridx)\n",
    "USQR = np.logspace(-9,-2, gridy) \n",
    "X,Y = np.meshgrid(M4,USQR)\n",
    "m4_grid, usqr_grid = X.flatten(), Y.flatten()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "combined-grass",
   "metadata": {},
   "source": [
    "### mu-pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "compound-orlando",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<multiprocess.pool.Pool state=RUN pool_size=8>"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# compute rate and lifetimes\n",
    "tot_cores = mp.cpu_count()\n",
    "pool = Pool(tot_cores)\n",
    "############\n",
    "\n",
    "get_event_partial = partial(hnl_tools.get_event_rate_mode, modes=['nu_e_mu', 'mu_pi'], flavor_struct=[0.0,1.0,0.0], detector = ps191)\n",
    "event_rate_weak_ps191 = pool.map(get_event_partial, ((m4, usqr) for m4, usqr in zip(m4_grid, usqr_grid)))\n",
    "\n",
    "############\n",
    "pool.terminate()\n",
    "pool.restart()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "honest-module",
   "metadata": {},
   "outputs": [],
   "source": [
    "event_rate_weak_ps191_e_mu = np.array(event_rate_weak_ps191)[:,0]\n",
    "event_rate_weak_ps191_mu_pi = np.array(event_rate_weak_ps191)[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "planned-google",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save to file\n",
    "table_PS191_minimal = np.array([m4_grid, usqr_grid, event_rate_weak_ps191_e_mu, event_rate_weak_ps191_mu_pi])\n",
    "np.save(\"limits_computed/PS191_minimal_Umu4__mode__others\",table_PS191_minimal)\n",
    "np.savetxt(\"limits_computed/PS191_minimal_Umu4__mode__others.dat\",table_PS191_minimal, header=\" m4(GeV) USQR event_rate\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "funny-tennis",
   "metadata": {},
   "source": [
    "### load it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "defined-continent",
   "metadata": {},
   "outputs": [],
   "source": [
    "table_ps191 = np.load(\"limits_computed/PS191_minimal_Umu4__mode__others.npy\", allow_pickle=True)\n",
    "m4_grid_ps191, usqr_grid_ps191, event_rate_weak_ps191_e_mu, event_rate_weak_ps191_mu_pi = table_ps191"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "id": "quantitative-diversity",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-69-2f0333389b14>:42: UserWarning: FixedFormatter should only be used together with FixedLocator\n",
      "  ax.set_xticklabels(np.array([20,40,60,80,100,150, 200,400]))\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUEAAADiCAYAAADQ8XXFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAArxUlEQVR4nO2deZwcZZ3/308Swg1DwiEIlTAgIIeEIdxylRNRwANJQFtWdxEmUsj+XBViXM91XTaux+pqwcRboBeTLK7CIpChuM8kI0dEzgGLQ+UIQxAIhOT5/fF9OlPT6Z6Z7unuqu7+vl+vfnV3nd95Zuozz/E9jLUWRVGUdmVC2gYoiqKkiYqgoihtjYqgoihtjYqgoihtjYqgoihtjYqgUjOMMbMTn3uNMQuMMR2J77PLnjzCtdK0Y7y2KNlnUtoGKC3FlKLPC6y1g+77YmttX5XXStOOcdlijOkBBoAOa+2Saq+j1A/tCSr1otNaO5C2EaRohzHmAmC5E91ZadigjI72BJWa44aexcKzqsyxnUjv7GD3eba19ptljusqcYn+ciJXDzsqtOUQYEPvzxjTZa3tL3VNJT1UBJV6MBNYVvhijOkeZQh6oTGmGxGshaUOcOJSaY+u5naMwxaAwSrOUeqMiqBSDzoY/sB3woaeWTcwxVpbEJkOa+0SY0wv0hMr16uruCc4gh1dzo6pwK9c72xMdlRoyzJnA7BBPJWMoSKo1BwnJhe4XhUJwesB+sosEDyGE6ky16y491XODmttvzFmSple4Yh2VGjLQuA0YwzA0kpsVxqHiqBSF8rMp/XB0NyY61GdDvQjgtHTIDsgIXT1ssOtSBf+AehcYEZREVQaRvGigOtRzXOfB4FyglUPWxYmPqdmh5I+RlNpKYrSzqifoKIobY2KoKIobY2KoKIobU3LLYwYY54HnhjtuF32zU3Z+9gLd7/v6jMffuFP179cf8uUlDHuPXOT4FOnvWvrfY7/1u4rr5376Et/vvvVtO0pQWbbrkKmWWt3KN7YcgsjxpgV1tqDRzvOD+ItgD8DdwDvjUKvtRpCGYYxZjKAtfaNtG0pxg/iE4BrgCOj0LsjbXuKyXLbVYIxZrm1dmbx9rYdDkeh9yrwz8AJ1ME/TVGU5qBtRdARAtcB3/GD+G1pG6MoSuNpaxGMQm89cCbwOrDYD+JtUzZJaU8eAM4BHk/bkHakbecEk7g5mauAh4EPRKH3aF2MU1KjVea10qBV2k7nBEcgCr1rkbnB3YCvp2yO0mb4QbyNH8RdfhBvlbYt7UjLuchUSxR6kR/EDwHbpG2Lkj38IP5PYEaFp90Thd6nx3DcEbjVYcRbQWkg2hMczpvA5LSNaFeMMZ3GmB5jTLd7zS4USTLGLHXHzHbFki4oOne2MWaxO29B4b2wD/gGsG3hOiXu2VMoxqS0F9oTHM6rwOZpG9HGLLDWzgEwxiwGzgYOsdYOGmMK+fv63LbiTC99wCxrbZ8xZq5773QJVPuAw4CX2DgP4ALgQqSY0mmUySg9xh6d0oSoCA7nZWD7tI1QAMn0POgSksLwim8vlDln0L0PqyMywnUK+zXXXxujIjicvwE6OZ0e89zQdRDpoQFSGwTodMlPO4E9Spw7M3HMTPde8BIoiJxfOCaR6n6eK4u5HBjUFPjth4rgcF4Gtk7biHalVNp6a+0897EgaAO4DNVFx/UltheOnZs45J/dcQcXnTdiUaUGcT9wBqCuWSmgIjgc7QkqDScKvWeAy9K2o13R1eHhvAxs7gex/nNQGoYfxFP8ID7GD2J1z0oBFcHhvOjep6ZqhdJuHALcBOyXtiHtiIrgcJ5w79PSNEJRlMahIjicJ9z79BRtUFLAOWR3p22H0nh07ms4T7j36SnaoKSAK7VZqhi70uI0hQg637EpwPJ6OrZGofeyH8RPU3mMqFIlzp+vmyHXmA5kjuxCYLG1dpb7/c8CHktGirjtpwO9bv9SJGqk4G+44TrASeXOc8d3Fz6Xs9UP4htLbF4UhV7oMpVfndwRhd5xY20HJT0yPxx2D0k/8l967iiH14JbgGP9IDajHqnUggXW2oXOz28urjfmembJULnBMqFyq9y5ne79sUSoXPF1kuctSxzf7T6nxe+BDwAPpWhD21L3nqALSl8A9BZ6cc5DfwDx/l/k/lBLYq0dcNfoIRFFUEduBj6MRCWo82pjaVioXPKeYzVupJ6dK9dQdv8o130W+G015yrjpxE9wZkk/ghdz67D/eddxBiELTFfM7tONiYp9Ag+3IB7KS5Uzg1Hy4XKzaSyULmZxdcBdi869xDXY8QlWygc13D8IN7BD+KT/CDeLo37tzsNySxtjOnF9QRdCqRBa+1Ct+8xa+0ebp5mGNbaJW57n/vPPmrWaGNMP3D4eOw97pMDv4UJM55e+cu3PXLrl18fz7WUzLCJe18LbAucBXw7PXOGOPrMlbMmbbr1VWvXrDrm1p8ddFfa9pQg2XbNzO1ZySw9leFDlykgglf8cvsHkP/2PUDJSWuXC265MWY5NXB0fvONV75rjNlp57ef9pHxXkvJJO9iKL5YaXMyvzqcWA0uuyrsepWFnuWK8dZC8IP4OuD3Eydt/uXjz/nTVS62U2liCvOD7m/jcvfKBH4QvwmwyWZT3sxiHY+itmtaEnPEw0ijJ/gCI09Up44rxH4mMmn+Wz+It0zXIkVR6kUaItiHW5Fzq76ZdFCNQu8eZHHkIOASP4gz706kKErl1P3BdgsbM4HTXTLLfmDArcb1ICnUM0kUelcBnwFOAb6vvoPjwxhzQaJuSG/ic3dxzZAS57ZyWNtyZJ7ygbQNaUfqPifoFjiWlNgGGe0FFvF9YB/gXOASIIurd81Cv3NHAZjiVv8HERHoGunEVg5ri0LvBSBK2452JfMLI2kThZ71g/hC4JPIg6oiWCWlojKcKHYgfnvdQJe19puJcDoQf9KZSGjcMqoLoTsM+KYxJm+tnZU4b9h9yjru5011vmQ5O+rowQ/inYFjgb4o9J6v6j5K1eg819h4Esk1OCNlO1qZQvhbwSl6AdJDHABOSwho1SF0lK82t+E+tfyBKuAA4L+BvVK6f1ujPcEx4HqD9yAB+Up9GCzeMEKyjMZXmxtDj05pTrQnOHaWAgf5Qbxb2oY0M4WhL3CwG4rC8PC3wvs85wTf5eoHF8Laqg6hI1FtLrFt2H1q+bMqzUFDwuYayVhC66rBD+JCQoUvRqH3jVpfX6kvxpjJkE2HXz+I3w1cCxwVhd7tadtTTJbbrhKMMcuzEjbXlESh9xhwA3Cm+gwqSuugD3Nl9CJDsp60DVFairuAQ5H6w0qD0eFwBThn6euQLDUHRKH3RD3uo9SeVhnSpUGrtJ0Oh2uAiyk+C7DAL/0g3jRlk5QWwA/iXf0gPtsP4p3StqUdURGskCj0/oQ4Th8NXKZCODLNEiqXcljevkgWpFKr3kqdURGsgij08khM8anALX4Q75KySVmm34VJLgNWuM9LEQflEbHWDjaq9kcj76VkCxXBKolC77tIYoUDgC+lbE5mKRcq5z4ekuwROn/AHvfqcPsWuN5jb3HP0W1fnDiu2xizILFvgbvO0pHOc9s2fC5L3tgRXj2J43qqDrNTGo6K4DiIQu9/gSuBU/wgnpyyOc1I3UPlmqTanJIiGjY3fn4CzEHmCb+fsi3NxmDxhkaFyiXoGGHfcMYaOpcbynSuZB/tCY6f65CexVf8IO5I2ZZMkmaoXDNUmwNuA94O3JPS/dsa9ROsAX4Qz0BqoPRGoXdOI++tjI2kr5sT5Z4Sw2ulBOonqIyKS8X/beCTfhCflbI5yuh0k6GMQH4QT/eD+DPqZZAOKoK1Yz7i+vEjP4jnayr+7OJKus5J244EeyH/RKenbEdboiJYI6LQexN4P5Ic898QMdxk5LMURUkbFcEaEoXeGuCjwL8Cn0BXCBUl86gI1pgo9GwUel9ChPDv/SA+JW2bFEUpj4pg/fgX3IqxH8Q7pm2MoiilURGsE1HorQU+BmwD/EwTLSgjcDOwK7AibUPaEfUTrDN+EAfAD4E7gNlR6D2TskltSav4uqVBq7Sd+gmmRBR6IVLK8R3ACj+Ij0rZJCVj+EG8hx/EX9YiXumgItgAotBbjGSjfgW4wQ/itOrbKtlkD+BrgIpgCqgINogo9FYiMa13IVmpj0jZJEVRUBFsKFHoDQIfBJ4CfuMH8TtSNUhRFBXBRhOF3gvAicAbwE1+EB+dskmK0taoCDaSvPkH8mZF9M5pewJHAn9F5gj/VV1oFCUdVAQbiwd0AadEoRcDRwCXAP8M9PtBrIV22pMI2Ba4O21D2hH1E2wkefMO4F7gOWBncnYdgB/E7wEuBdYAx0Wh92h6RrYmreLrlgat0nbqJ5gN7kfqXeyADIcBiELvGsAHNkWq1x2XinVKKvhBvLcfxP/hB/H0tG1pR1QEG0nOWuDX7tuHkrui0LsPOBZ4CYj8IP6GpuJqG6YBnwM0qWoKqAg2niXu/ZPkzUHJHVHoPYDUyPgp8AXgXj+Ij22wfYrSVqgINp67kAp1mwFXkDdTkzuj0HslCr2zgPchw+Pr/SA+VzNVK0p9UBFsNDIk/hRSX3c6cBl5M7H4sCj0rgJmANcAPwD+1w/iHRpnqKK0B00hgoUyja5KWPOTs2uAU4HngROAr5Y6LAq9l5GU/Z8B3gPc7wfxbO0VKkrtyLwIOuHrNcZckLYtNSVnY+B0YD3wRfLm5FKHRaG3Pgq97yJxx88Ai4Gr/CAurqOrNC9LgYlIujWlwdTdT9CJ2AKg11rb77b1IK4incAia+3gKOd3AgMjHZc4Prt+gqXIm88DFwKDwMHk7EC5Q/0gnoQMpb+OPDT/AnwnCr2m9t9qBK3i65YGrdJ2afoJzgSmJAzpBDqstX3AIkQgy5IQvilOPFuNBcBvgA7gSvLmreUOjELvzSj0/hN4O/A7RDz7Nf64ufGDeF8/iC/WiKF0mFRqozHmXcjDuQyYB0wFeqy18yu9gbW2zxiTrPE6G+n1YK0dNMZ0u3vOLnHuErd/uTt2DqNXcDOF/1zNwrXzOKv7APaaYNjXWu76S2jev8u5rBzhlGeBjxz9iT+cOHGTLb9njLn5+E8+/us1f3t6/p2XHf14o+xuMjYBMCZ706lHn3n/9EmbbjN37ZoX88ZMezJte0qQ2barBeV6gl2u2/h5pKj4KmRIWgumuusVmAIbCmIPe7n9y4FOJ5Ile43GmB5jzHJjzHJ3/abihAWsvvMR/PXruc0Y3vqWbblx8Ed0j3beLT/Z7+pnH73ywPXr3vgXzIRZm22924pjzn5o7i775jZabVYUpTQl5wRdT3CZtXa1+34q8Hlr7SFV3cSYXtycoDFmgbv2ErfvRWvtdlX/BBvfq7nmBJPkzWbAz5EFk3XAOeTsj8ZyqkvN/lOgG3gcGSr/QucLhSzPa/lB/G7gWuCoKPRuT9ueYrLcdpVQ6ZzgcuRBBMBa+z9Ir7AWvEBijlBJIK4zOeDfkYWPheTND8ibUcPnotB7Eng34lLzPDJt8KgfxGf6QZx5LwBFSYuSD4e19iVrh/dArLXX1+iefcgiQGHlt69G120NcnY9OTsfOBNJvHoucB15s/1op7rC71cChyF+hU8j0Sk3+0F8QB2tVsbHOuBVxF1KaTAjusgYY84GxuJDs8pae0WZa8xG5hX7kCHxgNs2iOTWWzgW15ex0tTD4WLy5gjgCuAtwBPAB8jZ+8Z6unOq/jjwLeQfzw+Bb7tchm1Fqwzp0qBV2q7ccFjzCWYdcZn5NXAI0lv4GDn7P5Vcwg/iqcgc4Zlu02LEv3BZLU3NMq3yIKdBq7SdimAzIwsmC4G/c1u+DnyVnK1o+OQHsQecB/QA2wA3A1+JQu/G2hmbTbL8ILuCW/OAr0Wh93Da9hST5barhJo4SxtjtqmdScqYkQWTjwOfReaNvoRkoNm6kstEoRdHoXc+Ut/2n5B6tzf4QbzUD+KDRj5bqSNvQRbERp33VWpPpauGPzbGzKiHIcoo5KwlZ78DvBeZT/0AcG2lQggQhd5qF3myJ5KcYQaw3A/in/tB3FUzmxWlCahUBOcBjxtjfPc6yxjzoVHPUmpHzl4HHAoUCjVdTd5sVc2lotBb45IzvA34PhLNs8IP4tv9IM75QdxUkTeKUg0ViaC19nFr7UuAAT6JZEEuG/Cv1ImcfQQ4Hini/k4k5niLai8Xhd5gFHr/BOwKfBoZll0GDPhBfIEfxB3jtllRMkqlc4KfM8ZcCxwEnG2tPcdae09dLFNGRrLN+MCfgeOA35I3m4/nkk4MvwfsgxSIfxAJVXzSD+LvaiGguvEGUoN6bdqGtCMVrQ4bYxYBve7rMmvtamPMNoXwuizQkqvDI5E3ewM3ATshoVcfdAspNcEtmHwG+DDyT/M3yND5pij0msa1oFVWONOgVdqu5i4yxpiDkPC32dbac8ZpX81oOxEEyJt9gRuRUp7/B5xKzr5ey1u42ORzEPeaqUj50P8C8lHovVLLe9WDVnmQ06BV2q5ufoLGmG3dPGEmaEsRBMibA4AbEIH6DTCHnK358MoP4s2BjwD/CBwIvAz8CvgZcEdWe4dZfpD9IJ6BJMj9vKs4mCmy3HaVUNOkqm5leLrzG2w/wckiOXs/kkHmRcR9po+8qXmSzij0XotC76fIvPDRSAnRjwC3AX/0g3ieH8RaP7cydkCqC2pikRSoqidojNkdeeAOBqwOhzNE3hwMXA3siITZfQH4r0qjSyrBD+KtgDlIWN47EYfuq4GLgGui0Es9MUCWezN+EHcjdUaOiULvlrTtKSbLbVcJ9RwOz8jSCnHbiyDgMs58D4lCAOmlfYKcfajet/aDeC/g74F/QCIhBhAx/GkUeqtGOLWuZPlB9oP4XUiCkWOj0Ls5bXuKyXLbVUJVIlhu5dcYcz5SOwQkO8mcrKwQqwgmyJv3AxcDOwNrgC8D3yFn19X71s7R+hSkMNQ73f0vR/wPb4xC781625Akyw+yH8THAxFwXBR6N6VtTzFZbrtKqFYELxytrojLQt1dTf2ReqAiWETebAd8F4k9Brgb+AdytmET8H4QHwgESM90K+A5ZC7xV8CtUejVXZSz/CD7QXwYcClwRhR6d6VtTzFZbrtKqFYEFyH/vQ2SV7DfWvtEvYysBSqCZcib9yKZaHZFnHO/DHyrEb3CAm5l+UQka/nJwOaIs/flwCXAPfVaXW6VBzkNWqXtqhXBYe4vzjewE3jRWhvVxdJxoiI4AnmzLfAfwNluy41IfsKGVzhziyknI07YJyIVzf6AiGHelQuoGa3yIKdBq7RdTRdG3Orw7lkUQhXBMSC9wp8hUSaDQA85uzgtc1zS19OQfIlHuM23Ij3EJVHo/XW898jyg+zyCX4P+GwUev1p21NMltuuEqrtCX5ohLT5ZfeliYrgGMmbHZH6Iye7Lb8AziNnX07PKHAFyD/sXvsj7jYRIohXRKH3YjXXzfKD7Afx0UiC21lR6GWu5k6W264SqnWW/qYx5iJjzIeMMdOL9qljZzOTs88ilekCZOX248A95M3haZoVhd5jUeh9Iwq9A4ADgH8Ddgd+DPzVD+LfujRfFedRzDCFnkhrVjfPOKP1BM9H/LxmIS4xuzNULW6ptfZbDbCxIrQnWAUSe3wZklx1HRLC9W/kbEPdWMrhCkYdjPQOT0cWd9Yic5pXAVdFoTdiSrcs92b8ID4S8eU8IQq969K2p5gst10l1GxO0C2ODFprH6+VcbVERbBK8mZT4F+Bz7kttwNnkMvW79nVUD4SCQ08GUn7BfAAkjziSuD2YrebLD/IfhAfgbT3e6PQuyZte4rJcttVghZaUsZG3nQj84O7IMkRvg6E5GwmM8X4QbwncBIiiMciq8zPIUkkrgCiKPRez/KD7GpC/wI4Lwq929K2p5gst10lVLswUqg7nJyrKHVC2brDjUZFsAbkzVTEp7BQOuFZJLnqxeTsq6nZNQp+EG+DFJ0/BRHGrYHVwFXr1r5y5aonb772/mvmVrWw0s60tQg2IyqCNSJvDHAC8DWkpgnAX5D6xQtrmbi1HvhBvCnwLkTIPwBsb61dY4zpQ4bMV0Wh90yaNjYLKoJNhopgjRExPBERw0K7PoOs2v641slb64EfxJPWrhk8duLkrT44YcKkk5AFPoDlwG8RUbw3rVyIfhDvjSxMnR+F3g1p2DASKoJNhopgnRAxfB8ihjPc1qeAbwA/JZftB6TwIB9/zp/WAvshP8v7gcOQ6Z4ngd8B1yHziA0bNvtBvB+wEjg9Cr1FjbrvWFERbDJUBOtM3kxAhpdfQ/z4QHqGFyPD5HFHd9SDcg+yH8Q7IfOH70OGz1sjDtp3I4J4HXB3FHp1K4LkFnceAT4Whd4l9bpPtagINhkqgg1CxPBDwFeQyA6QxAy/QpK4LkvLtFKM5UH2g3gTZP7z3e51KBJQsBrxSVyKiOIjtRw6u/otMXB2FHo/rtV1a4WKYJOhIthgZJh8PFJz5P0MeRLciRRiWpKFoXI1D7IfxNshZU1PQAIGprtdMSKIS4Hro9B7fjy2+UG8A7ICf14Uej8Yz7XqgYpgk6EimCJ5szsShncWElUEsqJ8MdBLzv4lJcvG/SC7qJU9EDGchYjjtojL2O8ZEsXbotCraOXcD+ItkfnI70eht6Qa++qJimCToSKYAfJmS+CjwHkMDZXXAouQ3uHd5Br7h1frB9kP4klIKGlBFI8AJiFx2Lcg4aV9SI7E1GusjAcVwSZDRTBDyFD5WGSo/AGGEnb8EXEJyTcqLK/eD7JL6HAssrjSzZD4v4BkwekD+kaLcc4iKoJNhopgRsmbaUjx9jOREpMFbkUEcTE5+0K9bt/oB9kP4p0ZEsRZSBgiwOPATUjqrJuBgSj0rB/E1wEro9D7TCPsqwQVwSZDRTDj5M0miDCcgYS3be72rEXKdF4KXFXriJQ0H2Q3n7g3IobvQuo1F1LRPYOI4ZHAE8DxWRs+qwg2GSqCTUTebAV8EBHEWQwNl1cjhZguBW6qRc3kLD3ILhPO24FjEq9CT3EV0ju+DbgDWB6F3mtp2FkgS203HlQElWyTN29B8gWewVB4HkhUyn8Dl5Kz91V7+Sw/yK6nuAQ4Dsl+czSwp9v9JrL6fDsiirfXuv7KaGS57SpBRVBpHvJmH2R1+QyGfPMA7kd6h/9daXGorD/IfhB/FfgSsFUUeq8538HDkWHyEYjjdmHq4GnED/MuJLJlRRR6f6uXbVlvu7HStCJojOkGOq21C8d4vIpgqyCry0ciYngaQ/NoFllcuARxxl492qWy/iD7QdwNfAL4dKnCUi6a5R1IexyJiGKn270eSSp7N9JjvA14qFZRLVlvu7HSzCI4G0mpDtBhrZ01yvEqgq1I3kxGcgWegUSmbOr2rAH+FxHE68qVBGiVBzmJH8TbA4cggngo0nMs/KNYxZAg3oXMLVZVRKtV2i41ETTGdCAJOXuttf1uWw9Su6QTWGStHRzpfGvtoDGmE8BaO1otCRXBVkfqJ89GSnQem9jzLDJ/+Evg90mH7GZ4kN3c4NRqw/ASq9BHIb3Fo9x3GOot3sXQMPoPUeiNWkemGdpuLKQpgt3AXOBCa22/E7PZ1tpvFgTSWjt3DNfpGcuQWEWwzcib6cj84d8x9MCDPPCXAJeRs082w4PsB/EPEbehXWvlJuNqOh+KpAw7zH0u9BZfBVYgeRULr0eL790MbTcWUh0OG2N6cT1BY8wFSKGmhW7fY9baPdywdxjW2iXumC5kXnDUuEpjTD8yLFDaiD13gtu+ysHbbclHJ03gdGPYHsBa7HrLjS+9yuU/v5nffPYyMpte/5iz/njqxE22yL/5xt9OueUn+11dj3tsvu10Djz50j0mbz71EDNx8qHGTJgJ5kBjzGYA1trVYPvt+nUr7Pq1K9auGbz/wRvPj1986tZ1iC9nM3N7VkRwAbAsIXAvWmu3G+X82YhwlixM7YbXPe7rDsDbame90mycfxKTvngK795yU86YYDjZGJk/tJbX1ln+77U3WPSLm7nmvF+QqazYu82YO2mPw+c/DPbBGy/e/cRG3ne3A896+6TJWx9sJmxysDETDgZzQKIH+Dp2/cPASmvX/WH9ujceWPvaC3948MYL/jT4zJ3ZXlQYTvOKYIX30uGwMkTedACz16/nYxMmcHRiz2rg18DlwPXkbCZ6OX4QzwP+HfhcFHrfTtGOyUj88wHr1619h5kwYT9jJu4L7JY47BXgD+610r0/DMTFJU/TxA/iicAnbrhoWk9WRLB4OKwiqNQdY8zkB7/FbnvvzCnAR4CuxO7nEWfly4FbahGhUi1OfC5FnKWPiEIv9d5qck7QD+JtkfIE+yfe9wd2TJyyFln4fCTxegwpYfBkPX0ai3EJa38JHHfDRdNWZEUEu4DuxMLIj6y1c2p4LxVBZSM2mtzPm70Q16uPICFsBZ4GFiORG7el0UN0abq2ikJv0C1svJpm6NwYs3LvgIjinsh01Nvc5z0ZcvIusAoRxLjEewz8dTzi74T6SOT3W1hrOO+Gi6adm9bq8GxgPpJKqNdaO1CY40P+Gy8cyUWmivupCCobUfZBFofsAxAx/DDDI1RWI4lSrwZ+R87+uRG2JvGD+ApnX09alejGszrs4qR3RhLS7uZennsVPneUOPVvSA/9eSQd2fPI7+N1pIzD60gW8y2BLdz7Tsjvb0+3rxCDfmEUeo82rbN0pagIKqUY04Msgngo0ns4Edi36Ih+RBCvRhLD1n3eyw9iH1iIiMilwIIo9FbW+75JGpSLMSmOOwLbA1OL3rcGJiOO8psivo+vInOTryJC+SckvPJ2JMv3hh60iqDS1lT1IIsP4nuRanQ+w4d1LwDXAtcDETn7RI1M3Qg/iLdA4or/Een1fDwKvV/W637FZNFP0DmGU0looIqg0taM+0HOm82R6JST3Gv3oiOeQDJIR8AN5Owz1dpaDjc/eA7QG4Xec34QH4M4Pl9Zz9XYLIpgNagIKm1NTR9kGTbvhVShOx5JgdVRdNRDDInijeTsuCrSlcIP4kXAHKSY1eVIhu4VtSwHCiqCTYeKoFKKuj7IeTMROBAZMvtIPsCtio66j0IvEW4mN/7FQLeK/H4kqcRJyHzZFVHonTreaydREWwyVASVUjT0QZYSAjMRQTweSWSwWeKIQpnOG93rlvGKoquRfCrwUhR6i91iwx0MlQK9RbPIqAgqbUyqD3LebIbEsx+PCONhwCaJI+ohirsDP3T32xRYBywDzo9C79ZKrqUi2GSoCCqlyNSDnDdbINmij3OvUqJ4DyKItwK3V1u43g/izRHH4YIAnxOF3r1+EJ8CfNrd4y5gWRR6z5W6RqbabhyoCCptTaYfZBHFwxkSxcMZLoogYWi3MZQo9YHx+Cn6QfwhJIihi6ECV48DXS5SxQNei0LvuUy3XQWoCCptTVM9yEOieCwyn3g4EhGRZDUy51cQxrvI2Ypjct3cYRfiJL4PcJarg5xHomj+bO26++z69fdNmLjJ3VHojZrOLms4n8KJN1w07U4VQaVtaSoRLCZvJiGhc4Vs0UchkRVJ1gP3IqJ4p3sNJLNrV4IfxIe5+xxo7foZYN5ujHkgCr0Zbv/FiDAnkyQ8EoXeS9Xcr144R/MfAatvuGjaIaVEcFLjzVIUpSKkbsrv3euHAOTNrgwVXToKOCjx+pQ783ny5m5kzu9OJNRvcCy3jEKvkIYfY8zk3WbMnbTnEV9IZnuagvRQz0hsuwp4H4AfxCES//sUkhzhKWAgCr0XKvjJq8YP4rcgxbnOQTKOfwmpx7IR2hNU2oKm7gmOhbzZEnnIj2Iolf6OJY58iCFRvAu4f7RMOSO1nVt42QPJGvNSFHqRS5rwINJb3TRx+A+j0PuUq5wXIcL4FPAcEoZ4VxR6K13+v6m4hAnlnL/dfbZD4op3AP4Yhd4LfhB/FImzBllg+nwUeteWmxPUnqCitAI5+wpDLjaFqJZpSG+tIIpdSK9ob+Bj7sw15M3vkaQD9yPJUe8nZ8fUY3MJCla6V2HbemAvNxe3PbArkiDhKXfI1khR+ZnABxnyoZzvruMhC0EA+EG8BngN+Kco9H7hB/GBwC1IHPXEhDlnIFEzdyA9vyui0HtgtJ9Be4JKW9DyPcGxIGVL34EIYkEcy5Wi+AtOEF97g0fWrOXR7bbkQeDpWiaddUK5JdLz+5vryU0BcohYbo6I5ObAkij0bvKD+K3AZ5HMMc+51/PAfVHolXUl0tVhpa1RESxD3kwBZiALL/sn3otXowu8hmSJfhRZDEm+11Qga42KoNLWqAhWQN5MQIbSBwD7v7mOfSYY9pgwgT0pPc9Y4HWGBPIJhub8nnbvz5Cza+pp+kioCCptjYpg9Qxru7zZhqG0+W8r+rzTGC63CngJWTkuvF4p8X0tkh0axJl7W2RFeqp7D8nZfIU/hy6MKIoyTnJ2NZJhu3+jfXmzNUMrxR7wVmRRpPC+CyJgUzY6t3Kuq8E1ABVBRVFqRc6+jLik3FNyvwyzt0fSjJV7bYksiBRWfQtD1ZcQN5pV7vVorcxWEVQUpTHIosmz7pUZJox+iKIoSuuiIqgoSlujIqgoSlujIqgoSlujCyOK0mrkzYeRFdSb03RObhZUBBWllZDECd8BdgZeJW9uAH4H/I6cHRjx3DZFh8OK0lpsBvwc8dXbAinF+QPgMfLmYfLmhPRMyybaE1SUViJnXwO+AHyBvNkZeA/wXuDdSCTHKxuOzZtzkRT+tyEFne51CVzbCo0dVtqCto8dlhT9hwIryNnX3bbfISJZ4BWkBymF4nN2CbRO22nssKK0M9LDu71o6/8DjkGyUb8TSYRQqGGyDbAEYOUCdtlnF75P3tzLUPLVR1ql16giqCjtSs4+DDwM/BiAvNkROBBJvPrHwmG7TeXAiRN4H65+iON18uaPiCDOI2f/7K5hqi3uNCqSFHYGsJqcfbBWl1URVBRFyNlngaXutYFH/sLyvXfmo1ttxn6IQB6A5Buc4V6fShx+NXmzL5JP8HH3/gTwV+BRcvaRYfeU1exN3GsSsrDTgVTKW+uOOR84GRnObwaEwLnj/4EFFUFFUUZk5hd5DlhibSJ/n+QV3B/Yw6XXKrAHkkbLQ4baSS4CAnf+IUixp3IeKnswVGfksMS1HkIEtWaoCCqKUjkifLez8Tzj/khRpd2B6Yn37ZEFlwLbMSSAb7rXWiQ79SAwOXHsfyJuP3eSs8/X6kcooKvDSlvQKiucaVCXtpNh8CTgzbrNIRahq8OKomQHEb4R6x03Co0YURSlrcl8T9AY0wGchgSED1hrN65toCiKUiXN0BPsAZZba5cgFeoVRVFqRt1F0BjTYYzpNcZ0Jbb1GGO63XvHKJcYTHwe7VhFUZSKaERPcCaJEnvGmE6gw1rbBywCFoxy/iJgpjGmm9qU6lMURdlA3ecErbV9xpg5iU2zcb07a+2gEzeMMbNLnLvEWjsILHTHLB/DLU1hSV9REmwCYIwZ7ThlY1q67dJYGJnKkCc4uN6dm/PbCNdz7HLHzStzTA8ydwiwLxs7cCoKwLZI/Vqlclqh7aaV2pj51WFr7QDDRbPUMQtJ9BZLOUQqijFmobW2Z/QjlWJaue3SWB1+AZ3bU9LhyrQNaGJatu3SEME+3CqvWxnuS8EGpQ2x1rbsg1xvWrnt6j4cdgseM4FBY8ygtbbfGNPpFkS6gLNrfMuFNb6eoigtTMslUFCUWuFGKjORf9b9zq1LGQPN1HYqgopSBud1sMi5cml2ogpoprZrhrA5RQE2RBr1uAikjnFea9RIJmvtQvcQdwBj8VHNLO7n6jbGLHBuZ+O5Vku1XcuIoDFmQeJzJWF5ShPg5pCXO3eopcCPxnnJSiKZ5lPGR7UZcD/bPPezLWP8P0tLtV1LiKD7JcxOfK4kLE9pDjqA093nfmSuqWrc38eqxKZhkUxAMpLpQprYrctaO2CtneW+dgIrxnm9lmq7lhBB5BdbcKgu+QtRmhsXQlnoUXQhQjgMN0y7IPG9y81NjYWpDH+wp7je53yk1zm3OsuzgxOlqa43Xbyvbdsu8xEjo2GM6XbxyYUHpGRYntJSzAXmFG90c1AL3cPcB8ws9cCPFdfjyeyEfqVYa5cU3NOKV2vbue2aWgTdfN+q0Y5TWgf3kM5xvfyNcA/zALAgMQQcCy0byVSYF3dt1of0zjYSqHZtu2YfDvcAna6bX3hv6l+IUh43xCqsOpac5nArllOAOcnh3Rho5UimHoYnJO4odVC7tl1T9wSttd8sfDbGzHfd/S6GJmY7aLJfiFIa93tdDKxyKZ36KPrdukWxDcO4wvAu+XeSOLbRkUxpshDodj/bHEpMJbRz27WEs7Rr/MWIG8BC90saRH4hC8sNnRRFUVpCBBVFUaql2ecEFUVRxoWKoKIobY2KoKIobY2KoKIobY2KoKIobY2KoNKSGGMucPGvi40xi0vs7zXGLE2mgyra322MecwY05vY1mmMWVGhI7GScVQElZbDOf4OWGv7gV4kmqgjsb/weak7ZiNc7GsvLhmH2zYAXFjKgVhpXlQElVaku6iOdS9wWuL7WNNJLWSonnVBPAdrYJ+SIZo6bE5pfhJhjoX60t0MpcnqApa4HlgldBR9XwRczwhFuNwQtx/oKvT0XIzy8kTWlWJxVVoA7QkqWWAA6HRD0wFEiPoYEsUNuDm+C9znTve9I7F/o1yDLmxylTs+mXuycM4FyPC5D+h3YZcFFpPxfHjK+FARVFLFCd8shnpps4Alic/FCTB63XaQNGoXFsWGb5QrL3HeXCTr+GDRvkOQRKBd7pobzncJBbpLiavSGqgIKlmgMyFMnYnhbyfSg+soOn6pS5oxs9zCRoIOkISiuBIMJVjqjul31ytOxdYHzK9iWK40ASqCShYYGOFzd0EgXW9sOdJr3GiI6oSxL/G9C5ifSBO/pJD2CelNHmKM6XK9vY5CRTY2TtTbixNKpfXQLDJK05BMC+/891YkU8CXy3+nKCOhPUGlKXBD4jmJoXEvG/fYBhtoktIiaE9QaQncMHZA5+2USlERVBSlrdHhsKIobY2KoKIobY2KoKIobY2KoKIobY2KoKIobc3/B3njQ4a8h6JqAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 319.68x221.055 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "######\n",
    "fig, ax = std_fig(ax_form= [0.15,0.16,0.82,0.76])\n",
    "\n",
    "# HNL mass span\n",
    "x = np.logspace(np.log10(20),np.log10((m_neutral_kaon - m_mu)*1e3), 1000, endpoint=True)\n",
    "\n",
    "############################################\n",
    "# our PS191 bound (our own calculation)\n",
    "Xi,Yi,Zi=interp_grid(x=m4_grid_ps191*1e3, y=usqr_grid_ps191 , z=event_rate_weak_ps191_e_mu, logx=False, logy=True)\n",
    "level_0 = contour(Xi, Yi, Zi, levels = [2.3], colors=['royalblue'], linewidths=[0], zorder=10) #NEW CONSTRAINT\n",
    "bound_e_mu = level_0.collections[0].get_paths()[0].vertices  #grab the coords of the bound\n",
    "ax.plot(bound_e_mu[:,0],bound_e_mu[:,1], c='royalblue',  lw=1.5, zorder=10, rasterized=True, label='our mu e')\n",
    "\n",
    "############################################\n",
    "# our PS191 bound (our own calculation)\n",
    "Xi,Yi,Zi=interp_grid(x=m4_grid_ps191*1e3, y=usqr_grid_ps191 , z=event_rate_weak_ps191_mu_pi, logx=False, logy=True)\n",
    "level_0 = contour(Xi, Yi, Zi, levels = [2.3], colors=['royalblue'], linewidths=[0], zorder=10) #NEW CONSTRAINT\n",
    "bound_mu_pi = level_0.collections[0].get_paths()[0].vertices  #grab the coords of the bound\n",
    "ax.plot(bound_mu_pi[:,0],bound_mu_pi[:,1], c='royalblue',  lw=1.5, zorder=10, ls='--', rasterized=True, label='our mu pi')\n",
    "\n",
    "\n",
    "mn_mu_e, umu42_mu_e = np.genfromtxt(\"digitized/PS_191/umu4SQR_nu_mu_e.dat\", unpack=True)\n",
    "ax.plot(mn_mu_e, umu42_mu_e, lw=2, c='orange', label = r'Their mu e')\n",
    "\n",
    "mn_mu_pi, umu42_mu_pi = np.genfromtxt(\"digitized/PS_191/umu4SQR_mu_pi.dat\", unpack=True)\n",
    "ax.plot(mn_mu_pi, umu42_mu_pi, lw=2, c='orange', label = r'Their mu pi',ls='--')\n",
    "\n",
    "####\n",
    "# labels\n",
    "ax.set_title(r\"$|U_{e N}| = |U_{\\tau N}| = 0$\",fontsize=9)\n",
    "ax.legend(loc='upper center', markerfirst=False, frameon=False, handlelength=1.6,\n",
    "          framealpha=0.9, borderpad=0.2, edgecolor='None', facecolor='white', ncol=1, fontsize=8)\n",
    "\n",
    "ax.set_xlim(100,np.max(x))\n",
    "ax.set_ylim(1e-9,1e-2)\n",
    "ax.set_xscale(\"log\")\n",
    "ax.set_yscale(\"log\")\n",
    "ax.set_ylabel(r\"$|U_{\\mu N}|^2$\")\n",
    "ax.set_xlabel(r\"$m_{N}$/MeV\")\n",
    "\n",
    "# ax.set_xticks(np.array([20, 40, 60, 80, 100, 150, 200, 400]))\n",
    "ax.set_xticklabels(np.array([20,40,60,80,100,150, 200,400]))\n",
    "ax.grid(axis='y', which='both',dashes=(6,1),alpha=0.5,c='black',lw=0.1)\n",
    "ax.grid(axis='x', which='both',dashes=(6,1),alpha=0.5,c='black',lw=0.1)\n",
    "fig.savefig('plots/weak_decays_other_channels.pdf', dpi=300)\n",
    "fig.savefig('plots/weak_decays_other_channels.png', dpi=300)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "important-identification",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-60-d2e750ed4181>:4: RuntimeWarning: invalid value encountered in true_divide\n",
      "  mine(x)/theirs(x)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([       nan,        nan,        nan,        nan,        nan,\n",
       "              nan,        nan,        nan, 0.        , 1.64621533,\n",
       "       2.53138869, 3.8231285 , 4.56395446, 5.46785053, 5.66819018,\n",
       "       5.54826838, 5.90767926, 6.13144258, 6.4013163 , 6.649935  ,\n",
       "       6.70816448, 6.61379745, 6.64945735, 6.78346581, 6.72325039,\n",
       "       6.83419812, 6.74360731, 6.76414981, 6.70054417, 6.60023642,\n",
       "       6.50813433, 6.41843674, 6.26492069, 6.05241853, 5.80088093,\n",
       "       0.        , 0.        , 0.        ,        nan,        nan])"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mine = interpolate.interp1d(bound_e_mu[:,0], bound_e_mu[:,1], bounds_error= False, fill_value=np.inf)\n",
    "theirs=interpolate.interp1d(mn_mu_e, umu42_mu_e, bounds_error= False, fill_value=np.inf)\n",
    "x = np.linspace(40,400,40 )\n",
    "mine(x)/theirs(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "id": "behind-sherman",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "<ipython-input-61-0336af014098>:4: RuntimeWarning: invalid value encountered in true_divide\n",
      "  mine(x)/theirs(x)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([       nan,        nan,        nan,        nan,        nan,\n",
       "              nan,        nan,        nan,        nan,        nan,\n",
       "              nan,        nan,        nan,        nan,        nan,\n",
       "              nan,        nan,        nan,        nan,        nan,\n",
       "              nan,        nan,        nan, 0.        , 0.        ,\n",
       "       7.5818476 , 7.60518136, 7.58605793, 7.55368233, 7.53876833,\n",
       "       7.59437052, 7.61599989, 7.37800347, 6.97712841, 6.44558865,\n",
       "       0.        , 0.        , 0.        ,        nan,        nan])"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "mine = interpolate.interp1d(bound_mu_pi[:,0], bound_mu_pi[:,1], bounds_error= False, fill_value=np.inf)\n",
    "theirs=interpolate.interp1d(mn_mu_pi, umu42_mu_pi, bounds_error= False, fill_value=np.inf)\n",
    "x = np.linspace(40, 400, 40)\n",
    "mine(x)/theirs(x)"
   ]
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}