1b.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>MJD[days]</th>\n",
       "      <th>log10(E/GeV)</th>\n",
       "      <th>AngErr[deg]</th>\n",
       "      <th>RA[deg]</th>\n",
       "      <th>Dec[deg]</th>\n",
       "      <th>Azimuth[deg]</th>\n",
       "      <th>Zenith[deg]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>54971.15869961</td>\n",
       "      <td>3.43</td>\n",
       "      <td>1.07</td>\n",
       "      <td>128.331</td>\n",
       "      <td>53.070</td>\n",
       "      <td>256.558</td>\n",
       "      <td>143.042</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>54971.15957109</td>\n",
       "      <td>3.37</td>\n",
       "      <td>1.19</td>\n",
       "      <td>223.406</td>\n",
       "      <td>1.740</td>\n",
       "      <td>161.854</td>\n",
       "      <td>91.693</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>54971.16164457</td>\n",
       "      <td>5.37</td>\n",
       "      <td>0.53</td>\n",
       "      <td>332.873</td>\n",
       "      <td>-24.509</td>\n",
       "      <td>53.126</td>\n",
       "      <td>65.540</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>54971.16429956</td>\n",
       "      <td>5.43</td>\n",
       "      <td>0.33</td>\n",
       "      <td>243.988</td>\n",
       "      <td>-74.954</td>\n",
       "      <td>142.781</td>\n",
       "      <td>15.016</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>54971.16634297</td>\n",
       "      <td>5.33</td>\n",
       "      <td>0.55</td>\n",
       "      <td>36.228</td>\n",
       "      <td>-22.282</td>\n",
       "      <td>351.498</td>\n",
       "      <td>67.767</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107006</th>\n",
       "      <td>55347.24975340</td>\n",
       "      <td>5.82</td>\n",
       "      <td>0.86</td>\n",
       "      <td>358.582</td>\n",
       "      <td>-39.515</td>\n",
       "      <td>69.822</td>\n",
       "      <td>50.546</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107007</th>\n",
       "      <td>55347.25588194</td>\n",
       "      <td>5.62</td>\n",
       "      <td>0.61</td>\n",
       "      <td>110.336</td>\n",
       "      <td>-40.464</td>\n",
       "      <td>320.340</td>\n",
       "      <td>49.518</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107008</th>\n",
       "      <td>55347.25789030</td>\n",
       "      <td>3.39</td>\n",
       "      <td>1.17</td>\n",
       "      <td>150.619</td>\n",
       "      <td>27.546</td>\n",
       "      <td>280.710</td>\n",
       "      <td>117.500</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107009</th>\n",
       "      <td>55347.26411744</td>\n",
       "      <td>3.39</td>\n",
       "      <td>0.50</td>\n",
       "      <td>213.782</td>\n",
       "      <td>-1.188</td>\n",
       "      <td>219.806</td>\n",
       "      <td>88.760</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>107010</th>\n",
       "      <td>55347.27512234</td>\n",
       "      <td>3.24</td>\n",
       "      <td>1.77</td>\n",
       "      <td>31.417</td>\n",
       "      <td>16.139</td>\n",
       "      <td>46.147</td>\n",
       "      <td>106.190</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>107011 rows × 7 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "             MJD[days] log10(E/GeV) AngErr[deg]  RA[deg] Dec[deg]  \\\n",
       "0       54971.15869961         3.43        1.07  128.331   53.070   \n",
       "1       54971.15957109         3.37        1.19  223.406    1.740   \n",
       "2       54971.16164457         5.37        0.53  332.873  -24.509   \n",
       "3       54971.16429956         5.43        0.33  243.988  -74.954   \n",
       "4       54971.16634297         5.33        0.55   36.228  -22.282   \n",
       "...                ...          ...         ...      ...      ...   \n",
       "107006  55347.24975340         5.82        0.86  358.582  -39.515   \n",
       "107007  55347.25588194         5.62        0.61  110.336  -40.464   \n",
       "107008  55347.25789030         3.39        1.17  150.619   27.546   \n",
       "107009  55347.26411744         3.39        0.50  213.782   -1.188   \n",
       "107010  55347.27512234         3.24        1.77   31.417   16.139   \n",
       "\n",
       "       Azimuth[deg] Zenith[deg]  \n",
       "0           256.558     143.042  \n",
       "1           161.854      91.693  \n",
       "2            53.126      65.540  \n",
       "3           142.781      15.016  \n",
       "4           351.498      67.767  \n",
       "...             ...         ...  \n",
       "107006       69.822      50.546  \n",
       "107007      320.340      49.518  \n",
       "107008      280.710     117.500  \n",
       "107009      219.806      88.760  \n",
       "107010       46.147     106.190  \n",
       "\n",
       "[107011 rows x 7 columns]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "with open(\"icecube_10year_ps/events/IC59_exp.csv\", 'r') as f1:\n",
    "    lines = f1.readlines()\n",
    "content=[]\n",
    "column=lines[0].split()\n",
    "column.pop(0)\n",
    "for line in lines[1:]:\n",
    "    \n",
    "    content.append(line.split())\n",
    "\n",
    "event = pd.DataFrame(content, columns = column)\n",
    "event"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>log10(E_nu/GeV)_min</th>\n",
       "      <th>log10(E_nu/GeV)_max</th>\n",
       "      <th>Dec_nu_min[deg]</th>\n",
       "      <th>Dec_nu_max[deg]</th>\n",
       "      <th>A_Eff[cm^2]</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.20</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-73.74</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.20</td>\n",
       "      <td>2.40</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-73.74</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.40</td>\n",
       "      <td>2.60</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-73.74</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.60</td>\n",
       "      <td>2.80</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-73.74</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.80</td>\n",
       "      <td>3.00</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-73.74</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1995</th>\n",
       "      <td>9.00</td>\n",
       "      <td>9.20</td>\n",
       "      <td>73.74</td>\n",
       "      <td>90.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1996</th>\n",
       "      <td>9.20</td>\n",
       "      <td>9.40</td>\n",
       "      <td>73.74</td>\n",
       "      <td>90.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1997</th>\n",
       "      <td>9.40</td>\n",
       "      <td>9.60</td>\n",
       "      <td>73.74</td>\n",
       "      <td>90.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1998</th>\n",
       "      <td>9.60</td>\n",
       "      <td>9.80</td>\n",
       "      <td>73.74</td>\n",
       "      <td>90.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1999</th>\n",
       "      <td>9.80</td>\n",
       "      <td>10.00</td>\n",
       "      <td>73.74</td>\n",
       "      <td>90.00</td>\n",
       "      <td>0.00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>2000 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     log10(E_nu/GeV)_min log10(E_nu/GeV)_max Dec_nu_min[deg] Dec_nu_max[deg]  \\\n",
       "0                   2.00                2.20          -90.00          -73.74   \n",
       "1                   2.20                2.40          -90.00          -73.74   \n",
       "2                   2.40                2.60          -90.00          -73.74   \n",
       "3                   2.60                2.80          -90.00          -73.74   \n",
       "4                   2.80                3.00          -90.00          -73.74   \n",
       "...                  ...                 ...             ...             ...   \n",
       "1995                9.00                9.20           73.74           90.00   \n",
       "1996                9.20                9.40           73.74           90.00   \n",
       "1997                9.40                9.60           73.74           90.00   \n",
       "1998                9.60                9.80           73.74           90.00   \n",
       "1999                9.80               10.00           73.74           90.00   \n",
       "\n",
       "     A_Eff[cm^2]  \n",
       "0           0.00  \n",
       "1           0.00  \n",
       "2           0.00  \n",
       "3           0.00  \n",
       "4           0.00  \n",
       "...          ...  \n",
       "1995        0.00  \n",
       "1996        0.00  \n",
       "1997        0.00  \n",
       "1998        0.00  \n",
       "1999        0.00  \n",
       "\n",
       "[2000 rows x 5 columns]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "with open(\"icecube_10year_ps/irfs/IC59_effectiveArea.csv\", 'r') as f2:\n",
    "    lines = f2.readlines()\n",
    "content=[]\n",
    "column=lines[0].split()\n",
    "column.pop(0)\n",
    "for line in lines[1:]:\n",
    "    \n",
    "    content.append(line.split())\n",
    "\n",
    "effa = pd.DataFrame(content, columns = column)\n",
    "effa"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>log10(E_nu/GeV)_min</th>\n",
       "      <th>log10(E_nu/GeV)_max</th>\n",
       "      <th>Dec_nu_min[deg]</th>\n",
       "      <th>Dec_nu_max[deg]</th>\n",
       "      <th>log10(E/GeV)_min</th>\n",
       "      <th>log10(E/GeV)_max</th>\n",
       "      <th>PSF_min[deg]</th>\n",
       "      <th>PSF_max[deg]</th>\n",
       "      <th>AngErr_min[deg]</th>\n",
       "      <th>AngErr_max[deg]</th>\n",
       "      <th>Fractional_Counts</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-10.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-10.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-10.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-10.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2.00</td>\n",
       "      <td>2.50</td>\n",
       "      <td>-90.00</td>\n",
       "      <td>-10.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.000e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.00e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369595</th>\n",
       "      <td>8.50</td>\n",
       "      <td>9.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>90.000</td>\n",
       "      <td>7.597</td>\n",
       "      <td>7.830</td>\n",
       "      <td>6.935e+01</td>\n",
       "      <td>1.095e+02</td>\n",
       "      <td>1.81e+00</td>\n",
       "      <td>2.10e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369596</th>\n",
       "      <td>8.50</td>\n",
       "      <td>9.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>90.000</td>\n",
       "      <td>7.597</td>\n",
       "      <td>7.830</td>\n",
       "      <td>6.935e+01</td>\n",
       "      <td>1.095e+02</td>\n",
       "      <td>2.10e+00</td>\n",
       "      <td>2.43e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369597</th>\n",
       "      <td>8.50</td>\n",
       "      <td>9.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>90.000</td>\n",
       "      <td>7.597</td>\n",
       "      <td>7.830</td>\n",
       "      <td>6.935e+01</td>\n",
       "      <td>1.095e+02</td>\n",
       "      <td>2.43e+00</td>\n",
       "      <td>2.81e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369598</th>\n",
       "      <td>8.50</td>\n",
       "      <td>9.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>90.000</td>\n",
       "      <td>7.597</td>\n",
       "      <td>7.830</td>\n",
       "      <td>6.935e+01</td>\n",
       "      <td>1.095e+02</td>\n",
       "      <td>2.81e+00</td>\n",
       "      <td>3.25e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>369599</th>\n",
       "      <td>8.50</td>\n",
       "      <td>9.00</td>\n",
       "      <td>10.00</td>\n",
       "      <td>90.000</td>\n",
       "      <td>7.597</td>\n",
       "      <td>7.830</td>\n",
       "      <td>6.935e+01</td>\n",
       "      <td>1.095e+02</td>\n",
       "      <td>3.25e+00</td>\n",
       "      <td>3.77e+00</td>\n",
       "      <td>0.0000e+00</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>369600 rows × 11 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "       log10(E_nu/GeV)_min log10(E_nu/GeV)_max Dec_nu_min[deg]  \\\n",
       "0                     2.00                2.50          -90.00   \n",
       "1                     2.00                2.50          -90.00   \n",
       "2                     2.00                2.50          -90.00   \n",
       "3                     2.00                2.50          -90.00   \n",
       "4                     2.00                2.50          -90.00   \n",
       "...                    ...                 ...             ...   \n",
       "369595                8.50                9.00           10.00   \n",
       "369596                8.50                9.00           10.00   \n",
       "369597                8.50                9.00           10.00   \n",
       "369598                8.50                9.00           10.00   \n",
       "369599                8.50                9.00           10.00   \n",
       "\n",
       "       Dec_nu_max[deg] log10(E/GeV)_min log10(E/GeV)_max PSF_min[deg]  \\\n",
       "0              -10.000            0.000            0.000    0.000e+00   \n",
       "1              -10.000            0.000            0.000    0.000e+00   \n",
       "2              -10.000            0.000            0.000    0.000e+00   \n",
       "3              -10.000            0.000            0.000    0.000e+00   \n",
       "4              -10.000            0.000            0.000    0.000e+00   \n",
       "...                ...              ...              ...          ...   \n",
       "369595          90.000            7.597            7.830    6.935e+01   \n",
       "369596          90.000            7.597            7.830    6.935e+01   \n",
       "369597          90.000            7.597            7.830    6.935e+01   \n",
       "369598          90.000            7.597            7.830    6.935e+01   \n",
       "369599          90.000            7.597            7.830    6.935e+01   \n",
       "\n",
       "       PSF_max[deg] AngErr_min[deg] AngErr_max[deg] Fractional_Counts  \n",
       "0         0.000e+00        0.00e+00        0.00e+00        0.0000e+00  \n",
       "1         0.000e+00        0.00e+00        0.00e+00        0.0000e+00  \n",
       "2         0.000e+00        0.00e+00        0.00e+00        0.0000e+00  \n",
       "3         0.000e+00        0.00e+00        0.00e+00        0.0000e+00  \n",
       "4         0.000e+00        0.00e+00        0.00e+00        0.0000e+00  \n",
       "...             ...             ...             ...               ...  \n",
       "369595    1.095e+02        1.81e+00        2.10e+00        0.0000e+00  \n",
       "369596    1.095e+02        2.10e+00        2.43e+00        0.0000e+00  \n",
       "369597    1.095e+02        2.43e+00        2.81e+00        0.0000e+00  \n",
       "369598    1.095e+02        2.81e+00        3.25e+00        0.0000e+00  \n",
       "369599    1.095e+02        3.25e+00        3.77e+00        0.0000e+00  \n",
       "\n",
       "[369600 rows x 11 columns]"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "with open(\"icecube_10year_ps/irfs/IC59_smearing.csv\", 'r') as f2:\n",
    "    lines = f2.readlines()\n",
    "content=[]\n",
    "column=lines[0].split()\n",
    "column.pop(0)\n",
    "for line in lines[1:]:\n",
    "    \n",
    "    content.append(line.split())\n",
    "\n",
    "smear = pd.DataFrame(content, columns = column)\n",
    "smear"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD4CAYAAAAO9oqkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAmu0lEQVR4nO3deXhV1dn+8e/DPMhMmANhCIOACAREsRTFIqAVtVqxVpHSolaptX21aAdbq6/11V8d2lpFRq0KqKhUUcQBBS2jIPMQGRMCCQYIEDI/vz+yoVEZQqadk9yf68qVc9bZe58ncCX32WuvtZe5OyIiUrlVCbsAEREJn8JAREQUBiIiojAQEREUBiIiAlQLu4Ciatq0qcfExIRdhohIRFmxYsU+d4/6ZnvEhkFMTAzLly8PuwwRkYhiZjtO1K5uIhERURiIiIjCQEREUBiIiAgKAxERQWEgIiIoDEREBIWBiEjEeG1FAoP+7yOyc/NK/NgRO+lMRKQy2bgnjV+/8gUAR7NzqV61ZD/Ln/ZoZjbFzJLNbG2BtkfNbKOZrTaz182sYYHX7jWzeDPbZGaXFmgfFrTFm9mEAu3tzWxJ0D7TzGqU4M8nIhLxViccYNTExcefl8aaZIWJlmnAsG+0zQd6uPs5wGbgXgAzOxsYBXQP9nnazKqaWVXgH8Bw4Gzg+mBbgEeAx929E7AfGFusn0hEpAJZtj2VHz23hLNqlm5HzmnDwN0/AVK/0faeu+cETxcDbYLHI4EZ7p7p7tuAeKB/8BXv7lvdPQuYAYw0MwMuBl4N9p8OXFm8H0lEpGJYuCWFGycvoVn9mrxy6/nUq1V6gVASnU4/Ad4JHrcGdhV4LSFoO1l7E+BAgWA51n5CZjbOzJab2fKUlJQSKF1EpHyav34vY6ctJ6ZJXWaOO5+WDWqX6vsVKwzM7LdADvBiyZRzau4+0d3j3D0uKupbd2AVEakQ5nyxm1v/tYJureozY9wAourVLPX3LPI5h5ndDFwODHE/fjkjEYgusFmboI2TtH8FNDSzasHZQcHtRUQqnVnLdvGb2avpF9OYKTf3K/VrBccU6czAzIYB9wBXuHt6gZfmAKPMrKaZtQdigaXAMiA2GDlUg/yLzHOCEPkIuCbYfzTwZtF+FBGRyDb1023c89pqvhMbxfQx/cssCKBwQ0tfBv4DdDGzBDMbC/wdqAfMN7NVZvYMgLuvA2YB64F3gdvdPTf41H8HMA/YAMwKtgX4DfArM4sn/xrC5BL9CUVEIsA/PornT/9ez6Xdm/PcTX2pXaNqmb7/aWPH3a8/QfNJ/2C7+0PAQydonwvMPUH7VvJHG4mIVDruzqPzNvH0gi+58txWPHZtL6qV8ISywtAMZBGRkOTlOQ+8tZ5pn23n+v7RPHhlT6pWsVBqURiIiIQgN8+5b/YaZi7fxU8Gtuf3l3cjf+pVOBQGIiJlLDs3j1/N+oJ/f7GbX1zcibu+1znUIACFgYhImcrIzuWOl1by/oa9TBjelVu/2zHskgCFgYhImUnPyuGWF1awcMs+HhjZnZvOjwm7pOMUBiIiZSAtI5ux05axYsd+Hr3mHK6Niz79TmVIYSAiUsr2H8li9NSlrN+dxlPX9+byc1qFXdK3KAxEREpR8qEMbpy0lG1fHeHZG/sypFvzsEs6IYWBiEgpSTxwlB9PWsLetAym3tyPgZ2ahl3SSSkMRERKwfZ9R7hh0hLSjmbzwtj+9G3XOOySTklhICJSwrbsPcQNk5aQnZvHy+MG0KN1g7BLOq2yvwGGiEgF9mn8Pq56+jMcmHnL+RERBKAzAxGREvPZl/u4YdISAN4afyExTeuGXFHh6cxARKQEvLV6NzdPWXb8eSQFASgMRESKbcqibYx/eSW9oiOjS+hEFAYiIkWUl+c8/M4GHnhrPUPPbs4LY88Lu6Qi0zUDEZEiyM7N4zevrmb2ykR+PKAtf7qiR2hrEZQEhYGIyBk6nJnDbf/Kv+Hc/wztzO0XdQr9FtTFpTAQETkDKYcy+cm0ZaxPSuORH/Tkun5twy6pRCgMREQKafu+I4yeupS9aRk8d1NfLu5aPu8zVBQKAxGRQlidcIAxU5eR585LPxtAn7aNwi6pRCkMRERO4+PNKdz2rxU0qlOD58f2p2PUWWGXVOIUBiIipzD78wTueXU1sc3rMX1MP5rVrxV2SaVCYSAicgLuzrOfbOUv72zkgo5NeObGvtSvVT3sskrNaSedmdkUM0s2s7UF2hqb2Xwz2xJ8bxS0m5k9ZWbxZrbazPoU2Gd0sP0WMxtdoL2vma0J9nnKIn18lohEvLw854G31vOXdzby/V6tmDqmX4UOAijcDORpwLBvtE0APnD3WOCD4DnAcCA2+BoH/BPywwO4HzgP6A/cfyxAgm1+VmC/b76XiEiZyczJZfyMlUz9dDs/GdieJ687l5rVqoZdVqk7bRi4+ydA6jeaRwLTg8fTgSsLtD/v+RYDDc2sJXApMN/dU919PzAfGBa8Vt/dF7u7A88XOJaISJlKy8hm9JSlvL06iftGdOX3l3ejSgTPKj4TRb1m0Nzdk4LHe4Bjg21bA7sKbJcQtJ2qPeEE7SdkZuPIP+OgbduKMdFDRMqHvWkZjJ6ylPjkwzxx3blc2fukf4oqpGLfqC74RO8lUEth3muiu8e5e1xUVFRZvKWIVALxyYe5+unP2JWazpSb+1W6IICih8HeoIuH4Hty0J4IRBfYrk3Qdqr2NidoFxEpEyt27OeaZz4jMyeXmbecz6DOlfODZlHDYA5wbETQaODNAu03BaOKBgAHg+6kecBQM2sUXDgeCswLXkszswHBKKKbChxLRKRUvb9+LzdMWkzD2tV57bYLImaJytJw2msGZvYyMBhoamYJ5I8K+gswy8zGAjuAHwabzwVGAPFAOjAGwN1TzezPwLFlgB5w92MXpX9O/oil2sA7wZeISKmasXQn972+hh6tGzDl5n40Patm2CWF6rRh4O7Xn+SlISfY1oHbT3KcKcCUE7QvB3qcrg4RkZLg7vzhzXW8sHgH3+0cxdM39KFuTc2/1b+AiFQamTm5DH9yIVtTjnBudEMmjY6jelUt+Aha9lJEKonUI1ncOGkpW1OOAPCDPq0VBAXozEBEKrz45EP8ZNpy9qRl0LhuDVKPZIVdUrmjWBSRCu2TzSlc9fRnpGflMGPcAEb0bBF2SeWSwkBEKqzn/7OdMdOW0bphbd64fWCFW5CmJKmbSEQqnJzcPP781nqm/2cHQ7o248nre3OWRgydkv51RKRCScvI5o6XVvLJ5hR+9p32TBjejaqV5GZzxaEwEJEKY+dX6fxk+jK27zvCX67uyaj+uqFlYSkMRKRCWLotlVteWE6ew/Nj+3NBx6ZhlxRRFAYiEvFeXZHAvbNXE92oDpNv7kf7pnXDLiniKAxEJGLl5TmPvreJfy74kgs6NuGfN/SlQZ2KvTxlaVEYiEhESs/K4a6Zq5i3bi8/Oq8tf7qiu2YUF4PCQEQiTtLBo/x0+nI2JKXxh8vPZszAGPLvgi9FpTAQkYiyOuEAP52+nPSsXCaP7sdFXZuFXVKFoDAQkYgxd00Sv5q1iiZ1a/LabefRpUW9sEuqMBQGIlLuuTv/+Ciex97bTN92jXj2xr6VfjGakqYwEJFyLSM7lwmvreaNVbu5qndrHr66J7WqVw27rApHYSAi5da+w5nc8sIKVuzYz92XduHngzvqQnEpURiISLm0ac8hxk5fxr7DmTx9Qx9G9GwZdkkVmsJARMqdJ9/fwuPvb6Z5/ZrMuuV8zmnTMOySKjyFgYiUG3l5zlMfbuGJ97cA8ObtF9KiQa2Qq6ocFAYiUi4cTM/mrlmr+HBj8vE2BUHZ0dxtEQndhqQ0vv/3RSzcksKfR3YPu5xKSWEgIqF6Y2UiVz39KZk5ucwYdz43nh8TdkmVUrHCwMzuMrN1ZrbWzF42s1pm1t7MlphZvJnNNLMawbY1g+fxwesxBY5zb9C+ycwuLebPJCIRICsnjz/OWccvZ67inDYN+ff4C+nbTmsUh6XIYWBmrYFfAHHu3gOoCowCHgEed/dOwH5gbLDLWGB/0P54sB1mdnawX3dgGPC0mWlGiUgFtjctgx89t5hpn21n7IXtefGn59Gsnq4PhKm43UTVgNpmVg2oAyQBFwOvBq9PB64MHo8MnhO8PsTyZ4+MBGa4e6a7bwPigf7FrEtEyqml21K5/G+LWLc7jaeu783vLz9bt54uB4r8P+DuicBjwE7yQ+AgsAI44O45wWYJQOvgcWtgV7BvTrB9k4LtJ9jna8xsnJktN7PlKSkpRS1dRELg7kxZtI0fPbeYs2pW4807BnJFr1ZhlyWB4nQTNSL/U317oBVQl/xunlLj7hPdPc7d46KiokrzrUSkBKVn5fDLmat44K31XNS1GW/eMZDOzXXH0fKkOPMMLgG2uXsKgJnNBgYCDc2sWvDpvw2QGGyfCEQDCUG3UgPgqwLtxxTcR0Qi3PZ9R7j1XyvYtPcQd1/ahdu+25EqVXR/ofKmOB11O4EBZlYn6PsfAqwHPgKuCbYZDbwZPJ4TPCd4/UN396B9VDDaqD0QCywtRl0iUk68v34v3//7IvakZTBtTH9uv6iTgqCcKvKZgbsvMbNXgc+BHGAlMBF4G5hhZg8GbZODXSYDL5hZPJBK/ggi3H2dmc0iP0hygNvdPbeodYlI+HLznCfe38zfPoyne6v6PPPjvkQ3rhN2WXIKxbodhbvfD9z/jeatnGA0kLtnANee5DgPAQ8VpxYRKR8OpGdx54xVfLw5hWv6tuHBK3to/YEIoHsTiUiJWZt4kNteXMGegxk8eGUPbjivrdYfiBAKAxEpEa+tSOC+19fQqE4NZt1yPr3bajZxJFEYiEixZOXk8ee31vPC4h0M6NCYv/+oj9YnjkAKAxEpsj0HM7jtxRWs3HmAcYM6cM+lXaim2cQRSWEgIkXy0cZk7n71C9KzcvnHj/pw2TlaljKSKQxE5Izk5jn/9+5Gnv1kKwDz7xpErGYTRzyFgYgU2p6DGdw1cxX/2foVAHHtGikIKgh17olIoXywYS/Dn/yEVbsOcF1c9Ol3kIiiMwMROaXMnFwenruRaZ9tp1vL+vzt+t7sT89i5vJdp99ZIobCQERO6suUw4x/aSXrk9IYMzCG3wzrSq3qVVm2PTXs0qSEKQxE5FvcnVdWJHD/m+uoVb0Kk0fHMaRb87DLklKkMBCRr0nLyOZ3r69lzhe7GdChMU9c15sWDbQkZUWnMBCR41btOsD4lz9n94EM/mdoZ24b3ImquuV0paAwEBHy8pyJC7fy2LxNNK9fi1m3DKBvu8ZhlyVlSGEgUsklH8rgVzO/YFH8Pkb0bMHDV59Dg9rVwy5LypjCQKQSW7ApmV/P+oIjWTk8fHVPRvWL1i2nKymFgUgllJWTx6PzNvLcwm10bVGPGdcP0EziSk5hIFLJbN93hPEvr2RN4kFuHNCO317WTSuRicJApDJ5fWUCv3t9LdWqVuHZG/tyafcWYZck5YTCQKQSOJyZwx/eWMvslYn0b9+YJ647l1YNa4ddlpQjCgORCm5NwkHGv/w5O1PT+eUlsYy/OFZzB+RbFAYiFVRenjPl02088u5Gmp5Vkxnjzqd/e80dkBNTGIhUQAn707nwkY8AuLR7cx75wTk0rFMj5KqkPFMYiFQg7s7MZbt48O0Nx9ue+XFfzR2Q0yrW4jZm1tDMXjWzjWa2wczON7PGZjbfzLYE3xsF25qZPWVm8Wa22sz6FDjO6GD7LWY2urg/lEhllHTwKDdPXcaE2Wvo2brB8XYFgRRGcVc6exJ41927Ar2ADcAE4AN3jwU+CJ4DDAdig69xwD8BzKwxcD9wHtAfuP9YgIjI6bk7r65IYOjjn7B0Wyp/uqI7L/70vLDLkghT5G4iM2sADAJuBnD3LCDLzEYCg4PNpgMLgN8AI4Hn3d2BxcFZRctg2/nunhocdz4wDHi5qLWJVBbJaRncO3sNH2xMpl9MIx69phcxTeuGXZZEoOJcM2gPpABTzawXsAK4E2ju7knBNnuAYytitAYKrpOXELSdrP1bzGwc+WcVtG3bthili0Q2d2fOF7v5w5vryMjO5feXn82YC2KooiGjUkTFCYNqQB9gvLsvMbMn+W+XEADu7mbmxSnwG8ebCEwEiIuLK7HjikSSlEOZ/O6NNcxbt5c+bRvy6LW96Bh1VthlSYQrThgkAAnuviR4/ir5YbDXzFq6e1LQDZQcvJ4IRBfYv03Qlsh/u5WOtS8oRl0iFdZbq3fz+zfWciQrl/tGdGXshR00gUxKRJEvILv7HmCXmXUJmoYA64E5wLERQaOBN4PHc4CbglFFA4CDQXfSPGComTUKLhwPDdpEJJB6JIvbX/qcO15aSdvGdXh7/IWMG9RRQSAlprjzDMYDL5pZDWArMIb8gJllZmOBHcAPg23nAiOAeCA92BZ3TzWzPwPLgu0eOHYxWUTg3bV7+N0bazh4NJu7L+3CLYM6UK1qcQcCinxdscLA3VcBcSd4acgJtnXg9pMcZwowpTi1iFQ0B9Kz+OOcdbyxajfdW9XnXz89j64t6oddllRQmoEsUg59sGEvE2avYf+RLO66pDM/v6gj1XU2IKVIYSBSjhw8ms0D/17Pa58n0LVFPaaN6Uf3Vg1Ov6NIMSkMRMqJBZuSmfDaGlIOZzL+4k6MvziWGtV0NiBlQ2EgErJDGdk89PYGZizbRWyzs5h4U1/OadMw7LKkklEYiITopSU7eeTdjRzKyOa2wR25c0is1iOWUCgMREKQciiT/527gddXJgIw++cX0Ket7s8o4VEYiJShvDzn5WU7eeSdjRzNzj3eriCQsOnqlEgZWb87jR888xm/fX0t3Vs14J07B4VdkshxOjMQKWWHM3N4Yv5mpn62nUZ1qvPEdecy8txWWnRGyhWFgUgpcXfmrdvLn/69jj1pGfyof1vuubQrDepUD7s0kW9RGIiUgl2p6dw/Zx0fbkymW8v6/OOGProuIOWawkCkBGXl5DFp0Vae+mALVcz43WXduPmCGN1YTso9hYFICVmy9St+98ZatiQfZlj3Fvzh+2fTqmHtsMsSKRSFgUgxpR7J4uG5G3hlRQKtG9Zm8ug4hnRrfvodRcoRhYFIEeXlOa+s2MXD72zkcEYOtw3uyC8ujqV2Dc0glsijMBApgk17DvHb19ewfMd++sc05sGretC5eb2wyxIpMoWByBlIz8rhyQ+2MHnhNurVqsaj15zDNX3baM6ARDyFgUghzV+/lz/OWUfigaNcFxfNhOFdaVS3RthliZQIhYHIaSQeOMof56xj/vq9dG5+Fq/cej79YhqHXZZIiVIYiJzE0axc7nt9zfE7i04Y3pWxF7bX8pNSISkMRL4hL89584tE/u/dTSQdzADgrz/sxdV92oRcmUjpURiIFLBk61c8NHcDqxMO0rN1A45k5pCWkUOjOro2IBWbzndFgG37jnDLC8u5buJiUg5l8vh1vXjz9oH0aaf7CUnlUOwzAzOrCiwHEt39cjNrD8wAmgArgBvdPcvMagLPA32Br4Dr3H17cIx7gbFALvALd59X3LpECuNAehZPfRDPC4u3U71qFf5naGfGXthBE8ek0imJbqI7gQ1A/eD5I8Dj7j7DzJ4h/4/8P4Pv+929k5mNCra7zszOBkYB3YFWwPtm1tndc7/5RiIlJSsnjxcW7+CpD7ZwKCOb6/pFc9f3OtOsXq2wSxMJRbG6icysDXAZMCl4bsDFwKvBJtOBK4PHI4PnBK8PCbYfCcxw90x33wbEA/2LU5fIybg7767dw9DHP+bPb63nnDYNmHvnd3j46nMUBFKpFffM4AngHuDYPPwmwAF3zwmeJwCtg8etgV0A7p5jZgeD7VsDiwscs+A+IiVmdcIBHnx7A0u3pRLb7CymjunH4M5Rmj0sQjHCwMwuB5LdfYWZDS6xik79nuOAcQBt27Yti7eUCmD3gaM8Nm8Ts1cm0qRuDR68sgej+kVrjQGRAopzZjAQuMLMRgC1yL9m8CTQ0MyqBWcHbYDEYPtEIBpIMLNqQAPyLyQfaz+m4D5f4+4TgYkAcXFxXozapRI4nJnDsx9/ycRPtuLAbYM78vPBHalXS8tOinxTkcPA3e8F7gUIzgz+x91vMLNXgGvIH1E0Gngz2GVO8Pw/wesfurub2RzgJTP7K/kXkGOBpUWtSyQ3z5m1fBf/773N7DucyRW9WnHPsC60aVQn7NJEyq3SmHT2G2CGmT0IrAQmB+2TgRfMLB5IJX8EEe6+zsxmAeuBHOB2jSSSolq4JYWH3t7Axj2H6NuuEc/d1JfeWntY5LRKJAzcfQGwIHi8lROMBnL3DODak+z/EPBQSdQildOWvYd4aO4GFmxKIbpxbZ6+oQ/De7TQxWGRQtLtKCSiJR44yuBHPyI716lXqxr3jejK6AtiqFlNk8ZEzoTCQCJS0sGj/OOjeGYu20V2bv5Ygo/vvojGWl9ApEgUBhJR9hzM4OkF8cxYugvHuTYumpeW7ARQEIgUg8JAIkJyWgZPL/iSl5buJC/PuTauDbdf1Ik2jeocDwMRKTqFgZRryYcyeGbBVl5csoOcPOcHfVoz/uJYohtrmKhISVIYSLmUciiTZz/+kn8t2UF2rnNV79aMv7gT7ZrUDbs0kQpJYSDlyleHM5n4yVae/88OMnNyubJ3a35xcSwxTRUCIqVJYSDlQuqRrCAEtpORncsVvVoxfkgsHaPOCrs0kUpBYSChOpCexXMLtzLt0+2kZ+fy/XNa8YshsXRqphAQKUsKAwnFwfRsJi3aytRPt3MkK4fLerbkziGxxDavd/qdRaTEKQykTB08ms3kRduYumgbhzJzGNGzBXcO6UyXFgoBkdN5YGR37pr5BfVrlfyfboWBlIm0jGymLtrOpEVbOZSRw7DuLbjzkli6tax/+p1FBICrerfhqt5tSuXYCgMpVXvTMvjre5uZuXwXAEPPbs6dl8TSvVWDkCsTkYIUBlIq1u9OY9Kirfz7i93H7x00qHMUE2+KC7kyETkRhYGUGHdnweYUJi/cxqL4fdSpUZUbzmtHwzrVeeL9LTSuoxXGRMorhYEUW0Z2Lm+uSmTSwm1sST5M8/o1uWdYF27o344Gdarz+sqEsEsUkdNQGEiRpR7J4oX/7OCFxdvZdziLri3q8f+u7cX3e7WiRjUtNi+V08s/G0CrhrXCLuOMKQzkjH2ZcpjJi7bx2ooEMnPyGNwlip99pwMXdGyilcWk0ju/Y5OwSygShYEUiruzZFsqkxZu5f0NydSoVoWre7dm7IXtNVFMpAJQGMgpZefmMXdNEpMWbmNN4kEa163BL4bEcuOAdkTVqxl2eSJSQhQGckJpGdnMWLqTaZ9uZ/fBDDo0rctDV/XgB33aUKu61heWyFWtSv71rKpVdF2rIIWBfM2u1HSmfrqdmct2ciQrlwEdGvPAyB5c3LUZVaroeoBEvrsu6cz89Xu5Nq50ZvJGKoWB4O68sSqRx+ZtJungUcyMy89pyU8v7EDPNpopLBVLgzrV+XTCxWGXUe4oDCqx1CNZzP48gRnLdhGffBiAcYM6cPMFMbRqWDvk6kSkLBU5DMwsGngeaA44MNHdnzSzxsBMIAbYDvzQ3fdb/pjDJ4ERQDpws7t/HhxrNPC74NAPuvv0otYlp5aX53z25Ve8vGwn763bQ3au06dtw+Ov3zeiW3jFiUhoinNmkAP82t0/N7N6wAozmw/cDHzg7n8xswnABOA3wHAgNvg6D/gncF4QHvcDceSHygozm+Pu+4tRm3zDnoMZvLpiFzOX72JX6lEa1qnOjQNiuK5fNF1a1CNmwtthlygiISpyGLh7EpAUPD5kZhuA1sBIYHCw2XRgAflhMBJ43t0dWGxmDc2sZbDtfHdPBQgCZRjwclFrk3w5uXl8tCmFmct28uHGZPIcLujYhLsv7crQs5trVJCUW0+OOpe0o9lhl1GplMg1AzOLAXoDS4DmQVAA7CG/Gwnyg2JXgd0SgraTtZ/ofcYB4wDatm1bEqVXSDu/Smfm8p28sjyB5EOZRNWrya3f7ch1/aJp10QLy0v5N/LcE/4JkFJU7DAws7OA14BfuntawdsRuLubmRf3PQocbyIwESAuLq7EjlsRZObkMm/dXmYu28mn8V9RxeCiLs0Y1b8tF3WJolpVjakWkZMrVhiYWXXyg+BFd58dNO81s5bunhR0AyUH7YlAdIHd2wRtify3W+lY+4Li1FWZbN57iBlLdzF7ZQIH0rNp06g2v/5eZ66Ni6ZFg8i7WZaIhKM4o4kMmAxscPe/FnhpDjAa+Evw/c0C7XeY2QzyLyAfDAJjHvC/ZtYo2G4ocG9R66oM0rNyeGt1EjOW7uTznQeoXtUY2r0Fo/pFM7BjU00Ok1IX3agOELk3ZZNvK86ZwUDgRmCNma0K2u4jPwRmmdlYYAfww+C1ueQPK40nf2jpGAB3TzWzPwPLgu0eOHYxWf7L3Zm3bg8fbEjmnbV7OJyZQ8eouvzusm5c1bs1Tc7SfYKk7LRoUIvtf7ks7DKkBBVnNNEi4GQfQYecYHsHbj/JsaYAU4paS0Xl7qzbncbba5J4a/VudqUeBeCavm0Y1S+avu0a6ZbRIlIiNAO5nHF3NiQd4u01u3l7dRLbv0qnahVjQIfG7Eo9ShWDx67tFXaZIlLBKAzKAXdn095DvL06ibdXJ7F13xGqVjEu6NiEW7/bkaHdW1C/VjU6/fYdnQlIodWoWoWs3Lywy5AIoTAI0Za9h3hrdRJvr0kiPvkwVQwGdGjCT7/TgUu7N//adYAc/VLLGfro7sGs3KmJ/FI4CoMyFp98OP8MYM1uNu89jBmc174xoy/owbDuLbRgjJSY1g1r01o3HJRCUhiUga0pxwIgiY17DmEG/WIa88DI7gzr0YJm9TQfQETCpTAoJdv3HeHtNfnXANYnpQEQ164R93//bEb0bEnz+goAyadpIVIeKAxKiLuzOuEgf/twC3vSMlibmB8Afdo25PeXn82Ini1o2UCn7PJ1r9x6PjG6X5SUAwqDYkjLyOaz+H18vDmFjzelsPtgBgAdmuZPBhves6X6bOWU+sU0DrsEEUBhcEby8vIngX28OZmPN6fw+c4D5OY59WpWY2Cnpuw+uAeAZ2/sS2zzeiFXKyJSeAqD09h3OJOFW/I/+S/cso+vjmQB0LN1A277bkcGdY6id9uGVK9ahUv++vHx5SOlYvhhXDQLNqVwjtaClgpOYfAN2bl5rNx54Pin/2N9/03q1mBQ5ygGdW7Kd2KjaKp7AVUKI3q21D14pFJQGAAJ+9P5ZPM+Pt6czGfxX3EoM4eqVYy+bRtx96VdGBQbRfdW9XU3UBGpsCplGGRk57JkWyofb0rhky0px7t2WjeszeW9WvHdzk25oFNT6teqHnKlIiJlo9KFQcGF32tUq8KADk0Y1S+awV2i6Bh1lu79E4FW/eF7+n8TKaZKFwbHTBvTj/PaN6F2DS0KH+ka1qkRdgkiEa/ShsHgLs3CLqHSiG2WP8z23OiG4RYiIidVacNAyk6P1g1Y88eh1NM1GJFyq0rYBUjloCAQKd8UBgLAuEEdqKPrJyKVlrqJBID7RnTjvhHdwi5DREKiM4MIUSUYOjm4c1TIlYhIRaQzgwhRpYrptggiUmp0ZlCCLunWHIDmDbRwjYhEFp0ZlKAJw7syYXjXsMsQETlj5ebMwMyGmdkmM4s3swml9T4XdYnijos6ldbhRUQiUrk4MzCzqsA/gO8BCcAyM5vj7utL+r2mjulf0ocUEYl45eXMoD8Q7+5b3T0LmAGMDLkmEZFKo7yEQWtgV4HnCUHb15jZODNbbmbLU1JSyqw4EZGKrryEQaG4+0R3j3P3uKgojbcXESkp5SUMEoHoAs/bBG0iIlIGyksYLANizay9mdUARgFzQq5JRKTSKBejidw9x8zuAOYBVYEp7r4u5LJERCqNchEGAO4+F5gbdh0iIpVReekmEhGREJm7h11DkZhZCrAjeNoU2BdiOcWh2sOh2stepNYNFav2du7+reGYERsGBZnZcnePC7uOolDt4VDtZS9S64bKUbu6iURERGEgIiIVJwwmhl1AMaj2cKj2shepdUMlqL1CXDMQEZHiqShnBiIiUgwKAxERiewwMLNoM/vIzNab2TozuzPsmgrLzGqZ2VIz+yKo/U9h13QmzKyqma00s7fCruVMmNl2M1tjZqvMbHnY9ZwJM2toZq+a2UYz22Bm54ddU2GYWZfg3/vYV5qZ/TLsugrLzO4KfkfXmtnLZhYRi5yb2Z1BzesK8+8d0dcMzKwl0NLdPzezesAK4MrSWCGtpJmZAXXd/bCZVQcWAXe6++KQSysUM/sVEAfUd/fLw66nsMxsOxDn7hE3gcjMpgML3X1ScEPHOu5+IOSyzkiwqmEicJ677zjd9mEzs9bk/26e7e5HzWwWMNfdp4Vb2amZWQ/yFwnrD2QB7wK3unv8yfaJ6DMDd09y98+Dx4eADZxgUZzyyPMdDp5WD74iIpnNrA1wGTAp7FoqCzNrAAwCJgO4e1akBUFgCPBlJARBAdWA2mZWDagD7A65nsLoBixx93R3zwE+Bq4+1Q4RHQYFmVkM0BtYEnIphRZ0tawCkoH57h4ptT8B3APkhVxHUTjwnpmtMLNxYRdzBtoDKcDUoHtukpnVDbuoIhgFvBx2EYXl7onAY8BOIAk46O7vhVtVoawFvmNmTcysDjCCr68Z8y0VIgzM7CzgNeCX7p4Wdj2F5e657n4u+Yv59A9O7co1M7scSHb3FWHXUkQXunsfYDhwu5kNCrugQqoG9AH+6e69gSPAhHBLOjNB19YVwCth11JYZtaI/PXY2wOtgLpm9uNwqzo9d98APAK8R34X0Sog91T7RHwYBP3trwEvuvvssOspiuB0/yNgWMilFMZA4Iqg730GcLGZ/Svckgov+KSHuycDr5PfpxoJEoCEAmePr5IfDpFkOPC5u+8Nu5AzcAmwzd1T3D0bmA1cEHJNheLuk929r7sPAvYDm0+1fUSHQXARdjKwwd3/GnY9Z8LMosysYfC4NvA9YGOoRRWCu9/r7m3cPYb8U/4P3b3cf1ICMLO6wUADgi6WoeSfTpd77r4H2GVmXYKmIUC5HyjxDdcTQV1EgZ3AADOrE/y9GUL+tclyz8yaBd/bkn+94KVTbV9uFrcpooHAjcCaoO8d4L5goZzyriUwPRhdUQWY5e4RNUwzAjUHXs//naYa8JK7vxtuSWdkPPBi0N2yFRgTcj2FFoTv94Bbwq7lTLj7EjN7FfgcyAFWEjm3pnjNzJoA2cDtpxtwENFDS0VEpGREdDeRiIiUDIWBiIgoDERERGEgIiIoDEREBIWBiIigMBAREeD/A3F3wEckhKdyAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "de = [float(smear['log10(E_nu/GeV)_max'][i]) - float((smear['log10(E_nu/GeV)_min'][i])) for i in range(0,len(smear['log10(E_nu/GeV)_max']))]\n",
    "dphi = [float(smear['Dec_nu_max[deg]'][i]) - float((smear['Dec_nu_min[deg]'][i])) for i in range(0,len(smear['Dec_nu_max[deg]']))]\n",
    "#print(de)\n",
    "#print(de)\n",
    "\n",
    "e = [((float(smear['log10(E_nu/GeV)_max'][i]) + float((smear['log10(E_nu/GeV)_min'][i])))/2.0) for i in range(0,len(smear['log10(E_nu/GeV)_max']))]\n",
    "e2 = [pow((float(smear['log10(E_nu/GeV)_max'][i]) + float((smear['log10(E_nu/GeV)_min'][i])))/2.0,2) for i in range(0,len(smear['log10(E_nu/GeV)_max']))]\n",
    "e2dphide = [e2[i] * dphi[i]/de[i] for i in range(0,len(smear['log10(E_nu/GeV)_max']))]\n",
    "\n",
    "plt.figure()\n",
    "plt.plot(e,e2dphide)\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
  },
  "kernelspec": {
   "display_name": "Python 3.10.4 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.4"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}