{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/hep/davide/miniconda3/envs/root_env/lib/ROOT.py:301: RuntimeWarning: numpy.dtype size changed, may indicate binary incompatibility\n",
" return _orig_ihook( name, *args, **kwds )\n"
]
}
],
"source": [
"import ROOT as r\n",
"import ctypes\n",
"import numpy as np\n",
"from array import array\n",
"import root_numpy as rn\n",
"import matplotlib.pyplot as plt\n",
"import pickle\n",
"import math\n",
"\n",
"from tools.data_processing import *\n",
"from tools.kinematics import *"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"masses = {\n",
" 'Ds': 1.968,\n",
" 'Dplus':1.869,\n",
" 'pi':0.140,\n",
" 'mu_plus': 0.105,\n",
" 'mu_minus': 0.105,\n",
" 'e_plus': 0.0005,\n",
" 'e_minus': 0.0005,\n",
"}"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"l_index =0\n",
"meson_index=0\n",
"mother_index = None\n",
"data_index = None "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# RETRIEVE DESIRED DATA"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"mother_ID=['Ds','Dplus']\n",
"meson_ID =['pi','X']\n",
"l_flv = ['e','mu']\n",
"\n",
"data_type = ['MC','data']\n",
"\n",
"def find_file_path(l_index=l_index, meson_index=meson_index, mother_index=mother_index, data_index=data_index): \n",
" if data_index == 0:\n",
" return \"/disk/lhcb_data/davide/Rphipi/\"+data_type[data_index]+\"/\"+mother_ID[mother_index]+\"_phipi_\"+l_flv[l_index]+l_flv[l_index]+\"/\"+mother_ID[mother_index]+\"_\"+l_flv[l_index]+l_flv[l_index]+meson_ID[meson_index]+\".root\"\n",
" else:\n",
" return \"/disk/lhcb_data/davide/Rphipi/\"+data_type[data_index]+\"/\"+mother_ID[mother_index]+\"_phipi_\"+l_flv[l_index]+l_flv[l_index]+\"/\"+mother_ID[mother_index]+\"_phi\"+meson_ID[meson_index]+\"_\"+l_flv[l_index]+l_flv[l_index]+\".root\"\n",
"mother_index=0\n",
"\n",
"data = r.TFile(find_file_path(l_index=l_index, mother_index=mother_index, data_index=1))\n",
"MC_Ds = r.TFile(find_file_path(l_index=l_index, mother_index=mother_index, data_index=0))\n",
"\n",
"tree_name_Ds = mother_ID[mother_index]+'_OfflineTree/DecayTree'\n",
"\n",
"t_data = data.Get(tree_name_Ds)\n",
"t_MC_Ds = MC_Ds.Get(tree_name_Ds)\n",
"\n",
"\n",
"mother_index=1\n",
"\n",
"MC_Dplus = r.TFile(find_file_path(l_index=l_index, meson_index=meson_index, mother_index=mother_index, data_index=0))\n",
"\n",
"tree_name_Dplus = mother_ID[mother_index]+'_OfflineTree/DecayTree'\n",
"t_MC_Dplus = MC_Dplus.Get(tree_name_Dplus)\n",
"\n",
"#Switch on only the branches that you need\n",
"t_data.SetBranchStatus(\"*\",0)\n",
"t_MC_Dplus.SetBranchStatus(\"*\",0)\n",
"t_MC_Ds.SetBranchStatus(\"*\",0)\n",
"\n",
"\n",
"for branch in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0):\n",
" t_MC_Ds.SetBranchStatus(branch, 1)\n",
"\n",
"for branch in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0):\n",
" t_MC_Dplus.SetBranchStatus(branch, 1)\n",
" \n",
"for branch in return_branches(data_index=1, mother_index=1, l_index=l_index, meson_index=0): \n",
" t_data.SetBranchStatus(branch, 1)\n",
"#Create a dictionary\n",
"\n",
"#dict ={'branch_name'=[branch_value[event]]}\n",
"\n",
"MC_Ds_tuple_dict = {}\n",
"branches_needed=return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0)\n",
"for branch in branches_needed:\n",
" \n",
" MC_Ds_tuple_dict[branch] = rn.root2array(\n",
" \n",
" filenames=find_file_path(l_index, mother_index=0, data_index=0),\n",
" treename = tree_name_Ds,\n",
" branches = branch,\n",
" start=0,\n",
" stop=t_MC_Ds.GetEntries(),\n",
" )\n",
" \n",
"MC_Dplus_tuple_dict = {}\n",
"branches_needed=return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0)\n",
"for branch in branches_needed:\n",
" \n",
" MC_Dplus_tuple_dict[branch] = rn.root2array(\n",
" \n",
" filenames=find_file_path(l_index, mother_index=1, data_index=0),\n",
" treename = tree_name_Dplus,\n",
" branches = branch,\n",
" start=0,\n",
" stop=t_MC_Dplus.GetEntries(),\n",
" )\n",
" \n",
"data_tuple_dict = {}\n",
"\n",
"branches_needed=return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0)\n",
"for branch in branches_needed:\n",
" \n",
" data_tuple_dict[branch] = rn.root2array(\n",
" \n",
" filenames=find_file_path(l_index,mother_index=0, data_index=1),\n",
" treename = tree_name_Ds,\n",
" branches = branch,\n",
" start=0,\n",
" stop=t_data.GetEntries(),\n",
" )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HLT and PID Preselection"
]
},
{
"cell_type": "code",
"execution_count": 101,
"metadata": {},
"outputs": [],
"source": [
"#HLT1 PRESELECTION\n",
"data_tuple_dict_presel_1={}\n",
"MC_Ds_tuple_dict_presel_1={}\n",
"MC_Dplus_tuple_dict_presel_1={}\n",
"\n",
"\n",
"for label in return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0): \n",
" data_tuple_dict_presel_1[label] = data_tuple_dict[label][data_tuple_dict[\"Ds_Hlt1TrackMVADecision_TOS\"]]\n",
" \n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_1[label] = MC_Ds_tuple_dict[label][MC_Ds_tuple_dict[\"Ds_Hlt1TrackMVADecision_TOS\"]]\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_1[label] = MC_Dplus_tuple_dict[label][MC_Dplus_tuple_dict[\"Dplus_Hlt1TrackMVADecision_TOS\"]]\n",
"\n",
"#RareCharm D2pi l l HLT2 PRESELECTION\n",
"\n",
"data_tuple_dict_presel_2={}\n",
"MC_Ds_tuple_dict_presel_2={}\n",
"MC_Dplus_tuple_dict_presel_2={}\n",
"\n",
"for label in return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0):\n",
" data_tuple_dict_presel_2[label] = data_tuple_dict_presel_1[label][data_tuple_dict_presel_1[\"Ds_Hlt2RareCharmD2Pi\"+l_flv[l_index].capitalize()+l_flv[l_index].capitalize()+\"OSDecision_TOS\"]]\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0):\n",
" MC_Ds_tuple_dict_presel_2[label] = MC_Ds_tuple_dict_presel_1[label][MC_Ds_tuple_dict_presel_1[\"Ds_Hlt2RareCharmD2Pi\"+l_flv[l_index].capitalize()+l_flv[l_index].capitalize()+\"OSDecision_TOS\"]]\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_2[label] = MC_Dplus_tuple_dict_presel_1[label][MC_Dplus_tuple_dict_presel_1[\"Dplus_Hlt2RareCharmD2Pi\"+l_flv[l_index].capitalize()+l_flv[l_index].capitalize()+\"OSDecision_TOS\"]]\n",
"\n",
"#PID preselection\n",
"\n",
"#MC_PID_indices=np.where(MC_tuple_dict_presel_2[l_flv[l_index]+\"_plus_MC15TuneV1_ProbNN\"+l_flv[l_index]]>0.4)\n",
"\n",
"data_PID_indices_plus=np.where(data_tuple_dict_presel_2[l_flv[l_index]+\"_plus_MC15TuneV1_ProbNN\"+l_flv[l_index]]>0.4)\n",
"data_PID_indices_minus=np.where(data_tuple_dict_presel_2[l_flv[l_index]+\"_minus_MC15TuneV1_ProbNN\"+l_flv[l_index]]>0.4)\n",
"data_PID_indices_pi=np.where(data_tuple_dict_presel_2[\"pi_MC15TuneV1_ProbNNpi\"]>0.4)\n",
"\n",
"data_PID_indices = np.intersect1d(data_PID_indices_plus,data_PID_indices_minus)\n",
"data_PID_indices = np.intersect1d(data_PID_indices,data_PID_indices_pi)\n",
"data_tuple_dict_presel_3={}\n",
"MC_Ds_tuple_dict_presel_3={}\n",
"MC_Dplus_tuple_dict_presel_3={}\n",
"\n",
"for label in return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0):\n",
" data_tuple_dict_presel_3[label] = data_tuple_dict_presel_2[label][data_PID_indices]\n",
" \n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_3[label] = MC_Ds_tuple_dict_presel_2[label]#[MC_PID_indices]\n",
" \n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0):\n",
" MC_Dplus_tuple_dict_presel_3[label] = MC_Dplus_tuple_dict_presel_2[label]#[MC_PID_indices]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# BKGCAT and True ID preselection for MC Ds"
]
},
{
"cell_type": "code",
"execution_count": 102,
"metadata": {},
"outputs": [],
"source": [
"#Getting rid of MC Background using BKG_CAT in MC Ds\n",
"#Keeping only low mass and ghost background\n",
"MC_Ds_indices=np.where(MC_Ds_tuple_dict_presel_3[\"Ds_BKGCAT\"]<65)[0]\n",
"\n",
"MC_Ds_tuple_dict_presel_4={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_4[label] = MC_Ds_tuple_dict_presel_3[label][MC_Ds_indices]\n",
" \n",
"#No reflection\n",
"MC_Ds_indices=np.where(MC_Ds_tuple_dict_presel_4[\"Ds_BKGCAT\"]!=30)\n",
"\n",
"MC_Ds_tuple_dict_presel_5={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_5[label] = MC_Ds_tuple_dict_presel_4[label][MC_Ds_indices]\n",
" \n",
"\n",
"#No partially reconstructed\n",
"MC_Ds_indices=np.where(MC_Ds_tuple_dict_presel_5[\"Ds_BKGCAT\"]!=40)\n",
"\n",
"MC_Ds_tuple_dict_presel_6={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_6[label] = MC_Ds_tuple_dict_presel_5[label][MC_Ds_indices]\n",
"\n",
"#Getting rid of MC Background matching mother IDs\n",
"MC_Ds_indices=np.where(np.abs(MC_Ds_tuple_dict_presel_6[\"pi_MC_MOTHER_ID\"])==431)[0]\n",
"\n",
"MC_Ds_tuple_dict_presel_7={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_7[label] = MC_Ds_tuple_dict_presel_6[label][MC_Ds_indices]\n",
" \n",
"MC_Ds_indices=np.where(np.abs(MC_Ds_tuple_dict_presel_7[l_flv[l_index]+\"_plus_MC_MOTHER_ID\"])==333)\n",
"\n",
"MC_Ds_tuple_dict_presel_8={}\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_8[label] = MC_Ds_tuple_dict_presel_7[label][MC_Ds_indices]"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# BKGCAT and True ID preselection for MC Dplus"
]
},
{
"cell_type": "code",
"execution_count": 103,
"metadata": {},
"outputs": [],
"source": [
"#Getting rid of MC Background using BKG_CAT in MC Ds\n",
"#Keeping only low mass and ghost background\n",
"MC_Dplus_indices=np.where(MC_Dplus_tuple_dict_presel_3[\"Dplus_BKGCAT\"]<65)[0]\n",
"\n",
"MC_Dplus_tuple_dict_presel_4={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_4[label] = MC_Dplus_tuple_dict_presel_3[label][MC_Dplus_indices]\n",
" \n",
"#No reflection\n",
"MC_Dplus_indices=np.where(MC_Dplus_tuple_dict_presel_4[\"Dplus_BKGCAT\"]!=30)\n",
"\n",
"MC_Dplus_tuple_dict_presel_5={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_5[label] = MC_Dplus_tuple_dict_presel_4[label][MC_Dplus_indices]\n",
" \n",
"\n",
"#No partially reconstructed\n",
"MC_Dplus_indices=np.where(MC_Dplus_tuple_dict_presel_5[\"Dplus_BKGCAT\"]!=40)\n",
"\n",
"MC_Dplus_tuple_dict_presel_6={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_6[label] = MC_Dplus_tuple_dict_presel_5[label][MC_Dplus_indices]\n",
"\n",
"#Getting rid of MC Background matching mother IDs\n",
"MC_Dplus_indices=np.where(np.abs(MC_Dplus_tuple_dict_presel_6[\"pi_MC_MOTHER_ID\"])==411)[0]\n",
"\n",
"MC_Dplus_tuple_dict_presel_7={}\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_7[label] = MC_Dplus_tuple_dict_presel_6[label][MC_Dplus_indices]\n",
" \n",
"MC_Dplus_indices=np.where(np.abs(MC_Dplus_tuple_dict_presel_7[l_flv[l_index]+\"_plus_MC_MOTHER_ID\"])==333)\n",
"\n",
"MC_Dplus_tuple_dict_presel_8={}\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_8[label] = MC_Dplus_tuple_dict_presel_7[label][MC_Dplus_indices]"
]
},
{
"cell_type": "code",
"execution_count": 104,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABCsAAAH4CAYAAABjWWd3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl41tWZ+P/3zSYlCQQ0CBgw1ori0mqhVFTcximUDm217oJliuXCSqfW0SpOF9pBv5bqqFgr0lKXKot23LCCFSkuM1iBioLgUsctIhTQIAkImJzfH3nIL2gCD5Llgbxf15Urz+d8zjmf+3ngunJy5yyRUkKSJEmSJClXtGruACRJkiRJkmozWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskLTTImJcRKSIeLWe+69m7o/7WPm3ImJuRJRFxKaIeCUi/isiejRJ4JIkabdTa9yRIqIqIt6PiAURcVVEdPsU/f0oIk5shFAlNSCTFZI+rQ+BAyKiX+3CiPgSUJK5X7v8OuAe4P+A4cBXgOuBfwJuboJ4JUnS7msdMAA4BjgbuI/q8cSSiOi7k339CDixQaOT1ODaNHcAknZbFcDfqB4wLKxVfjYwF6gZOETEUOASYGRK6fe16j4REZOpTlxIkiTV56OU0jO1rh+NiFuAJ4HpEXFISqmymWKT1AicWSFpV0wHzoyIAMh8PzNTXtsPgb99LFEBQEqpMqU0q9EjlSRJe5SUUhnVsyQ+B/wzQGZ56bKI2BwR6yLi/ojYf2ubiHgD2Bv4Wa2lJSdm7l0REc9FxMaIKI+IxyLi803+xiQBJisk7Zr7gH2B4zLXA4GiTDkAEdGW6imbs5s8OkmStKebB3wEHJ25zgd+AZwMfJvqxMSsiNj6e8+pVC8pmUL1spIBVM8UBehC9RLVrwBnAOXAnIjo1OjvQtInuAxE0qeWUiqLiNlUL/14KvN9dkppXWayBVQPEvYC3mqeKCVJ0p4qpfRhRKyh+o8npJRGbb0XEa2B+cBKqv+w8mRK6bmI+Ago/diyElJKP/pY278A7wLfAO5s7PciaVvOrJC0q6YDp0fEXsDpfHIJyFap6UKSJEktSM1fSDInjy2MiE1Uz7hYmbnVe4edRJwYEU9GxIZM241AYTZtJTU8kxWSdtVDVE+5vArIA2Z+7P5aYBPQq4njkiRJe7iIaE/1LM5VETEQuBdYRvVyjy8D/TNV2++gnwOpXrK6HjiH6mUlXwJW7aitpMbhMhBJuySlVBERD1O9iea9KaWKj93fEhH/AwwCftwcMUqSpD3WSVT/TjMf+BbwNvDtlFICiIjuWfYzlOoZGqenlDZm2rYGOjZ4xJKy4swKSQ3hFqpnVEyq5/4NQL+I+PbHb0REq4gY3JjBSZKkPU9EFAK/BP4OzAHaAWxNVGScU0fTzUDbj5W1ozpZUVWr7BvAZxoqXkk7x5kVknZZSmke1btx13d/ZkT8FzAlIo4FHqR6h+1DgNHAG3haiCRJql+biNh64kcB0Be4EOgADE4pVUbEHODCiLgBeADoB3ynjr5eAoZExKPABuBl4HHgGuC2iPgdcCBwBVDWiO9J0nY4s0JSk0gp/TtwFnAQMBV4DPh3qgcHFzZjaJIkKfd1onqpx/9SvS/F6cBdwBEppUUAKaX7gHHAecCfgMFUz474uCsy3+cAC4C+mT5GAycAjwAjgXOpPuZUUjOIbWdJSZIkSZIkNS9nVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKW2aO4CGts8++6SSkpLmDkOSpJyzaNGiNSmlouaOoyVwPCJJUt2yHY/sccmKkpISFi5c2NxhSJKUcyLizeaOoaVwPCJJUt2yHY+4DESSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKXvcBpuS1NSqqqpYs2YNZWVlVFZWNnc4asFat25NYWEh++yzD61a+fcISdoTbNmyhdLSUj788MPmDkXKSkONR0xWSNIuKi0tJSIoKSmhbdu2RERzh6QWKKXEli1bWLVqFaWlpfTq1au5Q5IkNYDS0lIKCgooKSlxjKGc15DjEf/sIkm7qKKigv3224927do5iFCziQjatWvHfvvtR0VFRXOHI0lqIB9++CF77723YwztFhpyPGKyQpIagFPulSv8vyhJex4TFdrdNMR4xBGNJEmSJEnKKSYrJEmNYunSpfTp04e8vDwmTpzY3OE0qpKSEubMmdPcYUiSpFoqKys599xzyc/Pp3///s0dTqMaN24cw4YNa+4wGpQbbEpSIxh5+4JG7X/KiC9lXbekpIQVK1awYsUK9tlnn5ryo446isWLF/P6669TUlICwLPPPsu4ceOYP38+W7Zs4aCDDmLUqFFceOGFOx3jhAkTGDx4MNdffz0AI0aMoLi4mPHjx+90Xy3RuHHj+Pvf/85dd93V3KFIknLI0GlDG7X/mefMzLpuSUkJq1atok2bNrRt25ZDDjmE888/n1GjRmW1DOCNN97ggAMOYMuWLbRp0/C/mj7xxBM8+eSTrFq1iry8PObNm8ewYcMoLS1t8GftiRr732dHnFkhSS3AAQccwLRp02qulyxZwoYNG7apM3/+fE4++WSGDBlCaWkp69ev53e/+x3z5s37VM8sLS3lsMMO25Wwt/HRRx81a3tJkvRJM2fOZP369axcuZJf/OIX/Nd//RcjR45s7rCA6rFISUkJeXl5DdJfQ4wlPOY+eyYrJKkFGD58OHfeeWfN9R133MH555+/TZ3LLruM733ve4wZM4a8vDwigr59+zJjxow6+3zttdcYOHAgnTt3pmPHjpx22mm89957AJx88sk88cQTjBkzhvz8fCZPnszdd9/NhAkTyM/PZ+jQ6r8KvfHGGwwZMoROnTrRrVs3rrnmmpr+x40bx+mnn86wYcMoLCzk9ttv/0QMI0aMYPTo0XzlK1+hoKCA/v378+qrr9bcjwhuvvlmDj74YHr37g3A4sWLGThwIAUFBfTq1Ys77rijpv4DDzzA5z73OfLy8ujevTu//OUva+7NmDGDgw8+mIKCAo466iieffbZrD77iooKRo8eTVFREQUFBQwYMICNGzcyb948iouLt6m7dTnJ7Nmzufrqq5kxYwb5+fl84QtfAGDSpEn07NmTvLw8evXqxR/+8IesYpAkqbG1a9eOU045hQceeIA77riDpUuXAtXJjCOOOIKCggK6du3K5ZdfTkoJgOOPPx6AwsJC8vPzmT9//nbHF3W56KKL6NGjB3l5eRx++OE1yzKnTJnCqFGjmD9/Pvn5+Vx22WV89atfZcWKFeTn55Ofn8+KFSuoqqriJz/5CT169KBjx44MHTqU1atXA9XjlIhgypQplJSU8E//9E+feP7Wn+dXX301Xbt2Zd999+V3v/tdzf0RI0Zw4YUXMmTIEAoKCvjLX/7Cxo0bufDCCykqKqKwsJDzzz+/5o9I7777LoMGDSI/P5/CwkKOOeYYqqqqauKpb9y0I1OnTuWQQw4hPz+fkpISHnnkEeCTS1lrLyep699n+fLlHHPMMeTl5dGlSxdOP/30rGPYWSYrJKkFOProo/nggw9Yvnw5lZWVTJ8+fZt1jRs2bGD+/PmcdtppO9Xv+PHjWbt2LW+++Sbr1q1j7NixAMydO5eBAwfy61//mvLyckaNGsV5553Hj370I8rLy5k5cyaVlZV89atf5bjjjmPt2rUsWLCAyZMnc//999f0//DDD3PuuedSVlZW7zrMadOmMX78eMrKyjjxxBM5++yzt7k/a9YsFixYwLJlyygrK2PQoEGMGjWKdevWMWvWLC655BIWLVoEwMiRI7ntttuoqKjg5ZdfZvDgwQA8/fTTjBkzhunTp7N+/XouvfRSvvGNb7Bx48YdfkYXXXQRr7/+OsuWLeODDz7gxhtv3OHU2MGDB3PllVdy1llnUV5ezvPPP09ZWRmXXXYZc+bMoaKigkWLFtGvX78dPl+SpKZ06KGH0rt3b5566ikAOnXqxL333sv69et54okn+MMf/sD06dMBePLJJwEoKyujvLycAQMGAPWPL+py/PHH89JLL7F+/XpGjx7NGWecwYYNGxg5ciSTJk1iwIABlJeX86tf/YpZs2bRo0cPysvLKS8vp0ePHlxzzTXMmTOH5557jrVr19KzZ08uuOCCbZ7xzDPP8NJLL/Hoo4/WGcPKlStZv3497777Lg899BAXX3wxixcvrrk/Y8YMfvGLX/DBBx8wcOBAfvCDH7By5Ur+/ve/s2LFCj744IOa93jttddywAEH8P7777N27Vquu+46IiKrcVN9/vKXv3DRRRcxefJkysvL+etf/8qBBx64w3Z1/fv8+Mc/5mtf+1rNbJpLL710h/18WiYrJKmF2Dq74rHHHqNPnz7st99+Nffef/99qqqqKCoqyrq/Aw88kBNOOIFWrVrRuXNnfvjDH9b8UMvG008/TUVFBVdeeSVt2rShZ8+efPe73+Wee+6pqXPssccyZMgQANq3b19nP1//+tfp378/rVu3Zty4cSxZsmSb2RWXX345HTt2pH379jz44IP07t2b4cOH06pVKw477DBOP/10/vjHPwKQl5fH8uXLWb9+PR07dqyZ0TBlyhRGjx7NUUcdBcB5551Hx44dd/h+P/zwQ6ZOncoNN9xAUVEREUH//v3Za6+9sv6ctmrbti2tW7dm2bJlbNy4kaKiIvr06bPT/UiS1Nj22WefmtkQxx9/PIcccggAffr04dxzz93uz8+dHV+cddZZdOzYkVatWjFmzBhat27NkiVLso71t7/9LePHj2ffffelbdu2/OQnP+Hhhx/eZrnsT3/6U9q3b1/vWKR169b89Kc/pXXr1nz5y1/mm9/8Jvfee2/N/VNPPZV+/foREUQEd955J9deey2dOnWiQ4cOXH755TXjn7y8PFauXMlbb71F69atGTBgABGR1bipPrfddhujRo2qmSmx7777cvDBB2f9GdWWl5fHW2+9xbvvvku7du04+uijP1U/2TBZIUktxPDhw5k6dSq33377J5aAdO7cmVatWrFmzZqs+ystLeW0006ja9eudOrUibPPPpvy8vKdar9ixQoKCwtrvq6++mref//9mjrdu3ffYT+1l1J06NCBLl26sGrVqjr7KC0t5a9//es2z7z77rtrBlT33HMPDz30EL169eLYY4+tGRyVlpZy3XXXbdPu7bff3uHntXbtWrZs2cJnP/vZ7D6U7cjLy2Pq1KncdNNNdO/enUGDBtVMsZUkKZesWbOGLl26APDUU09x7LHH0rlzZwoLC2tmXdZnZ8cX//mf/8nnPvc5OnbsSGFhIe+9995Oj0dOPfXUmp/vffr0oW3btqxdu7amzo7GI126dOEzn/lMzXVxcXG9Y5HVq1ezadMm+vbtW/PMwYMHs27dOgAuvfRSevXqxSmnnML+++/P+PHjSSllNW6qz7vvvtsgYxGAa665hs2bN/OlL32JQw45hFtvvbVB+q2LyQpJaiH2339/DjjgAB555JFPLPfo0KEDAwYM4L777su6vyuuuIL8/HxeffVV1q1bx/Tp02vWoNYlIra57tatG71796asrKzma/369cyePXun3tc777xT83rjxo2899577LvvvnXW7datG6eccso2zywvL6/5QXv00Ufz8MMPs3r1as444wzOPPPMmnbjxo3bpt2GDRs477zzthvb3nvvTbt27Xj99dc/ca9du3bb/NWmqqpqmwHHxz8vgCFDhjB37lxWrlzJEUcc8YlpqpIkNbfly5fzyiuvcNxxxwFwzjnncO6557Jq1SrKysoYM2ZMzXihrp91OzO+mDNnDjfffDMPP/ww69ato6ysjL333rve+nU9r1u3bsyZM2ebn/EffvghPXv2zPo9v/fee9ssDS0tLa13LLL33nvTtm1bXnnllZrnrVu3rmZM0LFjRyZOnMjrr7/OrFmzmDhxIo8++ugujZt69OhR51gEPjkeqZ2kqevz6tGjB7fddhsrVqzg97//Pf/2b//GSy+9tMMYPg2TFZLUgkyZMoW5c+fWuSv2hAkT+M1vfsMtt9xS80Pr+eef/8QeEFtVVFSw1157kZ+fz6pVq7j22mu3++wuXbrw5ptv1lyfcMIJVFVVcdNNN7F582ZSSrz88ss1+0dk66GHHmLhwoVUVlby85//nMMOO4yDDjqozrqnnnoqixcv5t5776WyspKqqiqee+45XnrpJTZv3sw999xDRUUFbdq0oaCgoOaH9AUXXMAtt9zC3/72N6B6ecef//xn1q9fv93Y2rdvzznnnMMll1zCmjVrSCmxYMECNm3aRJ8+fSgvL+dPf/oTVVVVTJgwgYqKim0+r7fffrtmU61Vq1Yxa9YsNm3aRLt27Wo2QZUkKRds2bKFuXPncuqppzJs2DCOOOIIoHq8kJ+fT7t27XjuueeYOnVqTZvCwkIigjfeeKOmbGfGFxUVFbRq1YrCwkIqKyuZMGHCdjfj7NKlC++//z4ffPBBTdmoUaP48Y9/zIoVK4DqpbFbN5/MVmVlJePHj6eyspJnn32WBx98sN6NJ9u3b8/w4cO59NJLa/5IsXLlyppNLmfPnl3zeRQUFNC6dWsiYpfGTSNGjGDy5Mn8z//8D1A9pnjllVcA+MIXvsD06dP56KOPeOGFF7ZZvlLXv88DDzzAypUrger9SFq1atVo4xGTFVKumHpW3V9SAzrwwAPr3ZTxmGOO4fHHH+ehhx6iR48e5OfnM2LECE466aQ6648bN45nnnmGgoICBg8eXHPCR31GjhzJwoULKSgo4Jvf/CZt2rTh0Ucf5fHHH6+Z6jls2LDtDjLqcvbZZzN27FgKCwuZM2dOzaZddenSpQuzZ89m0qRJdO7cmS5duvCDH/yg5q8hv/3tb9lvv/3o0KEDEydO5O677waqEyu/+tWv+Pa3v01+fj69evXilltuySq+m2++meLiYg4++GA6derExRdfTFVVFZ07d+bGG29k+PDhdO/enTZt2myzpOWMM85g48aNdOrUiS9+8Ys1A6Gtp4o89thjTJo0aac+K0na3QydNrTOL+WOoUOH1pz08R//8R98//vf57bbbqu5/+tf/5qxY8fSsWNHfvKTn2wzu7NTp05ccskl9OvXj8LCQp555pmdGl8MGTKEk08+mc9+9rPsv//+ANudEXH44Yfz9a9/neLiYgoLC1mxYgX/8R//wXHHHceXv/zlmhO/nnjiiZ36DLp160aHDh3o0aMH//Iv/8J1111Xs89VXX7961/TuXNn+vTpQ0FBAccff3zNPhsvvvgiAwcOpEOHDvTr14/vfOc7fOUrX9mlcdNJJ53ExIkTGTFiBPn5+Rx99NG89tprAFx11VW8+OKLdOrUibFjx9bMKoW6/32eeuopjjzySDp06MBXv/pVfvnLX37q/S92JLY3ZXd31K9fv7Rw4cLmDkPaefUlJs6t+9hI5Y7ly5e70WEzGTFiBMXFxYwfP765Q8kp9f2fjIhFKSWPEGkCjkekhlFfYmLmOTObOJLm4zgjt82bN49hw4ZRWlra3KHknF0djzizQpIkSZIk5RSTFZIkSZIkKae0ae4AJEn6tG6//fbmDkGSJLVgJ554oktAGokzKyRJkiRJUk4xWSFJkiRJOWxPOxRBe76tR6/vCpMVkiRJkpSj2rdvz9q1a01YaLeQUmLz5s2888475OXl7VJf7lkhSZIkSTmquLiY0tJSVq9e3dyhSFlp06YNnTp1Yp999tm1fhooHkmSJElSA2vbti0HHHBAc4chNTmXgUiSGsXSpUvp06cPeXl5TJw4sbnDaVQlJSXMmTOnucOQJEnaYzizQpIaw9SzGrf/c2dkXbWkpIQVK1awYsWKbabjHXXUUSxevJjXX3+dkpISAJ599lnGjRvH/Pnz2bJlCwcddBCjRo3iwgsv3OkQJ0yYwODBg7n++usBGDFiBMXFxYwfP36n+2qJxo0bx9///nfuuuuu5g5FkiSpyTmzQpJagAMOOIBp06bVXC9ZsoQNGzZsU2f+/PmcfPLJDBkyhNLSUtavX8/vfvc75s2b96meWVpaymGHHbYrYW/jo48+atb2kiRJajomKySpBRg+fDh33nlnzfUdd9zB+eefv02dyy67jO9973uMGTOGvLw8IoK+ffsyY0bdszhee+01Bg4cSOfOnenYsSOnnXYa7733HgAnn3wyTzzxBGPGjCE/P5/Jkydz9913M2HCBPLz8xk6dCgAb7zxBkOGDKFTp05069aNa665pqb/cePGcfrppzNs2DAKCwu5/fbbPxHDiBEjGD16NF/5ylcoKCigf//+vPrqqzX3I4Kbb76Zgw8+mN69ewOwePFiBg4cSEFBAb169eKOO+6oqf/AAw/wuc99jry8PLp3784vf/nLmnszZszg4IMPpqCggKOOOopnn302q8++oqKC0aNHU1RUREFBAQMGDGDjxo3MmzeP4uLibepuXU4ye/Zsrr76ambMmEF+fj5f+MIXAJg0aRI9e/YkLy+PXr168Yc//CGrGCRJknY3JiskqQU4+uij+eCDD1i+fDmVlZVMnz6dYcOG1dzfsGED8+fP57TTTtupfsePH8/atWt58803WbduHWPHjgVg7ty5DBw4kF//+teUl5czatQozjvvPH70ox9RXl7OzJkzqays5Ktf/SrHHXcca9euZcGCBUyePJn777+/pv+HH36Yc889l7Kysm3irW3atGmMHz+esrIyTjzxRM4+++xt7s+aNYsFCxawbNkyysrKGDRoEKNGjWLdunXMmjWLSy65hEWLFgEwcuRIbrvtNioqKnj55ZcZPHgwAE8//TRjxoxh+vTprF+/nksvvZRvfOMbbNy4cYef0UUXXcTrr7/OsmXL+OCDD7jxxhtp1Wr7P34HDx7MlVdeyVlnnUV5eTnPP/88ZWVlXHbZZcyZM4eKigoWLVpEv379dvh8SZKk3ZHJCklqIbbOrnjsscfo06cP++23X829999/n6qqKoqKirLu78ADD+SEE06gVatWdO7cmR/+8Ic8+eSTWbd/+umnqaio4Morr6RNmzb07NmT7373u9xzzz01dY499liGDBkCVJ8zX5evf/3r9O/fn9atWzNu3DiWLFmyzeyKyy+/nI4dO9K+fXsefPBBevfuzfDhw2nVqhWHHXYYp59+On/84x8ByMvLY/ny5axfv56OHTvWzGiYMmUKo0eP5qijjgLgvPPOo2PHjjt8vx9++CFTp07lhhtuoKioiIigf//+7LXXXll/Tlu1bduW1q1bs2zZMjZu3EhRURF9+vTZ6X4kSZJ2ByYrJKmFGD58OFOnTuX222//xBKQzp0706pVK9asWZN1f6WlpZx22ml07dqVTp06cfbZZ1NeXr5T7VesWEFhYWHN19VXX837779fU6d79+477Kf2UooOHTrQpUsXVq1aVWcfpaWl/PWvf93mmXfffXfN8pV77rmHhx56iF69enHsscfWJCNKS0u57rrrtmn39ttv7/DzWrt2LVu2bOGzn/1sdh/KduTl5TF16lRuuukmunfvzqBBg1i6dOku97unioj2EbEwIhZHxKsRcUNU6xIRj0XEkoj4c0R0rtVmbEQsj4ilETGoVnnfiHguIpZFxMSIiOZ5V5IktRwmKySphdh///054IADeOSRRz6x3KNDhw4MGDCA++67L+v+rrjiCvLz83n11VdZt24d06dPJ6VUb/2P/37XrVs3evfuTVlZWc3X+vXrmT179k69r3feeafm9caNG3nvvffYd99966zbrVs3TjnllG2eWV5ezq233gpUL5d5+OGHWb16NWeccQZnnnlmTbtx48Zt027Dhg2cd955241t7733pl27drz++uufuNeuXbttNjmtqqraJlFT1+/DQ4YMYe7cuaxcuZIjjjiCCy64YLvPb+E2ASeklI4EDgUGACcBPwdmpZSOAGZlromIvsC3gM8Dg4FbI2LrFJjbgAtSSocC+wOnNuUbkSSpJTJZIUktyJQpU5g7dy55eXmfuDdhwgR+85vfcMstt9T8Ev38889/Yg+IrSoqKthrr73Iz89n1apVXHvttdt9dpcuXXjzzTdrrk844QSqqqq46aab2Lx5MyklXn755Zr9I7L10EMPsXDhQiorK/n5z3/OYYcdxkEHHVRn3VNPPZXFixdz7733UllZSVVVFc899xwvvfQSmzdv5p577qGiooI2bdpQUFBQkzC44IILuOWWW/jb3/4GVC/v+POf/8z69eu3G1v79u0555xzuOSSS1izZg0pJRYsWMCmTZvo06cP5eXl/OlPf6KqqooJEyZQUVGxzef19ttvU1VVBcCqVauYNWsWmzZtol27djWboKpuqdrWD7Qt0Br4B/A1YOvOpHdlrsl8n5FS2pJSKgVeBPpHRC+gdUppUR1tJElSIzFZIUktyIEHHljvpozHHHMMjz/+OA899BA9evQgPz+fESNGcNJJJ9VZf9y4cTzzzDMUFBQwePDgmhM+6jNy5EgWLlxIQUEB3/zmN2nTpg2PPvoojz/+eM1SkmHDhtUsycjW2WefzdixYyksLGTOnDlMnz693rpdunRh9uzZTJo0ic6dO9OlSxd+8IMf1GyU+dvf/pb99tuPDh06MHHiRO6++26gOrHyq1/9im9/+9vk5+fTq1cvbrnllqziu/nmmykuLubggw+mU6dOXHzxxVRVVdG5c2duvPFGhg8fTvfu3WnTps02S1rOOOMMNm7cSKdOnfjiF79IZWUl48ePrzlV5LHHHmPSpEk79Vm1NBHROiIWU52kmJdSWgoUpZRWA2S+d81ULwbertW8NFNWX7kkSWpEsb0pu7ujfv36pYULFzZ3GNLOm3pW3eXn1n1spHLH8uXL3eiwmYwYMYLi4mLGjx/f3KHklPr+T0bEopRSiztCJCIKgUeBK4CHUkoFte6tTykVRMRkYG5KaXqm/FZgHvAm8NOU0uBM+QBgXEpp0MceQ0SMAkYB9OrVq2/tmUSStm/otO0nvD9u5jkzGykSSY0t2/GIMyskSdIeLaWS4EhlAAAgAElEQVRUBvwJOBpYHRFFAJnv/8hUKwV61mpWnCmrr7yu50xOKfVLKfXbmZN1JEnSJ5mskCRJe5yI2CciCjKvPwP8M7AUeAQYlqk2jOpNNsmUnxURbSOiGDgceDal9BZQFRFfzNQ7r1YbSZLUSNo0dwBSi1Pfcg9JO+32229v7hCUu3oAd2aOGW0PTEspzYyI/wVmRMR3gFXAmQAppYURcT/wAlAFjE4pbcr09a/A7yOiHTAX+O8mfi+SJLU4JiskSdIeJ6X0AnBkHeVrgVPqaXMVcFUd5Qvr6kuSJDUel4FIUgPYeryk1Nz8vyhJkvYEJiskaRfl5eXxzjvvsHnzZva0E5a0+0gpsXnzZt555x3y8vKaOxxJkqRd4jIQSdpFxcXFrFmzhjfffJOPPvqoucNRC9amTRs6derEPvvs09yhSJIk7RKTFZK0i1q1akXXrl3p2rVrc4ciSZIk7RFcBiJJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyinNkqyIiNYR8VxEPJy57hIRj0XEkoj4c0R0rlV3bEQsj4ilETGoOeKVJEmSJElNp7lmVvwAWF7r+ufArJTSEcCszDUR0Rf4FvB5YDBwa0Ts1cSxSpIkSZKkJtTkyYqIKAa+BvyuVvHXgD9kXt+Vud5aPiOltCWlVAq8CPRvqlglSZIkSVLTa46ZFTcAPwKqapUVpZRWA2S+d82UFwNv16pXmimTJEmSJEl7qCZNVkTEvwD/SCktauB+R0XEwohYuHr16obsWpIkSZIkNbGmnllxLPD1iHgDmA6cHBF3Aasjoggg8/0fmfqlQM9a7YszZdtIKU1OKfVLKfUrKipqzPglSZIkSVIja9JkRUppbEqpOKVUApwNzE0pDQMeAYZlqg2jepNNMuVnRUTbzF4XhwPPNmXMkiRJkiSpabVp7gAyfgbMiIjvAKuAMwFSSgsj4n7gBar3uBidUtrUfGFKkiRJkqTG1mzJipTSPGBe5vVa4JR66l0FXNVkgUmSJEmSpGbVHKeBSJIkSZIk1ctkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTTFZIkqQ9TkT0jIgnI2JpRLwSEZdnysdFxDsRsTjzNaRWm7ERsTzTZlCt8r4R8VxELIuIiRERzfGeJElqSdo0dwCSJEmNYAswJqX0QkQUAH+LiEcz965PKV1bu3JE9AW+BXwe2Bd4OiIOTiltAm4D/jWltCgiHgROBe5rsnciSVIL5MwKSZK0x0kprUwpvZB5vR54AdhvO02+BsxIKW1JKZUCLwL9I6IX0DqltChT765MXUmS1IhMVkiSpD1aRJQAXwKezhRdFBEvRcTdEbF3pqwYeLtWs9JMWX3ldT1nVEQsjIiFq1evbsB3IElSy2OyQpIk7bEiIh/4I3BxSmkdcDPwOeBQ4DVgYkM9K6U0OaXUL6XUr6ioqKG6lSSpRTJZIUmS9kgR0Rb4b2BaSuk+gJTS6pRSZUqpCphE9YwLqJ4x0bNW8+JMWX3lkiSpEZmskCRJe5zMiR1TgOUppetqlXetVe1bwLLM60eAsyKibUQUA4cDz6aU3gKqIuKLmXrnAbMa/Q1IktTCeRqIJEnaEx0LDAeWRMTiTNmVwLkR8XmgHfAWMBIgpbQwIu6neiPOKmB05iQQgH8Ffh8R7YC5VM/WkCRJjchkhSRJ2uOklJ4Goo5bj2ynzVXAVXWULwSObLjoJEnSjrgMRJIkSZIk5RSTFZIkSZIkKaeYrJAkSZIkSTnFZIUkSZIkScopJiskSZIkSVJOMVkhSZIkSZJyiskKSZIkSZKUU0xWSJIkSZKknNKmuQOQJEmStGuGThtaZ/nMc2Y2cSSS1DCcWSFJkiRJknKKyQpJkiRJkpRTXAYiSZIkqUnUt1xFkj7OmRWSJEmSJCmnmKyQJEmSJEk5xWUgkiRJknYrnn4i7fmcWSFJkiRJknKKyQpJkiRJkpRTTFZIkiRJkqScYrJCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOcVkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknJKm+YOQJIkSVJ2hk4b2twhSFKTcGaFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOWUJk1WRET7iFgYEYsj4tWIuCGqdYmIxyJiSUT8OSI612ozNiKWR8TSiBjUlPFKkiRJkqSm19QzKzYBJ6SUjgQOBQYAJwE/B2allI4AZmWuiYi+wLeAzwODgVsjYq8mjlmSJEmSJDWhJk1WpGoVmcu2QGvgH8DXgD9kyu/KXJP5PiOltCWlVAq8CPRvwpAlSZIkSVITa/I9KyKidUQspjpJMS+ltBQoSimtBsh875qpXgy8Xat5aabs432OyiwvWbh69erGfQOSJEmSJKlRtWnqB6aUKoEjI6IQeDQiTmqAPicDkwH69euXdrU/SZIkSZ/e0GlDmzsESbu5ZjsNJKVUBvwJOBpYHRFFAJnv/8hUKwV61mpWnCmTJEmSJEl7qKY+DWSfiCjIvP4M8M/AUuARYFim2jCqN9kkU35WRLSNiGLgcODZpoxZkiRJkiQ1raZeBtIDuDMiAmgPTEspzYyI/wVmRMR3gFXAmQAppYURcT/wAlAFjE4pbWrimCVJkiRJUhNq0mRFSukF4Mg6ytcCp9TT5irgqkYOTZIkSZIk5Yhm27NCkiRJkiSpLiYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFOy3mAzItoChwD7ZopWAi+nlLY0RmCSJEmSJKll2m6yIiLaUH2M6HDgBGAvIDK3E7ApIp4E/gDcY+JCkiRJyh1Dpw2ts3zmOTObOBJJ2jn1LgOJiO8A/wdcD7wLXAwMBPoAh2ZeXwysAK4DXsu0kSRJkiRJ+tS2N7PiQqqTEQ+mlCrrqfO/wOSIaA18E7gC+H3DhihJkiRJklqSepMVKaUvZdtJJpnx35kvSZIkSZKkT83TQCRJ0h4nInpGxJMRsTQiXomIyzPlXSLisYhYEhF/jojOtdqMjYjlmTaDapX3jYjnImJZREyMiKjrmZIkqeHscrIiMxjo1RDBSJIkNZAtwJiU0uFAX+CCiDgS+DkwK6V0BDArc01E9AW+BXweGAzcGhF7Zfq6DbggpXQosD9wapO+E0mSWqCGmFnxf8DrDdCPJElSg0gprUwpvZB5vR54AdgP+BrVp5gB3JW5JvN9RkppS0qpFHgR6J/5g0zrlNKiOtpIkqRGst2jS7M0kv//OFNJylkjb19QZ/mUEVlv0SNpNxQRJcCXgO8ARSml1QAppdUR0TVTrRiYW6tZaaasEni7jvK6njMKGAXQq5eTTiVJ2hW7nKxIKd3ZEIFIkiQ1tIjIB/4IXJxSWteY202klCYDkwH69euXGu1BkiS1AG6wKUmS9kgR0Zbqk8qmpZTuyxSvjoiizP0i4B+Z8lKgZ63mxZmy+solSVIjympmRUTcs6M6KaUzdz0cSZKkXZc5sWMKsDyldF2tW48Aw4DrM99n1SqfFBE3APsChwPPppQ2RURVRHwxpfQ34Dyq962QJEmNKNtlIEV1lOUBRwDvAS83WESSJEm77lhgOLAkIhZnyq4EfgbMiIjvAKuAMwFSSgsj4n6qN+KsAkanlDZl2v0r8PuIaEf1vhb/3XRvQ5KklimrZEVK6aS6yiNiX+B+qv86IUmSlBNSSk9T/wbgp9TT5irgqjrKFwJHNlx0kiRpR3Zpz4qU0irg/wETGiYcSZIkSZLU0jXEBptVQEkD9CNJkiRJkpT1BpuH1lHcGuhN9XTJxXXclyRJkiRJ2mnZbrC5FKjrvPAAXgRGNlhEkiRJkiSpRcs2WVHXBpsfAe+mlP6vAeORJEmSJEktXLangTzR2IFIe5ypZzV3BJIkSZK0W2qIDTYlSZIkSZIaTLbLQOoVEa8CrVJKBzZAPJIkSZJ2E0OnDW3uECTtoXY5WQE8iTM0JEmSJElSA9nlZEVKyZNAJEmSJElSg2mImRWSJEmSdiP1Ld+Yec7MJo5Ekuq2U8mKiDgO6A20//i9lNJvGiooSZIkbWtn9wbwl06pmokZafeUVbIiInoBjwCHAgmIzK1Uq5rJCkmSJEmStMuy3RjzBmAN0J3qRMWXgB7AZcBLVM+2kCRJkiRJ2mXZLgM5GTgfWJ25jpTSSuC/IiKonlUxqBHikyRJkqSseJSqtOfIdmbFZ4A1KaUqYDOwT617C4FjGjowSZIkSZLUMmWbrPg7sF/m9RLgrFr3vglsaMigJEmSJElSy5XtMpDZwD8B9wLXAPdExImZe72AcQ0emSRJkhqcJyNIknYHWSUrUkr/Xuv1f2cSFacB7YA5KaX7Gyc8SZIkNSeTG5Kk5pDtzIptpJSeAp5q4FgkSZIkNSOTU5JyRb3JiojomVJ6e2c6i4j9U0pv7npYkiRJ2t34i+6n4wkWkvRJ25tZsSQi7gUmp5QWbK+TiDgeuIDqzTY7NmB8kiRJ+hR29hdgf2GWJOWS7SUrjgJ+BjwdEauBvwLLgPeAAPYGDgO+DHQGpgGfb9RoJUmSJEnSHq/eZEVK6XVgRERcAZwD/DMwAtg3U2Ul8ALw/4AZKaVVjRuqJEmSJElqCXa4wWZKaSVwfeZLUlObelbd5efOaNo4JEn6lNzLYve3Jy4T8v+llNs+1WkgkiRJ0p7GX14lKXeYrJAkSdJO2xP/0i5Jyh0mKyRJkqTtcMaFJDW9Vs0dgCRJkiRJUm3OrJAkSdIeyaUqkrT7+tTJiogoAj4LLE0pVTRcSJIkSWoJXF4hSapPVsmKiLgKqEwp/TRzPRi4D9gLeC8iTkkpPd94YUqSJKmlMImh5uT/Pyk3ZLtnxbnAK7WufwX8BRgALAYmNHBckiRJkiSphco2WbEf8H8AEdELOAz4aUrpWeAG4OjGCU+SJEmSJLU02e5ZsQHolHl9AlCWUlqUuf6A6uUgkiRJUqNxer4ktRzZJiueBK6ICIBLgT/VuncI8GYDxyVJkiTtljyFRJJ2XbbJiu8D04EHgeeAsbXunQs83cBxSZIkSVkxOSBJe56skhUppTep3kyzLt8ANjZYRJIkSZIkqUXL9ujSucD3Ukov1XG7GzAJOLkhA5MkSZKkXLGzM3jcS0XaNdmeBnIi0LGeewVUb7opSZIkSZK0y7JNVgCkjxdERGvgeOC9BotIkiRJkiS1aPUuA4mInwE/zVwm4JnMaSB1uamB45IkSZIkSS3U9vaseARYAwQwEbgOeONjdTYDL6WUnmqU6CRJkqQc5SkkktR46k1WpJQWAAsAImI98KeU0pqmCkySJEmSJLVM2R5dekdjByJJkiRJkgTZH13aDrgC+AbQva52KaWuDRuaJEmSJElqibJKVgC3AsOA+4E/AZWNFpEkSZIkSWrRsk1WfAu4OKV0c2MGI0mSJEmSlG2y4kNgSWMGIkmSJEl7ivpOi5l5zswmjkTaPWWbrPgtcA7wZCPGIkmSJEktkskNaVvZJitWAedFxGxgDrDhY/dTSumWBo1MkiRJkiS1SNkmK27IfO8FfKWO+wkwWSFJkiRJ21HfDApJ28oqWZFSatXYgUiSJEmSJAGYhJAkSZIkSTkl62RFRJRExA0R8deI+HtEHJYp/35EHNt4IUqSJEmSpJYkq2RFRPQHXgAGAX8DDgD2ytzuClzeKNFJkiRJkqQWJ9uZFdcDs4BDge8DUevefODLDRyXJEmSJElqobJNVnwRuCWllKg++aO2MqCwQaOSJEnaRRHx+4j4R0QsrVU2LiLeiYjFma8hte6NjYjlEbE0IgbVKu8bEc9FxLKImBgR8fFnSZKkhpVtsmIdsE899w4AVjZMOJIkSQ3mdmBwHeXXp5SOzHw9AtUJCeBbwOczbW6NiK1LXm8DLkgpHQrsD5za6JFLktTCZZuseAj4RUQU1ypLEVEAXArc1+CRSZIk7YKU0pPAe1lW/xowI6W0JaVUCrwI9I+IXkDrlNKiTL27MnUlSVIjyjZZcTnwIfAqMCdTdiPwf0Al8NOGD02SJKlRXBQRL0XE3RGxd6asGHi7Vp3STFl95Z8QEaMiYmFELFy9enVjxC1JUouRVbIipfQ+cDQwBniX6oTFCuBK4NiU0vpGi1CSJKnh3Ax8jupNw18DJjZUxymlySmlfimlfkVFRQ3VrSRJLVKbbCumlDYDUzJfkiRJu52UUs2Uh4iYBMzLXJYCPWtVLc6U1VcuSU1i6LShdZbPPGdmE0ciNa2sZlZExJMRcWFE+GcCSZK024qIrrUuvwUsy7x+BDgrItpm9ug6HHg2pfQWUBURX8zUO4/q49wlSVIjynZmxSrgWuDGiHgCmAbcn1keIkmSlHMiYhpwIrBPRJQCPwNOiojPA+2At4CRACmlhRFxP/ACUAWMTiltynT1r8DvI6IdMBf47yZ9I5IktUBZJStSSmdERB7wdeBMqtd7/iYi5gAzgAfct0KSJOWSlNI5dRTXu5w1pXQVcFUd5QuBIxswNEmStAPZngZCSqkipTQtpXQq0BUYlbn1W2BlYwQnSZIkSZJanqyTFbVlZlG8BrwOfAB8Jpt2EdEzs//F0oh4JSIuz5R3iYjHImJJRPw5IjrXajM2IpZn2gz6NPFKkiRJkqTdx04lKyKif0RcFxFvAU8CJwA3Agdl2cUWYExK6XCgL3BBRBwJ/ByYlVI6gupNq36eeV5fqje/+jwwGLg1IvbamZglSZIkSdLuJas9KyLil8AZwP7Aq8BtwIyU0rLtNvyYlNJKMktGUkrrI+IFYD/ga8CXM9XuAp4B/i1TPiOltAUojYgXgf7AUzvzXEmSJEmStPvI9jSQM4B7gOkppcUN8eCIKAG+BHwHKNp67nlKaXWtY8WKqd51e6vSTNnH+xpFZg+NXr16NUR4kiRJkiSpmWR7GshnG/KhEZEP/BG4OKW0LiJ2qb+U0mRgMkC/fv3SrkcoSZIkSZKaS7YzK4iIjlTPghgA7A1cmFJ6NSLOBJZmuyQkItpSfT75tJTSfZni1RFRlJlVUQT8I1NeCvSs1bw4UyZJkiRJLdbQaUPrLJ95zswmjkRqHFltsBkRvYFlwE+AjsBJQEHm9peBn2bZT1B9vvnylNJ1tW49AgzLvB5G9SabW8vPioi2EVEMHA48m82zJEmSJEnS7inbmRUTgZeBbwAfAptr3XsSuDnLfo4FhgNLImLr3hdXAj8DZkTEd4BVwJkAKaWFEXE/8AJQBYxOKW3K8lmSJEmSJGk3lG2yYiDwjZRSeUS0/ti91UDXOtp8QkrpaaC+DSpOqafNVcBVWcYpSZIkSfr/2rv3aMnK8s7j3580YhAvXBpw7AZyR4MEBXTwikIUJR1jUJtGcGEgqBkxOqMixvGSBWMcSUyUNVEMStTQthpBWwFFIMEkGmkRuQheEiF0AnQjoAGR6zN/7H2g+nRVdx36nKp96nw/a+1Vp956q/azd1XXfvqpd79bmueGOg2EZjTFIwY89jjgx7MTjiRJkiRJWuiGHVlxPvC2JBfRFC4Aqh1l8TqauSUkqROOOeOScYcgSZIkaQsMW6x4M/BPwL8C5wEFvJVmwstHAivmJDpJkiRJkrTgDFWsqKrrk/wm8D+Bg2iKFrsDZwN/XlWeBiJp3ho0EuP0o/cfcSSSJEmSYPiRFVTVrTSXLv3fcxeOJEmSJEla6IYuVkiSJEmS5qdlK5f1bV+9YvWII5GGM+zVQCRJkiRJkkbCYoUkSZIkSeoUixWSJEmSJKlTLFZIkiRJkqROGbpYkeSVSR47l8FIkiRJkiTNZGTFx4DdANJ4R5Jd5yYsSZIkSZK0UA28dGmSc4HLgO+0S4BqH34Y8E7gi8CNcxyjJEmSJGkIgy5ROtP+XtJU4zawWAGcBzwZeBHwBJpCxalJLgIuYcPihSRJkiRJ0qwYWKyoqr+c+jvJNsCdwKXArwNH0RQqPpHkPOCrVXXeHMcqSZIkSZIWgIFzViR5fZJnJXlUVd3VNn+sqlbQFCwCrAS2A06d+1AlSZIkSdJCsKnTQH4b+GNgpyTX0YykODzJLwBXtH3OrapL5zhGSZIkSdII9ZvLwnksNEoDR1ZU1fOrahfg8cAf0oykOJhmLotbaIoXr01yUHuaiCRJkiRJ0hbb7KVLq+rGnvkojq2q7YH9aIoXS4EzgFvnLEJJkiRJkrSgbLZYMcDV7e3bqmopsO8sxSNJkiRJkha4Tc1ZsYGq6i1sFHAdcFf72NV9nyRJkiRJkjRDQxcrelXV/cAvznIskiRJkqSO6jfpJjjxpubGQz0NRJIkSZIkaU5YrJAkSZIkSZ1isUKSJEmSJHWKxQpJkiRJktQpFiskSZIkSVKnWKyQJEmSJEmdYrFCkiRJkiR1yqJxByBJkiRJmr+WrVzWt331itUjjkSTxJEVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVO8GogkSZIkadZ5lRBtCUdWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjrFS5dKkiRJkkbGS5pqGI6skCRJkiRJnWKxQpIkSZIkdYrFCkmSJEmS1CkWKyRJ0kRK8tEk65Jc2dO2Q5Lzk1yR5CtJtu957MQkVye5MskLetr3TfLtJN9N8oEkGfW2SJK00FiskCRJk+oM4JBpbe8Gzq2qJwHntvdJsi9wGLB3+5wPJ9mmfc7HgGOr6onA7sBL5j50SZIWNosVkiRpIlXVxcAt05oPBT7R/v3J9v5U+6qquqeq1gJXAU9NshuwVVV9q89zJEnSHLFYIUmSFpLFVbUeoL3duW1fAlzf029t2zaoXZIkzSGLFZIkSbMgyXFJ1iRZs379+nGHI0nSvGaxQpIkLSTrkywGaG/Xte1rgaU9/Za0bYPaN1JVp1XVflW13+LFi2c9cEmSFhKLFZIkaSE5Bziy/ftImkk2p9qXJ9k6yRJgL+CbVfXvwP1JntL2e0XPcyRJ0hxZNO4AJEmS5kKSlcCBwE5J1gLvbJdVSX4fuAl4OUBVrUlyFnA5cD/wmqq6q32pVwEfTfJw4ELg70a6IZK0QCxbuaxv++oVq0ccibrAYoUkSZpIVbViwEMHD+h/MnByn/Y1wD6zGJokSdoMixWSJEmSpM5yxMXC5JwVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVO8Goi0pc5cPu4IJEmSJGmiOLJCkiRJkiR1iiMrJEmSJEkTY9nKZX3bV69YPeJItCUsVkiat44545JxhyBJkqQxGVSU0GTwNBBJkiRJktQpFiskSZIkSVKnWKyQJEmSJEmdYrFCkiRJkiR1ihNsSpIkSZImnlcJmV8cWSFJkiRJkjrFkRXSfHXm8v7tR6wabRySJEmSNMscWSFJkiRJkjrFYoUkSZIkSeoUixWSJEmSJKlTLFZIkiRJkqROsVghSZIkSZI6xWKFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTlk07gAkSZIkSRqXZSuX9W1fvWL1iCNRL0dWSJIkSZKkTnFkhTQTZy4fdwSSJEmSNPEsVkiSJEmSNI2nh4zXSE8DSfLRJOuSXNnTtkOS85NckeQrSbbveezEJFcnuTLJC0YZqyRJkiRJGo9Rz1lxBnDItLZ3A+dW1ZOAc9v7JNkXOAzYu33Oh5NsM7pQJUmSJEnSOIy0WFFVFwO3TGs+FPhE+/cn2/tT7auq6p6qWgtcBTx1JIFKkiRJkqSx6cLVQBZX1XqA9nbntn0JcH1Pv7VtmyRJkiRJmmATMcFmkuOA4wB22223MUcjaVIcc8YlfdtPP3r/EUciSZIkLSxdGFmxPsligPZ2Xdu+Flja029J27aRqjqtqvarqv0WL148p8FKkiRJkqS51YWRFecARwLvb2/P7Wn/UJK/AHYB9gK+OZYIJUmSJEnCS5qOykiLFUlWAgcCOyVZC7yzXVYl+X3gJuDlAFW1JslZwOXA/cBrququUcYrSZIkSdIwLGLMrpEWK6pqxYCHDh7Q/2Tg5LmLSJIkSZIkdU0X5qyQJEmSJEl6gMUKSZIkSZLUKRYrJEmSJElSp3ThaiCSJEmSJE2kfhNvOunm5jmyQpIkSZIkdYrFCkmSJEmS1CkWKyRJkiRJUqdYrJAkSZIkSZ1isUKSJEmSJHWKxQpJkiRJktQpXrpUkiRJkqQR6nc5U/CSpr0cWSFJkiRJkjrFYoUkSZIkSeoUixWSJEmSJKlTnLNCkiRJkqQOW4hzXDiyQpIkSZIkdYrFCkmSJEmS1CmeBiJJkiRJUgcMOt1jIbJYIU2aM5f3bz9i1WjjkCRJkqSHyMJgFd0AABhYSURBVNNAJEnSgpPk2iRXJLksyZq2bYck57ftX0myfU//E5NcneTKJC8YX+SSJC0MFiskSdJC9dyq2qeq9mvvvxs4t6qeBJzb3ifJvsBhwN7AIcCHk2wzjoAlSVooLFZIkiQ1DgU+0f79yfb+VPuqqrqnqtYCVwFPHUN8kiQtGBYrJEnSQlTA1Ckfx7dti6tqPUB7u3PbvgS4vue5a9s2SZI0R5xgU5IkLUQHVNWNSXYGzktyzZa+YJLjgOMAdtttty19OUmSHrJBVxVZvWL1iCN56CxWSJKkBaeqbmxv1yX5LLA/sD7J4qpan2QxsK7tvhZY2vP0JW3b9Nc8DTgNYL/99qu5jF+SJJjsS516GogkSVpQkjwyybZTf9NMmvld4BzgyLbbkTSTbNK2L0+ydZIlwF7AN0cbtSRJC4sjKyRJ0kKzC3B2kgK2BVYBnwe+BqxK8vvATcDLAapqTZKzgMuB+4HXVNVdY4lckqQ50MXTRixWSJKkBaWq/o3mMqTT/Rg4eMBzTgZOnsu4JEmaa/PptBFPA5EkSZIkSZ1isUKSJEmSJHWKxQpJkiRJktQpFiskSZIkSVKnWKyQJEmSJEmd4tVAJHXGMWdc0rf99KP3H3EkkiRJksbJkRWSJEmSJKlTHFkhqfMGjbiQJEmSNJksVkj9nLl83BFIkiRJ0oLlaSCSJEmSJKlTLFZIkiRJkqRO8TQQSZohr1oiSZIkzS2LFdJCMWgejiNWjTYOSZIkSdoMTwORJEmSJEmdYrFCkiRJkiR1isUKSZIkSZLUKRYrJEmSJElSp1iskCRJkiRJnWKxQpIkSZIkdYrFCkmSJEmS1CkWKyRJkiRJUqdYrJAkSZIkSZ2yaNwBSFqYjjnjknGHIEmSJKmjHFkhSZIkSZI6xZEVWtjOXD7uCCRJkiRJ0ziyQpIkSZIkdYrFCkmSJEmS1CkWKyRJkiRJUqdYrJAkSZIkSZ1isUKSJEmSJHWKxQpJkiRJktQpFiskSZIkSVKnLBp3AJLG7Mzl/duPWDXaOCbAMWdc0rf99KP3H3EkkiRJ0vzmyApJkiRJktQpFiskSZIkSVKnWKyQJEmSJEmd4pwVkmbFoPkaJEmSJGmmHFkhSZIkSZI6xWKFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsWrgWhhOHP5uCOYfwbtsyNWjTYOSZIkSQuOxQpJmmODLut6+tH7jzgSSZIkaX7wNBBJkiRJktQpFiskSZIkSVKneBqIpJkZOP/Hm0YahiRJkqTJ5cgKSZIkSZLUKY6skKSOcUJOSZIkLXQWKzRZvETpnLvs+tv6th/P2/u2f3CXk+YynHltUFFiLl/bgockSZLmA4sVkvoaVJSQJEmSpLnmnBWSJEmSJKlTLFZIkiRJkqRO8TQQacIMOn1jn6WPHXEkjeNvGn4ui5n0lSRJkjS5LFZofnIizRlzDorJNZcTdUqSJEnjYLFCmqcsPkiSJEmaVBYrJM1bnjYiSZIkTaZ5UaxIcghwCrAV8DdV9adjDkmj4ukekqQOMBeRJGm0Ol+sSLIN8CHgWcCNwNeTfKWqLh1vZNLs6trEmHNt0KgIDTYbc1MMeo3Tj95/i197Nl9/ruOUZsJcRJKk0et8sQJ4GnBVVV0PkGQVcChggjAfDRopccSqWXn5Lv2Hf7bmlHBuipmbL6eHzLRg06X4Z1o46VrxoWvxqPPMRSRJGrH5UKxYAlzfc38tcOB4QtFG5rj4MC4zKRBM6siHLnEUxggM+Ld8zN1vmtPVXvbeF/R/YIaFmbksPsz0tS2ETCRzEUmSRixVNe4YNinJEcCzq+o17f0VwIFV9eqePscBx7V3fx343iyHsRNw8yy/pgZzf4+e+3y03N+j5f5+0O5VtXjcQcw3w+Qibbv5yMy4Pd02SdszSdsCbk/XuT2bN1Q+Mh9GVqwFlvbcX9K2PaCqTgNOm6sAkqypqv3m6vW1Iff36LnPR8v9PVrub82CzeYiYD4yU25Pt03S9kzStoDb03Vuz+x52DhWOkPfBPZKsiTJ1sBy4NwxxyRJkhYOcxFJkkas8yMrqurnSV4LfJmmuPLJqloz5rAkSdICYS4iSdLodb5YAVBV5wDnjDGEORvSqb7c36PnPh8t9/doub+1xTqQi8DkfZbdnm6bpO2ZpG0Bt6fr3J5Z0vkJNiVJkiRJ0sIyH+askCRJkiRJC4jFilaSjyZZl+TKAY+/Ocll7XJlkvuS7DDqOCfFEPt71yQXJPluku8nec2oY5w0Q+zzHZOc2+7zbybZa9QxTpIkS5Nc3H5ffD/JCX36JMkH2n3+7SRPGUesk2DI/b1nkq8nuSvJm8YRp7Q5m/uubvscmOSSJN9JcvEo45upScuvJi1/maTcYNKOu5N2XBtye45KckXb51tJOntFjSG358Xt9lze9nvhOGIdxjDb09N3/yT3JnnpnMflaSCNJM8Gbgc+XlWb/CJOsgx4Y1U9byTBTaDN7e8kJwFbV9UJSRYDPwAeV1V3jjjUiTHEPv8gcHNVvTvJnsDHquqAUcc5KZLsCuxcVZcneRRwKfCyqrqsp89hwCuB3wWeTLPPf3MsAc9zQ+7vnYHdafb3rVV1yniilQYb4rt6V+AC4KCqujHJTlV186jjHNak5VeTlr9MUm4wacfdSTuuDbk9TwOuqaqftP+xf09V7TOmkDdpyO3ZDrijqirJ3sAXq2q3MYW8ScNsT9tvK+B84OfAR6vqs3MZlyMrWlV1MXDLkN1XACvnMJyJN8T+Xgs8KkmA7YCbgbtGEdukGmKf7wlc2Pa9Btg5yeNHEdskqqobq+ry9u//Ai4Hpu/PQ2muKlBVdSmwKMnSEYc6EYbZ31W1rqouAe4ZQ4jSUIb4rj4c+HRV3dj272yhAiYvv5q0/GWScoNJO+5O2nFtyO35l6r6SXv3H6c/3iVDbs/t9eDIgEcCN442yuEN+e8H4Hjg74B1o4jLYsUMJdkWOITmTdLc+QjwROA/gSuAP6qq+8cb0sS7Avg9gCRPpanUd7L6O98k2QPYn+bA22sJcH3P/bVtm7bAJva3NAn2BB6X5Bvt8OI/GHdAs2GC8qtJy1/mZW4wacfdSTuuDbk9rwa+MIp4ttSmtifJS5JcA5wHvH60kT00g7anLVS+BPirUcVisWLmlgH/VFXD/kqgh+ZEmorefwP2AU5N8ujxhjTx3g3skuS7wAnAGsDzxLZQOwTws8Aben4t0Bxxf2sBeBjNcfEg4LnACV2eR2AGJiW/mrT8Zd7lBpN2HFiI25PkQOAY4C0jDO0h2dz2VNVZVbUnzXfcx5N0+v/fm9mevwBOGGUBdtGoVjRBDqfjQxQnxLOAk9qhUz9M8iOaXyq+Md6wJlf7hXTE1P0k/wZ8f3wRzX9Jtqb5lXBlVX2uT5e1wFIe/Fwvadv0EAyxv6VJcD1wQ1XdAdyR5B+AvYGBE3LOE5OSX01U/jLfcoNJO+5O2nFtmO1p53Y4HXhhVf14lPHN1Ezen6q6OMkiYBfghlHEN1NDbM9+wKeas9zYCXhRknur6uy5iqnTlZ2uSfIY4DnA58cdywLwrzS/GpFkF5oD/bXjDGjSJXlM+yVKkiOBb0/AL1xj056vfDpwdVX92YBu5wCvaPs/Bbi/qq4f0FebMOT+libBl4BnJlnUnjpxAHDNmGPaIhOWX01U/jKfcoNJO+5O2nFtmO1JshvwOeCoqupsUQyG3p5f7Pn7KcA2jGiuh5kaZnuq6herao+q2oNm9MUfzmWhArwayAOSrAQOpKkS3QS8E9gaoKo+1PY5Gjikqg4fT5STY3P7u52R9pM0E7tsBbyvqj4ynmgnwxD7/OnAGTSz+/4QOKaqbh1LsBMgyTOBr9Gc7zs1XO5ttOf6tvs8wKk0Q7nvBo6tqjVjCHfeG3J/70ozhPnRbZ/bgSdW1U9HH7HU35D5yJuBV7Xtp1fVn44l2CFMWn41afnLJOUGk3bcnbTj2pDb89fAYcB17eP3VlUnL1865PacSFsco/k39Maq+tqoYx3GMNszrf8ZNFc3mdOrgViskCRJkiRJneJpIJIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjrFYoUkSZIkSeoUixWSJEmSJKlTLFaok5K8K0n1LD9Jcl6S3xx3bFsiycPbbdtnxOvduV3vHrP4mqckuXYzfXrfx/uT3JrkkiQnt5fbmlNJrm3X/fY+jz2zJ7Y95jqWYSQ5McmX+7Q/J8nnk6xLck97+6UkhycZ+ns8yeokV2zi8VOT3JZkmyTLkvxrkoc/1O2RpPnOfGTW12s+svFj5iMbP24+IsBihbrtJ8AB7fIqYClwQZLtxxrVlnk4zTXMR5ocADu3691jxOuFB9/HpwOHA58DjgKuSLLvCNZ/e7ve6Va0j3VCkscCJwAnT2t/A3ARcB9wPHAQ8Drgp8Df0lwrflgrgb2SPLHP+rcCXgp8rqruqqrVwM+AY2e+NZI0UcxHZo/5yMbMRzZcj/mIHmCxQl12b1V9o10+BxwJ7Aj8zpjjGokkvzDuGGZJ7/v45ap6D7A3cAPwqfagNJe+CDwxyV5TDT0Hwi/M8bpn4ljghqq6eKohyVOAU4A/qarfq6pVVXVxVX26qlYAzwRunsE6Pk9zwF/R57HnArvQJBBTPgK8MUlmuC2SNEnMRyaD+chwzEfUGRYrNJ9MDRd7fG9jkh2SnJbkpiQ/T/LPSZ42rc9W7ZC27ye5ux229rfT+rwuyQ+S3JXkh0neOO3xdyW5OcmTk3yjXdc1SZ4/rd/Lk1zerueOJJcmmao2/1d7+7HeIX/tUklekeTjSW4DVrevV0le1y+WaW27J1nZxnh3ku8mOaodUji17y6aWu8M999jk5yZ5PYkNyT5443fnuFV1W3AW4BfAX5rmOck2SfJRe0+vbt9L/9oiKf+B/CPbPhrxvOA7eiTHCR5a5JvJ7mz3d7zk+w9rc/zkvxLu7/uTHJlkuU9j2/qMzDIUcDZ09qOB9YBJ/V7QlV9vaq+My22Y5Nc1X6Or0vylp7+d9B8rpZPfy2a/bMOuLCn7fM079F/30zskrSQmI9Mi2Vam/lIf+YjD/Y3H9FQLFZoPlna3q6fakiyDfBV4NnA64EXAWuBryb5pZ7nfhh4B/A3wMHAq4HeA+QbgA8Aq4Dnt/1OSfLWaTFsC/x12/fQNpbPJHlM+zpPoKkEf4nmoPc7wKeAR7fPf157exIPDim9oef130fzBf07DDgg9JNkZ+DrwJNohuQdBJxKM9zyBuAVbdf/0bPemey/T7Wv+VrglW3/fkMZZ+LvgXsZ4sCT5peHL9IM4fy9Npb3AVsPua6VbBjvCpqD5B19+u4AvJ/mc/AymqGZX+15j3egSSquoPkMvJDm8/XI9vHNfQb6bd8uNL/ufGPaQ88GLqyqe4fZyCRvBv5fu77fAv4c+JPeBKGN7VfTM+Q1ydY0+/XTVXXfVHtVXQfcCLxgmPVL0gJhPjKA+chmmY88yHxEm1dVLi6dW4B30QwnW9QuS2kODncAu/T0Owa4C9i9p20r4BrgA+39PWkSgT8YsK6t2nV9aFr7+2kORo/oiamAZ/T0+Y227cXt/VcA6zexXdu1/Y+e1r5H2/6pPs8p4HX99k/P/fcAtwGLB6x3r/Z1DpzWPsz+e3L73N/t6bMtcBNw7TDv4yYevwH4qyE+D49vY/iNGX6OrqUZtrgYuAfYn+Y83VuB3wV+u33dPTbx2XhE2/+Vbdsz2uc8csBzNvkZGPCcg9vX/NVp7XcC75nWlp5/F4uAh7Xtj6ZJZE6c1v/twI+BRe39qe1/X0+fqf3w9D6xnQ+cNZPtcXFxcZmUBfOR3sfMR8xHetvMR1zmfHFkhbpsR5ov9HuAfweeAxxaVTf19DmYpvr7H0kWJVlE8+X597TVeppz3+6jmfynnz3bdX12Wvunab5wn9TT9rOq+qee+9e0t49rb78D7JjkY0l+K8l2w2xojy/NsP+U5wHnVNX6zfbc0DD77+k0vzicM/WkqvoZ8JWHGGuvYc89vIk2kUjysrbyP7R2v1xI82vGIe16z+0bUHJgkouT/Ixmu+8EHgv8Wtvle23bmWlmqJ4+wdpD+QxMbc8t/cKfdv8wHvx3cQ/wf9v2A2h+TfnM1HvZvp8X0vw68+sAVXU3zaRiL08eOPdzOXAdza9h0/24Jz5JWojMR4ZnPrIJ5iPmI5oZixXqsp/QVJ7/O80wyftohvz12olmaNo905ZX0xzwaW//qz2g9TP15T79wDp1f4eetjt7O9SDQ9QWtfevpBm+tifNweeWJJ+ZwcHs1iH7TbcjzQF0pobZf9sDt7cHlV4zmUhpI0kewZBxVzPs8AU0++cM4MY05+k+dQar/BTwcuAI4OyquqtPTL8MnEdzLu8Kms/e/m2Mj2hjubmNZTuahPLmJF9J8qvt41vyGZieLP0nsGRa2wVtTPuz4ZDdndrbH7DhezmVzO7Y03clsBtwQPs+vJjmV7TpiUi/mCRpoTEfGZ75yOaZjzzIfESbtGjcAUibcG9VrWn//pcktwN/m+TMqvpq234LzZffG/o8f+rL/8fAo5JsOyBBmDogL57WPnW/X3V5oKo6Gzg7yaNpzh/8IPAh4CUzeZ0e99EM/+v1yGn3H2q1eZj9dxuwXZKHT0sQdurznJl4Ls13UL/q+Uaq6grgxe35jM8A3gt8Mcnjq+qeIV7iLJr34WU053b2s4zmYPjSqroTHjg/dYPzO6vqa8BBaWZIfy7wlzS/fD25fXymn4Eb29sd2TDpuhh4fpKtphLRqroVWNPG1vt+TH1On0//JPN7PX9fRJPwHE7zK9yj2HDW7V479MQnSQuR+UjDfATzEfMRjZIjKzSfrASuork+95QLgCcAP6iqNdOWqRmnL6Q5uPa7PBI0Qydvpqk+93opzbWjr9joGUOoqp9W1Srg74Cp60hPfZkPOxETNNXqqSF/tEPlDprW5wLghUl2pL9B6x1m//0zzUH8RT0xbEtzEHpI0lzD+73AD2km1BpaVd1TVX9PM1nTYoZMUqqZ8fu9NO/HoHU+nCY5uL+n7cVA38u2VdWdVXUOzSRnG10rfMBnoJ+pff0b09o/SJP0vW0Tz53ydZpf2nbt816uqaqpmd+nfoH7NE2idARwdU2bxbvHXsDlQ6xfkhYK8xHMR8xHBjIf0axxZIXmjaqqJP+H5teMZ7XV5I8DrwH+IckpwI9ozuc7gGZSob+squ8lOQ34YDtL9T/SVGcPq6pXVtV9SU4G/izN5be+CjwL+CPg7VX182FjTPIHwL7Al2kSjl+jSTLOarfh7iQ/Al6a5Cqag/bmvni/AByd5Fvt9h3LxgfE99MMSb0oyUk0CcUTge2q6s9ozrG9EzgqyU+B+9pfiYbZf5cmOR/4cHu+403Am2iG9A1jUZKpGbYf1e6f19JMinVIz9DVgdJcquskmmGOP6JJCt5CU50fuspeVe/YTJcLgD+luZTbXwO/DLyV5tecqVgOpTmgfoFmWORS4DjgH9rHN/kZGBDXuiSX0ZyP+7me9kuTvAn48yT70MwO/580v2Q9A9iVZhIrquq2JO8CPpTm8nBfo0l0ngA8p6qmXx5sJc2lyF7Chgn3A5Ls3q5jNs4HlqSJYD5iPoL5iPmIRqM6MMuni8v0hQGzNtP8IvF94NyetsfQDHu7nuaAtZ5mpu5nTXve24B/a/vcBHxi2msfT1NZv7vt98YhY3pgdmyaL/fz2hjupfki/wCwbU//Q2kOave2z92DB2ff/u0+r/8Y4DPAz2gOhG8H3j09FmB3moPHre02XAW8oufxV9FMWnRf809/Rvtve5pzLO9o9907aGa1vnaI97Ha5X6ag+wa4GSaivuwn4ddaA5m17X77Vaag+0vbeZ51wKnbOLxjWbfpjnQ/wfwc5rJvp7W+zo0536e3b6397b7629oZz4f5jMwIJb/BfxowGMH0iQj69v3aB3NBGOHA5nW90jgWzTJ4M/a/X3CgNf9Ubv9vzLg8dfT/JvIpmJ3cXFxmdQF8xGmbZ/5iPmI+YjLyJa0HwBJ0hi1Q1GvpfmF7YIxhwNAksuB06rq1HHHIkmS5p75iLrEOSskqQPqwXNY/3jcsQAkWUYzw/hHxh2LJEkaDfMRdYkjKySNXZKHsYniaTWXCpMkSZoz5iNStziyQlIXfJSNr63+wNJO0CRJkjSXzEekDnFkhaSxaw/+m7rk1+W14TXVJUmSZpX5iNQtFiskSZIkSVKneBqIJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6pT/DzyvenPpi0AqAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1296x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"if l_flv[l_index]=='mu':\n",
" lower_MC = 1.8\n",
" upper_MC = 2.1\n",
" lower_data = 1.0\n",
" upper_data = 3.0\n",
" \n",
"if l_flv[l_index]=='e':\n",
" lower_MC = 1.7\n",
" upper_MC = 2.15\n",
" lower_data = 1.6\n",
" upper_data = 2.4\n",
"\n",
"#plt.suptitle('Reconstructed D_s mass', fontsize=15)\n",
"plt.subplot(1,2,1)\n",
"label=\"Ds_ConsD_M\"\n",
"plt.title('MC', fontsize=15)\n",
"left_h_Ds=[MC_Ds_tuple_dict_presel_8[label][i][0]/1000 for i in range(len(MC_Ds_tuple_dict_presel_8[\"Ds_ConsD_M\"]))]\n",
"label=\"Dplus_ConsD_M\"\n",
"left_h_Dplus=[MC_Dplus_tuple_dict_presel_8[label][i][0]/1000 for i in range(len(MC_Dplus_tuple_dict_presel_8[\"Dplus_ConsD_M\"]))]\n",
"plt.hist(left_h_Ds,alpha=0.7,bins=70,range=(lower_MC,upper_MC), label='MC after presel cuts',density=False);\n",
"plt.hist(left_h_Dplus,alpha=0.7,bins=70,range=(lower_MC,upper_MC), label='MC after presel cuts',density=False);\n",
"plt.legend(fontsize='12')\n",
"plt.ylabel('# events (a.u.)', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"\n",
"label=\"Ds_ConsD_M\"\n",
"plt.subplot(1,2,2)\n",
"plt.title('Data', fontsize=15)\n",
"right_h=[data_tuple_dict_presel_3[label][i][0]/1000 for i in range(len(data_tuple_dict_presel_3[\"Ds_ConsD_M\"]))]\n",
"plt.hist(right_h,alpha=0.7,bins=70,range=(lower_data,upper_data), color='green',label='Data after presel cuts',density=False);\n",
"plt.legend(fontsize='12')\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(18,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+l_flv[l_index]+'_dataMC_after_presel.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": 105,
"metadata": {},
"outputs": [],
"source": [
"Dplus_bf_phicut=np.float(MC_Dplus_tuple_dict_presel_8[\"Dplus_ConsD_M\"].shape[0])\n",
"Ds_bf_phicut=np.float(MC_Ds_tuple_dict_presel_8[\"Ds_ConsD_M\"].shape[0])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# PHI MASS CUT"
]
},
{
"cell_type": "code",
"execution_count": 106,
"metadata": {},
"outputs": [],
"source": [
"if l_flv[l_index]=='mu':\n",
" lower_phi_mass = 980\n",
" upper_phi_mass = 1060\n",
" \n",
"if l_flv[l_index]=='e':\n",
" lower_phi_mass = 850\n",
" upper_phi_mass = 1100\n",
" \n",
"#Cut on phi mass\n",
"\n",
"MC_Dplus_indices=[]\n",
"MC_Ds_indices=[]\n",
"data_indices=[]\n",
" \n",
"for i in range(len(MC_Ds_tuple_dict_presel_8[\"Ds_ConsD_M\"])):\n",
"\n",
" phi_m = MC_Ds_tuple_dict_presel_3[\"phi_M\"][i]\n",
" #fixing a window on the phi mass\n",
" if lower_phi_mass<phi_m<upper_phi_mass:\n",
" MC_Ds_indices.append(i)\n",
" \n",
"for i in range(len(MC_Dplus_tuple_dict_presel_8[\"Dplus_ConsD_M\"])):\n",
"\n",
" phi_m = MC_Dplus_tuple_dict_presel_3[\"phi_M\"][i]\n",
" #fixing a window on the phi mass\n",
" if lower_phi_mass<phi_m<upper_phi_mass:\n",
" MC_Dplus_indices.append(i)\n",
" \n",
"for i in range(len(data_tuple_dict_presel_3[\"Ds_ConsD_M\"])):\n",
"\n",
" phi_m = data_tuple_dict_presel_3[\"phi_M\"][i]\n",
" #fixing a window on the phi mass\n",
" if lower_phi_mass<phi_m<upper_phi_mass:\n",
" data_indices.append(i)\n",
"\n",
"MC_Ds_tuple_dict_presel_9 ={}\n",
"MC_Dplus_tuple_dict_presel_9 ={}\n",
"data_tuple_dict_presel_4={}\n",
"\n",
"\n",
"for label in return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0): \n",
" \n",
" data_tuple_dict_presel_4[label] = data_tuple_dict_presel_3[label][data_indices]\n",
"\n",
"for label in return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0): \n",
" MC_Ds_tuple_dict_presel_9[label] = MC_Ds_tuple_dict_presel_8[label][MC_Ds_indices]\n",
"\n",
"for label in return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0): \n",
" MC_Dplus_tuple_dict_presel_9[label] = MC_Dplus_tuple_dict_presel_8[label][MC_Dplus_indices]"
]
},
{
"cell_type": "code",
"execution_count": 107,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABCsAAAH4CAYAAABjWWd3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl0ldW9+P/3h0lKEghoEDBgqBXFodWKKCpO9RZKL221zoLlFsvCSm+t1wlvB9qLfi3VqqiVYqkzg7ZqhQpWpDjcixWsKAqORWtEKKAgiQiY7N8fOeQXNIGDZDiQ92uts8559t7Pfj7nwFpn53P23k+klJAkSZIkScoVLZo6AEmSJEmSpJpMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJG23iBgTESkiXquj/rVM/ZhPlH87IuZExJqI2BARr0bEryOiW6MELkmSdjo1xh0pIioj4v2ImB8RV0ZEl8/Q36URcXwDhCqpHpmskPRZfQT0jIg+NQsj4nCgJFNfs/xa4F7gH8BQ4KvAdcBXgJsbIV5JkrTzWgv0A44CzgTup2o8sSgiDtvOvi4Fjq/X6CTVu1ZNHYCknVY58HeqBgwLapSfCcwBqgcOETEYuAgYnlL6fY22j0fERKoSF5IkSXX5OKX0dI3jRyLiFuAJYGpE7J9Sqmii2CQ1AGdWSNoRU4HTIyIAMs+nZ8pr+hHw908kKgBIKVWklGY2eKSSJGmXklJaQ9UsiS8A/waQWV66OCI2RsTaiHggIvbefE5EvAnsDvysxtKS4zN1l0fEcxGxPiLKIuLRiPhio78xSYDJCkk75n5gT+CYzHF/oChTDkBEtKZqyuasRo9OkiTt6uYCHwNHZo7zgV8AJwLfoSoxMTMiNv/dczJVS0omUbWspB9VM0UBOlG1RPWrwGlAGTA7Ijo0+LuQ9CkuA5H0maWU1kTELKqWfjyZeZ6VUlqbmWwBVYOE3YB/Nk2UkiRpV5VS+igiVlH14wkppRGb6yKiJTAPWE7VDytPpJSei4iPgdJPLCshpXTpJ879K/Au8E3gzoZ+L5K25MwKSTtqKnBqROwGnMqnl4BslhovJEmS1IxU/0KSufPYgojYQNWMi+WZql7b7CTi+Ih4IiI+zJy7HijM5lxJ9c9khaQd9RBVUy6vBPKA6Z+oXw1sAHo0clySJGkXFxFtqZrFuSIi+gP3AYupWu5xBNA307TtNvrZh6olq+uAs6haVnI4sGJb50pqGC4DkbRDUkrlETGDqk0070splX+iflNE/C8wAPhxU8QoSZJ2WSdQ9TfNPODbwNvAd1JKCSAiumbZz2CqZmicmlJanzm3JdC+3iOWlBVnVkiqD7dQNaNiQh311wN9IuI7n6yIiBYRMbAhg5MkSbueiCgEfgm8DswG2gBsTlRknFXLqRuB1p8oa0NVsqKyRtk3gc/VV7ySto8zKyTtsJTSXKp2466rfnpE/BqYFBFHA3+iaoft/YGRwJt4txBJklS3VhGx+Y4fBcBhwPlAO2BgSqkiImYD50fE9cCDQB/gu7X09TIwKCIeAT4EXgEeA64GbouI3wH7AJcDaxrwPUnaCmdWSGoUKaX/As4A9gUmA48C/0XV4OD8JgxNkiTlvg5ULfX4P6r2pTgVuBs4OKX0LEBK6X5gDHAO8GdgIFWzIz7p8szzbGA+cFimj5HAccDDwHDgbKpucyqpCcSWs6QkSZIkSZKaljMrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOWUVk0dQH3bY489UklJSVOHIUlSznn22WdXpZSKmjqO5sDxiCRJtct2PLLLJStKSkpYsGBBU4chSVLOiYi3mjqG5sLxiCRJtct2POIyEEmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmn7HIbbEpSY6usrGTVqlWsWbOGioqKpg5HzVjLli0pLCxkjz32oEULf4+QpF3Bpk2bKC0t5aOPPmrqUKSs1Nd4xGSFJO2g0tJSIoKSkhJat25NRDR1SGqGUkps2rSJFStWUFpaSo8ePZo6JElSPSgtLaWgoICSkhLHGMp59Tke8WcXSdpB5eXl7LXXXrRp08ZBhJpMRNCmTRv22msvysvLmzocSVI9+eijj9h9990dY2inUJ/jEZMVklQPnHKvXOH/RUna9Zio0M6mPsYjjmgkSZIkSVJOMVkhSWoQL774Ir179yYvL4/x48c3dTgNqqSkhNmzZzd1GJIkqYaKigrOPvts8vPz6du3b1OH06DGjBnDkCFDmjqMeuUGm5LUAIbfPr9B+5807PCs25aUlLBs2TKWLVvGHnvsUV1+6KGHsnDhQpYuXUpJSQkAzzzzDGPGjGHevHls2rSJfffdlxEjRnD++edvd4zjxo1j4MCBXHfddQAMGzaM4uJixo4du919NUdjxozh9ddf5+67727qUCRJOWTwlMEN2v/0s6Zn3bakpIQVK1bQqlUrWrduzf7778+5557LiBEjsloG8Oabb9KzZ082bdpEq1b1/6fp448/zhNPPMGKFSvIy8tj7ty5DBkyhNLS0nq/1q6oof99tsWZFZLUDPTs2ZMpU6ZUHy9atIgPP/xwizbz5s3jxBNPZNCgQZSWlrJu3Tp+97vfMXfu3M90zdLSUg488MAdCXsLH3/8cZOeL0mSPm369OmsW7eO5cuX84tf/IJf//rXDB8+vKnDAqrGIiUlJeTl5dVLf/UxlvA299kzWSFJzcDQoUO58847q4/vuOMOzj333C3aXHLJJXz/+99n1KhR5OXlEREcdthhTJs2rdY+33jjDfr370/Hjh1p3749p5xyCu+99x4AJ554Io8//jijRo0iPz+fiRMncs899zBu3Djy8/MZPLjqV6E333yTQYMG0aFDB7p06cLVV19d3f+YMWM49dRTGTJkCIWFhdx+++2fimHYsGGMHDmSr371qxQUFNC3b19ee+216vqI4Oabb2a//fajV69eACxcuJD+/ftTUFBAjx49uOOOO6rbP/jgg3zhC18gLy+Prl278stf/rK6btq0aey3334UFBRw6KGH8swzz2T12ZeXlzNy5EiKioooKCigX79+rF+/nrlz51JcXLxF283LSWbNmsVVV13FtGnTyM/P50tf+hIAEyZMoHv37uTl5dGjRw/uuuuurGKQJKmhtWnThpNOOokHH3yQO+64gxdffBGoSmYcfPDBFBQU0LlzZy677DJSSgAce+yxABQWFpKfn8+8efO2Or6ozQUXXEC3bt3Iy8vjoIMOql6WOWnSJEaMGMG8efPIz8/nkksu4Wtf+xrLli0jPz+f/Px8li1bRmVlJT/5yU/o1q0b7du3Z/DgwaxcuRKoGqdEBJMmTaKkpISvfOUrn7r+5u/zq666is6dO7Pnnnvyu9/9rrp+2LBhnH/++QwaNIiCggL++te/sn79es4//3yKioooLCzk3HPPrf4R6d1332XAgAHk5+dTWFjIUUcdRWVlZXU8dY2btmXy5Mnsv//+5OfnU1JSwsMPPwx8eilrzeUktf37LFmyhKOOOoq8vDw6derEqaeemnUM28tkhSQ1A0ceeSQffPABS5YsoaKigqlTp26xrvHDDz9k3rx5nHLKKdvV79ixY1m9ejVvvfUWa9euZfTo0QDMmTOH/v37c9NNN1FWVsaIESM455xzuPTSSykrK2P69OlUVFTwta99jWOOOYbVq1czf/58Jk6cyAMPPFDd/4wZMzj77LNZs2ZNneswp0yZwtixY1mzZg3HH388Z5555hb1M2fOZP78+SxevJg1a9YwYMAARowYwdq1a5k5cyYXXXQRzz77LADDhw/ntttuo7y8nFdeeYWBAwcC8NRTTzFq1CimTp3KunXruPjii/nmN7/J+vXrt/kZXXDBBSxdupTFixfzwQcfcMMNN2xzauzAgQO54oorOOOMMygrK+P5559nzZo1XHLJJcyePZvy8nKeffZZ+vTps83rS5LUmA444AB69erFk08+CUCHDh247777WLduHY8//jh33XUXU6dOBeCJJ54AYM2aNZSVldGvXz+g7vFFbY499lhefvll1q1bx8iRIznttNP48MMPGT58OBMmTKBfv36UlZXxq1/9ipkzZ9KtWzfKysooKyujW7duXH311cyePZvnnnuO1atX0717d84777wtrvH000/z8ssv88gjj9Qaw/Lly1m3bh3vvvsuDz30EBdeeCELFy6srp82bRq/+MUv+OCDD+jfvz8//OEPWb58Oa+//jrLli3jgw8+qH6P11xzDT179uT9999n9erVXHvttUREVuOmuvz1r3/lggsuYOLEiZSVlfG3v/2NffbZZ5vn1fbv8+Mf/5ivf/3r1bNpLr744m3281mZrJCkZmLz7IpHH32U3r17s9dee1XXvf/++1RWVlJUVJR1f/vssw/HHXccLVq0oGPHjvzoRz+q/lLLxlNPPUV5eTlXXHEFrVq1onv37nzve9/j3nvvrW5z9NFHM2jQIADatm1baz/f+MY36Nu3Ly1btmTMmDEsWrRoi9kVl112Ge3bt6dt27b86U9/olevXgwdOpQWLVpw4IEHcuqpp/KHP/wBgLy8PJYsWcK6deto37599YyGSZMmMXLkSA499FAAzjnnHNq3b7/N9/vRRx8xefJkrr/+eoqKiogI+vbty2677Zb157RZ69atadmyJYsXL2b9+vUUFRXRu3fv7e5HkqSGtscee1TPhjj22GPZf//9Aejduzdnn332Vr8/t3d8ccYZZ9C+fXtatGjBqFGjaNmyJYsWLco61ltvvZWxY8ey55570rp1a37yk58wY8aMLZbL/vSnP6Vt27Z1jkVatmzJT3/6U1q2bMkRRxzBt771Le67777q+pNPPpk+ffoQEUQEd955J9dccw0dOnSgXbt2XHbZZdXjn7y8PJYvX84///lPWrZsSb9+/YiIrMZNdbntttsYMWJE9UyJPffck/322y/rz6imvLw8/vnPf/Luu+/Spk0bjjzyyM/UTzZMVkhSMzF06FAmT57M7bff/qklIB07dqRFixasWrUq6/5KS0s55ZRT6Ny5Mx06dODMM8+krKxsu85ftmwZhYWF1Y+rrrqK999/v7pN165dt9lPzaUU7dq1o1OnTqxYsaLWPkpLS/nb3/62xTXvueee6gHVvffey0MPPUSPHj04+uijqwdHpaWlXHvttVuc9/bbb2/z81q9ejWbNm3i85//fHYfylbk5eUxefJkbrzxRrp27cqAAQOqp9jq0yKie0Q8EREvRsSrEXFZpnxMRLwTEQszj0E1zhkdEUsy5wyoUX5YRDwXEYsjYnxERFO8J0naWaxatYpOnToB8OSTT3L00UfTsWNHCgsLq2dd1mV7xxf/8z//wxe+8AXat29PYWEh77333naPR04++eTq7/fevXvTunVrVq9eXd1mW+ORTp068bnPfa76uLi4uM6xyMqVK9mwYQOHHXZY9TUHDhzI2rVrAbj44ovp0aMHJ510EnvvvTdjx44lpZTVuKku7777br2MRQCuvvpqNm7cyOGHH87+++/Pb3/723rptzYmKySpmdh7773p2bMnDz/88KeWe7Rr145+/fpx//33Z93f5ZdfTn5+Pq+99hpr165l6tSp1WtQa/PJv++6dOlCr169WLNmTfVj3bp1zJo1a7ve1zvvvFP9ev369bz33nvsueeetbbt0qULJ5100hbXLCsrq/6iPfLII5kxYwYrV67ktNNO4/TTT68+b8yYMVuc9+GHH3LOOedsNbbdd9+dNm3asHTp0k/VtWnTZotfbSorK7cYcNT29/CgQYOYM2cOy5cv5+CDD/7UNFVtYRMwKqV0EHAYcF5EHJKpuy6ldEjm8TBUJSSAbwNfBAYCv42IzVNgbgPOSykdAOwNnNyYb0SSdiZLlizh1Vdf5ZhjjgHgrLPO4uyzz2bFihWsWbOGUaNGVY8Xavuu257xxezZs7n55puZMWMGa9euZc2aNey+++51tq/tel26dGH27NlbfMd/9NFHdO/ePev3/N57722xNLS0tLTOscjuu+9O69atefXVV6uvt3bt2uoxQfv27Rk/fjxLly5l5syZjB8/nkceeWSHxk3dunWrdSwCnx6P1EzS1PZ5devWjdtuu41ly5bx+9//nv/8z//k5Zdf3mYMn4XJCklqRiZNmsScOXNq3RV73Lhx/OY3v+GWW26p/tJ6/vnnP7UHxGbl5eXstttu5Ofns2LFCq655pqtXrtTp0689dZb1cfHHXcclZWV3HjjjWzcuJGUEq+88kr1/hHZeuihh1iwYAEVFRX8/Oc/58ADD2Tfffette3JJ5/MwoULue+++6ioqKCyspLnnnuOl19+mY0bN3LvvfdSXl5Oq1atKCgoqP6SPu+887jlllv4+9//DlQt7/jLX/7CunXrthpb27ZtOeuss7joootYtWoVKSXmz5/Phg0b6N27N2VlZfz5z3+msrKScePGUV5evsXn9fbbb1dvqrVixQpmzpzJhg0baNOmTfUmqKpdSml5SumFzOt1wAvAXls55evAtJTSppRSKfAS0DciegAtU0qb/2PenWkrSaph06ZNzJkzh5NPPpkhQ4Zw8MEHA1Xjhfz8fNq0acNzzz3H5MmTq88pLCwkInjzzTery7ZnfFFeXk6LFi0oLCykoqKCcePGbXUzzk6dOvH+++/zwQcfVJeNGDGCH//4xyxbtgyoWhq7efPJbFVUVDB27FgqKip45pln+NOf/lTnxpNt27Zl6NChXHzxxdU/Uixfvrx6k8tZs2ZVfx4FBQW0bNmSiNihcdOwYcOYOHEi//u//wtUjSleffVVAL70pS8xdepUPv74Y1544YUtlq/U9u/z4IMPsnz5cqBqP5IWLVo02HjEZIWUKyafUftDqkf77LNPnZsyHnXUUTz22GM89NBDdOvWjfz8fIYNG8YJJ5xQa/sxY8bw9NNPU1BQwMCBA6vv8FGX4cOHs2DBAgoKCvjWt75Fq1ateOSRR3jssceqp3oOGTJkq4OM2px55pmMHj2awsJCZs+eXb1pV206derErFmzmDBhAh07dqRTp0788Ic/rP415NZbb2WvvfaiXbt2jB8/nnvuuQeoSqz86le/4jvf+Q75+fn06NGDW265Jav4br75ZoqLi9lvv/3o0KEDF154IZWVlXTs2JEbbriBoUOH0rVrV1q1arXFkpbTTjuN9evX06FDB7785S9XD4Q231Xk0UcfZcKECdv1WTVXEVECHA48lSm6ICJejoh7ImL3TFkx8HaN00ozZXWVS2oEg6cMrvWh3DF48ODqO33893//Nz/4wQ+47bbbqutvuukmRo8eTfv27fnJT36yxezODh06cNFFF9GnTx8KCwt5+umnt2t8MWjQIE488UQ+//nPs/feewNsdUbEQQcdxDe+8Q2Ki4spLCxk2bJl/Pd//zfHHHMMRxxxRPUdvx5//PHt+gy6dOlCu3bt6NatG//+7//OtddeW73PVW1uuukmOnbsSO/evSkoKODYY4+t3mfjpZdeon///rRr144+ffrw3e9+l69+9as7NG464YQTGD9+PMOGDSM/P58jjzySN954A4Arr7ySl156iQ4dOjB69OjqWaVQ+7/Pk08+ySGHHEK7du342te+xi9/+cvPvP/FtsTWpuzujPr06ZMWLFjQ1GFI26+uxMTZtd82UrljyZIlbnTYRIYNG0ZxcTFjx45t6lBySl3/JyPi2ZRSs7qFSETkA3OBq1JK90dEEfAekIAxwD4ppXMiYiIwJ6U0NXPebzPnvQX8NKU0MFPeDxiTUhpQy7VGACMAevTocVjNmUSSPpu6EhPTz5reyJE0HccZuW3u3LkMGTKE0tLSpg4l5+zoeMSZFZIkaZcUEa2BPwJTUkr3A6SUVqaUKlJKlcAEqmZcQNWMiZo/xxVnyuoq/5SU0sSUUp+UUp/tubOOJEn6NJMVkiRpl5O5Y8ckYElK6doa5Z1rNPs2sDjz+mHgjIhoHRHFwEHAMymlfwKVEfHlTLtzgJkN/gYkSWrmWjV1AJIkfVa33357U4eg3HU0MBRYFBELM2VXAGdHxBeBNsA/geEAKaUFEfEAVRtxVgIjU0obMuf9B/D7iGgDzKFqtoYkSRx//PEuAWkgJiskSdIuJ6X0FFDb9uR1bvGeUroSuLKW8gXAIZ8+Q5IkNRSXgUiSJElSDtvVboqgXd/mW6/vCJMVkiRJkpSj2rZty+rVq01YaKeQUmLjxo2888475OXl7VBfLgORJEmSpBxVXFxMaWkpK1eubOpQpKy0atWKDh06sMcee+xYP/UUjyRJkiSpnrVu3ZqePXs2dRhSo3MZiCSpQbz44ov07t2bvLw8xo8f39ThNKiSkhJmz57d1GFIkiTtMpxZIUkNYfIZDdv/2dOyblpSUsKyZctYtmzZFtPxDj30UBYuXMjSpUspKSkB4JlnnmHMmDHMmzePTZs2se+++zJixAjOP//87Q5x3LhxDBw4kOuuuw6AYcOGUVxczNixY7e7r+ZozJgxvP7669x9991NHYokSVKjc2aFJDUDPXv2ZMqUKdXHixYt4sMPP9yizbx58zjxxBMZNGgQpaWlrFu3jt/97nfMnTv3M12ztLSUAw88cEfC3sLHH3/cpOdLkiSp8ZiskKRmYOjQodx5553Vx3fccQfnnnvuFm0uueQSvv/97zNq1Cjy8vKICA477DCmTat9Fscbb7xB//796dixI+3bt+eUU07hvffeA+DEE0/k8ccfZ9SoUeTn5zNx4kTuuecexo0bR35+PoMHDwbgzTffZNCgQXTo0IEuXbpw9dVXV/c/ZswYTj31VIYMGUJhYSG33377p2IYNmwYI0eO5Ktf/SoFBQX07duX1157rbo+Irj55pvZb7/96NWrFwALFy6kf//+FBQU0KNHD+64447q9g8++CBf+MIXyMvLo2vXrvzyl7+srps2bRr77bcfBQUFHHrooTzzzDNZffbl5eWMHDmSoqIiCgoK6NevH+vXr2fu3LkUFxdv0XbzcpJZs2Zx1VVXMW3aNPLz8/nSl74EwIQJE+jevTt5eXn06NGDu+66K6sYJEmSdjYmKySpGTjyyCP54IMPWLJkCRUVFUydOpUhQ4ZU13/44YfMmzePU045Zbv6HTt2LKtXr+att95i7dq1jB49GoA5c+bQv39/brrpJsrKyhgxYgTnnHMOl156KWVlZUyfPp2Kigq+9rWvccwxx7B69Wrmz5/PxIkTeeCBB6r7nzFjBmeffTZr1qzZIt6apkyZwtixY1mzZg3HH388Z5555hb1M2fOZP78+SxevJg1a9YwYMAARowYwdq1a5k5cyYXXXQRzz77LADDhw/ntttuo7y8nFdeeYWBAwcC8NRTTzFq1CimTp3KunXruPjii/nmN7/J+vXrt/kZXXDBBSxdupTFixfzwQcfcMMNN9Cixda/fgcOHMgVV1zBGWecQVlZGc8//zxr1qzhkksuYfbs2ZSXl/Pss8/Sp0+fbV5fkiRpZ2SyQpKaic2zKx599FF69+7NXnvtVV33/vvvU1lZSVFRUdb97bPPPhx33HG0aNGCjh078qMf/Ygnnngi6/OfeuopysvLueKKK2jVqhXdu3fne9/7Hvfee291m6OPPppBgwYBVfeZr803vvEN+vbtS8uWLRkzZgyLFi3aYnbFZZddRvv27Wnbti1/+tOf6NWrF0OHDqVFixYceOCBnHrqqfzhD38AIC8vjyVLlrBu3Trat29fPaNh0qRJjBw5kkMPPRSAc845h/bt22/z/X700UdMnjyZ66+/nqKiIiKCvn37sttuu2X9OW3WunVrWrZsyeLFi1m/fj1FRUX07t17u/uRJEnaGZiskKRmYujQoUyePJnbb7/9U0tAOnbsSIsWLVi1alXW/ZWWlnLKKafQuXNnOnTowJlnnklZWdl2nb9s2TIKCwurH1dddRXvv/9+dZuuXbtus5+aSynatWtHp06dWLFiRa19lJaW8re//W2La95zzz3Vy1fuvfdeHnroIXr06MHRRx9dnYwoLS3l2muv3eK8t99+e5uf1+rVq9m0aROf//zns/tQtiIvL4/Jkydz44030rVrVwYMGMCLL764w/1KkiTlIu8GIknNxN57703Pnj15+OGHmTRp0hZ17dq1o1+/ftx///0cccQRWfV3+eWXk5+fz2uvvUaHDh2YMWMGI0eOrLN9RGxx3KVLF3r16sXixYu3/83U8M4771S/Xr9+Pe+99x577rlnrW27dOnCSSedxMMPP1xr/ZFHHsmMGTP4+OOPuemmmzj99NNZvnw5Xbp0YcyYMVx66aXbFdvuu+9OmzZtWLp0Kfvvv/8WdW3atNlik9PKysotEjWf/LwABg0axKBBg/joo4/48Y9/zHnnncfTTz+9XTFJUi4aPGVwU4cgKcc4s0KSmpFJkyYxZ84c8vLyPlU3btw4fvOb33DLLbdU/xH9/PPPf2oPiM3Ky8vZbbfdyM/PZ8WKFVxzzTVbvXanTp146623qo+PO+44KisrufHGG9m4cSMpJV555ZXq/SOy9dBDD7FgwQIqKir4+c9/zoEHHsi+++5ba9uTTz6ZhQsXct9991FRUUFlZSXPPfccL7/8Mhs3buTee++lvLycVq1aUVBQUJ0wOO+887jlllv4+9//DlQt7/jLX/7CunXrthpb27ZtOeuss7joootYtWoVKSXmz5/Phg0b6N27N2VlZfz5z3+msrKScePGUV5evsXn9fbbb1NZWQnAihUrmDlzJhs2bKBNmzbVm6BKkiTtikxWSFIzss8++9S5KeNRRx3FY489xkMPPUS3bt3Iz89n2LBhnHDCCbW2HzNmDE8//TQFBQUMHDiw+g4fdRk+fDgLFiygoKCAb33rW7Rq1YpHHnmExx57rHopyZAhQ6qXZGTrzDPPZPTo0RQWFjJ79mymTp1aZ9tOnToxa9YsJkyYQMeOHenUqRM//OEPqzfKvPXWW9lrr71o164d48eP55577gGqEiu/+tWv+M53vkN+fj49evTgllug6w2+AAAgAElEQVRuySq+m2++meLiYvbbbz86dOjAhRdeSGVlJR07duSGG25g6NChdO3alVatWm2xpOW0005j/fr1dOjQgS9/+ctUVFQwduzY6ruKPProo0yYMGG7PitJkqSdRaSUmjqGetWnT5+0YMGCpg5D2n6Tz6i9/Ozabxup3LFkyRI3Omwiw4YNo7i4mLFjxzZ1KDmlrv+TEfFsSslbiDQCxyPS9tneZSDTz5reQJFIamjZjkecWSFJkiRJknKKyQpJkiRJkpRTvBuIJGmndfvttzd1CJIkSWoAzqyQJEmSJEk5xZkVUmOrayNN7dQqKytp0cL8r5re5ludSpIk7cwcWUvSDsrLy+Odd95h48aN7Gp3WNLOI6XExo0beeedd8jLy2vqcCRJknaIMyskaQcVFxezatUq3nrrLT7++OOmDkfNWKtWrejQoQN77LFHU4ciSZK0Q0xWSNIOatGiBZ07d6Zz585NHYokSZK0S2iSZSAR0TIinouIGZnjThHxaEQsioi/RETHGm1HR8SSiHgxIgY0RbySJEmSJKnxNNWeFT8EltQ4/jkwM6V0MDAzc0xEHAZ8G/giMBD4bUTs1sixSpIkSZKkRtToyYqIKAa+DvyuRvHXgbsyr+/OHG8un5ZS2pRSKgVeAvo2VqySJEmSJKnxNcXMiuuBS4Ga91YrSimtBMg8b174XQy8XaNdaaZMkiRJkiTtoho1WRER/w78K6X0bD33OyIiFkTEgpUrV9Zn15IkSZIkqZE19syKo4FvRMSbwFTgxIi4G1gZEUUAmed/ZdqXAt1rnF+cKdtCSmliSqlPSqlPUVFRQ8YvSZIkSZIaWKMmK1JKo1NKxSmlEuBMYE5KaQjwMDAk02wIVZtskik/IyJaZ/a6OAh4pjFjliRJkiRJjatVUweQ8TNgWkR8F1gBnA6QUloQEQ8AL1C1x8XIlNKGpgtTkiRJkiQ1tCZLVqSU5gJzM69XAyfV0e5K4MpGC0ySJEmSJDWpprgbiCRJkiRJUp1MVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJO1yIqJ7RDwRES9GxKsRcVmmvFNEPBoRiyLiLxHRscY5oyNiSeacATXKD4uI5yJicUSMj4hoivckSVJzYrJCkiTtijYBo1JKBwGHAedFxCHAz4GZKaWDgZmZYyLiMODbwBeBgcBvI2K3TF+3AeellA4A9gZObtR3IklSM2SyQpIk7XJSSstTSi9kXq8DXgD2Ar4O3JVpdnfmmMzztJTSppRSKfAS0DciegAtU0rP1nKOJElqICYrJEnSLi0iSoDDgaeAopTSSoDMc+dMs2Lg7RqnlWbK6iqXJEkNyGSFJEnaZUVEPvAH4MKU0toGvtaIiFgQEQtWrlzZkJeSJGmXZ7JCkiTtkiKiNfBHYEpK6f5M8cqIKMrUFwH/ypSXAt1rnF6cKaur/FNSShNTSn1SSn2Kiorq741IktQMmayQJEm7nMwdOyYBS1JK19aoehgYknk9hKpNNjeXnxERrSOiGDgIeCal9E+gMiK+nGl3To1zJElSA2nV1AFIkiQ1gKOBocCiiFiYKbsC+BkwLSK+C6wATgdIKS2IiAeo2oizEhiZUtqQOe8/gN9HRBtgDlWzNSRJUgMyWSFJknY5KaWngKij+qQ6zrkSuLKW8gXAIfUXnSRJ2haXgUiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKSYrJEmSJElSTjFZIUmSJEmScorJCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOWUVk0dgCRJkqTmYfCUwU0dgqSdhDMrJEmSJElSTjFZIUmSJEmScorLQCRJkqSdXF3LK6afNb2RI5Gk+uHMCkmSJEmSlFNMVkiSJEmSpJxiskKSJEmSJOUUkxWSJEmSJCmnmKyQJEmSJEk5xWSFJEmSJEnKKY2arIiIthGxICIWRsRrEXF9VOkUEY9GxKKI+EtEdKxxzuiIWBIRL0bEgMaMV5IkSZIkNb7GnlmxATgupXQIcADQDzgB+DkwM6V0MDAzc0xEHAZ8G/giMBD4bUTs1sgxS5IkSZKkRtSoyYpUpTxz2BpoCfwL+DpwV6b87swxmedpKaVNKaVS4CWgbyOGLEmSJEmSGlmj71kRES0jYiFVSYq5KaUXgaKU0kqAzHPnTPNi4O0ap5dmyj7Z54jM8pIFK1eubNg3IEmSJEmSGlSjJytSShWZZSDFQP+IOKEe+pyYUuqTUupTVFS040FKkiRJkqQm02R3A0kprQH+DBwJrIyIIoDM878yzUqB7jVOK86USZIkSZKkXVRj3w1kj4goyLz+HPBvwIvAw8CQTLMhVG2ySab8jIhoHRHFwEHAM40ZsyRJkiRJalytGvl63YA7IyKAtsCUlNL0iPg/YFpEfBdYAZwOkFJaEBEPAC8AlcDIlNKGRo5ZkiRJkiQ1okZNVqSUXgAOqaV8NXBSHedcCVzZwKFJkiRJkqQc0WR7VkiSJEmSJNXGZIUkSZIkScopJiskSZIkSVJOaewNNiVJkiRphwyeMrjW8ulnTW/kSCQ1FGdWSJIkSZKknGKyQpIkSZIk5RSTFZIkSZIkKaeYrJAkSZIkSTnFZIUkSZIkScopWd8NJCJaA/sDe2aKlgOvpJQ2NURgkiRJkiSpedpqsiIiWgGnA0OB44DdgMhUJ2BDRDwB3AXca+JCkiRJkiTtqDqXgUTEd4F/ANcB7wIXAv2B3sABmdcXAsuAa4E3MudIkiRJkiR9ZlubWXE+VcmIP6WUKupo83/AxIhoCXwLuBz4ff2GKEn1Y/jt82stnzTs8EaORJKkz2bwlMFNHYIkNYo6kxUppaxH75lkxh8zD0mSJEmSpM/Mu4FIkiRJkqScssPJiojoHhE96iMYSZIkSZKkrG9duhX/oCrp0bIe+pIkSZIkSc1cfSQrhvP/385UkiRJkiRph+xwsiKldGd9BCJJkiRJkgRusClJkiRJknJMVjMrIuLebbVJKZ2+4+FIkiRJ2tkNnjK4qUOQtJPLdhlIUS1lecDBwHvAK/UWkSRJkiRJataySlaklE6orTwi9gQeAK6rz6AkSZIkSVLztUN7VqSUVgD/DxhXP+FIkiRJkqTmrj422KwESuqhH0mSJEmSpKw32DygluKWQC/gSmBhfQYlSZIkSZKar2w32HwRSLWUB/ASMLzeIpIkSZIkSc1atsmK2jbY/Bh4N6X0j3qMR5IkSZIkNXPZ3g3k8YYORJIkSZIkCepng01JkqScExG/j4h/RcSLNcrGRMQ7EbEw8xhUo250RCyJiBcjYkCN8sMi4rmIWBwR4yMiGvu9SJLU3OxwsiIiXouIN+ojGEmSpHp0OzCwlvLrUkqHZB4PQ1VCAvg28MXMOb+NiN0y7W8DzkspHQDsDZzc4JFLktTM1cfMiicyD0mSpJyRUnoCeC/L5l8HpqWUNqWUSqnaQLxvRPQAWqaUns20uzvTVpIkNaAdTlaklIanlP6jPoKRJElqBBdExMsRcU9E7J4pKwbertGmNFNWV/mnRMSIiFgQEQtWrlzZEHFLktRsZHs3EEnba/IZTR2BJOnTbgb+h6pbso8BxgPn1EfHKaWJwESAPn361HbLd0mSlKXtSlZExDFAL6DtJ+tSSr+pr6AkSZIaQkqpespDREwA5mYOS4HuNZoWZ8rqKpd2CoOnDK61fPpZ0xs5EknaPlklKzLrNR8GDqDql4jNu2DX/NXAZIUkScppEdE5pfSvzOG3gcWZ1w8DEyLiemBP4CDgmZTShoiojIgvp5T+TtUsjLsbPXBJkpqZbGdWXA+sAroC7wKHA+8AZwPDgW80SHSSJEmfUURMAY4H9oiIUuBnwAkR8UWgDfBPqsYxpJQWRMQDwAtAJTAypbQh09V/AL+PiDbAHOCPjfpGJElqhrJNVpwInAtsnjoZKaXlwK8z9xr/DTCgrpMlSZIaW0rprFqKJ22l/ZXAlbWULwAOqcfQJEnSNmR7N5DPAatSSpXARmCPGnULgKPqOzBJkiRJktQ8ZZuseB3YK/N6EVDzNgffAj6sz6AkSZIkSVLzlW2yYhbwlczrq4FzI2JpRCwF/hO4qSGCkyRJkiRJzU9We1aklP6rxus/RsTxwClUbU41O6X0QMOEJ0mSJEmSmptsN9jcQkrpSeDJeo5FkiRJkiSp7mUgEdF9ezuLiL13LBxJkiRJktTcbW3PikURcWtEHL6tTiLi2Ii4k6rNNyVJkiRJkj6zrS0DORT4GfBURKwE/gYsBt4DAtgdOBA4AugITAG+2KDRSpIkSZKkXV6dyYqU0lJgWERcDpwF/BswDNgz02Q58ALw/4BpKaUVDRuqJEmSJElqDra5wWZKaTlwXeYhSZIkSQAMnjK4qUOQtIva2p4VkiRJkiRJje4z3bpUUiOafEbt5WdPa9w4JEmSJKmROLNCkiRJkiTlFJMVkiRJkiQpp5iskCRJkiRJOeUzJysioigijoiIvPoMSJIkSZIkNW9ZJSsi4sqI+EWN44HAW8D/AW9GxJcaKD5JkiRJktTMZDuz4mzg1RrHvwL+CvQDFgLj6jkuSZIkSZLUTGWbrNgL+AdARPQADgR+mlJ6BrgeOLJhwpMkSZIkSc1NtsmKD4EOmdfHAWtSSs9mjj8AdqvvwCRJkiRJUvPUKst2TwCXRwTAxcCfa9TtT9X+FZIkSZIkSTss22TFD4CpwJ+A54DRNerOBp6q57gkSZIkNZDBUwbXWj79rOmNHIkk1S6rZEVK6S2qNtOszTeB9fUWkSRJkiRJatayvXXpnIjYv47qLsAj9ReSJEmSJElqzrLdYPN4oH0ddQVUbbopSZIkSZK0w7LdswIgfbIgIloCxwLv1VtEkiRJkpqEe1lIyhV1Jisi4mfATzOHCXg6czeQ2txYz3FJkiRJkqRmamszKx4GVgEBjAeuBd78RJuNwMsppScbJDpJkiRJktTs1JmsSCnNB+YDRMQ64M8ppVWNFZgkSZIkSWqesr116R0NHYgkSZIkSRJkmayIiDbA5cA3ga61nZdS6ly/oUmSJEmSpOYo27uB/BYYAjwA/BmoaLCIJEmSJElSs5ZtsuLbwIUppZsbMhhJkiRJkqQWWbb7CFjUkIFIkiRJkiRB9smKW4GzGjIQSZIkSZIkyH4ZyArgnIiYBcwGPvxEfUop3VKvkUmSJEmSpGYp22TF9ZnnHsBXa6lPgMkKSZIkSZK0w7JKVqSUsl0uIkmSJGkXM3jK4KYOQVIzYxJCkiRJkiTllKyTFRFREhHXR8TfIuL1iDgwU/6DiDi64UKUJEmSJEnNSVbJiojoC7wADAD+DvQEdstUdwYua5DoJEmSJElSs5PtzIrrgJnAAcAPgKhRNw84op7jkiRJkiRJzVS2yYovA7eklBJVd/6oaQ1QWK9RSZIkSZKkZivbZMVaYI866noCy+snHEmSJEmS1NxldetS4CHgFxHxNPBupixFRAFwMXB/QwQnSZIkSdmq7Rar08+a3gSRSNpR2c6suAz4CHgNmJ0puwH4B1AB/LT+Q5MkSZIkSc1RVsmKlNL7wJHAKKpmVswGlgFXAEenlNY1WISSJEmSJKlZyXYZCCmljcCkzEOSJEmSJKlBZDWzIiKeiIjzI6JoRy4WEd0zfb0YEa9GxGWZ8k4R8WhELIqIv0RExxrnjI6IJZlzBuzI9SVJkiRJUu7Lds+KFcA1wDuZpMJ3ayYUtsMmYFRK6SDgMOC8iDgE+DkwM6V0MDAzc0xEHAZ8G/giMBD4bUTs9hmuK0mSJEmSdhLZ7llxGtAZ+A5QBtwMvBsRMyJiaOauINn0szyl9ELm9TrgBWAv4OvAXZlmd2eOyTxPSyltSimVAi8BfbN6Z5IkSZIkaaeU7cwKUkrlKaUpKaWTqUpcjMhU3Qos394LR0QJcDjwFFCUUlqZuc7KTP8AxcDbNU4rzZRJkiRJkqRdVNbJipoysyLeAJYCHwCf257zIyIf+ANwYUpp7WeJ4RP9jYiIBRGxYOXKlTvanSRJkiRJakLblayIiL4RcW1E/BN4AjgOuAHYdzv6aA38EZiSUro/U7xy8+admed/ZcpLge41Ti/OlG0hpTQxpdQnpdSnqGiH9gCVJEmSJElNLKtbl0bEL4HTgL2B14DbqNpLYvH2XCwigqpbny5JKV1bo+phYAhwXeZ5Zo3yCRFxPbAncBDwzPZcU5IkSdrZDJ4yuKlDkKQmlVWygqpExb3A1JTSwh243tHAUGBRRGzu5wrgZ8C0iPguVXceOR0gpbQgIh6gaiPOSmBkSmnDDlxfkiRJkiTluKySFSmlz9fHxVJKTwFRR/VJdZxzJXBlfVxfkiRJkiTlvqz3rIiI9hFxYURMi4jZEbFvpvz0iDig4UKUJEmSJEnNSbZ7VvQC5lB1149ngBOAgkz1EcApwJkNEaAkSZIkSWpesp1ZMR54haoNNgez5VKOJ4Bj6jkuSZIkSZLUTGW7wWZ/4JsppbKIaPmJupVA5/oNS5IkSZIkNVfZzqz4CGhbR11XYHX9hCNJkiRJkpq7bJMVjwJXRERejbKUmWUxCni43iOTJEmSJEnNUrbLQC4B/hd4A5gFJOBy4CAgDzirQaKTpEYw/Pb5tZZPGnZ4I0ciSZLq2+Apg2stn37W9EaORNL2yCpZkVJ6OyK+BFwEfIWqpMXewIPAr1NKLgORlDPqSj5IkiRJ2jlkO7OClNL7wE8yD0mSJEmSpAaR7Z4VkiRJkiRJjcJkhSRJkiRJyikmKyRJkiRJUk4xWSFJkiRJknKKyQpJkiRJkpRTsk5WRMS5EVHYkMFIkiRJkiRtz8yK24AeAFHlpxHRpWHCkiRJkiRJzVWruioiYiawEHg+8wggZapbAD8DZgDLGzhGSZIkSZLUjGxtZsUsoCswmqpkRQJuiogxwEC2TF5IkiTllIj4fUT8KyJerFHWKSIejYhFEfGXiOhYo250RCyJiBcjYkCN8sMi4rmIWBwR4yMiGvu9SJLU3NSZrEgp3ZBSGpZS+hJQQFVy4u/AfsB4qhIVd0XENRExsFGilSRJyt7tVP3AUtPPgZkppYOBmZlj/r/27j5arrq+9/j7I0Et4hMQsNck0GdqKUUFLT6iUMXS1FrUEEQXForaq1bXtVKtq9UuuLZX2t6q61axKFVLjFpBo6Ai2KKtViK1PCg+tEJJCyQIaEHk8Xv/2PuQyWQmZ05yzsw+c96vtfaaM7/5zezv3jOZ/c13fvu3kzweOBY4uH3Ou5M8qH3O+4CTq+oxwP7A8xY+dEmSlrahxYokr07y1CQPrao72+b3VdVamoJFgHXAnsA7Fz5USZKk0VXVJcDNfc3HAB9o//5ge3+mfX1V3V1Vm4CrgCckWQXsVlVfHfAcSZK0QIbOWQH8GvAHwD5JrqUZSXFckh8Drmj7XFBVly1wjJIkSfNleVVtAaiqLUn2bdtXABf39NvUtt0LXDegfTtJTgFOAVi1atU8h61ptnrd6kmHIEmds6PTQJ5VVfsBjwZ+h2YkxVE0c1ncTFO8eEWSI3uGSUqSJC1JVXVmVR1aVYcuX7580uFIkrSozXrp0qq6oao+3d49uaoeCRxKU7xYSXM+6C0LFqEkSdL82ZJkOUB7u7lt30ST18xY0bYNa5ckSQto1mLFEN9ob99YVSuBx89TPJIkSQvpfOCE9u8TaCbZnGlfk2T3JCuAg4CvVNV/APcleVzb70U9z5EkSQtkR3NWbKOqegsbBVwL3Nk+9o2BT5IkSZqQJOuAI2jm39oE/FG7rE/yW8CNwAsBqmpjknOBy4H7gJf3TDD+UuC9SR5IM6/F3411QyRJWoJGLlb0qqr7gJ+Y51gkSZLmTXsFs0GOGtL/dOD0Ae0bgUPmMTRJHTBsYtMNazeMORJJg+zsaSCSJEmSJEkLwmKFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjrFYoUkSZIkSeoUixWSJEmSJKlTlk06AEmSJGkpWL1u9aRDkKRFw5EVkiRJkiSpUxxZIUmSJEmtYSNgNqzdMOZIpKXNkRWSJEmSJKlTLFZIkiRJkqROsVghSZIkSZI6xWKFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjpl2aQDkCRJkqbJ6nWrJx2CJC16jqyQJEmSJEmdYrFCkiRJkiR1iqeBSJIkSdIshp3es2HthjFHIi0NjqyQJEmSJEmdYrFCkiRJkiR1isUKSZIkSZLUKRYrJEmSJElSp1iskCRJkiRJnWKxQpIkSZIkdYrFCkmSJEmS1CkWKyRJkiRJUqdYrJAkSZIkSZ2ybNIBSJIkSdJitXrd6oHtG9ZuGHMk0nRxZIUkSZIkSeoUixWSJEmSJKlTPA1E0qJ10tmXTjoESZIkSQvAYoW0q85ZM+kIJEmSJGmqWKyQJEmSdsKwiRUlSbvOOSskSZIkSVKnWKyQJEmSJEmdYrFCkiRJkiR1isUKSZIkSZLUKRYrJEmSJElSp1iskCRJkiRJneKlSyVJkiRpng27tO2GtRvGHIm0OFmskCRJknZg2H86JUkLx9NAJEmSJElSp1iskCRJkiRJnWKxQpIkSZIkdYrFCkmSJEmS1ClOsCktVuesGdx+/PrxxiFJkiRJ88xihSRJkoRX/ZCkLrFYIUmSJEljMqwotmHthjFHInWbc1ZIkiRJkqROcWSFNBfD5omQJEmSJM2bsY6sSPLeJJuTXNnTtleSC5NckeSzSR7Z89gbknwjyZVJnj3OWCVJkiRJ0mSM+zSQs4Gj+9reAlxQVb8IXNDeJ8njgWOBg9vnvDvJg8YXqiRJkiRJmoSxFiuq6hLg5r7mY4APtH9/sL0/076+qu6uqk3AVcATxhKoJEmSJEmamC7MWbG8qrYAVNWWJPu27SuAi3v6bWrbtpPkFOAUgFWrVi1gqJKWkpPOvnRg+1knHjbmSCRJkqSlpQvFil1WVWcCZwIceuihNeFwJEmSJGlOvKSptK0uXLp0S5LlAO3t5rZ9E7Cyp9+Ktk2SJEmSJE2xLhQrzgdOaP8+gWaSzZn2NUl2T7ICOAj4ygTikyRJkiRJYzTW00CSrAOOAPZJsgn4o3ZZn+S3gBuBFwJU1cYk5wKXA/cBL6+qO8cZryRJkiRJGr+xFiuqau2Qh44a0v904PSFi0iSJEmSJHVNF04DkSRJkiRJup/FCkmSJEmS1CkWKyRJkiRJUqdYrJAkSUtOkmuSXJHka0k2tm17Jbmwbf9skkf29H9Dkm8kuTLJsycXuSRJS8NYJ9iUJEnqkGdU1U09998CXFBVf57kte39Vyd5PHAscDCwH/DFJD/nVcokjcPqdasHtm9Yu2HMkUjj5cgKSZKkxjHAB9q/P9jen2lfX1V3V9Um4CrgCROIT5KkJcORFZIkaSkq4MIky4Azq+odwPKq2gJQVVuS7Nv2XQFc3PPcTW2bFqlhv1RLkrrDYoUkSVqKDq+qG9qCxKeTXL2rL5jkFOAUgFWrVu3qy0mStKR5GogkSVpyquqG9nYz8FHgMGBLkuUA7e3mtvsmYGXP01e0bf2veWZVHVpVhy5fvnwhw5ckaepZrJAkSUtKkock2WPmb+Bo4OvA+cAJbbcTgAvav88H1iTZPckK4CDgK+ONWpKkpcXTQCRJ0lKzH3BekgL2ANYDHwe+AKxP8lvAjcALAapqY5JzgcuB+4CXeyUQSZIWlsUKSZK0pFTVv9NchrTf94CjhjzndOD0hYxLkiRt5WkgkiRJkiSpUyxWSJIkSZKkTvE0EEmSJE2l1etWTzoESdJOcmSFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTnHOCkmSJElaZIbNybJh7YYxRyItDEdWSJIkSZKkTrFYIUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVOcYFOaNuesGdx+/PrxxiFJkiRJO8mRFZIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjrFYoUkSZIkSeoUJ9iUJEnSorZ63epJhyBJmmeOrJAkSZIkSZ1isUKSJEmSJHWKxQpJkiRJktQpzlkhSZIkSVNi2BwuG9ZuGHMk0q5xZIUkSZIkSeoUixWSJEmSJKlTLFZIkiRJkqROcc4KSZ1x0tmXTjoESZKkqeRcFlpsHFkhSZIkSfBaNnoAABbOSURBVJI6xWKFJEmSJEnqFE8DkQY5Z82kI5AkSZKkJctihSTN0bC5Nc468bAxRyJJkiRNJ08DkSRJkiRJneLICkmSJC0Kw65mIEmaPhYrpKVi2Dwcx68fbxySJEmSNAtPA5EkSZIkSZ3iyApJkiRJWqKGnV61Ye2GMUcibcuRFZIkSZIkqVMsVkiSJEmSpE7xNBBJkiRJ0jYGnR7iqSEaJ0dWSJIkSZKkTrFYIUmSJEmSOsXTQCRNxElnXzrpECRJkiR1lMUKLW3nrJl0BJIkSZKkPhYrJEmS1CmDJvaTJC0tFiskSZIkSbMaVkj0KiFaCE6wKUmSJEmSOsVihSRJkiRJ6hSLFZIkSZIkqVMsVkiSJEmSpE6xWCFJkiRJkjrFYoUkSZIkSeoUixWSJEmSJKlTlk06AEkTds6awe3Hrx9vHJIkSZLUslghSfPkpLMvHdh+1omHjTkSSZIkaXGzWCFJkiRJ2mmr160e2L5h7YYxR6Jp4pwVkiRJkiSpUxxZIUmSpIkY9muspOngiAvtCosVkhbUsHkcJEmSJGkYixWS5oVFCUmSJEnzxTkrJEmSJElSp1iskCRJkiRJneJpIJK0wIadInPWiYeNORJJkiRpcbBYoaXhnDWTjmDxGbbPjl8/3jgkSZI0VbxKiEZhsUKSJEmSNHEWMdTLOSskSZIkSVKnOLJC0twMPaXmdWMNQ5IkSdL0cmSFJEmSJEnqFEdWSJIkaUENOw9dkkbhXBZLk8UKTRev+rHgvnbdrYMf2G+8cUyDYZc0HWYulzr1cqmSJGmpsrgxHSxWSJoXr7rxTQPb37HfaWOORJIkSUuBo7amm8UKSZIkSdLUc8TF4mKxQtJAQ0/3kCRJkqQFZrFC0oKay+khnkoiSZIkCSxWaLFyIs05c6SEJEmSpMXCYoW0SFl8kCR1jZPdSVqMnMuimyxWSFq0PG2kMddLoEqSJGl2cy1iWPSYX4uiWJHkaOAMYDfgb6rqTyYckiRJWkLMRbbnKApJS5Xff+PR+WJFkgcB7wKeCtwAfCnJZ6vqsslGJs2vYad1HLLyEWOOZDyGjYrQcPMxgmLYa5x14mG7/Nrz+foLHac0F+YikqRdMai44eiM2XW+WAE8Ebiqqq4DSLIeOAYwQViMhk2Mefz6ufUfokv/4Z+vOSWcm2J6LaXTWLpWfOhaPOo8cxFJ0rya6+iMpVjEWAzFihXAdT33NwFHTCYULRVdKnpo7qMwllIRYN4MKQyedNfrBrbP23/qhxYkB69XmpBO5CLzlag6fFmSpsdCf6dPshiSqprYykeR5HjgaVX18vb+WuCIqnpZT59TgFPauz8HfHOew9gHuGmeX1PDub/Hz30+Xu7v8XJ/b7V/VS2fdBCLzSi5SNtuPjI3bk+3TdP2TNO2gNvTdW7P7EbKRxbDyIpNwMqe+yvatvtV1ZnAmQsVQJKNVXXoQr2+tuX+Hj/3+Xi5v8fL/a15MGsuAuYjc+X2dNs0bc80bQu4PV3n9syfB0xipXP0FeCgJCuS7A6sAS6YcEySJGnpMBeRJGnMOj+yoqp+lOQVwGdoiisfrKqNEw5LkiQtEeYikiSNX+eLFQBVdT5w/gRDWLAhnRrI/T1+7vPxcn+Pl/tbu6wDuQhM32fZ7em2adqeadoWcHu6zu2ZJ52fYFOSJEmSJC0ti2HOCkmSJEmStIRYrGgleW+SzUmuHPL47yX5WrtcmeTeJHuNO85pMcL+flSSi5J8Pcm3krx83DFOmxH2+d5JLmj3+VeSHDTuGKdJkpVJLmm/L76V5NQBfZLk7e0+/5ckj5tErNNgxP19YJIvJbkzyesmEac0m9m+q9s+RyS5NMm/JrlknPHN1bTlV9OWv0xTbjBtx91pO66NuD0vTnJF2+erSTp7RY0Rt+e57fZc3vZ7ziRiHcUo29PT97Ak9yR5/oLH5WkgjSRPA24D3l9VO/wiTrIaeG1VPXMswU2h2fZ3ktOA3avq1CTLgW8DP15Vd4w51Kkxwj5/B3BTVb0lyYHA+6rq8HHHOS2SPArYt6ouT/JQ4DLgBVX1tZ4+xwIvAX4DeCzNPv+liQS8yI24v/cF9qfZ37dU1RmTiVYaboTv6kcBFwFHVtUNSfapqpvGHeeopi2/mrb8ZZpyg2k77k7bcW3E7XkicHVVfb/9j/1bq+qQCYW8QyNuz57A7VVVSQ4GPllVqyYU8g6Nsj1tv92AC4EfAe+tqo8uZFyOrGhV1SXAzSN2XwusW8Bwpt4I+3sT8NAkAfYEbgLuHEds02qEfX4gcHHb92pg3ySPHkds06iqbqiqy9u//xu4HOjfn8fQXFWgquoyYFmSlWMOdSqMsr+ranNVXQrcPYEQpZGM8F19HPDhqrqh7d/ZQgVMX341bfnLNOUG03bcnbbj2ojb889V9f327hf7H++SEbfntto6MuAhwA3jjXJ0I/77AXgV8HfA5nHEZbFijpLsARxN8yZp4bwHeAzwX8AVwO9W1X2TDWnqXQH8JkCSJ9BU6jtZ/V1skhwAHEZz4O21Ariu5/6mtk27YAf7W5oGBwI/nuTL7fDi3550QPNhivKractfFmVuMG3H3Wk7ro24PS8DPjGOeHbVjrYnyfOSXA18Gnj1eCPbOcO2py1UPg/4q3HFYrFi7lYD/1hVo/5KoJ3zBpqK3v8ADgHemeRhkw1p6r0F2C/J14FTgY2A54ntonYI4EeB1/T8WqAF4v7WEvAAmuPikcAzgFO7PI/AHExLfjVt+cuiyw2m7TiwFLcnyRHAScDrxxjaTplte6rq3Ko6kOY77v1JOv3/71m25/8Cp46zALtsXCuaIsfR8SGKU+KpwGnt0KnvJPkuzS8VX55sWNOr/UI6fuZ+kn8HvjW5iBa/JLvT/Eq4rqo+NqDLJmAlWz/XK9o27YQR9rc0Da4Drq+q24Hbk/wDcDAwdELORWJa8qupyl8WW24wbcfdaTuujbI97dwOZwHPqarvjTO+uZrL+1NVlyRZBuwHXD+O+OZqhO05FPhQc5Yb+wC/muSeqjpvoWLqdGWna5I8HHg68PFJx7IE/BvNr0Yk2Y/mQH/NJAOadkke3n6JkuQE4F+m4BeuiWnPVz4L+EZV/dmQbucDL2r7Pw64r6quG9JXOzDi/pamwaeApyRZ1p46cThw9YRj2iVTll9NVf6ymHKDaTvuTttxbZTtSbIK+Bjw4qrqbFEMRt6en+j5+3HAgxjTXA9zNcr2VNVPVNUBVXUAzeiL31nIQgV4NZD7JVkHHEFTJboR+CNgd4Cqelfb50Tg6Ko6bjJRTo/Z9nc7I+0HaSZ22Q14W1W9ZzLRTocR9vmTgLNpZvf9DnBSVd0ykWCnQJKnAF+gOd93ZrjcG2nP9W33eYB30gzlvgs4uao2TiDcRW/E/f0omiHMD2v73AY8pqp+MP6IpcFGzEd+D3hp235WVf3JRIIdwbTlV9OWv0xTbjBtx91pO66NuD1/DRwLXNs+fk9VdfLypSNuzxtoi2M0/4ZeW1VfGHesoxhle/r6n01zdZMFvRqIxQpJkiRJktQpngYiSZIkSZI6xWKFJEmSJEnqFIsVkiRJkiSpUyxWSJIkSZKkTrFYIUmSJEmSOsVihTopyZuTVM/y/SSfTvJLk45tVyR5YLtth4x5vfu26z1gHl/zjCTXzNKn9328L8ktSS5Ncnp7ua0FleSadt1vGvDYU3piO2ChYxlFkjck+cyA9qcn+XiSzUnubm8/leS4JCN/jyfZkOSKHTz+ziS3JnlQktVJ/i3JA3d2eyRpsTMfmff1mo9s/5j5yPaPm48IsFihbvs+cHi7vBRYCVyU5JETjWrXPJDmGuZjTQ6Afdv1HjDm9cLW9/FJwHHAx4AXA1ckefwY1n9bu95+a9vHOiHJI4BTgdP72l8DfB64F3gVcCTwSuAHwN/SXCt+VOuAg5I8ZsD6dwOeD3ysqu6sqg3AD4GT5741kjRVzEfmj/nI9sxHtl2P+YjuZ7FCXXZPVX25XT4GnADsDfz6hOMaiyQ/NukY5knv+/iZqnorcDBwPfCh9qC0kD4JPCbJQTMNPQfCTyzwuufiZOD6qrpkpiHJ44AzgD+uqt+sqvVVdUlVfbiq1gJPAW6awzo+TnPAXzvgsWcA+9EkEDPeA7w2Sea4LZI0TcxHpoP5yGjMR9QZFiu0mMwMF3t0b2OSvZKcmeTGJD9K8k9JntjXZ7d2SNu3ktzVDlv7274+r0zy7SR3JvlOktf2Pf7mJDcleWySL7frujrJs/r6vTDJ5e16bk9yWZKZavN/t7fv6x3y1y6V5EVJ3p/kVmBD+3qV5JWDYulr2z/JujbGu5J8PcmL2yGFM/vu8zPrneP+e0SSc5LcluT6JH+w/dszuqq6FXg98NPAr4zynCSHJPl8u0/vat/L3x3hqf8JfJFtf814JrAnA5KDJL+f5F+S3NFu74VJDu7r88wk/9zurzuSXJlkTc/jO/oMDPNi4Ly+tlcBm4HTBj2hqr5UVf/aF9vJSa5qP8fXJnl9T//baT5Xa/pfi2b/bAYu7mn7OM179MuzxC5JS4n5SF8sfW3mI4OZj2ztbz6ikVis0GKysr3dMtOQ5EHA54CnAa8GfhXYBHwuyU/2PPfdwB8CfwMcBbwM6D1AvgZ4O7AeeFbb74wkv98Xwx7AX7d9j2lj+UiSh7ev8/M0leBP0Rz0fh34EPCw9vnPbG9PY+uQ0ut7Xv9tNF/Qv86QA8IgSfYFvgT8Is2QvCOBd9IMt7weeFHb9X/2rHcu++9D7Wu+AnhJ23/QUMa5+HvgHkY48KT55eGTNEM4f7ON5W3A7iOuax3bxruW5iB5+4C+ewF/QfM5eAHN0MzP9bzHe9EkFVfQfAaeQ/P5ekj7+GyfgUHbtx/Nrztf7nvoacDFVXXPKBuZ5PeA/9eu71eAPwf+uDdBaGP7mfQMeU2yO81+/XBV3TvTXlXXAjcAzx5l/ZK0RJiPDGE+Mivzka3MRzS7qnJx6dwCvJlmONmydllJc3C4Hdivp99JwJ3A/j1tuwFXA29v7x9Ikwj89pB17dau61197X9BczB6cE9MBTy5p88vtG3Pbe+/CNiyg+3as+1/Yl/7AW37hwY8p4BXDto/PfffCtwKLB+y3oPa1zmir32U/ffY9rm/0dNnD+BG4JpR3scdPH498FcjfB4e3cbwC3P8HF1DM2xxOXA3cBjNebq3AL8B/Fr7ugfs4LPx4Lb/S9q2J7fPeciQ5+zwMzDkOUe1r/kzfe13AG/ta0vPv4tlwAPa9ofRJDJv6Ov/JuB7wLL2/sz2v62nz8x+eNKA2C4Ezp3L9ri4uLhMy4L5SO9j5iPmI71t5iMuC744skJdtjfNF/rdwH8ATweOqaobe/ocRVP9/c8ky5Iso/ny/Hvaaj3NuW/30kz+M8iB7bo+2tf+YZov3F/safthVf1jz/2r29sfb2//Fdg7yfuS/EqSPUfZ0B6fmmP/Gc8Ezq+qLbP23NYo++9JNL84nD/zpKr6IfDZnYy116jnHt5Im0gkeUFb+R9Zu18upvk14+h2vRcMDCg5IsklSX5Is913AI8Afrbt8s227Zw0M1T3T7C2M5+Bme25eVD4ffePZeu/i7uB/9O2H07za8pHZt7L9v28mObXmZ8DqKq7aCYVe2Fy/7mfa4BraX4N6/e9nvgkaSkyHxmd+cgOmI+Yj2huLFaoy75PU3n+ZZphkvfSDPnrtQ/N0LS7+5aX0RzwaW//uz2gDTLz5d5/YJ25v1dP2x29HWrrELVl7f0raYavHUhz8Lk5yUfmcDC7ZcR+/famOYDO1Sj775HAbe1BpddcJlLaTpIHM2Lc1Qw7fDbN/jkbuCHNebpPmMMqPwS8EDgeOK+q7hwQ008Bn6Y5l3ctzWfvsDbGB7ex3NTGsidNQnlTks8m+Zn28V35DPQnS/8FrOhru6iN6TC2HbK7T3v7bbZ9L2eS2b17+q4DVgGHt+/Dc2l+RetPRAbFJElLjfnI6MxHZmc+spX5iHZo2aQDkHbgnqra2P79z0luA/42yTlV9bm2/WaaL7/XDHj+zJf/94CHJtljSIIwc0Be3tc+c39QdXmoqjoPOC/Jw2jOH3wH8C7geXN5nR730gz/6/WQvvs7W20eZf/dCuyZ5IF9CcI+A54zF8+g+Q4aVD3fTlVdATy3PZ/xycCfAp9M8uiqunuElziX5n14Ac25nYOspjkYPr+q7oD7z0/d5vzOqvoCcGSaGdKfAfwlzS9fj20fn+tn4Ib2dm+2TbouAZ6VZLeZRLSqbgE2trH1vh8zn9NnMTjJ/GbP35+nSXiOo/kV7qFsO+t2r7164pOkpch8pGE+gvmI+YjGyZEVWkzWAVfRXJ97xkXAzwPfrqqNfcvMjNMX0xxcB10eCZqhkzfRVJ97PZ/m2tFXbPeMEVTVD6pqPfB3wMx1pGe+zEediAmaavXMkD/aoXJH9vW5CHhOkr0ZbNh6R9l//0RzEP/Vnhj2oDkI7ZQ01/D+U+A7NBNqjayq7q6qv6eZrGk5IyYp1cz4/ac078ewdT6QJjm4r6ftucDAy7ZV1R1VdT7NJGfbXSt8yGdgkJl9/Qt97e+gSfreuIPnzvgSzS9tjxrwXm6sqpmZ32d+gfswTaJ0PPCN6pvFu8dBwOUjrF+SlgrzEcxHzEeGMh/RvHFkhRaNqqok/5vm14ynttXk9wMvB/4hyRnAd2nO5zucZlKhv6yqbyY5E3hHO0v1F2mqs8dW1Uuq6t4kpwN/lubyW58Dngr8LvCmqvrRqDEm+W3g8cBnaBKOn6VJMs5tt+GuJN8Fnp/kKpqD9mxfvJ8ATkzy1Xb7Tmb7A+Jf0AxJ/XyS02gSiscAe1bVn9GcY3sH8OIkPwDubX8lGmX/XZbkQuDd7fmONwKvoxnSN4plSWZm2H5ou39eQTMp1tE9Q1eHSnOprtNohjl+lyYpeD1NdX7kKntV/eEsXS4C/oTmUm5/DfwU8Ps0v+bMxHIMzQH1EzTDIlcCpwD/0D6+w8/AkLg2J/kazfm4H+tpvyzJ64A/T3IIzezw/0XzS9aTgUfRTGJFVd2a5M3Au9JcHu4LNInOzwNPr6r+y4Oto7kU2fPYNuG+X5L923XMx/nAkjQVzEfMRzAfMR/ReFQHZvl0celfGDJrM80vEt8CLuhpezjNsLfraA5YW2hm6n5q3/PeCPx72+dG4AN9r/0qmsr6XW2/144Y0/2zY9N8uX+6jeEemi/ytwN79PQ/huagdk/73APYOvv2rw14/YcDHwF+SHMgfBPwlv5YgP1pDh63tNtwFfCinsdfSjNp0b3NP/057b9H0pxjeXu77/6QZlbra0Z4H6td7qM5yG4ETqepuI/6ediP5mB2bbvfbqE52P7kLM+7BjhjB49vN/s2zYH+P4Ef0Uz29cTe16E59/O89r29p91ff0M78/kon4Ehsfwv4LtDHjuCJhnZ0r5Hm2kmGDsOSF/fE4Cv0iSDP2z396lDXve77fb/9JDHX03zbyI7it3FxcVlWhfMR+jbPvMR8xHzEZexLWk/AJKkCWqHol5D8wvbRRMOB4AklwNnVtU7Jx2LJElaeOYj6hLnrJCkDqit57D+waRjAUiymmaG8fdMOhZJkjQe5iPqEkdWSJq4JA9gB8XTai4VJkmStGDMR6RucWSFpC54L9tfW/3+pZ2gSZIkaSGZj0gd4sgKSRPXHvx3dMmvy2vba6pLkiTNK/MRqVssVkiSJEmSpE7xNBBJkiRJktQpFiskSZIkSVKnWKyQJEmSJEmdYrFCkiRJkiR1isUKSZIkSZLUKf8frC6Su5uqxnIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1296x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"if l_flv[l_index]=='mu':\n",
" lower_MC = 1.75\n",
" upper_MC = 2.1\n",
" lower_data = 1.65\n",
" upper_data = 2.15\n",
" \n",
"if l_flv[l_index]=='e':\n",
" lower_MC = 1.7\n",
" upper_MC = 2.15\n",
" lower_data = 1.0\n",
" upper_data = 3.0\n",
" lower_data = 1.6\n",
" upper_data = 2.4\n",
"\n",
"#plt.suptitle('Reconstructed D_s mass', fontsize=15)\n",
"plt.subplot(1,2,1)\n",
"label=\"Ds_ConsD_M\"\n",
"plt.title('MC', fontsize=15)\n",
"left_h_Ds=[MC_Ds_tuple_dict_presel_9[label][i][0]/1000 for i in range(len(MC_Ds_tuple_dict_presel_9[\"Ds_ConsD_M\"]))]\n",
"label=\"Dplus_ConsD_M\"\n",
"left_h_Dplus=[MC_Dplus_tuple_dict_presel_9[label][i][0]/1000 for i in range(len(MC_Dplus_tuple_dict_presel_9[\"Dplus_ConsD_M\"]))]\n",
"plt.hist(left_h_Ds,alpha=0.7,bins=70,range=(lower_MC,upper_MC), label='MC after presel cuts',density=False);\n",
"plt.hist(left_h_Dplus,alpha=0.7,bins=70,range=(lower_MC,upper_MC), label='MC after presel cuts',density=False);\n",
"plt.legend(fontsize='12')\n",
"plt.ylabel('# events (a.u.)', fontsize=15)\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"\n",
"label=\"Ds_ConsD_M\"\n",
"plt.subplot(1,2,2)\n",
"plt.title('Data', fontsize=15)\n",
"right_h=[data_tuple_dict_presel_4[label][i][0]/1000 for i in range(len(data_tuple_dict_presel_4[\"Ds_ConsD_M\"]))]\n",
"plt.hist(right_h,alpha=0.7,bins=70,range=(lower_data,upper_data), color='green',label='Data after presel cuts',density=False);\n",
"plt.legend(fontsize='12')\n",
"plt.xlabel('Reconstructed D_s Mass (GeV)', fontsize=15)\n",
"fig=plt.gcf()\n",
"fig.set_size_inches(18,8)\n",
"plt.savefig('/home/hep/davide/Rphipi/'+l_flv[l_index]+'_dataMC_after_presel_and_phi_mass_cut.png', format='png', dpi=100)"
]
},
{
"cell_type": "code",
"execution_count": 108,
"metadata": {},
"outputs": [],
"source": [
"Dplus_af_phicut=np.float(MC_Dplus_tuple_dict_presel_9[\"Dplus_ConsD_M\"].shape[0])\n",
"Ds_af_phicut=np.float(MC_Ds_tuple_dict_presel_9[\"Ds_ConsD_M\"].shape[0])"
]
},
{
"cell_type": "code",
"execution_count": 109,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
" Phi mass cut efficiency on signal \n",
" Dplus eff 0.964, Ds eff 0.9657\n"
]
}
],
"source": [
"eff_phi_cut_Ds=Ds_af_phicut/Ds_bf_phicut\n",
"eff_phi_cut_Dplus=Dplus_af_phicut/Dplus_bf_phicut\n",
"print(\"\\n Phi mass cut efficiency on signal \\n Dplus eff {0:.4g}, Ds eff {1:.4g}\".format(eff_phi_cut_Dplus,eff_phi_cut_Ds))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Boost Cross Checks"
]
},
{
"cell_type": "code",
"execution_count": 110,
"metadata": {},
"outputs": [],
"source": [
"## pD = ppi + q in D rest frame\n",
"#\n",
"## Lambda D to lab:\n",
"#\n",
"#p_Dplus_lab, _, _, _, = get_4_momenta(MC_Dplus_tuple_dict_presel_9, 'Dplus')\n",
"#p_pi_lab, _, _, _,= get_4_momenta(MC_Dplus_tuple_dict_presel_9, 'pi')\n",
"#p_l_plus_lab, _, _, _, = get_4_momenta(MC_Dplus_tuple_dict_presel_9, l_flv[l_index]+'_plus')\n",
"#p_l_minus_lab, _, _, _, = get_4_momenta(MC_Dplus_tuple_dict_presel_9, l_flv[l_index]+'_minus')\n",
"#\n",
"#beta_D, beta_D_vec = get_boost_D_to_lab(MC_Dplus_tuple_dict_presel_9, mother_index=1)\n",
"#\n",
"#sigma_Es = []\n",
"#sigma_pxs = []\n",
"#sigma_pys = []\n",
"#sigma_pzs = []\n",
"#\n",
"#D_mass = []\n",
"#pi_mass =[]\n",
"#l_mass =[]\n",
"#\n",
"#for j in range(len(MC_Dplus_tuple_dict_presel_9[\"Dplus_ConsD_M\"])):\n",
"# \n",
"#\n",
"# beta_D=r.TVector3(beta_D_vec[:,j][0], beta_D_vec[:,j][1], beta_D_vec[:,j][2])\n",
"# lambda_lab_to_D=r.TLorentzRotation(-beta_D)\n",
"#\n",
"# p_Dplus_lab_r = get_TLorentzVector(p_Dplus_lab, j)\n",
"# p_pi_lab_r = get_TLorentzVector(p_pi_lab, j)\n",
"# p_l_plus_lab_r = get_TLorentzVector(p_l_plus_lab, j)\n",
"# p_l_minus_lab_r = get_TLorentzVector(p_l_minus_lab, j)\n",
"#\n",
"# #boost it\n",
"# p_Dplus_D = lambda_lab_to_D.VectorMultiplication(p_Dplus_lab_r)\n",
"# p_pi_D = lambda_lab_to_D.VectorMultiplication(p_pi_lab_r)\n",
"# p_l_plus_D = lambda_lab_to_D.VectorMultiplication(p_l_plus_lab_r)\n",
"# p_l_minus_D = lambda_lab_to_D.VectorMultiplication(p_l_minus_lab_r)\n",
"#\n",
"# #impose momentum conservation in D lab frame\n",
"# total_p=p_l_minus_D+p_l_minus_D+p_pi_D\n",
"# sigma_E, sigma_px, sigma_py, sigma_pz = [total_p.E()-p_Dplus_D.E(), total_p.Px()-p_Dplus_D.Px(), total_p.Py()-p_Dplus_D.Py(), total_p.Pz()-p_Dplus_D.Pz()]\n",
"# \n",
"# D_mass.append(p_Dplus_D.M2())\n",
"# pi_mass.append(p_pi_D.M2())\n",
"# \n",
"# l_mass.append(p_l_minus_D.M2())\n",
"# l_mass.append(p_l_plus_D.M2())\n",
"# \n",
"# sigma_Es.append(sigma_E)\n",
"# sigma_pxs.append(sigma_px)\n",
"# sigma_pys.append(sigma_py)\n",
"# sigma_pzs.append(sigma_pz) \n",
"# \n",
"#sigma_Es = np.array(sigma_Es)\n",
"#sigma_pxs = np.array(sigma_pxs)\n",
"#sigma_pys = np.array(sigma_pys)\n",
"#sigma_pzs = np.array(sigma_pzs)\n",
"#\n",
"#D_mass = np.sqrt(np.array(D_mass))\n",
"#pi_mass = np.sqrt(np.array(pi_mass))\n",
"#l_mass = np.sqrt(np.array(l_mass))\n",
"#\n",
"#plt.subplot(2,2,1)\n",
"#plt.hist(sigma_Es, bins=50);\n",
"#plt.subplot(2,2,2)\n",
"#plt.hist(sigma_pzs, bins=50);\n",
"#plt.subplot(2,2,3)\n",
"#plt.hist(sigma_pxs, bins=50);\n",
"#plt.subplot(2,2,4)\n",
"#plt.hist(sigma_pys, bins=50);\n",
"#fig = plt.gcf()\n",
"#\n",
"#\n",
"#fig.set_size_inches(16,8)\n",
"#print(D_mass.mean(), pi_mass.mean(), l_mass.mean())"
]
},
{
"cell_type": "code",
"execution_count": 111,
"metadata": {
"scrolled": false
},
"outputs": [],
"source": [
"#beta_q = get_boost_D_to_q(MC_Dplus_tuple_dict_presel_9, l_index=l_index, mother_index=1)\n",
"#\n",
"#sigma_Es = []\n",
"#sigma_pxs = []\n",
"#sigma_pys = []\n",
"#sigma_pzs = []\n",
"#\n",
"#D_mass = []\n",
"#pi_mass =[]\n",
"#l_mass =[]\n",
"#\n",
"#for j in range(len(MC_Dplus_tuple_dict_presel_9[\"Dplus_ConsD_M\"])):\n",
"# \n",
"# \n",
"# beta_D=r.TVector3(beta_D_vec[:,j][0], beta_D_vec[:,j][1], beta_D_vec[:,j][2])\n",
"# lambda_lab_to_D=r.TLorentzRotation(-beta_D)\n",
"#\n",
"# p_Dplus_lab_r = get_TLorentzVector(p_Dplus_lab, j)\n",
"# p_pi_lab_r = get_TLorentzVector(p_pi_lab, j)\n",
"# p_l_plus_lab_r = get_TLorentzVector(p_l_plus_lab, j)\n",
"# p_l_minus_lab_r = get_TLorentzVector(p_l_minus_lab, j)\n",
"#\n",
"# #boost it from lab to D\n",
"# p_Dplus_D = lambda_lab_to_D.VectorMultiplication(p_Dplus_lab_r)\n",
"# p_pi_D = lambda_lab_to_D.VectorMultiplication(p_pi_lab_r)\n",
"# p_l_plus_D = lambda_lab_to_D.VectorMultiplication(p_l_plus_lab_r)\n",
"# p_l_minus_D = lambda_lab_to_D.VectorMultiplication(p_l_minus_lab_r)\n",
"# \n",
"# beta_q_r=r.TVector3(0,0, beta_q[j])\n",
"# lambda_D_to_q=r.TLorentzRotation(beta_q_r)\n",
"#\n",
"# #boost it from D to q\n",
"# p_Dplus_q = lambda_D_to_q.VectorMultiplication(p_Dplus_D)\n",
"# p_pi_q = lambda_D_to_q.VectorMultiplication(p_pi_D)\n",
"# p_l_plus_q = lambda_D_to_q.VectorMultiplication(p_l_plus_D)\n",
"# p_l_minus_q = lambda_D_to_q.VectorMultiplication(p_l_minus_D)\n",
"#\n",
"# #impose momentum conservation in D lab frame\n",
"# total_p=p_l_minus_q+p_l_minus_q+p_pi_q\n",
"# sigma_E, sigma_px, sigma_py, sigma_pz = [total_p.E()-p_Dplus_q.E(), total_p.Px()-p_Dplus_q.Px(), total_p.Py()-p_Dplus_q.Py(), total_p.Pz()-p_Dplus_q.Pz()]\n",
"#\n",
"# sigma_Es.append(sigma_E)\n",
"# sigma_pxs.append(sigma_px)\n",
"# sigma_pys.append(sigma_py)\n",
"# sigma_pzs.append(sigma_pz) \n",
"#\n",
"# D_mass.append(p_Dplus_D.M2())\n",
"# pi_mass.append(p_pi_D.M2())\n",
"# \n",
"# l_mass.append(p_l_minus_D.M2())\n",
"# l_mass.append(p_l_plus_D.M2())\n",
"# \n",
"# \n",
"#sigma_Es = np.array(sigma_Es)\n",
"#sigma_pxs = np.array(sigma_pxs)\n",
"#sigma_pys = np.array(sigma_pys)\n",
"#sigma_pzs = np.array(sigma_pzs)\n",
"#\n",
"#D_mass = np.sqrt(np.array(D_mass))\n",
"#pi_mass = np.sqrt(np.array(pi_mass))\n",
"#l_mass = np.sqrt(np.array(l_mass))\n",
"#\n",
"#plt.subplot(2,2,1)\n",
"#plt.hist(sigma_Es, bins=50);\n",
"#plt.subplot(2,2,2)\n",
"#plt.hist(sigma_pzs, bins=50);\n",
"#plt.subplot(2,2,3)\n",
"#plt.hist(sigma_pxs, bins=50);\n",
"#plt.subplot(2,2,4)\n",
"#plt.hist(sigma_pys, bins=50);\n",
"#fig = plt.gcf()\n",
"#\n",
"#\n",
"#fig.set_size_inches(16,8)\n",
"#print(D_mass.mean(), pi_mass.mean(), l_mass.mean())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cos theta_l before cut on phi mass"
]
},
{
"cell_type": "code",
"execution_count": 112,
"metadata": {},
"outputs": [],
"source": [
"costheta_Dplus_b=get_costheta_list(MC_Dplus_tuple_dict_presel_8, l_index=l_index, mother_index=1)\n",
"costheta_Ds_b=get_costheta_list(MC_Ds_tuple_dict_presel_8, l_index=l_index, mother_index=0)\n",
"costheta_MC_b=np.concatenate((costheta_Dplus_b,costheta_Ds_b),axis=0)\n"
]
},
{
"cell_type": "code",
"execution_count": 113,
"metadata": {},
"outputs": [],
"source": [
"costheta_data_b=get_costheta_list(data_tuple_dict_presel_3, l_index=l_index, mother_index=0)"
]
},
{
"cell_type": "code",
"execution_count": 114,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAHVCAYAAAAJnF2uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X2UX1V9L/73xoSHQn9LK7QgA0a9S3lIckdMBEGQQF0gYLi2SMCHkhrkhovXQq/eiv0JlbY3dpmfWEsx0mKTSguh0EIsUEQIFaoYAgwhBLC1RAhUQQoqIiQk+/dHhrmZZCaZJJPZ8/B6rZWV+Z5z5nw/Z84k+7y/e599Sq01AAAA0MJOrQsAAABg7BJKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoRSgEAAGhGKAUAAKAZoRQAAIBmhFIAAACaGdfqjffcc886YcKEVm8PwChzzz33/LjWulfrOkYybTMAg2mgbXOzUDphwoQsXbq01dsDMMqUUn7QuoaRTtsMwGAaaNts+C4AAADNCKUAAAA0I5QCAADQTLN7SvuyZs2arFq1Ki+++GLrUthKu+66azo6OjJ+/PjWpQAAsIO5bmdD25sFhlUoXbVqVX75l385EyZMSCmldTkMUK01zzzzTFatWpU3vOENrcsBAGAHc93OKwYjCwyr4bsvvvhiXvva1/rFHmFKKXnta1/rkzIAgDHCdTuvGIwsMKxCaRK/2COU8wYAMLa4/uMV2/u7MOxCKQAAAGPHsLqndGOz5t89qPu7fObULW5TSskHP/jBXHHFFUmSl19+Ofvss08OPfTQ/OM//mOS5KabbsqFF16Yl156KWvWrMkxxxyTSy65pNd+5s+fn09+8pPp6OjI2rVrs88+++TCCy/M4YcfPqBau7q68uSTT+aEE07YyqP8v44++ujMnTs3U6ZM2eZ9AADAlozl6/b+vOpVr8qkSZNSSump9bzzzstOO22+X3CPPfbI888/v13vffvtt2fnnXfermOYMGFCli5dmj333HO7ahkIPaUb2X333bN8+fL84he/SJLccsst2XfffXvW33XXXTn33HNz9dVX5/77788DDzyQgw8+uM99zZgxI/fdd1+WLVuWz372s3n/+9+fhx56aEB1dHV15cYbb9z+AwIAgFFouFy3z5w5M7fffvsmy3fbbbd0dXXlvvvuy6233prbbrstn/3sZ7f+QLfB7bffnm9/+9tD8l6DQSjtwwknnJAbbrghSXLllVfm9NNP71l38cUX54ILLsiECROSrP8E5Oyzz97iPg877LCcc845ueyyy3r2c9BBB+W//tf/mlNPPbXXtqtXr84FF1yQhQsXprOzMwsXLswf/MEfZO7cuT3bTJw4MStXrszKlStzwAEH5MMf/nAOPvjgnHjiifn5z3++yfsvWrQohxxySCZOnJjp06fnpz/96Vb/XAAAYDhpfd0+UK95zWty2WWX5ZJLLkmtNfPnz8/JJ5+cY489Nm9605vy6U9/epPvuf3223PSSSf1vP7Yxz6W+fPnJ0k++clP9tR03nnn9fq+lStXZt68ebn44ovT2dmZO+64IzNnzsw111zTs80ee+zR8x5HHXVU3vve9+bNb35zZs6cmXXr1m1Sy1/8xV9k8uTJOeigg/Lbv/3bWbNmzTb9HPojlPbhtNNOy1VXXZUXX3wxy5Yty6GHHtqzbtmyZXnb2962Tfs95JBD8vDDDydJPv/5z+e+++7L/fffn8svv7zXdjvvvHMuuuiizJgxI11dXZkxY8Zm9/vII4/knHPOyYMPPpjXve51+dM//dNe63/0ox9lzpw5ueOOO7J8+fIcccQR+ZM/+ZNtOgYAABguWl+3b419990348ePz1NPPZUkWbJkSa699tqsWLEiX//61wfcs/nUU0/lxhtvzIMPPpj7778/F110Ua/1EyZMyOzZs3Peeeelq6srRx555Gb3t2TJklxyySV5+OGH8x//8R+56qqreq2///77c9111+Wee+7JihUrsttuu/WE48EilPZh8uTJWblyZa688srtuqdzY7XWXu/x4Q9/OH/913+93fvdb7/9cthhhyVJTj/99Nx555291t9xxx3513/91xxxxBHp7OzMggUL8sQTT2z3+wIAQEutrttvvvnmdHZ2prOzM4sWLcqZZ56Zzs7OXqF4S/t997vfnVe/+tXZZZdd8r73vW+Ta/j+/Mqv/ErGjx+fWbNm5ZprrsnOO++8lUfX29vf/va8/vWvz0477ZQZM2ZsUsctt9yS++67L1OnTk1nZ2duvfXWPPbYY9v1nhsTSvsxffr0fOITn+g1BCBJJk2alHvvvXeb9nnfffflwAMPTJLccMMNOfvss3P//fdn6tSpefnllzf7vTvttFOvrvQNnwO08RTMG7+uteY973lPurq60tXVlRUrVgz6pxsAANBCi+v24447rufaevr06fnLv/zLdHV15bvf/W6/+3zyySezdu3a/Oqv/mqSLV/D93f9P27cuHz3u9/NKaeckptuuinHH3/8Fo9nw32tW7cuq1ev7vd9+8oSs2bN6jneRx55JH/4h3+4xffcGkJpPz7ykY/kwgsvzKRJk3otP/fcc3PRRRflBz/4QZL1J3XevHlb3N+SJUty6aWX5qMf/WjWrVuXJ598MtOmTcvnPve5/PSnP81zzz3Xa/vddtstL7zwQs/rjo6Onn9UXV1defTRR3vWPfbYY1myZEmSZOHChTniiCN67evII4/M4sWLe2p+8cUX82//9m8D/VEAAMCw1fq6fSCee+65zJ49O+ecc05P6Lvlllvyk5/8JKtXr8511123yTV8R0dHHnzwwbz00kv52c9+lltvvTVJ8vzzz+f555/PCSeckC984Qt9Bu++ssQ999yTZH3I3vCe0CVLluSxxx5LrTV/93d/t0kd7373u3P11VfnP//zP5MkP/3pT/P4449v9c9gc4b1I2EGMhX0jtLR0ZGPf/zjmyw//PDDM3fu3JxyyilZvXp11q5dm2OPPbbPfSxcuDD/8i//knXr1mXvvffOwoULc+CBB2bNmjWZMWNGfv7zn2ft2rX5H//jf2wy1fK0adPyf/7P/8nkyZPz+7//+zn11FOzYMGCHHzwwTnssMPy5je/uWfbt7zlLfmzP/uz3HvvvXn961+fL3zhC732tffee+eyyy7L9OnTk6z/B/lHf/RH+S//5b9s748JAADG9HV7f37xi1+ks7Oz55EwH/jAB/K7v/u7Pevf/va35zd/8zfz6KOP5tRTT93k8S1vfOMbc/LJJ+eAAw7IW97ylrz1rW9Nsj4UTp8+PS+//HLWrl27ybV/krz3ve/Nb/zGb+Taa6/NJZdcktmzZ+eEE07IzTffnOOPPz677757z7ZTp07Nxz72sTz88MM5/PDDc9ppp/XaV2dnZ84///wceeSRGTduXHbaaad8+ctfzn777Tegn8NAlA3HNQ+lKVOm1KVLl/Za9tBDD/V0kzMwK1euzEknnZTly5e3LsX5A5oqpdxTa/Vg5u3QV9sM0BfXfdtn/vz5Wbp06SbPTB1qt99+e+bOndvzXNft0dfvxEDbZsN3AQAAaGZYD99lyyZMmDAsekkBAICBmTlzZmbOnNm6jBx99NE5+uijW5chlAIjz6z5d2+yrOW9LEAfFs/ZdNm084e+DgCGPcN3AQAAaEYoBQAAoBnDdwG2l2GKAADbbHiH0r4u9LbHAC4SSyn54Ac/mCuuuCJJ8vLLL2efffbJoYce2jNV8k033ZQLL7wwL730UtasWZNjjjlmk+mc58+fn09+8pPp6OjI2rVrs88+++TCCy/c5PlDGxus6aG/+MUv5qyzzsov/dIvbdP3D+b00AAAjHINrts39Ad/8AfZY4898olPfKLfba677rq8+c1vzkEHHTSo+93YxRdfnC9/+cuZOnVq/uZv/mbA3zfYVq5cmW9/+9v5wAc+0KyGgTJ8dyO77757li9fnl/84hdJkltuuSX77rtvz/q77ror5557bq6++urcf//9eeCBB3LwwQf3ua8ZM2bkvvvuy7Jly/LZz34273//+/PQQw8NyXF88YtfzAsvvDAk7wUAAMPdddddlxUrVuzw9/nKV76S2267bcCBdO3atTukjpUrV+Zv//Zvd8i+B5tQ2ocTTjghN9xwQ5LkyiuvzOmnn96z7uKLL84FF1yQCRMmJEle9apX5eyzz97iPg877LCcc845ueyyy5Ksn375d37ndzJ16tS85S1vyZ133rnJ98ycOTPXXHNNz+s99tgjSbJq1aocddRR6ezszMSJE/Otb32r1/d96UtfypNPPplp06Zl2rRpvb43Sa655pqeKahnzpyZ2bNn59BDD80b3/jG/P3f//0mdTz//PM5/fTTM3ny5Bx00EG5+uqrt3i8AACwI11wwQV505velHe961155JFHepZ/5StfydSpU3PQQQflpJNOys9+9rN8+9vfzqJFi/LJT34ynZ2d+f73v9/ndn25//778853vjNvfOMbe41mvOiiizJp0qQccMAB+dSnPpUkmT17dv793/8973nPe3LxxRfnxz/+cY477rhMnDgxhxxySO65554k63tgP/zhD+dd73pXfuu3fisvv/xyPvaxj2Xy5Mk54IAD8qd/+qd91vKVr3wlBx54YDo7O/PBD34wSf+Z4VOf+lTuuOOOdHZ25uKLL96On/SOJ5T24bTTTstVV12VF198McuWLcuhhx7as27ZsmV529vetk37PeSQQ/Lwww/3vH7ppZdy991356tf/Wo+8pGPDHg/f/u3f5sTTzwxXV1deeCBB3LIIYf0Wv/xj388r3vd67J48eIsXrx4i/t77LHHctddd+XWW2/N7NmzN+lhveCCC3LiiSdm2bJl+fa3v53zzz8/P/3pTwdcLwAADKbvfOc7+Yd/+IesWLEiN954Y+6++/8+Lu60007L3XffnRUrVqSzszNf+cpXcvjhh2f69On5/Oc/n66urrzpTW/qc7u+LFu2LLfddlvuueeezJkzJ4899lgWLVqUVatWZdmyZVmxYkWWL1+eW265JfPmzeu5Dj/vvPPy6U9/OkcddVSWL1+eiy++OB/60Id69rtixYrccsst+Zu/+Ztceuml2XvvvbNs2bLcf//9WbBgQb73ve/1quPee+/N3Llz853vfCddXV358z//883+jD73uc/lyCOPTFdXV84777zt+GnveMP7ntJGJk+enJUrV+bKK6/MCSecMGj7rbX2en3qqacmSY444oi8+OKLeeqppwa0n8MOOyyzZs3KL37xi7z3ve/d5pD8ilNOOSWllLzhDW/IAQcckOXLl/da/41vfCPf+MY3Mnfu3CTJmjVr8thjj2XixInb9b4AALAt7rjjjrzvfe/LLrvskl122SXTp0/vWbdkyZJ85jOfyQsvvJCf/exnOfbYY/vcx0C3O/nkk7Pzzjtn5513zrHHHpu77ror3/rWt/KNb3wjb33rW5OsH1n46KOPbvK9d955Z08v6rve9a48//zzefrpp5Mk06dPz84775xk/fX29773vZ4ez5/85Cf5/ve/nze/+c09+7r11ltz6qmn5tWvfnWS9Pw9Ggil/Zg+fXo+8YlP5Pbbb88zzzzTs3zSpEm59957c8ABB2z1Pu+7774ceOCBPa9LKb3Wb/x6p512yrp165Ik69aty+rVq5MkRx11VP75n/85N954Y84888yce+65OeOMMzb73hsG4hdffHGz77vx61prFi1alDe96U2bfQ8AABgKG1+vbuiMM87IzTffnEmTJmX+/Pm5/fbbt2u7vq6Va635zGc+k1mzZm3rIWT33Xfv+brWmi9/+cv9BuPN6S8zjCSG7/bjIx/5SC688MJMmjSp1/Jzzz03F110UX7wgx8kWX/i582bt8X9LVmyJJdeemk++tGP9ix75ZOQ73znO9l1112z11579fqejo6OnnHnN9xwQ9asWZMkefzxx7P33nvnzDPPzKxZs3oNV3jFbrvt1msY7mtf+9o89NBDqbXmuuuu67Xttddem1prVq5cmUceeWSTHtDjjjsuX/7yl3teb9yTCgAAQ+md73xnrrvuurz00kt54YUX8vWvf71n3UsvvZRf+7Vfy9q1a3tN9LPx9XF/221s0aJFWb16dZ577rnceuutOfTQQ3Pcccflr/7qr3omR/3Rj37U0wO6oSOPPDJXXXVVkvW9u3vssccm1/zJ+uvtyy67rGfSo0cffXSTW+qOPfbYXH311XnuueeSpOfv/jLDxsc7nA3vntKGz/nr6OjIxz/+8U2WH3744Zk7d25OOeWUrF69OmvXru33E42FCxfmX/7lX7Ju3brsvffeWbhwYa+e0vHjx+ftb397nnvuuXz1q1/d5Ptnz56dE044ITfffHOOP/74nk9Tbr311sydOzfjx4/PHnvs0ef3zpo1K9OmTcv++++fxYsXZ86cOTnuuOOy33775a1vfWuef/75Xsf6jne8I0899VQuvfTS7Lbbbr329Yd/+Ic5++yzc+CBB2bcuHHp6OjITTfdNLAfJAAAo98QX7e/4x3vyH/7b/8tBx10UDo6OjJ16tSedZ/97Gfztre9LR0dHZkyZUrPBEYzZszImWeemS984Qu59tpr+91uY5MmTcoxxxyTJ554Iueff37233//7L///lmxYkUOOeSQ7Lzzztlll11y5ZVXbhI4//iP/zgf+MAHcuWVV2bcuHH52te+1ud7nHPOOVm5cmUOPvjg7Lzzznn1q1/dK2gn6+en+V//63/lsMMOy6677ppJkybla1/7Wr+ZobOzM6tXr87EiRMza9asYX1fadn4PsehMmXKlLp06dJeyx566KFeoW00O/roozN37txMmTKlaR0zZ87MSSedlFNOOWW79zWWzh9tzZq/6eiAy2dO7WPLIdLXs9kafqg2VpVS7qm1tv1PdYTrq23eZv5dwKjmuo+N9fU7MdC22fBdAAAAmhnew3dHsf5upB5q8+fPb10CAAAwhg27ntJWw4nZPs4bAMDY4vqPV2zv78KwCqW77rprnnnmGb/gI0ytNc8880x23XXX1qUAADAEXLfzisHIAsNq+G5HR0dWrVrV53TKDG+77rprOjo6WpcBAMAQcN3OhrY3CwyrUDp+/Pi84Q1vaF0GAACwGa7bGUzDavguAAAAY4tQCgAAQDNCKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNDKvnlAK0Mmv+3Zssu3zm1AaVAACMLXpKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoRSgEAAGjG7LsAI5yZgwGAkWyLPaWllK+WUp4qpSzvZ30ppXyplLKilHJfKeWQwS8TANgWpZSVpZQHSildpZSlresBgI0NZPju/CTHb2b9byR5fZKDk8xK8lfbXxYAMIim1Vo7a61TWhcCABvbYiittX4ryX9uZpMTk1xR17s3ybhSyn6DVSAAAACj12BMdNSR5PENXq/qXgYAtFeT3NI9hPd/bryylHJWKWVpKWXp008/3aA8AMa6IZ3oqJRyVpKzkmT//fcfyrcGGJ4Wz+l7+bTzh7YORrN31Fp/WEr51ST/VEp5uNZ6yysra62XJbksSaZMmVJbFQnA2DUYPaWrkmw4XLeje9kmaq2X1Vqn1Fqn7LXXXoPw1gDA5tRaf9j991NJrkliamYAhpXB6Cm9McmHkvxd98y762qtj2/hewCa6OvxKTBalVJ2T1JrrS90f318ki80LgsAetliKC2lXJnk6CR7llJWJbkwyfgkqbXOS3JtkmmllBVJVif57R1WLQCwNX4tyXWllJrkl5IsTHJ925IAoLcthtJa6+lbWF+TnDNoFQEMA9OfXZAs/ubg7rS/+0dhB6m1/nuSya3rAIDNGYx7SgEAAGCbCKUAAAA0I5QCAADQjFAKAABAM0IpAAAAzQilAAAANCOUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNCKUAAAA0I5QCAADQjFAKAABAM0IpAAAAzYxrXQDAYJg1/+7WJQAAsA30lAIAANCMUAoAAEAzQikAAADNCKUAAAA0I5QCAADQjFAKAABAMx4JA4wsi+dk+rNP9Fq06DVnNCoGAIDtpacUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGbMvgswhK7vemKTZSd37tugEgCA4UFPKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNCKUAAAA045EwwJCbNf/uTZZdPnNqg0rWm/7sgsHf6eI5fT7+pS99bbfoB3c3/ZkAAAwVPaUAAAA0I5QCAADQjFAKAABAM+4pBcaUHXL/6A4w/dkFyeJv9l447fw2xQAA7EB6SgEAAGhGTykwbPU1S+/0Zwc2oy0AACODnlIAAACaEUoBAABoRigFAACgGaEUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoZ17oAgB1h+rMLWpcAAMAA6CkFAACgGaEUAACAZgzfBejH9V1PbLLs5M59B7wtAABbpqcUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGbMvgs0N/3ZBcnib/ax5teHvBYAAIaWnlIAAACaEUoBAABoRigFAACgGfeUAiPe9GcXtC4BAIBtNKCe0lLK8aWU5aWUh0opn+pj/d6llFtLKStKKd8rpcwe/FIBgK1VSnlVKeW+Uso/tq4FAPqyxVBaStklybwk70kyOckppZRDNtrsY0mW1loPSnJEks+VUnYb7GIBgK32O0keal0EAPRnIMN3D03yYK318SQppSxMcmKSezfYZlWSyaWUkmSPJD9O8tIg1wowplzf9USv14t+cHejShipSikdWd9m/3GS321cDgD0aSDDdzuSPL7B61Xdyzb0F0kOSvJkkgeS/E6tdd3GOyqlnFVKWVpKWfr0009vY8kAwAB9Mcn/TrJJm/wKbTMArQ3W7LvnJ1mW5HVJOpNcUkr5fzbeqNZ6Wa11Sq11yl577TVIbw0AbKyUclKSp2qt92xuO20zAK0NJJSuSrLfBq87updt6Mgkf1fX+7ckj2Z9zykA0MYRSaaXUlYmuSrJMaWUK9qWBACbGkgoXZJkYimlo5QyPsmMJDdttM33kxybJKWUX8v6QLpyEOsEALZCrfX8WmtHrXVCktOS3FZr/VDjsgBgE1uc6KjW+mIp5ewkN2d9iL2i1rr0lce+1FrnJbkoyRWllIeSvCrJZ2qtP9yBdQMAADAKDGT23dRab0xy40bL5m3w9Q+T/PrglgYADIZa6+1Jbm9cBgD0aUChFGBbzZrvMSYAAPRvsGbfBQAAgK0mlAIAANCMUAoAAEAzQikAAADNCKUAAAA0I5QCAADQjFAKAABAM0IpAAAAzQilAAAANCOUAgAA0My41gUA9Gf6swtalwAAwA6mpxQAAIBmhFIAAACaEUoBAABoxj2lwI6zeE6mP/tEr0WLXnNGo2IAABiO9JQCAADQjFAKAABAM0IpAAAAzbinFBhSnj0KAMCGhFKArXB91xNb3ggAgAEzfBcAAIBmhFIAAACaEUoBAABoRigFAACgGaEUAACAZsy+CzAKzZp/d5/LL585dYgrAQDYPD2lAAAANCOUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADRj9l2AEWz6sws2WbboNWc0qAQAYNsIpcBW6+txIx41AgDAtjB8FwAAgGaEUgAAAJoxfBdghOjr/lEAgJFOTykAAADNCKUAAAA0I5QCAADQjFAKAABAM0IpAAAAzQilAAAANCOUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNCKUAAAA0I5QCAADQjFAKAABAM0IpAAAAzQilAAAANCOUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQzrnUBALQ1a/7dmyy7fObUBpUAAGORnlIAAACa0VMKDI7Fc1pXAADACKSnFAAAgGaEUgAAAJoRSgEAAGhGKAUAAKCZAU10VEo5PsncJK9KsqDW+rk+tjk6yeeT7JzkJ7XWowaxTmCUu77ridYlwKhTStk1yZ1Z397vnuSGJOfVWmvTwgBgA1sMpaWUXZLMS3Jkkh8m+U4p5Ru11ns32GbvJH+e5Nha6w9LKXvuqIIBgAF7Kcm7aq0/L6WMz/qAOi3JbW3LAoD/ayDDdw9N8mCt9fFa65okC5OcuNE2pyW5utb6wySptf54cMsEALZWXe/n3S/HZ/2Ip6calgQAmxhIKO1I8vgGr1d1L9vQAUn2KaXcVUp5oJTy0b52VEo5q5SytJSy9Omnn962igGAASulvKqU0pX1YfT2WuvyjdZrmwFoarAmOtopSWeSY7N+WNDvlVImbrxRrfWyWuuUWuuUvfbaa5DeGgDoT611ba21M+s/UD6ylDJto/XaZgCaGkgoXZVkvw1ed3Qv29DjSW6utf68e+juPyeZPDglAgDbq9b6XNZPdHRY61oAYEMDCaVLkkwspXR0T5IwI8lNG21zQ5J3llLGlVJ+Kck7kjw8uKUCAFujlLJnKeWXu7/eLcm7kyzf/HcBwNDa4uy7tdYXSylnJ7k560PsFbXWpaWU2d3r59Va7y2l/FOSZVk/kcLlG87OCwA08bokf11KKUl2TXJlrfXrjWsCgF4G9JzSWuuNSW7caNm8jV5/PuufUwoADAO11mVZP+cDAAxbgzXREQAAAGw1oRQAAIBmhFIAAACaEUoBAABoRigFAACgGaEUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoZ17oAYJhYPGfTZdPOH/o6AAAYU/SUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNCKUAAAA0M651AQAMrunPLuhz+aLXnDHElQAAbJmeUgAAAJoRSgEAAGhGKAUAAKAZ95QC/Vs8Z9Nl087vc9Pru57YwcUAADAa6SkFAACgGaEUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJrxnFIYzbbiOaMAANCCnlIAAACaEUoBAACC0AFjAAAQ/klEQVRoRigFAACgGaEUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGbGtS4AgIYWz8n0Z5/otWjRa85oVAwAMBbpKQUAAKAZoRQAAIBmhFIAAACaEUoBAABoRigFAACgGaEUAACAZoRSAAAAmvGcUmDr9PFcSwAA2FZ6SgEAAGhGKAUAAKAZoRQAAIBmhFIAAACaEUoBAABoRigFAACgGY+EAZIk13f1/ZiXkzv3HeJKAAAYS/SUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQjlAIAANCMUAoAAEAzQikAAADNCKUAAAA0I5QCAADQjFAKAABAMwMKpaWU40spy0spD5VSPrWZ7aaWUl4upZwyeCUCANuilLJfKeVb3W3490opv9e6JgDY2BZDaSlllyTzkrwnyeQkp5RSDulju1cl+ZMk3xjsIgGAbbImycdqrROTvC3JmaWUzsY1AUAvA+kpPTTJg7XWx2uta5IsTHJiH9v9zyTXJnlqEOsDALZRrfWHtdZl3V//LMmyJPu2rQoAehs3gG06kjy+wetVSY7ecINSyr5J3pdkWpKp/e2olHJWkrOSZP/999/KUoEWru96onUJDKJZ8+/u9Xr6s87vWFFKmZD1bfRHNlqubQagqcGa6OiLSX6v1rpucxvVWi+rtU6ptU7Za6+9BumtAYDNKaXskeSaJOfWWn+y4TptMwCtDaSndFWS/TZ43dG9bENTklxVSkmSPZOcUEp5udZ63aBUCQBsk1LK+Ky/vebKWuvft64HADY2kFC6JMnEUkpHkh8lmZFk9oYb1Frf8MrXpZT5Sf5RIAWAtsr6T4svT/JQrfX/a10PAPRli8N3a60vJjk7yc1ZP0HCP9Ral5ZSZpdSZm/+uwGAho5I8uEkx5RSurr/nNC6KADY0EB6SlNrvTHJjRstm9fPtjO3vywAYHvVWu9MUlrXAQCbM1gTHQEAAMBWE0oBAABoZkDDd4ERYPGc1hUAAMBW01MKAABAM0IpAAAAzQilAAAANCOUAgAA0IyJjmCsMSESAADDiJ5SAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGZMdAQAADAW9TcB5rTzh7QMPaUAAAA0I5QCAADQjOG7MBJ51igAAKOEnlIAAACaEUoBAABoRigFAACgGaEUAACAZkx0BKPc9V1PbLLs5M59G1QCAACb0lMKAABAM0IpAAAAzRi+CzBGTH92QesSAAA2oacUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoRSgEAAGhGKAUAAKAZoRQAAIBmxrUuYDDMmn/3Jssunzm1QSUAAABsDT2lAAAANCOUAgAA0IxQCgAAQDNCKQAAAM0IpQAAADQjlAIAANDMqHgkDAA7Xl+P30o8ggsA2D56SgEAAGhGKAUAAKAZoRQAAIBmhFIAAACaEUoBAABoxuy7MIpc3/XEoG4HAAA7mp5SAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGaEUgAAAJoRSgEAAGjGI2FgOFs8p8/FHunCjjT92QW5/uIFvZYtes0ZjaoBAEY7PaUAAAA0I5QCAADQjFAKAABAM+4pBWCLpj+7oP+Vi7/Z+/W083dsMQDAqKKnFAAAgGaEUgAAAJoRSgEAAGhGKAUAAKAZoRQAAIBmhFIAAACaEUoBAABoxnNKYUdbPGfTZZ7jCAAASQYYSkspxyeZm+RVSRbUWj+30foPJ/nfSUqSl5L891rr0kGuFcak67ueaF0CMEKVUr6a5KQkT9VaJ7auBwD6ssXhu6WUXZLMS/KeJJOTnFJKOWSjzb6X5J3dDd7/m+QvB7tQAGCrzU9yfOsiAGBzBnJP6aFJHqy1Pl5rXZNkYZITN9yg1vrdWutPul/emWTfwS0TANhatdZvJfnP1nUAwOYMZPhuR5LHN3i9KsnRm9n+vydZ1NeKUspZSc5Kkv33339gFUJr7gmFrePfzIiibQagtUGd6KiUcnSSWUne2df6WutlSS5LkilTptTBfG8AYOtpm4HtNWv+3X0uv3zm1CGuhJFqIKF0VZL9Nnjd0b2sl1LK5CSXJ3lPrfWZwSkPAACA0WwgoXRJkomllI4kP0oyI8nsDTcopeyf5O+TfLjW+r1BrxIAAGAH6KunVy/v0NpiKK21vlhKOTvJzVk/MdIVtdalpZTZ3evnJbkgyWuTXFpKSZKXa61TdlzZAMCWlFKuzPp5IPYspaxKcmGt9fK2VQH01t/w374MNCwKmiPLgO4prbXemOTGjZbN2+DrM5OcObilAQDbo9Z6eusaAEYioXZoDepERwAAADuKsDg6CaUAAMCotzXDhBlaQikAADBiCZsjn1AKAAAMOkNtGSihFAZRn//5vr5BIQAAMEIIpQAAwJAw1Ja+CKUAADCG9RcUDbVlqAil0MLiOZsum3b+0NcBALAV9HSyIwilAAAAWzDiJ27qq1NkmNipdQEAAACMXUIpAAAAzRi+C8PF4jm5vuuJ1lXAVuvr9/bkzn0bVAIAw8OIH+o7xPSUAgAA0IyeUgAAGCPMnju4/DwHh55SAAAAmhFKAQAAaEYoBQAAoBn3lMI2Gug9BP3NqGt2UgAAEEoBAIA+mMRnBFs8p3UFW8XwXQAAAJrRUwob2o5PlaY/u2AQCwEAgLFBTykAAADNCKUAAAA0I5QCAADQjHtKGbH6mhHu8plTt3vbAVk8J9Of7ftRLwAAMCRG2Cy7/RFKAQBghBv0D+BhCAmlAAAAO5gPDvrnnlIAAACaEUoBAABoxvBdAAAYhfoaLgrDkVDKqDKS/vO9vsvsvQAA9Nbf9ezlr//mEFcydAzfBQAAoBk9pYx6059d0PeKxaP30yYAABgphFIAABghRtKtSmy7PjtVXr/v0BcyRIRSAAAYhgRQxgr3lAIAANCMnlJGlX7vHwUAgOFo8ZxeL6c/O/ae0CCUAjDoNn7k0aIf9D8E7fKZU3d0OQDAMCaUAgAw9DbqHUqSTDt/6OuAIdLniD5Pg0gilAIAwJDqawIjo0YYy4RSAAB2nL56RLdmW72nMOoJpQAAAINooJNvbjwHw1gllAIAMHyN4N7TrXnOqGeSjlye/rD9hFIAAEaWrRkSPEICLIxlQilswBAKABhlRnBPK+3o/RxaQinDX7+fhv76kJYBbLu+GvdFrzmjQSUAmyHAjkkCaHtCKQAAg2NrhtW2tL11Cq8jlgA6PAmlAADQnwEG2Osv/tgmy6bHqJChIGiOfEIpAADsINsTmMZKoBUqEUoZsfwHBgCMZltzrdNXgB3Ka6WBvv9YCdpsHaEUAABGuNYf1g/0/VvXOdr09eSIkzv3bVDJ9tmpdQEAAACMXUIpAAAAzRi+CwAAMMz1NVR3oNsN9yG9QikAAGyFgYYDYGAM3wUAAKAZPaUAAGy9xXNaV7DNRuLwRhjNhFJGBMNkAIAWXIPAjieUAgAw5gmf0I5QCgDAiNJfgOxrCK6wCcOfUAoAwJDamlA5GPsFhjehFAAYdH1OJDOtQSGMeIImbL8d9UHQYPFIGAAAAJrRU8qY5ZNXABhc2lZgW+gpBQAAoBk9pQAAAGNQn/f/N7jPdFSE0unPLuhj6dQhrwMAYKzYEUN1Df+FsWlAw3dLKceXUpaXUh4qpXyqj/WllPKlUsqKUsp9pZRDBr9UAGBrbKn9BoDhYIuhtJSyS5J5Sd6TZHKSU/oInb+R5PVJDk4yK8lfDXKdAMBWGGD7DQDNDaSn9NAkD9ZaH6+1rkmyMMmJG21zYpIr6nr3JhlXStlvkGsFAAZuIO03ADQ3kHtKO5I8vsHrVUmOHsA2Gy9LKeWsJGd1v3y+lPLI1hS7GXsm+XGvJb/754O062Fr02MeG8bicY/FY04c9xiw/v/pr/72oB3z6wdhH6PJQNpvbfPgGkP/fnsZi8c9Fo85GZvHPRaPOcmnh7RtHtKJjmqtlyW5bLD3W0pZWmudMtj7Hc7G4jEnY/O4x+IxJ467dR1DaSwe83CibR48Y/GYk7F53GPxmJOxedxj8ZiToT/ugQzfXZVkw6G4Hd3LtnYbAGDoaJsBGBEGEkqXJJlYSukopYxPMiPJTRttc2OSDyZJ9yQK62qtjwcAaGUg7TcANLfF4bu11hdLKWcnuTnrQ+wVtdalpZTZ3evnJbk2ybRSyookq5P89g6suS+DPuxoBBiLx5yMzeMei8ecOO6xZCwe8w7XX/s9hCWMxfM6Fo85GZvHPRaPORmbxz0WjzkZ4uMutdahfD8AAADoMZDhuwAAALBDCKUAAAA0MyJDaSnlC6WUh7r/3FBK2bOf7Y4vpSzv3u5TQ13nYCqlvL+U8mApZV0ppd/pmUspK0spD5RSukopQ3nv0A6xFcc9ms71r5RSbuk+j98opbymn+1Gxbne0rkr632plLKilHJf92RqI9oAjvnoUspPus9tVynlghZ1DqZSyldLKU+VUpb3s37UneexRtusbe5ju9F0rsdM2zwW2+VE29zP+qE717XWEfcnyTFJxnV//SdJvtjHNrskWZn10+GPT7I0ySGta9+OYz4wyVuS3J5kyma2W5lkz9b1DuVxj8Jz/WdJfrf76/OSfGm0nuuBnLskv5nk+iQlySFJ7m9d9xAc89FJ/rF1rYN83Ed1n7/l/awfVed5LP7RNmubR/m5HhNt81hsl7fiuLXNO/Bcj8ie0lrrbbXWl7tf3plk3z42OzTJg7XWx2uta5IsTHLiUNU42GqtD9VaH2ldx1Ab4HGPqnOd9bV/rfvrKzKyj2VLBnLuTsz6WUNrrfXeJONKKfttvKMRZLT9vg5IrfVbSf5zM5uMtvM85mibxw5t86hum8diu5yMvt/XARlObfOIDKUbOSvJoj6WdyTZ8Fmpq7qXjXY1ySvDS/5n62KGyGg713vVWp9Oku6/f7Wf7UbDuR7IuRtt53egx/OO7mFEt5VSOoemtKZG23ke67TNvY2G/6+31mg712OlbR6L7XKibe7PkJ3rLT6ntJVSyjeT7N3Hqt+vtV7fvc3vJ3k56z+xGvEGcswD8I5a6w9LKb+a5J9KKQ/XWm8ZvCoH3yAd94iyuWPeit2MuHPNgN2TZL9a6wullOOSXFdKeWOtdV3rwhjbtM29aJtHGW0zW6Bt3oGGbSittf765taXUs5I8t4kx9TuQc8bWZX148Jf0dG9bNja0jEPcB8/7P77qVLKNUmmJhnW/xkOwnGPqnNdSnm6lLJXrfXpUspeSZ7qZx8j7lz3YSDn7pVt7trMNiPJFo+51vqzDb6+uZSyOusvlJ4ckgrbGG3neVTSNm/zPkbc/9fa5t7GUNs8FtvlRNvcnyE71yNy+G4p5fgkv5fkvbXWF/rZbEmSiaWUjlLK+CQzktw0VDW2UErZvZTyS698neT4JCvaVjUkRtu5vjHJh7q//lD6OJZRdK4Hcu5uTPLBJOme9W1drfXxjFxbPObuC55Xvn5bkj3SzwXQKDLazvOYo23u2yj6/3prjbZzPVba5rHYLifa5v4M3bneUTMo7cg/Sf4t68c3d3X/mde9/HVJbtxguxOSPJjkoawfbtK89u045vdl/ScTLyX5UZKbNz7mJG9MsizJ/Un+NckfJSmta9/Rxz0Kz/Vrk3wzyQPdf//KaD7XfZ27JLOTzO7+uiT586xv2LuymRkuR8qfARzzx5Ms7/5zb5KjW9c8CMd8ZZL/SLKm+9/0rNF+nsfaH22ztnmUn+sx0zaPxXZ5gMetbd6B57p0vyEAAAAMuRE5fBcAAIDRQSgFAACgGaEUAACAZoRSAAAAmhFKAQAAaEYoBQAAoBmhFAAAgGb+f2aKV+adw3ywAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1152x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.subplot(1,2,1)\n",
"plt.hist(costheta_Ds_b,alpha=0.7,density=True,bins=70, label=\"MC Ds tuple\");\n",
"plt.hist(costheta_Dplus_b,alpha=0.5,density=True,bins=70, label=\"MC Dplus tuple\");\n",
"plt.legend()\n",
"plt.subplot(1,2,2)\n",
"plt.hist(costheta_MC_b,alpha=0.7,density=True,bins=70, label=\"MC Ds+Dplus tuple\");\n",
"plt.hist(costheta_data_b,alpha=0.5,density=True,bins=70, label=\"data before cut\");\n",
"plt.legend()\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(16,8)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Cos theta_l after cut on phi mass"
]
},
{
"cell_type": "code",
"execution_count": 115,
"metadata": {},
"outputs": [],
"source": [
"costheta_Dplus=get_costheta_list(MC_Dplus_tuple_dict_presel_9, l_index=l_index, mother_index=1)\n",
"costheta_Ds=get_costheta_list(MC_Ds_tuple_dict_presel_9, l_index=l_index, mother_index=0)\n",
"costheta_MC=np.concatenate((costheta_Dplus,costheta_Ds),axis=0)\n",
"costheta_data=get_costheta_list(data_tuple_dict_presel_4, l_index=l_index, mother_index=0)"
]
},
{
"cell_type": "code",
"execution_count": 116,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA6UAAAHVCAYAAAAJnF2uAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzt3X+YXVV9L/73QsKPQr9Xq2lRhhj1PsrvjpgIomgi7QMChmopiJZLaigNxWuhV2/F3oLS3mIfc8VSxJhebKhYfhS8kGq8gBgKVCiEEEIISOs1QrAKpaBQhASyvn/MZDqZH5mTzMzZc2Zer+eZJ+fsvfben7Oz96z5nLX2WqXWGgAAAGjCTk0HAAAAwNQlKQUAAKAxklIAAAAaIykFAACgMZJSAAAAGiMpBQAAoDGSUgAAABojKQUAAKAxklIAAAAas3NTB37Vq15VZ86c2dThAZhk7rnnnn+ttU5vOo5Opm4GYCy1Wjc3lpTOnDkzK1eubOrwAEwypZQfNB1Dp1M3AzCWWq2bdd8FAACgMZJSAAAAGiMpBQAAoDGNPVMK0KRNmzZlw4YNef7555sOhe202267paurK9OmTWs6lCnBvUJ/7j9gPEhKgSlpw4YN+fmf//nMnDkzpZSmw6FFtdY8+eST2bBhQ173utc1Hc6U4F5hC/cfMF5G7L5bSvlyKeXxUsraEcrNLqW8WEo5YezCAxgfzz//fF75ylf6I7vDlFLyyle+UqtdG7lX2ML9B4yXVp4pXZrk6G0VKKW8LMmfJblxDGICaAt/ZHcm/2/t55yzhWsBGA8jJqW11luT/NsIxf5rkmuTPD4WQQEAADA1jPqZ0lLK3knel2RuktkjlD09yelJMmPGjNEeGmDMLFh695ju79L52/x1mKSnxeFDH/pQLr/88iTJiy++mFe/+tU59NBD8/Wvfz1J8s1vfjPnnXdeXnjhhWzatCnvfve7c/HFF2+1n6VLl+bjH/94urq68tJLL+XVr351zjvvvBx++OEtxbp69er88Ic/zDHHHLOdn/I/zJkzJ4sWLcqsWbN2eB90hql8rwznZS97WQ466KCUUvpiPfvss7PTTtv+7n/PPffMs88+O6pj33LLLdlll11G9RlmzpyZlStX5lWvetWoYgHYUWMxJcznk/xBrXXzSAVrrUtqrbNqrbOmT58+BocG6Fx77LFH1q5dm5/97GdJkptuuil777133/o777wzZ511Vq6++urcd999uf/++3PAAQcMua+TTjop9957b9asWZNPf/rT+Y3f+I08+OCDLcWxevXqLF++fPQfCMbJRLlX5s+fn1tuuWXQ8t133z2rV6/Ovffem5tvvjnf/va38+lPf3r7P+gOuOWWW/Kd73ynLccCGC9jkZTOSnJlKWV9khOSXFJK+bUx2C/ApHfMMcfkG9/4RpLkiiuuyMknn9y37sILL8y5556bmTNnJulpjTnjjDNG3Odhhx2WM888M0uWLOnbz/77759f/uVfzoknnrhV2Y0bN+bcc8/NVVddle7u7lx11VX51Kc+lUWLFvWVOfDAA7N+/fqsX78+++67b0455ZQccMABOfbYY/Pv//7vg46/bNmyHHLIITnwwAMzb968/PSnP93u8wIDNX2vtOoVr3hFlixZkosvvji11ixdujTHH398jjzyyLzhDW/IJz/5yUHb3HLLLTnuuOP63n/kIx/J0qVLkyQf//jH+2I6++yzt9pu/fr1Wbx4cS688MJ0d3fntttuy/z583PNNdf0ldlzzz37jvHOd74z733ve/PGN74x8+fPz+bNg9sT/vIv/zIHH3xw9t9///zWb/1WNm3atEPnAWB7jDoprbW+rtY6s9Y6M8k1SX631nrdqCMDmAI+8IEP5Morr8zzzz+fNWvW5NBDD+1bt2bNmrzlLW/Zof0ecsgheeihh5Ikn/3sZ3Pvvffmvvvuy6WXXrpVuV122SXnn39+TjrppKxevTonnXTSNvf73e9+N2eeeWYeeOCBvOY1r8mf//mfb7X+xz/+cS644ILcdtttWbt2bd7+9rfnz/7sz3boM0B/Td8r22PvvffOtGnT8vjjPUNt3HXXXbn22muzbt26/N3f/V3LLZuPP/54li9fngceeCD33Xdfzj///K3Wz5w5MwsXLszZZ5+d1atX54gjjtjm/u66665cfPHFeeihh/Iv//IvufLKK7daf9999+W6667LPffck3Xr1mX33XfvS44BxlMrU8JckeSOJG8qpWwopSwopSwspSwc//AAJreDDz4469evzxVXXDGqZzoHqrVudYxTTjklf/3Xfz3q/e6zzz457LDDkiQnn3xybr/99q3W33bbbfmnf/qnvP3tb093d3cuu+yyPPbYY6M+LjR1r9xwww3p7u5Od3d3li1bltNOOy3d3d1bJcUj7fdXf/VX8/KXvzy77rpr3ve+9w26b4bzC7/wC5k2bVoWLFiQa665Jrvssst2frqtvfWtb81rX/va7LTTTjnppJMGxXHTTTfl3nvvzezZs9Pd3Z2bb745jzzyyKiOCdCKEQc6qrWePFKZfmXnjyoagClo3rx5+djHPpZbbrklTz75ZN/ygw46KKtWrcq+++673fu89957s99++yVJvvGNb+TWW2/N17/+9fzpn/5p1q5dm513Hv7X/0477bRVt77+cxIOnA5i4Ptaa97znvfkK1/5ynbHDCNp4l456qijctRRRyXpeaZ0/vz5mTNnzjb3+cMf/jAvvfRSfvEXfzHJyPfNcPfczjvvnH/8x3/MzTffnGuvvTZf+MIXsmLFim0eu/++Nm/enI0bNw573KHu3wULFuSP//iPt3kMgLE2Fs+UAjAKH/7wh3PeeefloIMO2mr5WWedlfPPPz8/+MEPkvT8gbl48eIR93fXXXflkksuyW//9m9n8+bN+eEPf5i5c+fmM5/5TH7605/m6aef3qr87rvvnueee67vfVdXV1atWpWkZxCk73//+33rHnnkkdx1111Jkquuuipvf/vbt9rXEUcckRUrVvTF/Pzzz+ef//mfWz0VsE1N3yutePrpp7Nw4cKceeaZfUnfTTfdlJ/85CfZuHFjrrvuukH3TVdXVx544IG88MILeeaZZ3LzzTcnSZ599tk8++yzOeaYY/K5z32u777sb6j795577knSk2T3fyb0rrvuyiOPPJJaa/72b/92UBy/+qu/mquvvjr/9m89MwH+9Kc/zaOPPrrd5wBge416ShiAyaCVaSnGS1dXVz760Y8OWn744Ydn0aJFOeGEE7Jx48a89NJLOfLII4fcx1VXXZV/+Id/yObNm7PXXnvlqquuyn777ZdNmzblpJNOyr//+7/npZdeyu/+7u8OmvZh7ty5+dM//dMcfPDB+cM//MOceOKJueyyy3LAAQfksMMOyxvf+Ma+sm9605vyF3/xF1m1alVe+9rX5nOf+9xW+9prr72yZMmSzJs3L0lPcvAnf/In+c//+T+P9jQxQUzle2U4P/vZz9Ld3d03JcwHP/jB/P7v/37f+re+9a359V//9Xz/+9/PiSeeOGj6lte//vU5/vjjs+++++ZNb3pT3vzmNyfpSQrnzZuXF198MS+99NKg+y1J3vve9+b9739/rr322lx88cVZuHBhjjnmmNxwww05+uijs8cee/SVnT17dj7ykY/koYceyuGHH54PfOADW+2ru7s755xzTo444ojsvPPO2WmnnfLFL34x++yzT0vnAWBHlf7PPLTTrFmz6sqVKxs5NsCDDz7Y12WP1qxfvz7HHXdc1q5d23QoQ/7/lVLuqbWaKHUUhqqb3Sujs3Tp0qxcuXLQnKntdsstt2TRokV987qOhmsCaFWrdbPuuwAAADRG910AWjJz5swJ0UoKnWTL4EhNmzNnzogDNAE0RVIKk9yCpXcPWtbkM2EAAEwQKy4Yevncc9oahu67AAAANEZSCgAAQGMkpQAAADTGM6UAyfDPVOyoFp7FKKXkQx/6UC6//PIkyYsvvphXv/rVOfTQQ/umbfjmN7+Z8847Ly+88EI2bdqUd7/73YOmlli6dGk+/vGPp6urKy+99FJe/epX57zzzhs0F+JAYzVVxec///mcfvrp+bmf+7kd2n4sp6qgDRq4V/r71Kc+lT333DMf+9jHhi1z3XXX5Y1vfGP233//Md3vQBdeeGG++MUvZvbs2fnqV7/a8nZjbf369fnOd76TD37wg43FADAaWkoBGrLHHntk7dq1+dnPfpYkuemmm7L33nv3rb/zzjtz1lln5eqrr859992X+++/PwcccMCQ+zrppJNy7733Zs2aNfn0pz+d3/iN38iDDz7Yls/x+c9/Ps8991xbjgWtuO6667Ju3bpxP86XvvSlfPvb3245IX3ppZfGJY7169fnb/7mb8Zl3wDtICkFaNAxxxyTb3zjG0mSK664IieffHLfugsvvDDnnntuZs6cmSR52cteljPOOGPEfR522GE588wzs2TJkiQ9U0H83u/9XmbPnp03velNuf322wdtM3/+/FxzzTV97/fcc88kyYYNG/LOd74z3d3dOfDAA3Prrbdutd1FF12UH/7wh5k7d27mzp271bZJcs011/RNhzF//vwsXLgwhx56aF7/+tfna1/72qA4nn322Zx88sk5+OCDs//+++fqq68e8fMy+Z177rl5wxvekHe961357ne/27f8S1/6UmbPnp39998/xx13XJ555pl85zvfybJly/Lxj3883d3d+d73vjdkuaHcd999ecc73pHXv/71W/UgOP/883PQQQdl3333zSc+8YkkycKFC/P//t//y3ve855ceOGF+dd//dccddRROfDAA3PIIYfknnvuSdLTAnvKKafkXe96V/7Lf/kvefHFF/ORj3wkBx98cPbdd9/8+Z//+ZCxfOlLX8p+++2X7u7ufOhDH0oy/H36iU98Irfddlu6u7tz4YUXjuJMAzRDUgrQoA984AO58sor8/zzz2fNmjU59NBD+9atWbMmb3nLW3Zov4ccckgeeuihvvcvvPBC7r777nz5y1/Ohz/84Zb38zd/8zc59thjs3r16tx///055JBDtlr/0Y9+NK95zWuyYsWKrFixYsT9PfLII7nzzjtz8803Z+HChYNaWM8999wce+yxWbNmTb7zne/knHPOyU9/+tOW42XyueOOO/J//s//ybp167J8+fLcffd/THP1gQ98IHfffXfWrVuX7u7ufOlLX8rhhx+eefPm5bOf/WxWr16dN7zhDUOWG8qaNWvy7W9/O/fcc08uuOCCPPLII1m2bFk2bNiQNWvWZN26dVm7dm1uuummLF68uO/aP/vss/PJT34y73znO7N27dpceOGF+c3f/M2+/a5bty433XRTvvrVr+aSSy7JXnvtlTVr1uS+++7LZZddlocffnirOFatWpVFixbljjvuyOrVq/OFL3xhm+foM5/5TI444oisXr06Z5999ijONkAzPFMK0KCDDz4469evzxVXXJFjjjlmzPZba93q/Yknnpgkefvb357nn38+jz/+eEv7Oeyww7JgwYL87Gc/y3vf+94dTpK3OOGEE1JKyete97rsu+++Wbt27Vbrb7zxxtx4441ZtGhRkmTTpk155JFHcuCBB47quHSu2267Le973/uy6667Ztddd828efP61t111135oz/6ozz33HN55plncuSRRw65j1bLHX/88dlll12yyy675Mgjj8ydd96ZW2+9NTfeeGPe/OY3J+lpzf/+978/aNvbb7+9rxX1Xe96V5599tk88cQTSZJ58+Zll112SdJzjT/88MN9LZ4/+clP8r3vfS9vfOMb+/Z1880358QTT8zLX/7yJOn7F2CykpQCNGzevHn52Mc+lltuuSVPPvlk3/KDDjooq1atyr777rvd+7z33nuz33779b0vpWy1fuD7nXbaKZs3b06SbN68ORs3bkySvPOd78zf//3fZ/ny5TnttNNy1lln5dRTT93msfsnxM8///w2jzvwfa01y5Ytyxve8IZtHoOpY+A10t+pp56aG264IQcddFCWLl2aW265ZVTlhro+a635oz/6oyxYsGBHP0L22GOPvte11nzxi18cNjHeluHuU4BOp/suQMM+/OEP57zzzstBBx201fKzzjor559/fn7wgx8k6fkjdPHixSPu76677soll1yS3/7t3+5btqVV5o477shuu+2W6dOnb7VNV1dX3zNw3/jGN7Jp06YkyaOPPpq99torp512WhYsWLBV18ktdt9996264b7yla/Mgw8+mFprrrvuuq3KXnvttam1Zv369fnud787qAX0qKOOyhe/+MW+9wNbUpl63vGOd+S6667LCy+8kOeeey5/93d/17fuhRdeyC/90i/lpZde2mqgn4HX5HDlBlq2bFk2btyYp59+OjfffHMOPfTQHHXUUfmrv/qrvgHJfvzjH/e1gPZ3xBFH5Morr0zS07q75557DrrPkp5rfMmSJX2DHn3/+98f1I39yCOPzNVXX52nn346Sfr+He4+Hfh5ATqNllKAZLunpRhLXV1d+ehHPzpo+eGHH55FixblhBNOyMaNG/PSSy8N27py1VVX5R/+4R+yefPm7LXXXrnqqqu2aimdNm1a3vrWt+bpp5/Ol7/85UHbL1y4MMccc0xuuOGGHH300X0tOzfffHMWLVqUadOmZc899xxy2wULFmTu3LmZMWNGVqxYkQsuuCBHHXVU9tlnn7z5zW/Os88+u9Vnfdvb3pbHH388l1xySXbfffet9vXHf/zHOeOMM7Lffvtl5513TldXV775zW+2diJpjzbfK29729vya7/2a9l///3T1dWV2bNn96379Kc/nbe85S3p6urKrFmz+gYwOumkk3Laaaflc5/7XK699tphyw100EEH5d3vfncee+yxnHPOOZkxY0ZmzJiRdevW5ZBDDskuu+ySXXfdNVdcccWghPN//s//mQ9+8IO54oorsvPOO+crX/nKkMc488wzs379+hxwwAHZZZdd8vKXv3yrRDvpeSb8v/23/5bDDjssu+22Ww466KB85StfGfY+7e7uzsaNG3PggQdmwYIFnisFOk4Z+NxRu8yaNauuXLmykWPDVLJg6eCWrUvnzx6i5NTy4IMPbpW0TWZz5szJokWLMmvWrEbjmD9/fo477riccMIJo97XUP9/pZR7aq3NfsgON1TdPJXuFVrjmoBJZLi5p8foC8hW62bddwEAAGiM7rsAk9xwg7q029KlS5sOAQCYgLSUAlNWU48vMDr+39rPOWcL1wIwHiSlwJS022675cknn/QHVoeptebJJ5/Mbrvt1nQoU4Z7hS3cf8B40X0XmJK6urqyYcOGIad2YGLbbbfd0tXV1XQYU4Z7hf7cf8B4kJQCU9K0adPyute9rukwYMJzrwAw3nTfBQAAoDGSUgAAABojKQUAAKAxklIAAAAaIykFAACgMZJSAAAAGmNKGGC7LFh696Bll86f3UAkAABMBlpKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAaDDlFL2KaXcWkpZW0p5uJTyB0OUmVNK+UkpZXXvz7lNxAoAI9m56QAAgO22KclHaq1rSik/n2RVKeWGWuvqAeVuq7Ue10B8ANCyEVtKSylfLqU8XkpZO8z6U0op9/d+W3tPKWXW2IcJAGxRa/1RrXVN7+tnkqxJsnezUQHAjmml++7SJEdvY/3DSd5Raz0wyf9I8r/HIC4AoAWllJlJZie5fYjVb+v90vjbpZTuYbY/vZSyspSy8oknnhjHSAFgaCMmpbXWW5P82zbW/2Ot9Se9b2+Pb2oBoC1KKXsmuSbJWf3q4i3uSbJP75fGf5bkulLKoHq/1rqk1jqr1jpr+vTp4x80AAww1gMd/U6SZcOt9G0sAIyNUsq0JNcmuaLW+rWB62utz9Ran+t9fUOSjUn2am+UADCyMUtKSylzkixI8t+HK+PbWAAYvVJKSXJpkgdrrf9rmDLT+71+S5I9kzzenggBoHVjMvpuKeXg9FSO76m1PjkW+wQAhvX2JKckub+UsmXE3U8mmZEktdbFSU4upZzeu25jkg/WWl9se6QAMIJRJ6WllBlJvpbklFrrw6MPCQDYllrr7UnKCGUuSnJReyICgB03YlJaSrkiyZwkryqlbEhyXpJpSd83secmeWWSS3p6E+XFWqtpYQAAABjRiElprfXkEdafluS0MYsIAACAKWOsR98FAACAlklKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAQAAaIykFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxkhKAQAAaIykFAAAgMbs3HQAwMSwYOndg5ZdOn92A5EAADCVaCkFAACgMZJSAAAAGiMpBQAAoDGSUgAAABojKQUAAKAxRt8FJoYVFwxeNvec9scBAEBbaSkFAACgMZJSAAAAGiMpBQAAoDGSUgAAABojKQUAAKAxklIAAAAaIykFAACgMZJSAAAAGiMpBQAAoDE7Nx0AwHZbccHgZXPPaX8cAACMmqQUmNIWLL170LJL589uIBIAgKlJ910AAAAaIykFAACgMZJSAAAAGiMpBQAAoDEGOgLab6jRcxs276nL/uPNim/1/GtEXwCAcaelFAAAgMZISgEAAGiM7rswSQw13yYAAEx0WkoBAABojKQUAACAxkhKAQAAaIykFAAAgMZISgEAAGjMiElpKeXLpZTHSylrh1lfSikXlVLWlVLuLaUcMvZhAgAAMBm1MiXM0iQXJ/nrYda/P8lrkxyQ5M1J/irJL49FcMDENO+py7ZesOJbydxzmglme624YKu38556rKFAAABIWmgprbXemuTftlHk2CSX1x6rkuxcStlnrAIEAABg8mqlpXQkXUke7fd+wxDLkiSllNOTnJ4kM2bMGINDAxPVgqV3D7n80vmzJ/WxhzNUTE3GAwAwUbR1oKNa65Ja66xa66zp06e389AAAABMQGORlG5I0r+7blfvMgAAANimsUhKlyf5UJL0jry7udY6qOsuAAAADDTiM6WllCuSzEnyqlLKhiTnJZmWJLXWxUmuTTK3lLIuycYkvzVu0QJTy4CRcgEAmHxGTEprrSePsL4mOXPMIgIAAGDKaOtARwAAANCfpBQAOkwpZZ9Syq2llLWllIdLKX8wRJlSSrmolLKulHJv77gPADDhjMU8pUCbDTcPZ1OuX/1Ylv1gYsXUDuYepUGbknyk1rqmlPLzSVaVUm6ota7uV+b9SV6b5IAkb07yV0l+uf2hAsC2aSkFgA5Ta/1RrXVN7+tnkqxJsveAYscmubz2WJVk51LKPgGACUZSCgAdrJQyM8nsJLcPWNWVpP8UbRt6lw3c/vRSyspSysonnnhivMIEgGFJSgGgQ5VS9kxyTZKzaq0/2ZF91FqX1Fpn1VpnTZ8+fWwDBIAWSEoBoAOVUqalZ67wK2qtXxuiyIYk/bvrdvUuA4AJRVIKAB2mlFKSXJrkwVrr/xqm2PIkH+otf0iSzbXWR4cpCwCNMfouAHSetyc5Jcn9pZQtI+5+MsmMJKm1Lk5PK+rcUsq6JBuT/FYTgQLASCSlANBhaq23JykjlKlJzmxPRACw4ySl0AatzitqjksAAKYaz5QCAADQGEkpAAAAjdF9FxgT8566bMjly15xapsjAQCgk2gpBQAAoDGSUgAAABojKQUAAKAxklIAAAAaIykFAACgMUbfhUlkqBFwp8zotysu6Hs576nHkozBZ++3zz5zzxndPgEA2IqWUgAAABojKQUAAKAxklIAAAAaIykFAACgMQY6AtgeQwyolEyhAaUAAMaYllIAAAAaIykFAACgMbrvwhS0YOndTYewQ65f/digZcd3793y9guW3r1Vl9vxOM72nNuhyl46f3bL2wMATAZaSgEAAGiMllJg0pr31GUT7vgGRAIA2JqWUgAAABojKQUAAKAxklIAAAAaIykFAACgMQY6goYMOQjPim8NXjb3nLYc2wA8AAA0QUspAAAAjdFSCrTVgqV3Z95Tj2217PjuvRuKBgCApmkpBQAAoDGSUgAAABojKQUAAKAxnimFsbTigsHLxmH03PEw7GjAo4x/yP0CAEAvLaUAAAA0RlIKAABAYySlAAAANEZSCgAAQGMkpQAAADRGUgoAAEBjJKUAAAA0RlIKAABAYySlAAAANEZSCgAAQGMkpQAAADSmpaS0lHJ0KWVtKeXBUsonhli/Vynl5lLKulLKw6WUhWMfKgAAAJPNziMVKKXsmmRxkiOS/CjJHaWUG2utq/oV+0iSlbXWPyilTE/yT6WUy2qtPxuXqGEcLVh696Bll86f3UAkAAAw+bXSUnpokgdqrY/WWjcluSrJsQPKbEjy86WUkmTPJP+a5IUxjRQAAIBJp5WktCvJo/3eb+hd1t9fJtk/yQ+T3J/k92qtmwfuqJRyeillZSll5RNPPLGDIQMAADBZjNVAR+ckWZPkNUm6k1xcSvn/BhaqtS6ptc6qtc6aPn36GB0aAACATtVKUrohyT793nf1LuvviCR/W3v8c5Lvp6flFAAAAIbVSlJ6V5IDSyldpZRpSU5K8s0BZb6X5MgkKaX8UnoS0vVjGCcAAACT0Iij79Zany+lnJHkhvQksZfXWldumfal1ro4yflJLi+lPJjkZUn+qNb6o3GMGwAAgElgxKQ0SWqty5MsH7Bscb/XP0ryK2MbGgAAAJNdS0kpwHi6fvVjg5Yd3713A5EAANBuYzX6LgAAAGw3SSkAAACNkZQCAADQGEkpAAAAjTHQEXSgeU9d1r6DrbhgwLEHD0oEAAA7SkspAAAAjZGUAgAA0BhJKQAAAI3xTCnQ0a5f7RlXAIBOJimFSa6tgyIBAMB20n0XAACAxkhKAQAAaIykFAA6TCnly6WUx0stsvu/AAAZLUlEQVQpa4dZP6eU8pNSyuren3PbHSMAtMozpQDQeZYmuTjJX2+jzG211uPaEw4A7DgtpQDQYWqttyb5t6bjAICxoKUUoI2GGg152StObSASpoC39XbvfTzJ79daVzcdEAAMRVIKY2ioOTOX/eDuBiIBprh7kuxTa32ulHJUkutKKa+vtW4eWLCUcnqS05NkxowZbQ4TAHTfBYBJp9b6TK31ud7XNyTZmGSvYcouqbXOqrXOmj59ejvDBIAkklIAmHRKKdP7vX5Lkj3T040XACYc3XcBoMOUUq5IMifJq0opG5Kcl2RaktRaFyc5ubdbbtLTSvrBWuuLTcQKACORlMJEt+KCpiMAJpha68kjrL8oyUVtCgcARkX3XQAAABojKQUAAKAxklIAAAAa45lSYEIaas7XqWDB0sHz2l46f3YDkQAAtIeWUgAAABojKQUAAKAxklIAAAAaIykFAACgMZJSAAAAGiMpBQAAoDGSUgAAABojKQUAAKAxOzcdADBxXb/6saZDYDssWHr3kMsvnT+7zZEAALROSykAAACNkZQCAADQGEkpAAAAjZGUAgAA0BhJKQAAAI2RlAIAANAYSSkAAACNkZQCAADQGEkpAAAAjZGUAgAA0BhJKQAAAI2RlAIAANAYSSkAAACNkZQCAADQGEkpAAAAjZGUAgAA0JidWylUSjk6yaIkL0tyWa31M0OUmZPks0l2SfKTWus7xzBOGLUFS+8etOzS+bMbiGR4169+rOkQGIb/GwCgo624oOkIhjViUlpK2TXJ4iRHJPlRkjtKKTfWWlf1K7NXki8kObLW+qNSyqvGK2AYK/OeuixZ8a2tF849p5lgAABgimql++6hSR6otT5aa92U5Kokxw4o84EkV9daf5QktdZ/HdswAQAAmIxaSUq7kjza7/2G3mX97Zvk1aWUO0sp95dSfnusAgQAAGDyaumZ0hbslOTgJEcm2T3JnaWUO2qta/sXKqWcnuT0JJkxY8YYHRoAAIBO1UpL6YYk+/R739W7rL9Hk9xQa/333q67f5+eJHUrtdYltdZZtdZZ06dP39GYAQAAmCRaSUrvSnJgKaWrlDItyUlJvjmgzDeSvKOUsnMp5eeSvC3JQ2MbKgAAAJPNiN13a63Pl1LOSHJDepLYy2utK0spC3vXL661riql/N8ka5JMS3Jp/9F5oZMNOUovAAAwJlp6prTWujzJ8gHLFg94/9n0zFMKU4J5KwEAYPRa6b4LAAAA40JSCgAAQGMkpQAAADRGUgoAAEBjJKUAAAA0RlIKAABAYySlAAAANEZSCgAAQGN2bjoAGGTFBYOXzT2n/XHAdpj31GXNHbz3npn31GN9i5a94tSmogEA2C5aSgEAAGiMpBQAAIDGSEoBAABojKQUAACAxhjoCMZZowPgAADABKelFAAAgMZISgEAAGiMpBQAAIDGSEoBAABojKQUAACAxhh9F6Bhw43QvOwVp7Y5EgCA9tNSCgAAQGMkpQAAADRG910AAIDJZMUFTUewXbSUAgAA0BhJKQAAAI3RfRdacP3qx5oOAQAAJiUtpQAAADRGSykAdJhSypeTHJfk8VrrgUOsL0n+PMmvJHkhyYJa66r2Rgm0w4Kldw+5/NL5s9scCew4LaUA0HmWJjl6G+vfn+S1SQ5IsiDJX7UhJgDYIVpK6VxDDXU995z2xwHQZrXWW0spM7dR5Ngkl9daa5JVpZSdSyn71FofbUuAALAdJKUAMPl0JemfgG4YYlmSpJRyepLTk2TGjBltCQ5gW4bqkqw78uSm+y4ATGG11iW11lm11lnTp09vOhwApiAtpQAw+WxIsk+SO3vfd/UuA5hQhhuoiRYN9ThbB9JSCgCTz/IkH0qSUsohSTZ7nhSAiUpLKQB0mFLKFUnmJHlVKWVDkvOSTEuSWuviJNcmmVtKWZdkY5LfaihUYIIb7fObnv9kLEhKmVy20YVh3lOPtTEQgPFTaz15hPU1yZltCgeYgEaTLJr7lHbTfRcAAIDGaCkFAIAxtj0D+GiBZKqTlAIAwATjWU2mEt13AQAAaIyWUgAAmALMCcpEJSkFAADGlYSYbZGUAkxwnisCACYzz5QCAADQGEkpAAAAjdF9l/ZYccHgZXPPaX8cAADAhCIpBQAAxsx4DGo02vEVWo3JmA3N0H0XAACAxmgpBQBgShuuFa1drWadMl3KRIvT6PSTh5ZSAAAAGqOlFPpbcUHmPfVY01FAkmTeU5cNWrbsFadu9358kwwATGSSUgAAGIWJ1q2VSWqo2SwmiZa675ZSji6lrC2lPFhK+cQ2ys0upbxYSjlh7EIEAABgshoxKS2l7JpkcZL3JDk4yQmllEOGKPeyJH+W5MaxDhIAAIDJqZWW0kOTPFBrfbTWuinJVUmOHaLcf01ybZLHxzA+AAAAJrFWktKuJI/2e7+hd1mfUsreSd6X5Ivb2lEp5fRSyspSysonnnhie2MFAABgkhmrgY4+n+QPaq2bSynDFqq1LkmyJElmzZpVx+jYAAww1Mi9yY6N3gsAMJ5aSUo3JNmn3/uu3mX9zUpyZW9C+qokx5RSXqy1XjcmUQIAADAptZKU3pXkwFJKV5IfJzkpycL+BWqtr9vyupSyNMnXJaSMxvWrt54rdNkPBg+1Pu+px3J8994jbgsAMBmYeobJasSktNb6fCnljCQ3pOcZ1MtrrStLKQt71y8e5xgBAACYpFp6prTWujzJ8gHLhkxGa63zRx8WAAAAU8FYDXQEQBv0DWC04lvNBgIAMEZamRIGAAAAxoWWUgAAJqWhBga6dP7stm0PtEZSCgAAMJGsuKDpCNpK910AAAAao6WU5kyxb4AAAIDBJKUAHej61Y81HQIAwJjQfRcAAIDGSEoBAABojKQUAACAxkhKAQAAaIyBjgAAAJpiRgpJKQAAU8eCpXc3HQIwgKQUAABgR7Xa0jn3nPGNo4NJShkXA7+FnPdUz5yKx3fv3UQ4wzLXIwAANEtSCgAA0ArPf44LSSkdYd5TlzUdAgCAZ1JhHJgSBgAAgMZoKQUAACaXAd1st4xvsuwVpzYRDSOQlAIAABPOUI9vSSonJ0kpAABAf+MxoJFBkoYlKQUAAFoy3OCTE6UFc8tAVFu66w6k9XVikpQCAMAkNNETyC3MsoCklO0y1DDol86f3UAkAADAZCApBQCAKW603Vpbbe0cj1ZaLa2dT1IKAAAMItmjXSSlAAAwhTSZbE7ERHermFZ8q7lAprCdmg4AAACAqUtLKQAA0BEmYksroycpBQCADiIxY7KRlAIA0PGGmrZuohrtSLcw2UhKaavrVz82aNnx3Xs3EAkwkHmIAXbccK2Xkk0YmaQUAIAJabjWz6nyhZluukwVklIAADpKO7vq6moL409SCgAADdMqylQmKQUAAMjoxz8xfsqO2anpAAAAAJi6tJQy9lZckHlPDf6WCJigVlyQJEPftyu+9R+v557TpoAAmtFqF1pdbWFsSUoBAJhSJJUwsUhKAaADlVKOTrIoycuSXFZr/cyA9XOSXJ/k+72LvlZrPb+tQQISYGiBpJTGDfVAOADDK6XsmmRxkiOS/CjJHaWUG2utqwYUva3WelzbA4Qd0M5pXoCJxUBHANB5Dk3yQK310VrrpiRXJTm24ZgAYIdoKZ2Megct2YoBSgAmk64kj/Z7vyHJnCHKva2UsjbJ40l+v9a6emCBUsrpSU5PkhkzZox9pNBGQ3WVXfaKUxuIBNgeklIAmJzuSbJPrfW5UspRSa4rpby+1rq5f6Fa65IkS5Jk1qxZtYE4ASa04R41M//o2JGUAkDn2ZBkn37vu3qX9am1PtPv9Q2llI1J9kryw7ZECNuwPc+Pjrb100BDMPF5phQAOs9dSQ4spXSVUqYlOSnJN/sXKKVM7/f6LUn2TE83XgCYULSUAkCHqbU+X0o5I8kN6fmC+fJa68pSysLe9YuTnNz7vGiSbEzywVrri81EDGNL6ydMLpJSWtM7eNK8p4boU7/iW20OBoBa6/IkywcsW9zv9UVJLmp3XEw+Q3W1vXT+7AYiASYr3XcBAABojJZSAACAcTLU6L1G7t2allIAAAAaIykFAACgMS113y2lHJ1kUZKXJbms1vqZAetPSfLfk5QkLyT5nVrryjGOFQCAScyoujA1jZiUllJ2TbI4yRFJfpTkjlLKjbXWVf2KPZzkHbXWn5RS3pPkfyfpHo+AAWhI7yjcLZl7zvjFAUwKElBgi1ZaSg9N8kCt9dEkKaVcleTYJH1Jaa31H/uVvz2JJ3cBABhyShmYDIYawIgd00pS2pXk0X7vNySZs43yv5Nk2VAreifxPj1JZsyY0VqEAAB0hiF7VPxK28MAOsuYDnRUSpmTZEF6ni8dpNa6pNY6q9Y6a/r06WN5aAAAADpQKy2lG5Ls0+99V++yrZRSDk5yaZL31FqfHJvwtt9wXUQunT97h7dvddsm9Y973lODuxIs+8F/rO+EzwMATHAtPmfu2VFgJK0kpXclObCU0pXkx0lOSrKwf4FSyowkX0tySq314TGPcrIY6pd304OBTMSYAACAKWPEpLTW+nwp5YwkN6Snu+/ltdaVpZSFvesXJzk3ySuTXFJKSZIXa62zxi9sAAAAJoOW5imttS5PsnzAssX9Xp+W5LSxDQ0AAIDJrqWkFACAyWMsxuAYOIbF8d1mBAR2zJiOvgsAAADbQ0spg7U4mh4AAMBoSUoBABi161cPnpIOoBWS0gmuXfOmqkhgahlpbmMAgHaRlAIA0GOIR3jmPfVYlr3i1AaCAaYKSSkAANs076nLmg4BmMQkpQBTiD8sAYCJxpQwAAAANEZSCgAAQGN03wUAIInR+IFmaCkFAACgMZJSAAAAGqP77ngZYp6vHSnXf1J7c4QB7dZqV77ju/feekHv77b+22/5HXbp/NljExwAMClISgEAJrEFS+8ecrkpooCJQvddAAAAGiMpBQAAoDG67wIATHK66kJnGGosh0HjNkxCktJhbPXLe8W3hi8495zxD2ZbegcT6T8g0lBa/jwAQOcaYgDFkf5GAGia7rsAAAA0RlIKAABAY6ZM992hhkM3Vx4AAECzpkxSCgDQlPH4cnzIfb52VLsEJqCpMPiRpLSD9A1WZKAiAABgkpCUAgB0oKGmebn+qQYCARglAx0BAADQGC2lAAATyRBzjTY+LzrAONJSCgAAQGMmR0tpv28U5z3VMzrVslecOqjYoGcvhhgwaMv2O3JsAIDRWLD07iH/Fjk+/t4A/sNkG5F3ciSlE9hQF0wy9EUzXNkdLdcuEy0eAOgIQ3yxvd1fjgNMApJSAGDKGo/5QwHYPpJSAIBxNtT0LXlt613t9EqCycU9vTUDHQEAANAYLaUAQEfolK62Q8U5r4E4gKmlkwc/kpQCAOyoIUfh/5W2hwHQySSlAABjaMjnRwEYlqQUAGAEQ3XJTUzhAkxs2zM9ZZMMdAQAAEBjpnRLaatDMbf60PD2DO3c5DDQk2kI6sn0WQAAYCqa0kkpADB1DNcFF4BmSUoBgClryEGJVnxriHJ65gCMF0kpANAW29NS2eT8o+16NMQjKAA9JKUAwJRgqhaAiUlSCgB0rhUXNB0BAKM0aZNS34YCTCx9v5cHPq8395z2B8O4G+2gQq1u3+oI+QBMXJM2KQUAOtNov1j2rCZAZ5GU7iAVHgAAwOhJSgGAMddqa+eyV5w6zpEAMNFJSgGAxhgDAoCdmg4AAACAqUtLKQAAwBQy0UYu11IKAABAY1pKSkspR5dS1pZSHiylfGKI9aWUclEpZV0p5d5SyiFjHyoAsIW6GYDJYsSktJSya5LFSd6T5OAkJwxRsb0/yWuTHJBkQZK/GuM4AYBe6mYAJpNWWkoPTfJArfXRWuumJFclOXZAmWOTXF57rEqycyllnzGOFQDooW4GYNJoZaCjriSP9nu/IcmcFsoMXJZSyulJTu99+2wp5bvbE+wwXpXkX8dgP03o1Ng7Ne6kc2Pv1LiTzo29U+NOOi72T255Mdq4Xzv6WDqGunl8dGrcSefGLu72End7dWrcST45VrG3VDe3dfTdWuuSJEvGcp+llJW11lljuc926dTYOzXupHNj79S4k86NvVPjTjo39k6Nu9Opm/9Dp8addG7s4m4vcbdXp8adtD/2VrrvbkjSv7tPV++y7S0DAIwNdTMAk0YrSeldSQ4spXSVUqYlOSnJNweUWZ7kQ0nSO9DC5lrrowEAxoO6GYBJY8Tuu7XW50spZyS5IT1J7OW11pWllIW96xcnuTbJ3FLKuiQbk/zWOMY80Jh2OWqzTo29U+NOOjf2To076dzYOzXupHNj79S4207dPG46Ne6kc2MXd3uJu706Ne6kzbGXWms7jwcAAAB9Wum+CwAAAONCUgoAAEBjOiIpLaX8RinlgVLK5lLKsEMTl1KOLqWsLaU8WEr5RL/lv1BKuamUcn8p5cZSyivaFPeIxy2lvKmUsrrfz09LKWf1rvtUKeWxfuuOaUfcrcbeW259b5nVpZSV27t9E3GXUvYppdzae608XEr5g37r2n7Oh7tu+60vpZSLSinrSin39g5Y0tK2Dcd9Su//w9pSyj39793hrpt2aSH2OaWUn/S7Ds5tdduG4/54v5jXllJeKqX8Qu+6xs55KeXLpZTHSylrh1k/Ia9xtq2om9taN7d6voa715s6360eu0yQurmF37MT9vdVC7FPyHq5hbjVyWMX88Stj2utE/4nyX5J3pTkliSzhimza5L16Rn+flqSlUkO6V33F0l+v/f12UkualPc23XcJC9L8qMkr+19/6kkH2vonLcUe+85f9VoP3s7406yV5KDe1//fJJ/StLdxDnf1nXbr8yvJ7k+SUlySJL7Wt224bgPTfKfel+/J8nqka6bCRT7nCRf35Ftm4x7QPn3Jvn2BDnn7+y9dtcOs37CXeN+Wvp/VTe393x3ZL3c6rEzAermFuuHCfn7qsXYJ1y93GLcc6JOHqu4J2x93BEtpbXWB2ut3x2h2KFJHqi1Plpr3ZTkqiTH9q47NslXel9f3m/5eNve4x6Z5Hu11h+Ma1StGe05m7DnvNb6o1rrmt7XzyRZk2TvNsU30Lau2y2OTc/ImrXWuirJzqWUfVrctrG4a63/WGv9Se/b29PcOR5oNOdtQp/zAU5OckVbIhtBrfXWJP+2jSIT8RpnBOrmtuvUermlY0+QurlT6+SWYp+g9bI6uY0mcn3cEUlpi7qS9J9/bUPvsiSZXmt9Ikl6//3FNsW0vcf9QAZfsGeWUh4qpXy1lPLK8QhyGK3GXpNs6ZLzX3dg+7G2XcctpcxMMjs9v5y3aOc539Z1O1KZVrYdL9t77N9Jsqzf++Gum3ZoNfa39XZT+XYppXs7tx0PLR+7lPJzSY5Oz5QgWzR5zkcyEa9xxoa6eex0ar283cdusG7u1Dp5W3ENZ6LUy+rkiVUnN3Z9jzhPabuUUr6Vnq4bA/1hrfX6dsfTqm3FvZ372SXJvCTn9Fv8hSR/nJ4L91NJLkrvROhjYYxif1ut9UellF9M8n9LKQ/VWm8amwiHNobnfM8k1yQ5q983h+N6zqeiUsqcJAuSvKPf4rZfN9vpniT71FqfK6UcleS6Usrrmw5qO7w3yT/UWvt/GzrRzzkTkLq5vXVzp9bLibq5k3RgvaxOngImTFJaa/2VUe5iQ3r6OW/R1bssSZ4opUyvtT5RSpme5PFRHqvPtuIupWzPcd+TZFWt9cf99v1Ev30tTs9zO2NmLGKvtf6o99/HSynXpOebzZsywc95KWVaer6xuqLW+rV++x7Xcz6EbV23A8vcOaDMtBa2HS+txJ1SysFJLk3ynlrrk1uWb+O6aYcRY+/tOrbl9Q2llI3p+WOrpc89Trbn2INadho+5yOZiNc4UTenzXVzp9bLYxX7BKibO7VO7h/XNo8/AetldfLEqpMbu74nU/fdu5IcWErp6v2ldlKSb/auW57kN3tf/2a/5eNte447qK9577cnW/x6knVjGt22jRh7KWWP3i4JKaXskZ6uCeta3X6ctBJ3Sc8v5Adrrf9rwLp2n/NtXbdbLE/vN8K9o6BtrrU+2uK2jcVdSpmR5GtJTqm1Ptxv+baum3ZoJfbp/V6/Jcme6fkjakKf8954/1OSd6VnoIIty5o+5yOZiNc4Y0PdPHY6tV5u6dgTpG7u1Do5rRx/gtbL6uSJVSc3d33XNo/6tCM/Sd6Xnmz8hSQ/TnJD7/LXJFner9wxSR5I8mB6uhZtWf7KJN9Kcn/vv7/QpriHPO4Qce+R5Mn0jojWb/nl6XnQ/6EkN6an60K7zvmIsSd5fW9896VnlLw/SVIm+jlPT3eV2hv76t6fY5o650Ndt0kWJlnY+7qkp+vSut5YZ21r2zZeIyPF/b+TPNXvHK8c6bqZQLF/NMna3p9VSeZ0wjnvfT8/yZUDtmv0nKfnj/p/SbIpPb/LF3TCNe5nxP9XdXMb6+ZW4t7Wvd7U+d6O2CdE3TzS79mJ/PuqhdgnZL3cQtzq5LGLecLWx1t+UQEAAEDbTabuuwAAAHQYSSkAAACNkZQCAADQGEkpAAAAjZGUAgAA0BhJKQAAAI2RlAIAANCY/x9KsJg+IGRBWQAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 1152x576 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.subplot(1,2,1)\n",
"plt.hist(costheta_Ds,alpha=0.7,density=True,bins=70, label=\"MC Ds tuple\");\n",
"plt.hist(costheta_Dplus,alpha=0.5,density=True,bins=70, label=\"MC Dplus tuple\");\n",
"plt.legend()\n",
"plt.subplot(1,2,2)\n",
"plt.hist(costheta_MC,alpha=0.7,density=True,bins=70, label=\"MC Ds+Dplus tuple\");\n",
"plt.hist(costheta_data,alpha=0.5,density=True,bins=70, label=\"data before cut\");\n",
"plt.legend()\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(16,8)"
]
},
{
"cell_type": "code",
"execution_count": 142,
"metadata": {},
"outputs": [],
"source": [
"MC_Dplus_tuple_dict_presel_9[\"cos_thetal\"]=costheta_Dplus\n",
"MC_Ds_tuple_dict_presel_9[\"cos_thetal\"]=costheta_Ds\n",
"data_tuple_dict_presel_4[\"cos_thetal\"]=costheta_data"
]
},
{
"cell_type": "code",
"execution_count": 143,
"metadata": {},
"outputs": [],
"source": [
"MC_Dplus_tuple_dict_presel_10={}\n",
"\n",
"\n",
"indices=np.where(MC_Dplus_tuple_dict_presel_9[\"cos_thetal\"]!=-2)[0]\n",
"labels = return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"\n",
"for label in labels: \n",
" MC_Dplus_tuple_dict_presel_10[label]=MC_Dplus_tuple_dict_presel_9[label][indices]\n",
" \n",
"\n",
"MC_Ds_tuple_dict_presel_10={}\n",
"\n",
"\n",
"indices=np.where(MC_Ds_tuple_dict_presel_9[\"cos_thetal\"]!=-2)[0]\n",
"labels = return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"\n",
"for label in labels: \n",
" MC_Ds_tuple_dict_presel_10[label]=MC_Ds_tuple_dict_presel_9[label][indices]\n",
" \n",
"data_tuple_dict_presel_5={}\n",
"\n",
"\n",
"indices=np.where(data_tuple_dict_presel_4[\"cos_thetal\"]!=-2)[0]\n",
"labels = return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"\n",
"for label in labels: \n",
" data_tuple_dict_presel_5[label]=data_tuple_dict_presel_4[label][indices] "
]
},
{
"cell_type": "code",
"execution_count": 144,
"metadata": {},
"outputs": [],
"source": [
"data_tuple_dict=data_tuple_dict_presel_5\n",
"MC_Ds_tuple_dict=MC_Ds_tuple_dict_presel_10\n",
"MC_Dplus_tuple_dict=MC_Dplus_tuple_dict_presel_10"
]
},
{
"cell_type": "code",
"execution_count": 150,
"metadata": {},
"outputs": [],
"source": [
"if l_flv[l_index]=='mu':\n",
" \n",
" lower_Dplus_mass_data=1830\n",
" lower_Dplus_mass_mc=lower_Dplus_mass_data\n",
" \n",
" upper_Dplus_mass=1910\n",
" upper_Dplus_mass_mc=upper_Dplus_mass_data\n",
" \n",
" lower_Ds_mass = 1930\n",
" lower_Ds_mass_mc=lower_Ds_mass_data\n",
" \n",
" upper_Ds_mass = 2010\n",
" upper_Ds_mass_mc=upper_Ds_mass_data\n",
" \n",
" \n",
"if l_flv[l_index]=='e':\n",
" \n",
" lower_Dplus_mass_data=1810\n",
" lower_Dplus_mass_mc=1790\n",
" \n",
" upper_Dplus_mass_data=1920\n",
" upper_Dplus_mass_mc=1950\n",
"\n",
" lower_Ds_mass_data = 1920\n",
" lower_Ds_mass_mc=1890\n",
" upper_Ds_mass_data = 2020\n",
" upper_Ds_mass_mc=2050\n",
"\n",
"#Retrieve mc signal and data bkg events\n",
"\n",
"data_bkg_indices=[]\n",
"MC_Ds_sig_indices=[]\n",
"MC_Dplus_sig_indices=[]\n",
"\n",
"\n",
"for i in range(len(data_tuple_dict[\"Ds_ConsD_M\"])):\n",
" \n",
" m = data_tuple_dict[\"Ds_ConsD_M\"][i]\n",
" \n",
" #selecting the out of signal regions\n",
" if 0<m<lower_Dplus_mass_data or upper_Dplus_mass_data < m < lower_Ds_mass_data or upper_Ds_mass_data < m:\n",
" data_bkg_indices.append(i)\n",
" \n",
"for i in range(len(MC_Ds_tuple_dict[\"Ds_ConsD_M\"])):\n",
" \n",
" m = MC_Ds_tuple_dict[\"Ds_ConsD_M\"][i]\n",
" \n",
" if lower_Ds_mass_mc< m <upper_Ds_mass_mc:\n",
" MC_Ds_sig_indices.append(i) \n",
" \n",
"for i in range(len(MC_Dplus_tuple_dict[\"Dplus_ConsD_M\"])):\n",
" \n",
" m = MC_Dplus_tuple_dict[\"Dplus_ConsD_M\"][i]\n",
" \n",
" #selecting the signal regions\n",
" if lower_Dplus_mass_mc< m <upper_Dplus_mass_mc:\n",
" MC_Dplus_sig_indices.append(i)\n",
" \n",
"\n",
"data_tuple_bkg={}\n",
"MC_Ds_tuple_sig ={}\n",
"MC_Dplus_tuple_sig ={}\n",
"\n",
"labels=return_branches(data_index=1, mother_index=0, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"\n",
"for label in labels: \n",
" data_tuple_bkg[label] = data_tuple_dict[label][data_bkg_indices]\n",
"\n",
"labels=return_branches(data_index=0, mother_index=0, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"for label in labels:\n",
" MC_Ds_tuple_sig[label] = MC_Ds_tuple_dict[label][MC_Ds_sig_indices]\n",
" \n",
"labels=return_branches(data_index=0, mother_index=1, l_index=l_index, meson_index=0)+['cos_thetal']\n",
"\n",
"for label in labels: \n",
" MC_Dplus_tuple_sig[label] = MC_Dplus_tuple_dict[label][MC_Dplus_sig_indices] "
]
},
{
"cell_type": "code",
"execution_count": 151,
"metadata": {},
"outputs": [],
"source": [
"with open('/disk/lhcb_data/davide/Rphipi/BDT_training/'+l_flv[l_index]+l_flv[l_index]+'/MC_for_BDT_training_Ds_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(MC_Ds_tuple_sig, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
" \n",
"with open('/disk/lhcb_data/davide/Rphipi/BDT_training/'+l_flv[l_index]+l_flv[l_index]+'/MC_for_BDT_training_Dplus_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(MC_Dplus_tuple_sig, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
"\n",
"with open('/disk/lhcb_data/davide/Rphipi/BDT_training/'+l_flv[l_index]+l_flv[l_index]+'/data_for_BDT_training_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(data_tuple_bkg, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
" \n",
"with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/'+l_flv[l_index]+l_flv[l_index]+'/MC_Ds_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(MC_Ds_tuple_dict, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
" \n",
"with open('/disk/lhcb_data/davide/Rphipi/MC/for_fit/'+l_flv[l_index]+l_flv[l_index]+'/MC_Dplus_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(MC_Dplus_tuple_dict, handle, protocol=pickle.HIGHEST_PROTOCOL)\n",
"\n",
"with open('/disk/lhcb_data/davide/Rphipi/data_for_BDT_selection/'+l_flv[l_index]+l_flv[l_index]+'/data_for_BDT_selection_'+l_flv[l_index]+l_flv[l_index]+'.pickle', 'wb') as handle:\n",
" pickle.dump(data_tuple_dict, handle, protocol=pickle.HIGHEST_PROTOCOL)"
]
},
{
"cell_type": "code",
"execution_count": 152,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"12244"
]
},
"execution_count": 152,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"data_tuple_bkg[\"Ds_ConsD_M\"].shape[0]"
]
},
{
"cell_type": "code",
"execution_count": 153,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2025"
]
},
"execution_count": 153,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MC_Dplus_tuple_sig[\"Dplus_ConsD_M\"].shape[0]"
]
},
{
"cell_type": "code",
"execution_count": 154,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1262"
]
},
"execution_count": 154,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MC_Ds_tuple_sig[\"Ds_ConsD_M\"].shape[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.15"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini