time of flight difference

particleFlightTime.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import pylab, numpy, math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# initialize nice display of equations\n",
    "import sympy\n",
    "from sympy import init_printing\n",
    "init_printing()\n",
    "\n",
    "# from sympy import init_session\n",
    "# init_session()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from sympy import Symbol, sin, cos\n",
    "from sympy.solvers import solve "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# radius of curvature R of a charged particle with momentum q and charge q\n",
    "# in a magnetic field B is\n",
    "# \n",
    "# R = p / (q * B)\n",
    "#\n",
    "# with p = 1 GeV/c, q = unit charge, B = 1 Tesla one gets\n",
    "#\n",
    "# R [m] = 3.33564095198 * p [GeV/c] / (q [unit charges] * B [Tesla])\n",
    "\n",
    "# R = Symbol('R') # radius of curvature \n",
    "p = Symbol('p', real = True) # particle momentum\n",
    "q = Symbol('q', real = True) # particle charge\n",
    "B = Symbol('B', real = True) # magnetic field\n",
    "\n",
    "# radius of curvature for particles with unit charge\n",
    "R = 3.33564095198 * p / B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "rEcal = 1.3 # approximate radius of the ECAL barrel in m (to be checked)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# trajectory of the charged particle as function of the angle phi\n",
    "# particle goes to positive x (phi = 0) axis initially \n",
    "# \n",
    "phi = Symbol('phi', real = True)\n",
    "\n",
    "# phi > 0 corresponds to bending upwards here\n",
    "# origin of the circle in this case is at (x,y) = (0, R)\n",
    "xpos = R * sin(phi)\n",
    "ypos = R - R * cos(phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1069b36d0>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAF/CAYAAACmKtU3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FNX+x/H3SUjoJARCD4QqnYAUFYEQIICKIAgiIAQV\nUK4KwhVQr1e5YkNERMTGj1As2LDQBFECSFd6lR4C0ntNO78/JokBUjbLzs7s7vf1PPMkszu7+8ls\nst/MOTPnKK01QgghBICf1QGEEELYhxQFIYQQGaQoCCGEyCBFQQghRAYpCkIIITJIURBCCJHB1KKg\nlJqqlDqmlNqSzf1BSqk5SqmNSqmtSqkYM/MIIYTImdlHCrFAhxzu/xewVWsdAUQC7yil8pmcSQgh\nRDZMLQpa6+XAmRw2SQWKpX1fDDiltU42M5MQQojsWf1f+SRgjlLqCFAU6GFxHiGE8GlWdzR3ANZr\nrcsBEcAHSqmiFmcSQgifZfWRQgzwBoDWeq9Saj9wG/BH5o2UUjJAkxBCOEFrrfKyvdVHCvFAWwCl\nVGmMgrAvqw211rZaXn75ZcszSCbvyiWZJJOrF2eYeqSglPoSaAWUVEodAl4GAgC01h8DrwLTlFKb\nAQWM0FqfNjOTqxw4cMDqCDeRTI6zYy7J5BjJZC5Ti4LW+uFc7v8baG9mBiGEEI6zuvnIY8XExFgd\n4SaSyXF2zCWZHCOZzKWcbXdyJ6WU9oScQghhJ0optId1NHusuLg4qyPcRDI5zo650jMppWSRJc+L\nq1h9SqrD1qyBfPn+WQoUgCJFjKVQIXDhPhHCcnJkLPLClUXBY5qPmjbVJCdDcjIkJcG1a3DxorFc\nvWoUhqAgKFnSWEJDja+lSkGFCsYSFmZ8LVLE6p9IiOwppaQoiDzJ7ncm7fY8VQyPKQo55UxJgcuX\n4cwZOHUKTp40lhMn4OhROHwYEhLg0CHja4ECUK2asVSv/s/XOnWgqFxPLSwmRUHklRSFW6C1USz2\n7oU9e2D3buPrrl2wc6dxZFGvnrHUrw+33w5Vq97cPBUXF0dkZKRLMrmKZHKcHXOlZ5KiIPLKlUXB\n5zqalTI++O+8Ex55BP73P/jiC/jzTzh/HhYuhJgYCAyEr76CqCgoUQKio+HFF+HHH42jECHEP6ZN\nm0aLFi2sjpGluLg4wsLCLHv9+Ph4ihYtmmuhtzpnOo/paHYHf3+oUcNYunb95/ajR2HdOmOZPNko\nJmFhkbRoAS1aQMuWRn+F1ez2ny/YMxPYM5cdM4m8Cw8PZ+rUqURFRQFQsWJFLly4YHEqx0lRcECZ\nMtCpk7GA0dm9aRMsXw6zZ8OQIUbHdnS0sbRqJZ3ZQjgrOTmZfPk876MpPbenN//5XPORK+TLBxcu\nxDF0KHz3HRw/Dp99ZhSPceOMr1FRMGEC7MtyeD9z2Pnce7uxYy47ZrrRoUOH6Nq1K6VKlaJkyZI8\n/fTT193/3HPPERISQpUqVfj5558zbo+NjaV27doUK1aMqlWr8sknn2TcFxcXR4UKFRg7dixly5bl\nscce4+rVq/Tr14+QkBBq167N2LFjr2taOXLkCN26daNUqVJUqVKF999/P+O+K1euEBMTQ0hICHXq\n1GHdunU5/kx+fn68//77VK1aldDQUEaMGJHxob53716ioqIoWbIkoaGh9OnTh3PnzmU8Njw8nLFj\nx9KgQQOKFClCr169iI+Pp1OnThQtWpRx48Zx4MAB/Pz8SE1NBeD06dP079+f8uXLExISwgMPPJBl\nrpx+RlNZPYqfgyP9abtZsmRJtvdduKD1jz9q/dhjWpcqpXXdulq/+KLW69ZpnZpqTSar2DGT1vbM\nlZ7Jjr/vWmudnJys69evr4cNG6YvX76sr169qlesWKG11jo2NlYHBAToKVOm6NTUVP3hhx/qcuXK\nZTx23rx5et++fVprrZcuXaoLFSqk169fr7U2fu58+fLpUaNG6cTERH3lyhU9cuRIHRkZqc+ePasT\nEhJ0vXr1dFhYmNZa65SUFN2oUSP96quv6qSkJL1v3z5dpUoVvXDhQq211iNHjtQtW7bUZ86c0YcO\nHdJ16tTJeGxWlFI6KipKnzlzRsfHx+saNWroKVOmaK213rNnj168eLFOTEzUJ06c0C1bttRDhw7N\neGylSpV0w4YNdUJCgr569arWWuvw8HD966+/Zmyzf/9+rZTSKSkpWmut77nnHt2zZ0999uxZnZSU\npJctW5axHypUqODQz3ij7H5n0m7P2+dtXh9gxWLXPxJHJCdrvWKF1iNGaF2tmtZVqmj9wgtab9pk\nboEQnsuuv+8rV67UoaGhGR9umcXGxupq1aplrF+6dEkrpfSxY8eyfK4uXbro9957T2ttfBgGBgbq\na9euZdxfpUoVvWjRooz1KVOmZHxgrl69WlesWPG653v99dd1//79Mx6b+cPzk08+yXhsVpRS120/\nefJk3aZNmyy3/f7773XDhg0z1sPDw3VsbOx12+RUFI4cOaL9/Pz02bNnb3ruzEUht5/xRq4sCp7X\ncOdh/P3hrruM5c03YcMGmDXL6J8oXBh69YK+faFiRauTCk/hqotX89rsfejQISpVqoSfX9atzmXK\nlMn4vlChQgBcvHiRUqVKsWDBAkaPHs3u3btJTU3l8uXL1K9fP2P70NBQAgMDM9aPHDlyXXNRhQoV\nMr4/ePAgR44coXjx4hm3paSk0LJlyywfW9GBP64btz9y5AgAx44dY8iQIfz+++9cuHCB1NRUQkJC\nsn1sbg4dOkRISAhBQUE5bpfbz2gm6VNwkjPtv0pBo0Ywdizs3w9TpsCRI9CwIbRrZ5wae/myezOZ\nzY6ZwJ65HM1kHOHf+pJXYWFhxMfHk5KSkqfHXbt2jW7dujFixAiOHz/OmTNnuOeee9JbAYCbh2ko\nW7Yshw4dyljP/H1YWBiVK1fmzJkzGcv58+eZO3duxmPj4+Mzts/8fXZu3L58+fIAvPDCC/j7+7N1\n61bOnTvHzJkzM/oGssue05ATYWFhnD59+rp+iey2y+lnNJMUBYv4+RlHD5MnG1dcP/44zJhhDMPx\nr3/B1q1WJxTies2aNaNs2bKMGjWKy5cvc/XqVVauXJnr4xITE0lMTKRkyZL4+fmxYMECFi1alONj\nevTowRtvvMHZs2c5fPgwkyZNyviwbdq0KUWLFmXs2LFcuXKFlJQUtm7dyh9//HHTYxMSEhzqoB03\nbhxnz57l0KFDTJw4kYceeggwjnQKFy5MsWLFOHz4MG+//Xauz1W6dGn27t2b5X1ly5alY8eODB48\nmLNnz5KUlMSyZctu2i63n9FMUhSc5MpzygsUgIcegp9/hs2b/zm9tVUr4wK6xET3Z3IVO2YCe+ay\nY6bM/Pz8mDNnDnv27KFixYqEhYXx9ddfA2Q5Umf6etGiRZk4cSI9evQgJCSEL7/8ks6dO2e5bbr/\n/ve/VKhQgcqVKxMdHU337t0zmpf8/f2ZO3cuGzdupEqVKoSGhjJw4EDOnz8PwMsvv0ylSpWoXLky\nHTp0oG/fvrkOGNe5c2duv/12GjZsyH333cejjz6a8Vzr168nKCiITp060a1bt1yf6/nnn2fMmDEU\nL16c8ePH3/TzzZw5k4CAAGrWrEnp0qWZOHHiTfsht5/RTD43zIWnSEoyrp6ePBl27IDBg+HJJ41B\n/oR38/Tz3M3w4Ycf8vXXX7NkyRKXP7efnx979uyhSpUqLn9ud5FhLmzA7DbpgAB48EH47TdYvBgO\nHjQG7Rs82BivyYpMzrBjJrBnLjtmssrRo0dZsWIFqamp7Nq1i/Hjx2d7Pr9wLSkKHqBOHaNTescO\nCAkx+iK6doX1661OJoQ5EhMTeeKJJyhWrBht2rShS5cuDB482JTXcuVcBN5Amo880KVLRpF4+22I\niICXXoJmzaxOJVxFmo9EXknzkY8rXNgYb2nPHrj3XujRw+iYXrXK6mRCCE8nRcFJdmj/LVDA6Hze\nvRu6d4fOnePo3Nlep7PaYT9lxY657JhJ+B4pCl4gMBAGDDAG5YuMhDZtjKuk9++3OpkQwtNIn4IX\nOn8exo+HSZOgf3/4z3+M+auFZ5A+BZFX0qcgclSsGLzyitGMdOYM3HYbfPKJMZe1EELkRIqCk+zY\n/ntjpjJljLOUFiwwxlVq1AjcHduO+wnsmcuOmXzFhx9+SOnSpSlWrBhnzpyxOo6lpCj4gIYNYckS\n+O9/oV8/YzrRY8esTiU8UXh4OIUKFaJo0aIZyzPPPJNx/99//81jjz1GuXLlKFasGLVq1eKVV17h\ncqaRHrXWVKlShTp16tz0/JGRkfzf//2fW36WdElJSQwfPpxff/2V8+fPXzcyqS8ytSgopaYqpY4p\npbbksE2kUmqDUmqrUirOzDyuZMdxanLKpBR06wbbtkHZslCvHnz4oflNSnbcT2DPXHbMdCOlFHPn\nzuXChQsZS/rYPadPn+bOO+/k2rVrrF69mvPnz/PLL79w9uzZ6waIW7ZsGSdOnGD//v03DfCW1RhK\nZkpJSeHo0aNcvXqVWrVque117czsI4VYoEN2dyqlgoEPgE5a67rAgybn8XlFihhDd//2m9GkdNdd\nRqEQ4laNHz+eoKAgPvvss4w5DCpUqMCECROoV69exnbTp0+nc+fOdOzYkenTpzv1WidPnuS+++6j\nePHilChR4rp5Bvz8/NiXaR7cmJgYXnrpJeDmqT8feeSRjGIQHBxM27ZtARgyZAgVK1YkKCiIxo0b\n8/vvv2c8X2pqKq+//jrVqlWjWLFiNG7cmISEBAB27txJu3btKFGiBDVr1uSbb75x6uezkqlFQWu9\nHMipga4X8J3WOiFt+5Nm5nElO7b/5iVT3bqwdCk89phxGuuYMcYgfFZmcic75rJjpqxkd2bU4sWL\n6dq1a46PvXz5Mt999x19+vShd+/ezJo1iyQnfvHeeecdwsLCOHnyJMePH+eNN97Idtsbjz6OHTvG\nmTNniI+PJzY2lm1p/xWdO3eOxYsXA8bQ1Zs2beLMmTP06tWL7t27k5g2XPE777zDrFmzWLBgAefP\nnyc2NpZChQpx6dIl2rVrR58+fThx4gSzZs1i8ODB7NixI88/n5Ws7lOoDoQopZYopf5QSj1icR6f\n4ucHAwcaYyitWAFNmxozwwmRHa01Xbp0oXjx4hlLeh/A6dOnKVu2bI6Pnz17NgUKFCA6Opp7772X\npKQk5s2bl+ccgYGB/P333xw4cAB/f3+aN2+ea+50fn5+jB49moCAAPLnz59lkevduzfFixfHz8+P\nYcOGce3aNXbt2gXAlClTeO2116hevToA9erVIyQkhLlz51K5cmX69euHn58fERERdO3a1eOOFqye\njjMAaAS0AQoBq5RSq7XW2YwDah92bP91NlNYGMyfb0zy0749PP00PP885HPBb4cd9xPYM5ejmdRo\n17S565fzfi2EUooff/yRqKiom+4rUaJExjSW2Zk+fTrdu3fHz8+P/Pnz07VrV6ZPn06XLl3ylOO5\n557jlVdeITo6GoCBAwcycuRIhx5749SfWRk3bhxTp07lyJEjKKU4f/48J08aDRkJCQlUrVr1pscc\nPHiQNWvWXNdRnZycTN++fR39sWzB6qJwCDiptb4CXFFKLQMaADcVhZiYGMLDwwGj7S8iIiLjjyj9\nsFvWb229X79I2rQxhsv46iv46adIqlSxTz5fWc+NMx/m7tC2bVu+//57Xn755Sw7ixMSEvjtt99Y\nt24d3333HUDGDG6nT5++ae7jnBQpUoRx48Yxbtw4tm3bRlRUFE2bNqV169YUKlTourOd/v777+vm\nUc6tI3v58uW8/fbb/PbbbxlnSIWEhGQcUYSFhbFnzx5q16593eMqVqxIq1atcp1VzkxxcXFMmzYN\nIOPzMs+01qYuQDiwJZv7agKLAX+MI4UtQO0sttN2s2TJEqsj3MRVmVJStB4/XuuSJbWeOlXr1FTr\nM7maHXOlZ7Lj73u68PBwvXjx4izvO336tA4PD9ePPPKIPnjwoNZa64SEBD1s2DC9efNm/frrr+va\ntWvrY8eOZSxHjx7VVapU0e+//77WWuvIyEj90Ucf6StXrmQsiYmJN73W3Llz9e7du3VqaqqOj4/X\nZcuW1XFxcVprrZs3b65HjRqlk5OT9YIFC3TBggX1Sy+9pLU29nGFChWue679+/drpZROSUnRWms9\nb948Xa5cOX306FF97do1PXr0aO3v769//fVXrbXWb7/9tq5fv37G62/atEmfOnVKX7hwQVeqVEnP\nnDlTJyYm6sTERL127Vq9Y8cOF+z5nGX3O5N2e54+s80+JfVLYCVwm1LqkFLqUaXUIKXUoLRP+p3A\nz8BmYA3wqdZ6u5mZRO78/ODZZ40zlN59Fx5+2Bg6QwiATp06XXedQrdu3QAoXrw4K1euJCAggGbN\nmlGsWDHatm1LcHAw1apVY8aMGQwePJhSpUplLKVLl+aJJ55gxowZGc//5JNPUqhQoYzlscceuynD\n7t27adeuHUWLFuWuu+7iX//6F61atQLgvffeY86cORQvXpwvvvjipsl5sjpSyHxbhw4d6NChAzVq\n1CA8PJyCBQtmnE0FMGzYMHr06EF0dDRBQUEMGDCAq1evUqRIERYtWsSsWbMoX748ZcuW5fnnn8/o\noPYUMvaRyNGVK0aB+PVX+OYbY/4GYS4Z+0jklYx9JNymYEH46CP43/+gXTv4+GOQzyshvJcUBSfZ\n8ZxyMzM9/LBx2urkycaw3FeuWJ/pVtgxlx0zCd8jRUE4rEYNY3a31FS4+26Ij7c6kRDC1aRPQeSZ\n1sZ8DePGwZdfGldEC9eRPgWRV9KnICylFAwfDjNnQs+e8MEHVicSQriKFAUn2bH9192Z2raFlSuN\nGd6GDs16xFU77iewZy47ZhK+R4qCuCVVqhiFYcsW6NIFLl60OpEQ4lZIn4JwiaQkePJJ+PNPmDsX\nype3OpHncud8AsJ7uKpPQYqCcBmt4c03jWsZFi405oYWQlhHOprdyI7tv1ZnUsoYXfXll40zktat\nsz5TduyYSzI5RjKZS4qCcLn+/eGTT+Dee43CIITwHNJ8JEyzYgV07WpcBZ02ZpoQwo2caT6yej4F\n4cWaNzf6Fjp2hGvXoFcvqxMJIXIjzUdOsmMboh0znT0bxy+/wHPPQWys1Wn+Ycd9JZkcI5nMJUcK\nwnR16xpzM7RtC4mJMGiQ1YmEENmRPgXhNnv3QuvW8Mor8OijVqcRwvtJn4KwtapVYfFiozDkzw+9\ne1udSAhxI+lTcJId2xA9IVONGvDLL/DvfxszuVnFE/aVHUgmx9gxk7PkSEG4Xe3a8PPP0L69MbPb\nffdZnUgIkU76FIRl1q41CsIPP8Bdd1mdRgjvI8NcCI/StCl89hk88ABs3Wp1GiEESFFwmh3bED0x\nU3Q0vPuucYHbwYPuyQSeua+sIJkcY8dMzpI+BWG5Xr3gxAno0MGYm6F4casTCeG7pE9B2Mazz8Km\nTUYndGCg1WmE8Hwyn4LwaCkpxgB6JUvClCnGUNxCCOdJR7Mb2bEN0dMz+fvD55/Dhg3w1lvmZQLP\n31fuIpkcY8dMzpI+BWErRYrAnDlwxx1Qs6Yx77MQwn2k+UjY0tq1xiQ9y5ZBrVpWpxHCM0nzkfAa\nTZvC2LHGkcK5c1anEcJ3mFoUlFJTlVLHlFJbctmuiVIqWSnV1cw8rmTHNkRvy9S/P7RrB336QGqq\n6zKB9+0rs0gmx9gxk7PMPlKIBTrktIFSyh94C/gZkPNNxHXefdc4UhgzxuokQvgG0/sUlFLhwByt\ndb1s7h8KJAJNgLla6++y2Eb6FHzYkSNw++3GmUlRUVanEcJzeFyfglKqPNAZ+DDtJvnkFzcpVw6m\nT4dHHoGjR61OI4R3s/qU1AnAKK21Vkopcmg+iomJITw8HIDg4GAiIiKIjIwE/mnPc+f6xo0bGTp0\nqGWvn9V6+m12yZM5y60+X2AgPPpoJH36wPPPx+HvL++fJ71/rlyfMGGC5X//N67b5fcpLi6OadOm\nAWR8XuaZ1trUBQgHtmRz3z5gf9pyATgG3J/FdtpulixZYnWEm3h7pqQkrVu10vrVV2/9ubx9X7mK\nZHKMHTNprXXaZ2eePrMt71PItF1s2nazs7hPm51TeIaEBGjUCObNgyZNrE4jhL3Zrk9BKfUlsBK4\nTSl1SCn1qFJqkFJqkJmvK7xXhQrw/vvGaaqXLlmdRgjvY2pR0Fo/rLUup7UO1FqHaa2naq0/1lp/\nnMW2/bM6SrCrzG2tduErmR56yLi47d//dv45fGVf3SrJ5Bg7ZnKWXNEsPNKkSbBgAcyfb3USIbyL\njH0kPNaSJdC3rzGVZ1CQ1WmEsB+ZT0H4nCeeMOZh+PRTq5MIYT+262j2ZnZsQ/TFTGPHwqJF8Msv\neXucL+4rZ0gmx9gxk7OkKAiPVqwYfPwxDBwIFy9anUYIzyfNR8Ir9O0LpUvD229bnUQI+5A+BeGz\njh+HunXh11+hXo6XSQrhO6RPwY3s2Iboy5lKlYL//Q+efNKxuRd8eV/lhWRyjB0zOUuKgvAaAwdC\nUhKkjQcmhHCCNB8Jr7JhA3ToALt2QXCw1WmEsJb0KQiBccRQtCi8847VSYSwlvQpuJEd2xAlk+HV\nV41Jef76K/ttZF85RjI5xo6ZnCVFQXid0qVhxAh47jmrkwjheTym+ej1Za+Tzy/fTUuAfwCFAwpT\nJLBIlkvhwML4Kal9vubaNahdGz75BNq0sTqNENZwpvnI6uk4HXYh8QLJqck3LYkpiVxKusTFxIvX\nLReuXeBC4gUSUxIJKRhCyUIlM5bQQqEZ35cvWp4KxSoQFhRGmSJlyOfnMbtE5CB/fnjjDRg5Etat\nA5WnPwshfJfHHCk4mzMxJZHTV05z8vJJTl4+yYlLJzK+P37pOEcuHiHhfAKHzh3i5OWTlCpcKqNI\nVAqqRLWQalQPqU61kGqEBYVlHHXExcVlzJFqF5Lpeqmp0LgxvPACPPigfXJlRzI5RjI5zquPFJwV\n6B9ImSJlKFOkTK7bJqUk8ffFvzOKxIGzB/jzyJ/M2jqLPaf3cOrKKSoHV6ZaSDUKHCrAvqB91C1V\nlzqhdSgcWNgNP43ICz8/42hhyBDo0gXyef1vuxC3zuuPFFzpUuIl9p3Zx+7Tu/nr1F9sP7GdLce3\nsOvkLsoVLUfdUnWpV6oe9UrXo0HpBlQvUV36MyymNURFGdN3PvaY1WmEcC+5TsEiyanJ7Dm9hy3H\ntrD1+Fa2HN/ChqMbOHPlDLeXu50m5ZoYS/kmhBULQ0kDt1utWgU9e8Lu3RAYaHUaIdxHioIbOdKG\neOLSCf448gfrjqwzlsPr0GialGvC3RXvpkXFFjQu15j8+fK7LZO72SVT+/bQvTs8/rixbpdcmUkm\nx0gmx0mfgs2EFg6lY/WOdKzeEQCtNQnnE1h7eC2/x//OMz8/w66Tu2hcrjEtK7WkRcUW3Bl2J0UC\ni1ic3Pu89JIxvHa/fhAQYHUaIexLjhQsdv7aeVYeWsmyg8tYHr+cDX9voH7p+kRXjSa6ajRNyzeV\n02RdJCrKKAr9+lmdRAj3kOYjL3Al6QorD61k0d5FLNq3iANnD9A6vHVGkahSvIrVET3WkiUwaBDs\n2AH+/lanEcJ8MvaRG5k11knBgIK0qdKGt9q9xYZBG9j5r510q9WNVQmraD61OdXfr87whcNZdnAZ\nyanJbsl0K+yUKTISSpaEH36wV650kskxkslcUhRsrnSR0vSu35vpXaZzZNgRvn7wa4rmL8rQn4dS\nZlwZ+v3Qj9k7ZnMxUSYozo1SMHy4jJ4qRE6k+ciDxZ+L56ddP/HTrp9YnbCalpVa0qNODzrf1pmg\nAkFWx7OllBSoUQNmzoS77rI6jRDmkj4FH3bu6jnm7Z7H19u+ZsmBJURVjuKhOg/RqUYnudr6BpMm\nGf0L331ndRIhzCV9Cm5ktzbEoAJBlDtVjh96/sDBoQfpfFtnpm+aTrnx5ej5bU++3/E915KvuT2X\n3fYTQP/+sHhxHHv3Wp3kenbcV5LJMXbM5CxTi4JSaqpS6phSaks29/dWSm1SSm1WSq1QStU3M4+v\nCC4QTExEDAt6L2DvM3tpHd6a99a8R/nx5Xl6/tOs/3s9vnzkVbgwdOwIH31kdRIh7MfU5iOlVAvg\nIjBDa10vi/vvBLZrrc8ppToAr2it78hiO2k+coH9Z/YzY9MMpm2aRtHAovSP6E+f+n0ILRxqdTS3\n27sX7rgD4uOhYEGr0whhDlv2KSilwoE5WRWFG7YrDmzRWlfI4j4pCi6UqlNZemApsRtj+WnXT7Su\n3JoBjQbQvmp7/P185wT+e+6Bhx6Si9mE9/L0PoXHgPlWh3CUHdsQHc3kp/xoXbk1Mx6YQfyz8dxb\n/V7+u+S/VH+/OmNXjOXk5ZNuz+RucXFxDB4MkydbneQfdtxXkskxdszkLFsUBaVUa+BRYKTVWXxN\nsfzFeLzR46wbsI5ZD85i+4ntVH+/Ov1+6MeahDVe3ffQsSMcPQobNlidRAj7sHxQnbTO5U+BDlrr\nM9ltFxMTQ3h4OADBwcFERERkjEqYXqXdvZ7Oqtc3Y71p+ab8+POP/LznZ3rF9yK4QDDt/NrRpnIb\n2rVpl+fni4yMtNXPl3nd3x9iYuDVV+N45hnr89hx3Y7vX/ptdsljp8+DuLg4pk2bBpDxeZlXlvYp\nKKUqAr8BfbTWq3N4DulTsECqTmXhnoW8u/pdtp3YxtNNn2bQ7YMoXrC41dFcZv9+aNoUEhKMeZ2F\n8Ca261NQSn0JrARuU0odUko9qpQapJQalLbJf4HiwIdKqQ1KqbVm5nGlG/87sANXZ/JTfnSs3pFF\njyxifq/57Di5g6oTqzJkwRD2n9lvSSZXSc9VuTLUrQtz51qbB+y5rySTY+yYyVmmFgWt9cNa63Ja\n60CtdZjWeqrW+mOt9cdp9z+utS6htW6YtjQ1M49wXoMyDZjeZTpbntxCwYCCNPm0Cd2/6c6fR/60\nOtot698fYmOtTiGEPcgwF8IpF65dYMr6KYxbNY6GZRryUsuXaFahmdWxnHLpEpQvb0zXGep7l2wI\nL2a75iPGb2bgAAAgAElEQVThvYrmL8qzdz7L3mf2ck/1e+j+TXfaf9ae3+N/tzpanqVf4SxjIQkh\nRcFpdmxDtCJTgXwFGNxkMHue2cODtR6k7/d9iZoexZL9SyzL5Igbc/XsCbNmWZMlnR33lWRyjB0z\nOUuKgnCJQP9ABtw+gF1P7aJvg74MnDuQtjPasvPETqujOaRDB9i8GQ4ftjqJENaSPgVhiqSUJGI3\nxvK/pf+jWYVmjGk9hlqhtayOlaOYGIiIgKFDrU4ihGtIn4KwjQD/AAbePpDdT+/mjvJ30HJaSx79\n8VHiz8VbHS1b3bvD7NlWpxDCWlIUnGTHNkQ7ZlqzYg3PNX+O3U/vplzRcjT8uCHDFg7jzJVsL153\ni6z2VZs2sGkTnHTd0E95Ysf3TzI5xo6ZnCVFQbhFcIFgxkSNYdvgbVxOukzND2rywdoPSE5Ntjpa\nhgIFoG1be1zIJoRVpE9BWGLLsS08u/BZ/r74N+Ojx9O+WnurIwEwYwb88IM0IwnvYMv5FFxBioJ3\n0loz96+5DF80nOolqvNO9DvULFnT0kwnT0LVqsboqTL5jvB00tHsRnZsQ/S0TEopOt3Wia2Dt9K2\ncltaxLZg2MJhXLh2wbJcJUtCvXqwfLnpEW7iae+fVSSTubItCkqpEAeWYHeGFd4p0D+QZ+98lu2D\nt3P26llqT67Nt9u/tWwuh+hoWLTIkpcWwnLZNh8ppa4BR3J5fD6tdZjLU92cRZqPfMjyg8t5ct6T\nhAWFManjJKqGVHXr669eDQMHGhezCeHJXNqnoJTaqLWOyOUFc93GFaQo+J6klCQmrJ7AWyve4plm\nzzCy+Ujy53PPhAfJyVCqFGzbBmXLuuUlhTCFq/sU7nDg8Y5s45Xs2IboTZkC/AN4rvlzrB+0ng1H\nN1D/o/osP+i6hv6ccuXLB1FR8MsvLns5h3jT+2cmyWSubIuC1vpq+vdKqeJKqQZKqUZKqduVUo1u\n3EYIM1QMqsj3D33PW23foud3PXlmwTNcTLxo+utGRYEX/Z0L4bBcT0lVSr0KxAD7gNT027XWrU1N\ndn0GaT4SnL5ymmELh7H04FKmdJpCmyptTHutLVuga1djjgUhPJUp1ykopf4C6mqtE28l3K2QoiAy\nW7B7AYPmDqJjtY6MbTeWoAJBLn+N1FTj9FTpVxCezKzrFLZizKMsMrFjG6KvZOpYvSNbntwCQN0P\n6/LL3rw3/ueWy88PWrRw7/UKvvL+3SrJZC5HisLrwHql1CKl1Jy05SezgwmRk6ACQXzc6WNiO8fy\n6E+PMvTnoVxJuuLS12jRApYtc+lTCmF7jjQfbQc+wjhiSO9T0FrrpSZny5xBmo9Etk5fOc2guYPY\ncWIHn3f9nAZlGrjkeVesgCFD4I8/XPJ0QridWX0K67TWTW4p2S2SoiByo7Vm5uaZDF80nJHNRzLs\nzmH4qVsbxeXyZaNf4cwZyO+eSySEcCmz+hSWK6XeUErdmXZKaqP0U1J9mR3bEH05k1KKvg36svbx\ntfyw8wfazGjDoXOHbilXoUJQvboxx4I7+PL7lxeSyVyOFIVGGBepvQ68k2kRwnYqF6/M0piltK3c\nlsafNmbeX/Nu6fmaNIF161wUTggPIENnC6+1/OByes3uRa+6vRgTNYYA/4A8P8fHH8OqVTBtmuvz\nCWE2lzYfKaXuc+AFc91GCKu0qNSC9QPXs/n4ZiKnR+bYnJSdhg3d13wkhB3k1Hw0LvOwFlkstwNv\nuCuo3dixDVEy3Sy0cCjzes2jU41ONPm0CfN3z89Trjp1YNcuY5A8s1m9r7IimRxjx0zOypfDfUfJ\nve/gLxdmEcIUfsqPUXePonlYc3rN7kXf+n2JUlEOPbZwYShXDvbsgZrWTgonhFuY2qeglJoK3Asc\n11rXy2abiUBH4DIQo7XekMU20qcgXOL4peM89O1DFMxXkM+7fk7xgrlfrP/AA9CrF3Tv7oaAQriQ\nHafjjAU6ZHenUuoeoJrWujowEPjQ5DzCx5UqXIpFfRZRo0QNmk5pyrbj23J9TL16xgB5QvgCU4uC\n1no5cCaHTe4HpqdtuwYIVkqVNjOTq9ixDVEyOSbAP4AuBbrwUsuXiJweyewds3PcvlYto1/BbHbc\nV5LJMXbM5Kyc+hTcoTyQ+ZSQBKACcMyaOMKX9G3Qlzqhdej6dVfW/72e0ZGj8ffzv2m7atVkCG3h\nOxwZ5qIwMAyoqLUeoJSqDtymtZ7r0AsoFQ7MyapPQSk1B3hTa70ibX0xMEJrvf6G7aRPQZjm+KXj\n9PimB8XyF+OLbl9QJLDIdfefOQOVKsG5c6Dy1DorhLWc6VNw5EghFvgTuCtt/QjwLeBQUcjFYSAs\n03qFtNtuEhMTQ3h4OADBwcFEREQQGRkJ/HPoJuuy7sz69nXbeTHsRWZdnEXL2JY8X+F5QguHZty/\naVMcSsHx45GULm19XlmX9ezW4+LimJZ2pWX652Weaa1zXIA/075uyHTbptwel2nbcGBLNvfdA8xP\n+/4OYHU222m7WbJkidURbiKZHJdVrtTUVP3G8jd0hfEV9Poj66+7r1kzrX//3f2ZrCaZHGPHTFpr\nnfbZ6dBndfriSEfzNaVUwfQVpVRV4JojBUcp9SWwErhNKXVIKfWoUmqQUmpQ2if9fGCfUmoP8DEw\n2JHnFcIMSilG3T2K8dHjif4smjm75mTcV6UK7N9vYTgh3MSRPoVo4EWgNvAL0BzjeoIl5sfLyKBz\nyymEK61JWMMDXz3AyOYjeabZM4wcqQgJgVGjrE4mhONMmU8h7YlLYjTvAKzRWp9wIp/TpCgIKxw4\ne4D7vriP6KrRVNo1jt1/+TFpktWphHCcKRevKaV+1Vqf1FrPTVtOKKV+dT6md0jv3LETyeQ4R3KF\nB4ezvP9y1h5ey2zdl/jDiZZncjfJ5Bg7ZnJWTqOkFlRKlQBClVIhmZZwjOsLhPB6xQsWZ9Eji9CB\n51la7n4uJV6yOpIQpsq2+UgpNRQYApTDOA013QXgE6212w6kpflIWO3goWRqjRhI/TbbmddrHiUK\nlbA6khC5MmuO5me01hNvKdktkqIgrHb1KhQL0gyb9wI/7vqBhX0WUjGootWxhMiRKX0KWuuJSqm6\nSqkeSqm+6YvzMb2DHdsQJZPj8pqrQAHIH6h4oekbDLp9EC1iW7D7lGvHvrDjvpJMjrFjJmflekWz\nUuoVoBVQB5iHMcz178AMU5MJYTMlS8LJkzD0jqEUDSxK6+mtWfTIImqH1rY6mhAu40jz0VagAbBe\na90gbRTTz7XWbd0RMC2DNB8JyzVpApMmQbNmxvpnmz9jxC8j+LnPz9QvXd/acEJkwayxj65orVOU\nUslKqSDgONePVySETwgJgdOn/1nvU78P+f3zEz0zmnm95nF7udutCyeEizgyzMU6pVRx4FPgD2AD\nxtAVPs2ObYiSyXHO5CpaFC5evP627nW68/F9H9Px846sOrTK7ZnMJpkcY8dMzsr1SEFrnT4e0UdK\nqYVAUa31ZnNjCWE/RYrcXBQAOtfsTKB/IJ1ndWb2Q7O5u+Ld7g8nhIs4OsxFeaASRhFRGCPvLTM5\nW+bXlz4FYbmnnoLbboOnn876/l/2/kLv2b2Z22suTcs3dW84IbJgSp+CUuot4CFgO5CS6S63FQUh\n7CC7I4V07aq2Y2rnqXT6shMLei+gUdlG7gsnhIs40qfwAMZMa/dorTulL2YHszs7tiFKJsc5k6tw\nYbiUyygX99W4j4/u/Yh7Pr+HLce2mJ7JbJLJMXbM5CxHzj7aCwTi4BwKQnirgIDciwLAA7UeIDEl\nkfaftee3fr9Rs2RN88MJ4SKOXKcwG+M6hV/5pzBorfUzJmfLnEH6FITlxo2Do0eNr46YsWkGL/72\nIkv6LaFaSDVzwwmRBbOuU/gpbclMPqGFz8mXD5KTHd++b4O+XE2+SvTMaH5/9HfKFS1nXjghXMSR\nsY+mZbFMd0c4O7NjG6JkcpwzufJaFAAG3j6QAY0G0OGzDpy9etblmcwmmRxjx0zOymk+hW/Svm7J\nYpHrFITP8ffPe1EAGHX3KFqHt+b+L+/nStIV1wcTwoVymk+hnNb6SNqkOjfRWh8wL9ZNWaRPQVhu\n0iTYsQM++CDvj03VqfSZ3YfLSZf5tse35PNzpOVWiFvj0qGztdZH0r4eyGq5xaxCeJzkZOMMJGf4\nKT+mdZnGleQrPDH3CeSfHGFXOTUfXVRKXchmOe/OkHZkxzZEyeQ4Z3IlJxv9Cs4K9A/kux7fsfnY\nZl5a8pJLMplNMjnGjpmcle2vuNa6CIBSagzGdJyfpd3VG2OKTiF8yq0WBYAigUWY12sed/7fnVQt\nXpX+Dfu7JpwQLuLIdQqbtdb1c7vNTNKnIOzg1VeNaTlfe+3Wn2vnyZ20mtaKWd1m0bpy61t/QiGy\nYMp0nMAlpVQfpZR/2tIbyGEEGCG80+XLxlAXrlCzZE2+7PYlPb/ryc6TO13zpEK4gCNFoRfQAziW\ntvRIu82n2bENUTI5zplcFy8acyq4SlTlKN5s8yb3fXEfJy+ftOW+kkyOsWMmZ+XYQqqU8gf+pbW+\n3015hLCtixeNkVJdqX/D/uw+vZsus7rw30r/de2TC+EER/oUVgN3WtmoL30Kwg66d4cePYyvrpSq\nU+n5bU/y58vPjC4zUCpPTcBCZMusPoWNwI9KqUeUUt3Slq4OBuqglNqplNqtlBqZxf0VlVJLlFLr\nlVKblFId8xJeCHe6cMF1fQqZpV/DsPX4ViasnuD6FxAiDxwpCgWA00AUcF/akut8CmlNT5OADkBt\n4GGlVK0bNvsPMEtr3QjoCUx2PLq17NiGKJkc50yu06chJMT1WQAKBRRiZPmRvLXiLX7b/5s5L+IE\nO75/kslcjszRHOPkczcF9qRf/ayUmgV0BnZk2iYVCEr7Phg47ORrCWG6kychNNS85y9TpAxfdPuC\nXt/1YvXjqwkPDjfvxYTIhiN9CmHARCB9NvJlwBCtdUIuj3sQaK+1HpC23gdoprV+OtM2ZYBFQHGg\nMNBGa70hi+eSPgVhuWLF4NAhCArKfdtb8e6qd5mxeQYrHl1BoYBC5r6Y8Gpm9SnEYsynUC5tmZN2\nW24c+RTvBcRqrcOAe/jnqmkhbOXaNePCtWLFzH+toXcMpU5oHQbMGSBjJAm3c+Si/VCtdeYiME0p\n9awDjzsMhGVaDwNuPLp4FGgPoLVerZQqoJQqqbU+eeOTxcTEEB4eDkBwcDARERFERkYC/7TnuXN9\n48aNDB061LLXz2o9/Ta75MmcxS550tfz+v6dPAklSkSilPnv39KlS+lTrA8v7H2BiWsm0uBqA7fv\nn/R1O75/EyZMsPzv/8Z1u3wexMXFMW3aNICMz8s801rnuAC/AY8A/hhFpA/wqwOPy4cxv3M4xhzP\nG4FaN2wzH+iX9n0t4HA2z6XtZsmSJVZHuIlkclxec/3xh9YNG5qTJd2Nmfae3qtDx4bqNQlrzH3h\nHNjx/ZNMjkv77Mz1cz7z4kifQjjwPnBH2k0rgae11vG5FZy0U0wnpBWU/9Nav6GUGg38obWek3Y2\n0qdAEYzmpue01ouzeB6dW04hzPTjj/B//wc/3Tgxrclm75jN8EXDWT9wPcULFnfviwuP50yfQq5F\nwQ6kKAirTZoE27fDZAtOmh6yYAgHzx3k+4e+lwvbRJ6Y0tGslApTSn2vlDqRtnynlKrgfEzvkLmt\n1S4kk+PymishAcLCct/uVmSXaWy7sRy+cJj31rxnboAs2PH9k0zmMvPsIyG8xqFDUMGif4Xy58vP\nVw9+xevLX2dNwhprQgif4UifwiatdYPcbjOTNB8JqzVvDm+8AS1bWpfh+x3f8+zCZ9n0xCaCCph8\nsYTwCmZdp3Aqbdwjf6VUvrSL0G46ZVQIb7Z7N1Svbm2GB2o9QMdqHXlqwVPWBhFezZGi8CjGHApH\ngb+B7oDPzyFoxzZEyeS4vOQ6d86YYKdMGfPygGOZxkWPY+3htczaOsvcMGns+P5JJnM5MvbRARwY\nAE8Ib7V3L1SrBnY48adwYGE+7/o593x+D83DmhMWZHLvt/A5jvQpzACe0VqfTVsvDryjtX7UDfnS\nM0ifgrDMV1/Bt9/CN99YneQfry17jV/3/8rivovxU44c8AtfZFafQv30ggCgtT4DNMprOCE81Y4d\nUKOG1SmuN+ruUSSmJPLuqnetjiK8jCNFQSmlQjKthGBcoezT7NiGKJkcl5dcW7ZAvXrmZUmXl0z+\nfv7MfGAmb654k83HNtsik7tIJnM5UhTeAVYppV5VSo0BVgFvmxtLCPvYutU9RSGvKhevzJtt3qT/\nj/1JTk22Oo7wEg4Nc6GUqoMx85oGftNabzc72A2vL30KwhJXrhizrZ0/DwEBVqe5mdaa9p+1J6py\nFKPuHmV1HGEzMvaREC7255/Qvz9sNq+F5pYdPHuQxp82ZlnMMmqF3jjjrfBlZnU0iyzYsQ1RMjnO\n0VwbNkADN1277+y+qhRcidGRo3nsp8dISU2xRSYzSSZzSVEQIgfr1kGTJlanyN0TjZ8gwD+AiWsm\nWh1FeDhHrlN4BpiZdiqqJaT5SFilUSP44AO4806rk+Ruz+k93DHlDlY/vppqIdWsjiNswKzmo9LA\nOqXU10qpDkoGdBc+4upV2LkTIiKsTuKYaiHVeP7u53ly3pMyt7NwWq5FQWv9IlADmArEALuVUq8r\npaqanM3W7NiGKJkc50iujRvhttugYEHz84Br9tWQO4Zw/NJxvtr21a0Hwp7vn2Qyl0N9ClrrVIwB\n8Y4BKUBx4FullFyvILzW6tXQrJnVKfImn18+Prz3Q4YvGs65q+esjiM8kCN9CkOAvsApYArwvdY6\nSSnlB+zWWpt+xCB9CsIKXbvCgw9Cr15WJ8m7AT8NoGBAQSZ2lI5nX2bKdQpKqdHAVK31wSzuq+2O\nC9mkKAh30xpKlYL1682fhtMMpy6fos7kOszvPZ9GZWWoMl9lSkez1vrlrApC2n1uvbLZTuzYhiiZ\nHJdbrp07oUgR9xYEV+6rEoVK8EabN3hy3pO3dO2CHd8/yWQuuU5BiCwsX27t1Juu0C+iH4H+gUxZ\nP8XqKMKDyDAXQmShZ0+IjoZH3TZriDk2Ht1Ih886sOupXTKvsw+SsY+EcIGUFChd2hjiwhP7E240\n4KcBBBcI5u1oOVnQ18jYR25kxzZEyeS4nHJt2GB0Mru7IJi1r8ZEjSF2Yyx7Tu/J82Pt+P5JJnNJ\nURDiBosWGU1H3qJ0kdI8d9dzPPfLc1ZHER5Amo+EuEFkJIwYAffcY3US17mafJXaH9Rmyv1TiKoc\nZXUc4SbSfCTELTp71rg2oVUrq5O4VoF8BXi73ds8u/BZlw+vLbyLqUUhbQC9nUqp3Uqpkdls00Mp\ntU0ptVUp9bmZeVzJjm2Ikslx2eWaP98oCIULuzcPmL+vutbqSlD+IKZvmu7wY+z4/kkmc5lWFJRS\n/sAkoANQG3hYKVXrhm2qA6OAu7TWdYEhZuURwhE//QSdO1udwhxKKd5s+yavxL3C1eSrVscRNmVa\nn4JS6k7gZa11h7T1UQBa6zczbTMW2Km1nprLc0mfgjBdYqJxKuqOHVCmjNVpzHP/l/cTVTmKoXcM\ntTqKMJnd+hTKA4cyrSek3ZZZdeA2pdTvSqlVSqn2JuYRIkdxcVCzpncXBIDXol7jjd/f4MK1C1ZH\nETZkZlFw5F/7fEA1oBXwMPCpUsojLru0YxuiZHJcVrm++cYYGdUq7tpX9UrXI7pqNONXjc91Wzu+\nf5LJXPlMfO7DQObLf8IwjhYySwDWaK1TgANKqb8wisSfNz5ZTEwM4eHhAAQHBxMREUFkZCTwzxvi\nzvWNGzda+vpZraezSx47r9/4/iUlwfffR/Lnn77x/nXM15Gn1z7N4CaD2bZumyU/r7PrGzdutFUe\nO30exMXFMW3aNICMz8u8MrNPIR+wC2gDHAHWAg9rrXdk2qZ92m0xSqmSwHqgwY3zQUufgjDb/Pnw\n2muwYoXVSdznqflPEeAXwLsd3rU6ijCJrfoUtNbJwFPAQmA78JXWeodSarRSqlPaNguBU0qpbcBv\nwL9vLAhCuMOsWfDQQ1ancK//tPwP0zdN5+8Lf1sdRdiIqdcpaK0XaK1v01pX01q/kXbby1rrOZm2\nGa61rqO1rq+1/trMPK504yG/HUgmx2XOdeUKzJkD3btblwfcv6/KFCnDI/UfYdzKcdluY8f3TzKZ\nS65oFj7vhx+gSRMoW9bqJO43ovkIYjfGcvzScaujCJuQsY+Ez4uOhv794eGHrU5ijcHzBlMsfzHe\nbPtm7hsLjyLzKQiRR4cOQUQEJCRAwYJWp7HGwbMHafRJI/566i9KFCphdRzhQrbqaPZ2dmxDlEyO\nS881Ywb06GGPgmDVvqoUXImuNbsyYfWEm+6z4/snmcwlRUH4rNRUiI2FmBirk1jv+RbP8+EfH3L2\n6lmrowiLSfOR8FmLFsHIkcZQ2SpPB9jeqffs3kSUjuC55jIZj7eQ5iMh8mDyZBg8WApCuuF3Dmfi\n2okkpiRaHUVYSIqCk+zYhiiZHPfVV3EsXw69elmd5B9W76tGZRtRPaQ6X2/753IhqzNlRTKZS4qC\n8Elz50KfPtZMpmNn/77r37yz6h2kudZ3SZ+C8DlXr0J4+D9DZYt/pOpU6k6uy/sd36dNlTZWxxG3\nSPoUhHDAzJnQuLEUhKz4KT+G3zmcd1a9Y3UUYREpCk6yYxuiZMpdaiq88w60aRNndZSb2GVf9a7f\nmw1HN7D9xHbbZMpMMplLioLwKfPnG/0IERFWJ7GvAvkKMLDRQCavm2x1FGEB6VMQPqV1axg40HfH\nOXJUwvkE6n9Yn4NDD1I0f1Gr4wgnSZ+CEDlYuRIOHIAHH7Q6if1VKFaB1pVb8/mWz62OItxMioKT\n7NiGKJly9uqr8PzzEBBgr1zp7JZpcOPBjP18rO1OT7XbfgJ7ZnKWFAXhE9atg23boF8/q5N4jqjK\nUSSnJvN7/O9WRxFuJH0Kwifcf78xb8JTT1mdxLO8t/o9Vh9ezZfdvrQ6inCCzKcgRBbWr4f77oO9\ne+0xRLYnOXv1LJXfq8zOf+2kdJHSVscReSQdzW5kxzZEyZS1F180lswFwQ65bmTHTBtXb6TzbZ35\nbPNnVkfJYMf9ZMdMzpKiILxaXBzs2gUDBlidxHP1j+hP7MZY23U4C3NI85HwWlrDXXfBv/5lDH4n\nnJOqU6n+fnW+evArGpdrbHUckQfSfCREJnPmwKVLcqHarfJTfvRr0I/YDbFWRxFuIEXBSXZsQ5RM\n/0hKMmZVe/118Pe/+X7ZV45Jz9SvQT9mbZvF1eSr1gbC3vvJG0hREF7po48gLAzuvdfqJN6hUnAl\nGpZpyI87f7Q6ijCZ9CkIr3P6tDEs9pIlUKeO1Wm8x+ebP2fm5pn83Odnq6MIB8l1CkIAQ4ZAcjJ8\n8IHVSbzL5aTLlHunHLuf3k1o4VCr4wgHSEezG9mxDVEywfbt8MUXMHp0ztvJvnJM5kyFAgrRsXpH\nvtvxnXWBsP9+8nRSFITX0BqefBJeeQVKlrQ6jXfqWacns7bOsjqGMJGpzUdKqQ7ABMAfmKK1fiub\n7boB3wCNtdbrs7hfmo9ErqZPh0mTYPXqrM84ErfuWvI1yr5Tlq2Dt1KuaDmr44hc2Kr5SCnlD0wC\nOgC1gYeVUrWy2K4oMARYbVYW4f1OnzZOQf3wQykIZsqfLz+da3bmm23fWB1FmMTM5qOmwB6t9QGt\ndRIwC+icxXavAm8C14A8VTQr2bEN0ZczvfCCMXlOYwcvuPXlfZUXWWV6qM5DzNpmXROSp+wnT2Vm\nUSgPHMq0npB2WwalVCOgvNZ6ftpN0kYk8mzpUpg7F8aMsTqJb2hTuQ17Tu/hwNkDVkcRJshn4nPn\n+AGvlPIDxgOZpz3J9kghJiaG8PBwAIKDg4mIiCAyMhL4p0q7ez2dVa/vCeuRkZGmPv/ly9CrVxyD\nB0NwcN4en85O+8tu61m9fyuWr6BpYlO+3/E9z975rNvzpd9mh/2TeT1zNqvyxMXFMW3aNICMz8u8\nMq2jWSl1B/CK1rpD2vrzQGp6Z7NSKgjYA1xMe0gZ4DTQ6cbOZuloFtkZNgyOHYPPZSpht5r711zG\nrRxHXEyc1VFEDmzV0Qz8AVRXSoUrpQKBh4Cf0u/UWp/TWodqrStrrStjdDTfVBDs6sb/DuzA1zKt\nWgVffgkTJ+b9sb62r5yVXaY2lduw4egGTl0+5d5AeNZ+8kSmFQWtdTLwFLAQ2A58pbXeoZQarZTq\nZNbrCt9w6ZIx3/L770OJElan8T0FAwoSVTmK+bvn576x8CgyzIXwSE88AZcvw4wZVifxXdM2TmPu\nX3P5tse3VkcR2bBb85EQppgzBxYuNI4ShHXurX4vv+z7xRbDaQvXkaLgJDu2IfpCpmPHYOBA4wgh\nKMj55/GFfeUKOWUKLRxK/dL1+W3/b+4LhOftJ08jRUF4jNRUiIkxlhYtrE4jAO6vcT9z/5prdQzh\nQtKnIDzGm28aTUdxcRAQYHUaAbD52Ga6ftWVPc/ssTqKyIL0KQivtXw5vPsuzJolBcFO6pWqx8XE\ni+w9vdfqKMJFpCg4yY5tiN6a6cQJ6NULpk41pth0BW/dV66WWyalFNFVo/ll3y/uCYRn7idPIkVB\n2FpKCvTpYxQFmW/ZnqKrRrNo7yKrYwgXkT4FYWsjR8IffxinoOYzc6Qu4bSjF49Sc1JNTo44ST4/\neZPsRPoUhFf56iv4+mvjqxQE+ypTpAzhweGsPbzW6ijCBaQoOMmObYjelGnTJnjqKZg925ypNb1p\nX5nJ0UzRVaP5Za97+hU8eT95AikKwnZOnIAHHoD33oOGDa1OIxwRGR7J0oNLrY4hXED6FIStXL0K\nUVHGIpPmeI6zV88S9m4Yp0acItA/0Oo4Io30KQiPpjX07w8VK8L//md1GpEXwQWCqRZSjT+P/Gl1\nFKEDpdwAABMESURBVHGLpCg4yY5tiJ6e6eWX4cABiI0FP5N/Mz19X7lLXjK1qNiC5fHLzQuTxtP3\nk91JURC2MHUqfPYZ/PgjFCxodRrhjJaVWrqlKAhzSZ+CsNxPP8GgQbB0KdSoYXUa4axjF49R84Oa\nnBpxCj8l/2/agfQpCI+zfDk8/rgx0J0UBM9WukhpShUuxZZjW6yOIm6BFAUn2bEN0dMybdkCDz4I\nX3wBjRu7LxN43r6ySl4z3RV2F6sTVpsTJo037Cc7k6IgLPHXX9ChA0ycCG3bWp1GuEqTck1Yd2Sd\n1THELZA+BeF2+/ZBZKRx2mlMjNVphCutO7yOx+c8zqYnNlkdRSB9CsIDxMdDmzbwwgtSELxR/dL1\n2XN6D5eTLlsdRThJioKT7NiGaPdMCQlGQRgyBJ54wrpMYP99ZRd5zZQ/X35qh9Zmw98bzAmEd+wn\nO5OiINziwAFo1co49XToUKvTCDNJv4Jnkz4FYbo9e4zO5OHD4emnrU4jzBa7IZbF+xfzedfPrY7i\n86RPQdjOzp3QurXRhyAFwTc0Kd+EP478YXUM4SQpCk6yYxui3TJt2ADNm8cxZgwMHGh1muvZbV+B\n92S6rcRtHDp3iCtJV1wfCO/ZT3YlRUGY4rffoH17o1O5Xz+r0wh3CvAPoHqJ6mw/sd3qKMIJ0qcg\nXO7bb2HwYGMqzchIq9MIK/Se3ZvoKtH0i5D/CKxkyz4FpVQHpdROpdRupdTILO4fppTappTapJRa\nrJSqaHYmYZ7Jk42jg0WLpCD4snql6rHluIyB5IlMLQpKKX9gEtABqA08rJSqdcNm64HbtdYNgG+B\nsWZmchU7tiFamSklBYYNM4atWL4cIiKsz5QTO+bypkx1S9U1rSh4036yI7OPFJoCe7TWB7TWScAs\noHPmDbTWcVrrq2mra4AKJmcSLnbxInTtanQsr1oFVapYnUhYrV6pejJaqocytU9BKfUg0F5rPSBt\nvQ/QTGud5cmJSqlJwBGt9es33C59CjZ1+DB06gQNGsDHH0OgTM8rAK01QW8GsX/IfkoUKmF1HJ9l\nxz4Fhz/J0wpGI+Bt8+IIV1q1Cpo1g+7djZnTpCCIdEopaoXWYtepXVZHEXmUz+TnPwyEZVoPAxJu\n3Egp1RZ4AWiZ1sx0k5iYGMLDwwEIDg4mIiKCyLSezPT2PHeub9y4kaFp4zVY8fpZraff5o7XmzcP\npk+PZOpUKFIkjqVLs97+xmzu3B/y/t36+q28f9VCqrHn9B4S9ya6NN+ECRMs//u/cd0uv09xcXFM\nmzYNIOPzMs+01qYtGEVnLxAOBAIbgVo3bNMQ2ANUzeF5tN0sWbLE6gg3cUema9e0fuIJrWvW1Hrn\nTntkcoYdc3lbppeXvKz/8+t/XBcmjbftJzOlfXbm6XPb9OsUlFIdgQmAP/B/Wus3lFKjgXVa67lK\nqV+AusDRtIcc1Fp3ueE5tNk5Re7i46FnTwgNhZkzoVgxqxMJO/ts82fM2z2PL7t9aXUUn+VMn4LZ\nzUdorRcAC2647eVM37czO4O4dXPnGnMpDx9uLH5yLbzIRfWQ6uw+tdvqGCKP5E/bSZnbWu3CjExJ\nSTBihHGF8nffwXPP5a0g2HE/gT1zeVum9D4FVx/le9t+shvTjxSE59q7F/r0gZAQWL8eSpa0OpHw\nJCUKlcBP+XHy8klCC4daHUc4SMY+EjfRGmJjYeRI+M9/jCGvpblIOKPJp02Y2GEid4bdaXUUn2TL\nPgXhWU6ehAEDYN8+WLIE6ta1OpHwZBWDKpJw/qaz0IWNyf9/TrJjG+KtZvr+e6hfH6pVg7VrXVMQ\n7LifwJ65vDFThaIVXF4UvHE/2YkcKQhOnICnnjLGLvrqK2jRwupEwluEBYXJkYKHkT4FH6a1UQSG\nDoW+fWH0aChY0OpUwpvM2jqL2Ttm83X3r62O4pOkT0E4bN8+4+ggPh5++gmaNrU6kfBGFYq5vvlI\nmEv6FJxkxzZERzJduwZjxhhFoFUr41RTMwuCHfcT2DOXN2Yyoyh4436yEzlS8CGLFxtHBzVrwp9/\nQqVKVicS3q5c0XIcvXiUlNQU/P38rY4jHCB9Cj7gr7/g3/+G7dth/Hi4/36rEwlfUurtUmx5cgul\ni5S2OorPseN8CsJCZ88aU2Q2bw4tW8K2bVIQxP+3d/ZBVpV1HP983ZVVDAFnHVLBwYIYVHIETdSa\nsIzwDR1NQTNTmzIrXypnSKdMU0fRknQcM9QkX8m3EOxFEF01AUtxAVEDlDUwZUDbFURe99cfz3Ov\nl+u+3F323HP28vvMnLnn5bnnfvc5Z8/vPM/veX6/8lPbs5Y169ekLcMpETcKnSSLfYg5TRs2wKRJ\nMGQIfPhhMAaXXAI1NelpyhpZ1FWpmrraKFRqPWUF9ylUEFu3wp13hqGlw4fD7Nk+I9lJH28pdC/c\np1ABbN0aIphefjl8+tNw7bVwuIeacTLC92Z8jxF7jeC8Q85LW8oOh89T2MHYsiVMPrv6aujdG266\nCUaPBnXoFnCcZKntWcvq9avTluGUiPsUOkmafYibNoUopkOHwm23wc03w9y5UFNTlzmDkNW+1izq\nqlRN7lPoXnhLoRvR2AiTJwcjMGQI3H57mICWNUPgOIXU9qxl/jvz05bhlIj7FLoBb70VuoamTIFj\njw3pMA8+OG1VjlMa016fxpT6KUwbPy1tKTscPk+hgmhuhpkz4aSTggGoqoIFC+Dee90gON2LT/X4\nFOs2rUtbhlMibhQ6SVJ9iO+/H2YdDxkSciMfd1wIWnfDDTBgQDqatocsaoJs6qpUTV1tFCq1nrKC\n+xQywJYtMGtWcB7PnAnHHw933w0jR7q/wOn+eEuhe+E+hRRZvBjuuScs/fvDOefAuHHQt2/ayhyn\n62hobGDUlFE0XNyQtpQdDp+n0A1YujTMLZg6FZqaYPz40Do44IC0lTlOMvTq0ctbCt0I9yl0klL7\nEM3glVfgmmtgxIiQ6nLVqjC/4K23gq+gqwxCFvs1s6gJsqmrUjW5T6F74S2FBNi8GZ5/PmQ0e+yx\n4DM48US4/vowr6Daa93ZgehR1YPNzZs9p0I3wX0KXYAZvPEGPPFE6Aqqq4PBg0OY6rFj4aCD3GHs\n7NjUXF3D2kvX0qOqR9pSdijcp1AmzGDZMnjuOXj2WXjmmRB6YvTo4CO44w7Yc8+0VTpOdqjeqZot\nzVvcKHQDEvUpSBoj6XVJSyVNaOF4jaQ/xePzJGUyQeS6deHB/+tfw6mnwt57wxFH1PHkk3DYYfD4\n47ByZRhSevrp6RmELPZrZlETZFNXJWvKGYWuoJLrKQskZhQkVQG3AGOA/YHTJQ0tKvYd4D0zGwxM\nAiYmpacUzODdd0Mu40mTwhDRAw+Efv1gwoTgGB47NgSfu+yyeu6/H84/PziKs9A9VF9fn7aET5BF\nTZBNXZWsqSuNQiXXUxZIsvvoC8AyM2sAkDQVOBF4raDMWOCXcf0RghFJnI8+gjffDF1AS5eGzyVL\nYNGiEF5i2LCwHH44XHBBWN95523P0dTUWA6pHaKx0TWVShZ1VbKmrjQKlVxPWSBJo7APsKJgeyVw\nWGtlzGyLpCZJe5jZ+x39MTP44ANYs+bjZfXqMPxz5cqwrFgRPhsbYeBAGDQoOISHDYOTTw6tgr32\nysZbv+NUEl1pFJxkSdIodOlwoa99LQztzC0bNoS+/tyyfj307Bn682trt/0cNAhGjQqxg/r3D91B\nVds5Mq6hoaEr/qwuxTWVThZ1VbKm6p2q2bx1c5ecq5LrKQskNiRV0kjgCjMbE7cvBZrNbGJBmb/H\nMvMkVQPvmNkn3LSSsjse1XEcJ8NkaUjqi8BgSQOB/wLjgNOLykwHvg3MA74BzG7pRB39oxzHcZzO\nkZhRiD6CHwFPAFXAnWb2mqQrgRfNbAZwJ3CPpKXAe8D4pPQ4juM47dMtZjQ7juM45SFTAfGyONmt\nBE0/kbRY0gJJT0raN21NBeVOkdQsaXgWNEk6LdbVK5LuS1uTpH0lPS1pfrx+x5RB0x8krZK0qI0y\nN0fNCyQlnmevPU2Svhm1LJT0vKTPp62poNyhkrZIOjkLmiSNkvRyvMfrktZUii5JvSXNkFQfdZ3d\n5gnNLBMLoYtpGTAQ2BmoB4YWlfkBcGtcHwdMzYCmUcAucf37WdAUy/UCngXmAMPT1gQMBuYDveN2\nbQY0TQbOi+tDgeVJaoq/8yXgYGBRK8ePBf4a1w8D5mVA0+EF121MFjQVXOOngMeBU9LWBPQBFgP9\n43ai93gHdF0GXJvTROiqr27tfFlqKeQnu5nZZiA32a2QscAf4/ojwFfT1mRmdWa2IW6+APRPW1Pk\nKuA6YCOQtKO+FE3fBW4xsyYAM1uTAU3NQO+43gd4O2FNmNlzwP/aKJK/x83sBaCPpH5pajKzubnr\nRnnu8VLqCeAC4GFgddJ6oCRNZwCPmNnKWD7pe7xUXc3A7nF9d0IUiVYnjWTJKLQ02W2f1srEP6pJ\n0h4payrkO8BfE9QDJWiK3UX7mFlOS9KOo1LqaTAwRNI/JM2V9PUMaLoCOFPSCuAvhIdM2rSkO/GH\ncAcoxz3eLpL2IRj538VdWXCODgb2iF2SL0r6VtqCIrcA+0v6L7AAuKitwlmKkpqFi1pMyZoknQkM\nB36cnBygHU2SdgJuJAz1ze9OVFFp9VQNDAK+DAwAnpU0rOANNA1NZwB3mdmkOK/mXiALOfCKr1cm\n/jckHQWcCxyZthbgt8DPzMwkieTv8VLYmfAM+CrQE5graZ6ZLU1XFmOA+WZ2lKTPArMkHWRma1sq\nnKWWwtuEh0WOAYS3pOIy+wLEyW69rRMhMbpYE5KOJvTbjY1dFUnSnqZehAdbnaTlwEhgesLO5lLq\naSUww8y2WoiHtYRgJNLUdC7wIICZzQN2kVSboKZSKNbdnzJ0a7VHdC7fTrjH2+vWKQcjgKnxHj8F\nuFXS2JQ1rQBmmtlHZvYewad3UMqaAM4GHgUwszeA5cCQ1gpnySjkJ7tJ6kFwJE8vKpOb7AZtTHYr\np6Y4OuQ24IQy9SG2qcnMmsxsTzPbz8z2I0wMPMHM5qelKTKN4JQnPng/B7yZsqb/AEdHTUMJAwbK\n0g/cBtOBsyAfFaDRzFalKSiOqHsUONPMlqWpJYeZfabgHn8YON/Miq9vuXkM+KKkKkk9CQMFXk1Z\nE2x7n/cjGITW//fK4R3vgBf9GODfhFEjl8Z9VxIeagA1hDe7pYSH3cAUNR0f12cB7wAvx2Va2vVU\nVPZpEh59VKom4DeE0RkLgdPS1kQYcfQPwsikl4Gjy6DpAcIM/02EN8tzgfOIo6BimVui5gVlunZt\nagLuIIxYyd3j/0xbU1HZu4CTs6AJuCTe44uAC5PWVOL124swiXhh1HVGW+fzyWuO4zhOnix1HzmO\n4zgp40bBcRzHyeNGwXEcx8njRsFxHMfJ40bBcRzHyeNGwXEcx8njRsFxHMfJ40bBcUpA0rr4ubek\nh7bjPBdL2rWD36mLeSGO78B3donx8zcmHDTSqTB88prjFCGpysy2Fu1ba2a9uuDcy4FDLMTGKfU7\nTwM/tU6EKom/N8KSjRHmVBDeUnAyS8yqtUAh495uMWvU/i2UOyuWq5d0d9w3UNJT+jgj3oB29k+R\ndJukecBESfvFEN8LJV1d8FsDcxmuJJ0t6VFJf5O0RNLEgnK3SvpX1HxF3HchsDfwtKTZcd9oSXMk\nvSTpQUm7tVYdBeeuk3RjPP+rkg6JOpZIumq7Kt1xyhGbwxdfOrsQkgXdQIgHNKGF4wcQ4hvtEbf7\nxM8ZwLfi+jnAn9vZP4UQjC7Xep5OCAAHIePf2rg+kJjhihB98g1CZNoaoIGQxwKgb/ysIsSfOjBu\nLy/QWgs8A+watycAv2jhb9wmflXczmXSupAQRbUf0IMQ+6ZvQdn87/niSymLtxScrPMrYDRwCHB9\nC8e/AjxosXvEzBrj/pHA/XH9XuCL7ew34CEzy/WnHkEINJYr1xqzzWytmW0kRMTM5Q0fJ+klQgrS\nA4BPtHCilv2BOZJeJkRHLTXHdy4i6CvAYjNbZWabCNEvE88T7lQuWUqy4zgtUQvsRnjj3hVYX3Tc\naD3BSkf3F5+7FDYWrG8FqiXtB/yU4DtoknQXsEsr359lZmdsx+82F2loJtSV43QKbyk4Wef3wM8J\nb/cTWzj+FHBqboSNpL5x/xxgfFz/JiHhSVv7i3m+qFypiNCd9CHwQYxff0zB8bV8nC/3BeDImA2L\n6DcZ3IHfakuD43QKNwpOZpF0FrDRzKYC1wGHShpVWMbMXgWuAZ6RVE/I2QAh3/I5khYQHuoXtbMf\ntk17eRHwQ0kLCc5ha6Gc8clUmWZmCwl5B14H7iPkbMgxGfi7pNlmtprgl3gg6plDGxmxWqFFDR08\nh+Pk8SGpjpNx4pDUS8zspU5814ekOh3CWwqOk33eB6Z0ZvIawW/YnJgyp+LwloLjOI6Tx1sKjuM4\nTh43Co7jOE4eNwqO4zhOHjcKjuM4Th43Co7jOE6e/wMrXTBLso5a9AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x105b65c50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# draw an example trajectory for a 1 GeV particle in a 3.8 Tesla field\n",
    "phiValues = numpy.linspace(0, numpy.pi, 101)\n",
    "\n",
    "# make sure we have a 'square' plot\n",
    "pylab.figure(figsize=(6,6))\n",
    "\n",
    "pylab.plot(\n",
    "  [ xpos.evalf(subs = {phi : phiVal, B : 3.8, p : 1}) for phiVal in phiValues ],\n",
    "  [ ypos.evalf(subs = {phi : phiVal, B : 3.8, p : 1}) for phiVal in phiValues ],\n",
    "    label = 'charged particle'\n",
    ")\n",
    "pylab.grid()\n",
    "pylab.xlabel('x coordinate [m]')\n",
    "pylab.ylabel('y coordinate [m]')\n",
    "\n",
    "# get the plot limits to make the aspect ratio 1:1\n",
    "maxSize = max(pylab.xlim()[1], pylab.ylim()[1])\n",
    "\n",
    "# draw the ECAL inner surface\n",
    "pylab.plot(\n",
    "    [ rEcal * math.cos(phiVal) for phiVal in numpy.linspace(0, math.pi, 101)],\n",
    "    [ rEcal * math.sin(phiVal) for phiVal in numpy.linspace(0, math.pi, 101)],\n",
    "    label = 'ECAL surface'\n",
    "          )\n",
    "\n",
    "# make sure the aspect ratio is 1:1\n",
    "pylab.xlim((0, maxSize)); pylab.ylim((0, maxSize))\n",
    "\n",
    "pylab.legend()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq8AAAAqCAYAAACdg7LfAAAABHNCSVQICAgIfAhkiAAAE4RJREFU\neJztnXm4XkV9xz/JTUhCNkISCQmRSxaWoGgCgWxASA1VChVQdpc81EKIgsiiAamGWCAVgwSKKEt5\nMShYUBGwodIqtVUQUsSCtoEHNQJNAxFRQSIQbv/4zvHMOfcsc971Jvl9nud97n1n5szynt/Mb5bf\nzIBhGIZhGIZhbCX063QGDMMwDMMwjEz6Ax8Ghrjvn+1gXgzDMAzDMAyjkKOAie7/rwP7dzAvfYb+\nnc6AYRiGYRiGkclk4CT3/1PEHVnDMAzDMAyjwxwITEAmjeOAWdtRnHkMAoa7/9cA41uYlmEYhmEY\nRp9kOHAH8OaKfqFhDwI+BiwD7gMOSfl/ALgI+BBwrudeA3qALcADwFTP7wDgKuCDwPXAlEC/ojiL\n8lJvnPWWvYy5wNIK4Q3DMAzDMLYJPoQ6Vj1AdwW/0HiGAZd5348H/oBmKgEWASvc/93Aq8AY930Z\nmskcl0prB2A9sKv7PhN4OMCvKM6ivNQbZyNlL2IE6vAahmEYhmFstxR1UEM6r3lh9wPeQLaaoI5X\nD+rI7QA8D+zuhZ/k/b8sJ42FwE+97/2Al4A9SvyK4izKS71xNlL2Ik4HBrrPOwKfaQVls8ptY0Cn\nEjYMwzCMChwEzAFGoiXUzwDf72iOjCweA2YDP3ffow1GTzr3MagDNw/tnL/XCzsU+CvgFdSB/Bzq\nRHYDv/bS6AF+A7wFzX7m+f2iIM6ivBSlVxRnvWUfijqoByOThGHADOCf3DOXA5cAXS5MJxgGHA1c\n4L4fj/I/FXi2Q3nKZWCnM2AY2xjHdjoDhtFkmiXTFwPTc/zKlmON6rRq5jXNamCl+/9kFz7qgA0D\nXgR2c99nE58/vxBYh05FugB1lHyedPEV+RXFWZSXeuNME1r29yP74QeI69M04L8y4uwURbPKbafo\nqKzj0A9qGEZz2BsZ6xvGtkIzZfoydAD75Ay/ScAnPL970aHtc5uUttEaTgU2AOe57791f9e6vy+h\n2csj3feHUYcI1FncE3irey59qdIwYFOJX1mceXmpN06fKmW/C3gNbQq7y/nvBoym71A0q9x28swG\nFqDRwVlNTGsWmhZ/BTU6Q9A0+GNNfj403ArgW8D/AH9Etix/iUZKz6TC7olsPF5xYV9x339fkN/z\n0Mz1ZSn3vZzbeiT8Y4DzgY2pcJOATwOb0Y7Goajx/r8GygGwI/AgGkVlEZpuNDrtQUd57OjyEvo+\nO82bkb3OLQFhrwTO9r53odH0hWj2J5STgVtz/G5Es07TXZw/RL8/yFZqHPBN4G+R/LUCq6PC6mg4\nzZTpzcAS4Muobm7x4uqU4ix7F0WEyGTou9waOcL9/TgwGL3vR5E8dnnhelBfZBbwL6i+bSY+Huo1\nVHdO954ZAIxCdfT1Ar+iOIvy8rM646y37L8F5iNZe935HY7sSpvJfOAbqM5dUfHZHuBH3velLo4f\nBzwbqm8b0rUj0ahix4AMhTIdjSYGe25fAH4HvL2Jz1dJpyf12QL8TUba04BfIVsrkBD+gnyDbZBN\ny8sZYUYiZfU+z+1C4HHUmEfsgUZ4/nLc+4Cf0HvAEVoOiHdM9uT4h6Y7GriTeCdm9Ow6eh9D0hcZ\niToOIWYx04CPZLjviZRsFR4hvuIvi270btKdKdC7exX4doX0qpj9WB0VVker0QqZvo3yVT9/ObZV\nlL2LIkJksooMlVGPiV9Vs4E/B94WGPYQZBMa7cg/Gg0+AL7r4gIYizpuE9Bg69NeHIuBH7j/BwD/\nS3wk12FI9sr8iuIsyksjcdZTdlycl3p+TxBuuhHKSeh9fanBeE5FqyTp2eksQvVtw7r270j2fJvB\n59EPdpzndqRzu6qJz1dJZz1wHWrgV5JdKQegxt7/PSaiHYMfLcjvdS7NZSn3S4DnSDZMO6MR2xme\n253ACyTNOoagUd6iOsqxD1IQNTSyy2uMQ9M9Hzgz4/kLkfz0da5GxvAhLEZG+lksI9y05kBUgYt4\nP3o3ebtJ/9P5jw1M83rCbQKtjgqro+G0SqZnoZmvPKooznoIfRd5hMpkFRkqo0pdPwW4FpXrNpId\nhiK/bwKfDAg7CQ1G04O2Ec5/InAD2rV+HeoYRhyFVkQ+gWbu3+T5/RnqdH0QuAl1akL8iuIsyks9\ncTZS9n9Fg7xTgFVktxPNYAqNbdY/AtVB0CREd0n4UH3bkK4diipTs+0sPoBGGId7btEIYEXmE/U9\nXyWd+wPSPRXNDOwUEDbiWOBEshXjE8DdGc88hkZkoNmd14CHMsKtA+5Jud1fIW8QH66cpkq616KG\nLM05VF+KaDeT0JJPKKvIV5JjkEF9V45/Op53lYS5Hi0vZs1k9UczAb8nOQNYRI3wkbvVUWF1NJxW\nyXQX6gCkbQihuuJslBrVO68hMllVhsqo0frfwmgdA9GgOUSXdJKiWeUsqujbhnTt8cRT463mcmTb\nUe/oIvT5vHD3B6RxH8UzAGmGoVEV9FaMw53bNRnP/TOxMfc4Fy7rCJiHkID73F8hf5DfGFdJ93QX\n9qvIDghkU/cI9dmGtZOVaCQdypUl/msoP3evC1W8stHuOuDfc/wWo9/8tJI4fGo0ptCsjsZYHU3S\napm+F80e+1RVnM2gRvXOa4hMVpWhMmpY53VrZg7x4LhVDEL2/PMIMwVLUzarnEUVfVtJ16YbnoXI\nqL7VTEbLImci+55WPV8Ubge0mWEM8fEPS9HMC2gEMA81JIegmaJhqIG4mGwj5QvItu2C+GDi32X4\nvYwEYBBaWvoDSZvAiPFoeW0AsVF3WTlCqZJuDd2uchJa9jgfeCdaUvKP9ng3kqm3oKWQuejKvVHI\n7u8sNEMBqkyfQkeHbEazMjcQ78w8CL3HjS4fOyH7oHVeescgZfYb9FuOQe/xw16YI5CNZR4DkczM\nRr/tKGRvsxZ1arakwv8Abb4pGl0uQAr89YIwu6BlqX9MuQ9y+T8XNQJV7WzrxepoEqujSVot0z8l\nqWAnodnI4alwIwvS7wShMllVhoxtl/3QBqqxaCVjTQNxFenAqagzeThwM2p3jwGWI9v4VchM5jC0\nCj8T1dEHXNw/p7ijmkWRvm2qrl1Lcoddszka2T88gRrzqjZLoc+HhHuKeMcqyFB+I/F1b2PQqOK/\n0QxBxHy0zLVvKr63Ex/eC71ndeY4t4sz8rLa+e3ivv8DMuT38z0OKb4ekjY7ZeVIUyN/JqFKusPR\n7Eg0+rrHyz9IGKOR1MNIYS70/B8i/n3moA7DfM//y8hWEGRftI6kXdw+wC+Jz4Xch94j1xNReSMm\nkDx4Os3uqKKeiTrPZxDPBl5D9k7JI0ka62dxE+VH+hxHPFO2wn2uQLs776O+TTY1qs/GWB21Ohrh\n19E0rZbpi5Ai7TQ1qs28VpHJKu8yJJ/dFcIb2x4hOhDUz/PdRiPZ/DZa2Yi4CunYeinSt03XtZuA\nE+rNaQV2AL7nMlKPfW3o80Xh0rYTXWgkHG0a2QU1IJvpba/1DEmbpP5oJOPbbaUV40EZbhG3Ob/x\n7nu00zAaSAxAI7NHXDi/LGXlSFMjvzGuku5SZM/2F2hE1oPOtIvs1BYC70UN8ybiXZQRX0VLa/1Q\nxzS94/jzaMYmOlsvaxPhKnT8CEhuHyc5MzOU5I7kuUipZDEa7Qb2Nyh8gVixTEBlTC9tz6Z4iW+w\ny2NZJ/BqZBuYNRNzCWpcFmb4FVGjfoVmdTTJ9lxH07RDppeQX1fbSY1qndcqMlnlXYbks7tCeGPb\nI0QHgsyYaim3X6L65tfp06i2QThNnr5tia59jfgYB5+3oeWORwM/NxUUKGK+y+DtAWEbeb5KOuuJ\nZ/sGuucezwj3IGqUB7nvZ5DcNQi9leCkDLeIu52fL3Sj0NEZV6AGbXc0YnqFcqXhlyNNjeLGOCTd\nj5HsbA5Fnc03iJd+d0VKa1+X3qGpdB5EAjzL+S/Pyc8Jzv+dGX5LnN/+qFPxHLJLXI1GdOlKdyz5\nszk3kFy27kdyuWOIS+vE1HN7IVnI4z3kL1P7/IT8UeVAdKD1BrIN1m8mux6+gDofWX77B+RpPlZH\nI7bnOpqmHTK9mOpmFRHN1Fc1qnVeq8gkVJehVtV1Y+snRAdCfuc1fV7zIiTL3XXmJ0/ftkTXvkr5\n7tF62JvemwSiq8XeQLNrzXg+NNz3yd5E8AwaMUc8R7ZB/b+5OHdFyzxXZ4RJK8GhLg+rMsJ+F9mo\nlPE0Ok4jIrQcPjWqb0Dw0+1HfMdzmmgDxj6e2xLUEPsN9o5ICNcQ7/rO24gUHbKetSnqr51ftFqw\nBzra5Gni2Y+TvPDHk23TPdDl8VzP7a0kbWVnuDjTS6X7km0jGXEH+cd/ROyE7HuKdvVvovoouEZ4\nw2N11OpoXh1N0w6Z/iQyZeg0Naq/ixCZLCItQyHUsJlXo1wHQn7nNe22iMY6r1n6tmm6Nn097Ivo\nTMNmMgKNgh8hee1fZIzbj+LjIUKfr5LOdHob/oPslfzZkAfJblgHoYb9edSp2hud2Rd9ohmPE933\nY9GGjx+TtH2LmEK8/J3HWHQ4sj87FVqORkinOxYppqybbb6EhMtfJjwUCbA/YlqAlm9vRUflQLwj\nOs0G9zfL/itaVtuIOkR/REtwE9Hh0l9zeYqU8vM56ezs8uxv/jqMpP3QKWgm6AGSjCJ/KWOky0fW\nLIzPwagu/keO/1RU1udQGZqN1VGro0V11KddMh0tLW6NhMhkHlky1Ex67LNNfiBMB7aTLH3bKl3L\nj8g+2LoRBiNzhCdJdowPQD96enQ9laStUOjzVdK5AykPn+kunH/zzUnINs3PTz/UyU/voPXpdnEt\nS7kvB54luRw02YVd4rmdjV7Sbp7beWi2xlc6oeXwqRELe5qQdPuhqwvnZTw/3IX1K8oGev8Oa1An\nYYD7PI2OIkpzDLJ/eZnsTTS3u2cHoFHiOSn/LvSuok0qe5It/FGZ/BHqtd7/+6FluTn05hhks5nF\nqcjusIzL0YxfXgf+dvTOzsjxz6NG2KjZ6miM1VHh11Gfdsn0PSQ313WKGvnvAnrXBQiXyVAZCs1n\nd8Vntgf6oz7Nx91nW2YR5ToQ2jfzmqVvW6VrWZmKqFlciiqqrxBWI5unAzy3Q1CDl+7IhD4fGm6m\nc48a8H5od/sPSTbq/VHv/zzP7Xg0ouhOF9JjKnrpl6Tcd0WC5N8UcSU6AN3fSHIREqZoN/J0NLuY\nfqGh5fCJNp5kXf8bmu4J6Cibbs9tBFJI7/Xc9nRp3eG5fQQ1zlM8t4VoeePdnttYYjuYU9wz6asu\nNxHbMS5Cs7j+jMfuJA8B74cU9Xh6swjNmAxB7z066/MwpMTT9pIRl5K/NPod4uOXilhL9n3zOyH7\noC0kryMMpUZ4w2N1VFgdza6jEe2Q6f7oHcwMSKfVFL2LvLoQKpOh7zKEGtZ5zeIo4pWUr7Nt2/8u\nolwHguzQ07v5n81wi8zy6r1OOk/fLqIJujY9qr6X8oNi6+FClOFb0bl1u6De9QySRvkbUQV/qs7n\nQ8M9jK5muwXZX4xADe1pJJfO3kA7da9Ajf5r6IXMIvsIiRHAt5BhMciuYwG6jvFO9CLnI4U5A9n3\n7YzsjP2zFFei5bnL0EzJSHS0UHoaPbQcb0IKdALxUS3r0dLfDcBXKqb7NVeWK106PWiEdw3JA9kP\ndf7XAV9Ev+cgpJQ2eOHuQ4pgGVK6G9xz0czLV9A94SuRUtuClMlhxArydZf+UtQR7kGK8j1eOj1I\n+R7syuBTQ7ufb0YG75NR53mDS+dFspmHrghMM86VNW9peBRqTHdCSuslYps40MzLaLQaciC6RrOV\nWB21OlpUR6F9Mr0/krFO2byGvou8uhAqk6Hv0qifyci++7PER9a1ui3tFGU68FDgc2iSYCbqVK5C\n56uPRwOsbnTE3SrifSbfQZsa8zZV55Gnb2s0V9cCaph+TfgdyYZRxC3k27x1ijmU25OdRfIe6zwm\nk98Qnk1ymbkT1LDZGKOYKnW0XTK9muQGDqOcGlbXsxhEbHO+huxVN6N1lOnbZujaP7EcjVIMo1Ge\npni3cae4k+SGoTR/HxjPNWQfLQdamknbOrabGWh3p2HkUaWOtkOm90CXJZjcVsPqejFzCbPVNppP\nkb5thq79E8PQ8kizTx0wti8im8KjOp2RDKagA9jzzuLMuts+zSHkm9hMofelC4bR16hSR9sh0wOR\n6Vo9964bRh4jkH2x0RmK9G2jurYXs4Fv0PsoLcMI4aPItqsH3bDRF3YNp5lH8vrGiN3RAe9FdKHG\nMO/4qE+hK0ANo69StY62Q6aXI/tjwwhhMGqHX0B27Evd50a00W6kC3c6GhgNJPu8cKP1ZOnbZuja\nTN4VELFhGL25m/JD/Q1ja6IdMp3eQGwYZQxGGyH3Srn/DN1YdgI613gT+Rd3GIZhGIZhGEZbeAfx\nRTcRXciW20wFDMMwDMMwjD7FCuKjy0Cz9zeiTUKduFXKaBN5G1YMwzAMwzD6MmvREUrfQzOuRyAT\ngXOIr502DMMwDMMwjI6zM+qg7pZyv4vwI5iMrRQ7TcAwDMMwjK2NBeikjGdS7r9C57oahmEYhmEY\nRp/hi+i6U59p6Oisk9ufHaOdmM2rYRiGYRhbC9NQ53Qxutr4Uec+GpgIrEI2sIZhGIZhGIZhGIZh\nGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhGIZhlPH/47VRKz2aEWgAAAAASUVORK5CYII=\n",
      "text/latex": [
       "$$\\left(- \\frac{3.33564095198 p}{B} \\cos{\\left (\\phi \\right )} + \\frac{3.33564095198 p}{B}\\right)^{2} - 1.69 + \\frac{11.126500560526 p^{2}}{B^{2}} \\sin^{2}{\\left (\\phi \\right )}$$"
      ],
      "text/plain": [
       "                                                  2                           \n",
       "~  -3.33564095198~(-1)~p~cos(~)   3.33564095198~p~           11.126500560526~p\n",
       "~- ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ + ~~~~~~~~~~~~~~~~  - 1.69 + ~~~~~~~~~~~~~~~~~\n",
       "~               B                        B       ~                        2   \n",
       "                                                                         B    \n",
       "\n",
       "2    2   \n",
       " ~sin (~)\n",
       "~~~~~~~~~\n",
       "         \n",
       "         "
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# the expression for which we want to find the root\n",
    "xpos * xpos + ypos * ypos - rEcal * rEcal"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# find the angle corresponding to the charged particles trajectory when\n",
    "# it hits the ECAL\n",
    "\n",
    "# this seems to give somewhat non-sensical expressions (with coefficients 10^22...) \n",
    "# \n",
    "# phiSolutions = solve(xpos * xpos + ypos * ypos - rEcal * rEcal , phi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# create a function with the expression to be solved\n",
    "func = sympy.lambdify([phi, B, p], xpos * xpos + ypos * ypos - rEcal * rEcal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.62070977555 -1.62070977555276\n"
     ]
    }
   ],
   "source": [
    "# test this function\n",
    "print func(0.1, 3.8, 3), (xpos * xpos + ypos * ypos - rEcal * rEcal).evalf(subs = { phi : 0.1, B: 3.8, p : 3})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.text.Text at 0x106ac6fd0>"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY0AAAEPCAYAAAC+35gCAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcFNW5//HP44KoRMe4r8HrboJBibhdFS8accMYo+IS\nMy6JMaLJdYnRmID5xUS9cYm7JiruoEENiCBEGYXEoKCDKK4xRFxC3FCRKMI8vz9ODU6Gnpnunuk+\ndXq+79erX3Z113R/pxr7mTpP1Slzd0RERIqxXOwAIiKSDhUNEREpmoqGiIgUTUVDRESKpqIhIiJF\nU9EQEZGiRS0aZnaTmc0zs1ltPD/AzD4ws6ez23nVzigiIp9bIfL73wxcCdzazjqPuvvgKuUREZF2\nRN3TcPcpwPsdrGbVyCIiIh3Le0/DgV3MrNHMHjSzbWMHEhHpzmIPT3XkKWATd19oZvsB9wNbRs4k\nItJt5bpouPtHLe6PN7NrzOyL7v5ey/XMTBNoiYiUwd1LagHkenjKzNY1M8vu9wesdcFo5u7J3oYN\nGxY9Q3fMrvzxb8of91aOqHsaZnYXsCewlpnNBYYBKwK4+/XAt4CTzWwxsBAYEitrJc2ZMyd2hLKl\nnB2UPzblT0/UouHuR3bw/NXA1VWKIyIiHcj18FR3UV9fHztC2VLODsofm/Knx8od18oTM/Na+D1E\npHa5Q1NTuK24Yuw0gZnhtdQI7y4aGhpiRyhbytlB+WNLOf+iRXDrrQ2MHw833AA/+xmccAIMHgy7\n7AJbbw0bbgirrx6KxHLLQY8e4fmU5fqQWxGR2Jqa4JVXYMYMmDUr3GbPhtdfhzXXhK98BTbeONx2\n2gnWWQfWWis894UvQK9esMoqsMIKoXCkTsNTIiItfPwxPP44TJkCf/kLTJ8e9hb69YM+fWC77WDb\nbeG//ivsOaSsnOEpFQ0R6dYWL4YnnoBJk2DiRJg5E/r2hT32gN12gx13DHsPtUg9jUSlPK6bcnZQ\n/thi5f/oI7jnHvj2t2G99eDkk2HBAhg+HN5+G6ZOhV/9Cg44oP2Ckfr2L4d6GiLSLSxYAGPHwt13\nw8MPw667hqb0BRfAJpvETpcODU+JSM1asiQUiNtuCwVj113h8MPh4INhjTVip4tPPQ0REeAf/4Cb\nbgq3ddeFY4+FIUNqtzdRLvU0EpXyuGjK2UH5Y+vK/E1NMGECHHhgONLp/fdh3Lhw9NNpp1WmYKS+\n/cuhnoaIJO3jj+Hmm+GKK2DVVeHUU0OTe+WVYyerTRqeEpEkzZsHV14J118fDo89/fTQszBdILpo\nGp4SkZr32mthb2KbbeC998KJeKNHh3MqVDAqT0UjB1IeF005Oyh/bKXkf+01OOkk2H77MPQ0ezZc\ncw1svnnl8nUk9e1fDhUNEcm1N9+EU04JxWLNNeGll+Dii8NJeVJ96mmISC7Nnx+Kw/XXw3HHwY9/\nrENmu5p6GiKSvEWL4PLLYcstQ7O7sRF+8xsVjLxQ0ciBlMdFU84Oyh9by/zucP/98OUvh8kDJ0+G\nG28MU47nVerbvxw6T0NEops9OxwRNW8eXH01fP3rsRNJW9TTEJFoPvwwzCx7++3hyncnnxwuViTV\noZ6GiCTBPcw2u802oXA891zY01DByD8VjRxIeVw05eyg/DH87W+w337wy1/CT37SwO9/D2uvHTtV\neVLc/p0VtWiY2U1mNs/MZrWzzhVm9rKZzTSz7auZT0S6zuLF4SionXaCgQPDNbf79ImdSkoVtadh\nZrsDC4Bb3X2Zfz5mtj8w1N33N7OdgN+6+84F1lNPQyTHZs6EE06Aurpw3sVmm8VOJJBgT8PdpwDv\nt7PKYOCWbN1pQJ2ZrVuNbCLSeZ99BuefD/vsAz/4QTiUVgUjbXnvaWwIzG2x/DqwUaQsFZPyuGjK\n2UH5K2nWrDAUNW0aPP00HH/8shMK5jl/MVLPX44UjlVovetUcByqvr6e3r17A1BXV0ffvn0ZMGAA\n8PkHm9flxsbGXOXRspY7s/zwww3ccw+MHj2Aiy6CTTdt4OWXYcMN85GvOy83NDQwYsQIgKXfl6WK\nfp6GmfUGxrbR07gOaHD3kdnyC8Ce7j6v1XrqaYjkwJw54dKqyy0HI0ZAmd9LUiXJ9TSKMAY4FsDM\ndgbmty4YIpIPd9wB/fvD4MHwyCMqGLUq9iG3dwF/AbYys7lmdryZnWRmJwG4+4PAq2b2CnA98IOI\ncSumefcxRSlnB+XvCh9+CMccE867eOghOPPMsKdRjDzk74zU85cjak/D3Y8sYp2h1cgiIqWbPh2G\nDIG99w7nXayySuxEUmnRexpdQT0Nkepyh8sugwsvDFfP+9a3YieScpTT00jh6CkRyZH33oP6+jAj\n7bRpsOmmsRNJNeW9Ed4tpDwumnJ2UP5SPfkk9OsXTtCbMqXzBUPbPz0qGiLSIfcwDLX//mH+qMsu\ngx49YqeSGNTTEJF2LVwI3/9+uOzq6NGwxRaxE0lXqcXzNEQkoldfhV13hSVL4PHHVTBERSMXUh4X\nTTk7KH97Jk6EXXYJs9PefjusumrXv4e2f3p09JSI/Ad3+L//C32Le+6BPfaInUjyRD0NEVlq4UI4\n8UR46SW47z7YeOPYiaSS1NMQkbK9/nrYqzALh9OqYEghKho5kPK4aMrZQfmbTZsWrn1x2GGhf7Hy\nyl3ysh3S9k+Pehoi3dxdd8EPfwg33ggHHRQ7jeSdehoi3VRTU7gU6y23wJgxsN12sRNJtWnuKREp\nyiefhPmj5swJQ1Prrhs7kaRCPY0cSHlcNOXs0D3zv/02DBwYDq2dPDluweiO2z91Khoi3ciLL4YT\n9vbcM/QyqtXwltqhnoZINzF1arjuxQUXhLO8RdTTEJGC7rkHTjklHE779a/HTiMp0/BUDqQ8Lppy\nduge+S+7DE4/PcwllbeC0R22f63RnoZIjWpqgrPOggkT4M9/hk02iZ1IaoF6GiI16NNPwyG1b7wB\nf/wjrLFG7ESSR5p7SkT48MNwhb1Fi8KQlAqGdCUVjRxIeVw05exQe/nnzYO99oItt4S774aePePk\nKlatbf/uIGrRMLNBZvaCmb1sZmcXeL7ezN42s6ez2/Excoqk4NVXYbfdYPDgcD3v5ZePnUhqUbSe\nhpktD7wI7A28ATwJHOnuz7dY5ztAP3c/rYPXUk9DurVnnoH99oPzzoOTT46dRlKRWk+jP/CKu89x\n98+AkcDBrdax7CYibZg6FfbZJxxaq4IhlRazaGwIzG2x/Hr2WEsOHGpmM83sHjPbqGrpqijlcdGU\ns0P6+S+8sIFDDoHbboPDD4+dpnSpb//U85cj5nkaxYwnjQXudPfPzOx7wC3AwEIr1tfX07t3bwDq\n6uro27cvAwYMAD7/YPO63NjYmKs8Wk5jed68AVx0EVxwQQM9egDkK5+W87fc0NDAiBEjAJZ+X5Yq\nZk9jZ2C4uw/Kls8Bmtz9ojbWXx54193rCjynnoZ0KzfcEK6FMWEC9OkTO42kKrWexnRgCzPrbWY9\ngCOAMS1XMLP1WiwOBmZXMZ9ILv3mN/DrX8Ojj6pgSPVFKxruvhgYCjxEKAaj3P15MzvfzJovOnma\nmT1rZo3ZuvVx0lZW8+5jilLODmnld4ef/xx+/3uYMgU23zyt/IUof3qizj3l7uOB8a0eG9bi/rnA\nudXOJZI37mHSwcmT4bHHYJ11YieS7kpzT4nk3JIl4VDaWbPgwQc1LYh0HV1PQ6TGLF4Mxx0Hr78O\nkyZBr16xE0l3p7mnciDlcdGUs0O+8y9aBEOGhGt6jxtXuGDkOX8xlD892tMQyaFPPoHDDoPllgtT\nm6+0UuxEIoF6GiI5s3AhHHIIrL463HEHrLhi7ERSq1I7T0NEWvn4YzjwQFh7bbjzThUMyR8VjRxI\neVw05eyQr/wffQSDBkHv3nDLLbBCEYPHecpfDuVPj4qGSA588AHsuy9su204eU/XwpC8Uk9DJLL5\n80PB+NrX4MorQ/NbpBoq0tMws63M7GEzey5b3s7Mzis3pIh87v33Ye+9YZdd4KqrVDAk/4r5J/o7\nwlQei7LlWcCRFUvUDaU8Lppydoib/733QsHYc89wASUr43Jj2v5xpZ6/HMUUjVXcfVrzQjYO9Fnl\nIonUvnffDQXjf/4nzFpbTsEQiaHDnoaZjQdOBe5x9+3N7FvACe6+XzUCFkM9DUlJc8H4+tfhwgtV\nMCSecnoaxRSNzYAbgF2B94G/A0e7+5wyc3Y5FQ1JRXPB2HffcE0MFQyJqSKNcHf/m7sPBNYCtnL3\n3fJUMGpByuOiKWeH6uavRMHQ9o8r9fzl6PD0ITMbRrietwFu2b90d/9FZaOJ1A7tYUitKGZ46kxC\n0QBYGTgQmO3ux1c4W9E0PCV51nyU1D77qIch+VKRnkaBN1kJmOjue5b0gxWkoiF51VwwBg6Eiy9W\nwZB8qdaEhasCG5bxc9KGlMdFU84Olc3//vvhCKm99qpcwdD2jyv1/OUopqcxq8XicsA6gPoZIu1o\nnhpk9911HobUlmJ6Gr1bLC4G5rl7rk7u0/CU5MmHH4Y9jJ12gssvV8GQ/OrSnoaZfbG9H3T390p5\no0pS0ZC8aJ7evG/fMJeUCobkWVf3NJ4CZrRzky6S8rhoytmha/MvWAD77w99+oTZaqtRMLT940o9\nfzna7Gm4e+9Kv7mZDQIuB5YHfu/uF7V6fiXgVmAH4F3gCHf/R6VziZTq44/hgANgq63gmms0W63U\nrqIOuTWzNYAtgJ7Nj7n7Y516Y7PlgReBvYE3gCeBI939+Rbr/AD4irv/wMyOAA5x9yEFXkvDUxLN\nv/8dLtG60UZw880qGJKOSl1P47vAY8BE4HzgIWB4OQFb6Q+84u5zssb6SODgVusMBm7J7o8GBnbB\n+4p0mU8+gW98A9ZbD266SQVDal8x/8R/SPiCn+PuewHbAx90wXtvCMxtsfw6y57/sXQdd18MfNBR\ngz5FKY+LppwdOpf/00/h0EOhri5c0zvGJVq78/bPg9Tzl6OYovGJu/8bwMx6uvsLwFZd8N5dOp5k\nfQ0bkN0GGVZvDG8YDoQPtuWHW395PVZv2PnZrT7u+hf+4cJc5dH6xa3f80Ljwf7G3SsbK54YP4/W\n1/odrd/Q0EB9fT319fUMHz6cchRznsb9wHGEPY6BhOnRV3D3/ct6x89fd2dguLsPypbPAZpaNsPN\nbEK2zl/NbAXgLXdfu8BrqachVbN4MQwZAosWwR/+AD16xE4kUp6Kzz1lZgOA1YAJ7r6og9U7eq0V\nCI3wgcCbwBMUboT3cfeTzWwI8A01wiWmJUvgmGPCGd/33w8rrRQ7kUj5KtUIv9LMdgVw9wZ3H9PZ\ngpG91mJgKKGxPhsY5e7Pm9n5ZnZQttqNwJpm9jLwI+AnnX3fPEp5XDTl7FBa/qYmOP54eOcduPfe\nfBSM7rT98yj1/OXocO4pwol855nZ1sC9wEh3n94Vb+7u44HxrR4b1uL+p8DhXfFeIp3R1AQnnQSv\nvQbjxsHKK8dOJBJH0cNTZrYm8E3gSGATd9+8ksFKoeEpqSR3GDoUZs6ECROgV6/YiUS6RjnDU8Xs\naTTbHNga+BJhOEmk5rnD6afDjBkwcaIKhkgxPY2Ls57CL4BngX7uflAHPyYlSHlcNOXs0H5+dzjn\nHHj0URg/HlZbrXq5ilXL2z8FqecvRzF7Gn8DdnH3dyodRiRPhg+HBx+EyZNhjTVipxHJh5Iv95pH\n6mlIV7vgArjjDmhogHXWiZ1GpDIq3dMQ6RZ+85swLcijj6pgiLSm6dVyIOVx0ZSzw7L5r7gCrr0W\nHnkE1l8/TqZS1Nr2T03q+ctR1J5GNo35ui3Xd/fXKhVKJIbrroNLLw17GBttFDuNSD4VM/fUqcAw\n4F/AkubH3b1PZaMVTz0N6aybboJhw0IPY7PNYqcRqY6KzD1lZn8D+rv7u50JV0kqGtIZt98OZ58d\njpLacsvYaUSqpyJzTwGvAR+WF0mKkfK4aMrZAYYNa+DHP4ZJk9IsGKlvf+VPTzE9jb8Dk81sHNA8\nUaG7+6WViyVSeffdB1deGYaktt02dhqRNBQzPDU8u9u8ohGKxvkVzFUSDU9JqcaOhRNPDHNJbb99\n7DQicVT0ehpm9gUAd/+ojGwVpaIhpZgwAb7zHXjgAdhxx9hpROKp1PU0+pjZ08BzwHNmNsPMvlJu\nSFlWyuOiqWX/05/g2GPDBZR23DG9/K0pf1yp5y9HMY3wG4DT3X0Td98EOCN7TCQpDQ1w1FEwejTs\nskvsNCJpKqanMdPdv9rRYzFpeEo6MnUqfPObcPfdMGBA7DQi+VCpuaf+bmY/A24jNMGPBl4tI59I\nFI8/HgrGnXeqYIh0VjHDU8cD6xAu9ToaWDt7TLpIyuOiec8+bRocfDDcdhvsvfeyz+c9f0eUP67U\n85ejwz0Nd38POLUKWUS61PTpcNBBcPPNsO++sdOI1IY2expm9lt3/6GZjS3wtLv74MpGK556GtLa\nU0/BfvvBDTeEPQ0RWVZX9zRuzf57SYHn9A0tufX006FgXHedCoZIV2uzp+HuM7K7fd29oeUN0Dm0\nXSjlcdG8ZZ85MxSMa66BQw7peP285S+V8seVev5yFNMI/06Bx+o786Zm9kUzm2RmL5nZRDOra2O9\nJWb2dHa7vzPvKbXvmWdC7+KKK+DQQ2OnEalN7fU0jgSOAnYHprR46gvAEncfWPabml0MvOPuF5vZ\n2cAa7v6TAut95O5fKOL11NPo5mbNgn32gd/+Fo44InYakTR06dxTZvYlYFPgQuBswjkaEKZJf8bd\nF3ci6AvAnu4+z8zWAxrcfesC66loSIeefTYUjEsvhSOPjJ1GJB1dOveUu/8j618cBTzRop/xPNDZ\ni2Gu6+7zsvvzCJeSLaSnmT1pZo+bWc22NFMeF42dvblgXHJJeQUjdv7OUv64Us9fjmLOCL8b2LXF\nchPwB+Br7f2QmU0C1ivw1E9bLri7m1lbuwmbuPtbZrYp8IiZzXL3gmej19fX07t3bwDq6uro27cv\nA7LTf5s/2LwuNzY25ipPKstrrTWAffaBE05oYIMNAPKVT8tazttyQ0MDI0aMAFj6fVmqYuaeanT3\nvq0e69TcU9nw1AB3/6eZrQ9MLjQ81epnbgYecPfRBZ7T8FQ303IP46ijYqcRSVOlLvf6Tsuhoez+\nO6WGa2UMnx+V9R1gmSOjzKzOzFbK7q8F7EaYnl26ueamtwqGSPUVUzS+D5xrZnPNbC7wE+CkTr7v\nhcA+ZvYS8D/ZMmbWz8x+l62zLfCkmTUCjwC/dvcXOvm+udS8+5iiamefOTMUjMsv75qCkfK2B+WP\nLfX85Shm7qlXgJ3MrFe2vKCzb5rNZ7XM9HHZCYXfze7/Bdius+8ltaOxEQYNCtf1Puyw2GlEuqdi\neho9gUOB3sDyfH6N8F9UPF2R1NOofTNmwAEHwNVX68Q9ka5Sqetp/BGYD8wAPiknmEhnPPFEmK1W\nkw+KxFdMT2NDdz/C3S9290uabxVP1o2kPC5a6ex//SsceCDceGNlCkbK2x6UP7bU85ejmKLxFzNT\nb0GqbsoUGDwYbrklFA4Ria+YnsbzwObA34FPs4fd3XNTSNTTqD0PPwxDhsDIkTCw7FnORKQ9lepp\n7FdmHpGyTJgAxx4Lo0fDHnvETiMiLRUzPNXUxk26SMrjol2d/f77Q8G4//7qFIyUtz0of2yp5y9H\nMXsaD/L5lfp6Ema+fRH4cqVCSfc0ciT86Ecwfjz06xc7jYgU0mFPY5kfMNsBOMXdT6hMpNKpp5G+\nm2+Gn/4UHnoI+vSJnUake6hUT+M/uPtTZrZTqT8n0parroKLL4bJk2GrrWKnEZH2dNjTMLMzWtzO\nMrO7gDeqkK3bSHlctLPZf/3rMI/UY4/FKRgpb3tQ/thSz1+OYvY0erW4vxh4AFhmenKRUriH4ag/\n/jEUjHA9DBHJu/Yu93qbu3/bzH7k7pdXOVdJ1NNIS1MTDB0apgeZMAHWWit2IpHuqat7Gv3MbAPg\neDO7tfWT2Uy1IiX57DM47jiYOxceeQRWWy12IhEpRXs9jeuAh4GtCJMVtrxNr3y07iPlcdFSsv/7\n32GG2vnzwx5GHgpGytselD+21POXo82i4e5XuPs2wM3uvmmr239VMaPUgA8+gH33DYXivvtg5ZVj\nJxKRcpR8nkYeqaeRb/PmhYKxxx7hSKnlipmHQEQqrlLXCBcp29/+BrvtBt/8Jvz2tyoYIqnT/8I5\nkPK4aHvZn3467F2ceSb8/OdgJf09Ux0pb3tQ/thSz1+Oks8IFynGI4+Eqc2vvVaXZxWpJeppSJcb\nORJOOw1GjYK99oqdRkTaUpW5p0Tac9llcOml8Kc/wXa5uUyXiHQV9TRyIOVx0ebsTU1wxhnwu9/B\nn/+cTsFIeduD8seWev5yRCkaZnaYmT1nZkuyqdbbWm+Qmb1gZi+b2dnVzCjF++QTOOIIePJJmDoV\nNtkkdiIRqZQoPQ0z25pw9b/rgTPc/akC6yxPuNjT3oRZdZ8EjnT35wusq55GJO++CwcfDBtvHK6J\n0bNn7EQiUqxkztNw9xfc/aUOVusPvOLuc9z9M2AkcHDl00mxXnkFdtklnIdxxx0qGCLdQZ57GhsC\nc1ssv549VnNSHBedOhX++7/hwAMbuOiidE/aS3Hbt6T8caWevxwVO3rKzCYB6xV46lx3H1vES5Q0\n3lRfX0/v3r0BqKuro2/fvgwYMAD4/IPN63JjY2Ou8nS0fN55DVx1FYwaNYCVVoqfR8ta1nJxyw0N\nDYwYMQJg6fdlqaKep2Fmk2m7p7EzMNzdB2XL5wBN7n5RgXXV06iCpiYYPhxuuw3GjNG1vEVSl+p5\nGm0Fng5sYWa9gTeBI4Ajq5RJWlm4EOrr4fXXYdo0WGed2IlEJIZYh9weYmZzgZ2BcWY2Pnt8AzMb\nB+Dui4GhwEPAbGBUoSOnakHz7mNevfEG7LknrLhimB6kZcHIe/aOKH9cyp+eKHsa7n4fcF+Bx98E\nDmixPB4YX8Vo0sq0aWHuqKFD4eyz8znpoIhUj+aekjbddhucfjrcdBMcdFDsNCLS1VLtaUjOfPYZ\nnHUWjBsHkyfDV74SO5GI5EWiR9fXljyNi779drjK3osvwhNPdFww8pS9HMofl/KnR0VDlnrySdhx\nR9hpJ3jgAVhjjdiJRCRv1NMQIMxO+9OfwvXXwyGHxE4jItWgnoaUbOHCcGTUtGlhapAtt4ydSETy\nTMNTORBrXPSll2DnneHTT0PRKKdgpD6mq/xxKX96VDS6qVGjwoSDQ4fC7bdDr16xE4lICtTT6GYW\nLoQf/SgcSjtqFOzQ5iWwRKTWJXM9DYlj9mzo3x8+/hhmzFDBEJHSqWjkQKXHRd3h2mvD/FGnnx6G\no1ZbrWteO/UxXeWPS/nTo6Onatzbb8MJJ4RJB6dOha22ip1IRFKmnkYNGzcOvvtdOOYY+OUvoUeP\n2IlEJE90noYAsGBBGIaaOBHuuisMS4mIdAX1NHKgK8dFGxrgq18Nkw4+80zlC0bqY7rKH5fyp0d7\nGjViwQI45xy47z647jo48MDYiUSkFqmnUQMmTYKTToLdd4fLL9dEgyJSHPU0upn33gu9i4aGcEjt\nfvvFTiQitU49jRwodVzUHW69FbbdNpxvMWtWvIKR+piu8sel/OnRnkZinn8eTj4ZPvoIxo4N178Q\nEakW9TQS8eGH8ItfwIgR8POfwymnwPLLx04lIinT3FM1qKkpDEVtsw28+y489xycdpoKhojEoaKR\nA22Ni06ZEi69evXVMHo03HwzrLtudbN1JPUxXeWPS/nTE6VomNlhZvacmS0xszbnWjWzOWb2jJk9\nbWZPVDNjTC++CIceCkcfDf/7v/D44+FiSSIisUXpaZjZ1kATcD1whrs/1cZ6fwf6uft7HbxeTfQ0\n3ngDzj8/nKB35plw6qmwyiqxU4lIrUqmp+HuL7j7S0WuXtIvlKJ588IeRZ8+UFcX9jTOPlsFQ0Ty\nJ+89DQcmmtl0M/tu7DBd7Z//hLPOgs03b2DJktDkvvhi+OIXYycrXupjusofl/Knp2LnaZjZJGC9\nAk+d6+5ji3yZ3dz9LTNbG5hkZi+4+5RCK9bX19O7d28A6urq6Nu3LwMGDAA+/2DzsjxyZAOjRsGj\njw7g6KPh9NMb2WsvWH/9fOTTspa1XJvLDQ0NjBgxAmDp92Wpop6nYWaTaaen0WrdYcACd7+kwHNJ\n9DRmzIBLLoGHHoITTwxTgOTtaCgR6T6S6Wm0UjCwma1iZl/I7q8KfB2YVc1gXWHx4nC47IAB8I1v\nQL9+8OqrcNFFKhgikp5Yh9weYmZzgZ2BcWY2Pnt8AzMbl622HjDFzBqBacAD7j4xRt5yvPFGuFre\nppuGmWe///1QLM44A1Zf/T/Xbd59TFHK2UH5Y1P+9ESZe8rd7wPuK/D4m8AB2f1Xgb5VjtYpixbB\n+PFw443hetyHHQZjxsD228dOJiLSNTT3VCe5w7RpcPvtMGoUbL01HHccHH449OoVJZKISFF0PY0q\nWbIEnngC/vCHcFtlFTjqqPDYppvGTiciUjl5aIQnYf58uPdeOP542GAD+N73wp7EuHEwezb87Gfl\nF4yUx0VTzg7KH5vyp0d7Gm346KMw59Njj8Gf/hROvNttNzjggM4VCBGRlKmnQWhgv/giTJ8ehpie\neCIs77BDuO72wIGw667Qs2cXhhYRiaycnka3KRpLlsBbb8HcufDKK/Dyy+H27LNh+UtfCudQ9O8f\nroa3ww4qEiJS27p10Rg2zFm8GD7+GBYsCMNL774L77wDb78N//oXrLUWbLwxbLYZbLEFbL45fPnL\n4QJHK68cL39DQ8PSU/5Tk3J2UP7YlD+ubn/0VM+esOaaoUHdq1eY+G/ttUOxWG896NEjdkIRkbTV\nzJ5GLfweIiLVlOrcUyIikggVjRxI+VjvlLOD8sem/OlR0RARkaKppyEi0k2ppyEiIhWlopEDKY+L\nppwdlD825U+PioaIiBRNPQ0RkW5KPQ0REakoFY0cSHlcNOXsoPyxKX96VDRERKRo6mmIiHRT6mmI\niEhFRSmIdVJTAAAIFUlEQVQaZvZ/Zva8mc00s3vNbPU21htkZi+Y2ctmdna1c1ZLyuOiKWcH5Y9N\n+dMTa09jIvBld/8q8BJwTusVzGx54CpgELAtcKSZbVPVlFXS2NgYO0LZUs4Oyh+b8qcnStFw90nu\n3pQtTgM2KrBaf+AVd5/j7p8BI4GDq5WxmubPnx87QtlSzg7KH5vypycPPY3jgQcLPL4hMLfF8uvZ\nYyIiEknFLvdqZpOA9Qo8da67j83W+SmwyN3vLLBetzkcas6cObEjlC3l7KD8sSl/eqIdcmtm9cB3\ngYHu/kmB53cGhrv7oGz5HKDJ3S8qsG63KTAiIl2p1ENuK7an0R4zGwScBexZqGBkpgNbmFlv4E3g\nCODIQiuW+kuLiEh5YvU0rgR6AZPM7GkzuwbAzDYws3EA7r4YGAo8BMwGRrn785HyiogINXJGuIiI\nVEcejp4qSkcn+pnZSmY2Knv+r2b2pRg521JE/nozezvb83razI6PkbMQM7vJzOaZ2ax21rki+91m\nmtn21czXkY7ym9kAM/ugxbY/r9oZ22NmG5vZZDN7zsyeNbPT2lgvl59BMfnz+hmYWU8zm2ZmjVn2\n4QXWye13T5H5S/vucffc34DlgVeA3sCKQCOwTat1fgBck90/AhgZO3eJ+b8DXBE7axv5dwe2B2a1\n8fz+wIPZ/Z2Av8bOXGL+AcCY2Dnbyb8e0De73wt4scC/n9x+BkXmz+1nAKyS/XcF4K/ATq2ez+13\nT5H5S/ruSWVPo5gT/QYDt2T3RwMDq5ivI8Xkt+yWO+4+BXi/nVWWbnt3nwbUmdm61chWjCLyQ063\nPYC7/9PdG7P7C4DngQ1arZbbz6DI/JDTz8DdF2Z3exD+6GtqtUqev3uKyV/Sd08qRaOYE/2WruOh\nif6BmX2xOvE6VEx+Bw7NhhbuMbNCZ8nnVaHfL6X8DuyS7cI/aGbbxg7Uluxowu0JMym0lMRn0E7+\n3H4GZracmTUC84CJ7v5kq1Xy/N1TTP6SvntSKRqpd+uLyT8W+JKH+bgm8flfLqlo/ZdKSp/ZU8Am\n7t6XcGTf/ZHzFGRmvYA/AD/M/mJfZpVWy7n6DDrIn9vPwN2bslwbATuZ2ZdjZypFEflL+u5JpWi8\nAWzcYnljwl9SrdfZBMDMVgBWd/f3qhOvQx3md/f3sqErgBuBflXK1hVa/34bZY8lwd0/at6Fd/fx\nwIp5+ksRwMxWJAx93O7uhb5Qc/0ZdJQ/hc/A3T8AJhMmUW0pz989S7WVv9TvnlSKxtIT/cysB6HZ\nNKbVOmMIDR2AbwEPVzFfRzrMb2Ytp1wZTDg3JRVjgGNh6Zn88919XtxIxTOzdc3Msvv9CYei5+Z/\n+izbjcBsd7+8jdVy+xkUkz+vn4GZrWVmddn9lYF9CD2ZlnL73VNM/lK/e6KcEV4qd19sZs0n+i0P\n3Ojuz5vZ+cB0D3NZ3QjcZmYvA+8CQ+Il/k9F5j/NzAYDiwn566MFbsXM7gL2BNYys7nAMEJDDXe/\n3t0fNLP9zewV4GPguHhpl9VRfsL/6Ceb2WJgITn6t5PZDTgGeMbMns4eO5fsr9sEPoMO85Pfz2B9\n4BYLl2pYjnCS8YOpfPdQXP6Svnt0cp+IiBQtleEpERHJARUNEREpmoqGiIgUTUVDRESKpqIhIiJF\nU9EQEZGiqWhIt2BmcwqdYWxmB1mBqeq7+L2Hm9kZJa7/eqFprEt83xFmdmh2/w4ze7d5WaRcSZzc\nJ9IFnAIzeWYnN42twnuXuv6l7n5p6yfMbIVsUryS3tfdjzazm8vIIvIftKchNSObpuUFM7vdzGZn\nM3au3GKVU81shpk9Y2ZbZT9Tb2ZXFnit/mb2FzN7ysz+bGZbtlj/XjMbb2YvmdlFLX7mBDN70cJF\nb37Xxutulv3sdDN7rDlHoV+nxc8MN7PbzGwq4ezeL2U/OyO77ZKtZ2Z2VbYNJgHrtPe6IuVQ0ZBa\nsyVwtbtvC3xIuEBOs7fdvR9wLXBm9lhbf3k/D+zu7jsQph35VYvnvgocDvQBjjCzDc1sA+A8wgWQ\ndgO2avXazfdvAE51968BZwHXFPl7bQ0MdPejgX8B+2S/yxDgimydQ7LffxvCPFS7tvP7iZRFw1NS\na+a6++PZ/duB04BLsuV7s/8+BXwzu9/WX951wK1mtjnhi7fl/ysPu/tHAGY2m3BFxrWBR919fvb4\nPYQv8KXMbFXCF/k92dx8EC6M0xEnXNXu0xY/c5WZfRVYAmyRPb4HcKeHuYHeMrNH2vn9RMqioiG1\npuVf1tZquflLdwkd/9v/f4TicIiFaz43FHidlq/V+i/6Ql/WywHvu3s51+9e2OL+/wJvufu3s4no\nPskeL9S30Z6GdCkNT0mt2SSbGhzgKGBKma+zGvBmdr+jGWMdeBLY08zqLFxT4VA+/8I2wuSgHwF/\nN7NvwdIexHZlZvtndv9YwszJAI8RhsuWM7P1gb3KeG2RdqloSK15ETglGzZandC/gGX7C17gfksX\nA782s6cIX8rtru/ubxL6Hk8AU4G/Ax8U+JmjgRMsXH7zWcL1C4rR8j2vAb6TvcZWwIIsw33Ay4Tr\nIdwC/KXI1xYpmqZGl5ph4frTY929T6T3X9XdP872NO4lXDflj2W8zjBggbtf0uHKpb3uCML2Gd2V\nryvdi/Y0pNbE/CtoeHaRoVnAq+UUjMwC4HudPbmvJTO7A9gd+HdXvaZ0T9rTEBGRomlPQ0REiqai\nISIiRVPREBGRoqloiIhI0VQ0RESkaCoaIiJStP8P3QGH9wRL3JgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x106c617d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot this function for a 1 GeV particle\n",
    "\n",
    "pylab.plot(\n",
    "    phiValues,\n",
    "   [ func(phiValue, 3.8, 1) for phiValue in phiValues]\n",
    ")\n",
    "# draw a horizontal line at zero\n",
    "pylab.plot(pylab.xlim(), [0,0], '--')\n",
    "\n",
    "pylab.grid()\n",
    "pylab.xlabel('phi angle [rad]')\n",
    "pylab.ylabel('function value')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.66759049589\n"
     ]
    }
   ],
   "source": [
    "# find the phi angle where the particle hits the ECAL in the above example numerically\n",
    "import scipy.optimize\n",
    "\n",
    "print scipy.optimize.brentq(func, \n",
    "                            0, math.pi, # bracketing range\n",
    "                            args = (3.8, 1))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "ename": "ValueError",
     "evalue": "\n\\frac{\\left\\lvert{p}\\right\\rvert}{\\sqrt{mP^{2} + p^{2}}}\n           ^\nExpected a delimiter (at char 11), (line:1, col:12)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/IPython/core/formatters.pyc\u001b[0m in \u001b[0;36m__call__\u001b[0;34m(self, obj)\u001b[0m\n\u001b[1;32m    335\u001b[0m                 \u001b[0;32mpass\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    336\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 337\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mprinter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    338\u001b[0m             \u001b[0;31m# Finally look for special method names\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    339\u001b[0m             \u001b[0mmethod\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_safe_get_formatter_method\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mprint_method\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/interactive/printing.pyc\u001b[0m in \u001b[0;36m_print_latex_png\u001b[0;34m(o)\u001b[0m\n\u001b[1;32m    125\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mlatex_mode\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0;34m'inline'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    126\u001b[0m                     \u001b[0ms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mlatex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m'inline'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 127\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0m_matplotlib_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    128\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    129\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_print_latex_matplotlib\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/sympy/interactive/printing.pyc\u001b[0m in \u001b[0;36m_matplotlib_wrapper\u001b[0;34m(o)\u001b[0m\n\u001b[1;32m     88\u001b[0m         \u001b[0mo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreplace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'\\operatorname'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m''\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m         \u001b[0mo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mo\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mreplace\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34mr'\\overline'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34mr'\\bar'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 90\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mlatex_to_png\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     91\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     92\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_can_print_latex\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mo\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/IPython/lib/latextools.pyc\u001b[0m in \u001b[0;36mlatex_to_png\u001b[0;34m(s, encode, backend, wrap)\u001b[0m\n\u001b[1;32m     85\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     86\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'No such backend {0}'\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mformat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbackend\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 87\u001b[0;31m     \u001b[0mbin_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwrap\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     88\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mencode\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mbin_data\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     89\u001b[0m         \u001b[0mbin_data\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mencodestring\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbin_data\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/IPython/lib/latextools.pyc\u001b[0m in \u001b[0;36mlatex_to_png_mpl\u001b[0;34m(s, wrap)\u001b[0m\n\u001b[1;32m    104\u001b[0m     \u001b[0mmt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmathtext\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMathTextParser\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'bitmap'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    105\u001b[0m     \u001b[0mf\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBytesIO\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 106\u001b[0;31m     \u001b[0mmt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_png\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m12\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    107\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgetvalue\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    108\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/mathtext.pyc\u001b[0m in \u001b[0;36mto_png\u001b[0;34m(self, filename, texstr, color, dpi, fontsize)\u001b[0m\n\u001b[1;32m   3092\u001b[0m         \u001b[0mimage\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpixels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3093\u001b[0m         \"\"\"\n\u001b[0;32m-> 3094\u001b[0;31m         \u001b[0mrgba\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdepth\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_rgba\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtexstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcolor\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfontsize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3095\u001b[0m         \u001b[0mnumrows\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumcols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtmp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mrgba\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3096\u001b[0m         \u001b[0m_png\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwrite_png\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mrgba\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtostring\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumcols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnumrows\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfilename\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/mathtext.pyc\u001b[0m in \u001b[0;36mto_rgba\u001b[0;34m(self, texstr, color, dpi, fontsize)\u001b[0m\n\u001b[1;32m   3057\u001b[0m             \u001b[0mimage\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mpixels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3058\u001b[0m         \"\"\"\n\u001b[0;32m-> 3059\u001b[0;31m         \u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdepth\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_mask\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtexstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfontsize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3060\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3061\u001b[0m         \u001b[0mr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmcolors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolorConverter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mto_rgb\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcolor\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/mathtext.pyc\u001b[0m in \u001b[0;36mto_mask\u001b[0;34m(self, texstr, dpi, fontsize)\u001b[0m\n\u001b[1;32m   3030\u001b[0m         \u001b[0;32massert\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_output\u001b[0m\u001b[0;34m==\u001b[0m\u001b[0;34m\"bitmap\"\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3031\u001b[0m         \u001b[0mprop\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mFontProperties\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msize\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mfontsize\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3032\u001b[0;31m         \u001b[0mftimage\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdepth\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtexstr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdpi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprop\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mprop\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3033\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3034\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mftimage\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_array\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/mathtext.pyc\u001b[0m in \u001b[0;36mparse\u001b[0;34m(self, s, dpi, prop)\u001b[0m\n\u001b[1;32m   3003\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m__class__\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_parser\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mParser\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3004\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3005\u001b[0;31m         \u001b[0mbox\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_parser\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mparse\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfont_output\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfontsize\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdpi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   3006\u001b[0m         \u001b[0mfont_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mset_canvas_size\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbox\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mwidth\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbox\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheight\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbox\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdepth\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3007\u001b[0m         \u001b[0mresult\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfont_output\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget_results\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbox\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/opt/local/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages/matplotlib/mathtext.pyc\u001b[0m in \u001b[0;36mparse\u001b[0;34m(self, s, fonts_object, fontsize, dpi)\u001b[0m\n\u001b[1;32m   2337\u001b[0m                         \u001b[0merr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mline\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2338\u001b[0m                         \u001b[0;34m\" \"\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0merr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcolumn\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m\"^\"\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2339\u001b[0;31m                         six.text_type(err)]))\n\u001b[0m\u001b[1;32m   2340\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_state_stack\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2341\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_em_width_cache\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m{\u001b[0m\u001b[0;34m}\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mValueError\u001b[0m: \n\\frac{\\left\\lvert{p}\\right\\rvert}{\\sqrt{mP^{2} + p^{2}}}\n           ^\nExpected a delimiter (at char 11), (line:1, col:12)"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$$\\frac{\\left\\lvert{p}\\right\\rvert}{\\sqrt{mP^{2} + p^{2}}}$$"
      ],
      "text/plain": [
       "     ~p~     \n",
       "~~~~~~~~~~~~~\n",
       "   __________\n",
       "  ~   2    2 \n",
       "~~  mP  + p  "
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# mass of the particle (in GeV/c^2)\n",
    "mP = Symbol('mP', real = True)\n",
    "\n",
    "beta = Symbol('beta')\n",
    "gamma = 1 / sympy.sqrt(1 - beta * beta)\n",
    "\n",
    "# velocity of the particle\n",
    "betaSol = solve(sympy.Eq(p, mP * gamma * beta), beta)[1]\n",
    "betaSol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAEoAAAAsCAYAAADCfS42AAAABHNCSVQICAgIfAhkiAAAA11JREFU\naIHt2l+IVVUUx/HPqJOZ/aFsLCqpSIqkTErFMNOaQeyhqEjCZDAqKIoiIiI0YYKCCnsLeouejAl6\nKSohDXrSEqLpP5UUREF/QNIX+2P1sPZl7lzvvXPPved2zuT5woW915yz1++s2f/OPouKjphdtIDE\nfDyIx3AEy3EnjuLb4mSVj1Gcgr24NdmW4OPCFJWU03AifsGcZFuPHwpT1MCsogUkfsMq7MNfybYe\n7xSmqIE501/yn7EWn6TyEG4SwSoFA0ULqGMPPsCnWImXMFGoohIyiJ+VZxU+hrLMUStETzpatJAy\nsxTvivnphoK1VFRUlIva9uCfQlVUVFRUNCfrK8zJONwPIf83bsbFRYsogqyvMJfgq34IKTtZA3Xc\nbiOyBGoZPuqXkLKT5eBuHV7sk44as/AA5qX6czPR3+N5NDINN2JRKr+Gq8rir9OhtwC/9iiqEy7C\nplQ+YPIhZoy/zXk00gFzxWcreBvnlMVfp3PUInzfo6hO+D39VuM9/FgWf50EarZ8j2iHcD+ewKvi\nZHMQF4rveDtwHZ7K0Wc7Tu3G33aTK0CN1ViTk6ga54vgL6izDeAg7hWBG8RIhjY3dqkls7+1+LDJ\nxY/K/+vI3Xi/wXam2ND+IRaOg7gsQ5tjXei4HYc68Vc/9CawE8PYXWcfdOzQOwNbRHCfxaWiC1+A\nN3C2eN05JIZSIyPYVVefi2eS/1H83frZMjNf9Jo1eFm82F+JtzCefplZhv119XOF8EbuEQH8Gnck\n20liYrw61Zfg8yb3DuAnPI/78BC26X14j7Ww55IA0jiZT4gV7nTRFTeIHtLIOBaKnrAz2a4QQd6b\n6svxWZN7LxdJGdtEik83nIWHTT0mukYketQ4jKfxOv7E4lSG80ydH7tiHLek8tY2120RXbnGdjxZ\nV9+F20TQ63lELMV5M9bmb+tM/YfvMFX7tDTbme8Wc8gJYii1Ylh8uKwxIvIHiOCsSuLuanJfPwLV\njmYJIGNZGmgVqOHUeLsHWmwyLWcgCdiX6kdSeRRvJttKMfFfL4bfhixCe+TapHGz2L9txHd5NHwA\nLyhXtst0bGphzyUBpNVL8R6xis2kg7pXWthzSQBpFeV5+BJf9NJ4CVgqzpiGxFD7plA1FRUVFRXF\n8i83CZaQNNxbFwAAAABJRU5ErkJggg==\n",
      "text/latex": [
       "$$\\frac{p}{\\sqrt{mP^{2} + p^{2}}}$$"
      ],
      "text/plain": [
       "      p      \n",
       "~~~~~~~~~~~~~\n",
       "   __________\n",
       "  ~   2    2 \n",
       "~~  mP  + p  "
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# calculate beta (the velocity) from the particle's mass an momentum\n",
    "#\n",
    "# E = m * gamma \n",
    "# p = m * gamma * beta\n",
    "#\n",
    "# beta = p / E\n",
    "# where E = m * m + p*p\n",
    "beta = p / sympy.sqrt(mP * mP + p * p)\n",
    "\n",
    "beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKoAAAA1CAYAAADRRnAUAAAABHNCSVQICAgIfAhkiAAACBJJREFU\neJzt3XmsXHUVwPHPo68LdtG2NLYWpGKxgWhttCy1paVQQapCU3wWWgtiqyAaI9AgYqHPuGs1cUmM\n0cS1uMQlQlwSQUXUVsAloKIiakEwghUriBRon3+cO859w5s782Z5703n900mc+/9bWfuPff3O79z\n7+8MiUQH0DPaAiQ6hnF4FR7FVHyqznKH4A04NNt/f+tFSyTKnIpF2falmFlnuZfjiGz7a3hhI40f\n0kihRFeyH1diUvZ5sM5yz8a52fZdykqbSLSNL+BPOGMYZSYKUwG+g2c00vC4RgqNAsdn3w9jNp6L\nv3ZJnWOFuWIE/hWuwmfxWB3l9mf5lmJAKGvTTMVX8cwm81ZLOwGXoB/fw/KK9POwFZtxWe74Z8SP\n3I+dODqXthgfwfn4JObXmVZUZ5EsjdZ5khgCN+HzWFnRXiPns4hqZarJPxWvHuJzepZ+NXqz7fOF\n7Vkv08S5rIeNmFyUYbNQoAHMq1FZUd5qaVPwntz+K/GIuFOJk/LebHueuAsPy/b7RQ81u6KtCdiN\nOdn+cbiljrSiOotkaabOf4gLDH34jzgnNHY+i6hWppb8RWxTnrkvxrI6y8GFGJ99VtXIe1G9lQ7n\nhBTlrUxbiAPCuCbusgGhsBPwAI7M5T8qt91fpY0X4ze5/R4x7D6rRlpRnUWyNFonYQaUeoo+PK5s\nu5UYzvmsh8oyteQvYgYuFyPNumHIsA7/Fjfqg+I8VGMcLhgqoXeog23idiwRxjjl2d+d2fHDhHIs\nEy6M7+byThZD5n/Fyd4uTvg87Mm1MaB8MmYXpP25oM4iWYraK6oTfp0rtxbvwENVzlW7mKdY/iL+\nqboPdLLoNU8S5s8UvADfxpezTxFLhCk0HfeI0eyjuLeUYSQVdQA/y+1fgQ/hl1ify7MD3xSTkNJk\n5OvYlaU/kO0fIxTqkYp2Sg7pojQFdZZMkaFkabTOA1n68ThNKPJ2reOlQuF+WiNfLfkbZa2wd/tE\nr7hDXNcviZG0iAuEf/bNeB0+gWPxOeG7xej5UV+Dv2FLtr83+741+35YXMyXZfu3iItP9MDPwfOy\ncpVP16aIYaYorVad1WRptM4SN+Od2fdNyjZqs6zHWXXkqyV/o1wrTJn52TYcrvZDgcPxdmFSDORk\n2y0m3v9nNBR1dfZ9ufCxzRMujwGD3WUDosc/Ef8STmbKd//j+J3BE5deMXzsrpFWVGeRLI3WeSL+\nrmwL/lAMjS/RGjbgLXXkK5K/GfaK37gLT2THThOenSJeJG7ufeKGKZlCZ+BHtRodymg/Hc+vM29R\n2nJhw5VmxmuEfQLfV3aFzBI/fq6467bl6rgIP8m2e3GfsgtmJX5RR1pRnUWyNFrn4qzOidn+auFJ\nqHSJtXsyVSR/s2zDu7PtWfiD2vIuUDZXVot5yxz82ODJ7KBeYwMuFid1jlCkm7O07ULjb6ojb7W0\no8Rd0ieG/C04Rzw33icu5MXCK7BO+N1uEzPGaThTGOsL8Vrh3jkgesDXi2HmLGHr7KmRVlSnAlka\nrfM+4Vs9RfQ8a4RvtuQaauR8FlGtTJH8zXK1+J2zhb/4Svy+Rpk9wrTaKDqsXuEyu0TY+IlESxmP\n+zX3pPPCosT0UkqiFRwn3G/7Gyw/SUxYE4m2sVCYSrcb3ssqeU5Q/CAgkUgkEolEd1F6EjBQmCuR\nSCQSiUQikRg5un1df/5FiG6k269/x7BGvIqXGON0+yPUBeItn8QYp9sVNbnlOoRuVtRF4pW3RAcw\nkmumxhon4+NtbmM0AoSloGQHGVeMQBstCRDWAW22nW4d+mdqfkFbPYxGgLAUlOwgYoORuYAtCRDW\nAW22nW61UY8QgQ7azb7ssxQ3ijVFB2ObbadTFXWSWBx4qVhwuDM7frRy1I69Qxc1TuNLJoZillhI\ntxVfEW+6jxdLo+8VCyNXivX8I8W0UWgzUYXSOpsFFcd/q3yBrlKe/ZZYKlaJtpIjhfLnAy70iOgl\nwwkQlqevCXkabTPRBlZ58rA2TgzpW7FCrFmvvFBbtD4u7CaDwxUR4XMGxPr9egKEVdLfoCzDCUrW\nMXTq0E8o4A9y+70ibtHP8QHRk14j4hddn8s33tBD/wwRFnIF3idiRk0TQRSuE+vVFwglqIwbtUoE\nUisxUQT6ukasWT+gtTQblCwxgtwqFPMcMYvfgQ8b3FsuMjj+51yhOEOxWSjxncpB254iJialaC7H\nCtMiT48I1/NBER3lTXib5s2L/oK0jWJmv1MEKCvJdluTbY5ZOvU1rxnKMUzzocevxd14Y7bfI4Kx\nHSOGwU2id7x/iDqnih50p3LImyWidy4FrT1PONTz9uNCEYXkaSIyXiM8XUQsyV+PZSK0TYmH8K5s\n+6niBrpHREJ5QsR6+rRyNMKDik4d+k/BHZ4cH/9uMVkqMSBcNCfjG0IhhlJSQhHWijXqJVbhhtz+\nevH/StOV/xVklbBPG1VSokd+a8WxftV71b3iNw03KFnH0qlPpiptQmLoWy96wDzXZ/kniF6oiFNV\nV9TpIm7UdSJsZr7MjfUK3kJWCFcY4SI7U+MTsDFPp/WoJWV8hRgW+7PjM4UT/2yDJ1iEol4mLmwt\nhZqvHL6xRyjArmz/0Wx7I74lgvKeLXr3x0QIycqbp50sFybHhkyWPvxlBNtPtIG78DGdZZOfW5DW\niqBkHUWnDv3D5QYxg++kF6W/WJDWbFCyjqNb7shDRbTlO0ZbkBawULxjOksM9X8cVWkSiUQikUgk\nEolEIpFIJJrif7EJHzkPABDHAAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$\\frac{1.11265005605311 \\cdot 10^{-8} p^{2} \\phi}{B \\sqrt{mP^{2} + p^{2}}}$$"
      ],
      "text/plain": [
       "                     2  \n",
       "1.11265005605311e-8~p ~~\n",
       "~~~~~~~~~~~~~~~~~~~~~~~~\n",
       "         __________     \n",
       "        ~   2    2      \n",
       "    B~~~  mP  + p       "
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# time of flight of the particle \n",
    "\n",
    "# speed of light\n",
    "from scipy.constants import c\n",
    "\n",
    "tof = beta * phi * R / c \n",
    "\n",
    "tof"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAKwAAAASCAYAAADG8TXdAAAABHNCSVQICAgIfAhkiAAABoZJREFU\naIHt2nusXFUVBvDfhbaX2vbW1mJbqKEPa33TaypqqxUfiRqNkkY0JhWbmAA1UUEpKopW0qIVqaih\nChpFS7RRIfVBYjU+qBGwPJSAgihioUZEBRGFEtpe/1j7ZPaZe+beM6czGM39ksnMWfvbe+2zz358\na51hAhP4H8LAOOVnYTI+VqOtF+I0PIKp6bMZtzTkwdOwMXEfTd8b8VAf/dblnYD3p/IF2IMP408N\n2nsNPo+bs3s9lJXvwbb0++P4Nm5PvEV4PbZjX1ZnAT6AEQziCalu1Tj3G8vEHNqb+jMHG/CXNt5i\nfAT7cRDT8D7cW8fJcfi3mCDjYRjfwVGZbRv+ieUNePBM3I2V6Xoe7mrrT6/91uU9Dz/AE9P1dOzG\nfVjYoL33igdZ9TmEV2Xc9vKDOFcZT8JOzM9si/BbLPX4YqZYSGsz2zm4FVMy2yL8DWsy21qxiCfV\ncXSpGJCNNbifStyTM9vrku0zDXiTxOCekdmegr/i3X30W5d3FZ6qjOHE29Ggvc+l+5uMIzL7Kny6\nzc9e8Wx24kIcbzQ24J0V9nOwpcLeT2wWCzmfdLPxGNZntp24X/n+p4rddl1hyAtzrMGPu+jUL8Wu\n8WBmm5G+H27AO0WsuMsy2z04WvkB9tpvXd5Lxfg8ua3uP/DKBu0dFPf3mJYUmI4PCtmR4y6cipPE\nznyz0ViMl1fY94tF8XjiZPwCBzLb/ULSFAt5Cl6L3ytLoUfEAn3jWA6m44vpd90dtgoXpE5W7QDj\n8X6I3/wX/Nbl3aKlH3PcK2RUL/xuw4oK+0/HqUdo5hF8DbOSbRA34bk16vcKM1I/Lq4o26W1kOcl\n3u4K3h6xQ3fEZixJv5tO2CXi+F7fgDcgVtbVWI1NuEgcGcN99NsNb5oY5BzHiPH6SQ/8rhJBWBWu\nEcHUhWLyXymC0xyDuD71589CC14ugrsqLBcyYwu24qtaJ8Hh4NmpD1VB+5VaAeGRYqHvqeDtS7xK\nHbtcDEaBbifsSfgs7kjtdMpCjMWbk/zehtMz+4kiO/CsPvltwstxvjjaV1aUddvejcrBW447hd4t\nsFZE2+0LaAa+rxWcfQ9zK9pbJ2TFsZntbLxhnD7Wwcrk+6MVZdtTWdGnL4mgKx+beUIijCjLL4Se\n/Ypy5NZ0h50idpqfi4i1G97c5He/EN059omB74ffw+EtFotp0xicuu29wthy6MiK64eVgzhC+35B\naMM/aO22z8k4LxK6eVVmO0EElaMmSAO8QOc5tCOVHZOujxYL+rR0PUmM502JN2q81uNlbbbD0bAn\npvrf7JI3OV3fWsG9TmjHwT74bcobTP3aOk47ddu7Quw23WBv+hQ4U0y6AtNExuKQcpC2S+jDC/DJ\nxDkVQ23tHy8CyF/V/Hw51Vus8xz6birLpccskYfdKibrcbhBSMTSqTRPHFntqDthn260mB/SyiNO\n75J3n2oBfnXiFvnFXvuty8uxHedV2Ju0Nxn/EnFEFXbjZxX2feJEIh7sA0I/tuP05PsZydcBXNLB\nVy8wTdxne2qOyLI8UKONe/CjduNaEZnvzD5XiZu7PV2vaa+UMCRWwAGtYK3obKGfZnbBI5Ltt1X4\nuk48mEl98NtN/wpsNDppf0r6btLeqmQ/c9SdBx4SO1g79ou8NXGUF8FMOwZEZD6sFZmf3cFXr3Cj\nCLDacbexA1RCJowoxzIdsVDnHXaplr48Suig34mEcIEVqf71XfLgLUKX5Rp2QOQ5v9Env930j1jk\nVTvrpQ3bg7elsk5ZhG+JoDRH8cKiWDgDIr324or6M8RuPCgW/YPiZUI75hsn99kFzhOvq/MjfYno\n8zsy2xniZF2Q2c5K/c3fFHbE0tRo+/G0WmzzuzLb+clh3qnt4nhb0YB3BK5NHS7wJpESWthHv3V5\nq1NfLm/77MDXG7RXYIMY87dXlMHzU/1i9xwQKahrlHfUN+PXymM1JHRzPhE/IWRG3r/Cxyy9wXyx\n0bw1s10kctl5gP8h/FEr2zEssh+lrEtVemVI/LliWXL2qBDcW4Q0WCZu8grlFbIOrxZH4FzxNuNc\nEflpwJstxPc0sVMNpJu6s2F7veT9XXnXzLFJWSbU9Uukvy4Tb8tu6ND+S/AuITeGxIPfrKVhC6zG\ne8TzGxHZhIuVXzxMEs/1WKEVB4W0uET5zdThYnnq4x1Ct88Wr9jzP+tMFbvxHHESzBR/Jrq2h/2Y\nwAQmMIEJTGAC/x/4D847Sxknz33+AAAAAElFTkSuQmCC\n",
      "text/latex": [
       "$$4.33633323758e-09$$"
      ],
      "text/plain": [
       "4.33633323758e-09"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# time of flight (in seconds) for a straight path (infinite momentum)\n",
    "rEcal / c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtMAAALOCAYAAABrvJ90AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XuYZXdZJ/rvL51uArnQcUgEIRBCJIADJEZBVLQzAQfl\nEoVRjExMI8MlGpCHy4nCUYijBtQBFQ0KOKh4ZkAGRFAicpAcLgISTJOQBIZbJBAcgRCSoKFz+Z0/\nahcpmq7qXbveVauq9ufzPPupWrt2r/361T/eLL977dZ7DwAAsHoHjT0AAABsVpZpAACYkWUaAABm\nZJkGAIAZWaYBAGBGlmkAAJjRoMt0a+0RrbWPttY+3lo7Z5nX/ERr7bLW2kdaa//PkPMAAEClNtR9\npltr25J8LMnDknwuyQeTnN57v2LJa749yeuSnNJ7/0pr7U699y8OMhAAABQb8sr0g5J8ovd+Ze/9\npiSvTXLaPq95cpLf671/JUks0gAAbCZDLtN3TXLVkuPPTp5b6tuTnNBae09r7X2ttf844DwAAFDq\n4AHPPU1/5OAkxyf5wSTHJHlXa+3+i1eqAQBgIxtymf5cFhbkRcdk4er0Up9N8oHe+y1Jrmyt/e8s\nLNcfWvqi1towxW4AAFii995W8/ohax4XJfn21tqxrbUdSR6f5M37vOZNSXYlSWvtTknuneRT+ztZ\n792j6HHmmWeOPsNWechSnhv5IU9ZbtSHPOW5UR+zGGyZ7r3fnOTsJG9LcnmS1/Xer2itndtae/Tk\nNW9L8qXW2mVJ/i7Jc3rvXx5qJgAAqDRkzSO99wuSXLDPcy/Y5/jZSZ495Bx8o2OPPXbsEbYMWdaS\nZy151pFlLXnWkue4fAPiHNq1a9fYI2wZsqwlz1ryrCPLWvKsJc9xWaYBAGBGlmkAAJjRYF8nXqm1\n1jfDnAAAbF6ttfQNdGs8AADY0izTc+jCCy8ce4QtQ5a15FlLnnVkWUueteQ5Lss0AADMSGcaAACi\nMw0AAOvKMj2HdKvqyLKWPGvJs44sa8mzljzHZZkGAIAZ6UwDAEB0pgEAYF1ZpueQblUdWdaSZy15\n1pFlLXnWkue4LNMAADAjnWkAAIjONAAArCvL9BzSraojy1ryrCXPOrKsJc9a8hyXZRoAAGakMw0A\nANGZBgCAdWWZnkO6VXVkWUueteRZR5a15FlLnuOyTAMAwIx0pgEAIDrTAACwrizTc0i3qo4sa8mz\nljzryLKWPGvJc1yWaQAAmJHONAAARGcaAADWlWV6DulW1ZFlLXnWkmcdWdaSZy15jssyDQAAM9KZ\nBgCA6EwDAMC6skzPId2qOrKsJc9a8qwjy1ryrCXPcVmmAQBgRjrTAAAQnWkAAFhXluk5pFtVR5a1\n5FlLnnVkWUueteQ5Lss0AADMSGcaAACiMw0AAOvKMj2HdKvqyLKWPGvJs44sa8mzljzHZZkGAIAZ\n6UwDAEB0pgEAYF1ZpueQblUdWdaSZy151pFlLXnWkue4LNMAADAjnWkAAIjONAAArCvL9BzSraoj\ny1ryrCXPOrKsJc9a8hyXZRoAAGakMw0AANGZBgCAdWWZnkO6VXVkWUueteRZR5a15FlLnuOyTAMA\nwIx0pgEAIDrTAACwrizTc0i3qo4sa8mzljzryLKWPGvJc1yWaQAAmJHONAAARGcaAADWlWV6DulW\n1ZFlLXnWkmcdWdaSZy15jssyDQAAM9KZBgCA6EwDAMC6skzPId2qOrKsJc9a8qwjy1ryrCXPcVmm\nAQBgRjrTAAAQnWkAAFhXluk5pFtVR5a15FlLnnVkWUueteQ5Lss0AADMSGcaAACiMw0AAOvKMj2H\ndKvqyLKWPGvJs44sa8mzljzHZZkGAIAZ6UwDAEB0pgEAYF1ZpueQblUdWdaSZy151pFlLXnWkue4\nLNMAADAjnWkAAIjONAAArCvL9BzSraojy1ryrCXPOrKsJc9a8hyXZRoAAGakMw0AANGZBgCAdWWZ\nnkO6VXVkWUueteRZR5a15FlLnuOyTAMAwIx0pgEAIDrTAACwrizTc0i3qo4sa8mzljzryLKWPGvJ\nc1yWaQAAmNFcdKbPOSd54QuT29++biYAALaWWTrTc7FMf8u3JB//ePLv/l3hUAAAbCk+gLiM7duT\nm24ae4qNQ7eqjixrybOWPOvIspY8a8lzXHOzTO/dO/YUAABsNXNR8zjuuOTtb0/uda/CoQAA2FLU\nPJah5gEAwBAs03NIt6qOLGvJs5Y868iyljxryXNclmkAAJjRXHSmH/Sg5GUvSx784MKhAADYUnSm\nl+HKNAAAQ7BMzyHdqjqyrCXPWvKsI8ta8qwlz3FZpgEAYEaDdqZba49I8ttJtiV5Ve/9xfv8fXeS\n30zy2clTL+u9//f9nGdNnekf+ZHk534ueeQjZz4FAABb3Cyd6YMHHGZbkt9L8rAkn0vywdbam3vv\nVyx5WU/yP3vvzxhqjsSVaQAAhjFkzeNBST7Re7+y935TktcmOW2f17TJY1CW6W+kW1VHlrXkWUue\ndWRZS5615DmuIZfpuya5asnxZyfPLdWTPK619uHW2utba3cbYhDLNAAAQxisM91ae1ySR/Tenzw5\n/s9JHtx7f/qS13xLkut77ze11p6S5PG991P3c641daZ/+qeTU09Nzjxz5lMAALDFbajOdBZ60scs\nOT4mt33QMEnSe79myeEfJfmN5U62e/fuHHvssUmSnTt35sQTT8yuXbuS3Pb/3lju+ItfvDAf+UiS\nTPd6x44dO3bs2LFjx1v/eM+ePbn22muTJFdeeWVmMeSV6YOTfCzJqUmuTvIPSU5f+gHE1tqde+//\nPPn9x5I8t/f+vfs515quTD/tackDH5icddbMp9hSLrzwwq//HxJrI8ta8qwlzzqyrCXPWvKss6Gu\nTPfeb26tnZ3kbVm4Nd4f9d6vaK2dm+Si3vtbkjyjtfaYJDcn+VKS3UPMojMNAMAQBr3PdJW1Xpl+\n1rOSu941efazC4cCAGBLmeXK9EFDDbOR7NiR7N079hQAAGw1c7FMq3l8o8UCPmsny1ryrCXPOrKs\nJc9a8hyXZRoAAGY0F53p885LvvKV5EUvKhwKAIAtRWd6Ga5MAwAwBMv0HNKtqiPLWvKsJc86sqwl\nz1ryHJdlGgAAZjQXnelXvjL5wAeSV72qcCgAALYUnelluDINAMAQLNNzSLeqjixrybOWPOvIspY8\na8lzXJZpAACY0Vx0pv/iL5I/+ZPkTW8qHAoAgC1FZ3oZrkwDADAEy/Qc0q2qI8ta8qwlzzqyrCXP\nWvIcl2UaAABmNBed6Xe/O/nFX0ze857CoQAA2FJ0ppfhyjQAAEOYm2V6796xp9g4dKvqyLKWPGvJ\ns44sa8mzljzHNRfL9I4drkwDAFBvLjrTH/1octppycc+VjgUAABbis70MnSmAQAYgmV6DulW1ZFl\nLXnWkmcdWdaSZy15jssyDQAAM5qLzvQXv5iccELypS8VDgUAwJaiM70MV6YBABiCZXoO6VbVkWUt\nedaSZx1Z1pJnLXmOyzINAAAzmovOdO/JQQclt96atFW1YAAAmBc608toLTn4YFenAQCoNRfLdKLq\nsZRuVR1Z1pJnLXnWkWUtedaS57gs0wAAMKO56EwnyZ3ulFxxRXLUUUVDAQCwpehMr8CVaQAAqlmm\n55BuVR1Z1pJnLXnWkWUtedaS57gs0wAAMKO56Uzf5z7JG9+Y3O9+RUMBALCl6EyvYMcOV6YBAKg1\nN8u0msdtdKvqyLKWPGvJs44sa8mzljzHZZkGAIAZzU1n+qEPTX7t15If+IGioQAA2FJ0plfgyjQA\nANUs03NIt6qOLGvJs5Y868iyljxryXNclmkAAJjR3HSmf+zHkjPOSB772KKhAADYUnSmV+DKNAAA\n1SzTc0i3qo4sa8mzljzryLKWPGvJc1yWaQAAmNHcdKaf/OTku787ecpTioYCAGBL0ZlegSvTAABU\ns0zPId2qOrKsJc9a8qwjy1ryrCXPcVmmAQBgRnPTmf7FX0wOPzx53vOKhgIAYEvRmV7B9u3J3r1j\nTwEAwFYyN8v0jh1qHot0q+rIspY8a8mzjixrybOWPMc1N8u0zjQAANXmpjP9kpckV12VvPSlRUMB\nALCl6EyvwJVpAACqWabnkG5VHVnWkmctedaRZS151pLnuCzTAAAwo7npTP/pnyZvf3vymtcUDQUA\nwJaiM70CV6YBAKhmmZ5DulV1ZFlLnrXkWUeWteRZS57jskwDAMCM5qYz/da3Ji97WXLBBUVDAQCw\npehMr8CVaQAAqlmm55BuVR1Z1pJnLXnWkWUtedaS57gs0wAAMKO56Uz/wz8kP/dzyQc/WDQUAABb\nis70ClyZBgCg2lwt03v3jj3FxqBbVUeWteRZS551ZFlLnrXkOa65WqZdmQYAoNLcdKY//enklFOS\nK6+smQkAgK1FZ3oFrkwDAFDNMj2HdKvqyLKWPGvJs44sa8mzljzHZZkGAIAZzU1n+vrrk7vcJbnh\nhqKhAADYUnSmV+DKNAAA1SzTc0i3qo4sa8mzljzryLKWPGvJc1xzs0xv27bw85Zbxp0DAICtY246\n00lyu9slX/lKcsghBUMBALCl6EwfgKoHAACVLNNzSLeqjixrybOWPOvIspY8a8lzXJZpAACY0Vx1\npu961+QDH0judreCoQAA2FJ0pg/AlWkAACpZpueQblUdWdaSZy151pFlLXnWkue4LNMAADCjuepM\nP/CByR//cXLSSWufCQCArWWWzvTBK5zs5CQH2mBv6r1fupo3HNOOHa5MAwBQZ6Wax4VJ/tsBHm8Z\neL5Sah4LdKvqyLKWPGvJs44sa8mzljzHteyV6SQX9d5PWekft9beWTzPoCzTAABUmqvO9MMelpxz\nTvLwhxcMBQDAljLIfaZba9/fWjts8vsZrbWXttbuMeuQY3JlGgCAStPcGu/lSb7aWntgkmcl+USS\nPx10qoFYphfoVtWRZS151pJnHVnWkmcteY5rmmX65knH4keT/H7v/feTHD7sWMOwTAMAUOmAnenW\n2ruS/E2SJyZ5aJIvJNnTe7//8ON9fYaSzvRP/mRy2mnJ6acXDAUAwJYySGc6yeOTfC3Jz/Te/znJ\nXZP81gzzjc6VaQAAKh1wme69f773/t967++eHH+m9/4nw49WzzK9QLeqjixrybOWPOvIspY8a8lz\nXNPczeNxrbWPt9aua61dP3lctx7DVbNMAwBQaZrO9CeTPKr3fsX6jLTfGUo602efnZxwQvL0pxcM\nBQDAljJUZ/qfx1ykK7kyDQBApWmW6Ytaa69rrZ0+qXw8rrX22MEnG4BleoFuVR1Z1pJnLXnWkWUt\nedaS57imWabvmORfk/xQkkdNHo+e5uSttUe01j466Vyfs8LrHtdau7W19p3TnHdWlmkAACodsDM9\n84lb25bkY0keluRzST6Y5PR9KyOttcOT/HWSg5Oc3Xv/x/2cq6Qz/cIXJr0n55675lMBALDFlHam\nW2tPmeINV3rNg5J8ovd+Ze/9piSvTXLafl73X5O8KAv3sl7V8Ku1fXuyd++Q7wAAwDxZqebxC621\nxy7pSS99PLa19rgkz1zh3981yVVLjj87ee7rJrWOu/be3zp5apjL5BNqHgt0q+rIspY8a8mzjixr\nybOWPMd18Ap/e1cO3I3+2xX+tuJi3Fo7KMlLkpy59OkDvN+a7NhhmQYAoM6yy3Tvffcaz/25JMcs\nOT4mC1enFx2e5DuSXNhaS5I7J3lza+3R++tN7969O8cee2ySZOfOnTnxxBOza9euJLf9F9mBjrdv\n35Wbbpr+9Vv1ePG5jTLPZj7etWvXhppnsx/LU56OHTt2vJ7He/bsybXXXpskufLKKzOLIT+AeHAW\nPoB4apKrk/xD9vMBxCWvf2eSZw/5AcRXvCL54AeTV75yzacCAGCLGepLW2bSe785ydlJ3pbk8iSv\n671f0Vo7t7U21a31qulML1j8LzPWTpa15FlLnnVkWUueteQ5rpU602vWe78gyQX7PPeCZV57ypCz\nJJZpAABqLVvzaK09e8lhz20fDuxJ0nt/ybCjfcMsJTWP170uecMbkj//84KhAADYUmapeax0Zfrw\nLCzOJyT57iRvzsJC/ags9J83HVemAQCotGxnuvf+wt77uVm4C8d39t6f3Xt/VpKTk9xjvQasZJle\noFtVR5a15FlLnnVkWUueteQ5rmk+gHh0kqUr6E2T5zYdyzQAAJUOeGu81trzkzw+yRuzUPP40Szc\nmePXhx/v6zOUdKbf8Y7k134t+bu/KxgKAIAtpboznSTpvf9aa+1vkjw0Cx3q3b33i2eccVSuTAMA\nUGnZmkdr7VsWH0k+neQ1Sf4syT9Nntt0LNMLdKvqyLKWPGvJs44sa8mzljzHtdKV6X/M5DZ4y7hn\n8SyDs0wDAFBpsK8Tr1TVmf7wh5MzzkguuaRgKAAAtpRBvk68tXZQa+2M1tovT47v3lp70KxDjsmV\naQAAKk1za7zzkzwkyU9Njm+YPLfpWKYX6FbVkWUtedaSZx1Z1pJnLXmO64B380jy4N77Sa21i5Ok\n935Na237wHMNYvv2ZO/esacAAGCrmOY+0x9I8r1JLpos1Ucl+dve+0nrMeBkhpLO9NVXJyefnHz+\n8wVDAQCwpQzSmU7ysiR/keTo1tqvJ3lvkvNmmG90ah4AAFQ64DLde/+zJOdkYYG+Oslpvfc/H3qw\nIVimF+hW1ZFlLXnWkmcdWdaSZy15jmvZznRr7Yje+3WTL2j5P0n+5+RPvbX2Lb33a9ZlwkKWaQAA\nKi3bmW6t/XXv/ZGttSvzzV/e0nvvxw093JJZSjrTe/cmhx5qoQYA4JvN0pleaZn+/t77e1prh/Te\nbyyZcEZVy3TvyUEHJbfemrRVxQQAwFZX/QHE35n8/PvZR9pYWku2bUtuvnnsScalW1VHlrXkWUue\ndWRZS5615Dmule4zfXNr7ZVJ7tZa+90kS7f03nt/xrCjDWOxN719U94pGwCAjWSlmsdRSU5N8uIk\nv5yFZbov/uy9/8m6DVlU80iSI45IrroqueMdS04HAMAWMUvNY9kr0733LyR5bWvto733PWueboNw\nRw8AAKpMc5/pLbNIJ5bpRLeqkixrybOWPOvIspY8a8lzXNN8A+KWYpkGAKDKsp3pjaSyM33cccnb\n357c614lpwMAYIuovjXe4knPa60dueT4yNbar84y4EbgyjQAAFWmqXn8cO/9y4sHk98fOdxIw7JM\n61ZVkmUtedaSZx1Z1pJnLXmOa5pl+qDW2iGLB6212yfZMdxIw7JMAwBQ5YCd6dbaOUkek+S/Z+Ee\n009M8ube+4uHH+/rM5R1ph/0oOR3fzf5nu8pOR0AAFtE6X2mF/XeX9xauyTJw7LwpS2/0nt/24wz\njs6VaQAAqkx1a7ze+wW992f33p+zmRfpJNmxwzKtW1VHlrXkWUuedWRZS5615DmuZZfp1tp7Jz9v\naK1dv8/juvUbsZYr0wAAVJm7+0w/8pHJWWclj3pUyekAANgihrrP9GumeW6zcGUaAIAq03Sm//3S\ng9bawUlOHmac4VmmdasqybKWPGvJs44sa8mzljzHtVJn+nmtteuT3H9pXzrJvyR587pNWMwyDQBA\nlWnuM/2i3vsvrNM8y81Q1pk+88zklFOS3btLTgcAwBYxSGc6yfNaa2e01n558iZ3b609aKYJNwBX\npgEAqDLNMn1+kock+anJ8Q2T5zYly7RuVSVZ1pJnLXnWkWUtedaS57gO+A2ISR7cez+ptXZxkvTe\nr2mtbR94rsFYpgEAqDJNZ/oDSb43yUWTpfqoJH/bez9pPQaczFDWmX72s5O73CV5znNKTgcAwBYx\nVGf6ZUn+IsnRrbVfT/LeJOfNMN+G4Mo0AABVDrhM997/LMk5WVigr05yWu/9z4cebCiWad2qSrKs\nJc9a8qwjy1ryrCXPcU3TmU6S/53kusnre2vt7r33zww31nAs0wAAVJmmM/30JC/Iwpe13LL4fO/9\n/sOO9g0zlHWmzzsv+cpXkhe9qOR0AABsEbN0pqe5Mv3MJCf03r8021gbiyvTAABUmeYDiJ/JQsVj\nS9i+Pdm7d+wpxqVbVUeWteRZS551ZFlLnrXkOa5prkx/Osk7W2t/nWRxDe2995cMN9ZwduxwZRoA\ngBrTdKZfOPl18YUtC8v0uQPOte8MZZ3pV70qed/7kj/6o5LTAQCwRQzSme69v3DmiTYgnWkAAKoc\nsDPdWntLa+3Nk5+Lv7+mtfbzrbVD1mPISpZp3apKsqwlz1ryrCPLWvKsJc9xTfMBxE8nuSHJK5K8\nMsn1k+N7T443Fcs0AABVpulMX9R7/679Pddau6z3/h2DTpjazvSb3pS8+tXJX/5lyekAANgiZulM\nT3Nl+tDW2j2WvMk9khw6Odx0N5lzZRoAgCrTLNPPTvLu1tqFrbULk7w7yXNba4cm+ZMhhxuCZVq3\nqpIsa8mzljzryLKWPGvJc1zT3M3jra21eye5TxZuj/ex3vuNkz//9pDDDcEyDQBAlQN2ppOktXb/\nJPdLckgm95vuvf/psKN9w/uXdabf857knHOS97635HQAAGwRg9xnevKlLT+Y5DuS/HWSH07yniTr\ntkxXcmUaAIAq03Sm/1OShyX5fO/9iUkemGTnoFMNyDKtW1VJlrXkWUuedWRZS5615DmuaZbpf+u9\n35Lk5tbaHZP8S5Jjhh1rOJZpAACqTHOf6fOTPD/J47NwZ4+vJrl4cpV6XVR2pj/60eS005KPfazk\ndAAAbBGDdKZ77z87+fUPWmtvS3J47/2SWQbcCFyZBgCgyjQ1j7TWHthaOy3JSUm+vbX22GHHGo5l\nWreqkixrybOWPOvIspY8a8lzXNPczePVSe6f5LIkty750xuHGmpI27cnezfd9zYCALARTdOZvjzJ\nd5SVlmdQ2Zn+4heTe987ueaaktMBALBFzNKZnqbm8f4sfGHLlrBjh5oHAAA1plmm/yTJ37fW/ndr\n7dLJwwcQNzHdqjqyrCXPWvKsI8ta8qwlz3EdsDOd5L8nOSPJR/KNnelNyTINAECVaTrT7+u9P2Sd\n5lluhrLOdO/JQQclt9yy8BMAAJKB7jOd5OLW2v9I8pYki/fB6L33TXk3j9Zuuzp9u9uNPQ0AAJvZ\nNNdm75Dka0l+KMmjJo9HDznU0Oa96qFbVUeWteRZS551ZFlLnrXkOa5pvgFx9zrMsa7mfZkGAKDG\nATvTG0FlZzpJjjoqueyy5Oijy04JAMAmN9R9prccV6YBAKhgmZ5DulV1ZFlLnrXkWUeWteRZS57j\nOuAy3Vrb2Vp7aWvtQ5PHf2ut3XE9hhvKvC/TAADUmOY+029McmkWvgmxZeELXB7Qe3/s8ON9fYbS\nzvR975u84Q3J/bbMl6QDALBWQ91n+l77LM4vbK19eHWjbSyuTAMAUGGazvS/tdYeunjQWvv+JP86\n3EjDm/dlWreqjixrybOWPOvIspY8a8lzXNNcmX5akj9d0pP+cpIzhxtpePO+TAMAUGPFznRrbVuS\nF/fen7O4TPfev7Jewy2Zo7Qz/dCHJr/2a8kP/EDZKQEA2OTKO9O991taa9/fFrbZdV+ih7J9e7J3\n79hTAACw2U3Tmd6T5C9ba2e01h43eazbnTyGMO81D92qOrKsJc9a8qwjy1ryrCXPcU3TmT4kyTVJ\n/sM+z7+xfpz1sWPHfC/TAADUOOB9pjeC6s70Yx+bPOEJyeMeV3ZKAAA2uVk609N8A+IJrbV3tNYu\nmxw/oLX2f8865EYw7zUPAABqTNOZfmWS5yVZ/MjepUlOH2yidTDvy7RuVR1Z1pJnLXnWkWUtedaS\n57imWabv0Hv/wOLBpG+xqVfReV+mAQCoccDOdGvtgiRPT/L63vtJrbX/lORJvfcfXo8BJzOUdqaf\n8pTk5JOTpz617JQAAGxy5feZnjg7ySuSnNBauzrJp5M8YYb5NgxXpgEAqHDAmkfv/ZO991OTHJ3k\nPr337+u9Xzn4ZAOa92Vat6qOLGvJs5Y868iyljxryXNc09zN49bW2ouTfLX3ft3kuX8cfLIBzfsy\nDQBAjWk605cmuSDJdyZ5fO/9S621i3vvJ63HgJMZSjvTz3tecuihyfOfX3ZKAAA2uUHuM53k5t77\n/5WFW+S9u7V28kzTbSCuTAMAUGGaZTpJ0nt/XZKfSPLHSY4baqD1MO/LtG5VHVnWkmctedaRZS15\n1pLnuKa5m8eTF3/pvX+ktfbQJKcNN9Lwtm9Prrtu7CkAANjsDtiZ3giqO9MveUly1VXJS19adkoA\nADa5oTrTW8681zwAAKgxt8v03r1jTzEe3ao6sqwlz1ryrCPLWvKsJc9xTXOf6XdM89xm4so0AAAV\nlu1Mt9Zun+QOSd6ZZNeSPx2R5G967/cZfLrbZintTL/mNcnb3pb82Z+VnRIAgE1uls70SnfzeGqS\nn0/ybUk+tOT565P83urH2zhcmQYAoMKyNY/e+2/33u+Z5Lm993sueTyg926Z3sR0q+rIspY8a8mz\njixrybOWPMd1wM507/13932utXbnaU7eWntEa+2jrbWPt9bO2c/fn9Zau6S1dnFr7d2ttftON/ba\nzPsyDQBAjZnuM91a++ve+yMP8JptST6W5GFJPpfkg0lO771fseQ1h/fer5/8/ugkP9t7/+H9nKu0\nM33BBcnv/E7yN39TdkoAADa58vtMt9YObq29c9/nD7RITzwoySd671f23m9K8trs882Ji4v0xGFJ\nbp3ivGvmyjQAABVWXKZ77zcnubW1tnOGc981yVVLjj87ee4btNZ+trX2iSQvTvKMGd5n1eZ9mdat\nqiPLWvKsJc86sqwlz1ryHNdKd/NY9NUkl7bW/jbJv06e6733Ay2+U/Uyeu/nJzm/tXZ6kl9Ksnua\nf7cW875MAwBQY5pl+o2Tx+Jy3DLdovy5JMcsOT4mC1enl/O6JC9f7o+7d+/OsccemyTZuXNnTjzx\nxOzatSvJbf9FNu3xJZdcmGuuSRZvn73af7/Zjxef2yjzbObjXbt2bah5NvuxPOXp2LFjx+t5vGfP\nnlx77bXsThUtAAAgAElEQVRJkiuvvDKzmOkDiFOduLWDs/ABxFOTXJ3kH/LNH0A8vvf+icnvj07y\nS733B+3nXKUfQNyzJznzzOTDHy47JQAAm1z5BxBXeKNzD/SaSd/67CRvS3J5ktf13q9orZ07WZyT\n5OzW2kdaaxcneWaSM2eZZ7Xmveax+F9mrJ0sa8mzljzryLKWPGvJc1zT1Dz256JpXtR7vyDJBfs8\n94Ilvz9zxvdfk3lfpgEAqHHAmkdr7ZNJ3p/k3Une3Xu/bD0G22eG0prHpz+dnHJKMmM1BgCALWio\nmsd3JHlFkn+X5Ldaa59srb1plgE3ClemAQCoMM0yfXOSm5LckoUvVflCkv8z5FBDm/dlWreqjixr\nybOWPOvIspY8a8lzXNN0pq9LcmmSlyR5Ve/9i8OONLzt25O9e8eeAgCAzW6azvRpSR6a5LuzcIX6\n75O8q/f+/w4/3tdnKO1MX399cuc7J1/9atkpAQDY5GbpTE99n+nW2n2S/EgWbmF3dO/9kNWPOJvq\nZfrGG5MjjnB1GgCA2wzyAcTW2hsmd/T43SR3SHJGkiNnG3FjWOxMD/R9NRueblUdWdaSZy151pFl\nLXnWkue4pulMvyjJxZMvYdkStm1LDjooueWW5OBZ77QNAMDcG+zrxCtV1zyS5JBDki9/Obn97UtP\nCwDAJrVuXye+Fcz77fEAAFg7y/Qc0q2qI8ta8qwlzzqyrCXPWvIc1zQfQDyotXZGa+2XJ8d3b609\naPjRhjXPyzQAADWmuc/0H2Thmw//Q+/9Pq21b0nyt73371qPASczlHem73a35H3vS445pvS0AABs\nUrN0pqe5l8WDe+8ntdYuTpLe+zWtte0zTbiBuDINAMBaTdOZ3tta27Z40Fo7KgtXqje1eV6mdavq\nyLKWPGvJs44sa8mzljzHNc0y/bIkf5Hk6Nbaryd5b5LzBp1qHczzMg0AQI2p7jPdWrtvklMnh+/o\nvV8x6FTf/P7lnekTT0xe/erkpJNKTwsAwCY1SGe6tfY9SS7vvf/e5PiI1tqDe+8fmHHODcGVaQAA\n1mqamscfJLl+yfFXJ89tavO8TOtW1ZFlLXnWkmcdWdaSZy15jmuqL21Z2rHovd+SZNsKL98U5nmZ\nBgCgxjT3mf6LJO9M8vIkLclZSU7pvf/o8ON9fYbyzvTDHpacc07y8IeXnhYAgE1qls70NFemn5bk\n+5J8Lslnk3xPkqesfryNZfv2ZO/esacAAGAzO+Ay3Xv/P733x/fej548Tu+9/8t6DDekea556FbV\nkWUtedaSZx1Z1pJnLXmOa5q7eRyd5MlJjl3y+t57/5kB5xrcjh3zu0wDAFBjms70+5K8K8mHcts3\nH/be+xsGnm3pDOWd6dNPTx796OSnfqr0tAAAbFKD3Gc6ye177+fMONOGNc81DwAAakzzAcS/aq09\ncvBJ1tk8L9O6VXVkWUueteRZR5a15FlLnuOaZpl+ZpK3tNZubK1dP3lcN/RgQ5vnZRoAgBoH7Exv\nBEN0pp/+9OTbvz15xjNKTwsAwCY1yH2mW2sHtdbOaK398uT47q21B8065EbhyjQAAGs1Tc3j/CQP\nSbJ434sbJs9tavO8TOtW1ZFlLXnWkmcdWdaSZy15jmuau3k8uPd+Umvt4iTpvV/TWts+8FyDm+dl\nGgCAGtPcZ/oDSb43yUWTpfqoJH/bez9pPQaczFDemT733OSWW5Jf+ZXS0wIAsEkN0plO8rIkf5Hk\n6Nbaryd5b5LzZphvQ3FlGgCAtVpxmW6tHZTk00nOycICfXWS03rvf74Osw1qnpdp3ao6sqwlz1ry\nrCPLWvKsJc9xrdiZ7r3f2lr7/d77iUmuWKeZ1sU8L9MAANSYpjP9W0nen+QN5cXlKQ3Rmf79308u\nuyw5f9PflwQAgApDdaafluTPk+z1DYgAAHCbAy7TvffDeu8H9d63994PnzyOWI/hhrR9e7J379hT\njEO3qo4sa8mzljzryLKWPGvJc1y+AREAAGY0TWf6D5LcmuQ/9N7v01r7lizcZ/q71mPAyQzlnenX\nvS75X/8ref3rS08LAMAmNUtnem6/AXHHDlemAQBYm2k+gLi3tbZt8WDyDYi3DjfS+pjnmoduVR1Z\n1pJnLXnWkWUtedaS57h8AyIAAMxo2c50a+243vunJr/fN8mpkz+9o/e+rl/gMkRn+u/+Lvmv/zV5\n5ztLTwsAwCZV3Zl+fZKTW2vv6L2fGt+ACAAA32Clmse21trzk5zQWntWa+3ZSx7PWq8BhzLPy7Ru\nVR1Z1pJnLXnWkWUtedaS57hWWqYfn+SWJNuSHJ7ksCWPw4cfbVjzvEwDAFBjpc70z/fef6e19su9\n919Z57n2naW8M33JJckTnpBcemnpaQEA2KRm6UyvdGX6ZyY/f2z2kTYuV6YBAFirlZbpy1trH89C\nZ/rSfR6XrNeAQ5nnZVq3qo4sa8mzljzryLKWPGvJc1zL3s2j9356a+3OSf42yaOTrOqS90Y3z8s0\nAAA1lu1MbyRDdKavvjo5+eTk858vPS0AAJtU6X2mW2uv773/eGttfx/R6733B6x6wg3ElWkAANZq\npc70z09+Pno/j8cMPNfg5nmZ1q2qI8ta8qwlzzqyrCXPWvIc10qd6asnP69ct2nW0Twv0wAA1Fjp\nPtM3JFmuqNx770cMNtU3z1Lemd67N7nDHZKbby49LQAAm1RpZ7r3ftjkpL+a5Ookfzb50xOSfNus\nQ24U27cnt9yS9J60LXWfEgAA1stKnelFj+m9n997v27yeHmS04YebGitJQcfPJ9VD92qOrKsJc9a\n8qwjy1ryrCXPcU2zTH+1tfafW2vbJo8nJLlh6MHWg940AABrccD7TLfW7pnkd5J87+Sp9yb5+fX8\nYOIQnekkueMdk3/6p2TnzvJTAwCwyZR2phf13j+dLXArvP1xZRoAgLWYpuaxZc3rMq1bVUeWteRZ\nS551ZFlLnrXkOS7L9Bwu0wAA1JimM31c7/1TB3puSEN1pu91r+Rtb0uOP7781AAAbDKzdKanuTL9\nhv089/rVvMlG5co0AABrsewy3Vq7b2vtcUnu2Fp7bGvtcZOfu5Mcsm4TDmhel2ndqjqyrCXPWvKs\nI8ta8qwlz3GtdDePeyd5dJI7Tn4uuj7Jk4ccar3M6zINAECNaTrT39t7//t1mme5GQbpTD/4wcnv\n/E7yPd9TfmoAADaZQe4zneQTrbXnJzl2yet77/1nVjnfhuPKNAAAazHNBxD/MskRSd6e5K+XPDa9\neV2mdavqyLKWPGvJs44sa8mzljzHNc2V6dv33s8ZfJIRzOsyDQBAjWk607+a5H2999GuRg/VmX7k\nI5Ozzkoe9ajyUwMAsMmUdqZbazckWdxgn9da25tk8Tpu770fMduYG8f27cnevWNPAQDAZrVsZ7r3\nfljv/fDJ46De+yFLjjf9Ip3Mb81Dt6qOLGvJs5Y868iyljxryXNcB+xMt9a+cz9PfyXJP/Xeb64f\naf3M6zINAECNaTrT709ycpJLJk/dP8llWfgyl7N6728bdMIM15nevTv5wR9MnvjE8lMDALDJzNKZ\nnubWeFcnObH3fnLv/eQkJyb5VJKHJ/mN1Y+5cbgyDQDAWkyzTJ/Qe79s8aD3fnmS+/TeP5nbPqC4\nKc3rMq1bVUeWteRZS551ZFlLnrXkOa5p7jN9WWvt5Ulem6Ql+Ykkl7fWbpfb7u6xKc3rMg0AQI1p\nOtN3SPKzSb5v8tR7k5yf5MYkh/berx90wgzXmX7Oc5Jv/dbkuc8tPzUAAJtM6X2mF/Xe/zXJb00e\n+xp8kR6SK9MAAKzFsp3p1trrJz8v3c/jkuX+3WYyr8u0blUdWdaSZy151pFlLXnWkue4Vroy/fOT\nn49ej0HGsH178rWvjT0FAACb1QE70xvBUJ3pF70o+fKXkxe/uPzUAABsMqWd6dbaDVn+1nd9K3yl\n+LzWPAAAqLHSfaYf0Hs/fJnHpl+kk/ldpnWr6siyljxrybOOLGvJs5Y8x7XSMr34AcR3rNMs625e\nl2kAAGos25lure3JwkJ9VpKXZOELWxb13vtLhh/v67MM0pl+1auS970v+aM/Kj81AACbzCyd6ZWu\nTP9kkluSbEtyeJLDljwOn3XIjcSVaQAA1mLZZbr3/tHe+4uSPKn3fu6+j3WccTA7diR79449xfrT\nraojy1ryrCXPOrKsJc9a8hzXSlemkyS997euxyBjOOyw5IYbxp4CAIDNaq7vM/3udye/+IvJe95T\nfmoAADaZ0s50a+3HJz+PW+tgG9WRRy58aQsAAMxipZrH8yY/37Aeg4xh587k2mvHnmL96VbVkWUt\nedaSZx1Z1pJnLXmOa9lvQEzypdba25Pcs7X2ln3+1nvvjxlwrnUxr8s0AAA1VrrP9I4k35nkz5I8\nKd98n+n/b/jxvj7LIJ3p3hfu6PHVry78BABgfs3SmT7gBxBba0f13r/QWjssSXrv637/i6GW6SQ5\n6qjkssuSo48e5PQAAGwS1V/asujOrbWLk1ye5PLW2odaa/9+pgk3oHmseuhW1ZFlLXnWkmcdWdaS\nZy15jmuaZfoVSZ7Ve7977/3uSZ49eW5LmMdlGgCAGtPUPD7ce3/ggZ4b0pA1j4c/PHnuc5Mf+qFB\nTg8AwCYxS81jpbt5LPp0a+2XkrwmCx9CfEKST80w34Z05JGuTAMAMJtpah4/k+ToJG/Mwj2nj5o8\ntyXs3Dl/X9yiW1VHlrXkWUuedWRZS5615DmuA16Z7r1fk+Tps75Ba+0RSX47ybYkr+q9v3ifvz8r\nC7feuznJF5L8TO/9M7O+32rpTAMAMKsDdqbXdPLWtiX5WJKHJflckg8mOb33fsWS1+xK8v7e+42t\ntacl2dV7/8l9zjNYZ/rXfz25/vrkvPMGOT0AAJvEULfGW4sHJflE7/3K3vtNSV6b5LSlL+i9X9h7\nv3Fy+IEkdxt4pm+gMw0AwKyGXqbvmuSqJcefnTy3nCcleeugE+1DZ5q1kGUtedaSZx1Z1pJnLXmO\n64DLdGvthNbaO1prl02OH9Ba+7+nPP/U3YzW2n/OwteX/+a0/6aCzjQAALOa5tZ4r0zy3CR/MDm+\nNMn/TPKrU/zbzyU5ZsnxMVm4Ov0NWmsPS/K8JD8wqYN8k927d+fYY49NkuzcuTMnnnhidu3aleS2\n/yKb5XjnzuQzn7kwF14427/fjMeLz22UeTbz8a5duzbUPJv9WJ7ydOzYseP1PN6zZ0+unVxVvfLK\nKzOLab605aLe+3e11i7uvZ80eW5P7/3EA568tYOz8AHEU5NcneQf8s0fQDwpyeuT/Mfe+yeXOc9g\nH0D86EeTH/3RhZ8AAMyvoT6A+IXW2vFL3uQ/Jfn8NCfvvd+c5Owkb0tyeZLX9d6vaK2d21p71ORl\nv5Hk0CT/q7V2cWvtTav5H2CtdKZZC1nWkmctedaRZS151pLnuKapeZyd5BVJ7tNauzrJp7PwLYhT\n6b1fkOSCfZ57wZLfHz7tuYaw2JnuPWmr+u8QAADm3dT3mW6tHZrkoN779cOOtN/3HqzmkSSHHLJw\ndfr2tx/sLQAA2OBmqXkc8Mp0a+3IJD+d5NgkB7eFy7e99/6MWYbciBbvNW2ZBgBgNabpTL81yT2S\nXJLkoiQfmjy2jHnrTetW1ZFlLXnWkmcdWdaSZy15jmuazvTteu/PGnySEbnXNAAAs5jm1njPSnJD\nkrck+dri8733a4Yd7RtmGLQz/cM/nDz96cmP/MhgbwEAwAY3SGc6Cwv0byR5fpJbJ8/1JMetbryN\na7EzDQAAqzFNZ/o5SY7vvd+j937PyWPLLNLJ/NU8dKvqyLKWPGvJs44sa8mzljzHNc0y/fEk/zb0\nIGOatw8gAgBQY5rO9JuSfEeSd+a2zvS63hpv6M70b/xG8oUvJL/5m4O9BQAAG9xQnek3TR5LDbfZ\njuDII5OPf3zsKQAA2GwOWPPovf/xfh5/sh7DrRedaWYly1ryrCXPOrKsJc9a8hzXslemW2uv773/\neGvt0v38uffeHzDgXOtKZxoAgFks25lurX1b7/3q1to9kuzbHem9938afLrbZhm0M/3BDyZnnZVc\ndNFgbwEAwAY3S2d62ZpH7/3qya8/23u/cukjyc+uYc4Nx32mAQCYxTS3xvuh/Ty3pb4rUGeaWcmy\nljxrybOOLGvJs5Y8x7VSZ/qsLFyBvtc+venDk7x36MHW0x3vuLBM9560VV3YBwBgnq3Umb5jkiOT\nvCjJObmtN3197/1L6zPe12cZtDOdJIcdlnz+88nhhw/6NgAAbFCl95nuvX8lyVeS/ORaB9sMFqse\nlmkAAKY1TWd6LszThxB1q+rIspY8a8mzjixrybOWPMdlmZ5wr2kAAFZr2c70RrIenelHPzp58pOT\nxzxm0LcBAGCDKr3P9LyZt9vjAQCwdpbpCZ1pZiHLWvKsJc86sqwlz1ryHJdlesKVaQAAVktneuIl\nL0muuip56UsHfRsAADYonek1cGUaAIDVskxP6EwzC1nWkmctedaRZS151pLnuCzTE65MAwCwWjrT\nExdfnDzxicmePYO+DQAAG5TO9Bq4Mg0AwGpZpid0ppmFLGvJs5Y868iyljxryXNclumJI45Ibrgh\nufXWsScBAGCz0JleYufO5NOfXrhKDQDAfNGZXiO9aQAAVsMyvcS89KZ1q+rIspY8a8mzjixrybOW\nPMdlmV7ClWkAAFZDZ3qJH/ux5Iwzksc+dvC3AgBgg9GZXiNXpgEAWA3L9BI606yWLGvJs5Y868iy\nljxryXNcluklXJkGAGA1dKaX+N3fTT7xiYWfAADMF53pNdq5M/nyl8eeAgCAzcIyvYTONKsly1ry\nrCXPOrKsJc9a8hyXZXoJnWkAAFZDZ3qJSy9NfuqnFn4CADBfdKbXSGcaAIDVsEwvoTPNasmyljxr\nybOOLGvJs5Y8x2WZXuLQQ5OvfS256aaxJwEAYDPQmd7Hne6UfPSjCz8BAJgfOtMF9KYBAJiWZXof\n89Cb1q2qI8ta8qwlzzqyrCXPWvIcl2V6H+41DQDAtHSm9/HjP578xE8s/AQAYH7oTBdwZRoAgGlZ\npvdx5JFb/wOIulV1ZFlLnrXkWUeWteRZS57jskzvw5VpAACmpTO9j/PPTz7ykYWfAADMD53pAq5M\nAwAwLcv0PnSmWQ1Z1pJnLXnWkWUtedaS57gs0/twZRoAgGnpTO/jiiuSxz524ScAAPNDZ7qAK9MA\nAEzLMr0PnWlWQ5a15FlLnnVkWUueteQ5Lsv0Pg45ZOHnjTeOOwcAABufzvR+3PnOyZ49Cz8BAJgP\nOtNF9KYBAJiGZXo/tnpvWreqjixrybOWPOvIspY8a8lzXJbp/XBlGgCAaehM78fppyePeczCTwAA\n5oPOdBFXpgEAmIZlej+2+jKtW1VHlrXkWUuedWRZS5615Dkuy/R+bPUPIAIAUENnej9e8YrkoosW\nfgIAMB90pots9ZoHAAA1LNP7sdWXad2qOrKsJc9a8qwjy1ryrCXPcVmm90NnGgCAaehM78fHP578\nyI8s/AQAYD7oTBfZ6jUPAABqWKb3Y3GZ3gQX7WeiW1VHlrXkWUuedWRZS5615Dkuy/R+bN+e3O52\nyVe/OvYkAABsZDrTy7jb3ZL3v3/hJwAAW5/OdCG9aQAADsQyvYydO5Nrrhl7imHoVtWRZS151pJn\nHVnWkmcteY7LMr2Mu989+cxnxp4CAICNTGd6Gb/0S8nBBycveMG6vi0AACPRmS50/PHJJz4x9hQA\nAGxklull3OteySc/OfYUw9CtqiPLWvKsJc86sqwlz1ryHJdlehn3upcr0wAArExnehm9J4cfnlx9\ndXLEEev61gAAjEBnulBryXHHbd2qBwAAa2eZXsHxx2/NZVq3qo4sa8mzljzryLKWPGvJc1yW6RVs\n5Q8hAgCwdjrTK/jDP0wuuih55SvX/a0BAFhnOtPFXJkGAGAllukVbNXb4+lW1ZFlLXnWkmcdWdaS\nZy15jssyvYJjjkn+5V+SG28cexIAADYinekDuPe9k7/8y+S+9x3l7QEAWCc60wPYqrfHAwBg7SzT\nB7AVe9O6VXVkWUueteRZR5a15FlLnuOyTB+AO3oAALAcnekD+Ku/Ss4/P3nrW0d5ewAA1onO9AC2\nYs0DAIAagy/TrbVHtNY+2lr7eGvtnP38/Qdaa//YWruptfa4oedZrXveM/nMZ5Kbbx57kjq6VXVk\nWUueteRZR5a15FlLnuMadJlurW1L8ntJHpHkfklOb63te5O5f0pyZpL/MeQsszrkkOToo5Orrhp7\nEgAANppBO9OttYckeUHv/RGT419Ikt77i/bz2lcn+ave+xv287fROtNJcsopyfOfnzzsYaONAADA\nwDZiZ/quSZZe0/3s5LlN5fjj9aYBAPhmQy/TG/9WIVPYarfH062qI8ta8qwlzzqyrCXPWvIc18ED\nn/9zSY5ZcnxMFq5OL2fZ5Xv37t059thjkyQ7d+7MiSeemF27diW57f+Ihjr+2tcuzPvfnyTr835D\nH+/Zs2dDzePYsWPHG/140UaZZ7MfL9oo82z240UbZZ7NdLxnz55ce+21SZIrr7wysxi6M31wko8l\nOTXJ1Un+Icnpvfcr9vPaP07ylo3Ymb744uTMM5NLLhltBAAABjZLZ3rwL21prf1wkt9Osi3JH/Xe\nz2utnZvkot77W1pr353kjUmOTHJjks/33u+/zzlGXaavuy65y12SG25I2qriBQBgs9iIH0BM7/2C\n3vsJvffje+/nTZ57Qe/9LZPfP9h7P6b3fljv/U77LtIbwRFHJIcemvzzP489SY19/99CzE6WteRZ\nS551ZFlLnrXkOa7Bl+mtYqt9CBEAgLUbvOZRYeyaR5KccUZy6qnJ7t2jjgEAwEA2ZM1jq3BlGgCA\nfVmmp7SVlmndqjqyrCXPWvKsI8ta8qwlz3FZpqfkWxABANiXzvSU/uVfkvveN/nSl0YdAwCAgehM\nD+ioo5Kbbkq+/OWxJwEAYKOwTE+pta3Tm9atqiPLWvKsJc86sqwlz1ryHJdlehX0pgEAWEpnehV+\n4ReSww9Pnv/8sScBAKCazvTA7nUvV6YBALiNZXoVjj9eZ5pvJMta8qwlzzqyrCXPWvIcl2V6FbbK\nBxABAKihM70Kt96aHHrowr2m73CHsacBAKCSzvTADjooOfbY5FOfGnsSAAA2Asv0Km2FDyHqVtWR\nZS151pJnHVnWkmcteY7LMr1KW+VDiAAArJ3O9Cq97GXJ5ZcnL3/52JMAAFBJZ3oduDINAMAiy/Qq\nnXBCctllyQa5UD4T3ao6sqwlz1ryrCPLWvKsJc9xWaZX6Z73XPjp6jQAADrTM3jCE5JTTkn+y38Z\nexIAAKroTK+TU05J3vnOsacAAGBslukZnHJKcuGFm7c3rVtVR5a15FlLnnVkWUueteQ5Lsv0DI47\nbuHbED/+8bEnAQBgTDrTM/rpn06+7/uSpz517EkAAKigM72Odu1aqHoAADC/LNMzWvwQ4ga7YD4V\n3ao6sqwlz1ryrCPLWvKsJc9xWaZndM97JoccknzsY2NPAgDAWHSm12D37uTBD07OOmvsSQAAWCud\n6XXmftMAAPPNMr0Gix9C3IAXzVekW1VHlrXkWUuedWRZS5615Dkuy/Qa3OMeyWGHJZdfPvYkAACM\nQWd6jZ70pOSkk5Kzzx57EgAA1kJnegTuNw0AML8s02t0yikLy/Stt449yfR0q+rIspY8a8mzjixr\nybOWPMdlmV6ju90tOfLI5LLLxp4EAID1pjNd4MlPTu5//+QZzxh7EgAAZqUzPRL3mwYAmE+W6QK7\ndiXvetfm6U3rVtWRZS151pJnHVnWkmcteY7LMl3g274tudOdkksuGXsSAADWk850kac9LbnPfZJn\nPnPsSQAAmIXO9Ih27dKbBgCYN5bpIou96VtuGXuSA9OtqiPLWvKsJc86sqwlz1ryHJdlusid77zw\n+PCHx54EAID1ojNd6Oyzk7vcJXn+88eeBACA1ZqlM22ZLvShDyWPe1zyyU8m27aNPQ0AAKvhA4gj\nO/nk5Fu/NbnggrEnWZluVR1Z1pJnLXnWkWUtedaS57gs08XOOit5+cvHngIAgPWg5lHs3/4tufvd\nkw98IDnuuLGnAQBgWmoeG8Dtb5+ceWbyh3849iQAAAzNMj2Apz0tefWrkxtvHHuS/dOtqiPLWvKs\nJc86sqwlz1ryHJdlegDHH59853cmr3/92JMAADAknemBvPnNyXnnJe9739iTAAAwDZ3pDeSRj0w+\n97nk4ovHngQAgKFYpgeybVvy1KduzNvk6VbVkWUtedaSZx1Z1pJnLXmOyzI9oCc9aaE3fe21Y08C\nAMAQdKYHdvrpyUMekjzjGWNPAgDASmbpTFumB/budydPfnJyxRVJW9X/agAAWE8+gLgBff/3J9u3\nJ+9859iT3Ea3qo4sa8mzljzryLKWPGvJc1yW6YG1lpx1VnL++WNPAgBANTWPdXDddclxxyXveldy\nv/uNPQ0AAPuj5rFBHXFEcu65yVOektx669jTAABQxTK9Tp72tOTmm5NXvWrsSXSrKsmyljxrybOO\nLGvJs5Y8x2WZXifbtuX/b+/Oo+woyzyOf3/pkJCNJaw5LAIxkWVwkhhACEgAwcA4ICOKjOKoiHCQ\ngcNxGcCZIRxXBhE4ZgYVUBkUwiIg6NHAsIwYCSGSsC8RbAZiWAwEskDWZ/5469KVy+3ue2+qu/p2\n/z7nvKfqvvVW1XufvEmern6rih/9CL72NXjxxbJ7Y2ZmZmZF8JzpXnb22dDeDjNnlt0TMzMzM8vz\nc6ZbwMqVsPfeMGMGHHlk2b0xMzMzswrfgNgChg+HH/wATjsNVqwopw+eW1Ucx7JYjmexHM/iOJbF\ncjyL5XiWy8l0CQ4/HKZMgenTy+6JmZmZmW0MT/Moycsvp+kes2bBhAll98bMzMzMPM2jhWy7LXz7\n2zZWeZsAABIcSURBVHDyybBuXdm9MTMzM7NmOJku0Wc/CyNGwKWX9u55PbeqOI5lsRzPYjmexXEs\ni+V4FsvxLNfgsjswkElw5ZVw4IEwdiwcc0zZPTIzMzOzRnjOdB8wb156TN5NN8FBB5XdGzMzM7OB\nyXOmW9TkyXDNNXDccfDII2X3xszMzMzq5WS6jzj8cLjkEjjqKHjuuZ49l+dWFcexLJbjWSzHsziO\nZbEcz2I5nuXynOk+5IQT0iPzjjgCZs+Grbcuu0dmZmZm1hXPme6Dzj0X7rwzlZEjy+6NmZmZ2cDQ\nzJxpJ9N9UAR8/vOwaBHceisMGVJ2j8zMzMz6P9+A2E9I8MMfwqhRcMghsHhxscf33KriOJbFcjyL\n5XgWx7EsluNZLMezXE6m+6jBg+G662DaNNhnH7jvvrJ7ZGZmZmbVPM2jBfzqV/C5z8HXvw6nnFJ2\nb8zMzMz6J8+Z7seefhqOPRamTIHvfx+GDi27R2ZmZmb9i+dM92Pjx8OcObBkCRx8cLo5sVmeW1Uc\nx7JYjmexHM/iOJbFcjyL5XiWy8l0Cxk1Cm68EY45BiZNgssvh3Xryu6VmZmZ2cDlaR4tasECOP10\nWLUKZsyA/fYru0dmZmZmrc3TPAaQCRPg3nvhjDPSXOqTToJXXim7V2ZmZmYDi5PpFibBiSfCk0/C\nFlvAXnulq9Rr13a9n+dWFcexLJbjWSzHsziOZbEcz2I5nuVyMt0PbLYZXHQR3H033HJLullxxgxY\nsaLsnpmZmZn1b54z3Q/94Q/w3e+maSCnnprmVm+3Xdm9MjMzM+vbPGfaADjgALjpppRU//WvsMce\n6WUvTz1Vds/MzMzM+hcn0/3YuHFw2WUpiR4zBj7wATjwQDjrrHtYsqTs3vUPnqdWLMezWI5ncRzL\nYjmexXI8y+VkegDYZhuYPh1eeAHOOQceegjGjoWjj4brr4c33yy7h2ZmZmatyXOmB6hly+Dmm+Hn\nP4e5c+Goo1L50Idg663L7p2ZmZlZ72tmzrSTaWPxYrjtNvjNb9ITQd7zHjjyyFQmT4a2trJ7aGZm\nZtbzfAOi1aV6btWYMfCFL6Qr1S+/DN/5Tnqs3kknpaeAHHtsejrInDmwenU5fe6rPE+tWI5nsRzP\n4jiWxXI8i+V4lmtw2R2wvmXIEDjkkFQuvDDNs773Xpg9O00JWbgQJk2CKVNg//3T+g47pBfImJmZ\nmQ00nuZhDXnjjXSFevbstJw/HyJg4sT0ivOJE1MZN87TQ8zMzKy1eM609bqINOd6/vxUFixIy7/8\nBd79bth99/Sc68py/HgYMaLsXpuZmZm9k5Npq8s999zD1KlTe/QcK1em51s/+SQ88UTHcuFCGD0a\ndtttwzJ2LLzrXWn+ditd0e6NWA4kjmexHM/iOJbFcjyL5XgWp5lk2nOmrUcMH94x5SNv3bp01frZ\nZzvKrFlp2d4OS5akhHqnnTYsY8ZsWIYNK+VrmZmZmW3AV6atT1m9GhYtguef37AsXpzKiy+m5dCh\nKanefnvYdtv0Ypptt91wfautUhk9Ot1YaWZmZtYVT/OwASECli7tSK5feSU90q96uWRJKq++Cptu\nmpLqSnK95ZawxRYdy8r65punstlmG5ZWmnpiZmZmzelzybSkacAlQBtwRURcULV9KPDfwCRgCXB8\nRDxX4zhOpgs00OZWRaQ3PuaT66VLU3nttY71SnnjjVRefz0tly1L00pGjYKRI9OyUlauvIexY6cy\nYgRvl5Ej2eDz8OG1y7Bh6Qq7HyvYYaCNzZ7meBbHsSyW41ksx7M4fWrOtKQ2YAbwQWAR8ICkWyPi\niVyzk4AlETFO0vHABcAneqpPlixYsGBA/aWTOq4w77pr4/uvX59eYrNsGSxfnpaV9euuW8C++05l\nxYqObYsX8/bnN99MN2NWlxUr0rY1a9JV82HDNiybbtp5GTq0Y1m9PmRI58taZZNN3rksM7kfaGOz\npzmexXEsi+V4FsvxLFdP3oC4L/CniGgHkDQTOAbIJ9NHA+dl678gJd/Ww5YuXVp2F1rKoEEdV6Kr\nzZu3lJNPbv7Y69bBW2+lxDpf3nqrdlm1qmNZWX/jjY7Pq1e/c3316tpl1aqUzK9Zkz5X1tvaOhLr\nrsrgwe9cHzy4o1R/7qy0tXUs77prKREdn/Pbqte7qsuXQYNqf87XV9dVb6uUtrbW+k2C/64Xx7Es\nluNZLMezXD2ZTO8APJ/7/AKwX2dtImKtpNcljY6IV3uwX2Z9Rltbx3SQviCiI6muVVavhrVr03r1\ncs2a9MNBvr6yXqmvLpVtleMuXw7PPZfW163rKPnP1ds6q1u3Lv1WodbnWvX5uny7/Lb161Oc8sl1\nPtmuLp1tl+qvz9fVWu9q+fTT8NBDtbdX11Wv16rrrm1n+9W7T9EFGm/T2efHHoMbbqi/fVd9yNfV\nWi+6Lv8DYK22XR2nnmU9x6tu++qr8MwznR+nnv5214d61rvbp5FjNNqHIttV/h3tiWNb93oymfYk\n5z6qvb297C70G/0tllLHFJAyfOYz7VxySTnnrlfEOxPt6qS7VptKfeVzfj2/rVZ9Z/tV71O97dJL\n2znxxHfWd3aM6vXuPldv62q/WuftyVL5s2qkTVef585tb6h9Z33I19VaL7ouf7tRrbZdHaeeZT3H\nq9X2pZfa+fWva2+rp7/d9aGe9e72aeQYjfah6HZr17ZzwQWN7dPdtnrVm6h3VV9vEt9M4l9vHzZG\nj92AKOn9wPSImJZ9PgdYn78JUdJvszZzJA0GFkfENjWO1TOdNDMzMzPL6TM3IALzgHGSdgH+AhwP\nnFDV5lbgn4A5wHHAnbUO1OiXMjMzMzPrDT2WTGdzoE8HZpEejXdlRDwh6XxgXkTcBlwJXC1pIenR\neH6Sh5mZmZm1jJZ4aYuZmZmZWV80qOwOdEXSNElPSloo6V/K7k+rk9Qu6WFJ8yXNLbs/rUbSjyW9\nJOmRXN1oSXdIelrS7ZK2KLOPraSTeE6X9EI2RudnL36ybkjaSdLdkh6T9KikM7J6j88mdBFPj88G\nSdpU0v2SFmSxnJ7V75rVL5Q0U9ImJXe1JXQRz59KejY3Nt9bcldbiqS2LG63ZZ8bGp99NpnOvfRl\nGrAncIKkPcrtVcsLYGpETIyIfcvuTAv6CWk85p0N3BER40lz/s/u9V61rlrxDOB72RidGBG/LaFf\nrWgNcFZE7AW8H/hi9u+lx2dzOounx2eDIuIt4JCImABMAKZJ2o/0kraLImIc8BrpJW7WjS7iGcCX\nc2Pz4VI72nrOBB6Ht59E19D47LPJNLmXvkTEGqDy0hfbOL6Zs0kRcS/pL1Xe0cBV2fpVwEd6tVMt\nrJN4gsdowyLixYhYkK0vJ70cawc8PpvSRTzB47NhEbEyWx0CbEJKWA4BbszqPTYb0Ek8wWOzKZJ2\nBI4CrqAjhg2Nz76cTNd66csOnbS1+gRwu6R5kjbivX2Ws11EvJStvwRsV2Zn+onTJT0k6UpPS2hc\n9gSlicD9eHxutFw852RVHp8NkjRI0gLSGLwdeAZYGhHZa5BYhP9/r1t1PCOiMm3zm9nY/J6kkt4W\n0JIuBr4CrAeQtBUNjs++nEz7zsjiTYmI9wFHkn5teVDZHepPIt3N63G7cS4DdiP9+nIxcFG53Wkt\nkkYCvwDOjIhl+W0en43L4nkjKZ7L8fhsSkSsz6Yl7Eh6E/LuJXeppVXHU9JewDkRsTuwDzAa8H1m\ndZD0YeDliJhPx1Xphq/w9+VkehGwU+7zTqSr09akiFicLV8BbiZNpbGN85Kk7QEkjQFeLrk/LS0i\nXo4M6VduHqN1ym6Q+QVwdUTcklV7fDYpF8+fVeLp8blxIuJ14G5gf2ALSZUcZEfS//nWgFw8p0XE\ni1ndatL9KB6b9TkAOFrSn4FrgUOBS2hwfPblZPrtl75kv644nvSSF2uCpOGSRmXrI4AjgEe63svq\nUHnxENnyli7aWjeyhK/iWDxG6yJJpOf2Px4R+Reye3w2obN4enw2TtLWlekwkoYBh5PmoN8NfCxr\n5rFZp87iWRmb2dj12KxTRJwbETtFxK6kd53cFRGfosHx2aefMy3pSNJPCJWXvny75C61LEm7kq5G\nQ3pZz88dz8ZIuhY4GNiaNFft34FfAtcDOwPtwMcjYmlZfWwlNeJ5HjCV9Cv0AP4MnJKb82udkHQg\n8DvgYTqmcpwDzMXjs2GdxPNc0lt8PT4bIGlv0g1cbaQLeNdFxDey/5NmkqYkPAh8KnvYgHWhi3je\nCWxDmqIwHzg1d6Oi1UHSwcCXIuLoRsdnn06mzczMzMz6sr48zcPMzMzMrE9zMm1mZmZm1iQn02Zm\nZmZmTXIybWZmZmbWJCfTZmZmZmZNcjJtZmZmZtYkJ9NmZk2QdIykPXKfz5d0WBftp0q6rZtjTpX0\nkyL72duyF229KenBXN12kq6R9IykeZL+IOkj3RznGUnjq+oukfRVSQdKelySX0xhZqVzMm1m1iBJ\ng0lvGduzUhcR50XEnRt56NIf/J99t431p4iYlB1PpLeH3RMRYyNiMulNYzt2c4yZWbtKvwYBHwWu\njYjfA0cW0E8zs43mZNrMBpzs6umTkn6WXeG8IXs1L5L+TdJcSY9I+mFun3skXSzpAeCrwN8DF0p6\nUNJukn4q6aNZ230kzZa0QNL9kkZWnX+EpB9n2x6UdHS2aTWwNGtzsKT5WXmw+hhZm09lx5gv6QdZ\nwomk5ZK+kZ3/PknbZvXbSLox+35zJR2Q1U+XdLWk3wNXZa8svkPSo5Iul9Quaavs6vuZufN/U9IZ\n3YT7UGBVRPyoUhER/xcRM7JjtEm6MOvPQ5K+kDW7Fjg+d5wPAM9FxPOV03dzXjOzXuFk2swGqvHA\nf0bEnsAbwGlZ/YyI2Dci9gaGSfpwVh/AJhGxT0R8C7gV+HJETIqIZ7PtIWkI6arqGRExATgMeLPq\n3F8D7oyI/UjJ5oWShkfEfRFxVtbmS8BpETEROLD6GNkUk48DB2Rt1gOfzDYPB+7Lzv874OSs/lLg\n4ojYFzgOuCJ3yN2BwyLik8B04H8i4m+AG0mvIw/gx8Cns/MPIiW7V3cT571Ir+PtzEnA0qxP+wIn\nS3pXRDwKrJf03qzdJ4BrujmXmVmvczJtZgPV8xFxX7b+M1LCCnCopDmSHiYlunvm9rmu6hjVV0cF\nvAdYHBF/BIiI5RGxrqrdEcDZkuYDdwNDgZ2q2swGLpb0z8CWNY5xGPA+YF52nEOBXbNtqyPi19n6\nH4FdsvUPAjOy9r8ERkkaQUqUb42IVVm7KaQfCIiIWcBr2fpzwBJJE7Lv8GBEvEbXNpi6ImlGdsV8\nbi4Wn876NAcYDYzLtl0LfEJSG3AMcEM35zIz63VFzI0zM2tF+SRPpKvKQ4H/AiZFxCJJ5wGb5tqt\n6OIYXdXV8g8RsbDTzkVcIOlXwN8BsyV9KCKeqmp2VUScW2P3Nbn19XT8Wy9gv4hYnW+cpjWzsuoY\nnU2juAL4LLAd6Up1dx4jzXUGICJOl7QVMC/X5vSIuKPGvjOB24H/BR6OiFfqOJ+ZWa/ylWkzG6h2\nlvT+bP0fgXtJiXOQrr6OBD5WtU8+wVwGbFa1PYCngDGSJgNIGpVdWc2bBbw911jSxOrOSRobEY9F\nxH8AD5CueOfdCRwnaZus/WhJO3f1hUmJaf68f9tJu9mkKSRIOgLYMrftZmAaMDn7Hl2KiLuATSWd\nmqsekVufBZxWufFR0nhJw7N9nwX+CnwHT/Ewsz7KybSZDVRPAV+U9DiwOXBZRLwOXA48CvwWuL9q\nn/xV55nAVyT9UdJubzeIWEOaS/x9SQtIyWIlSa/s/3VgE0kPS3oUOL9G/87MboJ8iHRj4m826EjE\nE8C/ArdnbW4Htq/Rz/x5zwAmZzf6PQac0sl3Ox84Inv03HHAi6QfHirf7y7g+oio9yr8R4CDJT0r\n6X7gp6SbOCFd6X4ceDA732Vs+FvTa0k/SNxU57nMzHqV6v+30Mysf5C0C3BbdpOhVcluolwXEesk\n7U+6UbPyqLtBpHnYx0XEMzX23YVeiK3/DM2sr/CVaTMbqHwloXM7Aw9kV9YvJXsaiKQ9gYWkJ328\nI5HOrAU2V+6lLUWTdBDpaSqeQ21mpfOVaTMzMzOzJvnKtJmZmZlZk5xMm5mZmZk1ycm0mZmZmVmT\nnEybmZmZmTXJybSZmZmZWZOcTJuZmZmZNen/AZygtJfgoEnpAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x10835fd50>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot time of flight as function of particle's energy\n",
    "\n",
    "# assume it's a charged pion\n",
    "pionMass = 139.57018e-3 # in GeV\n",
    "\n",
    "energies = pylab.linspace(1,40, 100)\n",
    "\n",
    "pylab.figure(figsize = (12,12))\n",
    "\n",
    "pylab.plot(\n",
    "  energies,\n",
    "    [ 1e9 * (# convert to nanoseconds\n",
    "     tof.evalf(subs = {\n",
    "                p : math.sqrt(energy**2 - pionMass**2),\n",
    "                B : 3.8, # in Tesla\n",
    "                mP : pionMass,\n",
    "                phi : scipy.optimize.brentq(func, \n",
    "                            0, math.pi, # bracketing range\n",
    "                            args = (3.8, \n",
    "                                   math.sqrt(energy**2 - pionMass**2)\n",
    "                                   ))\n",
    "                \n",
    "            })\n",
    "     # subtract zero magnetic field time of flight\n",
    "     - rEcal / c\n",
    "     )\n",
    "     \n",
    "     for energy in energies]\n",
    ")\n",
    "\n",
    "pylab.xlabel(\"particle's energy [GeV]\")\n",
    "pylab.ylabel(\"time of flight difference w.r.t zero magnetic field [ns]\")\n",
    "\n",
    "pylab.grid()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}