{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quenching and drifting examples\n",
    "In this notebook we apply the quenching and drifting stages to a track dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This is need so you can import larndsim without doing python setup.py install\n",
    "import os,sys,inspect\n",
    "currentdir = os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))\n",
    "parentdir = os.path.dirname(currentdir)\n",
    "sys.path.insert(0,parentdir) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from larndsim import consts\n",
    "consts.load_detector_properties(\"../larndsim/detector_properties/module0.yaml\",\"../larndsim/pixel_layouts/multi_tile_layout-2.2.16.yaml\")\n",
    "\n",
    "from larndsim import quenching, drifting\n",
    "import importlib\n",
    "importlib.reload(drifting)\n",
    "importlib.reload(quenching)\n",
    "\n",
    "from math import ceil\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import numpy as np\n",
    "import h5py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "with h5py.File('module0_corsika.h5', 'r') as f:\n",
    "    tracks = np.array(f['segments'])\n",
    "\n",
    "x_start = np.copy(tracks['x_start'] )\n",
    "x_end = np.copy(tracks['x_end'])\n",
    "x = np.copy(tracks['x'])\n",
    "\n",
    "tracks['x_start'] = np.copy(tracks['z_start'])\n",
    "tracks['x_end'] = np.copy(tracks['z_end'])\n",
    "tracks['x'] = np.copy(tracks['z'])\n",
    "\n",
    "tracks['z_start'] = x_start\n",
    "tracks['z_end'] = x_end\n",
    "tracks['z'] = x\n",
    "\n",
    "tracks = tracks[(tracks['dx']>0)&(tracks['dE']>1e-3)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks_box = np.copy(tracks)\n",
    "tracks_birks = np.copy(tracks)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "threadsperblock = 256\n",
    "blockspergrid = ceil(tracks.shape[0] / threadsperblock)\n",
    "\n",
    "quenching.quench[blockspergrid,threadsperblock](tracks_box, consts.box)\n",
    "quenching.quench[blockspergrid,threadsperblock](tracks_birks, consts.birks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recombination\n",
    "\n",
    "The number of ionized electrons must be corrected by the recombination effect. Two models have been implemented:\n",
    "\n",
    "- Birks (Amoruso et al. NIM A 523 (2004) 275)\n",
    "\\begin{equation}\n",
    "R_{Birks} = \\frac{A_b}{1+k_b/\\epsilon\\cdot dE/dx }\n",
    "\\end{equation}\n",
    "\n",
    "- Modified box (Baller, 2013 JINST 8 P08005)\n",
    "\\begin{equation}\n",
    "R_{box} = \\frac{\\log(\\alpha + \\beta \\frac{dE/dx}{E\\rho})}{\\beta/\\epsilon\\cdot dE/dx}\n",
    "\\end{equation}\n",
    "\n",
    "where $\\epsilon=E\\rho_{\\mathrm{LAr}}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEKCAYAAAAIO8L1AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8nklEQVR4nO3deXxV9Z34/9c7N3tIIGRjCUuAkAXCDi4oIItCddS2ttXWUaut9Te1bp12dLpo/bbT1um003GoHaq01r1Vq1RRqLggosguAmEJBJKwhSxkIUCW9++Pc5JcYgg3yb25Wd7Px+M8cs+5Z3nnPiDv+9lFVTHGGGNaCgl2AMYYY7onSxDGGGNaZQnCGGNMqyxBGGOMaZUlCGOMMa2yBGGMMaZVocEOoL0SExN15MiRwQ7DGGN6lI0bNx5X1aT2XNPjEsTIkSPZsGFDsMMwxpgeRUQOtPcaq2IyxhjTKksQxhhjWmUJwhhjTKsC2gYhIguB3wIe4HFV/UWL94cDTwID3HPuV9XlgYzJGNM1amtrKSws5NSpU8EOpU+JjIwkNTWVsLCwTt8rYAlCRDzAYmABUAisF5FlqrrD67QfAn9R1cdEJBtYDowMVEzGmK5TWFhIbGwsI0eORESCHU6foKqUlJRQWFhIWlpap+8XyCqmGcBeVd2nqmeA54FrWpyjQJz7uj9wKIDxGGO60KlTp0hISLDk0IVEhISEBL+V2gKZIIYCBV77he4xbw8BN4pIIU7p4Tut3UhEbheRDSKyobT4MDx7PRRtDETMxhg/suTQ9fz5mQe7kfoG4E+qmgp8DnhKRD4Tk6ouUdVpqjptYFgt7H4D/jAXnr4Odr0JtVbHaYz5LI/Hw6RJk5g4cSJTpkxh7dq1wQ7pLO+++y5XXXVVp88JlEA2UhcBw7z2U91j3m4DFgKo6ociEgkkAsfOede4ITDzavj4cdj7D2cLi4ExcyHzKki/HKIH+vc3Mcb0SFFRUWzZsgWAFStW8MADD/Dee+8FN6geJJAliPVAuoikiUg4cD2wrMU5B4F5ACKSBUQCxW3eNSQUFjwM92yDuT+EwROhthp2/h3+9i34zzHwp6vgo8egrN0DB40xvVRFRQXx8fGA05j7ve99j/Hjx5OTk8MLL7wAwN13383DDz8MOAll1qxZNDQ0nHWfhx56iJtvvplLL72UESNG8PLLL/P973+fnJwcFi5cSG1tLQCrVq1i8uTJ5OTkcOutt3L69GkA3nzzTTIzM5kyZQovv/xy032rq6u59dZbmTFjBpMnT+bVV18N+GdyPgErQahqnYjcCazA6cK6VFW3i8jDwAZVXQZ8F/iDiNyL02B9i/q6BmpMAsz6nrOVF8CuNyD3NTjwAeS/72xv3u+UKCbeACMuhthBAfptjTFtGXn/6wG5b/4vrmzz/ZqaGiZNmsSpU6c4fPgwb7/9NgAvv/wyW7ZsYevWrRw/fpzp06cza9Ysfv7znzN9+nQuvfRS7rrrLpYvX05IyGe/R+fl5fHOO++wY8cOLrroIl566SUeeeQRPv/5z/P666+zcOFCbrnlFlatWsXYsWO56aabeOyxx7jjjjv45je/ydtvv82YMWP4yle+0nTPn/3sZ8ydO5elS5dSXl7OjBkzmD9/vn8/sHYK6DgId0zD8hbHfuz1egcws9MPGjAMLrjd2WrKYM9bTrLYvQL2rHQ28cDYhTD1Fhg9Fzw9bhoqY0w7eVcxffjhh9x00018+umnrFmzhhtuuAGPx0NKSgqzZ89m/fr1XH311fzhD39g1qxZ/OY3v2H06NGt3nfRokWEhYWRk5NDfX09CxcuBCAnJ4f8/Hx27dpFWloaY8eOBeDmm29m8eLFzJkzh7S0NNLT0wG48cYbWbJkCQArV65k2bJl/OpXvwKcXmAHDx4M5MdzXr3vr2RUPEz4krNVH4eNf4IDa2H/e7DrdWcLi4aL7oSLvwORcee9pTGmc873Tb8rXHTRRRw/fpzi4rZrsbdt20ZCQgKHDp27131ERAQAISEhhIWFNfUcCgkJoa6urkPxqSovvfQSGRkZZx0/evRoh+7nD8HuxRRYMYkw61/hn1+G+3bC/IcgPg1qT8LqR+AXw2DxBbD+CThTHexojTEBlJubS319PQkJCVx66aW88MIL1NfXU1xczOrVq5kxYwYHDhzgv/7rv9i8eTNvvPEG69at69CzMjIyyM/PZ+/evQA89dRTzJ49m8zMTPLz88nLywPgueeea7rmiiuu4NFHH6Wxln3z5s2d/I07r/eVIM6lXzJccq+zHfwIVvwAijZAcS68fh+s/BHkXAfTvg6DJkCIJ9gRG2M6qbENApxv6E8++SQej4fPf/7zfPjhh0ycOBER4ZFHHiElJYUFCxbwq1/9iiFDhvDEE09wyy23sH79eiIjI9v13MjISP74xz/ypS99ibq6OqZPn84dd9xBREQES5Ys4corryQ6OppLL72UyspKAH70ox9xzz33MGHCBBoaGkhLS+O1117z90fSLuJrm3B3MW3aNPXbehB1p53eTx8vgQKvbwr9BjlVVBNvgJRx/nmWMX3Mzp07ycrKCnYYfVJrn72IbFTVae25T98pQbQmNMIpNeRcB8W74YP/dno/lR+EtY8626AcmPhV55x+ycGO2BhjukzfThDeksbCtb8DVSjcAFufhU9fgiPb4MgDsPKHMGY+TLweMj4HYe0rchpjTE9jCaIlERg23dkW/gJ2vwlbnnNGbO9Z4WwR/WHcNZB9DaTNsS6zxpheyf6ytSU0wkkC2ddAVbFTotj6HBzeApv+7Gz9UuCyf3eqoULDgx2xMcb4Te/u5upP/ZLgwjvgW+/B//chXPqvTnKoOgp/vxt+OQKe+gIUbQp2pMYY4xdWguiIlGxnu+zf4b1fwo5Xne6yeaucDZxR2+Ovg4xFENEvuPEaY0wHWAmiM0I8TpL49jr47m5nNtlGu9+El7/hTB74wo3utOQ1wYvVmD7oXNN9Hzp0iOuuu67Va4I5vba3kSNHcvz48U6f0xlWgvCX2BT44hPw+5lQshf6D3emJi/4yBlrsfPvzrTkw6bDkMkw+98gLCrYURvTq51ruu8hQ4bw4osvfub8jk6T0VtZCcKfwiLhjjVww/NOqeK2FfCdTTBkCkQnOtOS73sX1vwG/jMdXv027F4JdWeCHbkxvZ73dN/5+fmMHz8egD/96U9cffXVzJ07l3nz5p11zfr165k8eTJ5eXm89957TJo0iUmTJjF58uSmEdCN8vPzyczM5JZbbmHs2LF87Wtf46233mLmzJmkp6fz8ccfA1BaWsq1117LhAkTuPDCC/nkk08AKCkp4fLLL2fcuHF84xvfwHsQ89NPP82MGTOYNGkS3/rWt6ivrw/Y5+TNShD+FhbltDs0ShgNt7/jvC7Lh9fuRQ9/gpw8Dpufhs1Po+ExNGT8EyE5X0BGXWa9oUzv81D/AN33RJtvn2u675Y2bdrEJ598wsCBA3n33XcBWLt2Ld/5znd49dVXGT58OPfccw+LFy9m5syZVFVVtTr9xt69e/nrX//K0qVLmT59Os8++yxr1qxh2bJl/Md//AevvPIKDz74IJMnT+aVV17h7bff5qabbmLLli385Cc/4ZJLLuHHP/4xr7/+Ok888QTgjIp+4YUX+OCDDwgLC+Nf/uVfeOaZZ7jppps699n5wBJEAJ08U8e2whP8Y8dRPik6wanaej4p/AYAo6WIK0PWca1nDaPOHMGz7XnY9jwV9GND1EXkxs+lNGUmSQP6caq2gTkZSWQOjiUi1OaIMsZX55ruu6UFCxYwcGDzSpQ7d+7k9ttvZ+XKlQwZMgSAmTNnct999/G1r32NL3zhC6Smpn7mPmlpaeTk5AAwbtw45s2bh4g0TQMOsGbNGl566SUA5s6dS0lJCRUVFaxevbppAaErr7yyqbSzatUqNm7cyPTp0wEn6SUnd82sDpYg/KzqdB3Ltx3mpY2FrM8vpeEcU13l6VD+p/4L/DH0Oi7RrWTrbhbwMZkhBcyt+Qdza/5BeVEMqxqm8F79RH731jTqQiJIS4whe3AcE1L7M3n4AMYN6U9kmCUN082d55t+V2hruu+YmJiz9gcPHsypU6fYvHlzU4K4//77ufLKK1m+fDkzZ85kxYoVZGZmnnVd4zTg4Ez97T0teGemAb/55pv5+c9/3qHrOyOgCUJEFgK/xVlR7nFV/UWL938DXObuRgPJqjogkDEFykf7SvjGkxuoa2jgVK2zRKEnRMgeFMuJmlouGZPIF6emcqCkmoxBsSzfdoT46DBuuySNUI9TJVV9uo7CA5+iO16hf97fGVC5ly963ueLnvcBeLX+YpYfn8E/jk1g2VaneBsaImQMimXSsAFMGxnP9JEDSY2PDs6HYEw35j3d98mTJ9s8d8CAATzxxBMsWLCAmJgY5syZQ15eHjk5OeTk5LB+/Xpyc3M/kyB8cemll/LMM8/wox/9iHfffZfExETi4uKYNWsWzz77LD/84Q954403KCsrA2DevHlcc8013HvvvSQnJ1NaWkplZSUjRozo0OfQHgFLECLiARYDC4BCYL2ILHNXkQNAVe/1Ov87wORAxRMoZ+oaePTtPfzu3Tzq3eLCjLSBXDcllUU5g4iNDDvr/BlpTjF2QuqAz9wrJiKUmLGTYOwk0AeheBdsehI++h0A13jWco3H6ab3Sf/LeK12Os+VZbD9kLL9UAXPrGtefapfRCgPXT2Oi0cnMGSA9ZYyfdO5pvv2RUpKCq+99hqLFi1i6dKlPP3007zzzjuEhIQwbtw4Fi1adP6btOKhhx7i1ltvZcKECURHR/Pkk08C8OCDD3LDDTcwbtw4Lr74YoYPHw5AdnY2P/3pT7n88stpaGggLCyMxYsXd0mCCNh03yJyEfCQql7h7j8AoKqtlpNEZC3woKr+o637+nW6705QVf66sZDvv/hJ07HMQbE8duNU0hJj2riyg0r3wfa/waqHW317Tfr3eKlmOm8VQuWps4uyo5JiuGRMIjPHJHLx6ITPJC1jAsGm+w6enjDd91CgwGu/ELigtRNFZASQBrTaxUBEbgduB5qyajCVnzzDD1/5lNc+OQxAanwUv/7ypKbSQUAMHAWXfhcuvgtOFMCuN5yxFQc/BOCSPf/JJQiaOp2ilMt4sXoin55O5qN9pewrrmZfcTV//vAAoSHCtJHxzMlIZk5GEhkpsU3LJRpjjLfu0kh9PfCiqrbauVdVlwBLwClBdGVgLZ2oqeUr//cRu45WEh3u4e556Xx9ZhrhoV00pMQT5iSLi77tbMdy4a2HQOth33tI4cekFn7MPQAJ6dTP/Bx74mexonwYq/eWsPlgGR/tK+WjfaX84o1cBsVFMntsEnMykrgkPdFKF8aYJoFMEEXAMK/9VPdYa64Hvh3AWPzidF09dzy1kV1HKxmdFMMTN09nZCCqk9ojORO++rwbYJUzF1Tucmeqj5I9eNb+lkx+S2ZMMndnLKLqkstZXTeOt/dW8N7uYo5UnOKFDQW8sKGAMI9w4agEFmSnMD8rxdoujOnjAtkGEQrsBubhJIb1wFdVdXuL8zKBN4E09SGYYLVBNDQo9/1lC69sOURSbAR/+5eLu3dvofpap/opdznkvg4nmhuwCYuBMfNoyPgcu2Iv5u2Dtby76xgbD5Sd1S133JA4FmSnsCA7hezBcVYVZdpl586dZGZm2r+bLqaq5Obm+qUNIqBrUovI54D/xunmulRVfyYiDwMbVHWZe85DQKSq3u/LPYOVIH75Zi6PvZtHTLiHF751EeOHBmhkaCCowtFP3WTxGhxpblhHPDDiYshYRHnqXN46FstbO46yek8xJ8801/gNHRDF/KxkFmQP4oJRAwnz2Cwtpm379+8nNjaWhIQESxJdRFUpKSmhsrKStLS0s97rdgkiEIKRIJ766AA/euVTPCHC0lumM3tsUpc+3+/K3Ubu3NfgwAfQ4NXraeAoSL+CM6MX8EHtWFbuKuetnUcprjzddEpsZChzMpJZkJ3CnIwk4qzdwrSitraWwsJCTp06FexQ+pTIyEhSU1MJCzv7/6UliABYtfMo3/zzBhoUHrluAl+eNuz8F/UkNWWw5y1nKdW9bzn7jcL7wag5NKRfzvaYC3kjX3lr51F2H61qOiU0xKvdIjuFodZuYUy3ZAnCz45XnWb+r9+j/GQt98xP5575Y7vkuUFTXweF651ksXslHNt+9vuDJ0L6FRxOmc3rJYP4x87iz0wnkj24ud1i3BBrtzCmu7AE4Wd3PbeZZVsPccmYRJ66bUbf+2NXXgB7VjrbvvegzmvBo5gkGLOAqhFzead2HMv31PDe7rPbLYb0j2S+mywuSEvouq7AxpjPsAThR1sLyrlm8QdEhoWw8p7ZDE/oxj2WukJtDex/v7l04d0rSkJgyBTqRs5mW+RkXj42lBW5pRzzbreICGV2RhI5Q/vzpWnDGBhjU5ob05UsQfjRrX9az9u5x/jW7FE8sMimCziLqrMG9+4VTumi4GNoqG1+PywGHTGTQwkXsvJUFs/vj2HXsaqzbnHJmERrtzCmC1mC8JPG0kNUmIc1/3YZCf0izn9RX3a6yukNte9dyHsHinee/X7sYKqGXsLKmix+viuFYuLPetvGWxgTeJYg/MRKD51UcdhJFvvecX5WHT3r7frELPbFTefNk1n86dAQSs40d8cb0j+SeVlOyeLCUQNtgSRj/MQShB80lh6iwz28/30rPXSaKhzb4ZQs9r0D+R+c1ditnnDKEybzccgEni0ew/vVQ2lwl0qPCfcwa2wS87NSuCwz2dotjOkESxB+0Fh6uGP2aO5f1P7FQMx51J2GgnXN1VGHNgPN/wbrwvuzr99U3qzJ5K/l6RRoCgAhAlNHxDuli6wURifFWFWUMe1gCaKTth86wZX/s8ZKD13pZCnsX+2ULvLegfIDZ71dGZXKxyET+duJdFbXZVNBPwBGJkQzPyuFeVkpTB8ZT6hN/WFMm7rbehA9zlMfOn+cvjJ9mCWHrhI9EMZd62zgLIzUWB21fzWxNYXMo5B5oaChIRyMzGDlqSxWlWXz5JqxPL5mP/2jwpiT4VRFzbapP4zxGytBuCpO1XLBz1ZRU1vPqu/OZnRSP78/w7RTQz0c2gL73oa8d52qKa/utGdCItks41h5Kov3G3LYramEhoQwI20gczOTmZORbFVRxrisBNEJr209TE1tPReNSrDk0F2EeCB1qrPN+p7bnXZtU3VUePFOLmAjF4RtBKA8JJ4PajP4MD+L5/dl8dPXhzJ8YAyXZSQxJzOZi0YlEBlmvaKM8ZUlCNcrm521jL44NTXIkZhziugHYy93NnC60+5/r6lKakDVUa70fMSVno8AKCWODyszWfdxFj//KJuDnmFcNDqpqXQxbGAfHx1vzHlYFRNw+EQNF/38bSJCQ9jww/m27GZPpArH98CBNZC/xulOW3XkrFNKtR/rGrKatrrELC7LGsScjCSmjRhoc0WZXq3bVTGJyELgtzgLBj2uqr9o5ZwvAw/h9HXcqqpfDWRMrVm53RnINScjyZJDTyUCSWOdbdqtTsIoyTsrYQysPMQiz3oWedYDUFbRj/UfZvDWB9n8OnQ8SaOnMCdrELPHJjOof2SQfyFjgi9gCUJEPMBiYAFQCKwXkWWqusPrnHTgAWCmqpaJSHKg4mnLm5863zQXjh8UjMebQBCBxDHONvUWJ2GU7nOmBMlfg+avIb6iiMs9G7nc47RhnMiL5uM9mTyhWZwYsZCLpk0hNT6a8UP6ExVubRem72kzQbh/5P+sql/rwL1nAHtVdZ97r+eBa4AdXud8E1isqmUAqnqsA8/plNLqM6zbX0JoiDA3M6WrH2+6iggkjHa2KTchqs6YC7d0UbdvNf0rC1ng2cQCNkHRM+QVDKZIE/mjjKEqeRqJWZdwQVYaWYPiCAmxnlGm92szQahqvYiMEJFwVT3TznsPBQq89guBC1qcMxZARD7AqYZ6SFXfbOdzOuWtHUdpULgkPZH+UVa91GeIQPxIZ5t8o/MfoewAHFjL6dyVhOx+g9EcZjSHmcU2OP43GlYLu94bxoueLGoGzyAxaxbTJk4gJc6qo0zv5EsV0z7gAxFZBlQ3HlTVX/vp+enAHCAVWC0iOapa7n2SiNwO3A4wfPhwPzy22ZvbneqlRVa9ZOJHQPwIIibd4Kx/ceADOFlKTcFWavI+IK5sG1l6kCw9CIdWwCEoeiuBd8LHUTN4BsnZlzJu0kVERdogS9M7+JIg8twtBIhtx72LAO8FnFPdY94KgXWqWgvsF5HdOAljvfdJqroEWAJOL6Z2xNCm03X1rM07DsC8rKA0f5juKiwKxswHIGrCl4kCqK1BizZRunM1J/M+YGDpZoY2lDC0djUcXA0Hf0X1GxFsj8jk1KCpJGRewrAJs/H0Swzqr2JMR503QajqTwBEpJ+7X9X2FU3WA+kikoaTGK4HWvZQegW4AfijiCTiVDnt8/H+nbblYDmnahvISIklOdaqCcx5hEUhI2eSMHImCQANDdQe2UHB1rc5lfcB8aVbGNxwhHFntsLBrXBwKayEo2FDqUicTOyYi0nOnkVISrYzCNCYbu68CUJExgNPAQPd/ePATaq6va3rVLVORO4EVuC0LyxV1e0i8jCwQVWXue9dLiI7gHrge6pa0qnfqB3W5jmPumh0Qlc90vQmISGEDRnPqCHjgbsAKD9WyJ5N71C9dy3xpVsYW7+XlNoiUg4XweHX4H2oCYmmPD6HiLQLiR87E0md7sxJZUw3c96BciKyFviBqr7j7s8B/kNVLw54dK3w50C5L//+Qz7OL2XJP0/l8nHWBmH8r6D4BLlb11K5Zy0xxZvJrstlWEjxZ86riEnDM3wGMaMvgmEzICnTShnGrwI1UC6mMTkAqOq7IhLT7ui6mZoz9WwuKCNE4IJRVoIwgTEsqT/D5i+C+YtQVfJLTvLS9p2U5q4h4shGsupzmSD7iaveDzv3w84XAKjzRNEwaBLhI6bB0KnO1n+Y0/vKmC7iUy8mEfkRTjUTwI10YTtBoGw6WEZtvTJ+aJx1bzVdQkRIS4whbfY0mD0NVWX30Spe2HOYwp0f4zm0nnH1uUyUPIZRDEUfOpurPjoRT+o0GDLFTRpTrGrKBJQvCeJW4CfAyzjTYbwPfD2QQXWFzQfLAJg2wv6DmeAQETIGxZIxKBYuHUt9w9fYebiCN/NK2LZ7D7UHN5LZsIeJksfEkDziTx6H3W86m0vj05DGEsbQqTB4gtMDyxg/8CVBzFfVu7wPiMiXgL8GJqSusflgOQCThw8IahzGNPKECOOH9mf80P4waxS19QvYVnSCD/NK+OO+Eo4dzCW9djcTQpyEkSP7iSzbD2X74dMXAVDxICnZzQlj6FRrzzAd5ksj9SZVnXK+Y13FH43UqsrUn75FafUZ3v/+ZTbts+kR6uobyD1Sybr9pazfX8qm/cdIrNnHxJC8plLG2JBCPLT4Px0W45QsBk+CIZOcn4npljT6GL82UovIIuBzwFAR+R+vt+KAuo6F2D0cKDlJafUZEvuFkxpvxXHTM4R6QppKGLddkoaqkldczcf7S1mXX8qj+0spLS9jvOQzMSSPSe6WWlsMBz90tkZh0TAop0XSGAseWyLGNGvrX8MhYANwNbDR63glcG8ggwq0zQVO+8OkYfG2HKXpsUSEMcn9GJPcj69e4ExBU1h2kvX5pXy8v5Tf7C8lr7iaBE6QE7KfcZJPTsh+Jofmk1Jb7CzhWrCu+YahUTBo/NlJIykDPNaJo686Z4JQ1a3AVhH5G1CtqvXQNMNrj55sZtOBcsDaH0zvkxofTWp8NJ+f7KyMeLzqNOv3l7LxwCTWHixjSdEJamuUgVQwLiSfHNnP1PADTPTkk1h3BArXO1uj0EhIGeeVNCZCUhaEhgfl9zNdy5fy5EpgPtA4xUaUeywoA+X8YVvRCQAmDRsQ3ECMCbDEfhEsyhnMopzBAJyqrWdb0Qk2Hihj44F0XjhQxu+qnYmaB1DJuJB8JnnyuTi6gGz2E3+6CIo2OlsjTzgkZztVVIMmOD9TxkFkXDB+RRNAviSISO/5l1S1SkR6bKtuQ4Oy60glAFmD7R+06VsiwzxMHzmQ6SOd7t2Ng/echFHGpgND+N2xHBY736GIo4pxIQe4OKqAC6MOkl6/jwE1B+HwFmfzFp/mlTTGO6/jhtrgvh7MlwRRLSJTVHUTgIhMBWoCG1bgHCw9SU1tPSlxEQyMsWKy6duaBu8lxnDdVKda6kRNLZsPlrHpQBkbD5bxScEAPqwe1zTZfywnyQ45wOy4I0yLLGJ0w37iq/IIaexyu3NZ8wOi4s8uaQzKcRvDrV2jJ/AlQdwD/FVEDgECDAK+EsigAin3SAUAmYOs9GBMa/pHhTEnI5k5Gc4U+A0Nyr7jVWwpOMHWgnK2Fpaz8VAM68qzmq4JpY6s0CMsGHiM6ZFFjNH9DKzchaemDPavdrZGnnBnbIZ30hg0HiL7d/Wvas7Dl+m+14tIJpDhHtrlrt/QI+087FQvZQ5uz9IWxvRdISHCmORYxiTHNpUyTtXWs+NwBVsLytlSUM7WgnK2lYSy7Vgq0DhESsmIqmRhYjEXRDlJI6FyF57yfDjyibN56z8cUrKd9o2Ucc7PxHQrbQSRr52eM4BsIBKYIiKo6p8DF1bg7D7qJohBliCM6ajIMA9ThsczZXh807Gy6jN8UnSCLQedUsbWgnJ2VQu7CuKA0cAsADLjYWHicaZHHWJMw34Sq3bhKd4JJw46m9dUIoSEOV1tk7Od5JEy3nkdN8TaNrqAL+tBPIizJGg2sBxYBKwBemSCyC85CUBaYr8gR2JM7xIfE87ssUnMHpsEOA3ghWU1TSWMrYXlbCs6QW5ZA7lliUAiMAGA0QMjmJ1UyYUxR8gMKWDQqX2El+yEsnw4+qmzbfN6WGR/SB7nJo1xzuvkLOtJ5We+lCCuAyYCm1X16yKSAjwd2LACQ1UpKHUSxAibXsOYgBIRhg2MZtjAaP5p4hDAmS5k99EqthQ4yWL7oRPkHq4kr/Q0eaXhLGU4MByYydABUUxNC+XS/sfJCStieF0+0WW7nGRRUwYH1zqbt5bVVEkZkJAOYbZiZEf4kiBqVLVBROpEJA44xtlrTZ+TiCwEfouzotzjqvqLFu/fAvwnzWtV/6+qPu5r8O1VWn2GqtN1xEaGMiDa6jWN6WqhnhCyh8SRPaT5m35tfQN7jlbxadEJPj10gk+LTrDjcAVF5TUUlcMyIoBRwChS4hYxfnAcM5JqmR59mDF6gNiK3cjR7VC8q/VqKgmBgaOchvGkTKekkZQJCWMscZyHLwlig4gMAP6AM+VGFfBhm1fQNOJ6MbAAKATWi8gyVd3R4tQXVPXOdkXdQQfc0sPwgdE2xYYx3USYV9L4svvds75B2VdcxbaiE3xaVMGnh06w41AFRytOc7SimFW7wBmzm0lCzATGDb2NnGnRXBBXTpangMSqvcjxXCjOhdJ9ULLX2XJfa36wd+JoTBpJmU7DeGiPnizCb9qarG+mqn4A3Kuqp4Hfi8ibQJyqfnKu67zMAPaq6j73fs8D1wAtE0SXOei2P4xIsOolY7ozT4iQnhJLekosX3A7RTU0KAdKT55V0vi0qIKS6jOs3l3M6t3ON1KIIzp8OhmD5pI5NI7xk8OZGF3MKC0gunyvkzSO7XTGbLSaODxu4sjo84mjrRLE/wBTcUoLUwBUNb8d9x4KFHjtFwIXtHLeF0VkFrAbJxkVtHKOXxwoaSxB9PgVU43pc0JCmgf1NbZpNDaEbz/kJIttRSfIPeKUNDYfLG9a98URx9ABM8kavJDMjDiyk8MZH3GMoWfy8ZTsbpE49jhba4kjORMSM5wBf4npzs+I3tnppa0EUSsiS4DUFtN9A9ByEaEO+jvwnKqeFpFvAU8Cc1ueJCK3A7cDDB8+vMMPO1hqJQhjehPvhvCF4wc3HS+pOs2uI5XsPFJJ7uEKco9UsvtopduuUcNbO481nRsRGsfYlMvIHHQ1mVPiyE4KJzv8KP0r85yk0Vri4O9nBxI31E0YbtJIchNIv5Qe3R23rQRxFc4kfVdw9nTfviri7MbsVJobowFQ1RKv3ceBR1q7kaouAZaAs2BQB2IB4GCpM1eA9WAypndL6BfBxWMiuHhMYtOxuvoG8ktOknukgtzDleQeqWDnYSdpbCs60TSJZ6Pk2HgyB19B1qDryMyIJSsxnNEhhwkr2Q3HvbaSvVBR5Gz73jk7kIj+zaWMpLHNSSQ+rUesvdHWdN/HgedFZKc79Xd7rQfSRSQNJzFcD3zV+wQRGayqh93dq4GdHXiOz5qqmKwEYUyfE+oJaVo/46oJzcdP1NSy+6hT0tjhJo5dRyo5VnmaY5XFrN5d3HyPEGFMcgqZg8aQNfh6MibGkpEczaCGo8jxPW7S2AXH9zi9qk6VQ9EGZ/MWEua2c4x1uuEmpju9qhLGQPTArvlAfODLVBsdSQ6oap2I3AmswOnmulRVt4vIw8AGVV0G3CUiV+OsUFcK3NKRZ/mi5kw9xypPE+YRBve3VeSMMY7+UWFnzXALToN4QdlJdroJo7HEcaD0JLlHKsk9UskrWw41nR8bGUp6cixjU+YwJvkqxmbFMja5HymhlW7i8Eoax/c4XXGP73K2lqLim5NFwujm1wNHQXjXtp+ed03q7qaja1LvOlLJFf+9mrTEGN751zn+D8wY0+tVn65zShtu28bOw5XsPlZJ+cnWp6eLjQxlbEosY1P6kZ4cS3pKP8amxJIcUYeU5jnJ4vhuKMlr7lF1pqrVewFOW4d30mjcBgw/75xVfl2Turc56DUGwhhjOiImIpTJw+OZ7DUHlapyvOoMe446DeG7j1Wx92hVU+JoXGvDW5ybONJTMkhPnsqYCf0YndyPwbERhJwsbk4WJXubk0fpvua2Du/ZcQFCQiF+ZOslj9jBHW4o92Uupgjgi8BI7/NV9eEOPTFIDpS4DdTW/mCM8SMRISk2gqTYsxvFVZXiqtPsOVrF7qOV7DlW5SaRKk7U1LLhQBkbWiSOqDAPo5JiGJ3Uj9FJFzMqZQGjx/djVFIMkSEKJwrOLm00JpATBc37LYXFOAmjA3wpQbwKnMDpyXS6Q0/pBqwEYYzpSiJCcmwkybGRzGyZOCpPs+dYc+LYV1xFXnE1xZWn2X6ogu2HKlrcC4YOiHITxxBGJ6czOqsfo5P6kdgvHKk7dfaIce8kcrLks1Or+8iXBJGqqgs7dPdu5EDTKGobJGeMCR4RITkukuS4sxMHOD2qGpNFXnEVeceqyCuu4kDJSQrLaigsq+E9r15V4FRXjU7u5yaP8YxOuoDRmf0YPjCaME8InCx1ksdPprc7Vl8SxFoRyVHVbec/tfuyEoQxprvrHxX2mTYOcCY0PFh60k0YbvIormLvsSoqTtW1Mmrc6ZI7IiHaSRzJHRvp7UuCuAS4RUT241QxCaCqOqHty7qP+galsMwShDGmZwrzhLglhLP/0Dc2kDcmjLxjzcmjqLzGTSbVsONoh57rS4JY1KE7dyOHT9RQW68kx0YQFe4JdjjGGOMX3g3kF45KOOu9mjP17D/enDDu/WX77+/LQLkDIjIRuNQ99H5HB88Fi83iaozpa6LCPWetvXFvB+4Rcr4TRORu4Bkg2d2eFpHvdOBZQdO8DoQ1UBtjjK98qWK6DbhAVasBROSXOFOAPxrIwPzpgJUgjDGm3c5bgsBplK732q93j/UYjbO4WgO1Mcb4zpcSxB+BdSLyN3f/WuCJgEUUAE1dXK0EYYwxPvOlkfrXIvIuTndXgK+r6uaARuVHqtpcxWQlCGOM8Vlba1LHqWqFiAwE8t2t8b2Bqloa+PA6r/J0HZWn6ogO9zAwJjzY4RhjTI/RVgniWZxV5TYC3nOCi7s/KoBx+c2xilMApMRFIj146T9jjOlqba0od5X7M63rwvG/YxXO/IJJsRFBjsQYY3oWX8ZBrPLl2DmuXSgiu0Rkr4jc38Z5XxQRFZF2LWbhi2OVToJItgRhjDHt0lYbRCQQDSSKSDzNXVvjgKHnu7GIeIDFwAKgEFgvIstUdUeL82KBu4F1HfoNzuNYZXMVkzHGGN+1VYL4Fk77Q6b7s3F7FfhfH+49A9irqvtU9QzwPHBNK+f9P+CXwKl2xO2zoxVWgjDGmI44Z4JQ1d+67Q//qqqjVDXN3Saqqi8JYihQ4LVfSIuSh4hMAYap6usdCd4XTVVMcZYgjDGmPXwZB/GoiIwHsoFIr+N/7syDRSQE+DVwiw/n3g7cDjB8+PB2PaexF1NyrFUxGWNMe/jSSP0gzrxLjwKXAY8AV/tw7yJgmNd+qnusUSwwHnhXRPKBC4FlrTVUq+oSVZ2mqtOSkpJ8eHSzYmukNsaYDvFlLqbrgHnAEVX9OjAR6O/DdeuBdBFJE5Fw4HpgWeObqnpCVRNVdaSqjgQ+Aq5W1Q3t/SXacrSxBGGN1MYY0y6+JIgaVW0A6kQkDjjG2SWDVqlqHXAnsALYCfxFVbeLyMMi4ksJpNOqT9dRfaaeiNAQ4iJ9mXbKGGNMI1/+am4QkQHAH3B6MVXhTPd9Xqq6HFje4tiPz3HuHF/u2R6N1UtJsRE2itoYY9rJl0bqf3Ff/l5E3gTiVPWTwIblH6UnzwCQ0M/aH4wxpr18qncRkaHAiMbzRWSWqq4OZGD+UFrlJIiB0WFBjsQYY3qe8yYIdwW5rwA7aF44SIHunyDcEkS8zeJqjDHt5ksJ4logQ1VPBzgWvyurbixBWIIwxpj28qUX0z6gR9bRWAnCGGM6zpcSxElgizuDa1MpQlXvClhUftJUgrAEYYwx7eZLgliG1wC3nqS0uhawBGGMMR3hSzfXJ7sikEAoO2klCGOM6ai21oP4i6p+WUS2cfaSowCo6oSARuYHpW4VU7w1UhtjTLu1VYK42/15VVcEEgil1gZhjDEd1tZ6EIfdnwdwGqcnAhOA0+6xbq2uvoETNbWIQP+oHtkJyxhjgsqX6b6/AXwMfAFnZtePROTWQAfWWeU1TgP1gKgwPCE2D5MxxrSXL72YvgdMVtUSABFJANYCSwMZWGc1dnG1MRDGGNMxvgyUKwEqvfYr3WPdWqmNojbGmE5pqxfTfe7LvcA6EXkVpzfTNUC3n83VurgaY0zntFXFFOv+zHO3Rq8GLhz/KbEeTMYY0ynnTBCq+pPO3lxEFgK/BTzA46r6ixbv3wF8G2eW2CrgdlXd0dnngrVBGGNMZ/ky3fc04Ad4rQcB5x8oJyIeYDGwACgE1ovIshYJ4FlV/b17/tXAr4GF7f0lWtM4zUa8rQVhjDEd4ksvpmdwejJtAxrace8ZwF5V3QcgIs/jtF80JQhVrfA6P4ZWRmx3VMUpJ0HYGAhjjOkYXxJEsap2ZLK+oUCB134hcEHLk0Tk28B9QDgwt7UbicjtwO0Aw4cP9+nhlW6CiIu0BGGMMR3hSzfXB0XkcRG5QUS+0Lj5KwBVXayqo4F/A354jnOWqOo0VZ2WlJTk030rauoAiLUEYYwxHeJLCeLrQCbOokGNVUwKvHye64qAYV77qe6xc3keeMyHeHxSedotQUT5tOy2McaYFnz56zldVTM6cO/1QLqIpOEkhuuBr3qfICLpqrrH3b0S2IOfNJYgrIrJGGM6xpcEsVZEstvb/VRV60TkTmAFTjfXpaq6XUQeBja47Rp3ish8oBYoA25uZ/zn1NhIHRtpJQhjjOkIX/56Xoiz5Oh+nFldBVBf1oNQ1eXA8hbHfuz1+u7PXOQHqkrlKWuDMMaYzvAlQfhlXEJXOnmmnvoGJTIshPBQX9rhjTHGtHTev57u2g8DgH9ytwHdfT2IxtKDtT8YY0zH+bIexN04g+WS3e1pEflOoAPrjMb2hzgbJGeMMR3mSxXTbcAFqloNICK/BD4EHg1kYJ1RUWMN1MYY01m+VNALzmR6jerdY92WVTEZY0zn+fIV+48460H8zd2/FngiYBH5gXVxNcaYzjvvX1BV/bWIvAtc4h76uqpuDmhUnVTRWIKwNghjjOkwX6b7vhDYrqqb3P04EblAVdcFPLoOamyDsComY4zpOF/aIB7DWcynURV+nDMpEKyKyRhjOs+nRmpVbVqnQVUb8K3tImgqrYrJGGM6zZcEsU9E7hKRMHe7G9gX6MA6o7mKqVvnMWOM6dZ8SRB3ABfjzMjauOjP7YEMqrOa52GyBGGMMR3lSy+mYzhTdfcY1aedBNEvwqqYjDGmo3yZamOsiKwSkU/d/Qki0urKb91FVVOCsBKEMcZ0lC9VTH8AHsBZswFV/YRuXqKwKiZjjOk8XxJEtKp+3OJYXSCC8ZfqM054MVaCMMaYDvMlQRwXkdE461AjItcBh325uYgsFJFdIrJXRO5v5f37RGSHiHziVmONaFf0rVBVqk41JghPZ29njDF9li8J4tvA/wGZIlIE3IPTs6lNIuIBFgOLgGzgBhHJbnHaZmCauzrdi8AjvofeutN1DdQ1KOGeECJCLUEYY0xH+bJg0D5VnQ8kAZnAbJrnZWrLDGCve/0Z4Hngmhb3fkdVT7q7HwGp7Qm+NU0N1Nb+YIwxnXLOBOHOufSAiPyviCwATgI3A3uBL/tw76FAgdd+oXvsXG4D3jhHLLeLyAYR2VBcXNzmQxu7uFr1kjHGdE5bX7OfAspwFgf6JvADnHUgPq+qW/wZhIjcCEzDKZ18hqouAZYATJs2TVs7p9GJxsWCbAyEMcZ0SlsJYpSq5gCIyOM4DdPDVfWUj/cuAoZ57ae6x84iIvNxks9sVT3t473P6VB5DQBDBkR29lbGGNOntdUGUdv4QlXrgcJ2JAeA9UC6iKSJSDjO2Ill3ieIyGScBvCr3RHbnVZY5iSI1Phof9zOGGP6rLZKEBNFpMJ9LUCUuy+AqmpcWzdW1ToRuRNYAXiApaq6XUQeBjao6jLgP4F+wF9FBOCgql7dmV+oOUFEdeY2xhjT550zQahqp1t5VXU5sLzFsR97vZ7f2We0VFDqdIqyEoQxxnSOL+MgehQrQRhjjH/0qgShqhSWOSWIYVaCMMaYTulVCaL8ZC3VZ+qJjQglLsoGyhljTGf0qgRR4JYeUgdG4zZ6G2OM6aBelSCs/cEYY/ynlyWIxh5MliCMMaazelmCsEFyxhjjL70yQQyzEoQxxnRar0oQNkjOGGP8p9ckCGcMhFOCGGolCGOM6bRekyBKq89QU1tPXGQo/aNsqm9jjOmsXpMgmtofBlr1kjHG+EOvSxDWxdUYY/yj1ySIplHU1kBtjDF+0WsShA2SM8YY/+pFCcIGyRljjD8FNEGIyEIR2SUie0Xk/lbenyUim0SkTkSu68yzmhuprQRhjDH+ELAEISIeYDGwCMgGbhCR7BanHQRuAZ7tzLO814EYOsAShDHG+EMgF02YAexV1X0AIvI8cA2wo/EEVc1332vozIOOV53hVG0DA6LDiI20MRDGGOMPgaxiGgoUeO0XusfaTURuF5ENIrKhuLj4M+9bA7Uxxvhfj2ikVtUlqjpNVaclJSV95v3mSfqsgdoYY/wlkAmiCBjmtZ/qHvM7GyRnjDH+F8gEsR5IF5E0EQkHrgeWBeJBNkjOGGP8L2AJQlXrgDuBFcBO4C+qul1EHhaRqwFEZLqIFAJfAv5PRLZ35FlWgjDGGP8LZC8mVHU5sLzFsR97vV6PU/XUKYVWgjDGGL/rEY3UbVFViqwEYYwxftfjE0Rx5WlO1zUwMCacmIiAFoiMMaZP6fEJosBKD8YYExA9PkHYIDljjAmMXpAgbJCcMcYEQq9JEFaCMMYY/+oFCcK6uBpjTCD0ggRhJQhjjAmEHp0gGhq8x0BYCcIYY/ypRyeI4qrTnKlvILFfOFHhnmCHY4wxvUqPThAFpe4qclZ6MMYYv+vRCcLaH4wxJnB6eIKwQXLGGBMoPTxB2CA5Y4wJlF6RIKwEYYwx/tejE4StJGeMMYET0AQhIgtFZJeI7BWR+1t5P0JEXnDfXyciI329d32DcqjcShDGGBMoAUsQIuIBFgOLgGzgBhHJbnHabUCZqo4BfgP80tf7H6s8RW29khQbQWSYjYEwxhh/C+QKOzOAvaq6D0BEngeuAXZ4nXMN8JD7+kXgf0VEVFXPddOSqtNc8B9vUVfvnGKlB2OMCYxAVjENBQq89gvdY62eo6p1wAkgoeWNROR2EdkgIhsqq6o5WnGakuozAMwcnRiI2I0xps/rEWt0quoSYAnAlKlTdfkD8wDwhAhJsRHBDM0YY3qtQCaIImCY136qe6y1cwpFJBToD5S0ddMQEQb1j/RnnMYYY1oRyCqm9UC6iKSJSDhwPbCsxTnLgJvd19cBb7fV/mCMMabrBKwEoap1InInsALwAEtVdbuIPAxsUNVlwBPAUyKyFyjFSSLGGGO6gYC2QajqcmB5i2M/9np9CvhSIGMwxhjTMT16JLUxxpjAsQRhjDGmVZYgjDHGtMoShDHGmFZJT+tVKiKVwK5gx9FNJALHgx1EN2GfRTP7LJrZZ9EsQ1Vj23NBjxhJ3cIuVZ0W7CC6AxHZYJ+Fwz6LZvZZNLPPopmIbGjvNVbFZIwxplWWIIwxxrSqJyaIJcEOoBuxz6KZfRbN7LNoZp9Fs3Z/Fj2ukdoYY0zX6IklCGOMMV3AEoQxxphW9agEISILRWSXiOwVkfuDHU+wiMgwEXlHRHaIyHYRuTvYMQWTiHhEZLOIvBbsWIJNRAaIyIsikisiO0XkomDHFAwicq/7f+NTEXlORPrUIjIislREjonIp17HBorIP0Rkj/sz/nz36TEJQkQ8wGJgEZAN3CAi2cGNKmjqgO+qajZwIfDtPvxZANwN7Ax2EN3Eb4E3VTUTmEgf/FxEZChwFzBNVcfjLDfQ15YS+BOwsMWx+4FVqpoOrHL329RjEgQwA9irqvtU9QzwPHBNkGMKClU9rKqb3NeVOH8EWq733SeISCpwJfB4sGMJNhHpD8zCWWcFVT2jquVBDSp4QoEod6XKaOBQkOPpUqq6GmeNHW/XAE+6r58Erj3ffXpSghgKFHjtF9JH/yh6E5GRwGRgXZBDCZb/Br4PNAQ5ju4gDSgG/uhWuT0uIjHBDqqrqWoR8CvgIHAYOKGqK4MbVbeQoqqH3ddHgJTzXdCTEoRpQUT6AS8B96hqRbDj6WoichVwTFU3BjuWbiIUmAI8pqqTgWp8qEbobdy69WtwEuYQIEZEbgxuVN2Lu7Tzecc49KQEUQQM89pPdY/1SSIShpMcnlHVl4MdT5DMBK4WkXycKse5IvJ0cEMKqkKgUFUbS5Mv4iSMvmY+sF9Vi1W1FngZuDjIMXUHR0VkMID789j5LuhJCWI9kC4iaSISjtPotCzIMQWFiAhOPfNOVf11sOMJFlV9QFVTVXUkzr+Ht1W1z35TVNUjQIGIZLiH5gE7ghhSsBwELhSRaPf/yjz6YGN9K5YBN7uvbwZePd8FPWY2V1WtE5E7gRU4vRKWqur2IIcVLDOBfwa2icgW99i/u2uAm77tO8Az7peofcDXgxxPl1PVdSLyIrAJp8ffZvrYlBsi8hwwB0gUkULgQeAXwF9E5DbgAPDl897HptowxhjTmp5UxWSMMaYLWYIwxhjTKksQxhhjWmUJwhhjTKssQRhjjGmVJQhjjDGtsgRhejx3zqGrRORbInJYRLZ4bTnuOb8XkZltXe/Dc0aKSI3X2BNERL1Hb4tIqIgUn2/qcXe69itaHLtHRB7z2j9nzL4QkSj3MzgjIokdvY/puyxBmN5gMrAFyAF+qKqTvLZt7jkXAh+d53pf5KnqJK/9amC8iES5+wvwbQqY5/jsFNTXu8cbtRXzealqjRtrn5rJ1PiPJQjT44jIWBFZIyLbROQHwCBVLQQm0MofehHJAnaran1b17vf6he45/xURB71MaTlOFOOA9zA2X/kEZEbReRj99v8/7lrm7wIXOmOeG6clXcI8H7LmEXkJhH5RES2ishTbkkmV0T+JCK7ReQZEZkvIh+4i8HM8P3TNObcLEGYHkVEIoC/Afepag7OlO+57tvjcKa6bqxeut09vgh404frHwR+ICJfwylV3ONjWM8D17urlk3Aa+p19w/9V4CZ7rf5euBrqloKfOzGBk7p4S/aPLXBIuBNERkH/BCYq6oTcRZHAhgD/BeQ6W5fBS4B/hX4dx/jNqZNPWYuJmNc1wIbVPVjd387cEpEhgHFqjqhlWuuoHlOolavB2eRFXdyt/uAOY0ljvNR1U/cEsANOKUJb/OAqcB659ZE0TyLZmM106vuz9taifmLwF9V9bj7rFIRicOZrXQbgIhsx1kpTEVkGzDSl7iNOR9LEKanyQG813+YCrzrHv/MjJ0iEg0MUNXGevhzXY/boD0YKHFX6muPZTiL1MwBErxDAJ5U1QdaueZV4DciMgWIblzXwjtmN6m05rTX6wav/Qbs/7XxE6tiMj1NCTAeQESm4nxr34pTtZPbyvmXAe+c73p3fvxncBaaqRKRluv5ns9S4CdejeKNVgHXiUiy+8yBIjICQFWr3NiWcna7hXfMbwNfEpGExuvbGZcxHWYJwvQ0TwGT3K6m3wfKcdY8yAH+2av9YbO74l5T+0Mb1+fjLCrzXVXdCfw/nPYIn6lqoar+TyvHd+C0IawUkU+Af+CUUho9B0zk7ATRFLM7pf3PgPdEZCvQZ9f/MF3Ppvs2vZqIbAIucFcW6+y9RgKvqer4TgfW9nP8FrN7v3xgWmM7hjG+shKE6dVUdYq//tDi9EDq7z1QLhD8FXPjQDkgDKdtwph2sRKEMcaYVlkJwhhjTKssQRhjjGmVJQhjjDGtsgRhjDGmVZYgjDHGtMoShDHGmFZZgjDGGNMqSxDGGGNa9f8DczOK6iNgsgQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "\n",
    "recomb_birks = tracks_birks[\"n_electrons\"] / (consts.MeVToElectrons * tracks_birks[\"dE\"])\n",
    "recomb_box = tracks_box[\"n_electrons\"] / (consts.MeVToElectrons * tracks_box[\"dE\"])\n",
    "\n",
    "ax.plot(np.sort(tracks_box[\"dEdx\"]), recomb_box[np.argsort(tracks_box[\"dEdx\"])], label=\"Box model\", lw=2)\n",
    "ax.plot(np.sort(tracks_birks[\"dEdx\"]), recomb_birks[np.argsort(tracks_birks[\"dEdx\"])], label=\"Birks model\", lw=2)\n",
    "ax.set_xlabel(\"$dE/dx$ [MeV/cm]\")\n",
    "ax.set_ylabel(\"Recombination factor\")\n",
    "ax.set_xlim(0,10)\n",
    "_ = ax.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracks_quenched = np.copy(tracks)\n",
    "quenching.quench[blockspergrid,threadsperblock](tracks_quenched, consts.box)\n",
    "tracks_drifted = np.copy(tracks_quenched)\n",
    "drifting.drift[blockspergrid,threadsperblock](tracks_drifted)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Electron lifetime\n",
    "The electrons traveling towards the anode are absorbed by impurities present in the liquid argon. The decrease follows an exponential curve $N_e^{\\mathrm{anode}}/N_e=\\exp(-t_{\\mathrm{drift}}/e_{\\mathrm{lifetime}})$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAEKCAYAAADq59mMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0PklEQVR4nO3dd3xUdfb/8ddJJyEJLdTQeychgmUV1wqsiiJI1fW7609poqvYlrXh2rsClm2uUsSGvStW1pICofdeQwmEhLSZ8/tjLpplKSmT3JnJeT4e8/DOvXdm3nc35OTez9zPEVXFGGOMqYgwtwMYY4wJPlY8jDHGVJgVD2OMMRVmxcMYY0yFWfEwxhhTYVY8jDHGVFiE2wFqQqNGjbRNmzZuxzDGmKCSkZGxR1WTjrWtVhSPNm3akJ6e7nYMY4wJKiKy6Xjb7LKVMcaYCrPiYYwxpsKseBhjjKkwKx7GGGMqzJXiISL/FJHdIrL0ONtFRJ4RkbUiki0iqWW2/V5E1jiP39dcamOMMUe4debxEjDwBNsHAR2dx7XAcwAi0gC4G+gP9APuFpH61ZrUGGPM/3CleKjqN8C+E+wyBHhZfX4A6olIM+BC4DNV3aeq+4HPOHERqpIt+wq45t8/szuvsLo+whhjqoXXq0yak8lHS3ZUy/sH6phHC2BLmedbnXXHW/8/RORaEUkXkfScnJxKhbjv/eV8vmI3N81bjNdrfU+MMcHjhW/W8372Dqa+vZSDhSV+f/9ALR5VpqovqmqaqqYlJR3zBsmTuu/SHjSMi+K7tXt47ut1fk5ojDHVI2PTfh77dBUAjw3vRUJMpN8/I1CLxzagZZnnyc66462vFk0SYnj8it4APPHZatI3nuhKmzHGuO9AQQmT52bh8SrX/KYt53RpUi2fE6jF413gKudbV6cCB1R1B/AJcIGI1HcGyi9w1lWbszs35rqz2uHxKpPnZpFbUFydH2eMMZWmqtz2Zjbbcg/TOzmRWwd2qbbPcuurunOB/wCdRWSriPxRRMaJyDhnlw+B9cBa4G/ABABV3QfcB/zsPKY566rVlAs706dlPbYfKOTWN7Kxvu/GmEA068fNfLxsJ/HRETw7KpWoiOr7FS+14RdhWlqaVnVixC37Chj8zLfkFZZy7yXd+f3pbfwTzhhj/GD59oNcOvN7iku9PDsqhYt7N6/ye4pIhqqmHWtboF62CjgtG8Ty8OW9ALj/gxUs3XbA5UTGGONTUFzKpLmZFJd6GdWvpV8Kx8lY8aiAwT2bMaZ/K4o9Xq6fm8WholK3IxljDHe9s4z1Ofl0alKXuy7qXiOfacWjgu68qBtdmsazYU8+d719zNlVjDGmxszP2sobGVuJiQxj+uhU6kSF18jnWvGooJjIcKaPTqFOZDhvZW3jjYytbkcyxtRS63MOMXW+74/Yey7uTqcm8TX22VY8KqFD43juHeI7Nbzz7aWs3X3I5UTGmNqmqNTD9XOzKCj2cFGvZow4peXJX+RHVjwqaXjfZC7t05zDJR4mzcmksMTjdiRjTC3y4IcrWbb9IK0axPLg0J6ISI1+vhWPShIR/npZT9o0jGXlzjzu/2CF25GMMbXEp8t28tLCjUSGC9NHpxBfDdOPnIwVjyqoGx3B9NGpRIWH8coPm6pt9kpjjDliW+5hbnkjG4DbBnahV3I9V3JY8aiiHi0SuWOwbwqAW9/MZsu+ApcTGWNCVanHyw1zszhwuIRzujTmj79p61oWKx5+cPXpbTi/WxPyCkuZ/GoWJR6v25GMMSHoqc/XkL5pP00SonlseO8aH+coy4qHH4gIjw7rRfPEGLI25/L4p6vdjmSMCTHfr93DjK/WEibw9MgUGsRFuZrHioef1IuN4ulRKYSHCc9/vY6vV1euAZUxxhwtJ6+IG+ctQhUmn9uRU9s1dDuSFQ9/OqVNA/50XkcAbpq3iN0HrX2tMaZqvF7lptcWkZNXRP+2Dbj+nI5uRwKsePjd+LM7cEaHhuzNL+bGeYvwWPtaY0wVvPjter5ds4f6sZE8PdJ3dSMQWPHws/Aw4ckRfWhUN4qF6/by3Fdr3Y5kjAlSmZv389gnvnayj1/Rm6aJMS4n+pUVj2rQOD6Gx6/oA/ja1/5s7WuNMRV04HAJ18/JorSa28lWlhWPajKgUxLjBrTHqzB5bhb78619rTGmfFSV2512sr2quZ1sZVnxqEY3X9CJ1Fb12HGgkFusfa0xppxm/7iZj5bupG50BM+OSqnWdrKVFXiJQkhkeBhPj0whISaCz1fs4qWFG92OZIwJcCt2HGTa+8sBeHBoT1o3jHM50bFZ8ahmLRvE8sgwX/vaBz9cae1rjTHHlV9UyqQ5vnayI0+pmXaylWXFowYM7NGMK09tTbHHy6Q5mda+1hhzTHe+s5R1TjvZuy+umXaylWXFo4ZM/V1XujZLYOPeAqbOX2LjH8aY//J6+hbeytxGnchwZtRgO9nKsuJRQ460r42NCuedRdt53drXGmMca3blcdc7ywCYNqQ7HWuwnWxlWfGoQe2T6jJtSA8A7n5nGWt25bmcyBjjtsPFHibOyeRwiYehqS0Ynlaz7WQry4pHDRvWN5mhKS04XOL8wBRb+1pjarN73l3G6l2HaJ8Ux33OH5fBwIqHC+67tAftk+JYvesQ97y7zO04xhiXvJ21jXnpW4iOCGPGmFTioiPcjlRuVjxcEBcdwYwxqURHhDEvfQvzs2z8w5jaZl3OIf48fwkA91zSnS5NE1xOVDFWPFzSpWkC91zi+yre1PlLWbv7kMuJjDE1pbDEw8TZmRQUe7i4d3NGnhIc4xxlWfFw0chTWnJJ7+YUFHuYNCeTwhIb/zCmNrjv/eWs3JlHm4axPHBZD1fbyVaWFQ8XiQgPDO1J20ZxrNyZx73vLXc7kjGmmr2fvZ3ZP24mKjyM6aNTiY+JdDtSpVjxcFnd6Aimj/ZNfDb3p828u3i725GMMdVk4558bn/TN85x50Vd6dEi0eVElWfFIwB0b57IXRd1A+CON7PZsCff5UTGGH8rKvUwaa5veqLBPZsy9tTWbkeqEiseAWJM/1b8rlcz8ot9A2k2/mFMaPFNjHqQlg3q8NDlvYJynKMsKx4BQkR4aGhPWjeMZfmOg9z/wQq3Ixlj/OTjpTt4aeFGIsOFGaNTSQjScY6yrHgEkPiYSGaMTiUqPIxXftjE+9k2/mFMsNuyr4Bb3sgG4I5BXemVXM/dQH5ixSPA9GiRyNTfdQXg9jeXsGmvjX8YE6yKS71MmptFXmEp53drwv+d0cbtSH5jxSMAXXVaawZ2b8qholImzsmkqNTGP4wJRo98vJLFW3JpUa8Ojw4L/nGOsqx4BCAR4eFhvWjZoA5Ltx3kwQ9Xuh3JGFNBny/fxd+/20BEmPDs6BTqxUa5HcmvrHgEqMQ6kUwflUpkuPDSwo18vHSH25GMMeW0LfcwN7++GIBbLuxMaqv6LifyP9eKh4gMFJFVIrJWRG4/xvbWIvKFiGSLyFciklxm2yMiskxEVojIMxJK54Jl9G5Zj9sH+cY/bnkjmy37ClxOZIw5mRKPl8lzszhwuITfdk7i/53Zzu1I1cKV4iEi4cAMYBDQDRglIt2O2u0x4GVV7QVMAx50Xns6cAbQC+gBnAIMqKHoNe4PZ7Th/G5NyCssZdKcTIpLvW5HMsacwOOfriZj036aJsTw+BV9CAsLyb9tXTvz6AesVdX1qloMvAoMOWqfbsCXzvKCMtsViAGigGggEthV7YldIiI8Nqw3LerVYfHWAzz8sY1/GBOoFqzazfNfryPcGedoEBda4xxluVU8WgBbyjzf6qwrazEw1Fm+DIgXkYaq+h98xWSH8/hEVf/njjoRuVZE0kUkPScnx+8HUJMSYyOZPjqFiDDhH99t4NNlO92OZIw5ys4Dhdz8mm+c46bzO3FKmwYuJ6pegTxgPgUYICJZ+C5LbQM8ItIB6Aok4ys454jImUe/WFVfVNU0VU1LSkqqydzVIqVVfW4b2AWAKa8vZut+G/8wJlCUOuMc+/KLObNjI8YPaO92pGrnVvHYBpTtfpLsrPuFqm5X1aGqmgJMddbl4jsL+UFVD6nqIeAj4LQaSe2ya85sy7ldGnOwsJTr52ZR4rHxD2MCwdNfrOGnjftoHB/NkyNCd5yjLLeKx89ARxFpKyJRwEjg3bI7iEgjETmS7w7gn87yZnxnJBEiEonvrKRWTAQlIjw2vDfNE2PI2pzLo5+scjuSMbXed2v2MH3BWsIEnh6ZQqO60W5HqhGuFA9VLQUmAZ/g+8X/mqouE5FpInKJs9vZwCoRWQ00Ae531r8BrAOW4BsXWayq79VkfjfVj4vi2dEphIcJL36zni9Xhux3BYwJeLvzCrlxXhaqMPncjpzWvqHbkWqMqKrbGapdWlqapqenux3Dr577ah0Pf7ySerGRfDj5TJrXq+N2JGNqFY9XufIfP7Jw3V5Ob9+QV/7Yn/AQu1wlIhmqmnasbYE8YG5O4Lqz2nF25yRyC0qYbOMfxtS4Z79cw8J1e2lUN4qnRvQJucJxMlY8glRYmPD48N40TYghfdN+nvhstduRjKk1vl+7h6e/WIMIPDUihcYJMW5HqnFWPIJYw7rRPDPKN/7x3FfrWLBqt9uRjAl5uw8WcsOrvnGO68/pyG86NnI7kiuseAS5fm0bcNP5nQC4+bXF7DxQ6HIiY0JXqcfL9XOz2HOomNPbN+SGczu6Hck1VjxCwPgB7TmzYyP25RczeW4WpTb+YUy1eOrzNfy4YR9J8dE8NbL2jXOUZcUjBISFCU+O6EPj+Gh+2riPJz+38Q9j/O3r1TnM+Mp3P8czI1NoHF/7xjnKsuIRIho54x9hAjMW2PiHMf6048Bh/jRvEarwp/M61ar7OY7HikcIObVdQ26+oDMAN81bxPbcwy4nMib4HT1v1cTfdnA7UkCw4hFixg9oz4BOSewvKGHSnEy7/8OYKnrs09X8vHE/TRKieaqWzFtVHlY8QsyR8Y9miTFkbs7lEev/YUylfbly16/9OUal0rCWzFtVHlY8QlCDuKhf+n/87Vvr/2FMZWzLPcxNTn+Omy/oRL+2od2fo6KseISovq0b/Ff/D+t/bkz5FZd6mTQnk9wCXx/ycWeFfn+OirLiEcKuObMt53VtwsHCUibOyaSo1ON2JGOCwiMfryRrcy7NE2N4IoT7kFeFFY8QJuKb/yq5fh2ytx7gwQ9t/MOYk/l02U7+/t0GIsKEZ0enUj+E+5BXhRWPEJcYG8mM0alEhgsvLdzIB9k73I5kTMDasq+AKa/7xjluG9iFvq3ru5wocFnxqAV6t6zH1MFdAbjtzWw27sl3OZExgefIOMfBwlLO69qEa85s63akgGbFo5b4/eltGNyzKYeKSpkwO5PCEhv/MKasBz5cweKtB2hRrw6PD++NiI1znIgVj1pCRHjo8l60bhjL8h0Hufe95W5HMiZgfLRkBy8t3EhkuDBjTCqJsZFuRwp4VjxqkYQY3/hHVEQYc3/azNtZ29yOZIzrNu3N59Y3sgG4Y1BX+rSs526gIGHFo5bp0SKRey7uDsCf5y9h7e5DLicyxj2FJR4mzskkr6iUgd2b8n9ntHE7UtCw4lELjerXkiF9mlNQ7GHC7AwOF9v4h6md/vrBcpZuO0jLBnV4eFgvG+eoACsetZCI8MBlPWmfFMfqXYe4852lbkcypsa9t3g7s37YTFR4GDNH9yWxjo1zVIQVj1oqLjqCmWP6EhMZxhsZW3ktfYvbkYypMetzDnH7m75xjjsv6krP5ESXEwUfKx61WOem8dw3pAcAd72zlJU7D7qcyJjqV1jiYcLsTPKLPVzUqxljT23tdqSgVKXiISJPi8i/nOUL/BPJ1KThaS0Z3jeZwhIvE2Zncqio1O1IxlSre99bxsqdebRtFMeDQ3vaOEclVfXMwwtscJbPqeJ7GZdMG9KDzk3iWZ+Tz5/fWoKquh3JmGrxdtY25v60haiIMKaPTiE+xsY5KquqxaMASBSRSKCVH/IYF9SJCmfGmFRio8J5d/F25vy02e1Ixvjd2t2H+PP8JQDcc3F3uje3cY6qqGrxuBtYB8wA5lQ9jnFLh8Z1eXBoTwDufW85S7cdcDmRMf5zuNjDxNmZFBR7GNKnOaP6tXQ7UtArV/EQkWOOKKlqqarOVNVrVfV9/0YzNW1InxaM7t+K4lIvE+dkcrCwxO1IxvjFXe8sZdWuPNolxfHAZTbO4Q/lPfP4SET2iMh3IjJTRMaLyBkiklCt6UyNu+uibnRrlsCmvQXc9ka2jX+YoPda+hZez9hKTGQYM8ekEhcd4XakkFCu4qGq3YAWwPXAD0B74E5gpYhsONFrTXCJiQxn5phU4qMj+GjpTl5auNHtSMZU2vLtB7nzbd9NsNOG9KBLU/t711/KPeahqkWqmgXMB34EdgKFwOJqymZc0qZRHI8M6wX4pqletCXX3UDGVEJeYYnTftnLFWnJXJFm4xz+VN4xj84icpOIfAksBE4DZgNdVPXSasxnXDKoZzOuPr0NJR5l4uxMcguK3Y5kTLmpKre9mc2GPfl0aRrPNOdmWOM/5T3zWAGMAZ4H0lT1JlX9TFXtN0oI+/PgrvRuWY9tuYeZ8vpiG/8wQeNf32/kwyU7qRsdwXNj+xITGe52pJBT3uIxHt8Zx0Rgi4isEJHXROROEbm02tIZV0VFhDF9VAoJMRF8vmI3f/t2vduRjDmpjE37eeDDFQA8OqwXbRvFuZwoNJW3eGQDk1V1gKo2Bs4H/gUUA5dXVzjjvpYNYnn8ij4APPzxKtI37nM3kDEnsC+/mElzMin1Kn84oy2DejZzO1LIKm/xuArIEJFXReRqoFRVP1LVh1X1yuqLZwLB+d2acO1Z7fB4lUlzsthzqMjtSMb8D69XuXHeInYcKCS1VT1uH9TF7Ughrbxf1R2vqqnAPUB94CUR+Y+IPCAiZ4mIXVAMcbdc2Jm01vXZebCQG19dhMdr4x8msExfsJZvVufQIC6K6U67ZVN9KvS/rqquVNUnVXUgvokQvwOG4/vqrglhkeFhTB+dSsO4KL5bu4env1jjdiRjfvHdmj08+flqROCpEX1oXq+O25FCXqVLs6oeVtUPVfV6VU2r6OtFZKCIrBKRtSJy+zG2txaRL0QkW0S+EpHkMttaicinzsD9chFpU9njMOXXNDGGZ0alIALPfrmGr1btdjuSMew8UMgNr2ahCpPP6chZnZLcjlQrlPc+jxtE5G/O8p1V/VDnMtcMYBDQDRglIt2O2u0x4GVV7QVMAx4ss+1l4FFV7Qr0A+y3WA05o0MjbjqvE6rwp3mL2JZ72O1IphYr8XiZNCeTvfnF/KZDIyaf29HtSLVGec882gNH+pTG++Fz+wFrVXW9c6/Iq8CQo/bpBnzpLC84st0pMhGq+hmAqh5S1QI/ZDLlNPG3HTi7cxL7C0qYODuT4lKv25FMLfXoJ6tI37SfpgkxPD2yD+FhNuFhTSlv8VCgjoj0AJr74XNb8GsxAtjqrCtrMTDUWb4MiBeRhkAnIFdE3hKRLBF59FgD9iJyrYiki0h6Tk6OHyKbI8LChCev6EPzxBgWbcn95Tv1xtSkj5fu5MVv1hMeJkwfnULDutFuR6pVyls8HgcEuBK4o/ri/JcpwAARyQIGANsADxABnOlsPwVoB1x99ItV9UVVTVPVtKQkuwbqb/XjopgxJpXIcOGlhRt5P3u725FMLbJpbz63vO6bVu+OQV1Ia9PA5US1z0mLh4i8jm9aklOAPsDf/PC524Cys5QlO+t+oarbVXWoqqYAU511ufjOUhY5l7xKgbeBVD9kMhWU0qo+f/mdb6jqtjeyWZdzyOVEpjYoLPEwflYmeUWlXNi9CX/8TVu3I9VKJy0eqjocSAcuUNULgc/88Lk/Ax1FpK2IRAEjgXfL7iAijUTkSL47gH+WeW09ETlyOnEOsNwPmUwlXHVaay7q1Yz8Yg/jZ2VQUFzqdiQT4u59bxnLdxykdcNYHh3e2xo7uaS8l606AS1EpBm+y0RV4pwxTAI+wTfp4muqukxEponIJc5uZwOrRGQ10AS433mtB98lqy9EZAm+y2n+OBsylSAiPHR5L9onxbF61yH+Mn+pTaBoqs2bGVuZ+9MWoiJ8jZ0SYiLdjlRrSXn+oYtIZ+D/4Rs4/4eqrqzuYP6Ulpam6enpbscIaat35TFk+vccLvHwwGU9Gd2/lduRTIhZtTOPITO+o7DEy0NDezKyn/2MVTcRyTjefXzlnZ5klapOUdVbgq1wmJrRqUk8Dw7tCcA97y5jydYDLicyoeRQUSnjZ2dQWOLl8tRkRpxijZ3cVt6bBC9xpmB/VUSOvh/DGAAuTWnBmP6tKPZ4mTAngwMFJW5HMiHgSGOn9Tn5dG4Sz18v7WHjHAGgvGMeF6nqFao6EhhYnYFMcLvzom70bJHIln2Hufn1RXhtAkVTRS//ZxMfZO8gLiqcmWNTqRNl87AGgvIWjzrOXFOtAeusYo4rJjLcGcj0NZB60RpImSpYtCWXv37g+zLlw8N60T6prsuJzBHlLR734OsiOAnfPFPGHFfLBrE84TSQevSTVfywfq+7gUxQ2p9fzMTZmZR4lKtPb8NFvfwxuYXxl/IWjxH4buRLxoqHKYfzujVh/Nnt8XiV6+dmsTuv0O1IJoh4vcpNr/km3uzTsh5/HtzV7UjmKOUtHqKqo1V1lKqOrtZEJmTcfH4n+rdtQE5eEZPnZlHqsQkUTfk89/U6FqzKoV5sJDPGWGOnQFTuWXVFZLiIDBaRwdWayISMiPAwnh2dQlJ8ND+s38cTn612O5IJAgvX7eHxT1cB8OSIPrSwxk4BqbzF42ugDtDIeRhTLo3jY3h2VAphAjO/WscXK3a5HckEsF0HC5k8NwuvwqTfduC3nRu7HckcR3mLxyqgP/BbfNOGGFNup7ZryC0XdgF8DaS27LP2K+Z/lXq8XD83iz2Hijm9fUP+dH4ntyOZEyhv8bgG2A/cDfzXdy9F5GkR+ZezfIF/45lQcd1Z7Tiva2MOFpYyYXYmhSUetyOZAPPop6v4acM+GsdH8/TIFGvsFODKWzx2ATGAF98khWV5gQ3O8jl+ymVCTFiY8PjwPiTXr8OSbQe4732bCNn86uOlO3nh6yONnVJJirfGToGuvMVjNr6eHrfya2vYIwqARBGJBGymMnNcibGRPDemL1HhYcz+cTPzs7a6HckEgA17/ruxU7+21tgpGJS3eOSr6lpVnayq84/adjewDpgBzPFrOhNyeiYncvclvgZSf35rKat35bmcyLipoLiUca9kkFdUyqAeTa2xUxApb/H4SET2iMh3IjJTRMaLyBkikqCqpao6U1WvVdX3qzWtCQmj+7XispQWHC7xMG5WBoeKrIFUbaSqTJ2/lFW78mjXKI5HhvWyCQ+DSHmnZO8GtACuB34A2gN3AitFZMOJXmvM0USE+y/rQacmdVmfk89tb2ZbA6layHfpcht1IsN5bmxf4q2xU1Ap922bqlqkqlnAfOBHYCdQCCyupmwmhMVGRTBzTF/iosL5IHsH//p+o9uRTA1atCWXae/5vjTx4NCedG4a73IiU1Hl7efRWURuEpEvgYXAafgG0buo6qXVmM+EsA6N6/Lo8N4APPDhCtI37nM5kakJ+5wJD4s9Xq46rTWXprRwO5KphPKeeawAxuD7xlWaqt6kqp+panH1RTO1weCezbjmN20p9SoTZmeSk1fkdiRTjTxe5YZXs36Z8HDq72zCw2BV3uIxHt8Zx0Rgi4iscDoL3ikil1ZbOlMr3DaoC6e0qc/uvCKun5tpEyiGsGe+WMO3a/ZQPzaSmWNSiY6wxk7BqrwD5i+o6vWqOkBVGwPnA/8CioHLqzOgCX2R4WHMcG4M+2H9Ph51JsUzoWXBqt088+UaROCZUSk0twkPg1ql5jlW1a2q+pGqPqyqV/o7lKl9GifEMH2Ub0qKF75ez8dLd7odyfjRln0F/GneIlR9U/Wf2THJ7UimimySfBMw+rdryO0DfRMoTnl9MetzDrmcyPhDYYmHCbMzyS0o4ZwujZlwdge3Ixk/sOJhAso1Z7ZlUI+mHCoqZfysTAqK7QbCYDft/eUs2XaA5Pp1ePKKPoTZhIchwYqHCSgiwiPDetEuKY5Vu/KYOn+p3UAYxN7I2MqcHzcTFRHG82P7khhrNwKGCiseJuDEx0Ty/Ni+1IkMZ37WNmb9sMntSKYSlm8/yNT5SwC4b0h3erRIdDmR8ScrHiYgdWoSz0OX9wR8lz2yNu93OZGpiAOHSxg/O4OiUi9XpCUz4hSbcDvUWPEwAWtInxZcfXobSjy+Gwj3HrIbCIOBqjLl9cVs2ltAt2YJTBvSw+1IphpY8TAB7c+Du5Laqh47DhRyw6uL8Hht/CPQvfDNej5bvouEmAieH9uXmEi7ETAUWfEwAS0qIowZY1JpGBfFd2v38ORnq92OZE7gP+v28sjHKwF44oo+tGoY63IiU12seJiA1yyxDs+OSiFMYPqCtXyxYpfbkcwx7DpYyPVzM/EqTPxte87rdnTHahNKrHiYoHB6h0ZMubAzAH+at4jNewtcTmTKKvF4mTg7kz2HijmjQ0NuOr+z25FMNbPiYYLGuLPac17XJhwsLGXcrAwKSzxuRzKOhz5aSfqm/TRNiOHpkb5pZkxos+JhgkZYmPD4Fb1p3TCW5TsOcufbdgNhIPggewf/+G4DEWHCjDGpNKob7XYkUwOseJigklgnkufG9CUmMozXM7Yy7+ctbkeq1dbuPsStb/iaif7ld13p27q+y4lMTbHiYYJOt+YJ3H+p7wbCu95dxpKtB1xOVDvlF5UyflYG+cUeLu7dnN+f3sbtSKYGWfEwQenyvsmM7t+K4lIv42ZlkFtgTS1rkqpyx1tLWLP7EB0a1+WhoT0RsXGO2sSKhwlad1/cjV7JiWzLPcyN8xbhtRsIa8zL/9nEu4u3ExcVzvNj+xIXHeF2JFPDrHiYoBUdEc7MManUi43kq1U5PPvlWrcj1QoZm/bz1w+WA/DwsF50aFzX5UTGDa4VDxEZKCKrRGStiNx+jO2tReQLEckWka9EJPmo7QkislVEptdcahNokuvH8vTIFETgqS9W89Wq3W5HCmk5eUVMmJ1BiUf5vzPacFGv5m5HMi5xpXiISDgwAxgEdANGiUi3o3Z7DHhZVXsB04AHj9p+H/BNdWc1gW9ApyRuPLcTqnDjvEVs3W83EFaHUo+XSXMy2XWwiFPa1OfPg7u6Hcm4yK0zj37AWlVdr6rFwKvAkKP26QZ86SwvKLtdRPoCTYBPayCrCQLXn9OBszsnkVtQwoTZmXYDYTV4+OOV/LhhH0nx0cwYnUpkuF31rs3c+n+/BVD2C/pbnXVlLQaGOsuXAfEi0lBEwoDHgSkn+gARuVZE0kUkPScnx0+xTaAKCxOeGtGH5Pp1yN56gGnvL3c7Ukh5P3s7f/vWdyPgzDGpNE6IcTuScVkg/+kwBRggIlnAAGAb4AEmAB+q6tYTvVhVX1TVNFVNS0pKqv60xnX1YqN4bkxfoiLCmPPjZt7IOOGPiCmnNbvyuPWNbACm/q4rp7Rp4HIiEwjcKh7bgJZlnic7636hqttVdaiqpgBTnXW5wGnAJBHZiG9c5CoReagmQpvA1zM5kWmXdAdg6vwlLNtuNxBWRV5hCde9kkFBsYchfZpztd0IaBxuFY+fgY4i0lZEooCRwLtldxCRRs4lKoA7gH8CqOoYVW2lqm3wnZ28rKr/820tU3uN7NeKK9KSKSr1Mn5WJgcKStyOFJSOdARcvyefLk3jedBuBDRluFI8VLUUmAR8AqwAXlPVZSIyTUQucXY7G1glIqvxDY7f70ZWE5ymDelBjxYJbN5XwI3zsuwGwkp4/uv1fLJsF/FOR8DYKLsR0PxKasOspGlpaZqenu52DFPDtuwr4OLp35FbUMIN53bkT+d3cjtS0Ph+7R6u/MePeBX+flWaNXaqpUQkQ1XTjrUtkAfMjamSlg1ieca5gfDpL9bw5UrrQFge23IPc/3cLLzq+wq0FQ5zLFY8TEg7q1MSUy7wdbW78dVFbNqb73KiwFZY4mHCrAz25RdzVqckbjzPztbMsVnxMCFv/ID2nN/N14HwulcyOFxsNxAez73vLWfx1gMk16/D0yP6WEdAc1xWPEzIO9KBsF2jOFbuzOP2t7KtA+ExvPbzFub+tJmoiDCeH9uX+nFRbkcyAcyKh6kVEmIief7KvsRGhfPOou28tHCj25ECypKtB/jLO0sB+OulPejRItHlRCbQWfEwtUanJvE8MqwXAPd/sIKfN+5zOVFg2J9fzLhZGRSXehndvxVXpLU8+YtMrWfFw9QqF/VqzjW/aUupV5kwO5PdBwvdjuQqj1eZ/GoW23IP07tlPe6++OjJrY05Nisepta5fVAX+rdt4PSmyKS41Ot2JNc8+dlqvl2zh4ZxUTw3JpXoiHC3I5kgYcXD1DoR4WFMH51K04QY0jft54EPV7gdyRWfLtvJ9AVrCRN4dlQKzevVcTuSCSJWPEytlBQfzcyxqUSGCy8t3Mj8rNo1A++GPfnc/NpiAG4d2IXTOzRyOZEJNlY8TK2V2qo+d1/sm4H3jreWsHz7QZcT1YyC4lLGvZJBXlEpA7s35bqz2rkdyQQhKx6mVhvTvxXD+iZTWOJl3KyMkJ+BV1W57c0lrNqVR/ukOB4d3stmyjWVYsXD1Goi4tzXUDtm4P3X9xt5b/F24qLCeeHKvsTHRLodyQQpKx6m1ouJDOe5MX2pFxvJglU5PPPlGrcjVYufNuz75csBjwzrTYfG8S4nMsHMiocxhP4MvLsPFjJxTialXuXas9rxu17N3I5kgpwVD2McZ3VK4ubzO6EaWjPwFpd6mTA7k5y8Ik5t14BbL+zsdiQTAqx4GFPGhLM7cF7X0JqB94EPV5C+aT9NE2J4dlQqEeH2z95Unf0UGVNGWJjwxIjetHVm4L0jyGfgfWfRNl5auJHIcGHm2FSS4qPdjmRChBUPY46SEBPJ82P7UicynLcXbeffQToD74odB7ntzWwA7rq4O6mt6rucyIQSKx7GHEPnpr/OwPvXIJyB90BBCeNmZVBY4uXy1GTG9m/ldiQTYqx4GHMcF/duzh+DcAZer1e5cV4Wm/YW0K1ZAvdf1sNuBDR+Z8XDmBO4fVAX+jkz8E6ck0mJJ/Bn4H3q89UsWJVDvdhIXriyLzGRNlOu8T8rHsacQGR4GDNGp9IkIZqfN+7n/g8CewbeT5ft5Jkvf50pt2WDWLcjmRBlxcOYk0iKj2bmmL6/zMD7zqJtbkc6prW7D3GTM1PuLRd24cyOSS4nMqHMiocx5dC3dX3ucmbgve3NbFbsCKwZePMKS7julXQOFZUyuGdTxg2wmXJN9bLiYUw5je3fistTfTPwXvdKBrkFxW5HAnwD5FNeX8y6nHw6NanLo8N62wC5qXZWPIwpJxHh/st+nYH3hlcX4QmAGXif+3odnyzbRXxMBC9cmUZcdITbkUwtYMXDmAqIiQzn+bF9qR8byderc3jis1Wu5vlq1W4e+9SX4akRfWjbKM7VPKb2sOJhTAUl149lxuhUwgRmLFjHx0t3uJJj817f2Y8q3HheR87t2sSVHKZ2suJhTCWc3qERdwzqCsDNry1mza68Gv38guJSrn0lnQOHSziva2Mmn9OxRj/fGCsexlTSNWe25eLezckv9nDdKxkcLKyZFraqyu1vLmHlzjzaNorjiRF9CAuzAXJTs6x4GFNJIsLDl/ekS9N41u/J56Z5i2ukhe0/vtvAu04r2Rev7EuCtZI1LrDiYUwVxEZF8MKVfUmIieDzFbt49su11fp5C9ft4cGPVgLw2PDedGxirWSNO6x4GFNFrRvG8cwoXwvbp75YXW0tbLflHmbSnCw8XmX82e0Z1NNayRr3WPEwxg/O7tyYKRd0RhVueHURG/b4t4VtYYmH8bMy2JdfzJkdGzHlAmsla9xlxcMYP5lwdnsu7N6EvMJSrn05nfyiUr+8r6py59tLyd56gOT6dXhmZArhNkBuXGbFwxg/EREev6IPHRrXZc3uQ9zyxmK/tLCd9eNmXs/YSkxkGC9c2Zf6cVF+SGtM1VjxMMaP6kb7BtDjoyP4cMlOnv96fZXeL2PTPqa9twyAh4b2onvzRH/ENKbKrHgY42ftk+ryxIg+ADz6yUq+XZNTqffZfbCQcbMyKfEofzijLZemtPBjSmOqxrXiISIDRWSViKwVkduPsb21iHwhItki8pWIJDvr+4jIf0RkmbNtRM2nN+bEzu/WhMnndsSrcP3cLLbsK6jQ64tLvYyfnUlOXhH92zbgjsFdqimpMZXjSvEQkXBgBjAI6AaMEpFuR+32GPCyqvYCpgEPOusLgKtUtTswEHhKROrVSHBjKuDGcztyTpfG5BaUcO0rGRwu9pT7tfe9v5yMTftplhjDjDGpRIbbRQITWNz6iewHrFXV9apaDLwKDDlqn27Al87ygiPbVXW1qq5xlrcDuwFrmWYCTliY8OSIPrRpGMuKHQe5463scg2gv5a+hVd+2ERUeBjPje1Lo7rRNZDWmIpxq3i0ALaUeb7VWVfWYmCos3wZEC8iDcvuICL9gChg3dEfICLXiki6iKTn5FTumrMxVZVYJ5IXr0ojNiqctxdt51/fbzzh/tlbc/nL20sBuO/S7vRpWa/6QxpTCYF8LjwFGCAiWcAAYBvwy3m/iDQDXgH+T1W9R79YVV9U1TRVTUtKshMT455OTeJ5bHhvAO7/cAX/Wbf3mPvtOVTEuFcyKC71Mrp/K0ac0qomYxpTIW4Vj21AyzLPk511v1DV7ao6VFVTgKnOulwAEUkAPgCmquoPNZLYmCoY3LMZ4wa0x+NVJs3JZHvu4f/aXurx+tYfKCS1VT3uvvjoIUBjAotbxeNnoKOItBWRKGAk8G7ZHUSkkYgcyXcH8E9nfRQwH99g+hs1mNmYKrnlws6c2bERe/OLGTcrg8KSXwfQH/poJT+s30dSfDTPje1LdES4i0mNOTlXioeqlgKTgE+AFcBrqrpMRKaJyCXObmcDq0RkNdAEuN9ZfwVwFnC1iCxyHn1q9ACMqYTwMOGZkSkk169D9tYD3PXOUlSVdxZt4+/fbSAiTJg5JpUmCTFuRzXmpMQf0ycEurS0NE1PT3c7hjEALNt+gMufW0hhiZeoiDBKPF5UYdqQ7lx1Whu34xnzCxHJUNW0Y20L5AFzY0JS9+aJPDa8N3WjIygu9X3X4w9ntOXKU1u7nMyY8otwO4AxtdFFvZpzYfemeLyKCDbGYYKOFQ9jXBIZHkak1QwTpOyylTHGmAqz4mGMMabCrHgYY4ypMCsexhhjKsyKhzHGmAqz4mGMMabCrHgYY4ypsFoxPYmI5ACb/PBWjYA9fngft9lxBBY7jsBix/Gr1qp6zJ4WtaJ4+IuIpB9vnpdgYscRWOw4AosdR/nYZStjjDEVZsXDGGNMhVnxqJgX3Q7gJ3YcgcWOI7DYcZSDjXkYY4ypMDvzMMYYU2FWPIwxxlSYFY+TEJH7RCTb6ZX+qYg0d9aLiDwjImud7aluZz0REXlURFY6WeeLSL0y2+5wjmOViFzoYsyTEpHhIrJMRLwiknbUtqA5DgARGehkXSsit7udp7xE5J8isltElpZZ10BEPhORNc5/67uZsTxEpKWILBCR5c7P1A3O+qA6FhGJEZGfRGSxcxz3OuvbisiPzs/XPBGJ8usHq6o9TvAAEsosTwaed5YHAx8BApwK/Oh21pMcxwVAhLP8MPCws9wNWAxEA22BdUC423lPcBxdgc7AV0BamfXBdhzhTsZ2QJSTvZvbucqZ/SwgFVhaZt0jwO3O8u1Hfr4C+QE0A1Kd5XhgtfNzFFTH4vwOqussRwI/Or+TXgNGOuufB8b783PtzOMkVPVgmadxwJFvGAwBXlafH4B6ItKsxgOWk6p+qqqlztMfgGRneQjwqqoWqeoGYC3Qz42M5aGqK1R11TE2BdVx4Mu2VlXXq2ox8Cq+Ywh4qvoNsO+o1UOAfzvL/wYurclMlaGqO1Q101nOA1YALQiyY3F+Bx1ynkY6DwXOAd5w1vv9OKx4lIOI3C8iW4AxwF3O6hbAljK7bXXWBYM/4DtrguA+jrKC7TiCLe/JNFHVHc7yTqCJm2EqSkTaACn4/moPumMRkXARWQTsBj7Dd1abW+YPRr//fFnxAETkcxFZeozHEABVnaqqLYHZwCR30x7fyY7D2WcqUIrvWAJSeY7DBC71XScJmnsARKQu8CZw41FXGoLmWFTVo6p98F1R6Ad0qe7PjKjuDwgGqnpeOXedDXwI3A1sA1qW2ZbsrHPNyY5DRK4GLgLOdf5RQBAex3EE3HGcRLDlPZldItJMVXc4l293ux2oPEQkEl/hmK2qbzmrg/JYAFQ1V0QWAKfhu5Qe4Zx9+P3ny848TkJEOpZ5OgRY6Sy/C1zlfOvqVOBAmVPdgCMiA4FbgUtUtaDMpneBkSISLSJtgY7AT25krKJgO46fgY7ON2KigJH4jiFYvQv83ln+PfCOi1nKRUQE+AewQlWfKLMpqI5FRJKOfHtSROoA5+Mbv1kADHN28/9xuP1NgUB/4PurZCmQDbwHtNBfv+EwA9+1xSWU+eZPID7wDSBvARY5j+fLbJvqHMcqYJDbWU9yHJfhu35bBOwCPgnG43DyDsb3DZ91wFS381Qg91xgB1Di/H/xR6Ah8AWwBvgcaOB2znIcx2/wXZLKLvPvYnCwHQvQC8hyjmMpcJezvh2+P6DWAq8D0f78XJuexBhjTIXZZStjjDEVZsXDGGNMhVnxMMYYU2FWPIwxxlSYFQ9jjDEVZsXDGGNMhVnxMKaaiMg9IjLFWV5YyfeoJyITKvG6NiJy2JnvqEpEpI7TkqBYRBpV9f1MaLDiYYwfODMNHPffk6qeXsm3rgdUuHg41qlvvqMqUdXDzvtsr+p7mdBhxcOEJBH50vlreZGIFIrIFUdtv8ppjLVYRF5x1t1UZhLGG4/a/3+2OX/drxKRl/Hd2dtSRKaKyGoR+Q5f35Ejrz9U5jUrRORvTuOeT50pJRCRt0Ukw1l/rfPSh4D2znE86uw31mn+s0hEXhCR8HL+b/Jfx+xkWSkiLzmZZ4vIeSLyvfgaIQXylPbGbW7fWm8Pe1TnAxiPrylOeJl13fFNC9LIed4A6Itvmpk4oC6wDEhxth9zG9AG8AKnHrVfLJCAb1qIKc62Q85/2+Cb1biP8/w1YOyRHM5/6+ArRg2d/cs2XeqKb5qcSOf5TOCqYxz30a871jEfydIT3x+SGcA/8U29MwR4+6j33Hjk9fawh515mJAlIlcBg4Axquops+kc4HVV3QOgqvvwzXM0X1Xz1ddY5y3gTGf/E23bpL5mYDjr5qtqgfqm9j7eRIcbVHWRs5yB75c4wGQRWYyvWVdLfJM7Hu1cfEXqZ2c841x8cxidzLGO+UiWJarqxVcUv1BVxVcE2xzznYzBpmQ3IUpEhuNr3jVEVUuq8aPyK/GaojLLHqCOiJwNnAecpqoFIvIVEHOM1wrwb1W9oxKfe7Is3jLPvdjvB3MCduZhQo6IXIRvkHmoqhYeY5cvgeEi0tDZvwHwLXCpiMSKSBy+2Xu/dfY/0bayvnH2qyMi8cDFFYidCOx3CkcXfD2oAfLw9dc+4gtgmIg0PpJdRFqX4/2PdczGVJr9ZWFC0b/x9dj+3teygWdV9R9HNqrqMhG5H/haRDxAlqpeLSIv8WsPkL+rapazf+axtomvdSll3jdTROYBi/E1EPq5Apk/BsaJyAp8U8r/4LznXmcAeynwkareIiJ/AT51vt1VAkwENp3ozY91zMA9FchnzH+xKdmNCUFOYXtfVXv48T034utbs8df72mCl122MiY0eYBEf94kCETiGwsxxs48jDHGVJydeRhjjKkwKx7GGGMqzIqHMcaYCrPiYYwxpsKseBhjjKkwKx7GGGMqzIqHMcaYCrPiYYwxpsL+P6tkcxllTf/HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "drift_time = []\n",
    "for track in tracks_drifted: \n",
    "    drift_time.append(abs(track[\"z\"] - consts.tpc_borders[track['pixel_plane']][2][0])/consts.vdrift)\n",
    "    \n",
    "drift_time = np.array(drift_time)\n",
    "ax.plot(np.sort(tracks_quenched[\"z\"]), np.exp(-drift_time / consts.lifetime)[np.argsort(tracks_quenched['z'])], lw=2)\n",
    "ax.set_ylabel(\"$N_{e}^{\\mathrm{anode}}/N_{e}$\")\n",
    "_ = ax.set_xlabel(\"$z$ coordinate [cm]\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Diffusion\n",
    "The diffusion coefficients (longitudinal and transverse) are proportional to $\\sqrt{2t_{\\mathrm{drift}}}$ (which in turn is given by $(z_{\\mathrm{start}} - z_{\\mathrm{anode}})/v_{\\mathrm{drift}}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEMCAYAAADeYiHoAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABTf0lEQVR4nO3dd3hUVfrA8e+bTkgILbRQEnqvkY6ACNKkCYqKgBQbWNbVXcuuP+vay66CCtJUpIhKU0TpIDVA6C2hhk4gCSGkzvn9cS4xhpA6k5kk5/M882Tm3jv3vheSeefec857RCmFYRiGYWTm5uwADMMwDNdkEoRhGIaRJZMgDMMwjCyZBGEYhmFkySQIwzAMI0smQRiGYRhZcmiCEJHeInJIRCJE5IUs1nuLyDxr/RYRCc6wrrmIbBKRfSKyR0R8HBmrYRiG8VcOSxAi4g5MAvoAjYH7RaRxps3GAleUUnWBj4F3rfd6AN8CjymlmgDdgBRHxWoYhmHczMOB+24LRCiljgKIyFxgILA/wzYDgVet5wuAz0REgF7AbqXULgClVHROB6tYsaIKDg62W/CGYRglwfbt2y8ppQKzWufIBBEEnMrwOgpod6ttlFKpIhILVADqA0pElgOBwFyl1HvZHSw4OJiwsDB7xW4YhlEiiMiJW61zZIIoCA+gM3AbkACsFJHtSqmVGTcSkUeARwBq1qxZ6EEahmEUZ45spD4N1Mjwurq1LMttrHaHACAafbWxTil1SSmVAPwCtM58AKXUFKVUqFIqNDAwyyskwzAMI58cmSC2AfVEJEREvIDhwOJM2ywGRlnPhwKrlK4euBxoJiK+VuLoyl/bLgzDMAwHc9gtJqtNYSL6w94dmK6U2icirwNhSqnFwDTgGxGJAC6jkwhKqSsi8hE6ySjgF6XUz46K1TAMw7iZFJdy36Ghoco0UhuGYeSN1b4bmtU6M5LaMAzDyJKr9mIyjKLvWjTsmAnJCeDhDcGdoVZHZ0dlGLlmEoRhOMKBJbD0b3Dt4l+Xtx4Fvd4EnzLOicsw8sAkCMOwp2vRsOx52PuDfl2rE9TuBvEXYMcs/YhYCQP+B3V7ODVUw8iJSRCGYS/7F8PPz+qrBk9f6Pk6hI4FN6upL3QMLHoCzuyEb4eYqwnD5ZlGasMoqGvR8P3DMP8hnRxqdYbHN0Lb8X8mB4DKjWHsCujxCrh76auJyR30FYVhuCCTIAyjIPYvgkltYd+P4Fka+n4Ao5ZA+ZCst3f3gC5/h0fWQrVWEBelryYWPwWJcYUbu2HkwCQIw8iPa5esq4aRkHAJgrvA43/cfNVwK+ZqwigCTIIwjLw6sBQmtfvrVcPIxbe+ariVG1cTj667+Woi6apjYjeMPDAJwjByKzEWfnoc5j3451XDE1m0NeRVpUY3X0183glObLRf7IaRDyZBGEZuHFuvP7R3fQcePtD7XX3VUC7YPvvPeDVRpTnEnIAZfeH3VyA1yT7HMIw8MgnCMLKTkgi/vgSz+kPsKajaEh5dD+0fK9hVw61UagTjVkKX50AE/vgvTOkO5/ba/1iGkQOTIAzjVs7ugildYfMkEHfo+gKMWwGB9R17XA8v6PFvGLMcyteGC/tgSjfY8DHY0hx7bMPIwCQIw8gsLRXWvQ9T74CLB6FCXRj7O3R/Edw9Cy+OGm3hsQ16gJ0tBVa8CjP7weVjhReDUaKZBGEYGUVHwozesOpNsKVC20f1LaXqbZwTj1dp6P8xPLgA/KrAyU3wRWfYPguKSal+w3WZBGEYoD9sw6brD9+obeBfDR76Cfq+B16+zo4O6vWEJzZB40GQHA9LnoI5w3WNJ8NwEJMgDONaNMx9UFdfTUmAZvfq7qt17nB2ZH/lWx6GzYQhU8E7AA7/Cp93hCMrnB2ZUUyZBGGUbEfX6A/ZQz/rD92h0+GeqVCqnLMjy5oINLcSWHAXXftp9j3w64u6x5Vh2JFJEEbJlJoMv/0bvh4E8eegZgd4fAM0vcfZkeVOQHUYuQh6/B+4ecDmyfBVD7hw0NmRGcWISRBGyXMpAqb1hI3/A3GD7i/DqKVQtqazI8sbN3fo8iyM/U13hz2/V3fL3faVacA27MIkCKPkUAp2fA1fdoGz4TohPLwMuv5Dj2QuqoLa6J5WLUdAaiL8/HeY+4BuWzGMAjAJwigZrl+B70fB4iethuhheoxBzXbOjsw+vP1g0CQYOkO3pRz6RbetRK52dmRGEWYShFH8ndio6yjtXwRe/jB4CtzzFfgEODsy+2s6RLel1Oyo21a+GQQrXtOD/wwjj0yCMIovm02PiJ7ZD+JOQ1AoPLYeWtzn7Mgcq2xNGL1Ut62IG2z4SP8bxEY5OzKjiDEJwiie4i/ouRVWvQnKBp2egTG/5n3OhqLKzV23rYxaCv5V4dRmPQjw0K/OjswoQhyaIESkt4gcEpEIEXkhi/XeIjLPWr9FRIKt5cEicl1Ewq3HF46M0yhmjq3TH4ZHV4NvBXjwB+j5WuHWUXIVwZ10W0vdnrodZs59sPxl3c3XMHLgsAQhIu7AJKAP0Bi4X0QaZ9psLHBFKVUX+Bh4N8O6SKVUS+vxmKPiNIoRWxqseQe+Hgjx56GW9eFY705nR+ZcpSvCA/Oh5+u6Ku2mz3S9qSvHnR2Z4eIceQXRFohQSh1VSiUDc4GBmbYZCMyyni8AeoiIODAmo7i6ek4nhjVv6+6stz+vJ/QpU83ZkbkGNzfo9LS+zRZQA05vhy9uh/2LnR2Z4cIcmSCCgFMZXkdZy7LcRimVCsQCFax1ISKyU0TWikgXB8ZpFHWRq/QtpeProXQlXWTvjn8V7bENjlKjrZ61rkE/SIqF+Q/BL8+bMh1Glly1kfosUFMp1Qp4FvhORMpk3khEHhGRMBEJu3jxYqEHaTiZzQZr3oVvhuiaRCG361tKdbo7OzLX5lsehs+G3u+AmydsnWLdcjrh7MgMF+PIBHEaqJHhdXVrWZbbiIgHEABEK6WSlFLRAEqp7UAkcNM0XkqpKUqpUKVUaGBgoANOwXBZCZfhu3thzX/0624vwUMLwb+yU8MqMkSg/eO6TEfZmnBmpy7TYSrDGhk4MkFsA+qJSIiIeAHDgcw3PBcDo6znQ4FVSiklIoFWIzciUhuoBxx1YKxGUXJmJ3zZFSJ+h1LlYcQP0O2fumunkTdBreGRtVCvl+7lNHuobui32ZwdmeECHJYgrDaFicBy4AAwXym1T0ReF5EB1mbTgAoiEoG+lXSjK+ztwG4RCUc3Xj+mlLrsqFiNImT7LJh2F8SehGqt4dG1ULeHs6Mq2nzLw/3z9MA60A393w3TV2lGiSaqmFR9DA0NVWFhYc4Ow3CUlOvwy3Ow81v9OnSMvofu4e3cuIqbiJXwwzi4fhkCasJ9X0O1Vs6OynAgEdmulArNap2rNlIbxp8uH4NpvXRy8PCBQV/oeZpNcrC/uj30VVm11voqbVov2D7TlA8voUyCMFzboV914+m53VAuBMatgJb3Ozuq4q1sTT1eInQMpCXDkqdh0UR9FWeUKCZBGK7JZoPVb+vSEImx0KAvPLIGqjRzdmQlg4e3vkob9IW+agv/Vl9NxJx0dmRGITIJwnA9SVf1AK617+hqpD3+D+6bDaXKOjuykqfl/fqqrVyIvoqb0g2Ob3B2VEYhMQnCcC2Xj8JXPeHgUj3xzQPf62k13cyvqtNUaQbjV0GdOyAhWpc02TrVtEuUAOavznAdkatgSne4eAAqNoBHVptCe67Ct7xO1h2fBFuq7lG25ClITXJ2ZIYDmQRhOJ9SsGkSfHsPJMZA/T76tkaFOs6OzMjI3QN6vQlDpup2iR1fw6y74ep5Z0dmOIhJEIZzpSTCwsdh+Ut6Yp/bn4fh34HPTaW3DFfR/F7dy6lMEJzaotslTm93dlSGA5gEYThP3BmY0Qd2zQFPXxg2S1dhNe0Nrq9aK92rrGYHuHoGpveB8DnOjsqwM/OXaDjHqa36m+eZHbrf/djfoMkgZ0dl5IVfJT3nRuhYSEuChY/B7/9n6jgVIyZBGIVvzwKY2V/P+hbcBcavMeMbiioPL+j/kR4zIe7wxye6i3LyNWdHZtiBSRBG4VFKVwr9wfrGGTpWT+5TukLO7zVcW+gYXVXXO0B3UZ7RR99CNIo0kyCMwpGSCD8+oiuFihv0fhf6fQjuns6OzLCXOt3/HFR3dhdMvQPOhDs7KqMATIIwHO/aJT24as988CwNw+dA+8f0pDVG8RJYH8athJod4epZfSVxYKmzozLyySQIw7EuHoavesCpzbpb5Njl0KC3s6MyHKl0BRi5EFo8ACkJMG8EbPjEjLwugkyCMBzn6Br46k64chyqttTlGkxjdMng4Q2DJus6WihY8X+6ImxqsrMjM/LAJAjDMXZ8rUdGJ8VCo7vh4V/Av4qzozIKk4iuo3Xv1+BRSleE/W4YJMY5OzIjl0yCMOxLKVj9H1hs1ezp9AwM+xq8Sjs7MsNZGg+Eh3+G0oH6qnJGH4g76+yojFzIMUGISEhulhkGaSmweCKsfVf3VOr/MfR8zYyMNiCoDYz9HcrXgfN79a3HCwecHZWRg9z85f6QxbIF9g7EKOKS4mHO/da0oKV0PaXQMc6OynAl5UN0kqjeFuKiYPpdZm4JF+dxqxUi0hBoAgSIyJAMq8oAPo4OzChC4i/Ad/fCmZ3gWwEemA/Vs5wD3SjpSleAkYvgx/F6QN03g2Hwl9B0SM7vNQrdLRME0ADoD5QF7s6w/Cow3oExGUVJdCR8O0T3VCoXDCN+NGW6jex5+eqG62X/hG1TYcHDetR1x4nOjszI5JYJQim1CFgkIh2UUpsKMSajqIgK01cOCdG6G+uD3+sCboaREzd36Ps+BFTXXWB/exniTkOvt0yblQvJ7grihggReQkIzri9UsrcYC7JjqzQA6BSr0PdnjBsJnj7OTsqoygRgc7P6AGUCx+HzZMh4TIM/MyUYHERuUkQi4D1wAogzbHhGEXC3h/gx0fBlgItR8Ddn5g/aCP/mg/TbRNzR8DuuZAYC8NmgGcpZ0dW4uUmQfgqpf7p8EiMoiFsOix9FlDQYaKegtLUVDIKqs4dMGqxHlx5eJn+ef8c8AlwdmQlWm5u9i0Vkb752bmI9BaRQyISISIvZLHeW0TmWeu3iEhwpvU1RSReRJ7Lz/ENO1v/ESz9G6CgxysmORj2VT1UT2XqXxVO/GHNGXLR2VGVaLlJEE+jk0SiiMSJyFURyXGsvIi4A5OAPkBj4H4RaZxps7HAFaVUXeBj4N1M6z8CluUiRsORlILfX4GVrwEC/T6CLn83ycGwv0qNYMxyKF8bzu2GGb0h5qSzoyqxckwQSil/pZSbUspHKVXGep2bGeXbAhFKqaNKqWRgLjAw0zYDgVnW8wVADxH9qSMig4BjwL5cnovhCLY0WPIU/PFfcPOAe76C28Y6OyqjOCtXSyeJys0gOgKm94aLh5wdVYmUm1IbIiIjROTf1usaItI2F/sOAk5leB1lLctyG6VUKhALVBARP+CfwGs5xPaIiISJSNjFi+ZS1O5Sk2HBGF14z8NHz+PQbKizozJKAr9KMHop1Oygu7/O6ANndzs7qhInN7eYJgMdgAes1/HoW0eO9CrwsVIqPruNlFJTlFKhSqnQwMBAB4dUwqQmwfejYP9C8C6jpwat38vZURklSamyeuBl3Z56rM2s/nB6u7OjKlFykyDaKaUmAIkASqkrgFcu3ncaqJHhdXVrWZbbiIgHEABEA+2A90TkOPAM8JKImGGWhSXlOsx9EA79AqXK6d4ltTo6O6oiKT4plcvXkolLTEGZCXPyzssXhs+GBv1099dZA+HkZmdHVWLkpptritXgrABEJBCw5eJ924B6VuXX08Bw/rwKuWExMArYBAwFVin9V9TlxgYi8ioQr5T6LBfHNAoqOQHm3q/LMvtadXPMJD95opRi+4krfLX+GMv3n0ufSC24gi9jO4cwtE0NSnm5OzfIosTDG+6dpec03/cjfDMEHpgLIbc7O7JiLzcJ4n/AT0AlEXkL/UH+r5zepJRKtb71LwfcgelKqX0i8joQppRaDEwDvhGRCOAyOokYzpIUD3OGw/H1ULqSvnKo1MjZURUZqWk2lu87z9T1Rwk/FQOAh5vgV8qDxJQ0jkcn8O9F+/jw98OMaFeLkR1rUcnf1L3MFXdP3UHCwxt2zYHZw+C+2VDvTmdHVqxJbi57rcquPQABViqlXK6Qe2hoqAoLC3N2GEVXYpz+ozu1GfyqwKglegJ6I0fxSanM33aK6X8cI+rKdQDK+nr+JQmk2RTL951jyro/k4eXuxuDWlVjXJfa1K/s78QzKEJsNvj5b7B9Jrh76SRh2sYKRES2K6WyLL98ywQhImWUUnEiUj6r9Uqpy3aMscBMgiiA6zF65OrpMChTXV85mIqsOTobe52ZG4/z3ZaTXE1MBf68jXRPm+r4et18gX7j9tPU9Uf5bf/59NtP3RoEMr5LbTrWqYCY8SXZU0pXgt36pU4Sw+eYK4kCyG+CWKqU6i8ix7DaH26sApRSqrb9Q80/kyDyKTEWvh4EZ3ZAQE0YvUSX7TZuaf+ZOKauP8qSXWdItek/jduCyzGuS23ubFQZd7fcfcAfu3SN6RuO8f32UySm6Ga9xlXL8MjttenXvCqe7qaq6S0pBcv+AVungLu3LstRt4ezoyqS8pUgihqTIPIh6apu8IvaCmVrweifoWyNnN9XAiml2HLsMp+viWTtYT3mxk2gT7OqjO9Sm5Y1yuZ735evJTN78wlmbTrOpfhkAKoG+DC6YzD3t6tJGR9TCDFLSsHPf4ewaXqczgPzoHY3Z0dV5BQoQYjIYHTvoljrdVmgm1JqoZ3jLBCTIPIo+Rp8OxROboSAGvDwL1C2prOjcjk2m+L3A+f5Ym0kO0/GAFDK053hbWswplMINcr72u1YiSlpLAo/zdT1x4i4oIcA+ft48FD7WozpHEJFP2+7HavYyNgm4VEKHpxvejflUUETRLhSqmWmZTuVUq3sF2LBmQSRBynX9UQ/x9aBfzWdHMqHODsql5KcamNR+Gm+WBtJ5MVrAJTz9WR0xxBGdqhFudK5GQqUPzabYu3hi3yxNpItx3RTn7eHG/fdVoPxXWrbNSkVCzYbLH1aj/j39NUTVwV3dnZURUZBE8RupVTzTMv2KKVcqnO8SRC5lJKoxzlErgK/yvDwMtMgncG1pFTmbD3JtA3HOBubCEBQ2VKM7xLCvbfVyLLh2ZG2n7jC52siWHHgAgDubsKAFtV4vFsd0/MpI5sNFj8J4d+CZ2kYscAM7sylgiaI6UAMf5bXmACUV0qNtmOMBWYSRC6kJutZ4I4sB9+K+sohsIGzo3IJ0fFJzNp4nFmbThB7PQWA+pX9eKxrHe5uUc3pDcaHzl3li7WRLN51hjSrYfzORpV5onsdWtcs59TYXIYtDRZN0OMkvPxg5GKo3sbZUbm8giaI0sC/gRv9yH4H3lRKXbNrlAVkEkQO0lJ1baWDS6FUeV0IrXITZ0fldOdiE/lyXSRztp5M70l0W3A5Hutah+4NKuGWyx5JheXU5QSmrj/KvG2nSErV8bYLKc8T3etye72KpousLQ1+ehT2fA8+ZfUVcuXMswwYGZleTCWdzQaLJ0L4bD1D16ilULV5zu8rxqKuJPDF2kjmb4siOU1/0PZoWInHutXhtuAsh/64lItXk5i58RhfbzqRPgajSbUyPHlHXXo1ruJyia1QpaXA/JG6lphfZT0JUXmX6pXvUvI7DuITpdQzIrKEv46DAEApNcC+YRaMSRC3oBT89i/Y9JluwBu5CGrkplp78XQi+hqTV0fyw44oUm0KEejbtCoT76hLo6q5mebEtcQlpjB7s24zuRSfBEDDKv48eUc9+jQtwYkiJRG+G6Y7YpStCQ//CgGZZxswIP8JorVSaoeIdM1qvVJqrR1jLDCTIG5h3Qew6g1w89T9xEvoYKLIi/FMWhXBIusevpvAgBbVmNC9LvWKQWNvYkoa88NO8fmayPTG9XqV/HiyRz36Naua68F7xUpSPHw9UFcIqFhf324qXdHZUbmc/CaIlUqpHiLyrlLqnw6N0A5MgsjCtmnw87OAwNDp0HSIsyMqdIfOXeWz1REs3X0GpXQvoMGtgpjQvS4hFUs7Ozy7S0pN4/uwKD5fE8npGF0XqnZgaZ68oy53N6+GR0kbnZ1wWc9tfWEfVG2hb6/6FL0rRUfKb4LYD4xDV1x9AF1iI51Saoed4ywQkyAy2fsDLBgLKOj/CYQ+7OyICtW+M7F8ujKCX/edA8DTXRjapgZPdKtTIsYRJKfa+HFHFJ+tjkgvIBhcwZcJ3esyqFWQ03tlFaqr5/Xc1pePQkhXeHABeDhuHEtRk98EMRQYC3RGz+2QMUEopdQd9g60IEyCyCBiBXw3HGwp0OMV6PJ3Z0dUaA6ei+Pj3w+zfN95ALw83Bh+Ww0e61qHamVLOTm6wpeSZuOnnaeZtDqCE9EJANQs78uE7nUY0rp6yUkUV47DVz3h2gVoNgwGTwG3EnLuOchvguiklPpDRF5RSr3u0AjtwCQIy+kdMLMfpCRAh4nQ600oAV0fIy/G88mKI+m3krw93HiwXS0e7VqbymXMnAupaTYW7zrDZ6siOHpJ91CvWd6XZ+6sx8CWQSWjjeJMuP7bSI6Hjk9BrzecHZFLyG+C2K6UaiMiO5RSrR0aoR2YBMFfvyU1Hw6Dvyj2yeFE9DX+u/IIC3eexqb0HAv3t63BE93rmsSQhTSbYunuM/x35RGOWiVE6lby42931i8ZvZ4iV+l5T2yp0PsdaP+4syNyuvwmiM3AbmAgMC/zeqXUU/YMsqBKfIJIuAzTekH0EV3R8oHvi/V91qgrCXy6MoIFO6JIsyk83IRhoTWYeEddgkrgraS8Sk2zsTD8DJ+sOJzeRtG4ahn+3qs+dzSsVLwH3O2aqwfTleDOGxnlN0FURI+efhd4JfN6pdQsewZZUCU6QaQkwjeD4OQmqNQExizTA+KKofNxiXy2KoK5206SkqZwdxOGtAriqR71SkTjs70lp9qYH3aKT1cd4XycHkfRskZZnuvVgE51i/HkRRs+hhWv6gmHRvwIIV2cHZHTFLTURgul1C6HRGZHJTZB2GzwwxjY95OuzDpuRbEcEBSbkMLnayOZ8ccxklJtiMDAFtV4+s76xbK7amFLTElj9paTTF4dQfQ1PSdFu5DyPH9XA0KLwMjyPMs4K513gB5tXUJLchQ0QdQHPgcqK6WaikhzYIBS6k37h5p/JTZB/PYv2PgpeJexfsmLV32lxJQ0Zm08zuQ1kelF9Po0rcKzPesXiwFuruZaUiozNx7ny7WRxFklPLo1COSfvRsWyZHm2bKlwYKHYf8iPZvi+FXgF+jsqApdQRPEWuB54Msbc0CIyF6lVFO7R1oAJTJBbJ8FS54CNw8Y8UOxmk0rNc3GjztO8/GKw+kjgzvUrsALfRrSogCztxm5E3s9hWnrjzJtwzGuJachAoNbBfH3Xg2KVxtPynU9kO50GFRvC6OWgGfJ6txQ0ASxTSl1W8ZJgrKaRMjZSlyCOLZetzvYUmHAZ9D6IWdHZBdKKX7bf573lx9Kn1WtcdUy/LNPQ1Ot1Ami45P4dFUEs7ecICVN4eXhxuiOwUzoVpcA32IyFerV8/BVD4g9BU2Hwj1fFfvefxlllyByM1LkkojUwSrYZw2gO2vH+Iy8unwU5j+kk0OHicUmOWw5Gs09n2/k0W+2E3EhnhrlS/Hf4S1Z+mRnutYPNMnBCSr4efPqgCaseLYr/ZtXJTnVxpR1R+ny3iq+XBtJYkqas0MsOP/KcP9cPYfE3gWw9j1nR+QycnMFURuYAnQErgDHgAeVUiccH17ulZgriMRYPdbh0iGo10v/Yru5OzuqAjly/ipvLzvIqoN61rQKpb14qkc97m9bEy8PM9rVleyOiuGdZQfZGBkNQLUAH57t1YDBrYrBYLvDy2HOcFA2uGcaNBvq7IgKhV3mg7AmDnJTSl21Z3D2UiISRFoqzLlPl9IIbARjfyvShcei45P4eMVh5mw9RZpNUdrLnfG312Zcl9r4eRfu1J5G7iml58x+Z9lBDp7THwcNq/jzYt9GdK1fxBt5N02G5S+Ch4/u9FGtlbMjcriCtkEEAP8H3G4tWgu8rpSKzcWBewP/BdyBr5RS72Ra7w18DbQBooH7lFLHRaQt+qoFdA2oV5VSP2V3rBKRIH59ETZP1jPCjV8F5UOcHVG+JKWmMfOP43y2KoKrSam4CTzQribP3Fmfin7ezg7PyKU0m2LhztN89Pvh9Mqx3RsE8nK/xtSt5Ofk6PJJKVjyNOyYBQE14JE1xb5EeEETxA/AXuDGwLiHgBZKqWyHH4qIO3AY6AlEoQv+3a+U2p9hmyeA5kqpx0RkODBYKXWfiPgCyUqpVBGpCuwCqimlUm91vGKfIG6M/nTz1JP+BHdydkR5ppRi2d5zvL3sAKcu6w+Ubg0CealvI+qbLqtFVmJKGjM36oQfn5SKu5vwUPtaPN2jHuVKF8HR/KlJumZT1DYI7gIPLQT34ntFW9AEcVOPpdz0YhKRDuhv/ndZr18EUEq9nWGb5dY2m0TEAzgHBKoMQYlICLAZCCqxCeLsbpjWE1ITof/HEDrG2RHl2a5TMbz58362Hb8CQP3Kfrzcr3HRvyVhpLt4NYmPfj/MvG0nsSkIKOXJM3fWY0T7WkWvamzcWZjSFeLPQ/sJ0Ps/zo7IYQrai+m6iHTOsLNOwPVcvC8IOJXhdZS1LMttrA//WKCCdZx2IrIP2AM8llVyEJFHRCRMRMIuXryYi5CKoITLMG+ETg6tRkCbojWvw7nYRP42L5yBk/5g2/ErVCjtxVuDm/LLU11McihmAv29eXtIM35+qgsd61Qg9noKry3Zz12frGPVwfPktr3TJZSpCvd+ra/YN0+CXTeVoysRcnPd9Dgwy2qLAN2TabTDIrIopbYATUSkkXX8ZUqpxEzbTMFqqwgNDS1Cv325ZLPBj49AzAmo2hL6flhk+mcnpaYxfcNxPl11hITkNLzc3RjTOYQnutehjE8x6T9vZKlR1TLMHteOFQcu8NbP+zl68RpjZobRpV5FXunfuOiMgK/ZHvq8q2dlXPIUVGqoZ6UrQXJMEEqpcKCFiJSxXsflct+ngRoZXle3lmW1TZR1iykA3Vid8fgHRCQeaAoU03tIt7D2HYj4XTdK3/dNkRnhuebQBV5fsj993oHeTarwcr9GppheCSIi9Gxcma71A/l603H+u/II649cos9/1zOmcwhP9ahXNHqqhY6BMzth5zcwdwQ8uhZ8i2FtqlvI8RaTiPxHRMoqpeKUUnEiUk5EclOHaRtQT0RCRMQLGA4szrTNYmCU9XwosEoppaz3eFjHrwU0BI7n8pyKh0O/wtp3Qdxg6DQoW9PZEeXo1OUExn8dxugZ2zh66Rp1Akvzzdi2fPFQG5McSigvDzfGdanNmue68UC7mqQpxZR1R+nx4RoWhZ92/dtOItD3AwhqA7En4afH9JV9CZGbRur0EhsZluVqEiER6Qt8gu7mOl0p9ZaIvA6EKaUWi4gP8A3QCrgMDFdKHRWRh4AXgBTAhu5WuzC7YxWrRuqYk/BFZz0orghMGZqYksbnayL5Ym0kSak2Snu58/Sd9RjdMcQMdDP+YndUDP9etI9dp2IAaF+7PK8PbOr6vdhiTsIXXSAxBu58FTr/zdkR2U1BezHtBm5TSiVZr0uhP+BdqmxosUkQaSkwoy9EbYX6feD+OS7b7nCjbtIbS/enTzozqGU1XuzbyMzmZtySzaaYH3aKd389yJWEFNzdhIc7BvP0nfXwd+X2qUO/6oGq4q6L+hXBruZZKWgvptnAShEZKyJjgd/5c0yEYW+r39LJoUwQDJrsssnhdMx1xn8dxqPfbCfqynUaVvFn3iPt+WR4K5McjGy5uQnD29Zk9XPdGNG+Jjal+GrDMe74cC2Ld51x3dtODXpDp2dApcGCMXDtkrMjcrhcldqwRkTfab38XSm13KFR5UOxuIKIWAnfDtHtDqN/hlodnR3RTVLTbMzceJyPfj9MQnIaft4ePNerPiPa18KjqPV1N1zCnqhY/r1oL+HWbaeu9QN5c1BT12y3SkuFWXfDyY1Q7y54YJ7LfonLLbvUYnJ1RT5BXD0PX3SCaxeh+8vQ9R/Ojugmu07F8OKPe9h/Vndk69esKq/c3dhcMRgFZrMp5oWd4u1fDhCXmIqPpxvP9qzPmE4hrvfFIzYKPu+o2wj7vA/tHnF2RAViEoSrs6XBN4Ph2FoIuV0P7XehCq1XE1P4YPkhvt58AqUgqGwp3hjUhDsaVnZ2aEYxc+FqIm8sPcCSXWcAPRfI20Oaud4kUfsWwvejwN0bHlldpGdyLGgbhOFof/xXJwffijB4ikslh+X7znHnR2uZtekEbiI8enttfn/2dpMcDIeo5O/Dp/e3YsbDtxFUthT7z8YxePIfvLZkH/FJt6y0U/iaDILWIyEtCRaM1TPTFUPmCsLZzu2BKd3BlgIPLoB6PZ0dEQCX4pP4v8X7+Hm3nhuqZY2y/GdwMxpXK7rlxY2iJSE5lU9WHGHahmOk2RTVAnx4Y1BTejRykS8nydfgy64QfQRuGwf9PnR2RPlSoCsIEekkIr+LyGEROSoix0TkqP3DLIFSk6yBNyl6xKYLJAelFIt3naHXx+v4efdZfL3cefXuxvzweEeTHIxC5evlwUt9G7F4YieaVw/gTGwiY2eF8ey8cGITUpwdHniV1oNY3Txh21dwaJmzI7K73IyDOAj8DdgOpM8vqJSKvuWbnKBIXkGseBU2fAzlQuCxDeDt3Br6F+ISeXnhXn7ffx6ATnUr8M6Q5q7Zm8QoUdJsihl/HOOD3w6RmGKjkr83/xncjDsbu8DVxMZP4bd/gV9leGJzkSvFUdCBcluUUu0cEpkdFbkEcXILzOitJygZ86suDOYkSikWbI/ijaX7iUtMxd/bg5f7NeK+22qYeaANl3Ls0jWe/34XYSd02fghrYJ45e7GlPV14rwTNhvM7AsnN0Gze+Geqc6LJR8KmiDeQZfK+BFIurFcKbXDnkEWVJFKEMnXdCmNy0eh09PQ83WnhXIhLpF//rCb1Yd0ufTuDQJ5a3AzqpUt5bSYDCM7aTbFzI3HeX/5QRJTbLrMuLOvJqIj4fNOkHodhn8HDfs5L5Y8KmiCWJ3FYqWUusMewdlLkUoQP/9d37Os1FhPaejhnGk2l+05y0s/7eFKQgplfDx4dUATBrcKMlcNRpGQ+WpicKsgXh3QhIBSTirXsfkL+PWfULoSTNhSZG41mXEQriRylR7z4Oap55Wu2rzQQ4hLTOHVRfv4caeuvt6lXkXeH9qCKgFmwJtRtGS+mqgW4MOH97akQ50KhR+MzQaz+sOJP6DZMLjnq8KPIR8K2ospQEQ+ujFzm4h8mGHyICMvkuJh8dP6ebd/OiU5bIqMps8n6/lx52l8PN14fWATvh7T1iQHo0hydxPGdg5h2dO306JGWc7EJvLAV5t5e9kBklMLuSy3mxsM/Aw8fWHP93BgaeEe3wFyM1BuOnAVuNd6xAEzHBlUsbX6LV1Tvkpz6FS45YITU9J4c+l+HvhqM6djrtOiegA/P9WFkR2CzS0lo8gLqViaBY914Kk76iLAl2uPMmjSH0RcuFq4gZSvrcuBAyx9Rk8ZXITlpg0iXCnVMqdlzubyt5iiwuCrO3UhvvGroFrLQjt0xIWrTPxuJwfPXcXdTXjyjrpM6F636E0kbxi5sP3EZZ6ZF86py9fx9nDj5X6NeKh9rcL7ImSz6YJ+JzZAiwdg8OeFc9x8Kmipjesi0jnDzjoBxXNcuaOkpcDipwAFHScWWnJQStfdv/vTPzh47irBFXz54fGOPHNnfZMcjGKrTa3y/PJUF4a2qU5Sqo1XFu3j4ZnbiI5PyvnN9uDmBgP+p+s07foOjq4pnOM6QG4+JR4HJonIcRE5AXwGPObYsIqZLV/ChX1QLhi6vlAoh7yamMIz88L5x4LdXE9JY3CrIJY+1YWWrlb0zDAcwN/Hkw+GtWDyg60JKOXJmkMX6fe/DWw9Vki3fCrU+bMi88/PQWpy4RzXznJMEEqpcKVUC6A50Ewp1UoptcvxoRUTsadhzdv6eZ/3wcvxo5J3R8XQ/9MNLAo/QylPdz4Y1oKP72tZNCaJNww76tusKsue7kKbWuU4F5fI/VM3M2l1BDZbIfTe7PgUVKirazVtnuz44znALdsgRGSEUupbEXk2q/VKqY8cGlkeuWwbxPxRsH8hNOwPw2c79FBKKab/cZx3lh0gJU3RqGoZPnugFXUCnVvCwzCcLSXNxoe/HeaLtZEA3F4/kI/vbUEFPwePQboxCZhnaZi4DQKCHHu8fMhvG0Rp66f/LR5GTiJW6uTg6Qu933HooeKTUpn43U7eWLqflDTFqA61+OmJjiY5GAbg6e7GC30aMmP0bZTz9WTd4Yv0/d96thx1cEm5uj2g0d2Qcg1+/7djj+UAZqCco6Qmw+T2cDkS7nwNOj/jsENFXoznsW+2c+RCPH7eHnwwrDm9m1Z12PEMoyg7G3udJ7/bSdiJK7gJPHdXAx7vWsdxvZxiTsJnbXUZjlFL9KRgLqSgA+XeE5EyIuIpIitF5KKIjLB/mMXMtqk6OVSsD+2fcNhhlu87x8DP/uDIhXjqVfJj0cROJjkYRjaqBpRi7iPteaJbHWwK3vv1EE/M3uG4CYnK1oQuf9fPf/mH7tVYROSmF1MvpVQc0B84DtQFnndkUEVewmVY+65+3utN8LB/pck0m+L95Qd59JvtxCel0rdZFX6a0MncUjKMXPBwd+MfvRvy1chQ/L09WLb3HIMn/cGxS9ccc8COT+qy/hcPwNYpjjmGA+QmQdzo+tIP+F4pFevAeIqHte/pCc1DukK9Xnbf/ZVryYyesZVJqyNxE3ipb0MmPdDa9FIyjDy6s3FlFk3sRN1Kfhy5EM+Azzaw8sB5+x/I0wf6WF8aV78N8RfsfwwHyE2CWGpNGtQGWCkigUCiY8Mqwi5F6NtLCNz1Ftj5vua+M7Hc/dkG1h+5RPnSXnw7th2P3O7A+6eGUczVDvRj4YRO9G5ShauJqYydFcZ/VxzB7u2z9e+CendB8lVY49hOK/aSm3EQLwAdgVClVApwDRiYm52LSG8ROSQiESJy0wgxEfEWkXnW+i0iEmwt7yki20Vkj/XTpUqLZ2vF/4EtFVo9CFWa2XXXv+8/z9DPNxF1RddSWvpkZzrWrWjXYxhGSeTn7cHnI1rz/F0NEIGPVxzmqbnhJKak5fzmvOj5ui63s30mXDxs3307QG4aqUei2x8etJ4PBXK8byIi7sAkoA/QGLhfRBpn2mwscEUpVRf4GLCuwbgE3K2UagaMAr7J3ek42fENcHCp7vPc/V92261SiinrInnkmzCup6QxpFUQ8x7tYCb1MQw7EhEmdK/L9FG34eftwZJdZxg+ZTMXrtrxhkmlhtB6JKg0/WXSxeXmFtNtGR5dgFeBAbl4X1sgQil1VCmVDMzl5iuPgcAs6/kCoIeIiFJqp1LqjLV8H1BKRJwzq05u2Wyw/CX9vNPTUMY+PYmSU228+OMe/vPLQZSC5+9qwIf3tsDH090u+zcM46+6N6zEgsc7EFS2FOGnYhg8aSMHzsbZ7wDdXtJfIg/9Aic22W+/DpCbW0xPZniMB1oDuekqEwScyvA6ylqW5TZKqVQgFsg808c9wA6lVCFV2sqnfT/C2V3gX1UX5LOD2IQURk3fytxtp/D2cGPyg62Z0L2uaW8wDAdrWKUMCyd0olXNspyOuc7Qzzey6qCdGq/9K0OHCfr5jTI8Lio/JT2vASH2DiQrItIEfdvp0Vusf+TGREYXL14sjJCylpb65390txfAq3T22+fCqcsJDP78DzYdjSbQ35v5j3agbzMzvsEwCkugvzdzxrdnQItqXEtOY9ysMOZsPWmfnXd4ArzLwLG1cGKjffbpALlpg1giIoutx1LgEPBTLvZ9GqiR4XV1a1mW24iIBxAARFuvq1vHGamUiszqAEqpKUqpUKVUaGBgYC5CcpDdcyE6QvdzbvlggXd34Gwc93y+kaMXr9Gwij+LJnSihanCahiFzsfTnf8Ob8lTd9TFpuDFH/fwv5V26OFUqhy0f1w/d+EeTbfsOC8i3tZtnQ8yLE4FTiilonKx721APREJQSeC4cADmbZZjG6E3oRu/F6llFIiUhb4GXhBKfVHbk/GKdJS/hwU1/0lcC/YhOmbj0Yz/uswriam0r52eaaMDKWMj5MmYTcMAxHh2V4NqFTGh38v2stHvx/mwtVEXhvQFHe3Atzubf84bP78z6uIWh3tF7SdZHcFcaP1ZJxSaq31+COXyeFGm8JEYDlwAJivlNonIq+LyI1G7mlABRGJAJ4FbnSFnYgesf2KiIRbj0p5PLfCsXu+rrVSsT40vadAu/p171lGTt/K1UQ9Mnrmw21NcjAMFzGifS0+f7A1Xh5ufLv5JBO/21GwbrBF4Coiu3Lfe4H/AG+QRWkNpdSPjg0tb5xSrM+WBp/dpmsuDf4SWgzP966+23KSfy3cg03BQ+1r8eqAJgX7dmIYhkNsPhrN+FlhXE1KpV1IeaaOKsBV/vUr8ElzSIqDh3+FWh3sG2wu5LdY32Pobq1lgbszPfrbOcaiaf9CnRzK1oKmQ/O9m2kbjvHSTzo5PNuzPq8PNMnBMFxV+9oVmP9YByr5e7Pl2GUenLqFmIR8zhiX8SpiretdRWSXIKoqpR4HXlRKPZzpMaawAnRZNhuss5pnOv8N3PNXB+nLtZG8sXQ/AG8MbMJTPeqZbqyG4eIaVS3DD493pGZ5X/acjuX+qVvyP+d1+8d1j6aja1xuXER2CeJF66eZfzorh5fBhf3gXw1aZm57z51JqyN4e9lBRODtIc14qEOwfWM0DMNhapT3Zd6j7aldsTQHzsblf9T1X9oiXGtcRHYJIlpEfgNCMnRzTX8UVoAua8PH+menp8Ejb4O8lVJ8suIw7y8/hAi8d09z7m9b0wFBGobhSFUDSjH30fbUs6rBDv9yM+di85EkblxFHFsLZ8LtHmd+ZZcg+gGvoOsifZjFo+Q6tQ2itoFPWWj9UJ7eqpTiw98O88mKI7gJfHRvC4aF1sj5jYZhuKRK/j7MfaQ9jaqW4eilawyfsonzcXlMEqXK6RpNAJsn2z/IfLplglBKJSulNgMdM3RzTX8UYoyu58Z/YJvReR41/dHvh/lsdQTubsInw1sxuFV1+8dnGEahquDnzZzx7WgaVIbj0Qk8+FU+2iTaPqIrve79AeLOOibQPLplghCRT6yn080tpgxio2D/IhB3/R+aB1PWRfLpKp0c/je8FQNaVHNQkIZhFLayvl58PaYd9Sv7EXEhnhHTthKbkIfpRcvVgkZ36+kCtk11XKB5kN0tphsltj/A3GL609apulRvk0EQkLn24K3N2XqS//xyEID3hzanX3NTV8kwipvypb34dlw7QqyG65EztnI1MQ9Jor1VxC9sOiQnOCbIPMjuFtN26+daYD+wv8TfYkq+pif6AGj/RK7ftnT3GV76aQ8Arw1owpDW5raSYRRXlfx9mD2uHdXLlWLXqRge/WY7yam23L25RlsIaqMH0O2e69hAcyHbYn0i8qqIXEIX6DssIhdF5JXCCc0F7VkAiTEQFArVsxx4eJM1hy7wt3nhKAXP9arPqI7BDg3RMAznq1a2FN+Na09FP282Rkbzwo+7c1fgTwTaWSMLbnwZdaLs2iCeBToBtymlyiulygHtgE4i8rfCCtCl7LTuut02Nleb7zsTy4TZO0hJU4zvEsKE7nUdGJxhGK6kZgVfpo8OpZSnOz/uOM3HK47k7o2NBuheTWd3wZmdjg0yB9ldQTwE3K+UOnZjgVLqKDACGOnowFzOhYO6a6uXPzTOeUrus7HXGTszjGvJaQxsWY2X+jYyI6QNo4RpXr0snz3QCjeB/608wrxtuZhPwtMHWtyvn2+flf22DpZdgvBUSl3KvFApdREoeSVGb1w9NLsnx66t8UmpjJ0Zxrm4RNoGl+e9oc1NcjCMEqpHo8q8MagpAC/9tJe1h3MxuVnrUfrnnu8hKd6B0WUvuwSRXfWpfFamKqJSk2HXHP28VfYXT6lpNp78bgf7z8YRUrE0Xz7UBm8PM3+0YZRkD7arxePd6pBmUzzx7Xb2nYnN/g2VGkLNDpAcr8dFOEl2CaKFiMRl8bgKNCusAF3C4WWQEA2VGkNQ62w3fXvZQVYfukg5X09mjL6NcqW9CilIwzBc2fO9GqRPXzpm5jYu5DTaus1o/XPH1w6P7Vay6+bqrpQqk8XDXylVsm4x7bBuL7V6SPcyuIVF4aeZtuEYnu7ClJGhBFcs+NzUhmEUD25uwvvDmnNbcDnOxyXx2LfbSUrNZsKhRgPAyw9Oh0F0lrMuO1yOc1KXeHFnIHIluHlC8/tuudn+M3H884fdALxydxNuCy5fWBEahlFEeHu4M/nBNlQN8GHHyRjeXXbo1ht7+eqR1aDbIpzAJIic7PsJlA3q3wWlK2S5SUxCMo9+G0Ziio1hbaozop2pzGoYRtYC/b2Z/GBrPNyE6X8c4/f952+9cbNh+ufu+ZCbcRR2ZhJETvYt1D9vMd+0UornF+zm1OXrNAsK4I1BTU2PJcMwstWqZjn+0bsBAM8v2MWZmOtZbxjSFUpX0jNXntlRiBFqJkFkJzYKoraCRyl9BZGFb7ec5Pf95/H38WDyg63x8TQ9lgzDyNm4zrXp3iCQmIQUnpqzk9S0LMpxuHv8+eV0z4LCDRCTILK3f5H+Wb9XlmMfDp+/ypvWdKFvD2lGjfK+hRmdYRhFmJub8OG9LalcxpuwE1f45FYjrZsO0T8PLC3020wmQWRn30/6Z5PBN61KTEnjqTk7SUrV7Q79m5vS3YZh5E350l78d7geaT1pTQQbjtw0NlnXfvOrDLEn4dzuQo3PJIhbiY3SpTU8SkG9XjetfvfXgxw8d5XaFUvz6oAmTgjQMIzioH3tCjzdoz5KwTPzwrl8LdM4ZDc3aNBXPz+wtFBjMwniVo78pn/W7XHT7aVtxy8zc+NxPNyE/w5vRWlvDycEaBhGcTHxjrq0CynPpfgkXluy7+YNGvXXPw+aBOEajvyuf2a6ekhMSeOfP+xGKXi8Wx2aVQ9wQnCGYRQn7m7Ce0ObU8rTnUXhZ1iRuetr8O3gHQAX9hfqoDmTILKSmgRHrTmR6vX8y6r/rTzC0YvXqBNYmol3mPLdhmHYR60KpXnuLt319eWFe4i9nmEmOg8vfTcDIGJlocXk0AQhIr1F5JCIRIjIC1ms9xaRedb6LSISbC2vICKrRSReRD5zZIxZOrERUq5B5WZQ5s/G5/1n4vhy3VFE4L2hLUwRPsMw7Gp0x2Ba1yzL+bgkPvwt0yjrGwni6OpCi8dhCUJE3IFJQB+gMXC/iDTOtNlY4IpSqi7wMfCutTwR+DfwnKPiy1b67aU70xcppXhtyT7SbIqR7WvRplY5p4RmGEbx5e4mvD2kOe5uwrebT7D/TNyfK2t31z+PrYO0PMxzXQCOvIJoC0QopY4qpZKBuUDmmXYGAjdmxFgA9BARUUpdU0ptQCeKwnejgTpD+8PyfefYcuwy5Xw9ebZnA6eEZRhG8degij8Pta+FTcGri/f9OVVpQBBUbKBLgEdtK5RYHJkggoBTGV5HWcuy3EYplQrEAlkXPMqCiDwiImEiEnbxYi4m4ciNuLMQfUTPHFe9LaAbpt/65QAAz/asT4BvySpmaxhG4fpbz/pUKO3F1uOXWbzrzJ8r6tyhf0auKpQ4inQjtVJqilIqVCkVGhgYaJ+dnvhD/6zZXg9zB2b8cZxTl69Tv7If97c1hfgMw3CsgFKe6bWa3l12kMQUqyx4MUoQp4EaGV5Xt5ZluY2IeAABQLQDY8rZ8fX6Z3AnAC7FJzFpdQQA/+7fGA/3Ip1TDcMoIoa1qUHDKv6ciU1k9hZrLutaHUHc4Ew4JF9zeAyO/LTbBtQTkRAR8QKGA4szbbMYsCZfZSiwSikn1LTN6Lh1BVGrMwBfro0kPimVOxpWoks9O12lGIZh5MDNTXje6vY6aXUE8Ump4O0HlZuCSoPTjq/u6rAEYbUpTASWAweA+UqpfSLyuogMsDabBlQQkQjgWSC9K6yIHAc+AkaLSFQWPaDs79ol3f7g6QvVWnIpPolvNp8A4G931nf44Q3DMDK6o2El2tQqx+VryczYcEwvrNle/zy12eHHd+j9EqXUL0qp+kqpOkqpt6xlryilFlvPE5VSw5RSdZVSbZVSRzO8N1gpVV4p5aeUqq6U2u/IWIE/M3K1VuDuydR1R0lMsdGjYSUzYtowjEInIvy9l/5yOv2PYyQkp0KNdnrlqa0OP74pIpTR6e36Z7VWXLmWnH718PSd9ZwYVNGQkpJCVFQUiYnO6ZlsFF0+Pj5Ur14dT0/TOzArHWpXoGWNsoSfimH+tlOMbpIhQdhsupifg5gEkdGNGZuCWjNn20kSktO4vX4gzauXdWpYRUFUVBT+/v4EBwebGfWMXFNKER0dTVRUFCEhIc4OxyWJCI93q8Oj32xn6vpjPNi+G57+1eDqGYiOgEDH3f42XXJuUCr9FlNqlVZ8u0lfPYzpFOzEoIqOxMREKlSoYJKDkSciQoUKFcyVZw56NqpM3Up+nI65zq97z0G1lnrF+T0OPa5JEDfEnoKES1CqPL+f8eFMbCIhFUtzu+m5lGsmORj5YX5vcubmJozqUAuA77ac1D2ZAM6ZBFE4zu3VP6u2YKZ19TCqQy3c3Mwvr2EYzjewVRClPN3ZdDSac75Wu6hJEIXk0mEA4vzrsOXYZXy93LmnTXUnB2XkhZ+fn0P337dvX2JiYoiJiWHy5Ml5fv+aNWvo319P/LJ48WLeeeedfMeSm3N99dVX+eCDDwB45ZVXWLFiBQDr16+nSZMmtGzZkuvXr/P888/TpEkTnn/++TzHcePfxHC8Mj6eDGihq0svOG0VC3VwgjCN1Ddc0hOGb72qS0H1bloFfx/Tq8L40y+//ALA8ePHmTx5Mk888US+9zVgwAAGDBiQ84Z28vrrr6c/nz17Ni+++CIjRowAYMqUKVy+fBl397yXr7/xb2IUjgfa1WRe2Clm7LMxwcsfiT8P8RfBzzG3wk2CuMG6glh8Wk8vek9rc/WQX8Ev/OyQ/R5/p1+e3xMeHs5jjz1GQkICderUYfr06ZQrV45u3brRrl07Vq9eTUxMDNOmTaNLly4kJCQwevRo9u7dS4MGDThz5gyTJk0iNDSU4OBgwsLCeOGFF4iMjKRly5b07NmTfv368cEHH7B0qZ4OcuLEiYSGhjJ69Gh+/fVXnnnmGXx9fencuXN6XDNnziQsLIzPPvuM0aNHU6ZMGcLCwjh37hzvvfceQ4cOJT4+noEDB3LlyhVSUlJ48803GTgwc0Hkv3rrrbeYNWsWlSpVokaNGrRp0waA0aNH079/f2JiYpg/fz7Lly9n2bJlXL16lfj4eNq0acOLL77IsmXL6N+/P0OHDgX0lUp8fDxnz57lvvvuIy4ujtTUVD7//HO6dOmS/m9SsWJFPvroI6ZPnw7AuHHjeOaZZzh+/Dh9+vShc+fObNy4kaCgIBYtWkSpUqXy/H9pQPPqAdQOLM3Ri9eIDwrBP3q37snkoARhbjGB7sFkJYhNsRWoXMab9rVzXVTWcGEjR47k3XffZffu3TRr1ozXXnstfV1qaipbt27lk08+SV8+efJkypUrx/79+3njjTfYvn37Tft85513qFOnDuHh4bz//vu3PHZiYiLjx49nyZIlbN++nXPnzt1y27Nnz7JhwwaWLl3KCy/oggI+Pj789NNP7Nixg9WrV/P3v/+d7CrRbN++nblz5xIeHs4vv/zCtm03l4QeN24cAwYM4P3332f27NksXryYUqVKER4ezn333XfLfX/33XfcddddhIeHs2vXLlq2bHnTsWfMmMGWLVvYvHkzU6dOZefOnQAcOXKECRMmsG/fPsqWLcsPP/xwy+MY2RMR+jerCkBEamW98LLjpiA1VxAACdGQGEOSe2kuUpYRjSvjbhqn8y0/3/QdITY2lpiYGLp27QrAqFGjGDZsWPr6IUOGANCmTRuOHz8OwIYNG3j66acBaNq0Kc2bN8/38Q8ePEhISAj16ukGxREjRjBlypQstx00aBBubm40btyY8+f1fMRKKV566SXWrVuHm5sbp0+f5vz581SpUiXLfaxfv57Bgwfj6+sLYNdbWLfddhtjxowhJSWFQYMG3ZQgNmzYwODBgyldWl+BDxkyhPXr1zNgwABCQkLSt8/4b23kT7/m1fjfqgg2x5WjFegrCAcxVxCQfvVwjCBA6Nk46z9Ao3jx9vYGwN3dndTU1Hzvx8PDA5vNlv46P336b8QCpF8lzJ49m4sXL7J9+3bCw8OpXLmyw8cLZDwXm81GcnIyALfffjvr1q0jKCiI0aNH8/XXX+d6nxnPraD/1gbUr+xHncDSHEiybitFO+4KwiQISE8Q+5Ir4eftQfva5Z0ckGEPAQEBlCtXjvXrdQn3b775Jv1q4lY6derE/PnzAdi/fz979tzcS8Tf35+rV6+mv65Vqxb79+8nKSmJmJgYVq7Uk8o3bNiQ48ePExmp/4DnzJmTp/hjY2OpVKkSnp6erF69mhMnTmS7/e23387ChQu5fv06V69eZcmSJXk6HkBwcHD6bbXFixeTkqKntjxx4gSVK1dm/PjxjBs3jh07/lpJtEuXLixcuJCEhASuXbvGTz/9RJcuXfJ8fCNnIsIdDStxTOlbTVw+mv0bCsDcYgK4chyAE7YqdKpfAW+PvPfmMJwvISGB6tX/7Fzw7LPPMmvWrPRG6tq1azNjxoxs9/HEE08watQoGjduTMOGDWnSpAkBAX8t1FihQgU6depE06ZN6dOnD++//z733nsvTZs2JSQkhFatWgG6DWHKlCn069cPX19funTp8pfEkpMHH3yQu+++m2bNmhEaGkrDhg2z3b5169bcd999tGjRgkqVKnHbbbfl+lg3jB8/noEDB9KiRQt69+6dfstozZo1vP/++3h6euLn53fTFUTr1q0ZPXo0bdvqWRjHjRtHq1atzO0kB7m9fiDz11fSL65k/8WhIMTZ0y/YS2hoqAoLC8vfm398BHbP4/mUR2jQ+3HGdalt3+BKgAMHDtCoUSNnh1FgaWlppKSk4OPjQ2RkJHfeeSeHDh3Cy8vL2aEVa8Xl96ewJKak0fL15YS7jcRHUuDFKPD2z9e+RGS7Uio0q3XmCgIgVk90d0ZV4IFa5ZwcjOFMCQkJdO/enZSUFJRSTJ482SQHw+X4eLrTvHo5zp0uT7Cch6vn8p0gsmMSBJAWexp34LJbRZpUM/M+lGT+/v7k+0rUMApRq5plOX+6HMGch6tnoaL9pyUwjdSALeEyANWCquPlYf5JDMNwfa1qlOW8su54xJ11yDHMp6EtDY/kOGxKqF+rhrOjMQzDyJUWNcpyQZUFQF299SDMgjAJIjEWQRGHLw2rlXV2NIZhGLlSpYwPCe76lnhi3CWHHMMkCOv2Uozyo06gY6uBGoZh2IuI4OWnbzFdjTEJwiHUdStB4EeN8r5OjsbIj+joaFq2bEnLli2pUqUKQUFB6a9vjAQubkwpbwPAN0CPpk68Gu2Q/Zf4XkzXYy/iC8SJP2V8Svw/R5FUoUIFwsPDAf3B6efnx3PPPZe+PjU1FQ8P5//fOioOU8q75PIrWxHOgC0hxiH7d/5fjZPFXonGF0j19DdTH9rLqw7qKvxqbK43HT16ND4+PuzcuZNOnToxfPhwnn76aRITEylVqhQzZsygQYMGzJw5k8WLF5OQkEBkZCSDBw/mvffeIy0tjbFjxxIWFoaIMGbMGPr06cPIkSPZunUroOeFuPvuu9mzZw/bt2/n2WefJT4+nooVKzJz5kyqVq1Kt27daNmyJRs2bOD++++nZs2avPbaa7i7uxMQEMC6detIS0vjhRdeYM2aNSQlJTFhwgQeffTRm87JlPI2MvP317fFbSmOqdFV4hNEXHw8VQE3L/NLXdxERUWxceNG3N3diYuLY/369Xh4eLBixQpeeuml9LLT4eHh7Ny5E29vbxo0aMCTTz7JhQsXOH36NHv36qloY2JiKFu2LMnJyRw7doyQkBDmzZvHfffdR0pKCk8++SSLFi0iMDCQefPm8fLLL6d/oCYnJ6ePrWjWrBnLly8nKCgo/fbNtGnTCAgIYNu2bSQlJdGpUyd69epFSEhI+rlkLOWdmppK69at0xPEDePGjWPDhg03JYEbV1fLli3L8t/pRinvl19+mbS0NBISEv6yPmMpb6UU7dq1o2vXrpQrV44jR44wZ84cpk6dyr333ssPP/yQfvViOF6ANbOgSk1yyP5LfIK4nhAPgLuXaX+wmzx803ekYcOGpd9aiY2NZdSoURw5cgQRSS9CB9CjR4/0ekuNGzfmxIkTNGnShKNHj/Lkk0/Sr18/evXqBcC9997LvHnzeOGFF5g3bx7z5s3j0KFD7N27l549ewK6XEfVqlXT959xnoVOnToxevRo7r333vRy47/99hu7d+9mwYIF6bEeOXLkLwnClPI2suJbSv8+uNsc09bm0EZqEektIodEJEJEXshivbeIzLPWbxGR4AzrXrSWHxKRuxwVY1rydf3Ewzv7DY0i58YHGsC///1vunfvzt69e1myZMlfymZnVY66XLly7Nq1i27duvHFF18wbtw4QH/Yz58/n8OHDyMi1KtXD6UUTZo0ITw8nPDwcPbs2cNvv/2WZRxffPEFb775JqdOnaJNmzZER0ejlOLTTz9Nf/+xY8fSE5I9mVLexY+Xtw8A7qqIJQgRcQcmAX2AxsD9ItI402ZjgStKqbrAx8C71nsbA8OBJkBvYLK1P7u77FGZ1WktiPWr44jdGy4iNjaWoKAgQE/3mZNLly5hs9m45557ePPNN9PLW9epUwd3d3feeOON9CuDBg0acPHiRTZt2gRASkoK+/bty3K/kZGRtGvXjtdff53AwEBOnTrFXXfdxeeff55+VXP48GGuXbv2l/eZUt5GVrx8rARhS8lhy/xx5C2mtkCEUuoogIjMBQYC+zNsMxB41Xq+APhMdEvxQGCuUioJOCYiEdb+Ntk7yNKth7KmdFfahZg5IIqzf/zjH4waNYo333yTfv1ynvHu9OnTPPzww+nfuN9+++30dffddx/PP/88x44dA8DLy4sFCxbw1FNPERsbS2pqKs888wxNmjS5ab/PP/88R44cQSlFjx49aNGiBc2bN+f48eO0bt0apRSBgYEsXLjwL+8zpbyNrHh767ZTT+WYBOGwct8iMhTorZQaZ71+CGinlJqYYZu91jZR1utIoB06aWxWSn1rLZ8GLFNKLch0jEeARwBq1qzZJqcJVQzHMeWajYIwvz/5c/l8FOU/b8IV/Cn3alS+9lFsy30rpaYAU0DPB+HkcAzDMApVQMVqRE04lt5YbW+OTBCngYzV76pby7LaJkpEPIAAIDqX7zUMwyjR3N3dqB7ouNvjjuzFtA2oJyIhIuKFbnRenGmbxcAo6/lQYJXS97wWA8OtXk4hQD1gqwNjNeyguMxOaBQu83vjuhx2BaGUShWRicBywB2YrpTaJyKvA2FKqcXANOAbqxH6MjqJYG03H92gnQpMUEqlOSpWo+B8fHyIjo6mQoUKZkS6kWtKKaKjo/GxeuMYrsXMSW3YRUpKClFRUX8ZX2AYueHj40P16tXx9PR0diglUrFtpDZch6en519G/hqGUfSV+HLfhmEYRtZMgjAMwzCyZBKEYRiGkaVi00gtIhcBewylrgg4Zv6+wmXOw7WY83At5jz+VEspFZjVimKTIOxFRMJu1aJflJjzcC3mPFyLOY/cMbeYDMMwjCyZBGEYhmFkySSIm01xdgB2Ys7DtZjzcC3mPHLBtEEYhmEYWTJXEIZhGEaWTIIwDMMwsmQSBCAib4jIbhEJF5HfRKSatVxE5H8iEmGtb+3sWLMjIu+LyEEr1p9EpGyGdS9a53FIRO5yYpg5EpFhIrJPRGwiEpppXZE5DwAR6W3FGiEiLzg7nrwQkekicsGa+fHGsvIi8ruIHLF+lnNmjDkRkRoislpE9lu/U09by4vaefiIyFYR2WWdx2vW8hAR2WL9fs2zplawH6VUiX8AZTI8fwr4wnreF1gGCNAe2OLsWHM4j16Ah/X8XeBd63ljYBfgDYQAkYC7s+PN5jwaAQ2ANUBohuVF7TzcrRhrA15W7I2dHVce4r8daA3szbDsPeAF6/kLN37HXPUBVAVaW8/9gcPW71FROw8B/KznnsAW6zNpPjDcWv4F8Lg9j2uuIAClVFyGl6WBGy33A4GvlbYZKCsiVQs9wFxSSv2mlEq1Xm5Gz8QH+jzmKqWSlFLHgAigrTNizA2l1AGl1KEsVhWp80DHFqGUOqqUSgbmos+hSFBKrUPP05LRQGCW9XwWMKgwY8orpdRZpdQO6/lV4AAQRNE7D6WUirdeeloPBdwBLLCW2/08TIKwiMhbInIKeBB4xVocBJzKsFmUtawoGIO++oGifR4ZFbXzKGrx5kZlpdRZ6/k5oLIzg8kLEQkGWqG/fRe58xARdxEJBy4Av6OvTmMyfCm0++9XiUkQIrJCRPZm8RgIoJR6WSlVA5gNTHRutLeW03lY27yMnolvtvMizV5uzsNwbUrf1ygS/eRFxA/4AXgm0x2DInMeSqk0pVRL9J2BtkBDRx+zxEwYpJS6M5ebzgZ+Af4POA3UyLCuurXMaXI6DxEZDfQHeli/+FAEz+MWXO48clDU4s2N8yJSVSl11rrdesHZAeVERDzRyWG2UupHa3GRO48blFIxIrIa6IC+7e1hXUXY/ferxFxBZEdE6mV4ORA4aD1fDIy0ejO1B2IzXJa6HBHpDfwDGKCUSsiwajEwXES8RSQEqAdsdUaMBVTUzmMbUM/qaeKFnnN9sZNjKqjFwCjr+ShgkRNjyZHoCdKnAQeUUh9lWFXUziPwRq9EESkF9ES3p6wGhlqb2f88nN067woP9LeLvcBuYAkQpP7sOTAJfa9vDxl61LjiA91oewoItx5fZFj3snUeh4A+zo41h/MYjL6fmgScB5YXxfOw4u2L7jkTCbzs7HjyGPsc4CyQYv1/jAUqACuBI8AKoLyz48zhHDqjbx/tzvB30bcInkdzYKd1HnuBV6zltdFfkiKA7wFvex7XlNowDMMwsmRuMRmGYRhZMgnCMAzDyJJJEIZhGEaWTIIwDMMwsmQShGEYhpElkyAMwzCMLJkEYRRJIlJWRJ4orPfZk4i8KiLPZXi9MZ/7yfO5iEiwiFy3avoUiIiUEl0iP1lEKhZ0f4brMQnCKKrKAnn9cBSgfF7fVxDWKPxs/86UUh3zufuy5O9cIpWu6VMgSqnr1n7OFHRfhmsyCcJwCSKyyvo2Gi4iiSJyb4Z1pUXkZ2uylL0ich/wDlDH2v59a7uFIrLdmlDlEWtZsOgJe75Gj0Cdlvl9meIYKXrCpV0i8k2G5c9mKCj4THbLszhmDRF5WUQOi8gG9FwXGY8Zb73ngIhMteL/zSqpcGObm87tFv8GI0RPLBMuIl+KiHsu/u1vOmcrnoMiMtOKe7aI3Ckif4ieZMeVy6wb9uLsIeTmYR4ZH8Dj6ElQ3DMsuweYmuF1ABBMholsrOXlrZ+l0B/MFaztbEB7a91N78vw/iboshgVM+2vDbrUSmnAD9iHLht9q+WZj3ljO1+gDLoswnMZjhtvvScVaGktmw+MyMW5ZZzMpxG6VIyn9XoyMDLTOWZ+z63O+UY8zdBfJLcD09HlZwYCCzPs4/iN95tH8XqUmGquhusTkZFAH+AepVRahlV7gA9F5F1gqVJqvWQ9ReRTIjLYel4DXczvHHBC6QmfcnIH8L1S6hKAUurGZDmdgZ+UUtesOH8EuqA/LLNavjjTMbtY2yVY292qYN8xpVS49Xw7+kM6p3PLqAc6GW3Td9MoRc5VSm91zjfi2WPFvA9YqZRSIrInU2xGMWUShOESRGQYerKmgUqplIzrlFKHRc8H3hd4U0RWAl9nen834E6gg1IqQUTWAD7W6muOjT5L+TlmUobnaegP+JzOLSMBZimlXszHsXOKx5bhtQ3z2VEimDYIw+lEpD+6sXWIUioxi/XVgASl1LfA++h5kq+i5xi+IQC4Yn2ANkTP15uVzO/LaBUwTEQqWMctby1fDwwSEV8RKY2uNrs+m+WZrbO2KyUi/sDdtzj+rdzq3DKfy0pgqIhUuhG/iNTKYd+3OmfDMN8CDJcwCz338R/WrZFPlVLTMqxvBrwvIjZ06enHlVLRVoPpXvTUqv8CHhORA+hS4FneUsr8PqXU8xnW7RORt4C1IpKGLq88Wim1Q0Rm8ufcE18ppXYCZLVc9NSWGY+5Q0TmAbvQt3y25fHf59eszi2rcxGRfwG/WT2nUoAJwIlb7fhW55zH+IxiypT7NowSxEpeS5VSTe24z+PouVIu2Wufhmswt5gMo2RJAwLEjgPlAE90u4RRzJgrCMMwDCNL5grCMAzDyJJJEIZhGEaWTIIwDMMwsmQShGEYhpElkyAMwzCMLJkEYRiGYWTJJAjDMAwjSyZBGIZhGFn6fxoXSc7kKk6UAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1, 1)\n",
    "ax.plot(np.sort(tracks_quenched[\"z\"]), tracks_drifted[\"long_diff\"][np.argsort(tracks_quenched[\"z\"])], label=\"Longitudinal diffusion\", lw=2)\n",
    "ax.plot(np.sort(tracks_quenched[\"z\"]), tracks_drifted[\"tran_diff\"][np.argsort(tracks_quenched[\"z\"])], label=\"Transverse diffusion\", lw=2)\n",
    "ax.set_xlabel(\"$z$ start coordinate [cm]\")\n",
    "ax.set_ylabel(\"Diffusion coefficient\")\n",
    "_ = ax.legend()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "NumbaKernel",
   "language": "python",
   "name": "numbaenv"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
