DoseNet_TimeBinning (2).ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "a137cd0d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Net counts (tail-trained background): 25953.979995303238\n",
      "Net counts (linear-estimate background): 25770.174999999952\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_88416/3622507761.py:121: RuntimeWarning: overflow encountered in exp\n",
      "  return m * np.exp(-t * x)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(25953.979995303238, 25770.174999999952)"
      ]
     },
     "execution_count": 84,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# integration & subtraction methods\n",
    "# time binning in next cell\n",
    "\n",
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "import scipy.optimize\n",
    "import datetime\n",
    "from scipy.integrate import trapz\n",
    "get_ipython().run_line_magic('matplotlib', 'inline')\n",
    "\n",
    "\n",
    "\n",
    "def setup(filename):\n",
    "    df = pd.read_csv(filename)\n",
    "    df = df.fillna(method = 'ffill') # forward-filling all the NaNs\n",
    "    timedata = df.loc[:, 'deviceTime_unix'] # finding time data\n",
    "    timevals = df['deviceTime_unix'].to_list()\n",
    "    df = df.set_index('deviceTime_unix')\n",
    "    spectra = df.loc[:, '0':'1023'] \n",
    "    array = spectra.to_numpy()\n",
    "    spectrum = array.sum(axis = 0)\n",
    "    x_all = np.array(range(0, 1024))\n",
    "    y_all = spectrum\n",
    "    return x_all, y_all\n",
    "\n",
    "\n",
    "# integration method\n",
    "\n",
    "def arguments (**keywordargs): # seperating p0 and gaussian kwargs \n",
    "    p0 = keywordargs['p0']\n",
    "    gauss_guess = []\n",
    " \n",
    "    gauss_guess.append(keywordargs['amp'])\n",
    "    gauss_guess.append(keywordargs['x0'])\n",
    "    gauss_guess.append(keywordargs['sigma'])\n",
    "    return  p0, gauss_guess\n",
    "\n",
    "def gaussian(x, amp, x0, sigma):\n",
    "    return amp*np.exp(-(x-x0)**2/(2*sigma**2))\n",
    "\n",
    "def gaussian_function(x, a, b, amp, x0, sigma): \n",
    "    return (x*a + b + gaussian(x, amp, x0, sigma))\n",
    "   \n",
    "def NetCounts_integration(x, y, xbounds_left, xbounds_right, **kwargs):\n",
    "    '''\n",
    "    Returns net counts of given peak using Gaussian integration.\n",
    "    \n",
    "        Parameters:\n",
    "            x (array): x values\n",
    "            y (array): y values\n",
    "            xbounds_left (int): left bound of peak\n",
    "            xbounds_right (int): right bound of peak\n",
    "            Keyword arguments:\n",
    "            p0 (list): estimate values for exponential curve fit\n",
    "            amp (int): amplitude of peak\n",
    "            x0 (int): mean of distribution\n",
    "            sigma (float): sigma (variance) value\n",
    "            a (int), Optional: estimate value for base; defaults to 1\n",
    "            b (int), Optional: estimate value for base; defaults to 1\n",
    "            \n",
    "        Returns:\n",
    "            integrated (float): Net counts as determined by Gaussian integration\n",
    "    '''\n",
    "    \n",
    "    r_idx = np.where(x == xbounds_right) # finding index of xbounds_right value    \n",
    "    r_idx = int(r_idx[0][0])\n",
    "    \n",
    "    l_idx = np.where(x== xbounds_left) # finding index of xbounds_left\n",
    "    l_idx = int(l_idx[0][0])\n",
    "    \n",
    "    xpeak = x[l_idx:r_idx] \n",
    "    ypeak = y[l_idx:r_idx] \n",
    "    \n",
    "    p0, gauss_guess = arguments(**kwargs)\n",
    "  \n",
    "    popt, pcov = scipy.optimize.curve_fit(gaussian, xpeak, ypeak, gauss_guess)\n",
    "    print(popt)\n",
    "\n",
    "    gauss_vals = gaussian(xpeak, *popt) # gaussian fit values\n",
    "    plt.figure(figsize = (10, 6)) \n",
    "    plt.plot(xpeak, ypeak, label=\"data\")\n",
    "    plt.plot(xpeak, gauss_vals, '-y', label = 'Gaussian') \n",
    "    plt.legend()\n",
    "    plt.title(\"Gaussian Integration\")\n",
    "    plt.grid()\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "\n",
    "    # integrating\n",
    "    integrated = gauss_vals.sum()\n",
    "    print('Net counts:', integrated)\n",
    "    \n",
    "    return integrated\n",
    "\n",
    "\n",
    "params = {'p0':[1000, 0.05], 'offset':None, 'tail length':0.46}\n",
    "# Note: if curve fit is underestimating drastically, try reducing tail length\n",
    "#       if curve fit is overestimating drastically, try increasing tail length\n",
    "\n",
    "\n",
    "\n",
    "# Using background subtraction \n",
    "\n",
    "def seperate_args (**keywordargs): # seperating p0 and gaussian kwargs \n",
    "    p0 = keywordargs['p0']\n",
    "    offset = keywordargs['offset']\n",
    "    if offset == None:\n",
    "        offset = 0.1\n",
    "    tail_length = keywordargs['tail length']\n",
    "    if tail_length == None:\n",
    "        tail_length = 0.25\n",
    "    plot_all = keywordargs['plot_all']\n",
    "    return p0, offset, tail_length, plot_all\n",
    "\n",
    "def exponential_curve(x, m, t): \n",
    "    return m * np.exp(-t * x) \n",
    "         \n",
    "def NetCounts_subtraction(x, y, xbounds_left, xbounds_right, **kwargs):\n",
    "    '''\n",
    "    Uses two approaches to subtract background counts and return net counts of given peak.\n",
    "    \n",
    "        Parameters:\n",
    "            x (array): x values\n",
    "            y (array): y values\n",
    "            xbounds_left (int): left bound of peak\n",
    "            xbounds_right (int): right bound of peak\n",
    "            Keyword arguments:\n",
    "            p0 (list): estimate values for exponential curve fit\n",
    "            offset (float), Optional: percentage of offset to avoid noise early in data; defaults to 0.1\n",
    "            tail length (float), Optional: percentage of data to use as tail; defaults to 0.25\n",
    "            plot_all (boolean): option to generate visual plots\n",
    "        Returns:\n",
    "            tailBackground_subtracted(float): Net counts as determined by tail-trained background subtraction \n",
    "            straightBackground_subtracted(float): Net counts as determined by subtracting a straight line between bounds\n",
    "    '''\n",
    "    \n",
    "    p0, offset, tail_length, plot_all = seperate_args(**kwargs)\n",
    "    \n",
    "    x = np.array(x)\n",
    "    y = np.array(y)\n",
    "    r_idx = np.where(x == xbounds_right) # finding index of xbounds_right value\n",
    "    r_idx = r_idx[0][0]\n",
    "    l_idx = np.where(x== xbounds_left) # finding index of xbounds_left\n",
    "    l_idx = int(l_idx[0][0])\n",
    "    \n",
    "    xpeak = x[l_idx:r_idx] # peak only\n",
    "    ypeak = y[l_idx:r_idx]\n",
    "    \n",
    "    tail_range = range(int(len(x)*float(tail_length)))\n",
    "    xtail = []\n",
    "    for i in tail_range:\n",
    "        xtail.append(x[l_idx-i])\n",
    "    xtail.reverse()\n",
    "    idx = np.where(x==xtail[0])\n",
    "    tailstart = int(idx[0]) # finding index of tail beginning\n",
    "   \n",
    "\n",
    "    x_roi = x[tailstart:r_idx] # including left tail, until xbounds_right\n",
    "    y_roi = y[tailstart:r_idx]\n",
    "      \n",
    "    leftx = x[tailstart:l_idx] # left tail\n",
    "    lefty = y[tailstart:l_idx]\n",
    "\n",
    "    # training the curve on the left tail \n",
    "  \n",
    "    left_p, left_c = scipy.optimize.curve_fit(exponential_curve, leftx, lefty, p0, maxfev = 1000000)\n",
    "    \n",
    "    # didn't use right-tail or two-tail training because of sensor settling noise\n",
    "   \n",
    "    offset = len(x)*float(offset) # to avoid noise early in the data\n",
    "    offset = int(offset)\n",
    "    \n",
    "    left_tail = exponential_curve(x[offset:], left_p[0], left_p[1])\n",
    "    if plot_all == True:\n",
    "        plt.figure(figsize = (5, 3))\n",
    "        plt.plot(xpeak, ypeak, label=\"data\")\n",
    "        plt.plot(xpeak, left_tail[l_idx:r_idx], '--', label=\"left tail trained\") # peak with left tail trained curve\n",
    "        plt.plot(xpeak, ypeak - left_tail[l_idx:r_idx], '-r', label = 'subtracted') # peak with curve subtracted\n",
    "        plt.legend()\n",
    "        plt.title(\"Background fit (left tail trained)\")\n",
    "        plt.grid()\n",
    "        plt.legend()\n",
    "        plt.show()\n",
    "    \n",
    "    # area under curve, tail-trained background subtracted \n",
    "    tailBackground_subtracted = scipy.integrate.trapz((ypeak - left_tail[l_idx:r_idx]), xpeak)\n",
    "    print('Net counts (tail-trained background):', tailBackground_subtracted)\n",
    "    \n",
    "    # depending on data, left tail training may be inaccurate\n",
    "    # option for a straight line between bounds as background fit\n",
    "    x1 = [x[l_idx], x[r_idx]]\n",
    "    y1 = [y[l_idx], y[r_idx]]\n",
    "    coefs = np.polyfit(x1, y1, 1)\n",
    "    poly = np.poly1d(coefs) \n",
    "    straight_y = poly(xpeak)\n",
    "    if plot_all == True:\n",
    "        plt.figure(figsize = (5, 3))\n",
    "        plt.plot(xpeak, ypeak, label=\"data\")\n",
    "        plt.plot(xpeak, straight_y, 'y--', label ='straight line')\n",
    "        plt.plot(xpeak, ypeak-straight_y, '-g', label = 'subtracted')\n",
    "        plt.legend()\n",
    "        plt.title('Background fit (linear estimator)')\n",
    "        plt.show()\n",
    "\n",
    "    # area under curve, straight background subtracted\n",
    "    straightBackground_subtracted = scipy.integrate.trapz((ypeak-straight_y), xpeak)\n",
    "    print('Net counts (linear-estimate background):', straightBackground_subtracted)\n",
    "\n",
    "    return tailBackground_subtracted, straightBackground_subtracted\n",
    "params = {'p0':[1000, 0.05], 'offset':None, 'tail length':0.46, 'plot_all':False}\n",
    "\n",
    "NetCounts_subtraction(x, y, 800, 920, **params)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "24c485be",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "from timestamp 1489611059 to 1489641001\n",
      "Net counts (tail-trained background): 5139.168521989112\n",
      "Net counts (linear-estimate background): 1147.2166666666574\n",
      "from timestamp 1489641301 to 1489651204\n",
      "Net counts (tail-trained background): 1711.5024055418153\n",
      "Net counts (linear-estimate background): 528.3166666666662\n",
      "from timestamp 1489793103 to 1489822810\n",
      "Net counts (tail-trained background): 4780.261173221473\n",
      "Net counts (linear-estimate background): 475.49999999999386\n",
      "from timestamp 1489823110 to 1489823710\n",
      "Net counts (tail-trained background): 164.95199993438354\n",
      "Net counts (linear-estimate background): 107.96666666666655\n",
      "from timestamp 1489997152 to 1490026860\n",
      "Net counts (tail-trained background): 4953.352546500654\n",
      "Net counts (linear-estimate background): 769.8166666666646\n",
      "from timestamp 1490027160 to 1490056868\n",
      "Net counts (tail-trained background): 5233.290895384511\n",
      "Net counts (linear-estimate background): 1317.5833333333312\n",
      "from timestamp 1490057168 to 1490082975\n",
      "Net counts (tail-trained background): 4614.640835726445\n",
      "Net counts (linear-estimate background): 1365.5833333333317\n",
      "from timestamp 1490256354 to 1490286062\n",
      "Net counts (tail-trained background): 5234.6545141426395\n",
      "Net counts (linear-estimate background): 991.3833333333273\n",
      "from timestamp 1490286362 to 1490316070\n",
      "Net counts (tail-trained background): 4824.735761075474\n",
      "Net counts (linear-estimate background): 943.3999999999817\n",
      "from timestamp 1490316370 to 1490346368\n",
      "Net counts (tail-trained background): 4587.176349268807\n",
      "Net counts (linear-estimate background): 165.04999999996517\n",
      "from timestamp 1490346668 to 1490376374\n",
      "Net counts (tail-trained background): 4603.365266285287\n",
      "Net counts (linear-estimate background): 476.8499999999908\n",
      "from timestamp 1490376674 to 1490406381\n",
      "Net counts (tail-trained background): 4643.703336706852\n",
      "Net counts (linear-estimate background): 319.1833333333245\n",
      "from timestamp 1490406681 to 1490428587\n",
      "Net counts (tail-trained background): 3370.7985385776265\n",
      "Net counts (linear-estimate background): 107.56666666665524\n",
      "from timestamp 1490515573 to 1490545281\n",
      "Net counts (tail-trained background): 4752.330291229336\n",
      "Net counts (linear-estimate background): 183.59999999999377\n",
      "from timestamp 1490545581 to 1490575289\n",
      "Net counts (tail-trained background): 5119.372186548126\n",
      "Net counts (linear-estimate background): 959.016666666661\n",
      "from timestamp 1490575589 to 1490601397\n",
      "Net counts (tail-trained background): 4257.0055834634\n",
      "Net counts (linear-estimate background): 974.2499999999947\n",
      "from timestamp 1490774756 to 1490804464\n",
      "Net counts (tail-trained background): 4914.978782995089\n",
      "Net counts (linear-estimate background): 1302.133333333325\n",
      "from timestamp 1490804764 to 1490834472\n",
      "Net counts (tail-trained background): 4906.504404043858\n",
      "Net counts (linear-estimate background): 1208.0833333333126\n",
      "from timestamp 1490834773 to 1490860580\n",
      "Net counts (tail-trained background): 4204.242182179242\n",
      "Net counts (linear-estimate background): 692.1500000000007\n",
      "from timestamp 1491206754 to 1491236462\n",
      "Net counts (tail-trained background): 5262.109636383189\n",
      "Net counts (linear-estimate background): 1181.2666666666648\n",
      "from timestamp 1491236762 to 1491266471\n",
      "Net counts (tail-trained background): 5134.949881821076\n",
      "Net counts (linear-estimate background): 1134.0166666666555\n",
      "from timestamp 1491266771 to 1491292579\n",
      "Net counts (tail-trained background): 4294.0650571430115\n",
      "Net counts (linear-estimate background): 780.2166666666623\n",
      "from timestamp 1491379563 to 1491409266\n",
      "Net counts (tail-trained background): 5258.9893984014925\n",
      "Net counts (linear-estimate background): 1128.9666666666667\n",
      "from timestamp 1491409566 to 1491439271\n",
      "Net counts (tail-trained background): 4688.996149818265\n",
      "Net counts (linear-estimate background): 1004.599999999981\n",
      "from timestamp 1491439571 to 1491469559\n",
      "Net counts (tail-trained background): 5052.49250580886\n",
      "Net counts (linear-estimate background): 1063.6000000000058\n",
      "from timestamp 1491469859 to 1491499568\n",
      "Net counts (tail-trained background): 5022.579276142338\n",
      "Net counts (linear-estimate background): 754.7166666666641\n",
      "from timestamp 1491499868 to 1491529576\n",
      "Net counts (tail-trained background): 5139.133838833394\n",
      "Net counts (linear-estimate background): 772.6166666666566\n",
      "from timestamp 1491529876 to 1491559874\n",
      "Net counts (tail-trained background): 4615.990690297945\n",
      "Net counts (linear-estimate background): 554.1166666666631\n",
      "from timestamp 1491560174 to 1491589882\n",
      "Net counts (tail-trained background): 4899.034888656014\n",
      "Net counts (linear-estimate background): 663.5999999999867\n",
      "from timestamp 1491590182 to 1491619890\n",
      "Net counts (tail-trained background): 5209.066197441332\n",
      "Net counts (linear-estimate background): 1243.783333333312\n",
      "from timestamp 1491620191 to 1491650164\n",
      "Net counts (tail-trained background): 5116.856340692745\n",
      "Net counts (linear-estimate background): 1299.6833333333325\n",
      "from timestamp 1491650465 to 1491680172\n",
      "Net counts (tail-trained background): 5177.727784515874\n",
      "Net counts (linear-estimate background): 902.9666666666526\n",
      "from timestamp 1491680472 to 1491710180\n",
      "Net counts (tail-trained background): 4825.741088295829\n",
      "Net counts (linear-estimate background): 432.0666666666539\n",
      "from timestamp 1491710480 to 1491740464\n",
      "Net counts (tail-trained background): 5399.325187577698\n",
      "Net counts (linear-estimate background): 1525.6333333333241\n",
      "from timestamp 1491740764 to 1491770471\n",
      "Net counts (tail-trained background): 5237.399336124744\n",
      "Net counts (linear-estimate background): 1114.233333333319\n",
      "from timestamp 1491770771 to 1491791174\n",
      "Net counts (tail-trained background): 3567.2134795031534\n",
      "Net counts (linear-estimate background): 898.083333333332\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_88416/3622507761.py:121: RuntimeWarning: overflow encountered in exp\n",
      "  return m * np.exp(-t * x)\n",
      "/tmp/ipykernel_88416/3622507761.py:121: RuntimeWarning: overflow encountered in multiply\n",
      "  return m * np.exp(-t * x)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEWCAYAAACXGLsWAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAh0UlEQVR4nO3de5hV1X3/8feHi0C8K6NFwIApseGiqIPF2qiJVTGaYI1GjCbYXKgJRqM2ippGa+VXf6k1SuPlodGijTeKplKNiZeosSlRBjsJN1GMKBOIjFgiREEGv/1jr9HDeGb2MM65DOfzep55zj5rr7339+yB85211t5rKyIwMzPrSK9KB2BmZtXPycLMzHI5WZiZWS4nCzMzy+VkYWZmuZwszMwsl5OFWY2TdLOkv610HFbdnCys6kj6vKQGSRskrZb0kKQ/L8NxQ9Ifl/o4lSTpLEn/VVgWEWdHxN9XKibrGZwsrKpIugC4Dvh/wN7AvsCNwMQKhlUWkvpUOgaz9jhZWNWQtCtwJTA1Iu6LiD9ExOaI+M+I+Faq00/SdZJWpZ/rJPVL6973V3Nha0HSLEk3SHpQ0npJT0v6SFr387TJr1KL5jRJAyU9IGmdpNclPSWp6P8ZSddLWinpDUkLJH28YF1vSZdKejEdd4GkoQXxTZX0AvBCKvuqpOXpmHMl7ZPKJel7ktZI+r2kX0sandZ9StKStP/fSvqbIjF+DLgZOCx9xnUF5+WqtHyUpCZJF6XjrJZ0Utr/8ymmSwv22UvStPTZ1kqaLWmPbfrFW4/gZGHV5DCgP/CjDupcBowHxgIHAocC396GY5wO/B2wO7AcmA4QEUek9QdGxE4RcQ9wIdAE1JG1ci4F2psfZ36KaQ/gTuDfJfVP6y5Ix/0UsAvwJeDNgm1PAv4UGCnpk8A/AJ8DBgEvA3enescCRwAfBXYDTgPWpnW3AH8dETsDo4GftQ0wIpYCZwPz0mfcrZ3P8kdkv4fBwHeAfwHOBA4BPg58R9J+qe65Kf4jgX2A/wVuaGe/1oM5WVg12RN4LSJaOqhzBnBlRKyJiGayL/4vbMMx7ouIZ9Ix7iD7gm/PZrIv7A+nFs5T0c5kahHxw4hYGxEtEfFPQD9g/7T6K8C3I2JZZH4VEWsLNv+HiHg9It5Kn+/WiHg2IjYBl5C1BIaleHYG/gRQRCyNiNUFsY6UtEtE/G9EPLsN56TY554eEZvJEtVA4PqIWB8Ri4HFwAGp7l8Dl0VEU4r3CuAUd6ltf5wsrJqsBQbmfNHsQ/bXdquXU1ln/a5g+U1gpw7q/iNZ6+NhSb+RNK29ipIulLQ0dQ+tA3Yl+5IFGAq82MFxVhYsb/X5ImID2XkZHBE/A75P9pf7q5JmStolVf0sWcvlZUlPSjqsg+PlWRsRW9LyW+n11YL1b/Heefsw8KPUVbcOWApsIWuJ2XbEycKqyTxgI1m3RntWkX1Btdo3lQH8AfhQ6wpJf/RBgkl/SV8YEfsBnwYukHR023ppfOJisq6j3VP3zu8BpSorgY90dKiC5a0+n6QdyVpcv00xzYiIQ4BRZN1R30rl8yNiIrAX8B/A7E4cqzusBI6PiN0KfvpHxG+7+ThWYU4WVjUi4vdkfeQ3pEHVD0nqK+l4Sd9N1e4Cvi2pTtLAVP+Had2vgFGSxqbxgiu2MYRXgda+eCSdKOmPJQl4g+wv5i1FttsZaAGagT6SvkM2NtHqB8DfSxqRBqkPkLRnOzHcCfxV+gz9yK4KezoiVkgaJ+lPJfUlS4wbgS2SdpB0hqRdU9dRa6ztfcYhknbo3CnJdTMwXdKHAdLvZbu/cq0WOVlYVYmIa8kGhL9N9uW7EjiH7K9lgKuABuDXwELg2VRGRDxPdjXVo2RXFm11ZVQnXAHclrpUPgeMSPvaQNbquTEiniiy3U+Bh4DnybqQNrJ119K1ZH/pP0z2RX4LMKBYABHxGPC3wL3AarIWyaS0eheyweb/TcdZC1yT1n0BWCHpDbJB7DPb+Yw/Ixtz+J2k19qpsy2uB+aSddWtB35JNlhv2xn54UdmZpbHLQszM8tVsmQh6dZ0U8+iNuXfkLRM0uKCfmgkXZJuRFom6biC8kMkLUzrZqT+YzMzK6NStixmARMKCyR9gmzahgMiYhSpv1XSSLJ+2VFpmxsl9U6b3QRMIes/HtF2n2ZmVnolSxYR8XPg9TbFXwOuTjfvEBFrUvlE4O6I2BQRL5Fd236opEHALhExL90MdTsdX1ZpZmYlUO67LD8KfFzSdLIrRv4mIuaTTSvwy4J6Talsc1puW55r4MCBMWzYsO6I2cysZixYsOC1iKhrW17uZNGHbE6e8cA4YHaaY6bYOER0UF6UpClkXVbsu+++NDQ0fOCAzcxqiaSXi5WX+2qoJrK5eSIingHeIZsSoYlsSoRWQ8juZG1Ky23Li4qImRFRHxH1dXXvS4xmZtZF5U4W/wF8EkDSR4EdgNfIbuqZpGz66eFkA9nPpEnS1ksan66C+iJwf5ljNjOreSXrhpJ0F3AU2cRwTcDlwK3Arely2reByWngerGk2cASsmkTphZMZPY1siurBpDdJftQqWI2M7Pitts7uOvr66PtmMXmzZtpampi48aNFYrKuqJ///4MGTKEvn37VjoUs+2epAURUd+2vKbmnG9qamLnnXdm2LBh+N6+niEiWLt2LU1NTQwfPrzS4ZjVrJqa7mPjxo3sueeeThQ9iCT23HNPtwbNKqymkgXgRNED+XdmVnk1lyzMzGzb1dSYRVvDpj3YrftbcfUJuXV22mknNmzYwKpVqzj33HOZM2dOt8bQXVasWMF///d/8/nPfx6AhoYGbr/9dmbMmPGB9z1r1iyOPfZY9tlnW56GamaVVNPJopL22WefkieKlpYW+vTp2q94xYoV3Hnnne8mi/r6eurr33eBRJfMmjWL0aNHO1lYTeroj9TO/MFZKe6GqpAVK1YwevRoIPvyPPnkk5kwYQIjRozgoosuerfeww8/zGGHHcbBBx/MqaeeyoYNGwC48sorGTduHKNHj2bKlCm0XgJ91FFHcemll3LkkUdy/fXXb3XM5uZmPvvZzzJu3DjGjRvHL37xCwCefPJJxo4dy9ixYznooINYv34906ZN46mnnmLs2LF873vf44knnuDEE08E4IorrmDy5Mkce+yxDBs2jPvuu4+LLrqIMWPGMGHCBDZv3txujHPmzKGhoYEzzjiDsWPH8tZbb7FgwQKOPPJIDjnkEI477jhWr15d2pNvZtvMyaJKNDY2cs8997Bw4ULuueceVq5cyWuvvcZVV13Fo48+yrPPPkt9fT3XXnstAOeccw7z589n0aJFvPXWWzzwwAPv7mvdunU8+eSTXHjhhVsd47zzzuP8889n/vz53HvvvXzlK18B4JprruGGG26gsbGRp556igEDBnD11Vfz8Y9/nMbGRs4///z3xfviiy/y4IMPcv/993PmmWfyiU98goULFzJgwAAefPDBdmM85ZRTqK+v54477qCxsZE+ffrwjW98gzlz5rBgwQK+9KUvcdlll5XqNJtZF7kbqkocffTR7LrrrgCMHDmSl19+mXXr1rFkyRIOP/xwAN5++20OO+wwAB5//HG++93v8uabb/L6668zatQoPv3pTwNw2mmnFT3Go48+ypIlS959/8Ybb7B+/XoOP/xwLrjgAs444wxOPvlkhgwZUnT7Qscffzx9+/ZlzJgxbNmyhQkTsseMjBkzhhUrVuTG2GrZsmUsWrSIY445BoAtW7YwaNCgzp42MysTJ4sq0a9fv3eXe/fuTUtLCxHBMcccw1133bVV3Y0bN/L1r3+dhoYGhg4dyhVXXLHVfQg77rgjAJdddtm7f+U3NjbyzjvvMG/ePAYMGLDV/qZNm8YJJ5zAj3/8Y8aPH8+jjz7a6Xh79epF37593728tVevXrS0tOTG2CoiGDVqFPPmzevMaTKzCnE3VBUbP348v/jFL1i+fDkAb775Js8///y7X7oDBw5kw4YN7Q6UT58+ncbGRhobGwE49thj+f73v//u+tbyF198kTFjxnDxxRdTX1/Pc889x84778z69eu7HHtHMRbue//996e5ufndZLF582YWL17c5eOaWWnUdMuimq88AKirq2PWrFmcfvrpbNq0CYCrrrqKj370o3z1q19lzJgxDBs2jHHjxnVqfzNmzGDq1KkccMABtLS0cMQRR3DzzTdz3XXX8fjjj9O7d29GjhzJ8ccfT69evejTpw8HHnggZ511FgcddNA2xb7bbru1G+NZZ53F2WefzYABA5g3bx5z5szh3HPP5fe//z0tLS1885vfZNSoUdt0PDMrrZqaSHDp0qV87GMfq1BE9kH4d2fbi2q/dLa9iQTdDWVmZrmcLMzMLFfNJYvttdtte+bfmVnl1VSy6N+/P2vXrvWXTw/S+jyL/v37VzoUs5pWU1dDDRkyhKamJpqbmysdim2D1iflmVnl1FSy6Nu3r5+2ZmbWBSXrhpJ0q6Q1khYVWfc3kkLSwIKySyQtl7RM0nEF5YdIWpjWzZCfhGNmVnalHLOYBUxoWyhpKHAM8EpB2UhgEjAqbXOjpN5p9U3AFGBE+nnfPs3MrLRKliwi4ufA60VWfQ+4CCgcZZ4I3B0RmyLiJWA5cKikQcAuETEvslHp24GTShWzmZkVV9aroSR9BvhtRPyqzarBwMqC902pbHBablve3v6nSGqQ1OBBbDOz7lO2ZCHpQ8BlwHeKrS5SFh2UFxURMyOiPiLq6+rquhaomZm9TzmvhvoIMBz4VRqjHgI8K+lQshbD0IK6Q4BVqXxIkXIzMyujsrUsImJhROwVEcMiYhhZIjg4In4HzAUmSeonaTjZQPYzEbEaWC9pfLoK6ovA/eWK2czMMiVrWUi6CzgKGCipCbg8Im4pVjciFkuaDSwBWoCpEbElrf4a2ZVVA4CH0o91UrXPcGlmPUPJkkVEnJ6zflib99OB6UXqNQCjuzU4MzPbJjU1N5SZmXWNk4WZmeVysjAzs1xOFmZmlsvJwszMctXUFOXbG18Wa2bl4paFmZnlcsvCbBu5Rbf98u+2fW5ZmJlZLicLMzPL5WRhZma5nCzMzCyXk4WZmeVysjAzs1xOFmZmlsv3WZhZTfA9FB+MWxZmZpbLycLMzHKVLFlIulXSGkmLCsr+UdJzkn4t6UeSditYd4mk5ZKWSTquoPwQSQvTuhmSVKqYzcysuFK2LGYBE9qUPQKMjogDgOeBSwAkjQQmAaPSNjdK6p22uQmYAoxIP233aWZmJVayAe6I+LmkYW3KHi54+0vglLQ8Ebg7IjYBL0laDhwqaQWwS0TMA5B0O3AS8FCp4jb7oDyQatujSo5ZfIn3vvQHAysL1jWlssFpuW15UZKmSGqQ1NDc3NzN4ZqZ1a6KJAtJlwEtwB2tRUWqRQflRUXEzIioj4j6urq6Dx6omZkBFbjPQtJk4ETg6Iho/eJvAoYWVBsCrErlQ4qUm23X3JVl1aasLQtJE4CLgc9ExJsFq+YCkyT1kzScbCD7mYhYDayXND5dBfVF4P5yxmxmZiVsWUi6CzgKGCipCbic7OqnfsAj6QrYX0bE2RGxWNJsYAlZ99TUiNiSdvU1siurBpCNcXhw28yszEp5NdTpRYpv6aD+dGB6kfIGYHQ3hmZmZtvId3CbmVkuTyRYhAcXzcy25paFmZnlcsvCrIdyC9jKycnCrIY54VhnuRvKzMxyuWVhZlZlqrHF55aFmZnlcrIwM7NcThZmZpbLycLMzHI5WZiZWS5fDWVmH0g1Xrlj3c8tCzMzy+VkYWZmuZwszMwsl5OFmZnlcrIwM7NcJUsWkm6VtEbSooKyPSQ9IumF9Lp7wbpLJC2XtEzScQXlh0hamNbNUHp4t5mZlU8pWxazgAltyqYBj0XECOCx9B5JI4FJwKi0zY2SeqdtbgKmACPST9t9mplZiZXsPouI+LmkYW2KJwJHpeXbgCeAi1P53RGxCXhJ0nLgUEkrgF0iYh6ApNuBk4CHShW3mb3H91BYq3KPWewdEasB0uteqXwwsLKgXlMqG5yW25YXJWmKpAZJDc3Nzd0auJlZLauWO7iLjUNEB+VFRcRMYCZAfX19u/XsPf7L0cw6o9wti1clDQJIr2tSeRMwtKDeEGBVKh9SpNzMzMqo3MliLjA5LU8G7i8onySpn6ThZAPZz6SuqvWSxqeroL5YsI2ZmZVJybqhJN1FNpg9UFITcDlwNTBb0peBV4BTASJisaTZwBKgBZgaEVvSrr5GdmXVALKBbQ9uW8m4W86suFJeDXV6O6uObqf+dGB6kfIGYHQ3hmZmZtuoWga4zcrCLQezrvF0H2ZmlsstC8vlv8bNzC0LMzPL5WRhZma5nCzMzCzXNicLSbtLOqAUwZiZWXXq1AC3pCeAz6T6jUCzpCcj4oLShWZmVl1q+WKPzl4NtWtEvCHpK8C/RsTlkn5dysDMzKy4SiStznZD9UkT/30OeKAkkZiZWdXqbLL4O+CnwPKImC9pP+CF0oVlZmbVpLPdUKsj4t1B7Yj4jaRrSxSTmZlVmc62LP65k2VmZrYd6rBlIekw4M+AOkmFVz7tAvQuZWBmZlY98rqhdgB2SvV2Lih/AzilVEGZmVl16TBZRMSTwJOSZkXEy2WKyczMqkxnB7j7SZoJDCvcJiI+WYqgzMysunQ2Wfw7cDPwA2BLTl0zM9vOdPZqqJaIuCkinomIBa0/XT2opPMlLZa0SNJdkvpL2kPSI5JeSK+7F9S/RNJyScskHdfV45qZWdd0Nln8p6SvSxqUvtT3kLRHVw4oaTBwLlAfEaPJrqqaBEwDHouIEcBj6T2SRqb1o4AJwI2SfCWWmVkZdbYbanJ6/VZBWQD7fYDjDpC0GfgQsAq4BDgqrb8NeAK4GJgI3B0Rm4CXJC0HDgXmdfHYZma2jTqVLCJieHcdMCJ+K+ka4BXgLeDhiHhY0t4RsTrVWS1pr7TJYOCXBbtoSmXvI2kKMAVg33337a6QzcxqXmenKP9isfKIuH1bD5jGIiYCw4F1wL9LOrOjTYodup14ZgIzAerr64vWMTOzbdfZbqhxBcv9gaOBZ4FtThbAXwAvRUQzgKT7yO4Sf1XSoNSqGASsSfWbgKEF2w8h67YyM7My6Ww31DcK30vaFfi3Lh7zFWC8pA+RdUMdDTQAfyAbG7k6vd6f6s8F7kwTF+4DjACe6eKxzWpGLT+ox7pfZ1sWbb1J9qW9zSLiaUlzyFomLcD/kHUd7QTMlvRlsoRyaqq/WNJsYEmqPzUiauJeD/9nt1rhf+vVr7NjFv/Je+MEvYGPAbO7etCIuBy4vE3xJrJWRrH604HpXT2emZl9MJ1tWVxTsNwCvBwRTSWIx8y2M241bB86dVNemlDwObKZZ3cH3i5lUGZmVl06lSwkfY5sUPlUsudwPy3JU5SbmdWIznZDXQaMi4g1AJLqgEeBOaUKzMzMqkdn54bq1ZookrXbsK2ZmfVwnW1Z/ETST4G70vvTgB+XJiQzM6s2ec/g/mNg74j4lqSTgT8nm35jHnBHGeIzM7MqkNeVdB2wHiAi7ouICyLifLJWxXWlDc3MzKpFXrIYFhG/blsYEQ1kj1g1M7MakDdm0b+DdQO6MxCzWuIb1aynyWtZzJf01baFaf6mLj9W1czMepa8lsU3gR9JOoP3kkM9sAPwlyWMy8zMqkiHySIiXgX+TNIngNGp+MGI+FnJIzMzs6rR2edZPA48XuJYzMy6zONApeW7sM3MLJeThZmZ5XKyMDOzXE4WZmaWqyLJQtJukuZIek7SUkmHSdpD0iOSXkivuxfUv0TScknLJB1XiZjNzGpZpVoW1wM/iYg/AQ4ElgLTgMciYgTwWHqPpJHAJGAUMAG4UVLvikRtZlajyp4sJO0CHAHcAhARb0fEOmAicFuqdhtwUlqeCNwdEZsi4iVgOXBoOWM2M6t1lWhZ7Ac0A/8q6X8k/UDSjmRToa8GSK97pfqDgZUF2zelsveRNEVSg6SG5ubm0n0CM7MaU4lk0Qc4GLgpIg4C/kDqcmqHipRFsYoRMTMi6iOivq6u7oNHamZmQGWSRRPQFBFPp/dzyJLHq5IGAaTXNQX1hxZsPwRYVaZYzcyMCiSLiPgdsFLS/qnoaGAJMBeYnMomA/en5bnAJEn9JA0HRgDPlDFkM7Oa19lncHe3bwB3SNoB+A3wV2SJa3aa/vwV4FSAiFgsaTZZQmkBpkbElsqEbWZWmyqSLCKikWyq87aObqf+dGB6KWMyM7P2+Q5uMzPLValuKDOzTvP045XnloWZmeVysjAzs1zuhuoCN4nNrNa4ZWFmZrmcLMzMLJeThZmZ5XKyMDOzXE4WZmaWy8nCzMxy+dLZEvHltVvz+TDr2dyyMDOzXE4WZmaWy8nCzMxyOVmYmVkuJwszM8vlq6HMzLrR9nrlX8VaFpJ6S/ofSQ+k93tIekTSC+l194K6l0haLmmZpOMqFbOZWa2qZDfUecDSgvfTgMciYgTwWHqPpJHAJGAUMAG4UVLvMsdqZlbTKpIsJA0BTgB+UFA8EbgtLd8GnFRQfndEbIqIl4DlwKFlCtXMzKhcy+I64CLgnYKyvSNiNUB63SuVDwZWFtRrSmVmZlYmZU8Wkk4E1kTEgs5uUqQs2tn3FEkNkhqam5u7HKOZmW2tEi2Lw4HPSFoB3A18UtIPgVclDQJIr2tS/SZgaMH2Q4BVxXYcETMjoj4i6uvq6koVv5lZzSl7soiISyJiSEQMIxu4/llEnAnMBSanapOB+9PyXGCSpH6ShgMjgGfKHLaZWU2rpvssrgZmS/oy8ApwKkBELJY0G1gCtABTI2JL5cI0M6s9FU0WEfEE8ERaXgsc3U696cD0sgVmZmZb8XQfZmaWy8nCzMxyOVmYmVkuJwszM8vlZGFmZrmq6dJZq2Hb67TOZtsLtyzMzCyXk4WZmeVysjAzs1xOFmZmlsvJwszMcjlZmJlZLicLMzPL5WRhZma5nCzMzCyXk4WZmeVysjAzs1xOFmZmlqvsyULSUEmPS1oqabGk81L5HpIekfRCet29YJtLJC2XtEzSceWO2cys1lWiZdECXBgRHwPGA1MljQSmAY9FxAjgsfSetG4SMAqYANwoqXcF4jYzq1llTxYRsToink3L64GlwGBgInBbqnYbcFJangjcHRGbIuIlYDlwaFmDNjOrcRUds5A0DDgIeBrYOyJWQ5ZQgL1StcHAyoLNmlKZmZmVScWShaSdgHuBb0bEGx1VLVIW7exziqQGSQ3Nzc3dEaaZmVGhZCGpL1miuCMi7kvFr0oalNYPAtak8iZgaMHmQ4BVxfYbETMjoj4i6uvq6koTvJlZDSr7Y1UlCbgFWBoR1xasmgtMBq5Or/cXlN8p6VpgH2AE8Ez5Iraewo9mNSudSjyD+3DgC8BCSY2p7FKyJDFb0peBV4BTASJisaTZwBKyK6mmRsSWskdtZlbDyp4sIuK/KD4OAXB0O9tMB6aXLCj7QPwXvdn2z3dwm5lZLicLMzPL5WRhZma5nCzMzCyXk4WZmeVysjAzs1xOFmZmlqsSN+UZvjehK3zOzCrHLQszM8vlZGFmZrmcLMzMLJeThZmZ5XKyMDOzXE4WZmaWy8nCzMxyOVmYmVkuJwszM8vlZGFmZrmcLMzMLFePSRaSJkhaJmm5pGmVjsfMrJb0iGQhqTdwA3A8MBI4XdLIykZlZlY7ekSyAA4FlkfEbyLibeBuYGKFYzIzqxmKiErHkEvSKcCEiPhKev8F4E8j4pw29aYAU9Lb/YFl3XD4gcBr3bCfcugpsfaUOKHnxNpT4oSeE2tPiRO6N9YPR0Rd28Ke8jwLFSl7X5aLiJnAzG49sNQQEfXduc9S6Smx9pQ4oefE2lPihJ4Ta0+JE8oTa0/phmoChha8HwKsqlAsZmY1p6cki/nACEnDJe0ATALmVjgmM7Oa0SO6oSKiRdI5wE+B3sCtEbG4TIfv1m6tEuspsfaUOKHnxNpT4oSeE2tPiRPKEGuPGOA2M7PK6indUGZmVkFOFmZmlsvJogM9ZYoRSSskLZTUKKmh0vEUknSrpDWSFhWU7SHpEUkvpNfdKxljq3ZivULSb9O5bZT0qUrGmGIaKulxSUslLZZ0XiqvqvPaQZzVeE77S3pG0q9SrH+XyqvtnLYXZ8nPqccs2pGmGHkeOIbs0t35wOkRsaSigRUhaQVQHxFVdwORpCOADcDtETE6lX0XeD0irk5JePeIuLiScaa4isV6BbAhIq6pZGyFJA0CBkXEs5J2BhYAJwFnUUXntYM4P0f1nVMBO0bEBkl9gf8CzgNOprrOaXtxTqDE59Qti/Z5ipFuEBE/B15vUzwRuC0t30b2BVJx7cRadSJidUQ8m5bXA0uBwVTZee0gzqoTmQ3pbd/0E1TfOW0vzpJzsmjfYGBlwfsmqvQfOtk/loclLUhTnlS7vSNiNWRfKMBeFY4nzzmSfp26qaqiy6yVpGHAQcDTVPF5bRMnVOE5ldRbUiOwBngkIqrynLYTJ5T4nDpZtK9TU4xUicMj4mCyWXmnpu4U6x43AR8BxgKrgX+qaDQFJO0E3At8MyLeqHQ87SkSZ1We04jYEhFjyWaIOFTS6AqHVFQ7cZb8nDpZtK/HTDESEavS6xrgR2RdaNXs1dSf3dqvvabC8bQrIl5N/znfAf6FKjm3qb/6XuCOiLgvFVfdeS0WZ7We01YRsQ54gmwcoOrOaavCOMtxTp0s2tcjphiRtGMaPETSjsCxwKKOt6q4ucDktDwZuL+CsXSo9Ysi+Uuq4NymQc5bgKURcW3Bqqo6r+3FWaXntE7Sbml5APAXwHNU3zktGmc5zqmvhupAuvzsOt6bYmR6ZSN6P0n7kbUmIJu+5c5qilPSXcBRZFMovwpcDvwHMBvYF3gFODUiKj6w3E6sR5E17QNYAfx1ax92pUj6c+ApYCHwTiq+lGw8oGrOawdxnk71ndMDyAawe5P9ET07Iq6UtCfVdU7bi/PfKPE5dbIwM7Nc7oYyM7NcThZmZpbLycLMzHI5WZiZWS4nCzMzy+VkYdaGpD0LZu/8XcFsnhsk3ViC482SdEqR8npJM7r7eGZd0SMeq2pWThGxluya9YrOOhsRDUBVTTlvtcstC7NOknSUpAfS8hWSbpP0sLLniZws6bvKnivykzTNBZIOkfRkmuTxp23utC30F5KekvS8pBPbOd6tkp6Q9BtJ56byHSU9qOz5BosknVaGU2E1yMnCrOs+ApxANo31D4HHI2IM8BZwQkoY/wycEhGHALcC7d1dPww4Mu3vZkn9i9T5E+A4snl/Lk/7nwCsiogD0zM4ftJdH86skLuhzLruoYjYLGkh2fQLrV/UC8m+/PcHRgOPZNMk0ZtsRtBiZqdJ4F6Q9BuyxNDWgxGxCdgkaQ2wdzrWNZL+P/BARDzVPR/NbGtOFmZdtwkgIt6RtDnemzvnHbL/WwIWR8RhndhX23l3is3Ds6lgeQvQJyKel3QI8CngHyQ9HBFXbtOnMOsEd0OZlc4yoE7SYZBN1y1pVDt1T5XUS9JHgP3Strkk7QO8GRE/BK4BDu6GuM3exy0LsxKJiLfTJbEzJO1K9v/tOmBxkerLgCfJupbOjoiNqesqzxjgHyW9A2wGvtYdsZu15Vlnzcwsl7uhzMwsl5OFmZnlcrIwM7NcThZmZpbLycLMzHI5WZiZWS4nCzMzy/V/Q6OTdZH9C3kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZIAAAEWCAYAAABMoxE0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeSklEQVR4nO3de5zVdb3v8ddbQFAEQxqNi4kWpngjGbxv82RHKD1hhonHEsoiTbsfd2juUtue2OVumxYYhoGFuclMKbU0vPaQow5qAhJCijJCMGIoeCHAz/nj951cDmuGNfxmrTWLeT8fj3ms3/r+bp/fD2a95/f9XZYiAjMzs+21U7ULMDOz2uYgMTOzXBwkZmaWi4PEzMxycZCYmVkuDhIzM8vFQWJmrZJ0raR/q3Yd1rk5SKymSPrfkhokbZC0StKdko6rwHpD0nvLvZ5qkjRB0p8K2yLi3Ij4TrVqstrgILGaIelrwFXA/wX2At4NTAHGVLGsipDUvdo1mLXGQWI1QdLuwOXA+RFxS0S8GhGbIuK3EXFhmqanpKskrUw/V0nqmcZt9dd24VGGpBmSfizpdknrJT0s6T1p3ANplj+nI6EzJL1T0u8krZP0kqQHJRX9fZL0Q0krJL0iab6kfykY103SxZL+mtY7X9LeBfWdL2kpsDS1fU7SsrTOOZIGpnZJ+i9JayS9LOlJSQencR+R9FRa/guS/k+RGg8ErgWOTtu4rmC//HsaPkFSo6R/TetZJenUtPynU00XFyxzJ0mT0ratlTRb0h7t+oe3muAgsVpxNNAL+E0b03wTOAoYDhwGHAFc0o51nAlcBvQDlgFXAETE8Wn8YRGxW0T8N/B1oBGoIzs6uhho7XlDj6aa9gBuBH4lqVca97W03o8AfYHPAK8VzHsqcCQwTNIHge8CnwAGAM8BN6XpTgKOB/YH3gGcAaxN46YDn4+IPsDBwD0tC4yIxcC5wLy0je9oZVveRfbvMAj4FnAd8ElgBPAvwLck7Zem/VKq/wPAQODvwI9bWa7VMAeJ1Yr+wIsRsbmNac4CLo+INRHRRBYKn2rHOm6JiEfSOmaRffi3ZhPZh/k+6cjowWjlwXUR8YuIWBsRmyPiP4GewPvS6M8Cl0TEksj8OSLWFsz+3Yh4KSJeT9t3fUQ8FhEbgYvIjiCGpHr6AAcAiojFEbGqoNZhkvpGxN8j4rF27JNi231FRGwiC7F3Aj+MiPURsQhYBByapv088M2IaEz1XgqMdTfdjsdBYrViLfDObXwIDST7K73Zc6mtVH8rGH4N2K2Nab9PdtRyl6RnJE1qbUJJX5e0OHU5rQN2J/sABtgb+Gsb61lRMPy27YuIDWT7ZVBE3AP8iOwv/tWSpknqmyb9ONkRz3OS7pd0dBvr25a1EbElDb+eXlcXjH+dt/bbPsBvUvffOmAxsIXsCM52IA4SqxXzgDfIukpas5Lsw6vZu1MbwKvArs0jJL0rTzHpL/CvR8R+wP8CvibpxJbTpfMh3yDrjuqXuoxeBpQmWQG8p61VFQy/bfsk9SY7Unsh1XR1RIwADiLr4rowtT8aEWOAPYFbgdklrKsjrAA+HBHvKPjpFREvdPB6rMocJFYTIuJlsj75H6cTvLtK6iHpw5K+lyb7JXCJpDpJ70zT/yKN+zNwkKTh6fzEpe0sYTXQ3PePpFMkvVeSgFfI/tLeUmS+PsBmoAnoLulbZOdCmv0U+I6koemE+aGS+rdSw43Ap9M29CS7eu3hiFguaaSkIyX1IAvNN4AtknaWdJak3VN3VHOtrW3jYEk7l7ZLtula4ApJ+wCkf5cd/gq7rshBYjUjIn5AdnL6ErIP5hXABWR/ZQP8O9AAPAksAB5LbUTE02RXff2R7Aqot13BVYJLgZmpm+YTwNC0rA1kR0tTIuK+IvP9AbgTeJqsW+oN3t5d9QOyI4S7yD7kpwO7FCsgIuYC/wb8GlhFdiQzLo3uS3bi++9pPWuBK9O4TwHLJb1CdkL9k61s4z1k5zj+JunFVqZpjx8Cc8i6/9YD/4/swgHbwchfbGVmZnn4iMTMzHIpW5BIuj7dtLSwoG0PSXdLWppe+xWMuyjdaLVE0qiC9hGSFqRxV6c+6eabz/47tT+cLoE0M7MKK+cRyQxgdIu2ScDciBgKzE3vkTSMrK/3oDTPFEnd0jxTgYlkfdJDC5Z5DvD3iHgv8F/Af5RtS8zMrFVlC5KIeAB4qUXzGGBmGp7JW5dyjgFuioiNEfEs2fX5R0gaAPSNiHnpZq8bWszTvKybgRObj1bMzKxyKn2H6V7Nd9tGxCpJe6b2QWRXdDRrTG2b0nDL9uZ5VqRlbZb0Munu55YrlTSR7KiG3r17jzjggAM6bIPMzLqC+fPnvxgRdcXGdZZHFRQ7kog22tuaZ+vGiGnANID6+vpoaGjYnhrNzLosSc+1Nq7SV22tTt1VpNc1qb2R7FERzQaT3cXbmIZbtr9tnvTYjN3ZuivNzMzKrNJBMgcYn4bHA7cVtI9LV2LtS3ZS/ZHUDbZe0lHp/MfZLeZpXtZY4J7WHppnZmblU7auLUm/BE4ge9BeI/BtYDIwW9I5wPPA6QARsUjSbOApssdJnF/wYLjzyK4A24XsDuE7U/t04OeSlpEdiTTf4WtmZhXU5e5s9zkSs9q0adMmGhsbeeONN6pdyg6tV69eDB48mB49erytXdL8iKgvNk9nOdluZtamxsZG+vTpw5AhQ/CV/uUREaxdu5bGxkb23XffkufzI1LMrCa88cYb9O/f3yFSRpLo379/u4/6HCRmVjMcIuW3PfvYQWJmZrn4HImZ1aQhk27v0OUtn3xym+PXrVvHjTfeyBe+8IU2pzvmmGN46KGHWL58OaeccgoLFy7capoZM2Zw0kknMXBge74JGq699lp23XVXzj777HbNV8yECRM45ZRTGDt2bO5l+YjEzKwE69atY8qUKduc7qGHHtrmNDNmzGDlypVFx23Z0toXWMK5557bISHS0RwkZmYlmDRpEn/9618ZPnw4X/3qVznxxBM5/PDDOeSQQ7jtttv+Od1uu+3W5nJuvvlmGhoaOOussxg+fDivv/46Q4YM4fLLL+e4447jV7/6Fddddx0jR47ksMMO4+Mf/zivvfYaAJdeeilXXpl98eUJJ5zAN77xDY444gj2339/HnzwQSALogsvvJCRI0dy6KGH8pOf/ATIrsi64IILGDZsGCeffDJr1qwpXuB2cNeWmVkJJk+ezMKFC3niiSfYvHkzr732Gn379uXFF1/kqKOO4qMf/WhJJ6rHjh3Lj370I6688krq69+6LaNXr1786U/ZN0CvXbuWz33ucwBccsklTJ8+nS9+8YtbLWvz5s088sgj3HHHHVx22WX88Y9/ZPr06ey+++48+uijbNy4kWOPPZaTTjqJxx9/nCVLlrBgwQJWr17NsGHD+MxnPtMh+8ZBYmbWThHBxRdfzAMPPMBOO+3ECy+8wOrVq3nXu9613cs844wz/jm8cOFCLrnkEtatW8eGDRsYNWpU0XlOO+00AEaMGMHy5csBuOuuu3jyySe5+eabAXj55ZdZunQpDzzwAGeeeSbdunVj4MCBfPCDH9zuWltykJiZtdOsWbNoampi/vz59OjRgyFDhrR578WnP/1pHn/8cQYOHMgdd9xRdJrevXv/c3jChAnceuutHHbYYcyYMYP77ruv6Dw9e/YEoFu3bmzevBnIQu6aa67ZKnzuuOOOsl0+7XMkZmYl6NOnD+vXrweyv/L33HNPevTowb333stzz7X6hHUAfvazn/HEE0/8M0QKl1XM+vXrGTBgAJs2bWLWrFntqnPUqFFMnTqVTZs2AfD000/z6quvcvzxx3PTTTexZcsWVq1axb333tuu5bbFRyRmVpO2dbluR+vfvz/HHnssBx98MCNHjuQvf/kL9fX1DB8+nPZ+Wd6ECRM499xz2WWXXZg3b95W47/zne9w5JFHss8++3DIIYe0GTotffazn2X58uUcfvjhRAR1dXXceuutfOxjH+Oee+7hkEMOYf/99+cDH/hAu2puix/aaGY1YfHixRx44IHVLqNLKLav23poo7u2zMwsFweJmZnl4iAxs5rR1briq2F79rGDxMxqQq9evVi7dq3DpIyav4+kV69e7ZrPV22ZWU0YPHgwjY2NNDU1VbuUHVrzNyS2h4PEzGpCjx492vWtfVY57toyM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsFweJmZnl4iAxM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsFweJmZnlUpUgkfRVSYskLZT0S0m9JO0h6W5JS9Nrv4LpL5K0TNISSaMK2kdIWpDGXS1J1dgeM7OurOJBImkQ8CWgPiIOBroB44BJwNyIGArMTe+RNCyNPwgYDUyR1C0tbiowERiafkZXcFPMzIzqdW11B3aR1B3YFVgJjAFmpvEzgVPT8BjgpojYGBHPAsuAIyQNAPpGxLyICOCGgnnMzKxCKh4kEfECcCXwPLAKeDki7gL2iohVaZpVwJ5plkHAioJFNKa2QWm4ZftWJE2U1CCpoampqSM3x8ysy6tG11Y/sqOMfYGBQG9Jn2xrliJt0Ub71o0R0yKiPiLq6+rq2luymZm1oRpdWx8Cno2IpojYBNwCHAOsTt1VpNc1afpGYO+C+QeTdYU1puGW7WZmVkHVCJLngaMk7ZqusjoRWAzMAcanacYDt6XhOcA4ST0l7Ut2Uv2R1P21XtJRaTlnF8xjZmYV0r3SK4yIhyXdDDwGbAYeB6YBuwGzJZ1DFjanp+kXSZoNPJWmPz8itqTFnQfMAHYB7kw/ZmZWQcoueOo66uvro6GhodplmJnVFEnzI6K+2Djf2W5mZrk4SMzMLBcHiZmZ5eIgMTOzXBwkZmaWi4PEzMxycZCYmVkuDhIzM8vFQWJmZrk4SMzMLBcHiZmZ5eIgMTOzXBwkZmaWi4PEzMxycZCYmVkuDhIzM8vFQWJmZrk4SMzMLBcHiZmZ5eIgMTOzXBwkZmaWi4PEzMxycZCYmVkuDhIzM8vFQWJmZrk4SMzMLBcHiZmZ5eIgMTOzXBwkZmaWi4PEzMxycZCYmVkuDhIzM8vFQWJmZrk4SMzMLJeqBImkd0i6WdJfJC2WdLSkPSTdLWlpeu1XMP1FkpZJWiJpVEH7CEkL0rirJaka22Nm1pVV64jkh8DvI+IA4DBgMTAJmBsRQ4G56T2ShgHjgIOA0cAUSd3ScqYCE4Gh6Wd0JTfCzMyqECSS+gLHA9MBIuIfEbEOGAPMTJPNBE5Nw2OAmyJiY0Q8CywDjpA0AOgbEfMiIoAbCuYxM7MKqcYRyX5AE/AzSY9L+qmk3sBeEbEKIL3umaYfBKwomL8xtQ1Kwy3btyJpoqQGSQ1NTU0duzVmZl1cNYKkO3A4MDUi3g+8SurGakWx8x7RRvvWjRHTIqI+Iurr6uraW6+ZmbWhGkHSCDRGxMPp/c1kwbI6dVeRXtcUTL93wfyDgZWpfXCRdjMzq6CKB0lE/A1YIel9qelE4ClgDjA+tY0HbkvDc4BxknpK2pfspPojqftrvaSj0tVaZxfMY2ZmFdK9Suv9IjBL0s7AM8CnyUJttqRzgOeB0wEiYpGk2WRhsxk4PyK2pOWcB8wAdgHuTD9mZlZByi546jrq6+ujoaGh2mWYmdUUSfMjor7YON/ZbmZmubQ7SCT1k3RoOYoxM7PaU1KQSLpPUl9JewB/JrsH5AflLc3MzGpBqUcku0fEK8BpwM8iYgTwofKVZWZmtaLUIOme7u34BPC7MtZjZmY1ptQguQz4A7AsIh6VtB+wtHxlmZlZrSj1PpJVEfHPE+wR8YzPkZiZGZQeJNeQPcZkW21mZradhky6vdVxyyefXMFK2qfNIJF0NHAMUCfpawWj+gLdis9lZmZdybaOSHYGdkvT9SlofwUYW66izMysdrQZJBFxP3C/pBkR8VyFajKzNlSi+6NWu1isOko9R9JT0jRgSOE8EfHBchRl+XXEB4E/TMw6l876O1lqkPwKuBb4KbBlG9OamdUU/+GVT6lBsjkippa1khrQlf+j7Oh2pH9bb4tVWqlB8ltJXwB+A2xsboyIl8pSlZlZJ+JAa1upQdL8zYUXFrQFsF/HlmNmOwp/+HYdJQVJROxb7kLMrOtx2OwYSgoSSWcXa4+IGzq2HCuFf/nMrDMptWtrZMFwL+BE4DHAQWKWlBLw/iOg8rzPy6/Urq0vFr6XtDvw87JUZGZmNaXUI5KWXgOGdmQhVnv8l55Z6Xbk35dSz5H8luwqLcge1nggMLtcRe3IduT/TJ2Z97tZ+ZR6RHJlwfBm4LmIaCxDPWZmVmNKPUdyv6S9eOuku78d0ayL81GeNSu1a+sTwPeB+wAB10i6MCJuLmNtVka19CFQS7Va5fn/R/WV2rX1TWBkRKwBkFQH/BFwkFibtvVL7g8Bs9q3U6nTNYdIsrYd85qZ2Q6s1COS30v6A/DL9P4M4I7ylGRmZrVkW9/Z/l5gr4i4UNJpwHFk50jmAbMqUJ+ZmXVy2+qeugpYDxARt0TE1yLiq2RHI1eVtzQzM6sF2wqSIRHxZMvGiGgg+9pdMzPr4rZ1jqRXG+N26chCzKxyfLWcdaRtBcmjkj4XEdcVNko6B5hfvrLMOh9/+JoVt60g+QrwG0ln8VZw1AM7Ax8rY11mXZLDympRm0ESEauBYyT9D+Dg1Hx7RNxT9srMzKwmlHRTYUTcGxHXpJ8OCRFJ3SQ9Lul36f0eku6WtDS99iuY9iJJyyQtkTSqoH2EpAVp3NWS1BG1mZlZ6ap5d/qXgcUF7ycBcyNiKDA3vUfSMGAccBAwGpgiqVuaZyowkey7UYam8WZmVkFVCRJJg4GTgZ8WNI8BZqbhmcCpBe03RcTGiHgWWAYcIWkA0Dci5kVEkH3t76mYmVlFVeuI5CrgX4E3C9r2iohVAOl1z9Q+CFhRMF1jahuUhlu2b0XSREkNkhqampo6ZAPMzCxT8SCRdAqwJiJKvXy42HmPaKN968aIaRFRHxH1dXV1Ja7WzMxKsb3f2Z7HscBHJX2E7IbHvpJ+AayWNCAiVqVuq+anDTcCexfMPxhYmdoHF2k3M7MKqvgRSURcFBGDI2II2Un0eyLik8AcYHyabDxwWxqeA4yT1FPSvmQn1R9J3V/rJR2VrtY6u2AeMzOrkGockbRmMjA73TX/PHA6QEQskjQbeIrs++LPj4gtaZ7zgBlkj2u5M/2YmVkFVTVIIuI+sq/vJSLWAie2Mt0VwBVF2ht460ZJMzOrAn/LoZmZ5eIgMTOzXBwkZmaWS2c62W5mZjlV4wnSPiIxM7NcHCRmZpaLg8TMzHJxkJiZWS4+2W41z19Pa1ZdPiIxM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsFweJmZnl4iAxM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsFweJmZnl4iAxM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsFweJmZnl4iAxM7NcKh4kkvaWdK+kxZIWSfpyat9D0t2SlqbXfgXzXCRpmaQlkkYVtI+QtCCNu1qSKr09ZmZdXTWOSDYDX4+IA4GjgPMlDQMmAXMjYigwN70njRsHHASMBqZI6paWNRWYCAxNP6MruSFmZlaFIImIVRHxWBpeDywGBgFjgJlpspnAqWl4DHBTRGyMiGeBZcARkgYAfSNiXkQEcEPBPGZmViFVPUciaQjwfuBhYK+IWAVZ2AB7pskGASsKZmtMbYPScMv2YuuZKKlBUkNTU1OHboOZWVdXtSCRtBvwa+ArEfFKW5MWaYs22rdujJgWEfURUV9XV9f+Ys3MrFVVCRJJPchCZFZE3JKaV6fuKtLrmtTeCOxdMPtgYGVqH1yk3czMKqgaV20JmA4sjogfFIyaA4xPw+OB2wrax0nqKWlfspPqj6Tur/WSjkrLPLtgHjMzq5DuVVjnscCngAWSnkhtFwOTgdmSzgGeB04HiIhFkmYDT5Fd8XV+RGxJ850HzAB2Ae5MP2ZmVkEVD5KI+BPFz28AnNjKPFcAVxRpbwAO7rjqzMysvXxnu5mZ5VKNrq0d2pBJt7c6bvnkkytYiZlZZfiIxMzMcnGQmJlZLg4SMzPLxUFiZma5OEjMzCwXB4mZmeXiIDEzs1wcJGZmlouDxMzMcnGQmJlZLg4SMzPLxUFiZma5OEjMzCwXB4mZmeXiIDEzs1wcJGZmlouDxMzMcnGQmJlZLg4SMzPLxUFiZma5OEjMzCwXB4mZmeXiIDEzs1wcJGZmlouDxMzMcnGQmJlZLg4SMzPLxUFiZma5OEjMzCwXB4mZmeXiIDEzs1wcJGZmlkv3ahdgWxsy6fZWxy2ffHIFKzEz27aaPyKRNFrSEknLJE2qdj1mZl1NTQeJpG7Aj4EPA8OAMyUNq25VZmZdS00HCXAEsCwinomIfwA3AWOqXJOZWZeiiKh2DdtN0lhgdER8Nr3/FHBkRFzQYrqJwMT09n3Akg4q4Z3Aix20rHKqlTqhdmqtlTqhdmqtlTqhdmrtyDr3iYi6YiNq/WS7irRtlYwRMQ2Y1uErlxoior6jl9vRaqVOqJ1aa6VOqJ1aa6VOqJ1aK1VnrXdtNQJ7F7wfDKysUi1mZl1SrQfJo8BQSftK2hkYB8ypck1mZl1KTXdtRcRmSRcAfwC6AddHxKIKltDh3WVlUit1Qu3UWit1Qu3UWit1Qu3UWpE6a/pku5mZVV+td22ZmVmVOUjMzCwXB8l2qKXHskhaLmmBpCckNVS7nmaSrpe0RtLCgrY9JN0taWl67VfNGpu1Uuulkl5I+/UJSR+pZo2ppr0l3StpsaRFkr6c2jvdfm2j1k61XyX1kvSIpD+nOi9L7Z1xn7ZWa9n3qc+RtFN6LMvTwP8ku/z4UeDMiHiqqoW1QtJyoD4iOtXNU5KOBzYAN0TEwante8BLETE5BXS/iPhGNetMdRWr9VJgQ0RcWc3aCkkaAAyIiMck9QHmA6cCE+hk+7WNWj9BJ9qvkgT0jogNknoAfwK+DJxG59unrdU6mjLvUx+RtJ8fy9IBIuIB4KUWzWOAmWl4JtkHS9W1UmunExGrIuKxNLweWAwMohPu1zZq7VQisyG97ZF+gs65T1urtewcJO03CFhR8L6RTvgLUCCAuyTNT4+K6cz2iohVkH3QAHtWuZ5tuUDSk6nrq+pdG4UkDQHeDzxMJ9+vLWqFTrZfJXWT9ASwBrg7IjrtPm2lVijzPnWQtF9Jj2XpRI6NiMPJnpB8fuqmsfymAu8BhgOrgP+sajUFJO0G/Br4SkS8Uu162lKk1k63XyNiS0QMJ3tyxhGSDq5ySa1qpday71MHSfvV1GNZImJlel0D/Iasa66zWp36zpv70NdUuZ5WRcTq9Ev7JnAdnWS/pr7xXwOzIuKW1Nwp92uxWjvrfgWIiHXAfWTnHDrlPm1WWGsl9qmDpP1q5rEsknqnE5lI6g2cBCxse66qmgOMT8PjgduqWEubmj9Eko/RCfZrOtk6HVgcET8oGNXp9mtrtXa2/SqpTtI70vAuwIeAv9A592nRWiuxT33V1nZIl89dxVuPZbmiuhUVJ2k/sqMQyB6Hc2NnqVXSL4ETyB5zvRr4NnArMBt4N/A8cHpEVP0kdyu1nkDWVRDAcuDzzX3m1SLpOOBBYAHwZmq+mOzcQ6far23UeiadaL9KOpTsZHo3sj+8Z0fE5ZL60/n2aWu1/pwy71MHiZmZ5eKuLTMzy8VBYmZmuThIzMwsFweJmZnl4iAxM7NcHCRmJZLUv+AJqn8reKLqBklTyrC+GZLGFmmvl3R1R6/PbHvV9FftmlVSRKwlux6/qk//jYgGoNN8JYCZj0jMcpJ0gqTfpeFLJc2UdJey74I5TdL3lH0nzO/TY0GQNELS/elhmn9ocfdxoQ9JelDS05JOaWV910u6T9Izkr6U2ntLul3Zd1MslHRGBXaFdVEOErOO9x7gZLJHjf8CuDciDgFeB05OYXINMDYiRgDXA609cWAI8IG0vGsl9SoyzQHAKLJnKH07LX80sDIiDkvfofL7jto4s5bctWXW8e6MiE2SFpA9rqL5Q3wBWTC8DzgYuDt75BTdyJ7KWszs9LC9pZKeIQuNlm6PiI3ARklrgL3Suq6U9B/A7yLiwY7ZNLOtOUjMOt5GgIh4U9KmeOs5RG+S/c4JWBQRR5ewrJbPMCr2TKONBcNbgO4R8bSkEcBHgO9KuisiLm/XVpiVyF1bZpW3BKiTdDRkj1OXdFAr054uaSdJ7wH2S/Nuk6SBwGsR8QvgSuDwDqjbrCgfkZhVWET8I13We7Wk3cl+D68CFhWZfAlwP1l31bkR8UbqDtuWQ4DvS3oT2ASc1xG1mxXjp/+amVku7toyM7NcHCRmZpaLg8TMzHJxkJiZWS4OEjMzy8VBYmZmuThIzMwsl/8PsREdht9YcGEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def extract_time(filename, starttime, endtime, **kwargs):\n",
    "    '''\n",
    "    Creates time bins of data and finds net counts of each time bin\n",
    "    \n",
    "        Parameters:\n",
    "            filename: data file in csv form\n",
    "            starttime (int): starting unix timestamp\n",
    "            endtime (int): ending unix timestamp\n",
    "            Keyword arguments:\n",
    "            plot_all (boolean): option to generate visual plots for each individual time bin\n",
    "        Returns:\n",
    "        tail_counts (int): net counts calculated by left-tail training\n",
    "        straight_counts(int): net counts calculated by linear estimation\n",
    "    '''\n",
    "    \n",
    "    \n",
    "    df = pd.read_csv(filename)\n",
    "    df = df.fillna(method = 'ffill') # forward-filling all the NaNs\n",
    "    timevals = df['deviceTime_unix'].to_list() # finding time data\n",
    "    df = df.set_index('deviceTime_unix')\n",
    "    lst = [[]]\n",
    "    start = timevals[0]\n",
    "    for val in timevals:\n",
    "        if val-start < 30000:\n",
    "                lst[-1].append(val)\n",
    "        else:\n",
    "            lst.append([val])\n",
    "            start = val\n",
    "    num_bins = len(lst)\n",
    "    tail_counts = []\n",
    "    straight_counts = []\n",
    "    binrange = range(0, num_bins)\n",
    "    for i in range (0, num_bins):\n",
    "        binval1 = lst[i][0]\n",
    "        last_binval = lst[i][-1] \n",
    "        print('from timestamp',binval1, 'to',last_binval)\n",
    "        spectra = df.loc[binval1:last_binval, '0':'1023'] # slicing specific time interval of spectra\n",
    "        array = spectra.to_numpy()\n",
    "        spectrum = array.sum(axis = 0)\n",
    "        x_all = np.array(range(0, 1024))\n",
    "        tail, straight = NetCounts_subtraction(x_all, spectrum, 560, 590, **params)\n",
    "        tail_counts.append(tail)\n",
    "        straight_counts.append(straight)\n",
    "    plt.bar(binrange, straight_counts, label='linear-estimate')\n",
    "    plt.title('Counts across time')\n",
    "    plt.xlabel('Time bins')\n",
    "    plt.ylabel('Counts')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    plt.bar(binrange, tail_counts, label='tail-trained')\n",
    "    plt.title('Counts across time')\n",
    "    plt.xlabel('Time bins')\n",
    "    plt.ylabel('Counts')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    return(tail_counts, straight_counts)\n",
    "    \n",
    "params = {'p0':[1000, 0.05], 'offset':None, 'tail length':0.46, 'plot_all':False}\n",
    "\n",
    "extract_time('/home/akhila/Downloads/etch_roof_CesiumIodide.csv', 1490046365,1491791174)\n",
    "x,y=setup('/home/akhila/Downloads/etch_roof_CesiumIodide.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d9e7d040",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efc6465f",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a83cf1d5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Myenv",
   "language": "python",
   "name": "myenv"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}