{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dd0f9938",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from scipy.optimize import minimize\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "090aaa0d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAG4CAYAAABb+t1HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAAxOAAAMTgF/d4wjAABI+klEQVR4nO3de3yT5f3/8Xfa0lag1NaCWMKxZYXKQTk4N0FAQEGZCJQpiMgYiEKhICKK4ECpoFNAJ5WhP8XhrG7KJggylcHcdDIEGWo7y0nKQbAF1kLPh/v3B18y0qY0CUnuJH09H48+sPeV3LluYsm71+e6rttiGIYhAACAABJidgcAAABcRYABAAABhwADAAACDgEGAAAEHAIMAAAIOGFmd8CbqqurdezYMUVFRclisZjdHQAA4ATDMHTmzBnFx8crJMTxWEtQB5hjx46pdevWZncDAAC44fDhw7JarQ7bgjrAREVFSTr3F9CsWTOTewMAAJxRWFio1q1b2z7HHQnqAHO+bNSsWTMCDAAAAeZi0z+YxAsAAAIOAQYAAAQcAgwAAAg4BBgAABBwCDAAACDgEGAAAEDAIcAAAICAQ4ABAAABhwADAAACDgEGAAAEHAIMAAAIOAQYAAAQcAgwAAAg4BBgAABAwCHAAACAgEOACWJllVVa/lGOyiqrzO4KAAAeRYAJYuWV1Xp+y16VV1ab3RUAADyKAAMAAAIOAQYAAAQcAgwAAAg4BBgAABBwCDAAACDgEGAAAEDAIcAAAICAY2qA+eSTT2SxWGxfv/71r21tTzzxhF3bjh07bG179+7V2LFjlZqaqtTUVJWVlZnRfQAAYJIwM198zZo1Wr58ue37e++9V5JUWlqqXbt22dqaNGmi3r17S5KKiop066236oMPPlBiYqJefPFFpaWladWqVb6/AAAAYArTAkx2drasVqtmzpxZq23NmjVKTU3VoEGDarWtWrVKVqtViYmJkqSxY8eqefPmeuSRR9SuXTsv9xoAAPgD00pIzz77rJ566indeOON2rRpk+14dXW1li1bpjvuuEOjRo3St99+a/e89evXq0uXLrbvY2Nj1apVK23cuNFnfQcAAOYyLcCkpKTo2WefVVlZmW677TYtWbJEklRSUqL09HTNmTNHO3bsUM+ePbVlyxbb87KyshQXF2d3rpiYGOXk5NT5Wp07d5bVapXVatWyZcu8c0EAAMBnTCshDR06VEOHDlVaWpoWLFigxx9/XCNGjFCnTp00evRoSdKcOXM0fPhwTZw4UQcOHFBoaKgKCgoUGxtrd66IiAgVFxfX+VrZ2dlq1qyZV68HAAD4junLqC0Wi5588kl16dJFW7dutWtr3Lix1q5dqx9++ME2whIbG1tr1VFJSYliYmJ81mcAAGAu0wOMdC7E3HTTTQ6XQ7ds2VLJycm2toSEBOXn59s9Ji8vT0lJST7pKwAAMJ9fBBhJ+v777zVkyJBax6uqqhQeHq6uXbtKOjd35sI9YU6ePKmTJ0/qlltu8VlfAQCAuUwJMLm5uZo+fbpyc3MlSW+//bauvvpqderUSR9//LEWLlyooqIiVVZW6tFHH9XixYsVGhoqSZowYYL27dunEydOSJLeeOMNTZkyRVar1YxLAQAAJjAlwISFhelvf/ubkpOT1a9fP1VWVuqxxx6zta9evVrt2rXTiBEjNHr0aA0cONDWFhMTo3fffVczZ87U7NmzdfjwYbvN8PA/2ccK7f4EACBYWAzDMMzuhLcUFhYqOjpaBQUFDW4V0jdHCzRhzb+Ud6ZczaPCtWbCdbq6VbTZ3QIAoF7OfH77zRwYeNZHWSeUd6ZckpR3plwfZ/9gco8AAPAcAkyQGpx8pZpHhUuSmkeFa1DnFib3CAAAzyHABKmrW0Vr5ZgekqSVY3pQPgIABBUCTBDrHN/M7k8AAIIFAQYAAAQcAgwAAAg4BBgAABBwCDAAACDgEGAAAEDAIcAAAICAQ4ABAAABhwADAAACDgEGAAAEHAIMAAAIOAQYAAAQcAgwAAAg4BBgAABAwCHAAACAgEOAAQAAAYcAAwAAAg4BBgAABBwCDAAACDgEmCAWHhaitIEdFR7G2wwACC5hZncA3hMRFqpZg39kdjcAAPA4fjUHAAABhwADAAACDgEGAAAEHAIMAAAIOAQYAAAQcAgwAAAg4BBgAABAwCHAAACAgEOAAQAAAYcAAwAAAo6pAeaTTz6RxWKxff3617+WJBUXF2vy5MmaMWOGxo0bp0OHDtk9b/v27RozZoymTp2qefPmyTAMM7rvV8oqq7T8oxyVVVaZ3RUAALzOYpj46T9x4kR169bN9v29996rmJgYpaSk6Pbbb9f48eP11Vdf6e6779bOnTvVqFEjHT16VH379tXOnTsVExOjhx56SHFxcXrkkUdqnb+wsFDR0dEqKChQs2bNfHlpPnemtEJdF36orxberKjIRmZ3BwAAtznz+W3azRyzs7NltVo1c+ZMu+M7d+7Uxo0b9fvf/16S1LVrV1ksFmVmZmr8+PFaunSp+vXrp5iYGEnS3XffrRtvvFHTp09XkyZNfH0ZAADABKaVkJ599lk99dRTuvHGG7Vp0ybb8fXr1ysxMVERERG2Y927d9eGDRts7V26dLG1de3aVcXFxdq2bZvP+g4AAMxlWoBJSUnRs88+q7KyMt12221asmSJJCkrK0txcXF2j42JiVFOTo6KioqUm5tr1x4WFqamTZsqJyenztfq3LmzrFarrFarli1b5p0LAgAAPmNaCWno0KEaOnSo0tLStGDBAj3++OMaMWKECgoKFBsba/fYiIgIFRcXq6CgQJLqbK9LdnZ20M+BAQCgITF9GbXFYtGTTz6pLl26aOvWrYqNjVVZWZndY0pKShQTE2MLLnW1AwCAhsH0ACOdCzE33XSTysrKlJCQoPz8fLv2vLw8JSUlKTIyUvHx8XbtJSUlOnv2rJKSknzdbQAAYBK/CDCS9P3332vIkCFKSUnRnj177EZZsrKydNttt0k6N3dmx44dtrbs7GxFR0frhhtu8HmfAQCAOUwJMLm5uZo+fbpyc3MlSW+//bauvvpqderUSddee6369OmjzZs3S5J2796tkJAQjR49WpI0ffp0bd26VaWlpZKktWvXav78+YqMjDTjUgAAgAlMmcQbFhamv/3tb3rttdfUs2dP3XfffXrsscds7W+++aZmzZqlTz/9VPn5+dq0aZNCQ0MlSYmJicrIyNCkSZMUGxurFi1a6KGHHjLjMgAAgElM3YnX27y5E29ZZZUytu7X1AEJiggL9ei53cFOvACAYOHM57ffzIEJNOWV1Xp+y16VV1ab3RUAABocAgwAAAg4BBgAABBwCDAAACDgEGAAAEDAIcCYrKyySss/ylFZZZVXnwMAQDAhwJjMndVMrIACADR0BJggkX2s0O5PAACCGQHGi3xV6vnmaIGmZe6SJE3L3KVvjhZ49fUAADAbAcZDHIUVX5V6Pso6obwz5ZKkvDPl+jj7B6++HgAAZiPAeIiZ81IGJ1+p5lHhkqTmUeEa1LmFz/sAAIAvEWCCwNWtorVyTA9J0soxPXR1q2iTewQAgHcRYHzIm3NiOsc3s/sTAIBgRoDxIXfLTOz7AgCAPQJMAGDfFwAA7BFgGhBGcgAAwYIA42e8GTIYyQEABAsCjJu8tfMtIQMAgPoRYNzAzrcAAJiLAOMGdr4FAMBcBBg3sPMtAADmIsC4gZ1vAQAwFwHGTex8CwCAeQgwJnNnNZO3VkABABAoCDAmcmc1EyugAAAgwJjKndVMrIACAIAA41X1lXrcWc3k6RVQ3F4AABCICDAeUjOsOFPqcWc1k6dXQLHzLwAgEBFgPMBRWHG21OPOaiZWQAEAGjoCjAc4Ciu+3uwuPCxEaQM7KjzMP95SSlMAAG/yj0+7AOcorDgq9bi7/NmZ50WEhWrW4B8pIizUnUvwOEpTAABvIsB4QF3zUi4s9bi7/Jll0wAA1EaA8ZD65qW4u/yZZdMAANRGgPERd+fEeHIujad28GV+i7n4+wcAAozPOLv8uWbI8NSyaU+WopjfYi7+/gHATwLMpk2b1L59e7tjTzzxhCwWi+1rx44dtra9e/dq7NixSk1NVWpqqsrKynzdZbfUV2aqK2R4Ytk0pSgAQDAJM7sDp0+f1owZM2QYhu1YaWmpdu3apeXLl0uSmjRpot69e0uSioqKdOutt+qDDz5QYmKiXnzxRaWlpWnVqlWm9N+THIWMS92o7rzByVfq9/86pLwz5T5Z1g0AgDeZPgIzf/58jR8/3u7YmjVrlJqaqpkzZ2rmzJmaPHmyrW3VqlWyWq1KTEyUJI0dO1Yvv/yyvvvuO1922yu8uXeMp3fwBQDATKYGmD//+c/q3bu32rRpYztWXV2tZcuW6Y477tCoUaP07bff2j1n/fr16tKli+372NhYtWrVShs3bvRZv73F2yGDHXwBAMHCtACTn5+vjRs3asKECXbHS0pKlJ6erjlz5mjHjh3q2bOntmzZYmvPyspSXFyc3XNiYmKUk5NT52t17txZVqtVVqtVy5Yt80j/vbXzLSEDAID6mTYHZsGCBXryySdrHW/SpIlGjx4tSZozZ46GDx+uiRMn6sCBAwoNDVVBQYFiY2PtnhMREaHi4uI6Xys7O1vNmnk2EJzf+bYhKKusUsbW/Zo6IKHOnX6deQwAAJ5iygjM22+/rf79+6tly5YXfVzjxo21du1a/fDDD7YRltjY2FqrjkpKShQTE+O1/jZ0zizbZWkvAMCXTBmBWbVqlbZt26a77rrL7rjFYtHWrVvVv39/27GWLVsqOTnZFloSEhKUn59v97y8vDwlJSV5vd8AAMA/mDIC8/LLL+vLL7+0fS1atEhXXXWVvvzyS/Xq1cvusVVVVQoPD1fXrl0lSSkpKXZ7wpw8eVInT57ULbfc4tNrCBY1N84L1F1eA7XfAAD3mBJgEhMTdc0119i+2rRpo/DwcF1zzTX6/PPPtXDhQhUVFamyslKPPvqoFi9erNDQc/MqJkyYoH379unEiROSpDfeeENTpkyR1Wo141ICmqON8wK1FBSo/QYAuMf0fWAcWb16tdq1a6cRI0Zo9OjRGjhwoK0tJiZG7777rmbOnKnZs2fr8OHDtg3vApE7q5k8tQKK3XmDFyNSAIKd6TvxSudGVc4vpx40aJCOHTt20cf36tVLmZmZPuiZ97mzmslTK6C8uTvvhaWp6zpc4bHzwjnnR6Qm9W3PqjAAQckvR2DgG97aOM+TN45s6BhJAQDHCDBe5K3N7jzJGxvnUZpyX83AwtweAHDMfz9ZA4yjsHK+1NPQhvDruqcTown1I7AAgHMIMB7SUMOKI3WVpvhw9oyaS98BoCEiwPhQIJSUnOHMByj3dPIO5hcBwDmB/UkaYNwdpfGn4GPGByilp/9hfhEAnGP+JyLq5U/lKTM+QCk9/U9d84sAoKEhwMAlnv4AZXTFNd5a+g4AgYYA04B4ohTl6Q/QhjS64qmwxvwiACDA+B1vznfxVCnKmx+gznzIB+qoTc2wFqjXAQD+gADjZ/xpvst5zqw68tTSXmdGZIJl1CZYrgMAzECAwUU5s+qIpb0AAF8jwOCinFl1FKhLeynhAEDgIsDATs1SkDOrjvxhaa87JSxKOAAQuAgwsHFUCnJm1ZHZS3spYTGaBKDhIcDApq5SkDOrjsxc2utsCcubH/JmB4iao0ncLwlAsCPAwMYfSkHucNRvR4HCmyUjfypHMSIFoCEgwMDG7FKQuxz1258CxaVwZyQlUCdVA4ArCDANXM2N8wJ1l1ezN9fzlAsDi7sjKYE6kgYAriDANHC+3DgvUOdl+GoH3ZqB5ffbc90aSXF2JM3seTsAcCkIMPCJYJqX4a3yVM3Sj0VyeyTFmRGpYCmzAWiYCDDwipqlqUCZl2HmqETN0s/YH7cJyDlJAOALBBh4Rc3S1KXMy/Bl6cnMUQlHpR9Pzu2hZAQgmLgUYPLy8vT0009r0KBBuvLKK9WoUSOFh4erbdu2GjlypNauXavy8nJv9RUBzN0VTsFUenKGNycjUzICEEzCnHlQZWWlFixYoLfeeks9evRQz549dfvtt6tx48aSpLNnz+rQoUNau3atnnzyST311FNKSUnxascReBx9OF84unJdhytqPcdR6Wlin3be7ywAwK/VG2CKioo0bdo0jRgxQkuWLKn3hGVlZXr11Ve1fPlyzZo1yyOdhH+pOb/FXTVHV9ZMuE5trmhs95jByVfq9/86pLwz5QG/JLi+sOYsT/39A0Agq/dfwFdeeUW/+c1vNHz4cKdOGBERoQceeEBDhw7Vhx9+eMkdhP/x1NJrR6MrNee7OCo9BeJybEelMHevw5dL3wHAX9UbYO6//35FRUW5fOJOnTqpX79+bnUKDUPNib3t4xo7nO9yYekpUOfE1Axrb27PDcjrAAB/UW+AiYiIsPv+qaee0pIlS/TPf/5TkvTwww8rOjpaPXv21Ndff33R5wIXqjm6ciCvqN6l1oG6HLtmWDOkgLgOAPBXLhfRV65cqVtvvVU/+clP9P/+3//Tc889p1mzZum1117Tiy++6I0+IohdOLrizFJrT26T781SVM0VPzXD2t0/bsN2/wBwCVwOMPfcc4+6d++u4uJiPf744xo6dKgWLlyobt26KTqajbbgPmeWWnvqhpNmlKIuDGuBeuNMAPAXLgeY4uJiffPNN5oyZYr++9//6oUXXpAkHT9+XJmZmR7vIMzny1UvzuyD4om9UvyhFOXJPV8CcWIzAFwKlz+R5syZoxUrVqiwsFAbN25Uhw4dtHnzZj344IO6+uqrvdFHmMyZVS/+trTX0Qf6hcc8fcdmMwOEJ0eT2K0XQKBw+dPm2LFjevnll/Xee++pf//+kqQhQ4YoNTVV77zzjqf7hwDhT0t7HX2g1zwmyWMlHLNXRnlyNIndegEECpcDzOuvv+7weKdOnTRs2LBL7hDM5W8jKc66sN+OPtAdHXO2hFPf6IrZ5ShPjyYBQCBw6lYC2dnZevrpp5Wbm6ucnBz95z//qfWY3NxcVVe791vbpk2bNG3aNB08eFDSuXk2aWlpuuyyy3Tq1Cmlp6erbdu2tsdv375dK1asUExMjC6//HKlp6fLYrG49dqwd34kJdBc2O+6du91Z0ffQNgt+PyE4J+v/tw2mnSmtMKnfYB/KausUsbW/Zo6IMEvRkUBb3Dq1+zOnTvrt7/9rdq2baurrrpK/fr1s/vq37+/7r//fm3ZssXlDpw+fVozZsyQYRi2Y+PHj1ffvn31wgsvaO7cufrZz36miopz/yAfPXpUY8aMUUZGhjIyMlReXq6nn37a5ddF8HK0wsfdVT/u7hbsa+5MCGbib/CiFIiGwKkRGOncpnSvvfaaPv74Yw0aNMjhY86cOeNyB+bPn6/x48fr1VdflSTt3LlTGzdu1O9//3tJUteuXWWxWJSZmanx48dr6dKl6tevn2JiYiRJd999t2688UZNnz5dTZo0cfn1EZwcfaA78yFf835FNUdXau4WvGbCdbq6VbRX7yLtjvpKgY5GlljKDSCQuDzRYdCgQaqurtbx48eVm5tr+zp48KAWLlzo0rn+/Oc/q3fv3mrTpo3t2Pr165WYmGi3i2/37t21YcMGW3uXLl1sbV27dlVxcbG2bdvm6qUAdhxNxnVnt2BPcndOUn2Tqs2etwMAl8rlAJORkaHY2Fi1atVK7du3t30lJCRoxYoVTp8nPz9fGzdu1IQJE+yOZ2VlKS4uzu5YTEyMcnJyVFRUpNzcXLv2sLAwNW3aVDk5OXW+VufOnWW1WmW1WrVs2TKn+4iGpa4PdVd3C74UNQOLJ1d3eXMZOQD4mssB5pFHHtGcOXP09ddf68CBA3Zfs2bNcvo8CxYs0JNPPlnreEFBgWJjY+2ORUREqLi4WAUF55an1tVel+zsbB05ckRHjhzRgw8+6HQf4Vn+vsLJmQ91b8938dZydG8uIwcAMzg9B+a8H//4x/rlL3+pli1b1mpLS0tz6hxvv/22+vfv7/AcsbGxOnv2rN2xkpISxcTE2IJLWVmZw3b4N39f4eRoNY8j/jbfxRmORpcm9mknKbCuAwDOc/lX4bVr12rt2rUO2958802nzrFq1Srdddddslgsslgs+sUvfqFDhw7JYrHo8OHDys/Pt3t8Xl6ekpKSFBkZqfj4eLv2kpISnT17VklJSa5eClBLIIYTZ7hbMmJnXgD+yuUA06tXL82dO1ehoaG1vubNm+fUOV5++WV9+eWXtq9Fixbpqquu0pdffqmlS5dqz549dqMsWVlZuu222yRJKSkp2rFjh60tOztb0dHRuuGGG1y9FAQIfy89BQJ3S18sxwXgr1wuIU2fPl0xMTHq2LGjQkL+94FSXl6u1157zalzJCYm2n2/e/duhYeH65prrpEk9enTR5s3b9bw4cO1e/duhYSEaPTo0bbXv/nmm1VaWqrIyEitXbtW8+fPV2RkpKuXggDh76UnR2oux/YHwTq6BKBhcjnAPPDAAzIMQ9HR9r/BnT592mNlnDfffFOzZs3Sp59+qvz8fG3atEmhoecmNSYmJiojI0OTJk1SbGysWrRooYceesgjrwvfC8bRFfZYAQDvcznA7N69u9YxwzC0a9cuVVRU6OGHH3a5ExMmTLBbTh0XF1fnPBvp3M0jhwwZ4vLrwP/UHF0JhkDjaMLs1a2ia11bMFwrAJjF5QBz/g7UjsTHx7sVYIDznCkX+fsHf133Rqp5bYFYGnOE++4AMIPLnwDp6emqqqpSdXW17auyslKLFi3Srl27vNFHwI639krxFH+4N5K73Lk/EhN9AZjB5QDz4IMP1rrzc0hIiIYPH6677rrLYx0DAlkgTph1dCsF+D+WuqOhcrmEdOLEiVrHSktL9dZbb+mLL77wSKcAX/D3UpSvXWyzO/iv8yNgk/q299tRScAbXA4w7dq1qzUCYxiGGjdurCVLlnisY4C3BcscFE+pa+5OTcx58X/+uIwf8DSXf/VMTU3Vvn377O6BdPToURUUFGj69One6COAenhiNMnZuTvMefFvlALRULg8AvPYY4/pyiuvlCQdP35c4eHhtW6uCJjN0Qe6N0tGZpejPDWaFIhzd2CvrmX8QLBx+V/bK6+8Us8//7yuuuoqtWrVSs2bN1fXrl2VmZnpjf4BbnG0Usmbq5f8fWUUGg5373sFBBqXR2AWL16s559/Xr/4xS/UvXt3tWjRQqdPn9aGDRtUWFioKVOmeKOfgGnMHl0BXOHsXdWBQOdygPnoo4+0d+9eXX755XbHf/7zn+uBBx7wVL8AvxGsk30JZsGLUiAaApcDTO/evWuFl/MOHjx4qf0B/F6wfPAHazAD0DC4/C9weXm5PvzwQxmGIUmqqqrS3/72N40YMUKXXXaZxzsI+JIz4YT5LgBgPpcDzOLFi7Vs2TI1adJEzZs3V0REhG666SadOnVKv/3tb73RR8BnCCf23Lm1AAD4gsslpGbNmmnz5s364osv9Omnn0qSevXqpRtuuMHjnQPgOa6WvmruJ7JmwnVMCAXgN1wOMJK0b98+denSRb169VJRUZF2794twzBq7dALwH+4OueF/UQA+DOXS0irVq1SUlKSVq9eLUlq0qSJ4uLiNHr0aB05csTjHQRgDvYTAeDPXA4wb775pv74xz/qvvvusx1LSkrSzTffrEmTJnm0cwDM4+ytBZgnA8AMLgeYHj16aOTIkYqMjLQ7XlhYaJsTAyA41LefCPfdAWAWlwNMSEiI1q1bZ/u+qKhIGRkZWrRokfr27evRzgGBKlj2iqmPo3ky8C1GwNBQufyva3p6ujZs2KCmTZsqPj5e0dHRSk1NVd++ffXaa695o49AwGkoy7GZJ2MuRsDQkLkcYC677DK99tpr+uabb7R69WqtW7dO//nPf7Rp0ybbXaoBBK8Lf+N3dp4MvIMRMDRk9QaY9957z+Hxtm3batiwYbr99tv1ox/VXppZVlambdu2XXIHAfgPR7/xc98d8zAChoas3gATHx+vZ555RhUVFU6f9Ntvv9WvfvUrNrcDAkx9c3e8/Rt/WWWVln+Uo7LKKo+eN1gxAoaGrN6N7Hr37i2LxaIBAwaod+/euvHGG5WUlKSWLVuqcePGqqioUGFhoQ4fPqyvvvpKH330kaKjo7VixQo1atTIF9cAwEPq2+xucPKV+v2/DinvTPlFf+Mvq6xSxtb9mjogwaV5QOWV1Xp+y15N6ts+6OcPeQojYGionJoD06tXL23ZskXt2rXT/Pnz1aVLFzVv3lxNmjTR5ZdfrjZt2mjAgAF6//33lZqaqpdffllNmjTxdt8B+Jizv/GfDyLlldW+7B6ABsTpWwlEREQoLS1NaWlpOnr0qL7++mudPHlSl112mVq2bKkePXooIiLCm30F4Af4jd//NZRl/GjY3LoXUqtWrdSqVStP9wUA4AGu3vcKCETEcwA+4akJukz0BSARYAD4iKfmxTC/BoBEgAEAAAHI5QDz3//+t862FStWXEJXADQklIIAXAqXA8zUqVNrHSsvL9eCBQs0d+5cj3QKQPCjFATgUrgcYPbu3au77rpL+/fvlyS9//77Sk5O1urVq5WQkODxDgIwD8txAfgrl5dRn7+/0X333acjR45o+/btmjFjhhYsWHDR8hKAwMNyXAD+yuVfq5566ik9+OCD+sMf/qBGjRpp5MiRGjt2rKKiotS6dWunz5OXl6dhw4YpKipKvXr10r///W+79ieeeEIWi8X2tWPHDlvb3r17NXbsWKWmpio1NVVlZWWuXgYAAAhgLo/ALFmyRK1bt9abb76p0aNHq6qqSitWrNC7776ryspKLVmyxKnzPPXUU5oxY4Zmz56tadOmady4cfrqq68kSaWlpdq1a5eWL18uSWrSpIl69+4tSSoqKtKtt96qDz74QImJiXrxxReVlpamVatWuXopABBw3L3PFBBsXA4wgwcP1jvvvKOoqChJUmhoqGbPnq2srCyNHj3aqQBTWlqq2bNny2q1SpIWLVqkyZMn29rXrFmj1NRUDRo0qNZzV61aJavVqsTEREnS2LFj1bx5cz3yyCNq166dq5cD4BIxT8a33LnhJaEHwcjlf3FWrlxpCy8XatOmjd555x2nzhEZGWkLL5JUUVGhBx54QJJUXV2tZcuW6Y477tCoUaP07bff2j13/fr16tKli+372NhYtWrVShs3bnT1UgB4wPl5Mnww+i9WfCEYuRxgwsPDlZuba/d16NAhbdy4UX/7299c7sCRI0f0xhtv6LHHHpMklZSUKD09XXPmzNGOHTvUs2dPbdmyxfb4rKwsxcXF2Z0jJiZGOTk5db5G586dZbVaZbVatWzZMpf7CMA3so8V2v3pjIa+nwwjYGioXC4htWvXThaLxWFbQkKC7r//fqfPtXv3bs2fP18ffPCBBgwYoH/84x9q0qSJRo8eLUmaM2eOhg8frokTJ+rAgQMKDQ1VQUGBYmNj7c4TERGh4uLiOl8nOztbzZpx51zAVy4MItd1uMKp53xztEDTMndJkqZl7tKaCdfp6lbR9T7PnZJKMGGlGBoqlyN7Wlqa9u3bpwMHDti+9u7dq2nTpunPf/6zS+e65ppr9P7772vdunXauXOn3n33Xbv2xo0ba+3atfrhhx9sIyyxsbG1Vh2VlJQoJibG1UsB4Ib6fuOvGUS+OVrg1Hk/yjqhvDPlkqS8M+X6OPuHBj+6AqBuLgeYBQsWqH379mrbtq3tq0OHDrrnnnv0q1/9yq1OjBgxQv369dOxY8dqtbVs2VLJycm20JKQkKD8/Hy7x+Tl5SkpKcmt1wbgmvrmvDgKIs4YnHylmkeFS5KaR4VrUOcWzN0AUCeXA0zN8o107lYCO3bs0AcffOB2Ry6//HJ169at1vGqqiqFh4era9eukqSUlBS7PWFOnjypkydP6pZbbnH7tQF4jqMgItU/v+XqVtFaOaaHJGnlmB51lo/cmScDIPi4HGBCQkIUGhpq93XZZZcpNTVV06ZNc+ochYWFev3111VYeO4foD179qi8vFw333yzPv74Yy1cuFBFRUWqrKzUo48+qsWLFys09NxvexMmTNC+fft04sQJSdIbb7yhKVOm2K1qAmAeR0HE2bJS5/hmdn/W5G55CkDwcXkS75133qn77rtPISH/yz7h4eFq27at4uPjnTrHiRMntHDhQs2ZM0f9+/dXQkKC3nrrLVv76tWrtXLlSl1//fV6/PHHbZvYSedWHL377ruaOXOm4uPjZbFYbBveAfAPNYOIo7LSxD7tXD6vp84DIPC5HGDS09PVoUMHh22nTp1yWGKqqWPHjjp48KDDtkGDBjmcC3OhXr16KTMzs/7OAvALg5Ov1O//dUh5Z8ptZSV3Vio5Og+AhqneElLNPV/CwsJqHcvNzdXBgwe1ePFiX/QZQICpWVaS5FYpyNl5Mp5aveToPKyMAvxDvSMwgwYN0v79++s9kWEYslgsbBQHwKELy0r/7+8Ha5WCnNnzpeZ56uKpvWEcnaeh7zsD+It6A0xaWpqmT5+uffv22SbSOlJZWamVK1d6tHMAglMwlYK4zxBgjnoDzHPPPafZs2fb5r1UVFSoUaNGDh9LCQmAM86Xgn6++vOLloI8xZshgxEZwBz1zoHp06ePfv3rX9u+v9hS6c8//9wzvQIQ9JwpBXmKow3xAnUuC/vgAOfUOwKTlJSk3/3ud2rdurVCQ0P1/fff6+9//7sMw7B7XHFxsVauXKmbbrrJa50FAE8JxJETd+8X5c6KL8Df1RtgHn74YS1YsEDz58/XsWPHZBiGNm7c6PCxdd3kEQBw6Rztg1NfgHE39AD+rt4SUqNGjbR06VLl5uaqsrJSEydOVHV1da2vkpISzZgxwxd9BgCbhlRSqes2DRfj7r2pAH/n8q0EpkyZ4vB4RESEHn300UvuEAA4q6HdWsDZfXAu5E7oAQKBywGmV69edba1aMEPBgDPqW90pSGOLrg6+dmd0AMEApcDDAD4gjOjK4wuOMeXK74AXyHAAPAb4WEhShvYUeFhIU6Nrnh6dCFQl1YDDREBBoDHXRhEXBERFqpZg3+kiLBQp0dXPDm64Gi/GE8gGAGeR4AB4HEXBhF3BdPcDXeDEcEHqBsBBoDfakhzNxyFFW+NCAHBgAADwCfcLSu5w1N7wzg6j7f2nSGsAK4hwADwCU+UlaT6g5Cn9oZxdB5nz03pB/A+AgyAgFJfEPLU3jCOzuPomKMRGWdGUwg5wKUhwAAIKo5WL7lTCnJ0nprH2sc1dnu0h5IRcGkIMABM4a05MTVXL0lyqxTkaBVUzWMH8ooa3E7AgL8gwAAwhafmxDhy4eolZ0tB9Z3H0TF2AgbMQ4ABENScKQW5GzyCaa8aINAQYAAENWdKQZcSPBrSXjWAPyHAAAh69ZWCAokv99MB/FmY2R0AADjv/NwhVxB6EIz4vxkA/o+3dtk1mzcnTANmIcAACGieGl3w1A6+jgRrMALMRIABENA8NbrgqR18a7qUYETwAepGgAHgt3w5d8Nbe7q4G4y8OSIEBAMCDAC/5cu5G97a08XZYFRztMVbI0JAsCDAAMD/8dTS6gvDiDPByNFoC7v8AhdHgAEAD3IURhwFowtDjqPRFm/v8svdsBHoCDAA4EHOlH5qhpwOzZs4HG3x5mZ73A0bgY4AAwAe5Ezpp2bIOZhfzD2VABcRYAAEHU+tXnJ0nvrO7Uzpx1HICdRbG9REaQq+YlqAycvL07BhwxQVFaVevXrp3//+t62tuLhYkydP1owZMzRu3DgdOnTI7rnbt2/XmDFjNHXqVM2bN0+GYfi6+wD8mKdWLzk6jzPnri+MBPNdrClNwVdMCzBPPfWUZsyYofXr16u4uFjjxo2ztY0fP159+/bVCy+8oLlz5+pnP/uZKioqJElHjx7VmDFjlJGRoYyMDJWXl+vpp5826zIABABv7ifj7rmDZcQFMIspAaa0tFSzZ8/WzTffrAEDBmjRokU6fPiwJGnnzp3auHGj7rzzTklS165dZbFYlJmZKUlaunSp+vXrp5iYGEnS3XffrfT0dBUVFZlxKQACgKNRE0+FGu4zBJjDlAATGRkpq9Vq+76iokIPPPCAJGn9+vVKTExURESErb179+7asGGDrb1Lly62tq5du6q4uFjbtm3zTecBBAWCBxDYTJ/Ee+TIEb3xxht67LHHJElZWVmKi4uze0xMTIxycnJUVFSk3Nxcu/awsDA1bdpUOTk5db5G586dZbVaZbVatWzZMu9cCAAEMCbfItCEmfniu3fv1vz58/XBBx9owIAB+sc//qGCggLFxsbaPS4iIkLFxcUqKDh3L5C62uuSnZ2tZs2oMwNAXc5Pvp3Utz2jUggIpo7AXHPNNXr//fe1bt067dy5U++++65iY2NVVlZm97iSkhLFxMTYgktd7QAAoGEwvYQkSSNGjFC/fv107NgxJSQkKD8/3649Ly9PSUlJioyMVHx8vF17SUmJzp49q6SkJF93GwAaPEpPMItfBBhJuvzyy9WtWzelpKRoz549dqMsWVlZuu222yRJKSkp2rFjh60tOztb0dHRuuGGG3zeZwBo6Nj3BWYxJcAUFhbq9ddfV2HhuZuZ7dmzR+Xl5br55pt17bXXqk+fPtq8ebOkc/NkQkJCNHr0aEnS9OnTtXXrVpWWlkqS1q5dq/nz5ysyMtKMSwEQRLy5XwwAzzJlEu+JEye0cOFCzZkzR/3791dCQoLeeustW/ubb76pWbNm6dNPP1V+fr42bdqk0NBzk8oSExOVkZGhSZMmKTY2Vi1atNBDDz1kxmUACDLnl1Z7mrvByJuB6sK7YV/X4QqnnlNWWaWMrfs1dUACE31hOlMCTMeOHXXw4ME62+Pi4rR27do624cMGaIhQ4Z4o2sA4HHuBiNvBaqad8NeM+E6tbmicb3PY6US/AnjpADgB3xZvqp5N+yPs3/w+msCnkaAAQAPcyeM+HJnYEd3w/aUC0tTgDcRYADAw5wJI2ZOGHb2btiuLpGuWZr65miBZzoMOECAAQATmH0vJmfuhu3qEmlKU/AlAgwAwCO8WZoCaiLAAAA8MnfF2dIU4AkEGABo4Dw5d8WZ0pTELQhw6QgwANDAmTF3hVsQ4FIRYACggWPuCgIRAQYAGrhLmbvCvi8wCwEGAOD03JULOTt3hvku8AYCDADAKTVHW5ydO+PN+S6Eo4aLAAMAqJej0RZvz51xJpwwGbjhIsAAAOrlaLTF2/u+uBNOGJFpOAgwAIB61TXa4s7cGW9iRKbhIMAAABy6cM6Lt0dbGDmBqwgwANAA1Xc3bEdzXrw52uLMyAkhBxciwABAA1Tf3bD98c7SlIdwIQIMAKAWdueFvyPAAABqcXfOS32lKX9FeSrwBNb/YQAAn3Fnzkt9panz/O0WBJSnAg8BBgDgU87egsBTGF0JTgQYAIBP+XqCsKdGVwhC/oUAAwDw6dyVQJ0gTJnJvxBgAABOz11xR835Lp7cFM/f5tLAdwgwAACvqWu+izMThOsLJ+7OpaEUFBwIMAAAt9VXenJ3vosz4cTRuZ0ZkfFmKYhw5DsEGACA2+orPbk738WZ4FPz3O3jGvt0dZMjzJPxHQIMAMBrnJ3vUnPkxJngU/PcB/KK/O72B/AeAgwAwCnurlSqb76Lo3KRo+DjqDx04bk9vbqJcpB/I8AAAJzirZVKdZWLLgwnzsyJ8eTqJolykL8jwAAATOXMyImzk4Hduf0BAhMBBgBgKmdGTry9+R37yQQeAgwAwCFf7s5b38iJp8tDF/L1vZngGaYFmEOHDumWW25RVFSUrr32Wm3dutWu/YknnpDFYrF97dixw9a2d+9ejR07VqmpqUpNTVVZWZmvuw8AQc+bu/O6w93yUH2jK76+NxM8w5QAYxiGJk+erEGDBikjI0OGYWjYsGHav3+/JKm0tFS7du3S8uXLtXz5cq1evVq9e/eWJBUVFenWW2/VE088oRdffFGdOnVSWlqaGZcBAPBzzoyuOFueoszkX0wJMF9//bXmzp2rOXPm6J577tGHH34oi8Wi999/X5K0Zs0apaamaubMmZo5c6YmT55se+6qVatktVqVmJgoSRo7dqxefvllfffdd2ZcCgDAjzkzuuJMeYoyk/8xJcAkJSVp4MCBtu9btGih5ORkRUREqLq6WsuWLdMdd9yhUaNG6dtvv7V77vr169WlSxfb97GxsWrVqpU2btzos/4DAC6Nr+bXODu6Ul95ijKT/zElwISHh9c6dvr0aQ0bNkwlJSVKT0/XnDlztGPHDvXs2VNbtmyxPS4rK0txcXF2z42JiVFOTo7X+w0A8Axfza9xdkO8+lBm8j9+sQrpk08+0ciRI2W1WtWkSRONHj1av/rVr/Sf//xHP/nJTzRx4kRVVZ3bCbGgoECxsbF2z4+IiFBxcXGd5+/cubOsVqusVquWLVvm1WsBAPgXVzfEc4Qyk/8JM7sDFRUVyszM1IoVK2q1NW7cWGvXrlX79u2Vk5Ojzp07KzY2ttaqo5KSEsXExNT5GtnZ2WrWjE2NAKChc1QKcnZJtjtlJk8u94Y900dgnnvuOc2bN08REREO21u2bKnk5GRbaElISFB+fr7dY/Ly8pSUlOT1vgIAXOfN+S6untuVDfFcLQd5e7M92DM1wLzyyisaMmSIWrduLUkqKyuzlYrOq6qqUnh4uLp27SpJSklJsdsT5uTJkzp58qRuueUW33UcAOA0d+e7OBNOXD23sxviuVMO8uZme6jNtACTkZGhAwcO6Pjx49q8ebPWrVunSZMmacuWLVq4cKGKiopUWVmpRx99VIsXL1Zo6Ln/OSdMmKB9+/bpxIkTkqQ33nhDU6ZMkdVqNetSAABe4E7wcSb0OLMhnrurjrgXk++YMgfm1Vdf1bRp0yRJS5YssR2fPn26QkJCtHr1aq1cuVLXX3+9Hn/8cdsmdtK5FUfvvvuuZs6cqfj4eFksFi1fvtzn1wAA8D/nQ8+lGpx8pX7/r0PKO1NOOchPmRJgJk6cqIkTJ9bZfuzYsYs+v1evXsrMzPR0twAAJvHlfZeccb4c9PPVn9vKQWdKK9w6V1lllTK27tfUAQl+c1uGYOAf/6cAABo0Z8pFvg45zpSDyiqrtPyjHJVVVtX5mPLKaj2/Za/KK6sv6TywR4ABAAQEf7u5pORcOPHUeQg59ggwAAAEAE+FpWBBgAEANHj+NgcH9eOdAgA0eL4sT1EK8gwCDAAgqPnb6AqlIM8w/V5IAAB4k6f2hvG3INTQEWAAAHCCp4IQPIMYCQBAkGhI82sIMACABsWf7o7taQ1pfg0BBgDQoHhzxZGnzt2QRlLcRYABAMCLso8V2v3pDGdHUhpy0CHAAADggDPloPrCyTdHCzQtc5ckaVrmLn1ztMCt89SlIZWMaiLAAADgQH3lIGfCyUdZJ5R3plySlHemXB9n/+DWeST3Q46n+NtoDwEGAAA3OBNOBidfqeZR4ZKk5lHhGtS5hVvncTbkeJO/jfYQYAAAcIMz4eTqVtFaOaaHJGnlmB66ulW0W+dxJuQ0NAQYAADc4Ew4kaTO8c3s/qxZCvJUyGloCDAAALipZjipT12loPrO42xYcoajuSz+Nr/FGQQYAAB85FJKQa6Gpbo4msvib/NbnEGAAQDAR8woBQXi6IozCDAAAPiIJ0tBjjhaah2IoyvOIMAAAOBDnioF1eQPS619iQADAICfcXbTugsf58ml1s6UncwuTRFgAADwI86OpNR8XIfmTTw2v8aZspPZpSkCDAAAfsTZkZSajzuYX+zz+TVmIsAAAOBHnF2p5OhxDWl+DQEGAAA/4uxKJW+vaLqQP97KgAADAIAXhYeFKG1gR4WHOf+R6+xIirdGXGryx1sZEGAAAHCTM+EkIixUswb/SBFhoZd0Hk9xZy6LL0d7nEWAAQDATc6EE0+dxxMhp665LM6EGl+N9jiLAAMAQADwRFhyNJfFHyfoOoMAAwBAELtwdMXRXBZ/nKDrDAIMAAA+5Mv5LjVHVyTVmsvijxN0nUGAAQDAhzw1b8YZjkZXas5lcTRB15k5MWZvbEeAAQAgSDk7unJhqHFmTow/zJsxLcAcOnRIt9xyi6KionTttddq69attrbi4mJNnjxZM2bM0Lhx43To0CG7527fvl1jxozR1KlTNW/ePBmG4evuAwDg99xZ/uzMnBh/mDdjSoAxDEOTJ0/WoEGDlJGRIcMwNGzYMO3fv1+SNH78ePXt21cvvPCC5s6dq5/97GeqqKiQJB09elRjxoxRRkaGMjIyVF5erqefftqMywAAwK84ml/j6vJnZ0Zt/GHejCkB5uuvv9bcuXM1Z84c3XPPPfrwww9lsVj0/vvva+fOndq4caPuvPNOSVLXrl1lsViUmZkpSVq6dKn69eunmJgYSdLdd9+t9PR0FRUVmXEpAAD4DU/Mr3Fm1MYfNrYzJcAkJSVp4MCBtu9btGih5ORkRUREaP369UpMTFRERIStvXv37tqwYYMkaf369erSpYutrWvXriouLta2bdt81n8AALzJ2ZVK3lrR5Myojdkb25kSYMLDw2sdO336tIYNG6asrCzFxcXZtcXExCgnJ0dFRUXKzc21aw8LC1PTpk2Vk5NT5+t17txZVqtVVqtVy5Yt89yFAADgBc6OpPhyRZO/CTO7A5L0ySefaOTIkbJarSooKFBsbKxde0REhIqLi1VQcG6Wc13tdcnOzlazZv6x9TEAAIHIl/vXOMP0AFNRUaHMzEytWLFC0rlwcvbsWbvHlJSUKCYmxhZcysrKHLYDAADvOD/a4y9Mj1HPPfec5s2bZ5vzkpCQoPz8fLvH5OXlKSkpSZGRkYqPj7drLykp0dmzZ5WUlOTTfgMAAPOYGmBeeeUVDRkyRK1bt5Z0bmRl5MiR2rNnj90oS1ZWlm677TZJUkpKinbs2GFry87OVnR0tG644Qbfdh4AAJjGtACTkZGhAwcO6Pjx49q8ebPWrVunSZMm6dprr1WfPn20efNmSdLu3bsVEhKi0aNHS5KmT5+urVu3qrS0VJK0du1azZ8/X5GRkWZdCgAAAcPRXBZ/m9/iDFPmwLz66quaNm2aJGnJkiW249OnT1dISIjefPNNzZo1S59++qny8/O1adMmhYaem2GdmJiojIwMTZo0SbGxsWrRooUeeughMy4DAICA42gui7/Nb3GGKQFm4sSJmjhxYp3tcXFxWrt2bZ3tQ4YM0ZAhQ7zRNQAAgkogjq44w/RVSAAAwHsCcXTFGcEVxwAAQINAgAEAAHacKTuZXZqyGIZhmPLKPlBYWKjo6GgVFBSwEy8AAAHCmc9vRmAAAEDAIcAAAICAQ4ABAAABhwADAAACDgEGAAAEHAIMAAAIOAQYAAAQcAgwAAAg4BBgAABAwCHAAACAgEOAAQAAAYcAAy1btszsLuACvB/+hffDv/B++Bcz3w9u5ghZrVYdOXLE7G7g//B++BfeD//C++FfvPV+cDNHAAAQlMLM7oA3nR9cKiwsNLkn/s0wDP6O/Ajvh3/h/fAvvB/+xVvvx/lzXqxIFNQlpCNHjqh169ZmdwMAALjh8OHDslqtDtuCOsBUV1fr2LFjioqKksViMbs7AADACYZh6MyZM4qPj1dIiOPZLkEdYAAAQHBiEi8AAAg4BBgAABBwCDAAACDgBPUyavzPoUOHdN999+mzzz5TYmKili1bpgEDBkiSiouLlZaWpssuu0ynTp1Senq62rZta3KPG45NmzZp2rRpOnjwoCTeD7MVFBTo9ddfV2xsrFq1aqUBAwbwnpjg7Nmzmjt3rpo1a6bKykqVlZXp17/+tSIiIiRJVVVVevjhh1VSUqL//ve/mjt3rrp3725yr4NPfn6+nn76aRUVFSkjI8N2vL6fie3bt2vFihWKiYnR5ZdfrvT0dM8vpjEQ9Kqrq43BgwcbzzzzjPG73/3O6N69u9G4cWNj3759hmEYxqhRo4zXX3/dMAzD2LNnj9G1a1ejvLzczC43GKdOnTISEhKMtm3b2o7xfphn3759xl133WXk5+fbHec98b2pU6cazz33nO37adOmGfPnz7d9P3v2bGPRokWGYRjGiRMnjA4dOhinT5/2dTeDWnl5ufGnP/3J6Nixo3HvvffatV3sZ+LIkSNG+/btjVOnThmGce69WrJkicf7R4BpAPbs2WN8/PHHtu9PnDhhNGnSxFixYoXxxRdfGJGRkUZpaamtvVu3brb/MeFdU6dONRYtWmQLMLwf5jl16pTx05/+1Dh+/Ljdcd4TcyQnJxvvvPOO7ftVq1YZt99+u2EYhvH9998bjRo1Mr777jtb++23324LNPCssWPH2gWY+n4mUlNTjQkTJtjadu3aZTRt2tQ4e/asR/vFHJgGICkpSQMHDrR936JFCyUnJysiIkLr169XYmKibVhWkrp3764NGzaY0dUG5c9//rN69+6tNm3a2I7xfpjnV7/6lTp06KAlS5bo+uuv18MPP6zKykreE5N069ZNGRkZqq6uliR98cUXuv/++yVJH3zwgSIiIuxKFrwn3tOoUSO77+v7mVi/fr26dOlia+vatauKi4u1bds2j/aLANMAhIeH1zp2+vRpDRs2TFlZWYqLi7Nri4mJUU5Ojq+61yDl5+dr48aNmjBhgt1x3g9zFBcX65VXXlGPHj20fPlyrVmzRi+99JLmzp3Le2KSZcuW6cCBAxo/frw+++wz3XHHHRo6dKgkfk7MdrG//6KiIuXm5tq1h4WFqWnTph5/fwgwDdAnn3yikSNHymq1qqCgQLGxsXbtERERKi4uNql3DcOCBQv05JNP1jrO+2GO7du3q6SkRHfffbcsFos6deqkO++8Uy+99BLviUmuuuoqvf3224qOjtatt95qG4mR+Dkx28X+/gsKCiTJJ+8PAaaBqaioUGZmpp544glJ5/4nKysrs3tMSUmJYmJizOheg/D222+rf//+atmyZa023g9zHD9+XJLUpEkT27H+/furpKREVVVVvCcmyM7O1j//+U+tXLlSzz77rEaOHKm//OUvkvg5MdvF/v7PBxdfvD8EmAbmueee07x582y1y4SEBOXn59s9Ji8vT0lJSWZ0r0FYtWqV7rrrLlksFlksFv3iF7/QoUOHZLFYdPjwYd4PE0RFRUk6V1o9r0WLFpKka665hvfEBJMmTVLfvn1t/52amqqHHnpIEv9ume1if/+RkZGKj4+3ay8pKdHZs2c9/v4QYBqQV155RUOGDLHdobusrEwjR47Unj177NJyVlaWbrvtNrO6GfRefvllffnll7avRYsW6aqrrtKXX36ppUuX8n6Y4Cc/+YkaNWqkf/3rX7Zjp06d0rXXXqtx48bxnpjg3//+t938vV/+8pc6fPiwJGn48OE6deqU7XuJ98SXUlJSLvozkZKSoh07dtjasrOzFR0drRtuuMGj/SDANBAZGRk6cOCAjh8/rs2bN2vdunWaNGmSrr32WvXp00ebN2+WJO3evVshISEaPXq0yT0OXomJibrmmmtsX23atFF4eLiuueYa9e3bl/fDBFdccYXS0tL00ksv2Y796U9/0mOPPcbPiEmGDh2qDz74wPb9t99+q2HDhkk6Nzo2YcIErVu3TtK5EmBWVpamTp1qSl+DXVVVld0cpPp+JqZPn66tW7eqtLRUkrR27VrNnz9fkZGRHu0Xd6NuAF599VX98pe/rHV8+vTpeuGFF5Sfn69Zs2bpqquuUn5+vhYvXqz4+HgTetowrVmzRgsXLtR3330nSbwfJqmsrNTDDz+s4uJiNWvWTG3bttW0adMk8Z6YoaCgQA899JDi4+PVtGlTHT9+XI8//riio6MlnStLpKWlqVmzZsrPz9fDDz+s5ORkk3sdfN5++209/PDDCgkJ0TPPPGMLKfX9TGzevFlvvPGGYmNj1aJFC82fP9/jfSPAAACAgEMJCQAABBwCDAAACDgEGAAAEHAIMAAAIOAQYAAAQMAhwAAAgIBDgAEAAAGHAAMAAAIOAQZA0MnJydFrr71m2uu/+OKLtjtcA/AOAgyAoHLw4EGtXr1av/jFL+yOf/TRRxoxYoQmTJigKVOmaPLkyVq8eLEWLlxY7zlzcnJsdxDv0KGD/vKXv9i1f/jhh4qPj9f111+vrKwsTZkyRY8++qjd3a0BeBa3EgAQNAzD0KBBg5SZmakWLVrYjj/zzDN69dVXtXnzZrVr1852fPPmzVqwYIHdnXPrUllZqQ4dOqhTp0768MMPa7UPHz5cL730ku1+MDt37tSzzz6rzMzMS78wALUwAgMgaLz33ntq2rSpXXj561//qkceeUSvv/66XXiRpCFDhmjSpElOnTssLExTpkzRtm3bdOLECbu2kydPKjY21u5mdj179tSBAwe0e/dut68HQN0IMABqOX36tGbPnq17771X3bp10yuvvCJJ+vvf/67bb79dTzzxhH7729+qefPm6tGjh/Lz823P3bNnj6ZNm6aUlBT16tVLn332maqqqpSZmalevXrpT3/6k3784x/rxz/+sSRp27Ztuv/++/Xcc88pPT3ddp4vvvhCVqtVoaGhWrdunSSpqKhIo0aN0siRI1VSUlKr36tXr9aNN95od2zevHnq1auX7fVqmjJlit33q1atUmpqqn76059q3LhxKiwstLVNmjRJFotFL7/8st1zfve73+mee+6pde4bbrhBL730ksPXBXCJDACoYfz48cbx48cNwzCMLVu2GBaLxdixY4dx5MgRo0OHDka/fv2MTZs2Gbm5uUZ8fLyxdOlSwzAMo7i42Bg1apRRVVVlGIZhLFiwwGjRooXx3//+1/jLX/5iSDLGjRtnvPvuu8aDDz5onDhxwmjWrJnxn//8xzAMw5g9e7aRlJRkDBgwwCgqKjLefvttIyQkxMjPz7f1bcqUKcbhw4dr9bmiosKIjIw0NmzYYDt27NgxQ5IxefLkWo//+OOPjYkTJxp33nmnkZ6ebpw5c8Z47733jIyMDMMwDKOkpMTo2LGjkZqaave8u+66y7BarUZlZaXt2PDhw43q6upar7Fy5Uqjffv2zv2lA3BJmNkBCoB/2b9/v7Zt22ZbxVNRUaGbbrpJubm56tWrl1q1aqUbb7xRQ4cOlXRulGH//v2SpMzMTOXn5+uZZ56RJJ09e1Zdu3ZVfn6+Bg0aJEn6+c9/rp/97GcaOXKk1q1bp8LCQrVv316S9NOf/lSvvPKK9uzZo/DwcI0aNUqtWrXSG2+8obS0NFVVVen06dOyWq21+r1v3z6VlpaqZcuWtmOHDh2SJF1xxRW1Hj9w4ED94x//0MKFC/Xiiy+qadOmevbZZ9W7d28tXbpU0rkyUEVFhd3zpk2bprfeekvr16/XiBEj9Pnnn+u6666TxWKp9RqtWrXSwYMHVVRUpCZNmrjwLgCoDwEGgJ3s7Gw1adJEjzzyiO3YggULbP8dEhKikJD/VZ8bN26s8vJySVJWVpaSk5PtnltTVFSU7b/j4uIkSfn5+YqPj1fjxo1VUFCgoqIihYeHKzQ0VJMmTdKqVauUlpamv/71rxo4cKDD8548edLWn/OaNm0qSbXmrJx35ZVX2vUjKytLzz77rK677ro6+9+nTx9169ZNK1eu1IgRI/T666/b/f1c6Hxf8vLyCDCAhzEHBoCdiooKHTx40G7uh/S/gFDfc//973/bHausrKx1rvP69u2rUaNG6Y9//KMk6fPPP9fgwYMVExNje8ykSZO0b98+bd26VevWrVNKSorDc1122WW2PpzXuXNnNW/eXFu3bq2373X139F1P/DAA/rrX/+qnTt36syZM3aTdx0JDQ116vUBOI8AA8BOUlKSSktL9fzzz9uOHTp0SJ999lm9z+3UqZM+++wz/f3vf7cde+edd2wjNDVZLBZ169ZNX3/9tV544QVFRkbaJuyeFx8fr2HDhun5559XcXGxYmNjHZ7rfOnowrAUGhqquXPn6rvvvtMf/vAHp/q/YsUKlZWV2Y45WgY9btw4RUVFKSUlRT//+c/rPN/Zs2cVGhpqtyoKgGcQYADYSU5O1uDBg7Vw4UJNnz5dv/nNb5Samqqbb75Z0rlRiurqarvnVFVVSZLuvvtuxcXFafjw4UpPT9fixYv12WefKS4uTsb/bTl1YZjJycnRCy+8oOuvv15XXnmlOnbsqK+++sr22PPuv/9+vffeexoyZEid/Y6Pj1ebNm10+PBhu+OzZs3SvffeqwceeMAuWElSQUGB3fdpaWnKysrSgAEDlJGRobFjxyo5ObnWazVt2lT33HOPysrKdNttt9XZp0OHDqlLly6KiIio8zEA3GT2LGIA/uf77783br31VuOyyy4zrr/+etsqoQ0bNhhNmzY1evToYWzfvt3YunWr0a5dO6N169bGJ598YhiGYfzzn/80unXrZjRt2tS48847jcLCQqO0tNRYvHixIckYPHiw8c033xiGYRhnz541rrvuOiMuLs4IDw83JBmSjClTptj1p7q62rj66quNoqKii/Y7NTXVeOSRRxy2/eEPfzBuvvlmY+TIkcZ9991njBo1yrjzzjuNf/zjH3aPS09PN6644gqjVatWxm9/+9s6X+ubb74x5s2bd9H+TJo0yVi4cOFFHwPAPezEC8A0+/fv1zvvvKO5c+dKOreT7smTJ5Wamqq33nrL9rjy8nI9+OCDevHFF+s939ixY7V9+3av9ttZXbp00ccff2y3MgqAZ1BCAmCahx56SH369LF9b7FYFBcXpx49etg9bv369Rct1ZyXkJCgIUOGONzq39c2bdqk8ePHE14ALyHAADBNXl6e5s2bp88++0yFhYU6duyYnnvuOV177bWqrKzUvffeq8cee0yvv/76Ree/XOjxxx/X1q1blZeX5+Xe1+348eP65z//qTlz5pjWByDYUUICYJqjR49q3rx5+stf/qKCggL17NlTTzzxhG666SYVFxerW7duioqK0h/+8Ad17NjR6fNWV1fro48+0i233OLF3tftvffe0+233+5wczsAnkGAAQAAAYcSEgAACDgEGAAAEHAIMAAAIOAQYAAAQMAhwAAAgIBDgAEAAAHn/wNCt/Cz8ERkUQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def gauss(x, amplitude=1, mu=0, sig=1):    \n",
    "    return amplitude * np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))\n",
    "\n",
    "def linear(x,*params):\n",
    "    a, b = params\n",
    "    return a*x + b\n",
    "\n",
    "np.random.seed(43)\n",
    "energy = np.linspace(10,100,100)\n",
    "spectral_line1 = gauss(energy, *[50, 42.1, 3])\n",
    "spectral_line2 = gauss(energy, *[100, 68.7, 1])\n",
    "spectral_line3 = gauss(energy, *[100, 18.1, 0.5])\n",
    "counts = (linear(energy,-3, 500) + #continuum spectrum\n",
    "spectral_line1 + #line emission\n",
    "spectral_line2 + #line emission\n",
    "spectral_line3 + #line emission\n",
    "np.random.normal(0, 7, size=len(energy))) #statistical noise\n",
    "\n",
    "# for this example, we will say that the errors on our measurements have to \n",
    "# do with the counting errors, so we will use Poissonian statistics here:\n",
    "\n",
    "errors = np.sqrt(counts)\n",
    "\n",
    "# Note that this assumes that the response of the detector is uniform across \n",
    "# all energies; in reality this is not the case; These would get tied up \n",
    "# in **systematic** uncertainties and the instrument response function.\n",
    "\n",
    "plt.figure(figsize=[8,6])\n",
    "plt.errorbar(energy,counts,errors, fmt='.')\n",
    "plt.xlabel(\"energy (GeV)\")\n",
    "plt.ylabel(\"flux (counts)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "465e0d20",
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss_function(parameters, data_x, data_y, model, verbose=False):\n",
    "    '''\n",
    "    This is a really generic loss function.\n",
    "    It can take in any number of parameters, any generic model.\n",
    "\n",
    "    Notice that the parameters are passed *first*.\n",
    "    This is because of the way scipy's libraries need the function to be formatted.\n",
    "    '''\n",
    "\n",
    "    loss = pow(residuals(parameters, data_x, data_y, model), 2.).sum()\n",
    "    if verbose:\n",
    "        print(loss)\n",
    "    return loss\n",
    "\n",
    "#guess = np.array([50. , 312.])\n",
    "#minimize(loss_function, guess, args=(data_x, data_y, model))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "59380684",
   "metadata": {},
   "outputs": [],
   "source": [
    "def calc_chisquare(meas, sigma, fit):\n",
    "\n",
    "     diff = pow(meas-fit, 2.)\n",
    "     test_statistic = (diff / pow(sigma,2.)).sum()\n",
    "\n",
    "     return test_statistic"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "8f55a039",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Continuum slope               :  -2.97 +-   0.03\n",
      "continuum intercept           : 498.34 +-   1.98\n",
      "peak0_A                       : 112.75 +-   9.68\n",
      "peak0_mu                      :  18.10 +-   0.08\n",
      "peak0_sigma                   :   0.47 +-   0.05\n",
      "peak1_A                       :  56.56 +-   3.97\n",
      "peak1_mu                      :  42.04 +-   0.23\n",
      "peak1_sigma                   :   2.96 +-   0.24\n",
      "peak2_A                       :  98.08 +-   6.36\n",
      "peak2_mu                      :  68.56 +-   0.07\n",
      "peak2_sigma                   :   1.00 +-   0.07\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAG4CAYAAABb+t1HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAAxOAAAMTgF/d4wjAABo3UlEQVR4nO3de3zOdR/H8dd17XDtaLYYZuXc2Mi5s6QoRQeniiSJiDklKVEq0pHq1pJUuhUdVYoUIh0lh9S2ksgcwoY2dr52/e4/lut22TW7tl3btV17Px+PPWq/0/W5dsXefY8mwzAMRERERKoRs6cLEBERESktBRgRERGpdhRgREREpNpRgBEREZFqRwFGREREqh1fTxdQkWw2GwcOHCA0NBSTyeTpckRERMQFhmFw/PhxoqKiMJudt7V4dYA5cOAAZ599tqfLEBERkTLYu3cv0dHRTs95dYAJDQ0FCn8AtWrV8nA1IiIi4oqMjAzOPvts++9xZ7w6wJzsNqpVq5YCjIiISDVzpuEfGsQrIiIi1Y4CjIiIiFQ7CjAiIiJS7SjAiIiISLWjACMiIiLVjgKMiIiIVDsKMCIiIlLtKMCIiIhItaMAIyIiItWOAoyIiIhUOwowIiIiUu0owIiIiEi1owAjIiIi1Y4CjIiIiFQ7CjAiIiJS7SjAeLFcawFzV+8g11rg6VJERETcSgHGi+VZbTy/9g/yrDZPlyIiIuJWCjAiIiJS7SjAiIiISLWjACMiIiLVjgKMiIiIVDsKMCIiIlLtKMCIiIhItaMAIyIiItWOAoyIiIhUOwowIiIiUu0owIiIiEi1owAjIiIi1Y4CjIiIiFQ7CjAiIiJS7SjAiIiISLWjACMiIiLVjgKMiIiIVDsKMN4q8UN2fjoXgOQDGR4uRkRExL0UYLxU4if/YdRPDQAYs3QLifvTPVyRiIiI+yjAeKnVubEcIgKA1ON5rEk+7OGKRERE3MejAWbDhg2YTCb719NPP20/9+ijjzqc27Rpk/3cH3/8waBBg4iPjyc+Pp7c3FxPlF91FVjpYd1APY4CUDfUn+6tIj1clIiIiPv4evLFFy1axNy5c+3f33777QDk5OSwZcsW+7ng4GA6d+4MQGZmJtdeey2fffYZzZs3Z968eYwfP5758+dX/huoqnL+Ic5nDwtrLWHckT48ObAfcQ3DPF2ViIiI23gswCQnJxMdHc2ECROKnFu0aBHx8fF07969yLn58+cTHR1N8+bNARg0aBB169bl/vvvp3HjxhVcdTWRdRR8A2l1lpmO/+ykVVQtT1ckIiLiVh7rQnrmmWd4/PHHueyyy1i5cqX9uM1mY86cOdx4443069eP33//3eG+5cuX07p1a/v3ERERNGzYkBUrVlRa7VVe9lEIisAIjKA2JzxdjYiIiNt5LMD079+fZ555htzcXHr16sXs2bMByM7OZtasWUyePJlNmzbRsWNH1q5da78vKSmJOnXqODwrPDycHTt2FPtarVq1Ijo6mujoaObMmVMxb6gqyToKgREYgeFEmI57uhoRERG381gX0jXXXMM111zD+PHjmT59Og899BB9+vShZcuWDBgwAIDJkydzww03MGzYMHbt2oWPjw/p6elEREQ4PMtisZCVlVXsayUnJ1OrVg3qRsk+CoG1/22BSfF0NSIiIm7n8WnUJpOJxx57jNatW7Nu3TqHc0FBQSxevJjDhw/bW1giIiKKzDrKzs4mPDy80mqu8rKPFXYhBYQRZsr0dDUiIiJu5/EAA4Uh5oorrnA6Hbp+/frExsbazzVr1oy0tDSHa1JTU4mJiamUWquFnHQICAP/EILRFHMREfE+VSLAAPz999/07NmzyPGCggL8/f1p06YNUDh25tQ1YY4cOcKRI0e4+uqrK63WKi8vC/xDMPxDCDZle7oaERERt/NIgElJSWHs2LGkpBSOz3jnnXeIi4ujZcuWrFmzhhkzZpCZmYnVauWBBx5g5syZ+Pj4ADB06FB27tzJoUOHAHjzzTcZOXIk0dHRnngrVVN+JvgFYfiHEEKOp6sRERFxO48EGF9fX7766itiY2Pp2rUrVquVBx980H5+wYIFNG7cmD59+jBgwACuvPJK+7nw8HA++OADJkyYwKRJk9i7d6/DYnjCvy0wQeAfQpACjIiIeCGPzEKKiopi+/btTs91796dAwcOnPH+Tp06sXTp0ooozTvkZ4FfMIYlhGCTAoyIiHifKjMGRtwoLxP8gzD8gglBY2BERMT7KMB4o/ws8CvsQgow5YPN6umKRERE3EoBxhvlZYF/MIZ/yL/fay0YERHxLgow3ig/E/wCwT8YAFOe9kMSERHvogDjjfIKB/FiMnPCCFCAERERr6MA443y/51GDWRjKWyRERER8SIKMN7GMArHvPj9G2AMf0z5mokkIiLeRQHG2xTkgVFgH/9S2AKjACMiIt7FIwvZSQU6OePILwh/HzNhtWrhZ9NidiIi4l0UYLxNflbhP/2CsJjN1K9zFti0I7WIiHgXdSF5m7ws8A0E878frV/g/0ONiIiIl1CA8Tb5mfYZSMC/AUZjYERExLsowHibk2vAnOQXpBYYERHxOgow3uaUNWAAtcCIiIhXUoDxNqesAQP82wKjACMiIt5FAcbb5GfZ14ABNIhXRES8kgKMtynSAqMuJBER8T4KMN6myBgYDeIVERHvowDjbYrMQgosPCYiIuJFFGC8TZF1YDSIV0REvI8CjLdxOgZGLTAiIuJdFGC8TZ6zWUhqgREREe+iAONt8rUOjIiIeD8FGG+T52wlXnUhiYiId1GA8Tb5zvZCUguMiIh4FwUYb5PnbDfqLDAMz9UkIiLiZgow3iY/q+gYGAyw5nqsJBEREXdTgPE2zmYhgcbBiIiIV1GA8Tanz0LyPRlgNA5GRES8hwKMtzm9BcZsLgwxCjAiIuJFFGC8SUE+2PIdW2Dg34G8mZ6pSUREpAIowHiTvH9Div/pAUZTqUVExLsowHiTkwN1T10HBgoDTZ5aYERExHt4NMBs2LABk8lk/3r66acByMrKYsSIEYwbN47BgwezZ88eh/s2btzIwIEDGT16NFOnTsXQGieF8rLAxx98fB2PazVeERHxMr4lX1JxFi1axNy5c+3f33777QAMGTKE66+/niFDhvDLL79w3XXXsXnzZvz8/Ni/fz8DBw5k8+bNhIeHc++99/Lkk09y//33e+ptVB35mY4DeE/yC1YXkoiIeBWPBZjk5GSio6OZMGGCw/HNmzezYsUK3nrrLQDatGmDyWRi6dKlDBkyhCeeeIKuXbsSHh4OwK233spll13G2LFjCQ528su7JsnLKtp9BOpCEhERr+OxLqRnnnmGxx9/nMsuu4yVK1fajy9fvpzmzZtjsVjsx9q2bcsnn3xiP9+6dWv7uTZt2pCVlcX69esrrfYqKz+z6ABeUBeSiIh4HY8FmP79+/PMM8+Qm5tLr169mD17NgBJSUnUqVPH4drw8HB27NhBZmYmKSkpDud9fX0JCQlhx44dxb5Wq1atiI6OJjo6mjlz5lTMG6oK8jKLTqGGf7uQFGBERMR7eKwL6ZprruGaa65h/PjxTJ8+nYceeog+ffqQnp5ORESEw7UWi4WsrCzS09MBij1fnOTkZGrVquX+N1HV5GWBf0jR4/5BhedERES8hMenUZtMJh577DFat27NunXriIiIIDfXcePB7OxswsPD7cGluPM1XrFdSEFqgREREa/i8QADhSHmiiuuIDc3l2bNmpGWluZwPjU1lZiYGAICAoiKinI4n52dzYkTJ4iJianssquevKxiupAUYERExLtUiQAD8Pfff9OzZ0/69+/P9u3bHVpZkpKS6NWrF1A4dmbTpk32c8nJyYSFhXHJJZdUes1VTr66kEREpGbwSIBJSUlh7NixpKSkAPDOO+8QFxdHy5Ytad++PZdeeimrVq0CYNu2bZjNZgYMGADA2LFjWbduHTk5OQAsXryYadOmERAQ4Im3UrXknVAXkoiI1AgeCTC+vr589dVXxMbG0rVrV6xWKw8++KD9/JIlS3j//fe57777eOGFF1i5ciU+Pj4ANG/enISEBIYPH864ceM466yzuPfeez3xNkqUay1g7uod5FoLKuW1ft51AKtvYNGTCjAiIuJlPDILKSoqiu3btxd7vk6dOixevLjY8z179qRnz54VUZpb5VltPL/2D4Z3aYLF16fCX2vn/sO0bNGi6IeqLiQREfEyVWYMjJRfILkYzlbiVQuMiIh4GQUYLxJMjsbAiIhIjaAAU0aVOb7FVYGmXAxn06j9g9WFJCIiXkUBpoxOjm/Js9o8XYpdELnON3PUXkgiIuJlFGAqUUW32gSSi+GvMTAiIuL9FGAqUUW32gSdqQvJZgVrXoW8roiISGVTgPEiweQUs5XAv2vD5GdWbkEiIiIVRAHGTTw+qNcw/p1GXcwsJID87MqtSUREpIIowLiJxwf1WrPxMxWApVbRc2Yf8A3QTCQREfEaCjBewpR7HADDWYCBf2ciqQtJRES8gwJMFVPmrqjcDLINf/Dxc37eL1hdSCIi4jUUYKqYsnZFmfKOcwInGzme5B8EeWqBERER76AA4yVMucfJMJwM4D1Ji9mJiIgXUYDxEqbcklpgQtQCIyIiXkMBxlvkZnDcOEOAsYSSn51e5fZvEhERKQsFmAqUfCDD4Z8VyZSbwXHO0IVkqYUtK73K7d8kIiJSFgowFSRxfzpjlm4BYMzSLSTuT6/Q10tOzWeZ9dLiw1JALcit+CAlIiJSGRRgKsjqpEOkHi/ceyj1eB5rkg9X2Gsl7k9n5LbGrDY6Fx+WLLUwKcCIiIiXUIApq4Izb4zYI7YedUP9Aagb6k/3VpEVVsrqpEMczrMAZwhLAQowIiLiPRRgyiL3BH/95zouMv1K8p4DTi+JaxjGiwM7APDiwA7ENQyrsHJ6xNYj0rdwhlGxYclSy75ar4iISHWnAFMGiWkFDMu7l++N1oxd9DWJn85zOkW5VVQth39WlLiGYbxS/2Ni2FN8WAoIg7yyBRiPb1QpIiJyGgWYMliddIjULAOAQ0Y4axIPEPxyZ4b4fA7W3GLvczYryV3hoLXxO5Hm9OLDUjnGwHh8o0oREZHTKMCUQZHxLbfeR+7Vz3Crz1qCF1wIW9+CAqvDPcXNSnJXODBlHeWoEVr8BRoDIyIiXkQBpgyKjG+Jro313Gu5Ju8JcrtOhQ1PwUsX4fvbckwUBpPyzEoqsZXGVoAp+yjHzhRgNAtJRES8iAJMGZ0+viX5QAY2zPwSfjXE/wQX3o1lzVSW+0/DZ9daerSKLPOspBJbabL/wYTBUc7cAkPuCXugOknjW0REpDpSgHGDIt1DB7Og0zAyR25iecHFBC4fRdwXA3npCl+gAmYlZaZi+AaSg6X4ayy1MGEQQo7DYY1vERGR6kgBxg2K7R7yC+SVgt6cuHsLNLmMjusG85rfU8SZ/3JvAccPYIQ2OPM1/iEYJjO10IaOIiJS/SnAuEGJi9ZZQqHbA2SO2sxOoyFBb/aC94ZC2h9FnlWm/ZMyDmALjTrzNWYzRmAE4abyrwWjbicREfE0BRg3cHXROiPoLB633krmXRshoDa8dAmWleNpSCrg+v5JRQJExgGMWg1LrNMIqkMdU/kH8qrbSUREPE0Bxk1Ks2idUSsKrnsORn+PqSCXtZZ7sayeyuqfd7k0U6lIgEjfV3ILDIUB5iw0E0lERKo/BRhPOqsZOdfN54a8xzBl7KXHppFE+hcGmFLNVDq6C1vtxiVeZgTX4SxTxe6KLSIiUhkUYKqA341zyOm3mLhhL/JKg4+42+cjFsT9RlxdX9cekPo7tjoxJV6WaGtCkq1R6cbXiEdpvJGIiHNVIsCsXLmSJk2aOBx79NFHMZlM9q9NmzbZz/3xxx8MGjSI+Ph44uPjyc0tfvn+auXszjS9YyE/G80579AyeL4dbFxwxu0JyDwCmYexndXijI9O3J/OXcnt+MY474zja9xFv3jL5vSfm8YbiYg45/EAc+zYMcaNG4dhGPZjOTk5bNmyhblz5zJ37lwWLFhA586dAcjMzOTaa6/l0UcfZd68ebRs2ZLx48d7qvwK8Z2tNVlDPi8cJ/PTawQvuJD+Pl+BzVr04pTvoU5M4WaNZ7A66RCH8/yA0q8EXBb6xVs2+rmJiLjG4wFm2rRpDBkyxOHYokWLiI+PZ8KECUyYMIERI0bYz82fP5/o6GiaN28OwKBBg3jllVf466+/KrPsimcyQctecPe35HadylifDwl6tQskfgTGKb/cdq6BJl1KfFyP2HpEBhaGxDONr1HLSfWkz01EahqPBpiPPvqIzp07c84559iP2Ww25syZw4033ki/fv34/fffHe5Zvnw5rVu3tn8fERFBw4YNWbFiRaXVXanMPljjBnBl3jPkd7oLPptC0KLudDNvxXzoV/jlPWg3qMTHxDUMI+G6+sT7fMiLA9sXO9VbLQDVkz43EalpPBZg0tLSWLFiBUOHDnU4np2dzaxZs5g8eTKbNm2iY8eOrF271n4+KSmJOnXqONwTHh7Ojh07in2tVq1aER0dTXR0NHPmzHHr+6gsVnzJb38HjN9GfmxfnvJbQNDr3aDzcGjY0aVnxDRvzr1+7xEbluf2+tQCICIilcljAWb69Ok89thjRY4HBwczYMAAHn74YX777Tcuuugihg0bRkFB4S/G9PR0IiIiHO6xWCxkZWUV+1rJycns27ePffv2cc8997j3jZyBv6+Z8Ve2wN+3+B9zqVfe9Qsk/4J4OucmcGJSCvR4xOXXwj+YVKMWpvQUl9+DqzWqBaDyKCyKiHgowLzzzjtcfvnl1K9f/4zXBQUFsXjxYg4fPmxvYYmIiCgy6yg7O5vw8PAKq7esLL4+TOxxLhZfH6fnXV151zkT+AW6/FonpRj1MB/b5fKrlK9GqQgKiyIiHgow8+fP55ZbbrFPkb7jjjvYs2cPJpOJ9evXO1xbv359YmNj7aGlWbNmpKWlOVyTmppKTEzJ66BUNcVuAlmBkmyN8Dm43eXrPVGjiIhISTwSYF555RW2bt1q/3rkkUdo0KABW7dupVOnTg7XFhQU4O/vT5s2bQDo37+/w5owR44c4ciRI1x99dWV+h7cocRNIP9Vpg0ei7HdaIr54Da31ygiIlKZPBJgmjdvTrt27exf55xzDv7+/rRr144ffviBGTNmkJmZidVq5YEHHmDmzJn4+BR2jQwdOpSdO3dy6NAhAN58801GjhxJdHR0pb4Hl8aclMCVTSDd3YWzzdYcn4M/Q342UHI4cnWjSk/TuBARkZrF4+vAOLNgwQIaN25Mnz59GDBgAFdeeaX9XHh4OB988AETJkxg0qRJ7N27l7lz51Z6ja6OOSlJSZtAursL5w+jIUZwPdi51uVwVJqNKj1F40JERGoWFzfbqVhDhw61T6fu3r07Bw4cOOP1nTp1YunSpZVQmXuVpdWmR2w93vpxD6nH89zUhWMiv80tWL77D6vPObdIOBp2aeNyPr94udYCEtb9yehuzezBz9kxERGRklTJFhhvVZZWm4rowsnrdBdk7KdHxrIyjW8p65gcZ60kajkREZGyUIBxE3eMiSnO6V045R7UawmFW98nbu8SFoYv5iz+cTkcaVq1iIhUBQowbuKuMTElcVuAiGwJo76mVeNovrFMoP3WafDP3hJv07RqERGpChRgqhm3BoiAMHKveISr8p4EWx7M64RlzYNEUHzLjqZVi4hIVaAAU81URIDYa9Qj57r5MHwt5n/+YoNlAv5fPwE5RVt3qsu0am/hzjWARES8iQJMFVPSWJoKDRD1W5Pd/y2G5N2PT8q38Hxb+O4/9jVjTnJlWrV+8ZbNqT83jTcSESmeAkwV48pYmopclyX5QAZbjHPZ0u1N6LsQtr9L8ILzucXnS7BZXXpGeX7xOluQrqYsUnf6z+2tjSkudxcqMIpITaMAI3aOv0C3khjUGe76itwrZzLCZwXBr1wMv34AxpmnPJdnnE5Nnmp9+s/NBC51F6qlRkRqIgWYGqSk7imnwcNsxtryBq7Ke4q8C8fCF9MJev1KLjdvA8Ow33tqC4AG+pbN6T+3QRec49hdGGGDA9swHdsF/P9nr5lhIlITKcDUICV1T50peBTgQ37b22DsFvJbD2COXwKBb10HKT8UaQEAPD7Qtzp2Ozkb39Qqqhbnmf6kw9d3wlNNYfGN/PXyQIb5fMYfa/8LtgIFRhGpkRRgxM6lAcJ+AeSfP5rLcp+joFEXeLMfq9/5T5EWAE/vn1Rdu50cfm42G/4bZvO2/0wKGnSAiYkkDv6ZO3yf5LWCa7n72yASF9xJXITN44FRRKSyKcCIA1eDxwmCyOsyBcb/TI+zoR5HAYgM9lELgDvYrPDBnfglLeP6vMfIu+wBqBVV2F10ojAsHrLVYk1WC1jcl1YRhbed/NyqYwuUiEhpKMBI+QTXIW7ANObd2pEepk28VvAgcVsexnT87wp7yerwy7k8s6lM2AhYMQ4OJ5M1ZBU7jWj7uSLdRf3vgsBwAj+5CxM1b+CziNRcCjDiFi1btGC10ZkmIxZD9lGCXz6f+32XQPZRt7+Ws1/OVW0asSuzqYoLNBN938dn/yYY8jFG0FkO54p0851TF/q/hvnIn9zts7yC35WISNWhAONhFbkJpLuUpkbbWS3gpv+SNfgTYk17CJnfEb56GnJPVFh91XUasbOQ47tjJXf4fE52v8UQWs/pfUW6+QJrk33jK4z1/Qhz6m8VXreISFVQqt+aqampPPnkk3Tv3p169erh5+eHv78/jRo1om/fvixevJi8vLyKqtUrVdYmkOVRlhpt9dsxJP8Bsvsuhj8+hxfa4ffTAvzJd7jOHS0nXjON+MifBKwYw335d2Gr27JUt9rqt+ONgqsI+GwC2NRtJCLez6UAY7VaeeCBBzj//PP58ccf6dixIw8++CAvvfQSCQkJTJw4kUaNGrF48WJat27N+++/X9F1SzVR0OhSuHM1XPcCOzZ9zmCf1ez46h0osLqt5aQqTCMu97gcWwF8OIr81rfwme2CMj3iOWs/TCcOwi/vla0GEZFqxLekCzIzMxkzZgx9+vRh9uzZJT4wNzeX1157jblz5zJx4kS3FFnTVYdupjMymUgMvYQ7M/1ILchjxVfpvPbrjaxuMILU44X/CZ5sORl2aeNSP/7kuJCbFvxgn0Z8PCe/5Bvd6GR30PAuTcrWmvbt85B9jNzLp8O335Sphhws5F42lcAvH4Nm15bpGSIi1UWJvxEXLlzIf/7zH2644QaXHmixWLj77ru55ppr+OKLL8pdoFRsN1NlhSPH6b9hrKl7Gz1SnqeeT2HXUXlbTjy97kx5mA/9Al89BX1eBr+gcj3LGtcfAmvjt3lhkXPVYfaWiIirSvytNWrUKEJDQ0v94JYtW9K1a9cyFSWVp7LG4BTp5rmyJ3GTVvDSBccY5fMxC8MWEWdNrNAaTlcVNo70J5+AT0fDJeMgumP5H2gywxUP4b9xHgHkOpzS1GoR8SYlBhiLxeLw/eOPP87s2bP5/vvvAbjvvvsICwujY8eO/Prrr2e8V7xDWVptnK7y6xdIi6tG8GZBD1q2jIW3biLw3ZuJM+2uqNIdVIWNI8f7fgBmX7hssvse2qIHRkgDbvJZ775niohUMaXuN3jxxRe59tprueiii3j11Vd59tlnmThxIq+//jrz5s2riBrlNJ4eE1PWVpviunlOEETepffB+J+x1WnJ+/6PEPDhMEj93W01V0Xmv7cyzGcVOb3mgY+f+x5sMpF38QRG+n4KBZoVKCLeqdS/AW+77Tbatm1LVlYWDz30ENdccw0zZszgvPPOIyxMe7BUhuow9bpMgs9iW8tJdMpNINFoBPO7wEej+W3HDqDqLFLnFtZcAlaO4yXr9dgi49z/+HN7k2P445v4gdufLSJSFZQ6wGRlZZGYmMjIkSP5559/eOGFFwA4ePAgS5cudXuBUnOcnFadSRAjdnclsf86Eo8HMWbJVgDil/xUbRapK9GGZ8DkS0LB9RXzfLMPrxZci/9PL4NhVMxriIh4UKkDzOTJk3nuuefIyMhgxYoVNG3alFWrVnHPPfcQF+f+/5OUmqPIgnT7/Vld/y4OEw7A4RNW1nz6DuT848Eqy898cDt89x9yer2AteSVDMrsw4JLMKfvhZTvnZ7XrCQRqc5K/bfngQMHeOWVVxyO9ezZk1q1atG2bVu3FSbVW1nG6fSIrcdbP+4h9Xiew7Tqk8cig0x053tCXnqM0T49Ie9iCAh3+qxTV/g9v+lZTq/xBD+sBKwcB5eMw1avDXDmTS/LM94pmwDy296K/8aXodHFRc6Xe+0aEREPKvXfim+88YbT4y1btqR3797lLki8Q1nG6TibqXTqsXmDLyBu5Otk91lED58tBM/vDD/MB6vjdOGqvDfSRN9/V6nucq/T86cHlvKOd8rrMAx+Xwnp+8t0v4hIVeVSC0xycjJPPvkkKSkp7Nixg99+K7phXEpKCjbtwVLteXqGk7OZSqcfK2h8GX3yHiHpGhtBX8+G7+fhe/G9+BABON8bKa5hGLnWAhLW/cnobs080uLg89cGbvf5nJzrvyTY1x+sRVcLPhlYSuv0z+1kC1RSVgTnN+8OP70KlzxQvjcgIlKFuPRbqlWrVrz88ss0atSIBg0a0LVrV4evyy+/nFGjRrF27dqKrlcqWHWY4VT4y9nEr8GXwKhv4cqHsXz/HF/434dv8kf0aFnX6d5IZV3jpSxjRYpsUpl1lD8/nMktedNIzGtQqtd3xamfW5EWqKbDYMtiKKjc7RVERCqSy2NgLBYLr7/+OmvWrKF79+5Orzl+/LjbCpPqpbJabk7/5bxo6PnEnTeAzGa9WPj4NGZ9OZ244Od5qcsD9F9p/H/RvHI4faxISS05RWoc0hFWTeHOzLs5bATb687MtQKlH6dT0s+6SAtURnPi/ALx3fk54F/Kdy8iUjWV+rdN9+7dsdlsHDx4kJSUFPvX7t27mTFjRgWUKNVBZbXcOOseAsDHj6UFV5J5149w3i10+D6ed/0fpXWB+7cnKKklp0iNK99n9eFaHLYG248t2ZhS5nE6Jf2si2zbEFsfOtyG38+LXX4NEZGqrtQBJiEhgYiICBo2bEiTJk3sX82aNeO5556rgBJF/q/IL+fTN4D0C4SL48kctZlvCtoQ+N5AeLMfHNjqkRojA2x0P7yIHjcOdajbAOdBzA2cbtvQbjA+f22gIaluex0REU8qdYC5//77mTx5Mr/++iu7du1y+Jo4cWJF1Chi5/SXszOWUF4o6MuJu7dAZCy8dg0BH95BS1NKhdR16piXkzU24QCvmR4jbuBjxMW1caj71gvOOXMQK6cig6FrNcDa9Epu8v3Kra8jIuIppQ4wF1xwAXfeeSetWrWiUaNG9q/GjRszfvz4MhWxcuVKmjRpYv8+KyuLESNGMG7cOAYPHsyePXscrt+4cSMDBw5k9OjRTJ06FUMrjdYozmYqFTsuJDACrnoMxm3FCK7LR/7TCXxvEOz90W31OJu23TpnE59aHqTZDfdD8yuL1O1yEHOj/Ha3McBnPdiKH4ysxe1EpLoodYBZvHgxixc770tfsmRJqQs4duwY48aNcwghQ4YMoUuXLrzwwgtMmTKF6667jvz8whkU+/fvZ+DAgSQkJJCQkEBeXh5PPvlkqV9XvEuJY3BqNSD3qqfokvs8BXVbwuK+7Hr1Ds7hIMkHyrdOTJExL1+uJvDDO7g/fwTWltcVe19xm1tWlF+DOvOq9Rp+/6n42YKVvRu3iEhZlTrAdOrUiSlTpuDj41Pka+rUqaUuYNq0aQwZMsT+/ebNm1mxYgU333wzAG3atMFkMtn3WXriiSfo2rUr4eGFK7DeeuutzJo1i8zMzFK/tnheZa87k0pt8i5/iMRbvmfEof6kUJ+xr64m8fPXISej6PRnFziMefHLofueOWT3f4tPbEVXv/WUwlai7bxa0IvRX2RWqcX9RETKotRbCYwdO5bw8HBatGiB2fz/Xzp5eXm8/vrrpXrWRx99ROfOnR2OLV++nObNm2OxWOzH2rZtyyeffMKQIUNYvnw548aNs59r06YNWVlZrF+/nl69epX27UglchZWyrpwW3mt/jObw3l+ABwqqMWardthYw/i8+8DAv4/RduFrp24BqHMv/AfVny5jhvqpRM3+B2OB9QHvqjYN1EKp7YSHbYGs2bbnwy7so2Hq5KK4ulFG0UqQ6n/t/fuu+/m5ptvplu3bg6L2XXq1IknnnjC5eekpaWxYsUKhg4d6nA8KSmJOnXqOBwLDw9nx44dZGZmkpKS4nDe19eXkJAQduzYUexrtWrViujoaKKjo5kzZ47LNYp7VaVF8orMZrrtfla3eYrD+QHAv11BH74KWxZjTvsdfxwXgbOQh/nAFvjqKXihHe23TiODEJre+RrUPrvS309JTn2/9Xwy6W771sMVibs4G7ekrkCpCUrdArNt27YixwzDYMuWLeTn53Pfffe59Jzp06fz2GOPFTmenp5ORESEwzGLxUJWVhbp6YXN3sWdL05ycjK1alXOOAOpHk4Oor1pwQ//H0R7fjveSvrx340joXtUPmxbQtCBe/nNkgv/qQu+FoKtufwecBjj3drQpAv0eITMxlfz/mPreNjkmS0YSnLq+0240o+47a9wvPtdni5L3ECbckpNVeoAc/nllxd7LioqyqUA884773D55ZdTv379IuciIiI4ceKEw7Hs7GzCw8PtwSU3N9fpeZHSOH0Q7am/5OcNvpC4pr2AKZzIyuHKR99n3e0tCPY1yLZCtwU7WTtlAKGB/3Z15lT9ZfpPvs9zL7gaNt+Hz18bPFyRiEjZlfp/F2fNmkVBQQE2m83+ZbVaeeSRR9iyZYtLz5g/fz633HILJpMJk8nEHXfcwZ49ezCZTOzdu5e0tDSH61NTU4mJiSEgIICoqCiH89nZ2Zw4cYKYmJjSvhWRIpzODDL7cJhwbFEd4JwLsUV15DDhUEVbW0rk4wftBvHHd8uB0g1YFhGpKkr9N/A999yDyWRyfIjZzA033MAtt9zi0jNeeeUVtm7dav965JFHaNCgAVu3buWJJ55g+/btDq0sSUlJ9gG6/fv3Z9OmTfZzycnJhIWFcckll5T2rYjUWIkNb2LEH+cDpd/KQESkKih1F9KhQ4eKHMvJyeHtt9/mp59+cukZzZs3d/h+27Zt+Pv7065dOwAuvfRSVq1axQ033MC2bdswm80MGDAAKJwFddVVV5GTk0NAQACLFy9m2rRpBAQElPatSBVUkdOqT50iXZrNE73R6v1+HKawS/bkVgaVsZieVBQD36RlEBUH9Vt7uhiRSlHqANO4ceMiLTCGYRAUFMTs2bPdUtSSJUuYOHEi3377LWlpaaxcuRIfn8LBac2bNychIYHhw4cTERFBZGQk9957r1teVzzP2bRqd4QaZ7tYn3NWkEv3emPw6RFbjyXf/s7hHHOFbGUglety8zYClz8NtRrC+J89XY5IpSh1gImPj2fixIkOa8D4+/sTGRlpDxmlNXToUIfp1HXq1Cl2tV+Anj170rNnzzK9llQ/7lgrxtku1sMubVzifc6CT2auFajegSauYRjzBnfmy9cfonu3m9T6Us319/ma3IvvwfLrO7DrKzi7q6dLEqlwpf5f2gcffJAmTZrQqFEjLBYLoaGhNGjQoMzhRaQinN5qU+Iu1sU4Pfgs2ZhSZN+j6qpV9FmcZTpO231veboUKQ9bAZeZf8baoic07w6713u6IpFKUeoAU69ePZ5//nkaNGhAw4YNqVu3Lm3atLEv9S9SFZy+aF5ZN088PfgYUKQlpzp7u6Abvr+vgKyjni5Fysh8dCd+FGCr1xbOuQhSfvB0SSKVotQBZubMmcycOZPbbruN//73v6xatYrp06ezcuVKXn755YqoUcQtyrJ54unB59YLzilTS05lc3Xc0J9GQwoatIefnf8PiHanrvrMh7bzm3EOmH0gqj0cSgRDK/CK9yv1GJjVq1fzxx9/ULt2bYfjN910E3fffbe76hIpVmVvAHlq8AkN8Cu6gq+LNVZm3aUZN5Tfbgi+PzwPF44uMmBZq7xWfT6HfiHR1ojmAGc1g4I8TOkpni5LpMKV+m/Szp07FwkvJ+3evbu89YiUyNN7KrnSkuOsxtOPVXYQK4415jo4cZDEzRu8ZnxPTWI+nESS0bjwGx8/OKs55tTfPFqTSGUo9d+ceXl5fPHFFxiGAUBBQQFfffUVffr0ITAw0O0FingrTwcxO79A6DCE1d9871Xje2oK87Fd7DZO2Zalbkt8jhS/ua2ItyjTGJg5c+YQHBxM3bp1sVgsXHHFFRw9elRjYESqq/Pvosc/7xIZXBimzjS+R+NiqpCCfEwZ+0ixnfJZRTQhef8RQNtEiHcr9RiYWrVqsWrVKn766Se+/fZbADp16qSl/KVGqyrdQWVW+xziYluzgE3cuKWDfXzPcSebVGpcTBXyTwqJtibs5yz7uKVEoyl3JRWuTH5y7SKt8yPeqEx/2+7cuZPWrVszfvx4hg8fDmDvUhKpiapMd1B5XDiatn8uIIicUs3UEs9J/GMnw/InA2b7uKXVRyM5XBAMqCtQvFupA8z8+fOJiYlhwYIFAAQHB1OnTh0GDBjAvn373F6gSEWp9q0mpVTi+z37AmzhTejns6FyC5MyW/37UQ7bCsPmybDSo83Z1KNwXZ+qPNVfpLxK/Tf3kiVLeO+997jrrrvsx2JiYrjqqqvsrTEi1YGzVhNvDjUlthKZTOR1Gskwn8/A5vr4Fo2J8ZweQbuI9M8F/h9W4s5tzqt+TxNMdqkWbRSpbkr9t3SHDh3o27dvkd2fMzIy7GNiRKorr+gKKgdrqxvxMxXg+/snLt9zckxMnlWLp1W2uIIk5ncu7CKyhxVfC7G184gx7VVXoHi1UgcYs9nMsmXL7N9nZmaSkJDAI488QpcuXdxanIhUMh9/5luvw/+7OWBTIKnyMg7Q6uw6gOO6REataKJNaZ6qSqRSlDrAzJo1i08++YSQkBCioqIICwsjPj6eLl268Prrr1dEjSJSid4r6IopKw12rCrzM9StVEmO/40tpEGRw7aQBtQ1HfNAQSKVp9TTqAMDA3n99deZMWMGv/zyCzabjZYtW3Luua4tWy4iVVsu/uRdEE/AhqehUfcyPUNTrStBgRVOHMIIbQD843DKCKlPfZMmVYh3K7EF5uOPP3Z6vFGjRvTu3Zvrr7/eaXjJzc1l/fr15S5QRCrGmQYs57e7Hf7Zg8+uNR6oTFySeRgMG0ZIvSKnjND61FMLjHi5EgNMVFQUTz31FPn5RRe0Ks7vv//Oww8/rMXtRKqwMw5Y9g+GLvdiWf8YZjQWpkrK+BuC64KPf5FTthAFGPF+JXYhde7cGZPJRLdu3ejcuTOXXXYZMTEx1K9fn6CgIPLz88nIyGDv3r388ssvrF69mrCwMJ577jn8/Pwq4z2IuMSbp0hXiM53kvT1x8SYUuyrvEoVcvxvCC06/gUKu5AiUYAR7+bSGJhOnTqxdu1a5s+fz7Rp00hOTsZkMjlc4+/vz1VXXcWkSZO47LLLKqRYkfI42eJQWjU1+CQeyuHOrHgOG2bGLNnCoju0JH2VcvxvqBXl9JQR2oD6pmNYtUK6eDGXB/FaLBbGjx/P+PHj2b9/P7/++itHjhwhMDCQ+vXr06FDBywWS0XWKuIRpwefmhJoVicd4nBO4XtMPVG4yqsCTBWScaDYFhhbSH2CTXkcz82AwDqVXJhI5Sj1LCSAhg0b0rBhQ3fXIlItlLUlp7rpEVuPt37cQ+rxPOqZjtG9oevj4JzJtRaQsO5PRndrpplJ5XDy5zgu4wA+dZo5v8g/hEzDgvnEQaitACPeybv/F1JEyiyuYRgvDuwAwMutk4nbPB0Mg+QDGQD2f7rKXSv21vQ1Zk7+HI2MAxBa3/lFJhOHjHBMJw5WbnEilUgBRkSKdXJ112a9JsCBbSR+/RFjlm4BsO9+XJGchRVtXVDIdOJvCHU+BgbgkBGB6fjflViRSOVSgBERu2LH9wTUhmufYvWXq0k9ngf8f/fjiqSwUjzz8b+hlvMxMACHqY35xKFKrEikcpU6wPzzzz/FnnvuuefKUYqIeNoZ14Zp3Y8eTQOp55MJ/H/3Y3eq6d1DrgomG1PeiWIH8QKkGmGYFGDEi5U6wIwePbrIsby8PKZPn86UKVPcUpSIVE1xNz3EwloLaG/aYd/9uKxjYpxRi4tr6puOYvgGQGB4sdccNmpjylSAEe9V6gDzxx9/cMstt/Dnn38C8OmnnxIbG8uCBQto1qyYEfEi4h0CwmjWbwZv+c+mtfVXEvenV+qYGClUz3QMI6QBmExOu/38fc20j22JjwKMeLFSB5j169fz6quv8tBDD9G1a1f69+9P37592blzJ59//nlF1CgiVUhB9AXMst5KwEfDWL11h9MxMe5slZGi6nMU278zkJx1+1l8fbjmonYaAyNerdQB5vHHH+eee+7h3Xffxc/Pj759+zJo0CBCQ0M5++yzK6JGEali3irojvXca+mxczZ1Qwq3DDk5JkatMhWvvunYv7tQn0FIfThxCLQar3ipUgeY2bNns2rVKpYsWcKaNWtYvHgxa9euZfr06TzwwAMVUaOIVEG5PZ4krq4vC2u/gQ/59jExq5MOVepMpZqonuloYRfSmYTWg/wsyD1eOUWJVLJSB5gePXrw66+/MmDAAAB8fHyYNGkSAwcOZPny5W4vUESqKB8/GLCIuIA0FvrNoVVk4a7IPWLrUTe08N8rYqaSFLbA2EpqgQmoDT6WwlYYES9U6gDz4osvEhoaWuT4Oeecw/vvv++WokSkmvAPIrv/W4SZMgl8byDkHndYwfdkq4wzmjJddvVc6UIymQpbYRRgxEuVOsD4+/uTkpLi8LVnzx5WrFjBV1995fJzUlNT6d27N6GhoXTq1Imff/7Z4fyjjz6KyWSyf23atMl+7o8//mDQoEHEx8cTHx9Pbm5uad+GiLhLQBiD86YCJvjvDZCZZl/B9+Q/ndGU6bIrnIVUzDYCpwqpD8e1nYB4p1Jv5ti4cWNMJpPTc82aNWPUqFEuPefxxx9n3LhxTJo0iTFjxjB48GB++eUXAHJyctiyZQtz584FIDg4mM6dOwOQmZnJtddey2effUbz5s2ZN28e48ePZ/78+aV9KyLiJlkEkD1gKaErx8CCbpj7LfZ0Sd7LVkBd/iGnpDEwoBYY8WqlDjDjx49n3LhxmM3/b7wpKCjgueeeczm85OTkMGnSJKKjowF45JFHGDFihP38okWLiI+Pp3v37kXunT9/PtHR0TRv3hyAQYMGUbduXe6//34aN25c2rcjIu7ia4EBi2DD0wS9eS3Xmu8ErrKfPnVq9flNz/JMjV7AlJWKr8mGEVKv5ItD6qkFRrxWqbuQpk+fTpMmTWjUqJH9q2nTptx22208/PDDLj0jICDAHl4A8vPzufvuuwGw2WzMmTOHG2+8kX79+vH777873Lt8+XJat25t/z4iIoKGDRuyYsWK0r4VEXE3kwm63kdO75d4wm8hllX3QF6mW6dW1/Q1Zn77ay+P599C8mEXus5PTqUW8UKlDjARERFFjuXl5bFp0yY+++yzUhewb98+3nzzTR588EEAsrOzmTVrFpMnT2bTpk107NiRtWvX2q9PSkqiTp06Ds8IDw9nx44dxb5Gq1atiI6OJjo6mjlz5pS6RpGaqtjNHUtgPfdarsmdjTltByy4nNU//uyWqdXOglBNGgycuD+du1ccZUHB9a4FwVC1wIj3KnWAMZvN+Pj4OHwFBgYSHx/PmDFjSvWsbdu2MWrUKD777DO6detGbm4uwcHBDBgwgIcffpjffvuNiy66iGHDhlFQUPiXU3p6epEQZbFYyMrKKvZ1kpOT2bdvH/v27eOee+4p7VsWqbHOuLljCfZTl+xBH0GbAfTYPolI/8IWg9JMrT69tcXZGjM1aTDw6qRDpGYXjkF0KQiqBUa8WKkDzM0338yaNWv48ssv7V/ffPMN+/bt48knnyzVs9q1a8enn37KsmXL2Lx5Mx988IHD+aCgIBYvXszhw4ftLSwRERFFZh1lZ2cTHl78pmYi4iFmX+h6H3GjFrGw7jvcYf6Mly63FTu1+lTOWltq+hozPWLrEelfGOBcev9qgREvVuoAM2vWLLp160bXrl3tXxdddBFRUVEcPXq0TEX06dOHrl27cuDAgSLn6tevT2xsrD20NGvWjLS0NIdrUlNTiYmJKdNri0glqBdHkzvfwGryoeNXw+C9OzCl7yty2aktLs5aW1xdY8ZbxTUMY0GLjVxo+tW19x9SH3L+ITcns8Z0s0nNUWKAOX3NF19f3yLHUlJS2L17NzNnzixzIbVr1+a8884rcrygoAB/f3/atGkDQP/+/R3WhDly5AhHjhzh6quvLvNri4jryjouBrMPiwuuInPE9+AbQPDCi5ng+37hcvcUbXFpWjfYaWuLK2vMeLPWtt9paj7k2vsPrgMmM/n/HKwx3WxSc5T4N1D37t1p0qRJiV/NmjXj+eefd+lFMzIyeOONN8jIKPy/re3bt5OXl8dVV13FmjVrmDFjBpmZmVitVh544AFmzpyJj09hH/zQoUPZuXMnhw4V9uu++eabjBw50mFWk4hUnPKMiwEKV5Dt8xJZAz+iq3k7wQsuhO3vsTrpoEOLy+60rDK3tlTkwF5PDxo2nTjIIaO2axebfSA4ElOmxsGI9ylxHZjx48czduxYdu7caQ8RzlitVl588UWXXvTQoUPMmDGDyZMnc/nll9OsWTPefvtt+/kFCxbw4osvcuGFF/LQQw/ZF7GDwhlHH3zwARMmTCAqKgqTyWRf8E5Eqg9bVAf65s0gqWsGgasfood/a5YEDeVw1v9bXM45KwgofWvLyYG9w7s0KVXQyrUWkLDuT0Z3a1bsfWV9truYj//NQaPobNBihdbDfOIw4FdhNYl4QokB5tlnn2XSpEk0bdoUKFyzxc/P+R8EV7uQWrRowe7du52e6969u9OxMKfq1KkTS5cudem1RKTqMjBjbX0TnNeHuB9f5rV1M1lo7srgq/sQ1zCM4zn5lVqPp8NJifJzMOUc47BRikkLIfX/bYFRK7V4lxK7kC699FKefvpp+/dnmir9ww8/uKcqEak2yjwmxuEhQXDpRBrFf0SsOYWOq26AZSMxpe91S42e7vZxm+N/Y5h8SKMULVIhkZg0lVq8UIl/48TExPDf//6XdevWsWHDBv7++2++/vprNmzY4PC1atUqdeWI1EBlHRPjdEXdwAget95K5ojvAAh+5SIe9H0Tsss2w/EkZ2vFuBJqqlzwOX4QI6QeRmkmkIbWx3RCU6nF+5TYhXTfffcxffp0pk2bxoEDBzAMo9hl+4vb5FFEapaSWmVOn3G0aOj5DoN0jbCzoe/LZHUcRYuF4wmZ3wkuGQ8X3g34u6VGV7qLqlyXUsZ+jNAoSC3FPSH1MO/fWmEliXhKiTHez8+PJ554gpSUFKxWK8OGDcNmsxX5ys7OZty4cZVRs4hUcSW1yjhb48UZW73WDM2fQnbfN+D3z+D5dvhtXogf1gqrvUr7JwVbWCnHsoTWVxeSeKVSd1qPHDnS6XGLxcIDDzxQ7oJExPuVdkXdgkZdYPgauO45/LYu4kv/Sfj+8g7YqkjXTmVJ34ut1tmluydEAUa8U6kDTKdOnYo9FxlZs5b1FpGyKdOKuiYTtOxF1rCvmGPtj+WbJ+GlS/DdsRIwKrbgquKfvYXda6URWg9TVhpmtIideJdyTBsQESm7Mq+oa/bhQ1sXMu/6AToNw/L5vXzo/zA+e76ugCrLpsIG/6bvw1arlF1IIfUwGTbOooSdq0WqGQUYEamefPzhgrv46YZ1zMofyK53H4T/3gj7t3i6sorZIdswIL0MLTC+FoyAcCJN/7ivFpEqQAFGRKqtxP3pjHkvmZ+MVtxheoTE4AthUS8CPhxGM9N+T5fnXtnHIO9E6VtgAFtIPQUY8ToKMCJSKdyy4N1pHGYznchnTVg/GLsFIyiCFf5TsawcD/+4ZzE8j0vfCwFh+AeHlfrn+KtvKw4bYY5r7ohUcwowIlIpyrsJ5EmnBiGns5lqNSD36me4Ou9JTNYcmNcJy5oHiaBq/fIu9TiZ9H0Qdk6pf46J+9MZsb83iUZTxizdQuJ+jYUR76AAIyJVlrNWm1N/gRc3myn5QAZ7jPps7fQ0DF+L+dhuNlgm4P/1E5BTviDjdAXhMij1OJlje6B2Kce/UNhKdTg/EDjzmjsi1Y0CjIhUWa60Npw+m+n0VX4TC84me8ASbs+bgs+eb+D5tvDdPLDmFHlWSeGkyLP3p1fedgNH/4SzmpX6th6x9Yi0FC7858qaOyLVhQKMiFQZ7hgnU9wqv5uNGLJv/QT6LoDtbxP88vnc7LMObIW/3J2FE1eeXSEzjpw5shPOal7q2+IahvHSZXncbF7r+po7ItWAAoyIVBnuGCfjbFyMvWXl7+PQogfctYHcKx5hpM8nBC+8BBI/ZHXSwRK3NyjtCsJudWQXRJS+BQag1dl1Gev3cenX3BGpwhRgRMSrnD4uBijasmI280vYlVyZ9zTbzh0Pq6bSI3kqkYVDRYoGn3//WaYVhN0hP6dwFlIZWmCgcBp1Xf4pXEtGxEsowIiI1zl1XIyzbp+T3UUGPozcHEVi//XEdezKaz6P09+8npe6WwAnwYdyrCBcHsd2g18QhNYv0+1GcD0sJivk/OPeukQ8SAFGRLyas26fIqHmjwy4ZByN4pdztjmVjmtvYfW7L7q0Y3ZFOjlAOP/wDohoWrgfVFlYQsk0LJi1qaN4EQUYEfFqzrp9ih3LElCLudYBZI78iR7RVupxFIDIYJ9yjXcp60ylkwOEbQcTIbJVmV8f4LBRG1OmAox4DwUYEfGIiliZtzind/uUNJbFCK5L3ICHmHdrB3qYNvGabRpx2x6FjL/L9PrlnalkPpwI9VuX6d6TDhOOSS0w4kUUYETEI9y1Mm9ZuTKWpWWLc1ltdKbJnW9AZir8pwOWdTOozfEzPttdi92d5HM4Eeq3KdczUo0wzCcOuqUekapAAUZEpAS2OufCTf+FOz7DnPpb4aq+3zzldFVfV9aTcaa4bqYQsjD/sxvqlS/AHDbUAiPeRQFGRMRVUe3Ivult7sybjM9fG+D5tvhtnIeFPPslxS2kV5LiuplamVKwhdSDkLrlKn2/UQdTupdsbCmCAoyISKltMloWrurb7xX8kpbxlWUifltfB2uey4vdudrN1Nn8OwUNLyh3zXuMepj/+avczxGpKhRgRETKwmSC5t3JGrqWGfm34/fTKzCvE3Fpn/HiwLZA8YvdlaabqZP5dwqi3RVg9mgxO/EaCjAiUq1V5mwmZ5L/Ps4q2/ls7rUCLn8A1s2i48reXG3+kVYNQp3e43I3k81KJ/MOCs6+sNx17jXqYsrPLByMLOIFFGBEpFrz5Gwmx5aUn0msey3Ebya/w5086reIoDd6wM41RVo9XO1m8tn7A7n4YSvnAF6AHCyFY2mO7i73s0SqAgUYEZEyctqS4utPfodhdM2dS36rG+CD4QQuuZ5Opt/s97m6p5LvH5+xuqADmMr3V/XJVirCmxRuSyDiBRRgRETK6EwtKTlYyL9gLIz/mYJzLuV1/6cJfPcWOLANcGEdGmsuvonv86ntonLXebKVyhzRFI79Ve7niVQFCjAi4nUqa1yMSy0pAWHkdZnCZblzC9eTea0nvHs75iN/nPnhv7wHgbX5zhbnvoIjmsDRXe57nogHKcCIiNcp67gYV4LP6de4ujv1MWqRe8WjMG4LBEUQ9Ho3nvadj+mflKIXZx2FtY+Re+kUoIwbODpT51xI/a3k60SqAQUYEZF/uRJ8yj1ouFYU9J5L5vBvMJtsBC+8GN6/E5+dn1OPo5j3/wSL+8DZnfml9pWA+7YkIDIWUn8HW+k2lRSpihRgRMTreXyqtZNF64zajZmUP5qsO9ZBWEMC1jzIxoB4gt7pD40vJfHCOYxZuhUo3ZYEZxTRFAybxsGIV/BYgElNTaV3796EhobSqVMnfv75Z/u5rKwsRowYwbhx4xg8eDB79uxxuHfjxo0MHDiQ0aNHM3XqVAwtzCQiZ+Cs1aQs3UVlUdKidbazWkCPR8kc9RMtcv7LiYm74epZrN7xT5m2JDgjH1+oEwOHk8r/LBEP81iAefzxxxk3bhzLly8nKyuLwYMH288NGTKELl268MILLzBlyhSuu+468vPzAdi/fz8DBw4kISGBhIQE8vLyePLJJz31NkSkmqqU7iJcX7Qu+UAG+fiS/HfhTteurhVTavVi4XCye54l4kEeCTA5OTlMmjSJq666im7duvHII4+wd2/hJmObN29mxYoV3HzzzQC0adMGk8nE0qVLAXjiiSfo2rUr4eHhANx6663MmjWLzMxMT7wVEZEzciWIOGulcXWtmFKLjIWD293zLBEP8kiACQgIIDo62v59fn4+d999NwDLly+nefPmWCwW+/m2bdvyySef2M+3bt3afq5NmzZkZWWxfv36yileRKQUXAkixbXSuDrDqVSiO8G+n7QnklR7Hh/Eu2/fPt58800efPBBAJKSkqhTp47DNeHh4ezYsYPMzExSUlIczvv6+hISEsKOHTuKfY1WrVoRHR1NdHQ0c+bMqZg3IiJSjJKCSIV1FzkT1R5OHIb0vRX3GiKVwKMBZtu2bYwaNYrPPvuMbt26kZubS3p6OhEREQ7XWSwWsrKySE8vHPxW3PniJCcns2/fPvbt28c999zj/jciInKK0g7+rbDuIqfFBUP9NuTv+YG5q3eQa9WUaqmePBpg2rVrx6effsqyZcvYvHkzH3zwAREREeTm5jpcl52dTXh4uD24FHdeRKQqKMvg3wrpLirO2Rdg7PmB59f+QZ7VBkCutaDEQOPKNSKVxeNdSAB9+vSha9euHDhwgGbNmpGWluZwPjU1lZiYGAICAoiKinI4n52dzYkTJ4iJiansskVEqqdmV+D75xrg/+Ng8qw2h0DjjCvXKORIZakSAQagdu3anHfeefTv35/t27c7tLIkJSXRq1cvAPr378+mTZvs55KTkwkLC+OSSy6p9JpFRMrKo4vrNbkMU+YhmpkOuP3RroQcEXfwSIDJyMjgjTfeICOjcFXK7du3k5eXx1VXXUX79u259NJLWbVqFVA4TsZsNjNgwAAAxo4dy7p168jJyQFg8eLFTJs2jYCAAE+8FRGRMnHHGjNl5h9EQaMuXG3+qfJfW8RNfD3xoocOHWLGjBlMnjyZyy+/nGbNmvH222/bzy9ZsoSJEyfy7bffkpaWxsqVK/HxKfxD3rx5cxISEhg+fDgRERFERkZy7733euJtiIh4fJuCssqP60+/P2accTp1rrWAhHV/Mrpbs2KDlivXiFQEjwSYFi1asHv37mLP16lTh8WLFxd7vmfPnvTs2bMiShMRKZWTLSnVjfXcXoSbJuGT8i3EdHN6zcnuoOFdmhQbTly5RqQiVK//ZRAREffwDeDNgu74f/u0FrWTakkBRkSkhnrF2hufw0nw24pKf23NVpLyUoAREamhjhNETo/H4ZNxmDL2V+pra7aSlJcCjIhIFVUZA4StcQMg9gaC3u5LJMcq7HVE3E0BRkSkgpU1iFTaVOtrn8EafSGfWh4sHNQrUg0owIiIVDCPrvniCrMPudc8xwvWPgS+NxA+GQ/Zao2Rqk0BRkSkBko+kOHwT0wm3izoQeadGyB9HzzfDr+NL2Ihz4NVOtLAXzmVAoyISBVQmQviJe5PZ8zSLQCMWbqFxP3p9iCTlBUBt74P/V/FL/E91lruxffXd8FmKxp6KpkG/sqpFGBERKqAyuxmWp10iNTjhS0rqcfzWLIxxTHQHMiA5t3JuuNLns0fgGXD4yS+0I8xi3/4/zX7050+29MhR2oOBRgRkRqmR2w96ob6A1A31B8DHALNmuTDhReazHxo60LmXRtZXbs/qdkUveYUzlp2XKXuISktBRgRkRomrmEYLw7sAMCLAztw6wXnOASa7q0iHW/wtdCjR2/qhvgBUM90jO67n4K/f3ZocTm9ZcdZyCmOuoektBRgRERqoFZRtez/PD3QxDUMAxy7g+IahvHioI4A/Of2S4hrVI/EBcOJf+1LoLDFpWnd4DMHIRE3UoARERGHQAPOu4Ps1zSOhqtmsrrjSxy2BgGFLS679/3tNAhVJnVF1RwKMCIiUoQr3UE92rewt7hE+mbRffPdtPvtGWpz3B52Kpu6omoOBRgRESni9IG+zrqDTu16mjfsCuLuWog5bQcbLBPw/+ZpyD0OOG8VqY6zldS6U7UowIiISBHFjYs5nUPXU4O2ZN/0NnfmTcbnr/XwfDv4PoG8nCyHVhFXZytVtcCg1p2qRQFGREScOn1cjKs2GS3JvvVTuDEBtr1F8IILGOCzHmxWwLXuKVBgkDNTgBEREfczmeDcq2Hk1+Re/hBjfD4m6NUukPgRPVpFaraSlJsCjIiIVByzGWtsP7rnPU1+p5Gw6n7iVtzAS10LAMNjs5Wk+lOAERGRCmfFl/z2Q2HcVmjdl47fjGSJ3yzijB0uP8PZwN+qNk5GKo8CjIiIVB6/QLhkPCdGbeEn41yC3u4HSwdhTk0+423FDfzVOJmaSwFGREQqX0At5lhvInPUJgiLZvdrw7jK/CO/7fjd6eXl2aZAvJMCjIiIeIwRHEli2we5wzyTL2znE79kK4nvzIDjhxyuc2VdmuK40s2krqjqRwFGRKQG8vc1M/7KFvj7uv5roCz3uGJ10iFSMwuDwyEiWHPAH15oj/9XM6lFJuD6ujTOuNLNpK6o6kcBRkSkBrL4+jCxx7lYfH0A18LJ6fe4S5HWlZvGwO2f4HNgC19bxuP//fOQl1XmdWlOp9YW76AAIyIibgsnZWmlcdq6Et2R7IHLuDt/Ar47VsAL7fDb8hp+WM/4LFe2KHC1tUVBp2pTgBEREbdxFoRcCRXFta58Z2tN1pDPodez+G15lTX+9+L767tgKxoqXN2iwFUV1a2kYOQeCjAiIlJmJbW4lCdU2IPP38eh1XVkDdvAC9a+WDbMhvmXwm8rwTDs11f0TCV3bUCp8TbuoQAjIiJOuWNcTFlDhdPgY/bhA9tlZN71A3QcCp+MI+jNa7nQnASUb6ZSmeoRj1KAERERp9wxLqasoeKMwcfXAheMhHHbsDbrwQK/OQS+M4A4064yz1QqVz3iEQowIiJSYco6/dml4GMJIe/ie+iS+xwFka3h9Wtp/+NEmpoOlHumUpnqkUqlACMiIhXK2QDdkrqnXA0+yQcySCeEn2MmwtgtGIHhrPR/AMvK8ZC+z23voTzr0EjF8FiA2bNnD1dffTWhoaG0b9+edevWOZx/9NFHMZlM9q9NmzbZz/3xxx8MGjSI+Ph44uPjyc3NrezyRUSkHFzpnipp3Zci41KOB5F79TNcnfckJms2/KcjrJqKKSvNLTU7q0czijzHIwHGMAxGjBhB9+7dSUhIwDAMevfuzZ9//glATk4OW7ZsYe7cucydO5cFCxbQuXNnADIzM7n22mt59NFHmTdvHi1btmT8+PGeeBsiIuJBxY1L2WPUJ+f6BTB8DRzZSfD8TkzwfR9yj7u9Bs0o8hyPBJhff/2VKVOmMHnyZG677Ta++OILTCYTn376KQCLFi0iPj6eCRMmMGHCBEaMGGG/d/78+URHR9O8eXMABg0axCuvvMJff/3libciIiIeUuK4lPpt4NZ32XzFW+yz1eGvedfDd/PAmlPkWe6aIu2MWmkqhkcCTExMDFdeeaX9+8jISGJjY7FYLNhsNubMmcONN95Iv379+P13x51Jly9fTuvWre3fR0RE0LBhQ1asWFFp9YuISOU7fdyMK+NSEvenc/fafN63Xc6d1gdI3PQlwS+fz80+68BmtV9TkVOk1UpTMTwSYPz9/YscO3bsGL179yY7O5tZs2YxefJkNm3aRMeOHVm7dq39uqSkJOrUqeNwb3h4ODt27Cj29Vq1akV0dDTR0dHMmTPHfW9EREQqjbNxMyWNkzm1m+lwNqyJnc3W9o8SYOTwV0I/+HUZq5MOum2KtFpbKk+VmIW0YcMG+vbtS3R0NMHBwQwYMICHH36Y3377jYsuuohhw4ZRUFD4H0N6ejoREREO91ssFrKysop9fnJyMvv27WPfvn3cc889FfpeRESk6ji9m6lJ3RBGfV+bN2zXcOeJUSR++iI9kh4gMhD7NSe7osrSreRKa0tFdlfVJB4PMPn5+SxdupRHH320yLmgoCAWL17M4cOH7S0sERERRWYdZWdnEx4eXin1ioiI+5VlE0hXnN7NtCs18/8tMjk+rOnwH+I6d+NV39kMMK/jpSv9iWsYVmHdSlrR1308HmCeffZZpk6disVicXq+fv36xMbG2kNLs2bNSEtznBKXmppKTExMhdcqIiIVw127YTtzajdTkYG/cQ3h4rE0jl9OtPkIHb8cCEtuZvWmXypk5d3iZk6p66n0PBpgFi5cSM+ePTn77LMByM3NtXcVnVRQUIC/vz9t2rQBoH///g5rwhw5coQjR45w9dVXV17hIiJSLRU78NcSylxrfzJH/gQRTenx8wQi/QqHJrhz5d3iZk5poG/peSzAJCQksGvXLg4ePMiqVatYtmwZw4cPZ+3atcyYMYPMzEysVisPPPAAM2fOxMenMJUPHTqUnTt3cujQIQDefPNNRo4cSXR0tKfeioiIVBGudEWdaeCvEVwXes4mbvyHLDh3M2N9lvFKk6+ICy1+nGVpaEVf9/H1xIu+9tprjBkzBoDZs2fbj48dOxaz2cyCBQt48cUXufDCC3nooYfsi9hB4YyjDz74gAkTJhAVFYXJZGLu3LmV/h5ERMQ1FTW+xZmTXVHlVvtsmg14jEmPvM4409fwnw74d7iTMNqU+9ElzZwS13gkwAwbNoxhw4YVe/7AgQNnvL9Tp04sXbrU3WWJiEgFcFuo8IA/jYbk9HkNvyO/4LfmETYGv47vDylw8WjwD/Z0eTWaxwfxioiIOFOZLTclatgB8+0fEzB4Kb5/rILn28HGBSTvOwJoSrQnVIH/KkRERIqqyJlJZdbkssI9lq57jsTvPiX+tS+BM0+J1rovFUMBRkREpDRMJmjZi9Wtn+awrXAQburxPNZs+BoMw+FSrftScRRgRESk2vJkN1OP2Pr2KdGRFivddz9D0OKeXGROtF9T3LovUn4KMCIiUm15spvp1CnR826/lLh7V7A98kbamf5g16IRsH9zyTtml1NNXgBPAUZERKSMTp0SnZhq5a7EVrxUcCMjDlxP4sK7iftmDC9dexZQMeu+1OQF8BRgRETEq7nSzeSOriiHna/zLKzpPB9C6tFxxTU87Tuf2GCNf3EnBRgREfFqrnQzuaMrqkh3Udtm0OtZMu/8Gl9TAcGvXASfTcGUmVrkXldmKtXk7iJnFGBERKRGqaydr092FxnhTZiYP4asIV/APykEv9yJe3zfhZzCsOLqTKWa3F3kjAKMiIjUKK62tpQl6JxpmwBbZCwMXEr+oA/oE5FCyMsd4dvnWf3rfs1UKgMFGBEREScqaoaTf+OLOHvil5j6LYRfP6DH1ngiA6xAxcxUOp23dEUpwIiIiFQ2kwmad4cR64m7biyvhrzMUPNnzL/wGHENQsscMly5z1u6ohRgREREPMVshrg+NB79AVmmANpvexhevoyC3z/n+bU7Sh0yvCWcuEIBRkRExNPMvrxb0I3MkT9Cu4EEfBrPu/6P4rP3B09XVmUpwIiIiFQVvgFw0RgyR/3EtwWtCXzvFnhrAPy93dOVVTkKMCIiIlWNJZTnC/qROWoz1DkXXr2KgI9H0Nj0t8Nlntzp2tODgRVgREREyqiiN5M0gs6Cq2fB2J8w/INZ5X8/llX3QMaBCt3pujoMBlaAERERKSNXplq7JeSERZN7zXNcmzcbU046vNCB1Z8ucWn9mLK0lHg6nLhCAUZERKQCuWs9meQDGewyoth6/lwYtooeps3UMx0DoG5I8evHVIcwUhYKMCIiIm5SUV1KRbqLjCbE3bWQF69vwM3mL3nN/Bhxe5eANdetr1uVKcCIiIi4SUWt3nvqTtendhfFtL+Md2zdaHb9fbBlMcELLmCAz3qwWd36+lWRAoyIiEgVV2Sna4fuIhMFza+GUd+Q2/VB4n0+IujVLpD0MRhGkWd5cuaSOynAiIiIeNjpXU+nh4zidrp2YDbzS/jVdMt7lp+b3Q0rJ8Mr3fD56yv7JRU5c6myKcCIiIh42KldT8WFjDPtdA3/Dyc2fBi59RwS+2+A2BsI/OhO3vKbhfnA5mK7olyZqeTpdV9OpwAjIiJShRQXMkp9387jcOlETozazFajOUFL+9Ij5TkigwrH55zaFeXKTKWqNptJAUZERKQKOfN4lzLcFxDGM9abyRz5I3ENgnnNmMbV5h9JuC7KeVdUNaEAIyIiUoW4NN6lDPcZIfXg2qdpPPJtrvLZTKflVxSOkzl+qEx1enowsAKMiIhIFVPSeJfy3GfUbsSk/LvJun0NZByAF9rj/9UsapFpv6akcFIVBgMrwIiIiHix4sKIrW5LuOUtuH05Pvs3scEyAf8fXiDxr4MlhpOyjtNxJwUYERERL+VSS0l0J7IHfkh8/jh8f1vO6jdmlhhOyjpOx50UYERERKqpkrYucLmlxGTiG1sbsm5fTY8relDPXNhaExlgo3tMnSKXl3WcjjspwIiIiFRTJW1dUOqWEpOJuC438J/hVxBr2s2rgc8T90kv+P2zIqv6lnWcjrsowIiIiHipsraUtIqKIMloQuMxH0GH2+DjeILe7MUFpuQKrLZ0PBZg9uzZw9VXX01oaCjt27dn3bp19nNZWVmMGDGCcePGMXjwYPbs2eNw78aNGxk4cCCjR49m6tSpGE72ehAREfEmZd3pulwtJb4WuPBuGL8Na9MreMX/GQLfuQkObCv9s9zMIwHGMAxGjBhB9+7dSUhIwDAMevfuzZ9//gnAkCFD6NKlCy+88AJTpkzhuuuuIz8/H4D9+/czcOBAEhISSEhIIC8vjyeffNITb0NERKTSVNRO1+DCmi6WUH5uchedcxP4NbAjvNaTgI/upKnpgNtrcZVHAsyvv/7KlClTmDx5MrfddhtffPEFJpOJTz/9lM2bN7NixQpuvvlmANq0aYPJZGLp0qUAPPHEE3Tt2pXw8HAAbr31VmbNmkVmZmaxryciIiLOuTJT6eQ1uVgY8ccFJN70NUZAGB/5PwQ5NWghu5iYGK688kr795GRkcTGxmKxWFi+fDnNmzfHYrHYz7dt25ZPPvkEgOXLl9O6dWv7uTZt2pCVlcX69esrrX4REZHq6vSuKFdmKhW5JgVye87h0tznIKAGDeL19/cvcuzYsWP07t2bpKQk6tRxnLIVHh7Ojh07yMzMJCUlxeG8r68vISEh7Nixo9jXa9WqFdHR0URHRzNnzhz3vREREZEKUNbxLq7cd3pXlCszlYq7JoOQUtXnTlViFtKGDRvo27cv0dHRpKenExER4XDeYrGQlZVFenphs1Zx54uTnJzMvn372LdvH/fcc4/734CIiIgblXW8S1nuc2WmUlVY9+V0vp4uID8/n6VLl/Lcc88BheHkxIkTDtdkZ2cTHh5uDy65ublOz4uIiEjpuTJTydPrvpzO4y0wzz77LFOnTrWPeWnWrBlpaWkO16SmphITE0NAQABRUVEO57Ozszlx4gQxMTGVWreIiIh4jkcDzMKFC+nZsydnn302UNiy0rdvX7Zv3+7QypKUlESvXr0A6N+/P5s2bbKfS05OJiwsjEsuuaRyixcREanByjpOx108FmASEhLYtWsXBw8eZNWqVSxbtozhw4fTvn17Lr30UlatWgXAtm3bMJvNDBgwAICxY8eybt06cnJyAFi8eDHTpk0jICDAU29FRESkyipL0CjLYODK5pExMK+99hpjxowBYPbs2fbjY8eOxWw2s2TJEiZOnMi3335LWloaK1euxMen8AfUvHlzEhISGD58OBEREURGRnLvvfd64m2IiIhUeSeDRkXfU9k8EmCGDRvGsGHDij1fp04dFi9eXOz5nj170rNnz4ooTURERKoBjw/iFRERESktBRgRERGpdhRgREREpNpRgBEREanhPD0luiw8vhKviIiIeJYrs46qWshRgBEREZESVbWp1VUjRomIiIiUggKMiIiIVDsKMCIiIlLtKMCIiIhItaMAIyIiItWOAoyIiIhUOwowIiIiUu0owIiIiEi1owAjIiIi1Y4CjIiIiFQ7CjAiIiJS7SjACHPmzPF0CXIKfR5Viz6PqkWfR9Xiyc/DZBiG4bFXr2AZGRmEhYWRnp5OrVq1PF1OlRUdHc2+ffs8XYb8S59H1aLPo2rR51G1VNTn4crvb7XAiIiISLXj6+kCKtLJxqWMjAwPV1K1GYahn1EVos+jatHnUbXo86haKurzOPnMM3USeXUX0r59+zj77LM9XYaIiIiUwd69e4mOjnZ6zqsDjM1m48CBA4SGhmIymTxdjoiIiLjAMAyOHz9OVFQUZrPz0S5eHWBERETEO2kQr4iIiFQ7CjAiIiJS7SjAiIiISLXj1dOo5f/27NnDXXfdxXfffUfz5s2ZM2cO3bp1AyArK4vx48cTGBjI0aNHmTVrFo0aNfJwxTXHypUrGTNmDLt37wb0eXhaeno6b7zxBhERETRs2JBu3brpM/GAEydOMGXKFGrVqoXVaiU3N5enn34ai8UCQEFBAffddx/Z2dn8888/TJkyhbZt23q4au+TlpbGk08+SWZmJgkJCfbjJf2Z2LhxI8899xzh4eHUrl2bWbNmuX8yjSFez2azGT169DCeeuop47///a/Rtm1bIygoyNi5c6dhGIbRr18/44033jAMwzC2b99utGnTxsjLy/NkyTXG0aNHjWbNmhmNGjWyH9Pn4Tk7d+40brnlFiMtLc3huD6Tyjd69Gjj2WeftX8/ZswYY9q0afbvJ02aZDzyyCOGYRjGoUOHjKZNmxrHjh2r7DK9Wl5envHhhx8aLVq0MG6//XaHc2f6M7Fv3z6jSZMmxtGjRw3DKPysZs+e7fb6FGBqgO3btxtr1qyxf3/o0CEjODjYeO6554yffvrJCAgIMHJycuznzzvvPPt/mFKxRo8ebTzyyCP2AKPPw3OOHj1qXHzxxcbBgwcdjusz8YzY2Fjj/ffft38/f/584/rrrzcMwzD+/vtvw8/Pz/jrr7/s56+//np7oBH3GjRokEOAKenPRHx8vDF06FD7uS1bthghISHGiRMn3FqXxsDUADExMVx55ZX27yMjI4mNjcVisbB8+XKaN29ub5YFaNu2LZ988oknSq1RPvroIzp37sw555xjP6bPw3MefvhhmjZtyuzZs7nwwgu57777sFqt+kw85LzzziMhIQGbzQbATz/9xKhRowD47LPPsFgsDl0W+kwqjp+fn8P3Jf2ZWL58Oa1bt7afa9OmDVlZWaxfv96tdSnA1AD+/v5Fjh07dozevXuTlJREnTp1HM6Fh4ezY8eOyiqvRkpLS2PFihUMHTrU4bg+D8/Iyspi4cKFdOjQgblz57Jo0SJeeuklpkyZos/EQ+bMmcOuXbsYMmQI3333HTfeeCPXXHMNoD8nnnamn39mZiYpKSkO5319fQkJCXH756MAUwNt2LCBvn37Eh0dTXp6OhEREQ7nLRYLWVlZHqquZpg+fTqPPfZYkeP6PDxj48aNZGdnc+utt2IymWjZsiU333wzL730kj4TD2nQoAHvvPMOYWFhXHvttfaWGNCfE087088/PT0doFI+HwWYGiY/P5+lS5fy6KOPAoX/keXm5jpck52dTXh4uCfKqxHeeecdLr/8curXr1/knD4Pzzh48CAAwcHB9mOXX3452dnZFBQU6DPxgOTkZL7//ntefPFFnnnmGfr27cvnn38O6M+Jp53p538yuFTG56MAU8M8++yzTJ061d532axZM9LS0hyuSU1NJSYmxhPl1Qjz58/nlltuwWQyYTKZuOOOO9izZw8mk4m9e/fq8/CA0NBQoLBr9aTIyEgA2rVrp8/EA4YPH06XLl3s/x4fH8+9994L6O8tTzvTzz8gIICoqCiH89nZ2Zw4ccLtn48CTA2ycOFCevbsad+hOzc3l759+7J9+3aHtJyUlESvXr08VabXe+WVV9i6dav965FHHqFBgwZs3bqVJ554Qp+HB1x00UX4+fnx448/2o8dPXqU9u3bM3jwYH0mHvDzzz87jN+788472bt3LwA33HADR48etX8P+kwqU//+/c/4Z6J///5s2rTJfi45OZmwsDAuueQSt9ahAFNDJCQksGvXLg4ePMiqVatYtmwZw4cPp3379lx66aWsWrUKgG3btmE2mxkwYICHK/ZezZs3p127dvavc845B39/f9q1a0eXLl30eXjAWWedxfjx43nppZfsxz788EMefPBB/RnxkGuuuYbPPvvM/v3vv/9O7969gcLWsaFDh7Js2TKgsAswKSmJ0aNHe6RWb1dQUOAwBqmkPxNjx45l3bp15OTkALB48WKmTZtGQECAW+vSbtQ1wGuvvcadd95Z5PjYsWN54YUXSEtLY+LEiTRo0IC0tDRmzpxJVFSUByqtmRYtWsSMGTP466+/APR5eIjVauW+++4jKyuLWrVq0ahRI8aMGQPoM/GE9PR07r33XqKioggJCeHgwYM89NBDhIWFAYXdEuPHj6dWrVqkpaVx3333ERsb6+Gqvc8777zDfffdh9ls5qmnnrKHlJL+TKxatYo333yTiIgIIiMjmTZtmttrU4ARERGRakddSCIiIlLtKMCIiIhItaMAIyIiItWOAoyIiIhUOwowIiIiUu0owIiIiEi1owAjIiIi1Y4CjIiIiFQ7CjAi4nV27NjB66+/7rHXnzdvnn2HaxGpGAowIuJVdu/ezYIFC7jjjjscjq9evZo+ffowdOhQRo4cyYgRI5g5cyYzZswo8Zk7duyw7yDetGlTPv/8c4fzX3zxBVFRUVx44YUkJSUxcuRIHnjgAYfdrUXEvbSVgIh4DcMw6N69O0uXLiUyMtJ+/KmnnuK1115j1apVNG7c2H581apVTJ8+3WHn3OJYrVaaNm1Ky5Yt+eKLL4qcv+GGG3jppZfs+8Fs3ryZZ555hqVLl5b/jYlIEWqBERGv8fHHHxMSEuIQXr788kvuv/9+3njjDYfwAtCzZ0+GDx/u0rN9fX0ZOXIk69ev59ChQw7njhw5QkREhMNmdh07dmTXrl1s27atzO9HRIqnACMiRRw7doxJkyZx++23c95557Fw4UIAvv76a66//noeffRRXn75ZerWrUuHDh1IS0uz37t9+3bGjBlD//796dSpE9999x0FBQUsXbqUTp068eGHH3LBBRdwwQUXALB+/XpGjRrFs88+y6xZs+zP+emnn4iOjsbHx4dly5YBkJmZSb9+/ejbty/Z2dlF6l6wYAGXXXaZw7GpU6fSqVMn++udbuTIkQ7fz58/n/j4eC6++GIGDx5MRkaG/dzw4cMxmUy88sorDvf897//5bbbbivy7EsuuYSXXnrJ6euKSDkZIiKnGTJkiHHw4EHDMAxj7dq1hslkMjZt2mTs27fPaNq0qdG1a1dj5cqVRkpKihEVFWU88cQThmEYRlZWltGvXz+joKDAMAzDmD59uhEZGWn8888/xueff24AxuDBg40PPvjAuOeee4xDhw4ZtWrVMn777TfDMAxj0qRJRkxMjNGtWzcjMzPTeOeddwyz2WykpaXZaxs5cqSxd+/eIjXn5+cbAQEBxieffGI/duDAAQMwRowYUeT6NWvWGMOGDTNuvvlmY9asWcbx48eNjz/+2EhISDAMwzCys7ONFi1aGPHx8Q733XLLLUZ0dLRhtVrtx2644QbDZrMVeY0XX3zRaNKkiWs/dBEpFV9PBygRqVr+/PNP1q9fb5/Fk5+fzxVXXEFKSgqdOnWiYcOGXHbZZVxzzTVAYSvDn3/+CcDSpUtJS0vjqaeeAuDEiRO0adOGtLQ0unfvDsBNN93EddddR9++fVm2bBkZGRk0adIEgIsvvpiFCxeyfft2/P396devHw0bNuTNN99k/PjxFBQUcOzYMaKjo4vUvXPnTnJycqhfv7792J49ewA466yzilx/5ZVX8s033zBjxgzmzZtHSEgIzzzzDJ07d+aJJ54ACruB8vPzHe4bM2YMb7/9NsuXL6dPnz788MMPnH/++ZhMpiKv0bBhQ3bv3k1mZibBwcGl+BREpCQKMCLiIDk5meDgYO6//377senTp9v/3Ww2Yzb/v/c5KCiIvLw8AJKSkoiNjXW493ShoaH2f69Tpw4AaWlpREVFERQURHp6OpmZmfj7++Pj48Pw4cOZP38+48eP58svv+TKK690+twjR47Y6zkpJCQEoMiYlZPq1avnUEdSUhLPPPMM559/frH1X3rppZx33nm8+OKL9OnThzfeeMPh53Oqk7WkpqYqwIi4mcbAiIiD/Px8du/e7TD2A/4fEEq69+eff3Y4ZrVaizzrpC5dutCvXz/ee+89AH744Qd69OhBeHi4/Zrhw4ezc+dO1q1bx7Jly+jfv7/TZwUGBtprOKlVq1bUrVuXdevWlVh7cfU7e9933303X375JZs3b+b48eMOg3ed8fHxcen1RcR1CjAi4iAmJoacnByef/55+7E9e/bw3XfflXhvy5Yt+e677/j666/tx95//317C83pTCYT5513Hr/++isvvPACAQEB9gG7J0VFRdG7d2+ef/55srKyiIiIcPqsk11Hp4YlHx8fpkyZwl9//cW7777rUv3PPfccubm59mPOpkEPHjyY0NBQ+vfvz0033VTs806cOIGPj4/DrCgRcQ8FGBFxEBsbS48ePZgxYwZjx47lP//5D/Hx8Vx11VVAYSuFzWZzuKegoACAW2+9lTp16nDDDTcwa9YsZs6cyXfffUedOnUw/l1y6tQws2PHDl544QUuvPBC6tWrR4sWLfjll1/s1540atQoPv74Y3r27Fls3VFRUZxzzjns3bvX4fjEiRO5/fbbufvuux2CFUB6errD9+PHjycpKYlu3bqRkJDAoEGDiI2NLfJaISEh3HbbbeTm5tKrV69ia9qzZw+tW7fGYrEUe42IlJGnRxGLSNXz999/G9dee60RGBhoXHjhhfZZQp988okREhJidOjQwdi4caOxbt06o3HjxsbZZ59tbNiwwTAMw/j++++N8847zwgJCTFuvvlmIyMjw8jJyTFmzpxpAEaPHj2MxMREwzAM48SJE8b5559v1KlTx/D39zcAAzBGjhzpUI/NZjPi4uKMzMzMM9YdHx9v3H///U7Pvfvuu8ZVV11l9O3b17jrrruMfv36GTfffLPxzTffOFw3a9Ys46yzzjIaNmxovPzyy8W+VmJiojF16tQz1jN8+HBjxowZZ7xGRMpGK/GKiMf8+eefvP/++0yZMgUoXEn3yJEjxMfH8/bbb9uvy8vL45577mHevHklPm/QoEFs3LixQut2VevWrVmzZo3DzCgRcQ91IYmIx9x7771ceuml9u9NJhN16tShQ4cODtctX778jF01JzVr1oyePXs6Xeq/sq1cuZIhQ4YovIhUEAUYEfGY1NRUpk6dynfffUdGRgYHDhzg2WefpX379litVm6//XYefPBB3njjjTOOfznVQw89xLp160hNTa3g6ot38OBBvv/+eyZPnuyxGkS8nbqQRMRj9u/fz9SpU/n8889JT0+nY8eOPProo1xxxRVkZWVx3nnnERoayrvvvkuLFi1cfq7NZmP16tVcffXVFVh98T7++GOuv/56p4vbiYh7KMCIiIhItaMuJBEREal2FGBERESk2lGAERERkWpHAUZERESqHQUYERERqXYUYERERKTa+R+ncOjtZ++E9QAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def custom_model(x, *params):\n",
    "    return linear(x, *params[0:2]) + gauss(x, *params[2:5]) + gauss(x, *params[5:8]) + gauss(x, *params[8:11])\n",
    "\n",
    "guess = np.array([-10, 500, 200, 20, 3, 200, 40, 3, 200, 68, 3], dtype=float)\n",
    "\n",
    "fit_energy=np.linspace(energy[0], energy[-1], len(energy)*100)\n",
    "coeff, cov = curve_fit(custom_model, energy, counts, p0=guess, sigma=errors)\n",
    "coeff_error = np.sqrt(np.diag(cov))\n",
    "param_names = [\"Continuum slope\", \"continuum intercept\",\n",
    "               \"peak0_A\", \"peak0_mu\", \"peak0_sigma\",\n",
    "               \"peak1_A\", \"peak1_mu\", \"peak1_sigma\",\n",
    "               \"peak2_A\", \"peak2_mu\", \"peak2_sigma\"]\n",
    "\n",
    "for name, param,error in zip(param_names, coeff, coeff_error):\n",
    "    print(\"{0:30s}: {1:6.2f} +- {2:6.2f}\".format(name, param, error))\n",
    "\n",
    "plt.figure(figsize=[8,6])\n",
    "plt.errorbar(energy,counts,errors, fmt='.')\n",
    "plt.xlabel(\"energy (GeV)\")\n",
    "plt.ylabel(\"flux (counts)\")\n",
    "plt.plot(fit_energy, custom_model(fit_energy, *coeff));"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
