restrepo/lfv.ipynb

## lfv.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  Inert doublet model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Populating the interactive namespace from numpy and matplotlib\n"
     ]
    }
   ],
   "source": [
    "%pylab inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import os, sys, inspect\n",
    "import commands\n",
    "from hep import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define functions to change from general basis to physical basis "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Check $\\mu \\to e \\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "SM={'G_F': 1.166e-05, 'alpha_em': 0.0078125, 'm_e': 0.000511, \\\n",
    "                'm_mu': 0.1057, 'm_tau': 1.777, 'vev': 246.0, \\\n",
    "                'MW': 80.385, 'MZ': 91.1876, 'WH':4.21e-3,'m_pi0':134.9766E-3,\\\n",
    "                'm_pipm':139.57018E-3,'hbarc':1.973269631E-16,'V_ud':0.97425,\\\n",
    "                'f_pi':131E-3}\n",
    "def Fme(x,xmin=0.996,xmax=1.005,xfit=1.001):\n",
    "    \"\"\"Fixing near to one values\n",
    "     xmin: close to 1 from below\n",
    "     xmax: close to 1 from above\n",
    "     xfit:  optimized 1 limit\n",
    "    \"\"\"\n",
    "    x=np.asarray(x)\n",
    "    if x.shape:\n",
    "        x[np.logical_and(x>xmin,x<xmax)]=xfit\n",
    "    else:\n",
    "        if x>xmin and x<xmax:\n",
    "            x=xfit\n",
    "        \n",
    "    return (1.-6.*x+3.*x**2+2.*x**3-6*x**2*np.log(x))/(6.*(x-1.)**4)\n",
    "\n",
    "def LFV(SM,dp,yuk):\n",
    "    '''Oscar Notes with \\pi^2 -> \\pi\n",
    "    '''\n",
    "    import numpy as np\n",
    "    const=False\n",
    "    SM=pd.Series(SM)\n",
    "    y=np.matrix(yuk)\n",
    "    \n",
    "    dp=pd.Series(dp)\n",
    "    FMl=[]\n",
    "    for i in range(1,4):\n",
    "        dp['Mtrp%d' %i]=dp['Mtr0%d' %i] #degenerate fermions\n",
    "        FMl.append(1./(dp.MHC**2)*Fme(dp['Mtr0%d' %i]**2/dp.MHC**2)\\\n",
    "                   -1./(dp['Mtr0%d' %i]**2)*( Fme(dp.MH0**2/dp['Mtrp%d' %i]**2)\\\n",
    "                                            + Fme(dp.MA0**2/dp['Mtrp%d' %i]**2)   )  )\n",
    "    FM=np.matrix(np.diag(FMl))\n",
    "    \n",
    "    Brmueg  =(3.*SM.alpha_em/(4.*16.*np.pi*SM.G_F**2)*np.abs(y[:,0].T*FM*y[:,1].conjugate())**2)[0,0]\n",
    "    Brtaumug=(3.*SM.alpha_em/(4.*16.*np.pi*SM.G_F**2)*np.abs(y[:,1].T*FM*y[:,2].conjugate())**2)[0,0]\n",
    "    if (Brmueg<5.7e-13):\n",
    "       if (Brtaumug<4.5e-8):\n",
    "           const=True\n",
    "\n",
    "    return Brmueg,Brtaumug,const"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare Implementation: \n",
    "\n",
    "Obtain the required info for masses and Yukawas from the SPheno spc output file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def lfvin(spc):\n",
    "    dp=pd.Series()\n",
    "    y=np.zeros((3,3))\n",
    "    for i in range(1,4):\n",
    "        dp['Mtr0%d' %i]=eval(spc.blocks['MASS'].entries[6000056+i])\n",
    "\n",
    "    for i in range(3): \n",
    "        for j in range(3): \n",
    "            y[i,j]=eval(spc.blocks['YN'].entries[i+1,j+1])\n",
    "\n",
    "\n",
    "    dp['MH0']=eval(spc.blocks['MASS'].entries[35])\n",
    "    dp['MA0']=eval(spc.blocks['MASS'].entries[36])\n",
    "    dp['MHC']=eval(spc.blocks['MASS'].entries[37])\n",
    "    return dp,y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Edit the input file: \n",
    "\n",
    "Here input file is modified with the object structures in hep.py module. But similar steps can be done by editing the file by hand"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "''"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "commands.getoutput('cp LesHouches.in.FT_SD_3gen LesHouches.in.radinertRSIII')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<PySLHA Doc: 33 blocks, 24 decays, 0 xsections>"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a=hep(MODEL='radinertRSIII',spc_input_file='LesHouches.in.radinertRSIII')\n",
    "a.runSPheno() #creates SPheno.spc.radinertRSIII"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We now obtain the required data for the $\\mu \\to e \\gamma$ analytical calculations, and get the corresponding evaluation. See the definitions of the functions `lfvin` and `LFV` before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.8575984851104539e-25, 3.0751963715623622e-25, True)"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dp,y=lfvin(a.LHA_out_with_comments)\n",
    "LFV(SM,dp,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We check also for the result with the yukawas transposed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "LFV(SM,dp,y.transpose())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we compere with the numerical result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'4.32974619E-26  # BR(mu->e gamma)'"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.LHA_out_with_comments.blocks['FLAVORKITLFV'].entries[701]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we change some input value of the input file, and check the change"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4.0261471288342711e-16, 2.0913790619140509e-15, True)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a.LHA.blocks['YNIN'][1,2]='-1E-0  # Yn(1,2)'\n",
    "a.runSPheno()\n",
    "dp,y=lfvin(a.LHA_out_with_comments)\n",
    "LFV(SM,dp,y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scan upont the mass squared scalar potential parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "commands.getoutput('cp LesHouches.in.FT_SD_3gen LesHouches.in.radinertRSIII')\n",
    "a=hep(MODEL='radinertRSIII',spc_input_file='LesHouches.in.radinertRSIII')\n",
    "df=pd.DataFrame()\n",
    "a.LHA.blocks['SPHENOINPUT'].entries[55]='0               # Calculate one loop masses'\n",
    "dm_masses=np.logspace(np.log10(2E3),np.log10(1E6),100)\n",
    "for MHX in dm_masses:\n",
    "    if np.where(dm_masses==MHX)[0][0]%10==0: #find the index of the array entry\n",
    "        print np.where(dm_masses==MHX)[0][0]\n",
    "    a.LHA.blocks['MINPAR'][6]='%0.8E  # mEt2' %MHX\n",
    "    a.runSPheno()\n",
    "    dp,y=lfvin(a.LHA_out_with_comments)\n",
    "    a.Series=a.Series.append(pd.Series({'muegamma_analy':LFV(SM,dp,y)[0],'muegamma_analyT':LFV(SM,dp,y.transpose())[0],\\\n",
    "                                         'muegamma':eval(a.LHA_out_with_comments.blocks['FLAVORKITLFV'].entries[701])}))\n",
    "    df=df.append(a.Series,ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdd0119aed0>]"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xd4lUXax/HvgBQpoUhNkCZSpQhSVCChimBAmjRR1BVE\nhVV0XV0RgqhYQUUsoKDiIijVIAq8QGDp4NKWLiKBgFJCqEIgmfePOSEJBEg/J8nvc11znZwnc57M\n4SHnzjwz94yx1iIiIpLL2w0QERHfoIAgIiKAAoKIiHgoIIiICKCAICIiHgoIIiICKCCIiIhHsgKC\nMeYLY8yfxpjNlx1vZ4zZYYzZZYz5ZxKv62SMGW+M+dYY08ZzLNAYs8wY84kxpnn6vA0REUmr5PYQ\nJgH3JDxgjMkFfOQ5XgvoZYypnrCOtXaOtbY/MBB4IO4wcArIBxxIfdNFRCQ9JSsgWGuXA8cvO9wI\n2G2t3WetvQBMBTpd5RRDgXGecy2z1nYAXgReTVWrRUQk3aVlDCEA2J/g+QHPMYwxfY0xo40x/saY\nN4F51tqNl70+Csibhp8vIiLp6IaMOKm1djIw2RgzCGgF+BljqlhrxxtjOuNuMxXB3XISEREfkJaA\nEAGUT/C8nOfYJdbascDYy47NAmZd68TGGK24JyKSCtZak9rXpuSWkfGUOOuAKsaYCsaYvEBP4IfU\nNuRy1lqvl+HDh3v9XCl5XXLqXq/O1b6fkuPp+e/mC9cuq1y/lH7PV69dVrx+vvK7l1bJnXY6BVgJ\nVDXGhBtjHrHWxgCDgAXAVmCqtXZ7mlvkQ4KCgrx+rpS8Ljl1r1fnat9P6XFvS+92ZYXrl9Lv+eq1\ng6x3/bLL755Jj6iS3owx1hfbJckTEhJCSEiIt5shqaBrl7UZY7CZdMsoU4WEhBAWFubtZkgq+PJf\nnnJtunZZU1hYWLoEcvUQRESyiWzbQxARkcylgCAiIoAPBwSNIYiIJI/GEEREJBGNIYiISLpQQBAR\nEcCHA4LGEEREkkdjCCIikkhaxxAyZPnr7OD112HTJihcGAoVAj8/97Wf35WlSBFX/PzgBv2LikgW\npR7CVaxdC3v3wqlTicvJk+7xxAlXTp50Je7rAgWgaNHEpVgxKF48/jFhuekmV4oUgVw+ewNPRLKC\ntPYQFBDSUWwsnD4Nx4+7ABEV5b5OWCIj48uxY/HlzBkXGEqUcKVkSShVKv4xrpQu7UrRomBSfdlF\nJDvKtreMQkJCCAoKylKLbeXKFX8bKaWio12QOHLElaNH3ePhw7B1KyxZAn/+GV/OnYMyZeKLv/+V\nJSDABRkFDpHsLSwsLF0m4aiHkEX99ZcLDIcOJS4HD7oSEeHKX3+5wFCuHNx8c3ypUAHKl3ePhQt7\n+92ISHrQLSO5prNnXWDYvz++hIfHl337IH9+FxgqVoRKleJL5cruMX9+b78LEUkOBQRJE2vdGMbv\nv7tB9Ljy22+uhIe7sYsqVeDWW91j1aqu3HIL5Mvn7XcgInEUECRDXbzoehV79sCvv8Lu3a7s3Ol6\nF/7+UL06VKvmHmvWhBo13MC4iGQuBYSMcviweyxUCG68USOzSbhwwfUmduyIL9u3u5I3rwsOtWrB\nbbdB7drusWhRb7daJPvKtgFh+PDh3p1l9PjjMGeOm0caHQ0FC7rR14QlqSy1okXjM9UuT0bIITfj\nrXUD3Fu3uvK//8GWLe7r4sWhbt34Uq+eu/WkHAyR1IubZTRixIjsGRB8ql0XL7pEgcsz1E6cSPx1\n3GNUVHzmWlwyQlSU62UklaV2002Js9QSJiPcdBPkyePtf4F0ERvrehSbNsWXDRvcP0+dOlC/PjRo\n4Er16sr6FkmpbNtD8Ha7nnzySVasWEGhQoUulcKFCycqfn5+l0qRIkUuPcaVPAk/yK11yQORkVdm\nqSXMUDt2zCUgHDvmkhEiI91tq7gstdKl47PUypRxzxMmJBQo4L1/tFSKjISNG+G//3Vl/Xo3dbZe\nPWjYEBo1gsaN3Ywn3bkTuToFhAxy8OBBjhw5wunTpzl16tSlx7hy8uRJTp48eenrEydOXHqMiori\n5MmT5MuXj6JFi1K0aFGKFSt26bF48eKXHhOWEiVKcNNNN1GkSBFyxd1DiY11wSMuSy2uxGWo/fGH\nK3FJCXnzxmemlS3rkhDiSrlyrpQtC7lze/Xf93qiolxwWLsW1q2DNWvg/Hlo0iS+NG7sYqWIOAoI\nPspay5kzZ4iKiuL48eNERUURFRVFZGQkx48f5/jx40RGRhIZGcmxY8cSlbNnz14KECVKlKBkyZKU\nKlWKkiVLUrp0aUqVKkWpUqUoU6YMpUuXxs/PD2OM64VERSXOUEuYpRaXkHD0qOtZlC/vstTiMtTi\nkhEqVvTJT9qICFi92pWVK12volo1uPtuaNrUlYAAb7dSxHsUELKhCxcucOzYMY4cOcKRI0c4evQo\nhw8f5siRIxw+fJjDhw/z559/XioXLlygTJkyl4q/vz/+/v6ULVuWgIAA/P39CQgIoFixYi5wXLjg\nAsXlGWpx5fff3a2nuOy0uHLLLa6UK+cTPYzz510vYsUKWL7cFT8/aN4cgoIgMNC9BZGcQgFBOHv2\nLH/88QeHDh1KVA4ePEhERAQHDx7kwIEDXLhwgXLlynHzzTdTrlw5ypcvz80330z58uWpUKEC5cuX\np2DBgq6ncfhw4iy1PXviy9Gj7pP21lvjs9SqVnV/rpcp47Ub/bGxburr0qXxJV8+FxxatoRWrdSD\nkOwt2wYEr087zYZOnTpFREQE+/fv58CBA4SHhxMeHs7+/fvZt28f4eHhFC5cmIoVK1KxYkUqVapE\npUqVqFy5MpUrV6ZChQpuoPzsWRcY4rLUdu1ymWq7drk/26tXd6VGDZeMULOmCyCZ3Kuw1jVryRJX\nFi92k7dat4Y2baBFi9QtRCjiazTtVNKdtZY///yT33//nb17914qv/32G7/99hsHDx4kICCAKlWq\nXCpVq1alatWqVKpUyQWLY8fis9S2bYsvR4+6IHHbba7UqeOy1cqWzbQeRWysm+q6cKErq1e7XIh2\n7VypX1/5EJK1Zdsegi+2K6eLjo5m3759/Prrr+zevTtRiYiIoEKFClSvXv1SqVGjBjVq1KBIkSIu\nX2PbtvgstS1b3KezMYmz1G6/PdOSEP76C5Ytg/nz4eefXcxq1w7at4d77nGpIiJZiQKC+ITz58+z\nZ88eduzYwfbt2xM9Fi1alFq1alGrVi1uu+02brvtNmrVqkXBAgXcjKiEWWobNrjpRLVrx2eqNWzo\nbjtlcJDYtw/mzXNl6VIXm+67D4KDXYwS8XUKCOLTYmNjCQ8PZ+vWrWzdupX//e9/bNmyhZ07d+Lv\n70/dunUvlXr16lG+fHnMqVMuQPzyi8tSW78eDhxwvYjGjeMz1SpWzLDbTX/95cYcQkNh7lw3C7dT\nJ7j/fvejdWtJfJECgmRJFy9eZPfu3WzatOlS2bBhA9HR0dSrV4/69evToEED6tevT5UqVch16pQL\nEGvWuGy11avdqHGTJnDXXa7ccUeGrBdlrfvRs2e7EhnpAkPnzm4GUzZZWUSyAQUEyVb++OMPNmzY\nwIYNG/jll19Yv349J06coEGDBjRs2JDGjRvTqFEjAvz9Xf7EqlWurFjhllmtW9dlqDVr5jLWihdP\n9zbu2gWzZsGMGW5G7v33Q7dublqrgoN4kwKCZHtHjhxh/fr1rFu3jrVr17JmzRry5s1LkyZNuPPO\nO7nzzjtp0KAB+WNiXA9i+XL4z3/c15UquQy1wED35/xNN6Vr2/btc4Hh++/dfhGdO0PPnu7H+UDu\nnuQwCgiS41hr2bt3L6tXr2bVqlWsWrWK7du3U6dOHe6++26aNWvG3XffTYkiRdy9nqVLISzM9SIq\nV47PUmvePF03lN63D777DqZNc4ngPXpA797uTpYW5ZPMkG0DghLTJCXOnDnDunXrWL58OcuXL2fV\nqlUEBAQQGBhI8+bNCQwMxL9kSTdAvXgxLFrkxiJuv91lqbVt6z6502km086d8O238O9/uwHoBx90\nRUtpSEZQYprINVy8eJFNmzaxbNmyS6VEiRK0aNGCFi1a0LJlS0oWLOhuLS1cCAsWuJlMrVvDvfe6\nRAR//zS3w1oXdyZPdj2HGjXg4Yehe3dlSUv6y7Y9BF9sl2RdsbGxbN68mSVLlrB48WKWLVtGpUqV\naN26NW3atKFZs2YUiIqKz1JbuNBNa23fHjp0cFNd0zgoEB0NP/0EX33lOinBwfDYY+7OlaaxSnpQ\nQBBJhQsXLrB+/XoWLlzIwoUL2bhxI40bN6Zdu3a0a9eOWtWqYVavdllqc+e6xf7at4eOHd3tpYIF\n0/Tzjxxxt5O++MItDfXoo/DII+nSKZEcTAFBJB2cPHmSJUuWMH/+fH766ScuXrzIvffeS4cOHWjV\nqhWFjh51WWo//OBmLwUFufmmwcFuJ7tUstZtAPT5526mUmAgDBjgYo5mKUlKKSCIpDNrLbt27WLe\nvHnMmzePNWvWcNdddxEcHExwcDDl/fxcz2HWLDf20KABdOniShr+xD91CqZOhU8/dclv/fu7W0ql\nSqXjm5NsTQFBJIOdOnWKBQsWEBoayo8//ki5cuXo1KkT999/P3WrVsUsXOiSEebOdSu5du/uMtXK\nlk31z1y/Hj75BGbOdEMYTz/tlszQ9FW5FgUEkUwUExPDypUrmTNnDrNmzSI2NpYuXbrQtWtXmtx+\nO7kWLXLJCKGhbkprz57QtWuqE+KOH4dJk2DcOLf66uDBLr8hX750fmOSLSggiHiJtZYtW7Ywc+ZM\nZsyYQWRkJF27dqV79+7c3aABuebPd/eAfv7ZTSXq08cNShcokOKfFRvrZih98AFs3uzGGZ580m2N\nLRJHAUHER+zYsYPp06fz/fffc/ToUR544AF69OhB45o1MbNnu2lFa9e6weiHHnIjyKmYb7ptG3z4\noctr6NwZhgxxd6pEFBBEfND27duZNm0aU6dOJTo6ml69etGnTx9qFisGU6bA119DVJTLUuvXzy2p\nkUJHj7oB6HHj3N5CL7zgJj9pnCHnUkAQ8WHWWjZu3MiUKVP49ttvKVmyJH379qV3796U+eMPN0Aw\nZQrUquWmFHXtmuJbSufOwTffwLvvun0bXnzR9Rw0bTXnUUAQySJiYmJYunQpkydPZvbs2dx55530\n69ePju3akX/hQpeMsGaNG4geMMAt5Z0CsbEuTeLNN9201RdecHem8ubNoDckPietAcFnE+ZDQkII\nCwvzdjNE0k3u3Llp2bIlkyZNIiIigj59+jBhwgQCKlXiyUWL+OW117AbNriR4vvuc/NMJ050qczJ\nkCuXG55YtQomTHCJbrfcAmPHuh3gJPsKCwsjJCQkzedRD0HEy8LDw/nqq6+YOHEiRYoU4bHHHuPB\nXr0otno1fPaZ+4Tv2xcGDoSqVVN07nXr4PXXXcfj+efhiSfSvOqG+LBs20MQySnKly/PK6+8wp49\ne3j33XdZsWIFlapUod/06ax86SXs2rUu8aBpU2jXziXAxcYm69wNG7ptP+fPd3HlllvgnXfgzJkM\nflOSJamHIOKDjhw5wldffcX48eO58cYbGThwIH26dqXwTz+5e0CRkS5L7ZFHUrSO9v/+ByNGuE3l\nXnjB9RhuvDED34hkKvUQRLKhkiVL8vzzz7Njxw7ee+89FixYQIVq1Ri8fj07J09204pWrXI77jz7\nLPz+e7LOe9ttbmxh/nxYtgyqVHFLZERHZ+z7kaxBAUHEh+XKlYvWrVszc+ZMNm3ahJ+fH80DA2k7\nfDjzHnqI2P/+100jatDArWmxbl2yzlunjlubb84cV6pXd5v4xMRk8BsSn6ZbRiJZzPnz55k2bRof\nfPABp0+fZvDgwfTr2pWCU6fCmDEuye2FF9x4QzKz1JYtg5degpMnYdQot6CeEtyyHuUhiORQ1lpW\nrFjBmDFjWLp0KY8//jhPDxhAwIoV8PbbbrOFF1+EBx5I1l7R1ro1+V56CYoXd4PPTZpkwhuRdKMx\nBJEcyhhD06ZNmTFjBmvWrOHMmTPUrl+fRxctYtuUKfDWW25ti2rVYPx4OH/+Oudza+9t3uzGqrt1\nc7Fkz55MekPidQoIItnALbfcwocffsju3bupXLkyLVu1Ivjjj1n55ptuE+dZs9wI8kcfXTdLLXdu\nt6Xnzp1urKFRI3juObf0kmRvCggi2chNN93E0KFD2bt3L+3bt+fBBx8kaOhQFjz7LHbGDFi40CUj\nfPDBdQNDwYIwdChs3erGFqpVcwvpXbyYSW9GMp3GEESysYsXLzJt2jRGjRpFgQIFGDp0KMHlymFe\nfdXNSPrnP91enfnzX/dcmza5pbYPH3bxpGXLTHgDkiIaVBaR64qNjWX27NmMHDkSay3Dhg3j/vLl\nyTVihPukHzrUDRzkyXPN81jrMp+fe84tuT16NFSsmDnvQa5PAUFEks1ay9y5cwkJCSEmJobhw4dz\nf5kymGHDYO9eePVVt9rqdTbuOXfOLbf9/vsuYfqFF5LVyZAMpoAgIikWFxiGDx+OMYZXX32V9jfe\niPnXv9yn/RtvwL33XjcZYd8+11vYsMHt4tahQya9AUmSAoKIpJq1llmzZjFs2DAKFy7MqDfeIOjE\nCZeMUKaMy2do2PC651mwAJ56CmrXduMLN9+cCY2XKygPQURSzRhDly5d2LRpE08//TSP/e1v3PPJ\nJ/zy5ZfQu7fbYKFnT3c76RratoUtW9y4wu23w3vvaTZSVqSAICLkzp2bPn36sH37djp16kRw5870\nWryYvfPnQ82acMcdbkbSiRNXPUf+/DBsGKxeDT/95DoWyVxaSXxEsgKCMeYLY8yfxpjNlx1vZ4zZ\nYYzZZYz5ZxKv62SMGW+M+dYY08ZzzBhjXjPGfGiM6Zs+b0NE0kPevHl58skn2bVrFzVq1OCOwECe\niYwkctkyOHbMJSN8+uk1V8GrUsWlOzz3HAQHu8VYtf9C1pDcHsIk4J6EB4wxuYCPPMdrAb2MMdUT\n1rHWzrHW9gcGAg94DncCygHRwIHUN11EMkqhQoUYNmwY27Zt48KFC1QLCmJ0zZpE//ADTJ3q7gst\nXnzV1xsDDz7o9l+IjHTLbi9YkIlvQFIl2YPKxpgKQKi1to7neRNguLX2Xs/zFwFrrX0ride+C3xj\nrd3o6UlEWmsnGGO+t9Z2T6K+BpVFfMj27dv5xz/+wfbt23nn7bfpbC3mH/9wt5Leew/Kl7/m6+fP\nd5vxBAW53IVixTKn3TmNNweVA4D9CZ4f8BzDGNPXGDPaGONvjHkTmGet3eiptx847vlaq6+LZAE1\natRg7ty5fPbZZwwPCaHluHFsnjoVatVyvYXXXnPTVa/innvcoHPBgm4mUmhoJjZeki1DBpWttZOt\ntUOArkAroJsxpr/n27OAdsaYD4ClGfHzRSRjtG7dmg0bNtC9e3fadOzIk4cPc3zhQli/3q2Ed437\nQoUKubX1/v1vN67Qty8cP37V6uIFab1lFGKtbed5ftVbRilulDF2+PDhl54HBQURFBSU1tOKSDqK\njIzklVdeYfr06YwYMYLHy5Yl97PPuttIY8ZAQMBVX3vmjNuqYfZsmDDB7eUjKRcWFkZYWNil5yNG\njMicxDRjTEVcQKjteZ4b2InrARwC1gK9rLXbU9uYBD9LYwgiWURcDsO5c+f4ZPRo7liwwM1EGj4c\nBg5062lfxaJF8NhjLiC89567pSSplyljCMaYKcBKoKoxJtwY84i1NgYYBCwAtgJT0yMYiEjWUrdu\nXZYtW8agQYMIfuABBh49yonQUPjuO7jzTrd43lW0auW+/ddfbihCeQvelayAYK3tba31t9bms9aW\nt9ZO8hz/yVpbzVp7q7X2zfRsWEhISKKukIj4LmMMDz30ENu2bcMYQ40uXZjSvz+2f39o0wZefvmq\ng85Firg9fF57za2FNHKkspxTKiwsjJCQkDSfR2sZiUi6W7NmDf3796d06dKMHzGCiqNHu705P/8c\nmjW76usiIqBfPzfGMHmy28tHkk9rGYmIz2ncuDHr16+nTZs23BEczLuNGxPzxhtuXaRBg+D06SRf\nFxDgchYeeACaNIEvv3R7MEjmUA9BRDLUnj17GDBgAMePH+fL0aOpPXEiLFsGEydCixZXfd2WLS5+\n1KsHn3wCfn6Z2OgsKtv2EDSGIJI93HLLLSxcuJBBgwbRqnt3hlWowIUxY1wiwuDBV13oqHZtN8hc\nuDDUr68B52vRGIKIZDkHDx5k4MCB7Nmzh8kffMDtkybBmjXw9dduRtJVTJ8OTz7pxqYHD77uvj05\nljbIEZEsxVrL1KlTeeaZZ+jfvz/Datcmz6BBLmdh6FC44YYkX/fbb9CjB5Qr5+42aT2kK2XbW0Yi\nkj0ZY+jVqxcbN25k06ZN3PH662z95htYudLNQNqzJ8nXVa4My5dDhQrQoIFbLUPSl88GBI0hiGRv\nZcuWZc6cOQwZMoSg3r15q0ULYrt3d9OLvvgiyelF+fLB++/DO+9A+/YuIVo3EzSGICLZyL59+3j4\n4YeJiYnh25dfptyLL0KlSjB+PJQsmeRrdu2Cbt2gbl347DMoUCCTG+2DdMtIRLK8ChUqsHjxYjp1\n6sTtffsy5e9/h1tvveZGPFWruu06jXHj0Ve50yQpoB6CiPiUjRs30rt3b+rVq8f47t0p9NRT8Oij\nEBKS5ICztfDxxzBiBEya5Ja/yKnUQxCRbKVevXqsX7+eokWLUnvIENZPmOBGkIOC4MCVu+4aA089\n5ZbSHjDArYkUG5v57c4OfDYgaFBZJOcqUKAAH3/8MWPGjKHDo48yqlkzYjt0cHst/PRTkq+56y5Y\nuxZ+/NGNLZw6lcmN9iINKotIjrB//3769OlD/vz5mfb00xR76il48EG3LGoSt5DOn3fLJa1cCT/8\n4Kar5hS6ZSQi2drNN9/M4sWLadSoEbWffJIVY8e6W0ht28Kff15RP18+N+to4EDXa1iyxAuNzqLU\nQxCRLGP+/Pn069ePpwcO5F/R0Zgvv4Rp0+Duu5Osv2gR9O7txqMHDszUpnqFlq4QkRwlIiKCHj16\nUKRIEaY+9BCFBw+GV15xI8tJLHL0668QHOz26Rk9+qorY2QLumUkIjlKQEAAS5YsoXr16tR58UW2\nfPYZTJgADz8MZ89eUb9KFVi1CnbscIHhxAkvNDqL8NmAoFlGInI1efLk4b333uPdd9+lVf/+fPXE\nE26uadOmEB5+Rf2iRWHePDfAfNddsG+fFxqdgTTLSEQE2LFjB507d6Z5s2aMq1KFG8aMceMKzZsn\nWf+DD+Dtt13eQsOGmdzYDKZbRiKSo1WvXp21a9dy9Ngxms6cybExY6B7d5e+nIS//919q317FxQk\nngKCiGR5hQsXZvr06QQHB1P3+efZOG4cfPSR21XnwoUr6nfq5PLbnnoKxo71QoN9lG4ZiUi2Ehoa\nymOPPcbo4cN5cN48l6n23XdQvPgVdX//He69161/9PbbkCuL/4msaaciIpfZvn07HTt2pNN99/G2\nMeSaO9etaXHrrVfUjYx0PQZ/f7eTZ758XmhwOtEYgojIZWrUqMGaNWvYsHkzwTt3cvbpp91ubEuX\nXlG3eHFYuBBiYlxvISdPS/XZgKBppyKSFsWLF+fnn3+mQoUKNPzsMw69+y488AB8+eUVdfPndxOT\natSAwEA4dCjz25sWmnYqIpJMY8eOZdSoUcx77z3qvfwyPPQQDB9+RWaztfDGG24Hz/nzk7zD5NM0\nhiAikgzz5s3j4YcfZvzIkXSeOBFq1nRbdObNe0XdCRNcvPjxR7dpW1ahgCAikkybN28mODiYwY89\nxpD16zFnz8LMmeDnd0XdGTPcgnjTp181x83nKCCIiKRAREQEHTp0oEnDhnycKxe51qxxSQlly15R\nd9Ei6NnTDTtkha05NctIRCQFAgIC+M9//sPv+/dzX3g45zt1cgsc7dx5Rd1WrWDuXLel87RpXmhs\nJlNAEJEcp3DhwoSGhlLW35+m8+Zx8pln3PSitWuvqNu4sZuWOmSIG1vIzhQQRCRHypMnD59//jnt\n2rWj/tixHBo50t0XWrDgirp16kBYGLz+ulscL7vKxltFiIhcmzGGkSNHEhAQQIPhw1ny5ptU69vX\nfer37Jmo7q23ury2Vq3gr7/gxRe91OgMpIAgIjneE088QalSpWj2xBPMHTWKRs8/D1FR8MQTiepV\nqADLlrmgcPYsjBiR5CZtWZbPBoSQkBCCgoIICgrydlNEJAfo0qULxYoV474ePfh6+HDavf22CwqX\ndQX8/V1PoU0biI6GUaO8HxTCwsLSZWUHTTsVEUlgw4YNdOjQgfeee45eEye6fTeT+NQ/etQFhVat\n4J13vB8UQHkIIiLpbteuXbRp04YXH3+cgbNnw513unGFy9bHjoyEtm3dzp1jxng/KCggiIhkgPDw\ncNq0acPD99/PS8uXY6pXd0td5M6dqF5UlAsKd94J77/v3aCggCAikkEOHz5M27ZtubdZM97Ytg1T\nurTbNOGGxMOvcUHhrru821NQQBARyUCRkZG0a9eOJnXr8sH+/ZjChWHKFMiTJ1G9qCg3ptC0KYwe\n7Z2goIAgIpLBTp48SYcOHahesSLjIyMx+fLB1KlXrJR6/LgLCkFB3hlo1lpGIiIZzM/Pj59//pm9\nhw7xUKFCxMbEQNeubr/mBIoVc4nO//d/8PLLbn+FrEQBQUQkGQoWLEhoaChHTpygT548xN5wA3Tr\ndkVQKF7cBYTQUHj1VS81NpUUEEREkunGG29k9uzZnDp3jp5AbO7c0L27y1BLoEQJFxSmToW33vJO\nW1NDAUFEJAXy58/PjBkzOBcT44KCMUkGhdKlXVAYPx4++sg7bU0pBQQRkRTKly8f33//PX/FxNA7\nd25iY2Ohd2+4eDFRvYAAt8nOO+/AxIleamwKKCCIiKRCvnz5mD59Oqejo92Ywpkz0LcvxMQkqlex\nottP4ZVX4LvvvNPW5PLZgBASEpIuizWJiGSUfPnyMWPGDE6cO0e/QoWwR4647dViYxPVq1rV7dI5\naBDMm5eDNdVqAAALdUlEQVT+7QgLCyMkJCTN51EegohIGp07d47g4GAqlirF+AMHMDVrwscfX5GI\nsHo1dOwI06dD8+bp3w7lIYiIeFn+/PmZM2cOuyMiGFyxIvaXX+CFF65IRGjSxCU5d+sG//2vlxp7\nDQoIIiLpoECBAsydO5f//vorL9Wrh50/H0aOvKJe69bw6adw332we7cXGnoNPrtBjohIVlOoUCF+\n/PFHWrVqhV+LFvzr3/8GPz945plE9bp0gWPH3IJ4K1a4TXd8gQKCiEg6Klq0KPPnzycwMJAi99/P\nU2PGQNGi0K9fonqPP+6Cwj33uB3Yihf3TnsT0qCyiEgGOHToEM2aNWNE7970+fxzGDcOOndOVMda\neO45WLvWTU298ca0/Uytdioi4qP27t1L8+bN+WzAANp/+KFby6Jly0R1YmNd+sKZM2720Q1puG+j\ngCAi4sO2bt1Kq1atmPn3v3PXmDHw889Qv36iOtHRbpC5YkX47LPUL5utaaciIj6sVq1a/PDDD9w/\nZgxbBw92n/y//pqoTt68MGOGm4o6YoSXGooCgohIhmvUqBFTpkyhxYcfcuDxx930okOHEtUpXBh+\n/NHt0OmtdY90y0hEJJN89913DBkyhC09elBs8WI3vcjPL1GdnTshMBAmTYJ7703Z+TWGICKShYwb\nN473x4xh8913c+PBg65bcNlWnCtXQqdObrihQYPkn1tjCCIiWchTTz1Fz169aLF1Kxfz5k1yMby7\n7nKDyx07wr59mdc29RBERDKZtZa//e1vRB44wIyTJ8kVFASjRl1Rb8wY+OILWL7c5bZdj3oIIiJZ\njDGGTz/9lOgbbuCZSpWw06e7rdUu88wz0KIFdO16xYZsGdMuX/xLXD0EEckJzpw5Q4sWLejTuDF/\nnz7dTS+6bCQ5JsYlOJco4XoL18pRyJQegjHmC2PMn8aYzZcdb2eM2WGM2WWM+WcSr+tkjBlvjPnW\nGNPGc6ypMeYTY8wEY8zy1DZcRCSrK1iwIKGhoXzw44/MffRRePhh2LAhUZ3cud2S2Rs3wltvZWx7\nktVDMMY0BU4DX1tr63iO5QJ2Aa2Ag8A6oKe1dkcSry8KvGOtfTzBsU5AKWvthCTqq4cgIjnGzp07\nCQwMZMGAAdSZOBFWrYJy5RLViYhw+ymMGeP2U0hKpvQQrLXLgeOXHW4E7LbW7rPWXgCmAp2ucoqh\nwLjLjvUGpqSgrSIi2VK1atX4/vvvaf3JJxzq1g2Cg+H06UR1AgLghx9g4EC3GF5GSMugcgCwP8Hz\nA55jGGP6GmNGG2P8jTFvAvOstRvjKhpjbgairLVn0vDzRUSyjWbNmvHhhx/SZMYMztasCb16uQGE\nBG6/3Y0jdO4M+/df5URpkCGzjKy1k621Q4CuuFtK3Ywx/RNUeQyYlBE/W0Qkq+rZsyeP9+9Pyx07\niDl9Gp5//oo6HTu62UcdO7oVUtNTWjbIiQDKJ3heznPsEmvtWGDs5S+01oZc7+QhIfFVgoKCCAoK\nSl0rRUSykJdffplff/2VfkeP8vW8eZgaNaB//0R1nn8etm+Hdu3CaNkyLNWro14u2dNOjTEVgVBr\nbW3P89zATlwP4BCwFuhlrd2e5kZpUFlEcrDo6Gjatm1L+1tv5YXQUPj2W5eQkKiO25+5WTN4/XV3\nLLOmnU4BVgJVjTHhxphHrLUxwCBgAbAVmJoewUBEJKfLmzcvM2bMYPySJfz88MPQsyfs2nVZHbdk\n9pQpLl6kByWmiYj4qB07dhAYGMjyhx7i1tBQWLMGihRJVGfzZmjVCn76CRo2zKZLV4SEhBAWFubt\nZoiIeE316tX55ptvaDZ5MicaNUpy5lGdOjB4cBgtW4ak+eephyAi4uM++ugjPv/kE9aXLMkNjRrB\n229fUWfkSBg2TPshiIhka9ZannjiCc6EhzN51y5MSAj07XtZHciVS7eMRESyNWMMY8eOJfzMGca2\nbg1DhsD69Ze+HxYWxogRIWn/Ob74l7h6CCIiVzpy5AiNGjVicpcuNJ0+3a1hUbr0pe9rC00RkRxk\n06ZNtG7dmi1dulBm+3ZYtAjy5AG0QY6ISI5St25dPvroI+6aP5/zBQrAs8+m27kVEEREspgePXrw\nQM+edDt7FrtgAXz9dbqc12cDggaVRUSu7vXXX+diwYK83aQJYYMGETJgQJrPqTEEEZEsKjIykoYN\nG/LVfffR9IcfML//rjEEEZGcqHjx4sycOZPOU6ZwOB1WhFZAEBHJwurWrcv7779P02XL0nwuBQQR\nkSyuT58+tA8OTvN5fDYgaFBZRCR5wsLC8PPzS/N5NKgsIpJNKDFNRETShQKCiIgACggiIuLhswFB\ng8oiIskTFhZGSEhIms+jQWURkWxCg8oiIpIuFBBERARQQBAREQ8FBBERARQQRETEQwFBREQAHw4I\nykMQEUke5SGIiEgiykMQEZF0oYAgIiKAAoKIiHgoIIiICKCAICIiHgoIIiICKCCIiIiHzwYEJaaJ\niCSPEtNERCQRJaaJiEi6UEAQERFAAUFERDwUEEREBFBAEBERDwUEEREBFBBERMRDAUFERAAFBBER\n8VBAEBERQAFBREQ8fDYgaHE7EZHk0eJ2IiKSiBa3ExGRdKGAICIigAKCiIh4KCCIiAiggCAiIh4K\nCCIiAiggiIiIhwKCiIgACggiIuKhgCAiIoACgoiIeCggiIgIoIAgIiIeCggiIgIoIIiIiIcCgoiI\nAMkMCMaYL4wxfxpjNl92vJ0xZocxZpcx5p9JvK6TMWa8MeZbY0wbz7EAY8xMY8znSb1GRES8I7k9\nhEnAPQkPGGNyAR95jtcCehljqiesY62dY63tDwwEHvAcrgNMt9b+DaiXhraLj9LWp1mXrl3OlqyA\nYK1dDhy/7HAjYLe1dp+19gIwFeh0lVMMBcZ5vl4BDDDG/B/wc8qbLL5OHypZl65dzpaWMYQAYH+C\n5wc8xzDG9DXGjDbG+Btj3gTmWWs3euo9CrxsrW0N3JeGn5/h0vOXI7XnSsnrklP3enWu9v2UHve2\n9G5XVrh+Kf2er147yHrXL7v87mXIoLK1drK1dgjQFWgFdDPG9Pd8ezHwjDHmE2Dv1c4REhJyqXjr\nP64CQuqPe1tW+0BJbl0FhMw9n68HhLCwsESflWllrLXJq2hMBSDUWlvH87wJEGKtbed5/iJgrbVv\npblRxiSvUSIikoi11qT2tTekoK7xlDjrgCqeQHEI6An0Sm1DEkrLGxIRkdRJ7rTTKcBKoKoxJtwY\n84i1NgYYBCwAtgJTrbXbM66pIiKSkZJ9y0hERLI3ZSqLiAjgxYCQ0uxnY0x1Y8wnxphpxpjHMr/F\nEic1mevGmALGmHXGmPaZ21q5XCp+9wKNMcs8v3/NM7/FEicV184YY14zxnxojOl7vfN7s4eQouxn\na+0Oa+1A3OB120xuqySWmsz1fwLTMq2Fci0pvX4WOAXkw+Ubifek9Np1AsoB0STj2nktIKQm+9kY\nEwz86DkuXpLSa2eMaQ1sA46QeKaaeEFKr5+1dpm1tgPwIvBqpjZWEknF52Y1YIW19nngyeud39fG\nEK6a/QxgrQ211rYH+mVyu+T6rnXtgoDGQG/gb5nbLEmma/7ueUQBeTOtRZJc17p2B4gPIDHXO1FK\n8hC8yhgTCHQB8gNLvNwcSQFr7VAAY8xDwFEvN0dSyBjTGXc7ogju1oRkHTOBscaYZsDS61X2tYAQ\nAZRP8Lyc5xjW2qUk4w2J11z12sWx1n6dqS2SlLjW794sYJY3GiXJcq1r9xcp6JV7+5bRVbOfjTF5\ncQPIP3ilZXI9unZZm65f1pVh186b006V/ZxF6dplbbp+WVdGXztlKouICOD9W0YiIuIjFBBERARQ\nQBAREQ8FBBERARQQRETEQwFBREQABQQREfFQQBAREUABQUREPP4fZNwgeS5AzO4AAAAASUVORK5C\nYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdd013278d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(df.mEt2,df.muegamma,'k-')\n",
    "plt.loglog(df.mEt2,df.muegamma_analy/4.,'r-')\n",
    "plt.loglog(df.mEt2,df.muegamma_analyT/4.,'b-')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scan upont $Y_{12}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "commands.getoutput('cp LesHouches.in.FT_SD_3gen LesHouches.in.radinertRSIII')\n",
    "a=hep(MODEL='radinertRSIII',spc_input_file='LesHouches.in.radinertRSIII')\n",
    "df=pd.DataFrame()\n",
    "a.LHA.blocks['SPHENOINPUT'].entries[55]='0               # Calculate one loop masses'\n",
    "dm_masses=np.logspace(np.log10(1E-6),np.log10(1E-0),100)\n",
    "for MHX in dm_masses:\n",
    "    if np.where(dm_masses==MHX)[0][0]%10==0: #find the index of the array entry\n",
    "        print np.where(dm_masses==MHX)[0][0]\n",
    "    a.LHA.blocks['YNIN'][1,2]='-%0.8E  # Yn(1,2)' %MHX\n",
    "    a.runSPheno()\n",
    "    dp,y=lfvin(a.LHA_out_with_comments)\n",
    "    a.Series=a.Series.append(pd.Series({'muegamma_analy':LFV(SM,dp,y)[0],'muegamma_analyT':LFV(SM,dp,y.transpose())[0],\\\n",
    "                                         'muegamma':eval(a.LHA_out_with_comments.blocks['FLAVORKITLFV'].entries[701])}))\n",
    "    df=df.append(a.Series,ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fdd06413e50>]"
      ]
     },
     "execution_count": 100,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xl4VdXVx/HvYhJxhpYqaIFWFMVSqApaHNJWC2LVVgYB\nFQFtK31ffa0VUATBeWZMCAghSmRSK62gT+tEoEBtcUTDEJkpWhxQMhAy7vePe4LXEMhN7nDu8Ps8\nD8+Tu++956xtYlbW2ffsZc45REREGvkdgIiIxAclBBERAZQQRETEo4QgIiKAEoKIiHiUEEREBFBC\nEBERT0gJwcyyzGy3ma2tMd7bzDaYWb6Zja7lfR3MbLaZPRc0drGZrTCzTDO7KPwpiIhIJIRaIWQD\nvYIHzKwRkO6NdwYGmVmn4Nc457Y6526qcSwHFAJHAP9pSNAiIhJ5ISUE59xK4Ksaw92Bj51z251z\n5cBC4KoQjrXCOXc5cCdwXz3jFRGRKAlnDaEtsDPo8X+8MczsejObaGYnec9ZLe//GmgWxvlFRCSC\nmkTjoM65HCDHzFqaWSbQ1cxGO+ceNbPfELjMdByBS04iIhIHwkkIu4DvBz0+2Rs7wDm3BxhRY2wx\nsPhwBzYz7bgnItIAzrnarsiEpD6XjIxvX/pZA5xqZu3MrBkwEHipoYHU5JwL69/48ePDfl1tz9U1\nVvP56seHe008z+9wj0P5Oprzq+/c4mV+0freRWJ+ifSzmWrzC+V3S7hC/djpfGA1cJqZ7TCzYc65\nSuAW4FUgD1jonFsfdkQRkpaWFvbranuurrGaz1c/DjWeUMVqfod7fLh5hyuU49V3brWN+zG/aH3v\nahtPpvnV9+c12eYXk98t4WbLaPwLhJW8xo8f73cIUaX5Ja5knptzyT8/73dng3/36k5lH0T6r+l4\no/klrmSeGyT//MJlLgLXnSLNzFw8xiUiEs/MDBejRWUREUliSggiIgIoIYiIiEcJQUREACUEERHx\nKCGIiAighCAiIh4lBBERAfxpoWlm9oCZTTWz68OfgoiIRIIfLTSvIrBVdhlqoSkiEjdi3kITOB1Y\n5Zy7A/hDfYIVEZHo8aOF5k6+SS6VYZxfREQiKCqLys65HOfc7UBpcAtN7+nFQG8zmwIsj8b5RURS\ninNUTJ8e9mH8aKFZAtRcVzhIWloa7du3p3379qSlpWnbWhGRGnJzc8n9y1/InzePj7/8Muzj1Sch\nHLKFJvApgRaag8KOyJObmxupQ4mIJB/nuGDdOs6dNYtJZvSePZsbbrwxrEOG1A/Ba6GZBrQCdgPj\nnXPZZnYZMJnApacs59wjYUXzzfnUD0FE5FC2b6fwmmvYsXYt6eecw9gFC2jbtm3Y/RDUIEdEJFF4\nawWlo0YxyYxTpkxhyPDhmAVyQLgJIZw1BBERiZWgqiDj3HO5e/582rZtG9FTaOsKEZF45hwVGRkU\nn3kmUz76iHfS08nIzY14MgBVCCIi8SsGVUEwVQgiIvEmhlVBMFUIIiLxJMZVQTBVCCIi8cCnqiCY\nKgQREb95VcH2tWuZHuOqIJgqBBERvwRVBZM/+oh3fagKgqlCEBHxw/btFA4YwPYPP/S1KgimCkFE\nJJaCq4K8PN+rgmAhVQhmlgX8CtjtnOsSNN6bb+9l9GiN93UA7gaOdc4N8MYuAK71zn2Gc+6CSExE\nRCTuHWIPongR6uZ2FwBFwNzqhOC10MwHfgF8QmD304HOuQ21vP+56oQQNHYV0No5N6uW12svIxFJ\nHs5RkZlJ6ciRTDLj+1Oncv2wYQf2IIqUmOxl5Jxb6W1zHexAC00vkOoWmgclhEMYDAwPNVARkYQU\n51VBMD9aaGJmpwBfO+eKwzi/iEj88nYmDb6vYPry5XGbDCBKnzJyzuUAOWbWMriFZtAaw41AdjTO\nLSLiuwSqCoLFvIWmNz6hroOrhaaIJJygfgVTzPh+ejrTo7BWUC03N5fc3Fy2bdvGtm3bwj5eyA1y\nzKw9sMQ59yPvcWNgI4FF5U+BfwODnHPrww5Ki8oikmh83IOoWriLyiGtIXgtNFcDp5nZDjMb5pyr\nBG4BXgXygIWRSAYiIgmlxt3GfuxBFClqoSki0lBxsgdRtZhUCCIiEiTO9iCKFO1lJCJSH3FWFUSS\nKgQRkVAkaVUQTBWCiEhdkrgqCKYKQUTkUFKgKgimCkFEpDZB/QoyEuhu43CoQhARCVZLv4J434Mo\nUlQhiIhUS5G1gkNRhSAikmJrBYcS6tYVWWa228zW1hjvbWYbzCzfzEbX8r4OZjbbzJ4LGmtrZi96\n4we9R0QkprZvp/D889k4ciQjzzmH4Rs3MmT48KhtSBfPQq0QsoFewQNex7R0b7wzMMjMOgW/xjm3\n1Tl3U41jdQFe8Ma7NihqEZFw1dKvIBWrgmB+dExbBSwxs+FATj3jFREJX4L2K4g2PzqmDQfuds5d\nAvwqjPOLiNRP0FpBonQxiyU/Oqa9CdxjZtcCW6NxfhGRg6T4J4hCEfOOac65tUC/ug6ujmkiEhFB\nXcwmm9EuPZ2MKHYxiyV1TBMRCVXQ3capUBWoY5qISE213G2c6p8gCoU6polIcqlRFYxdsIA2bdr4\nHVVMqGOaiAgcVBW8l5FBRm5uyiSDSNBeRiKS+PQJoohQhSAiiavGHkTvTJumtYIwqEIQkcQUVBWk\nSr+CaFOFICKJpZadSXW3cWSoQhCRxBHcxezccxmrtYKIUoUgIvGvti5mWiuIOFUIIhLfVBXEjCoE\nEYlPqgpiThWCiMQf3VfgCz9aaJ5hZovMLMPM+oY/BRFJGupt7KuQ9jIyswuAImCuc66LN9YIyCew\n2+knwBpgoHPuoI5pZvacc26A9/XtwL+cc6vM7K/Ouatqeb32MhJJNaoKwhaTvYyccyuBr2oMH2ih\n6ZwrB6pbaNYlBxhoZo8BLesTrIgkIVUFcSOcNYTaWmh2h0ALTaAb8Lhz7lOCWmg65z4HbvEqjD+H\ncX4RSXQp1q8g3kXlU0bOuRzn3O1AaXALTQAza2dmM4FngMejcX4RiXOH2JlUycBffrTQ3A78vq6D\nq4WmSJJK4X4FkaYWmiKSmIJ6G080o/20aVw3dGhS9DaOF+EuKodUIXgtNNOAVma2AxjvnMs2s+oW\nmo2ALLXQFJFa1fgEkaqC+KQWmiISPTWqgnZTp3L9sGGqCqIkJhWCiEi9aa0g4WgvIxGJLPU2Tliq\nEEQkcrZtC6wVeDuTjlNVkFBUIYhI+JyjIj2d4s6dD1QF01UVJBxVCCISnm3bAmsFH33E9O7dGTt/\nvhJBglKFICINE1wVrFvH+9Onk7FsmZJBAlOFICL1p6ogKalCEJHQqSpIaqoQRCQ0qgqSXrQ7pl1l\nZk+Z2QIzu9Qba2FmT5vZTDMbHJlpiEjUqCpIGbHqmHY8gd4IvzWz64CvnHMvm9lC59zAWl6vrStE\n4sG2bRR4VUGmqoK4lygd08YC6d7XJ/NNY53K+oUrIjFRVXWgKpiybh0fqCpICdHumPYEcCvwinPu\nA+91OwkkhbUEdVITkTihqiBlRbtjWl8Cl5T6mdnvvKcXe48zgCXROL+INEBVFRXTplF85pmqClJU\ntDumTQOm1RjbBwyv6+DqmCYSQ0FVwYwePRg7fz4nnXSS31FJHdQxTUQip6qKiowMSkePZmKjRnRI\nT+faG25Qv4IEpY5pItIw3n0F21QViEcd00RSTXVVcOedPGnGD1QVJA11TBOR0G3dSsE11xyoCsap\nKpAg2stIJBVUVVExdeqBu43XZmaS8eabSgbyLaoQRJLd1q0UDBjAtrw8VQVyWKoQRJJVVRXl1VXB\n+vV8OGOGqgI5LFUIIsnIqwq25uUx87zzGDdvnhKB1EkVgkgyqaqifMoUijt3ZtL69eTNnEnGG28o\nGUhIVCGIJIutWyns358t69apKpAGUYUgkuiCqoInN25UVSANpgpBJJFt2RJYK1BVIBGgCkEkEVVX\nBWedxURVBRIhfrTQ7GBms83suchMQSTFbNlC4bnnknfXXYw6/3xu/vhjBl9/vbaekLDFvIVm0Nhz\nzrkBhzif9jISqcm7r6BszBgea9yYTpmZDLz2WiUCOSAmexk551aaWbsawwdaaHqBVLfQPCghEGih\nmdHQIEVS3ubNFA4YwKb165l1/vncM28eJ554ot9RSZIJZw2hthaabSHQQtPMJppZGzN7hEALzfdr\nvF9/1ojUpaqK8kmTKP7Rj3g8P5+Ns2aR8frrSgYSFVH5lJFzLgfI8fol/AI41sxOdc49ZWYtgQeB\nrmY22jn3aDRiEEl4mzdT0L8/mzdsYPZPf8q4Z59VIpCo8qOF5h5gRF0HVwtNSVlVVZRPnkzp2LE8\n0bgxZ86aRfrgwVorkIOohaZIMtu0iYIBA9i8YQNZF1zAuJwcvve97/kdlSSIcBeVQ/3Y6XxgNXCa\nme0ws2HOuUqguoVmHrBQLTRFGqiqivInn6S4Sxee+Phj8mfPZtrf/65kIDGlFpoiftu0iYL+/fl4\n40ayL7yQe3JyaN26td9RSQKKSYUgIlFQVUX5449T3KULj2/axOasLKb97W9KBuIb7WUk4oePP6ag\nf3/y8/N55qKLGDd3rhKB+E4VgkgsVVZS9thjFP/4xzy2ZQtbs7NVFUjcUIUgEiv5+RT068fGTZuY\nm5bGPc88w3e/+12/oxI5QBWCSLRVVlL+yCMUd+3KY9u2sS07m2mvvKJkIHFHFYJING3cSEG/fqzf\nvJn5P/85Y7OzlQgkbqlCEIkGryoo6taNh7dv5z/PPMOUpUuVDCSuqUIQibQNGyjo1491W7aw4Be/\nYFx2Nt/5znf8jkqkTqoQRCKlooKyBx6g6Cc/4aGdO/lk7lymLFmiZCAJQxWCSCSsW0dBv37kbdvG\nwksu4Z7sbFq1auV3VCL1EnKFEKE2mpfUMnZp+NMQ8UlFBWX33UfR2Wfz4K5d7H72Waa89JKSgSSk\n+ux2Go02mgeNeePay0ji30cfUdivH2t37OCFXr0Yl5VFy5Yt/Y5KUljM9jJyzq0EvqoxfKCNpnOu\nHKhuo1mb2tpoqrWmJJ6KCkrHj6fw3HN5YPdu9ixcyKTFi5UMJOGFu4ZQWxvN7hBoowl0A54AbqVG\nG83DtNYUiV8ffkhhv368v3MnL/bpwz2zZ3PCCSf4HZVIRETtU0bOuRzn3O1AXwKXlPqZ2e8Aglpr\nHhgTiWvl5ZTecw+F3btz7+efs3fRIib9+c9KBpJUwq0QGtpG86CxmtRCU+LGBx9Q2K8f7+zaxZIr\nrmDcU09x/PHH+x2ViH8tNCF2bTS1qCxxoayM0vHjKZ08mQlHHsmlzz7LZX36+B2VyCGFu6gccoXg\ntdFMA1qZ2Q5gvHMu27v88yqBy09ZaqMpSeHddyns359/f/opf/vNbxifmclxxx3nd1QiUaUWmiLB\nSkvZP3YsZenpTGjRgl7z5tGrd2+/oxIJScwqBJGk9/bbFPbrxz937+b1vn2ZMH06xx57rN9RicSM\n9jIS2b+fkj/+kb0XXsi44mIaL1nCY88+q2QgKUcVgqS2f/2Lwv79+cfnn/OPgQO5f+pUjjnmGL+j\nEvGFKgRJTSUllNxyC1+npXF3aSktXnmFh7OzlQwkpalCkNSzahWFAwawbM8e3rruOh6ePJmjjjrK\n76hEfKcKQVJHcTH7fv979lxyCXc5R8vXXuOhWbOUDEQ8qhAkJbhlyygaNIhX9+7lg+HDeezxx2nR\nooXfYYnEFSUESW6Fhey79VZKFi5kbKtW3LBsGX3PO8/vqETikhKCJC332msUDx7MkqIi8m++mUkP\nP0zz5s39DkskboW0hhChbmmXemOdzCzTzBaZ2Y2RmYZIkL17Kb7uOr648kruOOooOv7jH4yfNEnJ\nQKQOIW1dEaVuaQYsdM5dU8vrtXWFNIhbupR9Q4bw55IS/nPrrdxx//00a9bM77BEYiImHdMi3S3N\nzK4AXvbeIxK+PXso7tuX/w4YwJ9ataLLP//JmEcfVTIQqYdwPnZaW7e0thDolmZmE82sTW2d0Zxz\nS5xzfYChYZxfBICqF16g+Ac/YP4rr/DsyJFMW7eOrl27+h2WSMKJyqKycy4HyAnqjHasmZ3qnHvK\nzC4GrgaaA8uicX5JEZ99RuHQoXyVm8uDHTrwf88/z5lnnul3VCIJK5yE0NBuacuB5XUdXB3T5JCc\no2rePEpGjODpigrchAlMv+MOGjdu7HdkIjHlW8e0WHVL846tRWWp3a5dFF5/PZ+99RaPderEHYsW\n0bFjR7+jEokLMVlU9rqlrQZOM7MdZjbMOVcJVHdLyyPwiSF1S5PocI6Kp56i+PTTyXzrLV5/9FEy\n335byUAkgtQxTeLf1q0UDhzIzg8/ZGrXrty1YAHt2rXzOyqRuBOTCkHEF5WVVEycSHHnzkzKy2PN\ntGlkrlqlZCASJdq6QuLT+vUUXHMN+Zs2kf3TnzI2J4eTTjrJ76hEkpoqBIkv5eWUTZhA0dln89D2\n7WzNzib9tdeUDERiQBWCxI9336VwwADe+eQT/tKrF+Nmz6ZVq1Z+RyWSMlQhiP9KSth/220U9OzJ\nuK++Yt/zzzN58WIlA5EYU4Ug/lqxgqKBA1n29desGjiQ+6ZM4dhjj/U7KpGUpApB/FFQwL6hQ/my\nd29GmXH83//OI9nZSgYiPlKFIDHnXnqJfUOH8mJJCVtHjGDiQw+pV4FIHNCNaRI7u3dTNHw4e5ct\n476TT+bmRYvo1q2b31GJJI1wb0xThSDR5xyVc+ZQetttzKmooHLcODJGjaJJE/34icSTmLfQ9MZb\nmNkaM+sT/hQkrm3eTOH555N/663c1qkTfT78kD+OGaNkIBKHfGmhaWb3AoXAOufcK7W8XpeMEl1F\nBeWPPUbZAw/weKNGtJ88mRtuvJFA51QRiYaYXDJyzq00s5obyBxooekFUt1C86CEQFALTTO7BFhH\noEGOfjsko3feoXDgQNbu2sWin/2MMVlZnHjiiX5HJSJ1CKdur62FZncItNAEugFPALfy7RaaaUAL\noDOwj0BvZUkGxcWUjB5NWVYWE1q04BeLFjH1iiv8jkpEQhTzFprOubEAZjYE+CIa55fYcy+/TPHQ\nofytqIh3r7uO+yZO5JhjjvE7LBGph5i30Ax6bu7hDq4Wmgniv/+l6Le/pfCNN5hw4okMX7qUfj16\n+B2VSEpQC02JD1VVVMycSenIkcyurKRizBhuvfNOmjZt6ndkIikrJovKXgvNNKCVme0Axjvnsr1L\nQq8S+Phqllpopoi8PAoHD2bLpk1M//GPGZWTww9/+EO/oxKRMOlOZQndvn3sHzuW8sxMHmzWjC4Z\nGQy69lp9lFQkTuhOZYkJ98orFA8bxhsFBazs25e7pk6lZcuWfoclIhGkCkEOb9cuim66icLly5nQ\nujVD5s2jZ8+efkclIrVQhSDRUVFBxeTJlI0fzwznaHT33aSPGqVFY5EkpoQgB1u9msIhQ/jo00/J\nOe88Rs+ZQ7t2NW9UF5FkowY58o0vvqB48GD2XHIJY/bu5bN585j+xhtKBiIpQhWCQGUllTNnUjp6\nNDkVFXz2P//Do/fdR4sWLfyOTERiSAkh1f373xQOGUL+zp3MOOss/vTMM3Tq1MnvqETEB7pklKo+\n+4x9gwfzVVoa4z7/nC3Z2Tz11ltKBiIpTBVCqqmooGLaNMrGjePpykp233wzD9x/P0cffbTfkYmI\nz2LeMc3MLjazFWaWaWYXRWYaEpI336SwY0fWjB3LH7t145IPPuDeSZOUDEQECP2SUTbQK3jA65iW\n7o13BgaZ2beuNzjn/uqc+x0wAhhQPUygW9oRBHooSLRt3UrRL3/Jf6+4glElJXw+fz4zVqzgtNNO\n8zsyEYkjMe+Y5pxbAawws9bAROC6BsYudSkooHTCBCoyM5lixhFjxjB55EiOOOIIvyMTkTjkR8e0\nal8DzcI4vxxKZSVVc+ZQMnIkS0tLeevKKxk5aRJt2rTxOzIRiWMx75hmZr8hcJnpOAKXnCSSXn+d\nwt//no9372b6qady8+zZXHPOOX5HJSIJIOYd05xzi4HFdR1cHdPqKS+PohEjKHjnHe4/+mh+PmcO\ns/r319bUIklMHdPk23btYt+oUVS++CKPNm7Md+65hxH/939aJxBJQeqYlqq+/pqyBx+kIiODWc7x\n2fDh3H7ffbRq1crvyEQkQakfQqIpKQlsS/3QQ7xUWcmqX/6S2ydNokOHDn5HJiI+C7dCUEJIFGVl\nVGVlUTJ2LCv37+fFn/yEm6dOpVu3bn5HJiJxQgkh2VVU4HJy2HfnnbxfXMycDh0Ykp7OxRdf7Hdk\nIhJn1DEtWVVU4ObPp/iuu9hYUEDmSSdxdXY2sy+7TJ8cEpGoUEKIN+XluHnzKB4zhvzCQjK++10u\nf/ppZl19tRKBiESVEkK82L8/cHfxvfeyft8+Mlu35vL0dGb9+tc0aqRdykUk+pQQ/FZQQOWMGZQ+\n8ghrysuZe/LJXDFzJrOuvFKJQERiSovKfvnkE8qefJLKGTN41Ywlp5/OgIcf5tJLL9WlIRFpEH3K\nKNG89x77HnoIli7lWeC9tDSGTphAjx49/I5MRBKcEkIiKC+Hv/6VwocfpmzDBqZWVVE0aBA333UX\nHTt29Ds6EUkSSgjx7NNPqZgxg7L0dPIrKpjZuDEd/vQnbhoxgpYtW/odnYgkmZgkBDPLAn4F7HbO\ndQka7w1M5pu9jB6t8b6rgMuBY4A5zrnXLHCB/H7gWGCNt1V2zfMlbkKorIRXX6V48mQarVjBC2Ys\nP+ssLhs9miuvvJKmTZv6HaGIJKlYJYQLgCJgbnVC8Fpo5hPY7fQTYA0w0Dl3UMc0MzseeNw591sz\n+zXwa+AL4GXn3LJaXp94CWHTJiqysiibPZv/lJWRWVlJ0+uv54b//V86d+7sd3QikgJicqdyJFto\nAqcDq5xzs8zseeCghJAwvvwS9/zzFM2YgcvPZ55zvN25Mz+77TYe6tuXI4880u8IRURC5kcLzZ1A\nmfd1ZRjn98fevbB0KUVZWTRZvZo3mzXjLy1acModd3DdsGGM0K6jIpKgYt5Ck0C3tGlmdiGwPBrn\nj7gvvoClSymeO5cmq1fz7+bNWeAcRw0dSt+hQ5nZo4fuHRCRhOdHC80S4Ka6Du5rC03nYMMG3NKl\nFC1YQNN161jVvDkvOEfza6/l8kGDmJqWRpMmutFbRPyjFprRsmcP5OZStmQJFUuXsr+khFecI/e4\n4/jOgAH0ufpqevbsSePGjWMbl4hIiNRCs6H27IHVq6l8801KXn6Zptu28f5RR7G4uJgvzjmHM6++\nmsv69OHaTp10OUhEUkJq3JhWVQX5+fDWW5SvXEnpm2/SdNcu8o4+mr8VFbG5QwdaX345F196KRdd\ndBEtWrSI3LlFRGJEdyrXVFER+OX/wQdUrllDycqVNP3oIwqaNOHtxo1ZVlzMl6efTuteveiZlkbP\nnj054YQTIjsBEREfpG5CKCuDLVtg/XpYv57S99+n7L33aL59O3uOPJK8xo1ZXlTEru99j6bdu3PG\nxRfTo0cPunbtyhFHHBGbiYiIxFDyJoSqqsDHPXfuhK1bYetWKjdtoiQvD9u0iSO++II9LVqwqUkT\n3ispIc85Sk89laO6d6fT2WfTpUsXunbtytFHH+33dEREYiJpE0J506aUNWnC582bs8OMDaWlrN+/\nn8LWraFjR44+6yx+0KkTZ5xxBmeccQZt2rTR4q+IpLSkTQh/uOEGTvzhDznllFNo3749HTp0oG3b\ntvrsv4jIISRtQojHuERE4lm4CUFNe0VEBFBCEBERjxKCiIgASggiIuIJKSGYWZaZ7TaztTXGe5vZ\nBjPLN7PRtbzvKjN7yswWmNml3tgFZpZpZrPMbGVkpiEiIuEKtULIBnoFD3gtNNO98c7AIDPrFPwa\n59xfnXO/A0YAA7yxlc65EcBS4Jnwwk9Mubm5focQVZpf4krmuUHyzy9cISUE59xK4KsawwdaaDrn\nyoHqFpq1CW6hWW0wML8esSaNZP+h1PwSVzLPDZJ/fuEKZw2hthaabSHQQtPMJppZGzN7hG+30MTM\nTgG+ds4Vh3H+wwr1G3+419X2XF1jNZ+vfhzpH8RYze9wjw8373CFcrz6zq22cT/mF63vXW3jyTS/\n+v68Jtv8YvG7JSqLys65HOfc7UBfAg10+pnZ74JeciOBy1BRk8zftPocTwnh8OPJ9AultvFkmp8S\nQm6tz0VybvXpmNaOQMe0Lt7j84AJzrne3uM7AeecezTsoMx0m7KISANEvWOax7x/1dYAp3qJ4lNg\nIDCooYEEC2dCIiLSMKF+7HQ+sBo4zcx2mNkw51wlUN1CMw9YmFAtNEVE5FvicnM7ERGJPd2pLCIi\ngA8JIYy7ns3MHjCzqWZ2fewirp8w5nexma3w7uK+KHYR109D5+e9poWZrTGzPrGJtv7C+P518r53\ni8zsxthFXD+R3HUgHoUxvw5mNtvMnotdtA1T3zl6/989bWYzzWzwYQ/unIvpP+ACoCuwNmisEbAJ\naAc0Bd4HOtV436+Bp4EngJ/FOu4YzO8i4GVgDvADv+cR6fl5r7sXuAPo4/c8ojE/77UGLPJ7HlGc\n3/HALL/nEcX5Pef3HCI9R+A64HLv64WHO3bMKwTX8LueTwdWOefuAP4Q/UgbpqHzc86tcM5dDtwJ\n3BeTYBugofMzs0uAdcDnfPvTanEljJ9PzOwKAkl9YdQDbaBw5uepbdeBuBGB+cW9BszxZL65ibjy\ncMeOlzWEOu96Bj7hm/8Ih51UHArlru6TvOe+BprFOL5w1TW/SQQ+ktyDwJYlN8U8wvCE9P1zzi1x\nzvUBhvoQYzgavOtAgqjP/39x+8dKHQ45R+/rk72vDzu/uG9Q7JzLAXLM7EhgmpldCCz3OayICZrf\nb8ysF3AcgU0Dk0L1/Kofm9kQ4Av/IoqsoO/fxd7Nmc2BZT6HFTFB87uFwK4Dx5rZqc65p3wOLSKC\n5tfSzDKBrmY22kXgBts48iKQbmaXA0sO98J4SQi7gO8HPT7ZGzvAOVdC4v1lWS2U+S0GFscyqAiq\nc37VnHP0CQHKAAAArUlEQVRzYxJRZIXy/VtO4v6hEsr8pgHTYhlUBIUyvz0EdmVOVIeco3NuHzA8\nlIP4dcnokHc9m1kzAnc9v+RLZJGh+Wl+8UzzS+z5QbTm6MMK+XwC6wGlwA5gmDd+GbAR+Bi40++V\nfM1P89P8Eu9fss8v2nPUncoiIgLEz6eMRETEZ0oIIiICKCGIiIhHCUFERAAlBBER8SghiIgIoIQg\nIiIeJQQREQGUEERExPP/YbLMujO0wMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fdd067a46d0>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(np.abs(df.Yn12),df.muegamma,'k-')\n",
    "plt.loglog(np.abs(df.Yn12),df.muegamma_analy/4.,'r--')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Scan upont the right handed neutrino mass "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "10\n",
      "20\n",
      "30\n",
      "40\n",
      "50\n",
      "60\n",
      "70\n",
      "80\n",
      "90\n"
     ]
    }
   ],
   "source": [
    "commands.getoutput('cp LesHouches.in.FT_SD_3gen LesHouches.in.radinertRSIII')\n",
    "a=hep(MODEL='radinertRSIII',spc_input_file='LesHouches.in.radinertRSIII')\n",
    "df=pd.DataFrame()\n",
    "a.LHA.blocks['SPHENOINPUT'].entries[55]='0               # Calculate one loop masses'\n",
    "dm_masses=np.logspace(np.log10(100),np.log10(1E5),100)\n",
    "for MHX in dm_masses:\n",
    "    if np.where(dm_masses==MHX)[0][0]%10==0: #find the index of the array entry\n",
    "        print np.where(dm_masses==MHX)[0][0]\n",
    "    a.LHA.blocks['MTFIN'][1,1]='%0.8E  # MTF(1,1)' %MHX\n",
    "    a.LHA.blocks['MTFIN'][2,2]='%0.8E  # MTF(2,2)' %(MHX+1000)\n",
    "    a.LHA.blocks['MTFIN'][3,3]='%0.8E  # MTF(3,3)' %(MHX+2000)\n",
    "    a.runSPheno()\n",
    "    dp,y=lfvin(a.LHA_out_with_comments)\n",
    "    a.Series=a.Series.append(pd.Series({'muegamma_analy':LFV(SM,dp,y)[0],'muegamma_analyT':LFV(SM,dp,y.transpose())[0],\\\n",
    "                                         'muegamma':eval(a.LHA_out_with_comments.blocks['FLAVORKITLFV'].entries[701])}))\n",
    "    df=df.append(a.Series,ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f642c476c10>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEECAYAAAAoDUMLAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3XmcjvX+x/HXZ2xjaWzZwzgkB0mmiKMaUSFFkq2UrFGU\nnKISUqdoE6Lse0P8pGxHccw4KpkUFdmyq0gRsg7f3x9z60yyzMw9933d98z7+Xh4PFzf+7qu+63r\nMX3m2r4fc84hIiIS4XUAEREJDSoIIiICqCCIiIiPCoKIiAAqCCIi4qOCICIigAqCiIj4pKogmNl4\nM9trZl+fM97QzDaY2SYz63Oe7Zqa2RgzizOzW1OM5zGzRDNr7P8/QUREMoKl5sU0M6sLHAGmOOeq\n+cYigE1AfeAHIBFo7ZzbcJ7tCwCvOuc6+5afBw4D651zCzPo3yIiIn5I1RmCc24FcOCc4ZrAZufc\nDufcKWAG0PQCu+gHjAQwswbAeuBnwNITWkREMl52P7YtBexKsbyb5CKBmbUDrgVeA3oCC51za3zr\nxQJ5gCrAUWCBHxlERCSD+FMQLsg5NxWYamY9SL6kFGVmFZxzY5xz/QDM7AFgfyC+X0RE0s6fgrAH\nKJNi+Qrf2B+ccyOAEefb2Dk35UI7NjPNuCcikg7OuXRfik/LY6fGn6/5JwIVzKysmeUEWgMfpjfI\nucY89BBjr7iCzyIi+D1bNrZVqsSGxETOnDmDcy4ofwYMGOD5vtKyXWrWvdQ6F/o8LeMZ+d8tFI5d\nuBy/tH4WqscuHI9fqPzs+StVZwhm9i7J1/4Lm9lOYIBzbqLvktBHJBeW8c657/xO5NN5wgQADh48\nyJIPP2Tn9Om82rw5mNGoUSMaNWrELbfcwmWXXZZRX/kXsbGxnu8rLdulZt1LrXOhz9M67rWMzhUO\nxy+tn4XqsYPwO36Z5WcvVY+dBpuZufPlcs7x3XffsWjRIhYtWsTnn39O66pVefbgQfK0aEHRDh2g\nXDkPEktKAwcOZODAgV7HkHTQsQtvZobz45JRWBWEcx0+fJj4+fPZPX48BT/9lAYnT5JUoAAnb7uN\nEr16keP664OQVs4VHx8f0r99yoXp2IW3LF0QUnLO8dXq1Xw5ejQ2fz5rDh3i12bNaNq0KQ0bNiQq\nKipAaUVEQoMKwgXs2bOH+fPn88EHH7BixQrq1KlDs2bNuKdMGYrUrAmXX55BaUVEQoMKQiocPnyY\nxYsXM3fuXP4xezbtkpL4tWxZcrdtS5EuXaB06Qz7LhERr6ggpNGpU6dYvngxm0aOpFB8PLefOsWR\nEiU4PHcuf4+JCch3iogEgwqCH06fPs1ny5fzxdtv8/pnn5EvXz7uvfde7r33XqpWrYqZploSkfCh\ngpBBzpw5Q2JiIrNmzWLWrFnkyZOHR2+8kRb58lGsRw89zioiIU8FIQCcc6xatYr4d96h5KxZND5+\nnOPFi5PzgQco8sgjUKqUZ9lERC5EBSHAzpw5w2f//S9rhw6l0OLFNDp1iv+2bs21Q4ZQSoVBREKI\nCkIQJSUlsfyjj3gvLo73FiygevXq3Hfffdxzzz0UKFDA63giksWpIHjk+PHjLFy4kOnTp7NkyRJu\nb9CAl0+coPQjj5DzttsgWzavI4pIFqOCEAIOHjzI+3FxHH3tNW7csYPSkZEcbd6ckn36YFWqeB1P\nRLKIoBQEMxsPNAH2Ol9PZd94Q+BN/jfb6ZBztmsK3AFcBkxwzn1syc9yvgBEAYkuuZnOud8XVgUh\npe3bt/PR669j06fT9NAhdl97LUXmzKG0Xn4TkQALVkGoCxwBppwtCGYWAWwiuSPaDyT3R2jtnNtw\nnu0LAK865zqbWTOgGcnd0hY455adZ/2wLQhnOef4/NNPmTt6NGMXLCAmJob27dtz9913kzt3bq/j\niUgmFLRLRmZWFpiXoiDcQHJfhEa+5b6AO/cswffZa8A059waM+sD/OqcG2tms5xz955n/bAvCCkd\nO3aMDz/8kAkTJvDFF1/QqlUrepUvT4U778QqVvQ6nohkEv4WhLR0TDtXKWBXiuXdvjHMrJ2ZvWFm\nJc1sMLDQObfGt94u4IDv76f9+P6wkTt3blq1asXixYv56quvKFmyJPNefJFfq1RhT4UKHBk1Co4e\n9TqmiGRx/hSEC3LOTXXOPQHcQ/IlpRZm1sX38ftAQzMbBiQE4vtDWZkyZejXrx+P//ILX8+fz3sl\nSvDpY49xuGBBdt1zD+7MGa8jikgW5e8lo4HOuYa+5QteMkpzKDN38803Ex0dTXR0NLGxsZm6acev\nv/7K3OHD2ThuHO/nyUPnzp1p3749RYoU8TqaiISw+Ph44uPj2b59O9u3bychISFo9xCiSS4IV/uW\nswEbST4D+BFYBbTJiL7Kme0eQmo551i5ciWjR49m7ty5NG7cmK5du3JTrVpYZKTX8UQkxAXlHoKZ\nvQt8ClQ0s51m9pBz7jTQA/gIWAfMyIhikJWZGbVr12bSpEls27aNG264gW7dujGzZEl+KF+eo2PH\nwokTXscUkUxKL6aFOOccy//zHxL796fGqlVclyMHx9q0oVj//lC2rNfxRCSE6E3lLOSnn35izuDB\n5Bg3jruSklg+ejRN27QhZ86cXkcTkRCggpAFJSUl8cGcObz19tts2rSJrl270rVrV4oVK+Z1NBHx\nkJfvIYhHsmfPzj0tW7Js2TL+/e9/s2fPHipVqkS7du1YFxcHW7d6HVFEwpDOEDKJAwcOMG7cOPYN\nHswzhw5xtHp1ir/8Mtnq1we1AhXJEnTJSP4kKSmJ+TNnsnnAAO7esYOoIkXI168feTp1At1rEMnU\nVBDkgr5YtYqlffpw5Sef8EmXLjzauzfl1BtaJNPSPQS5oOtq1qTPsmXU3LqV7Hnzcv3119OyZUs+\n//xzr6OJSAjSGUIWcvjwYSZMmMDQoUMpU6YMQ269lVq33ELEP/7hdTQRyQC6ZCRplpSUxOzZs1nz\n9NM8smcPOaKjKTx4MDmaNYMInTSKhCsVBEk35xzLPv6Ylf/8J3du2ECJQoXIO3AguTt3Vk9okTCk\ngiAZ4qsvv2RR796U+ewzNvbuTc/HH9dsqyJhJliT2403s71m9vU54w3NbIOZbfJ1Qjt3u6ZmNsbM\n4szsVt9YKTObY2bjzreNeOPaGjV4Ztkyaq9bx8+//MJVV11Fz5492blzp9fRRCRIUnvBeCJwe8oB\nX0/lt3zjVYA2ZlYp5TrOuQ+cc12AbkBL33A1YLZzrhNQ3Y/sEgDly5fnnXfeYd26dURGRlK9enU6\ndOjArqlTYd8+r+OJSAClqiA451bwv7aXZ9UENjvndjjnTgEzgKYX2EU/YKTv758AXc1sCfDvtEeW\nYChRogSvvPIKW7ZsoWzZsszs2pUjpUuz//77Yc8er+OJSAB40VO5A/Csc64B0MSP75cgKFSoEAMG\nDKDr3r1MfvJJZr3/PofLlePnFi1g+3av44lIBvKip/J/gMfN7G1gWyC+XzLeZZddxiMvvsiDP//M\nu/37E7d4MZ/Wrs2qVau8jiYiGUQ9lSVdjh8/zvhx4xg8ZAhVq1Zl4MCB1KpVy+tYIlmKeipLSDlx\n4gQTJ07kpZdeomrVqgwYMIBalSpB/vxeRxPJctRTWTyVK1cuHn74YTZv3sxdd91F2xYt2Fe8OD83\nbQrbdEVQJJzoxTTJUCdOnGDq8OEceeEFOhw/zsm77uLyoUOhdGmvo4lkenpTWULS8ePHmfzGG5x8\n6SUeOnmS3wYMoNSzz3odSyRTU0GQkHb06FEmDh7MuJEjqdq4MQMHDqR8+fJexxLJlNQPQUJanjx5\neGTQIBK2bePKK6+kVq1aPPzww+zRy20iIUcFQYIiKiqK/v37s3HjRqKiorj66qsZ3LEjv//rX3D8\nuNfxRAQVBAmywoUL88orr/Dtt9/y+7FjLH/+eX4rXpzjo0ZBUpLX8USyNN1DEE9t3ryZqd270zgh\ngYqFCxP11ltkb94cLN2XQUWyLN1Ulkzhy9Wreb9zZ1qsX8/2N9/kzi5diFD3NpE0UUGQTGXp0qX0\n6dMHM+PVV1/VlCUiaaCCIJnOmTNneO+993jmmWf4+9//zpAhQ6hapYouI4lcgh47lUwnIiKC1q1b\n891333HrrbdSv359VlWowOGnnoLff/c6nkimpTMECXkHDx7k7aeeouLkydyaKxe5XnmFXJ07Q7Zs\nXkcTCSmh3lO5kpm9bWYzzayjbyyPmU0ys9Fm1ja9wSXrKFCgAE+PGUPMxo28Wrs2ax97jF/KleP0\n0qVeRxPJVFJ1hmBmdYEjwJQU/RAigE0kT3/9A5AItHbObTjP9kbybKitzOx+4IBzboGZzXDOtT7P\n+jpDkAv6fOVKPnzwQXIdO8aNkydTr149ryOJhISgnCH401PZzO4EFvg+B7iC/7XePJ2e0JK11brh\nBl7csIFKr71Ghw4daN68Od9//73XsUTCXqB7Kpdwzs1zzjUG2qdY7wrf3/XYiKSLmdGyZUu+++47\natasSa1atXjqqac49NtveuNZJJ0C3VO5opkNM7PRwDLfx3NI7rE8EpgXiO+XrCMyMpK+ffvy7bff\nsn//fjr/7W/8Gh3NmSVLvI4mEnay+7HtHqBMiuUrfGN/cM4lAAnnjB0FOlxq57GxseqpLKlWvHhx\nJkyYwBeJibzWrh3d77iDyNq1uXzSJIiO9jqeSECc21PZX+qpLJmOc473Jk9mR8+edDtxArp357KX\nX4bISK+jiQRUUN5U9vVUjgUKA3uBAc65iWbWCHiT5EtP451zg9Mb5JzvU0EQvx05coSRffty2bhx\n/P7cczz25JPkzJnT61giAaOpK0QuYcuWLfTq1YtNmzYxfPhwbr/9dq8jiQSECoJIKs2fP5/HHnuM\natWq8eabb1K2TBnNjySZiuYyEkmlJk2asG7dOmrUqMEN117LnrJlOTV9OuiXDxFAZwiSRW3bto0x\n7drxYGIiRapWpXBcHFSs6HUsEb/oDEEkHcqVK8fLK1awZcYM3t62jUNXX83hXr3g2DGvo4l4RmcI\nkuUdPXqUEX36UGnMGJK6dqXZ0KFk00yqEoZ0U1kkg2zYsIGHH36YI0eO8M4773Ddddd5HUkkTXTJ\nSCSDVKpUiWXLltGjRw+aNGlCz549OXTokNexRIJGBUEkBTPjwQcfZN26dfz+++9UqVKFNQ8/DN98\n43U0kYDTJSORi0hISODjVq148sABrHNnol59FXLn9jqWyHnpHoJIgJ04cYIRzzzDlW+9xc358xMV\nF0dE/fpexxL5CxUEkSBZv349k5o3p/e2bZx++mlKDhzodSSRP9FNZZEgqVy5MoPXr+fDl14idvhw\nBg0axMmTJ72OJZJhUl0QzGy8me01s6/PGW9oZhvMbJOZ9TnPdpXM7G0zm2lmHX1jTc1sjJnFmdmt\n/v8zRIIjIiKCzr17s3TtWj7//HNiYmJITEz0OpZIhkhLP4S6wBFginOumm8sAthEck+EH4BEoLVz\nbsN5tjdghnOuVYqxAsCrzrnO56yrS0YS8pxzxMXF8cQTT9CuXTsG9e5N7mLFNGGeeCZol4yccyuA\nA+cM1wQ2O+d2OOdOATOApucJeSewwPd5Sv2AkWlKLBIizIy2bdvy9ddfs2vXLuZfdRW/1qkDu3d7\nHU0kXfy9h1AK2JViebdvDDNrZ2ZvmFkJ59w851xjoP3ZFc1sMLDQObfGzwwinipatCgzZswgx9ix\njF+3jiNXXsmJUaM0i6qEHX96Kl+Uc24qMNXMbjazvkAksAzAzHqQfJkpyswqOOfGnLu9eipLuGnW\nsiW/NmjAyw8+SJtevSg5aRKFZs+GMmUuvbFIOnjWUxnAzMqS3Ff57D2EG4CBzrmGvuW+gHPODfEr\nlO4hSJib//77rGvfnrI1anDn/PnkzZvX60iSBQT7sVPz/TkrEahgZmXNLCfQGvgwvWFEMosmd99N\nl+3bWVi6NNdccw0rVqzwOpLIJaXlsdN3gU+Bima208wecs6dBnoAHwHrSH6K6LvARBUJLwULFmTK\nlCm8/vrrtGzZkieeeIJj6rcgIUxvKosEwf79+3n00UdZu3Yt7z/+OJWaN4ciRbyOJZmMpq4QCSMz\nZ85kW8eOPGJG5MSJ5GjRwutIkomoIIiEmR9//JE3mjen55dfku+OOyg4aRJERXkdSzIBFQSRMOSc\nY9KIEWR76ima5s3LZfPmEVGnjtexJMypIIiEse+//57RTZpwKF8+np0zh9KlS3sdScKYZjsVCWPl\ny5fn5W+/pczddxMTE8PMmTO9jiRZmM4QRELE6tWrue+++7juuusYOXIk+fPn9zqShBmdIYhkEjEx\nMXz55ZdERUVxzTXXsOWJJ2DvXq9jSRaiMwSREDR/3jw2t2lDJzNyT51K9mbNvI4kYUA3lUUyqZ9+\n+onXmjbliTVryNuqFfnHjIHISK9jSQhTQRDJxJxzvDN4MCUGDuTmokUpmJAAf/ub17EkRKkgiGQB\na9esYdodd3Cwbl1eHzuWKL3IJucRlJvKGdlP2Teex8wSzaxxeoOLZCXXVK/O85s3k61gQWrUqKE+\nzhIQqTpDyOh+ymb2PHAYWO+cW3ie9XWGIHIBs2fPpnv37jz55JP07t2biAg9LCjJgnKGkJH9lM2s\nAbAe+Jk/91YQkVRo0aIFiYmJzJ07l6a33cahIUPg9GmvY0km4M+vFuntpxwL1ALaAp38+H6RLKts\n2bIkJCRQ++qrWde/P79efz38+KPXsSTMBeRc0zk31Tn3BMnNdIaZ2Wh8/ZSdc/18n00Hxgbi+0Wy\nguzZs/PM0KEc+eADJm7ezKGKFTm9aJHXsSSMZfdj2z1Ayu7hV/jG/uCcSwASzrexc27KxXYeGxtL\ndHQ00dHRxMbGEhsb60dUkczr1oYNuXrzZp5v3Ji+zZoR+fDDXDZ0KOjeQqYXHx9PfHw827dvZ/v2\n7X7vL9WPnZpZNDDPOXe1bzkbsJHkm8o/AquANhnRQlM3lUXS7vTp0wx75hl+HzWKmrNmcXvDhl5H\nkiAL1mOn6qcsEuKyZcvGE0OGcOO8eXTo2JFnn32WpKQkr2NJGNGLaSKZ0L59+2jbti1JSUnExcVR\nokQJryNJEGi2UxH5i6JFi7J48WLq1atHTEwMn06fDvv2eR1LQpzOEEQyuSVLlvDxPffwLJBv/nwi\nbrzR60gSIDpDEJGLatCgAT3Xr2dQqVIcbNCA3wcNAv3CJeehMwSRLOLUqVO88sgj3DF5MmX+8Q8K\nzZ0LmiQvU9FspyKSJnNnzuRA+/aUbd2aehMmkDzVmGQGKggikmabNm3innvuISYmhlGjRpEnTx6v\nI0kG0D0EEUmzihUrsnLlSk6dOkWdOnXYunWr15EkBKggiGRRefPmZdq0aXTs2JHatWuzcOFCzZqa\nxemSkYiwYsUKurRsScKpUxT+v/8j4qabvI4k6aBLRiLit7p167J09WqGFC3KoVtv5djrr+vR1CxI\nZwgi8oeTJ0/yUocOtJ41i+KNG1Pg3Xchd26vY0kqhV1PZUv2opkNN7N26Q0uIhkvZ86cDJw2jS+G\nD2fZwoXsj4mBM2e8jiVBktpLRhOB21MO+Hoqv+UbrwK0MbNKKddxzm1wznUDWgO3+Yabktw74STJ\nXdZEJMTc37UrpVesoO2BAzw3YABnVBSyhKD3VAauAj5xzv0T6J7e4CISWNddfz3T1q5l+fLlNG3a\nlN9++83rSBJgXvRU3s3/iouecRMJYUWLFmXJkiWULVuWmjVrsmHDBq8jSQAFvacyMAdoaGbDuEB7\nTREJHTly5OCtt97iySef5KabbuLzF16An37yOpYEQNB7KjvnjgGdLrVz9VQWCS2dOnWicuXKLG3Y\nkEqvv07UkiXYddd5HStLU09lEfHU7t27GR4bS7/du4kcO5ac7fSwYKhQT2URCaorrriC57/5hpfq\n1WN/p04ceuwxPZqaSejFNBFJF+ccI557jqtfe42CU6dS/d57vY6U5WnqChHxhJnR88UXOTJrFrc9\n8ghxcXFeRxI/6QxBRPz2zTffcNddd3H//ffz/PPPExGh3zW9oAY5IhIS9u3bR/PmzSlevDhTpkxR\n0x0P6JKRiISEokWLsnTpUvLmzcu/qlblcK9eutkcZlQQRCTD5MqVi0mTJlG0TRs2jhrFgYYN4dgx\nr2NJKumSkYgExPtxcZx+6CFuKVuWQsuXQ7FiXkfK9HQPQURC1heJiSyvX5+O2bMT9d//YlWqeB0p\nU1NBEJGQtnPnTkbVrUvu2rV5Zto0cuTI4XWkTEsFQURC3qFDh2jdujWnT59m1qxZREVFeR0pU9JT\nRiIS8qKiovjwww/529/+Rt26ddm9W72xQpEKgogERfbs2Rk1ahQPPPAAtWvXZs2aNZCU5HUsSUGX\njEQk6GbPns37HTsyokwZCiUkQKFCXkfKFII12+l4M9trZl+fM97QzDaY2SYz63Oe7SqZ2dtmNtPM\nOvrGSpnZHDMbd75tRCTza9GiBY8sWMB7W7dyoGpVyIC5/MV/qb1kNBG4PeWAmUUAb/nGqwBtzKxS\nynWccxucc92A1sBtvuFqwGznXCeguh/ZRSSM1albl1u++ooRJ09yqFo13OrVXkfK8lJVEJxzK/hf\nH+SzagKbnXM7nHOngBlA03O3NbM7gQW+zwE+Abqa2RLg3+kNLiLhr2LFijy8fj3/KlaMw3XrkrRk\nideRsjR/biqXAnalWN7tG8PM2pnZG2ZWwjk3zznXGGjvW68D8KxzrgHQxI/vF5FMoGjRovRfs4YX\na9Sg06BBHDlyxOtIWZY/PZUvyDk3FZhqZjebWV8gEljm+/g/QH8zuw/YdqF9qKeySNaRN29eXkpI\noFu3bsTGxrJgwQKKaaqLS/Kyp3JZknsqV/Mt3wAMdM419C33BZxzbojfofSUkUiW5Jxj0KBBTJky\nhcWLF1OhQgWvI4UVf58ySssZgvn+nJUIVPAVih9JvnHcJr1BRETMjAEDBlCiRAluuukm5s+bR41q\n1UDTXQRFah87fRf4FKhoZjvN7CHn3GmgB/ARsA6Y4Zz7LnBRRSSr6NKlCyNHjuSdevX4pVYtOHzY\n60hZgl5ME5GQtXzZMrY3bkyTkiUptHIlFCnidaSQpsntRCRTW7tmDQk33ki7fPkouGoVlC7tdaSQ\npcntRCRTu6Z6de5Ys4ZRJ0/y29VX4zZs8DpSpqWCICIhr3z58nT49lteiIpi0BtvoCsIgaFLRiIS\nNg4cOEDjxo2pXLkyY8aMIVu2bF5HCim6ZCQiWUbBggX5+OOP2bFjB61bt+bkyZNeR8pUVBBEJKzk\ny5eP+fPnk5SURNOmTTl69KjXkTINFQQRCTuRkZHMmjWLyy+/nJmVK3Ns7FivI2UKKggiEpayZ8/O\n5MmT2X799Rx+9FGODhvmdaSwp4IgImErIiKCge+9x5g2bTj45JMcefllryOFtYDMdioiEixmxrMT\nJ/Jqnjy0HTAAd+oUl/Xv73WssKTHTkUkU3DOMbRXL64dP56/f/01xcuV8zpS0IV0T2XfOnnMLNHM\nGqdYnmRmo82sbXqDi4ikZGY88eabrHjqKWIbNeKHH37wOlLYCWhPZZ8+wMwUy82BWc65rsBdaU4s\nInIRzz33HO3btyc2NlZFIY1SdQ/BObfC1/cgpT96KgOY2dmeyn9MNGJmDYD1JHdMO+sK4OyZxul0\n5hYRuaC+ffsCcMstt7Bs2TJKlCjhcaLw4M9N5fP1VK4JyT2VgRpAFPAbyWcQR4GFvvXOFoV0X+sS\nEbmYvn37cvr0aerXq8eKQYMo1LKl15FCXkB7Kp9dNrMHgP2+xTnAW2Z2BzAvEN8vIgLw7LPPkuvw\nYY7cfz85vv+ey55+2utIIc2fgrAHKJNi+Qrf2F8456ak+PtRoMOldh4bG0t0dDTR0dHExsYSGxvr\nR1QRyar+OXgwbxw9SusBA7C8ecnXs6fXkTJMfHw88fHxbN++ne3bt/u9v1Q/dmpm0cA859zVvuVs\nwEagPsk9lVcBbTKijaYeOxWRjOScY0inTnSYOpV8b79Nno4dvY4UEMF67FQ9lUUkbJkZfcaNY3Tz\n5hzr1o3jc+d6HSkk6cU0Eckyzpw5w4C772bj4cNM+/e/yZkzp9eRMpR6KouIpEFSUhL33nsvuXLl\nYvr06ZmqyY4a5IiIpEH27NmJi4tj3759PProo2rHmYIKgohkOZGRkcydO5fExET6ayK8P+iSkYhk\nWT///DM31qnDlMqVqTl1KkRFeR3JL7qHICLih++3bOGTatVodNVVFFm1CnLk8DpSuukegoiIH8pX\nqEClJUv4ct06fr77bsjCv4yqIIhIllezTh2Spk1jz+LFHPjnP72O4xkVBBER4I6WLflywABODB/O\n0f/7P6/jeEL3EEREfJxzDGrdmnW//UbcggVh946C7iGIiGQQM+OZadP4+cQJns6CM6OqIIiIpJAj\nRw5mz57NnDlzmDJlyqU3yERSXRAyuK9yUzMbY2ZxZnarf/8EEZGMVbhwYT744AN69+7N2rVrvY4T\nNGk5Q8iwvsrOuQ+cc12AboDaGIlIyKlSpQrDhg3jvubNObJsmddxgiLVBcE5twI4cM7wH32VnXOn\ngLN9lf+Qoq/yz/y1ZWY/YGRaQ4uIBEPbtm1pW6MGpxs1wm3e7HWcgPO3hWZ6+iov8H0+GFjonFvj\nZwYRkYDpPW0ab1asyEOxsRTduhVy5fI6UsAEpKcyXLyvspn1ILnTWpSZVXDOjQlUDhERf+TKlYtW\n8fEkVqpETJcuFJ882etIAeNvQUhvX+URwIiL7Vg9lUUkVESXK8fqYcOI6NGDEw88QK769b2OBHjY\nUxmC11dZL6aJSKhxzvHqjTdy/fHj1PviC6/jnFfQXkxTX2URycrMjI4ffMD9P/7Iskz61JGmrhAR\nSYNFixbRrVs31q5dS/78+b2O8yfqhyAiEmSdO3cmT548DBs2zOsof6KCICISZPv376dy5cosW7aM\nKlWqeB3nD5rcTkQkyC6//HKee+453uzYEffTT17HyTA6QxARSYekpCSmFi/OLddcQ9mlS72OA+gM\nQUTEE9lRTB6VAAAHQ0lEQVSzZ6fchAlEJiRwYsUKr+NkCJ0hiIj44a0aNWhy8CDRW7d6HUVnCCIi\nXrp92jRy7NjBsc8+8zqK33SGICLip8l//zu1ChemkseXjnSGICLisehXXqHDvn2E+y+yKggiIn66\nqUkTDufKxccff+x1FL+oIIiI+MnMePzxx0PuzeW0SlVByMh+yhcaExEJZ23btiUxMZGNGzd6HSXd\nUnuGkGH9lC8yJiIStnLnzk2XLl0YMeKirV5CWqoKQkb2U75Ej2URkbDVvXt39k6axJEtW7yOki7+\ndExLbz/lWCDPOWMiImGvZMmSdM+Vi52TJ1P5hRe8jpNmAempfLF+ys65fueOiYhkFoevvJLL0vM+\ngnOwdSscOQInTkCZMlCsGFjwLqT4UxDS1U/5YmMpqaeyiISjyDp1yP3uu6lePykpiQULFvDJf//L\nk2PGcCxbNpKAYsePkwM4Fh3N6YQEChUt+pdtPeupHKx+yr59601lEQlL65Yto3SDBkQlJV30t/t9\n+/YxduxYRo8eTalSpbjrrrsoUKAA+fLlI0eOHOzcuZM9337LofXrmb1xI1dddRX169fn9ttv56ab\nbiJ79r/+Ph+UBjm+fsqxQGFgLzDAOTfRzBoBb5J8c3q8c25weoOc830qCCISlpKSkvgpZ04KfPUV\n+a655rzrHNy6lWqxsdx22210796dGjVqXHSfJ0+eZOXKlSxdupQFCxawa9cuBl53Hfft2EG+Hj2I\nuO8+yJdPHdNERELNkHLlqDd4MDVbtfrrh7t2cbBiRV5r04YXJ0xI1/63bNnC7Lg4do4bR7NffuHG\nM2fg3nvJM2WKCoKISCjp3r07V155Jb169frzB2fO8FO1aszYu5euO3eSO3duv77HOcfKlSt5b+hQ\n8s2bx4vHj/tVEALylJGISFYWExNDfHz8X8YPDBjAjo0buenTT/0uBpB8z6B27drUrl2bgwcP8mLB\ngn7tT3MZiYhksJiYGFavXv2nsTPff48NGcLqxx6jxvXXZ/h3FihQwO996JKRiEgGO3XqFAUKFGDv\n3r3ky5cPgM3NmvHJZ5/R7ocfyJYtW0C+V/0QRERCTI4cOahSpQpr1qz5Y6znvn3kf+mlgBWDjKCC\nICISAK2LF+fI6NEAfPHFF6zfs4c7H3zQ41QXp4IgIhIAVUqXpsR//gPA0KFD6dmz53lfJgsluocg\nIhIA3yxaRMk77+To1q1cU70627ZtI3/+/AH9Tr2YJiISgk6eOMGByEhGPPAAvxcqxNChQwP+nSoI\nIiIhaneuXCxKSqLBli2UK1cu4N/nb0EI7QtaIiJh7OsqVSj6++9BKQYZQWcIIiIBsn//frJly0ZB\nP98gTq2gvIdgZuPNbK+ZfX3OeEMz22Bmm8yszwW2zWNmiWbW2LdsZvaimQ33dVYTEcmULr/88qAV\ng4yQ2sdOJwK3pxwwswjgLd94FaCNmVU6z7Z9gJkplpuS3EznJMltNyWTOd8cLhIedOyytlQVBOfc\nCuDAOcM1gc3OuR3OuVPADJL/Z/8HM2sArAd+TjF8FfCJc+6fQPf0BpfQpf+phC8du6zNnxfTSgG7\nUizv9o1hZu3MbCjQBqgFtAU6p1jvbHE57cf3B1xG/nCkd19p2S41615qnQt9ntZxr2V0rnA4fmn9\nLFSPHYTf8cssP3sBeVPZOTfVOdfLOdfROfcEMB0Y6/t4DtDQzIYBCYH4/oyigpD+ca+F2/9QUruu\nCkJw95fVCkJaeiqXJbmncjXf8g3AQOdcQ99yX8A554b4HcpMjxiJiKRDsN5DMN+fsxKBCr5C8SPQ\nmuRLRH7z5x8kIiLpk9rHTt8FPgUqmtlOM3vIOXca6AF8BKwDZjjnvgtcVBERCaSQfDFNRESCT9Nf\ni4gI4GFBSOvbz2bW1MzGmFmcmd0a/MRyVjqOXSUze9vMZppZx+AnlpTSM/PAuTMOiHfS8fN3s5kt\n9/0M3nSxfXt5hpCmt5+dcx8457oA3YCWQc4qf5bWY7fBOdeN5AcPbgtyVvmr9Mw8cO6MA+KdtB4/\nBxwGcnGJ2SE8KwjpffsZ6AeMDEJEuYD0HDszuxNY4BsXD6X1+J0z44CeAPRYWo+fc265c+4OoC8w\n6GL7DrV7CBd8+xnAzAYDC51za87dUDx30WPnnJvnnGsMtA9yLkmdix2/WP4340Cn4MaSVLroz5/P\nQSDnxXYSNv0QzKwHUB+IMrMKzrkxXmeS1DGzm4HmQCSwzOM4kkbOuX4AZvYAsN/jOJJGZnY3yZeS\n8pN8WemCQq0g7AHKpFi+wjeGc24EMMKLUJIqFzt2CYT4NCVy4eN3lnNuSlATSVpc7OfvfeD91OzE\n60tGF3z72cxyknwT8kNPksml6NiFNx2/8BaQ4+flY6d6+zlM6diFNx2/8BbI46c3lUVEBPD+kpGI\niIQIFQQREQFUEERExEcFQUREABUEERHxUUEQERFABUFERHxUEEREBFBBEBERn/8HQ9wO8kFN+BwA\nAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f642d9f7350>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.loglog(np.abs(df.MTF11),df.muegamma,'k-')\n",
    "plt.loglog(np.abs(df.MTF11),df.muegamma_analy/4.,'r--')"
   ]
  },
  {
   "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}