cms_double_ads_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",
    "#|0.2|1->0.3|2->0.4|3->0.6|4->0.8|5->1.0|\n",
    "file_pp = '../../data/music_cms/cms_pp_vir0.5_R0.2_trkpT4.0_maxT250_pt200-1001_m1-7_merged_315769.dat'\n",
    "file_pp1 = '../../data/music_cms/cms_pp_vir0.5_R0.3_trkpT4.0_maxT250_pt200-1001_m1-7_merged_338726.dat'\n",
    "file_pp2 = '../../data/music_cms/cms_pp_vir0.5_R0.4_trkpT4.0_maxT250_pt200-1001_m1-7_merged_352683.dat'\n",
    "file_pp3 = '../../data/music_cms/cms_pp_vir0.5_R0.6_trkpT4.0_maxT250_pt200-1001_m1-7_merged_369925.dat'\n",
    "file_pp4 = '../../data/music_cms/cms_pp_vir0.5_R0.8_trkpT4.0_maxT250_pt200-1001_m1-7_merged_380962.dat'\n",
    "file_pp5 = '../../data/music_cms/cms_pp_vir0.5_R1.0_trkpT4.0_maxT250_pt200-1001_m1-7_merged_389140.dat'\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d70a0423",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please set the path for pbpb jet data file(AdS/CFT ONLY)######\n",
    "#|0.2|1->0.3|2->0.4|3->0.6|4->0.8|5->1.0|\n",
    "file_pbpb = '../../data/music_cms/ADS/R0.2/cms_ads_0-10_R0.2_s1-4_6461.dat'\n",
    "file_pbpb1 ='../../data/music_cms/ADS/R0.3/cms_ads_0-10_R0.3_s1-4_6461.dat'\n",
    "file_pbpb2 ='../../data/music_cms/ADS/R0.4/cms_ads_0-10_R0.4_s1-4_6461.dat'\n",
    "file_pbpb3 ='../../data/music_cms/ADS/R0.6/cms_ads_0-10_R0.6_s1-4_6461.dat'\n",
    "file_pbpb4 ='../../data/music_cms/ADS/R0.8/cms_ads_0-10_R0.8_s1-4_6461.dat'\n",
    "file_pbpb5 ='../../data/music_cms/ADS/R1.0/cms_ads_0-10_R1.0_s1-4_6461.dat'\n"
   ]
  },
  {
   "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'] = [7.5, 6.5]\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",
    "\n",
    "data_pp1 = np.loadtxt(file_pp1, delimiter=',')\n",
    "data_pbpb1 = np.loadtxt(file_pbpb1, delimiter=',')\n",
    "\n",
    "data_pp2 = np.loadtxt(file_pp2, delimiter=',')\n",
    "data_pbpb2 = np.loadtxt(file_pbpb2, delimiter=',')\n",
    "\n",
    "data_pp3 = np.loadtxt(file_pp3, delimiter=',')\n",
    "data_pbpb3 = np.loadtxt(file_pbpb3, delimiter=',')\n",
    "\n",
    "data_pp4 = np.loadtxt(file_pp4, delimiter=',')\n",
    "data_pbpb4 = np.loadtxt(file_pbpb4, delimiter=',')\n",
    "\n",
    "data_pp5 = np.loadtxt(file_pp5, delimiter=',')\n",
    "data_pbpb5 = np.loadtxt(file_pbpb5, delimiter=',')\n",
    "\n",
    "# Indices of the data array\n",
    "i_pp = data_pp[:,0] \n",
    "i_pbpb = data_pbpb[:,0] \n",
    "\n",
    "i_pp1 = data_pp1[:,0] \n",
    "i_pbpb1 = data_pbpb1[:,0]\n",
    "\n",
    "i_pp2 = data_pp2[:,0] \n",
    "i_pbpb2 = data_pbpb2[:,0]\n",
    "\n",
    "i_pp3 = data_pp3[:,0] \n",
    "i_pbpb3 = data_pbpb3[:,0]\n",
    "\n",
    "i_pp4 = data_pp4[:,0] \n",
    "i_pbpb4 = data_pbpb4[:,0]\n",
    "\n",
    "i_pp5 = data_pp5[:,0] \n",
    "i_pbpb5 = data_pbpb5[:,0]\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",
    "\n",
    "jet_id_pp1 = np.where(i_pp1 < 0.1)\n",
    "jet_id_pbpb1 = np.where(i_pbpb1 < 0.1)\n",
    "\n",
    "jet_id_pp2 = np.where(i_pp2 < 0.1)\n",
    "jet_id_pbpb2 = np.where(i_pbpb2 < 0.1)\n",
    "\n",
    "jet_id_pp3 = np.where(i_pp3 < 0.1)\n",
    "jet_id_pbpb3 = np.where(i_pbpb3 < 0.1)\n",
    "\n",
    "jet_id_pp4 = np.where(i_pp4 < 0.1)\n",
    "jet_id_pbpb4 = np.where(i_pbpb4 < 0.1)\n",
    "\n",
    "jet_id_pp5 = np.where(i_pp5 < 0.1)\n",
    "jet_id_pbpb5 = np.where(i_pbpb5 < 0.1)\n",
    "\n",
    "# Extract jets\n",
    "jets_pp = data_pp[jet_id_pp]\n",
    "jets_pbpb = data_pbpb[jet_id_pbpb]\n",
    "\n",
    "jets_pp1 = data_pp1[jet_id_pp1]\n",
    "jets_pbpb1 = data_pbpb1[jet_id_pbpb1]\n",
    "\n",
    "jets_pp2 = data_pp2[jet_id_pp2]\n",
    "jets_pbpb2 = data_pbpb2[jet_id_pbpb2]\n",
    "\n",
    "jets_pp3 = data_pp3[jet_id_pp3]\n",
    "jets_pbpb3 = data_pbpb3[jet_id_pbpb3]\n",
    "\n",
    "jets_pp4 = data_pp4[jet_id_pp4]\n",
    "jets_pbpb4 = data_pbpb4[jet_id_pbpb4]\n",
    "\n",
    "jets_pp5 = data_pp5[jet_id_pp5]\n",
    "jets_pbpb5 = data_pbpb5[jet_id_pbpb5]\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",
    "assoc_pp1= []\n",
    "for i in range(len(jet_id_pp1[0])-1):\n",
    "  chunck = data_pp1[jet_id_pp1[0][i]+1:jet_id_pp1[0][i+1]]\n",
    "  assoc_pp1.append(chunck)\n",
    "chunck = data_pp1[jet_id_pp1[0][-1]+1:]\n",
    "assoc_pp1.append(chunck)\n",
    "\n",
    "assoc_pp2= []\n",
    "for i in range(len(jet_id_pp2[0])-1):\n",
    "  chunck = data_pp2[jet_id_pp2[0][i]+1:jet_id_pp2[0][i+1]]\n",
    "  assoc_pp2.append(chunck)\n",
    "chunck = data_pp2[jet_id_pp2[0][-1]+1:]\n",
    "assoc_pp2.append(chunck)\n",
    "\n",
    "assoc_pp3= []\n",
    "for i in range(len(jet_id_pp3[0])-1):\n",
    "  chunck = data_pp3[jet_id_pp3[0][i]+1:jet_id_pp3[0][i+1]]\n",
    "  assoc_pp3.append(chunck)\n",
    "chunck = data_pp3[jet_id_pp3[0][-1]+1:]\n",
    "assoc_pp3.append(chunck)\n",
    "\n",
    "assoc_pp4= []\n",
    "for i in range(len(jet_id_pp4[0])-1):\n",
    "  chunck = data_pp4[jet_id_pp4[0][i]+1:jet_id_pp4[0][i+1]]\n",
    "  assoc_pp4.append(chunck)\n",
    "chunck = data_pp4[jet_id_pp4[0][-1]+1:]\n",
    "assoc_pp4.append(chunck)\n",
    "\n",
    "assoc_pp5= []\n",
    "for i in range(len(jet_id_pp5[0])-1):\n",
    "  chunck = data_pp5[jet_id_pp5[0][i]+1:jet_id_pp5[0][i+1]]\n",
    "  assoc_pp5.append(chunck)\n",
    "chunck = data_pp5[jet_id_pp5[0][-1]+1:]\n",
    "assoc_pp5.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)\n",
    "\n",
    "assoc_pbpb3= []\n",
    "for i in range(len(jet_id_pbpb3[0])-1):\n",
    "  chunck = data_pbpb3[jet_id_pbpb3[0][i]+1:jet_id_pbpb3[0][i+1]]\n",
    "  assoc_pbpb3.append(chunck)\n",
    "chunck = data_pbpb3[jet_id_pbpb3[0][-1]+1:]\n",
    "assoc_pbpb3.append(chunck)\n",
    "\n",
    "assoc_pbpb4= []\n",
    "for i in range(len(jet_id_pbpb4[0])-1):\n",
    "  chunck = data_pbpb4[jet_id_pbpb4[0][i]+1:jet_id_pbpb4[0][i+1]]\n",
    "  assoc_pbpb4.append(chunck)\n",
    "chunck = data_pbpb4[jet_id_pbpb4[0][-1]+1:]\n",
    "assoc_pbpb4.append(chunck)\n",
    "\n",
    "assoc_pbpb5= []\n",
    "for i in range(len(jet_id_pbpb5[0])-1):\n",
    "  chunck = data_pbpb5[jet_id_pbpb5[0][i]+1:jet_id_pbpb5[0][i+1]]\n",
    "  assoc_pbpb5.append(chunck)\n",
    "chunck = data_pbpb5[jet_id_pbpb5[0][-1]+1:]\n",
    "assoc_pbpb5.append(chunck)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a38aafd6",
   "metadata": {},
   "source": [
    "## Jet-$R_{\\mathrm{AA}}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "760626b2",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R=0.3/0.2= [1.08065017 1.04772615 1.09741014 1.08342547 1.0894024 ]\n",
      "R=0.4/0.2= [1.11901673 1.0882585  1.16052045 1.18320653 1.12058585]\n",
      "R=0.6/0.2= [1.17554052 1.15685342 1.18915696 1.21801265 1.25177887]\n",
      "R=0.8/0.2= [1.21069166 1.15772495 1.21599646 1.28016079 1.43576899]\n",
      "R=1.0/0.2= [1.21506895 1.17263609 1.22558194 1.26044169 1.48563593]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAHBCAYAAADjMSi8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABL/0lEQVR4nO3de3yU9Zn38c8FRCFRYAA1ja0HhBaLrXSZrYfaYgVRVNS2oVhr1dqS7nZta2J62EpCIK5P0Rjsbm3rYLdaa1sLtfqgYjF4aFX0AbbVroiKxYoESoQBlBCFcD1/3JMYksl5Tpn5vl+v+5XMffjd18ydzFzzu38Hc3dERERE2huU7gBEREQkMylJEBERkbiUJIiIiEhcShJEREQkLiUJIiIiEpeSBBEREYkr45IEM5tlZv/XzDab2dtmttbMvtCD4w41s5vNbJuZ7TGzB83suBSELCIikpUyLkkAyoC3gVLgQuAx4Fdm9o1ujvtP4EqgHCgGxgCPmNnQ5IUqIiKSvSzTBlMyszHu/ma7db8CTnP34zs55v3Aa8BV7v6L2LqjgY3A19399uRGLSIikn0yriahfYIQ82egqIvDpsd+3tumnM3Ak8CMxEUnIiKSOzIuSejEacDLXWyfALzh7m+3W/9ibJuIiIj00pB0B9AdM5sKXAxc1cVuIWBnnPXR2LYODjnkEN+3b1/r4/e9730UFXVVWSGZoKGhgSOOOCLdYQxIufjaZctzHijPI9PiTHc8qTp/fX09W7ZsaXm4z90PSVTZGZ0kxHon/Aq4393vSGTZeXl5vPvuu4ksUlIgHA6zZs2adIcxIOXia5ctz3mgPI9MizPd8aTj/GaWl8jyMvZ2g5mNApYDfwe+2M3uUWBEnPWh2LYOxowZ06/4JD1KSkrSHcKAlYuvXbY854HyPDItznTHk+7zJ0LG9W4AMLN8oA44iqBXw7Zu9r8KiAAj3H1Pm/UrgZ3u/rn2x4TDYc+kjFdERKS/zAx3t0SVl3E1CWY2BFgCjAfO7S5BiFkR+/mZNuUUAZ8kqI0QERHJBVu636XnMrFNwo+B84BvAaPNbHSbbX9293diNQS4+9TYzzfM7GfALWZmQANQRXCr4pepDF5ERCSN6hNZWCYmCS1jHvwwzrbjCQZNGhxn2zeBPUAtkA88AXzB3ZuSEKOIiEjWy7gkwd2P68E+Z8ZZ9w7BkM5liY9KREQk92RckpApdu/ezbZt22g7loKkXl5eHkceeSTDhw9PdygiIjlHSUIcu3fv5h//+AdHH300w4YNI2jmIKnm7uzdu5fNmzcDKFEQEUmxjOvdkAm2bdvG0UcfTX5+vhKENDIz8vPzOfroo9m2rSedXEREJJGUJMSxb98+hg0blu4wJGbYsGG67SOSQyy2SPopSeiEahAyh66FiEh6KElIBbNgERERGUCUJIiIiEhcShKSLBqNshioBhYvXkw0Gne+qaSoqqrCzFqXwsJCLrjgAp5//vl+l71u3TqmTp1Kfn4+RUVFVFZW0tzc3OPjDxw4QDgcxsx44IEH+h2PiIgkXs4mCQ0NDYTDYSKRSFLKd3eqKyoYW1TESuAdoK60lLFFRVRXVJCqibVGjBjBqlWrWLVqFbfccgsvv/wyZ599Njt27OhzmdFolGnTpmFm3H///VRWVnLzzTczb968Hpdx++2388Ybb/Q5BhEReU8kEiEcDgMkdopjd8/JZfLkyd6ZdevWdbqtpxbMnev/lJ/vm8G9zbIZ/J/y833B3Ln9Pkd35s2b56NHjz5o3apVqxzwu+++u8/l3nDDDT5y5EjftWtX67qFCxf6sGHDDlrXmR07dviYMWP89ttvd8CXLVvW7TGJuCYiMjC0vFFL7wFrPIGflTlbk5BM0WiU2poaljU2UtRuWxGwrLGR2poadu7cmfLYTj75ZAA2bdrU5zKWL1/OOeecc9DgRpdccgl79+7liSee6Pb4iooKPvGJTzB16tQ+xyAiIsmnJCGRYr0Ylo4axfSmpg4JQosiYHpTE0tCoZT3fHj99dcBOP7441vX7d+/v9vF29weWb9+PRMmTDio3GOOOYb8/HzWr1/f5fmff/55/vu//5uampoEPisREUkGJQlJsBUY380+42L7pULLB/2rr77K1VdfzaRJk7joootat+fl5XW73Hnnna37R6NRRo4c2eE8oVCo24aZ3/jGN7j66qsZN25cwp6fiIgkh+ZuSKTYt+3CxYupKy2FPXs63XVDQQHTFi2COXOSGtL27dvJy8trfTx69GhWr17NoYce2rpu9erV3ZbTtuahr37zm9/w0ksvsWzZsn6XJSIiyackIQmKi4v5zje/ST3EveWwGVjR3Mxts2YlPZYRI0ZQV1dHc3Mzzz33HOXl5Vx66aU89dRTDBoUVCRNmjSp23IGDx7c+nsoFGLXrl0d9olGo4RCobjH79u3j29/+9t897vf5cCBA+zcuZPdu3cDsGfPHt566y0OP/zwPjxDERFJFt1uSIJQKERZeTkz8/Opb7etHrgwP5+y8vK4VfaJNmTIEMLhMKeccgolJSXceuutPPPMMyxZsqR1n97ebpgwYUKHtgebNm2isbGxQ1uFFnv27OGNN96grKyMUChEKBRqbUR5ySWX8LGPfSwJz15ERPpDNQlJMnfBAgAm1tQwvamJcQS3GFY0N1NWVta6PdUuu+wyFi5cyMKFC5k9ezbQ+9sNM2bM4Kabbjro2/8999zDsGHDmDJlStzjDzvsMB577LGD1m3dupUvfOEL3HDDDZx11ll9fUoiIpIsiexPOZCWZI+T0GLHjh0eAV8AHolEPBqNJqzs7sQbJ8Hd/e6773bA6+rq+lTujh07vLCw0KdNm+aPPPKI33bbbV5QUODXXXfdQfudcMIJftVVV3VazsaNGzVOgoh0oHES+g6NkzCwhEIh5gAVwJw5c1Jyi6E7s2fPZvz48dx44419Oj4UCrFy5Uqam5uZOXMm8+bNo7S0lPnz5x+03/79+3s1VLOIiGQWc0/N8MCZJhwO+5o1a+Jue/HFFznxxBMTd7KWcRBy9LVOhIRfExHJWC0jx+gds/fMbK27hxNVntokpIKSAxERGYB0u0FERETiUpIgIiIicSlJEBERkbhyNkloaGggHA4TiUTSHYqIiEi/RCIRwuEwwJhElpuzDRePOOIIOuvdICIiMpCUlJRQUlKCmb2ZyHJztiZBREREuqYkIQXM3hsqQUREZKBQkiAiIiJxKUnIYlVVVZhZ61JYWMgFF1zA888/3++y161bx9SpU8nPz6eoqIjKysoeD8F877338s///M8MGzaM0aNHc+6557Jnz55+xyQiIomlJCHJotEosBioZvHixbHHqTNixAhWrVrFqlWruOWWW3j55Zc5++yz2bFjR5/LjEajTJs2DTPj/vvvp7Kykptvvpl58+Z1e+ztt9/OpZdeyowZM1i+fDm3334748ePZ//+/X2OR0REkiSRs0UNpCXZs0AeOHDA585d4EOHjnSY7XCdFxR83ocOHelz5y7wAwcO9Psc3Yk3C+SqVasc8LvvvrvP5d5www0+cuRI37VrV+u6hQsX+rBhww5a115DQ4MfdthhHolEen1OzQIpkjs0C2TfoVkgB4bKyuuprb2PpqYXgN8A17Nnzz00Nb1Abe19VFZen5a4Tj75ZAA2bdrU5zKWL1/OOeecw/Dhw1vXXXLJJezdu5cnnnii0+N++9vfAnDFFVf0+dwiIpI6ShKSIBqNUlNTS2PjMqCo3dYiGhuXUVNTy86dO1Me2+uvvw7A8ccf37pu//793S7eZpKq9evXM2HChIPKPeaYY8jPz2f9+vWdnvvZZ5/lQx/6ED/72c94//vfT15eHqeccgpPP/10gp+liIgkgpKEBGrp6jhq1FKamqbTMUFoUURT03RCoSUp6R7Z8kH/6quvcvXVVzNp0iQuuuii1u15eXndLnfeeWfr/tFolJEjR3Y4TygU6rLNxdatW3nppZe4/vrrWbhwIcuWLaOgoIBzzz2Xf/zjHwl9ziIi0n85O+Jicm0Fxnezz7jYfsm1fft28vLyWh+PHj2a1atXc+ihh7auW716dbfltK156Ct35+2332bJkiWce+65AJx++ukce+yx/OhHP6K6urrf5xARkcRRkpBALTXyixcXUlpaR1e9+goKNrBo0TTmzEluTCNGjKCuro7m5maee+45ysvLufTSS3nqqacYNCioSJo0aVK35QwePLj191AoxK5duzrsE41GCYVCnZYRCoUwM84888zWdcOHD2fy5MmsW7eu509KRERSQrcbkqC4uJjm5hVAfSd7bKa5eQWzZs1KeixDhgwhHA5zyimnUFJSwq233sozzzzDkiVLWvfp7e2GCRMmdGh7sGnTJhobGzu0VWjrxBNPbG0x25a7tyYsIiKSOfTOnAShUIjy8jLy82fSMVGoJz//QsrLy+Le10+2yy67jIkTJ7Jw4cLWdatXr+52mTlzZuv+M2bM4A9/+ANvvfVW67p77rmHYcOGMWXKlE7PfcEFFwDw2GOPta7btWsXa9eube11ISIiGSSR/SkH0pLacRI+7/D9jBgnwd397rvvdsDr6ur6VO6OHTu8sLDQp02b5o888ojfdtttXlBQ4Nddd91B+51wwgl+1VVXHbTuoosu8sLCQr/jjjv8gQce8E996lM+ZswY37FjR5fn1DgJIrlD4yT0HRonITEaGhoIh8NEIpGklG9mVFdXUF//N2AaMJRFi6axZctGqqsrsDTO+DR79mzGjx/PjTfe2KfjQ6EQK1eupLm5mZkzZzJv3jxKS0uZP3/+Qfvt37+/w1DNv/zlL7n44ospKyujuLiYvLw8Hn300S7bMoiISNcikQjhcBhgTCLLNW93fzhXhMNhX7NmTdxtL774IieeeGLCztWSD+ToS50Qib4mIpK5Wr5C6S2z98xsrbuHE1WeejekgJIDEREZiHL2doOIiIh0TUmCiIiIxKUkQUREROJSkiAiIiJxKUkQERGRuJQkiIiISFxKEkRERCSujEsSzGycmd1mZs+bWbOZPd7D48JmtsLMdsSWOjM7Jcnh9ojNN2x++kZYFBER6YuMSxKAicB5wEvAyz05wMw+ANQRDA71pdgyBHjEzI5NUpwiIiJZLROThGXu/gF3nwW80MNjzgcOBz7j7g+6+4PAZ4ACgoQjbaLRKKwFnoDFixcHj1OkqqoKM2tdCgsLueCCC3j++ef7Xfa6deuYOnUq+fn5FBUVUVlZ2WGehnjWrFnD9OnTGTVqFKNGjWLatGk8++yz/Y5HREQSL+OSBHc/0IfD8oD9wJ42696OrUtLPb+7U1FVQdExRfC3IJLSH5dSdEwRFVUVpGrOjBEjRrBq1SpWrVrFLbfcwssvv8zZZ5/Njh07+lxmNBpl2rRpmBn3338/lZWV3HzzzcybN6/L4zZt2sS0adPYv38/d911F3fddRf79+/n7LPP5u9//3uf4xERkeTIlrkbfgcsAG42s/+IrasEosCSdARUOb+S2p/X0lTSBMODdXvYA7uh9ue1AFRXVSc9jiFDhnDqqacCcOqpp3Lcccdx2mmn8fDDD3PppZf2qcyf/vSn7N27l3vvvZfhw4dz9tlns3v3bqqqqvjOd77D8OHD4x734IMP8tZbb/H73/+eESNGAHD66aczZswYHnroIf71X/+1b09SRESSIuNqEvrC3euBTwOfA/4RWz4LnOPuDfGOqa+vP6gqvqqqKmHxRKNRam6uobG4sTVBaDUcGosbqbm5hp07dybsnD118sknA8G3+r5avnw555xzzkHJwCWXXMLevXt54oknOj1u3759DBkyhIKCgtZ1hx12GEOGDElZzYqISLZpe2sZKEpk2VmRJJjZ+whqDNYCM2LLWuBBMzsm3jFFRUW4e+uSiCShpRfDqCtG0XRsU8cEocVwaDq2idDloZT3fHj99dcBOP7441vX7d+/v9ul7Yf4+vXrmTBhwkHlHnPMMeTn57N+/fpOz/25z32O/Px8rr32WrZt28a2bdsoLS0lFAoxa9asBD9TEZHcUFVV1fpZBtQnsuysSBKAbxO0Syh294fd/WGCWoVmoDzl0bwNjOpmn1Gx/VKg5YP+1Vdf5eqrr2bSpElcdNFFrdvz8vK6Xe68887W/aPRKCNHjuxwnlAo1GXDzKKiIh577DF+97vfcdRRR3HUUUdx77338oc//IEjjjgioc9ZRET6L1vaJEwAXnD3fS0r3P1dM3sBOCFVQfi84Nv24sWLKf1xadAGoRMFbxWw6OuLmDNnTlJj2r59O3l5ea2PR48ezerVqzn00ENb161evbrbctrWPPTVli1bmDVrFpMnT+b2228H4NZbb+X888/n6aef5phj4lb6iIhImmRLkvB34DwzO8Td3wUws0OBk4BlqQ6muLiYb5Z9E3YT/5bDbmje0JySKvYRI0ZQV1dHc3Mzzz33HOXl5Vx66aU89dRTDBoUVCRNmjSp23IGDx7c+nsoFGLXrl0d9olGo4RCoU7LuOmmm9i3bx9Lly5tTVzOOussxo8fT01NDf/5n//Zy2cnIiLJlHG3G8ws38yKzawYOBo4ouWxmeXH9tlgZj9rc9jtBI01fm9m55vZBcB9wPuASIqfAqFQiPJry8lfmh8kCm3thvyl+ZRfWx63yj7RhgwZQjgc5pRTTqGkpIRbb72VZ555hiVL3uv00dvbDRMmTOjQ9mDTpk00NjZ2aKvQ1vr165k4ceJBNRuHHHIIEydO5NVXX03gsxYRkUTIxJqEI+nYbbHl8fHAawRxt361dfe1ZnYuMA+4K7b6r8DZ7v5cUqPtxIJ5CwCoubkmaMQ4KrjF0LyhmbJry1q3p9pll13GwoULWbhwIbNnzwZ6f7thxowZ3HTTTbz11lscfvjhANxzzz0MGzaMKVOmdFrGsccey0MPPcS7777LIYccAsA777zD//7v/zJz5sz+PC0REUmGti38c2mZPHmyd2bdunWdbuutHTt2ODNxPo1HIhGPRqMJK7s78+bN89GjR3dYf/fddzvgdXV1fSp3x44dXlhY6NOmTfNHHnnEb7vtNi8oKPDrrrvuoP1OOOEEv+qqq1ofr1mzxocMGeLnnXeeP/DAA75s2TI/99xzfciQIf6Xv/yly3Mm8pqISGZreaOW3gPWeAI/KzPudkO2CYVCMBmYAnPmzEnJLYbuzJ49m/Hjx3PjjTf26fhQKMTKlStpbm5m5syZzJs3j9LSUubPn3/Qfvv37z9oqObJkyfz8MMP89Zbb/GlL32Jyy+/nMbGRh555JHW8RtERCRzmOfoIDbhcNjXrFkTd9uLL77IiSeemLBztYyD0NL7QXov0ddERDJXy8gxesfsPTNb6+7hRJWXiW0Sso6SAxERGYh0u0FERETiUpIgIiIicSlJEBERkbiUJIiIiEhcShJEREQkrpxNEhoaGgiHw0QiKR+1WUREJKEikQjhcBhgTCLLzdkukEcccQSdjZOQaOrzKyIiyVRSUkJJSQlm9mYiy83ZmgQRERHpmpKELFZVVYWZtS6FhYVccMEFPP/88/0qd8OGDXzta1/jox/9KIMHD+bMM8/s8bG7du3iy1/+MqFQiBEjRvDFL36R7du39yseERFJjpy93ZAq0WgUli6FrVtZXFhIcXFxMJ9DiowYMYKHH34YgNdee43KykrOPvtsXnzxRUaNGtWnMl944QUeeughTj31VPbt29erYz//+c/z8ssvc/vttzNo0CC++93vcvHFF/OnP/2pT7GIiEjyqCYhSdydiupqisaOhZUr4Z13KK2ro2jsWCqqq0nVnBlDhgzh1FNP5dRTT+WSSy7hF7/4Bdu2bWtNHPpi5syZbNq0iSVLljBx4sQeH7dq1SpWrFjBnXfeyec+9zk+85nP8Mtf/pInn3ySurq6PscjIiLJoZqEJKm8/npq77uPphdegKIiAPYA1NdTO3MmANUVFSmPq2W2xU2bNvW5jEGD+pZbLl++nKOOOopPfepTres+/vGPc/zxx7N8+XKmTZvW55hERCTxVJOQBNFolJraWhqXLWtNEFoVFdG4bBk1tbXs3Lkz5bG9/vrrABx//PGt6/bv39/tkoiaj/Xr1zNhwoQO60888UTWr1/f7/JFRCSxlCQkkMWWUUuX0jR9escEoUVREU3TpxNasqT1mGRq+aB/9dVXufrqq5k0aRIXXXRR6/a8vLxulzvvvLPfcUSjUUaOHNlhfSgUCtpuiIhIRtHthmTYuhXGj+96n3Hjgv2SbPv27eTl5bU+Hj16NKtXr+bQQw9tXbd69epuy2lb8yAiIrlBSUICtVTILy4spLSuLmiD0ImCDRtYNG0ac5Ic04gRI6irq6O5uZnnnnuO8vJyLr30Up566qnWtgWTJk3qtpzBgwf3O5ZQKERDQ0OH9dFoNKU9PkREpGd0uyEJiouLaV6xAurr4++weTPNK1Ywa9aspMcyZMgQwuEwp5xyCiUlJdx6660888wzLFmypHWfVN1umDBhQty2B521VRARkfRSTUIShEIhysvKqJ05s2Pjxfp68i+8kLKysrj355PtsssuY+HChSxcuJDZs2cDqbvdMGPGDKqrq3nyySc544wzAFizZg1/+9vfmDFjRr/LFxGRxFKSkCQL5s4FoGbixKAR47hxFGzYQPOKFZSVlbVuTzUz4/vf/z5f/OIXWblyJVOnTm2ZFKTHGhsbeeihhwDYvHkzu3fvZunSpQCcd9555OfnAzBu3DimTJnCz372MwBOO+00pk+fzuWXX05NTU3rYEpnnHGGuj+KiGQg3W5IEjOjuqKC+r/9DaZNg6FDWTRtGls2bqS6ogKzZPdp6Nzs2bMZP348N954Y5+O37ZtG7NmzWLWrFk888wzrFu3rvXxtm3bWvfbv38/zc3NBx17zz33MGXKFK666iouv/xyJk+ezO9///t+PR8REUkOS9XIf5kmHA57Z7NAvvjii5x44okJO5dmgey/RF8TEclces/sOzNb6+69qx7uQs7WJDQ0NBAOh4lEIkk/l6M/dhERSZ5IJNJy63hMIsvN2TYJRxxxBJ3VJIiIiAwkJSUllJSUYGZvJrLcnK1JEBERka4pSRAREZG4lCR0IlcbdGYiXQsRkfRQkhBHXl4ee/fuTXcYErN3796D5p8QEZHUUJIQx5FHHsnmzZtpbGzUt9g0cncaGxvZvHkzRx55ZLrDERHJOTnbu6Erw4cPB6C+vp59+/alOZrclpeXx1FHHdV6TUREJHWUJHRi+PDh+mASEZGcptsNIiIiEpeSBBEREYlLSYKIiIjEpSRBRERE4lKSICIiInEpSRAREZG4lCSIiIhIXDmbJDQ0NBAOh4lEIukORUREpF8ikQjhcBhgTCLLtVwddjgcDvuaNWvSHYaIiLRjsZ+5+enUP2a21t3DiSovZ2sSREREpGtKEkRERCQuJQkiIiISl5IEERERiUtJgoiIiMSlJEFERETiUpIgIiIicSlJEBERkbgyLkkws3FmdpuZPW9mzWb2eC+O/ayZrTazvWa23cweNrOCJIYrIiKStTIuSQAmAucBLwEv9/QgM/sq8CtgOTAD+CrwCjAkCTGKiIhkvUz8AF3m7vcDmNlSejAOtZmNARYB33D3xW02/T45IYqIiGS/jKtJcPcDfTjs87GfdyYyFhERkVyWcUlCH51CcHviK2b2hpntM7Nnzez0dAcmIiIyUGVLklAIfAiYC3wXmAnsAR42s6PiHVBfX4+ZtS5VVVUpC1ZERCRRqqqqWj/LgKJElp3RU0W3tElw9zO72W8FcDYww90fjq0bDvwd+JG7V7Q/RlNFi4hkJk0V3XeaKjq+KMHf0+MtK9x9N7AW+HCaYhIRERnQsiVJeJEg+bR26w3oS0NIERGRnJctScIDsZ+fbllhZiOAycBzaYlIRERkgMu4cRLMLJ9gMCWAo4HhZlYce/yQuzea2QbgCXf/CoC7rzGz+4Gfmdn3gDeB7wD7gFtT+wxERESyQ8YlCcCRwJJ261oeHw+8RhD34Hb7XAbcBNQC+cBTwFnuHk1apCIiIlks45IEd3+Njm0L2u9zXJx1bwP/GltERESkn7KlTYKIiIgkmJIEERERiUtJgoiIiMSlJEFERETiUpIgIiIicSlJEBERkbiUJIiIiEhcShJEREQkrpxNEhoaGgiHw0QikXSHIiIi0i+RSIRwOAwwJpHlmntuztgdDod9zZo16Q5DRETaaRlyNzc/nfrHzNa6ezhR5eVsTYKIiIh0TUmCiIiIxKUkQUREROJSkiAiIiJxKUkQERGRuJQkiIiISFxKEkRERCQuJQkiIiISl5IEERERiavbJMHMPmNm3zSzD7Vbf3XywhIREZF06zJJMLMfAN8CxgGPmNk1bTZflcS4REREJM2GdLP9fOBj7r7fzOYDS8zsaHf/Nu8Nry0iIiJZqLvbDYPcfT+Au28HzgWOM7Of9eBYERERGcC6+6DfYmb/1PLA3d8FZhNMznVSMgMTEZHcE41GYfFiqK5m8eLFwWNJm+6ShCuB+rYr3P2Au38V+GSygkqFhoYGwuEwkUgk3aGIiOQ8d6eiupqisWNh5Up45x1K6+ooGjuWiupq3DVxdFcikQjhcBhgTCLLtb6+8Gb2T+7+P4kMJpXC4bCvWbMm3WGIiAhQUV1N7X330bhsGRQVvbehvp78mTMpu/hiqisq0hfgAGFma909nLDy+pEkPOju5ycqkFRTkiAikhmi0ShFY8fS9MILBycILerrGTpxIls2bmTkyJEpj28gSXSS0F3vhk4N5ARBRHKDumANEEuXwvTp8RMEgKIimqZPJ7RkCcyZk9rYclx34yQUmNkCM3vBzHbFlhfMrNrMDk9VkCIiksW2boXx47veZ9y4YD9Jqe4aLt4NHAp8BjgqtnwGOCS2TUQkY7mWAbFECgspeOWVzi4jAAUbNhApLEx7rJm+JFqXbRLM7CV3/1Bvtw0EapMgIpIZum2TsHkzQ086SW0SeiDRbRK6q0l428zOiRPEucCeRAUhIiK5KxQKUV5WRv7MmVBff/DG+nryL7yQ8rIyJQhp0F3DxcuBn5rZ7cAbsXXvB14DrkhiXCIikkMWzJ0LQM3EiTRNnw7jxlGwYQPNK1ZQVlbWul1Sq0ddIM3sSILkwIA33P0fyQ4s2XS7QUQk80SjUUYtXQpbtxIpLGTWrFmqQeiFtHSBdPdtwLZEnVRERCSeUCjU2s1RnR3TT5M0iYiISFxKEkQk6Wy+YfM1tJHIQNPjJMHMbm77U0RERLJbb2oSPh37eVYyAhEREZHM0ue5G0RERBItGo0Gczls3criwkKKi4uDxoySFjnbJqGhoYFwOEwkEkl3KCIiOc/dqaiqoOiYIvhJKdTNo/THpRQdU0RFVQV9nbE4V0QiEcLhMMCYRJabkpoEMzPPsCt8xBFHoHESREQyQ+X8Smp/XktTSRMMD9btYQ/shtqf1wJQXVWdxggzW0lJCSUlJZjZm4ksN6k1CWY20cxuBF5P5nlERGTgikaj1NxcQ2NxY2uC0Go4NBY3UnNzDTt37kxHeDmtNzUJLYMpdTnaYmx0xi8SDOl8IrAAmNqn6EQkO8xLT0Wiul0OEGuBY+mYILQYDk3HNhG6PASTUxiX9Lwmwd3PbfuzPTO7xMweAv4HOIZgsKwt7n6Du7+ciGBFRCQLvQ2M6mafUbH9JKUS2SbhV8DjwBnu/hqAmWVUOwQRyS2ephoM6Z3FixdT+uPSoA1CJwreKmDR1xcxZ44Ga+6KVSW29iyRbRI+RlBp9Ccze9TMvgwMTmD5IiKShYqLi2ne0Ay7O9lhNzRvaGbWrFkpjUsSmCS4+3Pu/m2CO0s/IGiHEDKzpWamKysiInGFQiHKry0nf2l+x0RhN+Qvzaf82nLNBpkGfbrdYGYzgOvc/Yz229z9ALACWGFmBUAx8C/Akv4EKiIi2WvBvAUA1NxcQ9O4wTC8kYLd+TRvaKbs2rLW7ZJa1tXwBWb2CeB24AMEH/I/AO4iaEIyz93v7rJws4nAFcAX3P0DiQo6EcLhsGucBJHUaLlLqhYC0p1oNMqo2IiLkcJCZs2apRqEXjCzte4eTlR53dUk/BBYBPwRuAh4FvgpUOHu+zoJUF0gRUSkT0KhEMQaJ6qJYvp11yYhz90j7r7e3RcCu939e/EShER1gTSzcWZ2m5k9b2bNZvZ4L54PZjbIzNaYmZvZBb05VkRERN7TXU1Cc7vHXQ33mKgukBOB84BngLw+HP9V4P19OE5ERETa6K4mYVLs23yzmTW3eXwg9ritRHWBXObuH3D3WcALvTnQzELAfwDX9eG8IiIi0kaXSYK7D3L3wW2WQW1/tts3IV0gY70j+qoaeApY2Y8yREREhCTMAhmnC+QsoIQkd4E0s48CVwEfTeZ5REREckVSZ4F09z3ufoe7n53M88T8F/Ajd9/Qk53r6+sxs9alqqoqudGJiIgkQVVVVetnGVCUyLJ7lCSY2Sgz+5WZNZjZVjOLxO7/t93n42ZWaWZPJzLAHsZ3CfAh4PqeHlNUVIS7ty5KEkREZCCqqqpq/SwD6hNZdk9rEhYBnweeBpYBpwMPmNmRZvafZrYVWAVUAv1pU9BrZpYH3AQsBAaZ2Ujem3C0wMwOT2U8IiIi2aKnbRLOAb7n7jUAFtRpRIDVBKMxPkzQBXK5u29PRqBdKCDo8lgbW9r6DfAqMC7FMYmIiAx4PU0SjgSebHng7m5m1cBXgEp373E1fxK8DXy63bpC4NfA94FHUx6RiIhIFuhN74b24yJsjv18JEGxAGBm+QSDKQEcDQw3s+LY44fcvdHMNgBPuPtX3H0/wSBObcs4LvbrX9392UTGJyIikit6kyTcYmZrgRdjy6ux9e8mOKYj6dhdsuXx8cBrBHH3ZaAmERER6aGeJgk/BE4iaLx4JAdP5vajWI+GPwN/Adb3Z0Ck2JDO1s0+x/W3DBEREelaj5IEdy9t+T02y+NHYstHYz+/DgyL7dIIHJbYMEVERCTVej3iortvIxj2uHXo41hvh/EEScNJCYtORERE0qZHSYKZbQEejC2PuPvbbbd7MILDy7FlaaKDFBERkdTr6WBK3yJIKH4KvGlmj5jZNWY2PnmhiYiISDr1KElw99+6+5UE4w9MIRh58TJgvZm9bGaLzGxabPRDERERyQK9muDJA8+6+zx3DxOMY/AD4BjgXmC7md1rZlclIVYRERFJoX7NAunuW939v939c8Bo4LPA34HvJCI4ERERSZ+ETRXt7vvcvc7dS919Qtttsd4PGaWhoYFwOEwkEkl3KCIiIv0SiUQIh8MAYxJZrsWmlkwKM5sIXAF8wd0/kLQT9UE4HPY1a9akOwyRnNDyLSF57zaSTfT30ndmtjbWHCAhElaT0CI2fXSpmf0ZWAvsBKYm+jwiIiKSXAlLEszsEjN7CPgfgoaMc4At7n6Du7+cqPOIiIhIavR6xMUu/IpgNsYzYnMnYGaqLRIRERmgEnm74WMEtxf+ZGaPmtmX0UyNIiIiA1Yiezc85+7fBo4lGDthKhAys6VmNitR5xEREZHU6FOSYGYzzOzJeNvc/YC7r3D3y4D3AcuAf+lHjCIiIpIGXSYJZvYJM3vRzN42s5+b2YfM7P8B/wX8pAflHwdMBD7Y/1BFREQklbqrSfghsAgIA+uBZ4FHgRPd/e54B6gLpIiISHboLknIc/eIu69394XAbnf/nrvva7+jukCKiIhkl+66QDa3e/xmF/uqC6SIiEgW6S5JmGRmbRMFiz02gkkh23Zx/BjB9NF/MrNXgLtQF0gREZEBq8vbDe4+yN0Ht1kGtf3Zbl91gRQREckiCZ+7IU4XyAeAkkSfR0RERJIr4UlCW+6+x93vcPezk3keERERSbweJQlmNsrMfmVmDWa21cwiZhZqt8/HzazSzJ5OTqiJ1dDQQDgcJhKJpDsUyULGe9PdiogkWyQSIRwOA4xJZLnm3n0HBDO7E/gi8CCwDTgN2AV8BpgLfB44gmD672fc/YxEBpkM4XDY16xZk+4wJEu1JAjq3gPRaJRRS5fC1q1ECgspLi4mFAp1f6DkLP3/9J2ZrXX3cKLK6+nthnOA77n7Re4+B/gIsA5YDVxNMDbCFcBRAyFBEJHkc3cqqiooOqYIflIKdfMo/XEpRccUUVFVQU++oIhIevV0qugjgda5Gtzdzawa+ApQ6e7XJyM4ERm4KudXUvvzWppKmmB4sG4Pe2A31P68FoDqquo0Rigi3elNw8X2Ayttjv18JEGxiEiWiEaj1NxcQ2NxY2uC0Go4NBY3UnNzDTt37kxHeCLSQz2tSQC4xczWAi/Glldj699NeFQiknA2P4VNKdcSjJjSPkFoMRyajm0idHkIJicvDJ+nWxoi/dHTJOGHwEkEDRSP5OD2JD+K9Wj4M/AXYL27H0hkkCIywLwNjOpmn1Gx/UQkY/UoSXD30pbfzexIgoaLHwE+Gvv5dWBYbJdG4LDEhiki/ZXKb9WLFy+m9MelQRuEThS8VcCiry9izpw5KYtLRHqnR10guy3EzIDxBEnDSe5e1e9Ck0xdICWZcr0LVzQapeiYoqDRYh5BX6i3Cb4+fBjYB0MjQ9myaQsjR45MZ6iSgXL9/6c/Et0FsjdtEjrlQabxcmxZmogyRWTgCoVCXFt2LQtvXcT+dw7A2VPhtI/Cuufg1kcZcuggri0rVYIgkuESkiSIiLRngw/BisbCw3+AoqL3NtTXY+eegw0+JH3BScaKRqMQG3xrsQbfSruE3G4YiHS7QZIp16tLo9EoRWPH0vTCCwcnCC3q6xk6cSJbNm5UbYIAweBblZXXU1NTS9Pgc6BxHAX5r9DcvILy8jIWLJhLcGdbupKRtxtEJPOl9O116VKYPj1+ggBQVETT9OmEliyBJDZczNUkbSCqrLye2tr7aGp6AQj+bvbsAaintnYmANXVFWmLL1cldRZIEclRW7fC+PFd7zNuXLCf5LxoNEpNTS2NjctoSRDeU0Rj4zJqamo1+FYaKEkQyRGewiVSWEjBK690GU/Bhg1ECguTGoeZloGwjBq1lKam6XRMEFoU0dQ0nVBoSdpjzfQl0ZQkiEjCFRcX07xiBdTXx99h82aaV6xg1qxZqQ1MMtRWgl70XRkX209SKWeThIaGBsLhMJFIJN2hiGSdUChEeVkZ+TNndkwU6uvJv/BCysvKkt5o0V3LQFgikUIKCrqpeSrYQCRSmPZYM3W57bYIkyeHAcYk8n9IvRtEkqCl1i83/7sC7k7l9ddTU1tL0/TpMG4cBRs20LxiBeVlZSyYq9bqEohGoxQVjT2o0eLBNjN06Els2aLeMN1JdO+GnK1JSCeD1LY0F0kDM6O6ooL6v/0Npk2DoUNZNG0aWzZupLqiQgmCtAqFQpSXl5GfPxNof4uqnvz8CykvT37Nk3SkLpAiCabBYA4WCoVauzlqlgbpzIIFcwGoqZlI0+DpsXESNtDcvIKysrLW7ZJaqkkQSRB3p6K6mqKxY2HlSnjnHUrr6igaO5aK6mpy9daeSE+YGdXVFdTX/w0WTYP5Q1m0aBpbtmykulo1T+mimgSRBKm8/npq77vvoFEGg7Fg6qmdGRsMpiL3BoNRzYr0hmqeMosaLqZYNBplVOwNM6I3zKyhYYg70jC70ldq+Nt3arg4QB1UFf3ww9DUxNXLl6squpcsQ5dRS5cGLfh7MAxxumJMtYOG2d3zG/Dr2bPnHpqaXqC29j4qK69PQ1Qi0hu63ZAildXVLFy6hH1tvmm+C1Bfz8JzzwF3qisr0xqj9MPWrTDuhK73GTs2Z4YhbhlmN36XtpZhdidy7bXfyJmaFZGBSDUJKRCNRvnBTTexr/2UuQBFRex7+A/84KabsmJc8mR/a/UMXRYdfjg893zXwf/1r9wyfHjaYszFYXZFpH9ytiZhLSmsgr3zTpgypcuq6P1TphC680741reSFkY23NCw+Rn6zr8K+OPQYHTBeNd582Z4dCXXND/INTuvSXV0Man8C9AwuyLZIONqEsxsnJndZmbPm1mzmT3eg2P+2cx+bmYbzKzRzF4ys3lmNjQFIXfvscfgox/tep+TToJHH01NPJJ47wJFzXDuWXGHIWbG1GD7O2mJDkjtELGZMsyuiPRPxvVuMLOLgB8BzwAnAf9w9zO7OaYG+DhwF/AK8FGgGqhz98/FPabInK8lMPCu/Br44Hnwfx/sfJ8Lz4eXH4IvJC8Mn5f8a52rrZIXL17MNbdeQ+OYffDsYJh2Fkz4CKz/K9Q9Cqc0k/9mHrf82y3MmZP9Hbs0zK70R66+jyRCons3ZGKSMMjdD8R+XwqM6UGSMMbd32y3rgS4DTjO3f/e/phUdoG85ZZbKP3+v8OGVzuvih4/jlv+zw/4VhJvN6RCrv5zR6NRio4poqmkCfKAdcDbwGHAh4F9MDQylC2btuTMh2JFRTW1tffR2LiMgxOFevLzZ1JWdjHV1bk3boR0L1ffRxIh67tAtiQIvTzmzTir/xz72VnLqZS54oorGDzI4JxOqqLPncrgQcYVV1yRngCl30KhEOXXlpO/NB/2AZOBKbGf+yB/aT7l15bnTIIAwTC7ZWUXM3ToRCiYDXYdBQWzGTp0ImVlF2uYXZEBIONqEtrqaU1CJ8d+E1gEvM/dt7XfnurBlObOm8vCWxexv+kAnN2mKvqRRxkydBDf/bdSrp8/8PuN5/I3gAMHDjDlrLN48qkn4PjBcGQzbBsMG5s54xNTeOLRRxk0KOPy8qRrP4DYrFmzcipZkt7L5feR/sr6moREMLNCYC5wV7wEAaC+vh4za12qqqqSGlN1VTXfu7qMQ3FYWwe/Xsgha+s4FOd7V5dRXVWd1PNL8s2b9x/8z+q3YP8L8MpP4KkFwc/96/if1W8xb95/pDvEjJDJX0xEBqKqqqrWzzISXHuedTUJZnYIUAe8H5js7tF4+2XKsMzZ9q0qV78BdN9Qr56hQyfmVEM9d+f6ykpqa2r4xODBfLSxkVfz81nR3ExZeTlzFyzQsMwSV66+jyRComsSsmqcBAvecX4BTAQ+0VmCkEkyOUnLRJn7mbIU6NngQematibVf2rXV1ZyX20tLzQ1vfeq7NlDPTCzthaAimrVoIlksmy73XALcBFwkbuvT3MsB3F3KiqqKSoaC6UrYd47lJbWUVQ0loqK7Ji7IRqNwuLFUF3N4sWLg8c5Q4MHtRWNRqmtqWFZY2OcQZlhWWMjtTU1WTHKqEg2y5okwcz+HbgauMzdn0x3PO1l82Q3B01etXIlvPMOpXV1SZm8KpUDAg3EwYO6HFgoheMyLx01iultaxDaKQKmNzWxJBRKbiwi0i8Z1ybBzPKB82IPrwWGA/Nijx9y90Yz2wA84e5fiR1zKXA3cAfB2AhtveruDe3Pk8o2Cdl+v7qiupra++6jcdmyg8eBqK8nf+ZMyi6+mOqK7O4PPyAGD0rhh2Y1weCSXaW+1wFDgaT+ZWTY+5v0jNok9F0u9G44ElgSW04lGIqm5fGRsX2GAIPbHDM99vNKglH02y7nxzvJ2rWp+2KVzZPdRKNRamprOyYIAEVFNC5bRk1tbdZXK4dCIcrLy8jPnwm0GwuDevLzL6S8vCy9SWAKqy0KIxFeKSjoMpwNBQUURiJJrj4Rkf7IuCTB3V9zd+tkeS22z3HufmWbY67s4pg70vRU2siM+9WWhGXU0qU0TZ/e5eRVTdOnE1qyJCHny2QHDR50yGeBf+eQQz6bk4MHFRcXs6K5uUO61GIzsKK5mVmzZqUyLBHppYxLElJl8uTcu1+dFFu3wvhuEqBx44L9csAhvMshvpeP8CCf5Ad8hAc5xPdyCO+mO7SUCoVClJWXMzM/P069ClyYn09ZeW6NQCkyEOVskpBKxcXFNDevoGM1dIvNNDevSPq3Kk/CEikspOCVbhKgDRuIFBYm5HyZrKXL34vvvMPz777LH4E1777Li++8w321tVxfWZnuEFNq7oIFXFxWxsShQzm/oIDvmTG7oICJQ4dycVkZcxcsSHeIItKNjGu4mCqpHkwpWye7iUajjH7/+/FXXul08qpBH/wg2zdvzupvjdFolLFFRQePCdBGPTBx6FA2bsmdCZ5aZPsAYpJ4arjYdxpMaYBquR9dUzMx1ohxHAUFG2huXkFZWdmAvV+9c+dOvHlfMHnVHx7t0LuBc6dyoHkfu3btSswHQ4Z2a+t+KKX3uvylbaLoNH0hCIVCEJseO/snyRbJLkoSUsTMqK6uoKzsakaNWgpsZdGiacyadVvKvlXZ/CR8wN4LfGAQHPU3+NAJMK3N5FV1j8IpzdDUzHEXHgef7f/pMvWbRWY0TRURSSwlCSkWCoVo+T41Jxu+Vr0FHH0APnkAwvtg3UPw7ENwGPBvwDCgic6bY/RWht4eK1y8mLrSUtizp9N9NhQUMG3Roiy58CKSC9QmIQ1aasyz4aW//PLLueuZX8EXmzvf6e7BfOnUS/nFL36RusBSrLs2CZuBk3K0TQLoHrP0jv5e+i4XBlNKiYaGBsLhMJFIJN2hDGjz58+Hjc2wu5MddgMbm6nO8ol81OVPRNIpEokQDocBxiSyXNUkpEM2VSUAnzzzTJ584Y9wuQeDaLfYDfzCOGPip/jT44+nKbrUaTs18unNTXxkH2wsKNDUyOibofSO/l76LtE1CUoS0iHLkoQDBw4w5ayzePKpJ+D4wXBkM2wbDBubOeMTU3ji0UcZNCh3Kq2i0SijrhgFb0PkCxF1+UNv+tI7+nvpOyUJCaIkIfE2btzI2HnzYPNmvnT00VRXV3PsscemO6y0yNJL3Gd605fe0N9L32mcBMlYxx9/PMQaJ2ZvE0URkdyhJCHFotEowSgJQbe54uLiWLdIERGJRqMQG6FzcWGh3iPTLHduFKeZu1NdUcHYoiJWAu8AdaWljC0qorqigly97SMiAge/R55fWsq/z5un98gMoJqEFGmZ/OegfvR79lAPzKytBaAiy7sJioh0Ru+RmUkNF1Mglyb/aRn62efl5t9VCzVcPJgaoklXcuk9MtnUcDFR1q5N2WRBGTP5jz6xJNfk6LgUA03GvEdKB2qTkAKa/EdEpHN6j8xcuZskTJ4cfLNOwVIYifBKQUGX4WwoKKAwEkluLJIyjuHoW2zapeh/XEuWvEdmw5JguZskpFBxcTErmps7nQhxM7CiuZlZs2alMiwRkYyg98jMpSQhBXJp8h+vChYRkZ7KpffIgSZ3Gy6m2NwFCwCYWFPD9KYmxhFUn61obqasrKx1u4hILmr7Hnn64MF8pLGRjfn5eo9MM3WBTLFoNMrSUaOCERcj2TX5T/vnltMjpakP5EHUBVJ6KhqNMio24mKksDCr3iNTQRM8Jcixxx7rRxxxBCUlJZSUlKT25Fn2AdJ2iuRzYrUkr+T6FMlZdo37S0mC9Ib+XnovEokQiURYu3bt3939uESVm7NJgmaBTJzqigruq61lWWPjQf2c64GZ+flcXFaWeyOlZdk17i+96Utv6O+l71STkCBKEhJDI6V1IouucSLoTV96Q38vfacRF7NBuj44klDln/KR0vShKyKSMkoSpF80UlonlMyISBbQOAm5JBtGShMRkZRRkiD9opHSpEcMNEq1yMCjJEH6RSOliYhkL7VJkH7TaJIiItlJXSAlYbJ5NEnpH/UIld5QF8i+UxdIyVihtt0c5ySkw6OIiKSR2iSIiIhIXEoSREREJC4lCSIiIhKXkgQRERGJS0mCiIiIxJWzSUJDQwPhcJhIJJLuUERERPolEokQDocBxiSyXI2TIImlDvESh/4spDc0TkLfaZwEyWz6FBARyRpKEkQk6VzfDUUGpJxtkyAiIiJdU5IgIiIicSlJEBERkbiUJIiIiEhcShJEREQkLiUJIiIiEpeSBBEREYkr45IEMxtnZreZ2fNm1mxmj/fwuBFm9nMzi5rZLjO728xGJzlcERGRrJWJgylNBM4DngHyenHcb4EPAl8FDgALgfuATyY4PhERkZyQiUnCMne/H8DMltKDySrM7DRgOjDF3f8YW7cZeNbMprl7XTIDFpFuaLhukQEp4243uPuBPhw2A/hHS4IQK+f/ARtj20REZKCYb8EiaZeJNQl9MQFYH2f9i7FtHbz0Epx55sHrPv95+PrXobERzjuv4zFXXhksb74JxcUdt//rv8Ls2bBpE3zpSx23X3stzJwZnPtrX+u4fe5cmDYN/vIXuOaajttvuAFOPx2efhq+//2O22+5BSZNgro6uP76jttvuw0+9CFYtgxuvrnj9rvugg98AO65B37yk47bly6FMWPgjjuCpb2HHoL8fPjxj+G3v+24/fHHg581NfDAAwdvGzYMli8Pfq+uhpUrD94+ejT87nfB7//+77Bq1cHb3/9++OUvg9+vuSZ4Ddv64AehZVbwkhJ4+eWDt0+aFLx+AJddBm+8cfD2006D//N/gt8/9znYvv3g7VOnQkVF8PuMGbB378HbL7gAysuD39v/3YH+9vS3F/yuv73YytceC+J9TH97/fnbS4RsSRJCwM4466PA2HgHvPvuuzzxxOOtj4899jjguIQHJiIikkxVVVX88IcjgUkARYks2zyD7xW2tElw9zO72e8RYI+7X9xu/S+Bse5+evtjwuGwr1mzJoHRiohIIljsVoPPy9zPp0xlZmvdPZyo8jKuTUIfRYERcdaHYttERESkl7IlSVhP/LYHnbVVEBERkW5kS5KwHCg0szNaVphZmKA9wvK0RSUiIjKAZVzDRTPLJxhMCeBoYLiZtbSpfcjdG81sA/CEu38FwN1XmdkK4BdmVs57gyk9qTESRERE+ibjkgTgSGBJu3Utj48HXiOIe3C7fWYDi4D/JqgheQD4ZtKiFBERyXIZlyS4+2tAl6NouPtxcdbtBL4cW0RERKSfsqVNgoiIiCSYkgQRERGJS0mCiIiIxJVxbRJERCTHVcVGWpyX3jBENQkiIiLSCSUJIiIiEpeSBBEREYkrZ5OEhoYGwuEwkZaJ3kVERAaoSCRCOBwGGJPIcjN6quhk0lTRIiKZyWLD6eXox1O/aKpoERERSQklCSIiIhKXkgQRERGJS0mCiIiIxKUkQUREROJSkiAiIiJxKUkQERGRuJQkiIiISFxKEkRERCQuJQkiIiISl5IEERERiUtJgoiIiMSlJEFERETiUpIgIiIiceVsktDQ0EA4HCYSiaQ7FBERkX6JRCKEw2GAMYks1zxHJ+wOh8O+Zs2adIchIiLtmAU/c/TjqV/MbK27hxNVXs7WJIiIiEjXlCSIiIhIXEoSREREJC4lCSIiIhKXkgQRERGJS0mCiIiIxKUkQUREROJSkiAiIiJxKUkQERGRuJQkiIiISFxKEkRERCQuJQkiIiISl5IEERERiUtJgoiIiMSVs0lCQ0MD4XCYSCSS7lBERET6JRKJEA6HAcYkslzzHJ2wOxwO+5o1a9IdhoiItGMW/MzRj6d+MbO17h5OVHk5W5MgIiIiXVOSICIiInEpSRAREZG4lCSIiIhIXEoSREREJK4h6Q5ARESkLcdaf5P0Uk2CiIiIxKUkQUREROJSkiAiIiJxZWSSYGYfNrOVZtZoZvVmtsDMBvfguLCZrTCzHbGlzsxOSUXMIiIi2SbjkgQzCwF1BC1WLgIWANcC87s57gOx44YAX4otQ4BHzOzYZMYsIiKSjTKxd8O/AMOAz7r7boIP+eFAlZndGFsXz/nA4cBn3H0XgJk9DbwJnAf8JPmhi4iIpFVRIgvLuJoEYAbwh3bJwG8IEocpXRyXB+wH9rRZ93ZsnbXfuaGhof+RSspp1s6+y8XXLlue80B5HpkWZ7rjSdP535fIwjIxSZgArG+7wt1fBxpj2zrzu9g+N5vZkWZ2JLAIiAJL2u/85ptvJixgSZ10/9MPZLn42mXLcx4ozyPT4kx3POk+fyJk4u2GELAzzvpobFtc7l5vZp8GHgC+GVu9BTjH3TtUGzQ2Nr5rZoe0WbUFqO9r0JIyY8xMGV7f5OJrly3PeaA8j8TGaR0qgXsr3a9bqs5fxHs1CPsSWXAmJgl9YmbvI6gxWAt8Nbb634AHzez0WG1EK3c/NMUhioiIDCiZmCREgRFx1odi2zrzbYJ2CcXuvg/AzB4FXgHKea92QURERHogE9skrKdd24NY98Z82rVVaGcC8EJLggDg7u8CLwAnJCFOERGRrJaJScJy4BwzO7zNutnAXuCJLo77O3BS23YGZnYocBLwWhLiFBERyWqZmCT8FHgHuNfMpplZCVAF1LbtFmlmG8zsZ22Ou52g8cbvzex8M7sAuI+gMUfSm5ia2Qdio0S+aGYvmNmNZv1vdSPJp2vXd3rtBiZdt9xlZk+Y2XNm9ryZLY2NQ9SpjEsS3D0KTAUGA8sIRlpcBMxrt+uQ2D4tx60FziUYUOku4BcEtyjOdvfnkh85+4HvuvuJwMeAU4DPpuC80n+6dn2n125g0nXLXRe6+8nu/lHgdYL2fJ3KuCTBzGYBPwA+BDQTdEtc5+7Nbfdz9+Pc/cp261a6+6fcfVRsmeLuj3dxrnFmdlsso2o2s6727XI+CXff4u5rYr+/CzwPfKC3z38gM7NiM3vazLabWZOZvWRmc9t1NU3UubL22pnZ0Wb2tpm5mR2WhPKz9rVLBzMbYmbfM7NXzOwdM3vDzBYl4Ty6bjmqp9e+u+sO0GZE4kFAAcEUCJ3KxN4NZcBGoJT3hlT+lZmNcff/SvC5JsbKf4agZ0Rc9t58EusI5pM4AbiZIMmaG2f/0cDFwPQEx5vpRgOPAjcRjHXxcYJbRYXA1Qk+VzZfu5sIRgstSFL52fzapcMdwFkEtZ7rCT5sP5yE8+i65a5ur31vrruZPQT8M0HD/mu7PLO7Z9QCjImz7lfAxi6O+Rjw/k62nQcM7mTboDa/LwUe72S/fyfofjm8zbrvEIzwOLzdvocCjwHXpvu1zIQF+A+ChMF07Xr0en0K2EHQbdeBw7rYV69dmheCW5z7gA/34hhdNy29/Tvr9tr35rrHtg0GbgS+09W5M+52g7vHG53qz3Q9aUUVwURQY9quNLPLCdo1nNXJuQ70MKwezScRq9a5G/izu9/cw7Kz3Xagq9sNVejatY3hvwhmPu3JKG1V6LVLt6uAR919XS+OqULXTXqhh9e+V/MeeXAL/07g8q4KzbgkoROnAS93sf1Kgmz+4ZaWmmZ2EfDfwPfd/ZF+nr+n80ncBrxFd9U3Wc7MBptZvpmdQTCI1U88lrrGcSW6di3+heCb3a093P9K9Nql2ynAy2b2IzPbHbsXfK+ZdfWl5kp03STxur3uZhYys6Pa7PI54H+7KjTjkwQzm0pwv6zTTNeDHhHnEIzK+H/N7DzgHuBmd1+YgDC6nU/CzD4BfAUIA382s7+YWa6O8rgntvyJYGyLTlvP6toFYveFq4EybzMgWFf02mWEQoIP/UnAJcCXgckEXbHjdinUdZMk6cm8RyHggVgDyL8SjCP0ra4KzcSGi63M7DiC9gj3u/sdXe3r7lvMbBrwFPAgcLu7fzfpQb53/qeIMyV1jjqdoPvpx4FK4EfA1zvbWdcOCNpuPOPuD/XmIL12aWex5SJ33w5gZlsIkuOzgJXxDtJ1k3Rw978RNFjssYxNEsxsFMHoi38HvtjDw4YR3P8+ABxmZoN6cR+vK32dTyInufv/xH590oIZ0O40s5vd/dUuDsvZa2dmEwnubX/KzEbGVufHfo4ws2Z339tFETn72mWAKPC3lgQh5kngXYIeDnGThBhdN0mkpFz3jLzdYGb5BFM+HwJc4O6NPTjmOGAFQReRKcAF9Pzebnf6Op+EQEvCcHxnO+jaMZ6gW9Mqgn/mKO89/zcIGjPGpdcu7V4k/jdyI/jwj0vXTZIgKdc942oSzGwIwZTP44HT3X1bD445CngEeBWY5e57Yw2BHjKzqLt/v59hLQe+bWaHu/tbsXU9mU9C4BOxnxvjbdS1A4Jvnp9ut+5c4LsEXeL+Fu8gvXYZ4QFgfmwcl5YeKZ8iSPrijvSq6yZJkpzrnsy+nX1ZCOZZcIJW8ae2Ww7t5JgHgDV07AN8EUEr4nM6OS4fKI4tqwgGlmh5nN9mvxCwheAfexpQQjDYzfXpfr0yaQEeJujfP4NgYJb5sdfpN10co2sX//ldSffjJOi1S/91Gk4wtO0qYCZwKbAJeETXTUsC/866vfbJuu5pf/JxXozXYm+O8ZbjOjnmg8QZhCm27XQ6H8jnuJ6ei+D+4qMEWdkWgpbocQc9ydUl9pr8b+wPcyfBrYZvAHldHKNrF//5XUn3SYJeuwxYgHHAQwQ9eqIEIzCGdN20JGrp6bVPxnW3WMEiIiIiB8nIhosiIiKSfkoSREREJC4lCSIiIhKXkgQRERGJS0mCiIiIxKUkQUREROJSkiAiIiJxKUkQERGRuJQkiIiISFxKEkRynJn9pQf7VCU/EhHJNBqWWUS6ZWbu7vGmRBaRLKaaBJEcZ2be5vdPmtmfzGytmf3RzD5sZoti2/5iZn9KX6QikmqqSRDJcS21BGY2imBO+nPdPWpmnwBq3P001SSI5CbVJIhIi9OA8cBjsXYKtwJH96YAM3vczDy2XNNu28VmtsLMtpvZu2a22cyWmtm5vTzHMjP7axfbf2RmO83sUDOrahPP0t6cR0SUJIjIewYBT7v7pDbLMX0o5zGChOM3LStityx+B2wGvgpMA74HDAOWm9kJvSj/18BJZvbh9hvMbDBQDNzr7u8At8di+XMfnodIzlOSICItngYmm9nHAMxskJn9U2zbW2Y2vIfl7HD3Z9x9a6yci4BrgK+4+5fd/ffu/kd3v8vdzwcuBPb2Is77gUbgC3G2fRo4iiCRwN3fcPdngN29KF9EYpQkiGQxM9thZtVmdnesmn+Hmd1mZoe239fdtwOzgFvN7DngBeAzsc0/BP5fHxsuXgOsdvc74m1092XuXt8u7k+a2RNm1hiLe7GZHR7bfw+wDJgdp7hLgG3Ao32IU0TaUZIgkqXM7BggBJQCOwgSgFuAEuDfYvscEdsGgLs/6e6nu/vJ7n6iu1fE1le4+wR3/2QvYxhCUN2/ohfHfAKoA7YS3Dq4BjgP+Hmb3X4NjDezyW2OywM+C/zW3Zt7E6eIxDck3QGISNKcHPv5I3f/Xuz3R83sU8C5ZvYg8BBwQxJjGA0cCmxqu9LMDBjcZlWzv9fV6gcEbSNmt9l/M7DSzE5y9/8l6IWxk6DmYG1st3MIkqJfJ+F5iOQk1SSIZK+TCe7d39Ru/cvAKHd/yd1PcPebUxBL+77W1wL72iwtNRv5BDUPvzWzIS0L8GRsv8kA7v4ucC/w+VjCAcHth78Dq5L8XERyhpIEkex1MvDHWFuDtgqB+jj7J8N24B3g/e3W3wX8c2xpK0RQw/BjDk4i3gHygA+02ffXwDHAaWY2FLgI+E2bGgkR6SfdbhDJXicDBzU0jDX+OwtY1GbOhlHAYcDrscffdfc/JCIAd99vZquA6UBlm/X/AP4Ri6ntITsJah2qCG6FtNc2uXksVsYlwPuAw9GtBpGEUpIgkoXMrAA4gY41Bt+M/fyhu8+P7XslcKa7X5mkcG4B7jOzL7n7XV3t6O57zOwZ4EPuvqCbfZvN7LcEDTKPBl509+cSFbSIKEkQyVYfIbideIKZVRLUKEwnaAtwmbvvTFUg7n6/md0C3GFmnybovvgmQaPG6bHd3m5zyHcIGikeAJYCbxHcVjgfuM7dX26z76+BbxB01ZyXzOchkouUJIhkp5MJBhA6n6Dr4HXAS8Asd78/1cG4e6mZ/RH4OvAzglsDDQSNDM9z9+Vt9n0y1gNjPkHbhcEEDRIfJnaLos2+q8zsNeA4dKtBJOGUJIhkp5OBv7r788R6BKSQxXoktO3WiLv/Hvh9Twpw92eBHs3p4O7HdxHIIIIaFU1OJdIH6t0gkp1OBtJ1f/6zBD0SvpWm87dVSRDLp9IdiMhApJoEkSwTGzfgI8Av0nD6rxHcSoD3ekukUwR4IPb7jq52FJGOlCSIZJlYFX9PJ2NK9LlfSsd5OxObEyJVY0KIZB3TuCMiIiISj9okiIiISFxKEkRERCQuJQkiIiISl5IEERERiUtJgoiIiMSlJEFERETiUpIgIiIicSlJEBERkbj+P202rey4HhEwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 540x468 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",
    "\n",
    "jet_pt_pp1_in = jets_pp1[:,1]\n",
    "jet_pt_pbpb1_in = jets_pbpb1[:,1]\n",
    "\n",
    "jet_pt_pp2_in = jets_pp2[:,1]\n",
    "jet_pt_pbpb2_in = jets_pbpb2[:,1]\n",
    "\n",
    "jet_pt_pp3_in = jets_pp3[:,1]\n",
    "jet_pt_pbpb3_in = jets_pbpb3[:,1]\n",
    "\n",
    "jet_pt_pp4_in = jets_pp4[:,1]\n",
    "jet_pt_pbpb4_in = jets_pbpb4[:,1]\n",
    "\n",
    "jet_pt_pp5_in = jets_pp5[:,1]\n",
    "jet_pt_pbpb5_in = jets_pbpb5[:,1]\n",
    "\n",
    "\n",
    "bins_11 = [200,250,300,400,500,1000]\n",
    "pt_min = 200\n",
    "pt_max = 1000\n",
    "\n",
    "# Fill Histogram\n",
    "n_pp, pt = np.histogram(jet_pt_pp_in, bins= bins_11 )\n",
    "n_pbpb, pt = np.histogram(jet_pt_pbpb_in, bins= bins_11 )\n",
    "\n",
    "n_pp1, pt = np.histogram(jet_pt_pp1_in, bins= bins_11 )\n",
    "n_pbpb1, pt = np.histogram(jet_pt_pbpb1_in, bins= bins_11 )\n",
    "\n",
    "n_pp2, pt = np.histogram(jet_pt_pp2_in, bins= bins_11 )\n",
    "n_pbpb2, pt = np.histogram(jet_pt_pbpb2_in, bins= bins_11 )\n",
    "\n",
    "n_pp3, pt = np.histogram(jet_pt_pp3_in, bins= bins_11 )\n",
    "n_pbpb3, pt = np.histogram(jet_pt_pbpb3_in, bins= bins_11 )\n",
    "\n",
    "n_pp4, pt = np.histogram(jet_pt_pp4_in, bins= bins_11 )\n",
    "n_pbpb4, pt = np.histogram(jet_pt_pbpb4_in, bins= bins_11 )\n",
    "\n",
    "n_pp5, pt = np.histogram(jet_pt_pp5_in, bins= bins_11 )\n",
    "n_pbpb5, pt = np.histogram(jet_pt_pbpb5_in, bins= bins_11 )\n",
    "\n",
    "# Statistical Errors\n",
    "err_n_pp = np.sqrt(n_pp)\n",
    "err_n_pbpb = np.sqrt(n_pbpb)\n",
    "\n",
    "err_n_pp1 = np.sqrt(n_pp1)\n",
    "err_n_pbpb1 = np.sqrt(n_pbpb1)\n",
    "\n",
    "err_n_pp2 = np.sqrt(n_pp2)\n",
    "err_n_pbpb2 = np.sqrt(n_pbpb2)\n",
    "\n",
    "err_n_pp3 = np.sqrt(n_pp3)\n",
    "err_n_pbpb3 = np.sqrt(n_pbpb3)\n",
    "\n",
    "err_n_pp4 = np.sqrt(n_pp4)\n",
    "err_n_pbpb4 = np.sqrt(n_pbpb4)\n",
    "\n",
    "err_n_pp5 = np.sqrt(n_pp5)\n",
    "err_n_pbpb5 = np.sqrt(n_pbpb5)\n",
    "\n",
    "# bin width\n",
    "dpt = (pt[1:]-pt[:-1])\n",
    "# bin center\n",
    "pt = pt[0:-1] + 0.5*dpt\n",
    "\n",
    "\n",
    "n_ev_pp = 315769\n",
    "n_ev_pp1 = 338726 \n",
    "n_ev_pp2 = 352683 \n",
    "n_ev_pp3 = 369925\n",
    "n_ev_pp4 = 380962 \n",
    "n_ev_pp5 = 389140\n",
    "\n",
    "n_ev_pbpb = 6461\n",
    "sigma = 5.15\n",
    "d_eta = 2*2\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",
    "\n",
    "dn_dpt_pp1 = (n_pp1*sigma)/n_ev_pp1/dpt/d_eta\n",
    "dn_dpt_pbpb1 = (n_pbpb1*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "dn_dpt_pp2 = (n_pp2*sigma)/n_ev_pp2/dpt/d_eta\n",
    "dn_dpt_pbpb2 = (n_pbpb2*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "dn_dpt_pp3 = (n_pp3*sigma)/n_ev_pp3/dpt/d_eta\n",
    "dn_dpt_pbpb3 = (n_pbpb3*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "dn_dpt_pp4 = (n_pp4*sigma)/n_ev_pp4/dpt/d_eta\n",
    "dn_dpt_pbpb4 = (n_pbpb4*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "dn_dpt_pp5 = (n_pp5*sigma)/n_ev_pp5/dpt/d_eta\n",
    "dn_dpt_pbpb5 = (n_pbpb5*sigma)/n_ev_pbpb/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",
    "\n",
    "err_dn_dpt_pp1 = (err_n_pp1*sigma)/n_ev_pp1/dpt/d_eta\n",
    "err_dn_dpt_pbpb1 = (err_n_pbpb1*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "err_dn_dpt_pp2 = (err_n_pp2*sigma)/n_ev_pp2/dpt/d_eta\n",
    "err_dn_dpt_pbpb2 = (err_n_pbpb2*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "err_dn_dpt_pp3 = (err_n_pp3*sigma)/n_ev_pp3/dpt/d_eta\n",
    "err_dn_dpt_pbpb3 = (err_n_pbpb3*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "err_dn_dpt_pp4 = (err_n_pp4*sigma)/n_ev_pp4/dpt/d_eta\n",
    "err_dn_dpt_pbpb4 = (err_n_pbpb4*sigma)/n_ev_pbpb/dpt/d_eta\n",
    "\n",
    "err_dn_dpt_pp5 = (err_n_pp5*sigma)/n_ev_pp5/dpt/d_eta\n",
    "err_dn_dpt_pbpb5 = (err_n_pbpb5*sigma)/n_ev_pbpb/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_pp1\n",
    "error_raa1 = ratio_error(dn_dpt_pbpb1,err_dn_dpt_pbpb1,dn_dpt_pp1,err_dn_dpt_pp1)\n",
    "\n",
    "raa2 =  dn_dpt_pbpb2/dn_dpt_pp2\n",
    "error_raa2 = ratio_error(dn_dpt_pbpb2,err_dn_dpt_pbpb2,dn_dpt_pp2,err_dn_dpt_pp2)\n",
    "\n",
    "raa3 =  dn_dpt_pbpb3/dn_dpt_pp3\n",
    "error_raa3 = ratio_error(dn_dpt_pbpb3,err_dn_dpt_pbpb3,dn_dpt_pp3,err_dn_dpt_pp3)\n",
    "\n",
    "raa4 =  dn_dpt_pbpb4/dn_dpt_pp4\n",
    "error_raa4 = ratio_error(dn_dpt_pbpb4,err_dn_dpt_pbpb4,dn_dpt_pp4,err_dn_dpt_pp4)\n",
    "\n",
    "raa5 =  dn_dpt_pbpb5/dn_dpt_pp5\n",
    "error_raa5 = ratio_error(dn_dpt_pbpb5,err_dn_dpt_pbpb5,dn_dpt_pp5,err_dn_dpt_pp5)\n",
    "\n",
    "err_1 = ratio_error(raa1,error_raa1,raa,error_raa)\n",
    "\n",
    "err_2 = ratio_error(raa2,error_raa2,raa,error_raa)\n",
    "\n",
    "err_3 = ratio_error(raa3,error_raa3,raa,error_raa)\n",
    "\n",
    "err_4 = ratio_error(raa4,error_raa4,raa,error_raa)\n",
    "\n",
    "err_5 = ratio_error(raa5,error_raa5,raa,error_raa)\n",
    "\n",
    "\n",
    "##########\n",
    "\n",
    "plt.errorbar(pt, raa2/raa, ms=9,mfc='red',mec='black',fmt='o', label=r\"R=0.4\",\n",
    "             xerr=0.5*dpt, yerr=err_2, color='red')\n",
    "\n",
    "plt.errorbar(pt, raa3/raa, ms=9,mfc='blue',mec='black',fmt='o', label=r\"R=0.6\",\n",
    "             xerr=0.5*dpt, yerr=err_3, color='blue')\n",
    "\n",
    "plt.errorbar(pt, raa4/raa, ms=9,mfc='green',mec='black',fmt='o', label=r\"R=0.8\",\n",
    "             xerr=0.5*dpt, yerr=err_4, color='green')\n",
    "\n",
    "plt.errorbar(pt, raa5/raa, ms=9,mfc='cyan',mec='black',fmt='o', label=r\"R=1.0\",\n",
    "             xerr=0.5*dpt, yerr=err_5, 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{R}}_{\\mathrm{AA}}/R^{\\mathrm{R=0.2}}_{\\mathrm{AA}}$\")\n",
    "plt.xlim(pt_min,pt_max)\n",
    "plt.ylim(0.8,2.0)\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\"$JS(MATTER+LBT),Q_{\\mathrm{sw}}=2GeV$\")\n",
    "\n",
    "# save plot to the Desktop\n",
    "plt.tight_layout()\n",
    "plt.savefig('cms_double_pT_0-10')\n",
    "print(\"R=0.3/0.2=\", raa1/raa)\n",
    "print(\"R=0.4/0.2=\",raa2/raa)\n",
    "print(\"R=0.6/0.2=\",raa3/raa)\n",
    "print(\"R=0.8/0.2=\",raa4/raa)\n",
    "print(\"R=1.0/0.2=\",raa5/raa)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "a64d3017",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAGICAYAAADMEtD8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAjklEQVR4nO3dfXyU1Z3//9chBJIQY1IiCUFIFISgliJkVVwrYGzZQF1bQRJsfdjqimxXv229WbW1Egs/lKrYlYolbdVaWsJdV0VBVBBsK+6KrrZFsFIJtSCRQgAx4S75/P6YyTSTZJJrJjOZSfJ+8rgeYc51rnN9TmaST667c5yZISIiIm3rFe8AREREugIlTBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREP2k2YzrmrnHPPOud2O+eOOOfedM7N8LBdX+fcQ865j51znzrnnnfOFUQlahERkU7m5QjzFuAI8B3gX4FXgF87525uZ7tHgK8DtwHTgGzgJedcSsTRioiIxIlrb+AC51y2mf29WdmvgXFmdkaIbU4HqoDrzOwpf9kgYCfwTTP7WRRiFxER6TTtHmE2T5Z+/wfktbHZF/1ff9Oknd3A74CScAIUERFJBJHe9DMO+HMb6wuBv5nZkWbl2/zrREREupTe4W7gnCsGvgxc10a1LOBgK+U1/nUiIiJdSlgJ03+X66+BZ8zsyWgHk5aWZnV1dYHXAwcOJC/vH2d+9+3bx2mnneaprT179gRtG0o4bSZCXa/9imUMXut2136FWzcWfeuu/YpV3e7ar3Dr6rPYert79uzho48+anxZa2b9Wt3YzDwtwGfwnVL9XyCtnbo/BHa2Uv4osDXUdmPHjrW2tLe+KV/X2hdOm4lQ12u/YhmD17rdtV/h1o1F37prv2JVt7v2K9y6+iy23y6wxULkKE/XMJ1zacBzQB/gS2ZW284m24HBzrnmWbrQvy4iM2fOjHTTqLSZCHXDEe94u2u/wq0bjlh8b7tSv2JVt7v2K9y64dBnsSUvj5X0Bp4BzgcuMrP32230H4+VfN3MlvjL8vxlIR8rKSoqsi1btoQTf1sx0F7fuiL1q+vprn1Tv7qe7tq3aPbLOfemmRW1ts7LNcxFwGTgW0B/51z/Juv+z8yOOefWA5hZsf/r35xzPwd+5JxzwD6gHNgFLIm4J2GYPXt2Z+ym06lfXU937Zv61fV01751Vr+8HGFWAfkhVp9hZlXOuY0AZjahyXZ9gfuAa4A0YBPwH2a2M9S+onmEKSIiEq4OHWGaWYGHOhNaKTuGb1i9W9oPUUREJLGF/RxmPB0+fJiPP/6YEydOxDsUkS4tOTmZAQMGkJGREe9QRLqMLpMwDx8+THV1NYMGDSI1NRXfpVERCZeZUVdXx+7duwGUNEU86jLzYX788ccMGjSItLQ0JUuRDnDOkZaWxqBBg/j444/jHY5Il9FlEuaJEydITU2Ndxgi3UZqaqoub4iEocskTEBHliJRpJ8nkfB0qYQpIiISL0qYIiIiHiRUwty3bx9FRUVUVFTEO5SYW7VqFZdeeimZmZn07duX4cOHc8stt7Bnzx7Ad7rMOceSJS0HRlqyZElgfVO7du3immuuYciQIaSkpDB48GCuuOIKXn31Vc9xFRQUBNpuXHJzc1vUe/fddykuLiYtLY28vDzuuece6uvrg+qYGfPmzWPw4MGkpqZyySWX8Pbbb3uOBWDRokVMmTKF/v3745xj48aNrdaLZjxe2hKR7qWiooKioiKA7JCVQo3KHo+lrVHk33333TZHmO9KbrnlFuvVq5ddf/319uyzz9rGjRtt0aJFNmrUKPvyl79sZr7R99PT023KlCkttp88ebKlp6cHjdB/4MABy83NtfPOO8+WLFliGzZssMcff9y+/OUv2/333+85tvz8fLv66qtt8+bNgeXNN98MqnPgwAEbOHCgFRcX24svvmiPPfaYpaWl2fe+972gevPmzbOUlBRbuHChvfTSS1ZSUmL9+/e3jz76yHM8F1xwgV144YU2Y8YMA+yVV15pUSea8Xhtq7voTj9XItFAG7OVxD1JNl1inTBzcnIMCLnk5OR0eB/tefbZZw2wn//85y3WnTx50tasWWNmvoQ5Y8YMS05OtgMHDgTq7N+/35KTk+3qq68OSpgVFRXmnLPq6uoW7TY0NHiOLz8/32699dY268ybN88yMzPt0KFDgbL58+dbampqoKyurs4yMjLs3nvvDdQ5cuSIZWdnh5V86uvrzczsj3/8Y8iEGc14vLTVnShhigRrK2Em1CnZWKuuru7Q+mh4+OGHGTNmDNddd12LdUlJSZSUlARejxs3jry8PFatWhUoW7VqFXl5eYwbNy5o24MHD9KnTx8+85nPtGg32ndDrl27lkmTJgU98F5WVkZdXR2bNm0C4LXXXuPw4cNMnz49UKdfv35cfvnlrF271vO+evVq/yMazXi8tCUiPVOPSpjxduLECV577TX+5V/+xVN95xylpaUsXbo0ULZ06VLKyspa1B0zZgzHjh3jmmuu4c0336ShoSHiOH/+85/Tp08fTj31VKZNm8auXbuC1m/fvp3CwsKgsiFDhpCWlsb27dsDdZKSkjjrrLOC6o0cOTJQJ1qiGY+XtkSkZ1LC7ET79+/n2LFjDBkyxPM2ZWVlbNq0ierqavbu3cumTZtaTZjFxcV85zvfYdmyZRQVFZGZmcnUqVN5+eWXw4rxiiuuYNGiRaxfv54HHniAzZs38/nPf55Dhw4F6tTU1JCZmdli26ysLGpqagJ10tPTSUpKalGntraW48ePhxVXW6IZj5e2RKRnUsKMg3BOkZ533nkMGzaM5cuXs3z5coYPH87o0aNbrbtgwQL+/Oc/88ADDzBhwgReeOEFvvjFL/KTn/zE8/7+67/+ixkzZvD5z3+emTNnsm7dOvbs2cMTTzzhuQ0Rke5ICbMT9e/fn759+/LXv/41rO1KS0uprKyksrKS0tLSNusOGzaM2267jWeffZZdu3YxevRovvvd70Y8G/m5557LiBEjeOuttwJlWVlZQUecjWpqasjKygrUOXLkSIvHMWpqakhLS6NPnz4RxdOaaMbjpS0R6ZmUMDtRcnIy//zP/8y6devC2q6srIzNmzfz+uuvt3o6NpTs7Gy+8Y1vUFNT06FBtps/81lYWNjiet6HH35IbW1t4PpfYWEh9fX17NixI6hea9cIOyqa8XhpS0R6JiXMTvbtb3+bLVu28Itf/KLFuoaGBl544YUW5SNHjmTmzJnMnDkz5C/tffv2tVr+/vvv07dvX0499dSI4v3Tn/7E9u3bGTt2bKCspKSEdevW8cknnwTKli1bRmpqKuPHjwfgoosuIiMjgxUrVgTq1NbWsnr16qA7gaMhmvF4aUtEeqhQz5vEY+kJz2Ga+QYuSEpKshtuuMGee+4527hxo/3kJz+x0aNHBw1csHDhwpBtLFy4MOg5zAceeMBGjx5tCxYssPXr19sLL7xgt912myUlJdlNN93kKa7nnnvOysrKAgMfLFq0yPLy8uyMM84IegaxcZCEyy67zF566SVbvHix9evXr9WBAlJTU+3HP/6xvfzyyzZ58mTr37+/7d271/P36o033rAVK1bYggULDLDy8nJbsWKFvfHGGzGJx2tb3YWewxQJhgYuSDwrV660CRMmWEZGhiUnJ9tZZ51lt956a2DUmXAT5tatW+2b3/ymjRw50tLT0+3UU0+1MWPG2KJFi+zEiROeYnrnnXfs0ksvtezsbOvdu7fl5OTYtddea7t3725Rd+vWrTZx4kRLSUmx3Nxcu/vuu+3kyZNBdRoaGmzu3Lk2aNAgS0lJsYsvvtjeeustT7E0uvbaa1v94+baa6+NWTxe2uouutvPlUhHtZUwnUV4M0gsFBUV2ZYtW1pdt23bNkaOHNnJEYl0b/q5EgnmnHvTzIpaW5dQ1zB70uDrIiKSOLwMvt6788Jp32mnnUaoI0zpuJMnT4Zc16tXL0/D0EVTosUjIj1X442Vzrm/h6qj30g9RFVVFcnJySGX1sa27UnxiIi0J6GOMCV28vLyeOONN0Kuz84OPQVcLCRaPCIi7VHC7CH69OnTeH4+ISRaPCIi7dEpWREREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDA70ZNPPhmYKqvp0nyCZzNj3rx5DB48mNTUVC655BLefvvtFu29++67FBcXk5aWRl5eHvfcc0+L+R5bs2PHDm688UZGjRpFUlISEyZMaLVeNOPw2pYXjd+3JUuWtFi3ZMmSFtORRaLxvTpy5Ei7cTjn6NWrF3l5eZSWlrJz507A96xpa+9386WqqqpDsYpI59BjJXGwYcMGUlNTA6/PPPPMoPX3338/c+bM4YEHHqCwsJAFCxZw2WWX8ac//Ync3FzAN6HxZZddxtlnn80zzzzDX/7yF2699VYaGhqYO3dum/vfunUra9as4cILL+TEiRMh60UzDi9thSM9PZ3Kykq+9rWvBZUvXbqU9PT0NhNdNN16661MmzYNM2Pnzp3Mnj2bKVOm8Ic//IGBAweyefPmQN0PPviAr371qzz66KOMGTMmUD5w4MBOiVVEOijUqOzxWLr7bCVPPPGEAfbJJ5+ErFNXV2cZGRl27733BsqOHDli2dnZQVNMzZs3zzIzM4Om3Zo/f76lpqYGlbWmvr4+8P+pU6fa+PHjYxqH17a8AmzGjBmWnJxsBw4cCJTv37/fkpOT7eqrrw6ayaU148ePt9mzZ4dc7+W9opUZZZ5//nkDbOvWrS3q//GPfzTAXnnllTZj60zd4edKJJpoY7aSHnVKtiC3oM1TYwW5BfEOkddee43Dhw8zffr0QFm/fv24/PLLWbt2baBs7dq1TJo0iYyMjEBZWVkZdXV1bNq0qc19eBmjNZpxeG0rHOPGjSMvL49Vq1YFylatWkVeXh7jxo2LqM1oOOWUUwDaPHIXka4poRJmrGcr2VW9C2vj367qXTHZb3NDhw6ld+/ejBgxgsWLFwet2759O0lJSZx11llB5SNHjmT79u1B9QoLC4PqDBkyhLS0tKB6kYpmHF7bCodzjtLSUpYuXRooW7p0KWVlZRG1F6mGhgZOnjzJiRMn+POf/8zs2bM566yzOPfcczs1DhHpGC+zlSRUwmycrWTmzJnxDiUmBg4cyJw5c/jlL3/J6tWrufDCC5k1axYPP/xwoE5NTQ3p6ekkJSUFbZuVlUVtbS3Hjx8P1MvMzGyxj6ysLGpqajocazTj8NpWuMrKyti0aRPV1dXs3buXTZs2hUyY9fX1nDx5MrCYWSDZNS4NDQ1hx/Ctb32L5ORk+vTpw4gRI3j//fdZtWpVi76KSGKbOXNm42xZIWcr0U0/nWjSpElMmjQp8LqkpISjR48yd+5cvvWtb0V9Oqv6+nqsyQThSUlJHb57NJGcd955DBs2jOXLl2NmDB8+nNGjR/O73/2uRd2hQ4eya1fwGYRXX32VOXPmBF7Pnj2b8vLysGK4/fbbA6eaP/74YxYtWsTkyZN5/fXXGTRoUPidEpGEpYQZZ9OmTWP58uVUVVVx5plnkpWVxZEjR6ivrw86SqmpqSEtLY0+ffoAvqOzQ4cOtWivpqaGrKwsAIqLi4OuZ77yyishHyFpLppxeG0rEqWlpVRWVmJmlJaWhqy3evVqjh07Fnh94403Mnbs2KCzGXl5eWHvf8iQIUGDyBcXF3P66afz8MMP8+CDD4bdnogkLiXMOGs84mv8WlhYSH19PTt27GDEiBGBes2vFRYWFra4/vfhhx9SW1sbqLd48WI++eSTwPqm7bUnmnF4bSsSZWVlgaPExx9/PGS9z372s0GvTznlFPLy8qI+Y0rfvn0588wz2bZtW1TbFZH4S6hrmD3RypUryc7OJj8/H4CLLrqIjIwMVqxYEahTW1vL6tWrKSkpCZSVlJSwbt26oIS4bNkyUlNTGT9+POBLkEVFRYGl8Q5OL6IZh9e2IjFy5MjATOkdTb7RcPToUf7yl78wePDgeIciIlGmI8xONHXqVM4//3xGjRpFfX09y5YtY9myZTzyyCOB65cpKSnceeedzJkzh6ysrMBD/g0NDdx8882BtmbNmsUjjzzClVdeyR133MEHH3xAeXk5t9xyS9AjHq2pra1lzZo1AOzevZvDhw+zcuVKACZPnkxaWlpU4/DaVqSaj5QUTU8//TQpKSlBZf/0T/8U+AOnqqqK119/HfDd5f3oo49y6NAhrr/++pjFJCJxEuoBzXgssR64ID8n34CQS35Ofof30Za77rrLhg8fbqmpqZaSkmJjxoyxp556qkW9hoYGmzt3rg0aNMhSUlLs4osvtrfeeqtFva1bt9rEiRMtJSXFcnNz7e6777aTJ0+2G8fOnTtDfg927twZkzi8tuUFrQwY0NTChQujNnBBa8sTTzwRiKPp0r9/f5swYYJt2LCh1TY1cIFI4qONgQucNbmLMhTn3DDgdmAccA7wWzOb4GG7ImAe0Hih6C3ge2b2P63VLyoqMv9tvS1s27aNkSNHthuriHinnyuRYM65N82s1ZsbvF7DPAeYDLwH/NnjTgcDL+M77XuNf+kNvOScy/e4XxERkYTg9RrmajN7BsA5t5I2RkJoYgpwCvAVMzvk3/Y1fA+FTgYeCz9cERGR+PB0hGlm4Q+BAsnASeDTJmVH/GXd5+l5ERHpEWL5WMkqoBZ4yDk3wDk3AHgYqAFWtLmliIhIgolZwjSzPcBEYCpQ7V+uBCaZ2b7WttmzZ0/Q7CHhDlMmIiISrvLy8qYTz4cc8itmCdM5NxDfkeSbQIl/eRN43jk3pLVt8vLygm7hVcIUEZFYKy8vD+QdYE+oerEcuOB2fNcxp5nZCQDn3AbgfeA24P/FcN8iIiJRFctrmIXA1sZkCWBmx4GtwNAY7ldERCTqYpkwdwHnOucCU1E45/oC5wJVMdyviIhI1Hk6JeucS8P37CTAICDDOTfN/3qNmdU653YAm8yscRDNnwH/Bvy3c24RvkdJ/gMYCFREqwMiIiKdwesR5gB8N/CsAC4Ezm7yeoC/Tm8gMNmhmb0J/Au+wQt+CTwFpAFfMLN3ohF8V/Pkk08G3QXcuDQfPNzMmDdvHoMHDyY1NZVLLrmEt99+u0V77777LsXFxaSlpZGXl8c999xDfX19u3Hs2LGDG2+8kVGjRpGUlBRyjsxoxuG1LS8av29LlixpsW7JkiVN73aLWON7lZGRQV1dXYv1l112Gc45vv71r3doP42OHz9OeXl5i+9JVVUVzjmee+65sNrbuHEjzjn+9Kc/RSU+EaFnDb4eb40Dem/YsME2b94cWKqrq4PqzZs3z1JSUmzhwoX20ksvWUlJifXv398++uijQJ0DBw7YwIEDrbi42F588UV77LHHLC0tzb73ve+1G8fTTz9tp59+uk2bNs0KCwtt/PjxrdaLZhxe2vIKsPT0dJsyZUqLdZMnT7b09PR2B19vT+N7lZ6ebitWrAhat3fvXktKSrL09HS79tprO7SfRp988knQwO6Njh49aps3b7aampqw2jt06JBt3rzZamtr26zXHX6uRKKJNgZfj3uSbLrEfraSnHZmK8np8D7a0vhL+JNPPglZp66uzjIyMuzee+8NlB05csSys7ODktC8efMsMzPTDh06FCibP3++paamBpW1pr6+PvD/qVOntpowoxmH17a8AmzGjBmWnJxsBw4cCJTv37/fkpOT7eqrr47abCUzZsywqVOnBq175JFHbMSIETZ27NioJMza2tqQCTPWlDBFgrWVMHvUBNK7qqtDZ0v/+nh77bXXOHz4MNOnTw+U9evXj8svv5y1a9cGytauXcukSZOC5r4sKyujrq6OTZs2tbmPxrk3OysOr22FY9y4ceTl5bFq1apA2apVq8jLy2PcuHERtdmasrIy1qxZEzRBdmVlJaWlpS3qbt++nbKyMgYPHkxaWhrnnHMOP/rRj2ho+MfIko2nStetW8e//uu/kp6ezk033RSY3Psb3/hG4JRyVVVVq6dkCwoKuO2223j44Yc5/fTTycrKoqysjIMHD7bYj07JikRPj0qYiWLo0KH07t2bESNGsHjx4qB127dvJykpibPOOiuofOTIkWzfvj2oXmFhYVCdIUOGkJaWFlQvUtGMw2tb4XDOUVpaytKlSwNlS5cupaysLKL2QpkwYQKZmZk8/fTTAOzatYvNmzczY8aMFnV3797NiBEjWLRoEWvWrOGGG25g9uzZzJ8/v0Xd66+/ns997nM8++yzXH/99WzYsAGAu+++m82bN7N582YGDhwYMq7ly5ezfv16KioqmD9/Ps899xzf/e53o9NpEWlVLAcukGYGDhzInDlzOP/886mvr6eyspJZs2ZRW1vLd77zHQBqampIT08nKSkpaNusrCxqa2s5fvw4ffr0oaamhszMzBb7yMrKoqampsOxRjMOr22Fq6ysjIceeojq6mrMjE2bNrFgwQJ+97vftahbX1/fOIoH4LsU0dDQwMmTJwNlvXr1anH03atXL6666ioqKyu55pprWLZsGaNGjWrxRwJAcXExxcXFgfYvvvhiamtr+elPf8pdd90VVPeqq65izpw5gddHjhwBfH9MXXjhhe32PTk5maeffprevX0/wu+++y6VlZUsWrSo3W1FJDJKmJ1o0qRJTJo0KfC6pKSEo0ePMnfuXL71rW95OlUajuZJIikpqcN3jyaS8847j2HDhrF8+XLMjOHDhzN69OhWE+bQoUPZtWtXUNmrr74alLRmz57d6nCMZWVljB8/ngMHDlBZWRnyKPbo0aPcd999/OpXv+Kvf/0rJ04Exuzg5MmTgeQGMGXKlHC7G2TixIlB7Z199tl8/PHHnDhxguTk5A61LSKtU8KMs2nTprF8+XKqqqo488wzycrK4siRI9TX1wcdkdXU1JCWlhY4EsvKyuLQoUMt2qupqSErKwvwHfE0vZ75yiuvhHyEpLloxuG1rUiUlpZSWVmJmbV6XbHR6tWrOXbsWOD1jTfeyNixY5k5c2agLC+v9TGXG6+Xzps3j7fffpvf/OY3rda74447+NnPfsbs2bMZM2YMmZmZPPPMM8ydO5ejR4+Snp4eqJuTkxNuV4M0P6rv06cPZsaxY8eUMEViJKES5r59+ygqKmLmzJlBv8i6s8YjvsavhYWF1NfXs2PHDkaMGBGo1/xaYWFhYYvrfx9++CG1tbWBeosXLw66WaVpe+2JZhxe24pEWVlZ4Cjx8ccfD1nvs5/9bNDrU045hby8PIqKijztp7S0lAcffJDzzz+fgoKCVuusWLGCm2++mf/8z/8MlD3//POt1u1OR/oi3UFFRQUVFRUA2aHqJNRNP6eddhpbtmzpMckSYOXKlWRnZ5Ofnw/ARRddREZGBitW/GPK0NraWlavXk1JSUmgrKSkhHXr1gUlxGXLlpGamsr48eMBX4IsKioKLI13YnoRzTi8thWJkSNHBv7A6mjybcu1117L5Zdfzi233BKyTl1dHX379g28brxO7UXjUfbRo0c7FqiIRGTmzJls2bIF4O+h6iTUEWas5efk4Np4dCS/g6fJ2jN16lTOP/98Ro0aRX19PcuWLWPZsmU88sgjgeuXKSkp3HnnncyZM4esrCwKCwtZsGABDQ0N3HzzzYG2Zs2axSOPPMKVV17JHXfcwQcffEB5eTm33HJL0CMeramtrWXNmjWA787Ow4cPs3LlSgAmT55MWlpaVOPw2lakmo+UFAtnn3124E7ZUL7whS/w6KOPMmzYMD7zmc/w6KOPBp0GbkufPn0444wzWL58Oeeeey4pKSmMGjUqCpGLSLT0qIRZtXdvXPc/YsQIHn/8cT788EPMjLPPPpunnnqKa665JqjenXfeSUNDA/fddx/79++nqKiIl156Kei6V1ZWFuvXr+emm27i8ssvJzMzk+985zue5hD9+OOPueqqq4LKGl/v3LkzcMoxmnF4aaurW7hwIbNmzeI//uM/SE1N5dprr+UrX/mK5zMmP/nJT7jtttu47LLLOHbsGDt37oxxxCISDtf0Lsp4KyoqMv8hcQvbtm1j5MiRnRyRSPemnyuRYM65N82s1ZsbEuoapoiISKJSwhQREfFACVNERMQDJUwREREPlDBFREQ8UMIUERHxQAlTRETEAyVMERERD5QwRUREPEiohNk4W4l/xPhu58knn8Q512JpPhaqmTFv3jwGDx5Mamoql1xyCW+//XaL9t59912Ki4tJS0sjLy+Pe+65h/r6+nbj2LFjBzfeeCOjRo0iKSkp5JRf0YzDa1teNH7flixZ0mLdkiVLAuujZcKECTjnePnll1usq6qqCnov+/Xrx9ChQ/nqV7/Kb3/72xb1T5w4wYIFCzj33HNJS0sjOzubCy64gPvvvz9q8YpI+CoqKhpnLwo5WwlmljDL2LFjLZR333035Lqu4oknnjDANmzYYJs3bw4s1dXVQfXmzZtnKSkptnDhQnvppZespKTE+vfvbx999FGgzoEDB2zgwIFWXFxsL774oj322GOWlpZm3/ve99qN4+mnn7bTTz/dpk2bZoWFhTZ+/PhW60UzDi9teQVYenq6TZkypcW6yZMnW3p6uvk+2h33t7/9zXr16mWAXXfddS3W79y50wB78MEHbfPmzbZx40Z74oknbNKkSQZYeXl5UP2ZM2dav379bN68ebZ+/XpbtWqV3X777XbRRRdFJd5wdYefK5FoArZYiBwV9yTZdOkpCfOTTz4JWaeurs4yMjLs3nvvDZQdOXLEsrOzg5LQvHnzLDMz0w4dOhQomz9/vqWmpgaVtaa+vj7w/6lTp7aaMKMZh9e2vAJsxowZlpycbAcOHAiU79+/35KTk+3qq69uN2GOHz/eZs+e3e6+HnroIXPO2aWXXmqZmZl27NixoPWNCXP16tUttv3+979vgL3yyitmZvbpp59acnKy/fCHP2xRt6Ghod1YYqE7/FyJRFNbCTOhTsnGWm5uQaunRBuX3NyCeIfIa6+9xuHDh5k+fXqgrF+/flx++eWsXbs2ULZ27VomTZoUNJVXWVkZdXV1bNq0qc19NE4l1llxeG0rHOPGjSMvL49Vq1YFylatWkVeXh7jxo2LqM3WLF26lAsvvJA77riDgwcPhhXv7NmzycvLC5xy//TTTzlx4gS5ubkt6mpCaZHE16MSZnX1LsBCLr71sTd06FB69+7NiBEjWLx4cdC67du3k5SUxFlnnRVUPnLkSLZv3x5Ur/mEyUOGDCEtLS2oXqSiGYfXtsLhnKO0tJSlS5cGypYuXUpZWVlE7bVmx44dbNmyhbKyMoqLixkwYEDQ/tqTlJTEpZdeyuuvvw74JkgfPHgw5eXl/OY3vwmadFtEEl+PSpjxNnDgQObMmcMvf/lLVq9ezYUXXsisWbN4+OGHA3VqampIT08nKSkpaNusrCxqa2s5fvx4oF5mZmaLfWRlZVFTU9PhWKMZh9e2wlVWVsamTZuorq5m7969bNq0KWTCrK+v5+TJk4HFzGhoaAgqa2hoCNpm6dKl9OrVi+nTp5OUlMRVV13F6tWr+fTTTz3HePrpp1PdZNLyJ598kk8++YSpU6eSmZlJUVERDz74YMTfAxHpPEqYnWjSpEncfffdfPGLX6SkpIRf/OIXTJ8+nblz57b4ZR0NrSWJ7uS8885j2LBhLF++nOXLlzN8+HBGjx7dat2hQ4eSnJwcWF599VXmzJkTVPaDH/wgaJulS5cyfvz4wCnUsrIyamtreeaZZzzH2Px7fumll/KXv/yFyspKrrvuOvbv38/tt9/OpZdeGpPPgIhEjxJmnE2bNo0DBw5QVVUF+I66jhw50uKxjJqaGtLS0ujTp0+g3qFDh1q0V1NTQ1ZWFgDFxcVBCaG9a5tNRTMOr21ForS0lMrKSiorKyktLQ1Zb/Xq1bzxxhuBZcyYMdxwww1BZTNnzgzUf+edd9i2bRtf+tKXOHjwIAcPHuScc85h4MCBYZ2W3b17Nzk5OUFlp5xyCqWlpfz0pz/lgw8+4Pvf/z6///3vWb16dfjfABHpNL3jHUBP13izR+PXwsJC6uvr2bFjByNGjAjUa36tsLCwsMX1vw8//JDa2tpAvcWLFwddJ2vaXnuiGYfXtiJRVlbGnDlzAHj88cdD1vvsZz8b9PqUU04hLy+v8bmrFhqT4q233sqtt94atG7dunUcOHCAz3zmM23GdvLkSTZs2MAll1wSso5zjttvv505c+awfft2rrjiijbbFJH40RFmnK1cuZLs7Gzy8/MBuOiii8jIyGDFihWBOrW1taxevZqSkpJAWUlJCevWrQtKiMuWLSM1NZXx48cDvgRZVFQUWE455RTPcUUzDq9tRWLkyJHMnDmTmTNndjj5NjIzKisrmThxIq+88krQ8utf/5oTJ04E3Z0byg9+8AP27NnDrFmzAN+gBQcPHmxR7/333wdocSQqIolFR5idaOrUqZx//vmMGjWK+vp6li1bxrJly3jkkUcCj3qkpKRw5513MmfOHLKysigsLGTBggU0NDRw8803B9qaNWsWjzzyCFdeeSV33HEHH3zwAeXl5dxyyy1Bj3i0pra2ljVr1gC+U4aHDx9m5cqVAEyePJm0tLSoxuG1rUg1HympozZv3syuXbuYP39+q6Mg3XfffSxdupQbbrghUPbee++RnZ3N8ePH2blzJ5WVlbzwwguUl5cH/nA4dOgQw4cP59prr2XixImceuqpvPfee9x3330MGjSIr3zlK1Hth4hEWagHNOOxxHrggpyc/NDPlIDl5OR3eB9tueuuu2z48OGWmppqKSkpNmbMGHvqqada1GtoaLC5c+faoEGDLCUlxS6++GJ76623WtTbunWrTZw40VJSUiw3N9fuvvtuO3nyZLtxND5s39qyc+fOmMThtS0vAFu4cGHI9QsXLuzQwAU33XSTZWRkWG1tbavr58+fb7169bI9e/a0+F6mpKTYGWecYVdffbW9+uqrQdsdO3bM7rvvPvv85z9vAwYMsJSUFBs6dKjdeOON9uGHH7bd6RjRwAUiwWhj4AJnCXTnZFFRkW3ZsqXVddu2bWPkyJGdHJFI96afK5Fgzrk3zazVmxsS6hpmdx98XUREEpOXwdcT6hrmaaedRqgjTBERkVhpvHnQOff3UHUS6ghTREQkUSlhioiIeKCEKSIi4kGXSpiJdEevSFennyeR8HSZhJmcnExdXV28wxDpNurq6khOTo53GCJdRpdJmAMGDGD37t3U1tbqL2ORDjAzamtr2b17NwMGDIh3OCJdhqfHSpxzw4DbgXHAOcBvzWyCx22vBO4CzgVqgTeAqWbmfVJBCAyztmfPHk6cOBHOpiLSTHJyMjk5Oe0Ooygi/+D1OcxzgMnA64DnczjOuX8Dfgz8EF/CzQIuDWO/QTIyMvQDLiIiceE1ca02s2cAnHMraWMkhEbOuWzgYeBmM/tpk1X/HXaUIiIicebpGqaZRTIV/HT/119EsK2IiEhCieVNPxcA7wHXO+f+5pw74Zz7H+fcRTHcp4iISEzEMmHmAiOAu4E7gMuBT4EXnHOaKVdERLqUWCZMB6QD15vZr8zsBeDLQD1wU2sb7NmzB+dcYCkvL49heCIiIlBeXh7IO0BeqHphz4fZeNNPe4+VOOeWAVcBaWZ2tEn5y8AhM5vafJu25sMUERGJtXjNh7kN31Gmax4PEMlNRCIiInETy4T5nP/rxMYC59ypwFjgnRjuV0REJOq8jvSThm/gAoBBQIZzbpr/9Rozq3XO7QA2mdn1AGa2xTn3DPBz59ydwN+B/wROAI9GsxMiIiKx5nXgggHAimZlja/PAKr8bSU1q/M14AFgAZAG/B641MxqIglWREQkXjwlTDOrouW1yOZ1ClopOwL8u38RERHpsrrMbCUiIiLxpIQpIiLigRKmiIiIB0qYIiIiHihhioiIeKCEKSIi4kFCJcx9+/ZRVFRERUVFvEMREZEepKKigqKiIoDsUHXCHnw9ljT4uoiIxFO8Bl8XERHpNpQwRUREPFDCFBER8UAJU0RExAMlTBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDBFREQ8UMIUERHxQAlTRETEg4RKmJqtRERE2lOQW4BzLuRSkFsQdpuarURERLod5xxG6NzlcESa2zRbiYiISAcpYYqIiHighCkiIuKBEqaIiIgHSpgiIiIeKGGKiIh4oIQpIiJdSn5OPq6Nf/k5+THZb++YtCoiIhIjVXur4rJfHWGKiIh4oIQpIiLigRKmiIiIBwmVMDX4uoh0tlgM5C1djwZfFxFpRywH8pauR4Ovi4iIdJASpoiIiAdKmCIiIh4oYYqIiHighCkiIuKBEqaIiIgHnhKmc26Yc26xc+4Pzrl659zGcHbinOvlnNvinDPn3JciilREJAbiNZC3dD1eB18/B5gMvA4kR7CffwNOj2A7EZGYitdA3tL1eD0lu9rMBpvZVcDWcHbgnMsC/j/ge+EGJyIikig8JUwza+jAPuYAvwfWd6ANERGRuIrpfJjOuVHAdcCoWO5HREQk1mJ9l+xC4MdmtiPG+xEREYmpmCVM51wZMAKY63WbPXv2BM0SUF5eHqvwREREACgvLw/kHSAvVL2wZytxzq0Ess1sQht1koEPgAXAE/7iIcA7QBmwxsw+ab6dZisREZF4isdsJf3wPUayAKjxL+/411UC/xej/YqIiMRErG76OQJMbFaWCywFvgtsiNF+RUREYsJTwnTOpeEbuABgEJDhnJvmf73GzGqdczuATWZ2vZmdBDY2a6PA/98/mtn/dDhyERGRTuT1CHMAsKJZWePrM4Aqf1tJ0QlLREQksXhKmGZWBbh26hR0tA0REZFEpdlKREREPFDCFBER8UAJU0RExAMlTBEREQ+UMEVERDxQwhQREfFACVNERLqUgtzcoIk6mi8Fubkx2W9CJcx9+/ZRVFRERUVFvEMREZEEtau6GoOQy67q6rDbrKiooKioCCA7VJ2wZyuJJc1WIiIi7XHO0VbmckCkuS0es5WIiIh0K0qYIiIiHihhioiIeKCEKSIi4oESpoiIiAdKmCIi0qXk5+TgIOSSn5MTk/16nUBaREQkIVTt3RuX/eoIU0RExAMlTBEREQ+UMEWkR4vXuKTS9egapoj0aI3jkobiIhiXVLqnhDrC1ODrIiISDxp8XUSkHbEcyFu6Hg2+LiIi0kFKmCIiIh4oYYqIiHighCkiIuKBEqaI9GjxGpdUIpebW9Dms7O5uQUx2a+ewxSRHi1e45JK5Kqrd0Eb9zZXV7uY7FdHmCIiIh4oYYpEWUE7p4sKYnS6SERiS6dkRaJsV/UurI3TRS5Gp4tEJLZ0hCkiIuKBEqaIiIgHSpgiIiIeJFTC1GwlIiLSnpycfEI/Oev868Oj2UpE4sA3+0UbN/3gNPuFSILSbCUiIiIdpIQpEmX5Ofm4Nv7lR3C6SETiTwlTJMqq9lZhZiGXqr1V8Q5RmojXuKTS9WjgAhHp0eI1Lql0PTrCFBER8UAJU0RExANPCdM5N8w5t9g59wfnXL1zbqOHbf7JOfeEc26Hc67WOfeec262cy6lw1GLiIh0Mq/XMM8BJgOvA8ketykFhgLzgfeBUcAc/9ep4YUpIiISX14T5mozewbAObeSNkZCaOJ+M/t7k9cbnXNHgcXOuXwz2xVmrCIiInHj6ZSsmTWE23CzZNno//xf88JtT0REJJ46+6afcUAD8JdO3q+ISKtiMS6pdE+dljCdc7nA3cAvzezj1urs2bMn6IHh8vLyzgpPRHqove0MNLFXA010e+Xl5YG8QxtnQMMefL3xGqaZTQhjmz7Ay8DpwFgzq2mtngZfFxGReGpr8PWYj/TjfCn7KXx32v5zqGQpIiKSyDpjaLwfAVcAXzCz7Z2wPxERkaiLacJ0zt0F3ARMN7PfxXJfIiIiseQpYTrn0vANXAAwCMhwzk3zv15jZrXOuR3AJjO73r/N1cA84Elgt3PuwiZN/sXM9kWjAyIiIp3B6xHmAGBFs7LG12cAVf62kpqs/6L/69f9S1PfwJdIRUREugRPCdPMqvA9lNRWnYJmr79Oy0QpIiLSJWm2EhEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDBFREQ8SKiEuW/fPoqKiqioqIh3KCIi0oNUVFRQVFQEkB2qTtizlcSSZisREZF4amu2koQ6whTpDgpyc4PmdW2+FOTmxjtEEYlAZ8xWItKj7Kqupq3zNq66utNiEZHo0RFmgivILWjnaKUg3iGKiPQIOsJMcLuqd2FtHK+46jaH+BURkSjREaaIiIgHSpgiIiIeKGGKiIh4oIQpIiLigRKmSJTl5+TgIOSSn5MTx+hEJFK6S1Ykyqr27o13CCISA93mCLO7Pq+Yn5OPa+Nffk5+vEMUEekRus0RZnd9XrFqb1W8QxARERLsCFOzlYiISDz0qNlKnHNtH2HiSKS+iohI4tFsJSIiIh2khCkiIuKBEqaIiIgHSpgiUZbbziNOuV30ESeRnq7bPFaSn5Pf5qMjel5ROkt19S5o4wa06i76iJNIT9dtEqaeVxQRkVjSKVkREREPlDBFREQ8UMIUERHxQAlTRETEAyVMERERDxIqYWrwdekOcnLyCT19tPOvF5FE0qMGXxcREekoDb4uIiLSQUqYIiIiHihhiohIl5Kbm9vOeM25MdmvEqaIiHQp1dXVHVofKSVMERERD5QwE1xBO6ceCmJ06kFERIJ5SpjOuWHOucXOuT845+qdcxs9bneqc+4J51yNc+6Qc+5Xzrn+HYq4h9lVXY1ByGVXjE49iIhIMK/Te50DTAZeB5LDaH85MBz4N6ABmA88DXw+jDZERETizusp2dVmNtjMrgK2etnAOTcO+CJwrZmtMrP/Br4GXOycuyyycEOL111TIiLSM3hKmGbWEEHbJUC1mb3apJ3/BXb610VVvO6aEhGRniGWN/0UAttbKd/mXyciIhK2nJycDq2PlNdrmJHIAg62Ul4DnNnaBu+9BxMmBJdNnw7f/CbU1sLkyS23+frXfYsXzdsGuPVWuPxy375vvLHl+rvvhssug7ffhm9/u+X6efPgoovgtdfgu99tuf5HP4LRo+Hll2Hu3JbrFy+GESNg9Wp46CFv/WjNk0/6lubWrIG0NFi0CJYvb7l+40bf1wcfhOeeC16Xmgpr1/r+P2cOrF8fvL5/f1i1yvf/u+6CzZuD159+OixZ4vv/t7/t+x42NXw4NI6zP3Mm/PnPwetHj/Z9/wC+9jX429+C148bB/fd5/v/1Kmwf3/w+uJi+P73ff8vKYG6uuD1X/oS3Hab7/+tfTa8fvb+/neYNq3l+n//dygthQ8/hGuuabk+0T97v/wlDB4My5bBY4+1XL9yJWRnd/3PXlpaLnV1oc9ApaTkcMEFe4PK9NmL72fvT3/a2+5nDyL77LUllgkzbMePH2fTpo2B1/n5BUBBnKIRkZ6grWQJcPSoLud0d1VVVezaVdX4Mi9UvbBnK3HOrQSyzWxCO/WWA6eZ2cRm5c8DmNmU5tt0ZLaSFOc41sb6vsDRBJqZxauC3Nw2Hx3Jz8mhau/ekOtFpG3OuXbrJNKsThJbbc1WEssjzO20/vhIIb5HS6LqGL7nEkNp/0ciMSkZiogkhlje9LMWyHXOXdxY4Jwrwnf9cm0M9ysiIhJ1no4wnXNp+AYuABgEZDjnGi81rzGzWufcDmCTmV0PYGabnXMvAk85527jHwMX/M7MXo5qL0RERGLM6ynZAcCKZmWNr88AqvxtJTWrUwo8DDyO72j2OeD/RRKoiIhIPHlKmGZWRTuXAc2soJWyg8A3/IuIiEiXpdlKRKRHi9dD8NL1JNRzmCIinW2v7kQXj7rNEWZ+Tg4OQi75XfSvxNzcgnYGlS+Id4giIj1CtznC7K7PK1ZX76KtJ0yrq7vqE6bdV25ubpuD/efk5OioRqQL6jZHmCKJQjPniHRPSpgiIiIeJFTC3LdvH0VFRVQ0TiMgIiLSCSoqKigqKgLIDlUn7MHXY6kjg693V76BodseJTeR3kPRYN4iXVlbg68n1BGmiIhIouo2CVOPX4iISCx1m8dKuu/jF31pe1TCvp0ViIhIj9ZtjjC7r7amxfayXjqbhloT6Z66zRGmSKLQoAQi3ZOOMEVERDxQwhQREfFACTPB6XqYiEhi0DXMBKfrYSIiiaHbHGHm5OQTenIv518vIiISmW5zhLl3b1W8QxARkW6s2xxhioiIxFJCJUzNViIiIvGg2UpERETCoNlKREREOkgJU0RExAMlTBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDBFREQ8UMIUERHxQAlTRETEAyVMERERD5QwRUREPEiohKnZSkREJB40W4mIiEgYNFuJiIhIBylhioiIeKCEKSIi4oESpoiIiAdKmCIiIh54SpjOubOdc+udc7XOuT3OuR8455I8bFfknHvROXfAv7zsnLug42GLiIh0rnYTpnMuC3gZMOAK4AfArcC97Ww32L9db+Aa/9IbeMk5l9+xsEVERDpXbw91ZgGpwJVmdhhfwssAyp1zP/SXtWYKcArwFTM7BOCcew34OzAZeKzD0YuIiHQSL6dkS4B1zRJjJb4kOr6N7ZKBk8CnTcqO+MtcmHGKiIjElZeEWQhsb1pgZn8Fav3rQlnlr/OQc26Ac24A8DBQA6yILFwREZH48JIws4CDrZTX+Ne1ysz2ABOBqUC1f7kSmGRm+8KOVEREJI5i9liJc24gviPJN/Gd1i3x//9559yQ1rbZs2cPzrnAUl5eHqvwREREACgvLw/kHSAvVL12B193zn0MPGpm9zYr/xQoN7MHQmy3AN8R5VlmdsJf1gd4H3jGzP5f8200+LqIiMRTRwdf306za5X+R0bSaHZts5lCYGtjsgQws+PAVmCoh/2KiIgkDC8Jcy0wyTl3SpOyUqAO2NTGdruAc/1HlQA45/oC5wJV4YcqIiISP14S5k+AY8BvnHOXOedmAuXAgqaPmjjndjjnft5ku5/hOxf83865Kc65LwFPAwMBzRAtIiJdSrsJ08xqgGIgCViNb4Sfh4HZzar29tdp3O5N4F/wDV7wS+ApfKdxv2Bm70QjeBERkc7iZaQfzOxd4NJ26hS0UrYeWB9RZCIiIglEs5WIiIh4oIQpIiLigRKmiIiIB0qYIiIiHihhioiIeKCEKSIi4oESpoiIiAcJlTD37dtHUVERFRUaCEhERDpPRUUFRUVFANmh6rQ7W0ln0mwlIiISTx2drURERKTHU8IUERHxQAlTRETEAyVMERERD5QwRUREPFDCFBER8UAJU0RExAMlTBEREQ+6bcIsLy+PdwgxoX51Pd21b+pX19Nd+9Zp/TKzhFnGjh1rbVm8eHGb65vyda194bSZCHW99iuWMXit2137FW7dWPStu/YrVnW7a7/CravPYvvtAlssRI6Ke5JsurSXMNtb36zTnuqF02Yi1A3ngxHveLtrv8KtG4u+ddd+xapud+1XuHX1WWy/3bYSZkKNJeuc+wTfaeJ9wN9bqZIdorw1ecAeD/XCaTMR6nrtVyxj8Fq3u/Yr3Lqx6Ft37Ves6nbXfoVbV5/F0O1mA6cBDWZ2SmsbJlTCFBERSVTd9qYfERGRaFLCFBER8aDLJUzn3NnOufXOuVrn3B7n3A+cc0ntbHOOc+4Ff/1jzrm/Oud+5pwb2FlxtyeSfjXbvpdzbotzzpxzX4plrOGI8P0q8Pej+VLZWXF70ZH3zDl3pXPuDedcnXNuv//z2S/WMXsR4XtWHuI9M+fcXZ0Ve1sifb+cc0XOuRedcwf8y8vOuQs6I2avOtC3c/x9q3XO/d0595hzLr0zYvbCOTfMObfYOfcH51y9c26jx+1Odc494Zyrcc4dcs79yjnXv6Px9O5oA53JOZcFvAy8C1wBDAUewpf4725j01OBncBT+C4MnwHMBsY65/7JzE7GMu72dKBfTf0bcHpMAoxQFPp1G/D7Jq+93gAQcx3pm3Pu34AfAz8EbgeygEtJgJ/HDvTrZ8ALzcq+DNwBrI16oGGKtF/OucH+7d4CrvEX3w685Jz7rJntimXcXnSgb6cCG4A/A6VAf3yfyYH43rtEcA4wGXgdSA5ju+XAcHy/FxuA+cDTwOc7FE2o22cTcQHuAmqAjCZl/wnUNi3z2NYXAAPGdPV+4fuFuw+43t+nL8W7Tx3pF1CQSP2Ict+ygU+AG+Ldh2j2K0RbzwPb4t2nDr5fs4B64NQmZVn+sn+Pd7862Le7gMNAZpOyy/0/e0Xx7pc/nl5N/r8S2Ohhm3H+PlzSpOx8f9llHYmnq52SLQHWmdnhJmWVQCowPsy29vu/9olGYB3U0X7NwXcktj4GsXVENN+vRBNp36b7v/4iVoF1UFTeM//pry8AS6MbXsQi7VcycBL4tEnZEX+Zi3aQEYq0b6PxPXN4sEnZS/gSy5QoxxgRM2uIYLMSoNrMXm3Szv/iO8tY0pF4ulrCLAS2Ny0ws7/i+0uqsL2N/df5+jjnRgD3A28A/xuLQMMUcb+cc6OA6/Cdvkw0HXq/gCf81y0+cs4tcM6lxiLICEXatwuA94DrnXN/c86dcM79j3PuotiFGpaOvmeNpuJLNomSMCPt1yp/nYeccwOccwOAh/Ed0a2IUazhirRvKcDxZmUn8Z3CHBnNADtZi++H3zbC+wy30NUSZhZwsJXyGv+69qwBjuH7Zn4G3ym/SP6CibaO9Gsh8GMz2xHtoKIg0n4dAx7Fd4q5GFgM/Du+v5oTRaR9ywVG4Lu2dAe+U2CfAi8453KiHGMkOvoz1qgMeMvM3o9GUFEQUb/MbA8wEd8fANX+5Upgkpnti36YEYn0PdsBfM451/Ta4FggCd/vx64qWp/hFrpawuyom4EL8V28TwfWOudS4htS5JxzZfh++c6NdyzRZGYfmdlNZvasmW00s3LgFuBfnXOfi3N4HeXwffauN7NfmdkL+G6wqAduimdg0eK/+3w8iXN0GTF/X1YAb+I7nVfi///zzrkh8YwtCn6Kb2Sbhc65XOfcOcAifJ/FRDiQSDhdLWHW4Lvjtbks/7o2mdn7ZvY/ZrYEmAScB1wd3RAjEna//H8VPoDv7q9ezrlMIMO/up9zrtWhnTpZh96vZlb6v47tUETRE2nfavBdI9rYWOC/9vQmcHYU44tUNN6z6fj+MFgWraCiINJ+3Y7v1PI0M3vB/wfOVHxJJVEug0TUNzPbDswEZgAfAX/Ad4nqbWBv1KPsPNH8vROkqyXM7TQ7B+2/7TuN1s9Zh2S+28EPAGdGLbrIRdKvfvgeI1mA70NQA7zjX1cJ/F9MIg1P1N4vfEmm6dd4i7Rv2/Alk+Y3jDgS46/6aLxnZcDvzOzDKMfWEZH2qxDYamYnGgvM7DiwFd/jG4kg4vfMzB4HcoBR+MZjvQkYhu8xjq6qxffDL9S1Tc+6WsJcC0xqdvRUCtQBm8JpyH/jT398d07FWyT9OoLv2krTZYZ/3XeBr8Ym1LBE7f0Cpvm/vhmNwKIg0r495/86sbHA/zzcWP7xB088deg9c84V4LvskWinYyPt1y7gXOdc4G5651xf4FygKgZxRqJD75mZHTWzP5pZNfA1fHlheUwi7RxrgVzn3MWNBc65InwHRx17Jjjez9mE+UxOFr5TBy8Bl+E7nXAEmNus3g7g501eP4jvrtiv4PtF9U18H/YdQL+u2q9W2ikggZ5f7MD7VY7vwesr/dv9AN8P/6p49yka7xm+B6g/Aq7Fd/v+JnzP0WZ15X75y+8ETgDZ8e5LlD6LY/39ed7/Xn0J3y/dE8Dn4t2vDvYtA98lnSn4LlHd7+/X1+PdpyYxpuH7Y3kasBnfkX3j67RQn0VgHfCB/3fIl/Hdmf7bDscT729IBN/As/GNTlHn/5DMAZKa1akCnmzyugzfc4oH8N1qvd3/Czlhfqgj6VcrbRSQQAmzg+/XFuAQvtved+BLmn3j3Z9ovGf4bvp5DN+zwHX4Rmn5bLz7E43PIr7rXy/Euw9Rfr+KgVf9vz8O4PsDZ0K8+9PRvuG7rPOiv091+B6z+3K8+9Is5sbfaa0tBW28Z5nAE/julj0M/Doav+81vZeIiIgHXe0apoiISFwoYYqIiHighCkiIuKBEqaIiIgHSpgiIiIeKGGKiIh4oIQpIiLigRKmiIiIB0qYIiIiHvz/pz7VZ1uSbUEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 540x468 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure size 540x468 with 0 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R03 = raa1/raa\n",
    "R04 = raa2/raa\n",
    "R06 = raa3/raa\n",
    "R08 = raa4/raa\n",
    "R10 = raa5/raa\n",
    "\n",
    "\n",
    "#print(R03)\n",
    "#print(R04)\n",
    "#print(R06)\n",
    "#print(R08)\n",
    "#print(R10)\n",
    "\n",
    "R_bins = [0.3,0.4,0.6,0.8,1.0]\n",
    "\n",
    "\n",
    "R_CMS_300_400 = [1.019,0,0,0,0]\n",
    "R_CMS_400_500 = [1.038,1.048,1.103,1.12,0]\n",
    "R_CMS_500_1000 = [0.998,0.933,0.948,0.987,1.02]\n",
    "#R_CMS_500_1000_err = [0.01,0.01,0.01,0.01,0.01]\n",
    "\n",
    "R_200_250 = [R03[0], R04[0], R06[0], R08[0], R10[0]]\n",
    "R_250_300 = [R03[1], R04[1], R06[1], R08[1], R10[1]]\n",
    "R_300_400 = [R03[2], R04[2], R06[2], R08[2], R10[2]]\n",
    "R_400_500 = [R03[3], R04[3], R06[3], R08[3], R10[3]]\n",
    "R_500_1000 = [R03[4], R04[4], R06[4], R08[4], R10[4]]\n",
    "\n",
    "R_500_1000_M = [0.9*R03[4], 0.9*R04[4], 0.9*R06[4], 0.9*R08[4], 0.9*R10[4]]\n",
    "R_500_1000_A = [0.85*R03[4], 0.85*R04[4], 0.85*R06[4], 0.85*R08[4], 0.85*R10[4]]\n",
    "\n",
    "\n",
    "#plt.errorbar(R_bins, R_CMS_300_400,ms=8,mfc='black',mec='black', fmt='s', label=\"CMS_300_400\", xerr=None, yerr=None, color='black')\n",
    "#plt.errorbar(R_bins, R_CMS_400_500,ms=8,mfc='black',mec='black', fmt='s', label=\"CMS_400_500\", xerr=None, yerr=None, color='black')\n",
    "plt.errorbar(R_bins, R_CMS_500_1000,ms=8,mfc='black',mec='black', fmt='s', label=\"CMS_500_1000\", xerr=None, yerr=None, color='black')\n",
    "\n",
    "\n",
    "#plt.errorbar(R_bins, R_200_250,ms=8,mfc='yellow',mec='black', fmt='s', label=\"200-250 M+LBT\",xerr=None, yerr=None, color='black')\n",
    "#plt.errorbar(R_bins, R_250_300,ms=8,mfc='blue',mec='black', fmt='s', label=\"250-300  M+LBT\",xerr=None, yerr=None, color='black')\n",
    "#plt.errorbar(R_bins, R_300_400,ms=8,mfc='green',mec='black', fmt='s', label=\"300-400  M+LBT\",xerr=None, yerr=None, color='black')\n",
    "#plt.errorbar(R_bins, R_400_500,ms=8,mfc='red',mec='black', fmt='s', label=\"400-500  M+LBT\",xerr=None, yerr=None, color='black')\n",
    "plt.errorbar(R_bins, R_500_1000,ms=8,mfc='magenta',mec='black', fmt='s', label=\"500-1000  M+LBT\",xerr=None, yerr=None, color='black')\n",
    "plt.errorbar(R_bins, R_500_1000_M,ms=8,mfc='Red',mec='black', fmt='s', label=\"500-1000  M+Martini\",xerr=None, yerr=None, color='black')\n",
    "plt.errorbar(R_bins, R_500_1000_A,ms=8,mfc='Blue',mec='black', fmt='s', label=\"500-1000  M+ADS\",xerr=None, yerr=None, color='black')\n",
    "\n",
    "\n",
    "#plt.xscale('log')\n",
    "\n",
    "plt.legend(loc='upper left')\n",
    "plt.ylim(0.7,2.0)\n",
    "plt.axhline(1, color = \"blue\", linestyle=\"dashed\", linewidth=1.5) \n",
    "\n",
    "\n",
    "fig = plt.figure()\n",
    "\n",
    "\n",
    "\n",
    "err_JS_300_400 = [err_1[2],err_2[2],err_3[2],err_4[2],err_5[2]]\n",
    "err_JS_400_500 = [err_1[3],err_2[3],err_3[3],err_4[3],err_5[3]]\n",
    "err_JS_500_1000 = [err_1[4],err_2[4],err_3[4],err_4[4],err_5[4]]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a61d37f",
   "metadata": {},
   "source": [
    "# for martini and ads load this file with reco respectively, and get the values of R_x00_y00 (3 arrays) and input that array in R_x00_y00_M and R_x00_y00_A to plot all three models simultaneously."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "f384161e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "R_300_400_A= [1.097410140890613, 1.1605204479420221, 1.1891569592333242, 1.2159964562826286, 1.2255819350077892]\n",
      "R_400_500_A= [1.0834254707157707, 1.1832065256294306, 1.2180126541687064, 1.280160791901948, 1.2604416948018737]\n",
      "R_500_1000_A= [1.0894023986419905, 1.1205858468179841, 1.2517788733138633, 1.4357689871926322, 1.485635932020225]\n",
      "err_JS_300_400_A= [0.06931080402058196, 0.07199129946631468, 0.07284418231214698, 0.07368295058305094, 0.07376735984984974]\n",
      "err_JS_400_500_A= [0.16368694118847074, 0.17431058714420328, 0.1765823604999301, 0.18235590293651, 0.17894588462428027]\n",
      "err_JS_500_1000_A= [0.38765515035748915, 0.3935150661543935, 0.42512527616212037, 0.47116679185456184, 0.48071687800133]\n"
     ]
    }
   ],
   "source": [
    "print(\"R_300_400_A=\",R_300_400)\n",
    "print(\"R_400_500_A=\",R_400_500)\n",
    "print(\"R_500_1000_A=\",R_500_1000)\n",
    "print(\"err_JS_300_400_A=\",err_JS_300_400)\n",
    "print(\"err_JS_400_500_A=\",err_JS_400_500)\n",
    "print(\"err_JS_500_1000_A=\",err_JS_500_1000)\n",
    "\n",
    "\n",
    "\n"
   ]
  }
 ],
 "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
}