Newer
Older
HCAL_project / create_true_images.ipynb
@Davide Lancierini Davide Lancierini on 2 Dec 2018 23 KB First commit
{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import os\n",
    "import pickle\n",
    "import sys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_piplus = pd.read_csv('/disk/lhcb_data/davide/HCAL_project_full_event/csv/piplus/MCtrue_piplus.csv', index_col=False)\n",
    "df_piplus0 = pd.read_csv('/disk/lhcb_data/davide/HCAL_project_full_event/csv/piplus0/MCtrue_piplus0.csv', index_col=False)\n",
    "df_Kminus = pd.read_csv('/disk/lhcb_data/davide/HCAL_project_full_event/csv/Kminus/MCtrue_Kminus.csv', index_col=False)\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "true_events_piplus=df_piplus.to_dict(orient='list')\n",
    "true_events_piplus0=df_piplus0.to_dict(orient='list')\n",
    "true_events_Kminus=df_Kminus.to_dict(orient='list')\n",
    "\n",
    "for key in true_events_piplus:\n",
    "    true_events_piplus[key]=np.array(true_events_piplus[key])\n",
    "    true_events_piplus0[key]=np.array(true_events_piplus0[key])\n",
    "    true_events_Kminus[key]=np.array(true_events_Kminus[key])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#DETERMINE MAX AND MIN X,Y\n",
    "max_X=np.max(true_events_piplus['true_x'][np.where(true_events_piplus['region']>=0)])\n",
    "min_X=np.min(true_events_piplus['true_x'][np.where(true_events_piplus['region']>=0)])\n",
    "\n",
    "#width_X=np.ceil(max_X-min_X)\n",
    "width_X=8403.0\n",
    "\n",
    "max_Y=np.max(true_events_piplus['true_y'][np.where(true_events_piplus['region']>=0)])\n",
    "min_Y=np.min(true_events_piplus['true_y'][np.where(true_events_piplus['region']>=0)])\n",
    "\n",
    "#width_Y = np.ceil(max_Y-min_Y)\n",
    "width_Y=6812.0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#number of events\n",
    "batch_size=25000\n",
    "\n",
    "n_batches=2\n",
    "save_test_images = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "piplus = np.array([{} for i in range(n_batches)])\n",
    "piplus0 = np.array([{} for i in range(n_batches)])\n",
    "Kminus = np.array([{} for i in range(n_batches)])\n",
    "\n",
    "for i in range(n_batches):\n",
    "    \n",
    "    for key in true_events_piplus:\n",
    "        piplus[i][key]= true_events_piplus[key][i*batch_size:(i+1)*batch_size]\n",
    "    for key in true_events_piplus0:\n",
    "        piplus0[i][key]= true_events_piplus0[key][i*batch_size:(i+1)*batch_size]\n",
    "    for key in true_events_Kminus:\n",
    "        Kminus[i][key]= true_events_Kminus[key][i*batch_size:(i+1)*batch_size]\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reconverting to np.arrays\n",
    "def create_grid_and_normalise(particle, X_pixels=64, width_X=width_X, width_Y=width_Y):\n",
    "    \n",
    "    particle_norm=np.array([{} for i in range(n_batches)])\n",
    "\n",
    "    for i in range(n_batches):\n",
    "        particle_norm[i][\"true_x\"]=particle[i][\"true_x\"]/(width_X/2)\n",
    "        particle_norm[i][\"true_y\"]=particle[i][\"true_y\"]/(width_Y/2)\n",
    "    \n",
    "        particle_norm[i][\"true_ET\"]=particle[i][\"true_ET\"]\n",
    "        particle_norm[i][\"region\"]=particle[i][\"region\"]\n",
    " #   for i in range(n_batches):\n",
    " #       for j in range(particle[i][\"region\"].shape[0]):\n",
    " #           if particle[i][\"region\"][j]>=0:\n",
    "                \n",
    "    Y_pixels=np.int(np.floor((X_pixels*width_Y)/width_X))\n",
    "\n",
    "    x_values = np.linspace(-X_pixels//2-1,np.ceil(X_pixels/2)+1,X_pixels+3)\n",
    "    y_values = np.linspace(-Y_pixels//2-1,np.ceil(Y_pixels/2)+1,Y_pixels+3)\n",
    "\n",
    "    X_grid, Y_grid = np.meshgrid(x_values, y_values)\n",
    "    \n",
    "    return X_grid, Y_grid, particle_norm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_grid, Y_grid, piplus_norm = create_grid_and_normalise(piplus)\n",
    "_, _,  piplus0_norm= create_grid_and_normalise(piplus0)\n",
    "_, _, Kminus_norm = create_grid_and_normalise(Kminus)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "52"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_pixels = 64\n",
    "Y_pixels = np.int(np.ceil((X_pixels*width_Y)/width_X))\n",
    "Y_pixels"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "def make_pics_dict(X_grid, Y_grid, particle, n_batches=n_batches):\n",
    "    \n",
    "    pic = np.zeros(shape=(Y_pixels,X_pixels,1),dtype=np.float32)\n",
    "    pics_dict = {}\n",
    "        \n",
    "    for n in range(n_batches):\n",
    "        \n",
    "        batch_size = particle[n][\"region\"].shape[0]\n",
    "        pics_dict[n]=np.array([pic for i in range(0,batch_size)])\n",
    "        \n",
    "        for k in range(batch_size):\n",
    "            event=k\n",
    "            if particle[n][\"region\"][event]>=0: #and particle_1[n][\"region\"][event]>=0 and particle_2[n][\"region\"][event]>=0:\n",
    "                #event = n*batch_size+k\n",
    "                \n",
    "                #print(event)\n",
    "                for i in range(X_grid.shape[0]):\n",
    "                    for j in range(Y_grid.shape[1]):\n",
    "                                \n",
    "                        if X_grid[i,j] < particle[n][\"true_x\"][event]*X_pixels/2 < X_grid[i,j+1]:\n",
    "                            if Y_grid[i,j] < particle[n][\"true_x\"][event]*Y_pixels/2 < Y_grid[i+1,j]:\n",
    "                                \n",
    "                                x_primed=int(Y_pixels/2-Y_grid[i,j])\n",
    "                                y_primed=int(X_pixels/2+X_grid[i,j])\n",
    "                                \n",
    "                                #pics_dict[n][k][x_primed,y_primed,0]=1\n",
    "    \n",
    "                                pics_dict[n][k][x_primed-1,y_primed-1,0]=particle[n][\"true_ET\"][event]\n",
    "                                #total_pic[x_primed-1,y_primed-1,0]=1\n",
    "                                #total_pic[x_primed-1,y_primed-1,0]=true_events['true_ET'][event]\n",
    "    \n",
    "                            \n",
    "        print('Converted '+str(n+1)+'/'+str(n_batches)+' batches of ' +str(batch_size)+' images')\n",
    "    \n",
    "  #  \n",
    "  #  if N % batch_size != 0:\n",
    "  #      \n",
    "  #      print('Converting last batch of '+str(N%batch_size)+' images') \n",
    "  #  \n",
    "  #      for k in range(batch_size*n_batches,N):\n",
    "  #          #if true_events['region'][k]>=0:\n",
    "  #  \n",
    "  #              for i in range(X_grid.shape[0]):\n",
    "  #                  for j in range(Y_grid.shape[1]):\n",
    "  #                      if X_grid[i,j-1] < X_norm[n_batches][k]*X_pixels/2 < X_grid[i,j]:\n",
    "  #                          if Y_grid[i-1,j] < Y_norm[n_batches][k]*Y_pixels/2 < Y_grid[i,j]:\n",
    "  #                              \n",
    "  #                              x_primed=int(Y_pixels/2-Y_grid[i,j])\n",
    "  #                              y_primed=int(X_pixels/2+X_grid[i,j])\n",
    "  #                              \n",
    "  #                              pics_dict[n_batches][k-(N-1)][x_primed,y_primed,0]=true_events['true_ET'][event]\n",
    "  #                              #total_pic[x_primed,y_primed,0]=true_events['true_ET'][event]\n",
    "    return pics_dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 1/2 batches of 25000 images\n",
      "Converted 2/2 batches of 25000 images\n"
     ]
    }
   ],
   "source": [
    "pics_dict_piplus = make_pics_dict(X_grid, Y_grid, piplus_norm, n_batches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 1/2 batches of 25000 images\n",
      "Converted 2/2 batches of 25000 images\n"
     ]
    }
   ],
   "source": [
    "pics_dict_piplus0 = make_pics_dict(X_grid, Y_grid, piplus0_norm, n_batches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Converted 1/2 batches of 25000 images\n",
      "Converted 2/2 batches of 25000 images\n"
     ]
    }
   ],
   "source": [
    "pics_dict_Kminus = make_pics_dict(X_grid, Y_grid, Kminus_norm, n_batches)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(25000, 52, 64, 1)"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pics_dict_piplus[1].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7f1b3b509250>"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAAD8CAYAAAAfZJO2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAADJVJREFUeJzt3W+IZYV5x/Hvr7vrmpiKmqhsXammbFN9UdewGIOltBpbm4ToC1OUUPbFwr6xYGgg1RYKgb6Ib6J9UQpLtNkXadSapIqE2GWjlEJZHaMm6kZXrY3LWjdtlaSBbnfN0xf3aEe949z5d+c++P3AcO85cy7nYe/s13PPHM+mqpCkrn5pvQeQpJUwYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaW1HEklyV5JkkzyW5abWGkqRJZblX7CfZADwLXAkcBh4Brq+qpxd6zUnZXCdzyrL2J+m943/4Of9bxzLJthtXsJ9LgOeq6gWAJHcCVwMLRuxkTuFjuWIFu5T0XnCg9k+87Uo+Tp4DvDRv+fCw7i2S7E4yl2TuOMdWsDtJeqeVRGzcod47PptW1Z6q2lFVOzaxeQW7k6R3WknEDgPnzlveChxZ2TiStDQridgjwLYk5yc5CbgOuG91xpKkySz7xH5VnUjyx8ADwAbgjqp6atUmk6QJrOS3k1TVd4DvrNIskrRkXrEvqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyapNSMmqTUjJqk1IyaptUUjluSOJEeTPDlv3RlJ9iU5NDyevrZjStJ4kxyJfQ246m3rbgL2V9U2YP+wLElTt2jEquqfgP962+qrgb3D873ANas8lyRNZLnnxM6uqpcBhsezFtowye4kc0nmjnNsmbuTpPHW/MR+Ve2pqh1VtWMTm9d6d5LeY5YbsVeSbAEYHo+u3kiSNLnlRuw+YOfwfCdw7+qMI0lLM8klFt8A/gX4SJLDSXYBXwauTHIIuHJYlqSp27jYBlV1/QLfumKVZ5GkJfOKfUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrS0asSTnJnkwycEkTyW5cVh/RpJ9SQ4Nj6ev/biS9FaTHImdAL5QVRcAlwI3JLkQuAnYX1XbgP3DsiRN1aIRq6qXq+r7w/OfAQeBc4Crgb3DZnuBa9ZqSElayJLOiSU5D7gYOACcXVUvwyh0wFmrPZwkLWbiiCX5APBN4PNV9dMlvG53krkkc8c5tpwZJWlBE0UsySZGAft6VX1rWP1Kki3D97cAR8e9tqr2VNWOqtqxic2rMbMkvWmS304GuB04WFVfmfet+4Cdw/OdwL2rP54kvbuNE2xzGfBHwA+TPD6s+zPgy8DdSXYBPwY+uzYjStLCFo1YVf0zkAW+fcXqjiNJS+MV+5JaM2KSWjNiklozYpJaM2KSWpvkEgtpXT1w5PF3rPv9X9m+DpNoFnkkJqk1IyapNSMmqTUjJqk1T+xr5nkSX+/GIzFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrRkxSa0ZMUmtGTFJrS0asSQnJ3k4yRNJnkrypWH9+UkOJDmU5K4kJ639uJL0VpMciR0DLq+qi4DtwFVJLgVuAW6tqm3Aq8CutRtTksZbNGI18t/D4qbhq4DLgXuG9XuBa9ZkQk3Vhl//tbd8SbNuonNiSTYkeRw4CuwDngdeq6oTwyaHgXPWZkRJWthEEauq16tqO7AVuAS4YNxm416bZHeSuSRzxzm2/EklaYwl/Xayql4DHgIuBU5LsnH41lbgyAKv2VNVO6pqxyY2r2RWSXqHjYttkORM4HhVvZbkfcAnGJ3UfxC4FrgT2Ancu5aDajpef/b59R5BWpJFIwZsAfYm2cDoyO3uqro/ydPAnUn+EngMuH0N55SksRaNWFX9ALh4zPoXGJ0fk6R14xX7klozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWjNiklozYpJaM2KSWps4Ykk2JHksyf3D8vlJDiQ5lOSuJCet3ZiSNN5SjsRuBA7OW74FuLWqtgGvArtWczBJmsREEUuyFfgU8NVhOcDlwD3DJnuBa9ZiQEl6N5Meid0GfBH4xbD8QeC1qjoxLB8Gzlnl2SRpUYtGLMmngaNV9ej81WM2rQVevzvJXJK54xxb5piSNN7GCba5DPhMkk8CJwOnMjoyOy3JxuFobCtwZNyLq2oPsAfg1JwxNnSStFyLHolV1c1VtbWqzgOuA75XVZ8DHgSuHTbbCdy7ZlNK0gJWcp3YnwJ/kuQ5RufIbl+dkSRpcpN8nHxTVT0EPDQ8fwG4ZPVHkqTJecW+pNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1jZOslGSF4GfAa8DJ6pqR5IzgLuA84AXgT+sqlfXZkxJGm8pR2K/W1Xbq2rHsHwTsL+qtgH7h2VJmqqVfJy8Gtg7PN8LXLPycSRpaSaNWAH/mOTRJLuHdWdX1csAw+NZazGgJL2bic6JAZdV1ZEkZwH7kvxo0h0M0dsNcDLvX8aIkrSwiY7EqurI8HgU+DZwCfBKki0Aw+PRBV67p6p2VNWOTWxenaklabBoxJKckuSX33gO/B7wJHAfsHPYbCdw71oNKUkLmeTj5NnAt5O8sf3fVdV3kzwC3J1kF/Bj4LNrN6YkjbdoxKrqBeCiMev/E7hiLYaSpEl5xb6k1lJV09tZ8hPg34APAf8xtR0vn3Oung4zgnOutuXO+atVdeYkG041Ym/uNJmbd+X/zHLO1dNhRnDO1TaNOf04Kak1IyaptfWK2J512u9SOefq6TAjOOdqW/M51+WcmCStFj9OSmptqhFLclWSZ5I8l2Sm7j+W5I4kR5M8OW/dGUn2JTk0PJ6+zjOem+TBJAeTPJXkxhmd8+QkDyd5YpjzS8P685McGOa8K8lJ6znnMNOGJI8luX+GZ3wxyQ+TPJ5kblg3U+/5MNNpSe5J8qPhZ/Tj05hzahFLsgH4a+APgAuB65NcOK39T+BrwFVvWzdrN348AXyhqi4ALgVuGP4MZ23OY8DlVXURsB24KsmlwC3ArcOcrwK71nHGN9wIHJy3PIszQo+bkv4V8N2q+g1G/5fPQaYxZ1VN5Qv4OPDAvOWbgZuntf8JZzwPeHLe8jPAluH5FuCZ9Z7xbfPeC1w5y3MC7we+D3yM0UWPG8f9PKzTbFuHv1iXA/cDmbUZhzleBD70tnUz9Z4DpwL/ynCefZpzTvPj5DnAS/OWDw/rZtnM3vgxyXnAxcABZnDO4WPa44xu0bQPeB54rapODJvMwvt/G/BF4BfD8geZvRmhx01JPwz8BPjb4eP5V4e73qz5nNOMWMas81ejy5DkA8A3gc9X1U/Xe55xqur1qtrO6GjnEuCCcZtNd6r/l+TTwNGqenT+6jGbzsLP6GVV9VFGp2JuSPLb6z3QGBuBjwJ/U1UXAz9nSh9xpxmxw8C585a3AkemuP/lmOjGj9OUZBOjgH29qr41rJ65Od9QVa8BDzE6h3dakjfunLLe7/9lwGeGf8nrTkYfKW9jtmYEVnZT0ik6DByuqgPD8j2Morbmc04zYo8A24bf/pwEXMfoxoqzbKZu/JjRTd1uBw5W1VfmfWvW5jwzyWnD8/cBn2B0kvdB4Nphs3Wds6purqqtVXUeo5/F71XV55ihGaHPTUmr6t+Bl5J8ZFh1BfA005hzyif/Pgk8y+j8yJ+v54nIMbN9A3gZOM7ovyq7GJ0j2Q8cGh7PWOcZf4vRx5sfAI8PX5+cwTl/E3hsmPNJ4C+G9R8GHgaeA/4e2Lze7/sw1+8A98/ijMM8TwxfT73x92bW3vNhpu3A3PC+/wNw+jTm9Ip9Sa15xb6k1oyYpNaMmKTWjJik1oyYpNaMmKTWjJik1oyYpNb+D1TuIblof2KtAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "i=1\n",
    "plt.imshow(pics_dict_piplus[1][i].reshape(52,64)+pics_dict_piplus0[1][i].reshape(52,64)+pics_dict_Kminus[0][i].reshape(52,64))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "pics_dict={}\n",
    "\n",
    "pics_dict[0]=pics_dict_piplus[0]+pics_dict_piplus0[0]+pics_dict_Kminus[0]\n",
    "pics_dict[1]=pics_dict_piplus[1]+pics_dict_piplus0[1]+pics_dict_Kminus[1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "for i in range(n_batches):\n",
    "    with open('/disk/lhcb_data/davide/HCAL_project_full_event/true/sample'+str(i)+'.pickle', 'wb') as handle:\n",
    "        pickle.dump(pics_dict[i], handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "a=np.array([pics_dict_Kminus[0][i].sum()+pics_dict_piplus[0][i].sum()+pics_dict_piplus0[0][i].sum() for i in range(len(pics_dict_Kminus[0]))])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAD8CAYAAAB+UHOxAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAE21JREFUeJzt3X+sXOWd3/H3Zx3yQ5uomHChru3UdOuqIVVjoltAon+wJAsGVjUrLRWo3XgjVG8lUBMpbdfkHzZJkRypG5pIWSRvcGNW2RArP4oF3mVdkijNHwEMIQRDIrzEDbe28N0aSFBUKrPf/jGPm8G+9p17PfdefJ/3SxrNOd/znJnnkcf3M885Z2ZSVUiS+vNrS90BSdLSMAAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnXrLUnfgdM4///xat27dUndDks4qjz/++N9U1cRs7d7UAbBu3Tr27du31N2QpLNKkv85SjsPAUlSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqfe1J8EPlus2/rgG9YPbrt+iXoiSaNzBiBJnTIAJKlTBoAkdcoAkKROGQCS1KlZAyDJ25M8muSHSfYn+WSrfynJT5M82W4bWj1JPp/kQJKnknxg6LE2J3mu3TYv3LAkSbMZ5TLQ14CrqurVJOcA30vyF23bf6iqr53Q/lpgfbtdBtwNXJbkPOAOYBIo4PEku6vqpXEMRJI0N7POAGrg1bZ6TrvVaXbZBNzb9vs+cG6SVcA1wN6qOtr+6O8FNp5Z9yVJ8zXSOYAkK5I8CRxh8Ef8kbbpznaY564kb2u11cALQ7tPtdqp6pKkJTBSAFTV61W1AVgDXJrknwC3A/8Y+GfAecAftuaZ6SFOU3+DJFuS7Euyb3p6epTuSZLmYU5XAVXVy8B3gI1Vdbgd5nkN+K/Apa3ZFLB2aLc1wKHT1E98ju1VNVlVkxMTs/6ovSRpnka5Cmgiyblt+R3Ah4Aft+P6JAlwA/B022U38OF2NdDlwCtVdRh4CLg6ycokK4GrW02StARGuQpoFbAzyQoGgbGrqh5I8q0kEwwO7TwJ/NvWfg9wHXAA+CXwEYCqOprk08Bjrd2nquro+IYiSZqLWQOgqp4CLpmhftUp2hdw6ym27QB2zLGPkqQF4CeBJalTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkTo3yZXCao3VbH3zD+sFt1y9RTyTp1JwBSFKnnAHM0Ynv7iXpbOUMQJI6ZQBIUqcMAEnqlAEgSZ0yACSpU7MGQJK3J3k0yQ+T7E/yyVa/KMkjSZ5L8tUkb231t7X1A237uqHHur3Vf5LkmoUalCRpdqPMAF4Drqqq9wMbgI1JLgc+A9xVVeuBl4BbWvtbgJeq6h8Cd7V2JLkYuAl4H7AR+JMkK8Y5GEnS6GYNgBp4ta2e024FXAV8rdV3Aje05U1tnbb9g0nS6vdV1WtV9VPgAHDpWEYhSZqzkc4BJFmR5EngCLAX+Gvg5ao61ppMAavb8mrgBYC2/RXg3cP1GfaRJC2ykQKgql6vqg3AGgbv2t87U7N2n1NsO1X9DZJsSbIvyb7p6elRuidJmoc5XQVUVS8D3wEuB85NcvyrJNYAh9ryFLAWoG3/O8DR4foM+ww/x/aqmqyqyYmJibl0T5I0B6NcBTSR5Ny2/A7gQ8CzwLeB323NNgP3t+XdbZ22/VtVVa1+U7tK6CJgPfDouAYiSZqbUb4MbhWws12x82vArqp6IMkzwH1J/hPwA+Ce1v4e4M+SHGDwzv8mgKran2QX8AxwDLi1ql4f73AkSaOaNQCq6ingkhnqzzPDVTxV9X+AG0/xWHcCd869m5KkcfOTwJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqcMAEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tSsPwqfZC1wL/B3gb8FtlfV55L8EfBvgOnW9BNVtaftcztwC/A68O+q6qFW3wh8DlgBfLGqto13OG9O67Y+eFLt4Lbrl6AnkvQrswYAcAz4eFU9keRdwONJ9rZtd1XVfx5unORi4CbgfcDfA/57kn/UNn8B+C1gCngsye6qemYcA5Ekzc2sAVBVh4HDbfkXSZ4FVp9ml03AfVX1GvDTJAeAS9u2A1X1PECS+1pbA0CSlsCczgEkWQdcAjzSSrcleSrJjiQrW2018MLQblOtdqr6ic+xJcm+JPump6dP3CxJGpORAyDJO4GvAx+rqp8DdwO/AWxgMEP44+NNZ9i9TlN/Y6Fqe1VNVtXkxMTEqN2TJM3RKOcASHIOgz/+X66qbwBU1YtD2/8UeKCtTgFrh3ZfAxxqy6eqS5IW2awzgCQB7gGerarPDtVXDTX7HeDptrwbuCnJ25JcBKwHHgUeA9YnuSjJWxmcKN49nmFIkuZqlBnAFcDvAT9K8mSrfQK4OckGBodxDgJ/AFBV+5PsYnBy9xhwa1W9DpDkNuAhBpeB7qiq/WMciyRpDka5Cuh7zHz8fs9p9rkTuHOG+p7T7SdJWjx+EliSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6NdKPwvds3dYHl7oLkrQgnAFIUqdmDYAka5N8O8mzSfYn+Wirn5dkb5Ln2v3KVk+Szyc5kOSpJB8YeqzNrf1zSTYv3LAkSbMZ5RDQMeDjVfVEkncBjyfZC/w+8HBVbUuyFdgK/CFwLbC+3S4D7gYuS3IecAcwCVR7nN1V9dK4B3U2OPHQ0sFt1y9RTyT1atYZQFUdrqon2vIvgGeB1cAmYGdrthO4oS1vAu6tge8D5yZZBVwD7K2qo+2P/l5g41hHI0ka2ZzOASRZB1wCPAJcWFWHYRASwAWt2WrghaHdplrtVPUTn2NLkn1J9k1PT8+le5KkORg5AJK8E/g68LGq+vnpms5Qq9PU31io2l5Vk1U1OTExMWr3JElzNFIAJDmHwR//L1fVN1r5xXZoh3Z/pNWngLVDu68BDp2mLklaAqNcBRTgHuDZqvrs0KbdwPEreTYD9w/VP9yuBroceKUdInoIuDrJynbF0NWtJklaAqNcBXQF8HvAj5I82WqfALYBu5LcAvwMuLFt2wNcBxwAfgl8BKCqjib5NPBYa/epqjo6llFIkuZs1gCoqu8x8/F7gA/O0L6AW0/xWDuAHXPpoCRpYfhJYEnqlAEgSZ0yACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6pQBIEmdMgAkqVMGgCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnRrlN4G1CNZtffAN6we3Xb9EPZHUi1lnAEl2JDmS5Omh2h8l+V9Jnmy364a23Z7kQJKfJLlmqL6x1Q4k2Tr+oUiS5mKUQ0BfAjbOUL+rqja02x6AJBcDNwHva/v8SZIVSVYAXwCuBS4Gbm5tJUlLZNZDQFX13STrRny8TcB9VfUa8NMkB4BL27YDVfU8QJL7Wttn5txjSdJYnMlJ4NuSPNUOEa1stdXAC0NtplrtVPWTJNmSZF+SfdPT02fQPUnS6cw3AO4GfgPYABwG/rjVM0PbOk395GLV9qqarKrJiYmJeXZPkjSbeV0FVFUvHl9O8qfAA211Clg71HQNcKgtn6ouSVoC85oBJFk1tPo7wPErhHYDNyV5W5KLgPXAo8BjwPokFyV5K4MTxbvn321J0pmadQaQ5CvAlcD5SaaAO4Ark2xgcBjnIPAHAFW1P8kuBid3jwG3VtXr7XFuAx4CVgA7qmr/2EcjSRrZKFcB3TxD+Z7TtL8TuHOG+h5gz5x6J0laMH4VhCR1ygCQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6Na/fA9DCW7f1wZNqB7ddvwQ9kbRcOQOQpE4ZAJLUKQNAkjplAEhSpwwASeqUASBJnZo1AJLsSHIkydNDtfOS7E3yXLtf2epJ8vkkB5I8leQDQ/tsbu2fS7J5YYYjSRrVKDOALwEbT6htBR6uqvXAw20d4FpgfbttAe6GQWAAdwCXAZcCdxwPDUnS0pg1AKrqu8DRE8qbgJ1teSdww1D93hr4PnBuklXANcDeqjpaVS8Bezk5VCRJi2i+5wAurKrDAO3+glZfDbww1G6q1U5VlyQtkXGfBM4MtTpN/eQHSLYk2Zdk3/T09Fg7J0n6lfkGwIvt0A7t/kirTwFrh9qtAQ6dpn6SqtpeVZNVNTkxMTHP7kmSZjPfANgNHL+SZzNw/1D9w+1qoMuBV9ohooeAq5OsbCd/r241SdISmfXbQJN8BbgSOD/JFIOrebYBu5LcAvwMuLE13wNcBxwAfgl8BKCqjib5NPBYa/epqjrxxLIkaRHNGgBVdfMpNn1whrYF3HqKx9kB7JhT7yRJC8ZPAktSp/xBmLPIiT8S4w/ESDoTzgAkqVMGgCR1ykNAJ5jpt3glaTlyBiBJnTIAJKlTBoAkdcoAkKROGQCS1CkDQJI6ZQBIUqf8HMBZbKbPLPj1EJJG5QxAkjplAEhSpwwASeqUASBJnTIAJKlTBoAkdeqMLgNNchD4BfA6cKyqJpOcB3wVWAccBP5lVb2UJMDnGPxo/C+B36+qJ87k+XUyfzVM0qjGMQP4zaraUFWTbX0r8HBVrQcebusA1wLr220LcPcYnluSNE8LcQhoE7CzLe8Ebhiq31sD3wfOTbJqAZ5fkjSCMw2AAv4qyeNJtrTahVV1GKDdX9Dqq4EXhvadajVJ0hI406+CuKKqDiW5ANib5MenaZsZanVSo0GQbAF4z3vec4bdkySdyhkFQFUdavdHknwTuBR4McmqqjrcDvEcac2ngLVDu68BDs3wmNuB7QCTk5MnBYTmxpPCkk5l3oeAkvx6kncdXwauBp4GdgObW7PNwP1teTfw4QxcDrxy/FCRJGnxnckM4ELgm4OrO3kL8OdV9ZdJHgN2JbkF+BlwY2u/h8EloAcYXAb6kTN4bknSGZp3AFTV88D7Z6j/b+CDM9QLuHW+zydJGi8/CSxJnTIAJKlTBoAkdcoAkKROGQCS1Cl/FL4z/pC8pOOcAUhSp5wByK+LkDrlDECSOtX1DGCm4+GS1AtnAJLUKQNAkjrV9SEgzcyTwlIfDADNys8OSMuTh4AkqVPOADQvHiaSzn7OACSpU84ANBbOCKSzT1cB4Ae/Fo8njqU3v0UPgCQbgc8BK4AvVtW2xe6DloazBOnNZVEDIMkK4AvAbwFTwGNJdlfVM4vZD705jBIIhoa0cBZ7BnApcKCqngdIch+wCViQAPCQz9lllH8vA0Ean8UOgNXAC0PrU8Bli9wHLSMLFfLzmY28mc57GJQaxWIHQGao1RsaJFuALW311SQ/OYPnOx/4mzPY/2zluM9QPrN4bcZg1nEvUj8Wm6/zU/v7ozzQYgfAFLB2aH0NcGi4QVVtB7aP48mS7KuqyXE81tnEcffFcfdlnONe7A+CPQasT3JRkrcCNwG7F7kPkiQWeQZQVceS3AY8xOAy0B1VtX8x+yBJGlj0zwFU1R5gzyI93VgOJZ2FHHdfHHdfxjbuVNXsrSRJy45fBidJnVqWAZBkY5KfJDmQZOtS92chJdmR5EiSp4dq5yXZm+S5dr9yKfs4bknWJvl2kmeT7E/y0VZf1uMGSPL2JI8m+WEb+ydb/aIkj7Sxf7VdZLGsJFmR5AdJHmjry37MAEkOJvlRkieT7Gu1sbzWl10ADH3dxLXAxcDNSS5e2l4tqC8BG0+obQUerqr1wMNtfTk5Bny8qt4LXA7c2v6Nl/u4AV4Drqqq9wMbgI1JLgc+A9zVxv4ScMsS9nGhfBR4dmi9hzEf95tVtWHo8s+xvNaXXQAw9HUTVfV/geNfN7EsVdV3gaMnlDcBO9vyTuCGRe3UAquqw1X1RFv+BYM/CqtZ5uMGqIFX2+o57VbAVcDXWn3ZjT3JGuB64IttPSzzMc9iLK/15RgAM33dxOol6stSubCqDsPgjyVwwRL3Z8EkWQdcAjxCJ+Nuh0KeBI4Ae4G/Bl6uqmOtyXJ8zf8X4D8Cf9vW383yH/NxBfxVksfbNyXAmF7ry/H3AGb9ugktD0neCXwd+FhV/XzwpnD5q6rXgQ1JzgW+Cbx3pmaL26uFk+S3gSNV9XiSK4+XZ2i6bMZ8giuq6lCSC4C9SX48rgdejjOAWb9uogMvJlkF0O6PLHF/xi7JOQz++H+5qr7Ryst+3MOq6mXgOwzOg5yb5PgbuuX2mr8C+BdJDjI4pHsVgxnBch7z/1dVh9r9EQaBfyljeq0vxwDw6yYG493cljcD9y9hX8auHf+9B3i2qj47tGlZjxsgyUR750+SdwAfYnAO5NvA77Zmy2rsVXV7Va2pqnUM/j9/q6r+Fct4zMcl+fUk7zq+DFwNPM2YXuvL8oNgSa5j8A7h+NdN3LnEXVowSb4CXMngGwJfBO4A/huwC3gP8DPgxqo68UTxWSvJPwf+B/AjfnVM+BMMzgMs23EDJPmnDE76rWDwBm5XVX0qyT9g8O74POAHwL+uqteWrqcLox0C+vdV9ds9jLmN8Ztt9S3An1fVnUnezRhe68syACRJs1uOh4AkSSMwACSpUwaAJHXKAJCkThkAktQpA0CSOmUASFKnDABJ6tT/A7wp3R/F2ulpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.hist(a/1000,bins=70);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'r' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-37-39bb05a3998a>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTBrowser\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'r' is not defined"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}