{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"WARNING: The TensorFlow contrib module will not be included in TensorFlow 2.0.\n",
"For more information, please see:\n",
" * https://github.com/tensorflow/community/blob/master/rfcs/20180907-contrib-sunset.md\n",
" * https://github.com/tensorflow/addons\n",
"If you depend on functionality not listed there, please file an issue.\n",
"\n"
]
}
],
"source": [
"import numpy as np\n",
"from math import pi\n",
"import matplotlib.pyplot as plt\n",
"\n",
"import tensorflow as tf\n",
"import zfit\n",
"import os\n",
"\n",
"os.environ[\"CUDA_VISIBLE_DEVICES\"] = \"1\"\n",
"\n",
"from utils.kin_utils_edit import *\n",
"from utils.BtoKllgamma_fullydiff_utils_fix import *\n",
"\n",
"import pickle\n",
"ztf = zfit.ztf\n",
"ztyping = zfit.util.ztyping\n",
"ztypes = zfit.settings.ztypes\n",
"\n",
"zfit.run.check_numerics = False\n",
"zfit.settings.set_verbosity(10)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"#TASK='LOAD'\n",
"TASK='SAMPLE'\n",
"n_events=100000\n",
"\n",
"lepton_flv='mu'\n",
"\n",
"mB_mass=5280\n",
"mKst_mass=892\n",
"mgamma_mass = 0.1\n",
"\n",
"neutral_mesons=True"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"if lepton_flv=='e':\n",
" l_mass = 0.5\n",
"elif lepton_flv=='mu':\n",
" l_mass = 105\n",
"\n",
"mBbar_mass=mB_mass-mgamma_mass\n",
"\n",
"lower_q2= 4*l_mass*l_mass\n",
"upper_q2=(mBbar_mass-mKst_mass)*(mBbar_mass-mKst_mass)\n",
"\n",
"\n",
"upper_mBbar2 = mBbar_mass**2\n",
"lower_mBbar2 = (2*l_mass+mKst_mass)*(2*l_mass+mKst_mass)\n",
"\n",
"if neutral_mesons:\n",
" filename='MC_tuple_{0}k_B0toK0'+lepton_flv+lepton_flv+'_fulldiff_{1}MeVcutoff.pickle'.format(n_events/1000, mgamma_mass)\n",
"else:\n",
" filename='MC_tuple_{0}k_B+toK+'+lepton_flv+lepton_flv+'_fulldiff_{1}MeVcutoff.pickle'.format(n_events/1000, mgamma_mass)\n"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"class dGamma(zfit.pdf.ZPDF):\n",
" \n",
" _PARAMS = ['mB','mKst','ml','mgamma']\n",
" _N_OBS = 5\n",
" \n",
" def _unnormalized_pdf(self, x):\n",
" \n",
"\n",
" mB = self.params['mB'] \n",
" mKst = self.params['mKst']\n",
" ml = self.params['ml']\n",
" mgamma = self.params['mgamma']\n",
" \n",
" q2, mBbar2, cos_theta_l, cos_theta_gamma, phi_gamma, = x.unstack_x()\n",
" \n",
" pB_Bbar, pKst_Bbar, k_Bbar, q_Bbar, p1_Bbar, p2_Bbar, p1, p2 = momenta_map(q2,\n",
" cos_theta_l,#phi_l,\n",
" cos_theta_gamma,phi_gamma,\n",
" #cos_theta_Kst,phi_Kst,\n",
" mB,mBbar2,mKst,\n",
" ml, mgamma, n_events=n_events, mode='tf')\n",
" \"\"\"\n",
" Scalar products\n",
" \n",
" \"\"\"\n",
" \n",
" \n",
" \n",
" p1k=scalar_product_4(p1_Bbar, k_Bbar)\n",
" p2k=scalar_product_4(p2_Bbar, k_Bbar)\n",
" \n",
" Bk=scalar_product_4(pB_Bbar, k_Bbar)\n",
" Kk=scalar_product_4(pKst_Bbar, k_Bbar)\n",
" \n",
" p1p2 = scalar_product_4(p1,p2)\n",
" p1B = scalar_product_4(p1_Bbar,pB_Bbar)\n",
" p2B = scalar_product_4(p2_Bbar,pB_Bbar)\n",
" \n",
" p1K = scalar_product_4(p1_Bbar,pKst_Bbar)\n",
" p2K = scalar_product_4(p2_Bbar,pKst_Bbar)\n",
" \n",
" BK = scalar_product_4(pB_Bbar,pKst_Bbar)\n",
"\n",
"\n",
" \"\"\"\n",
" R_{ij} functions\n",
" \n",
" \"\"\"\n",
" R12= (p1p2)/((p1k)*(p2k))\n",
" \n",
" RBK=(BK)/((Bk)*(Kk))\n",
" \n",
" RB1=(p1B)/((p1k)*(Bk))\n",
" RB2=(p2B)/((p2k)*(Bk))\n",
" \n",
" RK1=(p1K)/((p1k)*(Kk))\n",
" RK2=(p2K)/((p2k)*(Kk))\n",
" #\n",
" R11= (ml**2)/((p1k)*(p1k))\n",
" R22= (ml**2)/((p2k)*(p2k))\n",
" \n",
" RKK=(mKst**2)/((Kk)*(Kk))\n",
" RBB=(mB**2)/((Bk)*(Kk))\n",
" \n",
" \"\"\"\n",
" Phasespace\n",
" \n",
" \n",
" \"\"\"\n",
" beta_l=dphi2_tf(ztf.sqrt(q2), ml, ml)\n",
" phasespace_term= tf.expand_dims(dphi2_tf(mB, ztf.sqrt(mBbar2), mgamma)*dphi2_tf(ztf.sqrt(mBbar2), mKst, ztf.sqrt(q2))*beta_l,axis=-1)\n",
" \n",
" \"\"\"\n",
" Total PDF\n",
" \n",
" \"\"\"\n",
" if neutral_mesons==True:\n",
" pdf=-(R11+R22-2*R12)*phasespace_term\n",
" else:\n",
" pdf=-(RBB+RKK+R11+R22 + (RB1-RB2+RK1-RK2) + 2*RBK -2*R12)*phasespace_term\n",
" \n",
" return pdf[:,0]\n"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /disk/lhcb_data/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/tensorflow/python/ops/resource_variable_ops.py:435: colocate_with (from tensorflow.python.framework.ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Colocations handled automatically by placer.\n"
]
}
],
"source": [
"mB_par = zfit.Parameter('mB', mB_mass)\n",
"mKst_par = zfit.Parameter('mKst', mKst_mass)\n",
"mgamma_par = zfit.Parameter('mgamma', mgamma_mass)\n",
"mlepton_par = zfit.Parameter('ml', l_mass)\n",
"mBbar_par = zfit.Parameter('mBbar', mBbar_mass)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"scrolled": true
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /disk/lhcb_data/davide/miniconda3/envs/zfit_env/lib/python3.6/site-packages/zfit/core/sample.py:163: to_int32 (from tensorflow.python.ops.math_ops) is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Use tf.cast instead.\n"
]
}
],
"source": [
"lower = ((lower_q2, lower_mBbar2, -1. , -1. ,-pi,),)\n",
"upper = ((upper_q2, upper_mBbar2, 1. , 1. , pi,),)\n",
"\n",
"\n",
"#lower = ((lower_q2, lower_mBbar2, -1., -pi, -1., -pi, -1., -pi,),)\n",
"#upper = ((upper_q2, upper_mBbar2, 1., pi, 1., pi, 1., pi,),)\n",
"#\n",
"#obs_list=[\"q2\", \"mBbar\", \"theta_l\", \"theta_gamma\", \"phi_gamma\"]\n",
"#obs_list=[\"q2\", \"mBbar2\", \"cos_theta_l\", \"phi_l\",\"cos_theta_gamma\", \"phi_gamma\",\"cos_theta_Kst\",\"phi_Kst\"]\n",
"obs_list=[\"q2\", \"mBbar2\", \"cos_theta_l\",\"cos_theta_gamma\", \"phi_gamma\"]\n",
"obs = zfit.Space(obs_list, limits=(lower,upper))\n",
"\n",
"pdf = dGamma(obs=obs, ml = mlepton_par, mB = mB_par, mKst = mKst_par, mgamma= mgamma_par)\n",
"\n",
"sampler = pdf.create_sampler(n=n_events)\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"sampling done!\n"
]
}
],
"source": [
"if TASK=='SAMPLE':\n",
" \n",
" for i in range(1):\n",
" sampler.resample()\n",
" sample = sampler.to_pandas()\n",
" print('sampling done!')\n",
" MC_tuple={}\n",
" \n",
" for j in range(sample.keys().shape[0]):\n",
" \n",
" MC_tuple[sample.keys()[j]]=np.array([sample[sample.keys()[j]][i] for i in range(n_events)]).reshape(n_events,1)\n",
" \n",
" "
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"pB_Bbar, pKst_Bbar, k_Bbar, q_Bbar, p1_Bbar, p2_Bbar, p1, p2=momenta_map(MC_tuple['q2'], \n",
" MC_tuple['cos_theta_l'],\n",
" #MC_tuple['phi_l'],\n",
" MC_tuple['cos_theta_gamma'], \n",
" MC_tuple['phi_gamma'],\n",
" #MC_tuple['cos_theta_Kst'], \n",
" #MC_tuple['phi_Kst'],\n",
" mB_mass, MC_tuple['mBbar2'], \n",
" mKst_mass, l_mass, mgamma_mass,\n",
" mode = 'numpy')"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"p1_Bbar_to_boost=p1_Bbar\n",
"p2_Bbar_to_boost=p2_Bbar\n",
"p1_q = lorentz_boost_np(p1_Bbar_to_boost, -q_Bbar/q_Bbar[:,0:1])\n",
"p2_q = lorentz_boost_np(p2_Bbar_to_boost, -q_Bbar/q_Bbar[:,0:1])\n",
"z=np.array([[0.,0.,0.,1.] for i in range(p1_Bbar.shape[0])])\n",
"\n",
"cos_theta_l_np_Bbar=get_costheta_l_np(p1_Bbar,z)\n",
"cos_theta_l_np_q=get_costheta_l_np(p1_q,z)\n",
"\n",
"\n",
"MC_tuple['pB_Bar']=pB_Bbar\n",
"MC_tuple['pKst_Bar']=pKst_Bbar\n",
"MC_tuple['k_Bar']=k_Bbar\n",
"MC_tuple['q_Bar']=q_Bbar\n",
"MC_tuple['p1_Bar']=p1_Bbar\n",
"MC_tuple['p2_Bar']=p2_Bbar\n",
"MC_tuple['p1']=p1\n",
"MC_tuple['p2']=p2\n",
"MC_tuple['cos_theta_l(3)']=cos_theta_l_np_q.reshape(n_events,1)\n",
"MC_tuple['cos_theta_l(2)']=cos_theta_l_np_Bbar.reshape(n_events,1)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"if TASK=='LOAD':\n",
" with open(filename,'rb') as f:\n",
" MC_tuple= pickle.load(f) "
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(104.99999999999996, 104.99999999999993)"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"np.sqrt(scalar_product_4_np(p1_Bbar,p1_Bbar)).mean(),np.sqrt(scalar_product_4_np(p1_q,p1_q)).mean()"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA3gAAANCCAYAAAA0j26FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzde3xU1b3///eHIdySkEi5gxChAUXkokGxtj1cvAChYD2Kl0rRtiot2GrbX8X6a9X+Dj9pT0F6FAVEqFQONNRUOQiKX49HSxUVAZGIXI4SLlIEKZAEEAjr+8dM4sxkkkxuc9nzej4e88jM3mvv+czO7DX7s9faa5tzTgAAAACA5Ncs3gEAAAAAABoHCR4AAAAAeAQJHgAAAAB4BAkeAAAAAHgECR4AAAAAeETzeAcQSfv27V1OTk68wwAANLH33nvvkHOuQ7zjSBb8PgJA6qjvb2RCJng5OTlav359vMMAADQxMyuOdwzJhN9HAEgd9f2NpIsmAAAAAHgECR4AAAAAeAQJHgAAAAB4BAkeAAAAAHgECR4AAAAAeAQJHgAAAAB4BAkeAAAAAHgECR4AAAAAeAQJHgAAAAB4RPN4B4C6y5n2YsjrXTPy4xQJAAAAgODj83gfm9OCBwAAAAAeQQsegCbnnFNJSYmOHTum48ePq7y8PN4hoZH4fD61adNGbdu2VWZmpsws3iEhAvZBb2L/AxAJCR6AJuWc02effaaysjK1a9dOnTt3ls/n40DEA5xzKi8vV2lpqQ4dOqQTJ06oY8eO/G8TDPugN7H/AagOCR6AJlVSUqKysjL17NlTPp8v3uGgEZmZmjdvruzsbGVmZqq4uFglJSVq27ZtvENDEPZBb2L/A1AdrsED0KSOHTumdu3acWDpcT6fT+3atdOxY8fiHQrCsA96H/sfgGAkeACa1PHjx5WRkRHvMBADGRkZOn78eLzDQBj2wdTA/gegAgkegCZVXl5Oy0GK8Pl8DN6RgNgHUwP7H4AKJHgAmhwX/acG/s+Ji/+N9/E/BlCBBA8AAAAAPIIEDwAAAAA8gtskoNHkTHsx5PWuGflxigQAAABITbTgAQAAAIBHRJXgmdkoM9tmZjvNbFqE+WZm/xGYv9nMLg6ad6+ZFZnZFjNbamatGvMDAAAAAAD8ak3wzMwnaY6k0ZL6SbrZzPqFFRstKTfwuFPSk4Flu0n6saQ851x/ST5JNzVa9B6WM+3FkEeqvDeA+pk+fbrMTAcOHIh3KEBKYh8EkCiiacG7VNJO59zHzrlTkpZJGh9WZrykxc5vnaRsM+sSmNdcUmszay6pjaRPGyl2AEDA5s2b1alTJ3Xq1EmStGPHDplZyKNTp04aM2aMNm7cGOdoAe9hHwSQKKIZZKWbpD1Br/dKuiyKMt2cc+vN7PeSdks6IWmNc25NpDcxszvlb/1Tjx49ooseACDJf3A5cODAytebNm2SJE2ZMkVDhw5VeXm5PvjgA82ZM0dXXnmltmzZoi5dulS3OgB1xD4IIFFEk+BFunOmi6aMmZ0jf+veeZKOSFpuZrc6556tUti5+ZLmS1JeXl74+gEA1Th58qR27NihcePGVU6raCG44447Qg46MzIy9PDDD+ull17S7bffHvNYAS9iHwSQSKLporlX0rlBr7urajfL6spcKekT59xB59xpSYWSvlb/cAEgta1du1YjR45Uenq6unbtqhkzZqioqEjl5eUhB5EbN25UixYt1K9f6CXTvXr1kiQdPXo0pnGnGjM718xeM7OtgYHGfhKhTLUDlCFxsQ8CSHTRtOC9KynXzM6TtE/+QVJuCSuzQtJUM1smf/fNo865/Wa2W9JQM2sjfxfNkZLWN1r0AJBCCgsLdeONN2rw4MF65JFHdOrUKc2aNUsrV66UJA0aNKiy7MaNG9W/f3+lpaWFrGP16tWSpKFDh8Yu8NR0RtLPnHMbzCxT0ntm9opz7sOgMsEDlF0m/wBl4ZdAIIGwDwJIBrUmeM65M2Y2VdLL8o+CudA5V2RmkwPz50paJWmMpJ2Sjku6PTDvbTP7i6QN8v/YbVSgGyYAIHrFxcWaNGmS8vPzVVhYqGbN/B0whg8frry8PLVq1Up9+/aVJO3fv18HDhzQiBEjdOjQITnntG/fPj311FNatmyZJk2axMFlE3PO7Ze0P/C8xMy2yn+9enCCVzlAmaR1ZpZtZl0CyyLBsA8CSBbRtODJObdK/iQueNrcoOdO0pRqln1Q0oMNiBGQpBpv2bBrRn4MIwFib/r06Tpx4oRmzpxZeWApSZdccomysrKUm5srn88n6ctrf5YuXaqlS5dWlm3Xrp2efPJJ3XXXXSHrPnLkiO68806tXr1amZmZ+sUvfqF77rknBp8qNZhZjqTBkt4OmxVxgDIFEsOg5RmELAGwDwJIFlEleACA+CkvL9dzzz2nMWPGqHfv3hHnh1/7I0nz5s1Tr169VFZWpuXLl2vJkiXat2+fzELHxZo6daq++OIL7du3T8XFxRo5cqT69u2r0aNHN+0HSwFmliHpOUn3OOeOhc+OsEiVQcYYhCz+2AcBJBMSPABIcHv27NHhw4c1ZMiQKvOKiopUWlpa5eDS5/Np4sSJat26tSRp3Lhx+uijj/Too4/qvvvuU0ZGhiRVHni+9957atu2rS666CLdcccdWrhwIQeXDWRmafInd0ucc4URikQziBkSAPsggGRCggcgIdTUBTeZNEV34bKysmrnzZ49W1LVwR1yc3MrDywlycx022236e6779bKlSt10003SZK2b9+us2fPqn///pVlBw0apMLCSPkIomX+JpqnJW11zs2qpljEAcpiFWM49sHqsQ8CSCbR3CYBqFbOtBcrHwCaRs+ePWVmeuWVV0Kmr1u3Ts8884wkacCAAZL8Q69/8sknIa0JFcaOHStJeuGFFyqnlZaWKisrK6Rcdna2SkpKGvUzpKArJE2UNMLMNgUeY8xscsUgZfJf2/6x/AOUPSXpR3GKFbVgHwSQTGjBA4AEl5GRoRtuuEEFBQW67rrrNGrUKG3btk3z589XZmamsrKyKg8QN23aJOdcxIPLnJwcXXDBBVq9erVOnz6ttLQ0ZWRk6Nix0EvDjh49qszMzJh8Nq9yzq1V5GvsgstUO0AZEgv7IIBkQgseACSBefPmaeLEiXrjjTf005/+VO+8844KCgqUnp4ecXCHitaEcKNHj9bRo0f1+uuvS5L69OkjM1NRUVFlmU2bNoV0FwPAPgggedCCl0Do5hh74duc2y0gUWVnZ2vx4sVVpu/evTvk9T333FPj8OozZ87UzJkzK1+np6fr+uuv1wMPPKA//elPKi4u1oIFC7Ro0aLGCx7wAPZBAMmCFjwASHFz5sxRWlqaunTpoquuukrTpk1j9D4ghtgHATQmWvAAIMVlZ2dr+fLl8Q4DSFnsgwAaEwleHCVjl8xkjBkAAABIFXTRBAAAAACPoAXPYxg0BAAAAEhdtOABAAAAgEfQggcACc6sxvtlx4T/ntxAamIfBJBMSPCQcIK7mdLFFODADog39kEAyYQumgAAAADgEbTgIe5quvWCl2/LUNNno+USsfLFF19oypQpevXVV3Xo0CH16NFDDzzwgG655ZZ4hwZ4HvsfgKZAgtdI6FYYX4weCtTPmTNn1LVrV7366qvKycnRm2++qfz8fJ133nm6/PLL4x0e4GnsfwCaAl00AcADpk+fLjPTgQMH6rRcenq6fvOb36hXr15q1qyZvv71r+uKK67Qm2++2USRAt5U3T548OBB/eMf/wh5nDp1ShL7H4CmQQseAHjA5s2b1alTJ3Xq1KlB6ykrK9P69ev1k5/8pJEiQ7Kg50PDVLcPDhkyRMXFxSHTXnvtNQ0bNqzKOtj/ADQGEjwACYGDy4bZvHmzBg4c2KB1nD17VrfddpuGDBmiq6++upEiA1JDdfvgkiVLdOLEiZBpkcqx/wFoLCR4iIl4XiPH9XnwupMnT2rHjh0aN25cvdfhnNPkyZP16aef6uWXX06I+34ByaKmffCKK66odXn2PwCNiWvwPC5n2ouVDwDJb+3atRo5cqTS09PVtWtXzZgxQ0VFRSovL6/SKvD3v/+9smz79u01ZcoUrV27VmamZ599trKcc05TpkzRpk2btHr1amVkZMT6YwFJI9p9cOrUqWrZsqX+/ve/a8yYMTrnnHOUnZ2tb33rW9qzZ09lOfY/AI2NFjwASBKFhYW68cYbNXjwYD3yyCM6deqUZs2apZUrV0qSBg0aVKXsJZdcohkzZuiLL77Qo48+qr/97W+SVOVAdN26dXr11VfVtm3b2H4oIInUZR/cvHmzWrZsqdGjR2vChAn63e9+p3feeUdPP/20rr/+er399tuS2P8AND4SPCBK3AoD8VRcXKxJkyYpPz9fhYWFatbM3wFj+PDhysvLU6tWrdS3b19J0t69e3Xbbbdp3LhxWr58eWXZb3zjGxo6dKhatGih888/v3K9TzzxhFq2bKlzzz238v1++ctf6pe//GWMPyWQuOqyD0r+BK+kpETPPfecrrvuOknSHXfcofLyci1atEgfffSRWrduzf4HoNGR4DUBEoHEFm131bp0a+V/jqY2ffp0nThxQjNnzqw8sJSkSy65RFlZWcrNzZXP55MkPfLIIzp58qRmzZoVUvayyy5Tq1atdP755ystLU2S1LNnTznnYvthgCRUl31w165dOnr0qCZMmFCZ3FUYPny4Fi1apOLiYl1zzTXsfwAaHdfgAUCCKy8v13PPPacxY8aod+/eEedXdLk8e/asCgoKNHbsWPXs2TPi+ho62iaQauqyD0rS+++/L0m6/fbbq5Q9e/asJHGtHYAmQ4IHAAluz549Onz4sIYMGVJlXlFRkUpLSysPLnfv3q1Dhw7p4osvrlJ2+/btOnnyJAkeUEd12QelLxO8SPvh22+/rbS0NPXv37/pAgaQ0kjwACDBlZWVVTtv9uzZkr4c3KGmsgsXLpRECx5QV3XZB6UvE7zmzUOvhDl8+LCeffZZjRo1SllZWZL8t1ho3ry5MjIyKh8jRozQvn37GvtjAEgRJHgpJPiWCeHXl3E7BT+2AxJRz549ZWZ65ZVXQqavW7dOzzzzjCRpwIABIWXfeOONkLJr166tPBCtKAsgOnXZB6UvE7zw/fD+++/XiRMn9NBDD4WUTUtL09GjR1VaWqp9+/bp888/1+9///sm+jQAvI5BVgAgwWVkZOiGG25QQUGBrrvuOo0aNUrbtm3T/PnzlZmZqaysrMrWgOCy3/nOdzRs2DBt2bJFy5cvV6dOnXTmzBm1b98+zp8ISC512QdLS0v18ccfa9CgQfre976noqIifeUrX9Hy5cv13//933r88cdDum5u2LBB559/fuUALVlZWerbt69Onz4dl88KIPmR4IVhNEQAiWjevHlq2bKlVq1apTVr1mjw4MEqKCjQXXfdVaXLZXDZVatW6corr9Qbb7yhr33ta8rLy4vTJwCSW7T74AcffCDnnP7t3/5N69ev15w5c3T48GENHDhQL7zwgsaNGxey3g0bNlRej1deXq7XXntNr7/+epXWQgCIFgke4oIukEDdZGdna/HixVWm7969O6qyBw8e1MGDB7n+DqinaPfBiu6ZF110kfLz8/Xggw/WuN6NGzeqqKhI//Vf/6WTJ09Kkp5++mm6UgOoNxK8eiJBAZBMPvjgA0kMsAI0tffff19t27ZVjx49ai17+vRpbdmyRS+99JKGDRums2fP6j//8z/1gx/8QKNHj1a7du1iEDEAr2GQFQBIAVu2bJFEggc0tffff18XXnhhVGWLiop06tSpymvymjVrpuHDh+vkyZM6cuRIU4YJwMNowUOToZUTSBxbtmxRq1at1KdPn3iHAniWc04ffPCBbr755qjKb9iwQX369FHbtm0lSZ9//rkeeOABXXzxxerVq1dThgrAw0jwUhgJWOyxzREv8+fP1/z58+MdBuBpZqaSkpKoy2/cuFE7d+5URkaGfD6fsrOzNXr0aK1ataoJowTgdSR4dcDBOQAAaCyPPfaYHnvssXiHAcBjuAYPAAAAADyCBA8AAAAAPIIumjFGN08AAAAATYUWPAAAAADwCBI8AAAAAPAIumgCSS642++uGflxjAQAAADxFlWCZ2ajJP1Bkk/SAufcjLD5Fpg/RtJxSbc55zYE5mVLWiCpvyQn6XvOubca7RMkOK6586Zo/6/1/f+HLxdt4lbf5Zqac07+agJe5pyLdwgAAKS8WrtomplP0hxJoyX1k3SzmfULKzZaUm7gcaekJ4Pm/UHSS8658yUNlLS1EeIGkCR8Pp/Ky8vjHQZioLy8XD6fL95hIAKSb+/jfwygQjTX4F0qaadz7mPn3ClJyySNDyszXtJi57dOUraZdTGztpK+KelpSXLOnXLOHWnE+AEkuDZt2qi0tDTeYSAGSktL1aZNm3iHgTCcZEkNnGABUCGaLprdJO0Jer1X0mVRlOkm6Yykg5IWmdlASe9J+olzriz8TczsTvlb/9SjR49o40cM0M3Ue2LZlbNt27Y6dOiQMjMzOfjwsPLych0+fFjt27ePdygJwcwWShor6TPnXP8I84dJekHSJ4FJhc653zRFLBUnWbKzs5ti9UgQnGABUCGaFrxIF86E9wOorkxzSRdLetI5N1hSmaRpkd7EOTffOZfnnMvr0KFDFGEBSAaZmZlKT09XcXGxjhw5ojNnztCVyCOcczpz5oyOHDmi4uJipaenKzMzM95hJYo/ShpVS5m/OecGBR5NktxJ/pMshw8fphXPwypOsLRt2zbeoQBIANG04O2VdG7Q6+6SPo2yjJO01zn3dmD6X1RNggfAm8xMHTt2VElJiY4dO6bPPvuMA00P8fl8atOmjdq3b6/MzEwG0wlwzr1hZjnxjkPyn2Q5ceKEiouL1a5dO2VkZMjn8/G/SnLOOZWXl6u0tFSHDx/mBAuAStEkeO9KyjWz8yTtk3STpFvCyqyQNNXMlsnfffOoc26/JJnZHjPr65zbJmmkpA8bLXoAScHM1LZtW84uA6EuN7P35T8h+nPnXFGkQg29hIGTLN7FCRYAkdSa4DnnzpjZVEkvy3+bhIXOuSIzmxyYP1fSKvlvkbBT/tsk3B60irslLTGzFpI+DpsHAEAq2iCpp3Ou1MzGSHpe/pGoq3DOzZc0X5Ly8vLq1b+ZkywAkDqiug+ec26V/Elc8LS5Qc+dpCnVLLtJUl4DYgQAwFOcc8eCnq8ysyfMrL1z7lA84wIAJL+oErxUlSyjRyZLnKi/4P9xoty8HED9mVlnSQecc87MLpV/0LPP4xwWAMADSPDgSSS9AOLJzJZKGiapvZntlfSgpDSpsgfM9ZJ+aGZnJJ2QdJNjeFkAQCMgwQMAoJE5526uZf7jkh6PUTgAgBQSzX3wAAAAAABJgAQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCUTSTBMP+o6H4DgEAAHgfLXgAAAAA4BEkeAAAAADgEXTRBBpZqnaFDP/cu2bkxykSAACA1EWCB6BJBCd8JHsAAACxQYIHIK5o+QMAAGg8XIMHAAAAAB6R8i14qXq9FAAAAADvSfkED0DNuJYOAAAgeZDgAYga18sBAAAkNq7BAwAAAACPoAUPSDJcNwoAAIDqpGSCxwEyUD32DwAAgORFF00AAAAA8IiUbMEDUkW0rXG02gEAAHgDCR6AeiMxBAAASCx00QQAAAAAjyDBAwAAAACPIMEDAAAAAI8gwQMAAAAAjyDBAwAAAACPIMEDAAAAAI/gNgkAElbwbRh2zciPYyQAAADJgRY8AAAAAPAIEjwAAAAA8Ai6aAJIKMHdMuMpPA66iAIAgGRACx4AAAAAeAQteABiLlFa6QAAALyGFjwAAAAA8AgSPAAAAADwCBI8AAAAAPAIrsED0OQa45o7RrUEAACoHS14AAAAAOARJHgAAAAA4BEkeAAAAADgEVEleGY2ysy2mdlOM5sWYb6Z2X8E5m82s4vD5vvMbKOZrWyswAEAAAAAoWodZMXMfJLmSLpK0l5J75rZCufch0HFRkvKDTwuk/Rk4G+Fn0jaKqltI8UNADEVPMgLA7wAAIBEFU0L3qWSdjrnPnbOnZK0TNL4sDLjJS12fuskZZtZF0kys+6S8iUtaMS4ASDp5Ex7sfIBbzOzhWb2mZltqWZ+jT1fAACor2gSvG6S9gS93huYFm2Z2ZJ+IelsTW9iZnea2XozW3/w4MEowgIAIGH9UdKoGuYH93y5U/6eLwAANFg0CZ5FmOaiKWNmYyV95px7r7Y3cc7Nd87lOefyOnToEEVYAAAkJufcG5IO11Ck2p4vAAA0RDQ3Ot8r6dyg190lfRplmesljTOzMZJaSWprZs86526tf8gAwDVxSHrV9XzZH17QzO6Uv5VPPXr0iElwAIDkFU0L3ruScs3sPDNrIekmSSvCyqyQ9N3ANQVDJR11zu13zt3vnOvunMsJLPffJHcAAETVO8Y/kR4uAIA6qLUFzzl3xsymSnpZkk/SQudckZlNDsyfK2mVpDGSdko6Lun2pgsZAGpG6x6SQDS9YwAAqLNoumjKObdK/iQueNrcoOdO0pRa1vE/kv6nzhECAOA9KyRNNbNl8t9W6Khzrkr3TAAA6iqqBA8AklX4LQlo0UMsmNlSScMktTezvZIelJQm0fMFANC0SPAAJD3uK4dE45y7uZb5tfZ8AQCgPkjwAKQUrs8DAABeFs0omgAAAACAJEALHoCUxfV5AADAa1IiweP6HAAAAACpgC6aAAAAAOARKdGCBwDxQO8BAAAQa7TgAQAAAIBHkOABAAAAgEeQ4AEAAACAR5DgAQAAAIBHkOABAAAAgEeQ4AEAAACAR5DgAQAAAIBHcB88AGhE3PsOAADEEy14AAAAAOARJHgAAAAA4BF00QSABqBLJgAASCS04AEAAACAR5DgAQAAAIBHkOABAAAAgEeQ4AEAAACAR5DgAQAAAIBHMIomAAQkyoiY4XHsmpEfp0gAAECyoQUPAAAAADyCBA8AAAAAPIIumgAQB3TDBAAATYEWPAAAAADwCFrwAKCOEmUwFgAAgHC04AEAAACAR3i2BY8z7ACSSX3rrODluI4PAADQggcAAAAAHuHZFjwA8CJ6JwAAgJrQggcAAAAAHkGCBwAAAAAeQYIHAAAAAB5BggcAAAAAHkGCBwAAAAAeQYIHAAAAAB5BggcAAAAAHkGCBwAAAAAeQYIHAAAAAB4RVYJnZqPMbJuZ7TSzaRHmm5n9R2D+ZjO7ODD9XDN7zcy2mlmRmf2ksT8AAAAAAMCv1gTPzHyS5kgaLamfpJvNrF9YsdGScgOPOyU9GZh+RtLPnHMXSBoqaUqEZQEATSBn2ouVD8ReFCdHh5nZUTPbFHj8Oh5xAgC8pXkUZS6VtNM597EkmdkySeMlfRhUZrykxc45J2mdmWWbWRfn3H5J+yXJOVdiZlsldQtbFgAATwk6OXqVpL2S3jWzFc658N+/vznnxsY8QACAZ0XTRbObpD1Br/cGptWpjJnlSBos6e1Ib2Jmd5rZejNbf/DgwSjCAgAgYVWeHHXOnZJUcXIUAIAmFU2CZxGmubqUMbMMSc9Jusc5dyzSmzjn5jvn8pxzeR06dIgiLAAAElY0J0cl6XIze9/MVpvZhZFWxAlQAEBdRNNFc6+kc4Ned5f0abRlzCxN/uRuiXOusP6hAgDqK/g6vF0z8uMYScqI5uToBkk9nXOlZjZG0vPyX8seupBz8yXNl6S8vLzwdQAAECKaBO9dSblmdp6kfZJuknRLWJkVkqYGrs+7TNJR59x+MzNJT0va6pyb1YhxAwDCMJhKQqn15Ghwjxbn3Coze8LM2jvnDsUoRgCAB9XaRdM5d0bSVEkvS9oqqcA5V2Rmk81scqDYKkkfS9op6SlJPwpMv0LSREkjgkYJG9PYHwIAgARTeXLUzFrIf3J0RXABM+scOBEqM7tU/t/kz2MeKQDAU6JpwZNzbpX8SVzwtLlBz52kKRGWW6vI3VQAAPAs59wZM6s4OeqTtLDi5Ghg/lxJ10v6oZmdkXRC0k2B31MAAOotqgQPAADUTRQnRx+X9His4wIAeFs0o2gCAAAAAJIACR4AAAAAeARdNAEgwTE6JgAAiBYJHgCkOO6RBwCAd9BFEwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPILbJAAAKoXfc4/bJgAAUFUi36OWFjwAAAAA8AgSPAAAAADwCBI8AAAAAPAIrsEDgBSTyNcNAACAhqEFDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8gkFWAADVCh6QhZueAwCQ+GjBAwAAAACPIMEDAAAAAI8gwQMAAAAAjyDBAwAAAACPIMEDAAAAAI8gwQMAAAAAj+A2CQAAAABQi+BbByUyWvAAAAAAwCNI8AAAAADAI+iiCQCISnjXlF0z8uMUCQAAqA4teAAAAADgEbTgAQAAAECYZBlUJRwteAAAAADgESR4AAAAAOARJHgAAAAA4BEkeAAAAADgEQyyAgAAAABK3oFVgtGCBwAAAAAeQQseAAAAgJTlhVa7YCR4AAAAAFKG1xK6cHTRBAAAAACPoAUPAAAAgKd4vZWuJiR4AAAAABJSKidq9RVVgmdmoyT9QZJP0gLn3Iyw+RaYP0bScUm3Oec2RLMsAABe1JDfTgBIViRk8VdrgmdmPklzJF0laa+kd81shXPuw6BioyXlBh6XSXpS0mVRLgsAgKc05Lcz1rECSF31TcZ2zchv8DrQdKJpwbtU0k7n3MeSZGbLJI2XFPwjNV7SYueck7TOzLLNrIuknCiWBQDAa+r92+mc29+UgaXywVhjH5QGr6+x1lnf92/q9w7/rDVJ5e9YquB/nNiiSfC6SdoT9Hqvqp5hjFSmW5TLSpLM7E5JdwZelprZtihii6S9pEP1XDaekjHuZIxZSs64iTl2kjHuuMRsv23Q4hUx92yUYBJPQ347QxK8CL+Pnyv5vqONrV7f+QZ+Z5t8fXXU3n4bu+9BnD9rdZKxvm5sbIME3AaNuL/U6zcymgTPIkxzUQkWjocAACAASURBVJaJZln/ROfmS5ofRTw1MrP1zrm8hq4n1pIx7mSMWUrOuIk5dpIxbmJOSA357QydEPb7mALbrlZsA7aBxDaQ2AYS2yCSaBK8vZLODXrdXdKnUZZpEcWyAAB4TUN+OwEAqLdobnT+rqRcMzvPzFpIuknSirAyKyR91/yGSjoauIYgmmUBAPCahvx2AgBQb7W24DnnzpjZVEkvyz/U80LnXJGZTQ7MnytplfzDPO+Uf6jn22tatkk+yZca3M0zTpIx7mSMWUrOuIk5dpIxbmJOMA357YyCp7ddlNgGbAOJbSCxDSS2QRXmH7wLAAAAAJDsoumiCQAAAABIAiR4AAAAAOARSZvgmdkoM9tmZjvNbFqE+WZm/xGYv9nMLo5HnGExnWtmr5nZVjMrMrOfRCgzzMyOmtmmwOPX8Yg1LKZdZvZBIJ71EeYn1LY2s75B22+TmR0zs3vCyiTEdjazhWb2mZltCZrWzsxeMbMdgb/nVLNsjftAjGP+dzP7KPD//6uZZVezbI3fpaZUTdwPmdm+oO/BmGqWTaRt/eegeHeZ2aZqlo3Ltq6unkv073UyYNvEtw6Jp4b8VnhFQ+pwL2ho3eoFNWyDlPkeRM05l3QP+S9Y/19JveS/FcP7kvqFlRkjabX89xkaKuntBIi7i6SLA88zJW2PEPcwSSvjHWtYTLskta9hfsJt67Dvyj8k9UzE7Szpm5IulrQlaNrvJE0LPJ8m6bfVfK4a94EYx3y1pOaB57+NFHM036U4xP2QpJ9H8R1KmG0dNn+mpF8n0raurp5L9O91oj/YNpXbIW51SJw/d71+K7z0qG8d7pVHQ+pWrzxq2AYp8z2I9pGsLXiXStrpnPvYOXdK0jJJ48PKjJe02Pmtk5RtZl1iHWgw59x+59yGwPMSSVsldYtnTI0k4bZ1kJGS/tc5VxzvQCJxzr0h6XDY5PGSngk8f0bStREWjWYfaBKRYnbOrXHOnQm8XCf//bwSSjXbOhoJta0rmJlJmiBpaSxiiVYN9VxCf6+TANsmhTXgt8IzGlCHe0ID61ZP8PBxdKNL1gSvm6Q9Qa/3quo/OJoycWNmOZIGS3o7wuzLzex9M1ttZhfGNLDInKQ1Zvaemd0ZYX4ib+ubVP0BcKJt5wqdXOBeWIG/HSOUSeRt/j35W3Qjqe27FA9TA11LF1bTtSVRt/U3JB1wzu2oZn7ct3VYPZfs3+t4Y9v4xf17nUCi2adSQW11uOfUo271nAjH0Sn3PahJsiZ4FmFa+P0eoikTF2aWIek5Sfc4546Fzd4gf3fCgZIek/R8rOOL4Arn3MWSRkuaYmbfDJufkNva/DcXHidpeYTZibid6yJRt/kDks5IWlJNkdq+S7H2pKTekgZJ2i9/l8dwCbmtJd2smlvv4rqta6nnql0swrRE2NaJgG3jl2h1COIrmjrcU+pZt3pKhG2Qct+D2iRrgrdX0rlBr7tL+rQeZWLOzNLk/1Iucc4Vhs93zh1zzpUGnq+SlGZm7WMcZnhMnwb+fibpr/J3FQqWkNta/gOADc65A+EzEnE7BzlQ0cU18PezCGUSbpub2SRJYyV9xzkX8cAziu9STDnnDjjnyp1zZyU9VU08ibitm0u6TtKfqysTz21dTT2XlN/rBMK2UeLVIXEWzT7laVHW4Z7RgLrVMyJtg1T7HkQjWRO8dyXlmtl5gVaamyStCCuzQtJ3zW+opKMVTdjxErhm5mlJW51zs6op0zlQTmZ2qfz/o89jF2WVeNLNLLPiufyDaWwJK5Zw2zqg2haORNvOYVZImhR4PknSCxHKRLMPxIyZjZJ0n6Rxzrnj1ZSJ5rsUU2HXin5bkeNJqG0dcKWkj5xzeyPNjOe2rqGeS7rvdYJJ+W2TiHVInEWzT3lalHW4JzSwbvWE6rZBKn0PohbrUV0a6yH/yI3b5R9V7IHAtMmSJgeem6Q5gfkfSMpLgJi/Ln+Xms2SNgUeY8LiniqpSP4R0tZJ+lqcY+4ViOX9QFzJsq3byJ+wZQVNS7jtLH8Cul/SafnP0H9f0lckvSppR+Bvu0DZrpJWBS1bZR+IY8w75b8+qOJ7PTc85uq+S3GO+0+B7+xm+X8kuyT6tg5M/2PFdzmobEJs6xrquYT+XifDI9W3TbzrkDh/9qh/K7z6qEsd7sVHXetWLz5q2AYp8z2I9mGBDQYAAAAASHLJ2kUTAAAAABCGBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8onm8A4ikffv2LicnJ95hAACa2HvvvXfIOdch3nEkC34fASB11Pc3MiETvJycHK1fvz7eYQAAmpiZFcc7hmTC7yMApI76/kbSRRMAAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyCBA8AAAAAPIIEDwAAAAA8ggQPAAAAADyiebwDAJDacqa9GPJ614z8OEUCAABQP8HHM/E+lmlQgmdmoyT9QZJP0gLn3IwIZYZJmi0pTdIh59y/NOQ9ATQ+55xKSkp07NgxHT9+XOXl5TF776fGdQl5vXXr1pi9NxrO5/OpTZs2atu2rTIzM2Vm8Q4JSBrxrHvhHdTDCFfvBM/MfJLmSLpK0l5J75rZCufch0FlsiU9IWmUc263mXVsaMAAGpdzTp999pnKysrUrl07de7cWT6fL2Y/EKf3Hgl5fUH37Ji8LxrOOafy8nKVlpbq0KFDOnHihDp27MjBBRCFeNe98AbqYUTSkGvwLpW00zn3sXPulKRlksaHlblFUqFzbrckOec+a8D7AWgCJSUlKisrU8+ePZWdna3mzZvzw4ComJmaN2+u7Oxs9ezZU2VlZSopKYl3WEBSoO5FY6AeRiQNSfC6SdoT9HpvYFqwPpLOMbP/MbP3zOy71a3MzO40s/Vmtv7gwYMNCAtAXRw7dkzt2rWTz+eLdygJbfPeI5UPVOXz+dSuXTsdO3Ys3qEASYG6F42NehgVGpLgRTrN5MJeN5d0iaR8SddI+pWZ9Ym0MufcfOdcnnMur0OHDg0IC0BdHD9+XBkZGfEOAx6QkZGh48ePxzsMIClQ96IpUA9DatggK3slnRv0urukTyOUOeScK5NUZmZvSBooaXsD3hfwvFiOLFleXs4ZZDQKn8/HIBFAlKh70RSohyE1rAXvXUm5ZnaembWQdJOkFWFlXpD0DTNrbmZtJF0miSHygBjJmfZi5aMmXPeR2hqr+ynfoy+Z2UIz+8zMttRSboiZlZvZ9bGKDYmDfQaNje8UpAYkeM65M5KmSnpZ/qStwDlXZGaTzWxyoMxWSS9J2izpHflvpVDjjx0AAB7wR0mjaioQGI36t/L/jgIA0CgadB8859wqSavCps0Ne/3vkv69Ie9TH4l0s0EAjS+4tWkAt1ZAgnHOvWFmObUUu1vSc5KGNHlAAICU0ZAumgAAoB7MrJukb0uaG0VZRpkGAEStQS14AJDowq8ro7Wv/tiWjWq2pPucc+W1XTPjnJsvab4k5eXlhY9WDQBACBI8AABiL0/SskBy117SGDM745x7Pr5hAQCSHQkekhbXWQJIVs658yqem9kfJa0kuQMANAauwQOABDR9+nSZmQ4cOBDvUFAPZrZU0luS+prZXjP7fvAo0wASG3UwkhkJHgAkoM2bN6tTp07q1KmTJGnHjh0ys5BHp06dNGbMGG3cuDHO0SKcc+5m51wX51yac667c+5p59zc8JGmA2Vvc879JR5xAoiMOhjJjC6aAJCANm/erIEDB1a+3rRpkyRpypQpGjp0qMrLy/XBBx9ozpw5uvLKK7VlyxZ16dIlXuECgKdQByOZkeABQII5efKkduzYoXHjxlVOqzhDfMcdd4QcdGRkZOjhhx/WSy+9pNtvvz3msQKA11AHI9nRRRMA4mjt2rUaOXKk0tPT1bVrV82YMUNFRUUqLy8POYjYuHGjWrRooX79+oUs36tXL0nS0aNHYxo3klfOtBdDHkAqow6GF9GCBwAN0JB7wxUWFurGG2/U4MGD9cgjj+jUqVOaNWuWVq5cKUkaNGhQZdmNGzeqf//+SktLC1nH6tWrJUlDhw6t70cAgJREHQyvIsGLo5rOnDLsP2IpOEnh5tWxUVxcrEmTJik/P1+FhYVq1szfoWL48OHKy8tTq1at1LdvX0nSq+99pAMHDujioV/XoUOH5JzTvn379NRTT2nZsmWaNGkSBxcAUAd1qYP379+vAwcOaMSIEdTBSAp00QSAOJg+fbpOnDihmTNnVh5YSNIll1yirKws9e/fXz6fT5L0UdFmSdLqF55Thw4d1LFjRw0ePFjLli3Tk08+qUWLFoWs+8iRI5owYYIyMzPVtWtXzZ49O3YfDACSQF3q4Irr75YuXUodjKRACx6QIrgxfOIoLy/Xc889pzFjxqh3794R5wdf+/HRFn+C96sZj+qbl/RXWVmZli9friVLlmjfvn0ys5Dlp06dqi+++EL79u1TcXGxRo4cqb59+2r06NFN+8EAIAnUtQ6uSPDmzZunXr16UQcj4ZHgAYiKVwZjSITkds+ePTp8+LCGDBlSZV5RUZFKS0tDE7yiD+Tz+TT2uht1aa5/GO5x48bpo48+0qOPPqr77rtPGRkZklR54PHee++pbdu2uuiii3THHXdo4cKFHFwAgOpeB2/cuFE+n08TJ05U69atJVEHI7GR4NUg/IA2EQ4MUw2tTkg20VzPWFZWVu3yFV15gi/u/2jLZvU4r7daBQ4sJMnMdNttt+nuu+/WypUrddNNN0mStm/frrNnz6p///6VZQcNGqTCwsL6fSAAccXJtcZX1zp448aNys3NrUzuJOpgJDauwQOAGOvZs6fMTK+88krI9HXr1umZZ56RJA0YMECSf+jtfXuK1eeCC6usZ+zYsZKkF154oXJaaWmpsrKyQsplZ2erpKSkUT8DACSrutbBn3zySUiLXgXqYCQqWvCSkNdbFr1ythKoTkZGhm644QYVFBTouuuu06hRo7Rt2zbNnz9fmZmZysrKqjxA2LRpk5xz6tuvf5X15OTk6IILLtDq1at1+vRppaWlKSMjQ8eOHQspd/ToUWVmZsbkswFAoqtPHRwpwaMORqKiBQ8A4mDevHmaOHGi3njjDf30pz/VO++8o4KCAqWnp0e8uD83QgueJI0ePVpHjx7V66+/Lknq06ePzExFRUWVZTZt2hTSXQgAUl1d6+CKFr1w1MFIRCnRguf1Fi+gKYXfyDtVNPXnzs7O1uLFi6tM3717d8jre+65RyOuv63a9cycOVMzZ86sfJ2enq7rr79eDzzwgP70pz+puLhYCxYsqDKMNwCksrrUwffcc0+166EORiKiBQ8APGbOnDlKS0tTly5ddNVVV2natGmM3gYAMUIdjHhLiRY8IFqM2gkvyM7O1vLly+MdBgCkJOpgxBsJnselSsJCN1w/BqgBAABIbXTRBAAAAACPoAUvSdAy4z20OgIAAKCx0YIHAAAAAB5BC14Ti0UrTTJcZ5dIrVXRbq/GiDmRPjeazsBzz4l3CE3OORfvEAAghJnFO4SYoy5GNEjwAKCB3t/zz4jTB3TPjnodwffdC1+upnmNoanXDwBNgWQHiIwErwlwvVxq4/8PAACAeGnQNXhmNsrMtpnZTjObFmH+MDM7amabAo9fN+T9ACBZTZ8+XWamAwcOxDsULV26VHl5eWrdurVycnK09I/z4x2S55jZQjP7zMy2VDP/O2a2OfB408wGxjpGIBUkct37+OOPxzskeFS9W/DMzCdpjqSrJO2V9K6ZrXDOfRhW9G/OubENiBFoEFrUkAg2b96sTp06qVOnTnGN47777tPMmTM1efJk3XXXXVq6dKlm/Oo+te/QSVfljw/prinRZbMB/ijpcUmLq5n/iaR/cc7908xGS5ov6bIYxQakjESue++++2517txZ119/fVxjg/c0pIvmpZJ2Ouc+liQzWyZpvKTwBA8AUt7mzZs1cGB8G2n+8pe/6He/+53+/Oc/a8KECZKkW2+9Vd26n6s/L35aV+WPj2t8XuKce8PMcmqY/2bQy3WSujd1TEAqSuS6t0ePHnriiSdI8NDoGpLgdZO0J+j1XkU++3i5mb0v6VNJP3fOFUVamZndKelOSerRo0cDwgKaRjKMVtqUVky9osbWnPCWn2B1Wa66so21/ng4efKkduzYoaHDr66Mp6kHSwl/j9OnT+tnP/uZ8vPzKw8wJKl169a6cOBgvf/eO40eD6L2fUmrq5vJ72NqS8Xfm8ZSUfeOGzcubjHUVPcOGTJEb775Zg1LA/XTkGvwIo1NGz6c0QZJPZ1zAyU9Jun56lbmnJvvnMtzzuV16NChAWEBQPxseOct3XHTeF3Wp5u6du2qGTNmqKioSOXl5ep7Qf+Qsnv27NGPf/xjffWrX9WQr3bWiIv76p7vf0dbt26tLLNx40b9fPJtGjYoV0O+2ll5eXlatWpVlff9/PPP9YcZD+vbI4bqsj7d9PULe2rw4MF64oknJEnPP/+8du/erbvvvrvKsmlpLVRaUtLIWwLRMLPh8id491VXht9HoHZr167VyJEjlZ6eXqXuDW/BC657W7Vqpc6dO+vaa6+tUvfecMMN6tixo1q1alVj3Xv//ferX79+Sk9PV3Z2dtR1b4sWLXTs2LFG3hJAw1rw9ko6N+h1d/lb6So5544FPV9lZk+YWXvn3KEGvG+DpXpLDBpfU1/nx3c2Nhra2vd/Vv+X7vvR99T3wgH68bRfq0Mbn2bNmqWVK1dKkvpeeFFl2bfeektjxoxRs2bNNHnyZKWd01n79+7Ri38t0PHjxyVJf/7znzVx4kR99fx++sHUn6qZz6eXC5fqW9/6ltasWaORI0dKkkpLSzV06FD98+gxXTvhO+qec56OHfmnire+rz17/B0tCgoKlJGRoYEDB+rQodAq+NjRI2qTntEk2wTVM7MBkhZIGu2c+zze8QDJqrCwUDfeeKMGDx6sRx55RKdOnQqpewcNGlRZNrzu7dWrl4qLi/Xss89WqXsvuugi/fKXv5TP59PChQurrXtLSkr0ve99T71799bhw4e1bt26qOref/7zn8rIiFz3Ag3RkATvXUm5ZnaepH2SbpJ0S3ABM+ss6YBzzpnZpfK3GPIjFgWvDwyS7AlLIv1/EimWVPbp3t361b0/0jdGXK1ZT/1JzZo104Du2Ro+fLjy8vLUqlUr5fTOlSR9fvAzTRg3Tjk5OXr55ZfVsWPHykTqrnt+ITNT4atvaeJ3v6uvj7hav5/7RzVv7q+u75t6h3JycjRjxozKg4yCggLt3LlTf37pDZ0flEQGd9F8/fXXVVpaqi5dukSMv88FF4a8Hj/sUv2/j8zSkMu/3ngbCZXMrIekQkkTnXPb4x0PkKyKi4s1adIk5efnq7CwUM2a+TunBde9ffv2lSQdOHBA48Lq3gq//vWvZWYqKirSd7/7XeXn52v58uWVde+tt95abd27cePGkCQyWG1174ABAxptWwAV6p3gOefOmNlUSS9L8kla6JwrMrPJgflzJV0v6YdmdkbSCUk3Oe5KCcCDFjw2S1+cPKGf/erfKg8wJOmSSy5RVlaWcnNz5fP5JElPzpqhI0eOaO3atSEHGJIqDyYWPDZTZqZfPTKrcpoknXPOORowYIC2b/8yJ/jnP/03Wv9g4/qQBK/Cvn37dPDgQU2aNEm33npryLwDBw7o1ltvDWld/OLkSe3Z9bH6nH9h+KoQJTNbKmmYpPZmtlfSg5LSpMrfx19L+oqkJ8xMks445/LiEy2QvKZPn64TJ05o5syZtda9Dz30UK11b8VtFebOnRt13fv2229HTPCiqXurSwyBhmjQjc6dc6skrQqbNjfo+ePyDxONgGRsbUn21javaervkNe75AV/vsYa6KS8vFz/Z9UKfX34VTo357yI8yuuASkvL9ealc/rmmuuqTyrHO7MmTP6n1de0jdHXqOvdOhYZb5zTunp6ZWvb7nlFs2bN0//dv9PtXDOoxpxzVhd/a1rNaD71ZKkTz75RJJ0+eWX68orrwzZBu999BdJUt7QL1vqdm7bqvYdOynrnHPquikQ4Jy7uZb5P5D0gxiFA3hSeXm5nnvuOY0ZM0a9e/eOOD+47i0oKKi17l2xYoXGjh0b8bYK1dW9kydP1iOPPKJvf/vbmjBhgi6//HJJVeveYEuXLpUkDRs2rO4fHKhFg250DgCQ/vHpPh098k9dOPDiKvOKiopUWlpaeZBRUfbSSy+tdn17d+/SieNlOu+rfarMO3v2rLZu3RoyaECXLl20ZcsWzV6wRJf/ywi9+PxyfffaazR16lRJ/utEJCkzM7PK+l59aaWap6Vp+NVjKqdt27qlSpdNAEg0e/bs0eHDhzVkyJAq88Lr3oqyNdW9H3/8scrKynTBBRdUmVdT3fv88/6TdkuWLNHXvva1qOrewsJCpaWlafx4/+1pTp8+rczMTE2bNi2k3A033KCHHnqoli0BhCLBSyE5014MecTy/eANm/ceCXnA78TxsmrnzZ49W9KXF/mXldY+Ylqgy55atGhZZd7KlSt1+PBhXXvttSHTW7RooeHXjNGvZ8zW6rfe14UDL9aCBQskfXlwURI2UuY/Pt2r19as0tjrbgxprduxtUh9wkb8BIBEU1YWfd0bzWiVFXVvy5Z1q3vHjx+vefPmadeuXRoyZEitde+ePXv0wgsvaOLEiWrXrp0kKS0tTatWrdKiRYsqy7311ltau3atfv7zn9caOxCMBA8AGqhr93NlZlr3t9dCpq9bt07PPPOMpC8vpO/c9Vz5fD69/PLLVdZz+tQpSVK3c3uqVavW2vjuupD5n/1jv+6++25dcMEF+td//VdJqjIqmxRIDJ1Tt27dJEn9+/dXWlqa1qxZU1mmvLxc/9+0e9W6TRv96GehZ4y3f7hFubTgAUhwPXv2lJnplVdeCZkeqe7t0aNHtXXvqUDde95556l169Zau3ZtyPxPP/00qrq3ZcuWclHUvXfddZfS09P18MMPhyw/dOhQHT58WEePHpUk/exnP9PDDz/MSJuoswZdgwc0Blr40Nhi3brYJj1DV+VfqzUr/6p775ioK4aN1K7/3am/Ln1GmZmZysrKUlZWllRyRG2zsjR+wndUuHSxLr3imxp+Tb58zZpp+0dF+nj7Nv2xcLWaN2+u2374Y8199LeaNvUHGvK1b+jA/k+1/NlFapnWXCtWrKi8+P/ee+/VO++8o/Hjx6tluy764uRJrVn5V23d8r6WLVsmScrKytKtt96qRYsW6Uc/+pG+0iNXLxYW6MPNm/SHhf+pTl26hXyeHR99SBdNAAkvIyNDN9xwgwoKCnTddddp1KhR2rZtm+bPnx9a90rKzs7W7bffrgULFmjEiBG69tpr5fP5tHnzZn344Yf629/+pubNm+sXv/iFHn74Yd1yyy0aPny49u7dq7lz58rn81Vb9/bu3VsnTpxQQUGBNmzYUG3dO3DgQD377LNav369VqxYoe7du4d8nrS0NPXo0UPbt2/Xrl27dOTIEX3/+9+P7UaFJ5DgAUAj+NWMR9WiZQut/e9X9NYbr+n8Cy9SQUGB7rrrrio32f3l9N+rV25fvVCwRP8x4zdq0aKFeuX20Q0Tv1dZ5s6f/D8yM/112Z/0yqoV6tips64ee60e+/f/X+3bt68sd+WVV+rgwYN69tlndejzz9XuK+3V76JBWvz8Gt3wrZGV5R5//HE1a9ZMS5cu1Znycl182df07H+9UqUr5oH9+3TyxInKWzogtYSfcGNgLSS6efPmqWXLllq1apXWrFmjwYMHV1v3PvHEE+rXr58WLVqk+++/Xy1bttQFF1ygH/7wh5VlfvWrX8nM9PTTT+svf/mLunbtqgkTJujBBx+stu79/PPP1bFjR11yySV66623Qq7zC657lyxZom9+85t6++23q709Qm5urrZs2aLp06dr9uzZlSOAAnVBggdJjTNSZjK2xCVSzIkUC+qubVaWps+eGzJtQPds7d69u0rZtLQ0TbzjR5p4x4+qXZ/P59Pke+/T5HvvC5nevn3oyJ+TJk3SpEmTJNXcctmmTRstWLBACxYsqLHctg+3qFefvhxUpBDqHiSz7OxsLV68uMr06uree++9V/fee2+16/P5fHrwwQf14IMP1vi+wXVvTYLr3mh89atf1UMPPaTevXtr7NixUS0DhCPBAwBU2rG1SLnc/w4A4iI3N1f/l717j7arru+9//6ySapJdhJpwi0QQBrRSLlukOOjBWppE7BGPaKIg4tHjfQRL8PaQZ7HUS/DwTDYQ0t5jIZIo2I7TPEQMWoQLedUSm3aBIGQgIGIhGxCSTBALqAhm+/zx15JV3b2Ze297nO9X2PskTUva83vmntm/9Zn/eb8zc2bN7N8+fJml6I21vEBr1ano/gNqKRWM5ZrET9w9SfrUIkkqRLd3d2ce+65nHnmmc0uRW2s4wOeDjaasGqwlRprYGir1c3aJUnN9+CDD3LGGQffU1UaDW+TIEmSJLWAtWvXDjkAi1Qpe/CkBmqVHs/yOr72tqM4oruJxQxQ3kNVae+UN12XJBXBXXfd1ewSVKFW+Uw3GHvwJEmSJKkg7MEbYLg0XotbCbS7Vv62Qq3L68YkSZIawx48SZIkSSoIe/AkjdlYr31rlWvmWqUOSZKk01gTzgAAIABJREFUWrEHTxKZ2ewSVAAeR5IkNZ89eGPUqteitWpdal2/2ZuQL0N0NbsUjUEr9UL29fXR1eVxJElSM9mDJ3W4h7b+htzzm2aXoQLYtWsXEyZMaHYZkiR1NHvwVEj2ZFbuXza9wNnHbeewI15JHOJ3Phqbvr4+tm/fzrRp05pdiiRJHc1Pc1KHW7dtD99/aDvbn36Sl3+zm3y5z2upVJHMZO/evTz33HNs2rSJiRMn0t3d3eyypLbh31rVmseUwB48tRF75epn+cO7eOSZPbz5uB3MPvwVvOLQaHZJADy885X7Hz/97ItNrORAldZVvt5I67ajcTtfSVdXFxMmTGDatGl0d3cT0RrHjtTqurq66Ovr49BD/Sim2vFaaIEBr6MZmFRu3bY9rNu2p9llHODxhRftfzy3DY/X8vqhPd/DcAa+P0mVmzBhArt27WLq1KnNLkUF4rXQAk/RlCSp5iJiaURsjYh1QyyPiLgxIjZGxNqIOKPRNaq5Jk+ezPbt2+nr62t2KSqIfddCT548udmlqMkMeJIk1d43gDnDLJ8LzCr9zAe+2oCa1EK6u7uZOHEimzZt4rnnnmPv3r1eP6VR81poDcZTNCVJqrHMvDsijh9mlXnALdn/iX5VREyNiKMy86mGFKimiwgOP/xwdu7cyY4dO9i6dau9eRoTr4XWQAY8SZIabwawuWy6tzTvoIAXEfPp7+Vj5syZDSlOjRERTJ482VPqJNWUp2hKktR4g33FPuj5eZm5JDN7MrNn+vTpdS5LktTuDHiSJDVeL3Bs2fQxwJYm1SJJKpCqAl5EzImIDaVRwBYMs95ZEdEXEe+qZnuSJBXECuDy0mia5wDPe/2dJKkWxnwNXkR0AYuAC+j/JnJ1RKzIzIcGWe864M5qCpUkqV1ExLeB84BpEdELfBYYB5CZi4GVwIXARuAF4P3NqVSSVDTVDLJyNrAxMx8DiIhl9I8K9tCA9T4K3AacVcW2JElqG5n53hGWJ/CRBpUjSeog1ZyiOdQIYPtFxAzgHcDikV4sIuZHxJqIWLNt27YqypIkSZKkzlRNwKtkBLAbgGsyc8QbuzhKmCRJkiRVp5pTNCsZAawHWFa64eI04MKI2JuZt1exXUmSJEnSIKoJeKuBWRFxAvAkcAlwafkKmXnCvscR8Q3gB4Y7SZIkSaqPMQe8zNwbEVfTPzpmF7A0M9dHxFWl5SNedydJkiRJqp1qevDIzJX0D/VcPm/QYJeZV1azLUmd5/gFP2x2CZIkSW2lqhudS5IkSZJaR1U9eJKkodkDqXZXfgw/vvCiJlYiSaqUAU+SJEmSRtAuX9x6iqYkSZIkFYQBT5IkSZIKwoAnSZIkSQVhwJMkSZKkgjDgSZIkSVJBGPAkSZIkqSAMeJIkSZJUEN4HT5IkjZo3QZek1mTAkyRJI2qXG/xKUqfzFE1JkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJUh1ExJyI2BARGyNiwSDLp0TE9yPigYhYHxHvb0adkqRiMeBJklRjEdEFLALmArOB90bE7AGrfQR4KDNPBc4Dro+I8Q0tVJJUOAY8SZJq72xgY2Y+lpl7gGXAvAHrJNAdEQFMArYDextbpiSpaLxNgiRJtTcD2Fw23Qu8YcA6XwZWAFuAbuA9mfnywBeKiPnAfICZM2fWpdhqDbyFgvfFk6TmMeBJklR7Mci8HDD9J8D9wB8CJwI/iYh/ycwdBzwpcwmwBKCnp2fga0iSaqj8C6t2/bLKUzQlSaq9XuDYsulj6O+pK/d+YHn22wj8Cnhtg+qTJBWUAU+SpNpbDcyKiBNKA6dcQv/pmOWeAN4CEBFHACcBjzW0SklS4XiKpiRJNZaZeyPiauBOoAtYmpnrI+Kq0vLFwBeAb0TEg/Sf0nlNZj7TtKIlSYVgwJMkqQ4ycyWwcsC8xWWPtwB/3Oi6JEmVGTiAVLsw4EmSJEkqlE4e3ddr8CRJkiSpIAx4kiRJklQQnqIpSZIkqdCKcH+7SlXVgxcRcyJiQ0RsjIgFgyyfFxFrI+L+iFgTEW+qZnuSJEmSpKGNuQcvIrqARcAF9N/QdXVErMjMh8pWuwtYkZkZEacAt+JNXCVJkiSpLqrpwTsb2JiZj2XmHmAZMK98hczclZlZmpwIJJIkSZKkuqjmGrwZwOay6V7gDQNXioh3AF8EDgeGPOE1IuYD8wFmzpxZRVmSJEmSVJl2vd/dUKrpwYtB5h3UQ5eZ383M1wJvB74w1Itl5pLM7MnMnunTp1dRliRJkiR1pmoCXi9wbNn0McCWoVbOzLuBEyNiWhXblCRJkiQNoZpTNFcDsyLiBOBJ4BLg0vIVIuL3gF+WBlk5AxgP/LqKbUqSpBbXScORS1KrGXPAy8y9EXE1cCfQBSzNzPURcVVp+WLgvwOXR8RLwIvAe8oGXZEkSZIk1VBVNzrPzJXAygHzFpc9vg64rpptSJIkSZIqU9WNziVJkiRJrcOAJ0mSJEkFYcCTJEmSpIIw4EmSJElSQVQ1yIokSZIktZPyW7kUkT14kiRJklQQ9uBJkiRJantF75mrlD14kiRJklQQBjxJkiRJKggDniRJdRARcyJiQ0RsjIgFQ6xzXkTcHxHrI+Knja5RklQ8XoMnSVKNRUQXsAi4AOgFVkfEisx8qGydqcBXgDmZ+UREHN6caltT+bU0jy+8qImVSFJ7sQdPkqTaOxvYmJmPZeYeYBkwb8A6lwLLM/MJgMzc2uAaJUkFZMCTJKn2ZgCby6Z7S/PKvQZ4VUT8c0TcGxGXD/ZCETE/ItZExJpt27bVqVxJUlF4iqYkSbUXg8zLAdOHAmcCbwFeCfxbRKzKzEcOeFLmEmAJQE9Pz8DXkKSO4unbIzPgSZJUe73AsWXTxwBbBlnnmczcDeyOiLuBU4FHkCSNyPveDc6AJ0lS7a0GZkXECcCTwCX0X3NX7nvAlyPiUGA88AbgbxpaZRP47bsk1ZcBT5KkGsvMvRFxNXAn0AUszcz1EXFVafnizHw4In4ErAVeBm7OzHXNq1qSVAQGPEmS6iAzVwIrB8xbPGD6r4C/amRdjTbcKVT25kkajH8bqmPAkyRJTee1NFLn8v9/bXmbBEmSJEkqCHvwJElSx/OUMElFYcCTJEmS1JI8fXP0PEVTkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkE4yIokSVKFHG1TncJjvX1VFfAiYg7wt0AXcHNmLhyw/H3ANaXJXcCfZeYD1WxTkiRpn4Ej7PlBVFKnG/MpmhHRBSwC5gKzgfdGxOwBq/0KODczTwG+ACwZ6/YkSZIkScOrpgfvbGBjZj4GEBHLgHnAQ/tWyMyfla2/Cjimiu1JkqQOZC+dJFWumkFWZgCby6Z7S/OG8gHgjqEWRsT8iFgTEWu2bdtWRVmSJEmS1Jmq6cGLQebloCtGnE9/wHvTUC+WmUsoncLZ09Mz6OtIkiQN7NEbapk9fVLrGu7/sapTTcDrBY4tmz4G2DJwpYg4BbgZmJuZv65ie5IkSZKkYVQT8FYDsyLiBOBJ4BLg0vIVImImsBy4LDMfqWJbkiRJNWPvgaSiGnPAy8y9EXE1cCf9t0lYmpnrI+Kq0vLFwGeA3wW+EhEAezOzp/qyJUmSpGJpxwGF2rHmoqvqPniZuRJYOWDe4rLHHwQ+WM02JEmSWpEfbKXheU1sc1QV8CRJkiRpJJ4W3TgGPEmSJKkg7FmWAU+SJKlF+OFc7cDeuNZmwJMkSR2hmR9KvRZJUqMc0uwCJEkqooiYExEbImJjRCwYZr2zIqIvIt7VyPqkTnP8gh/u/5GKzB48SZJqLCK6gEXABUAvsDoiVmTmQ4Osdx39txxSizIQHMweSal1GfAkSaq9s4GNmfkYQEQsA+YBDw1Y76PAbcBZjS1PraRVA6QhTmpPBjxJkmpvBrC5bLoXeEP5ChExA3gH8IcME/AiYj4wH2DmzJk1L7TIWuWau3q95lhCl4O4aCxGczy36hcWncSAJ0lS7cUg83LA9A3ANZnZFzHY6qUnZS4BlgD09PQMfA21IT8sS6onA54kSbXXCxxbNn0MsGXAOj3AslK4mwZcGBF7M/P2xpQotT97JKvjFwjFZMCTJKn2VgOzIuIE4EngEuDS8hUy84R9jyPiG8APDHftrVU/LLdqXaqdSn/HHgudwYAnSVKNZebeiLia/tExu4Clmbk+Iq4qLV/c1AIljUktBp4ZTa/jcNsz1GkoBjxJkuogM1cCKwfMGzTYZeaVjahJlfEDcWerx+ihHlNqJAOeJEmSWlYjrrNrhwDWDjWqNRjwJEmS1FQOliLVjgFPkiRJNVPvG6S3ak+WN4ZXqzDgSZIktaFWCRStGrjqrVPft1qfAU+SJKnNtUrYG0kje/cqHZ2yHryZvZrJgCdJkqRhdWoI6dT3rfZmwJMkSVJdjDUgGayksTPgSZIkqXChqmjvR6qUAU+SJEkNZwCT6uOQZhcgSZIkSaoNA54kSZIkFYSnaEqSJBWIQ/RLnc0ePEmSJEkqCAOeJEmSJBWEAU+SJEmSCqKqgBcRcyJiQ0RsjIgFgyx/bUT8W0T8NiI+Vc22JEmSJEnDG/MgKxHRBSwCLgB6gdURsSIzHypbbTvwMeDtVVUpSZIkSRpRNT14ZwMbM/OxzNwDLAPmla+QmVszczXwUhXbkSRJkiRVoJqANwPYXDbdW5o3JhExPyLWRMSabdu2VVGWJEmSJHWmagJeDDIvx/pimbkkM3sys2f69OlVlCVJkiRJnamagNcLHFs2fQywpbpyJEmSJEljVU3AWw3MiogTImI8cAmwojZlSZLU3ioYafp9EbG29POziDi1GXVKkoplzKNoZubeiLgauBPoApZm5vqIuKq0fHFEHAmsASYDL0fEJ4DZmbmjBrVLktSSKhxp+lfAuZn5bETMBZYAb2h8tZKkIhlzwAPIzJXAygHzFpc9/k/6T92UJKmT7B9pGiAi9o00vT/gZebPytZfhe2lJKkGqrrRuSRJGtRoR5r+AHDHYAscZVqSNBoGPEmSaq/ikaYj4nz6A941gy13lGlJ0mhUdYqmJEkaVEUjTUfEKcDNwNzM/HWDapMkFZg9eJIk1d6II01HxExgOXBZZj7ShBolSQVkD54kSTVWyUjTwGeA3wW+EhEAezOzp1k1S5KKwYAnSVIdVDDS9AeBDza6LklSsXmKpiRJkiQVhAFPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIKoKeBExJyI2RMTGiFgwyPKIiBtLy9dGxBnVbE+SpHZhGylJaoYxB7yI6AIWAXOB2cB7I2L2gNXmArNKP/OBr451e5IktQvbSElSs1TTg3c2sDEzH8vMPcAyYN6AdeYBt2S/VcDUiDiqim1KktQObCMlSU1xaBXPnQFsLpvuBd5QwTozgKcGvlhEzKf/G0yAXRGxoYraAKYBz1T5Go1irfXTTvVaa320U63QRvXGdTWp9bha1NKCatZGdnj7CO1Vr7XWRzvVCu1Vr7XWQY3aRxhjG1lNwItB5uUY1umfmbkEWFJFPQduOGJNZvbU6vXqyVrrp53qtdb6aKdaob3qbadam6BmbWQnt4/QXvVaa320U63QXvVaa300u9ZqTtHsBY4tmz4G2DKGdSRJKhrbSElSU1QT8FYDsyLihIgYD1wCrBiwzgrg8tJIYecAz2fmQadnSpJUMLaRkqSmGPMpmpm5NyKuBu4EuoClmbk+Iq4qLV8MrAQuBDYCLwDvr77kitXsdJYGsNb6aad6rbU+2qlWaK9626nWhmrxNrLdfm/tVK+11kc71QrtVa+11kdTa43MQS+JkyRJkiS1mapudC5JkiRJah0GPEmSJEkqiLYOeBFxcUSsj4iXI2LIoUgjYk5EbIiIjRGxoGz+YRHxk4h4tPTvq+pY64jbioiTIuL+sp8dEfGJ0rLPRcSTZcsubGatpfUej4gHS/WsGe3zG1VrRBwbEf8nIh4uHS8fL1tW9/061PFXtjwi4sbS8rURcUalz21Cre8r1bg2In4WEaeWLRv0eGhyvedFxPNlv9/PVPrcJtT6F2V1rouIvog4rLSsofs2IpZGxNaIWDfE8pY5ZjW4sH2sm0r3zVD/b1tw39pG1q7WlmkjK6jV9nFstbZH+5iZbfsDvA44CfhnoGeIdbqAXwKvBsYDDwCzS8u+BCwoPV4AXFfHWke1rVLd/wkcV5r+HPCpBu3XimoFHgemVfte610rcBRwRulxN/BI2TFQ1/063PFXts6FwB303xPrHODfK31uE2p9I/Cq0uO5+2od7nhocr3nAT8Yy3MbXeuA9f8U+N9N3Ld/AJwBrBtieUscs/4M+zu0fWxyvUP9v221fYttZC1rbYk2ssJaz8P2cSz1tkX72NY9eJn5cGZuGGG1s4GNmflYZu4BlgHzSsvmAd8sPf4m8Pb6VDqmbb0F+GVmbqpjTUOpdr+01H7NzKcy8+elxzuBh4EZdayp3HDH3z7zgFuy3ypgakQcVeFzG1prZv4sM58tTa6i/75dzVLN/mm5fTvAe4Fv17GeYWXm3cD2YVZplWNWQ7B9rCvbyNqxjawP28c6aZf2sa0DXoVmAJvLpnv5rz9cR2TpnkOlfw+vYx2j3dYlHHwAX13q7l1az1M6qLzWBH4cEfdGxPwxPL+RtQIQEccDpwP/Xja7nvt1uONvpHUqeW4tjXZ7H6D/W6p9hjoe6qXSev9bRDwQEXdExOtH+dxaqXh7ETEBmAPcVja70ft2JK1yzKo6to9jYxtZO7aR9WH72DwtcbyO+T54jRIR/wQcOciiT2fm9yp5iUHm1eXeEMPVOsrXGQ+8Dfh/ymZ/FfgC/bV/Abge+B9jq7Rmtf5fmbklIg4HfhIRvyh9s1FTNdyvk+j/o/CJzNxRml3T/TrYZgeZN/D4G2qdhh27I9Rx8IoR59PfeL2pbHZDjofyMgaZN7Den9N/Gteu0rUjtwOzKnxuLY1me38K/Gtmln9D2Oh9O5JWOWY7mu3jfjX/O24baRs5iHZqI20fm6cljteWD3iZ+UdVvkQvcGzZ9DHAltLjpyPiqMx8qtR9urWaDQ1Xa0SMZltzgZ9n5tNlr73/cUR8DfhBs2vNzC2lf7dGxHfp736+mxbcrxExjv6G6x8yc3nZa9d0vw5iuONvpHXGV/DcWqqkViLiFOBmYG5m/nrf/GGOh6bVW/YhhcxcGRFfiYhplTy30bWWOah3ogn7diStcsx2NNvH/a9d87/jtpH7X9s2cuQ6DtAibaTtY/O0xPHaCadorgZmRcQJpW/+LgFWlJatAK4oPb4CqOQbz7EazbYOOr+49Id5n3cAg47eUyMj1hoREyOie99j4I/Lamqp/RoRAfwd8HBm/vWAZfXer8Mdf/usAC6PfucAz5dOpankuQ2tNSJmAsuByzLzkbL5wx0Pzaz3yNLvn4g4m/6/eb+u5LmNrrVU4xTgXMqO4ybt25G0yjGr6tg+jo1tZO3YRjavVtvH+miN4zUbNOpMPX7o/2PTC/wWeBq4szT/aGBl2XoX0j8q1C/pP3Vl3/zfBe4CHi39e1gdax10W4PUOoH+/2BTBjz/W8CDwNrSAXFUM2ulfxSgB0o/61t5v9J/ikSW9t39pZ8LG7VfBzv+gKuAq0qPA1hUWv4gZSPeDXXs1nF/jlTrzcCzZftxzUjHQ5PrvbpUzwP0X/D+xlbdt6XpK4FlA57X8H1L/wfop4CX6P8b+4FWPWb9GfJ3aPvYxHqH+3/bavsW28ha1toybWQFtdo+jq3Wtmgfo7RBSZIkSVKb64RTNCVJkiSpIxjwJEmSJKkgDHiSJEmSVBAGPEmSJEkqCAOeJEmSJBWEAU+SJEmSCsKAJ0mSJEkFYcCTJEmSpIIw4EmSJElSQRjwJEmSJKkgDHiSJEmSVBAGPEmSJEkqCAOeJEmSJBWEAU+SJEmSCsKAJ0mSJEkFYcCTJEmSpIIw4EmSJElSQRjwJEmSJKkgDHiSJEmSVBAGPEmSJEkqCAOeJEmSJBWEAU+SJEmSCsKAJ0mSJEkFYcCTJEmSpIIw4EmSJElSQRjwJEmSJKkgDHiSJEmSVBAGPEmSaiwilkbE1ohYN8TyiIgbI2JjRKyNiDMaXaMkqZgMeJIk1d43gDnDLJ8LzCr9zAe+2oCaJEkdoKqAFxFzImJD6RvIBYMsPy8ino+I+0s/n6lme5IktYPMvBvYPswq84Bbst8qYGpEHNWY6iRJRXboWJ8YEV3AIuACoBdYHRErMvOhAav+S2a+dTSvPW3atDz++OPHWpokqU3ce++9z2Tm9GbX0QQzgM1l072leU8NXDEi5tPfy8fEiRPPfO1rX9uQAiVJzTXWNnLMAQ84G9iYmY8BRMQy+r+RHBjwRu34449nzZo11b6MJKnFRcSmZtfQJDHIvBxsxcxcAiwB6OnpSdtHSeoMY20jqzlFc6hvHwf6bxHxQETcERGvH+rFImJ+RKyJiDXbtm2roixJklpeL3Bs2fQxwJYm1SJJKpBqAl4l3z7+HDguM08F/j/g9qFeLDOXZGZPZvZMn96JZ+tIkjrICuDy0mia5wDPZ+ZBp2dKkjRa1ZyiOeK3j5m5o+zxyoj4SkRMy8xnqtiuJEktLSK+DZwHTIuIXuCzwDiAzFwMrAQuBDYCLwDvb06lkqSiqSbgrQZmRcQJwJPAJcCl5StExJHA05mZEXE2/T2Gv65im5IktbzMfO8IyxP4SIPKkSR1kDEHvMzcGxFXA3cCXcDSzFwfEVeVli8G3gX8WUTsBV4ELik1apIkSZKkGqumB4/MXEn/aSbl8xaXPf4y8OVqtiFJkiRJqkxVNzqXJEmSJLUOA54kSZIkFYQBT5IkSZIKwoAnSZIkSQVR1SArkmrn+AU/3P/48YUXNbESSZIktSsDnqS6y0x27tzJjh07eOGFF+jr62t2SaqRrq4uJkyYwOTJk+nu7iYiml2SJEkdzYAnqa4yk61bt7J7924OO+wwjjzySLq6ugwCBZCZ9PX1sWvXLp555hlefPFFDj/8cH+3kiQ1kQFPUl3t3LmT3bt3c9xxx9HV1dXsclRDEcGhhx7K1KlT6e7uZtOmTezcuZPJkyc3uzRJkjqWg6xIqqsdO3Zw2GGHGe4Krquri8MOO4wdO3Y0uxRJkjqaAU9SXb3wwgtMmjSp2WWoASZNmsQLL7zQ7DIkSepoBjxJddXX12fvXYfo6upyAB1JkprMgCep7hx0ozP4e5YkqfkMeJIkSZJUEAY8SZIkSSqIwt4m4fgFP9z/+PGFFzWxEkmSJElqDHvwJEmSJKkgDHiSJEmSVBAGPEkqgGuvvZaI4Omnn252KZIkqYkMeJJUAGvXruWII47giCOOAODRRx8lIg74OeKII7jwwgu57777mlytJEmql8IOsiJJnWTt2rWceuqp+6fvv/9+AD7ykY9wzjnn0NfXx4MPPsiiRYv4oz/6I9atW8dRRx3VrHIlSVKdGPAkqc395je/4dFHH+Vtb3vb/nn7euk+9KEPHRD8Jk2axOc//3l+9KMf8f73v7/htUqSpPryFE1JaiP33HMPb3nLW5g4cSJHH300CxcuZP369fT19R0Q5O677z7Gjx/P7NmzD3j+q1/9agCef/75htYtSZIawx48SWoTy5cv5z3veQ+nn346X/ziF9mzZw9//dd/zQ9+8AMATjvttP3r3nfffZx88smMGzfugNe44447ADjnnHMaV7gkSWoYA54ktYFNmzZxxRVXcNFFF7F8+XIOOaT/BIzzzz+fnp4eXvGKV3DSSScB8NRTT/H000/zh3/4hzzzzDNkJk8++SRf+9rXWLZsGVdccYUBT5KkgvIUTUlqA9deey0vvvgi119//f5wB3DmmWcyZcoUTj75ZLq6uoD/uv7u29/+NtOnT+fwww/n9NNPZ9myZXz1q1/l61//+gGv/dxzz/Hud7+b7u5ujj76aG644YbGvTFJklRT9uBJUovr6+vjtttu48ILL+TEE08cdPnA6+8AbrrpJl796leze/duvvOd7/AP//APPPnkk0TEAc+/+uqr+e1vf8uTTz7Jpk2beMtb3sJJJ53E3Llz6/vGJElSzRnwJKnFbd68me3bt3PWWWcdtGz9+vXs2rXroIDX1dXFZZddxitf+UoA3va2t/GLX/yCv/mbv+Gaa65h0qRJAPvD37333svkyZP5/d//fT70oQ+xdOlSA54kSW3IUzQlqcXt3r17yGX7TqccOMDKrFmz9oc7gIjgyiuvZPfu3fsHZQF45JFHePnllzn55JP3zzvttNNYt25dLd+CJElqEHvwJLWE4xf8sNkl1MTjCy+q+Wsed9xxRAQ/+clP+Mu//Mv981etWsU3v/lNAE455RSg//YHv/rVr3j3u9990Ou89a1v5aMf/Sjf+973uOSSSwDYtWsXU6ZMOWC9qVOnsnPnzpq/D0mSVH8GPElqcZMmTeLiiy/m1ltv5Z3vfCdz5sxhw4YNLFmyhO7ubqZMmbI/pN1///1k5gGnbO5z/PHH87rXvY477riDl156iXHjxjFp0iR27NhxwHrPP/883d3dDXlvkiSptjxFU5LawE033cRll13G3XffzSc/+Un+4z/+g1tvvZWJEycOOsDKvh69gebOncvzzz/PT3/6UwBe85rXEBGsX79+/zr333//AadsSpKk9mEPniS1galTp3LLLbccNP+JJ544YPoTn/gEn/jEJ4Z8neuvv57rr79+//TEiRN517vexac//Wm+9a1vsWnTJm6++eaSKeJlAAAa7UlEQVSDbqUAsLb3uQOmTzlm6mjfhiRJqjN78CSpwy1atIhx48Zx1FFHccEFF7BgwQJH0JQkqU3ZgydJHW7q1Kl85zvfaXYZkiSpBgx4UpMUZdRISZIktQ5P0ZQkSZKkgjDgSZIkSVJBeIqm1ECelil1joiYA/wt0AXcnJkLByyfAvw9MJP+9vh/ZubBw5dKkjQKVfXgRcSciNgQERsjYsEw650VEX0R8a5qtidJnSgimv6j0YmILmARMBeYDbw3ImYPWO0jwEOZeSpwHnB9RIxvaKGSpMIZcw9eWeN1AdALrI6IFZn50CDrXQfcWU2hktSpMrPZJWj0zgY2ZuZjABGxDJgHlLeRCXRHf4KeBGwH9ja6UElSsVTTg7e/8crMPcC+xmugjwK3AVur2JYkSe1kBrC5bLq3NK/cl4HXAVuAB4GPZ+bLA18oIuZHxJqIWLNt27Z61StJKohqAt6IjVdEzADeASwe6cVswCSpOtdeey0RwdNPPz3scv/GNsRg57UO7Ir9E+B+4GjgNODLETH5oCdlLsnMnszsmT59eu0rlSQVSjUBr5LG6wbgmszsG+nFbMAkqTpr167liCOO4Igjjhh0+bp16zj88MPxb2xD9ALHlk0fQ39PXbn3A8uz30bgV8BrG1SfJKmgqhlFs5LGqwdYVrpAfxpwYUTszczbq9iuJGkQa9eu5dRTTx1y+bp16zj55JMbWFFHWw3MiogTgCeBS4BLB6zzBPAW4F8i4gjgJOCxhlYpSSqcanrw9jdepVG/LgFWlK+QmSdk5vGZeTzwv4D/23AnSbX3m9/8hkcffZTTTjtt0OV79+7lkUceMeA1SGbuBa6mf4Cxh4FbM3N9RFwVEVeVVvsC8MaIeBC4i/4zXp5pTsWSpKIYcw9eZu6NiH2NVxewdF/jVVo+4nV3krTP4wsvanYJbeOee+7hs5/9LKtWrWLKlCl87GMf44ILLqCvr++gHrxVq1bx2c9+lnvuuYc9e/awaNEi/umf/onrrruOt771rU16B50hM1cCKwfMW1z2eAvwx42uS5JUbFXdBy8zV2bmazLzxMy8tjRv8WDhLjOvzMz/Vc32JKnTLV++nPPPP5+dO3fyxS9+kU9+8pPceOONfPzjHwc4oAfvH//xH3nzm9/M7t27mTt3LgCf/OQnefbZZ3nnO9/JunXrmvIeJElS/VQV8CRJjbNp0yauuOIKLrroIlatWsXHPvYxPvWpT/H973+ff/3Xf+UVr3gFJ510EgC//OUvufLKK7ngggv46U9/yowZMxg/fjzXXnstf/d3f8dLL73EN77xjea+IUmSVHPVDLIiSWqga6+9lhdffJHrr7+eQw75r+/nzjzzTKZMmcKsWbPo6uoC4Etf+hJ79uzhpptuoquri/vuu4+TTz6ZcePG8Qd/8AcA9Pb2NuV9SJKk+rEHT5LaQF9fH7fddhsXXnghJ5544qDL911/l5l897vf5Y1vfCPHHnssmckDDzzA6aefDsDu3bsBmDp1auPegCRJaggDntSCjl/wwwN+pM2bN7N9+3bOOuusg5atX7+eXbt27Q94mzdvZtu2bfsD3caNG9mxYwdnnHEGAPfeey/Q3/P30ksv0d3dzYIFCw54zYsvvpjPfe5zdXxHkiSpHgx4ktQG9vW6DeaGG24A/muAlWeffRaACRMmAHDfffcB7A983/72txk3bhx/+qd/yrhx41i5ciVf//rX97/ev/3bv3HPPffwqU99qvZvRJIk1ZUBT5LawHHHHUdE8JOf/OSA+atWreKb3/wmAKeccgoARx99NNB/OwWAn//85xxyyCGccsop3HXXXfz93/89H/7whznyyCMBOOecc9i+fTvPP/88AH/+53/O5z//eSZNmtSQ9yZJkmrHQVYkqQ1MmjSJiy++mFtvvZV3vvOdzJkzhw0bNrBkyRK6u7uZMmUKU6ZMAWD69Om8/e1v5/bbb+fyyy/nwQcf5FWvehWf//znufHGG3nzm9/Ml770pf2vPW7cOGbOnMkjjzzC448/znPPPccHPvCBZr1VSZJUBQOeJLWJm266id/5nd9h5cqV/PjHP+b000/n1ltv5cMf/vBBNzj/+te/zqRJk7j99tvZuXMnhxxyCLfddhuf/vSnueaaaxg/fvwB68+aNYt169Zx7bXXcsMNN+wfjVOSJLUXA54ktYmpU6dyyy23HDT/iSeeGHTdb33rW2zZsoUZM2awcOFC/uIv/mLI1/693/s9Pve5z3HiiSfy1re+taZ1S5KkxjHgSVKBrV+/HoDXv/71w643a9YsNm/ezPLlyxtRliRJqhMDniQVWKUBr7u7m3PPPZczzzzzgPlre5+rW22SJKn2DHiSVGDr169n4sSJzJw5c9j1HnzwQc444wwDnSRJbc7bJEhSgX3ta19j165dRMSw661du3b/bRYkSVL7sgdPksRdd90FeEqmJEntzh48SZIkSSoIA54kSZIkFYQBT5IkSZIKwoAnSZIkSQVhwJMkSZKkgjDgSZIkSVJBdMRtEo5f8MMDph9feFGTKlGnGXjsSZIkSfVkD54kSZIkFURH9OANVN6rYm+eVH+ZSUQ0uwzVWWY2uwRJkjqePXiS6qqrq4u+vr5ml6EG6Ovro6urq9llSJLU0Qx4kupqwoQJ7Nq1q9llqAF27drFhAkTml2GJEkdzYAnqa4mT57M9u3b7cUruL6+PrZv387kyZObXYokSR2tI6/Bk9Q43d3dvPjii2zatInDDjuMSZMm0dXV5TV5BZCZ9PX1sWvXLrZv387EiRPp7u5udlmSJHU0A56kuooIDj/8cHbu3MmOHTvYunWrvXkt7OlnX6x43XE7X0lXVxcTJkxg2rRpdHd3G9wlSWoyA55UY9777mARweTJkz19rw3MHcXx6yjEkiS1Hq/BkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJkiRJBeGNzqU2UH7zdG8uLbWHiJgD/C3QBdycmQsHWec84AZgHPBMZp7b0CIlSYVTVQ9eRMyJiA0RsTEiFgyyfF5ErI2I+yNiTUS8qZrtSZLUDiKiC1gEzAVmA++NiNkD1pkKfAV4W2a+Hri44YVKkgpnzAGvksYLuAs4NTNPA/4HcPNYtydJUhs5G9iYmY9l5h5gGTBvwDqXAssz8wmAzNza4BolSQVUTQ/eiI1XZu7KzCxNTgQSSZKKbwawuWy6tzSv3GuAV0XEP0fEvRFx+WAvFBHzS2fBrNm2bVudypUkFUU1Aa+SxouIeEdE/AL4If29eIOyAZMkFUgMMm/gl5yHAmcCFwF/AvxlRLzmoCdlLsnMnszsmT59eu0rlSQVSjUBr5LGi8z8bma+Fng78IWhXswGTJJUIL3AsWXTxwBbBlnnR5m5OzOfAe4GTm1QfZKkgqom4FXSeO2XmXcDJ0bEtCq2KUlSO1gNzIqIEyJiPHAJsGLAOt8D3hwRh0bEBOANwMMNrlOSVDDV3CZhf+MFPEl/43Vp+QoR8XvALzMzI+IMYDzw6yq2KUlSy8vMvRFxNXAn/bdJWJqZ6yPiqtLyxZn5cET8CFgLvEz/rRTWNa9qSVIRjDngVdJ4Af8duDwiXgJeBN5TNuiKJEmFlZkrgZUD5i0eMP1XwF81si5JUrFVdaPzkRqvzLwOuK6abUiSJEmSKlPVjc4lSZIkSa3DgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIKoaRbOIjl/ww/2PH194URMrkSRJkqTRsQdPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkF0/CAr5YOqSFIn8u+gJEnFYQ+eJEmSJBWEAU+SJEmSCsKAJ0mSJEkF0fHX4Em14DVMkiRJagX24EmSJElSQRjwJEmSJKkgDHiSJEmSVBAGPEmSJEkqCAdZkdrMwAFdHl94UZMqkSRJUquxB0+SJEmSCsKAJ0mSJEkFYcCTJEmSpIIw4EmSJElSQRjwJEmSJKkgHEVTKhBH2JQkSeps9uBJkiRJUkHYgye1uYG9dpIkSepcBjxJ6jB+KSBJUnF5iqYkSZIkFYQBT5IkSZIKwoAnSZIkSQXhNXjDcMh5SZIkSe3EHjxJkiRJKgh78KQxcBRCSZIktSJ78CRJkiSpIAx4kiRJklQQBjxJkuogIuZExIaI2BgRC4ZZ76yI6IuIdzWyPklSMVV1DV5EzAH+FugCbs7MhQOWvw+4pjS5C/izzHygmm1KkkbP60YbKyK6gEXABUAvsDoiVmTmQ4Osdx1wZ+OrlCQV0Zh78Moar7nAbOC9ETF7wGq/As7NzFOALwBLxro9SZLayNnAxsx8LDP3AMuAeYOs91HgNmBrI4uTJBVXNadojth4ZebPMvPZ0uQq4JgqtidJUruYAWwum+4tzdsvImYA7wAWD/dCETE/ItZExJpt27bVvFBJUrFUE/BGbLwG+ABwx1ALbcAkSQUSg8zLAdM3ANdkZt9wL5SZSzKzJzN7pk+fXrMCJUnFVM01eJU0Xv0rRpxPf8B701AvlplLKJ3C2dPTM+jrSJLUJnqBY8umjwG2DFinB1gWEQDTgAsjYm9m3t6YEiVJRVRNwKuk8SIiTgFuBuZm5q+r2J4kSe1iNTArIk4AngQuAS4tXyEzT9j3OCK+AfzAcCdJqlY1AW/ExisiZgLLgcsy85EqttUSykehe3zhRU2sRJLUyjJzb0RcTf/omF3A0sxcHxFXlZYPe92dJEljNeaAV2Hj9Rngd4GvlE5B2ZuZPdWXLUlSa8vMlcDKAfMGDXaZeWUjapIkFV9V98EbqfHKzA8CH6xmG5LGzl5nSZKkzlLNKJqSJEmSpBZiwJMkSZKkgjDgSZIkSVJBGPAkSZIkqSAMeJIkSZJUEFWNoilJal3lo6hKkqTOYMCTOtDAD/7eQkGSJKkYDHhj5AfkzmNviCRJklqdAa9GvKG0JEmSpGYz4Ekdwh5ISZKk4nMUTUmSJEkqCHvw6sDTNSVJkiQ1gz14kiRJklQQBjxJkiRJKggDniRJkiQVhNfgSWW8flLtzJFSJUmSPXiSJEmSVBAGPEmSJEkqCE/RbDBPAZRUS56WKUmSytmDJ0mSJEkFYQ+eNAR7RiRJktRuDHjqaIY4SZIkFYkBT5JazMAvHrxeV5IkVcqAJ0ltxF5nSZI0HANenflhTJIkSVKjOIqmJEmSJBWEPXjqOPaqHmy4feL1X5IkSe3DHjxJkiRJKgh78FR49tip3XkMS5KkStmDJ0mSJEkFYcCTJEmSpIIw4EmSJElSQXgNnqSKDbwWzBE2JUmSWosBT1JNGP4kSZKaz1M0JUmSJKkg7MFrIns8Rsf9JUmSJA3PgCdpWN6DTZIkqX0Y8FQI9u6p3XjMSpKkeqgq4EXEHOBvgS7g5sxcOGD5a4GvA2cAn87M/1nN9tR57D3qPOW/81YNPe1Qo5qvgjbyfcA1pcldwJ9l5gONrVKSVDRjDngR0QUsAi4AeoHVEbEiMx8qW2078DHg7VVVKUlSG6mwjfwVcG5mPhsRc4ElwBsaX60kqUiq6cE7G9iYmY8BRMQyYB6wv/HKzK3A1ojwK+4KDNdbZS+B1H48DbOjVdJG/qxs/VXAMQ2tUJJUSNUEvBnA5rLpXqr45jEi5gPzAWbOnFlFWZ2n3U8X80OwmqXRx56nHHeU0baRHwDuGGyB7aMkaTSqCXgxyLwc64tl5hL6T0+hp6dnzK9TVPUOce0eEqVWYYhTScVtZEScT3/Ae9Ngy20fJUmjUU3A6wWOLZs+BthSXTmSJBVCRW1kRJwC3AzMzcxfN6g2SVKBVRPwVgOzIuIE4EngEuDSmlSlQrBXUEOpxbHRKsdXrXrs7PkrnBHbyIiYCSwHLsvMRxpfoiSpiMYc8DJzb0RcDdxJ/xDQSzNzfURcVVq+OCKOBNYAk4GXI+ITwOzM3FGD2tUEXi8nSSOrpI0EPgP8LvCViADYm5k9zapZklQMVd0HLzNXAisHzFtc9vg/cVSwmvObfrWKWhyLrXQ8t0qvoIqhgjbyg8AHG12XJKnYqgp4am+t9MF6LNq9/qJr5GmY9ixLkiT1M+CpKQxnkiRJUu0Z8NQQBjpJkiSp/gx4BTPcqWpjDVntGM7asWYNz9EqJUmSRnZIswuQJEmSJNWGAU+SJEmSCsJTNHWQ0ZzC5uluqoTHiSRJUmPYgydJkiRJBWEPnlqOvT1qBd5bT5IktSMDnqRCc/RNSZLUSQx4BeeHUkmSJKlzeA2eJEmSJBWEAU+SJEmSCsKAJ0mSJEkFYcCTJEmSpIJwkBVJhePgQpIkqVPZgydJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKggDniRJkiQVhAFPkiRJkgrCgCdJkiRJBWHAkyRJkqSCMOBJkiRJUkEY8CRJkiSpIAx4kiRJklQQBjxJkiRJKoiqAl5EzImIDRGxMSIWDLI8IuLG0vK1EXFGNduTJKld2EZKkpphzAEvIrqARcBcYDbw3oiYPWC1ucCs0s984Ktj3Z4kSe3CNlKS1CzV9OCdDWzMzMcycw+wDJg3YJ15wC3ZbxUwNSKOqmKbkiS1A9tISVJTHFrFc2cAm8ume4E3VLDODOCpgS8WEfPp/wYTYFdEbKiitkpMA56p8zaaxffWnnxv7alj31tcV5NtHFeTV2k9NWsjB7SPv42IdbUttdCK/P+zHtxfo+P+Gh331+idNJYnVRPwYpB5OYZ1+mdmLgGWVFHPqETEmszsadT2Gsn31p58b+3J96Yh1KyNLG8f/Z2MjvtrdNxfo+P+Gh331+hFxJqxPK+aUzR7gWPLpo8BtoxhHUmSisY2UpLUFNUEvNXArIg4ISLGA5cAKwasswK4vDRS2DnA85l50OmZkiQVjG2kJKkpxnyKZmbujYirgTuBLmBpZq6PiKtKyxcDK4ELgY3AC8D7qy+5Zhp2OmgT+N7ak++tPfnedJA6tpH+TkbH/TU67q/RcX+Njvtr9Ma0zyJz0EviJEmSJEltpqobnUuSJEmSWocBT5IkSZIKomMDXkR8ISLWRsT9EfHjiDi62TXVSkT8VUT8ovT+vhsRU5tdUy1FxMURsT4iXo6Ith9uNyLmRMSGiNgYEQuaXU8tRcTSiNhaxPt2RcSxEfF/IuLh0vH48WbXVCsR8YqI+I+IeKD03j7f7Jo6zUh/F0oDs9xYWr42Is5oRp2tooL99b7SflobET+LiFObUWerqLTdiYizIqIvIt7VyPpaTSX7KyLOK32mXB8RP210ja2kgv+PUyLi+2VtTCuN0dFwI31WGtPf+8zsyB9gctnjjwGLm11TDd/bHwOHlh5fB1zX7Jpq/P5eR/+NH/8Z6Gl2PVW+ly7gl8CrgfHAA8DsZtdVw/f3B8AZwLpm11KH93YUcEbpcTfwSFF+d/Tfn21S6fE44N+Bc5pdV6f8VPJ3gf7BWe4o/a7OAf692XW3+P56I/Cq0uO57q+R253Sev+b/sGA3tXsult5fwFTgYeAmaXpw5tdd4vvr/9332dTYDqwHRjf7NqbuM+G/aw0lr/3HduDl5k7yiYnMsQN2NtRZv44M/eWJlfRf2+lwsjMhzNzQ7PrqJGzgY2Z+Vhm7gGWAfOaXFPNZObd9P/hLpzMfCozf156vBN4GJjR3KpqI/vtKk2OK/0U5m9kG6jk78I84JbS72oVMDUijmp0oS1ixP2VmT/LzGdLk4VrF0ep0nbno8BtwNZGFteCKtlflwLLM/P/b+9+QqyqwzCOfx+0wrBFZFRYMRJWi+gPBQUVVAqahOAuCI1qIyHUbqgBXUTgKoIkXGjUwooIS8kobCG1SIpgaAgjpEglNwa10NXo0+LcmCGm7u9Od+6555znsztz7+I979z7nvOe3597CsB2l3NWki8DV0kSsJLqPmGWjiq4Vxq43ne2wQOQ9Kqk08BTwM6641kiz1J1/TGeVgOn5x2foSVNQpdImgDuoRrpagVJyyRNU93cHbXdmnNrgJK6kNoxZ9BcPEe3r4t98yVpNbAF2DvCuMZVyefrVuBqScckfSdp28iiGz8l+dpDNRvrN2AGeMH2pdGE10gD1/tF/w5eE0j6Arh+gZembB+yPQVMSXoJ2AHsGmmA/0O/c+u9Z4rqiciBUcY2DCXn1xJa4G8ZKWkQSSupnnK/+I+ZAY1m+yJwd28N70eS7rDdurWUY6qkLqR2zCnOhaRHqRq8h5Y0ovFWkq/XgUnbF6tBlk4ryddy4F5gHbAC+FrScds/LXVwY6gkXxuAaeAx4BbgqKSv2nQNHbKB632rGzzb6wvf+i5whAY1eP3OTdLTwBPAOvcm8DbJAP+7pjsD3DTv+EaqJ1rRAJIuo2ruDtg+WHc8S8H2H5KOARuBNHijUVIXUjvmFOVC0p3APuBx27+PKLZxVJKv+4D3e83dKmCTpFnbH48mxLFS+n08Z/s8cF7Sl8BdVGuzu6YkX88Au3v3pycl/QLcDnwzmhAbZ+B639kpmpLWzjvcDPxYVyzDJmkjMAlstn2h7njiP30LrJW0RtLlwJPA4ZpjigK9tQP7gRO2X6s7nmGSdO3fu+9KWgGsp0U1sgFK6sJhYFtvd7UHgD9tnx11oGOib74k3QwcBLZ2dFRlvr75sr3G9oTtCeBD4PmONndQ9n08BDwsabmkK4H7qdZld1FJvk5RjXYi6TqqjfN+HmmUzTJwvW/1CF4fuyXdBlwCfgW21xzPMO0BrqAa8gY4brs15ydpC/AG1c5LRyRN295Qc1iLYntW0g7gc6qdp96y/UPNYQ2NpPeAR4BVks4Au2zvrzeqoXkQ2ArM9NaqAbxs+9MaYxqWG4B3JC2jehD4ge1Pao6pM/6tLkja3nt9L9XOhpuAk8AFqifinVSYr53ANcCbvevirO3G/8zOYhTmK3pK8mX7hKTPgO+p7iv3dXVKe+Hn6xXgbUkzVNMPJ22fqy3omi10r0S1udmi670aOHsvIiIiIiIiFtDZKZoRERERERFtkwYvIiIiIiKiJdLgRUREREREtEQavIiIiIiIiJZIgxcREREREdESafAiIiIiIiJaIg1eRERERERES/wFhu0CpfWJJB4AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x1440 with 6 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.subplot(4,2,1)\n",
"plt.hist(MC_tuple['q2']/1e6,density=True, label=r'$\\frac{dR_{0}}{dq^{2}}$', bins=100);\n",
"plt.legend(fontsize='25');\n",
"plt.subplot(4,2,2) \n",
"plt.hist(MC_tuple['mBbar2']/1e6,density=True,label=r'$\\frac{dR_{0}}{d\\bar{p}^{2}_{B}}$',bins=100);\n",
"plt.legend(fontsize='25');\n",
"plt.subplot(4,2,3) \n",
"plt.hist(MC_tuple['cos_theta_l'],density=True,label=r'$\\frac{dR_{0}}{dcos\\theta_{l}}$',bins=100);\n",
"plt.legend(fontsize='25');\n",
"plt.subplot(4,2,4) \n",
"plt.hist(MC_tuple['cos_theta_gamma'],density=True,label=r'$\\frac{dR_{0}}{dcos\\theta_{\\gamma}}$',bins=100);\n",
"plt.legend(fontsize='25');\n",
"plt.subplot(4,2,5) \n",
"plt.hist(MC_tuple['phi_gamma'],density=True,label=r'$\\frac{dR_{0}}{d\\phi_{\\gamma}}$',bins=100);\n",
"plt.legend(fontsize='25');\n",
"plt.subplot(4,2,6)\n",
"#plt.hist(MC_tuple['phi_l'],density=True,label=r'$\\frac{dR_{0}}{d\\phi_{l}}$',bins=100);\n",
"#plt.legend(fontsize='25');\n",
"#plt.subplot(4,2,7) \n",
"#plt.hist(MC_tuple['cos_theta_Kst'],density=True,label=r'$\\frac{dR_{0}}{dcos\\theta_{K}}$',bins=100);\n",
"#plt.legend(fontsize='25');\n",
"#plt.subplot(4,2,8)\n",
"#plt.hist(MC_tuple['phi_Kst'],density=True,label=r'$\\frac{dR_{0}}{d\\phi_{K}}$',bins=100);\n",
"#plt.legend(fontsize='25');\n",
"\n",
"\n",
"\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,20);"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"sample_for_hist=np.array([\n",
" \n",
" MC_tuple['q2'],\n",
" MC_tuple['mBbar2'],\n",
" MC_tuple['cos_theta_l'],\n",
" #MC_tuple['phi_l'],\n",
" MC_tuple['cos_theta_gamma'],\n",
" MC_tuple['phi_gamma'],\n",
" #MC_tuple['cos_theta_Kst'],\n",
" #MC_tuple['phi_Kst'],\n",
"]\n",
").reshape(len(obs_list),n_events).transpose()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Mass checks ok\n",
"B rest frame ok\n"
]
}
],
"source": [
"Kst_mass_check = (np.sqrt(scalar_product_4_np(pKst_Bbar,pKst_Bbar))).mean()-mKst_mass<10**-10\n",
"lepton_mass_check = (np.sqrt(scalar_product_4_np(p2_Bbar,p2_Bbar))).mean()-l_mass<10**-10\n",
"B_mass_check = (np.sqrt(scalar_product_4_np(pB_Bbar,pB_Bbar))).mean()-mB_mass<10**10 \n",
"photon_mass_check = (np.sqrt(scalar_product_4_np(k_Bbar,k_Bbar))).mean()-mgamma_mass<10**-7\n",
"\n",
"if Kst_mass_check and lepton_mass_check and B_mass_check and photon_mass_check:\n",
" print('Mass checks ok')\n",
"else:\n",
" print('Mass checks not good')\n",
" print('Kst mass', Kst_mass_check)\n",
" print('lepton mass', lepton_mass_check)\n",
" print('B mass', B_mass_check)\n",
" print('photon_mass_check', photon_mass_check)\n",
" \n",
"pB_Bbar_to_boost=pB_Bbar\n",
"pB = lorentz_boost_np(pB_Bbar_to_boost, -pB_Bbar/pB_Bbar[:,0:1])\n",
"\n",
"spatial_rest_check=pB[:,1].mean()<10**-10 and pB[:,2].mean()<10**-10 and pB[:,3].mean()<10**-10\n",
"Brest_mass_check=pB[:,0].mean()==mB_mass\n",
"\n",
"if spatial_rest_check and Brest_mass_check:\n",
" print('B rest frame ok')\n",
"else:\n",
" print('B rest frame not good')"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA30AAAEvCAYAAADxU6hsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3df5CdVZ3n8fd3OwYISUxCCGYCmOBGMaAiRILioK6lEtAKoLjAjmRSzIQomZGadTQ6S4kOrtFxGSvIEGItEgRhGFwlU0SzqRSs7g4gP5IBA5WiDSHEpJKQOJAQIKT57h/9JNPpdN/79K/c7ue+X1Vd997znPPc88CpW3w45zlPZCaSJEmSpGr6D43ugCRJkiRp4Bj6JEmSJKnCDH2SJEmSVGGGPkmSJEmqMEOfJEmSJFWYoU+SJEmSKmxYozvQH8aPH5+TJ09udDckSZIkqSEee+yxFzLz2K6OVSL0TZ48mUcffbTR3ZAkSZKkhoiI57o75vJOSZIkSaowQ58kSZIkVZihT5IkSZIqzNAnSZIkSRVm6JMkSZKkCjP0SZIkSVKFGfokSZIkqcIMfZIkSZJUYYY+SZIkSaowQ58kSZIkVdiwRnegkTKTXbt28dJLL7Fnzx7a2toa3SUNUS0tLYwYMYLRo0czatQoIqLRXZIkSZKAJg59mcm2bdt4+eWXGTduHG95y1toaWnxP9bVY5lJW1sbu3fv5oUXXuCVV15hwoQJjiVJkgapyQvu6/bYhoXnH8aelFerzzB4+63BoWlD365du3j55Zd561vfSktLS6O7oyEsIhg2bBhjxoxh1KhRPPfcc+zatYvRo0c3umuSJElS897T99JLLzFu3DgDn/pVS0sL48aN46WXXmp0VyRJkiSgiUPfnj17GDlyZKO7oQoaOXIke/bsaXQ3JEmSJKCJQ19bW5uzfBoQLS0tbgokSZKkQaNpQx/gRhsaEI4rSZIkDSZNHfokSZIkqeoMfZIkSZJUYYY+SZIkSaqwpn1OnyRJUlUM1oeND9Z+Sc3GmT5JkiRJqjBDnyrrX/7lX2hpaeG3v/1tr9rfcsstjB07lp07d/ZzzyRJkqTDp1Toi4hzI2JdRLRGxIIujkdELCqOPxERp9drGxHjImJlRDxTvI4tyv9LRKzp8PdGRJzWHxerweuZZ54hIrr9O/bYY3t8zi996Ut8+tOf5tRTTwXg8ccf56/+6q94z3vew+jRozn22GM555xz+PnPf95l+8svv5yxY8fyrW99q0/XJkmSJDVS3Xv6IqIFuBH4GLAJeCQilmXmUx2qzQSmFn8zgJuAGXXaLgBWZebCIgwuAL6SmXcAdxTf/S7g3sxc0z+Xq8FqzZr2f8V/+Zd/yfve975Djvc09N1///08+OCD/PrXvz5Q9t3vfpeVK1dy0UUX8fnPf55XX32VO++8kwsvvJBrrrmGb37zmwedY9iwYVx55ZVce+21XHPNNYwZM6YXVyZJkiQ1VpmNXM4EWjNzPUBE3AXMAjqGvlnAbZmZwEMRMSYiJgKTa7SdBXy4aL8UeAD4SqfvvhS4s8dXpSFn9erVAHzhC1/gHe94R5/Pd/PNNzN58mTOPvvsA2V/8Rd/wa233sqRRx55UNkHP/hBvv3tb3P11Vczbty4g85z2WWX8dWvfpXbb7+d+fPn97lfkiRJ0uFWZnnnJOD5Dp83FWVl6tRqe1xmbgEoXid08d3/GUNfU1i9ejWjR4/m7W9/e5/P9frrr/PP//zPfOITnyAiDpSfffbZBwU+gJaWFi666CL27dvHunXrDjnXCSecwLRp07jnnnv63C9JkiSpEcrM9EUXZVmyTpm2XX9pxAxgT2Z2uQtHRMwF5gKceOKJZU6pQWz16tWceuqp7Nix45Bjw4cPZ/To0aXP9dhjj7Fnz54ul4l2ZfPmzQBMmNDV/3eAGTNm8JOf/IS9e/cyfPjw0v2QJKmsWo82AB9vIKlvyoS+TcAJHT4fD2wuWWd4jbZbI2JiZm4ploJu63TOS6gxy5eZS4AlANOnTy8VJDU4bdmyha1bt7J169Yu792bPXs2t956a+nzPf300wCcdNJJdev+/ve/50c/+hEzZszgbW97W5d1TjrpJF599VXWr1/PySefXLofkiQNBgZKSWVC3yPA1IiYAvye9jB2Wac6y4D5xT17M4AXizC3vUbbZcBsYGHxeu/+k0XEfwAuBs7p7YVp6Nh/P983vvENPvCBDxxyvLsw1p3t27cDMHbs2Jr19uzZw4UXXsjevXv54Q9/2G29Y4455sB5DX2SJEkaauqGvszcFxHzgRVAC3BLZq6NiHnF8cXAcuA8oBXYA8yp1bY49ULg7oi4AthIe8jb7xxg0/4NYFRt+0Pfpz/9aU455ZR+O2/7vkJd27t3LxdddBGrV6/mnnvu4V3vele3dd944w2Ag+4PlCRJ9WcR+9K2ijOQzXjNGhzKzPSRmctpD3YdyxZ3eJ/AVWXbFuU7gI920+YB4KwyfRtoffkxG0wG84/I6tWrGTZsWL9s4gL//niHP/zhD10ef/311/nsZz/LypUrWbp0KbNmzap5vv3nGT9+fL/0T5IkSTqcSj2cXRpIq1evZurUqbzpTW/ql/O9853vBKC1tfWQY21tbVx22WXce++9LF68mD/5kz+pe77f/e53HHXUUaXuEZQkSZIGG0Of+sXf/u3fEhF85CMfYdeuXXzrW9/izDPPZOzYsRxxxBGcfPLJfPvb36atre2gdi+++CLPPvss06ZN67e+nHHGGYwYMYLf/OY3B5W/8cYbzJ49m3vuuYfvf//7/Pmf/3mp8z388MPMmDHDnTslSZI0JJVa3inV8/jjjwPtm5685z3v4dlnnyUiOProo9m7dy/r1q3ja1/7GmvXruX2228/0G7NmjVkJq+++upB5fsNHz6cz372sweVRQQf+tCHeOCBB7rsy5ve9CY+9alPsWLFCjLzwL14f/3Xf80dd9zB+9//fo455phDvu8DH/jAIbN5Gzdu5Omnn+bzn/98j/+ZSJIkSYOBoU/9Yn/o++lPf8qoUaNYtGgRc+bMYeTIkaxfv54/+7M/4/777+eOO+7gyiuv5I//+I+Bf9/E5b777uO++w69f/K9733vQaFv165dAEycOLFmf6688kr+8R//kV/96ld86EMfAtqf3wfw4IMP8uCDDx7S5kc/+tEhoe/OO+/kiCOOKLUMVJIkqS+qspeEBh9Dn/ps586dbNy4EYCRI0dy//33c8YZZxw4ftJJJ/FP//RPTJkyhV27drFs2bIDoe/qq6/m6quvLv1dDzzwABHB1772tZr1PvKRj/D+97+fG2644UDo625msDv79u3j5ptv5gtf+AJjxozpUVtJ0tDUjLsr1rrmKl6v1Iy8p099tn+WD9rv7esY+PY75phjDjyDb/363j+JY+XKlVxyySU1H7Gw39/93d/xs5/9jCeffLJX3/XjH/+YP/zhD/zN3/xNr9pLkiRJg4Ezfeqz/aFv7NixzJs3r9t69ZZklrFo0aLSdc8+++xDNo7piTlz5jBnzpxet5ckSdXTjLPBGvqc6VOf7b8v75Of/CRHHnlkt/V27NgBwHHHHXdY+iVJkiTJmT71g/0zfdOnT69Z7+GHHwbaN2eRJEkajNxMRVXkTJ/6ZPfu3Qcegj5+/Phu6z3wwANs27aNiOC888470Pbaa6/lk5/8JG95y1uICP70T//0cHRbkiRJahrO9KlP1qxZwxtvvAHA9u3bu633ne98B4DPfOYzTJo0CYAXXniBb3zjG0ycOJHp06d3+cgGSZKknnK2TjqYM33qk447d/7iF7/oss6SJUv45S9/yZFHHsk3v/nNA+UTJ05k06ZNbN68mXvuuWfA+ypJkiQ1I2f61Cf7Q9+YMWNYsWIF11xzDV/+8pcZNWoUW7Zs4frrr+f6668H4Oabb+bkk08+0PaII444MOsnSZJUBc4yajBypk99sn/nzmuvvZZTTz2V6667jje/+c28+c1v5o/+6I/43ve+R0Tw/e9/n8svv7zBvZUkSZKaj6FPvfbaa6/x1FNPAXDWWWfxq1/9ivnz53P88cfz6quvMmnSJD73uc/x+OOP88UvfrHBvZUkSZKak8s71WtPPPEE+/bto6WlhXe/+90cddRR3HDDDdxwww2N7pokSQOi1tK9gXwot0sGD+Y/D6lnDH11DOQP+FC3f2nnO97xDo466qgG90aSdLg1KgBJknrG5Z3qtf2buJx22mkN7okkSZKk7jjTp17bH/re+973NrgnkiRpIDRyGaVLOKX+40yfemXfvn08+eSTgDN9kiRJ0mDmTJ96ZdiwYbzyyit9Ps8PfvAD/u3f/o19+/YB7ZvDXHfddQCcc845nHPOOX3+DkmSJKmZGfrUUN/73vd47rnnDnxevXr1gQ1ivv71rxv6JEmS+oEbLzU3Q58aasOGDY3ugiRJklRp3tMnSZIkSRVm6JMkSZKkCnN5pyRJaho+BkBV5dhWLaVm+iLi3IhYFxGtEbGgi+MREYuK409ExOn12kbEuIhYGRHPFK9jOxx7d0Q8GBFrI+LJiDiyrxcqSZIkSc2obuiLiBbgRmAmMA24NCKmdao2E5ha/M0FbirRdgGwKjOnAquKz0TEMOB2YF5mngJ8GHi995coSZIkSc2rzEzfmUBrZq7PzL3AXcCsTnVmAbdlu4eAMRExsU7bWcDS4v1S4ILi/ceBJzLzXwEyc0dmtvXy+iRJkiSpqZUJfZOA5zt83lSUlalTq+1xmbkFoHidUJS/HciIWBERj0fEl8tciCRJkiTpUGU2cokuyrJknTJtu+rTB4H3AXuAVRHxWGauOugLI+bSvpSUE088sc4pJUlST7kxhCRVQ5nQtwk4ocPn44HNJesMr9F2a0RMzMwtxVLQbR3O9X8y8wWAiFgOnE77fX8HZOYSYAnA9OnT6wXJLmUmEV3lUqn3Mns1HCVJkhqi3v/g2bDw/MPUEw2UMss7HwGmRsSUiBgOXAIs61RnGXB5sYvnWcCLxZLNWm2XAbOL97OBe4v3K4B3R8SIYlOXDwFP9fL6utXS0kJbm7cKqv+1tbXR0tLS6G5IkiRJQImZvszcFxHzaQ9jLcAtmbk2IuYVxxcDy4HzgFbal2TOqdW2OPVC4O6IuALYCFxctPlDRFxPe2BMYHlm9vv6khEjRrB7927GjBnT36dWk9u9ezcjRoxodDckSZIkoOTD2TNzOe3BrmPZ4g7vE7iqbNuifAfw0W7a3E77YxsGzOjRo3nhhRcYNWqUszLqN21tbezcuZPx48c3uiuSJEkSUDL0VdGoUaN45ZVXeO655xg3bhwjR46kpaXFe/zUY5lJW1sbu3fvZufOnRx99NGMGjWq0d2SJEmSgCYOfRHBhAkT2LVrFy+99BLbtm3zHj/1WktLCyNGjGD8+PGMGjXK/3kgSZKkQaNpQx+0B7/Ro0czevToRndFkqQBUWtXPnfkkzTQ3Bl0cCize6ckSZIkaYhq6pk+SZLUGP7ff0k6fJzpkyRJkqQKM/RJkiRJUoW5vFOSJA06fdmApt7S0UYZrP2SVH2GPkmSNKQYniSpZ1zeKUmSJEkV5kyfJEnqd87GSdLg4UyfJEmSJFWYoU+SJEmSKszQJ0mSJEkVZuiTJEmSpAoz9EmSJElShbl7pyRJQ5i7ZEqS6nGmT5IkSZIqzNAnSZIkSRXm8k5JkqR+4FJbSYOVM32SJEmSVGGGPkmSJEmqMJd3SpLUpFyOKEnNwZk+SZIkSaowQ58kSZIkVVip0BcR50bEuohojYgFXRyPiFhUHH8iIk6v1zYixkXEyoh4pngdW5RPjohXImJN8be4Py5UkiRJkppR3dAXES3AjcBMYBpwaURM61RtJjC1+JsL3FSi7QJgVWZOBVYVn/f7XWaeVvzN6+3FSZIkSVKzKzPTdybQmpnrM3MvcBcwq1OdWcBt2e4hYExETKzTdhawtHi/FLigj9ciSZIkSeqkzO6dk4DnO3zeBMwoUWdSnbbHZeYWgMzcEhETOtSbEhGrgZeA/5aZvy7RT0mSKsldNiVJfVEm9EUXZVmyTpm2nW0BTszMHRFxBvDziDglM1866Asj5tK+lJQTTzyxziklSZIkqTmVCX2bgBM6fD4e2FyyzvAabbdGxMRilm8isA0gM18DXivePxYRvwPeDjza8QszcwmwBGD69On1gqQkSQOq1mzchoXnH8aeSJJ0sDKh7xFgakRMAX4PXAJc1qnOMmB+RNxF+/LNF4swt71G22XAbGBh8XovQEQcC+zMzLaIOIn2zWHW9+EaJUkqxeAmSaqiuqEvM/dFxHxgBdAC3JKZayNiXnF8MbAcOA9oBfYAc2q1LU69ELg7Iq4ANgIXF+XnAN+MiH1AGzAvM3f2y9VKkiRJUpMpM9NHZi6nPdh1LFvc4X0CV5VtW5TvAD7aRflPgZ+W6ZckSZIkqbZSD2eXJEmSJA1NpWb6JEnSwPGRDJKkgWTokyRVRr3w5GYskqRm5PJOSZIkSaowQ58kSZIkVZihT5IkSZIqzNAnSZIkSRXmRi6SJEmSesXdh4cGZ/okSZIkqcKc6ZMkaYD5f8IlSY3kTJ8kSZIkVZihT5IkSZIqzNAnSZIkSRXmPX2SJJXgfXmSmpW/f0OfM32SJEmSVGGGPkmSJEmqMEOfJEmSJFWY9/RJkiRJaoha9wtuWHj+YexJtTnTJ0mSJEkVZuiTJEmSpAoz9EmSJElShRn6JEmSJKnCDH2SJEmSVGHu3ilJTarWjmkweHdNq9dvSZJ0MEOfJGlQMdRJktS/XN4pSZIkSRVWKvRFxLkRsS4iWiNiQRfHIyIWFcefiIjT67WNiHERsTIinilex3Y654kRsTsivtSXC5QkSZKkZlY39EVEC3AjMBOYBlwaEdM6VZsJTC3+5gI3lWi7AFiVmVOBVcXnjv4e+EUvrkmSJEmSVCgz03cm0JqZ6zNzL3AXMKtTnVnAbdnuIWBMREys03YWsLR4vxS4YP/JIuICYD2wtpfXJUmSJEmi3EYuk4DnO3zeBMwoUWdSnbbHZeYWgMzcEhETACLiaOArwMcAl3ZKUgU1arMWN4mRJDWjMjN90UVZlqxTpm1n3wD+PjN31+xUxNyIeDQiHt2+fXudU0qSJElScyoz07cJOKHD5+OBzSXrDK/RdmtETCxm+SYC24ryGcBnIuK7wBjgjYh4NTN/0PELM3MJsARg+vTp9YKkJEmSJDWlMjN9jwBTI2JKRAwHLgGWdaqzDLi82MXzLODFYulmrbbLgNnF+9nAvQCZ+ceZOTkzJwPfB/5758AnSZIkSSqn7kxfZu6LiPnACqAFuCUz10bEvOL4YmA5cB7QCuwB5tRqW5x6IXB3RFwBbAQu7tcrkyRJkiSVWt5JZi6nPdh1LFvc4X0CV5VtW5TvAD5a53uvLdM/SZIkSVLXSoU+SZJ6wl0yJUkaPMrc0ydJkiRJGqKc6ZMk9ZgzeZIkDR3O9EmSJElShRn6JEmSJKnCDH2SJEmSVGHe0ydJkiRp0Kl3//iGhecfpp4Mfc70SZIkSVKFGfokSZIkqcIMfZIkSZJUYYY+SZIkSaowQ58kSZIkVZi7d0rSIObOZZIkqa+c6ZMkSZKkCjP0SZIkSVKFGfokSZIkqcK8p0+ShjDv+ZMkSfU40ydJkiRJFWbokyRJkqQKM/RJkiRJUoUZ+iRJkiSpwtzIRZIqrN5GL5Ikqfqc6ZMkSZKkCjP0SZIkSVKFubxTktQll4ZKklQNhj5J6iMfkC5JkgazUss7I+LciFgXEa0RsaCL4xERi4rjT0TE6fXaRsS4iFgZEc8Ur2OL8jMjYk3x968RcWF/XKgkSZIkNaO6oS8iWoAbgZnANODSiJjWqdpMYGrxNxe4qUTbBcCqzJwKrCo+A/wWmJ6ZpwHnAjdHhDOSkiRJktQLZWb6zgRaM3N9Zu4F7gJmdaozC7gt2z0EjImIiXXazgKWFu+XAhcAZOaezNxXlB8JZC+vTZIkSZKaXpnQNwl4vsPnTUVZmTq12h6XmVsAitcJ+ytFxIyIWAs8CczrEAIlSZIkST1QJvRFF2WdZ9+6q1Om7aEVMh/OzFOA9wFfjYgjD+lUxNyIeDQiHt2+fXu9U0qSJElSUyoT+jYBJ3T4fDywuWSdWm23FktAKV63df7izHwaeBk4tYtjSzJzemZOP/bYY0tchiRJkiQ1nzKh7xFgakRMiYjhwCXAsk51lgGXF7t4ngW8WCzZrNV2GTC7eD8buBegqDuseP9W4B3Aht5eoCRJkiQ1s7q7YmbmvoiYD6wAWoBbMnNtRMwrji8GlgPnAa3AHmBOrbbFqRcCd0fEFcBG4OKi/IPAgoh4HXgD+EJmvtAvVytJkiRJTabUoxAyczntwa5j2eIO7xO4qmzbonwH8NEuyn8M/LhMvyRJkiRJtZV6OLskSZIkaWgy9EmSJElShRn6JEmSJKnCDH2SJEmSVGGGPkmSJEmqMEOfJEmSJFWYoU+SJEmSKszQJ0mSJEkVZuiTJEmSpAoz9EmSJElShQ1rdAckNZfJC+6reXzDwvMb8t0D+b311PtnIkmS1BfO9EmSJElShTnTJ6kyBnLGzNk4SZI0VBn6JA0ZBi9JkrTfYL1tYzAy9EkSBkpJklRd3tMnSZIkSRXmTJ8kDTBnESVJUiM50ydJkiRJFWbokyRJkqQKM/RJkiRJUoV5T5+kHqt3j5rbJEuSJA0ehj5Jg4qbnkiSJPUvl3dKkiRJUoUZ+iRJkiSpwgx9kiRJklRhhj5JkiRJqrBSoS8izo2IdRHRGhELujgeEbGoOP5ERJxer21EjIuIlRHxTPE6tij/WEQ8FhFPFq//qT8uVJIkSZKaUd3QFxEtwI3ATGAacGlETOtUbSYwtfibC9xUou0CYFVmTgVWFZ8BXgA+lZnvAmYDP+711UmSJElSkyvzyIYzgdbMXA8QEXcBs4CnOtSZBdyWmQk8FBFjImIiMLlG21nAh4v2S4EHgK9k5uoO510LHBkRR2Tma726QkmHnY9dkCRJGjzKLO+cBDzf4fOmoqxMnVptj8vMLQDF64QuvvvTwOquAl9EzI2IRyPi0e3bt5e4DEmSJElqPmVCX3RRliXrlGnb9ZdGnAJ8B7iyq+OZuSQzp2fm9GOPPbbMKSVJkiSp6ZQJfZuAEzp8Ph7YXLJOrbZbiyWgFK/b9leKiOOBnwGXZ+bvSvRRkiRJktSFMqHvEWBqREyJiOHAJcCyTnWWAZcXu3ieBbxYLNms1XYZ7Ru1ULzeCxARY4D7gK9m5v/rw7VJkiRJUtOru5FLZu6LiPnACqAFuCUz10bEvOL4YmA5cB7QCuwB5tRqW5x6IXB3RFwBbAQuLsrnA/8RuCYirinKPp6ZB2YCJUmSJEnllNm9k8xcTnuw61i2uMP7BK4q27Yo3wF8tIvy64DryvRLkiRJklRbqYezS5IkSZKGplIzfZKqp9az9DYsPP8w9kSSJEkDyZk+SZIkSaowZ/qkIczZOkmSJNXjTJ8kSZIkVZihT5IkSZIqzOWdkg5Ra9moJEmShhZDn9RA9cKV9+VJkiSpr1zeKUmSJEkVZuiTJEmSpAoz9EmSJElShRn6JEmSJKnCDH2SJEmSVGGGPkmSJEmqMEOfJEmSJFWYoU+SJEmSKsyHs0sVVe/B75IkSWoOzvRJkiRJUoUZ+iRJkiSpwgx9kiRJklRh3tMnDWLelydJkqS+cqZPkiRJkirM0CdJkiRJFWbokyRJkqQKM/RJkiRJUoWVCn0RcW5ErIuI1ohY0MXxiIhFxfEnIuL0em0jYlxErIyIZ4rXsUX5MRFxf0Tsjogf9MdFSpIkSVKzqrt7Z0S0ADcCHwM2AY9ExLLMfKpDtZnA1OJvBnATMKNO2wXAqsxcWITBBcBXgFeBa4BTiz9pwNXbJXPDwvMH7NySJEnqXwP533ZDUZmZvjOB1sxcn5l7gbuAWZ3qzAJuy3YPAWMiYmKdtrOApcX7pcAFAJn5cmb+X9rDnyRJkiSpD8qEvknA8x0+byrKytSp1fa4zNwCULxOKN9tSZIkSVIZZUJfdFGWJeuUadsrETE3Ih6NiEe3b9/eH6eUJEmSpMqpe08f7bNzJ3T4fDywuWSd4TXabo2IiZm5pVgKuq0nHc/MJcASgOnTp/dLkJR6w3v2JEmSNJiVmel7BJgaEVMiYjhwCbCsU51lwOXFLp5nAS8WSzZrtV0GzC7ezwbu7eO1SJIkSZI6qTvTl5n7ImI+sAJoAW7JzLURMa84vhhYDpwHtAJ7gDm12hanXgjcHRFXABuBi/d/Z0RsAEYDwyPiAuDjnXYLlSRJkiSVUGZ5J5m5nPZg17FscYf3CVxVtm1RvgP4aDdtJpfplyRJkiSptlIPZ5ckSZIkDU2GPkmSJEmqsFLLO6X+Umunyw0Lzz+MPZEkSZKagzN9kiRJklRhhj5JkiRJqjBDnyRJkiRVmKFPkiRJkirMjVzUr2pt1NLIc7tJjCRJkpqVM32SJEmSVGGGPkmSJEmqMJd36hD1llEO1FLJgVwaKkmSJDUrZ/okSZIkqcKc6VNT6OssorOQkiRJGqqc6ZMkSZKkCjP0SZIkSVKFGfokSZIkqcIMfZIkSZJUYYY+SZIkSaowQ58kSZIkVZihT5IkSZIqzNAnSZIkSRVm6JMkSZKkChvW6A5o6Jm84L5Gd0GSJElSSYa+JmRokyRJkpqHoa+iDHaSJEmSwHv6JEmSJKnSSoW+iDg3ItZFRGtELOjieETEouL4ExFxer22ETEuIlZGxDPF69gOx75a1F8XEZ/o60VKkiRJUrOqG/oiogW4EZgJTAMujYhpnarNBKYWf3OBm0q0XQCsysypwKriM8XxS4BTgHOBfyjOI0mSJEnqoTL39J0JtGbmeoCIuAuYBTzVoc4s4LbMTOChiBgTEROByTXazgI+XLRfCjwAfKUovyszXwOejYjWog8P9v4yJUmSJKmcWvtjbFh4/mHsSf8os7xzEvB8h8+birIydWq1PS4ztwAUrxN68H2SJEmSpBLKzPRFF2VZsk6Ztr35PiJiLu1LSQF2R8S6OudthPHAC43uhCrBsaT+4DhSf3Acqb84ltQfej3PvogAAASnSURBVDWO4ju9/8K+tB1gb+3uQJnQtwk4ocPn44HNJesMr9F2a0RMzMwtxVLQbT34PjJzCbCkRP8bJiIezczpje6Hhj7HkvqD40j9wXGk/uJYUn9wHJVTZnnnI8DUiJgSEcNp32RlWac6y4DLi108zwJeLJZs1mq7DJhdvJ8N3Nuh/JKIOCIiptC+Ocxvenl9kiRJktTU6s70Zea+iJgPrABagFsyc21EzCuOLwaWA+cBrcAeYE6ttsWpFwJ3R8QVwEbg4qLN2oi4m/bNXvYBV2VmW39dsCRJkiQ1k2jfcFMDISLmFstQpT5xLKk/OI7UHxxH6i+OJfUHx1E5hj5JkiRJqrAy9/RJkiRJkoYoQ98AiYhzI2JdRLRGxIJG90eDW0RsiIgnI2JNRDxalI2LiJUR8UzxOrZD/a8WY2tdRHyicT1XI0XELRGxLSJ+26Gsx+MmIs4oxl9rRCyKiK4enaMK62YsXRsRvy9+l9ZExHkdjjmWdIiIOCEi7o+IpyNibUR8sSj3d0ml1RhH/ib1gaFvAEREC3AjMBOYBlwaEdMa2ysNAR/JzNM6bDu8AFiVmVOBVcVnirF0CXAKcC7wD8WYU/O5lfYx0FFvxs1NtD/3dGrx1/mcqr5b6frf+98Xv0unZeZycCyppn3Af83MdwJnAVcV48XfJfVEd+MI/E3qNUPfwDgTaM3M9Zm5F7gLmNXgPmnomQUsLd4vBS7oUH5XZr6Wmc/SvmvumQ3onxosM38F7OxU3KNxE+3PSR2dmQ9m+03et3VooybRzVjqjmNJXcrMLZn5ePF+F/A0MAl/l9QDNcZRdxxHJRj6BsYk4PkOnzdRe7BKCfzviHgsIuYWZccVz7ukeJ1QlDu+VEtPx82k4n3ncglgfkQ8USz/3L8kz7GkuiJiMvBe4GH8XVIvdRpH4G9Srxn6BkZX64XdJlW1nJ2Zp9O+JPiqiDinRl3Hl3qju3HjeFJ3bgLeBpwGbAH+R1HuWFJNETES+ClwdWa+VKtqF2WOJQFdjiN/k/rA0DcwNgEndPh8PLC5QX3REJCZm4vXbcDPaF+uubVYmkDxuq2o7vhSLT0dN5uK953L1eQyc2tmtmXmG8AP+fdl5I4ldSsi3kT7f6jfkZn/qyj2d0k90tU48jepbwx9A+MRYGpETImI4bTfXLqswX3SIBURR0fEqP3vgY8Dv6V9zMwuqs0G7i3eLwMuiYgjImIK7Tcm/+bw9lqDWI/GTbHUaldEnFXsanZ5hzZqYvv/I71wIe2/S+BYUjeKf+//E3g6M6/vcMjfJZXW3TjyN6lvhjW6A1WUmfsiYj6wAmgBbsnMtQ3ulgav44CfFbsIDwN+kpm/jIhHgLsj4gpgI3AxQGaujYi7gado3+Hqqsxsa0zX1UgRcSfwYWB8RGwCvg4spOfj5vO07954FPCL4k9NpJux9OGIOI325VAbgCvBsaSazgY+BzwZEWuKsq/h75J6prtxdKm/Sb0X7ZvZSJIkSZKqyOWdkiRJklRhhj5JkiRJqjBDnyRJkiRVmKFPkiRJkirM0CdJkiRJFWbokyRJkqQKM/RJkiRJUoUZ+iRJkiSpwv4//UAEXmLJt3QAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1_Bbar[:,0], label=r'$p_{1}^{E, (2)}$', range=(0,2600),density=True,bins=100);\n",
"plt.legend(fontsize='25');\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"136"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"(p1_Bbar[:,0]>2500).sum()"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA20AAAEvCAYAAADW/SmEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAf/UlEQVR4nO3df7RfZX0n+vfnnkgkJkgggBiggQ6KSFWukV8W1GF6dSptkCILXQqltqAVq7dTR36sFtvCEn+04xWmIzhFo1KUYqfQQUCkMDgLxOGHI7+Ga/hphGuCkUkQSEh47h/nG3oCJznJ+X5Pzs45r9daWee7n/3s/f2csPnmvM/z7GdXay0AAAB00/8x2QUAAACwcUIbAABAhwltAAAAHSa0AQAAdJjQBgAA0GFCGwAAQIfNmOwCkmTevHltwYIFk10GAADApLjtttseb63tMtq+ToS2BQsW5NZbb53sMgAAACZFVT28sX2mRwIAAHSY0AYAANBhQhsAAECHCW0AAAAdJrQBAAB0mNAGAADQYUIbAABAhwltAAAAHSa0AQAAdJjQBgAA0GEzJrsAAACYTK21rFq1KitXrsxTTz2VdevWTXZJbIOGhoYya9as7LDDDpkzZ06qamDnFtoA6IwFp125yf0PnfvOrVQJMF201rJs2bL88pe/zE477ZRXvOIVGRoaGugP3Ex9rbWsW7cuTz75ZB5//PE8/fTT2XXXXQd2HQltAABMW6tWrcovf/nL/Mqv/EqGhoYmuxy2UVWVGTNmZMcdd8ycOXPy8MMPZ9WqVdlhhx0Gcn73tAEAMG2tXLkyO+20k8DGwAwNDWWnnXbKypUrB3ZOoQ0AgGnrqaeeyuzZsye7DKaY2bNn56mnnhrY+YQ2AACmrXXr1hllY+CGhoYGuqCN0AYAwLRm0REGbdDXlNAGAADQYUIbAABAh40Z2qrqoqpaVlV3vaD9I1V1X1XdXVWfGdF+elUt6e17+0QUDQAAMF1sznPavpLk/CRfXd9QVW9LsijJ61prq6tq1177/kmOT/LaJK9M8t2qelVrzWPlAQAAxmHMkbbW2o1JVryg+UNJzm2tre71WdZrX5TkG6211a21B5MsSXLQAOsFAACYVsZ7T9urkhxeVbdU1X+rqjf12ucn+cmIfkt7bQAAAFvNTTfdlKGhodx1111jdx7FRRddlLlz52bFiheOX219mzM9cmPHzU1ySJI3Jbm0qvZJMtralm20E1TVyUlOTpK99tprnGUAwLZtwWlXbnL/Q+e+cytVAjC1/Mmf/El+53d+JwcccMDzbbfffnu+/vWv57rrrsuDDz6YmTNn5jWveU3++I//OEcfffQGx59wwgk5++yzc8455+Sv/uqvtnb5GxjvSNvSJP/Qhv0gyXNJ5vXa9xzRb48kj452gtbaha21ha21hbvssss4ywAAANjQ9ddfn5tvvjl/9Ed/tEH7Zz7zmSxevDgHHXRQPvOZz+TMM8/M6tWr8653vSt/9md/tkHfGTNm5JRTTsnf/M3f5Iknntia5b/IeEPbPyb510lSVa9Ksl2Sx5NckeT4qppZVXsn2TfJDwZRKAAAwOa44IILsmDBgrz5zW/eoP0jH/lIfvrTn+ZLX/pSPvjBD+ZjH/tYbrrpphxyyCH51Kc+9aKpkO9973uzevXqfP3rX9+a5b/ImNMjq+qSJG9NMq+qliY5K8lFSS7qPQZgTZITW2styd1VdWmSe5KsTfJhK0cCMN2NNQUSgOSJJ57I3LlzN7r/i1/8Yk455ZQxz/Pss8/mn/7pn/L+978/VRvevfXCEJckQ0NDOeaYY/L9738/9913Xw499NDn9+25557Zf//9c9lll+XUU0/dgu9msMYMba2192xk1/s20v+cJOf0UxQAADC9tNbyta99bYO2Z599NmeccUZWrlyZt771rZt1nttuuy1PPfVU3vSmN43duefRR4fv6Np1111ftO/ggw/O3/3d32XNmjXZbrvtNvucgzTehUgAAAAGZu7cuXnf+/5lXGjNmjU57rjjsnLlylx55ZV59atfvVnnuffee5Mk++yzz2b1/+lPf5ovf/nLOfjgg/Orv/qrL9q/zz775JlnnskDDzyQ/fbbb7POOWjjvacNAABgQqxZsybHHntsvvvd7+bb3/72Zo+yJcny5cuTZJNTLdd76qmn8q53vStr1qzJl770pVH77LzzzhucdzIYaQMAADpj9erVOfbYY3PDDTfkqquuyuGHHz6u8wwvubFxa9asyTHHHJM77rgjl112WX7t135t1H7PPfdckrzo/ritSWgDAIAxTJUFhbr+7MfVq1fnmGOOyfe+971cddVV+fVf//UtPsf6x4n94he/2GifZ599Nscdd1yuvfbaLF68OIsWLdpo3/XnmTdv3hbXMiimRwIAAJPumWeeydFHH53vfe97ufrqq8cV2JLkNa95TZJkyZIlo+5ft25d3vve9+byyy/PF7/4xQ3uoxvN/fffn+23336z75GbCEIbAAAwEH/5l3+Zqsrb3va2rFq1Kuecc04OOuigzJ07NzNnzsx+++2XT33qU1m3bsOngj3zzDNZtGhRbrrpplxzzTU57LDDxl3DG9/4xsyaNSs/+MGLHxf93HPP5cQTT8xll12Wz3/+8/mDP/iDMc93yy235OCDD560lSMT0yMBAIABuf3225MML97x+te/Pg8++GCqKi972cuyZs2a3HfffTnjjDNy9913b/DA6hNOOCHf+c53cuqpp+b+++/P/fff//y+7bbbLscdd9zz21WVt7zlLbnhhhtGreElL3lJfuu3fivXXHNNWmsb3Iv28Y9/PBdffHEOPfTQ7Lzzzi96aPZhhx22wYjaI488knvvvTcf+tCH+vp76ZfQBgAADMT60Patb30rc+bMyRe+8IWcdNJJmT17dh544IH8/u//fq6//vpcfPHFOeWUU3L44YentZarr746SXL++efn/PPP3+Ccb3jDG54PbatWrUqS7L777pus45RTTsk3v/nN3HjjjXnLW97yfPttt92WJLn55ptz8803v+i4L3/5yxuEtksuuSQzZ84ccwrlRBPaAACAvq1YsSKPPPJIkmT27Nm5/vrr88Y3vvH5/fvss0/+/u//PnvvvXdWrVqVK664IocffniqKitXrtys97jhhhtSVTnjjDM22e9tb3tbDj300Jx33nkbhLaNjc6NZu3atbngggvyh3/4h9lxxx03+7iJILQBTFGbWums66uHjddEre42Vf++AAZp/ShbMnxv28jAtt7OO++cww47LNdcc00eeOCBLX6Pa6+9Nscff/xGl+cf6bOf/WyOOOKI3HnnnZvV/4W+9rWv5Re/+EXOPPPMLT520IQ2AACgb+tD29y5c/PBD35wo/3Gmtq4KV/4whc2u++b3/zmFy14siVOOumknHTSSeM+fpCENgDo01R5fhNAP+64444kyVFHHZWXvvSlG+3385//PEmy2267bZW6pgKhDYCtSsABmJrWj7QtXLhwk/1uueWWJMmBBx444TVNFZ7TBgAA9OXJJ598/mHW8+bN22i/G264IcuWLUtV5Td/8zefP/aTn/xkjjrqqLziFa9IVeV3f/d3t0bZ2wyhDQAA6MsPf/jDPPfcc0mS5cuXb7Tfpz/96STJsccem/nz5ydJHn/88fz5n/95br/99jFH6aYr0yMBYIoaayqqVTGBQRm5cuRVV12Vj370oy/qc+GFF+bqq6/OS1/60vzFX/zF8+277757li5dmvnz5+eZZ57J9ttvv1Vq3pYYaQMAAPqyPrTtuOOOueaaa/Knf/qnzz8I+7HHHsvHP/7xfOhDH0qSXHDBBdlvv/2eP3bmzJnPj7oxOqENAADoy/qVIz/5yU/mgAMOyNlnn52Xv/zlefnLX55XvvKV+dznPpeqyuc///mccMIJk1zttkdoAwAAxm316tW55557kiSHHHJIbrzxxpx66qnZY4898swzz2T+/Pl5//vfn9tvv33UaZOMzT1tAADAuP3oRz/K2rVrMzQ0lNe97nXZfvvtc9555+W8886b7NKmDKENAADGYOGejVs/NfLVr361RUQmiOmRAADAuK1fhOQNb3jDJFcydY0Z2qrqoqpaVlV3jbLvT6qqVdW8EW2nV9WSqrqvqt4+6IIBAIDuWB/aDjzwwEmuZOranOmRX0lyfpKvjmysqj2T/EaSR0a07Z/k+CSvTfLKJN+tqle11tYNqmAA+uf5XQAMwtq1a3PnnXcmMdI2kcYMba21G6tqwSi7/kOSf5/k8hFti5J8o7W2OsmDVbUkyUFJbu6/VAAAoEtmzJiRp59+uu/znH/++XniiSeydu3aJMOLm5x99tlJkiOOOCJHHHFE3++xLRvXQiRV9dtJftpa+59VNXLX/CTfH7G9tNcGAH0ba4QQgG3T5z73uTz88MPPb99xxx3PL3By1llnCW1bekBVzUpyZpL/a7Tdo7S1jZzn5CQnJ8lee+21pWUAAABTxEMPPTTZJXTaeFaP/NUkeyf5n1X1UJI9ktxeVa/I8MjaniP67pHk0dFO0lq7sLW2sLW2cJdddhlHGQAAAFPfFo+0tdbuTLLr+u1ecFvYWnu8qq5I8ndV9dcZXohk3yQ/GFCtAMA2YlNTWS10A7BlNmfJ/0syvJDIq6tqaVV9YGN9W2t3J7k0yT1Jrk7yYStHAgAAjN/mrB75njH2L3jB9jlJzumvLAAAAJJxrh4JMJ14phlWrQRgMo1nIRIAAAC2EqENAACgw0yPBGCgTCUEgMES2gDYYoIZMJW01lJVk10GU0hrbaDnMz0SAIBpa2hoKOvWeUIVg7Vu3boMDQ0N7HxCGwAA09asWbPy5JNPTnYZTDFPPvlkZs2aNbDzCW0AAExbO+ywQ1asWGG0jYFZt25dVqxYkR122GFg53RPGwAA09acOXPy9NNP5+GHH85OO+2U2bNnZ2hoyD1ubJHWWtatW5cnn3wyK1asyMte9rLMmTNnYOcX2oApYao+AHtT39e2+j0BdElVZdddd82qVauycuXKLFu2zKgb4zI0NJRZs2Zl3rx5mTNnzkCDv9AGAMC0VlXZYYcdBjqdDQbJPW0AAAAdJrQBAAB0mNAGAADQYe5pA5hEYy2gAlPRVF04CGCiCG0A0GFWEAXA9EgAAIAOM9IGsI0ytZKpeg1M5OiikUtgW2SkDQAAoMOENgAAgA4T2gAAADpMaAMAAOiwMUNbVV1UVcuq6q4RbZ+tqv9VVT+qqv9SVTuO2Hd6VS2pqvuq6u0TVTgAAMB0sDmrR34lyflJvjqi7dokp7fW1lbVp5OcnuQTVbV/kuOTvDbJK5N8t6pe1VpbN9iyAYDJNFVXrgToojFH2lprNyZZ8YK277TW1vY2v59kj97rRUm+0Vpb3Vp7MMmSJAcNsF4AAIBpZRD3tP1ekqt6r+cn+cmIfUt7bQAAAIxDXw/Xrqozk6xNcvH6plG6tY0ce3KSk5Nkr7326qcMAGDATH8E6I5xh7aqOjHJUUmObK2tD2ZLk+w5otseSR4d7fjW2oVJLkyShQsXjhrsAJgcfmAHgO4Y1/TIqnpHkk8k+e3W2lMjdl2R5PiqmllVeyfZN8kP+i8TAABgehpzpK2qLkny1iTzqmppkrMyvFrkzCTXVlWSfL+19sHW2t1VdWmSezI8bfLDVo4EAAZlrFHgh85951aqBGDrGTO0tdbeM0rz326i/zlJzumnKAAAAIb1tRAJwFTg/i0AoMsGseQ/AAAAE0RoAwAA6DDTIwEmkKmXAEC/hDbgRTYVNKzMBlOHXyowWawCClvG9EgAAIAOE9oAAAA6zPRIAIAJ1s9UVFMFASNtAAAAHSa0AQAAdJjQBgAA0GFCGwAAQIdZiASYFjyPChhLP88O8xkDTCShDaBPflgDACaS6ZEAAAAdZqQNAIBtyqZmOHiuHVORkTYAAIAOE9oAAAA6zPRIAIDNYNEhYLIYaQMAAOgwoQ0AAKDDTI+EbdR0XDnL1CQAYDoac6Stqi6qqmVVddeItp2q6tqq+nHv69wR+06vqiVVdV9VvX2iCgcAAJgONmek7StJzk/y1RFtpyW5rrV2blWd1tv+RFXtn+T4JK9N8sok362qV7XW1g22bGBbNNZI2VQdIQS2HiPy9Gs6zmSh+8YcaWut3ZhkxQuaFyVZ3Hu9OMnRI9q/0Vpb3Vp7MMmSJAcNqFYAAIBpZ7wLkezWWnssSXpfd+21z0/ykxH9lvbaAAAAGIdBrx5Zo7S1UTtWnVxVt1bVrcuXLx9wGQAAAFPDeFeP/FlV7d5ae6yqdk+yrNe+NMmeI/rtkeTR0U7QWrswyYVJsnDhwlGDHTAxJvKeD/eTAAAM1nhH2q5IcmLv9YlJLh/RfnxVzayqvZPsm+QH/ZUIAAAwfY050lZVlyR5a5J5VbU0yVlJzk1yaVV9IMkjSd6dJK21u6vq0iT3JFmb5MNWjgQAGD8r7wJjhrbW2ns2suvIjfQ/J8k5/RQFAADAsEEvRAIAAMAAjXchEgAApjgPmoZuENqAzrDyJADAi5keCQAA0GFG2ugMUzAAAODFjLQBAAB0mJE2mERGFwEAGIuRNgAAgA4T2gAAADrM9EgAgClqrEepmIoP2wahDQBgG9bPMy49H3N6cA/9ts/0SAAAgA4z0gYTaCr+BnMqfk8AMAhGtJgoQhsDY958dwhWAEw0/+7D1mN6JAAAQIcZaYM+GNECgNH5NxIGx0gbAABAhwltAAAAHWZ6JHSUaSUAACRG2gAAADpNaAMAAOgw0yMBAGAb5paKqa+vkbaq+r+r6u6ququqLqmql1bVTlV1bVX9uPd17qCKBQAAmG7GHdqqan6SP0qysLV2QJKhJMcnOS3Jda21fZNc19sGAABgHPqdHjkjyfZV9WySWUkeTXJ6krf29i9OckOST/T5PkwBhu4BAGDLjXukrbX20ySfS/JIkseS/O/W2neS7NZae6zX57Ekuw6iUAAAgOmon+mRc5MsSrJ3klcmeVlVvW8Ljj+5qm6tqluXL18+3jIAAACmtH6mR/6bJA+21pYnSVX9Q5LDkvysqnZvrT1WVbsnWTbawa21C5NcmCQLFy5sfdQBE8q0TgCAwdvUz1gPnfvOrVhJ9/WzeuQjSQ6pqllVVUmOTHJvkiuSnNjrc2KSy/srEQAAYPoa90hba+2Wqrosye1J1ia5I8MjZ7OTXFpVH8hwsHv3IAplehtrtMtvYwCAbZmfddiUvlaPbK2dleSsFzSvzvCoG5PE//QAADB19PVwbQAAACZWv89pAwCAzrCA2JYxQ2vbILRNQ/2s1OODEACYrvwcxGQxPRIAAKDDjLSxAb9BAgCAbhHaOsrDBgfL3ycAsC3zi/XpzfRIAACADhPaAAAAOsz0yG2Q4XEAAJg+jLQBAAB0mNAGAADQYaZHMiWYMgoAsPVZoXvrMNIGAADQYUIbAABAhwltAAAAHeaeNgAAmGDb6v3322rdU42RNgAAgA4T2gAAADpMaAMAAOgwoQ0AAKDDhDYAAIAOs3okAACdYsVC2FBfoa2qdkzyn5MckKQl+b0k9yX5ZpIFSR5Kclxr7Rd9VQkTyD8MAAB0Wb/TI/+fJFe31vZL8vok9yY5Lcl1rbV9k1zX2wYAAGAcxh3aqmqHJEck+dskaa2taa09kWRRksW9bouTHN1vkQAAANNVP9Mj90myPMmXq+r1SW5L8tEku7XWHkuS1tpjVbXraAdX1clJTk6Svfbaq48ytk2m5AEAMJX5eXdw+pkeOSPJ/5nkP7XWDkzyy2zBVMjW2oWttYWttYW77LJLH2UAAABMXf2EtqVJlrbWbultX5bhEPezqto9SXpfl/VXIgAAwPQ17umRrbX/r6p+UlWvbq3dl+TIJPf0/pyY5Nze18sHUikAADAtjDW18qFz37mVKumGfp/T9pEkF1fVdkkeSHJShkfvLq2qDyR5JMm7+3wPAACAaauv0NZa+2GShaPsOrKf8wIAADCs3+e0AQAAMIGENgAAgA4T2gAAADpMaAMAAOgwoQ0AAKDDhDYAAIAOE9oAAAA6TGgDAADoMKENAACgw4Q2AACADhPaAAAAOkxoAwAA6DChDQAAoMNmTHYBU9mC066c7BIAAIBtnJE2AACADhPaAAAAOkxoAwAA6DChDQAAoMOENgAAgA4T2gAAADpMaAMAAOiwvkNbVQ1V1R1V9V972ztV1bVV9ePe17n9lwkAADA9DWKk7aNJ7h2xfVqS61pr+ya5rrcNAADAOPQV2qpqjyTvTPKfRzQvSrK493pxkqP7eQ8AAIDprN+Rts8n+fdJnhvRtltr7bEk6X3dtc/3AAAAmLbGHdqq6qgky1prt43z+JOr6taqunX58uXjLQMAAGBK62ek7c1JfruqHkryjST/uqq+nuRnVbV7kvS+Lhvt4Nbaha21ha21hbvssksfZQAAAExd4w5trbXTW2t7tNYWJDk+yT+31t6X5IokJ/a6nZjk8r6rBAAAmKYm4jlt5yb5jar6cZLf6G0DAAAwDjMGcZLW2g1Jbui9/nmSIwdxXgAAgOluIKFtulpw2pWTXQIAADDFTcT0SAAAAAZEaAMAAOgwoQ0AAKDDhDYAAIAOE9oAAAA6TGgDAADoMKENAACgw4Q2AACADhPaAAAAOkxoAwAA6DChDQAAoMOENgAAgA4T2gAAADpsxmQXAAAAsCUWnHblJvc/dO47t1IlW4eRNgAAgA4T2gAAADpMaAMAAOgwoQ0AAKDDhDYAAIAOE9oAAAA6TGgDAADosHGHtqras6qur6p7q+ruqvpor32nqrq2qn7c+zp3cOUCAABML/2MtK1N8u9aa69JckiSD1fV/klOS3Jda23fJNf1tgEAABiHcYe21tpjrbXbe69XJbk3yfwki5Is7nVbnOTofosEAACYrgZyT1tVLUhyYJJbkuzWWnssGQ52SXYdxHsAAABMR32HtqqaneRbST7WWlu5BcedXFW3VtWty5cv77cMAACAKamv0FZVL8lwYLu4tfYPveafVdXuvf27J1k22rGttQtbawtbawt32WWXfsoAAACYsvpZPbKS/G2Se1trfz1i1xVJTuy9PjHJ5eMvDwAAYHqb0cexb07y/iR3VtUPe21nJDk3yaVV9YEkjyR5d38lAgAATF/jDm2ttf+epDay+8jxnhcAAIB/MZDVIwEAAJgYQhsAAECHCW0AAAAdJrQBAAB0mNAGAADQYUIbAABAhwltAAAAHSa0AQAAdJjQBgAA0GFCGwAAQIcJbQAAAB0mtAEAAHSY0AYAANBhQhsAAECHCW0AAAAdJrQBAAB0mNAGAADQYUIbAABAhwltAAAAHSa0AQAAdJjQBgAA0GFCGwAAQIdNWGirqndU1X1VtaSqTpuo9wEAAJjKJiS0VdVQkv+Y5N8m2T/Je6pq/4l4LwAAgKlsxgSd96AkS1prDyRJVX0jyaIk90zQ+02IBaddOdklAAAA09xETY+cn+QnI7aX9toAAADYAhM10lajtLUNOlSdnOTk3uaTVXXfBNVCd81L8vhkF8E2z3XEoLiWGATXEYPgOupTfXqyKxiXX9nYjokKbUuT7Dlie48kj47s0Fq7MMmFE/T+bAOq6tbW2sLJroNtm+uIQXEtMQiuIwbBdcQLTdT0yP+RZN+q2ruqtktyfJIrJui9AAAApqwJGWlrra2tqlOTXJNkKMlFrbW7J+K9AAAAprKJmh6Z1tq3k3x7os7PlGB6LIPgOmJQXEsMguuIQXAdsYFqrY3dCwAAgEkxUfe0AQAAMABCGxOiqj5bVf+rqn5UVf+lqnYcse/0qlpSVfdV1dtHtL+xqu7s7ftCVVWvfWZVfbPXfktVLdj63xGTpareXVV3V9VzVbXwBftcS/Stqt7Ru4aWVNVpk10P3VNVF1XVsqq6a0TbTlV1bVX9uPd17oh9W/TZxNRXVXtW1fVVdW/v37SP9tpdR2wWoY2Jcm2SA1prr0vy/yY5PUmqav8Mryb62iTvSPI3VTXUO+Y/ZfjZffv2/ryj1/6BJL9orf2rJP8hybb55A3G664kxyS5cWSja4lB6F0z/zHJv02yf5L39K4tGOkr+ZfPkfVOS3Jda23fJNf1tsf72cTUtzbJv2utvSbJIUk+3LtWXEdsFqGNCdFa+05rbW1v8/sZflZfkixK8o3W2urW2oNJliQ5qKp2T7JDa+3mNnyj5VeTHD3imMW915clOdJvlaaP1tq9rbX7RtnlWmIQDkqypLX2QGttTZJvZPg6gee11m5MsuIFzSM/TxZnw8+ZLf1sYoprrT3WWru993pVknuTzI/riM0ktLE1/F6Sq3qv5yf5yYh9S3tt83uvX9i+wTG9IPi/k+w8gfWybXAtMQgbu45gLLu11h5Lhn8gT7Jrr308n01MI72p+QcmuSWuIzbThC35z9RXVd9N8opRdp3ZWru81+fMDE8JuHj9YaP0b5to39QxTBGbcy2Ndtgoba4ltpRrgkEbz2cT00RVzU7yrSQfa62t3MRkD9cRGxDaGLfW2r/Z1P6qOjHJUUmObP/ybImlSfYc0W2PJI/22vcYpX3kMUurakaSl+fF01TYho11LW2Ea4lB2Nh1BGP5WVXt3lp7rDdlbVmvfTyfTUwDVfWSDAe2i1tr/9Brdh2xWUyPZEJU1TuSfCLJb7fWnhqx64okx/dW8ds7wzfQ/qA3JWBVVR3Su8fohCSXjzjmxN7rY5P884gQyPTlWmIQ/keSfatq76raLsM3/l8xyTWxbRj5eXJiNvyc2dLPJqa43n/zv01yb2vtr0fsch2xWTxcmwlRVUuSzEzy817T91trH+ztOzPD97mtzfD0gKt67QszvELX9hm+B+4jrbVWVS9N8rUMz/9ekeT41toDW/HbYRJV1buSnJdklyRPJPlha+3tvX2uJfpWVb+Z5PNJhpJc1Fo7Z5JLomOq6pIkb00yL8nPkpyV5B+TXJpkrySPJHl3a21Fr/8WfTZtze+FyVFVv57ke0nuTPJcr/mMDN/X5jpiTEIbAABAh5keCQAA0GFCGwAAQIcJbQAAAB0mtAEAAHSY0AYAANBhQhsAAECHCW0AAAAdJrQBAAB02P8PwlnjEvmYaRcAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1_Bbar[:,3], label=r'$p_{1}^{z, (2)}$',range=(-2600,2600),bins=100);\n",
"plt.legend(fontsize='25');\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {},
"outputs": [],
"source": [
"#plt.hist(pKst_Bbar[:,1],density=True, label=r'$p_{K}^{y, (2)}$',bins=500);\n",
"#plt.legend(fontsize='25');\n",
"#fig = plt.gcf()\n",
"#fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAA30AAAEvCAYAAADxU6hsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nO3de3yV1Z3v8e+PSECRq4lcExENIhJkVxpErUQdIYaXxPFSFa0XVLDiYPF0BpwXNc6rKji1cgQEaw1q9Ujq5SA4wEFQU7HqCCEolxkQClKMKBJAEiEGss4f2ckkZN9y2zt59uf9euWV7Gc969m/jU9Tvqz1rGXOOQEAAAAAvKldrAsAAAAAALQcQh8AAAAAeBihDwAAAAA8jNAHAAAAAB5G6AMAAAAADyP0AQAAAICHnRTrAppDUlKS69+/f6zLAAAAAICYKCws/M45lxyozROhr3///lq3bl2sywAAAACAmDCzL4O1Mb0TAAAAADyM0AcAAAAAHkboAwAAAAAPI/QBAAAAgIcR+gAAAADAwwh9AAAAAOBhhD4AAAAA8DBCHwAAAAB4mCc2Z2+I8vJylZSU6PDhwzp+/HisywFiLiEhQZ07d1aPHj3UoUOHWJcDAACAZhZXoa+8vFy7d+9W9+7d1b9/f7Vv315mFuuygJhxzqmiokLff/+9du/erdTUVIIfAACAx8TV9M6SkhJ1795dSUlJSkxMJPAh7pmZEhMTlZSUpO7du6ukpCTWJQEAAKCZxdVI3+HDh9W/f/9YlwG0Sl26dNGuXbvUu3fvWJcCwAMK5qVLZfuCtmcmJktTN0axIgCIX3EV+o4fP6727dvHugygVWrfvj3PuQJoPmX7pPH5wdvzcqJXCwDEubia3imJKZ1AEPxvAwAAwJviLvQBAAAAQDwh9AEAAACAh8XVM30AALRFoRZFYUEUAEA4hD4AAFq7UIuisCAKACAMpnciag4cOKCePXtqx44dDep3/fXX66mnnmqhqgAAAABvI/QhYpdffrnMrN7X2LFjI+r/+OOPKzs7W2eddVbNsZkzZ+qnP/2punTpouTkZF199dXatGlTnX65ubl69NFHdejQoWb9PAAAAEA8YHonIlZUVKTHHntMEyZMqHP8lFNOCdv3hx9+0PPPP6+33367zvGCggLdd999+ulPfyrnnB5++GH9wz/8g7Zs2aIePXpIktLT0zVgwAC98sormjx5cvN9IACIAzwPCAAg9CEiO3bs0MGDBzVq1Cj16tWrwf2XL1+udu3a6eKLL65zfOXKlXVev/zyy+ratav++te/6uqrr645Pm7cOC1atIjQBwANxfOAABD3mN6JiBQWFiohIUE+n69R/desWaMLLrgg7Abghw8fVmVlpbp3717neEZGhj799FMdOXKkUe8PAAAAxCtCHyJSWFio48eP6/TTT9epp55a83XDDTdE1P/LL79U7969w573wAMPaNiwYRo5cmSd43369FFFRYWKi4sbVT8AAAAQr5jeWcvsVdtiXUJIU68cGLP3Liws1HXXXadZs2bVOd61a9eI+h85ckQ9e/YMec6DDz6oDz/8UB9++KESEhLqtJ188sk11wEAAAAQOUIfIlJUVKQZM2bo7LPPblT/pKQkHThwIGj71KlTlZ+fr/fff18DBgyo115SUiJJSk5ObtT7AwAAAPEqoumdZpZlZlvNbLuZTQ/QbmY2x9/+uZn9JFxfM+thZqvM7Av/9+7+4+3N7CUz22hm/2VmDzXHB0VdH330kcxMU6ZM0aJFi/Szn/1MXbt2VceOHZWRkaH333+/5tydO3eqpKSk0c/zSZLP59OWLVsCtj3wwAN69dVX9d5772nQoEEBz9m0aZP69OkTdrQQAAAAQF1hR/rMLEHSM5KulLRH0lozW+qcq/03+Kskpfm/RkhaIGlEmL7TJb3rnJvlD4PTJU2TdIOkDs65dDM7RdIWM1vknNvVPB8ZkrR+/XpJ0qpVqzR//nyNHTtWkyZNUlFRkVavXq3s7Gxt3bpVqampKiwslCT16tVLe/furXOdpKQknXRS+AHjMWPGaNq0adq/f79OO+20muOTJ0/Wyy+/rLfeekvdu3evuX71M4PV1qxZo6ysrCZ/bgBoKWyNED2h/qwl/rwB4ESRTO/MkLTdOfc3STKzfEk5kmqHvhxJf3LOOUmfmFk3M+stqX+IvjmSMv39X5JUoKrQ5yR1MrOTJJ0s6UdJ3zf+IyKQ6tBXXFysgoICXXLJJTVtU6ZM0dy5czVnzhw9+eSTNaHv3HPPrXMNM1NJSYm6desmSXrxxRd15513aufOnerfv3+dc9PT05WRkaH8/Pw62y7Mnz9fknTFFVfUOT83N1ePPPKIJOno0aNavHhxve0dAKBVYWuE6An1Zy3x5w0AJ4hkemdfSX+v9XqP/1gk54Tq29M597Uk+b+f7j/+hqQySV9L2i3pSedcyYlFmdlEM1tnZuv27Qv+r30IrDr0zZs3r07gk6S7775bkmqmY86cOVPOuXpflZWVNYFPqpoGOnjwYPXr1y/ge+bm5mrOnDk6fvx4zbFA13XO1QQ+ScrLy9OIESN04YUXNstnBwAAAOJJJKEv0MZqLsJzIul7ogxJxyX1kXSmpP9lZvVW9nDOPeecG+6cG87iHg1TXl6uLVu2KCUlRbfccku99urplxUVFQ267vLlyzVv3ryg0z2zsrI0efJk7dmzp0HXbd++vebOndugPgAAAACqRDK9c4+klFqv+0k6cbO0YOckhuj7jZn1ds597Z8K+q3/+HhJ/885VyHpWzP7q6Thkv4WQa2IwMaNG1VRUaExY8aoXbv6uX/Xrl2SpNTU1AZdd+3atWHPmTJlSoOuKUkTJ05scB8AAAAAVSIJfWslpZnZmZK+knSTqoJZbUsl3e9/Zm+EpEP+MLcvRN+lkm6XNMv/fYn/+G5Jl5vZK5JOkXShpP/dyM+HAKqndp743F21ZcuWSZJGjx4drZIAAHGmIClVeqJX0HYWYwGA5hM29DnnjpnZ/ZJWSkqQtNA5t9nM7vW3PytpuaRsSdsl/SDpzlB9/ZeeJek1M7tLVUHvBv/xZyS9IGmTqqaHvuCc+7w5PiyqVIe+gwcP1msrKSnRc889p379+iknp+pB+A8++KBmQZfi4mK98MILuuOOO6JZMgDAa3LCTNtnMRYAaDYRbc7unFuuqmBX+9iztX52kiaf2C9YX//x/ZKuCHC8VP8TANECqlfjfOONN/TII4+oU6dOkqTS0lKNHz9e+/fv15tvvqmOHTvWHB8yZIhuu+023XbbbTGrGwAAAEDDRRT64B3Hjh3Tpk2bdP755+vQoUMaOnSoxo0bp/Lyci1ZskTFxcWaOXOmrr322po+2dnZys7OliRG+AAANULul9eJRdYAoLUg9MWZzZs36+jRoxo5cqQefPBBTZkyRXl5eXLOacSIEXrhhRd4lg8AEJlw++UBAFoFQl+cqX6ez+fzKS0tTStWrIhxRQAAAABaEqGvlqlXDox1CS2udugDAAAA4H2RbM4OD1m/fr0SEhKUnp4e61IAAAAARAEjfXGksrJSn332mQYNGlSzMicAIDIhFy2RWLgEANBqEfriSLt27VRaWtrgfqWlpdq+fbukquC4e/dubdiwQT169FBqampzlwkArROLlgAA2iimdyKsdevWyefzyefz6ciRI8rNzZXP59PDDz8c69IAAAAAhMFIH8LKzMyUcy7WZQAAAABoBEb6AAAAAMDDCH0AAAAA4GFM7wQAAEBcmb1qW9C21rpvc6iapdZbN1oHRvoAAAAAwMMY6QMAoA0rSEqVnugV/AT2DwSAuEfoAwCgLcuZG+sKAACtHNM7AQAAAMDDGOkDAABxo2BeulS2L2h7ZmKyNHVjFCsCgJZH6AMAAPGjbJ80Pj94e15O9GoBgCgh9AEA4BdyFIgFURAGo4gAWitCHwCgVYnpX5zDjQK1kHCfmcDZeoT9h4FGjiI29b5vrfvOtda6gHhD6AMAtC7xOP0uVmEzzHYPjEwF0FL/reLxvgcQNYQ+AADiVbjtHggaAOAJbNkAAAAAAB5G6EPUHDhwQD179tSOHTsa1O/666/XU0891UJVAQAAAN7G9E5E7PLLL9f7779f73h2draWLVsWtv/jjz+u7OxsnXXWWTXHnnnmGf3hD3/Qrl27JEnnnXeeZsyYobFjx9ack5ubq1GjRumuu+5S165dm/5BAMQtFkwBAMQjQh8iVlRUpMcee0wTJkyoc/yUU04J2/eHH37Q888/r7fffrvO8X79+umJJ55QWlqaKisr9dJLL+maa65RYWGhhg4dKklKT0/XgAED9Morr2jy5MnN94EAxJ8YLZgCRGLJhq+CtmVGrwwAHkToQ0R27NihgwcPatSoUerVK/hKb8EsX75c7dq108UXX1zneE5O3UUCHnvsMS1YsEAff/xxTeiTpHHjxmnRokWEPgAhMZLXvMKt7tla/zxD1t1Ka25JobZNkNg6AYgHhD5EpLCwUAkJCfL5fI3qv2bNGl1wwQUys6DnHD9+XK+//rpKS0t10UUX1WnLyMjQo48+qiNHjujkk09uVA0A4gAjec0r3OqeLahJwS2GdcebcIGyKX29GEbj8TOjdSD0ISKFhYU6fvy4Tj/99DrHr7rqKr3++uth+3/55Zfq3bt3wLaNGzdq5MiROnr0qE499VQtXrxY6enpdc7p06ePKioqVFxcXOeZQACAR8UouLXV0U0ACIXQV9vsdOnQ7lhXEVzX1JhtkltYWKjrrrtOs2bNqltShAurHDlyRD179gzYds4552jDhg06ePCg3nzzTd1+++0qKCjQkCFDas6pHt07cuRIIz8BAAARYJQQgAcR+mo7tFt65FCsqwjukditXFlUVKQZM2bo7LPPblT/pKQkHThwIGBbYmJizXWHDx+utWvXavbs2crLy6s5p6SkRJKUnMy/sAJeEPLZO0ZSAABoVoS+OPXRRx/p4osv1j/90z9p5MiRmj9/vj7//HOVl5dr6NCheuKJJ3TZZZdJknbu3KmSkpJGP88nST6fTy+++GJE51ZWVqq8vLzOsU2bNqlPnz5BRwsBtDE8ewcAQNQQ+uLU+vXrJUmrVq3S/PnzNXbsWE2aNElFRUVavXq1srOztXXrVqWmpqqwsFCS1KtXL+3du7fOdZKSknTSSeFvozFjxmjatGnav3+/TjvttJrj06dP19ixY5WSkqLDhw/r1VdfVUFBQb19/9asWaOsrKymfmwAAIBWqykL4wChEPriVHXoKy4uVkFBgS655JKatilTpmju3LmaM2eOnnzyyZrQd+6559a5hpmppKRE3bp1kyS9+OKLuvPOO7Vz507179+/zrnp6enKyMhQfn5+nW0X9u7dq1tvvVV79+5V165dNXToUK1YsUJjxoypOefo0aNavHixVq5c2ax/BgAAnCjUXnmSlDOsb5QqiZ5QQYPVJAFvIPTFqerQN2/evDqBT5LuvvtuzZ07V1u2bJEkzZw5UzNnzgx7zZ07d2rw4MHq169fwPbc3Fw98MADuvfee5WQkCBJEU35zMvL04gRI3ThhReGPRcAAC8Kt6rohMpuWjh8aRQril9su4C2iNAXh8rLy7VlyxalpKTolltuqddePf2yoqKiQdddvny55s2bF3S6Z1ZWliZPnqw9e/bojDPOiPi67du319y5rKYGAIhjYVYV7ZqXE6VCALRFhL44tHHjRlVUVGjMmDFq165dvfZdu3ZJklJTUxt03bVr14Y9Z8qUKQ26piRNnDixwX0AoDVh7zeg7eC5OngRoS8OVU/tPPG5u2rVi6iMHj06WiUBgLex9xsQVQQ3oK76wzzwvOrQd/DgwXptJSUleu6559SvXz/l5FRNFfnggw80btw49e3bV2YW8dYLAAAAAGKP0BeHqlfjfOONN1RWVlZzvLS0VOPHj9f+/fv19NNPq2PHjjXHhwwZoqefflonn3xyTGoGAAAA0DhM74wzx44d06ZNm3T++efr0KFDGjp0qMaNG6fy8nItWbJExcXFmjlzpq699tqaPtnZ2crOzpYk3XHHHTGqHAAAoPVjailaI0b64szmzZt19OhRjRw5Uu+8844GDhyovLw8vfzyyzr33HO1cuVKTZ8+PdZlAgAAAGgmjPTFmern+Xw+n9LS0rRixYoYVwSgLSqYly6V7QvanpmYLE3dGMWKgOgItXl7S27cHm7T+MwWe+fWidE0oGEIfbV1TZUe6RrrKoLr2rAtFAKpHfoAoNHK9knj84O3s2dYXIhVAGqrQm3dUda+e5SrARBPCH21xcG/Sq9fv14JCQlKT0+PdSkAAMSXEFt3rA4zkhcrsRxRYzQPaD480xdHKisr9dlnn2nQoEE1K3MCAAAA8DZG+uJIu3btVFpa2uB+paWl2r59u6Sq4Lh7925t2LBBPXr0UGpq06ecAgAAAGg5hD6EtW7dOl122WU1r3Nzc5Wbm6vbb7+djdoBAADagFDTZadeOTCKlSAWCH0IKzMzU865WJcBAAAAoBF4pg8AAAAAPIyRPgBAswu1NL0kqVNy9IoBagm33x3QVrHaKUKJKPSZWZakpyUlSHreOTfrhHbzt2dL+kHSHc659aH6mlkPSX+W1F/SLkk/d84d8LcNlfQHSV0kVUr6qXPuaFM+KAAgikIsTd9UIQMlYRIAgHrChj4zS5D0jKQrJe2RtNbMljrnttQ67SpJaf6vEZIWSBoRpu90Se8652aZ2XT/62lmdpKkVyT9wjn3mZmdJqmimT4vAKCta8FACQCAF0Uy0pchabtz7m+SZGb5knIk1Q59OZL+5KpW+/jEzLqZWW9VjeIF65sjKdPf/yVJBZKmSRot6XPn3GeS5Jzb34TPBwAxVzAvXSrbF7Q9MzFZmroxihVFJmTdjKgBANBmRBL6+kr6e63Xe1Q1mhfunL5h+vZ0zn0tSc65r83sdP/xgZKcma2UlCwp3zn37xHUCQCtU9k+aXx+8Pa8nOjV0hDh6gYAAG1CJKHPAhw7cf3+YOdE0jdQTZdI+qmqng9818wKnXPv1nlDs4mSJkpig3AAAFoAi54gnIEbc9Sp4kDQdl9lNy0cvjSKFQEIJJLQt0dSSq3X/SQVR3hOYoi+35hZb/8oX29J39a61l+cc99Jkpktl/QTSXVCn3PuOUnPSdLw4cPZRA4AACDKOlUc0JLBvw/anll4XxSrQWOFW/mTzdvbvkhC31pJaWZ2pqSvJN0kafwJ5yyVdL//mb0Rkg75w9y+EH2XSrpd0iz/9yX+4ysl/YuZnSLpR0mjJM1u5OcDAADwvGVd+si3/tKAbYy2AQgb+pxzx8zsflWFsQRJC51zm83sXn/7s5KWq2q7hu2qmpJ5Z6i+/kvPkvSamd0labekG/x9DpjZU6oKm07Scufcsub6wAAAAF6zOm1G0DZG2wBEtE+fc265qoJd7WPP1vrZSZocaV//8f2SrgjS5xVVbdsAAAAAD+J5QCB6Igp9AAAAQHPieUAgegh9AAB4WKgVOHOG9Y1iJQDiEYvEtA7tYl0AAAAAAKDlMNKHqDlw4IAGDRqkjz76SGeddVbE/a6//npddNFFevDBB1uwOgBANIXbA5BRSABoPoz0IWKXX365zKze19ixYyPq//jjjys7O7tO4HvmmWc0dOhQdenSRV26dNHIkSO1bFndxVpzc3P16KOP6tChQ836eQAAiAfV2zkE+5qwblysSwTQwhjpQ8SKior02GOPacKECXWOn3LKKWH7/vDDD3r++ef19ttv1zner18/PfHEE0pLS1NlZaVeeuklXXPNNSosLNTQoUMlSenp6RowYIBeeeUVTZ4ccJFYAIDHNOVZxHCjiLESq7pCbecgSWVfPMoef4DHEfoQkR07dujgwYMaNWqUevXq1eD+y5cvV7t27XTxxRfXOZ6Tk1Pn9WOPPaYFCxbo448/rgl9kjRu3DgtWrSI0AcAaLWhrq1ijz/A+5jeiYgUFhYqISFBPp+vUf3XrFmjCy64QGYW9Jzjx48rPz9fpaWluuiii+q0ZWRk6NNPP9WRI0ca9f7NZdeuXTIz3XHHHTGto1prqwcAAACtDyN9iEhhYaGOHz+u008/vc7xq666Sq+//nrY/l9++aV69+4dsG3jxo0aOXKkjh49qlNPPVWLFy9Wenp6nXP69OmjiooKFRcXN2gRGABAbDAaBwCtB6GvloK/F8S6hJAyUzJj9t6FhYW67rrrNGvWrDrHu3btGlH/I0eOqGfPngHbzjnnHG3YsEEHDx7Um2++qdtvv10FBQUaMmRIzTknn3xyzXUAAAAARI7Qh4gUFRVpxowZOvvssxvVPykpSQcOHAjYlpiYWHPd4cOHa+3atZo9e7by8vJqzikpKZEkJScnN+r9gbaqYF66VLYvaHtmYrI0dWMUKwIAAG0Nz/TFqY8++khmpilTpmjRokX62c9+pq5du6pjx47KyMjQ+++/X3Puzp07VVJS0ujn+STJ5/Npy5YtEZ1bWVmp8vLyOsc2bdqkPn36BB0trO2aa66RmWnu3Ln12n7zm9/IzHT33XfXa/v000914403qm/fvurQoYN69+6t0aNH67XXXgv4Prt27dJNN92kpKQkdezYUcOHD9d//Md/BDz3P//zP3X99derV69eSkxMVEpKiiZNmqTi4uKgn6Oh9dRWWVmpKVOmyMx07bXX6ujRo2H7oJUq2yeNzw/+dWh3rCsEAACtHCN9cWr9+vWSpFWrVmn+/PkaO3asJk2apKKiIq1evVrZ2dnaunWrUlNTVVhYKEnq1auX9u7dW+c6SUlJOumk8LfRmDFjNG3aNO3fv1+nnXZazfHp06dr7NixSklJ0eHDh/Xqq6+qoKCg3l59a9asUVZWVkSfbeHChfL5fPrnf/5nXXLJJTVh9d1339Xjjz+uwYMHa86cOXX6/PGPf9Qvf/lLJSQkaNy4cUpLS9O3336rdevWaf78+fr5z39e5/wvv/xSGRkZGjBggH7xi1+opKREf/7zn5WTk6PVq1frsssuqzn3hRde0D333KMOHTpo3LhxSklJ0RdffFGzhcUnn3yi1NTUJtVT29GjR3XrrbfqzTff1OTJkzVnzhy1a8e/73hVQVKq9ETwFXUZCQTQFNV7/AVT1r57FKsB0FiEvjhVHfqKi4tVUFCgSy65pKZtypQpmjt3rubMmaMnn3yyJvSde+65da5hZiopKVG3bt0kSS+++KLuvPNO7dy5U/37969zbnp6ujIyMpSfn19n24W9e/fq1ltv1d69e9W1a1cNHTpUK1as0JgxY2rOOXr0qBYvXqyVK1dG9Nl69OihRYsWadSoUbrxxhu1fv16/fDDD7r11lvVoUMHvfbaa3X2FtyyZYvuu+8+denSRWvWrNF5551X53p79uyp9x4FBQV65JFHlJubW3Ns/PjxysrK0u9+97ua0Ldt2zZNmjRJ/fv311/+8hf17fs/e0u99957uvLKK/XAAw9o8eLFTaqnWklJiXJycvTXv/5Vs2bN0rRp0yL6M0MbllN/RLuOvJzQ7WjzWDAFLSncHn8A2gZCX5yqDn3z5s2rE/gk6e6779bcuXNrpmPOnDlTM2fODHvNnTt3avDgwerXr1/A9tzcXD3wwAO69957lZCQIKkqKIaTl5enESNG6MILLwx7brWLLrpIv/3tb/XQQw9p0qRJ2rdvn/bu3as//vGP9ULUggULdOzYMf3mN7+p1yYp4Oc544wzNGNG3f8jHDNmjFJTU/Xpp5/WuXZFRYWefvrpOoFPki6//HKNGzdOb7/9tg4fPqzOnTs3uh6pavQxKytLO3bs0Msvv6xbbrklyJ8OGiPUs3WMpgEAgNaM0BeHysvLtWXLFqWkpAQMBtXTLysqKhp03eXLl2vevHlBp3tmZWVp8uTJ2rNnj84444yIr9u+ffuAz+eFM23aNBUUFOjVV1+VJN18880Bn+X75JNPJFVtPxGpYcOG1QTX2lJSUvTxxx/XvK7++S9/+YvWrl1b7/xvv/1Wx48f17Zt23TBBRc0up6tW7dq5MiRKisr04oVK3TFFVdE3BcRqn62LhBG0wAAQCtG6ItDGzduVEVFhcaMGRPwWa9du3ZJUr3nzMIJFGpONGXKlAZdU5ImTpzY4D5S1fTTf/zHf6yZFvqrX/0q4HkHDx6UpHojcaFUT2k90UknnaTKysqa1/v375ck/e53vwt5vdLS0ibVs23bNpWUlGjYsGH6yU9+EnE/eF+oZ/4YoQSaF1NtAbRWhL44VD2188Tn7qpVL6IyevToaJXUIr744gv9+te/Vvfu3XXo0CHdfffd+vTTT9WxY8c651UHuK+++kqDBg1q1hqq9zE8dOiQunTpElGfxtRz9dVX65xzztG//uu/6oorrtA777yjpKSkxhUNbwn1zB8jlAAAxAWW9ItD1aGvekSptpKSEj333HPq16+fcnKq/kL4wQcfaNy4cerbt6/MLKLn8GKtvLxcN954o8rKypSfn6+HHnpIGzduDDjaV/2s4IoVK5q9juprr1mzpsF9GlrPQw89pNmzZ6uoqEiXXXaZvvnmmwb1BxB/lmz4KuQXAMAbCH1xqHo1zjfeeENlZWU1x0tLSzV+/Hjt379fTz/9dM2IWGlpqYYMGaKnn35aJ598ckxqbqhf//rXKioq0r/8y79o9OjR+rd/+zddfPHF+sMf/lBvn7tf/vKXOumkk/Tb3/424F6CoVbLDOf+++9X+/btNXXqVG3btq1e+48//lgvEDalnl/96ldasGCBNm/erFGjRoXcBxAAAADxgemdcebYsWPatGmTzj//fB06dEhDhw7VuHHjVF5eriVLlqi4uFgzZ87UtddeW9MnOztb2dnZkqQ77rgjRpVH7q233tK8efM0YsQIPfroo5KkhIQELVq0SMOGDdM999yj4cOHa8CAAZKkwYMHa/78+br33nvl8/mUk5OjtLQ07d+/X+vWrVPnzp3rbFbfEIMGDdLChQs1YcIEnXfeecrKytLAgQNVUVGh3bt3a82aNUpOTtZ///d/1/Rpaj333nuvOnbsqLvuukuXXnqp3nvvvQY/nwkAAADvIPTFmc2bN+vo0aMaOXKkHggwT38AABQ9SURBVHzwQU2ZMkV5eXlyzmnEiBF64YUX2vSzfLt379aECRPUtWtXLVq0qM5KoikpKVq4cKGuueYa3XTTTfrwww+VmJgoSbrnnns0ZMgQPfnkkyooKNBbb72lpKQkDR06NOCKnw1x66236vzzz9fvf/97vf/++3rnnXfUqVMn9enTR9dff71uvPHGen2aWs8dd9yhDh066LbbbqsJftUhF4GF2pJBktQpOXrFREm4jd29+JkBAIhHhL44U/08n8/nU1paWos8xxZLqampKikpCdqek5Mj51zAtpEjR+rNN98Mef3+/fsH7S9VbdoeSHp6eoOfhWxqPTfffLNuvvnmBr1nXAu1JYNXhdvYHQAAeAKhr5bMlMxYl9Diaoc+AAAAAN5H6Isz69evV0JCgtLT02NdCgAgQqykCQBoCkJfHKmsrNRnn32mQYMG1durDgDQNKGCWc6wvlGsBACAugh9caRdu3YqLS1tcL/S0lJt375dUlVw3L17tzZs2KAePXqwKiQATyG4AQC8iH36ENa6devk8/nk8/l05MgR5ebmyufz6eGHH451aQAAAADCYKQPYWVmZoZcsRIAACCQZV36yLf+0oBtZe27R7kaIH4R+gAAANAiVqfNiHUJAEToAwAg5lidEwDQkgh9AADPCBeeWIwFABCPWMgFAAAAADyMkT4AcaNgXrpUti9wY6fk6BYDAAAQJYQ+APGjbJ80Pj/WVQAA4BmzV22LdQmIQNxN72TrASAw/rcBAADgTXEV+hISElRRURHrMoBWqaKiQgkJCbEuAwAAAM0srqZ3du7cWd9//72SkpJiXQrQ6nz//ffq3LlzrMsAPIktGQAAsRRXI309evTQgQMH9N133+nHH39kOhvinnNOP/74o7777jsdOHBAPXr0iHVJAAAAaGZxNdLXoUMHpaamqqSkRLt27dLx48djXRIQcwkJCercubNSU1PVoUOHWJcDAACAZhZXoU+qCn69e/dW7969Y10KAKANYYomgHjFCp1tX9yFPgAAALR+y7r0kW/9pQHbfJXdtHD40ihXBLRdhD4AAAC0OqvTZgRtyyy8L4qVAG0foQ8AYqwgKVV6olfgxk7J0S0GAAB4DqEPAGItZ26sKwAAICZCPS849cqBUazE2+JqywYAAAAAiDeEPgAAAADwMEIfAAAAAHgYz/QBQJwKt+9czrC+UaqkYdgvDwCAhiH0AQBaFUIdAADNi+mdAAAAAOBhEYU+M8sys61mtt3MpgdoNzOb42//3Mx+Eq6vmfUws1Vm9oX/e/cTrplqZqVm9uumfEAAAAAAiGdhp3eaWYKkZyRdKWmPpLVmttQ5t6XWaVdJSvN/jZC0QNKIMH2nS3rXOTfLHwanS5pW65qzJa1o6gcEAABAfBm4MUedKg4EbfdVdtPC4UujWBEQW5E805chabtz7m+SZGb5knIk1Q59OZL+5Jxzkj4xs25m1ltS/xB9cyRl+vu/JKlA/tBnZtdI+puksiZ8NgAAAMShThUHtGTw74O2ZxbeF8VqgNiLJPT1lfT3Wq/3qGo0L9w5fcP07emc+1qSnHNfm9npkmRmnVQV/q6UxNROAPCgWC3WwiIxAIB4FEnoswDHXITnRNL3RP8mabZzrtQsUHf/G5pNlDRRklJTU8NcEkA8KJiXLpXtC35Cp+ToFQMAANBKRBL69khKqfW6n6TiCM9JDNH3GzPr7R/l6y3pW//xEZKuN7N/l9RNUqWZHXXOzav9hs655yQ9J0nDhw8PFyQBNKNw4SozMVmaujGKFfmV7ZPG50f/fQEAAFqxSELfWklpZnampK8k3SRp/AnnLJV0v/+ZvRGSDvnD3L4QfZdKul3SLP/3JZLknPtZ9UXN7BFJpScGPgAxFi5c5eUEbWq1gREAAMCjwoY+59wxM7tf0kpJCZIWOuc2m9m9/vZnJS2XlC1pu6QfJN0Zqq//0rMkvWZmd0naLemGZv1kAFqnJgTG1qogKVV6olfwE5hWCgAAYiiSkT4555arKtjVPvZsrZ+dpMmR9vUf3y/pijDv+0gk9QFATOXMjXUFrQ4LpgAA0HpEtDk7AAAAAKBtimikDwCA2hjJAwCg7WCkDwAAAAA8jJE+AM0u5MImLGoCAGiiZV36yLf+0qDtZe27R7EaoPUj9AFofixsAgBoQavTZsS6BETB7FXbQrZPvXJglCpp+5jeCQAAAAAexkgfgFYl1NRQNm4HAABoOEIfgNYl1NTQNrhxOwAAQKwR+gC0GSEXiJFYJAYAACAAQh+AtiMOF4gJtx9ezrC+UaoEAAC0VSzkAgAAAAAeRugDAAAAAA8j9AEAAACAh/FMHwC0YTzzBwANt6xLH/nWXxqwzVfZTQuHL41yRUDLIvQBAAAgrqxOmxG0LbPwvihWAkQH0zsBAAAAwMMIfQAAAADgYUzvBAAPC/fMHwAA8D5G+gAAAADAwwh9AAAAAOBhTO8EAATE1FAAALyB0AcATcReeQAAoDVjeicAAAAAeBihDwAAAAA8jNAHAAAAAB5G6AMAAAAADyP0AQAAAICHsXonEKcK5qVLZfsCtmUmJktTN0a5IgAAALQEQh8Qr8r2SePzA7fl5US3FgAAALQYpncCAAAAgIcR+gAAAADAwwh9AAAAAOBhhD4AAAAA8DAWcgFQT0FSqvREr+AndEqOXjEAAABoEkIfgPpy5sa6AgAAADQTQh/QioXaS09iPz0AAACER+gDWrNQe+lJ7KcHAACAsAh9QAyFG8kL9+xcqGfvGAUEAACAROgDYivcSF44oZ69a6WjgEs2fBWyPWdY35i8d0u+bzjh/kwAAACagtAHeFQ8rsDZkuGJYAYAANoqQh/gVR5cgZPgBQAAqs1etS1o29QrB0axktaPzdkBAAAAwMMY6QMAMYoIAAC8i9AHAC2MQAkAAGKJ0AcAAAD4LevSR771lwZt91V208LhS6NYEdB0hD4AAADAb3XajJDtmYX3RakSoPkQ+gA0WCz32gMAAEDDEPoAtCo8/wYAANC82LIBAAAAADyM0AcAAAAAHsb0TqCFFcxLl8r2BW7slBzdYgAAABB3Igp9ZpYl6WlJCZKed87NOqHd/O3Zkn6QdIdzbn2ovmbWQ9KfJfWXtEvSz51zB8zsSkmzJCVK+lHSPzvn3mvaxwRiqGyfND4/1lUAAAAgToWd3mlmCZKekXSVpMGSbjazwSecdpWkNP/XREkLIug7XdK7zrk0Se/6X0vSd5Kuds6lS7pd0suN/nQAAAAAEOciGenLkLTdOfc3STKzfEk5krbUOidH0p+cc07SJ2bWzcx6q2oUL1jfHEmZ/v4vSSqQNM05V1TrupsldTSzDs658kZ9QgBRxwqcAAAArUckC7n0lfT3Wq/3+I9Fck6ovj2dc19Lkv/76QHe+zpJRQQ+AAAAAGicSEb6LMAxF+E5kfQN/KZm50l6QtLoIO0TVTWVVKmpqZFcEgAAAADiTiQjfXskpdR63U9ScYTnhOr7jX8KqPzfv60+ycz6SVos6Tbn3I5ARTnnnnPODXfODU9OZgVEAAAAAAgkktC3VlKamZ1pZomSbpK09IRzlkq6zapcKOmQf8pmqL5LVbVQi/zfl0iSmXWTtEzSQ865vzbhswEAAABA3As7vdM5d8zM7pe0UlXbLix0zm02s3v97c9KWq6q7Rq2q2rLhjtD9fVfepak18zsLkm7Jd3gP36/pLMl/cbMfuM/Nto5VzMSCAAAAACITET79Dnnlqsq2NU+9mytn52kyZH29R/fL+mKAMcflfRoJHUBAAAAAEKLKPQB8J5Q2yrkDDtxgV4AAAC0VZE80wcAAAAAaKMIfQAAAADgYUzvBNowpmgCAAAgHEIfIKlgXrpUti9oe2ZisjR1Y6P6qhP7SAIA4BXLuvSRb/2lQdt9ld20cPiJu5sBsUXoA6Sq0DY+P3h7Xk7j+7ZBoUYQAQCIZ6vTZoRszyy8L0qVAJEj9AExFC5cMUUTAAAATcVCLgAAAADgYYQ+AAAAAPAwQh8AAAAAeBihDwAAAAA8jIVcgAgUJKVKT/QK3MiWDAAAAGjFCH1AJHLmxroCAAAAoFEIfYgL8biBOnvtAQAAQCL0wUNCBrtOyZ7bQB0AAACIBKEP3lG2j2AHAAAAnIDQB7RiTNEEAABAU7FlAwAAAAB4GKEPAAAAADyM0AcAAAAAHkboAwAAAAAPI/QBAAAAgIexeieg8Ktk5gzr22LXBgAAQPOavWpbyPapVw6MUiWtAyN9AAAAAOBhhD4AAAAA8DBCHwAAAAB4GM/0AU3EM3sAAABozRjpAwAAAAAPY6QPUVUwL10q2xewLTMxWZq6McoVAQAAAN5G6EN0le2TxucHbsvLiW4tAAAAQBwg9CHqgj0DlxndMgAAAIC4wDN9AAAAAOBhhD4AAAAA8DBCHwAAAAB4GM/0oc0ItfKnJKlTcvSKAQAAANoIQh+a1exV20K2+5py8bJ9WjL496HPCbJITM6wvk15ZwAAAKDNYnonAAAAAHgYI32oJ9xo3dQrB8bkvZs0SggAAADEKUb6AAAAAMDDCH0AAAAA4GFM70RcWBJkgZdo9QcAAABihZE+AAAAAPAwRvrQYKH2y5tQ2U0Lhy9t1HWXdekj3/pLg7aXte/eqOsCAAAA8YzQh4Yr2yeNzw/Y1DUvp9GXXZ02o9F9AQAAAATG9E4AAAAA8DBCHwAAAAB4GKEPAAAAADyM0AcAAAAAHsZCLqhn4MYcdao4ELS9rH13rQ6yb10ZK3ACAIA4Fmo1cl8TVjkHmoLQF4eWPXVuyFCn9t21ZPDvG3VtVuAEAADxLNTfhTIL74tiJcD/IPR51OxV24K2+SoONDrUAQAAAGhbInqmz8yyzGyrmW03s+kB2s3M5vjbPzezn4Tra2Y9zGyVmX3h/969VttD/vO3mtmYpn5IAAAAAIhXYUOfmSVIekbSVZIGS7rZzAafcNpVktL8XxMlLYig73RJ7zrn0iS9638tf/tNks6TlCVpvv86AAAAAIAGimSkL0PSdufc35xzP0rKl5Rzwjk5kv7kqnwiqZuZ9Q7TN0fSS/6fX5J0Ta3j+c65cufcTknb/dcBAAAAADRQJM/09ZX091qv90gaEcE5fcP07emc+1qSnHNfm9npta71SYBrAQAAAECLC7U+xtQrB0axkuYRSeizAMdchOdE0rcx7yczm6iqqaSSVGpmW8NcNxaSJH0X6yICuzXWBaBhWvG9hDaE+wjNgfsIzSXu7qX/LUk6J8ZVeE6j7qMHm/CGTenbws4I1hBJ6NsjKaXW636SiiM8JzFE32/MrLd/lK+3pG8b8H5yzj0n6bkI6o8ZM1vnnBse6zrQ9nEvoTlwH6E5cB+huXAvoTlwH0Umkmf61kpKM7MzzSxRVYusnLir5FJJt/lX8bxQ0iH/1M1QfZdKut3/8+2SltQ6fpOZdTCzM1W1OMynjfx8AAAAABDXwo70OeeOmdn9klZKSpC00Dm32czu9bc/K2m5pGxVLbryg6Q7Q/X1X3qWpNfM7C5JuyXd4O+z2cxek7RF0jFJk51zx5vrAwMAAABAPDHnwj1ih8Yys4n+aahAk3AvoTlwH6E5cB+huXAvoTlwH0WG0AcAAAAAHhbJM30AAAAAgDaK0NdCzCzLzLaa2XYzmx7retC6mdkuM9toZhvMbJ3/WA8zW2VmX/i/d691/kP+e2urmY2JXeWIJTNbaGbfmtmmWscafN+Y2QX++2+7mc0xs0Bb58DDgtxLj5jZV/7fSxvMLLtWG/cS6jGzFDN738z+y8w2m9kD/uP8XkLEQtxH/E5qAkJfCzCzBEnPSLpK0mBJN5vZ4NhWhTbgMufcsFrLDk+X9K5zLk3Su/7X8t9LN0k6T1KWpPn+ew7x50VV3QO1Nea+WaCqfU/T/F8nXhPe96IC/3ef7f+9NMw5t1ziXkJIxyT9L+fcuZIulDTZf7/wewkNEew+kvid1GiEvpaRIWm7c+5vzrkfJeVLyolxTWh7ciS95P/5JUnX1Dqe75wrd87tVNWquRkxqA8x5pz7QFLJCYcbdN9Y1T6pXZxzH7uqh7z/VKsP4kSQeykY7iUE5Jz72jm33v/zYUn/Jamv+L2EBghxHwXDfRQBQl/L6Cvp77Ve71HomxVwkt4xs0Izm+g/1tO/36X830/3H+f+QigNvW/6+n8+8TggSfeb2ef+6Z/VU/K4lxCWmfWX5JP0n+L3EhrphPtI4ndSoxH6Wkag+cIsk4pQLnbO/URVU4Inm9mlIc7l/kJjBLtvuJ8QzAJJZ0kaJulrSb/3H+deQkhmdqqkNyX9yjn3fahTAxzjXoKkgPcRv5OagNDXMvZISqn1up+k4hjVgjbAOVfs//6tpMWqmq75jX9qgvzfv/Wfzv2FUBp63+zx/3ziccQ559w3zrnjzrlKSX/U/0wj515CUGbWXlV/Uf8/zrn/6z/M7yU0SKD7iN9JTUPoaxlrJaWZ2Zlmlqiqh0uXxrgmtFJm1snMOlf/LGm0pE2qumdu9592u6Ql/p+XSrrJzDqY2ZmqejD50+hWjVasQfeNf6rVYTO70L+q2W21+iCOVf8l3e8fVfV7SeJeQhD+/+55kv7LOfdUrSZ+LyFiwe4jfic1zUmxLsCLnHPHzOx+SSslJUha6JzbHOOy0Hr1lLTYv4rwSZJedc79PzNbK+k1M7tL0m5JN0iSc26zmb0maYuqVria7Jw7HpvSEUtmtkhSpqQkM9sjKVfSLDX8vvmlqlZvPFnSCv8X4kiQeynTzIapajrULkmTJO4lhHSxpF9I2mhmG/zH/lX8XkLDBLuPbuZ3UuNZ1WI2AAAAAAAvYnonAAAAAHgYoQ8AAAAAPIzQBwAAAAAeRugDAAAAAA8j9AEAAACAhxH6AAAAAMDDCH0AAAAA4GGEPgAAAADwsP8PSUmagVYIAA8AAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 1080x360 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"plt.hist(p1_Bbar[:,0],alpha=0.5, density=True,label=r'$p_{1}^{E, (2)}$', range=(0,2600),bins=100);\n",
"plt.hist(p1[:,0],histtype=\"step\", density=True,label=r'$p_{1}^{E, (3)}$', range=(0,2600), bins=100);\n",
"plt.hist(p1_q[:,0],alpha=0.3,label=r'$p_{1}^{E, (3)}$ xcheck', density=True,range=(0,2600), bins=100);\n",
"plt.legend(fontsize='20');\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(15,5)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABIcAAAGbCAYAAABeXfDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjAsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+17YcXAAAgAElEQVR4nOzdeZxdVYEv+t8iJOnIjGEIgRBUhCtDQxuDXlSMASEiIt1ox+nStC0PG6/QXu+H4T2koRVD67Nbu0WbhwjefsLVblE/QjM0EQWnhjA7NA8VJcxDIkEChLDeH6nEqtSp1EnqVNU5tb/fzyef1Fl7qHXOOvvsU7+99lql1hoAAAAAmmmz8a4AAAAAAONHOAQAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoME2H+8KtDJ9+vQ6e/bs8a4GAAAAwISxZMmSx2qtO6xf3pXh0OzZs3PzzTePdzUAAAAAJoxSyq9blbutDAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA02ObjXQEAAABgsGeffTZPPPFEVqxYkdWrV493degykyZNylZbbZXtt98+U6dOHdG+hEMAAADQZZ599tn85je/yXbbbZfZs2dn8uTJKaWMd7XoErXWrFq1Kk8++WR+85vfZNasWSMKiNxWBgAAAF3miSeeyHbbbZfp06dnypQpgiEGKKVkypQpmT59erbbbrs88cQTI9qfcAgAAAC6zIoVK7L11luPdzXoAVtvvXVWrFgxon0IhwAAAKDLrF69OpMnTx7vatADJk+ePOIxqYRDAAAA0IXcSkY7OvE+EQ4BAAAANNiws5WVUi5K8pYkj9Ra922x/H8meXe//f2XJDvUWp8opdybZEWS1Umer7XO6VTFAQAAABi5dqayvzjJPyb5cquFtdZPJvlkkpRSjkryV7XW/sNkz6u1PjbCegIAAACMqp8/+GSeW/3CgLIpkzbL3jMm9uDgw4ZDtdbvlVJmt7m/dya5dCQVAgAAABgPz61+Ifvvuu2AsjuWLh+n2oydjo05VEp5UZIjkvxrv+Ka5JpSypJSygnDbH9CKeXmUsrNjz76aKeqBQAAAMAGdHJA6qOSfH+9W8oOrrX+UZIFSU4qpbx+qI1rrRfUWufUWufssMMOHawWAAAAwNhbtmxZdtppp/ziF79oe5tjjz02n/70p0exVoN1MhxamPVuKau1PtD3/yNJLk8yt4O/DwAAAKBrnXvuuXnzm9+cl770pUmST3ziE3nVq16VrbfeOjvssEOOOuqo3HXXXQO2Oeuss/Kxj30sv/3tb8esnh0Jh0op2yQ5JMk3+5VtUUrZau3PSd6U5K7WewAAAACYOJ5++ulceOGFed/73reu7Prrr89f/uVf5gc/+EEWL16czTffPIceemieeOL3N2Htt99+eclLXpJ//ud/HrO6tjOV/aVJ3pBkeillaZKzkkxOklrrF/pWOybJNbXW3/XbdKckl5dS1v6er9Rar+pc1QEAAAC605VXXpnNNtssBx988Lqyq6++esA6/+t//a9ss802+f73v5+jjjpqXflb3/rWXHrppTnppJPGpK7tzFb2zjbWuThrprzvX/bLJH+4qRUDAAAA6FU33HBDXvnKV6av00xLK1asyAsvvJDttttuQPncuXPzsY99LCtXrsy0adNGu6odHXMIAAAAYMJauXJlzjnnnOy99975gz/4g+y2224544wzsmrVqkHr/vrXv86MGTM2uL+TTz45BxxwQF7zmtcMKN9ll12yatWqPPDAAx2t/1CG7TkEAAAAdJeDFy3O/ctXjnc1hjRz22n5/mlvHO9qdNSDDz6Yww47LHfffXeOOeaYHH300fn2t7+dT3ziE3n88cfzT//0TwPWX7lyZXbaaach9/fhD384N954Y2688cZMmjRpwLK1vYVWrhybNhYOAQAAQI+5f/nK3LvoyPGuxpBmn3bFeFeho5577rkcddRRuffee/Od73xn3ThCZ555ZvbZZ59ceOGFOfvss7Pzzjuv22b69OlZtmxZy/391V/9VS677LJ85zvfyUte8pJBy9cOUL3DDjuMwrMZzG1lAAAAABvwqU99KkuWLMl55503YIDpLbfcMsccc0xeeOGF3HDDDQO2OfDAA/PTn/500L5OPvnkfOUrX8nixYuz9957t/x9d911V3bZZZcN9jzqJOEQAAAA0PWuv/76/Omf/ml23XXXTJ06NTNmzMjhhx+eb3zjG+vWqbXm4osvziGHHJIXv/jFmTZtWg488MBccsklQ+73hhtuyB//8R/npS99aV71sp2z4447Zu7cuTnjjDOSJM+sXJlPfvKTmTFjRk444YRB27/4xS9Okjz00EMDyg8//PD87Gc/y+OPP76u7KSTTsqXvvSlXHrppdluu+3y0EMP5aGHHspTTz01qE5HHHHExr9Im0g4BAAAAHS1U045JfPmzcvixYtz6KGH5sMf/nDmzZuXJUuW5Ac/+EGSNePzHHHEETn++OOzfPnyHHfccTn++OPz0EMP5c/+7M9y7rnnDtrvueeem9e//vVZsmRJ5s+fn/f+xV/mqKOOyrPPPpurrroqSbL4qm9n+fLlede73pXJkycP2sczzzyTJJkyZcqA8v322y9z587NZZddtq7s/PPPz4oVKzJ//vzMmDFj3b9PfepTA/Z3+eWX5/3vf//IX7g2GXMIAAAA6FpnnHFGPvOZz+RP/uRPcskll2SLLbZYt+ypp55a1zPnXe96V6655pp8/OMfX9frJ0nOPvvs7L333jnnnHPygQ98YN208Q8//HA++tGP5rWvfW2uu+66TJkyJXcsXZ79d902SfLYY48lSW5YfE2S5P77789f//VfD6rfv//7vydJdtttt0HLzjrrrJx88sk58cQTM2nSpNRah32+X/ziF3PQQQfl1a9+dTsvT0cIh0ZRq9HjJ+KI7QAAADAabr311px33nmZM2dOvvKVrwzqnbPllltmyy23zLe//e184xvfyNvf/vYBwVCyZlDnt7zlLfnyl7+cW265JfPnz0+S/PznP8/q1avz8pe/fNB+kzUDSifJrTf9KEkG9ABq5RWveMWgsiOOOCInnXRSli5dmt13372t5zx58uT8wz/8Q1vrdopwaBS1Gj1+oo3YDgAAAKPl05/+dF544YUsWrSoZYCz1gUXXJAkOfXUU1suXzsu0OrVq9eV7bPPPtlmm21y0UUX5dFHH8273/3uzNznoKSv51CS/O53v8uD9y/NvvvumzvvvHPQfp966qlsv/322WmnnTJ79uyWv/tDH/rQsM+zv1bjGo02Yw4BAAAAXenqq6/Odtttl3nz5m1wve9+97vZeuut88pXvrLl8gcffDBJMmvWrHVl06dPz4033phjjz021113XRYuXJg3HPCyLFiwILfcckuSNbeSJckuu+wyZP1WrVqVI488suXyXiEcAgAAALrOM888k0cffTS77757Ntts6PhixYoVefLJJ4fsubN69ep897vfzY477pi99tprwLJ99903X/va17Js2bJce+21mb/gqFx11VV505velGeffTbPPfdckmTq1Kkt9/2lL30pSfLnf/7nA8pPPfXUHHbYYe0+1XHntjIAAACg66wdvPmRRx7Z4HpTp07NZpttlmXLlrVcfvHFF+fBBx/MqaeemlJKy3WmTJmSQw89NDvuPSdPL38sN954Yx5++OHsvPPOSQZPU58kP/rRj3LllVdmwYIFmTt37oBlt912Ww444IBhn2O30HMIAAAA6DrTpk3LvvvumwceeCBf/epXBy2/++67s3r16kyZMiWvetWrct99962bOWyt6667LieffHJmz56d008/fV35rbfeml/84heD9vmbX/0yd911V2bNmpVdd90106dPz0v23CtLlizJHXfcsW69X//611m4cGG22WabnH/++YP2c/vtt+fAAw8cydMfU3oOAQAAAF1p0aJFeetb35qFCxfmkksuyT777JPly5fntttuy3333bduLKGPf/zjOfzww/OWt7wl73jHO7LLLrvk9ttvz9VXX53dd9891157bbbZZpt1+/3sZz+bSy65JHPnzs0+++yTHXfcMb/61a/yzW9+K6UkF1100bpb2d7/3/9HTv/QCZk/f37e85735He/+12++tWvppSSK664YtDtbA899FAefvjhnuo5JBwCAAAAutKRRx6Z66+/PosWLcoPf/jDXHPNNZk+fXr222+/nHzyyevWmz9/fq677rqcffbZ+frXv54k2WOPPXLmmWfmIx/5SLbaaqsB+z366KPz/PPP5z/+4z/yta99Lc8880x22WWXLHjbsTnvnDOz5557rlv3zce8PTO2npK//du/zec///lMnz4973jHO3LWWWdl5syZg+p86623Ztq0aYPGN+pmwiEAAADoMTO3nZbZp10x3tUY0sxtp3VsX6973evyute9btj1DjnkkCxevLitfb7tbW/L2972tkHldyxdnj37TWW/1nHHHZfjjjuurX3fdttt2W+//TJp0qS21u8GwiEAAADoMd8/7Y3jXQWG0GuDUScGpAYAAADoGOEQAAAAQEM9/fTTueeee3pqprJEOAQAAADQEbfffnuSZP/99x/nmmwc4RAAAABAB9x+++3Zc88986IXvWi8q7JRhEMAAAAAHXDiiSfm5z//+XhXY6MJhwAAAAAaTDgEAAAA0GDCIQAAAIAGEw4BAAAANJhwCAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoMGEQwAAAACjYNmyZdlpp53yi1/8ou1tjj322Hz6058exVoNJhwCAAAAGAXnnntu3vzmN+elL31pkuRzn/tc9t9//2y99dbZeuut85rXvCZXXHHFgG3OOuusfOxjH8tvf/vbMauncAgAAACgw55++ulceOGFed/73reubNddd815552XW265JTfffHPe+MY35m1ve1vuuOOOdevst99+eclLXpJ//ud/HrO6CocAAAAAOuzKK6/MZpttloMPPnhd2dFHH50FCxbkZS97WV7+8pfn4x//eLbaaqv88Ic/HLDtW9/61lx66aVjVlfhEAAAAECH3XDDDXnlK1+ZUkrL5atXr85ll12Wp556Kv/1v/7XAcvmzp2b//iP/8jKlSvHoqrCIQAAAIB2rFy5Muecc0723nvv/MEf/EF22223nHHGGVm1atWgdX/9619nxowZg8rvvPPObLnllpk6dWpOPPHEXH755dlvv/0GrLPLLrtk1apVeeCBB0btufS3+Zj8FgAAAKBzvvOJ8a7Bhs07fbxr0HEPPvhgDjvssNx999055phjcvTRR+fb3/52PvGJT+Txxx/PP/3TPw1Yf+XKldlpp50G7WevvfbKbbfdluXLl+df//Vfc9xxx+X666/Pvvvuu26dadOmrdvHWNBzCAAAAGADnnvuuRx11FG59957853vfCf/+3//75x33nn58Y9/nFmzZuXCCy/MQw89NGCb6dOnZ9myZYP2NWXKlLzsZS/LnDlz8olPfCIHHHBA/u7v/m7AOk888USSZIcddhi9J9XPsOFQKeWiUsojpZS7hlj+hlLKb0spt/X9+2i/ZUeUUv6zlHJPKeW0TlYcAAAAYCx86lOfypIlS3LeeecNGGB6yy23zDHHHJMXXnghN9xww4BtDjzwwPz0pz8ddt8vvPBCnn322QFld911V3bZZZeWPY9GQzs9hy5OcsQw69xQaz2g7985SVJKmZTkc0kWJHlFkneWUl4xksoCAAAAzXT99dfnT//0T7Prrrtm6tSpmTFjRg4//PB84xvfWLdOrTUXX3xxDjnkkLz4xS/OtGnTcuCBB+aSSy4Zcr833HBD/viP/zgvfelL86qX7Zwdd9wxc+fOzRlnnJEkeWblynzyk5/MjBkzcsIJJwza/sUvfnGSDOo5dPjhh+dnP/tZHn/88XVlp512Wm644Ybce++9ufPOO3P66afn+uuvz7vf/e5BdTriiOGimM4ZNhyqtX4vyRObsO+5Se6ptf6y1vpcksuSHL0J+wEAAAAa7JRTTsm8efOyePHiHHroofnwhz+cefPmZcmSJfnBD36QZM34PEcccUSOP/74LF++PMcdd1yOP/74PPTQQ/mzP/uznHvuuYP2e+655+b1r399lixZkvnz5+e9f/GXOeqoo/Lss8/mqquuSpIsvurbWb58ed71rndl8uTJg/bxzDPPJFlzu1h/++23X+bOnZvLLrtsXdlDDz2U97znPdlrr70yf/783HTTTfm3f/u3LFiwYMD+Lr/88rz//e8f+QvXpk4NSP2aUsrtSR5I8pFa60+SzExyX791liY5aKgdlFJOSHJCksyaNatD1QIAAAB62RlnnJHPfOYz+ZM/+ZNccskl2WKLLdYte+qpp9b1zHnXu96Va665Jh//+MfX9fpJkrPPPjt77713zjnnnHzgAx/IdtttlyR5+OGH89GPfjSvfe1rc91112XKlCm5Y+ny7L/rtkmSxx57LElyw+JrkiT3339//vqv/3pQ/f793/89SbLbbrsNWnbWWWfl5JNPzoknnphJkybl4osvHvb5fvGLX8xBBx2UV7/61W28Op3RiXDoliS711qfKqW8Ock3kuyZpLRYtw61k1rrBUkuSJI5c+YMuR4AAADQDLfeemvOO++8zJkzJ1/5ylcG9c7Zcssts+WWW+bb3/52vvGNb+Ttb3/7gGAoWTOo81ve8pZ8+ctfzi233JL58+cnSX7+859n9erVefnLXz5ov8maAaWT5NabfpQkA3oAtfKKVwweSeeII47ISSedlKVLl2b33Xdv6zlPnjw5//AP/9DWup0y4nCo1vpkv5+vLKWcX0qZnjU9hfrHZrtmTc8iAAAAgGF9+tOfzgsvvJBFixa1DHDWuuCCC5Ikp556asvla8cFWr169bqyffbZJ9tss00uuuiiPProo3n3u9+dmfsclPT1HEqS3/3ud3nw/qXZd999c+eddw7a71NPPZXtt98+O+20U2bPnt3yd3/oQx8a9nn212pco9E24qnsSyk7l1JK389z+/b5eJKbkuxZStmjlDIlycIk3xrp7wMAAACa4eqrr852222XefPmbXC97373u9l6663zyle+suXyBx98MMnAYWymT5+eG2+8Mccee2yuu+66LFy4MG844GVZsGBBbrnlliRrbiVLkl122WXI+q1atSpHHnnkRj+3btLOVPaXJvlhkr1KKUtLKe8rpZxYSjmxb5Vjk9zVN+bQZ5MsrGs8n+SDSa5O8rMkX+0biwgAAABgg5555pk8+uij2X333bPZZkPHFytWrMiTTz45ZM+d1atX57vf/W523HHH7LXXXgOW7bvvvvna176WZcuW5dprr838BUflqquuypve9KY8++yzee6555IkU6dObbnvL33pS0mSP//zPx9Qfuqpp+awww5r96mOu2FvK6u1vnOY5f+Y5B+HWHZlkis3rWoAAABAU9W6ZjjiRx55ZIPrTZ06NZtttlmWLVvWcvnFF1+cBx98MKeeemr6bnwaZMqUKTn00EOz495z8vTyx3LjjTfm4Ycfzs4775xk8DT1SfKjH/0oV155ZRYsWJC5c+cOWHbbbbflgAMOGPY5dosR31YGAAAA0GnTpk3LvvvumwceeCBf/epXBy2/++67s3r16kyZMiWvetWrct99962bOWyt6667LieffHJmz56d008/fV35rbfeml/84heD9vmbX/0yd911V2bNmpVdd90106dPz0v23CtLlizJHXfcsW69X//611m4cGG22WabnH/++YP2c/vtt+fAAw8cydMfU52ayh4AAACgoxYtWpS3vvWtWbhwYS655JLss88+Wb58eW677bbcd99968YS+vjHP57DDz88b3nLW/KOd7wju+yyS26//fZcffXV2X333XPttddmm222Wbffz372s7nkkksyd+7c7LPPPtlxxx3zq1/9Kt/85rdSSnLRRRetu5Xt/f/9f+T0D52Q+fPn5z3veU9+97vf5atf/WpKKbniiisG3c720EMP5eGHH+6pnkPCIQAAAKArHXnkkbn++uuzaNGi/PCHP8w111yT6dOnZ7/99svJJ5+8br358+fnuuuuy9lnn52vf/3rSZI99tgjZ555Zj7ykY9kq622GrDfo48+Oo+vWJk7bl2Su776tTz77DPZcaedc+Qxx+YTZ5+ZPffcc926bz7m7Zmx9ZT87d/+bT7/+c9n+vTpecc73pGzzjorM2fOHFTnW2+9NdOmTRs0vlE3Ew4BAABAr5l3+vDrTBCve93r8rrXvW7Y9Q455JAsXry4rX2+7W1vy0vmvCH795u2fkOOO+64HHfccW2te9ttt2W//fbLpEmT2lq/GxhzCAAAAKBDem0w6kQ4BAAAANAxwiEAAACAhnr66adzzz339NRMZYlwCAAAAKAjbr/99iTJ/vvvP8412TgGpB5jM7edltmnXTGo7PunvXGcagQAAAB0wu23354999wzL3rRi8a7KhtFODTGWoVA64dFAAAAQO858cQTc+KJJ453NTaa28oAAAAAGkw4BAAAANBgwiEAAACABhMOAQAAADSYcAgAAACgwYRDAAAA0IVqreNdBXpAJ94nwiEAAADoMpMmTcqqVavGuxr0gFWrVmXSpEkj2odwCAAAALrMVlttlSeffHK8q0EPePLJJ7PVVluNaB/CIQAAAOgy22+/fZYtW5bHHnsszz33nFvMGKDWmueeey6PPfZYli1blu23335E+9u8Q/UCAAAAOmTq1KmZNWtWnnjiidx7771ZvXr1eFdpwnl42cr8bMW0jq031iZNmpStttoqs2bNytSpU0e0L+EQAAAAdKGpU6dmxowZmTFjxnhXZUJacNoVuXfRkR1br5e5rQwAAACgwYRDAAAAAA0mHAIAAABoMOEQAAAAQIMJhwAAAAAaTDgEAAAA0GDCIQAAAIAGEw4BAAAANJhwCAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoMGEQwAAAAANJhwCAAAAaDDhEAAAAECDCYcAAAAAGmzYcKiUclEp5ZFSyl1DLH93KeWOvn8/KKX8Yb9l95ZS7iyl3FZKubmTFQcAAABg5NrpOXRxkiM2sPxXSQ6pte6f5G+SXLDe8nm11gNqrXM2rYoAAAAAjJbNh1uh1vq9UsrsDSz/Qb+HP0qy68irBQAAAMBY6PSYQ+9L8m/9Htck15RSlpRSTtjQhqWUE0opN5dSbn700Uc7XC0AAAAAWhm251C7SinzsiYcem2/4oNrrQ+UUnZMcm0p5ee11u+12r7WekH6bkmbM2dO7VS9AAAAABhaR3oOlVL2T3JhkqNrrY+vLa+1PtD3/yNJLk8ytxO/DwAAAIDOGHE4VEqZleTrSd5ba727X/kWpZSt1v6c5E1JWs54BgAAAMD4GPa2slLKpUnekGR6KWVpkrOSTE6SWusXknw0yYuTnF9KSZLn+2Ym2ynJ5X1lmyf5Sq31qlF4DgAAAABsonZmK3vnMMv/IslftCj/ZZI/3PSqAQAAADDaOj1bGQAAAAA9RDgEAAAA0GDCIQAAAIAGEw4BAAAANJhwCAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAbbfLwrQDJz22mZfdoVLcu/f9obx6FGAAAAQFMIh7rAUAFQq8AIAAAAoJPcVgYAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoMGEQwAAAAANJhwCAAAAaDDhEAAAAECDCYcAAAAAGkw4BAAAANBgwiEAAACABhMOAQAAADSYcAgAAACgwYRDAAAAAA0mHAIAAABosM3HuwIAAAAAG3LwosW5f/nKAWUzt52W75/2xnGq0cQiHAIAAAC62v3LV+beRUcOKJt92hXjVJuJRzg0AUhQAQAAgE0lHJoAJKgAAADAphIOAQAAAD1n5rbTWnaMcCfNxhMOAQAAAD1nqADInTQbz1T2AAAAAA2m59AE1ap7na51AAAATHRD3W7Waj3WEA5NUK1CIF3rAAAAmOh0ith4bisDAAAAaDDhEAAAAECDua2s4Q5etDj3L185oMzYRAAAANAcwqGGu3/5yty76MgBZcYmAgAAgOZwWxkAAABAgwmHAAAAABpMOAQAAADQYMOGQ6WUi0opj5RS7hpieSmlfLaUck8p5Y5Syh/1W3ZEKeU/+5ad1smKAwAAADBy7QxIfXGSf0zy5SGWL0iyZ9+/g5J8PslBpZRJST6X5LAkS5PcVEr5Vq31pyOtNAAAADAxDTWrNqNn2HCo1vq9UsrsDaxydJIv11prkh+VUrYtpcxIMjvJPbXWXyZJKeWyvnWFQwAAAEBLrWbVZnR1YsyhmUnu6/d4aV/ZUOUtlVJOKKXcXEq5+dFHH+1AtQAAAAAYTifCodKirG6gvKVa6wW11jm11jk77LBDB6oFAAAAwHDaGXNoOEuT7Nbv8a5JHkgyZYhyAAAAALpEJ3oOfSvJf+ubtezVSX5ba30wyU1J9iyl7FFKmZJkYd+6AAAAAHSJYXsOlVIuTfKGJNNLKUuTnJVkcpLUWr+Q5Mokb05yT5Knkxzft+z5UsoHk1ydZFKSi2qtPxmF5zBhzdx2WmafdsWgsu+f9saO7g8AAABornZmK3vnMMtrkpOGWHZl1oRHbIJWIdD64c5I9wcAAAA0WyduKwMAAACgRwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoMGGncoeAAAAYGMcvGhx7l++ctj1Zm47Ld8/7Y1jUCM2RDgEAAAAdNT9y1fm3kVHDrve7NOuGIPaMBy3lQEAAAA0mHAIAAAAoMGEQwAAAAANJhwCAAAAaDADUveYmdtOGzRg18xtp41TbQAAAGDT+Ru3OwiHeowp/gAAAJgo/I3bHdxWBgAAANBgwiEAAACABhMOAQAAADSYcAgAABRJm7sAABsESURBVACgwYRDAAAAAA0mHAIAAABoMOEQAAAAQIMJhwAAAAAaTDgEAAAA0GDCIQAAAIAGEw4BAAAANJhwCAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA02ObjXQEAAACgNxy8aHHuX75yQNnMbafl+6e9cZxqRCcIhwAAAIC23L98Ze5ddOSAstmnXTFOtaFT3FYGAAAA0GDCIQAAAIAGc1sZg8zcdtqgboHuIQUAAICJSTjEIK1CIPeQAgAAwMQkHAIAAAAGGWpmMiYe4RAAAAAwSKuZyZiYhENMaEMl3cZPAgAAgDXaCodKKUck+UySSUkurLUuWm/5/0zy7n77/C9Jdqi1PlFKuTfJiiSrkzxfa53TobrDsFol3cZPAgAAgN8bNhwqpUxK8rkkhyVZmuSmUsq3aq0/XbtOrfWTST7Zt/5RSf6q1vpEv93Mq7U+1tGaM6bMYAYAAAATUzs9h+YmuafW+sskKaVcluToJD8dYv13Jrm0M9WjW5jBDAAAACamzdpYZ2aS+/o9XtpXNkgp5UVJjkjyr/2Ka5JrSilLSiknDPVLSiknlFJuLqXc/Oijj7ZRLQAAAABGqp1wqLQoq0Ose1SS7693S9nBtdY/SrIgyUmllNe32rDWekGtdU6tdc4OO+zQRrUAAAAAGKl2bitbmmS3fo93TfLAEOsuzHq3lNVaH+j7/5FSyuVZc5va9za+qgAAAEC3GWqMWnpHO+HQTUn2LKXskeT+rAmA3rX+SqWUbZIckuQ9/cq2SLJZrXVF389vSnJOJyoOAAAAjD8TFfW+YcOhWuvzpZQPJrk6a6ayv6jW+pNSyol9y7/Qt+oxSa6ptf6u3+Y7Jbm8lLL2d32l1npVJ58A48cMZgAAAL3n4EWLc//ylQPK/C3XbO30HEqt9cokV65X9oX1Hl+c5OL1yn6Z5A9HVEO6lhnMAAAAes/9y1fm3kVHDijzt1yztTMgNQAAAAATVFs9h2Ci060SAACAphIOQXSrBAAAoLmEQ3SUQaoBAABGn7sf6CThEB1lkGoAAIDR5+4HOkk4xKjrtt5EQ9UHAACg27TqIZT4G4bOEg4x6rqtN5FulgAAQK9o1UMIOk04BAAAAKNorMYH6ra7NugdwiEAAAAYRWM1PlC33bVB7xAOjbbvfGJw2bzTx74eXaZVoj3UelJuAAAAGD3CIcZFu4GPlBsAAABGl3AIRsFY3VMMAACbzF0OXWeovyPGglmdm004BKNgrO4pBgAAJo7xnJnMhexmEw7Rk/TMAQAAgM4QDtGT9MwBAACAzthsvCsAAAAAwPgRDgEAAAA0mNvKYAhDjdZvXCMAYFSYOQoGMdYojA3hEAyh1QnHuEYAADB2jDUKY0M4BADjrdO9BfQ+YDx43wF0paHuiOhqziljTjgEAAAAE5Rb8GiHcAjGiDGMABrOVVDYdI4fgFElHIIxYgwjAAAAupFwCDaC3j8AQCO06qnTynj23tGbCKBjhEOwEfT+AQAAYKIRDsFE4yoa3ard9+ZovIcdF3SK9xLQC3rxs2oc63zwosW5f/nKAWXuDmiwlj0n9x/zaow14RAAAACNdf/ylbl30ZEDytwdQNMIh6ABXA0BAOiQXuyVM568XtAThEPQAK6GAAAAMBThUDeTsveEUzb/l+Q7d6xXOvHvSZ0o9KoCGEW+y4wOrysAHSYcAhpNryoAAKDphENsOletJpxWvWjO3OKbed9r9xi4onbubY5d+hvPWeTwusJE4VhurJH0Qj9zi2/m7/+vf1lv2/d2tH7QLuEQsE6rXjTrn7AAAIA1RtILfdAF2CSnzDO0AeNDOMT46fAVll4c+6d1nTMmV5rG5PUax6toX/ybE7Li2VUDyr625XuNJQTQyzp9XtHbY+LRpnTIzG2nDQp5jEu5kXrheGxVx4YSDgET0opnV+WU+S8fUPb3V68cYm0AAPi9ViGQcSmZyDYb7woAAAAAMH70HKKrterOubZ8fVtNnZy/v+7u9dY7qLMVGkm3Q10Wobv0QldnGG+9cO5yLAP0tl441zSAcIiutjH39BrQDQAAADaecAgSVx2BNXrxs6CJV9t6sZ2YeIY69rwXYcIayYQuQ05EA13CmEMAAAAADabnUK9p92qpq6rjqtNXBg5etDj3Lx8405apNDtjJFeAALqO8z90j9E4HsfiGB/J7xjP3qwtfvfBPz5o0HfoVlqNZzqhODfQBuEQ9ID7l6/MvYuOHFBmKk0AABhaq+/QQGvCIdrTxDEtus16bXDK5ncnae9k1/U9mVzN2DhjNc6FdqGXTPT3a6ef33i+XuP1u32XYQL44o2/yopnVw0o+9qPF7f1HazV97czt/hVy0ld2tHp/Y2riX4OGQvj9Rmr7TqmrXColHJEks8kmZTkwlrrovWWvyHJN5P8qq/o67XWc9rZFug9ejIBAIy9Fc+uyinzXz6g7O+vHv62qaT197e//7/+ZZPr0un9AeNr2HColDIpyeeSHJZkaZKbSinfqrX+dL1Vb6i1vmUTt6WbuLJGD2l11SpJztxi8jjUhgmp02MvuJq1RrePm0Fzed/AqJtQvY7oHf7O3aB2eg7NTXJPrfWXSVJKuSzJ0UnaCXhGsi3AsIa8l9xUoQAAXUmvI+g+7YRDM5Pc1+/x0iQHtVjvNaWU25M8kOQjtdafbMS2KaWckOSEJJk1a1Yb1QJaaXUv+lZTR78Xzcxtpw26tWxCXQFyJZlOaeJVq25/zmM1jle367JZhhr3+g+l248fNspQPZ57cRbaTo9p2avMfMtE0U44VFqU1fUe35Jk91rrU6WUNyf5RpI929x2TWGtFyS5IEnmzJnTch1geK3uRR8Lrb7QuAIEAPB7Q/V4NnYjMN7aCYeWJtmt3+Nds6Z30Dq11if7/XxlKeX8Usr0drZtpKZdAeq259tt9WmaFq//SGbeGA0TqhfURHm/61VAr5koxx5MBO2eQ5xrxsxIex2NyXe1Fu+HraZOzt9fd/egso5q4myS7XJuHVXthEM3JdmzlLJHkvuTLEzyrv4rlFJ2TvJwrbWWUuYm2SzJ40mWD7ctwEhm3hgNekEBAHSv8fqu1pMXCqFNw4ZDtdbnSykfTHJ11kxHf1Gt9SellBP7ln8hybFJPlBKeT7JyiQLa601ScttR+m5wITwxRt/lb+5euCVkJnbThu03lZTJ7fsgtzts3SN15hIXafbr8yMVLtXdlwBmnhG0qbeD9BZHT7XtDyH33jC4D+YJ9L5DDbBmVt8c1BY1cjvu+PF94lN0k7PodRar0xy5XplX+j38z8m+cd2twWGtuLZVUPMvvXjAQ/f99o98r55vTdL13iNiQQAjEzLnr7r3WID6GFEb2orHIKuIw3uPuN1X/Z46rbeP+N1XHTT73W1monAOQ4YBU2cXazlc/ZdYQ3nGtYjHAJGjasmAAAA3U84NBFM9NS3ieNXjFO9m3hFaUz0wvvQGEEj1+ZMfFtNndzR4HRUfod2nli0J8PxHhnaCF6bVp/PM7c9qO3tD160OPcvHzhBR6uxJX1/a6ih3pt6RrGJhEMAMErGYnwOY4AAdKdWn8+nzBs8y9ZQ7l++cvA4lEIgYJQIh4CuNHPbaYNmY5u57bSWU5cyOsai1wuMplbv4a/9eHHnP0ca1uui1euadNnnw3i2ifHQGmk0ZmNtt0fQUGM8dvJ4bLcXE9C7hENAV2r1x9v6YRGjS48Uel3L9/DVK4dYm3YNNeukzweabDxnY20VAnX6eNSLCSY+4RD0sol8tbrFcztl87uTHDl43W7S6TGyXG2G8dPtn7HdXj9+z+d7y54nyRC9gifye3uI53bK5ncLW+hOPr8aQzgEAACMqpY9T6JXMEC3EA7BCA11n3e369V6tzSRrzBCp4zGceLYg003XsfPWPzeXu39y4iNZPyjVtv+/dVXZOa200b0uwfpxXNXL9aZniMcghHqmsE3N1Kv1hsAgO40kvGPWm17yrz2A0XfbWFkhEMwFAl9d1qvXdZciRyfASDHhPchwPjrhc/iEdSxF2an7HgdR+H1orUJ1Vu9F/TC5xVdSTgEAAAN1guzU3ZTHcdzZrJe1E0hIzA04RCdJaluhFM2/xczarDReuHKdNt81vUsn189xHHWWHrmjK9Wr//MbQ8ap9rQlTo9Oy9dQTgEwJjopqu+AHQvPXPGV6vX/5R5bxyn2gBjRTjUJK1S2nmnj309gAllLHoETaheR+0apytrrV7rpAGvN0ALejEBTSEcAmBExqJHkF5HY2eoK/Zeb6CJ9GICmkI41HTu+aRBZm47LbNPu2JA2ZlbtLj610U9NvTWGFrL2U9uPKEnXy9tTzcYz/fhSH6346e3jcU4YL3wHjGjFzDehENAY3z/tBb3y3fRwLR6x2ycVl/qe/X10vZ0g/F8H47kdzt+GE4vvEe6KagCmmmz8a4AAAAAAONHzyEYR7oQb5wRv14T+DbKXugyD0DvGJXvKBP4PMwo8Z6BMSMcgnHkD/eN4/UaWi90mQegdzjnAjSLcAigg/TgYTSZUnl0tHpdh+J47l1N/Hxu4nMeiYn0Gat3OrCxhEMAHaQHD6PJlMqjY2NeV8dz72ri53MTn/NITKTPWAEgsLGEQwCjbKird+18cRvJtkNpt5dEE68wjsVVdldzh6aXAxNVN723u+0zaCT1GY1zJEBTCYcARtlIplwfjenaJ9KV0U4bi6vs/mgZml4OTFTd9N7uts+gkdRnNM6RAE0lHAI2qNuuME4Uo3GltBf1wlXfXqgjvavTPUq6qYfKUCbSuC50n4l0juwqZg2DCU84BGxQN/1BMZF0+kppr+qFq769UEd6V6d7lHRTD5Wh6L3IaJpI50iAsSQcAqBxhhp3ydVlRpMeMyPnNYTO0tMKWEs4BEDj6LnAePC+GzmvIXSWnlbAWsIhgDYYd4Ze0ur9OtR63sPQe3phbCkAeotwCKANxp2hl7T7B6L3MPSmXhhbCoDeIhyiu5gJARpvNGZyG4ur6b04FoqxJmDicDwDMBLCIQC6Sqdnchurq+m9OBaKW1Bg4nA8AzASwiEAOs4V7Imn28fd6sWeWxujiWPMjEYvwk4aqzbp9vd2t382ANAe4RAAHeePgomn28fd6sWeWxujiWPMdLoXYaeNVZt0+3u72z8bAGiPcAiAtukRRH96DAxtJL09On2cbczsdZ3W7b1eOm2o17rd59y01wuA7iEcAqBt/uinPz0GhjaS3h6dPs7G87jt9l4vnTbS17pprxcA3UM4BLCJ9KLpDaPRTtp+7OidxGjamF5V3nMATGTCIYBN5A+F3jAa7aTtx47eSYymdo9l7zkAJjrhEMAENdKxLwA2ppecHnUA0LuEQwATlN4twEhtzOeIzxwA6F3CIQCgY9odI8isTCOnpw5sOscPwEBthUOllCOSfCbJpCQX1loXrbf83UlO7Xv4VJIP1Fpv71t2b5IVSVYneb7WOqczVQcAuk27YwSZlWnk9NSBTef4ARho2HColDIpyeeSHJZkaZKbSinfqrX+tN9qv0pySK11WSllQZILkhzUb/m8WutjHaw3ANBQI7nir7cAAMBg7fQcmpvknlrrL5OklHJZkqOTrAuHaq0/6Lf+j5Ls2slKAgCsNZIr/noLAAAM1k44NDPJff0eL83AXkHre1+Sf+v3uCa5ppRSk/xTrfWCVhuVUk5IckKSzJo1q41qAQC9YCL11ml3TCVosol0zAM0RTvhUGlRVluuWMq8rAmHXtuv+OBa6wOllB2TXFtK+Xmt9XuDdrgmNLogSebMmdNy/wBA75lIwUm7YypBk02kYx6gKTZrY52lSXbr93jXJA+sv1IpZf8kFyY5utb6+NryWusDff8/kuTyrLlNDQAAAIAu0E44dFOSPUspe5RSpiRZmORb/VcopcxK8vUk76213t2vfItSylZrf07ypiR3daryAAAAAIzMsLeV1VqfL6V8MMnVWTOV/UW11p+UUk7sW/6FJB9N8uIk55dSkt9PWb9Tksv7yjZP8pVa61Wj8kwAAMaB8VXYFN43AHSTdsYcSq31yiRXrlf2hX4//0WSv2ix3S+T/OEI6wgA0LWMr8Km8L4BoJu0c1sZAAAAABOUcAgAAACgwYRDAAAAAA3W1phDAADQVAaPBmCiEw4BAMAGGDwagInObWUAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABrMgNQA0IXMjgQAwFgRDgFAFzI7EgAAY8VtZQAAAAANJhwCAAAAaDDhEAAAAECDCYcAAAAAGkw4BAAAANBgwiEAAACABhMOAQAAADSYcAgAAACgwYRDAAAAAA0mHAIAAABoMOEQAAAAQIMJhwAAAAAaTDgEAAAA0GDCIQAAAIAGEw4BAAAANJhwCAAAAKDBhEMAAAAADSYcAgAAAGgw4RAAAABAgwmHAAAAABpMOAQAAADQYMIhAAAAgAYTDgEAAAA0mHAIAAAAoMGEQwAAAAANJhwCAAAAaDDhEAAAAECDCYcAAAAAGkw4BAAAANBgwiEAAACABhMOAQAAADSYcAgAAACgwdoKh0opR5RS/rOUck8p5bQWy0sp5bN9y+8opfxRu9sCAAAAMH6GDYdKKZOSfC7JgiSvSPLOUsor1lttQZI9+/6dkOTzG7EtAAAAAOOknZ5Dc5PcU2v9Za31uSSXJTl6vXWOTvLlusaPkmxbSpnR5rYAAAAAjJPN21hnZpL7+j1emuSgNtaZ2ea2SZJSyglZ0+soSZ4qpfxnG3XrdtPLeXlsvCvBuJieaPuG0vbNpe2bS9s3k3ZvLm3fXNq+uSbS3/a7typsJxwqLcpqm+u0s+2awlovSHJBG/XpGaWUm2utc8a7How9bd9c2r65tH1zaftm0u7Npe2bS9s3VxPavp1waGmS3fo93jXJA22uM6WNbQEAAAAYJ+2MOXRTkj1LKXuUUqYkWZjkW+ut860k/61v1rJXJ/ltrfXBNrcFAAAAYJwM23Oo1vp8KeWDSa5OMinJRbXWn5RSTuxb/oUkVyZ5c5J7kjyd5PgNbTsqz6Q7Tajb5Ngo2r65tH1zafvm0vbNpN2bS9s3l7Zvrgnf9qXWlkMAAQAAANAA7dxWBgAAAMAEJRwCAAAAaDDh0AiVUt5eSvlJKeWFUsqQU9uVUo4opfxnKeWeUspp/cq3L6VcW0r5//r+325sas5ItdN2pZS9Sim39fv3ZCnllL5lf11Kub/fsjeP/bNgU7R73JZS7i2l3NnXvjdv7PZ0lzaP+d1KKd8ppfys79xwcr9ljvkeM9S5u9/yUkr5bN/yO0opf9TutnS3Ntr+3X1tfkcp5QellD/st6zlZz+9oY22f0Mp5bf9Pss/2u62dK822v1/9mvzu0opq0sp2/ctc8z3sFLKRaWUR0opdw2xvDHneuHQyN2V5I+TfG+oFUopk5J8LsmCJK9I8s5Syiv6Fp+W5Lpa655Jrut7TG8Ytu1qrf9Zaz2g1npAkldmzYDtl/db5e/WLq+1XjkmtaYTNua4ndfXvv3DY8d9b2qn3Z5P8j9qrf8lyauTnNTv8z5xzPeMYc7day1IsmffvxOSfH4jtqVLtdl+v0pySK11/yR/k8EDlbb67KfLbcSxe0O/z/JzNnJbukw7bVdr/WS/7/SnJ/lurfWJfqs45nvXxUmO2MDyxpzrhUMjVGv9Wa31P4dZbW6Se2qtv6y1PpfksiRH9y07OsklfT9fkuRto1NTRsHGtt38JL+otf56VGvFWBjpceu4703Dtlut9cFa6y19P69I8rMkM8eshnTShs7dax2d5Mt1jR8l2baUMqPNbelew7ZfrfUHtdZlfQ9/lGTXMa4jo2Mkx67jvndtbNu9M8mlY1IzRl2t9XtJntjAKo051wuHxsbMJPf1e7w0v/9jYada64PJmj8qkuw4xnVj021s2y3M4BPJB/u6J17k1qKe0m7b1yTXlFKWlFJO2ITt6S4b1W6llNlJDkzy437FjvnesaFz93DrtLMt3Wtj2+99Sf6t3+OhPvvpfu22/WtKKbeXUv6tlLLPRm5L92m77UopL8qaXib/2q/YMT+xNeZcv/l4V6AXlFL+PcnOLRb9n7XWb7azixZldWS1YixsqO03cj9Tkrw1a7qhrvX5rOmKXvv+/7+T/Pmm1ZRO61DbH1xrfaCUsmOSa0spP++7OkGX6uAxv2XWfHE8pdb6ZF+xY763tHPuHmod5/3e1nb7lVLmZU049Np+xT77e1c7bX9Lkt1rrU/1jR33jay53cRx37s2pu2OSvL99W4pc8xPbI051wuH2lBrPXSEu1iaZLd+j3dN8kDfzw+XUmbUWh/s6572yAh/Fx20obYvpWxM2y1Ickut9eF++173cynl/0ny7U7Umc7oRNvXWh/o+/+RUsrlWdP99Htx3HetTrR7KWVy1gRD/2+t9ev99u2Y7y0bOncPt86UNrale7XT9iml7J/kwiQLaq2Pry3fwGc/3W/Ytu8X+KfWemUp5fxSyvR2tqVrbUzbDboTwDE/4TXmXO+2srFxU5I9Syl79PUgWZjkW33LvpXkuL6fj0vSTk8kusPGtN2ge5P7/rhc65isGdyc3jBs25dStiilbLX25yRvyu/b2HHfm9pp95Lki0l+Vmv99HrLHPO9ZUPn7rW+leS/9c1k8uokv+275bCdbelew7ZfKWVWkq8neW+t9e5+5Rv67Kf7tdP2O/d91qeUMjdr/p56vJ1t6VpttV0pZZskh6Tf+d8x3wiNOdfrOTRCpZRjkvxDkh2SXFFKua3WengpZZckF9Za31xrfb6U8sEkVyeZlOSiWutP+naxKMlXSynvS/KbJG8fh6fBpmnZdv3bvu/xi5IcluT/WG/7vy2lHJA13Q/vbbGc7tVO2++U5PK+74+bJ/lKrfWqDW1P12un3Q9O8t4kd5ZSbuvb7oy+mckc8z1kqHN3KeXEvuVfSHJlkjcnuSdrZqM8fkPbjsPTYBO02fYfTfLiJOf3fc4/3zdL0YY+++lybbb9sUk+UEp5PsnKJAtr/f/bt2MbBGIgiKKzAR1QEJ1RDj1cPbRAvASXEJKe5r3ckbW29C3vJjH3F/Xnvifnw86xu5+f5Wb+4mbmleSR5D4z7yTPJLek766f8ywDAAAAoJFvZQAAAADFxCEAAACAYuIQAAAAQDFxCAAAAKCYOAQAAABQTBwCAAAAKCYOAQAAABT7AiGNeamCj7/aAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 1440x504 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"\n",
"plt.hist(cos_theta_l_np_Bbar,bins=200, histtype=\"step\", label=r'$cos\\theta_{l}^{(2)}$', density=True);\n",
"plt.hist(cos_theta_l_np_q,bins=200, alpha=0.5,label=r'$cos\\theta_{l}^{(3)}$',density=True);\n",
"plt.legend(fontsize='20');\n",
"fig = plt.gcf()\n",
"fig.set_size_inches(20,7)\n"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [],
"source": [
"if TASK=='SAMPLE':\n",
" \n",
" with open(filename,'wb') as handle:\n",
" pickle.dump(MC_tuple, handle)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"#bin_n=12\n",
"#sample_for_hist_reduced=sample_for_hist[0:20000,:]\n",
"#H, edges = np.histogramdd(sample_for_hist_reduced/n_events, bins=bin_n)\n",
"#\n",
"#x=np.array([H[i,:].sum() for i in range(bin_n)])\n",
"#plt.plot(x)\n",
"#\n",
"#y=np.array([H[:,i,:].sum()/n_events for i in range(bin_n)])\n",
"#plt.plot(y)\n",
"#\n",
"#z=np.array([H[:,:,i].sum()/n_events for i in range(bin_n)])\n",
"#plt.plot(z)\n",
"#\n",
"#plt.plot(np.array([H[:,:,:,i].sum() for i in range(bin_n)]))\n",
"#plt.plot(np.array([H[:,:,:,:,i].sum() for i in range(bin_n)]))\n",
"##plt.plot(np.array([H[:,:,:,:,:,i].sum() for i in range(bin_n)]))\n",
"#\n",
"#"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.7"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Davide Lancierini