Warning, /estarlight/production/Au_vs_Fe_test/Gold_Iron_Q2_test.ipynb is written in an unsupported language. File is not indexed.
0001 {
0002 "cells": [
0003 {
0004 "cell_type": "markdown",
0005 "metadata": {},
0006 "source": [
0007 "## Heavy ion comparison\n",
0008 "This notebook is motivated by the studies [_H. Mäntyasaari, et. al._](https://arxiv.org/abs/1712.02508) utilizing [_eSTARlight_](). The cross sections for coherent $\\rho$ vector-meson (photo)electro-production $\\sigma_{\\gamma +A \\rightarrow V.M. + A}(Q^2)$ are obtained for two heavy-ion targets (A=56,197) in a series of $Q^{2}$ intervals at eRHIC energies:\n",
0009 " - 18 GeV electron beam\n",
0010 " - 275 GeV/n ion beam\n",
0011 "\n",
0012 "The results are then plotted after taking into account the $A^{4/3}$ scaling.\n",
0013 "\n",
0014 "First step is to load relevant python modules and initialize some global variables (i.e. working directory)\n",
0015 "\n",
0016 "*Note*: Running the study as is set-up will take around an hour "
0017 ]
0018 },
0019 {
0020 "cell_type": "code",
0021 "execution_count": 4,
0022 "metadata": {
0023 "collapsed": true
0024 },
0025 "outputs": [],
0026 "source": [
0027 "import itertools\n",
0028 "import sys, string, os\n",
0029 "from math import log10, floor\n",
0030 "import matplotlib.pyplot as plt \n",
0031 "import seaborn as sns\n",
0032 "home_dir = '/Users/michaellomnitz/Documents/install_test/'\n",
0033 "unit_decode = { 'mb.':-3, 'microb.':-6, 'nanob.':-9, 'picob.':-12, 'femtob.':-15}"
0034 ]
0035 },
0036 {
0037 "cell_type": "markdown",
0038 "metadata": {},
0039 "source": [
0040 "Some simple funcitons used for utilies and run workflow:\n",
0041 "- round_sig: Round a number to a given number of significant figures\n",
0042 "- make_input_file: Takes template input file for _eSTARlight_ and sets the desired $Q^2$ range and target nucleii.\n",
0043 "- run_target_study: Iterates over target and $Q^2$ and calls _eSTARlight_ and saves the log file to desired location\n",
0044 "- decode_string: Takes relevant line from _eSTARlight_ and extract the cross-section and units (i.e. reads \"gamma+X --> VM+X 204.161 microb.\" and returns $204.161\\times10^{-6}$)\n",
0045 "- find_x_section: Iterates over produced log files and extracts cross-sections\n",
0046 "\n",
0047 "*Note*: It is assumed that home_dir contains eSTARlight and the compiled executable e_starlight"
0048 ]
0049 },
0050 {
0051 "cell_type": "code",
0052 "execution_count": 11,
0053 "metadata": {
0054 "collapsed": false
0055 },
0056 "outputs": [],
0057 "source": [
0058 "def round_sig(x, sig=2):\n",
0059 " return round(x, sig-int(floor(log10(abs(x))))-1)\n",
0060 "\n",
0061 "def make_input_file(nuclei_def, q2_range, template_loc):\n",
0062 " if len(nuclei_def) != 2 or len(q2_range) != 2:\n",
0063 " print('Something is awefully wrong')\n",
0064 " return\n",
0065 " #load template\n",
0066 " with open(template_loc,'r') as f:\n",
0067 " data = f.readlines()\n",
0068 " # -- Set prodiction id\n",
0069 " data[6] = 'BEAM_2_Z = '+str(nuclei_def[0])+'\\n' \n",
0070 " data[7] = 'BEAM_2_A = '+str(nuclei_def[1])+'\\n' \n",
0071 " data[37] = 'MIN_GAMMA_Q2 = '+str(q2_range[0])+'\\n' \n",
0072 " data[38] = 'MAX_GAMMA_Q2 = '+str(q2_range[1])+'\\n' \n",
0073 " with open(home_dir+'slight.in','w') as f:\n",
0074 " f.writelines(data)\n",
0075 " \n",
0076 "def run_target_study(pwd_to_template):\n",
0077 " print('Running ',len(nuclei), ' nuclei')\n",
0078 " for key in nuclei.keys():\n",
0079 " print( 'Starting work on ',key,' -- ',len(Q2_bins), 'Q2 bins in the study')\n",
0080 " for idx in range(len(Q2_bins)-1):\n",
0081 " this_bin = [Q2_bins[idx],Q2_bins[idx+1]]\n",
0082 " make_input_file(nuclei[key],this_bin, pwd_to_template)\n",
0083 " os.chdir( home_dir )\n",
0084 " os.system( home_dir+'e_starlight > log' )\n",
0085 " os.rename(home_dir+'log',home_dir+'eSTARlight/production/Au_vs_Fe_test/'+key+'_Q2_range_'+str(this_bin[0])+'_'\n",
0086 " +str(this_bin[1])+'.txt')\n",
0087 " print(idx,'/',len(Q2_bins),' done for ',key)\n",
0088 " \n",
0089 "def decode_string( line ):\n",
0090 " space_pos = [pos for pos,char in enumerate(line) if char == ' ']\n",
0091 " n_space = len(space_pos)\n",
0092 " num = float(line[space_pos[n_space-2]:space_pos[n_space-1]])\n",
0093 " #num = round_sig(num,2)\n",
0094 " units = line[space_pos[n_space-1]+1:-1]\n",
0095 " return num,units\n",
0096 " \n",
0097 "def find_x_section(nucleii):\n",
0098 " x_sec_to_ret = []\n",
0099 " for i in range(len(Q2_bins)-1): \n",
0100 " file_loc = home_dir+'eSTARlight/production/Au_vs_Fe_test/'+nucleii+'_Q2_range_'+str(Q2_bins[i])+'_'+str(Q2_bins[i+1])+'.txt'\n",
0101 " with open(file_loc,'r') as f:\n",
0102 " data = f.readlines() \n",
0103 " x_sec, units = decode_string(data[47])\n",
0104 " x_sec_to_ret.append(x_sec*10**unit_decode[units])\n",
0105 " #print (x_sec,'10**',unit_decode[units])\n",
0106 " #return x_sec, units\n",
0107 " return x_sec_to_ret"
0108 ]
0109 },
0110 {
0111 "cell_type": "markdown",
0112 "metadata": {},
0113 "source": [
0114 "### Run the study\n",
0115 "This procedure may take some time. No events are generated but the look-up tables for _eSTARlight_ have to be built for each target and $Q^2$ interval\n",
0116 "We define variables for this study:\n",
0117 "- $Q^2$ intervals\n",
0118 "- Nucleii definition: atomic number and nucleon number"
0119 ]
0120 },
0121 {
0122 "cell_type": "code",
0123 "execution_count": 15,
0124 "metadata": {
0125 "collapsed": true
0126 },
0127 "outputs": [],
0128 "source": [
0129 "#Q2_bins = [round((i+1)*0.01,2) for i in range(10)]\n",
0130 "Q2_bins = [round((i+1)*0.01,2) for i in range(9)] + [round((i+1)*0.1,1) for i in range(9)]\n",
0131 "Q2_bins+= [(i+1)*1 for i in range(9)] + [(i+1)*10 for i in range(9)] + [100,200]\n",
0132 "nuclei = {'Au':[79,197], 'Fe':[26,56]}"
0133 ]
0134 },
0135 {
0136 "cell_type": "code",
0137 "execution_count": 16,
0138 "metadata": {
0139 "collapsed": false
0140 },
0141 "outputs": [
0142 {
0143 "name": "stdout",
0144 "output_type": "stream",
0145 "text": [
0146 "Running 2 nuclei\n",
0147 "Starting work on Au -- 38 Q2 bins in the study\n",
0148 "0 / 38 done for Au\n",
0149 "1 / 38 done for Au\n",
0150 "2 / 38 done for Au\n",
0151 "3 / 38 done for Au\n",
0152 "4 / 38 done for Au\n",
0153 "5 / 38 done for Au\n",
0154 "6 / 38 done for Au\n",
0155 "7 / 38 done for Au\n",
0156 "8 / 38 done for Au\n",
0157 "9 / 38 done for Au\n",
0158 "10 / 38 done for Au\n",
0159 "11 / 38 done for Au\n",
0160 "12 / 38 done for Au\n",
0161 "13 / 38 done for Au\n",
0162 "14 / 38 done for Au\n",
0163 "15 / 38 done for Au\n",
0164 "16 / 38 done for Au\n",
0165 "17 / 38 done for Au\n",
0166 "18 / 38 done for Au\n",
0167 "19 / 38 done for Au\n",
0168 "20 / 38 done for Au\n",
0169 "21 / 38 done for Au\n",
0170 "22 / 38 done for Au\n",
0171 "23 / 38 done for Au\n",
0172 "24 / 38 done for Au\n",
0173 "25 / 38 done for Au\n",
0174 "26 / 38 done for Au\n",
0175 "27 / 38 done for Au\n",
0176 "28 / 38 done for Au\n",
0177 "29 / 38 done for Au\n",
0178 "30 / 38 done for Au\n",
0179 "31 / 38 done for Au\n",
0180 "32 / 38 done for Au\n",
0181 "33 / 38 done for Au\n",
0182 "34 / 38 done for Au\n",
0183 "35 / 38 done for Au\n",
0184 "36 / 38 done for Au\n",
0185 "Starting work on Fe -- 38 Q2 bins in the study\n",
0186 "0 / 38 done for Fe\n",
0187 "1 / 38 done for Fe\n",
0188 "2 / 38 done for Fe\n",
0189 "3 / 38 done for Fe\n",
0190 "4 / 38 done for Fe\n",
0191 "5 / 38 done for Fe\n",
0192 "6 / 38 done for Fe\n",
0193 "7 / 38 done for Fe\n",
0194 "8 / 38 done for Fe\n",
0195 "9 / 38 done for Fe\n",
0196 "10 / 38 done for Fe\n",
0197 "11 / 38 done for Fe\n",
0198 "12 / 38 done for Fe\n",
0199 "13 / 38 done for Fe\n",
0200 "14 / 38 done for Fe\n",
0201 "15 / 38 done for Fe\n",
0202 "16 / 38 done for Fe\n",
0203 "17 / 38 done for Fe\n",
0204 "18 / 38 done for Fe\n",
0205 "19 / 38 done for Fe\n",
0206 "20 / 38 done for Fe\n",
0207 "21 / 38 done for Fe\n",
0208 "22 / 38 done for Fe\n",
0209 "23 / 38 done for Fe\n",
0210 "24 / 38 done for Fe\n",
0211 "25 / 38 done for Fe\n",
0212 "26 / 38 done for Fe\n",
0213 "27 / 38 done for Fe\n",
0214 "28 / 38 done for Fe\n",
0215 "29 / 38 done for Fe\n",
0216 "30 / 38 done for Fe\n",
0217 "31 / 38 done for Fe\n",
0218 "32 / 38 done for Fe\n",
0219 "33 / 38 done for Fe\n",
0220 "34 / 38 done for Fe\n",
0221 "35 / 38 done for Fe\n",
0222 "36 / 38 done for Fe\n"
0223 ]
0224 }
0225 ],
0226 "source": [
0227 "#\n",
0228 "run_target_study(home_dir+'eSTARlight/production/templates/slight_Q2test_template.in')"
0229 ]
0230 },
0231 {
0232 "cell_type": "markdown",
0233 "metadata": {},
0234 "source": [
0235 "### Note\n",
0236 "Next cell is not always needed. In this case the study was done in parts and the array holding the full $Q^2$ range needed to be built "
0237 ]
0238 },
0239 {
0240 "cell_type": "code",
0241 "execution_count": null,
0242 "metadata": {
0243 "collapsed": false
0244 },
0245 "outputs": [],
0246 "source": [
0247 "temp = [ (i+1)*1. for i in range(1,100)]\n",
0248 "Q2_bins = [round((i+1)*0.01,2) for i in range(10-1)] + [round((i+1)*0.1,1) for i in range(10)] + [ (i+1)*1 for i in range(1,100)]\n",
0249 "print(Q2_bins)"
0250 ]
0251 },
0252 {
0253 "cell_type": "markdown",
0254 "metadata": {},
0255 "source": [
0256 "### Load the points"
0257 ]
0258 },
0259 {
0260 "cell_type": "code",
0261 "execution_count": 17,
0262 "metadata": {
0263 "collapsed": false
0264 },
0265 "outputs": [],
0266 "source": [
0267 "au_points = find_x_section('Au')\n",
0268 "fe_points = find_x_section('Fe')"
0269 ]
0270 },
0271 {
0272 "cell_type": "markdown",
0273 "metadata": {},
0274 "source": [
0275 "### Further curating the data\n",
0276 "Some minor reformatting before plotting the results. As mentioned, we are interested in $A_1^{-4/3}\\sigma_{A_1}$ / $A_2^{-4/3}\\sigma_{A_2}$ so the data points for Gold(Au) and Iron(Fe) are multiplied by their respective scaling factors. Secondly, the average $<Q^2>$ is obtained for each bin assuming a funcitonal form $Q^{-4}$ for the bin centers of each interval."
0277 ]
0278 },
0279 {
0280 "cell_type": "code",
0281 "execution_count": 18,
0282 "metadata": {
0283 "collapsed": false
0284 },
0285 "outputs": [],
0286 "source": [
0287 "import math as m\n",
0288 "au_scale = 197**(-4/3)\n",
0289 "fe_scale = 56**(-4/3)\n",
0290 "x_bins = []\n",
0291 "for i in range(len(Q2_bins)-1):\n",
0292 " num = m.log(Q2_bins[i+1]) - m.log(Q2_bins[i])\n",
0293 " denom = (Q2_bins[i]**(-1)-Q2_bins[i+1]**(-1))\n",
0294 " x_bins.append(num/denom)\n",
0295 "y_val = [ (au_scale*au_points[i])/(fe_scale*fe_points[i]) for i in range(len(au_points)) ]\n"
0296 ]
0297 },
0298 {
0299 "cell_type": "markdown",
0300 "metadata": {},
0301 "source": [
0302 "## Plot the figure\n",
0303 "Set axis labels, figure title and curve legend, then plot and save the figure"
0304 ]
0305 },
0306 {
0307 "cell_type": "code",
0308 "execution_count": 19,
0309 "metadata": {
0310 "collapsed": false
0311 },
0312 "outputs": [
0313 {
0314 "data": {
0315 "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf8AAAGOCAYAAACZjwraAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XlcVOX+B/DPrCzDsO8iq+KGC6JlroWZZi5dtdzimlpm\nmXYrS8ubdW/erJvV9ZplWlo/81pmplGWXbVuZlouiIKggIAi+w4DzHp+f6AoscgycGaYz/v18gVz\nzsw53+EBP/M855znSARBEEBEREQ2Qyp2AURERNS5GP5EREQ2huFPRERkYxj+RERENobhT0REZGMY\n/kRERDaG4U9ERGRjGP5EREQ2huFPRACA4uJiDBs2DBEREdDpdGKXYzZd9X0RtQfDn6iLWblyJTZs\n2NDq161btw6zZ8+GXq9HSkpKB1TWPl31fRGJgeFPZIGmTJmCSZMmobi4uN7ydevWYcqUKTAajWbd\nX1xcHM6cOYMlS5bAy8sL58+fb/U2Lly40KLnpaSkYMGCBRg0aBAGDhyI559/vsN65OZ4X0RdEcOf\nyAK98847yM3NxU8//VS3LCMjA5988gleeuklyGQys+3LaDTib3/7G55//nnI5XL06tWr1SGp0Wgw\nf/58vP76680+78yZM3jooYcwYsQIfPXVV3j11Vfx3Xff4ZNPPmnPW2iUOd4XUVfF8CeyQGFhYYiK\nikJqamrdsn/84x8YN24chg4datZ97dixAx4eHrjzzjsBAOHh4a0OSZVKha1bt2Lv3r34+9//jsbu\nF2YymfDiiy9i4cKFWLhwIUJCQjBlyhQMGzYMp06dMsdbqccc74uoq5KLXQARNS40NBSXLl0CABw6\ndAgnT57Ed9991+B5mzZtwgcffFD3WKfTQSKRYOvWrXXLHnvsMSxevLjBawsKCrBx40Z8+umndcvC\nw8Oxc+dOmEwmSKU3+geZmZm45557bln3jh07MHToUNx77731lsfFxeHy5cuYPXt2veVyuRwKhUK0\n90Vkixj+RBYqJCQEBw8ehE6nw9q1a/HYY4/B19e3wfNmzZpVL2jXrVsHHx8fxMTE1C1zcXFpdB+v\nv/46SktLMXXq1LplgiDAZDLh0qVL6NGjR91yf39/7N+/v8l68/Pz8fTTT6N///6Ijo5usP7YsWMI\nCgqCWq2ut6/ExEQsWLBAtPdFZIsY/kQWKjQ0FFevXsXGjRshkUgaDUgAcHV1haura91jlUoFFxcX\nBAUFNbv93377DYcPH8auXbvg4OBQt9xgMGDatGk4f/58vZBUKBQICwtrdFsajQZPPPEEIiMjsX79\neiiVygbPSUhIaHBi34EDB1BaWtpglKAz3xeRLWL4E1mokJAQGI1GbN68Ge+9916jgdpWer0ef//7\n3/HnP/8ZAwcObLDe398f58+fx5QpU1q0PZVKheXLl+POO+9sdAgfAM6fP4/Kykps3boV0dHROHXq\nFF5//XWsWLECPj4+7Xo/15n7fRF1VQx/Igvl4eEBFxcXDBw4EHfddZdZt/3xxx+joKAACxcubHR9\nSEgIEhMTW7XNcePGNbmuqKgIeXl52LhxI7Zs2YK3334b3bp1w4oVKzBjxoxW7ac5HfG+iLoiidDY\nablEJDqTyYSoqCi89tprjQ6LW5Off/4ZTzzxBOLi4pocGSCizsNTXoksVEZGBqqqqtCvXz+xS2m3\n68fZGfxEloHhT2ShkpKSoFar0b17d7FLabfz58+jT58+YpdBRNdw2J+IiMjGsOdPRERkYxj+RERE\nNobhT0REZGMY/kRERDaG4U9ERGRjGP5EREQ2huFPRERkYxj+RERENobhT0REZGMY/kRERDaG4U9E\nRGRjGP5EREQ2huFPRERkYxj+RERENobhT0REZGMY/kRERDaG4U9ERGRjGP5EREQ2huFPRERkYxj+\nRERENsYiwj8+Ph4xMTGNrquursasWbOQlpZWb3lRURHGjBnTYDkRERE1Ty52AVu2bMHXX38NBweH\nBuvOnTuHl19+GXl5efWW6/V6rF69Gvb29p1VJhERUZches8/MDAQGzZsaHSdTqfDxo0bERoaWm/5\nG2+8gVmzZsHb27szSiQiIupSRA//8ePHQy5vfAAiKioKfn5+9Zbt2bMH7u7uGDVqVKv2YzAY21wj\nERFRVyL6sH9rffnll5BIJDh27BiSkpKwYsUKvP/++/Dy8mr2dSUlVY0u9/JSo6CgoiNKpTZim1gm\ntovlYZtYJktpFy8vdZPrrC78d+zYUfd9TEwMXnnllVsGPxEREd0g+rD/H8XGxuLzzz8XuwwiIqIu\nSyIIgiB2EZ2hqSEYSxmeoRvYJpaJ7WJ52CaWyVLapblhf4vr+RMREVHHYvgTERHZGIY/ERGRjWH4\nExER2RiGPxERkY1h+BMREdkYhj8REZGNYfhbiB07PsHUqeOh1Wpb9PzExAQ8+eSiusc6nQ6vvLIK\nixY9jKefXoIrVy5j//5YPPnkIjz55CIsWvQwoqOHo6Ki9trTsrJS/POf/2j1/ouLi/D222/UW1Za\nWorXXvtbm97HggVz62q8vo0PPvgAjz02HwsWPIRvvtlb99zrNefkZOOee8bUve7JJxdh27YtLdpf\nS2i1WsyYMdks2zp+/Ffs27cHWq0WsbF7b/2CDmQymfDmm6/hscfm48knFyEr60qjz/vj75bRaMRr\nr/0Ns2bNwuOPL8SlS6kwGAx45ZVVWLx4AZ544hFkZmZ00rsgInOwuul9u6offvgOY8feg0OHfsDE\nic0Hz44dn+DAgf2wt79xG+TY2K/g4OCIzZs/xuXLGXjnnX/i7bffrdvWW2+9gfvumwK1unbShy1b\n3se0aQ+2ev/u7h5wdFQhLu4UIiOjAADHjv2C228f3ur3odVqIQgC3n13c92y06dPIi4uDu+//xFq\namqwc+f2unU31xwcHFLvdZZq2LDan0tOTjZiY/di8uT7W/S6PXu+wFdffYGqqirMnv0QZsyY1e5a\njhz5CTqdDh98sA0JCefw7rvv4PXX3673nMZ+t44ePQIA+Oyzz3DgwI/YvPk93HffFBiNRmzatBUn\nThzH5s0b8Y9/vNnuGumGvJIqfHM0A1cKKqE3mOr9c1YpEeSrRt9QT3g6KRDoq4bKXiF2yWRFGP7X\n7DqcihPJ+Wbd5tDe3ngwusctn3f69En4+wfg/vun4+9/X42JEyfjxx8P4ssvd9V73hNPLEPfvhHo\n1i0A//jHm3j11dV169LT0+uCJjAwGBkZ6XXrkpPPIz09Dc8+uwIAoNFUIinpPJYv79nk/gFg//5Y\nZGZm4PHHl0Kr1WLu3BnYvTsW48ZNwEcffVAX/idO/IZnnlnR5Haaei+CANTU1ODpp5fAaDRi0aIl\n+P334wgPD8eLLy6HRqPBkiVPNag5Jye70Z+jwWDAm2++hqysKzCZTHj00ccxePCQW/78AaCqqgp/\n//tfUVFRgW7dAprd3v79sTh27Ci02hpcvZqFuXPnYeLEybh8ORNr1/4NMpkcJpMJL7+8BqdOnUBm\nZgbKy8uQkZGObdu2IDMzA/fccy+GDx+JjIx0bNz4L7z55vq6Wn766RBOnPgN27b9B6WlpZg3bybu\nv39Gk3e/bO5n3LdvRN3js2fP4Pbb7wAARET0R3JyUoPtNPa7NXr0nRg+fCQAIC8vF05OanTvHgSj\n0QiTyQSNRtNsbdQ6JRVaxP6agSPx2TCaBCgVUijlMigVUtjbyaF2lKCoXIvfzufht/N5da/zdLFH\nkK8aQT7quq/OKqWI74QsGf9iLcA33+zD5Mn3IzAwGAqFAomJCbjrrrtx1113N/r8O+8c2yAAe/YM\nx6+/HsHo0XciMTEBhYUFMBqNkMlk+L//24YFC24M4yYmJiAwMKjZ/ffrF4GmBAeH4OzZMwBqA7Km\npgZOTk5Nbqep95KWlorZs2MwefL9uHLlMpYvX4ZBgwYjPb0Ya9asQ07OVaxY8Qz+858vG9SckZFe\nb2j65ZfX4JdffoaLiyteeGE1yspKsWTJInz66a4G+23M3r1fIiQkDI89tgSJiQk4ffokYmP3Nrk9\njaYSb7/9Lq5cuYwVK57GxImTceLEb+jTpx+eeOIpxMfHQaOprNv+n/+8AGlpqZg//1GcPn0SX321\nG8OHj8S3336NSZOm1qtl9+7P8fzzqyCXy+Hp6Qm5XA5BEFBdXY233nodCoUCkZFRuOeee+te09zv\ny3UajQYqlVPdY6lUCoPBUC+4G/vdAgC5XI4VK1bghx/+izVr3oCDgwNyc7MxZ86Ma4dj3mnRz5ma\npqnRY//xTBw6mQWdwQRfd0dMGx2KqF5ekEgk9Z4rCAIKympQWmXAuZR8ZOZWICO3AqcuFODUhYK6\n57mp7RDko0agjxOCfNUI9nWGq5OywfbI9jD8r3kwukeLeunmVl5ejmPHjqKkpBi7d38OjaYSe/Z8\njvz83Fv25G52331TkJmZjieeeAT9+w9Er169IZPJUFFRgcuXM+v1gEtLS+Hu7t7s/huG/41bQMhk\nMsjltb3b+Pg4DBgwsNntNNUr7dEjHAEBAZBIJAgMDIKLiwskEglGjhwJhUKBwMBgKJV2KC0tqVcz\n0Piwf1paKs6ejcP58wkAAKPRgNLSUri6utY9Z/Pm9+o+uKxf/z5kMhkA4MqVyxg+fAQAoF+/CMjl\n8ia3BwA9eoQDALy9faDT6QAAkyZNxY4dn+DZZ5dCpXLCY48tabStIiOj8M47/0RJSQl+//14vecZ\nDAakpl6s+6BTWFgIZ2cXKBQKHDr0A+68cyxGjhyN1atfqBf+Len5q1QqVFXduLW1IAit6rG/8cYb\nmD9/MRYtehhjxtyF2267A4sXP4m8vFw89dTj+OSTz2BnZ9fi7VGtco0Oh05l4dCpLFRpDXBT22HO\nyBCM6O8LmbTx07IkEgm8XR3Qr6ca4f61h/IEQUBJhRaZeRXIzK3A5bxKZOZV4ExqIc6kFta9Vu2o\nQLCvM6aODEGov3OnvEeyPAx/kf3ww35MmjS1bni7pqYGDzwwBU8++cwte3I3S04+j6io27Bs2bNI\nTj6PvLwcAEB8/GkMGTK03nPd3NzqTvxrav8lJSVQKpUoKqr9T+PCheS61wuCAJlMBqlUil9/PYKp\nU6c3u52meqVffbUbaWmpWL58JQoLC6DRaDB69F3Yt+8LTJo0A0VFhaipqYazs0u9mpsSFBQMb29v\n/PnPC6DV1uCTT7bC2bn+f26LFj3R6GtDQkKQkHAOo0bdiYsXk2EwGJrdXmM9p19++R8GDozEggWL\n8N//fo8dOz7BoEGDrz1fCkEw1b12/PiJ+Ne/3sRttw2rF8Dp6Zeg0Whw9WoW/Pz88cEH79Yd7y8o\nyEdYWO0HVOkfQqElPf/+/Qfi6NEjGDt2HBISziE0tGUfdr///lsUFOTjmWeWwd7eHlKpFM7OLpDJ\naut2dnaBwWCAyWRq0faoVl5xFQ78fhlHE3KhN5jg5KDAg3f1QPTgblAqZK3enkQigbuzPdyd7RHZ\n88Ztzss0OlzOqx0ZuJxbgcy8Cpy7VIQLV0qwdPoA9At2b2ar1FXxbH+Rxcbuw/jxE+se29vbY8yY\naMTGftWq7QQEBGLXrp147LH52LJlE5YufQYAcPlyJvz9u9V7br9+/ZGamnLL/d9++3Dk5ubg8ccX\n4vDhg1CpVABqe9gREf0B1PaYAwOD2vQ+Jk2aisrKCjz++EKsXv0CXnhhNUaMGIU+ffrg0Ufn4fnn\nn8Yzz6yATCarV3NTpk6dhszMDDz55CIsXrwAvr5+DUKy6ddOR3b2VTz++ELs2fMFFApFq7fXu3df\nfPjhJixbthj79u3B9Okz69a5ublBrzfgvff+DQCYOHEy/ve/ww2G/FNSLuCeeybglVdWYd68WfDx\n8cXUqdMAAF5e3sjPrz0v5foHidYYPfouKJVKLF68ABs2vI1ly55BeXkZXnzxuWZfN2ZMNC5evIC5\nc+fimWeWYtmyZ/Dgg3Nw8WIynnjiESxbthiLFi2Bg4NDs9uhWqlXy/DunnN4cfNx/HQmG65OSswd\nF443Hx+OCbcHtin4m+OiUqJ/qAcmDw/Gkmn98c/Hh2PptP4wmQSs/yK+3mECsh28pa+F3Hqxs735\n5muYOnUawsN7t/q17723HiNGjMHAgYM6oLKm26Q9NVuagoJ8rFnzMtavf7/e8vXr30JERH+MHXtP\ng9dUV1fjnXf+CaVSiQEDBtUb9u8Mtvq3Yi75JVXYfuACEjNKAADBvmrcOywIUeFekErbdgy+PW1y\nPqMYG748B73BhPkTe2NEf782bYcaspS/Fd7Slxp45JHF+Oqr3a1+XVFRITQaTYcFf3PaWrOl+d//\nDuPZZ5di4cLHGqxLSbmAnj3DG32dg4MDXnzxZSxf/kKnBz+1ndFkwnfHM/HSR78jMaMEfYPd8Pzs\nSLw0bwiG9vZuc/C3V99gdyyfPQgOdjJ89G0SDp5sfN4H6prY87eQT2h0A9vEMrFdWi8ztwLbvkvC\n5bxKODsqMGdcOIb29jbb2fbmaJOs/Eq89fkZlGl0+NOoEEwaHsyrAdrJUv5Wmuv584Q/IiIz0+qN\n2HckHQdOXIYgACP7++HB6B5wcrC8iXgCvJ2w8qHBWLfzDL46ko4qrQEP3tWDHwC6OIY/EZEZJaYX\n45Pvk1FYVgMvV3vMm9AbfS38jHofN0e88NBgvPX5GRz4/QqqagyYN6G3aIckqOMx/ImIzMBgNOHL\n/6XhwO9XIJVIcO/tgZgyMgR2Zj57v6O4O9tjxdzBeOfzeBw5m4NqnRGLJveFXMZTw7oihj8RUTuV\nVGjx/r4EpGaVwcfdEYun9EOQb9PHWy2Vs6MSz82OxL93x+Nkcj5qdAYs+VN/q/kAQy3Hj3RERO1w\nPqMYr2z7HalZZRja2xur5w2xyuC/ztFejqdnDkL/UA8kXCrG25+fQVWNQeyyyMwY/kREbWASBMQe\nTcdbn9WG49xx4Vg8tR8c7Kx/QNVOIcPS6f1xWx9vpGSV4Z87T6O8Sid2WWRG1v9bSkTUySqr9dgc\nm4iES8Vwd7bD4/dHIMzfReyyzEouk2LR5H6wV8rxc3w2Xv/0NJbPGgR3Z3uxSyMzYM+fiKgV0rLL\n8Mq235FwqRgRoe54Zf5tXS74r5NKJZg3oRcm3BaI3OIqrP30FPKKq279QrJ4DH8iohYQBAEHT17B\n65+eRkmFFn8aHYq/PDDQIq/dNyeJRIIH7grDtNGhKCrXYu2O07iSX3nrF5JFY/gTEd2CSRCw81AK\n/nMwBSp7OZ6dOQiThwdDaiMT4UgkEkwaHoy548JRrtHhjR2nkXq1TOyyqB0Y/kREzTAYTfjomyQc\nPJmFbp4qrH54qMVP2tNRxkYF4JFJfVCjM2LdZ3FIzCgWuyRqI4Y/EVETdHoj3vsqAccScxHq74wV\ncwfb/AlvwyP8sORPEbwlsJVj+BMRNaKqxoB3dsXjTGoh+gW7YfmsQV3++H5LRYZ74S8PDIRMKsX7\nexNw9FyO2CVRKzH8iYj+oFyjwz93nsaFK6UY0ssLy2YMhL2SV0bfjLcEtm4WEf7x8fGIiYlpdF11\ndTVmzZqFtLQ0AIBer8dzzz2HOXPmYMaMGTh06FBnlkpEXVxRWQ3W7jiNy3mVGD3QH4unRkAht4j/\nKi1OmL8LVswZDBeVEv85mILYo+mwkbvEWz3Rf6O3bNmCv/71r9BqtQ3WnTt3DnPnzsWVKzc+UX79\n9ddwdXXFf/7zH3z44Yd49dVXO7NcIurCsgs1eO3atewThwVh3oRevLPdLVy/JbCHsz2+OpKOXT+m\n8gOAFRA9/AMDA7Fhw4ZG1+l0OmzcuBGhoaF1yyZMmICnnnoKQO11tzIZbzhBRO2XnlOO13fUXsP/\nwF1hmHFnGO9p30LXbwns5+GIA79fwcffJcNk4gcASyb6Qazx48cjKyur0XVRUVENlqlUKgBAZWUl\nli1bhr/85S8t2o+bmyPk8sY/KHh5We9NOLoqtoll6qrtcja1AOs+i4NWZ8STDwzC+GFBYpfUYpbS\nJl5eary5bDRe3nIMR87mwCSR4Lm5UZDZ6C2BLaVdmiJ6+LdFTk4OlixZgjlz5mDy5Mktek1JSeNT\nUnp5qVFQUGHO8qid2CaWqau2y+mLBdi0LwEAsHhqBAaHuVvN+7TENnl6xkD8e3c8jsZnI9jbCWOj\nAsQuqdNZSrs09wHE6j6SFRYWYsGCBXjuuecwY8YMscshIit25Gw2Nn51DjKpFE89MBBDenuLXZLV\nc7SX4/E/9YdSIcW3xzKg0xvFLokaYXHhHxsbi88//7zJ9Zs2bUJ5eTnee+89xMTEICYmBjU1NZ1Y\nIRFZO0EQsP94JrbtT4ajnRzLZw9CPxudta8juKiUGDs4AKWVOvx0JlvscqgREsFGTstsagjGUoZn\n6Aa2iWXqKu1iEgTsOpyKH05cgbuzHZ55cBD8PVVil9UmltwmFVU6PL/pGOzkUryxeDjslLZzcral\ntEuXGvYnImorg9GED785jx9OXIGfhyNefCjKaoPf0qkdlRg3pDvKq/Q4HNf4Sd0kHoY/EdkErc6I\nf+8+i+OJeQjr5owXHoqy+Xn6O9r427rDwU6O745fRrXWIHY5dBOGPxF1eZXVerz5WRwS0osxIMwD\ny2dFcp7+TqCyV2D80O6orNbj0Cn2/i0Jw5+IurSishqs/fQULmWXY3iEL56c1h92Cts5/iy2cUO7\nQ2Uvx4HfL6Oqhr1/S8HwJ6Iu62pBJV779BRyiqow4bZALLivD+Q2OumMWBzs5JhweyA0NQb8cOKy\n2OXQNfwrIKIuKTWrrG663gfv6oEHo3tAyul6RTE2KgBqRwX+e/IKKqv1YpdDYPgTURcUn1qIdZ/F\noVprxML7+mDC7YFil2TT7JVy3Ht7EKq1Rhz4nb1/S8DwJ6Iu5ei5HGz48hwAYOn0/hjR30/kiggA\n7hrcDS4qJQ6ezEJ5lU7scmwew5+IuozvfsvER98mwcFOhuWzIjGwh6fYJdE1dgoZJt4RBK3eiO9/\nY+9fbAx/IrJ612ft++LHNLip7bBy7mD0CHARuyz6gzsH+cNNbYfDp7JQVqkVuxybxvAnIqtmMJqw\n9dskfP/75bpZ+7p5OYldFjVCIZdh0vBg6AwmfHs8U+xybBrDn4isVrXWgA1fnsOvCbkI9XfGyrmD\n4eHCWfss2agBfvBwtsdPcdkoqWDvXywMfyKySvml1fjH9lM4d6kIEaHueG5WJNSOSrHLoluQy6SY\nPCIYBqMJ3xzLELscm8XwJyKrk5xZglc/PoHsQg3GDemOp2YMsKm7xlm74RG+8HZ1wM9nslFYVi12\nOTaJ4U9EVuWnM1fx1udnUKMz4uF7e2P23T0hk/K/MmtyvfdvNAn45lce+xcD/2KIyCoYTSbs+OEi\n/u/7C3Cwk2P5rEEYPdBf7LKoje7o5wtfd0ccPZeD/FL2/jsbw5+ILJ6mRo93dsXj0OksdPNS4aV5\nQ9Ar0E3ssqgdpFIJpo4MgdEkIPaXdLHLsTkMfyKyaDlFGqz55CTOZ5RgUA9PvPhQFLxcHcQui8xg\naB9vdPNU4dfEXOQWV4ldjk1h+BORxUq4VIQ1/3cKeSXVmDgsCE9O7w8HO7nYZZGZSCW1vX9BAL5m\n779TMfyJyOIIgoD/nriCd76Ih95gwqOT+mLGnWG8K18XNLiXFwK9nfDb+TxcLdSIXY7NYPgTkUUx\nGE345Ptk7DyUAmdHJVbMjcQdEb5il0UdRCqR4P5RoRAA7GPvv9Nw/IyILEZ5lQ7v7TmHi1llCPJR\nY+n0/nB35ox9Xd3AHh4I8VPjZHI+LudVINBHLXZJXR57/kRkEbLyK7Hmk5O4mFWGIb29sfKhwQx+\nGyG51vsH2PvvLOz5E5Ho4lIKsDn2PLQ6I+4fGYLJI4Ih4fF9mxIR4o4e3VwQl1KIjNxyBPs6i11S\nl8aePxGJRhAEfHssA+9+eQ6CScAT90dgysgQBr8Nqu39hwAA9h5h77+jsedPRKLQG4zY9l0yjifm\nwU1th2XTByDIl8d6bVmfIDf06u6Ks2lFSLtahrBuLmKX1GWx509Ena60UovXd8TheGIewvydsXre\nEAY/QSKR4E+ja4/9f3uMc/53JPb8iahTZeSWY8OX51BSocXwCF/Mm9ALCjnvyEe1wru7IthXjfi0\nQpRUaOGmthO7pC6JPX8i6jQnkvPx+qenUVqhxQN3hWHhfX0Y/NTAqAF+EATgWGKu2KV0WQx/Iupw\nJkHA3iOX8P7eBEikEiydMQD33h7EE/uoUbf39YFCLsWRszkQBEHscrokiwj/+Ph4xMTENLquuroa\ns2bNQlpaGgDAZDJh9erVmDlzJmJiYpCZyeNCRJZMqzPi/b0J+PpoBjxd7LEqJgqDeniKXRZZMEd7\nBaLCvZBXXIXUq2Vil9MliR7+W7ZswV//+ldotdoG686dO4e5c+fiypUrdcsOHjwInU6Hzz//HM8+\n+yxef/31ziyXiFqhqKwGaz89hVMXCtCruytemjcEAV5OYpdFVmDkAD8AwJGzOSJX0jWJHv6BgYHY\nsGFDo+t0Oh02btyI0NDQumWnTp3CqFGjAACDBg1CQkJCp9RJRC0nCAKOnsvB6q2/43J+JUYP9Mez\nswZB7agUuzSyEr2D3ODhbI8TSfmo0RnELqfLEf1s//HjxyMrK6vRdVFRUQ2WVVZWwsnpRs9BJpPB\nYDBALm/+rbi5OULexIlFXl68xMjSsE0sU0vapaS8Bht3x+O3xFw42Mnw5AODcM/tgTy+30G68t/K\n+GFB+M8PF3Dhajnuvi1I7HJaxdLbRfTwby0nJydoNDdu+2gymW4Z/ABQUlLV6HIvLzUKCirMVh+1\nH9vEMrWkXX5PysP2AxegqTGgd6ArFkzsA09XBxQWVnZSlbalq/+tDApzx04A+4+mY2CIu9jltJil\ntEtzH0BEH/ZvrcGDB+Pnn38GAJw5cwbh4eEiV0RE5VU6vLc3AZv2JUJvNGHuuHAsnx0JT1cHsUsj\nK+bp4oA+wW5IySpDbnHjHThqG4vr+cfGxqKqqgozZ85sdP24ceNw9OhRzJo1C4Ig4LXXXuvkCono\nZqcuFGCjMsqIAAAgAElEQVT7gWSUV+nRI8AFC+/rAx83R7HLoi5i5AA/nM8owZGz2Xjgzh5il9Nl\nSAQbuYiyqSEYSxmeoRvYJpbpj+2iqdFjx38v4nhiHuQyKaaNDsU9Q7tDKuWx/c5iC38rOr0Rz7x7\nFAq5FOuWDIdMavkD1pbSLs0N+1tcz5+ILN/ZtEJs+y4ZZZU6hPipsfC+vvD3VIldFnVBSoUMw/r5\n4PDpqzh3qZhzRJgJw5+IWqyqxoDPDqfgl7M5kEklmDY6FPcOC7SK3hhZr1ED/HH49FX8cjaH4W8m\nDH8iapG4C/n412enUVyuRaCPEx65ry8CvDlhD3W8QB8ndPd2QnxqIco1OjirOF9EezH8iahZNToD\ndv2Yhp/irkImlWDKiGBMGh4MuYy9feocEokEIwf4YefBFBxLzMX42wLFLsnqMfyJqEkXLpfgo2+T\nUFhWg0BfNeZP6I0gX8uevIS6pjv6+eKLH1Nx5GwO7hnanZNGtRPDn4gayC2uwtdH03E8MQ8SCTBx\nWBAe+VN/lDYxWRZRR3NyUGBQTy+cTM5Hek4FQv2dxS7JqjH8iahOfkkVvj6agWOJuRAEINDbCTHj\neyGsmwsUTUyPTdRZRg3ww8nkfPxyNpvh304MfyJCSYUWe49cwtFzuTAJArp5qXD/yBBEhntByuFV\nshD9gt3hprbDb0l5mDm2J+wU/EDaVgx/IhtmMgn4Me4q9vychmqtEX4ejpg6MgRDensz9MniSKUS\njOjvi29+zcTpCwW4I8JX7JKsFsOfyEZl5Jbjk+8vIDO3Ao52cvx5fC+MHujPGfrIoo3o74dvfs3E\nkbPZDP92YPgT2ZhqrQFf/XwJh05nQRCAO/r54MHonnDhtdNkBXzcHNGruyuSL5civ7Qa3rx5VJsw\n/IlsRGFZNX48fRU/x2dDU2OAj7sj/nxPOPoEW8+tUomA2pv9XLhSiqNnc/Cn0aFil2OVGP5EXZgg\nCEjOLMHBU1k4k1oIQai9ZOpPo0Mx4bZAKOScqIesz5Be3tjx34s4mpCDqSNDeKiqDRj+RF2QVmfE\nr4m5OHwqC1cLNQCAIF817o4KwG19vHnZHlk1O6UMt/Xxwc/x2TifWYyIEA+xS7I6DH+iLiS/pAqH\nT1/FkbM5qNYaIJNKcHtfH4yNCkCYvzNnRaMuY9QAP/wcn41fzuYw/NuA4U9k5QRBQGJGMQ6dzMLZ\ntCIIAFxUSowbEow7I7vB1clO7BKJzC7U3xl+Ho44fbEAldV6ODkoxC7JqjD8iaxUtdaAXxNycfh0\nFnKKaqfdDfN3xtioAAzp7c0b71CXJpFIMGqAP3b9mIrfzudhbFSA2CVZFYY/kRUxmQQkXS7Br+dy\ncOpCAXQGE+QyCe7o54u7hwQgxI9TnpLtuCPCF7t/SsOR+GyGfyuZJfzfffddaDQaPPjggzh9+jSm\nT59ujs0S0TW5xVU4ei4HvybkoqRCCwDwdnXAiAF+GD3Qn9fok01yUSnRP9Qd8WlFyC2ugq+7o9gl\nWQ2zhL9CocCKFSuwefNmZGZmMvyJzEBTo8fvSfn49VwO0rLLAQAOdjKMHuiHEf390KObC0/gI5s3\nuJcX4tOKEJdSgHtvDxK7HKthlvDv168fAGDRokX4+OOPzbFJIpsjCAJyiqqQmFGMxPRinM8ogcFo\nggRAvxB3jIjwRWS4F29mQnSTgT08IZEAcRcLGf6tYJbwHzlyZN33Dz/8sDk2SWQTKqp0OJ9RUhf4\n14f0AcDfU4U7+vngjn6+cHe2F7FKIsvl7KhEz24uSMkqQ5lGx0NgLdSu8C8uLkZ6enrdv8zMTLz7\n7rvmqo2oy9EbTEi9WobE9GIkZhTjcm4FhGvrVPZy3NbHG32D3dEv2B0eLgx8opaIDPfCxawyxKcW\nYvRAf7HLsQptCv8HHngAzs7O8Pb2RnBwMA4cOIA1a9bA15d3WCK6mVZnxKXsMqRcLUNqVhkuZpVC\npzcBAGRSCXoFuqJfiDv6BrsjyEfNaUqJ2iCypyc+P5yKuIsFDP8WalP4z549G0eOHMGoUaNw7733\n4uTJk+jbt6+5ayOyOiUVWqRklSI1qzbwr+RVwiQIdev9PBzRL8QdESHuCO/uCnslr7Ylai9vN0d0\n81IhMaMENToD/65aoE0/oWnTpmHatGk4dOgQXnzxReTn58NoNEIm44lIZDtMJgFZBZVIvdarT8kq\nQ1F5Td16mVSCEH81enZzRY8AF/To5gJnHo8k6hCRPb3wza8ZSLhUjCG9vcUux+K1OvzT09MREhIC\nABg7dizGjh2LkydPYtWqVVCr1Vi1apXZiyQSmyAIKNfokFWoQdq1Xv2l7DJUa411z3FyUGBQD8+6\noA/xU/MGOkSdJLKnJ775NQNxKYUM/xZodfhPnDgRv/zyC0pLSxEcHAyZTIYhQ4ZgyJAhuHjxYkfU\nSNRp9AYT8kuqkFtchZyi+l+rtYZ6z/Vxd0RUuAt6BLigZ4ALfN0ded09kUiCfdVwU9vhbFohDEYT\np7e+hVaHvyAIGDt2LNzd3VFTU4MFCxbgkUceAQCEh4ebvUAicxMEAeVVeuQWaZBTXIXca+GeW1SF\ngrJq3HSIHkDt8L23mwP6BLnBz8MRoX7OCAtwgbMjh/CJLIVEIkFkT08cPn0VKVdK0SfYXeySLFqr\nw1+lUuHAgQPw9PREQUEBXn31VWzbtg3z58/viPqI2kxvMCG/tBq5RZp6PfjcoipU/aEXDwDOjgr0\n7OYCXw9H+Lqr4OvhCD8PR3i62EMmZS+CyNJF9vTC4dNXEZdSyPC/hVaHf/fu3eHp6QkA8PLywrp1\n6zB79ux2hX98fDzWrVuH7du311t++PBhbNy4EXK5HNOnT8eDDz4IvV6PlStX4urVq5BKpXj11VcR\nFhbW5n2T9aus1iOnSIPswvohX1DadC++V6Ar/DxU8HWvDXhfD0eo7HlLUCJr1ivQFQ52csSlFGD2\n3T15GK4ZrQ7/gIAAfPHFF3jggQcA1A61VFZWtrmALVu24Ouvv4aDg0O95Xq9HmvXrsXu3bvh4OCA\n2bNnIzo6GmfOnIHBYMBnn32Go0eP4l//+hc2bNjQ5v2TdRAEAWUaHbILNcgpqrr2tTbwy6v0DZ6v\ndlSgRzeX2mC/3ot3d4SnK3vxRF2VXCbFgDAP/HY+D5fzKhHkqxa7JIvV6vB/6aWXsHTpUuzYsQN9\n+vRBcnIyBg8e3OYCAgMDsWHDBjz//PP1lqelpSEwMBAuLi4AgKioKJw4cQLh4eEwGo0wmUyorKyE\nXM7rObsSkyAgt0iDxLRCZBdWIbvoesg3POFOAsDT1R4D/Jzh76mCn4djXW/eyYG9eCJbFNnTE7+d\nz0NcSgHDvxmtTk4fHx/s2rULcXFxSE5OxpgxYzB27Ng2FzB+/HhkZWU1WF5ZWQm1+kbDqVQqVFZW\nwtHREVevXsW9996LkpISbNq0qUX7cXNzhLyJy668vPgL0lKCIKBaa0BltR6aaj0qq/WorNJDU62r\n+76i6tr31XpU1xhgr5RB5aCAykEBp2tfJRIJqmpq11dpDaiq0SO/pBpXCyqh1Rnr7VMmlcDfS4UA\nbzUCfdQI8Kn96u+l4mQenYx/K5aHbVLfXWp7fPhNEs6lF+PRaQNFq8PS2+WW/3NqNBo4ODhAKpUi\nJSUFdnZ2CAwMRGRkJCIjIzusMCcnJ2g0mnp1qNVqfPzxxxg5ciSeffZZ5OTkYN68eYiNjYWdnV2z\n2yspqWp0uZeXGgUFFWat3ZoIggCdwYSKKh3KNXqUabQo1+hQptE1+Kqp1qNKa2hwHL05Uomk3gx3\nzVHIpfBzd0RwNxd4OCnh56GCv6cK3m4OjV62U1FWDdttuc5n638rloht0rjeQa5IuFSM8yn58HJ1\nuPULzMxS2qW5DyDNhv/69esRFxcHHx8fqNVq5OfnQ6VSwdfXF0899ZTZC71ZWFgYMjMzUVpaCkdH\nR5w8eRILFy5EWloaFIraIV0XFxcYDAYYjcZbbK1rEQQBBqMJ1TojtDojanRG1OgMN31f+7haZ0SN\ntvZrtdaAaq2h3uOaa1+NpubDWSIB1I5KuDjZwc9DBUd7OVT2cjjaKeBoL6/752Rf+1jloIDKXgGV\nvRxKhQx6gxFVWiOqavSoqjFAU2OAIAhwsJPX/lPKYG8nh6OdHFKpxGL+cIjIOg3u6YWES8U4k1KI\ncUO7i12ORWo2/I8fP46dO3fCaDRi4sSJOHDgAAAgJiamwwqKjY1FVVUVZs6ciZUrV2LhwoUQBAHT\np0+Hj48PHn74Ybz44ouYM2cO9Ho9nn76aTg6OnZYPe1lNJmg1Zmg1Rtr/+mMDb6v0Rmhu/a1bp3e\n2CDMa24K+5b2pv9IAsDeTgZ7pRwuKiV83B1gr5RD7aiAi0oJZ5Xypq92cFYpoXZQtOuGMwq5DC5y\nGW+1SUSdYlBPT/zfgQuISylg+DehRQdMZTIZli5d2mFFBAQEYNeuXQCAyZMn1y2Pjo5GdHR0veeq\nVCqsX7++w2ppiZwiDX6Ky0a1zgDdHwK9fpCbYDCa2r0/hVwKO4UM9koZ3J3tYa+s/d7u2ld7hfxa\noMuuPU8Oe6UMDna1yx2UtT3s66+R8vIXIurCXJ3sEOrvjAtXSlFZrecJwI1oNvxHjhyJs2fPon//\n/pg0aRIAQKfTYe7cuZ1SnKU6npiH/5680mD59ZC2U8jg6mQHpeJGINspZbBTSGGnkMNOeeN5dnWB\nXf/xzd9zmkoiotaJ7OmJS9nliE8txIj+fmKXY3GaDf8JEybg2LFj2LlzJwAgIiICw4YNw4QJEzql\nOEs1ZWQwhvbxhkImrRfWvBc7EZFlGBzuhS//dwlxKQz/xjQb/mFhYXWz55lMJiQkJODgwYPYvHkz\nxo0bh7vvvrtTirQ0MqkUAV5OYpdBRERN8PNQwcfdEQnpRdDpjVAqeIfNm7V4PFkqlWLAgAHw8PDA\nG2+8YbPBT0RE1mFwT0/o9CaczygRuxSL02zPf/78+XVzIwvXzi5PTk7G/v37sXXr1o6vjoiIqI0i\nw73w3W+XcTqlAIN6eopdjkVpNvzHjBmDCxcuYNKkSRgxYgQA4JFHHsGHH37YKcURERG1Vai/M5xV\nSsSnFsJkEnhe1k2aDf+HH34YBoMB+/btwzPPPIPo6Oi6EQAiIiJLJpVIMKiHJ36Oz0bq1TKEd3cV\nuySLcctj/tdvp/vWW29BoVBg6NChnVEXERFRuw0Orx3uj0spELkSy9Liu6JIJBKMHz++I2shIiIy\nqz5BbrBTyBCXUoiZ0T3FLsditGr2mC1btnRUHURERGankMvQN9gN+SXVyG/iBm+2qFXhf+TIkXqP\ntVotdDpdvWX5+fntr4qIiMhM+oW4AwAv+btJq8L/5pP9Pv74Yzz66KNYvHgxXn75ZdTU1AAAli9f\nbt4KiYiI2qFfcG34J2YUi1yJ5WjxMX8Addf8A8CBAwfqpv09evQoFi9ejFdeecWsxREREbWXt5sD\nPJztkJxZwkv+rmlV+N/c8zeZTDAYDJDL5RgxYgTCwsLwwgsvICMjw9w1EhERtZlEIkHfYHccOZuD\nzLwKhPg5i12S6Fo17H9zz/75559HScmN4ye+vr744IMP8Je//MVsxREREZlD3+Drx/059A+0Mvyv\n3+QHAKKiouDl5VVvvVKpxLRp08xTGRERkZn0CXYDACSmM/yBVg773ywvLw9JSUlITk5GUlIS1q9f\nb866iIiIzMbZUYlAbyekXi2DVm+EnY3f5a/F4b937966oL9w4QK0Wi2MRiOmTJmCyMjIjqyRiIio\n3fqGuONyfiVSrpQiItRD7HJE1eLwX7lyJe644w489NBDCA8PR/fu3XH33XdjzZo1HVkfERGRWfQN\ndsP3v13G+YwSmw//Fh/z3759O0pLS/Htt9/Czs4OUqm03qV/REREliw8wBVymZTX+6MV4T906FDs\n2bMHw4YNw5w5c/Dvf/8bJpOpI2sjIiIyG6VChp4BLriSX4lyje7WL+jCWnW2v0QiwaxZs7Bv3z5U\nVFSgpKQEe/fu7ajaiIiIzKrvtbP+z2fadu+/VeF/nVqtxqpVq/Dll18iNjYWU6dONXddREREZnfj\nen/bnue/zZf6AbXX/X/00Uc4fPiwueohIiLqMEE+aqjs5TifUQxBEGz23LU29fz/KDo62hybISIi\n6lBSqQR9gtxQXK5FXkm12OWIxizhT0REZC36XrvFry3P9sfwJyIim8J5/hn+RERkY7xdHeDlao/k\nyyUw2ugl6wx/IiKyOf2C3VGtNSIjp0LsUkTR4vCfNWsWVq5cif3796OsrKwjayIiIupQ14f+bXW2\nvxZf6vfZZ58hMzMTP//8M1asWIGamhrcdtttGD16NCIiIjqyRiIiIrPqHeQGCWqv958yIkTscjpd\nq4b9g4KCEBMTg02bNuGDDz5AREQE9u7di9mzZ7eriPj4eMTExDRYfvjwYUyfPh0zZ87Erl276pZ/\n8MEHmDlzJqZNm4YvvviiXfsmIiLb4+SgQJCvGmlXy1CjM4hdTqdr8yQ/dnZ2GD16NEaPHt2uArZs\n2YKvv/4aDg4O9Zbr9XqsXbsWu3fvhoODA2bPno3o6GikpaUhLi4OO3fuRHV1NbZu3dqu/RMRkW3q\nF+KOjNwKXLxSigFhnmKX06lEP+EvMDAQGzZsaLA8LS0NgYGBcHFxgVKpRFRUFE6cOIFffvkF4eHh\nWLJkCRYvXow777yz84smIiKr1zeodp7/xHTbm+q3XdP7msP48eORlZXVYHllZSXUanXdY5VKhcrK\nSpSUlCA7OxubNm1CVlYWHn/8cXz//fe3nKLRzc0Rcrms0XVeXupGl5N42CaWie1iedgmbefi6gjl\nl+dwMavU7D9HS28X0cO/KU5OTtBoNHWPNRoN1Go1XF1dERoaCqVSidDQUNjZ2aG4uBgeHh7Nbq+k\npKrR5V5eahQU2OalHpaKbWKZ2C6Wh23Sfj0DXJCYXoyU9EK4OtmZZZuW0i7NfQARfdi/KWFhYcjM\nzERpaSl0Oh1OnjyJyMhIREVF4ciRIxAEAXl5eaiuroarq6vY5RIRkRXqd+2SvyQbu8ufxfX8Y2Nj\nUVVVhZkzZ2LlypVYuHAhBEHA9OnT4ePjAx8fH5w4cQIzZsyAIAhYvXo1ZLLGh/OJiIia0yuwtvN4\nMasUd0T4ilxN55EIgiCIXURnaGoIxlKGZ+gGtollYrtYHrZJ+xlNJjz5zhF4uNhjzSO3m2WbltIu\nHT7sbzKZ6v7t3r3bHJskIiLqcDKpFKH+zsgu1KCyWi92OZ3GLMP+8+bNw+23135iio+Px4wZM8yx\nWSIiog4X3t0VSZklSMkqRWRPL7HL6RRmCf/ly5dj4MCBAIBLly6ZY5NERESdomeACwAgJauM4d8S\nxcXFSE9PR3p6On744QdkZmbi3XffNVdtREREHS7U3xlSiQQpWaVil9Jp2hT+DzzwAJydneHt7Y3g\n4GAcOHAAa9asga+v7ZwpSUREXYO9Uo4gXydk5FRApzdCqej6V5C16YS/2bNnw9nZGaNGjcKiRYvg\n7u6Ovn37wt3d3dz1ERERdbieAa4wmgSk55SLXUqnaFP4T5s2De+88w7s7Ozw4osvoqCgAEaj0dy1\nERERdYrrx/0vZpWJXEnnaPWwf3p6OkJCau99PHbsWIwdOxYnT57EqlWroFarsWrVKrMXSURE1JF6\nBNRO9pNyxTaO+7c6/CdOnIhffvkFpaWlCA4Ohkwmw5AhQzBkyBBcvHixI2okIiLqUC4qJXzcHZF6\ntQwmkwCptPmbxVm7Vg/7C4KAsWPH4tFHH8WoUaPw4Ycf1q0LDw83a3FERESdpWeAC2p0RmQVVIpd\nSodrdfirVCocPHgQhw8fxr59+3D27Fls27atI2ojIiLqNDdf79/VtTr8u3fvDk9PTwCAl5cX1q1b\nh2+++cbshREREXWm8O7XbvJjA8f9Wx3+AQEB+OKLL+oeSyQSVFZ2/SESIiLq2rxdHeCsUiIlqxRd\n/Z53rT7h76WXXsLSpUuxY8cO9OnTB8nJyRg8eHBH1EZERNRpJBIJega44NSFAhSW1cDL1UHskjpM\nq8Pfx8cHu3btQlxcHJKTkzFmzBiMHTu2I2ojIiLqVD0DXHHqQgEuXill+DcmMjISkZGR5qyFiIhI\nVOHdb5z0N6K/n8jVdJxbHvPXaDQwmUwAgJSUFFy+fLnDiyIiIhJDd28n2ClkXf4mP832/NevX4+4\nuDj4+PhArVYjPz8fKpUKvr6+eOqppzqrRiIiok4hk0oR1s0Z5zNKUFGlg9pRKXZJHaLZ8D9+/Dh2\n7twJo9GIiRMn4sCBAwCAmJiYTimOiIios4UHuOJ8RglSs8oQGe4ldjkdokWX+slkMixdurSjayEi\nIhKdLUz202z4jxw5EmfPnoUgCJg0aRIAQKfTYe7cuZ1SHBERUWcL9XeBVCLp0sf9mx32nzBhAo4d\nO4adO3cCACIiIjBs2DBMmDChU4ojIiLqbHZKGYJ8nZCRWwGt3gg7hUzsksyu2fAPCwtDWFgYAMBk\nMiEhIQEHDx7E5s2bMW7cONx9992dUiQREVFn6hngivScCqRnl6N3kJvY5Zhdi6f3lUqlGDBgADw8\nPPDGG28w+ImIqMvqGVA7z39XHfpvtuc/f/58SCS19zS+Ps9xcnIy9u/fj61bt3Z8dURERCK4ftLf\nxS560l+z4T9mzBhcuHABkyZNwogRIwAAjzzyCD788MNOKY6IiEgMziolvF0dkJ5dDkEQ6jrCXUWz\n4f/www/DYDBg3759eOaZZxAdHd3l73REREQEACH+zvjtfB7yS6rh4+4odjlmdctj/nK5HNOnT8db\nb70FhUKBoUOHdkZdREREogr1cwYAXMouF7kS82vxjX0kEgnGjx/fkbUQERFZjBD/a+GfU447InxF\nrsa8Wny2PwBs2bKlo+ogIiKyKIHeTpBJJUjP6Xo9/1aF/5EjR275nNTU1FYXER8f3+j9Ag4fPozp\n06dj5syZ2LVrV711RUVFGDNmDNLS0lq9PyIioltRKmQI8HbC5bwK6A0mscsxq1aFf0tO9nvttdda\nVcCWLVvw17/+FVqttt5yvV6PtWvXYuvWrdi+fTs+//xzFBYW1q1bvXo17O3tW7UvIiKi1gj1c4bB\nKCCroFLsUsyqxcf8AdS71GHw4MEYNGhQvUsgBEFAUlJSqwoIDAzEhg0b8Pzzz9dbnpaWhsDAQLi4\n1F5rGRUVhRMnTuDee+/FG2+8gVmzZmHz5s2t2hcREVFrhPo748e4q7iUXY6QaycAdgWtCv+be/7B\nwcH417/+BWfn+j+M+fPnt6qA8ePHIysrq8HyyspKqNXquscqlQqVlZXYs2cP3N3dMWrUqFaFv5ub\nI+Tyxudn9vJSN7qcxMM2sUxsF8vDNulYUf2Aj75NQnZxVat+1pbeLq0K/1deeaXu+23btkGlUjV4\nzrZt29pdFAA4OTlBo9HUPdZoNFCr1di+fTskEgmOHTuGpKQkrFixAu+//z68vJq/53JJSVWjy728\n1CgoqDBLzWQebBPLxHaxPGyTjqeUCHCwkyEpvbjFP2tLaZfmPoC0Kvyv3+QHQN1wfEcJCwtDZmYm\nSktL4ejoiJMnT2LhwoX17igYExODV1555ZbBT0RE1BZSiQTBvs5IyixBVY0ejvYKsUsyi1aF/83y\n8vKQlJSE5ORkJCUlYf369WYpKDY2FlVVVZg5cyZWrlyJhQsXQhAETJ8+HT4+PmbZBxERUUuF+teG\nf3pOBfqFuItdjlm0OPz37t1bF/QXLlyAVquF0WjElClTEBkZ2a4iAgIC6i7lmzx5ct3y6OhoREdH\nN/m67du3t2u/REREt3L9RL9LOeW2F/4rV67EHXfcgYceegjh4eHo3r077r77bqxZs6Yj6yMiIhJV\n6LWZ/tK70DS/Lb7Of/v27SgtLcW3334LOzs7SKXSLneXIyIioj9ydbKDm9oOl3LKu8zN7Voc/kOH\nDsWePXswbNgwzJkzB//+979hMnWtGY+IiIgaE+rnjHKNDsXl2ls/2Qq0aoY/iUSCWbNmYd++faio\nqEBJSQn27t3bUbURERFZhLqh/y4yz3+rwv86tVqNVatW4csvv0RsbCymTp1q7rqIiIgsRkgXu71v\nmy/1A2qvxf/oo49w+PBhc9VDRERkcYL91JBIas/47wra1PP/o+YuxyMiIrJ29ko5/D1VyMgth7EL\nnO9mlvAnIiLq6kL8nKHTm5Bd2Ph08daE4U9ERNQCXemkP4Y/ERFRC4R2oZP+GP5EREQt0M1LBaVc\nyvAnIiKyFTKpFEG+alwtrIRWZxS7nHZh+BMREbVQiJ8zBAHIyLXu3j/Dn4iIqIWC/dQAgMt5lSJX\n0j4MfyIiohYK9L4e/hUiV9I+DH8iIqIW8nV3hFIhRSZ7/kRERLZBKpWgu5cTcoo00Busd6Y/hj8R\nEVErBPqoYTQJuFpovb1/hj8REVErBPo4AbDuk/4Y/kRERK0Q6FN70l+mFZ/0x/AnIiJqhQAvFaQS\niVWf8c/wJyIiagWFXAZ/T0dcya+EySSIXU6bMPyJiIhaKdBHDZ3ehLwS67y9L8OfiIiolaz9uD/D\nn4iIqJWCrPyMf4Y/ERFRK3W38ml+Gf5ERESt5Ggvh5erPS7nVUIQrO+kP4Y/ERFRGwT6qFFZrUdJ\nhVbsUlqN4U9ERNQG1nzSH8OfiIioDaz5pD+GPxERURtc7/lb40l/FhH+8fHxiImJabD88OHDmD59\nOmbOnIldu3YBAPR6PZ577jnMmTMHM2bMwKFDhzq7XCIiIrg62cFZpbTK8JeLXcCWLVvw9ddfw8HB\nod5yvV6PtWvXYvfu3XBwcMDs2bMRHR2N//3vf3B1dcWbb76J0tJS3H///Rg7dqxI1RMRkS0L9HFC\nwqViVFbr4eSgELucFhO95x8YGIgNGzY0WJ6WlobAwEC4uLhAqVQiKioKJ06cwIQJE/DUU08BAARB\ngHkp5BEAAA0vSURBVEwm6+ySiYiIAABBVjr0L3rPf/z48cjKymqwvLKyEmq1uu6xSqVCZWUlVCpV\n3fply5bhL3/5S4v24+bmCLm88Q8KXl7qRpeTeNgmlontYnnYJuKK6OmFb49loqhSX68tLL1dRA//\npjg5OUGj0dQ91mg0dR8GcnJysGTJEsyZMweTJ09u0fZKmrj5gpeXGgUF1vWJratjm1gmtovlYZuI\nz9WhNkaTLhWiIMIHgOW0S3MfQEQf9m9KWFgYMjMzUVpaCp1Oh5MnTyIyMhKFhYVYsGABnnvuOcyY\nMUPsMomIyIZ5uTrAXimzumv9La7nHxsbi6qqKsycORMrV67EwoULIQgCpk+fDh8fH6xZswbl5eV4\n77338N577wGoPWnQ3t5e5MqJiMjWSCUSBHo7IeVqGbR6I+wU1nEemkSwxkmJ26CpIRhLGZ6hG9gm\nlontYnnYJpbhPwcv4uDJLDw3OxJ9gtwspl2sctifiIjIGgwM8wQAnLqQL3IlLcfwJyIiaofeQa5w\nclDg1IUCmEzWMZjO8CciImoHmVSKweGeKNPokJJVKnY5LcLwJyIiaqehvWsv8zuRbB1D/wx/IiKi\ndrp56N9oBUP/DH8iIqJ2unnoPym9SOxybonhT0REZAbXh/5/ic8WuZJbY/gTERGZwfWh/1/PZlv8\nWf8MfyIiIjO4PvRfUqG1+LP+Gf5ERERmYi1n/TP8iYiIzKR3kCvUjkqLn/CH4U9ERGQmMqkUd/T3\ns/gJfxj+REREZjRyoD8Ayx76Z/gTERGZ0YAenhY/1z/Dn4iIyIxkMikGh3tZ9NA/w5+IiMjMhvb2\nBmC5Q/8MfyIiIjOz9Nv8MvyJiIjMrHbCH8sd+mf4ExERdQBLHvpn+BMREXUASx76Z/gTERF1AEse\n+mf4ExERdRBLHfpn+BMREXUQSx36Z/gTERF1EEsd+mf4ExERdSBLHPpn+BMREXUgSxz6Z/gTERF1\nIEsc+mf4ExERdTBLG/pn+P9/e/caE9WZx3H8xz2sXNYo2H3hYNcWXKOJYnSTdkksLOFFxUQGGdx0\nUlJUmpiQbEgKTRGbXtDGsm9I1drEkNSm1lYbS3pLTNOypVpFBEOEJhid1K11SdTggHZQzr5oO5Vy\nb9d5DjzfzyvOwzDP/+SfyY/z5+QAAMB95rbRP+EPAMB95rbRvyvCv7OzU36/f9T6Z599Jq/XK5/P\np8OHD0uShoeHVVdXJ5/PJ7/fr0AgEOlyAQCYNjeN/mNNF/DGG2/ogw8+UGJi4oj1oaEh7dy5U++9\n954SExO1adMm5ebmqr29XaFQSO+88446Ojq0a9cu7d2711D1AABMzb2j/3/8PVPR0VHGajF+5e/x\neNTY2Dhq/cKFC/J4PEpNTVV8fLxWrVql06dP68yZM8rJyZEkrVixQl1dXZEuGQCAaXPT6N/4lX9B\nQYEuX748aj0YDCo5OTl8PGfOHAWDQQWDQSUlJYXXY2JidOfOHcXGTnwqc+f+QbGxMWN+Ly0tecx1\nmENP3Im+uA89cafx+vL3v2aopfM7dQVu6G+rPBGu6hfGw388SUlJGhgYCB8PDAwoOTl51Prw8PCk\nwS9J168Pjrmelpasvr6bv79g/N/QE3eiL+5DT9xpor786Y8JSkqM0787/qMNjy66r6P/iX4xND72\nH8/ixYsVCAR048YNhUIhtbW1aeXKlcrOzlZLS4skqaOjQ5mZmYYrBQBgan4e/fcbHv277sq/ublZ\ng4OD8vl8qqmpUXl5uRzHkdfr1YIFC5Sfn6/W1laVlpbKcRzV19ebLhkAgClbvSRdLZ3f6XTPf5Xl\nmWukhijHccw/bSACxhvBMDZzH3riTvTFfeiJO03Wl7vDw/pnY6uio6P0r22P3rfR/4wc+wMAMBu5\nYfRP+AMAEGGr/2L2gT+EPwAAEbbE8+MDf9oMPeuf8AcAIMJMj/4JfwAADDA5+if8AQAwwOTon/AH\nAMAAk6N/wh8AAENMjf4JfwAADDE1+if8AQAwxNTon/AHAMAgE6N/wh8AAIN+Hv13XbwWsT1d91/9\nAACwSUx0tLYULtXNwVDE9iT8AQAwbPmf50V0P8b+AABYhvAHAMAyhD8AAJYh/AEAsAzhDwCAZQh/\nAAAsQ/gDAGAZwh8AAMsQ/gAAWIbwBwDAMoQ/AACWiXIcxzFdBAAAiByu/AEAsAzhDwCAZQh/AAAs\nQ/gDAGAZwh8AAMsQ/gAAWIbwBwDAMoQ/AACWIfzHcOLECdXW1qqqqko9PT2my8E9Tpw4oeeee850\nGdZrb29XdXW1qqur1d/fb7oc3IPPiLu4NU8I/zHcunVLL774osrLy/Xll1+aLgc/CQQC6u7u1g8/\n/GC6FOsdPnxYL7zwgoqLi/XRRx+ZLgc/4TPiPm7NE8JfUlNTkyoqKlRRUaG9e/cqNzdXt27d0ptv\nvqkNGzaYLs9av+5LRkaGnnrqKdNlQdLdu3eVkJCgtLQ09fX1mS4HP+Ez4j5uzZNY0wW4QVlZmcrK\nysLH165d0+7du1VZWal58+aZK8xyv+4L3CMxMVGhUEh9fX2aP3++6XIA13Jrnsz6K//Ozk75/X5J\n0vDwsOrq6uTz+eT3+xUIBMb8mV27dqmvr08NDQ365JNPIlmuNX5LXxAZU+lNSUmJ6urqdOjQIa1f\nv95kudbgM+M+U+mJa/PEmcX279/vrFu3ztm4caPjOI7z6aefOtXV1Y7jOM7Zs2edp59+2mR51qIv\n7kVv3Im+uM9M78msvvL3eDxqbGwMH585c0Y5OTmSpBUrVqirq8tUaVajL+5Fb9yJvrjPTO/JrA7/\ngoICxcb+cltDMBhUUlJS+DgmJkZ37twxUZrV6It70Rt3oi/uM9N7MqvD/9eSkpI0MDAQPh4eHh7R\nPJhBX9yL3rgTfXGfmdYTq8I/OztbLS0tkqSOjg5lZmYarggSfXEzeuNO9MV9ZlpP3PtryX2Qn5+v\n1tZWlZaWynEc1dfXmy4Joi9uRm/cib64z0zrSZTjOI7pIgAAQORYNfYHAACEPwAA1iH8AQCwDOEP\nAIBlCH8AACxD+AMAYBnCHwAAyxD+AABYhvAHMKG2tjaVlJSotLRUBw4csG5/YDbiCX8AJnT16lXN\nnTtX8fHx8vv92r9/vxITE63ZH5iNrHq2P4DpW7BgQfjrmJgYRUdHdmBoen9gNuJTBEDd3d2qqKhQ\nYWGhCgsLVVlZqatXr454TWtrqzwejxISEiRJWVlZKiws1FdffSVJGhoaUmNjowoKCvT4449r/fr1\nqqysVG9v76T7b968WW+//faINcdxlJeXp1OnTo25vySVlZVpzZo1Onjw4O86f8A2hD9guba2Nm3b\ntk3l5eVqbm5Wc3OzHn74YVVWVoZf8/333+v1119XdXX1iJ89dOiQHnnkEUnSs88+q2+++Ubvvvuu\nPvzwQx07dkxFRUW6ePHipDV4vV69//77I9a+/vprRUdHa/Xq1ePu39TUpNzc3N966oC1CH/AYqFQ\nSM8884y2b9+uNWvWhNe3bt2qc+fO6cqVKwqFQqqpqdHzzz+vOXPmjPk+ly5d0vHjx/Xyyy8rJSVF\nkhQVFaW1a9cqPz8//LrOzk75/X4VFRWpqKhIn3/+uSQpLy9PgUBAFy5cCL/26NGjKioq0tDQ0KT7\nA5gewh+w2BdffKG4uDg99thjI9Z//tv64OCgmpub1dvbqx07dsjv94/6c4AknT9/XhkZGUpNTR13\nr/7+fu3YsUMNDQ06evSo9u3bp7q6OvX39ys+Pl6FhYU6cuSIJCkYDOr48ePasGHDlPYHMD3c8AdY\n7Pz581q2bNmo9e7ubsXFxWnhwoVavHixvF7vtN63t7dXVVVVun37tnJyclRbW6uzZ8/q8uXL2rJl\nS/h1UVFRCgQCWr58uYqLi7V582ZVVVXp448/VnZ2th544AF5vd5p7w9gYoQ/YLGUlBSFQqFR6wcO\nHJDP51N8fPyU3mfp0qUKBALq7+9XSkqKHnroIR07dkwHDx5UV1eXpB9v4MvKytJbb7015nssWbJE\n6enpamlp0ZEjR/Tkk0/+9hMDMCHG/oDF1q5dq5MnT6qnp0fSj3fsv/baa/r2229H3PA3mUWLFikv\nL0+1tbW6efNmeH1wcDD89cqVKxUIBHTy5Mnw2rlz53Tvo0a8Xq8aGxt16dIl5eXl/Z5TAzABrvwB\niz344INqaGjQ9u3bdfv2bV25ckUbN25UU1PTtG+u27lzp/bs2aPi4mLFxsYqJSVF6enp2rp1qyQp\nNTVVe/bs0e7du1VfX6+hoSEtXLhQ+/btU1RUlCRp3bp1euWVV1RSUjLlqQOA6eMJfwDCXnrpJV2/\nfl2vvvpqOJDHk5WVpfb2duN34NfU1GjZsmV64oknjNYBzCSM/QGE1dbWqqGhYdLgl6T58+dr06ZN\n4Yf8mFBWVqZTp07xuF9gmrjyBwDAMlz5AwBgGcIfAADLEP4AAFiG8AcAwDKEPwAAliH8AQCwDOEP\nAIBlCH8AACxD+AMAYBnCHwAAy/wPtopdQoilc9wAAAAASUVORK5CYII=\n",
0316 "text/plain": [
0317 "<matplotlib.figure.Figure at 0x119baa9b0>"
0318 ]
0319 },
0320 "metadata": {},
0321 "output_type": "display_data"
0322 }
0323 ],
0324 "source": [
0325 "fig, ax = plt.subplots()\n",
0326 "fig.suptitle(r'$ \\gamma+A \\rightarrow \\rho +A$')\n",
0327 "plt.xlabel('$Q^{2}[\\\\rm{GeV}^{2}]$')\n",
0328 "plt.ylabel('$A_1^{-4/3}\\sigma_{A_1}$ / $A_2^{-4/3}\\sigma_{A_2}$')\n",
0329 "ax.semilogx(x_bins, y_val)\n",
0330 "ax.legend([r'A=197(Au)/A=56(Fe - density $\\rho_0$ = 0.138'])\n",
0331 "plt.savefig(home_dir+'Au_Fe_scaling_rho.eps')\n",
0332 "plt.show()"
0333 ]
0334 },
0335 {
0336 "cell_type": "markdown",
0337 "metadata": {},
0338 "source": [
0339 "## Save data\n",
0340 "We now save the data points to a csv file for plotting in the future. More specifically, paper plots are mostly done in ROOT so we will plot this with similar style"
0341 ]
0342 },
0343 {
0344 "cell_type": "code",
0345 "execution_count": null,
0346 "metadata": {
0347 "collapsed": false
0348 },
0349 "outputs": [],
0350 "source": [
0351 "import csv\n",
0352 "data = [x_bins, y_val]\n",
0353 "with open(home_dir+'Au_Fe_scaling_jpsi.csv','w') as f:\n",
0354 " out = csv.writer(f, delimiter = ',')\n",
0355 " out.writerows(zip(*data))\n"
0356 ]
0357 },
0358 {
0359 "cell_type": "markdown",
0360 "metadata": {},
0361 "source": [
0362 "## Prepare the file \n",
0363 "This step prepares the csv file to plot in ROOT using the new procedure (03-01-2018)."
0364 ]
0365 },
0366 {
0367 "cell_type": "code",
0368 "execution_count": null,
0369 "metadata": {
0370 "collapsed": false
0371 },
0372 "outputs": [],
0373 "source": [
0374 "import csv\n",
0375 "import math as m\n",
0376 "particles = ['rho','J_psi']\n",
0377 "temp_dir = home_dir+'Au_vs_Fe_Tests/'\n",
0378 "au_scale = 197**(-4/3)\n",
0379 "fe_scale = 56**(-4/3)\n",
0380 "for part in particles: \n",
0381 " x_bin = []\n",
0382 " iron = []\n",
0383 " gold = []\n",
0384 " f_iron = csv.reader(open(temp_dir+part+'_Iron/estarlight_total_xsec.csv','r'))\n",
0385 " for row in f_iron:\n",
0386 " q0 = float(row[0])\n",
0387 " q1 = float(row[1])\n",
0388 " num = m.log(q1) - m.log(q0)\n",
0389 " denom = (q0**(-1)-q1**(-1))\n",
0390 " if q0 < 90:\n",
0391 " x_bin.append(num/denom)\n",
0392 " iron.append(float(row[2]))\n",
0393 " f_gold = csv.reader(open(temp_dir+part+'_Gold/estarlight_total_xsec.csv','r'))\n",
0394 " for row in f_gold:\n",
0395 " q0 = float(row[0])\n",
0396 " if q0 < 90: \n",
0397 " gold.append(float(row[2]))\n",
0398 " y_bin = [ (gold[i]*au_scale)/(iron[i]*fe_scale) for i in range(len(iron))]\n",
0399 " with open(temp_dir+'Au_Fe_scaling_'+part+'.csv','w') as f:\n",
0400 " out = csv.writer(f, delimiter = ',')\n",
0401 " out.writerows(zip(*[x_bin,y_bin]))\n",
0402 " \n",
0403 " \n",
0404 "#for row in f:\n",
0405 "\n",
0406 "\n",
0407 "\n",
0408 " "
0409 ]
0410 },
0411 {
0412 "cell_type": "code",
0413 "execution_count": null,
0414 "metadata": {
0415 "collapsed": false
0416 },
0417 "outputs": [],
0418 "source": []
0419 },
0420 {
0421 "cell_type": "code",
0422 "execution_count": null,
0423 "metadata": {
0424 "collapsed": false
0425 },
0426 "outputs": [],
0427 "source": []
0428 },
0429 {
0430 "cell_type": "code",
0431 "execution_count": null,
0432 "metadata": {
0433 "collapsed": false
0434 },
0435 "outputs": [],
0436 "source": []
0437 },
0438 {
0439 "cell_type": "code",
0440 "execution_count": null,
0441 "metadata": {
0442 "collapsed": false
0443 },
0444 "outputs": [],
0445 "source": []
0446 },
0447 {
0448 "cell_type": "code",
0449 "execution_count": null,
0450 "metadata": {
0451 "collapsed": false
0452 },
0453 "outputs": [],
0454 "source": []
0455 },
0456 {
0457 "cell_type": "code",
0458 "execution_count": null,
0459 "metadata": {
0460 "collapsed": true
0461 },
0462 "outputs": [],
0463 "source": []
0464 }
0465 ],
0466 "metadata": {
0467 "kernelspec": {
0468 "display_name": "Python 3",
0469 "language": "python",
0470 "name": "python3"
0471 },
0472 "language_info": {
0473 "codemirror_mode": {
0474 "name": "ipython",
0475 "version": 3
0476 },
0477 "file_extension": ".py",
0478 "mimetype": "text/x-python",
0479 "name": "python",
0480 "nbconvert_exporter": "python",
0481 "pygments_lexer": "ipython3",
0482 "version": "3.6.0"
0483 }
0484 },
0485 "nbformat": 4,
0486 "nbformat_minor": 2
0487 }
