{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Blast-wave fit to $\\pi$, $K$, $p$ spectra in Pb-Pb collisions at the LHC" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Transverse momentum spectra in high-energy Pb-Pb collisions can be described by a superposition of thermal sources which move with a velocity $\\beta$ towards the detector. The functional form is given by\n", "\n", "$$\n", "\\frac{1}{2 \\pi p_T}\\frac{dN}{dp_T dy} \\propto m_T \\int_{0}^{1} \\hat r d \\hat r I_0\\left(\\frac{p_T \\sinh\\rho(\\hat r)}{T}\\right) K_1\\left(\\frac{m_T \\cosh\\rho(\\hat r)}{T}\\right) \n", "$$\n", "\n", "where $I_0$ and $K_1$ are [modified Bessel functions](https://en.wikipedia.org/wiki/Bessel_function#Modified_Bessel_functions). Here $T$ is the (kinetic) freeze-out temperature. The transverse mass is defined as $m_T = \\sqrt{p_T^2 + m^2}$. The variable $\\rho$ is the transverse rapidity, i.e., $\\rho = \\mathrm{arctanh}\\,\\beta$ where $\\beta$ is the transverse velocity. The transverse velocity depends and the radial coordinate $r \\equiv \\hat r R$ as \n", "\n", "$$\n", "\\beta(\\hat r) = \\beta_s {\\hat r}^n \\equiv \\beta_s \\left( \\frac{r}{R} \\right)^n\n", "$$\n", "\n", "where $\\beta_s$ is the velocity at the surface of the fireball (maximum radial distance $R$).\n", "\n", "In this problem we perform a simultaneous fit to pion, kaon, and proton transverse momentum spectra in central (0-5\\%) Pb-Pb collisions at 5.02 TeV center-of-mass energy (per nucleon-nucleon pair) in order to extract the model parameters $T$, $\\beta_s$, and $n$. \n", "\n", "As a minimizer we use the numerical minimization software library [iminuit](https://pypi.org/project/iminuit/). You can install iminuit with ``pip install iminuit``. Check out the [iminuit documentation](https://iminuit.readthedocs.io/en/stable/index.html) for more information. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import scipy.integrate as integrate\n", "from scipy.special import k1\n", "from scipy.special import i0\n", "from iminuit import Minuit" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "# define blast-wave integrand (the term we have to integrate over radial coordinate r)\n", "def dndpt_blastwave_integrand(r, pt, m, Tkin, beta_s, n):\n", " \n", " # your code here\n", " R = 1\n", " beta = beta_s*(r/R)**n # beta_s*2/(n+2)\n", " rho = np.arctanh(beta)\n", " mt = np.sqrt(pt**2+m**2)\n", " return pt * r * mt * i0(pt*np.sinh(rho)/Tkin) * k1(mt*np.cosh(rho)/Tkin)\n", " " ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# blast-wave function (dn/dpT dy) \n", "def dndpt_blastwave_not_vectorized(pt, m, Tkin, beta_s, n):\n", " integral, integral_err = integrate.quad(lambda r: dndpt_blastwave_integrand(r, pt, m, Tkin, beta_s, n), 0., 1.)\n", " return integral\n", "\n", "# define a vectorized version of the blast-wave function that can take a numpy array of pt values as input\n", "dndpt_blastwave = np.vectorize(dndpt_blastwave_not_vectorized, excluded=['m', 'Tkin', 'beta_s', 'n'])" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### a) Plot the blast-wave $p_T$ spectrum ($dN/dp_T$) for protons for two differnt surface velocities $\\beta_s = 0.2$ and $\\beta_s = 0.8$ in the range $0 \\le p_T \\le 3\\,\\mathrm{GeV}/c$. Use $T = 0.1\\,\\mathrm{GeV}$ and $n = 1.$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAG1CAYAAADQqgGtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABz/klEQVR4nO3dd3wUdf7H8dem9wKkEAi9hpIAQgBFASMBkSIqYEEET08Fy3F6P7w721k4T4/z1Jx6eogNpQnYaVIUEDAQpCNFQEIKLYGEtN35/bGyEhMggd2dTfJ+Ph77IDs7O/PZ2dndN9/5zncshmEYiIiIiMgl8TK7ABEREZHaQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKpBZYv349vXv3Jjg4GIvFQkZGhtklndeTTz6JxWLhyJEj551v+vTpWCwWfvrpJ/cUJm5x5v2/GFXdJ2raZ0JqBx+zCxCpbVavXs2iRYt46KGHiIiIcPn6SktLuemmmwgICOBf//oXQUFBNG3a1OXrrQ3c/V6Je9T2z4T2W8+llioRJ1u9ejVPPfUUJ06ccMv69uzZw/79+3n44Ye5++67ue2224iMjHTLums6d79X4h61/TOh/dZzKVSJmKigoOCSl5GTkwNQpf+xOmN9Ip6uOp+JqnLlZ0efy9pDoUpcaujQoXTr1o0ZM2aQlJREYGAgzZo141//+tclLffkyZM89NBDNGvWDH9/f6Kjo7nmmmvYsGGDY54z/TZ27NjByJEjCQsLo379+jz44IMUFRWVW96hQ4cYP348MTEx+Pv706FDB6ZNm1ZhvYcOHeLOO+8kLi4Of39/mjdvzr333ktJSYljnY888ggAzZs3x2KxOPp/nKln27Zt3HLLLURGRnLFFVcAsH//fu677z7atm1LYGAg9evX56abbrpgv5E77riDq666CoCbbroJi8VC3759y73+ytYHsHHjRgYNGkRYWBghISFcffXVfPfdd+WWf2YZu3bt4rbbbiM8PJyoqCgee+wxDMPg4MGDDBs2jLCwMGJjY/nnP/9ZhXfvV0eOHLnge/NbVd1WF9pHzvdeVeaHH37AYrHwySefOKalp6djsVjo2rVruXkHDRpEcnJytWqeM2cOFouFFStWVFj3G2+8gcViYcuWLY5pVd1nf8tZ72lV9h+Ab7/9lu7duxMQEEDLli154403Kl3exb6e3zrfZ6KqdV/os1OZqn7fXOrnsir7bXVe4+7du7njjjuIiIggPDyccePGUVhY6JivKt+18iv1qRKX2rx5MwUFBUycOJGJEycSExPDW2+9xaRJk2jTpg2DBw++qOXec889zJkzh4kTJ5KQkMDRo0f59ttv2b59e4UfuJEjR9KsWTOmTJnCd999x8svv8zx48d59913AcjOzqZnz55YLBYmTpxIVFQUX375JXfeeSf5+fk89NBDAGRmZtKjRw9OnDjB3XffTbt27Th06BBz5syhsLAQPz8/RowYwa5du/jwww/517/+RYMGDQCIiopy1HPTTTfRunVrnnvuOQzDAOydalevXs3o0aNp3LgxP/30E6+99hp9+/Zl27ZtBAUFVbodfv/739OoUSOee+45HnjgAbp3705MTEy5eSpb39atW+nTpw9hYWH86U9/wtfXlzfeeIO+ffuyYsWKcoEAYNSoUbRv356///3vfP755zzzzDPUq1ePN954g/79+/P888/zwQcf8PDDD9O9e3euvPLKKr2PF3pvKlPVbXWhfaQq79XZOnbsSEREBCtXrmTo0KEAfPPNN3h5ebFp0yby8/MJCwvDZrOxevVq7r777mrVPHjwYEJCQpg1a5YjFJwxc+ZMOnToQMeOHYGq77PncynvaVX3n82bNzNgwACioqJ48sknKSsr44knnqiwjzrj9Zxxvs9Edff7yj47F1LVffpiP5cX2m+r+xpHjhxJ8+bNmTJlChs2bOCtt94iOjqa559/Hqjed60AhoiL5OfnGxaLxQgLCzO2b9/umJ6Tk2MEBgYaN99880UvOzw83JgwYcJ553niiScMwBg6dGi56ffdd58BGJs2bTIMwzDuvPNOo2HDhsaRI0fKzTd69GgjPDzcKCwsNAzDMG6//XbDy8vLWL9+fYV12Ww2x98vvPCCARj79u2rtJ7KXveZdZxtzZo1BmC8++67532dy5YtMwBj9uzZVV7f8OHDDT8/P2PPnj2OaZmZmUZoaKhx5ZVXVljG3Xff7ZhWVlZmNG7c2LBYLMbf//53x/Tjx48bgYGBxtixY89b79nLvdB78/bbb1fYllXdVlXZR871Xp3L4MGDjR49ejjujxgxwhgxYoTh7e1tfPnll4ZhGMaGDRsMwFiwYEG1a7755puN6Ohoo6yszDHt8OHDhpeXl/G3v/3NMa2q+2xlnPGeVnX/GT58uBEQEGDs37/fMW3btm2Gt7e3cfbPT3VeT2X7xG+d6zNR3f2+Ot9RVd2nnfG5PN9+W93XOH78+HLPv/7664369es77lflcyS/0uE/cZmtW7diGAaTJ0+mXbt2julRUVG0b9+egwcPXvSyIyIiWLt2LZmZmRecd8KECeXu33///QB88cUXGIbB3LlzGTJkCIZhcOTIEcctNTWVvLw8NmzYgM1mY/78+QwZMoTLLruswjqqc3r4PffcU2FaYGCg4+/S0lKOHj1Kq1atiIiIuORm9t+uz2q1smjRIoYPH06LFi0c0xs2bMgtt9zCt99+S35+frnn/O53v3P87e3tzWWXXYZhGNx5552O6REREbRt25a9e/dWubbzvTfnUtVtVZ19pKr69OnDhg0bHH1gvv32W6699lqSkpL45ptvAHvrlcViKXdIp6o1jxo1ipycHJYvX+6YNmfOHGw2G6NGjQKo8j57IRf7nlZ1/7FarSxcuJDhw4fTpEkTx3zt27cnNTXVcd9Zr+dCLma/r+yzeiFV3aed8bn8LWe8xj59+nD06FHHfK74HNVmClXiMps3bwZgzJgxlT4eHBzM//73P0JCQggJCcHPzw9fX1/H/QEDBpxz2f/4xz/YsmUL8fHx9OjRgyeffPKcP+atW7cud79ly5Z4eXnx008/kZuby4kTJ/jvf/9LVFRUudu4ceMAe6fX3Nxc8vPzHYdfLkXz5s0rTDt9+jSPP/448fHx+Pv706BBA6Kiojhx4gR5eXlOXV9ubi6FhYW0bdu2wrzt27fHZrNVCLxn/ygChIeHExAQ4Dj0cPb048ePV7m2870351LVbVWdfeRsJSUlZGVllbtZrVbA/oNTVlbGmjVr2LlzJzk5OfTp04crr7yyXKhKSEigXr161a554MCBhIeHM3PmTMe0mTNnkpSURJs2bQCqvM9eyMW+p1Xdf3Jzczl9+nSF9xgo91xnvZ4LuZj9vrLP6oVUdZ92xufyt5zx2T5zluSZ9/xiP0d1lfpUicts2bKFevXq0bhx43LTi4qK2LZtG/fffz933nmn43/Gd999NyEhIUydOvWCyx45ciR9+vRh3rx5LFq0iBdeeIHnn3+ejz/+mEGDBp33uWe3KtlsNgBuu+02xo4dW+n8nTt3rnJ/iqo4u9XijPvvv5+3336bhx56iF69ehEeHo7FYmH06NGOGp25vury9vau0jTgkrZVVVr8qrqtLnYfWb16Nf369Ss3bd++fTRr1ozLLruMgIAAVq5cSZMmTYiOjqZNmzb06dOH//znPxQXF/PNN99w/fXXX1TN/v7+DB8+nHnz5vGf//yH7OxsVq1axXPPPeeYp6r77IW46z29EGe9HldwxmfnXPu0M5btDBd6zy/lu7YuUqgSl9m8eXOlH9i3336boqIibrjhhnLTf/jhh2o1tzds2JD77ruP++67j5ycHLp27cqzzz5b4YP+448/lvtf4e7du7HZbDRr1oyoqChCQ0OxWq2kpKScc102m42wsLByZ1+dy8WMFD1nzhzGjh1b7kyroqIil4xDExUVRVBQEDt37qzw2I4dO/Dy8iI+Pt7p663M+d6bc6nOtrrQPlLZe5WYmMjixYvLTYuNjQXAz8+PHj168M0339CkSRP69OkD2FuwiouL+eCDD8jOzq7QUb86NY8aNYp33nmHpUuXsn37dgzDcBz6A6q8z7pKVfef4OBgAgMD+fHHHyvMd/Zz3fV63LXfX8w+Xd36zvUd46rXWNXvWtHhP3GhLVu2kJubW+5LNTc3lylTppCamlruLBTDMNiyZUuV/kdqtVorHBKLjo4mLi6O4uLiCvOnpaWVu//KK68A9tPevb29ueGGG5g7d26lgSk3NxcALy8vhg8fzqeffsr3339fYb6z/ycfHBwMUK1A5O3tXaE14JVXXnEcdnImb29vBgwYwIIFC8odksjOzmbGjBlcccUVhIWFOX29lTnfe3MuVdlWVd1HKnuvIiMjSUlJKXcLCAhwPN6nTx/Wrl3LsmXLHKGqQYMGtG/f3nHG1Jnp1an5jJSUFOrVq8fMmTOZOXMmPXr0KPcjXdV91lWquv94e3uTmprK/PnzOXDggGO+7du3s3DhQre/Hnft9xezT1e3vnN9xzj7NVb3u1bUUiUukp2dTW5uLp07d+a6665jwoQJnD59mrS0NKxWa4XxZ/bs2UNxcTEdOnS44LJPnjxJ48aNufHGG0lMTCQkJIQlS5awfv36SsfU2bdvH0OHDmXgwIGsWbOG999/n1tuuYXExEQA/v73v7Ns2TKSk5O56667SEhI4NixY2zYsIElS5Zw7NgxAJ577jkWLVrEVVddxd1330379u05fPgws2fP5ttvv3UMNNitWzcA/vKXvzB69Gh8fX0ZMmTIeV/Tddddx3vvvUd4eDgJCQmsWbOGJUuWUL9+/Qtuj4vxzDPPsHjxYq644gruu+8+fHx8eOONNyguLuYf//iHS9ZZmQu9N5Wpyraq6j5yrvfqzI9WZfr06cOzzz7LwYMHy4WnK6+8kjfeeINmzZpVOORdnffX19eXESNG8NFHH1FQUMCLL75YYZ6q7rOuUtX956mnnuKrr76iT58+3HfffZSVlfHKK6/QoUMHfvjhB7e/Hnfs9xezT1e3vvPtt858jdX9rhU0pIK4xuLFiw3AWLdunfG73/3OCA8PN8LCwoxRo0YZBw4cqDD/nDlzjISEhCotu7i42HjkkUeMxMREIzQ01AgODjYSExON//znP+XmO3PK8LZt24wbb7zRCA0NNSIjI42JEycap0+fLjdvdna2MWHCBCM+Pt7w9fU1YmNjjauvvtr473//W26+/fv3G7fffrsRFRVl+Pv7Gy1atDAmTJhgFBcXl5vv6aefNho1amR4eXk5Tn0+U09ubm6F13T8+HFj3LhxRoMGDYyQkBAjNTXV2LFjh9G0adMLDlFwoSEVKlufYdhP/U9NTTVCQkKMoKAgo1+/fsbq1aurtIyxY8cawcHBFZZ51VVXGR06dDhvvWcv90LvTWWnz1dlW1V1HzGMyt+r88nPzze8vb2N0NDQckMfvP/++wZgjBkzpsJzqvv+nvn8WCwW4+DBg5XWUdV99rec9Z5WZf8xDMNYsWKF0a1bN8PPz89o0aKF8frrrztquJjXcylDKlS17gt9dipT1X3aGZ9Lwzj/fnspr/Hs7Vudz5HYWQzDhT0Qpc566aWXePjhhykoKMDf3/+C8z/xxBOOAe2c5cknn+Spp54iNze3whlNIiLOpO8bAfWpEhfZvHkzLVq0qFKgAnsndbPO8BEREXEGhSpxiS1bttC+ffsqz69QJSIiNZ06qovTGYbB1q1bK4z1cz579uxxYUUiIiKupz5VIiIiIk6gw38iIiIiTqBQJSIiIuIE6lPlJjabjczMTEJDQy/qMiYiIiLifoZhcPLkSeLi4vDyOn9blEKVm2RmZrrtemoiIiLiXAcPHqxwtYTfUqhyk9DQUMD+prjrumoiIiJyafLz84mPj3f8jp+PQpWbnDnkFxYWplAlIiJSw1Sl6446qrtYWloaCQkJdO/e3exSRERExIU0TpWb5OfnEx4eTl5enlqqREREaojq/H6rpUpERETECRSqXEyH/0REROoGHf5zEx3+ExERqXl0+E9ERETEzRSqRERERJxAoUpERETECRSqXEwd1UVEROoGdVR3E3VUFxERqXnUUV1ERETEzRSqRERERJxAF1SWi2MthfxMyPvZ/m9xPpQWQkkhlJwCWxl4+YC3H3j72m8B4RDUAILq228hMRBUD6pwkUoRERFPp1Al52cYcGI/ZG789XZ0D5w8DIbt0pfvFwqRzSCyqf3fqLYQ0xGi24Nv4KUvX0RExE0UqlwsLS2NtLQ0rFar2aVUXeEx2L0Edn4Je5fD6WOVz+ftB+GNIayRvRXKLxh8g+z/evvaW7OspWArhbISKDoBhUeh4Ij939PHoOQkZG+2385m8YL6rSG2E8T3gPhke9jy1i4rIiKeSWf/uYnHn/1XeAx+mAnbP4UD34FxVgj08oXYjhDXxX6LToDweAiOAq9L6JZXWgQnDsDxn+D4Pji2F3K2QdaWyoOcXwg06gZNL4eW/e21KGSJiIgLVef3W6HKTTw2VB3eBOvehM1zoOz0r9OjE6DNQPstLgl8/N1Xk2HYDy9mbbEfbjy4Fn5eb++3dbaAcGh+pT1gtRkIYXHuq1FEROoEhSoP5FGhyjDgx0XwzT/tgeWMmI7QZQy0HWTv4+RJbFbI3QEH1sDeFbBvBRTllZ8nriu0GwztrrP3zVIHeBERuUQKVR7IY0JV9jZY+GfYu8x+38sHEoZB97ugSc+aE0RsVnsr1p6v4cfF9pYsztqV67eGjjdApxuhQWvTyhQRkZpNocoDmR6qTuXCsmdhwzv2s/a8/SD5Hug1AUJj3V+Ps53Mhl1fwo7P7Z3rrSW/PhbbGTrdZL+FNTStRBERqXkUqjyQqaFq63z45AEo/uVwWcIwSHkK6jV3bx3uUnwSdnwBm2fbW7LOdLq3eEGra6DLbfY+WD5+5tYpIiIeT6HKA5kSqspKYPFjsPZ1+/2GiZA6BZpd7p71e4KCo7Btvv3MxrP7jwXVh8Sb4bLxUL+laeWJiIhnU6jyQG4PVScOwuw74ND39vuXPwT9H6vbQxAc+RE2vg+bPoRT2b9Ob9nf3qesTSp4eZtXn4iIeByFKg9y9uCfu3btck+o2vM1zBkPp4/bhx24/g37GX1iZy2D3Yvh+2n2Tu5nOriHx0P330G3OyAwwsQCRUTEUyhUeSC3tVTtXgIf3mzvqB3XBW56x/OGR/AANpvBqZIyCrN247dxOmE7PsKn+AQAZd6B7G58PRlxN5Pr05BSq40Sq0FJmY1Sq40ym4FhGFhtBlbD+PWkQ8uZfyx4e4G3lxc+Xha8vSz4elvw8/HCz9vb/q+PFwG+XgT5eRPo602Ar/3fYH8fQvx9HP8G+HphqSlnZIqI1EIKVR7ILaFq3zfwwY1QVgTth8AN/3PvoJ0mK7XayM4v4nBeETn5xeSeLOLIqRJyTxZz5FQxxwtLOFFYyvHCEvJOl2I7a8/3p4Sh3qu50/tL2nkdBMBqWFho687rZUP4wTCn35WPl4WwQF9CA3wIC/AlLNCHiEA/woN8iQzyJTLIj4ggP+oH+1Hvl1uDEH8C/XQYU0TEGRSqPJDLQ9XBdfDucCgtgNapMOr9Wnd2m2EYHM4r4qejBRw8VsiBY4XsP1rIweOnOXziNLmniqnu3uzn7UVogA8hAT4E+/kQ7OfFZbYfuO7UXDqcXu+Yb29YD75rdAdZkZfh4+2Ft5cFiwW8LRa8LBYMDM40WhkG2H5pySqzGVhtNsqsBsVlNkqsNkrKbBSX2SgutXK61MrpEvu/hSVWCorL7LeSS7tWZJCfN1Gh/jQI8ScqxJ+oUH9iwvyJDgsgJiyAmDB/YsMCCA/0VUuYiMh5KFR5IJeGqsyN8M5Q+2VcWvSFm2eCb4Bz1+FGhmHw8/HTbD+cz86sk+zJPcWe3AL25J6i8AJhw8/bi9hwe2g4O1TUD/GnXrAvEUF+RAb5ERnkS1igLwG+52nRyd4Gq/5tH5rhzLAMjbvDlY9A6wEuHSjVZjMoLLVysqiUk0Vl5J8uJb+olLzTpZwoPHMr4cTpUo4VlDhuRwtKKCmzVXk9gb7eNAwPIPaXW6OIQBpFBBIXEUijSPvf591GIiK1nEKVB3JZqMreCtMH2zulN+kNt80Bv2DnLd/FDMNg/9FCMg6eIOPgCbZl5rM9K5+TRWWVzu/jZaFxZCBN6gfTtF4QTeoFEV/PHgIahgdSP9gPLy8nh53j+2H1y7DhPbAW26fFdYW+j0LrazxqFHrDMDhVXMbRUyXknirmyMlick8Vk5NfTM7JIrLzi8nOLyI7v4jjhaVVWmZUqD+NIwOJj7Rv6yb1gmhSL5im9YOIDQtw/vYWEfEgClUeyGWh6pupsPQpaHQZjJkHAR50seZKFJdZ2XQwj7V7j/L9/uNs+vkEJyr5cff1ttAqOpT2saG0jA6hVXQILaNCaFo/CF9vLxMqxz5q+5pXYf1bUFpon9aomz1ctUrxqHBVFUWlVkcftMN5p8k8UUTmidNknjjNoROnOXT89AUPQ/p5exFfL5DmDYJpVj+Ypg2CaV4/mOZRwTRU4BKRWkChygO59PBf+nT7KOmBkc5drhNYbQabfj7Byl25fLf3KBsPnKD4N4en/Hy86BAXRlJ8BB3jwkmIC6NlVAh+PiaFpws5lWtvuTo7XMX3hJQnoGlvc2tzIsMwyDtdysFjpzl4vJCfj9v7sR04dpoDRwv4+fhpymzn/voI8PWieYMQWjQIpkVUsCMYt4gKJsivDo+XJiI1ikKVi+zbt4/x48eTnZ2Nt7c33333HcHBVTvUZvq1/9zoyKlilu/MZcWuXL75MbdCS1SDEH+SW9SjR7N6dGkSQbvYMM8NUOdzKhdWvWQPV2VF9mmtroGrH4eGnU0tzR3KrDYO5xWx/2gh+44W8NORAvYfLWDvEfuJBKXWc3+1NIoIpHVMCK2jQ2gdE0rrX1ojQwN83fgKREQuTKHKRa666iqeeeYZ+vTpw7FjxwgLC8PHp2r/467toSo7v4ivtmTx5ZbDrNt3rNxwBWEBPvRpHcXlrRqQ3KIeLRoE164zzvIzYcU/YMO7v3Zo73iDPVxFNjO1NLOUWW0cPH6afUdOsfeXkwz25BSwO/cUxwpKzvm8RhGBtI0NpU1MKG1jQ2gX6+GtliJS6ylUucDWrVt58MEHWbJkyUU9vzaGqrzTpXy6KZP5Gw/x/f7j5R7r2CiMfm2juapNFEnxEfiY1Q/KnY7ugWXPwZY59vveftDjbrjyYY88NGuWYwUl7M45xY85J/kx+9d/c04WVzq/j5eFVtEhtIsNpV3DMNo3DCOhYRhRoXVnDDYRMU+dDFUrV67khRdeID09ncOHDzNv3jyGDx9ebp60tDReeOEFsrKySExM5JVXXqFHjx5VWv78+fOZPn06VquVQ4cOceONN/LnP/+5yvXVllBltRms2n2E2ek/s3BrVrnT97s2ieDaTg1J7RBLfL0gE6s02eFNsOgx2LfCfj8gAq76k/36grVs7DBnOlFYwq7sU+zMPsmurJPszDp53jNBG4T4kxBnD1gd4uy3ZvWD1TleRJyqOr/ftaa3aEFBAYmJiYwfP54RI0ZUeHzmzJlMmjSJ119/neTkZF566SVSU1PZuXMn0dHRACQlJVFWVvELfNGiRZSVlfHNN9+QkZFBdHQ0AwcOpHv37lxzzTUuf22eIO90KbPWH2T66p84dOK0Y3qbmBBu7NaYIYlxNAwPNLFCD9IwEW5fYL9k0KLHIHc7LPwzrP8fpD4LbQbWuDMF3SEiyI8ezevRo3k9xzTDMMjMK2LH4Xx2ZJ1k2+F8th/OZ9+RAo6cKmblrlxW7sp1zB/s501CXBgd4sLp1CicTo3DaRkVgreCloi4Qa1pqTqbxWKp0FKVnJxM9+7defXVVwGw2WzEx8dz//33M3ny5Asuc82aNTz55JMsXLgQgBdeeAGARx55pNL5i4uLKS7+9XBGfn4+8fHxNa6lam/uKaav/ok56T87Bt4MD/RlWFIcN3ZrTKdG4bWrf5SzWcsg4wP4+hkoyLFPa9EPBk6B6Pbm1laDFZaUsfOXkLU1037bcTi/wpmlYB/gtENcGJ0ah9O5cTidG0fQXC1aIlJFdbKl6nxKSkpIT0/n0UcfdUzz8vIiJSWFNWvWVGkZ3bt3Jycnh+PHjxMeHs7KlSv5/e9/f875p0yZwlNPPXXJtZtlV/ZJXlqyiy82ZzmmtYsNZdzlzRiW1EijbFeVtw90Gwsdrodv/gnf/Qf2LoPXLofud0K/P6u/1UUI8vOhS5NIujT5dduVWW3syS1gy6E8Nh/KY2tmHlsz8ykssfL9/uPl+v2F+vvQqXE4ifERJDaOICk+gtjwmnsVAhHxDHWipSozM5NGjRqxevVqevXq5ZjvT3/6EytWrGDt2rVVWu6XX37Jn/70JwzDYMCAAUydOvWc89bUlqq9uaf499If+WRTJoZhP0p1dbtoxl/enF4t66tV6lId22s/JLjjM/v9oAaQ8iQk3QpedaAzv5tZbQb7jpzih5/zfrmdYGtm5S1aMWH+JMVHkBQfSZcmEXRuHK7xtERELVWuMmjQIAYNGlSlef39/fH39yctLY20tDSs1ku7QK6r5Zws4oWvdjJ3w8+O4RAGdYzloZQ2tI0NNbe42qReCxj9AexdDl/8CY7shE8m2gdwHfwixHUxu8JaxdvLPjJ/q+hQRnRtDNhbtHZmn+SHn/PYdPAEm37OY2dWPtn5xSzcms3CrdkAeFmgbWwYXZtE0KVJJF2bRNC8tg0HIiJOVSdaqkpKSggKCmLOnDnl+lmNHTuWEydOsGDBApfX5Kln/5VZbbz/3X7+uWgXJ4vtnfRT2kfzUEobOjYKN7m6Wq6sBNa+Diueh5JTgMV+SLD/YxAYYXZ1dUphSRlbDuWTcfA4Gw/Yr0N5OK+ownyRQb50bRJJ16aRdGsaSWLjCAL9dChcpDZTS9Vv+Pn50a1bN5YuXeoIVTabjaVLlzJx4kRzizNR+v5jPDZ/K9sO5wPQuXE4Tw7tQNcm6uPjFj5+cPkD0OkmWPwYbJ5tH5192yf2juwdb9BZgm4S5OdT4czDw3mnyThwgg0HjrPhwAk2H8rjeGEpS3fksHSH/aQDHy8LHeLC6Na0Ht2bRdKtWSTRoeqbJVJX1ZqWqlOnTrF7924AunTpwtSpU+nXrx/16tWjSZMmzJw5k7Fjx/LGG2/Qo0cPXnrpJWbNmsWOHTuIiYlxWV1nH/7btWuXR7RUFZVaefbz7bz33X7AfjbfI6ltublHE516bqZ9K+GzSXD0R/v9Fn1h8FSo39LUssSupMzG1sw8Nhw4Qfr+Y6TvP052fsUBS5vWD+KypvXo0TyS7s3q6ZChSA1XJwf/XL58Of369aswfezYsUyfPh2AV1991TH4Z1JSEi+//DLJycluqc9TDv/tzDrJ/R9uYFf2KQBGXtaY/xvYjvohGp3aI5QVw6qXYeULYC0Gb3+46hHo/aAGDvUwhmFw6MRp0vcfZ/1Px/j+p+PszD7Jb79RG4T40b1ZPUdLWLvYMP3nRaQGqZOhylN5SkuVYRi8v/YAz3y2jeIyGw1C/Jk6MpEr20S5vRapgqN74IuHYc/X9vvRCTD0FWh8mbl1yXnlnS5lw4HjfP/TMdbvO07GzyfKXXUAIDTAh+7N6pHcvB7JLerTMS6sblzGSaSGUqjyQGa2VJ0qLuOPszIcZzVd1SaKf45MpIFapzybYdj7WX01GQqPAhb7tQSvfgz8dUZmTVBUamXzoTzW7TvGun32Q4anistftSHYz5tuzerRs0U9eraoT6dG4fgqZIl4DIUqD2RWqMo5WcS4t9ezNTMfX28L/zewHeMvb67RpGuSgqP2y9z88JH9fngTGPpvaNnf3Lqk2sqsNrYfPsnafUf5bu8x1v90jLzTpeXmCfbz5rJm9oDVq6VaskTMplDlgcwIVXtyTzF22jp+Pn6a+sF+/O+O7iTFR7hl3eICu5fCZw/BiQP2+0m3QeozGpG9BrPZDHZknQlZR1m77xgnCsuHrBB/+5mJvVvWp2eL+iQ0DNN/ikTcSKHKg5jVp2rDgePcOX09xwtLaVo/iHfH96Bp/WCXr1dcrPgUfP00rH0DMCAkFq6bCu0Gm12ZOIHNZrAz+yRr9hxlzd6jrN17lPyi8ocLI4N86dWyPr1aNuDylvV1dqGIiylUeSB3tlQt25HDvR+kU1Rqo3PjcKbd0V39p2qb/WvsI7EftQ8jQqeRMOh5CKp3/udJjWK1GWw/nM/qPUdYveco6/cdo6Ck/NUZGoYHcHmrBlzeqj6Xt2xAdJjGyRJxJoUqD+SuUJW+/xi3vLmW4jIbfdtGkXZLV4L968QYr3VP6WlYPgVWvwKGzd5qNeTf0Hag2ZWJi5Rabfzw8wlW7T7K6j1H2LD/BCXW8mcXtokJ4fJWDejTugHJzevr8y9yiRSqPJA7QtW+IwWM+M8qjheWcnW7aF4f001nEdUFB9fD/Ht/HTQ06VZIfU6XuqkDTpdY+X7/Mb7dfYRVu4+wNTO/3DhZvt4WujSJpE+rBlzRugGdG0dojCyRalKo8iDu6lN19FQxI15bzf6jhXRuHM5Hd/ckyE//Q60zSk/D18/AmjTAgLBGMPw/9lHZpc44XlDC6j1H+Xb3Eb7dncvBY6fLPR4e6MvlrepzRasormzTgMaRQSZVKlJzKFR5IFe2VBWVWrn5ze/YeOAEjSMDmXff5USFqg9VnXTgO3ur1bG99vvJ90DKk+AbaGpZYo4DRwv5Zncu3/54hG93H+Hkbzq9t2gQzJVt7AGrZ4v6+o+YSCUUqjyQq0KV1WZw3wfpLNyaTXigL3Pv7U2r6BCnLV9qoJICWPQYfP8/+/0GbeD616FRN3PrElOVWW38cCiPb3YdYeWPuWQcPIHV9uvXv5+3F92bR3Jl6yiuahtF25hQnVUogkKVR3JVqJq6aCcvf70bP28v3ruzB8kt6jtt2VLD7V4C8yfAqSyweEPfyXDFJPBWa4TYL6mzZs9RVv6Yy8pdufx8vPyhwpgwf65qE8VVbaK5onUDwgN9TapUxFwKVR7E1X2qDh4rZNz09TxwdWuGJsY5bblSSxQeg8//CFs/tt+PT4br34B6zc2tSzyKYRjsO1LAyl25rNiVy5q9Rykq/fWsQm8vC12bRHBVmyj6to2mQ1yYWrGkzlCo8kCu7FNVUmbDz0dn+ck5nLmG4Od/hOJ88AuBQf+ApFtAP4xSiaJSK+t/OsbynfaQtTvnVLnHo0L96ftLwOrTpgFhAWrFktpLocoDmXlBZRHAfnmbj38PB1bb7ycMs49rpcvcyAUcPFbIil25LN+Zy+o9Ryg8awBSHy8L3ZpG0q9dNP3bRdM6OkStWFKrKFR5IIUq8Qg2K6z6Nyx7FmxlENYYRvwXml1udmVSQxSXWVm/7zjLd+awbGcOe3ILyj3eKCKQfu2i6N8uml4tGhDo521SpSLOoVDlgRSqxKMc2gBzfwfH9oDFC/r8Ea76P/DWYRypngNHC1m+K4evd+SwZs9Rist+7Yvl7+NF75b16d8+hv7tomkUoaE9pOZRqPJAClXicYpPwZf/Bxnv2+837g43vAWRzUwtS2qu0yVW1uw9wtc7cvh6ew6ZeUXlHm8XG8rV7aPp3y6GpHiN7i41g0KVB3HXiOoiF23LXPj0D1CcB/5hMPRl6HC92VVJDWcYBjuzT7J0ew7LduSw4cBxzhoWi/rBfvRrF83V7aLp0yaKEF2jUDyUQpUHUkuVeLQTB2DOnfDzOvv9buNg4BSNxC5Oc6yghBW7cliyPYeVu3LLje7u5+1Fz5b1SWkfzdXtY3SYUDyKQpUHUqgSj2ctheVT4JupgAFR7eGmtyG6vdmVSS1TarWxft8xlmzPYemObPYfLSz3eELDMFLaR5OSEEOnRuE6m1BMpVDlgRSqpMbYswzm/R5OZYNPIFz7AnS5TWNaiUsYhsGe3FP2gLU9m/T95Q8TxoYFcPUvAat3y/r4++hsQnEvhSoPpFAlNcqpXHuw2rPUfr/zKBg8Ffx1XUlxrWMFJSzbkcOS7dms2JVbbkysYD9vrmobxTUJMfRvG0N4kM5WFddTqPJAClVS49hssOol+PoZMKxQvzXcNB1iO5pdmdQRRaVW1uw9ypJt2SzZnk12frHjMW8vC8nN63FNQgzXJMTQODLIxEqlNlOo8kAKVVJj7V8Dc8bDyUzwCYBBz0PXsTocKG5lsxlsPpTH4m3ZLN6Wzc7sk+Ue7xAXxjUJMQxIiKV9w1D1wxKnUajyQApVUqMVHIX598CPi+z3O4+G66aCX7C5dUmdtf9oAYu3ZbNoWzbf/3SsXD+s+HqBDEiIJbVDLN2aRmo8LLkkClUeSKFKajzH4cCnwbDZzw4c+Q5EtTW7Mqnjjp4q5usdOSzcms03P+aWG9W9frAfKe1jSO0Yw+WtGqiju1SbQpUH0eCfUuv89K19TKtTWeAbbL8oc+ebzK5KBIDCkjJW7jrCoq1ZLNmeTf5Z42EF+3nTr100qR1i6dcuWgOOSpUoVHkgtVRJrXIqB+beCftW2u93vwtSnwUff3PrEjlLqdXG2r3HWLg1i0Xbssp1dPfz8aJPqwakdowlpX0M9YL9TKxUPJlClQdSqJJax2a1Dxa68gX7/UaX2Q8Hhjc2ty6RSthsBpt+PsFXW7NYuCWLn84acNTLAsnN6zOok70fVkxYgImViqdRqPJAClVSa+1aCB/fDUUnILCe/aLMra42uyqRczIMg13Zp1i4NYuvtmSx7XB+uce7NolgUMeGDOwYS3w9DdVQ1ylUeSCFKqnVju+HWbfD4QzAAv3+DH0eBi8vsysTuaADRwv5authvtqSxYYDJ8o91rFRmCNgtYzS4Ld1kUKVB1KoklqvtAi++j9In26/32YgXP8GBEaYWZVItWTlFTlasNbuO1puqIa2MaEM6hTLtZ0a0jo6RGNh1REKVR5IoUrqjI3vw2eTwFoMkc1h1PsahV1qpKOnilm8LZsvt2Sxes8RSq2//ly2iArm2o4NGdQploSGYQpYtZhClQdSqJI6JTMDZo6BvAP2izIPfRk6jzS7KpGLlldYypLt2Xy55TArdx2hxPrrWFjN6gcxqFNDru3YkI6NFLBqG4UqD6RQJXVO4TGY+7tfL8qcfC8MeBq8dRFcqdlOFpXy9Y4cPv/hMCt2lR9sNL5eINd2bMi1nRrSuXG4AlYtoFDlAjt37mTUqFHl7n/44YcMHz68Ss9XqJI66bfDLjS9wn5R5pAoU8sScZaC4jK+3pHDl1sO8/WOHIpKfw1YjSMDubaTPWAlKmDVWApVLnbq1CmaNWvG/v37CQ6u2rXPFKqkTtv+Gcy7B0pOQlgjGPUeNOpmdlUiTlVYUsbynbl8vvkwX2/P4XSp1fFY48hABndSC1ZNpFDlYjNmzGDBggXMnDmzys9RqJI6L3cnfHQrHP0RvP3tF2TucpvZVYm4xOkSK8t35tgD1o4cCkt+E7A6N+S6TnHqg1UDVOf3u9YMIrNy5UqGDBlCXFwcFouF+fPnV5gnLS2NZs2aERAQQHJyMuvWrbuodc2aNavcoUARqYKotnDX19B2sP3MwAUT4ItHwFpqdmUiThfo582gTg159ZaupP/1Gl67tSvXdW5IkJ83Px8/zRsr9jLk1W+56oXlPP/VDrYcykNtHDVfrWmp+vLLL1m1ahXdunVjxIgRzJs3r1x/p5kzZ3L77bfz+uuvk5yczEsvvcTs2bPZuXMn0dHRACQlJVFWVlZh2YsWLSIuLg6wJ9ZWrVpx4MABAgLOfSmD4uJiiot/vc5Ufn4+8fHxaqkSsdnsfayWP2e/3/RyuOkd9bOSOuFMC9ZnlRwibN4gmMGdGnJdYkPaxoSqBctD1PnDfxaLpUKoSk5Opnv37rz66qsA2Gw24uPjuf/++5k8eXKVl/3ee++xcOFC3n///fPO9+STT/LUU09VmK5QJfKLHV/YL29TchLCGsPo9yGui9lVibhNYUkZy3bk8tkPmXy9I6fcWYQto4K5rnMcQxIb0io61MQqRaHqN6GqpKSEoKAg5syZUy5ojR07lhMnTrBgwYIqL3vIkCHcfffdDBky5LzzqaVKpApyd8FHN8PR3eATAENfhc43mV2ViNsVFJexZHs2n/1wmBU7c8uNg9UuNpTrOjdkSGIcTetX7eQocZ7qhCofN9VkqiNHjmC1WomJiSk3PSYmhh07dlR5OXl5eaxbt465c+decF5/f3/8/f1JS0sjLS0Nq9V6weeI1DlRbez9rObeBT8uhI9/B9mb4eonwMvb7OpE3CbY34dhSY0YltSI/KJSlmyzB6xvfsxlR9ZJdmSd5MVFu+jcOJzrOjdkcOc4GkUEml22/EadCFXOEh4eTnZ2drWeM2HCBCZMmOBIuiLyGwHhcPOH8PUz8O1UWPVvyN4GN7yl6wZKnRQW4MuIro0Z0bUxeYWlLNyaxac/ZLJ6z1F++DmPH37O47kvdtCtaSRDOjfk2s4NiQ49dx9fcZ86EaoaNGiAt7d3hUCUnZ1NbGysSVWJiIOXN6Q8ATEdYMFE2L0Y3kqxh60Grc2uTsQ04UG+jOwez8ju8Rw9VcyXW7L4dFMm6346Rvr+46TvP87fPttGzxb1GZIYx6COsUQE+Zlddp1Va4ZUOB8/Pz+6devG0qVLHdNsNhtLly6lV69eLl13WloaCQkJdO/e3aXrEakVOt0I47+yDxB69Ed482rYvcTsqkQ8Qv0Qf27r2ZSZv+/Fd49ezePXJdClSQQ2A1bvOcqjH2/msmeWMH76euZvPMSp4opns4tr1ZqO6qdOnWL37t0AdOnShalTp9KvXz/q1atHkyZNmDlzJmPHjuWNN96gR48evPTSS8yaNYsdO3ZU6GvlChr8U6QaTuXYL8h88DuweMGAZ6HnvaBTzEUqOHiskM9+OMwnmzLZfjjfMT3A14ur28UwJDGOvm2jCPBVP8WLUSfP/lu+fDn9+vWrMH3s2LFMnz4dgFdffZUXXniBrKwskpKSePnll0lOTnZLfQpVItVUVgyfTYKMX4Yv6XIbDJ4KPv7m1iXiwXbnnOSTTYf5dFMm+44UOKaH+vuQ2jGWIYlxXN6yPj7edeJAlVPUyVDlqc4++2/Xrl0KVSLVYRjw3X9g0V/BsEF8Txj1vgYKFbkAwzDYmpnPJ5sy+XRTJofzihyPNQjx49pODRmWFEfXJpEaZPQCFKo8kFqqRC7Bj0tgzjgozofwJvYO7LEdza5KpEaw2Qy+33+cTzYd4ovNWRwrKHE81igikKFJcQxLiqNdrH6bKqNQ5YEUqkQuUe4u+HAUHNsLfiH2IRfaDjK7KpEapdRq49vdR/g0I5OFW7MoOOtCz21jQhmaFMfQxDji6wWZWKVnUajyIDr8J+JEhcdg9ljYtxKwwDVPQe8H1IFd5CIUlVpZuj2HBRmHWP6bUdy7NY1kWFIc13ZqSIOQut2PUaHKA6mlSsRJrKXw5Z/g+2n2+0m3wnX/Ugd2kUuQd7qUhVuyWLDpEKv3HOVMMvD2snBFqwYM7xLHNQmxhPjXieEty1Go8kAKVSJOZBiw7k346v/sHdib9LZ3YA+ub3ZlIjVedn4Rn27K5JNNmfzwc55jeoCvF9ckxDI8KY4r20ThW0fOIFSo8kAKVSIusHsJzP6lA3tkM7hlFkS1NbsqkVpjb+4pPtmUyYKM8kM0RAb5cm2nhlzfpRHdmtbuMwgVqjyI+lSJuFjODpgxEk7sB/9wGDkdWvY3uyqRWsUwDH74OY8FGZl8+kMmuSeLHY81jgxkWFIcw5Ma0Tom1MQqXUOhygOppUrEhQqOwMzb4MAasHjDtf+A7r8zuyqRWqnMamPN3qPM32g/g/Dsy+EkNAzj+i6NGJoUR0xY7bjIs0KVB1KoEnGxsmL49EHY9KH9fs/7YMAz9os1i4hLnC6xsmR7tuMMwjKbPVJYLHB5ywYMS4pjYMdYQgN8Ta704ilUeSCFKhE3MAz45kX4+hn7/TYD4Yb/gX+IuXWJ1AHHCkr4fPNh5m88RPr+447pZzq4X98ljj6ta14Hd4UqD6I+VSIm2PIxzLsHrMUQ2wlungnhjcyuSqTOOHC0kAUZh5iXcYi9ub92cK8X7Md1ne0d3JPiI2pEB3eFKg+klioRNzu4Hj66GQpyISQWbpkJcUlmVyVSpxiGwZZD+czbeIhPNmVy5NSvHdybNwhmeFIjhneJo2n9YBOrPD+FKg+kUCViguP7YcYoyN0OvsFw4/90aRsRk5T9comc+RsPsXBrNqdLf71ETremkQzv0ojrOjUkMtjPxCorUqjyQApVIiYpyoNZY2HvMrB4QeoU6HmP2VWJ1GkFxWUs2pbFvI2ZfPtjLr/0b8fX20LfttGM6NKI/u2j8fcx/0QThSoPpFAlYiJrKXw+CTa8a7/f4/cwcIrODBTxADn5RXyyKZOPNxxi2+F8x/SwAB+uS4xjhMkDjCpUeSCFKhGTGQas+jcsecJ+X2cGinicnVkn+XjjzyzYmElWfpFjepN6QVzfpRHXd2lEswbu7X+lUOVBdPafiIfZOh/m/R7KiqBhov3SNqGxZlclImex2gy+23uUjzcc4qsthyko+bX/VdcmEVzftTFDOjckIsj1/a8UqjyQWqpEPMjB9fDhKCg8CuHx9mAVk2B2VSJSicKSMhZtzebjjYfK9b/y8/aif7toRnRtRN+20fj5uGb8K4UqD6RQJeJhju2FD26Co7vBPwxGvgst+5ldlYicR05+EQsyMvl44yG2n9X/KjLIl6GJcdzQrTGdG0c4dZ0KVR5IoUrEAxUeg49uhQOrwcsHhrwMXW41uyoRqYLth+3jX83beMhxgecr20Tx7vgeTl2PQpUHUqgS8VBlxTD/Ptgyx36/76Nw1f/ZL14mIh6vzGpj1Z6jfLzhZ1I7xHJtp4ZOXb5ClQdSqBLxYDYbfP00fDvVfj/pVrjuJfDxrEEIRcT9qvP7XbOuaigi4gpeXpDyhD1IWbwh4wOYcZN94FARkSpSqBIROeOycfZrBPoGw97lMG0Q5B0yuyoRqSEUqkREztb6Ghj/JYTEQM5WeCsFsreaXZWI1AAKVS6WlpZGQkIC3bt3N7sUEamqhonwuyXQoC2czIRpA2HvCrOrEhEPp47qbqKO6iI10Onj9iEX9q8CL18YlgaJo8yuSkTcSB3VRUScITASxsyDDiPAVgrz7oZv/mm/jqCIyG8oVImInI+Pv/3Cy73vt99f+jf44mGwWc//PBGpcxSqREQuxMsLBjwDA58HLLD+LZh1O5SeNrsyEfEgClUiIlXV8x4Y+Q54+8OOz+CdoVBw1OyqRMRDKFSJiFRHwjC4fT4ERMDP62DaADj+k8lFiYgnUKgSEamupr3hzkUQHg9Hd8P/BsDhTWZXJSImU6gSEbkYUW3twSqmI5zKhrcHw55lZlclIiZSqKqGf/3rX3To0IGEhAQeeOABNMSXSB0XFgfjvoBmfaDkJHxwE/wwy+yqRMQkClVVlJuby6uvvkp6ejqbN28mPT2d7777zuyyRMRsAeFw21zocL19LKuP74LVr5hdlYiYQKGqGsrKyigqKqK0tJTS0lKio6PNLklEPIGPP9wwDXreZ7+/6K+w8C9gs5lbl4i4Va0JVStXrmTIkCHExcVhsViYP39+hXnS0tJo1qwZAQEBJCcns27duiovPyoqiocffpgmTZoQFxdHSkoKLVu2dOIrEJEazcsLUp+Da56231/zqr3VqqzE3LpExG1qTagqKCggMTGRtLS0Sh+fOXMmkyZN4oknnmDDhg0kJiaSmppKTk6OY56kpCQ6duxY4ZaZmcnx48f57LPP+Omnnzh06BCrV69m5cqV56ynuLiY/Pz8cjcRqeUsFrj8Abj+DfDygS1zYMZNUHzS7MpExA1q5QWVLRYL8+bNY/jw4Y5pycnJdO/enVdffRUAm81GfHw8999/P5MnT77gMmfPns3y5csdoe2FF17AMAz+9Kc/VTr/k08+yVNPPVVhui6oLFJH7F4CM2+H0gJomAi3zoEQdRkQqWl0QeXfKCkpIT09nZSUFMc0Ly8vUlJSWLNmTZWWER8fz+rVqykqKsJqtbJ8+XLatm17zvkfffRR8vLyHLeDBw9e8usQkRqkVQrc8RkENbCPYfW/AXBsr9lViYgL1YlQdeTIEaxWKzExMeWmx8TEkJWVVaVl9OzZk2uvvZYuXbrQuXNnWrZsydChQ885v7+/P2FhYeVuIlLHNOpqH8sqoikc36dBQkVquToRqpzl2WefZfv27WzdupWXX34Zi8VyweekpaWRkJBA9+7d3VChiHic+i3hzsUQ0wkKcu2DhO5dYXZVIuICdSJUNWjQAG9vb7Kzs8tNz87OJjY21qXrnjBhAtu2bWP9+vUuXY+IeLDQGBj3+VmDhN4IW+eZXZWIOFmdCFV+fn5069aNpUuXOqbZbDaWLl1Kr169XLputVSJCGAfJPTWOfYLMltLYPY4WPem2VWJiBPVmlB16tQpMjIyyMjIAGDfvn1kZGRw4MABACZNmsSbb77JO++8w/bt27n33nspKChg3LhxLq1LLVUi4uAbADe+DZeNBwz44mFYNgVq30nYInWSj9kFOMv3339Pv379HPcnTZoEwNixY5k+fTqjRo0iNzeXxx9/nKysLJKSkvjqq68qdF4XEXEpL28YPBVCYmD5FFjxd3tfq2tfsD8mIjVWrRynypOkpaWRlpaG1Wpl165dGqdKRH61/i34/GHAsB8WHPGm/ZI3IuIxqjNOlUKVm1TnTRGROmTrfPvlbKwl0PwqGP0B+IeaXZWI/EKDf4qI1BQdhts7sPuFwL4VMP06KDhidlUichEUqlxMZ/+JyAW1uArGfgpB9eFwBkxLhRMHzK5KRKpJh//cRIf/ROSCjuyG94ZD3kEIbQhj5kF0e7OrEqnTdPhPRKQmatDKflmbqHZw8jBMGwgH15ldlYhUkUKViIgnCYuDcV9C4+5QdALeHQa7l5hdlYhUgUKVi6lPlYhUW1A9uH0BtLwaSgthxmjYMtfsqkTkAtSnyk3Up0pEqq2sBOb9HrZ+DFhg8IvQ/XdmVyVSp6hPlYhIbeDjBze8BZfdCRjw+R9hxQu6rI2Ih1KoEhHxZF7eMPifcOWf7PeXPQML/ww2m7l1iUgFClUupj5VInLJLBbo/xdInWK//91/YMEEsJaZW5eIlKM+VW6iPlUi4hQZH9oDlWGFttfCjW+Db4DZVck52Gw2SkpKzC5DzsPX1xdv73NfzLw6v98+zi5ORERcKOlmCAiH2XfAzi/ggxth9AwI0H/WPE1JSQn79u3DpkO1Hi8iIoLY2FgsFsslLcfpocrb2xur1ersxYqIyBntroUxH9uHWvjpG3h3KNw6F4Lrm12Z/MIwDA4fPoy3tzfx8fF4eam3jScyDIPCwkJycnIAaNiw4SUtz+mhSkcTRUTcoNkVcMen8P4NkLkR3h4IY+ZDeCOzKxOgrKyMwsJC4uLiCAoKMrscOY/AwEAAcnJyiI6OPu+hwAtxenT+bdPZPffc40iAIiLiRHFdYNxXENYIjuyyX4j56B6zqxJwHLHx8/MzuRKpijPBt7S09JKW4/L2yEGDBnHttdfy5JNPUlBQ4OrViYjULVFtYPxXUK+l/ULM01Iha7PZVckvLrWPjriHs94nl4eqYcOGsXbtWmJiYujduzevv/56neq0pyEVRMTlIprA+IUQ2wkKcuHtwXBgrdlVidQ5buk55+3tzeDBg/nDH/7AX//6VxISEvj000/dsWrTTZgwgW3btrF+/XqzSxGR2iwkCsZ+BvE9oTgP3hsOu5eaXZVIneLyUDVw4ECaNGnCLbfcwg8//MArr7zCBx98wPz583nooYdcvXoRkbojMALGzINWKb9ciHkUbJ1vdlUidYbLxqnas2cPLVu25O9//zudOnWq0Jv+f//7H+3atXPV6kVE6ia/IBj9IXx8F2ybD3PGQfFJ6DrG7MpEaj2XtVTdc889NG/enHvvvZf77ruPtLQ0vvnmG/Ly8hzzfPHFF65avYhI3eXjBzdOgy5jwLDBJxNhzX/MrkpqmOnTp5OQkEBQUBDt27fn888/N7skj+eyULV48WL27dvHkCFDyMnJ4dChQzzzzDPUq1ePVq1aAdCiRQtXrV5EpG7z8oahr0Cvifb7Cx+F5X8HjSUoVTB37lwmTpzIY489xpYtW0hNTeWee+4xpZa0tDSaNWtGQEAAycnJrFu37rzzW61WHnvsMZo3b05gYCAtW7bk6aefdss4mlU+/BcTE0ObNm3o1KkTHTt2dPwbGRl53ufNmjWLjIwMx/1FixbxwQcfXHTBIiJSRRYLDHgGAiJg2TOwfAoU5UPqs/bHRM5h6tSp/PGPf+Tmm28GYPDgwUyfPt3tdcycOZNJkybx+uuvk5yczEsvvURqaio7d+4kOjq60uc8//zzvPbaa7zzzjt06NCB77//nnHjxhEeHs4DDzzg0nqrHKoyMzPZuXMnW7ZsYcuWLSxevJht27Zx+vRpOnTowJdfflnp8wICAti2bRsJCQkADBgwgEcffdQ51YuIyPlZLHDVI+AfCl/9H3yXBsX5MOTf9tYscQvDMDhdas4l3AJ9vas1DtPJkyf57rvvmDp1qmPawoUL6dKliyvKO6+pU6dy1113MW7cOABef/11Pv/8c6ZNm8bkyZMrfc7q1asZNmwYgwcPBqBZs2Z8+OGHF2zhcoYqhypvb28SEhJISEhg5MiRrFmzhi+//JJ58+Zx5MiRcz7vf//7H6NGjaJv374kJSWxefPmOjUYWlpaGmlpaboeooiYq+c99mD1yUTY+B6UnILr/2vvfyUud7rUSsLjC01Z97a/pRLkV/Xz0jZt2oSXlxeJiYkUFhYyY8YMXn75ZebNm1ftdT/33HM899xz569v2zaaNGlSYXpJSQnp6enlGmK8vLxISUlhzZo151xe7969+e9//8uuXbto06YNmzZt4ttvvy0XEl2lylv5yJEjLFy4kM8//5yNGzfSrVs3Bg4cyNdff01UVNQ5n9ehQwfS09OZP38+mzdvpmnTpvzlL39xSvE1wYQJE5gwYQL5+fmEh4ebXY6I1GVdbgX/EJhzJ2ydByUFMPJd8A00uzLxIBkZGbRr14709HSuuOIKAEaMGMGgQYMAOHjwIA888AAHDx7kxIkTXHnllUybNq3SZd1zzz2MHDnyvOuLi4urdPqRI0ewWq3ExMSUmx4TE8OOHTvOubzJkyeTn59Pu3bt8Pb2xmq18uyzz3Lrrbeetw5nqFafqsTERB5++GHee++9al1w0M/Pj+DgYJ5++umLKlJERJwkYRjcHAwzb4MfF8H7N8ItH9lbscRlAn292fa3VNPWXR0ZGRl07dqVTp06sXbtWlatWsVf//pX/va3v/Hkk08yevRoHn/8cVJTUzEMg+3bt59zWfXq1aNevXqX+hKqZdasWXzwwQfMmDGDDh06kJGRwUMPPURcXBxjx4516bqrHKpeeOEFtm7dyr///W/+8Ic/EB8fT8eOHR23gQMHnvf5n3/+ueP4poiImKh1Coz5GD4YCfu/hXeHwa1zIMi9P351icViqdYhODNlZGQwZswYwsLC6NGjBz169GDnzp2sXbuWoqIi1q9fz+WXXw7YX9eZPtOVuZTDfw0aNMDb25vs7Oxy07Ozs4mNjT3n8h555BEmT57M6NGjAejUqRP79+9nypQpnhOqhg0bxqRJkxz39+3b5+i0/t57710wVC1evJi//OUvdOnShS5dutCyZcuLr1pERC5N095wx6fw3gg4lA7TB8OY+RAac8GnSu1VVlbG1q1bad++fbnpmzZt4rrrriMgIIDevXvTrl07RowYwfjx40lKSjrn8i7l8J+fnx/dunVj6dKlDB8+HACbzcbSpUuZOHHiOZdXWFiIl1f5EaO8vb3dct3hKoeqe+65h927dxMbG0vnzp0dt/vuu69KfYWuvvpqxo8fz8aNG5k2bRp79uzho48+uqTiRUTkEsR1gXFf2luqcrbB2wPh9gX2CzRLnbRjxw6Kior429/+RlRUFEFBQbz22mv89NNP3HnnnQB8/fXXrFixgrlz59K7d29WrVp1zjMDL/Xw36RJkxg7diyXXXYZPXr04KWXXqKgoMBxNiDAq6++yrx581i61H6tyyFDhvDss8/SpEkTOnTowMaNG5k6dSrjx4+/6DqqqsqhavHixYC9KW/9+vUcOnSITz75hKVLl9KsWTN279593uf37NmTli1b0rJlS2688cZLq1pERJwjuh2M/wreHQrH9sK0QTD2E6ivowl1UUZGBg0bNiQwMJA+ffoQHBzMFVdcwbJlyxyH3Ly8vOjXrx/9+vVjz549bNu2zWXDLYwaNYrc3Fwef/xxsrKySEpK4quvvirXef3IkSPs2bPHcf+VV17hscce47777iMnJ4e4uDh+//vf8/jjj7ukxrNZjGoOMZqUlFRhMM/333+fd999F8DR0/6MDRs28OKLL3LixAk6derkuHxNXXPm7L+8vDzCwsLMLkdEpLy8Q/YWq6M/QnA03D4fYjqYXVWNVVRUxL59+2jevDkBAQFml1NlDz/8MHv27Dnn8AkLFy6kf//++Pr6sm/fPvr378+yZcto1qyZewt1svO9X9X5/a72ZWrODOZ5xoABA9i6des55x81ahTXXXcdzz77LG3atGHEiBEsWrSouqsVERFXCm9kPxQY0wkKcux9rA6lm12VuFlGRgadO3c+5+OzZ8+mXbt2JCUlcdttt/Hmm2/W+EDlTNU+FaG6g3lGRkZyyy23ANClSxeGDx9O//79GTBgwMVXbZIXX3yRt99+G4vFwuTJk7ntttvMLklExHlCouyd1z+4CX5eD+8Mg1tn2Tu1S52wadMm7r333nM+/tZbb7mxmpqn2i1VZwbz7NOnDz/99BNNmzY95yVqAFq2bMk///lPx4UMIyIiLrpYM23evJkZM2aQnp7O+vXrefXVVzlx4oTZZYmIOFdgpP0swGZ9oOSk/ezA3UvNrkrcJDc3lxtuuMHsMmqsaocqsJ/mOHLkSJ5++mn+8Ic/nHdE9eLiYl577TWaNGnCwIED6dixIykpKRw6dOiiizbD9u3b6dWrFwEBAQQGBpKYmMhXX31ldlkiIs7nHwK3zobWA6DsNHw4GrZ/ZnZVIh7vokJVVZzpif/xxx+ze/duduzYwZNPPskf/vAHiouLufnmm2nVqpXT1rdy5UqGDBlCXFwcFouF+fPnV5gnLS2NZs2aERAQQHJycrUurtixY0eWL1/OiRMnOH78OMuXL69xwVBEpMp8A2HUB/YR2K0lMOt2+GG22VWJeDSXDe96rnGtRo0axd133+309RUUFJCYmMj48eMZMWJEhcdnzpzJpEmTeP3110lOTuall14iNTWVnTt3Eh0dDdjPbCwrK6vw3EWLFpGQkMADDzxA//79CQ8Pp2fPnue9VE9xcTHFxcWO+/n5+U54lSIibuTjBzdMA9+JsOlD+Pgue8tV19vNrkzEI1V7SIUL+e2QCmfGtWrfvj3p6elVHtfqUlgsFubNm+cYgRUgOTmZ7t278+qrrwL2UVnj4+O5//77mTx5crXX8bvf/Y7rr7/+nJfeefLJJ3nqqacqTNeQCiJS49hs8MXD8P3/7PcHPg897zG3Jg9XU4dUqKtMG1KhumbNmsW8efN47rnnWLhwIV988YXjmkHuUlJSQnp6OikpKY5pXl5epKSksGbNmiovJycnB4CdO3eybt06UlPPfXHMRx99lLy8PMft4MGDF/8CRETM5OUFg/8JvX65NMhX/wff/NPcmkQ8kMuv7nhmXKszF1wcMGAAjz76qKtXW86RI0ewWq3lRmAFiImJYceOHVVezrBhw8jLyyM4OJi3334bH59zbz5/f3/8/f0vumYREY9iscCAZ8AvBFb8HZb+DUoKoP9j9sdExPWhqrrjWnmy6rRqnZGWlkZaWlq5Q6IiIjWSxQL9HgW/IFj8uL21qqQQBk5RsBLBDYf/qjuulSs0aNAAb29vsrOzy03Pzs52XMvIVSZMmMC2bdtYv369S9cjIuI2lz8I175o/3vta/DZQ/Z+VyJ1nMtbquDXca1GjhzpjtVVuv5u3bqxdOlSR+d1m83G0qVLmThxokvXrZYqEamVetwFvkHwyURInw6lRTAsDbzd8rMi4pFc3lJ1xueff+7S5Z86dYqMjAzHxZ737dtHRkYGBw4cAGDSpEm8+eabvPPOO2zfvp17772XgoICxo0b59K61FIlIrVWl1thxJtg8YYfPoK546GsxOyqRExTa0LV999/T5cuXejSpQtgD1FdunTh8ccfB+wXdn7xxRd5/PHHSUpKIiMjg6+++qpC53UREamGTjfCqPfA2w+2LYBZY+ytVlLjTZ8+nYSEBIKCgmjfvr3Lf8drA7eFqsWLF/OXv/yFOXPmOEZbd6a+fftiGEaF2/Tp0x3zTJw4kf3791NcXMzatWtJTk52eh2/lZaWRkJCAt27d3f5ukRETNFuMNz8IfgEwK6v4MNR9jMDpcaaO3cuEydO5LHHHmPLli2kpqZyzz3mjE1W3auhTJkyhe7duxMaGkp0dDTDhw9n586dbqnVbaHq6quvZvz48QBMmzaN0aNHu2vVptLhPxGpE1qlwK1zwDcY9i6H92+E4pNmVyUXaerUqfzxj3/k5ptvpkWLFgwePJiTJ93/fp65GsoTTzzBhg0bSExMJDU11TFuZGVWrFjBhAkT+O6771i8eDGlpaUMGDCAggLXB32Xj6h+xvTp07njjjucuaoapTojsoqI1FgH18H7N0BxPjS6DG6bC4ERZlfldhVG6DYMKC00pxjfoGoNeXHy5EkiIiJYvXq144jOww8/THp6OsuWLXNVlZVyxtVQcnNziY6OZsWKFVx55ZWVzuOsEdVdeprGhg0bePHFFzlx4gSdOnVyFFyX6Ow/EalT4nvA2E/gvevh0PfwzhAYMx+C65tdmblKC+G5OHPW/edM8Auu8uybNm3Cy8uLxMRECgsLmTFjBi+//DLz5s2r9qqfe+45nnvuufPOs23bNpo0aVJh+pmroZw9YPjFXA0lLy8PgHr16lX5ORfLpaFq1KhRPPXUU7Rv354NGzYwYsQInn/+eQYMGODK1XqUCRMmMGHCBEfSFRGp9eK6wNjP4L3hkPUDTB8Mty+AUJ0YVBNkZGTQrl070tPTueKKKwAYMWIEgwYNAuDgwYM88MADHDx4kBMnTnDllVcybdq0Spd1zz33XHA4pbi4ysOmM66GYrPZeOihh7j88svp2LFjlZ5zKVwaqiIjI7nlllsA6NKlC8OHD6d///51KlSJiNRJsR3hji/g3aGQux2mXwu3fwLhjcyuzBy+QfYWI7PWXQ0ZGRl07dqVTp06sXbtWlatWsVf//pX/va3v/Hkk08yevRoHn/8cVJTUzEMg+3bt59zWfXq1XNLC9G5TJgwgS1btvDtt9+6ZX0u7ajesmVL/vnPf3Km21ZERIQrVyciIp4kqg2M+wLC4+Hobnh7EBzfb3ZV5rBY7IfgzLhV8xJCZ0JVWFgYPXr04A9/+ANjxoxh7dq1FBUVsX79ei6//PJfXpbFcW3fyjz33HOEhISc93ZmPMnfutSroUycOJHPPvuMZcuW0bhx42psgYvnlFC1du3aSqcXFxfz2muv0aRJEwYOHEjHjh1JSUnh0KFDzlhtjaAhFUSkTqvXwh6sIpvDif3w9rVw1PnD6ohzlJWVsXXrVtq3b19u+qZNm7jiiisICAigd+/etGvXjgceeMAx4Pa53HPPPY6Buc91O9fhv7OvhnLGmauh9OrV65zrNAyDiRMnMm/ePL7++mu39uV2ytl/TZo0cSTNys7+KygoYPPmzfzwww+OW2ZmJrt3777UVdcYOvtPROq0/Ex4dxgc2QUhsfbO7FFtza7KZc53Npkn27JlC506deLyyy/nlVdeISgoiNdee42ZM2eyceNGYmNjsdlsrFixgrlz5zJt2jRWrVrlGHjb2WbOnMnYsWN544036NGjBy+99BKzZs1ix44djr5Wr776KvPmzXOEr/vuu48ZM2awYMEC2rb9dR8LDw8nMDCw0vW4/ey/c3U0MwyDY8eOVfrYlVdeyaJFiwgODqZnz5707NmT48ePExkZWdXViohIbRAWB3d8Du8Oh5yt9har2xfY+16Jx8jIyKBhw4YEBgbSp08fgoODueKKK1i2bJnjkJuXlxf9+vWjX79+7Nmzh23btrksVI0aNYrc3Fwef/xxsrKySEpKqnA1lCNHjpQbVPy1114D7IOCn+3tt992+dBOVW6pqlevHu+99x4hISHlphuGwahRoxzHPM9uqQoNDaVx48ZYrVbat29PmzZtWLRoEZs2bXLyy/B8aqkSEQEKj9nPCjy8CQIj7cMtxCWZXJTz1dSWqocffpg9e/acc/iEhQsX0r9/f3x9fdm3bx/9+/dn2bJlNGvWzL2FOpnbW6r69u1LaGhopQNnde7cudLntGrVio0bN3L69Gm2bt3K9u3bz3scVEREarmgevazAN+/4ZdxrIbCmI+h8WVmVybYW6rOdEKvzOzZs7nvvvsIDQ0lODiYN998s8YHKmdy6YjqrVu3ZsGCBbRp0wYfH5eO3uCxzh78c9euXWqpEhEBKMqHGSPhwBrwC4FbZ0PT3mZX5TQ1taUqKiqK119/nRtuuMHsUtzKWS1VLg1VDRo0oGfPnuzYsYOAgAASEhLo0KEDTzzxhDNXWSPo8J+IyG+UFMCHo2HfSvtYSrfMhOaVX0akpqmpoaquqhGXqVm1apWj5/2pU6fYtm0bW7dudeUqRUSkpvALhltmwUe3wp6l8MFNMHoGtLra7MpELorLBv88ffp0uU7tISEh9OjRgx49erhqlSIiUtP4BtqDVJuBUFZkb7natdDsqkQuitNDlWEYzJ07l9atWzN48GA6d+5cbnDQMWPGOHuVIiJSk/kGwMj3oP0QsJbYW662f2p2VSLV5pKWqqeffpr09HQyMjJ4++23ufPOO5kxYwYATu7CJSIitYGPH9z4NnS8AWylMGssbPnY7KoumX7zagZnvU8u6VNVVlbmGJirW7durFy5kuuvv57du3djqeY1iGq6s8/+ExGR8/D2hRFvgrcfbPoQ5t4J1lJIHGV2ZdXm7e0NQElJyTlH8RbPUVhYCICvr+8lLcfpZ/95eXnRt29fXnrppXLjV5WUlDB27Fhmz55NWVmZM1dZI+jsPxGRKrJZ4dMHYeN7gAWGvQpdbjO7qmoxDIMDBw5QWlpKXFwcXl4u68Isl8AwDAoLC8nJySEiIoKGDRtWmMfUIRW8vLw4ePAg3t7elV5FetWqVecdWKy2UqgSEakGmw2++CN8P81+/7p/wWXjza2pmkpKSti3bx82m83sUuQCIiIiiI2NrfRomqlDKlxo56mLgUpERKrJywsGT7UfClz7Onz2B/uhwOTfm11Zlfn5+dG6dWtKSkrMLkXOw9fX13G49lLVzWHORUTE81ksMPDv9r5Wq1+BL/9kD1a9J5pdWZV5eXlp8M86RAd5RUTEc1kscM3T0OeP9vuL/gLf/svcmkTOQaFKREQ8m8UC/R+Dvo/a7y95Elb8w9SSRCqjUCUiIp7PYoG+k6H/X+33lz0Ly54DjQMlHkShysXS0tJISEige/fuZpciIlLzXfkIXPM3+98rnoelf1OwEo/h9CEVpHIaUkFExIm+ew2+mmz/u/f99n5XdWxwaXGP6vx+q6VKRERqnp73wrUv2v9e/Qp89aharMR0ClUiIlIz9bjLPigowNrX4ItHFKzEVApVIiJSc102Hoa+Alhg/Zv2QUI1grmYRKFKRERqtq63w/D/ABZIfxs+fUDBSkyhUCUiIjVf0i0w4r9g8bJfiPmTifYLM4u4kUKViIjUDp1Hwg1vgcUbMj6A+fcqWIlbKVSJiEjt0fEGuHEaePnADzPh47vBWmZ2VVJHKFRV4vrrrycyMpIbb7yxwmOfffYZbdu2pXXr1rz11lsmVCciIufVYTjcNN0erLbMgY9/Z78Qs4iLKVRV4sEHH+Tdd9+tML2srIxJkybx9ddfs3HjRl544QWOHj1qQoUiInJe7YfAyPfAyxe2zoM54xWsxOUUqirRt29fQkNDK0xft24dHTp0oFGjRoSEhDBo0CAWLVpkQoUiInJB7a6F0R+Atx9s/wRm3wFlJWZXJbVYjQtVK1euZMiQIcTFxWGxWJg/f36FedLS0mjWrBkBAQEkJyezbt06p6w7MzOTRo0aOe43atSIQ4cOOWXZIiLiAm1SYfQM8PaHHZ/BrNuhrNjsqqSWqnGhqqCggMTERNLS0ip9fObMmUyaNIknnniCDRs2kJiYSGpqKjk5OY55kpKS6NixY4VbZmamu16GiIi4S+tr4OYZ4BMAu76EmWOgtMjsqqQW8jG7gOoaNGgQgwYNOufjU6dO5a677mLcuHEAvP7663z++edMmzaNyZPtF9/MyMi4qHXHxcWVa5k6dOgQPXr0qHTe4uJiiot//d9Qfn7+Ra1TREScoFUK3PwRfHgz/LgQZt4Koz4A3wCzK5NapMa1VJ1PSUkJ6enppKSkOKZ5eXmRkpLCmjVrLnn5PXr0YMuWLRw6dIhTp07x5ZdfkpqaWum8U6ZMITw83HGLj4+/5PWLiMglaNkPbp0FPoGwewl8dDOUnja7KqlFalWoOnLkCFarlZiYmHLTY2JiyMrKqvJyUlJSuOmmm/jiiy9o3LixI5D5+Pjwz3/+k379+pGUlMQf//hH6tevX+kyHn30UfLy8hy3gwcPXvwLExER52h+Jdw6G3yDYM/X9parkkKzq5JaosYd/nOHJUuWnPOxoUOHMnTo0Asuw9/fH39/f9LS0khLS8Nq1ai+IiIeoXkfuHUOfHAT7F0GH462Hxr0CzK7MqnhalVLVYMGDfD29iY7O7vc9OzsbGJjY02pacKECWzbto3169ebsn4REalEs8vhtjngFwL7VsCMkVBSYHZVUsPVqlDl5+dHt27dWLp0qWOazWZj6dKl9OrVy5Sa0tLSSEhIoHv37qasX0REzqFpb7htrj1Y/fQNzBilYCWXpMaFqlOnTpGRkeE4g2/fvn1kZGRw4MABACZNmsSbb77JO++8w/bt27n33nspKChwnA3obmqpEhHxYE16wm0fg1+oPVh9cBMUnzK7KqmhLIZhGGYXUR3Lly+nX79+FaaPHTuW6dOnA/Dqq6/ywgsvkJWVRVJSEi+//DLJyclurrS8/Px8wsPDycvLIywszNRaRETkNw6uh/dHQHE+NOlt78zuH2J2VeIBqvP7XeNCVU1zdkf1Xbt2KVSJiHiqn9PhveuhOA+a9PolWFW8ZJnULQpVHkgtVSIiNcChdHj3l2AV39PemV3Bqk6rzu93jetTJSIi4jKNusHt8yEgHA5+B+/fAEW6IoZUjUKVi+nsPxGRGqZRV7h9AQREwMG1ClZSZTr85yY6/CciUsNkZsC7w6DoBDTubh9+ISDc7KrEzXT4T0RE5FLFJcHYTyAwEn5eD++NgKI8s6sSD6ZQJSIici4NE+H2X4LVoe8VrOS8FKpcTH2qRERquIadYeynZwWr6+H0CbOrEg+kPlVuoj5VIiI1XNZmeGconD4GcV1hzDwIjDC7KnEx9akSERFxtthOv/SxqgeZG+C94XD6uNlViQdRqBIREamq2E6/HAqsB5kbfzkUqGAldgpVIiIi1RHb0R6sguorWEk5ClUupo7qIiK1UGxH+1mBClZyFnVUdxN1VBcRqYWyt8I7Q6DwKMR1+aXzeqTZVYkTqaO6iIiIO8R0KH8o8N3harGqwxSqRERELsXZwepwhg4F1mEKVSIiIpfqty1WGiC0TlKoEhERcYYKwWq4glUdo1DlYjr7T0SkDlGwqtN09p+b6Ow/EZE6JGsLvDv0l7MCdUmbmkxn/4mIiJjpzDhWjkvaqI9VXaBQJSIi4gpnRl4/E6zeHwFFeWZXJS6kUCUiIuIqsR1/vQjzoXR7i5WCVa2lUCUiIuJKsZ1+CVaRvwQrtVjVVgpVIiIirhbb6Zc+VpFw6Ht4/wYoyje7KnEyhSoRERF3aNgZbl8AARHw83oFq1pIocrFNE6ViIg4NEy0HwoMiICf18EHN0LxSbOrEifROFVuonGqRETEITPDPo5VUR7EJ8Ntc8E/1OyqpBIap0pERMSTxSX9cigwHA6uhffVYlUbKFSJiIiYIa4LjJn/S7D6Dj4YCcWnzK5KLoFClYiIiFka/XIJG/9wOLAaZoyEkgKzq5KLpFAlIiJipkbdfglWYbB/FcwYpWBVQylUiYiImK3xL8HKLxR++uaXYFVodlVSTQpVIiIinqDxZeWD1YejofS02VVJNShUiYiIeIr47jDmY/ALgX0r4MObFaxqEIWqSlx//fVERkZy4403VusxERGRSxbfwz5ulW8w7F0GH90CpUVmVyVVoFBViQcffJB333232o+JiIg4RZOecNsce7Da8zXMvFXBqgZQqKpE3759CQ2tfGTb8z0mIiLiNE17w62zwTcIdi+BmbdBWbHZVcl51LhQtXLlSoYMGUJcXBwWi4X58+dXmCctLY1mzZoREBBAcnIy69atc3+hIiIil6rZ5XDLLPAJhN2LYdbtClYerMaFqoKCAhITE0lLS6v08ZkzZzJp0iSeeOIJNmzYQGJiIqmpqeTk5DjmSUpKomPHjhVumZmZ7noZIiIiVdO8D9wyE3wCYNdXMPsOKCsxuyqphI/ZBVTXoEGDGDRo0Dkfnzp1KnfddRfjxo0D4PXXX+fzzz9n2rRpTJ48GYCMjAyX11lcXExx8a//m8jPz3f5OkVEpJZqcRXc/JF9mIWdX8CccXDTdPD2NbsyOUuNa6k6n5KSEtLT00lJSXFM8/LyIiUlhTVr1ri1lilTphAeHu64xcfHu3X9IiJSy7TsB6NngLc/7PjMHqyspWZXJWepVaHqyJEjWK1WYmJiyk2PiYkhKyurystJSUnhpptu4osvvqBx48blAtn5Hjvbo48+Sl5enuN28ODBi3tRIiIiZ7S6+pdg5QfbP4W5vwNrmdlVyS9q3OE/d1iyZMlFPXY2f39//P39SUtLIy0tDavV6qzyRESkLmudAqM+sI9ftW0+eHnD9f8Fb/2km61WtVQ1aNAAb29vsrOzy03Pzs4mNjbWlJomTJjAtm3bWL9+vSnrFxGRWqjNABj1Hnj5wpa5MP8esOk/72arVaHKz8+Pbt26sXTpUsc0m83G0qVL6dWrlyk1paWlkZCQQPfu3U1Zv4iI1FJtB8HId8DLBzbPhvn3KViZrMaFqlOnTpGRkeE4g2/fvn1kZGRw4MABACZNmsSbb77JO++8w/bt27n33nspKChwnA3obmqpEhERl2k3GG58Gyze8MNH8Mn9YLOZXVWdVeMOwH7//ff069fPcX/SpEkAjB07lunTpzNq1Chyc3N5/PHHycrKIikpia+++qpC53UREZFaIWEo3Pg/mHMnZHwAFi8Y8jJ41bh2kxrPYhiGYXYRtdnZHdV37dpFXl4eYWFhZpclIiK1zeY58PFdYNig2x0w+F8KVk6Qn59PeHh4lX6/FarcpDpvioiIyEX5YRbM+709WHX/HVz7IlgsZldVo1Xn91sRVkREpLboPBKG/QewwPq34Mv/A7WduI1ClYvp7D8REXGrpJth2Kv2v9e9AQv/rGDlJjr85yY6/CciIm6V/g58+oD97973wzVP61DgRdDhPxERkbqu21gYPNX+9+pXYOnf1GLlYgpVIiIitVX3O2HQC/a/v50Ky6eYW08tp1DlYupTJSIipkq+G1J/CVMrnoflz5tbTy2mPlVuoj5VIiJiqtWvwKK/2v/u/xhc+bC59dQQ6lMlIiIi5fW+H1KetP/99dOw6t+mllMbKVSJiIjUFVf8Afr90lq1+HFYk2ZuPbWMQpWIiEhdctUjcNVk+98L/wxr/2tuPbWIQpWLqaO6iIh4nL6Toc8vfaq+fATW/8/cemoJdVR3E3VUFxERj2IYsORJWPWS/f6Ql+1jW0k56qguIiIi52ex2Duu95xgv//pg7DxA1NLqukUqkREROoqiwVSn4UevwcMWDABNs00u6oaS6FKRESkLrNYYNDzcNmdgAHz74HNc8yuqkZSqBIREanrLBa49kXoejsYNvj4btg63+yqahyFKhfT2X8iIlIjeHnBdf+GpFvBsMLcO2HH52ZXVaPo7D830dl/IiJSI9isMO8e2DwLvHxh1PvQdqDZVZlGZ/+JiIjIxfHyhuGvQYcRYCuFWWPgxyVmV1UjKFSJiIhIed4+MOK/0H4IWEvgo1tgzzKzq/J4ClUiIiJSkbcv3DAN2gwCazF8eDPs+8bsqjyaQpWIiIhUzscPRr4Dra6BstMwYxTsX2N2VR5LoUpERETOzcff3lm9RT8oLYAPboSD682uyiMpVImIiMj5+QbA6BnQrA+UnIL3R8ChDWZX5XEUqlxM41SJiEit4BcEt8yEJr2hOB/eGw6HN5ldlUfROFVuonGqRESkVig+Ce+NgJ/XQWA9uOMziOlgdlUuo3GqRERExDX8Q+G2ORDXFU4fg3eGQs4Os6vyCApVIiIiUj0B4TDmY4jtDIVH4N2hcGS32VWZTqFKREREqi8wEm5fANEd4FQ2vDMEju01uypTKVSJiIjIxQmqZw9WUe3gZKb9UODx/WZXZRqFKhEREbl4IVFw+ydQvxXkHbS3WOX9bHZVplCoEhERkUsTGgNjP4XI5nBiv73FKv+w2VW5nUKViIiIXLqwOHuwimgCx/bYO6+fyjG7KrdSqKrE9ddfT2RkJDfeeGO56QcPHqRv374kJCTQuXNnZs+ebVKFIiIiHigi3h6swhrBkV3w7jAoOGp2VW6jUFWJBx98kHfffbfCdB8fH1566SW2bdvGokWLeOihhygoKDChQhEREQ8V2cwerEJiIWcbvDcMCo+ZXZVbKFRVom/fvoSGhlaY3rBhQ5KSkgCIjY2lQYMGHDtWN3YUERGRKqvf0h6sgqMha7P9WoFFeWZX5XI1LlStXLmSIUOGEBcXh8ViYf78+RXmSUtLo1mzZgQEBJCcnMy6deucXkd6ejpWq5X4+HinL1tERKTGi2pjH24hqD5kboT3b7Bf4qYWq3GhqqCggMTERNLS0ip9fObMmUyaNIknnniCDRs2kJiYSGpqKjk5v3aWS0pKomPHjhVumZmZVarh2LFj3H777fz3v/91ymsSERGplWIS7MEqIAJ+Xg8fjISS2tttxsfsAqpr0KBBDBo06JyPT506lbvuuotx48YB8Prrr/P5558zbdo0Jk+eDEBGRsZFr7+4uJjhw4czefJkevfufd75iouLHffz8/Mvep0iIiI1VmwnGDPP3mn9wGr4cDTcMgt8A82uzOlqXEvV+ZSUlJCenk5KSopjmpeXFykpKaxZs+aSl28YBnfccQf9+/dnzJgx5513ypQphIeHO246TCgiInVWo65w28fgFwL7VsLM26Cs+MLPq2FqVag6cuQIVquVmJiYctNjYmLIysqq8nJSUlK46aab+OKLL2jcuLEjkK1atYqZM2cyf/58kpKSSEpKYvPmzZUu49FHHyUvL89xO3jw4MW/MBERkZouvjvcOgd8g2D3Epg1FspKzK7KqWrc4T93WLJkSaXTr7jiCmw2W5WW4e/vj7+/P2lpaaSlpWG1Wp1ZooiISM3TtBfc/BHMGAm7voS5d8KNb4N37YgjtaqlqkGDBnh7e5OdnV1uenZ2NrGxsabUNGHCBLZt28b69etNWb+IiIhHaXEVjP4AvP1g+ycw726w1Y6Gh1oVqvz8/OjWrRtLly51TLPZbCxdupRevXqZWJmIiIg4tEqBke+Clw9smQsLJkIVjwR5shoXqk6dOkVGRobjDL59+/aRkZHBgQMHAJg0aRJvvvkm77zzDtu3b+fee++loKDAcTagu6WlpZGQkED37t1NWb+IiIhHajvIfujP4g2bZsBnD9X4YGUxDMMwu4jqWL58Of369aswfezYsUyfPh2AV199lRdeeIGsrCySkpJ4+eWXSU5OdnOl5eXn5xMeHk5eXh5hYWGm1iIiIuIxNs+Bj+8CwwY97oZB/wCLxeyqHKrz+13jQlVNc3ZH9V27dilUiYiI/FbGhzD/XsCAXhNhwDMeE6wUqjyQWqpERETOI/0d+PQB+999/gj9H/OIYFWd3+8a16dKREREaqFuY+HaF+1/f/NPWPmCufVcBIUqERER8Qw97oIBz9r/XvYsfPuSqeVUl0KVi+nsPxERkWroPRGuftz+95In4LvXzK2nGtSnyk3Up0pERKQalj0HK563/z14KnS/05Qy1KdKREREara+j8LlD9r//nwSbHzf3HqqQKHKxXT4T0RE5CJYLJDyFCTfa7+/YCL8MNvcmi5Ah//cRIf/RERELoJh2Fuqvp9mH339prchYZjbVq/DfyIiIlI7WCxw7T8h6VYwrDBnPOz80uyqKqVQJSIiIp7NywuGvgIdbwRbGcy6HXYvMbuqChSqRERExPN5ecP1b0D7oWAtgY9uhb0rzK6qHIUqF1NHdRERESfx9oEb/gdtBkJZEXw4GvavMbsqB3VUdxN1VBcREXGS0iL46GbY8zX4hcLtC6BxN5esSh3VRUREpPbyDYBRH0CzPlByEt6/Hg5vMrsqhSoRERGpgfyC4OaPID4ZivLg3eGQvc3UkhSqREREpGbyD4FbZ0NcVzh9DOb93j6ulUkUqkRERKTmCgiHMR9Dm0Fw03T7uFYm8TFtzXVEWloaaWlpWK1Ws0sRERGpnQIj4ZaPzK5CZ/+5i87+ExERqXl09p+IiIiImylUiYiIiDiBQpWIiIiIEyhUiYiIiDiBQpWIiIiIEyhUiYiIiDiBQpWIiIiIEyhUuVhaWhoJCQl0797d7FJERETEhTT4p5to8E8REZGaR4N/ioiIiLiZQpWIiIiIEyhUiYiIiDiBj9kF1BVnuq7l5+ebXImIiIhU1Znf7ap0QVeocpOTJ08CEB8fb3IlIiIiUl0nT54kPDz8vPPo7D83sdlsZGZmEhoaisViKfdYfn4+8fHxHDx4UGcGVoO228XRdqs+bbOLo+12cbTdLo6rtpthGJw8eZK4uDi8vM7fa0otVW7i5eVF48aNzztPWFiYPkAXQdvt4mi7VZ+22cXRdrs42m4XxxXb7UItVGeoo7qIiIiIEyhUiYiIiDiBQpUH8Pf354knnsDf39/sUmoUbbeLo+1WfdpmF0fb7eJou10cT9hu6qguIiIi4gRqqRIRERFxAoUqERERESdQqBIRERFxAoUqERERESdQqBIRERFxAoUqN0lLS6NZs2YEBASQnJzMunXrzjv/7NmzadeuHQEBAXTq1IkvvvjCTZV6lupst+nTp2OxWMrdAgIC3Fit+VauXMmQIUOIi4vDYrEwf/78Cz5n+fLldO3aFX9/f1q1asX06dNdXqenqe52W758eYV9zWKxkJWV5Z6CPcCUKVPo3r07oaGhREdHM3z4cHbu3HnB59X177aL2W76boPXXnuNzp07O0ZL79WrF19++eV5n2PGvqZQ5QYzZ85k0qRJPPHEE2zYsIHExERSU1PJycmpdP7Vq1dz8803c+edd7Jx40aGDx/O8OHD2bJli5srN1d1txvYL09w+PBhx23//v1urNh8BQUFJCYmkpaWVqX59+3bx+DBg+nXrx8ZGRk89NBD/O53v2PhwoUurtSzVHe7nbFz585y+1t0dLSLKvQ8K1asYMKECXz33XcsXryY0tJSBgwYQEFBwTmfo++2i9tuoO+2xo0b8/e//5309HS+//57+vfvz7Bhw9i6dWul85u2rxnicj169DAmTJjguG+1Wo24uDhjypQplc4/cuRIY/DgweWmJScnG7///e9dWqenqe52e/vtt43w8HA3Vef5AGPevHnnnedPf/qT0aFDh3LTRo0aZaSmprqwMs9Wle22bNkyAzCOHz/ulppqgpycHAMwVqxYcc559N1WUVW2m77bKhcZGWm89dZblT5m1r6mlioXKykpIT09nZSUFMc0Ly8vUlJSWLNmTaXPWbNmTbn5AVJTU885f210MdsN4NSpUzRt2pT4+Pjz/i9G7LSvXZqkpCQaNmzINddcw6pVq8wux1R5eXkA1KtX75zzaH+rqCrbDfTddjar1cpHH31EQUEBvXr1qnQes/Y1hSoXO3LkCFarlZiYmHLTY2Jiztn/Iisrq1rz10YXs93atm3LtGnTWLBgAe+//z42m43evXvz888/u6PkGulc+1p+fj6nT582qSrP17BhQ15//XXmzp3L3LlziY+Pp2/fvmzYsMHs0kxhs9l46KGHuPzyy+nYseM559N3W3lV3W76brPbvHkzISEh+Pv7c8899zBv3jwSEhIqndesfc3HpUsXcaNevXqV+19L7969ad++PW+88QZPP/20iZVJbdO2bVvatm3ruN+7d2/27NnDv/71L9577z0TKzPHhAkT2LJlC99++63ZpdQoVd1u+m6za9u2LRkZGeTl5TFnzhzGjh3LihUrzhmszKCWKhdr0KAB3t7eZGdnl5uenZ1NbGxspc+JjY2t1vy10cVst9/y9fWlS5cu7N692xUl1grn2tfCwsIIDAw0qaqaqUePHnVyX5s4cSKfffYZy5Yto3HjxuedV99tv6rOdvutuvrd5ufnR6tWrejWrRtTpkwhMTGRf//735XOa9a+plDlYn5+fnTr1o2lS5c6ptlsNpYuXXrOY8G9evUqNz/A4sWLzzl/bXQx2+23rFYrmzdvpmHDhq4qs8bTvuY8GRkZdWpfMwyDiRMnMm/ePL7++muaN29+wedof7u47fZb+m6zs9lsFBcXV/qYafuaS7vBi2EYhvHRRx8Z/v7+xvTp041t27YZd999txEREWFkZWUZhmEYY8aMMSZPnuyYf9WqVYaPj4/x4osvGtu3bzeeeOIJw9fX19i8ebNZL8EU1d1uTz31lLFw4UJjz549Rnp6ujF69GgjICDA2Lp1q1kvwe1OnjxpbNy40di4caMBGFOnTjU2btxo7N+/3zAMw5g8ebIxZswYx/x79+41goKCjEceecTYvn27kZaWZnh7extfffWVWS/BFNXdbv/617+M+fPnGz/++KOxefNm48EHHzS8vLyMJUuWmPUS3O7ee+81wsPDjeXLlxuHDx923AoLCx3z6LutoovZbvpus38GV6xYYezbt8/44YcfjMmTJxsWi8VYtGiRYRies68pVLnJK6+8YjRp0sTw8/MzevToYXz33XeOx6666ipj7Nix5eafNWuW0aZNG8PPz8/o0KGD8fnnn7u5Ys9Qne320EMPOeaNiYkxrr32WmPDhg0mVG2eM6f6//Z2ZjuNHTvWuOqqqyo8JykpyfDz8zNatGhhvP32226v22zV3W7PP/+80bJlSyMgIMCoV6+e0bdvX+Prr782p3iTVLa9gHL7j77bKrqY7abvNsMYP3680bRpU8PPz8+Iiooyrr76akegMgzP2dcshmEYrm0LExEREan91KdKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAkUqkREREScQKFKRERExAn+HymrmVJVGBaWAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "mp = 0.938 # proton mass in GeV\n", "Tk = 0.1 # kinetic freeze-out temperature in GeV\n", "n = 1\n", "\n", "pt = np.linspace(0.1, 3., 100)\n", "\n", "# your code here\n", "plt.plot(pt,dndpt_blastwave(pt,mp,Tk,0.8,n),label = r\"$\\beta_S$ = 0.8\")\n", "plt.plot(pt,dndpt_blastwave(pt,mp,Tk,0.2,n),label = r\"$\\beta_S$ = 0.2\")\n", "\n", "# plt.plot(pt,dndpt_blastwave(pt,mp,Tk,0.,n),label = r\"$\\beta_S$ = 0\")\n", "plt.legend(loc='center right')\n", "plt.title(r\"$p_T$ spectra from blast-wave model for protons\")\n", "plt.ylabel(r\"$\\frac{1}{2 \\pi p_T}\\frac{dN}{dp_T dy}$\")\n", "plt.yscale('log')\n", "# plt.xlim((0.1,3))\n", "plt.show()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### b) Download the [ALICE data from hepdata](https://www.hepdata.net/download/submission/ins1759506/1/csv). Here is the corresponding [ALICE paper](https://arxiv.org/abs/1910.07678). Extract the data tables from the tar file, e.g. on the command line with ``tar xvf HEPData-ins1759506-v1-csv.tar``. Read the pion, kaon, and proton spectra using the code below:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# read pion, kaon, and proton pT spectra in 0-5% Pb-Pb collisions at 5.02 TeV\n", "# Hepdata https://www.hepdata.net/record/ins1759506\n", "pt = {}\n", "dndpt = {}\n", "dndpt_err = {}\n", "pt[\"pion\"], dndpt[\"pion\"], dndpt_err[\"pion\"] = np.loadtxt(\"HEPData-ins1759506-v1-csv/Table1.csv\", delimiter=',', usecols=(0,3,6), unpack=True, skiprows=13, max_rows=58)\n", "pt[\"kaon\"], dndpt[\"kaon\"], dndpt_err[\"kaon\"] = np.loadtxt(\"HEPData-ins1759506-v1-csv/Table3.csv\", delimiter=',', usecols=(0,3,6), unpack=True, skiprows=13, max_rows=53)\n", "pt[\"proton\"], dndpt[\"proton\"], dndpt_err[\"proton\"] = np.loadtxt(\"HEPData-ins1759506-v1-csv/Table5.csv\", delimiter=',', usecols=(0,3,6), unpack=True, skiprows=13, max_rows=51)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.11 , 0.13 , 0.15 , 0.17 , 0.19 , 0.225, 0.275, 0.325,\n", " 0.375, 0.425, 0.475, 0.525, 0.575, 0.625, 0.675, 0.725,\n", " 0.775, 0.825, 0.875, 0.925, 0.975, 1.05 , 1.15 , 1.25 ,\n", " 1.35 , 1.45 , 1.55 , 1.65 , 1.75 , 1.85 , 1.95 , 2.1 ,\n", " 2.3 , 2.5 , 2.7 , 2.9 , 3.1 , 3.3 , 3.5 , 3.7 ,\n", " 3.9 , 4.25 , 4.75 , 5.25 , 5.75 , 6.25 , 6.75 , 7.5 ,\n", " 8.5 , 9.5 , 10.5 , 11.5 , 12.5 , 13.5 , 14.5 , 15.5 ,\n", " 17. , 19. ])" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "pt[\"pion\"]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "# define fit ranges for each particle, define also particle masses\n", "pt_min = {\"pion\": 0.5, \"kaon\": 0.2, \"proton\": 0.3}\n", "pt_max = {\"pion\": 1.0, \"kaon\": 1.5, \"proton\": 3.0}\n", "mass = {\"pion\": 0.139, \"kaon\": 0.494, \"proton\": 0.938}" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "# prepare data in the selected pt ranges\n", "pt_meas = {}\n", "dndpt_meas = {}\n", "dndpt_meas_err = {}\n", "for particle in [\"pion\", \"kaon\", \"proton\"]:\n", " sel = (pt[particle] >= pt_min[particle]) & (pt[particle] <= pt_max[particle])\n", " pt_meas[particle] = pt[particle][sel]\n", " dndpt_meas[particle] = dndpt[particle][sel]\n", " dndpt_meas_err[particle] = dndpt_err[particle][sel]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We take the normalization of the blast-wave function as an arbitrary parameter that we adjust in the fit range for a given set of blast-wave parameters $T$, $\\beta_s$, and $n$. A simple least squares minimization gives the following formula for the normalization (factor in front of the blast-wave formula):\n", "\n", "$$\n", "A = \\frac{\\sum_i \\frac{y_i^\\mathrm{data} y_i^\\mathrm{BW}}{\\sigma_i^2}}{\\sum_i \\frac{ {y_i^\\mathrm{BW}}^2 }{\\sigma_i^2} }\n", "$$\n", "\n", "This is implemented in the function 'normalization' (see below)." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "# calculate normalization of the blastwave function for given blastwave parameters and particle species\n", "def normalization(particle, Tkin, beta_s, n):\n", " dndpt_bw = dndpt_blastwave(pt_meas[particle], mass[particle], Tkin, beta_s, n)\n", " err_square = dndpt_meas_err[particle]**2\n", " A_num = np.sum((dndpt_meas[particle] * dndpt_bw) / err_square)\n", " A_den = np.sum(dndpt_bw**2 / err_square)\n", " return A_num / A_den " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Next we implement a function that returns the $\\chi^2$ for a given particle species" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "# chi2 for a given particle species\n", "def chi2_particle(particle, Tkin, beta_s, n):\n", " \n", " # blast-wave prediction\n", " dndpt_bw = dndpt_blastwave(pt_meas[particle], mass[particle], Tkin, beta_s, n)\n", "\n", " # normalization\n", " A = normalization(particle, Tkin, beta_s, n)\n", "\n", " pulls = (dndpt_meas[particle] - A * dndpt_bw) / dndpt_meas_err[particle]\n", " return np.sum(pulls * pulls) " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "#### c) Implement the total $\\chi^2$ as the sum of the contributions from pions, kaons, and protons" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# total chi2 for all considered particles\n", "\n", "def chi2(Tkin, beta_s, n):\n", "\n", " # your code here\n", " return chi2_particle(\"pion\",Tkin,beta_s,n) + chi2_particle(\"kaon\",Tkin,beta_s,n) + chi2_particle(\"proton\",Tkin,beta_s,n)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we perform the $\\chi^2$ fit using iminuit" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Migrad
FCN = 79.55 Nfcn = 96
EDM = 2.94e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 Tkin 0.0896 0.0028 0 0.3
1 beta_s 0.9072 0.0025 0 1
2 n 0.735 0.013 0 3
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Tkin beta_s n
Tkin 7.65e-06 -5.38e-06 (-0.778) 3.34e-06 (0.092)
beta_s -5.38e-06 (-0.778) 6.24e-06 1.28e-05 (0.388)
n 3.34e-06 (0.092) 1.28e-05 (0.388) 0.000173
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 79.55 │ Nfcn = 96 │\n", "│ EDM = 2.94e-05 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ Tkin │ 0.0896 │ 0.0028 │ │ │ 0 │ 0.3 │ │\n", "│ 1 │ beta_s │ 0.9072 │ 0.0025 │ │ │ 0 │ 1 │ │\n", "│ 2 │ n │ 0.735 │ 0.013 │ │ │ 0 │ 3 │ │\n", "└───┴────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌────────┬───────────────────────────────┐\n", "│ │ Tkin beta_s n │\n", "├────────┼───────────────────────────────┤\n", "│ Tkin │ 7.65e-06 -5.38e-06 3.34e-06 │\n", "│ beta_s │ -5.38e-06 6.24e-06 1.28e-05 │\n", "│ n │ 3.34e-06 1.28e-05 0.000173 │\n", "└────────┴───────────────────────────────┘" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = Minuit(chi2, Tkin=0.1, beta_s=0.85, n=1.)\n", "m.errordef = Minuit.LEAST_SQUARES\n", "\n", "# define limits for Tkin, beta_s, and n\n", "m.limits = [(0, 0.3), (0, 1.), (0., 3.)]\n", "\n", "# perform the chi2 minimization and calculate the covariance matrix of the parameters\n", "m.migrad()\n", "m.hesse() " ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Now we plot the fit result along with the measured pion, kaon, and proton spectra" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAGlCAYAAADgRxw/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABomUlEQVR4nO3deVyU1f4H8M8wCigq7orOGOaCO5aKaZFYpllpRfwqK9P2DA1TW8y8Zl1Tq5u4oHYttTLMBFxupS02KJmlqZi7pahIuKGCK8vM9/fHYZBlBgaYYRY+79frCZl5nmfOPA0z3znne75HIyICIiIiIjfk5ewGEBEREVUUAxkiIiJyWwxkiIiIyG0xkCEiIiK3xUCGiIiI3BYDGSIiInJbDGSIiIjIbTGQISIiIrfFQIaIiIjcFgMZIiIiclsMZIiIiMhtOS2QuXLlCm644QZMmDDBWU0gIiIiN+e0QGbatGm45ZZbnPXwRERE5AGcEsj89ddfOHDgAAYPHuyMhyciIiIPUe5AZtOmTRgyZAhatGgBjUaD1atXl9gnJiYGgYGB8PX1Re/evbF169Yi90+YMAHTp0+vcKOJiIiIgAoEMpcvX0ZwcDBiYmIs3r9ixQqMGzcOU6ZMwY4dOxAcHIxBgwbh9OnTAIA1a9agffv2aN++feVaTkRERNWeRkSkwgdrNFi1ahUeeOCBgtt69+6NXr16Yd68eQAAk8kEvV6PMWPG4I033sDEiROxbNkyaLVaXLp0Cbm5uRg/fjz+9a9/WXyM7OxsZGdnF/xuMplw7tw5NGrUCBqNpqJNJyIioiokIrh48SJatGgBLy87ZrZIJQCQVatWFfyenZ0tWq22yG0iIk8++aQMHTq0xPFLliyR8ePHl/oYU6ZMEQDcuHHjxo0bNw/YUlNTKxN6lFADdnT27FkYjUY0a9asyO3NmjXDgQMHKnTOiRMnYty4cQW/Z2ZmolWrVkhNTUW9evUq1V4iIiKqGllZWdDr9ahbt65dz2vXQKa8Ro4cWeY+Pj4+8PHxKXF7vXr1GMgQERG5GXunhdh1+nXjxo2h1Wpx6tSpIrefOnUKzZs3t+dDEREREdk3kPH29kaPHj2wYcOGgttMJhM2bNiAPn36VOrcMTEx6NSpE3r16lXZZhIREZGHKPfQ0qVLl/D3338X/J6SkoLk5GQ0bNgQrVq1wrhx4zBixAj07NkTISEhiI6OxuXLl/HUU09VqqGRkZGIjIxEVlYW/P39K3UuIiIi8gzlDmT++OMP9O/fv+B3cyLuiBEjsHTpUjzyyCM4c+YM/vWvf+HkyZPo3r071q9fXyIBmIiIiKiyKlVHxhnMPTKZmZlM9iUiInITjvr8dtqikURERESV5TaBDJN9iYiIqDgOLREREZHDcWiJiIiIqBgGMkREROS2GMgQERGR23KbQIbJvkRERFQck32JiIjI4ZjsS0RERFQMAxkiIiJyWwxkiIiIyG0xkCEiIiK35TaBDGctERERUXGctUREREQOx1lLRERERMUwkCEiIiK3xUCGiIiI3BYDGSIiInJbDGSIiIjIbblNIMPp10RERFQcp18TERGRw3H6NREREVExDGSIiIjIbTGQISIiIrfFQIaIiIjcFgMZIiIiclsMZIiIiMhtMZAhIiIit+U2gQwL4hEREVFxLIhHREREDseCeERERETFMJAhIiIit8VAhoiIiNwWAxkiIiJyWwxkiIiIyG0xkCEiIiK3xUCGiIiI3BYDGSIiInJbDGSIiIjIbTGQISIiIrflNoEM11oiIiKi4rjWEhERETkc11oiIiIiKoaBDBEREbmtGs5ugLswGo1ISkpCeno6AgICEBoaCq1W6+xmERERVWsMZEphDl7WrFmDL7/8EmfOnCm4T6fTYfbs2bj//vsZ4BARETkJA5liSgteCktLS8NDDz2ERo0aISMjo+B2c4ATHh5e7sdkMERERFQ+nLVUSEJCAqKionDixIkKn0Oj0QAA4uLiCoKZ0gIVS49ZkWCIiIjIlTnq85uBDFSgMW3aNEyZMsUu59NoNNDpdEhJScGaNWusBioAEBERgeL/CywFQ0RERO6MgUw+e1+IhIQEvPzyy0hLS7ND64oaN24cZs2aZTFQEZESw1LF9zEHQwA49ERERG6NgUw+e1yIwnkw0dHR9m2gnU2dOhWLFi3i0BMREbk1BjL5Knsh7JEH42wceiIiInfDyr52kJCQgIiICIcHMRqNBg0bNnTY+c2x59ixY2E0Gh32OERERK6u2gQyRqMRUVFRJfJVKqpRo0YArveOmJl/j4qKssvjWCMiSE1NRVJSEoxGIxITE7F8+XIkJiYyuCEiomqj2gQySUlJle6JadKkCcaOHQuDwYBTp04hPj4eLVu2LLKPTqdDXFwcJk2aBJ1OVyLQMdNoNGjUqBE0Go3VfWyxZs0aBAYGon///njsscfQv39/BAYGIiEhocLnJCIichfVJkdm+fLleOyxx8r9eA0aNMCIESNw//33W5wtVFaNmIiICAAo0hNUOMcFQImcHb1ej2eeeQZvv/12udtb/PzMoSEiIldQ7ZN9Y2JiEBMTA6PRiEOHDpX7QiQmJqJ///7lesypU6di0qRJlZrqbCm5WK/XIzo6utSCeQAQGBiItLS0Cg2HFZ6+zanaRETkbNU+kDGr6IUwGo02BwaNGzfGxx9/bLfejIouQWCtR6c8DAYDwsLCuAwCERE5FWctVZJWqy2opmstJ8Xb2xtz5szByZMn7Toko9VqERYWhmHDhiEsLMzmACI8PBxxcXEl8nD0ej3Gjh1r0znS09ORkJDAPBoiIvJI1aZHxsxaHZl69eph48aN6N69u51aaj+WelOSkpJsGiobMmQIvvnmGy6DQERETsWhpXz2rOw7d+5cJCQkoHHjxti8eTPat29v59Y6TnmGyqxhHg0REVUVDi3ZkVarRZMmTfC///0PADBv3jy3CmKA0ofKzFO6Bw0aVOo5WIuGiIjcXbUMZIxGI5555hnk5uZiyJAhePjhh53dpAqxlkNjrmUzYsQIm87DWjREROSuquXQ0uzZszF27FjUrVsX+/btg06ns3Mrq5a1GUkVmXJuxhwaIiKyJ+bI5KvshTh69Cg6d+6MK1euYMGCBXjxxRcd0ErXUNk8GubQEBGRvTBHxg5EBC+88AKuXLmC22+/Hc8//7yzm+RQZeXRlKVwDg0REZErqlaBzBdffIEffvgBPj4+WLRoEby8PP/pl5ZHU55aNEwGJiIiV1RthpZycnKg1+tx+vRpTJ8+HW+88YYDW+l6KlOLZvTo0Vi9enWR2js6nQ6zZ89m/gwREdmEOTL5KnohVq9ejQcffBDNmzfH8ePHUbNmTQe20j1UJoeGycBERFQezJGppCVLlgAAhg8fziAmny05NN7e3haPNQc+Y8eO5TATERE5TbUIZE6dOoVvv/0WAPDUU085uTWupbQcmqlTpyInJ8fqsUwGJiIiZ6vh7AZUhWXLlsFoNKJ3797o2LGjs5vjcsLDw3H//feXyKH5+uuvbTr+6NGjXF2biIicwuMDGREpGFZib4x15hW6CwsICLDp2DFjxmDcuHE4f/58wW1MBiYioqrg8UNL27dvx969e+Hr64tHHnnE2c1xK6GhodDpdKXWnPHy8sKlS5eKBDEAkJaWhoiICC5zQEREDuXxgYy5N+bBBx9E/fr1ndsYN2NLMrC1a8pkYCIiqgoeHchcu3YNsbGxADisVFFlJQOfO3fO6rHFk4FZVI+IiOzNo3Nk1qxZgwsXLkCv1+OOO+5wdnPcVmWTgWfMmIGjR49i8uTJLKpHRER25dGBjHlYacSIEZxBU0mVSQb+/vvv8f3335e43ZxHw6J6RERUUVU+tHThwgX07NkT3bt3R5cuXbBo0SKHPE5aWhp+/PFHAMDIkSMd8hjVXVnJwBqNBo0bN7ZagJB5NEREVFlVHsjUrVsXmzZtQnJyMn7//Xe89957yMjIsPvjfP755zCZTAgNDUWbNm3sfn6yLRl4zJgxyM3NtXqOwnk0zKEhIqLyqvJARqvVonbt2gCA7OxsiEi51/mxxeeffw6ASb6OVloycFxcHNq1a2fTedasWYPAwED0798fjz32GPr374/AwEBO3yYiolKVO5DZtGkThgwZghYtWkCj0WD16tUl9omJiUFgYCB8fX3Ru3dvbN26tcj9Fy5cQHBwMHQ6HV599VU0bty4wk/AkoyMDBw4cACAmnZNjhUeHo6jR4/CYDAgNjYWBoMBKSkpCA8PtzmPJjo6ukgiMMBaNEREVLZyBzKXL19GcHAwYmJiLN6/YsUKjBs3DlOmTMGOHTsQHByMQYMG4fTp0wX71K9fH7t27UJKSgpiY2Nx6tSpij8DC3bu3AkAaNu2LWvHVBFzMvCwYcMQFhZWkFxtS1E9a5hDQ0REZSl3IDN48GD8+9//ttrT8dFHH+G5557DU089hU6dOmHhwoWoXbs2Fi9eXGLfZs2aITg4uNRFB7Ozs5GVlVVkK8uOHTsAADfddJONz4ocxZY8mtJwYUoiIiqNXXNkcnJysH37dgwYMOD6A3h5YcCAAdiyZQsAtRL1xYsXAQCZmZnYtGkTgoKCrJ5z+vTp8Pf3L9j0en2Z7TD3yNx8882VeTpkJ6Xl0YwdO9amc6SnpzMZmIiISrBrIHP27FkYjUY0a9asyO3NmjXDyZMnAQDHjh1DaGgogoODERoaijFjxqBr165Wzzlx4kRkZmYWbKmpqWW2w9wjw0DGdVjLo7n//vttOn7v3r1MBiYiohKqvCBeSEgIkpOTbd7fx8cHPj4+Nu9/8eJFHDp0CACHllyNpaJ65hyatLS0UmevTZs2rcRtLKhHRER27ZFp3LgxtFptieTdU6dOoXnz5vZ8KKt27doFQA1bNGnSpEoekyrOlhwaLy/LL1MmAxMRkV0DGW9vb/To0QMbNmwouM1kMmHDhg3o06dPpc4dExODTp06oVevXqXux0Rf91PWwpQmk8nqsSyoR0RUvZV7aOnSpUv4+++/C35PSUlBcnIyGjZsiFatWmHcuHEYMWIEevbsiZCQEERHR+Py5cuVLkwXGRmJyMhIZGVlwd/f3+p+zI9xT5VdmHLNmjUYPnw4F6UkIqpmyh3I/PHHH+jfv3/B7+PGjQOgFmZcunQpHnnkEZw5cwb/+te/cPLkSXTv3h3r168vkQDsKJyx5L4qszBldHR0iduYQ0NE5Pk04oj1ARzI3COTmZmJevXqFbnv2rVrqFOnDoxGI44fP27TVG1ybUajEYGBgWUmA1uj0Wig0+mQkpLCFdCJiJyotM/vyqjytZYcaffu3TAajWjcuDF0Op2zm0N2YO+CesyjISLyLG4TyNiS7Ft4WKkiJfHJNdmroF5CQgJr0RAReRi3CWQiIyOxb98+bNu2zeo+nLHkuSpbUO+vv/5CREQEF6YkIvIwHpUjExISgm3btmHFihV4+OGHndRCqkq25tD4+fnh8uXLFu9jHg0RkeMxR6YMubm5+PPPPwGwR6Y6sTWHxloQA3BhSiIid+YxgcyBAweQnZ2NunXrok2bNs5uDlWh0nJo4uPjMXnyZJvOw4UpiYjcT5WvtVRRMTExiImJsfrBUjg/xlpJe/Jc1grqabVaNGzYEO+++26Z5/jrr78QGBjIonpERG7EY3Jkxo4di9mzZyMqKspicTSqvmzJo2nUqBHOnTtX4n7z8BSL6hERVQ5zZMrApQnImtLyaMzOnz9vMcjhwpRERK7NIwIZk8nEpQmoVNbyaBo2bIhatWpxYUoiIjflEYHM4cOHcenSJfj6+qJDhw7Obg65KEu1aE6fPo05c+bYdPyaNWtYUI+IyMW4TbJvaczDSt26dUONGh7xlMhBLC1M2bZtW5uO5cKURESuxyN6ZMzDSqwfQxURGhoKnU5X6rIW1mbCMYeGiMi53CaQKW2tJSb6UmXYkgxsaw4NERFVLbcJZEpba2nXrl0A2CNDFWctGbh58+Zo0aKFTedIT093RNOIiKgUbhPIlObSpUsAgMaNGzu5JeTOLCUDnzhxAsuWLbPp+ICAAAe3kIiIimNmLFEhlpKBb7/9duh0OqsF9cyLToaGhlo8p9FotFhxmIiIKs8jemSIHKmsHBoRwciRIy0GJwkJCZyyTUTkQAxkiGxgLYfGPJvp3XffRUREBFJTUwvuS0hIQERERJG1m4DrU7YZzBARVR4DGSIbWcqhOXPmDKKioqDVahEfH48OHTpg5syZuHr1KqKiorjsARGRg7nNopGFV78+dOhQkUWn/Pz8cOXKFRw5cgStW7d2ckupOvrzzz8RGRmJX375BQCg1+uL9M5YYzAYSuTkEBF5omq/aGRp06+JnK1bt27YtGkTPv/8czRt2tSmIAbglG0iospym0CGyNVpNBoMHz4cBw8exIMPPmjTMZyyTURUOQxkiOysfv36WLlyJZo1a2Z1H41GA71eb3XKNhER2YaBDJEDaLVazJ8/v9T1m6Kjo1lPhoiokhjIEDmIecq2TqcrcZ+vry9SU1ORl5fnhJYREXkOt5m1ZGYp65mzlsiVFa7sm5mZiU8++QTbt28HAAQHB2P+/Pno27evTcezMjARuatqP2uJyF2Zlz0YNmwYXnzxRfz+++9YuHAhGjRogF27duHWW2/F008/jTNnzpQ4lpWBiYhKx0CGqIpptVq88MILOHjwIJ5++mkAwJIlSxAUFISPP/64oEgeKwMTEZXNbQKZmJgYdOrUCb169XJ2U4jsokmTJvj000+xefNmBAcH4/z583jxxRfRp08fbN26lZWBiYhswBwZIheQl5eH+fPnY/LkycjKyrL5OFYGJiJ3wRwZIg9Wo0YNvPzyyzhw4AAef/xxm49jZWAiqu4YyBC5kICAACxbtgyzZs2yeX8iouqMgQyRCxozZgxatmxp9X5WBiYiUhjIELkgrVaLOXPmQKPRWKwOLCKYNWuWzfVkjEYjEhMTsXz5ciQmJjJJmIg8BgMZIhdlrgxsrWfm448/xqFDh8o8D2vREJEn46wlIhdXuLJvo0aNsHnzZsycORPZ2dnw9vbGa6+9hjfffBO1atUqcay5Fk3xP3NzL09cXBzCw8Or5HkQUfXmqFlLDGSI3NDhw4cxZswYrFu3DgDQunVrzJkzB/fdd1/BPkajEYGBgSUK6plpNBrodDqkpKRwyQMicjhOvyaiAm3atMG3336LhIQE6PV6pKSkYMiQIXjggQdw7NgxAEBSUpLVIAZQeTapqalISkqqqmYTEdkdAxkiN6XRaPDggw9i//79eP3111GjRg2sWbMGHTt2xPTp05GammrTeViLhojcmdsEMlyigMgyPz8/zJgxA7t27UK/fv1w9epVvPnmm5g0aZJNx7MWDRG5M+bIEHkQEcGXX36J8ePH4/Tp06XuyxwZIqpKzJEhojJpNBo88cQTOHjwICIjI0vdDwCio6MZxBCRW2MgQ+SB6tevj3nz5uGPP/5AmzZtStyv0+nKNfWaBfWIyFUxkCEAgNEIJCYCy5ern/yc8gw9evTAwYMHMX/+fPj5+RXcftddd6Ffv342nYMF9YjIlTGQISQkAIGBQP/+wGOPqZ+Bgep2MwY67kur1WLUqFE4cuQIRo4cCQBYvHgxgoKC8Omnn8JkMlk91lxQr/g07rS0NERERDCYISKnY7JvNWE0AklJQHo6EBAAhIYCWq0KViIigOKvAvPyPnFx6mdUFFD4s0ynA2bPBswjE9bOT64nKSkJL730Evbs2QMA6NOnDxYsWIDg4OAi+7GgHhHZk6OSfSFuJjMzUwBIZmZmwW21a9cWAHLkyBEntsx1xceL6HQiKlxRm04n8vXXJW8vvGk0Io0aqZ+W7tNo1LmtnT8+3tnPnKzJycmRDz/8UPz8/ASAaLVaGTt2bJG/K4PBIADK3AwGg/OeCBG5DUuf3/bAoSUPZ+5xKf6lOi0NePjhkrcXJgJkZJTsrTHfBwDPP2/9/BERHJ5yVTVr1sT48eNx4MABREREwGg0Ijo6Gh06dMCKFSsgIjYXymNBPSJyJgYyHsBagGA0qiGh0gKRyrAl0Bk7VrWDeTiuSafTYeXKlVi/fj3atm2L9PR0PProoxg0aBDy8vJsOgcL6hGRMzFHxs0lJFjPX2nYUAUMzjZ1KvD228zDcXXXrl3D+++/j/feew/Z2dmoWbMmfH19cenSpRKrZwPMkSGi8uHq1/mqYyBT0UTdqCggOrrKm1tu/v5AVhYDHVdx+PBhjB49GuvXr7e6j7mgXnlq0RBR9cZAJl91C2Ss9bh89BEwbpz1HBeNBqhXD8jMrNzja7XOH+Jp2BA4f77ygQ7ZTkSQkJCAqKgopKWllbhfr9cjOjq6XEGM0WhEUlIS0tPTERAQgNDQUPbkEFUjnLWUrzrNWoqPtz5jyNpMo/Ju9euLTJgg0rRp0dv1+uszkswzlCy1wdqsJvM+DRrYr62VmVlllpcnYjCIxMaqn3l5zvq/6x6ysrJkwoQJ4uXlJQDE29tbnn/+ebly5Uq5zhMfHy86na7IbCedTifxnNpGVG04atYSAxkXYOnDNS+v9KnRtm7t21sPAAp/yJf2AW9perUtgY5GIzJ1qmMDGVsCHb1ePR9OE6+4P//8U2677baCIKRjx46SmJho07Hx8fGi0WiKBDEARKPRiEajYTBDVE0wkMnnaYGMtQ9XewUABkPpgYitKhromAMye/YiVWSbOtW2Xhv22FhnMplk6dKl0rhx44JgZPjw4XLy5Emrx+Tl5ZXoiSkezOj1esnjhSbyeAxk8nlSIOPIoaPCPREijv+ALivQqejwlL02rbbsa7VyJXtsbJGRkSEvvvhiQS+Lv7+/xMTEWAxGWFSPiMyqfUG8mJgYdOrUCb169XJ2U+zCnjVezEmvxX+Pjr4+a0erBcLCgGHD1E9751iWdv7wcJWQ27Jl0WN0OiA+Hvjvf0t/Ho0albyvvEpLWBYBUlOB//s/2wr7VXcNGzbEggUL8Ntvv+Hmm29GZmYmIiMj0bt3b/zxxx9F9mVRPSJyNLcJZCIjI7Fv3z5s27bN2U2xi6Sk0qvqlkWjAfR6YOVKywFCXJxrzdYJDweOHgUMBiA2Vv1MSVG3OzLQ0WiA+vUr3m5zUGku7MeifdeFhIRg69atmDdvHurVq4ft27cjJCQEkZGROH/+PADbi+WxqB4RVZhd+3eqgDsOLVkadomNLd8wUXXI63D1hOOpUzn0ZE16ero8/vjjBUNFTZs2lc8//1xyc3NFp9NZTPYFc2SIqhXmyORzt0Cmssm8lj48y5uo6ylcMeHY0hTv6uznn3+WDh06FAQq/fr1k+jo6IIZSsWDGM5aIqo+HBXIsCCeA5VWeVcE8PMDLl+2fKxGo4ZWUlLU76xYW7bSKvua/18ARf9/mP9fVEbh/1f8/wLk5OTgo48+wjvvvIOrV6+iRo0auO+++7Bt27YixfXKW1SPBfWI3Bsr++Zzl0DGaFSLIlYkD6ZwxVpXynNxd5aqJOv1wH/+o6okp6VVLqj57DPgySdL3l5dl044duwYXn75ZaxduxaAWqDyhRdewI033ogWLVqUKxAxVxk+Ueh/nk6nw+zZs7lEApGbYGXffO4ytGQw2DY0cffdIi1bcuioqlgbnipriritW/v2qlLypk0iubkswicisnbtWgkMDCwYUrrnnnvk8OHDNh/PgnpEnoFDS/lcrUfG2rft5cuBxx4r+/jYWODhh6vnN3ZXY63H5tlngSlTyj6++LpUdeoAly6V3K869rhduXIF7733Ht5//33k5ubC19cXkyZNwquvvgofHx+rxxmNRgQGBhbpiSmMK3ATuQ8OLeVzpUDG2oKOs2erhQ779y/7HAaDqrtCrsFSYAqoYUJrQ0/mHJldu4CffgLWrgW++Qa4cMH641TXvJoDBw4gMjISP//8MwCgffv2iImJwYABAyzun5iYiP42/CEZDAaE8Q+JyKU5KpBxmzoyrsacPGqtgNqZM0Dz5taPN9eBMX9QkmuwVNhPq1XBKVB68cEGDVRRvS++uL4itzXmInzffmvnJ+DiOnTogJ9++gmxsbFo3rw5Dh06hLvuugvDhg2zWBSPBfWIqCwMZCrAlqq8o0YB+TXBSrBUeZdcW2lF+ywNEZ0+bdt5H3wQGDQI+OQTICPDPm11dRqNBsOGDcOBAwcwZswYeHl54auvvkJQUBDmzJmDvLy8gn1ZUI+IysKhpQpITLRt2AgAOncGzp1TQxVmer0KYqpLfoQnsXUGUnleI2Y1agADBgCPPALcf7/q4SnPY7qrHTt2YNSoUdi6dSsAoHv37liwYAFuueWWghyZtLQ0WHqrYo4Mkftgjkw+VwhkbE3kHTpUldjXaDz7g4hKMk+/Lyuv5ocfgFWrgK+/BpKTr99fsyZw991AmzZqGYpC5VcK8rA8KRA2mUz45JNP8MYbbxQsb/Dss89ixowZ2LhxIyLyiwAVfrvS5HdtxsXFcQo2kRtgIJPPFQIZW79tM5G3eiutCB9Qckjq4EEVtHz9NbB7t/XzevKspzNnzuD111/HkiVLAACNGjXC+++/j3r16uGVV14pMnupvAX1ABbVI3ImBjL5XCGQsfXbdnWbkUIlWZvSXdbQ4u7dwK23AhcvWt9Hr/fc19gvv/yCUaNGYc+ePQCAvn37Yt68ecjMzKxwEMKiekTOxVlLLqTwLJbimMhLhZW26ndpMjJKD2IANetp+HBg3z67Nddl3HbbbdixYwc+/PBD+Pn54ddff0WvXr2wdu1a3HfffQgLCyt3EBMREVGiHk1aWhoiIiKQkJBg76dARFWEgUwFhYcDb79dcjqutVksVH1ZmtJdFltnEy9frhLKe/ZUwbWts6XcQc2aNTF+/HgcOHAAERERMBqNmDVrFjp06ICVK1daTP61xGg0IioqyuL+5tvGjh0LY+FqhkTkNhjIVNDvvwMzZ6qhpe7dgU8/tf3bNlFZbJ1N3Levmu20fTswdizQogVw330q1+baNYc2scrodDqsXLkS69atQ5s2bfDPP//g4Ycfxt13342//vqrzOOTkpKsVgYGVDCTmpqKpKQkezabiKoIA5kKOHgQuPde4MoVYOBAYMsW4Omnbf+2TVSW0FDVu1e8x8/MXFBx0ybVezNvHhASovK3vv1WLXvRogUQGQls3Qrk5akk9eXL1U937Hy4++67sWfPHkyZMgU+Pj744Ycf0KVLF0yZMgVXr161ehyL6hF5NgYy5fTPP6qAWUYG0KuXml7t6+vsVpGnsbWasFYLNG6sApbffwcOHADefFMFQefPA/PnA717A7VqqZl2jz2mfgYGqkRkd+Pr64u3334bu3fvxqBBg5CTk4N33nkHXbt2xfr16y0ew6J6RJ6NgUw5XLiganscOwa0a6e++dap4+xWkacqbzVhAAgKAqZNUwnGP/54fQmMQsVyAahZVBER7hnMAEC7du2wbt06rFy5Ei1btsThw4cxePBgREREIDU1tci+oaGh0Ol0BXVnitNoNNDr9QjleiFEbqnKA5nU1FSEhYWhU6dO6NatG1auXFnVTaiQa9dUtdXdu9UaSt9/DzRp4uxWkaer6KwnrVb1vKSkWN9HBBg5svSaNa5Mo9EgIiIC+/fvx/jx46HVahEfH4+OHTviww8/RG5uLgBAq9Vidn73VvFgxvx7dHQ068kQuSupYv/884/s3LlTRETS09OlRYsWcunSJZuPz8zMFACSmZlZcFvt2rUFgKSkpNi5tYrJJPLwwyKASN26IvnNJ3JpBoN6zdqy9e0r8umnIhcvOrvVFbdr1y7p27evABAA0rlzZ9m0aVPB/fHx8aLT6QruByB6vV7i4+Ntfoy8vDwxGAwSGxsrBoNB8vLyHPFUiDySpc9ve6jyQKa4bt26yfHjx23e3xmBzNy56s2+Zk2Rn392yEMQ2V1srG1BjEZz/d916og8/7zItm0qgHc3RqNRFi9eLI0aNSoIVkaMGCGnTp0SkcoFIpYCIZ1OV65AiKg6c5lAZuPGjXLfffdJQECAAJBVq1aV2GfevHlyww03iI+Pj4SEhMjvv/9u8Vx//PGHdO7cuVyPX9WBzI4dIt7e6k1+1iy7n57IYWztkYmLE5kxQ6Rt26K3d+8uMn++yIULInl56nyxseqnq3dEnD17Vp577rmCgKN+/fqyYMGCCvegxMfHi0ajKRLEABCNRiMajYbBDJENXCaQ+e6772TSpEmSkJBgMZD56quvxNvbWxYvXix79+6V5557TurXr1/wjcgsIyNDOnXqJJs3by7X41dlIJOVJdKunXpTHzLEPb+hUvWVlyei0xXtcSneE6PXXw9KTCaRxESRxx8X8fG5vp+3t0jt2kWP1elE3OGze8uWLdK9e/eCwKNXr16yffv2cp0jLy+vRE9M8WBGr9dzmImoDC4TyBQ52EIgExISIpGRkQW/G41GadGihUyfPr3gtmvXrkloaKh8/vnnZT7GtWvXJDMzs2BLTU2tskBm+PDrb9pnz9r11ERVIj5eBSzFgxnzbdaCkYwMkeho9dq3FgSVdrwryc3NldmzZ0vdunUFgHh5ecno0aPl/PnzNh1vMBisBjGFN4PB4NDnQeTuHBXI2HXWUk5ODrZv344BAwYU3Obl5YUBAwZgy5YtAFQVzZEjR+KOO+7A8OHDyzzn9OnT4e/vX7Dp9Xp7Ntmqzz4DvvgC8PJSs0UaNaqShyWyq4pM4QaAhg2B0aOtn9cc0rz0kusX16tRowZefvllHDx4EMOGDYPJZMK8efPQoUMHfPnll2UudcCCekSuza6BzNmzZ2E0GtGsWbMitzdr1gwnT54EAGzevBkrVqzA6tWr0b17d3Tv3h27S5n/OXHiRGRmZhZsxWtEOMKBA+oNGgCmTr1ei4PIHVV0CndSUtFVuy05dQro0gVYskRVunZlAQEBiI2NxU8//YSgoCCcOnUKTzzxBO68807s37+/1ONsPT8RVb0qryNz2223wWQyITk5uWDr2rWr1f19fHxQr169IpsjGY2q+umVK8AddwATJzr04YiqhCMXrjxwQC3R0bIlMG4ccOhQZVrqeHfeeSd27dqFadOmwdfXFwaDAcHBwXjzzTdxxUI0xoJ6RK7NroFM48aNodVqcerUqSK3nzp1Cs2bN7fnQznMf/8L7NwJNGgALFvGtZOo+rK1g+G559SSBxcuALNmqerCAwcCq1eXrCjsKnx8fPDmm29i3759uO+++5Cbm4vp06ejU6dOWLt2bZF9WVCPyLXZNZDx9vZGjx49sGHDhoLbTCYTNmzYgD59+tjzoRwiIwN46y3173fftf2NnMgT2bpw5YIFwOHDasmOe+9Vt//4I/Dgg0Dr1mrJhNOnVW+nqy1c2bp1a/zvf//DmjVr0KpVKxw7dgz3338/hg4diqNHjxbsFx4ejri4OLQslmyk0+kQFxeH8HIseW80GpGYmIjly5cjMTERRle4EETurLzZwRcvXpSdO3fKzp07BYB89NFHsnPnTjl27JiIqOnXPj4+snTpUtm3b588//zzUr9+fTl58mSlspLnzZsnHTt2lPbt2zts1tJLL6kUxq5dRXJzK3UqIo9QkVlPR46IvP66SKNG1/evUcP1p3BfunRJJk6cKDVr1hQAUqtWLZk2bZpcu3atYJ/KVvZlUT2qzlxm+rW1qYgjRowo2Gfu3LnSqlUr8fb2lpCQEPntt9/s1mBH1ZFJThbx8lJvsJxFSXRdfHzJadh6fdlByNWrIp99VrLQnqtP4d63b5+EhYUVvLcFBQXJTz/9VOnzsqgeVXeOCmQ0ImXMPXQxWVlZ8Pf3R2ZmZkHir5+fH65cuYKUlBQEBgaW+5wiQL9+apbGww8DK1bYudFEbs5oVH8f6elqyDU01Lb8MaNR5c+UNvspIABITXWtfDQRQWxsLMaPH1+Q8/foo4/io48+qtDsJKPRiMDAQJywciE0Gg10Oh1SUlKYa0Mey9Lntz1U+awlV7RihXqTrlUL+PBDZ7eGyPVUZNYTYNsU7vR0NUNw0yb1pcIVaDQaPP744zhw4ABGjx4NLy8vfPXVVwgKCsKcOXOQV84s5qSkJKtBDKACp9TUVCQlJVW26UTVjtsEMjExMejUqRN69epl1/NevgxMmKD+/eabKnmRiOzD1incmzapXtGbb1Y1aa5dc2y7bFW/fn3MnTsXW7duRa9evXDx4kVERUWhV69e+O2332w+D4vqETmO2wQykZGR2LdvH7Zt22bX806fDqSlqdkV5oCGiOzD1lGY++5TPaLJyaomjV4PTJqk/jZdQY8ePbBlyxYsXLgQ9evXR3JyMvr06YPnn38eGRkZZR7PonpEjuM2gYwjZGQAH32k/v2f/wC+vs5tD5GnsXUK9+rVaghq5kygVSvg7FngvfdUfs1jjwFbtzp/+rZWq8ULL7yAgwcPYuTIkQCARYsWoUOHDliyZAlMJpPVY1lUj8hxqnUg8/HHwNWrqjv7gQec3Roiz6PVAvm15EoEM+bfo6PVfg0bAq+9pmrSxMWpICgvTwUuvXsDtWsD/furwKZ/fxXkJCRU5bNRmjZtiiVLlmDTpk3o0qULzp49i6effhq33347/vzzT4vHsKgekeNU20AmOxuYO1f9e9w4698YiahyyrtwZY0awEMPqbyZ7dtV7gwA5OQU3S8tDYiIcE4wA6helh07duDDDz+En58fNm/ejJtvvhnjxo3DxYsXS+xvr6J6LKhHVFS1nX792WfAyJHqzfXIEcDb23FtJqKKTeF2l+nbJ06cwCuvvIK4uDgAQIsWLRAdHY2IiIgSPTBGoxFJSUlIT09HQEAAQkNDbe6JSUhIQFRUVJEZUDqdDrNnzy5XdWEiZ3DU9OtyF8RzltIq+9aqVatcBfFMJpFu3VRRrhkzHNRgIqo0g8FyMb3iW58+al+TybntXbdunbRp06ag2N1dd90lBw8etMu5WVCP3B0L4uWzFNHVrl0bV69etblHZsMGYMAANeZ+4oRaIJKIXM/y5SonxlY33aSGih9+2Hm9rNeuXcPMmTMxffp0ZGdnw9vbG6+//jomTpyIWrVqVeicLKhHnoAF8ezIPFPp6acZxBC5MltnIw8dqqZv79wJDB+uyinMmAGcP+/Y9lni6+uLKVOmYM+ePRg0aBBycnLw7rvvonPnzvjuu+8qdE4W1COyrtoFMvv3A999p5J7x451dmuIqDS2Tt9OSFB5MtOmAc2bA//8A0ycqO6LigJSUtT+VTmFu23btli3bh1WrlyJli1bIiUlBffeey/Cw8Nx/Pjxcp2LBfWIrKt2gUx0tPr5wANAmzbObAkRlaU807cbNVLVuY8eBZYuBbp2VZW758wB2rYF+vYFWrSo2incGo0GERER2L9/PyZMmACtVotVq1ahY8eOmDlzJnKKT8WyggX1iKyrVjkyZ86oYlvXrqnZE7fdVkWNJqJKSUhQPSuFR1f0ehXEWJusIwL89JNaP+2HHyzvYw6GLE0Dd4Q9e/Zg1KhR+OWXXwAAHTt2xPz58xEWFlbqceYcmbS0NFh6y2aODLmDap8jY4+1lhYsUEFMr17ArbfasXFE5FDh4aqnxWAAYmPVz5SU0oMPjQa46y41lNysmeV9zDHB2LFVUym4S5cu2LRpE5YuXYomTZpg//796N+/P4YPH46TJ09aPY4F9YhKYdc5UFXA0vQtW6Zf5+WJBASoqZrLl1dBQ4nIJdg6hXv16qpt17lz52TUqFEFU6r9/f1l3rx5kpeXZ/WY+Ph40el0RaZf6/V6Tr0mt8Dp1/kqOrT066+qF8bfHzh9mgXwiKoLW6dwe3sDL7wAvPKKmvVUVbZt24ZRo0Zh+/btAICbb74ZCxYsQEhIiMX9K1NQj8iZqv3QUmWtXq1+3ncfgxii6sTW/NecHLVsSdu2wLBhwI4djm2XWa9evfD7779j3rx58Pf3x44dO3DLLbdg1KhROHfuXIn9tVotwsLCMGzYMISFhTGIoWqvWgQyIsCqVerfXBySqHqxdQr3+vXAwIGAyQR89RXQo4fKsfnxx+u5NI6avq3VahEZGYmDBw9i+PDhEBEsXLgQQUFBWLp0qcUEXyJSqkUgs38/8PffgI8PMGiQs1tDRFXJ1incgwYB338PJCcDjz+ujvvpJxXc3HyzqhgcGOjY6dvNmjXD559/jsTERHTq1Alnz57FU089hdtvvx27d++23wMReZBqEciYh5XuvBOoW9epTSEiJyjPCtzBwcCyZcDhw8DLL6ulTJKTgVmzSi5e6agVuPv164edO3di5syZqF27Nn755RfcdNNNmDBhgsWVte2NK2yTO3GbZN+YmBjExMTAaDTi0KFD5Ur2DQkBtm0D/vtf4LnnqrjhROQyKrIC9+nTQLt2QFaW5fs1GhUQpaQ4ZgXu48eP45VXXkFCfrTUsmVLzJo1y+LK2vbAFbbJURyV7Os2gYxZeWctnTihxr81GlW2vHlzJzSaiNxWYqIaRiqLwQCUUdeuUr777juMGTMGR44cAQAMHDgQ8+bNQ7t27ez2GAkJCYiIiCiRk2MOmOLi4hjMUIVx1lIFrV2rfvbpwyCGiMrP1uWLpk4FDh50XDvuuece7NmzB5MnT4a3tzd++OEHdOnSBVOmTMHVq1crfX6j0YioqCiLicXm28aOHcthJnI5Hh/ImPNjOFuJiCrC1unbiYlAx44qZya/JIzd1apVC++88w727NmDgQMHIicnB++88w66dOlS4ZW1zbjCNrkrjw5kLlxQ3b0AAxkiqhhbpm83bapqVIkA8fFAz55qtpPBcH3qtj21a9cO69evL1hZ+8iRIxVeWduMK2yTu/LoQGbdOiAvD+jUSSXrEbklexQvqew5HFVAxQ3YMn17wQLgf/8Ddu8GnnhCHfPjj8Add6hh7TVrVH0awH6XsvDK2uPHjy+ysvb7779v88raZlxhm9yWXRc8qALlWWvp4YfVGipvvlnFjSTHuHTp+sI4ly45pw15eWrxnthY9bOUdXEqdYxZfLyITld0USCdTt1eVeewRxs8gKXLoNdbvgxHjoi89JKIj8/1fTt3FomKctyl/PPPP+W2224rWIOpU6dOkpiYaPPxeXl5otPpCtZ+Kr5pNBrR6/WlrgVFVBpHrbXksYHMtWsideqoN4qtW53QULK/qghkSgs64uJEWrYs+ikUECASEyNy4IDI/v0i+/aJ7Nkjsnu32qKjRZo3L3nMJ5+IZGSIZGWpF6vJVLIt8fEiGk3J1Q01GrXZ8ulX2XPYow22XFs3Ud6ncPKkyBtviNSrZ32xyopcSmtMJpMsXbpUmjRpUhCAPPHEE3Ly5Embjo+PjxeNRlMimDHfxsUpqTK4aGQ+W6dfr1sH3HMP0KIFkJoKeHn0IFo1kZWlVv0EgO++U0kIthbuMJmA8+fVFJQNG4CjR4EaNYAmTYDMTHXf7t3AH38A165dP87LS1VEy8lRmyP5+gK1aqmfvr6qdkBurvX969dXVdoaNgQaNFBb/frqZ+3a6jkHBpas4mZWVgEUo7FyxxeWkABERRU9l06nxmyqwXTejAzgxhurrhbN+fPnMWnSJCxcuBAiAn9/f0ybNg0vvvhimWszWaojo9frER0dzanXVCnVvo5MeQvivfCCKoA3ahQwf74TG072kZCgyqympV2/zfxBOGgQcOyYCk5SU1XBoH/+Ufv+848KXjIyHJvXodGoIEurvZ44ce7c9cSIqubrC9SpA5w9W/a+330HDB5c8nZ7FVBJSFBTeYq/1ZivU/HSutZUpJqdi3BWLZryrqxtxhW2yREcFch45NCSySTSooXqtv3+eyc1lMovM/N6f/t3313vt4+Lszy8Yc+tVq3Kn8NguP5cDAbbjlm/XuTCBTUGcfSoGp7auVPk7bdtO75rV5GQEJF27UQaNxbRaivWdj8/kTZtRG67TSWXjR0r8thjth0bG2v9/2leXsmkkOLjKnp92WM0bp6nExtr26W87z71UrCnvLw8iYmJEX9//4JhohdeeEEyMjLs+0BEZWCOTD5bApn0dPWm4OUlcvmykxpK5RMfXzL/xM9PJCjItiCmdm2RmjWL3taggcikSSJ//KHyUhwZCBX/QLf1k8taEGBrIFQ4eBJRUXxWlkhKisjChY5/zoAKurZuFTl1qmSuT0WfR/HXhr3ydJzE1ssAiPj6ioweLXLsmH3bcPLkSRk+fLgAKu+lcePGsnjxYjEajfZ9ICIrmCOTz5YcmR9+UKMNQUHAgQNObjBdV3hooHlzlX+xY4eahxofb//HMw9dvP02MGWK/c9fXOFxgcqOJZjzU9LS1OdbcbYkVdh6jt271RBUerrazENzqalqWCg7u+znYVa7tnrMG28EWrdWuUfLlpV9XGwsMGyY9edgjzwd8/mcMDxly/+KRo3UZdu6Vd1Wsybw5JPAG28Abdvary2bNm3CSy+9hL179wIAbr31VsyfPx/dunWz34MQWcChpXy29Mi8/776ZvN//+ekRlZHZU3n+PJLNfRRFT0Exb+1N2zo+McoPjxiHlKx1ptky5CKuSei+DkqMmupouewdrx5u/VWkVtuqXyP17p1lh/fHj06hZ+LE4enbPlfYTKJ/PSTSP/+1+/38hJ5/HGRvXvVeewx+SsnJ0c++OAD8fPzEwCi1WrllVdesfs3ZaLCOLSUz5ZA5okn1BvAv//tpEZ6otLePS19QLRooQr4vPaaSPv2VR/A2DtQadSo/AGBvQIRW4uXOOocth5/7ZrIwYMq72fBAvX/PiKi5JCftU2nE7nzTlWAZfZsleA2Z45tx5aWp2N+Di4wPFWe/xW//ipyzz1F973lFpGmTe0Xi6WmpkpERIQAargpICBAvvrqKzFZKgdAVEkMZPLZEsh066b+wNeudVIjPU1p32StfUC42tawYem9I7YEKhUJCOwRiNjjK3hlz1HZon6l9er4+VX+/29pPTL2Sjiu7HWo4Cm2bxcJDy+9+ZWNxdavXy9t27YtCGgGDBggBw4cqPgJiSxgIJOvrEAmO/v6F8CjR53YUE9R2jdZWz6EfH3tF4xYCjJsPXbqVPsEKlVd2ddTlHVtz54V2bxZZPFikddfF3ngAZGOHW2bhVW7tsi0aSLffKMyZB2RcGztOVTR0FRenkizZvaJxay5evWqTJ06VXx8fASA1KxZUyZNmiSXnTRjIi8vTwwGg8TGxorBYGBFYQ/AQCZfWYHMrl3qD9vf33KxVLLC0odtWd9kq2Izv0OvXFlyVpNeL/L117bnojgqUCHbVOTa5uSoYabyvGbq1RPp21fkxRdF5s+3fSp7acNTTh6asjUW+/nnyj/W33//LYMHDxZz70xgYKCsreLu7fj4eNHpdAVtACA6nY6Vhd0cA5l8ZQUyX3yh/qBDQ53YSHdj7Zvm1KlVH7SU9iFhrc5MeXJRGKi4J0uv0YAAVf9/2jSRYcNEunQRqVGj4q8/az0y9hyaMp+vnK9BW2fzt22r/jQq+yXOZDJJQkKC6PX6gkBi6NChJdazcwTzMgmFgxiAyyR4AgYy+coKZCZMUH/QkZFObKQ7KWvoyB5bkyZl95hY63Ep/KZV2lpL9shFIddmSwCQna3WuPrySxXkDB58vTpmaVvduuqYgwdFitdVcYGZU+WpQwOI9OwpsmZN5QOaS5cuyeuvvy41atQQAFKrVi2ZNm2aXLt2rXIntsK8cGXxIKZwMMOFK90XA5l8ZQUyAweqP+T//teJjXQXjh46Khyk2NJjYq3HpTzPh70tZMnSpba/buvWFbn9dpFx41Rw85//2HacA2dO2TKbv0UL1eTata/fHhysTlvZmnd79+6VsLCwgoCiffv28uOPP1bupBYYDAarQUzhzWBL0Egup9oHMvPmzZOOHTtK+/btSw1kzAsN//abExvriix9yJf3a56ld09bpyWzx4SczdJrsFkzkVGj1JTvW26pXHK6g2dO2TqCevq06oyqU+f6Pl26iKxYUbnY3mQyybJly6RZs2YFAcUjjzwiaWlpFT9pMbGxsTYFMrFlBY3kkqp9IGNWWo/Mtm3HCv6wL150YiOrWlk9EdbqvISFlS9osfbuaWuQwh4TcrayXoO5uWpoaulStU5Anz4iPj5l/300aaJmXlljp+Gp8nwfOHtWZPJklfts3rdjR9XJVJk/vfPnz8uYMWPEy8tLAEjdunXlo48+ktzc3IqfNB97ZDwbA5l8pQUyX3yRLoBKeKs2yhpzt0edl6lTOduHqq/cXNuHl4KCREaMUOtc7dp1/e+gsmtvFZKXnSc7Zxlk8+hY2TnLIHnZpf+tnTunJm7Vr1+0mV98oZ5aRe3YsUN69+5dEFx069ZNfvnll4qfUK7nyFhK9mWOjPtjIJOvtEBm0qQMAUQeesiJDaxKZY25m6cml9XTYsvUZQYqVN1Z+tJQv75aMfzGGy3/DdWpI3LHHdfLjVeyR6YytWwuXFDVzguv2NGunep8qmhAYzQaZdGiRdKwYcOCYGPkyJFy+vTpip1Qrs9aKh7McNaS+2Mgk6+0QOahhy4WdCB4PFvG3MuztlFlyugTVRelBfSnT6uifG+9pZZaqFu3fH9/tq69Ze3LiI1/q1lZIu+9p9LbzKdo00bVIszJKftpWnLmzBl59tlnC4KOBg0ayMKFCyvcc2Kpjoxer2cQ4+YYyOQrLZDp1OmaACKrVjmvfVWmsom6hbexY5mIS2RveXkif/4p8vHHIiNHlr2w5ogRIr//fj2aKH4ue9ayEZVHOHNm0e87N96o8p6LV0KwtYDxr7/+KsHBwQXBR69eveSPP/6w/ZoVecqs7OtpHBXIaERE4EYsLQNeu3ZtXL2aC2/vHOTkaHDkCNC6tZMb6mjLlwOPPWafcxkMQGgokJQEpKcDAQHqd63WPucnImXJEuDVV4GMDOv7+PkBffqov8HQUKB3b2DrVqB//7LPbzAAYWHlatLly8CCBcD77wNnzljeR6NRP+PigPDw0s+Xl5eH+fPnY/LkycjKyoJGo8GoUaPw73//Gw0aNChX28izWPr8tgcPCmRaA9iLunWBCxcALy+nNtO+jMaSQUZSkm1vbI0aAefOqS9WxWk0gE4HpKQwaCGqKoX/nhs1Anx9gd9+A375RW3nzxfdv2ZN9c3s0KGyzx0bCwwbVqFmZWUBrVoBmZmW7y/v20V6ejomTJiA2NhYAEDTpk3xwQcfYPjw4dCYIyOqVhwVyHjQx303AEDXrh4WxCQkAIGBKmh57DH1MzBQfXXS6a5/VSpOowH0emDhwuu/F78fAKKjGcQQVSWtVvWaDBsGDBwI3H478NprwNq1wNmzwO7dwPz5wKOPAi1bArm5tgUxgPqiU0E7dlgPYgD1XSg1VcVgtjUlAF9++SV+/vlndOjQAadPn8aIESMQFhaGPXv2VLidRMV50Ed+sPpvsJObYU8JCUBEBHDiRNHb09KARx4B6te33tMCqCAlIkL1B7dsWXQfnc62fmIiqjpeXkCXLsCoUWr4ODUVOHwY+OQToHbt0o+tWxe4ckWNFVVAenqhZsCIfkjEo1iOfkiEF4wF9y1frmIrW/Xv3x+7du3CjBkzULt2bWzatAndu3fHhAkTcPHixQq1lagIu2bcVAHryb7fCiCyYIETG2dP5Vk+oFatshN1OX2ayL1ZK+1bfKtZU+TWW1U1PINBxMZ1kczzBx5EvBxH0fee49DJg4gvkhRceJaTrY4dOyYPPvhgQTJwy5Yt5euvvxZTZReFIrfAZN981nNkDgHQYfNmoG9f57bRLhITbcuBiYkBXniBibpE1UFCAhAVVbSXNiBA9bxevAj8/DNw/HjRY2rVAm67DbjzTrXddJPF9wejEXixWQI+zogAIEW6601Qvbwjasfhe7/wgqTgG28EJk8GnngCqFHD9qfx3XffYcyYMThy5AgAYODAgZg3bx7atWtn+0nI7TDZN5+lC1Grlg7Xrp3Iv1/1sLoVS8m8X39t26ykSiT3EZEbsvR+YQ5MRFQ2rsEAbNigAptTp4oe36ABcMcdwF13AQMGAG3aFJz3SrNA+GacsJhzYIIG1xqpbN8F/9Vi5szrs5zatlUBzWOP2R7QXL16FTNmzMCMGTOQk5MDb29vvP7665g4cSJq1apVoUtjb0ajEUlJSUhPT0dAQABCQ0Oh5ZfECnNUIOMRQ0ve3oMEEGnVqpz9nK4gPt5y0YapU+1TCZSIqi+TSWTPHpHZs0WGDi268JJ5a91a5Lnn1DoG5XjPuXRJ5P33i9ahaddOLX1QnpHrQ4cOyaBBgwqGm1q3bi3/+9//HHM9ysFSUT6dTseifJXAgnj5LF2ImjXHCyAycOBlJ7asAkqr1Gke67ZjASwiquZyc0W2bBF55x2R228v/T3G2lZsPaiLF0VmzChaKTgoSO1m69uTyWSSuLi4IoHD/fffL0ePHnXARSibeZmEwkEMwGUSKstRgYxHzFoymdTU6w4dcpzcknIwGtVYt6WRPfNt1qYGcOo0EVVEjRrALbeocaCNG1WNqW++Ue9FN9xg2zmKTfGuUwd4/XU1ovXee0DDhsDBg2qYqWtXNUpuMpV+So1Gg4ceegj79+/Hq6++iho1amDNmjXo2LEjpk+fjpycqntvNxqNiIqKglh4bzbfNnbsWBiNxhL3k3O4TSATExODTp06oVevXiXuM5m6AAA6dXKjQCYpqeS0aksiI9VU6cI4dZqI7KFOHeDee9WXosOHy65DU6OGSjj+/nvg6tUid9WtC0ycqAKad99V1SH271eVIoKDgfj4sgOaOnXq4P3330dycjL69euHq1ev4s0330S3bt2wYcOGSj1VWyUlJeFEKe/NIoLU1FQk2VpQhxzObQKZyMhI7Nu3D9u2bStxn4hKVmvbthzFDZytcNGG0tx6K3D0qErei41VP1NSGMQQkX1ptcC8earH11qhzbw8YO5c4O67VdfL4MHAnDnA338X7FKvHvDWW+pt6+23AX9/YM8eNbHq5puBNWssd0QX1rlzZxgMBixbtgzNmjXDwYMHMWDAADz66KP4559/7PaULUm38b3Z1v3I8dwmkCmdSpP38XFyM8rD1gqcAQFFK4GGhXE4iYgcIzzccgFNvR74/HPVrfLcc6pX+No1YP16NSzVrp3aoqLUbdeuwd8fmDJFfe+aPFn12OzaBTzwANCrF/Ddd6UHNBqNBo8//jgOHDiAMWPGwMvLCytWrECHDh0wa9Ys5OXlOeQSBNj43mzrfuR4HjH9WqO5DMAPSUlpuO22lqWfwFUYjWphE2vfLrgOEhE5S2lTvAEVgezbB6xbp7akpKI5fbVqqZo1996rNr0eGRnAf/6jOnDMxYd79wbeeUfNBC/RCVSsDTvr1MGo0aPx+++/AwC6du2K+fPn47bbbrPzUzciMDAQaWlpFvNkNBoNdDodUlJSOBW7nDj9Op+lrGfgkgAiSUknnNiyctq2TaRJE+szkjSaktV5iYhcUVaWyKpVahp38XISgEjXriITJ4ps3iyn0/Pk1VeLFiS/7TaRn38udL74eDEVq2xu0unEuHKlLFq0SBo2bFgwk2jkyJFy6tQpuz4d86yl4jOXOGupcjj9Op9HBDIrVoj4+qo/0GbNRJo2LXuJASIid2AyiSQni0ybJtK3r4iXV9H3t8aNRZ58Us7992t5/cUL4uNz/a7+/UX2TYsXEzRiLBYMGaERE9QXvDNnzsizzz5bEGDUr19f5s+fL3l2LEdhqY6MXq9nEFMJXKIgn1sPLYkA06apAWNAdbkuW6YGj7nEABF5oowMlTfz7bdqGOrChev31aiBa33C8D8Zgkm/DcHhvFY4ikC0ROnVhWufUkPuv/32G1566SXs3LkTANCzZ0/Mnz/f4uzWimBlX/viEgX53DaQyc4Gnn1WBS4A8MorwAcfMGAhouojNxf49VdVu+Z//1MFZwo55Xcjml0+UuZpjD8ZoL0zTP3baMSCBQswadIkZGVlQaPR4IUXXsB7772HBg0aOOJZUAU5KpDxkFlLLu7MGZX4tmyZClwWLgQ++ohBDBFVLzVrAv36qS9xBw6oQObDD9VtWi0yLts29fRg4vWpz1qtFqNHj8bBgwfxxBNPQESwcOFCBAUFYenSpTCVVbyG3B4DGUc7dkwtx715syqosH69Wq2aiKi6a98eGD8eSEwETp9G+v2jbDosHSWnPjdv3hxffPEFDAYDOnXqhDNnzuCpp57C7bffjj///NPODSdXwkDGkf76S+W7/P03EBgIbNmiVpslIqKiGjaEdsxLSIUOJlguyGeCBsehx4mYNciYOs9idfSwsDAkJyfj/fffh5+fHzZv3oybb74Z48aNQ1ZWlqOfBTkBAxlH2bNHBTGpqUBQkErm7djR2a0iInJZoWFavNNoNgCUCGbMv/8bk/DE+Tlo9PYYQK9Hzk0hwPTpaqgqX82aNfHqq69i//79eOihh2A0GjFr1ix06NABX331lcX6MOS+GMg4wh9/qDHfU6fUIiObNpVcL4mIiIrQaoHB/w3H/yEOaSg6ceMEdPg/xKH9+KH4pO1M/IJbYYIG3snbgDffVF8UO3ZU6yPs2AGIQK/XIy4uDuvXr0fbtm2Rnp6OYcOG4a677sKBQoEPuTfOWrK3pCQ1rfriRVW2ct06gJnzREQ2S0gAXnnZiNZpSQhAOtIRgKO6UHw0W1uwzFxiIvDRayfRfNtaPIhVuBMb4I1C1YVvuEEtuRAeDvTti2s5Ofjggw/w3nvv4dq1a6hZsyYmTJiASZMmwc/PzynP0xaeNAWclX3zuXRBvK1bRerUUQWcwsJUtUsiIiq3vDwRg0EkNlb9tFTrzmQSWb9epGdPkXq4II8iVlbXeEiya9YuWoSveXORl14S2bBBDh88KPfee29BkbtWrVrJqlWrxGQyVfVTLJOlonw6nc5ti/KxIF4+l+2R2b9f5cRkZAB33KHqJNSq5Zy2EBFVIyJqVe3Jk1V6Yi1cQUTdH/BGu3h0PPw/aDIzr+/cuDHk/vuxRafD8MWLcSQ1FQBwzz33YO7cubjxxhud9CyKSkhIQERERIl8Hk3+olRxcXEIN3dPuQkWxMvnkoHM8ePArbeqDPpevYANG1S1XiIiqjImE7BihVp1+6+/1G2tmucgJnwDBl+Nh3btavVlM580aIDtej3e2bsX641GePn44M0338Rrr70GX1/fshfPdBDzwpUnLMzKAtx34UoWxHNVZ84AAweqIKZDB7U2PYMYIqIq5+UFDBumFub+9FOgVSvg+ElvDJk/GG0Nn+CzmSdh/P4n4MUXgaZNoTl/Hj3//BNrjUZk1KiBhdnZ+G3KFNzUuTN2Tp6symb07w889pj6GRioEngcLCkpyWoQAwAigtTUVCQlJTm8Le6AgUxlZGUBgwer6pR6PfDDD0Djxs5uFRFRtVajBvD008ChQ8C8eUDz5sDRo8DIZ2ug88t34uv+C2A68Y/KGI6MBJo3R928PIwE8B2AbUeOoPu//w0pHkykpQEREQ4PZtLT08veqRz7eToGMhWVl6ey4bdvV8HLjz+qYIaIiFyCj4+KUw4fBt5/H2jYUH3vfOQRoEeIFt9e6geZm19Yb+NGYPRomJo2RR0AmvytCHMmxtixatjJQQICSlYursx+no6BTEW99prKhfHzU8sOBAU5u0VERGRB7drAq68CKSnA22+r0f/kZOC++1R6o2GTFrj9dmDuXHjFxpZ+MhFV6DQx0WHtDQ0NhU6nK0jsLU6j0UCv1yM0NNRhbXAnDGQqYtkyYNYs9e/PPwd69HBue4iIqEz16qlE4JQU9V20Vi21cswddwB33QVs3Qrg9GnbTvbQQ8DLL6t19Oy8MKVWq8Xs2arCcfFgxvx7dHS0WyX6OpJTApkHH3wQDRo0QEREhDMevnJ27ACee079e9IkwM2mvxERVXeNGgEzZ6ohp8hItSj3Tz+pGqZvxdg2XGPMzALmzgVuu00lAb/6akFFYXsIDw9HXFwcWrYsOhNXp9O55dRrR3LK9OvExERcvHgRn332GeLi4sp1rFOnX585A/TsqaZb33MPsHZtlUzFIyIixzl6FHjnHeCzzwCYjDiKQLREGrxQ8uPRBA1OQIdH0BpvNtiNwdnZqHHlyvUd2rcHHn1UbXZYX4+VfcvmlB6ZsLAw1HW3Kcp5eSpD7PhxoG1b4MsvGcQQEXmAwEBg8WJg716gV28tolD6wpVjEY3dfu0w9Px51LlyBfPuuAPXhg4FfH3VVKl33gE6dQJuukllGR8/XuG2abVahIWFYdiwYQgLC3PbIMaRyh3IbNq0CUOGDEGLFi2g0WiwevXqEvvExMQgMDAQvr6+6N27N7Zu3WqPtjrX668DBgNQpw6wejVQv76zW0RERHbUoQMwYwawCuGIsLJwZQTisArh+PLL/+CZZ55BNoAxP/+MFklJWDxjBkyff67W26tRQ2UUv/66WvcpNBRYsAA4e9Ypz82TlTuQuXz5MoKDgxETE2Px/hUrVmDcuHGYMmUKduzYgeDgYAwaNAinbU2gckUbNgAffaT+/dlnQOfOzm0PERE5RGgooNMBqzXhCMRRhMGAYYhFGAxojRSsQjhatgTuu88fn3zyCX799Vd0794d58+fxzNjx+KWuXPxx9tvAydPAh9/DPTrB2g0wC+/AC+9pCoE33cfsHw5cPmys5+uZ6jMQk0AZNWqVUVuCwkJkcjIyILfjUajtGjRQqZPn15kP4PBIA899FCZj3Ht2jXJzMws2FJTU6t20cjMTJFWrdTCYy++aP/zExGRS4mPF9Fo1FZ47UnzVqeOyH/+I3Llito/NzdXZs+eLfXq1RMAotFoZNSoUXLu3Dm1Q2qqyIcfitx0U9ET+fmJPP64yLp1Irm5znvCVcRRi0baNUcmJycH27dvx4ABAwpu8/LywoABA7Bly5YKnXP69Onw9/cv2PRVXXRu3Dg1vtm6NfDBB1X72EREVOXCw4G4OKDYhCE0agS0aAFcugSMHw+0awcsWgSI1MDLL7+MgwcP4vHHH4eIYMGCBQgKCsLSpUthatFCHbBjh1o/4a23gBtvVD0yX36pKsS3bAlERak54IXn4BiNqmbN8uXqpwML8ZXGaDQiMTERy5cvR2JiIoxOaodFlYmCUKxHJi0tTQDIr7/+WmS/V199VUJCQgp+v/POO6Vx48ZSq1YtadmyZYn9C3Nqj8y336qoWaMR2bjRvucmIiKXlpcnYjCIxMaqn3l5quPk009F9PrrHStt24osXy5iNKrjDAaDdOrUSQAIALn11lslOTm56MlNJpEtW0RGjxZp3LhoT027diJTp4rMny+i0xW9T6dTXUZVKD4+XnQ6XcHzASA6nU7iy9kOR/XIOCWQqQxLF8IhgUxGhkhAgHrhvPKK/c5LRERu7+pVkehokSZNrscYwcEi33yjYpScnBx5//33xc/PTwCIVquVqKgoyx/iOTnqi/OwYSK1alkezzJv5jGvKgpm4uPjRaPRFAlikD98ptFoyhXMuMXQUuPGjaHVanHq1Kkit586dQrNmze350M53ssvq6Xbg4KAadOc3RoiInIhvr5qJOjIEeDdd1XV4F27VB5vaCiwZUtNvPrqq9i/fz8iIiJgNBoxe/ZsBAUFITY2FlJ4+KhmTVWbLDYWOHUKWLpULRRlSRWt9wSo4aSoqKiibS1ohuQ3Y6zTh5nsGsh4e3ujR48e2LBhQ8FtJpMJGzZsQJ8+fSp17piYGHTq1Am9evWqbDPLlpCgxi29vNQspVq1HP+YRETkdurUUSkvR46o4r6+vmrVgn79VGySkaHHypUrsX79erRr1w4nT57E448/jjvvvBP79+8vecK6ddV07exs6w9qXu9p4UK7VRK2JCkpCSeKrwBepBmC1NRUJCUlOawNtih3IHPp0iUkJycjOTkZAJCSkoLk5GQczy/4M27cOCxatAifffYZ9u/fj1GjRuHy5ct46qmnKtXQyMhI7Nu3D9u2bavUecp08aKaIgeo+f+9ezv28YiIyO01aqRq3/39N/DCC6qMzLp1qibeo48CN944CLt378a7774LX19fGAwGdOvWDa+//jouXbpU9GTp6bY96OjRqpLwu++q8sR2lm5jO2zdz2HKOxZlMBhKjJUBkBEjRhTsM3fuXGnVqpV4e3tLSEiI/Pbbb/YaCnN8jszkydezt65dq/z5iIio2vnrL5XyYk5t0WpFnn9e5MQJkSNHjsiQIUMKPj/1er3ExcWJyWRSBxsMpefJmDcfn6K/9+unMpHtlINi7fO++GYwGGw6n0sm+zqDQwOZEyeuJ1pVcVY4ERF5nuRkkXvvvR5r+PqKvPqqmk+ydu1aCQwMLAgIBg0aJIcOHVLTo3Q664VsNBo1berCBZHPPhO5886i+9aqperT/PCDOlcF5eXliU6ns5jsi/yEX71eL3k2PoZbJPu6vcmTgatXgVtvBR580NmtISIiNxccDHzzDZCUpBbKvnZNlSS78Ubgzz+HYOvWvXjrrbfg7e2N77//Hl26dMG/pk5F9vvvQ8Tyek8iAKKjAX9/4Mkn1dLdx44B772nJqhcvaryPAcOVAtJvfkmcPBguduu1Woxe7Zad0qjKdoO8+/R0dHOX//JrmGRA82bN086duwo7du3d0yPTHLy9YjWjkNhREREImpa9rffinTrdr3zpFkzkXnzRHbvPigDBw4s6O1o2vRFeRBxchxF68gcg17CEW990MBkUp9ho0aJNGhQtCenTx+Rjz9WPTnlYKmOjF6vd5k6MhoRB6Y8O4ClZcA1mssA/JCUlIbbbmtZ+gmsGTgQ+PFHtcL1V1/Zr8FERESFmEzAihVqEODwYXVb69bAO+8IvL3j8cor4/DPP5sBtIQXBKFIQgDSkY4AJCEUotFCpwNSUoBSO0Oys4G1a9Xs2/Xrr0/X9vVVow5PPQXccUcZJ1GMRiOSkpKQnp6OgIAAhIaGlrsnxtLntz0wkAGA778H7r4b8PYGDhxQrygiIiIHyskBPv0UeOcdtcYkAHTtCoSHZ2PqVCt1ZAoxGICwMBsfLD1dDTctWaKWSTDT64ERI9TWtm25n0N5OCqQYY6M0QhMmKD+PWYMgxgiIqoS3t7AqFFqyvZ776mUl927YVMQA9g+SxuAWnV7wgRgzx61ntNLLwH166t6NP/+t1o46vbbVTG+4tPBXRwDmaVL1f/YBg2ASZOc3RoiIqpm/PyAiRNVUb3XXlOFfm0REFCBB9NogF69gJgYFQl99RUwaJC6PSlJDTcFBADPPqsq+5kHbVxk8UpL3CaQcUhl35wcYMoU9e9//UsFM0RERE7QsCEwc6bqofHzK21PEzSaVGzd+h/k5uZW/AF9fVVe6Pr1wPHjajmetm1Vj8ynn6ppVh06AMOHqyGo/v2Bxx5TPwMDVRV8F1C9c2SWL1f/U5o3V1URra1tQUREVIUSEoCICEsrEJhveAjAKnTu3BkxMTHo16+ffR5YBPjlF2DxYmDlSuDyZcv7madjx8UB4eE2nZo5Mo6QPz8eL73EIIaIiFxGeLiKEXS64vdo4OsLDB06GQ0bBmLv3r0ICwvD8OHDcdKcMVwZGo1a9XLJEuDECesjFVW4eGVZqm8g8/vvavP2Bp5/3tmtISIiKiI8XA0WGAxqYewPPwR69gSuXdNg7dqboNEcRt++KwD4YNmyZQgKCsKcOXOQl5dnnwYkJwPnz1u/37x45aRJ1ntuqkD1DWTmzlU/H30UaNbMuW0hIiKyQKtVU6yHDQPGj1cTjuLi1FqRGRle+PXXhxEQkInWrScjK+sSoqKi0KtXL2zZsqXyD27rtKiZM4EWLdQUrO3bK/+45VQ9A5n0dODrr9W/X37ZuW0hIiKykUYDPPQQsHcv8N//qvghPd0HKSnvoGXL06hT52EkJyejb9++eOaZZ3DmzJmKP5it06ICAoCsLGDhQtVldPPNwPz5wIULFX/scnCbQMaus5YWLgRyc9WaSj16VP58REREVahGDeC554C//gJmzFAlYdLSGuHSpRVo1uwggFuwePFiBAUFYeHChTBWJI8lNFQl6RRbZ6mARqNmMx07BmzYoLqNvL2BnTuByEgV4Dz5pJrW7cB5RdVv1lJ2NtCqFXD6tOqV+b//c3yjiYiIHOj8eRXQzJmjFqYEgHr1DMjKigSwHz179sSCBQvQs2fP8p3YPH0KKBqMWJu1lJEBLFsGLFqkuo3MOnRA1hNPwP+ttzhrqdJWrFBBjE4HPPCAs1tDRERUaQ0aqFSVv/5Stey8vICsrP7QaPagZs3P8ccfJxESEoJRo0bh3Llztp84f/qUtCzaSSAtdZanXjdqBERFqRLFW7YATz8N1K6tlv956y07PNOSqlcgI1J0yrWt5ROJiIjcgE6nOkP27FHrQop4ITd3OLy8DkNkJhYuXIGgoCAsXrwYJpPJpnMmIByBchRhMGAYYhEGAwIlBQkopX6MRgPccosqrJeerlI6une3z5Ms/lDVamhp82ZVqdDXV00Za9y4ahpNRETkBL/9BrzxBrBxo/rdyysLJtN7AOagb9+bEBMTg+6lBBjWCvNVoB4eC+LZhbk35vHHGcQQEZHHu+UWVYfmu++Abt0Ak6kegBnQaA7j11874eabQ/Dyyy8jMzOzxLFGoxolstTd4UL18KpRIJOVBaxerf49erRTm0JERFRVNBpg8GA1meiLL9QySSIBABZB5E/MnZuK9u2D8MUXX6DwIE1Skirua425Hl5SksOfQqmqTyCzfr2ach0U5LBxOiIiIlfl5QU88YTKu42OVnm5QAcAq3D69Co8+eQi9O/fH3vzZxvZWg/P1v0cxW0CmUrXkVmzRv0cOtR+jSIiInIzPj5qyOjIEWDyZKB2bQHQB8AmbNw4AcHBT2DChAnw979i0/lsrZvnKNUj2Tc3F2jaVFUZ/OUXVQiPiIiIcPIkMHUqsGiRwGjUADAB+BxNmy5EXt7POH++FkRKFsXTaNQsqZQUtZRCWZjsWxlJSSqIadJEZT4RERERAKB5c2DBAmDfPk1+jVgvACNx+nQizp37ESKARlO0z8M8ayk62rYgxpGqRyCzdq36ed99zr/iRERELqh9e1Xw/vffgX79jAB8AdwP4ApEis5q0lmph+cMnh/IiDA/hoiIyEYhIYDBoMW6dUDHjtkA/ADUB3AGfn6xePvtX3DkiLhEEAN4WCBjcV2rPXuAo0dVEby77qrqJhEREbkdjQa4+25gzx4ffP65oEmTKwCa4PLlx/D22w3Ru/e7+Pvvw85uJgAPC2QsMvfG3HUX4Ofn3LYQERG5ES8vYPhwDY4fr40ZM7Lh63sFQCfs2PEvtG9/Es888ymuXr3q3DY69dGrgjk/hsNKREREFeLrC7z+ug9OnqyN55/PgJfXNYjcisWLn0HjxklYuDDRaW3z7EDmn3+AbdtUH9mQIc5uDRERkVvz9wc+/rgRjh3zwR13HAFgxJUrAzFq1G0IDPwev/+eWuVtcptApkIF8cy9Mb17A82aOaZhRERE1YxOp8GGDTdi69ZraNNmL4AaOHZsEG65pRHCwn7BmTPZVdYWtwlkIiMjsW/fPmzbts32g8yBzP33O6ZRRERE1VivXn74++/O+PzzFNSrtwdAbWzceBsCAi5j1KgDyMlxfBvcJpApt4sXgQ0b1L8ZyBARETnM8OGtcf58Z4wduxFa7V8wGhti4cIOaNjwFGJiMmAyOe6xPTeQ+eEHICcHaNsW6NDB2a0hIiLyaF5eGsya1Q+nTjXFHXesAJCOy5ebYfToRmjV6iR++inXMY/rkLO6AvO06/vvt1JghoiIiOytUSN/bNjwCLZsOQu9/r8AspCW1hwPPVTTIY/nmYGMCLBunfo3p10TERFVuVtu6YqjR5/F7NnfolatRQAckzDjmYFMWhpw9qxaV6l3b2e3hoiIqFry8vLCyy8PQ1paBB5++EPHPIZDzupsf/6pfnboAPj4OLctRERE1VyDBg2waNGbDjm3Zwcy3bo5tx1ERETkUAxkiIiIyG0xkCEiIiK35TaBjM1LFGRnAwcOqH8zkCEiIvJobhPI2LxEwf79gNEINGgAtGxZNY0jIiIip3CbQMZm5mGl4GAWwiMiIvJwnhvIcFiJiIjI4zGQISIiIrfFQIaIiIjclkcFMjXPnQFOnVK5MZ07O7s5RERE5GAeFcjUPrxf/aNdO6B2bec2hoiIiBzOowIZv7/zAxkOKxEREVULnhXIHGYgQ0REVJ14VCBTm4EMERFRteIxgUwN5KL20b/ULwxkiIiIqgWPCWTa4xC8cnOAunWBG25wdnOIiIioCnhMINMN+fVjunYFvDzmaREREVEpPOYTvyCQ4bASERFRtcFAhoiIiNyW2wQyMTEx6NSpE3r16mXxfgYyRERE1Y9GRMTZjSiPrKws+Pv7IzMzE/Xq1QMANNSk4hxaqR0uXAD8/Z3XQCIiIirB0ue3PbhNj0xpumIvAOBagJ5BDBERUTXiEYFMN+wBAFxu09HJLSEiIqKq5CGBjOqRucJAhoiIqFrxkECGPTJERETVkUcEMkFQSxNcCWzn5JYQERFRVfKIQKYG8gAAJt9aTm4JERERVSWPCGSIiIioemIgQ0RERG6LgQwRERG5LQYyRERE5LYYyBAREZHbYiBDREREbouBDBEREbktBjJERETkthjIEBERkdtiIENERERui4EMERERuS0GMkREROS2GMgQERGR22IgQ0RERG7LKYHMN998g6CgILRr1w6ffPKJM5pAREREHqBGVT9gXl4exo0bB4PBAH9/f/To0QMPPvggGjVqVNVNISIiIjdX5T0yW7duRefOndGyZUvUqVMHgwcPxg8//FDVzSAiIiIPUO5AZtOmTRgyZAhatGgBjUaD1atXl9gnJiYGgYGB8PX1Re/evbF169aC+/755x+0bNmy4PeWLVsiLS2tYq0nIiKiaq3cgczly5cRHByMmJgYi/evWLEC48aNw5QpU7Bjxw4EBwdj0KBBOH36dKUbS0RERFRYuXNkBg8ejMGDB1u9/6OPPsJzzz2Hp556CgCwcOFCfPvtt1i8eDHeeOMNtGjRokgPTFpaGkJCQqyeLzs7G9nZ2QW/Z2ZmAgCysrIKbrsEgQnApSuXi9xORERErsH8+Swi9j2xVAIAWbVqVcHv2dnZotVqi9wmIvLkk0/K0KFDRUQkNzdX2rZtKydOnJCLFy9K+/bt5ezZs1YfY8qUKQKAGzdu3Lhx4+YB2+HDhysTepRg11lLZ8+ehdFoRLNmzYrc3qxZMxw4cAAAUKNGDfznP/9B//79YTKZ8Nprr5U6Y2nixIkYN25cwe8XLlzADTfcgOPHj8Pf39+eza92srKyoNfrkZqainr16jm7OW6L19F+eC3th9fSPngd7SczMxOtWrVCw4YN7XreKp9+DQBDhw7F0KFDbdrXx8cHPj4+JW739/fni8pO6tWrx2tpB7yO9sNraT+8lvbB62g/Xl72nTBt17M1btwYWq0Wp06dKnL7qVOn0Lx5c3s+FBEREZF9Axlvb2/06NEDGzZsKLjNZDJhw4YN6NOnjz0fioiIiKj8Q0uXLl3C33//XfB7SkoKkpOT0bBhQ7Rq1Qrjxo3DiBEj0LNnT4SEhCA6OhqXL18umMVUWT4+PpgyZYrF4SYqH15L++B1tB9eS/vhtbQPXkf7cdS11OTPPrJZYmIi+vfvX+L2ESNGYOnSpQCAefPm4YMPPsDJkyfRvXt3zJkzB71797ZLg4mIiIjMyh3IEBEREbkKp6x+TURERGQPDGSIiIjIbTGQISIiIrflkoFMaatnW7Jy5Up06NABvr6+6Nq1K7777rsqaqnrK8+1XLp0KTQaTZHN19e3ClvrmmxZ8b24xMRE3HzzzfDx8UHbtm0LEuGru/Jey8TExBKvSY1Gg5MnT1ZNg13U9OnT0atXL9StWxdNmzbFAw88gIMHD5Z5HN8ri6rIdeT7pGULFixAt27dCgoH9unTB+vWrSv1GHu9Hl0ukCnv6tm//vorhg0bhmeeeQY7d+7EAw88gAceeAB79uyp4pa7noqsRF6vXj2kp6cXbMeOHavCFrumslZ8Ly4lJQX33nsv+vfvj+TkZIwdOxbPPvssvv/+ewe31PWV91qaHTx4sMjrsmnTpg5qoXvYuHEjIiMj8dtvv+HHH39Ebm4uBg4ciMuXL1s9hu+VJVXkOgJ8n7REp9NhxowZ2L59O/744w/ccccduP/++7F3716L+9v19WjXlZvsICQkRCIjIwt+NxqN0qJFC5k+fbrF/R9++GG59957i9zWu3dveeGFFxzaTndQ3mu5ZMkS8ff3r6LWuScAJRZFLe61116Tzp07F7ntkUcekUGDBjmwZe7HlmtpMBgEgJw/f75K2uSuTp8+LQBk48aNVvfhe2XZbLmOfJ+0XYMGDeSTTz6xeJ89X48u1SOTk5OD7du3Y8CAAQW3eXl5YcCAAdiyZYvFY7Zs2VJkfwAYNGiQ1f2ri4pcS0AVPLzhhhug1+tLjabJOr4m7a979+4ICAjAXXfdhc2bNzu7OS4nMzMTAEpdjI+vy7LZch0Bvk+WxWg04quvvsLly5etVvW35+vRpQKZ0lbPtjYmfvLkyXLtX11U5FoGBQVh8eLFWLNmDZYtWwaTyYS+ffvixIkTVdFkj2HtNZmVlYWrV686qVXuKSAgAAsXLkR8fDzi4+Oh1+sRFhaGHTt2OLtpLsNkMmHs2LG49dZb0aVLF6v78b2ydLZeR75PWrd7927UqVMHPj4+ePHFF7Fq1Sp06tTJ4r72fD06ZfVrck19+vQpEj337dsXHTt2xMcff4x3333XiS2j6iooKAhBQUEFv/ft2xeHDx/GrFmz8MUXXzixZa4jMjISe/bswS+//OLsprg1W68j3yetCwoKQnJyMjIzMxEXF4cRI0Zg48aNVoMZe3GpHpmKrJ7dvHlzrrZtgT1WIq9ZsyZuuummImtrUdmsvSbr1auHWrVqOalVniMkJISvyXyjR4/GN998A4PBAJ1OV+q+fK+0rjzXsTi+T17n7e2Ntm3bokePHpg+fTqCg4Mxe/Zsi/va8/XoUoFMRVbP7tOnT5H9AeDHH3+s9qtt22MlcqPRiN27dyMgIMBRzfRIfE06VnJycrV/TYoIRo8ejVWrVuHnn39G69atyzyGr8uSKnIdi+P7pHUmkwnZ2dkW77Pr67ECicgO9dVXX4mPj48sXbpU9u3bJ88//7zUr19fTp48KSIiw4cPlzfeeKNg/82bN0uNGjXkww8/lP3798uUKVOkZs2asnv3bmc9BZdR3ms5depU+f777+Xw4cOyfft2efTRR8XX11f27t3rrKfgEi5evCg7d+6UnTt3CgD56KOPZOfOnXLs2DEREXnjjTdk+PDhBfsfOXJEateuLa+++qrs379fYmJiRKvVyvr16531FFxGea/lrFmzZPXq1fLXX3/J7t27JSoqSry8vOSnn35y1lNwCaNGjRJ/f39JTEyU9PT0gu3KlSsF+/C9smwVuY58n7TsjTfekI0bN0pKSor8+eef8sYbb4hGo5EffvhBRBz7enS5QEZEZO7cudKqVSvx9vaWkJAQ+e233wru69evn4wYMaLI/l9//bW0b99evL29pXPnzvLtt99WcYtdV3mu5dixYwv2bdasmdxzzz2yY8cOJ7TatZinABffzNduxIgR0q9fvxLHdO/eXby9veXGG2+UJUuWVHm7XVF5r+XMmTOlTZs24uvrKw0bNpSwsDD5+eefndN4F2LpGgIo8jrje2XZKnId+T5p2dNPPy033HCDeHt7S5MmTeTOO+8sCGJEHPt65OrXRERE5LZcKkeGiIiIqDwYyBAREZHbYiBDREREbouBDBEREbktBjJERETkthjIEBERkdtiIENERERui4EMERERuS0GMkREROS2GMgQERGR22IgQ0RERG6LgQwRERG5rf8HkET8dRJa3ZMAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# measured pion, kaon, proton pt spectra\n", "plt.errorbar(pt[\"pion\"], dndpt[\"pion\"], yerr=dndpt_err[\"pion\"], fmt='o', color='black')\n", "plt.errorbar(pt[\"kaon\"], dndpt[\"kaon\"], yerr=dndpt_err[\"kaon\"], fmt='o', color='blue')\n", "plt.errorbar(pt[\"proton\"], dndpt[\"proton\"], yerr=dndpt_err[\"proton\"], fmt='o', color='red')\n", "plt.xlim(0.,3.)\n", "plt.ylim(1,1e4)\n", "plt.yscale(\"log\")\n", "\n", "# fit parameters\n", "Tkin = m.values['Tkin']\n", "beta_s = m.values['beta_s']\n", "n = m.values['n']\n", "\n", "# fit results\n", "ptv = np.linspace(0., 3., 100)\n", "dndpt_bw_pion = dndpt_blastwave(ptv, mass[\"pion\"], Tkin, beta_s, n);\n", "dndpt_bw_kaon = dndpt_blastwave(ptv, mass[\"kaon\"], Tkin, beta_s, n);\n", "dndpt_bw_proton = dndpt_blastwave(ptv, mass[\"proton\"], Tkin, beta_s, n);\n", "\n", "A_pion = normalization(\"pion\", Tkin, beta_s, n)\n", "A_kaon = normalization(\"kaon\", Tkin, beta_s, n)\n", "A_proton = normalization(\"proton\", Tkin, beta_s, n)\n", "\n", "plt.plot(ptv, A_pion * dndpt_bw_pion, color='black')\n", "plt.plot(ptv, A_kaon * dndpt_bw_kaon, color='blue')\n", "plt.plot(ptv, A_proton * dndpt_bw_proton, color='red')" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Finally we draw the $1\\sigma$ error ellipse for the $T$ and $\\beta_s$ with the aid of 'get_cov_ellipse'" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "def get_cov_ellipse(cov, centre, nstd, **kwargs):\n", " \"\"\"\n", " Return a matplotlib Ellipse patch representing the covariance matrix\n", " cov centred at centre and scaled by the factor nstd.\n", " \"\"\"\n", "\n", " # Find and sort eigenvalues and eigenvectors into descending order\n", " eigvals, eigvecs = np.linalg.eigh(cov)\n", " order = eigvals.argsort()[::-1]\n", " eigvals, eigvecs = eigvals[order], eigvecs[:, order]\n", "\n", " # The anti-clockwise angle to rotate our ellipse by \n", " vx, vy = eigvecs[:,0][0], eigvecs[:,0][1]\n", " theta = np.arctan2(vy, vx)\n", "\n", " # Width and height of ellipse to draw\n", " width, height = 2 * nstd * np.sqrt(eigvals)\n", " return Ellipse(xy=centre, width=width, height=height,\n", " angle=np.degrees(theta), **kwargs)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAG+CAYAAADfgAMRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy89olMNAAAACXBIWXMAAA9hAAAPYQGoP6dpAABdJ0lEQVR4nO3de1xUdf4/8NcMMAyCgAiCgxiKJuIFFAUxMy0KU/OaWpq3XG9pbbJp+pXUtJYurquZpdtqmrbrJS8ppYV4y0RMLmkqioiACCgqIHeY+fz+8MdZpxkQEDgwvJ6Px3nEnPM5n3mfz7ry8nNuCiGEABERERHJQil3AURERERNGcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMjKXuwAypNPpcPPmTTRv3hwKhULucoiIiKgKhBC4f/8+NBoNlMqqz3cxjDVAN2/ehJubm9xlEBERUQ2kpqaiTZs2VW7PMNYANW/eHMCD/zFtbW1lroaIiIiqIjc3F25ubtLv8apiGGuAyk9N2traMowRERE1MtW9xIgX8BMRERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGjTaMrVu3Du7u7lCr1fD398eZM2cqbFtaWorly5fDw8MDarUa3t7eOHTokF6bEydO4KWXXoJGo4FCocC+ffsM+pkyZQoUCoXeMmjQIL02d+/exYQJE2Brawt7e3tMmzYNeXl5tXLMREREZHoaZRjbsWMHgoODsXTpUsTExMDb2xtBQUG4deuW0fYhISHYsGED1q5di4sXL2LWrFkYOXIkYmNjpTb5+fnw9vbGunXrKv3uQYMGIT09XVr++9//6m2fMGECLly4gPDwcISFheHEiROYMWPG4x80ERERmSSFEELIXUR1+fv7o3fv3vj8888BADqdDm5ubnjzzTexcOFCg/YajQaLFy/GnDlzpHWjR4+GlZUVtm3bZtBeoVBg7969GDFihN76KVOmIDs72+isGQBcunQJXl5e+O2339CrVy8AwKFDhzB48GDcuHEDGo2mSseXm5sLOzs75OTkwNbWtkr7EBERkbxq+vu70c2MlZSUIDo6GoGBgdI6pVKJwMBAREZGGt2nuLgYarVab52VlRVOnjxZ7e8/duwYWrVqhU6dOmH27Nm4c+eOtC0yMhL29vZSEAOAwMBAKJVKREVFVdhncXExcnNz9RYiIiJqGhpdGMvKyoJWq4Wzs7PeemdnZ2RkZBjdJygoCKtWrUJCQgJ0Oh3Cw8OxZ88epKenV+u7Bw0ahG+++QYRERH4+OOPcfz4cbz44ovQarUAgIyMDLRq1UpvH3Nzczg4OFRYGwCEhobCzs5OWtzc3KpVFxERETVejS6M1cSaNWvQsWNHeHp6QqVSYe7cuZg6dSqUyuod/iuvvIJhw4ahW7duGDFiBMLCwvDbb7/h2LFjj1XfokWLkJOTIy2pqamP1R8RERE1Ho0ujDk6OsLMzAyZmZl66zMzM+Hi4mJ0HycnJ+zbtw/5+flITk5GfHw8bGxs0L59+8eqpX379nB0dMTVq1cBAC4uLgY3EZSVleHu3bsV1gYAlpaWsLW11VuIiIioaWh0YUylUsHX1xcRERHSOp1Oh4iICAQEBFS6r1qthqurK8rKyrB7924MHz78sWq5ceMG7ty5g9atWwMAAgICkJ2djejoaKnNkSNHoNPp4O/v/1jfRURERKbJXO4CaiI4OBiTJ09Gr1694Ofnh9WrVyM/Px9Tp04FAEyaNAmurq4IDQ0FAERFRSEtLQ0+Pj5IS0vDsmXLoNPpsGDBAqnPvLw8aYYLAJKSkhAXFwcHBwe0bdsWeXl5eP/99zF69Gi4uLggMTERCxYsQIcOHRAUFAQA6Ny5MwYNGoTp06dj/fr1KC0txdy5c/HKK69U+U5KIiIialoaZRgbN24cbt++jSVLliAjIwM+Pj44dOiQdFF/SkqK3vVgRUVFCAkJwbVr12BjY4PBgwdj69atsLe3l9qcPXsWAwcOlD4HBwcDACZPnozNmzfDzMwM586dw5YtW5CdnQ2NRoMXXngBK1asgKWlpbTft99+i7lz5+K5556DUqnE6NGj8dlnn9XxiBAREVFj1SifM2bq+JwxIiKixqfJPGeMiIiIyJQwjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRERGRjBjGiIiIiGTEMEZEREQkI4YxIiIiIhkxjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhk1GjD2Lp16+Du7g61Wg1/f3+cOXOmwralpaVYvnw5PDw8oFar4e3tjUOHDum1OXHiBF566SVoNBooFArs27fPoI93330X3bp1g7W1NTQaDSZNmoSbN2/qtXN3d4dCodBbPvroo1o7biIiIjItjTKM7dixA8HBwVi6dCliYmLg7e2NoKAg3Lp1y2j7kJAQbNiwAWvXrsXFixcxa9YsjBw5ErGxsVKb/Px8eHt7Y926dUb7KCgoQExMDN577z3ExMRgz549uHz5MoYNG2bQdvny5UhPT5eWN998s3YOnIiIiEyOQggh5C6iuvz9/dG7d298/vnnAACdTgc3Nze8+eabWLhwoUF7jUaDxYsXY86cOdK60aNHw8rKCtu2bTNor1AosHfvXowYMaLSOn777Tf4+fkhOTkZbdu2BfBgZuztt9/G22+/XePjy83NhZ2dHXJycmBra1vjfoiIiKj+1PT3d6ObGSspKUF0dDQCAwOldUqlEoGBgYiMjDS6T3FxMdRqtd46KysrnDx58rFqycnJgUKhgL29vd76jz76CC1btkSPHj3w6aefoqys7LG+h4iIiEyXudwFVFdWVha0Wi2cnZ311js7OyM+Pt7oPkFBQVi1ahX69+8PDw8PREREYM+ePdBqtTWuo6ioCO+++y5effVVvfT71ltvoWfPnnBwcMCpU6ewaNEipKenY9WqVRX2VVxcjOLiYulzbm5ujesiIiKixqXRhbGaWLNmDaZPnw5PT08oFAp4eHhg6tSp2LRpU436Ky0txdixYyGEwJdffqm3LTg4WPq5e/fuUKlUmDlzJkJDQ2FpaWm0v9DQULz//vs1qoWIiIgat0Z3mtLR0RFmZmbIzMzUW5+ZmQkXFxej+zg5OWHfvn3Iz89HcnIy4uPjYWNjg/bt21f7+8uDWHJyMsLDwx95Ttjf3x9lZWW4fv16hW0WLVqEnJwcaUlNTa12XURERNQ4NbowplKp4Ovri4iICGmdTqdDREQEAgICKt1XrVbD1dUVZWVl2L17N4YPH16t7y4PYgkJCTh8+DBatmz5yH3i4uKgVCrRqlWrCttYWlrC1tZWbyEiIqKmoVGepgwODsbkyZPRq1cv+Pn5YfXq1cjPz8fUqVMBAJMmTYKrqytCQ0MBAFFRUUhLS4OPjw/S0tKwbNky6HQ6LFiwQOozLy8PV69elT4nJSUhLi4ODg4OaNu2LUpLS/Hyyy8jJiYGYWFh0Gq1yMjIAAA4ODhApVIhMjISUVFRGDhwIJo3b47IyEjMmzcPr732Glq0aFGPI0RERESNRaMMY+PGjcPt27exZMkSZGRkwMfHB4cOHZIu6k9JSYFS+b9Jv6KiIoSEhODatWuwsbHB4MGDsXXrVr27IM+ePYuBAwdKn8uv/Zo8eTI2b96MtLQ07N+/HwDg4+OjV8/Ro0cxYMAAWFpaYvv27Vi2bBmKi4vRrl07zJs3T+86MiIiIqKHNcrnjJk6PmeMiIio8WkyzxkjIiIiMiUMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRERGRjBjGiIiIiGTEMEZEREQkI4YxIiIiIhkxjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyajRhrF169bB3d0darUa/v7+OHPmTIVtS0tLsXz5cnh4eECtVsPb2xuHDh3Sa3PixAm89NJL0Gg0UCgU2Ldvn0E/QggsWbIErVu3hpWVFQIDA5GQkKDX5u7du5gwYQJsbW1hb2+PadOmIS8vr1aOmYiIiExPowxjO3bsQHBwMJYuXYqYmBh4e3sjKCgIt27dMto+JCQEGzZswNq1a3Hx4kXMmjULI0eORGxsrNQmPz8f3t7eWLduXYXf+8knn+Czzz7D+vXrERUVBWtrawQFBaGoqEhqM2HCBFy4cAHh4eEICwvDiRMnMGPGjNo7eCIiIjItohHy8/MTc+bMkT5rtVqh0WhEaGio0fatW7cWn3/+ud66UaNGiQkTJhhtD0Ds3btXb51OpxMuLi7i008/ldZlZ2cLS0tL8d///lcIIcTFixcFAPHbb79JbQ4ePCgUCoVIS0ur8vHl5OQIACInJ6fK+xAREZG8avr7u9HNjJWUlCA6OhqBgYHSOqVSicDAQERGRhrdp7i4GGq1Wm+dlZUVTp48WeXvTUpKQkZGht732tnZwd/fX/reyMhI2Nvbo1evXlKbwMBAKJVKREVFVfm7iIiIqOlodGEsKysLWq0Wzs7OeuudnZ2RkZFhdJ+goCCsWrUKCQkJ0Ol0CA8Px549e5Cenl7l7y3vu7LvzcjIQKtWrfS2m5ubw8HBocLagAdhMTc3V28hIiKipqHRhbGaWLNmDTp27AhPT0+oVCrMnTsXU6dOhVLZMA4/NDQUdnZ20uLm5iZ3SURERFRPGkYaqQZHR0eYmZkhMzNTb31mZiZcXFyM7uPk5IR9+/YhPz8fycnJiI+Ph42NDdq3b1/l7y3vu7LvdXFxMbiJoKysDHfv3q2wNgBYtGgRcnJypCU1NbXKdREREVHj1ujCmEqlgq+vLyIiIqR1Op0OERERCAgIqHRftVoNV1dXlJWVYffu3Rg+fHiVv7ddu3ZwcXHR+97c3FxERUVJ3xsQEIDs7GxER0dLbY4cOQKdTgd/f/8K+7a0tIStra3eQkRERE2DudwF1ERwcDAmT56MXr16wc/PD6tXr0Z+fj6mTp0KAJg0aRJcXV0RGhoKAIiKikJaWhp8fHyQlpaGZcuWQafTYcGCBVKfeXl5uHr1qvQ5KSkJcXFxcHBwQNu2baFQKPD222/jgw8+QMeOHdGuXTu899570Gg0GDFiBACgc+fOGDRoEKZPn47169ejtLQUc+fOxSuvvAKNRlN/A0RERESNRqMMY+PGjcPt27exZMkSZGRkwMfHB4cOHZIurk9JSdG7HqyoqAghISG4du0abGxsMHjwYGzduhX29vZSm7Nnz2LgwIHS5+DgYADA5MmTsXnzZgDAggULkJ+fjxkzZiA7Oxv9+vXDoUOH9O7U/PbbbzF37lw899xzUCqVGD16ND777LM6HA0iIiJqzBRCCCF3EaQvNzcXdnZ2yMnJ4SnLxySEwN27d5GdnY38/HwUFBSgoKCgwp8LCwuhUCikRalUVvqzmZkZmjVrBmtra1hbW0s/V7TO3LxR/vuHiIiqoKa/v/mbgRotIQSuX7+Oa9eu4ebNmxUuJSUlcpcqUalUsLa2hoODA1q2bCktjo6ORn8u//zn5+QREZHpYBijRqG0tBTx8fGIjY1FbGws4uLiEBcXh+zsbLlLq5aSkhKUlJTg3r17SExMrPJ+zZo1qzCoPfxz69at4ebmBnt7eygUijo8EiIiqi0MY9QgCSEQFxeH/fv34+DBg4iLi0NxcXG1+mjZsiU0Gg00Gg0cHBwMThka+7l8BkoIAZ1OByGEwc8Pfy4rK9M73Vl+yvPh//55XV5eHu7evYt79+6hqlcJlH9HVR97Ym1tjTZt2qBNmzZwc3Mz+l8GNiKihoFhjBqU+Ph4bNmyBf/5z3+QkpJSaVuNRoMePXrAy8sLrq6uUvDSaDRo3bp1gz+1p9Vqce/ePdy5cwdZWVm4c+eOtFT2uays7JF95+fn4/Lly7h8+XKFbcoD28Mhzc3NDe7u7mjXrh3c3NxgaWlZm4dMRERG8AL+BqgpXsAfHR2N4OBgnDhxwuj2J598Ej169JAWHx8fg1dPNQVCCOTm5hoNardv38bNmzeRmpqKGzduIDU1Ffn5+TX+LoVCAY1GA3d3dymglf/s7u4ONzc3qFSqWjw6IqLGraa/vxnGGqCmFMZu376NxYsX49///rfeKTszMzMEBgZixIgRGDp0KNq0aSNjlY2TEEJ6o0N5OPvzf1NTU1FQUFCj/pVKJVxdXfUCmru7O9q3bw9PT084OzvzNCgRNSkMYyakqYSx6OhoDB48WO8VUk8++SRmz56N8ePHN8mZr/omhEB2drYUzlJSUnD9+nVpSUpKMnjFV1XZ2dnB09MTnTt3hqenp7S0b98eFhYWtXwkRETyYxgzIU0hjJ0+fRrPP/888vLyAADNmzfHsmXLMHfuXJ76amAKCgqQnJysF9AeDmy3b9+uVn8WFhbo0KGDFM7Kw1qnTp1M9s87ETUNDGMmxNTD2P379+Hj44Nr164BAJ566il89913lb5MnRqu/Px8vaCWmJiIS5cuIT4+HsnJydXqS6PRGJ1Nc3V15SlPImrw+NBXajTWrl0rBbGAgACEh4fDyspK5qqopqytrdGlSxd06dLFYFt+fj6uXLmC+Ph4abl06RKuXLli9FEl5Q/qPXLkiN56GxsbvXBWHtY6dOjAmVQiavQ4M9YAmfrMmLe3N86dOwfgwaMsOnXqJHNFVN+0Wi2Sk5OlcPZwWMvKyqpyP2ZmZmjfvr3BTJqnpydatGhRh0dARGSIpylNiCmHMSEEzMzMIIRA165dcf78eblLogYmKyvLYCYtPj4eSUlJVX5ILgA4OzvD09MTXl5e8PHxQY8ePdCtW7cG//w5Imq8GMZMiCmHsbKyMulOut69e+PMmTMyV0SNRWFhIRISEgyC2uXLl1FYWFilPszMzNC5c2eDZ9bZ29vXbfFE1CQwjJkQUw5jQgg88cQTSE1NhVqtxr179zhTQY9Fp9MhNTXV6CnPzMzMKvXRrl07vYDWo0cPtG7dmjcNEFG1MIyZEFMOYwAwZcoUbNmyBQCwbt06vPHGGzJXRKbq7t27OH/+vPSC+djYWFy8eBFarfaR+7Zq1cogoHl4eECpVNZD5UTUGDGMmRBTD2NnzpyBv78/AMDJyQmJiYlo3ry5zFVRU1FUVIQLFy7oBbTff/+9Sm8iaN68Oby9vfUCmpeXF+/oJCIADGMmxdTDGACMHTsWu3btAgBMnToVGzdu5Ckhko1Wq0VCQoJeQIuNjcWdO3ceua+FhQW6du0q3STQo0cPeHt78x8YRE0Qw5gJaQph7OrVq/Dy8kJpaSkA4G9/+xs+/fRTBjJqMIQQuHHjhkFAS0lJeeS+CoUCHTp0MDjNyVd8EZk2hjET0hTCGADs2LEDr776qvS4gg8++ACLFy+WuSqiyt25cwdxcXGIi4uTAlp8fDx0Ot0j99VoNOjRowd69uyJPn36oE+fPnBwcKiHqomoPjCMmZCmEsYA4KuvvsKMGTOkz2vWrMFbb70lY0VE1VdQUGBwo8C5c+eMvmXgzzw9PREQEICAgAD07dsXnTt35k0CRI0Uw5gJaUphDABWrlyJ+fPnS5/feOMNrFq1CpaWljJWRfR4ysrKEB8fb3CaMycnp9L97Ozs4O/vj759+yIgIAD+/v6ws7Orp6qJ6HEwjJmQphbGACAkJAQffvih9Ll3797YtWsXnnjiCRmrIqpdQghcv34dv/32GyIjIxEZGYmYmBjp2kljFAoFunTpojd79uSTT/L6SqIGiGHMhDTFMCaEwKZNmzBnzhzp1E6LFi2wbds2DB48WObqiOpOYWEhYmJicOrUKSmgZWRkVLqPg4MD+vTpI82e+fn5wcbGpp4qJqKKMIyZkKYYxsrFxsbi5ZdfxrVr16R1ISEhWLJkifQaJSJTVj57Vh7MTp06hd9//73SB9UqlUp0795dmj0LCAiAh4cHZ8+I6hnDmAlpymEMALKzszFlyhR8//330rru3btj/fr1CAgIkLEyInnk5+fj7NmzerNnWVlZle7j5OSkd2qzV69eaNasWT1VTNQ0MYyZkKYexoAHswP/+Mc/sHDhQr0ZgenTp+Ojjz7i4wCoSRNC4OrVq3qzZ3/88Uelj9cwNzeHt7e3dGozICAATzzxBGfPiGoRw5gJYRj7nzNnzmDWrFmIjY2V1jk6OmLlypWYNGkSf5EQ/X/379/HmTNnpNmz06dP4969e5Xu07p1a73ZM19fX97FTPQYGMZMCMOYvrKyMqxbtw7vvfce7t+/L61/5pln8MUXX8DLy0vG6ogaJp1Oh8uXL+vNnl28eLHSfaysrNC3b18MHDgQAwYMQO/evfneTaJqYBgzIQxjxqWlpWHevHnSOy2BB+8FnD17NkJCQuDk5CRjdUQNX3Z2NqKioqTZs6ioKOTm5lbYvlmzZujXrx8GDBiAgQMHwtfXlzfSEFWCYcyEMIxV7tChQ5gzZ47eHZc2NjZ45513EBwczBc0E1WRVqvFxYsXERkZiZMnT+LYsWNITU2tsL2NjQ369esnzZz17NkT5ubm9VgxUcPGMGZCGMYerbCwEB999BFWrlyJgoICab2TkxMWL16MmTNnQq1Wy1ghUeMjhEBSUhKOHj2KY8eO4ejRo0hLS6uwva2tLZ5++mlp5szHxwdmZmb1WDFRw8IwZkIYxqouPT0dy5cvx1dffaV316WrqysWL16M119/nRckE9VQ+V2b5cHs6NGjlT6Q1s7ODv3798ezzz6LwMBAdOnShTfZUJPCMGZCGMaqLyEhASEhIdi5c6fe+rZt2+K9997DpEmTeCEy0WMSQuDKlSt6M2e3bt2qsL2zszOee+45BAYG4rnnnkPbtm3rsVqi+scwZkIYxmouLi4OS5cuxf79+/XWazQavP3225gxYwZfukxUS4QQiI+Pl2bNjh07VunDaDt27CgFs4EDB/J5gWRyGMZMCMPY4zt79iyWLFmCgwcP6q1v3rw5Zs6cib/+9a9o06aNTNURmSadToeLFy8iIiIChw8fxrFjx5CXl2e0rUKhgK+vrzRz9tRTT8HKyqqeKyaqXQxjJoRhrPacPn0an3zyCfbt24eH/6hbWFhg/PjxeOedd9C1a1cZKyQyXaWlpThz5owUziIjI1FWVma0raWlJfr16yeFs549e/JmAGp0GMZMCMNY7bt8+TJWrVqFLVu2oLi4WG/b4MGDMX/+fDzzzDO82JioDuXl5eGXX37B4cOHcfjwYZw7d67Ctvb29ggMDMTQoUPx4osvolWrVvVYKVHNMIyZEIaxupOZmYm1a9fiiy++MHhVTK9evRAcHIxRo0bxDkyienDr1i0cOXIEhw8fRnh4OFJSUoy2UygU8Pf3x0svvYShQ4eiW7du/IcTNUgMYyaEYazu5eXlYdOmTVi1ahWSk5P1tjk5OWHatGmYOXMm3N3d5SmQqIkRQiAxMVE6pRkREVHhuzXd3NwwdOhQDB06FAMHDuS1ZtRgMIyZEIax+lNWVoZdu3bh008/1XsZOfDgX+ODBw/GG2+8gaCgIF6/QlSPtFotTp8+jbCwMISFheGPP/4w2s7Kyko6nTlkyBC4urrWc6VE/1PT39/KOqypzq1btw7u7u5Qq9Xw9/fHmTNnKmxbWlqK5cuXw8PDA2q1Gt7e3jh06FC1+rx+/ToUCoXR5eH3JRrbvn379to9eKoV5ubmePXVVxEdHY1jx45h3Lhx0rv3hBD44YcfMGTIEHTo0AEfffRRpc9UIqLaY2ZmhqeeegqhoaE4f/48kpKS8PnnnyMoKEjvmYGFhYU4cOAAZs6ciTZt2sDX1xdLly7Fb7/9Bp1OJ+MREFWDkMmlS5fEqVOnRGlpaY323759u1CpVGLTpk3iwoULYvr06cLe3l5kZmYabb9gwQKh0WjEDz/8IBITE8UXX3wh1Gq1iImJqXKfZWVlIj09XW95//33hY2Njbh//77UDwDx9ddf67UrLCys8rHl5OQIACInJ6dGY0OPJz09XXzwwQfCzc1NANBbLCwsxPjx48Uvv/widDqd3KUSNUn3798Xe/fuFdOmTRMuLi4G/z8tX5ydncXrr78u9uzZo/d3NFFdqenv7zoPY5999pkYM2aM3rq//OUvQqlUCqVSKTp37ixu3bpV7X79/PzEnDlzpM9arVZoNBoRGhpqtH3r1q3F559/rrdu1KhRYsKECTXuUwghfHx8xOuvv663DoDYu3dvdQ5HD8NYw1BWVia+//57MWjQIKFQKAz+ou/WrZtYs2aNuH37ttylEjVZWq1W/Pbbb2Lp0qXC19e3wmCmUqlEUFCQWLt2rUhKSpK7bDJRDTaM+fj4iFmzZkmfjx8/LhQKhRg/frwIDQ0V1tbWIjg4uFp9FhcXCzMzM4PAM2nSJDFs2DCj+zg4OIh///vfeusmTJggnnjiiRr3efbsWQFA/Prrr3rrAQiNRiNatmwpevfuLTZu3FitWRSGsYbn6tWrYv78+aJly5ZGZ8tGjhwpvv/+e1FSUiJ3qURN2o0bN8S//vUvMWzYMGFlZVVhOOvSpYt49913xS+//FLjMzREf9Zgw1iLFi3EZ599Jn1+6623hKurqxRO5s+fLzp27FitPtPS0gQAcerUKb318+fPF35+fkb3efXVV4WXl5e4cuWK0Gq14ueffxZWVlZCpVLVuM/Zs2eLzp07G6xfvny5OHnypIiJiREfffSRsLS0FGvWrKnweIqKikROTo60pKamMow1UIWFhWLr1q0iICDA6F/wTk5OYt68eeL333+Xu1SiJq+goED8+OOP4o033jB62UH54uDgIF577TWxfft2ce/ePbnLpkaswYYxKysrsXHjRulz165dxdSpU6XPGzduFFZWVtXqsybB6datW2L48OFCqVQKMzMz8eSTT4o33nhDqNXqGvVZUFAg7OzsxMqVKx9Z73vvvSfatGlT4falS5ca/QuCYaxh++OPP8T8+fMrvGalR48ePI1J1EDodDpx7tw58fe//10EBAQYvfQAgDAzMxMDBgwQK1euFPHx8XKXTY1MTcNYnd9N2aZNG+kpy8nJybhw4QIGDhwobc/KykKzZs2q1aejoyPMzMyQmZmptz4zMxMuLi5G93FycsK+ffuQn5+P5ORkxMfHw8bGBu3bt69Rn9999x0KCgowadKkR9br7++PGzduGDz5vdyiRYuQk5MjLampqY/sk+TXpUsXfPLJJ0hNTcUPP/yAMWPG6N3lFRsbi7/+9a/QaDQYNWoUDhw4gNLSUhkrJmq6FAoFunXrhkWLFuHUqVPIzMzEli1bMGbMGL1HEGi1Whw7dgzvvPMOPD094enpiaVLl+LixYsyVk+mrs7D2JgxY/Dll19i7ty5ePnll9G8eXMMHTpU2h4XFycFoqpSqVTw9fVFRESEtE6n0yEiIgIBAQGV7qtWq+Hq6oqysjLs3r0bw4cPr1GfGzduxLBhw+Dk5PTIeuPi4tCiRYsKn+puaWkJW1tbvYUaD3NzcwwePBg7d+5Eeno61q1bh969e0vbS0tLsXfvXgwbNgxt2rTBX//6V5w6dYq33RPJyMnJCZMmTcLOnTtx+/ZtREREYN68eejQoYNeu8uXL2P58uXo0qULunfvjg8//BBXr16VqWoyWXU0UycpKCgQkydPFi1atBDt2rUT+/btk7ZlZ2eLZs2aiUWLFlW73+3btwtLS0uxefNmcfHiRTFjxgxhb28vMjIyhBBCTJw4USxcuFBqf/r0abF7926RmJgoTpw4IZ599lnRrl07vesDHtVnuYSEBKFQKMTBgwcN6tq/f7/46quvxPnz50VCQoL44osvRLNmzcSSJUuqfGy8gN80POo0ppubm/jb3/4mfvvtNz4mg6gBuXz5svjHP/4h+vfvX+HpzJ49e4qPP/6Yd2aSngZ7zVhltFqtyM7OrvEdaGvXrhVt27YVKpVK+Pn5idOnT0vbnnnmGTF58mTp87Fjx0Tnzp2FpaWlaNmypZg4caJIS0urVp/lFi1aJNzc3IRWqzXYdvDgQeHj4yNsbGyEtbW18Pb2FuvXrzfatiIMY6altLRUhIWFiZdfflmoVCqjf7G3b99eLFq0SMTFxTGYETUgaWlpYvXq1RXetANA9OnTR/zzn/8UN27ckLtckllNf39X+3VIxcXF+Pnnn3Hr1i20bdsW/fr143vBahlfh2S6srOz8f3332PHjh0IDw9HWVmZQRtPT0+MGzcO48aNQ+fOnWWokoiMSU5Oxq5du7B9+3ZER0cbbFcoFOjXrx/GjRuHl19+Gc7OzjJUSXKql3dTZmVloV+/fkhISJDWqVQqvPbaa1ixYoV0ofvq1auxfv165OTkwM3NDc8//zymTZtW7WvDmiqGsabhzp072LNnD3bs2IGjR48avYasW7duUjD787UsRCSfq1evYufOndixY4d0k9rDlEolBg4ciHHjxmHUqFFo2bKlDFVSfauXMPb222/js88+AwB06tQJd+/exe3bt6FQKODm5oYjR45g9+7dWLhwIR7uVqFQwMzMDAsWLMAHH3xQjcNqmhjGmp7MzEzs3r0bO3bswC+//AJj/7f09fXFuHHjMHbsWDzxxBMyVElExly6dAk7duzAjh07EB8fb7Dd3NwcgYGBGDduHEaMGAF7e/v6L5LqRb2EMU9PTyQkJGDz5s2YOHEihBD46aefMG/ePFy5cgV+fn64cuUK7t27h+effx5jxoxBZmYm9u7di5iYGCgUCsydOxdr1qyp0UE2FQxjTVtaWhp27dqFHTt24PTp00bbBAQEYNy4cRgzZgw0Gk09V0hExgghcP78eezYsQPbt2/HtWvXDNqoVCoMGjQI48ePx/Dhw6FWq2WolOpKvYQxa2trlJWVoaioCAqFQlqflZWFgIAA6Q9eQEAAfvnlF702//3vfzF9+nQUFhbi119/RZ8+fapcZFPDMEblkpOTsXPnTmzfvh0xMTEG2xUKBZ5++mnpX9wMZkQNgxAC0dHR2LFjB3bu3ImUlBSDNvb29njllVcwZcoU+Pn56f3OpMapXsKYlZUVLCwskJuba7Dtq6++wsyZM6FQKPDtt9/ilVdeMWjz0Ucf4f/+7/8wdepUbNy4scpFNjUMY2TM1atXpVMh58+fN9qmZ8+eeOmllzB06FD07NkTSmWdP0qQiB5Bp9MhKioK27dvx65du5Cenm7QxtPTE1OmTMHEiRP5j6pGrF7CWIcOHZCUlISEhASDi/Fv374NZ2dnKBQKo9uBBxcsOzk5oWPHjrh8+XKVi2xqGMboUR51jQoAuLi4YMiQIRg6dCgCAwNhY2NTz1US0Z9ptVocP34cW7Zskd7k8jClUokXXngBU6ZM4WnMRqhewtj06dOxceNGjBkzBjt27DDYbmVlhZKSEuTn51f4B8je3h5lZWXIy8urcpFNDcMYVZUQAufOncOePXsQFhZm9FQm8OA6lYEDB2Lo0KEYOnQo3N3d67dQIjJw//59fPfdd9i8eTNOnDhhsJ2nMRufegljly5dQo8ePVBaWornn38ey5cvh5+fn7T9+vXrOH36tNFTlABQVlYGKysrmJmZoaioqMpFNjUMY1RTaWlp+PHHHxEWFobw8HAUFhYabde1a1cpmPXp0wdmZmb1XCkRPSwxMRHffPMNtmzZguTkZIPtnp6eeP311zF16lQ4OjrKUCFVRb2EMQDYsmUL/vKXv0jPRHJ2doa/vz969eoFX19f+Pr6Vvi+xu3bt2P8+PFwc3Mz+oeNHmAYo9pQWFiIY8eOISwsDAcOHKjwBfQtW7bEiy++iKFDhyIoKIi33RPJSKfT4fjx49i8ebPR05iWlpYYM2YMZs+ejYCAAM6WNTD1FsYAICoqCu+++67etOrDfyBcXV31wlmvXr0QFxeHsWPHIicnB1OmTOEF/JVgGKPaVn7LfVhYGMLCwnD69GmjzzIzMzPD008/Ld0E8OSTT8pQLREBjz6N6e3tjdmzZ2PChAm8JrSBqNcwVi41NRVHjx7Fb7/9hujoaJw7d04vxf85sQshoFQqsWHDBgwdOpSviqgAwxjVtdu3b+PgwYMICwvDoUOHcP/+faPtOnbsKJ3O7NevH1QqVT1XSkQAkJCQgA0bNuDrr7/G3bt39bY1b94ckyZNwuzZs9GlSxeZKiRApjD2ZzqdDpcuXUJMTAyio6MRExODuLg4g4v1y0Oai4sLevToIS2jRo2qrVIaNYYxqk8lJSU4efIkDhw4gAMHDiAxMdFoO1tbWwQFBWHo0KF48cUXK7wcgYjqTmFhIXbu3Ikvv/wSUVFRBtv79++P2bNnY9SoUfzHkwwaRBgzRgiBK1euSOEsOjoasbGxBs8qUyqVRl+a3BQxjJFcyv//Wn4685dffoFWqzVop1AoEBAQIM2ade3aldeuENWzmJgYfPnll/j2228NbtZp1aoVZs6ciblz56JVq1YyVdj0NNgwVpGrV69K4aw8oN25c0eOUhochjFqKLKzs/HTTz8hLCwMP/74o8HpkXJt27aVbgB45plnYGdnV8+VEjVd2dnZ+Oabb/DFF18YPMNTrVZjypQp+Nvf/oYOHTrIVGHT0ejCGFWMYYwaorKyMpw+fVqaNbtw4YLRdmZmZujduzcCAwPx3HPPISAgAJaWlvVcLVHTI4TAsWPH8MUXX2Dv3r16s9oKhQKjRo3C/Pnz4e/vL2OVpo1hzIQwjFFjkJSUhB9++AFhYWE4evQoSkpKjLazsrJC//798dxzzyEwMBDe3t58TRNRHUtJScHq1avx1VdfGVy33b9/f8yfPx+DBw/m/xdrGcOYCWEYo8YmLy8PR44cQUREBCIiIiqcNQMePNfs2WefRWBgIAIDA42+Oo2Iase9e/ewYcMGrFmzBhkZGXrbOnfujPnz52P8+PGcva4lDGMmhGGMGrubN2/iyJEjOHz4MA4fPoy0tLQK27q7u0vB7Nlnn+VdmkR1oLi4GNu2bcPKlSsN3mer0WiwaNEiTJ8+naHsMTGMmRCGMTIl5Xdolgezo0ePIicnp8L23t7eGDhwIAYOHIj+/fvzjQBEtUin0yEsLAyffvopTp48qbetTZs2CAkJwdSpU/lYjBpiGDMhDGNkysrKyhATEyOFs19//bXC680UCgV69OghhbN+/frxTk2iWhIZGYmPP/4Y33//vd56d3d3LFmyBJMmTeJ7a6uJYcyEMIxRU1JQUIBff/1VCmexsbFGX9UEPHgeoa+vLwYOHIgBAwagX79+aN68eT1XTGRaYmNjsXTpUhw4cEBvvZeXF/7+979j2LBhfI5gFTGMmRCGMWrK7t69ixMnTuDYsWM4evQozp07V2FbMzMz9OrVS5o5e+qpp2BtbV2P1RKZjjNnzmDp0qU4dOiQ3vq+ffvio48+wtNPPy1TZY0Hw5gJYRgj+p+srCycOHECR48exdGjRyu9U9Pc3Bx+fn7SzFlAQADDGVE1/fLLL3j33XcRGRmpt/61117DypUr+V7pSjCMmRCGMaKK3bp1C8ePH8fRo0dx7NgxXLp0qcK2ZmZm8Pb2RkBAAPr27YuAgAC4u7vzlAvRIwghsH//fvzf//0fLl68KK23t7dHaGgoZsyYwWeUGcEwZkIYxoiqLiMjQwpnR48exZUrVypt7+zsLAWzgIAA+Pr6wsrKqp6qJWpctFot/v3vf2PRokW4d++etN7Pzw/r169Hjx49ZKyu4WEYMyEMY0Q1d/PmTRw7dgzHjh1DZGQkLly4UOENAQBgYWGBHj166M2eubm51WPFRA3frVu3sGDBAmzZskVap1Qq8dZbb2H58uW8keb/YxgzIQxjRLUnJycHUVFRiIyMRGRkJE6fPl3pc86AB89bKp8569u3L3r06MHnLhEBOH78OGbPnq13eYBGo8HWrVvx7LPPylhZw8AwZkIYxojqjk6nw6VLlxAZGYlTp04hMjLS4Inkf2ZpaQlfX1+905utW7eup4qJGpaSkhL84x//wPLly1FUVATgwSzZihUrsHDhwiZ9LRnDmAlhGCOqX3fv3sXp06el2bOoqCiDlyv/mbu7u97sWffu3WFhYVFPFRPJLykpCTNnzkR4eLi0bvDgwdi6dSscHBxkrEw+DGMmhGGMSF5arRZ//PGH3uzZ1atXK93HysoKvXv31ps943s2ydTpdDp8+OGHWLp0qXRt5hNPPIFdu3ahd+/eMldX/xjGTAjDGFHDc/v2bWnmLDIyEmfOnEFhYWGl+3To0EFv9qxr1658vQyZpPDwcIwfPx5ZWVkAAJVKhQ0bNmDKlCnyFlbPGMZMCMMYUcNXWlqKc+fO6c2eXb9+vdJ9bGxs0KtXL/To0UNaPD09YW5uXj9FE9WhGzduYOzYsdLDYhUKBb799lu8+uqrMldWfxjGTAjDGFHjlJ6erjd7dvbsWRQXF1e6j1qtRrdu3fQCWrdu3dCsWbN6qpqo9pSUlODtt9/Gl19+CeDBWzH27t2LoUOHylxZ/WAYMyEMY0SmoaSkBLGxsVI4O3XqFG7cuPHI/ZRKJTw9PfUCmo+PT5O9KJoaFyEEZs2ahX/9618AHvyD4+DBgxgwYIC8hdUDhjETwjBGZLpu3bqF2NhYvSUhIaFK+z7xxBMGAa1NmzZ8vRM1OFqtFq+99hq2b98O4MEp+sjISHTt2lXmyuoWw5gJYRgjalru37+P33//XS+gXbhwAaWlpY/c19HRET4+PnohrWPHjrxRgGRXWlqKkSNH4ocffgAAvPDCC/jpp59krqpuMYyZEIYxIiopKcGFCxf0Atrvv//+yOefAYC1tTW6d++uF9C6du0KS0vLeqic6H8KCgrQpUsX6eaWn3/+Gc8//7y8RdUhhjETwjBGRMbodDokJiYanOa8devWI/c1NzeHl5eXwWlO/h1Dde0///kPJkyYAADo06ePdLelKWIYMyEMY0RUVUIIpKenGwS0pKSkKu3v4eEBb29vdO7cGZ07d4anpyc6deoEGxubOq6cmgqdToeuXbtK77PMyspCy5YtZa6qbtT09zcfbkNE1IgpFApoNBpoNBoMGTJEWp+dnY24uDi9gHbp0iVotVq9/RMTE5GYmGjQr5ubGzw9PaWlPKi5uLjwhgGqFqVSiRdffFEKY6dPn9b7s0pAo36b57p16+Du7g61Wg1/f3+cOXOmwralpaVYvnw5PDw8oFar4e3tjUOHDlW7zwEDBkChUOgts2bN0muTkpKCIUOGoFmzZmjVqhXmz5+PsrKy2jloIqIqsLe3x4ABAzBv3jx88803OH/+PO7fv48zZ85gw4YNmDVrFvz9/WFlZWV0/9TUVISHh2Pt2rWYM2cOnn32WWg0Gtjb28Pf3x9TpkzBRx99hH379iE+Pr5KNxtQ0/XwTNif/0FAjXhmbMeOHQgODsb69evh7++P1atXIygoCJcvX0arVq0M2oeEhGDbtm346quv4OnpiZ9++gkjR47EqVOn0KNHj2r1OX36dCxfvlz6/PDDGbVaLYYMGQIXFxecOnUK6enpmDRpEiwsLPD3v/+9DkeEiKhy5e/PfPidgWVlZUhKSkJ8fLy0XLp0CZcuXUJ2drZBH7m5uThz5ozBP1TNzc3RoUMHvVm08lOednZ2dX1o1MDFxcVJP7do0UK+QhqoRnvNmL+/P3r37o3PP/8cwINz0m5ubnjzzTexcOFCg/YajQaLFy/GnDlzpHWjR4+GlZUVtm3bVuU+BwwYAB8fH6xevdpoXQcPHsTQoUNx8+ZNODs7AwDWr1+Pd999F7dv34ZKpXrksfGaMSKSmxACt2/flsLZw2HtUa99+jONRmP0lKerqytPeTYBsbGx8Pf3R2lpKRwdHZGenm6yrwBrUteMlZSUIDo6GosWLZLWKZVKBAYGVniXRnFxMdRqtd46KysrnDx5stp9fvvtt9i2bRtcXFzw0ksv4b333pNmxyIjI9GtWzcpiAFAUFAQZs+ejQsXLkizcEREDZlCoUCrVq3QqlUr9O/fX29bQUEBrly5ojeTFh8fj8uXLxt9/dPNmzdx8+ZNHDlyRG+9jY0NOnXqpDeT5unpiQ4dOvAxHCYiPj4eI0aMkE5jz5o1y2SD2ONolCOSlZUFrVarF3gAwNnZGfHx8Ub3CQoKwqpVq9C/f394eHggIiICe/bskc5dV7XP8ePH44knnoBGo8G5c+fw7rvv4vLly9izZw8AICMjw2gf5duMKS4u1vsLLDc3tyrDQEQki2bNmsHHxwc+Pj5667VaLVJSUgxm0i5duoSsrCyDfvLy8hAdHY3o6Gi99eU3JbRr1w7u7u4Gi5ubW5XOMpB8ysrK8K9//QuLFy+WTnf36NED7733nryFNVCNMozVxJo1azB9+nR4enpCoVDAw8MDU6dOxaZNm6rVz4wZM6Sfu3XrhtatW+O5555DYmIiPDw8alRbaGgo3n///RrtS0TUUJiZmaFdu3Zo164dBg8erLctKysLly9fNjjtmZSUBJ1Op9dWCIG0tDSkpaVJZy8eplQq4erqajSolYc1CwuLOj1WMk6n0yE8PBwLFizAuXPnpPU+Pj748ccfGaIr0CjDmKOjI8zMzJCZmam3PjMzEy4uLkb3cXJywr59+1BUVIQ7d+5Ao9Fg4cKFaN++fY37BB5cZwYAV69ehYeHB1xcXAwubC3vs6J+Fi1ahODgYOlzbm4u3NzcKvxOIqLGxtHREY6Ojnjqqaf01hcVFSEhIUFvFi0xMRHXr1+v8GG2Op0OqampSE1NxS+//GKwvTystWnTBm5ubtJ/H/7Z2dmZr4yqJUIInD9/Hjt37sTWrVuRkpKit/21117DunXreA10JRplGFOpVPD19UVERARGjBgB4MH/OSMiIjB37txK91Wr1XB1dUVpaSl2796NsWPHPlaf5XeItG7dGgAQEBCADz/8ELdu3ZLuwAwPD4etrS28vLyM9mFpacnrI4ioSVKr1ejWrRu6detmsC0/Px8pKSlISkrC9evXpaX8s7FTn4B+WKvoOmJzc3NoNBqDwPbwf52cnDjDZoRWq8WVK1cQFxeHU6dO4cCBA0hOTjZo17NnT6xduxZ9+/aVocrGpVGGMQAIDg7G5MmT0atXL/j5+WH16tXIz8/H1KlTAQCTJk2Cq6srQkNDAQBRUVFIS0uDj48P0tLSsGzZMuh0OixYsKDKfSYmJuI///kPBg8ejJYtW+LcuXOYN28e+vfvj+7duwN48CJULy8vTJw4EZ988gkyMjIQEhKCOXPmMHAREVWDtbW19GYAY/Ly8pCcnGwQ0q5fv47k5OQKwxrw4JqmlJQUg1mcP7Ozs0PLli2lxdHRscLP5T9X9Oy2xkYIgZycHCQkJOg9QPjcuXMoKCgwuo9SqURQUBCmTJmC0aNHc/axihptGBs3bhxu376NJUuWICMjAz4+Pjh06JB0sXxKSgqUyv8907aoqAghISG4du0abGxsMHjwYGzduhX29vZV7lOlUuHw4cNSSHNzc8Po0aMREhIi9WFmZoawsDDMnj0bAQEBsLa2xuTJk/WeS0ZERI/PxsYGXbp0QZcuXYxuLyoqwo0bN3Djxg2kpqYa/W9lgQ0AcnJykJOTg2vXrlW5Lisrq0qDW/PmzWFtbY1mzZrB2tpa7+eH/1ubQUYIgeLiYhQUFKCgoAD5+fkoKChATk4O0tPTcfPmTaSlpUl3vpYvhYWFj+zbwsICAwYMwLBhw/Dyyy9XemkPGddonzNmyvicMSKi+lFYWIi0tLQKg9qdO3dw584d3L17F/X969LS0lIKZyqVCgqFAkqlUnr7S2U/l5SU6IWugoICgxslaqpdu3Z6L5zv168fH+z7/zWp54wRERHVBisrK3To0AEdOnSotJ1Wq0V2djbu3LmjF9L+/PnP2x7nNVHljz26d+9ejfuoKXt7e+mdp23atIG3tzd69OgBb29vvTNKVDsYxoiIiB7BzMxMOt345JNPVmkfIQTu379vENTy8vL0Zq2q8t+SkhIIISCEgE6nq/RnnU4HlUolne58+NTnn0+H2tjYoHXr1lLw0mg0aN26td5r/qjuMYwRERHVAYVCAVtbW9ja2qJdu3Zyl0MNmPLRTYiIiIiorjCMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRERGRjBjGiIiIiGTEMEZEREQkI4YxIiIiIhkxjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpJRow5j69atg7u7O9RqNfz9/XHmzJkK25aWlmL58uXw8PCAWq2Gt7c3Dh06VK0+7969izfffBOdOnWClZUV2rZti7feegs5OTl6fSgUCoNl+/bttXfgREREZDIabRjbsWMHgoODsXTpUsTExMDb2xtBQUG4deuW0fYhISHYsGED1q5di4sXL2LWrFkYOXIkYmNjq9znzZs3cfPmTaxcuRJ//PEHNm/ejEOHDmHatGkG3/f1118jPT1dWkaMGFEn40BERESNm0IIIeQuoib8/f3Ru3dvfP755wAAnU4HNzc3vPnmm1i4cKFBe41Gg8WLF2POnDnSutGjR8PKygrbtm2rUZ8AsGvXLrz22mvIz8+Hubk5gAczY3v37q1xAMvNzYWdnR1ycnJga2tboz6IiIioftX093ejnBkrKSlBdHQ0AgMDpXVKpRKBgYGIjIw0uk9xcTHUarXeOisrK5w8ebLGfQKQBrw8iJWbM2cOHB0d4efnh02bNqGRZl4iIiKqY+aPbtLwZGVlQavVwtnZWW+9s7Mz4uPjje4TFBSEVatWoX///vDw8EBERAT27NkDrVZb4z6zsrKwYsUKzJgxQ2/98uXL8eyzz6JZs2b4+eef8cYbbyAvLw9vvfWW0X6Ki4tRXFwsfc7Nza18AIiIiMhkNMowVhNr1qzB9OnT4enpCYVCAQ8PD0ydOhWbNm2qUX+5ubkYMmQIvLy8sGzZMr1t7733nvRzjx49kJ+fj08//bTCMBYaGor333+/RnUQERFR49YoT1M6OjrCzMwMmZmZeuszMzPh4uJidB8nJyfs27cP+fn5SE5ORnx8PGxsbNC+fftq93n//n0MGjQIzZs3x969e2FhYVFpvf7+/rhx44be7NfDFi1ahJycHGlJTU2ttD8iIiIyHY0yjKlUKvj6+iIiIkJap9PpEBERgYCAgEr3VavVcHV1RVlZGXbv3o3hw4dXq8/c3Fy88MILUKlU2L9/v8F1aMbExcWhRYsWsLS0NLrd0tIStra2egsRERE1DY32NGVwcDAmT56MXr16wc/PD6tXr0Z+fj6mTp0KAJg0aRJcXV0RGhoKAIiKikJaWhp8fHyQlpaGZcuWQafTYcGCBVXuszyIFRQUYNu2bcjNzZWu73JycoKZmRkOHDiAzMxM9OnTB2q1GuHh4fj73/+Od955p55HiIiIiBqDRhvGxo0bh9u3b2PJkiXIyMiAj48PDh06JF2An5KSAqXyfxN/RUVFCAkJwbVr12BjY4PBgwdj69atsLe3r3KfMTExiIqKAgB06NBBr56kpCS4u7vDwsIC69atw7x58yCEQIcOHbBq1SpMnz69jkeEiIiIGqNG+5wxU8bnjBERETU+Teo5Y0RERESmgmGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRERGRjBjGiIiIiGTEMEZEREQkI4YxIiIiIhkxjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZNeowtm7dOri7u0OtVsPf3x9nzpypsG1paSmWL18ODw8PqNVqeHt749ChQ9Xus6ioCHPmzEHLli1hY2OD0aNHIzMzU69NSkoKhgwZgmbNmqFVq1aYP38+ysrKauegiYiIyKQ02jC2Y8cOBAcHY+nSpYiJiYG3tzeCgoJw69Yto+1DQkKwYcMGrF27FhcvXsSsWbMwcuRIxMbGVqvPefPm4cCBA9i1axeOHz+OmzdvYtSoUdJ2rVaLIUOGoKSkBKdOncKWLVuwefNmLFmypO4Gg4iIiBov0Uj5+fmJOXPmSJ+1Wq3QaDQiNDTUaPvWrVuLzz//XG/dqFGjxIQJE6rcZ3Z2trCwsBC7du2S2ly6dEkAEJGRkUIIIX788UehVCpFRkaG1ObLL78Utra2ori4uErHlpOTIwCInJycKrUnIiIi+dX093ejnBkrKSlBdHQ0AgMDpXVKpRKBgYGIjIw0uk9xcTHUarXeOisrK5w8ebLKfUZHR6O0tFSvjaenJ9q2bSu1iYyMRLdu3eDs7Cy1CQoKQm5uLi5cuPCYR05ERESmplGGsaysLGi1Wr3AAwDOzs7IyMgwuk9QUBBWrVqFhIQE6HQ6hIeHY8+ePUhPT69ynxkZGVCpVLC3t6+0jbE+yrcZU1xcjNzcXL2FiIiImoZGGcZqYs2aNejYsSM8PT2hUqkwd+5cTJ06FUql/EMQGhoKOzs7aXFzc5O7JCIiIqon8ieRGnB0dISZmZnBXYyZmZlwcXExuo+TkxP27duH/Px8JCcnIz4+HjY2Nmjfvn2V+3RxcUFJSQmys7MrbWOsj/JtxixatAg5OTnSkpqaWoVRICIiIlPQKMOYSqWCr68vIiIipHU6nQ4REREICAiodF+1Wg1XV1eUlZVh9+7dGD58eJX79PX1hYWFhV6by5cvIyUlRWoTEBCA8+fP692BGR4eDltbW3h5eRmtydLSEra2tnoLERERNQ3mchdQU8HBwZg8eTJ69eoFPz8/rF69Gvn5+Zg6dSoAYNKkSXB1dUVoaCgAICoqCmlpafDx8UFaWhqWLVsGnU6HBQsWVLlPOzs7TJs2DcHBwXBwcICtrS3efPNNBAQEoE+fPgCAF154AV5eXpg4cSI++eQTZGRkICQkBHPmzIGlpWU9jxIRERE1dI02jI0bNw63b9/GkiVLkJGRAR8fHxw6dEi6WD4lJUXverCioiKEhITg2rVrsLGxweDBg7F161a9i/Ef1ScA/POf/4RSqcTo0aNRXFyMoKAgfPHFF9J2MzMzhIWFYfbs2QgICIC1tTUmT56M5cuX1/2gEBERUaOjEEIIuYsgfbm5ubCzs0NOTg5PWRIRETUSNf393SivGSMiIiIyFQxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRERGRjBjGiIiIiGTEMEZEREQkI4YxIiIiIhkxjBERERHJiGGMiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERychc7gLIkBACAJCbmytzJURERFRV5b+3y3+PVxXDWAN0//59AICbm5vMlRAREVF13b9/H3Z2dlVurxDVjW9U53Q6HW7evInmzZtDoVDIXY4scnNz4ebmhtTUVNja2spdToPAMTHEMTHEMTGO42KIY2LoccdECIH79+9Do9FAqaz6lWCcGWuAlEol2rRpI3cZDYKtrS3/kvgTjokhjokhjolxHBdDHBNDjzMm1ZkRK8cL+ImIiIhkxDBGREREJCOGMWqQLC0tsXTpUlhaWspdSoPBMTHEMTHEMTGO42KIY2JIrjHhBfxEREREMuLMGBEREZGMGMaIiIiIZMQwRkRERCQjhjGqE+vWrYO7uzvUajX8/f1x5syZStvv2rULnp6eUKvV6NatG3788Ue97Xl5eZg7dy7atGkDKysreHl5Yf369dL269evQ6FQGF127dpVJ8dYXfU9JgCQkZGBiRMnwsXFBdbW1ujZsyd2795d68dWU3KMSWJiIkaOHAknJyfY2tpi7NixyMzMrPVjexy1PS6ZmZmYMmUKNBoNmjVrhkGDBiEhIUGvTVFREebMmYOWLVvCxsYGo0ePblDjIseY/Otf/8KAAQNga2sLhUKB7Ozs2j6sx1LfY3L37l28+eab6NSpE6ysrNC2bVu89dZbyMnJqZPjqwk5/pzMnDkTHh4esLKygpOTE4YPH474+PjqFS6Iatn27duFSqUSmzZtEhcuXBDTp08X9vb2IjMz02j7X3/9VZiZmYlPPvlEXLx4UYSEhAgLCwtx/vx5qc306dOFh4eHOHr0qEhKShIbNmwQZmZm4vvvvxdCCFFWVibS09P1lvfff1/Y2NiI+/fv18txV0aOMRFCiOeff1707t1bREVFicTERLFixQqhVCpFTExMnR/zo8gxJnl5eaJ9+/Zi5MiR4ty5c+LcuXNi+PDhonfv3kKr1dbLcT9KbY+LTqcTffr0EU8//bQ4c+aMiI+PFzNmzBBt27YVeXl5Uj+zZs0Sbm5uIiIiQpw9e1b06dNH9O3bt16O+VHkGpN//vOfIjQ0VISGhgoA4t69e/VxuFUix5icP39ejBo1Suzfv19cvXpVREREiI4dO4rRo0fX23FXRq4/Jxs2bBDHjx8XSUlJIjo6Wrz00kvCzc1NlJWVVbl2hjGqdX5+fmLOnDnSZ61WKzQajQgNDTXafuzYsWLIkCF66/z9/cXMmTOlz126dBHLly/Xa9OzZ0+xePHiCuvw8fERr7/+ek0OodbJNSbW1tbim2++0Wvj4OAgvvrqqxofS22RY0x++uknoVQqRU5OjrQ9OztbKBQKER4e/tjHVBtqe1wuX74sAIg//vhDr08nJyfpz0F2drawsLAQu3btktpcunRJABCRkZG1dmw1JceYPOzo0aMNLozJPSbldu7cKVQqlSgtLX2cw6kVDWVMfv/9dwFAXL16tcq18zQl1aqSkhJER0cjMDBQWqdUKhEYGIjIyEij+0RGRuq1B4CgoCC99n379sX+/fuRlpYGIQSOHj2KK1eu4IUXXjDaZ3R0NOLi4jBt2rRaOKrHI+eY9O3bFzt27MDdu3eh0+mwfft2FBUVYcCAAbV7kNUk15gUFxdDoVDoPUNIrVZDqVTi5MmTtXmINVIX41JcXAzgwXE+3KelpaV0zNHR0SgtLdXrx9PTE23btq3we+uLXGPSkDWkMcnJyYGtrS3MzeV9u2JDGZP8/Hx8/fXXaNeuHdzc3KpcP8MY1aqsrCxotVo4OzvrrXd2dkZGRobRfTIyMh7Zfu3atfDy8kKbNm2gUqkwaNAgrFu3Dv379zfa58aNG9G5c2f07dv3MY/o8ck5Jjt37kRpaSlatmwJS0tLzJw5E3v37kWHDh1q8QirT64x6dOnD6ytrfHuu++ioKAA+fn5eOedd6DVapGenl7LR1l9dTEu5aFq0aJFuHfvHkpKSvDxxx/jxo0b0jFnZGRApVLB3t6+yt9bX+Qak4asoYxJVlYWVqxYgRkzZtTCUT0eucfkiy++gI2NDWxsbHDw4EGEh4dDpVJVuX6GMWoU1q5di9OnT2P//v2Ijo7GP/7xD8yZMweHDx82aFtYWIj//Oc/DWJWrC5VZUzee+89ZGdn4/Dhwzh79iyCg4MxduxYnD9/XsbK686jxsTJyQm7du3CgQMHYGNjAzs7O2RnZ6Nnz55QKk3zr0MLCwvs2bMHV65cgYODA5o1a4ajR4/ixRdfNNljfhSOiaHqjklubi6GDBkCLy8vLFu2rP4LrgfVGZMJEyYgNjYWx48fx5NPPomxY8eiqKioyt8l77wimRxHR0eYmZkZ3IWVmZkJFxcXo/u4uLhU2r6wsBD/93//h71792LIkCEAgO7duyMuLg4rV640mGb+7rvvUFBQgEmTJtXWYT0WucYkMTERn3/+Of744w906dIFAODt7Y1ffvkF69atM7jLsD7J+efkhRdeQGJiIrKysmBubg57e3u4uLigffv2tX2Y1VYX4wIAvr6+iIuLQ05ODkpKSuDk5AR/f3/06tVL6qOkpATZ2dl6s2OVfW99kWtMGjK5x+T+/fsYNGgQmjdvjr1798LCwqKWjqzm5B4TOzs72NnZoWPHjujTpw9atGiBvXv34tVXX61S/U3znwBUZ1QqFXx9fRERESGt0+l0iIiIQEBAgNF9AgIC9NoDQHh4uNS+tLQUpaWlBv8SMTMzg06nM+hv48aNGDZsGJycnB73cGqFXGNSUFAAAFUet/rUEP6cODo6wt7eHkeOHMGtW7cwbNiwxz2sx1YX4/IwOzs7ODk5ISEhAWfPnsXw4cMBPPiFY2FhodfP5cuXkZKSUuH31he5xqQhk3NMcnNz8cILL0ClUmH//v1611PJqSH9OREPbo6Urjmrkipf6k9URdu3bxeWlpZi8+bN4uLFi2LGjBnC3t5eZGRkCCGEmDhxoli4cKHU/tdffxXm5uZi5cqV4tKlS2Lp0qUGjyx45plnRJcuXcTRo0fFtWvXxNdffy3UarX44osv9L47ISFBKBQKcfDgwfo52CqSY0xKSkpEhw4dxNNPPy2ioqLE1atXxcqVK4VCoRA//PBD/Q6AEXL9Odm0aZOIjIwUV69eFVu3bhUODg4iODi4/g78EepiXHbu3CmOHj0qEhMTxb59+8QTTzwhRo0apfe9s2bNEm3bthVHjhwRZ8+eFQEBASIgIKB+DvoR5BqT9PR0ERsbK7766isBQJw4cULExsaKO3fu1M+BV0KOMcnJyRH+/v6iW7du4urVq3qPEqrOYxzqihxjkpiYKP7+97+Ls2fPiuTkZPHrr7+Kl156STg4OFT4SA1jGMaoTqxdu1a0bdtWqFQq4efnJ06fPi1te+aZZ8TkyZP12u/cuVM8+eSTQqVSiS5duhiEhfT0dDFlyhSh0WiEWq0WnTp1Ev/4xz+ETqfTa7do0SLh5ubWYJ4Z9TA5xuTKlSti1KhRolWrVqJZs2aie/fuBo+6kJMcY/Luu+8KZ2dnYWFhITp27Gj0z5Hcantc1qxZI9q0aSMsLCxE27ZtRUhIiCguLtZrU1hYKN544w3RokUL0axZMzFy5EiRnp5eZ8dYXXKMydKlSwUAg+Xrr7+uq8Oslvoek/JHfBhbkpKS6vJQq6y+xyQtLU28+OKLolWrVsLCwkK0adNGjB8/XsTHx1erboUQQlR9Ho2IiIiIahOvGSMiIiKSEcMYERERkYwYxoiIiIhkxDBGREREJCOGMSIiIiIZMYwRERERyYhhjIiIiEhGDGNEREREMmIYIyIiIpIRwxgRUR0KCQmBQqHAxx9/LHcpdWrQoEFQKBQ4cuSI3KUQNToMY0TUIEydOhUKhaJGi7+/v9zlG3Xjxg2sWrUKTk5OmDt3boXtdDod9u7di9dffx1eXl5o2bIlLCws0KJFC3Tt2hUTJ07Et99+i9zc3Fqtb/r06VAoFGjZsiWKi4urvF/Hjh2hUCgwbNgwad2yZcsAAO+88w50Ol2t1klk6hjGiKhBuH79OpydnQ0WGxsbqY2x7c7OzhgwYIB8hVdi8eLFKCwsxIIFC2BtbW20TVRUFLy8vDBq1Ch8/fXXuHTpEnJycmBra4vCwkJcuHAB27Ztw2uvvQY3Nzf885//rLX6pk2bBgC4e/cuvv/++yrtc/z4cVy9elVvfwDo06cPgoKCEBsbi23bttVajURNQnXfiE5EVJ+mT58uAAg3Nze5S6mWGzduCDMzM6FSqcSdO3eMttm7d69QqVQCgGjZsqVYsWKF+OOPP4ROp5PaZGZmiu+++04MHz5cKJVK4e/vX6t1enl5CQBi0KBBVWo/efJkAUA4OzuL0tJSg+MBILp06VKrNRKZOs6MEVGDFhMTAwDo2bOnzJVUz1dffQWtVovBgwfDwcHBYHt8fDxee+01lJSUoHv37jh37hxCQkLQpUsXKBQKqV2rVq0wevRo7Nu3D+fOnUNAQECt1lk+u/Xzzz8jLS2t0rb379/Hd999BwCYNGkSzM3N9baXH+uFCxfw66+/1mqdRKaMYYyIGqyysjL88ccfABpXGBNCYOPGjQCA8ePHG20TEhKC/Px8WFtbY+/evdBoNI/st0uXLo88TXn9+nW8/fbb6NKlC2xsbNCsWTN4enrir3/9K1JSUgzaT5w4ERYWFtDpdNi8eXOlfe/YsQP5+fkAgNdff91gu0qlwujRowEA//rXvx55PET0AMMYETVYFy5ckC4s79Gjh8zVVN0ff/yBGzduAACefvppg+3p6enYs2cPgAdhqH379rXyvd9++y08PT2xZs0aXLx4EWVlZQCAy5cv47PPPkPXrl3x888/6+3j5OQkXYj/qDD29ddfAwD69u0LT09Po2369+8PAPjpp58e51CImhSGMSJqsMpPUQKNa2bsxIkTAAA3Nze4uLgYbD969CiEEACgd0fi4wgPD8ekSZOg1WqxYMECJCUlobCwEPn5+YiPj8eYMWNw//59jBkzxmCGrPxU5dWrV6Xa/+zy5cs4deqUXntjyu9szczMRHx8fG0cGpHJYxgjogYrNjYWwIPrplxdXWWupuqioqIAAN7e3ka3X7x4UfrZx8fnsb9Pp9Nhzpw50Ol0WLduHT7++GO4u7tLj/7o1KkTdu7ciWHDhiE3NxerVq3S2z8oKAht2rQBAGzatMnod5Svt7GxwdixYyuspWPHjtIdsJGRkY99bERNAcMYETVY5TNjjzpFuXnzZigUChw7dqxW2j2umzdvAnhwCtCYO3fuSD8bu7gfeDBL5eLiYnQpn6Eqd+LECSQkJMDR0RF/+ctfKqxr0qRJAAxPISqVSkyZMgUA8N133yEvL09vu1arxdatWwEAY8eO1XvciDEtW7YE8L9xIKLKmT+6CRFR/dPpdPj9998BNK5TlABw+/ZtABUHraooKytDZmam0W0lJSV6n8vvXMzJyan0RoDy/ZKTkw22TZ06FR9++CHy8/OxY8cOvVORBw8eRHp6OoDKT1GWc3BwQHJysjQORFQ5zowRUYN05coVaYamtsLYxIkTUVhYKF1kXleKiooAAJaWlka3l88cAQ8euGqMp6cnhBDSkpSUVOH3lc9AlZaWIjMzs8Ll3r17AIDCwkKDPtq3by89PPfPpyrLP3t6eqJv374V1lHOysoKwP/GgYgqxzBGRA1S+fViQO2FMTMzM6jVaiiVdftXX3nYKg8/f+bl5SX9HBcX99jfp9VqATy4eP7hAFfZYkz5rNepU6dw5coVAA9m+cLCwgAYf5yFMeUB8+HQSUQVYxgjogap/HoxOzs7tGvXrtr763Q6zJo1C2ZmZvjiiy8AGL9mrHxdREQEli9fjrZt20KtVsPX17fCOwsfpfxasYpmvQYOHCg92HX//v01+o6Hld+xaez0Y3WMHj0a9vb2AP43G7Zt2zaUlpbC3NxcuubsUcqPu6Jr5ohIH8MYETVID1+8//AT6auipKQE48aNw6ZNm/Dtt9/ijTfeeOQ+CxcuRFhYGObNm4elS5fi+vXrGDZsGLKzs6tde/nM17Vr14xub926NUaNGgUA2Lp1a6WnIKviqaeeAgBkZGTg7NmzNe5HrVZLD6n95ptvoNVqpWeLDR06FM7Ozo/s4/79+8jKygIAdO7cuca1EDUlDGNE1CCVn76r7inKvLw8DBkyBAcPHsSBAwfwyiuvVHnfU6dOYd68eVi0aBE2btyInJwc/Pe//63W9wP/e/Dp77//Lj209s8++OADWFtbIz8/HyNGjHisOw8HDhyIDh06AADmzZtncIH/n1U0Ywf871Rleno6VqxYgfPnzwOo+inKs2fPQqfTwdzcXAqJRFQ5hjEianCuX78uBYbqPHk/KysLzz77LGJiYnD48GEEBQVVed9Zs2bpvWvxmWeeAfDgERPV9dRTT8Hc3BwlJSUVXhPm6emJbdu2QaVS4dy5c+jevTs++OADXLhwQe+artzcXBw6dAhvvvlmhd9nbm6O9evXw9zcHCdPnkT//v0RERGB0tJSqc21a9ewfv169O7dWzpta0zPnj2lZ5+tWLECwIOZvMGDB1fp2MufsdazZ89HPgKDiB5gGCOiBqemT95//fXXERMTg4iICPTp06da3+nu7q73uUWLFgAqn0WqiK2tLYYMGQKg8mvCRowYgePHj6NTp064c+cO3nvvPXTt2hUWFhZwdHSEnZ0d7Ozs8OKLLyIsLAzNmzfHihUrjB7bc889h127dqF58+aIiopCYGAgrK2t4ejoCLVaDQ8PD8yePRtnz5595Gnf8tkxnU4HAJg8eTLMzMyqdOzlx1vROzmJyBDDGBE1OOV3Upa/5LqqxowZAyEEVqxYIb2XsaoqChsV3Xn4KDNnzgQA/Oc//6m0jz59+uDixYvYvXs3pkyZAk9PT9ja2iInJwdKpRKdO3fGhAkTsGXLFqSnpyMkJARqtdpoXyNGjMDVq1exdOlS+Pn5wcbGBtnZ2bC0tIS3tzf+8pe/YO/evZg/f36ltU+YMEHvO6p6ivLatWuIjIyElZVVlS/2JyI+9JWIGqAVK1ZIp8iqY+LEiejXrx+mTZuGV155Bdu3b9c79VifgoKC4OHhgcTERPzyyy+VPttMqVRi1KhR0kX9j6NVq1ZYtmwZli1bVuM+WrRoYfRZZI+ybds2AMArr7wizSwS0aNxZoyITMrUqVPx73//G3v27MGrr75a7Rmy2qJUKqVA+dFHH8lSQ33Kz8/H2rVrYWlpiaVLl8pdDlGjwjBGRCbn9ddfx1dffYXdu3djwoQJ0kNR69srr7wCPz8/HDx4EGfOnJGlhvry+eefIysrC2+99RaeeOIJucshalR4mpKITNK0adOg0+kwc+ZMKBQKfPvtt/Veg0KhwIYNG7Bv3z6Tf0+jtbU1li1bhrffflvuUogaHYWo6dWpRERERPTYeJqSiIiISEYMY0REREQyYhgjIiIikhHDGBEREZGMGMaIiIiIZMQwRkRERCQjhjEiIiIiGTGMEREREcmIYYyIiIhIRgxjRERERDJiGCMiIiKSEcMYERERkYz+Hx2BH5Wp9yocAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "V = m.covariance\n", "# V_Tkin_betas = V[0:2, 0:2] # sub-matrix giving the covariance matrix for Tkin and beta \n", "V_Tkin_betas = np.array([[V[0,0], V[0,1]],[V[1,0],V[1,1]]])\n", "from matplotlib.patches import Ellipse\n", "fig, ax = plt.subplots()\n", "el = get_cov_ellipse(V_Tkin_betas, (Tkin, beta_s), 1, linewidth=2, fill=False)\n", "ax.add_artist(el)\n", "ax.set_xlim(Tkin - 0.04 * Tkin, Tkin + 0.04 * Tkin)\n", "ax.set_ylim(beta_s - 0.01 * beta_s, beta_s + 0.01 * beta_s)\n", "ax.set_xlabel(r\"$T_\\mathrm{kin}$ (GeV)\", fontsize=18)\n", "ax.set_ylabel(r\"$\\beta_\\mathrm{s}$\", fontsize=18)\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "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.10.12" } }, "nbformat": 4, "nbformat_minor": 4 }