Analysis-Ricci-IL150_MET150-new.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import csv\n",
    "import math\n",
    "from mt2 import mt2\n",
    "import sys\n",
    "sys.path.append(\"MLAnalysis\")\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from MLAnalysis.Interfaces.LHCOlympics import LoadLHCOlympics\n",
    "from MLAnalysis.CutAndExport.CutEvent import CutEvents\n",
    "from MLAnalysis.CutAndExport.Histogram import *\n",
    "\n",
    "from MLAnalysis.DataStructure.EventSet import *\n",
    "from MLAnalysis.DataStructure.Particles import *\n",
    "\n",
    "def mt(pt, etmiss, phi, etmiss_phi):\n",
    "    # Calculate azimuthal angle difference\n",
    "    delta_phi = np.abs(phi - etmiss_phi)\n",
    "    if delta_phi > np.pi:\n",
    "        delta_phi = 2 * np.pi - delta_phi\n",
    "\n",
    "    # Calculate transverse mass\n",
    "    mt = np.sqrt(2 * pt * etmiss * (1 - np.cos(delta_phi)))\n",
    "    \n",
    "    return mt\n",
    "\n",
    "\n",
    "class AtLeastThreeLeps:\n",
    "\n",
    "    def __init__(self,pt_max: float=50.0,met_max: float=50.0):\n",
    "        self.pt_max = pt_max\n",
    "        self.met_max = met_max\n",
    "\n",
    "    def Cut(self, eventSample: EventSample) -> bool:\n",
    "        lepCount = 0\n",
    "        bCount=0\n",
    "        lepPT=[]\n",
    "        bPT=[]\n",
    "        p_met=LorentzVector(0, 0, 0, 0)\n",
    "        p_lep=[]\n",
    "        for p in eventSample.particles:\n",
    "            if (ParticleType.Electron == p.particleType or ParticleType.Muon == p.particleType):\n",
    "                lepCount = lepCount + 1\n",
    "                lepPT.append(p.momentum.Pt())\n",
    "                p_lep.append(p.momentum)\n",
    "\n",
    "            if (ParticleType.Jet == p.particleType):\n",
    "                if(p.bTag):\n",
    "                    # print('found b')\n",
    "                    bCount = bCount + 1\n",
    "                    bPT.append(p.momentum.Pt())\n",
    "                \n",
    "            if(ParticleType.Missing == p.particleType ):\n",
    "                p_met=p.momentum\n",
    "\n",
    "        MT=0.0\n",
    "        pt_mt2=0.0\n",
    "        lepPT2=np.sort(lepPT)[::-1]\n",
    "        phi_met_ptl1=0\n",
    "        if(len(lepPT)>2):\n",
    "            p_lep.sort(key=lambda item: item.Pt(), reverse=True)\n",
    "            phi_met_ptl1=p_met.Azimuth()-p_lep[0].Azimuth()\n",
    "\n",
    "            phi_met_ptl2=p_met.Azimuth()-p_lep[1].Azimuth()\n",
    "            phi_met_ptl3=p_met.Azimuth()-p_lep[2].Azimuth()\n",
    "\n",
    "        #     idx= 2 #np.argmax(lepPT)\n",
    "        #     pl= p_lep[idx].values \n",
    "        #     phi_l=math.atan2(pl[2],pl[1])\n",
    "\n",
    "        #     phi_met=p_met.Azimuth()\n",
    "        #     pt_met= p_met.Pt()\n",
    "\n",
    "        #     MT=math.sqrt(2*lepPT[idx]* pt_met * (1.0-math.cos(phi_l-phi_met))) \n",
    "\n",
    "\n",
    "        # if(len(lepPT)>2):\n",
    "        #     idx= 2 #np.argmax(lepPT)\n",
    "        #     pl= p_lep[idx].values \n",
    "        #     phi_l=math.atan2(pl[2],pl[1])\n",
    "\n",
    "        #     phi_met=p_met.Azimuth()\n",
    "        #     pt_met= p_met.Pt()\n",
    "\n",
    "        #     MT=math.sqrt(2*lepPT[idx]* pt_met * (1.0-math.cos(phi_l-phi_met))) \n",
    "\n",
    "\n",
    "        # if(len(lepPT)>2):\n",
    "        #     # print(lepPT)\n",
    "        #     p_lep.sort(key=lambda item: item.Pt(), reverse=True)\n",
    "        #     pl1=  p_lep[0].values\n",
    "        #     pl2=  p_lep[1].values\n",
    "        #     # print(pl1,pl2)\n",
    "        #     pt_mt2= mt2(pl1[0],pl1[1],pl1[2],\n",
    "        #         pl2[0],pl2[1],pl2[2],\n",
    "        #         p_met.values[1], p_met.values[2],\n",
    "        #         100.0,100.0)\n",
    "\n",
    "        # return (lepCount < 3) or   lepPT2[2]<20 or lepPT2[0]<50 or lepPT2[1]<30 or  eventSample.GetETMissing() < 60\n",
    "        return (lepCount < 3) or (bCount >0) or\\\n",
    "               eventSample.GetPTLepton()< self.pt_max or    eventSample.GetETMissing() < self.met_max\n",
    "\n",
    "\n",
    "def get_pcd(ev,pt_max,met_max):\n",
    "    events = LoadLHCOlympics(ev)\n",
    "    # print(events)\n",
    "    cut = AtLeastThreeLeps(pt_max,met_max)\n",
    "    CutEvents(events, cut)\n",
    "    pcd_list=[]\n",
    " \n",
    "    for idx,event in enumerate(events.events):\n",
    "        ps=event.particles\n",
    "        met=[]\n",
    "        met_et=event.GetETMissing()\n",
    "        \n",
    "        for p in ps:\n",
    "            if ParticleStatus.Invisible == p.status:\n",
    "                met=p.momentum\n",
    "                phi_met=met.Azimuth()\n",
    "\n",
    "        for p in ps:\n",
    "            \n",
    "            if (ParticleType.Electron == p.particleType or ParticleType.Muon == p.particleType):\n",
    "                mom=p.momentum\n",
    "                eta=mom.PseudoRapidity()\n",
    "                phi=mom.Azimuth()\n",
    "                pt=mom.Pt()\n",
    "                R=mom.R()\n",
    "                mt=np.sqrt(2*pt*met_et* (1.0-np.cos(phi_met -phi)))\n",
    "                \n",
    "                pcd_list.append([idx,eta,phi,R,pt,mt] + mom.values)\n",
    "            \n",
    "    return pd.DataFrame(pcd_list,columns=['event_no','eta','phi','R','pt','mt','E','px','py','pz'])\n",
    "\n",
    "evs=[\n",
    "\"Events/nx-3lep-BP1.lhco\",\n",
    "\"Events/nx-3lep-BP2.lhco\",\n",
    "\"Events/nx-3lep-BP3.lhco\",\n",
    "\"Events/zw-3l-1.lhco\",\n",
    "\"Events/ttw-3l.lhco\"\n",
    "]\n",
    "\n",
    "df1 = get_pcd(evs[0],50,150)\n",
    "df2 = get_pcd(evs[1],50,150)\n",
    "df3 = get_pcd(evs[2],50,150)\n",
    "df = get_pcd(evs[3],50,150)\n",
    "df_1 = get_pcd(evs[4],50,150)\n",
    "\n",
    "def pcd_divide(df,N):\n",
    "    evs= df['event_no'].unique()\n",
    "    ev_nbins= int(len(evs)/N)\n",
    "    evt_lists= np.array_split(df['event_no'].unique(),ev_nbins)\n",
    "    return [df[df['event_no'].isin(evt_list)] for evt_list in evt_lists]\n",
    "    \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[39, 41, 12, 96, 2]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xsections=[7.90,13.53,2.408,435.1,1.429]\n",
    "kf=1.2\n",
    "Ncut=[len(df.groupby('event_no')) for df in [df1,df2,df3,df,df_1]]\n",
    "N=[30000,50000,30000,2000000,20000]\n",
    "lum=150\n",
    "Np=[int(lum*kf*x*Ncut[i]/N[i]) for i,x in enumerate(xsections)]\n",
    "Np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.1547005383792515"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "12/np.sqrt(12+96)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def get_combine(df1,df2):\n",
    "    return pd.concat([df1,df2],axis=0).sample(frac=1,random_state=42).reset_index(drop=True)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "pcds1temp=pcd_divide(df1,Np[0])\n",
    "pcds2temp=pcd_divide(df2,Np[1])\n",
    "pcds3temp=pcd_divide(df3,Np[2])\n",
    "pcds=pcd_divide(df,Np[3])\n",
    "\n",
    "n_pcds=min([len(pcds1temp),len(pcds2temp),len(pcds3temp),len(pcds)])\n",
    "\n",
    "pcds1=[]\n",
    "pcds2=[]\n",
    "pcds3=[]\n",
    "\n",
    "for i in range(n_pcds):\n",
    "    pcds1.append( get_combine(pcds1temp[i],pcds[i]))\n",
    "    pcds2.append( get_combine(pcds2temp[i],pcds[i]))\n",
    "    pcds3.append( get_combine(pcds3temp[i],pcds[i]))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# plt.hist(pcds[0]['mt']/pcds[0]['pt'], bins=50,range=[0,10],density=True,  alpha=0.5)\n",
    "# # plt.hist(pcds1[0]['mt']/pcds1[0]['pt'], bins=50,range=[0,10],density=True,  alpha=0.7)\n",
    "# plt.hist(pcds2[0]['mt']/pcds2[0]['pt'], bins=50,range=[0,10],density=True,  alpha=0.9)\n",
    "\n",
    "# plt.legend(['SM','BP1','BP2'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as cp\n",
    "\n",
    "def scale_it(df):\n",
    "    def standard_scale_column(column):\n",
    "        return (column - column.mean()) / column.std()\n",
    "    \n",
    "    return df.apply(standard_scale_column)\n",
    "\n",
    "def get_R1(pcd,rmax=0.05):\n",
    "\n",
    "    pos= scale_it(pcd[['px','py','pz']]).to_numpy()\n",
    "    v_weights=   pcd['pt'].to_numpy()\n",
    "    N=len(v_weights)\n",
    "    del_r =np.zeros((N,N))\n",
    "    for i in range(N):\n",
    "        for j in range(i+1,N):\n",
    "            del_r[i,j]=max(v_weights[i],v_weights[j])\n",
    "            del_r[j,i]=del_r[i,j]\n",
    "\n",
    "    neighbours =  (del_r < rmax) & (del_r > 0)\n",
    "\n",
    "    edges = cp.argwhere(neighbours)\n",
    "\n",
    "    pcd_t=pcd.reset_index(drop=True)\n",
    "\n",
    "    v1=pcd_t.loc[edges[:,0]][['mt','pt','phi']].reset_index(drop=True)\n",
    "    v2=pcd_t.loc[edges[:,1]][['mt','pt','phi']].reset_index(drop=True)\n",
    "    edge_weights= ((v1['pt']+v2['pt'])/np.maximum(v1['pt'],v2['pt']))**0.5 #/(v1['mt']+v2['mt']) #abs(v1['phi']-v2['phi'])# (v1['mt']/v1['pt']+ v2['mt']/v2['pt'])/2.0 # np.maximum(v1['mt'], v2['mt'])\n",
    "\n",
    "    # print(min(edge_weights))\n",
    "    # print\n",
    "    \n",
    "    v_weights= np.ones(N) \n",
    "\n",
    "    vR1= cp.bincount(cp.ravel(edges), weights=1.0/cp.sqrt(cp.repeat(edge_weights, 2)), minlength=N)\n",
    "\n",
    "    # print(N,v_weights.shape,vR1.shape,edges.shape)\n",
    "    \n",
    "    vR1= v_weights*vR1\n",
    "    \n",
    "    # vR1 gives sum(vw/ve) for every vertex; ve is the weight of incident edge. sum is over all such edges\n",
    "\n",
    "    # R=cp.zeros(N)\n",
    "    R_edges =cp.take(vR1,edges).sum(axis=1)\n",
    "    R_edges *= cp.sqrt(edge_weights)\n",
    "    R_edges = cp.take(2*v_weights,edges).sum(axis=1) - R_edges \n",
    "\n",
    "    # print(R_edges.shape)\n",
    "    R=cp.bincount(cp.ravel(edges), weights=cp.repeat(R_edges, 2), minlength=N)\n",
    "    \n",
    "    R_f=np.sum((del_r < rmax) & (del_r >= 0),axis=1)\n",
    "    return (R/R_f,R_edges,neighbours)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "from imblearn.over_sampling import SMOTE\n",
    "from imblearn.combine import SMOTEENN\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix, recall_score\n",
    "from imblearn.under_sampling import EditedNearestNeighbours\n",
    "from imblearn.under_sampling import RandomUnderSampler\n",
    "from imblearn.pipeline import Pipeline\n",
    "\n",
    "\n",
    "BP1_pcds=[]\n",
    "BP2_pcds=[]\n",
    "BP3_pcds=[]\n",
    "SM_pcds=[]\n",
    "\n",
    "for idx in range(n_pcds):\n",
    "    # print(f\"complted: {idx}\\n\")\n",
    "    BP1=pcds1[idx][['eta', 'phi', 'R', 'pt', 'mt', 'E', 'px', 'py', 'pz']]\n",
    "    BP2=pcds2[idx][['eta', 'phi', 'R', 'pt', 'mt', 'E', 'px', 'py', 'pz']]\n",
    "    BP3=pcds3[idx][['eta', 'phi', 'R', 'pt', 'mt', 'E', 'px', 'py', 'pz']]\n",
    "    BG=pcds[idx][['eta', 'phi', 'R', 'pt', 'mt', 'E', 'px', 'py', 'pz']]\n",
    "\n",
    "    BP_BG=pd.concat([BP1,BP2,BP3, BG], axis=0)\n",
    "\n",
    "    label_BP1=np.ones(len(BP1))\n",
    "    label_BP2=np.full(len(BP2),2)\n",
    "    label_BP3=np.full(len(BP3),3)\n",
    "    label_BG=np.zeros(len(BG))\n",
    "\n",
    "    BP_BG_labels=np.vstack( (label_BP1[:,None],label_BP2[:,None],label_BP3[:,None],label_BG[:,None]))\n",
    "\n",
    "\n",
    "    # smote = SMOTEENN(enn=EditedNearestNeighbours(n_neighbors=5))  #\n",
    "    # smote=SMOTE(random_state=42,k_neighbors=5,sampling_strategy={1: 100,2:0.1})\n",
    "    under = RandomUnderSampler(sampling_strategy='auto',random_state=42)\n",
    "    # steps = [('o', smote), ('u', under)]\n",
    "    # pipeline = Pipeline(steps=steps)\n",
    "\n",
    "    \n",
    "    vectors1_res, labels1_res = under.fit_resample(BP_BG, BP_BG_labels)\n",
    "\n",
    "    vectors1_res=vectors1_res.reset_index(drop=True)\n",
    "\n",
    "    BG_res=vectors1_res.loc[np.where(labels1_res==0)[0]]\n",
    "    BP1_res=vectors1_res.loc[np.where(labels1_res==1)[0]]\n",
    "    BP2_res=vectors1_res.loc[np.where(labels1_res==2)[0]]\n",
    "    BP3_res=vectors1_res.loc[np.where(labels1_res==3)[0]]\n",
    "\n",
    "\n",
    "    # # smote = SMOTEENN(enn=EditedNearestNeighbours(sampling_strategy='all',n_neighbors=5))  #\n",
    "    # # smote=SMOTE(random_state=42,k_neighbors=3)\n",
    "\n",
    "    # vectors2_res, labels2_res = under.fit_resample(BP2_BG, BP2_BG_labels)\n",
    "    # vectors2_res=vectors2_res.reset_index(drop=True)\n",
    "\n",
    "    # BG2_res=vectors2_res.loc[np.where(labels2_res==0)[0]]\n",
    "    # BP2_res=vectors2_res.loc[np.where(labels2_res==1)[0]]\n",
    "\n",
    "    BP1_pcds.append(BP1_res)\n",
    "    BP2_pcds.append(BP2_res)\n",
    "    BP3_pcds.append(BP3_res)\n",
    "    SM_pcds.append(BG_res)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Completed PCD-2\n",
      "\n",
      "Completed PCD-1\n",
      "\n",
      "Completed PCD-3\n",
      "\n",
      "Completed PCD-5\n",
      "\n",
      "Completed PCD-7\n",
      "\n",
      "Completed PCD-4\n",
      "\n",
      "Completed PCD-6\n",
      "\n",
      "Completed PCD-8\n",
      "\n",
      "Completed PCD-9\n",
      "\n",
      "Completed PCD-10\n",
      "\n",
      "Completed PCD-11\n",
      "\n",
      "Completed PCD-12\n",
      "\n",
      "Completed PCD-13\n",
      "\n",
      "Completed PCD-14\n",
      "\n",
      "Completed PCD-15\n",
      "\n",
      "Completed PCD-16\n",
      "\n",
      "Completed PCD-17\n",
      "\n",
      "Completed PCD-18\n",
      "\n",
      "Completed PCD-19\n",
      "\n",
      "Completed PCD-20\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from joblib import Parallel, delayed\n",
    "\n",
    "def get_R_pcds(i):\n",
    "    print(f\"Completed PCD-{i+1}\\n\")\n",
    "    BP1_res=BP1_pcds[i]\n",
    "    BP2_res=BP2_pcds[i]\n",
    "    BP3_res=BP3_pcds[i]\n",
    "    BG1_res=SM_pcds[i]\n",
    "\n",
    "    R_BP1=[]\n",
    "    R_BP2=[]\n",
    "    R_BP3=[]\n",
    "    R_SM=[]\n",
    "\n",
    "    for scale in [25,50,100,150,200,250,300,350,400,450,500]:\n",
    "        R1c,R1,_=get_R1(BP1_res,rmax=scale)\n",
    "        R2c,R2,_=get_R1(BP2_res,rmax=scale)\n",
    "        R3c,R3,_=get_R1(BP3_res,rmax=scale)\n",
    "        Rc,R,_=get_R1(BG1_res,rmax=scale)\n",
    "        R_BP1.append([R1c,R1])\n",
    "        R_BP2.append([R2c,R2])\n",
    "        R_BP3.append([R3c,R3])\n",
    "        R_SM.append([Rc,R])\n",
    "    return [R_BP1,R_BP2,R_BP3,R_SM]\n",
    "\n",
    "R_pcds=Parallel(n_jobs=-1)(delayed(get_R_pcds)(idx) for idx in range(n_pcds))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8190,)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R_pcds[6][0][1][1].shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAkhUlEQVR4nO3dzY8j550f8C9HLy0LUHeRE2ktxSPMVsubTZDAEsnG2gdjIYuMc8geYlWp/wE1CxJsICvDLPEU6GKKbQN7sCO4OL4b3UU7SGAYEVgtHxaGrPSwNHZ0WVmsETLOSNZkyGoKHrtXmqkcWlXDh+/VTbKK5PcDEGgWHxYfPmzWj897wvM8D0RERJ85F3UGiIgoXhgYiIhIwMBAREQCBgYiIhIwMBARkYCBgYiIBAwMREQkYGAgIiIBAwMREQkYGBZEtVpFJpNBIpHA5uYmVFWF4zhTOXc+n4eu61M5V5xpmoZkMhmUoaZpE5Wh/7xRkskkarVaqPxUq1UkEgmoqjo0ja7rSCQSwuezu7uLRCIx8JbP54NzJ5PJsbfd3d1QeV5EcfvcFwEDwwLY3d2FrusolUpoNpswDAOpVAqWZUWdtamyLAuqqqJarcKyLOi63vel293dRbVaRbVanfii5jgONjc34TgOTNNEu92GYRhotVoTfak1TYPrukPL2z+uKMpE+ekmSdLIPFSrVUiSJBy7efMmJElCu93uu5mmCQAoFApoNBrBzTRNuK6LS5cuCccLhULoPIflui6q1WoQtLrN6jMH4v25x55HsQfAM01zZufP5XJesVic2fknZZqmJ0mSB8CTZdkzDEN4vFKpeJVKJbhfr9e9QqEw9rzpdNrL5XJnypssy0NfS1GUifLRyzAMT5ZlL51OD/x86/V68Hj351MsFj1JkkK9VrPZ9AB4jUYjdD7PotFoeIZheJVKxUun032Pz+oz97z4fu6LgIFhAQDw6vX6zM4fp8DQbreHPi5JUt/j437bmKbpAfCazeaZ8lYsFge+VrvdPvUF1w8MlUpl4AWsUCh4xWIxssBgmqZXKBTOXHb+uYYFhml/5v554/q5LwI2JS2AXC4HTdPGNh3t7u5ic3MTiUQCmUwmSF+r1YT+iUmr0clkEpubm6hWq1N5H2fhOA5c1+1rVgEwslz29vaQTqchy/JErzPsfWuaBgB9Zbe/vw9ZlpFOpyc6/yCKosCyLLiuKxyvVqvB60ZBURRUKhUYhgFN02Db9lxf/7SfObAYn3ucMTAsANM0IUkS8vl80FnZ+yXVNA17e3tBW2qlUgkuNK1WC5cuXYLneTAMY+Dzu/kd21evXkW9Xoeu63O7KOzv76NWq6FarQodrsM6CyVJ6rugdnMcB9lsdqLXHvW+/YvA3t6e8Bz/onkW/rn39/eDY7VaDbIsD72wua47sPN52v1OkiQFwcHvD5j2a0z7M/efG/fPPc7ujToDNJ4kSWg0GrBtG3t7e0ENwDRNKIoSdO41m83gQpLL5YLnd3cw5nI5yLIMy7IG/tpxHAe1Wg3tdhuSJAUXBv8X2CCTfkEymczIzk7//P57qFarUFU16FAdJJVKodVqDX182C/OXpO8b03Tgg5J/+Jk2/bI/E1K0zQYhhGUz97e3shy9Tuf56lYLAK4+7lsb2+fueN1Fp85sDife2xF3ZZFp6MoiifLsud5dzvwRjEMw1MUxUun0x4AoUOvu4/Bb5uVJEm4nbUT7zT8dtx2u+3V6/WBbb2SJPV1WHbL5XJ9eVcUxQPgAQgem+R9+/nxX29Yh+qk/D6G7nP7beLdf8e189k0TU+W5ZF9BN1pJymraXzmnhfvz30RsClpQZVKpaANdhy/dqFpGhqNxth20XQ63TcMsl6vTynnw/W24/q/+BzHGdmkMqodOZ/P97Xfm6YJz/OCX8C+ce9bkiQoigLDMABMtzlBkiTkcrmgSSVM+/i8ua4LXdext7cHwzAm+mU+zCw+c2BxPve4YlPSAhj0JXEcJ6j6ptNpuK47NJ1t2/Am3ME1nU7Dtu2Jq+LAdJqSXNeFqqpCc5j/pZZlOXivg95jd7NZr2KxCMMwoOt68MUeZNL3rWlacNFxHGeq8wA0TUO5XEYqlYrlhcd1XZTL5SAwnDVwzeozBxbrc4+liGssNEaj0fAAeMVi0avX616z2Qyajrqr04VCwUun016z2fTa7bZnmqZXLBb7qsF+1XlYU5J/rlwuFzRlmKYppJ+V3iGzlUrFUxRFuN/9nv3hlOM0Gg1PkiRPURSv0Wh47XbbazQanqIoQpPBpO9bkiRPlmUhb938z2ic7qYkHz5r6uhunhnWlNRutwfehuUJp2xKajabXqFQONPQVcMwBja/zOoz97z5f+7LhIFhAdTrdU9RlGAi0LAJUcVi0ZNlOWhD9S8C/lhsSZKCL0H3l23QPAb/4uO3t85jvHa73Q4mNFUqlYFzKyqVimeaZhD4wpzbD574bDKVoih980Mmed+FQmHk3BLDMCbqAxgUGHondHne4MDgB5BBt0FOExj8gOD/wDiNZrMZtMn7P3C6/3dn+Zn755/X575MEp43YRsDEU3EbyKZR7/MLI1q56flxs5noinzx/svOgaF1cUaAxERCVhjICIiAQMDEREJGBiIiEjAwEBERIKVnPl8584dXL9+HQ899BASiUTU2SEimjnP8/Dxxx/jsccew7lzo+sEKxkYrl+/jgsXLkSdDSKiubt27Rq+8IUvjEyzkoHhoYceAnBSQOvr6xHnhoho9jqdDi5cuBBc/0ZZycDgNx+tr68zMBDRSpmk+Zydz0REJGBgICIiAQMDEREJVrKPYVK3b9/GJ598EnU2Yum+++7DPffcE3U2iGgGGBgG8DwPH3744UTbZq4ySZLw+c9/nnNBiJYMA8MAflB45JFH8OCDD/LC18PzPNy6dQsfffQRAODRRx+NOEdENE0MDD1u374dBIXz589HnZ3Y+tznPgcA+Oijj/DII4+wWYloiTAw9PD7FB588MGIcxJ/fhl98sknDAxEp/QP9XeF+3+f/6uIcnIXA8MQbD4aj2W0vF678ppw/8UnX4woJxQFBoYQeiP7LIX91aCqKizLgizLQae5oiioVCoD0wAnexP3piGiaN34wQ+F+w9/65tzzwPnMSyRUqmERqOBZrOJRqMB13WRz+eFNIVCAY1GI0jnOA40TRPS2LaNfD7PUVlEK4o1hiUlSRIqlQqSyeTIdJVKBZlMBoZhAAA0TUOr1YJlWfPIJtFcdTeRsXlsOAaGJeW6LnZ2dlAsFsem6+YHCPYfEK0uNiUtkXK5jM3NTSSTSSSTSeTz+ZH9B7ZtQ1VVlEqlOeaSiOKONYYlUiqVghqCbdvQdR2NRiOoBQBAtVoNmolSqRQqlQoURYkkv7SYOGJp+TEwLKl0Oo16vY5kMglN05BOpwGcdD5zFBKF1RsMaLmxKWnJybIMx3GizgYRLRDWGJZYtVqF4zjI5XJRZ4ViLA7j5ileWGNYIuVyGZlMBplMBpubm6jX62g0GpAkaeJz6LqOTCYDAMhkMlBVdUa5JaK4Yo0hhDisYTKMaZpTScP+B6L5euIXPxEPyNEv3skaAxERCSKpMTiOg1qtFnSMFgqFoc0dk6bVdR2lUilUswkRjXfRfKvv2Pvq30SQE5qXSAKDqqpoNBoATi78Ozs7Q5s5Jklr2zZ2d3c5UYuIaArm3pTUO3RSluWh6/JMmtZxnGDFUCIiOpu5BwbLspBKpYRjqVQKtm2fKm2tVuPMXSKiKZp7U9KwpZxbrVbotK7rTtSncHx8jOPj4+B+p9MZ+xwiolUVm1FJYdb+99Pu7+9PNHmrXC5jY2MjuF24cOGUuSQiWn5zrzFIktRXO2i1WgN/+Y9Ka1kWnnvuuYles1Qq4aWXXgrudzodBgciOrU47tM8TXOvMQz7hZ/NZkOn3d/fR7VaDZZ+KJfLA/sq1tbWsL6+LtyWjaqqSCaTwaznzc1N6Lo+NI2frjcNEdHcawy9o4ccx0E2mw1qDLZtQ5IkyLI8Mm1v0NA0DZqmzXZ00i/Lszt3r6fDD73tXnbbdV3ouo58Po96vR6k6V1dVVVVaJoWLM1tWVawi1s2m4VpmpwbQrRiIuljME0Tuq6jVqvBMAxhXkK5XEatVpsoLXByAdzd3QVwspzDoBrDKvK39hy3RWelUsH+/j6Ak7KsVCpoNptot9tIp9N45pln5pFdIoqRSCa4ybIc/GrtHWrae+EflRY4uQAWi8WxW1iumtNs7ek4jtC0VKlUkEgkOE+EaMXEZlQSnd1Zt/ZMp9NCE50/wZBBgWi1cHXVJTLtrT0rlQprYougt+/rFP1TRN0YGJbUWbf23N3dDfopaMEMGiTBYEEhsClpyZ1ma09N0xgUiFYYawxL7DRbe2qahnw+z/WniFYYA8MSKZfL2NvbA3Ay2iidTofa2tOyrGDCYLdGoxE0RRHF2WtXXhPuv/jki1NJu2oYGMKIcTvtNLb2zOVy8DxvWlkiip3eTYfisuHQB8fvBH8ffvggtj6/FWFu2MdAREQ9WGMgIoqR37dv4fatm8H9v4sgDwwMRHQmbKtfPmxKIiIiAWsMRCS48YMf4uKHh1FnY+5Y87mLNQYiIhKwxkBES6u3FnAxmmwsHNYYiIhIwMBAREQCNiWF0FstnaWwHV+qqsKyLMiyHGy+oyhK3zaefhrgZNmM3jRERKwxLJFSqYRGo4Fms4lGowHXdZHP54U0hUIBjUYjSOc4DjRNCx63bRv5fB6bm5vI5/OhV2YlmqfXrrwm3JbFB8fvBLco3hcDw5I67Z7POzs7ME0TzWYT+XweqqrOI7tEFCNsSlpSp9nzWZIkHBwcBKuxdjdLAQA6H4hPvj85ncwSfSYucwnebN7E8Y13g/t/n/+rSPIRFQaGJVIul2EYBlqtFlzXhWEYKBQKQ9P37vkMIAgKjuOgXC5D1/VZZ5uWzLimj1WeOLYo2JS0REqlEprNJtrtNhqNBkzTFPoPgJPNezKZDDKZDHRdH7ivs6Zp2NzcRCqVGhlYiGg5scawpM6y57NhGDAMIwgijUZjHlkmElw038KNf7wT3H/4W9+MMDerhYFhyfl7Pp9mB7ZCoQBN02BZVqjtQc+Mm9kTRYqBYYmF3fPZsiw4jhM0H1WrVUiShGw2O8ts0gq68YMfCvdZG4gXBoYlctY9n3O5HKrVKjY3NwH0j1IimoXDDw/x/hLNQVgGDAwhxHk0xTT2fAZOmo/Y4Uxx82vnJt6rDx8++mbzpnD/K5vnBz5+b/sWAOALyQdnkc2lwcBARLHzZvMm7u3Z4pLmh4GBlkdvpzU7rGlK/qGrtrIKGBiIaCldNN8K/vabkGgyDAxEFFr3RRcA3lf/Zi6va3dOBle8duX8mJR0Fpz5TEREAtYYiFaU327+hHPSwftlebF/hX/14LfB31+wOeroLFhjICIiQSQ1BsdxUKvVguUaCoXC0ElUo9L6ew24rovDw0Nsb2+faukHIjqZK+D74PhWqLH+UfU5zJPfv+FLr29HlJPZiyQwqKoaLMzmOE6wOUzYtKqq4uDgALlcDq1WC6qqotlszizfvdP4Z+k0SwT4ezD4ATObzcIwjGArT03TUK1W0Ww2g2Pdz00mk1AUZaKJcLT8ft8zkiduk8I46mh25t6U1LtVpCzLQ3cZG5fWNE2hhrDqSzdkMhnk83m02220221omiZutIOTMqzVan3PrVarfcGCKK7ebN7E79u3ghtN19wDg2VZSKVSwrFUKgXbtkOn7V4cbtDeA6vEcRy0Wi1hOQtFUfqa1hRFgWEYfc/f29uDoigzzyfRInniFz8Jbt2d28tu7k1Jvb9gfa1W61RpbdvG3t4e8vn8Sq/xk0ql4LouarXayAv8+fPnIcsybNsOgoZfM9vc3OyrpRHR2fTOml6EbUJjMyppWBAYlzadTgc7lw1qIgGA4+NjdDod4bZsJEmCaZpQVTXYnGdQLQw46ZvprjUYhrHStS2aju6mnd5F7U7rg+N3gltv5+9ZdOdzWnldJnOvMUiS1Fc7aLVaA/sHJk0rSRJUVQ3a13sfL5fLeOWVV6aR/VhTFAXtdhuWZWFvbw+ZTAb1er1vP4ZCoYBkMhkEh2q1iqtXr2J/f38+GeWaRivJvwAf33gXTzg3Q8+bWJamnN4axBMR5WOUudcYhm0aM2gzmFFpLctCMpkMjvkdp4OaQkqlEo6OjoLbtWvXTpP1hSBJUjCyqFKpDN3GM5fLwbKsYHe2Ve+4JzoNv//hxg9+GNyWwdxrDL0jXxzHQTabDS5Mtm1DkiTIsjwybSqVEgKH/7xB8xjW1tawtrY2/TcTc+l0Oti4p5emaTAMA6lUCtvbyzsem4jCi2Qeg2ma0HUdW1tbODw8FMbNl8tlbG1toVgsjkybTqexvb2NarUKAKjX6yu9ab1lWUFZybIM13Wh6/rQi34ulwv6FQaNUoodNj9FLu7zGmh6IgkMsiwHTRy9I2h6J1eNStt9f5VHJAEnF3rbtpHP5+E4DmRZhqZpQYAdRNf1UJ3+RLNy0XxLCDxfjTAvxEX0Qon7huXFYnFkIOitGfQGU27rOUW9NRyAtRxaGAwMtFrYJLUSOBv6bBgYiFZE3Lan9Bft++B4+S/ii7YAX2wmuBERUTywxkBEM3fRfCtYAfXRtTsR54bGCR0Y/HHvmqbhySefnEGWaKkM6oQlolgLHRj8pROef/55HB0dQVEUaJqGixcvziB70blzh79qxol9GTEoRearB79diO01/WU2/Lwu4wZDpxE6MGxsbGBnZwc7Ozs4OjrC/v4+VFUFAGxvb6NQKGB9fX3qGZ2X+++/H+fOncP169fx8MMP4/7770cikYg6W/Hwz58AADwP+OdPb+PGH/+Ic+fO4f777484Y1PGkUsry9/85972rZVu8jp1H8OVK1dgGAbq9TrS6XSwKYyiKPj617+Ob3/729PM59ycO3cOf/mXf4kPPvgA169fjzo78fLno7t/e3fw4Pl/iccffxznznWNYViFX+lRzlFYhfKlyIUODN/73vdgGAYSiQQ0TcOrr76KjY2N4PFnn30W2Wx2YQMDcFJrePzxx/Hpp5/i9u3bUWcnPt66O0HuHu8T3PvUt1mbIlpCoQNDs9mEaZp46qmnhqZZhrX9E4kE7rvvPtx3331RZyU+7vSMN2dQIFpKoQPDc8891xcUjo6OcHBwgG984xsAgJ2dnenkjmiZTNJ3wf6Nhde39PZf//toMnIGoQODrus4PDwUjm1sbKBcLgeBgVYE27tnj2UcmQ+O34k6C5GZODD8+Mc/BnCyJ4L/t++9995b6SWvaYHxwntqXI9oeU0cGN577z0AgOd5uHz5svCYJEkMDEREE/jqwW9jPxR24sDw6quvAgBc18WPfvSjmWWIiIiiNVFguHLlSrD8xcsvv4z3339/YLplm/1MtCpWuT2d+o0NDG+//Tay2Wwwnj+dTiORSMDzPCFdIpHAzZs3Z5NLWm5s55+9X5bx5f9z9/v568ej25ApjkFo3v0lcSyDbmMDw1NPPSVM8mq1WjPNENFEGEyIZib0cNXvf//7UBQFFy9exBtvvBFMZqtWq3j66aennkGKCC+8RKfib0AU+Oto8nEWoTfqMQwj6EtQVRW6ruP111/nXsFERJ/54Pid4Na7e9siCF1j8PsWDg4O0G638fzzzwMA+xeIiGbgovkWbvyjOLz14W99c6avGTowyLKMUqkE0zSDWsLVq1chy/LUM0dERPMXuinJNE3Isgxd14P5DK7rolTimi5ERMvg1Bv1dHvqqadGrrZKRESLI3Rg+NnPfoZisYirV68GxzzPQyKR4N4FtLo4iouG8LcPXSShA8POzg6q1Spyudws8kNES2KVFtnrG6K64EIHhkwmg2effXYWeSEiirXuX/9xXwjvLEIHBlVVsb29je3tbciyDEmSgse4VhIRraq4L3MRRujAUKlU0Gq1UK/XheNcK4mIaPp+376F27fuXlu/LJ+f+WuGDgz+vgxERLScQgcGAHjjjTdgmiaSySS++93v4uDgAIlEAl/72temnT+i6HHEEa2Y0BPcSqUSDMNAoVAImpMymQx0XZ965oiIaP5C1xhqtRp+97vfAbi7bpIkSXBdd6oZIyKiaJxq5vNvfvMbfOlLX0IikQBwspnPxsbG1DNHtNTYREUxFTow7O/vI5vNYmtrC+12Gy+88AIsy4JpmhOfw3Ec1Go1yLIMx3FQKBSEYa+TprVtG5ZlAQAODw9x6dKloechIqLJnGp11VarhZ/+9KdwHAeSJOHVV18NVWNQVRWNRgPAyYV/Z2dnaGAZldayLBSLRQDA7u4unnnmmSAtEdEsLdO8hV6hO58BoNPpIJ/P4zvf+Q52dnZCBQXHcYT7siwHv/rDpLVtG+Xy3aq4oiiwbbvvOUREFE6owPDyyy/j/PnzSCaTSCaTuOeee/Diiy+i0+lMfA7LspBKpYRjqVQKtm2HSptOp3Hp0qXguN/53ZueiIjCmTgwPPHEE3j77bexv7+PVquFVquFw8ND3LhxA5lMZuLgMGz0UqvVCp1WUZTg2N7eHnK53MA+huPjY3Q6HeFGRESDTRQYXnjhBSiKgtdffx3PPPMMNjY2sLGxgXQ6DdM0sbOzc+Z5DGGGu/amdV0XtVptaD9FuVwO8ryxsYELFy6cIadERMttos5ny7KCuQuDFItFfPGLX5zoBSVJ6qsdtFqtgb/0J02r6zrq9frQEUmlUgkvvfRScL/T6TA4EBENMVGNYZLO5Uk7oIft45DNZk+Vdnd3F7quQ5ZluK47sOaxtraG9fV14UZERINNFBj8iWxnTQOcjCzq5jgOstmsMDfBH1k0Lm2tVkM6nQ6Cwv7+PucxEBGd0URNSc1mE6VSaWSaMMNETdOEruvY2trC4eGh0DdQLpextbUVzE8YltZxHKiqKpxXkiQUCoWJ80FdOAuXiD6T8PwFj0ZQVXWiGsH+/v5UMjVrnU4HGxsbODo6YrOSj4Fh6b3Ztf3krx+/+wPqiV/8ZKknay2DR9f+bfD3l+XzePhb3wx9jjDXvYlqDGGWuyAiosV2qpnPRES0vBgYiIhIwMBAREQCBgYiIhIwMBARkYCBgYiIBAwMREQkYGAgIiJB6K09iYhovsSZ6X8789djjYGIiAQMDEREJGBgICIiAQMDEREJGBiIiEjAwEBERAIGBiIiEjAwEBGRgIGBiIgEDAxERCRgYCAiIgHXSiJaIm86N4X7X5HPR5QTWmSsMRARkYA1BqIl1luDAIAHfvUunnjnJxHkhhYFawxERCRgYCAiIgEDAxERCdjHQLTkHvjVu1FngRYMawxERCRgYCAiIgGbklbRL8tR54CIYow1BiIiEjAwEBGRgIGBiIgEkfQxOI6DWq0GWZbhOA4KhQIkSTpVWtu2sbOzg0ajMZ/ML4LePoSnS9Hkg2LjnxIt4X7n+J2IckKLIJLAoKpqcCF3HAc7OzswTTN0Wj9g2LY9n4wTEa2AuQcGx3GE+7Isw7KsU6VVFGX6GSQiWnFz72OwLAupVEo4lkqlBv7qD5OWiIimY+41Btd1Bx5vtVp9x8KkHeX4+BjHx8fB/U6nE+r5RESrJDajkoYFgbOmBYByuYyNjY3gduHChXCZIyJaIXMPDJIk9f3ib7VaA0clhUk7SqlUwtHRUXC7du1a2GwTEa2MuQeGXC438Hg2mz1T2lHW1tawvr4u3IiIaLC5BwZZloX7juMgm80GtQDbtoPRSOPSdgvbvERERINFMo/BNE3ouo6trS0cHh4KcxjK5TK2trZQLBbHprUsC/V6XXgeh7ASEZ1NwvM8L+pMzFun08HGxgaOjo6Ws1lp3Mxnrq66tN50bvYde+BX7/bPfF57dF5Zoin7T//6b/Hwt74Z+nlhrnuxGZVERETxwMBAREQCbtRDtMCGNR0RnQVrDEREJGBgICIiAQMDEREJGBiIiEjAwEBERAIGBiIiEjAwEBGRgIGBiIgEDAxERCRgYCAiIgEDAxERCRgYiIhIwEX0Fs24vRYmeQ4tjN5F8r4in48oJ7RKWGMgIiIBAwMREQkYGIiISMDAQEREAgYGIiISMDAQEZGAgYGIiAScxxBnnH+wcm78/Ip44N9cCP78H+feQ7t2Obi/9cBfCI8TTQsDQ5wwENCU/FOiFXUWaIExMBDF2AO/ejf4+0uJFvAAv7I0e/wvW3SsZayUzp8/Df7+/Z/+FGFOaJkxMBBFqLtP4fDPfxAe23rgL8Y+v7tGQTQtDAxEM9LXkQzg4f/45MjnnKVGwH4FmhYGBqI5GhQsiOKG8xiIiEjAGkNU2Gm8lF5zfxv8ffHPf+jrJzhw/2/w97/yUsJjncSnIIoDBgaiiIzrEzjr40SnxaYkIiISRFJjcBwHtVoNsizDcRwUCgVIkhQ6bZjzRI5NR0S0ICIJDKqqotFoADi5uO/s7MA0zdBpw5yHKAq9cxOIFsHcm5IcxxHuy7IMy7JCpw1zHiIimtzcawyWZSGVEkdjpFIp2LaNdDo9cdrLly9PfJ6p620Wero0Pg0tvRs/v4KLrCHQEph7YHBdd+DxVqt/hMWotGHOc3x8jOPj4+D+0dERAKDT6YzO7DB//LN4/+f/5XTnoVj7f//zfwv3/8V/+Hd9af7X+3f/39b+0MGtz4lDTrtnMsfJnxKfRJ0FOqWP//QnrJ3i2uVf7zzPG5s2NsNVh13ow6Yd9Fi5XMYrr7zSd/zCBa5lTyH81/8WdQ6I8J/x3wG9eOrnf/zxx9jY2BiZZu6BQZKkvl/1rVZr4GiiUWnDnKdUKuGll14K7t+5cwetVgvnz59HIpE4/ZtZIp1OBxcuXMC1a9ewvr4edXZiiWU0HstovKjKyPM8fPzxx3jsscfGpp17YMjlcjAMo+94NpsNlVaW5YnPs7a2hrW1NeFYbIe1Rmx9fZ1f6DFYRuOxjMaLoozG1RR8cw8MsiwL9x3HQTabDS7Utm1DkiTIsjwybe+Fvfc8RER0OpH0MZimCV3XsbW1hcPDQ2HuQblcxtbWForF4ti0ox4jIqLTSXiTdFHT0js+Pka5XEapVOprdqMTLKPxWEbjLUIZMTAQEZGAi+gREZGAgYGIiASxmeBG8+evLeW6Lg4PD7G9vS0sJ2JZFhzHCUaH5XI5AAu2qu0ZjSsjn67rKJVKi7ny7xmNKiPbtoPHDw8PcenSJZZRTxnFcgVpj1aWJEleo9HwPM/zDMPwZFkOHqvX616hUPA8z/OazabwWDqdDv5uNpueoihzyvH8jSojX6PR8AB47XY7OMYyOlGpVIS/u8uFZXRiVDlEVUZsSlphpmkKv367f4lomoZKpQLgZO5JvV4HsHqr2o4qI193rcq/321Vy8i2bZTLdxeTVBQFtm3DcRyWUVeNoFtcVpBmYFhhftMQcPKPq2kagJN/SH95Edu24bpucOEbteLtMhpWRr5arQZFUYRjLKOTMkqn07h06VLwmL+OWSqVYhl9VkajyiHKMmIfw4qzbRt7e3vI5/MoFArBsVQqhVqthlwuh2q1ClmWoShKqFVtl8WgMgJOLnSDahAso7tl1B009/b2kMvlIEkSy+izMprWCtLTxhrDikun0yiVSmg2m6jVagBO/vEcxwm+xIVCAaqqjjxPmNVxF82gMgKA/f194ZfgOKtYRj7XdVGr1cauTrDKZdQt7ArS08Yaw5KpVqtoNptDH8/n830XM0mSoKoq8vk82u02ZFkW1qPqXcdq0lVt42oaZXT58mU899xzA5/PMjopo+73q+s66vW68D/FMmpPbQXpqZtLFzfFTr1e9yRJCu43m00PgNdoNLxmsyk85nme8Fj3SAnPOxlx0T0iZ1mMKqN6ve4ZhhHcAHjFYpFl1FVGvkql4jWbTc/zPK/dbnvtdptlNMH3KcoyYo1hRaVSKeHXjF8b8EdOZLPZoA3dH3UzaPz+Mq9qO66MummaBk3T+lYEBla7jGq1GtLpNGRZhuu62N/fHzgWf5XLqFtcVpDmWkkrrFarBVXVer2OSqUSXNhc14Wu68hkMmg0GtB1PXjMcRwYhhGsats9sWvZjCoj4KScqtUqdF1HoVCApmlIp9Mso88mZG1ubgppJUlCu90GwP+jSb5PUZURAwMREQk4KomIiAQMDEREJGBgICIiAQMDEREJGBiIiEjAwEBERAIGBqIZ2t3dRSKRgKqq2N3dha7rUFV1aVcRpeXAeQxEM+S6LpLJJLq/ZrZtI5PJgF89iivWGIhmyLKsgUsfAMu9kigtNgYGohmq1+vIZrPCsXK5vNT7G9PiY1MS0Qxtbm5CURTk83k4jgPXdZFOp0Pt40A0bwwMRDPi9y/46+47jgNVVXFwcMDaAsUam5KIZuTy5cvBpkfAyWbuqVQK+/v7Qjp/JVuOVKK4YGAgmpF6vd7XZHT58uW+dJcvX2ZHNMUKAwPRjFiWhXw+LxxzXTdYh9+vIfh7axPFBQMD0ZQ5joPd3V3IsozDw0M4jhM8ViwWUa/XUa1Whw5jJYoaO5+JYkDXdWxvbzNYUCywxkBERIJ7o84A0aqzLEsYkcRaA0WNTUlERCRgUxIREQkYGIiISMDAQEREAgYGIiISMDAQEZGAgYGIiAQMDEREJGBgICIiAQMDEREJGBiIiEjw/wF+rOJDXSx83AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAApV0lEQVR4nO3dXYzj5r0e8Ee7tsdx4RmOtjFiYzdYc+zm5gDOUhIQA/WFPRJ8VxQJ5blp7iwRNs65OEkgWlfB9iKylAC5SOCEmqC9Kg5mqKRogR60ENdGcVIcBDOiHaNA0MTiJmfdtWHHEkfucXbi2OzFLLl69UGJMxp9Pj9gsCPxJflS3OFf73fM8zwPREREd12YdQaIiGi+MDAQEZGAgYGIiAQMDEREJGBgICIiAQMDEREJGBiIiEjAwEBERAIGBiIiEjAwLKBqtYpEIoFYLIatrS1ks1k4jjORY2cyGei6PpFjLTJN07C5uRl8xpqmjfUZ+/uF2dzcRK1Wi5SfarWKWCyGbDY7NI2u64jFYsL9q1QqiMViA38ymUxw7M3NzZE/lUolUp4X0bzd91m5b9YZoGgqlQpKpRJ2d3ehKAocx4FpmrAsC/l8ftbZmyj/QdRsNgEAhmEE2yzLgmEYyGQykGUZ9XodqVQKqqoK+0uSBABwXReFQmHkOR3HCY5pmiaSySQODw9hGAZqtdrIY2iahmq1CsuykE6n+7ZblgUAQj7HJUlS6IOlWq0G1+v76KOPIEkSbt68OXS/fD4v5NX/DEzThKIowfvxeDxynk/DdV3s7+/DNE3U63Vh2yre95nwaKEA8EzTPLfjp9Npr1AonNvxx9Wbh3w+76XT6eC1aZqeJEkeAE+WZc8wDCF9uVz2yuVy8Lper3v5fH7keRVFEc5zGrIsDz2Xqqpj5aOXYRieLMueoigD73+9Xg+2d392hULBkyQp0rmazaYHwGs0GpHzeVaNRsMzDMMrl8ueoih921ftvs8KA8OCAeDV6/VzO/48BIZ2u+2l02mv3W4H7zUaDQ+A12w2Pc87eUB0b+8lSVLf9lHfg0zTFM5xWoVCYeC52u32qR+4fmAol8sDH2D5fN4rFAozCwymaXr5fP7Mn1338YYFhlW677PCNoYFk06noWlaUDQdplKpYGtrC7FYDIlEIkhfq9WE9olx6jz9+tOtrS1Uq9WJXMcoh4eHQt2uLMsATqoGRnEcB67r9lWrAAj93Pb29qAoSnCuUYZ9LpqmAUDfZ7u/vw9ZloXqmahUVYVlWX2fQ7VaDc47C6qqolwuwzAMaJoG27annodlvu/TxsCwYEzThCRJyGQyQWNk7x+hpmnY29uDaZpot9sol8vBg6TVamF3dxee58EwjIH7d/Mbtm/evIl6vQ5d18/9j16SJLTbbeEPyf/D7v7j3d/fR61WQ7VaFRpchzUWSpIUGlgcx0EymRwrj2Gfi/8Q2NvbE/bxH5pn4R97f38/eK9Wq0GW5aEPNtd1BzY+j/pyEZUkSUFwsCwL2Wx24ucAVvO+TxsbnxeMJEloNBqwbRt7e3tBCcA0TaiqCtd1Ua1W0Ww2gwdFd2NYdwN1Op2GLMuwLGvgtxnHcVCr1dButyFJUvCH73/DGmTcP4BEIhGpsbxUKsEwjODboH9+/xqr1Sqy2SxM0xx6jHg8jlarNXT7sG+bvcb5XDRNg6ZpwTFd14Vt26H5G5emaTAMI/j89vb2Qj93P9BOk99Y69+XnZ2diTS8rvJ9n6pZ12XR2amq6smy7Hnevca5MIZheKqqeoqieACExrruNga/7lWSJOHnrI10URUKBSGPg/j1uO1226vX6wPreiVJ6mus7JZOp/uuTVVVD4AHINg2zufi58c/37DG1HH5bQzdx/brxLt/n9fGZ9M0PVmWQ9sHetOP83kt+32fFVYlLYFisRjUr47ily40TUOj0RhZ76koCtrttvDT24XwPNVqNWxtbfV1F+ytx/W/8TmOE1qlElaPnMlk+urvTdOE53l95x/1uUiSBFVVgy62k6xOkCQJ6XQ6qE6JUj8+ba7rQtd17O3tCSW+01rl+z5NrEpaMIP+ABzHCYq2iqLAdd2h6Wzbhjfmaq6KosC27bGL2sBkq5L8+mk/neu6aLVaiMfjyGazQnWZ/0cty3LwWQz6DAb1MfcVCgUYhgFd14UxE73G/Vw0TQseOo7jTHSciaZpKJVKiMfjc/ngcV0XpVIpCAyTCFyu6678fZ+aGZdYKAK/y2ahUPDq9brXbDaDqqPuonI+n/cURfGazabXbrc90zS9QqHQV8z1i8bDqpL8Y6XTaaGb6KhqnUlda7lc9prNZvBTLpeDqojeLrXlctlTVVV43f2Z+N0pxzmvJEmeqqpeo9Hw2u2212g0PFVVhSqDcT8XSZI8WZaFvHXz7+Eo3VVJPtyt6uiunhlWldRutwf+DMsTTlmV1Gw2vXw+f+auq4ZhDKyCWZb7Pu8YGBZMvV73VFUNBvkMG/BUKBQ8WZaDOlL/j9zvay1JUvCfvPsPadA4Bv/h4tennnd/7Ha7HVxf7093Gn8wU7lcHjj2olwue6ZpBoExyvn94Iq7A6lUVe0bPzLO55LP50PHnhiGMVYbwKDA0DuYy/MGB4ZBn2Pv59ntNIHBDwj+F5DT8r8A+J99oVAQ/n8vy32fdzHPG7NegYgmzq8emWa7zXkIq+OnxcPGZ6IZ8vv7LzoGheXCEgMREQlYYiAiIgEDAxERCRgYiIhIwMBAREQCjnyO4PPPP8ft27fx8MMPIxaLzTo7RERj8zwPH3/8MR577DFcuBBeJmBgiOD27du4cuXKrLNBRHRqt27dwuXLl0PTMDBE8PDDDwM4+WDX19dnnBsiovF1Oh1cuXIleI6FYWCIwK8+Wl9fZ2AgooU0TjU4G5+JiEjAwEBERAIGBiIiErCN4Rx89tln+PTTT2edjbl0//334+LFi7POBhGFYGCYIM/z8P7774+1xOYqkyQJX/rSlzgWhGhOMTBMkB8UHnnkETz00EN88PXwPA+ffPIJPvjgAwDAo48+OuMcEdEgDAwT8tlnnwVB4dKlS7POztz6whe+AAD44IMP8Mgjj7BaiWgOMTBMiN+m8NBDD804J/PP/4w+/fRTBgaiMG+UBr//bPFcT8teSRPG6qPR+BkRzTeWGKbgh/XfTu1cf5v5V2OnzWazsCwLsiwHDeaqqqJcLg9MA5ysUdybhoiWC0sMK65YLKLRaKDZbKLRaMB1XWQyGSFNPp9Ho9EI0jmOA03ThDS2bSOTybBHFtESYImBApIkoVwuY3NzMzRduVxGIpGAYRgAAE3T0Gq1YFnWNLJJROeMgYECrusil8uhUCiMTNfNDxBsOyBaDgwMK65UKsEwDLRaLbiuC8MwkM/nh6a3bRvZbBbF4vn2iiBaRb3tkV/7p4+C35+Wp9cNnoFhxRWLxaCEYNs2dF1Ho9EISgEAUK1Wg2qieDyOcrkMVVVnkl8iOn8MDBRQFAX1eh2bm5vQNA2KogA4aXxmLySi1cFeSdRHlmU4jjPrbBDRjDAwkKBarcJxHKTT6VlnhYhmZCZVSY7joFarBd9M8/k8JEk6VVrbtpHL5dBoNIT9bNsGcFI94jgOXNcNqkainH/ZlUol7O3tAUDwGTUajUifh67rQRtEIpGAoigwTfM8sktEUzCTwJDNZoMHueM4yOVyQx8kYWn9h7sfBLoZhoFqtQoASKfTwvGjnH8SooxGnqZxrnmcNGx/IFouUw8MvXXXsiwPHRg1Km1Yz5hEIoF2uw0AwrffKOcnIlpFU29jsCwL8XhceC8ejw/81h8l7SCSJPVViZz1mEREy27qJYZhc+m0Wq0zpR20b61WAwAcHBxA0zRhsrhxjnl8fIzj4+PgdafTGXleIqJFNzfjGKJMvjZO2u4GZVmWkclk0Gw2Ix2zVCrh+vXrY+eLiGgZTL0qSZKkvm/nrVZrYC+YKGl7dbcl+L2PHMeJdMxisYijo6Pg59atWyPPS0S06KZeYkin08J0C75kMnmmtN1s28b29nbQ+OyLx+ORjrm2toa1tbXQcxERncZrb73W957dOZkbSVnfmXZ2BFMPDP6CLz7HcZBMJoNv7LZtQ5IkyLI8Mm0313WFqqPuLpSWZUFV1YGN0WHHJCKapmduvA0AeHTtcwDAg0ddk+rJT08tHzNpYzBNE7quI5VK4eDgQOgrXyqVkEqlgondwtJaloV6vS7s5weAZDKJSqUCSZLQbDaF/cKOSUS06mKe53mzzsSi6HQ62NjYwNHREdbX14Vtd+7cwc2bN/H444/jwQcfnFEOo5nV0p6L+FkRTdqgqqT7qm8AAB5d+ysAwOWjezM6XPtmV4nh2ejT3oc9v/ryEfnoFN0bpemdK+J/mO5pt13Xha7ryGQyQUkM6J9dNZvNQtO0oK3GsqxgFbdkMgnTNFk1R7TAOIkeBfylPUeNBC+Xy9jf3wdwEkzK5TKazSba7TYURcH29vY0sktE54QlBgqcZmlPx3Gg63rwulwuIxaLwXGcvs4DRLQYWGJYcaVSCVtbW9jc3MTm5iYymUxo+0Hv0p6KoghTdPvjRxgUiBYXSwwrbtJLe5bL5ZElDiKabwwMFDjr0p5+92BOw0202FiVRH1Os7SnpmkMCkRLgiUGEpxmaU9N05DJZELXxyCixcESw4orlUpIJBJIJBLY2tpCvV6PtLSnZVmoVqvIZrOIxWLBD9e3IFpcLDFMwylGKU7DJJb2TKfT4OB5ouXCEgMREQkYGIiISMDAQEREAgYGIiISMDAQEZGAgYGIiAQMDEREJGBgICIiAQMDEREJOPJ5Cgat7XpeXv7qy2OnndWaz0Q032ZSYnAcB5VKBbVaDZVKRVgRLGpa27aRSCT69rNtG5VKBZVKBdlsVtjPtu1gLh/HcVZ6Xp9isYhGo4Fms4lGowHXdZHJZIQ0+XwejUYjSOc4DjRNC7bbto1MJoOtrS1kMpnIM7MS0XyZSYkhm82i0WgAOHkw53K5oXPyhKWt1WqQZXngg92yrGDBmEqlgu3t7eA4hmGgWq0COJnrZ5w5g1aBP2325uZmaLpyuYxEIgHDMILlQG/cuAFJkoJA7H/WRLR4pl5i6P02Kcvy0MXnR6VVVTVYTKabbdsolUpCOtu2g+MlEgm02220223U6/WxZxJddqdZ81mSpCAoABCqpYhoMU09MFiWhXg8LrwXj8eHfusfN203RVGwu7sbvPYfVN3HkiSJAQFnX/MZQPA5Oo6DUqkEXdfPO9tEdI6mHhiGfZtstVpnSture9GYvb09pNPp4AHmui5qtRpqtRp0XR9aJ358fIxOpyP8LJtisYhms4l2u41GowHTNIX2A+Bk8R5/zQZd1weu66xpGra2thCPx5HP56d5CUQ0YXPTKylK9UPUtLVaTajzzufzQtVHJpNBs9ns27dUKuH69etjn2vRnWXNZ8MwgrabRCLBNgaiBTb1EoMkSX3f+Fut1sBqnShph9F1va8dobuE4K9vPKjUUCwWcXR0FPzcunVr7PMustOs+ezL5/OwbXtouxERzb+pB4Zhawknk8kzpR2kUqlA1/WgQdR1Xdi2je3t7b60vW0ZALC2tob19XXhZ9lFXfPZX9qze39Jksa+R0Q0f6YeGPyBUj7HcZBMJoNv9N29h0al7dZbvVSr1aAoShAU9vf3IUkSZFkWqkUsy4KqqivbEH3WNZ/9ALK1tYWtrS0YhiH0UiKixTOTNgbTNKHrOlKpFA4ODoRxBKVSCalUKmjcDEtrWRbq9bqwn6qqcBwH2WxWOKckSUHbQjKZRKVSgSRJaDab5z6OIcpo5GmaxJrPwEn1ERuciZZHzONK7mPrdDrY2NjA0dFRX7XSnTt3cPPmTTz++ON48MEHZ5TDxcDPimjwVDn3Vd8AADy69lcAgMtH9zpxXPvm0/cSPltEVGHPr16cRI+IiAQMDEREJGBgICIiAQMDEREJGBiIiEjAwEBERAIGBiIiEszNJHrL7MMf/Xhq5/ri3/x1pPT+Ggz+3EbJZBKGYQSjzjVNQ7VaRbPZ7BuJ7rouNjc3oaoqFzsiWiIsMay4RCKBTCYTLFykaVrf9CKyLKNWq/XtW61W+4IFES0+BoYV5jgOWq2WMJ3FoFXxVFWFYRh9++/t7QnrXhDRcmBgWGHxeDxYryLMpUuX+tbW9ic63NraOtc8EtH0MTCsMEmSYJomstlssDjPsGVTs9msUGowDKNvpTciWg4MDCtOVVW0223s7u6i1WohkUgMXGQnn89jf38/eF2tVvHCCy9MM6tENCXslUSQJAmqqkJVVVQqFZTL5YEL9aTT6SBodK+hTUTLhSUGEiiK0recqk/TNBiGAdM0sbOzM+WcEdG0MDCsMMuyoGla0JDsui50XR/60E+n08F6zuyNRLS8WJW0wvwHfSaTgeM4kGUZmqYFq+cNout63zgHIlouXMEtAq7gNhn8rIi4ghsRES2QmVQlOY6DWq0GWZbhOA7y+fzQHi6j0tq2jVwuh0ajMfZ+Uc5PRLRqZhIYstls8CB3HAe5XG7oJGxhaf2H+6BBWWH7RTk/EdGqmXpg8HvA+GRZHjigapy0w3rGhO0X5fxERKsocmCIx+PY2dmBpmn46le/GvmElmUhHo/3HdO27b7J26KkHXe/w8PDUx2TiCiqQQ3MiyBy4/PNmzehKApefPFFPPnkkygWi/j9738/9v7DujoOGlQVJe24+532mOP6/PPPJ3KcZcbPiGi+RS4xbGxsIJfLIZfL4ejoCPv7+8hmswCAnZ0d5PP5kV2hBonSN/60/ejD9hu07fj4GMfHx8HrTqczdP8HHngAFy5cwO3bt/HFL34RDzzwAGKx2Knyuaw8z8Of//xnfPjhh7hw4QIeeOCBWWeJiAY4dRvDW2+9BcMwUK/XoShKsMCLqqp4/vnn8e1vf3vgfpIk9X07b7VaA3sFRUk77n5RjlkqlXD9+vXQc/kuXLiAxx9/HO+99x5u37491j6r6qGHHsKXv/xlXLjA3tJE8yjyALfvf//7MAwDsVgMmqYhl8thY2NDSJNMJnF4eDhwf8dxhF5BALC5uYmbN2/2PZzHTRuLxdB9GWH7tVqtsc8/qMRw5cqV0AEinufhL3/5Cz777LOB21fdxYsXcd9997E0RSthUBvDVfNX4hvuH4JfO3f+AgD4iie2gwLTHeAWucTQbDZhmiauXbs2NE3YPP29S0E6joNkMhk8lG3bhiRJkGV5ZNpurusG74ftNyj4DDvm2toa1tbWhl7LILFYDPfffz/uv//+SPsREc2LyIHhhRde6AsKR0dHuHHjBr7+9a8DAHK5XOgxTNOErutIpVI4ODgQxhCUSiWkUqlgvp6wtJZloV6vC/v5XVjD9gvbRkS06iJXJfkP03HfXyZRimJEtDqidEtdqqqkn/3sZwBOql78333vvPNO35QURES0mMYODO+88w6Ak8bV3oZlSZIYGIiIlsTYgeHVV18FcNLI+9Of/vTcMkRERLM1VmB46623gukvXnnllaEjna9evTqhbBERLa6+doQFMzIwvPnmm0gmk0G/fEVR+sYNACfdND/66KPzySUREU3NyMBw7do1YbDWpOYUIiJaVu+2Pxm67fLmQ1PMyelEHsfwgx/8AKqq4urVq3j99deDwWzVahXPPvvsxDNIRDSPPvzRj4Pfr74vdtV/d9qZmbDIgcEwDHznO98BcLLgTblcxnPPPYfnn38ev/vd7yaeQSKiedA7VqE3GIRZP37v3gt3JuujRRI5h37bwo0bN9But/Hiiy8CANsXiIjG4A9im2eRA4MsyygWizBNE/l8HsDJGg298xMREdFiihwYTNPE/v4+dF0P5kRyXRfFYvQh2kREi6y3kdmvMlr0CXNOvVBPt2vXroXOtkpERIsj8kopv/jFL/DEE0/g4sWLwc+FCxdw8eLF88gfERFNWeQSQy6XQ7VaRTqdPo/8EBHNFb9bapReSIsucmBIJBL4xje+cR55ISKaqUHTZ9/3m/85g5zMVuTAkM1msbOzg52dHciyLKx8xrmSiIgWX+TAUC6X0Wq1gpXTfJwriYhWgTBYbUlFDgz+ugxEREvr5j8Ev65CIOh1qrHZr7/+OkzTxObmJr73ve/hxo0biMVieO655yadPyIimrLI3VWLxSIMw0A+nw+qkxKJBHRdn3jmiIho+iKXGGq1WjBZnj9vkiRJcF137GM4joNarQZZluE4DvL5vNCIPW7asG21Wi3oUtt7bNu2AZysLeE4DlzXhaIoY+efiJbEGyXxtfv2bPIxZ0418vnXv/41nnrqKcRiMQAni/lsbGyMfYxsNhusEe04DnK5HEzTjJx21LZe5XIZhUIBhmGgWq0CANLp9NBzE9Hy+GH9t8HvdmcPAHC5c/L8eObQBQBcHbBf55zzNY8iB4b9/X0kk0mkUim022289NJLsCxr7Ier4zjCa1mWYVlW5LRh21zXhWmaUFU12F6pVFAoFACcVH21220A/aUJIlo9izDj6TSdanbVVquFn//853AcB5Ik4dVXXx27xGBZFuLxuPBePB6Hbdt91TlhaQ8PD4duk2VZCAq1Wk14DTAgEBENc6peSZ1OB5lMBuvr0ecQHNYWMWjJ0LC0Ydu6A4zrumi1WsK04K7rolarAQAODg6gadrAacOPj49xfHwcvO50VrFQSUSrJlKvpFdeeQWXLl3C5uYmNjc3cfHiRbz88ssTeWBGabwOS9u7Tdd1vPDCC8J7+XweqqpCVVXs7Owgk8kMPFapVMLGxkbwc+XKlbHzSES0qMYODE888QTefPNN7O/vo9VqodVq4eDgAB9++CESicTYwUGSpL7SQavVGli1E5Z2nOO4rgvLsvqO3d0+4fdo6m2zAE665h4dHQU/t27dGusaiWj+PPH3fxf8PHPjbTxz42089cs/4qlf/nHWWZs7Y1UlvfTSS1BVFa+++qrwvqIoME0TlUoFuq7jJz/5ychjpdNpGIbR934ymYyUVpblkcc5PDwc2FV1e3s7aHz29bZXAMDa2hrW1taGXgsRza/eCfHuO/7fM8rJ4hkrMFiWFYxdGKRQKODJJ58c64S9dfmO4yCZTAYPcNu2IUkSZFkOTTuoFNB9HP9YvQ98WZZRLpeFa1NVlY3RRIugd9xBmJ4xCas4tcVpjRUYxulxFGUcg2ma0HUdqVQKBwcHQlfXUqmEVCoVdC0NSxu2zdcbXCRJQjKZRKVSgSRJaDabHMdAtGD+0Rk9YeetC38SXvfXCdAwMc8fvhzCf/CeNc2i63Q62NjYwNHR0al6ZBHRGd0tMQwLDA/+r3uD2P5PrL+n46L4itcfxq598+l7L54tRj5mlOfXWCWGZrOJYjE8I4Mab4mIJiFoL7hbPeSXBv7N50/MKktLbazAsL29jWazOTINEdGkiFNYnJQQLnf+NCw5TdBYgYF18EQ0j/7rBXF9mKcWuPponpxq5DMR0aQF1UV3F8n5wL1XOrjclY7jDs4fAwMRzZ2r1nuIc2K7mWFgIKJT6R1A1u3lr7481jE+/NGPg9+vvn+3V6PL8QazxsBARGfTtT5yoH108u/dbpXdDcnd/t155YnOhIGBiCbiVlebQNH91ckv//HfDk3/7noC/+I3/SumrR+zCmnWGBiIaDwTXgbzcqeB9WM2JM8jBgYiCuW3JVyt/Q9hpbPusbm//tf/cuj+7EW0eBgYiOjM+PBfLgwMRHOgu3dOry/+zV+fyznFkcV7AIBnbpxUD13efCjYdvVczk7zjIGBaEKCKhfzV33bUl9KDX3Av/bWa/e6avbscx758/nTTHTzp6buvD/RU9OCYWAgWjLdAeAfm6OnpybqxcBAtAh6egS9drdH0NXDB/Bu+xNhW/cf9TMDDsUFa2gUBgaic9D7sP7sk4+Av70OAPiafEnYNqgaaZDuNQj8aafvaw/v888AQKfFwEAU0aiqmvt6gkIU3QHls09Ojv1O/bf42j+xSoimh4GBaAre61qI/j//Jtq+T/z93+HBI3FKCU4vTeeJgYFWWthEcIN6FwFi982zlA6I5hUDA9Gc6i5l/DNLCDRFMwkMjuOgVqtBlmU4joN8Pg9JkiKnDdtm2zYAQFEUOI4D13WhKErk8xNNCxuLaV7MJDBks1k0Gg0AJw/pXC43dPnQsLRh2wzDQLVaBQCk02nh+FHOT9Tbw8i3fvwe1ofs01l7dOg+RPNu6oHBcRzhtSzLsCwrctpRx0kkEmi32wAglAainJ/otBgAaJFdmPYJLctCPB4X3ovH40HVz7hpxzmOJEl9VURRzk9EtIqmXmJwXXfg+61Wf+NaWNpRx3FdF7VaDQBwcHAATdMgy3Kk8x8fH+P4+Dh43el0Bu5LRLRM5qZX0rAHdtS0/rbuBmVZlpHJZNBsNiMds1Qq4fr162Pni4hoGUy9KkmSpL5v561Wa2CvoLC0o47T3Zbg9z5yHCfS+YvFIo6OjoKfW7duRbhSWgruH04amQf8EC2rqQeGdDo98P1kMhkpbdg227axvb3dty0ej0c6/9raGtbX14UfWk7/2Pyo7+fd9ifCimVEq2LqVUmyLAuvHcdBMpkUxh9IkgRZlkPT9n7D794myzLK5XKwzbIsqKo6cj8iIppRG4NpmtB1HalUCgcHB8IYglKphFQqhUKhMDLtsG2SJCGZTKJSqUCSJDSbzbH2IyIiIOZ5njfrTCyKTqeDjY0NHB0dsVppSfhzJQ2aJfWZG2+zLYFm4itevO+9a998+t6LZ4uRjxnl+TX1NgYiIppvDAxERCSYm3EMRLPgT63N6bOJ7mGJgYiIBAwMREQkYGAgIiIBAwMREQkYGIiISMDAQEREAnZXpaXnj24e5Or0skG0MFhiICIiAQMDEREJGBiIiEjAwEBERAI2PhMBnF6bqAtLDEREJGCJgVZK74I8/qyqXHaJ6B6WGIiISMDAQEREgplUJTmOg1qtBlmW4TgO8vk8JEmKnDZsm23bsCwLAHBwcIDd3V1hGwAoigLHceC6LhRFOc9LJiJaGDMJDNlsFo1GA8DJwz2Xy8E0zchpw7ZZloVCoQAAqFQq2N7eDtIahoFqtQoASKfTQ89Ny8FfpQ3gSm1E45h6VZLjOMJrWZaDb/ZR0oZts20bpVIp2KaqKmzbDvZJJBJot9tot9uo1+tDSytERKto6oHBsizE43HhvXg8HlTvjJs2bJuiKNjd3Q3ed1032O6TJIkBgYhogKlXJfkP6V6tVitS2lHHUVU1eG9vbw/pdDoIBK7rolarAThpf9A0DbIs9x3r+PgYx8fHwetOpzPwnEREy2RuxjEMe9BHTdu7zQ8CfvsCAKGRWpZlZDIZNJvNvmOVSiVcv3597HwRES2DqVclSZLUVzpotVoDq3XC0o57HF3X+9oRutsn/B5NvW0WAFAsFnF0dBT83Lp1a8yrJCJaXFMvMaTTaRiG0fd+MpmMlFaW5ZHHqVQq0HUdsiwHJQnHcbC9vY12uy3s19teAQBra2tYW1sbeU00e1yMh2hypl5i6K3LdxwHyWRSGGPgf3sPSzvqOLVaDYqiBEFhf38/2K9cLgf7WZYFVVXZEE1EdNdM2hhM04Su60ilUjg4OBDGEZRKJaRSqWAMQljaYdscx0E2mxXOKUlS0LaQTCZRqVQgSRKazSbHMRARdYl5nufNOhOLotPpYGNjA0dHR1hf57Rr8yS0Kskf4Ob+AZ07f5lSjohO7ytef9X2tW8+fe/Fs8XIx4zy/OJcSUREJJib7qpEkzJ0au1jlhaIxsESAxERCRgYiIhIwKokWhp+IzNnUCU6G5YYiIhIwMBAREQCViXRQgkbrwAAcP/A3kdEZ8QSAxERCRgYiIhIwKokWki9g9iAk95IrEYiOjuWGIiISMDAQEREAlYl0ULxB7HF3+9fipWIJoMlBiIiEjAwEBGRgFVJNHeCQWw3/6Fv21X3vSnnhmj1MDDQTIwcwQzglvunvvfiXIGN6NwxMNC5GhYAguU2Abw7ZDbUp47/eC55IqJwDAw0c+vHrB4imiczCQyO46BWq0GWZTiOg3w+D0mSIqc9j210PoYttwkA4cuSE9G0zSQwZLNZNBoNACcP6VwuB9M0I6c9j210Pi53GsLrdVYTEc2tqQcGx3GE17Isw7KsyGnPYxuN9sP6b/veszt7APof/k/98o+I3/09DiJaFFMPDJZlIR4XHxPxeBy2bUNRlLHTHh4eTnxb7/mXxhulwe8/Wxx7n9fctwEAH7h/wrvrieD9Z268jWfu/s5SANFymHpgcF134PutVitS2vPY1uv4+BjHx8fB66OjIwBAp9MZeIxJ++NPDeH1/3379Uj7/y7WDk/wH/7L2Md68O6/TwJ4Ev9tYBqutEw0Gf/P+7Tvvc4/3+l6Ef0Z5D+3PM8bmXZueiUNe2BHTTvJbaVSCdevX+97/8qVKyG5IyI6B//pv3e9+PenPszHH3+MjY2N0DRTDwySJPV9O2+1WgN7BYWlPY9tvYrFIr71rW8Frz///HO0Wi1cunQJsVgMwEkUvnLlCm7duoX19dXqX7Oq176q1w2s7rUvw3V7noePP/4Yjz322Mi0Uw8M6XQahmH0vZ9MJiOllWV54tt6ra2tYW1tTXhvWLfW9fX1hf0Pc1areu2ret3A6l77ol/3qJKCb+qBQZZl4bXjOEgmk8ED17ZtSJIEWZZD0/Y+oCexjYiIZtTGYJomdF1HKpXCwcGBMIagVCohlUqhUCiMTHse24iIVp5HZ3Lnzh3vu9/9rnfnzp1ZZ2XqVvXaV/W6PW91r33VrjvmeWP0XSIiopXBhXqIiEjAwEBERIK5GeA272zbRi6XCybf8/nzLLmui4ODA+zs7ARTayzLDK/Drt227eD6Dw4OsLu7K/QuAwBFUeA4DlzXHetzmSenue5lv+dh25b5nodtW5Z7Lph1I8ciME3TazQa3qCPS5Ikr9FoeJ7neYZheLIsB9sURQl+bzabnqqqY22bJ2HXXi6Xhd+7rymfz3sAPABeOp322u12sG0Rrv20173s9zxs2zLf87Bty3DPezEwRDDoP0W9Xg9+Nwwj+I/QbDaF/xSedxJERm2bV73X3mg0hDw3m00PgNdsNj3PO/ks2u228HDw0y3StUe57mW/56O2Les9D9u2bPfcxzaGM0qn08HvpmlC0zQA4TPDhm1bFIqiYHd3N3jtzzXVfV2DBhQu+rWHXfey3/NxLOM9D7Os95xtDBNg2zb29vaQyWSQz+cBnH5m2EWiqmrw+97eHtLpdPBQcF0XtVoNwEk9vKZpkGV5Ka592HWvwj0Ps8z3fJhlvecMDBOgKApkWYau66jVasKDo9dpZ3+dZ/4DobtRrruRTZZlZDIZNJvN0GMsmkHXPSzdabYtmlW45+Na9Hu+soGhWq2G/qfNZDJCNdEokiQhm80ik8mg3W5PbIbX8zDpa9d1HfV6Xci/4zhBjxS/R4bjODO99vO+7lW654Oswj3vNc/3/Exm3cixSHo/rnq9PrAhstFoDG14arfbodvm1bD/KuVyOWhw9hseexto2+22B2Ahrz3Kda/KPR+0bRXu+aBty3bPfWx8jqi7GBiPx4VvG/7MsH7VUrfuWVxHzTA7r3qLwLVaLbhW13Wxv78fXF+5XA7SWZYFVVUX9tqjXHe3Zbznw7Yt+z0ftm0Z7zmwwlVJUViWhXq9DuDe7K+qqkJRFOzs7KBarQIA6vW6UN+8DDO8Drt2x3GQzWaFtJIkBfXMyWQSlUoFkiSh2Wwu3LWf5rqB5b7nYduW+Z6P2rYM97wXJ9EjIiIBq5KIiEjAwEBERAIGBiIiEjAwEBGRgIGBiIgEDAxERCRgYCCakUqlglgshmw2i0qlAl3Xkc1mF2L2TVpuHMdANCOu62JzcxPdf4K2bSORSIB/ljRLLDEQzYhlWcGkc70WYQZOWl4MDEQzUq/XkUwmhfdKpdLirAtMS4tVSUQzsrW1BVVVkclk4DgOXNeFoihnnv6a6KwYGIhmwG9f8Nfu8Cfnu3HjBksLNHOsSiKagcPDQ8iyLKx4Fo/Hsb+/L6RzXRe6rrOnEk0VAwPRDNTr9b4qo8PDw750h4eHbIimqWNgIJoBy7KQyWSE91zXDRZ38UsI6XSaVUs0dQwMRFPkOA4qlQpkWcbBwQEcxwm2FQoF1Ot1VKvVod1YiaaBjc9Ec07Xdezs7DBY0NSwxEBERAKu+Uw0xyzLEnoksdRA08CqJCIiErAqiYiIBAwMREQkYGAgIiIBAwMREQkYGIiISMDAQEREAgYGIiISMDAQEZGAgYGIiAQMDEREJPj//ah+ImLafYAAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.rcParams.update(\n",
    "{\n",
    "    'text.usetex': True,\n",
    "    'text.latex.preamble': r'\\usepackage{slashed}',\n",
    "})\n",
    "# plt.subplot(2, 2, 1)\n",
    "plt.figure(figsize=(4, 3)) \n",
    "\n",
    "# second last index is for the scale, second for type and first the PCD\n",
    "# 1 => Scale=50\n",
    "plt.hist(R_pcds[3][0][1][1], bins=40,density=True,  alpha=0.5)\n",
    "plt.hist(R_pcds[3][1][1][1], bins=40,density=True,  alpha=0.5)\n",
    "plt.hist(R_pcds[3][2][1][1], bins=40,density=True, alpha=0.5)\n",
    "plt.hist(R_pcds[3][3][1][1], bins=40,density=True, alpha=0.5)\n",
    "\n",
    "plt.legend(['BP1','BP2','BP3','SM'])\n",
    "\n",
    "plt.title(r'Scale $=50$ GeV, MET $>150$ GeV')\n",
    "plt.xlabel(r'$R_1$')\n",
    "plt.ylabel(r'Density')\n",
    "plt.tight_layout()\n",
    "\n",
    "# plt.savefig(f'plots/hist_150_R1_Scale_50GeV_MET150.pdf',bbox_inches='tight')\n",
    "\n",
    "\n",
    "plt.rcParams.update(\n",
    "{\n",
    "    'text.usetex': True,\n",
    "    'text.latex.preamble': r'\\usepackage{slashed}',\n",
    "})\n",
    "# plt.subplot(2, 2, 1)\n",
    "plt.figure(figsize=(4, 3)) \n",
    "\n",
    "plt.hist(R_pcds[3][0][5][1], bins=40,density=True,  alpha=0.5)\n",
    "plt.hist(R_pcds[3][1][5][1], bins=40,density=True,  alpha=0.5)\n",
    "plt.hist(R_pcds[3][2][5][1], bins=40,density=True, alpha=0.5)\n",
    "plt.hist(R_pcds[3][3][5][1], bins=40,density=True, alpha=0.5)\n",
    "\n",
    "plt.legend(['BP1','BP2','BP3','SM'])\n",
    "\n",
    "plt.title(r'Scale $=250$ GeV, MET $>150$ GeV')\n",
    "plt.xlabel(r'$R_1$')\n",
    "plt.ylabel(r'Density')\n",
    "plt.tight_layout()\n",
    "\n",
    "# plt.savefig(f'plots/hist_150_R1_Scale_250GeV_MET150.pdf',bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from imblearn.over_sampling import SMOTE\n",
    "from imblearn.combine import SMOTEENN\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.metrics import accuracy_score\n",
    "from sklearn.metrics import confusion_matrix, recall_score\n",
    "from imblearn.under_sampling import EditedNearestNeighbours\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.model_selection import cross_val_predict\n",
    "from sklearn.metrics import precision_recall_fscore_support\n",
    "\n",
    "\n",
    "def get_cl_scores(vector1,vector2,subset_size=None):\n",
    "\n",
    "    if(subset_size):\n",
    "        if(len(vector1) >=subset_size and len(vector2) >=subset_size ):\n",
    "            indices1 = np.random.choice(len(vector1), size=subset_size, replace=False)\n",
    "            vector1 = np.array(vector1)[indices1]\n",
    "\n",
    "            indices2 = np.random.choice(len(vector2), size=subset_size, replace=False)\n",
    "            vector2 = np.array(vector2)[indices2]\n",
    "        \n",
    "        else:\n",
    "            vector1 = np.array(vector1)\n",
    "            vector2 =np.array(vector2)\n",
    "\n",
    "    else:\n",
    "        vector1 = np.array(vector1)\n",
    "        vector2 =np.array(vector2)\n",
    "\n",
    "    vectors = np.hstack((vector1,  vector2)).reshape(-1, 1)\n",
    "\n",
    "    labels1=np.full(len(vector1),1)\n",
    "    labels2=np.full(len(vector2),2)\n",
    "    labels = np.hstack((labels1,labels2)).reshape(-1, 1)\n",
    "    \n",
    "    smote = SMOTE(random_state=42,k_neighbors=15) # SMOTEENN(enn=EditedNearestNeighbours(sampling_strategy='all'),random_state=42)  #SMOTE(random_state=42,k_neighbors=5)\n",
    "    vectors_res, labels_res = smote.fit_resample(vectors, labels)\n",
    "    \n",
    "    print(labels_res,np.sum(labels_res))\n",
    "    scaler = StandardScaler()\n",
    "    X_scaled = scaler.fit_transform(vectors_res)\n",
    "\n",
    "    # print(vectors.shape,labels.shape)\n",
    "    df_vectors=pd.DataFrame(vectors)\n",
    "    df_labels=pd.DataFrame(labels)\n",
    "   \n",
    "    # Choose a Model (Support Vector Machine)\n",
    "    model = SVC(kernel='rbf', C=0.5, random_state=42)\n",
    "\n",
    "    scores = cross_val_score(model,X_scaled,labels_res, cv=10, scoring='accuracy',n_jobs=-1)\n",
    "\n",
    "    return np.mean(scores)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "def get_acc_pcds(idx):\n",
    "    \n",
    "    BP1_res=BP1_pcds[idx]\n",
    "    BP2_res=BP2_pcds[idx]\n",
    "    BP3_res=BP3_pcds[idx]\n",
    "    BG1_res=SM_pcds[idx]\n",
    "    \n",
    "    R_BP1=R_pcds[idx][0]\n",
    "    R_BP2=R_pcds[idx][1]\n",
    "    R_BP3=R_pcds[idx][2]\n",
    "    R_SM=R_pcds[idx][3]\n",
    "\n",
    "    # vector1 = BP1_res[['px','py','pz','E','mt']].to_numpy()\n",
    "    # vector2 = BP2_res[['px','py','pz','E','mt']].to_numpy()\n",
    "    # vector3 = BP3_res[['px','py','pz','E','mt']].to_numpy()\n",
    "    # vector4 = BG1_res[['px','py','pz','E','mt']].to_numpy()\n",
    "\n",
    "    # normal_score_bp1=get_cl_scores(vector1,vector4,20000)\n",
    "    # normal_score_bp2=get_cl_scores(vector2,vector4,20000)\n",
    "    # normal_score_bp2=get_cl_scores(vector3,vector4,20000)\n",
    "\n",
    "    acc_R_BP1=[]\n",
    "    scales=[25,50,100,150,200,250,300,350,400,450,500]\n",
    "    for i in range(len(R_BP1)):\n",
    "        if(i==0): continue\n",
    "        print(f\"BP1, Filtration scale: {scales[i]} \\n\")\n",
    "\n",
    "        vector1 = R_BP1[i][1].to_numpy()\n",
    "        vector2 = R_SM[i][1].to_numpy()\n",
    "        print(vector1.shape,vector2.shape)\n",
    "        acc_R_BP1.append(get_cl_scores(vector1,vector2,2000))\n",
    "        \n",
    "        print(\"-------------------------------------------\\n\")\n",
    "\n",
    "\n",
    "    acc_R_BP2=[]\n",
    "    for i in range(len(R_BP2)):\n",
    "        print(f\"BP2, Filtration scale: {scales[i]} \\n\")\n",
    "\n",
    "        if(i==0): continue\n",
    "        \n",
    "        vector1 = R_BP2[i][1].to_numpy()\n",
    "        vector2 = R_SM[i][1].to_numpy()\n",
    "        acc_R_BP2.append(get_cl_scores(vector1,vector2,2000))\n",
    "        \n",
    "        print(\"-------------------------------------------\\n\")\n",
    "\n",
    "\n",
    "    acc_R_BP3=[]\n",
    "    for i in range(len(R_BP3)):\n",
    "        print(f\"BP3, Filtration scale: {scales[i]} \\n\")\n",
    "        \n",
    "        if(i==0): continue\n",
    "\n",
    "        vector1 = R_BP3[i][1].to_numpy()\n",
    "        vector2 = R_SM[i][1].to_numpy()\n",
    "        acc_R_BP3.append(get_cl_scores(vector1,vector2,2000))\n",
    "        \n",
    "        print(\"-------------------------------------------\\n\")\n",
    "\n",
    "\n",
    "    return [acc_R_BP1, acc_R_BP2,acc_R_BP3]\n",
    "   \n",
    "    \n",
    "acc_pcds=[]\n",
    "for idx in range(n_pcds):\n",
    "    print(f\"********************************** PCD ID: {idx} ********************************** \\n\")\n",
    "    acc_pcds.append(get_acc_pcds(idx)) \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>pcd_no</th>\n",
       "      <th>scale</th>\n",
       "      <th>BP1</th>\n",
       "      <th>BP2</th>\n",
       "      <th>BP3</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>0</td>\n",
       "      <td>50</td>\n",
       "      <td>0.99700</td>\n",
       "      <td>0.69425</td>\n",
       "      <td>0.73125</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>0</td>\n",
       "      <td>100</td>\n",
       "      <td>0.60275</td>\n",
       "      <td>0.68800</td>\n",
       "      <td>0.60450</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>0</td>\n",
       "      <td>150</td>\n",
       "      <td>0.62825</td>\n",
       "      <td>0.58700</td>\n",
       "      <td>0.61050</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0</td>\n",
       "      <td>200</td>\n",
       "      <td>0.58550</td>\n",
       "      <td>0.52275</td>\n",
       "      <td>0.48525</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0</td>\n",
       "      <td>250</td>\n",
       "      <td>0.51200</td>\n",
       "      <td>0.56300</td>\n",
       "      <td>0.50675</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>195</th>\n",
       "      <td>19</td>\n",
       "      <td>300</td>\n",
       "      <td>0.53650</td>\n",
       "      <td>0.49550</td>\n",
       "      <td>0.51875</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>196</th>\n",
       "      <td>19</td>\n",
       "      <td>350</td>\n",
       "      <td>0.51650</td>\n",
       "      <td>0.51350</td>\n",
       "      <td>0.50675</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>197</th>\n",
       "      <td>19</td>\n",
       "      <td>400</td>\n",
       "      <td>0.48550</td>\n",
       "      <td>0.54050</td>\n",
       "      <td>0.51675</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>198</th>\n",
       "      <td>19</td>\n",
       "      <td>450</td>\n",
       "      <td>0.51225</td>\n",
       "      <td>0.52550</td>\n",
       "      <td>0.50025</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>199</th>\n",
       "      <td>19</td>\n",
       "      <td>500</td>\n",
       "      <td>0.49400</td>\n",
       "      <td>0.52775</td>\n",
       "      <td>0.50525</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>200 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     pcd_no  scale      BP1      BP2      BP3\n",
       "0         0     50  0.99700  0.69425  0.73125\n",
       "1         0    100  0.60275  0.68800  0.60450\n",
       "2         0    150  0.62825  0.58700  0.61050\n",
       "3         0    200  0.58550  0.52275  0.48525\n",
       "4         0    250  0.51200  0.56300  0.50675\n",
       "..      ...    ...      ...      ...      ...\n",
       "195      19    300  0.53650  0.49550  0.51875\n",
       "196      19    350  0.51650  0.51350  0.50675\n",
       "197      19    400  0.48550  0.54050  0.51675\n",
       "198      19    450  0.51225  0.52550  0.50025\n",
       "199      19    500  0.49400  0.52775  0.50525\n",
       "\n",
       "[200 rows x 5 columns]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bp1_acc=[]  \n",
    "bp2_acc=[] \n",
    "bp3_acc=[] \n",
    "acc=[]\n",
    "scales=[50,100,150,200,250,300,350,400,450,500] \n",
    "for i in range(n_pcds):\n",
    "    temp1=[]\n",
    "    temp2=[]\n",
    "    temp3=[]\n",
    "\n",
    "    for j in range(0,10): # scales\n",
    "        temp1.append(acc_pcds[i][0][j])\n",
    "        temp2.append(acc_pcds[i][1][j])\n",
    "        temp3.append(acc_pcds[i][2][j])\n",
    "\n",
    "        temp=[i,scales[j],acc_pcds[i][0][j],acc_pcds[i][1][j],acc_pcds[i][2][j]]\n",
    "        acc.append(temp)\n",
    "\n",
    "    bp1_acc.append(temp1)\n",
    "    bp2_acc.append(temp2)\n",
    "    bp3_acc.append(temp3)\n",
    "    \n",
    "df_acc=pd.DataFrame(acc,columns=['pcd_no','scale','BP1','BP2','BP3'] )\n",
    "# df_acc['BP2']*=0.95\n",
    "\n",
    "df_acc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEiCAYAAAD9DXUdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB1q0lEQVR4nO39d3xc933nC79Pm94HAMEClgGpLksCQbnGjYC85e6zGxuQfZ3d3BKTiJ+UTRyJEDd7YyuNAmNvsonzRICcbF67z3VCAZI3z7N399qEZMeWm0hClK1GkRiADUSdGUyv59w/BjPEEANgQKKR+r31GmHmzG/O+c4BeD7n9/s2yTAMA4FAIBAI5pA32gCBQCAQbC6EMAgEAoGgAiEMAoFAIKhACINAIBAIKhDCIBAIBIIKhDAIBAKBoAIhDAKBQCCoQAiDQCAQCCoQwiBYF/r6+pAkic7OzkXHdHd3I0kS3d3d5W3Hjx9HkqSqj/b29vK+vV7vso/jx4+v+ffcaLq6uvB6vUiSRHNzM11dXQSDwZo/txRer5eBgYHVMlWwiVE32gDBewePx7PkhaWvrw+Px1OxbWZmBo/Hw8jIyKKfO3z4MG1tbeXXwWCQ9vZ2+vv7aWlpKW/3+Xw3b3yNRCIRnn/+efr7+zl58mTFe4ODg/T29tLe3k4gEODkyZMcOHCAjo6O8pjjx4+Xz0EkEuHIkSM1Hbf0nQOBAP39/bS2tnL69Gl6e3sZGBhYdj9dXV309fUxODhYcS7n2w5U2Cq4cxHCIFg3fD4fgUCAgYGBBReYwcFBfD7fAmEosdj2EoFAoOq2atvXiqGhIU6fPk0kEiEUCi14PxKJMDg4yMDAAIFAgO7u7gWiAEWhg+I56erqore3d9ljd3Z2lsWmRFtbW9WLfDVaWlrKolLtM729vTz++OM17Utw+yOWkgTrymc/+9mqF7r+/v4NvRsdGBioedllMVpaWjh8+PCSYjQyMoJhGAwPD5cFoMSxY8cqtrW1tdHX11eT7UNDQzUJyFJ0dHRUPV4kEimfH8F7AyEMgnWlo6ODwcFBIpFIxfa+vr4NvfB0dHTQ09NDb28vXV1dDA0Nrevxg8EgkUik6syotIyzGCdOnCjf8ddCyZ/Q3NxcIQSl83/jct/zzz9PIBCoWJYT3NkIYRCsK6ULzPPPP1/eVlpaWezCFolEqjqfl7tgrhSPx1MWh8HBQTo7O1f9GM8//zwDAwP09fVVONkXm6l4PJ4FInojwWCQ1tbWmo7f2dlJMBhkZGSEkydP0t3dXRbB0u/mxIkTFZ8piaXgvYMQBsG6c+O6+YkTJ5a88Hg8HgzDWPCodf38Zjhy5Aj9/f0Eg0E6OztXJRqnpaWFtrY2Ojo6OHz4MM3NzUtGaUHRL1PNXzGfxWYaNxIMBhkYGKC/vx+Px0MgEKCnp6dCCLq6uhgYGCiLUSQSYWhoSDid32MIYRCsO48//jhDQ0Plu+RqzujNwuHDh+nv7wegubl52bv3pbhxVvT4449XXISrsZwolPZ749JXZ2fngrDe0pg9e/aUQ3jnzxhKNgHlGV1fX9+KlqkEdwZCGATrjsfjoa2trbykspkvPJFIhO7ubk6cOEFvb29Nd+aLceOso7SvYDC45DLacuemvb19gd+mv78fwzAWhKm2tLQQDocrHvMjmTweDx0dHeUZnVhGem8ihEGwIXR1dXHixAn6+/s35YWnJAjd3d10dXUtGsa5kv2V1vfnb4PrMwmPx1PV17DccY8cOVIOf12KlpYWhoaGlp31lJzvg4ODBIPBBdFTgjsfkccg2BA6OjrK6+ulpZqlWOxidit38NUIBoP09PQAxUzsm5nJVFv+8Xg85Qt4ib6+Pjo6Osrf4ejRowwODpYvxAMDAzVflPv7+zl48CChUIijR48SCAQIBoMVQhMIBDh8+DCdnZ309vaWc0qCwWDFzKKtrQ2Px0NXV9emXeITrC1CGAQbRukCvNzFPRKJLFquYbValpcEoRSZdDOCU3LunjhxgqGhIbq7uysym48ePVpRlmNmZqZCFI8cOcLx48fLS06nTp2qOTehpaWFkZERuru7OXToEENDQ+Uoo/kzid7eXrq7u9m/fz8Ara2t5d/DfB5//HH6+vpuOTdCcHsiGav1L0sguI1Zap1fIHivIYRBIBAIBBUI57NAIBAIKhDCIBAIBIIKhDAIBAKBoAIhDAKBQCCoQAiDQCAQCCoQwiAQCASCCu7YBDdd1xkbG8PpdCJJ0kabIxAIBGuCYRjEYjG2bduGLK/Ovf4dKwxjY2M0NTVttBkCgUCwLly+fJkdO3asyr7uWGFwOp1AsZXiejSBX4pcLsd3vvMdHnvsMTRNE7YIW4QtwpZVIxQKsWfPnvI1bzW4Y4WhtHzkdDpxuVwbaksul8Nms+FyuTb8j0jYImwRttx5tgCrumQunM8CgUAgqEAIg0AgEAgqEMIgEAgEggruWB+DQCC4PdB1nWw2u+y4XC6Hqqqk02kKhcI6WLY5bNE0DUVR1vQYNyKEQSAQbBjZbJaRkRF0XV92rGEYNDY2cvny5Q3PTVpvWzweD42Njev2vYUwCASCDcEwDK5du4aiKDQ1NS2bnKXrOvF4HIfDsWqJXDfLetliGAbJZJLJyUkAtm7dumbHmo8QBoFAsCHk83mSySTbtm3DZrMtO7605GSxWDaFMKyXLVarFYDJyUkaGhrWZVlJOJ8FAsGGUFqbN5lMG2zJ5qcknKWchbVGCINAINhQNtpfcDuw3udICINAIBDUSGdnJ16vlwMHDvDII4+wb98+uru7q47Zv38/+/fvp7m5ecGYzY4QBoFAIFgBR48e5dSpU7z22mucOnWKSCRCe3t7xZjDhw9z5swZzpw5w/DwMMFgkK6urooxQ0NDtLe3E4lE1tH62tgQ5/PQ0BCHDh3izJkzS44LBoMMDAwQCAQIBoMcPnwYj8ezPkYKBALBMng8Hnp6evB6vUuO6+npYf/+/fT29gLQ1dVFKBRicHBwPcxcMesuDKUL/dDQ0LJjOzs7y+IRDAY5dOgQ/f39a22iQCDYAAzDIJVbPFlM13VS2QJqNr/qkUBWTbmpdfxIJEJXVxdHjhxZdtx8SgKxWf0r6y4MHR0dNY0LBoMVrwOBwKZVV4FAcOukcgXu+71vb8ix3/r9T2Ez1XY5PHbsGL29vczMzDA7O0tvby+HDx9edPzQ0BCdnZ0cPXp0tcxdczatj2FwcHBBHwWfz1fTTEMgEAjWiqNHj3L+/HlGR0c5deoU/f39C/wHfX19Zedzd3c3PT09y84qNhObNsFtMYdMKBSquj2TyZDJZMqvo9EoUIz7Xa/Y38UoHX+j7Zhvg7ClEmFLddbSllwuh2EY6LqOruuYFYk3vtK+6HjDMIjH4jicjlVfgjErUs1lOUoPgEceeYRvf/vb+P1+Dh06REtLC4ZhcOjQIZ555pmKzy62/9L3Xwpd1zEMg1wutyDBbS1+N5tWGBZjMcE4duwYTz/99ILt3/3ud2vKqlwPTp48udEmlBG2VEfYUp21sEVVVRobG4nH4zUV0QOwmhQKmdSq2xJL1zYun8+TTqeJxWLFz8393LVrF2+++SZ79+4ln8+TyWTKN6fLHjsWW9Znks1mSaVSfP/73yefz1e8l0wmazN+BWxaYfB4PAtmB6FQaNGopKNHj/KlL32p/DoajdLU1MQnPvEJ/H7/Wpq6LLlcjpMnT9Le3r4puj0JW4Qtm8GWdDrN5cuXcTgcWCyWZceXmt47nc4Nc9qqqorFYsHpdJZtee6557h48SL/4l/8C1wuF6qqYjaba+4cWUuXyXQ6jdVq5aMf/eiCczUzM3PT32cxNq0wtLW1lT3382ltba063mw2YzabF2zXNG3D/3GVELZUR9hSnTvdlkKhgCRJyLJcU5RRabml9JmNQJIknnnmGZ5//nkKhQKxWIyWlhbOnDlT9olKkrSsjd3d3eVgmgMHDtDS0rJkxKUsy0iSVPX3sBZ/IxsqDJFIpGIGMDQ0hMfjIRAIEAgEKsYGg0FaW1tFHoNAINgwShdvXdeJRqO4XK4FAlBLSH1PT8+a2LdarLvsDg4OltPDjx07xsDAQPm9G1/39/fT3d3NwMAAvb29IodBIBAI1oF1nzG0tbXR1tZWVTFvvPAHAoHyuFrzHwQCgUBwa2zaPAaBQCAQbAxCGAQCgUBQgRAGgUAgEFSwYmHw+Xx88Ytf5OzZs2tgjkAgEAg2mhULw8jICC0tLXzhC19g3759HD16lNHR0TUwTSAQCAQbwYqFwe12c+jQIU6fPs3p06cJBAJ0dnZy4MABvvrVr9acBi4QCASCzclN+xjOnj3LU089RU9PD3v27OGZZ55hz549dHR08LWvfW01bRQIBIJNgWjtuQh/8id/wt69e+ns7KS5uZkzZ87w/PPPc/DgQT7zmc/wne98h7/7u79bC1sFAoFgw1mN1p6Dg4M0Nzfj9Xo3ZXvPFQvD8PAw/f39nD9/nieeeAK3271gzI21yQUCgeBOpNTac7kmYj09PTz//PNAsRRQT08Pw8PDhMNhWlpaOHjw4HqYWzMrznxub29fUNnwhRdeQJIkPv3pTwNw6NCh1bFOIBAINjE309ozGAxWLC319PQgSRLBYHBBjbiNYsXC8NRTT5X7MJdoa2ujtbW1LAwCgUCwYgwDckv0WjD04vs5FaRVTsHSrFBjKe9bbe3Z0tJS8X6pjfFmEQW4CWGYmZlZUDvc7XavSU1wgUDwHiKXgj9/ZNG3JcBl6EirLQoAv/kamGpr6HX06FGeeOIJotEoFy5c4OjRo5w5c6aiTUBfX195ecnn89HT07NovbfN2PZzxcLQ2trKt771LX7xF3+xvO3ll1/eVGonEAgE60FLSwsnT57E6/XS1dVVng0cPny4ptLax48fL/spNhMrFoZnn32W1tZWjhw5QltbG8PDw2Xvu0AgENw0mrV4574IhqETjcZwuZyrP2vQrLf08UAgQDAYXLBMtBRdXV3s379/080W4CaEIRAIEAqFGBgYYGRkhLa2Nj7zmc+shW0CgeC9hCQtvZyj66DlQbPBBnVwq0ZfXx/BYJC2traaP9PV1UV7e/umbSdw0/0YNusXEggEgrXk2LFjnDhxYkFrz1q7Sw4ODtLX10dfX1/F9jNnzqxoxrGW3JQwvPzyywwNDVU4nP1+P0888cSqGSYQCASbjdVo7dnW1oZhGGtm42qw4vnYU089RUdHBxcuXKC/v59wOMzJkyc5derUWtgnEAgEgnVmxTOGF154gVAoBMCv/uqvcvz4cVwuF5/97GdX3TiBQCAQrD8rnjHMnwI1Nzfz0ksvASybEi4QCASC24MVC0NHRwcvvvgiUIzVffLJJzlw4AB79uxZdeMEAoFAsP6seCnpmWeeKT93u92cPHmSYDC46YpACQQCgeDmWPGM4atf/WrF6z179ghREAgEgjuIFQvDq6++yre+9a21sEUgEAgEm4AVLyUdPXqUxx9/vLx8ND+pY/fu3atomkAgEAg2ghULw6FDhzAMg7/6q7/ir/7qr8rbJUni/Pnzq2qcQCAQCNafFQvD6dOn18IOgUAg2PR0dnYyODhYrhknyzIdHR0V1VHnj4Fik54bx2x2NqQSVTAY5Pjx4wwMDHD8+PEl+52Wuh319fXR3d296XqjCgSC9xar0fN5aGiI9vZ2mpubaW9vLzfr2SyseMYgy/KC1p5QrLpa61JSZ2dnuUx3MBjk0KFDi9YXaW9vLxeoGhoaoru7u6IhhkAgEGwUpV4KXq93yXE9PT3s37+f3t5eIpEIhw4d4qWXXsLj8XD8+PGKa+JmYMXCEA6HK14bhkFfX9+CQlKLcaMyBgKBRbOmS9tLDu6Wlhb6+vqEMAgEdyCGYZDX84u+r+s6OT1HrpBDNlZ3sUOV1ao3vMtxMz2fPR5PWRSgeA3cbCshKxYGt9u9YNuRI0f41Kc+VVN11cHBQXw+X8U2n8/H0NDQgpKzi52samMXYzaZxe+vaahAINhA8nqe537+3KLvG4ZBJpPBbDbf1EV8KQ49eAhN0Woae6s9n+H6zW4wGOTYsWN0d3ffkv2rzYqFYXR0dMG2YDBY8xrZYhf7UmG++bS0tFTsd2hoaNGxmUyGTCZTfh2NRgF4dzyCx+vGaa7tl74W5HK5ip8bibClOsKW6qylLblcDsMw0HW9/FiqHHXpvbUoWa3rOrqkLzvOMAyeeuopnnjiCWKxGOfPn+fo0aOcPn2aZ599tjxmfs9nr9fLsWPH6OjoQNevH+NXf/VXee655zh48CBf+MIXKt6rZp9hGORyORRFqXhvLX43KxaGlpYWJEmq+OX4fL4Kx8rNUE0wAoEAPT099PX1lXMnSse7kWPHjvH0008v2B6+cJYfjL17S7atFidPntxoE8oIW6ojbKnOWtiiqiqNjY3E43Gy2SyGYfD4rsdX/Ti1kIwna5qF5PN50uk0sVgMgH379jEwMMDu3bv5pV/6JR566CHy+Ty//Mu/vOB6VLpZLXH8+HGOHz/O3/7t39LS0sL3vve9RY+bzWZJpVJ8//vfJ5+vXG5LJpM1fsvaWbEwVLtbXwkej2fBPkKh0KLdj44cOVKekZRa55XCwOZz9OhRvvSlL5VfR6NRmpqa2G7Pkmt8AF/DVu7a4sSkrn8gVi6X4+TJk7S3t6NpGzdzEbYIWzaTLel0msuXL+NwOLBYLMuONwyDWCyG0+lc9aWkWlFVFYvFgtPprLClubmZiYkJXC4XqqpiNptxuVw17fM3f/M3+e3f/m1effXVRduDptNprFYrH/3oRxecq/kN01aLFQvDiy++SCAQ4OGHHy5ve+GFF5AkiU9/+tPLfr6tra2q87i1tbXq+GAwWBaCkm+hmoiYzWbMZvOC7aZCgq36ZUbDVixmjbu2uFDkjfmj0jRtw/+hlxC2VEfYUp21sKVQKCBJErIs1xS8UlpqKX1mI5Akqfwovf7GN75BMBjkscceK0dtLmXj4OAgwWCw7Jfo6+vD4/Hw6KOPLvqZ0n6r/R7W4m9kxcLQ3d29IKzq4MGDHDhwoCZhuPFuPxgM0traWr7YDw0N4fF4yuP279/PyMgIHo+H3t7eFSeJ5E0utHycHcolLk2aMCkKgXr7ht1xCASC25tb7fnc1tZGX18fzc3NwMIopc3AioVhZmZmwRTJ4/GsaDrT399Pd3c3Bw4c4NSpUxU5DMeOHePAgQPl8K+enh4GBwcJhUJ0dnYuOtVaHImcvQFrfJytkongpIJJlWny2Va4H4FA8F5nNXo+QzEBbqlIpo1mxcLQ2trKt771LX7xF3+xvO3ll1+uuu6/GCWnMhQb/8znxpO6KidPUsjb6nGlxsjJFs5PqphVmQbX8uuaAoFA8F5jxcLw7LPP0traypEjR2hra2N4eLic+r2ZMRQTBbMLX+oiWcnCuQkZkyrjsZk22jSBQCDYVKxYGErFo1544YVypNBnPvOZtbBt1dE1O1Ihy5bMCJclE++Myzyw3Y3DvOLTIBAIBHcsNx2VNF8MVhKVtNEULF7UxATbcxcZiWmcG5e4f5sbi6Ys/2GBQCB4D7DimK/u7u4F/oSDBw9uupTupcjb6lHTIXYWLjMdSfDuRIxcYfmsR4FAsPqsRSbzncZ6n6MVC8NqRCVtOJJMzl6POXmNHcY1xsJJglMJdF38gQoE60WptEM2m91gSzY/pezm9cpr2ZCopE2BrJG3+LAnRtnqsjAyLWFSJHbXiRwHgWA9UFUVm83G1NQUmqYtm7Sm6zrZbJZ0Or1hCW7rbYthGCSTSSYnJ/F4PAvqJK0V75mopGoYqoWCbscVD5J3WxiekjFpCts91o02TSC445Ekia1btzIyMsLFixeXHW8YBqlUCqvVuuE3b+tti8fjobGxcc2PU+Kmo5IGBgYYGRm5raKSqqGbXMjJKTzxYbKOezk/LmFSZOqdC8trCASC1cVkMrFv376alpNyuRzf//73+ehHP7rhpULW0xZN09ZtplDipuM05yemvfzyy/T29nLixIlVMWq9yVvr0BLj1KdHuWpu5txEDJMi47Ztjjo1AsGdjCzLNRXRUxSFfD6PxWLZcGHYTLasBTe9ODY6OsrRo0fx+/03UaZikyFJ5GwNqIlxtuavkk5nOTcRJZldvJuUQCAQ3KmsaMYQjUbLrTWHh4eRJInnn3/+tl5KKiMXy2aYYpfY4bFwMbmFd8dj3LvNhVkVOQ4CgeC9Q00zhm984xt86lOfwuPx8Pd///ccOXKEUCjEI488cmeIwhylshnm2WG2K1HGo2nOT8TJixwHgUDwHqKmGcPhw4dpb28nGAyye/fu8vaNjgxYC0plM2yRd9nqe5ArYTCrMs31DuQN6uMgEAgE60lNM4YzZ87wyCOP0NLSwqc+9Sm+9a1vrbVdq8ZsZvkxN1KweJEKGRyz56k364xMJ7gcXv32eQKBQLAZqUkYHnnkEZ555hlCoRBHjhzh7/7u7/D5fAwNDfHd7353rW28JXrOqiSyK89oztvqUdIhXPEgTpPE+ck447PpNbBQIBAINhcrjko6ePAgzz//PKFQiBMnTnDs2DH8fj9f/OIX18K+W2Y0JvN//CBJOr9CcZgrm6ElxvClL2GWJc6NR5mJ38QURCAQCG4jbimXu6Ojg+985zsMDw9v2pIYVtXg51MF/uCHKfIrrYUka+QtXszREeqNaQo6nBuPEU3n1sZYgUAg2ASsSpEPj8fDk08+uRq7WnWefDiPWYFXr+U5/pMU+gqrFBqqlYJqwxx+ly1KjEQ2z7nxGKlsYY0sFggEgo1lYytRrQP3eAx+78NWFAm+eynP18+kV1zCVje7wdCxRN6l0ZxnJp7h3YkY2bwIYxUIBHced7wwADy6TeOpD1qRgP//hRx/+/OV+wny1jqUbAzr7AUabRLXZlNcmIpREKW6BQLBHcZ7QhgAPr5T4zdbi/VYvvlWlv53VigO88pm2GKj1NtNXJpJMjIdF41GBALBHcV7RhgA/qe9Jn7lfcWqqX1nM/yP4RU2CJEV8rY6TLFLOJJX8NnMjEwnuBJOrYG1AoFAsDG8p4QB4HP3mXn8HhMAf3Y6zQ8uryzCyFDMFMwuTJELOHNT2E0q5ydjTEZFjoNAILgzWLEwvPjii+zduxdFUcoPWZbXvV74rfCFh8z804CGbsCxH6c4M76yKqq6ZsdQNMzhd3FLCVRJ5p3xGOGEaFEoEAhuf1YsDIcOHaKnp4dQKFR+hMNhQqHQWth361RZ/5ckiX/bauGjTSo5Hb7ygyRvTa9MHIplM9JYQu/iM+vk8jrnJmLEM6JUt0AguL1ZsTDs37+fz3zmM7jd7gWPzcjWkQHQF16sFVniqQ9Y2d+okC7Av/9+kpHIynIT8rYGlHQIc/hd6u0Ks6kc58ajpHMix0EgENy+rFgYOjs7+exnP8uLL77I2bNnGR0dLT82I56pM2wZ+lOk/EIfgKZIfPkjNu6rU4hl4anvJRmLrSA3oVQ2I34Vc3SULQ4zU9FijkNOlOoWCAS3KStu7VlaRjp58mTFdkmSmJmZqWkfwWCQgYEBAoEAwWCQw4cP4/F4Fh07ODiIz+cjGAzS0dGxovIbuqxhm/45W1/9Yyb2/w4Fc+XMxqpK/OFHbTzxcoJgRKf7ewn+tM1OnbVGzZQ18lYf5mgQQ7PR4NrKWCSFSZG5a4tTlOoWCAS3HSsWhgsXLtzyQTs7Ozlz5gxQvPAfOnSI/v7+qmMHBgY4cuRI+XVXVxe9vb01H+vSvYdwXfkm5ugoW3/y+4y3Pkne3lgxxmmSOPYxG7/9UoKxuMHR7yX52idtuMy1iYOhWikUspjD76IrFuocHi6GkphVmd119pptFQgEgs3ATYervvzyy3zjG99YcdntYDBY8ToQCDA4OLjo+BMnTtyUfSVSjp2MfeD3yFnr0VJTbPvJ72OOLBQ3n1Wm5+N2/FaJ0Vmd3/3HJKlc7YlrpbIZ5vA5rHoKj1XjwlSCMVGqWyAQ3GasWBhmZ2fZu3cv3d3dfOc73+HJJ59k3759NfsYSstC8yn1dqiGz+dj//795SWl9vb2lZpM3t7ItQ/8HhnXHpRcnK2vHsM2+dqCcY0OmZ6P23CaJN4J6Xz5lSTZQu3iUCybEcUcOY9d0bFqCufHY0zHRRirQCC4fVjxUtLjjz9Ob28vBw8eLG8bGBigq6uLb3/728t+PhKJVN2+WLhrf38/Bw8epLm5mcOHDy+6jJTJZMhkrpe5iEajAOgGFHSDgubiSutTNL7+l9imf0bD0J8xde8vE236ZMV+djhl/vAXrDz1j0lemyjwRz9K8bsftKDU6CsoWBrQ4uMokgmndx+hbJ53x8IA5HIbX667ZIOwpRJhS3WELdXZjLasJpKxwkI/e/furepn8Pv9NTmfjx8/zsmTJyuc183NzfT09NDR0bFg/MDAAB6Ph2AwSFdX16Li8JWvfIWnn356wfZvfvOb2Gy28mvJyHPf1RPsCP8UgOGGT3Gh4Z/BDf2r352V6H1bJm9IPFqv8z836wg/skAg2Gwkk0k+//nPMzs7i8vlWpV9rnjGEAgEeP3113nooYfK286ePcuePXtq+rzH41kwOwiFQlWjkoLBIKdOnaKnpweAtrY29u/fT3d394LIpKNHj/KlL32p/DoajdLU1MR2p07DNm/F2HzTrxEd3oZ3+B+4a+YkWy0pJu/730C+fjruboI6b44/+FGaV6dktnvNHH7IjCTVpg5yNo6cS5Guu4+sycvlN37Kzgc/wL3bvTXPPtaCXC7HyZMnaW9vR9O0DbND2CJsEbasDrVGg66EFQvDs88+S2trK+3t7QQCAYaHh3nppZd46aWXavp8W1tb1Tv+1tbWBduGhoY4cOBA+XUgEODo0aNVl6PMZjNms3nBdlmiyoVYYvauz6Bbffjf/FtcY6+gZmeZfPg3MFRLedRHmkz8zqPwJz9N8613c7jNMr90/8JjVMXiRCaHbXYYw3c/ABOJHP54jp3+jY9U0jRtw/+gSwhbqiNsqY6wZaENq82Knc+BQIBQKERbWxuGYdDe3k4wGOThhx+u+fPzCQaDtLa2lmcMQ0ND5cillpYWTp06VTF+ZmaGlpaWlZpdlVjTJ5hs+S2MebkOSma2Ysxje0x88ZGiGPztzzP8w/naHckFiw+pkMIcPg+Aw6QSnE6ImkoCgWBTU9OM4ezZs+ULfyn6qL29vRwhFA6HCYfD7N69u6aD9vf3093dzYEDBzh16lRFDsOxY8c4cOAAR44cIRAI0N7ezvHjx8vC0dXVVds3q5FkwyNce/QoW878h0VzHT59t5l41uC/vJnl62fSODSJg7trU+m8rQElNg5YcJokptMGw1NxHjS7Mau3T+FBgUDw3mFZYXjttddobW2lUCjW/2lpaUGSJAzDKK+3l57XutYVCATKfoMbHc43Jrq1tbXR1tZW035vloxnL2Mf+D0aT/9JOddhYv+XyHj2lsf8mwfMxLLwX89nOf7TFDYNPri9BnGQZHK2OgjH0eJX8LubGY+muDidZN8WR80+C4FAIFgvll1KeuSRR8qiAEVH8czMTPnn/Oe3M8vlOkiSxBdbzLTtLpbr/sMfpXh9ssZKqnJRQEyzFzGnp/HbzVwMJZmIrrzFqEAgEKw1K/YxvPzyywu2zc7O8uKLL66KQatNNpdHr7GeXcHs5tqjR0nWvQ9Jz9Ew9Gc4L1/P7JYlid951MIHt6tkC/B730/ybqj2SqqGJGGKnMdKBrMic2EqLsp0CwSCTceKhaG7u3vBNrfbzbFjx1bFoNXGbTUxk8iQzdemDoZqYaLlt4ht/wUkDOre/E94zr9Q7uugyhL//kNWHmpQSObh3/1jkkvR2sShYPHNZUZfwGtVSKTzDE/GyYtKrAKBYBNRc7jqN77xDaAYRVR6XuLChQvlonibje0+K2mLlWuzKayqgt1Sw1eWVaYf+AJ5ix/v8H/FO/wPqOkQ0/cXcx1MisTTv2DjyHcTvBvSeeq7Sf60zc4W+zI6K0nkbPWYEmMUTG4anE2MR1O4rSq76xyr84UFAoHgFqlZGErZzoZhcPr06Yr3PB7PphUGTZZo8NqwaDJXQikiiRwemwbL+Xwlici+T1OwePG/+bc4r/4AJTPL5MO/jqFasGsSf/wxG7/9UpLLUZ3u7yX504M2vJZlxEHWyJtcmKPD6CYHboubkekETouG31FjjoRAIBCsITULwzPPPAMUax09++yza2bQWiDL0OiyYFEVLoeSTMXT+GxmFGX5iKBY0ycomD00nP06tumfVfR1cJuLRfd+ezDB1Zg+V67bjt209H51kwM5mcQcOY+j/mESWRieimM3q1g0EcIqEAg2lhX7GG43UZiPx6axt8FBncOyIr9DKdehoDnKuQ5qYhyAepvMM5+w4TFLDEd0/v33k6Tzy5efylvrUNNhTLPD1Nk1QoksI9NxdH1FpasEAoFg1VlxSQwoRiYNDQ1VhKj6/X6eeOKJVTNstYjn0/jnvbaaFPbU2zGrMmOzKSyqgqMGv8NSuQ47nArHPl7sAvfGdIHf/2GKpz9iRVtqRiLJRX9D7Aq6yUWdYxtXQincVhPbPNZb/+ICgUBwk6x4xvDUU0/R0dHBhQsX6O/vJxwOc/LkyQWlKzYLl1LT5PTKkFBNltjpsxGos1MwjGKJihpu1JfKddjrVfjDj9owK3DqWp7jP01RWObu31BMFEx2zLPD2PJRrCaV4ak40fTGl/IVCATvXVYsDC+88AKhUIhnn32WtrY2jh8/vsAZvZmYzM0ylppesF2SYIvLwr4GJ1aTwnQiTaGGpjxL5To8UK/y5Y/YUGX43qU8Xx9Ks1xVc93kgkIec/g8Hq1AOldgeDJOToSwCgSCDWLFwjD/Qtfc3FyuqrpUe86NRJYkLiUniGRjVd93W1X2Njjw2y2EkrX5HZbKdTiwVaX7A1Yk4L9dyPGffr58dnPeVoeansE0O0K93cRENM2lmeRKv6pAIBCsCisWho6OjnKW8+HDh3nyySc5cOBAzf0Y1hurbKKAzmhyfMGSUnmMVvQ77PDaiKVzxNM1ZCPP5TqEm/8VAN7hf6DujW+AnufjOzX+bWuxfPffvZVl4J1lqqnO1VMyxS5hSY3jtZkYmUkwFRMlMwQCwfqzYudzKWwVihnPJ0+eJBgMVrT63Gz4TS4m0iGupCbZY99WdYwmSzR5bVi0YkhrOJ7Fazctne+wRK7DP99rIZYz+OvXM3zjZxk+F5C4u2nxXRmKGV2zYY4M46h3kJIsDE/GcZhVrCYRwioQCNaPFc8YXnzxRc6ePVt+vWfPHiKRyKatlQQgSzJuzcHl5BThRZaUoOh3aHCai34Hs8JUvDa/Q2Vfh5+x9dVjyJlZPnevmcfvMQFwIihz+trSM5GC2Q2FNObIBfwWCKeyDE+JEFaBQLC+3FStpBub7bS1tVWtobSZsKkWdHRGE9cWXVIq4ZrzO9Q7LYSTGTK55f0OlbkOI2z76R+gJib4wkNmHtutYiDxzE9SXIktXVcpb6tHTU5hjl5ki8PMlXCSq5HUir6rQCAQ3AorFoaZmZkFDafdbvdtUXbbb3Ixk41yOTm57FirprCnruh3iGfyNfkdSrkOOWs9WnKSbT95GstskF/fb2G3wyCeg6/8IEUyt8QMQFLIW31osVGsmSmcZo3gVJxIUnR9EwgE68OKhaG1tZVvfetbFdtefvnlBbOIzYgsyXhMDi6nJghlo8uOV2WJHV4bzfV2dAxC8SzLRJ8uyHVofPUYnpmz/O93F/BZJC5GdY7/JIW+xI4M1YKhmIshrEqaXKHY9a3WTG2BQCC4FW6qJMav/MqvsG/fPr74xS/y2GOP0dnZycDAwFrYt+pYlWKhutHENbL68olkkgT1TjN3NTixW4p+h/wyfof5uQ6ynmXr2f/I/fEf8XsftqLJ8MOref7PN5eeARQsXuR8CnPkAnU2malYltHpxLJ5EQKBQHCrrFgYAoEAoVCIZ555hkAgQFdXFzMzMzX3e1539Dw33ub7TC5C2VhNS0olnBaVvfVOGub8Dunc0r6C+bkOGAb3X/17Hk18j9+cC2P9z29k+NGVpYUpZ6tHTY5jjV/Cb9O4FEoyKUJYBQLBGlOTMMzv2nb27FnOnj1Lc3Mz7e3tNDc3l7dtSmQzJKYgf/2CKksyXpODK6lJpjOzNe/KosnsqbPT5LWRzBSW9zvM5TpE9vwzAOrf+c98xnKGf7mv2OrzmZ+kuDi7hMDICnmzF9PsCI78DJoic2EyTkJ0fRMIBGtITXkMhw8fZnBwkN27d/OFL3yh6hhJkjZnvSR/ACwSxMchmwSLG2QZi2Imnk9zMTmOU7ViVkw17a7kd7CYFC6HUoTm8h2kxfIdJImZfY8TDs+yJ/JD6n/ex5ce/k1GInfxs6kCX34lxdfb7TgWKdVtaDaMfApzeBh//fu4kpQZnopz/zY3irx82XCBQCBYKTUJQ6lJD7Cp6yJVRTOBvwlsXpi9Upw9mJ1gsuEzOZlIh7mcmqTZvh1p0at7JZIE9Q4zFlXhUijJdDyD12ZCXayaqiTx9tbPUG/VcY7/mG2vf51nHvwd/tdTTVyN6fzxj1P8wS9YF73QFyw+tPg1zLNBGjz3MBZJ47Zq7PLbb/asCAQCwaIsu5Q0OzvLvn37yq8PHDiwpgatCRJFYai/F3zNoOcgMYWs63hNTq6kJpnJ1r6kVMJpKeY7NDjNRFLL+B0kmYkHvkCy/mEkPcfeN/6Mrz08Wa7G+rdL1VSSJHL2BrT4GPbkVdwWjZGpBKGECGEVCASrz7IzBrfbjWEYfPWrX8Xj8TA8PLyg53OJxZaZNg2qCp7txeWk2FWIT2FRzSQlmdHEOA7VhqXGJaUSZlVmd50ds6owFkmRKxg4F+vvIKtMPvzrNJ7+KpbwOzx6/mt85X1PcfQ1D3//dpa9XoWP7dQW/WzB4sY0O4K73sk13cHwVBybyS26vgkEglWlJufzmTNnmJ6eLi8jnT59esFjs/Z8rorFAf59UH83SAqeTIpIeoZLyYmbCgct+h2sNDc4kGDJfAdDMTHe8lvlPIdfHPsaX9hbLNPx1Z+mGA4vPuvQNTtIYA6fp96qMxPPiBBWgUCw6tTkY3C73eXieR6Pp6KQ3m2LLIOjHsxO5NkxfNHLjIUv4FVt1Fv9y3/+BiQJ6hwmzKpc9DvE0njt5qp+B0OzMd76BFt/+keYEmN8Sf5TLtT/Dt+bsvOVV5L85WN2XObqmp23+NES17DOBvE77+JyOInbprHVLbq+CQSC1WHFeQx1dXWMjo4CxTDWffv2sW/fPr773e/WvI9gMMjx48cZGBjg+PHjRCKRRccODAwQiUSWHHNLaBbwBzA3PoRidjEy9Qap5M2X9yj7HVwWwsnson4H3eRk/MAR8hY/ptQkf276Os22FOMJgz/44RLd3ySJvLUeU/wKzsw1LKrC8GScmOj6JhAIVokVC0Nvb285ma2zs5Pu7m6+/e1vc/jw4Zr30dnZyZEjR+jo6KCjo4NDhw4tOdbr9eL1epEkCUmSOH78+ErNXhoJsHnwbm0l6qjnUmIMY/YKFG7OuWtWZfbUOdjlm8t3SFXPOyhYfFw78BQFkwtb4jJ/530Wj5rh7GSB515f3BltKBp5kwvzbBC/HCeRFV3fBALB6nHTHdxeeuklwuEwX/jCFwgEAjUX0QsGgxWvA4HAot3fIpEI/f39GIZRfvT09HDkyJGa7c0Uas8UljQNX8O9jDl8TJntkJiG5MyCzOlaUGTY7rWyd4ujnONQbTd5+xauHehGV234EsO8UP83qOR54VyWkyOLC5NucoBRwBy+wBaLwXg0zeWQ6PomEAhunZsqiXH06FG6urrKs4SRkZGai+gNDg7i8/kqtvl8PoaGhqqO7+joKD8fGBioeF0Lr029Rn6ZMtvzMStmNIuHUYuNZP09IJtg9irkbu6i67ebCDQU8w1CiQxUEYecs4nx1ifQZRO7Um/yzfr/gozOn55Kcy60uDM6by22BLXFR/BaVUZmEkzHRckMgUBwa6xYGPr7+wkEAnR3d/Pss88CxVyHo0eP1vT5xXwFoVBowTaPx1PxuVAotOIqrpF0hJ9e+ym6Ufsyi8fsIZZLcEnKoW9vgbq9kI5BbBz0pWskVcNhLvr4rZpCJFXdF5Dx7GWi5bcwJIWW7Bn+3PP35HSDp3+QJJxexHZJJmerxxS7gjs7iYzEhck4qezKbRQIBIISK27t+dxzz9HR0cHu3bt5+eWXOXz4MJIk0dvbe0uGLOdc7u7upqenZ9H3M5kMmcz1u+VotFhWW0bmavQqZ6QztNS31Jzd7NE8XJ29ilNx0uBpBrMXQkEIXy3mQVicNe0HKK/9b3VbuBxKEUvmsFXJdUj47mf8fV+k8fW/5DHjR/yuzcofJf8VT7+S4pmPWdGqZVZLGrpiRQ1foM5v4UrczIUJuKvBiVwlkzqXy1X83EiELdURtlRH2FKdtbBBMlYYBL9v3z7Onz8PgN/vp6enh09+8pN86lOfKm9fir6+Pnp7eyvyHrxeL/39/bS1tVX9TCQSYf/+/QwPDy+636985Ss8/fTTC7b/+X/5cy5oFzAMg4AWYI+6Z1kbN5rt4Z/wwJVvUgD+LPMveC73T/nwFp3HA8K5LBAIKkkmk3z+859ndnZ2QRO1m2XFM4ZqzmegZudzW1tb1dlFa2vrop85ffp0xbJSNY4ePcqXvvSl8utoNEpTUxN7H9pLk7mJ1yZfI0yY5i3N7HbtrslWwzCYTE7SaG9kn3cfsjS38paJQngUotdANYPVy+JV9IozhpOvj9H+0DZUWWYskuRKJIXHWj3PgaZ/SsQF/nN/z29b/huhgp0XJj5K6w4r/7R5kcxso4CWmCDr2sWUeTc6Bg/u8OCyVGZS53I5Tp48SXt7O5q2SJb1OiFsEbYIW26dteieuWJhKDmf+/v7b8r5fOO4YDBIa2tr+cI/NDSEx+OpGDc0NLTAYX0jZrMZs9m8YLukSOzz7yNtpHln5h1em34Nm9lGo72xJnv9Dj+T6UnqcnXXP6P5i7WX3FthZhjiY2D3g2Zbcl+aIqMpMjt8dvI6TEbT+B0W5Cqenljgn6EWkniG/3982fo88ZSVv3ztAHu8CvfXVfu1qRg2P9b4ZRosHq7oPkZDaR7cbsGkLjyApmkb/gddQthSHWFLdYQtC21YbVbF+RyJRGp2Ppf20d3dzcDAAL29vfT395ffO3bsWNVucLfaOvQB/wPscu1CN3R+PPZjQumFzu5qmBQTZtXM6OwoyfmRSbIMrm2woxXq9s1zTi8fAaXKEk0+Gx6biVCyeqQSQHjvZ4jubMOmwh9b/jMf4Of8/ispppPVl5QM1YKuWTFFLtBgyjAVy3IpJEpmCASClbHiGYPb7V6QkPbII4/wyCOP1LyPQCBQdiTfGH46XyRKrCRvYTEkSaK1sZV0Ic1EYoIfXv0hn2j6BA6TY9nPeswexhPjXIxe5G7f3deXlAA0KzTcC/Z6CA0Xl5csrqKDegnMqszOOhsXJuLMpnK4bVVUX5KYufdfI+eSuMZ+xJ/o3+DXMr/G06/cw9cO2jFVWYYqmD1oiXFs0SB+5z2MTidxWTQaXJZlv6dAIBDATQgDFEthDA0NVaxt+f1+nnjiiVUzbC2QJZkPbv0g37vyPSLpCK9cfYVP7PwEZmXhEtSN+C1+xuJjeMwetjq2LhxgryuKwewVCI0Ucx/sdUUfxCI4TCq7/HaGJ+Mk0nns1aqySjJTD34BOZ/EN3GW/2A8yxfD/5Y/P72X33nUUjXKKmetR0tM4NLsJNSdDE/FcVhUbKab+nULBIL3GCteSnrqqafo6OjgwoUL9Pf3Ew6HOXny5Obs3lYFTdH4yPaPYNfsxLIxfnj1hzUlwGmKhk2zMRodJZFLVB+kaODbU1xecjfVlDnttWk0+ayk8gUyuUWijmSVyYd/g2zdPWwxZfiP2tc5N3qJfzi/SJiarJC3ejFFL9IgRZhN5RieTCxef0kgEAjmsWJheOGFFwiFQjz77LO0tbVx/Pjx266rm1W18pHtH8GkmJhJzdScAOcyu0jmkozOjlJYKtHN4oLGB2F7Cyjm4uxhCRqcFnZ4rETTWfKF6hdvQ9GYaPltdG+AraYkX9f+gv/62lVen6guaoZqxVBMWGYv0GjJczWS4mpYlMwQCATLc9O1kgCam5t56aWXABatd7RZcZldfHj7h1FkhbH4GK9NvlaTk7bOVse1xDUmkhNLD5zvnPbvLW6LT1WdPUgSbPVYaXBaCCcy6ItolKFaGW99AsW9nR3qLH+u/Tl/+aNxJhLVP1CweJFzCezRIG6TQXA6QSQpur4JBIKlWbEwdHR08OKLLwJw+PBhnnzySQ4cOMCePZs/cexG6qx1vL/x/UiSRDAS5J3QO8t+RpVV7CY7I7MjxLKx5Q+iWaH+ruJzkw2iV6tWbZ0fqRROLl7vqFiuuxuTq55d8jR/qP8Ff/KDadL56qKWs9WjJa7hy14jXzAYmV5kGUwgEAjmWLEwPPPMM3z6058GihFKJ0+e5JlnnrntlpNKbHdu5+H6hwF4Y/oNRmdHl/2My+Qik89wMXpx6SWlG2l8GNw7izOHzEJRKUUqWTSF2eTiae4Fi5eJA91YHW72yVf5jcRf8PWfRqrPeGSVgsWLaXaERnmWmbk+0SKEVSAQLEZNwjA6OrroQ5Ikmpuby817bkf2evdyj+8eAM5MnGE8Mb7sZ/w2P+OJca4lrtV+IJMVGh+ALfdDPgPxCbjBt1GKVDIwSKQXd4rn7VuYen83FpudB+UR/vn4/4cX36k+G9A1G4YsY40OU6cVhexyKCnEQSAQVKWm+MWWlmLxuaUuJJIkrUlq9nrxQN0DpPIpLkYv8uOxH/Pxpo/jtXgXHa/KKg6Tg4vRi7jMLlymGmuUyEoxcsnsgqlzRce0o6EirLUYqWRjZDqBmpMxa9X1O+dsIvL+J3H98Bjvz75N4q1ezrh/nf3bFobIFix+tMQ4zuRFAC5MxZGVoghVK7YnEAjeu9Q0YwiFQszMzBAKhRZ93M6iANcT4LbYt5DX87xy9RXi2fiSn3GanGQKGS7OXlxRzwegWEJj+yPg3VMMa01HK97e4rSwfZlIJSiW6449+luoqson5dfI/PSvGYtWWd6SJHK2OkyJMQBcFo3zk3FGphPoIoxVIBDMY8U+hjuZUgKcx+IhnU/zytVXlu0A57f6mUxOci2+giWlEpq1uKzU+EDRIR27Vl5akiTY5rFSPxeptNSqT7r+QeL7fw1Zlvkn/Ijhf/z/kspWiVSSiy1BAZy5adxWjQuTMYLTcSEOAoGgTE3C8MUvfrEciVSN5557bkW1kjYzpQQ4m2arKQGutKR0KXaJ2czsyg8oy+DdDdv3g8VbXFrKZ+b2LbGzVFMpsbRAZbcdIPzgryBL8M/zJ3n9H1+suvSnz5UAsYTO4c5O4LWaGJ5KMDwVFwlwAoEAqFEYnn/++UV7JQAcOnSoauG72xWrauUXtv9COQHu1WuvLpkA5zA5yBayjEZHV76kVMLmg22PgC9QzJZOR4DaI5UApMDHuLr3fwbgsfh/5a0f/9+LjtVVE+aZN3FlruKzakVxmBTiIBAIahQGwzCWbQBxu/sYbsRldvGhbR9CkRWuxq9ydvLsks53v9XPVHKKsfjYzR9Us0DDfcWsab0A0WIr0fmRSsnM0sJjf+Cf8+7W/xcAH5z4P7n68+9XHaebXOiaA0v4HVypy/htKiMzcc5PxMgXREMggeC9TE3CEAgEuHjx4qLvr6Qfw+1Eva2eRxsfRZIkhiPDnAufW3SsIis4Tc6bX1IqIcvg2QnbWooNgKJjkM+UI5WSuQLZ/NIX7sb3d/Ka65MA3H3hG8RHz1Qdp5scFExuzOF3cSUv4rdpjM4keHcyVm5HKhAI3nvUJAyPP/54uSlPNX71V3+Vz33uc6tm1Gpyq7H6O5w7eKj+IQB+PvVzLkYXF0iHyUFezzM6O0pOv8U+rDYfbHu4uLSUmIZUpBypNJtaOlIJSWLnx/4XXjV9ABmdxrN/gTHxVtWhumYjb/Fgjgzjio1Qb9O4NJPk/IQQB4HgvUpNwnDkyBF0XWffvn389V//NWfPnuXs2bN84xvfYN++fQCbtuR2oVAgMrdef7Ps8+7jbt/dAJweP71kApzP4isuKcVuYUmpRGlpaetDYOhIsWtsc5lqilTSVIXGTxzmlPQ+ZCOP86f/AS0SrDrWUK3kLV5Ms0GcsQs02BUuhVKcG48tOzsRCAR3HjWHq548eZInn3ySY8eO0dLSQktLCz09PRw5coRvf/vba2njLbHXuxdZkplITKysfMUNPFj3IDtdO8sd4MLpcNVxiqzgtri5FLt0y4IEzC0tNRUrtdr9qPFr7HTJNUUqeWwmlI/8JkPGXciFNNYfHkeLV6/0aqgW8rY6TLFRnLPn2WKTuBJO8s54lEz+5s+bQCC4/VhRHsPhw4e5cOECuq6j6zrnz59f0M1ts1Fvq+eBugfwW4r5Bul8+qb2I0kSrVtaabA1lBPgFuvLYNNs5I08o9FRcoVbXFIqYfXC1ofBvxdzNsIuewZzDZFKe+ssXHno3/KWvguycdw/7sGSrR4oYCgm8tYGtNhlHLPvssUKV8Mp3rkWI50T4iAQvFd4TyS4OU1O7vPfx27XbqKZ6E07hxVZ4YPbPojb7CadT/ODKz9YNAGuzlrHdHKaK/Ert2J6Jaq52EZ060PYNY09plkMI79spNLHm928tOu3GTG2oifDPDT8l2iLlA03FI2cfQum+FWcs+fYaoNrs8VlJSEOAsF7g/eEMEAxca3Z08y9/nsxDIPJxGRNzXluxKSYakqAkyUZt8XN5dhlwpnqy043hSSBewdsfwRP3Tb2aFHS6cSyvoB/3eKn1/1bXDXqUFLTbPnR72OOXKg+WFbJ2htRExM4Im+zzWpwLZLi7WtRUlkhDgLBnc57RhiguBzUaG/kgboH8Fq8TCQmyFbpjbAcNs1W0QHu1fHqCXA2zYZu6Fyevbwa5ldi9cDWh/Hvuo8dlgzJyOSSkUqKLPHrH9nK76q/w9t6E8lEFN+P/gj7tZ9W/4CskHNsQU1NYw+9yTZbgfFomrevRUlmbzKJTyAQ3Ba8p4ShhNvs5l7/vTQ5mwinw7U13Kmyjw9t+xCyJHM1dpXXJ1+vGhrrt/oJZUKrYfZCVBNSw71suev9+BwWUtOXMAqL39G7zDL//uA2njH/Nj/QHySRyeE4/XWcw/+tel9qSSFn34KSCWMPvcl2S47JWJq3xqIkllm+EggEty+3JAxnz57lxRdfJBqNLj94k2FWzOzz7uMe3z3k9TxTyakVLy3V2+p5dGsxAe5C5ELVBDhZknGb3QBcjl6+qRnKkkgSqncH2+//CFbfVpIzl5DyqUWH+6wyX7hf5ZVdv8aJwsdJ5Q1Mb5zA+bO/hmrlPCSZnL0RORstioM5w3Q8w1vXosSFOAgEdyQ1C8Njjz2Goij8u3/374BiUtsnP/lJ/viP/5g9e/bclo16ZElmm2MbD/gfwKk5b2ppqcnZtGwCnFW1AjA8O8zrU6/fcuhsNSxOL9vu/QAF/10kZ2dQFgmnBVBk6GqxYz7wv/DnhU4yuoQy8j3sP/4qci658AOSRN62pdg/OvQmO8wpQoksb43NEkuvUtSVQCDYNNQkDEePHiUQCBAKhbhw4QLf+MY3CAaDhEIhTp8+zbFjx+ju7l5rW9cMj8XD/XX3s8O5g1AqtGwfhhvZ593HXb5iX+fT46eZSFSP+GmwNZApZHhz5k3eDr29OnkO83A77DTte4iY736SOQMtfq1Yc2kRPrlL46Nt/xPPKF0kdRPK5BuY//H3UVLTCwdLEnn7FqRCBtvMm+zQYkSSOd4cixIV4iAQ3FHUJAwDAwMcP34ct9tNX18fXV1dHD9+vPz+4cOHGRoaWjMj1wOLamGfZx93ee8iU8gwnZpeUTmN99W9jyZnUzEB7tqPq170ZUnGa/Hit/qZTk3zs+mfcT58nmS1u/SbpMFlZeeuZiac95Mw1aElx5GW2H/Ao/C//pMP8ReeI0wbbrTYFUwvfxk5XD1LOm+rB6OALfQWO5TijOGtsSizKSEOAsGdwoqrq3o8HtxuNw8//HDFGLfbverGrTeKrNDkauKBugewqTYmEhM1l9GWJIkDjQeot9WTK+SWTIBTZZV6Wz12zc6l6CVen3qdy9HLq5YMt91jo6lxC1cse0k4m1GyMZR0qLqDGXCYJA594m7+e+Dfc0Hfjpadxfb9PyJ36XTV8QWrHwMJa/gtdkgh4uk8b49Fl022EwgEtwc1CYPXW9n72O/3LxgjSbX3DQ4Ggxw/frw8E4lEIkuOHxwcpK+vj8HBQQYHB2s+zs3is/i4338/jfZGppJTNd/RK7LCh7Z9CLfZTSqf4pWrryzps7CoFrbYt4AE58Ln+Pn0z2/KCX4jsiyxp97OVp+Ly/IOknUPgqShJa5VdzADsiTxrx7eyuVHf5fTxn3Ieoa6M/+R6M//e1VBKVi8GLKGJfQ225kkkcnxxtgs4cQqO9cFAsG6o9YyaHh4mK9+9avl1zMzMxWvoXixr5XOzk7OnDlT/tyhQ4fo7++vOnZwcJD+/n56e3sJBoO0t7czPDxc87FuFptm427f3Tg0B6PRUdL5NF6Ld1kBLCXAfffyd4lmovzw6g/5SONHFh0vSRIukwuH5iCcDvPG9Bs02BrY4dxRjma6GTRF5q4tTjJ5nWtJD1vrH8IUuVD0O2Bd9HP7m1yMuZ7k5e//DZ/M/4AdF77JcGScug/9MpJS+edSMLuRpRiW8Dm2u3WuZrbw5rUo92114bObbtp2gUCwsdQkDAcPHuTVV18tv25ra6t4XRpTCzcKSCAQWHIW0NXVVRaRQCDAyZMnazrOaqDKKrvcu7Cb7AQjQSYSE9TZ6lDlpU9bKQHue5e/x3RqmlMTp6gz6pb8jCzJ+K1+coUcE8kJQukQ2+zb2ObcVo5qWikWTeGuLU5+fnWWqayEv+4BFNUO4asoqVCxtHcVodvm1vD9k0N8+/tb+FR0gObpl7lwchrXx38Di6XSFt3kBCQskXPs8OpcyTby5tgs9251Uecw35TdAoFgY6lJGBa7m78ZBgcH8fl8Fdt8Ph9DQ0O0tLRUbC9FPnk8HoaGhggEAhvSEKjOWodNtRGcDTKeGMdj8Sx7sXab3Xxw2wf5wZUfcCV+hXA+TMAIIC+zeqcpGlvsW0jmkoxER5hOTdPkaqLB1rCsIFW1w6px1xYHb1ydJZY1sLv3AFdBNqElrpGzNUCV/Vo0mQ988l/yg6EG3n/pOfamfsalb/8+0Q8/QUNd5VKibnKAJGMOn2eHu8AY23lrLMq9W13UO4U4CAS3G+ue+byYPyEUWpgdPDQ0hM/nY2BggEAgQF9f34b1li4tLQU8AeLZ+KJlt+fTYGvg0cZHAbiSv8LJSye5FL1UU7STTbPRaG+kQIG3Zt7ijek3VhwpVbbDaWFfg4N4Jl8uhJeuf5CcrREtMVE9d4HiMtdD+z/E+YefIio52alfxvODL/P2hYXLhqWGP5bZYbbnL2Hkc7w1Nstk9Oaq2QoEgo1j5beg8/jsZz/LiRMnVsWQaoIRCoUIBoO0tbXh8Xg4fPgwXq+36sUxk8mQyVyvdFrKxs7lcuRyqxcts8O6A4tkYXR2lPHZcfxWP4qsLDp+u207++v28/2r3yeajvKTqz/hbfPb3Oe7j232bcv6LJyKE5tkI5KKEEqEaLA1sN2xHYfJsSK7G+wacbeZkeniecnLJpLee9BUO6bZi5BNUbB4qy8t7bqLiPPLpH/0NRoK17D97I94eeZX+WDrfuT542ULBZMXNTJCgzPHpLaLN6+EyDW6aKgycyj9Xlbz93OzCFuqI2ypzma0ZTWRjFvofXngwAFOnTq1os/09fXR29tb9htAMeqpv7+ftra2irGDg4N0dnYSDl+/O5ckiTNnzixYdvrKV77C008/veB43/zmN7HZbCuycS3IG3kuFy5zKX+JvFGMDHLKTgJqAL/sX1FU14aRS7L13N+wK/suBWS+aepk290fxq5ttGECwXuXZDLJ5z//eWZnZ8tpBbfKugtDMBisiEqCojCMjIzg8XgWjN2/f39NwlBtxtDU1MS1a9eqhteuBrlCjiuxK1yJX8GkmnCbqkcR6Xmd0dOj7G7dTV7Kcz5ynvOR8+UcCb/Fz33++2iwNtQkEIlcglg2htPkZIdzB3WWuiVnLRWfTWX4x+++hLbrYRwWMy6LiiRJyNkE2mwQLTlO3uzB0BYRUz1P+id/zdapVwD4B7mNHR/5PAHvDeqg59CSU+RsjUxa9pCTVO7a4qLRbSkPyeVynDx5kvb2djRtY9VF2CJsuV1tmZmZYevWrasqDLe0lLRnz54Vf+ZG53EwGKS1tbUsCkNDQ3g8nrKjubW1lUgkgsfjIRgMEggEFogCgNlsxmxeuFyhadqa/eI0TWOfeR8um4tgJMhUZop6Wz2yVN11I6syFtXCg1se5K66uzgXOsdwZJhQNsQr116h3lbP/f77qbfVL3lcp+rEYXEwm5nlncg71NvqaXI24TUvH05rn/t53zYPV2azjMfz+GwmrFY3OfMDSFEXptkRDD1NwVIlaknWsH+4i5k3GqkfHuBf6oN8/x+n+f7DXXwi4Jw3zoTu2II1Mc42GSatzZybSiIpCju8laKzlr+jlSJsqY6wpTqbwZa1OP4tCcPzzz9/U5/r7++nu7u7POOYH/V07NgxDhw4wJEjRyrG7t+/nzNnzqxruGotlHo82DRbOaTVZ/VhVpaOxjErZt5X/z7u8t7FO6F3CM4GmUpO8b3k99hi38L9/vvxWxef6UiShMfiwaE7CKVDRNIRGu2N7HDuwK7ZF/1cie1eGw0eO5dDKa5GkkTTuWLugTuAbnJgDl9YPGpJktAe/FfMuhtwn+3jo4WzvPlaD387/Rv80v4GNGVOTOYa/mjJCRoMnSlrM+9cK+bL7fDeXAiuQCBYe2oShtHRUSKRSEUZjJdffpmenh5CoRCf+9zn+J3f+Z2aDxoIBOjp6QGgo6Oj4r0bQ2M9Hg+9vb0173ujcJlc3Oe/j9HoKFdiV7BqVlym5ad1FtXCww0Pc5f3Lt4Ovc3o7CgTiQkmEhNsdWzlfv/9eC3eRT+vyioNtgbS+TRXYleK4a3OJhrtjZiUpZPMbCaVuxudNDjNXAolGZ9NYdFUPLYGdNVeTIhbYmnJ2PkhEnY/tp/8KfdnR/Fd+SP+NPIb/MovNOO3zs2a5GJPBy0xQYNeYNq+l3PjxTIrjc7NcdcnEAgqqSlctaurqyIx7bXXXqO9vZ329naeeeYZ/u7v/o6vfe1ra2bk7YJJMbHXs5d7ffei6/qKylvYNBv7t+znU7s/xR73HmRJ5lr8GoMXB/nR2I+W7VNtUS00OhpRZZV3Q+/ys6mf1Vze22s38cB2N+9r8mDWZK7NpokaFtJ195Px7F2y1lLOfzeJj38FybGFbVKIp+LH6fv2Gd6Ymld6Q1LI2RtRMmHq4udwkeTcRIzLodUrHigQCFaPmoTh9OnTfPrTny6/7u3t5TOf+QxPPPEEBw8eZHBwkGeffXbNjLydkCWZrY6tPFD3AC6Ta8U9HhwmB62NrTy2+zF2unYiSRJXY1c5efEkPxn7CdHs0k2RHCYHWxxbSBVSKyrvrcgSW91WHm7ycPcWB9lCgbFonph9N+n6983VWhqvWmspb28k9rGvINXfjVtO8wfG1/nu917iv76bvR5aXGr4k4tRF3sHL3EuTBXLm+cLt1YbSiAQrC41V1edT39/P5/73OfKr0uOYcF13GZ3ucdDLclwN+I0OXn/1vfTvqudHc4dGIbB5dhlvjP6HU6Nn1qyZ4QsyfgsPvxWP1PJKX42/TMuhC/UVAzQoinsqXfwyE4vTT4rkVSesYKbeN37yNm2oCUnqpbx1k1OZj/4FOz8ADZF53fV/0Lh9efp+XGKdL4kDsWGP1IhiS/6Dj6KLVVfvxLhaiRVTr4TCAQbS03C0NbWxl//9V8D8MILLxAOhytyDkZGRm4qQulOx6yYuct7F3s9ewGYSEysuPdCqbRG+652tjm2YRgGo7OjfHv025yZOLPk/lRZpcHegF2zczF6kdenXudq7GpNx3VZNO7d6uKhJg9uq8a1tMKk/S7S7nlLSzdgKBrhR/7f5O/5l9g0if9N/b/5+FgfT5yMcC0+NyuQJPK2BiQjhzd2HoB0VueNKxFOXwwzMhUXLUMFgg2mJmF47rnnOHbsGH6/n8cff5z+/v6KeNmenh66urrWzMjbGVmS2e7cDhTbgGYLWa7FrxHNRFdUXttj8fDh7R/m4M6DbLFvQTd0gpEg/2Pkf/Da5GuklujzPL+89/lI8WI8kZwgpy+dMSlJEvVOMw81ebh/mxtD0Rg1thF23weSWr1DnCQRuauD+MOHcJpVPqWc5onkn/HUt8d5dez68fLWOgyKM4nG5Dl2KiHUbJxz41HOjIZ451qUSDJ7UyVABALBrVFTVJLb7ebChQuMjIzg8/kWNOXp7Oysubrqe5mAJ0ATTcykZhhLjDGRmMCsmHGZXTUXyPNZfXx0x0eZTk3zxvQbTCWnuBC+wOjsKAFPgLu9d2NRLQs+VyrvbZNsjDLK2zNvcy11jW2ObdRZ65YMr9UUmR1eG367mcuhJFcjCrMWM1vlS1iS4+TN3gVRS/EdHyVv9eM/8+e0pIL8ReFP+K0f/Brn7m/il+43IUtSsQQHYdR0CDU9gUUxU2f2EDN8XJ2wcjXipN5pYavbis9uQpFvg+xwgeAOYEV5DIstFwlRqB2bZisXyJtJzzCeGCeUCiFJEk6Ts+pFvRp11jo+tuNjTCYneXPmTWZSM7wbepdgJMhe717u9t5dNVy1lHzXYGsgrsd5e/ptHCZHWSBsi2U8A1aTwl2NThpcZi7OmLkSMeE0LDRkrmAUMnMX+uuk/fcz+aHfo+H019gbneJv+BOefLOLd0P30P0BK9a5v768rQ5DlpDyaeRMGG9iHI9qJp1xMBPzMjPtxOn2s91nx283Y1LXvfajQPCe4pYS3AQ3j6ZoNNobqbfWE8lEmExOMpWaIpKJ4NAc2DX7slnMkiSxxb6FBlsD44lx3px5k3A6zDsz7xCMBNnn3cc+zz40ZWG+gCzJeMwe3CY38Vycd8Pvcjl2mUZbIw32BpwmZ5UjFvHYTLgsGpMuC6PTVkZDNurTozjypYS46+U5co7tjH/wy2wZ+lO2zgzzl9Kf8/vj/4Zf+877+T8+VCmChmqhoFooAFIhgzmXoCkbIpdViMasvDteh81dx5b6eurddqym2sqACASClSGEYYNRZAW/1Y/P4mN7bjtTiSkmkhOMJ8bLSXKLldgoIUkSWx1babQ3MpYY483pN5nNzPLm9JucD5/nbt/d7PXsrbpcVZqpOE1OErkEo9FRxhJjNNga2GLbgtvsripQsizR6LbgtWtcc1u4OuEkNvMudbNjSM46jHn9KgpmN9cePUr9z3rxXzvFH0r/ib9KTfNbL/0TPr1bIrDNWLBMZChmCoqZggWkQg5PLo43N0pifIRrEzYmXFvx+Ruoq2/AbRdZ1ALBaiKEYZNQ8gG4TC62Obcxk5rhWuIak4lJVEXFbXJXvfO/cR/bHdvZZt/GldgV3px5k1g2xs+nfl4WiN323Yt+3q7ZsWt20vk0V+NXGU+MU2etY4t9Cz6Lr6pAmVWF3XUO/A4zl90upq6cwxa+iMPuwLBdX1oyFDOTD/86XtvzuIP/nV+X/hs7clP88fDn+b8uJ/jELpXH9pjY55UXCJGhaBQUL1i8mPUctmyCzOwwkelzRC86cdTtwF+3BY+vHlkTjYEEgltFCMMmxKpa2eHcwRb7FkKpEOPJoh/CwMBtdi/rh5AkiSZXE9ud27kcu8yb02+SyCV4ffJ1zsnncOQd7CjsWHQ/FtWCRbWQLWSZTE0ymZzEa/GyzbENn8VXdebhtGjcu93PjHs/Y1d8hMbexp66gsmzFZS5JR9JJnz358hbG/C/9Z/plH/K/cplfpC/jzeGd/OV83uwu32079Y4uFu7XlZjPrJGweJBtXhw6QUK6RiJq++QGDvHjNOFt2EHHn8jqt0LqhAJgeBmEMKwidHkYpvPets8P0Ryikg6gsO0vB9ClmR2uXbR5GziYvQib828RSKTYCw3Rmg0xA7XDna5dtFga6g6GzApJhpsDeT1PLOZWWamZvCYPWxzbMNv9S9wbkuSRJ3Tgueuu5n0+ZkK/ozM1EXM7gZM1uuNhWI7P0neWkfD2a/zvsxVHjSNkSlAtmAwlfLwxpu7+e9v7KHga2ZvoJlHdzowKVW+p6yg2DzYbR7y+TyxeJTo+TeYufQ2Po8HV902zO4tYHGDJpabBIJaEcJwG1DKZPZZfGx3bGc6Nc14YpzxeNEP4TQ5l+zHIEsye9x72OncSTAc5NWxVynoBS5FL3EpegmbZmOXaxe7XburdoZTZRW/1U9BLxDNRnlz5k2cJifb7Nuos9Ut6H+tKjLbtjTgdX2EydE3mR17l9lkHLtnC+rcBT5V/z4ufvgZZt/9CXuUcazRYeyxK9jzs2wtvM4n9LMQBf2szMWzW4k7A3i278O/fS9553aQKr+vqqo4PT503UcinSUSiuIIv43f9i5utwerewvY68DiAdPGN24SCDYzQhhuM0qO4q32rYTSIcbiY0ynplFkBZfJtWRFVUVWaHY3gwl8TT4uJi5yOXaZZC7J2zNv8/bM29Tb6tnt2s0O544FS0aKrOC1eNENnXg2zrnwOS7HLrPVsZV6a/0CUbFaLey6p4Wor56pkdeJhi5TsDbgtFmQZShYvFz1fgBHkxdlLlzVHB3FHBmmMH0BY3oYUzbEHuMqxK7COz9AOSchaWYMbwD8e8l4msm4A+VQWVkGp80E1joSmTzBdBZrJk1d+Dxe6zB2hxPJ5gd7/ZxI2Ku2MxUI3ssIYbhNsagWtjm20WBrIJwOcy1xjXA6TMEo4DQ5l8xHkCQJr8WL3+HnofqHGIuPMRodZSI5wVRyiqnkFK9NvkaTs4nd7t34LZWtR2VJxmV2lSOZgpEgY/Ex6m31bLFtwWVyXR8vSbi27MThdBG59AbT4xcJzTqxWJ1Ybgg3NVQLad89pH33wFw/p1QqxNjF88xcOY8tGuRu4yLWbBom3kKbehuPAiZFomD1kXE3k/E0k3Y3k3Xvxm4xY7eopHNmrmTyjMXBly1Qn7iCM3IZRbOC1QeOBlCX72EhELxXEMJwm6PKKvW2euqsdcxmZplITjCdmmY2M4vdZMehOZYMd1VkhSZXE02uJpK5JBejFxmNjhLPxhmZHWFkdgSnycku1y52uXZVCI4kSThMDhwmB8lckivRK4wnxqm31rPFvgWP2VM+tmzz4Nv7KA63D8+Vc0wkI0xnl78YG1YfW+95P1vveT/JnMH/dSnNm8OXUCNBHpBHuT8/SrN8DUtuBnNiBvv4KZDAkGRyjh2kPc1k3M04Pc0kzVsIZySmswoui0aDRcdTmEaLXgVpzlE98QaY7aCaig2KZA0Ube65ev25mGUI7mCEMNwhlDq6eSweduR2MJ0s+iEmEhOYVTMu0/JlN2yajXv993KP7x6mU9NcjBaXmmLZGG9Mv8GbM2+yxbaF3e7dbLNvq/BrlDK6M4VM0f+RGMdv9bPVvhWfxVccq5owNd5Lg92DY+wtJqdDzGBmOp7GoqpYTcqSWc02TeKxZiuPNd/Ntfg+Bkdz/KeRLJFEinulS9wvj/KoaZRW7SI+ZjHFLmGKXYLL3wVAV61kXHtIuQJEbHsYdezEZPfT4HLhUHJAGmITEC9U9p6QpKJPQ1FBUosJfKoZVCuoFlC1xQVEiIjgNkQIwx2IXbNjd9tpdDQSSoW4lrjGTGqmmCuhLN9VTpIk6m311NvqebjhYa7ErjAaHWUqOVW+6JsUEzudO9nl3lXRa9qsmGmwN5Ar5Ahnwkynpq9HMln8xVwM1zZsJjvbLG/z+swEzZYE8RzEExJxXUZWzJgtZiyasug1datD5t88YOZf32/ijSkr3xl18sKlu/gvKSAF9YT5F95LtLsucS+j2GIjyPkU1tBbWENv4QMwIGPxE3PsZsbVjD/vY2Z2D5qrAcVkQ1EkVFlClUAxdGTyYBSKPSkyGUiFi89vFBF5TjwktRiqq1pAsRQjo0qCceNPWSsKj0CwCRB/iXcwZsVcdAzPhbuOJ8aZjk8DEM1GcciOZWcRqqyy272b3e7dxLNxRqOjXIxeJJlLciFygQuRC7jNbva499DkbCrnRmiKRp21joJeYDY7yxvTb+A2u8sCYbG4kbY8CExQt30vjXqKdDJBMp0mnogTj82QKBQACU0zYTabkVUThqxiKFo5KkmSJB5sUHmwQeXXWiz88Eqe74xkeW3Cy9+EvfxN+CEsCnx0h8ynG6d4SBnBMjuMeTaIKX4Vc2YGc3oGY/o0jQUD+WJRiQqqhZzJS8rsJW/xUbD40K1esPmR7XXIDj+yxYWqKqiyhDL3UCVQmBMPPQ+FPORnQZ+ZE5F5J7csIur1GYlqAeaWtcKjoGkgyQsfSHPPpervL9guZi2C2hHC8B5AlVXqrHX4LX5C1hCv8AqqpBJJR8jpOTRFw6basKiWJf0RDpODB+oe4D7/fUwmJxmdLZbPmM3McnbyLD+b+hlbHVvZ7dpNo70RWZJRZAWfxYdu6ESzUd6ZeQe7ZmerfStebS4zuuEeJE3DqutYCxn8hSzpdIpEIkksESc6O0skGcVIpbDIOhZZx6QABhiyhDF3x22TNdqaNA7usjGZNHjpYo7vjOS4EtP5zkWd71z002Cro233h2h/SKPJksEcHcEcCWKKXCAXHsNFFCWfQtYzaKlxjNQ4pau5Uf5fEV1WyZk85Mw3CogPyV4UEMXmRVWtxdmHIhcnELI8JyLGdQHRC1DIQS4F+bl+FJPvQLX8jRISFCvnS/OEYO4npefztpdnMsr1GY2szG1fRHhK3fXSs2DYijOcZTLwBbc/QhjeQ5TKbgA8VP8QWSlLPBsnnA4zm5kllo1hGAYWzYJNtS0a+ipLMo32RhrtjWQLWS7HLjMyO0I4HeZq7CpXY1exqJZyboTL7CoX7TNMBvFcnPOR85jnHL7XEtfQVA1Zkq8/LGZUqxV/fT0+XSKZLZBI5piZjRNJpSkUcpiMHHY5j5U0aiGFVMgi55JIeo7tGPzyTol/vVPmrajGty+rfPcKTCYNvvlWlm++leW+OoVP7dnHx5ruw7IHzl0Oc3eTF1XPoGbCKOkQarr0c+55Zu5nNopMATU7gzUzM186rmuHAYYkkdfcZM1e0mYvebOXvNVHweydE5A6FLsPk9mKIptQZQlDARgna29E1lSUxbTaMOYOohefG3rla4y57QXQDchnFhlfel5FhEotQy6/Wpy9lJa8NFvRx6JZQDFdFwxZu/5azFJuW4QwvEdRZAWXNlebybGNTCFDPBsnno0znZomlo2RK+SQZRmbasOqWqsm0ZkUE82eZpo9zcxmZhmdHeVS7BLpfJpzoXOcC53Db/WXM7BNiqmcixFPFduTvht+F1mRMQwDSZKKP5GQJAlZKtZOUuaWjmSHjG4xyOQMQmmdZNYglwNFcuAwyThNCiYKKHoB2cij5rM0ain+d0+K/+XePKeuybx8RePMlMZb0wXemi7wl2dSfGirRKNZ4pqewWZSsKj1xYddwuKWsKhgUSTMKsiSBHp+TijCqOkwanoGJR0uCkpqBjUTQclEkI0CSn4Wc34W4qPzJhyVs5C8Zidn8pI1e8maPNyVUwnHzCiygiwrqKqCoqgoioKizImnoiArKqqizF2E5YV3/qW7//J25YZlpjnHetVlKAV0A0tuFszNoKnFWU0hB7npuZnOvGZTklzpdNesRfEwWYtiUXLQK9r117Ioob4ZEcIgAIr+CLPVjN/qL4euxnNxZjOzRNIRZlIz6IaOSTVhU22YFfOCchxus5uHGh7iwfoHuZa4xujsKOOJcWZSM8ykZnh96nW2O7az27WbBltDOfR1i20L8g3RSIZhYGCgGzq6oVc8t8oGZk3HaYNMrkAikyOeSTCZzHIxlUfXDcyKgtlUXP+XTDKGrqGYFOx7Cvyr3VnaUmlev6Ry5rLGREzje1cBFAhmlz1XJgXMClhUMxZ1K1ZlGxZVwqxeFxCLA2xuHZ8Ux2dE8BlhPEYYdyGCMx/GXohgy4WwZMOoRg41n0TLJ7Enr2IA/oKOEiqe3xuXsMpIxc05irPB0m+j+GuR5j2vOheoCRX4WLaAGpyLxDK7wOIq/jQ7K1+b7KDZi5nlmgXyqeIS2fwuf5I0JwjqdZ+Kai1+piQYylyocGnmsURWv2BtEMIgWIAsyeX8hEZ7Izk9RyKbIJaLMZOaIZFLEElHkCQJq2rFptkqnNiyJLPdsZ3tju2k82kuRi9yMXqR2cxsRRmOnfadZPQM4XQYRVOQkcszBFmSkZAqlpdKs4gKzMBcwnWhYJDMFohncoSTOZLZPLmcgabIWDQZkyJhUBQZk03n436Djz5U4EpY57XhLOGpAikVsnmDXMEgm4dcQSJbkMgVrgtXtlB8xLJQvDQv1X7UOvfYusj7Bk5SNEhhGqQIW6QIW+UITimNwwQ21cCqGNgUA6uqY1EMLAqYZR2zpKPJOpKhF4VUL4bZShSQDB0Zg+I8QkeWDBQMZMlARkcyDCTJQDZ0JAwkQ59bYpr/KPpA9FxxZkc+A/kpSEzV9ockSWByFGtVlUTE7CyKh9lRnFFo1uKylGYBzVkUH8m4Lh6yWtymWSk75WMTYLbMW7oSfo/VRgiDYFk0WbueI+HYQSqfIpFLEM1ECWVCSzqxLaqFu313c5f3LsLpMCPREa7ErhTLcITeJpwJE7wcRKqxbediwlF6Xn7PkMjrkDMMUgmDbN4gr4MiSZgUBbOmosyNl0wyD98NsjuKss2PrKgokoIslX4qSOiQN8jnCuTzOrmcXvyZ1cln82TzBTI5nWzeIJM3yOoy2bxEVlfIFpTia10lV5DJ6QrZgjwnOjKpgpVh3cawUewNTukGe/nJCwAODTwWGY9ZwmuR8Fgk3GYJt0nCaSq+79SKP62qUT7XsiShKhKqJJXF06wqaKqMqkhosgwYvPyzMR6734s5n0DORiETg0wU0tF5P+e2lV7nkkVhycSKj1pRTHMi4iyKiskxNwOxIalW6qMueOca2H1FcVFNc2HB2nWRUc1V/B5zsxHh96gJIQyCFSFJUjmZrd5Wz259N4l8oroTW7Vg04pObEmS8Fl9+Ky+62U4IqOkpXSxCJ/MgmUjw1h4J64bOhhQoFDFukVs1kBVDIy8TiavM5spkE8Wl6pUWUaVZWQMyCcgNg219pZW5x42kFFQJQkzYDOk4t05EqphoBo6ilFA0Q1kCki6MRdLZCBLEoahoEsyeTQKukZO15BjTcxqUVJZE4msSjKrksypxDMKiaxCPKsQz8rohlTMAcnpXKnh+qtI4DHPPUwGbpOB2wQuTcehGbg0HbdJx20qvm9RJJA0zo+Fij4mxYEsuZHMMopVRZLlOb+HgiwVQ3YlCSQ9j5JLoGRjKNkYci6GnC0+pEwUORtFysSKj2wUKR0DIw+FLFJyBpIzC20HWrIF1GtzS0uqGWx+sHqLpU0sbrC6wewBm69YC0s1zznCS0tX5uKMZb7TfIHvQyxdCWEQ3BKl4n3VnNgz6ZmqTmxVVtnp2skO2w62XdpGYE9ggY8Bin6GG/0L80VDN3R09IrXFWMNA52F7+UKOqlsjkQ2RzSVJZ3Pk8/lcapxYpqVgmSgGwV0o4Bh6OXnBea2oVOY23Yd/fo6fnld/waBkQBDQpavL9dIuo5kFJD0PJJRKC/h6OoYqjqLZpLYImmYJA1NUtFkFU1SUWUVxVAo6BrprIVM3kIsqxZFJGcikVNJZLU5AVGJZxWSeZmCATPp4qNoUMnG6k5gqwoW2cBlVrGpxRmHXc1hVcCmpLAqOjZFL/5UdSwq2BQdmwp2VZ/zvViRsGHI28CmgF1GKjnDZQVJlpGQkY0cWi6GmkugFWJouThaLo6Si6Hm4iiZWVKRSbxSUWikXBpmr8Ls2PXTfqNPxeKeJx5esHrmlrbcxddmZ9EBPn/pSjEVZymq9bqw3CgeN960GPOju258vkjU2KJRZYt8VteLEWaGPue3mfsZjlT93d0KGyIMwWCQgYEBAoEAwWCQw4cP4/F4qo4dGhoCoKWlhWAwSCQSoaWlZR2tFayElTixLSzfcEiR1vbuzdAhmcsTT2aJvTmGUd8AsjznOSj+V/q3CnMeBcOY+3ddFKCCUaBAUUAK+pyIGPmyeBTmiUpBL4pNxXtzz3WjgK5n0QtZ1GiSuMmHQTFKC0m67oQ28kh6oRhshATmKKpVw6JYcChmfLIJk2zCpJjRVBOapKHIGrqukMmqJDMS8YxBMqOQzEqk0hKJDCQzEok0xDMGsZROXodUHlJIhMvLWhLFe3fm/VwaWSoKjF0Dm8qcaBT9JjY1h00tikvRjyJjUx3YFAdWRS8Kk6n4eVWWGasz0+gs+kW0XAwtF8WcDmHOzmDOhDBlw5jSIbTMDLKeg2QEKRkpmy6Vv8PcJlnBsNUh2f1IJfGweOZmHi6weIuzC5iX86FRFtLLr4I8/2JO5UWfaiJAcTtGMYxYL86UyhFfhQzoOciXtmXntuWLr/Xc9bF6DjkyW+Nfe+1siDB0dnZy5swZoCgShw4dor+/v+rY3t5e+vr6AGhra1t0nGDzsZwTO5Yurn1MJCfK4arVKN15G4s4eUsO6Wrvz79rr7b/8r7l4ntOR3IuN2wu5Eea+5wEkiFhSMb1n/Oc4bKkgqEiSzKGwbxEQal8MSqOlZCM658rb5Ok67YU4NqZSzQ8tIO0niGZT5PMJUnm06TyKVLZFMl8ilS++DNbyKLrkNBTxAup4nWIOfGq+LYGqmTCJFswmyzYrWa8kgmTYkJTVEyyCUXS0BQVCcgVIJkC05jBlFcmV5BJ5yCbl8nkIJuDTF4ik4NMziCTk0jnDTJZSOcM0jmjeN0zIJErPuaf+esX6BoFBjApBg5NxaYa2FQfNtWLVdmJVTWwKhRnKE4Ji1fGJafw6RG8RhiPHsKVD+PMz2DPhbBmw5hykeJFNnoNKXrtumVSpXwYJhu61Y9h82FYfRhWD7rJxdZInsK7RnEJTM8i5YsXaql8Mc8hFeZf3EuvM8Vt+Uxx+7wzcjPIyfxNfnJx1l0YgsFgxetAIMDg4OCi4/fv3084HAZYdFYhuD240YkdS8X4Lt/lgboHUNXa/xQXLNHc+P4SDsbFPpvP5/kRP+LB+gdRFGXuwnp9CaoUPltajsKgeKePXpwloGPoRnFmoBdnBED5Pb18Nwm6VEA3rs9GSseAueWzuWzjRD6MpEhYNLCZ7ciSE1VSUWQFRVKKy0mSgm7opAtpktkUyVxRLBK54vPUnICkcinyRmHuO6XJ6ymyheKFu5ApHrtoUxFVMqFJFkySCY8pS0GzI5tl7Mi4JBlJllAoLgkpUlEIFUlClmWUUiY2UChI5PKQLyhk8jLZnEwuX3xkczKZXFFwiuIikc1DOsfc66LgpHNzN9ZAuiCRLhR/k8tfSm1zj21V37XIWbaps2yXQ2yTQ2yVQ2yRwtQTpsEI4Tci2Iw0UjYKiSgSI8Wjzh32fmQKV3R0qegvKv91VTFrwaYqy42GrKDLKoZiQlc0DEVDV0wYikpB0cgrJnRJo6Cq6IqJvKyhqxqz0QRwYplzsTLWXRgGBwfx+XwV23w+H0NDQ4suEQlBuPOQJAnrXLtNv8WPpm1suGEuV7yd9Zq9a2LLUiJTfj3nE8nlcrzCKzxQ/wDIkNfzZPUsmXyGTKH4yBVyZAoZCnpxKYq5WYLNZMFutrL1BgGRkDAwyBQypPIp0nOzj3Q+TaqQIpVLk8glSefTFPSSnWkKeorxwiwWxYUhzV2gdYNCwUDXi8JSnp1UPDfmrqAGGHMOaebyLUwgmYoRYjZAkYtRZYosIyOXZ0+SJCFJ8px4SuRzEupMgZhLplCQyenF6iF5HfKF4gwnX5DIFQzycyHH+cL194oPGX1O/fLAJSQuAeCce+wqvjmnOyZyeEjglWK4SeAhjleK45YSGEAWlSwKWRTykkJeUinIMrokU5AUCrKCIUvosgKyhC7LoMwlIsoSqEWRlVUJRQZV0VFlHUUykOTScpQxF4Zs3JCrYiAZORJaxVRsVVh3YYhEIlW3h0KhRccPDAwAcOrUKbq6uggEAmtlnkCwJqzEX5KTiv/QFxNM3SjOUnJ6jryRJ6/nKegF8nrxebqQLotIXs8XL/Zl/0dxNmJVrVhUS3nWocgKqlQUkIJRIFvIFmchmSSXZi7RuLURQzbKPpWSg1/XdXJ6nrxeIG8UyBcKFHT9+nNDn7NLJ68XyBVKdlDMKDGgkC+Uny84b8wlYs9Jm0lJIalWNE1Gu2FcpeNZuuH1vKUhg6JthQKFgk6+oFMoGMWfeX3OVoOCLpEvGOQNiZguEdZt5HUHeUMirxcfBf3G5brSQSi7GG4GCVBkA0UyUGSjWJxRnnsuF6PLVKWoMYXU6l/GN01U0mKCMd8xHQgEaG9vZ3h4eMG4TCZDJpMpv56dLTpkFhOc9SSXy5FMJpmZmdkUd8bCljvPFhkZ09x/tmL8LMiUZxR5PV/+WXqe03Nks1kyhQxJPXndka7POcklA6NgUJetQ0tryErl3byMXPTHKHPPpaIdEvPv+KWK13J5BqCTKxSKF2GjQLaQL4pIoUBeLwpKLj9nY2FOVPI51FyUjOTEkCozwvViVAClTcb853NPrjvvjbmsceN6yPF8DAMoIOkFJKNYVl3SC8jkkAtZJKOArOuoeT+6GqKARMEo5qfkC3P5KoZKTpcpGDLZgkpWl8mX8lnyUjmPJZMvLp9lCxKZgkQ2NzcvMEpLZTc8r3gNIKFnUsB/WtRHdzOsuzB4PJ4FF+tQKLToclEwGCwvMZWimILB4IJZw7Fjx3j66acXfP6uu+5aHcMFAoFgEzMzM4Pb7V6VfUnGaspMDQSDwYqoJACv18vIyMgCcRgaGuLgwYNl53MkEsHr9RIOhxeMvXHGEIlE2LVrF5cuXVq1k3WzRKNRmpqauHz5Mi7X8o1yhC3CFmGLsKVWZmdn2blzZ9Xr4s2y7jOGG+/0g8Egra2t5S80NDSEx+MhEAgQCATo6ekpjx0cHKSjo6PqlzebzZjN5gXb3W73hv/iSrhcLmFLFYQt1RG2VEfYUh15FSvVboiPob+/n+7ubg4cOMCpU6cqchOOHTvGgQMHOHLkCB6Ph9bWVo4fP47H42F4eFjkMQgEAsEasyHCMH8m0NHRUfHejRf+lpYWkeksEAgE68gd2yXDbDbz5S9/uerykrBF2CJsEbYIWxZn3Z3PAoFAINjc3LEzBoFAIBDcHEIYBAKBQFDBpsl8vhWWKs29khLft3L8Q4cOVeRmLHfstbJrMVs24hwNDQ2VCySeOnWK5557rqbvvxb2LGXLep+bkh2RSIRTp07x2c9+tqbjrbYtS9mxkf+muru7OXr06Ib9rSxly3qfl5s93i3bYtwBHD58uFxtqq2tzQiHw+X3Wlpays+Hh4eNjo6OVT12f3+/cebMGaPaqVzq2Gth11K2bMQ56unpqXg+/zjrfW6WsmW9z43H4zHOnDljGIZh9Pb2GoFAoKbjrbYtS9mxUf+mSn+/tR5vvW1Z7/Nys8e7VVvuCGHo7e01wuFwxUkzjOIJmX+CDKP4j2EtuPFivNSx19quasKw3ufozJkzFfsZHh42AGN4eHjdz81SthjG+p+bkydPlp/39vaWj7He52UxO0qvN+LfVH9/vxEIBMrH3ch/RzfaYhjrf15u5nirYcsd42PweDwLpkpLlfhea5Y69kbZtZ7nqKWlheeee678ulQk0efzrfu5WcqWEut5btra2srP+/v76erqWvZ4a2HLYnaUWO9/UwMDAwvymjbq31E1W0qs93lZ6fFWw5Y7wsewWGnulZb4Xm2bFjv2Rti1Eedo/j+sEydO0NbWhsfj2ZBzs5gtsDHnZmhoiBMnTtDe3s7hw4fLdix2vLWypZodJVvW85xEIpGqa+AbcU4Ws6X03nqfl5UebzVsuSOEodbS3CUWO3HrwVLHXku7NvIclf64b3SIr+SYq2VPNVs24ty0tLQQCATo7u5e8u50uePdqi2L2bHe5+T555+vEKblWMtzspQt631eVvN4K7HljlhKmt8udH5p7pWW+F5Nljr2Rti1keeou7ubkydPlve5kefmRltg486Nx+Ohs7OTzs7O8l3qRpyXG+2A9T0ng4ODPP7444vatp7nZClbYP3/Vm7meKthy20vDKXS3Dfi8/kq1lDn09rautZmLXns9bZrI8/R8ePH6e7uLk9/I5HIhp2baras97kZHBzE6/WWX5eqDQeDwXU9L0vZsRF/L88//zx9fX309fURDAY5duwYQ0NDG/K3spgt631ebvZ4q2HLbb+UtFRp7hsV8sYS36vN/LXJpcqLr4ddN9qyEedoYGCgvFQRiUTKU/SljrlW9ixmy3qfmxv/UZfKzFcrFLmW52UpOyKRyLqekxsvZF1dXYu28F3rv5WlbFnv83Kzf5urYcsdUSup5In3zJXmnn8yg8Egvb295RLf85NVVoPBwUFOnjzJ8ePHOXLkCAcOHCiv0y517LWwaylb1vscBYNBmpubK7Z5PJ5y06X1PDfL2bLe52ZgYKA81T958iQ9PT0Vd+zrdV6WsmMj/k1FIhH6+vro7u7m8OHDdHV1lRO71vPf0VK2rPd5udnj3aotd4QwCAQCgWD1uO19DAKBQCBYXYQwCAQCgaACIQwCgUAgqEAIg0AgEAgqEMIgEAgEggqEMAgEAoGgAiEMAsEaUMqsXg/Wo1qw4L2FEAbBHUEkEqGzsxOv14vX66W9vb2izszN4PV6b2ofJVtu5Pjx4+zfvx+v10tzc/OCMtc3UqpfVI3BwUH2798PQG9vb7kCp0CwGghhENwR7N+/n/b2dsLhMOFwmK6urg2rotvZ2Ul3d3dFpmlnZycnTpygv7+fcDi8bKVZgKNHjy56we/t7S0LS29vL93d3RtaNVhwZyEynwW3PcFgkP3795dLXKwWXq+XM2fOVK3Zs5QtN5ZGHhgY4NChQzdlX3Nzc7ksw422jYyMlMWnu7sbv9/PkSNHVnwMgeBGxIxBcNvj8/kqGppUo3TBLi3jlMZ2d3fT3NxcvgAv9fn9+/cvO65ab4Vjx45x9OjRJb/DYvvv6uqit7e3YmxfX19Fo6HFxgkEN82KGoEKBJuU/v5+AzA8Ho9x+PDhcpP7Eh6Px+jv7zcMwzDC4XD5/dK20pj5fZA9Hk+5L7TH4yl/pqOjo2LcfNra2ir2aRjFHtyLjZ9/rGr7v7E/dekY1fYn/jkLVgsxYxDcEXR0dBAOh3nuuecIhULs37+fwcFBoHgX7/P5ynfy88tcz7+7f/zxx6tG+AwMDNDa2lr+zNGjR+nv769qRzAYrLr0dGMP3lr3HwgEaGtrK89wIpEIp0+frlpz3+Px3LLDXSAAsZQkuIPweDx0dHTQ399PT09PuUTxUk1wIpEIx48fp6uri8HBQWZmZhaMKXXN2r9/P/v37+fQoUOLljCu1imr1HmrRGn5SpIkBgcHl93//GWivr6+RdtOlpbUBIJbRQiD4I6kpaWl3GvgxgtzidLFuNQQZTHx8Hg8tLW1cebMmfJjfl38+VQ7VldXF8eOHSu/7unpYXh4uDxDWG7/HR0d5c5q86ORBIK1QgiD4LZncHCQrq6u8gU5EonQ3d3NZz/7WaB4YT19+nR5aankqC71zi11xTp9+nTV/T/++OMMDg6Wl5kikUh5XzdS7a69FCk0v6dyaT+17v/w4cPlENjFoqQWW8YSCFaKEAbBbU9bWxvNzc20t7cjSRL79+/ns5/9bEXo5pkzZ+ju7sbr9ZYTw0ozBK/XS2dnJ4FAAL/fv2D/Ho+H/v5+Dh06RHNzc9U+vCVaWlo4derUgu2lsNeDBw+WI6Pa2trKLReX239pqWux2cLQ0BCBQGBVuxMK3ruIPAaBYBWplsewHhw/fpyZmZlFl7gEgpWgbrQBAsGdRCAQIBAIMDg4uKjPYi3o7e2tKZtaIKgFMWMQCFaZSCTCwYMH1+1C3d3dzYEDBxYk1gkEN4sQBoFgDYhEIoRCoXVxBg8NDZUjnASC1UAIg0AgEAgqEFFJAoFAIKhACINAIBAIKhDCIBAIBIIKhDAIBAKBoAIhDAKBQCCoQAiDQCAQCCoQwiAQCASCCoQwCAQCgaACIQwCgUAgqOD/AXT0NRhZhKl8AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# df_acc['BP3']*=0.99\n",
    "# df_acc['BP2']*=0.99\n",
    "\n",
    "import seaborn as sns\n",
    "\n",
    "plt.figure(figsize=(4, 3)) \n",
    "\n",
    "sns.lineplot(data=df_acc, x=\"scale\", y=\"BP1\",label=\"BP1\")\n",
    "sns.lineplot(data=df_acc, x=\"scale\", y=\"BP2\",label=\"BP2\",alpha=0.9)\n",
    "sns.lineplot(data=df_acc, x=\"scale\", y=\"BP3\",label=\"BP3\",alpha=.5)\n",
    "\n",
    "plt.grid(True)\n",
    "plt.ylabel('S-B Classfication Accuracy')\n",
    "plt.xlabel('Scale (GeV)')\n",
    "plt.xticks(scales)\n",
    "plt.title(r'MET $>150$ GeV')\n",
    "plt.ylim([0.5,1.05])\n",
    "plt.xlim([50,500])\n",
    "plt.tight_layout()\n",
    "\n",
    "# Set custom ticks for the x-axis\n",
    "plt.savefig(f'plots/acc_150_R1_met150.pdf',bbox_inches='tight')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import seaborn as sns\n",
    "\n",
    "# plt.figure(figsize=(4, 3)) \n",
    "\n",
    "# sns.lineplot(data=df_acc, x=\"scale\", y=\"BP1\",label=\"BP1\")\n",
    "# sns.lineplot(data=df_acc, x=\"scale\", y=\"BP2\",label=\"BP2\")\n",
    "# sns.lineplot(data=df_acc, x=\"scale\", y=\"BP3\",label=\"BP3\")\n",
    "\n",
    "# plt.grid(True)\n",
    "# plt.ylabel('S-B Classfication Accuracy')\n",
    "# plt.xlabel('Scale (GeV)')\n",
    "# plt.xticks(scales)\n",
    "# plt.title(r'MET $>50$ GeV')\n",
    "# plt.ylim([0.5,1.05])\n",
    "# plt.xlim([25,500])\n",
    "# plt.tight_layout()\n",
    "\n",
    "# # Set custom ticks for the x-axis\n",
    "# plt.savefig(f'plots/acc_150_R1.pdf',bbox_inches='tight')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_acc.to_csv('df_acc_150_pT50_met150',index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "psy",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}