nuts_neutrino_experiments.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aca7bf6f-be63-40ef-94ff-d86196839c69",
   "metadata": {},
   "source": [
    "# NUTS experiments for neutrino NSI problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ba0bc164-77a6-4aa7-bf0a-c2f9a374b6a7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "import arviz as az\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import stats\n",
    "import pickle\n",
    "import pandas as pd\n",
    "import optax\n",
    "\n",
    "import json as js\n",
    "from tqdm import tqdm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ebea9a3c-458e-4a4d-bf79-a217792facda",
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.environ[\"CUDA_DEVICE_ORDER\"]=\"PCI_BUS_ID\"\n",
    "os.environ[\"CUDA_VISIBLE_DEVICES\"]=\"1\" #select GPU, -1 means use CPU\n",
    "\n",
    "from jax import config\n",
    "config.update(\"jax_enable_x64\", True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e8382b22-4332-4453-9e9a-a66d70e11f4c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "2024-05-19 18:36:49.602800: W external/xla/xla/service/gpu/nvptx_compiler.cc:679] The NVIDIA driver's CUDA version is 11.2 which is older than the ptxas CUDA version (11.8.89). Because the driver is older than the ptxas version, XLA is disabling parallel compilation, which may slow down compilation. You should update your NVIDIA driver or use the NVIDIA-provided CUDA forward compatibility packages.\n"
     ]
    }
   ],
   "source": [
    "from snuddd import config\n",
    "from snuddd.targets import Nucleus, electron_xe, Electron\n",
    "from snuddd.models import GeneralNSI\n",
    "from snuddd.jinterp import interp1d\n",
    "\n",
    "import numpyro\n",
    "from numpyro import optim\n",
    "from numpyro.diagnostics import hpdi\n",
    "import numpyro.distributions as dist\n",
    "from numpyro import handlers\n",
    "from numpyro import infer\n",
    "from numpyro.infer import MCMC, NUTS, autoguide, reparam\n",
    "\n",
    "import jax.numpy as jnp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6db32e11-fee5-40d9-a0fe-e6d3413e4e42",
   "metadata": {},
   "source": [
    "## Preparation for SNuDD "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6157ab0f-3cfa-4dfe-bda4-07cca4473d4d",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_nucleus = Nucleus(54, 132, mass=131.9041535 * config.u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "11862eba-36bd-4852-9b1c-42a7aa8e9384",
   "metadata": {},
   "outputs": [],
   "source": [
    "SM_matrix = np.array([[0, 0, 0],\n",
    "                      [0, 0, 0],\n",
    "                      [0, 0, 0]])\n",
    "\n",
    "SM_eta = 0\n",
    "SM_phi = 0\n",
    "\n",
    "sm_model = GeneralNSI(SM_matrix, SM_eta, SM_phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "cfde6386-4b0a-4b28-87c8-c27636989426",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_nucleus.update_model(sm_model)\n",
    "Xe_nucleus.prepare_density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "14828b57-a23c-4242-ac9c-2f540b647b46",
   "metadata": {},
   "outputs": [],
   "source": [
    "E_Rs = np.logspace(-1, 1, 100) / 1e6 # Recoil Energy in GeV \n",
    "\n",
    "sm_spec=Xe_nucleus.spectrum(E_Rs, nu = '8B')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "31809d25-9ca1-483a-a424-e97e77bd6669",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAHTCAYAAAD7+WZaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxfklEQVR4nO3deVhUZf8G8PvMwLAPuyAIiisoaoq4K2CKaZmaVrYQroW+Zcqrlm2mlZSWy1uSWQaaubRpWZaZu6L+3DDNFVMBAdn3dYbz+wOZGlkEzsDMwP25rrl0znnmPF+WcW6f85znCKIoiiAiIiIinZHpuwAiIiKi5oYBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIiIiIdIwBi4iIiEjHGLCIdEAQhDo9Dhw4gLfffhuCICA9PV3fZTdITEwM3n77bWRnZ2ttj46OhiAIuHnzZqP0W9/j//nnn5gyZQq8vLxgbm4Oa2tr9O7dG8uWLUNmZmaV49b2M2tInebm5rh161aV/YGBgfD19dXaduLECYwfPx6enp4wMzODi4sLBgwYgP/+97/16ruhavqZStUYvxNTp06FmZkZzp8/X2Xf+++/D0EQsHPnTp31R9RQJvougKg5OHbsmNbzd955B/v378e+ffu0tnft2rXeH9iGJiYmBosXL8bkyZNhZ2en2f7www/j2LFjaN26tf6Ku+vzzz/HrFmz0KVLF8yfPx9du3ZFWVkZTp06hbVr1+LYsWPYvn271muioqLg7e1d5Vhdu3ZtUA0lJSV444038NVXX9Xa7pdffsGjjz6KwMBALFu2DK1bt0ZycjJOnTqFrVu34qOPPmpQ//VR089Uqsb4nVi1ahX27t2L0NBQnDhxAqampgCA8+fPY9GiRZg8eTLGjBmjs/6IGooBi0gH+vfvr/Xc2dkZMpmsynZDVFhYCEtLS8nHcXZ2hrOzsw4qkubYsWOYOXMmRowYgR07dsDMzEyzb8SIEfjvf/+L3377rcrrfH190adPH53V8dBDD2Hz5s2YN28eevbsWWO7ZcuWwcvLC7t374aJyT//JE+aNAnLli1rUN+6+plK1Ri/E0qlEuvXr0dwcDDeffddLF68GGVlZQgJCYGLiwtWrVql0/6IGoqnCIn05M6dO3jqqadga2sLFxcXTJ06FTk5OVptrl27hqeffhqtWrWCmZkZfHx8sGbNmirHOnLkCB588EHY2NjA0tISAwcOxC+//FKlXeXpyTNnzmDixImwt7dHhw4d6tzf22+/jfnz5wMAvLy8tE6j1XQ66PLly3jqqafg4uICMzMzeHp64rnnnkNJSQkAIC4uDlOmTEGnTp1gaWkJd3d3jBkzptpTQHWxdOlSCIKAdevWaYWrSgqFAo8++mi9j7ts2TIIgoCrV69qtl27dg0mJibYtGlTlfYLFiyAo6MjXnnllVqPm5GRAScnJ61wVUkmu/8/0bX9TOvyva3tZ/rvr7Muv4f3qu53orLev/76676//zUZPnw4wsLCsHTpUpw+fRpvv/02zp07h/Xr18PW1larbUNrJ5KKAYtITyZMmIDOnTvj+++/x6uvvorNmzdj7ty5mv0XL16Ev78/Lly4gI8++gg///wzHn74YcyePRuLFy/WtDt48CCGDRuGnJwcrF+/Hlu2bIGNjQ3GjBmDbdu2Vdv3Y489ho4dO+Lbb7/F2rVr69zf9OnT8dJLLwEAfvjhBxw7dgzHjh1D7969q+3n3Llz8Pf3x/Hjx7FkyRL8+uuviIiIQElJCUpLSwEASUlJcHR0xPvvv4/ffvsNa9asgYmJCfr164crV67U63uqVquxb98++Pn5wcPDo96vValUWg+1Wq3ZX3m8+Ph4zbZOnTqhT58+2LFjR5Xj2djY4I033sDu3burnCr+twEDBuDEiROYPXs2Tpw4gbKysnrVXam6n2ldvrf3+5nW9fewvu73+38/y5cvh6enJyZOnIgPPvgAYWFhGDFihFabxqqdqE5EItK50NBQ0crKqtp9ixYtEgGIy5Yt09o+a9Ys0dzcXCwvLxdFURRHjhwptmnTRszJydFq9+KLL4rm5uZiZmamKIqi2L9/f7FVq1ZiXl6epo1KpRJ9fX3FNm3aaI73777feuutKnXVtb/ly5eLAMQbN25otYuKiqqyfdiwYaKdnZ2Ymppa7feiOiqVSiwtLRU7deokzp07t9bj3yslJUUEIE6aNKnO/VUet7qHXC7XtDt8+LAIQFy/fr3W68eNGyf26dOnyvFOnjwplpSUiO3btxf79Omj+TkEBASI3bp107RPT08XBw8erOnT1NRUHDhwoBgREaH1M61JbT/Te9X0va3pZyqKdf+9qE51P7O6/v7XxebNm0UAoqura7XfKym1E0nFESwiPbn3NFWPHj1QXFyM1NRUFBcXY+/evRg/fjwsLS21RlVGjx6N4uJiHD9+HAUFBThx4gQmTpwIa2trzbHkcjlCQkKQmJhY7SjQhAkTtJ7Xtb/6KCwsxMGDB/HEE0/UOg9HpVJh6dKl6Nq1KxQKBUxMTKBQKHDt2jVcunSpXn1KsXHjRpw8eVLrceLECc3+yhGshIQErdddunQJ7u7u1R5ToVDg3XffxalTp/DNN99U28bR0RGHDx/GyZMn8f7772Ps2LG4evUqFi5ciO7du9f5atN7f6aA9O9tY/xeVKrt978uysvL8fHHH0MmkyE1NRXnzp1rstqJ6oIBi0hPHB0dtZ5XzhcqKipCRkYGVCoVPv74Y5iammo9Ro8eDQBIT09HVlYWRFGs9iotNzc3ABVzfO51b/u69lcfWVlZUKvVaNOmTa3twsPD8eabb2LcuHHYuXMnTpw4gZMnT6Jnz54oKiqqV59OTk6wtLTEjRs36vU6APDx8UGfPn20Hn5+fpr97u7ukMlkWqcIv/32W1y5cgUhISE1HnfSpEno3bs3Xn/99VpP//Xp0wevvPIKvv32WyQlJWHu3Lm4efNmnSe6V/c7IPV72xi/F5Vq+/2viw8//BDHjh3D5s2b0alTJ0ydOlXrtY1ZO1Fd8CpCIgNkb2+vGYX6z3/+U22byvWdZDIZkpOTq+xPSkoCUBE67iUIQoP6qw8HBwfI5XIkJibW2m7Tpk147rnnsHTpUq3t6enp9V4yQC6X48EHH8Svv/6KxMTE+4a7+jAxMYGrq6smYMXFxWHmzJkYM2ZMtaNHlQRBwAcffIARI0Zg3bp1derL1NQUixYtwsqVK3HhwoU6vebenykg/XvbGL8XunDx4kW89dZbeO655/Dkk0+ibdu2GDRoEF5//XWsWLECgOHWTi0HAxaRAbK0tERQUBDOnj2LHj16QKFQ1Ni2X79++OGHH/Dhhx/CwsICQMXpk02bNqFNmzbo3LmzTvur60iDhYUFAgIC8O233+K9996rNugBFcHg3qv9fvnlF9y+fRsdO3a8b+33WrhwIXbt2oUZM2bgxx9/rPK1lJWV4bfffmvQWkkeHh5ISEhAcnIyRo8eja5du2LLli33fd3w4cMxYsQILFmypMrk++Tk5GpHnypP4VWORDZEXb+3Nf1M6/N70VRUKhVCQ0Ph5OSE1atXA6hYJiU8PBwrVqzAhAkTMGjQIIOsnVoWBiwiA7V69WoMHjwYQ4YMwcyZM9GuXTvk5eUhLi4OO3fu1FyZFhERgREjRiAoKAjz5s2DQqFAZGQkLly4gC1btlQ7siGlv+7du2vah4aGwtTUFF26dKn2mCtWrMDgwYPRr18/vPrqq+jYsSPu3LmDn376CZ999hlsbGzwyCOPIDo6Gt7e3ujRowdOnz6N5cuXN3j0acCAAfj0008xa9Ys+Pn5YebMmejWrRvKyspw9uxZrFu3Dr6+vlUC1oULF6BSqaocr0OHDpo5ZB4eHjh37hyGDBkCT09P7Ny5E1ZWVnWq64MPPoCfnx9SU1PRrVs3zfaRI0eiTZs2GDNmDLy9vVFeXo7Y2Fh89NFHsLa2xssvv9yg7wOAOn9va/qZ2tjY1Pn3oqlERETg1KlT+PXXX7VG4d555x3s3LkTU6dORWxsLCwsLAyudmph9D3Lnqg5qstVhGlpaVrbq7vi6saNG+LUqVNFd3d30dTUVHR2dhYHDhwovvvuu1qvPXz4sDhs2DDRyspKtLCwEPv37y/u3Lmzzn3Xt7+FCxeKbm5uokwmEwGI+/fvr/Eqv4sXL4qPP/646OjoKCoUCtHT01OcPHmyWFxcLIqiKGZlZYnTpk0TW7VqJVpaWoqDBw8WDx8+LAYEBIgBAQG1fn9qExsbK4aGhoqenp6iQqEQraysxF69eolvvfWW1lWNtV1FCED8/PPPNW3nzp0rAhCffvppTf33+vdVhPd6+umnRQBaVxFu27ZNfPrpp8VOnTqJ1tbWoqmpqejp6SmGhISIFy9evO/XWdvPtK7fW1Gs/mdaqa6/FzV9L6q7irAuv//3io2NFU1NTcUZM2ZUu//YsWOiTCbTukKyobUTSSWIoig2aaIjIiIiauZ4FSERERGRjjFgEREREekYAxYRERGRjjFgEREREekYAxYRERGRjjFgEREREekYFxrVgfLyciQlJcHGxqbOizoSERGRfomiiLy8PLi5uUEm0+2YEwOWDiQlJVW5/QUREREZh4SEBJ3euxRgwNIJGxsbABU/IKVSqedqiIiIqC5yc3Ph4eGh+RzXJQYsHag8LahUKhmwiIiIjExjTO/hJHciIiIiHWPAIiIiItIxBiwiIiIiHeMcLCIiA6RWq1FWVqbvMoiMmqmpKeRyuV76ZsAiIjIgoigiJSUF2dnZ+i6FqFmws7ODq6trk69TyYBFRGRAKsNVq1atYGlpycWLiRpIFEUUFhYiNTUVANC6desm7Z8Bi4jIQKjVak24cnR01Hc5REbPwsICAJCamopWrVo16elCTnInIjIQlXOuLC0t9VwJUfNR+X5q6jmNDFhERAaGpwWJdEdf7ycGLCIiIiIdY8AiIiIyEtHR0bCzszO6Y9fH5MmTMW7cuDq3P3DgAARBMLgrbxmwiIhIstTUVLzwwgvw9PSEmZkZXF1dMXLkSBw7dkzTpl27dhAEAcePH9d67Zw5cxAYGFiv/gIDAyEIAgRBgJmZGdzd3TFmzBj88MMPVdpWthMEAVZWVujUqRMmT56M06dP37efs2fP4pFHHkGrVq1gbm6Odu3a4cknn0R6ejoAw/1wbyyV38d7f4YlJSVwdHSEIAg4cOCAfoozMAxYREQk2YQJE3Du3Dls2LABV69exU8//YTAwEBkZmZqtTM3N8crr7yikz5nzJiB5ORkxMXF4fvvv0fXrl0xadIkPP/881XaRkVFITk5GX/99RfWrFmD/Px89OvXDxs3bqzx+KmpqRg+fDicnJywe/duXLp0CV9++SVat26NwsJCnXwN+lBaWirp9R4eHoiKitLatn37dlhbW0s6bnPDgEVERJJkZ2fjyJEj+OCDDxAUFIS2bduib9++WLhwIR5++GGtti+88AKOHz+OXbt21Xi8wMBAzJkzR2vbuHHjMHnyZK1tlpaWcHV1hYeHB/r3748PPvgAn332GT7//HP88ccfWm0rF5ts164dgoOD8d133+GZZ57Biy++iKysrGrriImJQW5uLr744gv06tULXl5eGDZsGFatWgVPT0/cvHkTQUFBAAB7e3sIgqCp8bfffsPgwYNhZ2cHR0dHPPLII7h+/brm2Ddv3oQgCPjhhx8QFBQES0tL9OzZU2vED6g4befp6QlLS0uMHz8eGRkZWvuvX7+OsWPHwsXFBdbW1vD396/ytbdr1w7vvvsuJk+eDFtbW8yYMaNOx65JaGgotm7diqKiIs22L7/8EqGhoVXanj9/HsOGDYOFhQUcHR3x/PPPIz8/X7NfrVYjPDxc831asGABRFHUOoYoili2bBnat28PCwsL9OzZE999912datUnBiwiIgMmiiIKS1V6edz7QVcTa2trWFtbY8eOHSgpKam1bbt27RAWFoaFCxeivLxcF98iLaGhobC3t6/2VOG95s6di7y8POzZs6fa/a6urlCpVNi+fXu13wsPDw98//33AIArV64gOTkZq1evBgAUFBQgPDwcJ0+exN69eyGTyTB+/PgqX/Prr7+OefPmITY2Fp07d8ZTTz0FlUoFADhx4gSmTp2KWbNmITY2FkFBQXj33Xe1Xp+fn4/Ro0fjjz/+wNmzZzFy5EiMGTMG8fHxWu2WL18OX19fnD59Gm+++Wadjl0TPz8/eHl5ab72hIQEHDp0CCEhIVrtCgsL8dBDD8He3h4nT57Et99+iz/++AMvvviips1HH32EL7/8EuvXr8eRI0eQmZmJ7du3ax3njTfeQFRUFD799FP89ddfmDt3Lp599lkcPHiwTvXqCxcaJSIyYEVlanR9a7de+r64ZCQsFff/mDAxMUF0dDRmzJiBtWvXonfv3ggICMCkSZPQo0ePKu0rPzC//vrrKh/KUslkMnTu3Bk3b968b1tvb28AqLFt//798dprr+Hpp59GWFgY+vbti2HDhuG5556Di4sL5HI5HBwcAACtWrXSmiA+YcIErWOtX78erVq1wsWLF+Hr66vZPm/ePM0o3+LFi9GtWzfExcXB29sbq1evxsiRI/Hqq68CADp37oyYmBj89ttvmtf37NkTPXv21Dx/9913sX37dvz0009aQWbYsGGYN2+e5vlbb71132PXZsqUKfjyyy/x7LPPIioqCqNHj4azs7NWm6+//hpFRUXYuHEjrKysAACffPIJxowZgw8++AAuLi5YtWoVFi5cqPl+rV27Frt3//P7XlBQgBUrVmDfvn0YMGAAAKB9+/Y4cuQIPvvsMwQEBNSpXn3gCBYREUk2YcIEJCUl4aeffsLIkSNx4MAB9O7dG9HR0VXaOjs7Y968eXjrrbckzweqjiiKdVr7qHJUqra27733HlJSUrB27Vp07doVa9euhbe3N86fP1/rsa9fv46nn34a7du3h1KphJeXFwBUGVn6dwCtvJVL5a1dLl26pAkVle59XlBQgAULFqBr166ws7ODtbU1Ll++XKWfPn36aD2vy7Fr8+yzz+LYsWP4+++/ER0djalTp1Zpc+nSJfTs2VMTrgBg0KBBKC8vx5UrV5CTk4Pk5GStfk1MTLRqvXjxIoqLizFixAjNSKm1tTU2btyodcrVENVrBKt9+/aSO5wzZw5mz54t+ThERC2BhakcF5eM1Fvf9WFubo4RI0ZgxIgReOuttzB9+nQsWrSoytwpAAgPD0dkZCQiIyOr7JPJZFVOydV1FW61Wo1r167B39//vm0vXboEAJrwUxNHR0c8/vjjePzxxxEREYFevXrhww8/xIYNG2p8zZgxY+Dh4YHPP/8cbm5uKC8vh6+vb5VAaWpqqvl7ZdCrPI1Yl1O08+fPx+7du/Hhhx+iY8eOsLCwwMSJE6v08++QU9dj16ZyXtm0adNQXFyMUaNGIS8vr0ofNYXXui7+Wfm9+OWXX+Du7q61z8zMrAGVN516jWDdvHkTWVlZEEWxQY9bt261mEtZiYh0QRAEWCpM9PKQugJ2165dUVBQUO0+a2trvPnmm3jvvfeQm5urtc/Z2RnJycma52q1GhcuXKhTnxs2bEBWVlaVU3TVWbVqFZRKJYYPH16nYwOAQqFAhw4dNF+XQqHQ1FgpIyMDly5dwhtvvIEHH3wQPj4+NU6kr03Xrl2rLIdw7/PDhw9j8uTJGD9+PLp37w5XV9c6nR6ty7HvZ+rUqThw4ACee+65au/x17VrV8TGxmr9Dhw9elRzGtfW1hatW7fW6lelUmktn9G1a1eYmZkhPj4eHTt21Hp4eHjUq96mVu85WHPnzsVbb73VoM5kMp6RJCJqbjIyMvD4449j6tSp6NGjB2xsbHDq1CksW7YMY8eOrfF1zz//PFauXIktW7agX79+mu3Dhg1DeHg4fvnlF3To0AErV66s9j/nhYWFSElJgUqlwu3bt/HDDz9g5cqVmDlzpubqvkrZ2dlISUlBSUkJrl69is8++ww7duzAxo0ba1xc8+eff8bWrVsxadIkdO7cGaIoYufOndi1a5dmmYK2bdtCEAT8/PPPGD16NCwsLGBvbw9HR0esW7cOrVu3Rnx8vGauU33Mnj0bAwcOxLJlyzBu3Dj8/vvvVeZIdezYET/88APGjBkDQRDw5ptv1unigboc+34eeughpKWlQalUVrv/mWeewaJFixAaGoq3334baWlpeOmllxASEgIXFxcAwMsvv4z3338fnTp1go+PD1asWKH1s7axscG8efMwd+5clJeXY/DgwcjNzUVMTAysra2rvXLRYIj1IAiCuHjx4vq8RKevN1Q5OTkiADEnJ0ffpRCRESsqKhIvXrwoFhUV6buUeikuLhZfffVVsXfv3qKtra1oaWkpdunSRXzjjTfEwsJCTbu2bduKK1eu1Hrt5s2bRQBiQECAZltpaak4c+ZM0cHBQWzVqpUYEREhjh07VgwNDdW0CQgIEAGIAESFQiG2bt1afOSRR8QffvihSn2V7QCI5ubmYocOHcTQ0FDx9OnTtX5d169fF2fMmCF27txZtLCwEO3s7ER/f38xKipKq92SJUtEV1dXURAETY179uwRfXx8RDMzM7FHjx7igQMHRADi9u3bRVEUxRs3bogAxLNnz2qOk5WVJQIQ9+/fr9m2fv16sU2bNqKFhYU4ZswY8cMPPxRtbW01+2/cuCEGBQWJFhYWooeHh/jJJ5+IAQEB4ssvv1zr970ux67Ov7+Ge1VX/59//ikGBQWJ5ubmooODgzhjxgwxLy9Ps7+srEx8+eWXRaVSKdrZ2Ynh4eHic889J44dO1bTpry8XFy9erXYpUsX0dTUVHR2dhZHjhwpHjx4UBRFUdy/f78IQMzKyqq2rtreV435+S2IYt1PxF6/fh0ODg6wt7dvUJiT+npDlZubC1tbW+Tk5NSY5ImI7qe4uBg3btyAl5cXzM3N9V0OUbNQ2/uqMT+/63WKsEOHDpI6k/p6IiIiImPASVFEREREOsaARURERKRjklZy37t3L/bt24eYmBgkJiYiPT0dlpaWcHZ2Rvfu3REQEIBHHnkErq6uuqq30ZiYmGhW1+3Tpw+++OILPVdERERExqreASs/Px//+9//8PnnnyM+Pl6zWJm5uTkcHBxQVFSECxcu4M8//8TXX38NExMTPProo5g7dy4GDRqk8y9AV+zs7BAbG6vvMoiIiKgZqNcpwrVr16Jjx4544403YGdnh3fffRf79u1Dbm4uCgsLkZiYiIyMDJSVleHy5cvYsGEDnnzySfz+++8YOnQoHnvsMdy4caOxvhYiomahMW6CTNRS6ev9VK9lGkxNTfHMM89g/vz56NatW507KSoqwpYtWxAREYGQkJAGL1Rak0OHDmH58uU4ffo0kpOTsX37dowbN06rTWRkJJYvX47k5GR069YNq1atwpAhQzT7FQoFunfvDgsLC7z33nv1uoEkl2kgIl0oLy/HtWvXIJfL4ezsDIVCIXk1daKWShRFlJaWIi0tDWq1Gp06daqy4LnBLNNw+fLlBi21YGFhgalTpyI0NBSJiYn1fv39FBQUoGfPnpgyZUq1t0fYtm0b5syZg8jISAwaNAifffYZRo0ahYsXL8LT0xNAxW2A3NzccOHCBTz88MM4f/48wxIRNSmZTAYvLy8kJycjKSlJ3+UQNQuWlpbw9PRs8rvJ1GsEC6i4T5CJiaS58Y1KEIQqI1j9+vVD79698emnn2q2+fj4YNy4cYiIiKhyjFGjRuGdd96pcvfxSiUlJSgpKdE8z83NhYeHB0ewiEgnRFGESqXSur8dEdWfXC6HiUnN99U0mBEsAHBzc8Nzzz2HqVOnomvXrjotpjGUlpbi9OnTVe4DFRwcjJiYGABAVlYWLC0tYWZmhsTERFy8eBHt27ev8ZgRERFYvHhxo9ZNRC2XIAgwNTWFqampvkshogaq93hZTk4OVqxYge7du2PgwIFYv3498vPzG6M2nUhPT4dardbcWLKSi4sLUlJSAACXLl1Cnz590LNnTzzyyCNYvXo1HBwcajzmwoULkZOTo3kkJCQ06tdARERExqXeI1jJycn46quvEBUVhePHj+PEiROYM2cOnnjiCUyZMgWDBw9ujDolu3d4UBRFzbaBAwfi/PnzdT6WmZkZzMzMdFofERERNR/1HsFycHDAyy+/jNjYWJw6dQozZ86EQqFAVFQUAgIC4OPjg+XLl+POnTuNUW+9OTk5QS6Xa0arKqWmplYZ1SIiIiLSBUlT6nv37o1PPvkEycnJ2Lx5Mx588EFcu3YNr776Kjw8PDB+/Hj8/PPPel3TRaFQwM/PD3v27NHavmfPHgwcOFBPVREREVFzppPLARUKBSZNmoRJkyYhMTERUVFRiI6Oxo8//oiffvoJrq6uuH37ti66qlZ+fj7i4uI0z2/cuIHY2Fg4ODjA09MT4eHhCAkJQZ8+fTBgwACsW7cO8fHxCAsLa7SaiIiIqOWq9zINdXXlyhVMnz4dR48ehSAIjXq58YEDBxAUFFRle2hoKKKjowFULDS6bNkyJCcnw9fXFytXrsTQoUN10j8XGiUiIjI+jfn5rdOAVVBQgG+++QZffvklYmJiIIoiLC0tMXHiRE3QaY4YsIiIiIyPQa2DVZ3Dhw/jyy+/xHfffYfCwkKIogh/f39MmzYNTz31FGxsbHTRDREREZFRaHDAun37NjZs2IDo6Ghcv34doijC0dER06dPx7Rp0+Dr66vLOomIiIiMRr0D1jfffIOoqCj88ccfUKvVkMlkCA4OxtSpUzFu3DiuPExEREQtXr0D1qRJkwAA7dq1w5QpUzBlyhS0adNG54URERERGasGBaxp06bhwQcfbIx6iIiIiIxevQPW5s2bG6MOIiIiomZD0kruAKBSqbBy5Ur07dsXSqUSJib/ZLbY2FjMmjULV69eldoNERERkdGQtExDUVERgoODERMTAycnJyiVShQUFGj2e3l5ISoqCg4ODnj33XclF0tERERkDCSNYC1duhRHjx5FREQEUlJSMH36dK39tra2CAgIwO7duyUVSURERGRMJAWsbdu2ITAwEAsWLIAgCBAEoUqb9u3bIz4+Xko3REREREZFUsCKj4+Hv79/rW2USiVycnKkdENERERkVCQFLBsbG6SlpdXa5vr163B2dpbSDREREZFRkRSw+vfvj507d9Y4QpWYmIhdu3Zh6NChUrohIiIiMiqSAtb8+fORmZmJ4cOHIyYmBiqVCgBQWFiIvXv3Ijg4GGVlZQgPD9dJsURERETGQNIyDUOHDsWaNWswe/ZsDBkyRLPdxsYGACCXyxEZGQk/Pz9pVRIREREZEUEURVHqQS5duoS1a9fixIkTyMzMhFKpRL9+/TBr1ix069ZNF3UatNzcXNja2iInJwdKpVLf5RAREVEdNObnt04C1v2Ul5dDJpO8aLzBYsAiIiIyPo35+S0p9XzxxRf3baNWq/H0009L6YaIiIjIqEgKWDNnzsSPP/5Y435RFPHMM8/g22+/ldINERERkVGRvEzDU089hcOHD1fZVxmuvvnmG4SFhUnphoiIiMioSApYP//8Mzp06IBHH30U58+f12wXRREhISHYunUrXnjhBaxZs0ZyoURERETGQlLAsrW1xe7du2Fra4uHHnoIN2/ehCiKePbZZ7F582Y8//zz+PTTT3VVKxEREZFRkHxpn5ubG37//XeUlpYiODgYTz31FLZs2YLp06dj7dq1uqiRiIiIyKjoZO2Ezp0749dff0VKSgq+/fZbTJs2DevWrdPFoYmIiIiMTr1Wcl+yZEmt+/v27YvY2Fi4u7trtRUEAW+++WbDKiQiIiIyMvVaaLShi4UKggC1Wt2g1xoDLjRKRERkfBrz87teI1j79+/XaedEREREzVG9AlZAQEBj1UFERETUbDTfGwQSERER6Um9RrBqEx8fDwBwdHSElZWVrg5LREREZHR0NoLVrl07tG/fHl9//bWuDklERERklHQ2guXp6QlBEGBjY6OrQxIREREZJZ0FrJs3b+rqUERERERGjZPciYiIiHRMUsCSy+V45plndFULERERUbMgKWAplUp4eHjoqhYiIiKiZkFSwOrbty/OnTunq1qIiIiImgVJAWvx4sXYt28fNmzYoKt6iIiIiIyepKsIf//9dwQGBmLq1Kn4+OOP0bdvX7i4uEAQBK12giDgzTfflFQoERERkbEQRFEUG/pimaxuA2CCIECtVje0G4PXmHfjJiIiosbRmJ/fkkaw9u/fr6s6iIiIiJoNSQErICBAV3UQERERNRtcaJSIiIhIx3R2qxy1Wo309HSUlJRUu9/T01NXXREREREZNMkB6/Tp03jttddw6NAhlJaWVttGEASoVCqpXREREREZBUkBKzY2FkOGDIGJiQmCg4Oxc+dO9OzZE66urjhz5gzS0tIQGBiItm3b6qpeIiIiIoMnaQ7WO++8AwA4ceIEfvzxRwDA+PHj8euvv+LmzZsICwvDhQsXsGjRIumVEhERERkJSQHryJEjePTRR+Hj46PZVrmsloWFBT755BO4ubnhtddek1YlERERkRGRFLBycnLQvn17zXNTU1Pk5+f/c3CZDIGBgdi7d6+UboiIiIiMiqSA1apVK2RlZWmeu7q64tq1a1ptiouLUVhYKKUbIiIiIqMiaZJ7165dceXKFc3zQYMGYceOHTh+/Dj69++PS5cu4ZtvvoG3t7fkQo1B+LazUCptYWYqh7mpDOamcpib/PN3C1M5zExlsDCVVzxXyKv83cJUDjMTGWQy4f4dEhERkUGSFLAefvhhzJ07F8nJyWjdujVeeeUVbN++HYMGDYKDgwOysrJQXl7eYuZg/X4xFTKz/Ps3rAPzu0HMUmECc1MZLBUmmhBmqagIZJaKiv0WpnJYmclhoTCB1b+2W5nd/VNhAkszOazuHuvem3ETERGRbkm62XNZWRkyMzNhb28PhUIBAIiJicF7772Hv//+G23btsVLL72Ehx9+WGcFG6LKm0Wu3fMnBIUlisvKUVymrvhTpUZxqbriz7vbi8rUKCpVo0RVjqK7+yqfNzZBQEXgUshhbfZP8LI2M4HV3Ye1mRzWZqawMpPDxrxymwlszE3+td0U1mYmkHOkjYiIjFRj3uxZUsBqLvLy8jBs2DCUlZVBrVZj9uzZmDFjRp1fr6sfkLpcRIlKjcLSisBVXHb373cDWWGpGoWlKhSVVf5djaJSlaZ9YakaBaUqFJWqUXB3X0GpGgUlFW0ag5VCDmtzE03gsjE3gdLcFDbm//xdaWEKpYUJbMz++buthSmU5qawVMg5okZERHrBgNXI1Go1SkpKYGlpicLCQvj6+uLkyZNwdHSs0+sb8wekK+XlIorKKgJYYYka+SUqFJSoUFCqQkFJRQir2FbRJr9EhfxilWa75lGsQl6xCqVq3Yy2mciEitBlfjd0WZjC9u7DztIUdhaKiueWprCzMIWdpQL2lhXPzUzkOqmBiIhapsb8/JZ8qxyVSoWPP/4YW7ZsweXLl1FYWKi5LU5sbCzWrVuHOXPmoHPnzpKLbSxyuRyWlpYAKq56VKvVaG65UyYTNKcAYSP9eCUqNfKLK0JXXrEKucVlmvCVW1yGvGIV8u7+Wfk8t6gMuXf/zCkqg6pchKpcRGZBKTILqr/NUm0sFXJN6HKwUsDO0hT2lgrYWyngYGla8aeVAvZ39ztYKWBuylBGRESNT1LAKioqQnBwMGJiYuDk5ASlUomCggLNfi8vL0RFRcHBwQHvvvuu5GJrcujQISxfvhynT59GcnIytm/fjnHjxmm1iYyMxPLly5GcnIxu3bph1apVGDJkiGZ/dnY2AgICcO3aNSxfvhxOTk6NVm9zYGYih5m1HI7WZg16vShWjKjlFqmQU1SG3OIy5BRWBK+cojJkF5Vpglh2YSmyi8qQXfjP83IRmtOkSTnFde7XUiGHg5UCjncDl6O1GRytFXCyqvjT0doMjlYKONuYwcFKAVO5pJVMiIiohZIUsJYuXYqjR4/i/fffx/z587F48WLN7XMAwNbWFgEBAdi9e3ejBqyCggL07NkTU6ZMwYQJE6rs37ZtG+bMmYPIyEgMGjQIn332GUaNGoWLFy/C09MTAGBnZ4dz587hzp07eOyxxzBx4kS4uLg0Ws0tnSAIsFSYwFJhAldb83q9trxcRF6JCtmFpcgqLENWYSmyC0uRWVB2989SZBeWIbOgFFmFpcgoKEVWQSlU5eLdUFaExKyiOvVlb2kKJ2szOFmbwdnmn0fl81Z3H/aWCi6tQUREGpIC1rZt2xAYGIgFCxYAQLWTldu3b4+zZ89K6ea+Ro0ahVGjRtW4f8WKFZg2bRqmT58OAFi1ahV2796NTz/9FBEREVptXVxc0KNHDxw6dAiPP/54tccrKSlBSUmJ5nlubq4OvgqqK5lM0MzTalu3aXIQRRG5xaq7pyNLkJFfEcQyCkqRkV+KjIISZBaUIi2vBBl3T1mqy8W7Aa4M11JrX37DVC7AydoMrZTmcLExg4vSHC7Ku8+V5nC9+1BamHBSPxFRCyApYMXHx2P8+PG1tlEqlcjJyZHSjSSlpaU4ffo0Xn31Va3tlac2AeDOnTuwsLCAUqlEbm4uDh06hJkzZ9Z4zIiICCxevLhR6ybdEoR/QpmXk9V925eXi8gqLEV6finS80uQnl+CtLx/Pe4+T82rCGZlahHJOcVIvs/pSnNTWUXYsjWHm60FXG3N0drOAq2V5mhtZw53OwvYWpgyhBERGTlJAcvGxgZpaWm1trl+/TqcnZ2ldCNJeno61Gp1ldN9Li4uSElJAQAkJiZi2rRpEEURoijixRdfRI8ePWo85sKFCxEeHq55npubCw8Pj8b5AkgvZDLh7vwsM3S5z1UBpapypOdXhK3U3GLcqfwztxgpuSW4k1OMlNxi5BSVobisHDczCnEzo+bbR1mYyuFmZw43Owu42VrA3d4C7nb//Nna1hwmnBtGRGTQJAWs/v37Y+fOncjJyYGtrW2V/YmJidi1a1eVCef6cO+IgCiKmm1+fn6IjY2t87HMzMxgZtawyd3U/ChMZBVhyM6i1nZFpWrcya0Y5ar8Mzmn6J8/s4uRUVCKojI1rqcV4HpaQbXHkcsEuCrN4eFggTb2lmhjbwEPe0t4OlrC08ESztZmnA9GRKRnkgLW/PnzERQUhOHDh2P16tWa5RkKCwtx7NgxvPTSSygrK9Ma7WlqTk5OkMvlmtGqSqmpqZzETk3KQiFHOycrtKvlFGVxmboicGUXITG7CEl3H7ezi3A7qwhJ2cUoVZdXPM8uApBZ5RhmJjJ4OFSELU8HS7R1tEQ7Ryt4OlrCw94SChOOfhERNTZJAWvo0KFYs2YNZs+erbXkgY1NxSkVuVyOyMhI+Pn5SatSAoVCAT8/P+zZs0drvtiePXswduxYvdVFVB1zUzm8nKxqnCdWXi4iPb8ECVmFSMwqQkJmIRIyi5CYXYj4zEIkZRejRFWOuNR8xFUzMV8mAO72FmjnWNFHO0creDlboYOTNdztLXjrIyIiHZG80GhYWBgCAgKwdu1anDhxApmZmVAqlejXrx9mzZqFbt266aLOWuXn5yMuLk7z/MaNG4iNjYWDgwM8PT0RHh6OkJAQ9OnTBwMGDMC6desQHx+PsLCwRq+NSJdkMgGtlOZopTSHX9uq+8vU5UjOLkZ8ZkXgupVRgFsZhbiZUYD4zEIUlqqRkFmEhMwiHL6WrvVahVyGto6WaO9shfbO1ujobI0OrazRwdkKNuamTfQVEhE1D83iVjkHDhxAUFBQle2hoaGIjo4GULHQ6LJly5CcnAxfX1+sXLkSQ4cO1Un/xnCrHCJRFJGWV4Ib6QW4mVGAmxmFuJFWgBvpBbiRUYDSWm427qI0Q8dW1ujUygadXKzR2cUGnVvZwNaSwYuIjBfvRWjgGLDI2JWXi7idXYS/0wvwd1o+rqfl43pqAeLS8pGWV1Lj61rZmKGLqw26uNigs6sNvF1t0KmVDSwUvCURERk+gw1Y48ePx/PPP4+HHnqoRa/bw4BFzVlOURmup+Uj7k4+rqXm4VpqPq7dyb87yb4qQQC8HK3g01oJn9Y28HZVoqubEq1tzVv0vxNEZHgMNmDJZDIIggB3d3dMnToVU6dO1dx6piVhwKKWKL9EhWt38nAlJQ9XKv9MyUNGDTfutrc0RVc3Jbq52aLb3T/bO1lxSQki0huDDViXLl3C559/jk2bNiE9PR0ymQwjRozAjBkz8Oijj8LERPIceqPAgEX0j9S8YlxOzsOl5Ny7jzzEpeVDXV71nxprMxN0dVOih7sturexRXd3W3g5WXGki4iahMEGrEplZWXYsWMH1q9fjz/++AOiKMLJyQmTJ0/GtGnT0LlzZ13UarAYsIhqV1ymxtU7efgrKRd/JeXgr6RcXEzKRUk1E+uV5ibo6WGHBzzs0LONHXp62MHZhgv7EpHuGXzA+reEhASsX78e0dHRSEhIAAAMGTIEM2bMwMSJE5vlCugMWET1p1KXIy4tH+cTc3D+dg7+TMzBxeTcaq9m9HCwQG9Pe/T2tEcvTzv4tFbClLcLIiKJjCpgAYBarcaOHTvw8ssvIykpqaIjQYCDgwNeeeUVhIeHQyZrPv84MmAR6UapqhxX7+QhNiEb5xKyEZuQjbi0fNz7r5SFqRwPeNjBv509/No5oLenHdfqIqJ6M5qAFRcXhy+++AIbN27EnTt3oFAoMGHCBISGhuLs2bNYs2YNEhMTER4ejuXLl+uqW71jwCJqPLnFZTiXkI0zt7JxJj4LZ+OzkFus0mojEwBvVyX6tXdAPy9H9PNygL2VQk8VE5GxMOiAVVJSgm+//RZffPEFDh8+DFEU0aVLF8yYMQOTJ0+Gg4ODVtsRI0bgypUruHPnjuTiDQUDFlHTKS8XcT0tHydvZuHUrUycupmF+MzCKu26uNigf3sHDOjghP7tHWBnycBFRNoa8/Nb0mV+L730Er7++mvk5ORAoVBg0qRJeP755xEQEFBtezMzM4wcORJHjx6V0i0RtWAymYBOLjbo5GKDp/tVLAuTmluMEzcyceJGBk78nYlrqfkVS0fcycOGY7cgCEA3NyUGdXDCwI5O6NvOgYuhElGjkrwOVufOnTWjVY6Ojvd9zdGjR/HHH39g0aJFDe3W4HAEi8iwZOSX4P9uZOLY3xmIuZ5R5cbXChMZ/NvZY0gnZwzp5AQfVyXX4yJqgQz2FOGBAwcQGBiow3KMEwMWkWFLzS3Gsb8zcDQuHUeupSMpp1hrv5O1GQI6OyPI2xlDOjnD1oIT5olaAoMNWFSBAYvIeIiiiL/TC3D4ahoOX0vHsb8zUFiq1uyXywT4edoj0NsZI3xc0LGVNRc+JWqmjCJgxcfHAwAcHR1hZWWli0MaDQYsIuNVqirHqVuZOHAlDfsvp+LaPacT2zpaYriPC4b7uMC/nT1MuP4WUbNhFAFLJpNBJpMhMjISzz//vC4OaTQYsIiaj4TMQhy4koq9l1MRE5eBUvU/C5/aWZpiuI8LHurmisGdnGBuyonyRMbMYK8i/DdPT08IggAbGxtdHZKIqMl5OFgiZEA7hAxoh4ISFQ5fS8Oei6nYd/kOsgrL8N3pRHx3OhFWCjmG+bhglK8rgrq04lWJRKSFc7B0gCNYRM2fulzEyZuZ+O1CCn67kIKU3H8mylsq5Bju44JHerRGQBdnmJkwbBEZA6M4RdiSMWARtSzl5SLOJWbjtwsp+OV8MhKzijT7bMxMENzNFY/1dkf/9o6Qc/kHIoPFgGXgGLCIWi5RFBGbkI2f/0zGL38ma41suSrNMfYBN4zv7Q5vV/7bQGRoDDpglZaWYseOHTh58iSys7OhVqurtBEEAevXr5fSjUFjwCIioGJk69StLOyIvY1f/kxGTlGZZl/X1ko80acNxj7gzvskEhkIgw1Yt27dwogRI3D9+nXUdhhBEKoNXs0FAxYR3atEpcb+y6nYfvY29l1ORZm64t9IhVyGEd1c8GQfDwzq6MRTiER6ZLAB67HHHsOOHTsQEhKCqVOnok2bNjAxqf7CxLZt2za4SEPHgEVEtckuLMWPsUnYdjIBF5NzNdvd7SzwVF8PPOHvgVY25nqskKhlMtiAZWdnB39/f+zZs0eXNRkdBiwiqqsLt3PwzakE7Dh7G7nFKgCAiUzASF9XPNuvLfq3d+DK8URNxGDXwSovL0evXr10VQsRUbPn624LX3dbvDbaB7vOJ2PT8Vs4E5+NX+5Oku/YyhqTB7bDY73dYanQ2VKFRNTEJI1gjRw5EgqFAjt37tRlTUaHI1hEJMXFpFxsOnELO87e1twX0dbCFE/19cRzA9rCzc5CzxUSNU8Ge4rw7NmzGDJkCKKjozFx4kRd1mVUGLCISBfyiitWio86ehPxmYUAKm4+PcrXFWEBHeDrbqvnComaF4MNWEuWLMHJkyexa9cuBAQEoFevXrC1rfoPgCAIePPNNyUVasgYsIhIl9TlIvZdTsWXR27g2N8Zmu1DOjkhLKADBnZw5DwtIh0w2IAlk9XtrvJcpoGIqGEuJuXi88N/46dzSVCXV/xz3d3dFjMDO+Chbq6QcZkHogYz2IB18ODBOrcNCAhoaDcGjwGLiBpbQmYhvjj8N7adSkBxWTkAoLOLNV5+sDNG+TJoETWEwQYsqsCARURNJSO/BNExNxEdcxN5d5d5YNAiahgGLAPHgEVETS2nqAxRR29g/ZEbWkFrXnAXjOjqwjlaRHVg8AErJiYG0dHRiI2N1RTZq1cvPPfccxg8eLAu6jRoDFhEpC/VBS2/tvZ4dZQ3/Ns56Lk6IsNm0AFr3rx5WLlypeZehDKZDOXlFfMDBEHAyy+/jBUrVkiv1IAxYBGRvuUUleGzg9fx5dEbmjlaD3q3woKHvNHF1UbP1REZpsb8/K7bZYA12LhxI1asWIEuXbpgy5YtSE5OhkqlQkpKCrZu3Qpvb2+sXr0aGzdu1FW9RERUDVsLUyx4yBsH5wfhqb6ekMsE7L2cilGrD+H17eeRkV+i7xKJWhRJI1gDBgxAUlISLly4ABubqv9Dys3NRffu3dG6dWscP35cUqGGjCNYRGRorqfl48PdV/DrhRQAgI25CeYM74znBrSFqVzS/62Jmg2DHcG6cOECJkyYUG24AgClUonHHnsMf/31l5RuiIionjo4W+PTZ/2w7fn+6NpaibxiFd75+SJGrjqE/VdS9V0eUbMn+b8x9xsA45UsRET606+9I3a+NBjvP9YdjlYK/J1WgClRJxH21Wmk5BTruzyiZktSwPL19cX333+P/Pz8avfn5eXh+++/R7du3aR0Q0REEshlAib19cT++YGYMcQLcpmA3/5KwfAVBxF19IZmhXgi0h1JASssLAyJiYkYMGAAvv/+e6SnpwMA0tPT8d1332HgwIFITEzEzJkzdVIsERE1nNLcFK8/3BU/vzQYvTztkF+iwuKdFzFuzVGcT8zRd3lEzYrkZRrmzp2L1atXa04F/nuZBlEU8dJLL2H16tXSKzVgnORORMamvFzElpPxeP/Xy8grVkEmAC8EdMCc4Z1gZiLXd3lETcKg18ECgCNHjiAqKgqxsbHIzc3VLDQaGhqKIUOG6KJOg8aARUTGKjWvGO/8fAk7zyUBALq42OCjJ3rC191Wz5URNT6DD1gtHQMWERm73y6kVKyXVVAKuUzAf4I64sWgjlCYcEkHar4MdpmG+Ph45Obm1tomLy8P8fHxUrohIqJG9pCvK36fOxQPd28NdbmI/+29hnFrjiIuNU/fpREZJUkBy8vL677zqyIjI+Hl5SWlGyIiagKO1mZY80xvfPJ0L9hbmuJici7GfHwU35xMuO+SPESkTVLAEkXxvm86vimJiIzLIz3csHvuUAzp5ISiMjUWfP8nZm+NRW5xmb5LIzIajX5yPTExscaV3omIyDC1sjHHhil98cpD3jCRCdh5LgkP/+8wYhOy9V0akVEwqe8LlixZovX8wIED1bZTq9VITEzE1q1b0a9fvwYVR0RE+iOTCZgZ2AH92jtg9pazSMgswsRPY7BoTFc8278t79RBVIt6X0Uok/0z6CUIwn1PAbq5uWH79u3w9/dvWIVGgFcRElFzl1NUhoU//Ild5ytuHv24Xxu8M84X5qZcM4uMl0Et03Dw4EEAFXOrhg0bhsmTJyM0NLRKO7lcDgcHB3h7e2uFsuaIAYuIWgJRFLHu0N/44LfLKBeBnh52WPtsb7S2tdB3aUQNYlAB698WL16MoKAgDB06VJc1GR0GLCJqSQ5fS8NLW84iu7AMTtYKRD7jh75eDvoui6jeDDZgUQUGLCJqaRIyCzFj4ylcTsmDiUxAxGPd8XgfD32XRVQvBrPQ6COPPILTp083qKOioiJ8+OGH+PTTTxv0eiIiMhweDpb4YdZAPNKjNVTlIuZ/9ydW/XGVS/MQ3VWvgJWQkIC+ffviwQcfRHR09H1XcQeAU6dOYc6cOWjbti3eeustODk5NbhYIiIyHJYKE/xvUi/8J6gDAGDVH9cw/7s/UaYu13NlRPpXr1OEoigiKioKS5YsQXx8PGQyGby9vdG7d2+4uLjA3t4eRUVFyMzMxLVr13Dq1Cnk5ORAJpPhiSeewHvvvYd27do14pejHzxFSEQt3eYT8XjzxwtQl4sY0skJkc/0ho25qb7LIqqVwc3BKi8vxy+//ILo6GgcPHgQmZmZVdrIZDL06NED48aNw/Tp0+Hm5qaTghvL+PHjceDAATz44IP47rvv6vVaBiwiImD/5VT8Z/MZFJaq4e1qg41T+6KV0lzfZRHVyOAC1r0uXbqExMREZGRkwMLCAs7OzujWrRtsbW11UWOT2L9/P/Lz87FhwwYGLCKiBjqfmIMp0SeRnl8CLycrfD29H9zsuIwDGabG/Pyu90ru1fHx8YGPj48uDqU3QUFBNa5KT0REddO9jS1+mDkQT31+HDfSC/DEZ8ewZUZ/eDhY6rs0oiZV7xVAo6KikJaW1hi1NNihQ4cwZswYuLm5QRAE7Nixo0qbyMhIeHl5wdzcHH5+fjh8+HDTF0pE1AJ4Olrim7ABaOdoicSsIjzx2TH8nZav77KImlS9A9a0adPg5uaGIUOG4KOPPsK1a9cao656KSgoQM+ePfHJJ59Uu3/btm2YM2cOXn/9dZw9exZDhgzBqFGjEB8f38SVEhG1DO52Ftj2wgB0cLZCck4xnlx3HNfu5Om7LKImU++AFRMTg3nz5iEzMxPz58+Ht7c3fHx88Nprr+H48eONUeN9jRo1Cu+++y4ee+yxavevWLEC06ZNw/Tp0+Hj44NVq1bBw8OjwWtylZSUIDc3V+tBRETaXJTm2PbCAHi72iAtrwRPrjuOqwxZ1ELUO2D1798fERER+Ouvv3Dt2jV88MEHcHZ2xrJlyzBo0CC0bt0aL7zwAnbt2oWSkpLGqLleSktLcfr0aQQHB2ttDw4ORkxMTIOOGRERAVtbW83Dw4OrFxMRVcfJ2gxbn++P7u62yCwoxbNfnEB8RqG+yyJqdJLuwtyhQwfMmzcPhw4dwp07d/DFF1+gX79++PrrrzFmzBg4OTlhwoQJ+Oqrr6pdyqEppKenQ61Ww8XFRWu7i4sLUlJSNM9HjhyJxx9/HLt27UKbNm1w8uTJGo+5cOFC5OTkaB4JCQmNVj8RkbGzs1Tgq2l90cXFBql5JXh2/Qmk5hbruyyiRiUpYP2bo6MjpkyZgh07diA9PR3bt2/HE088gaNHjyI0NBQuLi4IDAzUVXf1JgiC1nNRFLW27d69G2lpaSgsLERiYiL8/f1rPJaZmRmUSqXWg4iIalYZsjwdLBGfWYiQ9f+H7MJSfZdF1Gh0FrD+zdzcHI8++ijWr1+P5ORkHDlyBHPnztUaMWoqTk5OkMvlVfpOTU2tMqpFRESNp5XSHF9P74dWNma4cicPk6NOoqBEpe+yiBpFowQsURRx7do1JCYmQhAEDBw4EMuWLcPly5cbo7taKRQK+Pn5Yc+ePVrb9+zZg4EDBzZ5PURELZmHgyU2Te8HO0tTxCZkY8bGUyhRqfVdFpHOSQpYP/74I6ZOnYqsrCzNtps3b6J79+7w9vZG27Zt8cwzz6C8vHFv/Jmfn4/Y2FjExsYCAG7cuIHY2FjNMgzh4eH44osv8OWXX+LSpUuYO3cu4uPjERYW1qh1ERFRVZ1dbBA9pS+sFHLEXM/Awh/OQwc3FSEyKJJWcl+7di0SExNhb2+v2TZnzhxcvHgRw4YNQ0ZGBrZu3Yphw4Zh2rRpkoutyalTpxAUFKR5Hh4eDgAIDQ1FdHQ0nnzySWRkZGDJkiVITk6Gr68vdu3ahbZt2zZaTUREVLMHPOzw6bN+mBJ9Ej+cuY1OrWwwM7CDvssi0hlJ9yL09PTEiBEjsH79egBATk4OnJ2dMWHCBGzZsgVlZWXo1asXlEplg5dEMAa8FyERUcNsPHYTb/34FwQBWPusH0Z2c9V3SdSCNObnt6RThGlpaWjdurXm+ZEjR6BSqfDUU08BAExNTTFixAjExcVJq5KIiJql5wa0w3MD2kIUgbnbYvFXUo6+SyLSCUkBS6lUIiMjQ/P8wIEDkMlkGDJkiGabqakpCgoKpHRDRETN2FuPdMXgjk4oLFVjxoZTSM3jGllk/CQFLG9vb+zcuROZmZnIycnB1q1b0bt3b605Wbdu3eJyCEREVCMTuQxrnu6N9s5WSMopxvMbT6O4jFcWknGTFLBmz56NpKQkuLu7w8PDA0lJSVpX5qnVahw5cgQ9e/aUXCgRETVftpamWB/qD1uLiuUb3v3lor5LIpJEUsCaMGEC1qxZg27duqFz586IiIjA1KlTNfv37t2LwsJCPPTQQ5ILJSKi5s3LyQofP9ULALDpeDx+/jNJzxURNZykqwipAq8iJCLSnWW/XUbkgeuwMTPBz7MHo62jlb5LombKYK8ijI+PR25ubq1t8vLyNAt+EhER3U/4iM7o09YeeSUqvLj5LFd6J6MkKWB5eXlh9erVtbaJjIyEl5eXlG6IiKgFMZHL8L+nesHO0hTnb+fg/V+b/jZrRFJJCliiKN739gY8A0lERPXlZmeBDydWXCAVdfQmdv+VoueKiOqnUW72/G+JiYmwsbFp7G6IiKiZGd7VBdMHV5wBmf/tOdzOLtJzRUR1V+97ES5ZskTr+YEDB6ptp1arkZiYiK1bt6Jfv34NKo6IiFq2BQ954+TNTJxLzMHCH85jwxR/CIKg77KI7qveVxHKZP8MegmCcN9TgG5ubti+fTv8/f0bVqER4FWERESNJy41H6P/dxilqnIsn9gDj/fx0HdJ1Ew05ud3vUew9u/fD6BibtWwYcMwefJkhIaGVmknl8vh4OAAb29vrVBGRERUHx1bWWPu8M744LfLeOfniwjo7IxWSnN9l0VUK0nrYC1evBhBQUEYOnSoLmsyOhzBIiJqXCp1OcZHxuD87RyM6OqCdSF+PFVIkjXm5zcXGtUBBiwiosZ3KTkXYz4+AlW5iI+f6oUxPd30XRIZOYM6RXiv0tJS7NixAydPnkR2djbU6qoLwgmCgPXr10vtioiIWjCf1kr8J6gjVu+9hkU//YWBHRzhaG2m77KIqiVpBOvWrVsYMWIErl+/Xutkd0EQqg1ezQVHsIiImkapqhxjPj6CK3fy8GhPN/zv7r0LiRrCYEew5s6di7i4OISEhGDq1Klo06YNTEwkD4oRERFVS2Eiw/LHe2DcmqP46VwSxvdyR5B3K32XRVSFpBEsOzs7+Pv7Y8+ePbqsyehwBIuIqGm998tFfH74Bto7WeG3OUOhMOHV6lR/Bnuz5/LycvTqxeFZIiJqWrMf7AQnawX+Ti/AxmM39V0OURWSAtaAAQNw6dIlXdVCRERUJzbmppgX3AUAsHrvNWTkl+i5IiJtkgLW+++/j/379+O7777TVT1ERER18ngfD3RtrUResQor/7iq73KItEiakb5z504EBQXhySefREBAAHr16gVbW9sq7QRBwJtvvimlKyIiIi1ymYC3xnTFpHXHsflEPJ7t3xberpwHS4ZB0iT3ut4Ch8s0EBFRY5m56TR+vZCCQR0dsWlaP67wTnVmsMs0VN6XkIiISF9eG+2DvZdTcTQuA3su3kFwN1d9l0QkLWAFBAToqg4iIqIG8XCwxPTBXog8cB3v7bqEgC7OMDOR67ssauG4cAgRERm9WUEd4WxjhlsZhdh8Il7f5RBJD1gqlQorV65E3759oVQqtVZyj42NxaxZs3D1Kq/uICKixmNtZoI5wzsBANYevI7isuY775eMg6SAVVRUhKCgIMybNw+3bt2CUqnUuiehl5cXoqKisHHjRsmFEhER1eZxPw+42ZrjTm4Jtp1M0Hc51MJJClhLly7F0aNHERERgZSUFEyfPl1rv62tLQICArB7925JRRIREd2PwkSGWUEdAQCRB+I4ikV6JSlgbdu2DYGBgViwYAEEQaj20tj27dsjPp7nw4mIqPE93qcNWt8dxfrmFEexSH8kBaz4+Hj4+/vX2kapVCInJ0dKN0RERHViZiL/ZxRrP+dikf5IClg2NjZIS0urtc3169fh7OwspRsiIqI6e+LuKFZKbjFHsUhvJAWs/v37Y+fOnTWOUCUmJmLXrl0YOnSolG6IiIjqzMxEjlmBHQBUjGKVqDiKRU1PUsCaP38+MjMzMXz4cMTExEClUgEACgsLsXfvXgQHB6OsrAzh4eE6KZaIiKgunvD3gKvy7igWrygkPZB0L0IAWLt2LWbPnl3tvQblcjkiIyOrXF3Y3PBehEREhmfjsZt468e/0NrWHAfmB3J1d6qiMT+/JS80GhYWhnPnzuHFF1+Ev78/OnTogF69eiEsLAxnz55t9uGKiIgM0xN9KkaxknOK8cOZ2/ouh1oYySNYxBEsIiJD9cXhv/HuL5fQxcUGv80ZUu1yQtRyGewI1ogRI/DVV1+hoKBAV/UQERHpzON9PGCpkOPKnTwcu56h73KoBZEUsPbt24fJkyfDxcUFzz77LHbv3o3y8nJd1UZERCSJrYUpJvRuAwCIirmp32KoRZG80OjSpUvRvn17bN68GaNHj4a7uzv++9//4syZM7qqkYiIqMFCB7YDAPxx6Q7iMwr1Wwy1GJIClru7O1555RX8+eefiI2Nxdy5cyGXy7Fy5Ur4+/ujW7du+OCDD5CQwEtkiYhIPzq2ssbQzs4QRWDDsZv6LodaCJ1PchdFEXv37sWmTZuwfft25OfnQyaToaysTJfdGBROciciMmz7L6diSvRJ2JiZ4PhrD8LKzETfJZEBMNhJ7tURBAHDhw/H4sWL8dJLL8HExITzsoiISK8COjvDy8kKeSUqfH8mUd/lUAug04CVnZ2NdevWYejQoejQoQMiIiJgYWGByZMn67IbIiKiepHJBIQOaAsAiI65ifJyrlBEjUvyGGlpaSl27tyJTZs24ddff0VZWRnkcjkeeughhISEYOzYsTA3N9dFrURERA02sY8HPvz9Kv5OK8Cha2kI7NJK3yVRMyYpYM2YMQPff/89cnJyIIoi+vTpg5CQEEyaNAnOzs66qpGIiEgyazMTPN6nDaKO3kR0zE0GLGpUkgLW+vXr0a5dO/znP/9BSEgIOnfurKu6iIiIdC50QDtEx9zEgStpuJ6Wjw7O1vouiZopSQHr4MGDGDJkiK5qISIialTtnKwQ1KUV9l1OxTcnE7BwtI++S6JmStIkd4YrIiIyNk/6ewAAfjh7Gyo1r3KnxiH5KkKVSoWVK1eib9++UCqVMDH5Z1AsNjYWs2bNwtWrV6V2Q0REpBNBXVrBwUqBtLwSHLqWpu9yqJmSFLCKiooQFBSEefPm4datW1Aqlfj3uqVeXl6IiorCxo0bJRdKRESkCwoTGcY94A4A+O4018SixiEpYC1duhRHjx5FREQEUlJSMH36dK39tra2CAgIwO7duyUVSUREpEsT/SpuAP3HxVRkFZTquRpqjiQFrG3btiEwMBALFiyAIAgQBKFKm/bt2yM+Pl5KN01i/PjxsLe3x8SJE/VdChERNbKubkp0ba1EqbocO/9M0nc51AxJCljx8fHw9/evtY1SqUROTo6UbprE7NmzeSqTiKgFqRzF4mlCagySApaNjQ3S0mqfIHj9+nWjWHQ0KCgINjY2+i6DiIiayNgH3GAiE/BnYg6upOTpuxxqZiQFrP79+2Pnzp01jlAlJiZi165dGDp0qJRucOjQIYwZMwZubm4QBAE7duyo0iYyMhJeXl4wNzeHn58fDh8+LKlPIiJq3hytzTDMu2I1d94AmnRNUsCaP38+MjMzMXz4cMTExEClUgEACgsLsXfvXgQHB6OsrAzh4eGSiiwoKEDPnj3xySefVLt/27ZtmDNnDl5//XWcPXsWQ4YMwahRo7Tmfvn5+cHX17fKIymJ596JiFqqytOEP5zhmlikW5JWch86dCjWrFmD2bNnay06WnmqTS6XIzIyEn5+fpKKHDVqFEaNGlXj/hUrVmDatGmaqxhXrVqF3bt349NPP0VERAQA4PTp05Jq+LeSkhKUlJRonufm5urs2ERE1HSCvFvB0UqB9PwSHLyahgd9XPRdEjUTkhcaDQsLw7lz5/Diiy/C398fHTp0QK9evRAWFoazZ89WWbpB10pLS3H69GkEBwdrbQ8ODkZMTEyj9BkREQFbW1vNw8PDo1H6ISKixmUql2FcL66JRbonaQSrko+PD1avXq2LQ9Vbeno61Go1XFy0/9fh4uKClJSUOh9n5MiROHPmDAoKCtCmTRts3769xiskFy5cqHXaMzc3lyGLiMhITfRrg/VHbuCPS3eQVVAKeyuFvkuiZkAnAcsQ3LsGlyiK1a7LVZP6LIZqZmYGMzOzOrcnIiLD5dNaiW5uSvyVlIuf/0xCyIB2+i6JmgHJpwj1zcnJCXK5vMpoVWpqapVRLSIiouqMfcANALDrfN3PfBDVpl4B6+mnn8YPP/zQ4M6kvr46CoUCfn5+2LNnj9b2PXv2YODAgTrti4iImqdRvq0BACduZCA9v+Q+rYnur14Ba+vWrbhw4UKDO2vo6/Pz8xEbG4vY2FgAwI0bNxAbG6tZhiE8PBxffPEFvvzyS1y6dAlz585FfHw8wsLCGlwrERG1HB4OlujRxhblIrD7L45ikXT1noMVGxvb5LeUOXXqFIKCgjTPKyeYh4aGIjo6Gk8++SQyMjKwZMkSJCcnw9fXF7t27ULbtm2btE4iIjJeo7u3xp+JOdh1PhnP9OPnB0kjiKIo1rWxTFYx4FWfyeOVKrt5++238dZbb9X79YYsNzcXtra2yMnJgVKp1Hc5RETUAPEZhRi6fD9kAnDy9eFwtObFTM1dY35+12sEKyoqSnKHDzzwgORjEBER6ZqnoyV83ZW4cDsXv1+8g6f6euq7JDJi9QpYoaGhjVUHERGR3o3u3hoXbudi1/lkBiySxOiXaSAiItKV0XevJoy5noHMglI9V0PGjAGLiIjornZOVujaWgl1uYg9F3k1ITUcAxYREdG/PNyjYhTrFy46ShIwYBEREf3LKF9XAEBMXDqyC3makBqGAYuIiOhf2jtbw9vVBqpyEb9fvKPvcshIMWARERHd4+HuFacJd51P1nMlZKwYsIiIiO4x6m7AOhqXjpzCMj1XQ8ZIcsBSqVRYuXIl+vbtC6VSCROTf5bWio2NxaxZs3D16lWp3RARETWZjq2s0cXFBmVqEX9c4mlCqj9JAauoqAhBQUGYN28ebt26BaVSiX/fecfLywtRUVFNfu9CIiIiqYK7uQAA9l1O1XMlZIwkBaylS5fi6NGjiIiIQEpKCqZPn66139bWFgEBAdi9e7ekIomIiJpakHcrAMChq2koU5fruRoyNpIC1rZt2xAYGIgFCxZAEIRqbwLdvn17xMfHS+mGiIioyfVsYwcHKwXySlQ4fStL3+WQkZEUsOLj4+Hv719rG6VSiZycHCndEBERNTm5TEBgZ2cAwH6eJqR6khSwbGxskJaWVmub69evw9nZWUo3REREelF5mpDzsKi+JAWs/v37Y+fOnTWOUCUmJmLXrl0YOnSolG6IiIj0YmhnZ8hlAq6l5iMhs1Df5ZARkRSw5s+fj8zMTAwfPhwxMTFQqVQAgMLCQuzduxfBwcEoKytDeHi4ToolIiJqSrYWpvBraw8A2H+Fo1hUdyb3b1KzoUOHYs2aNZg9ezaGDBmi2W5jYwMAkMvliIyMhJ+fn7QqiYiI9GSYdyv8341M7LuciucGtNN3OWQkJAUsAAgLC0NAQADWrl2LEydOIDMzE0qlEv369cOsWbPQrVs3XdRJRESkF8O8W+H9Xy/j2PUMFJWqYaGQ67skMgKSAxYA+Pj4YPXq1bo4FBERkUHp1Moa7nYWuJ1dhJjr6XjQx0XfJZERkLxMQ25ubq1t8vLyuA4WEREZLUEQMIxXE1I9SQpYXl5e9x25ioyMhJeXl5RuiIiI9KoyYO2/nKp1SziimkgKWKIo3vcXjb+IRERk7Pq3d4SZiQxJOcW4cidP3+WQEZAUsOoiMTFRc1UhERGRMbJQyDGwgyMAniakuqn3JPclS5ZoPT9w4EC17dRqNRITE7F161b069evQcUREREZimHerbD/ShoOXE7DrMCO+i6HDJwg1vMcnkz2z6CXIAj3PQXo5uaG7du33/eehcYsNzcXtra2yMnJgVKp1Hc5RETUCBKzCjH4g/2QywSceWMEbC1N9V0SSdSYn9/1HsHav38/gIq5VcOGDcPkyZMRGhpapZ1cLoeDgwO8vb21QhkREZExamNvic4u1rh6Jx8Hr6Xh0Z5u+i6JDFi9A1ZAQIDm74sWLUJQUBDvNUhERC1CQGdnXL2Tj5i4dAYsqpWkoaVFixYxXBERUYsxsIMTAODY3xl6roQMnU5WcgcqJrWnp6ejpKSk2v2enp666oqIiEgv/L0cIJcJuJVRiNvZRXC3s9B3SWSgJAes06dP47XXXsOhQ4dQWlpabRtBEKBSqaR2RUREpFfWZibo0cYWZ+Ozcex6Bib6tdF3SWSgJJ0ijI2NxZAhQ3Ds2DEEBwdDFEX06NEDwcHBcHJygiiKCAgIQEhIiK7qJSIi0qsB7SvWw4q5nq7nSsiQSQpY77zzDgDgxIkT+PHHHwEA48ePx6+//oqbN28iLCwMFy5cwKJFi6RXSkREZAAq52Edv57Bu5VQjSQFrCNHjuDRRx+Fj4+PZlvlL5uFhQU++eQTuLm54bXXXpNWJRERkYHwa2sPU7mApJxi3Moo1Hc5ZKAkBaycnBy0b99e89zU1BT5+fn/HFwmQ2BgIPbu3SulGyIiIoNhoZCjl6c9AF5NSDWTFLBatWqFrKwszXNXV1dcu3ZNq01xcTEKC5nwiYio+fhnHhYDFlVPUsDq2rUrrly5onk+aNAg/P777zh+/DgA4NKlS/jmm2/g7e0trUoiIiIDUnnj52Och0U1kBSwHn74YRw6dAjJyckAgFdeeQWiKGLQoEFwdnZG9+7dkZ2dzTlYRETUrDzgaQczExnS80sQl5p//xdQiyMpYIWFheH27dtwdKxI8j179sTevXvx0EMPwcnJCcOHD8fOnTsxfvx4nRRLRERkCMxM5PBv5wCA87CoepIWGjU1NYWLi4vWtoEDB+KXX36RVBQREZGhG9DBEUfi0hETl4HnBrTTdzlkYCSNYMnlcjzzzDO6qoWIiMhoDLg7D+v4jQyUl3MeFmmTFLCUSiU8PDx0VQsREZHR6O5uCyuFHNmFZbiUkqvvcsjASApYffv2xblz53RVCxERkdEwlcvg73V3HhaXa6B7SApYixcvxr59+7BhwwZd1UNERGQ0/r1cA9G/SZrk/vvvvyMwMBBTp07Fxx9/jL59+8LFxQWCIGi1EwQBb775pqRCiYiIDM2A9hX3Jfy/G5lQqcthIpc0bkHNiCBKWCFNJqvbL5IgCFCr1Q3txuDl5ubC1tYWOTk5UCqV+i6HiIiaiLpcRK8lvyO3WIUd/xmEBzzs9F0S1UNjfn5LGsHav3+/ruogIiIyOnKZgH7tHbHn4h2c+DuDAYs0JAWsgIAAXdVBRERklPza2mPPxTuITcjWdylkQHiymIiISILKUauz8dl6rYMMCwMWERGRBD3a2EImACm5xUjJKdZ3OWQgGLCIiIgksFSYoItrxQTp2IQsPVdDhoIBi4iISCKeJqR7MWARERFJ1MvTDgBwlhPd6S5JASs+Ph65ubXffykvLw/x8fFSuiEiIjJove6OYJ1PzIFKXa7fYsggSApYXl5eWL16da1tIiMj4eXlJaWbRpeXlwd/f3888MAD6N69Oz7//HN9l0REREakg7M1bMxMUFSmxpU7efouhwyApHWwRFHE/RaCl7BQfJOxtLTEwYMHYWlpicLCQvj6+uKxxx6Do6OjvksjIiIjIJMJ6OlhhyNx6YhNyEY3N1t9l0R61uhzsBITE2FjY9PY3Ugil8thaWkJACguLoZarTaKYEhERIaDE93p3+o9grVkyRKt5wcOHKi2nVqtRmJiIrZu3Yp+/fo1qLhKhw4dwvLly3H69GkkJydj+/btGDdunFabyMhILF++HMnJyejWrRtWrVqFIUOG1LmP7OxsBAQE4Nq1a1i+fDmcnJwk1UxERC1L5UR3ruhOQAMC1ttvv635uyAIOHDgQI0hCwDc3NzwwQcfNKQ2jYKCAvTs2RNTpkzBhAkTquzftm0b5syZg8jISAwaNAifffYZRo0ahYsXL8LT0xMA4Ofnh5KSkiqv/f333+Hm5gY7OzucO3cOd+7cwWOPPYaJEyfCxcVFUt1ERNRyVI5gxaXmI6eoDLYWpvotiPRKEOt5LuzgwYMAKuZWDRs2DJMnT0ZoaGiVdnK5HA4ODvD29oZMprszkYIgVBnB6tevH3r37o1PP/1Us83Hxwfjxo1DREREvfuYOXMmhg0bhscff7za/SUlJVphLTc3Fx4eHo1yN24iIjIeQ5btQ0JmEb6a1hdDOjnruxy6j9zcXNja2jbK53e9R7D+fYPnRYsWISgoCEOHDtVpUfVRWlqK06dP49VXX9XaHhwcjJiYmDod486dO7CwsIBSqURubi4OHTqEmTNn1tg+IiICixcvllQ3ERE1P7087JGQWYTY+GwGrBZO0tDSokWL9BquACA9PR1qtbrK6TwXFxekpKTU6RiJiYkYOnQoevbsicGDB+PFF19Ejx49amy/cOFC5OTkaB4JCQmSvgYiImoeKk8Tch4WSVqmwZAIgqD1XBTFKttq4ufnh9jY2Dr3ZWZmBjMzs/qUR0RELcAD/1rRvT6fQ9T8SA5YpaWl2LFjB06ePIns7Gyo1eoqbQRBwPr166V2VS0nJyfI5fIqo1WpqamcpE5ERE2qm5sSCrkMmQWlSMgsgqejpb5LIj2RFLBu3bqFESNG4Pr167WuG9WYAUuhUMDPzw979uzB+PHjNdv37NmDsWPHNkqfRERE1TEzkcPHTYlzCdk4m5DFgNWCSQpYc+fORVxcHEJCQjB16lS0adMGJia6P+uYn5+PuLg4zfMbN24gNjYWDg4O8PT0RHh4OEJCQtCnTx8MGDAA69atQ3x8PMLCwnReCxERUW16edhVBKz4bIx9wF3f5ZCeSEpD+/btw4MPPogNGzboqp5qnTp1CkFBQZrn4eHhAIDQ0FBER0fjySefREZGBpYsWYLk5GT4+vpi165daNu2baPWRUREdK9ennaIjuFE95ZOUsAqLy9Hr169dFVLjQIDA+9765pZs2Zh1qxZjV4LERFRbSqvJLyYlIsSlRpmJnL9FkR6IWmZhgEDBuDSpUu6qoWIiMjoeTpYwsFKgVJ1OS4m5eq7HNITSQHr/fffx/79+/Hdd9/pqh4iIiKjJgiCZhTrHE8TtliSThHu3LkTQUFBePLJJxEQEIBevXrB1ta2SjtBEPDmm29K6YqIiMhodG2txL7LqbickqfvUkhP6n0vwn+r6z0GBUGodn2s5qIx72VERETGZ+e5JLy05Sx6edph+6xB+i6HamBQ9yL8t/379+uqDiIiomaji6sNAOBqSh7Ky0XIZFzRvaWRFLD+feNnIiIiquDlZAVTuYCCUjVuZxfBw4ELjrY0kia5ExERUVWmchk6OFsDAK5wHlaLJDlgqVQqrFy5En379oVSqdRayT02NhazZs3C1atXpXZDRERkVCpPE165w4DVEkk6RVhUVITg4GDExMTAyckJSqUSBQUFmv1eXl6IioqCg4MD3n33XcnFEhERGQtNwOIIVoskaQRr6dKlOHr0KCIiIpCSkoLp06dr7be1tUVAQAB2794tqUgiIiJj482A1aJJCljbtm1DYGAgFixYAEEQIAhVr5Jo37494uPjpXRDRERkdDq7VASs62n5KFWV67kaamqSAlZ8fDz8/f1rbaNUKpGTkyOlGyIiIqPjbmcBazMTqMpF3EgvuP8LqFmRFLBsbGyQlpZWa5vr16/D2dlZSjdERERGRxAEdHapuJLwcgrvSdjSSApY/fv3x86dO2scoUpMTMSuXbswdOhQKd0QEREZpS6uFauDcx5WyyMpYM2fPx+ZmZkYPnw4YmJioFKpAACFhYXYu3cvgoODUVZWhvDwcJ0US0REZEwqJ7pf5VINLY6kZRqGDh2KNWvWYPbs2RgyZIhmu41NxS+UXC5HZGQk/Pz8pFVJRERkhConuvOmzy2PpIAFAGFhYQgICMDatWtx4sQJZGZmQqlUol+/fpg1axa6deumizqJiIiMTuUIVmJWEfJLVLA2k/yxS0ZCJz9pHx8frF69WheHIiIiajbsrRRoZWOG1LwSXL2Th96e9vouiZqIpDlYI0aMwFdffaW1ejsRERH9gyu6t0ySAta+ffswefJkuLi44Nlnn8Xu3btRXs7F1IiIiCp1cWHAaokkLzS6dOlStG/fHps3b8bo0aPh7u6O//73vzhz5oyuaiQiIjJaHMFqmSQFLHd3d7zyyiv4888/ERsbi7lz50Iul2PlypXw9/dHt27d8MEHHyAhIUFX9RIRERkVTcC6kwdRFPVcDTUVQdTxT1sURezduxebNm3C9u3bkZ+fD5lMhrKyMl12Y1Byc3Nha2uLnJwcKJVKfZdDREQGpKhUja6LfoMoAidfHw5nGzN9l0R3Nebnt86vFxUEAcOHD0enTp3Qpk0bLF++XLMAKRERUUtjoZCjrYMlbmYU4kpKHgNWC6HTgJWdnY1vvvkGmzZtQkxMDERRhI2NDSZMmKDLboiIiIxKF1cb3MwoxOWUXAzu5KTvcqgJSA5YpaWl2LlzJzZt2oRff/0VZWVlkMvleOihhxASEoKxY8fC3NxcF7USEREZpS6uSuz+6w5vmdOCSApYM2bMwPfff4+cnByIoog+ffogJCQEkyZNgrOzs65qJCIiMmpcqqHlkRSw1q9fj3bt2uE///kPQkJC0LlzZ13VRURE1Gx00dz0OR/l5SJkMkHPFVFjkxSwDh48qHWTZyIiIqqqnaMlFCYyFJWpkZBViLaOVvouiRqZpHWwGK6IiIjuz0QuQ0dnawDAZZ4mbBHqHbB69+6NdevWaW3bvXs3wsPDq22/ePFimJjw7uFERNSyVZ4mjEvN13Ml1BTqHbBiY2ORkpKite348eNYvXp1ja/hyrVERNTSeTpYAgASMgv1XAk1BUmnCImIiKhuPCoDVhYDVkvAgEVERNQEPOwtAAAJmUV6roSaAgMWERFRE6gcwUrKLoK6nFNnmjsGLCIioibgojSHqVyAqlxEcg5HsZo7BiwiIqImIJcJcLfjacKWokHrJ2zatAnHjx/XPI+LiwMAjB49ukrbyn1EREQtnYeDJW5mFCIhsxADOjjquxxqRA0KWHFxcdUGp99++63a9oLAWwIQERHxSsKWo94B68aNG41RBxERUbPnYc+1sFqKegestm3bNkYdREREzZ6Hw905WFmcg9XccZI7ERFRE+EIVsvBgEVERNREKudgpeaVoLhMredqqDExYBERETURe0tTWCnkAIBEniZs1hiwiIiImoggCLySsIVgwCIiImpCmoDFeVjNGgMWERFRE+JE95aBAYuIiKgJaZZq4O1ymjUGLCIioiakGcHiHKxmjQGLiIioCXEOVsvAgEVERNSE2thXnCLMLVYhp6hMz9VQY2HAIiIiakJWZiZwtFIA4ChWc8aARURE1MR4mrD5Y8C6y8TEBA888AAeeOABTJ8+Xd/lEBFRM8bFRps/E30XYCjs7OwQGxur7zKIiKgF8LDnUg3NHUewiIiImhhHsJo/owhYhw4dwpgxY+Dm5gZBELBjx44qbSIjI+Hl5QVzc3P4+fnh8OHD9eojNzcXfn5+GDx4MA4ePKijyomIiKriau7Nn1GcIiwoKEDPnj0xZcoUTJgwocr+bdu2Yc6cOYiMjMSgQYPw2WefYdSoUbh48SI8PT0BAH5+figpKany2t9//x1ubm64efMm3NzccOHCBTz88MM4f/48lEplo39tRETU8lSu5p6YVQRRFCEIgp4rIl0TRFEU9V1EfQiCgO3bt2PcuHGabf369UPv3r3x6aefarb5+Phg3LhxiIiIqHcfo0aNwjvvvIM+ffpUu7+kpEQrrOXk5MDT0xMJCQkMZUREdF9l6nL4vbMH5SKw/78BcFaa67ukFik3NxceHh7Izs6Gra2tTo9tFCNYtSktLcXp06fx6quvam0PDg5GTExMnY6RlZUFS0tLmJmZITExERcvXkT79u1rbB8REYHFixdX2e7h4VG/4omIqMXruErfFVBGRgYD1r3S09OhVqvh4uKitd3FxQUpKSl1OsalS5fwwgsvQCaTQRAErF69Gg4ODjW2X7hwIcLDwzXPs7Oz0bZtW8THx+v8B9Tc+Pv74+TJk/ouo870WW9j963r40s9npTX1/e19Wlf+T9cjlDfH9/fhtM33991a195Bqq2z/yGMvqAVene89f1Oac9cOBAnD9/vs59mZmZwczMrMp2W1tb/gN8H3K53Ki+R/qst7H71vXxpR5Pyuvr+9qG9KVUKo3qd1cf+P42nL75/q5fXzKZ7q/5M4qrCGvj5OQEuVxeZbQqNTW1yqgW6d9//vMffZdQL/qst7H71vXxpR5Pyuvr+1pj+z00Fsb2feX7u+mO1xLf381mkrufnx8iIyM127p27YqxY8c2aJJ7feXm5sLW1hY5OTlG9b83Iro/vr+Jmq/GfH8bxSnC/Px8xMXFaZ7fuHEDsbGxcHBwgKenJ8LDwxESEoI+ffpgwIABWLduHeLj4xEWFtYk9ZmZmWHRokXVnjYkIuPG9zdR89WY72+jGME6cOAAgoKCqmwPDQ1FdHQ0gIqFRpctW4bk5GT4+vpi5cqVGDp0aBNXSkRERGQkAYuIiIjImBj9JHciIiIiQ8OARURERKRjDFhEREREOsaA1cTGjx8Pe3t7TJw4Ud+lEJFEP//8M7p06YJOnTrhiy++0Hc5RKRDUj+vOcm9ie3fvx/5+fnYsGEDvvvuO32XQ0QNpFKp0LVrV+zfvx9KpRK9e/fGiRMnGuWWG0TU9KR+XnMEq4kFBQXBxsZG32UQkUT/93//h27dusHd3R02NjYYPXo0du/ere+yiEhHpH5eM2D9y6FDhzBmzBi4ublBEATs2LGjSpvIyEh4eXnB3Nwcfn5+OHz4cNMXSkSSSX2/JyUlwd3dXfO8TZs2uH37dlOUTkT3YQif5wxY/1JQUICePXvik08+qXb/tm3bMGfOHLz++us4e/YshgwZglGjRiE+Pl7Txs/PD76+vlUeSUlJTfVlEFEdSH2/Vze7oq43mCeixqWLz3PJRKoWAHH79u1a2/r27SuGhYVpbfP29hZfffXVeh17//794oQJE6SWSEQ60pD3+9GjR8Vx48Zp9s2ePVv8+uuvG71WIqofKZ/nUj6vOYJVR6WlpTh9+jSCg4O1tgcHByMmJkZPVRFRY6jL+71v3764cOECbt++jby8POzatQsjR47UR7lEVA9N9XluFDd7NgTp6elQq9VwcXHR2u7i4oKUlJQ6H2fkyJE4c+YMCgoK0KZNG2zfvh3+/v66LpeIJKjL+93ExAQfffQRgoKCUF5ejgULFsDR0VEf5RJRPdT181zq5zUDVj3dO8dCFMV6zbvgVUZExuN+7/dHH30Ujz76aFOXRUQ6cL/3t9TPa54irCMnJyfI5fIqo1WpqalVUjARGTe+34mar6Z6fzNg1ZFCoYCfnx/27NmjtX3Pnj0YOHCgnqoiosbA9ztR89VU72+eIvyX/Px8xMXFaZ7fuHEDsbGxcHBwgKenJ8LDwxESEoI+ffpgwIABWLduHeLj4xEWFqbHqomoIfh+J2q+DOL93aBrD5up/fv3iwCqPEJDQzVt1qxZI7Zt21ZUKBRi7969xYMHD+qvYCJqML7fiZovQ3h/816ERERERDrGOVhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhE1CJFR0dDEATNY9KkSVr7Dxw4AEEQ8Pbbb+unwHvExcVp1duuXTt9l0REtWDAIiKD0qNHD60gUd0jLi5OZ/2NHTsWixYtwsSJE3V2zNo89dRTEAQBW7durbVdRkYGzMzM4OTkhNLSUjg4OGDRokVYtGgRbG1tm6RWImo4E30XQERUqbi4GJcuXUKrVq0wc+bMatvI5XJ06NBBZ32OGzcOkydP1tnx7mfatGnYunUroqKiqoya/dumTZtQWlqKkJAQKBQKODg4aEbToqOjm6ZYImowBiwiMhjnzp2DSqXC0KFDDebUnK49+OCDaNeuHf744w8kJCTAw8Oj2nZRUVEAKgIZERkfniIkIoNx5swZAECfPn30XEnNsrOzMWTIEMjlcqxbt05r36FDhzBmzBg4OTnBzMwMnTp1whtvvIHCwkJNG0EQMGXKFJSXl2PDhg3V9nH69GmcO3cOffv2ha+vb6N+PUTUOBiwiMhgnD59GoDhBqykpCQMGTIEJ0+exLfffovnn39es2/t2rUIDAxETEwMHnnkEcyePRvu7u547733MGLECJSWlmraTpkyBTKZDNHR0RBFsUo/HL0iMn48RUhEBqNyBOuXX37B4cOHq+z38PDQW+i4evUqgoODkZWVhd9++w2BgYGafRcvXsRLL72EBx54AH/88QccHBw0+95//30sXLgQH3/8Mf773/8CqPg6RowYgd27d+PQoUMICAjQtC8pKcHmzZthaWlZ6xwtIjJsDFhEZBBKS0tx4cIFAMDKlSurbTNlyhS9BKyTJ09i9OjRkMvlOHjwIB544AGt/Z999hlUKhX+97//aYUrAFiwYAFWrFiBLVu2aAIWAEydOhW7d+/Gl19+qRWwtm/fjqysLISGhkKpVDbq10VEjYcBi4gMwp9//omysjI8+eST913CoJKfn59m1MvMzAxdu3bF8uXL8eCDD+qsrsOHD+Ojjz6Ci4sLdu/eXe0VjMePHwcA/Pbbb/jjjz+q7Dc1NcXly5e1to0bNw6Ojo747rvv8Mknn8DGxgYA8OWXXwKoCGBEZLwYsIjIIFQGpe7du9epfVlZGS5cuIAvv/wSo0aNQl5eHt544w1MmDABqampUCgUOqnr7NmzyM/Px6hRo2pc3DMzMxMA8N5779X5uAqFAs8++yxWr16Nb775BtOmTUNCQgL27t2LTp06YejQoboon4j0hJPcicggVE5w79GjR53a//XXXygtLcXIkSPh6uqKTp064emnn0ZOTg6SkpJ0VteLL76IqVOn4ttvv0VISAjUanWVNpWn8nJzcyGKYo2Pe1We7qwctYqOjkZ5eTlHr4iaAQYsIjIIlSNYdQ1YZ86cQbt27eDm5gYAuHbtGt5//324uLigTZs2OqtLJpPhiy++wPTp07Fly5ZqQ1a/fv0A/HOqsK66d+8Of39/xMTE4PLly4iOjoZcLkdoaKjO6ici/WDAIiK9Kysrw/nz52Fra4u2bdvW6TVnzpxBfHw8rK2tYW5ujs6dO+PKlSv46quvYGKi29kPgiBg3bp1mDFjBrZs2YJnnnlGK2TNmjULJiYmeOmll5CQkFDl9dnZ2Th79my1x64cxZo+fTr+/vtvjB49Gq1bt9Zp/UTU9DgHi4j07q+//kJJSQns7e1rXMHdwsICr7zyiub52bNnMX/+fEyfPh1paWmYNWsWZs6ciREjRjRKjYIg4LPPPtOELVEU8fXXX8PExAS+vr6IjIzEzJkz0aVLF4wePRodOnRAbm4u/v77bxw8eBCTJ0/G2rVrqxz3qaeeQnh4OI4ePQqAa18RNRcMWESkd5Xzr1JSUrB48eJq2wwYMEATsMrLy3Hu3Dm8+uqr6NixIzp27IiIiAhMmjQJzz33HMzNzRulTkEQsHbtWshkMqxduxaiKGLz5s0wMTHBjBkz8MADD2DFihU4dOgQfvrpJ9ja2sLT0xNz586t8bSfUqnExIkTsXHjRri4uODhhx9ulNqJqGkJYnUzL4mIDNjly5fh4+OD69evo3379gAqTjPa29tj48aNeOyxx+57jOjoaEyZMgVRUVFNerNnXai8mvHmzZt6rYOIasY5WERkdM6cOQNLS0utZRNMTU0xdOhQ/Pjjj/U61pQpUyAIgsGvmh4XFwdBECAIAm7duqXvcojoPniKkIiMztmzZ+Ht7Q2ZTPv/iMOHD8fSpUuhVqshl8trPcYDDzyARYsWaZ4b+k2VHRwctOq1s7PTXzFEdF88RUhERESkYzxFSERERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOsaARURERKRjDFhEREREOvb/VuvznvMXQYoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(E_Rs*1e6, abs(sm_spec), label=\"SNuDD Standard Model\") # abs because of the residual imaginary part from the density matrix. All elements have zero imaginary part anyway.\n",
    "plt.legend()\n",
    "plt.title(r\"Theoretical CE$\\nu$NS rate in Xe\")\n",
    "plt.ylabel(r'${\\rm d}R/{\\rm d}E_R\\,\\,\\left[{\\rm ton}^{-1}\\,{\\rm yr}^{-1}\\,{\\rm keV}^{-1}\\right]$', size=14)\n",
    "plt.xlabel(r'$E_{R}\\,\\,\\left[{\\rm keV}\\right]$', size=14)\n",
    "\n",
    "plt.ylabel(r\"Event rate [events/(ton$\\cdot$year$\\cdot$keV)]\")\n",
    "\n",
    "plt.xlim(xmin=1e-1, xmax=1e1)\n",
    "plt.ylim(ymax=1e6, ymin=1e-5)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "41c813f9-a5d8-4da2-ad1e-9277757df7d7",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "\"\"\" IN SI UNIT \"\"\"\n",
    "\"\"\" Reduced Planck Constant in J·s \"\"\"\n",
    "HBAR = 1.054571800E-34\n",
    "\"\"\" Speed of Light in m/s \"\"\"\n",
    "C = 2.99792458E8\n",
    "\"\"\" Avogaro's Number in 1/mol \"\"\"\n",
    "A_V = 6.022140857E23\n",
    "\"\"\" Charge of Electron in C \"\"\"\n",
    "CH_E = 1.6021766208E-19\n",
    "\n",
    "\"\"\" keV in J \"\"\"\n",
    "KEV = 1e3 * CH_E\n",
    "\"\"\" GeV in J \"\"\"\n",
    "GEV = 1e9 * CH_E\n",
    "\n",
    "GeV = 1.0\n",
    "MeV = GeV / 1e3\n",
    "keV = GeV / 1e6\n",
    "eV = GeV / 1e9\n",
    "\n",
    "# GeV ^ -1\n",
    "m = 1.0 / (HBAR * C / GEV) / GeV\n",
    "cm = m / 1e2\n",
    "km = m * 1e3\n",
    "mm = m / 1e3\n",
    "um = m / 1e6\n",
    "nm = m / 1e9\n",
    "pm = m / 1e12\n",
    "fm = m / 1e15\n",
    "\n",
    "# GeV\n",
    "kg = 1.0 * C**2 / GEV * GeV\n",
    "ton = kg * 1e3\n",
    "g = kg * 1e-3\n",
    "\n",
    "# GeV ^ -1\n",
    "sec = 1.0 * C * m\n",
    "dy = 24.0 * 3600.0 * sec\n",
    "yr = 365.0 * dy\n",
    "\n",
    "# Fermi's const\n",
    "Gf = 1.1663787e-5 / GeV**2\n",
    "\n",
    "# Neutron mass\n",
    "mn = 0.9314940954 * GeV\n",
    "\n",
    "# Mixing angle\n",
    "sin2the = 0.2223"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "e3f0fd69-1d5a-4e29-aeb3-126b330ea26a",
   "metadata": {},
   "outputs": [],
   "source": [
    "NSI_matrix = np.array([[0.0, 0.0, 0.0],\n",
    "                       [0.0, 0.0, 0.0],\n",
    "                       [0.0, 0.0, 0.0]])\n",
    "\n",
    "NSI_eta = 0\n",
    "NSI_phi = 0\n",
    "\n",
    "NSI_model = GeneralNSI(NSI_matrix, NSI_eta, NSI_phi)\n",
    "\n",
    "Xe_nucleus.update_model(NSI_model)\n",
    "Xe_nucleus.prepare_density()\n",
    "\n",
    "NSI_spec=Xe_nucleus.spectrum(E_Rs, nu = '8B')\n",
    "acc_data = pd.read_csv('acceptance.csv').drop(columns='Unnamed: 0')\n",
    "\n",
    "accept = np.array(acc_data['acceptance'])\n",
    "energy = np.array(acc_data['energy'] * keV)\n",
    "accept = interp1d(energy, accept)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "b329ac42-973a-42ec-9baa-b538fac95553",
   "metadata": {},
   "outputs": [],
   "source": [
    "Xe_electron = Electron(Xe_nucleus, None, None)\n",
    "Xe_electron.update_model(NSI_model)\n",
    "Xe_electron.prepare_density()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "91bece62-534c-4037-8f1d-d72926aef461",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "90% UL B8 solar neutrino rate: 13.12 per ty\n"
     ]
    }
   ],
   "source": [
    "er = energy\n",
    "tot_rate = (np.abs(Xe_nucleus.spectrum(er, nu = '8B')) * accept(er).astype(\"float64\") * (yr**(-1) * keV**(-1) * ton**(-1)))\n",
    "\n",
    "delta_energy = er[1] - er[0]\n",
    "total_rate = np.sum(tot_rate) * delta_energy\n",
    "upper_lim_B8_flux = 5.25e6 * 2.3 * cm**-2 * sec**-1\n",
    "lower_lim_B8_flux = 5.25e6 * 0.9 * cm**-2 * sec**-1\n",
    "tot_flux = 5.46e6 * cm**-2 * sec**-1\n",
    "upper_lim_total_rate = total_rate / tot_flux * upper_lim_B8_flux\n",
    "lower_lim_total_rate = total_rate / tot_flux * lower_lim_B8_flux\n",
    "print(\"{:} {:.2f} per ty\".format(\"90% UL B8 solar neutrino rate:\", upper_lim_total_rate / (yr**(-1) * ton**(-1))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "82dddd3a-3dfb-4d4b-ba8e-907d7e05bae6",
   "metadata": {},
   "outputs": [],
   "source": [
    "pp_sum = jnp.real(jnp.sum(Xe_electron.spectrum(er_energies, nu = 'pp')))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86cd8232-980d-4354-bd17-6065798aa7bc",
   "metadata": {},
   "source": [
    "## Define the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "deff05a9-1996-46d9-88c3-1d4155ede6e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "import jax\n",
    "import jax.numpy as jnp\n",
    "\n",
    "import numpyro\n",
    "from numpyro import distributions as dist, infer\n",
    "from numpyro import handlers\n",
    "\n",
    "\n",
    "# @handlers.reparam(config={\"phi\": infer.reparam.CircularReparam()})\n",
    "def nsi_rate(rate = None):\n",
    "    eps_ee = numpyro.sample(\"eps_ee\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_em = numpyro.sample(\"eps_em\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_et = numpyro.sample(\"eps_et\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_mm = numpyro.sample(\"eps_mm\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_mt = numpyro.sample(\"eps_mt\", dist.Uniform(-5.0, 5.0))\n",
    "    eps_tt = numpyro.sample(\"eps_tt\", dist.Uniform(-5.0, 5.0))\n",
    "    sineta = numpyro.sample(\"sineta\", dist.Uniform(-1., 1.))\n",
    "    phi = numpyro.sample(\"phi\", dist.Uniform(-0.5 * jnp.pi, 0.5 * jnp.pi))\n",
    "    \n",
    "    NSI_matrix = jnp.asarray([[eps_ee, eps_em, eps_et],\n",
    "                       [eps_em, eps_mm, eps_mt],\n",
    "                       [eps_et, eps_mt, eps_tt]])\n",
    "    NSI_eta = jnp.arcsin(sineta)\n",
    "    NSI_phi = phi\n",
    "\n",
    "    NSI_model = GeneralNSI(NSI_matrix, NSI_eta, NSI_phi)\n",
    "    Xe_nucleus.update_model(NSI_model)\n",
    "    Xe_nucleus.prepare_density()\n",
    "\n",
    "    tot_rate = jnp.sum(jnp.abs(Xe_nucleus.spectrum(er, nu = '8B')) * accept(er) * (yr**(-1) * keV**(-1) * ton**(-1))) * delta_energy\n",
    "    tot_lambda = tot_rate * 0.6 * ton * yr + 5.38\n",
    "    \n",
    "    # Xe_electron = Electron(Xe_nucleus, None, None)\n",
    "    Xe_electron.update_model(NSI_model)\n",
    "    Xe_electron.prepare_density()\n",
    "    pp_sum_NSI = jnp.real(jnp.sum(Xe_electron.spectrum(er_energies, nu = 'pp')))\n",
    "    # pp_10kev_NSI = jnp.real(Xe_electron.spectrum(np.array([1.0e-5]), nu = 'pp'))\n",
    "    \n",
    "    measured_pp_ratio = numpyro.sample(\"measured_pp_ratio\", dist.TruncatedNormal(1.72, 1.72, low=0.0), obs=pp_sum_NSI/pp_sum)\n",
    "    \n",
    "    \n",
    "    numpyro.sample(\"rate\", dist.Poisson(tot_lambda), obs=rate)\n",
    "\n",
    "    \n",
    "    # numpyro.factor('rate_lh',jax.scipy.stats.poisson.logpmf(rate, tot_lambda))\n",
    "\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "01549039-5a3a-4ddd-8c86-50b98f030cfa",
   "metadata": {},
   "source": [
    "## Run samplers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "c65758c1-fafe-4ce7-a769-dd9e1aee47f7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [3:20:41<00:00,  1.72s/it, 31 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:17:27<00:00,  1.18s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:12:05<00:00,  1.13s/it, 31 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:15:45<00:00,  1.16s/it, 63 steps of size 6.00e-02. acc. prob=0.98]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 7h 28min 2s, sys: 2h 40min 59s, total: 10h 9min 1s\n",
      "Wall time: 10h 6min 35s\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    eps_ee      0.06      2.69      0.07     -4.32      4.32  11005.85      1.00\n",
      "    eps_em     -0.02      2.46     -0.01     -3.93      4.31  14296.77      1.00\n",
      "    eps_et      0.02      2.40      0.02     -3.93      4.14  12701.39      1.00\n",
      "    eps_mm      0.02      2.75      0.03     -4.23      4.52  15635.22      1.00\n",
      "    eps_mt      0.02      2.66      0.06     -4.34      4.23  12516.71      1.00\n",
      "    eps_tt     -0.02      2.75     -0.00     -4.57      4.23  12687.41      1.00\n",
      "       phi     -0.00      0.10      0.00     -0.14      0.14   4690.40      1.00\n",
      "    sineta     -0.57      0.04     -0.57     -0.63     -0.51   3074.12      1.00\n",
      "\n",
      "Number of divergences: 10\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(\n",
    "    infer.NUTS(nsi_rate, target_accept_prob=0.98, step_size=0.06, adapt_step_size=False),\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=4,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "time1 = time.time()\n",
    "%time sampler.run(jax.random.PRNGKey(11), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "time2 = time.time()\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "eefd47a5-9090-4e49-ab8e-347c7c28c1c2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/plain_step_fixed_size_0.98_5k_4.json'"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/plain_step_fixed_size_0.98_5k_4.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "14b78f77-4872-40d3-80d8-759e091b01ee",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:15:26<00:00,  1.16s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:16:17<00:00,  1.17s/it, 31 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:19:27<00:00,  1.20s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [3:25:50<00:00,  1.76s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|███████████████████████████████████████████| 7000/7000 [3:22:13<00:00,  1.73s/it, 127 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:59:25<00:00,  1.54s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 12h 17min 16s, sys: 4h 28min 22s, total: 16h 45min 38s\n",
      "Wall time: 16h 39min 34s\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    eps_ee      0.02      2.67      0.04     -4.14      4.51  16171.02      1.00\n",
      "    eps_em     -0.04      2.45     -0.04     -4.26      3.97  20274.39      1.00\n",
      "    eps_et      0.03      2.41      0.03     -4.06      4.07  20152.91      1.00\n",
      "    eps_mm      0.03      2.77      0.04     -4.29      4.51  22517.11      1.00\n",
      "    eps_mt      0.03      2.67      0.02     -4.23      4.41  19077.30      1.00\n",
      "    eps_tt     -0.05      2.72     -0.07     -4.29      4.41  19756.33      1.00\n",
      "       phi      0.00      0.10      0.00     -0.14      0.15   8561.14      1.00\n",
      "    sineta     -0.57      0.04     -0.57     -0.63     -0.52   7224.22      1.00\n",
      "\n",
      "Number of divergences: 4\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(\n",
    "    infer.NUTS(nsi_rate, target_accept_prob=0.98, step_size=0.06, adapt_step_size=False),\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=6,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "time1 = time.time()\n",
    "%time sampler.run(jax.random.PRNGKey(11), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "time2 = time.time()\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "fbe3b355-2ac7-4b8a-9ccc-8e610f4f4c15",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/plain_step_fixed_size_0.98_5k_6.json'"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/plain_step_fixed_size_0.98_5k_6.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "f74007cf-619f-4592-8f12-7bc5dc97a5b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:17:10<00:00,  1.18s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:08:54<00:00,  1.10s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:04:14<00:00,  1.06s/it, 63 steps of size 6.00e-02. acc. prob=0.99]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 6h 28min 33s, sys: 2h 17min 52s, total: 8h 46min 25s\n",
      "Wall time: 8h 45min 15s\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    eps_ee      0.03      2.64     -0.00     -4.12      4.46  11872.91      1.00\n",
      "    eps_em     -0.05      2.45     -0.06     -4.33      3.86  13718.94      1.00\n",
      "    eps_et     -0.02      2.43     -0.04     -4.05      4.14  13942.64      1.00\n",
      "    eps_mm      0.03      2.76      0.06     -4.21      4.57  15032.29      1.00\n",
      "    eps_mt     -0.00      2.68     -0.01     -4.47      4.18  13123.19      1.00\n",
      "    eps_tt     -0.05      2.76     -0.09     -4.71      4.07  13575.80      1.00\n",
      "       phi      0.00      0.10      0.00     -0.15      0.14   3235.14      1.00\n",
      "    sineta     -0.57      0.04     -0.57     -0.62     -0.52   5097.28      1.00\n",
      "\n",
      "Number of divergences: 13\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(\n",
    "    infer.NUTS(nsi_rate, target_accept_prob=0.98, step_size=0.06, adapt_step_size=False),\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=4,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "time1 = time.time()\n",
    "%time sampler.run(jax.random.PRNGKey(10), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "time2 = time.time()\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "355327de-3578-477a-9ace-9c233a301c2f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/plain_step_fixed_size_0.98_5k_4_1.json'"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/plain_step_fixed_size_0.98_5k_4_1.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "abc6573c-18b9-4161-95d6-2bbb13688a0c",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:16:11<00:00,  1.17s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:16:53<00:00,  1.17s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:17:42<00:00,  1.18s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:13:17<00:00,  1.14s/it, 63 steps of size 6.00e-02. acc. prob=0.98]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:15:29<00:00,  1.16s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:07:07<00:00,  1.09s/it, 63 steps of size 6.00e-02. acc. prob=0.99]"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 9h 59min 12s, sys: 3h 30min 2s, total: 13h 29min 14s\n",
      "Wall time: 13h 27min 30s\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    eps_ee     -0.02      2.66      0.02     -4.68      3.97  17042.27      1.00\n",
      "    eps_em     -0.03      2.46     -0.04     -4.15      4.09  21011.05      1.00\n",
      "    eps_et      0.02      2.44      0.01     -3.85      4.35  19984.46      1.00\n",
      "    eps_mm      0.02      2.75      0.03     -4.26      4.51  24833.45      1.00\n",
      "    eps_mt     -0.01      2.67     -0.02     -4.49      4.18  20305.25      1.00\n",
      "    eps_tt     -0.05      2.71     -0.09     -4.47      4.23  22079.37      1.00\n",
      "       phi      0.00      0.10      0.00     -0.15      0.14   6748.51      1.00\n",
      "    sineta     -0.57      0.04     -0.57     -0.63     -0.51   5734.20      1.00\n",
      "\n",
      "Number of divergences: 7\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(\n",
    "    infer.NUTS(nsi_rate, target_accept_prob=0.98, step_size=0.06, adapt_step_size=False),\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=6,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "time1 = time.time()\n",
    "%time sampler.run(jax.random.PRNGKey(9), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "time2 = time.time()\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "8414236d-5ca5-4a5b-897e-22a7d72bf68e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/plain_step_fixed_size_0.98_5k_6_1.json'"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "# az.to_json(arviz_obj, \"nrer_az/plain_step_fixed_size_0.98_5k_6_1.json\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "efe99e5e-29b8-4cea-a653-d6f3b26bcb84",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:22:06<00:00,  1.22s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:19:12<00:00,  1.19s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:15:06<00:00,  1.16s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:19:25<00:00,  1.20s/it, 63 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:19:06<00:00,  1.19s/it, 31 steps of size 6.00e-02. acc. prob=0.99]\n",
      "sample: 100%|████████████████████████████████████████████| 7000/7000 [2:11:48<00:00,  1.13s/it, 31 steps of size 6.00e-02. acc. prob=0.99]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 10h 28min 19s, sys: 3h 20min 50s, total: 13h 49min 9s\n",
      "Wall time: 13h 47min 42s\n",
      "\n",
      "                mean       std    median      5.0%     95.0%     n_eff     r_hat\n",
      "    eps_ee      0.03      2.67      0.02     -4.36      4.26  17304.76      1.00\n",
      "    eps_em     -0.00      2.45     -0.02     -4.13      4.07  21887.66      1.00\n",
      "    eps_et      0.01      2.42      0.02     -4.12      4.00  19493.99      1.00\n",
      "    eps_mm     -0.00      2.78      0.02     -4.53      4.30  22535.82      1.00\n",
      "    eps_mt      0.04      2.71      0.05     -4.40      4.32  18721.92      1.00\n",
      "    eps_tt     -0.05      2.73     -0.07     -4.41      4.31  20707.98      1.00\n",
      "       phi      0.00      0.09      0.00     -0.15      0.14   9002.24      1.00\n",
      "    sineta     -0.57      0.04     -0.57     -0.63     -0.52   9446.46      1.00\n",
      "\n",
      "Number of divergences: 3\n"
     ]
    }
   ],
   "source": [
    "sampler = infer.MCMC(\n",
    "    infer.NUTS(nsi_rate, target_accept_prob=0.98, step_size=0.06, adapt_step_size=False),\n",
    "    num_warmup=2000,\n",
    "    num_samples=5000,\n",
    "    num_chains=6,\n",
    "    progress_bar=True,\n",
    "    chain_method='sequential'\n",
    ")\n",
    "\n",
    "time1 = time.time()\n",
    "%time sampler.run(jax.random.PRNGKey(8), rate = 6, extra_fields = ['num_steps', 'accept_prob'])\n",
    "time2 = time.time()\n",
    "sampler.print_summary()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "093d54d6-f397-48bc-8121-ae3cf244a690",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'nrer_az/plain_step_fixed_size_0.98_5k_6_2.json'"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "arviz_obj = az.from_numpyro(sampler)\n",
    "az.to_json(arviz_obj, \"nrer_az/plain_step_fixed_size_0.98_5k_6_2.json\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:numpyro_gpu]",
   "language": "python",
   "name": "conda-env-numpyro_gpu-py"
  },
  "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.11.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}