{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Interaction probability" ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "This notebook demonstrates the verification of the **nuclear interaction module** in the `GT simulation` package.\n", "The goal is to simulate the propagation of cosmic ray protons through the interstellar medium (ISM) and calculate the probability of inelastic nuclear interactions as a function of energy.\n", "\n", "The interaction probability is governed by an exponential attenuation law. For a particle traversing a grammage $X$ (column density of matter in g/cm$^2$), the probability of undergoing an inelastic collision is given by:\n", "$$\n", "P_{\\text{int}}(X) = 1 - \\exp\\left(-\\frac{X}{\\Lambda}\\right)\n", "$$\n", "where $\\Lambda$ is the nuclear interaction length.\n", "\n", "Workflow overview:\n", "- Define the ISM properties (homogeneous hydrogen gas) and the primary particle spectrum (protons, 1 GeV -- 1 TeV, log-uniform distribution).\n", "- For each particle, assign a target path length (grammage $X$) based on the empirical Path Length Distribution.\n", "- Trace particles through the medium until they either interact or pass the assigned grammage.\n", "- Calculate the fraction of interacting particles in energy bins.\n", "- Compare the simulated interaction probability with the theoretical expectation for $\\Lambda = 55 \\text{ g/cm}^2$.\n", "\n", "Expected output:\n", "A plot showing the interaction probability vs. energy, where simulation data points should align with the theoretical curve, confirming the correct implementation of the physics module." ] }, { "cell_type": "code", "execution_count": 1, "id": "2", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from gtsimulation.Global import Units" ] }, { "cell_type": "markdown", "id": "3", "metadata": {}, "source": [ "## Path length distribution (grammage)\n", "\n", "We use a parameterized function to represent the mean escape length (grammage) traversed by cosmic rays in the Galaxy as a function of their rigidity $R$. This distribution corresponds to the \"Nested Leaky Box\" model or similar empirical fits (e.g., *Simon et al., 1998*).\n", "\n", "The function below calculates the grammage $X(R)$ in g/cm$^2$." ] }, { "cell_type": "code", "execution_count": 2, "id": "4", "metadata": {}, "outputs": [], "source": [ "def grammage(R, Rb=4.8362, coeffs=(2.2089, 0.4517, -0.2111, 0.0256), k=-0.6):\n", " a, b, c, d = coeffs\n", " ln_xb = np.log(Rb)\n", " A2 = np.exp(a + b*ln_xb + c*ln_xb**2 + d*ln_xb**3) / (Rb**k)\n", " return np.piecewise(R, [R <= Rb, R > Rb], [\n", " lambda t: np.exp(a + b*np.log(t) + c*np.log(t)**2 + d*np.log(t)**3),\n", " lambda t: A2 * t**k\n", " ])" ] }, { "cell_type": "code", "execution_count": 3, "id": "5", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGvCAYAAACTjDUBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABMLklEQVR4nO3de1xT9/0/8FfCVW4JAQRRFIL3O7fW2rsE7X2dAq7rvSvQdlu3dR0p+/72bbt9VwTd2q1bV9De1m6dEN26dr0R6NXaCgSsdyUBxSuXkISLckt+fyCZCGggCScJr+fjkcfMycnJm3lqXnyuIovFYgERERGRmxALXQARERHRWDC8EBERkVtheCEiIiK3wvBCREREboXhhYiIiNwKwwsRERG5FYYXIiIicisML0RERORWvIUuwBnMZjNOnjyJ4OBgiEQiocshIiIiG1gsFrS3tyM6Ohpi8ejtKx4ZXk6ePImYmBihyyAiIqJxaGxsxIwZM0Z93SPDS3BwMICBHz4kJETgaoiIiMgWJpMJMTEx1u/x0XhkeBnsKgoJCWF4ISIicjOXG/LBAbtERETkVhheiIiIyK0wvBAREZFbYXghIiIit8LwQkRERG6F4YWIiIjcCsMLERERuRWGFyIiInIrDC9ERETkVhheiIiIyK0wvBAREZFbYXghIiIit+KRGzOS++vq6UNLew/O9fWjp8+Mnn4zevrMEItECPD1QqCfNwL9vBDs54Mpvl5Cl0tERBNI0PCi0WiQlZWF6urqYcfVajUAoLKyEps3b4ZUKhWgQnKmju4+HDxlwoFTJuw/1Y7jbV04bTyH06ZzaD/XZ9M1vMUiKG+aj6zr5E6uloiIXIVg4UWlUkEul0Oj0Qx7Ta1WIzc3FwBQWFiI1NTUYQGH3E9LRzd21LXgiyMtqGzQ42hr1yXP9/cRI8DXGz5eIvh6i+HrJYbZAnR296Grpx+dPX3oM1vw3AcHsHi6BFfFh03QT0JEREISWSwWi6AFiES4sASNRoPU1FS0tbUBAHQ6HeLj46HVaiGX2/bbtclkgkQigdFoREhIiFPqJtscOGXCv3efxGeHmrH/lGnY61Eh/lgYHYIF04IRFx6EqBB/REn8EBnij2B/n0te22KxQLntW5RUHUdUiD8++Mm1CA30ddaPQkRETmbr97fLjXlJTEzE5s2brc8NBgMAQCaTCVQRjdVJw1m8U3sS79SewMHT7UNeWzgtBNfOCcfK2eFYOl1iV9gQiUR45o5FqDraBl1zJ36h2o3N9yVDJBLZ+yMQEZELc7nwAgDp6enWP2/duhUKheKSY166u7vR3d1tfW4yDf8Nn5zLYrFgR10rtnypw2eHmzHYmObrJcaq+VNx85IoXD07HOFBfg793ABfb7x4VwK+++evoD7QhL/uPIr7V8Y69DOIiMi1uGR4GWQwGKBSqS473iU/Px/PPvvsBFVFF+rpM+Pd3Sex5ct6HLigW+jKOBm+mzAdNy+eBknApbt/7LUoWoJf3jIfz7y7H799/wBSYmVYGM3uQiIiT+VyY14ulJOTA6VSedmxLiO1vMTExHDMixP19ZtRUnUcfyw/gtOmcwCAKT5eyEyegYeuicOssMAJrcdisSDrr1VQH2hCfEQg3v/JtfDz5hRqIiJ34rZjXgYVFhZag8vguJfRuo78/Pzg5+fY7ggamcViwUf7zqDwo4PQNXcCAKYG++GBq2Nx9xWznN7KMhqRSITC9GVY/fxn0DZ34tNDzVizKEqQWoiIyLlcYoXdwXAySKVSITEx0RpcSkpKuM6LC6g+2oZ1f/kKj7xVDV1zJ2SBvnj69oX4QnkjHrthtmDBZZAs0BffWT4dAPDh3tOC1kJERM4jWMuLWq1GWVkZgIExKykpKUhPT4dOp0NGRsaQc6VSKbKzs4UokwCYzvWi4IOD+Ns3xwAMdA89fG0csq+TX3Y680S7eXEUXvmyHuoDZ9DTZ4avt0vkcyIiciDBx7w4A9d5cZwP957G0//eizOmgTFF6UkzkLtmHqaG+Atc2cjMZguuzC9Hc3s3Xn8wBTfMmyp0SUREZCO3H/NCwmpu78b/+9cefLTvDAAgNiwAz61dgpXx4QJXdmlisQhrFkXira+P4cO9pxleiIg8ENvUaZgvjjTj5j98gY/2nYG3WITHbojHhz+9zuWDy6CbF08DAHy8/wz6+s0CV0NERI7Glhey6u034/dlh/HyZ1pYLMDcyCC8sD7B7dZMuTJOBmmAD/SdPdjVoHeb0EVERLZhywsBABr1Xcgs2om/fDoQXO6+cib+/aNr3C64AIC3lxhpCyIBAB9x1hERkcdheCHsqGvBbS9+iZpjBgT7e+OluxPx2+8ugb+P+y7ydvOSgTVePtx3Gmazx41JJyKa1NhtNMm9ubMBz7y7H/1mC5bFSPGnuxIQIwsQuiy7XT07HEF+3jhj6kZNowFJs0KFLomIiByELS+TVG+/Gf/vX3vwq3f2od9swXcTpmNr9gqPCC4A4OfthVXzB2Yafbj3lMDVEBGRIzG8TEKGrh7c/+ouvPX1MYhEgPKm+fh95jK37iYayc2L/9t15IHLGRERTVrsNppkzpjO4d5XvsHhMx0I9PXCC99LQNrCSKHLcorr50XA30eMRv1Z7DtpwuLpEqFLIiIiB2DLyyTS0NKJdX/5CofPdCAyxA+qR1d6bHABgABfb1w/NwIA9zoiIvIkDC+TxIFTJqS/vBPH284iNiwAqkdWYsE095sGPVaDC9a9v+cUu46IiDwEw8skUNWgR2bRTrR0dGPBtBCUPrLSYwbmXk7qgqnw9RZD19KJg6fbhS6HiIgcgOHFw32ja8W9r+xC+7k+pMSG4h/ZKxAR7Cd0WRMm2N8HN5zvOnrv25MCV0NERI7A8OLBqhr0ePD1Spzt7cd1cyPw14euhGSKj9BlTbhblw50Hf3nW3YdERF5AoYXD1VzrA0PvFaJrp5+XDsnHMX3JmGKr2dNhbaVYkEk/LzFaGjtwr6TJqHLISIiOzG8eKA9x42479Vd6Ojuwwq5DMX3JnvcGi5jEejnbV2w7j97uGAdEZG7Y3jxMPtPmnDPK99Yx7i8cn/KpG1xudBg19F7355k1xERkZtjePEgjfou3P/aLhjP9iJxphSvPXgFAv24DiEArJo/FVN8vNCoP4s9J4xCl0NERHZgePEQ+s6BJf+b27sxPyoYrz14BYIYXKwCfL2xasH5rqNv2XVEROTOGF48wNmefvzgjUroWjoxXToFbzx0xaScVXQ5ty0Z7DrirCMiInfG8OLm+vrN+NHfNag5ZoBkig/eeCgFkSH+Qpflkm6cPxUBvl44YTiL2kaD0OUQEdE4Mby4MYvFgl+9sxflB5vg5y3GK/cnY/bUYKHLcln+Pl5QLBjYy4ldR0RE7ovhxY0Vf67D27saIRYBf7wrAcmxMqFLcnmDs47e33MKZjO7joiI3BHDi5v65GATNnx4EADwv7ctxJpFUQJX5B6unxuBID9vnDSeQ01jm9DlEBHRODC8uKG6pnY8/nYNLBbgritm4v6VsUKX5Db8fbyweuFA19G/a7nXERGRO2J4cTOGrh784I0qtHf34Yo4GZ69YxFEIpHQZbmVO5ZHAxiYddTXbxa4GiIiGiuGFzfS22/GD/+uwdHWLswInYK/3J0IX2/+FY7V1bPDERboi9bOHuzQtgpdDhERjRG/+dzIb/9zADvqWhHg64XN9yUjLMhP6JLcko+X2Dpw953aEwJXQ0REY8Xw4ibe3X0Sr3/VAAB4fv1yLJgWImxBbu4757uOPtp7Gud6+wWuhoiIxoLhxQ3UNXXgqW3fAgAeuyGeM4scIHFmKGaETkFnTz/KDzQJXQ4REY0Bw4uL6+rpw2N/q0ZnTz9WyGV4Im2u0CV5BJFIhDuWDbS+sOuIiMi9MLy4MIvFgv/5514cPtOBiGA//PGuBHh78a/MUb6zfDoA4NNDzTB29QpcDRER2YrfhC7s7V2N+GfNCXiJRfjTXQmYGsw9ixxpXlQw5kcFo6ffjA/3cbsAIiJ3wfDiovadNOKZf+8DAPxizTxcKQ8TuCLPNLjmyztcsI6IyG0wvLigsz39ePztGvT0m6FYMBXZ18qFLslj3b50ILzs1LXijOmcwNUQEZEtGF5c0G/+sx/a5k5EhvihMH0ZxGKuoOssMbIAJM0KhcUyMB2diIhcH8OLi/lo32n8/ZtjEImA32cuhyzQV+iSPN532HVERORWGF5cyGnjOSjPr+eSfa0cV88OF7iiyeHWJdPgLRZhzwkj6prahS6HiIgug+HFRZjNFjxRUgtDVy8WTw/Bz1fPE7qkSSMsyA83zIsAAGzXcM0XIiJXx/DiIjZ/ocNX2lZM8fHCH76XwA0XJ9jaxBkAgH/WnIDZbBG4GiIiuhR+Q7qAg6dN2PTxIQDA07cvRHxEkMAVTT6r5k9FiL83ThnP4Wsdd5omInJlDC8C6+034+clu9Hbb4FiwVSsT4kRuqRJyd/HC7ed3y5gG7uOiIhcGsOLwP78SR32nTRBGuCD59YugUjEadFCWZswsF3Ah3tPoaunT+BqiIhoNAwvAtp7wog/VdQBAH79ncVc/l9gSbNCMVMWgM6efny874zQ5RAR0SgYXgTS3dePJ0t3o89swS1LonD70mlClzTpiUQirE0caH3ZpjkucDVERDQahheB/LH8CA6ebkdYoC9+853F7C5yEWsTBmYd7ahr4XYBREQuiuFFALsbDfjLp1oAwG+/uxhhQX4CV0SDZoYFIHlWKMwW4J1aDtwlInJFgoYXjUaDpKSkYcd1Oh0KCwuhUqlQWFgIg8Ew8cU5SW+/Gcpt38JsAe5YFo2bFrO7yNUMrvmyrfoELBau+UJE5GoECy8qlQrAQIC5WEZGBnJzc5Geno709HRkZWVNdHlOU/y5DgdPt0MW6Itn7lgkdDk0gluXTIOvtxiHzrRj/ymT0OUQEdFFBAsv6enpSExMHHZcp9MNeS6Xy6FWqyeqLKfSNXfgD+VHAAC/um0BN110UZIAH6QtiAQAqKo5cJeIyNW43JgXtVoNmUw25JhMJhuxhcadmM0W5G3fg54+M66bG4E7l08XuiS6hPTkga6jf9WcQE+fWeBqiIjoQi4XXkYb36LX60d9T3d3N0wm05CHqympasQ39XpM8fHCb+/k7CJXd92cCESF+KOtqxflB7jmCxGRK3G58DKaSw3azc/Ph0QisT5iYlxrif0m0zn89v0DAICfr56LGFmAwBXR5XiJ/7vmS0lVo8DVEBHRhVwuvEil0mGtLHq9HlKpdNT35OXlwWg0Wh+Nja71ZfPMu/vQfq4Py2ZI8ODVcUKXQzbKSB4IwZ8dbsZpI9d8ISJyFS4XXhQKxYjHk5OTR32Pn58fQkJChjxcxScHm/D+ntPwEouQv3YpvMTsLnIXceGBuCJWBrOFK+4SEbkSlwgvF3YJyeXyIa/pdDokJydfsuXFVZ3r7cfT/94HAPjBNXFYGO06oYpsk3F+4G5pVSPXfCEichGChRe1Wg2lUglgYMzK4LovAFBaWgqlUgmVSoWioiKUlpYKVaZdXvpUi2P6LkSF+OMnqXOELofG4ZYl0xDg64WG1i5UNrQJXQ4REQEQWTzw10mTyQSJRAKj0ShYF1J9SyfWvPA5evrMeOnuRNyyhCvpuqtc1W6UVB1HRtIMbMxYJnQ5REQey9bvb5foNvI0FosFT/97n3VNl5sXRwldEtkh8/zA3f/sOYWO7j6BqyEiIoYXJ/hg72l8frgZvt5i/PqORVzTxc0lzQqFPDwQXT39eP/bU0KXQ0Q06TG8OFhHdx9+/e5+AMAj18cjNjxQ4IrIXiKRyDpteivXfCEiEhzDi4O9WH4Ep03nMFMWgMduiBe6HHKQtYnT4SUWofpoG+qa2oUuh4hoUmN4cSBdcwde3VEPAHjmjoXw9/ESuCJylMgQf6yaPxUA8PYutr4QEQmJ4cWBfvPefvT2W3DjvAismh8pdDnkYN+/YiaAgQXrzvX2C1wNEdHkxfDiIJ8cbMInh5rh4yXCr25bKHQ55ATXzY1AtMQfhq5efLTvtNDlEBFNWgwvDtDTZ8av3xsYpPvg1XGQRwQJXBE5g5dYhPUpA60vf//mmMDVEBFNXgwvDvDajnrUt3QiPMgPP141W+hyyIkyU2ZALAK+qddD29whdDlERJMSw4udmtrP4cWKOgBA7k3zEOzvI3BF5EzTJFOsA3f/sYutL0REQmB4sVPhh4fQ0d2HZTMkSE+cIXQ5NAHuOj9wV1V9HN19HLhLRDTRGF7s8O1xA1TVxwEAT9+xCGIxV9KdDK6fG4FpEn+0dfXio31nhC6HiGjS8bb1RJPJBL1eP+YPkMlkgm2O6EwWiwX/958DAIA7l0cjcWaowBXRRPH2EiMzOQZ/KD+Ct785hjuWRQtdEhHRpGJzeHnuuedwxRVXYKybUFdVVSE/P3/Mhbm6j/adwa56Pfy8xfjFTfOFLocmWGZKDF6sOIKdulbomjs4w4yIaALZHF5mz56NtWvXjvkDxtNa4+p6+szY8MFAq0vWtXJMl04RuCKaaNOlU3DDvKmoONiEt3cdw//cyrV9iIgmik1jXsrLy/Hyyy/j0UcfRUNDAwBg27ZtePTRRy/73qysLLsKdEVvfn0UDa1dCA/ywyPcv2jSuvvKgYG7JVVccZeIaCLZFF7KyspQWlqK9PR05Obmora2FuvWrYNarXZ2fS7H0NWDP5YfAQD8fPVcBPnZ3HhFHuaGeVMxI3QKjGd78e7uk0KXQ0Q0adgUXlJSUhAXF4fU1FSUlJSgrKwMNTU1kEqlTi7P9fyxvA7Gs72YHxWMzOQYocshAXmJRbj7ylkABlrjiIhoYtg8VXrTpk3WP//iF7+AXq+HVqt1SlGuqr6lE29+3QAA+OUtC+DFqdGTXmbyDPh6ifHtcSN2NxqELoeIaFKwKbysW7cOCQkJQ46lpqaiurraKUW5qg0fHEBvvwU3zIvAdXMjhC6HXEBYkB9uXToNAFtfiIgmis0tL6mpqcOOxcXFDTtmMpnsq8hFWSwWpMTKIAv0xS9vWSB0OeRC7r1qoOvo3d0n0dbZI3A1RESez+4Vdk0mExoaGqwPpVLpiLpcjkgkwsPXyvHVU6swNzJY6HLIhSTESLEoOgTdfWaUVjcKXQ4RkccTWca66twFHnnkEajV6iEDd+vr69Ha2uqI2sbNZDJBIpHAaDR65Oq+5Hr+sesYntq+BzNlAfj0yRu4VQQR0TjY+v1tV8tLfHw86urqUFVVZX1s2LDBnksSuaU7lkcj2N8bx/Rd+OxIs9DlEBF5NLvCi0KhGHYsLS3NnksSuaUAX29kJA1MnX9rJwfuEhE5k10rrIWGhmLTpk2Qy+WQSqUwGAzYunUrtm7d6qj6iNzGPStm4tUd9ag41IRjrV2YGRYgdElERB7JrvCSm5sLg8EwZMxLTU2NvTURuSV5RBCunROOL4604I2dDfjVbdzviIjIGewKL2lpacP2Ltq2bZtdBRG5s4eujsMXR1pQUtmIn6Vx+wgiImewe8CuLceIJovr50ZAHh6I9u4+bNccF7ocIiKPZNevhVqtFkVFRUhJSQEwsJBbSUkJKisrHVIckbsRi0W4f2Usnv73Pry+owH3XDmL06aJiBzMrpaXoqIixMXFwWKxYHC5GDuWjSHyCOuSZiDYzxu6lk5OmyYicgK7Wl4KCgqGbRsw0vRposkkyM8bmSkxeOXLery2owE3zpsqdElERB7FrpaX5ORkbNq0ybqfUUVFBce8EAG4/6pYiETA54ebUdfUIXQ5REQexa7wUlJSgpaWFuvzVatWQa1W210UkbubGRYAxYJIAMDrX9ULXA0RkWexK7yEhYVhw4YN3D+IaAQPXh0LANhWfQLGrl5hiyEi8iB2hZddu3ahvb19yDHONCIacJU8DPOjgnG2tx9bq44JXQ4RkcewK7zk5OQgISEBa9aswfr16zFnzhzubUR0nkgkwkNXxwEAXt/RgN5+s8AVERF5BrvCS1xcHKqrq5Geno7k5GR8/PHHWLVqlaNqI3J7dyyPRniQH04az+H9PaeELoeIyCPYvXa5RCIZtkUAEQ3w9/HC/VfNwu/KDmPzFzrcsSwaIhEXrSMisofNLS9GoxGPPPII8vLyUFtb68SSiDzLPStmwd9HjL0nTNipaxW6HCIit2dzeJFIJHj55ZeRn5+PyspKPPLII0PWeCGikYUG+iIjKQYAsPlzncDVEBG5P5HFjvX8jUYjiouLodPpkJaWhrVr1zqytnEzmUyQSCQwGo2cxk0uoaGlEzf+7lNYLEDZz67DnMhgoUsiInI5tn5/2xVeLlRTU4OtW7dCJBIhLS1N0IG7DC/kinLerMJH+85gfXIMCtKXCl0OEZHLmfDwcqFt27ahrKwMs2fPxpNPPunoy18Wwwu5ouqjeqz7y074eonx5VM3Ymqwv9AlERG5FEHDyyCj0QiJROKsy4+K4YVc1Xdf2oGaYwb86MbZeHLNPKHLISJyKbZ+f9u1zsumTZtGPF5RUYEtW7YIElyIXFn2tXIAwFvfHEVXT5/A1RARuSe7wktCQgIaGhrQ0NBgPbZ582YUFRUhNDQUW7ZsGfe1dTodiouLoVKpUFhYCJ2OszTI/a1eFIWZsgAYunrxj12NQpdDROSW7AovpaWlUCgUSEtLs7bCqFQqFBQUYN26dQgNDR33tVUqFbKzs5Geno7c3FwUFBTYUyqRS/ASi5B93UDry5YvdOjp45YBRERjZVd4SUpKQl1dHY4cOQKpVAoA0Ov1iI2NBQC7VhLdunWrPaURuaz0pBnWLQPeqT0hdDlERG7HrvBiNBqtfzYYDACAtrY26zF7unpkMhmSkpKg0+mgVqu54SN5DH8fL/zgmoENG1/+TAuz2Wlj5omIPJLdGzPKZDKEhYWhtbUVGzduhEKhwKZNm1BbWwt7JjKVlpYCAOLj41FaWor09HR7SiVyKfesmIlgf29omzvx8f4zQpdDRORW7J4qPdj6cuHMovLycqjVauTn54/7uiqVClKpFDqdDjk5OcjOzkZRUdGI53Z3d6O7u9v63GQyISYmhlOlyaVt/Ogg/vyJFstipPjXYyu5YSMRTXoOnypdUVEx4nGJRDJsSnRqaqo1uIz2vkvR6XSorKyEQqFAdnY2tFotSkpKRu2Gys/Pt9YhkUgQExMz5s8kmmgPXh0HP28xdjcasFPLDRuJiGzlbeuJH3/8MeRy+ZgubrFYUFZWNuatAjQaDVJSUqzP5XI58vLyrONqLpaXl4cnnnjC+nyw5YXIlYUH+WF9Sgz+uvMo/vKZFitnhwtdEhGRW7A5vKxfvx7V1dVj/oDMzMwxvycxMRFFRUVDxrm0trYiMTFxxPP9/Pzg5+c35s8hElrWtXL87Ztj+OJIC/YcN2LJDC7sSER0OU7dHsAearUaGo3GOgVboVDY3PLD7QHInfxsay3+WXMCNy2Kwsv3JgldDhGRYFxibyOhMLyQOzl8ph2rn/8cAPDRT6/DvKhggSsiIhLGhOxtRET2mxsZjFuWRAEA/vRJncDVEBG5PoYXIhfwoxvnAADe+/Yk6po6BK6GiMi1MbwQuYCF0SFIWxgJiwX4M1tfiIguyeHh5cIdponIdo+vGmh9eaf2BBpaOgWuhojIddkdXmpra1FRUWF9KJVKR9RFNOksmSHBjfMiYGbrCxHRJdm8zstIMjMzYTAYrNOZAaCmpsbemogmrR+nzsEnh5qxveYEHk+dgxhZgNAlERG5HLvCS1paGrKysoYc27Ztm10FEU1miTNDce2ccHxxpAUvfapF/tolQpdERORy7Oo2io+Pt+kYEdnu8dSBsS+q6kacMJwVuBoiItczppaXLVu2DHne1taGoqIi6z5EFosFJSUlqKysdFyFRJNMSqwMV8nDsFPXij9VHEH+2qVCl0RE5FLG1PLy8ssvQ6vVoq2tDW1tbQCA5ORkWCwWDC7U64EL9hJNuCdWzwUAlFYdx7HWLoGrISJyLWNqeSkoKEBqauolz1EoFHYVREQDrS+DY1/+UH4Ev8tcJnRJREQuY0wtLxcHlwvXdDEajdi2bRtCQ0MdUhjRZPfz1fMAAP+sOQ5tM1fdJSIaZNeAXbVabf2zRCLBunXrhhwjovFbHiOFYsFUmC3AH9RHhC6HiMhljHmqtNFoRElJCUQiEcrKyoa9Xl1djYcfftghxRFNdj9Lmwv1gSa8++1J/PDG2dxxmogI42h5kUgkUCgUqKqqglarRV1d3ZBHbm6uM+okmpQWRUtw8+IoWCzAC+rDQpdDROQSRBY7pgeVl5dfdgCvEEwmEyQSCYxGI0JCQoQuh8guh8+0Y80Ln8NiAd778TVYPF0idElERE5h6/e3XWNeRtsKoKKiYtiaMEQ0PnMjg3H70mgAwO/L2PpCRGRXeElISEBDQ8OQWUebN29GUVERQkNDGWCIHOSnijnwEotQcbAJlQ16ocshIhKUXeGltLQUCoUCaWlp2LRpEwBApVKhoKAA69at47RpIgeRRwQhM3kGAKDgg4NcDJKIJjW7wktSUhLq6upw5MgR687Ser0esbGxAACRSGRvfUR03k9S58LPW4yqo20oP9AkdDlERIKxK7wYjUbrnw0GAwBYtw0AAJ1OZ8/liegCURJ/PHB1LABg40eH0G9m6wsRTU52hZe4uDjIZDKEhYWhtbUVGzduhEKhwKZNm1BbW8umbSIHe+z62Qjx98ahM+34V80JocshIhKEXVOlgf+2vkgk/52+WV5eDrVajfz8fPuqGydOlSZP9pdPtSj48CCmS6eg4snr4eftJXRJREQOMSFTpYGB0HJhcGloaEBqaqpgwYXI0z2wMhaRIX44YTiLv319TOhyiIgmnN3hpba2FhUVFdaHUql0RF1ENIopvl74SepcAMCfPqlD+7legSsiIppYY97b6EKZmZkwGAzWmUbA6AvXEZHjZCbPwJYvdNC1dKLoMx2eXDNP6JKIiCaMXeElLS0NWVlZQ45t27bNroKI6PK8vcTIvWkeHnlLgy1f6nD3ipmYJpkidFlERBPCrm6j+Ph4m44RkeOtWRSFK2JlONdrxqaPuG0AEU0edrW8aLVaFBUVISUlBQBgsVhQUlKCyspKhxRHRKMTiUT45a0LcOefd2B7zXE8eHUsN20koknBrpaXoqIixMXFwWKxWNd04douRBNneYwUdyyLhsUC/PY/B/jfHxFNCna1vBQUFCA1NXXIMYVCYVdBRDQ2v1gzDx/uO42dulZ8cqgJq+ZHCl0SEZFT2dXykpqaio0bN2L9+vUABhan45gXookVIwvAg+e3DXju/YPo6zcLWxARkZPZFV7y8vIglUqtrS2pqalQq9UOKYyIbPfDG2dDFuiLuqYO/KOyUehyiIicyq7wkpycjKysLMjlckfVQ0TjEOLvg58q5gAAni87DONZLlxHRJ7LrvBSX18PYGDWwyDONCISxl1XzER8RCBaO3vwYvkRocshInIau8JLQkICkpOTUVBQgLy8PKSkpCAtLc1RtRHRGPh4ifG/ty8CALz+VQPqmjoEroiIyDnsHrBbWlqKhIQEWCwWFBcXY9WqVY6qjYjG6Pq5EVAsmIo+swW/eW8/p04TkUcSWTzwXzdbt9Qm8kQNLZ1Ie/4z9PZb8OoDyZw6TURuw9bvb7t3lQYG9jPauHEjVq9ejZtuuskRlySicYoND8RD18QBAH7z3gH09HHqNBF5Foe2vBgMBiQnJ6Ours5RlxwXtrzQZNfR3YcbN32K5vZu5N08HznXc/0lInJ9E9ryMkgqlSI9Pd2RlySicQjy80bumnkAgBcr6tDUfk7gioiIHGdM4WX79u2XPWf27NnjLoaIHGdd4gwsmyFBR3cfNnxwUOhyiIgcZkzhpaysDO3t7TCZTKM+tFqts2olojEQi0V49juLIRIB2zUnsKteL3RJREQOMaYxL2KxeMiCdBezWCwQiUTo7+93SHHjxTEvRP+Vt30P3t51DPMig/He49fAx8uhvcVERA7jlDEv2dnZqKurg16vH/FRV1eHrKwsu4snIsfJXTMPoQE+OHSmHW981SB0OUREdhtTeMnJyUFcXBwkEsmID7lcjpycHGfVSkTjEBroi6dung9gYN+j00YO3iUi9zam8JKQkOCQc4hoYmUkxSBxphSdPf34v//sF7ocIiK7uHTnt1qtRnFxMdRqNdRqtdDlELktsViE39y5GGIR8N63p7CjrkXokoiIxs1lw4tarUZpaSmys7PZHUXkAIuiJbjvqlgAwK/e2YvuPmEH1hMRjZfLhpecnBwUFBQAAORyOcrKygSuiMj9PbF6LsKD/KBr7sRfPuWyBkTknlwyvOh0Ouj1ekilUmg0GhgMBsjlcqHLInJ7If4+eOaOhQCAlz7Roq6pXeCKiIjGzu7wsnHjRqxfvx4AUF5eDpPJZHdRGo0GMpkMKpUKcrkcxcXFUKlUo57f3d09bLE8IhrZrUum4cZ5EejpN+OX2/fCbPa4jeWJyMPZFV6eeuopSKVSKBQKAEBqaqpDBtbq9XrodDooFApIpVJkZ2cjIyNj1PPz8/OHTNmOiYmxuwYiTyUSDQzeneLjhV0NepRUNQpdEhHRmNgVXlJSUpCVleXwLh25XA6pVAqpVAoA1v/VaDQjnp+Xlwej0Wh9NDbyH2OiS5kRGoCfr54LAHju/QPcuJGI3Ipd4aW+vh4AhmwZUFlZaV9FwJjDkJ+fH0JCQoY8iOjSHlgZi8XTQ2A614dfv8u1X4jIfdgVXhISEpCcnIyCggLk5eUhJSUFaWlpdhcll8uRnJwMg8EAYGAAr1wuR2Jiot3XJqIB3l5ibFi71Lr2yycHm4QuiYjIJmPamHEk9fX1KCoqAgCsX7/eYSvsGgwGKJVKJCUlobq6Gkql0uYWGW7MSGS7/3tvP7Z8WY9pEn98/LPrEOzvI3RJRDRJ2fr9bXd4cUUML0S26+rpw00vfIFj+i7cdcVM5K9dInRJRDRJOWVX6Ut9GKcnE7mnAF9vFKxbCgB4e9cxbh1ARC7PrvBiNBqxevVqSKVShIaGYs2aNQwxRG7oqvgw3LNiJgBAue1bdHb3CVwREdHo7AovSqUSOTk5MJvN6O/vR1ZWFvLz8x1VGxFNoKduXoDp0ik43nYWGz86JHQ5RESjsiu8JCUlYd26ddbn6enpSE5OtrsoIpp4QX7e1vEur3/VgF31eoErIiIamV3hJSwsbNix0NBQ659ra2vtuTwRTbDr5kYgM3kGACBXtRtne7jzNBG5Hm973lxWVgadTmddAddgMECr1UKn0wEASktL8dFHH9ldJBFNnP+5dSE+O9yMhtYuFHx4EM/csUjokoiIhrCr5aWsrAwtLS2oq6tDXV0dWlpaIJFIrM/1ejY7E7kbyRQfFKYvAzDQfcTZR0TkauxqeSkqKkJqauqor5eXl9tzeSISyPVzI3D3lTPxt2+O4cnS3fjwp9dBMoWL1xGRa7Cr5eXi4FJRUYHt27eP+joRuY9f3rIAs8ICcMp4Ds++u0/ocoiIrOxqeQGA7du3W8e4WCwWVFVVYe3atXYXRkTCCvTzxu8yliGzaCe2a05g9cIo3LQ4SuiyiIjsa3l56qmn8PHHH2PXrl1oaWmBVqtFTk6Oo2ojIoElx8qQc308AOB//rkHLR3dAldERGRny0t8fDyysrJQX18PkUiE2NhYVFRUOKo2InIBP1XMwScHm3DwdDue2vYtNt+XDJFIJHRZRDSJ2dXyIpfLcfToUcTFxUGlUjmqJiJyIX7eXnh+/XL4eomhPtCEt745JnRJRDTJ2RVeDAYD5HI5TCYTWlpasGbNGhQVFTmqNiJyEQumhUB583wAwP+9tx+Hz7QLXBERTWYii8VicdTFysvLkZycDIlE4qhLjoutW2oTke3MZgseeL0Snx9uxvyoYPzrh1fD38dL6LKIyIPY+v1t967SmzZtGrKTNPvCiTyTWCzCpoylCAv0xcHT7Sj8kJs3EpEw7AovJSUlaGn57+qbqampUKvVdhdFRK5parA/CtOXAgBe3VGPTw81CVwREU1Gdm/MuGHDBnbNEE0iqQsicd9VswAAT5Z+i+Z2Tp8moollV3jZtWsX2tuHDtyrrKy0qyAicn2/vGUB5kYGoaWjG0+U1MJsdtjQOSKiy7IrvGRnZyMhIQFr1qzB+vXrMWfOHKSlpTmqNiJyUf4+XvjT9xPh7yPGF0da8NKndUKXRESTiN3rvFRXVyM9PR0pKSkoKyvDqlWrHFUbEbmwuZHB+PV3FgMAfl92GN/oWgWuiIgmC7u3BygtLUVmZiY+/vhjKJXKIRszEpFny0iagbUJ02G2AI//owat3D6AiCaAXeElJSUFDz/8MDZv3oykpCRs3boVra387YtoshCJRPjNnYsRHxGIM6ZuPFGym+NfiMjp7AovoaGhAICtW7di/fr1AACZTGZ/VUTkNgL9vPHnuxPh5y3GZ4eb8fLnWqFLIiIPZ1d40Wq1KC8vh1arxfLly1FfX4+2tjZH1UZEbmJ+VAieuWMRAGDTR4fwlbblMu8gIho/u8JLZmYmNBoNqqurYTKZUFxcDIPB4KDSiMidfC8lBmsTz49/ebsGp43nhC6JiDzUmMJLRUUFtmzZgi1btlj3H/jFL34BnU6HrVu3stWFaBITiUT47Z1LMD8qGC0dPXjsb9Xo6TMLXRYReaAxbcwYFhaG8vJyLF++fMTXDQYD4uPjBR+0y40ZiYRztLUTt734JdrP9eGBlbHW7iQiostxysaMWVlZ1uDS0NAw5AEAUqkUWVlZ4y6aiNzfrLBAvLB+OQDg9a8a8E7tCWELIiKPM6bwEhYWZv1zW1sbMjIyoFKphpwTHx/vmMqIyG2lLojEj1fNBgA8tW0PDp42XeYdRES2G1N4kUql1j8nJCQgMzMTTz75JGJjY63HRSKRo2ojIjf2U8VcXDsnHGd7+5H11yq0dfYIXRIReYgxhRedTof29naYTCaYTCaIRKIhz00mE7RarvFARICXWIQ/fi8BM2UBaNSfxY/e1qCvnwN4ich+YxqwKxaLh7SsWCyWEZ/39/c7tsox4oBdItdx6HQ7vvvSDnT19OOhq+Pwv7cvFLokInJRThmwm52djbq6Ouj1euj1erS1tVn/rNfrUVdXxwG7RDTEvKhg/D5zOQDg1R31KK1qFLYgInJ7Y2p5qampQUJCgt3nOBtbXohcz/Nlh/GH8iPw9RJja84KJMwMFbokInIxTml5sSWUCB1ciMg1/SR1DlYvjERPvxnZb1bjpOGs0CURkZuya3sAIiJbicUi/H79csyLDEZzezcefqMKnd19QpdFRG6I4YWIJkyQnzdeeSAZ4UG+2H/KhJ/8owb9Zpt7romIADC8ENEEmxEagOL7kuHrLYb6QBM2fHBA6JKIyM04PLwMbhVARDSaxJmh+F3GMgDA5i/q8fauYwJXRETuxNveC9TW1kKv11ufFxUVYevWrfZelog83O3LolHf0onflx3Gr/61FzGhAbhmTrjQZRGRG7ArvGRmZsJgMAzZNqCmpsbemohokvjxqtnQNXfgX7Un8chb1SjJuQoLo7m8ARFdml3hJS0tbdiidNu2bbOrICKaPEQiEQrSl+KU8Ry+qdfjwdd34Z+PXY1o6RShSyMiF2bXmJeRdpDmrtJENBZ+3l4ovi8ZcyODcMbUjQde2wXj2V6hyyIiF2ZXy4tWq0VRURFSUlIADOxtVFJSgsrKSocUR0STg2SKD1578AqsfWkHDp/pQM6bVXjjoSvg5+0ldGlE5ILsankpKipCXFwcLBYLBncZGMNuA0REVtOlU/DqAykI8vPG1zo9niz9FmauAUNEI7Cr5aWgoACpqalDjikUCrsKIqLJa1G0BH+5JxEPvlaJd3efRFigL56+feGQ3euJiOxqebk4uFRUVKC+vt6ugkaiVCphMBgcfl0icj3XzonApvNrwLz+VQNerKgTuCIicjVj2lV6JNu3b4dOpwMw0GVUVVXl0HVeNBoNkpKS0NbWNmRK9qVwV2ki9/fajno8++5+AMBv7lyMe1fMErgiInI2W7+/7eo2euqpp2AwGKDX6yGXy2EwGJCTk2PPJYfR6XSQy+UOvSYRub4Hr45DW2cP/lhRh/99Zy9CA3xw29JoocsiIhdgV3iJj49HVlYW6uvrIRKJEBsbi4qKCkfVBpVKhfT0dCiVSoddk4jcx8/S5kLf1YO3vj6Gn22tRbC/D66fGyF0WUQkMLvGvMjlchw9ehRxcXFQqVSOqgkAhq3ceynd3d0wmUxDHkTk/kQiEZ69YzFuWzoNvf0W5LxZha91rUKXRUQCsyu8GAwGyOVymEwmtLS0YM2aNSgqKnJIYSUlJTbPXMrPz4dEIrE+YmJiHFIDEQnPSyzC7zOXY9X8qTjXa8YPXq+E5lib0GURkYDsHrB7ofLyciQnJ0Mikdh1HbVajeTkZGvLS3x8PKqrq0dtienu7kZ3d7f1uclkQkxMDAfsEnmQc739+MEbldhR14oQf2/8PWsFFk+3798aInIttg7YtavlBQA2btyI9evXW587aj2GkpISFBcXo7i4GDqdDvn5+dBoNCOe6+fnh5CQkCEPIvIs/j5e2HxfMlJiQ2E614f7Xt2Fw2fahS6LiARgV8vLU089Zd3LaHCDxu3bt2Pt2rWOqe48kUgErVZr86wjTpUm8lzt53px95Zv8O1xIyKC/fCP7BWIjwgSuiwicoAJaXlJSUlBVlaW06YyGwwGFBYWAhhYzXe0lhcimjyC/X3w14euwPyoYDS3d+N7xV9D29whdFlENIHsCi+Dq+le2FXkyE0ZpVIpcnNzYbFYUFRUhMTERIddm4jclzTAF3/PWsEAQzRJ2RVeEhISkJycjIKCAuTl5SElJQVpaWmOqo2IaFSyQAYYosnK7r2NSktLkZCQAIvFguLiYqxatcpRtRERXdJIAaauiQGGyNM5dKq0q+CAXaLJRd/Zg+9v/hoHT7cjPMgXbz18JeZH8b99InfjlAG7mzZtuuw5W7ZsGcsliYjsNtgCs3BaCFo6evC94q/x7XGD0GURkZOMqeVFJpMhJSXlkudUVVWhtVXY5bvZ8kI0ORm7enH/a7tQ22hAsJ83XnswBcmxMqHLIiIbOWVX6dTUVISFhSEpKWnUczywF4qI3IQkwAdvPXwlfvB6Jb6p1+PeV3bhlfuTsXJ2uNClEZEDjSm8lJaWwmg0oqqqCsDAOi8XJyOZjL/lEJFwgvy88fqDVyD7zSp8caQFD7xeiT9/PxFpCyOFLo2IHMSuAbs1NTXQ6/UQiUQuNcuI3UZE1N3Xjx/9vQZl+8/ASyzChrVLkJHMTVuJXJmt398Om21UUVGBsrIypKWlCR5kGF6ICAD6+s14avseqKqPAwB+ect8ZF8XL3BVRDSaCduYsba2Fo8++ijS09NRVlYGnU5n7yWJiBzC20uMjelLkX3dwBYmz71/EBs+OMixeURubkxjXgY1NDSgtLQURUVFEIlEWLduHaqrqxEXF+fo+oiI7CISifDLWxZAFuiLDR8cxMufadHa0Y3n1i6Bj5fdv78RkQDG9F/uli1bkJKSgqSkJOh0OpSWluLIkSPYsGGDNbhs377dKYUSEdnjkevjUbhuKcQioLT6OB5+owod3X1Cl0VE4zCmMS9isRjp6elYv349pFLpkA0ZAaCtrQ0bNmxw6OaM48ExL0Q0GvX+M/jx2zU429uPRdEheO2BFEwN8Re6LCKCk9Z5yc7ORmFh4SX7i7du3TqWSxIRTSjFwki8nb0CP3i9EvtOmvDdl77CGw+lYPbUYKFLIyIbjanlpaamBgkJCXaf42xseSGiyzna2on7X92FhtYuSKb4oOjeJKyQhwldFtGk5pTZRraEEqGDCxGRLWaFBWLboyuRMFMK49le3PvKNyipahS6LCKyAYfaE9GkFRbkh7ezVuDWpdPQ229Brupb5H9wAGYzp1ITuTKGFyKa1Px9vPDi9xLw+KrZAICiz3R49G/V6OrhTCQiV8XwQkSTnlgswhOr5+GF9cvh6yXGR/vOIOPlnThhOCt0aUQ0AoYXIqLz7kyYjr9nXYmwQF/sO2nCHS9+iV31eqHLIqKLMLwQEV0gOVaGd350NRZOC0FrZw++v/lrvPX1UaHLIqILMLwQEV1kRmgAtj26ErctnYY+swX/7197kbd9D3r6zEKXRkRgeCEiGtEUXy+8eFcClDfNh0gEvL3rGNYX78QpI8fBEAmN4YWIaBQikQiP3hCPVx9IQYi/N2qOGXDbH7/EV3UtQpdGNKkxvBARXcaN86bivR9fax0Hc88r3+ClT+suuVUKETkPwwsRkQ1mhgVg+2MrkZE0A2YLUPjhIWS/WQ1jV6/QpRFNOgwvREQ28vfxQmH6UuSvXQJfLzHK9p/BLX/8AppjbUKXRjSpMLwQEY2BSCTCXVfMxPbHVmJWWABOGM4i8+WdKP5cy20FiCYIwwsR0Tgsni7Bez++Breen0793PsH8fBfq6Dv7BG6NCKPx/BCRDROwf4++NNdCfi/OxfD11uMioNNuOmFz/HlEc5GInImhhciIjuIRCLcs2IW/vnYSsRHBKKpvRv3vPINnnv/ABe1I3IShhciIgdYFC3Bez++Ft+/ciYAoPhzHb770g7UNXUIXBmR52F4ISJykCm+Xnjuu0tQdG8SQgN8sO+kCbe9+AXe+KqBg3mJHIjhhYjIwdYsisKHP70O18wOx7leM57+9z7c9+ounDRwawEiR2B4ISJygsgQf/z1oSvw7B2L4O8jxpd1LVjzwufYrjnOlXmJ7MTwQkTkJGKxCPevjMX7j1+L5TFStJ/rwxMlu5HzZjWaTOeELo/IbTG8EBE5mTwiCKpHrsKTq+fCx0uEj/efQdrzn2NbNVthiMaD4YWIaAJ4e4nxo1Vz8O8fXYMl0yUwnu3Fz0t348HXKzkWhmiMGF6IiCbQgmkh+OdjK5F70zz4eovx6aFmrH7+c/x1ZwP6OSOJyCYML0REE8zbS4zHbpiN9x+/Bokzpejo7sP/vrMP6S9/hYOnTUKXR+TyGF6IiAQye2owSh9ZiV9/ZxGC/LxRc8yA2/74JQo+PIhzvf1Cl0fkshheiIgE5CUW4b6rYqF+4nqsWRSJPrMFf/lUi7TnP0PFwTNCl0fkkhheiIhcQJTEH0X3JqPo3iREhfijUX8WD71ehay/VqFR3yV0eUQuheGFiMiFrFkUhfKfX4+c6+TwFotQtv8M0p7/DH+qOMKuJKLzRBYPXGTAZDJBIpHAaDQiJCRE6HKIiMbl8Jl2/Opfe/FNvR4AMFMWgP936wKkLYyESCQSuDoix7P1+5vhhYjIhVksFrxTexLPvX8ATe3dAIBr54TjV7ctxNzIYIGrI3Istw8vGo0GarUaAFBZWYnNmzdDKpXa9F6GFyLyNJ3dffjzJ3XY8kU9evrN8BKLcM+VM/ETxVzIAn2FLo/IIWz9/nbZMS9qtRq5ubnIzc1FSkoKUlNThS6JiEgwgX7eyL1pPtRPXI/VCyPRb7bgjZ1HcX3hJ3j5My3Hw9Ck4pItLxqNBqmpqWhrawMA6HQ6xMfHQ6vVQi6XX/b9bHkhIk+3o64Fv/3PAew/NbCo3XTpFOTeNA+3L42GWMzxMOSe3LrlJTExEZs3b7Y+NxgMAACZTCZQRUREruXq2eF478fXYFPGMkSF+OOE4Sx+8o9a3PHnL/H54WZu+EgezSVbXi6mVCqh0WhQVlY24uvd3d3o7u62PjeZTIiJiWHLCxFNCmd7+vHKlzq8/JkOHd19AICV8WFQ3jQfy2KkwhZHNAZuP2B3kMFgQFJSEqqrq0cdsPvMM8/g2WefHXac4YWIJpPWjm689KkWb+48ip5+MwDgpkVReGL1XM5MIrfgMeElJycHSqXykmNd2PJCRPRfx9u68IL6CLZrjsNsAUQi4Pal0fiJYg7iI4KELo9oVB4RXgoLC5Geng65XG4d92LLdGkO2CUiGljk7gX1Yby/5zQAQCwC7kyYjsdXzUFseKDA1REN59YDdgFApVIhMTHRGlxKSkpsXueFiIiAuZHBeOnuJPzn8WugWDAVZguwXXMCq373KX62tRZ1TR1Cl0g0Li7Z8jI4NfpCUqnUOnX6ctjyQkQ03O5GA55XH8anh5oBDHQn3bpkGn68ag7mRXFMDAnPI7qNxovhhYhodN8eN+DFijqU7T9jPaZYEIlHb4hH0qxQASujyY7hheGFiOiS9p804U+fHMEHe09j8JvgyjgZHr0hHtfPjeDmjzThGF4YXoiIbKJt7kDRZ1r8s+YEevsHvhLmRwUj61o5bl8WDV9vlx0eSR6G4YXhhYhoTE4Zz+KVL+rx913H0NUzsFdSZIgfHlgZh+9fOROSKT4CV0iejuGF4YWIaFyMXb34266jeH1HA5raB9bQCvD1QnrSDNy/MpZrxZDTMLwwvBAR2aW7rx/v7j6FzZ/rcOhMu/X4DfMi8MDKWFw3J4KbQJJDMbwwvBAROYTFYsFX2la8tqMe5QebrIN75eGBuHvFLKQnzoAkgF1KZD+GF4YXIiKHa2jpxF93HkVpVSPaz28C6e8jxp3Lp+OeFbOweLpE4ArJnTG8MLwQETlNZ3cf/lV7Am/uPIqDp//bpbR0hgR3XTETty+LRpCft4AVkjtieGF4ISJyOovFgqqjbXhz51F8sPeUdap1oK8X7lgeje+lzMTSGRKuGUM2YXhheCEimlCtHd3YrjmBt3cdg66l03p8XmQwMpJn4LsJ0xEW5CdgheTqGF4YXoiIBGGxWPBNvR7/2HUMH+w9je4+MwDAWyxC6oKpWJc4AzfMm8rF72gYhheGFyIiwRnP9uLd3SdRWtWI3ceN1uOhAT64fVk01ibOwDJ2K9F5DC8ML0RELuXQ6XaUVjXind0n0Xx+8TtgYMr17cui8Z3l0ZBzAbxJjeGF4YWIyCX19ZuxQ9uK7Zrj+GjfaZzrNVtfWzJdgu8sj8YtS6YhWjpFwCpJCAwvDC9ERC6vo7sPH+87jXdqT+LLuhb0m//7lZQ8KxS3Lp2GW5ZMQ2SIv4BV0kRheGF4ISJyK60d3Xh/zym8U3sSVUfbrMdFIiBllgw3LY7CmsVRmM4WGY/F8MLwQkTktk4bz+H9Pafwnz2nUH1BkAGAZTMkuGnxNKxeFMlNIj0MwwvDCxGRRzhpOIsP957Gh/tOo7JBjwu/teIjApG2MAppCyORECPlRpFujuGF4YWIyOM0t3fj4/2n8eHe0/ha12pd0RcAwoP8cOO8CKQumIpr5kRwewI3xPDC8EJE5NFM53rx2aFmfLz/DD492GTdKBIAfLxEWCEPww3zpuKGeRGQhwdyLRk3wPDC8EJENGn09JlR1aBH+cEmlB84g4bWriGvx8im4Ia5A0FmhTwMgWyVcUkMLwwvRESTlq65A+UHmvDp4SZU1rehp/+/a8n4eImQODMU182NwLVzwrE4WsKxMi6C4YXhhYiIAHR29+ErbSs+PdSEzw4343jb2SGvSwN8cJU8DCtnh2NlfBi7mATE8MLwQkREF7FYLDja2oUv6lrwxeFmfKVtRccFY2UAYJrEHyvkYbhKHoYV8jDEyKYwzEwQhheGFyIiuoy+fjN2Hzfiq7oW7NC2QHPUMKSLCQCiJf64Uh6GK+NkSImTsWXGiRheGF6IiGiMzvb0o/poG77WteJrXSt2HzcMmY4NAOFBvrgiToakWTIkzwrFwugQ+HiJBarYszC8MLwQEZGdunr6UH20Dbvq9dhVr0dNowE9fUNbZqb4eGFZjATJs2RInCXF8phQyAJ9BarYvTG8MLwQEZGDdff149vjRuyq16P6aBuqj7bBeLZ32Hlx4YFImClFQsxAmJkXFQxfb7bOXA7DC8MLERE5mdlsgba5A1Xng4zmWBt0zZ3DzvP1FmNRdAiWzZBiWYwES6ZLIQ8P5BTtizC8MLwQEZEADF09qGk0oOZoG2qPG7G70TBi60yQnzcWRYdgyXQJFk+XYPH0EMSFB8FrEgcahheGFyIicgGD07NrGw3YfdyAPceN2HvSiHO95mHnTvHxwoJpwVgULcHC6BAsnBaCeVHB8PfxEqDyicfwwvBCREQuqq/fDG1zJ3YfN2D/SRP2njBi30kTzvb2DztXLBoYQzN/WggWRAVjflQI5k8LxnSp560/w/DC8EJERG6k32xBfUsn9p00Yv8pE/afHHi0dvaMeH6QnzfmRgZhXlQw5kYGY15kMGZHBiEiyM9tQw3DC8MLERG5OYvFgqb2bhw4ZcLB0+04eP5/65o60Gce+etbMsUHc6YGYU5kEOIjgjB76sAjWjLF5QcIM7wwvBARkYfq6TOjvqUTh8604/Dp9vOBph1H9V0Y7Vt9io8X4sIDIY8IhDwiCPERgYgLH3gE+/tM7A8wCoYXhhciIppkzvX2Q9fciSNN7ThypgPa5g7UNXWgobVz2ErBFwoP8kNceABiwwIRGx6I2LBAzAoLwKywgAkNNgwvDC9EREQAgN5+M47pu6Br7kR9Swd0zZ0Dj5YOtHSMPKZmUHiQL2bKAqyPGFkAZoUFYl5kMCQBjg02tn5/ezv0U4mIiMjl+HiJER8xMAYGiBzymulcL462dEHX0oGGli40tHbiaGsnjrZ2obWzBy0dAw/NMcOQ9xWmL0VmcszE/RAXYHghIiKaxEL8fbBkhgRLZkiGvdZ+rhdHW7twTH/B4/zz2LBAAaodwPBCREREIwr29zm/+u/wYCMk7hJFREREboXhhYiIiNwKwwsRERG5FYYXIiIicisML0RERORWGF6IiIjIrTC8EBERkVtx2XVedDodVCoV5HI5dDodsrOzIZVKhS6LiIiIBOay4SUjIwPV1dUABoJMVlYWSktLBa6KiIiIhOaS3UY6nW7Ic7lcDrVaLVA1RERE5EpcMryo1WrIZLIhx2QyGTQajUAVERERkatwyW4jg8Ew4nG9Xj/i8e7ubnR3d1ufm0wmZ5RFRERELsAlW15GM1qoyc/Ph0QisT5iYoTZopuIiIiczyXDi1QqHdbKotfrR51tlJeXB6PRaH00NjZOQJVEREQkBJfsNlIoFCgqKhp2PDk5ecTz/fz84OfnZ31usVgAsPuIiIjInQx+bw9+j4/GJcOLXC4f8lyn0yE5OdnmdV7a29sBgN1HREREbqi9vR0SiWTU10WWy8Ubgeh0OhQVFSElJQWVlZXIy8uzObyYzWacPHkSwcHBEIlEI54zeN3RjPS6yWRCTEwMGhsbERISYvPP4gou9/O64ueM91pjfZ+t59ty3njuK8B9763JdF+N572Ourd4X7n+Z7njfXW5c4S4rywWC9rb2xEdHQ2xePSRLS7Z8gIMtL4UFBQAANLT08f0XrFYjBkzZlzyHC8vr0v+n36p10NCQtzqHwLg8j+vK37OeK811vfZer4t59lzXwHud29NpvtqPO911L3F+8r1P8sd76vLnSPUfXWpFpdBLjlgdyL88Ic/tOt1dzNRP48jP2e81xrr+2w935bzeF+5/ufYcy2h7i3eV67/We54X13uHFe+r1y228gVmUwmSCQSGI1Gt/othlwf7y1yBt5X5AyucF9N2paX8fDz88PTTz89ZGYTkSPw3iJn4H1FzuAK9xVbXoiIiMituOyAXXdSXFwMuVwOjUaD9PT0YVO9icbDYDAgPz8f69evR2JiotDlkIdQqVQAgMrKSqSlpUGhUAhcEXkKlUoFqVSKsrIy5OTkOPW7kN1GdtLpdNBqtVAoFMjNzYVSqRS6JPIQVVVVo26JQTQearUaOp0O6enpyMnJ4b9X5DAGgwGVlZVQKBRISUmxzhZ2FoaXi2g0GiQlJQ07rtPpUFhYCJVKhcLCQuuXilqtRnx8/JDziC421vsKGFhp2ta1jWhyGut9NfhL1uA5o61aTjTWe0sqlVoDy2DLizOx2+gCKpXK2v1zsYyMDFRXVwMY+MvLyspCaWkpDAbDkC8Y/qZMFxvPfUV0OfbeV0VFRU7/7Zjckz33llqthlQqdfovXgwvFxhtMbyLW1PkcjnUajWAgbTJwEKXMp77iuhy7LmvCgsLkZeXx/F5NCJ77i2FQgGZTIacnByUlZU5rUZ2G9lArVZDJpMNOSaTyaDRaJCcnIzW1lbrcQ6sJFtd6r4iGq/L3VdqtRoKhQKJiYnWwbtEtrjUvVVcXIzCwkIAA7/UO3sIBVtebDBay4per4dCoUBVVZV1INzmzZsntjhyW5e6r4CBfyguDDIMxmSLS91XOp0OGRkZkMvlMBgMUCgUY95+hSavS91bmZmZUKvVUKvVKCsrc3r3N8OLHQb/IrOzs4UthDzKhYMrOY2VHMVgMEAul6OtrU3oUsjDDI79HAzCE/HvFruNbCCVSq2/DQ/S6/WcCUJ24X1FzsD7ipzFle4thhcbjJYiOc2Q7MH7ipyB9xU5iyvdWwwvo7iwb+/iEfmD6yPwNxkaK95X5Ay8r8hZXPXe4piXCwwONAKA/Px8pKSkWPvwSktLoVQqkZKSgsrKSq7FQTbjfUXOwPuKnMUd7i1uzEhERERuhd1GRERE5FYYXoiIiMitMLwQERGRW2F4ISIiIrfC8EJERERuheGFiIiI3ArDCxFNChqNBkqlEmq12inXLy4uhlKpHHXzOiJyHIYXIhpCo9EgJycHIpEISqUShYWFKCwsRE5OzrBt7pOSkqBSqWy67qXOLS4uRmho6JBdtMdybVvodDrk5OQMW+J8MHQUFxdDpVJBrVajuLjY+rOqVCokJSVBJBKhsLBwyHsLCwsRGhqKnJwcZGdnIyUlZdjeL0TkBBYioou0tbVZAFja2tqsx8rKyixSqXTYsQufX8rlzlUoFJbq6upLnm/rZ42ktLTUotVqh31mQUHBkGPV1dUWAENqGTw20udf+P6RPoOIHI8tL0RkE4VCAYPBMKTbRaFQ2LyvyVjOHel8nU6HkpISm99/OYOtKLm5uUOOJyYmIjs7e9gxuVyO4uLiIcfVarV12XQimjgML0Rkk8EuncTEREE+v6CgwKHXy8/PR05OzoivZWRkDAtaOTk5KCoqGnJMo9EM26yOiJyPGzMS0ajUajWkUik0Gg1aW1uh1WqtX9YajQZZWVnW8R7AwPgQnU4HqVSK6upqZGRkQKPRQKFQDDtXo9Fg69atSElJAYAhY0UuvrZarUZVVZX1nJ6eHjz//PNITExEaWkpDAYDkpKSkJ6eblPI0el0MBgMowaPi8fFAEB2djaUSiV0Op31fdypmUgYDC9ENKrBrhuZTIasrCzk5eVZX0tMTMT69eutzw0GA7KystDW1gYAiI+Ph1KptAaBi8/NyMiAVqu1HsvPzx/12gqFAgqFAvHx8dbw4+vri+rqagADIUKpVA7r7nEkqVQKhUKBoqIiFBQUoLi4GJmZmU77PCIaHbuNiOiyEhMTkZycDKVS6ZDrlZSUDOt+kslkY7pGdna2dQyMTqdDcnKyze8dbDm5ePaUSqWCUqmESCRCTk7OsGnPOTk51nEvBoOBLS9EAmF4ISKbSKXSS66RIpVKkZ2djcLCQhQXFyMnJ8dp40EGQ0V2djaKi4uh0WjGPBYnNzd32BiWC7udcnJyhoWT9PR0GAwGFBcXc6wLkYAYXojIJvHx8daWigvXY7mwdSIsLAy5ubnIzs4eNovnwnMVCsWQawDDW0EuvvaFBkNUTk7OuAfyFhQUQK/XD5tBNFIdF0pPT4dSqeQsIyIBccwLEQ0xOJA2PT0dxcXFUCgU1unDZWVlKCwstLZIbN26FTKZDOnp6ZDL5dBqtYiPj7eOk8nIyEB2drb1mheeW1paCqVSibS0NGsXTH5+PgoKCmAwGIadPxhUBmsCBrp/EhMTxx0kqqurUVhYCKVSifj4eOsYnIKCglG7hPLy8tjqQiQwkcVisQhdBBG5P7VaDY1GY21x0el0UCqVI65q60gqlcqm8KJSqazrtTizFmd/BhGx24iIHKSsrGxISJHL5Vi/fv1lu2HGIycnxxqWhFp3hoiEw24jInKIgoICFBYWQq1WW1se9Hq9U6YvZ2RkWMfDMLwQTT7sNiKiSWFw3E1aWppTurGKi4uh1WqRl5fHKdRETsbwQkRERG6FY16IiIjIrTC8EBERkVtheCEiIiK3wvBCREREboXhhYiIiNwKwwsRERG5FYYXIiIicisML0RERORWGF6IiIjIrfx/q3pA/8AznbkAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "R_line = np.geomspace(1, 1e3, 100)\n", "\n", "fig = plt.figure()\n", "ax = fig.subplots()\n", "\n", "ax.plot(R_line, grammage(R_line))\n", "ax.set_xlabel(\"Rigidity [GV]\")\n", "ax.set_ylabel(\"Mean Escape Length [g/cm$^2$]\")\n", "ax.set_xscale(\"log\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "6", "metadata": {}, "source": [ "## Simulation settings\n", "\n", "Here we define the initial conditions and physical parameters:\n", "* Particles $10^5$ protons.\n", "* Energy Range 1 GeV to 1 TeV.\n", "* Medium: Interstellar hydrogen gas with number density $n = 1 \\text{ cm}^{-3}$." ] }, { "cell_type": "code", "execution_count": 4, "id": "7", "metadata": {}, "outputs": [], "source": [ "T_min = 1 * Units.GeV\n", "T_max = 1 * Units.TeV\n", "\n", "n_events = 100_000\n", "\n", "pdg_code = 2212 # proton\n", "n_ism = 1 # interstellar hydrogen gas density [cm^-3]\n", "m_H = 1.66e-24 # hydrogen mass [g]\n", "rho_ism = n_ism * m_H # mass density [g cm^-3]" ] }, { "cell_type": "markdown", "id": "8", "metadata": {}, "source": [ "We generate a sample of primary protons with energies distributed log-uniformly. For each proton, we calculate its rigidity $R$ and assign a corresponding grammage $X(R)$ that it attempts to traverse." ] }, { "cell_type": "code", "execution_count": 5, "id": "9", "metadata": {}, "outputs": [], "source": [ "def log_uniform(a, b, size=None):\n", " log_samples = np.random.uniform(np.log(a), np.log(b), size)\n", " return np.exp(log_samples)" ] }, { "cell_type": "code", "execution_count": 6, "id": "10", "metadata": {}, "outputs": [], "source": [ "T_range = log_uniform(T_min, T_max, n_events)\n", "R_range = np.sqrt((T_range / Units.GeV + 0.93827) ** 2 - 0.93827 ** 2) # rigidity [GV]\n", "X_range = grammage(R_range)" ] }, { "cell_type": "markdown", "id": "11", "metadata": {}, "source": [ "## Running the simulation\n", "\n", "We use `ProcessPoolExecutor` to parallelize the simulation.\n", "The dataset is split into `n_tasks` chunks. Each task simulates the transport of a batch of protons through the matter layer.\n", "\n", "**Note:** The `process_chunk_safe` function runs the interaction module for each particle. It returns a boolean flag indicating whether an interaction occurred (`True`) or the particle passed through without interaction (`False`)." ] }, { "cell_type": "code", "execution_count": 7, "id": "12", "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "100%|███████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 20/20 [00:33<00:00, 1.67s/it]\n" ] } ], "source": [ "from concurrent.futures import ProcessPoolExecutor, as_completed\n", "from tqdm import tqdm\n", "\n", "n_tasks = 20\n", "T_chunks = np.array_split(T_range, n_tasks)\n", "X_chunks = np.array_split(X_range, n_tasks)\n", "seeds = range(n_tasks)\n", "\n", "def process_chunk_safe(T_chunk, X_chunk, i_task):\n", " from gtsimulation.Interaction import NuclearInteraction\n", " sim = NuclearInteraction(seed=i_task)\n", " interaction_flag = np.zeros(T_chunk.size, dtype=np.bool)\n", " for i in range(T_chunk.size):\n", " primary, secondary = sim.run_matter_layer(\n", " pdg=pdg_code,\n", " energy=T_chunk[i],\n", " mass=X_chunk[i],\n", " density=rho_ism,\n", " element_name=['H'],\n", " element_abundance=[1.0]\n", " )\n", " interaction_flag[i] = (len(secondary) > 0)\n", " del sim\n", " return (i_task, interaction_flag)\n", "\n", "with ProcessPoolExecutor() as executor:\n", " futures = [\n", " executor.submit(process_chunk_safe, t, x, s) \n", " for t, x, s in zip(T_chunks, X_chunks, seeds)\n", " ]\n", " results = []\n", " for f in tqdm(as_completed(futures), total=len(futures)):\n", " results.append(f.result())" ] }, { "cell_type": "markdown", "id": "13", "metadata": {}, "source": [ "The results arrive out of order due to parallel processing. We sort them by task index to ensure they match the original energy array `T_range`." ] }, { "cell_type": "code", "execution_count": 8, "id": "14", "metadata": {}, "outputs": [], "source": [ "interaction_flag = np.hstack([item[1] for item in sorted(results, key=lambda x: x[0])])" ] }, { "cell_type": "markdown", "id": "15", "metadata": {}, "source": [ "## Analysis and Visualization\n", "\n", "We compute the interaction probability as the ratio of interacting particles to the total number of particles in logarithmic energy bins." ] }, { "cell_type": "code", "execution_count": 9, "id": "16", "metadata": {}, "outputs": [], "source": [ "bins = np.geomspace(T_min, T_max, 50) / Units.GeV\n", "\n", "interaction_count, _ = np.histogram(T_range / Units.GeV, bins=bins, weights=interaction_flag.astype('float'))\n", "total_count, _ = np.histogram(T_range / Units.GeV, bins=bins)\n", "\n", "interaction_fraction = np.divide(\n", " interaction_count,\n", " total_count,\n", " out=np.zeros_like(interaction_count),\n", " where=total_count != 0\n", ")" ] }, { "cell_type": "markdown", "id": "17", "metadata": {}, "source": [ "Finally, we plot the simulated data against the theoretical prediction for a constant nuclear interaction length $\\Lambda = 55 \\text{ g/cm}^2$." ] }, { "cell_type": "code", "execution_count": 10, "id": "18", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj8AAAGvCAYAAAC0IrTpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8ekN5oAAAACXBIWXMAAA9hAAAPYQGoP6dpAABZEUlEQVR4nO3de1ib9f3/8WegQOkBAmi1rdUS6rlqy8GzVttgp3MHFex0J+cE5pz7dh7I2Nn9NinMbW46V6i66TY3Sqw76DYlrXaeC8Sq9VzSam3V2kKgrS2lkN8ft0RSTklIcgN5Pa4rV5P7vnPfb7p78urn/hwsPp/Ph4iIiEicSDC7ABEREZFYUvgRERGRuKLwIyIiInFF4UdERETiisKPiIiIxBWFHxEREYkrCj8iIiISVxR+REREJK5MMLuA0aanp4dt27YxdepULBaL2eWIiIhIEHw+H7t27WLGjBkkJAzdtqPwc5Bt27Yxa9Yss8sQERGRMGzZsoUjjjhiyGMUfg4ydepUwPjLS0tLM7kaERERCUZHRwezZs3y/x4fisLPQXofdaWlpSn8iIiIjDHBdFlRh2cRERGJKwo/IiIiElcUfkRERCSuqM+PiIiMWt3d3XR1dZldhowCSUlJJCYmRuRcCj8iIjLq+Hw+3n//fbxer9mlyChitVo5/PDDRzwPn8KPiIiMOr3BZ9q0aUyaNEmTzsY5n8/HRx99xPbt2wGYPn36iM6n8CMiIqNKd3e3P/hkZWWZXY6MEqmpqQBs376dadOmjegRmCnhx+Px4HQ6sdlseDweSktLsVqtAx7rdrtxuVwANDY2smLFCv+xQ50nlGuIiMjo0dvHZ9KkSSZXIqNN7z3R1dU19sJPcXExzc3NgBFSSkpKqK+vH/BYl8tFeXk5ANXV1SxatMj/3aHOE8o1RERk9NGjrvjjdDoBo7GjsLAQu90esD9S90TMh7p7PJ6Azzabzd+yczC3201lZaX/c1FREW63G4/HM+R5QrmGiIiIhC7SndFdLhcej4eioiLKyspwOBwRPX9fMQ8/LpeLzMzMgG2ZmZm43e5+x+bm5rJixQr/596/6MzMzCHPE8o1REREosHhcPT7x3g0uN1u/+83j8cT8LtuqH0jFelwYrfb/U96PB4P+fn5ET1/XzEPP4MlxdbW1gG3FxUV+d/X1dVht9uxWq1DnieUa3R2dtLR0RHwEhERGana2lr/Y5xoqqmpIS8vD4vFQllZGTabLah9IxXNfrQ1NTVRbfkZNaO9hms+83q9OJ1Ofz+ecM4z0L7KykpuueWWICoUEREJjtvtJj8/n7q6On9rRrTk5eXR1tYG9A8kQ+0bidraWpYsWRKx8/VVXV1NRUVFRIPawWIefqxWa78WmNbW1mH/R3E4HDQ0NPiPG+o8oVyjoqKCG264wf+5o6ODWbNmBf8DSeT0dMOeHdDZYbz2dUDnLjiwz9jn6/74zx6YkPLxKxWSJkL6LDjkaLN/AhERwOjiUVNTQ05ODh6PJ6q/yGHoYBONFprm5mZKS0sjfl6Xy4Xdbic3Nxen0xnw9CeSYh5+7HY7NTU1/bYP9Wyvuroah8OBzWbzt94MdR6bzRb0NVJSUkhJSQnhJ5Cw9fRAx7uw4y3Y2QI7N0L7u7BrG+x6H3Z/YASbMO06YgEfzr+ejw4/1b8tY3IyM62pkaheRMzk80HXR7G/btIkCHOEkc1mw2az4XQ6o9r60/tkBIxRUn0fbw21byBOpxOPx4PVaqW5uZni4mLcbndA/V6vl5ycnIDveTweampqKCgooLW1lcsvv5ympiYcDof/mh6Ph4aGBqqqqvB6vbhcLlpaWvy/rz0eD8XFxf7f9Xa7ffyEn4P/0ns7NfUmU7fbjdVq9R/ndDrJzc31/2WsXLlywDl7+p5nqH0SIwc6YfursG09vPei8dr+qtGKMyQLpKTBxDRImWq8T5oIlkRISARLInsP+GhseZ9kXycpdDGRTo62bGXqu2uZ+u5anus5nt8euIRnek4kNWkCrhsXKACJjHVdH8GtM2J/3e9tg+TJIX3F7Xb7h2gXFRUN++jL6/UGjGweSFZW1qDn6Ps70WazUVhYSEtLy7D7BqqjpKTE/5gsJycHh8PRb7h5bW1tQKuP1+ulsLCQ5uZmrFYrDoeD2tpaysvLsdvtNDQ0+Keaqa+v94fB3NxccnJy8Hq9/t/7vdeONlP6/NTX1+NwOCgoKKCxsTFg/p3KykoKCgooLy/3p8C+rFar/y99qPMMtU+ioHMXbFkHbz9jvLY2Qff+/sclJEGmDbLmQFYOZBwFU2fA1MNh6nSYMs0IOUNo2drOV159ituXzGPOtCl0Ay0db3Poi7/H+mY9p/Mapye/xsacr2J/ZTFte/Yr/IhIzDQ1Nfl/Ty1ZsoTq6mr/L/iBWK1Wqqqqwr6ex+MhNzcXwN/C0vuobah94WppaQn4WVauXInNZvNvq6io8O/LysoKmKW7b+NG7+dgur5Eminhx2az+f+HPrhJq29Isdls+Hy+sM4z1D6JAJ8P3n8Z3noM3mqAdxuNPjl9pWbA9Hkw/RTjdfjJkDEbEiNz282ZNoW5M9ONDzNPhuN/D+0/gKdvh3W12DY9wKGcFpFriYjJkiYZrTBmXDdEfQfX5ObmYrVa/U8tBjs+3JYft9vNokWL+rWY9E7vMti+gfQ2LlRXV2O1Wgd8ROZ2uyksLOxXf9/wMlyQGQ1PYUbNaC8ZA7q7wLMWXv27EXh2vx+433okHHUWHHWm8WemLexn5WFLnwkX/QK2rSfh3XUsSXwcuDi2NYhI5FksIT9+MsNArSqXX3459fX1g4afkbT89P2HPhgdhouKivwtLIPtG8xQj9fAmHLm4FqLior69bPt7bg8kEhPjhgOhR8ZWm/geeUheP1h2Of9ZF/SZLAtgKMLYY7dCD+jxakl8O46rpywmtaeA2ZXIyJxwOVy+VtLevuq9q5K4HK5qK6ujnjHZ6vVSn5+vr+1pqWlxf8EZah9g2lpaSEnJwer1UpmZibFxcXDjurqHWTU29UEjEFJbreburo6/+e+kyzm5ub6Z3Suqqqiqqoqpi1CFt9Qz5XiUEdHB+np6bS3t5OWlmZ2OeZ5fwOs/wu8tBI+2vHJ9snT4ITPwnGfNlp3JsR+pNyGre1cfMdTPHz92Z889jrYgU4O3HY8E/bt5B17DUee/YXYFikiYdu3bx+bNm0iOzubiRMnml1O3HC5XAEjuzwej3+0lt1ux+l0YrVaB23RiYWh7o1Qfn+r5Uc+sa/dCDsv/MkYndVr8qFwwufghM8bj7SG6ZA8KkxIoe24L3Do+t+R+ep9oPAjIjKkhoaGgIkLbTYbS5Ys8S/R0dDQMOA0MmORwo/A9tdg3Qp48W/QtcfYlpAEx34K5n3JeKQVoU7KsdR6/JfIfOEupmx7Bj58Aw491uySRERGraqqKqqrq3G5XP5+S62trZSWlg44t89YNvZ+o0lk9PTAm/+B534Pm5/8ZPuhx0HeVXDS5TA5a9CvjwVdU2ayuieXCxKbofEeuKja7JJEREa1wfokWa3WqC/TEUsKP/HmQCe8VAdP/xZ2vmVssyQafXhOLYXZZ8d+hFYU3d99gRF+XvwrLPoRpEwxuyQRETGZwk+86NwFTfcaLT273jO2TUyH/K9Dwdch/Qhz64uSp3tOpDPdRkq7xwh9BV83uyQRETGZws94t38PNN4NT/8GPtppbJs6A874pvF4K2WqqeVFm48EWk/4MtOfvcX4e8i/ely1bImISOgUfsar/R8ZLT1P3w57PjS2ZdrgnBuN/jwTkk0tL5baji5ietMvjLXF3n4GZp9ldkkiEgTNxCIHi9Q9ofAz3nQfMObnefzWT2ZgzpgNCxxG6DF51NZW717a9gyw5lcfkV6JvSclHU6+HJr/COtqFH5ERrmkpCQAPvroI1JTtS6ffOKjjz4CPrlHwqXwM174fMY6Ww0/hg9fM7alHwkLboZTroDEkd0okbDVuxf7L9eyt6t7yONSkxIjvxL7qaVG+HntYWh/d9z2cRIZDxITE7FarWzfvh2ASZMmYdHj6rjm8/n46KOP2L59O1arlcTEkc03p/AzHrz3Ijz6/U+GrKdmwLnlRudeE2ZgHkzbnv3s7er2r8Y+kI3bd7O0bn3kV2I/7ESYfY7xd9R0rzHyS0RGrcMPPxzAH4BEwBhy33tvjITCz1i2Zyes+Sk03wf4IDEFTv8GnH0DpFrNrm5QAauxx9KppUb4af6jEQ6TNG2+yGhlsViYPn0606ZNo6ury+xyZBRISkoacYtPL4Wfsaj7gNF68fjPjCUpAOYWgf0nYJ1lammj2rEXQdoR0PEuvLIK5l1pdkUiMozExMSI/cIT6ZVgdgESoi3roOZc+M/NRvA57CT42n+g6B4Fn+EkTvhknp/na4x+UiIiEncUfsaKvV54+DtwTyFsf8Xo1/PpX0HZWmOxUQlO7leNx4PvrYd3G82uRkRETKDwM9r5fLDhQbizwHjUBTD/S3C922jFGAsrrI8mk7PgpGLj/fPjY3ViEREJjcLPaNa+FR5YAs6rYc92yDoarnoEPvc7mJRpdnVj12mlxp+v/h12vW9qKSIiEnsKP6ORzwfu++Gu0+GtRyExGc77Hlz7tLHwqIzM9FNg1unQcwCa/mB2NSIiEmMKP6ONdwv8+VL45/XQ2QEz86DsSTjPMarm7Bnzelt/mu41VroXEZG4ofAzWvh88MJf4K4zoGWN0Sm38Kdw9WMw7Tizqxt/jv+sscDrnu3wykNmVyMiIjGk8DMa7NkBdV+Cf3wT9u+CI06FbzwFZ/2f6WtxjVuJSZ8Me3/u9xr2LiISR/Sb1WxvPgb/uM5ogUhIgvO/Z4SeOB/FtXH77rD2hSTva/C/XxjD3rc8D0eeHpnziojIqKbwY5auvcZ6XE33GJ8PORYuW2F0xh2FYrUae8bkZFKTEllat37I41KTEsmYnDyia/mHvb/wJ3juLoUfEZE4ofBjhu2vGcPXt79qfD7tWrD/GJIiuJBnkIIJNTv37Ocbf2qOyWrsM62puG5cEJOgBcDp1xrh57WHjc7mmiVbRGTcU/iJJZ/PWFTzvxVwYC9MPhQuWQ5z7KaUs9W7F/sv1w4basAINvddfSpZg7S2RHI19pnW1Miu6D6Uvqu9N64wOpmLiMi4pvATK3u98K//MybWA7CdD5fUwNTDTCupbc9+9nZ1c/uSecyZNmXIYyPW0jIanX7tx6u93wcLHJA82eyKREQkihR+YmXNz4zgkzABFv4Qzvw2JIyOwXZzpk1h7sx0s8swzzGfgozZ0LYZXqqD/KvNrkhERKJI4SdWFn4fPnwd7D+BI/LNriZqYjJKK9ISEuHUMni0Ap5bbowCs1jMrkpERKJE4SdWUjPgqofNriJqYjpKKxrmfxEe/znseMOYZHLOIrMrEhGRKFH4kYiI+SitSJuYDvO+COtqjEkPFX5ERMYthR+JmJiO0oqG078B62phYwN8+AYceqzZFYmISBSY0uPW4/FQXV2N0+mkuroar9c75PFut5u8vLx+251OJ16vd8Dvu91u3G63/3q970UGlWmD4z5tvH/uLnNrERGRqDEl/BQXF1NeXk5RURFFRUWUlJQMeqzT6QQYMLwUFxeTkZFBRkYGFosFi8VCdXU1ADU1NeTl5WGxWCgrK8Nms0Xnh5Hx5YzrjD9f/Jux5pqIiIw7MQ8/Ho8n4LPNZsPlcg16fFFREbm5uf22e71e6uvr8fl8/ldVVRXl5eUA5OXl0dbWRltbGw0NDVit1oj+HDJOHXkGTJ8HB/ZB071mVyMiIlEQ8/DjcrnIzMwM2JaZmRnWY6mioiL/e6fTGfAZwGq1KvRIaCwWOONbxvt1K+BAp7n1iIhIxMU8/AzWv6e1tTWk8/QNNV6vl9bW1oBHW16vF6fTidPpxOFw9Gtx6tXZ2UlHR0fAS+LciZ+HqTNgz3Z42Wl2NSIiEmGjZrTXcJ2eh+JwOKiqqgrYVlpa6g9INpuNwsJCWlpa+n23srKSW265Jexry+g33OSK/YbfJybBaaXg+onR8XnelZr0UERkHIl5+LFarf1aeVpbW8N+POX1enG5XP2+7/F4/H2FbDYbHo8Hj8fTr+NzRUUFN9xwg/9zR0cHs2ZpZe/xIJSJF/utRp93Faythg82wKa1YDsvmqWKiEgMxTz82O12ampq+m3Pzw9vyYempqZ+wcftdrNo0SLa2toCth/c1wggJSWFlJSUsK4to1swEy8Ouhp9aoYx6WHjCnj2dwo/IiLjSMzDz8EtLx6Ph/z8fH+AcbvdWK3WAYeme73eAYPOwaHGZrMFPAZzuVwUFRWp83McGtHEi6dfC413w1uPadJDEZFxxJR5furr63E4HDidTmpqaqivr/fvq6ys9M/tA0ZwcTgcA+7rdXBQslqt5OfnU11dTW1tLY2NjQHXEAlKVs4nkx4+e6e5tYiISMRYfD6fz+wiRpOOjg7S09Npb28nLS3N7HKiasPWdi6+4ykevv5s5s5MN7scUwz7d/DOc3DvYkhMhqUbYOphsS9SRESGFcrvb1NafkTGjFmnwREF0L3f6P8jIiJjnsKPyFAsFjjzeuN9492wf4+59YiIyIgp/IgM57iLIWM27G2D9Q+YXY2IiIzQqJnkUCJvq3fvsMO8JQgJicaSF/++yej4nH+1sU1ERMYkhZ9xaqt3L/ZfrmVvV/eQx6UmJZIxOTlGVY1h866Ex38ObZvh9YfhhM+ZXZGIiIRJ4Wecatuzn71d3dy+ZB5zpk0Z9Lh+SzvIwJInQ0EJ/K8anv4tHP9ZLXkhIjJGKfyMc3OmTYnbYewRd2oJPP0b2NpkDIE/6gyzKxIRkTCow7NIsKZMg1O+YLx/5rfm1iIiImFT+BEJxZnXAxZ449+w/XWzqxERkTAo/IiE4pCjP1ny4pk7zK1FRETCovAjEqqzv2P8+VIdtG81txYREQmZwo9IqI7Ih6POgp4ueO4us6sREZEQKfyIhOOspcafzX+EvV4TCxERkVAp/IiE4+hCmHYC7N8NTfeYXY2IiIRA4UckHBYLnPV/xvvnlkPXPnPrERGRoCn8iIRr7mWQPgv2bIcX/2p2NSIiEiSFH5FwJSbBGdcZ75/5LfQMvY6aiIiMDgo/IiMx/8uQmgGtHnj1H2ZXIyIiQVD4ERmJlClwapnx/qlfgc9nbj0iIjIsLWwqAmzcvnvI/RmTk5lpTR1452llxmzP778MG13GSDARERm1FH4krmVMTiY1KZGldeuHPC41KRHXjQsGDkCTMiH/a/DsnfDkrxR+RERGOYWfMWqrdy9te/YPun+4lgwxzLSm4rpxwbB/l0vr1tO2Z//grT9nXAfrauGdZ+DtZ+GoM6JUsYiIjJTCzxi01bsX+y/Xsrdr6NFFqUmJZExOjlFVY9dMa+rgoSZYaTPglCvAfZ/R9+eo+sgUJyIiEafwMwa17dnP3q5ubl8yjznTpgx63JD9VCTyzvo/eOFP8NZj8N5LMP1ksysSEZEBKPyMYXOmTWHuzHSzy5BeWTlw4iWw4UF46tdQ/AezKxIRkQFoqLtIJJ19g/HnKw/Bjo3m1iIiIgNS+BGJpMPnwtGLAR88/WuzqxERkQEo/IhE2rk3GX+++DfwbjG3FhER6UfhRyTSZp0K2edCzwF4+nazqxERkYOow7NIhPSde2nyCd8ke9P/6HH/iTePLuXA5MMBjcATERkNFH5EIqD/3Es+6pOPoYA3eer+H/OzA18GhpkpWkREYkLhRyQCBpp7acqWH8B/v8LVKY9z1lW38saulOFnihYRkahTnx+RCOqde2nuzHRmn/ZZmDGfhO59HL/5/iEnpBQRkdgxJfx4PB6qq6txOp1UV1fj9XqHPN7tdpOXlzfgdrfb7T9n7/twriEScRYLnHuz8X7dChL3eU0tR0REDKY89iouLqa5uRkwQkpJSQn19QOvheR0OrHZbAHBpldNTQ21tbUA2O32gHOEcg2RqDnmQjhsLnywgaxX7gVONbsiEZG4F/Pw4/F4Aj7bbDZcLtegxxcVFQ26Ly8vj7a2NgCsVmvY1xCJmoQEOOdGcH6NrA33MoW5ZlckIhL3Yv7Yy+VykZmZGbAtMzNzwJadYFit1oDgE41riIzICZ+DQ44hcX8HX018zOxqRETiXszDz2B9b1pbW8M6l9PpxOl04nA4/C0+oVyjs7OTjo6OgJdIRCUkwrnlAJRMeISE/btMLkhEJL6NmqHu4XRILi0t9bf62Gw2CgsLaWlpCekalZWV3HLLLSFfWyQkcy+lc3Ul1vYWOl+5D7K/b3ZFIiJxK+YtP1artV8LTGtra79HV8Ho27fHZrPh8XjweDwhXaOiooL29nb/a8sWrcUkUZCQyPbcbwOQ9XItdKr1R0TELDEPP3a7fcDt+fn5IZ3H7XazaNGiftszMzNDukZKSgppaWkBL5FoaLd9lpae6Uzo9MK6WrPLERGJWzEPPzabLeCzx+MhPz/f3yrjdrv7jdbq1fexlc1mo6qqyv/Z5XJRVFSE1Wod9hoipkhI5I4Dlxjvn7lDrT8iIiYxpc9PfX09DoeDgoICGhsbA+bfqayspKCggPJyo4Ooy+WioaEhYF9vyMnPz6e6uhqr1UpLS0vAeYa6hohZ/tVzBlXp/yWl3WO0/pxzo9kliYjEHYvP5/OZXcRo0tHRQXp6Ou3t7aY9Auu7OvhANm7fzdK69Tx8/dnMnZkew8ri04at7Vx8x1ND/n2HcsyTn9rOrCeWQmomLH0JUqZGsXoRkfgQyu/vUTPaSwz9VwcfWGpSIhmTk2NUlURSe85nmfXSHdDaAutWwDk3mF2SiEhcUfgZZQZaHXwgGZOTtTJ4jG3cvjusff0kTIAF5fBQGTzzWyi4Biaqo72ISKwo/IxSvauDi/kyJieTmpTI0rr1Qx4XUmvc3CL4322w8y14vgYW3DzyQkVEJCgKPyLDmGlNxXXjgiH7YUGIrXGJE+C878KDX4dn74BTSyDVOvJiRURkWAo/IkGYaU2N/GPGEy81Wn8+fA2e/R0s1KzPIiKxEPN5fkTkYwkJcH6F8f6538NHoa9vJyIioVPLj0iMBXSOTjuXnKwTSd35Ch8+Ws0HpxphSB3aRUSiR+FHJEYG6zi9KOFT3JP8CpPX38tVz5/MDtJJTUrEdeMCBSARkSgIK/xUVFRQWVkZ6VpExrVBO077zuKjf6xm0ofr+W+Bm6ds32Fp3Xra9uxX+BERiYKwwk99fT05OTnk5+czb968CJckMn4N2nF68Q/hz5dxyGt/4riTSmJfmIhIHAkr/DQ3N5Oens6mTZtYtWoVAJdeemlECxOJKzmL4Mgz4J1nOXT9ncCnzK5IRGTcCmu0V3q6MflednY2O3fupLy8nCVLlrBq1So2b94cyfpE4oPFAgt/CEDmaw8wy/KByQWJiIxfYYWfJUuWcO2115KVlYXH46GhoYG6ujouvfRS2traWLNmTaTrFBn/Zp8FOYuw+A6wdMKDZlcjIjJuhRV+mpubycvLY+fOnVRWVpKdne3f19bWhtfrjVR9IvFlkdH6c0nC06S0vmFyMSIi41NY4ec73/kO11xzDQDt7e08+OCD/sddK1euxGq1Rqo+kfgyYz7t2ReRYPFxWNNtZlcjIjIuhRV+UlJS/O/T09O57LLLcLlcACxfvpyFCxdGpjqROLQ970a6fRbS3n4U3m02uxwRkXEn6NFe7e3trFy5EovFQkNDQ7/9zc3N/tYgEQlfZ8bRrOo+h+IJ/4PVt8BX/2l2SSIi40rQLT/p6enY7XaamppoaWlh48aNAa/y8vJo1ikSV24/cBk9CUmwaS14njC7HBGRcSWkeX6ys7NZvnw5q1evZtGiRQH7NMRdJHK2cihtx3+RrFf+CKt/CtkLjOHwIiIyYkGFn1WrVmG320lLSwNg06ZN3H333f79Xq+XhoYGHn300ehUKRKHts+7nqw362FrM7z2Tzjhc2aXJCIyLgT12OvWW2+lqanJ/3n58uW0tbX5Xz6fj507d0atSJF41D3pUDjzW8aH1T+F7i5zCxIRGSeCavnpG3wAVqxYwfz58wO22e32yFUlIoYzvgWN98DOjfDCnyD/arMrEhEZ88Ia6n5w8AGwqD+CSORNTIMFHw8meGIZ7N9jbj0iIuNAUC0/ffv3DKStrY2VK1fS2NgYkaJEpI+8r8GzvwPv2/DcXXDuzWZXJCIypgXV8nNwH5+DXwA+ny+qhYrErQnJsOhHxvunfgN71L9ORGQkgmr5qaqq6je0/WDq8yMSWRu37/7kQ8Yicg6ZS+qODez4z894/4yfkDE5mZnWVPMKFBEZo4IKP8MFH4CMjIwRFxMPtnr30rZn/6D7A37hSVzKmJxMalIiS+vWB2w/K+Ez/CV5A2kv38fnm09m54TpuG5coAAkIhKisOb5ObgPkOb5Cc5W717sv1zL3q7uIY9LTUokY3JyjKqS0WamNRXXjQsGCMlns/vfTzFl65P81dbAOS1fpG3PfoUfEZEQBRV+br31VqxWq3/B0uXLl7NkyZKAYzTPz/Da9uxnb1c3ty+Zx5xpUwY9To8zZKY1deB74OJboeZcZm19hFMspwFnx7w2EZGxTvP8mGDOtCnMnZludhkyFk0/GU65Al58gO8lPQC+r5pdkYjImDOieX46Ojro6OgI2CYiUbbwB/QkpnBawutMffsxs6sRERlzwgo/7e3tXHDBBVitVjIyMli8eLE/BIlIlKXPZMdJpQAc/vytWvZCRCREYYWfyspKHA4HPT09dHd3s2zZMlauXBnp2kRkEDtOuZYPfWmkdGyCpj+YXY6IyJgSVvgpKCgIGP4+f/588vLyIlaUiAytJ3kKvzlwmfHhiUrY125uQSIiY0hY4WegOX1CmefH4/FQXV2N0+mkuroar9c75PFut3vAcOV2u6murqa6upri4uKA87jdbtxut/96ve9Fxou/dZ/PPusc2NsKT/7K7HJERMaMoOf56auhoQG3243VagWMeX5sNhuzZ88O6qLFxcU0NzcDRjApKSmhvr5+wGOdTic2m23A8OJyuSgvNxZ9rK6uZtGiRf7z1tTUUFtbCxgj0QY7v8hYdYAJfHBqBUc99nV47vfGiu8ZR5ldlojIqBdUy095eTkNDQ2sW7eOdevWkZ6ezo4dO9i4cSMbN25kx44dQS9q6vF4Aj7bbDZcLtegxxcVFZGbm9tvu9vtprKyMuA4t9vtP39eXp5/7bGGhgZ/UBMZT3YdaYfsc6G7E1w/NrscEZExIaiWn5qamqCWuAiGy+UiMzMzYFtmZiZut3vAkDOY3NxcVqxY4f/c+8ir77kVeGTcs1hgcSXUnAOvPASnlsFRZ5hdlYjIqBZUy89wwWfNmjX9Ho0NZrD+Pa2trUF9v6+ioiL/+7q6Oux2e8CjOKfTidPpxOFw9Gtx6tXZ2emfr6jvvEUiY8bhcyH3K8b7/34XenrMrUdEZJQLquVnIKtWrfIHCp/PR1NTE5deemnYhQzX6Xm47zqdTn9/H4DS0lJ/ELLZbBQWFtLS0tLvu5WVldxyyy1hX1vETL0L4SYe/22OeclJ4nvreXftvXiPKQa0VIqIyEDCCj/f/e538Xq9tLa2YrPZ8Hq9lJWVBfVdq9Xar5WntbV1RI+oHA5Hv349Ho/H/xjNZrPh8XjweDzYbLaA71ZUVHDDDTf4P3d0dDBr1qywaxGJhYFWfi9N/AzfS/orSU/8Py5/NIOPmEhqUiLLv5xH1hAL5SogiUi8CSv85OTkUFJSwqZNm7BYLMyePZs1a9YE9V273U5NTU2/7fn5+eGUQnV1NQ6Hwx/CwAg+ixYtoq2tLeDYg/saAaSkpJCSkhLWtUXMMtDK75buAvbXP81hu97hiTPW89rx3+Ybf2rmq/euG/JcqUmJuG5coAAkInEjrPBjs9l4++23yc7O5rbbbuOmm24K6bt9eTwe8vPz/a02vUPoDz4OjMdbfVt3nE4nubm5/uCzcuVKSktLsdlsVFVV+Y9zuVwUFRWpA7SMKwOu/H7RrVD3Jaa9XMu0BSX9AtLBNm7fzdK69bTt2a/wIyJxI6zw0zuvT1tbGzt27GDx4sVYrVYWLlwY1Pfr6+txOBwUFBTQ2NgYMAdPZWUlBQUF/vl7XC4XDQ0NAfuKiorweDwUFxcHnNdqtfr7+uTn51NdXY3VaqWlpUXz/Eh8OO5imH0ObH4SGn7EzOI/KtSIiBzE4vP5fCM9yerVq8nPzyc9PT0SNZmqo6OD9PR02tvbSUtLi+i5N2xt5+I7nuLh689m7syx/3clo9T7L0PNueDrga8+DNnnDHqo7kkRGS9C+f0d1vIWfS/U0dHBokWLxkXwERkXDj/JmO0Z4D8O6D5gbj0iIqNMWOGnvb2dCy64AKvVSkZGBosXL9b8OCKjyfnfh9QM2P4KNN1rdjUiIqNKWOGnsrISh8NBT08P3d3dLFu2jJUrV0a6NhEJ16RMWPgD4/3jP4M9O82tR0RkFAkr/BQUFATM+jx//vwBV10XERPlfQ0OOwn2tcOan5pdjYjIqBFW+MnIyAhqm4iYKCERLqo23jffB9vWm1qOiMhoEdRQ94PX7WpoaPDPxwOfDH2fPXt2pOsTkZE46kyYWwQbnPCfcrj6UWMxVBGROBZU+CkvL6ewsNA/ois9PZ0dO3awY8cO/zE7d+4c0dpeIhIlhT+FN/4NW56HF/8K8640uyIREVMFFX5qamqGXdldREap9JmwoBxcP4HHfgjHXmiMBBMRiVNB9fkZKPh0dHRw9913c/fdd2uYu8hod/p1cMix8NEOWPNzs6sRETFVWB2eN23axMKFC3nsscd47LHHyMvLY/369REuTUQiZkIyXPQL433TPer8LCJxLay1vR588EGampoCtlVUVDBv3rxI1CQi0WBbAHMvgw0PwiM3wtcbzK5IRMQUYbX8ZGdn99uWn58/4mJEJMou+DkkT4WtTfDC/WZXIyJiirDCj8fj6bdt06ZNIy5GRKIsbTqcX2G8d/2ExH1t5tYjImKCsB572e12LrjgAv+szi6Xi6qqqogWJiJRcmoZvPAX2P4Kh62rBD5rdkUiIjEVVsvP/Pnzqampwefz4fP5qK2tZeHChZGuTUSiIXECfPqXAGS+8TfyLa+bXJCISGyF1fJTUFBARUUFy5Yti3Q9IhILR50BuV8B9/3cmnQPXd1fMbsiEZGYCavlp7S0tN9szmvWrIlIQSISI/ZbODAxi2MStpL1cq3Z1YiIxExYLT8Wi4Vrr72WnJwcbDYbra2t1NfX69GXyFgyKZP3Tv8hs55YyjT3b+CMKyDTZnZVIiJRF1bLz7Jly/D5fOzYsYN169axceNGWltbI12biERZ+5xLeKr7RBK6O+GRm8DnM7skEZGoC6vlZ6C1vlavXh2RgkQkhiwWfnjgalYnf4+EltXGBIgnFZldlYhIVIXV8tMbfDo6OvzremnhU5GxaZNvOh/Ou8748N8K2Ku5f0RkfAsr/LS3t3PBBRdgtVrJyMhg8eLFWtxUZAzbccq1cMgxsGc7NPzI7HJERKIqrPDjcDgoKyujp6eH7u5uSkpKqKysjHRtIhIjvsQU+MxvjA/u+2HT/8wtSEQkisIKP3l5eVx22WX+z0VFRVrbS2SsO+pMyP+68f6f34auvebWIyISJWGFn6ysrH7bMjIy/O/Xr18fdkEiYiL7T2DqDGjbBE+oNVdExqewRns1NDTg8XiwWq0AeL1eWlpa/Aue1tfX8+ijj0asSBGJkYlpxtIXf7sCnrkTTrwUZswzuyoRkYgKO/ykp6ezY8cO/7b09HQ2btwIoDl/RMay4y6CEy+BVx6Cf34LSh6HxCSzqxIRiZiIzfPTl+b8ERnjLqyGlsfh/Zfh2Tvh7O+YXZGISMSMaJ6fcPeLyCg3ZRosvtV4/3glfPimufWIiERQWOFHROLAvCshZxF0d8I/vgk93WZXJCISEQo/IjIwiwU++1tIngrvNsJzd5ldkYhIRCj8iMjg0o+AxT8z3q/5GezYaG49IiIRoPAjIkPL/SrYzoMD++Af1+nxl4iMeaaEH4/HQ3V1NU6nk+rqarxe75DHu91u8vLyQjpPqNcQkUFYLPDZOyB5Cmx5Dp6vMbsiEZERCWuo+0A2b97M7Nmzgzq2uLiY5uZmwAgpJSUl1NfXD3is0+nEZrPhdrtDOk8o1xCRYViPhAv+Hzz8HVj9UzhmMWTlmF2ViEhYwg4/69evD5jMsKamhrq6umG/1zsLdC+bzYbL5Rr0+KKiopDPE+o1RCQIeV+DV/4Om9bC36+Fr/0HEhLNrkpEJGRhhZ/LL78cr9frX94C4IUXXgjquy6Xi8zMzIBtmZmZuN1ucnNzg65hqPM0NTVF5Boi0ofFAp+7E+46E7Y8D8/8VpMfisiYFFb4KSwspKSkJGDbgw8+GNR3B+t7E+qSGEOdJ5RrdHZ20tnZ6f/c0dERUh0iccV6JFy4zOj4vObnMKcQDp9rdlUiIiEJq8NzTk7/Z/0DbQtFpDokD3WegfZVVlaSnp7uf82aNSsidYiMW/O+CMdcCD1d8NA34MB+sysSEQlJWC0/LS0t1NTUUFBQAIDP52PlypU0NjYO+12r1dqvBaa1tTXgEVowhjpPKNeoqKjghhtu8H/u6OhQABIZSu/kh3edDh+8DGuXwaIfmV2ViEjQwmr5qampITs7G5/Ph8/nA/D/ORy73T7g9vz8/JBqGOo8oVwjJSWFtLS0gJeIDGPKNLj418b7p34NW4x/+Gz17mXD1vYhX1u9e00sXEQkzJafqqqqfouXDhY4Dmaz2QI+ezwe8vPz/a0ybrcbq9Xa7zggoJP1UOc5uIXn4GuISASc8Dk46XJ4eSU8VMa2Kx7D/tsm9nYNPQlialIirhsXMNOaGqNCRUQChRV+Fi1aREdHBytXrgSM0V/z588P+vv19fU4HA4KCgpobGwMmH+nsrKSgoICysvLAWNUV0NDQ8C+3uHvQ51nqH0iEiEXVcPbT0NrCxNX/5C9XRdz+5J5zJk2ZcDDN27fzdK69bTt2a/wIyKmsfiCfV7Vx6ZNmyguLva3vrzwwgvU19czb968SNcXcx0dHaSnp9Pe3h7xR2AbtrZz8R1P8fD1ZzN3ZnpEzy0Sjt57cqjAApAxOXnwsLLpf3DfZwEfJftv4P+uWzro/a3/D4hItITy+zuslp8HH3yQpqamgG0VFRXjIvyIxJOMycmkJiWytG79kMcN+agq+1w483p45rcsS1rBhx9dCSjYiMjoFVb4yc7O7rct1A7LImK+mdZUXDcuoG3P4MPVg3pUtfAH7H1jNVk7XyF57U0w5+/GqDARkVEorPBz8PIRYDwKE5GxZ6Y1deT9byak8O75v2FW/YVMffcJWFcLp5VFpD4RkUgLK/zY7XYuuOAC/0rrLpeLqqqqiBYmImNLZ8Yx3HrgSn6adB889kOYfQ4cdoLZZYmI9BPWPD/z58+npqbGP89PbW0tCxcujHRtIjLG3N99AbtmnQ/dneC8Gro0p4+IjD5hr+qenZ3NsmXL/J83b97M7NmzI1GTiIxZFt499zaO//uF8OFr8Oj3PpkMUURklAgq/KxatQq73e4fOnb33XcH7Pd6vTQ0NPDoo49GvkIRGVO6Jx0Kl9bAny6BpnvBdp4xIaKIyCgR1GOvW2+9NWBo+/Lly2lra/O/fD4fO3fujFqRIjLG5CyEs5Ya7/95PXjfMbUcEZG+gmr5OXhOnxUrVvSb0TnY5S1EZGzauH13aPsX/gA2Pwlbm+HBErjqkShWJyISvLD6/GRkZPjft7e343K5/CO/RGR8CXYiRDAmQ8yYnGx8SEyCy+6BmnNhy3PG6u/HXh/dYkVEghBW+HG5XFxzzTUApKenc9lll3H33Xf7t4nI+BHMRIi9+i2DkZltdHh+8Ovwv9uYPPkUwhxkKiISMUGHn/b2dlauXInFYvEvNNpXc3Ozwo/IODWiiRBPKjLW/3Lfx6zH/49p/DSyxYmIhCjo8JOeno7dbqeqqoqWlpZ+S1z0rsIuItLPhVWw1c2ED17mjuQ7oOdCsysSkTgW0mOv7Oxsli9fzurVq1m0aFG0ahKR8SYpFYr/SHfNuZzW9Trbm38Fs35udlUiEqfCevien5/PbbfdRkdHBwBr1qzxvxcRGdAhc9h6jrEMzrT1d8Jb/R+fi4jEQljhZ+XKlezYscP/eeHChbhcrogVJSLjU0fOZ7n/QKHxYVUptL9rbkEiEpfCCj9ZWVksW7bMP+OziEiwfnbgS+w9ZC7sbYX6q+DA8KPIREQiKazws27dOnbt2hWwrbGxMSIFicj4tp8k3ln0e5iYDu82Gut/iYjEUFjhp6ysjPnz57N48WKWLFnC0UcfTWFhYaRrE5FxqivtKLj04zUCG1fAi38ztyARiSthhZ/s7Gyam5spKioiPz+fxx57jIULF0a6NhEZz465ABZ813j/r6Xw/gZTyxGR+BH2VKvp6emUlJRw8803s2nTJlatWhXJukQkHixwwBw7HNgLdV+CvV6zKxKROBDW8hYAq1atwuPxAODz+WhqauLSSy+NWGEiEgcSEuDSFVC7ANo2wUNl8IW/GttFRKIkrPDz3e9+F6/XS2trKzabDa/XS1lZWaRrE5FxKnAF+EQmnv97bP+8lIQ3/8v2h39C17kV4S+nISIyjLDCT05ODiUlJWzatAmLxcLs2bNZs2ZNpGsTkXFmqBXiL034Gr9KXs4092+4vnEC373xZgUgEYmKsMKPzWbj7bffJjs7m9tuu42bbrop0nWJyDg09ArxZ7Pj2S4O2XAPlQm/470tnwLrqTGvUUTGv7DCT1tbGzabjba2Nnbs2MHixYuxWq0a8SUiwxpyhfhLqtnd+hpTtj3DkQ3XQM4TMCkzluWJSBwIK/wUFRXR3d0NwLJly1i9ejX5+fkRLUxE4lDiBLYsuosp9xUyq+Ntdv3lq7z9qT9CQv//VGVMTtZjMREJS1jhp6CggIqKCv/oLq3wLiKRkpZ1ONf5buIB34+YuvV/PFNzPbce+GK/41KTEnHduEABSERCFlb4KS0t7Tesfc2aNXrsJSIjNtOayu9u/Co7XprCkWuuo3TCI1y08Dy8xy7xH7Nx+26W1q2nbc9+hR8RCVlY4cdisXDttdeSk5ODzWZj586dOJ1OhR8RiYiZ1lQ490vQvQXWLuOIp77HEbYTYfZZZpcmIuNAWDOJLVu2DJ/Px44dO1i3bh0tLS20trZGujYRiXcLHHDiJdDTZcwA3eoxuyIRGQfCavmpqanp189n9erVESlIRMQvIQE+dxe0bYZtL8ADX4BrGsyuSkTGuLBafhYtWsQvfvELliwxnsGvXr2agoKCiBYmIgJA8iRjyYup02HHG+C8GnoOmF2ViIxhYbX8VFRUYLPZsNvtgBGGVq1aFfTaXh6PB6fTic1mw+PxUFpaitVqDflYp9Ppr+Hg77vdbgByc3PxeDx4vV5yc3ND/2FFxHxp0+GKv8K9F8JGF9OTfgx86qBlMvrTcHgRGUhY4Sc/P5/LLrss7EddxcXFNDc3A0a4KSkpob6+PuRji4uL+x1fVVVFeXk5NTU11NbWAmC32wc9v4iMETPmw6W1sPIrZL32J65N9rG0zjLkVzQcXkQGElb42bRpE2CM+urV2NgYVMtP70rwvWw2Gy6XK+RjvV4v9fX1FBUV+fdXV1dTXl4OQF5eHm1tbUD/ViERGaNO+Cxc8DN47PuUJ/yFKy4+m47siwY8VMPhRWQwYYWf+fPnk5+fT1ZWFg0NDbhcLqqqqoL6rsvlIjMzcLr6zMxM3G53v8dSQx1rs9kCgo/T6Qz4DAo9IuPSGddB2yYsjXdz5BNLYfbRcIRmmBeR4IXd4bm+vp758+fj8/mora0Neo4fr9c74PaBhsoPdWzfYOP1emltbcVmswVsczqdOJ1OHA5Hv1YkERmjLBb4VBUcvRgO7IMHlkDrJrOrEpExJKyWn82bN5Odnc2yZctob2/H5XKRkZHB7Nmzwy5ksKATzLEOh6Nfy1PfjtE2m43CwkJaWlr6nauzs5POzk7/546OjqDrEBGTJE6AonvhDxfC+y/BX4rg6sdgcpbZlYnIGBBWy0/fPjrp6elcdtllg/bbOZjVau3XynNwS04ox3q9XlwuV7/v923p6R0pNlDrT2VlJenp6f7XrFmzgvo5RMRkKVPgypWQPgt2boQHLof9e8yuSkTGgKBbftrb21m5ciUWi4WGhv6TjDU3N3PNNdcMex673U5NTU2/7QOtCh/MsU1NTQMOc1+0aJG/w3Ovg/sPgTFs/4YbbvB/7ujoUAASGSvSpsOXHoR7LoCtTcYcQEv+YrQMiYgMIuiWn/T0dOx2O01NTbS0tLBx48aAV+8oq+H07ZcDRgtNfn6+P8C43W5/C81wx/Yef3CosdlsAY/BXC4XRUVFA7YupaSkkJaWFvASkTHk0GONFqAJE+HN/8LDS8HnM7sqERnFQvrnUXZ2NsuXL2f16tX9lrfYvHlz0Oepr6/H4XBQUFBAY2NjwBw8lZWVFBQU+MPUUMf2OjgkWa1W8vPzqa6uxmq10tLSonl+RMazI08z+gDVfQle+BOkzYBjrjO7KhEZpSw+X3j/RFq/fn1Af5yamhrq6uoiVphZOjo6SE9Pp729PeKtQBu2tnPxHU/x8PVnM3dmekTPLSJA073w8HcA2Hbm/+PMNTn6/5tInAjl93dYD8Yvv/xyvF5vwGOkF154IZxTiYhETv7VsOsDWLuMGc/8kM8mXAecbXZVIjLKhBV+CgsLKSkpCdj24IMPRqQgEZEROe+7sLcV1tXyy6TlbH0nH2YGt+6giMSHsIa65+TkBLVNRCTmPp4E0TvnEpIs3Rzp+ga8/YzZVYnIKBJWy09LSws1NTUUFBQA4PP5WLlyJY2NjREtTkQkLAkJvLvgNpre2IydF4xZoK96GKafYnZlIjIKhBV+ampqsNvt9O0rHWa/aRGR6EhI4rqu/+PpaXdyyM4mDtx3CZsurqMz45h+h2ZMTtbipyJxJKzwU1VV1W+ou91uj0hBIiKRkDE5mYSkVM7b+g0eSP45J+/bRFp9EZfv/xFv+w4PODY1KRHXjQsUgETiRFDhZ82aNQELlx4cfAAsFkvkqhIRGaGZ1lRcNy6gbc9+EvcVsO+RJRzW+joNWb9k08X1dE09AoCN23eztG49bXv2K/yIxImgwk9DQwMFBQVDPtqqq6tj3rx5kapLRGTEZlpTPw406XD1v+APF5G88y2OffRK+Np/jMkQRSTuBDXaq6qqCqvVSkZGxoAvq9VKdXV1tGsVEQnflGnw1X9Cxmxo2wz3fw52bze7KhExQVDhp7S0lI0bN9La2jrga+PGjf3m/RERGXXSZsBX/glpR8CON+G+z5C4d4fZVYlIjAX12Ku4uJjs7OxB96enp1NcXByxokREoibjKKMF6I8Xw4evk/3wErK4weyqRCSGgmr5GaiDczjHiIiMClk5xrw/U6cz0fsWDyT/XC1AInEkrBmeRUTGvKwcuOoRuiYdxrEJ75L9yBWw+0OzqxKRGFD4EZH4lZXDpotX8r4vg4ltb8B9n1EnaJE4oPAjInFtf3o2V+z/AV2TDoMPX4M/XAQd28wuS0SiSOFHROLeJt90Hj/zPvZPmQk732L/isW88forbNja7n9t9e41u0wRiZCwlrcQERkvMiYnk5qUSOnDrcyknL8k38rsXe8w5a+f4cr93/cvhaElMETGD4UfEYlrfZfBANi35ww6H7mCme0tuDKWsfmiv/JK13QtgSEyjij8iEjc+2QZDIB0mPFfuP9zJG1/laMfKcay+D5T6xORyFKfHxGRg02ZBlc9AjNyYW8r2Y98gdMTXjW7KhGJEIUfEZGBTMo0ZoKefQ6JXXu4L6mKqW83mF2ViESAwo+IyGBSpsIXnXQcdQEpli6ObCiFF+vMrkpERkjhR0RkKEkTece+nAe7z8bi64aHSuGZO82uSkRGQOFHRGQ4CRO4qesb7Jh7tfH5se/Do9+Hnh5z6xKRsCj8iIgEwUcC75/+Yyj8qbHh2TthVQkc2G9uYSISMg11FxEJ0sYP98Dsq7CeN5WZa2/GssHJ7rb3ecdeQ0/yVDImJ2seIJExQOFHRGQYvbNAL61b//GWwzg34UZ+n3Q7U7Y+BX+4iGv230R70jSWfzmPrMnJQ55LAUnEXAo/IiLDOHgWaMPZvPfh2Rz16NWcsPdt1mb8jCW7b+Cr93YPeS4tkyFiPoUfEZEgBM4C3bvxXDhqNfylmJQdb7Bq4k9559N3sXvW+QOeY+P23VomQ2QUUPgRERmJjKPg64/Byi+TsOl/zH70avj0bZB/tdmVicggNNpLRGSkUq3wxQfhlCvB1w0Pfwf+WwE9Qz8CExFzKPyIiETChGT4/F1w/g+Mz8/dBQ8sgX3t5tYlIv0o/IiIRIrFAgtuhuL7YEIqbGyAuwuh1WN2ZSLSh8KPiEiknfh5uPo/MHUG7HgDViyETU+aXZWIfMyUDs8ejwen04nNZsPj8VBaWorVag35WLfbDUBubi4ejwev10tubm7I1xARibgZ86FkDfztStjmhvs/R+bpPwLmmF2ZSNwzJfwUFxfT3NwMGCGlpKSE+vr6kI+tqamhtrYWALvdHnCOUK4hIhIVadPha/+Gf14PL9cz49kfc1vSuVgOFADpZlcnErdiHn48nsBn3zabDZfLFdaxeXl5tLW1AQS06oRyDRGRqEpKhUtXwPRT8DX8iKLE//HRw8Xwpb9C+kyzqxOJSzHv8+NyucjMzAzYlpmZ6X+EFeqxVqu13+OsUK4hIhJ1FguceT2bL/wTbb4pTPrwRahdAJufMrsykbgU8/Dj9XoH3N7a2hrysV6vF6fTidPpxOFw+Ft8QrlGZ2cnHR0dAS8RkWjYM/McPrP/Z+zNPB72fAj3fRae/g34fGaXJhJXRs0Mz4MFlqGO7duJ2WazUVhYSEtLS0jXqKys5JZbbgmhUhGR8L3rm4bnc3/nxOYfw0t/g4YfwZZ1xhxBE9UPSCQWYt7yY7Va+7XAtLa2DjgSa7hj+/bt6R3V5fF4QrpGRUUF7e3t/teWLVvC+8FERILkm5AKlyyHi38Nicnw+sNQex68/7LZpYnEhZiHH7vdPuD2/Pz8kI51u90sWrSo377MzMyQrpGSkkJaWlrAS0Qk6iwWY/2vqx+F9CONiRDvtkPTH/QYTCTKYh5+bDZbwGePx0N+fn7A3D29LTpDHWuz2aiqqvLvc7lcFBUV+fcNdQ0RkVFjZi6UrYWjL4AD++DhpfDg12Gf+h+KRIspfX7q6+txOBwUFBTQ2NgYMP9OZWUlBQUFlJeXD3ms1WolPz+f6upqrFYrLS0tAecZ6hoiIqPKpEy4og6evQNW/xQ2PAjbXoCiP8CMeWZXJzLuWHw+ta/21dHRQXp6Ou3t7RF/BLZhazsX3/EUD19/NnNnqmOjSLwJ6r8BW9aB82po32L0B7LfAqd9AxK0GpHIUEL5/T1qRnuJiAgw61Qo+x97ndeS6vkvPFrBrlce5d0Fv6R70qEhny5jcjIzralRKFRk7FL4EREZZbbuT8X+1lVc1jOdH0z4M1PffYJD/3w+N3eV8XjP/JDOlZqUiOvGBQpAIn0o/IiIjDJte/azt6uH/CU3s2XClRzx+Lc5pPU1/pD8C3aecBXvn1ZhDJcfxsbtu1lat562PfsVfkT6UPgRERml5kybwtEzT4XjnoDVt8Bzd5H16h/J+uBpuKQGjsgzu0SRMUk96ERERrukifCpSvjSKpg6HXa+BfcUwuO3QneX2dWJjDlq+RERibGN23eHt3/OIvjms/DITbDBCWur4M3/Gq1A046PQqUi45PCj4hIjGRMTiY1KZGldeuHPTY1KZGMyckD7MiAonvguIvg4RvgvReh5lxY4ICzlkKi/rMuMhz9v0REJEZmWlNx3biAtj37hz122CHqcy+DI8+Ef/0fvPUorPl/8Nq/jAVSDzsxglWLjD8KPyIiMTTTmhq5kVdp0+HKOnipDv7jgPfWQ80COPdmOPs7kbmGyDikDs8iImOZxQKnfAGuex6OvQh6uuCJW6HmXFI/cJtdnciopPAjIjIeTD0cvvAAXHYPTDoEPnwN2z8v4ScT/kjC/qE7WIvEG4UfEZHxwmKBk4rgW41wypVY8HHVhMc42rkIXn/E7OpERg2FHxGR8WZSJlzyezZd+Gfe7plG0p734G9XwgNfAO87ZlcnYjqFHxGRcWrPEeeyeH8VH57yTUiYAG/+B+48FZ76NRwYfsSZyHil8CMiMo7tI4UPTv0ufONpOOosOLAXXD+BmnPA84TZ5YmYQuFHRCQeTDsOrnoEPr/84w7Rr8P9n4OVXwHvFrOrE4kphR8RkXhhscC8K+D6Jji1FCwJ8Oo/4M4CeKIKuvaaXaFITGiSQxGRca7/WmEJcMoPSDniMmY88yMmv/+8MTfQC3+Gwp/AiZcaQUlknFL4EREZp4JbS+zbXJb8PFVpTia0vwPOq+H5WvjUrTAzL1alisSUwo+IyDgVzFpiG7fvZmmdhauLr+PETffB07fDludgxUI4eQks/CFYZ8WuaJEYUPgRERnHgl1L7K3WbnxHf4MJMz7PYY1VZLz1ILxUR88rf2fniVfRdeZ3mHH49BhULBJ9Cj8iInFs4Edjl3GSZT7fT/oLp/Mah75Ug/fFP9N+7k2kn/tNSJpoVrkiEaHwIyISxwZ/NHY2+L7M5i1ryHzm51h3bYQnb4EX74HzvgunXAGJ+hUiY5PuXBGRODfko7EjLmXDEedx3/JKbrX+k6SOd+Gf3zL6Bp3/fTjh85DwyawpW717h+xjBEZrUzCP4kSiReFHRESGlpBIffd5XHX5TZy4dSU8+SvYuRGcX4PDfwXnfQ+OvZCt7fuw/3Ite7u6hzxdalIiy7+cR9bk5EGPUUCSaFL4ERGRoPgmTIQzr4fcr8Jzd8Ezd8L7L8PfroAZ8+k66dvs7Urm9iXzmTNtyoDn2LlnP9/4UzNfvXfdkNdSQJJoUvgREZHQTEwz+v0UlMAzv4V1K2DbC8ze9jX+npyDdd8PmT3j84NOlDjc8PtQApLrxgUKQBIyhR8REQnP5CwovMVoDXr6N/SsW8G8Ay3w6FXw0ilw7s1w7KcD+gRBcMPvg5ufaD1te/Yr/EjIFH5ERGRkJh8CF/w/3rBdxf/++EOumfg4ie+9CHVfgkOPh3NvMjpGhzA6LNj5iUTCoYVNRUQkIrpTD6HywBd584pnjFaflDT48DV48OtwZ57xeEyLp8oooJYfEREJSv8FUgfe3z0xExb+AM74FjSugGfvgrbN8O+b4IllcPo3oOAaSM2IQdUi/Sn8iIjIkIJbINWQmpRIRu8IrVSr0QJ0+nXGivHP3AHt78Can8GTv4b5XzKCUKYtqvWLHEzhR0REhhTMAqm9Bhx+njwJTiuF/K/BKw/B07+BDzbAuhpYVwvHX2y0Es06bdARYiKRZEr48Xg8OJ1ObDYbHo+H0tJSrFZryMe63W5cLhcAjY2NrFixImAfQG5uLh6PB6/XS25ubrR/NBGRcSkiHZATk+Dky+GkYvA8Ac/eCRtd8Nq/jNeMXDjtG3Di52FCSiTKFhmQKeGnuLiY5uZmwAg3JSUl1NfXh3ysy+WivLwcgOrqahYtWuQ/tqamhtraWgDsdvug5xcRkRizWCDnfOO1/TV49nfw0krY5oaHSuGxH0D+1UZL0dTDza5WxqGYj/byeDwBn202m7/1JpRj3W43lZWV/n1FRUW43W7/d/Ly8mhra6OtrY2GhoZBW5ZERMRE046Hz90JN7wKC38IU2fAnu2wdhn8+kSovwo2PwU+n9mVyjgS8/DjcrnIzMwM2JaZmel/TBXssbm5uaxYscK/3ev1+vf3slqtCj0iImPB5EOM+YCWvgRFf4BZp0PPAaOP0B8/DXedDs/Xwr52syuVcSDmj716Q8rBWltbQz62qKjIv62urg673e4PO16vF6fTCRj9gcrKyrDZ+o8o6OzspLOz0/+5o6MjmB9DRESiITEJ5l5qvN5/GRrvMR6Jffg6/OdmaPgRzL2U1COLALUGSXhGzWivwYJOMMf2Bp3e/j5AQMdom81GYWEhLS0t/c5VWVnJLbfcEk7JIiISTYefBJ+53VhC48U6aLrHCEHr/0LO+r/waPIRpG34OmRcBZMyhzubiF/MH3tZrdZ+rTytra0DPp4K9liHw9GvX0/f/kK9I8UO7kMEUFFRQXt7u/+1ZcuW0H8oERGJnonpxlD5bz4HVz8Gp1xJT+JEjk14l+nP3gK/PBZWfgXeaoCebrOrlTEg5uHHbrcPuD0/Pz+sY6urq3E4HNhsNrxeL16vF7fbzaJFi/p97+D+QwApKSmkpaUFvEREZBSyWODI0+CS3/P6Fxv5YddV7D1kLnTvh1f/AX8pMjpJu34C2183u1oZxWIefg7ud+PxeMjPzw+Yn6e3hWa4Y51OJ7m5uf7gs3LlSqxWKzabjaqqKv/3XC4XRUVF6vwsIjJO9KSk86fuC2i55N/wjafgtGshNRN2vQdP/RruOg1qFsBzv4fdH5pdrowypvT5qa+vx+FwUFBQQGNjY8AcPJWVlRQUFPjn7xnsWI/HQ3FxccB5rVarv69Pfn4+1dXVWK1WWlpaNM+PiMh4dfhJcOEyKPwpvPkfo3/QW4/Ce+vhvfX4Hv0+u2eeQ/ucz9Fx1GJ6kqf0O8WAM1OHYat377AzYUfqWhI+i8+nyRP66ujoID09nfb29og/AtuwtZ2L73iKh68/m7kz0yN6bhGReDLsf0/37MTb+FfeefxeTrZ8Mthlny8JV08e/+o+gyd6TqETYx2y1KREXDcuGFEo2erdi/2Xa9nbNXS/o0hcS/oL5ff3qBntJSIiEqrBV5qfwMb0S1jamc2KT1s5ua0B68a/M7FjExcnPsfFic/RnTSFXUcV8nqWna+snULbnv0jCiRte/azt6ub25fMY860/q1LvfUurVs/4mvJyCj8iIjImBPsSvOpSYmccFIuh1nPAt+PjUdhLzvhlYdI7NiKdeNDnL7xIZpSUul5/FOQVwRzFkFS+MFkzrQpat0f5RR+RERkzAl2pfmA/jUWC8yYb7wK/x+82wivPETXy6tI++gD2PiQ8UqaBEcXwvGfNf6cqCAz3ij8iIjImDSileYTEoxh80eexhsnOfjJ7+5l2Qlvc+QHLpJ3bzWGzr/6D3oSktgz40x2HXUBluMvYvoR/VcKkLFH4UdEROJaxpSJvDLhROyvHAdcwEmWTVyYuI4LEpqYwzamvruWqe+uhae/z/5pJ5N8wqfhmMUwfZ7RmiRjjsKPiIjEtf6P0M4BvsI+4E1vC2lvP0bSW/8ho/VFkre/BNtfgicqYep047HY0RdA9oKI1hTMkHnQsPlwKfyIiEjcG/QR2sxcODGXDVuvZfEdj1B7+k5yvE8x5d3/kbjrPXDfD+778VkmcFhWLqWJOaS0HgIzTg27VSjYIfOgYfPhUvgREREZRsbkZPYkZXLpc+mAjWSu4PSEVzk/YT3nJawnO+EDDt2xju8lrYMH/wqPHg45C42X7TyYcmjA+QYfom/sG27IfO9xGjYfHoUfERGRYQw8uux8APYAb7ZvZsqWx8l8by0T330Wdr8PLz5gvAAOOwlsCzj8sDPJTOoKaoh+QXamQk2UKPyIiIgEYcjRZTNPgRNOAZZC1z7Y8hy0rIGNa+CDl/2vQ7iT5qQkPpoxjz0zzmTPjDP4aFouvgkTA06nvjzRpfAjIiISSUkTjUddtvOM9cZ2fwib1oLnCfCsxdL+DpM/aGTyB43wwm8gMQVmnQpHnQVHnQlHFECygk80KfyIiIhE05RD4aQi4+XzQdsm2PQkbH7S+HP3+8b7zU8axyckGR2tjzzDeM06FSZlmvszjDMKPyIiIrFisUCmzXjlfdUIQzs3GsHn7Wfh7aehYytsed54PX278b1Dj4NZp33yysox9ccY6xR+REREzGKxwCFHG6/8q40w5H0bNj9t9Bt65znY8SZ8+Lrxct9nfC81kyMPnc91iVm0b/iIV/efTk/ywCuZB9N/KJh5hcZTPySFHxERkdHCYoGM2cZr/heNbXt2GK1A7zwLWxph2wuwt5W0d1ZzcxLw7Ep6nrGw0TeD9T1zeMln48WeHF73HUkXE4adCyjYeYXG05xCCj8iIiKj2eRD4LhPGy+AA/vh/Zdhy/N8tOk5Jmxzk7x7C8dYtnJMwlYuZy0APQnJtKcdy8M7DoPmt+GEM2Da8ZCYFHD6tj37h51XKNg5hcbKzNQKPyIiImPJhGQ4Ig+OyGPSGd80tu3eDlub4d0m2OaGrW4S9nnJ8L7Mlye8DE+64EkgMRkOOxGmnwKHnwyHn4ylZxYAc6ZNYe7M8FewH0szUyv8iIiIjHVTpsGxFxov8I8q27Lhaf792H/44qw2prS+Ap3txmOzbS/4v3qCJQFX8uEcsjoPsucbEzIediKkzQhpiY5gWpBgdMxMrfAjIiIy3nw8qqw9J4vKA5mcdfHZzJ2RZgyz37beeGz2/kvw3ktY9mxnTsI28GwDz78+OUdqBkw7AaadQEbybPIs+0noPAkYunVopC1IsaDwIyIiMs59spZYFmQsMl7HG1u2vLOJv/3r31SdZeHwvRvhgw2w4y3Y22YMvX/7aWYCD6YA999irGZ/6HFG/6FDj4NDj4VDjgESzfnhwqDwIyIiMk5lTE4mNSkxiLXEcuk+awH0Pobq2mcMsd/+KnzwCrveeYldW15ihqUVdr1nvDyPB5zj2NRD+WvSoUx/KheOPNEIRIfMgfQjISEhSj9heBR+RERExqmBF2Ttr9/oq6SJMP1k4wW8vbWdi+94ijsvzeH4CduY2PYmKW1vkuLdSErbWyTv2UbS3g85I/FDeO1VeK3PyRNTjEkZs3Igaw7WhJkcb+kCzo78DxwkhR8REZFxbMgFWYPU24L0rVUtH2+Z9fFrEQCT2UuOZRsnJL3H909LYuruTcajs9YW6O40WpC2vwrAEcA1E84BvjiimkZC4UdERESGFEoL0tS+QaunG7zvwM4WYxmPnRvZve11Xtg8h2OjXPNQFH5ERERkWGG1ICUkQma28TraDsDmre38+Y6n+EIUagy6LBOvLSIiIhJzCj8iIiISVxR+REREJK4o/IiIiEhcUfgRERGRuKLwIyIiInFF4UdERETiisKPiIiIxBVTJjn0eDw4nU5sNhsej4fS0lKsVmvIx4a7T0REROKXKeGnuLiY5uZmwAgpJSUl1NfXh3xsuPtEREQkfsX8sZfH4wn4bLPZcLlcIR8b7j4RERGJbzEPPy6Xi8zMzIBtmZmZuN3ukI4Nd5+IiIjEt5iHH6/XO+D21tbWkI4Nd5+IiIjEt1GzqvtggSXUY0Pd19nZSWdnp/9ze3s7AB0dHUHXE6zduzro6fyI3bs66OiwRPz8IiIio120fhf2/t72+XzDHhvz8GO1Wvu1wLS2tg44EmuoY8Pdd7DKykpuueWWfttnzZoV5E8UujNuj9qpRURExoRo/S7ctWsX6enpQx5j8QUTkSLI4/EEjMQCyMjIYNOmTf3CyVDHtra2hrXv4Gsc3PLT09NDa2srWVlZWCwDJ9KCggIaGxsH/RkH29/R0cGsWbPYsmULaWlpg35/tBnu5x2N1xnJuUL9brDHh3vfDLdf91XsrjUW76vhjtF9NTquFe65Rut9NdT+aN1bPp+PXbt2MWPGDBIShu7VE/OWH5vNFvDZ4/GQn5/vDyVutxur1YrNZhvy2IGCUjD7DpaSkkJKSkrAtuHmA0pMTBzyf7Dh9qelpY2p/5gM9/OMxuuM5FyhfjfY40d63+i+Mv9aY/G+Gu4Y3Vej41rhnmu03lfB7I/GvTVci08vU/r81NfX43A4/Kmw7/w7lZWVFBQUUF5ePuyx4e4bqeuuu25E+8eaWP08kbzOSM4V6neDPX6k943uK/OvNRbvq+GO0X01Oq4V7rlG630VyrXMEPPHXvGso6OD9PR02tvbx9S/pGR0030l0aD7SqJlNNxbWtsrhlJSUvjxj3/c7zGbyEjovpJo0H0l0TIa7i21/IiIiEhcGTXz/MSr2tpabDYbbreboqKifp28RcLh9XqprKxkyZIl5Obmml2OjCNOpxOAxsZGCgsLsdvtJlck44HT6cRqtdLQ0EBZWVnUfxfqsZeJPB4PLS0t2O12ysvLcTgcZpck40RTU1NIE4eKBMPlcuHxeCgqKqKsrEz/zZKI8Hq9NDY2YrfbKSgooKqqKurXVPiJILfbTV5eXr/tHo+H6upqnE4n1dXV/l9KLpeLnJycgONEDhbqfQVgt9uHnbJBJNR7q/cfar3H5Ofnx7JcGSNCva+sVqs/8PS2/ESbHntFiNPp9D++OljfCRc9Hg8lJSXU19fj9XoDfkHpX+pysHDuK5FgjPTeqqmpicm/0GVsGcl95XK5BpyrLxoUfiKkqKhowO0Ht+bYbDZcLhdgpF0FHhlKOPeVSDBGcm9VV1dTUVGhPorSz0juK7vdTmZmJmVlZTQ0NEStRtBjr6hzuVxkZmYGbMvMzMTtdpOfn8/OnTv929UxVYI11H0lMhLD3Vsulwu73U5ubq6/87PIcIa6r2pra6murgaMRoFYdAFRy0+UDday09rait1up6mpyd+JcMWKFbEtTsasoe4rMP5D0zcIKVhLsIa6t3rXW7TZbHi9Xux2+6D/0hfpa6j76vLLL8flcuFyuWhoaIjJ43uFH5P03gilpaXmFiLjSt+OqRqCLJHk9Xqx2Wy0tbWZXYqMI719X3tDdKz+u6XHXlFmtVr9/xrv1draqpE4MiK6ryRadG9JNIy2+0rhJ8oGS7EaIiojoftKokX3lkTDaLuvFH6ioO+zzYNHQ/TOjaF/RUmodF9JtOjekmgYzfeV+vxESG9HLYDKykoKCgr8zzDr6+txOBwUFBTQ2NiouVgkaLqvJFp0b0k0jJX7SgubioiISFzRYy8RERGJKwo/IiIiElcUfkRERCSuKPyIiIhIXFH4ERERkbii8CMiIiJxReFHRORjbrcbh8OBy+WK6XVra2txOByDLv4oIpGl8CMiEeN2uykrK8NiseBwOKitraW6uhqHw0FGRkbMQ0WoPB4PZWVl/abi7w0ntbW1OJ1OXC4XtbW1eDyeIc/ndDrJy8vDYrFQXV0dsK+6upqMjAzKysooLS2loKCg39pHIhIdmuRQRCLK4/GQk5NDW1tbwNT1brebpqYmSktLzStuGE6nk9zc3ICp+AsLCyksLKS8vNy/ze12k5eXR3NzM7m5uUOes/fYg/8+wAhAvecd6NoiEh1q+RGRiMrMzBxw+3AhYTTqba3pG3zA+FmCDXG9gaa2tjZgu8vl8k/7LyKxpfAjIlHldrv9j4cuv/xyk6sJTWVlJWVlZQPuKy4uDnpRxrKyMmpqagK2ud1utfKImEQLm4pIVPS2dNTV1fkXMLRarbhcLhwOB2VlZdhsNjweDw0NDQGLHLpcLn84aGxspKqqCpfLRVlZGQ6HA4Camhqam5txOp14PB6sVivNzc0UFxf7v+twOMjNzaW+vh6v10teXh5FRUVUVVUNW7/H48Hr9Q4aUA7uFzRQzb1KS0txOBx4PB7/+bRKuoh5FH5EJCpKS0sH/AVvt9ux2+0Bgae+vh63201ubi4ejweHw0FzczMAra2t/r4xdrud5uZmampqyMzMxOv1UlJSQltbGwA5OTk4HA5/MGltbfWfx2q14nA4otLnaKiae69tt9upqamhqqqK2traMdcKJjKe6LGXiERV334tvaO9srKyKCgo8G+3Wq3+kU69wcblcvmPb2xs9B+Xk5PT77yDKS0tZeXKlYARUPLz84Ouu7eF5uARXU6nE4fDgcVioaysDK/XO2TNvcrKyvytYV6vVy0/IiZSy4+IRFXfx0bDDQ3vlZubG/BYqW9rTd/zWa1WSktLqa6uxmq1+h+l9VVaWkptbS2ZmZkhdzAuLy+npqYm4HtFRUUUFRVRXV1NWVmZP8QMVXPv94qLi6mtrVVfHxGTqeVHRCJqsLlqvF6v/7HQUJYsWdJvPqC+nw8+f1ZWFuXl5ZSWlvYblQVGi0swfXwGUlVVRWtra7+RWgeHuOFq7lVUVITD4dAoLxGTqeVHRCLG7Xb7RzVVVlb6H1G1tLRQW1tLRUUFbreburo6wOj/4/F4/N+z2Wzk5uZSVVWFw+HwPxqz2+3+R0q9nYp7W1laWlrIycnBarWSmZlJcXFxv5ai3NzcsANHc3Ozf6LGnJwcWlpaACMY9W31Gajmg1VUVKjVR2QU0CSHIjJm9Yah3haf3o7HB8/S7HQ6gwo/Zk40qEkORWJHj71EZMxqaGgICDk2m40lS5b4l6noDUdjcYJFEYkePfYSkTGrqqqK6upqXC6Xv8WktbWV0tJSXC6Xf6FQhR8R6UuPvUREPtbbH6mwsHDAPjvRUltbS0tLCxUVFRoCLxIDCj8iIiISV9TnR0REROKKwo+IiIjEFYUfERERiSsKPyIiIhJXFH5EREQkrij8iIiISFxR+BEREZG4ovAjIiIicUXhR0REROLK/wfReLeeTnFJxgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig = plt.figure()\n", "ax = fig.subplots()\n", "\n", "ax.stairs(interaction_fraction, bins)\n", "\n", "L = 55 # nuclear interaction length [g cm^-2]\n", "T_line = np.geomspace(T_min, T_max, 100)\n", "R_line = np.sqrt((T_line / Units.GeV + 0.93827) ** 2 - 0.93827 ** 2) # rigidity [GV]\n", "ax.plot(T_line / Units.GeV, 1 - np.exp(-grammage(R_line) / L),\n", " label=rf'$\\Lambda$ = {L} g/cm$^2$')\n", "\n", "ax.set_xscale('log')\n", "ax.set_xlabel('Energy [GeV]')\n", "ax.set_ylabel('Interaction probability')\n", "ax.legend()\n", "\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3.12 (GTsimulation)", "language": "python", "name": "gtsimkernel" }, "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }