cms_double_martini_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'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d70a0423",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Please set the path for pbpb jet data file(MARTINI 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/MARTINI/R0.2/cms_martini_0-10_R0.2_s1-6_12790.dat'\n",
    "file_pbpb1 ='../../data/music_cms/MARTINI/R0.3/cms_martini_0-10_R0.3_s1-6_12790.dat'\n",
    "file_pbpb2 ='../../data/music_cms/MARTINI/R0.4/cms_martini_0-10_R0.4_s1-6_12790.dat'\n",
    "file_pbpb3 ='../../data/music_cms/MARTINI/R0.6/cms_martini_0-10_R0.6_s1-6_12790.dat'\n",
    "file_pbpb4 ='../../data/music_cms/MARTINI/R0.8/cms_martini_0-10_R0.8_s1-6_12790.dat'\n",
    "file_pbpb5 ='../../data/music_cms/MARTINI/R1.0/cms_martini_0-10_R1.0_s1-6_12790.dat'\n",
    "\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.08555746 1.09252599 1.08742146 1.10790295 1.09855704]\n",
      "R=0.4/0.2= [1.13838423 1.14330314 1.13712018 1.14173331 1.14726646]\n",
      "R=0.6/0.2= [1.20679015 1.18408026 1.19239073 1.16678523 1.30287189]\n",
      "R=0.8/0.2= [1.25126776 1.20022524 1.20173558 1.2077353  1.35372505]\n",
      "R=1.0/0.2= [1.26828245 1.21144635 1.21224454 1.29235556 1.29698375]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgkAAAHBCAYAAADjMSi8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABKuElEQVR4nO3de3xcVbn/8c/TNtAm0HZ6IwTlUlopVqXaOXJRQSGUawU1pYiIihLP8XDUhCoeadrQ8ON3CqWtHtHDBC+IqNiK8Cu0WFIuR6Bwmh4FpS13pDStDXSg0DTQJs/vjz0Jabpzn1tmvu/Xa7+S2XvttZ+ZPck8s/Zaa5u7IyIiItLZkEwHICIiItlJSYKIiIiEUpIgIiIioZQkiIiISCglCSIiIhJKSYKIiIiEyrokwcxmmdn/M7MtZvaWma03s8/3Yr8DzewGM9tuZrvM7B4zOzINIYuIiOSkrEsSgErgLaAC+DTwAPBrM/u3Hvb7IfBlYA5QBowD7jOz4akLVUREJHdZtk2mZGbj3P3VTut+DZzo7kd1sc97gJeAS939l4l1hwEvAt9w95tTG7WIiEjuybqWhM4JQsKfgZJudpuR+HlHh3q2AA8DZyUvOhERkfyRdUlCF04Enulm+xTgFXd/q9P6jYltIiIi0kfDMh1AT8zsNOB84NJuikWA10PWxxPb9nPAAQf4nj172h8feuihlJR011gh2aCxsZHx48dnOoxBKR9fu1x5zoPleWRbnJmOJ13Hb2hoYOvWrW0P97j7AcmqO6uThMTohF8Dd7n7L5JZd0FBAe+8804yq5Q0iEaj1NfXZzqMQSkfX7tcec6D5XlkW5yZjicTxzezgmTWl7WXG8xsDLAK+DvwhR6Kx4FRIesjiW37GTdu3IDik8woLy/PdAiDVj6+drnynAfL88i2ODMdT6aPnwxZN7oBwMwKgTrgEIJRDdt7KH8pEANGufuuDuvXAK+7++c67xONRj2bMl4REZGBMjPc3ZJVX9a1JJjZMGAZMBk4s6cEIWF14udnOtRTAnyCoDVCREQkH2ztuUjvZWOfhB8DZwPfAsaa2dgO2/7s7m8nWghw99MSP18xs58CS83MgEagmuBSxa/SGbyIiEgGNSSzsmxMEtrmPPhByLajCCZNGhqy7ZvALmAxUAg8BHze3ZtTEKOIiEjOy7okwd2P7EWZT4ase5tgSufK5EclIiKSf7IuScgWO3fuZPv27XScS0HSr6CggAkTJjBy5MhMhyIikneUJITYuXMn//jHPzjssMMYMWIEQTcHSTd3Z/fu3WzZsgVAiYKISJpl3eiGbLB9+3YOO+wwCgsLlSBkkJlRWFjIYYcdxvbtvRnkIiIiyaQkIcSePXsYMWJEpsOQhBEjRuiyj4hIBihJ6IJaELKHzoWISGYoSUgHs2AREREZRJQkiIiISCglCSkWj8epBWqA2tpa4vHQ+02lRHV1NWbWvhQXF3Puuefy5JNPDrjuDRs2cNppp1FYWEhJSQnz5s2jpaWl1/u3trYSjUYxM+6+++4BxyMiIsmXt0lCY2Mj0WiUWCyWkvrdnZqqKiaWlLAGeBuoq6hgYkkJNVVVpOvGWqNGjWLt2rWsXbuWpUuX8swzz3D66aezY8eOftcZj8cpLS3FzLjrrruYN28eN9xwA/Pnz+91HTfffDOvvPJKv2MQEZF3xWIxotEoQHJvcezueblMnz7du7Jhw4Yut/XWgrlz/SOFhb4F3DssW8A/UljoC+bOHfAxejJ//nwfO3bsPuvWrl3rgN922239rvfaa6/10aNH+xtvvNG+buHChT5ixIh91nVlx44dPm7cOL/55psd8BUrVvS4TzLOiYhIrgPqPYmflXnbkpBK8XicxYsWsaKpiZJO20qAFU1NLF60iNdffz3tsR133HEAbN68ud91rFq1ijPOOGOfyY0uvPBCdu/ezUMPPdTj/lVVVXzsYx/jtNNO63cMIiKSekoSkikximH5mDHMaG7eL0FoUwLMaG5mWSSS9pEPL7/8MgBHHXVU+7q9e/f2uHiHyyObNm1iypQp+9R7+OGHU1hYyKZNm7o9/pNPPsnPfvYzFi1alMRnJSIiqaAkIQW2AZN7KDMpUS4d2j7on3/+eS6//HKmTZvGeeed1769oKCgx+WWW25pLx+Pxxk9evR+x4lEIj12zPy3f/s3Lr/8ciZNmpS05yciIqmhezckU+LbdnFtLXUVFbBrV5dFnysqonTJErjsspSG9Nprr1FQUND+eOzYsaxbt44DDzywfd26det6rKdjy0N//fa3v+Xpp59mxYoVA65LRERST0lCCpSVlfHdb36TBgi95LAFWN3Swk2zZqU8llGjRlFXV0dLSwtPPPEEc+bM4aKLLuKRRx5hyJCgIWnatGk91jN06ND23yORCG+88cZ+ZeLxOJFIJHT/PXv28J3vfIcrr7yS1tZWXn/9dXbu3AnArl27ePPNNzn44IP78QxFRCRVdLkhBSKRCJVz5jCzsJCGTtsagE8XFlI5Z05ok32yDRs2jGg0yvHHH095eTk33ngjjz32GMuWLWsv09fLDVOmTNmv78HmzZtpamrar69Cm127dvHKK69QWVlJJBIhEom0d6K88MIL+fCHP5yCZy8iIgOhloQUmbtgAQBTFy1iRnMzkwguMaxuaaGysrJ9e7pdfPHFLFy4kIULFzJ79myg75cbzjrrLK6//vp9vv3ffvvtjBgxglNOOSV0/4MOOogHHnhgn3Xbtm3j85//PNdeey2nnnpqf5+SiIikSjLHUw6mJdXzJLTZsWOHx8AXgMdiMY/H40mruydh8yS4u992220OeF1dXb/q3bFjhxcXF3tpaanfd999ftNNN3lRUZFfddVV+5Q7+uij/dJLL+2ynhdffFHzJIiIJBGaJ2FwiUQiXAZUAZdddllaLjH0ZPbs2UyePJnrrruuX/tHIhHWrFlDS0sLM2fOZP78+VRUVHD11VfvU27v3r19mqpZRESyi7mnZ3rgbBONRr2+vj5028aNGzn22GOTd7C2eRDy9LVOhqSfExGRHGRm6909mqz61CchHZQciIjIIKTLDSIiIhJKSYKIiIiEUpIgIiIiofI2SWhsbCQajRKLxTIdioiIyIDEYjGi0SjAuGTWm7cdF8ePH09XoxtEREQGk/LycsrLyzGzV5NZb962JIiIiEj3lCSkgdm7UyWIiIgMFkoSREREJJSShBxWXV2NmbUvxcXFnHvuuTz55JMDrnvDhg2cdtppFBYWUlJSwrx583o9BfMdd9zBP/3TPzFixAjGjh3LmWeeya5duwYck4iIJJeShBSLx+NALVBDbW1t4nH6jBo1irVr17J27VqWLl3KM888w+mnn86OHTv6XWc8Hqe0tBQz46677mLevHnccMMNzJ8/v8d9b775Zi666CLOOussVq1axc0338zkyZPZu3dvv+MREZEUSebdogbTkuq7QLa2tvrcuQt8+PDRDrMdrvKiogt8+PDRPnfuAm9tbR3wMXoSdhfItWvXOuC33XZbv+u99tprffTo0f7GG2+0r1u4cKGPGDFin3WdNTY2+kEHHeSxWKzPx9RdIEVEeobuAjk4zJt3DYsX30lz81PAb4Fr2LXrdpqbn2Lx4juZN++ajMR13HHHAbB58+Z+17Fq1SrOOOMMRo4c2b7uwgsvZPfu3Tz00ENd7ve73/0OgC996Uv9PraIiKSPkoQUiMfjLFq0mKamFUBJp60lNDWtYNGixbz++utpj+3ll18G4Kijjmpft3fv3h4X73CTqk2bNjFlypR96j388MMpLCxk06ZNXR778ccf55hjjuGnP/0p73nPeygoKOD444/n0UcfTfKzFBGRZFCSkERtQx3HjFlOc/MM9k8Q2pTQ3DyDSGRZWoZHtn3QP//881x++eVMmzaN8847r317QUFBj8stt9zSXj4ejzN69Oj9jhOJRLrtc7Ft2zaefvpprrnmGhYuXMiKFSsoKirizDPP5B//+EdSn7OIiAxc3s64mFrbgMk9lJmUKJdar732GgUFBe2Px44dy7p16zjwwAPb161bt67Hejq2PPSXu/PWW2+xbNkyzjzzTABOOukkjjjiCH70ox9RU1Mz4GOIiEjyKElIorYW+draYioq6uhuVF9R0XMsWVLKZZelNqZRo0ZRV1dHS0sLTzzxBHPmzOGiiy7ikUceYciQoCFp2rRpPdYzdOjQ9t8jkQhvvPHGfmXi8TiRSKTLOiKRCGbGJz/5yfZ1I0eOZPr06WzYsKH3T0pERNJClxtSoKysjJaW1UBDFyW20NKymlmzZqU8lmHDhhGNRjn++OMpLy/nxhtv5LHHHmPZsmXtZfp6uWHKlCn79T3YvHkzTU1N+/VV6OjYY49t7zHbkbu3JywiIpI99J85BSKRCHPmVFJYOJP9E4UGCgs/zZw5laHX9VPt4osvZurUqSxcuLB93bp163pcZs6c2V7+rLPO4o9//CNvvvlm+7rbb7+dESNGcMopp3R57HPPPReABx54oH3dG2+8wfr169tHXYiISBZJ5njKwbSkd56ECxy+nxXzJLi733bbbQ54XV1dv+rdsWOHFxcXe2lpqd93331+0003eVFRkV911VX7lDv66KP90ksv3Wfdeeed58XFxf6LX/zC7777bj/55JN93LhxvmPHjm6PqXkSRER6huZJSI7Gxkai0SixWCwl9ZsZNTVVNDS8AJQCw1mypJStW1+kpqYKy+Adn2bPns3kyZO57rrr+rV/JBJhzZo1tLS0MHPmTObPn09FRQVXX331PuX27t2731TNv/rVrzj//POprKykrKyMgoIC7r///m77MoiISPdisRjRaBRgXDLrNe90fThfRKNRr6+vD922ceNGjj322KQdqy0fyNOXOimSfU5ERHKRma1392iy6tPohjRQciAiIoNR3l5uEBERke4pSRAREZFQShJEREQklJIEERERCaUkQUREREIpSRAREZFQShJEREQkVNYlCWY2ycxuMrMnzazFzB7s5X5RM1ttZjsSS52ZHZ/icHvFrjbs6szNsCgiItIfWZckAFOBs4GngWd6s4OZvReoI5gc6ouJZRhwn5kdkaI4RUREclo2Jgkr3P297j4LeKqX+5wDHAx8xt3vcfd7gM8ARQQJR8bE43FYDzwEtbW1weM0qa6uxszal+LiYs4991yefPLJAde9YcMGTjvtNAoLCykpKWHevHn73achTH19PTNmzGDMmDGMGTOG0tJSHn/88QHHIyIiyZd1SYK7t/ZjtwJgL7Crw7q3Eusy0s7v7lRVV1FyeAm8EERS8eMKSg4voaq6inTdM2PUqFGsXbuWtWvXsnTpUp555hlOP/10duzY0e864/E4paWlmBl33XUX8+bN44YbbmD+/Pnd7rd582ZKS0vZu3cvt956K7feeit79+7l9NNP5+9//3u/4xERkdTIlXs3/B5YANxgZv8nsW4eEAeWZSKgeVfPY/HPF9Nc3gwjg3W72AU7YfHPFwNQU12T8jiGDRvGCSecAMAJJ5zAkUceyYknnsi9997LRRdd1K86/+u//ovdu3dzxx13MHLkSE4//XR27txJdXU13/3udxk5cmTofvfccw9vvvkmf/jDHxg1ahQAJ510EuPGjWPlypX8y7/8S/+epIiIpETWtST0h7s3AJ8CPgf8I7F8FjjD3RvD9mloaNinKb66ujpp8cTjcRbdsIimsqb2BKHdSGgqa2LRDYt4/fXXk3bM3jruuOOA4Ft9f61atYozzjhjn2TgwgsvZPfu3Tz00ENd7rdnzx6GDRtGUVFR+7qDDjqIYcOGpa1lRUQk13S8tAyUJLPunEgSzOxQghaD9cBZiWU9cI+ZHR62T0lJCe7eviQjSWgbxTDmS2NoPqJ5/wShzUhoPqKZyCWRtI98ePnllwE46qij2tft3bu3x6Xjh/imTZuYMmXKPvUefvjhFBYWsmnTpi6P/bnPfY7CwkKuuOIKtm/fzvbt26moqCASiTBr1qwkP1MRkfxQXV3d/lkGNCSz7pxIEoDvEPRLKHP3e939XoJWhRZgTtqjeQsY00OZMYlyadD2Qf/8889z+eWXM23aNM4777z27QUFBT0ut9xyS3v5eDzO6NGj9ztOJBLptmNmSUkJDzzwAL///e855JBDOOSQQ7jjjjv44x//yPjx45P6nEVEZOBypU/CFOApd9/TtsLd3zGzp4Cj0xWEzw++bdfW1lLx44qgD0IXit4sYsk3lnDZZZelNKbXXnuNgoKC9sdjx45l3bp1HHjgge3r1q1b12M9HVse+mvr1q3MmjWL6dOnc/PNNwNw4403cs455/Doo49y+OGhjT4iIpIhuZIk/B0428wOcPd3AMzsQOADwIp0B1NWVsY3K78JOwm/5LATWp5rSUsT+6hRo6irq6OlpYUnnniCOXPmcNFFF/HII48wZEjQkDRt2rQe6xk6dGj775FIhDfeeGO/MvF4nEgk0mUd119/PXv27GH58uXticupp57K5MmTWbRoET/84Q/7+OxERCSVsu5yg5kVmlmZmZUBhwHj2x6bWWGizHNm9tMOu91M0FnjD2Z2jpmdC9wJHArE0vwUiEQizLliDoXLC4NEoaOdULi8kDlXzAltsk+2YcOGEY1GOf744ykvL+fGG2/kscceY9mydwd99PVyw5QpU/bre7B582aampr266vQ0aZNm5g6deo+LRsHHHAAU6dO5fnnn0/isxYRkWTIxpaECew/bLHt8VHASwRxt3+1dff1ZnYmMB+4NbH6r8Dp7v5ESqPtwoL5CwBYdMOioBPjmOASQ8tzLVReUdm+Pd0uvvhiFi5cyMKFC5k9ezbQ98sNZ511Ftdffz1vvvkmBx98MAC33347I0aM4JRTTumyjiOOOIKVK1fyzjvvcMABBwDw9ttv87e//Y2ZM2cO5GmJiEgqdOzhn0/L9OnTvSsbNmzocltf7dixw5mJ8yk8Fot5PB5PWt09mT9/vo8dO3a/9bfddpsDXldX1696d+zY4cXFxV5aWur33Xef33TTTV5UVORXXXXVPuWOPvpov/TSS9sf19fX+7Bhw/zss8/2u+++21esWOFnnnmmDxs2zP/yl790e8xknhMRkVwF1HsSPyuz7nJDrolEIjAdOAUuu+yytFxi6Mns2bOZPHky1113Xb/2j0QirFmzhpaWFmbOnMn8+fOpqKjg6quv3qfc3r1795mqefr06dx77728+eabfPGLX+SSSy6hqamJ++67r33+BhERyR7meTqJTTQa9fr6+tBtGzdu5Nhjj03asdrmQWgb/SB9l+xzIiKSi8xsvbtHk1VfNvZJyDlKDkREZDDS5QYREREJpSRBREREQilJEBERkVBKEkRERCSUkgQREREJlbdJQmNjI9FolFgs7bM2i4iIJFUsFiMajQKMS2a9eTsEcvz48XQ1T0KyWeKnBkKKiEgqlJeXU15ejpm9msx687YlQURERLqnJCGHVVdXY2btS3FxMeeeey5PPvnkgOp97rnn+PrXv86HPvQhhg4dyic/+cle7/vGG2/wla98hUgkwqhRo/jCF77Aa6+9NqB4REQkNfL2ckO6xONxWL4ctm2jtriYsrKy4H4OaTJq1CjuvfdeAF566SXmzZvH6aefzsaNGxkzZky/6nzqqadYuXIlJ5xwAnv27OnTvhdccAHPPPMMN998M0OGDOHKK6/k/PPP509/+lO/YhERkdRRS0KKuDtVNTWUTJwIa9bA229TUVdHycSJVNXUkK57ZgwbNowTTjiBE044gQsvvJBf/vKXbN++vT1x6I+ZM2eyefNmli1bxtSpU3u939q1a1m9ejW33HILn/vc5/jMZz7Dr371Kx5++GHq6ur6HY+IiKSGWhJSZN4117D4zjtpfuopKCkBYBdAQwOLZ84EoKaqKu1xtd1tcfPmzf2uY8iQ/uWWq1at4pBDDuHkk09uX/fRj36Uo446ilWrVlFaWtrvmEREJPnUkpAC8XicRYsX07RiRXuC0K6khKYVK1i0eDGvv/562mN7+eWXATjqqKPa1+3du7fHJRktH5s2bWLKlCn7rT/22GPZtGnTgOsXEZHkUpKQRJZYxixfTvOMGfsnCG1KSmieMYPIsmXt+6RS2wf9888/z+WXX860adM477zz2rcXFBT0uNxyyy0DjiMejzN69Oj91kcikaDvhoiIZBVdbkiFbdtg8uTuy0yaFJRLsddee42CgoL2x2PHjmXdunUceOCB7evWrVvXYz0dWx5ERCQ/KElIorYG+driYirq6oI+CF0oeu45lpSWclmKYxo1ahR1dXW0tLTwxBNPMGfOHC666CIeeeSR9r4F06ZN67GeoUOHDjiWSCRCY2Pjfuvj8XhaR3yIiEjv6HJDCpSVldGyejU0NIQX2LKFltWrmTVrVspjGTZsGNFolOOPP57y8nJuvPFGHnvsMZYtW9ZeJl2XG6ZMmRLa96CrvgoiIpJZaklIgUgkwpzKShbPnLl/58WGBgo//WkqKytDr8+n2sUXX8zChQtZuHAhs2fPBtJ3ueGss86ipqaGhx9+mI9//OMA1NfX88ILL3DWWWcNuH4REUkuJQkpsmDuXAAWTZ0adGKcNImi556jZfVqKisr27enm5nx/e9/ny984QusWbOG0047re2mIL3W1NTEypUrAdiyZQs7d+5k+fLlAJx99tkUFhYCMGnSJE455RR++tOfAnDiiScyY8YMLrnkEhYtWtQ+mdLHP/5xDX8UEclCutyQImZGTVUVDS+8AKWlMHw4S0pL2frii9RUVWGW6jENXZs9ezaTJ0/muuuu69f+27dvZ9asWcyaNYvHHnuMDRs2tD/evn17e7m9e/fS0tKyz7633347p5xyCpdeeimXXHIJ06dP5w9/+MOAno+IiKSGpWvmv2wTjUa9q7tAbty4kWOPPTZpx9JdIAcu2edERCQXmdl6d+9b83A38rYlobGxkWg0SiwWS/mxHCUIIiKSOrFYrO3S8bhk1pu3fRLGjx9PVy0JIiIig0l5eTnl5eWY2avJrDdvWxJERESke0oSREREJJSShC7ka4fObKRzISKSGUoSQhQUFLB79+5MhyEJu3fv3uf+EyIikh5KEkJMmDCBLVu20NTUpG+xGeTuNDU1sWXLFiZMmJDpcERE8k7ejm7ozsiRIwFoaGhgz549GY4mvxUUFHDIIYe0nxMREUkfJQldGDlypD6YREQkr+lyg4iIiIRSkiAiIiKhlCSIiIhIKCUJIiIiEkpJgoiIiIRSkiAiIiKhlCSIiIhIqLxNEhobG4lGo8RisUyHIiIiHdjVhl1tmQ5jUInFYkSjUYBxyazX8nXa4Wg06vX19ZkOQ0REOmlLEHx+fn4+DYSZrXf3aLLqy9uWBBEREemekgQREREJpSRBREREQilJEBERkVBKEkRERCSUkgQREREJpSRBREREQilJEBERkVBZlySY2SQzu8nMnjSzFjN7sA/7ftbM1pnZbjN7zczuNbOiFIYrIiKSs7IuSQCmAmcDTwPP9HYnM/sa8GtgFXAW8DXgWWBYCmIUERHJedn4AbrC3e8CMLPl9GIeajMbBywB/s3dazts+kNqQhQREcl9WdeS4O6t/djtgsTPW5IZi4iISD7LuiShn44nuDzxVTN7xcz2mNnjZnZSpgMTEREZrHIlSSgGjgHmAlcCM4FdwL1mdkjYDg0NDZhZ+1JdXZ22YEVERJKlurq6/bMMKElm3bmSJBhwEPBVd7/N3e8FzgdagMvDdigpKcHd2xclCSIiMhhVV1e3f5YBDcmsO1eShDjgwINtK9x9J7AeeH+GYhIRERnUciVJ2EjQmmCd1hvQn46QIiIieS9XkoS7Ez8/1bbCzEYB04EnMhKRiIjIIJd18ySYWSHBZEoAhwEjzaws8XiluzeZ2XPAQ+7+VQB3rzezu4Cfmtn3gFeB7wJ7gBvT+wxERERyQ9YlCcAEYFmndW2PjwJeIoh7aKcyFwPXA4uBQuAR4FR3j6csUhERkRyWdUmCu7/E/n0LOpc5MmTdW8C/JBYRERms5numI5CEXOmTICIiIkmmJEFERERCKUkQERGRUEoSREREJJSSBBEREQmlJEFERERCKUkQERGRUEoSREREJFTeJgmNjY1Eo1FisVimQxERERmQWCxGNBoFGJfMei1x/+m8E41Gvb6+PtNhiIhIJ21T7ubnp9PAmNl6d48mq768bUkQERGR7ilJEBERkVBKEkRERCSUkgQREREJpSRBREREQilJEBERkVBKEkRERCSUkgQREREJpSRBREREQvWYJJjZZ8zsm2Z2TKf1l6cuLBEREcm0bpMEM/sP4FvAJOA+M/t2h82XpjAuERERybBhPWw/B/iwu+81s6uBZWZ2mLt/h3en1xYREZEc1NPlhiHuvhfA3V8DzgSONLOf9mJfERERGcR6+qDfamYfaXvg7u8AswluzvWBVAYmIiIimdVTkvBloKHjCndvdfevAZ9IVVDp0NjYSDQaJRaLZToUERGRAYnFYkSjUYBxyazX3Pt3x24z+4i7/28yg0mnaDTq9fX1mQ5DREQ6aevw1r9Pp/xmZuvdPZqs+gbSr6AmWUGIiIhI9ul3kuDu5yQzEBEREckuPc2TUGRmC8zsKTN7I7E8ZWY1ZnZwuoIUERGR9OupJeE24EDgM8AhieUzwAGJbSIiIpKjeppM6Vh3P7/TumeAK83s6dSEJCIiItmgp5aEt8zsjM4rzexMYFdqQhIREZFs0FNLwiXAf5nZzcAriXXvAV4CvpTCuERERCTDuk0S3P0p4BNmNoEgOTDgFXf/RzqCExERkczpqSUBAHffDmxPcSwiIiKSRXSTJhEREQmlJEFEUs7QveVFBqNeJwlmdkPHnyIiIpLb+tKS8KnEz1NTEYiIiIhkF11uEBERkVB5myQ0NjYSjUaJxWKZDkVERGRAYrEY0WgUYFwy6zX33t2x28z+190/YmZ/dvcP9+kgZua9PVCaRKNRr6+vz3QYInmhrdNiVv0TkKyl90v/mdl6d48mq76UtiSY2VQzuw54OZXHERERkeTrS5LQNplSt7MtmtkEM6swsz8D64HXgdP6F56IiIhkSq+TBHc/s+PPzszsQjNbCfwvcDhwGbDV3a9192eSEayIiIikT6+mZe6lXwMPAh9395cAzEyXlERERAapZPZJ+DDB5YU/mdn9ZvYVYGgS6xcREZE0SlqS4O5PuPt3gCOA/yDohxAxs+VmNitZxxEREZH06FeSYGZnmdnDYdvcvdXdV7v7xcChwArgnwcQo4iIiGRAt0mCmX3MzDaa2Vtm9nMzO8bM/gf4T+Anvaj/SGAq8L6BhyoiIiLp1FNLwg+AJUAU2AQ8DtwPHOvut4XtoCGQIiIiuaGnJKHA3WPuvsndFwI73f177r6nc8FkDYE0s0lmdpOZPWlmLWb2YB+eD2Y2xMzqzczN7Ny+7CsiIiLv6mkIZEunx692UzZZQyCnAmcDjwEF/dj/a8B7+rGfiIiIdNBTS8K0xLf5FjNr6fC4NfG4o2QNgVzh7u9191nAU33Z0cwiwP8BrurHcUVERKSDbpMEdx/i7kM7LEM6/uxUNilDIN29tT9PJKEGeARYM4A6REREhOTOuAi0f8ivBlabWREwCygHliX7WB2Z2YeAS4EPpfI4IiIi+SKld4F0913u/gt3Pz2Vx0n4T+BH7v5cbwo3NDRgZu1LdXV1aqMTERFJgerq6vbPMqAkmXX3KkkwszFm9mszazSzbWYWS1z/71jmo2Y2z8weTWaAvYzvQuAY4Jre7lNSUoK7ty9KEkREZDCqrq5u/ywDGpJZd29bEpYAFwCPEsygeBJwd2JOhB+a2TZgLTAPGEifgj4zswLgemAhMMTMRgMjE5uLzOzgdMYjIiKSK3rbJ+EM4HvuvgjAgjaNGLAOeC9wL8EQyFXu/loqAu1GEcGQx8WJpaPfAs8Dk9Ick4iIyKDX2yRhAtB+rwZ3dzOrAb4KzHP3Xjfzp8BbwKc6rSsGfgN8n2CGSBEREemjvoxu6DwvwpbEz/uSFAsAZlZIMJkSwGHASDMrSzxe6e5NZvYc8JC7f9Xd9xJM4tSxjiMTv/7V3R9PZnwiIiL5oi9JwlIzWw9sTCzPJ9a/k+SYJrD/cMm2x0cBLxHE3Z+JmkRERKSXepsk/AD4AEHnxQlAx+mWf5QY0fBn4C/ApoFMiJSY0tl6KHPkQOsQERGR7vUqSXD3irbfzWwC8MHE8qHEz28AIxJFmoCDkhumiIiIpFufZ1x09+0E0x63T32cGO0wmSBp+EDSohMRkbwSj8dh+XLYto3a4mLKysqIRCI97ygp0dvJlLaa2c1m9hkz26+VwAPPuPtyd69OepQiIpLT3J2q6ipKDi+Bn1RA3XwqflxByeElVFVXtU0UJGnW25aEbxGMOPgvYJSZ/Qm4B7jH3Z9NVXAiIpIf5l09j8U/X0xzeXP7dHi72AU7YfHPgylwaqprMhhhfrK+ZGeJywofJUgYziG4PfTzJBIGgmGJe1IQZ9JFo1Gvr6/PdBgieaGtF7G+C0qYeDxOyeEl+yQI+9gJw2PD2bp5K6NHj053eIOKma1392iy6utTnwQPMorHE8t8Myvm3YThjkSAdcDd7v6zZAUpItIfdrUGOQ0K64EjCE8QCNY3H9FM5JIITE9jXDKwu0C6+zZ3/5m7fw4YC3wW+Dvw3WQEJyIieeAtYEwPZcYkykla9elyQ78PYmaeZb1OjjjiCB8/fjzl5eWUl5dnOhyRnKbLDdKd2tpaKn5cwa7zd3VZpujOIpZ8YwmXXXZZGiMbPGKxGLFYjPXr1/+9p7mE+iKlSYKZTQW+BHze3d+bsgP1g/okiKSPkgTpjvokJE+y+yQM6HJDmMTtoyvM7M8EV5peB05L9nFERCQ3RCIR5lwxh8LlhbCz08adULi8kDlXzFGCkAF9nkypK2Z2IXAJwYRKy4DLgGXufm2yjiEiIrlpwfwFACy6YRHNk4bCyCaKdhbS8lwLlVdUtm+X9Era5QYzayW4G+OliXsnYGYvuPvEpBwgyXS5QSR9dLlBeisejzMmMeNirLiYWbNmqQWhDzI6BLIHHwYuBv5kZs8Ct6I7NYqISB9EIhFIdE5UF8XMS1qfBHd/wt2/QzDa9T8I+iFEzGy5mc1K1nFEREQkPfqVJJjZWWb2cNg2d29199XufjFwKLAC+OcBxCgiIiIZ0G2SYGYfM7ONZvaWmf3czI4xs/8B/hP4SS/qPxKYCrxv4KGKiIhIOvXUkvADYAkQBTYRTMd8P3Csu98WtoOGQIqIiOSGnjouFrh7LPH7JjP7V3f/XlhBDYEUEZGBisfjkBjdUFtcTFlZWdCZUTKip5aElk6PX+2m7K+B4cDH3b3C3evRiCcREekFd6eqqoaSkolQsQbmv01FRR0lJROpqqohy2b2zxs9JQnTzKylbenwuDXxuKMPE1xe+JOZ3W9mX0FDIEVEpBfmzbuGxYvvpLn5Edh1GviB7NpVSnPzIyxefCfz5l2T6RDzUtLv3WBmQ4BSgksPnwZWA7e7+7KkHmiANJmSSHp0nhxHzcfSWTwe59BDj+Lt1nIYWgunnQEfmgRPPAv3r4aWyzhwSIxt217SxEo9yPp7N4QMgbwb0G0WRfKMu1NVU0PJxImwZg28/TYVdXWUTJxIVU16mo/NtAyGZcyY5bzdeiRMXAPPPwV3/xauvQbuuT14PHENb7ceSSSyLOOxZvuSbMmccXE/7r4L+EViEZE8Mu+aa1h85500P/UUlJQAsAugoYHFM2cCUFNVlbkAJYu8CEOfh/ufbn+vtCspgftXwNHHwJ4XMxNeHuvV5QYzGwP8CDidoDPj/wOudPd4hzIfBc4EznT3k1ITbvIcccQRPn78eMrLyykvV0OHSDLF43FKJk7cJ0HYR0MDw6dOZeuLL6a0+TgFX6wkFS65BHa8AXff1XWZc8+DMaPgl79MX1yDkdnf3f3IZFXX25aEJcAFwD3AduAk4G4z+wwwN7FtPMFohseSFVwqjR8/HvVJkHyS1g/M5cthxozwBAGgpITmGTOILFvWPk+/5LH3vCdYuvOBY2Go+sL3QnejEPust30SzgC+5+7nuftlwAeBDcA64HLgf4EvAYe4+8eTGWAuMvQNR3Lctm0weXL3ZSZNCsqlkGsZFEvsqKMo2Lixq9MIQMHTTxM78siMx5rtS7L1NkmYALTfq8GDaxQ1wHuBee5+trv/yt1fS0GMIpIEaf2nX1xM0bPPdhtP0XPPESsuHlT/MCU1ysrKGPrgg9DQEF5gyxaGPvggs2bpXoHp1pfRDZ3nRdiS+HlfkmIRkRxRVlZGy+rV3f7Tb1m9Wv/0BQhuDz2nspLCmTP3f880NFD46U8zp7JSwx8zoC+jG5aa2XpgY2J5PrH+naRHlcM05ajkg7Z/+otnzqRpxYp9+yYk/ulX6p++dLBg7lwAFk2dSvOMGTBpEkXPPUfL6tVUVla2b5f06u3ohiXABwj6Ikzg3ZY8A9YCjwJ/Bv4CbHL31lQEm0zpnkzJ3Zl3zTUsWryY5jPOCP4Ann2WltWrmZP4A7BUDHIVyZB93vOd/unrPS9d6Tz51qxZs5RM9kGyJ1Pq84yLZjaBIFn4IMHNnD4IvB8YkSjS5O4HJSvAVEl3klBVU8PiO+8M/1Y1cyaV55+vMeOSk/RPX/qqLXVUv5K+y3iSEFpJ8HVgMkHS8AF3rx5wpSmWziQhW8aMi2SK/ulLX+j90n9ZOS2zB55x9+WDIUGA4E5UlqZlzPLlQXNrL8aMpzIOSR+95iKSC5J+7wYJkSVjxkVERPoib5OE6WjMuIiISHfyNklIp7KyMt5eubLbMeNvr1ypMeMiIpJVlCSkiXsLnHFq6EQhnHlasF1ERCSLJGV0w2BkJeZ8PU0HW08w9dShBfD4UCg9FaZ8EDb9Feruh+NbYOseOJrgOkiK+Pz8PNeZoN7Z+9LrIX2h90v/JXt0Q19mXJT+egsYC3xiD0T3wIaV8PhKOAj4V4IZJtYkyg1y+uMWEckdeZskHF5wOONXjKe8vJzy8vKUHqu2tpaKH1ewi11BQhDSWlD0ZhFLvrGEy3TbXBER6aNYLEYsFgMYl8x68/ZyQ9onUzq8hObyZhgZUmAnDI8NZ+vmrYN+MiW1JAT0OuxLr4f0hd4v/ZeVkylJ9yKRCHOumEPh8kLY2WnjTihcXsicK+YM+gRBRERyS95ebki3BfMXALDohkU0H9EMY4JLDC3PtVB5RWX7dhERkWyhyw1pFo/HGfOlMfAWxD4fy7mb3aiZMKDXYV96PaQv9H7pP41uGOQikUh7x0V1UhQRkWymPgkiIiISSi0JmVCdaESbn9kwki0ej8Py5bBtG7XFxZSVlQUtJ3lGr4OI5Aq1JKRZPB4HaoEaamtrE48HN3enqqaGkokTYc0aePttKurqKJk4kaqaGvKl34teh3DxeBxqa6Emd97zInnD3fNymT59uqdTa2urz527wIcPH+0w2+EqLyq6wIcPH+1z5y7w1tbWtMaTTHMXLPDCj3zE2bJl3xd5yxYv/MhHfO6CBRmOMD30OuyrtbXV586f68MPGu58uMg52bxoWpEPP2i4z50/d1C/5yW12v52pO+Aek/iZ6VGN6RJVVUNixffSVPTCqCkw5YGCgtnUll5PjU1VWmLJ1ni8TglEyfS/NRTUFKyf4GGBoZPncrWF1/MqVEcnel12F9VdRWLf76YprKmfScRS8wNUvmVSmqqazIWn2QvjW7ov2SPbsi6JMHMJgHfAU4EpgJ/cvdP9rDPPwHfAD5B8Am8Gfg1sNDdm8P2SfuMiyUTaW5+in0ThDYNDB8+la1bU/sBYj0X6bvaWqirg9tv77rM7NlQWgq5PJpjELwO6fxLz6dZRiX5lCT0Xz4MgZwKnA08BhT0cp/ZBPdQXAg8C3wIqEn8/FzYDuvXg6XkUzPMcmAG4QkCQAnNzTOIRJYBKfwAScVf3LZtMGlS92WOPjool8u2bYNJR3dfZuLEjL4OdnXa3vDBnU+PIDxBIFjffEQzkUsiuvOpSBbLxo6LK9z9ve4+C3iql/v8h7uf7O617v6gu/+QoDXis2Z2ROpC7a1twOQeykxKlEsdT8ESKy6mYMOGbo9bsHEjseLilBw/W5YlBx8MTzzZ7evAX//K0pEjMxZjWr0FjOmhzBhy4s6nIrks61oS3L21H/u8GrL6z4mfJcDfO2+cPh3S1SWhtraYioo6du3qukxR0XMsWVI66FrkS0tL2fPNb0JDQ/i1+C1b2HPvvcxYujQpx0vrt+G+WAv89/BuXwfuX8O3W+7h269/O93RAen9Vr3PnU+7oDufimS/rOuT0JGZLQfG9dQnoYt9vwksAQ519+2dt2dXn4QtDB/+gdT3SUjFB+x6YH0BjJ0I996/7wdkQwOceSq89gJM35PSZuWMewj4ewEUdfM67HoBjtgDp2QmxHQmCeqTIAOhPgn9lw99EgbMzIqBucCtYQkCQENDA9ahU8L8+fOprq5OSTyRSIQ5cypZvHhmF6MbPk1lZeXg/Gf5FjBxDxz4AhxzNJSeClM+CJv+CnX3w/EtMGpP0pqVs/Ua89KlS6m4sQImdPM6/H0PS89byre+9a1Mh5tybXc+7XZ0wxWD9D0vkmWqq6u5+uqr2x521fmtX3KuJcHMDgDqgPcA0909HlYu3UMg3Z15865h0aLFNDfPACZRVPQcLS2rmTOnkgUL5u6TtAwW7c3K5++C3cAGgoTgIOD9wAgoujN5zcrZ+xIthWFXwDdbg+62nV+HPcAPh8DexUBmkoR0/6m7O/Ounhfc+XTSUBjZRNHOQlqea2HOFXNYMH/BoHzPS+qpJaH/cn4IZEd9TRIs+I/zG+B04GPuvqmrshm9C+SY5cA2YrHiQX8XyHg8zoSSCez9xt4um5WH/XgYjVsbk/I8s/czpQaG/BEif4Yv7f/NmVsKIf5haD0DyMx8GJn6U4/H44xJTFMdKx7873lJPSUJ/afLDd1bCpwHnN5dgpBJwRz+wTfqXOmv5a1D4ZYC+NLukA/HEXhrn/uidn2sLP2vUVtbzLe/XUJT/FPww0Vw1FCY0ATbC+HFFmitpHD40yxdWpwz5723IpFI+5s9z566yKCXMy0JZvbvwDXABe7++57KZ6olAd79Npzulz4138JrgftgyDEwJOzDcQ60Pk3QuDPwj4hsfbvu2zl1BMHcGNuAYmAWsCstnVOzlb4ZSl/o/dJ/Od+SYGaFBJMpARwGjDSzssTjle7eZGbPAQ+5+1cT+1wEXAv8AthiZid0qPJ5d29MT/S9k60fdP2zDXgftNZAayU8uxye7fjhOBq4ilTPAZFp+3dO7ZgQDfLOqQOku2KKDF7ZOJnSBGBZYjmBoNtX2+MJiTLDgKEd9pmR+PllghHrHZdzUh7xIOGe/CUWK6ao6NnEEdoupVQlfo4GgjkgYrHipBwvmy1YMJfKyvMZPnwqFM0Gu4qiotkMHz6VysrzWbBgbqZDTCt3p6qqhpKSiVCxBua/TUVFHSUlE6mqyt+7YooMJll9uSGVMnm5IZdkyxwQ2UQd9QK5elMzST1dbui/vBrdkEpKEpJHHwb7y/d/ctlyUzMZnPL972cgcr5Pggw+bc3oixZNpbnlJNjzQYqKXqSlZTWVlZV518yerdI7fDQ7bmqWp9+BRJImG/skyCBjZtTUVNHQ8AKcuRI+tZAlS0rZuvVFamqqNGFOXsqOm5qJyMAoSZCkiUQSt/09BS677DI1I2eZVHRc7V2H1nDJ7NA6GDu6Srh4PA61tVBTQ21tbfBYMkZJQiaYZfPUgSIDVlZWRkvLaqChixJbaGlZzaxZs9IZlmQxd6emqoqJJSWcU1HBv8+fT11FBRNLSqipqtJomAxRnwRJmng8DtUxYBu1JbUaD5/HcvqmZpIS18ybx52LF/NUc/O775Zdu2gAZi5eDEBVTU2mwstbGt2QZvF4nOVjxgRz8cViOfFB6u5cM28eixct4mPNQ/kQTTxfVMjqlhYq58xh7oL8upFP5yGQuXCO+6O1tZVTTynloYfXMpTTaOEDDOVvtLCGUz5+Ivc/VMeQIWrMlOBvZmJJyb4JQgcNwNThw3lxq24t3pNkj27A3fNyOfzww3369Ol+0003eTq0trb6grlzffTw4T4b/CrwC4qKfPTw4b5g7lxvbW1NSxypsGDuXP9IYaFv6XRJeAv4RwoLfcHcuZkOMS06nuNzCvB/z6Fz3B9t74unwGPgCxI/N6TzfZHerhha+rnEwC/oocwFiXKZjjXbF+Al9+R9ViatosG2TJ8+vc//bwYiVz9Id+zY4aOHD9/veXV8fqOHD/d4PJ7pUFMuV89xf2TN+yIL/mlr6XlZQPDFqbsy30+Uy3Ss2b4A9e7J+6zU5YY0yJqmtBQ0+dcCdcDt3ZSZDZSSpNHwWfp+zZpz3J00XvJJ+/uiK1n6fpF91dbWUldRwe27dnVZZnZREaVLlnBZvt1GtY+Sfbkhfy8Irl//7iiDFC/Lx4xhRhcfHhB06ZrR3MyySCS1saRA2kfDp+mcDdpznObz3xXNkiB9UVZWxuqWlm7GwsDqlhaNhsmA/E0S0ihr/mGmoHGrOBbj2aKibg/7XFERxbEkXU3MUllzjruTxkbPtL8vBtn7RfYViUSonDOHmYWF+yUKDcCnCwupnDNHnRYzIH+ThOnT9Q8zCdL+DSCNH3SD8hxnyQemvhlKX81dsIDzKyuZOnw45xQV8T0zZhcVMXX4cM6vrGTuggWZDjEvqU9CGvR0vXoL8IFMX68egJqqKu5cvJgVTU2dRsPDzMJCzq+szPnxzbl+jvtD7wvpD91FdWA0BDJJi0Y3JE/HoX8XEPRCzsehf7l8jvuj4/vi7KIiv9IsL98X0ndt/6il79DohuRI92RK7u9OODSjuZlJBM3PuTThUOeJovLtG0A+nOP+0DdD6au2v5L8/HQamGS3JChJSLOc/yBt+xDM0/cV5ME57ge7Onhf+Pz8fV9I7ylJ6D8lCUmSqSQByO0P0lx+bn2h12EfShKkL5Qk9J/mSRAREZG00F0gM0HfLkVEZBBQS4KIiIiEUpIgIiIioXS5QZJLl1JERHJG3rYkNDY2Eo1GicVimQ5FRERkQGKxGNFoFGBcMuvVEEiRVNAQyH1oCKT0hYZA9p+GQIqIiEhaKEkQERGRUEoSREREJJSSBBEREQmlJEFERERCaZ4EkVTQqAYRyQFqSRAREZFQShJEREQklJIEERERCaUkQUREREIpSRAREZFQShJEREQklJIEERERCaUkQURERELlbZLQ2NhINBolFotlOhQREZEBicViRKNRgHHJrNc8T2eGi0ajXl9fn+kwRPKCWfAzT//dSB8l3i7o7dJ3Zrbe3aPJqi9vWxJERESke0oSREREJJSSBBEREQmlJEFERERCKUkQERGRUEoSREREJJSSBBEREQmlJEFERERCKUkQERGRUFmXJJjZJDO7ycyeNLMWM3uwl/uNMrOfm1nczN4ws9vMbGyKwxURkSRzM7xtmk7JqGGZDiDEVOBs4DGgoA/7/Q54H/A1oBVYCNwJfCLJ8YmIiOSFbEwSVrj7XQBmtpxe3KzCzE4EZgCnuPt/J9ZtAR43s1J3r0tlwCLSPdds/CKDUtZdbnD31n7sdhbwj7YEIVHP/wAvJraJiIhIH2VjS0J/TAE2hazfmNi2n6efhk9+ct91F1wA3/gGNDXB2Wfvv8+Xvxwsr74KZWX7b/+Xf4HZs2HzZvjiF/fffsUVMHNmcOyvf33/7XPnQmkp/OUv8O1v77/92mvhpJPg0Ufh+9/ff/vSpTBtGtTVwTXX7L/9ppvgmGNgxQq44Yb9t996K7z3vXD77fCTn+y/fflyGDcOfvGLYOls5UooLIQf/xh+97v9tz/4YPBz0SK4++59t40YAatWBb/X1MCaNftuHzsWfv/74Pd//3dYu3bf7e95D/zqV8Hv3/528Bp29L73QdtdwcvL4Zln9t0+bVrw+gFcfDG88sq+2088Ef7v/w1+/9zn4LXX9t1+2mlQVRX8ftZZsHv3vtvPPRfmzAl+7/y+gzx47/E+juEZvff03ttP+HvvgUTA+r83kPdeMuRKkhABXg9ZHwcmhu3wzjvv8NBDD7Y/PuKII4Ejkx6YiIhIKlVXV/ODH4wGpgGUJLNu8yy+wXtbnwR3/2QP5e4Ddrn7+Z3W/wqY6O4ndd4nGo16fX19EqMVkS619VTP4v83kkX0fuk3M1vv7tFk1Zd1fRL6KQ6MClkfSWwTERGRPsqVJGET4X0PuuqrICIiIj3IlSRhFVBsZh9vW2FmUYL+CKsyFpWIiMgglnUdF82skGAyJYDDgJFm1tandqW7N5nZc8BD7v5VAHdfa2argV+a2RzenUzpYc2RIJIFdG1ZZFDKuiQBmAAs67Su7fFRwEsEcQ/tVGY2sAT4GUELyd3AN1MWpYiISI7LuiTB3V8Cup20292PDFn3OvCVxCIiIiIDlCt9EkRERCTJlCSIiIhIKCUJIiIiEkpJgoiIiIRSkiAiIiKhlCSIiIhIKCUJIiIiEipvk4TGxkai0Sixthu9i4iIDFKxWIxoNAowLpn1ZvWtolNJt4oWEclSulV0v+lW0SIiIpIWShJEREQklJIEERERCaUkQUREREIpSRAREZFQShJEREQklJIEERERCaUkQUREREIpSRAREZFQShJEREQklJIEERERCaUkQUREREIpSRAREZFQShJEREQkVN4mCY2NjUSjUWKxWKZDERERGZBYLEY0GgUYl8x6zfP0ft3RaNTr6+szHYaIiHRmFvzM08+ngTCz9e4eTVZ9eduSICIiIt1TkiAiIiKhlCSIiIhIKCUJIiIiEmpYpgMQERHZhzosZg21JIiIiEgoJQkiIiISSkmCiIiIhFKSICIiIqGUJIiIiEgoJQkiIiISSkmCiIiIhFKSICIiIqGUJIiIiEiovE0SGhsbiUajxGKxTIciIiIyILFYjGg0CjAumfWa5+n0l9Fo1Ovr6zMdhoiISNKY2Xp3jyarvrxtSRAREZHuKUkQERGRUEoSREREJJSSBBEREQmlJEFERERCKUkQERGRUEoSREREJJSSBBEREQmlJEFERERCZWWSYGbvN7M1ZtZkZg1mtsDMhvZiv6iZrTazHYmlzsyOT0fMIiIiuSbrkgQziwB1gAPnAQuAK4Cre9jvvYn9hgFfTCzDgPvM7IhUxiwiIpKLhmU6gBD/DIwAPuvuOwk+5EcC1WZ2XWJdmHOAg4HPuPsbAGb2KPAqcDbwk9SHLiIiklElyaws61oSgLOAP3ZKBn5LkDic0s1+BcBeYFeHdW8l1lnnwo2NjQOPVNJOd+3sv3x87XLlOQ+W55FtcWY6ngwd/9BkVpaNScIUYFPHFe7+MtCU2NaV3yfK3GBmE8xsArAEiAPLOhd+9dVXkxawpE+m/+gHs3x87XLlOQ+W55FtcWY6nkwfPxmy8XJDBHg9ZH08sS2UuzeY2aeAu4FvJlZvBc5w9/2aDZqamt4xswM6rNoKNPQ3aEmbcWamDK9/8vG1y5XnPFieR7bFmel40nX8Et5tQdiTzIqzMUnoFzM7lKDFYD3wtcTqfwXuMbOTEq0R7dz9wDSHKCIiMqhkY5IQB0aFrI8ktnXlOwT9EsrcfQ+Amd0PPAvM4d3WBREREemFbOyTsIlOfQ8SwxsL6dRXoZMpwFNtCQKAu78DPAUcnYI4RUREclo2JgmrgDPM7OAO62YDu4GHutnv78AHOvYzMLMDgQ8AL6UgThERkZyWjUnCfwFvA3eYWamZlQPVwOKOwyLN7Dkz+2mH/W4m6LzxBzM7x8zOBe4k6MyR8i6mZvbexCyRG83sKTO7zsz2G3op2Ufnrv/02g1OOm/5y8weMrMnzOxJM1uemIeoS1mXJLh7HDgNGAqsIJhpcQkwv1PRYYkybfutB84kmFDpVuCXBJcoTnf3J1IfOXuBK939WODDwPHAZ9NwXBk4nbv+02s3OOm85a9Pu/tx7v4h4GWC/nxdyrokwcxmAf8BHAO0EAxL3ODuLR3LufuR7v7lTuvWuPvJ7j4msZzi7g92c6xJZnZTIqNqMbPuynZ7Pwl33+ru9Ynf3wGeBN7b1+c/mJlZmZk9amavmVmzmT1tZnM7DTVN1rFy9tyZ2WFm9paZuZkdlIL6c/a1ywQzG2Zm3zOzZ83sbTN7xcyWpOA4Om95qrfnvqfzDtBhRuIhQBHBLRC6lI2jGyqBF4EK3p1S+ddmNs7d/zPJx5qaqP8xgpERoezd+0lsILifxNHADQRJ1tyQ8mOB84EZSY43240F7geuJ5jr4qMEl4qKgcuTfKxcPnfXE8wWWpSi+nP5tcuEXwCnErR6biL4sH1/Co6j85a/ejz3fTnvZrYS+CeCjv1XdHtkd8+qBRgXsu7XwIvd7PNh4D1dbDsbGNrFtiEdfl8OPNhFuX8nGH45ssO67xLM8DiyU9kDgQeAKzL9WmbDAvwfgoTBdO569XqdDOwgGLbrwEHdlNVrl+GF4BLnHuD9fdhH501LX99nPZ77vpz3xLahwHXAd7s7dtZdbnD3sNmp/kz3N62oJrgR1LiOK83sEoJ+Dad2cazWXobVq/tJJJp1bgP+7O439LLuXPca0N3lhmp07jrG8J8Edz7tzSxt1ei1y7RLgfvdfUMf9qlG5036oJfnvk/3PfLgEv4twCXdVZp1SUIXTgSe6Wb7lwmy+Xvbemqa2XnAz4Dvu/t9Azx+b+8ncRPwJj013+Q4MxtqZoVm9nGCSax+4onUNcSX0blr888E3+xu7GX5L6PXLtOOB54xsx+Z2c7EteA7zKy7LzVfRudNkq/H825mETM7pEORzwF/667SrE8SzOw0gutlXWa6HoyIOINgVsb/Z2ZnA7cDN7j7wiSE0eP9JMzsY8BXgSjwZzP7i5nl6yyPuxLLnwjmtuiy96zOXSBxXbgGqPQOE4J1R69dVigm+NCfBlwIfAWYTjAUO3RIoc6bpEhv7nsUAe5OdID8K8E8Qt/qrtJs7LjYzsyOJOiPcJe7/6K7su6+1cxKgUeAe4Cb3f3KlAf57vEfIeSW1HnqJILhpx8F5gE/Ar7RVWGdOyDou/GYu6/sy0567TLOEst57v4agJltJUiOTwXWhO2k8yaZ4O4vEHRY7LWsTRLMbAzB7It/B77Qy91GEFz/bgUOMrMhfbiO153+3k8iL7n7/yZ+fdiCO6DdYmY3uPvz3eyWt+fOzKYSXNs+2cxGJ1YXJn6OMrMWd9/dTRV5+9plgTjwQluCkPAw8A7BCIfQJCFB502SKSXnPSsvN5hZIcEtnw8AznX3pl7scySwmmCIyCnAufT+2m5P+ns/CYG2hOGorgro3DGZYFjTWoI/5jjvPv9XCDozhtJrl3EbCf9GbgQf/qF03iQFUnLes64lwcyGEdzyeTJwkrtv78U+hwD3Ac8Ds9x9d6Ij0Eozi7v79wcY1irgO2Z2sLu/mVjXm/tJCHws8fPFsI06d0DwzfNTndadCVxJMCTuhbCd9NplhbuBqxPzuLSNSDmZIOkLnelV501SJDXnPZVjO/uzENxnwQl6xZ/QaTmwi33uBurZfwzweQS9iM/oYr9CoCyxrCWYWKLtcWGHchFgK8EfdilQTjDZzTWZfr2yaQHuJRjffxbBxCxXJ16n33azj85d+PP7Mj3Pk6DXLvPnaSTB1LZrgZnARcBm4D6dNy1JfJ/1eO5Tdd4z/uRDXoyXEv8cw5Yju9jnfYRMwpTYdhJdT+RzZG+PRXB98X6CrGwrQU/00ElP8nVJvCZ/S7wxXye41PBvQEE3++jchT+/L9NzkqDXLgsWYBKwkmBET5xgBsaIzpuWZC29PfepOO+WqFhERERkH1nZcVFEREQyT0mCiIiIhFKSICIiIqGUJIiIiEgoJQkiIiISSkmCiIiIhFKSICIiIqGUJIiIiEgoJQkiIiISSkmCSJ4zs7/0okx16iMRkWyjaZlFpEdm5u4edktkEclhakkQyXNm5h1+/4SZ/cnM1pvZf5vZ+81sSWLbX8zsT5mLVETSTS0JInmurZXAzMYQ3JP+THePm9nHgEXufqJaEkTyk1oSRKTNicBk4IFEP4UbgcP6UoGZPWhmnli+3Wnb+Wa22sxeM7N3zGyLmS03szP7eIwVZvbXbrb/yMxeN7MDzay6QzzL+3IcEVGSICLvGgI86u7TOiyH96OeBwgSjt+2rUhcsvg9sAX4GlAKfA8YAawys6P7UP9vgA+Y2fs7bzCzoUAZcIe7vw3cnIjlz/14HiJ5T0mCiLR5FJhuZh8GMLMhZvaRxLY3zWxkL+vZ4e6Pufu2RD3nAd8GvuruX3H3P7j7f7v7re5+DvBpYHcf4rwLaAI+H7LtU8AhBIkE7v6Kuz8G7OxD/SKSoCRBJIeZ2Q4zqzGz2xLN/DvM7CYzO7BzWXd/DZgF3GhmTwBPAZ9JbP4B8D/97Lj4bWCdu/8ibKO7r3D3hk5xf8LMHjKzpkTctWZ2cKL8LmAFMDukuguB7cD9/YhTRDpRkiCSo8zscCACVAA7CBKApUA58K+JMuMT2wBw94fd/SR3P87dj3X3qsT6Knef4u6f6GMMwwia+1f3YZ+PAXXANoJLB98GzgZ+3qHYb4DJZja9w34FwGeB37l7S1/iFJFwwzIdgIikzHGJnz9y9+8lfr/fzE4GzjSze4CVwLUpjGEscCCwueNKMzNgaIdVLf7uUKv/IOgbMbtD+S3AGjP7gLv/jWAUxusELQfrE8XOIEiKfpOC5yGSl9SSIJK7jiO4dn99p/XPAGPc/Wl3P9rdb0hDLJ3HWl8B7OmwtLVsFBK0PPzOzIa1LcDDiXLTAdz9HeAO4IJEwgHB5Ye/A2tT/FxE8oaSBJHcdRzw34m+Bh0VAw0h5VPhNeBt4D2d1t8K/FNi6ShC0MLwY/ZNIt4GCoD3dij7G+Bw4EQzGw6cB/y2Q4uEiAyQLjeI5K7jgH06GiY6/50KLOlwz4YxwEHAy4nHV7r7H5MRgLvvNbO1wAxgXof1/wD+kYip4y6vE7Q6VBNcCumsY3LzQKKOC4FDgYPRpQaRpFKSIJKDzKwIOJr9Wwy+mfj5A3e/OlH2y8An3f3LKQpnKXCnmX3R3W/trqC77zKzx4Bj3H1BD2VbzOx3BB0yDwM2uvsTyQpaRJQkiOSqDxJcTjzazOYRtCjMIOgLcLG7v56uQNz9LjNbCvzCzD5FMHzxVYJOjTMSxd7qsMt3CToptgLLgTcJLiucA1zl7s90KPsb4N8IhmrOT+XzEMlHShJEctNxBBMInUMwdPAq4Glglrvfle5g3L3CzP4b+AbwU4JLA40EnQzPdvdVHco+nBiBcTVB34WhBB0S7yVxiaJD2bVm9hJwJLrUIJJ0ShJEctNxwF/d/UkSIwLSyBIjEjoOa8Td/wD8oTcVuPvjQK/u6eDuR3UTyBCCFhXdnEqkHzS6QSQ3HQdk6vr8ZwlGJHwrQ8fvaB5BLCdnOhCRwUgtCSI5JjFvwAeBX2bg8F8nuJQA746WyKQYcHfi9x3dFRSR/SlJEMkxiSb+3t6MKdnHfjoTx+1K4p4Q6ZoTQiTnmOYdERERkTDqkyAiIiKhlCSIiIhIKCUJIiIiEkpJgoiIiIRSkiAiIiKhlCSIiIhIKCUJIiIiEkpJgoiIiIT6/24x60rPdxneAAAAAElFTkSuQmCC\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_ev_pp = 348300\n",
    "#n_ev_pp1 = 438421\n",
    "#n_ev_pp2 = 491662\n",
    "#n_ev_pp3 = 527578\n",
    "#n_ev_pp4 = 459913\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",
    "n_ev_pbpb = 12790\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": "iVBORw0KGgoAAAANSUhEUgAAAcwAAAGICAYAAADMEtD8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABAYklEQVR4nO3dfXyU1Zn/8c8hBJIQY1IiCUEICkJQSxFSFdcKGFs2UNdWkARbX7Zake3qtvVh1dZKLPxQqmJXKpa0VWtpCU9dFQVRQbCt2BVdbYtgpRJqQSKF8GTCU3L9/pjJNJNkknsmM5lJ8n3zul9hzn3uM9fJzOSa++kcZ2aIiIhI63rEOwAREZHOQAlTRETEAyVMERERD5QwRUREPFDCFBER8UAJU0RExIM2E6Zz7irn3LPOuV3OuSPOuTedc9M9bNfbOfeQc+5j59wnzrnnnXODoxK1iIhIB/Oyh3kLcAT4DvBvwCvAr51zN7ex3SPA14DbgKlANvCScy4l4mhFRETixLU1cIFzLtvM/tGk7NfAWDM7I8Q2pwOVwHVm9pS/bACwA/immf0sCrGLiIh0mDb3MJsmS7//A/Ja2ewL/p+/adTOLuB3QHE4AYqIiCSCSC/6GQv8pZX1BcDfzexIk/Kt/nUiIiKdSs9wN3DOFQFfAq5rpVoWcKCF8mr/OhERkU4lrITpv8r118AzZvZktINJS0uz2trawOP+/fuTl/fPI7979+7ltNNO89TW7t27g7YNJZw2E6Gu137FMgavdbtqv8KtG4u+ddV+xapuV+1XuHX1Xmy53d27d/PRRx81PKwxsz4tbmxmnhbgU/gOqf4vkNZG3R8CO1oofxTYEmq7MWPGWGvaWt+Yr2ttC6fNRKjrtV+xjMFr3a7ar3DrxqJvXbVfsarbVfsVbl29F9tuF9hsIXKUp3OYzrk04DmgF/BFM6tpY5NtwEDnXNMsXeBfF5EZM2ZEumlU2kyEuuGId7xdtV/h1g1HLH63nalfsarbVfsVbt1w6L3YnJfbSnoCzwDnAxeZ2fttNvrP20q+ZmaL/WV5/rKQt5UUFhba5s2bw4m/tRhoq2+dkfrV+XTVvqlfnU9X7Vs0++Wce9PMClta5+Uc5kJgEvAtoK9zrm+jdf9nZsecc+sAzKzI//PvzrmfAz9yzjlgL1AG7AQWR9yTMMyaNasjnqbDqV+dT1ftm/rV+XTVvnVUv7zsYVYC+SFWn2Fmlc65DQBmNr7Rdr2B+4BrgDRgI/AfZrYj1HNFcw9TREQkXO3awzSzwR7qjG+h7Bi+YfVuaTtEERGRxBb2fZjxdOjQIT7++GNOnDgR71BEOrXk5GT69etHRkZGvEMR6TQ6TcI8dOgQVVVVDBgwgNTUVHynRkUkXGZGbW0tu3btAlDSFPGo08yH+fHHHzNgwADS0tKULEXawTlHWloaAwYM4OOPP453OCKdRqdJmCdOnCA1NTXeYYh0GampqTq9IRKGTpMwAe1ZikSRPk8i4elUCVNERCRelDBFREQ8SKiEuXfvXgoLCykvL493KDG3cuVKLr30UjIzM+nduzfDhg3jlltuYffu3YDvcJlzjsWLmw+MtHjx4sD6xnbu3Mk111zDoEGDSElJYeDAgVxxxRW8+uqrnuMaPHhwoO2GJTc3t1m9d999l6KiItLS0sjLy+Oee+6hrq4uqI6ZMXfuXAYOHEhqaiqXXHIJb7/9tudYABYuXMjkyZPp27cvzjk2bNjQYr1oxuOlLRHpWsrLyyksLATIDlkp1Kjs8VhaG0X+3XffbXWE+c7klltusR49etj1119vzz77rG3YsMEWLlxoI0eOtC996Utm5ht9Pz093SZPntxs+0mTJll6enrQCP379++33NxcO++882zx4sW2fv16e/zxx+1LX/qS3X///Z5jy8/Pt6uvvto2bdoUWN58882gOvv377f+/ftbUVGRvfjii/bYY49ZWlqafe973wuqN3fuXEtJSbEFCxbYSy+9ZMXFxda3b1/76KOPPMdzwQUX2IUXXmjTp083wF555ZVmdaIZj9e2uoqu9LkSiQZama0k7kmy8RLrhJmTk2NAyCUnJ6fdz9GWZ5991gD7+c9/3mzdyZMnbfXq1WbmS5jTp0+35ORk279/f6DOvn37LDk52a6++uqghFleXm7OOauqqmrWbn19vef48vPz7dZbb221zty5cy0zM9MOHjwYKJs3b56lpqYGympray0jI8PuvffeQJ0jR45YdnZ2WMmnrq7OzMz+9Kc/hUyY0YzHS1tdiRKmSLDWEmZCHZKNtaqqqnatj4aHH36Y0aNHc9111zVbl5SURHFxceDx2LFjycvLY+XKlYGylStXkpeXx9ixY4O2PXDgAL169eJTn/pUs3ajfTXkmjVrmDhxYtAN76WlpdTW1rJx40YAXnvtNQ4dOsS0adMCdfr06cPll1/OmjVrPD9Xjx5tv0WjGY+XtkSke+pWCTPeTpw4wWuvvca//uu/eqrvnKOkpIQlS5YEypYsWUJpaWmzuqNHj+bYsWNcc801vPnmm9TX10cc589//nN69erFqaeeytSpU9m5c2fQ+m3btlFQUBBUNmjQINLS0ti2bVugTlJSEmeddVZQvREjRgTqREs04/HSloh0T0qYHWjfvn0cO3aMQYMGed6mtLSUjRs3UlVVxZ49e9i4cWOLCbOoqIjvfOc7LF26lMLCQjIzM5kyZQovv/xyWDFeccUVLFy4kHXr1vHAAw+wadMmPve5z3Hw4MFAnerqajIzM5ttm5WVRXV1daBOeno6SUlJzerU1NRw/PjxsOJqTTTj8dKWiHRPSphxEM4h0vPOO4+hQ4eybNkyli1bxrBhwxg1alSLdefPn89f/vIXHnjgAcaPH88LL7zAF77wBX7yk594fr7//u//Zvr06Xzuc59jxowZrF27lt27d/PEE094bkNEpCtSwuxAffv2pXfv3vztb38La7uSkhIqKiqoqKigpKSk1bpDhw7ltttu49lnn2Xnzp2MGjWK7373uxHPRn7uuecyfPhw3nrrrUBZVlZW0B5ng+rqarKysgJ1jhw50ux2jOrqatLS0ujVq1dE8bQkmvF4aUtEuiclzA6UnJzMv/zLv7B27dqwtistLWXTpk28/vrrLR6ODSU7O5uvf/3rVFdXt2uQ7ab3fBYUFDQ7n/fhhx9SU1MTOP9XUFBAXV0d27dvD6rX0jnC9opmPF7aEpHuSQmzg337299m8+bN/OIXv2i2rr6+nhdeeKFZ+YgRI5gxYwYzZswI+Ud77969LZa///779O7dm1NPPTWieP/85z+zbds2xowZEygrLi5m7dq1HD58OFC2dOlSUlNTGTduHAAXXXQRGRkZLF++PFCnpqaGVatWBV0JHA3RjMdLWyLSTYW63yQeS3e4D9PMN3BBUlKS3XDDDfbcc8/Zhg0b7Cc/+YmNGjUqaOCCBQsWhGxjwYIFQfdhPvDAAzZq1CibP3++rVu3zl544QW77bbbLCkpyW666SZPcT333HNWWloaGPhg4cKFlpeXZ2eccUbQPYgNgyRcdtll9tJLL9miRYusT58+LQ4UkJqaaj/+8Y/t5ZdftkmTJlnfvn1tz549nn9Xb7zxhi1fvtzmz59vgJWVldny5cvtjTfeiEk8XtvqKnQfpkgwNHBB4lmxYoWNHz/eMjIyLDk52c466yy79dZbA6POhJswt2zZYt/85jdtxIgRlp6ebqeeeqqNHj3aFi5caCdOnPAU0zvvvGOXXnqpZWdnW8+ePS0nJ8euvfZa27VrV7O6W7ZssQkTJlhKSorl5uba3XffbSdPngyqU19fb3PmzLEBAwZYSkqKXXzxxfbWW295iqXBtdde2+KXm2uvvTZm8Xhpq6voap8rkfZqLWE6i/BikFgoLCy0zZs3t7hu69atjBgxooMjEuna9LkSCeace9PMCltal1DnMLvT4OsiIpI4vAy+3rPjwmnbaaedRqg9TGm/kydPhlzXo0cPT8PQRVOixSMi3VfDhZXOuX+EqqO/SN1EZWUlycnJIZeWxrbtTvGIiLQlofYwJXby8vJ44403Qq7Pzg49BVwsJFo8IiJtUcLsJnr16tVwfD4hJFo8IiJt0SFZERERD5QwRUREPFDCFBER8UAJU0RExAMlTBEREQ+UMDvQk08+GZgqq/HSdIJnM2Pu3LkMHDiQ1NRULrnkEt5+++1m7b377rsUFRWRlpZGXl4e99xzT7P5Hluyfft2brzxRkaOHElSUhLjx49vsV404/DalhcNv7fFixc3W7d48eJm05FFouG1OnLkSJtxOOfo0aMHeXl5lJSUsGPHDsB3r2lLr3fTpbKysl2xikjH0G0lcbB+/XpSU1MDj88888yg9ffffz+zZ8/mgQceoKCggPnz53PZZZfx5z//mdzcXMA3ofFll13G2WefzTPPPMNf//pXbr31Vurr65kzZ06rz79lyxZWr17NhRdeyIkTJ0LWi2YcXtoKR3p6OhUVFXz1q18NKl+yZAnp6emtJrpouvXWW5k6dSpmxo4dO5g1axaTJ0/mj3/8I/3792fTpk2Buh988AFf+cpXePTRRxk9enSgvH///h0Sq4i0U6hR2eOxdPXZSp544gkD7PDhwyHr1NbWWkZGht17772BsiNHjlh2dnbQFFNz5861zMzMoGm35s2bZ6mpqUFlLamrqwv8f8qUKTZu3LiYxuG1La8Amz59uiUnJ9v+/fsD5fv27bPk5GS7+uqrg2Zyacm4ceNs1qxZIdd7ea1oYUaZ559/3gDbsmVLs/p/+tOfDLBXXnml1dg6Ulf4XIlEE63MVtKtDskOzh3c6qGxwbmD4x0ir732GocOHWLatGmBsj59+nD55ZezZs2aQNmaNWuYOHEiGRkZgbLS0lJqa2vZuHFjq8/hZYzWaMbhta1wjB07lry8PFauXBkoW7lyJXl5eYwdOzaiNqPhlFNOAWh1z11EOqeESpixnq1kZ9VOrJV/O6t2xuR5mxoyZAg9e/Zk+PDhLFq0KGjdtm3bSEpK4qyzzgoqHzFiBNu2bQuqV1BQEFRn0KBBpKWlBdWLVDTj8NpWOJxzlJSUsGTJkkDZkiVLKC0tjai9SNXX13Py5ElOnDjBX/7yF2bNmsVZZ53Fueee26FxiEj7eJmtJKESZsNsJTNmzIh3KDHRv39/Zs+ezS9/+UtWrVrFhRdeyMyZM3n44YcDdaqrq0lPTycpKSlo26ysLGpqajh+/HigXmZmZrPnyMrKorq6ut2xRjMOr22Fq7S0lI0bN1JVVcWePXvYuHFjyIRZV1fHyZMnA4uZBZJdw1JfXx92DN/61rdITk6mV69eDB8+nPfff5+VK1c266uIJLYZM2Y0zJYVcrYSXfTTgSZOnMjEiRMDj4uLizl69Chz5szhW9/6VtSns6qrq8MaTRCelJTU7qtHE8l5553H0KFDWbZsGWbGsGHDGDVqFL/73e+a1R0yZAg7dwYfQXj11VeZPXt24PGsWbMoKysLK4bbb789cKj5448/ZuHChUyaNInXX3+dAQMGhN8pEUlYSphxNnXqVJYtW0ZlZSVnnnkmWVlZHDlyhLq6uqC9lOrqatLS0ujVqxfg2zs7ePBgs/aqq6vJysoCoKioKOh85iuvvBLyFpKmohmH17YiUVJSQkVFBWZGSUlJyHqrVq3i2LFjgcc33ngjY8aMCTqakZeXF/bzDxo0KGgQ+aKiIk4//XQefvhhHnzwwbDbE5HEpYQZZw17fA0/CwoKqKurY/v27QwfPjxQr+m5woKCgmbn/z788ENqamoC9RYtWsThw4cD6xu315ZoxuG1rUiUlpYG9hIff/zxkPU+/elPBz0+5ZRTyMvLi/qMKb179+bMM89k69atUW1XROIvoc5hdkcrVqwgOzub/Px8AC666CIyMjJYvnx5oE5NTQ2rVq2iuLg4UFZcXMzatWuDEuLSpUtJTU1l3LhxgC9BFhYWBpaGKzi9iGYcXtuKxIgRIwIzpbc3+UbD0aNH+etf/8rAgQPjHYqIRJn2MDvQlClTOP/88xk5ciR1dXUsXbqUpUuX8sgjjwTOX6akpHDnnXcye/ZssrKyAjf519fXc/PNNwfamjlzJo888ghXXnkld9xxBx988AFlZWXccsstQbd4tKSmpobVq1cDsGvXLg4dOsSKFSsAmDRpEmlpaVGNw2tbkWo6UlI0Pf3006SkpASVffaznw18wamsrOT1118HfFd5P/rooxw8eJDrr78+ZjGJSJyEukEzHkusBy7Iz8k3IOSSn5Pf7udozV133WXDhg2z1NRUS0lJsdGjR9tTTz3VrF59fb3NmTPHBgwYYCkpKXbxxRfbW2+91azeli1bbMKECZaSkmK5ubl2991328mTJ9uMY8eOHSF/Bzt27IhJHF7b8oIWBgxobMGCBVEbuKCl5YknngjE0Xjp27evjR8/3tavX99imxq4QCTx0crABc4aXUUZinNuKHA7MBY4B/itmY33sF0hMBdoOFH0FvA9M/tDS/ULCwvNf1lvM1u3bmXEiBFtxioi3ulzJRLMOfemmbV4cYPXc5jnAJOA94C/eHzSgcDL+A77XuNfegIvOefyPT6viIhIQvB6DnOVmT0D4JxbQSsjITQyGTgF+LKZHfRv+xq+m0InAY+FH66IiEh8eNrDNLPwh0CBZOAk8EmjsiP+sq5z97yIiHQLsbytZCVQAzzknOvnnOsHPAxUA8tb3VJERCTBxCxhmtluYAIwBajyL1cCE81sb0vb7N69O2j2kHCHKRMREQlXWVlZ44nnQw75FbOE6Zzrj29P8k2g2L+8CTzvnBvU0jZ5eXlBl/AqYYqISKyVlZUF8g6wO1S9WA5ccDu+85hTzewEgHNuPfA+cBvwnzF8bhERkaiK5TnMAmBLQ7IEMLPjwBZgSAyfV0REJOpimTB3Auc65wJTUTjnegPnApUxfF4REZGo83RI1jmXhu/eSYABQIZzbqr/8Wozq3HObQc2mlnDIJo/A74B/I9zbiG+W0n+A+gPlEerAyIiIh3B6x5mP3wX8CwHLgTObvS4n79OTyAw2aGZvQn8K77BC34JPAWkAZ83s3eiEXxn8+STTwZdBdywNB083MyYO3cuAwcOJDU1lUsuuYS33367WXvvvvsuRUVFpKWlkZeXxz333ENdXV2bcWzfvp0bb7yRkSNHkpSUFHKOzGjG4bUtLxp+b4sXL262bvHixY2vdotYw2uVkZFBbW1ts/WXXXYZzjm+9rWvtet5Ghw/fpyysrJmv5PKykqcczz33HNhtbdhwwacc/z5z3+OSnwiQvcafD3eGgb0Xr9+vW3atCmwVFVVBdWbO3eupaSk2IIFC+yll16y4uJi69u3r3300UeBOvv377f+/ftbUVGRvfjii/bYY49ZWlqafe9732szjqefftpOP/10mzp1qhUUFNi4ceNarBfNOLy05RVg6enpNnny5GbrJk2aZOnp6W0Ovt6WhtcqPT3dli9fHrRuz549lpSUZOnp6Xbttde263kaHD58OGhg9wZHjx61TZs2WXV1dVjtHTx40DZt2mQ1NTWt1usKnyuRaKKVwdfjniQbL7GfrSSnjdlKctr9HK1p+CN8+PDhkHVqa2stIyPD7r333kDZkSNHLDs7OygJzZ071zIzM+3gwYOBsnnz5llqampQWUvq6uoC/58yZUqLCTOacXhtyyvApk+fbsnJybZ///5A+b59+yw5OdmuvvrqqM1WMn36dJsyZUrQukceecSGDx9uY8aMiUrCrKmpCZkwY00JUyRYawmzW00gvbOqKnS29K+Pt9dee41Dhw4xbdq0QFmfPn24/PLLWbNmTaBszZo1TJw4MWjuy9LSUmpra9m4cWOrz9Ew92ZHxeG1rXCMHTuWvLw8Vq5cGShbuXIleXl5jB07NqI2W1JaWsrq1auDJsiuqKigpKSkWd1t27ZRWlrKwIEDSUtL45xzzuFHP/oR9fX/HFmy4VDp2rVr+bd/+zfS09O56aabApN7f/3rXw8cUq6srGzxkOzgwYO57bbbePjhhzn99NPJysqitLSUAwcONHseHZIViZ5ulTATxZAhQ+jZsyfDhw9n0aJFQeu2bdtGUlISZ511VlD5iBEj2LZtW1C9goKCoDqDBg0iLS0tqF6kohmH17bC4ZyjpKSEJUuWBMqWLFlCaWlpRO2FMn78eDIzM3n66acB2LlzJ5s2bWL69OnN6u7atYvhw4ezcOFCVq9ezQ033MCsWbOYN29es7rXX389n/nMZ3j22We5/vrrWb9+PQB33303mzZtYtOmTfTv3z9kXMuWLWPdunWUl5czb948nnvuOb773e9Gp9Mi0qJYDlwgTfTv35/Zs2dz/vnnU1dXR0VFBTNnzqSmpobvfOc7AFRXV5Oenk5SUlLQtllZWdTU1HD8+HF69epFdXU1mZmZzZ4jKyuL6urqdscazTi8thWu0tJSHnroIaqqqjAzNm7cyPz58/nd737XrG5dXV3DKB6A71REfX09J0+eDJT16NGj2d53jx49uOqqq6ioqOCaa65h6dKljBw5stmXBICioiKKiooC7V988cXU1NTw05/+lLvuuiuo7lVXXcXs2bMDj48cOQL4vkxdeOGFbfY9OTmZp59+mp49fR/hd999l4qKChYuXNjmtiISGSXMDjRx4kQmTpwYeFxcXMzRo0eZM2cO3/rWtzwdKg1H0ySRlJTU7qtHE8l5553H0KFDWbZsGWbGsGHDGDVqVIsJc8iQIezcuTOo7NVXXw1KWrNmzWpxOMbS0lLGjRvH/v37qaioCLkXe/ToUe677z5+9atf8be//Y0TJwJjdnDy5MlAcgOYPHlyuN0NMmHChKD2zj77bD7++GNOnDhBcnJyu9oWkZYpYcbZ1KlTWbZsGZWVlZx55plkZWVx5MgR6urqgvbIqqurSUtLC+yJZWVlcfDgwWbtVVdXk5WVBfj2eBqfz3zllVdC3kLSVDTj8NpWJEpKSqioqMDMWjyv2GDVqlUcO3Ys8PjGG29kzJgxzJgxI1CWl9fymMsN50vnzp3L22+/zW9+85sW691xxx387Gc/Y9asWYwePZrMzEyeeeYZ5syZw9GjR0lPTw/UzcnJCberQZru1ffq1Qsz49ixY0qYIjGSUAlz7969FBYWMmPGjKA/ZF1Zwx5fw8+CggLq6urYvn07w4cPD9Rreq6woKCg2fm/Dz/8kJqamkC9RYsWBV2s0ri9tkQzDq9tRaK0tDSwl/j444+HrPfpT3866PEpp5xCXl4ehYWFnp6npKSEBx98kPPPP5/Bgwe3WGf58uXcfPPN/Nd//Veg7Pnnn2+xblfa0xfpCsrLyykvLwfIDlUnoS76Oe2009i8eXO3SZYAK1asIDs7m/z8fAAuuugiMjIyWL78n1OG1tTUsGrVKoqLiwNlxcXFrF27NighLl26lNTUVMaNGwf4EmRhYWFgabgS04toxuG1rUiMGDEi8AWrvcm3Nddeey2XX345t9xyS8g6tbW19O7dO/C44Ty1Fw172UePHm1foCISkRkzZrB582aAf4Sqk1B7mLGWn5ODa+XWkfx2HiZry5QpUzj//PMZOXIkdXV1LF26lKVLl/LII48Ezl+mpKRw5513Mnv2bLKysigoKGD+/PnU19dz8803B9qaOXMmjzzyCFdeeSV33HEHH3zwAWVlZdxyyy1Bt3i0pKamhtWrVwO+KzsPHTrEihUrAJg0aRJpaWlRjcNrW5FqOlJSLJx99tmBK2VD+fznP8+jjz7K0KFD+dSnPsWjjz4adBi4Nb169eKMM85g2bJlnHvuuaSkpDBy5MgoRC4i0dKtEmblnj1xff7hw4fz+OOP8+GHH2JmnH322Tz11FNcc801QfXuvPNO6uvrue+++9i3bx+FhYW89NJLQee9srKyWLduHTfddBOXX345mZmZfOc73/E0h+jHH3/MVVddFVTW8HjHjh2BQ47RjMNLW53dggULmDlzJv/xH/9Bamoq1157LV/+8pc9HzH5yU9+wm233cZll13GsWPH2LFjR4wjFpFwuMZXUcZbYWGh+XeJm9m6dSsjRozo4IhEujZ9rkSCOefeNLMWL25IqHOYIiIiiUoJU0RExAMlTBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDBFREQ8UMIUERHxIKESZsNsJf4R47ucJ598Eudcs6XpWKhmxty5cxk4cCCpqalccsklvP32283ae/fddykqKiItLY28vDzuuece6urq2oxj+/bt3HjjjYwcOZKkpKSQU35FMw6vbXnR8HtbvHhxs3WLFy8OrI+W8ePH45zj5ZdfbrausrIy6LXs06cPQ4YM4Stf+Qq//e1vm9U/ceIE8+fP59xzzyUtLY3s7GwuuOAC7r///qjFKyLhKy8vb5i9KORsJZhZwixjxoyxUN59992Q6zqLJ554wgBbv369bdq0KbBUVVUF1Zs7d66lpKTYggUL7KWXXrLi4mLr27evffTRR4E6+/fvt/79+1tRUZG9+OKL9thjj1laWpp973vfazOOp59+2k4//XSbOnWqFRQU2Lhx41qsF804vLTlFWDp6ek2efLkZusmTZpk6enp5ntrt9/f//5369GjhwF23XXXNVu/Y8cOA+zBBx+0TZs22YYNG+yJJ56wiRMnGmBlZWVB9WfMmGF9+vSxuXPn2rp162zlypV2++2320UXXRSVeMPVFT5XItEEbLYQOSruSbLx0l0S5uHDh0PWqa2ttYyMDLv33nsDZUeOHLHs7OygJDR37lzLzMy0gwcPBsrmzZtnqampQWUtqaurC/x/ypQpLSbMaMbhtS2vAJs+fbolJyfb/v37A+X79u2z5ORku/rqq9tMmOPGjbNZs2a1+VwPPfSQOefs0ksvtczMTDt27FjQ+oaEuWrVqmbbfv/73zfAXnnlFTMz++STTyw5Odl++MMfNqtbX1/fZiyx0BU+VyLR1FrCTKhDsrGWmzu4xUOiDUtu7uB4h8hrr73GoUOHmDZtWqCsT58+XH755axZsyZQtmbNGiZOnBg0lVdpaSm1tbVs3Lix1edomEqso+Lw2lY4xo4dS15eHitXrgyUrVy5kry8PMaOHRtRmy1ZsmQJF154IXfccQcHDhwIK95Zs2aRl5cXOOT+ySefcOLECXJzc5vV1YTSIomvWyXMqqqdgIVcfOtjb8iQIfTs2ZPhw4ezaNGioHXbtm0jKSmJs846K6h8xIgRbNu2Lahe0wmTBw0aRFpaWlC9SEUzDq9thcM5R0lJCUuWLAmULVmyhNLS0ojaa8n27dvZvHkzpaWlFBUV0a9fv6Dna0tSUhKXXnopr7/+OuCbIH3gwIGUlZXxm9/8JmjSbRFJfN0qYcZb//79mT17Nr/85S9ZtWoVF154ITNnzuThhx8O1KmuriY9PZ2kpKSgbbOysqipqeH48eOBepmZmc2eIysri+rq6nbHGs04vLYVrtLSUjZu3EhVVRV79uxh48aNIRNmXV0dJ0+eDCxmRn19fVBZfX190DZLliyhR48eTJs2jaSkJK666ipWrVrFJ5984jnG008/napGk5Y/+eSTHD58mClTppCZmUlhYSEPPvhgxL8DEek4SpgdaOLEidx999184QtfoLi4mF/84hdMmzaNOXPmNPtjHQ0tJYmu5LzzzmPo0KEsW7aMZcuWMWzYMEaNGtVi3SFDhpCcnBxYXn31VWbPnh1U9oMf/CBomyVLljBu3LjAIdTS0lJqamp45plnPMfY9Hd+6aWX8te//pWKigquu+469u3bx+23386ll14ak/eAiESPEmacTZ06lf3791NZWQn49rqOHDnS7LaM6upq0tLS6NWrV6DewYMHm7VXXV1NVlYWAEVFRUEJoa1zm41FMw6vbUWipKSEiooKKioqKCkpCVlv1apVvPHGG4Fl9OjR3HDDDUFlM2bMCNR/55132Lp1K1/84hc5cOAABw4c4JxzzqF///5hHZbdtWsXOTk5QWWnnHIKJSUl/PSnP+WDDz7g+9//Pr///e9ZtWpV+L8AEekwPeMdQHfXcLFHw8+CggLq6urYvn07w4cPD9Rreq6woKCg2fm/Dz/8kJqamkC9RYsWBZ0na9xeW6IZh9e2IlFaWsrs2bMBePzxx0PW+/SnPx30+JRTTiEvL6/hvqtmGpLirbfeyq233hq0bu3atezfv59PfepTrcZ28uRJ1q9fzyWXXBKyjnOO22+/ndmzZ7Nt2zauuOKKVtsUkfjRHmacrVixguzsbPLz8wG46KKLyMjIYPny5YE6NTU1rFq1iuLi4kBZcXExa9euDUqIS5cuJTU1lXHjxgG+BFlYWBhYTjnlFM9xRTMOr21FYsSIEcyYMYMZM2a0O/k2MDMqKiqYMGECr7zyStDy61//mhMnTgRdnRvKD37wA3bv3s3MmTMB36AFBw4caFbv/fffB2i2JyoiiUV7mB1oypQpnH/++YwcOZK6ujqWLl3K0qVLeeSRRwK3eqSkpHDnnXcye/ZssrKyKCgoYP78+dTX13PzzTcH2po5cyaPPPIIV155JXfccQcffPABZWVl3HLLLUG3eLSkpqaG1atXA75DhocOHWLFihUATJo0ibS0tKjG4bWtSDUdKam9Nm3axM6dO5k3b16LoyDdd999LFmyhBtuuCFQ9t5775Gdnc3x48fZsWMHFRUVvPDCC5SVlQW+OBw8eJBhw4Zx7bXXMmHCBE499VTee+897rvvPgYMGMCXv/zlqPZDRKIs1A2a8VhiPXBBTk5+6HtKwHJy8tv9HK256667bNiwYZaammopKSk2evRoe+qpp5rVq6+vtzlz5tiAAQMsJSXFLr74Ynvrrbea1duyZYtNmDDBUlJSLDc31+6++247efJkm3E03Gzf0rJjx46YxOG1LS8AW7BgQcj1CxYsaNfABTfddJNlZGRYTU1Ni+vnzZtnPXr0sN27dzf7XaakpNgZZ5xhV199tb366qtB2x07dszuu+8++9znPmf9+vWzlJQUGzJkiN1444324Ycftt7pGNHABSLBaGXgAmcJdOVkYWGhbd68ucV1W7duZcSIER0ckUjXps+VSDDn3Jtm1uLFDQl1DrOrD74uIiKJycvg6wl1DvO0004j1B6miIhIrDRcPOic+0eoOgm1hykiIpKolDBFREQ8UMIUERHxoFMlzES6oleks9PnSSQ8nSZhJicnU1tbG+8wRLqM2tpakpOT4x2GSKfRaRJmv3792LVrFzU1NfpmLNIOZkZNTQ27du2iX79+8Q5HpNPwdFuJc24ocDswFjgH+K2Zjfe47ZXAXcC5QA3wBjDFzLxPKgiBYdZ2797NiRMnwtlURJpITk4mJyenzWEUReSfvN6HeQ4wCXgd8HwMxzn3DeDHwA/xJdws4NIwnjdIRkaGPuAiIhIXXhPXKjN7BsA5t4JWRkJo4JzLBh4GbjaznzZa9T9hRykiIhJnns5hmlkkU8FP8//8RQTbioiIJJRYXvRzAfAecL1z7u/OuRPOuT845y6K4XOKiIjERCwTZi4wHLgbuAO4HPgEeME5p5lyRUSkU4llwnRAOnC9mf3KzF4AvgTUATe1tMHu3btxzgWWsrKyGIYnIiICZWVlgbwD5IWqF/Z8mA0X/bR1W4lzbilwFZBmZkcblb8MHDSzKU23aW0+TBERkViL13yYW/HtZbqm8QCRXEQkIiISN7FMmM/5f05oKHDOnQqMAd6J4fOKiIhEndeRftLwDVwAMADIcM5N9T9ebWY1zrntwEYzux7AzDY7554Bfu6cuxP4B/BfwAng0Wh2QkREJNa8DlzQD1jepKzh8RlApb+tpCZ1vgo8AMwH0oDfA5eaWXUkwYqIiMSLp4RpZpU0PxfZtM7gFsqOAP/uX0RERDqtTjNbiYiISDwpYYqIiHighCkiIuKBEqaIiIgHSpgiIiIeKGGKiIh4kFAJc+/evRQWFlJeXh7vUEREpBspLy+nsLAQIDtUnbAHX48lDb4uIiLxFK/B10VERLoMJUwREREPlDBFREQ8UMIUERHxQAlTRETEAyVMERERD5QwRUREPFDCFBER8UAJU0RExAMlTBEREQ+UMEVERDxQwhQREfEgoRKmZisREZF40GwlIiIiYdBsJSIiIu2khCkiIuKBEqaIiIgHSpgiIiIeKGGKiIh4oIQpIiLigRKmiIiIB0qYIiIiHihhioiIeKCEKSIi4oESpoiIiAcJlTA1+LqIiMSDBl8XEREJgwZfFxERaSclTBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPPCVM59xQ59wi59wfnXN1zrkN4TyJc66Hc26zc86cc1+MKFIREZE46umx3jnAJOB1IDmC5/kGcHoE24mIiCQEr4dkV5nZQDO7CtgSzhM457KA/wd8L9zgREREEoWnhGlm9e14jtnA74F17WhDREQkrrweko2Ic24kcB0wMpbPIyIiEmuxvkp2AfBjM9se4+cRERGJqZglTOdcKTAcmON1m927d+OcCyxlZWWxCk9ERASAsrKyQN4B8kLVC3u2EufcCiDbzMa3UicZ+ACYDzzhLx4EvAOUAqvN7HDT7TRbiYiIxFM8Zivpg+82kvlAtX95x7+uAvi/GD2viIhITMTqop8jwIQmZbnAEuC7wPoYPa+IiEhMeEqYzrk0fAMXAAwAMpxzU/2PV5tZjXNuO7DRzK43s5PAhiZtDPb/909m9od2Ry4iItKBvO5h9gOWNylreHwGUOlvKyk6YYmIiCQWTwnTzCoB10adwe1tQ0REJFFpthIREREPlDBFREQ8UMIUkW5tcO7goAFTmi6DcwfHO0RJEDEdS1ZEJNHtrNqJEXoAF1elSy/ER3uYIiIiHihhioiIeKCEKSIi4kFCJcy9e/dSWFhIeXl5vEMREZFupLy8nMLCQoDsUHXCnq0kljRbiYh0NOdc6xf94Eikv5MSW/GYrURERKRLUcIUkW4tPycf18q//Jz8eIcoCUL3YYpIt1a5pzLeIUgnoT1MERERD5QwRUREPFDCFBER8UAJUyTKNJi3SNeki35EokyDeYt0TdrDFBGRTiVeR3G0hykiIp1KvI7iaA9TRETEg4RKmBp8XURE4kGDr4vEgQbzFomtWH7GNPi6iIhIOylhikSZBvMW6Zp0laxIlGkwb5GuSXuYIiLSqcTrKI72MEVEpFOJ11Ec7WGKiIh4oIQpIiLigRKmiIiIB0qYIiIiHihhioiIeKCEKSIi4oESpoiIiAcJlTA1W4mIiMSDZisREREJg2YrERERaSclTBEREQ+UMEWkWxucm4tzLuQyODc33iFKgtDg6yLSre2sqqK1KzlcVVWHxSKJTXuYIiIiHihhioiIeOApYTrnhjrnFjnn/uicq3PObfCwzWedc08457Y752qcc+8552Y551LaHbWIiEgH83oO8xxgEvA6kOxxmxJgCDAPeB8YCcz2/5wSXpgiIiLx5TVhrjKzZwCccytoZSSERu43s380erzBOXcUWOScyzeznWHGKiIiEjeeDsmaWX24DTdJlg3+z/8zL9z2RERE4qmjL/oZC9QDf+3g5xURaVF+Tg4OQi75OTlxjE4SSYclTOdcLnA38Esz+7ilOrt37w66YbisrKyjwhORbqpyzx7MLORSuWdPvEOUGCsrKwvkHVo5Ahr24OsN5zDNbHwY2/QCXgZOB8aYWXVL9TT4uoiIxFNrg6/HfKQf50vZT+G70vZfQiVLERGRRNYRQ+P9CLgC+LyZbeuA5xMREYm6mJ7DdM7dBdwEfNXMfhfL5xJJFBrMW6Rr8rSH6ZxLwzdwAcAAIMM5N9X/eLWZ1TjntgMbzex6/zZXA3OBJ4FdzrkLGzX5VzPbG40OiCQaDeYtEluDc3PZ2crnKD8nJyYXa3k9JNsPWN6krOHxGUClv62kRuu/4P/5Nf/S2NfxJVIREZGwxOtLqaeEaWaV+G5Jaq3O4CaPv0bzRCkiItIpabYSERERD5QwRUREPFDCFJFuLTd3cKtXNefmDo53iJIglDAT3OA2PsyD9WEWaZeqqp2AhVx860U6ZuACaYedVTuxVq4Hc1WtXoslcZCfk9PqVXoazFukc1LCFIkyDdYtElvx+lKaUIdk9+7dS2FhIeXl5fEORUREElQsZpgpLy+nsLAQIDtUnbBnK4klzVbSnHOu9UOyOBLpNRTpbHzzQ7R6G7w+Y91Ia7OVJNQepoiISKJSwhQREfFACVNEurWcnHx8I3+2vPjWi+gq2YSXn5Pf6q0j+fowi7TLnj2V8Q5BOgklzARXqQ+ziEhC0CFZERERD5QwRaJMY5OKdE1KmCJRprFJRWIrXl9Ku0zC1CDlIiLdQ7y+lHaZi340SLmIiMRSl9nDFBERiSUlTBEREQ8SKmFqthIREYmHbjVbiWb1kESh2S9EYiuWnzHNViLSgTQ2qUjX1GWukhVJFBqbVCS2cnLyqWrlzodYfSntMglTg5SLiHQP8fpS2mUSpgYpFxGRWNI5TBEREQ+UMEVERDxQwhQREfFACVNERMQDJUwREREPlDBFREQ8UMIUERHxIKESpgZfFxGReOhWg693VYNzc9lZVRVyfX5ODpV79nRgRCIiXVdrg693mZF+uqqdVVWtj8nfSjIVEZHoSahDsiIiIolKCVNERDqV3NxcnHMhl9zc3Jg8rxKmiIh0KlVtnIpqa32klDBFREQ8UMIUERHxwFPCdM4Ndc4tcs790TlX55zb4HG7U51zTzjnqp1zB51zv3LO9W1XxCIiInHgdQ/zHGAS8B7wlzDaXwaMB74BfA34LPB0GNt3e/k5OTgIueTn5MQxOhGR7sPrfZirzOwZAOfcCloZCaGBc24s8AVgnJm96i/bBfzBOXeZmb0cYcwtys3NbfVEb05ODns64Q3+GpRARCQxeNrDNLP6CNouBqoakqW/nf8FdvjXRVW8rpoSEZHuIZYX/RQA21oo3+pfJyIiEracNk5FtbU+UrEcGi8LONBCeTVwZksbvPcejB8fXDZtGnzzm1BTA5MmNd/ma1/zLb2BY60E05vmbQPceitcfrnvuW+8sfn6u++Gyy6Dt9+Gb3+7+fq5c+Gii+C11+C7322+/kc/glGj4OWXYc6c5usXLYLhw2HVKnjooebrf/lLGDgQli6Fxx5rvn7FCsjOhief9C1NrV4NaWmwcCEsW9Z8/YYNvp8PPgjPPRe8LjUV1qzx/X/2bFi3Lnh9376wcqXv/3fdBZs2Ba8//XRYvNj3/29/2/c7bGzYMGgYZ3/GDPhLk7Pjo0b5fn8AX/0q/P3vwevHjoX77vP9f8oU2LcveH1REXz/+77/FxdDbW3w+i9+EW67zff/lt4bXt97//gHTJ3afP2//zuUlMCHH8I11zRfr/ee72e833tpabnU1oY+ApWSksMFFwSfGtF7L77vvT//eU+b7z2I7L3XmoQaS/b48eNs3Lgh8Dg/fzAw2NO2x6D1MVcjD0tEurDWkiXA0aM6ndPVVVZWsnNnZcPDvFD1wp6tpOGiHzMb30a9ZcBpZjahSfnzAGY2uek27ZmtxDnXZsJMpJlZRCQxONf212n97eg+WputJJbnMLfR8rnKUOc2RUREElYsE+YaINc5d3FDgXOuEN/5yzUxfN4uJTd3cBuDDA+Od4giIt2Cp3OYzrk0fAMXAAwAMpxzDaeaV5tZjXNuO7DRzK4HMLNNzrkXgaecc7cB9cA84HfRvgezK6uq2klrZ2erqnR2VkSkI3i96KcfsLxJWcPjM4BKf1tJTeqUAA8Dj+Pbm30O+M9IAhUREYknTwnTzCpp40JTMxvcQtkB4Ov+RUREpNPqMrOVaMxVEYlEvG6Cl84noe7DbA+NuSoikeiMY0xLfHSZPUwREZFYUsIUibLc3Nw2bgXKjXeIIhIBJcwEl5OTT+gzs86/XhKJZs4R6Zq6zDnMrmrPnsp4hyAiIiTYHubevXspLCykvGEaARERkQ5QXl5OYWEhQHaoOmEPvh5L7Rl8PTd3sH9UnJbl5OR3yr213NzcVg/h5eTk6Cq/BKPBvEU6r9YGX+8yh2S76hByOh8mIpIYEuqQrIiISKJSwhQREfFACVMkyjTUmkjX1GXOYYokCl2EJdI1aQ9TRETEAyVMERERD7pMwuyqQ8jpfJiISGLoMucwO+OgBF7ofJiISGLoMnuYIiIisaSEKSIi4oESpoiIiAcJlTA1W4mIiMRDt5qtREREpL1am60kofYwRUREEpUSpoiIiAdKmCIiIh4oYYqIiHighCkiIuKBEqaIiIgHSpgiIiIeKGGKiIh4oIQpIiLigRKmiIiIB0qYIiIiHihhioiIeJBQCVOzlYiISDxothIREZEwaLYSERGRdlLCFBER8UAJU0RExAMlTBEREQ+UMEVERDzwlDCdc2c759Y552qcc7udcz9wziV52K7QOfeic26/f3nZOXdB+8MWERHpWG0mTOdcFvAyYMAVwA+AW4F729huoH+7nsA1/qUn8JJzLr99YYuIiHSsnh7qzARSgSvN7BC+hJcBlDnnfugva8lk4BTgy2Z2EMA59xrwD2AS8Fi7oxcREekgXg7JFgNrmyTGCnxJdFwr2yUDJ4FPGpUd8Ze5MOMUERGJKy8JswDY1rjAzP4G1PjXhbLSX+ch51w/51w/4GGgGlgeWbgiIiLx4SVhZgEHWiiv9q9rkZntBiYAU4Aq/3IlMNHM9oYdqYiISBzF7LYS51x/fHuSb+I7rFvs///zzrlBLW2ze/dunHOBpaysLFbhiYiIAFBWVhbIO0BeqHptDr7unPsYeNTM7m1S/glQZmYPhNhuPr49yrPM7IS/rBfwPvCMmf1n0200+LqIiMRTewdf30aTc5X+W0bSaHJus4kCYEtDsgQws+PAFmCIh+cVERFJGF4S5hpgonPulEZlJUAtsLGV7XYC5/r3KgFwzvUGzgUqww9VREQkfrwkzJ8Ax4DfOOcuc87NAMqA+Y1vNXHObXfO/bzRdj/Ddyz4f5xzk51zXwSeBvoDmiFaREQ6lTYTpplVA0VAErAK3wg/DwOzmlTt6a/TsN2bwL/iG7zgl8BT+A7jft7M3olG8CIiIh3Fy0g/mNm7wKVt1BncQtk6YF1EkYmIiCQQzVYiIiLigRKmiIiIB0qYIiIiHihhioiIeKCEKSIi4oESpoiIiAdKmCIiIh4kVMLcu3cvhYWFlJdrICAREek45eXlFBYWAmSHqtPmbCUdSbOViIhIPLV3thIREZFuTwlTRETEAyVMERERD5QwRUREPFDCFBER8UAJU0RExAMlTBEREQ+UMEVERDzosgmzrKws3iHEhPrV+XTVvqlfnU9X7VuH9cvMEmYZM2aMtWbRokWtrm/M17W2hdNmItT12q9YxuC1blftV7h1Y9G3rtqvWNXtqv0Kt67ei223C2y2EDkq7kmy8dJWwmxrfZNOe6oXTpuJUDecN0a84+2q/Qq3biz61lX7Fau6XbVf4dbVe7HtdltLmAk1lqxz7jC+w8R7gX+0UCU7RHlL8oDdHuqF02Yi1PXar1jG4LVuV+1XuHVj0beu2q9Y1e2q/Qq3rt6LodvNBk4D6s3slJY2TKiEKSIikqi67EU/IiIi0aSEKSIi4kGnS5jOubOdc+ucczXOud3OuR8455La2OYc59wL/vrHnHN/c879zDnXv6Pibksk/WqyfQ/n3GbnnDnnvhjLWMMR4es12N+PpktFR8XtRXteM+fclc65N5xztc65ff73Z59Yx+xFhK9ZWYjXzJxzd3VU7K2J9PVyzhU65150zu33Ly875y7oiJi9akffzvH3rcY59w/n3GPOufSOiNkL59xQ59wi59wfnXN1zrkNHrc71Tn3hHOu2jl30Dn3K+dc3/bG07O9DXQk51wW8DLwLnAFMAR4CF/iv7uVTU8FdgBP4TsxfAYwCxjjnPusmZ2MZdxtaUe/GvsGcHpMAoxQFPp1G/D7Ro+9XgAQc+3pm3PuG8CPgR8CtwNZwKUkwOexHf36GfBCk7IvAXcAa6IeaJgi7ZdzbqB/u7eAa/zFtwMvOec+bWY7Yxm3F+3o26nAeuAvQAnQF997sj++1y4RnANMAl4HksPYbhkwDN/fxXpgHvA08Ll2RRPq8tlEXIC7gGogo1HZfwE1jcs8tvV5wIDRnb1f+P7g7gWu9/fpi/HuU3v6BQxOpH5EuW/ZwGHghnj3IZr9CtHW88DWePepna/XTKAOOLVRWZa/7N/j3a929u0u4BCQ2ajscv9nrzDe/fLH06PR/1cAGzxsM9bfh0salZ3vL7usPfF0tkOyxcBaMzvUqKwCSAXGhdnWPv/PXtEIrJ3a26/Z+PbE1sUgtvaI5uuVaCLt2zT/z1/EKrB2ispr5j/89XlgSXTDi1ik/UoGTgKfNCo74i9z0Q4yQpH2bRS+ew4PNCp7CV9imRzlGCNiZvURbFYMVJnZq43a+V98RxmL2xNPZ0uYBcC2xgVm9jd836QK2trYf56vl3NuOHA/8Abwv7EINEwR98s5NxK4Dt/hy0TTrtcLeMJ/3uIj59x851xqLIKMUKR9uwB4D7jeOfd359wJ59wfnHMXxS7UsLT3NWswBV+ySZSEGWm/VvrrPOSc6+ec6wc8jG+PbnmMYg1XpH1LAY43KTuJ7xDmiGgG2MGa/T78thLee7iZzpYws4ADLZRX+9e1ZTVwDN8v81P4DvlF8g0m2trTrwXAj81se7SDioJI+3UMeBTfIeYiYBHw7/i+NSeKSPuWCwzHd27pDnyHwD4BXnDO5UQ5xki09zPWoBR4y8zej0ZQURBRv8xsNzAB3xeAKv9yJTDRzPZGP8yIRPqabQc+45xrfG5wDJCE7+9jZxWt93AznS1httfNwIX4Tt6nA2uccynxDSlyzrlSfH9858Q7lmgys4/M7CYze9bMNphZGXAL8G/Ouc/EObz2cvjee9eb2a/M7AV8F1jUATfFM7Bo8V99Po7E2buMmL8vy4E38R3OK/b//3nn3KB4xhYFP8U3ss0C51yuc+4cYCG+92Ii7EgknM6WMKvxXfHaVJZ/XavM7H0z+4OZLQYmAucBV0c3xIiE3S//t8IH8F391cM5lwlk+Ff3cc61OLRTB2vX69XECv/PMe2KKHoi7Vs1vnNEGxoK/Oee3gTOjmJ8kYrGazYN3xeDpdEKKgoi7dft+A4tTzWzF/xfcKbgSyqJchokor6Z2TZgBjAd+Aj4I75TVG8De6IeZceJ5t+dIJ0tYW6jyTFo/2XfabR8zDok810Ovh84M2rRRS6SfvXBdxvJfHxvgmrgHf+6CuD/YhJpeKL2euFLMo1/xlukfduKL5k0vWDEkRjf6qPxmpUCvzOzD6McW3tE2q8CYIuZnWgoMLPjwBZ8t28kgohfMzN7HMgBRuIbj/UmYCi+2zg6q2a/D79Q5zY962wJcw0wscneUwlQC2wMpyH/hT998V05FW+R9OsIvnMrjZfp/nXfBb4Sm1DDErXXC5jq//lmNAKLgkj79pz/54SGAv/9cGP45xeeeGrXa+acG4zvtEeiHY6NtF87gXOdc4Gr6Z1zvYFzgcoYxBmJdr1mZnbUzP5kZlXAV/HlhWUxibRjrAFynXMXNxQ45wrx7Ry1757geN9nE+Y9OVn4Dh28BFyG73DCEWBOk3rbgZ83evwgvqtiv4zvD9U38b3ZtwN9Omu/WmhnMAl0/2I7Xq8yfDdeX+nf7gf4Pvwr492naLxm+G6g/gi4Ft/l+xvx3Ueb1Zn75S+/EzgBZMe7L1F6L47x9+d5/2v1RXx/dE8An4l3v9rZtwx8p3Qm4ztFdb+/X1+Ld58axZiG78vyVGATvj37hsdpod6LwFrgA//fkC/huzL9t+2OJ96/kAh+gWfjG52i1v8mmQ0kNalTCTzZ6HEpvvsU9+O71Hqb/w9ywnyoI+lXC20MJoESZjtfr83AQXyXvW/HlzR7x7s/0XjN8F308xi+e4Fr8Y3S8ul49yca70V8579eiHcfovx6FQGv+v9+7Mf3BWd8vPvT3r7hO63zor9Ptfhus/tSvPvSJOaGv2ktLYNbec0ygSfwXS17CPh1NP7ea3ovERERDzrbOUwREZG4UMIUERHxQAlTRETEAyVMERERD5QwRUREPFDCFBER8UAJU0RExAMlTBEREQ+UMEVERDz4/7N0FCdHhxUAAAAAAElFTkSuQmCC\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_M= [1.087421457060924, 1.1371201834649416, 1.1923907263828304, 1.2017355836199048, 1.2122445381874856]\n",
      "R_400_500_M= [1.107902948996912, 1.1417333072053066, 1.1667852266008045, 1.2077353007835328, 1.2923555587444864]\n",
      "R_500_1000_M= [1.0985570406473855, 1.1472664622184126, 1.3028718885511639, 1.3537250450673393, 1.296983750176387]\n",
      "err_JS_300_400_M= [0.05722454528682157, 0.058917362987257485, 0.06069459864960222, 0.060693234914417186, 0.06080044276875342]\n",
      "err_JS_400_500_M= [0.15511130554633093, 0.15791887747022149, 0.15901250213085974, 0.16216272340527998, 0.17006620500934602]\n",
      "err_JS_500_1000_M= [0.25694207865218277, 0.26393293540568064, 0.2890884673570552, 0.295898961361187, 0.2834177949359136]\n"
     ]
    }
   ],
   "source": [
    "print(\"R_300_400_M=\",R_300_400)\n",
    "print(\"R_400_500_M=\",R_400_500)\n",
    "print(\"R_500_1000_M=\",R_500_1000)\n",
    "print(\"err_JS_300_400_M=\",err_JS_300_400)\n",
    "print(\"err_JS_400_500_M=\",err_JS_400_500)\n",
    "print(\"err_JS_500_1000_M=\",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
}