diff --git a/notebooks/leastSquaresFits_minimize.ipynb b/notebooks/leastSquaresFits_minimize.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1eea5b6a50fe27311b0c333f9aad331f6d199bd4 --- /dev/null +++ b/notebooks/leastSquaresFits_minimize.ipynb @@ -0,0 +1,305 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Least squares Fit\n", + "\n", + "A couple of examples to show how to use: \n", + "- scipy.optimize.minimize\n", + "\n", + "https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html \n", + "\n", + "(Minimization of scalar function of one or more variables. It can handle constraints)\n", + "\n", + "See also:\n", + "https://matthew-brett.github.io/cfd2019/chapters/08/using_minimize" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from __future__ import print_function\n", + "import numpy as np\n", + "%matplotlib inline\n", + "import matplotlib.pyplot as plt\n", + "from scipy.optimize import minimize\n", + "from scipy.optimize import Bounds\n", + "from scipy.stats import norm, chi2, lognorm" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fit a gaussian\n", + "Consider a Gaussian. We are measuring some feature which has a Gaussian distribution in $x$. This could be an inhomogeneous spectral line for $x=E$ the energy of emitted photons. We are interested in the resonance frequency and the linewidth, i. e. we want to estimate them form our observations." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "def gaussian_parent(x, mu, sigma):\n", + " return norm.pdf(x, mu, sigma) \n", + "\n", + "def gaussian_sample(mu, sigma, sample_size):\n", + " return norm.rvs(mu, sigma, sample_size)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEYCAYAAABiECzgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8WklEQVR4nO3dd3xUZdr/8c+VAglFuggGBTFUkQhItaAYSIBQBESUYnlEV7Hg6oq77j7K+uy6+lPEhoKLoKv0IiV0BRVBmhEDSBUlgHSRnpBcvz/mhB1DCBOYycnMXO/X67ySOec+c74TwpUz95xz36KqGGOMKXoRbgcwxphwZQXYGGNcYgXYGGNcYgXYGGNcYgXYGGNcEuV2gKJQuXJlrVmzptsxjDFhavXq1ftVtUre9WFRgGvWrMmqVavcjmGMCVMi8lN+660LwhhjXGIF2BhjXGIF2BhjXBIWfcDGGP/IysoiIyODkydPuh2lWIqJiSEuLo7o6Gif2lsBNsb4LCMjg7Jly1KzZk1ExO04xYqqcuDAATIyMqhVq5ZP+1gXhDHGZydPnqRSpUpWfPMhIlSqVKlQ7w6sABtjCsWK77kV9mdjBdgYY1xifcDG5GPNmjV8+eWXfPPNN+zdu5cmTZrQokULEhMTueSSS9yOZ4pQWloau3btomPHjn5/bivAxnj55ZdfePzxx5k4cSIAcXFxXHrppQwfPpzMzEwuu+wy3nzzTXr06GFvxYPQ6dOniYoqXNlLS0tj1apVASnA1gVhjGPChAnUr1+f6dOnM3ToUHbu3MmOHTtYvXo1v/32G59//jnVqlWjV69edOvWjV9//dXtyGFp+/bt1KtXjwEDBnDttdfSs2dPjh8/ztChQ7n++uu55pprGDhwILmz/bRt25Y///nP3HzzzQwfPpzVq1dz880307RpUzp06MDu3bvPtHvmmWdo3rw5derU4csvvyQzM5O//e1vTJgwgYSEBCZMmODX12JnwMYAn3zyCX379qVVq1aMHj2aunXr/m57yZIladu2LStWrGD48OE8++yzJCcns2DBAsqUKeNSanc98cQTpKWl+fU5ExISeP3118/bbuPGjfz73/+mTZs23HfffbzzzjsMGjSIv/3tbwD069ePWbNmkZKSAsCvv/7KkiVLyMrK4uabb+bTTz+lSpUqTJgwgb/85S+MHj0a8Jwhr1ixgtTUVF544QUWLlzI0KFDWbVqFW+99ZZfXytYATaG6dOn079/f26++WZSU1OJjY09Z9uoqCj++Mc/ctVVV9GrVy+6dOnC7NmzC9zH+F+NGjVo06YNAH379uWNN96gVq1avPzyyxw/fpyDBw/SsGHDMwW4d+/egKdwp6enk5iYCEB2djbVqlU787y33347AE2bNmX79u0Bfx0BLcAikgQMByKB91X1pTzbxdneETgO3KOqa0QkBvgCKOlknKyq/+vs8zzwALDPeZo/q2pqIF+HCV1ffvklvXv35vrrr2fGjBk+F9Lu3bszZswY+vfvT58+fZg2bVrY9Qn7cqYaKHl/1iLCww8/zKpVq6hRowbPP//8767HLV26NOC5WaJhw4YsW7Ys3+ctWbIkAJGRkZw+fTpA6f8rYH3AIhIJvA0kAw2APiLSIE+zZCDeWQYCI5z1p4BbVbUxkAAkiUhLr/2GqWqCs1jxNRfk8OHD9O3blyuuuILU1FTKli1bqP379u3La6+9xqeffsq7774boJQmPz///POZIjpu3DhuuOEGACpXrszRo0eZPHlyvvvVrVuXffv2ndk3KyuLdevWFXissmXLcuTIET+m/69AfgjXHNiiqttUNRMYD3TN06Yr8KF6LAfKi0g15/FRp020s2gAs5ow9Pjjj5ORkcFHH31EhQoVLvg5OnTowFNPPcWmTZv8nNCcS/369Rk7dizXXnstBw8e5A9/+AMPPPAAjRo1olu3blx//fX57leiRAkmT57MM888Q+PGjUlISODrr78u8Fi33HIL69evD8iHcKhqQBagJ55uh9zH/YC38rSZBdzg9XgR0Mz5PhJIA44C//Jq8zywHVgLjAYqnOP4A4FVwKorrrhCjfE2ZcoUBfS555676OfauXOnVqhQQZs3b65ZWVl+SFd8rV+/3u0I+uOPP2rDhg3djnFO+f2MgFWaT50K5Blwfh1iec9iz9lGVbNVNQGIA5qLyDXO9hFAbTxdE7uBV/M7uKqOVNVmqtqsSpWzZgIxYezXX3/lwQcfpEmTJmc+Nb8Y1atXZ8SIEaxYsYJhw4b5IaEJF4EswBlADa/HccCuwrZR1V+BxUCS83iPU5xzgFF4ujqM8dk///lPDhw4wPvvv+/zsIHn07t3bzp37syLL77I/v37/fKcJn81a9YkPT3d7Rh+EcgCvBKIF5FaIlICuBOYkafNDKC/eLQEDqvqbhGpIiLlAUQkFrgN+MF5XM1r/+5AaPxLmCLx008/MXz4cPr168d1113n1+f+17/+xdGjR3nxxRf9+rwmdAXsMjRVPS0ig4B5ePpzR6vqOhF5yNn+LpCK5xK0LXguQ7vX2b0aMNa5kiICmKiqs5xtL4tIAp6uiu3Ag4F6DSY0DFuwieGLNp95fNkTU1jirB+cWMdvx2nQoAH3338/77zzDo8++ii1a9f223ObEJVfx3CoLU2bNvWx+9yEssQ/f6CADhkyJGDH2LVrl5YuXVp79eoVsGO4qTh8CFfcFeZDOLsTzoSN9BkjKVG6HEOGDAnYMapVq8ZTTz3FCy+8QFpaGgkJCQE7VnGX951Hrsfbxfv1nUcws8F4TFhI+24tJyrGc/l9b/Lest2czs4J2LGeeOIJypQpwyuvvBKwYwSDwYl12P5SJ1rUqkiLWhXZ/lIntr/UKaiKb9u2bVm1alXAnt/OgE1YeHTUfMo260JOdAyjl/6IAE8n1fPrMbzP+Co9Mp6lQM0hs+2Mz5yTnQGbkLdjxw62/BZBRHQMACezcli69YDfj+N9xteo/Gl2vdadLsfnhXXxPZ2dw88Hj5O+8zAvz/3BL+88jh07RqdOnWjcuDHXXHMNEyZMKHAoysGDB3PTTTdRv359Vq5cye233058fDzPPfcccO7hLfOaP38+rVq1okmTJvTq1YujR4+e1aawrACbkDds2DBO/bwWce4DiomOoE3tSgE9ZqkKl3L33Xfz/vvvh/V1wa8t2MSe305yLDOb0Ut/ZNiCi79de+7cuVSvXp3vvvuO9PR0kpKSGDRoECtXriQ9PZ0TJ04wa9asM+1LlCjBF198wUMPPUTXrl15++23SU9PZ8yYMRw44PlDvHHjRgYOHMjatWu55JJLeOedd353zP379/Piiy+ycOFC1qxZQ7NmzXjttdcu+rVYATYh7dChQ4wcOZKONXKoVi6W0iUiub9NrYCdlXqf8VVpdz8nTp7i7bffDsixgsHXWw+Q49z/6q93Ho0aNWLhwoU888wzfPnll5QrV47PP/+cFi1a0KhRIz777LPfDbDTpUuXM/s1bNiQatWqUbJkSa666ip27NgBnD285VdfffW7Yy5fvpz169fTpk0bEhISGDt2LD/99NNFvxbrAzYhbcyYMRw7doynn/oj/1h+nBr4v+/XW+4ZX47CjI1HaHrP//LuuyP485//7Le77oJJ69qV+C7jV1T9986jTp06rF69mtTUVJ599lnat2/P22+/fc6hKHOHmIyIiDjzfe7j3CEn8xve0puqkpiYyLhx4y46vzc7AzYhS1V57733aNWqFY0bNy6SY+Y94ytzVRN++eUXZszIexNoeHgysQ7VLonx6zuPXbt2UapUKfr27ctTTz3FmjVrgPMPRVmQcw1vmatly5YsXbqULVu2AHD8+HG/jH5nZ8AmZC1ZsoSNGzcy4JWJ1Bwy+8z63O8DcXVC3jO+5KY12VajBu+99x49evTw67GCQVRkBDUqlvLrO4/vv/+ep59+moiICKKjoxkxYgTTp0+nUaNG1KxZ85xDURYkd3jLBx98kPj4eP7whz/8bnuVKlUYM2YMffr04dSpUwC8+OKL1Klzkb8/+d2dEWqL3QkXnu68804tX768Hj9+vMiOmXU6W1v9Y6E2+OscfXnOBs06na1Dhw5VQDdv3lxkOQLlQu6Eu+Pdr/WOd78OQBr/8PfwlsVlOEpjXLN3716mTJnCgAEDinS+ttwzvmsuL8fTSfWIiozg/vvvJzIyklGjRhVZjuJg2IJN1Bwym29+PMg3Px6k5pDZ1Bwy2y9XQoQK64IwIWnMmDFkZWXx4IPuj9VUvXp1unTpwgcffMDf//53SpQo4XakIjE4sU5QXAPt5vCWVoBNyFFVRo0axY033kj9+vWL7Lh5xz7w7mt+8MEHmTZtGtOnT+eOO+4oskyBoKphNwGpr1QLN3OaFWATcpYvX86WLVv4y1/+UqTHLeiMLyfnauLi4vjoo4+CugDHxMRw4MABKlWqZEU4D1XlwIEDxMTE+LyPFWATcj766CNiY2OL1VUHERER3HXXXbz22mvs27ePYJ0mKy4ujoyMDPbt2+d2lGIpJiaGuLg4n9tbATYhJTMzkwkTJtCtW7dCTzMfaP369ePll19m/PjxPProo27HuSDR0dHUqlXL7Rghw66CMCElNTWVgwcP0rdvX7ejnOWaa64hISGBjz76yO0oppiwAmxCykcffcSll15K+/bt3Y6Sr379+rFy5Uo2btzodhRTDFgBNiHj0KFDzJo1iz59+hAVVTx71/r06UNERISdBRvACrAJIRMnTiQzM5N+/fq5HeWcqlWrRmJiIv/5z3/IyQncrBwmOAS0AItIkohsFJEtInLWRFzOdPRvONvXikgTZ32MiKwQke9EZJ2IvOC1T0URWSAim52vFQL5GkzwmDBhAnXq1KFJkyZuRynQXXfdxU8//cSKFSvcjmJcFrAC7Ewp/zaQDDQA+ohIgzzNkoF4ZxkIjHDWnwJuVdXGQAKQJCItnW1DgEWqGg8sch6bMLdnzx6WLFnCHXfcUeyvT+3SpQvR0dFMmjTJ7SjGZYE8A24ObFHVbaqaCYwHuuZp0xX40BmvYjlQXkSqOY9z5/uIdhb12mes8/1YoFsAX4MJElOnTiUnJycobnIoX748HTp0YNKkSdYNEeYCWYAvB3Z4Pc5w1vnURkQiRSQN2AssUNVvnDZVVXU3gPP1Uv9HN8Fm0qRJ1K1bl2uuucbtKD7p1asXO3bssG6IMBfIApzf+8C8N0qfs42qZqtqAhAHNBeRQv3PEpGBIrJKRFbZXTuhLZi6H3J16dKFEiVKMHHiRLejGBcFsgBnADW8HscBuwrbRlV/BRYDSc6qPSJSDcD5uje/g6vqSFVtpqrNgvW2T+Ob3O6HXr16uR3FZ+XLl6d9+/ZMnjzZuiHCWCAL8EogXkRqiUgJ4E4g77wsM4D+ztUQLYHDqrpbRKqISHkAEYkFbgN+8NpngPP9AODTAL4GEwQmTZpEvXr1gqb7Idcdd9zBjh07+Oabb87f2ISkgBVgVT0NDALmARuAiaq6TkQeEpGHnGapwDZgCzAKeNhZXw34XETW4inkC1Q1d57pl4BEEdkMJDqPTZjau3cvS5YsoVevXkHT/ZArtxvCroYIX1LY8SuDUbNmzXTVqlVuxzAB8P777/PAAw+QlpZWZBNv+lOnTp1Yv34927ZtC7o/IMZ3IrJaVZvlXW93wpmgNn36dGrVqsW1117rdpQL0r17d7Zv387atWvdjmJcYAXYBK0jR46wcOFCunXrFrRnjykpKYgI06dPdzuKcYEVYBO05s6dy6lTp+jWrZvbUS5Y1apVad26tRXgMGUF2ASt6dOnU7lyZdq0aeN2lIvSvXt30tLS2L59u9tRTBGzAmyCUmZmJrNnz6ZLly5ERka6HeeidO3quUPfzoLDjxVgE5QWL17M4cOHg7r7IdfVV1/NNddcYwU4DFkBNkFp+vTplC5dmttuu83tKH7RvXt3vvzyS/bv3+92FFOErACboKOqzJgxgw4dOhAbG+t2HL/o2rUrOTk5pKamuh3FFCErwCbopKWlsXPnTlJSUtyO4jdNmjShevXqzJw50+0opghZATZBZ9asWYgIHTt2dDuK34gInTp1Yt68eWRmZrodxxQRK8Am6MycOZMWLVpw6aWhNRR0SkoKR44c4YsvvnA7iikiVoBNUPnll19YuXIlnTt3djuK37Vr146YmBjrhggjVoBNUJk9ezZASPX/5ipVqhTt2rVj5syZhMMgWcYKsAkys2bNokaNGjRq1MjtKAHRuXNnfvzxRzZs2OB2FFMErACboHHy5Enmz59/ZgCbUJTbtWLdEOHBCrAJGosXL+b48eMh2f2QKy4ujuuuu84KcJiwAmyCxuzZs4mNjaVt27ZuRwmoTp06sWzZMg4ePOh2FBNgVoBNUFBVUlNTz1wpEMo6duxITk4O8+fPdzuKCTArwCYobNq0iW3btoXUzRfn0rx5cypVqmS3JYcBK8AmKOQWo3AowJGRkSQlJTFnzhybsj7EBbQAi0iSiGwUkS0iMiSf7SIibzjb14pIE2d9DRH5XEQ2iMg6EXnca5/nRWSniKQ5S+j/jzSkpqbSsGFDrrzySrejFImOHTuyf/9+bDLZ0BawAiwikcDbQDLQAOgjIg3yNEsG4p1lIDDCWX8a+KOq1gdaAo/k2XeYqiY4i71PC3FHjx5lyZIlYXH2m6tDhw6IiHVDhLhAngE3B7ao6jZVzQTGA13ztOkKfKgey4HyIlJNVXer6hoAVT0CbAAuD2BWU4wtWrSIrKyssCrAlSpVomXLllaAQ1wgC/DlwA6vxxmcXUTP20ZEagLXAd94rR7kdFmMFpEKfktsiqXU1FTKli0b9HO/FVbHjh1ZuXIle/bscTuKCZBAFuD8blXKe4N7gW1EpAwwBXhCVX9zVo8AagMJwG7g1XwPLjJQRFaJyKp9+/YVMropLnIvP2vfvj3R0dFuxylSuWf8c+fOdTmJCZRAFuAMoIbX4zhgl69tRCQaT/H9WFWn5jZQ1T2qmq2qOcAoPF0dZ1HVkaraTFWbValS5aJfjHFHeno6GRkZJCcnux2lyCUkJFC1alXmzJnjdhQTIIEswCuBeBGpJSIlgDuBGXnazAD6O1dDtAQOq+pu8dzo/29gg6q+5r2DiFTzetgdSA/cSzBuyz37S0pKcjlJ0YuIiCApKYn58+eTnZ3tdhwTAAErwKp6GhgEzMPzIdpEVV0nIg+JyENOs1RgG7AFz9nsw876NkA/4NZ8Ljd7WUS+F5G1wC3A4EC9BuO+OXPm0KhRIy6/PDw/g01KSuLQoUOsWLHC7SgmAKIC+eTOJWKpeda96/W9Ao/ks99X5N8/jKr283NMU0wdOXKEr776isGDw/dvbPv27YmIiGDu3Lm0atXK7TjGz+xOOFNsffbZZ2RlZYVl90OuihUr0qJFC+sHDlFWgE2xNWfOHMqUKRN2l5/llZSUxKpVq7CreUKPFWBTLKkqc+bMoV27dpQoUcLtOK5KTk5GVW10tBBkBdgUSz/88AM///xzWF5+llfTpk2pXLmyXQ8cgqwAm2Ipt88znPt/c0VERNChQwfmzZtno6OFGCvApliaO3cu9erVC5vRz84nOTmZffv2sXr1arejGD+yAmyKnWPHjrFkyRLrfvDSvn17RIR58+a5HcX4kRVgU+wsWbKEzMxMK8BeqlSpQtOmTa0fOMRYATbFzty5c4mNjeXGG290O0qxkpSUxLJlyzh06JDbUYyfWAE2xc7cuXO55ZZbQn7yzcJKSkoiJyeHRYsWuR3F+IkVYFOsbN26lc2bN9vVD/lo0aIF5cqVs26IEHLeAiwidURkkYikO4+vFZHnAh/NhKPcD5msAJ8tKiqKxMRE5s6di2cYFRPsfDkDHgU8C2QBqOpaPENLGuN3c+fO5aqrruLqq692O0qxlJSUxM6dO1m3bp3bUYwf+FKAS6lq3rHwTgcijAlvp06d4rPPPiMpKQnPkNAmrw4dOgA2S0ao8KUA7xeR2jhTBYlITzxTARnjV0uXLuXYsWNniow5W1xcHNdcc42NjhYifCnAjwDvAfVEZCfwBPBQgXsYcwHmzJlDdHQ0t956q9tRirWkpCS+/PJLjh496nYUc5F8KcCqqrcBVYB6qnqDj/sZU6BhCzZRc8jsM8ukyLZUf3Iao5blnTrQeEtKSiIrK4vPP//c7SjmIvlSSKcAqOoxVT3irJscuEgmXAxOrMP2lzrRolZFGpU/zU//6sygyhsYnFjH7WjF2g033EDp0qWtHzgEnHNKIhGpBzQEyonI7V6bLgHsCnnjF6ezc/j54HH27TtB+Zv6k9je+n/Pp2TJkjTo/TSzyzZj9pDZv9v2eLt4+wMWRAqaE64u0BkoD6R4rT8CPBDATCaMvLZgE3t+O0lOibJccn1X5u2KpvG1bqcq/u5pWoVHHulM8gvjKVv1CiY8aPPFBaNzFmBV/RT4VERaqeqyIsxkwsjXWw+Q49xTIFEl+XrrAXcDBYncG1V2r1tO2apXuJzGXChf+oC/FZFHROQdERmdu/jy5CKSJCIbRWSLiAzJZ7uIyBvO9rUi0sRZX0NEPheRDSKyTkQe99qnoogsEJHNztcKPr9aU+y0rl0JUc8g49GitKldyeVEweGqq64ivm499p2OJX3nYV6e+wOns22w9mDjSwH+CLgM6AAsAeLwdEMUSEQigbeBZKAB0EdEGuRplgzEO8tAYISz/jTwR1WtD7QEHvHadwiwSFXjgUXOYxOknkysQ+TeHzi1ayP9W9aw/stCuLLTw0TUbMqxzGxGL/2RYQs2uR3JFJIvBfhqVf0rcExVxwKdgEY+7Ncc2KKq21Q1ExgPdM3TpivwoXosB8qLSDVV3a2qawCcKy82AJd77TPW+X4s0M2HLKaYioqM4OBnozm97EP+2rUxUZF2haOvsipeRUS05/Pwk1k5LLXum6Djy297lvP1VxG5BigH1PRhv8uBHV6PM/hvEfW5jYjUBK4DvnFWVVXV3QDO10vzO7iIDBSRVSKyyqbzLp5yrwO+pM//o3SP/ztzPbCdyfkmsfEV5GSdAiAmOsK6b4JQQVdB5Brp9LM+B8wAygB/9WG//G7mzzuEU4FtRKQMnuuQn1DV33w45n+fRHUkMBKgWbNmNnRUMTQ4sQ4Vdy3jnnvuYc2aNVx33XVuRwoqzyQ3ZMQHn5Bd7nIeueM2674JQgWeAYtIBPCbqh5S1S9U9SpVvVRV3/PhuTOAGl6P44C8tzids42IROMpvh+r6lSvNntEpJrTphqw14csppiaM2cOl112GQkJCW5HCTpRkRFUiTrJrg8e5fb4aOu+CUIF/oupag4w6AKfeyUQLyK1RKQEniEsZ+RpMwPo71wN0RI4rKq7xTMU1r+BDar6Wj77DHC+HwB8eoH5jMtOnz7N/PnzbfSzC5DbfXOo1m1c+cws2o36wbpvgpAvXRALROQpYAJwLHelqh4saCdVPS0ig4B5QCQwWlXXichDzvZ3gVSgI7AFOA7c6+zeBugHfC8iac66P6tqKvASMFFE7gd+Bnr58kJN8bNixQoOHTpkk29egMGJdc50OVx99dXUrVuX2bNnn2cvU9z4UoDvc74+4rVOgavOt6NTMFPzrHvX63vN87y5678i//5hVPUA0O68qU2xN2fOHCIiIkhMTHQ7SlBLTk7m3//+NydPnrR59ILMeTuNVLVWPst5i68x5zNnzhxatWpFhQp2L83FSE5O5sSJEyxZssTtKKaQrNfeuOKXX35h9erV1v3gB23btiUmJsYGaQ9CVoCNK3In3+zYsaPLSYJfqVKlaNu2rRXgIGQF2LgiNTWVyy67jMaNG7sdJSQkJyezadMmtm7d6nYUUwi+TEu/yJd1xvgqKyuLefPm0bFjRyIi7BzAHzp16gRgV0IEmXP+9otIjIhUBCqLSAVnFLKKzq3B1YssoQk5y5Yt4/Dhw2eKhrl4tWvXpm7duqSmpp6/sSk2Cjr9eBBYDdRzvuYun+IZ5cyYCzJ79myio6O57bbb3I4SUjp16sTixYs5duzY+RubYuGcBVhVh6tqLeAp5xbk3EvQGqvqW0WY0YSY2bNnc+ONN3LJJZe4HSWkdOzYkVOnTrFokfUQBgtfrgN+U0Rai8hdItI/dymKcCb0/PTTT6xbt866HwLgxhtvpGzZstYPHETOeyeciHwE1AbSgGxntQIfBi6WCVW5fZRWgP2vRIkSJCYmkpqaiqra+BpBwJdbkZsBDZzbho25KLNnz6Z27drUqWNDJwZCp06dmDp1Kt9//z3XXmuzmxZ3vlwDlI5nSiJjLsqJEyf47LPP6Nixo52dBUjunYXWDREcfCnAlYH1IjJPRGbkLoEOZkLPZ599xokTJ+jcubPbUUJWtWrVaNq0KbNmzXI7ivGBL10Qzwc6hAkPM2bMoEyZMtx8881uRwlpKSkpvPDCC+zdu5dLL813xi5TTPhyFcSS/JaiCGdCh6oya9YsOnToQMmSJd2OE9JSUlJQVbspIwj4civyERH5zVlOiki2iBRqfjZj1qxZw65du0hJSXE7Ssi77rrruPzyy5k5c6bbUcx5nLcLQlXLej8WkW54ppw3xmczZ85ERGz0syIgInTu3JmPP/6YU6dO2TuOYqzQI6Go6nTgVv9HMaFs5syZtGrViipVqrgdJSykpKRw9OhRFi9e7HYUUwBfuiBu91p6ishLnD29vDHntHPnTtasWWPdD0Xo1ltvJTY21rohijlfzoBTvJYOwBGgayBDmdCSe0mUFeCiExsbS2JiIjNnzsTuoSq+fLkK4l6v5QFV/T9V3evLk4tIkohsFJEtIjIkn+0iIm8429eKSBOvbaNFZK+IpOfZ53kR2Skiac5inYrF3MyZM6lVqxYNGjRwO0pYSUlJ4eeff2bt2rVuRzHn4EsXRJyITHOK4R4RmSIicT7sF4ln2MpkoAHQR0Ty/g9MBuKdZSAwwmvbGCDpHE8/TFUTnMWutSnGjhw5wsKFC+nWrZvd/VbEUlJSEBE+/fRTt6OYc/ClC+IDYAaeQdgvB2Y6686nObBFVbepaiYwnrO7LroCH6rHcqC8iFQDUNUvgIO+vQxTXM2bN49Tp07RrVs3t6OEnapVq9K6dWumT5/udhRzDr4U4Cqq+oGqnnaWMYAvH2VfDuzwepzhrCtsm/wMcrosRotIvnOai8hAEVklIqv27dvnw1OaQJg2bRqVK1emdevWbkcJS926dePbb79l+/btbkcx+fClAO8Xkb4iEuksfYEDPuyX3/vNvJ8G+NImrxF4hsdMAHYDr+bXSFVHqmozVW1mlz65IzMzk9mzZ5OSkkJUlC93vRt/69rV86bTuiGKJ18K8H3AHcAveApeT2fd+WQANbwexwG7LqDN76jqHlXNVtUcYBR2U0ixtWTJEg4fPmzdDy6Kj4+nYcOG1g1RTPlyFcTPqtpFVauo6qWq2k1Vf/LhuVcC8SJSS0RKAHfi6Uv2NgPo71wN0RI4rKq7C3rS3D5iR3c8w2WaYmj69OmUKlWKxMREt6OEte7du/PFF19w4IAvb1xNUfLlKoixIlLe63EFERl9vv1U9TQwCJgHbAAmquo6EXlIRB5ymqUC24AteM5mH/Y6zjhgGVBXRDJE5H5n08si8r2IrAVuAQb78DpNEcvJyWH69OkkJSURGxvrdpyw1q1bN3JycuymjGLIl465a1X119wHqnpIRK7z5cmdS8RS86x71+t7BR45x759zrG+ny/HNu5atWoVu3btsu6HYqBJkybExcUxbdo07rnnHrfjGC++9AFHeF9pICIV8a1wmzA2depUoqKibO63YkBE6N69O/Pnz+fo0aNuxzFefCnArwJfi8jfRWQo8DXwcmBjmWCmqkyePJl27dpRsWJFt+MYoGfPnpw8edKmKipmfPkQ7kOgB7AH2AfcrqofBTqYCV5paWls3bqVnj17uh3FONq0aUPVqlWZPHmy21GMF5+6ElR1PbA+wFlMiJg8eTKRkZHW/1uMREZGcvvttzN27FiOHTtG6dKl3Y5kuIDxgI0piKoyadIkbrnlFipXrux2HOOlV69eHD9+nLlz57odxTisABu/Sk9PZ/Pmzdb9UAzdeOONVKlShUmTJrkdxTisABu/mjx5MhEREdb9UAxFRUXRvXt3Zs2axYkTJ9yOY7ACbPxs8uTJ3HTTTVStWtXtKCYfPXv25NixY9YNUUxYATZ+8/3337N+/Xp69erldhRzDm3btqVy5cpMmDDB7SgGK8DGj8aNG0dkZKQV4GIsOjqanj17MmPGDLspoxiwAmz8QlUZP348t912m818XMz16dOHEydOMGNG3rGxTFGzAmz84ptvvuHHH3+kT598h/AwxcgNN9xAXFwc48ePdztK2LMCbPxi3LhxlCxZ0q5+CAIRERH07t2buXPncvCgzfrlJivA5qJlZ2czceJEOnbsSLly5dyOY3zQp08fsrKymDp1qttRwpoVYHPRFi9ezC+//GLdD0GkSZMmxMfHM27cOLejhDUrwOaiffLJJ5QpU4bOnTu7HcX4SETo06cPn3/+OTt37nQ7TtiyAmwuyvHjx5k0aRI9e/a0mS+CTN++fVFVPv74Y7ejhC0rwOaiTJ8+nSNHjjBgwAC3o5hCio+Pp3Xr1owdOxbP5DSmqFkBNhdl7NixXHnlldx0001uRzEXoH///qxfv541a9a4HSUsWQE2F2znzp0sXLiQfv36ERFhv0rBqHfv3pQsWZKxY8e6HSUsBfR/jYgkichGEdkiIkPy2S4i8oazfa2INPHaNlpE9opIep59KorIAhHZ7HytkPd5TdH4z3/+Q05ODv3793c7irlA5cuXp2vXrnzyySdkZma6HSfsBKwAi0gk8DaQDDQA+ohIgzzNkoF4ZxkIjPDaNgZIyuephwCLVDUeWOQ8NkVMVRk7diytW7cmPj7e7TjmIgwYMIADBw6Qmpp6/sbGrwJ5Btwc2KKq21Q1ExgPdM3TpivwoXosB8qLSDUAVf0CyO82na5A7vulsUC3QIQ3BVu5ciUbNmyws98Q0L59ey677DLGjBnjdpSwE8gCfDmww+txhrOusG3yqqqquwGcr5fm10hEBorIKhFZtW/fvkIFN+c3atQoSpUqxZ133ul2FHORoqKi6N+/P7NmzWL37t1uxwkrgSzAks+6vNe6+NLmgqjqSFVtpqrNbHQu/zpy5Ajjxo2jd+/edutxiPif//kfsrOz+eCDD9yOElYCWYAzgBpej+OAXRfQJq89ud0Uzte9F5nTFNK4ceM4duwYAwcOdDuK8ZP4+HhuvfVWRo0aRU5OjttxwkYgC/BKIF5EaolICeBOIO8ApDOA/s7VEC2Bw7ndCwWYAeRe9T8A+NSfoc35jRw5kkaNGtGiRQu3oxg/euCBB9i+fTsLFy50O0rYCFgBVtXTwCBgHrABmKiq60TkIRF5yGmWCmwDtgCjgIdz9xeRccAyoK6IZIjI/c6ml4BEEdkMJDqPTRFZvXo1q1evZuDAgYjk14NkglX37t2pVKkSI0eOdDtK2IgK5JOraiqeIuu97l2v7xV45Bz75ju0lqoeANr5MaYphFGjRhETE0Pfvn3djmL8rGTJktxzzz0MHz6cPXv22MSqRcBuXzI+O3z4MB9//DG9e/emfPnybscxAfDAAw9w+vRpRo0a5XaUsGAF2PhszJgxHD16lEcffdTtKCZA6tatS/v27RkxYgRZWVluxwl5VoCNT3JycnjzzTdp3bo1TZs2dTuOCaDHHnuMXbt22WwZRcAKsPHJnDlz2Lp1K4899pjbUUyAJScnU7t2bd544w23o4Q8K8DGJ2+88QbVq1fn9ttvdzuKCbCIiAgeffRRvv76a1atWuV2nJBmBdic14YNG5g/fz4PP/ww0dHRbscxReCee+6hTJkyvPnmm25HCWlWgM15vf7665QsWdLufAsj5cqV45577mHcuHE2Z1wAWQE2Bdq9ezdjxozh3nvvxcbUCC9PPvkkOTk5vP76625HCVkSDnNBNWvWTK0vy3fDFmxi+KLNZ61/vF08gxPruJDIuOXuu+9mxowZ/Pzzz1SoYHMfXCgRWa2qzfKutzNgc5bBiXXY8n/JXFoqgszdm6n+6/ds+b9kK75h6E9/+hNHjx5lxIgR529sCi2gtyKb4PXagk3sPZpJiWrx7I/0nBU/nVTP7VimiH22N5Yrn5nFyN9g5JDZZ9bbuyH/sDNgk6+vNu+DCM/f58xsWLr1gMuJjBsGJ9ah81XR5GSdBCAmOoJH2ta24usnVoBNvqIP/vi7/3RtaldyOZFxS0ZWaSKiYwA4mZVjf4z9yAqwOcu/Zqez+kSl3/2ne3vxVoYt2ORyMuOG1rUrIeoZpD2KHPtj7EfWB2zOot/P4qd/PcuyZcto2bKl23GMy55MrMP0b3ey48fNZO34joF/edntSCHDLkMzv3P48GFq1apFq1atmD179vl3MCHNLkn0j3NdhmZnwOZ3Xn/9dQ4dOsTQoUPdjmKKgcGJdX5XaJOTk1m5ciX3/3mbi6lCh/UBmzN27drFK6+8Qo8ePWzISZOvF198kQMHDvDPf/7T7SghwQqwOeO5554jMzOTl16yafZM/po2bUrfvn0ZNmwY27dvdztO0LMCbABYs2YNY8aM4fHHH+fqq692O44pxv75z38SERHBM88843aUoBfQAiwiSSKyUUS2iMiQfLaLiLzhbF8rIk3Ot6+IPC8iO0UkzVk6BvI1hANV5cknn6RSpUo899xzbscxxVxcXBx/+tOfmDhxIkuXLnU7TlALWAEWkUjgbSAZaAD0EZEGeZolA/HOMhAY4eO+w1Q1wVlSMRdlypQpLFmyhKFDh1KuXDm345gg8PTTT1O9enUef/xxsrOz3Y4TtAJ5Btwc2KKq21Q1ExgPdM3TpivwoXosB8qLSDUf9zV+cPjwYR577DESEhJ44IEH3I5jgkTp0qV59dVXWb16NW+99ZbbcYJWIAvw5cAOr8cZzjpf2pxv30FOl8VoEcl3jDwRGSgiq0Rk1b59+y70NYS8Z599lj179jBq1CiiouyqROO73r17k5SUxHPPPceOHTvOv4M5SyALsOSzLu9dH+dqU9C+I4DaQAKwG3g1v4Or6khVbaaqzWwg8fwtW7aMd999l8cee4xmzc66RtyYAokI77zzDjk5OTzyyCOEw01d/hbIApwB1PB6HAfs8rHNOfdV1T2qmq2qOcAoPN0VppBOnTrFwIEDiYuL4+9//7vbcUyQqlWrFi+88AIzZ85k0qRJbscJOoEswCuBeBGpJSIlgDuBGXnazAD6O1dDtAQOq+rugvZ1+ohzdQfSA/gaQtZf/vIX0tPTee+99yhTpozbcUwQe+KJJ2jWrBl/+MMfbP64QgpYAVbV08AgYB6wAZioqutE5CERechplgpsA7bgOZt9uKB9nX1eFpHvRWQtcAswOFCvIVQtWrSIV199lYcffpjk5GS345ggFxUVxccff8zJkye59957ycnJcTtS0LDBeMLMwYMHufbaaylTpgxr1qyhVKlSbkcyIeK9997joYceYtiwYTzxxBNuxylWbE44Q05ODvfddx979uzh448/tuJr/GrgwIGkpKTwzDPPYCc8vrECHEb+8Y9/8Omnn/LKK6/YYDvG70SE0aNHc9lll3H77bezd+9etyMVe1aAw8ScOXP429/+xt13383jjz/udhwToipXrszUqVPZu3cvvXv35vTp025HKtasAIeBjRs3ctddd9G4cWNGjhyJSH6XWRvjH02bNuW9995j8eLF/PGPf7Trgwtgtz6FuF27dtGhQwdKlCjB1KlTrd/XFIkBAwaQlpbG66+/To0aNXjqqafcjlQsWQEOYYcPHyY5OZkDBw6wZMkSatWq5XYkE0ZeffVVdu3axdNPP81ll11G37593Y5U7FgBDlFHjx4lJSWFDRs2MHv2bJo0aXL+nYzxo4iICD788EP27dvHvffeyyWXXEKXLl3cjlWsWB9wCPrtt9/o0KEDX3/9Nf/5z39ITEx0O5IJUyVLlmTatGk0adKEHj16MGXKFLcjFStWgEPMoUOHSExMZMWKFUyYMIE77rjD7UgmzJUrV4758+fTvHlzevfuzfjx492OVGxYAQ4hW7dupXXr1nz77bdMmTKFHj16uB3JGMBThOfOncsNN9zAXXfdxSuvvGJXR2AFOGR89dVXtGjRgr1797JgwQLrazPFTtmyZZkzZw69evXiT3/6Ew888ACZmZlux3KVFeAgp6oMHz6cW2+9lYoVK7J8+XJuvvlmt2MZk6/Y2FjGjRvHc889x7///W9uueUWfvrpJ7djucYKcBA5nZ3Dy3N/oNvbS3l57g/s2befrl278sQTT5CUlMTy5cuJj493O6YxBYqIiODvf/8748eP5/vvvychIYGpU6ee9ft9Ojv0R1Wz0dCCxLAFmxi+aPNZ648sm8D/9mjGo48+ane4maCzdetW7rzzTjaXrEP5G+46a/vj7eIZnFjHhWT+ZaOhBbnBiXVIqFH+d+uiDmfw2Rt/5LHHHrPia4JS7dq1Wbp0KU/cFk/m7k2/25ZQo3xIFN+CWAEOEn+ftoa0Hb/+bt3pcnF8tjfWnUDG+EmJEiV4/vnn6Zt4PXra+VDudCaX6sGQv1LC7oQr5rZu3cprr73GBx98QNbpbG565GVirmxM23qXMTixDlGR9jfUBL9hCzYxMf0wElXCsyKqBPMzoMEdT/N8z+vp0aNHaM7araohvzRt2lSDSVZWls6YMUO7du2qERERWqJECb3vvvt006ZNbkczpkhkZmbqqFGjtE6dOgrolVdeqUOHDtUdO3a4He2CAKs0n9pkH8IVE9nZ2Xz55ZdMmTKFyZMn88svv1C1alXuu+8+Bg0aRPXq1d2OaEyRy8nJYcaMGbz55pt89tlnREREcOutt9KjRw+6d+9O1apV3Y7ok3N9CGcF2CWqypYtW/jiiy+YP38+Cxcu5ODBg8TExJCcnMyAAQPo2LEj0dHRbkc1pljYtm0bH3zwARMmTGDz5s2ICE2bNqV9+/bcdtttNG/enNKlS7sdM1+uFGARSQKGA5HA+6r6Up7t4mzvCBwH7lHVNQXtKyIVgQlATWA7cIeqHiooh9sFOCsrix9//JH09HTS0tL49ttvWb58Ofv37wegevXqtG/fno4dO5KcnGzTxBtTAFUlPT2dadOmMX/+fJYvX052djaRkZE0btyYZs2akZCQQOPGjalXrx4VK1b83f6ns3N4bcEmvt56gNa1K/FkEXyWUuQFWEQigU1AIpABrAT6qOp6rzYdgUfxFOAWwHBVbVHQviLyMnBQVV8SkSFABVV9pqAsgSzAWVlZ/LJ3H8MXbWF1xlHioo9R//RWMnb8zPbt29m2bRtbt249MzVLREQEdevWpUWLFrRu3Zo2bdpQv359u4zMmAt0+PBhli5dyrJly/j666/59ttvOXTov+dkVapU4eqrr6ZmzZr8dsWNpHPFWc/xeLurGZxYN2AZ3SjArYDnVbWD8/hZAFX9p1eb94DFqjrOebwRaIvn7DbffXPbqOpuEanm7F/gT66wBfiTTz5h7dq1nDx5khMnTnD8+HGOHz/O0aNHOXLkCL/99huHDx/m119/JbJxl3wvIM9cM43qB9OoWbMmderUoW7dutSvX59GjRoRG2uXjhkTKKrKjh07+O6779i4cSMbN25k27ZtbN++nUOXt+aS1neetc9vX48ncsNcKlSoQNmyZbnkkksoU6YMpUuXpnTp0sTGxhITE0NsbCyDBg0qdN/zuQpwIK/ruBzY4fU4A89Z7vnaXH6efauq6m4Apwhf6s/QAFOmTGHWrFnExMQQExNDqVKlzvxDlC1blqpVq1K+fHkqVKjAtlL1Sctn3sGnn/5TyF9EbkxxJCJcccUVXHHFFaSkpPxuW3Z2NkOnpzFuzS9kZkOUKI1j9lP3hhocatiTQ4cOceTIEY4cOUJGRgbHjx/n2LFjnDhx4swJ2d133+23D/8CWYDze0+d93T7XG182bfgg4sMBAYCXHHF2W85CjJ58uRCdQm8PPcHRi/9kZNZOcRER3B/m1pWfI0phiIjI/lrt+soU3oTS7ceoE3tSgxO7OhzH7C/ewwCWYAzgBpej+OAXT62KVHAvntEpJpXF8Te/A6uqiOBkeDpgihM8ML2xz6ZWAcBr39QK77GFFdRkRE8nVSPpy9gX39/VhPIArwSiBeRWsBO4E4gb2fpDGCQiIzH08Vw2Cms+wrYdwYwAHjJ+fppAF+DTy7mH9QYE74CVoBV9bSIDALm4bmUbLSqrhORh5zt7wKpeK6A2ILnMrR7C9rXeeqXgIkicj/wM9ArUK/BGGMCyW7EMMaYALPhKI0xppixAmyMMS6xAmyMMS6xAmyMMS6xAmyMMS4Ji6sgnOuKCzv3dWVgfwDiWIbgOr5lsAz+OP6Vqlol78qwKMAXQkRW5XfZiGUIr+NbBssQyONbF4QxxrjECrAxxrjECvC5jXQ7AJahOBwfLEMuy+Dn41sfsDHGuMTOgI0xxiVWgI0xxiVhU4BFZLSI7BWRdK91z4vIThFJc5aOXtueFZEtIrJRRDp4rW8qIt87296QQozQXJgMIlJJRD4XkaMi8lae5ymqDIkisto51moRudWFDM291n0nIt2LOoPX9iucf4+nLjZDIX8GNUXkhNf6d934GYjItSKyTETWOceMKcoMInK317o0EckRkYQizhAtImOdY20QZ77KC86gqmGxADcBTYB0r3XPA0/l07YB8B1QEqgFbAUinW0rgFZ4pk2aAyQHKENp4AbgIeCtPNuKKsN1QHXn+2uAnS5kKAVEOd/nzoASVZQZvLZPASZ5t7nQDIX8GdT0bufS70IUsBZo7DyuVNT/J/Ls1wjY5sLP4S5gvNfv5nag5oVmCJszYFX9AjjoY/OueH7Ip1T1RzwDxjcXzxRIl6jqMvX8xD8EugUig6oeU9WvgJPe64s4w7eqmjsV1DogRkRKFnGG46qaO+1pDM7cgEWZwTleN2Abnp9D7roLzlDY458jU1H+DNoDa1X1O2ffA6qaXdT/Dl76ALmzqRdlBgVKi0gUEAtkAr9daIawKcAFGCQia523IRWcdQXN1pyRz/pAZDgXtzL0AL5V1VNFnUFEWojIOuB74CGnIBdZBhEpDTwDvJCnbSAynOvfoZaIfCsiS0TkxgAe/1wZ6gAqIvNEZI2I/MmFDN564xTgIs4wGTgG7MYzI8//U9WDF5oh3AvwCKA2kIDnB/qqsz5gszUXIsO5FHkGEWkI/At40I0MqvqNqjYErgeedfoeizLDC8AwVT2ap72/M5zr+LuBK1T1OuBJ4BMRuSQAxy8oQxSeLrG7na/dRaRdEWcAPH+QgeOqmttnW5QZmgPZQHU83ZN/FJGrLjRDICflLPZUdU/u9yIyCpjlPDzXbM0Zzvd51wciw7kUaQYRiQOmAf1VdasbGbzabBCRY3j6o4syQwugp4i8DJQHckTkJJ4+Yb9lONfxnXcdp5zvV4vIVjxnpEX5M8gAlqjqfmdbKp5+0/8UYYZcd/Lfs9/cbEWV4S5grqpmAXtFZCnQDPjyQjKE9Rmw02+TqzuQ+xd1BnCn099ZC4gHVqjqbuCIiLR0PuHsz0XOylxAhnwVZQYRKQ/MBp5V1aUuZajl9LchIlcCdYHtRZlBVW9U1ZqqWhN4HfiHqr7l7wwF/AyqiEik8/1VeH4ftxXx7+M84FoRKeX8e9wMrC/q/xMiEoFnIt7xueuKOMPPwK3iURpoCfxwwRl8/aQw2Bc8fzF3A1l4/mLeD3yEp19xLZ6iW82r/V/wXP2wEa9PM/H8tUt3tr2FczdhgDJsx/PhwFGnfYOizAA8h6e/K81rubSIM/TD88FXGrAG6ObGv4XXfs/z+6sgLihDIX8GPZyfwXfOzyDFpd/Hvk6OdOBllzK0BZbn8zxF9ftYBs+VMOuA9cDTF5PBbkU2xhiXhHUXhDHGuMkKsDHGuMQKsDHGuMQKsDHGuMQKsDHGuMQKsDHGuMQKsDHGuMQKsAlZIvKgiOyW348h2yhPm1hngJtIPxxvsXiNHe2se0JE3hGREiLyRe5dfcaAFWAT2q4FnlPVBK/l+zxt7gOmqmq2H443Ds84Bd7uBMapaiawCM8oXsYAVoBNaGuE5xbmgtyNc8++eGae+EFE3heRdBH5WERuE5GlIrJZRJrn7iQifUVkhXNW/Z5zBj0Z6CwiJXOfD8+oWV85u013jmcMYAXYhLaGwAde3Q8DvTeKSAngKlXd7rX6amA4nrPnenhGv7oBeAr4s7NffTxnsm1UNQHP8IR3q+oBPLMiJDnPdScwQf97v386niE1jQHCfDhKE7pEpAawV1WvLaBZZeDXPOt+zO2mEM8g8ItUVUXkezxTAwG0A5oCKz0DXxGLZ6ok+G83xKfO1/tyn1g9M0hkikhZVT1yES/PhAgrwCZUXQv8cJ42J/BMc+TtlNf3OV6Pc/jv/xcBxqrqs5xtOvCaiDQBYlV1TZ7tJckzzZQJX9YFYUJVI85TgFX1EBDpzLBRGIvwDNB+KYCIVHTGKkY9s2YsBkbz+0HDEZFKwD71DOZtjBVgE7IaAf28+n+/FZEy+bSbj6eP12equh7PWMnzRWQtsADPjM25xgGN8Ro03HELkFqYY5nQZuMBm7AmItcBT6pqvyI41lQ8s4tsDPSxTHCwM2AT1lT1W+Bzf9yIURDniovpVnyNNzsDNsYYl9gZsDHGuMQKsDHGuMQKsDHGuMQKsDHGuMQKsDHGuMQKsDHGuOT/Aw8JzLK1reASAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<Figure size 360x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Create sample\n", + "## SAMPLE SIZE\n", + "sample_size = 5000\n", + "##################\n", + "\n", + "# Prepare toy data\n", + "mu = 1540 # True values that we will try to estimate\n", + "sigma = 11 # using a least-squares fit\n", + "\n", + "x_arr = np.linspace(1500, 1580, 101)\n", + "bins = 10\n", + "sample = gaussian_sample(mu, sigma, sample_size)\n", + "hist = np.histogram(sample, bins=bins, range=(1500, 1580))\n", + "bin_width = np.diff(hist[1])[0]\n", + "normalization = bin_width * sample_size\n", + "x = hist[1][:-1]+bin_width/2\n", + "y = hist[0]/normalization\n", + "y_errors = np.sqrt( (np.sqrt(hist[0]) / normalization) **2) #sqrt(x^2) just to take the positive value\n", + "\n", + "# Plot our toy measurement results\n", + "plt.figure(figsize=(5, 4))\n", + "plt.xlabel(r'$E$ (meV)')\n", + "plt.ylabel('count rate')\n", + "plt.plot(x_arr, gaussian_parent(x_arr, mu, sigma), '-', color='black', label='parent')\n", + "plt.errorbar(x, y, yerr=y_errors, fmt='.', ms=7, capsize=3, label='sample')\n", + "plt.legend()\n", + "plt.tight_layout()\n", + "\n", + "# Save data\n", + "data = np.vstack((x, y, y_errors))\n", + "np.savetxt('data', data)\n", + "np.savetxt('sample', sample)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Load data from disk. Format (3,12): (x, y, y_error) x N \n", + "data = np.loadtxt('data')\n", + "x = data[0, :]\n", + "y = data[1, :]\n", + "y_error = data[2, :]\n", + "# The sample used to generate\n", + "sample = np.loadtxt('sample')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": { + "scrolled": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbgAAAEoCAYAAAAqrOTwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABucklEQVR4nO2dd3hU1daH3zWT3ggBAoTeO9KLdAERBCk2sIGKWK792it69V5791NREbCgoCIqTUSQIi30TugQQkJ6z2Rm9vfHDJrLBVOYyZkz2e/zzHNmzpzym5PJ+c1ee+21RSmFRqPRaDT+hsVoARqNRqPReANtcBqNRqPxS7TBaTQajcYv0Qan0Wg0Gr9EG5xGo9Fo/BJtcBqNRqPxS0o1OBGZLiIpIrKzxLpXRWSviGwXkXkiEn2efS8TkX0ickBEHiuxPkZElopIgntZ3SOfRqPRaDQaN2Vpwc0ALjtr3VKgvVKqI7AfePzsnUTECrwPDAfaAhNEpK377ceAZUqpFsAy92uNRqPRaDxGqQanlFoJpJ+17hellN39ch1Q/xy79gAOKKUOKaVswNfAaPd7o4GZ7uczgTHll67RaDQazfnxRB/cLcCic6yvBxwv8fqEex1AbaVUEoB7GesBHRqNRqPR/EnAhewsIk8CduDLc719jnXlrgsmIlOAKQDh4eFdW7duXd5DaM7GmQGORLBEgbWh0Wr8Eqe7BN6ZpVIKVWKJAoVCnWMpAsHWAAIDrOf8J9J4CGcKOFLAGgsW/RvbCDZt2pSqlKrlreNX2OBEZCIwEhiszl3Q8gTQoMTr+sBJ9/NkEamrlEoSkbpAyvnOo5SaBkwD6Natm4qPj6+oZA2g8mahcl6AoOuR6u8hEmq0JFNTUFzMvuRUdp9KYc+pFHadSiEhJQ2bw1Gh450xtCLAabXSOCaapjVjaFYzhqY1XMvGNaoTFhTosc9QVVHKicp+Agq+RyJuRSLuNFpSlUNEjnrz+BUyOBG5DHgUGKCUyj/PZhuBFiLSBEgExgPXud/7EZgIvORezq+IDk35ULkfoXJfh+AhSPRbiAQZLclUZBUUuo3sNLtPpbD7VAqH0zL+bKVVCwmmTZ1Yru9+ETFhoQRYrARYLQRYLAS6l2fWBVos7vf+d5siu4PDaRkcSk3nYGo6+5JPs3TvgT/PA1CvWhRNalb/L+NrWjOGmLBQRHS7ryyIWCDqRZSyo3LfBEsMEnat0bI0HqRUgxOR2cBAoKaInACexZU1GQwsdf8zrVNK3SEiccAnSqkRSim7iNwNLAGswHSl1C73YV8C5ojIrcAx4GoPfy5NCZRSqNy3IO8DCBmFVHsJEd0C+DuSs3P/NLEzhpaYlf3n+7UjI2hbJ5bL2rSgbZ1Y2tSJJa5apMfMpUuDuP96bbPbOZqeycHUdA6lZXAwNZ3DqenMObaDgmL7n9tFh4bQvWF9hrdtycCWTQgP0j9i/g4RK1R7CeXMQGU/B9aGSHBvo2VpPISYabocHaIsP0opVM5/IH8GhF6NRD3v+qfW/A+Jmdn8vHMvP+3cS8LpNMAVMmwUU522dWrRpk4sbevE0rZOLWLCw4wV68apFKeyc1zGl5pOwuk0ViQc5nRuHiEBAQxo0YThbVsyoHkTHdb8G5QzB5U+HhzJSI25SEAToyVVCURkk1Kqm9eOrw3Of3H1MTwLBd9A2EQk8gkdvjqLzIJCFu/ez0879xJ/LBGAzvXrclmblnSoV5tWsbWICDZXK8jhdLL5+EkW7d7Pkj0JpOblExoYwKAWTRnetiX9mzchJPCC8sv8EmU/jkq/GiQSqTEHsej6E95GG1wJtMGVHVfL7QXI/xzC70AiHtDm5qaw2M5v+w/y0869rDpwhGKnk2Y1YxjVvjUj27emQfVqRkv0GA6nk/hjiSzcvZ9f9iSQnl9AWFAgl7RoyvB2rejXrBHBAdrszqBsm1HpN0JQF6T6p7qf2stogyuBNriyo3L/z9XvFnYzEvlYlTc3h9PJ+iPH+WnnXpbsOUCezUatiHBGtm/FFe3b0KZOLb+/Rnank41HT/zZssssKCQ8KIjBrVwtu75NGxGkzQ5V8CMq6yEIvQqJetHvvxdGog2uBNrgyobKn+0KTYaMRqq97MoWq4Iopdh9KoUfd+xlwa59nM7NIyI4iEtbt+CKDq3p0ag+VkvVvDbFDgfrj7jMbuneBLIKi4gMDubSNs2Z3LsbTWvGGC3RUJw5b0Pe+0jko0j4rUbL8Vu0wZVAG1zpqMLFqMz7IHgAEv1+lcyWzCksYvambczbtodDaekEWiwMaNGEUe1bM7BFU93/dBbFDgdrDx9n0e79LN69n0K7nSs6tOYf/XrRMCbaaHmGoJQTlfUAFC52/R+FDDFakl+iDa4E2uD+HmXbiEqfBIEdkJjPqtwgbofTyZwtO3h7xVoy8gvo3rAeozq0ZliblkSHhhgtzxSk5eXzyR/xfBm/FYdTcWOPTvyjXy8iQ4KNllbpKFWASr8B7AeRGt8iAc2NluR3aIMrgTa486Psx1BpV4El2p0BFm20pEpl7eFj/PuX39mfkkr3hvV47NIBtK9b22hZpiUlJ5d3fl/Lt1t2EhMexgOD+nBlp3ZYqlh/lHIko9LGgEQjNb5DLL4xPMRf0AZXAm1w58Y1hudacJx2mVsVGsNzPCOTl5euYum+A9SLjuKRwf0Y1qaFTgzwEDuTknlh8XK2nEiiXd1Ynho26H8Gofs7qugPVMbN7iIJr+rvlgfRBlcCbXD/i1J2VMbtYFvrSmuuIlUYcotsfLR6A5+t30ygxcLtfXtwc68uOuXdCyil+HnnPl5dtorknFxGtW/Nw4P7UTsqwmhplYbKfQ+V+w4S9S9dzsuDeNvg9N3A5Kicl8C2ylWhpAqYm1Mpfti+mzd+W83p3HzGdGzDg5f0pXZk1bnZVjYiwqgOrbmkVVOmrdnI9LWb+HXfAW7v04NbenetGj8qwu8C22ZU9r9cfdyBbUvfR2M4ugVnYv4cDhA2CUvUE0bL8Tqbjify7yW/szMpmU716vLksIF0rFfHaFlVjuMZmby0dCW/7jtI/egoHhs6gCGtmvl96E4501Gpo0GCkRrzEEuk0ZJMjw5RlkAb3F+oojWojMkQ3BeJ/tCv60uezMrmtWWrWbBrH7UjI3hocF9Gtm9d5RIefI0/Dh3l37/8TsLpNHo3acCTlw6kRWxNo2V5FWXb5MqsDB6MRL/r96bubbTBlUAbnAtlP4RKuwasdZCYrxGLf4bnCoqL+eSPeD75Ix6F4tbe3bjt4u66aLAPYXc6mR2/jXd+X0tekY3bLu7OPQN7E+DHA+hV3qeonJddtV3DJxktx9ToPjjNf6GcmaiMO0ACXC03PzW35fsPMXXRMk5l5zKibUseGtyPetFRRsvSnEWAxcKNPTpzebtWvLpsFR+u2cDWxCTeGDeCGj4y44LHCbsFbPGonFcgsBMS1MloRZrzoFtwJkKpYlTGrWDbhMTMQoK6Gi3J49gcDl5ftpoZ6zfTKrYmzwwfRLeG9Y2WpSkj323dxXOLlhEdGso7V42kU/26RkvyCsqZhUobC8qB1PxBzzxQQbzdgvPfOIKfoZRyTchoW4dUe9Evze14RhbXzfiGGes3c2P3Tnx76wRtbibjyk7t+HrSeAKtFm6YOYcvNm7FTD+iy4pYqiHRb4MzFZX5MEo5jZakOQfa4MxCwddQMAfCb0dCxxitxuOsO3Kcqz79iiNpmbx71UieumyQrmxvUtrWjeX7ydfTp1kj/rV4Oc8t+g2H0/8MQAI7IFFPgG0l5M80Wo7mHGiDMwGqeBcq+wUI6odEPGC0HI/z9abt3Prl99QID+P7yddxaZsWRkvSXCDVQkP44NrRTO7djdmbtnPvtz9TWGw3WpbnCb0Oggeict9GOU4ZrUZzFtrgfBzlzEZl3guWGkj0q3419Y3d6eSFJct5duEyLm7akG9uHl9lq9f7IxYRHh7Sj6eGDWTZvoNM+uJbMvILjJblUUQEiXwKlN2VdKLxKfznbumHuPrdngZHEhL9FmLxnzm6sgsLuf3rH/h8w1Ym9ezCh9eOrpIV66sCN/bozNtXjWRXUgoTZnzD8YwsoyV5FAloCOG3QeHPKNtGo+VoSqANzpcpXACFi5CIe5CgLkar8RhH0zO5dvo3rDt8nBdGDuHxSwdU2YlHqwrD2rTgsxuuJD0vn/Gffc2upGSjJXkUiZgCljhU9vMo5YehWJOi7yo+inKcRmU/D4EdXb8O/YR1R45zzfTZpOfnM/36cVzduYPRkjSVRLeG9Zg96VqCAqzcMGsuqw8eNVqSxxAJRaIeB/s+yJ9ttByNG21wPsifoUmVj1R7GRH/yCb8ZvNfySRzb5lAz8YNjJakqWSa1arB1zePp0F0NW7/+gd+2LbbaEmeI/hSCOrtSjhxphutRoM2ON+k8Aco+g2JfAAJaGa0mgvG7nTy4pIVPLNgGb2b6GSSqk7tyAi+nHgN3RvV49Efl/Dh6g1+MVZORJAo1w9TlfO60XI0lMHgRGS6iKSIyM4S664WkV0i4hSRc45CF5FWIrK1xCNbRO53vzdVRBJLvDfCY5/I5CjHKdeQgMAuEDbJaDkXTE5hEbd//QOzNmxhYs/OfDheJ5NoIDIkmGkTxjKqfWveXL7Gb8bKSUBzCLsRCr5FFe8wWk6VpyyxrxnAe8CsEut2AuOAj863k1JqH9AJQFyl7hOBeSU2eVMp9Vr55Po3SilU1pOA3R2aNPcMAUfTM7nj6/kcy8jkhZFDdH+b5r8Islp5Zcxl1I6M4JO18aTk5PHWlSNMP8BfIu5BFf7k6kOP+cavhvaYjVKvvFJqJZB+1ro9bgMrK4OBg0op/+lV9gYFc12Tl0Y8hAQ0MlrNBbE7KUUnk2hK5b/Gyu0/yMM/LDZ9S04sEUjkw1C8DQrmlb6DxmtU1k+L8cDZqUV3i8h2dwj0vJVKRWSKiMSLSPzp06e9q9JAlOO0a3buoJ4Qdr3Rci6IE5lZTJ49j9DAQOboZBJNGbixR2ceGdKPxXsSeH/lOqPlXDghoyGwMyr3DZTyr8HtZsLrBiciQcAVwNwSqz8AmuEKYSYB5+2RVUpNU0p1U0p1q1WrljelGorKeRlUERL1vKlDGtmFhdw++wdsdgefXjeWRjqZRFNGbunVlbEd2/L+qvX8nnDYaDkXhIgFiXwInKch/0uj5VRZKuNOOhzYrJT6c2SnUipZKeVQrhLcHwM9KkGHz6JsG6DwRwifjAQ0MVpOhSl2OLjv2wUcSc/k3atH0qxWDaMlaUyEiPDsiEtoXbsWD/+wyPQVTySoOwT1Q+VOQzlzjZZTJakMg5vAWeFJESk5SdRYXEkrVRKlil3T4FjqIRF3GC2nwiilmLrwN/44fIznLx9C7yYNjZakMSGhgYG8e9VInAru+/ZniuzmrgoikfeDyoT8GQYrqZqUZZjAbGAt0EpETojIrSIyVkROAL2BBSKyxL1tnIgsLLFvGDAU+P6sw74iIjtEZDswCPC/EvllJX8W2BOQqCcRCTVaTYX5+I94vt26kzv69uDKTu2MlqMxMQ1jonllzDB2nUrhX4uXGy3ngpDADhA8BJU3HeXMNFpOlaPUfFyl1ITzvPU/6UFKqZPAiBKv84H/iVMppW4sh0a/RTlOoXLfg+CBEDzYaDkVZtHu/bz+22oub9eK+wZebLQcjR9wSctm3NGnBx+u2UCn+nW5qlN7oyVVGIm4H5U2CpX3iatfTlNpmDebwQ9QOf8BZUcin0JEjJZTIbaeSOLR+YvpXL8u/7niUiwm/Rwa3+Pegb3p3aQBzy38jd1JKUbLqTAS2BJCRkL+5yiH/2aC+yLa4AxCFa1xzxRwu2u6DRNyPCOLO7+ZT+3ICP7vmisINvkAXY1vYbVYeGPsCGLCQ7nn25/IKig0WlKFkYh7QNlQeeetjaHxAvqOZABK2VxVDqwNTTtTQFZBIVNm/4DD6eSj8WOICQ8zWpLhFNvsJOxKZGf8EfLzChEEEUBcS3Ev4a/X53ovLCKYDt2b0LBZrGlb9p4iJjyMt68cyQ0z5/DI/MV8cO1oU0YJJKAxKnQs5M9Ghd+KWOuWvpPmgtEGZwR508FxGKn+MSLmq8toczi499ufOZ6RyafXj6NpTf+ZiLU82Isd7N95gu0bDrF9wyF2bzlKUUExABarBZRCKSpcSDimViSdejWj08XN6dSrObXqVPOkfNPQqX5dHr90AM8vXs6HqzdwV7+eRkuqEBLxD1TBfFTu+0i1F4yWUyXQBlfJKMdJVO7/QfBQJHiA0XLKjWs4wDLWHTnOy1cMq1JVSpRSHN53io0r97F9wyF2bT7yp6E1blmHYVd2o2OPZrTv1phq1cPPub9ym15J81P89+vMtFy2rTvIlrUH2LQmgd9+2gpA/Sa16Ny7Gd36taLTxc0JCqo6/77XdbuILSeSeGfFH3SMq0PfZuYrZSfWeqiwa92tuNtMX47PDIiZpqno1q2bio+PN1rGBeHMfAgKlyC1liDWOKPllJtvNm/nmQXLuKtfzyqRMVlUWMy29QfZsGIvG37fy+kk1+Djxi3r0LFHEzp2b0r77k3OaWiewOl0cjQhmS1rD7B17UF2xB+mMN9GWEQwPQa2ps/Q9nTr25KQsCCvnN+XyLcVc8302WQWFLLwzpuICgkxWlK5UY4U1OnBEDoGS7V/GS3HcERkk1LqnDPSeOT42uAqD1W8F5U2GsInY4l82Gg55WZv8mmumT6b7g3r8/F1Y03ZF1IW8nIKWbloOxt+38uWtQcoKigmJCyILhe3oMfA1nTv34qYWpGGaLPZ7Gxbd5A1v+xk7bLdZGfmExwSSLd+LelzaXt6DGhNeKT5bvxlZWdSMtd8OpuxF7XjxVFDjZZTIZxZT0HBfCT2d8RSNcP7Z9AGVwKzG5wz43awbUJqLUMs5upPybPZuPKTr8gtsjF/yg3U8MOkEqfTydJ5m5nx5hIy03KJrRtNz0Gt6TmoDR26NyEoONBoif+Fw+5gR/wR1vyykz9+3UX66RwCAq1cOq4bt/zzMr81uteWreLjP+L57PpxXNzUfGE+ZT+ASh2BRNyPRNxltBxD0QZXAjMbnLLFo9KvQyL+iUTcbrSccqGU4pH5i/l55z5m3HClX/a77dl6jA9e/JGEnYm07dyIyY+MoPVFDUyTxeh0Otm77TjL5m9h8dwNVK8ZyV1PX8HFQ/yvqkxhsZ3R076g2Ongp9tvJDzIfOFZZ/qtYN+D1FqBqx591cTbBqfHwVUCSilUzmtgiYXwm4yWU26+37abH3fs5R/9e/mduaWfzuH1x+fy4IQPSEvO5uFXruG1L2+nTaeGpjE3AIvFQtvOjbhn6hje+PpOoqqH8a97vuCF+74kPSXbaHkeJSQwgBdHDeVkZjZvLl9jtJwKIeE3gzMVCn82Wopfow2uMihaAcWbkYi7TVdvMiEllecX/Uavxg24s6//TPpQbLPz7fSVTB7+OisWbOOa2wbw8aJ/csmozqYytnPRqkMD3pl7N5MeGMaGFXuZMvJNFs/dWOHhCr5It4b1uL57J77YsJVNxxONllN+gvpAQAtU3md+9XfxNXSI0sso5XAllqgipOZCRHyrH+fvyLcVc/WnX5FRUMj8KTdQK8I7mYKVTfyqfXz0nwWcOHyaHgNaMeWxkdRrXNNoWV7hxOHTvPPsPHZsPEyH7k2497mx1G/iH/Mq5hbZGPXRLIIDApg/5QbTVdJR+XNR2U8i1Wciwb2NlmMIOkRpdgoXg30/EnGfqcwN4OVfV3IwNZ3XxlzmF+Z28lgaz/1jFk9PmYFyOnnug4k89+EkvzU3cI2de2nGZO57fhyH9iZx15h3+GbaCuzFDqOlXTARwUH86/KhHE7LMOcs4KFXgCUGlf+Z0Ur8Fm1wXkQphcr7EKxNIWRE6Tv4EPHHTvD1pu1M6tXFlJlqJSnIK2LGW0u4feSbbFt3kFv+eRn/9+P99BjY2mhplYLFYuGyq7sz7ecH6DGwNTPeXMK9V7/Hvh3HjZZ2wfRt1ohxF7Xlkz/i2ZWUXPoOPoRIMIReB0UrUPZDRsvxS7TBeZOiFWDfh0RMQcQ8l9pmt/PMgmXUqxbFvQPMO5hbKcWKBVu57fI3+OajFQwY3pGPF/2TqycPqFJVQM4QExvFU29fz9Pv3kB2Rj4Pjv+Aj19eQGG+zWhpF8RjQwcQEx7K1IW/ma4/S8KuAwJQ+XOMluKXmOeuazJcrbcPwFofQkYZLadcTPsjnoOp6Tw74hLCgswVVj3DicOneeTGabz80DdUrxHBa1/ezkMvX0ON2CijpRnOxUPa8dHPD3DZ1T34fsZq7rjiLXZsPGy0rApTLTSEBwb1ZfvJUyzavd9oOeVCrDUheAgUfI9S5v6h4Ytog/MWxduheCsSdrOp+t5O5+bxyR8buaxNCwY0b2K0nAqxd9sx7r3qPY4eSObe58by1px/0K5LY6Nl+RThkSHcM3UMr34+hYAAC4/f8gk/f7XWdC2gM4zp2IZmNWP4aI35skUl7EpQmVC02mgpfoc2OC+hCmaDhEHoWKOllIsPV2+g2OHkwUv6Gi2lQhzel8TTt88gumYE7/9wH8Ov6YHVqr/m56N9tya8PfduuvZpwfv/+pF3p/6Aw+E0Wla5sVos3NyrC3uTT7P+6Amj5ZSPoItBqqEKFxutxO/Q//leQDmzoGABhIxCLBFGyykzxzOy+GbTdq7s1I5GMdFGyyk3J4+m8uTk6QSHBPKf6ZOr7PQy5SU8MoRn3r+Ja24bwKI5G5jx5hKjJVWIUe3bUD0slBnrNhstpVyIBELIYChapsOUHkYbnDcomAcUIWETjFZSLt5fuQ6LRfhHv15GSyk3p5MyefyWT3HYnfz701upXa+60ZJMhdVq4eYHL+Py8T359tOVLP1hk9GSyk1IYAATunZkRcIhjqRlGC2nXEjwMFA5YDPhcAcfRhuch1FKofJnQ+BFSGBbo+WUmQOn05i/Yw/Xd+9E7SjztDoBMtNzeeLW6eRmF/DCJ7fQsFms0ZJMyx1PjOKiXs1455l57N5y1Gg55ea6bhcRYLUya8MWo6WUj+A+IOE6TOlhtMF5GtsG12zdJmu9vb3iD0IDA7nt4u5GSykXeTmFPDX5M1JOZvDcBxNp0a6e0ZJMTUCglSffvI7YutH8654vSE40V0uoVkQ4I9u14vttu8gqKDRaTpkRCYLgwVD4K0oVGy3Hb9AG52FUwVcg1Uw1sHv7yVP8svcAt/TuSkyYeWplFhbYePbOmRxJOMVT79xA+27mzPr0NSKjw5j6wU0U2+w8949ZFOQVGS2pXEzq2YWCYjtztuwwWkq5kJBhrmxK2wajpfgNpRqciEwXkRQR2Vli3dUisktEnCJy3jpiInJERHaIyFYRiS+xPkZElopIgnvpFx0mypEKhUshdCwi5pmL663la6geFsqknl2MllJmim12Xrj3C/ZsOcojr1xL9/6tjJbkVzRoGsvjb07gaEIyrz46B6fTPJmVrevUomfjBnyxcSvFDhOVJAvuBxKGKjRnko8vUpZyDjOA94BZJdbtBMYBH5Vh/0FKqdSz1j0GLFNKvSQij7lfP1qGY/k2BXMBOxI23mglZWbdkeOsOXSMx4cOICLYHPNSORxOXnnkGzatTuD+f42j//CORksCoKigiPULNrPq+/XkZeUTFBJIYHCgaxkU+OfrwOAAgkKC/nrvv9YFEFMnmlbdmxs+q0HXPi2Z8tjlfPjvn5n1zlIm3T/MUD3lYVLPztz5zY8s3XuAEe3M8eNHJAQVPACKlqLUs4hYjZZkeko1OKXUShFpfNa6PcCF/AOOBga6n88EVmByg1PK4Sq3E9QLCWhqtJwyoZTizd/WUCcqggndfMMkSsPpdPLOM/NYvWQntz06gmFXGdtn6LA72Lp8J8u+WsWa7zeQn1NA9drVqNWgJrZCG8VF9j+XxUXF2AqLKS4qvY+lUdv6jLpzGENvGkBYpHFh4ytuuJgjCcl889EKGjaL5ZJRnQ3TUh4GtmhK45hoZqzfbBqDA5CQy1CFi6B4EwT5z/RURuHtgnwK+EVEFPCRUmqae31tpVQSgFIqSUTMn/ZmWwXORCTsEaOVlJnlCYfYmpjECyOHmGaqkc/eWMIv38cz4c5LGDepn2E6Uk+mM+/thSydtYKM5CzCokLpf1UvBl3Xj4sGtsVqPf+vb6UUxTaX4ZU0vTPLQ9uP8dMHS3jvnk/59PEvGXHbEK5/6koiq1d+dquIcNdTV5B4JJW3nvqe+k1q0bJ9/UrXUV4sItzUozPPL17OlhMn6Vw/zmhJZSOoPxCCKlyMaIO7YLx9V+ujlDrpNrClIrJXKbWyPAcQkSnAFICGDRt6Q6NHUAXzQaJddeVMwmfrNlMvOoqxF7UzWkqZWLtsN99+upIR1/bkxnuMuc7H9yUy59UfWfbFShx2B71Hd2fw9f3pOaIzQSFlC/GKCEHBgQQFn7uEW8uuzbjs5kHs3ZDAD+8tYt7bC1g663cm/Ws8IyYPxhpQuaGrwKAAnnr7eu4a+w5vP/0978z9R6VrqAhjLmrLq8tW8cP2PaYxOLGEo4IvhqJy3SY158GrWZRKqZPuZQowDzjzkyRZROoCuJcpf3OMaUqpbkqpbrVq+eZEjUoVQNFyCBlmmrqTe5NPs+HoCa7vdhEBFt9Ppj2dlMmbT35LszZx3P7EyErvn9qzPoGpV77KrW0f4LevVnHZrYOZsf9dpn73MP3G9SyzuZWH1j1a8Nise/m/Ta/QpEND3rnrY+7s+ghbfqv87MCo6uHc8cQoDu1N4scv1lb6+StCeFAQg1o05Zc9CdhNlCQjQb3BcQzlMOFM5T6G1+5sIhIuIpFnngOX4kpOAfgRmOh+PhGY7y0dlULRClD5iImGBnwVv43gACtXdmpvtJRScdgdvPLwN9iLHTz+xoRKm+pGKcXGxVt46JKp3Nv7CbYt38WEx8fyxZEPuPf9ydRtWrtSdDS7qDGvLnuWZ759iIKcAh4Z8jxTr3yVpEOVO/9Zn6Ht6N6/FZ+/u5TTp7Iq9dwVZXjblqTnF7DRTPUpg9yVhIrWG6vDDyjLMIHZwFqglYicEJFbRWSsiJwAegMLRGSJe9s4EVno3rU2sFpEtgEbgAVKqTPD9F8ChopIAjDU/dq0qIKFYKlpmk7hrIJCftyxh1Ht2xAd6vvDGb764Dd2bjrC3c+OqZTZtx12B7/NXs0dXR7miRH/JjEhidtfu4kvj37AzS9MoHps5de4FBH6jevJp7vf4uYXJrDpl23c2vZ+Pn3iK/JzCipNw51PXYHD4eSj//xcKee8UPo3b0JYUCALd+0zWkrZCWjhmunb9ofRSkxPWbIoz1eSY945tj0JjHA/PwRcdJ5jpgGDyy7Td1HOXFcLLuwq06T1fr9tFwXFdm7ofs4/j0+xfcMhZn+wnMGjO3PJFd7P4CvML+L5q19n46ItNGhdj39+eheDr+9LoI/MixcUEsR1T4zj0kkDmf7EV3z90jx+mbGcW/9zPUNu7I/Fy+Hmug1imHDnJcx86xc2rNjr87OihwQGcEnLZvyy9wDPDL+EwL9J/vEVRCyooF5gW4dSyvDhImbG9ztffJ2i5UCRacKTTqX4Mn4bXRvE0aaObyevZmXk8crD3xDXMIZ/PD3a6+fLy87niREvEr94K3e/eyuf7HyDy24e5DPmVpKacTE8MuNu3ln7b2Ib1eLVm9/n3t5PsHut91sqV97cjwbNYvm/F36ksMD3q98Pb9uSzIJC1h8xT5hSgnqBMwUc5p2I1hfQBneBqMKFYKkNgV2NllImVh44zPGMLG7o3sloKX+LUoo3Hp9LVkYej71xHaHhwV49X3ZaDo8OfZ7df+zn8S/vY/Q/LvN6a8gTtOnZgrfXvMCjs+4hNTGd+/o8xYxnvvbqOQODArj7mdEkJ2bw9YfLvXouT9CvWSPCg4LMNdt30MWupc0cCT2+iu//B/swypntSucNGY6IOS7lFxu3USsinKGtmxst5W/5+at1bPh9H5MfGUHztt5N8U5LyuCfg57l0PZjPPvdQwwa38er5/M0FouFITf057O9bzNs0iC+fOE7vn97gVfP2bFHUwaP7sx3n63i6IHKTXYpL8EBAQxu1ZSlexOwmaV0l7UBWOJQRXr6nAvBHHdlX6XoV6DYNOHJI2kZrDp4hPFdO/p0X0ROZj6z3l1K597NueL63l49V/LR0zw44BlOHU7hxQWP03vUeUur+jyhEaE88PHt9B3Xkw8emMHyr9d49XyTHx5BSFgQ7z03H6WUV891oYxo24qswiLWHj5mtJQyISIQ3Ats61HKPEMcfA1tcBeAK3uyHgT6frIGwJfx2wi0WLi2Swejpfwtsz/8jfycQm577HKvdrCf2H+SB/o/TXZqDi8vfYbOl/j2dSkLVquVx7+4lw792/DKxHfZ/Ot2r50rukYEtzx4GTvjD/PrfN+eRbtP04ZEBgez2ERhSgnq7ZpdwL7XaCmmRRtcBVHOTLD9AaHDTZHlVFBczPfbdjGsbUtqRYQbLee8JCdm8NNX6xg6titNWtbx2nlOHUnhgf7PUFxYzGvLp9K2V0uvnauyCQoJ4vkfHqVB63pMHfcqJ/af9Nq5hl3VjTadGvLJK4vIy/Xd+deCAgIY0qoZS/ceNE+Y8sx4ON0PV2G0wVUU21rAjpikNNfqg0fJLbJxVSffLsu1aO4GnA4n193lvVEk9mI7/77uLWyFNl7//XmaXdTYa+cyiojocP698AnEInzw4AyvncdisXDbIyPIzshjxc/bvHYeTzCwRRNyiorYc+q8hZN8CrHWBmt9VLH3WuH+jja4CqKK1oBEQKA5qvAv3XuA6NAQujfy3UK5DruDX+dtpmvflsTGRXvtPLOmzmHPugQe+Oh2Grb23xnAa9arwY3PXM2GhVtYv9B7IcTWnRrStHVdFs5Z79N9cV0auJKVNh/3XovW4wS0gWIdoqwo2uAqiu0PCOqJiO9X4S92OFiecIhBLZr6dN3J+FX7SUvJZtiV3kv02PLbDr5+6Qcuu+USBl5rrmzJijD67sto0CqOj/45k2Jb6dP0VAQRYfg1PTi0J4n9O313rFlsZAQNqlczlcFJYBtwHEE584yWYkp8927nwyj7MXCcQILMcYPcePQE2YVFDPHxoQFLvounWky416pjZKVm89KN71K/ZV3uevtmr5zD1wgMCuSONyZxfN9J5r+3uPQdKsigUZ0ICQvil+/ivXYOT9Clfhybjp/06ZbmfxHQFlBgN1GpMR9CG1xFOFMjLvhiY3WUkV/3HSQ0MIC+TRsZLeW8pJ/OYcPvexkypguBXiimrJTitVv+j5y0HJ6YfT+h4b5fg9NT9BjemR4jOvP583PJSPFOkeTwiBB69G/F2mV7cPpw5f6uDeNIy8vnWIY5ikUT2Na1LN5trA6Tog2uAqiiNWCpC9YmRkspFadS/LrvIH2bNiYk0HfDqcvmb8Zhd3otPDn/vcWs+3kTt71yI807+f7fzdPc8fpEivJtzHhqttfO0fOSNmSk5rB/h++GKbvUP9MPZ5KpaCy1QaJRdm1wFUEbXDlRygG2dRB8sSmGB+w8mUxyTi5DWjczWsp5UUqx5Lt42nZpRIOmnq+PeXDbEaY9PIuel3dhzD3DPX58M9CgVT3G3DOcRZ/+xoEt3qlv2L1fKyxWC+uW7/HK8T1Bs1o1iAoJZpNJ+uFExNWK04kmFUIbXHmx7wKVZZr+t1/3HcAqwqAWTY2Wcl52bT5K4pFUr7TeCvIKeXHCW0TWiOSh6XeZ4keJt7jh6auoVjOS/7v/M6/0QUVGh9GhW2PW/+a7BmcRoUuDOFMlmhDQBuz7UMo7SUL+jDa48lLkLn8U5N0SUp5i6d6D9Ghcn2o+PO/bkm83EhoWRL9hnq8k8sH9Mzix7ySPfX4P0bUqfx43XyIiOpybX5jAjlV7WDnXO4OHew5qw5GEZJKOp3vl+J6gS/04Dqamk5FfOfPoXSgS2BYoBvsho6WYDm1w5UTZ1kJAa8Raw2gppXI4LYNDaekMaeW72ZN5uYWsWrKDAZdf5PEZA1Z8s4ZFny5j/GNj/KIMlycYdssgmnduwrRHPqcwv8jjx+85qA0A637z3T6jM+PhtpwwSSsu0HVN0f1w5UYbXDlQqhBsm8Ak4ck1h44CMKB5Y2OF/A2rl+ykqKDY4+HJ7PQc3rpjGm16teCmqdd49Nhmxmq1ctdbN5NyLJUf3/f8sIG4hjVo1Lw2G1b4bp9Rh7g6BFosbDmeZLSUsmFtAoSgin039OuraIMrD8X7gGIkyPszS3uC3adSqBEeRoPq0UZLOS87Nh4iukYErTo28Ohxl89eQ15WPve8N5kAH84eNYIO/drQslsz1v7knTFrrS9q4NNT6IQEBtAwJpqj6RlGSykTIlYIaAiO40ZLMR3a4MrDmRBBoG/XczzD/pRUWtTy7VDq3m3HadOpoceTP36ZuYKmFzWiRRffTa4xks6XtGfPugQK8jxfIDmuUU0yUnN9uvhyXLUoTmRlGy2j7FjrgsMkIVUfQhtcOVDFu0CqgcW7E3B6AqdSHDidRsvYmkZLOS/ZGXkkHkml9UWebb2lHDvN/viDDLm+v0eP6090HtwBh93BztWeDyXGNXL9qEo6mubxY3uKetFRnMzKMVpG2bHUA4dJQqo+hDa48lC8GwLbmiLVPDEzi4JiOy1jfbcFt889ILj1RQ09etxtK1wt7a6XmmOePiNo16c1AYFWti7b4fFj13Mb3MljPmxw1aLIyC8gz2YzWkqZEGtdUJm6JmU50QZXRpQqdtWDC2hrtJQysS/FdXNpUct3W3B7tx3DYhFatPNsRf9tK3YRGRNB4/aebRn6EyFhwbTp3ZIty3d6/Nh1G7gMLtGHW3D1o6MAOJlpkjCl1R01cupWXHnQBldW7AeBYsQk/W8JKakANPfhPri9247TuGUdjw8P2Pb7LjoOaIvFh2dO8AU6X9KBA5sPk53u2VBdSFgQNWKjOHk01aPH9SRx1VwGl2iWfrgzBqf74cqFvgOUlTPFTs+MSfFxEk6nUS86iojgIKOlnBOn08m+7cc93v+WfPQ0pw6ncNEAc/wQMZLOl7RHKcX23z0/viquUQ1O+nALrp67BZdothacNrhyUarBich0EUkRkZ0l1l0tIrtExCki5xzAJCINRGS5iOxxb3tfifemikiiiGx1P0Z45uN4D2XfDRIG1sZGSykT+1NSaenDrbcTh1PJyyn0eP/bmZv1RQO1wZVGqx7NCQkPZutvng9T1mtU06f74GqGhxEcYDVPJqWlFmBF6USTclGWFtwM4LKz1u0ExgEr/2Y/O/BPpVQboBfwDxEp2YH1plKqk/uxsByajaF4l6uCiViNVlIqNoeDw2kZPp1BuXfbMQCPt+B0/1vZCQwKpEO/Nmz5zfOJJnGNapCZ5rtDBUSEuGpRpumDEwkASx1wmGQWBB+hVINTSq0E0s9at0cp9bcz8CmlkpRSm93Pc4A9gGezCSoJpZxg32ua8OSRtAzsTqdPJ5js2XqMiKgQ6jX2rEbd/1Y+Og1qz7E9iaQleXbQc1xDdyalL4cpq0WZpw8O3GPhdAuuPFTKXUBEGgOdgfUlVt8tItvdIdDqlaGjwjiOg8pDTJJBeeD0mQxK3w1RJuxMpEX7+h41orSkDE4dTqFjP3P8nXyBM6HcXWs8Ox6urtvgTvlw0eW4apHmGgtnrauzKMuJ1w1ORCKA74D7lVJnfi59ADQDOgFJwOt/s/8UEYkXkfjTp097W+65cbgncAzw3RmxS5KWlw9AbFSEwUrOT05WPjG1Ij16zIIcV3X4arWiPHpcfyYoJNArxw0IdIXyHQ7fnd07NCiQIrvdaBllxxIFzlyjVZgKrxqciATiMrcvlVLfn1mvlEpWSjmUUk7gY6DH+Y6hlJqmlOqmlOpWq1Ytb8o9P2fCApY6xpy/nGQWuPo9qoV4Nv3ekxQVFhMS6tkMzzM3U4tVhyfLSl6260dBWFSYR4/rsDuAv4zOF7FaLDicvmvA/4OEgdIDvcuD1+4E4ir38SmwRyn1xlnv1S3xciyupBXf5UxYwGoeg4sKCcbqw/1QhQU2gj1scMrpmsTTYvH9SjO+Ql6Wq7UfXs2zBme3u4wjIMB3v4MBYsHu9PzEr95CJBwo1hOfloOyDBOYDawFWonICRG5VUTGisgJoDewQESWuLeNE5EzGZF9gBuBS84xHOAVEdkhItuBQcADnv5gnkQ5ksBSCxHfHFN2NpkFhUT78ASnTqeTooJiQkI9Gx5z6hZcufGawRW7WnDWAN9uwTnN1oIDUOaYqNUXKHUeEaXUhPO8Ne8c254ERrifrwbO+VNaKXVjOTQaj+OUaVpvAJkFBUSHhhot47zYilz9HsEhnv3BcOZmJboFV2by/wxRevb78meI0ocNLsAiOJRCKWWK+rJ/GVw+oPuZy4L+qVsWnMlgqW20ijKTmV9INR9uwRUWuArchoTpFpzReDtEafXhEOWZDF6HMkmY8k+D0/1wZcV3v32+hCMZrOYxuKxC3w5RFhW4+hA8nWTyVx+c/lqXlbysPCxWCyFhnk1IMkOSSYC7pW+aRJP/asFpyoK+E5SCUgWgshGTteCiw3zX4M604DwdotRZlOUnP7uA8GphHg/RnemD82WDO5OEZTebwTm1wZUVfScoDUeKa2mSFpzd6SSnqMin++CKCs+04LwTotR9cGUnLzufcA/3v0HJLErfNzjTtOAsugVXXrTBlYYz2bU0SQsu2z0GzpdDlIX57hacl0KUVt2CKzN5WfmEebj/Df4KUfpyH1zAny04s/XBaYMrK7777fMVHO7qKRaDBpmXk+zCIgAig313kPef4SsP3/zOtNyKi/Q4obJy+ngaUTGer3hzJgwdGFRqorbxmCbJxN3S1sMEyow2uNI482WyhBuro4ycCbs4le+GXWrWqQZASlKmR4/bpENDLBZh74YDHj2uv3JsbyIHthymx4iunj/2wRQCgwKoHRft8WN7ivS8fASI8uFox3+h3GXFxAQ/GnwEbXClUuRe+m6LqCTB7j6PIneIyBep0yAGi0U8Xmk+PCqMZp0as2PVHo8e119ZOnMFFquFwdf39fixD+1JonGL2j490Pt0bh4x4WF/hip9HmVzPzFHwQlfwCR/WQNRboMTsxic69edLxeRDQoKoFbdaBKPpnr82B36tWXPuv0U23SY8u9wOBws/fx3egzvTEwdz07moZTi0L4kmrSqW/rGBnI6N59aEeaIzLhwG5xJKir5AtrgSuNPgzPHlyrIBC04cE2I6Y25wtr3a4OtsJj98Yc8fmx/YvOvO0g7mcHQmwZ4/NgZp3PISs+jaWvfrv5zOjfPXAZ3pgVnkh/bvoA2uFJQf4YFvDOtiKcJsroMzubwbYOr16gmJ4+mojzcwd+hX2sAHaYshaWzVhBZPZxeo7p5/NgH97qKkzdtHefxY3sSl8F5PoPUa5z5sa1DlGVGG1ypFAHB5qhVhyvJJNBiwebDIUpwteByswvJzvRsynN0rWo0aF2PHat2e/S4/kRuZh5r5m1g0IS+BAV7/ofb4X0ug2vSyndbcE6lSMvTIUp/Rxtcaagi04UEggICfD9E6Z7x+eQRz/fDdezXhl1r9uHw8VasUfw+5w9shcVcOmmQV45/aG8SsXHRRHhhALmnyMwvwO50msvglDa48qINrjRMaHDBAVafTjIBqNe4JgCJXuqHy8vK5/COYx4/tj/wy6zfadS2Pi27NvXK8Q/tTaJpa99OMEnJdRUsrmlKgzPX/chItMGVhikNzvdbcLXrVXcNFTjmeYPr2L8NAFuW+fY8ukZwYv9Jdv+xj0snDvRK2L2osJjEI6k+b3Cpua7QuKlacOhhAuVFG1yp2DDbFyo4wOrzfXCBQQHE1qvOSS8MFYhtWIsWXZvy6+e/ezyJxez8MnMFFosw+Ib+Xjn+kYRTOJ2Kpj4/RCAXgFgzGZzJMrp9AW1wpaGKTVc5IDQwkFybrfQNDaZh01oc2JXoFRMafutgDm0/ypIZKzx+bLOSejKdnz/8he7DO1OjrmfHvp1h96ajADRvV88rx/cUiVk5CFAr0kwG56oza7aIkpFogysNCfnri2US6kVHcSIjy2gZpdJzUBsSj6ZxeN8pjx97xOTBdLqkPe/c9TH7Nx30+PHNhlKK12/9P2yFxdzx+kSvnWfFwm00bxtH7XreMVBPsfNkMk1rxhAaaI7hPwDKmQIEgkQbLcU0aIMrDQk1XXHThtWjOZ6ZhdPHw3MXD22HxWph5eLtHj+2NcDKk7PvJzo2iuevep2s1GyPn8NM/PzhL8Qv2cZtr9xI/ZbeGZ928mgq+3ecYMDlF3nl+J5CKcWOk6foEGeOGUL+xD3xslmGLPkC2uBKQ8JMNz1Fg+rVKLI7SMnJNVrK3xIdE8FFPZuyatEOr4Qpo2tV49lvHyI9KYN/X/dWlR02cCIhiWkPf07XSy/iiruGee08KxZsQ0QYMMK3De5Udi6pefl0jPPdcXrnxJkCllijVZgKbXCl4W7BmSlZoVFMNADHTBCm7H9ZR04eS+Pg7pNeOX6r7s255/3JbP51BzOe/sYr5/BlHHYHr0x8l8DgAB769E6v/fp3Op389uMW2nVtTC33bBG+yvaTrpB4h3omMzh3C05TdrTBlYJIKODkrxRd36dhddcN5mh6prFCysDFQ9piDbCwcvEOr51j+K2DGTF5MF+/NI/V89Z77Ty+yNcv/8CedQnc8/5t1KxXw2vnWb98L4lH07h8fA+vncNTbE88RaDFQuvYmkZLKR/OFNNMvOwraIMrDRPOolu3WhQBFgvHMzKNllIqUdXD6dy7OasWb/dqK/kf795K6x7NeXXS+xzbm+i18/gSCZsP8flzcxl47cUMGt/Ha+dRSjH3k9+pU786/YZ18Np5PMWOpGRa16lFUIB5sqOVMxdUHmLVIcryoA2uNEw4i26AxUK96ChThCgB+l3WgVMnMti/84TXzhEUHMjTc/9JUEggU8e9Sn6Oef6eFcFWaOPlm94lOjaKe96f7NVz7dp0hD1bjzFuUj+fnv8NXDUod55MpoPp+t+SXUvdgisXpRqciEwXkRQR2Vli3dUisktEnCJy3nLkInKZiOwTkQMi8liJ9TEislREEtxL380pNmELDlyJJsdMEKIE6D24HQGBVlZ5MUwJENugJk9+/QCJ+0/y2i3vm6pftbxMf3I2R3ef4J+f3kVUTKRXzzX305VEVQ9n6DjPzwzuaQ6nppNns9HRjBmUoJNMyklZWnAzgMvOWrcTGAesPN9OImIF3geGA22BCSLS1v32Y8AypVQLYJn7tW9iwhYcQKPq0RzLyDLFTTyyWqgrTOmlbMqSdBrUnskv38iq79Yz59UfvXouo9i2Yhffv7WAUXcOo/uwTl4915H9p9iwYi9XXN+bkFDfr7BxJsHEfBmUboPTSSblolSDU0qtBNLPWrdHKbWvlF17AAeUUoeUa1K1r4HR7vdGAzPdz2cCY8ojulL5swVnLoNrWD2anKIiMgvMMUi9//COpCRlsnfbca+f66oHRzLgmt5Mf+JLNi7Z6vXzVSa5mXm8Muk94prX4bZXbvD6+b6dvpLg0EBGXdfL6+fyBDtOJhMeFESTmjFGSykfjhTXUrfgyoU3++DqASXvVifc6wBqK6WSANzL8/7VRGSKiMSLSPzp06e9Jva8WNzhHac5+rPO0KSmK+q7L8XztR69Qe/BbQkODWT+52u8fi4R4Z+f3EnDtvV5dvTLLJmx3OvnrAx2rtnLQ5dMJe1kBo/MvJvQ8BCvnm/X5iP89uNWRlzTg6jq5ih5te7IcS6qVweLyQZLK8cJkGjEYqIJWn0Abxrcub5B5Y4/KaWmKaW6KaW61apVywOyyonFXTTWkVT5574AOteviwAbj3ovccOThEeGMPamvvy+cDsHdnk/yzE0IpTXVzxHu76tee2W/2Paw7NMOxD8+L5Epl75Kg/0e5qM5Cyem/cwbXu19Oo5C/KKeP3xudSuF80N9wz16rk8xcHTaRxMTWdIq2ZGSyk/9n0Q4N2/qT/iTYM7ATQo8bo+cGY0b7KI1AVwL1O8qOPCsMQAwSiTGVxUSAht6sSywSQGB3DVrf2JrBbKjLeWVMr5omIi+c+iJ7nirmHMff0nnhn9MnnZ5kkmSj+Vwdt3TmNy+wfZ8usOJv1rPDP2v0PPy72f7PHpa4s4dTyDB/99NWHh5ij+u2TvAQCGtG5usJLyoZQT7PshsLXRUkyHNw1uI9BCRJqISBAwHjjTq/8jcKbi60Rgvhd1XBAiAtY4cHqn0oY36d6oHtsSk3x+6pwzhEeGcO3tg9i0OoFt6yunQHJAYAD3vDeZe9+fTPySbdx38ZOcPOj54s+epCC3gM+fm8vEFvew6NPfGHXHpcxIeJfrn7zS62FJgE1r9rPg6/WMndiHDt2beP18nuKXPQl0rl+X2pERRkspH47joPKRgFZGKzEdZRkmMBtYC7QSkRMicquIjBWRE0BvYIGILHFvGyciCwGUUnbgbmAJsAeYo5Ta5T7sS8BQEUkAhrpf+y7WOHCYpyV0hh4N61Nkd/yZOWYGRl3Xi5p1qvHZG4srNQN01J3DePmXp0k/lcndPR/njx83+lwGqsPu4OePljKxxT3Mem4O3Yd35pNdb3L3u7dSPbZyymPlZBXw5pPf0aBZLBPvv7RSzukJjmdksif5NMPatDBaSvmx73UtdQuu3JQ6lF8pNeE8b807x7YngRElXi8EFp5juzRgcNllGoy1ARQuNlpFueneqD5WEVYdPEq3hvWNllMmgoIDueHuIbz11Hf8sXQXfS5tX2nn7jSoPe+t/w9Tx73Ks2NeoXG7Bgya0JdBE/pQt4lx6dlKKdb+GM8nj3/J8b2JtO/bmqnzHvF6P9u5+PDFH8lIzeWZ924kKNg8U80s2eMKTw41WXgSQBXvBSwQYD7tRqMrmZQBsTYAlYly5hgtpVxUCw2hS4M4ViQcNlpKuRgyujMNmtbiszcWU5hfuTVA45rV4YPNr3Dv+5MJjw7js6dmc1Ozu7mv71PMf38xGSmVm027e91+HhzwDM+OfQWU4rl5j/DG788bYm5rftnJbz9tZcIdg2jZ3hw/mM6wZE8C7erGUj/atwtBnxP7PrA2ctfF1ZQHbXBlIaCha+nw/hgtTzOwRVP2Jp/mVLZ5zNkaYOXOp67g5LF0Pvz3T5V/fquVUXcO461VL/D5ofe59d/XUZBTwHv3fMr4elN4YsSLLP38d6+W+zqRkMTz17zOfRc/SWJCEvd9MIWPd7zBxaO7GzIfWEZqDu9M/YEW7eox/vZBlX7+CyEpK4ftJ08xrLUJw5PgClHq8GSFME+1USOxupNBHcchsO3fb+tjDGzRhFeXrWJFwmHGd+1otJwy07l3c66ZMoBvPlrBRb2aMWhkJ0N01Gkcy/jHxjL+sbEc3nGU32avYfns1bwy8T2CQgLpfUU3LpnQj26XdapwyE4pRU56Lmkn00lNTGfdz5tYMO1XAoMDuGnqNVz14EhCI4z79a6U4t2pP1CQV8Q/X7qagEDfrjd5Nr/sTQDgUhP2vylnLjhOIKFXGy3FlGiDKwtnDM5+zFgdFaBZzRjqR0eZzuAAbrx7CDs2HubdZ+fRsn196jU2dnqTJh0acWuHRtz8wnh2r93Pb1+tZuXcP/h9zloiosPpd2UvLrmuLx36t8FqdZlAQW4BqYnppJ3McD9cJpaW5Hqe5n6v2PZXpqvFauHy24ZwwzNXEVPH+DKtv87fzNplu5n8yAgaNTdfqahf9h6gZWxNmtQw/lqWG7u7YFSAbsFVBG1wZUAskSiJRjmOn3P0ui8jIgxs0ZRvt+yksNhOSKB5/uTWACuPvTaef4x9h/88OJvXvridkDDj6x1aLBba92lN+z6tueutSWz+dQe/zV7F8q9Xs+jTZcTUrU54VChpJzPOGcYMiwylRlx1asRVp12f1tSIi3G/di3rNa9D9drRlf/BzkFyYgYfvvgT7bs2ZsxN3ptyx1uczs1j07FE7u5vjlJi/8MZgwvUQwQqgnnudkYT0NQ12NKEDG3dnC82bmXxnv2M6WiuEGututH88z9X8/zdn/PMHTN47oOJhPrQwOKAwAB6DO9Mj+GdKcwvYt1P8az63jWpardhnagRF0PNejF/GlqNuBjCIs2RLJCdkcfTUz5DRHjwP1djtZqvy/67rbtQwPC25qwComxbQar/VVFJUy60wZWVwE6Q/wVK2XCNWzcPPRvVp1nNGD7fsJXRHdoYkqRwIfQc1IaHX7mWVx+dw9O3z+D5jyb5ZPWMkLBgBl7bh4HXmq+lczaF+TaeuWMmp05k8OInt1C3gcmKEwM2h4MvNm6lb9NGNKvlvdnMvYVSCmxrIbiX6f5nfQXz/SQzCAnqDNigeI/RUsqNiHBD907sTEpmW6J5Bn2XZODlF/HYa+PZu+0YT02eTl6OOWZJMCP2Ygf/fuArEnae4NFXrzVVtZKSLNq1n9O5eUzs2cVoKRXDcRicyUhQb6OVmBZtcGUlsJNrWbzFUBkVZXTHNkQGB/P5BnPqB9fM34+/cR37d57gycmfkpttrimMzIBSiref+Z6NK/fxj2dGV+pAe0+ilGLm+s00qxlDv2aNjJZTMWzrXEttcBVGG1wZEWttsMShTGpw4UFBXNmpHYv3JJCck2u0nArTZ2g7nnr7eg7uSeKJWz4lJ9M8xZHNwGdvLOHXHzZzw92DGXFtT6PlVJj4Y4nsOpXCxJ6dTRveU0VrwRIH1oZGSzEt2uDKQ1AnsG01WkWFub7bRTicTr7ZtN1oKRdEr0va8vQ7N3B4/ykeu+UTsjLyjJZkeuzFDj5+eQFzP/mdEdf25Lq7zFNJ71zMWL+Z6NAQRncwV1LVGZRygm297n+7QLTBlQMJ7AzOJJTDnP1YDWOiGdiiCd9s3mGaGQbOR4+BrZn6fzdx/OBpHpv0CZnp5m2VGs3pU1k8OvFjvp+xmpETenHX01eY+qZ6LD2TZfsOMr5rR1MNi/kv7HtBZer+twtEG1x5+LMfbquRKi6IG7p3JjUvn0W7E4yWcsF07duS5z6YSNKxNB6b+AkZqeYpR+YrbFq9n7vHvcvhfUk8+tp4/vHMaFMOByjJrA1bCLBYuL7bRUZLqTi2ta5lkEnH7/kI5v4mVzaBbYBglM2c/XAAfZo2pGmNGD7faN7PUJLOFzfnuQ8ncioxnUcnfkx6SrbRkkyBw+Fk1tu/8PSUGcTUjOCdb+9m4OUmNgQ32YWFfLd1F5e3b0Ws2eZ9K4EqWgvWpq6+f02F0QZXDkSCILCdqVtwriEDF7HjZDLbEs01S/n5uKhnM/710c2cPpXFIxM/JjW5civ+m4300zk8eeunzP5wOUPGduHNb+6ifpNaRsvyCHO37CS/uJhJZh0aAChlg+J4CNbhyQtFG1x5CewMxTtdX0KTMrpjWyKCg5i53j9acQAdujfhhY9vJv10Do/cNI0Th08bLcknWb98D3ePe5e9247zwItX8uCLVxESaq7CBeej2D2wu0ej+rSpE2u0nIpTvN01g7cOT14w2uDKiQT1AIrBtsFoKRUmIjiI8V06snDXPnYnpRgtx2O069KYFz+5hbzsQu69+n3mzVxNYYF5f4h4kpSTmTx/9+dMvWsWkdVCefObu7h0XDejZXmUrzdt52RWDpMvNvfnUoW/AoF6/JsH0AZXXoJ7g4S5v4Tm5fa+3YkOC+U/S393lQTyE9p0ash7399Dyw71mfbSAiYNeYU5H/9Ofl6R0dIMwV7sYO4nvzNl5BtsXpPAzQ9exnvf30OTlnWMluZRsgsLeW/lOno3aUD/Zo2NllNhlHJC4SII7otYooyWY3q0wZUTkRAI6gdFv7q+jCYlKiSEewb0ZsPREyzbd9BoOR6lVt1oXvpsMq9+PoXmbeP47I3FTBz8Ml++/ys5WVWn+snO+MPcPe5dpr++mM69m/PRzw9wzW0DCAwyaer83/DB6g1kFRTy6JD+ph7iQPFWcCYhISOMVuIX+N83vRKQkCGooiVQvAOCzJt5dm2XDny1cSuvLFtF/xZNCLKaayLL0mjfrQkvdGvCvh3H+eajFXzx3jK++2w1V1zfmzET+xAdY94su78jMz2X6a8tZum8TcTWjebZ92+k1yXmHPBcFnYmJTNz3WbGdWpn7r43QBUuAIIh2NwD7X0FMVN4qlu3bio+Pt5oGShnFiqlN4TfgiXyIaPlXBArDxzhttnzeGxof27u1dVoOV7l8L4kvv5oBasW7yAoOIAR1/bgylv6UyPWP0JB2Rl5LF+wjS/e+5WCvCLGTerHdXde4hNz6HkLm93OuE++IquwkJ9vv4lqoSFGS6owSjlQp/tDYGcs1d8zWk6lICKblFJe6zTVBldBnOmTwHEKS63FRku5YCZ/NY+tJ5L45e6biQkzx1xlF8LxQyl8M20Fy3/ehsUiDLuyG1dPHkDteuab8Tk3u4C1y3azctF2tqw9gMPupH23Jtz9zGgatfD/MVRv/raGD9dsYNr4MQxoYc5ZD86gitajMm5Eqr2FhFaNEKU2uBL4ksGpvC9QOc8jNRchAc2MlnNBJKSkMnraF4zv2pFnhl9itJxKI+l4Ot9++ju/fL8JpRSXXNGZa28bSL3GNY2W9rfk5xWx7rfdrFy0g02r92MvdlC7XnX6XdaBAcM70qxtnLn7ocrI9pOnuHb614zp2Jb/XHGp0XIuGGfWM1D4I1LrD8QSZrScSsFwgxOR6cBIIEUp1d69Lgb4BmgMHAGuUUplnLVfK/c2Z2gKPKOUektEpgK3AWcGKz2hlFpYmlifMjhHEur0ACTiISRiitFyLpipC5cxZ/MOfrz9RpqbcHLIC+H0qSy+m76SRXM2YC920H94R8ZO6kuzNnE+U7aqsMDGhhV7WbloOxtX7sNWZKdG7Sj6X9aB/sM70qpjgyphamcostsZ+/GX5Nls/HT7jUSFmDc0CaCUHZXSB4IvxhL9ptFyKg1fMLj+QC4wq4TBvQKkK6VeEpHHgOpKqUf/5hhWIBHoqZQ66ja4XKXUa+UR60sGB+BMHQcSgKXGHKOlXDDpeflc+v4MOjeoy8cTxhotxxAyUnP4fsZqFsxeR0G+jeDQQJq1iaNFu3quR/v61Gtcs9JMz1ZUTPyq/axctJ11y/dQVFBM9ZoR9B3mMrW2nRtisfiGAVc2r/66ik/WxvPJdWPpZ+JhAWdQRatRGbcg0e8jIUONllNpeNvgSs2iVEqtFJHGZ60eDQx0P58JrADOa3DAYOCgUupo+SX6LhIyFJX7JsqRbPqacTHhYdzZrwev/LqKVQeP+MVNo7xUrxnJrQ8N5+pb+7Nx5T4SdiWSsCuRxd9uZP7nfwAQGhbkMr32LsNr0a4ecY1qXLDRKKXIzysiN6uAownJrFy0nbW/7SY/t4io6uEMvqIz/Yd3pH23Jj7TqjSKrSeSmL5uE9d0bu8331NVuAAkAoL7Gy3FryhTH5zb4H4u0YLLVEpFl3g/Qyl13h56d5hzs1LqPffrqcAkIBuIB/55dojzXPhaC04VJ6DSLkeinkXCrjdazgVjs9u5/MNZBFmtzLvteoIC9CgScBUmPnHoNPt3nSBhZyIHdiVycM9JbEWuKYdCw4Np2roujZrH0rB5beo2iCEkNIj83ELycgrJzSkkP6eQPPfrsx85WfnkZBXgdPw1rjIiKoSLh7Sn/4gOdOrZDGuAfw3hqCi5RTau+vQriux2frr9RiKCg42WdMEoZUOlXAzBl2CJfsVoOZWK4SFKt4jGVNDgRCQIOAm0U0olu9fVBlIBBfwLqKuUuuU8+08BpgA0bNiw69GjvtMIVEqh0kYBQVhqfm+0HI/we8Jhpnz9A7f27sojQ/SvyfPhsDs4djCF/TsTSdh5gsP7T3HsQDK52YXn3Sc4JJCwiBAiokIIiwghPDKEiMgQIqqFERkdSmS1MCKrhVKrTjTtuzX2ywHZF4JTKe6e8xMrEg7x2Q1X0rNxA6MleQRV8AMq6xGk+gwk+GKj5VQqhocoz0OyiNRVSiWJSF3g7woaDsfVeks+s6LkcxH5GPj5fDsrpaYB08DVgqugXq8gIhA6HpXzPKp4BxLYwWhJF8yAFk24tksHpq/dxIDmTfzmJuJprAFWmrSqS5NWdRl2pev/UylF+ukcUhIzKCqyEx4ZQnhEMOFRoYSFB2vDukDe+30ty/Yf5KlhA/3qe6nyvwBrU1170gtUNJj/IzDR/XwiMP9vtp0AzC65wm2KZxgL7KygDuMJHQ0SisqfXfq2JuGxoQNoFBPNo/OXkF14/haJ5r8REWrERtGmcyM69Wrm7p+rSbXq4drcLpAlexJ4f9V6ruzUjhu6dzJajsdQtm1QvB0Ju75KZcFWFqUanIjMBtYCrUTkhIjcCrwEDBWRBGCo+zUiEiciC0vsG+Z+/+z43SsiskNEtgODgAc88mkMQCyREDISChegnP4x2WZYUCCvjh1OSk4uzyxY5lfFmDXmY2/yaR6dv5jO9esydfglfmUEKu9DkGoQWjUzl71NWbIoJ5znrf8plqaUOgmMKPE6H/ifQVVKqRvLodHnkbAJqIK5UDAfwv3jo3WMq8P9gy7m9d/W0LxmDHcP0OETTeWTnl/AXXN+JCokhHeuGulXiU+qeD8ULUMi7kEs/lkX1Wiqdr6xh5DA9hDQHlXwjV+1dm67uDtjO7bl3ZXr+GH7bqPlaKoYxQ4H9337M6dz8njvmlHERvqXCai8aSBhEOYfP4p9EW1wHkLCJoB9PxRvNlqKxxARnh85hF6NG/DUT0tZf+S40ZI0VYiXlq5kw9ETvDByKB3j/Gv+OmU/BoU/Q+gExBJttBy/RRucpwi5HCTCr5JNAIKsVt69eiSNYqK5e+5PHDydZrQkTRVg7pYdfLFxK7f07srojm2MluNxVN7HQAASfrPRUvwabXAeQixhEDoGChejnOlGy/EoUSEhfDRhDEFWK1O+/oG0vHyjJWn8mE3HE3lu4W/0bdqIhy7pa7Qcj6Mcp6Dgewi9CrGae/46X0cbnAeR0PGADQrmGS3F49SPrsYH40eTmpvPnd/Mp7DYbrQkjR+SlJXDvXN/Ji46ijfGjcDqh7U2Vd5ngBMJv9VoKX6P/317DEQCW0JgV1T+bJRyGC3H43SMq8Pr40awPfEU93/3Mza7NjmN58gpLOKuOT9SUGzng2uuMPXkpedDOVKh4GsIGYkE+M9gdV9FG5yHkfBJ4DgGhQuMluIVhrRqxtQRg1mecJg75/xIQXGx0ZI0fkBukY3bZs8jISWVt668nGZ+OmWTyn0LlB2J+IfRUqoE2uA8TfBQCGiFyn3fL1txAOO7duTFUUNZc/AoU2b/QJ7NZrQkjYk5mZXN9TPnsD3xFK+PG0H/5o2NluQVVPFuKJgLYTcgAY2NllMl0AbnYUQsSMTd4Djst604gKs6tefVMcPZdCyRW7/8npzCIqMlaUzItsQkrv50Nicys/howhiGtWlhtCSvoJRC5fwHpJpuvVUi2uC8QRVoxQGM6tCaN6+8nJ0nk5n0xXdkFui6lZqys3DXPm6cNZeQwAC+uXm838ztdk6KfgXbeiTyfsQSZbSaKoM2OC9QVVpxAMPatODdq0exLyWVm2bN1UMINKWilOK9let44PuFtK9bm7m3TKC5n/a5gXu+t5yXIaAFhF5jtJwqhTY4b1FFWnEAg1o25aPxozmansmNs+aSnJNrtCSNj1JYbOeheYt49/e1jOnYhhk3XElMeJjRsrxL/ufgOIZEPoGI/9TSNAPa4LxEVWrFAfRp2oiPrxtLUnYON86ay8ks/5hZQeM5TufmcdPnc/l51z7+eUkfXrpimF8VTz4XypGGyn0fggciwX2MllPl0AbnTapQKw6gR6P6TL9+HGl5+dwwcy7HMzKNlqTxEfYmn+aa6bPZl5zKu1eNZEqfHn417c35ULlvgypEIh81WkqVRBucF/mvVpwfVjc5F53rxzHzxqvIs9m4fuZc9iWnGi1JYzDL9x9iwoxvsDudfDnpGi7100zJs1HFe6FgDoRdjwQ0M1pOlUQbnLcJvhQCO6Ny30A5q0bfVPu6tZl541U4lZNrps9m7padfjWNkKZsKKX4bN0m7vxmPk1qVOfbW66jfd3aRsuqFP4aFhClhwUYiDY4LyMiSNST4ExF5b5jtJxKo3XtWvxw2w10aRDHUz8v5Z65P5Oik0+qDLlFNh7+YTEvLV3Jpa2b8+XEa6gd5V/zuf0tRYvBttY9mWm00WqqLNrgKgEJ7Aih10H+TJRto9FyKo2aEeF8ct1YHh7cj5UHDzPig1nM3bJDt+b8nJ1JyYz75EsW7NrHvQN689ZVIwkNDDRaVqWhnOmo7OcgoD2ETTBaTpVGzHSz6datm4qPjzdaRoVQzjxU2hWAQmr8hFjCjZZUqRxJy+DpBb+y4egJejVuwL8uH0LDmGijZWk8iFKKmeu38NqyVdSICOP1scPp1rC+0bIqHWfm/VC4FKkxz1WAXXNeRGSTUqqbt46vW3CVhFjCkWovgSMRlfOq0XIqncY1qjPzxqt4/vLB7ExKZtRHnzN97SbsTqfR0jQe4ODpNCZ/NY//LP2d/s2b8MNtN1RJc1OFS6BwIRJxtzY3H0C34CoZZ/Z/IP8zpPpnVXZcTHJ2LlMXLeO3/YdoX7c2L44aSuvatYyWpakAaXn5vPv7WuZs3kFoUCAPDurDdd0uqhJDAM5GOdNRqZeDpQ5SYw4iVScsW1G83YLTBlfJKFWISh0DKh+puQCxRBotyRCUUizavZ9/LV5OdmERt13cnbv69fD7gb/+QmGxnVkbNvPh6o0UFhczvmtH7u7fy/+rkvwNzswHoXAJUuM7JLC10XJMgbcNTt9NKhmREKj2Mip9PCrnRVfYsgoiIoxo14reTRry0tLf+WD1en7Zm8ALI4fSpUGc0fI058GpFAt27uON5as5mZXDoBZNeXhIP5rVjDFamqGowl+g8Gck4j5tbj6EbsEZhDPnTcj7AIn+AAkZbLQcw1l54AjPLvyVpKwcbujeiQcu6UN4UJDRsjQliD92gpeWrmTHyWTa1onl0aH96dVYz0qtnBnu0GQtpMa3OjRZDgwPUYrIdGAkkKKUau9eFwN8AzQGjgDXKKUyzrHvESAHcAD2Mx+krPufjT8ZnFI2VNqV4ExzhyqrGy3JcHKLbLy1fA1fbNxK3WqRPH3ZIAa1aFol+3N8iSNpGby2bDVL9x2gdmQEDwy6mNEd22LRfxfXgO6sB6DwF3doso3RkkyFLxhcfyAXmFXC4F4B0pVSL4nIY0B1pdT/FFtzG1w3pVTqWevLtP/Z+JPBAajiPai0qyBkKJbot4yW4zNsPn6SJ39ayqG0dNrWiWVy724Ma9uCAItO+q1MMgsKeX/lOmbHbyPQauW2Pt25uVeXKjWmrTRU/leo7KlIxANIxJ1GyzEdhhucW0Rj4OcSBrcPGKiUShKRusAKpVSrc+x3hHMbXJn2Pxt/MzgAlfsBKvdNpNpbSOgIo+X4DDa7nfk79vLp2ngOp2XQoHo1bu7VlXEXtdU3WC9js9v5Mn4b/7dqPblFNq7q1I57B15MrYiqNXazNFTxdlTaBAi+GIn+CBH9A6y8+KrBZSqloku8n6GU+p8Ym4gcBjIABXyklJpWnv3Pxi8NTtlR6ePBftQV4ghoaLQkn8KpFMv2HeTjPzayLfEUUSHBjOnYlvFdO1b5xAZPk2ez8cO23Uz7YyOnsnPp26wRjwzuT6vaNY2W5nMoxylU2jWAFak5T5fjqiBmN7g4pdRJEYkFlgL3KKVWlsfgRGQKMAWgYcOGXY8ePVqez2cKlP0YKm0cWOshNb5GJNRoST6HUopNxxP5Kn47v+xJoNjppGfjBkzo2pEhrZoRaLUaLdGUFNntrDxwhAW79rF8/yEK7Xa6NIjjnv69uLhpI6Pl+STKmYtKvw4cx5GYL5HAtkZLMi2+anDlDjGKyFQgVyn1mg5R/i+q6HdUxhQIuQKp9opOrPgbUnPz+G7rLr7ZvIPErGxqRYRxVacOXNulA3WrVc1xheWh2OFg7eHjLNy1j6X7DpBbZCMmLJTL2rZkZPtWdKkfp79/50EpOyrjdrD9gVSfhgT3M1qSqfHVcXA/AhOBl9zL+WdvICLhgEUpleN+finwfFn3r2pI8ACIuNc1QWJgRwi/0WhJPkvNiHBu79uDyRd3Y9XBI3wVv50PV6/nozUbGNiiCRO6XkTfZo10ll8JnEoRfyyRBbv2sWRPAhn5BUQGB3Np6xZc3q4VvZo00Ek8paCUchVRtq1Col7Q5mYCypJFORsYCNQEkoFngR+AOUBD4BhwtVIqXUTigE+UUiNEpClwZpbPAOArpdSL7mPWONf+pYn15xYcgFJOVOZdULQSiZmFBHnth43fcTwjizmbd/Dt1p2k5xdQKyKcUR1a06tRA7o2rEdEcNUbU6eUYsfJZH7etZdFu/eTkpNHaGAAg1s2Y0S7VvRr1khXjikHKvcjVO7rEH4HlsgHjZbjF/hEiNJX8HeDA1DObNf4OJXvqkZujTVakqmw2e0s3XuQn3buZfWhoxQ7HFhFaB9Xm56NGtCzcQO6NIgjLMg/MzGVUuxLSWXhrn0s2LWPE5nZBFqtDGjemMvbtWJgi6Z++9m9iSr4GZX1IISMQqq9pkO4HkIbXAmqgsEBqOL9qPSrIaCNqyUnVa/14QkKi+1sOXGSdUeOs/7IcXacTMbudBJosdCxXh16NnYZXuf6dQk2aUsmPb+APadS2JWUwp5TKexISuZ4RhZWEXo3bcjIdq0Y0qo5kSHBRks1Lcq2EZU+CQI7ITGf6f9HD6INrgRVxeAAVMECV4WEsBuwRD1jtBy/IM9mY/Pxk6w7fJz1R4+zKykFp1IEWa10ql+Xno0b0KtxAzrWq0OQj2VlKqVIyc1jd1IKu06lsDsphT3JKZzMyvlzm/rRUbStE0vvJg25rE2LKl342FMo+yFU2rVgiUFqfKOHA3gYbXAlqEoGB+DMfgnypyPVXkZCxxotx+/IKSwi/lgi648eZ93h4+xNPo0CQgMDaFe3NnWiIqkdGU5sZASxEeHUjoogNiKC2pHhXu27UkpxIjOb3W4j23XK1TpLzcsHQIAmNWJoW6cWbevG0rZOLG3qxBIdGuI1TVUR5UhzRVJUARIzBwnQdTc9ja9mUWoqAYl8CGXfjcp6CiyxVXb+OG8RGRLMoJZNGdSyKeAqTbXx6AnWHznOrlMpbD1xkpScPGwOx//sGx0aQmxkBLUjI4iNDP9zGRvhWgLk24rJtxVTUFxMnq2YfJuNgmI7+TYb+TbXuoLiYvd2rnX5xcWk5uaTU1QEQIDFQvNaNejfvDFt68TStm4srWvX0oWovYxSBajM28GRitT4QpubSdEtOB9HObNR6de7BpVWn4kEXWS0pCqFUorMgkJScvJIyc0lOSfX9TynxPPcXFJz83GW8X8pwGIhLCjQ9QgMIjQogPCgIPfrQKLDQmlduyZt68TSMramafsHzYpSDlTmvVD0KxL9PhIyxGhJfotuwVVxxBIF1T9FpU9AZdwC1T9BgjobLavKICJUDwuleljo35assjudpOXlk+I2PYsIoW4TCw8MdD93mZiv9e9p/kIpJyr7SShaikQ+pc3N5GiDMwFijYWYz1HpE1EZN0P0h0hwL6NlaUoQYLFQ2x2y1JgTl7k9AwXfIxH3IOE3GS1Jc4Ho0gUmQaxxSMyXYI1DZdyGKlpptCSNxm9QSqFynoeCORB+J4TfbbQkjQfQBmcixBqLxHwBAc1QGXeiCn8xWpJGY3pc5vZvyP8KwicjEffrgdx+gjY4kyGWGCRmFgS2Q2Xehyr4yWhJGo1pUaoYlf0U5M+EsElIxMPa3PwIbXAmRCxRSPXpENQNlfUQKn+u0ZI0GtOhnLmojDugYC6E34VEPq7Nzc/QBmdSxBKBVJ8GQX1R2U+i8mYZLUmjMQ3Kcco1p5vtDyTqRSyROizpj2iDMzEioUj1DyB4CCrnBVTuR0ZL0mh8HlW8zzUbt+O4a063sKuNlqTxEtrgTI5IEBL9NoSMROW+jjPnbcw0eF+jqUxU0RpU+gTA6ZqNW8/p5tfocXB+gEggVHsVJSGQ9z5K5UPkYzrkotGUQOV/70ooCWjmarlZ6xotSeNltMH5CSJWiHoBJaGQ/xlKFULUM671Gk0VRiknKvddyHsfgi5Got9FLJFGy9JUAtrg/AgRC0Q+5W7JfYxypkC11xGLnjZFUzVRzhxU1iNQtAxCxiLV/qXnc6tC6D44P0NEsEQ+jEQ+DUXLUenXoxwnjZal0VQ6yn4AlXYVFK1AIp9Aqr2kza2KoQ3OT5HwG5HoD8BxBJV6BapwsdGSNJpKQxUucpmbykZiZiLhk3SfdBVEG5wfIyGDkBrzIaAxKvNenFlPoJx5RsvSaLyGUnacOa+gMu+DgJZIjXlIUA+jZWkMQhucnyMBDZGY2RB+BxR8h0obhyreabQsjcbjKGe6a0qpvE8g7Dok5gvEWsdoWRoD0QZXBRAJxBL5IFJ9Fqh8VNq1qLxPUcpptDSNxiMo2zZU6liwbUaiXsISNVX3t2m0wVUlJLgnUvMnCB6EynkZlXELypFstCyNpsIoVYQz5w1U+njAgtT4BgkbZ7QsjY+gDa6KIZZo1zigqBegeAsqdRSqcJnRsjSacqNsW1GpYyDvQwgdg9T8AQlsZ7QsjQ9RqsGJyHQRSRGRnSXWxYjIUhFJcC+rn2O/BiKyXET2iMguEbmvxHtTRSRRRLa6HyM895E0pSEiSNg1SI3vXROoZt6JM3uqa3C4RuPjKFWAM/s/qPRrQeUj1T/FUu0/iKWa0dI0PkZZWnAzgMvOWvcYsEwp1QJY5n59Nnbgn0qpNkAv4B8i0rbE+28qpTq5HwvLL11zoUhAM6TGHAi7BfK/QqVdiSrea7Qsjea8qKL1qNRRkP8ZhE5Aai7Q9SQ156VUg1NKrQTSz1o9Gpjpfj4TGHOO/ZKUUpvdz3OAPUC9CxGr8TwiQViiHnPNL+fMRKVdhcqbgVIOo6VpNH+inNk4s6aiMm4EFFL9cyzVpiKWCKOlaXyYivbB1VZKJYHLyIDYv9tYRBoDnYH1JVbfLSLb3SHQ/wlxlth3iojEi0j86dOnKyhXUxoS3NedgNIXlfNvVwUU+wGjZWmqOEo5UfnfoVIvhYKvXbNu1/gJCe5ptDSNCfB6komIRADfAfcrpbLdqz8AmgGdgCTg9fPtr5SappTqppTqVqtWLW/LrdKIJQaJ/gCp9irYD6FSR6Ny30M5842WpqmCqOJdqPQJqOzHwdoIqfE9lqgndG1VTZmpqMEli0hdAPcy5VwbiUggLnP7Uin1/Zn1SqlkpZRDuQZifQzoUgM+goggoaORmgshZAgq9x1U6hBU3hcoZTNanqYKoJxZrqSntCvBcRSJegmJmY0Eti11X42mJBU1uB+Bie7nE4H5Z28grsJvnwJ7lFJvnPVeyYmYxgK6tIaPIdaaWKLfRmK+BmsTVM7zqNTLUAXzdf+cxiu4wpHfok5fCvlfQ9j1SM1fkLBxrpkyNJpyUpZhArOBtUArETkhIrcCLwFDRSQBGOp+jYjEiciZjMg+wI3AJecYDvCKiOwQke3AIOABz34sjaeQoC6ukkfVPwGJQmU9jEobjSr8Tc8crvEISilXdmT6eFT2ExDQxB2OfBqxRBktT2NixEw3qW7duqn4+HijZVRZlHJC4SJU7tvgOAKBnZHIf+pitpoKoZQC20pU7gdQvBkstZDIhyBkjK78X0UQkU1KqW7eOr6e8FRTZkQsEHo5hFwKBd+7ElDSb0AF9XMZne4j0ZQBpZxQtBSV+yHYd4ElDol8BsKuQiTEaHkaP0IbnKbciARC2LUQOhryv0TlfoRKG4MKGYFE3I8ENDZaosYHUcoOhQtReR+C/YArMzLq3xB6hS6MrPEK2uA0FUYkBMJvhdBrUHnTIf8zVOESVOiVSMTdeqoSDYAr+7bgB1TeNHAcc83TVu0NCBmOiNVoeRo/Rhuc5oIRSyQSeR8q7HrXr/P82a5sy7AbkPCbEevf1gHQ+ClKFUL+HFTeJ+A8BQHtkej/g+BLdFakplLQBqfxGGKtiUQ9hQqbhMp919Wiy5+JCh6MhF0JQX0R0V85f0c5MyD/a1T+LHCmQWA3pNq/IaiPTh7RVCr6bqPxOBJQH4l+GWW/C5X/tSshpWgJWGJRoaOR0HFIQDOjZWo8jLIfQOV9AQXfA4UQ1A+JuAMJ6m60NE0VRQ8T0HgdpWxQtAJVMA+KVgAOCOyEhI6DkBF6rJOJUfZjULgAVbgQ7PuAIFfSSNgkJLCl0fI0Po63hwlog9NUKsqRCoXzUQXfgz0BCIaQoUjolRDUSycdmADlSHSNhyxYCHZ3EaLAzkjICAgZiVhrGCtQYxr0ODiNXyHWmq7My7BbwL7TZXQFP6MKfwZLXVToGHcIs5HRUjUlUI5kKFyMKlwAxVtdKwPaI5GPurIhrXGG6tNozoVuwWkMR6kiKFzmMjvbasAJgV2R0NEQ3A+x6mkEjUA5UqFoiaulVhwPKAho7W6pjUACGhotUWNydAtO4/eIBEPoCCR0BMpxCgp+RBV8h8p+BgBlbQzBfZCgPq4wpp7k0msoZwYU/uLqU7OtB5xgbYZE3ONqqenkII2J0Aan8SnEWgcipkD4beA4CEVrULY1UDAPlf8lYEUFXoQE94GgPhDYUQ89uACUMwNsm1G2eFcrrXgXYAdrIwi/3dVaC2ip0/s1pkTfGTQ+iYhAQHMIaI6ET3RlYhZvRRWtBtsaVO57wLsgEaigXm7D6wvWhvpm/DcoRyLY4t2GtslVMguAQAjsAOGTkZBhENBWX0eN6dEGpzEFIkEQ1MM9c8GD7pbHOlTRGlcrr+hX14bW+qigi12GF9ABrHFVtmqGUk6XgRW7Dc22CZxJrjclAgK7IBGjIKibuyUcbKxgjcbDaIPTmBKxVHf1CYUMd0274jjqatkVrXEV9C2Y494wDBXQDAJaIAEtwP0cS5xftVCUcoIzBexHoXg7qjgebJtBZbk2sMRCUDckcLLL0AJa6iEZGr9HG5zG9LjCmY0hoDESdr2ran3xTrDvRdkPuMbbFa1yZWn+uVM4KqApWJsiAU0hoKnL/KwNfbKyvVIKVAY4ToAjEewnXOFGx4m/Htj+2sHaBEIuRYK6QWBXsDbwK0PXaMqCNjiN3yESAEGdIKgTJW/pypnpCtnZD6DsCa7ntvWowvkltrKirA3A2gAsUSCRYIlEJAoskSBRYIlwLd3vuZ6HlslAlHICNlBFoArdyyKgqMS6XHCcMbDEv0xN5Z/1QauBtb6rrzJ4kCtt39oAAtsilpgLvo4ajdnRBqepMogl2hWeC+p2lvHluWYotx9C2Q+B45DLVIpPgMoGZzaK4lKObkX9aXiRgHIbl+0sI7OVcpySgiNcBmZt6CpUbK0H1nrudfUQS2Q5r4BGU7XQBqep8oglHCztILAd52uDKVUEzmxQOeDMcRvfX0ulSq7LASwgwUCwa3nmQbArmUNC3OuCgJD/3kbC3QZWrdKugUbjj2iD02jKgEgwWGsBtc79fuXK0Wg0ZaBq5k9rNBqNxu/RBqfRaDQav0QbnEaj0Wj8klINTkSmi0iKiOwssS5GRJaKSIJ7Wf08+14mIvtE5ICIPFbe/TUajUajqShlacHNAC47a91jwDKlVAtgmfv1fyGuMgnvA8OBtsAEEWlb1v01Go1Go7kQSjU4pdRKIP2s1aOBme7nM4Ex59i1B3BAKXVIKWUDvnbvV9b9NRqNRqOpMBXtg6utlEoCcC9jz7FNPeB4idcn3OvKur9Go9FoNBXGm+PgzjU0qNzTh4vIFGCK+2VRyb5Ak1ETSDVaRAUwq24wr3az6gbzajerbjC39lbePHhFDS5ZROoqpZJEpC6Qco5tTgANSryuD5wsx/4AKKWmAdMARCTem9ObexOzajerbjCvdrPqBvNqN6tuML92bx6/oiHKH4GJ7ucTgfnn2GYj0EJEmoirPPt4935l3V+j0Wg0mgpTlmECs4G1QCsROSEitwIvAUNFJAEY6n6NiMSJyEIApZQduBtYAuwB5iildrkPe879NRqNRqPxFKWGKJVSE87z1uBzbHsSGFHi9UJg4Tm2SzvX/mVgWgX28RXMqt2susG82s2qG8yr3ay6QWs/L6JUufM+NBqNRqPxeXSpLo1Go9H4JZVqcOcp+zVVRBJFZKv7McK9vkeJddtEZGyJfSaIyA4R2S4ii0Wk5nnO97i7TNg+ERlmBt0i0lhECkoc48OK6vaw9mvduneJyCt/cz5fu+al6jbympd4v6GI5IrIQyXWdXV/Xw6IyDsi554y3IhrfqG6ffiavygix0Ukt5TzeeSaV7Z2T153T+gWkTARWSAie93/o+fNx6jQNVdKVdoD6A90AXaWWDcVeOgc24YBAe7nZ4YSBLgfKUBN93uvAFPPsX9bYBsQDDQBDgJWE+huXPI8PnLNawDHgFru92YCg01wzcuq27BrXuL974C5JbcBNgC9cY0pXQQM95Vr7gHdvnrNe7m/Q7l/s5/HrrkB2j123T2h2/2/O8j9PAhY5cnveaW24NS5y36db9t85crEBAjhr0Hi4n6Eu38ZRvHX+LqSjAa+VkoVKaUOAwdwlQ/zdd0exUPamwL7lVKn3a9/Ba48xyF87ZqXVbdHKY92ABEZAxwCdpVYVxeIUkqtVa7/8Fmcu6SdIdfcA7o9iie0u4+zTrmrLP0NHrvm7nNWpnaP4Qnd7v/d5e7nNmAzrjHTZ1Oha+4rfXB3u8NI06XEzAIi0lNEdgE7gDuUUnalVDFwp3vdSVzO/uk5jvl3pcJ8WTdAExHZIiK/i0g/D2sut3ZcX6bW7vBGAK4bVoNzHNOnrnk5dINB11xEwoFHgefO2rYerut3hvNdS0OuuQd0g+9d87JSGdccvKMdvH/dK6RbRKKBUbgK8J9Nha65LxjcB0AzoBOQBLx+5g2l1HqlVDugO/C4iISISCAuo+gMxAHbgcfPcVyPlAozQHcS0FAp1Rl4EPhKRKI8qLvc2pVSGW7t3+AKIRwB7PwvPnXNy6HbyGv+HPCmUursfpOyXkujrvmF6vbFa15WvH3NwXvavX3dK6Tb/QN0NvCOUurQuTY5x7pSr7k3a1GWCaVU8pnnIvIx8PM5ttkjInlAe9wfVCl10L3PHM493c7flQrzWd1KqSKgyP18k4gcBFoCHitpUwHt8Uqpn4Cf3PtMARznOLSvXfMy6Tb4mvcErhJXAkw04BSRQlx9FSVDNee7lkZd8wvS7YvXXCn1XhkP7dVrDt7T7u3rfgG6pwEJSqm3znPoil3z0jrpPP3grE5OoG6J5w/girOCqyPxTOJAI/eHqYmr9ZPEX4kD/wJeP8d52vHfnZKHuLCO4MrSXeuMTlx9SIlAjJHX3P061r2sDmwFWvr6NS+HbsOu+Vn7TOW/kwY24kocOJOsMcJXrrkHdPvkNS+x/u8SNTx6zStZu0evu4e+Ly/g+mFk8fQ1r/AfpIIXYzaum3wxLke+FfgcV7/Jdlw1Kuu6t70RV2fkVlwdj2NKHOcOXOW/tuP6dV7Dvf4K4PkS2z2JK9tmH+fIzPFF3biSIHa5/5ibgVE+cs1nA7vdj/El1vv6NS9Vt5HX/Kz9pvLf//jdgJ3u6/kefxVmMPyaX6huH77mr7j3d7qXU715zStbuyevuyd042qJKVz3xa3ux2RPXXNdyUSj0Wg0fokvJJloNBqNRuNxtMFpNBqNxi/RBqfRaDQav0QbnEaj0Wj8Em1wGo1Go/FLtMFpNBqNxi/RBqfRaDQav0QbnEaj0Wj8kv8HZ3FH1iwqM8oAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# The objective function to be minimized.\n", + "\n", + "# chi2 definition\n", + "def chi2(args):\n", + " mu = args[0]\n", + " sigma = args[1]\n", + " \n", + " model = np.array([norm.pdf(e, mu, sigma) for e in x])\n", + "\n", + " chisq = np.sum(((y - model)/y_error)**2)\n", + " return chisq\n", + "\n", + "scanMu = np.linspace(1538.,1542.0,25)\n", + "scanSigma = np.linspace(10,12,25)\n", + "Z_binned = [[chi2((a,b)) for a in scanMu] for b in scanSigma]\n", + "fig1, ax2 = plt.subplots(constrained_layout=True)\n", + "\n", + "p2 = ax2.contour(scanMu, scanSigma, Z_binned, [10, 20, 50, 100])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "mu = 1539.8898 +- 0.1124\n", + "sigma = 11.2214 +- 0.0808\n", + " fun: 3.0017805277361016\n", + " hess_inv: array([[0.01264244, 0.00026046],\n", + " [0.00026046, 0.00652938]])\n", + " jac: array([ 3.87430191e-07, -9.83476639e-07])\n", + " message: 'Optimization terminated successfully.'\n", + " nfev: 39\n", + " nit: 11\n", + " njev: 13\n", + " status: 0\n", + " success: True\n", + " x: array([1539.88978469, 11.22137155])\n" + ] + } + ], + "source": [ + "# Perform the fit minimizing the least squares\n", + "\n", + "res = minimize(chi2, [1545, 9])\n", + "err = np.sqrt(np.diag(res.hess_inv))\n", + "print('mu = {:1.4f} +- {:1.4f}'.format(res.x[0], err[0]))\n", + "print('sigma = {:1.4f} +- {:1.4f}'.format(res.x[1], err[1]))\n", + "print (res)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot the result" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "def model_function(x, args):\n", + " mu, sigma = args[0:2]\n", + " return norm.pdf(x, mu, sigma)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAWAAAAEYCAYAAABiECzgAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABK+UlEQVR4nO3deVzVVf748ddhxwUVRMUFUBAlRc3QNBuzcktNm3QmyqXFxszsWzY2ozPNTPu0jVrZYpaJaWZj4/JzSzPLLM19AxUBWRURlX2/nN8f915CRL3AXVjez8fjPuR+PudzPm8+wptzz+d8zlFaa4QQQtifk6MDEEKIxkoSsBBCOIgkYCGEcBBJwEII4SCSgIUQwkFcHB2ANbVu3VoHBgY6OgwhhLjCgQMHMrTWvpW3N6gEHBgYyP79+x0dhhBCXEEplVjVdumCEEIIB5EELIQQDiIJWAghHKRB9QEL0ZiUlJSQkpJCYWGho0MRJh4eHnTs2BFXV1eLyksCFqKeSklJoXnz5gQGBqKUcnQ4jZ7WmosXL5KSkkLnzp0tOka6IISopwoLC/Hx8ZHkW0copfDx8anWJxJJwELUY5J865bq/n9IAhZCCAeRPmAhrCgzM5M9e/aQnp6Oj48PrVu3JjQ0FC8vL0eHJuogaQELYQX5+fmsWbOG9957j3379lFYWMixY8fYsmULH3/8MbGxsY4O0SZuu+22G5Z5/PHHiY6OBuD111+v9vHNmjWrWXBVSEhIoGfPnlarr7ZUQ1oRIzw8XMujyMLeioqK+OKLL0hLS6Nfv34MGDCAFi1aoLUmIyODb775hvPnz3PHHXdwxx13WK3f9sSJE4SGhlqlLntp1qwZubm5Nj/mWhISEhgzZgzHjx+3Sn1Vqer/RSl1QGsdXrmsdEEIUQulpaWsWrWKs2fP8sADD9CtW7fyfUopfH19mTp1Khs3buTHH3/ExcWF22+/3epxbNmyhbS0NKvW2a5dO0aOHHndMubk+MMPP/Diiy/SunVrjh8/zi233MLy5ctRSjFkyBDeeecdVq9eTUFBAX369KFHjx6sWLGi/Pjc3FzGjRvH5cuXKSkp4dVXX2XcuHE3jPGBBx7g4YcfZtSoUQA88sgj3Hvvvdxyyy1MnjyZvLw8ABYuXHhVa3vp0qXs37+fhQsXAjBmzBhmz57NkCFD2Lp1K//6178oKioiKCiIzz//nGbNmjFnzhzWr1+Pi4sLw4cP55133qnJpS0nCViIGtJas2bNGs6cOcN99913RfKtyNXVlXHjxlFSUsKOHTsIDAykY8eOdo7W9g4dOkRUVBTt27dn0KBB/Pzzz1f8sXnjjTdYuHAhhw8fvupYDw8P1qxZg5eXFxkZGQwYMICxY8fe8NNCREQEq1atYtSoURQXF7N9+3Y++ugjtNZs27YNDw8PTp8+zYMPPmjxRF0ZGRm8+uqrfPfddzRt2pQ333yTefPmMXPmTNasWcPJkydRSpGZmVmdy1MlScBC1FBUVBTR0dHcfffd9O7d+7pllVLce++9pKam8s033/DEE0/g4eFhtVhu1FK1h/79+5f/YenTpw8JCQkWt/a11vztb39j586dODk5kZqayvnz52nXrt11j7vnnnv4v//7P4qKitiyZQuDBw/G09OTrKwsZs6cyeHDh3F2diYmJsbi72PPnj1ER0czaNAgAIqLixk4cCBeXl54eHjw+OOPM3r0aMaMGWNxndciN+GEqIGioiK2bt2Kn5+fRTeSwNjKGz9+PFlZWWzatMnGEdqfu7t7+dfOzs6UlpZafOyKFSu4cOECBw4c4PDhw7Rt29aiBxo8PDwYMmQI3377LatWrSIiIgKA+fPn07ZtW44cOcL+/fspLi6+6lgXFxfKysrK35vPp7Vm2LBhHD58mMOHDxMdHc1nn32Gi4sLe/fuZfz48axdu9Yqf/QkAQtRAz/++CM5OTmMGjUKJyfLf406derE4MGDOXbsGCkpKTaMsG5ydXWlpKTkqu1ZWVm0adMGV1dXduzYQWJildPnVikiIoLPP/+cn376iREjRpTX5+fnh5OTE1988QUGg+Gq4wIDAzl8+DBlZWUkJyezd+9eAAYMGMDPP/9cPnIlPz+fmJgYcnNzycrKYtSoUSxYsKDKrpTqkgQsRDWlp6ezZ88e+vbtW6O+3Ntuu42mTZvy/fff2yC6um3atGn06tWLiRMnXrF94sSJ7N+/n/DwcFasWEH37t0trnP48OHs3LmToUOH4ubmBsCMGTOIjIxkwIABxMTE0LRp06uOGzRoEJ07dyYsLIzZs2fTt29fAHx9fVm6dCkPPvggvXr1YsCAAZw8eZKcnBzGjBlDr169uOOOO5g/f34troSRDEMTopq+/PJLkpOTefrpp2nSpEmN6tizZw/ffvstU6ZMsXjilsrq4zC0xqA6w9CkBSxENaSnp3P69GkGDBhQ4+QLEB4eTvPmzdmxYwcNqREkqkdGQQhRDbt378bFxYV+/frVqh4XFxcGDx7Mxo0biY2NpWvXrlaKsGE6duwYkydPvmKbu7s7v/76q4Misg5JwEJYKCcnh2PHjnHzzTfXqvVrdvPNN7Nr1y5++eUXScA3EBYWZpWbXnWNdEEIYaG9e/diMBgYOHCgVepzdnYmPDychIQEMjIyrFKnqF8kAQthgeLiYvbv309oaCje3t5Wq7dPnz44OTlx4MABq9Up6g9JwEJY4NixYxQWFlqt9WvWrFkzQkNDOXz4cJXjY0XDJglYCAscOXIEX19fm8zhEB4eTmFhIVFRUVav29bee+89QkNDadWqFW+88QYAa9euLZ9+UlyfJGAhbuDSpUskJyfTu3dvmywBFBAQgI+PT73shvjwww/ZtGkTly9fZs6cOYAk4OqwaQJWSo1USp1SSsUqpeZUsV8ppd4z7T+qlOpbab+zUuqQUmqDLeMU4nqOHDmCUoqwsDCb1K+UIjw8nJSUFKtPKWlL06dPJz4+nrFjxzJ//nxmzpzJL7/8wvr163n++efp06cPcXFxjg6zTrPZMDSllDPwATAMSAH2KaXWa60r/mm8B+hqet0KfGT61+wZ4AQg67kIh9Bac+TIEbp06WLTZYV69erFtm3bOHbs2A1nAKvKs88+a/VhWn369GHBggXX3P/xxx+zZcsWduzYwYYNxjbSbbfdxtixYxkzZgwTJkywajwNkS1bwP2BWK11vNa6GPgKqDzD8jhgmTbaA7RUSvkBKKU6AqOBT20YoxDXlZiYSFZW1g2nm6ytJk2a0KVLF6KiouTJuEbElg9idACSK7xP4crW7bXKdADOAQuAvwDNbReiENd35MgR3NzcqjU5TE316NGDdevWkZqaWu2bfddrqYq6y5Yt4KruVlT+015lGaXUGCBda33DuxJKqWlKqf1Kqf0XLlyoSZxCVKmkpITo6GhuuukmXF1dbX6+7t274+zsbNP1yuyhefPm5OTkODqMesGWCTgF6FThfUfgrIVlBgFjlVIJGLsu7lJKLa/qJFrrT7TW4VrrcF9fX2vFLgTx8fEUFxfbbRVdDw8PgoODiY6OrtfdEBEREbz99tvcfPPNchPuBmyZgPcBXZVSnZVSbkAEsL5SmfXAFNNoiAFAltb6nNZ6rta6o9Y60HTc91rrSTaMVYirnDx5Eg8PDwIDA+12zh49epCTk0NSUpLdzlkbCQkJtG7dmkceeaR8cctBgwYRHR3NoUOHCAoKcnCEdZvNErDWuhSYCXyLcSTD11rrKKXUdKXUdFOxTUA8EAssBmbYKh4hqsNgMHDq1ClCQkJwdna223m7deuGi4tLve+GEJax6WxoWutNGJNsxW0fV/haA0/doI4fgB9sEJ4Q15SYmEhBQYFdbr5V5ObmRkhICNHR0dxzzz3VWu5I1D/yvytEFU6ePImLiwvBwcF2P3doaCj5+fmkpqba/dzCviQBC1GJ1pqTJ08SHBxsl9EPlQUHB+Pk5MSpU6fsfm5hX5KAhagkNTWVnJwch6235uHhgb+/PzExMQ45v7AfWRFDiEpOnDiBk5OTQ1epCAkJYevWrVy+fJlWrVrVur7522J4d/vpq7Y/c3dXZg0LqXX9omakBSxEJTExMQQEBODp6emwGLp161YeizXMGhZCwhujubWzN7d29ibhjdEkvDG6XiXfIUOG0NBWPZcELEQFmZmZZGRkOHyNNm9vb1q3bi3dEA2cJGAhKoiNjQVwyOiHyrp27UpCQgJFRUVWqa/UUEbSpXyOp2bx1paTlBrKal1nXl4eo0ePpnfv3vTs2ZNVq1bx8ssv069fP3r27Mm0adPKn+obMmQIs2bNYvDgwYSGhrJv3z7uv/9+unbtygsvvAAYH+zo3r07Dz/8ML169WLChAnk5+dfdd6tW7cycOBA+vbtyx/+8Adyc3Nr/b04giRgISqIjY2lRYsWtG7d2tGh0K1bN8rKyqz2OO+8bTGczy4kr9jAkp/PMH9b7VvXW7ZsoX379hw5coTjx48zcuRIZs6cyb59+zh+/DgFBQXlU1WCcZzzzp07mT59OuPGjeODDz7g+PHjLF26lIsXLwJw6tQppk2bxtGjR/Hy8uLDDz+84pwZGRm8+uqrfPfddxw8eJDw8HDmzZtX6+/FESQBC2FSWlpKfHw8Xbt2tcnKF9XVqVMnPD09rdYN8UvcRcpMU0wUlpTxc9zFWtcZFhbGd999x1//+ld++uknWrRowY4dO7j11lsJCwvj+++/v2KppbFjx5Yf16NHD/z8/HB3d6dLly4kJxsnRuzUqRODBg0CYNKkSezateuKc+7Zs4fo6GgGDRpEnz59iIyMJDExsdbfiyPIKAghTJKSkigpKakT3Q8ATk5OdOnShbi4OLTWtf6jcFuQD0dSMtEaPFydGBTkU+sYQ0JCOHDgAJs2bWLu3LkMHz6cDz74gP3799OpUydefPFFCgsLy8u7u7sDxu/N/LX5fWlpKcBV32fl91prhg0bxsqVK2sdv6NJC1gIk9jYWJydnencubOjQykXFBREbm4u1phq9blhIfh5edDUzZmpgzpbZQTE2bNnadKkCZMmTWL27NkcPHgQgNatW5Obm8vq1aurXWdSUhK7d+8GYOXKldx+++1X7B8wYAA///xzeX99fn5+vb1ZKS1gIUxiY2MJCAjAzc3N0aGU69KlCwBxcXG0adOmVnW5ODvRybsJnYDnR1pnjotjx47x/PPP4+TkhKurKx999BFr164lLCyMwMBA+vXrV+06Q0NDiYyM5IknnqBr1648+eSTV+z39fVl6dKlPPjgg+U3KF999VVCQurPkDozVZ/nHa0sPDxcN7RxgsI+srKyWLBgAcOHD2fgwIGODucKH3zwAS1atGDSpCtnZD1x4kS1n9Z7YJGxZbnqibr1PZolJCQwZsyYej0bXFX/L0qpA1rr8MplpQUsBL8NPwsKCiIuLo7Lly/Tvn172rRpg4uLY39NunTpwsGDByktLa1xLJWfhAucsxGQJ+EcTRKwEMC6detYs2YNH374Ienp6eXbnZ2dGT58OP/3f//H8OHDHTI9ZFBQEHv37iUpKam8S6K6Zg0LqReJNjAwsF63fqtLErBo1I4dO8Zf/vIXtmzZgre3N2PGjGHQoEG0bduWc+fOERcXx/Lly7nnnnvo3r07S5YssXsXRWBgIE5OTsTFxV2VgK0xOkJYT3W7dGUUhGi0li5dSt++fdm9ezfDhg1j8+bNREZGMm3aNMaNG8f06dN5++23SUxMZMWKFRQVFTF48GDeeecdyspq/xSZpdzc3PD397/qgQwPDw8uXrxYr9ePa0i01ly8eBEPDw+Lj5EWsGh0tNa8+OKLvPzyywwdOpTnnnuOvXv3XnP1Czc3Nx566CFGjRrF1KlTef7559m9ezdfffWV3eYL7tKlC99//z25ubk0a9YMgI4dO5KSkmKVIWrCOjw8POjYsaPF5SUBi0ZFa80TTzzB4sWLeeyxx/j444/5+uuvad26NV5eXtc9tmXLlqxevZp58+Yxe/ZsJk+ezIoVK+yyZlxQUBDff/898fHx9OrVCwBXV9c6NWZZVJ90QYhG5fXXX2fx4sXMnTuXTz/9FCcnJxITEy1OZEop/vznP/P222+zatUqnnjiCbt0R7Rr1w4PDw8SEhJsfi5hP9ICFo3GqlWreOGFF5g0aRKvvfYaSilSUlIoKSmpdkty9uzZZGdn88orr+Dn58crr7xio6iNnJycCAgI4MyZMzY9j7AvScCiQao87rXoXAxpK/5Kl57hfPrpp+UjB86cOYNSisDAwGqf46WXXiIlJYXXXnuNwYMHM2zYMGuFX6XAwEBOnTpFZmYmLVu2tOm5hH1IF4RokCquAHFLe0+a7/4I/w5+/Lpj8xWTwJw5cwY/P78arX6hlGLhwoWEhoYyadIk0tLSrPktXMXcSpduiIZDErBo8I7+70NOnz5NZGTkFfP8FhcXk5KSUqsbWU2aNOHrr78mJyeHiRMnYjAYrBFyldq0aUOTJk0kATcgkoBFg1VqKONETCy5/gMZNnsht/9u8BX7k5OTKSsrq1UCnr8thtFfJOB5x+N8//33tLnnaQLnbLTKZOeVmbtKzpw5I2N/GwhJwKLB+veGY2TSFPf23Uly73JVUkxISMDJyYlOnTrV+Bzmro67x0XQtns/Sn9dwZ5nb7HZY7+BgYFkZ2dz+fJlm9Qv7EsSsGiw1vwShZOrsb+3sPTqFSASExNp3769VaafVErR98HZFBUV8dxzz9W6vmsxt9ZlNETDIAlYNEgxMTEk7t0GhhLg6hUgiouLSU1NJSAgoNbnMi92mVjqxbA/L2TlV6vYtm1breutio+PD82aNZN+4AZCErBokGbNmkXJwTX4NnWpcgUIc/9vTYafVVZxscvTqgOdxz7NU089RUlJSa3rrkz6gRsWScCiwdm0aRObNm3in/94gS5+PvTs0ILnR3bHxfm3H/fExESUUrXq/zW7YrHL0jL8w4dy+vRpPvvss1rXXZXAwEDy8vLKVxEW9ZckYNGgGAwG/vznP+PbMZD3zwfz65lL/HrmEoFzNl4xOiEhIYH27dtfMSa4pm4L8sE8I6SHqxOjwoO5/fbbeemll8jLy6t1/ZWZu03q60rA4jfyJJxoUFauXMnJkyf573//y4QJ91VZpqSkhNTUVKvN6/vcsBDWHkolq6CER24LZNawEG73epNBgwbx7rvv8re//c0q5zHz8fGhadOmJCUlccstt1i1bmFf0gIWDUZpaSkvvfQSvXv35v77779mOXP/rzVuwMFvi11W7Oq47bbbGDt2LG+++abVuwqUUgQEBJCQkCD9wPWctIBFg7Fs2TJiY2NZt27ddZcOSkhIQCmFv79/rc95vbXWXn/9dcLCwnjrrbd48803a32uigICAoiOjiYrK0vmhajHJAGLBqG4uJiXX36Zfv36ce+99163bGJiIn5+flbp/73RWmsRERF8+OGH/PWvf8Xb27vW5zOr2A8sCbj+ki4I0SAsXbqUxMREXn755euukVZaWmq18b+WmDt3Lrm5ubz//vtWrbdNmzZ4eHjIjbh6ThKwqPcMBgNvv/02/fr1Y8SIEdctm5qaisFgsFsCDgsLY+zYsbz77rvk5ORYrV5zP7Ak4PpNErCo99auXUtsbCx/+ctfbrhCsDlhWWP8r6X+/ve/c/nyZRYtWmTVegMCArh06ZJVE7uwL0nAol7TWvPmm28SHBzM73//+xuWT0pKwtfXlyZNmtghOqP+/fszdOhQ/vOf/1BYWGi1emU8cP1n0wSslBqplDqllIpVSs2pYr9SSr1n2n9UKdXXtN1DKbVXKXVEKRWllHrJlnGK+mvnzp3s27eP2bNn33BxzLKyMpKTk+3W/VDR3LlzSUtL44svvrBane3atcPNzU0ScD1mswSslHIGPgDuAW4CHlRK3VSp2D1AV9NrGvCRaXsRcJfWujfQBxiplBpgq1hF/fXWW2/Rpk0bpkyZcsOyaWlpFBcXW2X4WXXdeeed9OnTh3fffddqY3fNU2kmJSVZpT5hf7ZsAfcHYrXW8VrrYuArYFylMuOAZdpoD9BSKeVnep9rKuNqesmIc3GF6OhoNm3axNNPP23RkkLmROWIFrBSimeffZaoqCi2b99utXr9/f1JT0+noKDAanUK+7FlAu4AJFd4n2LaZlEZpZSzUuowkA5s01r/WtVJlFLTlFL7lVL7L1y4YK3YRT3w/vvv4+7uzvTp0y0qbx4z6+XlZePIqhYREUHbtm1ZsGCB1eo0/zGRVnD9ZMsEXNXt6Mqt2GuW0VobtNZ9gI5Af6VUz6pOorX+RGsdrrUO9/X1rU28oh7JzMxk2bJlPPTQQ1es83YtWmuSkpIc0vo1c3d3Z8aMGWzcuJGYGOssWdShQwecnZ0lAddTtkzAKUDFsT4dgbPVLaO1zgR+AEZaPUJRb33++efk5+fz9NNPW1Q+IyOD/Px8h/T/VjR9+nTc3Nx47733rFKfi4sL7du3lwRcT9kyAe8DuiqlOiul3IAIYH2lMuuBKabREAOALK31OaWUr1KqJYBSyhMYCpy0YayiHjEYDCxcuJBBgwZx8803W3SMI/t/K2rTpg0TJ05k6dKlZGVlWaVOf39/zp49a5MJ4IVt2SwBa61LgZnAt8AJ4GutdZRSarpSytxptwmIB2KBxcAM03Y/YIdS6ijGRL5Na73BVrGK+mXz5s3Ex8db3PoFYwJu2rSpVedjqKkZM2aQl5fH8uXLrVKfv78/ZWVlpKamWqU+YT82nYxHa70JY5KtuO3jCl9r4KkqjjsKWNa0EY3OwoULad++/XWnnKwsKSkJf3//Gz4pZw/h4eGEh4fz0UcfMWPGjFrHZH6qLzEx0SpLLAn7kSfhRL0SHx/Pt99+y7Rp03B1dbXomOzsbDIzMx3e/1vRk08+SVRUFLt27ap1XZ6enrRp00b6geshScCiXlm8eDFOTk5MnTrV4mPMiakuJeCIiAhatGjBxx9/fOPCFvD39yclJYWysjKr1CfsQxKwqDeKi4tZsmQJY8aMoWPHjhYfl5SUhJubG+3atbNhdNXTpEkTHn74YVavXo01xq8HBARQXFxMWlqaFaIT9iIJWNQb69evJz09nWnTplXruKSkJDp27HjdVTIcYfr06eV/VGrL3LqXeSHql7r1EynEdXzyySd06tSJkSMtHxJeWFjI+fPn61T3g1loaCh33HEHn376aa3nh/Dy8qJly5YkJyffuLCoMyQBi3ohLi6Obdu28fjjj99w1rOKzAmpLiZggKlTpxIbG8vOnTtrXZe/vz9JSUmyUGc9IglY1AufffZZtW++gbH7wcnJiQ4dKk9DUjeMHz+eFi1a8Omnn9a6Ln9/f/Ly8rh06ZIVIhP2IAlY1HmlpaVERkYyatSoaifSpKQk/Pz8cHNzs1F0tdOkSRMeeughVq9eTWZmZq3qMrfyZTha/XHDBKyUClFKbVdKHTe976WUesH2oQlhtHXrVs6ePctjjz1WrePMC3DW1e4Hs6lTp1JYWMiXX35Zq3pat26Np6enJOB6xJIW8GJgLlAC5U+pRdgyKCEqWrJkCb6+vowePbpax509exaDwVDnE3Dfvn3p06cPn332Wa3qUUqV9wOL+sGSBNxEa7230rZSWwQjRGUXLlxg/fr1TJ48udrdCOZEZM8FOGtCKcXUqVM5ePAghw4dqlVdnTp14tKlS+Tm5t64sHA4SxJwhlIqCNM8vUqpCcA5m0YlhMmKFSsoKSnh0UcfrfaxycnJtG7dmqZNm9ogMuuaOHEibm5uLF26tFb1yATt9YslCfgpYBHQXSmVCjwLWLYEgRC1oLVmyZIl9O/fn549q5yP/7rHJiUl1fnWr1mrVq247777WLFiBcXFxTWux8/PDxcXF0nA9YQlCVhrrYcCvkB3rfXtFh4nRK0cPHiQY8eO1aj1e+HCBQoLCx0+/291PPLII1y8eJGNGzfWuA5nZ2c6dOggD2TUE5Yk0m8AtNZ5Wusc07bVtgtJCKOlS5fi7u5ORET17/maH8mt6zfgKho2bBh+fn617obw9/fn3LlzFBUVWScwYTPXnA9YKdUd6AG0UEpVnHjVC/CwdWCicSsuLmblypWMGzeOli1bVvv45ORkmjdvXqNjHcXFxYVJkyYxb948zp8/T9u2bWtUj7+/P1prUlJSCAoKsnKUwpqu1wLuBowBWgL3Vnj1Bf5k88hEo7Zp0yYuXrzIww8/XKPjExMT68wE7NXx8MMPYzAYajUmuFOnTiilpB+4HrhmC1hrvQ5Yp5QaqLXebceYhCAyMpJ27doxfPjwah+blZVFdnZ2vbkBV1GPHj3o168fS5cuZdasWTWqw93dnbZt20oCrgcs6QM+pJR6Sin1oVJqifll88hEo5WRkcHGjRuZOHEiLi7VXzXL3P9bn27AVfTII49w9OhRjhw5UuM6zBO0GwwGK0YmrM2SBPwF0A4YAfyIcen4nOseIUQtrFy5kpKSkhp3P5gnYG/Tpo2VI7OPBx54AFdXV5YtW1bjOgICAigtLeXcORmyX5dZkoCDtdb/APK01pHAaCDMtmGJxiwyMpI+ffoQFlazHzPz+N+6NgG7pXx8fBg9ejRffvklpaU1e+hUJuapHyz5CS0x/ZuplOoJtAACbRaRaNROnDjBgQMHmDJlSo2Oz8/P58KFC/W2+8FsypQppKWl8d1339Xo+GbNmuHt7S0JuI6zJAF/opRqBbwArAeigTdtGpVotL744gucnZ158MEHa3R8XVyAsyZGjRpFq1atatUNIRO0133XTcBKKScgW2t9WWu9U2vdRWvdRmu9yE7xiUakrKyM5cuXM3z48BovoJmUlFT+NFh9Zn4AZe3atWRnZ9eoDn9/fwoKCsjIyLBydMJarpuAtdZlwEw7xSIauR9//JHk5GQmT55c4zqSkpLo0KFDjUZP1DVTpkyhoKCAb775pkbHy0KddZ8lXRDblFKzlVKdlFLe5pfNIxONzhdffEHz5s0ZN25cjY4vLi7m7Nmz9b77wezWW2+la9euNe6G8Pb2pmnTptIPXIdZkoAfwzgj2k7ggOm135ZBicYnPz+f1atXM378eJo0aVKjOlJSUtBa1/sbcGZKKSZNmlT+yaAmxwcEBEgCrsNumIC11p2reHWxR3Ci8Vi/fj05OTm16n5ITExEKVUvn4C7lokTJ6K1rvGjyf7+/mRlZdV6vTlhG/VzoKRocJYvX07Hjh0ZMmRIjetISkqiXbt2uLu7Wy8wBwsKCuK2227jiy++qNFoBvOnAekHrpskAQuHS09PZ8uWLUycOLHGD08YDAZSUlIaTP9vRZMmTSIqKqpGjya3bdsWDw8PScB1lCRg4XCrVq3CYDAwadKkGtdx9uxZSktLG0z/b0V//OMfcXV1Zfny5dU+1rxQpyTgusmSZem3W7JNiJpavnw5ffr0qfayQxXVxwnYLeXj48OoUaP48ssvazS5TkBAAJcuXSInR6ZwqWuumYCVUh6m4WatlVKtKgxBCwTa2y1C0aDFxMSwd+/eWrV+wZiA68sCnDUxefJkzp07x/bt1W/7yEKdddf1WsBPYBxy1p3fhp8dANYBH9g+NNEYLF++HCcnpxo/egzGJ+iSkpIaZPeD2ejRo2nZsmWNuiHatWuHq6urdEPUQddMwFrrd7XWnYHZpkeQzUPQemutF9oxRtFAaa1Zvnw5d999N+3b1/xDVVpaGsXFxQQGBlovuDrGw8ODP/zhD/zvf/8jLy+vWsc6OzvTqVMnScB10A2f19Rav6+Uug3jDGguFbbXfJYQ0SjN3xbDu9tPl78vTDnB+TNnGDD+iVrVm5CQANTfCdgtNWnSJBYvXszatWuZOHFitY4NCAhgx44dFBQU4OnpaaMIRXVZchPuC+Ad4Hagn+kVbuO4RAM0a1gICW+M5tbO3tza2ZvfN4/D09OTRf+cUat6ExMT8fb2pnnz5laKtG66/fbbCQgIqFE3hIwHrpssmbEkHLhJy5x2wgpKDWUkXconO7+YS/Fl3Pf7+2uVOM39v6GhoVaMsm56d3ssl9vfypZvV9Np5nKcm7UC4Jm7uzJrWMh1jzVPUJSYmEj37t3tEa6wgCUJ+DjGJYlkbRNRa/O2xXA+u5AyDa49h+MVVLuPw+fPn6ewsLBB9/+azRoWwsiOL3LTTV/TOn0/hxYusPhYFxcXOnbsWN5dI+oGSx7EaA1EK6W+VUqtN78sqVwpNVIpdUopFauUmlPFfqWUes+0/6hSqq9peyel1A6l1AmlVJRS6pnqfVuirvol7iJlps9STq4enCvzqlV99X0BzurqGtKNdmNmccErhLe2nKTUUGbxsYGBgaSlpVFQUGDDCEV1WNICfrEmFSulnDEOVxsGpAD7lFLrtdbRFYrdA3Q1vW4FPjL9Wwr8WWt9UCnVHDiglNpW6VhRD90W5MOR5Ew04KQN3B7culb1JSYm0qpVK1q0aGGdAOu4edticA+9A5xc+OyneBTw/EjLuhTMnxKkG6LusGQ2tB+rellQd38gVmsdr7UuBr4CKk/0Og5Ypo32AC2VUn5a63Na64Om8+cAJ4D6vcSBAOC5YSF4ZCdSdPYk94V63bDv8nq01iQmJjaa1i8YP0HgZGw3FRk0P8ddtPhYcz+wdEPUHZaMgshRSmWbXoVKKYNSypI1UjoAFScxTeHqJHrDMqYn724Gfr1GfNOUUvuVUvsvXLhgQVjCUeZviyH475s58/8+4OKmd/nmRC7Bf9/M/G0xNaovPT2dgoKCRpWAbwvyQSnj17q0iNu6WL42gouLC506dZIEXIdY0gJurrX2Mr08gPGAJQ9iqKqqq04ZpVQz4BvgWa11lUlfa/2J1jpcax3u6+trQVjCUWYNC+H7aaEUpUTz0p+fJPHNMSS8MbrGreAzZ84A0LlzZ2uGWac9NywEPy8PXEtyydq7lr5u1bs3HhAQwPnz56UfuI6o9mxoWuu1wF0WFE0BKs6M3RE4a2kZpZQrxuS7Qmv9v+rGKeqm5cuXo5Sq9oMEVUlISMDb27vR9P+aP0GczSqkCFeyf13Nfc/+u1qfIMx/rGQ8cN1ww5twSqn7K7x1wjgu2JIxwfuArkqpzkAqEAE8VKnMemCmUuorjDffsrTW55RSCvgMOKG1nmfBuUQ9oLVm2bJl3HnnnbWetaysrIyEhAR69OhhpejqvlnDQq74tDA180G+/vprnhjU0eI62rdvj4uLC2fOnJEbcXWAJS3geyu8RgA5XH0z7Spa61KMKyp/i/Em2tda6yil1HSl1HRTsU1APBALLAbMj0QNAiYDdymlDpteoyz/tkRdtGfPHuLi4mq17JDZuXPnKCoqalTdD5VNnjyZ3Nxc1q1bZ/ExLi4uMj9wHWLJXBCP1rRyrfUmjEm24raPK3ytMS74Wfm4XVTdPyzqsWXLluHp6cn48eNrXZe5/7cxPIBxLYMHD8bf359ly5ZVazY587wQ+fn5NV4AVViHJaMgOiql1iil0pVS55VS3yilLP/MIwRQVFTEV199xf331+7RY7OEhAR8fX1p1qyZFaKrn5ycnJg0aRJbt27l3DnLb8aZPzWY/4gJx7GkC+JzjH217TEOEft/pm1CWGzDhg1kZmZapfvBYDCQlJTUqLsfzKZMmUJZWRkrVqyw+JgOHTrg5uYmCbgOsCQB+2qtP9dal5peSwEZ7yWqJTIyEj8/P+6+++5a15WamkpJSYkkYKBbt24MGDCAyMhIi1dNdnJyIjAwkPj4eBtHJ27EkgScoZSapJRyNr0mAZY/fiMavfPnz7Np0yYmT56Mi4slT79fn7nl1pgewLiehx9+mOPHj3Po0CGLj+nSpQuXL18mMzPTdoGJG7IkAT8G/BFIwzgj2gTTNiEsYl5M8uGHH7ZKfWfOnMHPz08mFjd54IEHcHd3JzIy0uJjzJ8epBXsWJY8CZektR6rtfbVWrfRWt+ntZYxLMIiWms+//xz+vXrx0033VTr+oqLi0lOTpbuhwpatWrF2LFj+fLLLykuLrboGPMNTOkHdixLRkFEKqVaVnjfSim1xKZRiQbj8OHDHDt2jEceecQq9SUmJlJWVkZQUJBV6msoHn74YTIyMti8ebNF5ZVSdOnShfj4eIv7joX1WdIF0UtrnWl+o7W+jHFyHCFuKDIyEjc3NyIiIqxSX1xcXPnDBOI3I0aMoF27dnz+ueUDlDp37kx+fj7p6ek2jExcjyUJ2Ekp1cr8RinljWXzCItGrri4mBUrVjBu3Di8vS2ftet64uLiCAgIsMrNvIbExcWFKVOmsGHDBs6fP2/RMdIP7HiWJOD/AL8opV5RSr0M/AK8ZduwREOwceNGMjIyrHbzLTs7m4yMDLp06WKV+hqaRx99FIPBYPGinS1atMDHx0cSsANZchNuGcYpKM8DF4D7tdZf2DowUf999tln+Pn5MWLECKvUFxcXByD9v9fQvXt3Bg4cyJIlSyzu1+3SpQuJiYmUlpbaODpRFYumo9RaR2utF2qt35dlgYQlzp49y+bNm3nkkUes1l0QHx9P06ZNadOmjVXqa4geffRRoqOj2bdvn0Xlg4ODKSkpISkpycaRiapUez5gISwRGRlJWVkZjz1mnSHjWmvi4+MJCgpCKZmn6VoeeOABPD09WbLEsoFKgYGBODs7Exsba+PIRFUkAQur01qzZMkS7rjjDoKDg61SZ1paGvn5+dL/ewNeXl5MmDCBlStXkp+ff8Pybm5u+Pv7l3fvCPuSBCysbufOncTGxlqt9Qu/9f9KAr6xqVOnkp2dzX//+1+LygcHB5Oenk52tiVLPQprkgQsrO6zzz4rb4lZS2xsLO3atbPKVJYN3eDBgwkJCWHx4sUWlTd/SpFuCPuTBCysKjMzk9WrV/Pggw9abbLvwsJCkpKSrNad0dAppXj88cf5+eefiY6+8T1zX19fmjdvLt0QDiAJWFjV8uXLKSgoYNq0aVarMy4uDq01Xbt2tVqdDd3DDz+Mq6srn3766Q3LKqUIDg4mLi6OsrIyO0QnzCQBC6vRWrNo0SLCw8Pp27ev1eo9ffo0np6edOwoC7FYqk2bNowbN45ly5ZRVFR0w/LBwcEUFRWRmppqh+iEmSRgYTV79uzh+PHjVm39aq2JjY0lKCgIJyf5ca2OP/3pT1y8eJE1a9bcsGyXLl1QSnH69Gk7RCbM5CdaWM2iRYto1qxZtRaIvJGzZ8+Sl5cn3Q81MHToUAIDA/nkk09uWNbDwwN/f39iYmLsEJkwkwQsrOLy5cusWrWKiRMnWnWhTHOLTB4/rj4nJyemTZvGjh07OHny5A3Lh4SEcP78eVklw44kAQurWL58OYWFhTzxxBNWrTc2NpaOHTvStGlTq9bbWEydOhVXV1c+/vjjG5bt1q0bgLSC7UgSsKg1rTUfffQR/fr14+abrTdVdF5eHqmpqTL8rBbatGnDhAkTWLp0KXl5edct6+Pjg4+PjyRgO5IELGptx44dnDhxgpkzZ1q1XnMiCAkJsWq9jc2MGTPIysriq6++umHZkJAQEhISLBo5IWpPErCotYULF9K6dWv++Mc/WrXekydP0qJFC9q1a2fVehubQYMG0bNnTz788MMbTlMZEhKCwWCQhzLsRBKwqJWkpCTWrVvH448/joeHh9XqLS4uJi4uju7du8vsZ7WklGLGjBkcPHiQvXv3Xresv78/Hh4e0g1hJ5KARa0sWrQIgOnTp1u13tjYWAwGA927d7dqvY3VpEmTaN68Oe+///51yzk5OdG1a1dOnz4tT8XZgSRgUWNFRUUsXryYe++9l4CAAKvWffLkSTw9PWXxTStp3rw5jz32GF9//TXnzp27btlu3bqRn59PcnKynaJrvCQBixpbtWoVFy5c4KmnnrJqvQaDgZiYGLp16yZPv1nR008/TWlpKR999NF1ywUHB+Pi4mLRRD6iduSnW9SI1pp58+bRo0cPhg4datW6zXfhpfvBuoKCghgzZgwff/wxhYWF1yzn7u5OcHAwJ06csHhtOVEzkoBFjfzwww8cOXKEWbNmWf0m2cmTJ3F1dZXJ123gmWee4cKFCzcckhYaGkpOTg4pKSl2iqxxkgQsamTevHn4+voyceJEq9ZbVlbGyZMnCQoKwtXV1ap1C7jrrrvo2bMnCxYsuG7rNiQkBGdnZ+mGsDFJwKLaYmJi2LBhAzNmzLDq0DMwDmvLzc2lR48eVq1XGCmleOaZZzhy5Ag7duy4ZjkPDw+CgoKkG8LGJAGLaluwYAHu7u48+eSTVq/7+PHjuLq6ytNvNjRp0iTatm3LW2+9dd1yoaGhZGVlcfbsWTtF1vhIAhbVcuHCBZYuXcrEiRNp27atVes2GAycOHGCbt264ebmZtW6xW88PDx45pln+Pbbbzly5Mg1y5lHoUg3hO1IAhbV8t5771FYWMjs2bOtXveZM2fIz8+X7gc7mD59Os2aNePtt9++ZhlPT0+6dOlCdHS0dEPYiCRgYbHs7GwWLlzIfffdR2hoqNXrj4qKKh8CJWyrVatWTJs2ja+++orExMRrluvZsyeZmZnyUIaN2DQBK6VGKqVOKaVilVJzqtivlFLvmfYfVUr1rbBviVIqXSl13JYxCsstWrSIzMxM5s6da/W6S0tLOXHiBKGhobi4uFi9fnG1Z599FqUU8+bNu2aZ0NBQXF1dOXr0qB0jazxsloCVUs7AB8A9wE3Ag0qpmyoVuwfoanpNAyo+orMUGGmr+ET1FBYWMm/ePO6++2769etn9fpjY2MpKiqS7gc76tSpExMnTmTx4sWkp6dXWcbNzY3u3bsTFRVFaWmpnSNs+GzZAu4PxGqt47XWxcBXwLhKZcYBy7TRHqClUsoPQGu9E7hkw/hENURGRpKWlmaT1i/AsWPHaNKkCZ07d7ZJ/aJqc+fOpaioiHfeeeeaZcLCwigsLCQ2NtaOkTUOtkzAHYCKHUcppm3VLSMcrLi4mDfeeIP+/ftz1113Wb3+/Px8Tp06RVhYGM7OzlavX1xbt27diIiI4IMPPuDChQtVlgkKCqJp06bSDWEDtkzAVT2fWvlWqiVlrn8SpaYppfYrpfZf6wdI1E5kZCQJCQm8+OKLNpmb99ixYxgMBqsuZyQs98ILL1BQUMD8+fOr3O/k5ETPnj2JiYmhoKDAztE1bLZMwClApwrvOwKVR3RbUua6tNafaK3Dtdbhvr6+NQpUXFtxcTGvvvoqAwYMYORI23TJHz58GD8/P6uPKxaWCQ0N5Y9//CPvv/8+ly5V3evXq1cvDAaDjAm2Mlsm4H1AV6VUZ6WUGxABrK9UZj0wxTQaYgCQpbW+/mSlwq6WLFlCUlISL730kk1av2lpaaSlpdGnTx+r1y0s98ILL5Cbm8t//vOfKvf7+fnh6+vLwYMH7RxZw2azBKy1LgVmAt8CJ4CvtdZRSqnpSinz8gmbgHggFlgMzDAfr5RaCewGuimlUpRSU20Vq6haUVERr732GoMGDWLYsGE2OcehQ4dwdnYmLCzMJvULy/Ts2ZOIiAgWLFhAWlraVfuVUtxyyy2cPXv2hhO6C8vZdByw1nqT1jpEax2ktX7NtO1jrfXHpq+11vop0/4wrfX+Csc+qLX201q7aq07aq0/s2Ws4mqLFi0iJSXFZq1fg8HAsWPH6NatG56enlavX1TPK6+8Ut7lVJVevXrh4uLCgQMH7BxZwyVPwokqZWVl8fLLLzN06FCbjHwAOHHiBAUFBdL9UEcEBwczdepUPvnkE+Lj46/a7+npSY8ePTh27BjFxcUOiLDhkQQsqvTmm29y8eJF3nrrLZutSrx3715atWoljx7XIf/85z9xdnbmX//6V5X7+/btS3FxMcePywOq1iAJWFwlJSWF+fPnM2nSJJsNDTt37hzJycn069dPlp2vQ9q3b88zzzzDihUrOHz48FX7O3XqhK+vr3RDWIkkYHGVf/7zn5SVlfHKK6/Y7Bx79+7F1dVVxv7WQXPmzMHHx4dnn332qlnQKt6Mk3mCa08SsLjCwYMHWbp0KTNnziQwMNAm58jPz+f48eP06tXL6itqiNpr2bIlr7zyCj/++CPffPPNVft79+6Nm5sbe/bscUB0DYskYFGurKyMmTNn4uvryz/+8Q+bnefQoUOUlpbaZFIfYR2PP/44YWFhPP/881etoOzh4UHfvn05fvw4WVlZDoqwYZAE3MjN3xZD4JyNBM7ZSMCzX3HCtSuG8Il8vq/q2bFqy2AwsG/fPgIDA+XJtzrMxcWFBQsWkJCQUOV0lQMGDADg119/tXdoDYok4EZu1rAQZgwJQgHOTVrQov/v+cvMPzFrmG3WZDO3msy/wKJumr8thse2FtD6vrm8f8oT/2e/JnDORuZviwGgRYsW9OjRg4MHD1JUVOTgaOsvScCCX+Iu/jYDkrMrv8TbZhZQrTW7du2ibdu2suhmHWf+w9ys2224+4Xg7ObGjCFBV/xhHjhwIEVFRfJ4ci1IAhYEeBRSVmLs5/NwdWJQkI9NznPy5EkyMjK4/fbbZehZPWD8w2z8f9JOrvy/vTFX7G/fvj0BAQHs2bNHJmuvIUnAjdzbG4/z0ZzHyT2yDV1moLCkjA9+iCv/qGktWmt++uknvL29uemmygujiLrotiAfyv9OGoo5s2cLmZmZV5T53e9+R3Z2trSCa0gW32rkcvf+l5KMBL6YdS+jR4+22Xni4uI4d+4c9957L05O8ne/rpu/LYYPf4grf19akEvqt4sZEXGeX7f8NjStS5cu+Pv7s2vXLm6++WZcXV0dEW69JQm4ETt69Civv/46EydOtGny1VqzY8cOvLy86N27t83OI6xn1rCQq27E/s3rBP/+979Zv349Y8eOBYwPZtx5551ERkZy4MABublaTdIUaaQKCgp46KGH8Pb2ZsGCBTY9V3R0NGfPnuXOO++UJYfqsRdffJHevXvzpz/96YrliwIDA+ncuTO7du2SSXqqSRJwI/WXv/yFqKgoIiMjad26tc3OYzAY2L59O23atKFXr142O4+wPTc3N5YvX05mZiZPPPHEFY8pDxkyhLy8PPbu3evACOsfScCN0IYNG1i4cCHPPvuszZYZMtu/fz+XL19m6NCh0vfbAPTs2ZPXXnuNNWvWsGjRovLt/v7+hISE8NNPP5GTk+PACOsX+Y1oZFJSUnjsscfo1asX//73v216rsLCQnbu3Ennzp1lyskG5LnnnmPEiBE888wzV8yKNmLEiPJPPMIykoAbkaKiIiZMmEBBQQFfffWVzSfC2bFjB/n5+QwdOlTG/TYgTk5OLF++nDZt2vCHP/yBy5cvA+Dt7c2AAQM4cuQIKSkpDo6yfpAE3Ig8/fTT/PrrryxdupTQ0FCbnis1NZW9e/fSr18/2rdvb9NzCftr3bo1//3vf0lOTmbKlCmUlZUBxnHBzZo1Y/PmzVdNZSmuJgm4kfjkk09YvHgxc+fOZfz48TY9V1lZGRs2bKB58+Y2W85ION6AAQNYsGABGzZsYM6cOQC4u7szbNgwzp49y759+xwcYd0nCbgR+Pbbb3nqqacYMWKETSdZN9uzZw9paWmMHDlS5vtt4GbMmMFTTz3F22+/zaeffgpAWFgYQUFBfPfdd1y6ZJt5RRoKScAN3MGDBxk/fjw9e/bk66+/tvk43AsXLrBjxw5CQkJs3s0hHE8pxYIFCxg5ciRPPvkk3333HUopxo4di5OTE+vWrSvvnhBXkwTcgJ05c4ZRo0bh4+PDxo0b8fLysun5SktL+eabb3Bzc2PMmDFy462RcHFxYdWqVYSGhnLfffexe/duvLy8GDlyJElJSTJn8HVIAm6gEhISuPPOOykuLmbz5s12uRG2detWzp8/z3333Ufz5s1tfj5Rd3h5ebF161b8/Py45557OHToEL1796Zbt25s376d1NRUR4dYJ0kCboASEhIYMmQIWVlZfPfdd3aZfezEiRPs27ePAQMG0LVrV5ufT9Q97dq1Y/v27bRo0YLhw4dz9OhRxo4dS7Nmzfj666/Jy8tzdIh1jiTgBiY2NpY777yzPPn27dvX5udMS0tj7dq1+Pn5cffdd9v8fKLu8vf3Z/v27Xh4eDBkyBAOHz7MAw88QH5+PqtXr5b+4EokATcge/fuZeDAgeTk5LBt2zZuueUWm58zJyeHlStX4u7uTkREBC4uMsFeYxccHMyuXbvw9fVl6NChHD16lNGjR5OQkMCWLVtkfHAFkoAbiA0bNnDnnXfSvHlzfvnlF8LDw21+zuLiYlauXFk+s5qtb/KJ+iMgIICffvqJbt26MWbMGHbv3s2AAQPYt28fO3fudHR4dYYk4HqurKyMl19+mbFjx9K9e3d2795tl/XWiouLWbFiBWlpaUyYMIF27drZ/Jyifmnbti0//vgjI0aMYMaMGaxZs4YePXrwww8/yMgIE0nA9dilS5e49957+de//sXEiRP56aef7LLUe1FREcuXLyc5OZn7779fFtgU1+Tl5cW6dev461//yqJFi5g3bx6tWrViy5YtMnUlkoDrrS1bthAWFsa2bdv44IMPWLZsGU2aNLH5efPy8li+fDmpqalMmDCBnj172vycon5zdnbmjTfeYNWqVURFRfGPf/yDjIwMNm/ezPfff9+o+4Tljkk9k5mZyZw5c1i0aBE33XQT69evt8vNNoDz58+zcuVK8vLymDBhgjzpJqrlj3/8I/3792fixIksXLiQW2+9laysLHJychg9enSjvIHb+L7jeqLUUMa8bTH8EneR24J8ePbuYFYs/4I5c+Zw4cIFZs+ezSuvvGK3uRaioqJYt24dHh4ePProozLDmaiRwMBAfvzxR/7zn//w0ksvcfToUaKjo0lNTSUiIgJvb29Hh2hXqiE1/8PDw/X+/fsdHUatzd8Ww7vbT5e/12UGco9/z6XN7zJw4EAWLlxol/G9YJxUffPmzRw9epQOHTrwwAMPyFNuwiri4+N56qmn2LJlCz6tfek24c807dKXu3t25M/Du+Hi3HB6SJVSB7TWVw1NkgRcR933wc8cTs4sf68zzvDKYC8mTZpkl6V9tNacOHGCLVu2kJuby+DBg/nd734ni2oKq9Ja88Trn/K/XUdpctMQnFw9UGWluKgyZgztcdXKzPXVtRKwdEHUMdnZ2Uz9x3x+OpVB0x7GH0hdZgBvfy76dbdL8k1NTWXr1q0kJSXRpk0bIiIipMtB2IRSik/+/ifOL9zFkZQsALSTCwXnY4latZEzwX+nc+fODo7SdiQB1wElJSX88MMPLFu2jG+++YaCggJuviWcnn53kOXuxe1dfZk1LMSmH8m01sTHx/PLL78QHx9P06ZNGTNmDDfffLMspilsav62mPLkC8YuN50Rz2cbPmPp0qUMHDiQP/3pT0yYMMEuI33sSbogHOTixYts376dDRs2sGHDBi5fvkyLFi2IiIjgkUce4dZbb7XLdI7Z2dkcPXqUI0eOkJGRQbNmzbj11lvp168f7u7uNj+/EGC86Tx/Www/x11kUJAPs4aFsPfXPfz73/9mx44d5OXl4e7uzpAhQxg/fjyjR4+uV5/KpA/YgbTWJCYmsmfPHnbv3s2uXbs4dOgQWmtatWrFvffey/3338/w4cPx9PS0aSxlZWWkp6dz+vRpTp8+TXJyMgCdOnWib9++9OzZs1EOBxJ118WLF/n8889ZvXo1UVFR5ObmAtC5c2fuuOMO7rrrLsLDwwkJCamz9ygckoCVUiOBdwFn4FOt9RuV9ivT/lFAPvCI1vqgJcdWxdYJuPLQsOcqdQuUlJSQlJREfHw8p06d4tSpUxw/fpzDhw+TmZkJQJMmTejXrx933303w4YNIzw83GYJT2tNZmYm58+fJy0tjdTUVJKTkykqKgKgffv2hISEEBYW1uiG/4j6x2AwcObMGTZt2sSWLVuIjY0lKSmp/Oe5SZMmhIaG0qtXL3r16kXXrl0JDg4mMDDwik9zN/o9tgW7J2CllDMQAwwDUoB9wINa6+gKZUYBT2NMwLcC72qtb7Xk2KpYMwFrrSkoKCAnJ4fs7GyysrL4bN8FvksyUKIVzhjoXHSGlsk/kZKSQmpqKmfPnr1iuj0vLy9CQ0Pp06cPvXv3pn///vTq1QtXV9dax1ZUVERhYSEFBQXk5eWRl5dXHmt2djaXLl3i8uXLlJaWlh/n6+uLv78/nTp1IigoiGbNmtUqDiEcRWvN+fPnOX36NL/88guHDh0iPj6e9PR00tPTKSgouKK8t7c3fn5+0KQVlwLvxjXgZpSzC7rMQGlWOn/o0YInR/TCy8uL5s2b06RJE6ve+3DEKIj+QKzWOt4UwFfAOKBiEh0HLNPGvwJ7lFItlVJ+QKAFx9ba+PHjiYqKoqSkhJKSEoqKiiguLqawsJCioqKrHpFsN/kd3Nt3B8CAM9GXyij9eTcd/doSGhrKHXfcQbt27WjXrh0dOnTA29sbpRRaa7TWpKens23btvL3WmvKysrK/y0rK8NgMFBWVkZpaSkGg4HS0lJKS0spKSmhuLi4/HUt7u7ueHl54e3tTXBwMD4+PrRt25Y2bdrg5uZmzcsnhMMopcp/1373u9+htSYnJ4fz58+Tnp5ObGwsp0+fJiEhgbNnz5KZmUl2djY5F2Jwuel+lLMx9SknZ8oKsnj3mT/xbqX63dzcyl/u7u64urrSrl079uzZY7Xvw5YJuAOQXOF9CsZW7o3KdLDwWACUUtOAaWCcDLo6srOzMRgMuLm54enpiaurKy4uLri6ul5x8T08PHB1dSOlBaRqA2XKGWcM9O3YnH5PTsfZ6cqbZUVFRcTHx3PmzBlzjAA4OTld9bVSCicnp/KXs7MzTk5OuLi44OzsjIeHR3lM5rjc3d3x8PAofzVt2pSmTZvSrFkzuXEmGiWlFF5eXnh5edG1a1cGDRpUvk9rTX5+Prm5ueTl5fHRz6msO5FNcRm4KE3fIB96vPCCMUHn5JCfn09BQUF5Q6ywsLC8kWbt7kJbJuCqbuFX7u+4VhlLjjVu1PoT4BMwdkFUJ8Bt27ZVp3gVd2rHNKindYRoiJRS5Y0UgH8HBNLmit/jUQ77PbZlAk4BOlV43xE4a2EZNwuOtTsXZyeeH9md5x0diBCixurS77Et0/4+oKtSqrNSyg2IANZXKrMemKKMBgBZWutzFh4rhBD1ms1awFrrUqXUTOBbjEPJlmito5RS0037PwY2YRwBEYtxGNqj1zvWVrEKIYQjyIMYQghhY9cahiZ3kIQQwkEkAQshhINIAhZCCAeRBCyEEA4iCVgIIRykQY2CUEpdABKreVhrIMMG4dS3GEDiqGsxQN2Ioy7EAHUjjprGEKC19q28sUEl4JpQSu2vanhIY4tB4qh7MdSVOOpCDHUlDmvHIF0QQgjhIJKAhRDCQSQBm2ZSc7C6EANIHBXVhRigbsRRF2KAuhGHVWNo9H3AQgjhKNICFkIIB5EELIQQDtLgErBSaolSKl0pdbzCtheVUqlKqcOm16gK++YqpWKVUqeUUiMqbL9FKXXMtO89ZV5LyAZxKKV8lFI7lFK5SqmFleqpcRzVjGGYUuqA6VwHlFJ3Oeha9K+w7YhS6vf2vhYV9vub/k9mO+haBCqlCips/9gR10Ip1UsptVspFWU6p4cDrsXECtsOK6XKlFJ97HktlFKuSqlI07lOKKXmVjimZjFUXCCyIbyAwUBf4HiFbS8Cs6soexNwBHAHOgNxgLNp315gIMblkTYD99gwjqbA7cB0YGGlfTWOo5ox3Ay0N33dE0i1Rgw1iKMJ4GL62g9Ir/DeLteiwv5vgP9WLGPnaxFYsZyDfi5cgKNAb9N7H0f8jlQ6LgyId8C1eAj4qsLPaQIQWJsYGlwLWGu9E7hkYfFxGC9okdb6DMaJ4fsr48rMXlrr3dp4dZcB99kqDq11ntZ6F1BYcXtt46hmDIe01uZln6IAD6WUuwOuRb7WutT01gPTWoD2vBam890HxGO8FuZtdr0W14nNntdiOHBUa33EdOxFrbXBwdfiQWAl2P1aaKCpUsoF8ASKgezaxNDgEvB1zFRKHTV95Ghl2na9VZlTqthuqziuxVZx3CiG8cAhrXWRDWO4ZhxKqVuVUlHAMWC6KSHb7VoopZoCfwVeqlTWET8XnZVSh5RSPyqlfmfjOKqKIQTQSqlvlVIHlVJ/sXEM14qjogcwJWAbxlFVDKuBPOAckAS8o7W+VJsYGksC/ggIAvpgvHj/MW2v9arMVorjWmwRx3VjUEr1AN4EnrBhDNeNQ2v9q9a6B9APmGvqc7TntXgJmK+1zq1U3t7X4hzgr7W+GXgO+FIp5WWjOK4VgwvG7rGJpn9/r5S620YxXC8OwPjHGcjXWpv7bO15LfoDBqA9xi7LPyulutQmBluuilxnaK3Pm79WSi0GNpjeXmtV5hTT15W32yqOa7F6HNeLQSnVEVgDTNFax9kqhhvFUaHMCaVUHsY+aXtei1uBCUqpt4CWQJlSqhBjn7DdroXpE0iR6esDSqk4jC1Se16LFOBHrXWGad8mjH2my60dww3iMIvgt9avOT57XYuHgC1a6xIgXSn1MxAO/FTTGBpFC9jUR2P2e8D813M9EGHq6+wMdAX2auPKzDlKqQGmu5lTgHU2jKNKtojjWjEopVoCG4G5WuufbRnDDeLobOpjQykVAHQDEux5LbTWv9NaB2qtA4EFwOta64UOuBa+Siln09ddMP58xtvzWmBcGLeXUqqJ6f/lDiDaEb8jSikn4A/AV+Ztdr4WScBdyqgpMAA4WasYLL1bWF9eGP86ngNKMP51nAp8gbE/8SjGpOtXofzfMY5+OEWFO5cY/7IdN+1biOmpQRvGkYDxZkCuqfxNtY2jOjEAL2Ds3zpc4dXG3tcCmIzxxtdh4CBwnzX+T6r7/1HhuBe5chSEPa/FeNO1OGK6Fvc64loAk0xxHAfecuDvyBBgTxX12Ot3pBnGUTFRQDTwfG1jkEeRhRDCQRpFF4QQQtRFkoCFEMJBJAELIYSDSAIWQggHkQQshBAOIglYCCEcRBKwEEI4iCRg0WAopZ5QSp1TV84bG1apjKdpYhtnK5zvB1VhDmnTtmeVUh8qpdyUUjvNT/UJURVJwKIh6QW8oLXuU+F1rFKZx4D/aa0NVjjfSoxzE1QUAazUWhcD2zHO3CVElSQBi4YkDOMjzNczEdNz+sq44sRJpdSnSqnjSqkVSqmhSqmflVKnlVL9zQcppSYppfaaWtWLTC3o1cAYpZS7uT6MM2XtMh221nQ+IaokCVg0JD2Azyt0P0yruFMp5QZ00VonVNgcDLyLsfXcHeOMV7cDs4G/mY4LxdiSHaS17oNxSsKJWuuLGFdCGGmqKwJYpX97vv84xik1haiS9E+JBkEp1QlI11r3uk6x1kBmpW1nzN0UyjgJ/HattVZKHcO4JBDA3cAtwD7jZFd4YlwqCX7rhlhn+vcxc8XauHJEsVKqudY6pxbfnmigJAGLhqIXcPIGZQowLnNUUVGFr8sqvC/jt98PBURqredytbXAPKVUX8BTa32w0n53Ki01JYSZdEGIhiKMGyRgrfVlwNm0wkZ1bMc4QXsbAKWUt2muYrRx1YwfgCVcOVE4Sikf4II2TuAtxFUkAYuGIgyYXKH/95BSqlkV5bZi7OO1mNY6GuN8yVuVUkeBbRhXbDZbCfSmwkThJncCm6pzLtG4yHzAolFRSt0MPKe1nmyHc/0P4wojp2x9LlE/SQtYNCpa60PADms8iHE9phEXayX5iuuRFrAQQjiItICFEMJBJAELIYSDSAIWQggHkQQshBAOIglYCCEcRBKwEEI4yP8HZgUhDhuOVGAAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 360x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "x_arr = np.linspace(1500, 1580, 101)\n", + "\n", + "plt.figure(figsize=(5, 4))\n", + "plt.xlabel(r'$E$ (meV)')\n", + "plt.ylabel('count rate')\n", + "plt.errorbar(x, y, yerr=y_error, fmt='.', ms=7, capsize=3, label='sample')\n", + "plt.plot(x_arr, model_function(x_arr,[1545, 9]), '-', color='grey', label='initial_values')\n", + "plt.plot(x_arr, model_function(x_arr,[1540.1525, 11.1390]), '-', color='black', label='fit')\n", + "plt.legend()\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "What happens if you reduce significantly the number of events or increase the number of bins ? (mind empty bins!)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "hide_input": false, + "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.8.5" + }, + "toc": { + "base_numbering": 1, + "nav_menu": {}, + "number_sections": true, + "sideBar": true, + "skip_h1_title": false, + "title_cell": "Table of Contents", + "title_sidebar": "Contents", + "toc_cell": false, + "toc_position": { + "height": "calc(100% - 180px)", + "left": "10px", + "top": "150px", + "width": "225.438px" + }, + "toc_section_display": true, + "toc_window_display": true + } + }, + "nbformat": 4, + "nbformat_minor": 2 +}