cms_final_lbt_mart_ads_r0.4_0-10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "88c93f3a",
   "metadata": {},
   "source": [
    "## Jet Analysis(Nuclear Modification Factor)( $R_{\\mathrm{AA}}$)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b06e4546",
   "metadata": {},
   "source": [
    "# Preparations\n",
    "Please set path of jet data files for both pp and PbPb.\n",
    "Also if you use jet cone size different from $R=0.4$, please change the value jetR below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1dfb1c9a",
   "metadata": {},
   "source": [
    "# Loading all the files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "dcfa051a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set File Paths\n",
    "# Please set the path for pp jet data file\n",
    "file_pp = '../../data/music_cms/cms_pp_vir0.5_R0.4_trkpT4.0_maxT250_pt200-1001_m1-7_merged_352683.dat'\n",
    "\n",
    "# Please set the value of jet cone size you used in the jet reconstruction\n",
    "jetR = 0.4\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d70a0423",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please set the path for pbpb jet data file\n",
    "file_pbpb = '../../data/music_cms/LBT/R0.4/cms_music_0-10_R0.4_s1-7_35k.dat'\n",
    "file_pbpb1 = '../../data/music_cms/MARTINI/R0.4/cms_martini_0-10_R0.4_s1-6_12790.dat'\n",
    "file_pbpb2 = '../../data/music_cms/ADS/R0.4/cms_ads_0-10_R0.4_s1-4_6461.dat'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "fcee61ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "# preperations\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# pi and 2pi \n",
    "pi = 3.1415926\n",
    "twopi = 2.0*pi\n",
    "\n",
    "# define plot style\n",
    "width = 0.05\n",
    "plotMarkerSize = 8\n",
    "labelfontsize = 15\n",
    "import matplotlib as mpl\n",
    "mpl.rcParams['figure.figsize'] = [8, 6]\n",
    "mpl.rcParams['lines.linewidth'] = 2\n",
    "mpl.rcParams['xtick.top'] = True\n",
    "mpl.rcParams['xtick.labelsize'] = 15\n",
    "mpl.rcParams['xtick.major.width'] = 1.0\n",
    "mpl.rcParams['xtick.minor.width'] = 0.8\n",
    "mpl.rcParams['xtick.minor.visible'] = True\n",
    "mpl.rcParams['xtick.direction'] = \"in\"\n",
    "mpl.rcParams['ytick.right'] = True\n",
    "mpl.rcParams['ytick.labelsize'] = 15\n",
    "mpl.rcParams['ytick.major.width'] = 1.0\n",
    "mpl.rcParams['ytick.minor.width'] = 0.8\n",
    "mpl.rcParams['ytick.minor.visible'] = True\n",
    "mpl.rcParams['ytick.direction'] = \"in\"\n",
    "mpl.rcParams['legend.fontsize'] = 15\n",
    "mpl.rcParams['legend.numpoints'] = 1\n",
    "mpl.rcParams['font.size'] = 15\n",
    "mpl.rcParams['savefig.format'] = \"pdf\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "707d8005",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ratio_error(v1,e1,v2,e2):\n",
    "  #v1, e1: numerator value and error\n",
    "  #v2, e2: denominator value and error  \n",
    "  error1 = e1/v2\n",
    "  error2 = (e2/v2)*(v1/v2)\n",
    "  error = np.sqrt(error1*error1+error2*error2)\n",
    "  return error"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fe7826ff",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load files\n",
    "data_pp = np.loadtxt(file_pp, delimiter=',')\n",
    "data_pbpb = np.loadtxt(file_pbpb, delimiter=',')\n",
    "data_pbpb1 = np.loadtxt(file_pbpb1, delimiter=',')\n",
    "data_pbpb2 = np.loadtxt(file_pbpb2, delimiter=',')\n",
    "\n",
    "# Indices of the data array\n",
    "i_pp = data_pp[:,0] \n",
    "i_pbpb = data_pbpb[:,0] \n",
    "i_pbpb1 = data_pbpb1[:,0]\n",
    "i_pbpb2 = data_pbpb2[:,0]\n",
    "\n",
    "\n",
    "# Get Indices of jets in the data array\n",
    "jet_id_pp = np.where(i_pp < 0.1)\n",
    "jet_id_pbpb = np.where(i_pbpb < 0.1)\n",
    "jet_id_pbpb1 = np.where(i_pbpb1 < 0.1)\n",
    "jet_id_pbpb2 = np.where(i_pbpb2 < 0.1)\n",
    "\n",
    "# Extract jets\n",
    "jets_pp = data_pp[jet_id_pp]\n",
    "jets_pbpb = data_pbpb[jet_id_pbpb]\n",
    "jets_pbpb1 = data_pbpb1[jet_id_pbpb1]\n",
    "jets_pbpb2 = data_pbpb2[jet_id_pbpb2]\n",
    "\n",
    "# Extract associated charged particles for pp\n",
    "assoc_pp= []\n",
    "for i in range(len(jet_id_pp[0])-1):\n",
    "  chunck = data_pp[jet_id_pp[0][i]+1:jet_id_pp[0][i+1]]\n",
    "  assoc_pp.append(chunck)\n",
    "chunck = data_pp[jet_id_pp[0][-1]+1:]\n",
    "assoc_pp.append(chunck)\n",
    "\n",
    "# Extract associated charged particles for pbpb\n",
    "assoc_pbpb= []\n",
    "for i in range(len(jet_id_pbpb[0])-1):\n",
    "  chunck = data_pbpb[jet_id_pbpb[0][i]+1:jet_id_pbpb[0][i+1]]\n",
    "  assoc_pbpb.append(chunck)\n",
    "chunck = data_pbpb[jet_id_pbpb[0][-1]+1:]\n",
    "assoc_pbpb.append(chunck)\n",
    "\n",
    "assoc_pbpb1= []\n",
    "for i in range(len(jet_id_pbpb1[0])-1):\n",
    "  chunck = data_pbpb1[jet_id_pbpb1[0][i]+1:jet_id_pbpb1[0][i+1]]\n",
    "  assoc_pbpb1.append(chunck)\n",
    "chunck = data_pbpb1[jet_id_pbpb1[0][-1]+1:]\n",
    "assoc_pbpb1.append(chunck)\n",
    "\n",
    "assoc_pbpb2= []\n",
    "for i in range(len(jet_id_pbpb2[0])-1):\n",
    "  chunck = data_pbpb2[jet_id_pbpb2[0][i]+1:jet_id_pbpb2[0][i+1]]\n",
    "  assoc_pbpb2.append(chunck)\n",
    "chunck = data_pbpb2[jet_id_pbpb2[0][-1]+1:]\n",
    "assoc_pbpb2.append(chunck)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a38aafd6",
   "metadata": {},
   "source": [
    "## Jet-$R_{\\mathrm{AA}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "760626b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi4AAAGdCAYAAAA1/PiZAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACHLElEQVR4nO3deVxU1fvA8c8RUQRTweWLuO+S5Urm0qJpbqWWoWKrLWqLv0q01Y2vmm2G2mJKX0srAxUtl7TUck0rtVJzN6XMFRO0RJTl+f0xzMQww+oAzvC8X6/7wjn33HPPvQPOM2e7RkRQSimllHIHpYq7AkoppZRSeaWBi1JKKaXchgYuSimllHIbGrgopZRSym1o4KKUUkopt1G6uCvgicqXLy9NmzYt7mqofIqPj6dq1arFXQ2Pp/c5f/R+WZSU++Bp17l9+/YkEfFzZZlGp0O7np+fn1y4cKG4q6HyKSQkhG3bthV3NTye3uf80ftlUVLug6ddpzEGETGuLFO7igpBlSpVirsKqgCGDh1a3FUoEfQ+54/eL4uSch9KynVeCW1xKQQhISHiSRGzUkopVRDa4qKUUkopd3LC1QVq4KKUUkqpwnLc1QVq4KKUUkopt6GBi1JKKaXchq7jUkzOnz/P6dOnSUlJKe6qKKXywNvbm2rVqlGhQoXiropSJZoGLoUgPj6ekJAQhg4d6nRq2/nz5zl16hQ1atSgXLlyGOPSAddKKRcTES5evMixY8cANHhRKhdRUVFERUUBuHx9EJ0OXQhymw596NAhgoKC8PX1LcJaKaWuVFJSEsePH6dhw4bFXRWl3IIxZruIhLiyTB3jUgxSUlIoV65ccVdDKZVP5cqV0+5dpYqZBi7FRLuHlHI/+nerVPHTwMUDGGP0P1SllFIlgscELsaYhsaYWcaYncaYNGPMujwcU9cYI062GCd5+xpjdhljko0xe4wxAwvlQpRSSimVLY8JXIBmQC9gP3Agn8eOAtpn2sZk3mmMuQlYBKwFegJfAtHGmG5XWOcCsbawZG1pyS69MC1atIjbbruNSpUqUbZsWRo3bkx4eDjHj/+7WKK1Pp9++qnD8Z9++qnT+v7+++888MAD1K5dGx8fH2rVqkXfvn3ZsGFDjvXJfP2//PKLLb1u3bqMGjXK9joiIiLbh2EOHjyYkJB/x5KtW7fO4d5at8cee8yWr1OnTrb00qVL06BBA0aMGMH58+dzrLOIMHnyZGrVqkW5cuW45ZZb7OoOMGfOHKfnnzlzpkN5ixcvpmHDhqSlpQGwdetWHn74YRo2bIivry9NmjThv//9L8nJyQ7Hfvfdd9x44434+PhQr1493n77bbv9+/fv56mnniI4OBhfX1/q16/PM888Q2Jiol2+Q4cOMWzYMJo3b46XlxedOnVyONexY8coX748hw8fLtA5AGJiYmjdujXly5enRo0aPPjgg3a/e1mNGDECY4zd7wJYfh+s9zQ0NDTb45VSxc+TpkMvE5ElAMaYWPI3BWu/iHyfw/6xwAYReTrj9VpjTDNgHLCqQLX1ACNHjmTatGk8/PDDjBgxggoVKrBnzx5mzpzJkSNH+Pzzz215y5cvT0xMDPfff79dGdHR0ZQvX55//vnHlpaQkEC7du2oXr06r776KkFBQcTFxbF06VK2bNnCLbfckmu9QkNDady4sUuvd968edSvX98urVq1anavO3fuzOTJk0lNTWXr1q2MHTuWo0ePEhsbm225r732GhMnTuTNN9+kadOmREZG0rVrV3799VcCAwPt8n777bd2A7uz1ic9PZ1x48bx3HPP4eXlBcD8+fP57bffeOGFF2jUqBE7d+5k7Nix7Ny5k0WLFtmOPXToEN27d+fOO+/k1Vdf5ccffyQ8PBxfX19bgLZ69Wq+++47nnjiCZo3b87hw4cZM2YMW7Zs4fvvv6dUKct3od27d7NixQratWuX7WDWGjVqMHDgQCZMmMCcOXNs6Xk9x9KlSxk0aBBPPfUUb775JidOnGDMmDHccccdbN++3ZbPas+ePcyePdvpVObHHnuMHj168OSTT2b7PimlrhIi4nEbEAusy0O+uoAAd+aQpyxwGXg8S/qDQBpQMesxbdq0kZzs2bMnx/25yaizgGSzWfYXpqVLlwogs2fPdtiXmpoqK1assKvvoEGDxNvbW86ePWtL/+uvv8Tb21vuvfdeu/pGRUWJMUZOnTrlUHZ6enqO9QLknXfecUivU6eOjBw50vZ6/PjxUrlyZadlPPTQQ5L5PVy7dq0AsmvXrhzPfeutt8o999xjlzZx4kQB5PTp006PuXjxolSoUEH++9//2tL++ecfqVKliowePdqW9tFHHwkgf//9d451WLVqlfj4+Mj58+dtafHx8Q75Zs2aJYDExcXZ0oYOHSqNGjWSlJQUW9oTTzwhNWvWtN33M2fOOLwHX3/9tQCybt06W1paWprt3/fcc4/ceuutTuu7ceNGKVOmjJw5c8aWltdzDBw4UFq3bm2Xb8mSJQI4/Ru77bbbZMyYMQ6/C5k5ew+zutK/X6VKEmCbuPgz3pO6iq7ERxnjYk4YYyKNMZnnKjcAvIF9WY7Zi6WrzbVf693E1KlTad26NY888ojDPi8vL3r27GmX1r59e4KCguy+4S9atIigoCDat29vlzcxMZEyZcoQEBDgULa7DUJu06YNAHFxcU73b968mfPnzzNgwABbmp+fH71792blypX5Pt/cuXPp1q0b11xzjS3NWZdYq1atAOy6VVauXEm/fv0oXfrfhtiwsDD+/PNPfv31VwAqV67s8B44Kytra0d2OnbsSEBAADEx/w4ry+s5UlJSqFixol2+SpUqAVi/XNjExsayb98+XnzxxTzVSyl19Srpgcsl4D3gUaALMAt4Asg8ONc/42dilmMTsuy3OX78uN04hIiICFfWudilpKSwefNmevTokedjjDEMHDiQ6OhoW1p0dDRhYWEOeVu3bs2lS5d44IEH2L59O+np6S6ptzOpqakOW9YPPau0tLQ857WyBixZu3ys9u3bh5eXF40aNbJLDw4OZt++rLEyNGjQgNKlS9OkSRNmzZrlsP/bb7+lQ4cOOdYJYMuWLZQqVYoGDRoAcOHCBY4ePUrTpk0d6mGtZ05lAQXqmjPG0K5dO9asWZNrfbOe45FHHmHjxo18/PHHnD9/ngMHDjBmzBhuu+02rr32Wlu+ixcvMnLkSF577TX8/PzyXUelVP5kHjMGBLm6/BIduIjICREZLiJLRWSdiEQA4UAfY0yLgpYbFBRk16zlaYHLX3/9xaVLl6hdu3a+jgsLC2P9+vWcOnWKkydPsn79eqeBS5cuXRgxYgTz588nJCSESpUqcc899+T64ZZff/31F97e3g7bxx9/7DR/y5YtHfLOnTvXLo+IkJqayqVLl9i4cSOvvPIKISEh1KxZ02mZCQkJlC9f3jYexcrf35+kpCQuX74MQPXq1Zk4cSKffPIJy5Yto127djz++ONMnTrVdszx48c5ceIE1113XY7XffLkSSZNmsQDDzxgG6NjHfhqbbHIXA9rPZ1JSkrihRde4NZbb7W1LuVXixYt+PHHH7Pdn9057rjjDubMmcPQoUOpWLEiTZo0IS0tza5VD+DVV1+levXqDuOrlFKFIyIiIvOwiuxHyxeQJw3OdZVYYAbQBtjBvy0rFbPks7a0OP8fvQTIb7dNq1ataNiwIQsWLEBEaNy4MS1btmTTpk0OeSMjI3nyySf54osv2LBhA1999RWff/45M2bM4PHHH3dJ/StWrOg0GPrvf//LiRMnHNJjYmJsLRRW9erVs3u9ePFivL29ba87duzIhx9+eMVdXN27d6d79+621z179iQ5OZlJkybxzDPPUKpUKU6ePAk47xqyunz5MgMGDKB8+fJ2QU9BiAiPPvoop0+f5ssvvyxwOVWqVOH06dOIiMN9yukca9eu5fHHH+eZZ56hZ8+enDp1ioiICO6++27WrFmDl5cXR44cYcqUKaxdu9btuhmVUs5p4OJIsvz8DUgBmgLrM+VrCqST/6nXbq9y5cqULVuWP/74I9/HDhw4kJiYGESEgQNzXgqnYcOGjBo1ilGjRnHmzBm6devGyy+/zLBhw1zyIVS6dGm7ac9WlStXdhq4NGvWLNfWjNtuu43XX3+d0qVLU6dOHVuLRXb8/f35559/SEtLs2t1SUhIwNfXlzJlymR7bGhoKAsWLCAuLo769evbpjeXLVvWaX4R4cEHH2T37t189913dnWztrScO3fO7hhrS4uz63jhhRf4/PPPWb16tcPspvwoW7asrestc9CX2zlGjhxJnz59eP31121pLVu2pGnTpixZsoR+/frx4osv0rNnT5o0aWJrVUpPT+fSpUskJiZSsWJFDWiUcjMluqsoG9ZFHLYDiMglLOu39M+SbyCwRUTOUWxMNlvh8vb2pmPHjnz99df5PjYsLMw2rdVZN1F2qlSpwsMPP0xCQgKnT5/O93mLir+/PyEhIbRs2TLXoAWgadOmpKWlcejQIbv0ffv2OYw3ycr6gWv9aR3M7Gy9E4Bnn32WJUuWsGTJEoey/fz8qFWrlsNYFuvrrPmnTp3KlClT+Pjjj7n55ptzrGduEhMTKV++vEPQkts59u3bR8uWLe3SmjRpQrly5fjtt98Ay5owixcvxt/f37YdPXqUd999F39/f9vTnpVS7sNjAhdjjK8xJtQYEwrUAKpaXxtjfDPyHDLGzM50TIQx5i1jTD9jTFdjzARgKrBYRHZmKn4i0MkYM80Y08kY8waWxe4mFN0VXl2effZZtm3b5jDGAyzfaL/66iunxwUHBzN06FCGDh2a7QdzfHy80/SDBw9StmxZh5kk7qxDhw5UqFCBhQsX2tKSkpJYtmyZw8ysrGJjY6lSpQp16tQBLIvslSlThiNHjjjkffXVV3n33Xf59NNPuemmm5yW17NnTz7//HPbwnVgWQOmVq1adi1N8+bNY+TIkURGRtrNhiqouLg4h4G9eTlHnTp1+Omnn+zS9u7dy8WLF6lbty4A//vf/1i7dq3d9p///IcBAwawdu1aqlatesX1V0oVLU/qKqoGLMySZn1dD4jDcr2ZR0Huw7Jq7mNAOeAP4E3glcyFiMimjIBoEpZZR0eAe0WkWBafyzqTxfqNO7cZLq7Uu3dvwsPDefTRR/nuu+/o27cv5cuXZ9++fcycOZO6detmO+vI2Wqvmc2dO5d58+bx4IMP0qJFC1JSUlizZg0zZszgiSeewMfHp0B1vtIugZ07d9otlAeWcTLWmTcF4ePjw4svvsjEiRPx9/e3LUCXnp7O//3f/9ny3XPPPbRt25bmzZuTlpbG/PnzmT9/Pm+//bZt6rGPjw9t2rRh+/btPPzww7ZjP/vsM15++WUGDx5MjRo1+P77f9dabNCgge3D+7nnnmPevHk88MADDBkyhK1btzJr1izef/99271bv349Dz/8MN26daNdu3Z2ZdWsWdM2CDkpKYkVK1YAlhVyz58/b1uEr1evXvj6+tqO27ZtGx07drS9zus5Hn/8cUaMGEFQUJBtjMuECROoW7cuvXr1AnDaFWhdidnZar5KKTfg6oVhdCv8BeiyoggWnMtObGysdOrUSSpUqCDe3t7SqFEjGTlypJw4ccKufs4WhbN655137Oq/e/duefLJJyU4OFjKly8vFStWlNatW8uMGTPsFkdzJrtzVatWTV5++WXb64IsQOds69Kliy1fXhYvcyY9PV0mTZokNWrUEB8fH7npppvkp59+ssvz0ksvSePGjaVcuXLi4+MjrVu3lo8//tihrDfffFMaNGjgcD3Z1f+jjz6yy7tx40a54YYbpGzZslKnTh2ZPn263f7x48dnW9b48eNt+Y4cOZJtviNHjtjynT59Wry8vOwWlsvrOdLT02XGjBly/fXXi6+vrwQFBcmAAQPkt99+y/F+6wJ0ShUdCmEBOiNF+C29pAgJCZFt27Zlu3/v3r1X9C09q+JocblaGWOYPn06Tz75pG0htXPnzlG5cmXee+89hg0bVsw1LFynTp2idu3abNq0iRtuuKG4q5OrWbNmMWXKFA4cOFCsg2TT09NJT0+nS5cuVK1aNcdHNLj671cpT2aM2S4ijk2fV8BjxrhcTeLj4wkJCSEqKqpIzmeNQpXFM888g7e3Nz/88ANr165l2LBheHl52boPPNl//vMfHnvsMaZPn17cVcmViDB9+nRGjx5d7DN7JkyYgLe3d64P8VRK5U1UVJS1qzY/zw3ME08a43LVqFq1Kjm1uKjCs3XrVtu//fz8uOWWWwgODmbZsmXUqlWrGGtWdMaOHcvs2bMdplhfbU6ePMl9993HAw88UNxVYejQodx5550ATh81oZTKH+skDGPMGVeXrV1FhaCou4qUUkVH/36VyjvtKlJKKaVUiaaBi1JKKaXchgYuSimllHIbGrgopZRSym1o4OIJjLFsSimllIfTwEUppZRSbkMDl0JQ1AvQFYeIiAiMMTRq1Mjp/kaNGmGMISIiwmFfvXr1MMbYPRG5bt26GGNy3PKSZ86cOcTFxWW7v2vXrrZzDh482JZeqlQpatasyaBBg4iLi7vi+zN48GCnz8mxWrdunV29vL29qV+/Ps8//zwXLlwAYM6cOXm6J0opdbXRBejcTFEuQJeQkEAscBII/OADQkND8ff3L5Jz+/j4cOTIEbZt22b3Ib1161bi4uKcPgxxy5YttsAgOjqasWPHAvD5559z6dIlW74ePXoQGhrKY489Zkvz8vKye3KxszwNGjSwffBPmTLF7uF9gMOTpZs2bcpHH31Eeno6+/btY/To0fTq1YtffvmFMmXK5PeW5Nu8efOoX78+KSkp/PTTT4wePZrExESioqK444472LJliy1vbGwsb731ll1a2bJlC72OSimVX4W5AJ0GLm5KRJg0bhyRU6bQHWgIrBkxgueffprwUaMYM2FCoS+j7ufnR+vWrYmJibELXGJiYrjtttvYvn27wzHR0dH4+flx3XXX2QUurVq1sstXunRpatasSbt27bI9f3Z5rIFLkyZNcjzeeg3WPB06dMDX15dBgwaxbds2OnTokO1xxhjWrl17xU8Ybt68Oddddx0AN998M8eOHWPOnDlERUVRtWpV25ObAVswnNs1KaWUJ9OuIjc1adw4voiMZHdyMjHAJGD+hQvsTk7mi8hIJo0bVyT1CAsLY8GCBbZnJYkICxYsICwszCFvWloaCxYsoE+fPjzyyCPs3buXHTt2FEk986pFixYAHD16tFjOf80115CSklIs51ZKKXeggYsbSkhIIHLKFJYlJRGUZV8QsCwpicgpU0hMTCz0uvTr149Tp06xadMmADZu3Eh8fDz9+vVzyLt27VpOnTpFWFgYoaGheHt7Ex0dXWh1S09PJzU11W5LT0/P8Zg//vgDsIzDKQppaWmkpqZy8eJFvvvuO2bNmsXdd99dJOdWSil3pIGLu7BOeTaG2IAAuiUnOwQtVkFAt+RkFvr7F/o06UqVKtGjRw9iYmIASzdRjx49HMaSgKWbyJo/ICCAbt26ERMTU2hPtu7bty/e3t5224QJExzypaamkpKSws6dO3nxxRfp0aMHbdu2dciTeYN/gw7rVpDraNmyJd7e3vj6+nLTTTdRr149t3iys1JKFRcNXNzQScD5XJ5/NczIVxTCwsKIjY3l0qVLxMbGOu0munz5MosXL+buu++2DXoNCwvj999/txts6kpTp05l69atdtvQoUPt8mzfvh1vb2/KlClDixYtOH/+vEMrUFxcnEMABNC1a1e7tPXr1+e7jjExMWzdupUffviB6Ohozpw5w1133ZVry5BSSpVUOjjXXWT6Nh/4wQesGTECMgahOnPIz4+uU6fCkCGFXrU+ffrw2GOPMXr0aC5cuEDv3r0d8qxcuZLExER69epl68Lq1KkTZcuWJTo6OseBsAXVsGHDHKckAwQHB/Pxxx+TkpLCpk2bePnllxk2bBjz58+35QkKCmLr1q12x91www3MnDmTNm3a2NKaNGmS7zo2a9bMNji3bdu2NG7cmDZt2rBixQruvPPOfJenlFKeTgMXNxQaGsrzTz/NcXDaXXQMWJWWxqz+/YukPn5+ftx5551MnTqV/v374+fn55DH2orR30mdFi5cyLRp0/Dy8ir0umbl6+trC27at29PcnIy48aNIzw8nBtvvBGAMmXKOA2AmjRpkmtglF/BwcEA7N27VwMXpZRyQruKCkFhL0Dn7+9P+KhR9Pb15XiWfceBPr6+hI8aRaVKlQrl/M488cQT9O7dm8cff9xh34ULF1i2bBmDBg1i7dq1dltkZCSnTp3i22+/LbK65mTkyJFUqVKF119/vVjO/+uvvwJQq1atYjm/Ukq5gi5A52aKYgG6MRmDTJtNmUK35GQaYukeWpWWRnh4uG1/UenUqVO2a5osWbKEpKQknnnmGVsrhlXHjh155ZVXiI6O5vbbb3dpnfbv30+VKvZ/Mz4+PrRs2TLbY3x9fRkxYgRjx47l4MGD2a4MnBcJCQnExsY6pPfq1cv27507d/LPP/+Qnp7O4cOHmThxIrVr16Znz54FPq9SShU3XYBOOTDGMHbiRIaHhxMbEMBJoOvUqczq379IW1ryIjo6mkaNGjkELQDe3t4MGDCAzz77jPfff9+lK8GOGjXKIa1BgwZ2jxpwZvjw4bz55ptMmTKFWbNmFfj8hw8fdto1duTIEdu/77vvPsDyfgYFBXHzzTczadIkp7OylFJKgSmsqaglWUhIiOTU4rJ3717bWAaXsE551vdSqULn8r9fpTyYMWa7iLh0MKC2uHgCDViUUkqVEDo4VymllFJuQwMXpZRSSrkNDVyUUkop5TY0cFFKKaWU29DApRAU9gJ0Siml1NVMF6BzM0WxAJ1SSil1tSrMBei0xUUppZRSbkMDF49gMjallFLKs3lM4GKMaWiMmWWM2WmMSTPGrMvDMTcYYz4yxhwyxiQZY/YbY8YbY3yy5IswxoiTrUehXZBSSimlHHhM4AI0A3oB+4EDeTxmINAAeD3j2PeAcGCek7zngPZZti1XVuUrl5CQwAcfwMSJ8MEHH5CQkFAk542IiMAYk+1DCBs1aoQxhoiICId99erVwxhj98ygunXrYozJcctLnjlz5hAXF5ft/q5du9rOOXjwYFt6qVKlqFmzJoMGDSIuLu6K74+17Mzns7p48SIVKlSw1TerCxcu4Ofnh6+vL3///bfD/jlz5thdU5UqVejcuTMbNmxwuK7stsGDB9vue+ZnOlnf1+7duzucNzQ01O5BmuvWrcMYY3uitVJKFQVPGpy7TESWABhjYsnbSObXRCTzwKF1xphkYJYxpo6I/J5pX6qIfO/C+l4REWHSpHFERk6he3do2BDWrBnB888/TXj4KMaMmYAxhdt95OPjw5EjR9i2bZt19DgAW7duJS4uDh8fH4djtmzZYgsMoqOjGTt2LACff/45ly5dsuXr0aMHoaGhPPbYY7Y0Ly8v0tLScszToEEDLly4AMCUKVPo2LGj3fmzPrywadOmfPTRR6Snp7Nv3z5Gjx5Nr169+OWXXyhTpkx+b4md8uXLs27dOk6dOsV//vMfW/ry5cvJ6RlhS5cuJSkpCYAvvviCBx54wGm+b7/9lnLlynHy5EkmT55M9+7d2blzJ2PHjuXxxx+35Xv55ZdJTExkxowZtrSqVavmWPdVq1axdetWbrjhhjxdq1JKFRWPCVxEJL0Axzgb7fxzxs8g4Hcn+68KkyaN44svItm9O5mgIGvqBY4fh969IwEYO3ZiodbBz8+P1q1bExMTYxe4xMTEcNttt7F9+3aHY6Kjo/Hz8+O6666zC1xatWpll6906dLUrFmTdu3aZXv+7PJYA5cmTZrkeLz1Gqx5OnTogK+vL4MGDWLbtm106NAh2+OMMaxdu9auBSKrJk2a8Pfff7Nw4UKGDx9uS4+JiaFPnz589tlnTo+Ljo6mfv36iAjR0dHZBi433HAD5cuXB6Bt27bUrFmThQsX8vLLL9OgQQNbvoCAANLT03O9F5nz16hRg1deeYUvvvgiT8copVRR8aSuIldpD6QDv2VJr2SMOWOMSTHG/GyM6VcMdQMs3UORkVNYtiwpU9BiERQEy5YlERk5hcTExEKvS1hYGAsWLLC1IIgICxYsICwszCFvWloaCxYsoE+fPjzyyCPs3buXHTt2FHod86NFixYAHD161CXlDRw4kJiYGNvrv//+mxUrVji9P2B5b7/++msGDhxIWFgYq1ev5syZ3GcTBgUFUbVqVZfU2xjD6NGjWbp0Kbt27bri8pRSypU0cMnEGBMIjAE+EZHTmXYdAp4H+gP3AMeBRdkFL8ePH7cbT+BsnEcBamfbYmMD6NYt2SFosQoKgm7dklm40J/Cnm3Ur18/Tp06xaZNmwDYuHEj8fHx9OvneGvWrl3LqVOnCAsLIzQ0FG9vb6Kjowutbunp6aSmptpt6ek5N8z98ccfgGUcjisMGjSIzZs328r9/PPP8ff359Zbb3Waf9GiRVy+fJmwsDAGDRpEamoqsbGxuZ7nwoULnD171mX17t+/P40aNeKVV15xSXlKqZLDOlYuY7hCNp9UBaeBSwZjTBlgAfAPMCLzPhH5VEQiRWStiCwF7gS+B8Y5KysoKAgRsW2uCVz+dfIkZDMm1qZhQ0u+wlapUiV69Ohha1WIiYmhR48eDmNJwNIFYs0fEBBAt27diImJyXG8x5Xo27cv3t7edtuECRMc8qWmppKSksLOnTt58cUX6dGjB23btnXIk3kDSwtS5jRn1xEcHMz111/P/PnzAcv9GTBgAKVKOf/Ti46OJjg4mObNm3P99dfTrFmzbIM76/mPHTvGsGHDqF69Og8//HC+7lF2SpUqxUsvvcTChQs5cCCvY92VUsoSuFg//7B80XcpDVwAYwkLPyZjZpKI5Dg1RyzvxmKguTHGqwiqCIhtCwyM4uBBvxxzHzrkR2BgVMYxhSssLIzY2FguXbpEbGys026Qy5cvs3jxYu6++27boNewsDB+//13tmwpnMlZU6dOZevWrXbb0KFD7fJs374db29vypQpQ4sWLTh//rxDoBAXF+cQAAF07drVLm39+vVO6xEWFkZMTAx//fUXa9asybab6MSJE6xbt85uf1hYGBs3buTPP/90yF+pUiW8vb2pWbMmixYtIjY2NtdBt/lx//33U7t2bV599VWXlamUUldKAxeLaUBfoK+I7MvjMdZIosiFhoayalUax7OJY48dg1Wr0ujfv3+R1KdPnz78888/jB49mgsXLtC7d2+HPCtXriQxMZFevXqRmJhIYmIinTp1omzZsoXWXdSwYUNCQkLstqAs/WvBwcFs3bqVzZs388Ybb/DHH38wbNgwuzxBQUEOARDAzJkz7dLatGnjtB5hYWH89NNPvPrqq9SoUSPbQbILFiwgPT2dHj162O5Rz549ERFbi01mGzZs4Mcff+TTTz8lICCAsLAw28BkVyhdujTPP/88n376Kb//ftWOU1dKlTAeM6uooIwxLwHDgQEisimPxxgsY112iEhabvldzd/fn/DwUfTuHekwQPf4cejTx5fw8HAqVapUJPXx8/PjzjvvZOrUqfTv3x8/P8fWIGtw4iyYWrhwIdOmTcPLq4garzLx9fW1zYhq3749ycnJjBs3jvDwcG688UYAypQpYzdryqpJkyZO07OqV68ebdu2ZerUqTz33HPZ5rPeI+t5s+4bOXKkXVqrVq0oX748N9xwA7Vq1eLWW2/l3Xff5YUXXsi1Tnn1yCOPMGnSJF5//XWXlamUUlfCYwIXY4wvlkXkAGoAFYwxoRmvV4hIkjHmELBeRB7NOOZeYDIwBzhmjMn8Vfg3EYnPyLceWATsA/yAIcCNwF2FelE5GDPGMlajWbMpdOuWTMOGlu6hVavSCA8Pt+0vKk888QSXLl2yWz/E6sKFCyxbtoxBgwY5dNX8/PPPhIeH8+2333L77bcXVXWzNXLkSN5++21ef/11Fi9e7NJyP/vsMx588EGn+w8fPswPP/zAiBEj6NOnj92+lStX8sYbb3Dw4MFsF/y75ZZb6NWrF9OmTePZZ5+lbNmyLql32bJlGTVqFC+99BJt2rSxdZMppVRx8ZjABagGLMySZn1dD4jDcr2Zv9Z3y/g5OGPL7GEsAQ1YZhU9C1THMlX6J+AOEVl5pZUuKGMMY8dOZPjwcGJjAzh5Erp2ncqsWf2LrKUls06dOmW7psmSJUtISkrimWeecWhN6NixI6+88grR0dEuD1z2799PlSr26xD6+PjQsmXLbI/x9fVlxIgRjB07NsdAIb8GDBjAgAEDst0fExNDqVKlGDVqlEN31rXXXktkZCTR0dGMG+d0PDgAo0ePpmPHjnz66ac8+uijLqk3wLBhw5g8eTKbN2/OdjaUUkoVFY8Z4yIicSJistniMvLUFZHBmY4ZnMMxczLle1RE6otIORHxE5GbizNoyczf358hQ2DsWBgyZEixBC25iY6OplGjRk67QLy9vRkwYACLFy+2WznXFUaNGkX79u3tttDQ0FyPGz58OBUqVGDKlCkurU9OoqOj6dKli0PQAlCtWjVuv/32XMcCdejQgc6dOzNlyhSXztSyBnNKKXU1MIU1FbUkCwkJkW3btmW7f+/evQQHB7vwjNa1WvS9VKqwuf7vVynPZYzZLiK5DwbMB0/qKirBNGBRSilVMnhMV9HVJD4+npCQEKKiooq7KkoppVSRi4qKss66zMsDj/NFW1wKQdWqVcmpq0gppZTyZEOHDmXo0KEYY3J/2Fo+aYuLUkoppdyGBi5KKaWUchsauCillFLKbWjgopRSSim3oYGLBzD8u5KLUkop5ck0cFFKKaWU29DARRVIREQExphsn+XTqFEjjDFEREQ47KtXrx7GGA4dOmRLq1u3LsaYHLe85JkzZw5xcXHZ7u/atavtnIMHD7allypVipo1azJo0CDi4uJcdp8uXLiAn58fvr6+/P3333k6Zvny5Rhj7Orx+++/88ADD1C7dm18fHyoVasWffv2ZcOGDQ7HL1iwgMDAQLtl/xctWsRtt91GpUqVKFu2LI0bNyY8PJzjx4/b8mR3zzZt2pTrfc9aX6WUKiy6jkshsC5AZ53HXpgSEhIgNhZOnuSDwEBCQ0Px9/cv1HNa+fj4cOTIEbZt22ZdaAiArVu3EhcXh4+Pj8MxW7ZssX3ARUdHM3bsWAA+//xzu2cV9ejRg9DQUB577DFbmpeXF2lpaTnmadCgARcuXABgypQpdOzY0e78FStWtHvdtGlTPvroI9LT09m3bx+jR4+mV69e/PLLL5QpUya/t8TB0qVLSUpKAuCLL77ggQceyHcZCQkJtGvXjurVq/Pqq68SFBREXFwcS5cuZcuWLdxyyy12+b/88kt69eqFMZYOxJEjRzJt2jQefvhhRowYQYUKFdizZw8zZ87kyJEjfP7557ZjR44c6fA8p+DgYLZs2WJ7ffjwYe677z7ee+89WrdubUuvXr16vq9NKeWZoqKirIuwunwBOkRENxdvbdq0kZzs2bMnx/15kZ6eLmMmTBCfSpWEgQOF0aPFb8AA8alUScZMmCDp6elXfI6cjB8/XipXriy33367jBw50m5feHi4dOvWTSpXrizjx4+32/d///d/4ufnJzfeeKMEBwdnW76zY/Oa58iRIwLIsmXLcjz+oYcekqzvVXR0tADy3Xff5XgsIGvXrs0xj4hI7969pX79+lKvXj3p2bNnrvlFRJYtWyaAHDlyREREoqKixBgjp06dcsib9X1OS0uTqlWrSmxsrIiILF26VACZPXu2w7GpqamyYsUKu2t65513cq3frl278nz9nsgVf79KlRTANnHxZ6x2FbmpcZMmEfnFFyTv3g0xMTBpEhfmzyd5924iv/iCcZMmFUk9wsLCWLBgga1bQkRYsGABYWFhDnnT0tJYsGABffr04ZFHHmHv3r3s2LGjSOqZVy1atADg6NGjV1xWQkICX3/9NQMHDiQsLIzVq1dz5oz9IpIiQkREBNWqVeOaa67hwQcf5Pz583Z5EhMTKVOmDAEBAQ7nsLaqWG3dupXExERuv/12AKZOnUrr1q155JFHHI718vKiZ8+eV3qZSilVpDRwcUMJCQlMiYwkadkyCAqy3xkURNKyZUyJjCQxMbHQ69KvXz9OnTrFpk2bANi4cSPx8fH069fPIe/atWs5deoUYWFhhIaG4u3tTXR0dKHVLT09ndTUVLstPT09x2P++OMPwDIO50otWrSIy5cvExYWxqBBg0hNTSU2NtYuz9tvv82ECRMYOnQosbGxlCtXjueff94uT+vWrbl06RIPPPAA27dvz/EavvzyS26++WYqVKhASkoKmzdvpkePHnmuc9Z7lrlrTimlrgYauLgJk2kLiI0luVs3x6DFKiiI5G7d8F+4sNCnSVeqVIkePXoQExMDQExMDD169HAYSwKWMS3W/AEBAXTr1o2YmBi7QaSu1LdvX7y9ve22CRMmOORLTU0lJSWFnTt38uKLL9KjRw/atm3rkCfzBpYWpMxpWa8jOjqa4OBgmjdvzvXXX0+zZs3sArW0tDRef/11hg0bxqRJk+jevTuzZs0iODjYrpwuXbowYsQI5s+fT0hICJUqVeKee+5hzZo1Dtfy5ZdfcscddwDw119/cenSJWrXrp3ne/bMM8/Y3a9bb701z8cqpVRR0MDFHZ08CdnM5rFp2NCSrwiEhYURGxvLpUuXiI2NddpNdPnyZRYvXszdd99tG/QaFhbG77//bjfw05WmTp3K1q1b7basg6W3b9+Ot7c3ZcqUoUWLFpw/f96hFSguLs4hAALo2rWrXdr69ettx5w4cYJ169bZ3YuwsDA2btzIn3/+CVi6o06cOEHfvn3tzuestSoyMpIDBw7w5ptv0qlTJ7766iu6devGzJkz7c75888/2wIXq6zdSTl57rnn7O7X7Nmz83ysUkoVBQ1c3IRk2qICA/E7eDDH/H6HDhEVGEjhtGXY69OnD//88w+jR4/mwoUL9O7d2yHPypUrSUxMpFevXiQmJpKYmEinTp0oW7ZsoXUXNWzYkJCQELstKEsrVXBwMFu3bmXz5s288cYb/PHHHwwbNswuT1BQkEMABDBz5ky7tDZt2tiOWbBgAenp6fTo0cN2vT179kREmD9/PgAnMwLLatWq2Z0v6+vM1zNq1CiWLl3K77//TsuWLXn55ZdtLT0rVqygfv36NGnSBIDKlStTtmxZW/dXXtSuXdvuflnLUkqpq4UGLm4oNDSUtFWrINMaHHaOHSNt1Sr69+9fJPXx8/PjzjvvZOrUqfTu3Rs/Pz+HPNbgpH///vj7++Pv70+tWrW4dOkSCxcuLLaxFL6+voSEhNC+fXuee+45xo0bx4IFC/jhhx9secqUKeMQAAE0adLELu2aa66xHWO93htvvNF2vdbjrPsCAwMBOH36tF2dsr52pkqVKjz88MMkJCTY8mfuJgLw9vamY8eOfP311/m+L0opdbXSwMUN+fv7Myo8HN/evR2Dl+PH8e3Th1Hh4VSqVKnI6vTEE0/Qu3dvHn/8cYd9Fy5cYNmyZQwaNIi1a9fabZGRkZw6dYpvv/22yOqak5EjR1KlShVef/31Apdx+PBhfvjhB0aMGOFwvc8//zzbt2/n4MGD1KpVi8DAQJYsWWJ3/OLFi+1ex8fHOz3PwYMHKVu2LBUrVuTy5cusXr3aoZvo2WefZdu2bcydO9fh+PT0dL766qsCX6dSShUHXYCuEBTFAnQTxowBYEqzZpaBug0b4nfoEGmrVhEeHm7bX1Q6depEp06dnO5bsmQJSUlJPPPMM9x44412+zp27Mgrr7xCdHS0bQqvq+zfv58qVezXPvLx8aFly5bZHuPr68uIESMYO3YsBw8ezHZl4JzExMRQqlQpRo0a5dA1de211xIZGUl0dDTjxo3j+eefZ9SoUVSpUoWbb76ZRYsWsXfvXrtj5s6dy7x583jwwQdp0aIFKSkprFmzhhkzZvDEE0/g4+PDmjVrEBGHwbS9e/cmPDycRx99lO+++46+fftSvnx59u3bx8yZM6lbt26+Zh0ppVRe6AJ0brYVxQJ0VmfPnhWiooQJEyQqKkoSEhJcVnZOrAvQ5cS6QNydd94pjRo1yjbfE088IRUrVpTk5GSHY/NSflbWBeicbQ0aNLDlc7YAnYjIuXPnpFKlSjJ06NBsz00OC7Bdd911cvvtt2d7bM+ePaVp06YikrGQ4JgxUqVKFSlfvrzce++9Mm/ePLsF6Hbv3i1PPvmkBAcHS/ny5aVixYrSunVrmTFjhqSkpIiIyLPPPit9+/bN9pyxsbHSqVMnqVChgnh7e0ujRo1k5MiRcuLECbtr0gXocqcL0CmVdxTCAnRGCmkqakkWEhIi27Zty3b/3r17Haa8XgnrnBF9J0uuxo0b89xzzzFkyJDirorHc/Xfr1KezBizXURCcs+Zd9pV5AE0YFEHDhwo7ioopVSR0MG5SimllHIbGrgopZRSym1o4KKUUkopt6GBSzHRQdFKuR/9u1Wq+GngUgy8vb25ePFicVdDKZVPFy9etD2rSilVPDRwKQTWBegyFt9xUK1aNY4dO0ZSUpJ+g1PKDYgISUlJHDt2LNtnSSml/hUVFWV9zInLF6DTdVwKQW7ruACcP3+e06dPk5KSUkS1UkpdCW9vb6pVq0aFChWKuypKuQ1dx8WDVKhQQf8DVEoppfJJu4qUUkop5TY8JnAxxjQ0xswyxuw0xqQZY9bl8biKxpiPjDEJxphzxph5xpjKTvL1NcbsMsYkG2P2GGMGuvwilFJKKZUjjwlcgGZAL2A/kJ/1zxcAnYDHgMHADcAXmTMYY24CFgFrgZ7Al0C0MabbFdZZKaWUUvngSWNclonIEgBjTCx5GMlsjGkPdANuFZENGWnHgB+MMV1FZE1G1rHABhF5OuP1WmNMM2AcsMrF16GUUkqpbHhMi4uIpBfgsJ7AKWvQklHOj8CRjH0YY8oCnbG0zGQWA7Q3xlQsWI2VUkoplV+e1OJSEE2BfU7S92bsA2gAeDvJtxdL4NcY2Jp5x/790KmTfeYBA+DJJyEpCXr1cjzh4MGW7cwZCA113P/EEzBwIBw9Cg884Lh/5Ejo3dty7mHDHPePGQNdu8Ivv8CzzzrunzwZOnSAzZvh5Zcd90+bBi1bwpo1MGmS4/5Zs6BJE1i2DN56y3H/J59ArVowfz68/77j/thYqFIF5syxbFmtWAG+vjBjBizIGkIC69ZZfk6ZAsuX2+8rVw5WrrT8e+JE+OYb+/2VK8OiRZZ/v/QSbNliv79mTfj0U8u/n33Wcg8za9wYrEv2DB0KWR/U3LKl5f4B3H8//Pmn/f727eHVVy3/vuce+Osv+/1dusDYsZZ/9+wJWdcuvPNOGDXK8u+sv3egv3v6u2f5t/7uOe7X3z3Lvwv7d8/VSnrg4g8kOklPAOpnyoOTfAlZ9ttcvnyZ9evX2V7XqVMXqFvQOiqllFJuIy4ujvXrv2D69BEAQa4u3yMXoLOOcRGRTrnkWw1cEJG7sqR/CtQXkQ7GmI7AJqCViPySKU9D4CDQXUTsxrnkZQE6pZRSytMVxgJ0HjPGpYASAGdjVPz5t0XF+jNrPv8s+5VSSilVyEp64LKPf8eyZJZ57MtvQIqTfE2BdPI39VoppZRSV6CkBy4rgcCMdVoAMMaEYBnfshJARC5hWb+lf5ZjBwJbRORcEdVVKaWUKvE8ZnCuMcYXywJ0ADWACsYY6zj1FSKSZIw5BKwXkUcBRGSLMWYV8LExZhSWFpTXgU2Z1nABmAisM8ZMw7I4Xa+MrUchX5ZSSimlMvGYwAWoBizMkmZ9XQ+Iw3K9XlnyDASmAh9iaYFaDjydOYOIbMoIgiYBT2BZ5+XerINylVJKKVW4PHJWUXHTWUVKKaWUzipSSimlVAmngUshiI+PJyQkhCjrsoJKKaVUCRIVFUVISAjk4bmB+aVdRYVAu4qUUkop7SpSSimlVAmngYtSSiml3IYGLkoppZRyGxq4KKWUUsptaOCilFJKKbehgYtSSiml3IYGLkoppZRyGxq4KKWUUsptaOBSCHTlXKWUUiWZrpzrZnTlXKWUUkpXzlVKKaVUCaeBi1JKKaXchgYuSimllHIbGrgopZRSym1o4KKUUkopt6GBi1JKKaXchgYuSimlcmSMwRhT3NVQCtDApVDoAnRKKaVKMl2Azs3oAnRKKU9ibW3RzwuVX7oAnVJKKaVKNA1clFJKKeU2NHBRSimllNvQwEUppZRSbkMDF6WUUkq5DQ1clFJKKeU2NHBRSimllNvQwKUQ6AJ0SimlSjJdgM7NlJQF6KwLgOtvkFKeTRegUwWlC9AppZRSqkTTwEUppZRSbsOjAhdjzLXGmG+MMUnGmOPGmAnGGK9cjokwxkg220uZ8s3JJk/Twr8ypZRSSgGULu4KuIoxxh9YA+wB+gINgLewBGdjcjj0f8BXWdLuAl4AVmZJ3wc8nCUtrkAVVkoppVS+eUzgAjwOlAP6ich5YLUxpgIQYYx5IyPNgYj8CfyZOc0YMxbYJyK/ZMl+QUS+d33VlVLq6mEdjJtbug7WVcXBk7qKegJfZwlQYrAEM7fmtRBjTGXgdiDatdVTSiml1JXypMClKZauHBsR+QNIytiXV/cA3jgPXK41xpw3xlwyxmwyxuQ5IFJKKfcj2WxKFR9PClz8gUQn6QkZ+/IqDPhJRA5mSf8ZGAn0Bu4DvLB0R7XNWsDx48cxxti2iIiIfJzePSQkJMAHH8DEiXzwwQeW10oppUq8iIgI2+cfEOTq8j0pcLlixpjqWLqVHFpbRGS6iLwvIutFJBboAhwDXs6aNygoCBGxbZ4UuIgIYydOJKh+ffjmG7h0iRFr1hBUvz5jJ07UPm+llCrhIiIibJ9/wHFXl+9Jg3MTgIpO0v0z9uXFACwLws7PLaOIJBljVmBpgSkxxk2aROQXX5C8ezcEWQLpCwDHjxPZ23IrJo4dW3wVVEop5dE8qcVlH1nGshhjagG+ZBn7koMwYJOIHM1j/hLV4ZuQkMCUyEiSli2zBS02QUEkLVvGlMhIEhMTi6V+SimlPJ8nBS4rge7GmGsypQ0ELgLrczvYGFMXaEceZxMZY8oBdwDbs+7bjqXZxtO2gIAAkjt0cAxarIKCSO7QAf+FC4u9rgXZlFJKXf08KXCZCVwCFhtjuhpjhgIRQGTmKdLGmEPGmNlOjg8DUoGFWXcYYyoaYzYaY4YZY7oYYwYCa7EMOppcCNdy9WrRIuf9118PJ08WTV2UUoVMw3x19fGYwEVEErAMmPUClgH/BaYC47NkLZ2RJ6sw4BsROeNk3yUgHssKvCuAKCwzmG4VEYfHQLch+0mE7rwBsGOHk9uTya5dRAUGFntdC3x9SimlrmpGZ4G4XkhIiGzb5hDPuD1jDPj4wG+/Oe8uOnYMGjYk4cQJKlWqVOT1U0oVDuuKufp5ofLLGLNdREJcWWaRtriY7NaRVu4jLQ1u6w3Hs8xwO34cunSBtDQNWpRSShWaPAcuxhiHAa7O0rI5tpkx5g3gj3zUzW3Fx8cTEhJCVFRUcVfF9VJS4PBd0KAZ3DEQXhxt+dmgGRw+bNmvlFKqRIuKiiIkJASgiqvLznNXkTHmJxFpnel1eWCriARnk78alhVmHwSCgQlArIgcuOJaX+U8uqsIsIwISQBigZNAINAf6wLF2pyslGfRriJVUIXRVZTrAnTGmPFkDHA1xqRnJAvwNzDdSf4wLMFKcywzdIYAC0WkZM2+8Xj+WN5apZRSqujkGriIyH+B/xpjIkUkPA9lfgasA24SkTgAY4yG6R5FhyopVZJoQ4u6muRncO4HxpjNxpgjAMaYVsaY/zrJ1wrLGmwbjTHfGmMexvn0Y6WUUkqpfMlP4PI+MIp/n8D8CxCaNZOI7BCR54A6wGtY1lbxN8bEGmP6X1FtVbHK/ODIs2fPQlQUTJhAVFQUZ8+ezfxQLaWUUqpQ5Cdw8RORzdYXYvmEupxdZhFJF5FVInI/UB3LonCPF7im6qogIkwcO5b6QUHcMWIEL40fz5oRI6gfFMTEsWM1cFFKKVWo8vN06IvGmIpkLDJqjGkJJOXx2LpAM6Bxfiqnrj6Txo3ji8hIdicnY1uC7sIFjgO9IyMBGDtxYnFVTymllIfLz3Tom4ApQCPgOywr2/fP3AqTJb9Oh/YwCQkJ1A8Ksg9aMjkONPPx4YiunKuUh8m8FIJSeVcY06HzteS/McYfaI+li2mziJx1kifrdOh5WKZD13NJjd1AiDGODzDyAB8Aa4D5OeQZCHTFTSdKazeXUtnQwEUVTLGs45JZxoMMV+SSTadDe6iTWJrbctIwI59SSqmSKyoqyrp6vMtXzs11cK4x5qeMn+nGmDQnW7wx5tlMh+h06DZtLN/ePWwLjIrioJ9fjpd+yM+PwKioYq9rgTallFIuMXToUDKGTJxxddlX/HRoY0xVYIuINMySXgpLr8GDQB9gFTBfRBZe0QndQEkd43IMuE7HuCjlURISEoiNDeDkSQgMjCI0NBR/f//irpZyE1fl06FFJB7o5SQ963To5cDQKz2fKj7+/v6EjxpFb19fsjwbmuNAH19fwkeN0qBFKQ+QeemDb0bApfHo0gfqqpCvMS7ZyW2mkIhcAOZkbMqNjZkwAYBmU6bQLTmZhli6h1alpREeHm7br5Ryb7r0gbpa5buryBgTALwL3A6kAUuBFzIG7lrztAV6AD1EpIPrqusePLWrKDNL83EsJ0+eJDAwkP79+2tLi1IeInO3cDnsnwMfClxElz5QeVPs06EzKjEXy/osXwKnsUyPPgfcDYwBBgBVscyb+15EbnJlhd1BSQhcPJLRh0cqBZalD1YD13tDpBd07wINm8PBHbDqWwhPg10plm+vbrn0gSoyBop3OnSG7sCLIjIFwBhjgChgK1AL+ArLlOiVIvKXqyqqlFKqaJwEznjDF/Vh97cQlGk0/vHj0Ps2qHAYTqYUWxVVCVaQwKUasMn6QkTEGDMReBQYJyKTXFU5pYqUDjZUCoBrpk3jx5dGcCBL0AKW18u+hcYN4K43p8EzzxRLHZWbKISW7IIOzk3L8vpYxs/VV1AXjxEfH09ISAhDhw5l6FCdSOU+tKtIKatOXRyDFqugIMt+eDZjUypbLl+ArqCByzRjzHZgb8b2W0Z6tk+LLkmqVq2KjnFRSrmrv/+GFs1zznP99XD+fNHUR7k1ly9AV5DAZTpwHZZBuNWwf3jFu8aYzcDPwC/APhFJv9JKKlU0tKtIKYDAwA9Ys2YEcCHbPIcP+9G161R0eK7Kmetbsq9o5dyMJ0Bfn7E1z/h5LVAuI0uSiJS/0kq6G51VpJRyZwkJCdSvH8Tu3clOu4uOHYPrrvPhyBGdDq1yVuwPWcxKRE4D32RsgG2WUSMsgcx1V1Q7pZRSRc7f35/w8FH07h3JsmVJDrOK+vTxJTw8XIMWVSxcsnJuZmJpwjmQscW6unyllFKFb8yYjFWym02hQzcvrm+YxJFDvqxalbFK9hhdJVsVjyt+yKJypF1FSilPkZCQQEBsLJw8SZSukq3y6arrKlJKKeXZ/P39YYhlAK4Ow1VXgyt+OrRSSimlVFHRwKUQWBegi4qKKu6qKKWUUkUuKiqKkJAQKIQF6HSMSyHQMS5KKU9iXYlDPy1UfhXGGBdtcVFKKaWU29DARSmllFJuw6MCF2PMtcaYb4wxScaY48aYCcYYr1yOqWuMESdbjJO8fY0xu4wxycaYPcaYgYV3NUoppZTKymOmQxtj/IE1wB6gL9AAeAtLcDYmD0WMAr7L9NruwVDGmJuARcAM4GmgFxBtjEkQkVVXfAFKKaWUypXHBC7A41iekdRPRM4Dq40xFYAIY8wbGWk52S8i3+ewfyywQUSezni91hjTDBgHaOCilPJICQkJkLEA3QeBgYSGhlrWdlGqmHhSV1FP4OssAUoMlmDm1isp2BhTFugMLMiyKwZob4ypeCXlK6XU1UZEGDtxIkH168M338ClS4xYs4ag+vUZO3EiOiNVFRdPanFpCnybOUFE/jDGJGXsW5bL8R8ZYwKA00A0MFpELmbsawB4A/uyHLMXS/DXGNh6ZdVXSqmrx7hJk4j84guSd+/G+pTFCwDHjxPZuzcAE8eOLb4KqhLLk1pc/IFEJ+kJGfuycwl4D3gU6ALMAp7A0pqSuWyclJ+QZT8Ax48fxxhj2yIiIvJQfaWUujokJCQwJTKSpGXLsHs0NEBQEEnLljElMpLExMRiqZ+6ukVERNg+/4Cg3PLnlye1uBSIiJwAhmdKWmeMOQXMMMa0EJEd+S0zKCiI48ePu6yOqmiY3LMoVTLExkK3bo5Bi1VQEMnduuG/cKHtOUZK2UREWDYAY1z+YehJLS4JgLOxJv782zKSV7EZP9tkKhsn5ftn2a+UUu7v5Elo1CjnPA0bWvIpVcQ8KXDZh2Usi40xphbgi+PYlNxIlp+/ASlZy894nQ4cyGf56iokuummGwJEBQbid/AgOfE7dIiowMBir6tuV/dWGDwpcFkJdDfGXJMpbSBwEVifz7JCM35uBxCRS8BaoH+WfAOBLSJyLv/VVUqpq1NoaChpq1ZBdl3ex46RtmoV/ftn/S9RqcLnSYHLTCwDbRcbY7oaY4YCEUBk5inSxphDxpjZmV5HGGPeMsb0yzhuAjAVWCwiOzOVPxHoZIyZZozpZIx5A8sidBOK4NqUUqrI+Pv7Myo8HN/evR2Dl+PH8e3Th1Hh4VSqVKlY6qdKNo8ZnCsiCcaYLsC7WKY+J2IJQCKyZC0NZH4MwD4sq+Y+hmXNlz+AN4FXspS/yRgTCkzCMuvoCHCvrpqrlPJEE8ZYFhyf0qwZyd26QcOG+B06RNqqVYSHh9v2K1XUjC4i5HohISGybdu24q6GUkpdsYSEBAIyVs6NCgykf//+2tKi8swYs11EQlxZpse0uCillHI9f39/25RnnfisrgaeNMblqhEfH09ISAhRUVHFXRWllFKqyEVFRRESEgJQxdVla1dRISg5XUXWJdv0d0gpT6Z/6aqgCqOrSFtclFJKKeU2dIyLKpCEhATrk+4JDPxAH3WvlFKqSGiLi8oXEWHixLHUrx9kfdI9a9aMoH79ICZOHKuPuldKKVWotMVF5cukSeP44otIdu9OzvT8tQscPw69e0cCMHbsxGKrn1JKKc+mg3MLgacOzk1ISKB+/aAsQcu/jh+HZs18OHLkhK7zoJQH0cG5qqB0cK7b2I7lT92zttjYALp1cx60AAQFQbduySxc6F/sdS3YppRS6mqngYvKM33SvVJKqeKmgUshiI+vTUhIG6KiZlH8DxV33RYYGMXBg345XvuhQ34EBkYVe12vngewK6VUyaML0LmZkjrG5dgxuO46HeOilKfRMS6qoHSMiypW/v7+hIePondvX2dPuqdPH1/Cw0dp0KKUUqrQ6HRolS9jxkwAoFmzKXTrlkzDhpbuoVWr0ggPD7ftV0oppQqDdhUVAk/tKsrMsnJuQMbKuVH6qHulPJh2FamCKoyuIm1xUQXi7+9vfdI9+rB7pTybBizqaqKBi7oC+t+ZUkqpoqWDc5VSSinlNjRwUUoppZTb0MBFKaWUUm5DA5dCEB8fT0hICFFRUcVdFaWUUqrI6cq5bqYkTIdWSimlcqMr5yqllFKqRNPARSmllFJuQwMXpZRSSrkNDVyUUkop5TY0cFFKKaWU29DARSmllFJuQwMXpZRSSrkNDVwKgS5Ap5RSqiTTBejcjC5Ap5RSSukCdEoppZQq4TRwUUoppZTb8KjAxRhzrTHmG2NMkjHmuDFmgjHGK5djbjDGfGSMOZRx3H5jzHhjjE+WfBHGGHGy9Sjcq1JKKaWUVeniroCrGGP8gTXAHqAv0AB4C0twNiaHQwdm5H0dOAg0ByZm/LwnS95zQNZAZe+V1l0ppZRSeeMxgQvwOFAO6Cci54HVxpgKQIQx5o2MNGdeE5EzmV6vM8YkA7OMMXVE5PdM+1JF5PvCqb5SSimlcuNJXUU9ga+zBCgxWIKZW7M7KEvQYvVzxs8g11VPKaWUUlfKkwKXpsC+zAki8geQlLEvP9oD6cBvWdIrGWPOGGNSjDE/G2P6Fbi2SimllMo3Twpc/IFEJ+kJGfvyxBgTiGVMzCcicjrTrkPA80B/LGNfjgOLnAUvx48fxxhj2yIiIvJ8EUoppZQ7i4iIsH3+UQg9Fx6zAJ0xJgV4TkSmZUn/E/hYRF7OQxllsAzwrQm0EZGEHPIaYDNQTkRaZt6nC9AppZRSugBdbhKAik7S/TP25SgjEPkYaAb0yiloARBLxLcYaJ7blGullFJKuYYnzSraR5axLMaYWoAvWca+ZGMalmnUt4tIXvIDSMamlFJKqSLgSS0uK4HuxphrMqUNBC4C63M60BjzEjAcuF9ENuXlZBktNPcAO0QkrWBVVkoppVR+eFKLy0zgaWCxMeZ1oD4QAURmniJtjDkErBeRRzNe3wtMBuYAx4wx7TKV+ZuIxGfkWw8swtJ64wcMAW4E7irUq1JKKaWUjccELiKSYIzpArwLLMMyw2gqluAls9JA5jEp3TJ+Ds7YMnsYS0ADlllFzwLVsUyV/gm4Q0RWXnntlVJKKZUXHjOr6Gqis4qUUkopnVWklFJKqRJOA5dCEB8fT0hICFFRUcVdFaWUUqrIRUVFERISAlDF1WVrV1Eh0K4ipZRSSruKlFJKKVXCaeCilFJKKbehgYtSSiml3IYGLkoppZRyGxq4KKWUUsptaOCilFJKKbehgYtSSiml3IYGLoVAF6BTSilVkukCdG5GF6BTSimldAE6pZRSSpVwGrgopZRSym1o4KKUUkopt6GBi1JKKaXchgYuSimllHIbGrgopZRSym1o4KKUUkopt6GBSyHQBeiUUkqVZLoAnZvRBeiUUkopXYBOKaWUUiWcBi5KKaWUchsauCillFLKbWjgopRSSim3oYGLUkoppdyGBi5KKaWUchsauCillFLKbWjgopRSSim3oYFLIdCVc5VSSpVkunKum9GVc5VSSildOVcppZRSJZwGLkoppZRyGx4VuBhjrjXGfGOMSTLGHDfGTDDGeOXhuIrGmI+MMQnGmHPGmHnGmMpO8vU1xuwyxiQbY/YYYwYWzpWo4hYXF4cxpsg3pZRSOStd3BVwFWOMP7AG2AP0BRoAb2EJzsbkcvgCoDHwGJAOvA58AdycqfybgEXADOBpoBcQbYxJEJFVrrwWVfwWLFjA5cuX8fb2Lu6qKKWUysSTWlweB8oB/URktYjMBP4LhBtjKmR3kDGmPdANeEhEFonI58D9wE3GmK6Zso4FNojI0yKyVkSeA74CxhWkshEREXbftIOCgrjnnnv47bffbHkGDx5sHZVd6OVYiQgtW7Zk7ty5trQ9e/bQpUsXfH19CQoKYty4caSlpeVa1qFDhxg2bBjNmzfHy8uLTp06ZZs3L+fYt28fN954IxUrViQsLIx//vnHbv+GDRuoUaOGQ/rw4cN59NFH83D1/3L3oGXOnDlOW3RmzpyZ67F5eS8WLlxInz59qFGjBuXLl6dNmzZER0fb9v/444+5ti6VKlWKv//+O0/X07lzZ3r37p3t/ttuu41GjRqRmpqap/KUUu7LY1pcgJ7A1yJyPlNaDJbWk1uBZTkcd0pENlgTRORHY8yRjH1rjDFlgc5YWloyiwE+MsZUFJFz+a1wxYoV+eqrrwA4fPgwY8eOpUuXLuzevRs/P78iLwcsLQ1nz57l3nvvBSAhIYGuXbty7bXXsmTJEn777TdGjhxJeno6kyZNyrGs3bt3s2LFCtq1a0dKSkq2+fJ6jsGDB9OwYUMmTJjAc889x+TJk5k8eTIA6enpPPPMM7z66quUL1/ervxRo0bRtGlTXnrpJRo2bJjrPdi5cyfNmzfPNZ87+PbbbylXrpztdf369XPMn9f3IjIyknr16jF16lSqVKnCihUruPfeezlz5gz/93//R8OGDdmyZYst/9tvv83y5ctZterfxklvb2+uueaaPF1HcHAwq1evdrpv7dq1rF27ls8++4zSpT3pvzSllFMi4hEbcBqIcJJ+AXguh+MWAOucpH8JfJnx72sBATplyXNDRvoNmdPbtGkjuRk/frxUrlzZLm3jxo0CyIIFC0RE5KGHHpLcynJVOVYdOnSQl19+2fZ68uTJUqlSJTl37pwt7fXXX5dy5crZpTmTlpZm+/c999wjt956q9N8eTnH33//LYCcPn1aRERiYmIkJCTElv+DDz6QG264QdLT052eo0uXLhIeHp5jfa3efPNNuXjxol3aypUrpWPHjlKxYkW55pprpE2bNrJq1ao8lVccPvroIwHk77//ztdxeX2/4+PjHY4dNGiQ1K1b12m5vXr1kvbt2+erLpm9/fbb4uXlJZcuXXLY17FjR2nRokW2771SqvgA28TFn/ee1FXkDyQ6SU/I2Hclx1l/Zs2XkGU/AMePH7drEo+IiMjh9P9q06YNYBkYmtkXX3xB06ZN8fHx4aabbmLPnj2FUs6hQ4fYvHkzoaGhtrSVK1fSvXt3KlT4t7ctLCyMixcvsn79+hzrUapU3n698nKOy5cvA9haD3x9fW1p58+fZ8yYMUyfPj3bAa733HMP8+bNIz09Pdf6JCcn4+PjY3u9bt067rrrLm6//XYWL17MZ599RteuXT3y231e3+8qVRzXlGrVqhXHjx93Wu6uXbu47rrrnO5LTk5m/PjxNGzYkHLlyhESEsKmTZvs8gQHB5OWlsahQ4fs0r/66iu+++47XnnlFR3crNRVIvMQBiDI1eV7UuBy1QgKCrKLDvMauFgDjcDAQFva77//Tnh4OGPHjuWzzz7j3LlzdO/eneTkZJeX88033+Dn50eLFi1safv27aNp06Z25deuXRtfX1/27duXp+vKTV7OERAQQN26dXnnnXc4e/Zs5lUZmThxIl27dqV9+/bZnqNDhw6cOnWKXbt25ViXgwcP0rhxY7u0efPmcccddzB+/Hhuu+027rzzTl577TU6d+5ckMvNExEhNTU11y03DRo0oHTp0jRp0oRZs2blmv9K3u8tW7Y43DuAc+fOcfToUa6//nqHfampqdxxxx188MEHPPfccyxdupQGDRrQq1cvTp06ZcsXHBwMwP79++2OHzt2LB07duSOO+7I9dqUUkUjIiIic8+E828zV8CTvjImABWdpPvzb8tIdsdVzeU468+s5ftn2Z9v1g+fw4cP8+STT3LNNdfQteu/Y4LPnDnDkiVL6NChA2BpTWnQoAFz5szh8ccfd2k527dvJzg42K6lJCEhgUqVKjnU29/fn4SEAl+2nbyeY8aMGfTv35+XX36ZRo0a8d5773Ho0CH+97//5RqQNGvWDC8vL3788Ue7wOyff/6xGxPz5Zdf8thjj9kdW65cOb7++mumTZtGWFiYXUBYWObOncvDDz+caz7JZuXr6tWrM3HiRNq2bUtaWhoxMTE8/vjjJCUlMWLEiGzLK+j7/c033/DFF1/w4YcfOuz79ddfAZy2uLz11lv8+OOPbN++3Rb0dOrUidq1axMbG8tTTz0FQI0aNahQoYJd8LRkyRK2bdvGhg0bHMpVSnkuT2px2QfYfVU0xtQCfDP25fm4DE0zHfcbkOIkX1Ms06cPFKC+/PXXX3h7e+Pt7U2TJk04fPgw8+fPp3r16rY81apVswUbAHXq1KFNmzb8+OOPLi/n5MmTTrsAcpO1dSAvM44KomfPnpw+fZr9+/ezd+9eateuTXh4OCNGjKBmzZq899571K5dm9q1azNjxgy7Y0uXLk2lSpU4efKkLW3atGn06tXLLl/WQAYsLTqDBw9mwoQJ1KhRg+7du7N7927bfmMMgwcPtr1et26d3Qyq3PY707t3b7Zu3Zrrlp3u3bszZswYunXrRs+ePZk7dy4DBgxg0qRJeeouy4+4uDjuvfde+vbta3edVtbAJWuLS3p6Om+99RYPPvgg9evXt/3+GGNo2LAhf/zxh13+pk2b2lpcRIRx48bRs2dPbr75ZpRSJYcntbisBJ4zxlwjItY5lgOBi0BOgzFWAmONMTeJyCYAY0wIUD9jHyJyyRizFugPZG5vHwhskQLMKALLbKA1a9ZgjCEwMJCgoCCHfvpq1ao5HFetWjVOnDjh8nKSk5Px9fW1y+Pv78+5c46Xl5CQgL+/pcFp/fr1dt0mt956K+vWrcvhyu3l5RxWvr6+tm/mq1evZseOHcyfP58dO3YwduxYNm/eDED79u256aab7GYHlS1b1q5r7Nlnn+WTTz7h+PHjBAUF8ccff1C7dm2HelSsWJF3332X6dOn88033/D444/z6KOP8v3339vybNy4kYMHD9KoUSOn15jb/qwCAgKoWNFZA2LBhYaGsmDBAuLi4rKdXZSf9wLg7Nmz9OzZkzp16jBv3jynZe7atYv//Oc/DkHxrl27iI+PZ8aMGQ6BJkCPHj3sXgcHB9taXBYuXMiuXbuYM2eO03MqpTyXJwUuM7FMV15sjHkdS+ARAURKpinSxphDwHoReRRARLYYY1YBHxtjRvHvAnSbRGRNpvInAuuMMdOwLE7XK2Oz/981H0qXLp3r+iqnT592mtasWTOXlxMQEGDXIgGWb7lZxzYcPXqUpKQk21iINm3a2H37z+sU1/ycI6u0tDRGjBjBG2+8Qbly5Vi3bh233XabLX+XLl1Yv369XeCSmJhIQECAXTmhoaEsXryY4cOHs3z5cgYNGpRtPb28vOjWrRvdu3cn60M0n3/+eSZOnMjHH3/s9Njc9md1pV1FzliD2ZwGsebnvUhKSuLOO+/k8uXLLF++3CHotfr111+djm+xBs1r1qxxGqTVqVPH7nVwcDBLliwhLS2NiIgIBgwYQKtWrbK9FqWUZ/KYwEVEEowxXYB3sazZkghMxRK8ZFYayPoYgIEZeT/E0n22nCxrtojIJmNMKDAJeAI4Atwrhbxq7unTp9m8ebOtm+ePP/7gp59+ytOHWn7LadKkid3aG2DpnnnzzTf5+++/bQHJ/PnzKVeuHLfeeitgCVTyusCdM3k5R1bvv/8+/v7+DBz471MXkpKSbP++cOGC3Yd6fHw8SUlJDoNHQ0NDGTp0KMOHD+fs2bMOrQrDhg3D29ubTp06UblyZTZu3Mjs2bOJioqyy/fII4/w5ptvcuCA817D3PZnZe0qcqXY2FiqVKniEBBkltf3IjU1lf79+3Pw4EE2b97stEXP6tdff+WBBx5wSLd2ZVaqVMk2Ey4nwcHBJCYmEhkZycGDB1myZEmuxyilPJCr51frVvB1XLJ66KGHpEqVKlKvXj2ZN2+eLF68WK6//noJCgqyrTPiqnJERL7++mu7tVJERM6ePSuBgYHStWtXWb16tcyaNUv8/Pxk9OjRuV7jhQsXZOHChbJw4UJp166dXHvttbbXFy5cKPA5/vrrL6lSpYps377dlvbTTz9JmTJlZPbs2fK///1PypQpI7/88ott/4oVK8QYI2fOnHEor3Xr1rJz506Jiopy2PfWW2/JDTfcYFu/pV27dvLFF1/Y5bH8GVnWkrnvvvtk7dq1dmvW5La/MPTr109ee+01WbFihSxbtkzuv/9+AeTtt9+25Zk7d654eXlJXFycLS2v78WQIUMEkOnTp8uWLVvstuTkZFu+48ePCyCzZ892qGNKSoq0bNlSGjduLHPnzpX169dLbGysPP/88zJnzhyH/AcOHBBAfHx8ZMiQIa64TUqpQkYhrONS7B/ynri5MnBp06aNLFq0SBo1aiRlypSRDh06yK5du1xejojIpUuXJCAgQD7++GO79N27d0vnzp3Fx8dHAgMDZcyYMZKamprrNR45ckSwLNDnsB05cqTA5/i///s/eeSRRxzSp02bJoGBgRIYGGj3AS0i8vTTT0unTp2cljdhwgRp27atnDx5MtdrcsYamFy+fFkaNmwo77//vtPAJbv9heGll16Sxo0bS7ly5cTHx0dat27t8L5aF6kryHtRp06dPL231mD4hx9+cFrPo0ePyr333ivVq1cXHx8fqVevntx3332yf/9+h7ypqalSpkwZ8fHxkaNHjxbsxiilipQGLm6y1a5dW9q0aSOzZs3K9s28Wj399NPSq1ev4q6GS6WmpkqNGjXkk08+cbp/z549cvPNNxe4fGtgIiIye/ZsqV27ttPAJbv9SinlaWbNmiVt2rQRIE5c/BlrRPI+uE/lTUhIiGQdvOku/vzzTxo3bswvv/zidDExdxQTE8PYsWPZu3dvtqvd7t69226gcn4YY7D+HaWmptK0aVNq1qxpm1mV236llPJUxpjtIlLwQZBOeMzgXOUaNWvW5MMPP+TEiRMeE7iICLNnz85xif6CBi3W8q1Kly7tsCx9bvuVUkrlnba4FAJ3bnFRSimlXKUwWlw8aeVcpZRSSnk4DVyUUkop5TY0cFFKKaWU29DAxQMcOHCAiIgILl++bJf+4YcfUrFiRduTo10pJiaGsmXLkpKS4rIyR4wYgTHGtv3nP//hwQcf5K+//nLZOTLbs2cPXbp0wdfXl6CgIMaNG5fvB0QeO3aM8uXLY4zhn3/+ydexIkLLli2ZO3euLS0uLg5jDHFxcfkqq6AK83wLFy6kT58+1KhRg/Lly9OmTRuio6Md8g0fPpxHH33U5edXSnkmDVw8wMqVK5kxYwZlypSxS+/atSvfffddjrNpCmrHjh1ce+21eHt7u6zMXbt20a5dO7Zs2cJ3333HyJEjmTdvHk899ZTLzmGVkJBA165dMcawZMkSxo0bx1tvvcX48ePzVc5zzz3n8DTpvFqwYAFnz57l3nvvLdDxV7vIyEjKly/P1KlTWbp0KZ07d+bee+/lnXfescs3atQo5s2bp7OtlFJ54+qFYXQr+gXoHnnkEencuXORnMuqZ8+e8vDDD7u0zGrVqslTTz1ll9a7d28JCAhw6XlERCZPniyVKlWSc+fO2dJef/11KVeunF1aTtavXy/+/v7y5ptvCiB///13vurQoUMHefnll+3SrKsNZ13NtrC44nyHDh1ymh4fH++QNmjQIKlbt65DepcuXSQ8PLzAdVBKXV0KcwE6bXEpBFWrVmXbtm0MHTo0x3zz5s3jlltuoXLlypQvX5727dvbPeQwPT2d8uXLM3PmTMaMGUPNmjWpWLEiQ4YMIT09HYBq1arx4YcfsnbtWlsXy7fffms79oMPPiiUa/zll19o2bKl7fW3335L5cqVefbZZ/Pd3QJw6tQpTp8+TXBwsF16tWrV8PLK+kzMK7dy5Uq6d+9OhQoVbGlhYWFcvHiR9evX53p8Wloa//d//8e4ceOoUqVKvs9/6NAhNm/eTGhoaK55W7Zs6fR36f777+fmm2/O97ldITk5mXnz5nHbbbfRqFEjp3mc3ZdWrVpx/Phxh/R77rmHefPm2X6vlVLubejQoWQsC3LG1WVr4FKMdu/ezUMPPcTChQtZuHAh1atX56677rKNVTl8+DAXLlxgypQpXLx4kY8++ohRo0bxv//9jyVLliAiLF26FB8fH0aPHs2WLVvYsmULHTp0sB3bvHlzl9c7Pj6eEydO2AKXd955hzvvvJPXXnuNadOmFSjQ2LVrFwBNmza1paWnp7NlyxZ69+7t9BgRITU1NdfNmX379tmdC6B27dr4+vqyb9++XOs7c+ZMLl26VOBurG+++QY/Pz9atGiRa94bb7yRrOsC/fjjj0RHRxMZGVmg8xfUzz//zPDhw6levTqPPPIIlStX5ssvv8zz8Vu2bHG6sGGHDh04deqU7fdAKaWy5eomHN3y9pDFrFJSUmTfvn0CyJ49e0REZPHixQLIW2+9ZZe3evXq8vrrr4vIv039W7ZsscuzePFiKVWqlPzzzz/5rktuVq9eLcYYiY+Pl8cee0wCAgJk7dq1V1RmZGSkABIXFycpKSny+++/y2OPPSbt2rVz+kRnkX8fEpjb5kzp0qVl6tSpDuk1atSQl156Kce6njlzRvz9/eXLL7+0q0d+uoqGDBkiISEhDunOum5mz54t3t7edk/y7tChg9x33315Pl928tJVlJiYKO+99560bt1aAGnZsqVMnz5d/vrrr3yda82aNWKMkY8++shhX0pKinh5eTl9QrdSyn1RCA9Z1CX/i8nFixd55513+OSTTzh8+DBJSUm2fb6+vgDs3LmTChUqMHz4cNs+ESExMdHWDL9z506MMVx//fV25e/cuZP69evj5+fncO5z585x4sSJXOuYtUXC6pdffiEgIIB+/foRHx/PDz/8QMOGDW37IyIiiIiIyLX8rPUFqFu3ri2tZs2abN26lcqVKzs9pnfv3mzdujVf53GF0aNH065dO3r16lXgMk6ePJnnLqZ27dqRkpLCL7/8Qrt27Zg/fz4//fST3Qydli1b8ssvv+RYTkHel6+++op+/fpRrlw57rvvPmbPnm3XRZhXcXFx3HvvvfTt25fBgwc77C9dujSVKlXi5MmT+S5bKVWyaOBSDESEXr16cfDgQcLDw7n++uvx9/fnww8/ZO7cudSqVQuwdJ/cfPPNdrOFfvvtNy5evGh7tk52AcquXbuy7SZauHAhQ4YMyVM9ndmxYwciwsaNG1mzZo1d0ALw3//+N98fkNYZRe+88w4pKSmsXr2a8ePHM2nSJN59912nxwQEBFCxYsV8ncfK39+fc+fOOaQnJCTg7++f7XG7d+/mww8/ZMOGDSQmJgLYgs5z587h5eVFuXLlcj1/cnKyLUDNTdOmTalQoQJbt26lZcuWvPDCC4SHh1O7dm1bntyCFijY+1K2bFl8fX25ePEi586dIzExERHBGJPnMs6ePUvPnj2pU6cO8+bNy/FcycnJ+aqfUqrk0TEuxWDDhg2sW7eORYsWER4ezu23305ISAh79+7l+uuvp1Qpy9uya9cuhzEQO3bsoFSpUlx33XWAJXBxNk5i165dDq0wVo899liemuOys2PHDp5++mnuv/9+Hn74YeLj4237RowYAVhaAPI6cDQ9PZ09e/bQoUMHQkJCaN++PePGjaNHjx4sWLAg2wGbc+fOxdvbO9fNmaZNmzqMZTl69ChJSUnZtjQBHDx4kJSUFNq3b4+/vz/+/v62cS41a9bk//7v//J0zQEBAbbAJzelSpWibdu2bN26lalTp5KcnMyLL75olydzILFx40Zuvvlm2rRpwy233MKePXsK9L4AdO7cmWPHjjF79myOHTvGbbfdRoMGDZgwYQK///57rscnJSVx5513cvnyZZYvX55jsJaYmEhAQECe66aUKqFc3fekW+5jXD799FMB5OzZs7a0lStXCiBDhw4VEZGkpCQpVaqULFiwwO7YcePGScOGDW2vr7/+eocptdZjFy1alGM9CiI5OVm8vb1l8eLFkpSUJK1atZJbb71VUlJSbHnIZlxJdqxjez755BO79I8//tjp+B2rM2fOyNatW3PdnJk8ebL4+/vL+fPnbWlvvvlmrtOh4+PjZe3atXbbCy+8IICsWLFC9u3bl6drHjdunNSuXdshPbsxJ6NHj5aaNWvKNddc43SavfWe//XXX9K2bVvb79amTZukXbt2dnnycr7sHD58WMaMGSM1a9YUY4x06dLF4X2zSklJkV69ekmVKlXkwIEDOZZ7+vRpAWTp0qV5qodSyj1QCGNciv1D3hO33AKXPXv2SKlSpeSBBx6QNWvWyKuvvipNmzYVY4y89957IiKydetWAWT//v12x951113Sr18/2+uOHTvKTTfdJOvXr5ctW7ZIenq67diDBw/mWI+C+Omnn+zKjouLkypVqsjw4cNteZx9QK5du1YAp4N4Fy5cKIDs2rXLLj0+Pl5KlSol48ePd+k1iIicPXtWAgMDpWvXrrJ69WqZNWuW+Pn5yejRo+3yzZ07V7y8vCQuLi7bspwNzs3pekVEvv76awHk9OnTdunZBRJLly4VQK677jpJTU11KM96z5cvXy7+/v7SokUL21arVi27PHk5X25SU1Nl+fLlctddd0np0qWd5hkyZIgAMn36dNmyZYvdlpycbJd3xYoVYozJdiC2Uso9aeDiJlteFqB7//33JTAwUCpUqCADBgyQjRs3CiCbNm0SEZEPP/xQfH19JS0tze64+vXrS0REhO31xo0b5dprr5XSpUvLf/7zH9uxfn5+kp6enu35C+qjjz5yKPubb74RLy8v22wRZx+QX375pQCye/duh33jxo0THx8fu1Ybq/bt28uNN97ougvIZPfu3dK5c2fx8fGRwMBAGTNmjENQYA1Kcvpgdxa45HS9IiKXLl2SgIAA+fjjj+3SswsktmzZIoB89dVXTsuz3vOlS5fKHXfckWOevJwvP06ePOk0vU6dOtnO9Mp6vqefflo6depU4Doopa4uhbkAXbF/yHviVpDp0J7kmmuucehuGTduXIn6YMrL9T799NPSq1cvu7TsAonevXtLz549sy3LGpScOXNGAgMD5aeffhIRkbS0NNm+fbuIOH9finqlXmdSU1OlRo0a2XY5KaXcV2G0uOisIuVyzzzzDG3btqVq1aps3LgRgM2bNxMeHl7MNSs6ebne5557jsaNG3PgwAGni7IlJyezc+dOYmNj+eabb9ixY0eu561cuTILFy7kqaee4sKFC1y+fJnQ0FBat27t9H25GixcuJBy5coRFhZW3FVRSrkBYwmIlCuFhIRI1pVOlXImJiaG6tWrc+uttwKW9U7q1avHkSNH2L9/Pz179qRevXpMnz6dO++802kZ8fHxNG3atEBP0c58vsxr6BSl6OhoatSowS233FIs51dKFR5jzHYRCXFlmTod+iowbNgwypUrR1xcnEeer7AcPXqULl26EBwcTLNmzXj++ecpykD83XffpWnTpnTu3BmgQAuzhYWF2YKWrLp37056ejq//fZbtkHL/v37adeuHS+//HK+z321GDRokAYtSqk80xaXQpCfFpdt27YRGhrKQw89xI4dO/jiiy8KtW55PV9qaiqlS1/dPYknTpzg2LFjhISEcPnyZW6//Xaefvpp7rnnniI5f3BwMJ9//nmO674opVRJpi0uHkZEGD58OFFRUYwfP54zZ87w9ddfF9v5jDGMHz+e1q1b88orr9h1Haxbt45OnTrZ5Z08eTKtW7cmODiYH374Icdzb9iwgbZt2+Ln54cxhvj4eFq2bMncuXMd8qakpDB16lTatGlD+fLlqVChAt27d3c4R/Xq1QkJsfw9lClThubNm3P06FHb/uHDh/Poo4/m5xbl2YMPPshvv/1Gv379eO6554B/F4HbtWsXDRs25OzZs4Dl0QQzZ84slHrs2bOHLl264OvrS1BQEOPGjcv16dyxsbF06NCBypUr4+PjQ5MmTZg0aZLt4Z5Wc+bMsT1xPPNWkGvJ63t6NVi4cCF9+vShRo0alC9fnjZt2tg9XiGztLQ03n//fdq1a8c111yDj48PjRs3Zvjw4Rw4cCBf5+3cuXO2DxQFbE/izu7BoUqVGK4e7aub+84qAmxTuI8cOSJ16tSx7Vu7dq3ceuutdnmtU3k/++wzuf3227Mt99KlS1KlShW57777ZMOGDfLDDz9ITEyM1KpVSy5fvmyXNyEhQW644QapXLmyTJgwQVatWiXR0dHStm1bKVOmjKxfv97pOc6cOSM1a9a0PaDSeg1ly5YtlPVsRCzTfTPPxiHTdOOpU6dKaGiovPvuu9KnT59COf/Zs2elevXq0qVLF1m1apW8//774uvr67AWTVYzZ86U0aNHy+LFi+Xbb7+V1157TXx8fOSpp56yy2ed5v3tt9/arcFy6tSpfNWzoO9pcWnXrp0MGjRI5s+fL998842MHDlSAHn77bft8v3zzz9yyy23SLly5WTkyJHy5ZdfyjfffCNvvPGG1K5d2+lDNHPyxBNP2C0umdm3334rgHz22WcFvi6ligO6jot7bO4cuFifJp2XwMW6iNihQ4ekQYMG2ZZrXWwt83ofHTp0cFjxNz09XTp16iS1atWSw4cP2+27ePGi1KpVS5zd2+TkZOnUqZNMmTLFYV+XLl0kPDw8+4u+AjkFLunp6dKlSxcJCgqS+Pj4Qjn/5MmTpVKlSnZTnF9//fVcV/915uWXX5aKFSvarc9TkKdeZ1XQ97Q4OXu/Bg0aJHXr1rVLu/POO8Xf31927NjhkP+ff/6RuXPn5uu8b7/9tnh5ecmlS5cc9nXs2FFatGhRKGszKVWYCiNw0a6iQhAfH09ISAhRUVE55jt37hzDhg0jICCASpUqER4ezquvvurw0EJXyOu5rA9rLF26tN0zgpw9/K5s2bIAeHl5Zdt8fdNNN9G9e3cAAgMDMcawYcMGNm/eTGhoqF3e2bNns27dOmbPnk29evXs9vn4+HD//fezfft2u2cjpaWlcd9999GqVStGjhzpcP577rmHefPmZfu8I2fOnj3LsGHDqFKlCgEBAbzwwgu2a3nnnXfyVMbff//N0aNHKV26dJ6fSZRfK1eupHv37lSoUMGWFhYWxsWLF1m/fn2+yqpcubJDV5ErFOQ9LW7OntrdqlUrjh8/bnu9aNEili9fTlRUlNOHmfr5+fHggw/apSUnJzN+/HgaNmxIuXLlCAkJYdOmTbb9wcHBpKWlcejQIbvjvvrqK7777jteeeWVfD3cUqniFBUVZe3Kd/yDulKujoR0y1uLy+XLl6Vdu3Zy/fXXy/z582X58uXSsmVLqV69ut2S/lbp6emSkpKS63Yl5yJTi8Hly5elYsWKcvToUUlPT5eHH37YocXFKmvrTGY///yz3H777XLTTTfJli1b5Pvvv5eZM2eKn5+fw6rATZo0kdatW2d7z2bMmCGA3TfcRx99VAYPHpztN9FffvlFAPnll1+yLTezhIQEadKkibRp00aWLFkir7zyigAyY8YMCQwMlIsXL9ry5tTicv/998ukSZNk6dKl0rZtW6fvzZW8pyIiVatWdfo4BF9fX3njjTdyvdbU1FS5cOGCbNy4UerUqePQMmVtcalWrZp4eXlJ48aNZebMmbmWm1lB3tMrdaX31Zm7775brrvuOtvr9u3b273OTUpKitx2221SvXp1mTlzpqxatUoGDBgg11xzja0l8s8//xRAFi9ebHdsSEiIdOzYMV/1VepqgXYVuceWl8Bl8uTJUqVKFfnrr79sad99950Adkv6W1k/RHLbruRcWY+fNWuW1KtXT2688UZ57rnnChS4iIi0aNFCnn/+edvrIUOGOPT/79+/XwB56623si1nwoQJds9J2rRpk+35Pdbn8kyfPt3umJSUFPHy8pKoqKhsy81s1KhRUrFiRbvugipVqkjFihUduqKyC1yio6Plpptusj0+4PHHH3foFhO5svdURKR06dIydepUh/QaNWrISy+9lOu1li1b1naOBx980CGQ/Oqrr2TixIny9ddfy4oVK+TBBx8UQCIjI3MtW6Rg76krXOl9zWrNmjVijLE90uLkyZNijHEIGrMGTJmD6ddee03Kly9v9+yxy5cvS2BgoLz77ru2tAoVKsjkyZNtr7/44gsBZMOGDfm8C0pdHQojcLm657t6qPT0dKZNm8ZTTz1FQECALb1+/foAtGjRwuGY3r17s3Xr1kI9l+V37F9Dhw5l6NChTsvNnLdu3brZrgmTmprK3r17bTNvAE6ePOnQHP/rr78CcO2112Z7Lfv27cPHx8c226ljx44Odc6qdOnSVKpUiZMnT+aYDyzXNHfuXIYMGWJXv4oVK3Lu3Dkef/xxu/xZr9lal7CwMLtVYN9//32n5yvoe+oqmzdvJikpiR9//JEJEyYwfPhwZsyYYdvfvXt3WzcfQM+ePUlOTmbSpEk888wzlCqVc09zQd5TV3DlfY2Li+Pee++lb9++DB48GIADBw4gIg7r9rz44ou88cYbttfz589nwIABpKen89Zbb/Hggw9Sv359W7eqMYaGDRvyxx9/2I5p2rQp+/fvByy/T+PGjaNnz57cfPPNLrkepTyBBi7F4Ndff+X06dN07drVLv3PP/8EnAcuAQEBVKxYsUjO5Ur79u3j8uXLduMAkpOT8fX1tcv3999/A87HFwBcuHCBpUuX0qNHj3yvL1O2bFmnY3Sy2r9/P/Hx8dx+++126WlpaTzzzDO28T+uUtD31Mrf359z5845pCckJODv75/r8a1btwYsY3eqVKnCQw89xMiRI2nQoEG2x4SGhrJgwQLi4uJswW92CvM9zcmV3lers2fP0rNnT+rUqcO8efPs0gG7sUUAjz/+OP3792fdunU899xztGrVCrBMj4+Pj2fGjBl2gaFVjx49bP8ODg5m3759gGVa9q5du5gzZ84VX4tSnsSjBucaY4YYYw4aY5KNMduNMV3ycMwwY8xqY8wpY8w5Y8x3xphuTvLFGWMky5b713gnjh07BkDNmjXt0tesWUOFChWcfvucO3cu3t7euW6uOJcr7dy5E29vb7tF2gICAhwGrFavXh3Abh2WzGbMmME///zjdABubhITE+1am7Jz5MgRAOrUqWNL27x5M3FxcQVaFTc3BX1PrZo2bWr7kLM6evQoSUlJ+V4UzxrEWO9BdqyDQ/MySLQw39OcXOl9BUhKSuLOO+/k8uXLLF++3C7Qrlq1KuB4XfXq1SMkJISkpCSuueYa28D3EydOAJa/ua1btzpsmVs1g4OD2b9/P2lpaURERDBgwABbAKSUsvCYFhdjzCBgJhABbAIeBpYbY24QkV9zOHQ08BXwHnABuB/4yhhzl4gszZL3MyDztJICTcOoXLkygN3zYc6cOUNkZCTXX3+90w+FgjZ/5+dcxhi7rpe4uDg6derE9u3b6dLFEgOePXuWf/75h9q1awPw+uuv06NHj2y7bHbs2EHTpk3tPiiaNGnCli1b7PJ16NCBChUq8OGHH3L33Xfb7Vu/fj3jxo1jyJAh3HTTTfm6/vj4eJKSkpw+xDArLy8v2zWCpan++eeft/3b1a60S6Nnz568+eab/P3331xzzTWApXuiXLly2T5GIDvfffcdgMPMn6xiY2OpUqWKXXCXnfy+p0lJSdx3330cPnyYtLQ0+vbtS9WqVTlz5gyTJk3if//7HyNGjCAhIQEvLy9q1arFgQMHHFrvrvS+pqam0r9/fw4ePMjmzZupVq2a3f42bdoQEBDAhx9+yEMPPeRw/Pbt22nZsqXtb8sawFWqVIk2bdrkeO7g4GASExOJjIzk4MGDLFmypMDXoZTHcvWgmeLagP3Ah5lelwJ2AZ/mclwVJ2mbgbVZ0uKAKXmpS26Dc5OTkyUoKEiaN28uy5cvl5iYGGnVqpUEBATIk08+meOx+ZWfc5FlwKKzQbcfffSRPPTQQzkel1mPHj3kvvvus0uzruty+vRpu/T3339fABkwYIAsXbpUli9fLk8//bSULVtWBgwY4LBYXV6sWLFCjDFy5swZW9ratWsFkLVr19rlPXPmjPj4+MhNN90kK1eulMGDB0twcLA0atRI7r33XruBuFeDs2fPSmBgoHTt2lVWr14ts2bNEj8/P7sF6ObOnSteXl4SFxdnS+vevbu8+eabsmLFCvn6669l3Lhx4ufnJwMHDrQrv1+/fvLaa6/JihUrZNmyZXL//fc7LMSW3b20ys97umjRIhk6dKjd9f3888+2GTX33XeftG3bVrZs2SK7du0qtJk2Q4YMEUCmT59ut/Deli1bbGsX/e9//xNAunTpItHR0bJu3TqJiYmRhx9+WHx8fGTkyJG28lJSUqRly5bSuHFjmTt3rqxfv15iY2Pl+eeflzlz5tid+8CBAwKIj4+PDBkypFCuT6mihM4qyjb4qI9lpkCPLOnjgL8KUN57wP4saS4LXERENm/eLM2bN5dy5crJjTfeKOvXr5eKFSvme7ppXuT1XIURuAQFBcnrr79ul3bp0iUJCAiwrbyb2aJFi6Rdu3bi6+srgJQpU0bmzZuXh6t07umnn5ZOnTrZpX355ZcCyO7dux3yz58/X+rWrSve3t7SunVrOXz4sMyZM0f8/Pxs13HhwgW56667pHnz5tKsWTOnM4aKyu7du6Vz587i4+MjgYGBMmbMGNtsJpF/Z9hkDrrGjBkjzZo1Ez8/P6lYsaK0atVK3n77bYfA8KWXXpLGjRtLuXLlxMfHR1q3bu3wnuV0L63y+p4eOnRIGjZsKKNGjZIvv/xSUlNTJS0tTQIDA+XChQvSvHlz+fTTT2Xy5Mnyzjvv5LpCcEHVqVMn25lIme/jkiVLpGPHjuLr6ytly5aVevXqyaBBg+Tzzz93mKF19OhRuffee6V69eri4+Mj9erVk/vuu89ulpGIZYp6mTJlxMfHR44ePVoo16dUUdLAJftAo1fGfyx1s6T3z0ivms/yfgKWZUmLAxKwdA+dA2KBOs6OL8hKoIcPHxZAtmzZku9jXXUuwDatuEWLFhIcHHzFgUt2nn76aenVq1eOeb766ish06MF8is1NVVq1Kghn3zyiV36uHHjHIKZ/HDWMlBS5fde5vaeJiYmysKFC2XgwIHSrVs3ERHp27evvP/++9K3b185fvy43H777XLPPffI6tWrXXINSqnCo4FL9oHGfRkBSqUs6V0z0hvno6xHMo7pnCV9OjAIuBkYChwD/gAqZi2jevXqdt/SnC0SltXixYulVKlStiX3C1N25yqMFpfsHD16VMqVK+fwjTOre++9VypXruzQrZQX0dHR0rBhQ4fFxrp27SpLly7Nd3lWzloGSqqC3Mvs3tOjR4/KhQsXRETk2LFjUqlSJRGxPPepdu3atjVrWrVqJUFBQba8Sqmry/jx4zN/Bh4XF3/mX7WDc40xFYHqueUTkX255cnHOdtgGXw7XUTWZjnPM5lebjTGbAZ+wTIIeFrmvEFBQXbLg+fFjh07qF+/vsun3Bb3ubJTs2ZNPvzwQ06cOJHjwNnM01DzS0SYPXu2w1Tb1atXF7hMgAYNGrBt2zZWr17Nxx9/zPTp0wv1qd5Xs4Lcy+ze0127dvHiiy/aBom/++67AHTq1IkRI0bYnk7etm1bfv31V4dBuUqpq0NERAQREREAGGPy92GYB0bE9bMlXMEY8xjwQW75RMQYY3oBX2LpKvo9Uxn9gQVANRHJ8WEoxpj6WAbl/gjcLSJpeajjbmC7iNg9lCQkJES2bduW2+FXnexmFWVeaG3OnDmsW7fObm2JrMd5uj///JOAgAB8fX05fvw4zZo1IyEhobirpZRSVx1jzHYRCXFlmVdti4uI/A/4Xx6zW1tdmgK/Z0pvCpzNQ9BSDfg649iwvAQt1mpmbCVe5nViOnfuzCeffFKMtSlc2bUMKKWUKnxXbYtLfhlj9gMbReSxjNelsHTl7BSR+3M4rjywDqgIdMgtyMl03HUZ5YeLyNuZ97lri4tSSinlSiWqxaUAIoBPjTFxwHfAQ0Aj4F5rBmPMrcA3QBcRWZ+RvBhoDgwGGhhjbOudi8j3GcfdgWVhuuXAcSwtOWOwDM6dU3iXpJRSSqnMPCZwEZHojNaTF4CxwG7gTrFfNdcAXhk/rawPpnE2YtCa7yhQDcsg3ErAX1hW231ZRM676BKUUkopTxPk6gI96llFIvKBiDQUkbIi0lpEvsmyf52IGBFZlynNZLdlyrNTRLqISFUR8RaRQBEZLCJOR0vHx+ept0ldZaKiooq7CiWC3uf80ftlUVLugwdeZ66zg/PLowKXq8WZM2eKuwqqADzwP4yrkt7n/NH7ZVFS7kNJuc4r4TGDc68mxpi/sTw7SbmXKoBGnYVP73P+6P2yKCn3wdOuM1hEXLpomAYuSimllHIb2lWklFJKKbehgYtSSiml3IYGLkoppZRyGxq4uAljTC1jzDfGmL3GmN3GmDeMMSb3I1Vx0vetaOh9zh+9XxYl5T5czddpjFlvjNlhjNlpjIk1xlTI7RgNXNxHKvCCiAQDrYAbgX7FWyWVB/q+FQ29z/mj98uipNyHq/k6+4hICxFpjmU1+udyO6BEBy7GmBrGmH+MMZKx6q6ry29ojJmVEUmmGWPW5ZD32oyIOMkYc9wYM8EY42XdLyInRGRbxr8vAzuBWq6u89XKGDM4433Kuj1eCOcq0e+bMaa0MeZFY8xBY8wlY8yfxpiphXAet7/Pxph12fxeijGmvYvPdVXfL2NMmDHmp4z/U48ZYz42xrh81dTivg/GmLsyzn3JGHPEGBPuokvLep6r9v3Oa91yq1dGfc5l5C0F+JGHBxd7zJL/BfQm8A+Wm1UYmgG9gO8B7+wyGWP8gTXAHqAv0AB4C0tgOcZJ/srAXUA3l9f46ncbcDHT68OFcI6S/r7NwXKf/4vlyeu1gGsL4TyecJ+fBLI2bU/A8q12q4vPddXeL2NMHyAaeA/LN+bqwCTgS2NMGxFJL2jZThTbfTDGdMTyfLsPgVFYWi5eN8aki8i0Al1N9q7a9zsvdctPvYwxK4AbsDyqZ2SuZxeRErkBtwBnsfzyCVA+h7ytgJrZ7OsFeGWzr1Smf8cC67LJ9xKQAFTIlPY8kJQ5LSO9LLAWGFnc97CI36/Bub1P+r655D73AFKAa/U+F+j+lcn4f+X9knS/gBhge5a0Phl/s8Gech+Ar4GNWdLeynjPy3jKdebh/c61bvmpV8Y+L+AN4Pnczl8iu4oymqrewfLNKC8rFEYAq40xVbKU8yCwDMu3UweS928ZPYGvxf6BjTFAOeDWLPWeB/wsIm/lseySLAJ93/LrEeBbEdmTj2Mi0Pts1QPwx9L6kJ0IPO9+eQPnsqQlWk+VzTERuN99aAmszpK2Cst7nl3XYATud505ymPd8lSvTGWmAXOBB3MruEQGLsDjWCLP9/KYfzCWb6FfmYwRz8aYvliaC18Wkay/yPnVFEuTvI2I/IElMm2aKXkW8Dd5aUrzXL8ZY1KNMfuNMcNyyTsYfd/y60bggDHmXWPM+Yy+6cW5jFUYjN5nqzDgT2BjDnkG43n360PgZmPMg8aYCsaYxli6inIKggfjfvfBB7icJc36OjibYwbjftfpCrnWyxjjb4z5T6Ys9wC/5lZwiQtcMvr3JgLhIpKSl2NEJAHojiWqXmqM6QXMB94SkdddUC1//v12kllCxj5r3+qjQAjwszHmF2PM0y44t7s4AYwFHgB6Y+lbnWmMGZHdAfq+FUgglv9oW2L5EH4YaAN8bozz6ZN6ny2MMb5YukcWSEbbtzOeeL9E5EssvzdRWFpe9mNp+r8nh2Pc8T4cwjIWI7O2GT8DnB3gptfpCrnWK+Pn8oxBvruA64Bnciu4JA7OfQX4XkRW5OcgETlhjOkKfAd8CfxPRF4ojApmc/7vyL7J1eOJyNdY+petVhpjfIAxxpjp2TVd6vuWbyZj6ysifwEYY04A67E0aX/j7CC9z4AloPYj524iwPPulzGmMzATmA6sBP6DpYvkc2NM14xuAGf1cLf7MBPLF6YhWMZ2tAWss4qy7T5xw+ssEiJyGMdAMFclqsXFGNMMSx/+BGNMJWNMJcA3Y3dFY0y5XIooh2XwXTpQPmP6liskABWdpPtn7FPOxWL5llM3l3z6vuVdArDLGrRk2ISlOTy3mUUl/T6HAYckY9ppHnjS/XoLWCoiL4jIOhGZj2XmSicsM0py4k734UPg/YztLJYZRhMz9p3M5Vh3uk5XKLR6lajABWiEZRDZFiw3LoF/x7n8iWXArlPGmLpYBmF9j2Vg0Z3kfYxMbvZh3xeJMaYWlqBqn9MjFPw73z/bZnl93/JtL86/mRly+EZZ0u+zMaYilsGIuba2ZOSvi2fdr6bAL5kTRGQ/lqULGmR3kLvdBxFJE5HhQFWgOZaWpe8zdn+f3XHudp0uUmj1KmldRZuAzlnSegAvYJmW5nRNkIzBQ6uB34D+InIxY3DVCmNMgoi8fIX1Wgk8Z4y5RkT+zkgbiOWPfv0Vlu3JQrHMCvvd2U593wpkOfBfY0wVEbHOuLsFS8C/w9kBep8BuBvLgP9cAxcPvV+/A60zJxhjgrG0MsQ5O8Cd70PGuJUEAGPMk8BmEXH6YezO13mFCq9euc2X9vSNPKwPguU/8204zonvi2W0ePdsjvPF8uEaiqWVZ3em176Z8vljGXy6GugKDMWyMN6k4r4/V8sGLMISYPbE8m3lk4z37f/0fXPpfa6AZdntLVjGbNwLHAVW633O8b59BfySx7wed7+wDKhMx9Jl1BW4D8sA3SOAn6fcB6AdlrW/umJZMn8hcB5oXsLe71zrVpj1KpSLcqeNvAUujYEq2ezrAJhs9tXNKNvZVjdL3muBb7FEoyew9Js6XZioJG7A5Iz/CJMy7tF24IFcjtH3rWD3uiGwAriA5VvlHMBf73O2118l4wPoxTzm97j7haUr8QksS8lfAI5hmTlT35PuA5YZdluxfACfxzLQ9voS+H7nqW6FVS+TUbhSSiml1FWvpA3OVUoppZQb08BFKaWUUm5DAxellFJKuQ0NXJRSSinlNjRwUUoppZTb0MBFKaWUUm5DAxellFJKuQ0NXJRSSinlNjRwUUoppZTb0MBFKVUsjDG/5CFPROHXRCnlTnTJf6XUVcsYIyJiirseSqmrh7a4KKWKhTFGMv37ZmPMRmPMdmPMBmPMtcaYqRn7fjHGbCy+miqlriba4qKUKhbW1hRjTACwEughIgnGmI7AFBFpry0uSqmstMVFKVXc2gONgLUZ417eA2rkpwBjzDpjjGRsz2bZd5cxZpUx5i9jzGVjzDFjTKwxpkc+z7HMGLMrh/3vGmMSjTFljTERmeoTm5/zKKVypoGLUqq4lQI2i0jLTFvtApSzFksQFGNNyOhuWgQcAx4DugIvAuWAlcaYBvkoPxq4zhhzbdYdxhgvIBRYLCKXgP9l1OXnAlyHUioHGrgopYrbZqCNMaYVgDGmlDGmdca+v40xFfJYzlkR+V5ETmaU0xd4FnhURB4Wkc9FZIOIfCIidwB9gIv5qOcSIAkY5GRfZ+A/WIIbRORPEfkeOJ+P8pVSeaCBi1LK5YwxZ40xE40x8zK6aM4aY2YZY8pmzSsifwH9gfeMMTuA3cDdGbunAz8WcHDus8BWEZnjbKeILBOR41nqfbMxZr0xJimj3h8YY67JyH8BWAYMdFJcGHAa+LYA9VRK5YMGLkoplzLG1Ab8gRHAWSxByTRgKPBURp6qGfsAEJFNItJBRFqISLCIjM1IHysiTUXk5nzWoTSWrppV+TimI7AGOIml2+dZoBfwUaZs0UAjY0ybTMd5A/2ABSKSlp96KqXyr3RxV0Ap5XFaZPx8V0RezPj3t8aYW4AexpgvgRXA5EKsQ2WgLHA0c6IxxgBemZLS5N+pla9hGWszMFP+Y8A3xpjrRORXLLOfErG0sGzPyNYdS6AWXQjXoZTKQltclFKu1gLLWJA3s6QfAAJEZL+INBCRt4qgLlnXexgJpGTarC1AvlhaaBYYY0pbN2BTRr42ACJyGVgMDMgIgsDSdfQ7sKWQr0UphQYuSinXawFsyBi7klkgcNxJ/sLwF3AJqJkl/RPghowtM38sLTEzsA9sLgHeQK1MeaOB2kB7Y4wP0BeIydRyo5QqRNpVpJRytRaA3WDajAGutwFTMz2jKAAoD/yR8foFEfnaFRUQkVRjzBagGzAuU/op4FRGnTIfkoildSYCSzdWVpkDrrUZZYQB1YFr0G4ipYqMBi5KKZcxxvgBDXBsWXk64+d0EflvRt7BQCcRGVxI1ZkGfGGMeUBEPskpo4hcMMZ8DzQRkQm55E0zxizAMui4BrBXRHa4qtJKqZxp4KKUcqXrsXRBNzDGjMPS8tINy9iS+0UksagqIiJLjDHTgDnGmM5YpjKfwTJwt1tGtn8yHfI8loG46UAs8DeWLqE7gNEiciBT3mjg/7BM2x5fmNehlLKngYtSypVaYFl07Q4s04hHA/uB/iKypKgrIyIjjDEbgCeB2Vi6deKxDKTtJSIrM+XdlDHz6b9YxsJ4YRl0+xUZ3UuZ8m4xxsQBddFuIqWKlAYuSilXagHsEpGdZMzEKUImYyZQ5inOiMjnwOd5KUBEfgDy9AwjEamXQ0VKYWl50gdEKuViOqtIKeVKLYDiGu/RD8tMoGeK6fyZjcNSl1uKuyJKeRptcVFKuUTGuibXAx8Xw+mHYekGgn9nKRWnKGB5xr/P5pRRKZU/GrgopVwio3smrw9EdPW59xfHebOT8QykolqzRqkSxeiaSUoppZRyFzrGRSmllFJuQwMXpZRSSrkNDVyUUkop5TY0cFFKKaWU29DARSmllFJuQwMXpZRSSrkNDVyUUkop5TY0cFFKKaWU2/h/hi5cw0m045gAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Prepare arrays of Jet-pT\n",
    "jet_pt_pp_in = jets_pp[:,1]\n",
    "jet_pt_pbpb_in = jets_pbpb[:,1]\n",
    "jet_pt_pbpb1_in = jets_pbpb1[:,1]\n",
    "jet_pt_pbpb2_in = jets_pbpb2[:,1]\n",
    "\n",
    "bins_13 = [400,500,1000]\n",
    "y_axis13= [0.704,0.742]\n",
    "yerror_13 =[0.086,0.111]\n",
    "\n",
    "pt_min = 400\n",
    "pt_max = 1000\n",
    "\n",
    "n_ev_pp = 352683\n",
    "n_ev_pbpb = 35000\n",
    "n_ev_pbpb1 = 12790\n",
    "n_ev_pbpb2 = 6461\n",
    "sigma = 5.15\n",
    "d_eta = 2*2\n",
    "\n",
    "# Fill Histogram\n",
    "n_pp, pt = np.histogram(jet_pt_pp_in, bins= bins_13 )\n",
    "n_pbpb, pt = np.histogram(jet_pt_pbpb_in, bins= bins_13 )\n",
    "n_pbpb1, pt = np.histogram(jet_pt_pbpb1_in, bins= bins_13 )\n",
    "n_pbpb2, pt = np.histogram(jet_pt_pbpb2_in, bins= bins_13 )\n",
    "\n",
    "# Statistical Errors\n",
    "err_n_pp = np.sqrt(n_pp)\n",
    "err_n_pbpb = np.sqrt(n_pbpb)\n",
    "err_n_pbpb1 = np.sqrt(n_pbpb1)\n",
    "err_n_pbpb2 = np.sqrt(n_pbpb2)\n",
    "\n",
    "# bin width\n",
    "dpt = (pt[1:]-pt[:-1])\n",
    "# bin center\n",
    "pt = pt[0:-1] + 0.5*dpt\n",
    "\n",
    "# Jet Spectrum\n",
    "dn_dpt_pp = (n_pp*sigma)/n_ev_pp/dpt/d_eta\n",
    "dn_dpt_pbpb = (n_pbpb*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "dn_dpt_pbpb1 = (n_pbpb1*sigma)/n_ev_pbpb1/dpt/d_eta\n",
    "dn_dpt_pbpb2 = (n_pbpb2*sigma)/n_ev_pbpb2/dpt/d_eta\n",
    "\n",
    "# Errors \n",
    "err_dn_dpt_pp = (err_n_pp*sigma)/n_ev_pp/dpt/d_eta\n",
    "err_dn_dpt_pbpb = (err_n_pbpb*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "err_dn_dpt_pbpb1 = (err_n_pbpb1*sigma)/n_ev_pbpb1/dpt/d_eta\n",
    "err_dn_dpt_pbpb2 = (err_n_pbpb2*sigma)/n_ev_pbpb2/dpt/d_eta\n",
    "\n",
    "# Generate Plots\n",
    "fig = plt.figure()\n",
    "\n",
    "# Calculate RAA and error\n",
    "raa =  dn_dpt_pbpb/dn_dpt_pp\n",
    "error_raa = ratio_error(dn_dpt_pbpb,err_dn_dpt_pbpb,dn_dpt_pp,err_dn_dpt_pp)\n",
    "\n",
    "raa1 =  dn_dpt_pbpb1/dn_dpt_pp\n",
    "error_raa1 = ratio_error(dn_dpt_pbpb1,err_dn_dpt_pbpb1,dn_dpt_pp,err_dn_dpt_pp)\n",
    "\n",
    "raa2 =  dn_dpt_pbpb2/dn_dpt_pp\n",
    "error_raa2 = ratio_error(dn_dpt_pbpb2,err_dn_dpt_pbpb2,dn_dpt_pp,err_dn_dpt_pp)\n",
    "\n",
    "\n",
    "plt.errorbar(pt, y_axis13,ms=8,mfc='blue',mec='black',mew=2,fmt= 's',label=\"CMS [JHEP 05(2021)284]\",\n",
    "            xerr=None, yerr=yerror_13, color='black')\n",
    "\n",
    "plt.errorbar(pt, raa, ms=9,mfc='red',mec='black',fmt='o', label=r\"MATTER+LBT\",\n",
    "             xerr=0.5*dpt, yerr=error_raa, color='red')\n",
    "\n",
    "plt.errorbar(pt, raa1, ms=9,mfc='yellow',mec='black',fmt='o', label=r\"MATTER+MARTINI\",\n",
    "             xerr=0.5*dpt, yerr=error_raa1, color='yellow')\n",
    "\n",
    "plt.errorbar(pt, raa2, ms=9,mfc='cyan',mec='black',fmt='o', label=r\"MATTER+AdS/CFT\",\n",
    "             xerr=0.5*dpt, yerr=error_raa2, color='cyan')\n",
    "\n",
    "\n",
    "\n",
    "#axes setting\n",
    "plt.legend(loc='upper left')\n",
    "plt.xscale('log')\n",
    "plt.xlabel(r\"$p^{\\mathrm{jet}}_{\\mathrm{T}}$ [GeV]\")\n",
    "plt.ylabel(r\"$R^{\\mathrm{jet}}_{\\mathrm{AA}}$\")\n",
    "plt.xlim(pt_min,pt_max)\n",
    "plt.ylim(-0.3,1.6)\n",
    "plt.axhline(1, color = \"blue\", linestyle=\"dashed\", linewidth=1.5) \n",
    "plt.text(pt_min+5,-0.05, r\"PbPb(0-10%)\"r\"$\\sqrt{s_{\\mathrm{NN}}}=5.02 TeV$\")\n",
    "plt.text(pt_min+5,-0.15,r'$anti-k_{\\mathrm{t}},R=0.4,(|y_{\\mathrm{jet}}|<2)$')\n",
    "plt.text(pt_min+5,-0.25,r\"$\\hat{q}=\\hat{q}^{\\mathrm{run}}_{\\mathrm{HTL}}f(Q^{\\mathrm{2}}),\\alpha^{\\mathrm{fix}}_{\\mathrm{s}}=0.35,Q_{\\mathrm{sw}}=2GeV$\")\n",
    "\n",
    "# save plot to the Desktop\n",
    "plt.tight_layout()\n",
    "plt.savefig('Raa(0-10)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55bbf083",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "009adc34a3a8452e02639bd458ec0b5a4261dfbc6c9492726e6a0ec46c9767d3"
  },
  "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.6.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}