{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pylhe\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'particles': [{'status': -1.0, 'e': 1.0, 'mother1': 0.0, 'mother2': 0.0, 'pz': 1.0, 'px': 0.0, 'py': 0.0, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': -1.0, 'id': -11.0}, {'status': -1.0, 'e': 1.0, 'mother1': 0.0, 'mother2': 0.0, 'pz': -1.0, 'px': -0.0, 'py': -0.0, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': 1.0, 'id': 11.0}, {'status': 1.0, 'e': 1.0, 'mother1': 1.0, 'mother2': 2.0, 'pz': -0.32835755731, 'px': 0.019159631178, 'py': -0.94435916001, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': 1.0, 'id': -13.0}, {'status': 1.0, 'e': 1.0, 'mother1': 1.0, 'mother2': 2.0, 'pz': 0.32835755731, 'px': -0.019159631178, 'py': 0.94435916001, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': -1.0, 'id': 13.0}], 'eventinfo': {'scale': 2.0, 'weight': 23213.001, 'pid': 1.0, 'nparticles': 4.0, 'aqed': 0.007546771, 'aqcd': 0.2980729}}\n"
     ]
    }
   ],
   "source": [
    "filepath='myMGTest2/Events/run_01/unweighted_events.lhe'\n",
    "for i, e in enumerate(pylhe.readLHE(filepath)):\n",
    "    if i == 0:\n",
    "        firstevent = e\n",
    "        print e"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'status': 1.0, 'e': 1.0, 'mother1': 1.0, 'mother2': 2.0, 'pz': -0.32835755731, 'px': 0.019159631178, 'py': -0.94435916001, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': 1.0, 'id': -13.0}\n",
      "{'status': 1.0, 'e': 1.0, 'mother1': 1.0, 'mother2': 2.0, 'pz': 0.32835755731, 'px': -0.019159631178, 'py': 0.94435916001, 'm': 0.0, 'color1': 0.0, 'color2': 0.0, 'lifetime': 0.0, 'spin': -1.0, 'id': 13.0}\n"
     ]
    }
   ],
   "source": [
    "for particle in firstevent['particles'][2:]:\n",
    "    print particle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-13.0\n"
     ]
    }
   ],
   "source": [
    "print firstevent['particles'][2]['id']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#make histogram of cos theta\n",
    "mupluscth = np.zeros(10000);\n",
    "muminuscth = np.zeros(10000);\n",
    "\n",
    "for i, e in enumerate(pylhe.readLHE(filepath)):\n",
    "    mupluscth[i] = e['particles'][2]['pz']/np.sqrt(e['particles'][2]['px']**2 + e['particles'][2]['py']**2 + e['particles'][2]['pz']**2)\n",
    "    muminuscth[i] = e['particles'][3]['pz']/np.sqrt(e['particles'][3]['px']**2 + e['particles'][3]['py']**2 + e['particles'][3]['pz']**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAIABJREFUeJzt3Xl8FPX9x/HXZzcHioIoeHGrqI1HPSLWYgUVEVEBxQO8QFFqFW9/CtVipfWu2mpRQcVbEfGKguKJtooWVDwAgYgHIOUQ0QpCsruf3x+ztGtMyAY2mWTyfj4e+2Bn5ru7bybJJ998Z+Y75u6IiEi0xMIOICIiuafiLiISQSruIiIRpOIuIhJBKu4iIhGk4i4iEkEq7iIiEaTiLiISQSruIiIRlBfWB7ds2dI7dOgQ1seLiDRI77333nJ3b1Vdu9CKe4cOHZg+fXpYHy8i0iCZ2ZfZtNOwjIhIBKm4i4hEUFbF3cx6mtkcMys1s2GVbG9nZq+b2Qdm9pGZ9cp9VBERyVa1xd3M4sAo4AigCBhgZkUVml0JjHf3vYH+wB25DioiItnLpufeGSh19/nuXgaMA/pUaONAs/Tz5sDXuYsoIiI1lc3ZMq2BBRnLC4H9K7T5I/CSmZ0HNAW65ySdiIhskFwdUB0A3O/ubYBewENm9rP3NrMhZjbdzKYvW7YsRx8tIiIVZVPcFwFtM5bbpNdlGgyMB3D3qUAToGXFN3L3Me5e7O7FrVpVew6+iIhsoGyK+zSgk5l1NLMCggOmJRXafAUcCmBmvyAo7uqai4iEpNoxd3dPmNlQYDIQB8a6+0wzGwlMd/cS4BLgbjO7iODg6iAP4c7bHYZNrHT9F9cfWcdJRETCldX0A+4+CZhUYd2IjOezgC65jSYiIhtKV6iKiESQiruISASpuIuIRJCKu4hIBKm4i4hEkIq7iEgEqbiLiERQwyvuSz+FF38PqVTYSURE6q2GV9znvw7vjIK3bg07iYhIvdXwivv+Z8Pu/eDVP0Hpq2GnERHJzrxX4LVrIJWsk4/LavqBesUMet8OS2fDk4NhyBRo0SHkUCLSWFU1pxVkzGu14vOgXjVvAwdeBAWb1nquhlfcgQ4jptDezqSk4EoW3tqbfmV/ZA2FYccSEfm5stXw+CnB8xMfrpPCDg1xWCbtS9+WC8rP5Rf2Fdfk30swGaWISD3iDs+dD0tmQr97YcuOdfbRDba4A0xJ7c1fE/3oF/8np8VfCjuOiMhPXH3lefDxE9xUfjwd7l1Lh2ET1zuMk0sNurgD3J7syyvJvflD3sMU26dhxxERAaBL7GOuyHuEycli7kj2rvPPb/DF3Ylxcfk5LPBW3FFwG9uwIuxIItLItbMljMq/jc98ey4u/x0eQqlt8MUd4Hua8tvyi9mUNYwpuIVCysKOJCKNVFN+5O78m3GMs8ovYRWbhJIjEsUdYJ634cLyc9nDPuem/NHoAKuI1DUjxa35d7Cjfc255efzlW8TWpasiruZ9TSzOWZWambDKtl+q5nNSD/mmtnK3Eet3iupfflL4gR6x6dyTrziPbxFRGrXhXlP0SP+Hn9OnMLbqd1DzVLtee5mFgdGAYcBC4FpZlaSvm8qAO5+UUb784C9ayFrVu5I9maX2AIuzRvPXG/DK6l9w4oiIo1Iz9i/uCDvKcYnunJ/8vCw42TVc+8MlLr7fHcvA8YBfdbTfgDwWC7CbRjjsvIhfOwd+Wv+KHaxr8KLIiKNwq72Fbfk38l7qU5cmTgDsLAjZVXcWwMLMpYXptf9jJm1BzoCr1WxfYiZTTez6cuWLatp1qytpYAhZReziibck38zrPqm1j5LRBq3VqzknoK/8B1NObvsQsrIDzsSkPsDqv2BCe5e6cw47j7G3YvdvbhVq1Y5/uifWsKWDCm7mK1tJYw/DRI6g0ZEcqz8R+4uuJkt+Q9nll3CMlqEnei/sinui4C2Gctt0usq059Qh2R+6kPficvKz4Iv/wnPXxRcCiwikgupFDx9NnvafC4sP4eZXndTC2Qjm+I+DehkZh3NrICggP/sVBQz2xVoAUzNbcSN82zqQDjoMpjxMPzj5rDjiEhUvP5nmPUM1yZO4qXUfmGn+Zlqi7u7J4ChwGRgNjDe3Wea2Ugzy7ymtj8wzr0edo8P/j3scTy89if4eELYaUSkofvgkaCzuO8g7kn2CjtNpbKa8tfdJwGTKqwbUWH5j7mLlWNm0GcUfLcQnjknmFO53a/CTiUiDdHn/4DnLoAdukGvv8Bb9XPSwshcoVqtvELo/2hQ2B8bAN98FnYiEWlolpcGc7Nv2RGOfwDi9ePMmMo0nuIOsOmWcPITwfNHT4DVmmRMRLL0wzJ49HiIxeGk8bDJFmEnWq/GVdwBttox6MGv/ArGnQyJtWEnEpH6rmxV0CH8/msYMK5Ob7qxoRpfcQdofwD0vRO+ehuePjs4pUlEpDLJBDxxOiyeAceNhbadw06UlQZ5D9Wc2OM4+H4RvDwCNtsGel4XHHgVEVnHHSZeBPMmw1G3wq5Hhp0oa42iuFd9W6sd+aLbOfDOHbD5tnDghXWaS0TquTdugPcfhN9cCsVnhJ2mRhpFca+aQY9r4Icl8MpVQQ9+rwFhhxKR+uC9B2DKdfDLk+CQK8NOU2ONvLgDsVgw/r5qOZQMhaatoFP3sFOJSJjmTg6mLNnxUOh9W4Mcsm2cB1QryiuEEx+GrYtg/Kmw8L2wE4lIWL56B8YPhG13hxMerNfnsq+Pivs6TZrByROCnvujxwcXK4hI47L4I3jkBGi2fVAPCjcLO9EGa/TDMhUPtnaw85lQcDUtH+wDZ7wIW7St4pUiEinffAYPHxsU9NOehc22DjvRRlHPvYIvfDtOKxsGa/8DD/aG/ywJO5KI1LbvFsGDfcBTcOozkejUNfqee2VmeQc4ZQI82BceOgYGPR9MXSAi0bNqOTzUF9Z8BwOfg1Y7V9qs6lOq6yf13KvStjP0fwS+mQePHBf05EUkWtZ8Dw/3C6YjGTAOtt8r7EQ5o+JehQ7DJtLh7tWc9eNQEgs/4O0/92CXYU+HHUtEcqVsdTBD7JJP4ISHoEOXsBPllIZlqvFyqphLys/m1vw7uSP/b5A8qsGeGiUiaeU/wrgBwfxSx94NO/f476aGNvxSFRX3LDybOpCmibVcm38vTDgjmDyoQoGv6hvii+sbzlwUIo1C+ZpgTvb5bwQXMO5xXNiJakVWwzJm1tPM5phZqZkNq6LNCWY2y8xmmtmjuY0ZvkeThzKy/FSYXQJPnhnMFCciDUuiDMafBqWvQO/bIz3dSLXF3cziwCjgCKAIGGBmRRXadAKGA13cfTcgkjNwjU0eEcxFM+sZeOosFXiRhiRZDhNO/98Mj/ucGnaiWpXNsExnoNTd5wOY2TigDzAro81ZwCh3/xbA3ZfmOmi98euhwbmwL/8hmG/imDEQ1+iWSL2WTMCTg+HT54P7nhafEZmx9apkU5VaAwsylhcC+1doszOAmb0FxIE/uvuLOUlYH3U5Pyjwr1wFFoNjRoedSESqkkzA00Ng1rNw+LXQ+aywE9WJXHU584BOQDegDfCmme3h7iszG5nZEGAIQLt27XL00SE58MKgwL96NViMGEeT0pmlIvVLsjzosc96FrpfDQecG3aiOpNNNVoEZF6L2ya9LtNCoMTdy939c2AuQbH/CXcf4+7F7l7cqlWrDc1cf/zmYjjkD/DR49ySfwdxkmEnEpF1EmuDg6freuyN7GY82RT3aUAnM+toZgVAf6CkQptnCHrtmFlLgmGa+TnMWX8ddCkcOoK+8bcZlX8bBZSHnUhEyn+EcSfBnEnBGHsj6rGvU21xd/cEMBSYDMwGxrv7TDMbaWa9080mA9+Y2SzgdeD/3P2b2gpd7/zmEq4uP5We8WmMyb+FQsrCTiTSeJWtgkdPgNJXg9MdG8kYe0VZjbm7+yRgUoV1IzKeO3Bx+tEo3Zc8gtU04bq8e7g//0bOLL+EVWwSdiyRxmXN90FhX/BucKLDL08MO1FodA5fDj2ePJgfvYBb8u/koYLrGFR2WdiRRBqsGl/1vXoFPHI8LJ4RXEW+2zG1mK7+U3HPsZJUF9aWF3B7/m08VnANrDoUmrYMO5ZItH3/NTx0LKz4LLg13q7/+wUQ9fPZq6Jz92rB5NR+nFV+KTva1zC2ZzCdqIjUjuWlcG8P+G4hnPLkTwp7Y6biXkveSP2SU8uGw6qlwTfekplhRxKJnq8/gLE9grNjBj0PHQ8KO1G9oeJei6b5rnB6+kLdsUfAF2+FG0gkSua/AfcfBflN4YzJkbrRRi6ouNe2bYpg8EvBzXYfOgZmPxd2IpGGb+YzwR3SmreFwZOh5U5hJ6p3dEC1hjbo4MwW7YKexaMnBFfMHXkzFJ+R+3AijcCIK87nqrwH+cA7MXjBRXx37QfAB7p3QgUq7rUs85fBJpzDqPzbOOT5i+A//4Zuw4OZJUWkWjFSXJH3CIPzXuDl5L6cX34uP9Lkv9sb61kxVVFxr0M/0oQh5Rdzrd/LCW/cACs+hz5/h7zCsKOJ5Ext3JWsCWv5W/4oDo9P577E4fwpcaom6quG9k4dS5DHZYkhwYRjH4+HB/vAqsYzU4NITbXkO8YV/InDYu9xdfmpXJ0YqMKeBe2hUFgw4dhx98Gi9+GeQ2HZ3LBDidQ7O9oini4YwS62kN+WX8R9ySPCjtRgqLiHafdjYdBEWPsfuLd7cGqXiADQJfYxTxVcRRMr48SyP/ByqjjsSA2KinvY2u4HZ70Km28HDx8L7z8YdiKRcLnD1Dt4MP96FvtWHFM2ko98x7BTNTgq7vVBiw7BufAdD4KS82DS/wV3kBFpbMrXwDPnwOThvJwq5tiyq1noEbixTwh0tkxIKjujIM7pXJ63CUP+NYZ333mTc8suYDnNgY0700CkQfh+MTx+CiyaDt2G87sXd8PV/9xgKu71SJI41yZO5pNUB27Iv5uSwis4u+yi9f5JWhunnYnUuYXTYdzJwfGnEx6Cot74izpvfWPo12I9VJLqQr+yP5IixhMFIzkurgOtElHuwXGm+3oF13uc+TIU9a7+dVKtrIq7mfU0szlmVmpmwyrZPsjMlpnZjPTjzNxHbVxmeQeOXvtnpqd25i/5o2HipcENf0Wiomx1ML5ech60PwCGTIFtdgs7VWRUW9zNLA6MAo4AioABZlZUSdPH3X2v9OOeHOdslL6lGaeVD+PuRC+YdjeMPTy4qlWkoVs+L7i+48PHoOvlcMpTsOmWYaeKlGx67p2BUnef7+5lwDigT+3GknWSxLkmcQqc+Ah8Mx9Gd4VZJWHHEtlwnzwJY7rBD0uCm2sc/HuIxcNOFTnZFPfWwIKM5YXpdRX1M7OPzGyCmbXNSTr5n18cBWe/CVvtCONPhRcu1zCNNCyJtcHw4oQzguGX3/4Ddjo07FSRlasDqs8BHdx9T+Bl4IHKGpnZEDObbmbTly1blqOPbkRadAimDv7VOfDuXXBvD9rakrBTiVRrR1sE93QPhhcPGBpcmd28sj6i5Eo2p0IuAjJ74m3S6/7L3TNnvroHuLGyN3L3McAYgOLiYq9R0kbup6c8HkiPWCE3fT2aiQW/56ryQTydOhDQ9MFS3zgnx1/lyryH4bvNoP+jusdpHcmmuE8DOplZR4Ki3h84KbOBmW3n7ovTi72B2TlNKT/zUmo/ZpV14Jb8O7i14E4OTb7PleVnsJLNw44mjUR186dvyffckD+Gw+Lv82ZyDw46ZwJsvm0dpZNqh2XcPQEMBSYTFO3x7j7TzEaa2boTUs83s5lm9iFwPjCotgLL/yz0VvQv+wM3lPenR2w6kwsv5zexj8KOJULX2IdMLrycg2IfcXX5qQwsv1yFvY5ldYWqu08CJlVYNyLj+XBgeG6jSTZSxLgz2Zs3U3tya/4oHiq4nvsTPaDsYCjYNOx40shsyhouz3uMgXkvMyfVhlPLh/Optws7VqOk6QciYqZ34Oiya7gs73EG571A6Z/35bLyIbzvO/+knaYlkNry69gn3JB3N61tOWMTPbkh0Z+1FIQdq9HS9AMRspYC/pQ4lZPLhtPEyphQcDVX5T3ApqwJO5pE2Oas5tq8e3i04FrKyOOEsj8wMnGaCnvIVNwj6K3UHhy+9gYeSPZgYPwlXiq8jK6xD8OOJRHULTaDyYWXcWL8de5KHEWvsuuY7ruGHUvQsExkrWITrk4M5LnkAdyYP4YHCm7gyeSBsPoAXeYtG++Hpdycfwf94v9kbqo155RfyAzfKexUkkE994h733emV9l1/C1xDL1jU+Hv+8EHj0AqFXY0aYhSSfjX3XB7MUfHpnJboi9HlV2rwl4PqefeCJSRz62J45mU3J/JbZ+CZ8+B9+6DXjfB9ntX+3rNGS9AMOf6xIth8YewQzd6zj6K+b591i/X91HdUs+9EZnj7YLpC/reBd9+CWMOhucuhNUrwo4m9dnqFfDcBcH0AT8shePGwqnP1KiwS91Tz72xicVgrwGway+Ycj28OxpmPQOH/AH2GQjxjf+WUA8tIhJljLzqYs7Pe4rN+JH7kkfw12X9WPVwEypc9iL1kHrujVWT5tDzOjj7n7D1bsGf23d1gTkvBHfHkcbLHWY+A6M6MyL/IT5K7cCRZddyTeIUVrFJ2OkkSyrujd02RTDoeTjhQUiWw2P9g1ueLZgWdjIJw4J/wb094ImBkL8Jp5Vdzmnlw4MhPWlQNCwjYAZFfdjpQegff50LvnyKVvd2Z1KyMzclTgS2Czuh5EiVQ2YX7wBTroNZz8Jm20Dv22Gvk3nz9y+Glkk2joq7/FeCPB5OHsZTyd9wZnwSQ/Kep0fBdJ5K/oZRyT586T+d+Ek/lA3fTraQC/KegjvehYKm0HUY/Po8KNws7GiykVTc5WdW04TbksfyaPJQzsl7lpPir9Iv/ibPprowKtGHzzzaN1lY3y+tqBwU3tEWcX7e0xwdm8pqCuHAi4KirgvcIkPFXaq0nOaMTJzGnYnenJU3kVPir9C34C2eT/2Kvyf6Mtd1N8WGZjf7giF5z3N0bCo/UsDo5FGMSRzJB90HhB1NckzFXaq1jC24NnEydyWO5sy8SZwWf4nehVN5Jbk39yZ7MTVVhO4ClRu1chqpO5S+Am/fzsTCN/jBmzAmeSRjEkexgmYb/r5Sr6m4S9ZW0IwbE/0ZkziSQXmTOSX+Co/Fr2F2qh1jkz0pSf461JkAdX59BYm18PET8PbfYdls2Hx7risfwGPJQ/iepmGnk1qm4i41tpLN+WviOO5M9KZ3/G3OiL/ATfljuDxvHI8ku/NI4lCW0iLsmI3Xt1/A+w/C+w/BqqWwze5wzGjY7VhGX/lypS/RwfHoyaq4m1lP4G9AHLjH3a+vol0/YAKwn7tPz1lKqZfWUsATyW48kezKAbFZDI5P4oK8pxgaf5rXU3sxPtmN11J7k1AfovYly4ML0N67Hz57LTi9tdPhsP8Q2OHgYFkalWp/6swsDowCDgMWAtPMrMTdZ1VotzlwAfBubQSV+syYmtqNqandaJ/4NyfGp9Av/ibd4x+wzJvzZPIgxie7hh2yXslZT3nZXPhoHHzwMPywBJq1hm7DYO9ToHmb3HyGNEjZdKk6A6XuPh/AzMYBfYBZFdr9CbgB+L+cJpQG5UvflhsT/bk5cTzdYjM4MT6FM+MTOTvvObjncdi9HxT1hWa6MGqDff81fPJkMJ6++EOwGHTqAfueDjt1z8n8QNLwZfNd0BpYkLG8ENg/s4GZ7QO0dfeJZqbiXo/V1dhqkjivpvbl1dS+tOJbjo3/k+GJT+DFYfDicGjfBXY/BnY5UoU+Gyu/gk8nwuzn4cu3AIft94HDr4Pdj4XNt632LaRx2ehf8WYWA24BBmXRdggwBKBdO81V0VgsowWjk0cz+ouj2dEWcVTsHY7+fCo7fXkJTLwkmFN+l16wyxHBwT+ND2Ok2M2+oFvsQw6PT4O/fhFs2LoIul4OexwPLXWDDKlaNsV9EZB5tUqb9Lp1Ngd2B6ZY8EO5LVBiZr0rHlR19zHAGIDi4mJNPdgIfeat+VuyH39LHksnW0T32PtcHpsPr18Lr18Dm28HO3QLHh27Vtqrj+qZHa34lgNjn9A1/iEHxj6hpX0PwHupTnDYSNj1KNhqx5BTSkORTXGfBnQys44ERb0/cNK6je7+HdBy3bKZTQEu1dkysn7GPG/DvGQbLj/zyOAmEHMnB2d6zHsJPnwsaNZyF2h/ALTdH9p0jlBxc9rbEjrHPmU/m8N+sU/pGFsCwHJvxj9Se/Bmck/+kdqT5TTniy6Vn6uvc/ulKtUWd3dPmNlQYDLBqZBj3X2mmY0Eprt7SW2HlEZgs61hn1ODRyoFSz6B+VPg8zfgk6eDU/wANt2KsfltmekdmJVqzyxvz1e+Nb6e2atzWQA36L1SSVj5JSz+CBbPgK9n8EHhNFrYDwB865sxLbULj5R3Z2qqiFnefr3/n43JKY1HVmPu7j6JCrdecfcRVbTttvGxpFGLxWC7PYNHl/ODYr98Lix4Fxa8y/bvv8lBsY/Iywtu8v2DN+Ez357PfVs+T23H574dX/g2LPatWE6zjS6UWUml4Id/w8oF8N0CWDEfls2B5XNg+TxIrEn/3/Jg6yJeTO7Hx74D01K7UOrb101GaVR0zpTUf7EYbL1r8Nh3ID3fmUghZXSyhRTFvqTIvmQHW0xxbC69Y1OJ2f8O55R7nCW0YLFvyTfenJXelJVsxkrfjN9f8SprPJ81FLCWfNZSQJIYv445MZwYKWKk2IQymtoaNmEtTVlDM1vNlnzPVhY8WvId/HkgpMp/mrt5O15f0YJSP4RSb83sVDvmeFvWfhHeFA3SeKi4S4O0lgI+8R34JLnDT9YXUkY7W0p7W8K2toLt7Bu2tRVszwo62L/ZIvYDLfiBQiuv4p2rV+5xvqEZK7wZy70ZX7E17Q7oDFu0hS3aQ/O2wfOCppyu4REJiYq7RMpaCoIDtb6+qzOdJpTRnFUUWjmFlNOEMppQRgwnhZHCcIwUMX6kgFU04UcvZBVNWEMBFWfB7HtY7R7A1Bi61JSKuzRCxhoKWUMh6IRciSgdxRERiSD13EVyQMMmUt+o5y4iEkEq7iIiEaTiLiISQSruIiIRpAOqEjodjBTJPfXcRUQiSMVdRCSCVNxFRCJIxV1EJIJU3EVEIkjFXUQkglTcRUQiKKvibmY9zWyOmZWa2bBKtp9tZh+b2Qwz+6eZFeU+qoiIZKva4m5mcWAUcARQBAyopHg/6u57uPtewI3ALTlPKiIiWcum594ZKHX3+e5eBowD+mQ2cPfvMxabolsgiIiEKpvpB1oDCzKWFwL7V2xkZucCFwMFwCGVvZGZDQGGALRr166mWUVEJEs5O6Dq7qPcfUfgcuDKKtqMcfdidy9u1apVrj5aREQqyKa4LwLaZiy3Sa+ryjig78aEEhGRjZNNcZ8GdDKzjmZWAPQHSjIbmFmnjMUjgXm5iygiIjVV7Zi7uyfMbCgwGYgDY919ppmNBKa7ewkw1My6A+XAt8DA2gwtIiLrl9V87u4+CZhUYd2IjOcX5DiXiIhsBF2hKiISQSruIiIRpOIuIhJBKu4iIhGk4i4iEkEq7iIiEaTiLiISQSruIiIRpOIuIhJBKu4iIhGk4i4iEkEq7iIiEaTiLiISQSruIiIRpOIuIhJBKu4iIhGk4i4iEkFZFXcz62lmc8ys1MyGVbL9YjObZWYfmdmrZtY+91FFRCRb1RZ3M4sDo4AjgCJggJkVVWj2AVDs7nsCE4Abcx1URESyl03PvTNQ6u7z3b0MGAf0yWzg7q+7++r04jtAm9zGFBGRmsimuLcGFmQsL0yvq8pg4IXKNpjZEDObbmbTly1bln1KERGpkZweUDWzU4Bi4KbKtrv7GHcvdvfiVq1a5fKjRUQkQ14WbRYBbTOW26TX/YSZdQeuALq6+9rcxBMRkQ2RTc99GtDJzDqaWQHQHyjJbGBmewOjgd7uvjT3MUVEpCaqLe7ungCGApOB2cB4d59pZiPNrHe62U3AZsATZjbDzEqqeDsREakD2QzL4O6TgEkV1o3IeN49x7lERGQj6ApVEZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCMqquJtZTzObY2alZjasku0Hmdn7ZpYws+NyH1NERGqi2uJuZnFgFHAEUAQMMLOiCs2+AgYBj+Y6oIiI1Fw291DtDJS6+3wAMxsH9AFmrWvg7l+kt6VqIaOIiNRQNsMyrYEFGcsL0+tERKSeqtMDqmY2xMymm9n0ZcuW1eVHi4g0KtkU90VA24zlNul1NebuY9y92N2LW7VqtSFvISIiWcimuE8DOplZRzMrAPoDJbUbS0RENka1xd3dE8BQYDIwGxjv7jPNbKSZ9QYws/3MbCFwPDDazGbWZmgREVm/bM6Wwd0nAZMqrBuR8XwawXCNiIjUA7pCVUQkglTcRUQiSMVdRCSCVNxFRCJIxV1EJIJU3EVEIkjFXUQkglTcRUQiSMVdRCSCVNxFRCJIxV1EJIJU3EVEIkjFXUQkglTcRUQiSMVdRCSCVNxFRCJIxV1EJIKyKu5m1tPM5phZqZkNq2R7oZk9nt7+rpl1yHVQERHJXrXF3cziwCjgCKAIGGBmRRWaDQa+dfedgFuBG3IdVEREspdNz70zUOru8929DBgH9KnQpg/wQPr5BOBQM7PcxRQRkZrIpri3BhZkLC9Mr6u0jbsngO+ArXIRUEREai6vLj/MzIYAQ9KLP5jZnA18q5bA8tykyinlqhnlqrn6mk25asBu2Khc7bNplE1xXwS0zVhuk15XWZuFZpYHNAe+qfhG7j4GGJNNsPUxs+nuXryx75NrylUzylVz9TWbctVMXeTKZlhmGtDJzDqaWQHQHyip0KYEGJh+fhzwmrt77mKKiEhNVNtzd/eEmQ0FJgNxYKy7zzSzkcB0dy8B7gUeMrNSYAXBLwAREQlJVmPu7j4JmFRh3YiM52uA43Mbbb02eminlihXzShXzdXXbMqqPHZPAAAFXUlEQVRVM7WeyzR6IiISPZp+QEQkguptcTez481sppmlzKzKo8pVTY2QPgD8bnr94+mDwbnItaWZvWxm89L/tqikzcFmNiPjscbM+qa33W9mn2ds26uucqXbJTM+uyRjfZj7ay8zm5r+en9kZidmbMvp/tqYqTTMbHh6/RwzO3xjcmxArovNbFZ6/7xqZu0ztlX6Na2jXIPMbFnG55+ZsW1g+us+z8wGVnxtLee6NSPTXDNbmbGtNvfXWDNbamafVLHdzOy2dO6PzGyfjG253V/uXi8fwC+AXYApQHEVbeLAZ8AOQAHwIVCU3jYe6J9+fhfwuxzluhEYln4+DLihmvZbEhxk3jS9fD9wXC3sr6xyAT9UsT60/QXsDHRKP98eWAxskev9tb7vl4w25wB3pZ/3Bx5PPy9Kty8EOqbfJ16HuQ7O+B763bpc6/ua1lGuQcDfK3ntlsD89L8t0s9b1FWuCu3PIzgRpFb3V/q9DwL2AT6pYnsv4AXAgF8B79bW/qq3PXd3n+3u1V3kVOnUCGZmwCEEUyFAMDVC3xxFy5xqIZv3PQ54wd1X5+jzq1LTXP8V9v5y97nuPi/9/GtgKdAqR5+faWOm0ugDjHP3te7+OVCafr86yeXur2d8D71DcL1Jbctmf1XlcOBld1/h7t8CLwM9Q8o1AHgsR5+9Xu7+JkFnrip9gAc98A6whZltRy3sr3pb3LNU1dQIWwErPZgKIXN9Lmzj7ovTz/8NbFNN+/78/BvrmvSfZLeaWWEd52piZtPN7J11Q0XUo/1lZp0JemOfZazO1f7amKk0snltbebKNJig97dOZV/TuszVL/31mWBm6y54rBf7Kz181RF4LWN1be2vbFSVPef7q06nH6jIzF4Btq1k0xXu/mxd51lnfbkyF9zdzazK043Sv5H3ILhGYJ3hBEWugOB0qMuBkXWYq727LzKzHYDXzOxjggK2wXK8vx4CBrp7Kr16g/dXFJnZKUAx0DVj9c++pu7+WeXvkHPPAY+5+1oz+y3BXz2H1NFnZ6M/MMHdkxnrwtxfdSbU4u7u3TfyLaqaGuEbgj938tK9r8qmTNigXGa2xMy2c/fF6WK0dD1vdQLwtLuXZ7z3ul7sWjO7D7i0LnO5+6L0v/PNbAqwN/AkIe8vM2sGTCT4xf5Oxntv8P6qxMZMpZHNa2szF2bWneAXZld3X7tufRVf01wUq2pzuXvmNCP3EBxjWffabhVeOyUHmbLKlaE/cG7milrcX9moKnvO91dDH5apdGoED45QvE4w3g3B1Ai5+ksgc6qF6t73Z2N96QK3bpy7L1DpUfXayGVmLdYNa5hZS6ALMCvs/ZX+2j1NMBY5ocK2XO6vjZlKowTob8HZNB2BTsC/NiJLjXKZ2d7AaKC3uy/NWF/p17QOc22XsdgbmJ1+Phnokc7XAujBT/+CrdVc6Wy7EhycnJqxrjb3VzZKgNPSZ838Cvgu3YHJ/f7K9dHiXD2AYwjGndYCS4DJ6fXbA5My2vUC5hL85r0iY/0OBD98pcATQGGOcm0FvArMA14BtkyvLwbuyWjXgeC3cazC618DPiYoUg8Dm9VVLuDX6c/+MP3v4Pqwv4BTgHJgRsZjr9rYX5V9vxAM8/ROP2+S/v+XpvfHDhmvvSL9ujnAETn+fq8u1yvpn4N1+6ekuq9pHeW6DpiZ/vzXgV0zXntGej+WAqfXZa708h+B6yu8rrb312MEZ3uVE9SvwcDZwNnp7UZw86PP0p9fnPHanO4vXaEqIhJBDX1YRkREKqHiLiISQSruIiIRpOIuIhJBKu4iIhGk4i4iEkEq7iIiEaTiLiISQf8PvBZqMsv7OpQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi41LCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvSM8oowAAIABJREFUeJzt3Xd8VHX2//HXySQBxUa10FUsrPK1RKxrb6ACCiiorCiKjVVBV9F1LfizF77uLhbsHSmWLKAoirqy4oIVAZGIhWAhotgQMpM5vz/u4HeMCZnAJHdy834+HvPI3Hs/M/fNnXBy53Pv/Vxzd0REJFrywg4gIiLZp+IuIhJBKu4iIhGk4i4iEkEq7iIiEaTiLiISQSruIiIRpOIuIhJBKu4iIhGUH9aKW7Vq5Z06dQpr9SIiDdJbb731jbu3rqldaMW9U6dOzJkzJ6zVi4g0SGb2WSbt1C0jIhJBKu4iIhGUUXE3syPNbKGZlZjZyCqWdzCzGWb2jpm9b2Y9sx9VREQyVWNxN7MYMAboAXQFBppZ10rNLgfGu/uuwADgjmwHFRGRzGWy594dKHH3xe5eDowDeldq48AmqeebAl9kL6KIiNRWJmfLtAWWpE2XAntWanMV8IKZ/RloBhyalXQiIrJOsnVAdSDwoLu3A3oCj5jZ797bzIaa2Rwzm1NWVpalVYuISGWZFPelQPu06XapeemGAOMB3P0NoCnQqvIbuftYdy9y96LWrWs8B19ERNZRJsV9NtDFzDqbWSHBAdPiSm0+Bw4BMLMdCYq7ds1FREJSY5+7uyfMbBgwDYgB97v7PDMbBcxx92LgQuAeMxtOcHB1sOvO2yLSCHQaOaXaZZ/ecFQ9JvmtjIYfcPepwNRK865Iez4f2De70UREZF3pClURkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhqeMV92Yfw/GWQTIadREQkZzW84r54BswaAzNHh51ERCRnNbzivudZsFNfeOkaKHkp7DQiIpn56AWYcV299TpkNPxATjGDXv+AZQtg0hAY+go07wRUP8ZDmOM7iIiw/GN46nTYrAPsNxzyNqjzVTa8PXeAwmZwwqPgSXjyZChfGXYiEZHf6TRyCjuOnMSC23vz3S8V7PfZEDr97eV6WXfDLO4ALbeB4+6Frz6AycNBg1CKSM5xbii4l+2tlPPiwyj1NvW25obXLcNvu1/Oi/VlxPvjuOKtJsAR4YUSEalkSOw5esf+w03xE/h3slu9rrvh7rmn/KOiD9MrduVv+Y9SZB+GHUdEBID98uZyWf5jPF+xB3dU9Kr39Tf44u7kMSJ+Dku8NXcU/p3N+TbsSCLSyHW0r/hnwd9Z5O0YET8bsHrP0OCLO8APNOPM+Ag2ZBVjC2+jCeVhRxKRRmojVnJvwa04xhnxEaykaSg5IlHcARZ5Oy6In8vO9gk3F9xNcLc/EZH6YyQZXXAHne1LzomfzxLfPLQsGRV3MzvSzBaaWYmZjaxi+Wgzezf1+MjMVmQ/as2mJ3fnlsTx9Iq9wTmxyvfwFhGpW8PzJ3JY7G1GJQbxRvIPoWap8WwZM4sBY4DDgFJgtpkVp+6bCoC7D09r/2dg1zrImpE7Knqxfd4SLsofz0fejunJ3cOKIiKNSM+8WZyX/wzjEgfycMXhYcfJaM+9O1Di7ovdvRwYB/ReS/uBwBPZCLdujIvjQ5nrnfnfgjFsb5+HF0VEGoWu9im3FNzNnOR2XJE4lTAOoFaWSXFvCyxJmy5NzfsdM+sIdAaqvATLzIaa2Rwzm1NWVlbbrBlbTSFDy0fwM025t+BW+Hl5na1LRBq31qzgnsJbWUEzziofTjkFYUcCsn9AdQAw0d0rqlro7mPdvcjdi1q3bp3lVf/W17RgaPkI2tgKGP8nSOgMGhHJsvKV3FN4C835iTPKL+IbNg070a8yKe5LgfZp0+1S86oygFC7ZH7rPd+Wi+NnwGeva4gCEcmuZBKePpNu9gnnxYcxzzuFneg3Minus4EuZtbZzAoJCvjvTkUxsx2A5sAb2Y24fp5N7gf7XwzvPgr/vjXsOCISFS+PggXFXJs4KSdP3KixuLt7AhgGTAMWAOPdfZ6ZjTKz9GtqBwDj3HNw9/igy2Dn/vDyNTB3YthpRKShe/sReH00FJ3GfRU9wk5TpYwGDnP3qcDUSvOuqDR9VfZiZZkZ9B4D35fCM+fApu2gw15hpxKRhuiT12DyBbD1QdDjJnj9hbATVSkyV6jWKL8JDHg8KOxPDAwGzxcRqY1vFgX3kGi5LfR/EGK5cWZMVRpPcQfYsAWcNCF4/vjxsFKDjIlIhn4qg8f6QV4BnPgkbLBZ2InWqnEVdwhu8jHgcVjxefAXOLE67EQikutW/wSP94cfv4aB4369tWcua3zFHaDj3tDnTvhsJjxzdr3dsFZEGqCKOEwYDF++B/0fgPZ7hJ0oIw3yTkxZsXM/+GEpvHgFbLQ5HHFdcOBVRGQN9+DgacmLcPT/wva5eWZMVRpvcQfY5zz48SuYdQdsvAXse37YiUQkl7xyPbzzKOx/MZ0mtoGJU2p+TY5onN0ya5jB4dfCTn2DPfj3xoWdSERyxZwH4NUbYZeTg2tlGpjGvecOkJcX9L///A08ey5s2Aq6HBp2KhEJ08LnYcoI2PYwOOZ/G2SXbePec18jvwmc8Ci02RHGD4LSt8JOJCJh+ew/wQHULbrl/Lnsa6PivkbTTeCkSdCsdXDK0zclYScSkfr25Xvw+AmwaVs4aSI02SjsROtM3TLpNt4cBj0N9x0OD/eG056HzdrX/DoRafi+KYFHjmPpqgL6f38eX/y//4adaL1oz72yltsEBX71j/Bwr+CiBRGJtu9Lgx064OTyy/iCViEHWn+NYs+908jqT1/69Iajfj9zy27BMAWP9IFHjoXBk4OhC0Qken4qg4f7wOofYPBkPrm9NOxEWaE99+p02DMYpmD5omA8idU/hp1IRLJt1ffw6HHBnvuJ42HL/wk7UdaouK/NNgcFR8u/eDcYSTL+S9iJRCRLdhw5iTevO5z4l/MYvPLPdLrz27V+y29oGkW3zNpU92H+2l2zw1Fw7F3w1FAYfwoMeKzBnholIinxX7i34FaKbCHnx4fxSnKXsBNlXaMv7hnpdjyU/xTch3XiadDvfhV4kRxT447aGvFVMO5E9s6bz4j42UxO7l0P6epfRt0yZnakmS00sxIzG1lNm+PNbL6ZzTOzx7MbMwcUnQZHXA8LimHS6VCRCDuRiNRWYnVwoeLHL3NJ4gyeSe4XdqI6U+Oeu5nFgDHAYUApMNvMit19flqbLsClwL7u/p2ZtamrwKHa+xzwJLzwV7A8OO4eiOnLj0iDkCgPrjxd9AIcczsTJrQOO1GdyqQydQdK3H0xgJmNA3oD89PanAGMcffvANx9WbaD1rfqv+INCwr8i38Lxps4dqwKvEiuq4jDpNNg4VToeQvsPhgmROfgaVUyqUptgSVp06XAnpXabAdgZjOBGHCVuz+flYS5aN/zggI//cpgD/7YuyEvFnYqEalKRQKeOgMW/AuOvAG6nxF2onqRrV3OfKALcCDQDnjNzHZ29xXpjcxsKDAUoEOHDlladUj2uyAo8C9dHRT4PneqwIvkmAISMPHU4FjZYdfAXmeHHaneZFLclwLpA6y0S81LVwq86e5x4BMz+4ig2M9Ob+TuY4GxAEVFRb6uoXPGH0cEBf7la4Kffe5SF41IjmhCOWMKbocF7wQnQ+x9TtiR6lUmZ8vMBrqYWWczKwQGAMWV2jxDsNeOmbUi6KZZnMWcuWv/i+CQK2DuBJhwim64LZIDmrKaewpu5dDYO3DUbY2usEMGxd3dE8AwYBqwABjv7vPMbJSZ9Uo1mwYsN7P5wAzgL+6+vK5C55w/Xhj05X04GcadqCtZRUK0Iat4oOBm9sv7gL/Eh8IeQ8KOFIqM+hDcfSowtdK8K9KeOzAi9Wic9jobCjaEf50Pj/WHgU9Ak43DTiXSqGzMSh4ovIldrIQL4udQnNyXm8MOFRJ1EGfT7qcEBf7pM4PRJE+aABs0r7JpxlfTiUhGNuUnHiq8gT/YZwyLn8fzye5hRwqVinst1ViUu/WHgqYw4VR46BgY9Aw0a/hjQ4vktO+XMqHwajraMs6KX8BLyd3DThQ6jQpZF3Y8BgaOg28Wwf1HworPw04kEl3fLIL7j2BL+5ZT4peosKeouNeVLocGd3T6eVlw276v54WdSCR6lr4N9x8BiVUMKP8bs5Jdw06UM1Tc61LHfeDU1IW69/eAT2eGm0ckSj6eEXR9FjaD06YxzzuFnSinqLjXtc27wpAXYKM2wUHWBf8KO5FIwzfvGXj8eNisA5z2QnDvY/kNHVCtY2sOwDbnQu4vvIVu4wbxt8RpwCHhBhNpqN68G567BNrvCSeOq/aMtDWidHel2lBxz5KafoG+YxNOLL+MMQV/57qC+9jcvmN0oi9g9RNQpKFLVsALl8OsO3ihYnfOW3Qmq67+T9ipcpaKez36haYMjY/gOr+P8/OfooN9zSXxoZSz7nd10vny0iiU/wyTzoCFU7gv0YNrEyeRVK/yWqm417ME+VycGMpnvjl/KRhPW/uGM8uH8x2bhB1NJDf9+DU8cQJ8+R70uIlrnm4XdqIGQX/6QmGMqejDsPI/8z+2mKcLr2QbqzzQpoiw7EO491AoWwgnPAZ7nhl2ogZDxT1Ek5N7M7D8r2xkv/BU4ZXsnadz4UV+9fHLwTUiFavh1KmwQ8+wEzUoKu4he9u3o0/5KL725jxccAO8/XDYkUTC5Q5vjIFH+8Km7eD06bDVrmGnanBU3HNAqbehb/nVvJHsCsV/hql/Ce75KNLYxFfBM2fDtMtgh6OCa0Q2a+B3bQuJDqjmiB/ZkFPjF/Pxfm/CG/+Erz6A4x8KLn4SibA1Z3y14TvGFt7GLnkfc2u8H/98pw/+zqshp2u4VNxzSAUxOOJa2HKXYA9+7IFwwiPQVgMhScOxLqfn7mIl3F14G81YxdDy4byQ3KOu4jUa6pbJRd36w5BpYLFgTJp3Hgs7kUjdcOeE2AyeLBzFKi/kuPKrVdizJKPibmZHmtlCMysxs5FVLB9sZmVm9m7qcXr2ozYyW/4PDH0FOuwJz54DUy7S/VklWsp/hmfO5saCe3gzuSO9y6/hI28fdqrIqLFbxsxiwBjgMKAUmG1mxe4+v1LTJ919WB1kbLyatYSTn4bpVwb98EvnQL8HoEXnsJOJrJ+yj2D8n6DsQ0bH+/KPimN1xWmWZbI1uwMl7r7Y3cuBcUDvuo0lv4rlB/3wJzwKyxfD3QfA/OKwU4msu7kT4Z6DgnsdDHqK2yv6qrDXgUy2aFtgSdp0aWpeZX3N7H0zm2hm+m6VbTseA2e9FgxtOn5QMCqeummkIUmshikXwqQhsPlOcNbrsM3BYaeKrGz9ufwX0MnduwEvAg9V1cjMhprZHDObU1ZWlqVVNyLNO8Fp02Cvc+DNu+C+w2lvX4edSqRG21op3HMIzL4X9h4GgyfDJluFHSvSMjkVcimQvifeLjXvV+6+PG3yXuCmqt7I3ccCYwGKioq8VkklkF8IR14PHfeFZ89hSuFlXBkfzNPJ/dDwwZJ7nJNj07k8/1GWf9WUS+IXMn3G7jDjhbCDRV4mxX020MXMOhMU9QHAiekNzGxLd/8yNdkLWJDVlPJ7Ox4NW+zMh6P7M7rwTg6peJvL46exgo0BDQUs4WvJ99xUMJZDYu/wakU3LoqfRRmbhR2r0aixuLt7wsyGAdOAGHC/u88zs1HAHHcvBs4zs15AAvgWGFyHmWWN5h0ZUP43zoxNZnj+BPZospCL4mfx72S3sJNJI3dg3rvcXHAXm/ALV8X/xEMVh+M6aFqvMrpC1d2nAlMrzbsi7fmlwKXZjSaZSJLHnRW9eC3ZjdEFY3ik8AYeTBzODYmBrKJJ2PGksVn9I6PyH+BP+S/yYbI9J8cvY6FrbJgwaPiBiJjnnTim/Fouzn+SIfnPsV/eB1wcH8rbvl3Y0SSiKnf97Zs3lxsL7uHk2HLuTfTg5sQJrKYwpHSi70kRsppCrkkM4qTyS2lq5UwsvJor8x9iQ1aFHU0ibGNWcn3+PTxWeD2rvYD+5Vfw/xKDVNhDpuIeQTOTO3PE6ht5uOIwTom9wAtNLuaAvPfCjiURdGDeu0xrcjHHx17hrsQx9Cy/nrd8+7BjCeqWiayf2YCrEoMprtiHmwrG8lDhjUyq2A9W7g0btgg7njR0Py3jtoI7OC72OguT7TgrPpz3fZuwU0ka7blH3Nu+HT3Lr+f2xLH0ynsD/rlHMMpkMhl2NGmIkhXw5lj4RxFH5c3i9sSxHFN+rQp7DtKeewNR3XnrmSingNGJ/jxXsSfPt38qGGXyrQeh582w1S7ZCyk5b72uf1gyG6aMgK/eh60P4sgFR/OJb5nlhJIt2nNvRD70DsHwBX3uhO8+CW4GMnk4rPw27GiSy1Z+C8XnwX2Hws9l0P9BGPS0CnuO0557I9PpsueATdmY67kgfxKnzH6QH2eP55bE8Vx7zS2QFws7otRCnV6JnFgdjAXz6k2w+sdgTJgDR0KTjdf/vaXOac+9kfqRDbkmMYie5dez0NtzbcH9cOc+sPD54O7z0ni5w7ynYUz34EbVbXcPRnA84loV9gZExb2R+8jbM6D8cs4qvwAq4vDECfDgUVA6J+xoEoYl/4X7DocJg6GgGZw8CQY9BZt3DTuZ1JK6ZXLM+hw4XXfG88nucO7l8PZD8MqNcO8hsGMvOORKOt2ysMpXaRCy6Dj80ru4IH8SPWP/5WvfjFsSQ5n0+f4s3vbQsKPJOlJxl/8TK4A9ToduA+CNMfCfv8OHU7gx/4/cUdGLz3yLsBNKhjLdSdjWSrkg/yl65r3JzzRldLwvYyuO4hea1nFCqWsq7vJ7TTaCAy+BolPh37fSe9Z99Iu9yjPJfRmT6MNij8ZNFhrzsMjb2FLOz3+Ko/NmsZIm3FHRi3sSR/E9G4UdTbJExV2qt1Eb6HEjf3y1G2fkT+Hk2HSOLZzJ5ORe/CNxbNjpZB38wT7hzPzJHJ03i18o5K6KY7gn0ZPv2CTsaJJlKu5SozI247rESdydOJrT86cyKPYivZq8AY/PgL3Ohs77g637XaAa8x50/XAOzHuPM2JT2Dc2j5+8KWMrjmZs4ii+VVGPLBV3ydhyNuXGxEDuThzNqfnTOL/0NXj4ueBmx3udDTv1g4LG2Vebi3+gConTOzaTM2JT2C5vKV96C66LD2RcxcH8QLPQckn9UHGXX2V6EG4FGzM60Y/zh98JcyfArDvh2XNh+lVQNAR2Hwyb6OrFsLSzZQyMvczxsVdpbd8zP9mR4eVnMzm5N3H9l280MvqkzexI4HaC2+zd6+43VNOuLzAR2MPddaJ01BU0hd0Gwa4nwyevBkX+1RvgtZuhy+HBsi6HB2fhSJ3KJ8GheW9zYuwl9o/NpcKNl5O78mDFEcxM7oRunt741FjczSwGjAEOA0qB2WZW7O7zK7XbGDgfeLMugkoOM4OtDwweyz+Gdx6Bdx+Hj56DZm1gl4Gw65+g1ba1ettc7OrINdvYUo6L/Zv+sddoYytY6i25Ld6P8RUH8BUtw44nIcpkz707UOLuiwHMbBzQG5hfqd01wI3AX7KaUBqWltvAoVfBQZdDyYvw9sPwn3/CzNuh3R6wU1/o2kfdNuvjhy9g7kSYO4GXmrz/6176ExUH80pyF5K68FzIrLi3BZakTZcCe6Y3MLPdgPbuPsXMVNwFYvmwfY/g8eNX8N4TMHcSPD8Snr8UOu4LOx0L22svPCMrPocFk+HDyfDZfwCHrXbj6vggJlfsTRmb1clqw7liWrJhvY+umFkecBswOIO2Q4GhAB066I7ojcbGW8B+w4NH2Ucw7yn4YBJMuRCmXMizhVszvWI3pid3Z4F3QP3DBDdT+eo9WPQiLPhXMIY68GGyPc8nj+XZin35ZLG+/Uj1MinuS4H2adPtUvPW2BjYCXjFgnOdtwCKzaxX5YOq7j4WGAtQVFSkoQcbo9bbBcPGHnAJlH0IC6eSePEJhudP4kKbyJfegpnJnZhZ8QdmJndiGc3DTrxeanXc4Mev4OMZ8PFLwc+V3wTz23XnuvhApiX30BAQkrFMivtsoIuZdSYo6gOAE9csdPfvgVZrps3sFeAinS3TeGVU0MygzY7QZkf6TtmBVnzPQbF3OCDvfQ7Oe5t+sdcAWJRsy+zk9ryV3I63vAufRqW4ucO3i+HzN+CzN+CzmcENVACatYZtD4FtDoFtDoKN2jBW3SNSSzUWd3dPmNkwYBrBqZD3u/s8MxsFzHH34roOKbkpm/2x37ApEyoOZELFgRhJutrn7JP3AfvmzePo2CxOzH8ZgOW+MTy2N2zRDbbYOXg07wx51R9EzGbOdXmvPJK0szL+YJ+yc94n7GyL4aZz4ZfvggYbtIAOe8MeQ6DTH4N/21r+PSKZyKjP3d2nAlMrzbuimrYHrn8sacycPOZ5J+ZVdOKeiqMxkmxrX7Bb3iKKbCH9vy+FkpfAK4IXFG4ErbpAy23/79GiM2zSLtgLrgdGkjasoK19Qzsro6N9zbZ5X9DFlrK1fUFTiwMQ9xgLvT3seAxstSt02AdabadiLlmny9Uk5zl5LPJ2LKpox5McxF8+hyaU08VK6Zr3GV0Tn7H1ki/Z/+c3g1MESTuck5fP600240tvwXLflBXejBVszApvxg80Y5UXsIpCVlPAagqpII8TL/uAPJw8kuSRZAPKaWar2IDVnBVbxSa2khb8QEsLHq34ni3sWwqt4je5lyRbU+JbMTP5B0q8LfOTHVno7SmngE976SwhqVsq7tIgraaQD3xrPqjY+td5nw4/CuKrgr7rbxcH54P/8AVvvvpftuJbOtlXNM/7kc34iSaWWOd1xz3GcjbhW9+Eb3wTPqcNXyZbUuqtWeqtKPVWlHprVtEkG/9UkXWi4i71pl7OmS5o+uuB2jUunF55vU5TytmUn2licZoQpynlNKWcPJwkRhLDMZLk8QuF/ExTfvEm/ExTVlGITteUXKfiLo2QsYomwZ61TsiViNJRHBGRCNKeu0SGLpUX+T/acxcRiSAVdxGRCFJxFxGJIBV3EZEI0gFVkRDo4K/UNe25i4hEkIq7iEgEqbiLiESQiruISASpuIuIRJCKu4hIBKm4i4hEUEbF3cyONLOFZlZiZiOrWH6Wmc01s3fN7HUz65r9qCIikqkai7uZxYAxQA+gKzCwiuL9uLvv7O67ADcBt2U9qYiIZCyTPffuQIm7L3b3cmAc0Du9gbv/kDbZDN0CQUQkVJkMP9AWWJI2XQrsWbmRmZ0LjAAKgYOreiMzGwoMBejQoUNts4qISIaydkDV3ce4+zbAJcDl1bQZ6+5F7l7UunXrbK1aREQqyaS4LwXap023S82rzjigz/qEEhGR9ZNJcZ8NdDGzzmZWCAwAitMbmFmXtMmjgEXZiygiIrVVY5+7uyfMbBgwDYgB97v7PDMbBcxx92JgmJkdCsSB74BT6jK0iIisXUbjubv7VGBqpXlXpD0/P8u5RERkPegKVRGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQjKqLib2ZFmttDMSsxsZBXLR5jZfDN738xeMrOO2Y8qIiKZqrG4m1kMGAP0ALoCA82sa6Vm7wBF7t4NmAjclO2gIiKSuUz23LsDJe6+2N3LgXFA7/QG7j7D3VemJmcB7bIbU0REaiOT4t4WWJI2XZqaV50hwHNVLTCzoWY2x8zmlJWVZZ5SRERqJasHVM3sZKAIuLmq5e4+1t2L3L2odevW2Vy1iIikyc+gzVKgfdp0u9S83zCzQ4G/Age4++rsxBMRkXWRyZ77bKCLmXU2s0JgAFCc3sDMdgXuBnq5+7LsxxQRkdqosbi7ewIYBkwDFgDj3X2emY0ys16pZjcDGwETzOxdMyuu5u1ERKQeZNItg7tPBaZWmndF2vNDs5xLRETWg65QFRGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiaCMiruZHWlmC82sxMxGVrF8fzN728wSZtYv+zFFRKQ2aizuZhYDxgA9gK7AQDPrWqnZ58Bg4PFsBxQRkdrL5B6q3YESd18MYGbjgN7A/DUN3P3T1LJkHWQUEZFayqRbpi2wJG26NDVPRERyVL0eUDWzoWY2x8zmlJWV1eeqRUQalUyK+1Kgfdp0u9S8WnP3se5e5O5FrVu3Xpe3EBGRDGRS3GcDXcyss5kVAgOA4rqNJSIi66PG4u7uCWAYMA1YAIx393lmNsrMegGY2R5mVgr0B+42s3l1GVpERNYuk7NlcPepwNRK865Iez6boLtGRERygK5QFRGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiSAVdxGRCFJxFxGJIBV3EZEIUnEXEYkgFXcRkQhScRcRiaCMiruZHWlmC82sxMxGVrG8iZk9mVr+ppl1ynZQERHJXI3F3cxiwBigB9AVGGhmXSs1GwJ85+7bAqOBG7MdVEREMpfJnnt3oMTdF7t7OTAO6F2pTW/godTzicAhZmbZiykiIrWRSXFvCyxJmy5NzauyjbsngO+BltkIKCIitZdfnyszs6HA0NTkT2a2cB3fqhXwTXZSZZVy1Y5y1V6uZlOuWrAb1ytXx0waZVLclwLt06bbpeZV1abUzPKBTYHlld/I3ccCYzMJtjZmNsfdi9b3fbJNuWpHuWovV7MpV+3UR65MumVmA13MrLOZFQIDgOJKbYqBU1LP+wEvu7tnL6aIiNRGjXvu7p4ws2HANCAG3O/u88xsFDDH3YuB+4BHzKwE+JbgD4CIiIQkoz53d58KTK0074q056uA/tmNtlbr3bVTR5SrdpTO+i1IAAAFYElEQVSr9nI1m3LVTp3nMvWeiIhEj4YfEBGJoJwt7mbW38zmmVnSzKo9qlzd0AipA8BvpuY/mToYnI1cLczsRTNblPrZvIo2B5nZu2mPVWbWJ7XsQTP7JG3ZLvWVK9WuIm3dxWnzw9xeu5jZG6nP+30zOyFtWVa31/oMpWFml6bmLzSzI9YnxzrkGmFm81Pb5yUz65i2rMrPtJ5yDTazsrT1n5627JTU577IzE6p/No6zjU6LdNHZrYibVldbq/7zWyZmX1QzXIzs7+ncr9vZrulLcvu9nL3nHwAOwLbA68ARdW0iQEfA1sDhcB7QNfUsvHAgNTzu4Czs5TrJmBk6vlI4MYa2rcgOMi8YWr6QaBfHWyvjHIBP1UzP7TtBWwHdEk93wr4Etgs29trbb8vaW3OAe5KPR8APJl63jXVvgnQOfU+sXrMdVDa79DZa3Kt7TOtp1yDgX9W8doWwOLUz+ap583rK1el9n8mOBGkTrdX6r33B3YDPqhmeU/gOcCAvYA362p75eyeu7svcPeaLnKqcmgEMzPgYIKhECAYGqFPlqKlD7WQyfv2A55z95VZWn91apvrV2FvL3f/yN0XpZ5/ASwDWmdp/enWZyiN3sA4d1/t7p8AJan3q5dc7j4j7XdoFsH1JnUtk+1VnSOAF939W3f/DngRODKkXAOBJ7K07rVy99cIduaq0xt42AOzgM3MbEvqYHvlbHHPUHVDI7QEVngwFEL6/GzY3N2/TD3/Cti8hvYD+P0v1rWpr2SjzaxJPedqamZzzGzWmq4icmh7mVl3gr2xj9NmZ2t7rc9QGpm8ti5zpRtCsPe3RlWfaX3m6pv6fCaa2ZoLHnNie6W6rzoDL6fNrqvtlYnqsmd9e9Xr8AOVmdl0YIsqFv3V3Z+t7zxrrC1X+oS7u5lVe7pR6i/yzgTXCKxxKUGRKyQ4HeoSYFQ95uro7kvNbGvgZTObS1DA1lmWt9cjwCnunkzNXuftFUVmdjJQBByQNvt3n6m7f1z1O2Tdv4An3H21mZ1J8K3n4HpadyYGABPdvSJtXpjbq96EWtzd/dD1fIvqhkZYTvB1Jz+191XVkAnrlMvMvjazLd39y1QxWraWtzoeeNrd42nvvWYvdrWZPQBcVJ+53H1p6udiM3sF2BWYRMjby8w2AaYQ/GGflfbe67y9qrA+Q2lk8tq6zIWZHUrwB/MAd1+9Zn41n2k2ilWNudw9fZiRewmOsax57YGVXvtKFjJllCvNAODc9Bl1uL0yUV32rG+vht4tU+XQCB4coZhB0N8NwdAI2fomkD7UQk3v+7u+vlSBW9PP3Qeo8qh6XeQys+ZrujXMrBWwLzA/7O2V+uyeJuiLnFhpWTa31/oMpVEMDLDgbJrOQBfgv+uRpVa5zGxX4G6gl7svS5tf5Wdaj7m2TJvsBSxIPZ8GHJ7K1xw4nN9+g63TXKlsOxAcnHwjbV5dbq9MFAN/Sp01sxfwfWoHJvvbK9tHi7P1AI4l6HdaDXwNTEvN3wqYmtauJ/ARwV/ev6bN35rgP18JMAFokqVcLYGXgEXAdKBFan4RcG9au04Ef43zKr3+ZWAuQZF6FNiovnIB+6TW/V7q55Bc2F7AyUAceDftsUtdbK+qfl8Iunl6pZ43Tf37S1LbY+u01/419bqFQI8s/77XlGt66v/Bmu1TXNNnWk+5rgfmpdY/A9gh7bWnpbZjCXBqfeZKTV8F3FDpdXW9vZ4gONsrTlC/hgBnAWellhvBzY8+Tq2/KO21Wd1eukJVRCSCGnq3jIiIVEHFXUQkglTcRUQiSMVdRCSCVNxFRCJIxV1EJIJU3EVEIkjFXUQkgv4/WKxdv0SokwoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plot histograms, see that we recover (1 + cos^2 theta)\n",
    "\n",
    "xs=np.arange(-1,1,0.01)\n",
    "plt.hist(mupluscth,bins=50,density=True)\n",
    "plt.plot(xs,(3/8.)*(1+xs**2))\n",
    "plt.show()\n",
    "\n",
    "xs=np.arange(-1,1,0.01)\n",
    "plt.hist(muminuscth,bins=50,density=True)\n",
    "plt.plot(xs,(3/8.)*(1+xs**2))\n",
    "plt.show()"
   ]
  },
  {
   "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.17"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}