{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Fitting a mixture of templates and parametric models\n", "\n", "The class `iminuit.cost.Template` supports fitting a mixture of templates and parametric models. This is useful if some components have a simple shape like a Gaussian peak, while other components are complicated and need to be estimated from simulation or a control measurement.\n", "\n", "In this notebook, we demonstrate this usage. Our data consists of a Gaussian peak and exponential background. We fit the Gaussian peak with a parametric model and use a template to describe the exponential background." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "from iminuit import Minuit\n", "from iminuit.cost import Template, ExtendedBinnedNLL\n", "from numba_stats import norm, truncexpon\n", "import numpy as np\n", "from matplotlib import pyplot as plt" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "We generate the toy data." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAp7klEQVR4nO3deXRUZZ7G8aeyB0gqbNm6gwSGVUGQJUagsSFtRpaGkW6lzSAIDT2QMEKOsjSb7MthE0QYFgE9YLrpaRhGaGwI2iiEgDEoCkSUINiaoGOTsHQSktz5w2O1RdKShKrUW8n3c06dQ269de/vvgTr8X3fe6/NsixLAAAABvHxdAEAAAC3I6AAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIzj5+kCaqK8vFxffPGFQkJCZLPZPF0OAACoAsuydO3aNUVHR8vH54fHSLwyoHzxxReKiYnxdBkAAKAGLl++rB//+Mc/2MYrA0pISIikb08wNDTUw9UAAICqKCwsVExMjON7/Id4ZUD5blonNDSUgAIAgJepyvIMFskCAADjEFAAAIBxCCgAAMA4XrkGBQBgnrKyMt26dcvTZcCDfH195efn55JbgBBQAAB37fr16/r8889lWZanS4GHNWjQQFFRUQoICLir/RBQAAB3paysTJ9//rkaNGig5s2bcwPNesqyLJWUlOirr75Sbm6u2rRpc8ebsf0QAgoA4K7cunVLlmWpefPmCg4O9nQ58KDg4GD5+/vrs88+U0lJiYKCgmq8r2pHmyNHjmjw4MGKjo6WzWbTnj17nN63LEuzZ89WVFSUgoODlZCQoPPnzzu1+eabb5SUlKTQ0FCFhYVpzJgxun79eo1PAgDgeYycQNJdjZo47ae6H7hx44buv/9+rVu3rtL3ly1bpjVr1mjDhg3KzMxUw4YNlZiYqKKiIkebpKQkffTRRzp48KBef/11HTlyROPGjav5WQAAgDql2gHl0Ucf1YIFC/Rv//ZvFd6zLEurV6/WzJkzNWTIEHXu3FmvvPKKvvjiC8dIy9mzZ3XgwAFt3rxZcXFx6t27t9auXau0tDR98cUXd31CAABUxcMPP6xJkya5bf+jRo3S0KFD3bZ/T7h48aJsNptOnTrl9mO5dA1Kbm6u8vLylJCQ4Nhmt9sVFxenjIwMDR8+XBkZGQoLC1P37t0dbRISEuTj46PMzMxKg09xcbGKi4sdPxcWFrqybACAOzxvr+XjFdTu8eBWLr1RW15eniQpIiLCaXtERITjvby8PIWHhzu97+fnpyZNmjja3G7x4sWy2+2OF08yBgDURyUlJZ4uodZ4xZ1kp0+froKCAsfr8uXLni4JAFAHlJaWKiUlRXa7Xc2aNdOsWbMc93J59dVX1b17d4WEhCgyMlJPPvmkrly54vT5jz76SIMGDVJoaKhCQkLUp08fffrpp5Ue6+TJk2revLmWLl3q2LZgwQKFh4crJCREv/71rzVt2jR16dLF8f5300QLFy5UdHS02rVrJ0k6ffq0+vXrp+DgYDVt2lTjxo1zutiksumroUOHatSoUY6fW7ZsqUWLFmn06NEKCQlRixYttHHjRqfPnDhxQl27dlVQUJC6d++u7OzsKvft3XJpQImMjJQk5efnO23Pz893vBcZGVnhL7i0tFTffPONo83tAgMDHU8u5gnGAABX2b59u/z8/HTixAm98MILWrlypTZv3izp28un58+fr/fff1979uzRxYsXnb7g//rXv+onP/mJAgMDdfjwYWVlZWn06NEqLS2tcJzDhw/rZz/7mRYuXKipU6dKknbs2KGFCxdq6dKlysrKUosWLbR+/foKn01PT1dOTo7jwpIbN24oMTFRjRs31smTJ7Vr1y4dOnRIKSkp1T7/FStWOILHhAkTNH78eOXk5Ej69uZ7gwYNUseOHZWVlaXnn39ezz77bLWPUVMuXYMSGxuryMhIpaenOxJgYWGhMjMzNX78eElSfHy8rl69qqysLHXr1k3St39x5eXliouLc2U5ALxUy2n7Kmy7uGSgBypBXRcTE6NVq1bJZrOpXbt2On36tFatWqWxY8dq9OjRjnatWrXSmjVr1KNHD12/fl2NGjXSunXrZLfblZaWJn9/f0lS27ZtKxxj9+7deuqpp7R582Y98cQTju1r167VmDFj9PTTT0uSZs+erT//+c8VbrvRsGFDbd682XFn1k2bNqmoqEivvPKKGjZsKEl68cUXNXjwYC1durTCMosfMmDAAE2YMEGSNHXqVK1atUpvvvmm2rVrp507d6q8vFxbtmxRUFCQ7r33Xn3++eeO73N3q/YIyvXr13Xq1CnHCt7c3FydOnVKly5dks1m06RJk7RgwQLt3btXp0+f1lNPPaXo6GjHSuYOHTroX//1XzV27FidOHFCR48eVUpKioYPH67o6GhXnhsAAD/owQcfdLp/S3x8vM6fP6+ysjJlZWVp8ODBatGihUJCQtS3b19J0qVLlyRJp06dUp8+fRzhpDKZmZn65S9/qVdffdUpnEhSTk6Oevbs6bTt9p8lqVOnTk63jT979qzuv/9+RziRpF69eqm8vNwx+lFVnTt3dvzZZrM5zXKcPXtWnTt3drrZWnx8fLX2fzeqHVDeffddde3aVV27dpUkpaamqmvXrpo9e7YkacqUKZo4caLGjRvnSJoHDhxwOsEdO3aoffv26t+/vwYMGKDevXtXmPcCAMBTioqKlJiYqNDQUO3YsUMnT57U7t27Jf1joWpV7prbunVrtW/fXi+//HKNH6T4/SBSVT4+PhWei1TZ8W8PVzabTeXl5dU+njtUO6A8/PDDsiyrwmvbtm2Svj25efPmKS8vT0VFRTp06FCFIa8mTZpo586dunbtmgoKCvTyyy+rUaNGLjkhAACqKjMz0+nn48ePq02bNjp37pz+7//+T0uWLFGfPn3Uvn37CusnO3furLfffvsHg0ezZs10+PBhffLJJ3r88ced2rZr104nT550an/7z5Xp0KGD3n//fd24ccOx7ejRo/Lx8XEsom3evLm+/PJLx/tlZWX68MMP77jv24/zwQcfON1o9fjx49Xax93wiqt4AABwh0uXLik1NVU5OTl67bXXtHbtWj3zzDNq0aKFAgICtHbtWl24cEF79+7V/PnznT6bkpKiwsJCDR8+XO+++67Onz+vV199tcI0S3h4uA4fPqxz587pV7/6lWMR7cSJE7VlyxZt375d58+f14IFC/TBBx/c8ZEBSUlJCgoK0siRI/Xhhx/qzTff1MSJEzVixAjH+pN+/fpp37592rdvn86dO6fx48fr6tWr1eqbJ598UjabTWPHjtWZM2e0f/9+LV++vFr7uBsEFABAvfXUU0/p73//u3r27Knk5GQ988wzGjdunJo3b65t27Zp165d6tixo5YsWVLhy7lp06Y6fPiwrl+/rr59+6pbt27atGlTpWtSIiMjdfjwYZ0+fVpJSUkqKytTUlKSpk+frmeffVYPPPCAcnNzNWrUqDs+YK9BgwZ644039M0336hHjx76xS9+of79++vFF190tBk9erRGjhypp556Sn379lWrVq3005/+tFp906hRI/3v//6vTp8+ra5du2rGjBlOl0i7m826fZLKCxQWFsput6ugoIBLjoE6iKt4vEtRUZFyc3MVGxt7V0+vhfSzn/1MkZGRevXVVz1dSo390O9Ddb6/XXqZMQAAqJqbN29qw4YNSkxMlK+vr1577TUdOnRIBw8e9HRpRiCgAADgATabTfv379fChQtVVFSkdu3a6b//+7+dnmdXnxFQAADwgODgYB06dMjTZRiLRbIAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAcJtRo0Zp6NChtX7c559/Xl26dKn145qI+6AAANyiskcWuJMrH4fwwgsvyAufBFOnEFAAALiN3W73dAn1HlM8AIB66w9/+IM6deqk4OBgNW3aVAkJCbpx40aFKZ5r164pKSlJDRs2VFRUlFatWqWHH35YkyZNcrRp2bKlFi1apNGjRyskJEQtWrTQxo0bnY43depUtW3bVg0aNFCrVq00a9Ys3bp1q5bO1rsQUAAA9dKXX36pX/3qVxo9erTOnj2rt956S4899lilUzupqak6evSo9u7dq4MHD+rtt9/We++9V6HdihUr1L17d2VnZ2vChAkaP368cnJyHO+HhIRo27ZtOnPmjF544QVt2rRJq1atcut5eiumeAAA9dKXX36p0tJSPfbYY7rnnnskSZ06darQ7tq1a9q+fbt27typ/v37S5K2bt2q6OjoCm0HDBigCRMmSPp2tGTVqlV688031a5dO0nSzJkzHW1btmypZ599VmlpaZoyZYrLz8/bEVAAAPXS/fffr/79+6tTp05KTEzUI488ol/84hdq3LixU7sLFy7o1q1b6tmzp2Ob3W53hI7v69y5s+PPNptNkZGRunLlimPb7373O61Zs0affvqprl+/rtLSUoWGhrrh7LwfUzwAgHrJ19dXBw8e1J/+9Cd17NhRa9euVbt27ZSbm1vjffr7+zv9bLPZVF5eLknKyMhQUlKSBgwYoNdff13Z2dmaMWOGSkpK7uo86ioCCgCg3rLZbOrVq5fmzp2r7OxsBQQEaPfu3U5tWrVqJX9/f508edKxraCgQB9//HG1jnXs2DHdc889mjFjhrp37642bdros88+c8l51EVM8QAA6qXMzEylp6frkUceUXh4uDIzM/XVV1+pQ4cO+uCDDxztQkJCNHLkSD333HNq0qSJwsPDNWfOHPn4+Mhms1X5eG3atNGlS5eUlpamHj16aN++fRXCEP6BERQAQL0UGhqqI0eOaMCAAWrbtq1mzpypFStW6NFHH63QduXKlYqPj9egQYOUkJCgXr16qUOHDgoKCqry8X7+859r8uTJSklJUZcuXXTs2DHNmjXLladUp9gsL7xVXmFhoex2uwoKClhcBNRBld2B1JV3CYVrFRUVKTc3V7GxsdX6wvZmN27c0I9+9COtWLFCY8aM8XQ5Rvmh34fqfH8zxQMAwB1kZ2fr3Llz6tmzpwoKCjRv3jxJ0pAhQzxcWd1FQAEAoAqWL1+unJwcBQQEqFu3bnr77bfVrFkzT5dVZxFQAAC4g65duyorK8vTZdQrLJIFAADGIaAAAADjEFAAAC7hhReFwg1c9XtAQAEA3BVfX19J4pbtkCTdvHlTUsXb/lcXi2QBAHfFz89PDRo00FdffSV/f3/5+PD/vvWRZVm6efOmrly5orCwMEdwrSkCCgDgrthsNkVFRSk3N5dny0BhYWGKjIy86/0QUAAAdy0gIEBt2rRhmqee8/f3v+uRk+8QUAAALuHj41NvbnUP9yOgAKhVPGcHQFWwkgkAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADCOywNKWVmZZs2apdjYWAUHB6t169aaP3++LMtytLEsS7Nnz1ZUVJSCg4OVkJCg8+fPu7oUAADgpVweUJYuXar169frxRdf1NmzZ7V06VItW7ZMa9eudbRZtmyZ1qxZow0bNigzM1MNGzZUYmKiioqKXF0OAADwQn6u3uGxY8c0ZMgQDRw4UJLUsmVLvfbaazpx4oSkb0dPVq9erZkzZ2rIkCGSpFdeeUURERHas2ePhg8f7uqSAACAl3H5CMpDDz2k9PR0ffzxx5Kk999/X++8844effRRSVJubq7y8vKUkJDg+IzdbldcXJwyMjIq3WdxcbEKCwudXgAAoO5y+QjKtGnTVFhYqPbt28vX11dlZWVauHChkpKSJEl5eXmSpIiICKfPRUREON673eLFizV37lxXlwoAAAzl8hGU3//+99qxY4d27typ9957T9u3b9fy5cu1ffv2Gu9z+vTpKigocLwuX77swooBAIBpXD6C8txzz2natGmOtSSdOnXSZ599psWLF2vkyJGKjIyUJOXn5ysqKsrxufz8fHXp0qXSfQYGBiowMNDVpQIAAEO5PKDcvHlTPj7OAzO+vr4qLy+XJMXGxioyMlLp6emOQFJYWKjMzEyNHz/e1eUA8AItp+2rdpuLSwa6qxwABnB5QBk8eLAWLlyoFi1a6N5771V2drZWrlyp0aNHS5JsNpsmTZqkBQsWqE2bNoqNjdWsWbMUHR2toUOHurocAADghVweUNauXatZs2ZpwoQJunLliqKjo/Wb3/xGs2fPdrSZMmWKbty4oXHjxunq1avq3bu3Dhw4oKCgIFeXAwAAvJDN+v4tXr1EYWGh7Ha7CgoKFBoa6ulyAFRDVaZzqoIpHsD7VOf7m2fxAAAA4xBQAACAcQgoAADAOAQUAABgHJdfxQMA3+eqRbEA6hdGUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjMNlxgDqrMouceYZPoB3YAQFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHF4mjEAr8STioG6jREUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjOPn6QIAeIeW0/ZV2HZxycA7tqlNnj4+ANdhBAUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcdwSUP7617/q3//939W0aVMFBwerU6dOevfddx3vW5al2bNnKyoqSsHBwUpISND58+fdUQoAAPBCLg8of/vb39SrVy/5+/vrT3/6k86cOaMVK1aocePGjjbLli3TmjVrtGHDBmVmZqphw4ZKTExUUVGRq8sBAABeyM/VO1y6dKliYmK0detWx7bY2FjHny3L0urVqzVz5kwNGTJEkvTKK68oIiJCe/bs0fDhw11dEgAA8DIuH0HZu3evunfvrl/+8pcKDw9X165dtWnTJsf7ubm5ysvLU0JCgmOb3W5XXFycMjIyKt1ncXGxCgsLnV4AAKDucnlAuXDhgtavX682bdrojTfe0Pjx4/Wf//mf2r59uyQpLy9PkhQREeH0uYiICMd7t1u8eLHsdrvjFRMT4+qyAQCAQVweUMrLy/XAAw9o0aJF6tq1q8aNG6exY8dqw4YNNd7n9OnTVVBQ4HhdvnzZhRUDAADTuDygREVFqWPHjk7bOnTooEuXLkmSIiMjJUn5+flObfLz8x3v3S4wMFChoaFOLwAAUHe5PKD06tVLOTk5Tts+/vhj3XPPPZK+XTAbGRmp9PR0x/uFhYXKzMxUfHy8q8sBAABeyOVX8UyePFkPPfSQFi1apMcff1wnTpzQxo0btXHjRkmSzWbTpEmTtGDBArVp00axsbGaNWuWoqOjNXToUFeXAwAAvJDLA0qPHj20e/duTZ8+XfPmzVNsbKxWr16tpKQkR5spU6boxo0bGjdunK5evarevXvrwIEDCgoKcnU5AADAC7k8oEjSoEGDNGjQoH/6vs1m07x58zRv3jx3HB4AAHg5nsUDAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBw/TxcAwHu1nLbP0yUAqKMYQQEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwjp+nCwBgppbT9nm6BAD1GCMoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGcXtAWbJkiWw2myZNmuTYVlRUpOTkZDVt2lSNGjXSsGHDlJ+f7+5SAACAl3DrfVBOnjyp//qv/1Lnzp2dtk+ePFn79u3Trl27ZLfblZKSoscee0xHjx51ZzkAUOH+LheXDPRQJQB+iNtGUK5fv66kpCRt2rRJjRs3dmwvKCjQli1btHLlSvXr10/dunXT1q1bdezYMR0/ftxd5QAAAC/itoCSnJysgQMHKiEhwWl7VlaWbt265bS9ffv2atGihTIyMirdV3FxsQoLC51eAACg7nLLFE9aWpree+89nTx5ssJ7eXl5CggIUFhYmNP2iIgI5eXlVbq/xYsXa+7cue4oFQAAGMjlIyiXL1/WM888ox07digoKMgl+5w+fboKCgocr8uXL7tkvwAAwEwuDyhZWVm6cuWKHnjgAfn5+cnPz09/+ctftGbNGvn5+SkiIkIlJSW6evWq0+fy8/MVGRlZ6T4DAwMVGhrq9AIAAHWXy6d4+vfvr9OnTztte/rpp9W+fXtNnTpVMTEx8vf3V3p6uoYNGyZJysnJ0aVLlxQfH+/qcgAAgBdyeUAJCQnRfffd57StYcOGatq0qWP7mDFjlJqaqiZNmig0NFQTJ05UfHy8HnzwQVeXAwAAvJBb74Pyz6xatUo+Pj4aNmyYiouLlZiYqJdeeskTpQAAAAPZLMuyPF1EdRUWFsput6ugoID1KICb3H5Ds7qKG7UBtac63988iwcAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYByP3KgNgFnqyz1PAHgPRlAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh4ACAACMQ0ABAADGIaAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADCOn6cLAABPajltX4VtF5cM9EAlAL6PERQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGIeAAgAAjENAAQAAxuFpxkA9VNkTfAHAJIygAAAA4xBQAACAcQgoAADAOAQUAABgHAIKAAAwDgEFAAAYh8uMAeAOKrss++KSgR6oBKg/GEEBAADGIaAAAADjuDygLF68WD169FBISIjCw8M1dOhQ5eTkOLUpKipScnKymjZtqkaNGmnYsGHKz893dSkAAMBLuXwNyl/+8hclJyerR48eKi0t1W9/+1s98sgjOnPmjBo2bChJmjx5svbt26ddu3bJbrcrJSVFjz32mI4ePerqcoB6j9vaA/BGLg8oBw4ccPp527ZtCg8PV1ZWln7yk5+ooKBAW7Zs0c6dO9WvXz9J0tatW9WhQwcdP35cDz74oKtLAgAAXsbta1AKCgokSU2aNJEkZWVl6datW0pISHC0ad++vVq0aKGMjIxK91FcXKzCwkKnFwAAqLvceplxeXm5Jk2apF69eum+++6TJOXl5SkgIEBhYWFObSMiIpSXl1fpfhYvXqy5c+e6s1QAcGBaDPA8t46gJCcn68MPP1RaWtpd7Wf69OkqKChwvC5fvuyiCgEAgIncNoKSkpKi119/XUeOHNGPf/xjx/bIyEiVlJTo6tWrTqMo+fn5ioyMrHRfgYGBCgwMdFepAADAMC4fQbEsSykpKdq9e7cOHz6s2NhYp/e7desmf39/paenO7bl5OTo0qVLio+Pd3U5AADAC7l8BCU5OVk7d+7U//zP/ygkJMSxrsRutys4OFh2u11jxoxRamqqmjRpotDQUE2cOFHx8fFcwQMAACS5IaCsX79ekvTwww87bd+6datGjRolSVq1apV8fHw0bNgwFRcXKzExUS+99JKrSwEAAF7K5QHFsqw7tgkKCtK6deu0bt06Vx8eAADUATzNGKhDuDy29tze1zzdGHAtHhYIAACMQ0ABAADGIaAAAADjsAbFVZ63V7KtoPbrAACgDmAEBQAAGIeAAgAAjENAAQAAxiGgAAAA4xBQAACAcQgoAADAOFxmDBiKW6l7l8oeM8DfGVBzjKAAAADjEFAAAIBxCCgAAMA4BBQAAGAcAgoAADAOAQUAABiHy4wrc/uTiV31VGJXPfGYJycDAOo4RlAAAIBxCCgAAMA4BBQAAGAc1qB4mrvWu6DO4VbqAOoTRlAAAIBxCCgAAMA4TPG4U2WXA3sbLo32eheDnqywrWXRTg9UAgBVxwgKAAAwDgEFAAAYh4ACAACMwxoUb+CqtSzeeEmzN9bsIrevHfH0uhHWssChHv+7RO1hBAUAABiHgAIAAIzDFE9dVRcucfZSt9/xtbKpkZoMiVc6xTKNKRYAdRMjKAAAwDgEFAAAYBwCCgAAMA5rUOqKmqw5ceft52tzDUxVjuWu86pkv5WuOfFCdeU8PKnCeqRKnj5dlTa1yvRLiHlsRr3BCAoAADAOAQUAABiHKR54F3dOHdXi0La7hvVrOi1Tm9M5nr4jrSfv0Hv737tUST1VuHS8xr8vtTl9Y9pUkWn14I4YQQEAAMYhoAAAAOMQUAAAgHFYg1IV9em28a66XNmdn3PFsTw8/1xxHQbz4d9X0zUxt68n8cZLpatWsxfeDqAqx68r60Jq89YHdRgjKAAAwDgEFAAAYBwCCgAAMA5rUEzj6Tlh/INpt+v3At5wH5aqqMm9Ujx+Dqav56gj9zC647Fr+1im/T27ECMoAADAOAQUAABgHKZ4UH11ZDrCJTzcF56eUvH4tEYt8cbzrPy2+jXYkaenMEzjDTVWhelTgmIEBQAAGIiAAgAAjENAAQAAxmENCuqnujKPDCfeuFbEXWq1L0z79+Sqx2+4c12GN/SZh9elMIICAACMQ0ABAADGYYoHQL3mldNCpk0PeKO60odVmZrx0nNlBAUAABiHgAIAAIzj0YCybt06tWzZUkFBQYqLi9OJEyc8WQ4AADCEx9ag/O53v1Nqaqo2bNiguLg4rV69WomJicrJyVF4eLinygIAwLt56ZqT23lsBGXlypUaO3asnn76aXXs2FEbNmxQgwYN9PLLL3uqJAAAYAiPjKCUlJQoKytL06dPd2zz8fFRQkKCMjIyKrQvLi5WcXGx4+eCgm9XKBcWFrqnwGLLPfsFAHiPyr5j6tP3gxu+Y7/73rasO/ejRwLK119/rbKyMkVERDhtj4iI0Llz5yq0X7x4sebOnVthe0xMjNtqBADUc0vqxlRJjbnx/K9duya7/Yf37xX3QZk+fbpSU1MdP5eXl+ubb75R06ZNZbPZXHqswsJCxcTE6PLlywoNDXXpvvEP9HPtoJ9rB/1cO+jn2uOuvrYsS9euXVN0dPQd23okoDRr1ky+vr7Kz8932p6fn6/IyMgK7QMDAxUYGOi0LSwszJ0lKjQ0lH8AtYB+rh30c+2gn2sH/Vx73NHXdxo5+Y5HFskGBASoW7duSk9Pd2wrLy9Xenq64uPjPVESAAAwiMemeFJTUzVy5Eh1795dPXv21OrVq3Xjxg09/fTTnioJAAAYwmMB5YknntBXX32l2bNnKy8vT126dNGBAwcqLJytbYGBgZozZ06FKSW4Fv1cO+jn2kE/1w76ufaY0Nc2qyrX+gAAANQinsUDAACMQ0ABAADGIaAAAADjEFAAAIBx6mVAWbdunVq2bKmgoCDFxcXpxIkTP9h+165dat++vYKCgtSpUyft37+/lir1btXp502bNqlPnz5q3LixGjdurISEhDv+veBb1f19/k5aWppsNpuGDh3q3gLriOr289WrV5WcnKyoqCgFBgaqbdu2/LejCqrbz6tXr1a7du0UHBysmJgYTZ48WUVFRbVUrXc6cuSIBg8erOjoaNlsNu3Zs+eOn3nrrbf0wAMPKDAwUP/yL/+ibdu2ub1OWfVMWlqaFRAQYL388svWRx99ZI0dO9YKCwuz8vPzK21/9OhRy9fX11q2bJl15swZa+bMmZa/v791+vTpWq7cu1S3n5988klr3bp1VnZ2tnX27Flr1KhRlt1utz7//PNarty7VLefv5Obm2v96Ec/svr06WMNGTKkdor1YtXt5+LiYqt79+7WgAEDrHfeecfKzc213nrrLevUqVO1XLl3qW4/79ixwwoMDLR27Nhh5ebmWm+88YYVFRVlTZ48uZYr9y779++3ZsyYYf3xj3+0JFm7d+/+wfYXLlywGjRoYKWmplpnzpyx1q5da/n6+loHDhxwa531LqD07NnTSk5OdvxcVlZmRUdHW4sXL660/eOPP24NHDjQaVtcXJz1m9/8xq11ervq9vPtSktLrZCQEGv79u3uKrFOqEk/l5aWWg899JC1efNma+TIkQSUKqhuP69fv95q1aqVVVJSUlsl1gnV7efk5GSrX79+TttSU1OtXr16ubXOuqQqAWXKlCnWvffe67TtiSeesBITE91YmWXVqymekpISZWVlKSEhwbHNx8dHCQkJysjIqPQzGRkZTu0lKTEx8Z+2R836+XY3b97UrVu31KRJE3eV6fVq2s/z5s1TeHi4xowZUxtler2a9PPevXsVHx+v5ORkRURE6L777tOiRYtUVlZWW2V7nZr080MPPaSsrCzHNNCFCxe0f/9+DRgwoFZqri889T3oFU8zdpWvv/5aZWVlFe5WGxERoXPnzlX6mby8vErb5+Xlua1Ob1eTfr7d1KlTFR0dXeEfBf6hJv38zjvvaMuWLTp16lQtVFg31KSfL1y4oMOHDyspKUn79+/XJ598ogkTJujWrVuaM2dObZTtdWrSz08++aS+/vpr9e7dW5ZlqbS0VP/xH/+h3/72t7VRcr3xz74HCwsL9fe//13BwcFuOW69GkGBd1iyZInS0tK0e/duBQUFebqcOuPatWsaMWKENm3apGbNmnm6nDqtvLxc4eHh2rhxo7p166YnnnhCM2bM0IYNGzxdWp3y1ltvadGiRXrppZf03nvv6Y9//KP27dun+fPne7o0uEC9GkFp1qyZfH19lZ+f77Q9Pz9fkZGRlX4mMjKyWu1Rs37+zvLly7VkyRIdOnRInTt3dmeZXq+6/fzpp5/q4sWLGjx4sGNbeXm5JMnPz085OTlq3bq1e4v2QjX5fY6KipK/v798fX0d2zp06KC8vDyVlJQoICDArTV7o5r086xZszRixAj9+te/liR16tRJN27c0Lhx4zRjxgz5+PD/4K7wz74HQ0ND3TZ6ItWzEZSAgAB169ZN6enpjm3l5eVKT09XfHx8pZ+Jj493ai9JBw8e/KftUbN+lqRly5Zp/vz5OnDggLp3714bpXq16vZz+/btdfr0aZ06dcrx+vnPf66f/vSnOnXqlGJiYmqzfK9Rk9/nXr166ZNPPnEEQEn6+OOPFRUVRTj5J2rSzzdv3qwQQr4LhRaPmXMZj30PunUJroHS0tKswMBAa9u2bdaZM2escePGWWFhYVZeXp5lWZY1YsQIa9q0aY72R48etfz8/Kzly5dbZ8+etebMmcNlxlVQ3X5esmSJFRAQYP3hD3+wvvzyS8fr2rVrnjoFr1Ddfr4dV/FUTXX7+dKlS1ZISIiVkpJi5eTkWK+//roVHh5uLViwwFOn4BWq289z5syxQkJCrNdee826cOGC9ec//9lq3bq19fjjj3vqFLzCtWvXrOzsbCs7O9uSZK1cudLKzs62PvvsM8uyLGvatGnWiBEjHO2/u8z4ueees86ePWutW7eOy4zdZe3atVaLFi2sgIAAq2fPntbx48cd7/Xt29caOXKkU/vf//73Vtu2ba2AgADr3nvvtfbt21fLFXun6vTzPffcY0mq8JozZ07tF+5lqvv7/H0ElKqrbj8fO3bMiouLswIDA61WrVpZCxcutEpLS2u5au9TnX6+deuW9fzzz1utW7e2goKCrJiYGGvChAnW3/72t9ov3Iu8+eablf739ru+HTlypNW3b98Kn+nSpYsVEBBgtWrVytq6davb67RZFuNgAADALPVqDQoAAPAOBBQAAGAcAgoAADAOAQUAABiHgAIAAIxDQAEAAMYhoAAAAOMQUAAAgHEIKAAAwDgEFAAAYBwCCgAAMA4BBQAAGOf/AbM+raqHP0BfAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "rng = np.random.default_rng(1)\n", "\n", "s = rng.normal(0.5, 0.05, size=1000)\n", "b = rng.exponential(1, size=1000)\n", "b = b[b < 1]\n", "\n", "ns, xe = np.histogram(s, bins=100, range=(0, 1))\n", "nb, _ = np.histogram(b, bins=xe)\n", "n = ns + nb\n", "\n", "plt.stairs(nb, xe, color=\"C1\", fill=True, label=\"background\")\n", "plt.stairs(n, xe, baseline=nb, color=\"C0\", fill=True, label=\"signal\")\n", "plt.legend();" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Fitting a parametric component and a template\n", "\n", "We now model the peaking component parametrically with a Gaussian. A template fit is an extended binned fit, so we need to provide a scaled cumulative density function like for `iminuit.cost.ExtendedBinnedNLL`. To obtain a background template, we generate more samples from the exponential distribution and make a histogram." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Migrad
FCN = 87.01 (chi2/ndof = 0.9) Nfcn = 194
EDM = 4.93e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4951 0.0020 0 1
2 x0_sigma 0.0484 0.0018 0
3 x1 630 40 0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0_n x0_mu x0_sigma x1
x0_n 1.43e+03 0.000324 (0.004) 0.0185 (0.274) -441 (-0.284)
x0_mu 0.000324 (0.004) 3.87e-06 -5.26e-08 (-0.015) -0.000347 (-0.004)
x0_sigma 0.0185 (0.274) -5.26e-08 (-0.015) 3.21e-06 -0.0185 (-0.251)
x1 -441 (-0.284) -0.000347 (-0.004) -0.0185 (-0.251) 1.69e+03
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 87.01 (chi2/ndof = 0.9) │ Nfcn = 194 │\n", "│ EDM = 4.93e-05 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬──────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼──────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", "│ 1 │ x0_mu │ 0.4951 │ 0.0020 │ │ │ 0 │ 1 │ │\n", "│ 2 │ x0_sigma │ 0.0484 │ 0.0018 │ │ │ 0 │ │ │\n", "│ 3 │ x1 │ 630 │ 40 │ │ │ 0 │ │ │\n", "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌──────────┬─────────────────────────────────────────┐\n", "│ │ x0_n x0_mu x0_sigma x1 │\n", "├──────────┼─────────────────────────────────────────┤\n", "│ x0_n │ 1.43e+03 0.000324 0.0185 -441 │\n", "│ x0_mu │ 0.000324 3.87e-06 -5.26e-08 -0.000347 │\n", "│ x0_sigma │ 0.0185 -5.26e-08 3.21e-06 -0.0185 │\n", "│ x1 │ -441 -0.000347 -0.0185 1.69e+03 │\n", "└──────────┴─────────────────────────────────────────┘" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# signal model: scaled cdf of a normal distribution\n", "def signal(xe, n, mu, sigma):\n", " return n * norm.cdf(xe, mu, sigma)\n", "\n", "# background template: histogram of MC simulation\n", "rng = np.random.default_rng(2)\n", "b2 = rng.exponential(1, size=1000)\n", "b2 = b2[b2 < 1]\n", "template = np.histogram(b2, bins=xe)[0]\n", "\n", "# fit\n", "c = Template(n, xe, (signal, template))\n", "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1=100)\n", "m.limits[\"x0_n\", \"x1\", \"x0_sigma\"] = (0, None)\n", "m.limits[\"x0_mu\"] = (0, 1)\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABQI0lEQVR4nO3de3xT9f0/8NdpoLfQGwV6SxFUQLk4b4O1UJzan3ytKHxLFScw5pigolDwBjqVOcdNxCIT8MIUHeq0i6isYxMUqAPBIW4KDmH0q2164dqWpi2U9Pz+qOc0l5PkJD1JTprX8/HI49FzcnLyzmmSvvv5fN6fjyCKoggiIiIiHYkKdQBEREREzpigEBERke4wQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7PUIdgD/a29tRXV2NhIQECIIQ6nCIiIhIBVEUcebMGWRmZiIqynMbSVgmKNXV1cjOzg51GEREROSHyspKmEwmj8eEZYKSkJAAoOMFJiYmhjgaIiIiUqOxsRHZ2dny33FPwjJBkbp1EhMTmaAQERGFGTXDMzhIloiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEi3bFarRAEAYIgwGq1hjocIgoBJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdnxOUnTt34uabb0ZmZiYEQcCmTZsc7hdFEU888QQyMjIQFxeH/Px8HD582OGYU6dOYcqUKUhMTERycjJmzJiBpqamLr0QIiIi6j58TlCsVit+9KMf4YUXXlC8f/ny5Xj++eexbt067NmzB0ajEePGjUNra6t8zJQpU3DgwAF89NFH2Lx5M3bu3ImZM2f6/yqIiIioWxFEURT9frAg4L333sPEiRMBdLSeZGZm4oEHHsCDDz4IAGhoaEBaWhpee+013H777fjmm28wdOhQfP7557j66qsBAFu2bEFBQQGqqqqQmZnp9XkbGxuRlJSEhoYGJCYm+hs+EemU9BkHgLKyMtxwww0wGAwhjoqIusqXv9+ajkGpqKhAbW0t8vPz5X1JSUkYNWoUdu/eDQDYvXs3kpOT5eQEAPLz8xEVFYU9e/Yonvfs2bNobGx0uBFR92Q2mzF06FB5u6CgAAMGDIDZbA5hVEQUbJomKLW1tQCAtLQ0h/1paWnyfbW1tejXr5/D/T169EDv3r3lY5wtWbIESUlJ8i07O1vLsIlIJ8xmM4qKimCxWBz2WywWFBUVMUkhiiBhUcWzcOFCNDQ0yLfKyspQh0REGrPZbJg7dy6Uep2lfcXFxbDZbMEOjYhCQNMEJT09HQBQV1fnsL+urk6+Lz09HceOHXO4//z58zh16pR8jLOYmBgkJiY63IioeykvL0dVVZXb+0VRRGVlJcrLy4MYFRGFiqYJysCBA5Geno5t27bJ+xobG7Fnzx7k5OQAAHJyclBfX499+/bJx3z88cdob2/HqFGjtAyHiMJITU2NpscRUXjr4esDmpqacOTIEXm7oqICX375JXr37o3+/fujuLgYTz/9NAYNGoSBAwfi8ccfR2Zmplzpc+mll+J//ud/cNddd2HdunVoa2vDfffdh9tvv11VBQ8RdU8ZGRmaHkdE4c3nMuPt27fj2muvddk/ffp0vPbaaxBFEU8++SReeukl1NfXY8yYMVizZg0GDx4sH3vq1Cncd999+PDDDxEVFYVJkybh+eefR69evVTFwDJjou7HZrNhwIABsFgsiuNQBEGAyWRCRUUFS46JwpQvf7+7NA9KqDBBIeqepCoeAA5JiiAIAIDS0lIUFhaGJDYi6rqQzYNCRNQVhYWFKC0tdenuNZlMTE6IIgxbUIhIdziTLFH3xBYUItItq9UKQRAgCAKsVqviMfbJyNixY5mcEEUgJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3WGCQkRERLrDBIWIiIh0hwkKERER6Q4TFCIiItIdJihERESkOz1CHQARkbPq+hb55wPVDYiPP48UYzSykuNCGBURBRMTFCLSFUt9C65/doe8XbR2N6KiYxHX04CtD1zDJIUoQrCLh4h05bT1nOL+ljab2/uIqPthgkJERES6wwSFiIiIdIcJChEREekOB8kSUchZ6lvk8SVHjjWFOBoi0gMmKEQUVDabTf55586dGDYyD2Of2eHhEUQUidjFQ0RBYzabMeSSS+XtgoICXDl0MJoP7QphVESkR2xBIaKgMJvNKCoqgiiKDvtPHqsBNi1G34mPIn5ILgAgKjoWFzyyORRhEpFOsAWFiALOZrNh7ty5LsmJvVPbXoLYbnN7PxFFFiYoRBRw5eXlqKqq8niM7cwJnK06EKSIiEjvmKAQUcDV1NSoOq7urUfRfq5V8b72c60YYUqGIAiwWq1ahkdEOsQEhYgCLiMjI9QhEFGYYYJCRAGXl5cHk8kEQRBCHQoRhQkmKEQUcAaDAatWrQp1GEQURpigEFFQFBYWYuWLG2Doleqw35DQB6njHwxRVESkV5wHhYiCJr/gFmQciEHVqskAgL5FixA38AqI59twMsSxEZG+sAWFiIJKiDLIP8dmD3fYJiKSMEEhooCyWq0QBAGCIKC5WZvy4APVDbDUt2hyLiLSJyYoRKQ7cT09t6oUrd2N/Gd3MEkh6sY4BoWIdKelzYaSyZfj4n69AABHjjVhzhufuRxz2noOWclxoQiRiAKMCQoRhZxzi0lcTwN+PLA3kw+iCMYEhYhCrqWtc5HA0ntykNU3hckJUYRjgkJEujIsMwlGI5MTokjHQbJEFDQ19coLAXobFEtEkYctKEQUNKebzzlsl96Tg/h4I2LQhsHLtX8+q9WKXr06Bto2NTXBaDRq/yREFBBMUIgoZDq6c4ywWrWZH4WIug928RAREZHuMEEhIiIi3WGCQkRERLrDBIWIAqrabjr6o8ebQhgJEYUTJihEFDCW+hZc/+wOefu1Xd+FMBoiCidMUIgoYE5bz3k/iIhIgeYJis1mw+OPP46BAwciLi4OF110EX77299CFEX5GFEU8cQTTyAjIwNxcXHIz8/H4cOHtQ6FiIiIwpTmCcqyZcuwdu1a/P73v8c333yDZcuWYfny5Vi9erV8zPLly/H8889j3bp12LNnD4xGI8aNG4fWVuVZJomIiCiyaJ6g7Nq1CxMmTMBNN92EAQMGoKioCDfccAP27t0LoKP1pKSkBL/+9a8xYcIEXHbZZXj99ddRXV2NTZs2aR0OEYUBm61zscCdO3c6bEvE9s59rZVfO2wTUfejeYKSm5uLbdu24dtvvwUA/Otf/8Knn36KG2+8EQBQUVGB2tpa5Ofny49JSkrCqFGjsHv3bsVznj17Fo2NjQ43IuoezGYzhg4dKm8XFBRgwIABMJvN8r6tZR+gZv298vbx0kWwrJuBrWUfBDVWIgoezae6X7BgARobG3HJJZfAYDDAZrPhd7/7HaZMmQIAqK2tBQCkpaU5PC4tLU2+z9mSJUvwm9/8RutQiSjE3n//fUydOtVhjBoAWCwWFBUVobS0FAAwf9Z0l2NsZ05g3qzp6J9qRGFhYdBiJqLg0LwF5Z133sHGjRvx5ptv4osvvsCGDRuwYsUKbNiwwe9zLly4EA0NDfKtsrJSw4iJKFQefvhhl8QDgLxv7ty5mDt3ruIxHQcCxcXFil1CRBTeNG9Beeihh7BgwQLcfvvtAIARI0bgu+++w5IlSzB9+nSkp6cDAOrq6pCRkSE/rq6uDpdffrniOWNiYhATE6N1qEQUYhaLxe19oiiiqqrKyxlEVFZW4u9//7vcjUxE3YPmLSjNzc2IinI8rcFgQHt7OwBg4MCBSE9Px7Zt2+T7GxsbsWfPHuTk5GgdDhFFAHfdw0QUvjRvQbn55pvxu9/9Dv3798ewYcOwf/9+rFy5Er/85S8BAIIgoLi4GE8//TQGDRqEgQMH4vHHH0dmZiYmTpyodThEpDNR0bH4qqoeJw7vx7XXXqvJOaWWWSLqPjRPUFavXo3HH38c9957L44dO4bMzEzMmjULTzzxhHzMww8/DKvVipkzZ6K+vh5jxozBli1bEBsbq3U4RBRizuXBcQOvAADk5eXBZDLBYrEojjERBAFZWVkA4PYYyejRozWOmohCTRA9fep1qrGxEUlJSWhoaEBiYmKowyEiN0pefgMPzp8HW9NJeZ8hoQ9WPLsSxXdNg9lsRlFREQA4JCCCIACAXMWjdIy9Pd9aEB9vRIoxGlnJcfJ+q9WKXr16AQCamppgNBo1fHVE5Ctf/n5zLR4iCgiz2Yz5s6Y7JCdAR3nw/FnTYTabUVhYiNLSUmRmZjocYzKZUFpaisLCQrfHGHqlyj8Xrd2N8as/Rf6zO2CxWz2ZiMIXExQi0pzNZvNcHozO8uDCwkIcPHhQ3l9WVoaKigqHuU2cj1nzxrvImLHG5ZwtbTYuUEjUTTBBISLNlZeXeywRFsWO8uDy8nIAHZV+krFjxzpsS+z3XTUqF0KU6zFE1H1oPkiWiCKXpb4Fp63nsPfAf1UdX1NTE+CIiChcMUEhIk1Y6lsweunHAIDW74+peoz9ZI2B4LwI4Q033KDYOkNE+sMuHiLShP3YjxjTMBgS+rg9VhAEZGdnIy8vL2DxqFmEkIj0iwkKEWlOiDKg9/Uzle/7oYS4pKQkYK0ZUvmy81T60iKETFKI9I8JChEFRPyQXPSd+KhDOTDgWEIcCJ4qiKR9XGCQSP+YoBBRwMQPyXUoB1YqIdbaF3t2+VRBRET6xASFiALKvhzYXQmxlo4fq1N1HCuIiPSNCQoRdSt9+6WpOi7QFURE1DVMUIioW7lyVC5MJpM8GNdZMCqIiKjrmKAQUbdiMBiwatUqAHBJUoJRQURE2mCCQkTdjppFCIlI3ziTLBF1S4WFhcjPz0dSUhKAjgoiziRLFD7YgkJEqlitVgiCAEEQYLVaFY9pP9eK75aNx3fLxqP9XGuQI+x4/hGmZDlGNYsQEpE+MUEhIiIi3WGCQkRhJyU+GrF2jSGtlV9DbOfMsETdCcegEFHY+XzH33Fm40Py9vHSRUhMTUPPnGkhjIqItMQEhYjCztSpU13W2mk8eQzYvCJEERGR1tjFQ0RhR2khQEBpHxGFKyYoRKQJS32L4v64nqycISLfsYuHiDRx2nrOYbv0nhzExxsRgzYMXh6ioIgobDFBIaKAGJaZBKPR6HbOFCIiT9jFQ0Rhx91CgETUfTBBIaKQMxqNEEURoijCaDR6PObPf/4zAKUkhUkLUXfCBIWIwoq7hQBT09KROv7BEEVFRFpjgkJEYaewsBAHDx6Ut8vKyvD8e/9A/KCfhDAqItISExQiCkvOCwFGcSFAom6FCQoRERHpDhMUIiIi0h3Og0JEASVV3xAR+YItKERERKQ7bEEhIr9Z6lvkKe6PHGsKcTRE1J0wQSEiv1jqWzB66cehDoOIuil28RCRX5wXByQi0hITFCJSxWazyT/v3LnTYVsitns/JpDsn//19/6KypPsdiIKV0xQiMgrs9mMoUOHytsFBQUY95MRaD60S97XfGgXatbf63DMgAEDYDabgxLjgV0fOTz/vdNuxYjhw/HS6296XOOHiPRJEMOw/q+xsRFJSUloaGhAYmJiqMMh6tbMZjOKiopcSoUFQYAoiug78VEAwPFNi10eKy3oV1paisLCQk3jslqt6NWrFwBg48aNmDp1qmI5syAIAXl+IvKdL3+/maAQkVs2mw0DBgxAVVWV22OieqVCAGBrOql4vyAIMJlMqKiocJievqvsE5SsrCxYLJagPj8R+c6Xv9/s4iEit8rLyz0mJwDQ3nTSbXICAKIoorKyEuXl5VqHJ3OXnATr+YlIe0xQiMitmpoaXZ4rHJ+fiHzDBIWI3MrIyNDlucLx+YnIN0xQiMitvLw8mEwmebCrkqheqTD0SnV7vyAIyM7ORl5eXiBCBNAxBsVdjMF4fiLSHhMUInLLYDBg1apVAOCSAEjbqfmz0Dt/luLjpWNKSkoCOkB1+fLlHu8P9PMTkfaYoBCRR4WFhSgtLUVmZqbD/rSMTPSd+Cjih+Qifkgu+k581KUlxWQyBaXEd8KECVj54gaX5zck9MHKFzewxJgoDDFBISKvCgsLcfDgQXm7rKwMW3b/G/FDcuV98UNykTFjjcMxFRUVQUsO8gtucXj+vkWLkHX3euQX3BKU5ycibTFBISJV7LtIxo4dq9hlIkR5PyaQ7J8/Nnu4wzYRhRcmKERERKQ7PUIdABGRP4xGo+PU9vUNoQuGiDTHFhQiIiLSHSYoROQXS32Ly764nhzzQUTaYBcPEfnltPUcAOBXYwZi4hVZAIAYtGGw5ylJiIhUCUgLisViwdSpU5Gamoq4uDiMGDEC//znP+X7RVHEE088gYyMDMTFxSE/Px+HDx8ORChEFEDt51rx+M3DMMKUjIHJPZCZHBeyWFKM0YotOEotPUSkf5onKKdPn8bo0aPRs2dP/PWvf8XBgwfx7LPPIiUlRT5m+fLleP7557Fu3Trs2bMHRqMR48aNQ2trq9bhEFGEyEqOw+Y5Y+TtX+ReAKCzpYeIwovmXTzLli1DdnY2Xn31VXnfwIED5Z9FUURJSQl+/etfY8KECQCA119/HWlpadi0aRNuv/12rUMioghh34JzYd9eAOpCFwwRdYnmLSgffPABrr76atx6663o168frrjiCrz88svy/RUVFaitrUV+fr68LykpCaNGjcLu3bsVz3n27Fk0NjY63IiIiKj70jxBOXr0KNauXYtBgwbhb3/7G+655x7MmTMHGzZsAADU1tYCANLS0hwel5aWJt/nbMmSJUhKSpJv2dnZWodNREREOqJ5gtLe3o4rr7wSixcvxhVXXIGZM2firrvuwrp16/w+58KFC9HQ0CDfKisrNYyYiIiI9EbzBCUjIwNDhw512HfppZfi+++/BwCkp6cDAOrqHPuG6+rq5PucxcTEIDEx0eFGRERE3ZfmCcro0aNx6NAhh33ffvstLrigY0T9wIEDkZ6ejm3btsn3NzY2Ys+ePcjJydE6HCIiIgpDmlfxzJs3D7m5uVi8eDFuu+027N27Fy+99BJeeuklAIAgCCguLsbTTz+NQYMGYeDAgXj88ceRmZmJiRMnah0OERERhSHNE5Qf//jHeO+997Bw4UI89dRTGDhwIEpKSjBlyhT5mIcffhhWqxUzZ85EfX09xowZgy1btiA2NlbrcIiIiCgMBWSq+/Hjx2P8+PFu7xcEAU899RSeeuqpQDw9ERERhTkuFkhERES6wwSFiIiIdIcJChEREekOExQiUsVoNEIURYiiCKPR6PcxgWT//LFx8UF/fiLSDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiv4ntNvnnnTt3wmazeTiaiEg9JihE5Je9n/wVNevvlbcLCgowYMAAmM3mEEbl6sixJnxtaZBvlvqWUIdERCoIoiiKoQ7CV42NjUhKSkJDQwMSExNDHQ5RxDGbzZg0qQiA49eHIAgAgNLSUhQWFoYgsk5/O1CLWW/sc9kf19OArQ9cg6zkuBBERRTZfPn7zRYUIvKJzWbD3Llz4ZycAID0/05xcXHIu3vcJSAtbTactp4LcjRE5CsmKETkk/LyclRVVbm9XxRFVFZWory8PIhREVF3wwSFiHxSU1Oj6XFEREqYoBCRTzIyMjQ9johICRMUIvJJXl4eTCYTAEHxfkEQkJ2djby8vOAGRkTdChMUIvKJwWDAqlWrFO+TqnhKSkpgMBiCGRYRdTNMUIhIFUt9izyXyOBR1+P2hc/B0CvV4RiTyaSLEmMiCn89Qh0AEemfpb4Fo5d+7LT3YmTMWIOqVZMBAGVlZbjhhhvYckJEmmALChF55W7eECGqMxkZO3YskxMi0gwTFCIiItIdJihERESkO0xQiIiISHeYoBAREZHuMEEhom4pxRiNuJ7Kg3Yt9S1BjoaIfMUyYyLqlrKS47D1gWscKpA27bfglU8ruJoxURhggkJE3VZWchyykuPk7a8tDSGMhoh8wS4eIlJktVohCAIEQUBzszXU4RBRhGGCQkRERLrDBIWIvKqpb1Xc724QKhFRV3EMChF5dbq5Y1Dpr8YMxMQrsuT9MWjD4OWhioqIujMmKESk2sX9emF4VpK8bbVybAoRBQa7eIiIiEh3mKAQERGR7jBBIaKI0n6uFT8bdQEEQWAXFZGOMUEhIiIi3WGCQkRERLrDKh4i8pvRaIQoiqEOg4i6IbagEBERke4wQSEiIiLdYYJCREREusMEhYgi1oHqBljqW0IdBhEpYIJCRIqq7f5wHz3eFMJItJNijHbYLlq7G/nP7mCSQqRDTFCIyIWlvgXXPfOxvL3uT3+B2G5z+QMfbrKS41z2tbTZcNp6LgTREJEnTFCIyMW775aiZv298vbx0kWwrJuBb3Z9FMKoiCiSMEEhIgdmsxnzZ02Hremkw37bmROYP2s6zGZziCIjokjCBIWIZDabDXPnzvU4+VpxcTFsNlsQoyKiSMQEhYhk5eXlqKqqcnu/KIqorKxEeXl5EKMiokjEBIWIZDU1NZoeR0TkLyYoRCTLyMjQ9DgiIn8xQSEiWV5eHkwmEwRBULxfEARkZ2cjLy8vyJFpR2zvHD/TWvm1wzYR6QcTFCKSGQwGrFq1yuMxJSUlMBgMQYpIW1vLPlAsn95a9kEIoyIiJUxQiAhWqxWCIEAQBIwbNw4rX9wAQ69Uh2MMCX2w8sUNKCwsDFGUXcPyaaLwwgSFiFzkF9yCjBlr5O2+RYuQdfd65BfcEsKo/MfyaaLwE/AEZenSpRAEAcXFxfK+1tZWzJ49G6mpqejVqxcmTZqEurq6QIdCRD4Qojq7cWKzhztshxuWTxOFn4AmKJ9//jlefPFFXHbZZQ77582bhw8//BDvvvsuduzYgerq6rBtNiYi/WP5NFH4CViC0tTUhClTpuDll19GSkqKvL+hoQHr16/HypUrcd111+Gqq67Cq6++il27duGzzz4LVDhEFMFYPk0UfgKWoMyePRs33XQT8vPzHfbv27cPbW1tDvsvueQS9O/fH7t371Y819mzZ9HY2OhwI6LAsdS3uOyL62kI29WMvZVPA+FfPk3U3QQkQXn77bfxxRdfYMmSJS731dbWIjo6GsnJyQ7709LSUFtbq3i+JUuWICkpSb5lZ2cHImwi+sFp6zmH7dJ7crD1gWuQlRwXooi6xr582l2SEs7l00TdkeYJSmVlJebOnYuNGzciNjZWk3MuXLgQDQ0N8q2yslKT8xKROsMyk8I2OZEUFhaitLQUmZmZDvsNCX0wb+lajoMj0hnNE5R9+/bh2LFjuPLKK9GjRw/06NEDO3bswPPPP48ePXogLS0N586dQ319vcPj6urqkJ6ernjOmJgYJCYmOtyIiHxVWFiIgwcPytuPPLcBWXevx8hrbwxhVESkpIfWJ7z++uvx1VdfOey78847cckll+CRRx5BdnY2evbsiW3btmHSpEkAgEOHDuH7779HTk6O1uEQETmw78a55IqREGr/G8JoiMgdzROUhIQEDB8+3GGf0WhEamqqvH/GjBmYP38+evfujcTERNx///3IycnBT37yE63DISIiojCkeYKixnPPPYeoqChMmjQJZ8+exbhx47BmzRrvDyQiIqKIEJQEZfv27Q7bsbGxeOGFF/DCCy8E4+mJiIgozHAtHiIiItIdJihE5LBI3s6dO9HORfOIKMSYoBBFOLPZjKFDh8rbBQUFuH/iaDQf5tITRBQ6IRkkS0T6YDabUVRUBFEUHfafOlYLbF4RoqiIiNiCQhSxbDYb5s6d65KcdFDaR0QUPGxBIYpQ5eXlqKqq8nrcJ598AqPRGISIQufIsSZ8bWmQt1OM0WE/tT9RuGOCQhShampqND0uHKXEd6zO/MqnFXjl0wp5f1xPQ1gvjkjUHbCLhyhCZWRkaHpcOMpIVl7QtKXN5rKiMxEFFxMUogiVl5cHk8kEQRDcHpOemYW8vLwgRkVE1IEJClGEMhgMWLVqlcdjHlm0xGFxPSKiYGGCQhTBCgsLsfLFDTD0SnXYb0jog74TH0V+wS0hioyIIh0HyRJFuPyCW5BxIAZVqyYDAPoWLULcwCsgRHXPlhOj0SiXVttX7hCRvrAFhYgckpHY7OHdNjkhovDBBIWIiIh0hwkKERER6Q4TFCKKaO3nWvHdsvH4btl4tJ9rDXU4RPQDJihERESkO0xQiIiISHeYoBAREZHuMEEhIiIi3eFEbUQRyFLfIi+Gd+RYU4ijCZ0UYzTieirP+WKpb8HwrKQgR0REEiYoRBHGUt+C0Us/9npcXE8DUozRQYgodLKS47B5zhgMXt6xXXpPDv5+qB6vfFrB1YyJQowJClGE8faHt/SeHMTHG5FijEZWclyQogqdTLvXOCwzCUdPnw9hNEQkYYJCRA6GZSbBaDSGOgwiinAcJEtERES6wwSFiBAVHYuvquohiiJbT4hIF5igEBERke4wQSGiiGaz2eSfd+7ciXa7bSIKHSYoRBSxzGYzhg4dKm8XFBTg/omj0XxoVwijIiKACQpRROIKvh3JSVFRESwWi8P+U8dqcXzTYuz95K8hioyIACYoRBSBbDYb5s6dC1EUFe7t2Pf6c79x6P4houBigkJEEae8vBxVVVUejzlZV4PXzX/F15YGWOpbghQZEUk4URsRRZyamhpVxz38xk4Y9wmI62nA1geuiYiZdYn0gi0oRBRxMjIyVB1n6JUCAGhps3FtHqIgY4JCRBEnLy8PJpMJgiC4PcaQ0AcxpmFBjIqI7DFBIaKIYzAYsGrVKgBwSVKk7d7Xz4QQZQh6bETUgQkKEUWkwsJClJaWIjMz02F/WkYm+k58FPFDcuV97edaMcKUDEEQYLVagx0qUURigkJEEauwsBAHDx6Ut8vKyrBl978dkhMiCg0mKEQU0QyGzm6csWPHOmwTUegwQSGKMEpzesT1NCDFGB2CaIiIlHEeFKJuzmq1olevXgCApqYmnLaeg9jeOUPqwitF3HTjGM7xQUS6whYUogiz95O/omb9vfL2vdNuxejLL4XZbA5hVEREjpigEEWQ999/H88tuAe2ppMO+y0WC4qKipikEJFuMEEhiiAPP/wwpMXw7EmL5hUXF3OBPHC1ZyI9YIJCFEEsFovb+0RRRGVlJcrLy4MYUfg5UN3ABQSJgoCDZInIgdqF9CJV0drdiIqO5QKCRAHGFhQicqB2Ib3uKsUYjbiejnOh2Fc9tVZ+DbHdxgUEiQKMCQpRBMnKygKgvECeIAjIzs5GXl5ecIPSmazkOGyeM0beLupT41D1dLx0ESzrZqD50K5QhEcUMZigEEWQ5cuXK+6XFsgrKSnhTKoAMu26bZ779RyXqifbmRM4vmkxtpZ9EOzQiCIGExSibs6+KiclJQXFi1+AoVeqwzEmkwmlpaUoLCwMdnghZzQaIYoiRFGE0Wh0PUB0rXqSLFu0kFVPRAHCBIWom7FarRAEAYIg4M0338TQoUPl+woKCvB6yW+RdM10eV9ZWRkqKioiMjnpqtpqC3r06MEVjokCQBBFD/8e6FRjYyOSkpLQ0NCAxMTEUIdDpCv2U9sLggDXj7gA+7lQmpqalFsOIpj9NVSD15BIHV/+frMFhagbU/7/I+z+J9G9A9WcF4VIa5onKEuWLMGPf/xjJCQkoF+/fpg4cSIOHTrkcExraytmz56N1NRU9OrVC5MmTUJdXZ3WoRAR+cXXcSVFa3cj/9kdTFKINKR5grJjxw7Mnj0bn332GT766CO0tbXhhhtucOijnTdvHj788EO8++672LFjB6qrq8O+/9u+35/90UThy2w2O4zbUYvzohBpS/OZZLds2eKw/dprr6Ffv37Yt28fxo4di4aGBqxfvx5vvvkmrrvuOgDAq6++iksvvRSfffYZfvKTn2gdEhGRKmazGUVFRW66xjoZEvog+Zpf4OTmFUGKjCjyBHyq+4aGBgBA7969AQD79u1DW1sb8vPz5WMuueQS9O/fH7t372aCQqQlQfBYJkudbDYb5s6d6zU56TPpCcRfeBXE82046fFIIuqKgCYo7e3tKC4uxujRozF8+HAAQG1tLaKjo5GcnOxwbFpaGmpraxXPc/bsWZw9e1bebmxsDFjMRBSZysvLUVVV5fU4QYiCEGWAiLYgREUUuQJaxTN79mx8/fXXePvtt7t0niVLliApKUm+ZWdnaxQhUfe29PmXXCZlMyT0Qer4B0MUkX6pXSTRZq0PbCBEBCCACcp9992HzZs345NPPoHJZJL3p6en49y5c6ivr3c4vq6uDunp6YrnWrhwIRoaGuRbZWVloMImCnv2FSiJyclIv/P38nbfokXIuns94gexK9WZ2kUSDcZkAMoLCBKRdjRPUERRxH333Yf33nsPH3/8MQYOHOhw/1VXXYWePXti27Zt8r5Dhw7h+++/R05OjuI5Y2JikJiY6HAjIlfOFSj3TrsVta/eJ2/HZg+HEMW1dpTk5eXBZDLJ6xK5E5N1KZoP7VJcQJBr8xBpR/MEZfbs2fjjH/+IN998EwkJCaitrUVtbS1aWjrmB0hKSsKMGTMwf/58fPLJJ9i3bx/uvPNO5OTkcIAskR1fS9elChSLxeKw33mhOwCIio7FV1X17tefiUAGgwGrVq0CAI9JSst/P8fxTYsVFxCcN/PnnGqASCOaD5Jdu3YtAOCnP/2pw/5XX30Vv/jFLwAAzz33HKKiojBp0iScPXsW48aNw5o1a7QOJajsm9V37tyJG264QXerwlrqW7zO05BijEaW3UquFB7UVqCwG8KzwsJClJaWYs6cOQ6JXkZmFmqqO7YN/9zo9TxcQJCo67gWjwbMZrPLF5rJZMKqVat0MwGdpb4Fo5d+7PW4uJ4GbH3gGiYpOmC/Hoy3tV62b9+Oa6+91us5+xYtQvxFVwMANt8/BsOzkrQJtpuRvmOAjsUUR48eLW+rUVZWhhtvvDFQ4RGFLa7FE0TumtUtFguKiopgNptDFJkjtTNccjbM8MQKFG3Zt36OHTvW59ZQd1MmEJF6TFC6wFOzurSvuLiYzb0UcL5WoFBguatIJCL1Aj6TbDhyHqvhblyGt4mdRFFEZWUlysvLXcbkkCulMTIcE6NOXl4e0jIyUVdT7fG4mKxLgxRR92I0GiGKImw2GwYMGACLxeJxvM/o0aODGB1R98QExYnSWA134zLUNqurPS6SuRsjwzEx6hgMBiz4zVLMm/lzj8dJJcZxPQ1IMUYHI7RuRar0KSoq8njcf+qaEH9GZIJN1AXs4nGiNP7C3bgMtc3qao+LZO7GvXBMjHr5Bbeg78RHXWeOtdsuvScHm+8fw6SvCwoLC7HyxQ0er3PR2t0Yv/pT5D+7A5b6lmCHSNQtMEFR0H6uFd8tG4/vlo1H+7lWt8d5m9hJEARkZ2cjLy8PgO/zWrhz2HJCPs/ew9X42tLAL8FuyLl0Xc1YpvghuciY0Vmy37doEbLu+YM858nIQZkYnpXE5KSL8gtucbnO9jP2SjPLMsEm8h8TlC7wNLGTtF1SUqLpfCiW+hZc/+wOeZv/qXVPzjPCFhQUYMCAAaqqwuxniuXMsYFjf13bW5scZuyVZpZtPrQrFKERdQtMULpImtgpMzPTYb/JZEJpaanm86CwK6T786V03VLfgq8tDfja0oAjx5qCHSr94OTmFYozyx7ftJjT3xP5iYNkNVBYWIj8/HyHiZ3UziTry2Rc1P15K10XBAGTJk0CAHxbdRz/b/WeYIcYEaSqHXdSjNGI66muZWrZooW4/5d36G5maSK9Y4Kika5O7ETu2bcMKFVFdKcp/NWUrktON7PFLFSykuOwec4YDF7u/djaagunGiDyAxMU0r3iP30p/+xcdtzdpvBnSXr4yPThvcTfK5HvOAZFgf2CatJofC34U5WhJFDxhQPnsTbdbQp/lqR3T/y9EvmOCYqTrWUfoGb9vfK2NBrfn4Fu9mXFb775pmJVxvvvv69JfIGsFtCqPFqr83RnakrXSYc8/F7SM7PkqQaISD0mKHbMZjPmz5quOBp//qzpXVr4b+rUqYpVGVOnTtUkPlYLdA/2pevUPTyyaAnHpBH5gQnKDzxVT0i6svCfpwUFtYpv2aKFXJiwG3A7U2lCHyx5/iWvj4+KjsUFj2zGBY9sRlR0bKDCjHhSpY8oinjOze+r78RHkV9wS4giJApvTFB+4MvCf1qyTzg8jUvxFh/QWS1A4U9pptKsu9djbP44ed++PbsiavyRnrn7fcUPycWRY03yXDWcTJFIPVbx/EAPC/8VFBQgIzMLDy9a4vBfVwzacO2116o6x94D/0WfQVe47D9yrAnt51pR+VzHImfZ80rd/nftXNab3NOXV6EttTF3R84zwrYc3oOJ194t77t32q0wJPRB7+tnIn5ILoCOaqWWNsekhQsDBoe7GXydq9DWTbsKqXa/j3ApgbfH+ZsoGJig/KCrC/95m9hJrZpqC+bN/Dn6TnxU/qMTI7apfvzKfxzDmspPuxSD8xfqB3df3aXzUdc1H/4MJzevcNkvjT+S3i8tbTaUTL4cF/frJR8Tjn8Au6uWNhum/2Gvw75wKYEnCjZ28fzA14X/1OjKeJBT216Sm+8d/iP2UC1gSOiDGNMwv59TSUubDSfOdDZLd6U8Wqsy60hUv/1Vj/fbv18u7tcLw7OS5Bv/8OlbuJTAEwUbE5QfqKme8GXhP+fF3nxlO3MCZ6sO+PSY3tfP1HxhuOZDuzDx2lHyti+L1tnryuJ34c5+vRx/xyI4V2653O/H+4WISM+YoNjxVD2x8sUNHhf+s/8jVPLyG4qLvfnK1nTaZd+8p593Wy0gdQkBHWM3vls2Ht8tG4/2c61+PX/zoV04vmkxjtU6jrtRWrTOE18Wv1PDftChHhfIs5/v5bDlBEYv/RgFz27FCFMyRpiSUfDs1oCsPt1Wf0zT8xEBygl2NQf7UhBwDIqT/IJbkHEgBlWrJgPoGI0fN/AK5Bdc4/Yx9tOti+02WNbN12Q8iqFXinxOyR/3n0T6nb+HZfXPHOJT03LiPAOtp8eJ7Tac2qZc0iotWldcXIwJEyZ4bFVSs/hdcXExPvz0S9Ux24+R0Tt36+VIzfpadr8YjMmanYsIcL+URLSt85+enTt3ql4clcgXTFAUuBuN7459//HZqgOwnTnh9Tn69OmDkydPuk1kpPEkzYd24dTWF+X9x0sXObSgqIkPgPJ5nCpA7Hl7HfZl154WQVNbvv3Fnl0AHMfX+BpzdyPNZ9KR9M7w+r6Kybo0SJGRsxRjNIxGIy54ZHOoQ9GU0tiY5kO7UGX3uSwoKIDJZMKqVas8tjIT+YoJisaUumWUTJ48GWvWrIEgCIpJSu/rZ6Ll8B4c37RY4TlcxyM4l+Pa81YBIrEv41X7OtyVXduXIapx/FgdgPTOmH/oXnIXs32XllIpshblyUqrJMegDYNNfQFoW17p/Fz2XVdClAG9r5+peD3saT3+iNTLSo7D1geucfkdznnjM6/vw/ZzrRhhSgagj5Jd+8/unm8du2XdfS6l7trS0tKgJCnhVuas9F3C6jrvmKBoTOqW8eamm27Cddddhzlz5jiMzZBaCOIGjYJl3Qyv53E3UZf9/tNuumoAdFQF/ZAgLS+6DAs++LYjDpWvQ6tF0Pr2SwMqO3721L0kObXtJcQNGhWwP8rumrbVlHzbVyd1TKYWrdhV5e257MUPyUXfiY/i1NYXHRLUKGNvtFtPKZ6XgisrOa5b/8HRqts30rj7fLO83DsOktVYjGkYDAl9vB43evRoFBYW4uDBg/K+NW+8K88+qbar6KzlG5d9zYc/c1hQsL2l0f0J7FpvLuzb+V+It9fhS9l1VlaW1/LtK0d1dtmoee2BrlpxV/bpPAmaM+dqpXun3Yqq1VNR/WJnsikt8PjGW+/ga0sDPq84pSqm+CG5DrOVJubchii763q8dBGq183A7o/LVJ2PyBe+dPtSJ0/fJSwv94wJipN+CTGK+9VWXEjN8d5I/2HY/6dx1ahcuUVAbReLzVrvsu/k5hVey1K9UfM61JZdL1++vOOcTkmKtF1SUoLaM50fVNWvXeVxweKuWqm99QzaWx2rjWxnTmD5QzNx3X3LfRr0a99i1Lj7HZx3Xjiy6SRmTZ8SEeXbFFxd7fYl8hUTFCf9EmOx46Gfytu/yL0AAPB5xSnV81hIzfEu5cB22weqPZfrqe1ikSo3tFiTxXltF7evQ0XZdeXJzj/IDe3RWLH2VfRNS3c4JiMzS+6zlv6T+NWYgVg+bayqeE98+Azaz7W6dJ+I7TbFffaUSicPW07I5cHNzVbF57Q/j/1kc2oWc1RyfNNi2FqVn8sf0vN3ZWFL0o7z+9DWanUp/3d+TzU2NsrvQ6u1471hX7ou7fOHmvO4dlF2bPva7atVzN5iVDvpoz/xKH1PdHVOo2BwjluPMarBMSgKBmakyl/0fztQi9f/WYdXPq3AK59WyMd46z+MH5KLmAt+5FKuLP0HfNv6/QDcj2mQulg8NakaEvogbuAVLtUu/lJa28Xd6/BUdv3yG2/hnnvvdzlv0tifA39ZKZ8nYfDVGHXddQ6PvbhfL9x61Y14wmSCxWLx+ge/+fBnDrOsHi9dhKjYBACiwz5DQh9sHbQSw++a5vf4EufrbF+90Lt3b6+LObpz1vIN4i/SbjkBtRVWFFhbyz5Azfp58rZzBR6g/J7KysoKWozOzGYz5syZI2/bfyfEDRrl8TtJEASYTCafZtvWIsZAVRGpGRsG6G8siVLceotRLbageOHuF6qm/1BNubK7MQ1quljsK3262qUjkapkmg/tcohF4q2s2Ww2Y9b0KS7x2M6cwKkfkhPpPK025f5ZNbP6SpS6s9x1qcyfNR1ms9mv8SVS9YLzc0nVC++//76qeJUoddNpgU3toWM2mzF/1nTXz4HddvPhzxTfU9XV1UGJ0Zm7LkrpO6Hl8B6330n23bWBHCCr9aSPnqgdH6K3sSRKsegtRrWYoOiYpy6Wh595CR+tegBRe14LyHPbr+2ilr/dHEo8zeqbOv5Bv8/rT9eHt+oFANi4caPfMQVqgjWtKqzIN2o/B/Xb/6C43/5xweqmUxOzVDmn9J1kMpkCXmLsbdJHgF2b3Q0TFBWUpo2X5i7wNF5BjajoWLz3RRW+qqpHdZPrB8+5cqNv0SJk3b0e0352G04f/TfqagLz35btzAm0VOxXdaylvgV7D1ejR48ePnVz2F/D1pZmAEBrS7PcT5z70+sVX3vcRT/27cX8QOr62LVTudnWPiGz73tXU71w/PhxpKSmuq1W8hxXu+pkMCo6Fv0fel+zCivSnrfJCSW2Ju/VW2vXroXNZtNsoU2l81itVlWfXalyzvk7qaysDBUVFS7JifNzfX+yyWX8xt7D1fLnfe/hao9jPNRO+uiuiiiYi5Xaj3eR/j5osfxIpOEYFB3wVsXhrosl0E34aroepP7OQHzgUuKjER8TLW/HZg9Hy+E9XR5vc+L4MQAXOOxzHgtg3/cu2rzPfQIA5weOgXjS966eE39+SnGG3LieBsVuJ08TtwWrqZ3c0/JzuWDBAixZusxhEXN/x1y4G7uxbNky1eeQKnnsv5PGjh3r8l5Tei6l97j990bR2t2KE9lJ4yfUXlel44I1boW0xRaUAJKmKr/gkc1+zWTqjS9N+P50jUhdD55ehz/9mmk/W6zqemQmx2HznDHydlGfGk3G2/Tp289h2934Eqnvve2UulaqeDfN31GxCYiK9TyrrtLYn5Y2G0omX47N949ByeTLHZ/LTfdfMJrayTOtu9Ya6k+j/rRjia+WC3ZOnTpVdSxqKnm8jWWxf4+rIY2fUHtdnY8L5rgV0hYTFBW8lbE6l+cGS15eHkwmk9duhT6TnuiYAG7QT+R9niZPk/jS9aD2OGmNIefH/Gf/XojtNrQ7NcOmJXS2oGxc7Xmqd7Xa2ztfl5pZa8/8a4tLIuBMel3Ozd+Pr34Dpvv/iMxZ61XF5jz25+J+vTA8KwkX93NNcNQ2tVNwqf1cdoUoihABzL5/Dv71/SmPpaTexm6oHTNm/9l1R+1YFvvPn0RpOgB7eXl5SMvIdHu/UtdmqMatKJVq+/JaqQO7eLxQKhV0LmNVKs/1ly9ryEjVLkVFRR7PGdf/MghRBojo7KpYvnw5pk6d6nYtIMB91wPQuVbMkWNNPpU5975+JoQog8tjls2bjqjYBMx5sbOp2Lnk8litNk3n902fLL+uqLheXmetbVfRYiO9LsCx+fu1o3GIijaonpLfduYEvn9mgsvvPsUYrdjlY98FptTUTsFn/7n09PnqMlFEbbUF4xa8iNj+lyFGbMO3y/8XgOP6NGrHxHhj/x53R81zSe/x1PEPukwR4O37Zvr8RVj+kHIlkTQGpbW1VfVr16okv7nZCkFIBtAxYP7hhx+W77t32q1epz4A1K/XE27rEHUFExQPpFJB5y+Y9tYzLscqLWKnFamLRYlU7fLg/HkOXRSGXqkuXRZR0bH4qqoew7OSAACxsbEuawE5c/e6pHEz7hYPc2b/xePuMe2tZ2B1GspSXd0xiG7u3LkoKSnx+jwAkJraMY/NqVPuByFKryvh6gmqzglA/kJ1uM5BWl1ZaTE64IfFC5cH9KnJD4WFhSgtLfX6+dKCNC7EXZl8V8fEKL3Hnb9L/HkuTwuYuvu+ATIV16RS+r7zJR4txw1NnTpV9d+M+bOmo3+qEaOuu5Hr9ShgF48b/pbM+lOeay/W7h8Utc2At95ahAvv7myN6Fu0yKHpXxLX04AUY+d/3M5rAXmi9LrUdI8Adl1MQ3JVP0Z+Dh/LeJcuXYq6ujpUVFR4PxiA9cAnqmOJiu2F9Dt/L29LVUWBSE6UfvdZyXEYnpXkcMuM0C+ucODL56srpHEh7mY57sqYGPvPrhpajb/x9D2qVNlo/7n057V3NW77Lh1f/2YUFxfjRKNy95zS/CXBrEYKNSYobvjbLNqVReyaD+3C0XWz5O3jpYtg3XC310FlWclx+Etx58yumxffhT/P7hxcWnpPDjbfP0YxE1fbJaD0utQuaCgIUXLTsNrH2LMv4/XmnnvugcFgUP262lsaERWXqOrY46WLUPvqffK2t0nrJHE9fe92kRYU3Fr2gc+PJf2wfx8GYkyKNC6k+dAuhwVCCwoKkJXdHyUvv4GUCy9DWkamX89v/9n1RJpavSvPZc/b96h9TO2tTQ6fy4KCAmT3v0DVa9eiJL/50C5MvHaUX4+Vupg+/Ns2t8ccOdZZnv3yG285LEZq/1qlY6TuME/nCZfp79nF40ZXmvz8WcTOXbfHqeO1EJ2aPJ1bQgA4/Cc9LDMJRqNR875v6XXZj5NR9Ti7cuWuLPA36baf4ZW1v/d4jD9jMIzDrsWZf6orD1ZqRlYaN2TfJdfSZsPz036Cix+ox9ayDzBv1nSHVaTdPpddE7C7ga+B+D1T4DiPSenqGBX72aSd1dVUY97Mn/t9bqDjsxvX0wCr1eryHrfUt2B4VpLL1Orto34BUUW3rzd1bz2K7HmlAOBxXJ5SV1FNtQXzZv4cfSc+CniJx9eSfPvPe+r4BxWf31erPvwcxqHKy4d46063f62eWrqcp7NwN2ZJT9iC4kZXmvyWTxuLzfePcbg5l4na8zZTqSAIMOzdgPfvdd8SEgxKTcmqHmc3U6raBceUTLltEl56/U30S3f83WRkdm3tEnflwd74UokgVeMU3zUNzynMkOsJZ8fsHv74xz8iM9OxCiUtI9Oh/D91/IOqytQNCX3Qd+KjiBs0SnWXaWLObTD06u1TzAZjMlrabFj6v53VOz9NOgGxvbPrwbkLwl0JvD9aK79G+/k2h21fvn88zX6rZtFTJfbPf9qH7mpPvH0vquka93V4QfPZzt+bXruKmKC44U+poNRc+PPCG13GCiiViUrUzFRaW23B6aP/xvCspIAlJ55eq7umZDVisi7t/PmHRRB9jUtqhr1r2s9w+NB/5PvKysrwn2+89POreF3O/dpqnLV843I9pK4ZT91y+QW3qH4ub7NjUviYMGGCw5iUsrIybNn9b4fy//hBP3EZX+Fcpm4/9smXLtPG3e84NNwl9/aeQMRkXYrmQ7sw77bORT3fWHQ3LOtmYO8nf3X7OH8+T0qOly6C5fdTHLfXzUDz4c9UPd7d7LfSNcwvuMWneJw/7+0tjT49Xoma8m01v2dfhhcodQkOGDBAd3PCMEFxw5cF6+z5M4On2m4PT91OUlO/KIo+NdVJj/vzn/8MwH2S4u/ChH0nPgpDbGc8ahZBtKc0M6r99R07diwSExNdXrv96/KUYrorD1aj9f/2e5zgzV2SkmJ0nCFXDS78F56cP5fO712DweAyEaLSzNHuZpP2tcu03dpZ2XbnA4s8Htt34qM4+92/cHzTYpcSf9uZE3huwT0oefkNt2MefP08uSW2uzz3yc0rVFfgKc1+K11D53EZ9jfn1+VuQseuUlO+rfb3rOY4bwuf6ilJYYLigbsF69w1ufrTXAio7/YI5OJvUlmkcxN0ZpYJz730us8LE0pN0Ep9ou6agKNiE5CU7HgtujozqrdFB49vWuz32hjeKoDcNblmOc2QqwYX/iMlXekyffdYP7ddH+q6j0Q8+MB8zH1rH4DgrzWjtgLvxIfPuI2n+E9fYvzqT+VbwbNbMcKUjBGmZMx5o7OVxtfqQ2dRsQmIT3AsyZau+/FNix3WeLO/htL2iQ+fUfU8XekqkhLpSZMmobGx6y1DWuAgWS/yC25BxoEYVK2aDKCjaTBu4BVoP9fqsi+/QHmQkzdSt4e7JjxBEGAymQK++FthYSHy8/ORlNTxQSorK8MNN9wAg8GA7du3q16YULoenv4riB+Si5gLfuRyDd/+5RXIvbS/y/Pb83VgqKffoaS18mvEDbwCFzyyGWK7DZZ1Mzw2qUbFJXpt3pWaXGP7X+Zy36CsPhBFETabDQMGDIDFYlF8TcH63VNwuLx36xtcjlGa98jdXEjevjs8aTv2X7efQyHKgNbv/626WyG2/2UuY7F8+Tz5Q6rAU9PNYh+PJ0qvQYgy+FV9KOkz6QnEX3iVw3fbmjfexe/2nodl9c8cnsv5+aMzhsjb3l6rVl1FQMcilQ8++GDIJ35kC4oXzs3x3ppcPZ3HXbmpp26PYC/+ptQEDfjWxRCbPVxVF4bSNXT3/F3l/Fy2is/djh1x+H246fIyDrtW1fN6a3K170p07l7jwn/kja9dpvY8dX3Y36/mPJ7GYnUlRm/Ufg7VjA1Teg3VPzymK9WH0kze9p/hxvp6h9Lo46WLULV6KqpfnOEYs934G2+JmH1Xkbu/NWpfx4IFC3QxJoUJihfOzfGl9+R4rMjxdJ6tD1zjtrJH74u/+drFYL/QndLNn2uolebDn8FS+rTHsSPS76NfWrrDMXLX1SB18x6oaYJ3172ml9896Zu/VTPe3ptqu4/aTlu8jsVyF2NXK318qcDzNDbM3biM82dO4Pj7SzB5sG+dDd7iWTBnpstztbeeQXur03gep/E3is+l0J3ubqFRX7oE9TAmhV08KkjN8ZKvLQ0ep593Jys5zqECx3l9lfghuUi++EocXtFRY++uiyOQ3HWfSFVN7roigI4PStbd6yFEGRDX04AfD+ztseJI6RrGxwdmXg/puaTmZk+k0sT4Ibl4a8ndcrPsG+9swtNf9kCrraMp2FvTupomV4mn7jWKLEprLjmztVodu2UGjUKGXVeNL90Bzp9D6fnVdB9F9UpF05d/8xir9HlKHZ7ntjvJn24g6TUIUQaH8yb3TkX9KfcDWaV4hCiDw3V0RwDw/p/e8Pr9J3HXxR0fb8T58+cxYMAATdZGAly7iuyfW5raAHB8v6jtFpOmuCguLsaECRNC8l3EBCWElNZXiUEbBv8w74+eFn/zujChIGDFsyvlcThKi1zpgS/lerH9L3O4/v9bkI9r/1+U/PvaOmhlx1pNgOLEa2pG59sLVPcW6Ze7RSA3/HIkUn+YjPGk9RzufmOffIzzQpvSwnPJ1/yi87zXz/Q4gZin96b03/fF/XrJ73FAeQr3hB/9Dxr+4XkZCunzJPS/DM/cdgV+9kNx5ObFdyE+vrPCz9vnyYHd982m/Ra89PE38l2ekhP7eNpbmlQtciqKIqqqqvCb3/wGixYt8jq5nqcuf60WbpR0dBUtkLftF1wEOlr+nRe89aU0WqvFFP3FBCXEnFtVAN/XcggWtwsTJvTBimdXoviHVTn1zNdyPecWHSMg/76G3zUN/VONLgvC9Ujsi5Tr7vI4868SzgobeZT+SVFK7qVjtpZ9gHnLl7j8AZdKb+2b+qN6RLsuqqdycUvpv29373HpPKKtzcNZ7OL74fM0/IJ+bt/j7p4LQpRDV4fz943Uov3Wnu8g/vcfuOOOO7zG03x4j+rZoyWDBg1SXADSlwVDtZwuoLi4GAvnzHS5nlJX1taCS/BtqlFxwVtfhWqaAyYo5JNbby3Cmm/j5W6ovkWL0Hvw1bj11uu8PDJ07P9LVdsHq/Y4pa6ZYSPz0Nja+R+xXluTSB+U/klROiY9IRo3PvWox9YF++4LTxU6EqXWG6VFRe3f42veeBdL/x0tV/qo4e/n6Y13NuE3n53FkZWdr8HT943asXK+LBJqf+6f/vSnbq+F2nNoZePGjR4Tj98smIcoiJr80xOqaQ6YoJBPpIUJpW6ozYvvQlbfFF3/Abb/L9Vmy8G47b/Hsdoaj2Np1I4dAVy7ZvS4pgWFPzXdA86l7d6qDe27cyRKCbX9e/yqUbkQvt4PQF2Zc1c+T/9bkI9RuS0YvLJj29P3zZFjTbj0ssuQ2i8DJ4+5/49f7RgMiXOpv7troYaasXxq4unTpw+OHz/u8ThvXV1OZwWgHE96ZlbIpjlgFY8flEqG1TbjdwfSoGFRFDFyUKZfyUmwr2FWchyGZyXhR/17Y83vVwPwPGuuNNDXl64ZX2fxJfKF2mZ2+25M51lqlcjdOT/clD7P9u9x+3EjakqIu/p58vZ9I53zlU8rMGHNbkTl3un5/CpLk+3Zl/q7uxZq+DtDuT0RwJQpU7wep1bHjLzuk6VHFi0J2Xg4tqD4QW2/MbkXymsolfU69yVnZpnw0JOL5fU5+DslPVHbzL582lj8ONd1puIjx5pcVrTVglRCbPvHH3DqeK28Pz0zC48sWhLwz5PzOb3Fk5Scgl/epm78iSGhD5Y/swK9e/fGW2+9hYyMDOTl5XXpD7b0/TP7vvtRazf5ZUrvVAgQcepU53IEUVFRaG93HH8z57GnccuoISgpKfE7Bnvxg0Yh1jTM7XglX9cr0lJIE5QXXngBzzzzDGpra/GjH/0Iq1evxsiRI0MZkmpq+o3Js1Bew8LCQkyYMAHl5eWoqanR5IuHKJC8dQ9IXRE/L7xR8X2sVDHkT6ul0nlSh+fhb68sxNGv/qmLz1P8kFy8v+oBnD76b5d4bDYbHsvMQk21xe3jhdgE9J3wCHq0NWPFU485HGsymbBq1SoUFha6rcJypjSuR+n7B4DDvtzcXOzatQs1NTX4pj4KGyriMPLay5F3VRZMJhOqqizw1PrhjX2pdtygUR1Vjk2nYeiVIu8PJUEMUdnAn/70J/z85z/HunXrMGrUKJSUlODdd9/FoUOH0K9fP4+PbWxsRFJSEhoaGpCYmBikiImIQstsNsul/vZf3VJ3pbeJ/Sz1LZq0Wmp1Hi18bWnA+NWfuuzffP8YeR4QZ9J1dC5plq7jyhc3AIBiBYzztXa+Fkq0uD5v7/0eC8xfYWnhCNw+sj/MZjMmTSqCPwmKVCrtbr00e56uoz98+fsdsjEoK1euxF133YU777wTQ4cOxbp16xAfH48//OEPoQqJiEjXpO6BrKwsh/1qZx2WxmJ5Gm+ihlbnCRXpOprcXMf7f3kHnn3qUcWWKmlfcXExbDaby7VQugXi+hQWFmLe0rUwJPTx+bFpGZkuyYnS9PihHlsZki6ec+fOYd++fVi4cKG8LyoqCvn5+di9e7fL8WfPnsXZs2fl7YaGjgW29LLiIhFRsOTn5+Pf//43du3ahdraWqSnpyM3NxcGgyEivxObzjSi/Wyz4v7GRuWB8IDn67hlyxaPFVPSBGZbtmwJWoVLc9MZtJ9txof/PILmpjMAgLqEQci4czXuuKAFFyecx8KFC3HypPvqnZSUFLz22mtIvegy/OyVzx2um/UssLRwBC7s2znwNzk+GglRbWhsVDffjRrSe1RV540YAhaLRQQg7tq1y2H/Qw89JI4cOdLl+CeffFJERzsWb7zxxhtvvPEW5rfKykqvuUJYVPEsXLgQ8+fPl7fb29tx6tQppKamui0V9VdjYyOys7NRWVnJ8S0BxOscHLzOwcHrHBy8zsETqGstiiLOnDnjsjiqkpAkKH369IHBYEBdXZ3D/rq6OqSnp7scHxMTg5iYGId9ycnJgQwRiYmJ/AAEAa9zcPA6Bwevc3DwOgdPIK61NBOvNyEZJBsdHY2rrroK27Ztk/e1t7dj27ZtyMnJCUVIREREpCMh6+KZP38+pk+fjquvvhojR45ESUkJrFYr7rzT8yyARERE1P2FLEGZPHkyjh8/jieeeAK1tbW4/PLLsWXLFqSlpYUqJAAd3UlPPvmkS5cSaYvXOTh4nYOD1zk4eJ2DRw/XOmQTtRERERG5w8UCiYiISHeYoBAREZHuMEEhIiIi3WGCQkRERLoTkQnKCy+8gAEDBiA2NhajRo3C3r17PR7/7rvv4pJLLkFsbCxGjBiBsrKyIEUa3ny5zi+//DLy8vKQkpKClJQU5Ofne/29UAdf38+St99+G4IgYOLEiYENsJvw9TrX19dj9uzZyMjIQExMDAYPHszvDhV8vc4lJSUYMmQI4uLikJ2djXnz5qG1tTVI0YannTt34uabb0ZmZiYEQcCmTZu8Pmb79u248sorERMTg4svvhivvfZawOMMyVo8ofT222+L0dHR4h/+8AfxwIED4l133SUmJyeLdXV1isf/4x//EA0Gg7h8+XLx4MGD4q9//WuxZ8+e4ldffRXkyMOLr9f5jjvuEF944QVx//794jfffCP+4he/EJOSksSqqqogRx5efL3OkoqKCjErK0vMy8sTJ0yYEJxgw5iv1/ns2bPi1VdfLRYUFIiffvqpWFFRIW7fvl388ssvgxx5ePH1Om/cuFGMiYkRN27cKFZUVIh/+9vfxIyMDHHevHlBjjy8lJWViY899phoNptFAOJ7773n8fijR4+K8fHx4vz588WDBw+Kq1evFg0Gg7hly5aAxhlxCcrIkSPF2bNny9s2m03MzMwUlyxZonj8bbfdJt50000O+0aNGiXOmjUroHGGO1+vs7Pz58+LCQkJ4oYNGwIVYrfgz3U+f/68mJubK77yyivi9OnTmaCo4Ot1Xrt2rXjhhReK586dC1aI3YKv13n27Nnidddd57Bv/vz54ujRowMaZ3eiJkF5+OGHxWHDhjnsmzx5sjhu3LgARiaKEdXFc+7cOezbtw/5+fnyvqioKOTn52P37t2Kj9m9e7fD8QAwbtw4t8eTf9fZWXNzM9ra2tC7d+9AhRn2/L3OTz31FPr164cZM2YEI8yw5891/uCDD5CTk4PZs2cjLS0Nw4cPx+LFi2Gz2YIVdtjx5zrn5uZi3759cjfQ0aNHUVZWhoKCgqDEHClC9XcwLFYz1sqJEydgs9lcZqtNS0vDf/7zH8XH1NbWKh5fW1sbsDjDnT/X2dkjjzyCzMxMlw8FdfLnOn/66adYv349vvzyyyBE2D34c52PHj2Kjz/+GFOmTEFZWRmOHDmCe++9F21tbXjyySeDEXbY8ec633HHHThx4gTGjBkDURRx/vx53H333Xj00UeDEXLEcPd3sLGxES0tLYiLiwvI80ZUCwqFh6VLl+Ltt9/Ge++9h9jY2FCH022cOXMG06ZNw8svv4w+ffqEOpxurb29Hf369cNLL72Eq666CpMnT8Zjjz2GdevWhTq0bmX79u1YvHgx1qxZgy+++AJmsxl/+ctf8Nvf/jbUoZEGIqoFpU+fPjAYDKirq3PYX1dXh/T0dMXHpKen+3Q8+XedJStWrMDSpUuxdetWXHbZZYEMM+z5ep3/+9//4v/+7/9w8803y/va29sBAD169MChQ4dw0UUXBTboMOTP+zkjIwM9e/aEwWCQ91166aWora3FuXPnEB0dHdCYw5E/1/nxxx/HtGnT8Ktf/QoAMGLECFitVsycOROPPfYYoqL4P7gW3P0dTExMDFjrCRBhLSjR0dG46qqrsG3bNnlfe3s7tm3bhpycHMXH5OTkOBwPAB999JHb48m/6wwAy5cvx29/+1ts2bIFV199dTBCDWu+XudLLrkEX331Fb788kv5dsstt+Daa6/Fl19+iezs7GCGHzb8eT+PHj0aR44ckRNAAPj222+RkZHB5MQNf65zc3OzSxIiJYUil5nTTMj+DgZ0CK4Ovf3222JMTIz42muviQcPHhRnzpwpJicni7W1taIoiuK0adPEBQsWyMf/4x//EHv06CGuWLFC/Oabb8Qnn3ySZcYq+Hqdly5dKkZHR4ulpaViTU2NfDtz5kyoXkJY8PU6O2MVjzq+Xufvv/9eTEhIEO+77z7x0KFD4ubNm8V+/fqJTz/9dKheQljw9To/+eSTYkJCgvjWW2+JR48eFf/+97+LF110kXjbbbeF6iWEhTNnzoj79+8X9+/fLwIQV65cKe7fv1/87rvvRFEUxQULFojTpk2Tj5fKjB966CHxm2++EV944QWWGQfK6tWrxf79+4vR0dHiyJEjxc8++0y+75prrhGnT5/ucPw777wjDh48WIyOjhaHDRsm/uUvfwlyxOHJl+t8wQUXiABcbk8++WTwAw8zvr6f7TFBUc/X67xr1y5x1KhRYkxMjHjhhReKv/vd78Tz588HOerw48t1bmtrExctWiRedNFFYmxsrJidnS3ee++94unTp4MfeBj55JNPFL9vpWs7ffp08ZprrnF5zOWXXy5GR0eLF154ofjqq68GPE5BFNkORkRERPoSUWNQiIiIKDwwQSEiIiLdYYJCREREusMEhYiIiHSHCQoRERHpDhMUIiIi0h0mKERERKQ7TFCIiIhId5igEBERke4wQSEiIiLdYYJCREREusMEhYiIiHTn/wN1ymumgQVVkgAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m.visualize()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The fit succeeded and the statistical uncertainty in the template is propagated into the result. You can play with this demo and see what happens if you increase the statistic of the template.\n", "\n", "Note: the parameters of a parametric components are prefixed with `xn_` where `n` is the index of the component. This is to avoid name clashes between the parameter names of individual components and for clarity which parameter belongs to which component." ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## Extreme case: Fitting two parametric components\n", "\n", "Although this is not recommended, you can also use the `Template` class with two parametric components and no template. If you are in that situation, however, it is simpler and more efficient to use `iminuit.cost.ExtendedBinnedNLL`. The following snipped is therefore just a proof that `iminuit.cost.Template` handles this limiting case as well." ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Migrad
FCN = 92.87 (chi2/ndof = 1.0) Nfcn = 190
EDM = 1.67e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4964 0.0019 0 1
2 x0_sigma 0.0487 0.0016 0
3 x1_n 629 29 0
4 x1_mu 0.97 0.14 0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0_n x0_mu x0_sigma x1_n x1_mu
x0_n 1.3e+03 0.000424 (0.006) 0.0123 (0.209) -264 (-0.252) -0.369 (-0.076)
x0_mu 0.000424 (0.006) 3.44e-06 -9.8e-08 (-0.032) 0.000293 (0.005) -1.64e-05 (-0.065)
x0_sigma 0.0123 (0.209) -9.8e-08 (-0.032) 2.64e-06 -0.0129 (-0.271) -1.64e-05 (-0.075)
x1_n -264 (-0.252) 0.000293 (0.005) -0.0129 (-0.271) 851 0.337 (0.085)
x1_mu -0.369 (-0.076) -1.64e-05 (-0.065) -1.64e-05 (-0.075) 0.337 (0.085) 0.0183
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 92.87 (chi2/ndof = 1.0) │ Nfcn = 190 │\n", "│ EDM = 1.67e-05 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬──────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼──────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", "│ 1 │ x0_mu │ 0.4964 │ 0.0019 │ │ │ 0 │ 1 │ │\n", "│ 2 │ x0_sigma │ 0.0487 │ 0.0016 │ │ │ 0 │ │ │\n", "│ 3 │ x1_n │ 629 │ 29 │ │ │ 0 │ │ │\n", "│ 4 │ x1_mu │ 0.97 │ 0.14 │ │ │ 0 │ │ │\n", "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌──────────┬───────────────────────────────────────────────────┐\n", "│ │ x0_n x0_mu x0_sigma x1_n x1_mu │\n", "├──────────┼───────────────────────────────────────────────────┤\n", "│ x0_n │ 1.3e+03 0.000424 0.0123 -264 -0.369 │\n", "│ x0_mu │ 0.000424 3.44e-06 -9.8e-08 0.000293 -1.64e-05 │\n", "│ x0_sigma │ 0.0123 -9.8e-08 2.64e-06 -0.0129 -1.64e-05 │\n", "│ x1_n │ -264 0.000293 -0.0129 851 0.337 │\n", "│ x1_mu │ -0.369 -1.64e-05 -1.64e-05 0.337 0.0183 │\n", "└──────────┴───────────────────────────────────────────────────┘" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# signal model: scaled cdf of a normal distribution\n", "def signal(xe, n, mu, sigma):\n", " return n * norm.cdf(xe, mu, sigma)\n", "\n", "# background model: scaled cdf a an exponential distribution\n", "def background(xe, n, mu):\n", " return n * truncexpon.cdf(xe, xe[0], xe[-1], 0, mu)\n", "\n", "# fit\n", "c = Template(n, xe, (signal, background))\n", "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1_n=100, x1_mu=1)\n", "m.limits[\"x0_n\", \"x1_n\", \"x0_sigma\", \"x1_mu\"] = (0, None)\n", "m.limits[\"x0_mu\"] = (0, 1)\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAABAjklEQVR4nO3de3RU5aH+8WdnkAQTEgyQGxMBOVqE2uMFoRFDpc2SU6onnBCv4KEu6xV7AE9FqIottYIsj1wsytF6PYBUQpBKo7bFIKlSUNTz8wBi1aghkKClJCQKSLJ/f+CMM8lMZs/Mnpk9yfez1l4re887e97Zc3vyvvvdr2GapikAAAAHSUl0BQAAADoioAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMchoAAAAMfplegKRKK9vV379u1T3759ZRhGoqsDAAAsME1Thw8fVkFBgVJSum4jScqAsm/fPhUWFia6GgAAIAJ1dXVyu91dlknKgNK3b19JJ55gZmZmgmsDAACsaG5uVmFhofd3vCtJGVA83TqZmZkEFAAAkoyV0zM4SRYAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQWA47S2tsowDBmGodbW1kRXB0ACEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjEFAAAIDjhB1QtmzZoksvvVQFBQUyDEPPP/+83+2maWrevHnKz89Xnz59VFJSor/97W9+ZQ4ePKgpU6YoMzNT/fr103XXXaeWlpaonggAAOg+wg4ora2t+ud//mctX7484O2LFi3SsmXLtGLFCm3btk3p6emaMGGCjhw54i0zZcoU7dy5U3/605+0ceNGbdmyRTfccEPkzwIAAHQrhmmaZsR3NgytX79ekyZNknSi9aSgoED/+Z//qZ/97GeSpKamJuXm5uqpp57SlVdeqd27d2vEiBF64403NGrUKEnSSy+9pIkTJ2rv3r0qKCgI+bjNzc3KyspSU1OTMjMzI60+AIfyfMYlqaqqShdffLFcLleCawUgWuH8ftt6Dkptba0aGhpUUlLi3ZaVlaUxY8Zo69atkqStW7eqX79+3nAiSSUlJUpJSdG2bdsC7vfo0aNqbm72WwB0T5WVlRoxYoR3feLEiRoyZIgqKysTWCsA8WZrQGloaJAk5ebm+m3Pzc313tbQ0KCcnBy/23v16qXs7GxvmY4WLFigrKws71JYWGhntQE4RGVlpcrLy1VfX++3vb6+XuXl5YQUoAdJilE8c+fOVVNTk3epq6tLdJUA2KytrU0zZsxQoF5nz7aZM2eqra0t3lUDkAC2BpS8vDxJUmNjo9/2xsZG7215eXk6cOCA3+3Hjx/XwYMHvWU6Sk1NVWZmpt8CoHupqanR3r17g95umqbq6upUU1MTx1oBSBRbA8rQoUOVl5enTZs2ebc1Nzdr27ZtKioqkiQVFRXp0KFD2rFjh7fMK6+8ovb2do0ZM8bO6gBIIvv377e1HIDk1ivcO7S0tOiDDz7wrtfW1uqdd95Rdna2Tj31VM2cOVP33nuvTj/9dA0dOlR33323CgoKvCN9zjzzTP3Lv/yLrr/+eq1YsUJfffWVbr31Vl155ZWWRvAA6J7y8/NtLQcguYU9zHjz5s0aP358p+3Tpk3TU089JdM0dc899+jRRx/VoUOHdOGFF+rhhx/WGWec4S178OBB3XrrrXrhhReUkpKiyZMna9myZcrIyLBUB4YZA91PW1ubhgwZovr6+oDnoRiGIbfbrdraWoYcA0kqnN/vqK6DkigEFKB78ozikeQXUgzDkCRVVFSorKwsIXUDEL2EXQcFAKJRVlamioqKTt29brebcAL0MAQUAHHV2toqwzBkGIZaW1s73V5WVqZdu3Z516uqqlRbW0s4AXoYAgoAx/E9x2TcuHGccwL0QAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAHHV1tbm/XvLli1+6wDgQUABEDeVlZUaMWKEd33ixIkaMmSIKisrE1grAE5EQAEQF5WVlSovL1d9fb3f9vr6epWXl/uFlPT0dJmmKdM0lZ6eHu+qAnAAAgqAmGtra9OMGTNkmman2zzbZs6cSXcPAC8CCoCYq6mp0d69e4Pebpqm6urqVFNTE8daAXAyAgqAmNu/f7+lcuPHj1dra2vA21pbW2UYhgzDCFoGQPdBQAEQc/n5+YmuAoAkQ0ABEHPFxcVyu90yDCNwAZ/ta9+s0xsfH4xTzQA4FQEFQMy5XC4tXbpUkjqHFMOQfE6enbdhpy5bsZWQAvRwBBQAcVFWVqaKigoVFBT4be+fk6f+l/zMuz7ujAGSpA8PtMS1fgCchYACIG7Kysq0a9cu73pVVZWWrX9NJ5/+Xe+2kjNzE1E1AA7TK9EVANCzuFwu79/jxo3TCzv/HrDcyzsb/NYLMoKcvwKgWyKgAIip1tZWZWRkSJJaWkJ325w28MSVY6v3fKbqPZ95t7cfOxKbCgJwJAIKgLja8ck3J7+ufbNOr37U7Hf7eYOztfamIr9zUF7e2aBN79bFrY4AEo+AAiBudnxyUFN/u927Pm/DTqX0TutU7vwh2Tp/SLbfNgIK0LMQUADEzUef+V8Bdn7pSKX1OVkFGYa+tzhBlQLgSAQUAAlz2ahCpaenc+l6AJ0wzBgAADgOAQVAt8UEg0DyIqAAAADHIaAAAADHIaAAAADHIaAAAADHIaAAAADHIaAAAADHIaAASDpr36zTGx8fDF0QQNKyPaC0tbXp7rvv1tChQ9WnTx8NGzZMv/rVr2SapreMaZqaN2+e8vPz1adPH5WUlOhvf/ub3VUB0E0My8nwW5+3YacuW7GVkAJ0Y7YHlPvvv1+PPPKIfvOb32j37t26//77tWjRIj300EPeMosWLdKyZcu0YsUKbdu2Tenp6ZowYYKOHGE6dQCdnT8kWyt/Mtq7Pu6MAZLkN+MxgO7F9oDy+uuvq7S0VD/60Y80ZMgQlZeX6+KLL9b27SdmMDVNU0uWLNFdd92l0tJSfec739Ezzzyjffv26fnnn7e7OgCSQFtbm/fvLVu2+K17nO3O8v7tPvKJzPbOZQB0H7YHlAsuuECbNm3S+++/L0n63//9X/3lL3/RD3/4Q0lSbW2tGhoaVFJS4r1PVlaWxowZo61btwbc59GjR9Xc3Oy3AEgOOz75phvmz7sbO91eWVmpESNGeNcnTpyoIUOGqLKyMmiZ+2dNU/2K67S9+sUY1RpAotk+m/GcOXPU3Nys4cOHy+Vyqa2tTb/+9a81ZcoUSVJDQ4MkKTc31+9+ubm53ts6WrBggX75y1/aXVUAMfbGxwc19bfbvetb3v/c7/YNGzZo6tSpfueoSVJ9fb3Ky8tVUVEhSSovL+9Upu3w51o852ZdePpAlZWVxegZAEgU21tQnnvuOa1atUqrV6/WW2+9paeffloPPPCAnn766Yj3OXfuXDU1NXmXuro6G2sMIFY6niMyv3Sk37kks2fP7hQ8JHm3zZgxQzNmzAhYxmPmzJkBu4QAJDfbW1Buv/12zZkzR1deeaUk6ayzztInn3yiBQsWaNq0acrLy5MkNTY2Kj8/33u/xsZGnX322QH3mZqaqtTUVLurCiDOLhtV6LdeX18ftKxpmtq7d2+IPZqqq6vTH//4R283MoDuwfYWlC+++EIpKf67dblcam9vlyQNHTpUeXl52rRpk/f25uZmbdu2TUVFRXZXB0APEKx7GEDysj2gXHrppfr1r3+tP/zhD/r444+1fv16Pfjgg/q3f/s3SZJhGJo5c6buvfde/f73v9e7776rf//3f1dBQYEmTZpkd3UAJJjvaJstW7YoLS1NpmmqurratsfwtMwC6D5s7+J56KGHdPfdd+uWW27RgQMHVFBQoBtvvFHz5s3zlpk9e7ZaW1t1ww036NChQ7rwwgv10ksvKS0tze7qAEig7dUvav/jd3rXJ06cKLfbraVLl6q0tFRut1v19fUBzzExDEODBg2SpKBlPMaOHWt/5QEklGF29al3qObmZmVlZampqUmZmZmJrg6AACorKzV5crkk/68YwzAkyW+EjiS/AGKljK+Wlhalp6d32t7a2qqMjIwuywCIn3B+v5mLB4Dt2traNGPGDHUMJ9I3IWPmzJkqLS1VRUWFCgoK/Mq43W5VVFSorKxMZWVlAcu4MvrHrP4AEo+AAsB2NTU1XY7AMc0To29qampUVlamXbt2eW+rqqpSbW2t37VNOpa5Y/HTyr/u4dhUHoAjEFAA2G7//v1hlXO5XN5t48aN81v38N02/JzRMlI6lwHQfRBQANjO9xpHdpSLlJU5fgA4EwEFgO2Ki4vldrslGQFvNwxDhYWFKi4ujlkdrMzxA8C5CCgAbOdyubR06dKAt3lG6CxZsiRgV44dKisrVV5e3ulKtZ45fggpgPMRUADERFlZmWYtfKTTaBvfETqx4BlB1NUcP8zfAzgfAQVAzIwe/0O/0TaBRujYLZwRRACci4ACIKZ8R9sEG6Fjp3BHEAFwJgIKgG7FKSOIAETH9rl4ACDe1r5Zp7Q+J2tYToZ3BFFXc/y43e6YjiACED0CCoCkc9pA/zl15m3YqZTeJyYbXXtTkZYuXary8nIZhhFwjp9YjiACYA+6eAAknfMGZ2vlT0Z71+eXjtT4bw2UJH14oCXo/D2xHkEEwD4EFABJ6bzB2d6/LxtVqAkj8/xutzLHDwDnIqAAsKS1tVWGYcgwDLW2tkZcJpbajx3RVWMGex/fyhw/AJyJgAIAAByHgALANjs+Oej9e+2bdXp5Z0PMHqvjRIDtXBkW6FYYxQPAFm98fFBTf7vdu+47ssZuGzZs0OzZs73rEydOVHZOvozRU2LyeADij4ACwBYfHmjxW59fOlJpfU5WQYah7y2297GmTp3a6RonBw80SBsfsPeBACQMAQVATFw2qlDp6ekxOVk20AXYpEDbACQrzkEBAACOQ0ABAACOQ0ABAACOQ0ABkHDp6ekyTVOmaSo9PT1omePHj8vtdnvn1AHQfRFQACQNl8ulpUuXSlKAkEJoAboTAgqApBJsIsD+uXnqf8nPElQrAHYjoABIOoEmAly2/jWdfPp3E1grAHYioABISh0nAkxhIkCgWyGgAAAAx+FKsgBiyjNCBwDCQQsKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAEva2tq8f2/ZssVv3cNsD10mlhL9+ADsQ0ABEFJlZaVGjBjhXZ84caKGDBmiyspK77bt1S9q/+O3dFkmlgI9/siRI7Vu3bou5/gB4EyGmYTj/5qbm5WVlaWmpiZlZmYmujpAt1ZZWany8vJOQ4U9c+FUVFRIkiZPLpcUvExZWZmt9WptbVVGRoYkadWqVZoyZWpcHx9A+ML5/SagAAiqra1NQ4YM0d69ewPebhiGBg0aJEldlnG73aqtrfW7+mu0fAPKoEGDVF9fH9fHBxC+cH6/6eIBEFRNTU3Q4CFJpmlq7969IcvU1dWppqYmFlWUpKDhJF6PD8B+BBQAQe3fv9+R+0rGxwcQHgIKgKDy8/Mdua9kfHwA4SGgAAiquLhYbrfbe7JpR57zO9xut6TgZQoLC1VcXByzeg4aNChoHaXYPz4A+xFQAATlcrm0dOlSSeoUADzrS5cu9ZbpyFNmyZIlMT1BddGiRQHr6BHrxwdgPwIKgC6VlZWpoqJCBQUFftvdbrd3+G5ZWZlmLXxEroz+QcvEUmlpacA6uvoO0KyFjzDEGEhCDDMGYInncydJVVVVuvjii/VWXZM+PNAiSXp5Z4P+/L8fa+/SK/zKxKrlwneYcUtLi9LT0/3qeMfip/Xsvn66v/xsXTn61JjUAUB4wvn97hWnOgFIcr5BY9y4cXqrrkmXrdjqV8ZI8S8T724V38cbfs5oGQ0fxvXxAdiHgAIgIp6Wk/HfGqgJI/MkSQUZhr63OJG1AtBdEFAARGXCyDxvF0pra2vcHjc9Pb3T5fcBdB+cJAsAAByHgAIAABwnJgGlvr5eU6dOVf/+/dWnTx+dddZZevPNN723m6apefPmKT8/X3369FFJSYn+9re/xaIqAGKo/dgRXTVmsAzDiGv3DoDuz/aA8o9//ENjx47VSSedpBdffFG7du3Sf/3Xf+mUU07xllm0aJGWLVumFStWaNu2bUpPT9eECRN05MgRu6sDAACSkO0nyd5///0qLCzUk08+6d02dOhQ79+maWrJkiW66667VFpaKkl65plnlJubq+eff15XXnml3VUCAABJxvYWlN///vcaNWqULrvsMuXk5Oicc87RY4895r29trZWDQ0NKikp8W7LysrSmDFjtHXr1kC71NGjR9Xc3Oy3AACA7sv2gPLRRx/pkUce0emnn66XX35ZN998s/7jP/5DTz/9tCSpoaFBkpSbm+t3v9zcXO9tHS1YsEBZWVnepbCw0O5qAwAAB7E9oLS3t+vcc8/Vfffdp3POOUc33HCDrr/+eq1YsSLifc6dO1dNTU3epa6uzsYaAwAAp7E9oOTn52vEiBF+284880x9+umnkqS8vBNXnGxsbPQr09jY6L2to9TUVGVmZvotAACg+7I9oIwdO1Z79uzx2/b+++9r8ODBkk6cMJuXl6dNmzZ5b29ubta2bdtUVFRkd3UAAEASsn0Uz6xZs3TBBRfovvvu0+WXX67t27fr0Ucf1aOPPipJMgxDM2fO1L333qvTTz9dQ4cO1d13362CggJNmjTJ7uoAAIAkZHtAOf/887V+/XrNnTtX8+fP19ChQ7VkyRJNmTLFW2b27NlqbW3VDTfcoEOHDunCCy/USy+9pLS0NLurAwAAklBMJgu85JJLdMkllwS93TAMzZ8/X/Pnz4/FwwMAgCTHXDwAAMBxCCgAAMBxCCgAAMBxCCgALElPT5dpmjJNU+np6RGXiSXfx0/rc3LcHx+AfQgoAADAcQgoAADAcQgoAADAcWJyHRQAPYPZ3ub9e8uWLbr44ovlcrkSWKPOXt7pP0v6sJwMnT8kO0G1AWAVAQVARLZXv6j9j9/pXZ84caLcbreWLl2qsrKyBNbshGE5GZKk6j2fqXrPZ363rb2piJACOJxhmqaZ6EqEq7m5WVlZWWpqamJmYyABKisrNXlyuST/rw/DMCRJFRUVjggpb3x8UB8eaPGuv7yzQdV7PtPCsrN05ehTE1gzoGcK5/ebFhQAYWlra9OMGTPUMZxIkmma3glBS0tLE97dc/6Q7E4tJR1bUwA4EyfJAghLTU2N9u7dG/R20zRVV1enmpqaONYKQHdDQAEQlv3799taDgACIaAACEt+fr6t5QAgEM5BAWCJ54TT9tTBys7J18EDgVtIDMOQ2+1WcXFxnGsIoDshoAAI6Y2PD+qyFVu9664LrpWev69TOc8oniVLliT8BFkAyY2AAiAkz1Dd8d8aqAkj8ySdpe3fPVUrF/9SnzV+05Lidru1ZMkSRwwxBpDcCCgALJswMs97/ZArR9+oX9x8lbKysiRJVVVVjrySLIDkxEmyACLmG0bGjRtHOAFgGwIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKgIBaW1tlGIYMw9CRL79IdHUA9DAEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAAA4DgEFAA9SvuxI7pqzGAZhqHW1tZEVwdAEAQUAADgOAQUAADgOL0SXQEAySs9PV2maSa6GgC6IVpQAACA4xBQAACA4xBQAPRYa9+s0xsfH0x0NQAEQEAB0GMMy8nwW5+3YacuW7GVkAI4EAEFQI9x/pBsrfzJaO/6uDMGSJI+PNCSqCoBCIKAAiCgtrY279/vvb1dZntbF6WTx3mDs71/l5yZm8CaAOgKAQVAJ5WVlRoxYoR3/f5Z01S/4jptr34xgbUC0JMQUAD4qaysVHl5uerr6/22tx3+XIvn3KzKysoE1QxAT0JAAeDV1tamGTNmdHnxtZkzZ/p1/wBALBBQAHjV1NRo7969XZQwVVdXp5qamrjVCUDPREAB4LV//35bywFApAgoALzy8/NtLedE3XV0EtDdEFAAeBUXF8vtdksyAt5uGIYKCwtVXFwc34rZhNFJQPIgoABQa2urDMNQr169dP/99wcsYxgnQsuSJUvkcrniWT1bMDoJSC4EFAB+SktLNWvhI3Jl9Pfb7na7VVFRobKysgTVLHKMTgKST8wDysKFC2UYhmbOnOndduTIEU2fPl39+/dXRkaGJk+erMbGxlhXBYBFo8f/UPnXPexdr6qqUm1tbVKGE4nRSUAyimlAeeONN/Tf//3f+s53vuO3fdasWXrhhRe0du1avfrqq9q3b1/SfvEB3ZWR8k03zrhx45KyW8eD0UlA8olZQGlpadGUKVP02GOP6ZRTTvFub2pq0uOPP64HH3xQ3//+93XeeefpySef1Ouvv66//vWvsaoOgB6sJ4xOArqbmAWU6dOn60c/+pFKSkr8tu/YsUNfffWV3/bhw4fr1FNP1datWwPu6+jRo2pubvZbAMAqz+gkz4m+nSX36CSgO4pJQFmzZo3eeustLViwoNNtDQ0N6t27t/r16+e3PTc3Vw0NDQH3t2DBAmVlZXmXwsLCWFQbQDflcrm0dOlSSQoaUpJ1dBLQXdkeUOrq6jRjxgytWrVKaWlptuxz7ty5ampq8i51dXW27BdAz1FWVqaKigoVFBT4bXf1HaBZCx/hPDjAYWwPKDt27NCBAwd07rnnqlevXurVq5deffVVLVu2TL169VJubq6OHTumQ4cO+d2vsbFReXl5AfeZmpqqzMxMvwUAwlVWVqZdu3Z51+9Y/LQG3fS4Ro//YQJrBSCQXnbv8Ac/+IHeffddv23XXnuthg8frjvuuEOFhYU66aSTtGnTJk2ePFmStGfPHn366acqKiqyuzoA4Me3G2f4OaNlNHyYwNoACMb2gNK3b199+9vf9tuWnp6u/v37e7dfd911uu2225Sdna3MzEz99Kc/VVFRkb773e/aXR0AAJCEbA8oVixevFgpKSmaPHmyjh49qgkTJujhhx8OfUcAANAjxCWgbN682W89LS1Ny5cv1/Lly+Px8AAAIMkwFw8AvzlotmzZonbmpAGQYAQUoIerrKzUiBEjvOsTJ07UTyeN1Rd/48rOABInIeegAHCGyspKlZeXd5rl9+CBBmnjAwmqFQDQggL0WG1tbZoxY0ancHJCoG0AED8EFKCHqqmp0d69e0OWq66uVnp6ehxqBADfoIsH6KH2799va7lk9vJO/3nAhuVk6Pwh2QmqDQCJgAL0WPn5+baWS0anDTzRMlS95zNV7/nM77a1NxURUoAEIqAAPVRxcbHcbrfq6+uDnIci9c/NV3FxcZxrFj/nDc7W2puK9OGBFu+2l3c2qHrPZ/rwQAsBBUggAgrQQ7lcLi1dulTl5eUyDKNDSDEkmfr3Wff4zV3THZ0/JLtTEOnYmgIg/jhJFujBysrKVFFRoQE5/jOJn5ydo4GTfs4svwAShhYUoIcrPPcipV65RFp6hSRpYPkv1GfoOTJSXBqWk5HYysVAenp60C4tAM5BQAF6uA8PtMhI+aYbZ9GtVyitz8mMZAGQUAQUAH4uG1XIdU8AJBznoAAAAMchoADo0VpbW2UYhgzDUGtra6KrA+BrBBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQAAOA4BBQASumdpme3fSLTNJmHB4AjEFAA9GhtbW3ev7ds2aJ2n3UAiUNAAdBjVVZWasSIEd71iRMn6qeTxuqLPa8nsFYAJAIK0CP5TpB35MsvEl2dhKisrFR5ebnq6+v9th880KDPnr9P26tfTFDNAEgEFAA9UFtbm2bMmCHTNAPcemLbM4t/6df9AyC+CCgAepyamhrt3bu3yzJ/b9yvmpqaONUIQEe9El0BAIi3/fv3Wyq3ruZdNZx8moblZOj8IdkxrhUAXwQUAD1Ofn6+pXJrd7fohdZ3T/x9UxEhBYgjungA9DjFxcVyu90yDCNwAcNQ/9x8PfjTKzX+WwMlSR8eaIljDQEQUIAeaMcnB71//3l3YwJrkhgul0tLly6VpE4hxTAMGZIeffg3urpoqCaMzEtADQEQUIAe5o2PD2rqb7d717e8/7kkaVhORqKqlBBlZWWqqKhQQUGB33a3262KigqVlZV5t7UfO6KrxgyWYRhqbW2Nd1WBHolzUIAepmNXxfzSkRo5OKdHnl9RVlamkpISZWVlSZKqqqp08cUXy+VyJbhmAAgoQA932ajCHj3/jm8YGTduHOEEcAi6eIAeyGz3n3+GC5IBcBoCCtDN+V7WvrW1VdurX9T+x2/x3j5x4kQNGTJElZWVCawlAPgjoAA9yIYNG7R4zs1qa/m73/b6+nqVl5cTUgA4BgEF6EFmz54tz1wzvjxz0sycOZPuHgCOQEABepCOM/f6Mk1TdXV1zD8jZnsGnICAAsCP1XlqACCWCCgA/FidpwYAYomAAvQggwYNkhR4/hnDMFRYWKji4uL4VioJMCwbiD8CCtCDLFq0KOB2z3w0S5Ys4UJlHby55Y8MywYSgIAC9CClpaWatfARuTL6+20PNP9MT5Geni7TNGWaZsAr6i6fN5Nh2UACEFCAbs63O2LLli0aNe5i5V/3sHdbVVWVamtre2Q4sYZh2UAiEFCAbsZ3iOzq1as1YsQI720TJ07UTyeN1ZcfvuHdxvwzkfEMy+7VqxczHAMxwGSBQDc2depU73/7HgcPNEgbH0hQjQDAGlpQgG6sYzj5emvc65Fs6LYBEs/2gLJgwQKdf/756tu3r3JycjRp0iTt2bPHr8yRI0c0ffp09e/fXxkZGZo8ebIaGxvtrkpcdZyQDUByqqys9OsWA5AYtgeUV199VdOnT9df//pX/elPf9JXX32liy++2O9He9asWXrhhRe0du1avfrqq9q3bx8n6AFIuMrKSpWXl3c5JQCA+DDMwG3Atvnss8+Uk5OjV199VePGjVNTU5MGDhyo1atXq7y8XJL03nvv6cwzz9TWrVv13e9+N+Q+m5ublZWVpaamJmVmZsay+pa1trYqIyNDktTS0hJwuCIQD77vRSt4v57Q1tamIUOGaO/evSHL9s/Nl86for9/fS4PxxCwJpzf75ifg9LU1CRJys7OliTt2LFDX331lUpKSrxlhg8frlNPPVVbt24NuI+jR4+qubnZbwFgReCrxsoIsr0Hq6mpsRROFi5cqGXrX9PJp4f+ZwpA5GIaUNrb2zVz5kyNHTtW3/72tyVJDQ0N6t27t/r16+dXNjc3Vw0NDQH3s2DBAmVlZXmXwsLCWFYb6NYMwwgWW3o0q5Mk5uTkKIVh2UDMxTSgTJ8+Xf/3f/+nNWvWRLWfuXPnqqmpybvU1dXZVEOg+/EdgTLp2lvlysj2u93tdmvlypXxrpbjWZ0kMS8vTxLz8wCxFrOAcuutt2rjxo2qrq6W2+32bs/Ly9OxY8d06NAhv/KNjY3eD35HqampyszM9FsAdNZxBMrzTz4k37PMPFeNLS0tTUDtnK24uFhut9s7L1EwY8eO1fbqF5mfB4gx2wOKaZq69dZbtX79er3yyisaOnSo3+3nnXeeTjrpJG3atMm7bc+ePfr0009VVFRkd3WApBXu0PVgI1DaWw96//ZcNTbU/DM9kcvl0tKlSyWpU0jxXd+4caMWz7k54Pw8kydP5lIDgE1sDyjTp0/XypUrtXr1avXt21cNDQ1qaGjQl19+KUnKysrSddddp9tuu03V1dXasWOHrr32WhUVFVkaweNUHec7obkX8dTW1qYZM2YEuTCbfzkEV1ZWpoqKChUUFPhtHzRokPfv2bNnq6v5eSSOM2AH2wPKI488oqamJl100UXKz8/3Lr/73e+8ZRYvXqxLLrlEkydP1rhx45SXl5fUTaMdm9Vp7kW8WR2B8tprr8WhNsmtrKxMu3bt8q5XVVVp586d3nUr10jhOAPRs30uHiuXVUlLS9Py5cu1fPlyux8+7jzN6h2ft2c69p46hT3iy+oIlGAj5eDPd/LEcePGhX1/jjMQPebiiUJXzepMx454CncECmKL4wxEj4AShVDN6p7p2GtqauJYK/RE4YxAQfg8JxUfP37861GJHGcg1ggoUbDarG61HBCprkagdCyHyPke5458jzvHGYgeASUAq8M7rTare8rZNeMxMyf3DOGODAs2AsWV0T8m9eupysrKNGvhI52Oq+9xZyQfED0CShRCNasbhqHCwkIVFxfHuWZIdpGODCs89yLNf+Zl73rJrMUadPMTenbbJ1zzxEajx/9Q+dc97F2fdO1Pdez4N4GEkXxA9AgoUbByYaclS5bQ3IuwBLvgmmdkWLAfvTc+PqjLVmzVLza+5922R4UyUlwalmN9dmOENiwnQ0bKN5/r5598SJ81+o/cCfV6AegaASVKwZrV3W43Q4wRtmhGhn14oEWSNO6MAd5t80tHau1NRTp/SHan8ojc+UOytfIno7ssw0g+IDoEFBsEurBTbW2tpXDC+STwZXVkWK9evYK+X0rOzPX+fdmoQsJJBKxMBXDe4NDHlZF8QOQIKDbpeGEnunUQCUZ8dU+8rkD4CCgBxGpeHbv2y7w/3ZfVkWFILryuQPgIKB3YOa+Ob/fN6tWrA+53w4YNCaufVQyPjh8rI8PgRIzkA+xGQPER6egJK6ZOnRpwv1OnTnVE/eAMVkaGITkwkg+IDgHla7GeVyfYfn23d9Vdw7w/PUdXI8OeeeYZ73qw90tan5NDnuCJ6PmeSBvowm2M5AOiQ0D5mhPm1emqu8YJ9UP8BBoZ9uCDD2rOnDnebVwMzDk6XrgtnJF8AAIjoHzNKfPqBOquaW1t1fjx4y3dv6v6JeM5IMlYZ7v4dgv84x//0OWXX96pe29vfb0mTy7Xbff/t17e2dBxF4gj3wu3dfeRfD35c4n4IaB8Ldx5dTqyct0EK6LtrmG0QPc0e/bsgN17Mk1Jppb9+i69svtEQOGqsQC6AwLK12Ixr04056sE665JxLw/DI9OvI4tJx21Hf5c04Z+yVVjHWLtm3Vas/1TvfHxwURXBUhaBJSv2T2vTsfhwJEI1l1jpX52NcHaNaw5EcOjEyXQsY9Hk/jQtC8IJwnSsdVq3oadmlP5ri5bsZWQAkSIgOLDrnl1gg0HDleg7pqVK1fGbd4fu4Y1Mzw6PvLy8hJdhR6r49w880tHavy3Bkr6Zo4kAOEhoHQQzbw6UtfDgX0NGDDAcneNb1fIKaeconfffTei+oXTxWLXsOZo90O30AmDBg0KeR2UsWPHxqk2CMR3bp7LRhVqwsjuGxj5XCIeCCgBRDOvTqjhwB433nijpNDdNYG6Rs4666yw6xduF4tdw5qj2U9P6hYKxPfE62XLlknq+mJt3XnUSDKw60R5p+vpn0vEDwHFZlaHIQ8bNixkd1KwrpF9+/Z12l9X5zhs2LChyy6WQPeLdti1pz6RDo8Op1sonud8JGp4ZVlZmRY+/KT6Dcj1297x4mBwlvZjR3TVmMFdvl+cNmS3q/o4pbvWaccMsUFAsZnVYb55eXlddidZ6RqRgo8U8t1uZT+RPg+7hjX77ifZr5obqPk72ibxNz4+qIc/HqD0q5d6tw0s/4Xyrv1NVPsFrEr2zyWSDwElgJ37mrx/e4YLhlo8Z+qHGq7s4TlfIFh3ktWuotdee63Ttg0bNvg1wX7++edB7x8spNg57Lqr8ycC7SeZr5obqPk7NzdXQ4cO9dsWbpO450TL7w3/pgXl8u+comPP3RbVfgGrkvlzieTUK9EVcKLyR7Z6/563YadSeqdZut9PLhyqf8rJ0GW33qXFc27WiRlOfQPAN+uVb+9TWp9DOvLlF95b175Zp5GDc3T+kGzLXSwNDZ2vHjp16tSQJ+mG4hl27ekC8t1fuMOuFy1apKlTp1rej1Ou6hsuT/N3x2P/97//vVNZT5N4uKOvSs7M1aqv/374npmdHivS/SL2fD/fyShZP5dIXrSgBFBxc5H37/mlI7Ww7Kwul59ceOK/49/+pVZzKt9V5T/cGjhprlwZ/l9Evuue6yTM27DTb9tlK7bq3o27tPuQtZdmT5NLa7Z/qu0ffebdFmk46dhFEM2w644jj5577jnL+7HabXT11VertbXVti4VK/3awfZrdfSWh+dkysmTJ6u5uTlgmTc+PuhtoQt0GXua2p1tWE6GzPZvXoPZv/mdypb8udN7rON7qrm5OWbnVEXzHg+32zeW54nE6vMNhzGTUFNTkynJbGpqSnRVvLbX/t18dtsnfsvjm/7vxHXIJfOOxU+bq17/qFOZJzfv9pa5e+0b5uA7NpqD79honnr7BtPVd4D3tkCLq+8A89TbN5gDJ/3cdGX077Ks1SUnr8Bct26d33PzHG9JZlVVlXn8+PEuj8W6devMQYMG+e3X7Xab//M//2NpP8ePHzfdbrdpGEbAOvpuX7VqVafH6t+/v5mdnd3p8Ts+r45aWlq85VtaWiw/r3Xr1pnV1dURH/OqqqqA7yfPe8F3eXXnp5b3W11d3eXzRWytW7fOHJib7/+Z9fmctrS0BHxP+a573oeh3ptWRfMet/K5LCws9H6u7apzOHWM5rkjPsL5/SagxJCVD0THMr5BZ9bCFaZkfL34fhmc2DZr4QqfMtGHE9/9z1q4ImCIenLz7k4ha3vt373PZ926dQG/wAzD8Nse6gvCs5+O++q4LdiXZbDH7+pLrKvXK9TzmjlzZsTH+4knnuhUl2e3fWIOvmOj+eMntvkdZ986hlpWr17d5TFG7AR7v/guq1atCvqeSkRACfUeX7duXcjPpe/nKxaBwEodI3nuiB8Cis0CvbEjCR+RCPTfQmFhod9/NPaGE5//Smb8zhx8x0azcFbFN489qyLgf/Z3r33DUvDxDTodg0/HH+BArSOFhYV+LTGRLMHeN8FaikIdZ8MwzIEDB0Zcn4ULF5pbPzjgF/p+/MQ2c/AdG81nt33SqZ5WW2toQUkMq5/LAQO6biGVZG7cuNE8fvx42K2YwQTaj5XQ69s60tV3UqjH6iic70grn0PfFhwrzz1WIv3N6AkIKDZLZEAxzeAfrGi6Fawsdyx+OmQLyq9e2NkpxFhZCmdVBAw+vkHnyc27O3WTzbzvYfOUgXlRPa9AXSp2dN8MHDjQcotOx8XVd4A5cNLPOwU/39Ypj3Cb2hFfdn8uI+2y7CjYe3zVqlWW6+IJvaF+7K12w4TzHRlNMI+0WyhSBJTgCCg2S3RACbaf1atXW/5iCdQsG2rxdD2Eeh7ba//uF2KsLB1bUDzn33QMLb7r/S/5mS1f+Dfe/YBfyAraTRZm983MmTMjOs6+i2/XWqBw4hFOUzviK5zPZaRLuK+z1a7XUIun2zDariKPcL4jrR7Xjl2b0XQLRYqAEhwBJUaPJ33z30KimwtNM7z/1AoLC8P6T0k60fVg9Xn5Ho9QS7D/wrZ+cKBTi4nverQtJ54la9w089TbN1g8GdkwM0/JtrTfux9eY85auKJTPVPS+popaRmh9xFmy4fVpnbEV6xbNn1/YK28X+zsCvZ8doN9J4XbDRPO92gkLSjRdgtFKtG/GU5GQLFRoB8Bu5pcAwkneYdq6vcsnn5s3/0MGjTI0n9OVp5XoGMU6svA6nG1st9IluycfHPWwhXm3Q+vsWV/nhFVg+/YaLpn/M67fWD5L8xTb99g3rnmr2Htz+p/WHzpOY/Vz6VdS6jQYEdgsjpCJ5zHCnSOWVffN1aPq+/vQrzO1wp17pydvxnJ3hJDQLGJlTPxfT/AdjQXhvvm66qpv6v9BLtfuM/L6jEKNBIgnOMazSgZ+/bb9YiqYKOewh19E84XT7J/WXVXVj9fdiyhul2i7XIK57st3G7nSL9vrH7fRdotFC7fYx/r34xk/8wTUGwQSbOoHc2FkfxHHOwEsFBv4khaPiI9RtGMPAp3lEyg/1bs2O+dd95pFhQUBHxevkK1gllZrL72yf5l1Z1Z/XxFu4Q6cTXaFpRwug1j0Vpj5bj6rkfy3KNtQQmnizuc5xrqsZKx1ZSAYoNoPmiRvtmjOdO845vWd72rH65Izh0J9xj5nssSzXG1Mkom0LBMO/brWXwDSjhDJyP5D4uLTyW/SH+4wv1xi+YCa1Y/u1bY2b3V1feo73EN9I9DpBeXi4QdQdTqb0a8RyPFAgHFBtE0i0bSXBjtmeaR/lBFcuGvcFsDfC9EFs1xDTZKJtSFrSLdb7AvtK4eK9IL0AV7LEbkJLeOwTTUhQdl9YfdQpdpoIsaxmokny87u7daWloiCvyRXlwu3Nc02MX2wl2s/GbYPRopUf/chPP7zVw8QVidd8KO+ybLNOae5xVuPfLy8jrtIxKlpaUB5wYaNGhQxPvsar+B+L5GvvPweISaFyTY/EZdPZYTXntEb+XKlQHno1q5cqV3fdXKlRqY6/8ZOblvP6WkZfhtc2X01xVzFutI/jm6/uZbg353+G7/+c9/Hvbnz/PZjdV7PJQtW7bo2LFjXT52qO/NYJ9vK3OKBeL7+OHMv9WVUK9LLH4jIpnPKO5iFJJiKp7noISTjiNtLrSjn9SOFhQrzaCRNGf6vk52HFcr3VnhPq+O+7WyVFVVRdzkGu5jcVXY5NTxcxnqSq7Byrzy/2q92675xQrviLHcq+4L633Ub0Cu9+++WaeYoabJaGpqitt7PNiSkpLS6bFjcXE5K+w+r8jqb4bd59IksquILh6bxKs5Pl5nmnfFSjOoXaOaoj2udjQ3h9qvlcXTNRTJeyHcx2Jene4h0gs8dtzmmbPr1vnLIv5x7PrCh6Hm+uo8eq2riVDtXMIdgWfl4nJWhPP9F+13ZEd2/kYk4sJ1vggoNrJ6vY5oLpDllLlV7Jz3J9TxiOa4RjIUO9ScPh1bYqwsXY0ACvWfUbgBhRaU7sGugOIRzUnnj2/6v4AXFvRMuxDOjOqBrv/ju27/Yv0CilLniz4GmvT02W2f+LVUhTMfV6gl1DWegl3M7eDBg7Z+T1h9Hp7BBrFAQInR43X1RormxXTS3Cp2DFW0ejziebXFjvtdu3ZtVNPKWx2ebPULwwmvPZJPNKNmurqi87PbPrF8EcNgV1D2XAxx1esfmdk5+ab9s67LTOmTablssPmuPMvAST83XRn9Az6HaC7oePt/PWGuev0jv4A06dqfmqcMzPUrl5HZz0zP7Of//Dp0bwVb7O4q8v0+tBsBxWbxmFfBKXOr2HGxp2guMharM8utnHkfzpn/4TYvd8Uprz2SU6SjZuy6wJuVrk4rF1iLZPnhldeZgS+gGGgJ3i0VvCvrxNJ3VGnEdfRMhBruhKrWl9DdbZ4lvC7B2Hz/EFBsFq8fUifMrWJHU7LTA0pXJ7mFOiHY83rE46Q15tVBMMePHzerq6vN1atXm9XV1QFbBK2+N6P9zFvt6gx1QclIlurq6rBOXO3Y0nD8+HHzz3/+c4gLO4bXnSTJzM75pjXJ050Uq3NyUjL6m1ljp5gDLr3dzL3qPm+Xm+9y6u0bzNyr7jMzi64Ie/92t+ASUJKYU68S2F26IsINFsFej1gcD6e+9nCWYD/0vudUDRgwIOr3Ziy6OiP9PIV6DpGMigt3RI7VCzp29dm1eyLJQF1Fnm4p39ahaCdatfMcOAJKEnPylUG7Q1dEuGfDd/V6xPICUE577eEMVi/K5unGjPa9aXdXZySfp2DPNRYj8KzcJ9ixsPLZjXZepEDHoqvjY9fIIztHERJQEDPJ3hVB1wySVagRGB1bFex6bya6q9Plcll+DnZ0S3X1HKI9pna2oIR6DgMGDLA0J5mdr58VBBTEVDJ3RdA1g2QVSRiw672ZyK5O36G2oZ6Dnd1SwZ5DNMfUjrmKInkO0TwW56CEiYCCaHSHrir0PE64oGMgTvs82dUtFavnEM1cRZE8h1BLqK6rRI7iYS4e9DjB5gqJdG4OIB6szqMTzXxXkXDa5ylUfUpLSy3vy+1263e/+52ys7P17LPPavPmzVHPWROsfv3791d2drbfNpfL1ak+4T6HUOyeq8hWtkajMP3mN78xBw8ebKamppqjR482t23bZul+tKDADnTNIJk4fSSd0z5P0Y4Y2rhxo7l27dpO5/3YdQEzKxeqDNa9ZVdXkV1dV5E8b0d38axZs8bs3bu3+cQTT5g7d+40r7/+erNfv35mY2NjyPsSUAD0RE7rTklW0cw95pRjbUdXUSKeQ1IElNGjR5vTp0/3rre1tZkFBQXmggULQt6XgAKgp2LkmD2imXss0a1VXT0HK0si3y/h/H4bpmmaFnqCbHXs2DGdfPLJqqio0KRJk7zbp02bpkOHDmnDhg1+5Y8ePaqjR49615ubm1VYWKimpiZlZmbGq9oA4AhtbW2qqanR/v37lZ+fr+Li4k7nKyC0YMdx8+bNGj9+fMj7V1dX66KLLop9RbvgeQ719fWaNWuWPv/8cwX7Wc/OztZzzz2niy66KGHvl+bmZmVlZVn6/e4Vpzr5+fzzz9XW1qbc3Fy/7bm5uXrvvfc6lV+wYIF++ctfxqt6AOBoLpcr4T+M3UGw47h//35L97daLpZ8n0OfPn1UXl4uwzD8QophGJKkxx57TD/4wQ8SUc2IJMUonrlz56qpqcm71NXVJbpKAIBuyqkjpkLxjBAaNGiQ33ZHjMiJQEJaUAYMGCCXy6XGxka/7Y2NjcrLy+tUPjU1VampqfGqHgCgBysuLpbb7VZ9fX3A7hLDMOR2u1VcXJyA2nWtrKxMpaWl3aILMCEtKL1799Z5552nTZs2ebe1t7dr06ZNKioqSkSVAACQdKLbZOnSpZK+6R7x8KwvWbLEsT/6nm6fq666KqHnm0QrYV08t912mx577DE9/fTT2r17t26++Wa1trbq2muvTVSVAACQ1P26S5JRQrp4JOmKK67QZ599pnnz5qmhoUFnn322XnrppU4nzgIAkAjdqbskGSVkmHG0whmmBAAAnCGc3++kGMUDAAB6FgIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwHAIKAABwnIRd6j4anovfNjc3J7gmAADAKs/vtpWL2CdlQDl8+LAkqbCwMME1AQAA4Tp8+LCysrK6LJOUc/G0t7dr37596tu3b6epsKPV3NyswsJC1dXVMc9PDHGc44PjHB8c5/jgOMdPrI61aZo6fPiwCgoKlJLS9VkmSdmCkpKSIrfbHdPHyMzM5AMQBxzn+OA4xwfHOT44zvETi2MdquXEg5NkAQCA4xBQAACA4xBQOkhNTdU999yj1NTURFelW+M4xwfHOT44zvHBcY4fJxzrpDxJFgAAdG+0oAAAAMchoAAAAMchoAAAAMchoAAAAMfpkQFl+fLlGjJkiNLS0jRmzBht3769y/Jr167V8OHDlZaWprPOOktVVVVxqmlyC+c4P/bYYyouLtYpp5yiU045RSUlJSFfF5wQ7vvZY82aNTIMQ5MmTYptBbuJcI/zoUOHNH36dOXn5ys1NVVnnHEG3x0WhHuclyxZom9961vq06ePCgsLNWvWLB05ciROtU1OW7Zs0aWXXqqCggIZhqHnn38+5H02b96sc889V6mpqfqnf/onPfXUUzGvp8weZs2aNWbv3r3NJ554wty5c6d5/fXXm/369TMbGxsDln/ttddMl8tlLlq0yNy1a5d51113mSeddJL57rvvxrnmySXc43z11Veby5cvN99++21z9+7d5o9//GMzKyvL3Lt3b5xrnlzCPc4etbW15qBBg8zi4mKztLQ0PpVNYuEe56NHj5qjRo0yJ06caP7lL38xa2trzc2bN5vvvPNOnGueXMI9zqtWrTJTU1PNVatWmbW1tebLL79s5ufnm7NmzYpzzZNLVVWVeeedd5qVlZWmJHP9+vVdlv/oo4/Mk08+2bztttvMXbt2mQ899JDpcrnMl156Kab17HEBZfTo0eb06dO9621tbWZBQYG5YMGCgOUvv/xy80c/+pHftjFjxpg33nhjTOuZ7MI9zh0dP37c7Nu3r/n000/HqordQiTH+fjx4+YFF1xg/va3vzWnTZtGQLEg3OP8yCOPmKeddpp57NixeFWxWwj3OE+fPt38/ve/77fttttuM8eOHRvTenYnVgLK7NmzzZEjR/ptu+KKK8wJEybEsGam2aO6eI4dO6YdO3aopKTEuy0lJUUlJSXaunVrwPts3brVr7wkTZgwIWh5RHacO/riiy/01VdfKTs7O1bVTHqRHuf58+crJydH1113XTyqmfQiOc6///3vVVRUpOnTpys3N1ff/va3dd9996mtrS1e1U46kRznCy64QDt27PB2A3300UeqqqrSxIkT41LnniJRv4NJOVlgpD7//HO1tbUpNzfXb3tubq7ee++9gPdpaGgIWL6hoSFm9Ux2kRznju644w4VFBR0+lDgG5Ec57/85S96/PHH9c4778Shht1DJMf5o48+0iuvvKIpU6aoqqpKH3zwgW655RZ99dVXuueee+JR7aQTyXG++uqr9fnnn+vCCy+UaZo6fvy4brrpJv385z+PR5V7jGC/g83Nzfryyy/Vp0+fmDxuj2pBQXJYuHCh1qxZo/Xr1ystLS3R1ek2Dh8+rGuuuUaPPfaYBgwYkOjqdGvt7e3KycnRo48+qvPOO09XXHGF7rzzTq1YsSLRVetWNm/erPvuu08PP/yw3nrrLVVWVuoPf/iDfvWrXyW6arBBj2pBGTBggFwulxobG/22NzY2Ki8vL+B98vLywiqPyI6zxwMPPKCFCxfqz3/+s77zne/EsppJL9zj/OGHH+rjjz/WpZde6t3W3t4uSerVq5f27NmjYcOGxbbSSSiS93N+fr5OOukkuVwu77YzzzxTDQ0NOnbsmHr37h3TOiejSI7z3XffrWuuuUY/+clPJElnnXWWWltbdcMNN+jOO+9USgr/g9sh2O9gZmZmzFpPpB7WgtK7d2+dd9552rRpk3dbe3u7Nm3apKKiooD3KSoq8isvSX/605+Clkdkx1mSFi1apF/96ld66aWXNGrUqHhUNamFe5yHDx+ud999V++88453+dd//VeNHz9e77zzjgoLC+NZ/aQRyft57Nix+uCDD7wBUJLef/995efnE06CiOQ4f/HFF51CiCcUmkwzZ5uE/Q7G9BRcB1qzZo2ZmppqPvXUU+auXbvMG264wezXr5/Z0NBgmqZpXnPNNeacOXO85V977TWzV69e5gMPPGDu3r3bvOeeexhmbEG4x3nhwoVm7969zYqKCnP//v3e5fDhw4l6Ckkh3OPcEaN4rAn3OH/66adm3759zVtvvdXcs2ePuXHjRjMnJ8e89957E/UUkkK4x/mee+4x+/btaz777LPmRx99ZP7xj380hw0bZl5++eWJegpJ4fDhw+bbb79tvv3226Yk88EHHzTffvtt85NPPjFN0zTnzJljXnPNNd7ynmHGt99+u7l7925z+fLlDDOOlYceesg89dRTzd69e5ujR482//rXv3pv+973vmdOmzbNr/xzzz1nnnHGGWbv3r3NkSNHmn/4wx/iXOPkFM5xHjx4sCmp03LPPffEv+JJJtz3sy8CinXhHufXX3/dHDNmjJmammqedtpp5q9//Wvz+PHjca518gnnOH/11VfmL37xC3PYsGFmWlqaWVhYaN5yyy3mP/7xj/hXPIlUV1cH/L71HNtp06aZ3/ve9zrd5+yzzzZ79+5tnnbaaeaTTz4Z83oapkk7GAAAcJYedQ4KAABIDgQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOAQUAADgOP8fyoeT2U0GM4sAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m.visualize()" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The result is identical when we fit with `iminuit.cost.ExtendedBinnedNLL`." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Migrad
FCN = 92.87 (chi2/ndof = 1.0) Nfcn = 190
EDM = 1.67e-05 (Goal: 0.0002)
Valid Minimum No Parameters at limit
Below EDM threshold (goal x 10) Below call limit
Covariance Hesse ok Accurate Pos. def. Not forced
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
Name Value Hesse Error Minos Error- Minos Error+ Limit- Limit+ Fixed
0 x0_n 990 40 0
1 x0_mu 0.4964 0.0019 0 1
2 x0_sigma 0.0487 0.0016 0
3 x1_n 629 29 0
4 x1_mu 0.97 0.14 0
\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
x0_n x0_mu x0_sigma x1_n x1_mu
x0_n 1.3e+03 0.000424 (0.006) 0.0123 (0.209) -264 (-0.252) -0.369 (-0.076)
x0_mu 0.000424 (0.006) 3.44e-06 -9.8e-08 (-0.032) 0.000293 (0.005) -1.64e-05 (-0.065)
x0_sigma 0.0123 (0.209) -9.8e-08 (-0.032) 2.64e-06 -0.0129 (-0.271) -1.64e-05 (-0.075)
x1_n -264 (-0.252) 0.000293 (0.005) -0.0129 (-0.271) 851 0.337 (0.085)
x1_mu -0.369 (-0.076) -1.64e-05 (-0.065) -1.64e-05 (-0.075) 0.337 (0.085) 0.0183
" ], "text/plain": [ "┌─────────────────────────────────────────────────────────────────────────┐\n", "│ Migrad │\n", "├──────────────────────────────────┬──────────────────────────────────────┤\n", "│ FCN = 92.87 (chi2/ndof = 1.0) │ Nfcn = 190 │\n", "│ EDM = 1.67e-05 (Goal: 0.0002) │ │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Valid Minimum │ No Parameters at limit │\n", "├──────────────────────────────────┼──────────────────────────────────────┤\n", "│ Below EDM threshold (goal x 10) │ Below call limit │\n", "├───────────────┬──────────────────┼───────────┬─────────────┬────────────┤\n", "│ Covariance │ Hesse ok │ Accurate │ Pos. def. │ Not forced │\n", "└───────────────┴──────────────────┴───────────┴─────────────┴────────────┘\n", "┌───┬──────────┬───────────┬───────────┬────────────┬────────────┬─────────┬─────────┬───────┐\n", "│ │ Name │ Value │ Hesse Err │ Minos Err- │ Minos Err+ │ Limit- │ Limit+ │ Fixed │\n", "├───┼──────────┼───────────┼───────────┼────────────┼────────────┼─────────┼─────────┼───────┤\n", "│ 0 │ x0_n │ 990 │ 40 │ │ │ 0 │ │ │\n", "│ 1 │ x0_mu │ 0.4964 │ 0.0019 │ │ │ 0 │ 1 │ │\n", "│ 2 │ x0_sigma │ 0.0487 │ 0.0016 │ │ │ 0 │ │ │\n", "│ 3 │ x1_n │ 629 │ 29 │ │ │ 0 │ │ │\n", "│ 4 │ x1_mu │ 0.97 │ 0.14 │ │ │ 0 │ │ │\n", "└───┴──────────┴───────────┴───────────┴────────────┴────────────┴─────────┴─────────┴───────┘\n", "┌──────────┬───────────────────────────────────────────────────┐\n", "│ │ x0_n x0_mu x0_sigma x1_n x1_mu │\n", "├──────────┼───────────────────────────────────────────────────┤\n", "│ x0_n │ 1.3e+03 0.000424 0.0123 -264 -0.369 │\n", "│ x0_mu │ 0.000424 3.44e-06 -9.8e-08 0.000293 -1.64e-05 │\n", "│ x0_sigma │ 0.0123 -9.8e-08 2.64e-06 -0.0129 -1.64e-05 │\n", "│ x1_n │ -264 0.000293 -0.0129 851 0.337 │\n", "│ x1_mu │ -0.369 -1.64e-05 -1.64e-05 0.337 0.0183 │\n", "└──────────┴───────────────────────────────────────────────────┘" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "def total(xe, x0_n, x0_mu, x0_sigma, x1_n, x1_mu):\n", " return signal(xe, x0_n, x0_mu, x0_sigma) + background(xe, x1_n, x1_mu)\n", "\n", "c = ExtendedBinnedNLL(n, xe, total)\n", "m = Minuit(c, x0_n=500, x0_mu=0.5, x0_sigma=0.05, x1_n=100, x1_mu=1)\n", "m.limits[\"x0_n\", \"x1_n\", \"x0_sigma\", \"x1_mu\"] = (0, None)\n", "m.limits[\"x0_mu\"] = (0, 1)\n", "m.migrad()" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAGdCAYAAAA44ojeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/P9b71AAAACXBIWXMAAA9hAAAPYQGoP6dpAAA/SklEQVR4nO3de3xU9Z3/8fdkMAnkxiWQCwmEuhQBr9WKgNBa8oAu4sIvZNUKFl1W7IpdAusFUJS1FpS6FlSU4lK1D6VU0xEqsrQWq6SC2KW69VZqK9Ukksg1IYMJODm/P3DGTDIzOTNzZuZM5vV8PObxyDnznTOf880M+fC9OgzDMAQAAGAjaYkOAAAAoDMSFAAAYDskKAAAwHZIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANhOr0QHEIn29nZ98sknysnJkcPhSHQ4AADABMMwdPz4cRUXFystLXQbSVImKJ988olKS0sTHQYAAIhAbW2tSkpKQpZJygQlJydH0ukbzM3NTXA0AADAjObmZpWWlvr+joeSlAmKt1snNzeXBAUAgCRjZngGg2QBAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAC243a75XA45HA45Ha7Ex0OgAQgQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUAAAgO2QoAAAANsJO0HZuXOnrrjiChUXF8vhcGjz5s1+zxuGobvuuktFRUXq3bu3ysvL9cEHH/iVOXLkiGbNmqXc3Fz17dtXc+fOVUtLS1Q3AgAAeo6wExS3263zzjtPa9euDfj8qlWr9NBDD2ndunXas2ePsrKyNGXKFLW2tvrKzJo1S++++65eeuklbd26VTt37tS8efMivwsAANCjOAzDMCJ+scOh559/XjNmzJB0uvWkuLhY//Ef/6FbbrlFktTU1KSCggI9+eSTuvrqq/X+++9r1KhR+sMf/qCLLrpIkrR9+3ZNnTpVdXV1Ki4u7vZ9m5ublZeXp6amJuXm5kYaPgCb8n7HJWnbtm2aPHmynE5ngqMCEK1w/n5bOgZl//79amhoUHl5ue9cXl6exowZo927d0uSdu/erb59+/qSE0kqLy9XWlqa9uzZE/C6bW1tam5u9nsA6JlcLpdGjRrlO546darKysrkcrkSGBWAeLM0QWloaJAkFRQU+J0vKCjwPdfQ0KBBgwb5Pd+rVy/179/fV6azlStXKi8vz/coLS21MmwANuFyuVRZWan6+nq/8/X19aqsrCRJAVJIUsziWbJkiZqamnyP2traRIcEwGIej0cLFixQoF5n77mqqip5PJ54hwYgASxNUAoLCyVJjY2NfucbGxt9zxUWFurTTz/1e/7zzz/XkSNHfGU6y8jIUG5urt8DQM9SU1Ojurq6oM8bhqHa2lrV1NTEMSoAiWJpgjJs2DAVFhZqx44dvnPNzc3as2ePxo4dK0kaO3asjh07pr179/rKvPzyy2pvb9eYMWOsDAdAEjlw4ICl5QAkt17hvqClpUV//etffcf79+/XW2+9pf79+2vIkCGqqqrSvffeq+HDh2vYsGFatmyZiouLfTN9Ro4cqW9/+9u64YYbtG7dOp06dUo333yzrr76alMzeAD0TEVFRZaWA5Dcwp5m/Morr+iyyy7rcn7OnDl68sknZRiG7r77bq1fv17Hjh3TpZdeqkcffVRf/epXfWWPHDmim2++WS+88ILS0tI0c+ZMPfTQQ8rOzjYVA9OMgZ7H4/GorKxM9fX1AcehOBwOlZSUaP/+/Uw5BpJUOH+/o1oHJVFIUICeyTuLR5JfkuJwOCRJ1dXVqqioSEhsAKKXsHVQACAaFRUVqq6u7tLdW1JSQnICpBhaUADYDivJAj0TLSgAbMvtdsvhcMjhcMjtdgcs0zEZmThxIskJkIJIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwHRIUAABgOyQoAADAdkhQAACA7ZCgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUADElcfj8f28c+dOv2MA8CJBARA3LpdLo0aN8h1PnTpVZWVlcrlcCYwKgB2RoACIC5fLpcrKStXX1/udr6+vV2VlpV+SkpWVJcMwZBiGsrKy4h0qABsgQQEQcx6PRwsWLJBhGF2e856rqqqiuweADwkKgJirqalRXV1d0OcNw1Btba1qamriGBUAOyNBARBzBw4cMFXusssuk9vtDvic2+2Ww+GQw+EIWgZAz0GCAiDmioqKEh0CgCRDggIg5iZMmKCSkhI5HI6Azwc7DyB1kaAAiDmn06k1a9ZI6pqMOBwOv8GzI5dtj2tsAOyJBAVAXFRUVKi6ulrFxcV+50tKSjRg2i0JigqAXZGgAIibiooKvffee77jbdu2af/+/eoz/JIERgXAjkhQAMSV0+n0/Txx4kS/YwDw6pXoAAD0bG63W9nZ2ZKklpYWU68pW/xil3PvLvumlWEBsDlaUAAAgO2QoAAAANuhiwdAXHWcRjxy2XalpWcmMBoAdkULCgAAsB0SFAAAYDskKAAAwHZIUAD0WOyADCQvEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYjuUJisfj0bJlyzRs2DD17t1bZ555pn7wgx/IMAxfGcMwdNddd6moqEi9e/dWeXm5PvjgA6tDAdCDdN7DB0DPZnmCcv/99+uxxx7TI488ovfff1/333+/Vq1apYcffthXZtWqVXrooYe0bt067dmzR1lZWZoyZYpaW1utDgcAACQhyxOUXbt2afr06br88stVVlamyspKTZ48WW+88Yak060nq1ev1p133qnp06fr3HPP1c9+9jN98skn2rx5s9XhAEgCRrvH93Nr7Tt+x8HKeDxdywDoOSxPUMaNG6cdO3boL3/5iyTp//7v//T73/9e//iP/yhJ2r9/vxoaGlReXu57TV5ensaMGaPdu3cHvGZbW5uam5v9HgB6hhP7dunAhpt8xwerl6t+3Vyd2LcrZJmysjK5XK64xgogfixPUBYvXqyrr75aZ511ls444wxdcMEFqqqq0qxZsyRJDQ0NkqSCggK/1xUUFPie62zlypXKy8vzPUpLS60OG0CMhBo7cuKD13Vw8wp5Wg77nfccP6SDm1foxL5dOrFvV8Ay9fX1qqysJEkBeijLE5Rnn31WzzzzjDZu3Kg//vGPeuqpp/TAAw/oqaeeiviaS5YsUVNTk+9RW1trYcQAEuXYK0+EfP7wb3+iI7/9ScDnvAPvq6qq6O4BeqBeVl/w1ltv9bWiSNI555yjjz76SCtXrtScOXNUWFgoSWpsbFRRUZHvdY2NjTr//PMDXjMjI0MZGRlWhwogwTq3inTW3s3zhmGotrZWv/nNb3zdyAB6BstbUE6cOKG0NP/LOp1Otbe3S5KGDRumwsJC7dixw/d8c3Oz9uzZo7Fjx1odDoAUEKx7GEDysrwF5YorrtAPf/hDDRkyRKNHj9abb76pBx98UP/yL/8iSXI4HKqqqtK9996r4cOHa9iwYVq2bJmKi4s1Y8YMq8MBYDNp6ZkaevtWtX78JzX+fKkl1/S2zALoOSxPUB5++GEtW7ZMN910kz799FMVFxfrxhtv1F133eUrc9ttt8ntdmvevHk6duyYLr30Um3fvl2ZmZlWhwMgwTpPD+497AI50pzKKBktZ06+PMcPBX1tWvYAOdR9V9D48eOtCheATTiMjku8Jonm5mbl5eWpqalJubm5iQ4HQBAul0tXzpnnl2A4c/LVf9I89RkxzjdDJ5iBM063sAQq43A4fANlW1palJWV1aWM2+1WdnZ2yDIA4iecv9/sxQMgJlwulyorK0NOIe4zYpwGzlgqZ/YAvzLOnHwNnLFUfUaMC1pm8ODBMb8HAIlDggLAch6PRwsWLFCoBtojO9bLaPeoz4hxKpr7qO/8wMrlGvy9DeozYpzvXKAy7777bmyCB2ALJCgALFdTU6O6urqQZTzHD6mt7nSS4Uhz+s5nlp7td+zVuYzT2bUMgJ6DBAWA5Q4cOGCqnKflaIwjAZCsSFAAWK7jIoyhOLP7xTSOjivM7ty5kxVngSRCggLAchMmTFBJSYkcDkfQMs6cfGWUjI5ZDC6XS6NGjfIdT506lQ0GgSRCggLAck6nU2vWrAlZpv+keQHHmljBO4Oovr7e7zwbDALJgwQFQExUVFSouro65BTiWAg1g4gNBoHkQYICIGYqKiq6nUJste5mEHk3GKypqYlZDACiR4ICIKbMTCG2ktkZRGbLAUgMEhQAPYrZGURmywFIDBIUAD1KdzOIHA6HSktLNWHChDhHBiAcJCgAktLIZdv9fi5b/KIk/xlEnZMU7/Hq1atZiRawORIUAD2OdwZRcXGx3/mSkhJVV1eroqIiQZEBMKtXogMAgFioqKhQeXm58vLyJEnbtm3T5MmTaTkBkgQtKABMcbvdcjgccjgccrvdEZeJpc7v3zEZmThxIskJkERIUAAAgO2QoABISkb7lyvBtta+43cMIPkxBgWAZTrPrElLz4zJ+5z44HUde+UJ3/HB6uVy5uRry/mh9/8BkDxIUAAkncNbH+hyznP8kGbPnp2AaADEAl08AADAdkhQAPQYgXYwBpCcSFAAAIDtkKAAAADbIUEB0GME2yAQQPJhFg+AhEtLz9TQ27eaKnNi3y4d3LwiTpEBSBRaUAAklT4jxmngjKVyZg/wO+/MydfTTz+doKgAWI0WFABJp8+IccoYep7q1lwlSRpYuVy9h12g6dMnJTgyAFahBQVAUnKkfbnxX2bp2X7HAJIfCQoAALAdEhQAAGA7jEEBEFNmZugAQGe0oAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAUzwej+/nnTt3+h17Ge1fnmutfcfvOB7MxAggOZCgAOiWy+XSqFGjfMdTp05VWVmZXC6XX5kDG27yHR+sXq76dXN1Yt+uuMR4Yt+uLjGOHj1av/zlL2UYhrKysuISBwBrOAzDMBIdRLiam5uVl5enpqYm5ebmJjocoEdzuVyqrKxU538qvDsHV1dXS1LAMl4DZyxVnxHjLI2r/WSran9cKUkaMO0WHd76QJcyHWOsqKiw9P0BhC+cv98kKACC8ng8KisrU11dXcDnHQ6HBg8eLElBy0inN/Ib/L0Nli5H3zFBcWYPkKflcNAYS0pKtH//fjmdLIcPJFI4f7/p4gEQVE1NTcjEwzAM1dXVhSwjSZ7jh9RW967V4X15/SDJiXQ6xtraWtXU1MTs/QFYjwQFQFAHDhyw7FqelqOWXSsSVt4LgNgjQQEQVFFRkWXXcmb3s+xakbDyXgDEHgkKgKAmTJigkpIS32DTzrzjO0KVkU6PQckoGR2rMOXMHhD0OYfDodLSUk2YMCFm7w/AeiQoAIJyOp1as2aNJHVJQLzHa9as8ZUJpv+keZYOkO2s7zevD3jeG+Pq1asZIAskGRIUACFVVFSourpaxcXFfudLSkp803e9ZTq3ZDhz8mMyxbizPsMv0cAZS7u8f8cYASSXXokOAID9VVRUqLy8XHl5eZKkbdu2afLkyX6tEhUVFSra+bnq1lwlSRpYuVy9h10Q05aTjvqMGKeMoef5vf/+TXfScgIkKVpQAJjS8Q/9xIkTA/7h75iMZJaeHbfkJNj7k5wAyYsWFAARK1v8YqJDANBDkaAASEpp6ZkaevvWRIcBIEbo4gEAALZDggIAAGyHBAUAANhOTBKU+vp6zZ49WwMGDFDv3r11zjnn6H//9399zxuGobvuuktFRUXq3bu3ysvL9cEHH8QiFAAx1H6yVR/dP00f3T9N7SdbEx0OgB7E8gTl6NGjGj9+vM444wz9z//8j9577z3913/9l/r1+3IfjlWrVumhhx7SunXrtGfPHmVlZWnKlClqbeUfOAAAEINZPPfff79KS0v1xBNP+M4NGzbM97NhGFq9erXuvPNOTZ8+XZL0s5/9TAUFBdq8ebOuvvpqq0MCAABJxvIWlF/96le66KKL9M///M8aNGiQLrjgAj3++OO+5/fv36+GhgaVl5f7zuXl5WnMmDHavXt3wGu2tbWpubnZ7wEAAHouyxOUDz/8UI899piGDx+uX//61/q3f/s3/fu//7ueeuopSVJDQ4MkqaCgwO91BQUFvuc6W7lypfLy8nyP0tJSq8MGAAA2YnmC0t7erq997WtasWKFLrjgAs2bN0833HCD1q1bF/E1lyxZoqamJt+jtrbWwogBAIDdWJ6gFBUVadSoUX7nRo4cqY8//liSVFhYKElqbGz0K9PY2Oh7rrOMjAzl5ub6PQAAQM9leYIyfvx47du3z+/cX/7yFw0dOlTS6QGzhYWF2rFjh+/55uZm7dmzR2PHjrU6HAAAkIQsn8WzcOFCjRs3TitWrNCVV16pN954Q+vXr9f69eslSQ6HQ1VVVbr33ns1fPhwDRs2TMuWLVNxcbFmzJhhdTgAACAJWZ6gfP3rX9fzzz+vJUuW6J577tGwYcO0evVqzZo1y1fmtttuk9vt1rx583Ts2DFdeuml2r59uzIzM60OBwAAJKGY7GY8bdo0TZs2LejzDodD99xzj+65555YvD0AAEhy7MUDAABshwQFAADYDgkKAACwnZiMQQHQ82RlZckwjJBl0tIzNfT2rXGKyH7vD8A6tKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKgIgZ7R7fz6217/gdA0A0HEZ3w/JtqLm5WXl5eWpqamJnYyBBXC6XrpwzT56Ww75zzpx89Z80T31GjEtgZKH9/b7LEx0CkLLC+ftNCwqAsLlcLlVWVvolJ5LkOX5IBzev0Il9uxIUGYCeggQFQFg8Ho8WLFgQck2UIzvW090DICokKADCUlNTo7q6upBlPMcPqa3u3ThFBKAnIkEBEJYDBw6YKudpORrjSAD0ZCQoAMJSVFRkqpwzu1+MIwHQk5GgAAjLhAkTVFJSIofDEbSMMydfGSWj4xgVgJ6GzQIBmFK2+EXfz20XfVdG3YqgZftPmidHmjMeYQHooWhBARC2PiPGaeCMpXJmD/A778zJ18AZS229DgqA5EALCoCI9BkxThlDz1PdmqskSQMrl6v3sAtoOQFgCVpQAESsYzKSWXo2yQkAy5CgAAAA2yFBAQAAtkOCAgAAbIcEBQAA2A4JCgAAsB0SFAAAYDskKAAAwHZIUAAE5Ha75XA45HA45Ha7Ex0OgBRDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAASClMnwaSAwkKAACwHRIUAABgO70SHQCA5JWWnqmht29NdBgAeiBaUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAVAShm5bHvAnwHYCwkKgIA8Ho/v5507d8po94QoDQDWIkEB0IXL5dKoUaN8x1OnTlX9urk6sW9XAqMCkEpIUAD4cblcqqysVH19vd95z/FDOrh5BUkKgLggQQHg4/F4tGDBAhmGEbTMkR3r6e4BEHMkKAB8ampqVFdXF7KM5/ghtdW9G6eIAKQqEhQAPgcOHDBVztNyNMaRAEh1JCgAfIqKikyVc2b3i3EkAFIdCQoAnwkTJqikpEQOhyNoGWdOvjJKRscxKmt1HD/TWvuO33RqAPZBggLAx+l0as2aNSHL9J80T440Z5wistaJfbt0YMNNvuOD1ctVVlYml8uVwKgABEKCAkBut1sOh0MOh0NTpkxRdXW1nNkD/Mo4c/I1cMZS9RkxLkFRRufEvl06uHmFPC2H/c7X19ersrKSJAWwGRIUAF1UVFSoaO6jvuOBlcs1+HsbkjY5Mdo9OrJjfeDnvphSXVVVRXcPYCMxT1Duu+8+ORwOVVVV+c61trZq/vz5GjBggLKzszVz5kw1NjbGOhQAYejYjZNZenbSdutIUlvdu/IcPxT0ecMwVFtbq5qamjhGBSCUmCYof/jDH/STn/xE5557rt/5hQsX6oUXXtBzzz2nV199VZ988okqKipiGQqAFGZ2WrTZadYAYi9mCUpLS4tmzZqlxx9/XP36fTklsampSRs2bNCDDz6ob33rW7rwwgv1xBNPaNeuXXr99ddjFQ6AFGZ2WrTZadYAYi9mCcr8+fN1+eWXq7y83O/83r17derUKb/zZ511loYMGaLdu3cHvFZbW5uam5v9HgBgVkbJaDlz8oM+73A4VFpaqgkTJsQxKgChxCRB2bRpk/74xz9q5cqVXZ5raGhQenq6+vbt63e+oKBADQ0NAa+3cuVK5eXl+R6lpaWxCBtAD+VIc6r/pHmBn/tizZfVq1fL6UzecTZAT2N5glJbW6sFCxbomWeeUWZmpiXXXLJkiZqamnyP2tpaS64LIHX0GTFOA2cs7TJ9uqSkRNXV1YyDA2yml9UX3Lt3rz799FN97Wtf853zeDzauXOnHnnkEf3617/WyZMndezYMb9WlMbGRhUWFga8ZkZGhjIyMqwOFUCK6TNinDKGnqe6NVdJOj19ev+mO2k5AWzI8gRl0qRJevvtt/3OXX/99TrrrLN0++23q7S0VGeccYZ27NihmTNnSpL27dunjz/+WGPHjrU6HADw03n6NMkJYE+WJyg5OTk6++yz/c5lZWVpwIABvvNz587VokWL1L9/f+Xm5ur73/++xo4dq0suucTqcAAAQBKyPEEx48c//rHS0tI0c+ZMtbW1acqUKXr00Ue7fyEAAEgJcUlQXnnlFb/jzMxMrV27VmvXro3H2wMAgCTDXjwAAMB2SFAA+G2St3PnTjbNA5BwJChAinO5XBo1apTveOrUqSorK9OJD9h6AkDiJGSQLAB7cLlcqqyslGEYfufr6+tl1D2QoKgAgBYUIGV5PB4tWLCgS3IiKeA5AIgnEhQgRdXU1Kiurq7bcgXfWaG0dGu2rQAAs0hQgBR14MABU+U8LUdjHAkAdMUYFCBFFRUVmSrnzO4X40gSq2zxi13O/f2+yxMQCYCOaEEBUtSECRNUUlIih8MRtIwzJ18ZJaPjGBUAnEaCAqQop9OpNWvWSFKXJMV73H/SPL/N9QAgXkhQgBRWUVGh6upqFRcX+50vKSnRwBlL1WfEuARFBiDVMQYFSHEVFRVasPNzac1VkqSBlcvlGHaB+vTQlpO09EwNvX1rosMA0A1aUAD4deNklp5Ntw6AhCNBAQAAtkOCAgAAbIcEBUBKaz/Zqo/un6aP7p+m9pOtiQ4HwBdIUAAAgO2QoAAAANshQQEAALZDggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtsJsxAHb4BWA7tKAAAADbIUEBkNKMdo/v59bad/yOASQOCQqAlHVi3y4d2HCT7/hg9XLVr5srl8uVwKgASCQoQEpyu91yOBxyOBxyu92JDichTuzbpYObV8jTctjvvOf4IVVWVpKkAAlGggIg5RjtHh3ZsT5kmaqqKnk8dPcAiUKCAiDltNW9K8/xQ0GfNwxDtbW1qqmpiWNUADoiQQGQcjwtR02VO3DgQIwjARAMCQqAlOPM7meq3MIX/q6yxS/GOBoAgZCgAEg5GSWj5czJD1nGmZOvjJLRcYoIQGckKABSjiPNqf6T5oUs03/SPDnSnHGKCEBnJChAChq5bHvAn1NJnxHjNHDGUjmzB/idd+bka+CMpeozYpzvHNOygfhjLx4AKavPiHHKGHqe6tZcJUkaWLlcvYddQMsJYAO0oABIaR2TkczSs0lOAJsgQQEAALZDggL0cIHGT7BBHgC7YwwKkGJcLleXDfKcOfnqP2me38BQAEgkWlCAFLJlyxZVVlYG3CDv4OYVOrFvV4IiAwB/JChACrnttttkGEbQ54/sWE93j6T2k6366P5p+uj+aWo/2ZrocICURIICpJD6+vqQz3uOH1Jb3btxigYAgiNBAeDH7EZ6ABBLJCgA/JjdSC+VeDxfdnvt3LnT7xhAbJCgAClk8ODBcjgcQZ9ng7yuTnzwukaNGuU7njp1qsrKyuRyuRIYFdDzkaAAKWTVqlUhn2eDvK4Ob32gy9id+vp6VVZWkqQAMUSCAvRwHbsj+vXrp2effdbUBnmpIi09U0Nv36qht29VWnqmqdd4Z0JVVVXR3QPECAkK0MN0XDl248aNXbonFi5cqLxvzPGdG1i5XIO/tyElk5NoGIah2tpa9erVix2OgRhgJVmgB5s9e3aXdU/q6+tl1D3oO2aDPAB2RAsK0IMFWpQt1EJtAGAXlicoK1eu1Ne//nXl5ORo0KBBmjFjhvbt2+dXprW1VfPnz9eAAQOUnZ2tmTNnqrGx0epQACAiZlbTDTUbCkD0LE9QXn31Vc2fP1+vv/66XnrpJZ06dUqTJ0/266NduHChXnjhBT333HN69dVX9cknn6iiosLqUOIq0I6xAJLPiX27/DZTBJAYlo9B2b59u9/xk08+qUGDBmnv3r2aOHGimpqatGHDBm3cuFHf+ta3JElPPPGERo4cqddff12XXHKJ1SEBgCkn9u3Swc0rui3nzMnXz9at0axZs+IQFZCaYj4GpampSZLUv39/SdLevXt16tQplZeX+8qcddZZGjJkiHbv3h3rcICUQjeEeUa7R0d2rO+2XP7MuzT4exs0ffr0OEQFpK6YzuJpb29XVVWVxo8fr7PPPluS1NDQoPT0dPXt29evbEFBgRoaGgJep62tTW1tbb7j5ubmmMUMIDW11b0rz/FD3ZZzONKY9QTEQUxbUObPn6933nlHmzZtiuo6K1euVF5enu9RWlpqUYRAz/b0008HXJRtwLRbEhSRfZndJNHjPhbbQABIimGCcvPNN2vr1q363e9+p5KSEt/5wsJCnTx5UseOHfMr39jYqMLCwoDXWrJkiZqamnyP2traWIUNJL3OK8cWXv+I79i3KNtwxnp1ZnaTRGdWX0lsIAjEmuUJimEYuvnmm/X888/r5Zdf1rBhw/yev/DCC3XGGWdox44dvnP79u3Txx9/rLFjxwa8ZkZGhnJzc/0eALpyuVxdVo5teOJm3zGLsgWXUTJazpz87ssNHqkT+3axgSAQY5aPQZk/f742btyoLVu2KCcnxzeuJC8vT71791ZeXp7mzp2rRYsWqX///srNzdX3v/99jR07lhk8QAdut1vZ2dmSpJaWFmVlZYUs73K5VFlZ2WUhNk/L4S5lvfvP4EuONKf6T5rX7Syez/72Bx3e+kCX8/X19Zo5c6Ykc78vAKFZ3oLy2GOPqampSd/85jdVVFTke/ziF7/wlfnxj3+sadOmaebMmZo4caIKCwuT/n8eNPcikTwejxYsWNDtKrFmFiBLZX1GjNPAGUu7jtvpcHzslScCvrZj3fP9B6IXky6eQI/rrrvOVyYzM1Nr167VkSNH5Ha75XK5go4/SQaBmtVp7kU81dTUqK6urttybfXvxyGa5NZnxDgVzX3UdzywcrnfcaAWqc5ee+21mMQGpBL24omSt1m9vr7e73x9fb0qKytJUhAXBw4cMFWOGSjmdBynE8m4nWBLJgAwjwQlCqGa1b3nqqqqaO5FzBUVFZkq552BgthK5hZhwC5IUKLQXbO6YRiqra1VTU1NHKNCKpowYYJKSkq6XTk2Y/DIOEXUs3gHFQ+5dYupmT7jx4+PQ1RAz0aCEgWzzepmywGRcjqdWrNmjaTQy9szxTg63pk+Utd67njsdFLPQLRIUKJgtlndbDkgGhUVFaqurlZxcbHf+c4zUhAd70yfzvU8ePDgBEUE9EwkKAG43W45HA45HA653e6g5bprVnc4HCotLdWECRPCuq5V8SG5RTJ1vaKiQu+9957veNu2bRr8bz/V0Nu3aujtW5WWnhmTWFNNnxHjutTzn/70J98xSw0A0SNBiUKoZnXv8erVq2nuRdginbpetvhFnb38Jd/xjS+doFsnRjp+r48ePapzzjnHd8xSA0D0SFCiFKxZvaSkRNXV1aqoqEhQZEhWTF1PDiOXbff9PGvWLH5fgMVIUCwQqFl9//79ppITumvQkZmp6zNnzuTzEmPeWTvRdIux1AAQHRIUi3Rs7p04cSLdOoiImanrSB4sNQBEjgQFsBGmpPdM/F6B8JGgBBCrjf+sui4bE/ZcTEnvmfi9AuEjQenEyo3/Oo4v2bhxY8DrbtmyJWHxmcX06PgxM3UdyaPzUgMAzCNB6SCWsydmz54d8LqzZ8+2RXywBzNT12FPLDUAWIsE5Qux3vgv1HXtEB/sI9TU9aeffrrb11sxAwXd61jPgVaWZakBIDokKF9I1MZ/HROOUONJ2JgwtQSbuj5t2jTfOcYf2UeglWXNLjUAIDASlC/YYeO/YONJ3G63LrvsMlPXCBVfMo4BScaYrdJ56vqWLVviPv4I5qXSUgOp/L1E/JCgfCHajf+ysrJkGIYMw1BWVlbEcUQ7noTZAj3Tli1bAo4/qqur08yZMzXo/92hssUvJig6ALAeCcoXwt34z4xImt+7G09iZXxmMT068W677baQY5aO7Fgvo536BNBzkKB8weqN/zpPBw5Hd+NJ4jlbwKppzYmYHt2TdG456cxz/JDa6t6NUzQIpOPePCOXbadFC4gSCUoHVm38F2w6cLgCjSd5+umnTcVnRR+xVdOaU216dKC6j0ef/aljn8bkugCQCCQonUSz8Z8UejpwuLzjSTp2hfTr109vv/12RPGF08Vi1bTmaK9Dt5B5zqy+iQ4BKYLvJeKBBCWAaEbjdzcd2Cs/P9/UeJJAXSPnnHNO2PGF28Vi1bTmaK6T6t1C3oHXn3/+ecjxUV4Zg0fGKTIEkirrz6T69xLxQ4JiMbPTkK+66ipJoceTBJu58cknn3S5XqguhGDX8XaxBHpdtNOuvfFEOj06nG6heHapJGJ6ZajxUR050nrutNZkZebzYrcpu6HisUt3rd3qDLFBgmIxs9N8L7/88pDjXaZPn95t14gUfKZQx/NmrtNZtNOuw9XxOsm+am6g5u9om8QrKiqUP32J0rL6+53veNxa+w4zeRAzyf69RPIhQbFYd9OVvcaPHx9yvIvZrqLXXnuty7nOC3odOnQo6OuDJSlWTrsePHhwWNdJ5lVzAzV/FxQUaNiwYX7nImkS7zNinIrmPuo7zh17pTpW68Hq5apfN1cn9u2K/AaAIJL5e4nkRIISQKDpgp0fwZjd7M07biTYeBezXSwNDQ1dzgXamDBcVk67XrVqVVjXscOqvpEI1vx9+PBhHTlyxO9cpE3iHbtxmnc/K0+L/3U9xw/p4OYVJCmwXLJ+L5G8SFAiFChp8T4WvZERsDl+8ODBpq9vtuuksLBQkn+3QqQziDp3PUQz7brzzKNnn33W9HXM3vs111wjt9ttWZeKmX7tYNcNd/aWd9XhmTNnqrm52dRrwsHCbfZw1h1f/mdm2HUPqLm5uctnrPNnKlAZq8ZcRPMZD7fbN5bjRGL1/Ya9kKAEYMVo/M7N8QMrlyvtO2t91x39g1dUtvjFoK01c7Y1y5mTH/T63q6RyZMnR7UoXEeBuh4imXYdqJtj4cKFuu+++0xdx0z3kleg/Wms6lIxc1/e65rtkgskUDddtFi4LfFO7NulAxtu8h0frF7e5Xsa6DNlxXc5UqE+47FYbdvqGNGzkKDEUMfm+MzSs8OaZeFIc6r/pHlBnzcMQ60XXquiyrs0c+bMqLt0vAJ1PYQz7TrUKP/vfve7pq5jtptMCtydZWWXild3sxe2bNkS0XWlwN10UtdWunB5Wo5GHBOic2LfLh3cvEKelsN+5zvOwAtnll48mPmMW7nadixiJEnpWUhQbKzPiHEaOGOpnNkD/M47c/I1cMZS9R4+Rkd2rLf0Pb1dD1ded6OG3varsF5rZpS/WaG6l55++umwrxvNLAMz9/XMM8+Edc2OvN10VnNm94vJdRGa0e4J+r3s+Bm69dZbI56lZzWzM3SmT59uyWrbsYyRWUQ9BwmKCe0nW/XR/dP00f3T1H6yNei5SHTXnRSoq2jw9zaoz4hxaqt7V57jwWfoRMNz/JA+2/9myG6ojo8hi36pXr16dTvKv6Pu1i+ZMmVKwO6ladOmRXRP3lkGv/nNbwLfc5B+bTOzFw4ePKiBAwd2O3srkPb2dtP/qKalZ2rIrVtCdv9Jp5PYjJLRYceC6Jn9XpppKXnssccsmabuFeg6brfb1HfXO0PHbLevmZjDGRcS7SyieK5+m6jtLnoaEpQkEKyrKNZN+B73sZhde+Sy7QETn87nzl7+ku/4xpdOdBlzEolAXSqh+rXNzkqYNWuWpNCLqQUybdo0ZfQr1KD/d4epLp3uuv8kqf+keSzcliBWfi8XL15s2ZiqYJ/xcLonvd+F7rp9YzFOJJpZRIxbSU4kKDEU66Wvw2nCd+bka8C0W8K7/hd7u1h9HwXfWRHRdU588Lol4206d6l016/9wQcfmLpusObvAQMGqH///kFedVq404O76/7rM2KcqevAelZ3rVkxpirUZ3z27NmmYzEzkydW40QiXTyScSvJiwTFhI7TNb2rdQY6F28ZJaO7beqXpPyZd53uFhp+ie9c5z9sgRhGu+n7MluuY9dDuPV67JUnTL1Hd7674XUNve1XKlv8oobe9itded2NQfu1DcPQf/7ooZD11XH2QkVFhRyVD/qeG1i5XFn/8lP1ufYxU7GFMz04VPcfEsfs9zIaVm7YaWYcl9kZOuGOEwmn2yWSWUSJGrcSi9WkUxEJSjcCTRWse3i2PvnJXL9zVq3gGc7YFjNN/ZLUe8i5XZr7+37z+m5fd+iX95i6r851FIq368FsvXYs03lGRKQ63peZ8QLtLYdDvrd3RtWZd5zupgrUJWe2u8Vz/JA+/tF00+Oaopkphtgw+72MVucxF8HGOEQzBV4Kb4aO2XEivXr10saNG8Pqduludp/32q2tX3534rX6bce6D3RfVi59kEpjWUhQQgg2VbC99bjaW1v8zsVyBc9QXSxBm/oD/I+/43WyR38z4Os66+6+gtVRZx27HsKq1y/K5Fw0PeT1O0rLzFFaZnbIMr77+mCP6esOmHYLXSowJdj3Mha6G5sR7cqugWboeHfaNgxDWVlZEb1XoCUCuut2CTa7L9gimIlY/TZeSx+kAhKUIEJNFQwl2hU8I+k6CtTU3/HY7OtCCXRfZuvI18U0YlzE9ep+93emyuVNnKOS7z+t4hs3WHpdSUrLzFbh9Y/4jmPZpcLGf8kvnO9XNLxjLqJdATaQrVu3drswY6BYzIi02yXQLKI//elPvuNoVr+NVCQreUfaxZRKXUUkKEFEOoU3mhU8A3V7mO06irSpP5yuh873ZbaOHI403/tEWq/tnzUrrXdut+VyLpgaVpeK2etKp38fDU/c7DuOZZcKG//1DLHscus45iKaFWBDSUtLC2vhtWjey8tMt0vHmI4ePapzzjnHdxzv1W+jWck73C6mVJuNRIISRDRTBSN5bbBuj0i7jmIxg8h7X95xMo0/X2rudR2mK0dTr1mjL+u2TCR/EMxc1ytQV1bncUOB6t57buAMc3Ummfvdx3qmGOyp47iQYCvSersQOq5xEm7i4J2Sb3bcQ6hxIuG67LLL5Ha7u33vUF1FsVj9tvN4k0B1Hy4zXUxWz0ZKhrEsJChBRDNVMNzXmun2sMPmb977CjcO73TljteIRJ/hY0yPt7Hiut3x1kM43XKRjE2ww+8e0Qs2hqnj9P9AZQKNqfKOC5k+fbrpGTpLly4Nuysj0Gak3XUrBBsnEomdO3fq5MmTId873qvfdnz/cDYIDaW730ssZiMlQ1cRCUoQkU4VjGQFTzPdHone/M17X+HM2PHKGDzyy5+jrNdIx9uEe10z2urfj6hbLtz3SvTvHtboM/ySwNPCO0z/D1Sm85iqgZXL5bj6ES16I0ODr11leobOD3/4Q7/j/Pz8bls5xo8fH1G3QudxIpGaOnWq8vPz/Y7NLi4Xyeq33elcF4cORbeSt9kuJqtnIyVLVxEJShCRThWMZAVPs90eocpF2tRvtuuh/6R5+uyDPaZm7HQ0cMZSOTO/HOVvRb12Hm/jzMyKuEsl2HXNaP37mxF3y8XqMwJ76fy9DDRWzEwZq1aT7ri8vnHJdZICd304HA798pe/1G9/+9uIuxWs2jSwvb29y3vPnj1bVVVVpl5vdvXb7gTrYolUOF1MVs5GSqaF60hQQgjWHB+oyTWa6aZmuz1iufmb1RsThqqPeNVrqPfyNq0f3Lwi4v2UupsBZGXXDBv/IZBoPhe9z/y68qcvUVqW/yrHadkDlD99iRa+3qvbRQy9m4p69+PyJjdDFv3Sb9sKK4W7Qec111xjeoxFsHEZobpYzAi0mvTgwYNlGIZmzpwZdL8e7/E111xj6n2i7SryxtPc3GzyzmKrV6IDsLs+I8YpY+h5qltzlaTTzau9h12g9pOtXc5FOmLf2+0RqpsnHpu/BbtXR5pTrR//yfTsGzP1EU29ev/HGe19dUxKWmvfUe9hF2jo7VtltHtUv25uyPtN652r9s9Cf4m9XTOZQ87t+vov7sHMe7HxX89h5rMbqEyw15n5tyOYU5/+LervfMfPeOexWOF8n8Ll3aDTzPdQkoZd94AObLqz2wSj87iMyZMny+l0RrXY3datW/Xtb39bbrdbeXl5kk53MV1yySW+pMX7Xp3f/5JLvuz+y8/P1+HDhwPeg8PhUElJSdRdRV6PPfaYbrnlFstawSJFgmJCOE2ukV6//6R5Orh5RdAy8dr8zYqmZLP1Eet6DfVen32wR0d++xPfuYPVy+XMyVf/SfPUZ8S4bn8fWaMv0/H/7b4fvLt6s9PvHsnHzOcnGO9nM9rvvKflqE7s2xXV9ylSZr+HB6uXK6PfI754Aul8D1OnTvXdg+E5FXGM/7ajVWk12/3+Q3Ttoy8r84Yb/N5rwIABfsnH1KlTlZb2ZSdHsPEusegqWrx4sR555BGtWbMm7HE6ViJBsQlvV8SR3/7Eb0xDxy95IvWkLoYTH7yuw1sf6HLeO3bE26UU6veR1jvb1D+MZurN7r972Fuwz093uvtsmv3Onzpar6bfb+xy3tT3KXtAVFtY9Bk+Rpklo03de+d4OvIu8xDsNXnjZ4UVV3f3Fejfn8OHAyxh0Gn8TSBp2QPUf9I8LXojQ4veCL4LuiS1fvz3bq/n5R2TEslMJ6s4DCvmSMVZc3Oz8vLy1NTUpNxcc4tshSPUVvex5ml1W9Z1ZCWzXRGDv7fBFvEGE+59BPt9xKI+7Pq7h/0E+qx07B7truvDzGfTVFdn9gA5FHqfLKu+T2avG869d3xdKGbu0yvUd9fq7q6BlcuVXjRC9Q9/J+R7h1M/HXm7jvbv329Zd084f78ZJGszdt38zczsm2Toigh3Snew30cs6sOuv3vYS7Cp7Z/97Q++c/0s+Gya+YznnPdtUy0XVnyfOgs2A8/s2LBwlkxobzms7POnmCob6rsb6UraQeNqbfFb3TrQMged79NsciJZt5lipGhBCSCRLSh25+2nTdauCPd7r+rQCz/qtlz+Fbcqa9Q3ui2X7PWB5BKsK8KrY/eFVZ/NUNcxPKdi/n2SI00yvuzqCHUPZr/fORdNN9VF21H+FbfK4Twjqjo1G58Z3d2Dd4kFK8b+bNy4Ud/5zneivo4U3t9vxqAgLKFG/ScDq6d0J3t9IHmYXXG69/AxcqQ5LftsdjfTx4xovk9mujDCfZ9wNgnteO3MIedGVadWjuXr7h4ObX9YDlnT/hDtZoqRoosHYUvmrggzK9mGO603mesDySOSFaet+mwGu048vk9pvc4I+N6BmIknnDEYXh3vIZo6jXQl7c7M3IPRelztrS1RvY8VmylG9f508XRFF0/PFk4zOWAXVndPWsVu36fu4omke8fKe+guPjMiuYdIruVwOCyfxcMgWSCE7lbNJTmBHdlhxelA7PZ96jae4WNMX8uZk6/86YuV1jtb7vdeVevHf4p6dehwVtKWw/9PdCT30G08wTZhzclP6BRjKcFjUNauXasf/ehHamho0HnnnaeHH35YF198cSJDQoroM2Kceg8fc7rZvOWonNn9lFEymu4Z2JZdVpwOxG7fp1DxGO2ebuvRkZmjgdNvV3tri46+/N9+Za0YBB8sPkl+59KLz9LJT/4c0T2Y4f28ONKcAeOpqPinqK4frYR18fziF7/Qd7/7Xa1bt05jxozR6tWr9dxzz2nfvn0aNGhQyNfSxQMgFdmtOyVZmalHKfQMmETXtRVdRd3dw9/vuzyq6weSFF08Dz74oG644QZdf/31GjVqlNatW6c+ffropz/9aaJCAgBb83UPdBpoSfdkeLqrRzObo1q5GWgkgt2DGcnyeUlIF8/Jkye1d+9eLVmyxHcuLS1N5eXl2r17d5fybW1tamtr8x03NTVJUsx2XGxvOxGT6wJAtDLLzlfR9Q+r7ZP35Wk5Jmd2X2UUj5Qjzcm/XWEIVY+f7d9rasbUZ/v3KrPk7DhF3FXHe/j8+BEde/UpGa3Hg78gI1v5l1cp84tune4+L7H4G+u9ppnOm4QkKIcOHZLH41FBQYHf+YKCAv35z3/uUn7lypX6z//8zy7nS0tLYxYjAAChHKy+J9EhhKetRYdc95ounrc6dqEcP37ct7tzMEmxUNuSJUu0aNEi33F7e7uOHDmiAQMG+HZytEpzc7NKS0tVW1sbk/EtOI16jg/qOT6o5/ignuMnVnVtGIaOHz+u4uLibssmJEHJz8+X0+lUY2Oj3/nGxkYVFhZ2KZ+RkaGMjAy/c3379o1liMrNzeULEAfUc3xQz/FBPccH9Rw/sajr7lpOvBIySDY9PV0XXnihduzY4TvX3t6uHTt2aOzYsYkICQAA2EjCungWLVqkOXPm6KKLLtLFF1+s1atXy+126/rrr09USAAAwCYSlqBcddVVOnjwoO666y41NDTo/PPP1/bt27sMnI23jIwM3X333V26lGAt6jk+qOf4oJ7jg3qOHzvUdVLuxQMAAHo29uIBAAC2Q4ICAABshwQFAADYDgkKAACwnZRMUNauXauysjJlZmZqzJgxeuONN0KWf+6553TWWWcpMzNT55xzjrZt2xanSJNbOPX8+OOPa8KECerXr5/69eun8vLybn8vOC3cz7PXpk2b5HA4NGPGjNgG2EOEW8/Hjh3T/PnzVVRUpIyMDH31q1/l3w4Twq3n1atXa8SIEerdu7dKS0u1cOFCtba2xina5LRz505dccUVKi4ulsPh0ObNm7t9zSuvvKKvfe1rysjI0D/8wz/oySefjHmcMlLMpk2bjPT0dOOnP/2p8e677xo33HCD0bdvX6OxsTFg+ddee81wOp3GqlWrjPfee8+48847jTPOOMN4++234xx5cgm3nq+55hpj7dq1xptvvmm8//77xnXXXWfk5eUZdXV1cY48uYRbz1779+83Bg8ebEyYMMGYPn16fIJNYuHWc1tbm3HRRRcZU6dONX7/+98b+/fvN1555RXjrbfeinPkySXcen7mmWeMjIwM45lnnjH2799v/PrXvzaKioqMhQsXxjny5LJt2zbjjjvuMFwulyHJeP7550OW//DDD40+ffoYixYtMt577z3j4YcfNpxOp7F9+/aYxplyCcrFF19szJ8/33fs8XiM4uJiY+XKlQHLX3nllcbll1/ud27MmDHGjTfeGNM4k1249dzZ559/buTk5BhPPfVUrELsESKp588//9wYN26c8d///d/GnDlzSFBMCLeeH3vsMeMrX/mKcfLkyXiF2COEW8/z5883vvWtb/mdW7RokTF+/PiYxtmTmElQbrvtNmP06NF+56666ipjypQpMYzMMFKqi+fkyZPau3evysvLfefS0tJUXl6u3bt3B3zN7t27/cpL0pQpU4KWR2T13NmJEyd06tQp9e/fP1ZhJr1I6/mee+7RoEGDNHfu3HiEmfQiqedf/epXGjt2rObPn6+CggKdffbZWrFihTweT7zCTjqR1PO4ceO0d+9eXzfQhx9+qG3btmnq1KlxiTlVJOrvYFLsZmyVQ4cOyePxdFmttqCgQH/+858DvqahoSFg+YaGhpjFmewiqefObr/9dhUXF3f5UuBLkdTz73//e23YsEFvvfVWHCLsGSKp5w8//FAvv/yyZs2apW3btumvf/2rbrrpJp06dUp33313PMJOOpHU8zXXXKNDhw7p0ksvlWEY+vzzz/W9731PS5cujUfIKSPY38Hm5mZ99tln6t27d0zeN6VaUJAc7rvvPm3atEnPP/+8MjMzEx1Oj3H8+HFde+21evzxx5Wfn5/ocHq09vZ2DRo0SOvXr9eFF16oq666SnfccYfWrVuX6NB6lFdeeUUrVqzQo48+qj/+8Y9yuVx68cUX9YMf/CDRocECKdWCkp+fL6fTqcbGRr/zjY2NKiwsDPiawsLCsMojsnr2euCBB3Tffffpt7/9rc4999xYhpn0wq3nv/3tb/r73/+uK664wneuvb1dktSrVy/t27dPZ555ZmyDTkKRfJ6Liop0xhlnyOl0+s6NHDlSDQ0NOnnypNLT02MaczKKpJ6XLVuma6+9Vv/6r/8qSTrnnHPkdrs1b9483XHHHUpL4//gVgj2dzA3NzdmrSdSirWgpKen68ILL9SOHTt859rb27Vjxw6NHTs24GvGjh3rV16SXnrppaDlEVk9S9KqVav0gx/8QNu3b9dFF10Uj1CTWrj1fNZZZ+ntt9/WW2+95Xv80z/9ky677DK99dZbKi0tjWf4SSOSz/P48eP117/+1ZcAStJf/vIXFRUVkZwEEUk9nzhxoksS4k0KDbaZs0zC/g7GdAiuDW3atMnIyMgwnnzySeO9994z5s2bZ/Tt29doaGgwDMMwrr32WmPx4sW+8q+99prRq1cv44EHHjDef/994+6772aasQnh1vN9991npKenG9XV1caBAwd8j+PHjyfqFpJCuPXcGbN4zAm3nj/++GMjJyfHuPnmm419+/YZW7duNQYNGmTce++9ibqFpBBuPd99991GTk6O8fOf/9z48MMPjd/85jfGmWeeaVx55ZWJuoWkcPz4cePNN9803nzzTUOS8eCDDxpvvvmm8dFHHxmGYRiLFy82rr32Wl957zTjW2+91Xj//feNtWvXMs04Vh5++GFjyJAhRnp6unHxxRcbr7/+uu+5b3zjG8acOXP8yj/77LPGV7/6VSM9Pd0YPXq08eKLL8Y54uQUTj0PHTrUkNTlcffdd8c/8CQT7ue5IxIU88Kt5127dhljxowxMjIyjK985SvGD3/4Q+Pzzz+Pc9TJJ5x6PnXqlLF8+XLjzDPPNDIzM43S0lLjpptuMo4ePRr/wJPI7373u4D/3nrrds6cOcY3vvGNLq85//zzjfT0dOMrX/mK8cQTT8Q8Todh0A4GAADsJaXGoAAAgORAggIAAGyHBAUAANgOCQoAALAdEhQAAGA7JCgAAMB2SFAAAIDtkKAAAADbIUEBAAC2Q4ICAABshwQFAADYDgkKAACwnf8PnGX4Ai7lLKIAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "m.visualize()" ] } ], "metadata": { "kernelspec": { "display_name": "venv", "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.9.16" }, "orig_nbformat": 4, "vscode": { "interpreter": { "hash": "bdbf20ff2e92a3ae3002db8b02bd1dd1b287e934c884beb29a73dced9dbd0fa3" } } }, "nbformat": 4, "nbformat_minor": 2 }